基于能量的混凝土弹塑性损伤本构模型

基于能量的混凝土弹塑性损伤本构模型

弹塑性损伤本构的数值实现对于实际存在的许多复杂非线性问题，获得准确的分析解是不可能的。相反，它只能通过算术和改进的数值计算方法（例如有限技术）来寻求类似的解。因此,数值算法是本构模型的重要组成部分,其计算效率的高低极大地影响着本构模型的应用。在非线性有限元分析中,应力和应变之间的非线性关系导致平衡方程是应变(因而也是节点位移)的非线性方程,一般需要采用迭代法加以求解。同时,由于变形历史取决于材料的非线性本构关系,应采用增量分析方法以跟踪位移、应变和外部作用引起的应力的发展。利用增量法求解问题一般采取三种方式:Newton-Raphson方法、修正Newton-Raphson方法和准Newton-Raphson方法。对于存在严重非线性特别是经常涉及材料软化段处理的结构非线性分析,采用切线刚度矩阵的Newton-Raphson方法仍然是最为有效的数值方法。因此,材料非线性本构关系的数值实现实际上包括两方面的内容,即应力更新算法以及非线性方程组的求解算法。为保证计算结果的可信性,迭代过程中采用的应力更新算法必须是无条件稳定的。考虑到不同的计算目的需要,求解非线性方程组宜于采用Newton-Raphson方法,因此,一般还需要给出和应力更新算法一致的切线模量(即算法一致性切线模量)。针对本文第一部分建议的基于能量的混凝土弹塑性损伤本构模型,通过将损伤演化和塑性变形解耦,本文建立了该模型的弹性预测-塑性修正-损伤修正数值分析框架,给出了无条件稳定的应力更新算法和相应的算法一致性切线模量。在弹性预测-塑性修正过程中,利用谱分解回映算法建立了塑性流动因子和硬化参数的统一迭代格式,大大减小了有效应力更新的计算量。根据上述数值算法,作者编制了非线性有限元分析程序,对单调和低周反复荷载作用下的混凝土材料和结构试验进行了数值模拟,分析结果表明:建议的弹塑性损伤本构模型可以较好地描述混凝土材料的典型非线性行为,其数值方法是有效的。上述结论为进一步在结构非线性分析中应用本文建议模型奠定了基础。1应力释放率的修正根据初边值问题的数值解法,所建议本构模型的数值实现需要将时间离散为一系列的增量步[tn,tn+1]∈R+(n=0,1,2,……)进行,即在应变历史t→ε≡∇su(t)及增量步n的{σn,d±n,εpn,κn}的值已知的条件下,更新增量步n+1的状态变量{σn+1,d±n+1,εpn+1,κn+1}。根据运算符分解的概念,模型的应力更新过程可以分为如下弹性预测、塑性修正和损伤修正三个步骤,即(a)弹性预测(Ⅰ)应变更新:在应变场un+1给定的情况下,Gauss积分点处的应变相应更新为εn+1=εn+Δε(1)(Ⅱ)弹性试算应力εp,trialn+1=εpn；κtrialn+1=κn；ˉσtrialn+1=C0∶(εn+1-εpn),dtrial±n+1=d±n(2)(b)塑性修正(Ⅲ)检查屈服条件F(ˉσtrialn+1‚κn){＜0elastic⇒ˉσtrialn+1=ˉσn+1‚转步骤(Ⅴ)＞0plastic⇒回映修正(3)(Ⅳ)塑性回映修正:当有效应力处于塑性加载状态,弹性试算应力和塑性硬化函数都必须“回映”到屈服面上,常用的方法是最近点投影或者截平面算法。一旦加上有效应力空间的塑性一致性条件限制后,状态变量将修正为{ˉσn+1,d±n,εpn+1,κn+1}(4)(c)损伤修正(Ⅴ)计算损伤能释放率:将有效应力张量ˉσn+1分解为ˉσ+n+1和ˉσ-n+1,然后根据以下计算损伤能释放率Y+n+1=√E0(σ+n+1∶Λ0∶σn+1)(5a)Y-n+1=α(¯Ι1)n+1+√3(¯J2)n+1(5b)损伤变量d±n+1更新为d±n+1={d±nG±(r±n+1)若Y±n+1≤r±n其他(6)相应的,损伤阀值r±n+1更新为r±n+1=max{r±0,maxτ∈[0,n+1]Y±τ}(7)(Ⅵ)应力张量σn+1更新为σn+1=(1-d+n+1)ˉσ+n+1+(1-d-n+1)ˉσ-n+1(8)可以看出,建议模型中的损伤变量仅为有效应力的函数,因此上述应力更新可以解耦为弹性预测-塑性修正步中更新{ˉσ,εp,κ}(此时,损伤变量dnn+1保持上一步的值d±n不变)和损伤修正步中更新σ两步来进行。后者较为直接简单,故本文中将主要给出有效应力张量的更新算法。2应力状态的回射式解析非线性计算力学中,一般采用回映算法更新{ˉσn+1,εpn+1,κn+1}。已经有专门的专著介绍采用J2塑性屈服准则的径向回映算法,采用其他屈服准则的塑性力学模型数值算法也已经发展得比较完善。对于本文建议模型中采用的屈服准则,Lee和Fenvas在其塑性-损伤模型中给出了谱分解格式的回映算法。然而,对于平面应力状态,在求解塑性流动因子的迭代过程中,该方法需要插入一级非常复杂的迭代步骤以求解中间变量,大大地增加了数值计算工作量,其算法稳定性和收敛性也无法完全得到保障。同时,在塑性硬化参数的迭代求解过程中,并未考虑塑性一致性函数残量的影响,这可能存在一定问题。针对屈服函数是主应力函数的特点,通过对初始体积模量的简单修正,本文发展了一类适用于任意应力状态的统一谱分解格式回映算法。利用偏量空间和球量空间力学性能的解耦,可以推导证明更新后的有效应力张量ˉσDn+1和试算有效应力张量ˉσtrialn+1具有以下简单关系ˉσn+1=c1ˉσtrialn+1+c21(9)式中:1为二阶单位张量;系数c1和c2分别表示为塑性流动因子增量Δλp的函数c1=1-Δλp2G0∥ˉstrialn+1∥(10a)c2=Δλp((¯Ι1)trialn+13∥ˉstrialn+1∥G0-3αpΚ0)(10b)其中,ˉs为有效应力张量ˉσ的偏量部分;K0和G0分别为材料的初始体积模量和剪切模量;对于平面应力状态,K0与弹性模量E0和泊松比ν0之间的关系为Κ0=1+2ν03(1-ν20)E0(11)于是,如下关系式成立ˉσtrialn+1ˉσn+1=(c1ˉσtrialn+1+c21)ˉσn+1=ˉσn+1(c1ˉσtrialn+1+c21)=ˉσn+1ˉσtrialn+1(12)根据张量和矩阵运算基本原理,(12)式表明对称张量ˉσtrialn+1和ˉσn+1具有相同的特征向量;同时,其特征值矩阵(或有效主应力张量)之间的关系为^ˉσn+1=c1^ˉσtrialn+1+c21(13)由于上述特征值矩阵之间的简单关系,基于后退欧拉法,可以在有效主应力空间内建立塑性流动因子和硬化参数的Newton-Raphson迭代格式,以确定有效应力张量和塑性变形,可以看出,上述谱分解回映算法无需每次重新计算有效应力张量的特征向量,而仅需根据(13)式更新其特征值,大大减少了计算量,提高了建议模型的计算效率。3应力率重塑张量为有效切线刚度退化张量n+1增量步的应力σ(除标明外,本节所有物理量均为增量步n+1更新后的值,为书写简便,省略下标n+1)更新完成后,可以给出与连续体有效切线刚度Cep相对应的算法一致性有效切线模量ˉCalg=dˉσdε。建议模型的算法一致性切线模量按照如下方法得到。应力更新式(8)两边对应变微分可以得到dσdε=(Ι-ω)∶ˉCalg-[ˉσ+d(d+)dε+ˉσ-d(d-)dε](14)式中:四阶对称张量ω为损伤演化引起的有效切线刚度退化张量ω=d+Q++d-Q-(15)Q+和Q-为有效应力率张量˙ˉσ的正、负投影张量,为ˉσ的特征值ˆˉσi和特征向量pi的函数Q+=∑iΗ(^ˉσi)(Ρii⊗Ρii)+23∑i,j>1〈^ˉσi〉-〈^ˉσj〉^ˉσi-^ˉσj(Ρij⊗Ρij)(16a)Q-=Ι-Q+(16b)其中,H(x)为Heaviside函数;〈x〉=xH(x)=(x+|x|)/2为Macaulay函数;Pij为二阶对称张量,表示为:Ρij=12(pi⊗pj+pj⊗pi)(17)注意到˙r±=˙Y±,根据损伤变量的表达式(见文中(39)式)可以得到d(d±)dε=h±dr±dε=h±dY±dε(18)式中:h±称为损伤演化函数,表示为由式(5)可知如下关系成立ˉσ+d(d+)dε=R+∶ˉCalg(20a)ˉσ-d(d-)dε=R-∶ˉCalg(20b)式中:R+和R-为四阶张量,分别表示为R+=h+E02Y+ˉσ+⊗[ˉσ∶Λ0∶Q++ˉσ+∶Λ0](21a)R-=h-[ˉσ⊗(α1+32√3¯J2ˉs)](21b)于是,将(20)式代入到(14)式中,可以给出建议模型的算法一致性切线模量为dσdε=(Ι-ω-R)∶ˉCalg(22)式中:R=R++R-。4混凝土结构非线性分析的数值试验方案为验证建议模型及其数值算法的有效性,作者编写了相应的非线性有限元分析程序,对单调和反复加载情况下的普通素混凝土和高强素混凝土试验进行了数值模拟,并通过一个三点缺口素混凝土简支梁的数值分析说明了建议模型在混凝土结构非线性分析中的应用。除非特别指明,在所有的数值模拟中,均采用二维平面应力等参单元;加载方式按照试验原型施加。混凝土的通用材料参数分别取为:ν0=0.20‚f-b0/f-0=1.16‚αp=0.20‚ˉf+y=ft‚ˉf-y=fc。4.1单轴拉倒测试4.1.1模型参数选取文献进行了混凝土的单轴受拉试验,得到了应力-应变曲线的下降段。试验实测弹性模量为E0=3.8×104MPa,抗拉强度ft=3.40MPa。在数值模拟过程中,模型参数分别取为:E0=3.8×104MPa,f+0=3.68MPa,A+=0.0,B+=0.683,α+E=19.0。数值模拟应力-应变曲线和试验结果对比如图1a所示。4.1.2材料模型参数第二个单轴受拉对比试验是Gopalarantam-Shah试验。该试验实测混凝土的弹性模量为E0=3.1×104MPa,抗拉强度ft=3.48MPa。在数值模拟过程中,模型参数分别取为:E0=3.1×104MPa,f+0=3.57MPa,A+=0.0,B+=0.518,α+E=19.0。数值模拟应力-应变曲线和试验结果对比如图1b所示。从图1可以看出,建议模型可以很好的描述混凝土材料在单轴受拉应力状态下包括软化段在内的应力-应变全过程曲线。4.2单轴压4.2.1弹性模量验证文献进行了两批高强混凝土的单轴受压试验,其中,A类试件的抗压强度稍低,实测弹性模量E0=3.8×104MPa;B类试件的抗压强度稍高,实测弹性模量E0=3.9×104MPa。在数值模拟过程中,A类试件模型参数分别取为:E0=3.8×104MPa;f-0=43.0MPa;A-=1.8;B-=0.75;α-E=19.0;B类试件模型参数分别取为:E0=3.9×104MPa,f-0=62.0MPa;A-=2.0,B-=1.10;α-E=19.0。数值模拟应力-应变曲线和试验结果对比如图2a和2b所示。4.2.2数值模拟结果Karsan和Jirsa进行了普通混凝土的单轴受压试验,实测弹性模量为E0=3.17×104MPa,抗压强度为fc=27.6MPa。在数值模拟过程中,模型参数分别取为:E0=3.17×104MPa;f-0=10.2MPa;A-=1.0;B-=0.16;α-E=19.0。数值模拟应力-应变曲线和试验结果对比如图2c所示。从图2可以看出,无论是对于普通混凝土或是高强混凝土,建议模型均能很好的模拟单轴受压应力状态下包括软化段在内的应力-应变全过程曲线:在确定性本构意义上说,其吻合程度是令人满意的。4.3应力-应变耦合过程利用建议模型对著名的Kupferetal.的双轴试验进行了数值分析,给出了三种应力状态下[应力比分别为(-1)∶0,(-1)∶(-1)和(-1)∶(-0.52)]的应力-应变全过程曲线,其他双轴应力组合下计算得到的混凝土抗压强度结果如表1所示。正如已经指出的:无需引入“强度包络线”、“拉压软化系数”等经验型参数,本文建议模型即可很好的模拟双轴受压强度提高和“拉压软化效应”。4.4低每周重复负荷4.4.1受拉强度验算Taylor在其单轴重复受拉试验中的实测弹性模量为E0=3.17×104Mpa,受拉强度ft=3.47MPa。数值模拟中的模型参数分别取为:E0=3.17×104MPa,f+0=3.50MPa;A+=0.0;B+=8.5;α+E=0.1。数值模拟应力-应变曲线和试验结果对比如图3a所示。4.4.2模型模拟结果分析Karsan和Jirsa单轴重复受压试验的实测弹性模量E0=3.0×104MPa,单轴受压强度fc=27.6MPa。数值模拟过程中模型参数分别取为:E0=3.0×104MPa,f-0=15.0MPa,A-=2.1,B-=0.395,α-E=1.0。图3b给出了数值模拟应力-应变曲线和试验结果对比。从图3可以看出,无论是受拉还是受压情况下,建议模型可以很好的描述塑性变形的影响,刚度退化以及峰值点后的强度软化也与试验结果吻合良好。但需要指出的是,由于弹性卸载-再加载假定,建议模型不能模拟试验中出现的卸载-再加载次滞回环,对于地震动等动力作用下的耗能计算,这将导致一定的局限性。4.4.3受拉受力分析图4给出了混凝土在单轴低周反复作用下的数值模拟结果:先受拉至超过混凝土单轴抗拉强度后(OAB)卸载(BC),然后重新受压加载后(CDE)卸载(EF),然后进行第二次受拉加载后(FGH)卸载(HI),最后进行第二次受压加载后(IJK)卸载(KL)。在数值模拟过程中,模型参数分别取为:E0=3.0×104MPa;f+0=ft=3.0MPa,A+=0.0,B+=1.0,α+E=0.5,fc=30MPa,f-0=13.3MPa,A-=1.0,B-=0.20,α-E=2.0。图4中受拉卸载后受压加载(B-C-D-E和H-I-J-K)曲线表明建议模型可以很好地描述受压后由于受拉裂缝闭合导致的刚度恢复(即单边效应),这对于混凝土结构在地震动作用等反复荷载作用下的非线性分析是相当重要的。同时,卸载(BC和EF)段刚度小于材料的初始刚度(OA)表明建议模型可以较好地描述损伤引起的刚度退化。同样,无论是受拉还是受压,建议模型均可以较好的模拟卸载后可能存在的不可恢复变形。4.5试验验证和分析Peterson曾对集中荷载作用下的素混凝土缺口简支梁进行了详细的试验研究,实测材料属性分别为:弹性模量E0=30GPa,泊松比ν0=0.2,密度为ρ0=2400kg/m3,单轴受拉强度ft=3.33MPa。该试验的典型特征是材料的非线性行为主要由I型受拉裂缝控制,为建议本构模型的结构非线性分析提供了很好的验证标准,不同的研究者对此试验也进行了广泛的理论研究。由于结构和荷载对称性,取半边结构进行非线性分析。经过网格收敛性分析,整个结构共划分为280个四节点、一次减缩积分、平面应力等参单元。考虑到混凝土开裂后可能出现的不