5. Variational formulation

Insofar as the constraints provided by the law of behavior are Eulerian, a variational formulation written on the current configuration (Eulerian) and not on the initial configuration is chosen, i.e.:

\(\underset{{F}_{\text{int}}\text{.}\delta v}{\underset{\underbrace{}}{\underset{\Omega }{\int }\sigma {\nabla }_{x}\delta vd\Omega }}=\underset{{F}_{\text{ext}}\text{.}\delta v}{\underset{\underbrace{}}{\underset{\Omega }{\int }\rho f\delta vd\Omega +\underset{\partial {\Omega }^{f}}{\int }{t}^{d}\delta v\text{dS}}}\) \(\delta v\) Kinematically admissible

Here we are only interested in the work of internal forces and their variation in the perspective of a resolution by Newton’s method. We will find in [bib 4] the demonstration of the expressions presented.

5.1. Case of the continuous medium

Here we rewrite the work of domestic efforts in index form, i.e.:

\({F}_{\text{int}}\text{.}\delta v=\underset{\Omega }{\int }{\sigma }_{\text{ij}}\frac{\partial \delta {v}_{i}}{\partial {x}_{j}}d\Omega\)

We also need to express the variation in internal forces in the current configuration \(\Omega\) i.e.:

\(\begin{array}{}\begin{array}{}\delta {F}_{\text{int}}\text{.}\delta u\text{.}\delta v=\underset{\Omega }{\int }\left[{\sigma }_{\text{ij}}\frac{\partial \delta {u}_{p}}{\partial {x}_{p}}-{\sigma }_{\text{ik}}\frac{\partial \delta {u}_{j}}{\partial {x}_{k}}\right]\frac{\partial \delta {v}_{i}}{\partial {x}_{j}}d\Omega \text{}\text{rigidité géométrique}\text{}\\ \text{+}\underset{\Omega }{\int }\left[\frac{\partial {\sigma }_{\text{ij}}}{\partial \Delta {F}_{\text{pq}}}\frac{\partial \delta {u}_{p}}{\partial {x}_{q}^{-}}\right]\frac{\partial \delta {v}_{i}}{\partial {x}_{j}}d\Omega \text{}\text{rigidité de comportement}\text{}\end{array}\\ \end{array}\)

where \({x}^{-}\) are the coordinates of a point on the \({\Omega }^{-}\) configuration.

5.2. Finite element discretization

We discretize \(u\) movements and \(v\) virtual movements by finite elements. The notations are as follows, adopting the convention of summation of repeated indices:

\({u}_{i}(x)={N}^{n}(x){U}_{i}^{n}\frac{\partial {u}_{i}}{\partial {x}_{j}}={D}_{j}^{n}(x){U}_{i}^{n}\frac{\partial {u}_{i}}{\partial {x}_{j}^{-}}={D}_{j}^{-n}(x){U}_{i}^{n}\)

where:

\({N}^{n}(x)\) is the form function associated with node \(n\)

\({U}_{i}^{n}\), the i component of the node displacement of node \(n\)

\({D}_{j}^{n}(x)\), the components of the gradient of the shape functions on the \(\Omega\) configuration

\({D}_{j}^{-n}(x)\), the components of the gradient of the shape functions on the \({O}^{-}\) configuration

For the vector of internal forces, we obtain:

\(({F}_{\text{int}}{)}_{i}^{n}=\underset{\Omega }{\int }{\sigma }_{\text{ij}}{D}_{j}^{n}d\Omega\)

and for the tangent matrix, which is not symmetric in principle:

\(\begin{array}{}{K}_{i\text{}p}^{nm}=\underset{\Omega }{\int }\left[{D}_{p}^{m}{\sigma }_{\text{ij}}{D}_{j}^{n}-{D}_{k}^{m}{\sigma }_{\text{ik}}{D}_{p}^{n}\right]d\Omega \\ \text{+}\underset{\Omega }{\int }\left[{D}_{q}^{-n}\frac{\partial {\sigma }_{\text{ij}}}{\partial \Delta {F}_{\text{pq}}}{D}_{j}^{n}\right]d\Omega \end{array}\)

In the case of two-dimensional modeling (plane deformation), the expressions of the internal force vector and the tangent matrix are identical except that the indices corresponding to the components vary from 1 to 2 only.

In the case of axisymmetric modeling, by numbering the axes in order \((r,z,\theta ),\) the vector of internal forces is written:

\(({F}_{\text{int}}^{\text{axi}}{)}_{\alpha }^{n}=\underset{\Omega }{\int }\left[{\sigma }_{\alpha \beta }{D}_{\beta }^{n}+{\sigma }_{\text{33}}\frac{{N}^{n}}{r}{\delta }_{\alpha 1}\right]d\Omega\), \(\alpha \in \left\{\mathrm{1,2}\right\}\), \(\beta \in \left\{\mathrm{1,2}\right\}\)

and the tangent matrix:

\(\left[{K}^{\text{axi}}\right]=\left[K\right]+\left[{K}^{\text{corr}}\right]\)

with:

\({\left[{K}_{{}^{(1)}}^{\text{corr}}\right]}_{1\beta }^{\mathrm{nm}}=\underset{\Omega }{\int }\frac{{N}^{n}}{r}{\sigma }_{\beta \gamma }{D}_{\gamma }^{m}d\Omega +\underset{\Omega }{\int }\frac{{N}^{n}}{{r}^{-}}\frac{\partial {\sigma }_{\beta \gamma }}{\partial \Delta {F}_{\text{33}}}{D}_{\gamma }^{m}d\Omega\)

\({\left[{K}_{{}^{(2)}}^{\text{corr}}\right]}_{\alpha 1}^{\mathrm{nm}}=\underset{\Omega }{\int }{D}_{\alpha }^{n}{\sigma }_{\text{33}}\frac{{N}^{m}}{r}d\Omega +\underset{\Omega }{\int }{{D}_{\gamma }^{\text{-}}}^{n}\frac{\partial {\sigma }_{\text{33}}}{\partial \Delta {F}_{\alpha \gamma }}\frac{{N}^{m}}{r}d\Omega\)

\({\left[{K}_{{}^{(3)}}^{\text{corr}}\right]}_{11}^{\mathrm{nm}}=\underset{\Omega }{\int }\frac{{N}^{n}}{{r}^{-}}\frac{\partial {\sigma }_{\text{33}}}{\partial \Delta {F}_{\text{33}}}\frac{{N}^{m}}{r}\)

From an algorithmic point of view, the tangent elementary matrix \(K\) is not a priory symmetric. The overall resolution will therefore be done by default with a non-symmetric solver. However, it is possible to symmetrize the global tangent matrix before resolution (keyword SOLVEUR), which saves calculation time but may degrade convergence.

5.3. Expression of the tangent matrix of behavior

One gives here the expression for the tangent matrix (option FULL_MECA during Newton iterations, option RIGI_MECA_TANG for the first iteration). This is obtained by linearizing the system of equations that govern the law of behavior. Here we give the final result of this linearization. Details of this calculation can be found in [bib4].

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}{}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\\ +\frac{{J}^{-}}{J}\left[\text{KJ}-\frac{3}{2}K\alpha (T-{T}_{\text{ref}})(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}\)

\(H\) and \(H\Delta \stackrel{ˉ}{F}\) are given by:

In the elastic case (\(f<0\)):

\({H}_{\text{ijkl}}=\mu ({\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\mu }{3}{\delta }_{\text{ij}}\Delta {\stackrel{ˉ}{F}}_{\text{kp}}{\stackrel{ˉ}{b}}_{\text{lp}}^{e-}\)

and

\(H\Delta \stackrel{ˉ}{F}=2\mu {\tilde{\stackrel{ˉ}{b}}}^{\text{eTr}}\)

In the case of plastic (\(f=0\)) or viscoplastic:

\(\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}})\\ \text{}-2\mu \left[\frac{{\delta }_{\text{ij}}}{\mathrm{3a}}+\frac{\stackrel{ˉ}{R}\Delta p{\tilde{\tau }}_{\text{ij}}}{{\tau }_{\text{eq}}(\stackrel{ˉ}{R}+\mu \text{tr}{\stackrel{ˉ}{b}}^{\text{eTr}})}\right]\Delta {\stackrel{ˉ}{F}}_{\text{kp}}{\stackrel{ˉ}{b}}_{\text{lp}}^{e-}\\ \text{}+\frac{3{\mu }^{2}\text{tr}{\stackrel{ˉ}{b}}^{\text{eTr}}(\stackrel{ˉ}{R}\Delta p-{\tau }_{\text{eq}})}{a{\tau }_{\text{eq}}^{3}(\stackrel{ˉ}{R}+\text{}\mathrm{mu}\text{tr}{\stackrel{ˉ}{b}}^{\text{eTr}})}{\tilde{\tau }}_{\text{ij}}{\tilde{\tau }}_{\text{kq}}\Delta {\stackrel{ˉ}{F}}_{\text{qp}}{\stackrel{ˉ}{b}}_{\text{lp}}^{e-}\end{array}\)

and

\(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}\Delta p\tilde{\tau }}{{\tau }_{\text{eq}}(\stackrel{ˉ}{R}+\mu \text{tr}{\stackrel{ˉ}{b}}^{\text{eTr}})}\right]+\frac{3{\mu }^{2}\text{tr}{\stackrel{ˉ}{b}}^{\text{eTr}}(\stackrel{ˉ}{R}\Delta p-{\tau }_{\text{eq}})}{a{\tau }_{\text{eq}}^{3}(\stackrel{ˉ}{R}+\mu \text{tr}{\stackrel{ˉ}{b}}^{\text{eTr}})}(\tilde{\tau }:{\stackrel{ˉ}{b}}^{\text{eTr}})\tilde{\tau }\)

Where \(a=\frac{{\tau }_{\text{eq}}^{\text{Tr}}}{{\tau }_{\text{eq}}}\)

and \(\stackrel{ˉ}{R}=\text{R'}(p)+\underset{\text{uniquement cas visqueux}}{\underset{\underbrace{}}{{\sigma }_{0}\times {(1+{(\frac{\Delta p}{{\dot{\varepsilon }}_{0}\Delta t})}^{2})}^{\frac{-1}{2}}\times \frac{1}{m{({\dot{\varepsilon }}_{0}\Delta t)}^{\frac{1}{m}}}\times {(\Delta p)}^{-1}}}\),

\({R}^{\text{'}}(p)\) being the derivative of isotropic work hardening compared to the cumulative plastic deformation \(p\).

• For option RIGI_MECA_TANG, these are the same expressions as those given for FULL_MECA but with \(\Delta p=0\) and with all the variables and coefficients of the material taken at time \({t}^{-}\) (in principle, in the viscous case, you should take the expressions of FULL_MECA in the elastic case, take the expressions of in the elastic case, all the variables being taken at the moment \({t}^{-}\)). In particular, we will have \(\Delta \stackrel{ˉ}{F}=\text{Id}\).