3. Digital integration

3.1. Discretized equations

The time discretization of behavioral equations is based on an implicit Euler schema, i.e. the various variables of the problem are expressed at the final moment of the time step. The expression of the stress at the end of the time step thus obeys (), the deformation and the damage being also expressed at the end of the time step. There is still the delicate part which consists in calculating the evolution of the damage.

Integrating the damage evolution equation is a non-linear scalar problem. It results from the temporal discretization of the Kuhn and Tucker conditions (), by indicating by \(n\) the quantities at the beginning of the time step and by omitting any index for the one at the end of the time step:

(3.1)\[ b- {b} _ {n} ≥0\ text {;}\ text {;} Q (\ mathrm {\ epsilon}) -\ mathrm {\ varphi} (b) ≤0\ text {;} (b- {b} _ {n})\ text {;} 0\ text {;};} (b) _ {n})\ left [Q (\ mathrm {\ epsilon}) -\ mathrm {\ varphi} (b)\ right] =0\]

It is then a question of determining the value of \(b\), \(\mathrm{\epsilon }\) being known; this is called the integration of the law of behavior.

It should also be noted that deformation \(\mathrm{\epsilon }\) corresponds to the mechanical part of the deformations, that is to say the total deformation from which the shrinkage deformations, in particular thermal deformations (expansion), have been removed. In the terminology specific to the Aster code (Code_Aster), it is a law expressed in mechanical deformation.

3.2. Solving the algorithmic evolution equation

To treat the alternative (), we use the method described in [Lorentz, 2020] where we also demonstrated the existence and the uniqueness of the \(b\) solution. We start by evaluating function \(Q(\mathrm{\epsilon })\); this involves solving the scalar equation (). The reader is referred to document [R5.03.28] for the resolution method, the function \(Q\) being equal to one factor less than the function \(\mathrm{\Gamma }\) defined in this fascicle.

We then proceed with a bungee pull assuming that \(b={b}_{n}\). If \(Q(\mathrm{\epsilon })-\mathrm{\varphi }({b}_{n})⩽0\), then it’s the solution; it’s about the elastic regime. If not, we look at the ultimate value \(b=1\). If \(Q(\mathrm{\epsilon })-\mathrm{\varphi }(1)⩾0\), then it’s the solution; it’s the saturated diet. Finally, the last case, the damaging regime, solution \(b\) is the root of the following equation, between \({b}_{n}\) and 1:

(3.2)\[ \ mathrm {\ varphi} (b)\ text {-} Q (\ mathrm {\ epsilon})\ text {=} 0\]

The member on the left is ascending with \(b\) and has an opposite sign at both terminals. The resolution is achieved by a Newton method with controlled terminals.

The precision with which the two scalar equations are solved, that for evaluating \(Q(\mathrm{\epsilon })\) and (), is indirectly set by the user, so that the error made on the stress is less than the tensile limit \({\sigma }_{c}\) multiplied by the parameter RESI_INTE_RELA of STAT_NON_LINE. We can refer to [Lorentz,2020] for the calculation of the expected precisions for these two scalar equations.

3.3. Tangent matrix

The tangent matrix is obtained by varying the stress in () with respect to the deformation. To do this, we start by noting \({E}_{\text{+}}(\mathrm{\epsilon })\) and \({E}_{\text{-}}(\mathrm{\epsilon })\) the parts positive and negative of the elasticity operator:

(3.3)\[ {E} _ {\ text {+}} (\ mathrm {\ epsilon}) (\ mathrm {\ epsilon}))\ text {=}\ mathrm {\ epsilon}} (\ mathrm {\ epsilon}) ⟩} __} _ {\ text {+}} _ {\ text {+}}} +2\ mathrm {\ mu} {+}}\ text {;} {E} _ {\ text {-}}} (\ mathrm {\ epsilon})\ text {\ epsilon})\ text {=}\ mathrm {\ lambda} {⇒\ mathit {tr}} (\ mathrm {\ epsilon}} (\ mathrm {\ epsilon}) ⟩}}} _ {\ text {-}} +2\ mathrm {\ mu} {\ mathrm {\ tr} (\ tr}} (\ mathrm {\ epsilon}) (\ mathrm {\ epsilon})} ⟩} _ {\ text {-}}\]

By differentiating (), it then comes:

(3.4)\[ \ mathrm {\ delta}\ mathrm {\ sigma}\ text {\ sigma}\ text {=}\ left [B (b)\ frac {{\ mathit {dE}}}} {d\ mathrm {\ epsilon}}\ text {\ epsilon}} (\ mathrm {\ epsilon}})\ mathrm {\ epsilon}} (\ mathrm {\ epsilon}})\ mathrm {\ epsilon}} (\ mathrm {\ epsilon}})\ mathrm {\ epsilon}} (\ mathrm {\ epsilon}})\ mathrm {\ epsilon}} (\ mathrm {\ epsilon}})\ mathrm {\ epsilon}} {-}}} {d\ mathrm {\ epsilon}}} (\ mathrm {\ epsilon}}} (\ mathrm {\ epsilon}}} (\ mathrm {\ epsilon}\ text {\ epsilon}}\ text {\ epsilon}}\ text {\ epsilon})\ mathrm {\ epsilon}\ mathrm {\ epsilon}\ mathrm {\ epsilon}\ text {\ epsilon}\ text {\ epsilon}\ text {\ epsilon}\ text {\ epsilon}\ text {+}}} (\ mathrm {\ epsilon})\ frac {\ epsilon}\ text {\ epsilon}\ text {\ epsilon}\ text}} {\ mathit {db}}\ mathrm {\ delta} b\]

In elastic or saturated regimes, the last term is zero (no change in \(b\)). In a damaging regime, the variation in \(b\) is obtained by differentiating equation (). This leads to:

(3.5)\[ \ mathrm {\ delta} b\ text {=} {=} {\ left (\ frac {d\ mathrm {\ varphi}} {\ mathit {db}} (b)\ right)} ^ {-1}\ frac {\ mathit {dQ}} {\ frac {\ mathit {dQ}} {\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\ mathit {dQ}} {d\]

The tricky part then is to derive \(\mathrm{\Gamma }\), \({E}_{\text{+}}\), and \({E}_{\text{-}}\) from \(\mathrm{\epsilon }\). For this purpose, refer to document [R5.03.28] in which these operations are detailed. Although these quantities depend on the natural deformations, it can be noted that the derivation methodology that is implemented does not require the derivative of associated eigenvectors which, as we recall, are not defined when the eigenvalues are not distinct.