【论文阅读】1998-StefanieReese-A Theory of Finite Viscoelasticity and Numerical Aspects-Sec3

3-Integration of the Evolution Equation

3.1 Operator split : isotropic

• Evolution eqn:

$$\begin{array}{r}\boxed{-\frac{1}{2}{\mathcal{L}}_v{b}_e\cdot b_e^{-1}={\mathcal{V}}^{-1}:{\tau }_{NEQ}}\end{array}$$

• rewrite it as the potential form:

$$\begin{array}{c}\begin{array}{r}-\frac{1}{2}{\mathcal{L}}_v{b}_e\cdot b_e^{-1}=\frac{1}{\tau }\frac{\partial }{\partial {\tau }_{NEQ}}(\frac{\tau }{2}{\tau }_{NEQ}:{\mathcal{V}}^{-1}:{\tau }_{NEQ})=\frac{1}{\tau }\frac{\partial {\Phi }_{vis}}{\partial {\tau }_{NEQ}}\end{array}\end{array}$$

• ${\Phi }_{vis}$: interpreted as creep potential
• Algorithm of the integration of evolution equations
• The main features of algorithms is to carry out an operator split of material derivative of $b_e$ into an elastic predictor $E$ and an inelastic corrector $I$:

$$\begin{array}{rl}{\dot{b}}_e& =\dot{\overline{(F\cdot C_i^{-1}\cdot F^T)}}\\ & =\dot{F}{F}_i^{-1}{F}_i^{-T}{F}^T+F\dot{\overline{(C_i^{-1})}}{F}^T+F{F}_i^{-1}{F}_i^{-T}{\dot{F}}^T\\ & =\dot{F}{F}^{-1}{F}_e{F}_e^T{F}^{-T}{F}^T+F{F}^{-1}{F}_e{F}_e^T{F}^{-T}{\dot{F}}^T+F\dot{\overline{(C_i^{-1})}}{F}^T\\ & =l{b}_e+b_{e}l+F\dot{\overline{(C_i^{-1})}}{F}^T\end{array}$$

$$\begin{array}{c}\begin{array}{r}{\dot{b}}_e=\underset{E}{\underbrace{l\cdot b_e+b_e\cdot l}}+\underset{I}{\underbrace{F\cdot \dot{\overline{(C_i^{-1})}}\cdot F^T}}\end{array}\end{array}$$

In elastic predictor step the material time derivative of $C_i^{-1}$ is set to zero

$$\begin{array}{c}\begin{array}{r}(C_i^{-1}{)}_{trial}=(C_i^{-1}{)}_{t=t_{n-1}}\Rightarrow (b_e{)}_{trial}=(F{)}_{t=t_n}\cdot (C_i^{-1}{)}_{t=t_{n-1}}\cdot (F{)}_{t=t_n}^T\end{array}\end{array}$$

• Analogy to elastoplasticity, the state calculated with the elastic predictor step is called the trial state

In inelastic corrector step, the spatial velocity gradient $l$ is set to zero, leads ${\mathcal{L}}_v{b}_e={\dot{b}}_e$

by:
$${\dot{b}}_e=\underset{I}{\underbrace{F\cdot \dot{\overline{(C_i^{-1})}}\cdot F^T}}={\mathcal{L}}_v{b}_e$$

$$\begin{array}{c}\begin{array}{r}{\mathcal{L}}_v{b}_e={\dot{b}}_e\end{array}\end{array}$$

by:
$$\begin{array}{r}-\frac{1}{2}{\mathcal{L}}_v{b}_e\cdot b_e^{-1}={\mathcal{V}}^{-1}:{\tau }_{NEQ}\end{array}$$

$$\begin{array}{c}\begin{array}{r}{\dot{b}}_e=-(2{\mathcal{V}}^{-1}:{\tau }_{NEQ})\cdot b_e\end{array}\end{array}$$

• Solving the differential equation by exponential mapping:

$$\begin{array}{c}\begin{array}{r}{b}_e=\exp [-2{\int }_{t_{n-1}}^t{\mathcal{V}}^{-1}:{\tau }_{NEQ}dt]\cdot b_e\end{array}\end{array}$$

$$\begin{array}{c}\begin{array}{r}(b_e{)}_{t=t_n}\approx \exp [-2\underset{\Delta t}{\underbrace{(t_n-t_{n-1})}}({\mathcal{V}}^{-1}:{\tau }_{NEQ}{)}_{t=t_n}]\cdot (b_e{)}_{trial}\end{array}\end{array}$$

• approximate expression is first order accurate

For isotropy case

• ${\tau }_{NEQ}$ is commutes with $b_e$ and also with $(b_e{)}_{trial}$
• Prove $b_e$ and $(b_e{)}_{trial}$ share the same set of eigenvectors in the case of isotropy
• Note that:

$$\begin{array}{c}\begin{array}{r}{b}_e=\underset{\mathsf{a}}{\underbrace{\exp [-2\Delta t{\mathcal{V}}^{-1}:{\tau }_{NEQ}]}}\cdot (b_e{)}_{trial}\end{array}\end{array}$$

• By assumption of isotropy of ${\Psi }_{NEQ}$ that ${\tau }_{NEQ}$ has the same eigenspace as $b_e$ , (isotropic tensor function representation theorem)
• Since ${\mathcal{V}}^{-1}$ is isotropic, ${\mathcal{V}}^{-1}:{\tau }_{NEQ}$ shares the same eigenspace as $b_e$
• $\mathsf{a}$ has the same Eigen space as $b_e$
• Because of the properties of exponential map on symmetric tensor, the existence of inverse of $\mathsf{a}$ is guaranteed

$$\begin{array}{c}\begin{array}{r}{\mathsf{a}}^{-1}\cdot b_e=(b_e{)}_{trial}\end{array}\end{array}$$

• ${\mathsf{a}}^{-1}\cdot b_e$ and $b_e$ shares the same eigenspace
• Thus, $b_e$ and $(b_e{)}_{trial}$ shares the same eigenspace
• Solve (8) in terms of principal values, finite viscoelastic evolution is:

by:

$$\begin{gathered} {\mathcal{V}}^{-1}=\frac{1}{2{\eta }_D}({\text{1}}^4-\frac{1}{3}\text{1}\otimes \text{1})+\frac{1}{9{\eta }_V}\text{1}\otimes \text{1} \\ {\mathcal{V}}^{-1}:{\tau }_{NEQ}=\frac{1}{2{\eta }_D}dev[{\tau }_{NEQ}]+\frac{1}{9{\eta }_V}{\tau }_{NEQ}:\text{1} \end{gathered}$$

$$\begin{array}{c}\begin{array}{r}{\lambda }_{Ae}^2=\exp [-\Delta t(\frac{1}{{\eta }_D}dev[{\tau }_A]+\frac{2}{9{\eta }_V}{\tau }_{NEQ}:\text{1})]({\lambda }_{Ae}^2{)}_{trial}\end{array}\end{array}$$

• $dev[{\tau }_A]$ are the principal values of $dev[{\tau }_{NEQ}]$
• ${\lambda }_{Ae}^2$ are the principal values of $(b_e{)}_{t=t_n}$
• $({\lambda }_{Ae}^2{)}_{trial}$ are the principal values of $(b_e{)}_{trial}$
• Take the Logarithm of both sides:

$$\begin{array}{c}\begin{array}{r}\boxed{{\epsilon }_{Ae}=-\Delta t(\frac{1}{2{\eta }_D}dev[{\tau }_A]+\frac{1}{9{\eta }_V}{\tau }_{NEQ}:\text{1})+({\epsilon }_{Ae}{)}_{trial}}\end{array}\end{array}$$

• ${\epsilon }_{Ae}=\ln {\lambda }_{Ae}^2$ : are the elastic logarithmic principal stretches
• (11) furnishes an implicit nonlinear equation for determining the principal values of $b_e$

3.2 Linearized case

$$\begin{array}{c}\begin{array}{r}\boxed{-\frac{1}{2}({\mathcal{L}}_v{b}_e\cdot b_e^{-1}{)}_{lin}=\frac{1}{2\tau }(b_e-1)}\end{array}\end{array}$$

• in inelastic corrector step:${\mathcal{L}}_v{b}_e={\dot{b}}_e$

$${\dot{b}}_e=-\frac{1}{\tau }(b_e-1)\cdot b_e$$

• Solving the differential equation by exponential mapping:

$$\begin{gathered} b_e=\exp [-\frac{1}{\tau }{\int }_{t_{n-1}}^t(b_e-1)dt]\cdot (b_e{)}_{trial} \\ (b_e{)}_{t=t_n}\approx \exp [-\frac{1}{\tau }\underset{\Delta t}{\underbrace{(t_n-t_{n-1})}}((b_e{)}_{t=t_n}-1)]\cdot (b_e{)}_{trial} \end{gathered}$$

• For isotropic case

$$\begin{array}{c}\begin{array}{r}{\lambda }_{Ae}^2=\exp [-\frac{\Delta t}{\tau }({\lambda }_{Ae}^2-1)]({\lambda }_{Ae}^2{)}_{trial}\end{array}\end{array}$$

• hence:

$$\begin{array}{c}\begin{array}{r}{\epsilon }_{Ae}=-\frac{\Delta t}{\tau }\frac{1}{2}({\lambda }_{Ae}^2-1)+({\epsilon }_{Ae}{)}_{trial}\end{array}\end{array}$$

• linearizes ${\epsilon }_{Ae}=\ln {\lambda }_{Ae}$ with respect to ${\lambda }_{Ae}^2$ around thermodynamic equilibrium ${\lambda }_{Ae}^2=1$:

$$\begin{array}{c}\begin{array}{rl}({\epsilon }_{Ae}{)}_{lin}& =[\frac{\partial (\ln {\lambda }_{Ae})}{\partial ({\lambda }_{Ae}^2)}{]}_{EQ}({\lambda }_{Ae}^2-1)\\ & =[\frac{\partial (\frac{1}{2}\ln {\lambda }_{Ae}^2)}{\partial ({\lambda }_{Ae}^2)}{]}_{EQ}({\lambda }_{Ae}^2-1)\\ & =[\frac{1}{2}\frac{1}{{\lambda }_{Ae}^2}{]}_{{\lambda }_{Ae}^2=1}({\lambda }_{Ae}^2-1)\\ & =\frac{1}{2}({\lambda }_{Ae}^2-1)\end{array}\end{array}$$

• evolution equation (12) is valid for small perturbation away from thermodynamic equilibrium, no restriction for (15)
• substitue (15) into (14):

$$\begin{array}{c}\begin{array}{r}\boxed{{\epsilon }_{Ae}=-\frac{\Delta t}{\tau }{\epsilon }_{Ae}+({\epsilon }_{Ae}{)}_{trial}}\end{array}\end{array}$$

• Solving for ${\epsilon }_{Ae}$:

$$\begin{array}{c}\begin{array}{r}{\epsilon }_{Ae}=(1+\frac{\Delta t}{\tau }{)}^{-1}({\epsilon }_{Ae}{)}_{trial}\end{array}\end{array}$$

• is explicit update formula for the logarithmic principal values of $b_e$

3.3 Local Newton Iteration

for isotropic case

• non-equilibrium free energy ${\Psi }_{NEQ}$ can be represented as a function of principal values $b_{Ae}={\lambda }_{Ae}^2$ of ${\mathbf{b}}_e$:

$$\begin{array}{c}\begin{array}{r}{\Psi }_{NEQ}={\Psi }_{NEQ}(b_{1e},b_{2e},b_{3e})\end{array}\end{array}$$

• Split the non-equilobrioum Kirchhoff stress tensor and the non-equilibrium part of free energy into volumetric and deviatoric parts:

$$\begin{array}{c}\begin{array}{r}{\Psi }_{NEQ}=\widehat{\overline{({\Psi }_{NEQ}{)}_D}}({\overline{b}}_{1e},{\overline{b}}_{2e},{\overline{b}}_{3e})+\widehat{\overline{({\Psi }_{NEQ}{)}_V}}(J_e)\end{array}\end{array}$$

• Here, $J_e$ denotes the determinant of $F_e$:

$$\begin{array}{c}\begin{array}{r}{J}_e={\lambda }_{1e}{\lambda }_{2e}{\lambda }_{3e}\end{array}\end{array}$$

• And the principal values of ${\overline{b}}_e$:

$$\begin{array}{c}\begin{array}{r}{\overline{b}}_{Ae}=J_e^{-(2/3)}{b}_{Ae}=J_e^{-(2/3)}{\lambda }_{Ae}^2\end{array}\end{array}$$

• where:

$$\begin{array}{c}\begin{array}{r}{\overline{b}}_e=h_e^{-(2/3)}{b}_e\end{array}\end{array}$$

Ogden Models

• ${\Psi }_{NEQ}$ split into the deviatoric and volumetric part :

$$\begin{array}{c}\begin{array}{r}{\Psi }_{NEQ}=\underset{({\Psi }_{NEQ}{)}_D}{\underbrace{\sum _{r=1}^3\frac{({\mu }_m{)}_r}{({\alpha }_m{)}_r}({\overline{b}}_{1e}^{({\alpha }_m{)}_r/2}+{\overline{b}}_{2e}^{({\alpha }_m{)}_r/2}+{\overline{b}}_{3e}^{({\alpha }_m{)}_r/2}-3)}}+\underset{({\Psi }_{NEQ}{)}_V}{\underbrace{\frac{K_m}{4}(J_e^2-2\ln J_e-1)}}\end{array}\end{array}$$

• $({\mu }_m{)}_r$ and $({\alpha }_m{)}_r$ are elasticity constant
• the shear modulus ${\mu }_m$ :

$$\begin{array}{c}\begin{array}{r}{\mu }_m=\sum _{r=1}^3\frac{1}{2}({\mu }_m{)}_r({\alpha }_m{)}_r\end{array}\end{array}$$

• For rubberlike material, we have $K_m\approx 1000({\mu }_m{)}_r$, which means that the material behavior is nearly incompressible
• Non-equilibrium part of the Kirchhoff stress tensor:

By:

$$A:(u\otimes v)=u\cdot Av=(u\otimes v):A$$

$$\begin{array}{rl}d{b}_e& =d(\sum _{A=1}^3{\lambda }_{Ae}^2{n}_A\otimes n_A)\\ & =\sum _{A=1}^3(2{\lambda }_{Ae}d{\lambda }_{Ae}{n}_A\otimes n_A+{\lambda }_{Ae}^2(d{n}_A\otimes n_A+n_A\otimes d{n}_A))\\ \Rightarrow & n_A\cdot d{b}_e\cdot n_A=2{\lambda }_{Ae}d{\lambda }_{Ae}\\ \because & n_A\cdot d{b}_e\cdot n_A=d{b}_e:(n_A\otimes n_A)\\ \Rightarrow & d{b}_e:(n_A\otimes n_A)=2{\lambda }_{Ae}d{\lambda }_{Ae}\\ \Rightarrow & \frac{\partial b_e}{\partial {\lambda }_{Ae}}d{\lambda }_{Ae}:(n_A\otimes n_A)=2{\lambda }_{Ae}d{\lambda }_{Ae}\\ \Rightarrow & \frac{1}{2{\lambda }_{Ae}}\frac{\partial b_e}{\partial {\lambda }_{Ae}}:(n_A\otimes n_A)=1\\ \because & 1=\frac{\partial b_e}{\partial {\lambda }_{Ae}}:\frac{\partial {\lambda }_{Ae}}{\partial b_e}\\ \Rightarrow & \frac{\partial {\lambda }_{Ae}}{\partial b_e}=\frac{1}{2{\lambda }_{Ae}}(n_A\otimes n_A)\end{array}$$

$$\begin{array}{c}\begin{array}{rl}{\tau }_{NEQ}& =2\frac{\partial {\Psi }_{NEQ}}{\partial b_e}\cdot b_e\\ & =2\sum _{A=1}^3\frac{\partial {\Psi }_{NEQ}}{\partial {\lambda }_{Ae}}\frac{\partial {\lambda }_{Ae}}{\partial b_e}\cdot {\lambda }_{Ae}^2{n}_A\otimes n_A\\ & =2\sum _{A=1}^3\frac{\partial {\Psi }_{NEQ}}{\partial {\lambda }_{Ae}}\frac{1}{2{\lambda }_{Ae}}(n_A\otimes n_A)\cdot {\lambda }_{Ae}^2{n}_A\otimes n_A\\ & =\sum _{A=1}^3\frac{\partial {\Psi }_{NEQ}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}{n}_A\otimes n_A\\ & =\sum _{A=1}^3{\tau }_A{n}_A\otimes n_A\end{array}\end{array}$$

• ${\tau }_A$ is the principal Kirchhoff stress :

By: (19), (20), (21)

$$\frac{\partial J_e}{\partial {\lambda }_{Ae}}=\frac{\partial ({\lambda }_{1e}{\lambda }_{2e}{\lambda }_{3e})}{\partial {\lambda }_{Ae}}=\frac{J_e}{{\lambda }_{Ae}}$$

$$\begin{array}{c}\begin{array}{rl}{\tau }_A& =\frac{\partial {\Psi }_{NEQ}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}\\ & =\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}+\frac{\partial ({\Psi }_{NEQ}{)}_V}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}\\ & =\sum _{B=1}^3\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_{Be}}\frac{\partial {\overline{b}}_{Be}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}+\frac{\partial ({\Psi }_{NEQ}{)}_V}{\partial J_e}\frac{\partial J_e}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}\\ & =\underset{dev{\tau }_A}{\underbrace{\sum _{B=1}^3\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_{Be}}\frac{\partial {\overline{b}}_{Be}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}}}+\underset{\frac{1}{3}{\tau }_{NEQ}:\text{1}}{\underbrace{\frac{\partial ({\Psi }_{NEQ}{)}_V}{\partial J_e}{J}_e}}\end{array}\end{array}$$

• Verify:

By: (22)
$$\begin{gathered} \frac{\partial \det A}{\partial A}=\det A{A}^{-1} \\ \Rightarrow \frac{\partial J_e^2}{\partial b_e}=J_e^2{b}_e^{-1}\Rightarrow \frac{\partial J_e}{\partial b_e}=\frac{J_e}{2}{b}_e^{-1} \\ \Rightarrow \frac{\partial J_e^{-2/3}}{\partial b_e}=-\frac{2}{3}{J}_e^{-5/3}\frac{\partial J_e}{\partial b_e}=-\frac{1}{3}{J}_e^{-2/3}{b}_e^{-1} \end{gathered}$$

$$\begin{array}{rl}& \frac{\partial (\Phi A)}{\partial C}=A\otimes \frac{\partial \Phi }{\partial C}+\Phi \frac{\partial A}{\partial C}\\ & \\ \Rightarrow & \frac{\partial {\overline{b}}_e}{\partial b_e}=\frac{\partial (J_e^{-2/3}{b}_e)}{\partial b_e}=J_e^{-2/3}\frac{\partial b_e}{\partial b_e}+b_e\otimes \frac{\partial J^{-2/3}}{\partial b_e}\\ & =J_e^{-2/3}{\text{1}}^4-\frac{1}{3}{J}_e^{-2/3}{b}_e\otimes b_e^{-1}\\ & =J_e^{-2/3}({\text{1}}^4-\frac{1}{3}{b}_e\otimes b_e^{-1})\\ \Rightarrow & 2\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial b_e}\cdot b_e=2\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_e}:\frac{\partial {\overline{b}}_e}{\partial b_e}\cdot b_e\\ & =2\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_e}\cdot J_e^{2/3}{\overline{b}}_e:(J_e^{-2/3}({\text{1}}^4-\frac{1}{3}{b}_e\otimes b_e^{-1}))\\ & =2\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_e}\cdot {\overline{b}}_e:({\text{1}}^4-\frac{1}{3}{b}_e\otimes b_e^{-1})\\ & =dev[2\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_e}\cdot {\overline{b}}_e]\\ & =dev[{\tau }_A]\end{array}$$

$$\begin{array}{r}\Rightarrow \sum _{B=1}^3\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_{Be}}\frac{\partial {\overline{b}}_{Be}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}=2\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_e}\cdot {\overline{b}}_e=dev[{\tau }_A]\end{array}$$

$$\begin{gathered} \because {\tau }_{NEQ}:\text{1}={\tau }_A:\text{1} \\ \Rightarrow \frac{\partial ({\Psi }_{NEQ}{)}_V}{\partial J_e}{J}_e=\frac{1}{3}{\tau }_{NEQ}:\text{1} \end{gathered}$$

• (11) is nonlinear equation, solved by means of Newton iteration, Local Newton Iteration
• (16) is the linearized evolution equation, is explicitly solved

Remark 4:

• ${\Psi }_{NEQ}$ represent strain energy in a single Maxwell element
• To model viscous of certain material it’s necessary to use several Maxwell elements parallel
• In this case, solve an evolution equation for each Maxwell element and additional contribution to non-equilibrium part of Kirchhoff stress tensor

Local Newton Iteration

• [1] Non-linear equation, from (11)

$$\begin{array}{c}\begin{array}{r}{r}_A={\epsilon }_{Ae}+\Delta t(\frac{1}{2{\eta }_D}dev[{\tau }_A]+\frac{1}{9{\eta }_V}{\tau }_{NEQ}:\text{1})-({\epsilon }_{Ae}{)}_{trial}=0\end{array}\end{array}$$

• [2] Linearize around ${\epsilon }_{Ae}=({\epsilon }_{Ae}{)}_k$

$$\begin{array}{c}\begin{array}{r}{r}_A\approx \underset{({\tilde{r}}_A{)}_k}{\underbrace{r_A{|}_{({\epsilon }_{Ae}{)}_k}}}+\sum _{B=1}^3\underset{(K_{AB}{)}_k}{\underbrace{\frac{\partial r_A}{\partial {\epsilon }_{Be}}{|}_{({\epsilon }_{Ae}{)}_k}}}(\Delta {\epsilon }_{Be}{)}_k=0\end{array}\end{array}$$

Calculation of $K_{AB}$:

$$K_{AB}=\frac{\partial r_A}{\partial {\epsilon }_{Be}}$$

First, principal stress (26) has to be determined:
from (23):
$$\begin{gathered} \frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_{Be}}=\sum _{r=1}^3\frac{({\mu }_m{)}_r}{2}{\overline{b}}_{Be}^{[({\alpha }_m{)}_r/2-1]} \\ \frac{\partial ({\Psi }_{NEQ}{)}_V}{\partial J_e}=\frac{K_m}{2}(J_e-\frac{1}{J_e}) \end{gathered}$$
from (21):

$$\begin{gathered} \frac{\partial {\overline{b}}_{Be}}{\partial {\lambda }_{Ae}}=\frac{\partial (J_e^{-2/3}{\lambda }_{Be}^2)}{\partial {\lambda }_{Ae}}=-\frac{2}{3}{J}_e^{-2/3}\frac{{\lambda }_{Be}^2}{{\lambda }_{Ae}} \\ \frac{\partial {\overline{b}}_{Ae}}{\partial {\lambda }_{Ae}}=\frac{\partial (J_e^{-2/3}{\lambda }_{Ae}^2)}{\partial {\lambda }_{Ae}}=\frac{4}{3}{J}_e^{-2/3}{\lambda }_{Ae} \end{gathered}$$

from (26) & (21): $A\ne B,A\ne C,B\ne C$

$$\begin{array}{rl}{\tau }_A& =\sum _{B=1}^3\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_{Be}}\frac{\partial {\overline{b}}_{Be}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}+\frac{\partial ({\Psi }_{NEQ}{)}_V}{\partial J_e}{J}_e\\ & =\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_{Ae}}\frac{\partial {\overline{b}}_{Ae}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}+\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_{Be}}\frac{\partial {\overline{b}}_{Be}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}+\frac{\partial ({\Psi }_{NEQ}{)}_D}{\partial {\overline{b}}_{Ce}}\frac{\partial {\overline{b}}_{Ce}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}\\ & +\frac{\partial ({\Psi }_{NEQ}{)}_V}{\partial J_e}{J}_e\\ & =\underset{dev[{\tau }_A]}{\underbrace{\sum _{r=1}^3({\mu }_m{)}_r(\frac{2}{3}{\overline{b}}_{Ae}^{[({\alpha }_m{)}_r/2]}-\frac{1}{3}{\overline{b}}_{Be}^{[({\alpha }_m{)}_r/2]}-\frac{1}{3}{\overline{b}}_{Ce}^{[({\alpha }_m{)}_r/2]})}}+\underset{\frac{1}{3}{\tau }_{NEQ}:\text{1}}{\underbrace{\frac{K_m}{2}(J_e^2-1)}}\end{array}$$

By means of the relation:
from ${\lambda }_{Ae}=\exp [{\epsilon }_{Ae}]$

$$\begin{gathered} \frac{\partial {\overline{b}}_{Ae}}{\partial {\epsilon }_{Ae}}=\frac{\partial {\overline{b}}_{Ae}}{\partial {\lambda }_{Ae}}\frac{\partial {\lambda }_{Ae}}{\partial {\epsilon }_{Ae}}=\frac{\partial {\overline{b}}_{Ae}}{\partial {\lambda }_{Ae}}{\lambda }_{Ae}=\frac{4}{3}{\overline{b}}_{Ae} \\ \frac{\partial {\overline{b}}_{Ae}}{\partial {\epsilon }_{Be}}=\frac{\partial {\overline{b}}_{Ae}}{\partial {\lambda }_{Be}}\frac{\partial {\lambda }_{Be}}{\partial {\epsilon }_{Be}}=\frac{\partial {\overline{b}}_{Ae}}{\partial {\lambda }_{Be}}{\lambda }_{Be}=-\frac{2}{3}{\overline{b}}_{Ae}.(A\ne B) \\ \frac{\partial J_e}{\partial {\epsilon }_{Be}}=\frac{\partial J_e}{\partial {\lambda }_{Be}}\frac{\partial {\lambda }_{Be}}{\partial {\epsilon }_{Be}}=\frac{J_e}{{\lambda }_{Be}}{\lambda }_{Be}=J_e \end{gathered}$$

We arrive: $A\ne B,A\ne C,B\ne C$

$$\begin{array}{rl}& \frac{\partial dev[{\tau }_A]}{\partial {\epsilon }_{Ae}}=\sum _{r=1}^3({\mu }_m{)}_r({\alpha }_m{)}_r(\frac{4}{9}{\overline{b}}_{Ae}^{[({\alpha }_m{)}_r/2]}+\frac{1}{9}{\overline{b}}_{Be}^{[({\alpha }_m{)}_r/2]}+\frac{1}{9}{\overline{b}}_{Ce}^{[({\alpha }_m{)}_r/2]})\\ & \frac{\partial dev[{\tau }_A]}{\partial {\epsilon }_{Be}}=\sum _{r=1}^3({\mu }_m{)}_r({\alpha }_m{)}_r(-\frac{2}{9}{\overline{b}}_{Ae}^{[({\alpha }_m{)}_r/2]}-\frac{2}{9}{\overline{b}}_{Be}^{[({\alpha }_m{)}_r/2]}+\frac{1}{9}{\overline{b}}_{Ce}^{[({\alpha }_m{)}_r/2]})\\ & \frac{\partial (\frac{1}{3}{\tau }_{NEQ}:\text{1})}{\partial {\epsilon }_{Ae}}=K_m{J}_e^2\end{array}$$

Assume ${\eta }_D$ and ${\eta }_V$ to be constant

$$\begin{array}{r}{K}_{AB}=\frac{\partial r_A}{\partial {\epsilon }_{Ae}}={\delta }_{AB}+\Delta t\frac{1}{{\eta }_D}\frac{\partial dev[{\tau }_A]}{\partial {\epsilon }_{Be}}-\Delta t\frac{1}{3{\eta }_V}{K}_m{J}_e^2\end{array}$$

• [3] Solve for $(\Delta {\epsilon }_{Be}{)}_k$:

$$\begin{array}{c}\begin{array}{r}\sum _{B=1}^3(K_{AB}{)}_k(\Delta {\epsilon }_{Be}{)}_k=-({\tilde{r}}_A{)}_k\end{array}\end{array}$$

• [4] Update:

$$\begin{array}{c}\begin{array}{r}({\epsilon }_{Ae}{)}_{k+1}=({\epsilon }_{Ae}{)}_k+(\Delta {\epsilon }_{Ae}{)}_k,k\leftarrow k+1\end{array}\end{array}$$

• [5] Repeat until $||r_A||<tol$