高压水工隧洞透水衬砌渗流-应力-损伤耦合分析方法研究

2018-04-17周亚峰文喜雨伍鹤皋

水利学报 2018年3期

周　利，苏　凯，周亚峰，文喜雨，伍鹤皋

（1.武汉大学水资源与水电工程科学国家重点实验室，湖北武汉　430072；2.长江勘测规划设计研究院，湖北武汉　430010；3.长江岩土工程总公司，湖北武汉　430010）

1　研究背景

随着社会经济与工程技术的发展，水工隧洞无论在数量上还是在规模上都达到了前所未有的高度，带来了巨大的社会与经济效益，同时也面临诸多问题。例如：隧洞在内水压力作用下，衬砌混凝土受拉开裂，并在施工缺陷以及温度微裂纹等影响下，形成性态不一的宏观裂缝，进而导致内水外渗［1-2］，此类衬砌结构均称为透水衬砌。在高内压下衬砌环向拉应力超过混凝土抗拉强度后，衬砌混凝土损伤破坏，表现为微裂纹萌生、扩展并贯通，直至产生宏观裂缝，洞内高压水沿裂缝通道外渗使得渗透水压增大，在影响结构应力场分布的同时又加剧了损伤，如此相互影响，是一个复杂的多场耦合过程。卞康等［3］采用增量变弹性损伤刚度计算迭代的方法分析了衬砌水压致裂过程中的耦合效应；Simanjuntak等［4］探讨了衬砌开裂过程中水力耦合效应对隧洞渗流量的影响；Bian等［5］在分析实际工程时，提出了一种渗流-应力耦合方法；Fahimifar等［6］在平面轴对称假定下，提出了一种解析方法以分析衬砌开裂过程中的耦合效应。然而从我国进行的6次较大规模的隧洞衬砌直接压水试验结果来看（详见表1），当隧洞内压增大到某一数值时，裂缝出现，此后无新的裂缝产生，只是缝宽随内压增大，裂缝呈少而宽分布［7］，而这一现象在现有的研究中却鲜有体现。考虑到水工隧洞衬砌裂缝分布的这一特点，刘秀珍［8］按照弹性地基曲梁结构推导了圆形有压隧洞衬砌裂缝宽度的计算方法；沈威［9］采用室内水工隧洞结构模型试验研究了衬砌开裂过程，同样得出了裂缝少而宽的分布特征。

表1　隧洞压水试验参数及结果［7］

衬砌裂缝出现后，洞内高压水势必会沿着裂缝向外流动，使得围岩变形增大，当衬砌因钢筋的约束作用径向变形与围岩不一致时，两者可能出现脱离现象，单独承担自身范围内的水压力［10-11］，体现了衬砌与围岩典型的有条件联合承载机理［12-13］。此外，工程监测资料显示，隧洞发电运行时，衬砌钢筋应力普遍小于设计预期，并没有达到既安全又经济的效果，说明现有的衬砌设计方法仍存在一定的问题，以致高估了钢筋应力值。例如：在广蓄电站高压水工隧洞充水加压初期，钢筋应力随着内压增大；衬砌开裂后，钢筋应力值变化明显放缓甚至出现下降的趋势，最大值仅为49.5 MPa［14］。类似的情况也发生在周宁水电站水工隧洞充排水试验中［15］。然而目前的衬砌设计理论、解析及数值方法尚不能很好地解释这一现象［11］。

鉴于此，本文基于混凝土塑性损伤模型和透水衬砌理论，考虑在高内压作用下衬砌损伤开裂、渗流与应力的耦合作用以及衬砌与围岩有条件联合承载特性，基于大型通用有限元软件平台ABAQUS进行二次开发，通过实用程序GETVRM获取材料损伤，利用子程序USDFLD实现材料渗透系数随损伤的动态更新，提出并实现了透水衬砌渗流-应力-损伤耦合算法，模拟分析了水工隧洞初次充水加压过程，以此研究衬砌损伤开裂过程中的多场耦合作用，探讨钢筋应力的演化特征，为水工隧洞衬砌设计问题提供一定参考。

2　理论基础

水工隧洞水压力由衬砌和围岩承担。本文采用弹塑性本构关系模拟围岩，以Mohr-Coulomb准则作为岩体的破坏准则，由于研究重点在于透水衬砌的耦合分析，限于篇幅，在此不做赘述。对于衬砌，本文采用间接耦合的方法对隧洞初次充水期衬砌进行渗流-应力-损伤耦合分析，考虑渗流场通过等效节点荷载影响应力场分布，在结构自身所承担的荷载以及等效节点荷载的作用下计算应力场后即可获取结构应力、应变以及材料损伤，并建立损伤与渗流之间的耦合关系，根据损伤实时调整材料渗透系数并返回渗流场进行计算，直至计算收敛、加载结束，详见图1。

图1　渗流-应力-损伤耦合关系

2.1渗流场与等效节点荷载的数值解法假设水和材料不可压缩，水在材料孔隙间的渗流满足微分方程［16-17］：

式中：θ为含水率；分别是x、y、z方向上的渗透系数；t为时间变量；水力势H=p/γ+z，其中p为孔隙水压力，z为位置水头，γ为水容重。

计算空间域经过离散、插值和积分，应用Galerkin加权余量法及格林公式可得求解渗流场的矩阵方程：

式中：［Ks］为总渗透矩阵；{h}为节点水头列向量；{A}为对渗流边界积分得到的节点荷载。

渗流场计算完成后，水力梯度所产生的水荷载以渗透体积力的方式作用在结构内部，按式（3）换算成作用在模型节点上的等效荷载［3］：

式中：［N］为插值函数；Ω为节点积分域。

2.2应力场的计算将计算空间域离散化，得到渗流与应力的耦合方程：

式中：［Km］为结构刚度矩阵；{u}为节点位移列向量；{F}为结构的节点荷载。

2.3混凝土塑性损伤模型及钢筋模型相较于预裂缝模型与分布裂缝模型，塑性损伤模型能够较好地模拟混凝土材料软化过程、卸载刚度退化及反复加载过程［18-19］，并且不必预先假定混凝土开裂位置，程序会根据结构的受力状况判定并计算结构各个部位的损伤开裂情况，亦不会出现大面积开裂的情形，在工程分析中应用广泛［20］。多轴应力状态下，混凝土应力σ计算如下：

式中：d为单元损伤；为初始弹性刚度；ε为单元应变；εpl为塑性应变。

单元损伤d由等效塑性应变表示如下：

式中：t、c代表拉、压状态。

在ABAQUS中，混凝土损伤后，应力、损伤与塑性应变的关系可以按照列表形式输入，本文衬砌材料采用C25混凝土，抗拉、抗压强度取标准值1.78 MPa、16.7 MPa，结合《混凝土结构设计规范》（GB50010-2010）附录C确定相关参数，具体应力、损伤与塑性应变之间的关系如图2所示。

采用埋藏式模型模拟钢筋，在迭代计算中钢筋由线性杆单元模拟，选取理想弹塑性本构关系，应力σs与应变εs关系如下：

图2　塑性损伤演化曲线

式中：Es为钢筋弹性模量；fy为屈服应力。

2.4损伤开裂与渗流的耦合关系采用塑性损伤模型时，损伤将导致混凝土结构刚度矩阵的退化，以此模拟混凝土损伤开裂对结构受力状态的影响，同时混凝土损伤开裂对结构的渗透特性亦会带来显著的影响。材料单元由损伤相与未损伤相两部分组成，单元渗透系数计算如下［21］：

式中：km为未损伤相渗透系数；kd为损伤相渗透系数；损伤相的塑性体积应变为单元塑性体积应变；若单元同时出现拉损伤与压损伤，分别计算对应状态下的单元渗透系数，k取较大值。

考虑到脆性材料出现宏观裂缝后，其渗透特性较开裂前会出现突增现象，因此引入突跳因子ξ以计算损伤相渗透系数［22］：

式中：对于压剪损伤ξ=100，对于拉损伤

3　技术路线

在耦合计算时，首先根据式（1）—（3）计算渗流场与等效节点荷载，然后根据式（4）计算应力场，并依据塑性损伤演化曲线获取材料损伤，根据式（11）和式（12）调整材料渗透系数后重新计算渗流场，如此循环直至计算收敛、加载结束。

为提高计算效率并保证计算精度［3］，本文内压分3个阶段进行施加以模拟隧洞初次充水运行工况：（1）通过试算确定混凝土即将出现损伤时所对应的内压，该内压采取一次性加载的方式进行施加；（2）衬砌损伤后，内压水头增量以较小值逐级缓慢施加；（3）贯穿缝出现后，内水外渗通道形成，适当增大内压水头增量。每施加一次内压水头增量，须保证在该级内压作用下耦合迭代完全，渗流场达到稳定，以保证耦合精度。

本文基于ABAQUS平台，通过Fortran语言调用子程序USDFLD进行二次开发，使用USDFLD时，必须调用实用程序GETVRM以获取材料点的数据，以此完成耦合迭代过程，具体步骤如下（详见图3）：

①施加上覆压力以模拟初始应力场的影响，计算隧洞开挖后应力场，考虑隧洞开挖后，应力释放完全后进行衬砌支护，即衬砌不承担开挖释放荷载，计算衬砌支护后应力场；

②衬砌内壁施加增量水头计算渗流场，将各单元渗透体积力转换为等效节点荷载；

③将等效节点荷载施加于模型节点计算得到应力、应变以及材料损伤，调用实用程序GETVRM获取单元损伤值，对损伤单元根据式（11）和式（12）调用子程序USDFLD调整其渗透系数；

④根据新的渗透系数计算渗流场，检查渗流场是否稳定。当前后两次渗流场计算结果满足式（13）时，渗流场达到稳定状态；如不满足，则继续求解等效节点荷载施加于模型节点进行耦合迭代，直至渗流场稳定；

式中：Pn+1、Pn为迭代n+1和n次后节点孔压；δ为迭代控制阀值，为提高计算效率并保证计算精度，将δ取为 1%［3］。

⑤施加下一级增量水头，在同一级内压水头下重复②—④步，直至加载完成、计算结束。

图3　耦合迭代流程图

4　耦合算例分析

4.1模型参数及边界条件如图4为一竖井段圆形隧洞轴对称模型，X、Y轴平行于隧洞断面，位于水平面内，竖直方向为Z轴，向上为正，衬砌内半径r1为4.2 m，外半径r2为5 m，围岩沿径向取30倍开挖洞径［23］，为300 m，模型厚10 m，衬砌与围岩均由八节点实体单元模拟，共计9865个节点和7888个单元，衬砌与围岩之间共节点。在衬砌内外两侧铺设直径为28 mm的钢筋，钢筋到隧洞中心距离R1、R2分别为4.3 m、4.9 m。模型材料参数如表2所示。

图4　有限元模型

隧洞埋深100 m，沿Z轴负向施加100 m的岩体上覆压力于模型上表面以模拟初始应力场的影响。在初始时刻，衬砌及围岩处于饱和状态，结构内部为饱和渗流，衬砌及围岩孔隙率设为1.5%、10%。采用增量水头法计算渗流场及应力场［24］，忽略重力对增量水头影响范围即水力半径的影响，针对增量水头的影响范围将围岩外圈设为零水头边界［25］。模型顶面以及底面设定为不透水边界，衬砌内壁设为定水头边界，计算过程中增大衬砌内壁的定水头边界得到增量水头以模拟充水加压过程。计算渗流场时，模型施加全位移约束；计算应力场时，模型底部和外圈节点施加法向位移约束。不考虑模型自身以及洞内水体的重力，只考虑衬砌及围岩中水体的重力。内压水头从0.5 MPa逐级加载至1.1 MPa，每级增量水头控制在0.025～0.1 MPa，保证每一级内压水头下耦合迭代完全。

表2　材料参数

4.2衬砌损伤开裂特征通过试算，内压小于0.5 MPa时，衬砌不会出现损伤。内压增至0.65 MPa时，损伤分布均匀，从内向外损伤值由大变小，最大值为0.108，详见图5（a）。在0.67 MPa下迭代4次达到稳定，衬砌环向出现3条贯穿缝，等间距分布于0°（隧洞右腰）、120°和240°位置（逆时针），裂缝分布少而宽，详见图5（b）。此后直至加载完毕，不再出现新的裂缝，裂缝处损伤值随内压持续增大至0.999，且裂缝间混凝土中部出现拉损伤，原因在于裂缝出现后，衬砌结构类似于弹性地基（曲）梁，使得裂缝间衬砌外侧出现反弯点，拉应力增大，出现明显的拉损伤区［26］，详见图5（c）。通过耦合迭代计算得到的衬砌损伤开裂特征与隧洞压水试验结果（表1）吻合较好［7］，说明本文建立的耦合模型能够较好地反映衬砌的开裂过程。

图5　衬砌损伤开裂情况

4.3渗透孔压演化特征为避开模型顶面以及底面边界条件的影响，选取模型中部位置对渗透孔压演化特征进行分析。内压为0.5 MPa时，衬砌未进入塑性状态，无损伤，衬砌承担了大部分的水压力，围岩内孔压相对较小，不超过0.115 MPa，详见图6（a）。在内压达到0.67 MPa，衬砌开裂前，衬砌内侧出现损伤，损伤区域孔压较大，围岩内孔压有所增大，最大值为0.257 MPa，详见图6（b）；3条贯穿缝出现后，内水外渗，渗流场出现明显的重分布现象，围岩内孔压升高，在衬砌环向0°、120°及240°裂缝位置渗流场变化明显，孔压梯度明显降低，裂缝处围岩承受的水压接近内压水头，孔压自裂缝位置向围岩深部逐渐消散，渗流场分布与衬砌损伤演化情况相对应，基本满足等间距分布，详见图 6（c）。

图6　隧洞渗流场分布（单位：Pa）

随着内压的增大，衬砌损伤峰值与渗透系数峰值同步变化，表示材料损伤与渗透系数实现了同步更新。裂缝处混凝土损伤与渗透系数在衬砌开裂时出现陡增现象，此后随内压持续增大，加载完成时，渗透系数相较于初始1×10-9m/s提高了约3个量级，裂缝成为内水外渗的主要通道，详见图7。

4.4钢筋应力演化特征充水加压过程中，裂缝处钢筋应力随内压增大，而未开裂位置钢筋全程处于低应力状态，应力值小于14 MPa。当内压小于0.67 MPa时，钢筋应力缓慢增长；贯穿缝形成后，裂缝位置混凝土应力释放，大部分荷载转由钢筋承担，裂缝处钢筋应力突增至50.120 MPa；当内压从0.67 MPa增至0.9 MPa时，裂缝处钢筋应力缓慢增至62.770 MPa，当内压达到1.0 MPa时，裂缝处钢筋应力进一步增大，在内压1.1 MPa迭代稳定后裂缝处钢筋应力达到228.800 MPa，详见图8。原因在于贯穿缝出现后，裂缝处的钢筋所发挥的作用类似于锚栓，单独承担拉力，而未开裂位置钢筋与衬砌混凝土共同承担衬砌环向拉力，混凝土承担大部分拉力，钢筋所承担的拉力较小。

图7　衬砌损伤峰值及渗透系数峰值随内压的演化

图8　钢筋应力演化特征

5　衬砌与围岩有条件联合承载机理研究

衬砌开裂前，透水性远弱于围岩，内部水力梯度较大，承担大部分的水荷载，而围岩内部水力梯度较小，承担的水荷载较小，因此衬砌沿径向朝外变形的趋势更大，衬砌与围岩处于压紧状态；衬砌开裂后，透水性显著增强，内部水力梯度大幅降低，所承担的水荷载明显减小，此时围岩的透水性相对较弱，内部水力梯度相对较大，承担大部分的水荷载，因此围岩沿径向朝外变形的趋势更为显著，围岩拉着衬砌向外变形，同时衬砌受到内置钢筋的约束，在衬砌与围岩交界面上产生拉应力，当拉应力超过交界面粘结强度时，衬砌与围岩将出现脱离现象，表现出有条件联合承载特性［10，27-29］。

在原有模型的基础上，对计算模型进行局部修改，考虑衬砌与围岩之间不共节点，在衬砌与围岩之间建立面-面接触单元以模拟衬砌与围岩有条件联合承载特性，接触单元切向摩擦系数取为0.5，局部衬砌-接触单元-围岩模型如图9所示。此时衬砌与围岩接触面只存在法向力与切向摩擦力，当法向力为压力时，两者处于压紧状态，一旦法向受拉，两者脱离，不再传递法向应力，也不再满足位移连续条件，即默认衬砌与围岩交界面粘结强度为零，接触压力与间隙的关系如图10所示。在应力场计算中，接触面的作用主要是传递法向应力；在渗流场计算中，接触可以保证孔隙流体在不同材料之间的渗透作用。

5.1衬砌与围岩的接触状态内压增加到0.625 MPa前，衬砌未开裂，衬砌渗透特性远弱于围岩，衬砌承担了大部分水荷载，在0.625 MPa内压下迭代前，10 m厚度范围内衬砌内壁内水外渗总流量仅为1.109 L/min，此时衬砌与围岩不存在脱离现象，详见图11（a）；在0.625 MPa内压下迭代完成后，衬砌环向在290°位置（衬砌右腰为0°位置，逆时针旋转）形成一条径向贯穿裂缝，同时衬砌除贯穿缝外的其他部位的损伤值在迭代过程中有所增大，渗透性增强，两部分因素的作用下使得衬砌内壁内水外渗总流量陡增至6.850 L/min。此后随着内压增大，不再出现新的裂缝，衬砌损伤开裂特征与隧洞压水试验一致（表1）。由于裂缝处水力梯度的降低以及衬砌向外变形受到内置钢筋的约束，衬砌与围岩首先在290°开裂位置出现脱离现象，产生的空隙由渗透水充满，脱离范围在迭代过程中逐步扩大，迭代11次后达到稳定状态，脱离范围达75%，详见图11（b）。随着内压增至0.635 MPa，迭代4次后达到稳定状态，衬砌与围岩已经完全脱离，详见图11（c）。

图9　衬砌-接触单元-围岩模型

图10　接触压力与间隙的关系

图11　衬砌与围岩接触状态（深色为脱离部分）

5.2钢筋应力演化特征引入接触单元后，衬砌钢筋应力表现出明显的差异。衬砌开裂前钢筋应力较小，内压达到0.625 MPa时，衬砌开裂，钢筋应力呈整体上升趋势，且存在较大波动，迭代过程中裂缝处钢筋应力最大值达到327.500 MPa，此后随迭代过程钢筋应力逐渐趋于稳定，并在水压力缓慢增加的过程中，钢筋应力基本稳定并呈缓慢降低的趋势。在1.1 MPa内压作用下，钢筋应力最大值为56.560 MPa，整体量值较小，详见图12。出现此现象的原因在于衬砌开裂后，衬砌与围岩逐步脱离，围岩成为水荷载的主要承担对象，而衬砌承担的水荷载逐渐减小，导致钢筋应力量值较小。同时与完全联合承载情况不同，完全联合承载时衬砌在开裂之后仍然会随着围岩向外扩张变形，导致钢筋应力持续增大。相比之下衬砌与围岩相互脱离，即有条件联合承载极有可能是隧洞发电运行过程中钢筋应力普遍小于设计预期的原因，并且钢筋应力的这一演化特征和量值很好地解释了相关工程的监测结果［14-15，30］。

图12　有条件联合承载时钢筋应力演化特征

6　结语

本文基于混凝土塑性损伤模型与透水衬砌理论，通过等效节点荷载实现了渗流场对应力场的作用，建立了损伤开裂与渗流之间的耦合关系，以此提出了透水衬砌渗流-应力-损伤耦合算法。基于ABAQUS平台进行二次开发，通过Fortran语言调用实用程序GETVRM获取材料点的损伤值，利用子程序USDFLD实现了材料渗透系数随损伤的动态更新，完成了隧洞初次充水期间衬砌损伤开裂的耦合迭代过程。研究了衬砌损伤开裂、渗透孔压以及钢筋应力的演化特征，并采用面-面接触单元探讨了衬砌与围岩有条件联合承载机理，可以发现：（1）本文提出的耦合算法所得到的衬砌损伤开裂特征与衬砌压水试验结果吻合，表现为衬砌在一定内压下开裂，此后不再出现新的裂缝，只在裂缝间衬砌外侧出现拉损伤，较好地反映了水工隧洞初次充水加压工况中衬砌损伤开裂的过程。（2）衬砌开裂前后，渗流场表现出明显的差异。开裂后，裂缝处孔压梯度明显降低，渗透系数呈现数量级的增大，裂缝成为内水外渗的主要通道，这与工程普遍认知是相符的，说明本文所提出的耦合算法是可行的，具有一定的推广性。（3）完全联合承载时，衬砌与围岩之间径向变形协调，裂缝处钢筋应力持续增大，未开裂区钢筋全程维持在低应力水平；考虑有条件联合承载特性后，衬砌开裂之后与围岩逐步脱离，围岩成为承担水荷载的主体，衬砌钢筋应力出现应力回缩现象，在加载完毕时应力量值维持在较小水平，揭示了有条件联合承载是隧洞发电运行过程中钢筋应力普遍小于设计预期的原因，解释了相关工程的监测结果，说明本文所提出的耦合算法可以为隧洞衬砌设计问题提供一定参考。

参考文献：

［1］杨林德，丁文其.渗水高压引水隧洞衬砌的设计研究［J］.岩石力学与工程学报，1997，16（2）：17-22.

［2］杜小凯，任青文，陈伟 .有压引水隧洞内水外渗作用研究［J］.四川大学学报（工程科学版），2008，40（5）：63-68.

［3］卞康，肖明.水工隧洞衬砌水压致裂过程的渗流-损伤-应力耦合分析［J］.岩石力学与工程学报，2010，29（S2）：3769-3776.

［4］SIMANJUNTAK T D Y F，MARENCE M，MYNETT A E，et al.Mechanical-hydraulic interaction in the lining cracking process of pressure tunnels［J］.Int.J.Hydropow.Dams，2013，20（5）：98-105.

［5］BIAN K，XIAO M，CHEN J.Study on coupled seepage and stress fields in the concrete lining of the underground pipe with high water pressure［J］.Tunnelling and Underground Space Technology，2009，24：287-295.

［6］FAHIMIFAR A，MOHAMMAD R Z.A theoretical solution for analysis of tunnels below groundwater considering the hydraulic-mechanical coupling［J］.Tunnelling and Underground Space Technology，2009，24：634-646.

［7］张有天.水工隧洞建设的经验和教训（上）［J］.贵州水力发电，2001，15（4）：76-84.

［8］刘秀珍.关于圆形有压水工隧洞衬砌裂缝计算的探讨［J］.水力发电，1983，12：39-46.

［9］沈威.圆形水工压力隧洞衬砌变形特性与限裂设计研究［D］.大连：大连理工大学，2012.

［10］苏凯，伍鹤皋，周创兵 .内水压力下水工隧洞衬砌与围岩承载特性研究［J］.岩土力学，2010，31（8）：2407-2412.

［11］ZHOU Y，SU K，WU H.Hydro-mechanical interaction analysis of high pressure hydraulic tunnel［J］.Tunnelling and Underground Space Technology，2015，47：28-34.

［12］GRAZIANI A，BOLDINI D.Influence of hydro-mechanical coupling on tunnel response in clays［J］.Journal of Geotechnical&Geoenvironmental Engineering，2012，138（3）：415-418.

［13］YOO C.Interaction between tunneling and groundwater-numerical investigation using three dimensional stress-pore pressure coupled analysis［J］.Journal of Geotechnical&Geoenvironmental Engineering，2005，131（2）：240-250.

［14］叶冀升.广蓄电站水工高压隧洞设计施工的若干问题［J］.水力发电学报，1998（2）：39-50.

［15］李煊明.周宁水电站透水衬砌高压钢筋混凝土岔管研究［J］.水电站设计，2007，23（2）：15-17.

［16］吴梦喜.饱和—非饱和土中渗流Richards方程有限元算法［J］.水利学报，2009，40（10）：1274-1279.

［17］纪佑军，刘建军，程林松 .考虑流-固耦合的隧道开挖数值模拟［J］.岩土力学，2011，32（4）：1229-1233.

［18］LUBLINER J，OLIVER J，OLLER S，et al.A plastic-damage model for concrete［J］.International Journal of Solids&Structures，1989，25（3）：299-326.

［19］LEE J，FENVES G L.Plastic-damage model for cyclic loading of concrete structures［J］.Journal of Engineering Mechanics，1998，124（8）：892-900.

［20］苏凯，张伟，伍鹤皋，等.考虑摩擦接触特性的钢衬钢筋混凝土管道承载机理研究［J］.水利学报，2016，47（8）：1070-1078.

［21］贾善坡，陈卫忠，于洪丹，等.泥岩隧道施工过程中渗流场与应力场全耦合损伤模型研究［J］.岩土力学，2009，30（1）：19-26.

［22］周亚峰.水工隧洞渗流-应力-损伤/开裂耦合理论及方法研究［D］.武汉：武汉大学，2016.

［23］周亚峰，苏凯，伍鹤皋.水工隧洞钢筋混凝土衬砌外水压力取值方法研究［J］.岩土力学，2014，35（S2）：198-203.

［24］张有天，张武功.隧洞水荷载的静力计算［J］.水利学报，1980（3）：52-62.

［25］潘家铮，唐家骥，张有天，等.隧洞水荷载静力计算（学术讨论）［J］.水利学报，1981（5）：75-81.