孔词，肖明

(武汉大学水资源与水电工程科学国家重点实验室，湖北武汉，430072)

随着我国西部大型水电站建设，水电站地下洞室水工混凝土结构向着体形复杂化，地质条件多样化以及受力条件复杂化的方向发展[1]。高效、合理地计算地下洞室复杂水工混凝土结构的配筋是保证水工混凝土结构安全稳定、经济环保设计的关键。

对承受内水压力为主的有压隧洞，现行电力行业“水工隧洞设计规范”的衬砌结构配筋基本是按照“边值法”和“公式法”进行[2]。

1)“边值法”假定混凝土衬砌在内水荷载作用下不开裂，按照“厚壁圆筒”理论计算结构配筋，当水头较高、围岩弹性抗力系数较低时，衬砌结构将发生开裂，此时“边值法”计算结果往往偏大。

2)“公式法”假设混凝土衬砌已经开裂，混凝土不承担拉力，这和实际的水工隧洞钢筋混凝土衬砌的受力状况相符合，但“公式法”一律采用钢筋允许应力代替实际应力，导致低水头有压隧洞的配筋常出现负值，同时“公式法”用弹性抗力表征围岩的力学性能，不能反映围岩进入塑性以后的真实承载情况[3]。

梁倩等[4−5]针对“边值法”和“公式法”计算衬砌配筋的缺陷，提出了相应的改进配筋方法，但均将围岩作为弹性介质，没有真实地反映围岩和衬砌结构的联合承载机理。随着对衬砌结构承载机理的深入研究，人们认为衬砌不是主要承载体，大部分荷载应由围岩承担[6−8]。目前，在混凝土结构设计领域，对于形状简单，受力明确的混凝土结构可以按照规范[9−10]来计算结构配筋；但对于结构型式和受力条件复杂的水工混凝土结构，目前缺乏相应的设计规范指导混凝土结构的配筋计算[11]。刘霞等[12]将有限元和GES0 算法相结合，寻找钢筋混凝土D 区构件的优化解，建立了D 区平面构件的拉压杆计算模型并给出了钢筋的布置方案。黄博豪等[13]根据应力图形法对结构进行初步配筋，然后将钢筋加入混凝土结构后进行非线性计算，并开发了一套非线性有限单元法配筋系统，但该方法不满足现行设计规范的要求，在工程领域没有被广泛采纳。

可见，不论是DL/T5195—2004“水工隧洞设计规范”还是SL 191—2008“水工混凝土结构设计规范”，均不能很好地指导考虑围岩承载的大型复杂地下洞室复杂水工混凝土结构的配筋；而采用钢筋混凝土非线性有限元计算复杂水工混凝土结构的配筋不满足现行规范的要求。

鉴于此，本文作者首先以弹塑性损伤本构模型[14]模拟围岩的力学特性，采用线弹性本构关系计算水工混凝土结构的应力和应变，提出考虑围岩与结构联合承载的三维非线性有限元计算方法；其次，提出适宜编程的水工结构任意复杂断面轮廓线及相应配筋截面的提取方法，并以结构应力为基础构建分段插值函数，得到配筋截面上的应力分布图形，依据截面的应力分布图形按照现有规范对“杆系结构”和“非杆系结构”进行配筋计算；最后，开发复杂水工结构的自动化配筋软件，并将其应用于拉哇水电站复杂闸室结构的配筋实践中，验证了该方法的合理性和高效性。

1 三维有限元计算及其应力的处理

考虑围岩对内水压力的承载作用，本文作者采用三维弹塑性损伤有限元模拟洞室的开挖和支护效应，将混凝土材料视为弹性介质，将围岩作为弹塑性损伤介质，通过非线性迭代计算，得到混凝土结构的应力和应变。

1.1 洞室围岩和混凝土衬砌结构的三维有限元计算方法

岩体作为一种弹塑性介质，在洞室开挖后，其应力状态将发生改变，易产生卸荷，导致开裂破坏[15]。本文采用Z-P屈服准则作为岩体的弹塑性判别准则，考虑岩体的损伤破坏特性，岩体的应力和应变关系为

式中：为考虑损伤后的应力；Dep为弹塑性应力矩阵；εij为岩体的应变张量；sij为应力偏张量；δij为Kronecker 符号；D为岩体的损伤变量，可由下式计算得到

式中：eD=R为岩体的损伤常数。

洞室开挖以后，为方便迭代计算，将洞室开挖荷载R

分解为弹性荷载和塑性荷载2个部分

式中：σ0为岩体的初始地应力；Re和Rp分别为弹性荷载和塑性荷载；s为弹性系数。对弹性荷载Re进行一次加载计算；对塑性荷载Rp进行分级加载和分级迭代计算[16]进而模拟洞室围岩的开挖效应，并采用文献[17−18]中方法模拟锚杆、锚索的支护效应。

考虑到洞室开挖完毕后一段时间，混凝土衬砌结构才施加上去，为反映围岩变形对衬砌结构的影响，将式(4)中塑性荷载分解为2个部分

式中：φ为开挖释放系数，可根据岩体的开挖释放特性确定[15]。

衬砌施加前，采用空气单元模拟衬砌单元，并对整个结构施加Rp1进行迭代计算。衬砌施加后，将衬砌空气单元的力学参数用混凝土材料的力学参数替代，重新形成刚度矩阵，并将Rp2和内水压力施加到结构上进行迭代计算，进而得到混凝土结构的应力和应变。

1.2 混凝土衬砌结构有限元计算的应力处理

有限元计算的位移在单元边界处能够保证连续，但应力通常是不连续的，且单元内部应力精度高于单元边界应力精度[19]。有限元计算的单元应力通常为积分点上的应力，为便于后续的断面应力插值计算，本文采用绕节点平均计算单元节点应力，并通过区域峰值应力修正单元节点应力。

根据有限元计算的单元形心处的主应力及其方向，采用斜截面公式σen求出单元形心点处沿截面法向的正应力σen，按照下式计算得到单元节点的截面法向正应力σnn。

式中：H为与节点相关的单元数目；(σen)i为与节点相关的第i个单元形心点处的截面法向正应力。

采用绕节点平均得到单元节点应力的方法原理简单，操作方便，但存在应力均化现象，无法反映实际的峰值应力[20]。为了反映峰值应力，引入应力修正系数λ，根据下式对单元8 个节点应力进行一次修正。

式中：λi为第i个单元的应力修正系数；(σen)max表示与第i个单元的8 个节点相关的所有单元形心点处的截面法向正应力的最大值；(σnn)max表示第i个单元的8个节点通过绕节点平均后计算的截面法向正应力的最大值。

2 结构断面轮廓及配筋截面提取方法

复杂混凝土结构断面轮廓以及配筋截面的提取是数值配筋的关键，针对8节点六面体单元，可按下列步骤进行。

1)通过程序自动计算切割面和六面体单元十二条棱边的交点个数以及每个交点的全局坐标，然后判断交点连接方式，使交点连接形成凸多边形。

2)对所有相关单元进行循环遍历，得到切面切割有限元模型的所有线条信息，对所有线条进行去重操作，可得到切面切割有限元模型的轮廓线。

3)判断轮廓线的数目，对封闭构件和非封闭构件进行自动识别，进而根据轮廓线对配筋截面进行提取。

4)通过斜率突变准则对角点截面进行提取，从而反映断面转角处的应力集中现象，使配筋结果更安全可靠。

2.1 复杂断面轮廓提取方法

为减小计算量，复杂断面轮廓线的提取分3步。

首先，确定全局坐标系下切割面的平面方程

式中：a，b，c和d为空间平面一般式方程的系数，可根据实际配筋需求通过点法式或三点式方程进行确定。根据式(8)和整体有限元模型进行单元初选，将单元的8个节点坐标代入式(8)中得到8个数值，若所有数值均大于0 或者均小于0，说明该单元与平面没有交点，否则至少有一个交点，将与平面有交点的单元记录到一个列表中。对所有单元进行遍历，从而实现单元初选。

其次，对经过单元初选得到的单元，分别求取单元十二条棱边的参数方程，

式中：α，β和γ分别表示单元某条棱边线段与全局坐标系的xg，yg和zg轴的夹角；(xg0，yg0，z0g)为线段其中一个端点的全局坐标；L0为线段的长度。将式(9)代入式(8)进行求解，得到切割面和单元所有棱边的交点个数和对应的坐标，然后判断交点的连接方式，使其形成凸多边形。任意平面切割六面体形成的凸多边形类型见图1所示[21]。

图1 六面体和任意平面相交的切割平面Fig.1 Cutting planes crossed by arbitrary planes for hexahedron

最后，对所有单元进行循环初选计算，将确定的每个初选单元和切面的交点连接起来，得到切面切割有限元模型的所有线条信息，再进行线条的去重操作，便得到了断面轮廓线。

2.2 配筋截面的提取方法

提取断面轮廓线后，可以判断封闭轮廓线的数目，确定切割面截取的构件是封闭构件还是非封闭构件。若封闭轮廓线数目为2(封闭轮廓线数目大于2的构件本文不进行考虑)，则为封闭构件；若封闭轮廓线数目为1，则为非封闭构件(如图2所示)。

对于具有2组封闭轮廓线的封闭构件，可以通过比较轮廓线沿某一方向(如xg方向)的最大宽度，区分出内轮廓和外轮廓(如图2(a)所示)。取内轮廓所有线段的中垂线和外轮廓相交进而得到若干配筋截面。

对于只有1组封闭轮廓线的非封闭构件，需分割该组封闭轮廓线，形成一条较长的非封闭轮廓线1 和一条较短的非封闭轮廓线2，然后取较短轮廓线2上所有线段的中垂线和较长轮廓线1相交进行配筋截面的提取(如图2(b)所示)。

根据开口断面的突变特点，按以下步骤挑选出非封闭构件的分割点：

1)计算每个初选单元的相邻单元数目和单元信息，若对应于某个初选单元的相邻单元数目为1，则该单元位于开口断面突变处，如图2(b)中I和II所示。

图2 配筋截面提取Fig.2 Extraction of reinforcement sections

2)求解切割面和单元I、单元II 在开口突变处的交点信息，进而挑选出分割点①和分割点②。

2.3 角点截面的提取及配筋截面优化

鉴于六面体单元剖分的特点，通过上述方法常无法提取角点截面并对其进行配筋计算，而断面角点处有应力集中现象，往往是最危险截面，更应该进行配筋计算分析。为提取断面角点处的配筋截面，本文以斜率突变为准则自动识别所有角点，并判断角点连接方式，实现角点处配筋截面的自动提取。

图3所示为角点识别与配筋截面优化。依次计算断面2条轮廓线上线段的斜率，若2条相邻线段的斜率差值大于一阈值δ0，则将2条线段的公共点视为一角点，遍历所有线段得到所有角点

图3 角点识别与配筋截面优化Fig.3 Corner recognition and reinforcement section optimization

式中：Pin和Pex分别表示内侧轮廓和外侧轮廓的角点集合。以内轮廓中某一角点为初始点，向外轮廓角点集合中寻找使2个角点间距离最小的角点作为终点，连接起始点和终点，构成角点截面，对内轮廓角点集合进行循环遍历，进而得到所有角点截面。为便于结果显示，控制相邻配筋截面的间距不小于一阈值δ1，优化所有配筋截面。

3 结构配筋截面的内力计算方法

地下洞室混凝土结构受力条件复杂，沿配筋截面的应力分布常表现出明显的非线性特性，尤其是在断面转角处。董正中等[22−23]在研究中假定两插值点间应力为线性分布，无法反应截面应力分布的非线性特性。为此，本文通过分段二次插值原理构建截面正应力分布函数，既能反映截面应力分布的非线性特性，又能克服等距高次插值容易产生龙格现象[24]的缺陷，同时减小计算的复杂度。

图4所示为截面应力分布图形。对每一个配筋截面建立局部坐标系，坐标原点为配筋截面内轮廓一侧的端点，x轴和截面重合，y轴和截面外法线平行。在配筋截面上等距选取n(n为奇数)个应力插值点，可以通过插值原理构造出(n- 1)/2 个分段二次函数，从而用分段函数描述截面法向正应力的分布情况，图4中以9个应力插值点为例。

图4 截面应力分布图形Fig.4 Sectional stress distribution graph

首先，由插值点的全局坐标(xg，yg，zg)根据等参单元逆变换算法[25]计算应力插值点的自然坐标(ξ，η，ζ)。

然后，根据下式计算插值点的截面法向正应力

式中：Nj为第i个单元的第j个节点的形函数；(σnn)j为第i单元第j个节点的截面法向正应力，由式(6)计算得到；λi为第i个单元的节点应力修正系数。

每个应力插值点的局部坐标为

式中：L为截面线段的长度；(σn)k为第k个插值点经过应力修正以后的截面法向正应力。

对于第l个分段函数，构造插值函数的3 个插值点分别为P(2l-1)，P(2l)和P(2l+1)，将其分别记为p1(x1,y1)，p2(x2,y2)和p3(x3,y3)。这3 个插值点通过拉格朗日插值可得到分段函数的表达式，各分段函数的系数为

式中：al，bl和cl为第l个分段插值函数的系数。

构造出截面上的应力分段函数后，对函数进行积分便可得到截面的轴力。对于二次多项式，可以直接应用辛普森公式精确求积，将每个分段函数求积结果相加可得整个截面上的轴力，按照下式求取截面的轴力N

式中：n为插值点个数；yk表示第k个插值点的截面法向正应力。截面弯矩通过应力分段函数对截面中心取矩得到，对其中某一分段函数，应力对截面中心的矩计算如下

式中：Ml为第l个分段插值函数对截面中心的矩；fl(x)为第l个应力分段函数。将所有分段函数对截面中心的矩相加便可得到截面弯矩。

4 数值配筋方法及自动化配筋软件平台

地下洞室中的混凝土结构形体复杂、外形尺寸大、空间整体性强，不宜简化为“杆系结构”计算内力[26]。混凝土衬砌结构在内水压力作用下，极易发生开裂破坏，但此时外荷载产生的内力与弹性方法计算所得结果差别不大[27]。为此，本文考虑围岩和混凝土结构的联合承载特性，建立包含围岩和混凝土结构的三维有限元模型，对围岩进行非线性计算，对混凝土衬砌进行线弹性计算，得到混凝土结构的应力和应变。基于截面应力分布图形，按照规范[9−10]进行结构配筋。

根据混凝土结构的形状尺寸和受力特性，混凝土结构可划分为“杆系结构”和“非杆系结构”，“杆系结构”和“非杆系结构”的配筋计算分别采用“结构内力法”和“弹性应力图形法”，配筋计算流程如图5所示。

图5 自动化配筋流程Fig.5 Flow of automatic reinforcement

实际工程中的部分混凝土结构，如水闸底板、船闸闸室、矩形渡槽的槽身等结构可以简化为“杆系结构”计算配筋；但像水工结构中的各种深梁、短梁以及泄水建筑物的闸墩等结构，无法简化为“杆系结构”计算配筋[28]。为此，本文开发了“杆系结构”和“非杆系结构”的配筋计算程序，对于可简化为“杆系结构”的混凝土构件，可计算出截面内力，按照偏心受拉、压承载力计算截面配筋；对于大体积混凝土结构等“非杆系结构”，根据截面应力分布图形按照“弹性应力图形法”计算截面配筋。

为了实现复杂混凝土结构内力、配筋计算以及结果绘制的自动化，本文基于Python 编程语言和AutoCAD绘图平台，开发自动化配筋软件平台，通过一次性处理多个断面，大幅提高实际工程中配筋计算效率。自动化配筋软件平台包括GUI 界面、主计算程序模块和批量绘图模块3个部分，图6所示为软件平台的功能结构图。

图6 软件平台功能结构图Fig.6 Software platform function structure diagram

GUI 图形用户界面基于wxPython 库开发，用户可交互式输入计算文件路径、配筋计算参数以及绘图参数，计算成功后点击开始绘图按钮，程序将自动关联到电脑的AutoCAD并新建一个DWG文件，并将配筋结果绘制于DWG文件中，状态栏中将显示程序计算及绘图的结果。

5 工程实例分析

5.1 有限元模型和计算参数

以拉哇水电站导流隧洞闸室结构为例，建立了包含闸室结构和围岩的三维整体有限元计算模型如图7(a)所示，闸室结构有限元模型如图7(b)所示。

图7 有限元模型和分析断面Fig.7 Finite element model and analysis sections

模型岩体材料分为3 类，混凝土设计强度为C35，钢筋选用HRB400，材料的物理力学参数如表1所示。计算工况为弧门挡水工况，计算荷载包括衬砌结构自重，上游有压段内水压力(水头86.24 m)以及弧门水推力。其中，内水的水头按面荷载施加在有压段衬砌结构内表面，弧门水推力按集中荷载施加在闸墩大梁结构的支铰处。

表1 材料的物理力学参数Table 1 Physical mechanical parameters of materials

5.2 闸室结构配筋计算分析

闸室上游侧有压隧洞承受较大内水压力，混凝土结构在承受内压作用下极易发生破坏；弧形闸门在弧门挡水期承受巨大的水推力，闸墩大梁混凝土结构的稳定性控制了闸室结构的整体安全稳定性。故选取闸室上游侧1-1断面和闸墩大梁结构2-2 断面作为分析断面(如图7(b)所示)，计算分析混凝土结构的配筋。

闸室上游侧1-1 断面可简化为“杆系结构”，按照“结构内力法”计算结构配筋。配筋截面是以衬砌厚度为高，宽度为1 m 的矩形截面。1-1 断面配筋结果见图8。

由图8(a)可知：1-1 断面顶拱的应力分布比较均匀且数值相对较小；而断面底角处应力集中现象明显，沿衬砌厚度方向应力变化剧烈，最大拉应力超过C35混凝土的极限抗拉强度，混凝土发生开裂。由图8(b)～(d)可见：计算的断面内力图分布规律和有限元计算的应力分布规律吻合，闸室上游侧由于承受较大的内水压力，导致衬砌结构承受环向拉力，其中轴力在洞室转角处由于应力集中的缘故明显增大，最大值达2 657.16 kN；弯矩在顶拱和洞室转角处表现为外弯，其余部分表现为内弯，其中在底角处，弯矩显著增大，达−433.13 kN·m。截面的应力分布在顶拱处基本呈现出线性分布规律，在边墙和底板处存在一定的非线性分布特性，在转角处，截面应力分布的非线性特性显著加剧。衬砌结构的配筋率在边墙和底板处呈现出外侧较大、内侧较小的规律，而顶拱和转角处则相反，尤其在底角部位，内侧配筋率显著大于外侧配筋率。若采用对称配筋，则建议配筋率取0.77%；若采用非对称配筋，则建议衬砌内侧配筋率取0.77%，衬砌外侧配筋率取0.48%。

图8 1-1断面配筋结果Fig.8 Reinforcement results of 1-1 section

2-2 断面闸墩大梁混凝土结构体积和截面厚度较大，应采用“弹性应力图形法”计算截面配筋。目前规范中没有明确规定大体积混凝土结构的最小配筋率，只规定厚度大于2.5 m的卧置在地基上以承受竖向荷载为主的底板单位宽度内的配筋面积不小于2 500 mm2[9−10]。为此，本文的闸墩大梁结构的单位宽度内的最小配筋面积暂取2 500 mm2。2-2断面配筋结果如图9所示。

由图9(a)可见：在弧门推力作用下，大梁结构主要表现为受压，结构的最大压应力为−3.00 MPa，出现在大梁结构的支铰附近，结构压应力未超过混凝土的极限抗压强度。由图9(b)可见：大梁背部出现大量拉应力区，拉应力由梁体表面向梁体内部逐渐减小，最大拉应力为0.45 MPa。由图9(c)可见：拉应力分布图形和大梁结构的第三主应力分布规律相符合，拉应力图形面积在大梁背部中点部位达到最大值614.10 kN/m。由图9(d)可见：大梁结构单位宽度内的配筋面积均为2 500 mm2/m，大梁结构配筋可按2 500 mm2/m进行，但部分截面的受拉区向梁体内部延伸较长，尤其在大梁顶部，受拉区几乎贯穿整个大梁结构，为较危险区域，工程中应重点关注。

图9 2-2断面配筋结果Fig.9 Reinforcement results of 2-2 section

6 结论

1)在内水压力作用下，有压段衬砌结构承受较大环向拉应力，衬砌结构极易发生开裂破坏，最大配筋率高达0.77%；衬砌截面应力在顶拱处呈现线性分布规律，在底板和边墙处存在一定的非线性分布特性，尤其在衬砌转角处，截面应力分布的非线性特性显著加剧；衬砌边墙和底板的配筋率呈现出外侧较大、内侧较小的规律，而顶拱和转角处则相反，尤其在底角部位，内侧配筋率显著大于外侧配筋率；若采取非对称配筋，则建议取内侧配筋率为0.77%，外侧配筋率为0.48%，若采取对称配筋，则建议配筋率取0.77%。

2)在弧门水推力作用下，大梁大体积混凝土结构以受压为主，但压应力尚未超过混凝土的抗压强度；大梁背部出现拉应力区，拉应力由梁体表面朝梁体内部逐渐减小，拉应力图形面积均不大，大梁结构配筋可按照2 500 mm2/m进行，但大梁顶部拉应力区延伸较长，几乎贯穿整个结构，需重点关注。

3)本文所提出的数值配筋方法对各种复杂水工结构均具有较好的适应性，能够反映衬砌结构在复杂荷载条件下的非线性应力分布特性，能对“杆系结构”和“非杆系结构”进行配筋计算和自动绘图，解决了人工难以处理的复杂水工结构的配筋难题，为地下工程中的复杂水工混凝土结构的配筋设计提供了一个有效手段。