单晶塑性abaqus-umat 自学1
本文以翻译黄永刚老师的umat文章为主， 这个颜色的为补充材料和个人的理解。
摘要
包含单晶塑性的abaqus用户材料子程序。本文回顾有限元形式的单晶弹塑性和黏弹塑性，包含小变形理论，以及严格的有限变形和有限旋转理论。单晶的非弹性变形由于晶体滑移，假定遵从Schmid准则。提出多种多样的分切应力与剪切变形的自硬化和潜硬化关系，并且包含在子程序内。
1. 简介
有限元程序ABAQUS已经广泛应用在固体的变形与应力分析。除了自身大量的本构模型，ABAQUS还提供用户子程序接口。应力、应变以及求解依赖的状态变量以增量形式求出。当umat子程序被访问时，提供增量步开始的状态（应力以及求解依赖的状态变量）和应变和时间增量。umat子程序执行二段程序：更新应力和状态变量使其等于增量步结束时对应的值，提供材料的Jacobian矩阵， ∂ Δ σ / ∂ Δ ε \partial \Delta\sigma/\partial \Delta\varepsilon ∂Δσ/∂Δε，为Newton-Rhapson迭代提供相应的本构模型。
本文的主要目标是提供单晶塑性的本构的框架。运动学部分严格满足有限变形。塑性变形仅仅考虑位错滑移，不考虑扩散变形，孪晶和晶界滑移。Schmid应力，或滑移系上分切应力假定是滑移的驱动力。完整的子程序用于单晶和双晶的变形和断裂分析。
2. 单晶弹塑性本构模型回顾
2.1运动学
晶体力学的运动学理论的先驱是Taylor（1938），精准的力学理论是Hill（1966），Rice（1971）和Hill and Rice（1972）搭建的。下面是关于Asaro and Rice（1977）和Asaro（1983）简单的理论总结。
晶体材料的弹性变形和旋转是由晶格承担的。单晶的非弹性变形源于晶体滑移。材料流动是通过位错在晶格上的运动。总的变形梯度
F
\bm{F}
F:
F
=
F
∗
⋅
F
P
\bm{F=F^*\cdot F^P}
F=F∗⋅FP
其中
F
P
F^P
FP 是中间构型的塑性剪切，
F
∗
F^*
F∗ 为晶格拉伸和旋转。假设弹性特征不受滑移的影响，应力由
F
∗
F^*
F∗ 单独决定。
F
P
F^P
FP 的变化率与
α
\alpha
α 滑移系的滑移率
γ
˙
(
α
)
\dot{\gamma}^{(\alpha)}
γ˙(α) 有关。
F
P
˙
⋅
F
P
−
1
˙
=
∑
α
γ
˙
(
α
)
s
(
α
)
m
(
α
)
\dot{\bm{F}^P}\cdot\dot{\bm{F}^{P-1}}=\sum_{\alpha}\dot{\gamma}^{(\alpha)}{s}^{(\alpha)}{m}^{(\alpha)}
FP˙⋅FP−1˙=α∑γ˙(α)s(α)m(α)
对所有激活的滑移系求和，
s
(
α
)
{s}^{(\alpha)}
s(α) 和
m
(
α
)
{m}^{(\alpha)}
m(α) 分别是参考构型下滑移方向和滑移面的法向。
补充：单滑移变形梯度描述
速度梯度张量
L
\bm L
L
L
=
F
˙
⋅
F
−
1
=
F
˙
∗
⋅
F
P
⋅
F
P
−
1
⋅
F
∗
−
1
+
F
∗
⋅
F
˙
P
⋅
F
P
−
1
⋅
F
∗
−
1
=
F
˙
∗
⋅
F
∗
−
1
+
F
∗
⋅
F
˙
P
⋅
F
P
−
1
⋅
F
∗
−
1
\begin{aligned} \bm{L}&=\bm{\dot{F}}\cdot\bm{F^{-1}}=\bm{\dot F^*\cdot F^P}\cdot\bm{ F^{P-1}\cdot F^{*-1}}+\bm{ F^*\cdot \dot F^P}\cdot\bm{ F^{P-1}\cdot F^{*-1}}\\ &=\bm{\dot F^*}\cdot\bm{ F^{*-1}}+\bm{ F^*\cdot \dot F^P}\cdot\bm{ F^{P-1}\cdot F^{*-1}}\\ \end{aligned}
L=F˙⋅F−1=F˙∗⋅FP⋅FP−1⋅F∗−1+F∗⋅F˙P⋅FP−1⋅F∗−1=F˙∗⋅F∗−1+F∗⋅F˙P⋅FP−1⋅F∗−1
求
F
P
−
1
\bm{ F^{P-1}}
FP−1 由于：
[
I
−
γ
(
α
)
s
(
α
)
m
(
α
)
]
⋅
[
I
+
γ
(
α
)
s
(
α
)
m
(
α
)
]
=
I
−
γ
(
α
)
2
s
(
α
)
m
(
α
)
⋅
s
(
α
)
m
(
α
)
\begin{aligned} &[\bm{I}-{\gamma}^{(\alpha)}{s}^{(\alpha)}{m}^{(\alpha)}]\cdot[\bm{I}+{\gamma}^{(\alpha)}{s}^{(\alpha)}{m}^{(\alpha)}]\\\quad\\ =&\bm{I}-{\gamma}^{(\alpha)^2}{s}^{(\alpha)}{m}^{(\alpha)}\cdot{s}^{(\alpha)}{m}^{(\alpha)} \end{aligned}
=[I−γ(α)s(α)m(α)]⋅[I+γ(α)s(α)m(α)]I−γ(α)2s(α)m(α)⋅s(α)m(α)
滑移方向和滑移面法向垂直：
m
(
α
)
⋅
s
(
α
)
=
0
{m}^{(\alpha)}\cdot{s}^{(\alpha)}=0
m(α)⋅s(α)=0
因此：
F
P
−
1
=
I
−
γ
(
α
)
s
(
α
)
m
(
α
)
\bm{ F^{P-1}}=\bm{I}-{\gamma}^{(\alpha)}{s}^{(\alpha)}{m}^{(\alpha)}
FP−1=I−γ(α)s(α)m(α)
F
P
˙
⋅
F
P
−
1
=
[
γ
˙
(
α
)
s
(
α
)
m
(
α
)
]
⋅
[
I
−
γ
(
α
)
s
(
α
)
m
(
α
)
]
=
γ
˙
(
α
)
s
(
α
)
m
(
α
)
\bm{ \dot{F^P} \cdot F^{P-1}}=[\dot{\gamma}^{(\alpha)}{s}^{(\alpha)}{m}^{(\alpha)}]\cdot[\bm{I}-{\gamma}^{(\alpha)}{s}^{(\alpha)}{m}^{(\alpha)}]=\dot{\gamma}^{(\alpha)}{s}^{(\alpha)}{m}^{(\alpha)}
FP˙⋅FP−1=[γ˙(α)s(α)m(α)]⋅[I−γ(α)s(α)m(α)]=γ˙(α)s(α)m(α)
为了方便定义
s
∗
(
α
)
s^{*(\alpha)}
s∗(α) 为
α
\alpha
α 滑移系在变形后的构型中的滑移方向：
s
∗
(
α
)
=
F
∗
⋅
s
(
α
)
s^{*(\alpha)}=F^* \cdot s^{(\alpha)}
s∗(α)=F∗⋅s(α)
定义
m
∗
(
α
)
m^{*(\alpha)}
m∗(α) 为
α
\alpha
α 滑移系在变形后的构型中的滑移面的法向：
m
∗
(
α
)
=
m
(
α
)
⋅
F
∗
−
1
m^{*(\alpha)}= m^{(\alpha)}\cdot F^{*-1}
m∗(α)=m(α)⋅F∗−1
因此：
L
=
F
˙
∗
⋅
F
∗
−
1
+
γ
˙
(
α
)
F
∗
⋅
s
(
α
)
m
(
α
)
⋅
F
∗
−
1
=
F
˙
∗
⋅
F
∗
−
1
+
γ
˙
(
α
)
s
∗
(
α
)
m
∗
(
α
)
=
L
∗
+
L
P
\begin{aligned} \bm{L}&=\bm{\dot F^*}\cdot\bm{ F^{*-1}}+\dot{\gamma}^{(\alpha)}\bm{ F^*\cdot }{s}^{(\alpha)}{m}^{(\alpha)}\bm{ \cdot F^{*-1}}\\\quad\\ &=\bm{\dot F^*}\cdot\bm{ F^{*-1}}+\dot{\gamma}^{(\alpha)}\bm{s}^{*(\alpha)}\bm{m}^{*(\alpha)}\\\quad\\ &=\bm{L^*}+\bm{L^P} \end{aligned}
L=F˙∗⋅F∗−1+γ˙(α)F∗⋅s(α)m(α)⋅F∗−1=F˙∗⋅F∗−1+γ˙(α)s∗(α)m∗(α)=L∗+LP
当前状态的速度梯度矩阵：
L
=
F
˙
⋅
F
−
1
=
D
+
Ω
\bm{L}=\bm{\dot{F}\cdot F^{-1}}=\bm{D}+\bm{\Omega}
L=F˙⋅F−1=D+Ω
其中，
D
D
D 为对称的拉伸率，
Ω
\Omega
Ω 为反对称的自旋张量。可分解为晶格部分（上标*）和塑性部分（上标P）：
D
=
D
∗
+
D
P
Ω
=
Ω
∗
+
Ω
P
\bm{D=D^*+D^P}\quad \bm{\Omega=\Omega^*+\Omega^P}
D=D∗+DPΩ=Ω∗+ΩP
满足：
D
∗
+
Ω
∗
=
F
˙
∗
⋅
F
∗
−
1
D
P
+
Ω
P
=
∑
α
γ
˙
(
α
)
s
∗
(
α
)
m
∗
(
α
)
\bm{D^*+\Omega^*}=\bm{\dot F^*}\cdot\bm{ F^{*-1}}\quad \bm{D^P+\Omega^P}=\sum_{\alpha}\dot{\gamma}^{(\alpha)}\bm{s}^{*(\alpha)}\bm{m}^{*(\alpha)}
D∗+Ω∗=F˙∗⋅F∗−1DP+ΩP=α∑γ˙(α)s∗(α)m∗(α)
2.2 本构
依据 Hill and Rice（1972），存在弹性势能，
Φ
=
Φ
(
F
∗
)
\Phi=\Phi(F^*)
Φ=Φ(F∗) ，保证对称的晶格变形率张量
D
∗
\bm{D^*}
D∗, 和Cauchy应力的Jaumann客观应变率满足如下关系：
σ
∇
∗
+
σ
(
I
:
D
∗
)
=
L
:
D
∗
(2.2.1)
\overset{\nabla}{\bm{\sigma}}^*+\bm{\sigma}\left( \bm{I:D^*}\right)=\bm{L:D^*}\tag{2.2.1}
σ∇∗+σ(I:D∗)=L:D∗(2.2.1)
其中
I
I
I是二阶单位张量，
L
\bm{L}
L 是弹性张量具有如下对称性
L
i
j
k
l
=
L
j
i
k
l
=
L
i
j
l
k
=
L
k
l
i
j
L_{ijkl}=L_{jikl}=L_{ijlk}=L_{klij}
Lijkl=Ljikl=Lijlk=Lklij，Jaumann率
σ
∇
∗
\overset{\nabla}{\bm{\sigma}}^*
σ∇∗ 是晶格坐标轴上的共转应力率，Jaumann率
σ
∇
\overset{\nabla}{\bm{\sigma}}
σ∇ 是材料坐标轴上的共转应力率，二者之间的关系：
σ
∇
∗
=
σ
∇
+
(
Ω
−
Ω
∗
)
⋅
σ
−
σ
⋅
(
Ω
−
Ω
∗
)
(2.2.2)
\overset{\nabla}{\bm{\sigma}}^*=\overset{\nabla}{\bm{\sigma}}+(\Omega-\Omega^*)\cdot \sigma- \sigma\cdot(\Omega-\Omega^*)\tag{2.2.2}
σ∇∗=σ∇+(Ω−Ω∗)⋅σ−σ⋅(Ω−Ω∗)(2.2.2)
其中，
σ
∇
=
σ
˙
−
Ω
⋅
σ
+
σ
⋅
Ω
\overset{\nabla}{\bm{\sigma}}=\dot{\bm{\sigma}}-\Omega\cdot \sigma+\sigma\cdot\Omega
σ∇=σ˙−Ω⋅σ+σ⋅Ω
补充：
在Lagrange坐标下，变形率
D
∗
D^*
D∗ 和 Kirchhoff应力
τ
∗
=
J
σ
∗
\tau^*=J\sigma^*
τ∗=Jσ∗ 共轭
τ
∇
∗
=
L
:
D
∗
\overset{\nabla}{\tau}^*=\bm{L:D^*}
τ∇∗=L:D∗
其中：
τ
∇
∗
=
τ
∗
˙
−
Ω
⋅
τ
∗
+
τ
∗
⋅
Ω
=
J
˙
σ
∗
+
J
σ
∗
˙
−
J
Ω
⋅
σ
∗
+
J
σ
∗
⋅
Ω
=
J
σ
∗
(
I
:
D
)
+
J
σ
∗
˙
−
J
Ω
⋅
σ
∗
+
J
σ
∗
⋅
Ω
=
J
[
σ
∗
(
I
:
D
)
+
σ
∗
∇
]
\begin{aligned} \overset{\nabla}{\tau}^*&=\dot{\tau^*}-\Omega\cdot \tau^*+\tau^*\cdot\Omega\\ &=\dot{J}\sigma^*+J\dot{\sigma^*}-J\Omega\cdot \sigma^*+J\sigma^*\cdot\Omega\\ &=J\sigma^*({I:D})+J\dot{\sigma^*}-J\Omega\cdot \sigma^*+J\sigma^*\cdot\Omega\\ &=J\left[\sigma^*({I:D})+\overset{\nabla}{\sigma^*}\right] \end{aligned}
τ∇∗=τ∗˙−Ω⋅τ∗+τ∗⋅Ω=J˙σ∗+Jσ∗˙−JΩ⋅σ∗+Jσ∗⋅Ω=Jσ∗(I:D)+Jσ∗˙−JΩ⋅σ∗+Jσ∗⋅Ω=J[σ∗(I:D)+σ∗∇]
考虑到从初始构型到中间构型只发生滑移，忽略体积变化
J
=
det
F
≈
1
J=\det{F}\approx 1
J=detF≈1
因此：
σ
∇
∗
+
σ
(
I
:
D
∗
)
=
L
:
D
∗
\overset{\nabla}{\bm{\sigma}}^*+\bm{\sigma}\left( \bm{I:D^*}\right)=\bm{L:D^*}
σ∇∗+σ(I:D∗)=L:D∗
晶体滑移遵循Schmid准则，任意滑移系
α
\alpha
α的滑移率
γ
˙
(
α
)
\dot{\gamma}^{(\alpha)}
γ˙(α)假设只与当前应力
σ
\sigma
σ 的Schmid应力
τ
(
α
)
{\tau}^{(\alpha)}
τ(α)。Schmid应力是忽略晶格畸变后的分切应力。Asaro and Rice（1977）讨论了现有有限弹性畸变的一些可能的一般化。本文基于Rice（1971），热力学应力与滑移共轭。定义：
τ
(
α
)
=
m
∗
(
α
)
⋅
ρ
0
ρ
σ
⋅
s
∗
(
α
)
(2.2.3)
\tau^{(\alpha)}=m^{*(\alpha)}\cdot\frac{\rho_0}{\rho}\sigma\cdot s^{*(\alpha)}\tag{2.2.3}
τ(α)=m∗(α)⋅ρρ0σ⋅s∗(α)(2.2.3)
其中
ρ
0
{\rho_0}
ρ0 和
ρ
{\rho}
ρ 参考构型和当前构型的密度；Hill and Rice(1972)认为
τ
(
α
)
\tau^{(\alpha)}\,
τ(α)等于随晶格旋转的Kirchhoff应力的最大剪切分量。Schmid应力的变化率如下：
τ
˙
(
α
)
=
m
∗
(
α
)
⋅
[
σ
∇
∗
+
σ
(
I
:
D
∗
)
−
D
∗
⋅
σ
+
σ
⋅
D
∗
]
⋅
s
∗
(
α
)
(2.2.4)
\dot\tau^{(\alpha)}=m^{*(\alpha)}\cdot \left[ \overset{\nabla}{\bm{\sigma}}^*+\bm{\sigma}\left( \bm{I:D^*}\right)-\bm{D^*\cdot\sigma+\sigma\cdot D^*}\right]\cdot s^{*(\alpha)}\tag{2.2.4}
τ˙(α)=m∗(α)⋅[σ∇∗+σ(I:D∗)−D∗⋅σ+σ⋅D∗]⋅s∗(α)(2.2.4)
补充：
式（2.2.3）中：
ρ
0
ρ
=
d
v
d
V
=
J
\frac{\rho_0}{\rho}=\frac{dv}{dV}=J
ρρ0=dVdv=J
因此：
τ
˙
(
α
)
=
m
˙
∗
(
α
)
⋅
J
σ
⋅
s
∗
(
α
)
+
m
∗
(
α
)
⋅
(
J
σ
)
˙
⋅
s
∗
(
α
)
+
m
∗
(
α
)
⋅
J
σ
⋅
s
˙
∗
(
α
)
\dot\tau^{(\alpha)}=\dot m^{*(\alpha)}\cdot J\sigma\cdot s^{*(\alpha)}+m^{*(\alpha)}\cdot \dot{(J\sigma)}\cdot s^{*(\alpha)}+m^{*(\alpha)}\cdot J\sigma\cdot \dot s^{*(\alpha)}
τ˙(α)=m˙∗(α)⋅Jσ⋅s∗(α)+m∗(α)⋅(Jσ)˙⋅s∗(α)+m∗(α)⋅Jσ⋅s˙∗(α)
其中：
m
∗
(
α
)
⋅
(
J
σ
)
˙
⋅
s
∗
(
α
)
=
m
∗
(
α
)
⋅
[
(
J
σ
)
∇
+
Ω
∗
⋅
σ
−
σ
⋅
Ω
∗
]
⋅
s
∗
(
α
)
=
m
∗
(
α
)
⋅
[
σ
∇
∗
+
σ
(
I
:
D
∗
)
−
Ω
∗
⋅
σ
+
σ
⋅
Ω
∗
]
⋅
s
∗
(
α
)
\begin{aligned} m^{*(\alpha)}\cdot \dot{(J\sigma)}\cdot s^{*(\alpha)}&=m^{*(\alpha)}\cdot[ \overset{\nabla}{(J\sigma)}+\Omega^*\cdot \sigma- \sigma\cdot\Omega^*]\cdot s^{*(\alpha)}\\&=m^{*(\alpha)}\cdot[\overset{\nabla}{\bm{\sigma}}^*+\bm{\sigma}\left( \bm{I:D^*}\right)-\Omega^*\cdot \sigma+ \sigma\cdot\Omega^*]\cdot s^{*(\alpha)} \end{aligned}
m∗(α)⋅(Jσ)˙⋅s∗(α)=m∗(α)⋅[(Jσ)∇+Ω∗⋅σ−σ⋅Ω∗]⋅s∗(α)=m∗(α)⋅[σ∇∗+σ(I:D∗)−Ω∗⋅σ+σ⋅Ω∗]⋅s∗(α)
忽略体积变化
J
=
det
F
≈
1
J=\det{F}\approx 1
J=detF≈1
m
˙
∗
(
α
)
⋅
σ
⋅
s
∗
(
α
)
=
m
(
α
)
⋅
F
˙
∗
−
1
⋅
σ
⋅
s
∗
(
α
)
=
m
∗
(
α
)
⋅
F
∗
⋅
F
˙
∗
−
1
⋅
σ
⋅
s
∗
(
α
)
=
m
∗
(
α
)
⋅
−
L
∗
⋅
σ
⋅
s
∗
(
α
)
\begin{aligned} \dot m^{*(\alpha)}\cdot \sigma\cdot s^{*(\alpha)}&=m^{(\alpha)}\cdot \dot{F}^{*-1}\cdot \sigma\cdot s^{*(\alpha)}\\&=m^{*(\alpha)}\cdot{F^*}\cdot \dot{F}^{*-1}\cdot \sigma\cdot s^{*(\alpha)} \\&=m^{*(\alpha)}\cdot -L^* \cdot \sigma\cdot s^{*(\alpha)} \end{aligned}
m˙∗(α)⋅σ⋅s∗(α)=m(α)⋅F˙∗−1⋅σ⋅s∗(α)=m∗(α)⋅F∗⋅F˙∗−1⋅σ⋅s∗(α)=m∗(α)⋅−L∗⋅σ⋅s∗(α)
注：
I
˙
=
(
F
∗
⋅
F
∗
−
1
)
˙
=
F
∗
˙
⋅
F
∗
−
1
+
F
∗
⋅
F
˙
∗
−
1
=
0
\dot\bm{I}=\dot\bm{(F^*\cdot F^{*-1})}=\dot\bm{F^*}\cdot \bm{F^{*-1}}+\bm{F^*}\cdot \bm{\dot F^{*-1}}=\bm 0
I˙=(F∗⋅F∗−1)˙=F∗˙⋅F∗−1+F∗⋅F˙∗−1=0
其中：
F
∗
˙
⋅
F
∗
−
1
=
L
∗
\dot\bm{F^*}\cdot \bm{F^{*-1}}=\bm{L^*}
F∗˙⋅F∗−1=L∗ ,因此
F
∗
⋅
F
˙
∗
−
1
=
−
L
∗
\bm{F^*}\cdot \bm{\dot F^{*-1}}=\bm{-L^*}
F∗⋅F˙∗−1=−L∗
同理：
m
∗
(
α
)
⋅
σ
⋅
s
˙
∗
(
α
)
=
m
(
∗
α
)
⋅
σ
⋅
F
∗
˙
⋅
s
(
α
)
=
m
∗
(
α
)
⋅
σ
⋅
F
∗
˙
⋅
F
∗
−
1
s
∗
(
α
)
=
m
∗
(
α
)
⋅
σ
⋅
L
∗
⋅
s
∗
(
α
)
\begin{aligned} m^{*(\alpha)}\cdot \sigma\cdot \dot s^{*(\alpha)}&=m^{(*\alpha)}\cdot \sigma\cdot\dot{F^*}\cdot s^{(\alpha)}\\&=m^{*(\alpha)}\cdot \sigma\cdot\dot{F^*}\cdot F^{*-1} s^{*(\alpha)} \\&=m^{*(\alpha)}\cdot \sigma\cdot L ^*\cdot s^{*(\alpha)} \end{aligned}
m∗(α)⋅σ⋅s˙∗(α)=m(∗α)⋅σ⋅F∗˙⋅s(α)=m∗(α)⋅σ⋅F∗˙⋅F∗−1s∗(α)=m∗(α)⋅σ⋅L∗⋅s∗(α)
因此：
τ
˙
(
α
)
=
m
˙
∗
(
α
)
⋅
[
σ
∇
∗
+
σ
(
I
:
D
∗
)
+
Ω
∗
⋅
σ
−
σ
⋅
Ω
∗
−
L
∗
⋅
σ
+
σ
⋅
L
∗
]
⋅
s
∗
(
α
)
=
m
˙
∗
(
α
)
⋅
[
σ
∇
∗
+
σ
(
I
:
D
∗
)
−
D
∗
⋅
σ
+
σ
⋅
D
∗
]
⋅
s
∗
(
α
)
\begin{aligned} \dot\tau^{(\alpha)}&=\dot m^{*(\alpha)}\cdot [\overset{\nabla}{\bm{\sigma}}^*+\bm{\sigma}\left( \bm{I:D^*}\right)+\Omega^*\cdot \sigma- \sigma\cdot\Omega^*-L ^{*}\cdot \sigma+\sigma\cdot L ^*]\cdot s^{*(\alpha)}\\ &=\dot m^{*(\alpha)}\cdot [\overset{\nabla}{\bm{\sigma}}^*+\bm{\sigma}\left( \bm{I:D^*}\right)-D^{*}\cdot \sigma+\sigma\cdot D^*]\cdot s^{*(\alpha)} \end{aligned}
τ˙(α)=m˙∗(α)⋅[σ∇∗+σ(I:D∗)+Ω∗⋅σ−σ⋅Ω∗−L∗⋅σ+σ⋅L∗]⋅s∗(α)=m˙∗(α)⋅[σ∇∗+σ(I:D∗)−D∗⋅σ+σ⋅D∗]⋅s∗(α)
2.3 率相关的晶体材料硬化
率无关的塑性是率相关塑性的极限。基于Schmid准则，
α
\alpha
α 滑移系的滑移率
γ
˙
(
α
)
\dot\gamma^{(\alpha)}
γ˙(α) 依赖于相应的分切应力
γ
˙
(
α
)
=
a
˙
(
α
)
f
(
α
)
(
τ
(
α
)
g
(
α
)
)
(2.3.1)
\dot\gamma^{(\alpha)}=\dot{a}^{(\alpha)}f^{(\alpha)}{\left(\frac{\tau^{(\alpha)}}{g^{(\alpha)}}\right)} \tag{2.3.1}
γ˙(α)=a˙(α)f(α)(g(α)τ(α))(2.3.1)
其中，
a
˙
(
α
)
\dot{a}^{(\alpha)}
a˙(α) 为参考应变率，
g
(
α
)
g^{(\alpha)}
g(α) 是滑移系当前的强度，
f
(
α
)
f^{(\alpha)}
f(α) 无量纲的描述应力和应变率的函数。Hutchinson（1976）利用简单幂律形式表示多晶蠕变：
f ( α ) ( x ) = x ∣ x ∣ n − 1 (2.3.1a) f^{(\alpha)}(x)=x|x|^{n-1}\tag{2.3.1a} f(α)(x)=x∣x∣n−1(2.3.1a)
其中，
n
n
n 为率敏感指数，
n
→
∞
n\to\infty
n→∞ 近似应变率无关。
应变硬化以增量形式给出：
g
(
α
)
=
∑
β
h
α
β
γ
˙
(
β
)
(2.3.2)
g^{(\alpha)}=\sum_{\beta}h_{\alpha\beta}\dot{\gamma}^{(\beta)}\tag{2.3.2}
g(α)=β∑hαβγ˙(β)(2.3.2)
其中
h
α
β
h_{\alpha\beta}
hαβ 为滑移的硬化模量，在所有激活的滑移系上求和。
h
α
α
h_{\alpha\alpha}
hαα称为自硬化系数，
h
α
β
h_{\alpha\beta}
hαβ
(
α
≠
β
)
({\alpha \ne \beta})
(α=β) 为潜硬化系数。
Peirce，Asaro and Needleman（1982）和Asaro（1983a，b）利用简单形式自硬化系数：
h
α
α
=
h
(
γ
)
=
h
0
sech
2
∣
h
0
γ
τ
s
−
τ
0
∣
(
α
不
求
和
)
(2.3.3a)
h_{\alpha\alpha}=h(\gamma)=h_0\operatorname{sech}^2\large{\left| \frac{h_0\gamma}{\tau_s-\tau_0} \right|}\quad (\alpha 不求和)\tag{2.3.3a}
hαα=h(γ)=h0sech2∣∣∣∣∣τs−τ0h0γ∣∣∣∣∣(α不求和)(2.3.3a)
其中
h
0
h_0
h0 为初始硬化模量，
τ
0
\tau_0
τ0屈服强度，等于当前强度的初始值
g
(
α
)
(
0
)
g^{(\alpha)}(0)
g(α)(0)，
γ
\gamma
γ 是所有滑移系上的Taylor累计剪切应变：
γ
=
∑
α
∫
0
t
∣
γ
˙
(
α
)
∣
d
t
(2.3.3b)
\gamma=\sum_\alpha\int_{0}^{t} \left|\dot{\gamma}^{(\alpha)}\right|dt\tag{2.3.3b}
γ=α∑∫0t∣∣∣γ˙(α)∣∣∣dt(2.3.3b)
潜硬化模量：
h
α
β
=
q
h
(
γ
)
(
α
≠
β
)
(2.3.3c)
h_{\alpha\beta}=qh(\gamma)\quad(\alpha \ne \beta)\tag{2.3.3c}
hαβ=qh(γ)(α=β)(2.3.3c)
q
q
q 为常数。这些表达式忽略Bauschinger效应。
Bassani and Wu(1991)三阶段硬化：
h
α
α
=
{
(
h
0
−
h
s
)
sech
2
∣
(
h
0
−
h
s
)
γ
τ
s
−
τ
0
∣
+
h
s
}
G
(
γ
(
β
)
;
β
≠
α
)
(
α
不
求
和
)
(2.3.4a)
h_{\alpha\alpha}=\left\{(h_0-h_s)\operatorname{sech}^2\large{\left| \frac{(h_0-h_s)\gamma}{\tau_s-\tau_0} \right|}+h_s\right\}G(\gamma^{(\beta)};\beta\ne\alpha)\quad (\alpha 不求和)\tag{2.3.4a}
hαα={(h0−hs)sech2∣∣∣∣∣τs−τ0(h0−hs)γ∣∣∣∣∣+hs}G(γ(β);β=α)(α不求和)(2.3.4a)
h
β
α
=
q
h
α
α
β
≠
α
(2.3.4b)
h_{\beta\alpha}=qh_{\alpha\alpha}\quad\beta\ne\alpha\tag{2.3.4b}
hβα=qhααβ=α(2.3.4b)
h
s
h_s
hs 为易滑移阶段硬化模量，
G
G
G 是交滑移相关的函数
G
(
γ
(
β
)
;
β
≠
α
)
=
1
+
∑
β
≠
α
f
α
β
tanh
(
γ
(
β
)
γ
0
)
(2.3.4c)
G(\gamma^{(\beta)};\beta\ne\alpha)=1+\sum_{\beta\ne\alpha}f_{\alpha\beta}\tanh\left(\frac{\gamma^{(\beta)}}{\gamma_0}\right)\tag{2.3.4c}
G(γ(β);β=α)=1+β=α∑fαβtanh(γ0γ(β))(2.3.4c)
γ
0
\gamma_0
γ0 是相互作用滑移系的滑移总量，
f
α
β
f_{\alpha\beta}
fαβ 相互作用的强度。例子，共面滑移的相互作用小于非共面滑移。
上述表达中没有明确的屈服，只要分切应力非零就会发生塑性剪切，然而，对于较大率敏感指数n（n≥50）,当分切应力小于 τ 0 \tau_0 τ0 塑性剪切率远远小于参考应变率 a ˙ \dot a a˙。 所以激活的滑移系没必要区分 ( s ∗ ( α ) , m ∗ ( α ) ) (s^{*(\alpha)},m^{*(\alpha)}) (s∗(α),m∗(α)) 和 ( − s ∗ ( α ) , m ∗ ( α ) ) (-s^{*(\alpha)},m^{*(\alpha)}) (−s∗(α),m∗(α))，所以当 τ ∗ ( α ) \tau^{*(\alpha)} τ∗(α)为负允许 γ ˙ ∗ ( α ) \dot\gamma^{*(\alpha)} γ˙∗(α)取负值。
3. 向前梯度时间积分方案和增量方程
本文两种时间积分方案。第一种方案在应力，应变和状态变量如剪切应变、分切应力、滑移当前的强度增量之间存在线性关系，这种假设在3.1-3.3中。应力和状态变量在时间增量步开始前进行更新。第二种方案是用Newton-Rhapson迭代求解非线性增量方程，详见3.4节。隐式时间积分方案，应力和状态变量在时间增量结束时进行计算。
3.1 向前梯度时间积分方案
率相关固体的切线模量方法被Peice,Shih and Needleman(1984)用在子程序中。我们定义在时间增量
Δ
t
\Delta t
Δt 内
α
\alpha
α 滑移系的应变增量为
γ
(
α
)
\gamma^{(\alpha)}
γ(α):
Δ
γ
(
α
)
=
γ
(
α
)
(
t
+
Δ
t
)
−
γ
(
α
)
(
t
)
(3.1.1)
\Delta\gamma^{(\alpha)}=\gamma^{(\alpha)}(t+\Delta t)-\gamma^{(\alpha)}(t)\tag{3.1.1}
Δγ(α)=γ(α)(t+Δt)−γ(α)(t)(3.1.1)
采用线性差值：
Δ
γ
(
α
)
=
Δ
t
[
(
1
−
θ
)
γ
˙
t
(
α
)
+
θ
γ
˙
t
+
Δ
t
(
α
)
]
(3.1.2)
\Delta\gamma^{(\alpha)}=\Delta t\left[(1-\theta)\dot\gamma^{(\alpha)}_{t}+\theta\dot\gamma^{(\alpha)}_{t+\Delta t}\right]\tag{3.1.2}
Δγ(α)=Δt[(1−θ)γ˙t(α)+θγ˙t+Δt(α)](3.1.2)
参数
θ
\theta
θ 取值范围0~1，取0为简单的Euler积分。
θ
\theta
θ 取值建议0.5~1（Peirce et al，1984）。
滑移率
γ
˙
(
α
)
\dot\gamma^{(\alpha)}
γ˙(α) 是分切应力
τ
(
α
)
\tau^{(\alpha)}
τ(α) 和当前滑移系强度
g
(
α
)
g^{(\alpha)}
g(α)的函数。滑移率的Taylor展开:
γ
˙
t
+
Δ
t
(
α
)
=
γ
˙
t
(
α
)
+
∂
γ
˙
t
(
α
)
∂
τ
(
α
)
Δ
τ
(
α
)
+
+
∂
γ
˙
t
(
α
)
∂
g
(
α
)
Δ
g
(
α
)
(3.1.3)
\dot\gamma^{(\alpha)}_{t+\Delta t}=\dot\gamma^{(\alpha)}_{t}+\frac{\partial \dot\gamma^{(\alpha)}_{t}}{\partial \tau^{(\alpha)}}\Delta\tau^{(\alpha)}++\frac{\partial \dot\gamma^{(\alpha)}_{t}}{\partial g^{(\alpha)}}\Delta g^{(\alpha)}\tag{3.1.3}
γ˙t+Δt(α)=γ˙t(α)+∂τ(α)∂γ˙t(α)Δτ(α)++∂g(α)∂γ˙t(α)Δg(α)(3.1.3)
联立(3.1.1)~(3.1.3)给出以下增量表达式：
Δ
γ
(
α
)
=
Δ
t
[
γ
˙
t
(
α
)
+
θ
∂
γ
˙
t
(
α
)
∂
τ
(
α
)
Δ
τ
(
α
)
+
θ
∂
γ
˙
t
(
α
)
∂
g
(
α
)
Δ
g
(
α
)
]
(3.1.4)
\Delta\gamma^{(\alpha)}=\Delta t \left[\dot\gamma^{(\alpha)}_{t}+\theta\,\frac{\partial \dot\gamma^{(\alpha)}_{t}}{\partial \tau^{(\alpha)}}\Delta\tau^{(\alpha)}+\theta\,\frac{\partial \dot\gamma^{(\alpha)}_{t}}{\partial g^{(\alpha)}}\Delta g^{(\alpha)}\right]\tag{3.1.4}
Δγ(α)=Δt[γ˙t(α)+θ∂τ(α)∂γ˙t(α)Δτ(α)+θ∂g(α)∂γ˙t(α)Δg(α)](3.1.4)
3.2 增量表达式
在这一部分中，已知时间增量为 Δ t \Delta t\, Δt应变增量为 Δ ε i j \Delta\varepsilon_{ij}\, Δεij，推导出所有滑移系上滑移应变增量 Δ γ ( α ) \Delta \gamma^{(\alpha)} Δγ(α)，分切应力增量 Δ τ ( α ) \Delta \tau^{(\alpha)}\, Δτ(α)和当前强度 Δ g ( α ) \Delta g^{(\alpha)}\, Δg(α)的关系。共转应力增量 Δ σ i j = σ ∇ i j Δ t \Delta \sigma_{ij}=\overset{\nabla}{\sigma}_{ij}\,\Delta t\, Δσij=σ∇ijΔt通过应变增量 Δ ε i j \Delta\varepsilon_{ij}\, Δεij表示。应力增量的更新和有限变形的ABAQUS有限元程序相同。
为了方便引入每个滑移系的“Schmid因子”
μ
i
j
(
α
)
\mu_{ij}^{(\alpha)}\,
μij(α)和张量
ω
i
j
(
α
)
\omega_{ij}^{(\alpha)}\,
ωij(α):
μ
i
j
(
α
)
=
1
2
[
s
i
∗
(
α
)
m
j
∗
(
α
)
+
s
j
∗
(
α
)
m
i
∗
(
α
)
]
(3.2.1a)
\mu_{ij}^{(\alpha)}=\frac{1}{2}\left[s_i^{*(\alpha)}m_j^{*(\alpha)}+s_j^{*(\alpha)}m_i^{*(\alpha)}\right] \tag{3.2.1a}
μij(α)=21[si∗(α)mj∗(α)+sj∗(α)mi∗(α)](3.2.1a)
ω
i
j
(
α
)
=
1
2
[
s
i
∗
(
α
)
m
j
∗
(
α
)
−
s
j
∗
(
α
)
m
i
∗
(
α
)
]
(3.2.1b)
\omega_{ij}^{(\alpha)}=\frac{1}{2}\left[s_i^{*(\alpha)}m_j^{*(\alpha)}-s_j^{*(\alpha)}m_i^{*(\alpha)}\right] \tag{3.2.1b}
ωij(α)=21[si∗(α)mj∗(α)−sj∗(α)mi∗(α)](3.2.1b)
张量
ω
i
j
(
α
)
\omega_{ij}^{(\alpha)}\,
ωij(α)和旋转张量
Ω
\Omega\,
Ω和
Ω
∗
\Omega^*\,
Ω∗的关系：
Ω
i
j
−
Ω
i
j
∗
=
∑
α
ω
i
j
(
α
)
γ
˙
(
α
)
(3.2.1c)
\Omega_{ij}-\Omega^*_{ij}=\sum_\alpha\omega_{ij}^{(\alpha)}\dot\gamma^{(\alpha)}\tag{3.2.1c}
Ωij−Ωij∗=α∑ωij(α)γ˙(α)(3.2.1c)
通过晶体滑移硬化方程(2.3.2)，得出当前硬化函数增量的表达式如下：
Δ
g
(
α
)
=
∑
β
h
α
β
Δ
γ
(
β
)
(3.2.2)
\Delta g^{(\alpha)}=\sum_\beta h_{\alpha\beta}\Delta\gamma^{(\beta)}\tag{3.2.2}
Δg(α)=β∑hαβΔγ(β)(3.2.2)
联立方程（2.2.4），弹性本构（2.2.1）以及将应变增量分解为晶格变形（2.1.4）和塑性变形（2.1.5）得出应变增量
Δ
ε
i
j
\Delta \varepsilon_{ij}\,
Δεij与分切应力
Δ
τ
(
α
)
\Delta\tau^{(\alpha)}
Δτ(α)的关系：
Δ
τ
(
α
)
=
[
L
i
j
k
l
μ
k
l
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
]
:
[
Δ
ε
i
j
−
∑
β
μ
i
j
(
β
)
Δ
γ
(
β
)
]
(3.2.3)
\Delta\tau^{(\alpha)}=\left[L_{ijkl}\mu^{(\alpha)}_{kl}+\omega^{(\alpha)}_{ik}\sigma_{jk}+\omega^{(\alpha)}_{jk}\sigma_{ik}\right]:\left[\Delta \varepsilon_{ij}-\sum_\beta\mu_{ij}^{(\beta)}\Delta\gamma^{(\beta)}\right]\tag{3.2.3}
Δτ(α)=[Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik]:⎣⎡Δεij−β∑μij(β)Δγ(β)⎦⎤(3.2.3)
其中
L
i
j
k
l
L_{ijkl}\,
Lijkl是弹性模量矩阵。通过方程（2.2.2）给出共转应力增量
Δ
σ
i
j
\Delta\sigma_{ij}\,
Δσij：
Δ
σ
i
j
=
L
i
j
k
l
Δ
ε
k
l
−
σ
i
j
Δ
ε
k
k
−
∑
α
[
L
i
j
k
l
μ
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
]
⋅
Δ
γ
(
α
)
(3.2.4)
\Delta\sigma_{ij}=L_{ijkl}\Delta\varepsilon_{kl}-\sigma_{ij}\Delta\varepsilon_{kk}-\sum_\alpha \left[L_{ijkl}\mu^{(\alpha)}+\omega^{(\alpha)}_{ik}\sigma_{jk}+\omega^{(\alpha)}_{jk}\sigma_{ik}\right]\cdot\Delta\gamma^{(\alpha)}\tag{3.2.4}
Δσij=LijklΔεkl−σijΔεkk−α∑[Lijklμ(α)+ωik(α)σjk+ωjk(α)σik]⋅Δγ(α)(3.2.4)
补充：(2.2.4)可改写为
τ
˙
(
α
)
=
m
∗
(
α
)
s
∗
(
α
)
:
[
σ
∇
∗
+
σ
(
I
:
D
∗
)
−
D
∗
⋅
σ
+
σ
⋅
D
∗
]
\dot\tau^{(\alpha)}=\bm{{m^{*(\alpha)} s^{*(\alpha)}}}:\left[ \overset{\nabla}{\bm{\sigma}}^*+\bm{\sigma}\left( \bm{I:D^*}\right)-\bm{D^*\cdot\sigma+\sigma\cdot D^*}\right]
τ˙(α)=m∗(α)s∗(α):[σ∇∗+σ(I:D∗)−D∗⋅σ+σ⋅D∗]
m
∗
(
α
)
s
∗
(
α
)
\bm{{m^{*(\alpha)} s^{*(\alpha)}}}\,
m∗(α)s∗(α)为向量
m
∗
(
α
)
\bm{{m^{*(\alpha)}}}\,
m∗(α)和向量
s
∗
(
α
)
\bm{{s^{*(\alpha)}}}\,
s∗(α)的并积，写成分量形式为
m
i
∗
(
α
)
s
j
∗
(
α
)
m_i^{*(\alpha)}s_j^{*(\alpha)}
mi∗(α)sj∗(α)。
由（3.2.1）可知：
m
i
∗
(
α
)
s
j
∗
(
α
)
=
μ
i
j
(
α
)
−
ω
i
j
(
α
)
m_i^{*(\alpha)}s_j^{*(\alpha)}=\mu_{ij}^{(\alpha)}-\omega_{ij}^{(\alpha)}
mi∗(α)sj∗(α)=μij(α)−ωij(α)
代入上式，可得：
τ
˙
(
α
)
=
m
i
∗
(
α
)
s
j
∗
(
α
)
:
[
L
i
j
k
l
:
D
k
l
∗
−
D
i
k
∗
⋅
σ
k
j
+
σ
i
k
⋅
D
k
j
∗
]
=
(
μ
i
j
(
α
)
−
ω
i
j
(
α
)
)
:
[
L
i
j
k
l
:
D
k
l
∗
−
D
i
k
∗
⋅
σ
k
j
+
σ
i
k
⋅
D
k
j
∗
]
\begin{aligned} \dot\tau^{(\alpha)}&=\bm{{m^{*(\alpha)}_i s^{*(\alpha)}_j}}:\left[L_{ijkl} :D^*_{kl}-{D_{ik}^*\cdot\sigma_{kj}+\sigma_{ik}\cdot D^*_{kj}}\right]\\ &=(\mu_{ij}^{(\alpha)}-\omega_{ij}^{(\alpha)} ):\left[L_{ijkl} :D^*_{kl}-{D_{ik}^*\cdot\sigma_{kj}+\sigma_{ik}\cdot D^*_{kj}}\right]\\ \end{aligned}
τ˙(α)=mi∗(α)sj∗(α):[Lijkl:Dkl∗−Dik∗⋅σkj+σik⋅Dkj∗]=(μij(α)−ωij(α)):[Lijkl:Dkl∗−Dik∗⋅σkj+σik⋅Dkj∗]
其中：
μ
i
j
(
α
)
:
L
i
j
k
l
:
D
k
l
∗
=
μ
k
l
(
α
)
:
L
k
l
i
j
:
D
i
j
∗
\mu_{ij}^{(\alpha)}:L_{ijkl} :D^*_{kl}=\mu_{kl}^{(\alpha)}:L_{klij} :D^*_{ij}
μij(α):Lijkl:Dkl∗=μkl(α):Lklij:Dij∗
又因为：
L
i
j
k
l
=
L
k
l
i
j
L_{ijkl}=L_{klij}\,
Lijkl=Lklij
μ
i
j
(
α
)
:
L
i
j
k
l
:
D
k
l
∗
=
(
L
i
j
k
l
:
μ
k
l
(
α
)
)
:
D
i
j
∗
\mu_{ij}^{(\alpha)}:L_{ijkl} :D^*_{kl}=(L_{ijkl} :\mu_{kl}^{(\alpha)}):D^*_{ij}
μij(α):Lijkl:Dkl∗=(Lijkl:μkl(α)):Dij∗
其中：
ω
i
j
(
α
)
:
L
i
j
k
l
:
D
k
l
∗
=
ω
k
l
(
α
)
:
L
k
l
i
j
:
D
i
j
∗
\omega_{ij}^{(\alpha)}:L_{ijkl} :D^*_{kl}=\omega_{kl}^{(\alpha)}:L_{klij} :D^*_{ij}
ωij(α):Lijkl:Dkl∗=ωkl(α):Lklij:Dij∗
又因为：
L
i
j
k
l
=
L
k
l
i
j
=
L
i
j
l
k
L_{ijkl}=L_{klij}=L_{ijlk}\,
Lijkl=Lklij=Lijlk且反对称张量
ω
k
l
=
−
ω
l
k
\omega_{kl}=-\omega_{lk}
ωkl=−ωlk
ω
i
j
(
α
)
:
L
i
j
k
l
:
D
k
l
∗
=
(
L
i
j
k
l
:
ω
k
l
(
α
)
)
:
D
i
j
∗
=
(
L
i
j
l
k
:
ω
l
k
(
α
)
)
:
D
i
j
∗
=
−
(
L
i
j
k
l
:
ω
k
l
(
α
)
)
:
D
i
j
∗
\begin{aligned} \omega_{ij}^{(\alpha)}:L_{ijkl} :D^*_{kl}&=(L_{ijkl} :\omega_{kl}^{(\alpha)}):D^*_{ij}=(L_{ijlk} :\omega_{lk}^{(\alpha)}):D^*_{ij}\\&=-(L_{ijkl} :\omega_{kl}^{(\alpha)}):D^*_{ij} \end{aligned}
ωij(α):Lijkl:Dkl∗=(Lijkl:ωkl(α)):Dij∗=(Lijlk:ωlk(α)):Dij∗=−(Lijkl:ωkl(α)):Dij∗
(
L
i
j
k
l
:
ω
k
l
(
α
)
)
:
D
i
j
∗
=
−
(
L
i
j
k
l
:
ω
k
l
(
α
)
)
:
D
i
j
∗
(L_{ijkl} :\omega_{kl}^{(\alpha)}):D^*_{ij}=-(L_{ijkl} :\omega_{kl}^{(\alpha)}):D^*_{ij}
(Lijkl:ωkl(α)):Dij∗=−(Lijkl:ωkl(α)):Dij∗
因此：
ω
i
j
(
α
)
:
L
i
j
k
l
:
D
k
l
∗
=
0
\omega_{ij}^{(\alpha)}:L_{ijkl} :D^*_{kl}=0
ωij(α):Lijkl:Dkl∗=0
继续计算共转项：
(
μ
i
j
(
α
)
−
ω
i
j
(
α
)
)
:
[
−
D
i
k
∗
⋅
σ
k
j
+
σ
i
k
⋅
D
k
j
∗
]
(\mu_{ij}^{(\alpha)}-\omega_{ij}^{(\alpha)} ):\left[-{D_{ik}^*\cdot\sigma_{kj}+\sigma_{ik}\cdot D^*_{kj}}\right]
(μij(α)−ωij(α)):[−Dik∗⋅σkj+σik⋅Dkj∗]
考虑到
μ
i
j
(
α
)
D
i
k
∗
σ
k
l
\mu_{ij}^{(\alpha)}\;D_{ik}^*\;\sigma_{kl}\,
μij(α)Dik∗σkl均为对称常量，考虑对称性可得
−
μ
i
j
(
α
)
:
(
D
i
k
∗
⋅
σ
k
j
)
=
−
μ
i
k
(
α
)
⋅
σ
k
j
:
D
i
j
∗
-\mu_{ij}^{(\alpha)}:(D_{ik}^*\cdot\sigma_{kj})=-\mu_{ik}^{(\alpha)}\cdot\sigma_{kj}:D_{ij}^*
−μij(α):(Dik∗⋅σkj)=−μik(α)⋅σkj:Dij∗
μ
i
j
(
α
)
:
(
σ
i
k
⋅
D
k
j
∗
)
=
μ
i
k
(
α
)
⋅
σ
k
j
:
D
i
j
∗
\mu_{ij}^{(\alpha)}:(\sigma_{ik}\cdot D_{kj}^*)=\mu_{ik}^{(\alpha)}\cdot\sigma_{kj}:D_{ij}^*
μij(α):(σik⋅Dkj∗)=μik(α)⋅σkj:Dij∗
μ
i
j
(
α
)
:
[
−
D
i
k
∗
⋅
σ
k
j
+
σ
i
k
⋅
D
k
j
∗
]
=
0
\mu_{ij}^{(\alpha)}:\left[-{D_{ik}^*\cdot\sigma_{kj}+\sigma_{ik}\cdot D^*_{kj}}\right]=0
μij(α):[−Dik∗⋅σkj+σik⋅Dkj∗]=0
张量
ω
i
j
(
α
)
\omega_{ij}^{(\alpha)}\,
ωij(α)为反对称张量：
ω
i
j
(
α
)
:
(
D
i
k
∗
⋅
σ
k
j
)
=
ω
i
k
(
α
)
⋅
σ
k
j
:
D
i
j
∗
\omega_{ij}^{(\alpha)}:(D_{ik}^*\cdot\sigma_{kj})=\omega_{ik}^{(\alpha)}\cdot\sigma_{kj}:D_{ij}^*
ωij(α):(Dik∗⋅σkj)=ωik(α)⋅σkj:Dij∗
−
ω
i
j
(
α
)
:
(
σ
i
k
⋅
D
k
j
∗
)
=
−
σ
i
k
⋅
ω
k
j
(
α
)
:
D
i
j
∗
=
ω
j
k
(
α
)
⋅
σ
i
k
:
D
i
j
∗
-\omega_{ij}^{(\alpha)}:(\sigma_{ik}\cdot D_{kj}^*)=-\sigma_{ik}\cdot\omega_{kj}^{(\alpha)}: D_{ij}^*=\omega_{jk}^{(\alpha)}\cdot\sigma_{ik}: D_{ij}^*
−ωij(α):(σik⋅Dkj∗)=−σik⋅ωkj(α):Dij∗=ωjk(α)⋅σik:Dij∗
−
ω
i
j
(
α
)
:
[
−
D
i
k
∗
⋅
σ
k
j
+
σ
i
k
⋅
D
k
j
∗
]
=
[
ω
i
k
(
α
)
⋅
σ
k
j
+
ω
j
k
(
α
)
⋅
σ
i
k
]
:
D
i
j
∗
-\omega_{ij}^{(\alpha)}:\left[-{D_{ik}^*\cdot\sigma_{kj}+\sigma_{ik}\cdot D^*_{kj}}\right] =[\omega_{ik}^{(\alpha)}\cdot\sigma_{kj}+\omega_{jk}^{(\alpha)}\cdot\sigma_{ik}]: D_{ij}^*
−ωij(α):[−Dik∗⋅σkj+σik⋅Dkj∗]=[ωik(α)⋅σkj+ωjk(α)⋅σik]:Dij∗
弹性应变增量=总应变-塑性应变：
Δ
ε
i
j
e
=
Δ
ε
i
j
−
Δ
ε
i
j
P
=
Δ
ε
i
j
−
∑
β
μ
i
j
(
β
)
Δ
γ
(
β
)
=
D
∗
Δ
t
\begin{aligned} \Delta\varepsilon_{ij}^e&=\Delta\varepsilon_{ij}-\Delta\varepsilon_{ij}^P\\ &=\Delta\varepsilon_{ij}-\sum_\beta\mu_{ij}^{(\beta)}\Delta\gamma^{(\beta)}\\ &=D^*\,\Delta t \end{aligned}
Δεije=Δεij−ΔεijP=Δεij−β∑μij(β)Δγ(β)=D∗Δt
由此得出（3.2.3）。
共转的应力增量
Δ
σ
i
j
=
σ
∇
⋅
Δ
t
\Delta\sigma_{ij}=\overset{\nabla}{\sigma}\cdot\Delta t\,
Δσij=σ∇⋅Δt，结合（2.2.2）和（3.2.1c）可得：
Δ
σ
i
j
=
[
σ
∇
i
j
∗
+
σ
i
k
⋅
ω
k
j
(
α
)
γ
˙
(
α
)
−
ω
i
k
(
α
)
⋅
σ
k
j
γ
˙
(
α
)
]
Δ
t
=
L
i
j
k
l
:
(
Δ
ε
i
j
−
Δ
ε
i
j
P
)
−
σ
i
j
Δ
ε
k
k
+
∑
α
[
σ
i
k
⋅
ω
k
j
(
α
)
−
ω
i
k
(
α
)
⋅
σ
k
j
]
Δ
γ
(
α
)
=
L
i
j
k
l
:
Δ
ε
i
j
−
σ
i
j
Δ
ε
k
k
−
∑
α
[
L
i
j
k
l
:
μ
k
l
+
ω
j
k
(
α
)
⋅
σ
k
i
+
ω
i
k
(
α
)
⋅
σ
k
j
]
Δ
γ
(
α
)
\begin{aligned} \Delta\sigma_{ij}&=[\overset{\nabla }{\sigma}^*_{ij}+\sigma_{ik}\cdot\omega^{(\alpha)}_{kj}\dot\gamma^{(\alpha)} -\omega^{(\alpha)}_{ik} \cdot\sigma_{kj}\dot\gamma^{(\alpha)} ]\Delta t\\ &=L_{ijkl}:(\Delta\varepsilon_{ij}-\Delta\varepsilon^P_{ij}) -\sigma_{ij}\Delta\varepsilon_{kk}+\sum_\alpha[\sigma_{ik}\cdot \omega^{(\alpha)}_{kj} -\omega^{(\alpha)}_{ik}\cdot\sigma_{kj}] \Delta\gamma^{(\alpha)}\\ &=L_{ijkl}:\Delta\varepsilon_{ij}-\sigma_{ij}\Delta\varepsilon_{kk}-\sum_\alpha[L_{ijkl}: \mu_{kl} +\omega^{(\alpha)}_{jk} \cdot\sigma_{ki} +\omega^{(\alpha)}_{ik}\cdot\sigma_{kj}] \Delta\gamma^{(\alpha)}\\ \end{aligned}
Δσij=[σ∇ij∗+σik⋅ωkj(α)γ˙(α)−ωik(α)⋅σkjγ˙(α)]Δt=Lijkl:(Δεij−ΔεijP)−σijΔεkk+α∑[σik⋅ωkj(α)−ωik(α)⋅σkj]Δγ(α)=Lijkl:Δεij−σijΔεkk−α∑[Lijkl:μkl+ωjk(α)⋅σki+ωik(α)⋅σkj]Δγ(α)
利用线性差值函数，通过给定的应变增量
Δ
ε
i
j
\Delta \varepsilon_{ij}\,
Δεij求出特定滑移系上的剪切应变
Δ
γ
(
α
)
\Delta\gamma^{(\alpha)}
Δγ(α)，将增量关系（3.2.2）和（3.2.3）代入（3.1.4）：
∑
β
{
δ
α
β
+
θ
Δ
t
∂
γ
˙
(
α
)
∂
τ
(
α
)
[
L
i
j
k
l
μ
k
l
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
]
μ
i
j
(
β
)
−
θ
Δ
t
∂
γ
˙
(
α
)
∂
g
(
α
)
h
α
β
sign
(
γ
˙
(
α
)
)
}
Δ
γ
(
β
)
=
γ
˙
(
α
)
Δ
t
+
θ
Δ
t
∂
γ
˙
(
α
)
∂
τ
(
α
)
[
L
i
j
k
l
μ
k
l
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
]
Δ
ε
i
j
(3.2.5)
\begin{aligned} &\sum_{\beta}\left\{\delta_{\alpha\beta} + \theta\,\Delta t \frac{\partial \dot\gamma^{(\alpha)}}{\partial \tau^{(\alpha)}} \left[L_{ijkl} \mu_{kl}^{(\alpha)} +\omega_{ik}^{(\alpha)} \sigma_{jk} +\omega_{jk}^{(\alpha)} \sigma_{ik} \right]\mu_{ij}^{(\beta)} -\theta\,\Delta t \frac{\partial \dot\gamma^{(\alpha)}}{\partial g^{(\alpha)}} h_{\alpha\beta}\,\operatorname{sign}\left( \dot\gamma^{(\alpha)} \right) \right\} \Delta \gamma^{(\beta)}\\ &=\dot\gamma^{(\alpha)}\Delta t+\theta\,\Delta t \frac{\partial \dot\gamma^{(\alpha)}}{\partial \tau^{(\alpha)}} \left[L_{ijkl} \mu_{kl}^{(\alpha)} +\omega_{ik}^{(\alpha)} \sigma_{jk} +\omega_{jk}^{(\alpha)} \sigma_{ik} \right]\Delta \varepsilon_{ij}\tag{3.2.5} \end{aligned}
β∑{δαβ+θΔt∂τ(α)∂γ˙(α)[Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik]μij(β)−θΔt∂g(α)∂γ˙(α)hαβsign(γ˙(α))}Δγ(β)=γ˙(α)Δt+θΔt∂τ(α)∂γ˙(α)[Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik]Δεij(3.2.5)
其中
δ
α
β
\delta_{\alpha\beta}
δαβ为Kronecker delta。一旦通过应变增量
Δ
ε
i
j
\Delta \varepsilon_{ij}\,
Δεij求出特定滑移系上的剪切应变
Δ
γ
(
α
)
\Delta\gamma^{(\alpha)}
Δγ(α)，其他增量通过（3.2.2）-（3.2.4）求出。
补充：式（3.2.5）中的 ∂ γ ˙ ( α ) ∂ τ ( α ) \large{\frac{\partial \dot\gamma^{(\alpha)}}{\partial \tau^{(\alpha)}}}\, ∂τ(α)∂γ˙(α)和 ∂ γ ˙ ( α ) ∂ g ( α ) \large{\frac{\partial \dot\gamma^{(\alpha)}}{\partial g^{(\alpha)}}}\, ∂g(α)∂γ˙(α)是通过2.3部分的率相关的滑移本构求出的。
3.3 晶格旋转
在晶体变形过程中，晶格经受畸变和旋转；然而，第2章中当所有的率方程都是建立在旋转坐标系上（2.1.2ab）晶体旋转没有明显包含在本构方程中（Asaro and Rice,1977 和 Asaro,1983b 也是如此）。在变形后的构型中，晶格变形和旋转通过相互垂直的向量，
s
∗
(
α
)
s^{*(\alpha)}\,
s∗(α)滑移方向向量和
m
∗
(
α
)
m^{*(\alpha)}\,
m∗(α)滑移面法向向量表示。对方程（2.1.2 a，b）微分：
s
˙
∗
(
α
)
=
(
D
∗
+
Ω
∗
)
⋅
s
∗
(
α
)
(3.3.1a)
\dot s^{*(\alpha)}=(D^*+\Omega^*)\cdot s^{*(\alpha)}\tag{3.3.1a}
s˙∗(α)=(D∗+Ω∗)⋅s∗(α)(3.3.1a)
m
˙
∗
(
α
)
=
−
m
∗
(
α
)
⋅
(
D
∗
+
Ω
∗
)
(3.3.1b)
\dot m^{*(\alpha)}=-m^{*(\alpha)}\cdot(D^*+\Omega^*)\tag{3.3.1b}
m˙∗(α)=−m∗(α)⋅(D∗+Ω∗)(3.3.1b)
详细推导可见(2.2.4)的补充推导
相应的增量形式用应变增量
Δ
ε
i
j
\Delta\varepsilon_{ij}\,
Δεij和所有滑移系上的滑移增量
Δ
γ
(
α
)
\Delta \gamma^{(\alpha)}
Δγ(α)表示：
Δ
s
i
∗
(
α
)
=
{
Δ
ε
i
j
+
Ω
i
j
Δ
t
−
∑
β
[
μ
i
j
(
β
)
+
ω
i
j
(
β
)
]
Δ
γ
(
β
)
}
s
j
∗
(
α
)
(3.3.2a)
\Delta s^{*(\alpha)}_i=\left\{\Delta\varepsilon_{ij} +\Omega_{ij}\Delta t - \sum_\beta \left[ \mu_{ij}^{(\beta)} + \omega_{ij}^{(\beta)} \right] \Delta \gamma^{(\beta)} \right\}s^{*(\alpha)}_j \tag{3.3.2a}
Δsi∗(α)=⎩⎨⎧Δεij+ΩijΔt−β∑[μij(β)+ωij(β)]Δγ(β)⎭⎬⎫sj∗(α)(3.3.2a)
Δ
m
i
∗
(
α
)
=
−
m
j
∗
(
α
)
{
Δ
ε
j
i
+
Ω
j
i
Δ
t
−
∑
β
[
μ
j
i
(
β
)
+
ω
j
i
(
β
)
]
Δ
γ
(
β
)
}
(3.3.2b)
\Delta m^{*(\alpha)}_i=-m^{*(\alpha)}_j \left\{\Delta\varepsilon_{ji} +\Omega_{ji}\Delta t - \sum_\beta \left[ \mu_{ji}^{(\beta)} + \omega_{ji}^{(\beta)} \right] \Delta \gamma^{(\beta)} \right\}\tag{3.3.2b}
Δmi∗(α)=−mj∗(α)⎩⎨⎧Δεji+ΩjiΔt−β∑[μji(β)+ωji(β)]Δγ(β)⎭⎬⎫(3.3.2b)
Δ
s
i
∗
(
α
)
\Delta s^{*(\alpha)}_i\,
Δsi∗(α)和
Δ
m
i
∗
(
α
)
\Delta m^{*(\alpha)}_i\,
Δmi∗(α)每一步都要更新，以得到在式（3.2.1a,b）中的定义的当前构型的“Schmid 因子”
μ
i
j
(
α
)
\mu_{ij}^{(\alpha)}
μij(α) 和张量
ω
i
j
(
α
)
\omega_{ij}^{(\alpha)}
ωij(α)。
3.4 非线性增量方程
在这一节中，滑移系的滑移应变
γ
(
α
)
\gamma^{(\alpha)}
γ(α)增量方程（3.1.2）不是通过Taylor展开式（3.1.3）求解的。在3.2和3.3节中，除了方程（3.2.5）所有增量方程在这一节中适用，且在时间增量步末尾计算应力和状态变量增量。滑移系上的滑移应变增量
Δ
γ
(
β
)
\Delta \gamma^{(\beta)}
Δγ(β)的求解方法从线性方程（3.2.5）换成了非线性方程，非线性方程是将滑移率的通用表达式（2.3.1）代入增量方程（3.1.2）：
Δ
γ
(
α
)
−
(
1
−
θ
)
Δ
t
γ
˙
t
(
α
)
−
θ
Δ
t
a
˙
(
α
)
f
(
α
)
(
τ
t
α
+
Δ
τ
α
g
t
α
+
Δ
g
α
)
=
0
(3.4.1)
\Delta \gamma^{(\alpha)} -(1-\theta)\,\Delta t\, \dot \gamma_t^{(\alpha)}-\theta\,\Delta t \, \dot a^{(\alpha)} f^{(\alpha)} \left(\frac{\tau_t^{\alpha}+\Delta \tau^{\alpha}}{g_t^{\alpha}+\Delta g^{\alpha}} \right)=0 \tag{3.4.1}
Δγ(α)−(1−θ)Δtγ˙t(α)−θΔta˙(α)f(α)(gtα+Δgατtα+Δτα)=0(3.4.1)
其中，
Δ
γ
(
α
)
\Delta \gamma^{(\alpha)}
Δγ(α)与分切应力增量
Δ
τ
(
α
)
\Delta \tau^{(\alpha)}
Δτ(α)的关系式（3.2.3）以及与当前强度
Δ
g
(
α
)
\Delta g^{(\alpha)}
Δg(α)的关系式为非线性的。利用Newton-Rhapson迭代法求解以上非线性方程，线性解（3.2.5）为初始试探值。所有其他增量的求解采取相同的积分过程。
4. ABAQUS子程序
4.1 UMAT
为了在abaqus实现第3章中向前梯度积分的单晶塑性本构，编写umat子程序。子程序中包含小变形理论和有限变形和有限旋转理论。附录A讨论UMAT子程序的input文件，单晶圆棒拉伸的input文件实例在附录C中。
单晶UMAT子程序的解相关状态变量包含所有滑移系上的当前强度 g ( α ) g^{(\alpha)} g(α)，剪切应变 γ ( α ) \gamma^{(\alpha)} γ(α)， 分切应力 τ ( α ) \tau^{(\alpha)} τ(α)，滑移面法向 m ( ∗ α ) m^{(*\alpha)} m(∗α)，滑移方向 s ( ∗ α ) s^{(*\alpha)} s(∗α)以及累计滑移量 γ \gamma γ。解相关状态变量的输出形式详见附录B。应力，应变和状态变量通过ABAQUS求解。当调用子程序时，主程序提供了增量步起始时的状态(应力，解相关的状态变量)和（预估的）应变增量以及时间增量。UMAT子程序执行2个函数：在时间步末尾更新应力以及解相关状态变量，为本构模型提供材料的Jacobian矩阵 ∂ Δ σ / ∂ Δ ε \partial \Delta \sigma/\partial \Delta \varepsilon ∂Δσ/∂Δε。Jacobian矩阵依赖于第三章向前梯度时间积分，因为单晶模型是率相关的，在子程序进行数值积分。
子程序提供两种应力和解相关状态变量的更新方案，一种是3.1-3.3节线性方案的在增量步起始通过线性方程进行更新，或者是3.4节非线性方案通过Newton-Rhapson迭代在增量末端进行更新。当增量方程是稳定的非线性方案允许较长的时间增量步。在Newton-Rhapson迭代方案中，对Jacobian矩阵 ∂ Δ σ / ∂ Δ ε \partial \Delta \sigma / \partial \Delta \varepsilon ∂Δσ/∂Δε 进行简化，忽略滑移面法向和滑移方向增量对应变增量的导数 ∂ Δ m ∗ ∂ Δ ε \large\frac{\partial \Delta m^*}{\partial \Delta \varepsilon} ∂Δε∂Δm∗ 和 ∂ Δ s ∗ ∂ Δ ε \large\frac{\partial \Delta s^*}{\partial \Delta \varepsilon} ∂Δε∂Δs∗。当不考虑晶格旋转简化不产生误差。如果考虑晶格旋转，误差是弹性应变增量的量级，与1比为 O ( D ∗ Δ t ) O(D^* \Delta t) O(D∗Δt)。
在当前版本abaqus主程序和umat子程序的接口中，没有提供方程（3.3.2a,b）中的增量
Ω
i
j
Δ
t
\Omega_{ij} \Delta t \,
ΩijΔt，可以通过旋转增量矩阵
Δ
R
\Delta R
ΔR求出：
Δ
R
=
(
I
−
1
2
Ω
Δ
t
)
−
1
⋅
(
I
+
1
2
Ω
Δ
t
)
(4.1.1)
\bm{\Delta R=\left(I-\frac{1}{2} \Omega \Delta t \right)^{-1}\cdot\left(I+\frac{1}{2} \Omega \Delta t \right)} \tag{4.1.1}
ΔR=(I−21ΩΔt)−1⋅(I+21ΩΔt)(4.1.1)
(abaqus 理论手册，1989)
虽然当前版本的UMAT是针对立方晶体的，但是也可以用于非立方晶体，将在下一节讨论。对每一种立方晶体子程序可以接受三组滑移系。根据观察BCC晶体有三组激活的滑移系{110}<111>，{121}<111>，{123}<111>，FCC有一组滑移系{111}<110>。
在umat子程序的主程序下包含如下函数子程序，F，DFDX，HSELF，HLATNT，GSLP0，DHSELF和DHLATN。这些子程序表征晶体滑移和硬化。子程序F通过式（2.3.1）求出增量步开始时的滑移率
γ
˙
(
α
)
\dot \gamma^{(\alpha)}
γ˙(α)，子程序DFDX提供F的导数
d
γ
˙
(
α
)
d
τ
(
α
)
/
g
(
α
)
\large\frac{d \dot \gamma^{(\alpha)}}{d \tau^{(\alpha)}/g^{(\alpha)}}
dτ(α)/g(α)dγ˙(α)。子程序HSELF和HLATNT提供自硬化和潜硬化系数。定义的具体形式详见附录A。
子程序GSLP0提供滑移系初始强度
g
(
α
)
(
0
)
g^{(\alpha)}(0)
g(α)(0)，定义为屈服强度
τ
0
\tau_0
τ0。子程序DHSELF和DHLATN只在Newton-Rhapson迭代方法中用到，给出自硬化和潜硬化模量对滑移的导数，
d
h
α
η
d
γ
(
β
)
\large \frac{dh_{\alpha\eta}}{d \gamma^{(\beta)}}
dγ(β)dhαη。
相同滑移物理参数相同。
更换不同的滑移率方程，只需更改子程序F和DFDX即可。
用户需要提供给UMAT子程序的7组数据：
（1）立方晶体弹性模量；
（2）潜在激活的滑移系的数量；
\quad\quad
典型的滑移面，如，BCC晶体的（110）
\quad\quad
典型的滑移方向，如，BCC晶体的[111]
（3）晶体的初始方向；
（4）滑移率方程中的常数，如，参考应变率
a
˙
(
α
)
\dot a^{(\alpha)}
a˙(α)，以及幂指数n；
（5）自硬化和潜硬化模量；
（6）向前梯度时间积分参数
θ
\theta
θ;且参数NLGEOM表示是用小变形理论还是有限变形和有限旋转理论进行分析；
（7）迭代方法所需参数；
输入文件的详细结构见附录A。
UMAT主程序中的8个子过程，ROTATION，SLIPSYS，STRAINRATE，LATENTHARDEN，GSLIPINIT，ITERATION，LUDCMP和LUBKSB。前5个子过程与UMAT主程序的关系如图1。子程序ROTATION定义立方晶体在整体坐标下的初始方向，SLIPSYS是在参考坐标下所有滑移系（滑移方向和滑移面法向）。子程序STRAINRATE，调用函数子程序F和DFDX，计算增量步开始时的各个滑移系的滑移率。如果需要迭代求解函数子程序F也要被UMAT主程序调用。子程序LATENTHARDEN调用函数子程序HSELF和HLATNT，已生成硬化矩阵。子程序GSLPINIT调用函数子程序GSLP0，计算各个滑移系在参考状态下的当前强度。子程序ITERATION，调用函数子程序DHSELF和DHLATN，提供Newton-Rhapson迭代方法。最后两个子程序用于求解线性方程；LUDCMP做LU分解，LUBKSB完成反向代入。
LU分解求解线性方程组
1.高斯消元法（Gaussian elimination）最常用求解线性方程组的方法，核心是将系数矩阵化为三角矩阵；
2.反向代入法(back substitution)：
方程组：
y
=
[
U
]
x
y=[U]{x}
y=[U]x，其中系数矩阵
U
U
U 为上三角矩阵
[
u
11
u
12
u
13
u
14
0
u
22
u
23
u
24
0
0
u
33
u
34
0
0
0
u
44
]
{
x
1
x
2
x
3
x
4
}
=
{
u
11
x
1
+
u
12
x
2
+
u
13
x
3
+
u
14
x
4
u
22
x
2
+
u
23
x
3
+
u
24
x
4
u
33
x
3
+
u
44
x
4
u
44
x
4
}
=
{
y
1
y
2
y
3
y
4
}
\left[ {\begin{matrix} u_{11}& u_{12}& u_{13}& u_{14} \\ 0& u_{22}& u_{23}& u_{24} \\ 0& 0& u_{33}& u_{34} \\0& 0&0& u_{44} \end{matrix}}\right] \left\{ {\begin{matrix} x_1 \\ x_2\\ x_3 \\ x_4\end{matrix}}\right\} =\left\{ {\begin{aligned} &u_{11}x_1+u_{12}x_2+u_{13}x_3+u_{14}x_4 \\ &u_{22}x_{2}+u_{23}x_3+u_{24}x_4\\ &u_{33}x_3+u_{44}x_4 \\ &u_{44}x_4\end{aligned}}\right\} =\left\{ {\begin{matrix} y_1 \\ y_2\\ y_3 \\ y_4\end{matrix}}\right\}
⎣⎢⎢⎡u11000u12u2200u13u23u330u14u24u34u44⎦⎥⎥⎤⎩⎪⎪⎨⎪⎪⎧x1x2x3x4⎭⎪⎪⎬⎪⎪⎫=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧u11x1+u12x2+u13x3+u14x4u22x2+u23x3+u24x4u33x3+u44x4u44x4⎭⎪⎪⎪⎪⎬⎪⎪⎪⎪⎫=⎩⎪⎪⎨⎪⎪⎧y1y2y3y4⎭⎪⎪⎬⎪⎪⎫
⟹
{
x
4
=
y
4
/
u
44
x
3
=
(
y
3
−
u
34
x
4
)
/
u
33
x
2
=
(
y
2
−
u
23
x
3
−
u
24
x
4
)
/
u
22
x
1
=
(
y
1
−
u
12
x
2
−
u
13
x
3
−
u
14
x
4
)
/
u
11
\Longrightarrow \begin{cases} x_4=y_4/u_{44} \\x_3=(y_3-u_{34}x_4)/u_{33} \\x_2=(y_2-u_{23}x_3-u_{24}x_4)/u_{22} \\x_1=(y_1-u_{12}x_2-u_{13}x_3-u_{14}x_4)/u_{11} \end{cases}
⟹⎩⎪⎪⎪⎨⎪⎪⎪⎧x4=y4/u44x3=(y3−u34x4)/u33x2=(y2−u23x3−u24x4)/u22x1=(y1−u12x2−u13x3−u14x4)/u11
3. LU分解
3.1 定义
把一个方阵A分解为下三角阵和上三角乘积，即为LU分解
A
=
L
U
A=LU
A=LU
3.2 LU分解的意义
求解线性系统
y
=
A
x
y=Ax
y=Ax，
y
=
L
U
x
y=LUx
y=LUx，记
z
=
U
x
z=Ux
z=Ux，通过
y
=
L
z
y=Lz
y=Lz 求出
z
z
z，最后通过
z
=
U
x
z=Ux
z=Ux，求出
x
x
x。
3.3 如何获得L和U
利用消元法：
[
A
11
A
12
A
13
A
14
A
21
A
22
A
23
A
24
A
31
A
32
A
33
A
34
A
41
A
42
A
43
A
44
]
→
L
1
[
A
11
A
12
A
13
A
14
0
A
22
−
A
12
A
21
A
11
A
23
−
A
13
A
21
A
11
A
24
−
A
14
A
21
A
11
0
A
32
−
A
12
A
31
A
11
A
33
−
A
13
A
31
A
11
A
34
−
A
13
A
31
A
11
0
A
42
−
A
12
A
41
A
11
A
43
−
A
13
A
41
A
11
A
44
−
A
13
A
41
A
11
]
\left[ {\begin{matrix} A_{11}& A_{12}& A_{13}& A_{14} \\ A_{21}& A_{22}& A_{23}& A_{24} \\ A_{31}& A_{32}& A_{33}& A_{34} \\ A_{41}& A_{42}& A_{43}& A_{44} \end{matrix}}\right] \overset{L_1}{\to } \left[ {\begin{matrix} A_{11}& A_{12}& A_{13}& A_{14} \\ 0& A_{22}-A_{12}\large \frac{A_{21}}{A_{11}}& A_{23}-A_{13}\large \frac{A_{21}}{A_{11}}& A_{24}-A_{14}\large \frac{A_{21}}{A_{11}} \\ 0& A_{32}-A_{12}\large \frac{A_{31}}{A_{11}}& A_{33}-A_{13}\large \frac{A_{31}}{A_{11}}& A_{34} -A_{13}\large \frac{A_{31}}{A_{11}}\\ 0& A_{42}-A_{12}\large \frac{A_{41}}{A_{11}}& A_{43}-A_{13}\large \frac{A_{41}}{A_{11}}& A_{44}-A_{13}\large \frac{A_{41}}{A_{11}} \end{matrix}}\right]
⎣⎢⎢⎡A11A21A31A41A12A22A32A42A13A23A33A43A14A24A34A44⎦⎥⎥⎤→L1⎣⎢⎢⎢⎡A11000A12A22−A12A11A21A32−A12A11A31A42−A12A11A41A13A23−A13A11A21A33−A13A11A31A43−A13A11A41A14A24−A14A11A21A34−A13A11A31A44−A13A11A41⎦⎥⎥⎥⎤
A
A
1
=
L
1
A
A\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad A^1=L_1A
AA1=L1A
L
1
=
[
1
0
0
0
−
A
21
A
11
1
0
0
−
A
31
A
11
0
1
0
−
A
41
A
11
0
0
1
]
=
[
1
0
0
0
−
L
2
,
1
1
0
0
−
L
3
,
1
0
1
0
−
L
4
,
1
0
0
1
]
L_1=\left[ {\begin{matrix} 1& 0& 0&0\\ -\frac{A_{21}}{A_{11}}& 1& 0&0\\ -\frac{A_{31}}{A_{11}}& 0& 1&0\\ -\frac{A_{41}}{A_{11}}& 0& 0&1 \end{matrix}}\right]=\left[ {\begin{matrix} 1& 0& 0&0\\ -L_{2,1}& 1& 0&0\\ -L_{3,1}& 0& 1&0\\ -L_{4,1}& 0& 0&1 \end{matrix}}\right]
L1=⎣⎢⎢⎡1−A11A21−A11A31−A11A41010000100001⎦⎥⎥⎤=⎣⎢⎢⎡1−L2,1−L3,1−L4,1010000100001⎦⎥⎥⎤
同理：
[
A
11
A
12
A
13
A
14
0
A
22
1
A
23
1
A
24
1
0
A
32
1
A
33
1
A
34
1
0
A
42
1
A
43
1
A
44
1
]
→
L
2
[
A
11
A
12
A
13
A
14
0
A
22
1
A
23
1
A
24
1
0
0
A
33
1
−
A
23
1
A
32
1
A
22
1
A
43
1
−
A
24
1
A
32
1
A
22
1
0
0
A
43
1
−
A
23
1
A
42
1
A
22
1
A
44
1
−
A
24
1
A
42
1
A
22
1
]
\left[ {\begin{matrix} A_{11}& A_{12}& A_{13}& A_{14} \\ 0& A_{22}^1& A_{23}^1& A_{24}^1 \\ 0& A_{32}^1& A_{33}^1& A_{34} ^1\\ 0& A_{42}^1& A^1_{43}& A_{44}^1 \end{matrix}}\right] \overset{L_2}{\to } \left[ {\begin{matrix} A_{11}& A_{12}& A_{13}& A_{14} \\ 0& A_{22}^1& A_{23}^1& A_{24}^1 \\ 0&0& A^1_{33} -A^1_{23}\large \frac{A^1_{32}}{A^1_{22}}& A^1_{43} -A^1_{24}\large \frac{A^1_{32}}{A^1_{22}}\\ 0&0& A^1_{43}-A^1_{23}\large \frac{A^1_{42}}{A^1_{22}} & A^1_{44}-A^1_{24}\large \frac{A^1_{42}}{A^1_{22}} \end{matrix}}\right]
⎣⎢⎢⎡A11000A12A221A321A421A13A231A331A431A14A241A341A441⎦⎥⎥⎤→L2⎣⎢⎢⎢⎢⎡A11000A12A22100A13A231A331−A231A221A321A431−A231A221A421A14A241A431−A241A221A321A441−A241A221A421⎦⎥⎥⎥⎥⎤
A
1
A
2
=
L
2
A
A^1\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad A^2=L_2A
A1A2=L2A
L
2
=
[
1
0
0
0
0
1
0
0
0
−
A
32
1
A
22
1
1
0
0
−
A
42
1
A
22
1
0
1
]
=
[
1
0
0
0
0
1
0
0
0
−
L
3
,
2
1
0
0
−
L
4
,
2
0
1
]
L_2=\left[ {\begin{matrix} 1& 0& 0&0\\ 0& 1& 0&0\\ 0&\large{-\frac{A^1_{32}}{A^1_{22}}}& 1&0\\ 0&-\large{\frac{A^1_{42}}{A^1_{22}}}& 0&1 \end{matrix}}\right]=\left[ {\begin{matrix} 1& 0& 0&0\\ 0& 1& 0&0\\0& -L_{3,2}& 1&0\\ 0&-L_{4,2}& 0&1 \end{matrix}}\right]
L2=⎣⎢⎢⎢⎢⎡100001−A221A321−A221A42100100001⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎡100001−L3,2−L4,200100001⎦⎥⎥⎤
继续：
[
A
11
A
12
A
13
A
14
0
A
22
1
A
23
1
A
24
1
0
0
A
33
2
A
34
2
0
0
A
43
2
A
44
2
]
→
L
3
[
A
11
A
12
A
13
A
14
0
A
22
1
A
23
1
A
24
1
0
0
A
33
2
A
34
2
0
0
0
A
44
2
−
A
34
2
A
43
2
A
33
2
]
\left[ {\begin{matrix} A_{11}& A_{12}& A_{13}& A_{14} \\ 0& A_{22}^1& A_{23}^1& A_{24}^1 \\ 0& 0& A_{33}^2& A_{34} ^2\\ 0&0& A_{43}^2& A_{44}^2 \end{matrix}}\right] \overset{L_3}{\to } \left[ {\begin{matrix} A_{11}& A_{12}& A_{13}& A_{14} \\ 0& A_{22}^1& A_{23}^1& A_{24}^1 \\ 0& 0& A_{33}^2& A_{34} ^2\\ 0&0& 0& A^2_{44}-A^2_{34}\large \frac{A^2_{43}}{A^2_{33}} \end{matrix}}\right]
⎣⎢⎢⎡A11000A12A22100A13A231A332A432A14A241A342A442⎦⎥⎥⎤→L3⎣⎢⎢⎢⎡A11000A12A22100A13A231A3320A14A241A342A442−A342A332A432⎦⎥⎥⎥⎤
A
2
U
=
A
3
=
L
3
A
A^2\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad U=A^3=L_3A
A2U=A3=L3A
L
3
=
[
1
0
0
0
0
1
0
0
0
0
1
0
0
0
−
A
43
2
A
33
2
1
]
=
[
1
0
0
0
0
1
0
0
0
0
1
0
0
0
−
L
4
,
3
1
]
L_3=\left[ {\begin{matrix} 1& 0& 0&0\\ 0& 1& 0&0\\ 0&0& 1&0\\ 0&0&-\large{\frac{A^2_{43}}{A^2_{33}}}& 1 \end{matrix}}\right] =\left[ {\begin{matrix} 1& 0& 0&0\\ 0& 1& 0&0\\ 0& 0 & 1&0\\ 0& 0&-L_{4,3}&1 \end{matrix}}\right]
L3=⎣⎢⎢⎢⎡10000100001−A332A4320001⎦⎥⎥⎥⎤=⎣⎢⎢⎡10000100001−L4,30001⎦⎥⎥⎤
因此：
A
=
L
3
L
2
L
1
U
→
U
=
L
1
−
1
L
2
−
1
L
3
−
1
A
A=L_3L_2L_1U \to U=L_1^{-1}L_2^{-1}L_3^{-1}A
A=L3L2L1U→U=L1−1L2−1L3−1A；易证
L
i
−
1
=
−
L
i
L_i^{-1}=-L_i
Li−1=−Li
Fortran程序： 输入变量：矩阵
A
A
A，矩阵维数
N
N
N；输出变量：下三角矩阵
L
L
L，其逆矩阵
L
−
1
L^{-1}
L−1（L_1），上三角矩阵
U
U
U
Program Main
Implicit none
Integer,parameter:: N=4
Integer I
Real A(N,N),L(N,N),L_1(N,N),L_2(N,N)
DATA A/2,2,3,4,1,3,3,4,1,2,3,5,1,2,3,6/
CALL LU(N,A,L,L_1)
DO I=1,N
write(*,*) A(I,:)
ENDDO
write(*,*) "***********"
DO I=1,N
write(*,*) L(I,:)
ENDDO
write(*,*) "***********"
DO I=1,N
write(*,*) L_1(I,:)
ENDDO
End Program Main
Subroutine LU(N,A,L,L_1)
Integer I,J,K
Real A(N,N),L(N,N),L_1(N,N),L_2(N,N)
L=0
L_1=0
L_2=0
DO I=1,N
L(I,I)=1
L_1(I,I)=1
L_2(I,I)=1
ENDDO
DO I=1,N-1
DO J=I+1,N
L(J,I)=A(J,I)/A(I,I)
DO K=1,I
A(J,K)=0.0
END DO
DO K=I+1,N
A(J,K)=A(J,K)-A(I,K)*L(J,I)
END DO
END DO
END DO
DO J=2,N
L_1(J,1)=L(J,1)
ENDDO
DO I=2,N-1
L_2=0
DO K=1,N
L_2(K,K)=1
ENDDO
DO J=I+1,N
L_2(J,I)=L(J,I)
ENDDO
L_1= matmul(L_2,L_1 )
ENDDO
Return
End Subroutine LU
4.2 针对UMAT子程序修改ABAQUS输入文件
用户自定义材料UMAT作为材料参数写入ABAQUS输入文件。当输入文件的单元定义好后，接下来的步骤应该包含弹塑性单晶响应：
（1）在输入文件中定义单晶固体 * MATERIAL卡下必需定义 *USER MATERIAL卡。 *USER MATERIAL卡下存在两类参数，CONSTANTS（常数）和UNSYMM（应力应变采用非对称时启用）。在本模型中采用了大量第一类参量。在当前版本的UMAT中，第4.1节中定义的7组共160个变量。更多细节见附录A。
（2）接下来是*DEPVAR卡用户必需提供解相关状态变量的个数。这个数量等于9倍的独立滑移系个数NSLPTL加5，即9*NSLPTL+5。正如2.3节讨论的，滑移系
(
−
s
∗
(
α
)
,
m
∗
(
α
)
)
(-s^{*(\alpha)},m^{*(\alpha)})
(−s∗(α),m∗(α)) 和滑移
(
s
∗
(
α
)
,
m
∗
(
α
)
)
(s^{*(\alpha)},m^{*(\alpha)})
(s∗(α),m∗(α)) 在立方晶体中不认为是独立的。每个滑移系有9个独立变量，分别是滑移系当前强度
g
(
α
)
g^{(\alpha)}
g(α),剪切应变
γ
(
α
)
\gamma^{(\alpha)}
γ(α)，分切应力
τ
(
α
)
\tau^{(\alpha)}
τ(α)，滑移面法向
m
(
α
)
m^{(\alpha)}
m(α)，滑移方向
s
(
α
)
s^{(\alpha)}
s(α)。所有滑移系的累计滑移应变
γ
\gamma
γ 也被认为是解相关状态变量。对于FCC金属晶体解相关状态变量为113个（=9*12+5）。
（3）在UMAT子程序之后必需定义*USER SUBROUTINE卡。
（4）为包含单晶有限应变和有限旋转，用户必需将参数NLGEOM设为非零值。因此同时，用户必需在.INP文件的*STEP卡内设置几何非线性。
输入文件中的7组材料参数共160个分布在20个数据卡中，每8个参数一个卡：三个卡记录晶体弹性刚度，四个卡记录潜在的激活滑移系，两个初始晶体方向，三个卡记录参考滑移率，六个记录自硬化系数和潜硬化系数，一个记录积分参数，还有一个记录积分方案。进一步详情见附录A。
5. 改进和提升
只要不改变输入数据卡片数据的分布，修改UMAT子程序相对简单。接下来给出几条改进意见：
5.1 非立方晶体
改变ROTATION和SLIPSYS，必需包含非立方晶体的长宽比和相对方向的非正交性。例如，正交晶体的[110]方向不一定垂直于[-110]方向，也不一定于(-111)垂直。
5.2 其他的滑移和硬化模型
UMAT子程序很容易包含其他的单晶模型。为了方便，讨论一种非幂律的滑移率表达式，只需简单改变子函数F和DFDX。然而子程序STRAINRATE也要修改成更一般的形式：
γ
˙
(
α
)
=
a
˙
(
α
)
f
(
α
)
(
τ
(
α
)
,
g
(
α
)
,
T
)
\dot\gamma^{(\alpha)}=\dot a^{(\alpha)}f^{(\alpha)}(\tau^{(\alpha)},g^{(\alpha)},T)
γ˙(α)=a˙(α)f(α)(τ(α),g(α),T)其中
T
T
T 为其他内变量。
5.3非schmid效应
来自Asaro and Rice（1977）