0%

相场损伤模型

基于相场法的损伤模型主要用于研究裂纹的扩展演化问题。

裂纹的扩散拓扑

图1 尖锐裂纹与扩散拓扑裂纹:

图2 尖锐裂纹与扩散裂纹模型:

尖锐裂纹可以采用一个辅助的场变量 d 来描述,即

d(x)={1x=00x0

显然,d=0 对应未断裂状态,d=1 则对应完全断裂状态。将这个辅助场变量 d(x) 称为裂纹相场。裂纹相场与连续损伤力学相关,因为连续损伤力学中的标量损伤度 d 在宏观层面上描述了微裂纹和微孔洞的均匀化发展。为了改善式(1)中相场的非连续性,可以采用指数模型,即

d(x)=e|x|/l

d(0)=1 , d(±)=0

式(2)即表示了裂纹的扩散拓扑(正则化),其函数图像如图2(b)所示。裂纹的扩散拓扑程度由长度尺度参数 l 控制,当 l0 时,式(2)等价于式(1)。

能量耗散密度函数

考虑如下的齐次微分方程

d(x)l2d(x)=0

注意到式(2)是上式的一个解,而该微分方程是下面变分问题的欧拉方程

mind(x) I(d)=12V(d2+l2d2)dv

这里引入泛函的目的有两个:1、定义能量耗散密度函数,2、推导裂纹相场的控制方程

dv=Γdx ,于是有

I(d=e|x|/l)=Γ+e2|x|/ldx=lΓ+0e2x/ld(2x/l)=lΓ

由此引入能量耗散函数 Γl 与能量耗散密度函数 γ

Γl(d)=1lI(d)=Vγ dv

γ(d)=12l(d2+l2d2)

能量耗散函数 Γl 是关于裂纹相场 d(x) 的泛函。将式(6)中的一维形式推广到多维,即

Γl(d,d)=Vγ dv

γ(d,d)=12l(d2+l2dd)

此时裂纹相场 d(x) 是关于空间位置坐标 x 的函数。

裂纹相场的控制方程

求式(7)中能量耗散函数的变分问题,即

mind(x) Γl(d)=12lV(d2+l2dd)dv

取能量耗散函数的变分,即

δΓl(d,d)=V(γdδd+γdδd)dv=1lV[dδd+l2duδdv]dv=1l[Vdδddv+l2dnδd|VVl22d=Δd(d)δddv分部积分]=1lV(dl2Δd)δddv+ldnδd|V

分部积分公式:uvdx=uvuvdx

导数的变分等于变分的导数:δd=δd

显然变分 δΓl=0 的充分条件为

是否为必要条件?我不知道

{dl2Δd=0 in Vdn=0 on V

式中,n 为自然边界 V 的外法线矢量。同时裂纹相场还需满足强制边界条件

d(x)=1on Γ

式中,Γ 表示一维的开裂点、二维的裂纹线、三维的裂纹面。

综上,裂纹相场的控制方程为

{dl2Δd=0 in Vdn=0 on Vd(x)=1 on Γ

实际上,式(12)给定了椭圆型的偏微分方程、Neumann边界条件以及Dirichlet边界条件,该控制方程描述了裂纹相场的定常状态

MATLAB的pdetool工具箱支持求解二维的椭圆型、抛物型、双曲型以及特征值模式的偏微分方程,求解方法为有限元法。在pdetool工具箱中椭圆型偏微分方程定义为

aucΔu=f

按照式(12),令 a=1/l2c=1f=0 即可。图3展示了求解域为 2×2 的正方形,裂纹长度为 1,长度尺度参数 l=0.25 的相场解。

图3 裂纹相场:

pdetool工具箱生成的代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
function pde_phase_field_damage_line
[fig,ax]=pdeinit;
pdetool('appl_cb',1);set(ax,'DataAspectRatio',[1 1 1]);
set(ax,'PlotBoxAspectRatio',[730.45 486 486]);
set(ax,'XLimMode','auto');set(ax,'YLim',[-1 1]);
set(ax,'XTickMode','auto');set(ax,'YTickMode','auto');

% Geometry description:
pderect([0 1 1 -1],'R1');pderect([-1 0 1 0.001],'R2');
pderect([-1 0 -0.0010000000000000009 -1],'R3');
set(findobj(get(fig,'Children'),'Tag','PDEEval'),'String','R1+R2+R3')

% Boundary conditions:
pdetool('changemode',0)
pdesetbd(12,'neu',1,'0','0');pdesetbd(11,'neu',1,'0','0');
pdesetbd(10,'neu',1,'0','0');pdesetbd(9,'neu',1,'0','0');
pdesetbd(7,'dir',1,'1','1');pdesetbd(5,'neu',1,'0','0');
pdesetbd(4,'neu',1,'0','0');pdesetbd(3,'dir',1,'1','1');
pdesetbd(2,'dir',1,'1','1');pdesetbd(1,'neu',1,'0','0');

% Mesh generation:
setappdata(fig,'Hgrad',1.3);
setappdata(fig,'refinemethod','regular');
setappdata(fig,'jiggle',char('on','mean',''));
setappdata(fig,'MesherVersion','preR2013a');
pdetool('initmesh');pdetool('refine');

% PDE coefficients:
pdeseteq(1,'1','16','0','1','0:10','0','0','[0 100]');
setappdata(fig,'currparam',['1.0';'10 ';'0 ';'1.0']);

% Solve parameters:
setappdata(fig,'solveparam',char('0','6570','10','pdeadworst',...
'0.5','longest','0','1E-4','','fixed','Inf'));

% Plotflags and user data strings:
setappdata(fig,'plotflags',[1 1 1 1 1 1 7 1 0 0 0 1 1 0 0 0 0 1]);
setappdata(fig,'colstring','');setappdata(fig,'arrowstring','');
setappdata(fig,'deformstring','');setappdata(fig,'heightstring','');

% Solve PDE:
pdetool('solve')

裂纹扩展的变分原理

考虑由裂纹相场引起的能量耗散,构造总势能泛函,当泛函的变分为零时,泛函取得极值,此时连续体处于平衡状态,即最小势能原理

连续体的总势能由应变能、动能、耗散势能、载荷势能组成,即

Π(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)=1g(1)=0k 为正的小量,其作用为防止 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)=gclV(dl2Δd)δddv+gcldnδd|V

动能与载荷势能的变分

动能的定义

K(˙u)=V12ρ˙u˙udv

载荷势能的定义

W(u)=Vbudv+Vtuds

载荷势能并不等于载荷所作的功

《弹性体的应变能密度》中给出动能与载荷势能的变分,即

δK=Vρ¨uδudv

δW=V(σ+b)δudv+Vσ:δεdv


耦合方程

于是由将式(20)、(23)、(25)代入式(15)可得总势能变分

δΠ= δE+δD+δKδW= V(ρ¨uσb)δudv+V[gcl(dl2Δd)+g(d)ψ+]δddv+  V{[g(d)+k]εψ++εψσ}:δεdv+gcldnδd|V=0

显然变分 δΠ=0 的充分条件为

{σ+b=ρ¨ugcl(dl2Δd)=g(d)ψ+σ=[g(d)+k]εψ++εψ

为了体现损伤的不可逆性,即 ˙d0,引入能量历史最大函数

H(ε)=maxt[0,tf]ψ+(ε;t)

耗散势变分 δD0,仅需 g(d)ψ+0 即可,而 ψ+ 肯定是不小于零

所以引入 H(ε),不是很懂

于是,将式(27)改写为

{σ+b=ρ¨ugcl(dl2Δd)=g(d)Hσ=[g(d)+k]εψ++εψ

式(29a)是耦合损伤的平衡方程与裂纹相场方程。相应的自然边界条件为

{nσ=t on Vdn=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)=(1d)2

应变能密度的拉压分解

首先定义两个算符

{x+=(x+|x|)/2x=(x|x|)/2

为了方便展示,再将应变能密度函数给出

ψ(ε)=με:ε+12λtr2(ε)


分解形式1

文献[2]中将应变张量 ε 进行谱分解,即

ε=3i=1εinini

式中,εini 分别为应变张量 ε 的第 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+niniε=3i=1εinini

显然,对于式(34)所定义拉伸、压缩应变张量,恒有 ε=ε++ε

进一步可以定义拉伸、压缩应变能密度

{ψ+(ε)=με+:ε++12λtr(ε)2+ψ(ε)=με:ε+12λtr(ε)2

证明1

还需证明式(35)满足式(18b)。不难证明

tr2(ε)=tr(ε)2++tr(ε)2

因此还需证明

ε:ε=ε+:ε++ε:ε

这里先证明两个子结论:令一阶张量 mn 满足 m=n=1 ,且 mn=0 ,于是有

(nn):(nn)=ninjninj=i[n2i(jn2j)]=1

(nn):(mm)=ninjmimj=i[nimi(jnjmj)]=0

于是根据式(36),可得

ε:ε=(3i=1εinini):(3i=1εinini)=3i=1ε2i

ε+:ε+=(3i=1εi+nini):(3i=1εi+nini)=3i=1εi2+

ε:ε=(3i=1εinini):(3i=1εinini)=3i=1εi2

进一步可得

ε+:ε++ε:ε=3i=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

相应的拉伸、压缩应变能密度定义为

ψ±=μ(ε12±+ε22±+ε32±)+12λε1+ε2+ε32±

事实上,分解形式2与分解形式1是一致的

证明2

由式(41)所定义的拉伸、压缩应变能密度,显然满足式(18b)。

偏导数2

εψ±=2μ3i=1εi±nini+λε1+ε2+ε3±I


分解形式3

文献[3]中将应变张量 ε 分解为球形部分和偏斜部分,即

{εv=13tr(ε)Iεd=εεv

于是,在上述分解的基础上定义拉伸应变张量 ε+ 与压缩应变张量 ε

{ε+=ε+v+εdε=εv

式中,ε+v=13tr(ε)+I  ,  εv=13tr(ε)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.