基于相场法的损伤模型主要用于研究裂纹的扩展演化问题。
裂纹的扩散拓扑
图1 尖锐裂纹与扩散拓扑裂纹:
图2 尖锐裂纹与扩散裂纹模型:
尖锐裂纹可以采用一个辅助的场变量 d 来描述,即
d(x)={1x=00x≠0
显然,d=0 对应未断裂状态,d=1 则对应完全断裂状态。将这个辅助场变量 d(x) 称为裂纹相场。裂纹相场与连续损伤力学相关,因为连续损伤力学中的标量损伤度 d 在宏观层面上描述了微裂纹和微孔洞的均匀化发展。为了改善式(1)中相场的非连续性,可以采用指数模型,即
d(x)=e−|x|/l
d(0)=1 , d(±∞)=0
式(2)即表示了裂纹的扩散拓扑(正则化),其函数图像如图2(b)所示。裂纹的扩散拓扑程度由长度尺度参数 l 控制,当 l→0 时,式(2)等价于式(1)。
能量耗散密度函数
考虑如下的齐次微分方程
d(x)−l2d″(x)=0
注意到式(2)是上式的一个解,而该微分方程是下面变分问题的欧拉方程
mind(x) I(d)=12∫V(d2+l2d′2)dv
这里引入泛函的目的有两个:1、定义能量耗散密度函数,2、推导裂纹相场的控制方程
令 dv=Γdx ,于是有
I(d=e−|x|/l)=Γ∫+∞−∞e−2|x|/ldx=lΓ∫+∞0e−2x/ld(−2x/l)=lΓ
由此引入能量耗散函数 Γl 与能量耗散密度函数 γ
Γl(d)=1lI(d)=∫Vγ dv
γ(d)=12l(d2+l2d′2)
能量耗散函数 Γl 是关于裂纹相场 d(x) 的泛函。将式(6)中的一维形式推广到多维,即
Γl(d,∇d)=∫Vγ dv
γ(d,∇d)=12l(d2+l2∇d⋅∇d)
此时裂纹相场 d(x) 是关于空间位置坐标 x 的函数。
裂纹相场的控制方程
求式(7)中能量耗散函数的变分问题,即
mind(x) Γl(d)=12l∫V(d2+l2∇d⋅∇d)dv
取能量耗散函数的变分,即
δΓl(d,∇d)=∫V(∂γ∂dδd+∂γ∂∇d⋅δ∇d)dv=1l∫V[dδd+l2∇d⏟u⋅δ∇d⏟v′]dv=1l[∫Vdδddv+l2∇d⋅nδd|∂V−∫Vl2∇2d=Δd⏞∇⋅(∇d)δddv⏟分部积分]=1l∫V(d−l2Δd)δddv+l∇d⋅nδd|∂V
分部积分公式:∫uv′dx=uv−∫u′vdx
导数的变分等于变分的导数:δ∇d=∇δd
显然变分 δΓl=0 的充分条件为
是否为必要条件?我不知道
{d−l2Δd=0 in V∇d⋅n=0 on ∂V
式中,n 为自然边界 ∂V 的外法线矢量。同时裂纹相场还需满足强制边界条件
d(x)=1on Γ
式中,Γ 表示一维的开裂点、二维的裂纹线、三维的裂纹面。
综上,裂纹相场的控制方程为
{d−l2Δd=0 in V∇d⋅n=0 on ∂Vd(x)=1 on Γ
实际上,式(12)给定了椭圆型的偏微分方程、Neumann边界条件以及Dirichlet边界条件,该控制方程描述了裂纹相场的定常状态。
MATLAB的pdetool
工具箱支持求解二维的椭圆型、抛物型、双曲型以及特征值模式的偏微分方程,求解方法为有限元法。在pdetool
工具箱中椭圆型偏微分方程定义为
au−cΔu=f
按照式(12),令 a=1/l2、c=1、f=0 即可。图3展示了求解域为 2×2 的正方形,裂纹长度为 1,长度尺度参数 l=0.25 的相场解。
图3 裂纹相场:
由pdetool
工具箱生成的代码如下:
1 | function pde_phase_field_damage_line |
裂纹扩展的变分原理
考虑由裂纹相场引起的能量耗散,构造总势能泛函,当泛函的变分为零时,泛函取得极值,此时连续体处于平衡状态,即最小势能原理。
连续体的总势能由应变能、动能、耗散势能、载荷势能组成,即
Π(u,˙u,ε,d,∇d)⏟总势能=E(ε,d)⏟应变能+D(d,∇d)⏟耗散势能+K(˙u)⏟动能−W(u)⏟载荷势能
于是总势能的变分为
δΠ=δE+δD+δK−δW=0
应变能的变分
对于没有裂纹(损伤)的各向同性线弹性材料,其应变能密度为
ψ(ε)=με:ε+12λtr2(ε)
其应力为
σ=∂ψ(ε)∂ε=2με+λtr(ε)I
式(17)的证明较为容易,将式(16)展开后再对 εij 求导即可
由于裂纹的存在,必然导致材料力学性能的退化,于是可以定义退化后的应变能密度
˜ψ(ε,d)=[g(d)+k]ψ+(ε)+ψ−(ε)
式中,g(d) 为退化函数,需满足 g(0)=1 、g(1)=0 ;k 为正的小量,其作用为防止 d=1 时能量的完全退化导致刚度矩阵奇异;ψ+(ε) 、ψ−(ε) 分别表示由拉伸、压缩主导的应变能密度,且需要满足下式
ψ(ε)=ψ+(ε)+ψ−(ε)
式(18a)表明仅考虑拉伸引起的损伤(裂纹开裂),这也反映了各向异性的性能退化。退化函数的选取,以及应变能密度的拉伸、压缩分解将在后续进行讨论。
将式(18a)在求解域 V 上积分,可得连续体的应变能
E(ε,d)=∫V˜ψ(ε,d)dv
于是应变能的变分为
δE(ε,d)=∫V(∂˜ψ∂ε:δε+∂˜ψ∂dδd)dv=∫V{[g(d)+k]∂εψ++∂εψ−}:δεdv+∫Vg′(d)ψ+δddv
耗散势的变分
耗散势通过对耗散势密度函数在求解域 V 上积分得到,即
D(d,∇d)=∫Vφ(d,∇d)dv
式中,耗散势密度函数的定义为
φ(d,∇d)=gcγ(d,∇d)
式中,gc 为材料参数(为了统一量纲)。于是由式(7)可得
D(d,∇d)=gcΓl(d,∇d)
于是根据式(9)可得
δD(d,∇d)=gcδΓl(d,∇d)=gcl∫V(d−l2Δd)δddv+gcl∇d⋅nδd|∂V
动能与载荷势能的变分
动能的定义
K(˙u)=∫V12ρ˙u⋅˙udv
载荷势能的定义
W(u)=∫Vb⋅udv+∮∂Vt⋅uds
载荷势能并不等于载荷所作的功
在《弹性体的应变能密度》中给出动能与载荷势能的变分,即
δK=∫Vρ¨u⋅δudv
δW=∫V(∇⋅σ+b)⋅δudv+∫Vσ:δεdv
耦合方程
于是由将式(20)、(23)、(25)代入式(15)可得总势能变分
δΠ= δE+δD+δK−δW= ∫V(ρ¨u−∇⋅σ−b)⋅δudv+∫V[gcl(d−l2Δd)+g′(d)ψ+]δddv+ ∫V{[g(d)+k]∂εψ++∂εψ−−σ}:δεdv+gcl∇d⋅nδd|∂V=0
显然变分 δΠ=0 的充分条件为
{∇⋅σ+b=ρ¨ugcl(d−l2Δd)=−g′(d)ψ+σ=[g(d)+k]∂εψ++∂εψ−
为了体现损伤的不可逆性,即 ˙d≥0,引入能量历史最大函数
H(ε)=maxt∈[0,tf]ψ+(ε;t)
耗散势变分 δD≥0,仅需 −g′(d)ψ+≥0 即可,而 ψ+ 肯定是不小于零
所以引入 H(ε),不是很懂
于是,将式(27)改写为
{∇⋅σ+b=ρ¨ugcl(d−l2Δd)=−g′(d)Hσ=[g(d)+k]∂εψ++∂εψ−
式(29a)是耦合损伤的平衡方程与裂纹相场方程。相应的自然边界条件为
{n⋅σ=t on ∂V∇d⋅n=0 on ∂V
式(25b)的推导过程中使用了应力边界条件 n⋅σ=t
应变能密度的退化
这里具体讨论式(18a)中的退化函数 g(d) 以及拉伸、压缩应变能密度 ψ+(ε) 、ψ−(ε) 。
退化函数
退化函数 g(d) 首先需满足
{g(0)=1g(1)=0
这两种情况分别对应无裂纹时应变能没有退化、开裂时拉伸应变能完全退化。
考虑到式(27)或式(29b)中的第二个式子,当 d=1 时,为了保证弹性驱动力 ∂d˜ψ=g′(d)ψ+ 收敛到一个有限值1,于是令
g′(1)=0
文献[1]里这样解释的,我是没看懂
于是在式(30)的条件下选取退化函数,例如
g(d)=(1−d)2
应变能密度的拉压分解
首先定义两个算符
{⟨x⟩+=(x+|x|)/2⟨x⟩−=(x−|x|)/2
为了方便展示,再将应变能密度函数给出
ψ(ε)=με:ε+12λtr2(ε)
分解形式1
文献[2]中将应变张量 ε 进行谱分解,即
ε=3∑i=1εini⊗ni
式中,εi、ni 分别为应变张量 ε 的第 i 个主应变、主方向,且 ‖ni‖=1 。特别地,令 ni 均为列向量,记 Q=[n1,n2,n3] ,则有
diag([ε1,ε2,ε3])=QT⋅ε⋅Q
实对称矩阵的谱分解性质
式中,diag([ε1,ε2,ε3]) 表示主对角线元素为 [ε1,ε2,ε3] 的对角阵。事实上 Q 为正交矩阵,式(33b)表示应变张量 ε 经正交矩阵 Q 旋转后得到主应变微体。
于是,在谱分解的基础上定义拉伸应变张量 ε+ 与压缩应变张量 ε−
{ε+=∑3i=1⟨εi⟩+ni⊗niε−=∑3i=1⟨εi⟩−ni⊗ni
显然,对于式(34)所定义拉伸、压缩应变张量,恒有 ε=ε++ε− 。
进一步可以定义拉伸、压缩应变能密度
{ψ+(ε)=με+:ε++12λ⟨tr(ε)⟩2+ψ−(ε)=με−:ε−+12λ⟨tr(ε)⟩2−
证明1
还需证明式(35)满足式(18b)。不难证明
tr2(ε)=⟨tr(ε)⟩2++⟨tr(ε)⟩2−
因此还需证明
ε:ε=ε+:ε++ε−:ε−
这里先证明两个子结论:令一阶张量 m、n 满足 ‖m‖=‖n‖=1 ,且 m⋅n=0 ,于是有
(n⊗n):(n⊗n)=ninjninj=∑i[n2i(∑jn2j)]=1
(n⊗n):(m⊗m)=ninjmimj=∑i[nimi(∑jnjmj)]=0
于是根据式(36),可得
ε:ε=(3∑i=1εini⊗ni):(3∑i=1εini⊗ni)=3∑i=1ε2i
ε+:ε+=(3∑i=1⟨εi⟩+ni⊗ni):(3∑i=1⟨εi⟩+ni⊗ni)=3∑i=1⟨εi⟩2+
ε−:ε−=(3∑i=1⟨εi⟩−ni⊗ni):(3∑i=1⟨εi⟩−ni⊗ni)=3∑i=1⟨εi⟩2−
进一步可得
ε+:ε++ε−:ε−=3∑i=1ε2i
由式(37a)、(38)即可证明 ε:ε=ε+:ε++ε−:ε− 。
综上,由式(35)所定义的拉伸、压缩应变能密度满足式(18b)。
需要注意,直接将式(34)所定义的拉伸、压缩应变 ε+、ε− 代入式(16),令 {ψ+=ψ(ε+)=με+:ε++12λtr2(ε+)ψ−=ψ(ε−)=με−:ε−+12λtr2(ε−)
这显然不能满足式(18b)。
偏导数1
由式(27)或式(29b)中的第三个式子可知,还需 ψ+、ψ− 对 ε 的偏导数 ∂εψ+ 、∂εψ− 。于是由式(35)可得
{∂εψ+=2με++λ⟨tr(ε)⟩+I∂εψ−=2με−+λ⟨tr(ε)⟩−I
这个不知道怎么证明
分解形式2
文献[1]中采用主应变形式的应变能密度,即
ψ(ε1,ε2,ε3)=μ(ε21+ε22+ε23)+12λ(ε1+ε2+ε3)2
相应的拉伸、压缩应变能密度定义为
ψ±=μ(⟨ε1⟩2±+⟨ε2⟩2±+⟨ε3⟩2±)+12λ⟨ε1+ε2+ε3⟩2±
事实上,分解形式2与分解形式1是一致的
证明2
由式(41)所定义的拉伸、压缩应变能密度,显然满足式(18b)。
偏导数2
∂εψ±=2μ3∑i=1⟨εi⟩±ni⊗ni+λ⟨ε1+ε2+ε3⟩±I
分解形式3
文献[3]中将应变张量 ε 分解为球形部分和偏斜部分,即
{εv=13tr(ε)Iεd=ε−εv
于是,在上述分解的基础上定义拉伸应变张量 ε+ 与压缩应变张量 ε−
{ε+=ε+v+εdε−=ε−v
式中,ε+v=13⟨tr(ε)⟩+I , ε−v=13⟨tr(ε)⟩−I 。直接将式(44)所定义的拉伸、压缩应变 ε+、ε− 代入式(16),令
{ψ+=ψ(ε+)=με+:ε++12λtr2(ε+)ψ−=ψ(ε−)=με−:ε−+12λtr2(ε−)
证明3
当 tr(ε)>0 时,ε+=ε,ε−=0 。于是 ψ+=ψ(ε),ψ−=0,显然满足式(18b)。
当 tr(ε)<0 时,ε+=εd,ε−=εv。于是 ψ+=ψ(εd),ψ−=ψ(εv) 。
在《弹性体的应变能密度》中给出畸变能密度 ψd 与体积改变能密度 ψv
{ψd=ψ(εd)=μ[ε:ε−13tr2(ε)]ψv=ψ(εv)=16(2μ+3λ)tr2(ε)
且 ψd+ψv=ψ,于是满足式(18b)。
综上,由式(45)所定义的拉伸、压缩应变能密度满足式(18b)。
偏导数3
ψd、ψv 对 ε 的偏导数 ∂εψd、∂εψv
{∂εψd=μ[2ε−23tr(ε)I]=2μεd∂εψv=13(2μ+3λ)tr(ε)I=(2μ+3λ)εv
于是有
{∂εψ+=13(2μ+3λ)⟨tr(ε)⟩+I+2μεd∂εψ−=13(2μ+3λ)⟨tr(ε)⟩−I
参考文献
- [1] Miehe C , Welschinger F , Hofacker M . Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[J]. International Journal for Numerical Methods in Engineering, 2010, 83(10):1273-1311.
- [2] Hofacker M , Miehe C . Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation[J]. International Journal of Fracture, 2012, 178(1-2):p.113-129.
- [3] Amor H , Marigo J J , Maurini C . Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[J]. Journal of the Mechanics & Physics of Solids, 2009, 57(8):1209-1229.