FCC单晶塑性abaqus-vumat

• 0. 整体思路
• 1. 弹性各向同性
• 2. 12个滑移系
• 3. 旋转
• 4. 率相关的滑移硬化本构
• 5. 求解12个滑移量
• 6. 应力更新
• 7. 更新状态变量
• 8.材料参数

以黄永刚老师的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​+ΔσV​i=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​=⎣⎢⎢⎢⎢⎢⎢⎡​2G00000​02G0000​002G000​0002G00​00002G0​000002G​⎦⎥⎥⎥⎥⎥⎥⎤​

自定义材料参数：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φ2​0​sinφ2​cosφ2​0​00​1​⎦⎤​⎣⎡​100​0cosΦ−sinΦ​0sinΦcosΦ​⎦⎤​⎣⎡​cosφ1​−sinφ1​0​sinφ1​cosφ1​0​00​1​⎦⎤​(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φ1​cosφ2​−sinφ1​sinφ2​cosΦ−cosφ1​sinφ2​−sinφ1​cosφ2​cosΦsinφ1​sinΦ​sinφ1​cosφ2​+cosφ1​sinφ2​cosΦ−sinφ1​sinφ2​+cosφ1​cosφ2​cosΦ−cosφ1​sinΦ​sinφ2​sinΦcosφ2​sinΦ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 ​1​3 ​1​3 ​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 ​1​2 ​−1​0>⋅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(α)​=⎣⎡​s1​m1​s2​m1​s3​m1​​s1​m2​s2​m2​s3​m2​​s1​m3​s2​m3​s3​m3​​⎦⎤​=⎣⎡​uhvhwh​ukvkwk​ulvlwl​⎦⎤​

μ 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(α)​=⎣⎡​s1​m1​2s1​m2​+s2​m1​​2s1​m3​+s3​m1​​​2s1​m2​+s2​m1​​s2​m2​2s2​m3​+s3​m2​​​2s1​m3​+s3​m1​​2s2​m3​+s3​m2​​s3​m3​​⎦⎤​

ω 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(α)​=⎣⎡​02s2​m1​−s1​m2​​2s3​m1​−s1​m3​​​2s1​m2​−s2​m1​​02s3​m2​−s2​m3​​​2s1​m3​−s3​m1​​2s2​m3​−s3​m2​​0​⎦⎤​

展开：
μ 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(α)​=⎣⎢⎢⎢⎢⎢⎢⎡​s1​m1​s2​m2​s3​m3​21​(s1​m2​+s2​m1​)21​(s2​m3​+s3​m2​)21​(s1​m3​+s3​m1​)​⎦⎥⎥⎥⎥⎥⎥⎤​ωI(α)​=⎣⎡​21​(s1​m2​−s2​m1​)21​(s2​m3​−s3​m2​)21​(s1​m3​−s3​m1​)​⎦⎤​

通过状态变量 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=⎣⎡​F1​F7​F6​​F4​F2​F8​​F9​F5​F3​​⎦⎤​U=⎣⎡​U1​U4​U6​​U4​U2​U5​​U6​U5​U3​​⎦⎤​
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)=U1​U2​U3​+2U4​U5​U6​−U1​U5​U5​−U3​U4​U4​−U2​U6​U6​
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∗=⎣⎡​U2​U3​−U5​U5​U5​U6​−U3​U4​U4​U5​−U2​U6​​U5​U6​−U3​U4​U1​U3​−U6​U6​U4​U6​−U1​U5​​U4​U5​−U2​U6​U4​U6​−U1​U5​U1​U2​−U4​U4​​⎦⎤​

		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(γ)=h0​sech2∣∣∣∣∣​τs​−τ0​h0​γ​∣∣∣∣∣​(α不求和)(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 τs​props(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​​ω1​0−ω2​​ω3​ω2​0​⎦⎤​⎣⎡​σ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μ4​2(−ω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μ5​2(−ω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​ΔεV​I=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
     

标签：abaqus,NB,Delta,stretchOld,FCC,stretchNew,alpha,vumat,NSLPTL
来源： https://blog.csdn.net/Niklas4088/article/details/112713431