4. Integrating the behavioral relationship

4.1. Establishment of the scalar equation for the implicit schema and with constant elastic coefficients

\({\mathrm{\epsilon }}_{\text{tot}}\) refers to the total deformation at time \(t+\Delta t\) and by \(\mathrm{\Delta }{\mathrm{\epsilon }}_{\text{tot}}\) to the total deformation variation during the current time step.

We call \({\mathrm{\epsilon }}_{\mathrm{o}}\) the deformation at time \(t+\Delta t\) resulting from thermal expansion and anelastic deformations (which may include enlargement deformations, cf. [§3.2]). So we have:

\(\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{o}}=\left[\alpha \left(t+\mathrm{\Delta }t\right)\left(T\left(t+\mathrm{\Delta }t\right)-{T}_{\text{ref}}\right)-\alpha \left(t\right)\left(T\left(t\right)-{T}_{\text{ref}}\right)\right]\mathrm{Id}+{\mathrm{\epsilon }}_{\mathrm{a}}\left(t+\mathrm{\Delta }t\right)-{\mathrm{\epsilon }}_{\mathrm{a}}\left(t\right)\)

where \(\mathrm{Id}\) is the 2nd order identity tensor in dimension 3.

We pose \(\mathrm{\Delta }\mathrm{\epsilon }=\mathrm{\Delta }{\mathrm{\epsilon }}_{\text{tot}}-\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{o}}\)

Since we assume here that \(\mu\) is constant, we have the following relationship between the deviators of \(\Delta \sigma\) and \(\Delta \varepsilon\):

\(\mathrm{\Delta }\stackrel{~}{\mathrm{\sigma }}=2\mu \left(\mathrm{\Delta }\stackrel{~}{\mathrm{\epsilon }}-\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}\right)\) eq 4.1-1

However, the law of flow is written, implicitly:

\(\frac{\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}}{\mathrm{\Delta }t}=\frac{3}{2}g\left({\sigma }_{\text{eq}},{\lambda }^{-}+{\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}},T\right)\frac{\stackrel{~}{\mathrm{\sigma }}}{{\sigma }_{\text{eq}}}\) eq 4.1-2

So, by eliminating \(\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}\) between [éq 4.1-1] and [éq 4.1-2], we have:

\(\begin{array}{c}2\mu \mathrm{\Delta }\stackrel{~}{\mathrm{\epsilon }}=\mathrm{\Delta }\stackrel{~}{\mathrm{\sigma }}+3\mu \mathrm{\Delta }t\mathrm{.}g\left({\sigma }_{\text{eq}},{\lambda }^{-}+{\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}},T\right)\frac{\stackrel{~}{\mathrm{\sigma }}}{{\sigma }_{\text{eq}}}\\ \left({\stackrel{~}{\mathrm{\sigma }}}^{-}+2\mu \mathrm{\Delta }\stackrel{~}{\mathrm{\epsilon }}\right)=\left(1+3\mu \mathrm{\Delta }t\frac{g\left({\sigma }_{\text{eq}},{\lambda }^{-}+{\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}},T\right)}{{\sigma }_{\text{eq}}}\right)\stackrel{~}{\mathrm{\sigma }}\end{array}\) eq 4.1-3

By setting \({\stackrel{~}{\mathrm{\sigma }}}^{e}={\stackrel{~}{\mathrm{\sigma }}}^{-}+2\mu \mathrm{\Delta }\stackrel{~}{\mathrm{\epsilon }}\), we therefore have:

\({\sigma }_{\text{eq}}^{e}=3\mu \mathrm{\Delta }t\mathrm{.}g\left({\sigma }_{\text{eq}},{\lambda }^{-}+{\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}},T\right)+{\sigma }_{\text{eq}}\) eq 4.1-4

However, according to [éq 4.1-2]:

\({\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}}=\mathrm{\Delta }t\mathrm{.}g\left({\sigma }_{\text{eq}},{\lambda }^{-}+{\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}},T\right)\)

From where:

\(\begin{array}{c}{\sigma }_{\text{eq}}^{e}=3\mu {\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}}+{\sigma }_{\text{eq}}\\ {\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}}=\frac{1}{3\mu }\left({\sigma }_{\text{eq}}^{e}-{\sigma }_{\text{eq}}\right)\end{array}\)

Substituting this last expression into [éq 4.1-4], we have:

\({\sigma }_{\text{eq}}^{e}=3\mu \mathrm{\Delta }t\mathrm{.}g\left({\sigma }_{\text{eq}},{\lambda }^{-}+\frac{1}{3\mu }\left({\sigma }_{\text{eq}}^{e}-{\sigma }_{\text{eq}}\right),T\right)+{\sigma }_{\text{eq}}\)

If we ask, \({\sigma }_{\text{eq}}^{e},{\lambda }^{-},T\) and \(\Delta t\) being known:

\(f\left(x\right)=3\mu \mathrm{\Delta }t\mathrm{.}g\left(x,{\lambda }^{-}+\frac{1}{3\mu }\left({\sigma }_{\text{eq}}^{e}-x\right),T\right)+x-{\sigma }_{\text{eq}}^{e}\)

we can then calculate the quantity \({\sigma }_{\text{eq}}={\left({\sigma }^{-}+\mathrm{\Delta }\sigma \right)}_{\text{eq}}\) as being the solution of the scalar equation: \(f\left(x\right)=0\) or \(x={\sigma }_{\text{eq}}\), convention adopted for the following paragraphs.

In the case of a Lemaître law with threshold, LEMA_IRRA_SEUIL, the preceding equations are only useful once the threshold has been crossed. In fact, below the threshold, the behavior is elastic.

The threshold is discretized implicitly:

\(D({\sigma }^{-}+\mathrm{\Delta }\sigma )=\frac{1}{S}\underset{0}{\overset{t}{\int }}{\left({\sigma }^{-}+\mathrm{\Delta }\sigma \right)}_{\text{eq}}(u)\text{du}\)

In the same way as for the integration of the Von-Mises elasto-plastic behavior laws, two cases are then distinguished.

\(D({\sigma }^{-}+\mathrm{\Delta }\sigma )\le 1\) then \(\phantom{\rule{0.5em}{0ex}}g\left({\sigma }^{-}+\mathrm{\Delta }\sigma ,\lambda ,T\right)=0\)

\(D({\sigma }^{-}+\mathrm{\Delta }\sigma )>1\) then \(\phantom{\rule{0.5em}{0ex}}g\left({\sigma }^{-}+\mathrm{\Delta }\sigma ,\lambda ,T\right)=A\text{.}\left(\frac{2}{\sqrt{3}}{\sigma }_{\text{eq}}\right)\phi\)

This results from the equation above:

\(\phantom{\rule{0.5em}{0ex}}g\left({\sigma }^{-}+\mathrm{\Delta }\sigma ,\lambda ,T\right)\ne A\text{.}\left(\frac{2}{\sqrt{3}}{\sigma }_{\text{eq}}\right)\phi\) involves \(D({\sigma }^{-}+\mathrm{\Delta }\sigma )\le 1\)

But \(g\) can only take the value \(0\) or \(A\text{.}(\frac{2}{\sqrt{3}}{\sigma }_{\text{eq}})\varphi\) so

\(g({\sigma }^{-}+\Delta \sigma ,\lambda ,T)\ne A\text{.}(\frac{2}{\sqrt{3}}{\sigma }_{\text{eq}})\varphi\) involves \(D({\sigma }^{-}+\Delta \sigma )\le 1\) and \(\Delta \sigma =A\text{.}\Delta \varepsilon\)

Let’s be \(g({\sigma }^{-}+\Delta \sigma ,\lambda ,T)=A\text{.}(\frac{2}{\sqrt{3}}{\sigma }_{\text{eq}})\varphi\) then \(D({\sigma }^{-}+A\text{.}\Delta \varepsilon )>1\)

The threshold crossing criterion can therefore be written \(D({\sigma }^{-}+A\text{.}\Delta \varepsilon )>1\)

Gold \(D({\sigma }^{-}+A\text{.}\Delta \varepsilon )=\frac{1}{S}\underset{0}{\overset{t}{\int }}{({\sigma }^{-}+A\text{.}\Delta \varepsilon )}_{\text{eq}}(u)\text{du}\)

By discretizing time, we then have:

\(D({\sigma }^{-}+A\text{.}\Delta \varepsilon )=\frac{1}{S}\underset{0}{\overset{{t}^{-}}{\int }}{({\sigma }^{-}+A\text{.}\Delta \varepsilon )}_{\text{eq}}(u)\text{du}+\frac{1}{S}\underset{{t}^{-}}{\overset{t}{\int }}{({\sigma }^{-}+A\text{.}\Delta \varepsilon )}_{\text{eq}}(u)\text{du}\)

\(D({\sigma }^{-}+A\text{.}\Delta \varepsilon )=\frac{1}{S}({D}^{-}S+\frac{{\sigma }_{\text{eq}}^{-}+{({\sigma }^{-}+A\text{.}\Delta \varepsilon )}_{\text{eq}}}{2}(t-{t}^{-}))\)

4.2. Solving the scalar equation: routine principle ZEROF2

It is easy to show that, if the conditions required in paragraph [§3] on the characteristics of materials are met, the function \(f\) is strictly increasing and the equation \(f(x)=0\) admits a unique solution.

If \({\sigma }_{\text{eq}}^{e}=0\), then the solution is \(x=0\). If not, we have: \(f(0)=-{\sigma }_{\text{eq}}^{e}<0\)

The problem therefore consists in finding the solution of equation \(f(x)=0\) for any function \(f\), knowing that this solution exists, that \(f(0)<0\) and that \(f\) is strictly increasing.

The algorithm adopted in ZEROF2 is as follows:

• we start from \({a}_{0}=0\) and \({b}_{0}={x}_{\mathrm{ap}}\) where \({x}_{\mathrm{ap}}\) is an approximation of the solution. If necessary (that is to say if \(f({b}_{0})<0\)), we reduce ourselves by the secant method (\({z}_{n}=\frac{{a}_{n}f({b}_{n})-{b}_{n}f({a}_{n})}{f({b}_{n})-f({a}_{n})}\) then \({a}_{n+1}={b}_{n}\) and \({b}_{n+1}={z}_{n}\)) in one or more iterations in case \(f(a)<0\) and \(f(b)>0\):

(In the case of the figure above, this first sentence was done in one iteration: \({a}_{1}={b}_{0}\) and \(f({b}_{1})>0\)).

• we calculate \({N}_{d}\) = integer part \((\sqrt{{N}_{\text{max}}})\) where \({N}_{\text{max}}\) is the maximum number of iterations we have given ourselves. The equation is then solved by the secant method, using however the dichotomy method each time \(n\) is a multiple of \({N}_{d}\):

If \({N}_{d}\) divides \(n\)

\({z}_{n}=\frac{{a}_{n}+{b}_{n}}{2}\)

If not

\({z}_{n}=\frac{{a}_{n}f({b}_{n})-{b}_{n}f({a}_{n})}{f({b}_{n})-f({a}_{n})}\)

Finsi

\(n=n+1\)

So \(\mid f(z)\mid >\varepsilon\)

So \(f(z)<0\)

\({a}_{n+1}={z}_{n}\) \({b}_{n+1}={b}_{n}\)

If not

\({a}_{n+1}={a}_{n}\) \({b}_{n+1}={z}_{n}\)

Finsi

go to 1)

If not

The solution is: \(x={z}_{n}\to \text{FIN}\)

Finsi

This second part of the algorithm makes it possible to process in a reasonable number of iterations cases where \(f\) is very strongly non-linear, while the secant method would have converged too slowly. These cases of strong non-linearity are encountered in particular with the law of LEMAITRE, for values of \(\frac{n}{m}\) that are large.

4.3. Calculation of the stress at the end of the current time step

According to [éq 4.1-3], if \(x\) is the solution to the scalar equation, by asking:

\(b(x,{\sigma }_{\text{eq}}^{e})=\frac{1}{1+3\mu \Delta t\frac{g(x,{\lambda }^{-}+\frac{1}{3\mu }({\sigma }_{\text{eq}}^{e}-x),T)}{x}}=\frac{x}{{\sigma }_{\text{eq}}^{e}}\)

we have:

\(\stackrel{~}{\mathrm{\sigma }}=b\left(x,{\sigma }_{\text{eq}}^{e}\right){\stackrel{~}{\mathrm{\sigma }}}^{e}\) eq 4.3-1

In the case where \({\sigma }_{\text{eq}}^{e}=0\), which is equivalent according to the scalar equation to \(x=0\), \(b\) is extended by continuity. To do this, we ask \(y(x)={\lambda }^{-}+\frac{1}{3\mu }({\sigma }_{\text{eq}}^{e}-x)\) and \(G(x)=g(x\text{,y}(x),T)\). The derivative of \(G\) is expressed in terms of the partial derivatives of \(g\) at point \((x\text{,y}(x),T)\):

\(\text{G'}\left(x\right)=\frac{\partial g}{\partial x}\left(x\text{,y}\left(x\right),T\right)-\frac{1}{3\mu }\frac{\partial g}{\partial y}\left(x\text{,y}\left(x\right),T\right)\)

Extending \(b\) by continuity then gives:

\(b\left(\mathrm{0,0}\right)=\frac{1}{1+3\mu \mathrm{\Delta }t\text{G'}\left(0\right)}\)

and we have, always in the case where \({\sigma }_{\text{eq}}^{e}=0\), \(\tilde{\sigma }=0\).

Once we have calculated \(\stackrel{~}{\mathrm{\sigma }}\), we get \(\mathrm{\sigma }\) by the relationship (\(K\) is assumed to be constant here):

\(\mathrm{\sigma }={\mathrm{\sigma }}^{-}+\mathrm{\Delta }\mathrm{\sigma }=\stackrel{~}{\mathrm{\sigma }}+\left(\frac{1}{3}\text{Tr}\left({\mathrm{\sigma }}^{-}\right)+K\mathrm{.}\text{Tr}\left(\mathrm{\Delta }\mathrm{\epsilon }\right)\right)\mathrm{.}\mathrm{Id}\) eq 4.3-2

The internal variables are as follows: \(V1\) : cumulative plastic deformation, \(V2\): cumulative viscoelastic energy dissipation density.

4.4. Semi-implicit schema

With an implicit numerical scheme [éq 4.1-2], in the case, for example, where \(g\) does not depend on \(\lambda\), only the value of the constraint at the end of the time step intervenes when calculating \(\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}\). This can result in significant numerical errors if the constraint varies greatly over time (see [bib2]).

To remedy this and improve the resolution, the flow law is semi-implicitly discretized:

\(\frac{\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}}{\mathrm{\Delta }t}=\frac{3}{2}g\left({\left({\sigma }^{-}+\frac{\mathrm{\Delta }\sigma }{2}\right)}_{\text{eq}},{\lambda }^{-}+\frac{{\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}}}{2},{T}^{-}+\frac{\mathrm{\Delta }T}{2}\right)\frac{\left({\stackrel{~}{\mathrm{\sigma }}}^{-}+\frac{\mathrm{\Delta }\stackrel{~}{\mathrm{\sigma }}}{2}\right)}{{\left({\sigma }^{-}+\frac{\mathrm{\Delta }\sigma }{2}\right)}_{\text{eq}}}\) eq 4.4-1

To transform what was programmed previously in the most economical way (using the implicit formulation [éq 4.1-2]), simply divide each member of equation [éq 4.4-1] by 2:

\(\frac{\left(\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}/2\right)}{\mathrm{\Delta }t}=\frac{3}{2}\frac{g\left({\left({\sigma }^{-}+\frac{\mathrm{\Delta }\sigma }{2}\right)}_{\text{eq}},{\lambda }^{-}+\frac{{\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}}}{2},{T}^{-}+\frac{\mathrm{\Delta }T}{2}\right)}{2}\frac{\left({\stackrel{~}{\mathrm{\sigma }}}^{-}+\frac{\mathrm{\Delta }\stackrel{~}{\mathrm{\sigma }}}{2}\right)}{{\left({\sigma }^{-}+\frac{\mathrm{\Delta }\sigma }{2}\right)}_{\text{eq}}}\)

and to do the same thing with the relationship [éq 4.1-1]:

\(\frac{\mathrm{\Delta }\stackrel{~}{\mathrm{\sigma }}}{2}=2\mu \left(\frac{\mathrm{\Delta }\stackrel{~}{\mathrm{\epsilon }}}{2}-\frac{\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}}{2}\right)\)

We note that this system is of the same form as that constituted by the equations [éq 4.1-1] and [éq 4.1-2], the data being \(\frac{\mathrm{\Delta }\mathrm{\epsilon }}{2}\) instead of \(\mathrm{\Delta }\mathrm{\epsilon }\), the unknowns being respectively \(\frac{\mathrm{\Delta }\mathrm{\sigma }}{2}\) and \(\frac{\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}}{2}\) instead of \(\mathrm{\Delta }\mathrm{\sigma }\) and \(\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}\) and the function \(\frac{g}{2}\) replacing the function \(g\).

We can therefore use the resolution of paragraphs [§4.1] to [§4.3] as well as the corresponding algorithm by introducing \(\frac{\mathrm{\Delta }\mathrm{\epsilon }}{2}\) and dividing the function \(g\) by \(2\). It then remains to multiply the results \(\frac{\mathrm{\Delta }\mathrm{\sigma }}{2}\) and \(\frac{\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}}{2}\) by \(2\) to obtain the stress and viscous deformation increments calculated by the semi-implicit diagram (the \(\mathrm{\Delta }\mathrm{\sigma }\) and the \(\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}\) of equation [éq 4.4-1]).

It will be noted that the calculation of the tangent operator is not affected by this modification of the numerical diagram. Indeed, we obviously have:

\(\frac{\partial \mathrm{\Delta }\mathrm{\sigma }}{\partial \mathrm{\Delta }\mathrm{\epsilon }}=\frac{\partial \left(\frac{\mathrm{\Delta }\mathrm{\sigma }}{2}\right)}{\partial \left(\frac{\mathrm{\Delta }\mathrm{\epsilon }}{2}\right)}\)

4.5. Taking into account the variation of elastic coefficients as a function of temperature

We have, if \(\mathrm{A}\) is the elasticity tensor:

\(\mathrm{\Delta }\mathrm{\epsilon }=\mathrm{\Delta }{\mathrm{\epsilon }}_{\mathrm{v}}+\mathrm{\Delta }\left({\mathrm{A}}^{-1}\mathrm{.}\mathrm{\sigma }\right)\)

with:

\(\mathrm{\Delta }\left({\mathrm{A}}^{-1}\mathrm{.}\mathrm{\sigma }\right)={\mathrm{A}}^{-1}\left({T}^{-}+\mathrm{\Delta }T\right)\mathrm{.}\left({\mathrm{\sigma }}^{-}+\mathrm{\Delta }\mathrm{\sigma }\right)-{\mathrm{A}}^{-1}\left({T}^{-}\right)\mathrm{.}{\mathrm{\sigma }}^{-}\)

This is translated into the equations in [§4.4] by:

\(2\mu \left(\frac{\mathrm{\Delta }\stackrel{~}{\mathrm{\epsilon }}}{2}\right)-\left({\stackrel{~}{\mathrm{\sigma }}}^{-}+\frac{\mathrm{\Delta }\stackrel{~}{\mathrm{\sigma }}}{2}\right)=3\mu \mathrm{\Delta }t\frac{g\left({\left({\sigma }^{-}+\frac{\mathrm{\Delta }\sigma }{2}\right)}_{\text{eq}},{\lambda }^{-}+\frac{{\left(\mathrm{\Delta }{\epsilon }_{v}\right)}_{\text{eq}}}{2},{T}^{-}+\frac{\mathrm{\Delta }T}{2}\right)}{2}\frac{\left({\stackrel{~}{\mathrm{\sigma }}}^{-}+\frac{\mathrm{\Delta }\mathrm{\sigma }}{2}\right)}{{\left({\sigma }^{-}+\frac{\mathrm{\Delta }\sigma }{2}\right)}_{\text{eq}}}-{\stackrel{~}{\mathrm{\sigma }}}^{-}\mathrm{.}\left(\frac{2{\mu }^{-}+2\mu }{4{\mu }^{-}}\right)\)

By posing:

\({\stackrel{~}{\mathrm{\sigma }}}^{e}=\left(\frac{2{\mu }^{-}+2\mu }{4{\mu }^{-}}\right){\stackrel{~}{\mathrm{\sigma }}}^{-}+2\mu \left(\frac{\mathrm{\Delta }\stackrel{~}{\mathrm{\epsilon }}}{2}\right)\)

and

\(\text{Tr}\left({\mathrm{\sigma }}^{\mathrm{e}}\right)=\left(\frac{3{K}^{-}+3K}{6{K}^{-}}\right)\text{Tr}\left({\mathrm{\sigma }}^{-}\right)+3K\mathrm{.}\text{Tr}\left(\frac{\mathrm{\Delta }\mathrm{\epsilon }}{2}\right)\)

we are exactly back to the previous case [§4.4].