2. Variational formulation of the damage problem#
2.1. Case of a generic law#
Two equivalent approaches can be used to describe the process of damaging a fragile isotropic material. On the one hand, it is possible to derive the law of damage within the framework of generalized standard description. In this case it is necessary to define the free energy of the system, as well as the dissipation potential. The flow rule then establishes the evolution of the internal variables.
As for the description of damage we only need a scalar variable, the previous description is simplified and can be reduced to a variational problem under the stress of increasing damage [bib2].
To define a law of behavior with a damage gradient [R5.04.01] it is therefore sufficient to express the total free energy density (elastic+dissipation) as a function of deformation tensor \(\varepsilon\) and damage variable \(0\le a\le 1\). The spatial distribution of the damage is then given by a field \(a(x)\). The free energy density generally takes the following form:
\(\Phi (\varepsilon ,a)=A(a)w(\varepsilon )+\omega (a)+c/2{(\nabla a)}^{2}\) eq 2.1-1
Here \(c\) is the nonlocality parameter (C_ GRAD_VARI) \(w(\varepsilon )\) the elastic deformation energy, \(\omega (a)\) the dissipation energy, and \(A(a)\) the stiffness function. \(a=0\) corresponds to the healthy material and \(a=1\) corresponds to the completely damaged material: \(A(1)=\mathrm{0,}A(0)=1\). The evolution problem is now a simple Helmholtz free energy minimization problem \(F\equiv \int \Phi (\varepsilon ,a)d\Omega\) under \(\dot{a}\ge 0\) constraint [1] _ .
\({\mathrm{min}}_{(\varepsilon ,a)}F(\varepsilon ,a),\text{où}F(\varepsilon ,a)=\int \left[A(a)\varepsilon :E:\varepsilon +\omega (a)+c/2{(\nabla a)}^{2}\right]d\Omega\)
where we replaced \(w(\varepsilon )=\varepsilon :E:\varepsilon /2\) using the Hooke tensor definition. Two equations derive from the variational minimization problem: \(\delta F(\varepsilon ,a)/\delta \varepsilon =0\) [2] _ and \(\delta F(\varepsilon ,a)/\delta a\ge 0\). The inequality in the second equation is related to the presence of imposed stress. These two equations must be satisfied everywhere in integration domain \(\Omega\). They are completed by a Kuhn-Tucker \(\dot{a}\cdot \delta F(\varepsilon ,a)/\delta a=0\) coherence equation. On edges \(\partial \Omega\) we get an additional normality condition \(\nabla a\cdot n=0\), where \(n\) is a normal-vector. Finally, the damage variable and its gradient must be continuous within the integration domain to achieve the minimum of the functional in question (see [bib2,4] for more details).
2.2. Behavioral relationships#
The link between the variational formulation and the usual laws of evolution is direct. The condition of the material is characterized by deformation \(\varepsilon\) and damage \(a\), ranging between 0 and 1. The stress-strain relationship is defined, which remains elastic, and the stiffness is affected by the damage:
\(\sigma =\delta F(\varepsilon ,a)/\delta \varepsilon =A(a)E:\varepsilon\) eq 2.2-1
with \(E\) the Hooke tensor. The evolution of damage, which is always increasing, is governed by the following threshold function:
\(f(\varepsilon ,a)=-\delta \Phi (\varepsilon ,a)/\delta a=-\frac{1}{2}A\text{'}(a)\varepsilon :E:\varepsilon -\omega \text{'}(a)+c\Delta a\) eq 2.2-2
The coherence condition then takes on its usual form:
\(f(\varepsilon ,a)\le 0\dot{a}\ge 0\dot{a}f(\varepsilon ,a)=0\) eq 2.2-3
Two particularities of this formulation are noted. First, the threshold function is non-local because of the presence of the damaging Laplacian. Then, the absence of flow conditions is justified by the double role of damage \(a\), on the one hand it is presented as an internal evolution variable, on the other hand it fulfills the mission of the Lagrange parameter \(\lambda \equiv a\).
We also see the advantage of presenting damage laws in their variational form. It suffices to describe the total free energy density (eq.2.1-1), which includes dissipation, to completely define the law of evolution.
2.3. Identifying the parameters for law ENDO_SCALAIRE#
In law ENDO_SCALAIRE the stiffness and dissipation functions are chosen as follows:
\(\omega (a)=\mathrm{ka},A(a)={(\frac{1-a}{1+\gamma a})}^{2}\)
The parameters of this law of behavior are then five in number. On the one hand, the Young’s modulus \(E\) and the Poisson’s ratio \(\nu\) which determine the Hooke tensor by:
\({E}^{-1}\cdot \sigma =\frac{1+\nu }{E}\sigma -\frac{\nu }{E}(\text{tr}\sigma )\text{Id}\) eq 2.2-1
On the other hand, \(k,\gamma ,c\) that define the softening behavior, as well as the characteristic width of the damage band. The latter can be adjusted to the macroscopic parameters using the one-dimensional model, which admits a semi-analytical solution (bib 6,7). By noting the peak stress by \({\sigma }_{y}\), the energy of the Griffith break by \({G}_{f}\), and the size of the area damaged at break by \(D\) we get:
\(k=\frac{3{G}_{f}}{4D},c=\frac{3}{8}D{G}_{f},\gamma =\frac{3E{G}_{f}}{4{\sigma }_{y}^{2}D}-1\)
Numerical tests have shown that to avoid the presence of snap-backs in the 1D force-displacement response, it would be necessary to have \(\gamma \ge 2.8\). For this reason, the choice was made to simplify the entry of model data, we fill in not the complete set of macroscopic parameters \({\sigma }_{y},{G}_{f},D\), but directly the parameters of the model \(\gamma ,c\) and the peak constraint \({\sigma }_{y}\), given under the keywords factors ENDO_SCALAIRE (GAMMA, SY) and NON_LOCAL (C_ GRAD_VARI) of the operator. DEFI_MATERIAU. As for \(E\) and \(\nu\), they are classically given under the keyword factor ELAS or ELAS_FO. The following reasoning is only valid in the strict sense for 1D modeling, but may be useful for unsophisticated users. If the parameters \(E,\nu ,{G}_{f},{\sigma }_{y}\) are defined in advance, the user can vary the parameter \(D\) in order to satisfy the condition of absence of local snap-backs \(\gamma \ge 2.8\). He must then ensure that the size of the system in question is greater than the width of the damage band \(D\).
Example of concrete under traction |
\(\begin{array}{}E=30\text{GPa,}\nu =0.2\\ {G}_{f}=100\text{N/m}\\ {\sigma }_{y}=3\text{MPa}\end{array}\) |
|
The width of the damage belt must be chosen while respecting \(\gamma \ge 2.8\) <=> \(D\le 66\text{mm}\)
2.4. Integration of the law of behavior locally#
Here we present the method of integrating law ENDO_SCALAIRE in its local version (\(c=0\)), so that the user can make a generalization for the non-local case, which is generic and relies entirely on the algorithm presented in the doc. [R5.04.01].
The temporal discretization of equations [éq 2.2-1] to [éq 2.1-3] over a time step
is carried out by an implicit Euler diagram. Integrating the law of behavior in time consists in determining the state of stress and damage of the solution of the following nonlinear system:
\(\sigma =A(a)E:\varepsilon\) eq 2.4-1
\({f}_{\text{loc}}(\epsilon ,a)\le 0a-{a}^{-}\ge 0(a-{a}^{-})\cdot {f}_{\mathrm{loc}}(\epsilon ,a)=0\) eq 2.4-2
where variables without indices correspond to the final time step \(t\), such as deformation \(\varepsilon\); the state of the material at the start of time step \(({\varepsilon }^{\text{-}},{a}^{\text{-}})\) is indicated by the index « - ». The local threshold function is given by (eq. 2.2-2): \({f}_{\text{loc}}(\varepsilon ,a)=(1+\gamma )\frac{(1-a)}{{(1+\gamma a)}^{3}}\varepsilon :E:\varepsilon -k\)
A resolution method was proposed by [bib3]. It starts by examining the solution without evolution of the damage (also called elastic test) and then, if necessary, makes a correction to check the consistency condition. In the present case, the existence and the uniqueness of the solution guarantee the correct functioning of the method. Let us consider the elastic test:
\(a={a}^{-}\) solution if \({f}^{\text{el}}(\varepsilon )\equiv {f}_{\text{loc}}(\varepsilon ,{a}^{-})\le 0\) eq 2.4-3
Otherwise, the damage is obtained by solving \({f}_{\text{loc}}(\varepsilon ,a)=0\) (3rd order polynomial).
\((1\mathrm{-}a)(1+\gamma )\varepsilon \mathrm{:}E\mathrm{:}\varepsilon \mathrm{/}k\mathrm{=}{(1+\gamma a)}^{3}\) eq 2.4-4
It is the largest root that is chosen among the three existing.
It is still necessary to ensure that the damage does not exceed the value 1. In fact, when \(a=1\), the stiffness of the material point in question is cancelled out \(A(1)=0\). Insofar as no technique for removing « broken » finite elements is implemented (a technique that is possibly delicate when the finite elements have several Gauss points), zero pivots may appear in the stiffness matrix. This is why a numerical threshold of elastic residual stiffness is introduced for the tangent matrix, which can be entered under the keyword factor COEF_RIGI_MINI of the DEFI_MATERIAU operator. This dimensionless value is a multiplier of the elastic modulus of an isotropic linear model. To maintain reasonable conditioning of the stiffness matrix, we choose the value by default \(\mathrm{min}A(a)={10}^{-5}\).
An indicator \(\chi\), stored in the second internal variable, then specifies the behavior during the current time step:
elastic behavior (deformation energy below the threshold)
evolution of damage
saturated damage (\(a=1\)).
2.5. Integration of the law of behavior in non-local#
Here we only present the method of integrating law ENDO_SCALAIRE in its local version (\(c=0\)), because the generalization for the non-local case is generic and is based entirely on the algorithm presented in the doc. [R5.04.01]. Note that for the non-local version the threshold function is offset, so we obtain a polynomial of order 4 to be solved. As for the constraint, it is given by [éq 2.4-1] in all cases.
2.6. Description of internal variables#
There are three internal variables:
\(\mathrm{VI}(1)\) damage \(a\)
\(\mathrm{VI}(2)\) indicator \(\chi\)
\(\mathrm{VI}(3)\) residual stiffness \(1-A(a)\)