杨云浩，陈鸿杰，王 伟

（1.中国水电顾问集团成都勘测设计研究院，成都 610072；2.河海大学a.岩土力学与堤坝工程教育部重点实验室；b.岩土工程研究所，南京 210098）

弹塑性各向异性损伤模型的FLAC3D开发与数值验证

杨云浩1，陈鸿杰2a，2b，王 伟2a，2b

（1.中国水电顾问集团成都勘测设计研究院，成都 610072；2.河海大学a.岩土力学与堤坝工程教育部重点实验室；b.岩土工程研究所，南京 210098）

探讨适用于硬岩的各向异性损伤模型在FLAC3D中的实现方法。首先，基于不可逆热力学框架下的弹塑性损伤本构基本理论，结合弹塑性力学塑性势与非关联塑性流动法则，导出了弹塑性各向异性损伤本构模型的增量形式及损伤演化速率的计算公式，在此基础上，采用复合M C准则并考虑黏聚力随损伤的逐步弱化，详细给出了该损伤模型在FLAC3D中的实现流程并探讨了相关细节处理方法，编写了相应的DLL文件以供调用。进而，通过圆柱试样的压缩数值试验验证其正确性与合理性，结果表明该损伤模型能够较好地描述岩石材料的非线性力学行为；以一个概化的地下洞室分层开挖实例初步检验了模型的应用效果，计算所得洞周损伤较大部位与实际工程中直接或间接观测到的围岩损伤情况较为一致，由此表明以拉应变度量的损伤及其演化方程能够合理地反映硬岩的损伤本质。

不可逆热力学；弹塑性；各向异性损伤；FLAC3D；本构模型

1 研究背景

在岩体地下工程领域，施工期围岩的稳定性长期以来都是采用常规弹塑性本构模型进行分析。对于硬岩地下洞室，开挖卸荷引起的原岩应力重分布导致围岩内微裂纹大量萌生，微裂纹的逐步扩展使得围岩的刚度和强度属性发生劣化［1－2］，从而在宏观的应力应变曲线上表现出非线性特征，而传统的弹塑性模型通常不能反映这些特征，相应的稳定性分析成果的可靠性也就大打折扣。

近年来，不少学者引入损伤的概念来描述岩石受载过程中变形和强度属性的劣化并建立了相应的损伤本构模型。文献［3］研究了弹脆塑性损伤本构模型，但未明确损伤的启动与演化准则。文献［4－6］所提出的弹塑性损伤本构模型将损伤变量取为标量，以应变能或等效塑性应变作为损伤演化参量。然而，文献［7－8］指出：张拉性裂纹的出现是低围压条件下硬岩损伤的主要形式，其本质是拉裂纹萌生、扩展直至裂纹相互贯通的结果。因而，以应变能或等效塑性应变作为损伤演化参量比较适合于黏土岩和泥岩等软岩，却不能很好地反映硬岩损伤的本质。文献［1］提出的损伤本构模型将体积应变作为损伤参量，体现了硬岩损伤细观上的裂纹萌生扩展导致宏观上出现体胀这一典型特征，但该模型只能考虑各向同性损伤。文献［9］通过大理岩三轴卸荷试验发现试样主破裂面周边常有较多的近垂直于卸荷方向的微张裂隙，由此表明试样轴向与径向的损伤必然是不同的。因此，对于硬岩，有必要考虑损伤的各向异性。文献［10－11］分别针对花岗岩和混凝土所建立的损伤模型考虑各向异性损伤，将二阶损伤张量的损伤主值表达为拉应变的幂函数，但该模型仅考虑弹性损伤且其损伤演化方程在理论上不够严谨。文献［12］所提出的弹塑性各向异性损伤模型将二阶损伤张量与应变张量的正锥直接关联并假定二者共轴。该模型建立在不可逆热力学基础上，理论依据严谨，所采用的损伤准则及损伤演化方程也能很好地反映硬岩损伤的物理机制，模型参数少且物理意义较明确，有良好的应用前景。

借助大型商业软件提供的二次开发接口，编制可供其调用的用户自定义本构模块，是将最新理论研究成果快速应用于工程实践的有效途径。目前，在硬岩各向异性损伤模型的FLAC3D二次开发方面所开展的工作还较少。鉴于此，本文在文献［12］的理论公式基础上，首先导出了弹塑性各向异性损伤本构模型的增量形式及二阶损伤张量增量的具体表达式，进而，采用复合M C准则并考虑黏聚力随损伤的逐步弱化，详细给出了该本构模型在FLAC3D下的实现流程并探讨了相关细节处理方法。

2 弹塑性各向异性损伤本构模型理论

2.1 基于Helmholtz自由能势的弹塑性各向异性损伤本构模型

考虑损伤的各向异性，Dragon等人采用了如下形式的Helmholtz自由能势［12］，即

式中：λ，μ为拉梅弹性常数；εe为弹性应变张量；D为二阶损伤张量。符号“tr”表示对张量求迹；“：”表示张量的双点乘；等式右边后2项为损伤引起的对自由能的修正项，a1，a2为待定参数，用于描述损伤所造成的弹性属性的劣化，其单位与λ，μ相同。

求ψe（εe，D）关于εe的偏导，可得式（2）所示弹塑性损伤本构关系［12］

式中：C（D）为弹性损伤刚度阵，是一个四阶张量，其各分量如式（4）所示；ε，εp分别为总应变和塑性应变张量；δij为Kronecker记号。

2.2 本构模型的增量形式

为便于数值实现，对式（3）的第一式两边关于虚拟时标求导，得到增量形式本构方程为

通过引入塑性屈服准则函数fp（σ，γp）和塑性势函数g，可将式（5）表达为更为具体的形式，简要推导如下。

首先，由塑性一致性条件：f·p（σ，γp）＝0，有

式中γp为广义塑性剪应变，作为塑性硬（软）化参量，计算公式如下：

式中λp为塑性乘子，由式（9）可得

由式（10）可见，如果未发生损伤的演化（此时该式分子部分的第二项为零），则λp退化为常规弹塑性下的形式。

式中符号“”表示并矢运算。

与常规弹塑性本构关系式对比可见，式（11）等号右边第一部分即为常规的弹塑性应力修正量，后两部分则为损伤相关应力修正量。同样，若不存在损伤演化，该式退变为常规的弹塑性本构方程式。

2.3 损伤准则函数与损伤演化方程

由上述可见，弹塑性各向异性损伤本构模型的关键是损伤的度量，也即损伤准则及损伤演化方程的建立。在热力学框架内，损伤演化法则是由和损伤张量相关的共轭损伤力空间内的损伤耗散势所决定的［13］。文献［13］给出了损伤准则的一般形式，文献［14］进一步给出了一个较为具体的形式，但因过于复杂，数值实现较为困难。文献［12］基于室内试验观测，将损伤与拉应变直接关联，建立了如下损伤准则函数，即

进一步假定损伤演化率张量D·与ε＋共轴，则可得如下损伤演化方程

损伤乘子λd的具体形式由损伤一致性条件求得，即

式（12）和式（13）所示损伤准则与损伤演化方程形式简单、易于计算，而且物理意义也较为明确，在一定程度上反映了硬岩受载过程中损伤发生与发展的内在物理机制。鉴于此，本文采用该损伤准则与损伤演化方程来描述硬岩受载过程中的各向异性损伤。

3 弹塑性各向异性损伤本构的FLAC3D二次开发

3.1 屈服准则的选取

考虑工程上的实用性，采用复合M C准则及相应的非关联塑性流动法则（如式（15）所示），同时不考虑损伤对塑性屈服的影响。

式中：fs，ft分别为剪切、拉伸屈服函数；gs，gt分别为非关联的剪切、拉伸屈服塑性势函数；，ψ，σt分别为内摩擦角、剪胀角和单轴抗拉强度；cini，cres分别为黏聚力的初始值和残余值。需注意的是，为了考虑黏聚力随损伤发展的逐步弱化，此处借鉴了文［1］的做法：将黏聚力c看作是广义塑性剪应变γp的函数，即，随γp的增长，c呈指数函数形式从初值降至残余值，参数b控制c的下降速度。

约定σ1≤σ2≤σ3（即拉为正、压为负），则在（σ1，σ3）平面上，复合M C准则屈服面曲线如图1［15］所示。图中：h＝0为拉、剪屈服区分界线；剪切屈服区定义为fs＜0且h＜0；拉伸屈服区定义为ft＞0且h＞0；弹性区为fs＞0且ft＜0。

图1 复合M C准则屈服面及屈服区定义［15］Fig.1 Composite M ohr Coulomb failure criterion and domains used in the definition of the flow rule［15］

3.2 弹塑性各向异性损伤本构的数值实现流程

采用前述理论公式，借鉴弹塑性模型的增量式迭代计算格式，建立了如图2所示的弹塑性各向异性损伤模型的数值实现流程。

3.3 模型关键变量的计算方法

由图2可见，在弹塑性各向异性损伤本构数值实现时，在每个计算时步都要更新总应变和弹性应变，塑性应变在单元进入塑性屈服之后要更新，而损伤张量在有损伤的进一步演化时需进行更新。借鉴FLAC3D内置应变软化模型中塑性硬化参量的计算方法，对总应变、弹性应变、塑性应变及损伤的计算采用逐时步累加的方法求取。实际上，FLAC3D处理每个单元时是对其内各个子单元逐一处理的，即主计算程序传给本构程序的应变增量是子单元的应变增量，本构程序输出的也是子单元的应力张量。为求得计算时步内各关键变量的增量，本文采取了体积平均的处理方法。

以时步内总应变增量为例，其计算公式为

式中：VZone为单元的体积；VSubZone（k）为单元内第k个子单元的体积；ΔεSijubZone（k）为第k个子单元的应变增量；Overlayes为单元的体积被使用的次数，取值为1或2；TotSubZones为单元内包含的子单元总数。

图2 弹塑性各向异性损伤本构在FLAC3D内的计算循环流程Fig.2 Com putation procedures of elastoplastic anisotropic damage constitutivemodel in FLAC3D

图3 数值试验获得的不同围压下的应力－应变曲线Fig.3 Stress strain curves under different confining pressures obtained from numerical compression tests

图4 不同围压下数值试样内某单元损伤演化过程Fig.4 Damage evolution under different confining pressures obtained from numerical compression tests

表1 弹塑性各向异性损伤本构模型参数取值Table 1 Parameters used in elastoplastic anisotropic damage constitutivemodel

4 弹塑性各向异性损伤本构模型的数值验证

按前述数值实现流程以及相关考虑，编制了弹塑性各向异性损伤本构的动态链接库（DLL）文件以供FLAC3D调用。为检验数值实现的正确性，进行圆柱试样压缩的数值模拟实验。试样上下两端施加方向相反的速度边界条件模拟轴向加载，试样径向施加均匀应力边界条件以模拟围压施加。数值试验采用的损伤相关力学参数来自文献［12］（如表1所示），得到的不同围压下的应力应变曲线如图3所示。

由应力应变曲线可见，在应力达到峰值前有明显的非线性段，这是由于损伤导致材料刚度逐渐弱化的结果，这一点与真实应力应变曲线峰前段特征比较吻合。同时也可看到，由于采用的复合M C屈服准则考虑了黏聚力的劣化，所以模型亦能在一定程度上模拟岩石类材料的峰后应变软化力学行为。

图4为不同围压条件下试样内某单元的损伤主值随加载过程的变化（D2与D1的量值及演化过程基本相同，此处略），可见：在弹性损伤阶段，损伤值增加平缓；塑性屈服后，损伤值呈加速增长趋势。

分析可知，圆柱试样压缩试验条件下，按FLAC3D的约定，应变以拉为正、压为负，则有ε1＝ε2＞0＞ε3，于是有ε＋1＝ε＋2，ε＋3＝0（ε＋的含义见2.3节所述），由于假定损伤张量D与拉应变张量ε＋共轴，因此，理论上应有D1＝D2，D3＝0。对比理论分析与数值模拟的结果，可见二者基本吻合，这与D3不严格为零的原因与数值计算的精度有关。进一步考察损伤的分布（如图5所示），可见：较大的损伤出现在与试样轴向呈一定夹角的“X”形的窄条带上，与真实岩样三轴压缩试验的最终破裂面形状也较为吻合。

综上所述，本文弹塑性各向异性损伤本构模型的数值实现是正确的。

图5 试样内损伤分布云图Fig.5 Contours of damage obtained from numerical triaxial com pression test

5 弹塑性各向异性损伤模型初步应用效果

建立一高跨比为2.5的城门洞型地下洞室分层开挖概化模型（为节省篇幅，模型图略去）。施加侧压系数为3.0的初始地应力，大主应力方向垂直于洞室轴向，最大值为15.0 MPa。洞室分2层开挖，应用弹塑性各向异性损伤模型模拟洞室开挖后的围岩力学响应，围岩力学参数取表1所给值。

分层开挖过程中，洞周围岩损伤演化情况如图6所示（因损伤主值D2较小、D3近似为零，此处略去其图）。本例地应力场侧压系数及洞室高跨比均较大，因而边墙部位拉应变大于压应变，而顶拱部位则刚好相反，故而边墙部位的损伤远大于顶拱部位，且损伤程度随开挖的进行而进一步增大。

图6 分层开挖过程中的围岩损伤演化Fig.6 Evolution of damages of surrounding rock during staged excavation

此外，损伤的发生与进一步演化必将弱化围岩的刚度，从而导致围岩变形较不考虑损伤影响时有所增加。为进行对比，在保持相同的变形与强度参数值的前提下，将损伤的启动门槛值参数调高，使损伤始终不发生，进行不考虑损伤影响的开挖模拟计算。

图7所示为特征点处的水平位移对比，可见：由于拱脚部位为压应力集中区，拉应变量值很小，损伤也小，因此考虑损伤与不考虑损伤情况下的位移基本相当；边墙部位为卸荷松弛区，拉应变量值大，损伤相应也大，故而考虑损伤情况下的水平位移较不考虑损伤的情况有显著增加。

图7 关键点水平位移对比Fig.7 Horizontal disp lacements of four key points under damage and non damage conditions

6 结 语

（1）对于开挖于完整性较好的硬岩岩体内的地下洞室，围岩应力与变形仿真分析结果的可靠性在很大程度上取决于所采用的本构模型。常规的弹塑性模型因其难以很好地描述硬岩受载过程表现出的非线性力学行为，故而有必要引入弹（塑）性损伤本构模型开展仿真分析。通过FLAC3D的二次开发接口编写自定义损伤模型是开展考虑损伤效应的围岩稳定分析的一条捷径。

（2）弹塑性各向异性损伤本构的应力增量计算公式与常规弹塑性的计算式相比，增加了2个损伤修正项；在无损伤演化发生时，则在形式上可退化为常规弹塑性本构模型。

（3）由于考虑了损伤的各向异性，需要在每个时步内计算单元的应变张量主值及其主向以便更新拉应变张量，因而应用该模型进行计算耗时较多，今后还需在模型的执行效率方面进行改进，使其可以有效地应用到大规模问题的仿真分析中。

［1］ SALARIM R，SAEB S，WILLAM K J，etal.A Coupled Elastoplastic Damage Model for Geomaterials［J］.Com puter Methods in Applied Mechanics and Engineering，2004，193（27／29）：2625－2643.

［2］ GOLSHANIA，YOSHIAKIO，MASANOBU O，et al.AMicromechanical Model for Brittle Failure of Rock and Its Relation to Crack Growth Observed in Triaxial Compres sion Tests of Granite［J］.Mechanics of Materials，2006，38（4）：287－303.

［3］ 刘洪永，程远平，赵长春，等.采动煤岩体弹脆塑性损伤本构模型及应用［J］.岩石力学与工程学报，2010，

29（2）：358－365.（LIU Hong yong，CHENG Yuan ping，ZHAO Chang chun，et al.Constitutive Model for Elasto Brittle Plastic Damage of Coal Rock Mass Due to Mining and Its Application［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（2）：358－365.（in Chinese））

［4］ SHAO JF，JIA Y，KONDO D，et al.A Coupled Elasto plastic Damage Model for Semi brittle Materials and Ex tension to Unsaturated Conditions［J］.Mechanics of Ma terials，2006，38（3）：218－232.

［5］ 贾善坡，陈卫忠，于洪丹，等.泥岩隧道施工过程中渗流场与应力场全耦合损伤模型研究［J］.岩土力学，2009，30（1）：19－26.（JIA Shan po，CHENWei zhong，YU Hong dan，et al.Research on Seepage Stress Cou pling Damage Model of Boom Clay during Tunneling［J］.Rock and Soil Mechanics，2009，30（1）：19－26.（in Chinese））

［6］ 贾善坡，陈卫忠，于洪丹，等.泥岩弹塑性损伤本构模型及其参数辨识［J］.岩土力学，2009，30（12）：3607－3614.（JIA Shan po，CHEN Wei zhong，YU Hong dan，et al.Parameter Identification of New Elastoplastic Dam age Constitutive Model for Claystone［J］.Rock and Soil Mechanics，2009，30（12）：3607－3614.（in Chinese））［7］ MITAIM S，DETOURNAY E.Damage around a Cylindri cal Opening in a Brittle Rock Mass［J］.International Journal of Rock Mechanics and Mining Sciences，2004，41（8）：1447－1457.

［8］ CAIM，KAISER P K，TASAKA Y.Generalized Crack Initiation and Crack Damage Stress Thresholds of Brittle Rock Masses near Underground Excavations［J］.Interna tional Journal of Rock Mechanics and Mining Sciences，2004，41（5）：833－847.

［9］ 黄润秋，黄 达.高地应力条件下卸荷速率对锦屏大理岩力学特性影响规律试验研究［J］.岩石力学与工程学报，2010，29（1）：21－33.（HUANG Run qiu，HUANG Da.Experimental Research on Affection Laws of Unloading Rates on Mechanical Properties of Jinping Mar ble under High Geostress［J］.Chinese Journal of Rock Mechanics and Engineering，2010，29（1）：21－33.（in Chinese））

［10］周维垣，剡公瑞，杨若琼.岩体弹脆性损伤本构模型及工程应用［J］.岩土工程学报，1998，20（5）：54－57.

（ZHOU Wei yuan，YAN Gong rui，YANG Ruo qiong.Elasto brittle Damage Model for Rockmass Based on Field Tests in Laxiwa Arch Dam Site［J］.Chinese Journal of Geotechnical Engineering，1998，20（5）：54－57.（in Chinese））

［11］高政国，黄达海，赵国藩.碾压混凝土的正交异性损伤本构模型研究［J］.水利学报，2001，（5）：58－64.（GAO Zheng guo，HUANG Da hai，ZHAO Guo fan.An Orthotropic Damage Constitutive Model for RCC［J］.Journal of Hydraulic Engineering，2001，（5）：58－64.（in Chinese））

［12］CHIARELLIA S，SHAO J F，HOTEIT N.Modeling of Elastoplastic Damage Behavior of a Claystone［J］.Inter national Journal of Plasticity，2003，19（1）：23－45.

［13］JU JW.On Energy based Coupled Elastoplastic Damage Theories：Constitutive Modeling and Computational As pects［J］.International Journal of Solids and Structures，1989，25（7）：803－833.

［14］MURAKAMI S，KAMIYA K.Constitutive and Damage Evolution Equations of Elastic brittle Materials Based on Irreversible Thermodynamics［J］.International Journal of Mechanical Sciences，1997，39（4）：473－486.

［15］Itasca Consulting Group Inc.FLAC3DUser’s Manual：Theory and Background in FLAC3D［M］.Minneapolis：Itasca Consulting Group，2005.

（编辑：姜小兰）

Implementation and Numerical Verification of Elastoplastic Anisotropic Damage Constitutive M odel in FLAC3D

YANG Yun hao1，CHEN Hong jie2，3，WANGWei2，3
（1.Hydro China Chengdu Engineering Corporation，Chengdu 610072，China；2.Key Laboratory of Ministry of Education for Geomechanics and Embankment Engineering，Hobai University，Nanjing 210098，China；3.Research Institute of Geotechnical Engineering，Hohai University，Nanjing 210098，China）

The implementation of an anisotropic damage constitutivemodel in FLAC3Dfor hard rock is presented in this paper.An incremental form of elastoplastic anisotropic damage constitutive equation is first deduced using the plastic mechanics theory and the elastoplastic damage constitutive theory within the framework of thermodynamics of irreversible process.Meanwhile the concrete form of evolution rate of damage tensor is determined by the consis tence condition of damage criterion.On the basis of these formulas，the detailed flowchart of numerical calculation by using the elastoplastic anisotropic damagemodel and M C yield criterion is established and some details about

the numerical calculation are also discussed.A DLL（Dynamic Link Library）file corresponding to the elastoplastic anisotropic damagemodel iswritten to be called by the FLAC3Dsoftware.The uniaxial and triaxial compression nu merical tests of cylinder sample are then conducted to test the validity of this user defined model.The results show that thismodel is capable of describing the nonlinear behavior of hard rock.The validity of elastoplastic anisotropic damagemodel in practice is further tested by applying it to a numerical simulation of staged excavation of a fictitious underground opening.Result show that the highly damaged zone around the opening determined by the simulation is in agreementwith the actual observation in project area.The result also show that using the tensile strain tensor to measure damage and establish damage criterion grasps the nature of hard rock damage.

thermodynamics of irreversible process；elastoplastic；anisotropic damage；FLAC3D；constitutivemodel

TU452

A

1001－5485（2013）12－0048－06

10.3969／j.issn.1001－5485.2013.12.009

2012－12－26；

2013－01－24

杨云浩（1977－），男，山西黎城人，工程师，博士，从事岩石力学与地下工程方面的研究，（电话）15928690461（电子信箱）haoyun yang2001＠163.com。