混凝土结构在爆炸、冲击等强动载作用下具有显著的应变率效应^[1-2]。Bischoff等^[3]通过分析大量试验数据指出，虽然受试验方法和设备条件等因素影响，不同学者得出的结果存在一定差异，但学者们普遍认为：混凝土的抗拉强度和抗压强度均随着应变率的升高而增加，相对于抗压强度，抗拉强度增加趋势更明显，应变率敏感性更显著。考虑到爆炸冲击试验耗费高且数据离散性大，近年来，众多学者借助数值模拟方法来研究强动载作用下的混凝土结构性能。材料的动态本构模型是数值模拟的基础，ABAQUS程序提供的混凝土损伤塑性（concrete damage plasticity，CDP）模型可用于单向加载、循环加载及动态加载等情况，具有较高的精度和较好的收敛性^[4]，目前已得到了广泛应用。CDP模型最早由Lubliner等^[5]提出，Lee等^[6]对拉压损伤指标进行了改进，通过非弹性应变与受压损伤因子、开裂应变与受拉损伤因子的关系曲线定义了混凝土材料塑性变形阶段的损伤演化。杨宏明等^[7]建立了二维混凝土凸骨料细观模型，运用CDP模型分析了不同应变率下混凝土试件的裂纹扩展及走向，结果表明，随着应变率的升高，二维混凝土试件的裂纹数量明显增加，裂纹扩展速度加快，走向趋于复杂。周禹辛等^[8]基于ABAQUS程序对不同强度等级混凝土冻融后的动态力学性能进行了模拟，结果表明，相同冻融循环次数下，随着加载速率的降低，混凝土的强度逐渐下降。李敏等^[9]指出现有CDP模型仅能考虑混凝土强度的率相关性，无法模拟不同应变率的拉压损伤程度，利用其对爆炸冲击等强动载作用下的混凝土结构动力响应进行分析时误差较大。阎石等^[10]假定爆炸作用下混凝土应变率为100 s⁻¹，通过动态强度放大系数计算得到混凝土强度并输入到CDP模型中，对混凝土板的破坏模式进行了数值分析。实际上，结构在响应过程中应变率是不断变化的，这种简化计算局限于应变率恒定为100 s⁻¹的理想情况，应用于典型强动载作用下混凝土结构的动力响应分析时误差较大。

由于ABAQUS程序的CDP模型中输入的材料参数仅能考虑混凝土强度的应变率效应，不能包含拉压损伤的率效应，已有的研究成果一般只能采用忽略损伤或仅考虑单一加载条件下的混凝土强化和损伤特性，导致计算误差很大。因此，迫切需要解决CDP模型中损伤因子的应变率相关计算难题。本文中，拟通过定义新的应变率场变量，编制VUSDFLD用户子程序，开发能够考虑损伤因子率相关性的改进的CDP（modified CDP，MCDP）模型，开展混凝土材料拉压损伤与应变率相关研究，并对MCDP模型进行对比验证和讨论。

1.   ABAQUS程序中的CDP模型

ABAQUS程序中的CDP模型基于弹塑性理论框架开发，并引入了损伤力学参数。屈服面函数主要包括硬化变量$\tilde \varepsilon _{\text{c}}^{{\text{pl}}}$和$\tilde \varepsilon _{\text{t}}^{{\text{pl}}}$，这2个变量分别表示压缩和拉伸等效塑性应变，如图1所示。CDP模型采用损伤因子描述混凝土材料由损伤引起的刚度退化和导致的拉压屈服强度改变。

图  1  混凝土塑性损伤模型的单轴应力-应变曲线

Figure  1.  Uniaxial stress-strain curves of CDP model

下载: 全尺寸图片 幻灯片

在单轴受压状态，混凝土材料在达到初始屈服应力${\sigma _{{\text{c0}}}}$后先硬化后软化，直至压碎破坏；在单轴受拉状态，混凝土材料在拉伸屈服后表现为软化，直至开裂。在CDP模型中，用户输入的应变为受压非弹性应变$\tilde \varepsilon _{\text{c}}^{{\text{in}}}$和受拉开裂应变$\tilde \varepsilon _{\text{t}}^{{\text{ck}}}$，ABAQUS程序自动转化为塑性应变^[11]，压缩塑性应变$\tilde \varepsilon _{\text{c}}^{{\text{pl}}}$与受压非弹性应变$\tilde \varepsilon _{\text{c}}^{{\text{in}}}$、拉伸塑性应变$\tilde \varepsilon _{\text{t}}^{{\text{pl}}}$与受拉开裂应变$\tilde \varepsilon _{\text{t}}^{{\text{ck}}}$的关系为：

式中：${\sigma _{\text{c}}}$为压应力，${\varepsilon _{\text{c}}}$为压应变，${E_0}$为初始弹性模量，${d_{\text{c}}}$为压缩损伤因子，$\varepsilon _{{\text{0c}}}^{{\text{el}}}$为弹性压应变，${\sigma _{\text{t}}}$为拉应力，${\varepsilon _{\text{t}}}$为拉应变，${d_{\text{t}}}$为拉伸损伤因子，$\varepsilon _{{\text{0t}}}^{{\text{el}}}$为弹性拉应变。

爆炸荷载作用下混凝土结构的材料应变率可以达到10²～10³ s^–1[3]，混凝土抗拉强度和抗压强度随应变率升高而增加。在CDP模型中，输入混凝土的真实应力和非弹性应变数据时，可以通过勾选增加应变率相关的选项，实现在模型中考虑混凝土强度的应变率效应。而在输入损伤参数的界面中则没有应变率选项，无法实现混凝土材料拉压损伤与应变率相关。实际上，不同应变率下的真实应力-真实应变曲线对应不同的损伤参数，只有在CDP模型的损伤参数中充分考虑损伤的应变率相关特性，才能较准确地模拟混凝土结构在爆炸荷载作用下的损伤破坏行为。

2.   考虑损伤因子率相关性的MCDP模型

2.1   不同应变率下的混凝土应力-应变曲线

CDP模型中通过引入动态增长因子（dynamic increase factors，DIF），即混凝土的动态强度与静态强度的比值，来表征材料强度的率相关性。因此，可以依据欧洲-国际混凝土协会（CEB）规范^[12]中的DIF与应变率的关系模型，计算不同应变率下的混凝土动态强度；并结合GB 50011—2010《建筑抗震设计规范》^[13]附录C中的混凝土单轴应力-应变曲线（见图2），用动态强度值替换静态强度值；利用式(1)和(2)，可分别计算CDP模型中需要输入的混凝土材料在不同恒定应变率下的真实应力-受压非弹性应变和真实应力-受拉开裂应变全曲线。

图  2  混凝土单轴应力-应变曲线

Figure  2.  Uniaxial stress-strain curve of concrete

下载: 全尺寸图片 幻灯片

文献[12]建议的混凝土抗压动态增长因子${\theta _{\text{c}}}$的计算公式为：

式中：${\sigma _{{\text{cd}}}}$为动态抗压强度；${\sigma _{{\text{cs}}}}$为静态抗压强度；${\dot \varepsilon _{{\text{stat}}}}$为准静态应变率，取为3×10⁻⁵ s⁻¹；$\lg {\gamma} = 6.156\alpha - 0.49$，$\alpha = {(5 + 3{\sigma _{{\text{cu}}}}/4)^{ - 1}} = {(5 + 9{\sigma _{{\text{cs}}}}{\text{/}}{\sigma _{{\text{c0}}}})^{ - 1}}$，${\sigma _{{\text{cu}}}}$为立方体混凝土抗压强度（MPa），${\sigma _{{\text{c0}}}}$=10 MPa为参考值。

文献[12]建议的混凝土抗拉动态增长因子${\theta _{\text{t}}}$的计算公式为：

式中：${f_{{\text{td}}}}$为动态抗拉强度；${f_{{\text{ts}}}}$为静态抗拉强度；准静态应变率${\dot \varepsilon _{{\text{stat}}}}$取为3×10⁻⁶ s⁻¹；$\lg \eta = 7.11\delta - 2.33$，$\delta = 1/(10 + 6{f_{{\text{cs}}}}{{{f_{0}^{ - 1}}}})$，${f_{{\text{cs}}}}$为静态抗压强度（MPa），${f_{\text{0}}}$=10 MPa为参考值。

2.2   损伤因子计算方法

损伤因子的计算参见文献[11]，采用基于高斯积分求解的经典损伤理论法，通过能量面积比确定损伤因子。根据Najar的损伤理论^[14]，假设能量耗散过程不可逆，建立混凝土损伤变量的计算公式：

式中：${W_0} = {E_0}{\varepsilon ^2}/2$为无损状态下的应变能；${W_\varepsilon }$为损伤状态下的应变能，在Najar损伤理论中，${W_\varepsilon }$按线性简化求得，取值为$\sigma \varepsilon /2$，如图3中的阴影部分所示。

图  3  Najar线性损伤塑性模型

Figure  3.  Najar’s linear damage plastic model

下载: 全尺寸图片 幻灯片

由图3可知，采用简化公式计算的损伤状态应变能小于实际损伤的应变能。因此，本研究采用高斯积分方法求得${W_\varepsilon } = \int {\sigma {\rm{d}}\varepsilon } = \int {f(\varepsilon ){\rm{d}}\varepsilon }$，代入式(5)得到损伤因子：

式中：$f(\varepsilon )$由不同应变率下的混凝土应力-应变曲线确定。

2.3   损伤因子应变率相关在ABAQUS程序中的实现

在CDP模型中，尚无法考虑混凝土损伤因子的应变率相关性，MCDP模型采用VUSDFLD用户子程序定义了不同应变率下的模型参数。为了在程序中实现通过式(1)和(2)计算出的不同应变率下等效塑性应变对应不同的非弹性应变、开裂应变和损伤参数，需要在CDP模型的基础上创建一个描述等效塑性应变率的场变量。通过与ABAQUS主求解程序对应的Fortran程序接口实现应力、应变及对应的损伤因子随场变量变化。在子程序中更新材料点的场变量，则主求解程序中的模型参数也将不断被更新。ABAQUS中用户子程序开发流程如图4所示，调用VUSDFLD子程序进行有限元分析的步骤如下。

图  4  ABAQUS二次开发流程图

Figure  4.  Flow chart of secondary development of ABAQUS

下载: 全尺寸图片 幻灯片

(1) 在CAE模块中创建部件生成几何模型，并划分网格。

(2) 在输入文件中，编辑关键字*INITIAL CONDITIONS，TYPE=FIELD激活子程序。

(3) 定义场变量个数，利用Fortran语言编写场变量计算程序。创建参数，给部件赋材料属性。

(4) 装配部件，创建分析步，设置输出变量。

(5) 创建约束，施加荷载。

(6) 创建分析作业，输入子程序.for文件路径，提交分析。

(7) 在VIEW模块中进行可视化后处理。

参考ABAQUS使用说明中的VUSDFLD子程序通用格式^[15]，使用Fortran语言编写应变率场变量程序，计算框架如图5所示。通过调用内置程序VGETVRM得到当前增量步开始时材料积分点的等效塑性应变ε_i(t)，同时利用内置程序SetField得到上一增量步开始时材料积分点的等效塑性应变ε_i(t-∆t)，取两者差值的绝对值，得到当前增量步开始时材料积分点的等效塑性应变的增量。当前增量步的应变率为：

图  5  VUSDFLD子程序计算流程图

Figure  5.  Flow chart of the user subroutine VUSDFLD

下载: 全尺寸图片 幻灯片

式中：∆t为时间增量。以增量形式进行编程，在VUSDFLD子程序中将应变率赋值给场变量，通过更新材料积分点场变量，主求解程序调用不同应变率下的非弹性应变、开裂应变及损伤因子数据进行有限元计算。

3.   MCDP模型的讨论

3.1   有限元模型的建立

为了验证和讨论所建立的MCDP模型，在ABAQUS中取一个八节点减缩积分单元C3D8R进行分析，模拟C30等级混凝土在单向受压下的动态性能，设定应变率范围为10⁻⁵～10² s⁻¹。如图6所示，对混凝土施加Y方向位移约束条件，利用ABAQUS/Explicit求解器进行求解。混凝土材料模型采用MCDP模型，在CDP模型基础上增加描述等效塑性应变率的场变量（Field 1），实现应力、应变及对应的损伤因子随场变量变化定义。原CDP模型参数不变，文献[16]给出了其余模型参数推荐值，如表1所示。

图  6  有限元模型

Figure  6.  Finite element model

下载: 全尺寸图片 幻灯片

表  1  MCDP模型参数

Table  1.  Parameters of the MCDP model

膨胀角/ （°）	流动势 偏移量	双轴与单轴抗压 强度之比	不变量 应力比	黏性系数
30	0.1	1.16	0.6667	0.0005

下载: 导出CSV 

| 显示表格

采用混凝土轴心强度标准值确定MCDP模型塑性损伤因子，计算简单且具有较高精度^[17]。因此，选用C30混凝土轴心抗压强度标准值20.1 MPa和轴心抗拉强度标准值2.01 MPa分别作为静态抗压强度和抗拉强度。根据前文计算方法得到不同应变率下混凝土的动态抗压和抗拉强度，如表2所示。输入的混凝土应力-应变曲线及对应的损伤参数如图7～8所示。

表  2  不同应变率下的C30混凝土动态强度

Table  2.  Dynamic strengths of C30 concrete under different strain rates

应变率/s−1	抗压动态增长因子	动态抗压强度/MPa	抗拉动态增长因子	动态抗拉强度/MPa
10−5	1.023	20.56	1.138	2.29
10−4	1.055	21.20	1.175	2.36
10−3	1.169	23.49	1.307	2.63
10−2	1.295	26.02	1.453	2.92
10−1	1.434	28.82	1.616	3.25
1	1.588	31.93	1.796	3.61
10	1.760	35.37	1.997	4.01
100	2.775	55.78	3.154	6.34

下载: 导出CSV 

| 显示表格

图  7  混凝土的动态应力-应变关系曲线

Figure  7.  Concrete’s dynamic stress-strain curves

下载: 全尺寸图片 幻灯片

图  8  不同应变率下的损伤因子

Figure  8.  Damage factors under different strain rates

下载: 全尺寸图片 幻灯片

3.2   MCDP模型的退化验证

受试验装置和试验手段的限制，目前尚没有可靠可用的混凝土恒定高应变率加载的试验数据。因此，为了验证MCDP模型和VUSDFLD子程序的正确性，拟对MCDP模型的应变率进行退化处理，模拟静态加载工况，此时应变率设定为恒定10⁻⁵ s⁻¹，并将计算结果与CDP模型静态加载计算数据进行对比。图9给出了采用MCDP模型和CDP模型计算得到的真实应力-真实应变曲线的对比结果，可以发现，经过退化后MCDP模型计算的静态加载结果与原CDP模型吻合较好。

图  9  两种模型得到的真实应力-真实应变曲线的对比

Figure  9.  Comparison between true stress-true strain curves calculated by two models

下载: 全尺寸图片 幻灯片

3.3   恒定应变率单轴压缩动态加载

在恒定应变率加载条件下，分别选用MCDP模型和CDP模型进行计算，将混凝土单轴压缩动态应力-应变全曲线计算结果进行对比分析。由于CDP模型仅考虑混凝土强度的应变率效应，仅可在用户界面输入不同应变率的真实应力-非弹性应变数据，无法计入应变率相关的拉压损伤。而MCDP模型能有效地模拟损伤的应变率相关行为，可在子程序中输入不同恒定应变率下C30混凝土的真实应力-非弹性应变曲线和拉压损伤因子数据（见图7～8）供主求解程序调用。2种混凝土材料模型的压缩模拟结果的对比如图10所示，可以看出，2种计算模型描述的混凝土抗压强度均随应变率的升高而增加，二者模拟曲线的上升段较吻合。不同之处在于：应变率越高，CDP模型描述的混凝土压缩应力越显著高于MCDP模型，这主要是由于CDP模型计算混凝土动态力学行为时未考虑损伤，产生了误差，且应变率越高，误差越大。而MCDP模型能够充分考虑损伤的率相关性，在中高应变率下，应力幅值明显低于CDP模型。如图10(b)所示，在100 s⁻¹恒定应变率下，MCDP模型计算得到的混凝土真实应力-真实应变曲线与公式计算结果吻合，而CDP模型的误差高达20.6%，下降段曲线MCDP模型的计算结果也显著优于CDP模型的计算结果，更好地体现了材料的损伤软化过程。

图  10  恒定应变率单轴压缩动态加载的应力-应变曲线

Figure  10.  Stress-strain curves under uniaxial compression with constant strain rate

下载: 全尺寸图片 幻灯片

3.4   变应变率单轴压缩动态加载

为验证MCDP模型在变应变率下计算的可靠性，分别构造了4种工况分段进行变应变率加载模拟和对比验证，加载工况如表3所示。

表  3  变应变率加载工况参数

Table  3.  Parameters for variational strain-rate cases

工况	模型	时段1/ms	应变率1/s−1	时段2/ms	应变率2/s−1	时段3/ms	应变率3/s−1
A	MCDP	0～1.5	1	>1.5	100	—	—
	CDP	0～1.5	1	>1.5	100	—	—
B	MCDP	0～0.15	10	>0.15	100	—	—
	CDP	0～0.15	10	>0.15	100	—	—
C	MCDP	0～0.13	10	0.13～0.15	100	>0.15	10
	CDP	0～0.13	10	0.13～0.15	100	>0.15	10
D	MCDP	0～0.005	100	0.005～0.055	10	>0.055	100
	CDP	0～0.005	100	0.005～0.055	10	>0.055	100

下载: 导出CSV 

| 显示表格

工况A和工况B采用不同的应变率，分2段进行加载，混凝土变应变率单轴压缩动态加载的应力-应变曲线如图11所示。由图11可以看出，2种模型的初始段曲线与恒应变率加载对应的曲线基本一致，当应变率突然升高至100 s⁻¹时，计算得到的应力-应变曲线明显上升。其中，工况B中CDP模型计算得到的峰值应力达64 MPa，上升段斜率大于MCDP模型，而采用MCDP模型计算得到的峰值应力为55.27 MPa。CDP模型描述的变应变率情况下的混凝土压缩应力均显著高于MCDP模型，充分体现了结构抗强动载作用计算时，在混凝土本构模型中考虑损伤应变率相关的重要性，验证了MCDP模型对于应变率突然升高情况下计算结果的可靠性。

图  11  工况A和B的变应变率单轴压缩动态加载应力-应变曲线

Figure  11.  Stress-strain curves of case A and case B under uniaxial compression with varying strain rate

下载: 全尺寸图片 幻灯片

工况C和工况D分3段采用不同的应变率进行加载，其应变和应变率时程曲线如图12(a)和(b)所示。对混凝土的单轴动态压缩力学性能的预测结果如图12(c)和(d)所示：MCDP模型描述的变应变率下的混凝土压缩应力均显著优于CDP模型，也验证了MCDP模型对于应变率连续突然升高和降低的情况下计算结果的可靠性。

图  12  工况C和D的预测动态压缩应力-应变曲线

Figure  12.  Predicted dynamic compressive stress-strain curves of case A and case D

下载: 全尺寸图片 幻灯片

3.5   爆炸荷载作用下钢筋混凝土梁的响应计算

以Bin等^[18]试验中的钢筋混凝土（reinforced concrete，RC）梁B2为分析对象，采用MCDP模型模拟爆炸荷载作用下RC梁的动态响应。试验梁有限元模型如图13所示，梁截面尺寸为0.2 m×0.2 m，长2.5 m；混凝土强度等级为C40；纵向受拉钢筋为4B20（HRB335级），配筋率$\rho$=1.90%；箍筋为8@150（HRB235级），配筋率$\rho$=0.35%。参照文献[18]，目标梁上的爆炸荷载可简化为分布三角形荷载，如图14所示，荷载峰值分别为：p₁= 4.64 MPa，p₂= 11.27 MPa，p₃= 28.92 MPa。

图  13  混凝土梁的有限元模型

Figure  13.  Finite element model of concrete beam

下载: 全尺寸图片 幻灯片

图  14  梁爆炸荷载简化模型

Figure  14.  Simplified blast load model of beam

下载: 全尺寸图片 幻灯片

梁跨中位移的计算结果与试验结果的对比如图15所示。采用MCDP模型计算得到的跨中位移曲线与试验结果基本一致，跨中位移最大值与试验值的误差较小；而采用原CDP模型计算得到的跨中位移最大值与试验值吻合较差，误差高达38.0%。这是由于CDP模型仅能考虑混凝土强度的率相关性，不能模拟损伤的率相关效应，输入的材料参数不能包含损伤参数^[9]。

图  15  跨中位移计算结果与试验结果的对比

Figure  15.  Comparison between simulation results and test result of mid-span displacement

下载: 全尺寸图片 幻灯片

4.   结　论

定义了新的应变率场变量，基于VUSDFLD子程序建立了考虑损伤因子率相关性的改进的CDP模型（MCDP模型），实现了损伤因子与应变率相关，给出了恒定应变率和应变率突变情况下的混凝土材料单轴压缩应力-应变曲线计算结果，并与CDP模型进行了对比验证，主要结论如下。

(1) 对MCDP模型的应变率进行退化处理，模拟静态加载工况，计算结果与原CDP模型吻合较好，验证了MCDP模型和VUSDFLD子程序的准确性。

(2) 高应变率下，由于未能考虑损伤参数的应变率相关，CDP模型描述的混凝土应力要显著高于MCDP模型，且应变率越高，误差越大，100 s⁻¹时误差高达20.6%。而MCDP模型同时考虑了混凝土强度和损伤的应变率相关性，误差较小，能够更好地体现材料的损伤软化过程。

(3) 采用MCDP模型计算得到混凝土材料变应变率动态压缩加载过程的应力-应变关系曲线，模拟实际受强动载过程中应变率连续突然升高和降低的情况，计算结果同样显著优于CDP模型，进一步验证了提出的MCDP模型的可靠性。

(4) 采用MCDP模型模拟了爆炸荷载作用下RC梁的动态响应，计算结果与试验结果吻合较好，为更准确地预测爆炸冲击等强动载作用下混凝土结构的响应和破坏提供了可靠的技术途径。