FCC单晶塑性abaqus-vumat
以黄永刚老师的umat文章为基础，构建单晶塑性的vumat。本文是之前“单晶塑性abaqus-umat 自学”的后续，大量引用前文已推导出的结论，用（umat）代指上文。求解过程针对FCC。
0. 整体思路
先从最简单的入手：1) 弹性各向同性；2)选用向前梯度插值，先不要Newton迭代法求解；3) 滑移只考虑Schmid准则；4) 求解方程不是通过LUP分解和反向代入法，而是采用求逆矩阵的方式（逆矩阵求解方法见上一篇文章）。求逆矩阵的方法写起来相对简洁，但计算量大。
1. 弹性各向同性
弹性采用各向同性，且分析过程中将体积变化和畸变分开考虑，动态分析中可根据状态方程(EOS)单独给出P-V关系。
体变： Δ σ V = K Δ ε v (1.1a) \Delta\sigma_V=K\Delta\varepsilon_v\tag{1.1a} ΔσV=KΔεv(1.1a)
等待偏应力更新之后，再在主应力上叠加体应力
σ
i
=
σ
i
+
Δ
σ
V
i
=
1
,
2
,
3
(1.1b)
\sigma_i=\sigma_i+\Delta\sigma_V\quad i=1,2,3\tag{1.1b}
σi=σi+ΔσVi=1,2,3(1.1b)
客观应力率：(后续推导中不特意说明就是偏应力和偏应变)
σ
∇
∗
=
L
:
D
∗
(1.2a)
\overset{\nabla}{\bm{\sigma}}^*=\bm{L:D^*}\tag{1.2a}
σ∇∗=L:D∗(1.2a)
σ ˙ − Ω ∗ ⋅ σ + σ ⋅ Ω ∗ = L : D ∗ (1.2b) \bm{ \dot\sigma-\Omega^*\cdot\sigma+\sigma\cdot\Omega^*}=\bm{L:D^*}\tag{1.2b} σ˙−Ω∗⋅σ+σ⋅Ω∗=L:D∗(1.2b)
方程两边同乘以
Δ
t
\Delta t
Δt
Δ
σ
−
Ω
∗
⋅
σ
Δ
t
+
σ
⋅
Ω
∗
Δ
t
=
L
:
Δ
ε
∗
(1.2c)
\bm{ \Delta\sigma-\Omega^*\cdot\sigma\Delta t+\sigma\cdot\Omega^*\Delta t}=\bm{L:\Delta \varepsilon^*}\tag{1.2c}
Δσ−Ω∗⋅σΔt+σ⋅Ω∗Δt=L:Δε∗(1.2c)
由于各向同性假设，2个弹性系数 E μ E\, \mu Eμ，且vumat给出的应变增量的剪切分量不是工程应变，因此相应的弹性矩阵为（偏应变与偏应力）： L i j k l = 2 G δ i k δ j l = [ 2 G 0 0 0 0 0 0 2 G 0 0 0 0 0 0 2 G 0 0 0 0 0 0 2 G 0 0 0 0 0 0 2 G 0 0 0 0 0 0 2 G ] {L_{ijkl}}=2G\delta_{ik}\delta_{jl}=\left[\begin{matrix} 2G&0&0&0&0&0\\0&2G&0&0&0&0\\0&0&2G&0&0&0\\0&0&0&2G&0&0\\0&0&0&0&2G&0\\0&0&0&0&0&2G \end{matrix} \right] Lijkl=2Gδikδjl=⎣⎢⎢⎢⎢⎢⎢⎡2G0000002G0000002G0000002G0000002G0000002G⎦⎥⎥⎥⎥⎥⎥⎤
自定义材料参数：props(1)=杨氏模量
E
E
E，props(2)=泊松比
μ
\mu
μ.
2. 12个滑移系
滑移面
{
111
}
\{111\}
{111}：滑移方向
<
1
1
ˉ
0
>
<
10
1
ˉ
>
<
01
1
ˉ
>
<1\bar10>\,<10\bar1>\,<01\bar1>
<11ˉ0><101ˉ><011ˉ>
滑移面
{
1
1
ˉ
1
}
\{1\bar11\}
{11ˉ1}：滑移方向
<
110
>
<
10
1
ˉ
>
<
011
>
<110>\,<10\bar1>\,<011>
<110><101ˉ><011>
滑移面
{
11
1
ˉ
}
\{11\bar1\}
{111ˉ}：滑移方向
<
1
1
ˉ
0
>
<
101
>
<
011
>
<1\bar10>\,<101>\,<011>
<11ˉ0><101><011>
滑移面
{
1
1
ˉ
1
ˉ
}
\{1\bar1\bar1\}
{11ˉ1ˉ}：滑移方向
<
110
>
<
101
>
<
01
1
ˉ
>
<110>\,<101>\,<01\bar1>
<110><101><011ˉ>
初始方向：三个欧拉角
φ
1
\varphi_1
φ1，
Φ
\Phi
Φ，
φ
2
\varphi_2
φ2，通过props(10) props(11) props(12)赋值
经过欧拉角转动后的晶体取向
g
g
g 为：
g
=
[
cos
φ
2
sin
φ
2
0
−
sin
φ
2
cos
φ
2
0
0
0
1
]
[
1
0
0
0
cos
Φ
sin
Φ
0
−
sin
Φ
cos
Φ
]
[
cos
φ
1
sin
φ
1
0
−
sin
φ
1
cos
φ
1
0
0
0
1
]
(2.1)
g=\left[ \begin{matrix} \cos \varphi_2 & \sin \varphi_2 & 0\\ -\sin \varphi_2 &\cos \varphi_2 & 0\\ 0&0&&1\end{matrix} \right] \left[ \begin{matrix}1&0&0\\ 0&\cos \varPhi &\sin \varPhi\\ 0&-\sin \varPhi &\cos \varPhi \end{matrix} \right] \left[ \begin{matrix} \cos \varphi_1 & \sin \varphi_1 & 0\\ -\sin \varphi_1 &\cos \varphi_1 & 0\\ 0&0&&1\end{matrix} \right]\tag{2.1}
g=⎣⎡cosφ2−sinφ20sinφ2cosφ20001⎦⎤⎣⎡1000cosΦ−sinΦ0sinΦcosΦ⎦⎤⎣⎡cosφ1−sinφ10sinφ1cosφ10001⎦⎤(2.1)
=
[
cos
φ
1
cos
φ
2
−
sin
φ
1
sin
φ
2
cos
Φ
sin
φ
1
cos
φ
2
+
cos
φ
1
sin
φ
2
cos
Φ
sin
φ
2
sin
Φ
−
cos
φ
1
sin
φ
2
−
sin
φ
1
cos
φ
2
cos
Φ
−
sin
φ
1
sin
φ
2
+
cos
φ
1
cos
φ
2
cos
Φ
cos
φ
2
sin
Φ
sin
φ
1
sin
Φ
−
cos
φ
1
sin
Φ
cos
Φ
]
R
D
:
<
u
v
w
>
N
D
:
{
h
k
l
}
=\left[ \begin{matrix} \cos \varphi_1 \cos \varphi_2- \sin \varphi_1 \sin \varphi_2\cos \varPhi & \sin \varphi_1 \cos \varphi_2+ \cos \varphi_1 \sin \varphi_2\cos \varPhi & \sin\varphi_2\sin\varPhi\\ -\cos \varphi_1 \sin \varphi_2- \sin \varphi_1 \cos \varphi_2\cos \varPhi & -\sin \varphi_1 \sin \varphi_2+ \cos \varphi_1 \cos \varphi_2\cos \varPhi & \cos\varphi_2\sin\varPhi\\ \sin \varphi_1 \sin \varPhi &-\cos\varphi_1\sin\varPhi&\cos\varPhi \end{matrix} \right]\\\quad\\ \qquad\qquad\qquad RD:<uvw>\qquad\qquad\qquad\quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad ND:\{hkl\}
=⎣⎡cosφ1cosφ2−sinφ1sinφ2cosΦ−cosφ1sinφ2−sinφ1cosφ2cosΦsinφ1sinΦsinφ1cosφ2+cosφ1sinφ2cosΦ−sinφ1sinφ2+cosφ1cosφ2cosΦ−cosφ1sinΦsinφ2sinΦcosφ2sinΦcosΦ⎦⎤RD:<uvw>ND:{hkl}
当初始欧拉角确定后，计算滑移系在整体坐标中表示
当前滑移系：
滑移面
m
=
{
h
k
l
}
=
{
1
3
1
3
1
3
}
⋅
g
(2.2a)
m=\{h\,k\,l\}=\{\frac{1}{\sqrt{3}}\,\frac{1}{\sqrt{3}}\,\frac{1}{\sqrt{3}}\} \cdot g\tag{2.2a}
m={hkl}={3
13
13
1}⋅g(2.2a)
滑移方向
s
=
<
u
v
w
>
=
<
1
2
−
1
2
0
>
⋅
g
(2.2b)
s=<u\,v\,w>=<\frac{1}{\sqrt{2}}\,\frac{-1}{\sqrt{2}}\,0>\cdot g\tag{2.2b}
s=<uvw>=<2
12
−10>⋅g(2.2b)
同理计算其他滑移系，与此同时对滑移系进行归一化。
计算出12个滑移系的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 ( α ) ] (2.3a) \mu_{ij}^{(\alpha)}=\frac{1}{2}\left[s_i^{(\alpha)}m_j^{(\alpha)}+s_j^{(\alpha)}m_i^{(\alpha)}\right] \tag{2.3a} μij(α)=21[si(α)mj(α)+sj(α)mi(α)](2.3a)
ω i j ( α ) = 1 2 [ s i ( α ) m j ( α ) − s j ( α ) m i ( α ) ] (2.3b) \omega_{ij}^{(\alpha)}=\frac{1}{2}\left[s_i^{(\alpha)}m_j^{(\alpha)}-s_j^{(\alpha)}m_i^{(\alpha)}\right] \tag{2.3b} ωij(α)=21[si(α)mj(α)−sj(α)mi(α)](2.3b)
其中: s i ( α ) m j ( α ) = [ s 1 m 1 s 1 m 2 s 1 m 3 s 2 m 1 s 2 m 2 s 2 m 3 s 3 m 1 s 3 m 2 s 3 m 3 ] = [ u h u k u l v h v k v l w h w k w l ] s_i^{(\alpha)}m_j^{(\alpha)}=\left[\begin{matrix} s_1\,m_1&s_1\,m_2&s_1\,m_3\\s_2\,m_1&s_2\,m_2&s_2\,m_3\\s_3\,m_1&s_3\,m_2&s_3\,m_3\\ \end{matrix}\right]=\left[\begin{matrix} uh&uk&ul\\vh&vk&vl\\wh&wk&wl \end{matrix}\right] si(α)mj(α)=⎣⎡s1m1s2m1s3m1s1m2s2m2s3m2s1m3s2m3s3m3⎦⎤=⎣⎡uhvhwhukvkwkulvlwl⎦⎤
μ i j ( α ) = [ s 1 m 1 s 1 m 2 + s 2 m 1 2 s 1 m 3 + s 3 m 1 2 s 1 m 2 + s 2 m 1 2 s 2 m 2 s 2 m 3 + s 3 m 2 2 s 1 m 3 + s 3 m 1 2 s 2 m 3 + s 3 m 2 2 s 3 m 3 ] \mu_{ij}^{(\alpha)}=\left[\begin{matrix} s_1\,m_1&\frac{s_1m_2+s_2m_1}{2}&\frac{s_1m_3+s_3m_1}{2}\\ \frac{s_1m_2+s_2m_1}{2}&s_2\,m_2&\frac{s_2m_3+s_3m_2}{2}\\ \frac{s_1m_3+s_3m_1}{2}&\frac{s_2m_3+s_3m_2}{2}&s_3\,m_3\\ \end{matrix}\right] μij(α)=⎣⎡s1m12s1m2+s2m12s1m3+s3m12s1m2+s2m1s2m22s2m3+s3m22s1m3+s3m12s2m3+s3m2s3m3⎦⎤
ω i j ( α ) = [ 0 s 1 m 2 − s 2 m 1 2 s 1 m 3 − s 3 m 1 2 s 2 m 1 − s 1 m 2 2 0 s 2 m 3 − s 3 m 2 2 s 3 m 1 − s 1 m 3 2 s 3 m 2 − s 2 m 3 2 0 ] \omega_{ij}^{(\alpha)}=\left[\begin{matrix} 0&\frac{s_1m_2-s_2m_1}{2}&\frac{s_1m_3-s_3m_1}{2}\\ \frac{s_2m_1-s_1m_2}{2}&0&\frac{s_2m_3-s_3m_2}{2}\\ \frac{s_3m_1-s_1m_3}{2}&\frac{s_3m_2-s_2m_3}{2}&0\\ \end{matrix}\right] ωij(α)=⎣⎡02s2m1−s1m22s3m1−s1m32s1m2−s2m102s3m2−s2m32s1m3−s3m12s2m3−s3m20⎦⎤
展开：
μ
I
(
α
)
=
[
s
1
m
1
s
2
m
2
s
3
m
3
1
2
(
s
1
m
2
+
s
2
m
1
)
1
2
(
s
2
m
3
+
s
3
m
2
)
1
2
(
s
1
m
3
+
s
3
m
1
)
]
ω
I
(
α
)
=
[
1
2
(
s
1
m
2
−
s
2
m
1
)
1
2
(
s
2
m
3
−
s
3
m
2
)
1
2
(
s
1
m
3
−
s
3
m
1
)
]
\mu_{I}^{(\alpha)}=\left[\begin{matrix} s_1\,m_1\\s_2\,m_2\\s_3\,m_3\\\frac{1}{2}{(s_1m_2+s_2m_1)}\\\frac{1}{2}{(s_2m_3+s_3m_2)}\\\frac{1}{2}{(s_1m_3+s_3m_1)}\\\end{matrix}\right]\quad \omega_{I}^{(\alpha)}=\left[\begin{matrix} \frac{1}{2}{(s_1m_2-s_2m_1)}\\\frac{1}{2}{(s_2m_3-s_3m_2)}\\\frac{1}{2}{(s_1m_3-s_3m_1)}\\\end{matrix}\right]\quad
μI(α)=⎣⎢⎢⎢⎢⎢⎢⎡s1m1s2m2s3m321(s1m2+s2m1)21(s2m3+s3m2)21(s1m3+s3m1)⎦⎥⎥⎥⎥⎥⎥⎤ωI(α)=⎣⎡21(s1m2−s2m1)21(s2m3−s3m2)21(s1m3−s3m1)⎦⎤
通过状态变量 STATEV(3*NSLPTL+1)~STATEV(6*NSLPTL)记录滑移面法向
m
m
m ，STATEV(6*NSLPTL+1)~STATEV(9*NSLPTL)记录相对应的滑移方向
s
s
s，NSLPTL为滑移系的数量（状态变量的顺序尽量和黄老师Umat中保持一致）。
针对第一步未开始变形时状态变量为了零，需要为滑移系赋初始值：
C***** 滑移系回归一化赋值后应约1，太小就是还没有赋值
If (StateOld(1) .EQ. 0.0)
C*********组装g矩阵
g(1,1)=cos(PROPS(10))*cos(PROPS(12))-sin(PROPS(10))*sin(PROPS(12))*cos(PROPS(11))
g(2,1)=-cos(PROPS(10))*sin(PROPS(11))-sin(PROPS(10))*cos(PROPS(12))*cos(PROPS(11))
g(3,1)=sin(PROPS(10))*sin(PROPS(11))
g(1,2)=sin(PROPS(10))*cos(PROPS(12))+cos(PROPS(10))*sin(PROPS(12))*cos(PROPS(11))
g(2,2)=-sin(PROPS(10))*sin(PROPS(12))+cos(PROPS(10))*cos(PROPS(12))*cos(PROPS(11))
g(3,2)=-cos(PROPS(10))*sin(PROPS(11))
g(1,3)=sin(PROPS(12))*sin(PROPS(11))
g(2,3)=cos(PROPS(12))*sin(PROPS(11))
g(3,3)=cos(PROPS(11))
C*********初始滑移面
Do I=1,3
StateOld(3*NSLPTL+I)=(1.0/Sqrt(3))*(g(1,I)+g(2,I)+g(3,I))
StateOld(3*NSLPTL+I+3)=StateOld(3*NSLPTL+I)
StateOld(3*NSLPTL+I+6)=StateOld(3*NSLPTL+I)
StateOld(3*NSLPTL+I+9)=(1.0/Sqrt(3))*(g(1,I)-g(2,I)+g(3,I))
StateOld(3*NSLPTL+I+12)=StateOld(3*NSLPTL+I)
StateOld(3*NSLPTL+I+15)=StateOld(3*NSLPTL+I)
StateOld(3*NSLPTL+I+18)=(1.0/Sqrt(3))*(g(1,I)+g(2,I)-g(3,I))
StateOld(3*NSLPTL+I+21)=StateOld(3*NSLPTL+I)
StateOld(3*NSLPTL+I+24)=StateOld(3*NSLPTL+I)
StateOld(3*NSLPTL+I+27)=(1.0/Sqrt(3))*(g(1,I)-g(2,I)-g(3,I))
StateOld(3*NSLPTL+I+30)=StateOld(3*NSLPTL+I)
StateOld(3*NSLPTL+I+33)=StateOld(3*NSLPTL+I)
End Do
C*********初始滑移方向
Do I=1,3
StateOld(6*NSLPTL+I)=(1.0/Sqrt(2))*(g(1,I)-g(2,I))
StateOld(6*NSLPTL+I+3)=(1.0/Sqrt(2))*(g(1,I)-g(3,I))
StateOld(6*NSLPTL+I+6)=(1.0/Sqrt(2))*(g(2,I)-g(3,I))
StateOld(6*NSLPTL+I+9)=(1.0/Sqrt(2))*(g(1,I)+g(2,I))
StateOld(6*NSLPTL+I+12)=(1.0/Sqrt(2))*(g(1,I)-g(3,I))
StateOld(6*NSLPTL+I+15)=(1.0/Sqrt(2))*(g(2,I)+g(3,I))
StateOld(6*NSLPTL+I+18)=(1.0/Sqrt(2))*(g(1,I)-g(2,I))
StateOld(6*NSLPTL+I+21)=(1.0/Sqrt(2))*(g(1,I)+g(3,I))
StateOld(6*NSLPTL+I+24)=(1.0/Sqrt(2))*(g(2,I)+g(3,I))
StateOld(6*NSLPTL+I+27)=(1.0/Sqrt(2))*(g(1,I)+g(2,I))
StateOld(6*NSLPTL+I+30)=(1.0/Sqrt(2))*(g(1,I)+g(3,I))
StateOld(6*NSLPTL+I+33)=(1.0/Sqrt(2))*(g(2,I)-g(3,I))
End Do
End If
C*****组装“schmid”因子mu(12,6)和旋转张量omega(12,3)
Do I=1,NSLPTL
mu(I,1)=StateOld(3*NSLPTL+3*I-2)*StateOld(6*NSLPTL+3*I-2)
mu(I,2)=StateOld(3*NSLPTL+3*I-1)*StateOld(6*NSLPTL+3*I-1)
mu(I,3)=StateOld(3*NSLPTL+3*I)*StateOld(6*NSLPTL+3*I)
mu(I,4)=0.5*(StateOld(3*NSLPTL+3*I-2)*StateOld(6*NSLPTL+3*I-1)
1 +StateOld(3*NSLPTL+3*I-1)*StateOld(6*NSLPTL+3*I-2))
mu(I,5)=0.5*(StateOld(3*NSLPTL+3*I)*StateOld(6*NSLPTL+3*I-1)
1 +StateOld(3*NSLPTL+3*I-1)*StateOld(6*NSLPTL+3*I))
mu(I,6)=0.5*(StateOld(3*NSLPTL+3*I)*StateOld(6*NSLPTL+3*I-2)
1 +StateOld(3*NSLPTL+3*I-2)*StateOld(6*NSLPTL+3*I))
omega(I,1)=0.5*(StateOld(3*NSLPTL+3*I-1)*StateOld(6*NSLPTL+3*I-2)
1 -StateOld(3*NSLPTL+3*I-2)*StateOld(6*NSLPTL+3*I-1))
omega(I,2)=0.5*(StateOld(3*NSLPTL+3*I)*StateOld(6*NSLPTL+3*I-1)
1 -StateOld(3*NSLPTL+3*I-1)*StateOld(6*NSLPTL+3*I))
omega(I,3)=0.5*(StateOld(3*NSLPTL+3*I)*StateOld(6*NSLPTL+3*I-2)
1 -StateOld(3*NSLPTL+3*I-2)*StateOld(6*NSLPTL+3*I))
End Do
3. 旋转
在umat中给出了
Δ
R
\Delta R
ΔR DROT, 需要通过
F
=
R
U
→
F
U
−
1
=
R
F=RU\to FU^{-1}=R
F=RU→FU−1=R求出新老旧两个
R
R
R，通过做差求出
Δ
R
\Delta R
ΔR：
Δ
R
=
R
N
e
w
−
R
O
l
d
(3.1)
\Delta R=R^{New}-R^{Old}\tag{3.1}
ΔR=RNew−ROld(3.1)
F
=
[
F
1
F
4
F
9
F
7
F
2
F
5
F
6
F
8
F
3
]
U
=
[
U
1
U
4
U
6
U
4
U
2
U
5
U
6
U
5
U
3
]
F=\left[\begin{matrix} F_1&F_4&F_9\\F_7&F_2&F_5\\F_6&F_8&F_3 \end{matrix}\right] \quad U=\left[\begin{matrix} U_1&U_4&U_6\\U_4&U_2&U_5\\U_6&U_5&U_3 \end{matrix}\right]
F=⎣⎡F1F7F6F4F2F8F9F5F3⎦⎤U=⎣⎡U1U4U6U4U2U5U6U5U3⎦⎤
det
(
U
)
=
U
1
U
2
U
3
+
2
U
4
U
5
U
6
−
U
1
U
5
U
5
−
U
3
U
4
U
4
−
U
2
U
6
U
6
\det(U)=U_1U_2U_3+2U_4U_5U_6-U_1U_5U_5-U_3U_4U_4-U_2U_6U_6
det(U)=U1U2U3+2U4U5U6−U1U5U5−U3U4U4−U2U6U6
U
U
U 的伴随
U
∗
U^*
U∗
U
∗
=
[
U
2
U
3
−
U
5
U
5
U
5
U
6
−
U
3
U
4
U
4
U
5
−
U
2
U
6
U
5
U
6
−
U
3
U
4
U
1
U
3
−
U
6
U
6
U
4
U
6
−
U
1
U
5
U
4
U
5
−
U
2
U
6
U
4
U
6
−
U
1
U
5
U
1
U
2
−
U
4
U
4
]
U^*=\left[\begin{matrix} U_2U_3-U_5U_5& U_5U_6-U_3U_4& U_4U_5-U_2U_6\\ U_5U_6-U_3U_4& U_1U_3-U_6U_6& U_4U_6-U_1U_5\\ U_4U_5-U_2U_6& U_4U_6-U_1U_5& U_1U_2-U_4U_4 \end{matrix}\right]
U∗=⎣⎡U2U3−U5U5U5U6−U3U4U4U5−U2U6U5U6−U3U4U1U3−U6U6U4U6−U1U5U4U5−U2U6U4U6−U1U5U1U2−U4U4⎦⎤
DetOld=stretchOld(NB,1)*stretchOld(NB,2)*stretchOld(NB,3)
1 +2.0*stretchOld(NB,4)*stretchOld(NB,5)*stretchOld(NB,6)
2 -stretchOld(NB,1)*stretchOld(NB,5)*stretchOld(NB,5)
3 -stretchOld(NB,3)*stretchOld(NB,4)*stretchOld(NB,4)
4 -stretchOld(NB,2)*stretchOld(NB,6)*stretchOld(NB,6)
ROld(1,1)=(1.0/DetOld)*(defgradOld(NB,1)*
1 (stretchOld(NB,2)*stretchOld(NB,3)-stretchOld(NB,5)*stretchOld(NB,5))
2 +defgradOld(NB,4)*
3 (stretchOld(NB,5)*stretchOld(NB,6)-stretchOld(NB,3)*stretchOld(NB,4))
4 +defgradOld(NB,9)*
5 (stretchOld(NB,4)*stretchOld(NB,5)-stretchOld(NB,2)*stretchOld(NB,6)))
ROld(2,1)=(1.0/DetOld)*(defgradOld(NB,7)*
1 (stretchOld(NB,2)*stretchOld(NB,3)-stretchOld(NB,5)*stretchOld(NB,5))
2 +defgradOld(NB,2)*
3 (stretchOld(NB,5)*stretchOld(NB,6)-stretchOld(NB,3)*stretchOld(NB,4))
4 +defgradOld(NB,5)*
5 (stretchOld(NB,4)*stretchOld(NB,5)-stretchOld(NB,2)*stretchOld(NB,6)))
ROld(3,1)=(1.0/DetOld)*(defgradOld(NB,6)*
1 (stretchOld(NB,2)*stretchOld(NB,3)-stretchOld(NB,5)*stretchOld(NB,5))
2 +defgradOld(NB,8)*
3 (stretchOld(NB,5)*stretchOld(NB,6)-stretchOld(NB,3)*stretchOld(NB,4))
4 +defgradOld(NB,3)*
5 (stretchOld(NB,4)*stretchOld(NB,5)-stretchOld(NB,2)*stretchOld(NB,6)))
ROld(1,2)=(1.0/DetOld)*(defgradOld(NB,1)*
1 (stretchOld(NB,5)*stretchOld(NB,6)-stretchOld(NB,3)*stretchOld(NB,4))
2 +defgradOld(NB,4)*
3 (stretchOld(NB,1)*stretchOld(NB,3)-stretchOld(NB,6)*stretchOld(NB,6))
4 +defgradOld(NB,9)*
5 (stretchOld(NB,4)*stretchOld(NB,6)-stretchOld(NB,1)*stretchOld(NB,5)))
ROld(2,2)=(1.0/DetOld)*(defgradOld(NB,7)*
1 (stretchOld(NB,5)*stretchOld(NB,6)-stretchOld(NB,3)*stretchOld(NB,4))
2 +defgradOld(NB,2)*
3 (stretchOld(NB,1)*stretchOld(NB,3)-stretchOld(NB,6)*stretchOld(NB,6))
4 +defgradOld(NB,5)*
5 (stretchOld(NB,4)*stretchOld(NB,6)-stretchOld(NB,1)*stretchOld(NB,5)))
ROld(3,2)=(1.0/DetOld)*(defgradOld(NB,6)*
1 (stretchOld(NB,5)*stretchOld(NB,6)-stretchOld(NB,3)*stretchOld(NB,4))
2 +defgradOld(NB,8)*
3 (stretchOld(NB,1)*stretchOld(NB,3)-stretchOld(NB,6)*stretchOld(NB,6))
4 +defgradOld(NB,3)*
5 (stretchOld(NB,4)*stretchOld(NB,6)-stretchOld(NB,1)*stretchOld(NB,5)))
ROld(1,3)=(1.0/DetOld)*(defgradOld(NB,1)*
1 (stretchOld(NB,4)*stretchOld(NB,5)-stretchOld(NB,2)*stretchOld(NB,6))
2 +defgradOld(NB,4)*
3 (stretchOld(NB,4)*stretchOld(NB,6)-stretchOld(NB,1)*stretchOld(NB,5))
4 +defgradOld(NB,9)*
5 (stretchOld(NB,1)*stretchOld(NB,2)-stretchOld(NB,4)*stretchOld(NB,4)))
ROld(2,3)=(1.0/DetOld)*(defgradOld(NB,7)*
1 (stretchOld(NB,4)*stretchOld(NB,5)-stretchOld(NB,2)*stretchOld(NB,6))
2 +defgradOld(NB,2)*
3 (stretchOld(NB,4)*stretchOld(NB,6)-stretchOld(NB,1)*stretchOld(NB,5))
4 +defgradOld(NB,5)*
5 (stretchOld(NB,1)*stretchOld(NB,2)-stretchOld(NB,4)*stretchOld(NB,4)))
ROld(3,3)=(1.0/DetOld)*(defgradOld(NB,6)*
1 (stretchOld(NB,4)*stretchOld(NB,5)-stretchOld(NB,2)*stretchOld(NB,6))
2 +defgradOld(NB,8)*
3 (stretchOld(NB,4)*stretchOld(NB,6)-stretchOld(NB,1)*stretchOld(NB,5))
4 +defgradOld(NB,3)*
5 (stretchOld(NB,1)*stretchOld(NB,2)-stretchOld(NB,4)*stretchOld(NB,4)))
Write(*,*)"ROld=",ROld(1,:)
Write(*,*)"ROld=",ROld(2,:)
Write(*,*)"ROld=",ROld(3,:)
DetNew=stretchNew(NB,1)*stretchNew(NB,2)*stretchNew(NB,3)
1 +2.0*stretchNew(NB,4)*stretchNew(NB,5)*stretchNew(NB,6)
2 -stretchNew(NB,1)*stretchNew(NB,5)*stretchNew(NB,5)
3 -stretchNew(NB,3)*stretchNew(NB,4)*stretchNew(NB,4)
4 -stretchNew(NB,2)*stretchNew(NB,6)*stretchNew(NB,6)
RNew(1,1)=(1.0/DetNew)*(defgradNew(NB,1)*
1 (stretchNew(NB,2)*stretchNew(NB,3)-stretchNew(NB,5)*stretchNew(NB,5))
2 +defgradNew(NB,4)*
3 (stretchNew(NB,5)*stretchNew(NB,6)-stretchNew(NB,3)*stretchNew(NB,4))
4 +defgradNew(NB,9)*
5 (stretchNew(NB,4)*stretchNew(NB,5)-stretchNew(NB,2)*stretchNew(NB,6)))
RNew(2,1)=(1.0/DetNew)*(defgradNew(NB,7)*
1 (stretchNew(NB,2)*stretchNew(NB,3)-stretchNew(NB,5)*stretchNew(NB,5))
2 +defgradNew(NB,2)*
3 (stretchNew(NB,5)*stretchNew(NB,6)-stretchNew(NB,3)*stretchNew(NB,4))
4 +defgradNew(NB,5)*
5 (stretchNew(NB,4)*stretchNew(NB,5)-stretchNew(NB,2)*stretchNew(NB,6)))
RNew(3,1)=(1.0/DetNew)*(defgradNew(NB,6)*
1 (stretchNew(NB,2)*stretchNew(NB,3)-stretchNew(NB,5)*stretchNew(NB,5))
2 +defgradNew(NB,8)*
3 (stretchNew(NB,5)*stretchNew(NB,6)-stretchNew(NB,3)*stretchNew(NB,4))
4 +defgradNew(NB,3)*
5 (stretchNew(NB,4)*stretchNew(NB,5)-stretchNew(NB,2)*stretchNew(NB,6)))
RNew(1,2)=(1.0/DetNew)*(defgradNew(NB,1)*
1 (stretchNew(NB,5)*stretchNew(NB,6)-stretchNew(NB,3)*stretchNew(NB,4))
2 +defgradNew(NB,4)*
3 (stretchNew(NB,1)*stretchNew(NB,3)-stretchNew(NB,6)*stretchNew(NB,6))
4 +defgradNew(NB,9)*
5 (stretchNew(NB,4)*stretchNew(NB,6)-stretchNew(NB,1)*stretchNew(NB,5)))
RNew(2,2)=(1.0/DetNew)*(defgradNew(NB,7)*
1 (stretchNew(NB,5)*stretchNew(NB,6)-stretchNew(NB,3)*stretchNew(NB,4))
2 +defgradNew(NB,2)*
3 (stretchNew(NB,1)*stretchNew(NB,3)-stretchNew(NB,6)*stretchNew(NB,6))
4 +defgradNew(NB,5)*
5 (stretchNew(NB,4)*stretchNew(NB,6)-stretchNew(NB,1)*stretchNew(NB,5)))
RNew(3,2)=(1.0/DetNew)*(defgradNew(NB,6)*
1 (stretchNew(NB,5)*stretchNew(NB,6)-stretchNew(NB,3)*stretchNew(NB,4))
2 +defgradNew(NB,8)*
3 (stretchNew(NB,1)*stretchNew(NB,3)-stretchNew(NB,6)*stretchNew(NB,6))
4 +defgradNew(NB,3)*
5 (stretchNew(NB,4)*stretchNew(NB,6)-stretchNew(NB,1)*stretchNew(NB,5)))
RNew(1,3)=(1.0/DetNew)*(defgradNew(NB,1)*
1 (stretchNew(NB,4)*stretchNew(NB,5)-stretchNew(NB,2)*stretchNew(NB,6))
2 +defgradNew(NB,4)*
3 (stretchNew(NB,4)*stretchNew(NB,6)-stretchNew(NB,1)*stretchNew(NB,5))
4 +defgradNew(NB,9)*
5 (stretchNew(NB,1)*stretchNew(NB,2)-stretchNew(NB,4)*stretchNew(NB,4)))
RNew(2,3)=(1.0/DetNew)*(defgradNew(NB,7)*
1 (stretchNew(NB,4)*stretchNew(NB,5)-stretchNew(NB,2)*stretchNew(NB,6))
2 +defgradNew(NB,2)*
3 (stretchNew(NB,4)*stretchNew(NB,6)-stretchNew(NB,1)*stretchNew(NB,5))
4 +defgradNew(NB,5)*
5 (stretchNew(NB,1)*stretchNew(NB,2)-stretchNew(NB,4)*stretchNew(NB,4)))
RNew(3,3)=(1.0/DetNew)*(defgradNew(NB,6)*
1 (stretchNew(NB,4)*stretchNew(NB,5)-stretchNew(NB,2)*stretchNew(NB,6))
2 +defgradNew(NB,8)*
3 (stretchNew(NB,4)*stretchNew(NB,6)-stretchNew(NB,1)*stretchNew(NB,5))
4 +defgradNew(NB,3)*
5 (stretchNew(NB,1)*stretchNew(NB,2)-stretchNew(NB,4)*stretchNew(NB,4)))
Write(*,*)"RNew=",RNew(1,:)
Write(*,*)"RNew=",RNew(2,:)
Write(*,*)"RNew=",RNew(3,:)
do I=1,3
do J=1,3
RInc(I,J)=RNew(I,J)-ROld(I,J)
end do
end do
Write(*,*)"RInc=",RInc(1,:)
Write(*,*)"RInc=",RInc(2,:)
Write(*,*)"RInc=",RInc(3,:)
再通过
Δ
R
\Delta R\,\,
ΔRDROT，求出
Ω
Δ
t
\Omega\Delta t\,\,
ΩΔtDSPIN：
Δ
R
=
(
I
−
1
2
Ω
Δ
t
)
−
1
⋅
(
I
+
1
2
Ω
Δ
t
)
(3.2a)
{\Delta R=\left(I-\frac{1}{2} \Omega \Delta t \right)^{-1}\cdot\left(I+\frac{1}{2} \Omega \Delta t \right)}\tag{3.2a}
ΔR=(I−21ΩΔt)−1⋅(I+21ΩΔt)(3.2a)
Ω
Δ
t
=
2
(
Δ
R
−
I
)
(
Δ
R
+
I
)
−
1
(3.2b)
{\Omega \Delta t =2\left(\Delta R-I \right)\left(\Delta R+I \right)^{-1}} \tag{3.2b}
ΩΔt=2(ΔR−I)(ΔR+I)−1(3.2b)
Δ
R
=
[
Δ
R
11
Δ
R
12
Δ
R
13
Δ
R
21
Δ
R
22
Δ
R
23
Δ
R
31
Δ
R
32
Δ
R
33
]
\Delta R=\left[\begin{matrix} \Delta R_{11}&\Delta R_{12}&\Delta R_{13}\\\Delta R_{21}&\Delta R_{22}&\Delta R_{23}\\\Delta R_{31}&\Delta R_{32}&\Delta R_{33} \end{matrix}\right]
ΔR=⎣⎡ΔR11ΔR21ΔR31ΔR12ΔR22ΔR32ΔR13ΔR23ΔR33⎦⎤
Δ
R
+
I
=
[
Δ
R
11
+
1
Δ
R
12
Δ
R
13
Δ
R
21
Δ
R
22
+
1
Δ
R
23
Δ
R
31
Δ
R
32
Δ
R
33
+
1
]
\Delta R+I=\left[\begin{matrix} \Delta R_{11}+1&\Delta R_{12}&\Delta R_{13}\\\Delta R_{21}&\Delta R_{22}+1&\Delta R_{23}\\\Delta R_{31}&\Delta R_{32}&\Delta R_{33}+1 \end{matrix}\right]
ΔR+I=⎣⎡ΔR11+1ΔR21ΔR31ΔR12ΔR22+1ΔR32ΔR13ΔR23ΔR33+1⎦⎤
(
Δ
R
+
I
)
=
[
Δ
R
11
+
1
Δ
R
12
Δ
R
13
Δ
R
21
Δ
R
22
+
1
Δ
R
23
Δ
R
31
Δ
R
32
Δ
R
33
+
1
]
(\Delta R+I)^{}=\left[\begin{matrix} \Delta R_{11}+1&\Delta R_{12}&\Delta R_{13}\\\Delta R_{21}&\Delta R_{22}+1&\Delta R_{23}\\\Delta R_{31}&\Delta R_{32}&\Delta R_{33}+1 \end{matrix}\right]
(ΔR+I)=⎣⎡ΔR11+1ΔR21ΔR31ΔR12ΔR22+1ΔR32ΔR13ΔR23ΔR33+1⎦⎤
4. 率相关的滑移硬化本构
基于Schmid准则，
α
\alpha
α 滑移系的滑移率
γ
˙
(
α
)
\dot\gamma^{(\alpha)}
γ˙(α) 由分切应力
τ
(
α
)
\tau^{(\alpha)}
τ(α)决定：
γ
˙
(
α
)
=
a
˙
f
(
α
)
(
τ
(
α
)
g
(
α
)
)
(4.1a)
\dot\gamma^{(\alpha)}=\dot a f^{(\alpha)}\left(\frac{\tau^{(\alpha)}}{g^{(\alpha)}} \right)\tag{4.1a}
γ˙(α)=a˙f(α)(g(α)τ(α))(4.1a)
其中， f ( α ) ( x ) = x ∣ x ∣ n − 1 (4.1b) f^{(\alpha)}(x)=x|x|^{n-1}\tag{4.1b} f(α)(x)=x∣x∣n−1(4.1b)
应变硬化： g ˙ ( α ) = ∑ β h α β γ ˙ ( β ) \dot g^{(\alpha)}=\sum_\beta h_{\alpha\beta}\dot\gamma^{(\beta)} g˙(α)=β∑hαβγ˙(β)
其中自硬化系数：
h
α
α
=
h
(
γ
)
=
h
0
sech
2
∣
h
0
γ
τ
s
−
τ
0
∣
(
α
不
求
和
)
(4.2a)
h_{\alpha\alpha}=h(\gamma)=h_0\operatorname{sech}^2\large{\left| \frac{h_0\gamma}{\tau_s-\tau_0} \right|}\quad (\alpha 不求和)\tag{4.2a}
hαα=h(γ)=h0sech2∣∣∣∣∣τs−τ0h0γ∣∣∣∣∣(α不求和)(4.2a)
其中
h
0
h_0
h0 为初始硬化模量，
τ
0
\tau_0
τ0屈服强度，等于当前强度的初始值
g
(
α
)
(
0
)
g^{(\alpha)}(0)
g(α)(0)，
γ
\gamma
γ 是所有滑移系上的Taylor累计剪切应变：
γ
=
∑
α
∫
0
t
∣
γ
˙
(
α
)
∣
d
t
(4.2b)
\gamma=\sum_\alpha\int_{0}^{t} \left|\dot{\gamma}^{(\alpha)}\right|dt\tag{4.2b}
γ=α∑∫0t∣∣∣γ˙(α)∣∣∣dt(4.2b)
潜硬化模量：
h
α
β
=
q
h
(
γ
)
(
α
≠
β
)
(4.2c)
h_{\alpha\beta}=qh(\gamma)\quad(\alpha \ne \beta)\tag{4.2c}
hαβ=qh(γ)(α=β)(4.2c)
自定义材料参数：props(3)=参考滑移率
a
˙
\dot a
a˙，props(4)=滑移率指数
n
n
n，props(5)=参考硬化系数
h
0
h_0
h0，props(6)=初始屈服强度
τ
0
\tau_0
τ0，props(7)=break-through突破强度
τ
s
\tau_s
τsprops(8)=潜硬化系数
q
q
q。
5. 求解12个滑移量
向前梯度插值求解黄老师umat中的公式（3.2.5）：
∑
β
{
δ
α
β
+
θ
Δ
t
∂
γ
˙
(
α
)
∂
τ
(
α
)
[
L
i
j
k
l
μ
k
l
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
]
μ
i
j
(
β
)
−
θ
Δ
t
∂
γ
˙
(
α
)
∂
g
(
α
)
h
α
β
sign
(
γ
˙
(
β
)
)
}
Δ
γ
(
β
)
=
γ
˙
t
(
α
)
Δ
t
+
θ
Δ
t
∂
γ
˙
(
α
)
∂
τ
(
α
)
[
L
i
j
k
l
μ
k
l
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
]
Δ
ε
i
j
转
化
为
⟹
A
A
⋅
Δ
γ
=
Y
(5.1)
\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^{(\beta)} \right) \right\} \Delta \gamma^{(\beta)}\\ &=\dot\gamma_t^{(\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{5.1}\\\quad\\ &\qquad 转化为\implies \bm{AA\cdot \Delta \gamma=Y} \end{aligned}
β∑{δαβ+θΔt∂τ(α)∂γ˙(α)[Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik]μij(β)−θΔt∂g(α)∂γ˙(α)hαβsign(γ˙(β))}Δγ(β)=γ˙t(α)Δt+θΔt∂τ(α)∂γ˙(α)[Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik]Δεij转化为⟹AA⋅Δγ=Y(5.1)props(9)=梯度差值系数
θ
\theta
θ
求解
∂
γ
˙
(
α
)
∂
τ
(
α
)
\Large\frac{\partial \dot\gamma^{(\alpha)}}{\partial \tau^{(\alpha)}}
∂τ(α)∂γ˙(α) 和
∂
γ
˙
(
α
)
∂
g
(
α
)
\Large\frac{\partial \dot\gamma^{(\alpha)}}{\partial g^{(\alpha)}}
∂g(α)∂γ˙(α)，对（4.1） 求导：
∂
γ
˙
(
α
)
∂
f
(
α
)
=
a
˙
d
f
(
α
)
d
x
=
n
∣
x
∣
n
−
1
其
中
x
=
τ
(
α
)
g
(
α
)
(5.2a)
\frac{\partial\dot\gamma^{(\alpha)}}{\partial f^{(\alpha)}}=\dot a \quad \frac{d f^{(\alpha)}}{dx}=n|x|^{n-1} \quad 其中\, x=\frac{\tau^{(\alpha)}}{g^{(\alpha)}} \tag{5.2a}
∂f(α)∂γ˙(α)=a˙dxdf(α)=n∣x∣n−1其中x=g(α)τ(α)(5.2a)
∂ x ∂ τ ( α ) = 1 g ( α ) ∂ x ∂ g ( α ) = − τ ( α ) g ( α ) 2 (5.2b) \frac{\partial x}{\partial \tau^{(\alpha)}}=\frac{1}{g^{(\alpha)}} \quad \frac{\partial x}{\partial g^{(\alpha)}}=\frac{-\tau^{(\alpha)}}{g^{(\alpha)2}} \tag{5.2b} ∂τ(α)∂x=g(α)1∂g(α)∂x=g(α)2−τ(α)(5.2b)
∂ γ ˙ ( α ) ∂ τ ( α ) = n a ˙ ∣ τ ( α ) g ( α ) ∣ n − 1 1 g ( α ) (5.2c) \frac{\partial \dot\gamma^{(\alpha)}}{\partial \tau^{(\alpha)}}=n\dot a \left|\frac{\tau^{(\alpha)}}{g^{(\alpha)}} \right|^{n-1} \frac{1}{g^{(\alpha)}} \tag{5.2c} ∂τ(α)∂γ˙(α)=na˙∣∣∣∣g(α)τ(α)∣∣∣∣n−1g(α)1(5.2c)
∂ γ ˙ ( α ) ∂ g ( α ) = n a ˙ ∣ τ ( α ) g ( α ) ∣ n − 1 − τ ( α ) g ( α ) 2 (5.2d) \frac{\partial \dot\gamma^{(\alpha)}}{\partial g^{(\alpha)}}=n\dot a \left|\frac{\tau^{(\alpha)}}{g^{(\alpha)}} \right|^{n-1} \frac{-{\tau^{(\alpha)}}}{g^{(\alpha)2}} \tag{5.2d} ∂g(α)∂γ˙(α)=na˙∣∣∣∣g(α)τ(α)∣∣∣∣n−1g(α)2−τ(α)(5.2d)
设
A
i
j
(
α
)
=
L
i
j
k
l
μ
k
l
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
A_{ij}^{(\alpha)}=L_{ijkl} \mu_{kl}^{(\alpha)} +\omega_{ik}^{(\alpha)} \sigma_{jk} +\omega_{jk}^{(\alpha)} \sigma_{ik}
Aij(α)=Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik：
L
i
j
k
l
μ
k
l
(
α
)
=
2
G
μ
i
j
(
α
)
L_{ijkl} \mu_{kl}^{(\alpha)}=2G\mu_{ij}^{(\alpha)}
Lijklμkl(α)=2Gμij(α)
ω i k ( α ) σ j k = [ 0 ω 1 ω 3 − ω 1 0 ω 2 − ω 3 − ω 2 0 ] [ σ 1 σ 4 σ 6 σ 4 σ 2 σ 5 σ 6 σ 5 σ 3 ] = [ ω 1 σ 4 + ω 3 σ 6 ω 1 σ 2 + ω 3 σ 5 ω 1 σ 5 + ω 3 σ 3 − ω 1 σ 1 + ω 2 σ 6 − ω 1 σ 4 + ω 2 σ 5 − ω 1 σ 6 + ω 2 σ 3 − ω 3 σ 1 − ω 2 σ 4 − ω 3 σ 4 − ω 2 σ 2 − ω 3 σ 6 − ω 2 σ 5 ] \begin{aligned} \omega_{ik}^{(\alpha)} \sigma_{jk} &=\left[\begin{matrix} 0&\omega_1&\omega_3\\-\omega_1&0&\omega_2\\-\omega_3&-\omega_2&0 \end{matrix} \right]\left[\begin{matrix} \sigma_1&\sigma_4&\sigma_6\\ \sigma_4&\sigma_2&\sigma_5\\ \sigma_6&\sigma_5&\sigma_3 \end{matrix} \right]\\ &=\left[\begin{matrix} \omega_1 \sigma_4+\omega_3 \sigma_6&\omega_1 \sigma_2+\omega_3 \sigma_5&\omega_1 \sigma_5+\omega_3 \sigma_3\\ -\omega_1 \sigma_1+\omega_2 \sigma_6&-\omega_1 \sigma_4+\omega_2 \sigma_5&-\omega_1\sigma_6+\omega_2 \sigma_3\\ -\omega_3 \sigma_1-\omega_2 \sigma_4&-\omega_3 \sigma_4-\omega_2 \sigma_2&-\omega_3 \sigma_6-\omega_2 \sigma_5 \end{matrix} \right]\end{aligned} ωik(α)σjk=⎣⎡0−ω1−ω3ω10−ω2ω3ω20⎦⎤⎣⎡σ1σ4σ6σ4σ2σ5σ6σ5σ3⎦⎤=⎣⎡ω1σ4+ω3σ6−ω1σ1+ω2σ6−ω3σ1−ω2σ4ω1σ2+ω3σ5−ω1σ4+ω2σ5−ω3σ4−ω2σ2ω1σ5+ω3σ3−ω1σ6+ω2σ3−ω3σ6−ω2σ5⎦⎤
ω j k ( α ) σ i k = σ i k ω j k ( α ) = σ i k ω k j ( α ) T = ( ω i k ( α ) σ k j ) T \omega_{jk}^{(\alpha)} \sigma_{ik}=\sigma_{ik}\omega_{jk}^{(\alpha)}=\sigma_{ik}\omega_{kj}^{(\alpha)T}=(\omega_{ik}^{(\alpha)} \sigma_{kj})^T ωjk(α)σik=σikωjk(α)=σikωkj(α)T=(ωik(α)σkj)T
A i j ( α ) = L i j k l μ k l ( α ) + ω i k ( α ) σ j k + ω j k ( α ) σ i k = [ 2 ( ω 1 σ 4 + ω 3 σ 6 ) + 2 G μ 1 ω 1 σ 2 + ω 3 σ 5 − ω 1 σ 1 + ω 2 σ 6 + 2 G μ 4 ω 1 σ 5 + ω 3 σ 3 − ω 3 σ 1 − ω 2 σ 4 + 2 G μ 6 . . 2 ( − ω 1 σ 4 + ω 2 σ 5 ) + 2 G μ 2 − ω 1 σ 6 + ω 2 σ 3 − ω 3 σ 4 − ω 2 σ 2 + 2 G μ 5 . . . . 2 ( − ω 3 σ 6 − ω 2 σ 5 ) + 2 G μ 3 ] \begin{aligned} &A_{ij}^{(\alpha)}=L_{ijkl} \mu_{kl}^{(\alpha)} +\omega_{ik}^{(\alpha)} \sigma_{jk} +\omega_{jk}^{(\alpha)} \sigma_{ik}\\ &=\left[\begin{matrix} 2(\omega_1 \sigma_4+\omega_3 \sigma_6) +2G\mu_1& \omega_1 \sigma_2+\omega_3 \sigma_5-\omega_1 \sigma_1+\omega_2 \sigma_6+2G\mu_4& \omega_1\sigma_5+\omega_3 \sigma_3-\omega_3 \sigma_1-\omega_2 \sigma_4+2G\mu_6\\ ..&2(-\omega_1 \sigma_4+\omega_2\sigma_5)+2G\mu_2&-\omega_1\sigma_6+\omega_2 \sigma_3-\omega_3 \sigma_4-\omega_2 \sigma_2+2G\mu_5\\ ..&..&2(-\omega_3 \sigma_6-\omega_2 \sigma_5)+2G\mu_3 \end{matrix} \right] \end{aligned} Aij(α)=Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik=⎣⎡2(ω1σ4+ω3σ6)+2Gμ1....ω1σ2+ω3σ5−ω1σ1+ω2σ6+2Gμ42(−ω1σ4+ω2σ5)+2Gμ2..ω1σ5+ω3σ3−ω3σ1−ω2σ4+2Gμ6−ω1σ6+ω2σ3−ω3σ4−ω2σ2+2Gμ52(−ω3σ6−ω2σ5)+2Gμ3⎦⎤
滑移系的当前强度
g
(
α
)
g^{(\alpha)}
g(α) STATE(1)~STATE(NSLPTL)，滑移率
γ
˙
α
\dot\gamma^{\alpha}
γ˙α STATE(1*NSLPTL+1)~STATE(2*NSLPTL)和当前的分切应力
τ
(
α
)
\tau^{(\alpha)}
τ(α)STATE(2*NSLPTL+1)~STATE(3*NSLPTL)，Taylor累计剪切应变STATE(10*NSPLTL+1)通过状态变量传递回来
伪程序：
NSLPTL=12 ! 12个滑移系
C*****计算 DgammaDtau(12) 和 DgammaDg(12), 引入中间值 DFDX（12）
Do I=1,NSLPTL
DFDX(I)=dtArray(1)*props(9)*props(3)*props(4)*(Abs(STATEV(2*NSLPTL+I)/STATEV(I))**(props(4)-1.0))
DgammaDtau(I)=DFDX(I)/STATEV(I)
DgammaDg(I)=DFDX(I)*(-STATEV(I+2*NSLPTL))/(STATEV(I)**2)
End Do
C*****计算自硬化系数Haa,潜硬化系数Hab
Haa=props(5)*(1.0/(cosh(props(5)*STATEV(10*NSPLTL+1)/(props(7)-props(6)))))**2.0
Hab=props(8)*Haa
C*****计算A(NSLPTL,6)
Do I=1,NSLPTL
A(I,1)=2.0*(omage(I,1)*stressOld(nblock,4)+omage(I,3)*stressOld(nblock,6)+EG*mu(I,1))
A(I,2)=2.0*(-omage(I,1)*stressOld(nblock,4)+omage(I,2)*stressOld(nblock,5)+EG*mu(I,2))
A(I,3)=2.0*(-omage(I,3)*stressOld(nblock,6)-omage(I,2)*stressOld(nblock,5)+EG*mu(I,3))
A(I,4)=2.0*EG*mu(I,4)+omage(I,1)*stressOld(nblock,2)+omage(I,3)*stressOld(nblock,5)
1 -omage(I,1)*stressOld(nblock,1)+omage(I,2)*stressOld(nblock,6)
A(I,5)=2.0*EG*mu(I,5)-omage(I,1)*stressOld(nblock,6)+omage(I,2)*stressOld(nblock,3)
1 -omage(I,3)*stressOld(nblock,4)-omage(I,2)*stressOld(nblock,2)
A(I,6)=2.0*EG*mu(I,6)+omage(I,1)*stressOld(nblock,5)+omage(I,3)*stressOld(nblock,3)
1 -omage(I,3)*stressOld(nblock,1)-omage(I,2)*stressOld(nblock,4)
End Do
C*****计算偏应变
strainIncV=0.0 !体积应变
Do I=1,ndir
strainIncV=strainIncV+strainInc(nblock,I)
End Do
Do I=1,ndir
strainIncD(I)=strainInc(nblock,I)-strainIncV
End Do
Do I=1+ndir,ndir+nshr
strainIncD(I)=strainInc(nblock,I)
End Do
C*****组装AA(NSLPTL,NSLPTL)和 Y(NSLPTL)
Do I=1,NSLPTL
Y(I)=0.0
Do K=1,3
Y(I)=Y(I)+A(I,K)*strainIncD(K)
End Do
Do K=4,6
Y(I)=Y(I)+2.0*A(I,K)*strainIncD(K)
End Do
Y(I)=dtArray(1)*STATEV(NSLPTL+I)+DgammaDtau(I)*Y(I)
Do J=1,NSLPTL
AA(I,J)=0.0
Do K=1,3
AA(I,J)=AA(I,J)+A(I,K)*mu(J,K)
End Do
Do K=4,6
AA(I,J)=AA(I,J)+2.0*A(I,K)*mu(J,K)
End Do
If (I .EQ. J) Then
AA(I,I)=1.0+DgammaDtau(I)*AA(I,J)-DgammaDg(I)*Haa*SIGN(1.0,STATEV(2*NSLPTL+I))
Else
AA(I,J)=DgammaDtau(I)*AA(I,J)-DgammaDg(I)*Hab*SIGN(1.0,STATEV(2*NSLPTL+I))
End If
End Do
End Do
C*****求解Dgamma
Call matInv(AA,NSLPTL,Det)
Dgamma=matmul(AA,Y)
6. 应力更新
黄老师umat中的式（3.2.4）：
Δ
σ
i
j
=
L
i
j
k
l
Δ
ε
k
l
−
σ
i
j
Δ
ε
k
k
−
∑
α
[
L
i
j
k
l
μ
k
l
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
]
⋅
Δ
γ
(
α
)
\Delta\sigma_{ij}=L_{ijkl}\Delta\varepsilon_{kl}-\sigma_{ij}\Delta\varepsilon_{kk}-\sum_\alpha \left[L_{ijkl}\mu^{(\alpha)}_{kl}+\omega^{(\alpha)}_{ik}\sigma_{jk}+\omega^{(\alpha)}_{jk}\sigma_{ik}\right]\cdot\Delta\gamma^{(\alpha)}
Δσij=LijklΔεkl−σijΔεkk−α∑[Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik]⋅Δγ(α)
改写：
Δ
σ
I
=
2
G
Δ
ε
I
−
∑
α
A
I
(
α
)
⋅
Δ
γ
(
α
)
\Delta\sigma_{I}=2G\Delta\varepsilon_{I}-\sum_\alpha A^{(\alpha)}_{I}\cdot\Delta\gamma^{(\alpha)}
ΔσI=2GΔεI−α∑AI(α)⋅Δγ(α)
叠加体应力：
Δ
σ
I
=
Δ
σ
I
+
E
3
(
1
−
2
μ
)
Δ
ε
V
I
=
1
,
N
D
I
\Delta\sigma_{I}=\Delta\sigma_{I}+\frac{E}{3(1-2\mu)}\Delta \varepsilon_V \qquad I=1,\mathrm{NDI}
ΔσI=ΔσI+3(1−2μ)EΔεVI=1,NDI
Do I=1,6
stressP=0.0
Do J=1,NSLPTL
stressP=stressP+A(J,I)*Dgamma(J)
End Do
stressNew(nblock,I)=stressOld(nblock,I)+2.0*G*strainIncD(I)-stressP
End Do
Do I=1,3
stressNew(nblock,I)=stressNew(nblock,I)+(props(1)/(3.0*(1.0-2.0*props(2))))*strainIncV
End Do
7. 更新状态变量
整理一下用到的状态变量：
| 状态变量 | 序号 |
|---|---|
| 滑移系的当前强度 g ( α ) g^{(\alpha)} g(α) | STATE(1)~STATE(NSLPTL) |
| 滑移率 γ ˙ α \dot\gamma^{\alpha} γ˙α | STATE(NSLPTL+1)~STATE(2*NSLPTL) |
| 当前的分切应力 τ ( α ) \tau^{(\alpha)} τ(α) | STATE(2*NSLPTL+1)~STATE(3*NSLPTL) |
| 滑移面法向 m m m | STATE(3*NSLPTL+1)~STATE(6*NSLPTL) |
| 滑移方向 s s s | STATE(6*NSLPTL+1)~STATE(9*NSLPTL) |
| 累积滑移量 ∑ ∣ Δ γ ( α ) ∣ \sum\begin{vmatrix}\Delta\gamma^{(\alpha)}\end{vmatrix} ∑∣∣Δγ(α)∣∣ | STATE(9*NSLPTL+1)~STATE(10*NSLPTL) |
| Taylor累积滑移量 ∑ ∣ Δ γ ∣ \sum\begin{vmatrix}\Delta\gamma\end{vmatrix} ∑∣∣Δγ∣∣ | STATE(10*NSLPTL+1) |
滑移系的当前强度
g
(
α
)
g^{(\alpha)}
g(α) 更新：
Δ
g
(
α
)
=
∑
β
h
α
β
Δ
γ
(
β
)
\Delta g^{(\alpha)}=\sum_\beta h_{\alpha\beta}\Delta\gamma^{(\beta)}
Δg(α)=β∑hαβΔγ(β)
Do I=1,NSLPTL
Do J=1,NSLPTL
If (I .EQ. J) Then
StaveNew(I)=StaveOld(I)+Haa*Dgamma(J)
Else
StaveNew(I)=StaveOld(I)+Hab*Dgamma(J)
End If
End Do
End Do
滑移率 γ ˙ α \dot\gamma^{\alpha} γ˙α 更新 γ ˙ α = Δ γ ( α ) Δ t \dot\gamma^{\alpha}=\Large \frac{\Delta \gamma^{(\alpha)}}{\Delta t} γ˙α=ΔtΔγ(α):
Do I=1,NSLPTL
StaveNew(NSLPTL+I)=Dgamma(I)/dtArray(1)
End Do
分切应力
τ
(
α
)
\tau^{(\alpha)}
τ(α) 更新：
Δ
τ
(
α
)
=
[
L
i
j
k
l
μ
k
l
(
α
)
+
ω
i
k
(
α
)
σ
j
k
+
ω
j
k
(
α
)
σ
i
k
]
:
[
Δ
ε
i
j
−
∑
β
μ
i
j
(
β
)
Δ
γ
(
β
)
]
\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]
Δτ(α)=[Lijklμkl(α)+ωik(α)σjk+ωjk(α)σik]:⎣⎡Δεij−β∑μij(β)Δγ(β)⎦⎤
改写：
Δ
τ
(
α
)
=
A
I
(
α
)
⋅
[
Δ
ε
I
−
∑
β
μ
I
(
β
)
Δ
γ
(
β
)
]
\Delta\tau^{(\alpha)}=A^{(\alpha)}_I \cdot\left[\Delta \varepsilon_{I}-\sum_\beta\mu_{I}^{(\beta)}\Delta\gamma^{(\beta)}\right]
Δτ(α)=AI(α)⋅⎣⎡ΔεI−β∑μI(β)Δγ(β)⎦⎤
Do I=1,6
strainIncP(I)=0.0
Do J=1,NSLPTL
strainIncP(I)=strainIncP(I)+A(J,I)*(mu(J,I)*Dgamma(J))
End Do
End Do
Do I=1,NSLPTL
DTau=0.0
Do J=1,6
Dtau=Dtau+A(I,J)*(strainIncD(J)-strainIncP(J))
End Do
StateNew(2*NSLPTL+I)=StateOld(2*NSLPTL+I)+Dtau
End Do
更新滑移系：
Δ
m
i
∗
(
α
)
=
−
m
j
∗
(
α
)
{
Δ
ε
j
i
+
Ω
j
i
Δ
t
−
∑
β
[
μ
j
i
(
β
)
+
ω
j
i
(
β
)
]
Δ
γ
(
β
)
}
\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\}
Δmi∗(α)=−mj∗(α)⎩⎨⎧Δεji+ΩjiΔt−β∑[μji(β)+ωji(β)]Δγ(β)⎭⎬⎫
Δ s i ∗ ( α ) = { Δ ε i j + Ω i j Δ t − ∑ β [ μ i j ( β ) + ω i j ( β ) ] Δ γ ( β ) } s j ∗ ( α ) \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 Δsi∗(α)=⎩⎨⎧Δεij+ΩijΔt−β∑[μij(β)+ωij(β)]Δγ(β)⎭⎬⎫sj∗(α)
设：
B
(
3
,
3
)
=
Δ
ε
i
j
+
Ω
i
j
Δ
t
−
∑
β
[
μ
i
j
(
β
)
+
ω
i
j
(
β
)
]
Δ
γ
(
β
)
B(3,3)=\Delta\varepsilon_{ij} +\Omega_{ij}\Delta t - \sum_\beta \left[ \mu_{ij}^{(\beta)} + \omega_{ij}^{(\beta)} \right] \Delta \gamma^{(\beta)}
B(3,3)=Δεij+ΩijΔt−∑β[μij(β)+ωij(β)]Δγ(β)
首先：
B
=
[
Δ
ε
1
Δ
ε
4
Δ
ε
6
Δ
ε
4
Δ
ε
2
Δ
ε
5
Δ
ε
6
Δ
ε
5
Δ
ε
3
]
+
[
Ω
11
Ω
12
Ω
13
Ω
21
Ω
22
Ω
23
Ω
31
Ω
32
Ω
33
]
Δ
t
B=\left[\begin{matrix} \Delta\varepsilon_{1} & \Delta\varepsilon_{4} & \Delta\varepsilon_{6} \\ \Delta\varepsilon_{4} & \Delta\varepsilon_{2} & \Delta\varepsilon_{5} \\ \Delta\varepsilon_{6} & \Delta\varepsilon_{5} & \Delta\varepsilon_{3} \end{matrix} \right] +\left[\begin{matrix} \Omega_{11} &\Omega_{12} & \Omega_{13} \\ \Omega_{21} &\Omega_{22} &\Omega_{23} \\ \Omega_{31} & \Omega_{32} &\Omega_{33} \end{matrix} \right]\Delta t
B=⎣⎡Δε1Δε4Δε6Δε4Δε2Δε5Δε6Δε5Δε3⎦⎤+⎣⎡Ω11Ω21Ω31Ω12Ω22Ω32Ω13Ω23Ω33⎦⎤Δt
然后：
B
=
B
−
∑
β
[
μ
1
(
β
)
μ
4
(
β
)
+
ω
1
(
β
)
μ
6
(
β
)
+
ω
3
(
β
)
μ
4
(
β
)
−
ω
1
(
β
)
μ
2
(
β
)
μ
5
(
β
)
+
ω
2
(
β
)
μ
6
(
β
)
−
ω
3
(
β
)
μ
5
(
β
)
−
ω
2
(
β
)
μ
3
(
β
)
]
Δ
γ
(
β
)
B=B-\sum_\beta \left[ \begin{matrix} \mu_{1}^{(\beta)} & \mu_{4}^{(\beta)} +\omega_1^{(\beta)} & \mu_{6} ^{(\beta)} +\omega_3^{(\beta)} \\ \mu_{4}^{(\beta)} -\omega_1^{(\beta)} & \mu_{2}^{(\beta)} & \mu_{5}^{(\beta)} +\omega_2^{(\beta)} \\ \mu_{6}^{(\beta)} -\omega_3^{(\beta)} & \mu_{5}^{(\beta)} -\omega_2^{(\beta)} & \mu_{3}^{(\beta)} \end{matrix} \right] \Delta \gamma^{(\beta)}
B=B−β∑⎣⎢⎡μ1(β)μ4(β)−ω1(β)μ6(β)−ω3(β)μ4(β)+ω1(β)μ2(β)μ5(β)−ω2(β)μ6(β)+ω3(β)μ5(β)+ω2(β)μ3(β)⎦⎥⎤Δγ(β)
Do I=1,3
B(I,J)=strainIncD(I)+DSPIN(I,I)
End Do
B(1,2)=strainIncD(4)+DSPIN(1,2)
B(2,3)=strainIncD(5)+DSPIN(2,3)
B(3,1)=strainIncD(6)+DSPIN(3,1)
B(2,1)=strainIncD(4)+DSPIN(2,1)
B(3,2)=strainIncD(5)+DSPIN(3,2)
B(1,3)=strainIncD(6)+DSPIN(1,3)
Do K=1,NSLPTL
Do I=1,3
B(I,J)=B(I,J)-Dgamma(K)*mu(K,1)
End Do
B(1,2)=B(1,2)-Dgamma(K)*(mu(K,4)+omega(K,1))
B(2,3)=B(2,3)-Dgamma(K)*(mu(K,5)+omega(K,2))
B(3,1)=B(3,1)-Dgamma(K)*(mu(K,6)-omega(K,3))
B(2,1)=B(2,1)-Dgamma(K)*(mu(K,4)-omega(K,1))
B(3,2)=B(3,2)-Dgamma(K)*(mu(K,5)-omega(K,2))
B(1,3)=B(1,3)-Dgamma(K)*(mu(K,6)+omega(K,3))
End Do
Do K=1,NSLPTL
Do I=1,3
Term1(I)=StateOld(3*(NSLPTL-1+K)+I)
End Do
Term2=matmul(Term1,B)
Do I=1,3
StateNew(3*(NSLPTL-1+K)+I)=StateOld(3*(NSLPTL-1+K)+I)-Term2(I)
End Do
Do I=1,3
Term1(I)=StateOld(6*(NSLPTL-1+K)+I)
End Do
Term2=matmul(B,Term1)
Do I=1,3
StateNew(6*(NSLPTL-1+K)+I)=StateOld(6*(NSLPTL-1+K)+I)+Term2(I)
End Do
End Do
更新每个滑移系的累积滑移量
∑
∣
Δ
γ
(
α
)
∣
\sum\begin{vmatrix}\Delta\gamma^{(\alpha)}\end{vmatrix}
∑∣∣Δγ(α)∣∣ STATE(9*NSLPTL+1)~STATE(10*NSLPTL)和更新Taylor累积滑移量
∑
∣
Δ
γ
∣
\sum\begin{vmatrix}\Delta\gamma\end{vmatrix}
∑∣∣Δγ∣∣ STATE(10*NSLPTL+1)
StateNew(10*NSLPTL+1)=0.0
Do I=1,NSLPTL
StateNew(9*NSLPTL+I)=StateOld(9*NSLPTL+I)+Abs(Dgamma(I))
StateNew(10*NSLPTL+1)=StateNew(10*NSLPTL+1)+Abs(Dgamma(I))
End Do
StateNew(10*NSLPTL+1)=StateOld(10*NSLPTL+1)+StateNew(10*NSLPTL+1)
8.材料参数
| 材料参数 | 序号 |
|---|---|
| 弹性模量 E E E | props(1) |
| 泊松比 μ \mu μ | props(2) |
| 参考滑移率 a ˙ \dot a a˙ | props(3) |
| 滑移率硬化指数 n n n | props(4) |
| 参考滑移硬化率 h 0 h_0 h0 | props(5) |
| 屈服剪切应力 τ 0 \tau_0 τ0 | props(6) |
| 突破剪切应力 τ s \tau_s τs | props(7) |
| 潜硬化系数 q q q | props(8) |
| 向前插值系数 θ \theta θ | props(9) |
| 初始欧拉角 φ 1 \varphi_1 φ1 | props(10) |
| 初始欧拉角 Φ \Phi Φ | props(11) |
| 初始欧拉角 φ 2 \varphi_2 φ2 | props(12) |
求出各个滑移系的滑移量
Δ
γ
(
α
)
\Delta \gamma^{(\alpha)}
Δγ(α) Dgamma，计算出塑性应变增量
Δ
ε
P
\Delta \varepsilon^P
ΔεP strainIncP(6)和塑性变形梯度增量
L
P
Δ
t
L^P \Delta t
LPΔtdefgradP (9)
Δ ε P = ∑ β μ i j ( β ) Δ γ ( β ) \Delta \varepsilon^P=\sum_\beta\mu_{ij}^{(\beta)}\Delta\gamma^{(\beta)} ΔεP=β∑μij(β)Δγ(β)
C****计算塑性应变增量
Do I=1,6
strainInP(I)=0.0
Do J=1,NSLPTL
strainInP(I)=strainInP(I)+mu(J,I)*Dgamma(J)
End Do
End Do
L P Δ t = ∑ β [ μ i j ( β ) + ω i j ( β ) ] Δ γ ( β ) L^P\Delta t=\sum_\beta \left[ \mu_{ij}^{(\beta)} + \omega_{ij}^{(\beta)} \right] \Delta \gamma^{(\beta)} LPΔt=β∑[μij(β)+ωij(β)]Δγ(β)
C****计算变形梯度增量
Do I=1,3
defgradP(I)=strainInP(I)
End Do
defgradP(4)=strainInP(4)
defgradP(7)=defgradP(7)
Do I=1,NSLPTL
defgradP(4)=defgradP(4)+omega(I,1)*Dgamma(I)
defgradP(7)=defgradP(7)-omega(I,1)*Dgamma(I)
End Do
defgradP(5)=strainInP(5)
defgradP(8)=defgradP(8)
Do I=1,NSLPTL
defgradP(5)=defgradP(5)+omega(I,2)*Dgamma(I)
defgradP(8)=defgradP(8)-omega(I,2)*Dgamma(I)
End Do
defgradP(6)=strainInP(6)
defgradP(9)=defgradP(9)
Do I=1,NSLPTL
defgradP(9)=defgradP(9)+omega(I,3)*Dgamma(I)
defgradP(6)=defgradP(6)-omega(I,3)*Dgamma(I)
End Do