基于扩展有限元的混凝土重力坝水力劈裂分析

石艳柯， 霍俊怡， 张多新， 石路杨

(华北水利水电大学土木与交通学院， 郑州 450045)

随着南水北调等大型引水、调水工程的建设，一大批水工建筑物在高水头、大埋深等恶劣的水文地质环境条件下服役，其安全性的评估显得十分重要。混凝土坝上游常年受到高水压及淤泥力等复杂的荷载作用，且在初始裂纹面上受到水压力的作用，致使大坝初始裂纹极易起裂，进而扩展，严重时会形成贯穿裂纹导致大坝承载力降低，影响到坝体整体性和安全性。水力劈裂是指由于水压的抬高，在岩(土)体中引起裂纹产生和扩展的一种物理现象[1]，宁文祥等[2]通过页岩油储层水力压裂裂缝形态的试验，在油气藏的勘探中提出了的水力压裂技术；王志荣等[3]针对水力压裂条件下煤层气井初始压裂缝转向问题，研究了水力压裂近井区裂缝转向扩展机理。混凝土坝因体量较大，裂纹在水力劈裂作用下的分析常采用数值方法。水力劈裂作用下裂纹扩展的数值模拟涉及非线性和移动边界问题，扩展有限元(extend finite element method， XFEM)在求解不连续问题上独有优势，相对于诸多水力劈裂模拟方法(如数值流形法、有限元法、边界元法、扩展有限元法和无网格法)能够较好地解决这一问题。中外诸多学者[4-6]采用不同的数值模拟方法对混凝土坝进行了开裂研究，对单裂纹和多裂纹的扩展路径及裂纹影响因素进行了对比分析，得到了混凝土坝的断裂破坏全过程。甘磊等[7-9]考虑不同初始裂纹位置及长度等因素，对混凝土坝裂纹扩展进行了详细的分析，结果表明缝内水压作用对重力坝建基面抗滑稳定不利。董玉文等[10]和高景泉等[11]通过预置混凝土坝踵位置处裂纹，对裂纹面上不同分布水压力载荷作用进行开裂分析，并对重力坝水力劈裂的部位、位移、应力分布特性进行了分析。靳旭等[12]主要考虑坝体单一裂纹的存在对重力坝应力分布的影响，其结果表明基岩弹性模量对裂纹扩展影响较小，并可适当地稀疏网格提高工作效率，但在实际工程中裂纹往往是多条存在且在坝踵位置处极易开裂，且水力压裂对裂纹扩展影响较大，在重力坝开裂分析中应给予考虑，但并未考虑。王新萍等[13]基于扩展有限元对煤层转向压裂进行数值模拟，对裂纹轨迹进行了精确的预测，并将数值模拟结果与试验结果进行比较，证明了该算法的准确性和可靠性。刘金龙[14]以Cohesive单元为基础建立储层裂缝水力压裂数值模型，研究了裂缝在储层中的扩展规律及其扩展过程中缝间的相互干扰作用，参考该方法可简便且精确地对人为控制裂缝的扩展。郭秉乔等[15-16]采用扩展有限元模型和黏聚力裂纹模型对含有多个初始裂纹的混凝土重力坝进行扩展分析，结果表明随着裂纹的逐渐扩展，混凝土坝整体的刚度和承载力减低。石路杨等[17-18]建立了地下隧洞水力劈裂的扩展有限元法模型，分析了水压力对地下隧洞的稳定性的影响和隧洞在内水压力和裂纹内水压的共同作用下的破坏过程。

有限元软件ABAQUS功能的日益完善，子程序XFEM已经可以解决多数工程中的断裂问题，已能够有效并较为精准地对裂纹扩展进行数值模拟。文献[13-14]借助ABAQUS软件实现了水力压裂裂纹数值分析，郭秉乔等[15-16]也借助该软件内的子程序对重力坝进行裂纹扩展分析，并得到较理想的结果。

重力坝上游面多数为直断面，选取的宁利碾压混凝土重力坝上游面为直断面与折断面的组合面，上述文献中对该结构形式的重力坝裂纹扩展分析较少，文献[19]裂纹扩展分析考虑该结构形式的黄登重力坝，但仅考虑坝体裂纹的扩展，未对坝踵位置处裂纹进行分析。现建立二维混凝土重力坝水力劈裂数值分析模型，基于扩展有限元基本原理和最大周向拉应力强度因子理论[20]，采用ABAQUS软件子程序XFEM对带初始裂纹的混凝土坝进行裂纹扩展分析。通过在坝踵预置不同长度和位置的初始裂纹，分析裂纹的开裂路径与长度、角度及裂纹面水压力之间的关系。此外，在坝体预置多条初始裂纹，对不同裂纹长度扩展的路径及裂纹尖端应力进行分析，并选取坝体关键点及裂纹相对张开口位移进行分析，为后续的该类型混凝土坝裂纹扩展分析提供一定的参考依据。

1 扩展有限元的基本理论

1.1 扩展有限元位移模式

扩展有限元法(XFEM)是在常规有限元分析框架下，在裂纹单元引入加强函数(例如表述裂纹贯穿和裂尖单元两侧位移不连续的加强函数)来表示裂纹的存在。在XFEM中，裂纹的位移场如式(1)所示，含裂纹的扩展有限元计算网格如图1所示。

(1)

(2)

图1 XFEM加强模式Fig.1 XFEM enrichment mode

式中：I、K、J、M分别为常规有限元节点集、裂纹贯穿单元节点集、裂纹尖端单元节点集和裂纹交叉单元节点集；NI、NK、NJ、NM分别为相应节点集的插值函数；uI为常规有限元节点的自由度，uK、uJ、uM为节点附加自由度；H(x)为裂纹贯穿单元的附加函数，来扩展有限元与常规有限元位移相协调，同时可以表征裂纹贯穿位置的位移不连续的性性质；φJ(x)为裂尖单元的附加函数，能够精确地表征裂尖位置的位移场；JM(x)为多裂纹交叉位置处单元的附加函数，以实现主次裂纹的汇合;r、θ为裂尖在局部极坐标系下的半径、夹角。式(1)中等号右侧第一项表示常规单元，第二项表示裂纹贯穿单元，第三项表示裂纹交叉单元，第四项表示裂纹尖端单元。

1.2 水力劈裂的控制方程

对于含有附加函数的富集单元，裂纹面上同一点的相对位移w(即裂纹开裂宽度)可根据式(1)求得

(3)

根据虚功原理，结合构造的位移模式，带裂纹结构水力劈裂问题的虚功方程[12]表示为

(4)

将式(1)代入虚功方程式(4)中，并结合式(3)即可得到水力劈裂问题的扩展有限元支配方程为

Kd=f

蓝宝石英文名称为Sapphire，源于拉丁文Spphins，意思是蓝色。象征着稳健、端庄和聪慧。自古以来，人们迷醉于蓝宝石内所蕴含的安静和强烈的力量。据说，蓝宝石是距离神灵最近的宝石，它湛蓝的颜色在任何时候都不会改变，甚至影响到蓝色天空的形成。

(5)

d={uI，uK，uJ}

(6)

(7)

(8)

(9)

2 计算模型及参数

宁利水库位于云南省丽江市宁蒗县，由枢纽工程和输水渠道组成，枢纽主要建筑物由拦河大坝、上游拦砂坝和导流隧洞灌溉干渠组成。大坝为碾压混凝土重力坝，最大坝84.0 m，坝轴线长252.0 m。依托于此工程，选取大坝标准断面建立二维平面模型进行裂纹扩展分析。坝体高84 m，坝顶宽8 m，坝体总宽71.76 m，上游迎水面由垂直面与坝坡比1：0.2的折断面组成，下游坝面坝坡比为1∶0.76，沿坝轴线分别取上游、下游取92 m和60 m，土体垂直向下取100 m，坝体断面如图2所示。

模型的边界条件及荷载条件为：在坝基左右边界施加水平向约束，坝基底部边界施加竖向约束，考虑的荷载有坝体与坝基的自重、上游库水压力及裂纹面上的水压力P，其中上游水压力以超载系数为2.5的梯度线荷载施加，因库水面与起始裂纹面的水头差一定，故将裂纹面上水压力视为均布荷载，以垂直于裂纹面的均布线荷载施加。

根据文献[11，15]所研究的结果，在坝踵处预置裂纹只会不断向基岩底部扩展，坝体内预置的裂纹只会向下游面扩展，并未形成贯穿坝体与坝基的裂纹，因此将坝基和坝体视为一个整体进行裂纹扩展分析，未考虑坝体与坝基的接触滑移对裂纹的影响。坝体与坝基均采用有限元软件ABAQUS中平面应变单元CPE4R模拟，裂纹通过预置的方式布置，并和坝体与坝基组装形成一个整体，使用最大周向拉应力强度因子理论[20]，借助于软件中的子程序XFEM对重力坝进行裂纹扩展分析，研究裂纹的数量、位置及裂纹面水压力对裂纹扩展的影响，方便进一步对重力坝进行安全评估。宁利重力坝为碾压混凝土重力坝，为标准的“金包银”，需要考虑不同材料的分层，但因预置初始裂纹长度较长贯穿外部材料层，裂尖深入内部材料层，故仍将坝体混凝土材料视为各向同性。因宁利重力坝仍然处于待建中，其碾压混凝土的断裂能参数并未及时获得，因此参照文献[19]中相似的重力坝材料参数，给出表1所示的坝体混凝土与坝基基岩的材料参数。

L为裂纹长度；θ为预置裂纹与水平夹角；H为坝顶；G为坝踵； J为坝趾图2 重力坝裂纹预置示意图Fig.2 Pre-crack diagram of the dam

表1 基岩和坝体混凝土的材料参数Table 1 Material parameters of the bedrock and the dam

3 数值结果分析

3.1 坝踵裂纹扩展分析

由于混凝土重力坝的独有结构形式，重力坝的坝踵位置处易产生大应力区域，据此在坝踵位置通过预置不同裂纹长度L(2、5、8、10 m)和角度θ(15°、30°、45°、60°)的初始裂纹进行开裂分析，裂纹预置位置如图2所示。

L=5 m，θ=30°时，考虑裂纹面不同水压力大小作用影响，裂纹的具体扩展路径如图3所示。当不考虑裂纹面水压力作用即P=0时，由于上游水体的超载作用和坝体自重的作用，导致初始裂纹逐步向坝基底部扩展，并随着裂纹面水压力的增加，裂纹开裂深度逐渐增加，P=104Pa和P=105Pa，裂纹进入稳态扩展，当裂纹面水压力进一步增加(P=106Pa)，裂纹失稳扩展，直至形成了贯穿裂纹，致使坝体承载力失效，影响到了坝体的安全度。

图4给出了初始裂纹与水平向夹角θ=30°且裂纹面水压力P=105Pa时不同初始裂纹长度L的扩展路径，当θ=30°时，裂纹起裂方向与初始裂纹夹角约为48°，裂纹的扩展方向与初始裂纹长度无关。

图5给出了初始L=5 m且P=105Pa时不同θ的扩展路径，随着θ的增加，裂纹起裂方向与初始裂纹的夹角逐渐降低，当θ=60°时，由于裂纹面的水压力与坝体自重、超载静水压力的等效作用力正好与初始裂纹的预置方向垂直，并依据最大主应力起裂准则，此时裂纹起裂方向与初始裂纹的夹角趋于0°即沿初始裂纹方向起裂；当θ=45°时，裂纹的扩展深度最低，且扩展长度最短，主要是此时裂纹面的水压力抵消了部分坝体自重影响。

图3 不同裂纹面水压力的裂纹扩展路径(初始L=5 m)Fig.3 Crack propagation path under different water pressure on the crack surface (initial L = 5 m)

图4 不同初始裂纹长度扩展路径(初始θ=30°)Fig.4 Crack propagation path under different initial crack length (initial θ= 30°)

图5 不同初始裂纹角度扩展路径(初始L=5 m)Fig.5 Crack propagation path under different initial crack angles (initial L = 5 m)

3.2 坝体裂纹扩展分析

图3～图5已经对待建混凝土重力坝进行了坝踵位置的单裂纹扩展分析，得到了坝踵位置处预置裂纹的不同长度、位置及裂纹面水压作用下扩展规律。为了进一步对坝体裂纹扩展规律分析，考虑自重、上游超载水压作用和裂纹面水压力作用，通过预置多条初始裂纹A、B、C，图6给出了不同裂纹长度多条初始裂纹下扩展的等效mises应力，表2为不同位置处关键点位移数值。

由图6可看出，随着预置裂纹B的长度增加，预置裂纹尖端mises应力值逐渐增大，裂纹B开裂深度加大，不断地向下游坝基面方向扩展。裂纹A与C尖端虽出现应力集中，但裂纹面驱动力(主要由裂纹面处的剪应力提供)不足以支撑裂纹进行扩展。

表2给出预置裂纹B不同长度下的坝体坝顶H、坝踵G、坝趾J的位移值及裂纹最大张开位移数值可以看出随着裂纹长度的增加，坝顶位移及裂纹张开度也相应增加，但坝踵G和坝趾J位置处的位移受预置裂纹B的长度影响较小。

图6 多裂纹扩展模式下等效应力云图Fig.6 Equivalent stress diagram of multiple crack propagation

表2 不同裂纹长度下坝体关键点位移对比Table 2 Displacement comparison of key points of the dam under different crack lengths

图7所示为B裂纹的张开位移变化趋势，当B裂纹初始长度为5 m时，其最大张开位移Δux、Δuy、Δusum数值相差较小，均小于1 mm。随着裂纹B的初始长度增加，裂纹最大张开位移也随着增加，竖向张开位移Δuy增幅较大，水平向最大张开位移Δux增幅较小，说明垂直于裂纹面的法向力主要控制水力劈裂裂纹的张开位移，裂纹面的剪切力对裂纹张开位移贡献较小，裂纹初始长度对裂纹张开口位移有着显著的影响。

图7 裂纹B的最大张开位移Fig.7 Maximum opening displacement of crack B

4 结论

基于扩展有限元XFEM基本原理，选取某待建混凝土重力坝，考虑不同裂纹长度、位置及裂纹面水压力进行坝踵及坝体裂纹的水力劈裂数值分析，为混凝土重力坝劈裂分析提供一定的参考依据。考虑重力坝自重、上游超载静水压力和裂纹面水压力荷载作用，对坝踵和坝体预置初始裂纹水力压裂分析，可以得到如下结论。

(1)坝踵预置裂纹长度L=5 m，角度θ=30°时，随着裂纹面水压力数值的增加，由裂缝面剪应力提供的驱动力也随之增大，裂纹不断向坝基底部边界扩展，重力坝承载力逐渐降低，影响到结构整体安全性。

(2)坝踵位置处裂纹面水压力P=105Pa，角度θ=30°时，不同长度的预置裂纹起裂与预置裂纹夹角基本相同，且最终扩展深度也基本相同，此时裂纹面的剪应力大小主要由结构自重及超载静水压力提供；坝踵处预置与水平向不同夹角θ的初始长度L=5 m的裂纹，当裂纹夹角θ=45°时，裂纹扩展深度最浅，θ=60°时与初始裂纹夹角趋于0°，主要是复杂荷载在平行于裂纹方向的分量先增大后减小。

(3)随着折断面处初始裂纹长度的增加，裂纹尖端应力值逐渐增大，且裂纹面驱动力不足以支撑裂纹A、C扩展，裂纹B的扩展深度随之加深，不断向下游坝基面扩展。

(4)坝体预置裂纹B随着裂纹长度的增加，坝顶位移及裂纹张开度也相应增加，但坝踵和坝趾位置处的位移受预置裂纹B的长度影响较小。