复杂隧洞衬砌结构的三维有限元分析及结构配筋

董正中，肖 明，杨 莹(武汉大学水资源与水电工程科学国家重点实验室，武汉 430072)

0 引 言

随着我国大型水电站建设的展开，水电站隧洞衬砌结构向着体形大型化，地质条件多样化，受力条件复杂化的方向发展[1]。衬砌自身的结构形式也日益复杂。如何分析衬砌结构的力学特性，以及如何计算衬砌结构的配筋，使之同时满足安全性和经济性的要求，是一个重要的课题。压力隧洞的传统计算方法[2,3]假定围岩及衬砌材料为连续各向同性的弹性介质，忽视了围岩及衬砌材料的非线性性质的不利影响，并且是将围岩与衬砌结构分开进行计算，从而导致隧洞混凝土衬砌厚度及配筋量偏大。因此，有必要对考虑隧洞围岩、衬砌联合承载及材料非线性情形下的衬砌结构进行分析计算。基于此，本文首先采用三维弹塑性损伤有限元分析方法，对复杂受力条件下衬砌结构进行计算，研究探讨了各工况下隧洞衬砌结构的受力特点；其次利用有限元插值应力场在衬砌截面上积分计算内力，基于内力结果在配筋截面上进行配筋计算，以保证衬砌配筋满足要求，并提出优化建议。

该方法基于三维有限元计算，综合考虑了初始地应力场，开挖卸荷，复杂地质条件及地下水位等多种因素，研究了衬砌结构的稳定特性，并直接由有限元结果插值进行配筋计算。为在工程实际中选择经济合理的衬砌形式以及提高隧洞运行的安全度提供依据。相较于传统衬砌结构设计评价方法更高效，精确，可为类似工程提供借鉴与参考。

1 衬砌结构三维有限元分析的基本方法

本文对衬砌结构采用弹性模型进行计算；对围岩单元则采用弹塑性模型[4]，屈服准则采用Ziekiewicz-Pande准则，以增量变塑性刚度法进行迭代计算[5]。

衬砌结构是在隧洞开挖完毕后施加上去的 。为了反映隧洞开挖变形对衬砌作用的影响，可以洞室开挖释放的荷载:

(1)

改写成：

(2)

式中：σ0是开挖单元的初始地应力场，包括自重应力场和构造应力场；γ为岩体容重；α为荷载分配系数。

衬砌结构施加前，作用于结构的荷载为

(3)

衬砌结构施加后，作用于结构的荷载为：

(4)

式中荷载分配系数α的取值与岩性、地应力和支护时机等因素相关，工程中一般靠经验确定[6]。因此，本文尝试给出一种确定取值的数值计算方法。

弹塑性有限元计算中，开挖后围岩单元的应力状态按下式可分为弹性和塑性，即：

(5)

式中：F为屈服函数；σ为开挖后的围岩应力；σ0为岩体的初始应力；Δσ为开挖引起的应力增量。当岩体开挖后，若围岩单元进入塑性状态，则一定存在一个临界应力状态，满足：

F({σ})=F({σ0}+β{Δσ})=0

(6)

式中：β为单元的弹性系数。令p=1-β，p称为塑性系数，其大小反映了总应力增量Δσ中塑性荷载的比例。由式(6)可以看出，初始地应力场下已屈服的岩体单元p=1；开挖后仍处于弹性的单元p=0。

根据式(6)，采用Zienkiewicz-Pande屈服准则[7]，通过一次开挖计算，令全部开挖释放荷载作用于围岩，求出所有围岩单元的塑性系数p；塑性系数p的分布在一定程度上反映了围岩的承载条件，因此，式(2)中的荷载分配系数 可根据洞室顶拱和边墙单元的塑性系数p分布按式(7)综合确定。

α=(1-η)p

(7)

其中，η为支护时机滞后系数。视工程实际情况，对于自承能力较差的岩体，数值分析时可取用较小η值，使围岩单独承担的荷载减少；对于自承能力较强的围岩，数值分析时可适当放大η取值，使围岩承担部分塑性荷载。一般地，岩性越差，地应力越大，则η越小；反之，η越大。

洞室开挖荷载释放完毕后，再根据运行期引水隧洞的水位分布情况，计算隧洞衬砌受内、外水压力作用的受力特征 。

2 衬砌结构的配筋计算

本文基于应力图形法，编制了地下洞室衬砌结构的配筋计算程序。并可以自动生成配筋截面：对于衬砌结构有限元模型，先生成模型的拓扑关系，即结点、单元线、单元面、单元间的相互包含关系；再对模型进行消隐，建立表面单元线表，这里定义表面单元面仅有一相关单元或相关单元仅有一衬砌单元，表面单元线包含在表面单元面中，以一平面截取表面单元线，生成内外轮廓交点，依序连接交点，即形成内外轮廓线，即可获得配筋截面。

2.1 应力修正与内力修正

由于有限元计算的特点，单元内部有应力均化现象，从而导致配筋截面上计算弯矩偏小。引入应力修正系数r。首先对结点应力作估值，可由与该结点相关的单元应力取平均值求得。分别计算单元应力及结点应力在问题域局部的极值，令其为σeij、σnij，取

r=‖σeij‖/‖σnij‖

(8)

这里‖ ‖表示应力张量的度量，取为应力不变量的函数。一般有r>1。将r与结点应力相乘即得修正结点应力。

有时沿配筋截面仅有少数几层单元，这可能导致计算弯矩偏小甚至方向相反。引入应变修正方法，假设在配筋截面方向上，应力依线性分布，设应力梯度为沿配筋截面应力的斜率，则应力梯度可由配筋截面两端点应变计算获得。假设配筋截面两端点沿所取投影方向正应变分别为ε1和ε2，则可得配筋截面上应力梯度为：

k=γ(ε1-ε2)/L

(9)

式中：L为配筋截面长度；L依结构形状和受力情况取值。

2.2 配筋截面的内力及配筋计算

单个截面的内力计算过程如下：

(1)设置应力插值点，取配筋截面分点作为插值点，份数与单元层数成比例，一般为其2倍或以上[8]。

(2)由上述插值方法计算各插值点应力，并根据上述应力修正方法对应力进行修正。

(3) 由柯西公式σn=ninjσij计算各插值点沿某一方向正应力。通过数值积分计算配筋截面上内力，这里采用复合积分公式。依据沿配筋截面方向单元层数对截面上内力作修正[9,10]。

以一长圆筒为例。其内径r=5 m，外径R=6 m，受外压q1=0.5 MPa，内压q2=0.7 MPa，按弹性本构计算，则由解析解公式可得配筋截面上轴力N=500 kN，弯矩N=16.56 kN·m。建立该圆筒有限元模型，沿径向取5层单元，对其进行有限元分析。图1左右分别为经由插值及数值积分所得的轴力及弯矩结果。可以看出，与解析解相比，轴力误差在4%以内，弯矩误差在6%以内。若考虑应力修正取应力张量度量为第一应力不变量可得应力修正系数为1.01，从而使结果更趋近于解析解。

图1 圆筒模型截面轴力弯矩图Fig.1 Axial force and bending moment graph of cylinder model

截面上配筋参照《水工混凝土结构设计规范》(DL/T 5057-2009)。如果截面上应力分布接近线性，按正截面承载力方法计算配筋；若偏离线性较大，按拉应力图形面积计算，即As=KT/fy，这里T取为拉应力面积， 为钢筋强度设计值，K为承载力安全系数。裂缝宽度验算依据具体工程情况参照相应规范计算。

3 工程实例分析

3.1 工程概况

本文结合某水电站引水隧洞衬砌结构进行分析计算。电站装机容量2 400 MW。电站单机引用流量较大(Q=621.4 m3/s)，机组采用单管单机供水形式，对应4台机组。每条引水道由进口渐变段(矩形→圆形)、上平段(有压引水隧洞)、渐缩段(有压隧洞段→地下压力钢管段)、上弯段、斜井段、下弯段、下平段等组成。引水道进口底板高程1 575.00 m，出口中心高程为1 494.80 m。由于引水隧洞条数较多、洞径较大(最大开挖断面16.8 m，居已建发电引水隧洞最大开挖直径的前茅)、相邻隧洞间的岩柱厚度较小(略大于一倍开挖洞径)，且围岩地质条件较复杂，因而有必要对隧洞衬砌结构受力特性进行分析，并进行结构配筋计算以满足限裂要求。

3.2 计算参数及有限元模型

计算网格一共剖分了362 544个8节点空间等参单元，节点总数376 960个。由四条引水隧洞进口段始端建至有压段末端。有限元模型参见图2，根据引水隧洞地址剖面图，共划分了三类岩体单元。衬砌采用C25混凝土，衬砌结构参见图3。计算所取的岩体和衬砌参数见表1及表2。

图3 衬砌单元模型Fig.3 Lining element model

表1 地下洞室围岩参数Tab.1 Surrounding rock parameters of cavern

表2 衬砌材料物理力学参数设计值Tab.2 Design values of mechanical parameters of lining

模型X轴垂直于进水口段轴线方向，Y轴沿进水口段水的流向，Z轴与大地坐标重合，指向上为正，范围由1 433.6 m高程至地表。边界条件为：底部全约束，四周法向约束，顶部不约束。

3.3 计算工况

本文以四条隧洞均过流为最不利工况展示衬砌受力及配筋成果。

(1)内水压力由正常蓄水位1 618 m确定。水击压力，4号压力管道末端的最大水击压力值为36.8 m，其余3个引水道的压力管道末端的水击压力值也按4号引水道取，为36.8 m，水击压力分布按线性规律计算。

(2)外水压力由水库蓄水后地下水位线确定，水库蓄水后的地下水位线由渗流场分析求得[11,12]。根据计算所得的岩体渗流场节点水头，在衬砌结构上施加渗透体积力。

(3)计算使用自行编制的地下洞室三维有限元分析以及渗流分析程序。

3.4 计算结果分析

(1)衬砌应力及位移分布规律。进水口渐变段1-1以及上弯段2-2两段截面衬砌第三主应力分布规律参见图4及图5。两段典型截面变形示意图参见图6及图7。

受内水压力作用，隧洞典型断面第一主应力均为压应力，分布在-1.14～-0.19 MPa，应力矢量沿洞周径向。

图4 衬砌1-1断面第三主应力分布图Fig.4 3rd principal stress of 1-1 section of lining

图5 衬砌2-2断面第三主应力分布图Fig.5 3rd principal stress of 2-2 section of lining

图6 衬砌1-1断面变形示意图(单位：mm)Fig.6 Displacement of 1-1 section of lining

图7 衬砌2-2断面变形示意图(单位：mm)Fig.7 Displacement of 2-2 section of lining

断面第三主应力为拉应力，应力矢量沿洞周切向。1-1断面由于断面形式，四周拐角处出现应力集中，拉应力值达1.5 MPa；2-2断面由于高程比1-1断面高程低，内水压力较大，断面各部位第三主应力量值均大于1-1断面第三主应力。4号洞由于洞径最大，洞周围岩参数较其他三条隧洞洞周围岩参数低，故而衬砌承担较多围岩荷载，断面各部位拉应力量值大于其他隧洞拉应力值，两侧腰部拉应力值达2.5 MPa。计算结果表明，渐缩段以及下游衬砌拉应力超出混凝土抗拉强度设计值，应重点加强配筋。

1-1断面位移分布在0.2～0.7 mm。2-2断面各部位位移值相对于1-1断面位移略大，分布在0.3～1.2 mm。变形指向洞外，各断面底部位移大于隧洞顶部及两侧位移。其中圆形断面洞周变形较渐变段洞周变形均匀。

(2)隧洞衬砌配筋计算限于篇幅，本文展示4号隧洞1-1断面和2-2断面的截面内力以及配筋分布图，分布参见图8和图9。

图8 4号隧洞1-1断面配筋示意图Fig.8 Reinforcement drawing of 1-1 section of lining

图9 4号隧洞2-2断面配筋示意图Fig.9 Reinforcement drawing of 2-2 section of lining

1-1断面受截面形状影响，角点处弯矩与边上弯矩方向相反，由于内水压力，边上弯矩向内弯，近角点处弯矩向外弯。其余各断面弯矩基本向外。

各洞基本按小偏心受拉配筋，并利用裂缝宽度限值验算配筋率。按所处环境类别，取裂缝宽度限值为0.3 mm，先由单宽配筋截面上内力进行配筋计算，结果表明：1-1断面按构造配筋，2-2断面裂缝宽度均未超出限值。若对于其他断面出现裂缝宽度超出限值，需对这些截面通过裂缝宽度验算求得配筋率。

4 结 论

(1)考虑隧洞围岩、衬砌联合承载及材料非线性情形下对衬砌结构进行分析计算。为了反映隧洞开挖荷载对衬砌作用的影响，提出了“荷载分配系数”取值的一种数值计算方法，为类似工程计算提供参考。

(2)基于有限元分析结果，实现了地下洞室隧洞衬砌结构的配筋。采用应力图形法配筋，配筋截面自动生成，能适应各种截面形式，同时应力分析与配筋过程相对独立，具有较高的效率，对圆筒结构的验算表面，该方法具有较高的精度。

(3)采用应力修匀计算结点应力，再计算插值点应力。特别地，该方法对于网络稀疏的情况具有很好的精度。对某水电工程引水隧洞的计算实例表明，衬砌结构有限元分析以及配筋计算方法成果合理，可以用于实际工程的配筋设计中。

□

[1] 卞 康, 肖 明, 任 祎. 复杂调压井结构的三维有限元分析计算[J]. 中国农村水利水电, 2008,(3):115-117.

[2] 河海大学, 武汉大学, 大连理工大学,等. 水工钢筋混凝土结构学[M]. 北京：中国水利水电出版社, 2009.

[3] 吕西林, 金国芳, 吴晓涵. 钢筋混凝土结构非线性有限元理论与应用[M]. 上海：同济大学出版社, 1997.

[4] 肖 明. 地下高压岔管洞室稳定及衬砌损伤开裂分析[J]. 武汉水利电力大学学报, 1995,28(6):594-599.

[5] 俞裕泰, 肖 明. 大型地下洞室围岩稳定三维弹塑性有限元分析[J]. 岩石力学与工程学报, 1987,6(1):47-56.

[6] 肖 明. 地下洞室施工开挖三维动态过程数值模拟分析[J]. 岩土工程学报, 2000,22(4):421-425.

[7] 王勖成, 邵 敏. 有限单元法基本原理和数值方法[M]. 北京：清华大学出版社有限公司, 1997.

[8] 张雨霆, 肖 明, 熊兆平. 三维空间离散点数据场的插值方法[J]. 武汉大学学报: 工学版, 2008,41(4):34-37.

[9] 钱向东, 任青文, 赵 引. 一种高效的等参有限元逆变换算法 Ξ[J]. 计算力学学报, 1998,15:4.

[10] 李春光, 郑 宏, 葛修润, 等. 六面体单元等参逆变换的一种迭代解法[J]. 岩土力学, 2004,25(7):1-050.

[11] 张 巍, 肖 明. 地下工程复杂渗流场数值模拟研究与应用[J]. 岩土工程技术, 2005,19(1):13-16.

[12] 张 巍, 肖 明. 地下工程渗流断层数值模拟的隐式复合材料单元法研究[J]. 岩土工程学报, 2005,27(10):1 203-1 206.