5. Digital integration#
5.1. Implicit solving behavior equations in MFront#
The modified Cam-Clay behavior model is integrated implicitly (directive \(\texttt{@DSL Implicit}\)) via the MFront [HMPS15] _ tool. This integration intervenes on the total deformation increment imposed \(\Delta \boldsymbol{\epsilon}\) at the current moment \(t_{n+1}\) (with \(\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 \(t_n\). Elastic prediction assumes the increment of the total elastic deformation tensor (block \(\texttt{@Predictor}\)). If the plasticity criterion \(f_{n+1}\) presented equations_optimalite -2 is not verified, the plastic correction must ensure a solution positioned on the surface of the reversibility domain (block \(\texttt{@Integrator}\)). In this situation, the residue of the system of equations to be solved during integration is noted \(\boldsymbol{r}\) and the vector of unknowns is noted \(\boldsymbol{x}\). They
are written as:
- begin {Bmatrix}boldsymbol {r} _1\ r_2\ r_2\ r_3end {Bmatrix} =begin {Bmatrix} -Deltaboldsymbol {epsilon} +Deltaboldsymbol {epsilon} _1_1} +Deltaboldsymbol {epsilon} ^e +Deltaboldsymbol {epsilon} ^e +Deltaboldsymbol {epsilon} ^e +Deltaboldsymbol {epsilon} ^e +Deltaboldsymbol {epsilon} ^e +Deltaboldsymbol {epsilon} ^e +Deltaboldsymbol {epsilon} ^e +Deltaboldsymbol {epsigma} _ {n+1}}\ --Deltaxi+Deltaxi+Deltalambdacfrac {partial f_ {n+1}} {partial p_ {c, n+1}}}\ f_ {n+1}}}\ f_ {n+1} /K=0end {Bmatrix}},qquad
begin {Bmatrix}boldsymbol {x} _1\ x_2\ x_3end {Bmatrix} =begin {Bmatrix}Deltaboldsymbol {epsilon} ^e\Delta {epsilon} ^e\Deltax_2\ x_3end {Bmatrix} :label: non_lineary_system
In the previous system, \(\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 \(r_2 =0\) reflects compliance with the law of normality on the scalar work hardening variable. Equation \(r_3 =0\) expresses that the solution is located on the surface of the reversibility domain at time \(t_{n+1}\) (the equation is dimensioned by the compressibility module \(K\)). The iterative Newton-Raphson method is used to solve this system of equations \(\boldsymbol{r}=\boldsymbol{0}\), this one being non-linear. To do this, we use the calculation of the Jacobian matrix \(J\) at each iteration of the schema, expressed as the following block matrix:
The obtaining of each block is detailed in the appendix by analytical calculations. The value of the Newton-Raphson stop criterion is taken to \(10^{-14}\) (this value is specified in the \(\texttt{@Epsilon}\) directive).
5.2. Calculation of the algorithmic tangent operator#
The algorithmic tangent operator \(\mathbb{D}\), calculated once the behavioral equations have been solved (block \(\texttt{@TangentOperator}\)), is obtained by a chain rule of derivatives:
mathbb {D} =frac {partialDeltapartialDeltaboldsymbol {sigma}} {partialDeltaboldsymbol {epsilon}} =mathbb {C} _ {n+1}:frac {partialDeltadeltaboldsymbol {epsilon}} ^e} {partialDeltaboldsymbol {epsilon}} {partialDeltaboldsymbol {epsilon}} {partialDeltaboldsymbol {epsilon}} {partialDeltaboldsymbol {epsilon}}
The first term \(\mathbb{C}_{n+1}\) refers to the fourth-order elasticity tensor. This results from the expression of the elastic energy potential energie_elastique as:
mathbb {C} _ {n+1} =Upsilon (boldsymbol {epsilon} ^e_ {n+1})left (3Kmathbb {J} +2mumathbb {K} - 2kappamathbb {K} - 2kappamathbb {K} - 2kappamuepsilon} - 2kappamuepsilon} - 2kappamuepsilon} - 2kappamuepsilon} - 2kappamumuleft (boldsymbol {I}otimesboldsymbol {epsilon} _ {d, n+1} {K} - 2kappamubb {K} - 2kappamumuleft +boldsymbol {epsilon} _ {d, n+1} ^eotimesboldsymbol {I}right) +frac {left (2kappamuright) ^2} {K} {K}boldsymbol {epsilon}\ boldsymbol {epsilon} _ {d, n+1} _ {d, n+1} _ {d, n+1} ^eotimesboldsymbol {epsilon} _ {d, n+1} ^eright)
with \(\mathbb{J}\) and \(\mathbb{K}\) the fourth-order orthogonal projectors, such as for any symmetric second-order tensor \(\boldsymbol{a}\):
mathbb {J}:boldsymbol {a} =frac {mathrm {tr} (boldsymbol {a})} {3}boldsymbol {I}qquadtext {and}qquadmathbb {a} =boldsymbol {a})} {boldsymbol {a} -frac {mathrm {text {and}}qquadtext {and}qquadmathbb {K}:boldsymbol {a} =boldsymbol {a} -frac {mathrm {tr} {tr} (boldsymbol {a}) {a})} {3}boldsymbol {I}
Note that \(\mathbb{C}_{n+1}\) is not constant due to the nonlinearity of elasticity when \(\kappa>0\).
The second term \(\cfrac{\partial \Delta\boldsymbol{\epsilon}^e}{\partial \Delta\boldsymbol{\epsilon}}\) is written as:
cfrac {partialDeltaboldsymbol {boldsymbol {epsilon} ^e} {partialDeltaboldsymbol {epsilon}} = J^ {-1}} _ {eel}
where \(J^{-1}_{eel}\) represents the first upper left block of the inverse of the Jacobian matrix in jacobienne. The proof of this calculation can be found in reference [HMPS15] _.
5.3. Internal variables#
The internal variables, whose increments constitute the unknowns of the non-linear system systeme_non_lineaire to be solved (directive \(\texttt{@StateVariable}\)), are summarized in the following table.
A few additional auxiliary variables (directive \(\texttt{@AuxiliaryStateVariable}\)) appear in the internal variables field in the sense of Code_Aster. These are summarized in the table below.