3.3　坝基稳定分析的弹塑性损伤演化本构模型

大多数材料与结构，在宏观裂纹出现之前，已经产生了微观裂纹和微观空洞，人们将材料与结构中的这些微观缺陷的出现和扩展称之为损伤。连续损伤力学或者简称损伤力学就是为了研究这一现象而形成的一个力学新分支（谢和平，2004）。

损伤理论是对经典力学的渗透，在一定程度上弥补了其微观和细观研究的不足。这一学科的研究和发展正方兴未艾，国内外也有许多相关的教材、专著和文献出版。但是到目前为止，如何定义损伤变量、解释损伤变量和进行损伤变量的测定一直是损伤力学学科中最有争议的问题之一。

损伤力学中由于引入了损伤变量，损伤变量自身的演变规律同时又对材料的力学行为产生影响，因此损伤力学的核心问题主要有三个：

（1）定义适当的损伤变量来表述微观空隙的宏观力学效果；描述材料中损伤状态的场变量称之为损伤变量，属于本构理论中的内部状态变量（张全胜，2003）。不同的损伤过程，可以选取不同的损伤变量，即使同一损伤过程，也可以选取不同的损伤变量。

（2）建立描述损伤变量演变规律的方程，即建立损伤演化方程；材料内部的损伤是随着外界因素的变化而变化的，如荷载、温度或腐蚀作用。损伤演化方程是为了描述损伤发展而建立的变化方程。

（3）建立描述有损伤的材料力学行为的物性方程，在这个方程里包含有损伤变量。与普通的物性方程不同在于，损伤物性方程包含了反映材料缺陷和裂纹变化的损伤变量，表现为本构方程中弹塑性模量的减小，有效应力或有效应变的增大。

3.3.1　弹塑性损伤演化的一般格式

经典的弹性损伤理论是基于应变等效原理的，其应力—应变关系为：

式中　σij（i，j＝1，2，3）——Cauchy应力分量；

εkl（k，l＝1，2，3）——应变分量。

式中　uk——相对于直角坐标x＝［x1，x2，x3］T的位移分量；

Cijkl——弹性常数。

式中　δij——Kronecher符号；

λ、μ——Lame常数。

将式（3.3-2）和式（3.3-3）代入平衡方程可得：

利用弹性张量的对称性，可以得到对称的二阶偏微分方程：

由于岩土体介质具有显著的非线性特征，J.C.Simo等假设材料损伤直接相关于总应变历史，并基于应变等效性假设和有效应力的概念，提出了应变基各向同性的弹塑性损伤模型。

在均热和小变形条件下，自由能密度函数为：

式中　ψ——材料的各向同性的连续性；

σr——塑性松弛应力张量；

ε——总应变张量；

p——表征一组塑性内变量的矢量；

Φp——与p和σr相关的塑性势函数；

Φ0——无损弹性应变能密度函数。

式（3.3-7）对时间求导数：

结合Clausius-Duhem不等式，并考虑在均热时，gi＝0，则有；

对于弹塑性损伤材料的弹性阶段，增量型弹性损伤本构关系：

其中，弹性损伤切线模量：

式中　Ωd——损伤耗散势；

γ——等效应变，；

R——损伤扩展力，。

本构关系中塑性响应的表征，假设在有效应力空间进行，设屈服函数为：

依据应力应变的正交关系：

式中　——初始弹性应力。

因此弹塑性损伤域用f≤0表示。增量型弹塑性损伤的应力应变关系为：

其中，弹塑性损伤切线模量

3.3.2　损伤场的计算策略

结构损伤分析的过程首先是选择合适的损伤变量，确定其演化方程，然后与结构的平衡、几何、物理条件一起，构成岩体结构损伤的定解问题或者变分问题，用有限元离散结构，求解结构的应力、应变场的损伤场。最后再根据相应的损伤临界条件，判断岩土体的损伤程度和安全程度。

大多数情况下，损伤是和岩土体的状态相关的，即损伤是耦合的。与无损的岩石力学理论相比，由于损伤变量的引入，含损伤结构定解问题的方程数目需要增加，除了本构关系、平衡方程、几何关系、初始条件和边界条件以外，还需引入损伤的演化方程，因而使得损伤岩体的有限元分析变得更加复杂困难。于是，在可能的情况下，往往通过忽略损伤场与应力应变场的相互作用，来减轻损伤耦合计算带来的困难，通常的求解方法如下。

1.全解耦方法

在某些条件下，可以认为损伤对结构中的应力、应变场的影响很小，可以忽略。基于此，可以首先不考虑损伤，利用无损材料的本构关系、平衡方程，通过普通的有限元计算获得无损的应力应变场，然后代入损伤演化方程，得到对应的损伤场随时间或载荷的变化历史，进而可以进行相关的工程稳定性分析。

全解耦方法相对比较简单，无需重新编制相应的有限元程序，由于损伤而增加的工作量很小。同时，由于应力分析过程中没有考虑损伤耦合的影响，使得全解耦方法的结果往往偏于保守，在损伤影响较大时，计算结果误差很大。

2.全耦合方法

事实上，当结构中的损伤积累到一定程度，将导致结构弹性模量等材料参数和力学性能的变化，造成应力和应变的重新分布。全耦合方法即是在应力应变场的计算中计入损伤的影响。采用含损伤的本构关系，补充损伤演化方程，直接得到所有场变量的分布。用全耦合的方法进行稳定性分析是严格和准确的，但相应的工作量也大幅度的增加。而且，当前的商业有限元程序一般不包括全耦合损伤本构关系，因此，需要发展相应的计算有限元软件或进行二次开发。到目前为止，仅极少数很简单的问题得到了全耦合的解析解。

3.半解耦方法

半解耦方法是介于全解耦方法和全耦合方法之间的一种方法。通常的一种做法是在本构关系中引入损伤，而在平衡方程中不考虑损伤的影响。此方法的工作量比全耦合方法小，而解的精度比全解耦方法高。

3.3.3　损伤场的半解耦计算

弹塑性耦合损伤场的数值计算，主要困难在于迭代中刚度退化引起的非线性、不可逆变形引起的非线性，引入损伤造成的刚度矩阵的非对称性，而弹塑性-损伤耦合的本构模型，使得迭代计算的非线性程度更加严重。

考虑损伤耦合的情况下，应变增量表示为：

则弹塑性损伤耦合本构关系可表示为：

耦合损伤问题迭代计算的基本过程与普通的弹塑性迭代格式是基本一致的：

（1）通过弹性计算，确定初始应力应变σ0，ε0，达到初始平衡状态；

（2）施加足够小的增量步，通过屈服函数，判断单元应力状态是否进入塑性阶段，并得到；

（3）通过损伤势函数Ωd，确定，以及相应的。

为了减轻计算中的非线性困难，本书主要介绍采用半解耦方法建立弹塑性损伤本构关系。

在弹塑性本构方程中引入损伤变量，认为材料的弹性变形不产生损伤，进入塑性的同时进入损伤演化过程，选择Drucker-Prager准则作为屈服准则和损伤阈值的判据。损伤演化方程在唐春安（2003）的基础上进行了改进，以进入塑性屈服作为损伤起始点，而不是采用特定应变值作为损伤阈值。随着损伤破坏程度的增加，损伤系数ω也不断加大。对于一维损伤变量，按下式确定：

式中　λ——单元的残余强度系数；

εco——单元的最大压缩主应力达到其单轴抗压强度时对应的最大压缩主应变；

εto——极限拉伸应变；

f——屈服函数。

按照Mazars把一维损伤本构推广到三维的方法，可以把该本构关系推广到三维应力状态，用等效应变分别代替两式中的ε

式中，ε1，ε2，ε3分别为主应变；＜＞是一个函数，定义为：

假定单元初始状态是弹性的，随着单元应力的增加，当应力应变状态满足Drucker-Prager屈服准则时，材料单元进入塑性，同时进行应力状态判断：假如应力状态位于应力空间的拉伸区域时，单元按照拉伸损伤演化方程进行损伤演化，进入刚度退化处理；假如应力状态位于压剪区域时，单元按照剪切损伤演化方程进入损伤，同样进行刚度退化处理。也就是说Drucker-Prager屈服准则同时作为损伤阈值判据，根据塑性屈服点是拉断屈服形式还是压剪屈服形式来选择不同的损伤演化方程。

在半解耦损伤有限元计算中，其基本迭代格式与弹塑性有限元计算相同。在迭代过程中，损伤只是降低了材料的刚度，没有产生新的非线性。单元进入损伤后，进行相应的刚度退化处理，将单元的刚度矩阵修正为单元的损伤刚度矩阵［HD］再进入有限元的迭代计算。

弹塑性损伤有限元迭代的基本格式如下：

每迭代一次后，根据单元的应变状态按式判断单元是否进入损伤破坏，对于进入损伤破坏的单元应按式（3.3-21）和式（3.3-22）计算单元的损伤变量，并据此修正损伤应力矩阵和损伤刚度，再进行下一步迭代。