止于至善

0%

相场损伤模型

发表于 更新于 分类于

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

裂纹的扩散拓扑

图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
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)=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.