5. Numerical formulation

For the variational formulation, it is the same as that given in the note [R5.03.21] and which refers to the law of behavior with isotropic work hardening and the Von Mises criterion in large deformations. We only recall that this is an Eulerian formulation, with geometry updated at each increment and at each iteration, and that behavioral rigidity and geometric rigidity are taken into account.

We now present the numerical integration of the law of behavior and give the expression for the tangent matrix (options FULL_MECA and RIGI_MECA_TANG).

5.1. Expression of the discretized model

Knowing the stress \({\sigma }^{-}\), the cumulative plastic deformation \({p}^{-}\), the elastic deformations \({e}^{-}\), the displacements \({u}^{-}\) and \(\Delta u\), we seek to determine \((\sigma ,p,e)\).

Since the displacements are known, the gradients of the transformation from \({\Omega }_{0}\) to \({\Omega }^{-}\), noted \({F}^{-}\), and from \({\Omega }^{-}\) to \(\Omega\), noted \(\Delta F\), are known.

To integrate this behavior model, an implicit Euler diagram is used, since porosity is an explicit function of the deformation via equation 4.3-1, which is therefore known when the behavior is integrated.

Once discretized, the following system is then obtained:

\(F=\Delta F{F}^{-}\) eq 5.1-1

\(J=\text{det}F\) eq 5.1-2

\(J\sigma =\tau\) eq 5.1-3

\(\tau =sb{}^{e}\text{}\) eq 5.1-4

\({b}^{e}=\text{Id}-2e\) eq 5.1-5

• Equations of state:

\(s=-\left[2\mu \tilde{e}+K\text{tr}e\text{Id}+\mathrm{3K}\alpha \Delta T\text{Id}\right]\) eq. 5.1-6

\(A=-R(p)\) eq 5.1-7

Subsequently, the flow laws and the plasticity criterion are expressed directly as a function of the elastic deformation tensor \(e\).

• Flow laws

\(\begin{array}{}{D}^{p}\simeq -\frac{1}{2}F{\dot{G}}^{p}{F}^{T}=-\frac{1}{2\Delta t}\left[\underset{{b}^{e}}{\underset{\underbrace{}}{{\text{FG}}^{p}{F}^{T}}}-\Delta F\underset{{b}^{e-}}{\underset{\underbrace{}}{{F}^{-}{G}^{\text{p-}}{F}^{-T}}}\Delta {F}^{T}\right]\\ \text{}=-\frac{1}{2\Delta t}\left[\text{Id}-2e-\Delta F\left\{\text{Id}-2{e}^{-}\right\}\Delta {F}^{T}\right]\\ \text{}=(e-\underset{{e}^{\mathrm{Tr}}}{\underset{\underbrace{}}{\frac{1}{2}\left[\text{Id}-\Delta F\left\{\text{Id}-2{e}^{-}\right\}\Delta {F}^{T}\right]}})/\Delta t=(e-{e}^{\text{Tr}})/\Delta t\end{array}\) eq 5.1-8

Taking the trace and deviatory parts of the equation [éq 4.1-5], we get:

\(\text{tr}e-\text{tr}{e}^{\mathrm{Tr}}=\Delta pDf\text{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(\frac{-K\text{tr}e}{{\sigma }_{1}})\) eq 5.1-9

\(\tilde{e}=\{\begin{array}{}\tilde{{e}^{\text{Tr}}}-\frac{3}{2}\Delta p\frac{\tilde{e}}{{e}_{\text{eq}}}\text{}\text{si solution régulière}\\ 0\text{et}\mathrm{\Delta p}\ge \frac{2}{3}(\tilde{e}-{\tilde{e}}^{\text{Tr}}{)}_{\text{eq}}\text{}\text{si solution singulière}\end{array}\) eq 5.1-10

• Consistency requirements

\(\begin{array}{}F=\{\begin{array}{cc}2\mu {e}_{\mathrm{eq}}+{\sigma }_{1}Df\mathrm{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})\mathrm{exp}(\frac{-K\mathrm{tr}e}{{\sigma }_{1}})-R-{\sigma }_{y}& \text{si solution régulière}\\ {\sigma }_{1}Df\mathrm{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})\mathrm{exp}(\frac{-K\mathrm{tr}e}{{\sigma }_{1}})-R-{\sigma }_{y}& \text{si solution singulière}\end{array}\\ \text{avec}F\le 0\Delta p\ge 0F\Delta p=0\end{array}\) eq 5.1-11

5.2. Nonlinear system resolution

Integrating the law of behavior therefore boils down to solving the [éq 5.1-9], [éq 5.1-10] and [éq 5.1-11] system. We will see that this resolution is reduced to that of a single scalar equation, whose unknown \(x\) is the increment of the trace of the elastic deformations:

\(x=\text{tr}e-\text{tr}{e}^{\mathrm{Tr}}\) eq 5.2-1

Thanks to this choice, whether the solution is elastic or plastic, reached at a singular point or not, the equation [éq 5.1-9] relating to the trace of the elastic increment is always valid and makes it possible to express the cumulative plastic deformation increment:

\(\begin{array}{}\text{tr}e-\text{tr}{e}^{\mathrm{Tr}}=\Delta p\underset{G}{\underset{\underbrace{}}{Df\mathrm{exp}(\frac{-K\text{tr}{e}^{\mathrm{Tr}}}{{\sigma }_{1}})\mathrm{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})}}\mathrm{exp}(\frac{-K(\text{tr}e-\text{tr}{e}^{\mathrm{Tr}})}{{\sigma }_{1}})\\ \Rightarrow \Delta p(x)=\frac{1}{G}x\text{exp}\frac{Kx}{{\sigma }_{1}}\end{array}\) eq 5.2-2

This equation shows us that we can look for \(x\ge 0\) to guarantee positive cumulative plastic deformation and that the elastic solution is obtained for \(x=0\). It is also noted that the cumulative plastic deformation increment is a continuous and strictly increasing function of \(x\). With these remarks, if we denote by S the term [éq 5.2-3] in the plasticity criterion, it is then, again, a continuous and strictly increasing function of \(x\):

\(F=2\mu {e}_{\text{eq}}-S(x)\text{}\text{avec}\text{}S(x)=-{\sigma }_{1}\text{G exp}(-\frac{\text{Kx}}{{\sigma }_{1}})+R(p(x))+{\sigma }_{y}\) eq 5.2-3

At this stage, the resolution process is divided into two stages.

5.2.1. Examining singular points

Such a singular point is characterized by [éq 5.1-10] (low) and [éq 5.1-11] (low), so in particular by \(S(x)=0\). Because of the properties of \(S\), this equation has at most one positive solution, let’s say \({x}^{S}\) that exists if and only if \(S(0)\le 0\). Knowledge of \({x}^{S}\) makes it possible to deduce the elastic deformation tensor \(e\), the cumulative plastic deformation \(p\) as well as the thermodynamic forces \(s\) and \(R\).

Finally, this singular point will be a solution if the inequality in [éq 5.1-11] (down) is true, that is to say if:

\({\mathrm{\Delta p}}^{s}\ge \frac{2}{3}({\tilde{e}}^{s}-{\tilde{e}}^{\text{Tr}}{)}_{\text{eq}}\) eq 5.2.1-1

5.2.2. Regular solution

The flow equation [éq 5.1-10] (top), which determines the deviatory part of the elastic deformation tensor, makes it possible to deduce a scalar equation that is a function of the cumulative plastic deformation increment:

\(\tilde{e}-{\tilde{e}}^{\text{Tr}}=-\frac{3}{2}\mathrm{\Delta p}\frac{\tilde{e}}{{e}_{\text{eq}}}\text{}\Rightarrow \text{}\{\begin{array}{}{e}_{\text{eq}}={e}_{\text{eq}}^{\text{Tr}}-\frac{3}{2}\Delta p\\ \tilde{e}=\frac{{e}_{\text{eq}}}{{e}_{\text{eq}}^{\text{Tr}}}{\tilde{e}}^{\text{Tr}}\end{array}\) eq 5.2.2-1

We note that due to the positivity of \({e}_{\text{eq}}\), the legal value of \(\Delta p\) is limited:

\(\Delta p\le \frac{2}{3}{e}_{\text{eq}}^{\text{Tr}}\) eq 5.2.2-2

The consistency condition now determines \(x\):

\(F=2\mu {e}_{\text{eq}}^{\text{Tr}}-S(x)-3\mu \Delta p\le 0\) eq 5.2.2-3

Given this expression, the increase in the legal value of \(\Delta p\) is reduced to the only condition \(S(x)\ge 0\) or, equivalently, to \(x\ge {x}^{S}\).

The elastic solution is obtained for \(x\) = 0. It is the solution to the problem if and only if:

\(F(0)=2\mu {e}_{\text{eq}}^{\text{Tr}}-S(0)<0\) eq 5.2.2-4

If not, we must then solve:

\(F(x)=2\mu {e}_{\text{eq}}^{\text{Tr}}-S(x)-\frac{3\mu }{G}x\text{exp}(\frac{\text{Kx}}{{\sigma }_{1}})=0\text{}\text{avec}\text{}\{\begin{array}{}x>{x}^{s}\text{}\text{si}{x}^{s}\text{existe}\\ x>0\text{}\text{sinon}\end{array}\) eq 5.2.2-5

This function is continuous and strictly decreasing and tends to \(-\infty\) with \(x\). It therefore admits at most one solution. The demonstration of the existence of this solution is immediate. Indeed, it is enough to prove that \(F\) is positive on the lower bound of the search interval.

When \({x}^{S}\) does not exist, \(F(0)>0\) since the solution is not elastic.

When \({x}^{S}\) exists, the function is equal to:

\(F({x}^{s})=2\mu {e}_{\text{eq}}^{\text{Tr}}-3\mu \Delta {p}^{s}>0\iff \Delta {p}^{s}<\frac{2}{3}{e}_{\text{eq}}^{\text{Tr}}\) eq 5.2.2-6

This condition is true since the singular solution has been rejected.

5.3. Calculation process

The procedure for solving all the equations of the model is as follows:

1. We are looking for whether the solution is elastic.

• calculation of \(F(0)\)

• if \(F(0)<0\), the solution of the problem is the elastic solution \({x}^{\text{Sol}}=0\)

• otherwise we go to 2)

  • • • • • • 2. If \(S(0)>0\), the solution is plastic and regular

• we go to 4)

3. If \(S(0)<0\), we’re looking for whether the solution is unique

• we solve \(S({x}^{s})=0\)

• if \({x}^{s}\) verifies the inequality \({\mathrm{\Delta p}}^{s}\ge \frac{2}{3}({\tilde{e}}^{s}-{\tilde{e}}^{\text{Tr}}{)}_{\text{eq}}\), then the solution is singular \({x}^{\text{Sol}}={x}^{s}\)

• otherwise, \({x}^{s}\) is a lower bound to solve \(F(x)=0\), we go to 4)

4. The solution is plastic and regular

• we solve \(F(x)=0\)

5.4. Resolution

To solve the two equations \(S(x)=0\) and \(F(x)=0\), we use a Newton method with controlled bounds coupled with dichothomy when Newton gives a solution outside the interval of the two terminals. We now present the determination of the limits for each of the preceding cases (points 2), 3) and 4) of the preceding paragraph).

Upper and lower limits in the case where the function S is strictly positive at the origin ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~

We solve:

\(\{\begin{array}{}F(x)=0\\ F(0)>0\end{array}\iff \{\begin{array}{}\underset{{f}_{1}}{\underset{\underbrace{}}{2\mu {e}_{\text{eq}}^{\text{Tr}}-S(x)}}=\underset{3\mu \Delta p}{\underset{\underbrace{}}{\frac{3\mu }{G}x\text{exp}(\frac{\text{Kx}}{{\sigma }_{1}})}}\\ {f}_{1}(0)>0\end{array}\) eq 5.4.1-1

where the function \(\mathrm{\Delta p}(x)\) is continuous, strictly increasing, and zero at the origin and the function \({f}_{1}(x)\) is continuous, strictly decreasing, and strictly positive at the origin (see [Figure 5.4.1-a]).

We ask:

\({f}_{1}=\underset{{f}_{2}}{\underset{\underbrace{}}{2\mu {e}_{\text{eq}}^{\text{Tr}}-R(x)-{\sigma }_{y}}}+{\sigma }_{1}\text{Gexp}(-\frac{\text{Kx}}{{\sigma }_{1}})\text{}\text{alors}\text{}{f}_{2}(x)<{f}_{1}(x)\text{}\forall x\ge 0\) eq 5.4.1-2

where function \({f}_{2}(x)\) is continuous, strictly decreasing. In this case, solving the equations:

\({f}_{2}(\Delta {p}^{\text{inf}})=3\mu \Delta {p}^{\text{inf}}\) and \({x}^{\text{inf}}\text{exp}(\frac{{\text{Kx}}^{\text{inf}}}{{\sigma }_{1}})=G\Delta {p}^{\text{inf}}\) eq 5.4.1-3

to derive from this successively \(\mathrm{\Delta p}\) then \(x\) gives a lower bound \({x}^{\text{Inf}}\) which corresponds to the solution of the isotropic work hardening model and Von Mises criterion. If \({f}_{2}(0)<0\), the lower bound is taken to be zero: \({x}^{\text{inf}}=0\).

The upper bound \({x}^{\text{Sup}}\) is such that:

\({x}^{\text{Sup}}\text{exp}(\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma }_{1}})=\frac{G}{3\mu }{f}_{1}({x}^{\text{Inf}})\) eq 5.4.1-4

The equation of type \(x\text{exp}(\frac{\text{Kx}}{{\sigma }_{1}})=\text{constante}\) is solved by a Newton method.

Figure 5.4.1-a: graphical representation of the upper and lower boundaries

5.4.1. Upper and lower bounds in the case where the function S is negative or zero at the origin

The system to be solved is as follows:

\(\{\begin{array}{}S(x)=0\\ S(0)<0\end{array}\text{}\iff \text{}\{\begin{array}{}R({p}^{-}+\frac{x}{G}\text{exp}(\frac{\text{Kx}}{{\sigma }_{1}}))+{\sigma }_{y}={\sigma }_{1}\text{Gexp}(-\frac{\text{Kx}}{{\sigma }_{1}})\\ R({p}^{-})+{\sigma }_{y}<{\sigma }_{1}G\end{array}\) eq 5.4.2-1

The left part is a continuous function, strictly increasing from \(x\) and strictly positive at the origin, the part on the right is a continuous function, strictly decreasing from \(x\) and strictly positive at the origin. Using the properties of these two functions, a graphical representation (cf. [Figure 5.4.2-a]) of these functions shows that the upper bound \({x}^{\text{Sup}}\) is such that:

\({\sigma }_{1}\text{G exp}(-\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma }_{1}})=R({p}^{-})+{\sigma }_{y}\text{}\iff \text{}{x}^{\text{Sup}}=\frac{{\sigma }_{1}}{K}\text{log}(\frac{{\sigma }_{1}G}{R({p}^{-})+{\sigma }_{y}})\) eq 5.4.2-2

The lower bound \({x}^{\text{Inf}}\) is such that:

\(\begin{array}{}{\sigma }_{1}\text{G exp}(-\frac{{\text{Kx}}^{\text{Inf}}}{{\sigma }_{1}})=R({p}^{-}+\frac{{x}^{\text{Sup}}}{G}\text{exp}(\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma }_{1}}))+{\sigma }_{y}\\ \iff {x}^{\text{Inf}}={\langle \frac{{\sigma }_{1}}{K}\text{log}(\frac{{\sigma }_{1}G}{R({p}^{-}+\frac{{x}^{\text{Sup}}}{G}\text{exp}(\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma }_{1}}))+{\sigma }_{y}})\rangle }^{+}\end{array}\) eq 5.4.2-3

Figure 5.4.2-a: graphical representation of the upper and lower boundaries

Upper and lower bounds in the case where the function \(S\) is strictly negative at the origin and \({x}^{s}\) is not a solution ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The following system is solved:

\(\{\begin{array}{}F(x)=0\\ S(0)<0\\ S({x}^{s})=0\end{array}\text{}\iff \text{}\{\begin{array}{}\underset{{f}_{1}}{\underset{\underbrace{}}{2\mu {e}_{\text{eq}}^{\text{Tr}}-S(x)}}=\underset{3\mu \Delta p}{\underset{\underbrace{}}{\frac{3\mu }{G}x\text{exp}(\frac{\text{Kx}}{{\sigma }_{1}})}}\\ {f}_{1}(0)>0\\ 2\mu {e}_{\text{eq}}^{\text{Tr}}=\frac{3\mu }{G}{x}^{s}\text{exp}(\frac{{\text{Kx}}^{s}}{{\sigma }_{1}})\end{array}\) eq 5.4.3-1

Solution \({x}^{\text{Sol}}\) is such as \(S({x}^{\text{Sol}})>0\).

For the lower bound, we take \({x}^{\text{Inf}}={x}^{s}\). Given the properties of the functions \({f}_{1}\) (strictly decreasing) and \(3\mu \Delta p(x)\) (strictly decreasing) and (strictly increasing), the upper bound \({x}^{\text{Sup}}\) is such that (cf. [Figure 5.4.3-a]):

\({x}^{\text{Sup}}\text{exp}(\frac{{\text{Kx}}^{\text{Sup}}}{{\sigma }_{1}})=\frac{\mathrm{2G}}{3}{e}_{\text{eq}}^{\text{Tr}}\) eq 5.4.3-2

This equation is solved by a Newton method.

Figure 5.4.3-a: graphical representation of the upper and lower boundaries

5.5. Volume correction a posteriori

In the case of the Rousselier model, the change in volume plays a crucial role, so that the error made by approximation \(s=\tau\) in the flow equations can lead to a fairly accurate change in porosity, cf. [Bib2]. By following the proposal of this reference, we will correct the elastic deformation trace a posteriori (i.e. after having calculated all the quantities).

In fact, in the absence of approximation, the hydrostatic part of the flow equation leads to:

\(\frac{d}{\mathrm{dt}}(\mathrm{ln}{J}^{p})=\dot{p}Df\mathrm{exp}(\frac{\text{tr}\tau }{3{\sigma }_{1}})\) (eq. 5.5-1)

After discretization in time, we then obtain a corrected proposal for the changes in plastic and elastic volume:

\({J}_{\mathrm{corr}}^{p}={J}^{p\text{-}}\text{exp}(\text{tr}e-\text{tr}{e}^{\mathrm{Tr}});{J}_{\mathrm{corr}}^{e}=\frac{J}{{J}_{\mathrm{corr}}^{p}}\) (eq. 5.5-2)

We will then look for a new trace of elastic deformation such that the change in elastic volume corresponds to the value corrected above:

\({e}_{\mathrm{corr}}=\tilde{e}+t\mathrm{Id}\text{où}t\text{tel que :}\text{det}(\mathrm{Id}-2{e}_{\mathrm{corr}})={{J}_{\mathrm{corr}}^{e}}^{2}\) (eq. 5.5-3)

This leads to a polynomial of degree 3 in \(t\), whose solution closest to \(e\) will be chosen.

5.6. Behavior tangent matrix expression

One gives here the expression for the tangent matrix (option FULL_MECA during Newton iterations, option RIGI_MECA_TANG for the first iteration).

For option FULL_MECA, this is obtained by linearizing the system of equations that govern the law of behavior. The main lines of this linearization are given below.

For option RIGI_MECA_TANG, these are the same expressions as those given for FULL_MECA but with \(\Delta p=0\). In particular, we will have \(\Delta F=\text{Id}\).

The law of behavior can be put in the following general form:

\(\sigma =\sigma (\tau ,\Delta F)\) eq 5.6-1

\(\tau =\tau (e)\) eq 5.6-2

\(e=e({e}^{\mathrm{Tr}},f)\) eq 5.6-3

\({e}^{\text{Tr}}={e}^{\text{Tr}}(\Delta F)\) eq 5.6-4

The linearization of this system gives:

\(\delta \sigma =(\frac{\partial \sigma }{\partial \tau }:\frac{\partial \tau }{\partial e}:(\frac{\partial e}{\partial {e}^{\mathrm{Tr}}}:\frac{\partial {e}^{\mathrm{Tr}}}{\partial \Delta F}+\frac{\partial e}{\partial f}\frac{\partial f}{\partial J}\otimes \frac{\partial J}{\partial \Delta F})+\frac{\partial \sigma }{\partial \Delta F}):\delta \Delta F=H:\delta \Delta F\) eq 5.6-5

where \(H\) is the tangent matrix. Next, the five terms of the previous equation are determined separately.

In the linearization of the system, the tensor \(C\) defined below and the following two equations will often be used:

\({\mathrm{\delta a}}_{\text{ij}}=\frac{1}{2}({\delta }_{\text{ik}}{\delta }_{\text{jl}}+{\delta }_{\text{jk}}{\delta }_{\text{il}}){\mathrm{\delta a}}_{\text{kl}}\) eq 5.6-6

\(\delta {a}_{\text{pp}}=\delta {}_{\mathrm{kl}}\delta {a}_{\text{kl}}\) eq 5.5-7

\({C}_{\text{ijkl}}=\frac{1}{2}({\delta }_{\text{ik}}{\delta }_{\text{jl}}+{\delta }_{\text{jk}}{\delta }_{\text{il}})\) eq 5.6-8

• Calculation of \(\frac{\partial \sigma }{\partial \tau }\) and \(\frac{\partial \sigma }{\partial \Delta F}\)

The linearization of the relationship between the Cauchy constraint \(\sigma\) and the Kirchhoff constraint \(\tau\) gives:

\(J\sigma =\tau \text{}\iff \text{}\delta \sigma =\frac{1}{J}\delta \tau -(\frac{\sigma }{J}\otimes \frac{\partial J}{\partial \Delta F}):\delta \Delta F\) eq 5.6-9

Using the relationship [éq 5.6-6], we get for \(\frac{\partial \sigma }{\partial \tau }\):

\(\frac{\partial \sigma }{\partial \tau }=C\) eq 5.6-10

and for \(\frac{\partial \sigma }{\partial \Delta F}\):

\(\frac{\partial \sigma }{\partial \Delta F}=-\frac{\sigma }{J}\otimes \frac{\partial J}{\partial \Delta F}\) eq 5.6-11

with

\(\begin{array}{}\frac{\partial J}{\partial {F}_{\text{11}}}={\mathrm{\Delta F}}_{\text{22}}{\mathrm{\Delta F}}_{\text{33}}-{\mathrm{\Delta F}}_{\text{23}}{\mathrm{\Delta F}}_{\text{32}}\\ \frac{\partial J}{\partial {F}_{\text{22}}}={\mathrm{\Delta F}}_{\text{11}}{\mathrm{\Delta F}}_{\text{33}}-{\mathrm{\Delta F}}_{\text{13}}{\mathrm{\Delta F}}_{\text{31}}\\ \frac{\partial J}{\partial {F}_{\text{33}}}={\mathrm{\Delta F}}_{\text{11}}{\mathrm{\Delta F}}_{\text{22}}-{\mathrm{\Delta F}}_{\text{12}}{\mathrm{\Delta F}}_{\text{21}}\\ \frac{\partial J}{\partial {F}_{\text{12}}}={\mathrm{\Delta F}}_{\text{31}}{\mathrm{\Delta F}}_{\text{23}}-{\mathrm{\Delta F}}_{\text{33}}{\mathrm{\Delta F}}_{\text{21}}\text{}\frac{\partial J}{\partial {F}_{\text{21}}}={\mathrm{\Delta F}}_{\text{13}}{\mathrm{\Delta F}}_{\text{32}}-{\mathrm{\Delta F}}_{\text{33}}{\mathrm{\Delta F}}_{\text{12}}\\ \frac{\partial J}{\partial {F}_{\text{13}}}={\mathrm{\Delta F}}_{\text{21}}{\mathrm{\Delta F}}_{\text{32}}-{\mathrm{\Delta F}}_{\text{22}}{\mathrm{\Delta F}}_{\text{31}}\text{}\frac{\partial J}{\partial {F}_{\text{31}}}={\mathrm{\Delta F}}_{\text{12}}{\mathrm{\Delta F}}_{\text{23}}-{\mathrm{\Delta F}}_{\text{22}}{\mathrm{\Delta F}}_{\text{13}}\\ \frac{\partial J}{\partial {F}_{\text{23}}}={\mathrm{\Delta F}}_{\text{31}}{\mathrm{\Delta F}}_{\text{12}}-{\mathrm{\Delta F}}_{\text{11}}{\mathrm{\Delta F}}_{\text{32}}\text{}\frac{\partial J}{\partial {F}_{\text{32}}}={\mathrm{\Delta F}}_{\text{13}}{\mathrm{\Delta F}}_{\text{21}}-{\mathrm{\Delta F}}_{\text{11}}{\mathrm{\Delta F}}_{\text{23}}\end{array}\) eq 5.6-12

• Calculation of \(\frac{\partial \tau }{\partial e}\)

The relationship between Kirchhoff stress \(\tau\) and the elastic deformation tensor \(e\) is given by:

\(\begin{array}{}\tau =s{b}^{e}=-2\mu e-\lambda \text{Tr}e\text{Id}+4\mu ee+2\lambda (\text{tr}e)e\\ \text{-}\mathrm{3K}\alpha \Delta T\text{Id}+\mathrm{6K}\alpha \Delta Te\end{array}\) eq 5.6-13

After linearization, we obtain:

\(\delta \tau =2(\lambda \text{tr}e-\mu +\mathrm{3K}\alpha \Delta T)\delta e+\lambda (2e-\text{Id})\text{Tr}\delta e+4\mu (\delta ee+e\delta e)\) eq 5.6-14

From where

\(\frac{\partial {\tau }_{\text{ij}}}{\partial {e}_{\text{kl}}}=2(\lambda \text{tr}e-\mu +\mathrm{3K}\alpha \Delta T){C}_{\text{ijkl}}+\lambda (2{e}_{\text{ij}}-{\delta }_{\text{ij}}){\delta }_{\text{kl}}+2\mu ({\delta }_{\text{ik}}{e}_{\text{lj}}+{\delta }_{\text{il}}{e}_{\text{kj}}+{e}_{\text{il}}{\delta }_{\text{kj}}+{e}_{\text{ik}}{\delta }_{\text{jl}})\) eq 5.6-15

• Calculation of \(\frac{\partial {e}^{\text{Tr}}}{\partial \Delta F}\)

The relationship between the elastic strain tensor \({e}^{\text{Tr}}\) and the transformation gradient increment \(\Delta F\) is written as:

\({e}^{\text{Tr}}=\frac{1}{2}(\text{Id}-\Delta F{b}^{e-}\Delta {F}^{T})\) eq 5.6-16

Its linearization gives:

\(\frac{\partial {e}_{\text{ij}}^{\text{Tr}}}{\partial \Delta {F}_{\text{kl}}}=-\frac{1}{2}({\delta }_{\text{ik}}\Delta {F}_{\text{jp}}{b}_{\text{pl}}^{\text{e-}}+\Delta {F}_{\text{ip}}{b}_{\text{pl}}^{\text{e-}}{\delta }_{\text{jk}})\) eq 5.6-17

• Calculation of \(\frac{\partial e}{\partial {e}^{\text{Tr}}}\)

Elastic case

In the elastic case, the calculation of \(\frac{\partial e}{\partial {e}^{\text{Tr}}}\) is immediate since \(\delta e=\delta {e}^{\text{Tr}}\) from where

\(\frac{\partial e}{\partial {e}^{\text{Tr}}}=C\) eq 5.6-18

Plastic case — Regular solution

To determine \(\frac{\partial e}{\partial {e}^{\text{Tr}}}\), we operate in two steps. Using the discretized flow law, we first calculate \(\delta e\) as a function of \(\delta {e}^{\text{Tr}}\) and \(\delta \Delta p\). Then the consistency condition makes it possible to deduce \(\delta \Delta p\) as a function of \(\delta {e}^{\text{Tr}}\). These two steps are detailed below.

The deviatoric part of the discretized flow law is written as:

\(\tilde{e}-{\tilde{e}}^{\text{Tr}}=-\frac{3}{2}\Delta p\frac{\tilde{e}}{{e}_{\mathrm{eq}}}\) eq 5.6-19

After linearization, we obtain:

\(\underset{1/\alpha }{\underset{\underbrace{}}{(1+\frac{3}{2}\frac{\Delta p}{{e}_{\text{eq}}})}}\delta \tilde{e}=\delta {\tilde{e}}^{\text{Tr}}-\frac{3}{2}\frac{\tilde{e}}{{e}_{\text{eq}}}\delta \Delta p+\frac{9}{4}\Delta p\frac{\tilde{e}}{{e}_{\text{eq}}^{3}}(\tilde{e}:\delta \tilde{e})\) eq 5.6-20

To find \(\tilde{e}:\delta \tilde{e}\), we contract the equation [éq 5.6-20] with \(\tilde{e}\) which gives:

\(\tilde{e}:\delta \tilde{e}=\tilde{e}:\delta {\tilde{e}}^{\text{Tr}}-{e}_{\text{eq}}\delta \Delta p\) eq 5.6-21

From where

\(\delta \tilde{e}=\underset{{A}_{1}}{\underset{\underbrace{}}{\alpha \left[\frac{9\Delta p}{4{e}_{\text{eq}}^{3}}\tilde{e}\otimes \tilde{e}+C\right]}}:\delta {\tilde{e}}^{\text{Tr}}\underset{{A}_{2}}{\underset{\underbrace{}}{-\frac{3}{2}\frac{\tilde{e}}{{e}_{\text{eq}}}}}\mathrm{\delta \Delta p}\) eq 5.6-22

For the trace part of the discretized flow law, we have:

\(\text{Tr}e-\text{Tr}{e}^{\text{Tr}}=Df\Delta p\text{exp}(\frac{3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(-\frac{K}{{\sigma }_{1}}\text{Tr}e)\) eq 5.6-23

Which gives, by asking \({k}_{1}=1+\frac{DfK\Delta p}{{\sigma }_{1}}\text{exp}(\frac{3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(-\frac{K}{{\sigma }_{1}}\text{Tr}e)\):

\(\text{tr}\delta e=\begin{array}{}\underset{{\alpha }_{1}}{\underset{\underbrace{}}{\frac{1}{{k}_{1}}}}\text{tr}\delta {e}^{\mathrm{Tr}}\\ +\underset{{\alpha }_{2}}{\underset{\underbrace{}}{\frac{Df\text{exp}(\frac{-K}{{\sigma }_{1}}\text{tr}e)\text{exp}(\frac{3K\alpha \Delta T}{{\sigma }_{1}})}{{k}_{1}}}}\delta \Delta p\\ +\underset{{\beta }_{1}}{\underset{\underbrace{}}{\frac{\Delta p\text{exp}(\frac{-K}{{\sigma }_{1}}\text{tr}e)\text{exp}(\frac{3K\alpha \Delta T}{{\sigma }_{1}})}{{k}_{1}}}}D\delta f\end{array}\) eq 5.6-24

In the plastic case, the consistency condition is equal to:

\(2\mu {e}_{\text{eq}}+Df{\sigma }_{1}\text{exp}(-\frac{3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(-\frac{K}{{\sigma }_{1}}\text{Tr}e)-R-{\sigma }_{y}=0\) eq 5.6-25

From where

\(\begin{array}{}\frac{3\mu }{{e}_{\mathrm{eq}}}(\tilde{e}:\delta \tilde{e})+D{\sigma }_{1}\text{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(\frac{-K\text{tr}e}{{\sigma }_{1}})\delta f\\ -DfK\text{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(\frac{-K\text{tr}e}{{\sigma }_{1}})\text{tr}\delta e-h\delta \Delta p=0\end{array}\) eq 5.6-26

By injecting the relationship [éq 5.6-21] into the equation above, we then obtain, by setting

\({k}_{2}=3\mu +h+{\alpha }_{2}DfK\text{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(-\frac{K}{{\sigma }_{1}}\text{Tr}e)\):

\(\delta \Delta p=\begin{array}{}\frac{3\mu }{{e}_{\mathrm{eq}}}\underset{{\alpha }_{3}}{\underset{\underbrace{}}{\frac{1}{{k}_{2}}}}\tilde{e}:\delta {\tilde{e}}^{\mathrm{Tr}}-\underset{{\alpha }_{4}}{\underset{\underbrace{}}{\frac{{\alpha }_{1}DfK\text{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(-\frac{K}{{\sigma }_{1}}\text{Tr}e)}{{k}_{2}}}}\text{tr}\delta {e}^{\mathrm{Tr}}\\ \underset{{\beta }_{2}}{\underset{\underbrace{}}{+\frac{\text{exp}(\frac{-3K\alpha \Delta T}{{\sigma }_{1}})\text{exp}(\frac{-K}{{\sigma }_{1}}\text{Tr}e)({\sigma }_{1}-{\beta }_{1}DfK)}{{k}_{2}}}}D\delta f\end{array}\) eq 5.6-27

By replacing \(\delta \Delta p\) with its value obtained above in equations [éq 5.6–22] and [éq 5.6-24], we get:

\(\begin{array}{}\delta e=\underset{\mathrm{ddvetr}}{\underset{\underbrace{}}{\left[{A}_{1}+({A}_{2}+\frac{1}{3}{\alpha }_{2}\mathrm{Id})\otimes (\frac{3\mu {\alpha }_{3}}{{e}_{\mathrm{eq}}}\tilde{e})\right]}}:\delta \tilde{{e}^{\mathrm{Tr}}}+\underset{\mathrm{dtretr}}{\underset{\underbrace{}}{\left[\frac{1}{3}{\alpha }_{1}\mathrm{Id}+{\alpha }_{4}({A}_{2}+\frac{1}{3}{\alpha }_{2}\mathrm{Id})\right]}}\text{tr}\delta {e}^{\mathrm{Tr}}\\ +\underset{\mathrm{dedf}}{\underset{\underbrace{}}{\left[\frac{1}{3}{\beta }_{1}\mathrm{Id}+{\beta }_{2}({A}_{2}+\frac{1}{3}{\alpha }_{2}\mathrm{Id})\right]D}}\delta f\end{array}\) eq 5.6-28

From where

\(\frac{\partial e}{\partial {e}^{\text{Tr}}}=\text{ddvetr}+(\text{dtretr}-\frac{1}{3}\text{ddvetr}:\text{Id})\otimes \text{Id}\) eq 5.6-29

Plastic case — Singular solution

The approach is identical to that used previously.

For the discretized flow law, we obtain:

\(\tilde{e}=0\text{}\iff \text{}\delta \tilde{e}=0\) eq 5.6-30

for the deviatory part and for the trace part, the relationship is identical to that found for the regular solution.

\(\text{tr}\delta e={\alpha }_{1}\text{tr}\delta {e}^{\mathrm{Tr}}+{\alpha }_{2}\delta \Delta p+{\beta }_{1}D\delta f\) eq 5.6-31

where \({\alpha }_{1}\), \({\alpha }_{2}\) and \({\beta }_{1}\) have the same definitions as in the previous paragraph.

The consistency condition then makes it possible to find \(\delta \Delta p\) according to \(\partial {e}^{\text{Tr}}\).

\(Df{\sigma }_{1}\text{exp}(\frac{\mathrm{3K}\alpha \Delta T}{{\sigma }_{1}})\text{exp}(-\frac{K}{{\sigma }_{1}}\text{Tr}e)-R-{\sigma }_{y}=0\) eq 5.6-32

or after linearization:

\(\delta \Delta p={\alpha }_{4}\text{tr}\delta {e}^{\mathrm{Tr}}+{\beta }_{2}D\delta f\) eq 5.6-33

or finally:

\(\delta e=\underset{\mathrm{dtretr}}{\underset{\underbrace{}}{\frac{1}{3}\left[{\alpha }_{1}+{\alpha }_{4}{\alpha }_{2}\right]\mathrm{Id}\text{tr}\delta {e}^{\mathrm{Tr}}}}+\underset{\mathrm{dedf}}{\underset{\underbrace{}}{\frac{1}{3}\left[{\beta }_{1}+{\beta }_{2}{\alpha }_{2}\right]\mathrm{Id}D\delta f}}\) eq 5.6-34

From where

\(\frac{\partial e}{\partial {e}^{\text{Tr}}}=\text{dtretr}\otimes \text{Id}\) eq 5.6-35

• Calculation of \(\frac{\partial f}{\partial J}\)

Given the relationship 4.3-1, the derivative is simply written:

\(\{\begin{array}{cc}\frac{\partial f}{\partial J}=D\frac{1-{f}_{0}}{{J}^{2}}& \text{si}f>{f}_{0}\\ \frac{\partial f}{\partial J}=0& \text{si}f={f}_{0}\end{array}\) eq 5.6-36

Note:

The tangent matrix is not altered by the volume correction because the volume correction is performed a posteriori, i.e. after the constraints have been calculated.