一种新的混凝土各向异性弹塑性损伤本构模型及其数值实施

焦延涛，程立平

(1. 华北水利水电大学水利学院，郑州 450045；2. 平顶山学院化学与环境工程学院，平顶山 467000)

目前，国内外关于混凝土塑性-损伤模型的研究成果较多，但对于混凝土这种准脆性材料，由于模型的复杂性，对其进行非线性有限元计算分析仍然是相当具有挑战性的任务。

对于混凝土类材料，其非线性特征主要有弹塑性、粘性、损伤、断裂等。而损伤对材料的非线性行为起着很大的作用，特别是在受拉情况下，但混凝土在受压情况下，在宏观裂缝形成之前，除了破坏外，还表现出明显的塑性性能。由于损伤和塑性变形均对混凝土的非线性行为均具有贡献，因此，单纯的采用弹性-损伤模型来描述混凝土的非线性行为是不够的[1]。

目前，现有的弹塑性损伤本构模型，大多是以热力学为基础，通过连续损伤力学建立相关的热力学势，按照经典塑性力学的方法，用正交流动法则建立热力学基础的损伤演化方程。文献[2 - 13]就是在热力学框架下通过建立耗散势，从而建立了混凝土材料的弹塑性损伤模型。另外，近年来由吴建营等[14-17]提出并发展起来的统一相场损伤理论，也属于此类模型的范畴。这种方法建立的损伤演化方程符合严密的力学推理，有坚实的理论基础，算法的鲁棒性也较好，而且模型也能较为完整的描述混凝土材料在复杂应力状态下的非线性行为，尽管此类模型的优点很多，但损伤耗散势的建立是一项较为困难的工作，对研究者的力学功底要求较高，而且损伤耗散势的试验验证也有一定的难度。考虑到建立损伤耗散势的难度，一种将基于应变的损伤演化方程和基于有效应力的塑性方程结合起来的弹塑性损伤本构模型受到了一些学者的欢迎，文献[18 - 27]即通过这种方法构建了混凝土类材料的弹塑性损伤本构模型。通过这种方法建立的弹塑性损伤本构模型，不需要通过迭代就可以得到损伤变量的值，这样就简化了非线性计算过程，在数值计算时也能获得稳定的解。但由于考虑各向异性损伤的复杂性，大多作者均只采用了一个各向同性的损伤变量，这并不能反映混凝土材料在复杂应力状态下损伤各向异性的特性，也不能反映混凝土材料拉压损伤不等性的特点。

有鉴于此，本文借鉴前人的研究成果，采用两种不同的损伤演化方程，即由双曲型函数表示的拉伸损伤演化方程和由指数型函数表示的压缩损伤演化方程，并由应变来驱动损伤的演化以及控制损伤的各向异性特性。如文献[19, 28]所述，损伤主要是一种应变控制的现象，因此通过应变来驱动损伤的演化是合理可行的。

由于在有效应力空间中建立的塑性损伤模型，其塑性部分在数值计算上更稳定(见文献[2 - 4])，因此本文在有效应力空间中建立了该模型的塑性屈服准则、不同的运动硬化准则和各向同性流动准则。此外，由于通过使用应变等效假设，假定有效构型中的总应变，包括弹性和塑性应变，与名义构型中的总应变相等，因此有效应力和损伤演化的计算可以通过解耦的算法来实现。即计算过程中，应变的计算无需考虑损伤对其的影响，而有效应力的计算在有效构型中进行，也无需考虑损伤对其的影响。如此以来，有效应力以及应变的计算即可按照经典塑性理论的方法来实现，从而可简化本构方程的推导，为数值实现提供优势。另外，采用了解耦的算法后，模型的损伤部分则可进行独立计算，从而可显式实现。

另外，如文献[29]所述，混凝土材料的应力-应变曲线分为上升段及下降段。在应力-应变曲线的上升段，刚度矩阵为正定的，进行弹塑性计算分析时，采用常规的迭代方法可以获得满意的解答，但在应力-应变曲线的下降段，若仍然采用常规的迭代方法，则计算不能收敛，这是因为此时刚度矩阵不再是正定的，可称之为“负刚度”。针对此问题，国内外学者提出了许多相关的方法，以克服这种迭代不收敛情况。目前，常用的这些方法有逐步搜索法、虚假刚性弹簧法、位移控制法、强制迭代法、硬化刚度法、弧长法等。尽管采用上述几类方法，可克服混凝土应力-应变曲线下降段“负刚度”的影响，从而提升计算效率，但深入研究这几类方法后，可发现这几类方法大多通过控制荷载增量，来实现消除“负刚度”的影响。而通过控制荷载增量的方法，则需重组整体刚度矩阵或者分解刚度矩阵，从而会增加模型编程实现的难度。这对于一些希望利用通用有限元软件ABAQUS 的UMAT 用户子程序，或者ANSYS 软件的USERMAT 用户子程序来二次开发实现混凝土弹塑性损伤模型的一些研究者来说，无疑也增加了其编程难度，因为不管采用UMAT 或USERMAT 用户子程序来定义材料本构，均只需定义出材料的应力-应变关系，即Jacobian矩阵，软件即会自动组合总的刚度矩阵，而无需重组或分解刚度矩阵。显然采用上述几种方法会增加程序二次开发实现的难度，为此本文给出了一种与本文模型相适应的迭代算法，以配合ABAQUS 软件的UMAT 用户子程序使用。

综上所述，本文的目的旨在建立一个相对简单的、能反映混凝土拉压不等性及各向异性的特点，且能便于编程实现的混凝土塑性损伤本构模型。值得一提的是，本文工作只适用于小应变条件下的混凝土损伤-断裂分析。

1 损伤变量定义

1.1 各向同性损伤

1) 各向同性损伤变量的定义为了应用连续介质损伤力学的基本原理，首先要对名义构型和有效构型的概念进行解释。名义构型也可以称为损伤构型，而有效构型是一种虚构的存在，也可称其为未损伤构型[20]。如图1(a)所示，在损伤构型中，所有类型的损伤，如空隙和裂纹等都包含在杆件中；而在有效构型中，所有类型的损伤都从杆件中移除，如图1(b)所示。

图1 单轴拉伸荷载作用下混凝土柱有效构型与名义构型对比Fig. 1 The comparison of the nominal and effective configurations of cylindrical bar under uniaxial tension

根据文献[30]所述，可以得到名义应力和有效应力之间的以下关系：

其中：

式中：A、AD和分别为如图1 所示的杆件截面的整个面积、总损伤面积和有效面积；d为标量损伤量，其值在0～1 之间变化；带上划线的物理量为有效构型中的物理量；不带上划线的物理量为损伤构型中的物理量。

式(1)所示的关系只适用于单轴应力状态，当双轴或三轴应力状态下，采用各向同性损伤模型时，有效应力张量与名义应力张量之间的关系可表示为：

理论上，损伤变量可以由式(2)确定，但由于混凝土材料的不均匀性，通过物理试验的方法测定材料的损伤面积是一项复杂的工作[31]。为此文献[32]提出了应变等效假设，以期通过间接手段确定损伤变量。通过对文献[3 - 4]的分析，利用应变等效假设，可以得到以下关系：

利用广义Hook 定律，考虑塑性变形的应力-应变关系可定义为：

式中：E为四阶损伤弹性刚度张量；为四阶未损伤弹性刚度张量。对于各向同性线性弹性材料，可由以下公式给出：

根据式(5)，各向同性损伤变量的间接形式可表示为：

对于一维情况，上述方程可改写为：

式中，E为材料损伤后的弹性模量，而式(8)也即式(2)的等价形式。

2) 应变张量的分解

由于混凝土在拉伸和压缩状态下表现出损伤不同的特性，为了充分描述由于拉伸、压缩和循环荷载造成的混凝土拉压损伤不同的这一特性，使用谱分解技术将应变张量分解为正负两部分[28]。式(9)中，上标“+”和“-”分别表示拉伸和压缩物理量。将总应变张量进行分解，则可得：

式中， ε+和 ε-分别为总应变张量的拉伸和压缩部分。

对称二阶应变张量的谱分解形式可表示为：

总应变张量的正负部分可以表示为：

式中，H(x) 为阶跃函数(当x＞0 时，H(x)=1；当x＜0 时，H(x)=0)。

将式(11)代入式(12)，可得到如下表达式：

式中，I为四阶单位张量，而P+可以表示为：

由于弹性应变张量的归一化特征向量与总应变张量的归一化特征向量相同，因此，弹性应变张量的正负部分可用下式表示为：

根据式(15)可得，弹性应变张量可以表示为：

而根据式(5)和式(15)则可得，两种构型中的正应力张量可表示为：

负应力张量则可表示为：

根据式(16)～式(18)，可得出以下关系式：

根据文献[3]所述，如果假设式(3)所示的有效应力与名义应力的关系，对拉伸和压缩状态均有效，则有：

从而名义应力张量可以表示为：

式中，d+和d-分别为拉、压各向同性损伤变量。值得注意的是，损伤变量不能简单地分解为d+和d-，即d≠d++d-。

1.2 各向异性损伤

各向同性损伤假设在损伤演化过程中混凝土材料强度和刚度的退化在不同方向上相同，这与实际是不符的。因此，为了准确描述混凝土损伤各向异性的特性，本文在假设损伤主轴与应变主轴重合的基础上，采用了二阶对称张量ω来表示损伤。由于ω是一个二阶对称张量，因此，其也可用谱分解表示为：

式中：vi(i=1,2,3)为损伤张量的归一化特征向量；为第i个主损伤量。

由于假设损伤变量的主轴与应变张量的主轴重合，因此，式(22)中的向量vi与式(10)中的向量ni可视为相同，则式(22)可重写为：

由于Ei是个未知量，本文中将由基于应变的损伤演化方程确定。关于损伤演化方程的详细描述将在第2.2 节中给出。

另外，由于定义了二阶损伤张量，因此需对式(3)进行扩展，即：

如文献[20]所述，由式(25)得到的有效应力张量并不是对称的。而非对称的有效应力张量，不管是在本构方程的推导或是在有限元的编程计算分析过程中均是不便于应用的。为此，许多学者提出了不同的对称化方法。本文采用文献[33]采用的方法，对有效应力张量进行对称化处理。通过使用该方法，式(25)可重写为：

或：

其中：

基于谱分解的概念，二阶损伤张量可以表示为：

结合式(19)和式(32)，总的名义应力张量可表示为：

从式(33)可以看出，损伤影响张量M可以表示为：

且：

2 弹塑性损伤模型

2.1 塑性部分

由于混凝土在拉压荷载作用下呈现出不同的材料性能，本文采用文献[34]提出的能同时考虑拉、压塑性的屈服准则。另外，对于混凝土类材料，当计算其应力的下降段时，会出现“负刚度”，从而影响收敛速度，为此本文将模型的塑性部分建立在有效构型中。文献[2]也表明：将塑性损伤模型的塑性部分建立在有效构型中，可与前文所述的应变等效假定相互配合，从而使模型的推导过程大大简化，而数值实施时，其算法也更高效、更稳定。屈服准则在有效构型中可表示如下：

式中，H+和H-定义为：

式中，符号 〈〉表示Macauley 括号，定义为：

和：

式中，Gp为塑性势函数，其表达式将在下文给出。

结合式(40)和式(47)，式(38)可表示为：

参考文献[3 - 4, 34]，对于混凝土类材料，压缩状态下塑性变形阶段会产生明显的体积改变，而采用非关联的塑性流动准则可很好的重现这一现象，因此本文选用非关联流动法则来确定塑性应变张量，即：

本文中选用Drucker-Prager 屈服函数作为塑性势函数Gp，则：

式中， αp为膨胀系数，对于混凝土类材料，其取值范围一般在0.2～0.3，本文取0.2。对式(53)进行求导，则塑性流动方向 ∂Gp/可表示为：

2.2 损伤部分

由于损伤是一个不可逆过程，本文通过损伤加载函数、加载卸载条件和损伤变量的演化规律来描述模型的损伤部分。

损伤加载函数和加卸载条件可以表示为：

由于损伤和塑性变形都是导致混凝土非线性响应的原因，则在单轴荷载作用下和

i可表示为[35]：

1) 拉伸损伤演化方程

式中，d0为初始损伤。参考文献[37]的取值，本文取其值为0.063。

2) 压缩损伤演化方程

参考文献[3]，混凝土弹性刚度的退化在拉伸和压缩状态下表现出不同的特性(如图2 所示)，为此需定义出不同的压缩损伤演化规律。根据文献[38]所述，对于混凝土类材料，在单轴压缩载荷作用下，其塑性行为比损伤产生的早，即当等效应力达到压缩屈服强度时，材料中开始出现塑性行为，而当等效应力达到峰值应力时，则出现新的损伤(即压缩状态下≠。基于此，单轴压缩载荷下本文采用指数软化规律，其峰值后的应力-应变关系可表示如下：

图2 单轴荷载作用下混凝土力学特性Fig. 2 Concrete behavior under uniaxial

3 数值实施

根据本文引言部分所述，采用解耦算法，模型的塑性部分和损伤部分可以分开单独进行求解。另外，如引言部分所述，对于混凝土类材料，在应力-应变曲线的下降段，由于“负刚度”的存在，计算很难收敛。本文将模型的塑性部分建立在有效构型中，材料屈服后，则根据式(36)的屈服准则，经过迭代计算，屈服点后拉伸有效应力-应变曲线如图2(a)中虚线部分所示，屈服点后压缩有效应力-应变曲线则如图2(b)中虚线部分所示。屈服点前，拉伸/压缩有效应力-应变曲线与名义应力-应变曲线重合。从图2 可以看出，采用将模型的塑性部分建立在有效构型中的方法，有效应力不存在下降段，而通过有效应力计算Jacobian 矩阵，可有效避免“负刚度”的出现，从而大大提高收敛速度。但研究混凝土的损伤-断裂状态，实际上研究的是其损伤构型中的力学性态，需要获得的应力-应变曲线实际上是名义应力-应变曲线，当有效应力计算出以后，再根据损伤量的计算结果，即可根据式(33)计算出名义应力。而由于采用应变等效假设，有效构型中的应变与损伤构型中应变相等，从而名义应力-应变曲线即可获得。名义应力-应变曲线的下降段如图2的实线所示，为避免“负刚度”的影响，本文将名义应力作为状态变量，不参与Jacobian 矩阵的计算。

基于ABAQUS 平台，采用UMAT 子程序对本文模型进行二次开发，具体的模型数值求解原理如下。

3.1 塑性部分

本文采用向后欧拉隐式算法进行塑性数值积分，其包括两个方面：弹性预测和塑性修正。增量步n和n-1之间的未知变量可以更新如下：

在 增 量 步n-1 给 出 一 组 变 量(ε(n-1),,) 和 一个应变增量变量 Δε=Δtε˙，则式(64)为一个关于变量 (ε(n+1),,))的非线性方程组。更新变量来自前一个时间歩结束时的收敛值。非线性方程组的求解通过Newton-Raphson 迭代算法来进行。有效应力及相关变量更新过程如下：

1) 初始化：

2) 弹性预测：

3) 检查屈服及收敛条件：

4) 初始迭代：

5) 牛顿迭代：

① 应力迭代

② 检查收敛条件

如果 |l|＜TORL且＜TOLF，进入第6)步；否则进入③。

③ 更新内变量：

④k=k+1，返回①。

6) 更新最终迭代结果：

7) 结束。

3.2 损伤量及名义应力计算

采用解耦算法，本文中损伤量的计算采用显式算法。损伤量获得后，结合有效应力的计算结果根据式(33)计算名义应力。同时为了避免“负刚度”的影响，将名义应力做为状态变量进行更新，不参与刚度矩阵的迭代计算，具体算法流程如下：

4 模型验证

为了验证本文所提出的塑性损伤模型的适用性和有效性，本节将给出不同加载条件下由本模型计算的结果与混凝土试件试验结果的对比，具体如下所述。

4.1 单轴压缩加载-卸载试验

由于塑性本构方程是在有效构型中定义的，因此本文中塑性材料参数的确定采用文献[4]给出的方法。该方法可概括为：首先，通过连接图3所示的每个卸载点(A点～E点)和重新加载点(A'点～E'点)，可以确定每个循环的损伤弹性模量E。进而可由方程=(/E)σ求得有效应力。有效应力得出后，绘制其与相应应变的关系曲线，然后可根据式(49)和式(50)确定出塑性材料参数。对于损伤材料参数，由于本文损伤演化方程是由应力-应变方程推导而来，因此，其可由式(60)和式(62)与单轴拉、压应力-应变曲线拟合确定。

图3 文献[39]试验所得有效构型及名义构型中应力-应变曲线Fig. 3 Experimental stress-strain curves in the effective and nominal configurations from the experimental results of literature[39]

此外，如文献[37]所述，通常由于混凝土试件中初始微裂隙的存在，因此由单轴试验确定的初始弹性模量E0并不是完全无损的弹性模量，其没有考虑初始损伤d0的影响。因此，根据式(8)，可以将初始完全未损伤的弹性模量定义为=E0/(1-d0)。

使用上述方法，根据文献[39]提供的单轴压缩加载-卸载试验曲线拟合的压缩塑性和损伤材料参数数如表1 所示。数值模型计算结果与试验结果的对比如图4 所示。

表1 根据文献[39]试验结果得到的材料参数Table 1 Material constants identified from the experimental results[39]

图4 本文模型计算的单轴压缩加载-卸载结果与文献[39]试验结果对比图Fig. 4 The model responses in uniaxial loading-unloading compression compared to experimental results presented in literature[39]

如图4 所示，模型计算结果能很好地描述峰后区域，但单轴抗压强度的计算结果有点偏高。其主要原因是在试验过程中等效应力在达到单轴抗压强度前，新的损伤已经开始。因此，物理试验时材料的实际损伤是略大于初始损伤d0的，但模型计算时单轴抗压强度是由公式=(1-d0)来获得的，而的值是根据前述方法由试验曲线拟合得到，其值一定的情况下，(1-d0)的值比实际值大，从而即导致了模型计算结果比试验结果偏大。

4.2 单轴拉伸加载-卸载试验

按照第4.1 节中提出的相同方法，由文献[40]的单轴拉伸加载-卸载试验拟合得到的材料参数如表2 所示。将文献[40]单轴拉伸加载-卸载试验结果与模型计算果进行对比，如图5 所示。从图5可以看出，模型计算结果与试验结果吻合良好。

表2 根据文献[40]试验结果得到的材料参数表Table 2 Material constants identified from the experimental results of literature[40]

图5 本文模型计算的单轴拉伸加载-卸载结果与文献[40]试验结果对比图Fig. 5 The model responses in uniaxial loading-unloading tension compared to experimental results presented in literature[40]

4.3 单轴单调压缩试验

单轴单调压缩加载情况，根据文献[39]提供的试验数据拟合的材料参数如表3 所示。模型计算结果与试验结果的对比如图6 所示。

表3 单轴单调压缩试验所用材料参数表Table 3 Material properties used for the monotonic uniaxial compressive test

图6 本文模型计算的单轴单调压缩加载结果与文献[39]试验结果对比图Fig. 6 The model responses in monotonic uniaxial compression compared to experimental results presented in literature[39]

对于单轴单调加载情况，由于第4.1 节采用的方法已不再适用，因此单轴单调加载情况，采用origin 软件根据最佳近似原则对试验数据进行拟合得到材料参数，下文单轴单调拉伸加载情况，材料参数也根据此方法得到。从图6 可以看出，模型计算结果基本能反应实际情况。

4.4 单轴单调拉伸试验

单轴单调拉伸加载情况，根据文献[41]提供的试验数据拟合得到的材料参数如表4 所示。模型计算的应力-应变曲线与试验的应力-应变曲线对比如图7 所示。从图7 可以看出，模型计算结果基本能反应实际情况。

图7 本文模型计算的单轴单调拉伸加载结果与文献[41]试验结果对比图Fig. 7 The model responses in monotonic uniaxial tension compared to experimental results presented in literature[41]

表4 单轴单调拉伸加载试验所用材料参数表Table 4 Material parameters used for the monotonic uniaxial tensile test

4.5 双轴单调压缩加载试验

将文献[42]的双轴压缩试验结果与本文模型的计算结果进行对比，如图8 所示。材料参数如表5所示。

表5 双轴压缩加载试验所用材料参数表Table 5 Material constants used for the biaxial compressive test

从图8 可以看出，当应力比为σ1/σ2=-1/0时，采用本文模型得到的计算结果与试验结果吻合很好，但当应力比为 σ1/σ2=-1/-1 和σ1/σ2=-1/-0.52时，应力-应变曲线的上升段能与试验曲线基本吻合，而下降段则与试验曲线存在一定的差异。究其原因，可能为试验过程中，当应力比为σ1/σ2=-1/0时，材料只有一向受压，另外两向无约束作用，这样受压时，材料的微缺陷与内部裂纹很容易就会扩展，从而导致有效承载面积的减小，也即损伤量的增大；当保持 σ1与单向受压时相同，逐渐增加 σ2，此时，与单向受压时不同，材料在另一向由于受到约束作用，内部的微缺陷及微裂纹不容易扩展，甚至随着 σ2的增大会有闭合的趋势，这就使得材料的有效承载面积减小很慢甚至会增大，从而使损伤量保持不变或减小。而多轴加载状态下，采用本文模型计算时，虽然考虑了损伤各向异性的特点，但是计算三个主应变方向的损伤量，各向实际上均是按照单轴加载状态进行计算的，并没有考虑三个方向之间相互作用对材料损伤的影响。因此，多轴加载状态下，各轴之间相互作用对损伤的影响有待进行更深入的研究。

图8 本文模型计算的单轴及双轴单调压缩加载结果与文献[42]试验结果对比图Fig. 8 The model responses in monotonic uniaxial and biaxial compressive loading compared to experimental results reported by literature[42]

另外，需要说明的是本文利用谱分解技术将应变张量分解为正负两部分，并将对应的应力张量分解为正负两部分。通过这种处理，不仅可以考虑拉压损伤不同的特点，也可方便地考虑泊松比的影响。如图8(c)所示，如果双轴应力比为σ1/σ2=-1/-0.52 ， ε1和 ε2均为压缩应变，考虑泊松比影响，则 ε3为伸长应变。而伸长应变可能导致拉伸损伤，因此在双轴压缩试验中，考虑拉伸材料常数是必不可少的。

5 双边开口四点弯曲梁加载试验

5.1 本文模型计算结果分析

为验证本文模型的有效性，本节采用本文提出的混凝土弹塑性各向异性损伤本构模型对文献[43]所做的混凝土双边开口四点弯曲梁(DEN)试件的断裂破坏情况进行数值模拟。试件厚37.5 mm，截面为400 mm×150 mm 的长方形，槽的尺寸为25 mm×5 mm。试验采用的加载系统为固定支座加载系统，具体如图9 所示。

图9 文献[43]中双边开口混凝土梁试件采用的固定支座加载系统Fig. 9 The fixed loading supports of the DEN specimen in literature[43]

基于ABAQUS 软件的UMAT 子程序，对本文模型进行二次开发，计算过程中采用8 节点六面体线性非协调模式单元(C3D8I 单元)。为模拟固定支座对混凝土梁试件的线性约束作用，在有限元模型中，对支撑处节点的除竖向平动自由度外的其余自由度均约束。本文所用混凝土双边开口梁试件三维有限元网格如图10 所示。

图10 双边开口混凝土梁试件三维有限元网格Fig. 10 3D mesh of the DEN specimen

由于采用力加载方式在数值模拟过程中容易导致计算不收敛，因此本文计算时，采用位移加载方式在图9 所示的加载点F1R 及F2 处分别施加-1 mm 和1 mm 的强制位移荷载，而在F2R 和F1 两个加载点处分别施加-1/15 mm 和1/15 mm 的强制位移荷载。根据第4.3 节所述材料参数识别方法，得到的材料参数如表6 所示。

表6 双边开口混凝土梁试件断裂数值模拟材料参数Table 6 Material parameters used for the DEN specimen fracture simulation

图11～图13 为试件完全破坏时，采用本文模型得到的双边开口梁试件3 个主应变方向上损伤区分布情况。

图14 为文献[43]双边开口梁试验裂纹最终扩展路径。对比图11～图13 的损伤区分布可知，采用本文弹塑性各向异性损伤本构模型得到的损伤区与图14 中所示的文献[43]试验所得的试件断裂区基本吻合。图11～图13 中3 个主应变方向上损伤值存在大于1 小于0 的情况，分析原因可能为ABAQUS 软件在计算过程中将积分点处损伤值外插至节点处，然后进行加权平均计算，从而造成了损伤值不在0～1 范围内这一现象。

图11 第1 主应变方向对应的损伤区分布Fig. 11 Distribution of damage zone corresponding to the first principal strain direction

图12 第2 主应变方向对应的损伤区分布Fig. 12 Distribution of damage zone corresponding to the second principal strain direction

图13 第3 主应变方向对应的损伤区分布Fig. 13 Distribution of damage zone corresponding to the third principal strain direction

图14 文献[43]双边开口梁试验裂纹扩展路径Fig. 14 Crack modes reported in literature[43]

5.2 ABAQUS 自带塑性损伤模型(CDP 模型)计算结果分析

为进一步验证本文模型的有效性，本节采用ABAQUS 软件自带的CDP 模型对5.1 节所述的双边开口四点弯曲梁进行损伤计算。采用CDP 模型计算时，有限元模型、加载条件等与5.1 节完全相同。采用CDP 模型计算时所用材料参数，根据DEN 试件混凝土的强度等级并参考文献[44]所给出的方法进行计算确定，具体材料参数如表7 和表8 所示。

表7 CDP 模型计算时所用混凝土弹塑性参数Table 7 Elasto-plastic mechanic parameters of concrete used by CDP model

表8 CDP 模型计算时所用混凝土损伤参数Table 8 Damage parameters of concrete used by CDP model

由CDP 模型计算的DEN 试件断裂时的损伤云图，如图15 和图16 所示。

从图15～图16 可以看出，由CDP 模型计算的拉损伤区分布，能基本反映试件中裂纹的扩展趋势，但存在明显的分叉裂纹，与实际情况存在一定出入，而且压损伤区分布占较大一块区域，裂纹路径不明显，与实际的裂纹分布情况也存在较大差异。

图15 由CDP 模型计算的拉损伤云图Fig. 15 Distribution of damage in tension calculated by CDP model

图16 由CDP 模型计算的压损伤云图Fig. 16 Distribution of damage in compression calculated by CDP model

5.3 本文模型与CDP 模型计算结果的对比分析

对比本文模型及CDP 模型计算的裂纹路径，可以看出，本文模型计算的裂纹扩展路径更接近试验结果，而且本文模型可以考虑损伤各向异性的特点，而CDP 模型只能考虑混凝土材料拉压损伤不等的特点。

另外，采用本文提出的材料本构模型及数值计算方法，只需20 min 左右即可完成整个模型的计算；而在有限元模型及加载条件完全相同的情况下，采用ABAQUS 软件自带的CDP 本构模型则需约70 min 才能完成最终的计算。两者相比，本文模型可大大提高计算效率。

经对比分析可知，本文模型能反映混凝土损伤各向异性的特点，损伤分布情况也与试验裂纹路径基本吻合，而且也可大大提高计算效率，节约计算机时，因此模型可用于混凝土类材料的损伤计算分析。

6 结论

将损伤力学与经典塑性理论相结合，提出了一种新的塑性各向异性损伤本构模型来模拟混凝土的非线性行为。根据本文的工作，可以得出以下结论：

1) 本文建立了确定各向异性损伤变量的两种不同的损伤演化方程。该方法不仅考虑了混凝土在拉压荷载作用下的不同损伤机理，而且可考虑混凝土损伤各向异性的特点。

2) 在运用应变等效假设的基础上，通过解耦算法将塑性与损伤分离开来，可大大简化本构模型的推导，也可提高模型的计算效率。

3) 通过与大量试验数据的比较，结果表明：该模型能有效地模拟混凝土在不同加载条件下的非线性特性。

4) 对DEN 试件的数值模拟表明：该各向异性塑性-损伤本构模型与ABAQUS 软件自带的CDP本构模型相比，其不仅能反映混凝土损伤各向异性的特点，而且其计算的损伤路径更符合实际情况；从计算效率来说，其效率相比CDP 模型也更高。