前言

繼續塑性本構，隨動硬化塑性。
記錄自己學習abaqus軟體umat子程式的t過程，本文主要參考了《非線性本構關係在ABAQUS中的實現》第四章和技術鄰的視訊課程“非線性各向同性強化彈塑性umat子程式教程”

塑性力學增量形式實現

計算過程中，體應變和體應力是彈性關係$$\Delta {\sigma }_v=K\cdot \Delta {\epsilon }_v\text{\hspace{1em}(1)}$$

$K$為體積模量，$K=\frac{E}{3(1-\nu )}$ , $\Delta {\epsilon }_v=\Delta {\epsilon }_{11}+\Delta {\epsilon }_{22}+\Delta {\epsilon }_{33}$

所以在下面的討論中只考慮偏應力張量和偏應力張量

試應力： $${\mathbf{\sigma }}^{\mathbf{t}{\mathbf{r}}^{\prime }}(\mathbf{t})={\mathbf{\sigma }}^{\prime }(\mathbf{t})+{\mathbb{C}}^{\prime }:\mathbf{\Delta }{\mathbf{\epsilon }}^{\prime }\text{\hspace{1em}(2)}$$

線性隨動硬化塑性的屈服函式:$$f={\left[\frac{3}{2}({\mathbf{\sigma }}^{\prime }-{\mathbf{\alpha }}^{\prime }):({\mathbf{\sigma }}^{\prime }-{\mathbf{\alpha }}^{\prime })\right]}^{\frac{1}{2}}-{\sigma }_y={\hat{\sigma }}_{Mises}-{\sigma }_y\text{\hspace{1em}(3)}$$

其中，${\mathbf{\alpha }}^{\prime }$為背應力，三維問題6個分量；${\sigma }_y$初始屈服強度，${\hat{\sigma }}_{Mises}$ 有效應力的Mises等效應力。

線性隨動硬化模型：Mises等效背應力與等效塑性增量成正比
$$\begin{gathered} \Delta {\alpha }_{Mises}=c\Delta {\bar{\epsilon }}^p \\ \mathbf{\Delta }\mathbf{\alpha }=\frac{2}{3}c\Delta {\mathbf{\epsilon }}^{\mathbf{p}}\text{\hspace{1em}(4)} \end{gathered}$$
$c$隨動硬化係數。

${\mathbf{\alpha }}^{\prime }$背應力作為狀態變數存在STATEV(NSTATV)中，提取${\mathbf{\alpha }}^{\prime }$背應力並計算試應力的有效屈服應力 ${\sigma }_{Mises}^{tr}$

判斷 ${\sigma }_{Mises}^{tr}$ 與 ${\sigma }_y$ 關係:

若 ${\sigma }_{Mises}^{tr}＜{\sigma }_y$ ： $\mathbf{\sigma }(\mathbf{t}+\mathbf{\Delta }\mathbf{t})={\mathbf{\sigma }}^{\mathbf{t}\mathbf{r}}(\mathbf{t})$

一致切線剛度矩陣，為彈性剛度矩陣
$$\mathbf{D}\mathbf{D}\mathbf{S}\mathbf{D}\mathbf{D}\mathbf{E}(i,j)=\left[\begin{array}{cccccc}2G+\lambda & \lambda & \lambda & 0& 0& 0\\ \lambda & 2G+\lambda & \lambda & 0& 0& 0\\ \lambda & \lambda & 2G+\lambda & 0& 0& 0\\ 0& 0& 0& G& 0& 0\\ 0& 0& 0& 0& G& 0\\ 0& 0& 0& 0& 0& G\end{array}\right]\text{\hspace{1em}(5)}$$

若 ${\sigma }_{Mises}^{tr}\ge {\sigma }_y$ ，計算新增的 $\Delta {\bar{\epsilon }}^p$
$$\Delta {\bar{\epsilon }}^p=\frac{{\sigma }_{Mises}^{tr}-{\sigma }_y}{(c+3G)}$$
根據正交流動法則：
$$\mathbf{\Delta }{\mathbf{\epsilon }}^{\mathbf{p}}=\Delta {\bar{\epsilon }}^p\frac{3}{2}\frac{{\hat{\mathbf{\sigma }}}^{\prime }}{{\hat{\sigma }}_{Mises}}=\Delta {\bar{\epsilon }}^p\frac{3}{2}\frac{{\hat{\mathbf{\sigma }}}^{\mathbf{t}{\mathbf{r}}^{\prime }}}{{\hat{\sigma }}_{Mises}^{tr}}$$
代入(4)更新背應力，STATEV(1)……STATEV(NTENS)

應力增量=彈性剛度矩陣×彈性應變增量
$$\mathbf{\Delta }\mathbf{\sigma }={\mathbf{D}}_{\mathbf{e}}：\mathbf{\Delta }{\mathbf{\epsilon }}^{\mathbf{e}}$$
彈性應變=總應變-塑性應變增量
$$\mathbf{\Delta }\mathbf{\sigma }={\mathbf{D}}_{\mathbf{e}}：(\mathbf{\Delta }\mathbf{\epsilon }-\mathbf{\Delta }{\mathbf{\epsilon }}^{\mathbf{p}})\text{\hspace{1em}(6)}$$

屈服後一致切線剛度矩陣：

將 $\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}}$ 用 $\mathbf{d}\mathbf{\epsilon }$ 表示出來，再代回式(6)就可求出屈服後的DDSDDE(NTENS,NTENS)

一致性條件
$$\{\begin{array}{l}f({\mathbf{\sigma }}^{\prime },\mathbf{\alpha })=0\\ f({\mathbf{\sigma }}^{\prime }+d{\mathbf{\sigma }}^{\prime },\mathbf{\alpha }+d\mathbf{\alpha })=0\end{array}⟹\frac{\partial f}{\partial {\mathbf{\sigma }}^{\prime }}\mathbf{d}{\mathbf{\sigma }}^{\prime }+\frac{\partial f}{\partial \mathbf{\alpha }}\mathbf{d}\mathbf{\alpha }=0\text{\hspace{1em}(7)}$$
其中：$$\{\begin{array}{l}\frac{\partial f}{\partial {\mathbf{\sigma }}^{\prime }}=\frac{\partial {\left[\frac{3}{2}({\mathbf{\sigma }}^{\prime }-\mathbf{\alpha }):({\mathbf{\sigma }}^{\prime }-\mathbf{\alpha })\right]}^{\frac{1}{2}}}{\partial {\mathbf{\sigma }}^{\prime }}=\frac{6}{{\hat{\sigma }}_{Mises}}({\mathbf{\sigma }}^{\prime }-\mathbf{\alpha })\\ \\ \frac{\partial f}{\partial \mathbf{\alpha }}=-\frac{\partial f}{\partial {\mathbf{\sigma }}^{\prime }}\end{array}\text{\hspace{1em}(8)}$$
(4)、(6)和(8)代入(7)：
$$\begin{array}{rl}& \frac{\partial f}{\partial {\mathbf{\sigma }}^{\prime }}\mathbf{d}{\mathbf{\sigma }}^{\prime }=c\frac{\partial f}{\partial {\mathbf{\sigma }}^{\prime }}\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}}\\ & \\ & \mathbf{d}{\mathbf{\sigma }}^{\prime }=c\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}}\\ & \\ & {\mathbf{D}}_{\mathbf{e}}^{\prime }：(\mathbf{d}{\mathbf{\epsilon }}^{\prime }-\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}})=c\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}}\\ & \\ & {\mathbf{D}}_{\mathbf{e}}^{\prime }\cdot \mathbf{d}{\mathbf{\epsilon }}^{\prime }=c\mathbf{I}：\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}}+{\mathbf{D}}_{\mathbf{e}}^{\prime }：\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}}\\ & \\ & {\mathbf{D}}_{\mathbf{e}}^{\prime }：\mathbf{d}{\mathbf{\epsilon }}^{\prime }=(c\mathbf{I}+{\mathbf{D}}_{\mathbf{e}}^{\prime }):\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}}\\ & \\ & {\mathbf{D}}_{\mathbf{e}}^{\prime }:(\mathbf{I}-\frac{1}{3}\hat{\mathbf{I}}):\mathbf{d}\mathbf{\epsilon }=(c\mathbf{I}+{\mathbf{D}}_{\mathbf{e}}):\mathbf{d}{\mathbf{\epsilon }}^{\mathbf{p}}\end{array}\text{\hspace{1em}(9)}$$
其中：$\mathbf{I}={\delta }_{ik}{\delta }_{jl}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}{\vec{\mathbf{e}}}_{\mathbf{k}}{\vec{\mathbf{e}}}_{\mathbf{l}}$ 四階單位張量，$\hat{\mathbf{I}}={\delta }_{ij}{\delta }_{kl}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}{\vec{\mathbf{e}}}_{\mathbf{k}}{\vec{\mathbf{e}}}_{\mathbf{l}}$

（注：$\mathbf{I}:\mathbf{A}={\delta }_{ik}{\delta }_{jl}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}{\vec{\mathbf{e}}}_{\mathbf{k}}{\vec{\mathbf{e}}}_{\mathbf{l}}:A_{nm}{\vec{\mathbf{e}}}_{\mathbf{n}}{\vec{\mathbf{e}}}_{\mathbf{m}}={\delta }_{ik}{\delta }_{jl}{\delta }_{kn}{\delta }_{lm}{A}_{nm}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}={\delta }_{ik}{\delta }_{jl}{A}_{kl}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}=A_{ij}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}$）
（注：$\hat{\mathbf{I}}:\mathbf{A}={\delta }_{ij}{\delta }_{kl}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}{\vec{\mathbf{e}}}_{\mathbf{k}}{\vec{\mathbf{e}}}_{\mathbf{l}}:A_{nm}{\vec{\mathbf{e}}}_{\mathbf{n}}{\vec{\mathbf{e}}}_{\mathbf{m}}={\delta }_{ij}{\delta }_{kl}{\delta }_{kn}{\delta }_{lm}{A}_{nm}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}={\delta }_{ij}{A}_{kk}{\vec{\mathbf{e}}}_{\mathbf{i}}{\vec{\mathbf{e}}}_{\mathbf{j}}$）

再代入(6):
$$\begin{array}{rl}\mathbf{d}{\mathbf{\sigma }}^{\prime }& =\left\{{\mathbf{D}}_{\mathbf{e}}^{\prime }-{\mathbf{D}}_{\mathbf{e}}^{\prime }:(c\mathbf{I}+{\mathbf{D}}_{\mathbf{e}}^{\prime }{)}^{-1}:{\mathbf{D}}_{\mathbf{e}}^{\prime }\right\}:(\mathbf{I}-\frac{1}{3}\hat{\mathbf{I}}):\mathbf{d}\mathbf{\epsilon }\\ \mathbf{d}{\mathbf{\sigma }}_{\mathbf{V}}& =K\hat{\mathbf{I}}:\mathbf{d}\mathbf{\epsilon }\end{array}\text{\hspace{1em}(10)}$$
因此：$$\mathbf{d}\mathbf{\sigma }=\left\{\left[{\mathbf{D}}_{\mathbf{e}}^{\prime }-{\mathbf{D}}_{\mathbf{e}}^{\prime }:(c\mathbf{I}+{\mathbf{D}}_{\mathbf{e}}^{\prime }{)}^{-1}:{\mathbf{D}}_{\mathbf{e}}^{\prime }\right]:(\mathbf{I}-\frac{1}{3}\hat{\mathbf{I}})+K\hat{\mathbf{I}}\right\}:\mathbf{d}\mathbf{\epsilon }\text{\hspace{1em}(11)}$$

推導過程保持張量計法
$$\mathbf{A}=\left[\begin{array}{ccccccccc}{A}_{11,11}& A_{11,22}& A_{11,33}& A_{11,12}& A_{11,23}& A_{11,13}& A_{11,21}& A_{11,32}& A_{11,31}\\ A_{22,11}& A_{22,22}& A_{22,33}& A_{22,12}& A_{22,23}& A_{22,13}& A_{22,21}& A_{22,32}& A_{22,31}\\ A_{33,11}& A_{33,22}& A_{33,33}& A_{33,12}& A_{33,23}& A_{33,13}& A_{33,21}& A_{33,32}& A_{33,31}\\ A_{12,11}& A_{12,22}& A_{12,33}& A_{12,12}& A_{12,23}& A_{12,13}& A_{12,21}& A_{12,32}& A_{12,31}\\ A_{23,11}& A_{23,22}& A_{23,33}& A_{23,12}& A_{23,23}& A_{23,13}& A_{23,21}& A_{23,32}& A_{23,31}\\ A_{13,11}& A_{13,22}& A_{13,33}& A_{13,12}& A_{13,23}& A_{13,13}& A_{13,21}& A_{13,32}& A_{13,31}\\ A_{21,11}& A_{21,22}& A_{21,33}& A_{21,12}& A_{21,23}& A_{21,13}& A_{21,21}& A_{21,32}& A_{21,31}\\ A_{32,11}& A_{32,22}& A_{32,33}& A_{32,12}& A_{32,23}& A_{32,13}& A_{32,21}& A_{32,32}& A_{32,31}\\ A_{31,11}& A_{31,22}& A_{31,33}& A_{31,12}& A_{31,23}& A_{31,13}& A_{31,21}& A_{31,32}& A_{31,31}\end{array}\right]$$

$${\mathbf{D}}_{\mathbf{e}}^{\prime }=\left[\begin{array}{ccccccccc}2G& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 2G& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 2G& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& G& 0& 0& G& 0& 0\\ 0& 0& 0& 0& G& 0& 0& G& 0\\ 0& 0& 0& 0& 0& G& 0& 0& G\\ 0& 0& 0& G& 0& 0& G& 0& 0\\ 0& 0& 0& 0& G& 0& 0& G& 0\\ 0& 0& 0& 0& 0& G& 0& 0& G\end{array}\right]$$
$${\mathbf{D}}_{\mathbf{e}}^{\prime }+\mathbf{c}\mathbf{I}=\left[\begin{array}{ccccccccc}2G+c& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 2G+c& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 2G+c& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& G+c& 0& 0& G& 0& 0\\ 0& 0& 0& 0& G+c& 0& 0& G& 0\\ 0& 0& 0& 0& 0& G+c& 0& 0& G\\ 0& 0& 0& G& 0& 0& G+c& 0& 0\\ 0& 0& 0& 0& G& 0& 0& G+c& 0\\ 0& 0& 0& 0& 0& G& 0& 0& G+c\end{array}\right]$$
$${\left[{\mathbf{D}}_{\mathbf{e}}^{\prime }+\mathbf{c}\mathbf{I}\right]}^{-1}=\left[\begin{array}{ccccccccc}\frac{1}{2G+c}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \frac{1}{2G+c}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{1}{2G+c}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{G+c}{(2G+c)c}& 0& 0& \frac{-G}{(2G+c)c}& 0& 0\\ 0& 0& 0& 0& \frac{G+c}{(2G+c)c}& 0& 0& \frac{-G}{(2G+c)c}& 0\\ 0& 0& 0& 0& 0& \frac{G+c}{(2G+c)c}& 0& 0& \frac{-G}{(2G+c)c}\\ 0& 0& 0& \frac{-G}{(2G+c)c}& 0& 0& \frac{G+c}{(2G+c)c}& 0& 0\\ 0& 0& 0& 0& \frac{-G}{(2G+c)c}& 0& 0& \frac{G+c}{(2G+c)c}& 0\\ 0& 0& 0& 0& 0& \frac{-G}{(2G+c)c}& 0& 0& \frac{G+c}{(2G+c)c}\end{array}\right]$$
$$\begin{array}{rl}& {\mathbf{D}}_{\mathbf{e}}^{\prime }:(c\mathbf{I}+{\mathbf{D}}_{\mathbf{e}}^{\prime }{)}^{-1}:{\mathbf{D}}_{\mathbf{e}}^{\prime }=\\ & \left[\begin{array}{ccccccccc}\frac{4G^2}{2G+c}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \frac{4G^2}{2G+c}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{4G^2}{2G+c}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{2G^2}{(2G+c)}& 0& 0& \frac{2G^2}{(2G+c)}& 0& 0\\ 0& 0& 0& 0& \frac{2G^2}{(2G+c)}& 0& 0& \frac{2G^2}{(2G+c)}& 0\\ 0& 0& 0& 0& 0& \frac{2G^2}{(2G+c)}& 0& 0& \frac{2G^2}{(2G+c)}\\ 0& 0& 0& \frac{2G^2}{(2G+c)}& 0& 0& \frac{2G^2}{(2G+c)}& 0& 0\\ 0& 0& 0& 0& \frac{2G^2}{(2G+c)}& 0& 0& \frac{2G^2}{(2G+c)}& 0\\ 0& 0& 0& 0& 0& \frac{2G^2}{(2G+c)}& 0& 0& \frac{2G^2}{(2G+c)}\end{array}\right]\end{array}$$
$$\begin{array}{rl}& {\mathbf{D}}_{\mathbf{e}}^{\prime }-{\mathbf{D}}_{\mathbf{e}}^{\prime }:(c\mathbf{I}+{\mathbf{D}}_{\mathbf{e}}^{\prime }{)}^{-1}:{\mathbf{D}}_{\mathbf{e}}^{\prime }=\\ & \left[\begin{array}{ccccccccc}\frac{2Gc}{2G+c}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \frac{2Gc}{2G+c}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \frac{2Gc}{2G+c}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}& 0& 0\\ 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}& 0\\ 0& 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}\\ 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}& 0& 0\\ 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}& 0\\ 0& 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}\end{array}\right]\end{array}$$
$$\begin{array}{rl}& \left[{\mathbf{D}}_{\mathbf{e}}^{\prime }-{\mathbf{D}}_{\mathbf{e}}^{\prime }:(c\mathbf{I}+{\mathbf{D}}_{\mathbf{e}}^{\prime }{)}^{-1}:{\mathbf{D}}_{\mathbf{e}}^{\prime }\right]:(\mathbf{I}-\frac{1}{3}\hat{\mathbf{I}})+K\hat{\mathbf{I}}=\\ & \left[\begin{array}{ccccccccc}K+\frac{4Gc}{3(2G+c)}& K+\frac{-2Gc}{3(2G+c)}& K+\frac{-2Gc}{3(2G+c)}& 0& 0& 0& 0& 0& 0\\ K+\frac{-2Gc}{3(2G+c)}& K+\frac{4Gc}{3(2G+c)}& K+\frac{-2Gc}{3(2G+c)}& 0& 0& 0& 0& 0& 0\\ K+\frac{-2Gc}{3(2G+c)}& K+\frac{-2Gc}{3(2G+c)}& K+\frac{4Gc}{3(2G+c)}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}& 0& 0\\ 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}& 0\\ 0& 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}\\ 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}& 0& 0\\ 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}& 0\\ 0& 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0& \frac{Gc}{(2G+c)}\end{array}\right]\end{array}$$
Voigt表記下的一致剛度矩陣（注：應變為工程應變${\gamma }_{xy}=2{\epsilon }_{xy}$）:
$$\left[\begin{array}{cccccc}K+\frac{4Gc}{3(2G+c)}& K+\frac{-2Gc}{3(2G+c)}& K+\frac{-2Gc}{3(2G+c)}& 0& 0& 0\\ K+\frac{-2Gc}{3(2G+c)}& K+\frac{4Gc}{3(2G+c)}& K+\frac{-2Gc}{3(2G+c)}& 0& 0& 0\\ K+\frac{-2Gc}{3(2G+c)}& K+\frac{-2Gc}{3(2G+c)}& K+\frac{4Gc}{3(2G+c)}& 0& 0& 0\\ 0& 0& 0& \frac{Gc}{(2G+c)}& 0& 0\\ 0& 0& 0& 0& \frac{Gc}{(2G+c)}& 0\\ 0& 0& 0& 0& 0& \frac{Gc}{(2G+c)}\end{array}\right]$$

umat 程式

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
     1 RPL,DDSDDT,DRPLDE,DRPLDT,
     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),
     4 JSTEP(4),STRESS_TR(NTENS),STRESS_VALID(NTENS),DSTRAN_P(NTENS)

      EE=PROPS(1)  !modulus
      EMU=PROPS(2)   !Possion's ritio
      SGM_Y=PROPS(3)  !yield stress
      EC=PROPS(4)  !hardening 
      LMD=EE*EMU/((1.0+EMU)*(1.0-2.0*EMU))
      EG=EE/(2.0*(1.0+EMU))
      EK=EE/(3.0*(1.0-2.0*EMU))

      DDSDDE=0.0
      DO I=1,NDI
        DDSDDE(I,I)=2*EG
      ENDDO

      DO I=NDI+1,NTENS
        DDSDDE(I,I)=EG
      ENDDO

      DO I=1,NDI
        DO J=1,NDI
          DDSDDE(I,J)=DDSDDE(I,J)+LMD
        ENDDO
      ENDDO
      DO I=1,NTENS
        STRESS_TR(I)=STRESS(I)
      ENDDO
      DO I=1,NTENS
        DO J=1,NTENS
          STRESS_TR(I)=STRESS_TR(I)+ DDSDDE(I,J)*DSTRAN(J)
        ENDDO
      ENDDO
C     計算有效屈服應力
C     STATEV(1) ~ STATEV(NTENS) 為背應力  
      DO I=1,NTENS
        STRESS_VALID(I)=STRESS_TR(I)-STATEV(I)
      ENDDO
      STRESS_V=0.0
      DO I=1,NDI
        STRESS_V=STRESS_V+STRESS_VALID(I)
      ENDDO
      STRESS_V=STRESS_V/3.0
      STRESS_MISES=0.0
      DO I=1,NDI
        STRESS_MISES=STRESS_MISES+(STRESS_VALID(I)-STRESS_V)**2.0
      ENDDO
      DO I=1+NDI,NTENS
        STRESS_MISES=STRESS_MISES+2.0*STRESS_VALID(I)**2.0
      ENDDO
      STRESS_MISES=SQRT(1.5*STRESS_MISES)
C     判斷是否屈服
      IF (STRESS_MISES .LE. SGM_Y) THEN
        DO I=1,NTENS
          STRESS(I)=STRESS_TR(I)
        ENDDO
      ELSE   !屈服
      DP=(STRESS_MISES-SGM_Y)/(EC+3.0*EG)
      DO I=1,NDI
        DSTRAN_P(I)=DP*1.5*(STRESS_VALID(I)-STRESS_V)/STRESS_MISES
      ENDDO
      DO I=1+NDI,NTENS
        DSTRAN_P(I)=DP*1.5*STRESS_VALID(I)/STRESS_MISES
      ENDDO
      DO I=1,NTENS     !應力更新
        DO J=1,NTENS
          STRESS(I)=STRESS(I)+ DDSDDE(I,J)*(DSTRAN(J)-DSTRAN_P(J))
        ENDDO
      ENDDO
      DO I=1,NTENS    !背應力更新  
        STATEV(I)=STATEV(I)+DSTRAN_P(I)*EC
      ENDDO
      ENDIF
      
      RETURN
      END

結果