1. Definitions and constitutive equations
Consider a small-strain material update from \(t_n\) to \(t_{n+1}=t_n+\Delta t\). The total strain increment \(\Delta\boldsymbol\varepsilon\) is prescribed. The previous stress and internal variables are held fixed during the update. Elasticity is isotropic, and temperature is fixed while differentiating the mechanical update.
1.1 Stress, strain, and tensor operators
The deviatoric stress and von Mises equivalent stress are
\[\boldsymbol\sigma'
=\boldsymbol\sigma-\frac13\operatorname{tr}(\boldsymbol\sigma)\boldsymbol{1}_2,
\qquad
\bar\sigma
=\sqrt{\frac32\,\boldsymbol\sigma':\boldsymbol\sigma'}
=\sqrt{\frac32}\,\|\boldsymbol\sigma'\|.
\]
Here \(\boldsymbol{1}_2\) is the second-order identity, and the symmetric fourth-order identity is
\[(\mathbb{1}_4)_{ijkl}
=\frac12(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}).
\]
Define the deviatoric projector
\[\mathbb{1}^{\mathrm{dev}}_4
=\mathbb{1}_4-\frac13\boldsymbol{1}_2\otimes\boldsymbol{1}_2,
\qquad
\boldsymbol e
=\mathbb{1}^{\mathrm{dev}}_4:\boldsymbol\varepsilon.
\]
For any symmetric tensor \(\boldsymbol a\),
\[\mathbb{1}^{\mathrm{dev}}_4:\boldsymbol a
=\boldsymbol a-\frac13\operatorname{tr}(\boldsymbol a)\boldsymbol{1}_2.
\]
The elastic stiffness is
\[\boxed{
\mathbb D^e
=\lambda\,\boldsymbol{1}_2\otimes\boldsymbol{1}_2+2G\mathbb{1}_4
=K\,\boldsymbol{1}_2\otimes\boldsymbol{1}_2
+2G\mathbb{1}^{\mathrm{dev}}_4,
}
\qquad
\lambda=K-\frac23G.
\]
The shear modulus may also be denoted \(\mu\); throughout, \(\mu=G\).
The unit deviatoric flow direction is
\[\boxed{
\boldsymbol n
=\frac{\boldsymbol\sigma'}{\|\boldsymbol\sigma'\|}
=\sqrt{\frac32}\frac{\boldsymbol\sigma'}{\bar\sigma}.
}
\]
It satisfies \(\boldsymbol n:\boldsymbol n=1\) and \(\operatorname{tr}\boldsymbol n=0\). The bold symbol \(\boldsymbol n\) is distinct from the scalar material exponent \(n\).
1.2 Flow rule
Let \(s>0\) be the scalar deformation resistance. The equivalent viscoplastic strain rate is
\[\boxed{
\dot{\bar e}^{\mathrm{vp}}
=A\exp\!\left(-\frac{Q}{k_BT}\right)
\left[
\sinh\!\left(\xi\frac{\|\boldsymbol\sigma'\|}{s}\right)
\right]^{1/m}
=:f.
}
\]
The thermal factor uses an activation energy \(Q\) expressed consistently with \(k_B\). At fixed temperature, define
\[A_T=A\exp\!\left(-\frac{Q}{k_BT}\right),
\qquad
X=\xi\frac{\|\boldsymbol\sigma'\|}{s}
=\xi\sqrt{\frac23}\frac{\bar\sigma}{s}.
\]
Then
\[\boxed{f(\bar\sigma,s)=A_T(\sinh X)^{1/m}.}
\]
The normalization of the tensor strain rate is
\[\dot{\bar e}^{\mathrm{vp}}
=\sqrt{\frac23\,\dot{\boldsymbol e}^{\mathrm{vp}}:
\dot{\boldsymbol e}^{\mathrm{vp}}}
=\sqrt{\frac23}\,\|\dot{\boldsymbol e}^{\mathrm{vp}}\|.
\]
Since \(\dot{\boldsymbol e}^{\mathrm{vp}}\) has direction \(\boldsymbol n\),
\[\dot{\boldsymbol e}^{\mathrm{vp}}
=\|\dot{\boldsymbol e}^{\mathrm{vp}}\|\boldsymbol n
=\sqrt{\frac32}\,\dot{\bar e}^{\mathrm{vp}}\boldsymbol n
=\sqrt{\frac32}\,f\boldsymbol n.
\]
Thus the scalar rate \(f\) and the norm of the tensor rate differ by \(\sqrt{3/2}\).
1.3 Hardening rule
The saturation resistance is
\[s^*
=\hat s\left[
\frac{\dot{\bar e}^{\mathrm{vp}}}{A}
\exp\!\left(\frac{Q}{k_BT}\right)
\right]^n
=\hat s\left(\frac{f}{A_T}\right)^n.
\]
Substituting the flow rule also gives
\[\boxed{
s^*=\hat s(\sinh X)^{n/m}.
}
\]
The resistance evolution is
\[\boxed{
\dot s
=h_0\left|1-\frac{s}{s^*}\right|^a
\operatorname{sgn}\!\left(1-\frac{s}{s^*}\right)f
=:g.
}
\]
Introduce
\[z=1-\frac{s}{s^*},
\qquad
H=h_0|z|^a\operatorname{sgn}z.
\]
For \(z\ne0\),
\[H=h_0|z|^{a-1}z,
\qquad g=Hf.
\]
Both \(f\) and \(s^*\) depend on \(\bar\sigma\) and \(s\), so \(H\) and \(g\) inherit both dependencies. At \(z=0\), define \(H=0\) by continuity for \(a>0\).
2. Backward Euler stress update
2.1 Equivalent and tensor strain increments
Backward Euler integration gives
\[\Delta\bar e^{\mathrm{vp}}=f_{n+1}\Delta t,
\qquad
\boxed{
\Delta\boldsymbol e^{\mathrm{vp}}
=\sqrt{\frac32}\,\Delta\bar e^{\mathrm{vp}}\boldsymbol n_{n+1}
=\sqrt{\frac32}\,f_{n+1}\Delta t\,\boldsymbol n_{n+1}.
}
\]
Equivalently, write
\[\Delta\boldsymbol e^{\mathrm{vp}}
=\Delta\lambda\,\boldsymbol n_{n+1},
\qquad
\Delta\lambda
=\sqrt{\frac32}\,\Delta\bar e^{\mathrm{vp}}.
\]
The plastic multiplier \(\Delta\lambda\) is the norm of the tensor increment. It is related to the integral
\[\Delta\boldsymbol e^{\mathrm{vp}}
=\int_{t_n}^{t_{n+1}}\dot{\boldsymbol e}^{\mathrm{vp}}\,dt
=\int_{t_n}^{t_{n+1}}\dot\lambda\,\boldsymbol n\,dt
\approx\Delta\lambda\,\boldsymbol n_{n+1}.
\]
2.2 Elastic trial stress
The elastic predictor is
\[\boxed{
\boldsymbol\sigma^{\mathrm{tr}}
=\boldsymbol\sigma_n+\mathbb D^e:\Delta\boldsymbol\varepsilon.
}
\]
Taking its deviatoric part,
\[\boxed{
\boldsymbol\sigma'^{\mathrm{tr}}
=\boldsymbol\sigma'_n+2G\Delta\boldsymbol e,
\qquad
\Delta\boldsymbol e
=\mathbb{1}^{\mathrm{dev}}_4:\Delta\boldsymbol\varepsilon.
}
\]
Its equivalent magnitude is
\[\bar\sigma^{\mathrm{tr}}
=\sqrt{\frac32}\,\|\boldsymbol\sigma'^{\mathrm{tr}}\|.
\]
The strain decomposition gives
\[\begin{aligned}
\boldsymbol\sigma_{n+1}
&=\boldsymbol\sigma_n+
\mathbb D^e:(\Delta\boldsymbol\varepsilon-\Delta\boldsymbol e^{\mathrm{vp}})\\
&=\boldsymbol\sigma^{\mathrm{tr}}
-\mathbb D^e:\Delta\boldsymbol e^{\mathrm{vp}}.
\end{aligned}
\]
Because \(\operatorname{tr}(\Delta\boldsymbol e^{\mathrm{vp}})=0\),
\[(\boldsymbol{1}_2\otimes\boldsymbol{1}_2):
\Delta\boldsymbol e^{\mathrm{vp}}=\boldsymbol0,
\qquad
\mathbb{1}^{\mathrm{dev}}_4:\Delta\boldsymbol e^{\mathrm{vp}}
=\Delta\boldsymbol e^{\mathrm{vp}}.
\]
Therefore,
\[\boxed{
\boldsymbol\sigma_{n+1}
=\boldsymbol\sigma^{\mathrm{tr}}-2G\Delta\boldsymbol e^{\mathrm{vp}}.
}
\]
2.3 Hydrostatic stress and collinearity
Taking the trace shows
\[\operatorname{tr}\boldsymbol\sigma_{n+1}
=\operatorname{tr}\boldsymbol\sigma^{\mathrm{tr}}.
\]
Define the mean stress by \(p=\operatorname{tr}\boldsymbol\sigma/3\). Then \(p_{n+1}=p^{\mathrm{tr}}\). Removing the hydrostatic part gives
\[\boldsymbol\sigma'_{n+1}
=\boldsymbol\sigma'^{\mathrm{tr}}
-2G\Delta\lambda\,\boldsymbol n_{n+1}.
\]
Use \(\boldsymbol n_{n+1} =\boldsymbol\sigma'_{n+1}/\|\boldsymbol\sigma'_{n+1}\|\):
\[\left(1+\frac{2G\Delta\lambda}{\|\boldsymbol\sigma'_{n+1}\|}\right)
\boldsymbol\sigma'_{n+1}
=\boldsymbol\sigma'^{\mathrm{tr}}.
\]
For nonzero corrected deviatoric stress and \(\Delta\lambda\ge0\), the multiplier is positive. Hence the corrected and trial deviatoric stresses are collinear and have the same direction:
\[\boxed{
\boldsymbol n_{n+1}
=\frac{\boldsymbol\sigma'^{\mathrm{tr}}}
{\|\boldsymbol\sigma'^{\mathrm{tr}}\|}
=\sqrt{\frac32}
\frac{\boldsymbol\sigma'^{\mathrm{tr}}}{\bar\sigma^{\mathrm{tr}}}.
}
\]
Taking norms yields
\[\|\boldsymbol\sigma'_{n+1}\|
=\|\boldsymbol\sigma'^{\mathrm{tr}}\|-2G\Delta\lambda.
\]
Multiplying by \(\sqrt{3/2}\) and inserting
\(\Delta\lambda=\sqrt{3/2}f_{n+1}\Delta t\) gives the scalar stress equation
\[\boxed{
\bar\sigma_{n+1}
=\bar\sigma^{\mathrm{tr}}-3G\Delta t\,f_{n+1}.
}
\]
3. Coupled local residual equations
The two unknowns are \(\bar\sigma_{n+1}\) and \(s_{n+1}\). Suppress the subscript \(n+1\) on quantities evaluated at the end of the increment:
\[\boxed{
\begin{aligned}
R_1(\bar\sigma,s)
&=\bar\sigma-\bar\sigma^{\mathrm{tr}}+3G\Delta t\,f(\bar\sigma,s)=0,\\
R_2(\bar\sigma,s)
&=s-s_n-\Delta t\,g(\bar\sigma,s)=0.
\end{aligned}
}
\]
The previous resistance \(s_n\) and the trial stress \(\bar\sigma^{\mathrm{tr}}\) are fixed during the local solve.
After convergence,
\[\boxed{
\begin{aligned}
\Delta\bar e^{\mathrm{vp}}&=f\Delta t,\\
\Delta\boldsymbol e^{\mathrm{vp}}
&=\sqrt{\frac32}f\Delta t\,\boldsymbol n,\\
\boldsymbol\sigma_{n+1}
&=\boldsymbol\sigma^{\mathrm{tr}}
-2G\sqrt{\frac32}f\Delta t\,\boldsymbol n\\
&=p^{\mathrm{tr}}\boldsymbol{1}_2
+\frac{\bar\sigma}{\bar\sigma^{\mathrm{tr}}}
\boldsymbol\sigma'^{\mathrm{tr}}.
\end{aligned}
}
\]
The last form applies for \(\bar\sigma^{\mathrm{tr}}>0\).
4. Local Newton Jacobian
Let
\[J=\frac{\partial(R_1,R_2)}{\partial(\bar\sigma,s)}
=\begin{bmatrix}
\frac{\partial R_1}{\partial \bar \sigma}&\frac{\partial R_1}{\partial s}\\
\frac{\partial R_2}{\partial \bar \sigma}&\frac{\partial R_2}{\partial s}\\
\end{bmatrix}
=\begin{bmatrix}
J_{11}&J_{12}\\J_{21}&J_{22}
\end{bmatrix}.
\]
All partial derivatives in this section hold temperature fixed.
4.1 Derivatives of the flow rate
From \(X=\xi\sqrt{2/3}\,\bar\sigma/s\),
\[X_{\bar\sigma}=\frac{\xi}{s}\sqrt{\frac23},
\qquad
X_s=-\frac{X}{s}.
\]
Differentiate \(f=A_T(\sinh X)^{1/m}\):
\[\begin{aligned}
f_{\bar\sigma}
&=\frac{A_T}{m}(\sinh X)^{1/m-1}\cosh X\,
\frac{\xi}{s}\sqrt{\frac23}\\
&=\sqrt{\frac23}\frac{\xi}{ms}
\frac{\cosh X}{\sinh X}\,f.
\end{aligned}
\]
Thus
\[\boxed{
f_{\bar\sigma}
=\sqrt{\frac23}\frac{\xi}{ms}\coth X\,f,
\qquad
f_s=-\frac{\xi\|\boldsymbol\sigma'\|}{ms^2}\coth X\,f
=-\frac{X}{ms}\coth X\,f.
}
\]
Here \(\coth X=1/\tanh X\).
4.2 Derivatives of the saturation resistance
Since \(s^*=\hat s A_T^{-n}f^n\),
\[\boxed{
s^*_{\bar\sigma}=\frac{ns^*}{f}f_{\bar\sigma},
\qquad
s^*_s=\frac{ns^*}{f}f_s.
}
\]
For \(z=1-s/s^*\),
\[\begin{aligned}
z_{\bar\sigma}
&=\frac{s}{(s^*)^2}s^*_{\bar\sigma}
=\frac{ns}{s^*f}f_{\bar\sigma},\\
z_s
&=-\frac1{s^*}+\frac{s}{(s^*)^2}s^*_s
=-\frac1{s^*}+\frac{ns}{s^*f}f_s.
\end{aligned}
\]
4.3 Derivatives of the hardening rate
For \(z\ne0\),
\[\frac{d}{dz}\left(|z|^a\operatorname{sgn}z\right)
=a|z|^{a-1}.
\]
Therefore,
\[\boxed{
\begin{aligned}
H_{\bar\sigma}
&=h_0a|z|^{a-1}\frac{ns}{s^*f}f_{\bar\sigma},\\
H_s
&=h_0a|z|^{a-1}
\left(-\frac1{s^*}+\frac{ns}{s^*f}f_s\right).
\end{aligned}
}
\]
Since \(H=h_0|z|^{a-1}(s^*-s)/s^*\), an equivalent form for \(s^*\ne s\) is
\[\boxed{
\begin{aligned}
H_{\bar\sigma}
&=\frac{an sH}{f(s^*-s)}f_{\bar\sigma}
=\frac{an\xi\sqrt{2/3}}{m(s^*-s)}H\coth X,\\
H_s
&=\frac{aH}{s^*-s}
\left(-1+\frac{ns}{f}f_s\right).
\end{aligned}
}
\]
Finally, \(g=Hf\) implies
\[g_{\bar\sigma}=Hf_{\bar\sigma}+fH_{\bar\sigma},
\qquad
g_s=Hf_s+fH_s.
\]
4.4 Four entries of the Jacobian
Differentiating the residuals gives
\[\boxed{
\begin{aligned}
J_{11}&=1+3G\Delta t\,f_{\bar\sigma},\\
J_{12}&=3G\Delta t\,f_s,\\
J_{21}&=-\Delta t\left(Hf_{\bar\sigma}+fH_{\bar\sigma}\right),\\
J_{22}&=1-\Delta t\left(Hf_s+fH_s\right).
\end{aligned}
}
\]
For \(s^*\ne s\), the last two entries can be expanded as
\[\begin{aligned}
J_{21}
&=-\Delta t\left[
Hf_{\bar\sigma}
+\frac{an\xi\sqrt{2/3}}{m(s^*-s)}
Hf\coth X
\right],\\
J_{22}
&=1-\Delta t\left[
Hf_s+\frac{aH}{s^*-s}(nsf_s-f)
\right].
\end{aligned}
\]
At \(z=0\), the derivative \(H_z\) is \(0\) for \(a>1\) and \(h_0\) for \(a=1\). Use these limits in \(H_y=H_z z_y\). For \(0<a<1\), \(H\) is not differentiable there. The quotient expressions containing \(s^*-s\) apply away from that point.
5. Initialization, Newton corrections, and convergence
5.1 Initial stress estimate
The scalar relation
\(\bar\sigma^{\mathrm{tr}}-\bar\sigma=3G\Delta t\,f\)
provides a rate scale. Use
\[f_{\mathrm{guess}}
=\frac{0.998\,\bar\sigma^{\mathrm{tr}}}{3G\Delta t}.
\]
Choose a positive resistance estimate \(s_0\); a natural continuation value is \(s_0=s_n\). To invert the flow law at this resistance, set
\[y=\left(\frac{f_{\mathrm{guess}}}{A_T}\right)^m,
\qquad
\sinh\!\left(\xi\frac{\|\boldsymbol\sigma'\|_0}{s_0}\right)=y.
\]
The inverse hyperbolic sine identity is
\[\operatorname{asinh}y=\ln\left(y+\sqrt{y^2+1}\right).
\]
Hence
\[\boxed{
\|\boldsymbol\sigma'\|_0
=\frac{s_0}{\xi}\ln\left(y+\sqrt{y^2+1}\right),
\qquad
\bar\sigma^{(0)}
=\sqrt{\frac32}\,\|\boldsymbol\sigma'\|_0,
\qquad s^{(0)}=s_0.
}
\]
The factor \(0.998\) selects an initial rate. Newton iteration subsequently enforces both residual equations. The reference resistance value \(3.909\times10^8\) may be used when it belongs to the chosen material data and stress-unit system.
For zero previous deviatoric stress, the predictor also gives a strain-based interpretation:
\[\boldsymbol\sigma'^{\mathrm{tr}}=2G\Delta\boldsymbol e,
\qquad
\|\boldsymbol\sigma'^{\mathrm{tr}}\|=2G\|\Delta\boldsymbol e\|.
\]
Define the equivalent deviatoric strain-increment magnitude
\(\Delta\bar e=\sqrt{2/3}\|\Delta\boldsymbol e\|\). Then
\[\bar\sigma^{\mathrm{tr}}
=\sqrt{\frac32}\,2G\|\Delta\boldsymbol e\|
=3G\Delta\bar e.
\]
This relates the \(3G\) stress scale to the equivalent strain convention.
5.2 Newton correction
At iterate \(k\),
\[J
\begin{bmatrix}\delta\bar\sigma\\\delta s\end{bmatrix}
=-\begin{bmatrix}R_1\\R_2\end{bmatrix}.
\]
Define
\[D=\det J=J_{11}J_{22}-J_{12}J_{21}.
\]
For \(D\ne0\),
\[J^{-1}
=\frac1D\begin{bmatrix}
J_{22}&-J_{12}\\-J_{21}&J_{11}
\end{bmatrix}.
\]
Thus
\[\boxed{
\delta\bar\sigma=-\frac{J_{22}R_1-J_{12}R_2}{D},
\qquad
\delta s=-\frac{-J_{21}R_1+J_{11}R_2}{D}.
}
\]
5.3 Reference rate and admissibility
A numerical reference rate \(f_{\mathrm{ref}}=10\) defines a stress-to-resistance threshold. Let
\[y_{\mathrm{ref}}
=\left(\frac{f_{\mathrm{ref}}}{A_T}\right)^m,
\qquad
r_{\mathrm{ref}}
=\frac1\xi\operatorname{asinh}(y_{\mathrm{ref}}).
\]
Then
\[f\le f_{\mathrm{ref}}
\quad\Longleftrightarrow\quad
\frac{\|\boldsymbol\sigma'\|}{s}\le r_{\mathrm{ref}}
\quad\Longleftrightarrow\quad
\frac{\bar\sigma}{s}\le\sqrt{\frac32}\,r_{\mathrm{ref}}.
\]
The reference rate has the same inverse-time units as \(f\). It is an iteration safeguard; its value should accommodate the intended loading regime.
With the current iterate held fixed, construct a candidate
\[\bar\sigma_c=\bar\sigma^{(k)}+\eta\,\delta\bar\sigma,
\qquad
s_c=s^{(k)}+\eta\,\delta s,
\qquad \eta=1.
\]
If the candidate violates \(\bar\sigma_c\ge0\), \(s_c>0\), or the selected reference-rate bound, halve \(\eta\) and reconstruct the candidate from the same current iterate. Both corrections are scaled together:
\[\eta\leftarrow\frac12\eta.
\]
Allow at most 20 such cutbacks for a correction. If no admissible candidate is obtained, report failure of the material update. If the candidate is admissible, accept it and continue the local iteration.
5.4 Convergence criteria
Use the initial scales \(\|\boldsymbol\sigma'\|_0\) and \(s_0\) for the correction checks:
\[\boxed{
\left|\frac{\eta\,\delta\bar\sigma}{\|\boldsymbol\sigma'\|_0}\right|<10^{-9},
\qquad
\left|\frac{\eta\,\delta s}{s_0}\right|<10^{-9}.
}
\]
Also require the residuals evaluated at the accepted state to be small, for example
\[\frac{|R_1|}{\bar\sigma^{\mathrm{tr}}}<10^{-9},
\qquad
\frac{|R_2|}{s_0}<10^{-9},
\]
for positive reference scales. This ensures that small corrections correspond to a solution of the local equations. The formulas below concern the smooth branch with positive trial deviatoric stress, positive resistance, and a nonsingular local Jacobian.
6. Consistent tangent: dependence of the flow rate on strain
The consistent material tangent differentiates the converged discrete stress update:
\[\mathbb C^{\mathrm{ep}}
=\frac{\partial\boldsymbol\sigma_{n+1}}
{\partial\boldsymbol\varepsilon_{n+1}}
=\frac{\partial\Delta\boldsymbol\sigma}
{\partial\Delta\boldsymbol\varepsilon}.
\]
These derivatives are equal because the previous state is held fixed. In the following differentiation, write \(\boldsymbol\varepsilon\) for the variable end-of-increment strain.
From
\[\boldsymbol\sigma
=\boldsymbol\sigma^{\mathrm{tr}}-2G\Delta\boldsymbol e^{\mathrm{vp}},
\qquad
\Delta\boldsymbol e^{\mathrm{vp}}
=\sqrt{\frac32}\Delta t\,f\boldsymbol n,
\]
the product rule gives
\[\boxed{
\mathbb C^{\mathrm{ep}}
=\mathbb D^e
-2G\sqrt{\frac32}\Delta t
\left[
\boldsymbol n\otimes
\frac{\partial f}{\partial\boldsymbol\varepsilon}
+f\frac{\partial\boldsymbol n}{\partial\boldsymbol\varepsilon}
\right].
}
\]
The two required derivatives are obtained separately.
6.1 Differentiate the local residuals
At the converged state, \(R_1=R_2=0\). Their differentials satisfy
\[\begin{aligned}
dR_1
&=J_{11}\,d\bar\sigma+J_{12}\,ds
-d\bar\sigma^{\mathrm{tr}}=0,\\
dR_2
&=J_{21}\,d\bar\sigma+J_{22}\,ds=0.
\end{aligned}
\]
Therefore,
\[J\begin{bmatrix}d\bar\sigma\\ds\end{bmatrix}
=\begin{bmatrix}d\bar\sigma^{\mathrm{tr}}\\0\end{bmatrix},
\]
and
\[\boxed{
\begin{bmatrix}d\bar\sigma\\ds\end{bmatrix}
=\frac1D
\begin{bmatrix}J_{22}\\-J_{21}\end{bmatrix}
d\bar\sigma^{\mathrm{tr}}.
}
\]
The strain sensitivities of the local unknowns are
\[\frac{\partial\bar\sigma}{\partial\boldsymbol\varepsilon}
=\frac{J_{22}}D
\frac{\partial\bar\sigma^{\mathrm{tr}}}{\partial\boldsymbol\varepsilon},
\qquad
\frac{\partial s}{\partial\boldsymbol\varepsilon}
=-\frac{J_{21}}D
\frac{\partial\bar\sigma^{\mathrm{tr}}}{\partial\boldsymbol\varepsilon}.
\]
6.2 Apply the chain rule to the flow rate
Since \(f=f(\bar\sigma,s)\),
\[\begin{aligned}
\frac{\partial f}{\partial\boldsymbol\varepsilon}
&=f_{\bar\sigma}
\frac{\partial\bar\sigma}{\partial\boldsymbol\varepsilon}
+f_s\frac{\partial s}{\partial\boldsymbol\varepsilon}\\
&=\frac{f_{\bar\sigma}J_{22}-f_sJ_{21}}D
\frac{\partial\bar\sigma^{\mathrm{tr}}}{\partial\boldsymbol\varepsilon}.
\end{aligned}
\]
Define the scalar sensitivity
\[\boxed{
\beta
=\frac{f_{\bar\sigma}J_{22}-f_sJ_{21}}D
=\frac{df}{d\bar\sigma^{\mathrm{tr}}},
}
\]
where the last derivative follows the converged coupled solution.
The determinant admits the factorizations
\[D=J_{22}\left(J_{11}-\frac{J_{12}J_{21}}{J_{22}}\right)
=J_{21}\left(\frac{J_{11}J_{22}}{J_{21}}-J_{12}\right).
\]
Where the individual denominators are nonzero,
\[\boxed{
\beta
=\frac{f_{\bar\sigma}}
{J_{11}-J_{12}J_{21}/J_{22}}
+\frac{f_s}
{J_{12}-J_{11}J_{22}/J_{21}}.
}
\]
For \(J_{21}=0\) and \(D\ne0\),
\[\boxed{\beta=\frac{f_{\bar\sigma}}{J_{11}}.}
\]
The determinant form also applies when one of the factorized expressions is undefined.
7. Derivatives of the trial magnitude and flow direction
7.1 Derivative of the trial deviatoric stress
The predictor is
\[\boldsymbol\sigma'^{\mathrm{tr}}
=\boldsymbol\sigma'_n
+2G\mathbb{1}^{\mathrm{dev}}_4:\Delta\boldsymbol\varepsilon.
\]
Hence
\[\boxed{
\frac{\partial\boldsymbol\sigma'^{\mathrm{tr}}}
{\partial\boldsymbol\varepsilon}
=2G\mathbb{1}^{\mathrm{dev}}_4.
}
\]
7.2 Derivative of the equivalent trial stress
Introduce the scalar squared norm
\[w=\boldsymbol\sigma'^{\mathrm{tr}}:\boldsymbol\sigma'^{\mathrm{tr}},
\qquad
\bar\sigma^{\mathrm{tr}}=\sqrt{\frac32}\,w^{1/2}.
\]
First,
\[dw=2\boldsymbol\sigma'^{\mathrm{tr}}:
d\boldsymbol\sigma'^{\mathrm{tr}}.
\]
Then
\[\begin{aligned}
d\bar\sigma^{\mathrm{tr}}
&=\sqrt{\frac32}\frac12w^{-1/2}\,dw\\
&=\sqrt{\frac32}\,w^{-1/2}
\boldsymbol\sigma'^{\mathrm{tr}}:
\left(2G\mathbb{1}^{\mathrm{dev}}_4:d\boldsymbol\varepsilon\right).
\end{aligned}
\]
Because the trial deviatoric stress is trace-free,
\[\boldsymbol\sigma'^{\mathrm{tr}}:\mathbb{1}^{\mathrm{dev}}_4
=\boldsymbol\sigma'^{\mathrm{tr}},
\qquad
w^{-1/2}=\sqrt{\frac32}\frac1{\bar\sigma^{\mathrm{tr}}}.
\]
Consequently,
\[\boxed{
\frac{\partial\bar\sigma^{\mathrm{tr}}}
{\partial\boldsymbol\varepsilon}
=\frac{3G}{\bar\sigma^{\mathrm{tr}}}\boldsymbol\sigma'^{\mathrm{tr}}
=3G\sqrt{\frac23}\,\boldsymbol n
=\sqrt6\,G\boldsymbol n.
}
\]
Combining this with Section 6,
\[\boxed{
\frac{\partial f}{\partial\boldsymbol\varepsilon}
=\beta\,\frac{3G}{\bar\sigma^{\mathrm{tr}}}
\boldsymbol\sigma'^{\mathrm{tr}}
=\beta\sqrt6\,G\boldsymbol n.
}
\]
7.3 Derivative of the unit flow direction
Since
\[\boldsymbol n
=\sqrt{\frac32}
\frac{\boldsymbol\sigma'^{\mathrm{tr}}}{\bar\sigma^{\mathrm{tr}}},
\]
the quotient rule gives
\[\begin{aligned}
\frac{\partial\boldsymbol n}{\partial\boldsymbol\varepsilon}
&=\sqrt{\frac32}\left[
\frac1{\bar\sigma^{\mathrm{tr}}}
\frac{\partial\boldsymbol\sigma'^{\mathrm{tr}}}
{\partial\boldsymbol\varepsilon}
-\frac1{(\bar\sigma^{\mathrm{tr}})^2}
\boldsymbol\sigma'^{\mathrm{tr}}\otimes
\frac{\partial\bar\sigma^{\mathrm{tr}}}
{\partial\boldsymbol\varepsilon}
\right]\\
&=\sqrt{\frac32}\left[
\frac{2G}{\bar\sigma^{\mathrm{tr}}}\mathbb{1}^{\mathrm{dev}}_4
-\frac{3G}{(\bar\sigma^{\mathrm{tr}})^3}
\boldsymbol\sigma'^{\mathrm{tr}}\otimes\boldsymbol\sigma'^{\mathrm{tr}}
\right].
\end{aligned}
\]
Use
\[\boldsymbol\sigma'^{\mathrm{tr}}
=\sqrt{\frac23}\,\bar\sigma^{\mathrm{tr}}\boldsymbol n,
\]
so that
\[\boldsymbol\sigma'^{\mathrm{tr}}\otimes\boldsymbol\sigma'^{\mathrm{tr}}
=\frac23(\bar\sigma^{\mathrm{tr}})^2
\boldsymbol n\otimes\boldsymbol n.
\]
Therefore,
\[\boxed{
\frac{\partial\boldsymbol n}{\partial\boldsymbol\varepsilon}
=\sqrt{\frac32}\frac{2G}{\bar\sigma^{\mathrm{tr}}}
\left(\mathbb{1}^{\mathrm{dev}}_4-\boldsymbol n\otimes\boldsymbol n\right)
=\frac{2G}{\|\boldsymbol\sigma'^{\mathrm{tr}}\|}
\left(\mathbb{1}^{\mathrm{dev}}_4-\boldsymbol n\otimes\boldsymbol n\right).
}
\]
The projector in parentheses removes both volumetric changes and changes parallel to \(\boldsymbol n\); only changes in direction contribute to \(d\boldsymbol n\).
8. Collecting the consistent tangent
Insert the two derivatives into
\[\mathbb C^{\mathrm{ep}}
=\mathbb D^e
-2G\sqrt{\frac32}\Delta t
\left[
\boldsymbol n\otimes\frac{\partial f}{\partial\boldsymbol\varepsilon}
+f\frac{\partial\boldsymbol n}{\partial\boldsymbol\varepsilon}
\right].
\]
The flow-rate term contributes
\[-2G\sqrt{\frac32}\Delta t\,
\boldsymbol n\otimes(\beta\sqrt6\,G\boldsymbol n)
=-6G^2\Delta t\,\beta\,\boldsymbol n\otimes\boldsymbol n.
\]
The direction term contributes
\[\begin{aligned}
&-2G\sqrt{\frac32}\Delta t\,f
\sqrt{\frac32}\frac{2G}{\bar\sigma^{\mathrm{tr}}}
\left(\mathbb{1}^{\mathrm{dev}}_4-\boldsymbol n\otimes\boldsymbol n\right)\\
&\qquad
=-\frac{6G^2\Delta t\,f}{\bar\sigma^{\mathrm{tr}}}
\left(\mathbb{1}^{\mathrm{dev}}_4-\boldsymbol n\otimes\boldsymbol n\right).
\end{aligned}
\]
Thus
\[\begin{aligned}
\mathbb C^{\mathrm{ep}}
={}&K\boldsymbol{1}_2\otimes\boldsymbol{1}_2
+2G\mathbb{1}^{\mathrm{dev}}_4\\
&-\frac{6G^2\Delta t\,f}{\bar\sigma^{\mathrm{tr}}}
\mathbb{1}^{\mathrm{dev}}_4\\
&+\left(
\frac{6G^2\Delta t\,f}{\bar\sigma^{\mathrm{tr}}}
-6G^2\Delta t\,\beta
\right)\boldsymbol n\otimes\boldsymbol n.
\end{aligned}
\]
8.1 Coefficient of the deviatoric projector
Define
\[\boxed{
2G_{\mathrm{eff}}
=2G\left(1-\frac{3G\Delta t\,f}{\bar\sigma^{\mathrm{tr}}}\right)
=2G\left(1-\frac{3G\Delta\bar e^{\mathrm{vp}}}
{\bar\sigma^{\mathrm{tr}}}\right).
}
\]
At the converged solution, \(R_1=0\) gives
\[\boxed{
G_{\mathrm{eff}}=G\frac{\bar\sigma}{\bar\sigma^{\mathrm{tr}}}.
}
\]
8.2 Coefficient of the flow-direction dyad
Define
\[\begin{aligned}
\alpha
&=6G^2\Delta t
\left(\frac{f}{\bar\sigma^{\mathrm{tr}}}-\beta\right)\\
&=2G\left[
\frac{3G\Delta\bar e^{\mathrm{vp}}}{\bar\sigma^{\mathrm{tr}}}
-3G\Delta t\,\beta
\right].
\end{aligned}
\]
In terms of the local Jacobian,
\[\boxed{
\alpha
=6G^2\Delta t
\left[
\frac{f}{\bar\sigma^{\mathrm{tr}}}
-\frac{f_{\bar\sigma}J_{22}-f_sJ_{21}}
{J_{11}J_{22}-J_{12}J_{21}}
\right].
}
\]
The factor \(\beta\) includes the response of both local unknowns to the trial stress.
8.3 Complete material tangent
The final fourth-order tensor is
\[\boxed{
\mathbb C^{\mathrm{ep}}
=K\boldsymbol{1}_2\otimes\boldsymbol{1}_2
+2G_{\mathrm{eff}}
\left(\mathbb{1}_4-\frac13\boldsymbol{1}_2\otimes\boldsymbol{1}_2\right)
+\alpha\,\boldsymbol n\otimes\boldsymbol n.
}
\]
Equivalently,
\[\boxed{
\mathbb C^{\mathrm{ep}}
=\left(K-\frac23G_{\mathrm{eff}}\right)
\boldsymbol{1}_2\otimes\boldsymbol{1}_2
+2G_{\mathrm{eff}}\mathbb{1}_4
+\alpha\,\boldsymbol n\otimes\boldsymbol n.
}
\]
The scalar coefficient definitions can be stored as
\[\begin{aligned}
v_{1}&=3G\Delta t,\\
v_{2}&=v_{1}\beta,\\
\alpha
&=2G\left(
\frac{v_{1}f}{\bar\sigma^{\mathrm{tr}}}-v_{2}
\right).
\end{aligned}
\]
A useful consistency identity follows from \(J_{11}=1+3G\Delta t f_{\bar\sigma}\) and \(J_{12}=3G\Delta t f_s\):
\[1-3G\Delta t\,\beta
=\frac{J_{22}}D.
\]
Therefore,
\[2G_{\mathrm{eff}}+\alpha=2G\frac{J_{22}}D.
\]
The tangent acts on a deviatoric perturbation parallel to \(\boldsymbol n\) with this coefficient; for a deviatoric perturbation orthogonal to \(\boldsymbol n\), its coefficient is \(2G_{\mathrm{eff}}\).
9. Matrix representation with engineering shear strains
9.1 Stress and strain vectors
Use the component order \((11,22,33,12,23,13)\):
\[\boldsymbol\sigma_V
=\begin{bmatrix}
\sigma_{11}&\sigma_{22}&\sigma_{33}&\sigma_{12}&\sigma_{23}&\sigma_{13}
\end{bmatrix}^{\mathsf T},
\]
\[\boldsymbol\varepsilon_V^{\mathrm{eng}}
=\begin{bmatrix}
\varepsilon_{11}&\varepsilon_{22}&\varepsilon_{33}&
\gamma_{12}&\gamma_{23}&\gamma_{13}
\end{bmatrix}^{\mathsf T},
\qquad \gamma_{ij}=2\varepsilon_{ij}\quad(i\ne j).
\]
The strain components follow
\[\varepsilon_{ii}=\frac{\partial u_i}{\partial x_i}
\quad\text{(no sum)},
\qquad
\gamma_{ij}
=\frac{\partial u_i}{\partial x_j}
+\frac{\partial u_j}{\partial x_i}.
\]
The incremental relation is
\[d\boldsymbol\sigma_V
=C_V^{\mathrm{ep}}\,d\boldsymbol\varepsilon_V^{\mathrm{eng}}.
\]
For tensor-component vectors
\(a_V=(a_{11},a_{22},a_{33},a_{12},a_{23},a_{13})^{\mathsf T}\),
\[\boldsymbol a:\boldsymbol b=a_V^{\mathsf T}Wb_V,
\qquad W=\operatorname{diag}(1,1,1,2,2,2).
\]
In particular, the flow direction is normalized by
\[n_V^{\mathsf T}Wn_V=1,
\qquad
n_V=\begin{bmatrix}n_{11}&n_{22}&n_{33}&n_{12}&n_{23}&n_{13}\end{bmatrix}^{\mathsf T}.
\]
The matrix mapping engineering strain components to the tensor components of the deviatoric strain is
\[P^{\mathrm{eng}}
=\begin{bmatrix}
2/3&-1/3&-1/3&0&0&0\\
-1/3& 2/3&-1/3&0&0&0\\
-1/3&-1/3& 2/3&0&0&0\\
0&0&0&1/2&0&0\\
0&0&0&0&1/2&0\\
0&0&0&0&0&1/2
\end{bmatrix}.
\]
Let \(v=(1,1,1,0,0,0)^{\mathsf T}\). The volumetric term is \(Kvv^{\mathsf T}\), whose normal \(3\times3\) block has every entry equal to \(K\). Consequently,
\[D_V^e=Kvv^{\mathsf T}+2GP^{\mathrm{eng}}.
\]
Its shear diagonal entries are \(G\).
For the rank-one contribution,
\[C^\alpha_{ijkl}=\alpha n_{ij}n_{kl},
\qquad
d\sigma^\alpha_{ij}
=\alpha n_{ij}(n_{kl}d\varepsilon_{kl}).
\]
Expand the contraction:
\[\begin{aligned}
n_{kl}d\varepsilon_{kl}
={}&n_{11}d\varepsilon_{11}
+n_{22}d\varepsilon_{22}
+n_{33}d\varepsilon_{33}\\
&+2n_{12}d\varepsilon_{12}
+2n_{23}d\varepsilon_{23}
+2n_{13}d\varepsilon_{13}\\
={}&n_V^{\mathsf T}d\boldsymbol\varepsilon_V^{\mathrm{eng}}.
\end{aligned}
\]
For example, the two off-diagonal contributions combine as
\[C^\alpha_{ij12}d\varepsilon_{12}
+C^\alpha_{ij21}d\varepsilon_{21}
=2\alpha n_{ij}n_{12}d\varepsilon_{12}
=\alpha n_{ij}n_{12}d\gamma_{12}.
\]
Thus
\[\boxed{C_V^\alpha=\alpha n_Vn_V^{\mathsf T}.}
\]
There is no extra shear factor in this outer product when the input vector contains engineering shear strains.
9.4 Complete matrix
Combining the volumetric, deviatoric, and dyadic terms,
\[\boxed{
C_V^{\mathrm{ep}}
=Kvv^{\mathsf T}+2G_{\mathrm{eff}}P^{\mathrm{eng}}
+\alpha n_Vn_V^{\mathsf T}.
}
\]
Explicitly,
\[C_V^{\mathrm{ep}}
=
\begin{bmatrix}
K+\frac43G_{\mathrm{eff}}&K-\frac23G_{\mathrm{eff}}&K-\frac23G_{\mathrm{eff}}&0&0&0\\
K-\frac23G_{\mathrm{eff}}&K+\frac43G_{\mathrm{eff}}&K-\frac23G_{\mathrm{eff}}&0&0&0\\
K-\frac23G_{\mathrm{eff}}&K-\frac23G_{\mathrm{eff}}&K+\frac43G_{\mathrm{eff}}&0&0&0\\
0&0&0&G_{\mathrm{eff}}&0&0\\
0&0&0&0&G_{\mathrm{eff}}&0\\
0&0&0&0&0&G_{\mathrm{eff}}
\end{bmatrix}
+\alpha n_Vn_V^{\mathsf T}.
\]
If tensor shear strains are stored instead, define
\[\boldsymbol\varepsilon_V^{\mathrm{ten}}
=(\varepsilon_{11},\varepsilon_{22},\varepsilon_{33},
\varepsilon_{12},\varepsilon_{23},\varepsilon_{13})^{\mathsf T}.
\]
Since \(\boldsymbol\varepsilon_V^{\mathrm{eng}} =W\boldsymbol\varepsilon_V^{\mathrm{ten}}\), the corresponding matrix is
\[C_V^{\mathrm{ten}}=C_V^{\mathrm{ep}}W.
\]
Its shear columns are twice those of the engineering-shear matrix. In particular, the isotropic shear diagonal becomes \(2G_{\mathrm{eff}}\), and the dyadic term becomes \(\alpha n_Vn_V^{\mathsf T}W\).
10. Material-update sequence
- Elastic predictor. Compute \(\boldsymbol\sigma^{\mathrm{tr}}\), \(\boldsymbol\sigma'^{\mathrm{tr}}\), \(\bar\sigma^{\mathrm{tr}}\), and the unit direction \(\boldsymbol n\).
- Initial estimate. Choose \(s_0>0\), obtain the stress estimate from the inverse flow law, and establish the convergence scales and reference-rate bound.
- Local equations. At each Newton iterate, evaluate \(f\), \(s^*\), \(H\), \(g\), both residuals, and all four Jacobian entries.
- Newton correction. Solve for \((\delta\bar\sigma,\delta s)\), apply the common cutback factor if required, and reevaluate the accepted candidate.
- Convergence. Require both correction and residual convergence. Report failure if admissibility or local convergence cannot be achieved.
- State update. Store the converged resistance, accumulate \(\Delta\bar e^{\mathrm{vp}}=f\Delta t\), and reconstruct \(\boldsymbol\sigma_{n+1}\) using the unchanged trial mean stress and corrected deviatoric magnitude.
- Consistent tangent. Evaluate the Jacobian at the converged state, calculate \(\beta\), \(G_{\mathrm{eff}}\), and \(\alpha\), and assemble \(\mathbb C^{\mathrm{ep}}\) or its engineering-shear matrix.
The closed-form derivatives assume positive stress and resistance wherever a denominator requires them. At zero trial deviatoric stress, the normalized direction is undefined; the stress update and tangent must be evaluated through the appropriate constitutive limits whenever those limits exist.
浙公网安备 33010602011771号