衬砌水力劈裂已成为深埋压力隧洞发生事故的重大隐患，对水电建设和水资源开发构成了重大威胁，压力隧洞衬砌水力劈裂模拟已引起越来越多国内外学者的关注。众多学者从数学与力学角度出发，通过引入断裂力学、损伤力学等理论，建立了压力隧洞渗流–应力耦合分析模型。张巍等^[1]考虑渗流场与应力场、损伤场的关系，建立了应力–损伤–渗透系数方程，采用3维弹塑性有限元程序，对大型洞室围岩开展应力–损伤–渗流耦合分析。Zhou等^[2]认为压力隧洞衬砌与围岩为有条件联合承载，同时假定钢筋与衬砌混凝土非协调变形，采用不均匀系数计算开裂衬砌的二次渗透系数。陈卫忠等^[3]建立了弹塑性耦合损伤本构模型，并考虑渗透系数随围岩变形的动态演化，对锦屏Ⅱ级水电站引水隧洞运行期衬砌与围岩开展受力变形研究。胡云进等^[4]基于传统黏结单元模型，提出带有渗流分析功能的黏结单元模型，采用比奥固结理论，建立了压力隧洞渗流–应力–开裂耦合模型。Dadashi等^[5]分析了内水压力对衬砌渗透系数的影响，提出渗流场与应力场直接耦合模型，表明钢筋混凝土隧洞内水外渗不可避免。此外，杜小凯^[6]、Sousa^[7]、Zareifard^[8]等也对压力隧洞水力劈裂进行了研究，取得了一定成果。上述研究在建立压力隧洞水力劈裂分析模型时均假设隧洞衬砌为均质体，忽略了混凝土衬砌存在的初始缺陷，无法模拟出压力隧洞投入运行后衬砌的随机开裂位置和随机开裂路径，与实际不符。混凝土衬砌在施工、养护过程中由于施工质量、温度作用等原因，在投入运行前不可避免地存在一些随机分布的初始缺陷，而这些初始缺陷决定了衬砌的水力劈裂位置和劈裂路径^[9]。因此，要准确模拟内水压作用下压力隧洞衬砌的随机水力劈裂过程，必须合理考虑混凝土衬砌的随机初始缺陷。

针对目前研究的不足，假定混凝土衬砌为两相介质，将施工、温度等原因导致的初始缺陷视为弱相，采用孔压黏结单元模拟，并按Weibull分布对黏结单元力学参数进行随机赋值，其余部分视为强相，采用实体单元模拟；按体力理论施加内水压力，建立压力隧洞衬砌随机水力劈裂分析模型。最后，与某大尺度压力隧洞结构模型试验结果^[10]进行比较，验证了分析模型的合理可靠性。

压力隧洞混凝土衬砌由于施工质量、温度作用等原因，在投入运行前其内部存在大量随机分布的初始缺陷。初始缺陷相较于骨料与砂浆力学性能明显较差，是混凝土衬砌内部的薄弱环节。根据最弱环原理，在内水压作用下率先在这些薄弱环节处发生破坏，并失稳扩展形成混凝土衬砌贯穿裂缝^[11]。因此，将混凝土衬砌中的骨料与砂浆假定为强相，初始缺陷作为弱项，弱相基质分布于强项基质周围，将二者分开考虑。

衬砌开裂主要发生在初始缺陷等弱相基质中，而压力隧洞衬砌水力劈裂其实质是渗流–应力耦合问题。为合理描述混凝土衬砌开裂力学特性及渗流–劈裂耦合关系，基于传统渗流–应力耦合分析模型引入具有渗流分析能力的孔压黏结单元（图1）模拟混凝土衬砌中的弱相，孔压黏结单元在荷载作用下主要沿厚度方向发生法向变形。考虑到实际工程中压力隧洞衬砌主要沿径向发生开裂，故为方便设置黏结单元，在围岩与衬砌的实体单元间沿径向插入孔压黏结单元（图2）。图2中，1～3等为实体单元，9、10等为孔压黏结单元。衬砌与围岩连接处通过绑定约束模拟，确保连接处各节点各自由度变量协调。

图1(Fig. 1)

图1 孔压黏结单元 Fig. 1 Cohesive element of pore pressure

图2(Fig. 2)

图2 黏结单元径向插入 Fig. 2 Radial insertion of cohesive element

孔压黏结单元采用线弹性损伤力学模型（式（1）～（6）），应力–开度关系如图3所示。传统孔压黏结单元默认单元在受拉或受压时采用相同的初始法向刚度，对于压力隧洞水力劈裂分析，由于初始几何厚度为0，容易导致收敛性问题。Bandis等^[12]采用双曲线模型描述岩石节理的法向刚度与位移间的关系。该模型随着单元受压厚度趋近于0时法向刚度接近无穷大，以有效避免接触面两侧单元的嵌入。基于此思想，在传统孔压黏结单元^[4]基础上做改进，提出新型孔压黏结单元。假定孔压单元不考虑法向变形与剪切变形的耦合效应，发生张拉变形时仍采用与传统孔压黏结单元一致的本构关系，但在孔压单元受压出现压应变时，其法向刚度采用一个大值，以抵制黏结单元两侧单元嵌入，且受压刚度不考虑损伤效应，其余参数不变。在兼顾计算精度与收敛性的前提下，经过多次验算，取黏结单元受压时法向刚度值为受拉初始刚度的100倍（式（3））。黏结单元受拉达到抗拉强度之前处于弹性状态（式（1）），张拉应力达到抗拉强度以后考虑材料损伤对力学性质的影响，以损伤指标D^[13]（式（2））描述黏结单元损伤性状。退化后刚度由式（3）和（4）计算，相应的应力根据式（5）与（6）求得。

图3(Fig. 3)

图3 黏结单元应力–开度关系 Fig. 3 Relationship of stress–opening of cohesive element

式中， ${T_{\rm{n}}}$ 、 ${T_{\rm{s}}}$ 、 ${T_{\rm{t}}}$ 分别为法向与切向名义应力， $T_{\rm{n}}^{\rm{0}}$ 、 $T_{\rm{s}}^{\rm{0}}$ 、 $T_{\rm{t}}^{\rm{0}}$ 分别为法向与切向的名义应力峰值， $\mathop \delta \nolimits_{{\rm{m}},{\rm{max}}} $ 为加载过程最大有效相对位移， $\mathop \delta \nolimits_{{\rm{mf}}} $ 为黏结单元完全破坏时的有效相对位移， ${\delta _{{\rm{m0}}}}$ 为黏结单元起裂时有效相对位移， ${K_{{\rm{nn}}}}$ 、 ${K_{{\rm{ss}}}}$ 分别为法向和切向初始刚度， ${\delta _{\rm{n}}}$ 、 ${\delta _{\rm{s}}}$ 分别为法向与切向张开位移。

孔压黏结单元含6个节点（图1），其中1～4节点同时包含位移自由度与孔压自由度，5、6节点仅含孔压自由度。孔压单元渗流分析的同时考虑黏结单元内水体流动及其与周围实体单元间的水交换（图4），即孔压黏结单元中同时发生切向流（式（7））与法向流（式（10））。孔压单元切向流服从立方定律（式（8）），考虑单元开度对单元渗透性能影响（式（9））；法向流通过式（10）描述，以法向滤失系数近似模拟。

图4(Fig. 4)

图4 黏结单元中的渗流 Fig. 4 Seepage in the cohesive element

式中， $q$ 为单宽流量， $d$ 为黏结单元裂缝开度， ${k_{\rm{i}}}$ 为渗透参数， $\nabla p$ 为黏结单元切向孔隙压力梯度， $\mu $ 为流体动力黏滞系数， ${t_{{\rm{curr}}}}$ 为黏结单元当前几何厚度， ${t_{{\rm{orig}}}}$ 为黏结单元初始几何厚度， ${g_{{\rm{init}}}}$ 为黏结单元初始间隙， ${q_{\rm{t}}}$ 、 ${q_{\rm{b}}}$ 为黏结单元法向渗流量， ${c_{\rm{t}}}$ 、 ${c_{\rm{b}}}$ 为黏结单元法向滤失系数， ${p_{\rm{t}}}$ 、 ${p_{\rm{i}}}$ 、 ${p_{\rm{b}}}$ 分别为黏结单元顶面、中面（5、6两节点所在面）、底面的孔隙压力。

衬砌开裂最先发生于初始缺陷等弱相基质中。初始缺陷的随机分布特征决定混凝土衬砌的开裂位置与开裂路径。因此，假定强相（实体单元）的各项力学参数取为混凝土力学参数的各项均值，重点考虑衬砌内部弱相的随机特征对混凝土衬砌开裂的影响。

假定强相实体单元周围均匀分布弱相基质，不考虑初始缺陷的具体几何形态，通过假定弱相力学参数服从某种概率统计分布以描述弱相在衬砌中的力学随机性。假定材料为各向同性，根据需要对弱相黏结单元的断裂能、抗拉强度或弹性模量等参数进行随机赋值。由于混凝土材料力学参数的非负性，假定弱项材料力学性质满足Weibull分布^[14]，其分布密度函数通过式（11）定义：

式中： $u$ 为满足该分布参数（如强度、断裂能、弹模等）的变量； ${u_{\rm{0}}}$ 为与所有单元参数平均值有关的参数，但其数值并不是该参数的平均值；形状参数 $m$ 为均质度，定义了Weibull分布密度函数的形状，反映了参数的离散程度，其值越大，表明材料参数越均匀。

相应的随机变量 $u$ 的数学期望与方差根据式（12）～（13）计算：

式中， $\Gamma \left( {1 + 1/m} \right)$ 等为 $\Gamma $ 函数，可通过查表获得。

考虑到非高斯场直接生成算法的复杂性，通过ZMNL (零记忆非线性转换)变换^[13]，基于高斯场生成非高斯场（图5）。即假定高斯场X与非高斯场Y拥有相同的累积概率分布函数，由高斯场的分布函数求逆获得非高斯场样本：

图5(Fig. 5)

图5 高斯场向非高斯场转换示意图 Fig. 5 Diagram of transformation from gaussian field to non-gaussian field

式中， $W$ 为非高斯场分布函数， $\phi $ 为高斯场分布函数。

渗流–应力–劈裂耦合模型基于Biot固结理论^[15]建立模型力学平衡方程与渗流方程，并联立求解：

式中， ${\nabla ^2}$ 为拉普拉斯算子， ${\varepsilon _{\rm{v}}}$ 为固体体积应变， $\lambda $ 、 $G$ 为拉梅系数， $u$ 、 $v$ 为 $x$ 、 $y$ 方向上的位移分量， ${k_x}$ 、 ${k_{\rm{y}}}$ 为 $x$ 、 $y$ 方向上的渗透系数， $p$ 为孔隙压力， $\gamma $ 为固体颗粒重度， $\;{\rho _{\rm{w}}}$ 为水的密度， $\alpha $ 、 $Q$ 为Biot固结系数， $g$ 为重力加速度， $t$ 为时间。

压力隧洞衬砌与围岩采取不同的本构模型。隧洞衬砌由于弱相采用孔压黏结单元模拟，已考虑衬砌开裂。因此，强相采用实体单元模拟时，可采用较简单的线弹性模型。数值分析中假设钢筋跟衬砌混凝土位移协调，采用埋藏式钢筋模型，通过钢筋与混凝土的单元刚度矩阵组合获得单元综合刚度矩阵。隧洞围岩黏结单元采用线弹性损伤本构模型，抗拉强度取围岩均值。围岩实体单元本构模型选用线性扩展Drucker–Prager模型模拟，采用增量变刚度法进行迭代计算。

压力隧洞衬砌水力劈裂蒙特卡洛分析具体流程为：

1）ABAQUS生成初始计算模型，输出初始inp文件。

2）MATLAB读取初始inp文件，插入黏结裂缝单元。

3）利用MATLAB对模型随机赋值，调用ABAQUS进行数值计算。

a. 令 $i=1$ 。

b. MATLAB生成随机力学参数，对黏结单元随机赋值。

c. 对ABAQUS提交新的inp文件，进行第 $i$ 次数值计算， $i=i+1$ 。

d. 若 $i>M$ ，跳出循环；若 $i \le M$ ，跳至b。

4）MATLAB提取计算结果，进行蒙特卡洛分析。

以某大比尺压力隧洞结构模型试验为对象进行模拟分析。模型试验整体轮廓长6.5 m，宽5.5 m，高5.7 m。隧洞洞径0.8 m，衬砌厚度0.06 m，其中衬砌外设灌浆圈0.5 m。模型外部整体用钢板包裹。围岩按C10混凝土模拟，衬砌采用C20，环向钢筋采用 $\phi 8$ ，钢筋保护层厚度取7 mm。参照设计院提供的相关参数及文献[10]的数据，材料参数取值见表1。

表1(Tab. 1)

表1 数值模型材料参数 Tab. 1 Material parameters of numerical model

表1 数值模型材料参数 Tab. 1 Material parameters of numerical model 材料 弹性模量/ GPa 泊松比 $\mu $ 断裂能/ (N·m–1) 抗拉强度/ MPa 初始渗透系数/ (m·s–1) 混凝土 C10 175 0.168 100 1.00 5.27×10–8 C20 255 0.168 100 1.54 5.00×10–10 钢筋 2 000 0.300 — 235.00 —

根据文献[16]可知，混凝土骨料断裂能与水泥砂浆的比值约为3；文献[17]指出骨料–砂浆界面强度与弹性模量约为砂浆的40%～80%。因此，混凝土衬砌强相均取为材料均值，弱项材料参数仅考虑断裂能为随机值，均值取为25 N/m。本文根据模型施工质量取均质度为10。为提高计算收敛性，计算过程中黏结单元的弹性模量取100 GPa。利用ABAQUS软件进行隐式分析。针对大型物理模拟试验条件，考虑混凝土硬化后发生的体积收缩效应，顶部钢板与浇筑混凝土将产生间隙，因此，模型顶部边界设为自由面边界，底部边界设为固定边界，左右两侧采用水平方向约束。考虑到模型四周有钢板，模型四周视为不透水边界。模型边界及网格划分如图6所示，计算模型中衬砌与围岩通过绑定约束实现接触位置的位移与孔压协调。模型共划分10 220个单元，其中衬砌部分实体单元960个，孔压黏结单元960个，围岩与灌浆圈部分实体单元5 440个，孔压黏结单元2 720个，钢筋划分为140个truss单元。模型忽略自重影响，仅考虑内水压力，根据模型试验加载过程，内水压力随时间线性增加至1.8 MPa。黏结单元断裂能随机赋值采用Monte–Carlo方法。

图6(Fig. 6)

图6 数值模型与边界条件 Fig. 6 Numerical model and boundary

随机选取一个样本计算结果，分析压力隧洞模型在高内水压作用下的破坏过程。内水压力变化导致的压力隧洞孔压分布与裂缝开度云图如图7所示。图7表明：压力隧洞在受力初期（0.36 MPa）衬砌未发生水力开裂，内水压力主要由衬砌承担；压力水头沿衬砌径向衰减明显，沿环向均匀分布。随着内水压力增大，衬砌内部逐渐出现非均匀分布的微裂缝，当内水压达到1 MPa时，衬砌部分孔压黏结单元张开，表现为沿环向多处位置出现径向微裂缝区，但衬砌尚未出现贯穿裂缝。由水压分布云图可知，此时内水压力仍主要作用于衬砌。随着内水压继续增大（1.2 MPa），衬砌在45°、200°、225°、315°这4个位置发生贯穿裂缝。当内水压力为1.8 MPa时，225°位置处裂缝张开度最大达到0.037 7 mm，此时压力隧洞衬砌失去阻水作用，内水压力大部分由隧洞围岩承担。在加载后期，衬砌内部随机分布的微裂缝大部分未继续扩展甚至部分裂缝发生闭合，仅4条主裂缝在内水压作用下继续扩展，这与实际工程中压力隧洞衬砌开裂情况相符。分析其原因，主要由于某条主裂缝形成贯通，导致该主裂缝附近单元应力释放，向裂缝两侧发生回缩，一定范围内钢筋与混凝土单元间黏结作用力减小，导致混凝土环向应力减小，主裂缝附近微裂纹不增长或发生闭合；随着与裂缝间距离的增加，钢筋混凝土衬砌回复到开裂前应力状态。若此处应力大于混凝土抗拉强度，将出现第2条裂缝。此外，主裂缝的贯通导致高压水外渗，衬砌外侧水压力增大，内外水头差减小，混凝土衬砌向外发生位移的趋势得以减缓，且主裂缝位置开裂面受水压作用，有向裂缝两侧变形的趋势。图8为模型试验与数值模拟结果的裂缝分布比较，二者都在理论上易在裂缝的90°、270°等位置之外出现裂缝，验证了沿衬砌径向插入黏结单元并给黏结单元力学参数随机赋值可模拟衬砌随机水力劈裂，同时也说明数值模型考虑衬砌力学性质分布随机的必要性。

图7(Fig. 7)

图7 水压分布与裂缝张开度 Fig. 7 Water pressure distribution and fracture opening

图8(Fig. 8)

图8 衬砌裂缝分布 Fig. 8 Crack distribution of lining

图9为压力隧洞衬砌中某个钢筋单元的应力与内水压关系曲线及其均匀化曲线。由9图可知，钢筋应力变化可分为4个阶段。加载初期（0～1 MPa）由于衬砌未发生开裂，钢筋与内水压力表现为线性关系。此时，内水压力作为体力施加于钢筋混凝土衬砌，与以面力形式施加效果类似。随着内水压增加（0.8～1.2 MPa），钢筋与内水压力关系不再具有线性关系，主要由于衬砌部分单元发生破坏，隧洞衬砌水压分布发生变化。当内水压力达到1.2 MPa时，衬砌发生开裂，相较于衬砌完好状态，此时损伤开裂区钢筋应力出现急剧增长。当内水压达到1.7 MPa时，钢筋应力达到最大值约50 MPa。在加载后期（1.7 MPa以后），钢筋应力出现回缩，这主要是由于衬砌贯穿裂缝发生内水外渗，内外水压差减小，衬砌上渗透梯度降低，内水压力更多的由围岩承担。通过分析计算，在加载后期围岩承载达到内水压力的近60%，说明内水外渗改善了衬砌的受力状况。图10给出了钢筋单元应力均值曲线与物理模型试验最大裂缝处的钢筋应力实测结果曲线。分析图10可知两者变化趋势有较好的一致性，再次验证了本文模型模拟压力隧洞水力劈裂的可靠性。

图9(Fig. 9)

图9 样本内水压–钢筋应力关系 Fig. 9 Relationship of water pressure–reinforcement stress in sample

图10(Fig. 10)

图10 模型试验与数值模拟的钢筋应力值对比 Fig. 10 Comparison of reinforcement stress between model test and numerical simulation

根据数值结果分析，钢筋混凝土衬砌受力响应与普通钢筋混凝土结构完全不同。围岩与钢筋混凝土衬砌共同承受隧洞内水压力。在内水压作用下，衬砌混凝土弱相介质发生水力开裂导致内水外渗，衬砌外侧水压增加，衬砌水力梯度降低，内水压力被转移至围岩。相较于普通钢筋混凝土结构，混凝土衬砌裂缝分布呈非均匀且稀疏特性。因此，在压力隧洞内水外渗分析中须以体力形式考虑内水压力。此外，在高内水压作用下衬砌损伤对衬砌渗透系数影响明显，分析中有必要考虑渗流–应力–损伤耦合效应。由分析可知，混凝土衬砌不可避免地发生开裂，在衬砌开裂状态下，钢筋应力远未达到屈服应力。因此，在压力隧洞衬砌设计时，不必过分强调限裂设计，应将设计重点放在围岩及衬砌与围岩连结处。衬砌施工与保养应尽量避免温度效应、人为因素等导致的初始缺陷，以减少贯通裂缝数量与隧洞内水外渗量。压力隧洞衬砌开裂分析应考虑衬砌初始缺陷的随机性，通过对施工质量的调查合理选取随机参数。

将压力隧洞混凝土衬砌视为两相介质，施工、温度等原因引起的初始缺陷作为弱相，采用孔压黏结单元模拟，并按Weibull分布对黏结单元力学参数进行随机赋值；其余部分视为强相，采用实体单元模拟，实体单元力学参数选用混凝土相关参数；内水压力按体力施加，建立压力隧洞衬砌随机水力劈裂分析模型。以某压力隧洞结构模型试验为例，进行了内水压作用下的隧洞衬砌水力劈裂分析。计算结果表明压力隧洞衬砌在高内水压作用下会沿隧洞环向产生随机分布的多条裂缝。随着内水压的增长，衬砌中仅有2～4条裂缝继续扩展形成主裂缝，其余裂缝并不扩展甚至发生闭合。随着衬砌裂缝贯通，内水压力转而由围岩承担，衬砌钢筋应力会减小。数值模拟结果与实际工程裂缝分布及模型试验结果具有较好的一致性。此外，数值模拟的钢筋应力均匀化曲线与模型试验中最大裂缝处的钢筋应力曲线变化趋势比较吻合。因此，本文所建立的压力隧洞衬砌随机水力劈裂分析模型是合理可靠的，为压力隧洞内水外渗耦合分析提供了新的思路。