5. Numerical formulation#
For the variational formulation, it is the same as that given in note [R5.03.21] and which refers to the law of large deformations behavior. 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. Integration of different behavioral relationships#
In the case of incremental behavior, keyword factor COMPORTEMENT, knowing the stress tensor \({\sigma }^{\text{-}}\), the internal variables \({r}_{k}^{-}\), the trace divided by three of the elastic deformation tensor \(\frac{1}{3}\text{tr}\stackrel{ˉ}{{b}^{e}}\), the displacements \({u}^{-}\) and \(\Delta u\), the displacements and, the temperatures \({T}^{-}\) and \(T\), and the proportions of different metallurgical phases \({Z}_{k}\), \(\Delta {Z}_{k}\), the aim is to determine \((\sigma ,{r}_{k},\frac{1}{3}\text{tr}\stackrel{ˉ}{{b}^{e}})\).
Since the displacements are known, the gradients of the transformation from \({O}_{0}\) to \({O}^{-}\), noted \({F}^{-}\), and from \({O}^{-}\) to \(\Omega (t)\), noted \(\Delta F\), are known.
We will then ask:
\(\mathrm{DA}=\sum _{i=1}^{4}{K}_{i}{F}_{i}^{\text{'}}\langle {\mathrm{DZ}}_{i}\rangle\), \(\Delta {G}_{\gamma }=\frac{\sum _{i=1}^{4}\langle -\Delta {Z}_{i}\rangle ({\theta }_{i\gamma }{r}_{i}^{-}-{r}_{\gamma }^{-})}{{Z}_{\gamma }}\) and
\(\Delta {G}_{i}=\frac{\langle \Delta {Z}_{i}\rangle ({\theta }_{\gamma i}{r}_{\gamma }^{-}-{r}_{i}^{-})}{{Z}_{i}}\) \((i=\mathrm{1,}4)\)
The implicit discretization of the law gives:
\(F=\Delta {\text{F F}}^{-}\)
\(J=\text{det}F\)
\(\stackrel{ˉ}{F}={J}^{-1/3}F\)
\({\stackrel{ˉ}{b}}^{e}={J}^{-2/3}{b}^{e}\)
\(J\sigma =\tau\)
\(\tilde{\tau }=\mu {\tilde{\stackrel{ˉ}{b}}}^{e}\)
\(\text{tr}\tau =\frac{\mathrm{3K}}{2}({J}^{2}-1)-\frac{\mathrm{9K}}{2}{\varepsilon }^{\text{th}}(J+\frac{1}{J})\)
\({\varepsilon }^{\text{th}}={Z}_{\gamma }\left[{\alpha }_{\gamma }(T-{T}_{\text{ref}})-(1-{Z}_{\gamma }^{r})\Delta {\varepsilon }_{f\gamma }^{T\text{ref}}\right]+(\sum _{i=1}^{4}{Z}_{i})\left[{\alpha }_{f}(T-{T}_{\text{ref}})+{Z}_{\gamma }^{r}\Delta {\varepsilon }_{f\gamma }^{T\text{ref}}\right]\)
\(f={\tau }_{\text{eq}}-(1-\stackrel{ˉ}{f}){R}_{\gamma }-\frac{\stackrel{ˉ}{f}}{Z}\sum _{i=1}^{4}{Z}_{i}{R}_{i}-{\sigma }_{y}\)
\({\stackrel{ˉ}{b}}^{e}=\stackrel{ˉ}{F}{G}^{p}{\stackrel{ˉ}{F}}^{T}=\stackrel{ˉ}{F}{G}^{p-}{\stackrel{ˉ}{F}}^{T}-\frac{\Delta p\text{tr}{\stackrel{ˉ}{b}}^{e}}{{\tau }_{\text{eq}}}\tilde{\tau }-\Delta A\text{tr}\stackrel{ˉ}{{b}^{e}}\tilde{\tau }\)
If \({Z}_{\gamma }>0\) then \(\Delta {r}_{\gamma }=\Delta p+\Delta {G}_{\gamma }-\underset{\text{uniquement en viscosité}}{\underset{\underbrace{}}{\Delta t({\text{Cr}}_{\text{moy}}^{-}{)}^{m}}}\), otherwise \({r}_{\gamma }^{-}=0\) and \(\Delta {r}_{\gamma }=0\)
If \({Z}_{i}>0\), \(\Delta {r}_{i}=\Delta p+\Delta {G}_{i}-\underset{\text{uniquement en viscosité}}{\underset{\underbrace{}}{\Delta t({\text{Cr}}_{\text{moy}}^{-}{)}^{m}}}\), otherwise \({r}_{i}^{-}=0\) and \(\Delta {r}_{i}=0\)
In the resolution of this system, only the deviatoric constraint \(\tilde{\tau }\) is unknown because the trace of \(\tau\) is a function only of \(J\) (known).
We introduce \({\tau }^{\text{Tr}}\), the Kirchhoff tensor which results from an elastic prediction (\(\mathrm{Tr}\): trial, in English test):
\({\tilde{\tau }}^{\text{Tr}}=\mu {\tilde{\stackrel{ˉ}{b}}}^{\text{eTr}}\)
where
\({\stackrel{ˉ}{b}}^{\text{eTr}}=\stackrel{ˉ}{F}{G}^{p-}{\stackrel{ˉ}{F}}^{T}=\Delta \stackrel{ˉ}{F}{\stackrel{ˉ}{b}}^{e-}\Delta {\stackrel{ˉ}{F}}^{T}\), \(\Delta \stackrel{ˉ}{F}=(\Delta J{)}^{-1/3}\Delta F\), and \(\Delta J=\text{det}(\Delta F)\)
We obtain \({\stackrel{ˉ}{b}}^{e-}\) from the stresses \({\tau }^{-}\) by the thermoelastic stress-deformation relationship and from the trace of the elastic deformation tensor.
\({\stackrel{ˉ}{b}}^{e-}=\frac{{\tilde{\tau }}^{-}}{{\mu }^{-}}+\frac{1}{3}\text{tr}\stackrel{ˉ}{{b}^{e-}}\)
For the Kirchhoff tensor \(\tilde{\tau }\), we obtain:
\(\tilde{\tau }=\mu {\stackrel{ˉ}{\tilde{b}}}^{\text{eTr}}-\mu \Delta p\frac{\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}}{{\tau }_{\text{eq}}}\tilde{\tau }-\mu \Delta \text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}\tilde{\tau }\)
If \(f<0\), then we have \(\Delta p=0\) and:
\(\tilde{\tau }=\frac{{\tilde{\tau }}^{\text{Tr}}}{1+\mu \Delta \text{A tr}\stackrel{ˉ}{{b}^{\text{eTr}}}}\)
otherwise we get:
\(\text{tr}\stackrel{ˉ}{{b}^{e}}=\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}\)
\(\begin{array}{}1+\mu \Delta p\frac{{\text{tr}\stackrel{ˉ}{b}}^{\text{eTr}}}{{\tau }_{\text{eq}}}\\ \\ \tilde{\tau }\left[+\mu \Delta \text{tr}\stackrel{ˉ}{b}\text{eTr}\right]={\tilde{\tau }}^{\text{Tr}}\end{array}\)
By calculating the equivalent stress, we obtain the following scalar equation in \(\Delta p\):
\({\tau }_{\text{eq}}+\mu \Delta p\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}+\mu \Delta A{\tau }_{\text{eq}}\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}={\tau }_{\text{eq}}^{\text{Tr}}\)
Expression of \({\tau }_{\text{eq}}\):
In plasticity: \({\tau }_{\text{eq}}={\sigma }_{y}+{R}^{\text{'}}\mathrm{Dp}+D({r}^{-}\text{;T,Z})\)
with
\({R}^{\text{'}}=(1-\stackrel{ˉ}{f}){R}_{0\gamma }+\frac{\stackrel{ˉ}{f}}{Z}\sum _{i=1}^{4}{Z}_{i}{R}_{\mathrm{0i}}\)
and \(D({r}^{-}\text{;T,Z})=[1-\stackrel{ˉ}{f}]{R}_{\gamma }({r}_{\gamma }^{-}+\Delta {G}_{\gamma })+\frac{\stackrel{ˉ}{f}}{Z}\sum _{i=1}^{4}{Z}_{i}{R}_{i}({r}_{i}^{-}+\Delta {G}_{i})\)
In viscosity:
\({\tau }_{\text{eq}}={\sigma }_{y}+{R}^{\text{'}}\mathrm{\Delta p}+D({r}^{-}\text{;T,Z})+(1-\stackrel{ˉ}{f}(Z)){\eta }_{\gamma }(\Delta \text{p/}\mathrm{\Delta t}{)}^{1/{n}_{\gamma }}+\frac{\stackrel{ˉ}{f}}{Z}\sum _{i=1}^{4}{Z}_{i}{\eta }_{i}(\Delta \text{p/}\mathrm{\Delta t}{)}^{1/{n}_{i}}\)
with
\(\begin{array}{}D({r}^{-}\text{;T,Z})=[1-\stackrel{ˉ}{f}]{R}_{\gamma }({r}_{\gamma }^{-}+{\mathrm{DG}}_{\gamma }-\Delta t({\text{Cr}}_{\text{moy}}^{-}{)}^{m})\\ \text{+}\frac{\stackrel{ˉ}{f}}{Z}\sum _{i=1}^{4}{Z}_{i}{R}_{i}({r}_{i}^{-}+\Delta {G}_{i}-\Delta t({\text{Cr}}_{\text{moy}}^{-}{)}^{m})\end{array}\)
\(\Delta p\) check:
\((1-\stackrel{ˉ}{f}(Z)){\eta }_{\gamma }(\mathrm{\Delta p}/\mathrm{\Delta t}{)}^{1/{n}_{\gamma }}+\frac{\stackrel{ˉ}{f}}{Z}\sum _{i=1}^{4}{Z}_{i}{\eta }_{i}(\mathrm{\Delta p}/\mathrm{\Delta t}{)}^{1/{n}_{i}}=\frac{{\tau }_{\text{eq}}^{\text{Tr}}-\mathrm{\mu \Delta p}\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}}{1+\mathrm{\mu \Delta A}\text{'tr}\stackrel{ˉ}{{b}^{\text{eTr}}}}-D({r}^{-}\text{;T,Z})-{\sigma }_{y}-{R}^{\text{'}}\mathrm{\Delta p}\)
The resolution is done in Code_Aster by a secant method with search interval [bib4].
Note:
In the case of non-linear isotropic work hardening, the work-hardening slops \({R}_{\mathrm{0k}}\) and the workworking \({R}_{k}\) in the expressions of \({R}^{\text{'}}\) and \(D({r}^{-};T,Z)\) correspond to the variables \({r}_{k}\) taken at the instant \(t\) , i.e. \({r}_{k}={r}_{k}^{-}+\Delta {G}_{k}+\Delta p-\Delta t({\text{Cr}}_{\text{moy}}^{-}{)}^{m}\). Now, as we do not know the value of these variables \({r}_{k}\) a priory of time, we solve the equation in \(\Delta p\) by taking the slopes \({R}_{\mathrm{0k}}\) and the workings \({R}_{k}\) for the quantities \({r}_{k}^{-}+\Delta {G}_{k}-\Delta t({\text{Cr}}_{\text{moy}}^{-}{)}^{m}\). Once the equation has been solved in \(\Delta p\) , we check, for each phase, that we are in the correct interval when calculating the work hardening and the slope. Otherwise, for the phase or phases concerned, the following interval is taken and the equation is solved again in \(\mathrm{Dp}\) . This process is continued until the right interval is found for all phases.
We then find for the constraint deviator:
\(\tilde{\tau }=\frac{1}{1+\mu \Delta A\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}}\left[1-\mu \frac{\Delta p\text{tr}{\stackrel{ˉ}{b}}^{\text{eTr}}}{{\tau }_{\text{eq}}^{\text{Tr}}}\right]{\tilde{\tau }}^{\text{Tr}}\)
Once the cumulative plastic deformation, the stress tensor and the tangent matrix have been calculated, a correction is made on the trace of the elastic deformation tensor \({\stackrel{ˉ}{b}}^{e}\) to take into account the plastic incompressibility, which is not maintained with the simplification made on the flow law [éq 4.3.1]. This correction is carried out by using a relationship between the invariants of \({\stackrel{ˉ}{b}}^{e}\) and \({\tilde{\stackrel{ˉ}{b}}}^{e}\) and by exploiting the plastic incompressibility condition \({J}^{p}=1\) (or equivalently \(\text{det}\stackrel{ˉ}{{b}^{e}}=1\)). This relationship is written as:
\({x}^{3}-{\stackrel{ˉ}{J}}_{2}^{e}x-(1-{\stackrel{ˉ}{J}}_{3}^{e})=0\)
with \({\stackrel{ˉ}{J}}_{2}^{e}=\frac{1}{2}({\stackrel{ˉ}{\tilde{b}}}^{e}{)}_{\text{eq}}^{2}=\frac{({\tau }_{\text{eq}}{)}^{2}}{2(\mu {)}^{2}}\), \({\stackrel{ˉ}{J}}_{3}^{e}=\text{det}\stackrel{ˉ}{\tilde{{b}^{e}}}=\text{det}\frac{\tilde{\tau }}{\mu }\), and \(x=\frac{1}{3}\text{tr}\stackrel{ˉ}{{b}^{e}}\)
Solving this third-degree equation makes it possible to obtain \(\text{tr}\stackrel{ˉ}{{b}^{e}}\) and therefore the thermoelastic deformation \({\stackrel{ˉ}{b}}^{e-}\) at the following time step. In the case where this equation admits several solutions, we take the solution that is closest to the solution of the previous time step. This is also why we store in an internal variable \(\frac{1}{3}\text{tr}\stackrel{ˉ}{{b}^{e}}\).
Note:
In the case where transformation plasticity is not taken into account, the expressions obtained are the same by taking \(\Delta A=0\).
In the case where it is the restoration of work hardening that is neglected then we also have the same expressions but taking all \(\theta\) equal to 1.
5.2. Tangent matrix expression#
Here we only give the expressions of the tangent matrix (option FULL_MECA during Newton iterations, option RIGI_MECA_TANG for the first iteration). For the hypotheses concerning the metallurgical part, they are the same as those in the document [R4.04.02]. For the major deformations part, we will find in the appendix to [bib1], the details of the linearization of the law of behavior.
We ask:
\(J=\text{det}F\), \({J}^{-}=\text{det}{F}^{-}\), and \(\Delta J=\text{det}\Delta F\)
For option FULL_MECA, we have:
\(\begin{array}{}\stackrel{ˉ}{A}=\frac{\partial \sigma }{\partial \Delta F}=\frac{(\Delta J{)}^{-1/3}}{J}H-\frac{1}{3J\Delta J}(H\Delta \stackrel{ˉ}{F})\otimes B-\frac{{J}^{-}}{{J}^{2}}\tau \otimes B\\ \text{+}\frac{{J}^{-}}{J}\left[\text{KJ}-\frac{3}{2}K{\varepsilon }^{\text{th}}(1-{J}^{-2})\right]\text{Id}\otimes B\end{array}\)
where \(B\) is equal to:
\(\begin{array}{}{B}_{\text{11}}=\Delta {F}_{\text{22}}\Delta {F}_{\text{33}}-\Delta {F}_{\text{23}}\Delta {F}_{\text{32}}\\ {B}_{\text{22}}=\Delta {F}_{\text{11}}\Delta {F}_{\text{33}}-\Delta {F}_{\text{13}}\Delta {F}_{\text{31}}\\ {B}_{\text{33}}=\Delta {F}_{\text{11}}\Delta {F}_{\text{22}}-\Delta {F}_{\text{12}}\Delta {F}_{\text{21}}\\ {B}_{\text{12}}=\Delta {F}_{\text{31}}\Delta {F}_{\text{23}}-\Delta {F}_{\text{33}}\Delta {F}_{\text{21}}\\ {B}_{\text{21}}=\Delta {F}_{\text{13}}\Delta {F}_{\text{32}}-\Delta {F}_{\text{33}}\Delta {F}_{\text{12}}\\ {B}_{\text{13}}=\Delta {F}_{\text{21}}\Delta {F}_{\text{32}}-\Delta {F}_{\text{22}}\Delta {F}_{\text{31}}\\ {B}_{\text{31}}=\Delta {F}_{\text{12}}\Delta {F}_{\text{23}}-\Delta {F}_{\text{22}}\Delta {F}_{\text{13}}\\ {B}_{\text{23}}=\Delta {F}_{\text{31}}\Delta {F}_{\text{12}}-\Delta {F}_{\text{11}}\Delta {F}_{\text{32}}\\ {B}_{\text{32}}=\Delta {F}_{\text{13}}\Delta {F}_{\text{21}}-\Delta {F}_{\text{11}}\Delta {F}_{\text{23}}\end{array}\)
and where \(H\) and \(H\Delta \stackrel{ˉ}{F}\) are given by:
In the elastic case (\(f<0\)):
\({H}_{\text{ijkl}}=\frac{\mu }{(1+\mu \Delta A\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}})}({\delta }_{\text{ik}}{\stackrel{ˉ}{b}}_{\text{lp}}^{e-}\Delta {\stackrel{ˉ}{F}}_{\text{jp}}+\Delta {\stackrel{ˉ}{F}}_{\text{ip}}{\stackrel{ˉ}{b}}_{\text{pl}}^{e-}{\delta }_{\text{jk}}-\frac{2}{3}{\delta }_{\text{ij}}\Delta {\stackrel{ˉ}{F}}_{\text{kp}}{\stackrel{ˉ}{b}}_{\text{lp}}^{e-}-2\Delta A{\tilde{\tau }}_{\text{ij}}\Delta {\stackrel{ˉ}{F}}_{\text{kp}}{\stackrel{ˉ}{b}}_{\text{pl}}^{e-})\)
and
\(H\Delta \stackrel{ˉ}{F}=\frac{2\mu }{(1+\mu \Delta A\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}})}({\tilde{\stackrel{ˉ}{b}}}^{\text{eTr}}-\Delta A\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}\tilde{\tau })\)
If not in plastic or viscoplastic filler, we have:
\(\begin{array}{}{H}_{\text{ijkl}}=\frac{\mu }{a}({\delta }_{\text{ik}}{\stackrel{ˉ}{b}}_{\text{lp}}^{e-}\Delta {\stackrel{ˉ}{F}}_{\text{jp}}+\Delta {\stackrel{ˉ}{F}}_{\text{ip}}{\stackrel{ˉ}{b}}_{\text{pl}}^{e-}{\delta }_{\text{jk}})\\ -2\mu \left[\frac{{\delta }_{\text{ij}}}{\mathrm{3a}}+\frac{{\stackrel{ˉ}{R}}^{\text{'}}({\tau }_{\text{eq}}\Delta A+\Delta p){\tilde{\tau }}_{\text{ij}}}{{\tau }_{\text{eq}}({\stackrel{ˉ}{R}}^{\text{'}}+\mu \text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}(1+{\stackrel{ˉ}{R}}^{\text{'}}\Delta A))}\right]\Delta {\stackrel{ˉ}{F}}_{\text{kp}}{\stackrel{ˉ}{b}}_{\text{lp}}^{e-}\\ +\frac{3{\mu }^{2}\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}({\stackrel{ˉ}{R}}^{\text{'}}\Delta p-{\tau }_{\text{eq}})}{a{\tau }_{\text{eq}}^{3}({\stackrel{ˉ}{R}}^{\text{'}}+\mu \text{tr}{\stackrel{ˉ}{b}}^{\text{eTr}}(1+{\stackrel{ˉ}{R}}^{\text{'}}\Delta A))}{\tilde{\tau }}_{\text{ij}}{\tilde{\tau }}_{\text{kq}}\Delta {\stackrel{ˉ}{F}}_{\text{qp}}{\stackrel{ˉ}{b}}_{\text{lp}}^{e-}\end{array}\)
and
\(\begin{array}{}H\Delta \stackrel{ˉ}{F}=\frac{2\mu }{a}{\stackrel{ˉ}{b}}^{\text{eTr}}-2\mu \text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}\left[\frac{\text{Id}}{\mathrm{3a}}+\frac{{\stackrel{ˉ}{R}}^{\text{'}}({\tau }_{\text{eq}}\Delta A+\Delta p)\tilde{\tau }}{{\tau }_{\text{eq}}({\stackrel{ˉ}{R}}^{\text{'}}+\mu \text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}(1+{\stackrel{ˉ}{R}}^{\text{'}}\Delta A))}\right]\\ +\frac{3{\mu }^{2}\text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}({\stackrel{ˉ}{R}}^{\text{'}}\Delta p-{\tau }_{\text{eq}})}{a{\tau }_{\text{eq}}^{3}({\stackrel{ˉ}{R}}^{\text{'}}+\mu \text{tr}\stackrel{ˉ}{{b}^{\text{eTr}}}(1+{\stackrel{ˉ}{R}}^{\text{'}}\Delta A))}(\tilde{\tau }:{\stackrel{ˉ}{b}}^{\text{eTr}})\tilde{\tau }\end{array}\)
with
\({\stackrel{ˉ}{R}}^{\text{'}}=(1-\stackrel{ˉ}{f}){R}_{\mathrm{0\gamma }}+\frac{\stackrel{ˉ}{f}}{Z}\sum _{i=\mathrm{1,4}}{Z}_{i}{R}_{\mathrm{0i}}+\underset{\text{uniquement en viscosite}}{\underset{\underbrace{}}{(1-\stackrel{ˉ}{f}(Z)){\eta }_{\gamma }(\mathrm{\Delta p}/\mathrm{\Delta t}{)}^{(1-{n}_{\gamma })/{n}_{\gamma }}/{n}_{\gamma }\mathrm{\Delta t}+\frac{\stackrel{ˉ}{f}}{Z}\sum _{i=1}^{4}{Z}_{i}{\eta }_{i}(\mathrm{\Delta p}/\mathrm{\Delta t}{)}^{(1-{n}_{i})/{n}_{i}}/{n}_{i}\mathrm{\Delta t}}}\)
\(a=\frac{{\tau }_{\text{eq}}^{\text{Tr}}}{{\tau }_{\text{eq}}}\)
For option RIGI_MECA_TANG
for the plastic model: these are the same expressions as those given for FULL_MECA but with \(\Delta p=0\) and \(\Delta A=0\), all the variables and coefficients of the material being taken at time \({t}^{-}\). In particular, we will have \(\Delta \stackrel{ˉ}{F}=\text{Id}\).
for the viscous model: we only take the expressions of FULL_MECA in the elastic case with \(\Delta A=0\), all the variables being taken at the time \({t}^{-}\).