基于一类非局部宏−微观损伤模型的裂纹模拟1)

2020-06-10卢广达陈建兵

力学学报 2020年3期

卢广达陈建兵

(同济大学土木工程防灾国家重点实验室,土木工程学院,上海 200092)

引言

工程结构的损伤、破坏和解体往往是微裂纹形核、萌生、发展和传播扩散的结果.因此,准确把握固体材料与结构中裂纹的产生和发展过程至关重要[1].自20 世纪20 年代初Griffith 的开创工作[2]以来,人们对此进行了大量的研究,并先后发展了线弹性断裂力学[3-6]、弹塑性断裂力学[7-11]和损伤力学[12-18].然而,固体中的裂纹萌生与传播的模拟,迄今仍面临一系列挑战性问题.例如,经典断裂力学难以解决裂纹形核以及扩展方向等问题[19].20 世纪末,通过引入能量释放率最小化等准则和裂纹面密度连续化思想,基于变分原理发展了脆性材料的相场理论,为裂纹扩展模拟开辟了新的道路[20-21].最近,结合损伤力学和相场理论,Wu 等[22-25]提出了一类统一相场理论框架,对准脆性材料的静力与动力裂纹扩展分析问题取得了重要突破.这一途径,从计算角度来看,属于弥撒型裂纹模型的范畴.在传统上,人们认为断裂力学与损伤力学是分别处理和分析裂纹出现之后与产生之前的两种力学理论.统一相场模型的成功,从一个侧面表明,这一认识存在局限性.事实上,二者之间是可以相容而相得益彰的.与此同时,人们意识到,仅仅在宏观连续介质力学尺寸不足以把握材料损伤萌生与裂纹扩展的物理−力学本质,要从根本上解决上述问题,需要结合微细观机理的分析,由此发展了宏微观断裂力学[1].另一方面,人们逐步认识到非局部化性质对克服经典局部损伤力学中网格敏感性问题的根本性意义[26].鉴于裂纹本质上是不连续的,20 世纪末以来人们提出了一系列新的途径,包括内聚裂纹模型[27]、扩展有限元方法[28]和无网格方法[29].自2000 年以来,Silling 提出peridynamics(近场动力学[30]或毗域动力学[31])理论[32],将固体介质的相互作用采用积分和离散形式表达,给出了一类新的、具有非局部化性质的离散裂纹萌生和扩展模拟方法,并在一系列问题中取得了相当的成功[33-37].

虽然如此,统一相场理论需要额外求解相场演化方程,从而使得理论和求解的复杂程度较高.另一方面,近场动力学方法中本构模型的确定以及对准脆性裂纹的模拟尚存在困难.为此,结合统一相场理论和近场动力学的基本思想,针对准脆性材料的裂纹模拟,最近提出了一类新的非局部宏−微观损伤模型[38].在此基础上,本文进一步修正微观损伤准则,并对具有强非线性回弹特性的典型裂纹扩展问题进行分析,验证了所改进模型的有效性.

1 宏−微观损伤模型的基本理论

从微细观层次来看,裂纹的萌生和发展从几何上表现为出现了不连续性,从宏观上虽依然可认为是连续体,但与原初固体材料相比,出现了宏观损伤.这一损伤将引起材料自由能的变化,使得本构关系中出现非线性,因而导致材料表现出非线性行为.本节将给出这一基本路径的理论刻画[38].

1.1 拓扑损伤

考虑连续固体B,其边界为∂B.该连续体中的物质点记为x=(x1,x2,x3)∈B.考察物质点对(x,x′),其中∈B,ξ=x′−x为两点之间的空间差,则为相应的距离,这里为向量的长度.连续体发生变形后,物质点x和x′的位移分别记为u=u(x)和u′=u(x′).

根据变形几何(如图1 所示),可以定义物质点对(x,x′)之间的变形为

其中,ν=为物质点对的单位方向向量.由此,可定义物质点对的伸长量为

式中,H(·)为Heaviside 阶跃函数,当λ0 时,H(λ)=0,否则H(λ)=1.

图1 物质点对的变形Fig.1 Deformation of material point pair

可以合理地假设,当伸长量超过给定阈值¯λ 时,物质点对之间的联系(物质键) 开始弱化,直至完全分离,即物质键断裂.这一过程显然与单调非减的加载历史参数有关.为此,可定义加载历史参数为即物质键变形历史的最大超越伸长量,其中t为时间变量(对静力问题应理解为虚拟时间).

注记1物质键[38]的观念受启发于近场动力学方法[32-34],类似的概念也见于分子动力学[39]和虚内键模型[40]中.需要指出的是,在本文中是临界伸长量,而不是文献[38]中的临界伸长率.

不难验证,历史最大超越伸长量有dtκ0,这里,后续符号类似.因此,对于任意物质键(x,x′) 可以定义一个刻画其联系强弱程度的量ω(x,x′,t).不失一般性,可以对其进行归一化,即ω ∈[0,1].它是历史最大超越伸长量κ 的单调非减函数,即

注记2从物理意义上看,ω 是微细观损伤的度量,其具体形式可由分子动力学的相互作用势函数通过计算加以确定,或者根据理性假设给出.在本文中,假设为如下形式

其中γ>0 为模型参数,它表征微细观损伤的发展速度.不难注意到,是以γ 为参数的、Heaviside 阶跃函数H(κ)的极限函数列,而Heaviside 阶跃函数表征了不连续性.因此,从这个意义上说,微细观损伤刻画了固体在物质键层次的不连续程度,这一点与李杰等[18]认为损伤变量应打破物质空间连续性的观点是一致的.

事实上,连续体中某物质点处的损伤,往往不仅取决于该点本身,还必然与其周边一定邻域范围内的物质点相互作用有关.换言之,其损伤的变化,本质上是非局部化的.记存在影响的邻域(影响域)半径为ℓ,即ℓ=diamDℓ(x)/2,这里Dℓ(x)表示物质点x的影响域.根据上述表述,可以进一步提炼出损伤的数学语言:

定义1设连续体B⊂Rd(d=1,2,3)及其位移场u⊂Rd,若给定影响域半径ℓ > 0,存在阈值¯λ > 0,使得物质点x∈B有

则称连续体在x处发生损伤,此时微细观损伤ω>0.它从数学上刻画了位移场u在空间上关于方向ν的不连续程度,因此不妨称这一定义为¯λ −ℓ 语言.

作为对照,这里给出位移场u作为d元d维向量值函数其连续性的ϵ −δ 语言[41].

定义2设u:B⊂Rd→Rd,若对任意ϵ > 0,存在δ>0,使得x∈B有

其中Dδ(x)=,则称向量值函数u在x处连续.

注记3定义1 和定义2 恰好构成关于函数不连续性和连续性的一组对偶表述,而这要得益于本文采用物质键伸长量表征微细观损伤的改进.需要说明的是,这一损伤的−ℓ 语言是文献[38]中尚未认识到的,并区别于经典连续损伤力学理论[13,17]中损伤变量只表征材料刚度退化的观念.

对∀x∈B物质点给定的影响域Dℓ(x),设任意两物质点之间的影响函数为ϕ(x,x′),并满足如下条件

从宏观上看,对于固体B内任意物质点x,其损伤应是该物质点所连接的物质键之微细观损伤的加权平均.因此,对∀x∈B定义宏观层次损伤为

即根据微细观层次的物质键(x,x′) 的断裂程度在局部空间上进行几何上的加权平均,它连续地刻画了位移场在一点附近的整体不连续性.不难注意到,损伤区域JΩ={x∈B:Ω(x)>0}构成了裂纹拓扑J(其中dimJ=d−1)的非局部化邻域表达.因此,从这个意义上,不妨称之为拓扑损伤.

注记4拓扑损伤Ω 虽然在形式上与光滑处理的非局部损伤[42]相同,但二者存在着本质的区别:前者积分项内的损伤变量是建立在点对(两点)意义上的,而后者则是基于一点的.需要指出的是,基于一点的非局部损伤已经被证明存在应力锁死的问题[43],而拓扑损伤则不存在此问题[38].特别地,当微细观损伤ω 取为Heaviside 阶跃函数时,则式(10)退化为近场动力学的损伤表示[33].

注记5事实上,式(10) 以积分形式定义的拓扑损伤与统一相场理论[22]中含梯度算子的椭圆泛函[44]

的作用相同,即采用非局部化思路将裂纹拓扑表示为含一定带宽的弥散型界面,从而使得分析过程无需人为预设裂纹扩展路径.其中,式(11)的参数φ 和α 分别为相场函数和裂纹几何函数.

进一步地,可以证明以下命题:

命题1设裂纹拓扑J⊂Rd,以及Ω 和Θ 分别如式(10)和式(11)所定义,则

式中,ϑ(d)=Υ(d)/2,其中Υ(d) 表示d维单位球域Dℓ=1(0)的表面积测度.

需要说明的是,上式第一个等式(即椭圆泛函逼近)已经由Braides[44]给出了证明,而第二个等式(拓扑损伤逼近) 的具体证明因限于本文的篇幅而考虑另文给出.

1.2 能量退化因子

注意到,连续体的自由能密度必然随着拓扑损伤的发展(即裂纹萌生和扩展)而退化.设能量退化因子为g,则有损连续体的自由能为

式中小应变张量ε=∇Su,其中∇S(·) 表示对称梯度算子,ψ0为无损材料的自由能密度函数,即

其中E为四阶弹性张量.

为简化计,在本文中暂不考虑塑性的影响.因此,这里弹性应变和总应变相同,后文对二者将不再加以区分.

显然,能量退化因子与拓扑损伤有关,因而它是拓扑损伤的函数g=g(Ω),表征了物质点储能能力的退化,称为能量退化因子[45]且满足如下要求[22]

其中,g(0)=1 表示物质点处于线弹性无损状态;g(1)=0 则是完全损伤的刻画,即丧失了全部的储能能力; 此外,dΩg< 0 表示损伤总是耗能的; 而dΩg(1)=0 则保证当Ω=1 时,损伤驱动力∂Ωψ=dΩg(Ω)ψ0趋于零,即当一点的损伤发展完全后,其驱动力应随之消失.

若从能量意义上说,不妨称Φ=1 −g(Ω) 为基于能量的损伤.事实上,基于拓扑损伤的能量退化因子还可以进一步从物理上论证,g是[0,1] 上的凸函数[38],其定性的物理解释如下.

由于拓扑损伤Ω 考虑的是影响域内所有方向上物质键损伤的几何“平均”,而各个方向上的物质键并非同时地发生损伤,因为物质键的伸长量在各个方向上的分布一般是不均匀的.在初始阶段,只有位移场在空间变化剧烈方向上那些少数(相对于所有方向上的物质键总数而言) 超过阈值的物质键率先进入损伤状态,此时拓扑损伤Ω 在数值上非常小,但这些率先发生损伤的物质键因变形剧烈而储存了较大比重的变形能.因此,在Ω=0 附近处基于能量的损伤Φ 要大于拓扑损伤Ω,即随着其他方向物质键损伤的发展及其能量的释放,逐渐地一点邻域内所能释放的变形能比重越来越小,即能量退化因子g随着拓扑损伤Ω 的增大而趋于平缓.

根据上述解释,可以清楚地看到:在本模型中能量退化因子是拓扑损伤的非线性函数,不同于经典连续损伤力学[13,17]1 −Ω 的线性形式.借鉴于统一相场模型[22]的结果并加以适当修正,可取能量退化因子为如下凸函数[38]

其中,p和q均为退化函数.一些代表性取值的函数曲线如图2 所示.

图2 能量退化因子Fig.2 The energetic degradation factor

2 控制方程

2.1 本构方程

根据Clausius-Duhem 不等式[17],有能量耗散率dtΞ 不等式

其中,σ为Cauchy 应力张量.将式(13) 代入上述不等式并稍作整理,得

由前述拓扑损伤即能量退化因子的定义和性质可知,上述积分项内第二项必然不大于零.因此,该不等式成立的充要条件是

式中,σ0为有效应力张量.

注记6由此可见,能量退化因子表征了材料的物性关系,式(20)与传统损伤力学理论[18]的本构关系在形式上完全相同.其中,本文式(17)退化因子给出的是准脆性材料本构,p和q的取值刻画了材料的宏观脆性程度,可认为是物性参数,应根据标准试件进行标定.

2.2 平衡方程与边界条件

平衡方程和边界条件与常规连续介质力学的基本表达相同.考虑图3 所示固体B的参考构形,容易列出平衡方程为

其中,t为内力向量,为体力密度向量(当不计体力时,=0).

图3 固体的构形及其边界Fig.3 The solid configuration and its boundary

根据Cauchy 假定,即t=σ·nP,其中nP为子集P的边界外法向向量.由式(20)所表述的本构关系,利用散度定理并注意到P的任意性,可得

而在边界上则分别有

其中,nB为固体的边界外法向向量,和分别为给定的边界面力和边界位移.

注记7得益于式(10)所定义的拓扑损伤关于位移场是Lipschitz 连续的(具体证明见文献[38]),本文模型仍然可以采用连续介质力学的Cauchy 假定建立平衡方程.这一点区别于近场动力学积分形式的平衡方程中关于应力散度项的低阶逼近[46]

其中,f和T分别为近场动力学内力密度和内力状态向量.尽管通过引入高阶微分算子[47]可以一定程度提高其逼近精度以及规避因内力作用方式重构而带来的零能模式问题[37],但也增加了问题求解的复杂性和计算量.此外,近场动力学方法本来的计算效率就低于有限元方法[48].

3 数值求解方法

3.1 有限元离散方程

上述问题,可以直接采用常规有限单元法进行离散和求解.式(22)和式(23)这一边值问题经过弱形式逼近后,可得单元的平衡方程为

式中,Ue为单元结点位移向量,Ke为单元刚度矩阵,Re为单元结点力向量,其中

这里Ae为单元区域,∂tAe为边界单元的力边界,Be和Ψe分别为几何矩阵和形函数矩阵,其表达式可在一般有限元专著中找到,兹不赘述.此外,De为弹性矩阵,对于平面应力问题有

其中,E为弹性模量,µ为泊松比.

根据单元组集规则,可得结构整体平衡方程为

其中,整体刚度矩阵K、整体等效载荷向量R和整体位移向量U分别为

式中,Λe为单元转换矩阵,Ne为单元数目.

3.2 单元刚度矩阵的计算

在式(26)中,为了计算单元刚度矩阵,需要计算给定积分点的能量退化因子g.这时,应通过离散化方式计算式(10) 中的拓扑损伤Ω,然后将其代入能量退化因子式(17)中.此时,拓扑损伤可根据离散化后的各单元变形计算.对于常应变单元,有

其中,xe和xe′分别为单元Ae和Ae′的形心坐标,而=vol[Dℓ(xe)∩Ae′].若考虑采用等参单元进行离散,则二者为相应的高斯点及其权重体积.

注意到,能量退化因子是位移的函数.由式(26)可知,单元刚度矩阵也是位移的函数,因而由式(30)给出的整体刚度矩阵亦为位移的函数.因此,平衡方程(29)是一个非线性方程.在本文中,采用改进局部弧长法进行求解.限于篇幅,兹不赘述.

注记8在求解中,拓扑损伤Ω 是整体节点位移向量U的“显式”函数,即一旦求得整体节点位移向量U已知,直接将其代入式(31)就可以得到拓扑损伤Ω,整个过程自始至终只需要求解(29) 的平衡方程,无需额外的待求解方程.这一点比相场模型简单而高效,因为相场模型中相场函数和位移场均是待求未知量,需要采用联立或交替的策略求解平衡方程和额外的相场演化方程[49].

4 典型实例分析

以Ingraffea 和Grigoriu 所做的有机玻璃剪切梁试件[50]为例,来验证和分析本模型的有效性.注意到,该试件以其强非线性回弹特性而著称,是检验数值模型优良性的经典算例.近些年来,经过该算例检验的数值模型有基于混合(应变/位移)元的各向同性损伤模型[51]和统一相场模型[23]等.分析中试件按平面应力问题处理,其材料参数按文献[51]取:弹性模量E=3.102 GPa,泊松比µ=0.35.宏−微观损伤模型中的参数取:影响域半径ℓ=4 mm,临界伸长量=1.38×10−2mm,γ=2.0×105,以及p=4 和q=4.

4.1 不含圆孔情况

注意到,经典的局部损伤力学存在网格敏感性[18]:一旦所分析的固体存在或出现初始裂纹时,当计算点的间距加密时,基于经典局部损伤力学的分析结果,将会出现即使是很微小的外力作用就可以引起严重的损伤(当网格非常密集时,更甚者会出现不需要任何外力就可以使固体发生破坏的问题) 这一违背物理真实的现象.

为了考察本文损伤模型的网格敏感性问题,在有限元离散时采用3 种不同疏密程度的三角形网格进行划分.为此,首先考虑不含圆孔的情况,试件形状和尺寸如图4 所示.

图4 剪切梁基本尺寸(切口宽度为2 mm)Fig.4 Size of the shear beam(The notch width=2 mm)

具体地,考虑在切口一定范围的区域进行单元加密处理,3 种不同疏密程度的单元总数目分别为17 308,23 564 和55 262,单元分布情况见图5.

图5 有限元网格划分Fig.5 Finite element meshes

图6 裂纹发展模式Fig.6 The crack development modes

图6 裂纹发展模式(续)Fig.6 The crack development modes(continued)

图7 不同网格的载荷−位移曲线Fig.7 Load-deformation curves for different meshes

采用本文方法进行分析和求解.3 种不同网格划分情况下的裂纹扩展模式如图6 所示.与图6(d)的试验结果对比可见,3 种不同网格的分析结果均与试验结果相符,且随着网格加密裂纹扩展路径趋于试验结果.图7 进一步给出了外加载荷−切口张开位移和外加载荷−加载点位移曲线.特别值得指出的是,图7(b)表现出一般数值方法很难处理的强非线性回弹性质.从图中可见,3 种不同网格的数值分析结果具有良好的一致性.可见采用本文建议方法分析的结果克服了经典局部损伤力学所面临的网格敏感性问题.

图8 是载荷−切口张开位移和载荷−加载点位移曲线与试验结果的对比.从图8(a)可见,数值模拟的载荷−切口张开位移曲线与试验曲线的吻合良好.在图8(b)中,虽然未能完美吻合,但载荷−加载点位移曲线的最大载荷与最终加载点位移与实验值的接近,且总体的定性特征一致.

图8 数值与试验结果比较Fig.8 Comparisons between the numerical and experimental results

图9 和图10 分别给出了本文建议方法与基于混合元的各向同性损伤模型[51]和统一相场模型[23]分析所获得的载荷−变形曲线与裂纹模式的对比,可见三者较为一致.

图9 与其他数值模型的载荷−位移曲线比较Fig.9 The load-deformation curves compared with other numerical models

图10 与其他数值模型破坏模式的比较Fig.10 The comparisons of failure modes with other numerical models

图11 典型载荷步损伤与位移云图Fig.11 Damage and displacement field at typical loading steps

图11 典型载荷步损伤与位移云图(续)Fig.11 Damage and displacement field at typical loading steps(continued)

图11 是图8 的载荷−变形曲线中几个典型阶段的损伤云图和位移云图演化情况.从图中可见,裂纹随着加载过程而不断发展,同时位移场在裂纹两侧出现显著的不连续性.特别值得注意的是,当外加载荷已经达到峰值的时候,裂纹才刚刚具有肉眼可见的发展(图11(a)),从图11(b)亦可见,此时位移场尚未出现显著的不连续性.当裂纹明显萌生和发展时(图11(c),图11(e),图11(g)),试件的载荷−变形曲线已经进入显著的下降段(图8(a))或回弹段(图8(b)).换言之,对于此类问题,当出现肉眼可见的裂纹时,裂纹已经进入不稳定发展阶段.这与人们从实际工程的损伤和断裂分析中获得的经验完全相符.

4.2 含圆孔情况

考虑含圆孔情况,试件形状与尺寸如图12 所示.类似地采用3 种不同疏密程度的三角形网格进行划分,单元数目和网格划分具体情况见图13.

图12 剪切梁基本尺寸(切口宽度为2 mm)Fig.12 Size of the shear beam(the notch width=2 mm)

图13 有限元网格划分Fig.13 Finite element meshes

图14 裂纹发展模式Fig.14 The crack development modes

采用本文建议的方法,得到3 种不同的网格划分获得的裂纹最终发展模式如图14 所示.与图14(d)的试验结果对比表明,3 种网格均能良好地捕捉裂纹发展过程,并给出了试验未展示的裂纹后续扩展部分.图15 给出了载荷−变形曲线的分析结果.同样地,在图15(b) 中也表现出强非线性回弹特性.由此可见,采用不同的网格时,无论是载荷−裂纹张开位移曲线、还是载荷−加载点位移曲线都具有很好的一致性,不存在经典局部损伤力学中的网格敏感性问题.

图15 不同网格的载荷−位移曲线Fig.15 Load-deformation curves for different meshes

图16 是采用网格C 时的载荷−变形曲线.与图8 类似,此时载荷−加载点位移曲线出现了明显的回弹现象.

图17 是与图16 中标注的典型阶段相应的损伤云图与位移云图.与前述情况类似地可见,当外载荷达到峰值时,试件中才刚刚出现可见裂纹的扩展.而且,裂纹两侧位移场出现显著的不连续性.有意思的是,当存在圆孔时,裂纹的发展会被圆孔所诱导(图17(c)),并最终穿透圆孔继续发展(图17(e)).值得注意的是,在被圆孔俘虏后,当裂纹进一步穿透圆孔时,出现了载荷的小幅上升(见图16),然后继续下降,直至完全破坏并基本丧失承载力.