FCC单晶塑性abaqus-vumat

以黄永刚老师的umat文章为基础，构建单晶塑性的vumat。本文是之前“单晶塑性abaqus-umat 自学”的后续，大量引用前文已推导出的结论，用（umat）代指上文。求解过程针对FCC。

0. 整体思路

先从最简单的入手：1) 弹性各向同性；2)选用向前梯度插值，先不要Newton迭代法求解；3) 滑移只考虑Schmid准则；4) 求解方程不是通过LUP分解和反向代入法，而是采用求逆矩阵的方式（逆矩阵求解方法见上一篇文章）。求逆矩阵的方法写起来相对简洁，但计算量大。

---

1. 弹性各向同性

弹性采用各向同性，且分析过程中将体积变化和畸变分开考虑，动态分析中可根据状态方程(EOS)单独给出P-V关系。

体变： $$\Delta {\sigma }_V=K\Delta {\epsilon }_v\text{\hspace{1em}(1.1a)}$$

等待偏应力更新之后，再在主应力上叠加体应力
$${\sigma }_i={\sigma }_i+\Delta {\sigma }_{V}i=1,2,3\text{\hspace{1em}(1.1b)}$$

客观应力率：(后续推导中不特意说明就是偏应力和偏应变)
$${\stackrel{\nabla }{\sigma }}^{\ast }=\mathbf{L}:{\mathbf{D}}^{\ast }\text{\hspace{1em}(1.2a)}$$

$$\dot{\mathbf{\sigma }}-{\mathbf{\Omega }}^{\ast }\cdot \mathbf{\sigma }+\mathbf{\sigma }\cdot {\mathbf{\Omega }}^{\ast }=\mathbf{L}:{\mathbf{D}}^{\ast }\text{\hspace{1em}(1.2b)}$$

方程两边同乘以$\Delta t$
$$\mathbf{\Delta }\mathbf{\sigma }-{\mathbf{\Omega }}^{\ast }\cdot \mathbf{\sigma }\mathbf{\Delta }\mathbf{t}+\mathbf{\sigma }\cdot {\mathbf{\Omega }}^{\ast }\mathbf{\Delta }\mathbf{t}=\mathbf{L}:\mathbf{\Delta }{\mathbf{\epsilon }}^{\ast }\text{\hspace{1em}(1.2c)}$$

由于各向同性假设，2个弹性系数 $E\mu$，且vumat给出的应变增量的剪切分量不是工程应变，因此相应的弹性矩阵为（偏应变与偏应力）：$$L_{ijkl}=2G{\delta }_{ik}{\delta }_{jl}=\left[\begin{array}{cccccc}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{array}\right]$$

自定义材料参数：props(1)=杨氏模量 $E$，props(2)=泊松比 $\mu$.

---

2. 12个滑移系

滑移面 { 111 } \{111\} {111}：滑移方向$<1\bar{1}0><10\bar{1}><01\bar{1}>$
滑移面 { 1 1 ˉ 1 } \{1\bar11\} {11ˉ1}：滑移方向$<110><10\bar{1}><011>$
滑移面 { 11 1 ˉ } \{11\bar1\} {111ˉ}：滑移方向$<1\bar{1}0><101><011>$
滑移面 { 1 1 ˉ 1 ˉ } \{1\bar1\bar1\} {11ˉ1ˉ}：滑移方向$<110><101><01\bar{1}>$
初始方向：三个欧拉角 ${\phi }_1$，$\Phi$，${\phi }_2$，通过props(10) props(11) props(12)赋值
经过欧拉角转动后的晶体取向 $g$ 为：
$$g=\left[\begin{array}{cccc}\cos {\phi }_2& \sin {\phi }_2& 0& \\ -\sin {\phi }_2& \cos {\phi }_2& 0& \\ 0& 0& & 1\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \Phi & \sin \Phi \\ 0& -\sin \Phi & \cos \Phi \end{array}\right]\left[\begin{array}{cccc}\cos {\phi }_1& \sin {\phi }_1& 0& \\ -\sin {\phi }_1& \cos {\phi }_1& 0& \\ 0& 0& & 1\end{array}\right]\text{\hspace{1em}(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=<uvw>=<\frac{1}{\sqrt{2}}\frac{-1}{\sqrt{2}}0>\cdot g\text{\hspace{1em}(2.2b)}$$
同理计算其他滑移系，与此同时对滑移系进行归一化。

计算出12个滑移系的Schmid因子”${\mu }_{ij}^{(\alpha )}$和张量${\omega }_{ij}^{(\alpha )}$:

$${\mu }_{ij}^{(\alpha )}=\frac{1}{2}\left[s_i^{(\alpha )}{m}_j^{(\alpha )}+s_j^{(\alpha )}{m}_i^{(\alpha )}\right]\text{\hspace{1em}(2.3a)}$$

$${\omega }_{ij}^{(\alpha )}=\frac{1}{2}\left[s_i^{(\alpha )}{m}_j^{(\alpha )}-s_j^{(\alpha )}{m}_i^{(\alpha )}\right]\text{\hspace{1em}(2.3b)}$$

其中:$$s_i^{(\alpha )}{m}_j^{(\alpha )}=\left[\begin{array}{ccc}{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{array}\right]=\left[\begin{array}{ccc}uh& uk& ul\\ vh& vk& vl\\ wh& wk& wl\end{array}\right]$$

$${\mu }_{ij}^{(\alpha )}=\left[\begin{array}{ccc}{s}_1{m}_1& \frac{s_1{m}_2+s_2{m}_1}{2}& \frac{s_1{m}_3+s_3{m}_1}{2}\\ \frac{s_1{m}_2+s_2{m}_1}{2}& s_2{m}_2& \frac{s_2{m}_3+s_3{m}_2}{2}\\ \frac{s_1{m}_3+s_3{m}_1}{2}& \frac{s_2{m}_3+s_3{m}_2}{2}& s_3{m}_3\end{array}\right]$$

$${\omega }_{ij}^{(\alpha )}=\left[\begin{array}{ccc}0& \frac{s_1{m}_2-s_2{m}_1}{2}& \frac{s_1{m}_3-s_3{m}_1}{2}\\ \frac{s_2{m}_1-s_1{m}_2}{2}& 0& \frac{s_2{m}_3-s_3{m}_2}{2}\\ \frac{s_3{m}_1-s_1{m}_3}{2}& \frac{s_3{m}_2-s_2{m}_3}{2}& 0\end{array}\right]$$

展开：
$${\mu }_I^{(\alpha )}=\left[\begin{array}{c}{s}_1{m}_1\\ s_2{m}_2\\ s_3{m}_3\\ \frac{1}{2}(s_1{m}_2+s_2{m}_1)\\ \frac{1}{2}(s_2{m}_3+s_3{m}_2)\\ \frac{1}{2}(s_1{m}_3+s_3{m}_1)\end{array}\right]{\omega }_I^{(\alpha )}=\left[\begin{array}{c}\frac{1}{2}(s_1{m}_2-s_2{m}_1)\\ \frac{1}{2}(s_2{m}_3-s_3{m}_2)\\ \frac{1}{2}(s_1{m}_3-s_3{m}_1)\end{array}\right]$$

通过状态变量 STATEV(3*NSLPTL+1)~STATEV(6*NSLPTL)记录滑移面法向 $m$ ，STATEV(6*NSLPTL+1)~STATEV(9*NSLPTL)记录相对应的滑移方向 $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中给出了$\Delta R$ DROT, 需要通过$F=RU\to F{U}^{-1}=R$求出新老旧两个$R$，通过做差求出$\Delta R$：
$$\Delta R=R^{New}-R^{Old}\text{\hspace{1em}(3.1)}$$
$$F=\left[\begin{array}{ccc}{F}_1& F_4& F_9\\ F_7& F_2& F_5\\ F_6& F_8& F_3\end{array}\right]U=\left[\begin{array}{ccc}{U}_1& U_4& U_6\\ U_4& U_2& U_5\\ U_6& U_5& U_3\end{array}\right]$$
$$\det (U)=U_1{U}_2{U}_3+2U_4{U}_5{U}_6-U_1{U}_5{U}_5-U_3{U}_4{U}_4-U_2{U}_6{U}_6$$
$U$ 的伴随 $U^{\ast }$
$$U^{\ast }=\left[\begin{array}{ccc}{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\end{array}\right]$$

		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,:)

再通过$\Delta R$DROT，求出 $\Omega \Delta t$DSPIN：
$$\Delta R={\left(I-\frac{1}{2}\Omega \Delta t\right)}^{-1}\cdot \left(I+\frac{1}{2}\Omega \Delta t\right)\text{\hspace{1em}(3.2a)}$$
$$\Omega \Delta t=2\left(\Delta R-I\right){\left(\Delta R+I\right)}^{-1}\text{\hspace{1em}(3.2b)}$$
$$\Delta R=\left[\begin{array}{ccc}\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{array}\right]$$
$$\Delta R+I=\left[\begin{array}{ccc}\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{array}\right]$$
$$(\Delta R+I{)}^{}=\left[\begin{array}{ccc}\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{array}\right]$$

---

4. 率相关的滑移硬化本构

基于Schmid准则，$\alpha$ 滑移系的滑移率 ${\dot{\gamma }}^{(\alpha )}$ 由分切应力${\tau }^{(\alpha )}$决定：
$${\dot{\gamma }}^{(\alpha )}=\dot{a}{f}^{(\alpha )}\left(\frac{{\tau }^{(\alpha )}}{g^{(\alpha )}}\right)\text{\hspace{1em}(4.1a)}$$

其中，$$f^{(\alpha )}(x)=x|x{|}^{n-1}\text{\hspace{1em}(4.1b)}$$

应变硬化： $${\dot{g}}^{(\alpha )}=\sum _{\beta }{h}_{\alpha \beta }{\dot{\gamma }}^{(\beta )}$$

其中自硬化系数： $$h_{\alpha \alpha }=h(\gamma )=h_0{\mathrm{sech}}^2\left|\frac{h_0\gamma }{{\tau }_s-{\tau }_0}\right|(\alpha 不求和)\text{\hspace{1em}(4.2a)}$$

其中 $h_0$ 为初始硬化模量，${\tau }_0$屈服强度，等于当前强度的初始值 $g^{(\alpha )}(0)$，$\gamma$ 是所有滑移系上的Taylor累计剪切应变：
$$\gamma =\sum _{\alpha }{\int }_0^t\left|{\dot{\gamma }}^{(\alpha )}\right|dt\text{\hspace{1em}(4.2b)}$$
潜硬化模量：
$$h_{\alpha \beta }=qh(\gamma )(\alpha \ne \beta )\text{\hspace{1em}(4.2c)}$$

自定义材料参数：props(3)=参考滑移率 $\dot{a}$，props(4)=滑移率指数 $n$，props(5)=参考硬化系数$h_0$，props(6)=初始屈服强度${\tau }_0$，props(7)=break-through突破强度${\tau }_s$props(8)=潜硬化系数 $q$。

---

5. 求解12个滑移量

向前梯度插值求解黄老师umat中的公式（3.2.5）：
$$\begin{array}{rl}& \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 }\mathrm{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 {\epsilon }_{ij}\\ & \\ & 转化为⟹\mathbf{A}\mathbf{A}\cdot \mathbf{\Delta }\mathbf{\gamma }=\mathbf{Y}\end{array}\text{\hspace{1em}(5.1)}$$
props(9)=梯度差值系数 $\theta$

求解 $\frac{\partial {\dot{\gamma }}^{(\alpha )}}{\partial {\tau }^{(\alpha )}}$ 和 $\frac{\partial {\dot{\gamma }}^{(\alpha )}}{\partial g^{(\alpha )}}$，对（4.1） 求导：
$$\frac{\partial {\dot{\gamma }}^{(\alpha )}}{\partial f^{(\alpha )}}=\dot{a}\frac{d{f}^{(\alpha )}}{dx}=n|x{|}^{n-1}其中x=\frac{{\tau }^{(\alpha )}}{g^{(\alpha )}}\text{\hspace{1em}(5.2a)}$$

$$\frac{\partial x}{\partial {\tau }^{(\alpha )}}=\frac{1}{g^{(\alpha )}}\frac{\partial x}{\partial g^{(\alpha )}}=\frac{-{\tau }^{(\alpha )}}{g^{(\alpha )2}}\text{\hspace{1em}(5.2b)}$$

$$\frac{\partial {\dot{\gamma }}^{(\alpha )}}{\partial {\tau }^{(\alpha )}}=n\dot{a}{\left|\frac{{\tau }^{(\alpha )}}{g^{(\alpha )}}\right|}^{n-1}\frac{1}{g^{(\alpha )}}\text{\hspace{1em}(5.2c)}$$

$$\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}}\text{\hspace{1em}(5.2d)}$$

设 $A_{ij}^{(\alpha )}=L_{ijkl}{\mu }_{kl}^{(\alpha )}+{\omega }_{ik}^{(\alpha )}{\sigma }_{jk}+{\omega }_{jk}^{(\alpha )}{\sigma }_{ik}$：
$$L_{ijkl}{\mu }_{kl}^{(\alpha )}=2G{\mu }_{ij}^{(\alpha )}$$

$$\begin{array}{rl}{\omega }_{ik}^{(\alpha )}{\sigma }_{jk}& =\left[\begin{array}{ccc}0& {\omega }_1& {\omega }_3\\ -{\omega }_1& 0& {\omega }_2\\ -{\omega }_3& -{\omega }_2& 0\end{array}\right]\left[\begin{array}{ccc}{\sigma }_1& {\sigma }_4& {\sigma }_6\\ {\sigma }_4& {\sigma }_2& {\sigma }_5\\ {\sigma }_6& {\sigma }_5& {\sigma }_3\end{array}\right]\\ & =\left[\begin{array}{ccc}{\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{array}\right]\end{array}$$

$${\omega }_{jk}^{(\alpha )}{\sigma }_{ik}={\sigma }_{ik}{\omega }_{jk}^{(\alpha )}={\sigma }_{ik}{\omega }_{kj}^{(\alpha )T}=({\omega }_{ik}^{(\alpha )}{\sigma }_{kj}{)}^T$$

$$\begin{array}{rl}& A_{ij}^{(\alpha )}=L_{ijkl}{\mu }_{kl}^{(\alpha )}+{\omega }_{ik}^{(\alpha )}{\sigma }_{jk}+{\omega }_{jk}^{(\alpha )}{\sigma }_{ik}\\ & =\left[\begin{array}{ccc}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{array}\right]\end{array}$$

滑移系的当前强度 $g^{(\alpha )}$ 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）：
$$\Delta {\sigma }_{ij}=L_{ijkl}\Delta {\epsilon }_{kl}-{\sigma }_{ij}\Delta {\epsilon }_{kk}-\sum _{\alpha }\left[L_{ijkl}{\mu }_{kl}^{(\alpha )}+{\omega }_{ik}^{(\alpha )}{\sigma }_{jk}+{\omega }_{jk}^{(\alpha )}{\sigma }_{ik}\right]\cdot \Delta {\gamma }^{(\alpha )}$$
改写：
$$\Delta {\sigma }_I=2G\Delta {\epsilon }_I-\sum _{\alpha }{A}_I^{(\alpha )}\cdot \Delta {\gamma }^{(\alpha )}$$
叠加体应力：
$$\Delta {\sigma }_I=\Delta {\sigma }_I+\frac{E}{3(1-2\mu )}\Delta {\epsilon }_{V}I=1,\mathrm{N}\mathrm{D}\mathrm{I}$$

      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^{(\alpha )}$	STATE(1)~STATE(NSLPTL)
滑移率 ${\dot{\gamma }}^{\alpha }$	STATE(NSLPTL+1)~STATE(2*NSLPTL)
当前的分切应力 ${\tau }^{(\alpha )}$	STATE(2*NSLPTL+1)~STATE(3*NSLPTL)
滑移面法向 $m$	STATE(3*NSLPTL+1)~STATE(6*NSLPTL)
滑移方向 $s$	STATE(6*NSLPTL+1)~STATE(9*NSLPTL)
累积滑移量 $\sum \left|\begin{array}{c}\Delta {\gamma }^{(\alpha )}\end{array}\right|$	STATE(9*NSLPTL+1)~STATE(10*NSLPTL)
Taylor累积滑移量 $\sum \left|\begin{array}{c}\Delta \gamma \end{array}\right|$	STATE(10*NSLPTL+1)

滑移系的当前强度 $g^{(\alpha )}$ 更新：
$$\Delta g^{(\alpha )}=\sum _{\beta }{h}_{\alpha \beta }\Delta {\gamma }^{(\beta )}$$

      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 }$ 更新 ${\dot{\gamma }}^{\alpha }=\frac{\Delta {\gamma }^{(\alpha )}}{\Delta t}$:

      Do I=1,NSLPTL
      	StaveNew(NSLPTL+I)=Dgamma(I)/dtArray(1)
      End Do

分切应力 ${\tau }^{(\alpha )}$ 更新：
$$\Delta {\tau }^{(\alpha )}=\left[L_{ijkl}{\mu }_{kl}^{(\alpha )}+{\omega }_{ik}^{(\alpha )}{\sigma }_{jk}+{\omega }_{jk}^{(\alpha )}{\sigma }_{ik}\right]:\left[\Delta {\epsilon }_{ij}-\sum _{\beta }{\mu }_{ij}^{(\beta )}\Delta {\gamma }^{(\beta )}\right]$$
改写：
$$\Delta {\tau }^{(\alpha )}=A_I^{(\alpha )}\cdot \left[\Delta {\epsilon }_I-\sum _{\beta }{\mu }_I^{(\beta )}\Delta {\gamma }^{(\beta )}\right]$$

      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

更新滑移系：
$$\Delta m_i^{\ast (\alpha )}=-m_j^{\ast (\alpha )}\left\{\Delta {\epsilon }_{ji}+{\Omega }_{ji}\Delta t-\sum _{\beta }\left[{\mu }_{ji}^{(\beta )}+{\omega }_{ji}^{(\beta )}\right]\Delta {\gamma }^{(\beta )}\right\}$$

$$\Delta s_i^{\ast (\alpha )}=\left\{\Delta {\epsilon }_{ij}+{\Omega }_{ij}\Delta t-\sum _{\beta }\left[{\mu }_{ij}^{(\beta )}+{\omega }_{ij}^{(\beta )}\right]\Delta {\gamma }^{(\beta )}\right\}{s}_j^{\ast (\alpha )}$$

设：$B(3,3)=\Delta {\epsilon }_{ij}+{\Omega }_{ij}\Delta t-{\sum }_{\beta }\left[{\mu }_{ij}^{(\beta )}+{\omega }_{ij}^{(\beta )}\right]\Delta {\gamma }^{(\beta )}$
首先：
$$B=\left[\begin{array}{ccc}\Delta {\epsilon }_1& \Delta {\epsilon }_4& \Delta {\epsilon }_6\\ \Delta {\epsilon }_4& \Delta {\epsilon }_2& \Delta {\epsilon }_5\\ \Delta {\epsilon }_6& \Delta {\epsilon }_5& \Delta {\epsilon }_3\end{array}\right]+\left[\begin{array}{ccc}{\Omega }_{11}& {\Omega }_{12}& {\Omega }_{13}\\ {\Omega }_{21}& {\Omega }_{22}& {\Omega }_{23}\\ {\Omega }_{31}& {\Omega }_{32}& {\Omega }_{33}\end{array}\right]\Delta t$$
然后：
$$B=B-\sum _{\beta }\left[\begin{array}{ccc}{\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{array}\right]\Delta {\gamma }^{(\beta )}$$

      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 \left|\begin{array}{c}\Delta {\gamma }^{(\alpha )}\end{array}\right|$ STATE(9*NSLPTL+1)~STATE(10*NSLPTL)和更新Taylor累积滑移量 $\sum \left|\begin{array}{c}\Delta \gamma \end{array}\right|$ 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$	props(1)
泊松比 $\mu$	props(2)
参考滑移率 $\dot{a}$	props(3)
滑移率硬化指数$n$	props(4)
参考滑移硬化率$h_0$	props(5)
屈服剪切应力 ${\tau }_0$	props(6)
突破剪切应力 ${\tau }_s$	props(7)
潜硬化系数 $q$	props(8)
向前插值系数 $\theta$	props(9)
初始欧拉角 ${\phi }_1$	props(10)
初始欧拉角 $\Phi$	props(11)
初始欧拉角 ${\phi }_2$	props(12)

求出各个滑移系的滑移量$\Delta {\gamma }^{(\alpha )}$ Dgamma，计算出塑性应变增量 $\Delta {\epsilon }^P$ strainIncP(6)和塑性变形梯度增量 $L^P\Delta t$defgradP (9)

$$\Delta {\epsilon }^P=\sum _{\beta }{\mu }_{ij}^{(\beta )}\Delta {\gamma }^{(\beta )}$$

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\Delta t=\sum _{\beta }\left[{\mu }_{ij}^{(\beta )}+{\omega }_{ij}^{(\beta )}\right]\Delta {\gamma }^{(\beta )}$$

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