Digital integration
=====================

Implicit solving behavior equations in MFront
----------------------------------------------------------------

The modified Cam-Clay behavior model is integrated implicitly (directive :math:`\texttt{@DSL Implicit}`) via the MFront [HMPS15] _ tool. This integration intervenes on the total deformation increment imposed :math:`\Delta \boldsymbol{\epsilon}` at the current moment :math:`t_{n+1}` (with :math:`\boldsymbol{\epsilon}_{n+1}=\boldsymbol{\epsilon}_n+\Delta \boldsymbol{\epsilon}`) at each integration point, knowing the state of the internal variables and of the stress tensor at the previous instant :math:`t_n`. Elastic prediction assumes the increment of the total elastic deformation tensor (block :math:`\texttt{@Predictor}`). If the plasticity criterion :math:`f_{n+1}` presented :eq:`equations_optimalite` -2 is not verified, the plastic correction must ensure a solution positioned on the surface of the reversibility domain (block :math:`\texttt{@Integrator}`). In this situation, the residue of the system of equations to be solved during integration is noted :math:`\boldsymbol{r}` and the vector of unknowns is noted :math:`\boldsymbol{x}`. They
are written as:

.. math::
\ begin {Bmatrix}\ boldsymbol {r} _1\\ r_2\\ r_2\\ r_3\ end {Bmatrix} =\ begin {Bmatrix} -\ Delta\ boldsymbol {\ epsilon} +\ Delta\ boldsymbol {\ epsilon} _1\ _1} +\ Delta\ boldsymbol {\ epsilon} ^e +\ Delta\ boldsymbol {\ epsilon} ^e +\ Delta\ boldsymbol {\ epsilon} ^e +\ Delta\ boldsymbol {\ epsilon} ^e +\ Delta\ boldsymbol {\ epsilon} ^e +\ Delta\ boldsymbol {\ epsilon} ^e +\ Delta\ boldsymbol {\ epsilon} ^e +\ Delta\ boldsymbol {\ ep\ sigma} _ {n+1}}\\ -\ -\ Delta\ xi+\ Delta\ xi+\ Delta\ lambda\ cfrac {\ partial f_ {n+1}} {\ partial p_ {c, n+1}}}\\ f_ {n+1}}}\\ f_ {n+1} /K=0\ end {Bmatrix}},\ qquad 
  \ begin {Bmatrix}\ boldsymbol {x} _1\\ x_2\\ x_3\ end {Bmatrix} =\ begin {Bmatrix}\ Delta\ boldsymbol {\ epsilon} ^e\\\ Delta {\ epsilon} ^e\\\ Delta\ x_2\\ x_3\ end {Bmatrix}
  :label: non_lineary_system
  
In the previous system, :math:`\boldsymbol{r}_1=\boldsymbol{0}` satisfies the additive decomposition of the total deformation tensor and the normal flow law on the plastic deformation tensor. The second :math:`r_2 =0` reflects compliance with the law of normality on the scalar work hardening variable. Equation :math:`r_3 =0` expresses that the solution is located on the surface of the reversibility domain at time :math:`t_{n+1}` (the equation is dimensioned by the compressibility module :math:`K`). The iterative Newton-Raphson method is used to solve this system of equations :math:`\boldsymbol{r}=\boldsymbol{0}`, this one being non-linear. To do this, we use the calculation of the Jacobian matrix :math:`J` at each iteration of the schema, expressed as the following block matrix:

.. math::
 J =\ begin {Bmatrix}\ cfrac {\ partial\ boldsymbol {r} _1} {\ partial\ Delta\ boldsymbol {\ epsilon} ^e} &\ cfrac {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1} {\ partial\ boldsymbol {r} _1 da}\\\ cfrac {\ partial r_2} {\ partial\ partial r_2} {\ partial\ delta r_2} {\ partial\ delta\ xi} {\ partial\ delta\ xi} &\ cfrac {\ partial r_2} {\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ delta r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial r_2} {\ partial\ partial rsymbol {\ epsilon} ^e} &\ cfrac {\ partial r_3} {\ partial\ Delta\ xi} &\ cfrac {\ partial r_3} {\ partial\ Delta\ lambda}\ end {Bmatrix}
 :label: Jacobian
 
The obtaining of each block is detailed in the appendix by analytical calculations. The value of the Newton-Raphson stop criterion is taken to :math:`10^{-14}` (this value is specified in the :math:`\texttt{@Epsilon}` directive).

Calculation of the algorithmic tangent operator
-------------------------------------------

The algorithmic tangent operator :math:`\mathbb{D}`, calculated once the behavioral equations have been solved (block :math:`\texttt{@TangentOperator}`), is obtained by a chain rule of derivatives:

.. math::
\ mathbb {D} =\ frac {\ partial\ Delta\ partial\ Delta\ boldsymbol {\ sigma}} {\ partial\ Delta\ boldsymbol {\ epsilon}} =\ mathbb {C} _ {n+1}:\ frac {\ partial\ Delta\ delta\ boldsymbol {\ epsilon}} ^e} {\ partial\ Delta\ boldsymbol {\ epsilon}} {\ partial\ Delta\ boldsymbol {\ epsilon}} {\ partial\ Delta\ boldsymbol {\ epsilon}} {\ partial\ Delta\ boldsymbol {\ epsilon}} 
  
The first term :math:`\mathbb{C}_{n+1}` refers to the fourth-order elasticity tensor. This results from the expression of the elastic energy potential :eq:`energie_elastique` as:

.. math::
\ mathbb {C} _ {n+1} =\ Upsilon (\ boldsymbol {\ epsilon} ^e_ {n+1})\ left (3K\ mathbb {J} +2\ mu\ mathbb {K} - 2\ kappa\ mathbb {K} - 2\ kappa\ mathbb {K} - 2\ kappa\ mu\ epsilon} - 2\ kappa\ mu\ epsilon} - 2\ kappa\ mu\ epsilon} - 2\ kappa\ mu\ epsilon} - 2\ kappa\ mu\ mu\ left (\ boldsymbol {I}\ otimes\ boldsymbol {\ epsilon} _ {d, n+1} {K} - 2\ kappa\ mu\ bb {K} - 2\ kappa\ mu\ mu\ left +\ boldsymbol {\ epsilon} _ {d, n+1} ^e\ otimes\ boldsymbol {I}\ right) +\ frac {\ left (2\ kappa\ mu\ right) ^2} {K} {K}\ boldsymbol {\ epsilon}\\ boldsymbol {\ epsilon} _ {d, n+1} _ {d, n+1} _ {d, n+1} ^e\ otimes\ boldsymbol {\ epsilon} _ {d, n+1} ^e\ right)
  
with :math:`\mathbb{J}` and :math:`\mathbb{K}` the fourth-order orthogonal projectors, such as for any symmetric second-order tensor :math:`\boldsymbol{a}`:

.. math::
\ mathbb {J}:\ boldsymbol {a} =\ frac {\ mathrm {tr} (\ boldsymbol {a})} {3}\ boldsymbol {I}\ qquad\ text {and}\ qquad\ mathbb {a} =\ boldsymbol {a})} {\ boldsymbol {a} -\ frac {\ mathrm {text {and}}\ qquad\ text {and}\ qquad\ mathbb {K}:\ boldsymbol {a} =\ boldsymbol {a} -\ frac {\ mathrm {tr} {tr} (boldsymbol {a}) {a})} {3}\ boldsymbol {I}
  
Note that :math:`\mathbb{C}_{n+1}` is not constant due to the nonlinearity of elasticity when :math:`\kappa>0`. 

The second term :math:`\cfrac{\partial \Delta\boldsymbol{\epsilon}^e}{\partial \Delta\boldsymbol{\epsilon}}` is written as:

.. math::
\ cfrac {\ partial\ Delta\ boldsymbol {\ boldsymbol {\ epsilon} ^e} {\ partial\ Delta\ boldsymbol {\ epsilon}} = J^ {-1}} _ {eel}

where :math:`J^{-1}_{eel}` represents the first upper left block of the inverse of the Jacobian matrix in :eq:`jacobienne`. The proof of this calculation can be found in reference [HMPS15] _.


Internal variables
------------------

The internal variables, whose increments constitute the unknowns of the non-linear system :eq:`systeme_non_lineaire` to be solved (directive :math:`\texttt{@StateVariable}`), are summarized in the following table.

.. list-table:: Variables internes du modèle.

 *-**Appellation**
   - **Definition**
   - **Symbol**
   - **2D components**
   - **3D components**
 *-*eel*
   - Tensor of elastic deformations
   - :math:`\boldsymbol{\epsilon}^e`
   - V1-V4
   - V1-V6
 *-*lam*
   - Cumulative plastic deformation
   - :math:`\lambda`
   - V5
   - V7
 *-*epv*
   - Volume plastic deformation
   - :math:`\xi`
   - V6
   - V8
   
A few additional auxiliary variables (directive :math:`\texttt{@AuxiliaryStateVariable}`) appear in the internal variables field in the sense of Code_Aster. These are summarized in the table below.

.. list-table:: Variables auxiliaires du modèle.

 *-**Appellation**
   - **Definition**
   - **Symbol**
   - **2D components**
   - **3D components**
 *-*oldyIeldSize*
   - Ray at the previous moment :math:`t_n` of the reversibility domain
   - :math:`R(\xi_n)`
   - V7
   - V9
 *-*yieldSizeGrowth*
   - Positive increase in the reversibility range over the :math:`[t_n,t_{n+1}]` interval
   - :math:`H_n \approx \frac{\langle R(\xi_{n+1})-R(\xi_n)\rangle_+}{\Delta\lambda}`
   - V8
   - V10
 *-*PlasticIndicator*
   - 0 if elastic, 1 if plastic 
   - :math:`Ip`
   - V9
   - V11
 *-*Dissipation*
   - Volume density of energy dissipated
   - :math:`\overline{D}=\int_0^{t_{n+1}}\left(\boldsymbol{\sigma}:\dot{\boldsymbol{\epsilon}}-\dot{\psi}\right)\mathrm{d}\tau`
   - V10
   - V12