2. Models of elasto-plasticity-damage in general

2.1. Challenges and difficulties

In modeling the behavior of solids, we note that the most important non-linear phenomena to take into account, during a structure analysis under extreme load, concern the presence and evolution of anelastic deformation (plasticity) and the weakening of stiffness (damage). The relative importance of these two phenomena can vary drastically from one material to another. Thus, most often only plasticity is considered when modeling metals while neglecting damage, while for a concrete-type material, damage is considered more important and plasticity is not taken into account. However, this simplification is often not justifiable and it is necessary to model both phenomena, plasticity and damage, at the same time.

For some time now the thermodynamic framework for modeling both phenomena has been well established (e.g. see [bib1], [bib2]) and there is no obstacle to extending it to elasto-plastic-damageable couplings. As for the numerical resolution of these models within the framework of the finite element method (MEF), there are very robust algorithms in particular for plasticity (e.g. see [bib13]) and also for damage (e.g. see [bib14]) provided that the localization phenomenon does not become predominant. On the other hand, the numerical resolution of a coupled model is much less obvious.

Numerical difficulties concerning the algorithm for solving a coupled model arise mainly in two aspects:

• satisfaction of the conditions (Kuhn-Tucker) of the plasticity and damage thresholds,

• construction of the tangent matrix.

The first difficulty is due to the fact that plastic and damageable internal variables act very differently on the stress level, making the system of equations difficult to solve. The second point is a simple consequence of the increased complexity of the coupled model compared to single-phenomenon models.

Recently, an algorithm for solving an elasto-plastic-damagable model based on strict partitioning into a plastic part and a damaging part was proposed in [bib3] and [bib4]. It is about introducing an additional iterative process in order to satisfy plastic eligibility and damaging admissibility simultaneously. This iterative process is driven by what is called damage deformation, a variable that, unlike plastic deformation, is not a state variable, but rather represents a means of coupling between different models. Thus, in the model introduced in [bib3], the choice of the elasto-plastic model being free, it was nevertheless necessary to modify the damage formulation, which in this case only provided for the modeling of non-softening behaviors. In [bib4] we demonstrated the interest of integrating this type of law in the framework of mixed finite elements.

In this document we generalize the approach presented in [bib3] in order to be able to couple practically any elasto-plastic model with any damage model. The advantage of such a platform actually goes beyond the elasto-plastic-damageable coupling.

2.2. Thermodynamic framework

In this part we recall the usual approach to build models of damage and elasto-plasticity, based on thermodynamics. The combination of the two types of model can be done in a fairly natural way, provided that the coupling is done at the level of the calculation of the constraints. The entire approach presented here is carried out under the assumption that the internal variables of damage and plasticity are not directly coupled to each other. In addition, we are dealing with small deformations and without thermal phenomena.

2.2.1. Plasticity

Free energy \(\psi\) can be written as the sum of elastic energy and so-called stored energy, due to plasticity:

(2.52)\[ \ psi (\ varepsilon, {\ varepsilon} ^ {p},\ alpha)\ mathrm {=} {\ psi} ^ {e} (\ varepsilon\ mathrm {-} {-} {\ varepsilon} {\ varepsilon} ^ {p}) + {\ psi} ^ {p} (\ alpha)\]

where \(\varepsilon\) is total deformation, \({\varepsilon }^{p}\) is plastic deformation, and \(\alpha\) is a vector of internal variables. The essential assumption we make in () is that elastic energy depends only on elastic deformation, \({\psi }^{e}\mathrm{=}{\psi }^{e}({\varepsilon }^{e})\), where \({\varepsilon }^{e}\mathrm{=}\varepsilon \mathrm{-}{\varepsilon }^{p}\).

In the second step, the plastic dissipation rate is defined and it is ensured that it always remains positive,

(2.52)\[ {D} ^ {p}\ mathrm {=}\ sigma\ mathrm {:}\ dot {\ varepsilon}\ mathrm {-}\ rho\ dot {\ psi}\ mathrm {=}\ mathrm {=} (\ sigma\ mathrm {=}} (\ sigma\ mathrm {-} -}\ rho\ psilon}\ mathrm {-}\ mathrm {-}\ rho\ frac {\ mathrm {\ partial} {\ psi} {\ psi}\ mathrm {=}}\ mathrm {=}} (\ sigma\ mathrm {-}\ rho\ psilon}\ mathrm {-}\ rho\ psilon}\ mathrm {-}\ mathrm {\ partial} {\ varepsilon} ^ {e}})\ mathrm {:} {\ dot {\ varepsilon}} ^ {e} +\ sigma\ mathrm {:} {\ dot {\ varepsilon}}\ dot {\ varepsilon}} {\ dot {\ varepsilon}}} {\ dot {\ alpha}\ dot {\ alpha}\ dot {\ alpha}\ dot {\ alpha}\ dot {\ alpha}\ dot {\ alpha}\ dot {\ alpha}\ dot {\ alpha}\ dot {\ alpha}\ dot {\ alpha}\ dot {\ varepsilon}})\ mathrm {\ psilon}}}\ mathrm {-} Age} 0\]

Where

(2.52)\[ A=\ rho\ frac {\ partial {\ psi} ^ {p}} {\ partial\ alpha}\]

are the thermodynamic forces associated with internal variables and the symbol \(\ast\) refers to the appropriate product operator between \(A\) and \(\dot{\alpha }\).

We assume that the dissipation does not change, i.e. \(\dot{D}=0\), when the plastic variables do not change, i.e. \({\dot{\epsilon }}^{p}=\dot{\alpha }=0\), which allows us to define the constraint as:

(2.52)\[ \ sigma =\ rho\ frac {\ partial {\ psi} ^ {e}} {\ partial {\ epsilon} ^ {e}}\]

In this document, we are content with materials that have linear and isotropic behavior. In this case, the elastic contribution of free energy is given by:

(2.52)\[ \ rho\ psi ({\ epsilon} ^ {e}) =\ frac {1} {e}) =\ frac {1} {2} {\ epsilon} ^ {e}\ mathrm {:} {\ epsilon} ^ {e}\ mathrm {:} {\ epsilon} ^ {e}\]

Where \({D}^{e}\) is the standard isotropic elasticity tensor. From the two previous equations, we deduce the expression for the stress tensor:

(2.52)\[ \ sigma = {D} ^ {e}\ mathrm {:} {\ epsilon} ^ {e}\]

Another essential ingredient in an elasto-plastic model is the definition of the flow threshold function, which depends on the stress tensor and the thermodynamic forces associated with the internal variables. This function is negative when the deformations are purely elastic and reaches zero when plastic flow is imminent:

(2.52)\[ {\ Phi} ^ {p} = {\ Phi} ^ {p}\ left (\ sigma, A\ right)\ le 0\]

The complete characterization of the general plasticity model requires the definition of the laws of evolution of the internal variables, i.e. the variables associated with the phenomena of dissipation. In our case, the internal variables are the plastic strain tensor and the set of internal variables \(\alpha\). To establish the laws of evolution of these internal variables, the principle of maximum plastic dissipation is often used. However, this is not essential in our approach. We will use more general writing:

(2.52)\[ \ begin {array} {c} {\ dot {\ epsilon}}} ^ {epsilon}}} ^ {p} = {\ dot {\ gamma}} N= {\ dot {\ gamma}} ^ {p}\ frac {\ partial\ psi}} {\ partial\ psi}} {\ partial\ psi}} {\ partial\ psi}} {\ partial\ psi} {\ psi} {\ psi} {\ partial\ psi} {\ psi} {\ partial\ psi} {\ psi} {\ partial\ psi} {\ psi} {\ partial\ psi} {\ psi} {\ partial\ psi} {\ psi} {\ partial\ psi} {\ psi} {\ partial\ psi} {\ partial\ psi} {\ psi} {\ partial\ psi} {\ partial\ psi} {\ sigma}} ^ {p}\ frac {\ partial\ psi} {\ partial A}\ end {array}\]

\(\dot{\gamma }\) is the plastic multiplier, \(N\) is called the flow vector, and \(H\) is the work hardening module. In addition, any solution must meet the conditions of plastic admissibility, which are:

(2.52)\[ {\ Phi} ^ {p}\ le\ mathrm {0,}\ dot {\ gamma}\ ge\ mathrm {0,} {\ Phi} ^ {p}\ dot {\ gamma} =0\]

2.2.2. Damaging

For damage models it is more difficult to have a construction as general as for that of plasticity. Consequently, we will deal with a simple case that is fairly representative of the desired characteristics of the model. Thus, free energy is defined as:

(2.52)\[ \ psi\ left (\ epsilon, d\ right) = {\ psi} ^ {e}\ left (\ epsilon\ right)\ zeta\ left (d\ right)\]

where \(\zeta\) is the damage function and \(d\) is the damage variable, which is a scalar but can also be a tensor. Likewise, the \(\zeta\) function may in general depend on \(\epsilon\), in particular on the sign of its trace to distinguish behaviors that are typically very different in traction and compression.

Using (), damage dissipation is defined as,

(2.52)\[ {D} ^ {d} =\ sigma\ mathrm {:}\ dot {\ epsilon} -\ rho\ dot {\ psi} =\ left (\ sigma -\ rho\ frac {\ partial\ psi} {\ partial\ psi} {\ partial\ psi}} {\ partial\ psi} {\ right)\ mathrm {:}\ dot {\ epsilon} -\ rho\ frac {\ partial\ psi} {\ partial\ psi} {\ right)\ mathrm {:}\ dot {\ epsilon} -\ rho\ frac {\ partial\ psi} {partial d}\ dot {d}\ ge 0\]

The dissipation does not change, i.e. \({\dot{D}}^{d}=0\), when the damage variable does not change, i.e. \(\dot{d}=0\), which allows us to define the stress as:

(2.52)\[ \ sigma =\ rho\ frac {\ partial\ psi} {\ partial\ psi} {\ partial\ epsilon} {\ partial\ epsilon}} {\ partial\ epsilon}\ zeta\ left (d\ right)} {\ left (d\ right) = {D} ^ {right) = {D} ^ {e}\ left (d\ right)\ left (d\ right)\ mathrm {:} {\ epsilon}\ zeta\ left (d\ right) = {D} ^ {e}\ left (d\ right)\ left (d\ right)\ mathrm {:} {\ epsilon}\ zeta\ left (d\ right)\]

Where \({D}^{e}(d)\) represents the standard isotropic elasticity tensor that changes as a function of damage. Assuming isotropic damage, taking \(\zeta \mathrm{=}1\mathrm{-}d\) the constraint becomes:

(2.52)\[ \ sigma\ mathrm {=} (1\ mathrm {-} d) {D} ^ {e}\ mathrm {:} {\ varepsilon} ^ {e}\]

As for plasticity, a damage threshold function is introduced:

(2.52)\[ {\ Phi} ^ {d}\ mathrm {=} {\ Phi} ^ {d} (\ varepsilon, Y)\ mathrm {\ the} 0\]

Where \(Y\) is the dual variable to \(d\). This time, the threshold is a direct function of the total deformation, \(\varepsilon\), and the aim is to satisfy the conditions for the admissibility of damage:

(2.52)\[ {\ Phi} ^ {d}\ mathrm {\ le}\ mathrm {\ le}\ mathrm {\ le}\ mathrm {\ ge}\ mathrm {0,} {\ Phi} {\ Phi} ^ {d}\ dot {d}\ mathrm {=} 0\]

We must also formulate the \(d\) law of evolution:

(2.52)\[ \ dot {d} =\ dot {d}\ left (Y\ right)\]

2.2.3. Coupling

Finally, damage and plasticity formulations are combined to arrive at a coupled model. The free energy of the coupled model is therefore written as:

(2.52)\[ \ psi (\ varepsilon, {\ varepsilon} ^ {p},\ alpha, d)\ mathrm {=} {\ psi} ^ {e} (\ varepsilon\ mathrm {-} {-} {\ varepsilon} {\ varepsilon} {p})\ zeta (d) + {\ psi} ^ {p} (\ alpha)\]

Dissipation is written as:

(2.52)\[ D\ mathrm {=}\ sigma\ mathrm {::}\ dot {\ varepsilon}\ mathrm {-}\ rho\ dot {\ psi}\ mathrm {=} (\ sigma\ mathrm {-} (\ sigma\ mathrm {- -}}\ rho\ frac {\ frac {\ frac} {\ partial} {\ partial} {\ partial} {\ partial}} {\ mathrm {\ partial}} {\ mathrm {\ partial}} {\ partial}} {\ mathrm {\ partial}} {\ partial} {\ partial}} {\ mathrm {\ partial} {\ partial}} {\ mathrm {\ partial} {\ partial}} {\ mathrm {\ partial} {\ partial} {varepsilon} ^ {e}})\ mathrm {:} {:} {\ dot {\ varepsilon}} ^ {e} +\ sigma\ mathrm {:} {\ dot {\ varepsilon}}}} ^ {p}} ^ {p}\ mathrm {-}\ mathrm {-}\ mathrm {-}\ frac {-} A\ mathrm {-} A\ mathrm {-}\ mathrm {-} mathrm {\ partial} {\ psi} ^ {d}} {\ mathrm {\ partial} d}\ dot {d}\ mathrm {\ ge} 0\]

The stress tensor is equal to:

(2.52)\[ \ sigma\ mathrm {=}\ rho\ frac {\ mathrm {\ partial} {\ mathrm {\ partial}} {\ mathrm {\ partial} {\ varepsilon} ^ {e}}\ mathrm {=}\ mathrm {=}\ mathrm {\ partial} {\ partial}} {\ mathrm {\ partial}} {\ mathrm {\ partial}} {\ mathrm {\ partial}} {\ mathrm {\ partial}} {\ mathrm {\ partial}} {\ varepsilon} ^ {e}}\ zeta (d)\ mathrm {=} {D} ^ {e} (d)\ mathrm {:} (\ varepsilon\ mathrm {-} {\ varepsilon} ^ {p})\]

If the elasticity is isotropic linear and the damage is isotropic:

(2.52)\[ \ sigma\ mathrm {=} (1\ mathrm {-} d) {D} ^ {e}\ mathrm {:} (\ varepsilon\ mathrm {-} {-} {\ varepsilon} ^ {p})\]

We then have two threshold functions: one for plastic flow and the other for damage:

(2.52)\[\begin{split} \ begin {array} {c} {\ Phi} ^ {p} ^ {p}\ mathrm {=} {\ Phi} ^ {p} (\ sigma, A)\ mathrm {\ le} 0\\ {\ Phi} ^ {d} ^ {d} ^ {d}\ mathrm {d}\ mathrm {\ the} 0\ end {array}\end{split}\]

The laws of evolution of the damage variables are given by:

(2.52)\[\begin{split} \ begin {array} {c} {\ dot {\ varepsilon}}} ^ {p}} ^ {p}\ mathrm {=}} ^ {p} N\ mathrm {=} {\ dot {\ gamma}} {\ dot {\ gamma}}}} ^ {\ gamma}}} ^ {p}\ frac {\ mathrm {\ partial}}\\ dot\ alpha}\ mathrm {=} {\ dot {\ gamma}}} ^ {gamma}} ^ {p} H\ mathrm {=}\ mathrm {-} {\ dot {\ gamma}} ^ {p}\ frac {\ dot {\ gamma}}\ frac {\ mathrm {\ gamma}}\\ gamma}} {\ dot {d}} ^ {p}\ frac {\ mathrm {\ gamma}}\ frac {\ mathrm {\ gamma}}\ frac {\ mathrm {\ gamma}}\ frac {\ mathrm {\ partial}}\ frac {\ mathrm {\ partial}\\ partial}\ psi} {\ mathrm {\ partial}\ psi} {\ mathrm {\ partial}\ psi} {Y)\ end {array}\end{split}\]

Eligibility requirements are:

(2.52)\[ \ begin {array} {c} {\ Phi} ^ {p}\ mathrm {\ le}\ mathrm {\ le}\ mathrm {\ gamma}\ mathrm {\ ge}\ mathrm {0,} {\ Phi} ^ {p}\ mathrm {\ le}\ {p}\ dot {\ le}\ mathrm {\ le}\ mathrm {0,}\ dot {d}\ mathrm {\ ge}\ mathrm {\ ge}\ mathrm {0,} {\ phi} ^ {d}\ dot {d}\ mathrm {=} 0\ end {array}\]

2.3. Digital integration

2.3.1. Reminder of the classical resolution of nonlinear problems using finite elements

In the proposed method, we add an additional iterative process compared to a classical calculation using MEF, which we recall in this chapter. In an overall resolution of a structural problem, the state of equilibrium is sought, defined by:

(2.52)\[ {f} ^ {\ text {int}}} (u (t)) = {f} ^ {\ mathit {ext}} (t)\]

Where \({f}^{\text{int}}\) is the vector of internal forces and \({f}^{\mathit{ext}}\) is the vector of external forces (the load). The first is a function of the displacement vector, the main variable of the system, and the second a function of the parameter*t* (pseudo-time), defining the concept of trajectory both for loading and for the evolution of the displacement field. In the context of MEF, internal strength is determined by:

(2.52)\[ {f} ^ {\ text {int}} (u)\ mathrm {=} {\ mathrm {\ int}} _ {\ Omega} {B} ^ {T}\ sigma (\ varepsilon) d\ Omega\]

where \(\Omega\) defines the domain in question, \(\varepsilon\) the strain tensor, the strain tensor, \(\sigma\) the stress tensor, and \(B\) is the classical matrix corresponding to the discretized symmetric gradient operator that allows the transition from deformations to displacements.

As we considered the hypothesis of small deformations, tensor \(\varepsilon\) is defined as:

(2.52)\[ {\ epsilon} _ {\ mathit {ij}}\ mathrm {=}\ frac {1} {2} (\ frac {\ mathrm {\ partial} {u} _ {i}}} {\ mathrm {\ partial}}} {\ partial}} {\ partial}}} {\ mathrm {\ partial} {u}}} {i}}} {\ mathrm {\ partial}}} {\ mathrm {\ partial} {x} _ {i}})\]

Using the definition of \(B\), we have:

(2.52)\[ {\ varepsilon} _ {\ mathit {ij}}\ mathrm {=} {B} _ {\ mathit {ijk}} {u} _ {k}\]

After discretizing the field of travel with MEF. When the behavior of the material is non-linear, Newton’s method is most often used to solve (). To do this, we need to linearize (), which gives us:

(2.52)\[ {f} ^ {\ text {int}, (k+1)}\ mathrm {=} {f} ^ {\ text {int}, (k)} +\ frac {\ mathrm {\ partial} {f} {f}} ^ {\ text {int} {f} ^ {f} ^ {f} ^ {f} ^ {f} ^ {f} ^ {k)}}\ {u} ^ {f} ^ {(k)} ^ {k)}\]

Where the index \(k\) represents the current iteration number. At convergence we must have:

(2.52)\[ {f} ^ {\ text {int}, (k)}\ stackrel {k\ to\ mathrm {\ infty}} {\ to} {f} ^ {\ mathit {ext}} ^ {\ mathit {ext}} (t)\]

In order to be able to calculate the displacement correction from one iteration to another, \(\Delta {u}^{(k)}\), we must determine:

(2.52)\[ K\ mathrm {=}\ frac {\ mathrm {\ partial} {f} {f} ^ {\ text {int}, (k)}} {\ mathrm {\ partial} {u} ^ {(k)}}}\ mathrm {=}}\ mathrm {=} {=}\ underset {\ omega} {\ omega} {\ int}} {B} ^ {T} {D} {D}}}}\ mathrm {=}}\ mathrm {=}}\ mathrm {=}}\ mathrm {=}}\ underset {\ omega} {\ int}} {B} ^ {T} {D} {D}}}}\ mathrm {=}} {\ mathit} {epd}} Bd\ Omega\]

Where \({D}^{\mathit{epd}}\) is the coherent tangent module for the law of behavior used:

(2.52)\[ {D} ^ {\ mathit {epd}}\ mathrm {=}\ frac {\ mathrm {\ partial}\ sigma (\ varepsilon)} {\ mathrm {\ partial}\ varepsilon}\]

The method detailed below is positioned at the level of a material point (or even at the Gauss point) and intervenes at the level of the calculation of \(\sigma\) and \({D}^{\mathit{epd}}\) from the deformation \(\epsilon\), given for each call of the global process of Newton’s iterations.

2.3.2. Approach to the coupled problem

This document does not detail the numerical integration over time for each of the models. We simply assume that we know how to do it and this with an implicit Euler diagram (back). Let’s say, that in an incremental resolution by the finite element method, for each Newton iteration of the global iterative process, we want to calculate, for each Gauss point, the constraints \({\sigma }_{n+1}\) and the internal variables \({\epsilon }_{n+1}^{p}\), \({\alpha }_{n+1}\) and \({d}_{n+1}\) from the values of the total deformation, and from the values of the total deformation, \({\epsilon }_{n+1}\), while assuming that we converged at the previous moment and that we know \({\sigma }_{n}\), \({\epsilon }_{n}^{p}\), \({\alpha }_{n}\), and \({d}_{n}\). The indices \(n\) or \(n+1\) indicate that the value of the variable corresponds to time \({t}_{n}\) or \({t}_{n+1}\), respectively.

Thus it is assumed that the two models can be integrated in decoupled mode and that we are therefore able to write:

(2.52)\[ {\ sigma} _ {n+1} ^ {p} = {\ sigma} ^ {p}\ left ({\ stackrel {} {\ epsilon}}} _ {n+1}\ right)\]

and

(2.52)\[ {\ sigma} _ {n+1} ^ {d} = {\ sigma} ^ {d}\ left ({\ stackrel {} {\ epsilon}}} _ {n+1}\ right)\]

Exponents \(d\) or \(p\) in \({\sigma }^{d}\) and \({\sigma }^{p}\) mean that the stress is calculated with the damage model and with the plasticity model, respectively. \(\stackrel{̃}{\epsilon }\) is a deformation-type variable.

In the algorithm presented below, an additional iterative process is introduced so that the constraints \({\sigma }^{d}\) and \({\sigma }^{p}\) converge on a single value. It is similar to the one proposed in [bib3] except that here the construction of the damage model is more general, which requires considering a variable other than the damage deformation (used in [bib3]) on which one is iterating.

2.3.3. Construction of the nonlinear system

Compared to a classical resolution of the law of behavior, two stress values must be converged here by introducing another unknown such as deformation. This partitioned system is more easily constructed by introducing the following decomposition of the total deformation,

(2.52)\[ \ epsilon = {\ epsilon} ^ {e} + {\ epsilon} ^ {p} + {\ epsilon} ^ {d}\]

Where we introduced damage deformation, which allows us to write the elasticity relationship as:

(2.52)\[ \ sigma = {D} ^ {e}\ mathrm {:}\ left (\ epsilon - {\ epsilon} ^ {p} - {\ epsilon} ^ {d}\ right)\]

or even:

(2.52)\[ {\ epsilon} ^ {d}\ left (\ sigma,\ epsilon, {\ epsilon, {\ epsilon} ^ {p}\ right) =\ epsilon - {\ epsilon} ^ {p} - {p} - {({D}} - {D} ^ {e})} ^ {e})} ^ {-1}\ mathrm {:}\ sigma\]

We can summarize the equations to be solved in three groups:

Plasticity modulus

(2.52)\[ {\ sigma} _ {n+1} ^ {p} = {\ sigma} ^ {p} {\ left ({\ epsilon} _ {n+1} - {\ epsilon} _ {n+1} _ {n+1}} ^ {d}\ right)} _ {e} {\ epsilon} _ {n+1} ^ {p} = {\ epsilon} ^ {p}\ left ({\ epsilon} _ {n+1} - {\ epsilon} _ {n+1} _ {n+1} = {n+1}} ^ {d}\ right)\]

The index \(e\) means that the elastoplastic law is calculated with the elasticity operator \({D}^{e}\)

Damage module

(2.52)\[ {\ sigma} _ {n+1} ^ {d} = {\ sigma} ^ {d} {\ left ({\ epsilon} _ {n+1} - {\ epsilon} _ {n+1} _ {n+1}} ^ {p}\ right)} _ {e}\]

Coupling module

(2.52)\[ {\ epsilon} _ {n+1} ^ {d} = {\ epsilon} = {\ epsilon} _ {n+1} ^ {p} - {(D} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ {e})} ^ { {\ sigma} _ {n+1} ^ {p} = {\ sigma} _ {n+1} ^ {d}\]

In [bib3] we chose \({\epsilon }_{n+1}^{d}\) as the additional unknown, which is not possible in our case, since \({\sigma }^{d}\) is not a function of \({\epsilon }^{d}\) but of \({\stackrel{̄}{\epsilon }}^{\mathit{ed}}=\epsilon -{\epsilon }^{p}\), the variable that is finally used to drive the iterative processes described below. Moreover, it is also worth \({\stackrel{̄}{\epsilon }}^{\mathit{ed}}={\epsilon }^{e}+{\epsilon }^{d}=\epsilon -{\epsilon }^{p}\).

2.3.4. Resolution algorithm

To solve the non-linear system, defined in the equations (to) by assigning them to the various modules (plasticity, damage and coupling), the classical Newton method is applied. In order to linearize the system, the following derivatives must be defined:

(2.52)\[ {D} ^ {\ mathit {ep}}}\ left ({\ stackrel {tante} {\ epsilon}} ^ {\ mathit {ed}}\ right) =\ frac {\ partial {\ sigma} _ {n+1}} _ {n+1}} ^ {n+1} ^ {n+1} ^ {n+1} ^ {n+1} ^ {\ mathit {ed}}}} {n+1} ^ {n + 1} ^ {p}}} {\ partial {\ stackrel {e}}} {D} ^ {d}\ left ({\ stackrel {tante} {\ epsilon}}} ^ {\ mathit {ed}}\ right) =\ frac {\ partial {\ sigma} _ {n+1}} ^ {d}} ^ {d}} {\ d}} {\ d}}} {\ partial {\ stackrel {ed}}} {\ mathit {ed}}}} \ frac {\ partial {\ epsilon} _ {n+1} ^ {d}} {\ partial {\ stackrel {]} {\ epsilon}} ^ {\ mathit {ed}}}}} =1- {\ left ({D} ({D} ^ {e}\ right)} ^ {e}\ right)} ^ {e}\ right)} ^ {e}\ right)} ^ {-1} {D} ^ {d}\ left ({\ stackrel {!}} _ {n+1} ^ {\ mathit {ed}}\ right)\]

The non-linear system (à) is first brought back to the last of the equations, i.e. to \({\sigma }_{n+1}^{p}={\sigma }_{n+1}^{d}\), since we can express all the unknowns with respect to the value of \({\stackrel{̄}{\epsilon }}_{n+1}^{\mathit{ed}}\). We can therefore define the residue as:

(2.52)\[ {R} _ {\ sigma} ^ {(k)}} = {\ sigma} ^ {k)} = {\ sigma} ^ {p}\ left ({\ epsilon} _ {n+1} ^ {d, (k)}\ right) - {\ sigma} ^ {d}}\ right) - {\ sigma} ^ {d}\ right) - {\ sigma} ^ {d}}\ right) - {\ sigma} ^ {d}}\ right) - {\ sigma} ^ {d}}\ right) - {\ sigma} ^ {d}}\ right) - {\ sigma} ^ {d}}\ right) - {\ sigma} ^ {d, (k)}\ right) - {\ sigma} ^ {d, (k)}\ right) - {\ sigma} ^ {d}} Mathit {ed}, (k)}\ right)\]

where the exponent \(k\) indicates the current iteration and where:

(2.52)\[ {\ epsilon} _ {n+1} ^ {d, (k)} = {\ epsilon} _ {n+1} - {\ epsilon} _ {n+1} ^ {p, (k)} - {\ left ({D}, (k)}} - {\ left ({D}, (k)} {d}} - {d}} - {\ left ({D} ^ {d}\ right)} ^ {e}\ right)} ^ {e}\ right)} ^ {e}\ right)} ^ {e}\ right)} ^ {e}\ right)} ^ {e}\ right)} ^ {d}\ left ({\ stackrel} {d}\ left ({\ stackrel} {d}\ left ({\ stackrel}\ left {]] {\ epsilon}} _ {n+1} ^ {\ mathit {ed}, (k)}\ right)\]

Applying Newton’s method, we calculate the new approximation for \({\stackrel{̄}{\epsilon }}_{n+1}^{\mathit{ed},(k+1)}\) as:

(2.52)\[ {\ stackrel {tante} {\ epsilon}}} _ {n+1}}} _ {n+1}} _ {n+1}} _ {n+1} ^ {n+1} ^ {\ mathit {ed}, (k)} ^ {\ mathit {ed}, (k)} ^ {\ mathit {ed}, (k)} ^ {d}} +\ delta {\ stackrel {! it {ed}, (k)}\]

Where

(2.52)\[ \ Delta {\ stackrel {Belgium} {\ epsilon}}} _ {n+1}}} _ {\ epsilon}} =- {\ left (\ frac {\ partial {R} _ {\ sigma}} {\ sigma} ^ {\ epsilon}} {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {\ sigma} ^ {(k)}}}} {\ partial {(k)}}} {\ partial {\ stackrel {!\ right)} ^ {-1} {R} {R} _ {\ sigma} ^ {(k)}\]

The corresponding tangent module can be expressed entirely by the decoupled modules ():

(2.52)\[ \ left (\ frac {\ partial {R} _ {\ sigma}} {\ sigma} ^ {(k)}} {\ partial {\ stackrel {]} {\ epsilon}} _ {n+1} ^ {\ mathit {ed}}}}\ right) =- {\ sigma}}}\ right) =- {D} ^ {D} ^ {D} ^ {\ mathit {ep}, (k)}\ frac {\ partial {\ epsilon} ^ {\ mathit {epsilon}} ^ {\ mathit {ed}}}}\ right) =- {D} ^ {D} ^ {D} ^ {\ mathit {ep}, (k)}\ frac {\ partial {\ epsilon}} ^ {\ mathit {\ epsilon}} ^ {d, (k)}} {\ partial {\ stackrel {e} {\ stackrel {ounty}} {\ epsilon}}} _ {n+1} ^ {d, (k)}} - {D} ^ {d, (k)}} =-\ left (k)} =-\ left (k)} =-\ left (1)} =-\ left (1- {\ left ({D} ^ {e}\ right)} ^ {e}\ right)} ^ {-1} {d} ^ {d, (k)}\ right)\ right)\]

We recall here that the total deformation \({\epsilon }_{n+1}\) , resulting from the global process of Newton iterations, is constant for the determining local iterative process \({\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed}}\). The algorithm for calculating the variables in iteration \(k+1\) is performed as follows:

0. Initialization: \({\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k\mathrm{=}0)}\mathrm{=}{\varepsilon }_{n+1}\mathrm{-}{\varepsilon }_{n}^{p}\)

1. Loop on the \(k\mathrm{=}\mathrm{0,}\mathrm{...},{N}_{\mathit{max}}\)

2. Calculation of the increment:

\({\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k+1)}\mathrm{=}{\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k)}+\Delta {\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k)}\) where \(\Delta {\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k)}\mathrm{=}{({D}^{d,(k)}+{D}^{\mathit{ep},(k)}(1\mathrm{-}{({D}^{e})}^{\mathrm{-}1}{D}^{d,(k)}))}^{\mathrm{-}1}{R}_{\sigma }^{(k)}\)

3. \({\sigma }_{n+1}^{d,(k+1)}\mathrm{=}{\sigma }^{d}{({\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k+1)})}_{e}\)

4. \({\varepsilon }_{n+1}^{d,(k+1)}\mathrm{=}{\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k+1)}\mathrm{-}{\varepsilon }_{n+1}^{\mathit{el},(k+1)}\) where \({\varepsilon }_{n+1}^{\mathit{el},(k+1)}\mathrm{=}{({D}^{e})}^{\mathrm{-}1}\mathrm{:}{\sigma }^{d}({\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k+1)})\)

5. \({\sigma }_{n+1}^{p,(k+1)}\mathrm{=}{\sigma }^{p}{({\varepsilon }_{n+1}\mathrm{-}{\varepsilon }_{n+1}^{d,(k+1)})}_{e}\)

6. \({R}_{\sigma }^{(k+1)}\mathrm{=}{\sigma }_{n+1}^{p,(k+1)}\mathrm{-}{\sigma }_{n+1}^{d,(k+1)}\)

7. If \(\mathrm{\mid }\Delta {\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed},(k)}{R}_{\sigma }^{(k+1)}\mathrm{\mid }>\mathit{tol}\) resumes the loop at 1), otherwise FIN

**Note 1: the only places where the physical specificities of the coupled models are taken into account are points 3) and 5), where the plasticity and damage modules are used, respectively.*

**Note 2: The value of* \({\varepsilon }_{n+1}\) is considered constant in the process of balancing the constraints and when changing of \({\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed}}\) we actually vary the value of \({\varepsilon }_{n+1}^{p}\) . This plastic deformation is in principle determined by the elasto-plasticity module as a function of \({\varepsilon }_{n+1}\) (see ()). In reality, \({\stackrel{ˉ}{\varepsilon }}_{n+1}^{\mathit{ed}}\) is considered here as a variable and its link with \({\varepsilon }_{n+1}^{p}\) is only restored at the convergence of the iterative process. So:

(2.52)\[ {\ stackrel {} {\ varepsilon}}} _ {n+1}}} _ {n+1}} _ {n+1}\ mathrm {-} {n+1}\ mathrm {-} {-} {\ varepsilon} ^ {p} ({\ varepsilon} _ {n+1}\ mathrm {-}} {\ varepsilon} _ {n+1} ^ {d, (k)})\]

2.3.5. Construction of the coherent tangent module

Finally, in order to obtain optimal performance in the global iterative process, the coherent tangent module must be calculated, once the process described above has been converged, i.e. \({\sigma }_{n+1}={\sigma }_{n+1}^{p}={\sigma }_{n+1}^{d}\):

(2.52)\[ {D} _ {n+1} ^ {\ mathit {epd}}\ mathrm {=}\ frac {\ mathrm {\ partial} {\ sigma} _ {n+1}} {\ mathrm {\ partial}} {\ mathrm {\ partial}} {\ partial} {\ partial} {\ varepsilon} _ {n+1}}\]

First we observe that if we derive () with respect to \({\sigma }_{n+1}^{p}\) we can write:

(2.52)\[ I= {D} _ {n+1} ^ {\ mathit {ep}}}\ left ({D} _ {n+1} ^ {\ mathit {epd}} -\ frac {\ partial {\ epsilon} _ {n+1}} _ {n+1}} _ {n+1} {\ partial {\ epsilon}} _ {n+1} {\ epsilon} _ {n+1} {\ epsilon} _ {n+1} {\ epsilon} _ {n+1} {\ epsilon} _ {n+1} {\ epsilon} _ {n+1} {\ epsilon} _ {n+1} {\ epsilon} _ {n+1} {\ epsil\]

and that in the same way, by deriving () with respect to \({\sigma }_{n+1}^{d}\) we obtain:

(2.52)\[ I= {D} _ {n+1} ^ {d}\ left ({D} _ {n+1} ^ {\ mathit {epd}}} -\ frac {\ partial {\ epsilon} _ {n+1}} ^ {d}} ^ {d}} ^ {p}}} {\ p}} {\ left}} {\ left}} {\ p}} {\ left}} {\ partial {\ sigma} _ {n+1} ^ {d}}\ right)\]

with, at the end of the iterative process \({\sigma }_{n+1}={\sigma }_{n+1}^{p}={\sigma }_{n+1}^{d}\), we finally get:

(2.52)\[ \ frac {\ partial {\ epsilon} _ {n+1}} ^ {d}} {\ partial {\ sigma} _ {n+1}} = {({D} _ {n+1} ^ {\ mathit {epd}} ^ {epd}})})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}})} ^ {epd}}\]

By then using \({\epsilon }_{n+1}={\epsilon }_{n+1}^{e}+{\epsilon }_{n+1}^{p}+{\epsilon }_{n+1}^{d}\) and therefore that:

(2.52)\[ \ frac {\ partial {\ epsilon} _ {n+1}}} {\ partial {\ sigma} _ {n+1}} =\ frac {\ partial {\ epsilon} _ {n+1} ^ {e}}} {\ partial {\ sigma}}} {\ partial {\ sigma}} _ {n+1}} _ {n+1}} ^ {e}}} {\ partial {e}}} {\ partial {\ sigma} _ {n+1}} ^ {e}}} {\ partial {\ sigma}} _ {n+1}} ^ {e}}} {\ partial {\ sigma}} _ {n+1}} ^ {e}}} {\ partial {\ sigma}} {\ sigma} _ {n+1}} ^ {e}}} {\ partial {\ sigma}} {sigma} _ {n+1}} +\ frac {\ partial {\ epsilon} _ {n+1} ^ {d}} {\ partial {\ sigma} _ {n+1}}\]

we obtain the coherent tangent module directly as:

(2.52)\[ {D} _ {n+1} ^ {\ mathit {epd}}} = {\ left ({\ left ({D} _ {n+1} ^ {\ mathit {epd}}\ right)} ^ {-1}}}} ^ {-1}} + {\ left ({D} -1}} + {\ left ({D} -1}} + {\ left ({D} -1}} + {\ left ({D} -1}} + {\ left ({D} _ {n+1}} ^ {n+1} ^ {\ mathit {epd}}}\ right)} ^ {-1} + {\ left ({D} -1}} + {\ left ({D} -1}}} + {\ left ({D} -1}}} + {\ left ({D} -1}}} +1} ^ {\ mathit {epd}}\ right)}\ right)} ^ {-1}\ right)} ^ {-1}\]

2.4. Coupler limitations

There are a few conditions to be met for the use of this coupling method concerning both physical and numerical aspects:

• Distinct internal variables:

As specified in the introduction to the thermodynamic framework, the two coupled models should not share the same internal variables. The method is designed to be able to combine models representing decoupled physical phenomena, which is considered to be the case in most target applications.

• Existence of coherent tangent modules:

For the approach to work, it is essential that coherent tangent modules are available for both elasto-plastic and damageable modules. Under « coherent tangent module » we mean that it is the derivative coherent with the temporal discretization used of the stress tensor at the given moment in relation to the deformation tensor at the same instant.

• Non-singularity of tangent modules:

It is necessary several times to invert the tangent matrices of individual modules. If the behavior of a module has a flat non-linear slope (perfect elasto-plasticity for example) it is necessary to give it a non-zero value. This non-zero value must be chosen such that it is large enough for the algorithm to work and small enough for it to be consistent with the physical model being studied.