各向异性及多边损伤模型的理论模型

各向异性及多边损伤模型的理论模型

混凝土材料本构建模一般来说，材料结构的微缺陷（如金属和合金材料中的微孔、微切割带、高聚物材料中的微孔和银纹等）的产生、发育和融合是材料最终破坏的主要因素。对于混凝土和其他准刚性材料，在初始状态上可以认为是各向同性。在外部负荷的作用下，在某些情况下，在嘉靖和维护过程中形成的初始微缺陷沿着最大的尽量服从方向发展，形成了一条微裂缝。同时，随着应力历史的变化，这些微裂缝的发展方向也会发生一定程度的旋转和一些微裂缝，导致材料的明显的不同向异性。所谓的破坏引起的不同方向的异向异性。此外，当微裂缝通过弹性作用时，材料的刚性变化很大，柔度也会增加。相反，在压力作用下，微裂缝关闭（即损伤受抑制），材料的刚性基本恢复。这种由微裂缝打开、闭合和重新开口过程引起的刚度变化，即所谓的“单边效应”。上述微裂缝（损伤）引起的材料的各向异性和一侧效应对混凝土材料和结构的变形和浮力产生了非常重要的影响，如重复加载和非比例加载。自从Kachanov最早提出损伤力学的基本思想并用于分析金属蠕变问题以来,损伤力学已经逐渐发展成为包括连续损伤力学和细观损伤力学等在内的一门相对独立的固体力学分支.其中,尽管近年来细观损伤力学发展非常迅猛,但通过引入标量、矢量、二阶或四阶甚至更高阶张量形式的损伤内变量,连续损伤模型能够在宏观上反映微裂缝演化对材料非线性行为的影响,在混凝土材料本构建模领域得到了非常广泛的应用[14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].从物理意义上讲,“损伤”表征为材料的刚度退化,因此,采用刚度张量(或柔度张量,或其增量)作为损伤变量建立损伤模型更为合理.基于类似的考虑,文献引入了退化应变率的概念,并为众多学者所采用,随后Carol等归纳了弹性退化(损伤)理论;Govindjee等和Meschke等分别将该理论推广至多面弹性损伤和弹塑性损伤模型.在连续损伤力学基本理论框架内,通过将某些二阶张量如应力、应变或有效应力等进行正、负分解并引入相应的四阶投影算子,已有部分混凝土损伤模型[4,16,17,18,20,26,27,28,29]能够描述单边效应的影响.然而,研究表明:对于各向同性损伤,已有模型往往可以给出令人满意的预测结果,而在描述各向异性损伤时,这些模型会给出损伤卸载能量耗散不为零等热力学不相容结论.此外,由于投影算子的引入,已有模型的数值算法过于复杂或无法线性化而计算效率十分低下,大大限制了模型的应用.基于不可逆热力学原理和内变量理论,本文直接以材料柔度张量的增量作为损伤内变量,建立了热力学相容的各向异性单边损伤模型的一般形式,给出了热力学相容的投影算子,推导了模型的率本构关系并发展了其数值实现算法和算法一致性模量.在此基础上,建立了合理的损伤准则和演化法则,并将上述模型应用于混凝土材料.数值模拟结果初步验证了模型的有效性.建议模型完全构筑于热力学基础上,无需引入应变等效或应变能等效等经验性假定.需要说明的是,模型不考虑动态载荷下混凝土材料的应变率效应,其率相关推广将在后续研究中给出.此外,由于不可恢复变形可以通过塑性力学方法加以描述,本文将集中讨论微裂缝演化(损伤)对材料非线性行为的影响,因此,建议模型隶属于弹性损伤范畴,但推广至弹塑性损伤模型不存在原则性的困难.为方便后文表述,二阶和四阶单位张量分别记为1和Ⅰ,此外,引入以下张量运算符其中,A和B均为二阶对称张量,X和Y均为四阶张量,i,j,k和l为分量下标.1热容量单级弹性损伤模型1.1微裂缝演化的应力-应变耦合关系传统的连续损伤模型中,一般选取某些标量或二阶张量作为损伤内变量,其选择具有较大的随意性.同时,为了描述损伤演化对材料宏观力学性能的影响,通常还需要引入诸如应变等效或应变能等效等假定.Krajcinovic从原子、微细观和宏观等尺度综合分析了材料细观结构非均匀特性和微裂缝非线性演化特征对材料宏观力学行为的影响,结果表明:从一般意义上讲,作为沟通微细观和宏观尺度的桥梁,选取材料的柔度或刚度张量作为损伤内变量更为合理.对于初始各向同性材料,其线弹性柔度张量C0和刚度张量S0分别表示为式中,E0和v0分别表示材料的弹性模量和泊松比.出现损伤后,微裂缝演化导致材料刚度S发生退化,而柔度张量C则相应增加.若将损伤变量统一标记为B★(可以为标量、二阶张量或四阶张量,下标“★”表示张量的阶数),则可以给出如下关系式式中,Λ表示柔度张量的增量,为损伤内变量B★的非线性函数.注意到微裂缝演化与拉应力的大小和方向密切相关,将应力张量σ分解为其正、负分量之和的形式,即其中,正、负应力分量σ+和σ-分别表示为式中,σ(i)为应力σ的第i个特征值,相应的特征向量为n(i);<·>为McAuley尖括号,定义为<x>=max(0,x).引入四阶对称投影算子P±,将式(4)写为对于弹性损伤材料,其卸载路径直线指向原点,而且在整个卸载—再加载(再加载指再加载至先前的卸载点,后文同)的过程中,材料的损伤状态保持为当前的值不变,即=0.此时,弹性损伤材料的Gibbs自由能势ψ依赖于当前的应力σ和损伤状态B★,即其中,柔度张量的增量Λ表示为式中,Λ+和Λ-即为损伤内变量B★,其物理意义为分别在σ+和σ-单独作用下,微裂缝演化引起的材料柔度增量.局限于纯力学范畴(等温、绝热状态),任意可能的变形过程均需满足不可逆热力学第二定律,即如下Clausius-Duhem不等式将式(6)对时间微分并代入到式(8)给出由的任意性,可以给出如下割线应力-应变关系式将式(10)两边对时间微分,可以得到式中,为阻止微裂缝继续发展的弹性应变率即损伤不进一步发展(保持上一步的柔度C或刚度S不变)的前提下,应力的增量引起的应变,而则为损伤应变率即在当前应力水平下,由于刚度退化(或柔度增加)而导致的应变,也称为退化应变率.其中,根据式(7)率张量Λ表示为上述损伤应变率最早由Hueckel和Maier给出,随后被文献采用.1.2边界条件12上述本构关系的推导过程中,忽略了其中的关键步骤:既然依据热力学基本原理给出,所得到的本构关系成立与否需要验证其是否满足不可逆热力学的基本要求.对于弹性损伤模型而言,式(10)成立的前提是式(6)中定义的函数ψ确实是有效的Gibbs自由能势.上述要求实际上等效为:当材料损伤保持恒定不变即=0时,如下伪能量耗散率应该为零,即当式(10)成立时,式(13)给出的伪能量耗散率简化为对于Green弹性材料,当不考虑单边效应时,在整个损伤卸载—再加载过程中损伤状态不变即意味着其柔度张量保持当前值不变,即,此时上述热力学要求实际上等效于本构关系式(10).然而,对于上述单边弹性损伤模型中引入了投影算子P±且其表达式并不唯一,从式(14)可以看出=0并不一定可以保证.显然,当各向同性损伤或当P±保持固定不变即=0时(如比例加载等),的条件自动得到满足;而对于一般情况,投影算子P±满足能够保证热力学相容条件=0.综合考虑式(5)和式(15),可以得出结论:只有同时满足如下条件的投影算子才是热力学相容的,式中,的率张量.附录表明,热力学相容的投影算子唯一地表示为如下表达式式中,H(σ(i))为σ(i)的Heaviside函数,二阶对称张量N(ij)表示为1.3损伤演化法则对于任何可能的变形过程,根据不等式(9)和式(10)可知,本构关系必须满足如下耗散不等式式中,-Y±为与损伤变量Λ±功共轭的热力学广义力即损伤能释放率.从式(19)标准的二项式结构可以看出,应该为四阶对称正定张量.需要指出的是,上述式(14)和式(19)两个看起来类似的条件实际上是完全独立的:前者限制了卸载一再加载过程中材料损伤行为,而后者则规定了损伤加载演化法则应该满足的条件.构建合理的损伤变量演化法则是发展连续损伤模型的最大难点之一.由于混凝土非线性行为的复杂性,对于一般的各向异性损伤模型,给出高阶损伤变量的演化法则更加困难.已有连续损伤模型的演化法则往往只能基于唯象学给出,而目前细观损伤理论微裂缝的方向和大小分布函数往往依赖于人为假定,且难以有效地考虑软化段大量微裂缝的相互作用.基于上述考虑,本文类似于塑性流动法则,给出如下损伤演化法则的一般形式式中,λ±≥0为损伤演化因子;损伤演化方向ψ±为正定的四阶对称张量,以保证损伤耗散式(19)为非负值.上述损伤演化法则最早由Simo和Ju提出,并得到了较为广泛的应用.由于热力学相容的投影算子满足关系式(15),相应的损伤应变率可以简单地表示为式中,损伤应变的演化方向Γ±为需要指出的是,一旦确定了柔度增量的演化法则,损伤应变的演化法则即可根据式(21)给出;反之则不成立.损伤演化因子λ±与损伤状态有关:当处于损伤加载时,λ±>0,而λ±=0则表明损伤卸载或处于弹性区域内,因此,必须首先根据损伤准则确定相应的损伤状态.本文采用如下损伤准则和如下损伤加载/卸载条件式中,损伤加载函数f±与当前的应力状态有关;损伤阈值r±为累积损伤λ±的函数.当处于损伤加载即时,可以通过损伤一致性条件给出,即其中,Υ±为损伤加载函数f±的应力梯度,而h±则表示硬化/软化函数.1.4弹性损伤模型的建立为了确定材料的失效模式和局部化状态,往往需要给出材料的连续体切线模量.此外,连续体切线模量的推导过程往往对于发展本构模型的数值算法有一定的启发作用.对于建议模型,考虑到式(11)和(21),可以给出如下率本构关系式中,损伤演化因子±和相应的连续体切线模量Stan可以分为以下几种情况求得:(1)当F+<0,F-<0时(2)当,F-<0时(3)当时(4)当时其中,系数Hij和Δ分别表示为从上述率本构关系和切线刚度张量清楚地看出,由于引入了损伤应变率且采用了与塑性流动法则类似的损伤演化法则即式(20),只需要将经典的塑性力学基本公式中的初始刚度张量S0替换为割线刚度张量S,即可给出建议的弹性损伤模型.然而,此处的割线刚度张量S为内变量,其演化法则根据关系式S=C-1和式(20)给出的演化法则加以描述.2实现模型数值的算法2.1应用效果及改进的双面性损伤模型的健全本文建议的各向异性单边损伤模型率本构方程的数值积分实际上是对状态变量进行更新的过程.若整体非线性方程采用基于应变的有限元进行求解,一般将时间离散为若干子增量步[tn,tn+1]R+(n=0,1,…),初始条件为,且应变历史即t→ε≡▽su(t)已知.根据无条件稳定的后退欧拉方法,时间步tn+1的状态变量按下式更新式中,为弹性试算应力,表示为于是,式(33a)～(33d)等效为如下残量形式可以看出,式(34)和(36)实际上构成了模型数值实现的弹性预测——多面损伤修正方法:在损伤加载状态,弹性预测应力和其它状态变量通过多面损伤修正“回映”到当前损伤面.典型的回映算法有最近点投影方法和切平面方法两种.由于可以保证二阶收敛速度并能够给出与算法一致的线性化方程,本文采用最近点投影方法.相应地,将式(36)线性化为(此处,为了书写简便,省略了下标n+1,下同)根据式(37a),δσ可以用δ(Δλ±)表示为式中,四阶对称张量Π称为弹性算法模量,表示为将式(38)代入式(37b)中,可以得到当两个损伤面都处于损伤加载时,式(40)可以写成如下矩阵形式式中,矩阵J的分量Jij表示为因此,第k+1次迭代后的Δλ±可以求解为其中(·)(k+1)和(·)(k)分别表示第k+1次迭代和第k次迭代后变量的值.当仅有损伤面F+或F-处于损伤加载状态时,式(40)退化为单面损伤的形式,于是得到Δλ+和Δλ-或根据式(43)～(45)计算出损伤演化因子的增量δ(Δλ±)后,其它物理量可以相应地进行更新,不再赘述.在已有各向异性单边损伤模型中,几乎都无法给出相应的率本构关系和连续体切线模量,其非线性有限元数值实现算法则采用显式算法,或只能够采用效率较低的准Newton-Raphson算法或拟NewtonRaphson算法,而这些方法对于混凝土此类强非线性材料而言往往不够强大.Mahnken等发展了Ortiz模型的隐式求解算法,然而其过程异常繁琐,所得结果也极其复杂,效率不高而极大地限制了其应用.本文建议模型不仅在理论上解决了已有模型热力学不相容的理论缺陷,同时还给出上述无条件稳定的隐式数值实现算法.2.2计算力学中的结论为了保证建议模型在有限元数值实现时具有二阶收敛速度,在非线性方程组的迭代求解过程中,还必须使用与上述最近点投影算法一致的切线模量dσn+1/dεn+1.为此,将式(33)线性化,可以得到综合考虑式(46a)和(46b),将得到如下关系式于是,将式(47)代入到式(46c),可以将dλ±用dσn+1表示,并将其结果代回到式(47),即可给出与算法一致的切线模量dσn+1/dεn+1表示为:(1)当时(2)当时(3)当时(4)当时式中,detJ为矩阵J的行列式即可以看出,将第1.4节给出的连续体切线模量Stan中的割线刚度S替代为弹性算法模量Π,即得到了与算法一致的切线模量dσn+1/dεn+1.因此,当Δε→0也即时,算法模量趋近于割线刚度张量即Πn+1→Sn.相应的,与算法一致的切线模量退化为连续体切向刚度张量,即dσn+1/dεn+1→.上述结论与计算力学中的结论是一致的.3混凝土应用3.1损伤阈值的确定和工作原理在缺乏试验数据的情况下,通常很难直接给出柔度增量Λ±完全各向异性损伤的演化法则.受到各向同性和正交各向异性损伤演化法则的启发,本文建议如下柔度增量的演化方向张量相应地,材料柔度增加(或刚度退化)引起的损伤应变率表示为式中,E0和v0分别为材料的弹性模量和泊松比;为σ±的单位张量.可以看出,当v0分别取为v0=0和v0=-1时,式(54)分别退化为Ortiz和Losi等提出的演化法则.需要指出的是,Ψ±为正定的四阶对称张量,满足损伤耗散为非负值的热力学要求.研究表明,对于率无关材料,可以根据损伤耗散不等式(19)给出如下损伤加载函数f±表达式上述损伤加载函数f±没有考虑损伤演化对材料泊松比的影响.试验表明,拉应力对混凝土泊松比的影响不大,而在双轴或三轴受压应力下其泊松比会逐渐增加,最终破坏时将超过0.5甚至接近于1.0.基于上述事实,本文引入参数v±以反映材料泊松比的变化,并考虑到损伤阈值r±归一化的需要,模型建议采用如下损伤加载函数式中,四阶对称张量Θ±表示为与损伤演化张量Ψ±相同的形式即相应的,损伤面的应力梯度Υ±分别表示为本文取参数v+为混凝土的泊松比v0,即v+=v0,此时Θ+=E0Ψ+,则建议的损伤加载函数f+与式(56)是等效的.若将单轴受拉和等双轴受拉强度分别标记为ft和fbt,则损伤阈值r+的最大值表示为一般来说,混凝土的泊松比v0在0.15～0.25之间,由此可以得到fbt/ft的值在0.767～0.816之间,这与试验结果非常吻合.例如,泊松比v0=0.20时,相应的,fbt=0.79ft.类似的,将单轴受压和等双轴受压强度分别标记为fc和fbc(均为正值表示),则对应于上述强度的损伤阈值表示为对于典型的混凝土材料,强度比值fbc/fc一般在1.10～1.20之间,因此,v-的取值范围为0.587～0.653.本文中,取fbc/fc=1.16,相应的,v-=0.6384.将建议模型和上述损伤加载函数应用到一维应力状态下,即可得到适合于混凝土材料的损伤阈值函数r±.本文中,r+和r-的函数形式分别取为相应的,硬化/软化函数h+和h-分别表示为式中,则分别代表单轴受拉和单轴受压的线弹性极限强度,可以表示为单轴抗拉强度ft和单轴抗压强度fc的如下函数分别为控制单轴受拉和单轴受压曲线的硬化/软化形状的模型参数,与材料属性断裂能和单元特征长度有关.为了获得与网格大小无关的计算结果,当本文采用裂缝带理论对软化段进行规则化,即假定单轴受拉和单轴受拉应力-应变曲线与应变轴所围的面积应该分别等于单位特征长度lch的抗拉断裂能/lch和抗压断裂能/lch.将上述理论应用于建议的混凝土各向异性损伤模型,模型参数与材料的抗拉断裂能和抗压断裂能之间的关系表示为式中,a+和a-分别定义为不同a+和a-取值条件下,归一化的损伤阈值r+/ft和r-/fc的演化曲线如图1所示.3.2抗拉受力分析众所周知,混凝土是一种典型的具有初始微缺陷(微裂缝和微空洞)的准脆性材料.试验观察表明,在很低的应力水平下这些微缺陷即开始发展、融合、贯通和成核,而且,其损伤演化与应力历史密切相关且具有本质意义上的各向异性特性.本节中,将建议模型应用于混凝土典型非线性行为的数值模拟,其中,采用的材料属性分别为E0=3.0×104MPa,v0=0.20,ft=3.0MPa,fc=30.0MPa以及fbc/fc=1.16.首先考虑单轴受拉(σ1>0,σ2=σ3=0).此时,除了分量=λ+/E0>0外,柔度增量Λ+和Λ-的其它分量均为零.于是,利用建议模型,可以得到加载方向的应力-应变关系为式中,单轴受拉应力状态下,相对应变表示为其次,对于单轴受压(σ1=σ2=0,σ3<0),除了分量=λ-/E0>0外,Λ+和Λ-的其它分量为零.此时,加载方向的单轴受压应力-应变关系曲线可以表示为式中,参数与a-的关系由式(65)给出,单轴受压应力状态下的相对应变表示为求出上述单轴受拉和单轴受压应力-应变曲线与应变轴包围的面积,根据裂缝带理论即可得到参数a±的表达式(66).不同参数值a±给出的单轴受拉和单轴受压应力-应变关系曲线如图2.图3(a)给出了单轴反复拉压应力状态下加载方向x的应力σ1=σxx与该方向的应变值ε1=εxx的关系曲线(O→A→B→O→C→D→O→B→E→O).图中模型参数取为a+=0.50和a-=6.00.从中可以看出,建议模型能够较好地描述微裂缝演化引起的刚度退化和单边效应.最后,将建议模型应用于混凝土的双轴应力状态下(σ3=0),得到的强度包络线如图3(b),图中还给出了Kupfer等试验结果.数值模拟结果表明建议模型能够较好地描述双轴受压时的强度提高和双轴拉压应力状态下强度的降低.4数值模拟结果已有的正交各向异性或完全各向异性损伤模型存在热力学不相容、数值算法复杂且效率低等缺点.基于连续损伤力学和不可逆热力学原理,本文直接以材料柔度张量的增量作为损伤内变量,建立了各向异性单边损伤模型的一般形式,并给出了热力学相容的投影算子,推导了模型的率本构关系.在此基础上,发展了高效、稳定的Newton-Raphson数值实现算法及其算法一致性模量,并给出了适用于混