具体实施方式

以下结合说明书附图和具体实施例对本发明做进一步详细说明。

实施例：适用于高压、高应变率下的混凝土损伤本构模型构建方法，详细如下描述。

图1是具体实施例中强度面示意图。图中1是屈服强度面，2是最大强度面，3是残余强度面。当前加载面为分段函数，当λ≤λ_m时，通过插值屈服强度面和最大强度面描述混凝土材料应变硬化现象，当λ＞λ_m时，通过插值最大强度面和残余强度面描述混凝土材料应变软化(损伤)现象。即

其中，Y(σ_ij,η)为当前加载面，σ_ij是应力张量，η是强度面比例因子，λ是等效塑性应变，λ_m是峰值应力时等效塑性应变，J₂是第二应力偏量不变量，DIF为动态增强因子，σ_y是屈服强度面，σ_m是最大强度面，σ_r是残余强度面，r'是当前子午线与压缩子午线比值，即

其中，θ是Lode角，即

其中，ψ是拉压子午线比值，s_ij是偏应力张量。

引入修正系数来修正高围压段屈服强度面形式，屈服强度面由多段函数组成，即

其中，r_p＝P/T，T和f_c分别是混凝土单轴抗拉强度和单轴抗压强度，f_yc是混凝土抗压屈服强度，a_1y、a_2y分别是材料常数，a_1y＝0.9088，a_2y＝0.075/f_c，N为修正系数，N的建议取值0.45～2，当N取1时与实验值吻合较好。

当r_p≤0时，峰值压力可达-T，通过插值三轴拉伸点(-T,0)和在零压力时的0.45倍最大强度面(0,1.35T)；

当时，由于缺乏有效的试验数据，对初始屈服点(f_yc/3,f_yc)和在零压力时的0.45倍最大强度面(0,1.35T)线性插值，其中f_yc＝0.45f_c；

当时，根据三轴压缩试验数据定义，σ_y＝0.45σ_m；

当r_p＞Nf_c时，屈服强度面与压力无关。

最大强度面由多段函数构成，即

其中，a₁和a₂是材料常数，由三轴压缩试验数据确定，a₁＝0.5857，a₂＝0.025/f_c。是拉压子午线比值，即

拉压子午线比值由离散点插值决定。由于最大强度面为多段函数定义，当r_p＞f_c时，可根据无侧限单轴压缩和三轴压缩实验数据定义，当0＜r_p≤f_c/3时，没有有效的实验数据，由于压力低于单轴抗压强度，故采用拉伸子午线除以拉压子午线比值的方法。

对残余强度面进行了修正，弥补了混凝土靶体损伤后材料偏软的缺陷。即

其中，a₃是材料常数，a₃＝0.5。

其中，n、m、k是材料参数， m＝0.29，k＝1.86；当λ≤λ_m时，η从0到1变化，此时由屈服强度面向最大强度面变化，捕捉应变硬化现象；当λ＞λ_m时，η从1到0变化，此时由最大强度面向残余强度面变化，捕捉应变软化现象。

图2是本发明实施例损伤机制示意图，本实施例通过引入拉伸损伤指标、剪切损伤指标和静水压缩损伤指标以及剪切损伤对拉伸损伤影响因子和静水压缩损伤对其他应力状态影响因子，确定等效塑性应变。本发明引入了由混凝土内部孔隙坍塌引起的静水压缩损伤，并引入了与应力状态相关的变量，考虑各损伤之间相互影响关系。即

λ＝1-(1-δ_sλ_t)(1-λ_s)(1-δ_hλ_h) (9)

其中，λ_t、λ_s、λ_h分别是拉伸损伤、剪切损伤和静水压缩损伤指标，δ_s是剪切损伤对拉伸损伤影响因子，δ_h是静水压缩损伤对其他应力状态影响因子。即

拉伸损伤指标、剪切损伤指标和静水压缩损伤指标如下所示，

其中，是等效塑性应变增量，dμ_p是塑性体积应变增量，b₁、b₂、b₃是损伤常数，b₁＝1.6，b₂＝1.0，b₃＝1.5。

通过4个简单工况对本发明中损伤机制做出解释。

(1)当混凝土处在受拉状态，r_p＜0，只有λ_t累积，当处于压缩状态时，λ_s和λ_h同时存在，λ_h是由孔隙坍塌引起的塑性体积应变累积而成，当混凝土处于高围压状态，δ_h＝0，λ_h对外没有影响；当混凝土处于静水压力状态，无偏量产生，λ_s＝0，孔隙坍塌后仅存在体积硬化，与现有实验结果一致。

(2)当混凝土应力状态由三轴压缩向非三轴压缩变化时，δ_h＝1，λ_h对λ有贡献，随着三向围压增大，λ_h越大，与现有实验结果一致。

(3)当混凝土从受拉状态变为受压状态，假设之前没有剪切损伤和静水压缩损伤，λ_s＝0，δ_h＝0，所以，λ＝1-(1-δ_sλ_t)，拉伸时，r_p＜0，λ_s＝1，λ＝λ_t，当受压时，r_p≥0，δ_s＝0，λ＝λ_s，显然，拉伸损伤对后续压缩行为没有影响。

(4)当混凝土从受压状态变为受拉状态，假定之前没有拉伸损伤累积，λ_t＝0，所以只有剪切损伤和静水压缩损伤存在，λ＝1-(1-λ_s)(1-λ_h)，因此，前期的剪切损伤和静水压缩损伤对后续的拉伸性能有很大影响。当应力状态由静水压缩状态变为拉伸时，λ＝λ_h反应了孔隙坍塌对混凝土抗拉强度的折减。

为了便于结果输出和后处理，统一了输出的各损伤分量，即

其中，D_i是输出的各损伤分量，λ_i是不同应力状态下累积的等效塑性应变，t表示拉伸应力状态，s表示剪切应力状态，h表示静水压缩应力状态，ρ_c为材料损伤参数，

图3是本发明考虑应变率迟滞效应与未考虑应变率突变效应的对比示意图，本实施例引入时间松弛因子确定当前时间步的结束时刻的等效应变率。1是当在t时刻应变率突变时，应力突然增加；2是应力随应变率的突增而增加缓慢。混凝土是一种率敏感材料，对准确预测强动载下材料的破坏具有重要意义。本发明通过利用动态增强因子DIF径向放大强度面实现应变率效应的考虑，现有应变率效应无法捕捉应变率突变实验下观测到的应力迟滞响应，故引入了时间松弛因子，解决这一问题。即

其中，是瞬时应变率，是是当前时间步开始时刻的等效应变率，η₀是时间松弛因子，η₀＝0.1ms，Δt是时间步，n和n+1是一个时间步的开始与结束。

图4是本发明实施例具体实施的流程图，包括：

A1：基于流体弹塑性力学进行应力更新。

用户可根据弹塑性力学中屈服准则、流动法则和一致性条件进行应力更新推导。

图5是实施例中应力更新过程流程图。具体的当前加载面为

流动势函数为

其中，ω为塑性流动关联因子，当ω＝0时为非关联流动(J2流动)，当ω＝1时为完全关联，当ω介于0和1之间为部分关联。

将式(16)带入式(17)，得

流动法则为

其中，du为一致性参数。

将式(18)代入式(19)，得

其中，Y_,p＝dY(σ_ij,η)/dP。

等效塑性应变增量为

将式(20)代入式(21)，得

等效塑性应变为

λ＝1-(1-δ_sλ_t)(1-λ_s)(1-δ_hλ_h) (23)

建立修正的等效塑性应变和等效塑性应变增量之间的关系：

将式(23)代入式(24)，得(考虑应变率)

假定

将式(25)代入式(26)，得

应变增量分解为弹性和塑性两部分：

将式(20)代入式(28)，得

式(29)两侧同时乘以弹性逆矩阵C^-1

式(30)两侧同时乘以屈服面对应力的梯度

应满足一致性条件(Kuhn-Tucker准则)

其中，Y,λ＝dY(σ_ij,η)/dλ。

将式(26)、(27)、(32)代入式(31)，得

注意到：

则：

将式(34)代入式(33)，得

式(35)提取公因式，移项，得

注意到，在时间步开始的时候泰勒展开另外，则一致性参数：

当前强度面更新

Y_n+1＝Y^*+Y_,pdp^p+Y_,λdλ,dp^p＝KωY_,pdu (38)

更新应力及内变量

A2：基于Fortran语言和Visual Studio进行编程并编译生成LS-DYNA有限元软件所需求解器。

将上述理论推导过程通过Fortran语言表达并编译成可适用于R5版本的LS-DYNA有限元软件的ls971.exe求解器，并将其拷贝至LS-DYNA软件安装目录下的program文件夹下。

A3：基于本发明进行用户自定义工况数值计算并进行结果后处理。

本发明适用于评估侵彻爆炸强动载作用下混凝土材料的损伤破坏，同样适用于混凝土结构的静载力学性能分析。首先，用户要采用国际单位制(m-kg-s)建立模型关键字文件(.k)，然后去更改混凝土材料模型及其参数，最后选择求解器并递交模型关键字文件进行求解即可，结果后处理不变。

下表为本发明实施例中关键字参数输入格式，本发明实施方式已内嵌模量(包括弹性模量、体积模量和剪切模量)、拉伸强度计算公式和状态方程10组参数，本领域技术人员也可以根据经验输入。下表中前两行(MID～IG,IVECT～IEOS)参数选取已固定，其中，MID为材料模型编号，RO为材料密度，MT为用户自定义材料类型，LMC为材料参数数量，NHV为历史变量数量，IORTHO为定义材料各向同性和各向异性属性的参数，IBULK为体积模量参数位置，IG为剪切模量参数位置，IVECT为用户提供向量化子程序的参数，IFAIL为控制壳单元和实体单元材料失效准则的参数，ITHERM为控制计算单元温度的参数，IHYPER为控制计算变形梯度的参数，IEOS为控制状态方程的参数。

表1实施例中关键字参数输入格式

其中，f_c：单轴抗压强度；E：弹性模量，K：体积模量，K＝E/3/(1-2μ)；G：剪切模量，G＝E/2/(1+μ)；T：单轴抗拉强度，ω：塑性流动关联因子，0～1，建议0.5，可根据具体工况标定；frac_t：控制单元拉伸删除阈值，根据实验标定；frac_c：控制单元压缩删除阈值，本领域技术人员可根据实验标定；K₁～K₁₀：体积模量；μ₁～μ₁₀：体积应变；P₁～P₁₀：压力；rate：应变率效应开关，rate＜0时，不考虑应变率效应，rate≥0时，考虑应变率。

(2)本发明的关键点

1.引入三个强度面即屈服强度面、最大强度面和残余强度面，通过线性插值屈服强度面和最大强度面以达到描述混凝土材料应变硬化的目的，通过插值最大强度面和残余强度面以到达描述应变软化(损伤)的目的；

2.根据已有三轴压缩试验数据显示屈服强度面在高围压下与压力关联不大，将高围压下的屈服强度面拉平，用以描述高围压下混凝土材料无明显屈服点的现象；

3.引入基于微裂纹发展引起的拉伸损伤、剪切损伤和基于孔隙坍塌与压实引起的静水压缩损伤，健全了混凝土损伤机制的表达，并引入了拉伸损伤和剪切损伤之间的关联因子，同时在静水压缩损伤引入与应力状态相关的关联因子，以达到表征剪切损伤对拉伸损伤的影响和静水压缩损伤对其他应力状态影响的目的；

4.为解决应力随应变率突变而突变的问题，考虑应变率时引入时间松弛因子，以达到反映应变率突变时应力迟滞效应的目的。

图6是本发明技术方案，引入修正系数来修正高围压段屈服强度面形式；引入拉伸损伤指标、剪切损伤指标和静水压缩损伤指标以及剪切损伤对拉伸损伤影响因子和静水压缩损伤对其他应力状态影响因子，确定等效塑性应变；引入时间松弛因子确定当前时间步的结束时刻的等效应变率。

图7是利用本发明进行的单单元三轴压缩计算结果，本实施例引入修正系数来修正高围压段屈服强度面形式，引入拉伸损伤指标、剪切损伤指标和静水压缩损伤指标以及剪切损伤对拉伸损伤影响因子和静水压缩损伤对其他应力状态影响因子，确定等效塑性应变，通过计算所体现的混凝土高围压下的力学性能。

1是三轴压缩实验结果，2是利用本发明计算的混凝土三轴压缩下的应力应变曲线，3是轴向应力应变曲线，4是应力与径向应变曲线，5是应力与体积应变曲线。可以看出，计算结果与实验结果吻合较好，本发明能较好的考虑混凝土材料三轴力学性能。

图8是利用本发明进行的单单元复杂应力状态计算结果，对应本实施例中引入拉伸损伤指标、剪切损伤指标和静水压缩损伤指标以及剪切损伤对拉伸损伤影响因子和静水压缩损伤对其他应力状态影响因子，通过计算所体现的混凝土高围压下静水压缩损伤累积对混凝土单轴受力状态的影响。复杂应力状态，即混凝土先承受静水压缩，然后释放静水压缩施加单轴压缩。1是低强度混凝土先静水压缩后单轴压缩的实验值，2是利用本发明计算的低强度混凝土先静水压缩后单轴压缩后与未受静水压缩后单轴压缩强度的比值，2是高强度混凝土先静水压缩后单轴压缩的实验值，4是利用本发明计算的高强度混凝土先静水压缩后单轴压缩后与未受静水压缩后单轴压缩强度的比值。可以看出，计算结果与实验值吻合较好，本发明能较好考虑静水压缩损伤对后续应力状态的影响，预测混凝土复杂应力状态性能。

图9是利用本发明进行的单单元动态应变率计算结果，对应本实施例引入时间松弛因子确定当前时间步的结束时刻的等效应变率通过计算所体现的混凝土应变率突变引起的应力迟滞效应。1是应变率为5s^-1单轴拉伸应力应变曲线，2是应变率为10s^-1单轴拉伸应力应变曲线，3是应变率为20s^-1单轴拉伸应力应变曲线，4是应变率由10s^-1突变至5s^-1单轴拉伸应力应变曲线，5是应变率由10s^-1突变至20s^-1单轴拉伸应力应变曲线。可以看出，本发明在计算应变率突变时，应力不发生突变，而是缓慢变化，能较好的反应实验过程中应变率突变引起的应力迟滞效应。

图10是利用本发明进行的混凝土靶体弹体侵彻计算结果，本实施例采用以上技术方案，通过混凝土试件层次的侵彻计算，展示动态拉伸损伤、动态剪切损伤和动态静水压缩损伤的演化过程。图10(1)是利用本发明计算的有限厚度混凝土靶体侵彻结果，并注明了数值计算和试验结果(1中括号内)对比，可以看出，本发明在计算混凝土靶侵彻过程中能很好地预测开坑和震塌现象；图10(2)是利用本发明计算的侵彻过程中靶体的拉伸损伤；图10(3)是利用本发明计算的侵彻过程中靶体的剪切损伤；图10(4)是利用本发明计算的侵彻过程中靶体的静水压缩损伤。可以看出，本发明能较好的预测侵彻过程靶体的损伤破坏情况。

图11是利用本发明进行的混凝土靶体接触爆炸计算结果，本实施例采用以上技术方案，通过混凝土试件层次的爆炸计算，展示利用本发明计算混凝土遭受爆炸荷载作用下的损伤破坏情况。图11(1)是利用本发明针对混凝土薄板的锥形装药接触爆炸计算结果；图11(2)是混凝土薄板的锥形装药接触爆炸试验结果；图11(3)是利用本发明针对混凝土厚板的锥形装药接触爆炸计算结果；图11(4)是混凝土厚板的锥形装药接触爆炸试验结果。可以看出，本发明能较好的预测爆炸作用下混凝土毁伤破坏效应。

以上结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。