一种模拟准脆性断裂问题的相场法能量退化法则

2022-07-27庄梦如朱潇潇

中国农村水利水电 2022年7期

张巍，庄梦如，朱潇潇，沈蔚，陈敏

（淮安市水利勘测设计研究院有限公司，江苏淮安 223005）

0 引言

工程材料的裂缝形成与扩展是材料失效和结构破坏较为普遍的形式［1］，尤其是对水利工程中的防渗材料，如混凝土防渗体、土质防渗体、防渗衬砌等，一旦出现裂缝，可能会威胁工程结构的安全。如美国Teton大坝［2］由于高水压引起的裂缝扩展，导致溃坝事故；澳大利亚Kölnbrein 拱坝［3］由于坝踵裂缝出现了严重的渗漏，严重影响了工程效益。混凝土、岩石、非饱和黏性土及水泥砂浆等水利工程中常用防渗材料的裂缝扩展多为准脆性断裂或脆性断裂，量化这些裂缝的断裂过程，对预防裂缝引起的材料失效具有重要的意义。

在材料准脆性断裂模拟方面，相场法基于断裂能量变分原理，裂缝沿着使固体系统总能量最小的路径扩展，使得其可以较好地捕捉裂缝的开裂扩展。Miehe［4］等人给出了材料脆性断裂的相场法热力学框架，并成功应用于裂纹的静态、动态及多场耦合问题，此后相场法模拟准脆性断裂问题得到了广泛重视与蓬勃发展［1］。在相场法基本理论中，相场变量是区分材料的断裂区和未损伤区域的重要指标，为了模拟断裂区材料的力学性能变化，需要采用能量退化法则对相场变量与材料的力学特性进行耦合，材料的力学响应主要由能量退化法则控制［5］。因此，能量退化法则对相场法模拟材料的断裂行为十分重要。

Bourdin［6］提出了一个二次幂函数型式的能量退化法则，由于其函数形式简单、无需参数、易于编程实现，而被广泛使用，后发展为了一种单参数幂函数型式的能量退化法则［7］。为了进一步提高相场法在预测与裂纹形核以及现有裂纹扩展相关的临界载荷时的精度。Sargado［8］等运用加权思想，采用幂函数型式的能量退化法则对指数型式的能量退化法则进行修正。Wang［9］等基于广义线性软化法则，提出了一种适用于黏结模型的自动校准的能量退化法则。Lu［10］等从细观损伤角度出发，运用加权平均思想，建立了一种指数型式的材料宏观能量退化法则。Borden［11］提出了一种单参数多次幂函数型式的能量退化法则，其可较好地捕捉材料在断裂破坏前的线弹性行为。吴建营［1］在研究适用于脆性断裂和准脆性破坏的相场正则化内聚裂缝模型时，提出了基于多项式和幂函数型式的能量退化法则。Zivkovic［12］采用等效塑性应变，修正了二次幂函数能量退化法则的次数，形成了考虑塑性变形的能量退化法则。

在水利工程及岩土工程领域，现多采用二次幂函数型式的能量退化法则。刘国威［13，14］等研究了岩石的双平行翼型裂缝相交和动态裂缝三维曲面扩展等问题以及动力水力压裂作用下裂缝动态扩展。李鹏飞［15］等分析了包含不同岩桥倾角的预制双裂隙岩石类材料在单轴压缩作用下的损伤和破坏过程。易良平［16］等研究了煤砂互层中水力裂缝纵向延伸影响因素。Zhou［17］等分析了类岩石材料的压剪耦合脆性断裂，研究了不同地应力状态下岩石的水力压裂行为［18］。Bilgen［19］等模拟了岩石材料的巴西劈拉试验。Zhang［20］等研究了不同试样半径及不同初始预制裂缝宽度对直裂缝半圆弯曲岩石试样（NSCB）的裂缝扩展路径及断裂韧度的影响。Zhang［21］等在考虑拉压不等损伤的基础上，分析了单裂隙人造岩石的初始裂缝扩展路径与位移-荷载曲线。Reinoso［22］等结合黏结模型，模拟了岩石的巴西劈拉试验及含裂缝岩石的拉伸与压缩断裂行为。此外，侯越［23］等采用了三次多项式型式的刚度法则，分析了砂浆裂纹相互作用失效的行为。Bilgen［24］等采用Borden 能量退化法则，研究了不同参数对岩石张拉断裂行为的影响。Fei［25］等提出了拉伸-剪切双相场模型，采用统一相场黏结模型［26］，研究了岩石裂缝的张拉-剪切断裂耦合行为。Wang［27］等基于统一相场黏结模型框架，建立了复合型断裂的相场模型，分析了岩石类材料的断裂行为。

本文介绍了相场法模拟准脆性断裂问题的基本原理，详细描述了幂函数能量退化法则与Borden 能量退化法则，提出了一种模拟准脆性断裂问题的相场法耦合能量退化法则，编制了二维相场法计算程序，研究了耦合能量退化法则对单边缺口方形板拉伸断裂行为的影响，进而模拟了直裂缝半圆弯曲岩石试样断裂行为。

1 相场法简介

1.1 相场法基本理论

对于准脆性材料断裂的模拟而言，相场模型可以分为AT2、AT1和内聚裂缝模型等［1］，本文主要介绍AT2相场模型，其总势能可以采用式（1）表示。

式中：ψε(ε，φ) 为Helmholtz 自由能，ψε(ε，φ)=g(c)ψε(ε)；ψε(ε)为弹性应变能密度；ε为弹性应变张量；φ为相场变量，表征了尖锐裂缝弥散为有限裂缝带的拓扑形式，其取值范围为［0，1］，φ=1 表示裂缝，φ=0 表示材料完好；g(c)为能量退化法则函数，c=1-φ；Gc为临界能量释放率；l0裂缝尺度或者特征长度；Γ代表了非连续缺陷；Ω为区域外边界；u为位移场；b和t分别为体力与面力，如图1所示。

图1 裂缝几何正则化Fig.1 Schematic graph of fracture geometric regularization

总势能对位移场和相场变分，同时引入自由能历史最大状态变量H=max[ψε(ε)]作为裂缝扩展的不可逆条件，可得相场模型的控制方程。

式中：g＇(c)=∂g(c)/∂φ；σ为Cauchy 应力张量。控制方程的边界条件可按式（3）进行表示。

1.2 能量退化法则

在相场法基本理论中，能量退化法则耦合了相场与位移场，对材料断裂行为的模拟至关重要。能量退化法则g(c)是随着自变量c单调递增的函数，且满足以下条件［1，4］。

现有的AT2 相场模型中，能量退化法则多为幂函数能量退化准则以及Borden 能量退化法则。幂函数能量退化法则表达式如式（6）所示。

式中：m为幂函数能量退化法则的参数，m≥0。幂函数能量退化法则函数曲线均在g(c)=c2的下侧，且为下凹曲线。Borden能量退化法则可按式（7）进行表述。

式中：s为能量退化法则的参数，0＜s≤3。Borden 能量退化法则函数曲线为介于与g(c)=3c2-2c3（不包括）与g(c)=c3之间的曲线。

根据式（6）和式（7），s=2 时的Borden 能量退化法则与m=0 时的幂函数能量退化法则相同，均为g(c)=c2。参数s=2 的函数曲线为Borden 能量退化法则函数曲线上下限之内的一条曲线，这表明Borden 能量退化法则函数曲线包括了幂函数能量退化曲线不能涉及的一块区域，该区域为0＜s＜2的Borden 能量退化函数曲线。同时，s=3 时的Borden 能量退化法则与m=1 时的幂函数能量退化法则相同，均为g(c)=c3。参数s=3 的曲线为Borden 能量退化法则函数曲线的下限，而m=1 的曲线为幂函数能量退化法则函数曲线上限之下的一条曲线，这表明幂函数能量退化法则函数曲线包括了Borden 能量退化法则函数曲线不能涉及的一块区域，该区域为m＞1 的Borden 能量退化函数曲线。虽然幂函数与Borden 能量退化法则函数覆盖域可以相互弥补，但是能量退化曲线的线型并不能较好融合，这导致材料的损伤演化路径较为单一，无法全面反映材料的力学性质。

2 耦合能量退化法则

为更为全面地模拟材料的断裂行为，通过权重参数对幂函数与Borden 能量退化法则进行耦合，形成的耦合能量退化法则，如式（8）所示。

式中：w为权重参数，0 ≤w≤1。可以证明式（8）在0＜s≤3、m≥0 的情况下，满足能量退化法则的式（4）和式（5）的要求。当w=0时，式（8）为幂函数能量退化法则；当w=1时，式（8）为Borden 能量退化法则。根据式（8），结合式（4）和式（5），可知参数s和m的取值范围并不限于0＜s≤3 和m≥0，为分析方便，本文仅研究0＜s≤3和m≥0下的耦合能量退化法则。

图2显示了耦合能量退化法则参数s=0.1、m=4 时，不同参数w下，耦合能量退化法则的函数曲线。由图2可知，耦合能量退化法则函数曲线涵盖了幂函数与Borden 能量退化法则的所有区域，且函数曲线线型更加丰富，可用于较为全面地反映不同材料的能量退化行为。

图2 耦合能量退化法则函数曲线（s=0.1、m=4）Fig.2 coupled degradation function curves（s=0.1、m=4）

3 数值实现

3.1 有限元离散

位移场和相场问题的残差和可由式（9）和式（10）表示。

式中：g(c，k)=(1-k)g(c) +k；k仅为了数值计算需要，0＜k≪1；Nui和Nφi为形函数；i表示每个单元的第i个节点。采用Newton-Raphson迭代求解，刚度矩阵Kuu和Kφφ可按式（11）和式（12）计算。

根据式（9）～（12）可对位移场与相场进行交错耦合求解。

3.2 程序验证

本文编制了二维相场法计算程序，为了便于与其他学者的研究进行比较，耦合能量法则参数取w=1 和s=2 或w=0 和m=0，此时耦合能量退化法则退化为二次幂函数能量退化法则。

Hesch［28］、Liu［29］和Zhou［30］对单边缺口正方形板进行了拉伸断裂模拟，如图3所示，模型底部固定约束，两侧水平约束，在顶部施加拉伸位移。计算模型参数：弹性模量E为210 GPa、泊松比v为0.3、临界能量释放率Gc为27 kJ/m2、特征长度l0为0.015 mm、k为1×10-7。单元特征长度为0.005 mm。分两步进行加载计算，第一步加载步长为5.0×10-6mm，加载至0.005 mm，第二步加载步长为1.0×10-6mm，直至方形板完全发生断裂。

图3 单边缺口方形板拉伸断裂几何模型与边界条件Fig.3 Geometry and boundary conditions for tension fracture

图4（a）与4（b）显示了裂缝的扩展过程，裂缝沿着与加载方向垂直的方向进行扩展，扩展路径与Liu 等人的模拟较为接近。图4（c）示显示了荷载-位移曲线，就荷载峰值而言，本文结果与Hesch 的结果较为接近，在达到荷载峰值前，本文得到的荷载-位移曲线与Hesch、Liu 和Zhou 的研究基本一致，荷载超过峰值后，本文与Hesch、Liu 和Zhou 所得到的荷载-位移曲线线型基本相同。因此，本文编制的程序对裂缝扩展的模拟以及荷载-位移的捕捉是合理且准确的，可用于分析材料的静态断裂问题。

图4 裂缝扩展路径及荷载-位移曲线Fig.4 Fracture propagation and load-displacement curves of tensile fracture

4 耦合能量退化法则对单边缺口方形板拉伸断裂行为的影响

本节以3.2 节中单边缺口方形板拉伸断裂为例，研究了耦合能量退化法则参数对荷载峰值和荷载-位移曲线的影响。

图5（a）显示了当s=0.1与m=4时，不同w值对单边缺口方形板拉伸断裂行为的荷载-位移曲线的影响，其中w=0与w=1.0分别对应了m=4 的幂函数能量退化法则及s=0.1 的Borden 能量退化法则。可以看出，随着参数w的增大，荷载峰值及其对应的位移值都逐渐增加，且w=1.0 时的荷载峰值是w=0 时荷载峰值的2.45倍，这表明参数w对荷载峰值影响较大。同时，在荷载峰值前，荷载-位移曲线的非线性程度随着参数w的增加而减小，当w=1.0 时，荷载-位移曲线表现为线性，这说明了随着参数w的增加，材料的断裂行为由准脆性断裂向脆性断裂转变，即参数w影响了材料的断裂进程。

图5（b）显示了当w=0.5 与m=4 时，不同s值下的单边缺口方形板拉伸断裂行为的荷载-位移曲线，随着参数s的增加，荷载峰值变小。在w=0.5与m=4情况下，位移-荷载曲线均表现出非线性，随着参数s的增加，荷载峰值前后的荷载-位移曲线的非线性程度增大。这说明参数s影响了材料的断裂进程，可以反映材料的准脆性断裂行为。

图5（c）显示了当w=0.5 与s=0.1 时，不同m值下的单边缺口方形板拉伸断裂行为的荷载-位移曲线。随着参数m的增加，荷载峰值有所降低，但相对图5（a）与图5（b）而言，荷载峰值下降幅度较小。同时，随着参数m的增加，峰值荷载前的荷载-位移曲线非线性程度稍有增加，但是相对于参数w和s对荷载-位移曲线非线性的影响程度而言，参数m对荷载-位移曲线非线性行为影响较小。这表明参数w和s对单边缺口方形板拉伸断裂行为的影响大于参数m对单边缺口方形板拉伸断裂行为的影响。

图5 单边缺口方形板拉伸断裂的荷载-位移曲线Fig.5 Load-displacement curves of tensile fracture

综上可知，耦合能量退化法则可以反映受拉作用下单边缺口方形板的准脆性断裂以及脆性断裂行为，尤其是参数w和s对单边缺口方形板拉伸断裂的进程影响较大，它们对荷载峰值影响也较大。

5 直裂缝半圆弯曲岩石试样断裂模拟

直裂缝半圆弯曲（NSCB）试样为国际岩石力学和岩石工程学会（ISRM）建议的测量岩石断裂性能的方法之一［31］，已广泛应用于岩石、黏性土、混凝土及砂浆材料的断裂性能的测试。本节模拟了NSCB试样的I型断裂试验，主要分析耦合能量退化法则各参数对裂缝扩展路径及荷载-位移曲线的影响。由于NSCB 试验，需要对试样施加压力，则应对应变能进行分解，本文采用Miehe［32］等人提出的方法对对应变能进行分解。

直裂缝半圆弯曲试样几何模型（见图6）与计算参数参考Zhou 等人［33］对岩石材料的研究，选取弹性模量E为92 GPa、泊松比v为0.21、临界能量释放率Gc为9.6 J/m2、材料特征长度l0为0.40 mm、计算参数k为1.0×10-9。模型单元的特征长度为0.083 mm。

图6 直裂缝半圆弯曲试样I型断裂几何模型与边界条件Fig.6 Geometry and boundary conditions for NSCB

图7显示了不同耦合能量退化法则参数下的NSCB 试样I型裂缝扩展的荷载-位移曲线。当s=0.1 与m=4 时，随着参数w的增加，荷载峰值及其对应的位移均会增大，参数w=1.0时荷载峰值是参数w=0 时荷载峰值的2.43 倍。同时，当参数w的增量相同时，随着参数w的增加，荷载峰值增量先减小后增加。此外，随着参数w的增加，荷载峰值前的荷载-位移曲线由非线性特征变为线性，这表明随着参数w的增加，初始裂缝的开裂扩展会由准脆性变为脆性，即耦合能量退化法则中参数w可以反映NSCB试样的准脆性断裂及脆性断裂。

图7 NSCB试样I型断裂的荷载-位移曲线Fig.7 Load-displacement curves of NSCB

当w=0.5 与m=4 时，随着参数s的增加，NSCB 试样I 型裂缝扩展的荷载峰值减小，参数s=0.1 时荷载峰值时s=3.0 时荷载峰值的1.38 倍。同时，当参数s=3.0 时，荷载峰值前的荷载-位移曲线表现为非线性，即NSCB 试样发生准脆性断裂。当参数s=0.1 时，荷载峰值前的荷载-位移曲线表现为线性，即NSCB 试样发生脆性断裂。随着参数s的增加，荷载峰值前的荷载-位移曲线由线性逐渐转变为非线性，这表明参数s影响了NSCB试样的I型断裂进程。

当w=0.5 与s=0.1 时，随着参数m的增加，NSCB 试样I 型裂缝扩展的荷载峰值与荷载-位移曲线几乎不变，这表明参数m对NSCB试样I型裂缝扩展的影响较小。

综上可知，耦合能量退化法则可以反映NSCB试样I型裂缝扩展的准脆性行为和脆性行为，通过改变耦合能量退化法则的参数可以控制NSCB 试样I型裂缝的扩展进程，尤其是参数w和参数s对荷载峰值与荷载-位移曲线影响都较大。

图8、9分别显示了当s=0.1、m=4、w=0.2和s=0.1、m=4、w=0.8时NSCB 试样I 型裂缝扩展在荷载峰值与荷载降为0 的相场变量，其可以反映裂缝扩展情况。在荷载峰值时，两种参数下NSCB 试样的相场区域较小，s=0.1、m=4、w=0.2 时NSCB 试样的相场区域小于s=0.1、m=4、w=0.2时NSCB 试样的相场区域，且s=0.1、m=4、w=0.2 时NSCB 试样的相场变量最大值0.79 大于s=0.1、m=4、w=0.2 时NSCB 试样的相场变量最大值0.59。在荷载降为0时，两种参数下NSCB试样的相场变量最大值达到1.0，裂缝沿着初始裂缝直线扩展，且s=0.1、m=4、w=0.2时NSCB 试样的相场区域宽度大于s=0.1、m=4、w=0.8 时NSCB 试样的相场区域宽度，而两种参数下NSCB 试样的相场变量最大值区域基本一致。

图8 s=0.1、m=4、w=0.2时的NSCB试样I型裂缝扩展Fig.8 Fracture propagation of NSCB（s=0.1、m=4、w=0.2）

图9 s=0.1、m=4、w=0.8时的NSCB试样I型裂缝扩展Fig.9 Fracture propagation of NSCB（s=0.1、m=4、w=0.8）

耦合能量退化法则各参数较难通过试验进行直接测量，可根据试验数据，通过反分析得到。通过上述分析可知，参数w和s对NSCB试样I型裂缝扩展影响较大，参数反演过程中，可着重考虑这个两个参数。同时，为对NSCB试样I型裂缝扩展进行合理地模拟，精确反映荷载-位移曲线的型式，耦合能量退化法则中引入参数w和s的十分必要。