# Continuum damage modelling in geomechanics

Gilles Pijaudier-Cabot, Ludovic Jason, Antonio Huerta, Jean-François Dubé

► To cite this version:

Gilles Pijaudier-Cabot, Ludovic Jason, Antonio Huerta, Jean-François Dubé. Continuum damage modelling in geomechanics. Felix Darve, Ioannis Vardoulakis. Degradations and Instabilities in Geomaterials, 461, Springer, pp.77-105, 2004, International Centre for Mechanical Sciences book series (CISM), 978-3-211-21936-2. 10.1007/978-3-7091-2768-1_3. hal-01570542

HAL Id: hal-01570542

https://hal.science/hal-01570542

Submitted on 31 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Continuum damage modelling in geomechanics

Gilles Pijaudier-Cabot¹, Ludovic Jason¹*, Antonio Huerta², and Jean-François Dubé³

¹R&DO – GeM, Ecole Centrale de Nantes, France

²Laboratori de Càlcul Numèric, Universitat Politècnica de Catalunya, Barcelona, Spain

³LMGC/UMR5508, Université de Montpellier II, France

*also at EDF R&D, Clamart, France

**Abstract:** Continuum damage mechanics is a framework for describing the variations of the elastic properties of a material due to microstructural degradations. This chapter presents the application of this theory to the modelling of concrete. Several constitutive relations are devised, including incremental, explicit, and non local damage models. A general framework for damage induced anisotropy is also presented. Coupled damage and plasticity modelling is discussed. Finally, the issue of the experimental determination of the internal length in non local models is tackled.

## 3.1 Introduction

This chapter is concerned with the presentation of several continuum damage models and their implementation in non linear finite element analyses. Continuum damage means that the mechanical effects of progressive microcracking, void nucleation and growth are represented by a set of state variables which act on the elastic and/or plastic behaviour of the material at the macroscopic level. First, we will deal with the scalar damage model and present different versions of this type of constitutive relations. After having recalled how the local integration of the constitutive relation can be performed in the same spirit as for plasticity-based models, we will deal with more specific applications: cases where the evolution of damage is specified in an integrated form, where damage induced inelastic strains are introduced, where crack closure effects are modelled, and finally where damage is coupled to plasticity. The following section deals with damage induced anisotropy. A general framework for capturing directionality of damage is recalled. Coupling with plasticity and crack closure effects are also discussed. Non local damage is developed in the fourth section. We present integral and gradient models. In particular, a general scheme for the implementation of the gradient model, implemented into the finite element code “Code\_Aster”, is described. The fifth section deals with several computational issues related to damage such as the implementation of the model within secant or tangent algorithms. The chapter concludes on the problem of the calibration of the model parameters in damage models including the internal length. Inverse analysis of size effect tests on geometrically similar specimens is used for the experimental determination of these parameters.
## 3.2 Isotropic Damage Model

When damage is assumed to be isotropic, it is considered that it produces a degradation of the elastic stiffness of the material through a variation of the Young's modulus:

$$ \varepsilon_{ij} = \frac{1 + \nu_0}{E_0(1-d)} \sigma_{ij} - \frac{\nu_0}{E_0(1-d)} [\sigma_{kk} \delta_{ij}] \quad (1) $$

where $\nu_0$, $E_0$ and $\delta_{ij}$ are the Poisson's ratio and Young's modulus of the undamaged isotropic material and the Kronecker symbol respectively. The elastic (i.e. free) energy per unit mass of material is:

$$ \rho\psi = \frac{1}{2}(1-d)\varepsilon_{ij}C_{ijkl}^{0}\varepsilon_{kl} \quad (2) $$

where $C^0$ is the stiffness tensor of the undamaged material. This energy is assumed to be the state potential. The damage energy release rate is defined as the variable associated to the damage state variable in the state potential:

$$ Y = -\rho \frac{\partial \psi}{\partial d} = \frac{1}{2} \varepsilon_{ij} C_{ijkl}^{0} \varepsilon_{kl} \quad (3) $$

with the rate of dissipated energy:

$$ \dot{\phi} = -\frac{\partial \rho \psi}{\partial d} \dot{d} \quad (4) $$

For an isotropic damage model, this equation reduces to:

$$ \dot{\phi} = Y\dot{d} \quad (5) $$

The damage energy release rate is always positive and thus the rate of damage must be positive in order to comply with the second principle of thermodynamics.

The rules governing the evolution of damage can be defined following the same principles as those in elasto-plasticity. In fact, the general formalism is that of generalised standard materials where the loading function is assumed to be a pseudo potential of dissipation (Lemaitre, 1992). The evolution of damage requires the definition of a loading function:

$$ f(Y, d) = Y - Y_0 - Z \quad (6) $$

where $Y_0$ is a parameter which defines the threshold of damage, and Z is the hardening-softening controlling variable. The evolution law is prescribed according to the normality rule:
$$ \dot{d} = \dot{\lambda} \frac{\partial f}{\partial Y} \qquad (7) $$

with the Kuhn-Tucker conditions $\dot{\lambda} \ge 0, f \le 0$, and $\dot{\lambda}f = 0$. The evolutionary equations are completed by the definition of a hardening function:

$$ \dot{\lambda} = H \dot{Z} \qquad (8) $$

where $H$ can be regarded as the equivalent of a hardening-softening modulus in generalised plasticity. In most applications, this modulus is not constant. The major difference with plasticity is that the loading function and the evolution equation for damage are provided in an explicit form because these quantities are function of the total strain.

### 3.2.1 Damage Model with Crack Closure

Crack closure effects are of importance when the material is subjected to alternated loads, for instance reinforced concrete structures subjected to earthquakes. The fundamental property of the scalar damage model is that the time derivative of damage is constrained to be positive or zero (second principle of thermodynamics). During load reversals, however, micro cracks are closing progressively and the tangent stiffness of the material should increase while damage is constant (Fig. 1).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic response of a material subjected to uniaxial compression after being damaged in tension.</figcaption></figure>

Within isotropic damage modelling, one solution is to introduce two damage scalars ($d_t, d_c$), instead of one, in order to separate the mechanical effect of micro cracking as a function of the
type of loading. La Borderie (1991) developed the corresponding constitutive relations and applied it to concrete. The elastic energy is:

$$ \rho\psi = \frac{1}{2} \left[ \frac{1}{E_0(1-d_t)} \langle \sigma \rangle_{ij}^{+} \langle \sigma \rangle_{ij}^{+} + \frac{1}{E_0(1-d_c)} \langle \sigma \rangle_{ij}^{-} \langle \sigma \rangle_{ij}^{-} + \frac{\nu_0}{E_0} (\sigma_{ij}\sigma_{ij} - \sigma_{kk}^2) \right] + \frac{\beta_1 d_t}{E_0(1-d_t)} f(\sigma) + \frac{\beta_2 d_c}{E_0(1-d_c)} \sigma_{kk} \quad (9) $$

where $\langle \sigma \rangle^+$ is the positive part of the stress, i.e. the stress tensor expressed in its principal directions with the positive principal stresses only, and $\langle \sigma \rangle^- = \sigma - \langle \sigma \rangle^+$.

Inelastic strains due to the growth of damage are also taken into account in this model with a function $f(\sigma)$ introduced in order to represent the progressive vanishing of inelastic strains due to damage in tension when the material is subjected to compression:

$$ \begin{aligned} f(\sigma) &= \sigma_{kk} \quad \text{when } \sigma_{kk} \in ]0, \infty] \quad f(\sigma) \\ f(\sigma) &= (1 + \frac{\sigma_{kk}}{2\sigma_c})\sigma_{kk} \quad \text{when } \sigma_{kk} \in ]-\sigma_c, 0] \\ f(\sigma) &= -\frac{\sigma_c}{2}\sigma_{kk} \quad \text{when } \sigma_{kk} \in [-\infty, -\sigma_c] \end{aligned} \quad (10) $$

where $\sigma_c$ is the crack closure stress (material parameter). The stress-strain relation reads:

$$ \begin{aligned} \varepsilon_{ij} = \frac{1}{E_0(1-d_t)} \langle \sigma \rangle_{ij}^{+} + \frac{1}{E_0(1-d_c)} \langle \sigma \rangle_{ij}^{-} \\ + \frac{\nu_0}{E_0} [\sigma_{ij} - \sigma_{kk} \delta_{ij}] + \frac{\beta_1 d_t}{(1-d_t)} \frac{\partial f(\sigma)}{\partial \sigma_{ij}} + \frac{\beta_2 d_c}{(1-d_c)} \delta_{ij} \end{aligned} \quad (11) $$

Since the model has two damage variables, the evolutionary equations include two yield functions. These functions are constructed on the same principles as the function used for the one scalar model, except that there are two energy release rates associated to the two damage variables respectively:

$$ Y_t = -\rho \frac{\partial \psi}{\partial d_t}, \quad Y_c = -\rho \frac{\partial \psi}{\partial d_c} \quad (12) $$

The two loading functions are:

$$ f_t(Y_t, d_t) = Y_t - Y_{0t} - Z_t, \quad f_c(Y_c, d_c) = Y_c - Y_{0c} - Z_c \quad (13) $$

and an associated model can be constructed with two independent hardening functions similar to Eq. (8).
Note that this constitutive relation cannot be integrated explicitly. The stress decomposition into positive and negative parts is not known a priori. Therefore, the conditions of evolution of damage are not known for each integration step. Same as in classical plasticity, the integration of the constitutive relations is performed with an iterative scheme: the equations of the system are the state laws and the evolution laws recast as follows in order to obtain a well-conditioned algebraic system:

$$
\left\{
\begin{array}{l}
\displaystyle \epsilon_{ij} = \frac{1}{E_0(1-d_t)} \langle \sigma \rangle_{ij}^{+} + \frac{1}{E_0(1-d_c)} \langle \sigma \rangle_{ij}^{-} \\
\displaystyle + \frac{\nu_0}{E_0} [\sigma_{ij} - \sigma_{kk} \delta_{ij}] + \frac{\beta_1 d_t}{(1-d_t)} \frac{\partial f(\sigma)}{\partial \sigma_{ij}} + \frac{\beta_2 d_c}{(1-d_c)} \delta_{ij} \\
\\
\displaystyle Y_t - Z_t = 0 \\
\\
\displaystyle Y_c - Z_c = 0 \\
\\
\displaystyle \frac{1}{(1-d_t)} = F_t(Y_t, Z_t) \\
\\
\displaystyle \frac{1}{(1-d_c)} = F_c(Y_c, Z_c)
\end{array}
\right.
\qquad (14)
$$

Here, $F_t, F_c$ are calculated from the integrated equation of evolution of damage instead of an incremental one. It does not change the general picture of the integration process, however. The system is written in a matrix form:

$$
[J]\vec{x} = 0 \quad \text{with} \quad \vec{x}^T = \begin{Bmatrix} \sigma & Y_t & Y_c & \frac{1}{1-d_t} & \frac{1}{1-d_c} \end{Bmatrix} \qquad (15)
$$

where matrix [J] can be computed and programmed using formal calculus. A full Newton algorithm is implemented in order to solve this non-linear system. Note that there are two sorts of non-linearities: the first one is due to the evolution of damage (i.e. whether damage grows or not), the second is due to crack closure effects. For constant values of the damage variables, the material is orthotropic in the coordinate system of the principal stresses.

3.2.2 Integrated Damage Model

As we saw above, the evolution of damage is very often related to the state of strain.
Moreover, it is defined in an explicit, integrated way, which is easier to handle. The damage
loading function in Eq. (6) can be redefined as:

$$
f(\tilde{\epsilon}, \kappa) = \tilde{\epsilon} - \kappa \qquad (16)
$$
where $\tilde{\epsilon}$ is a positive equivalent measure of strain and $\kappa$ is a threshold value. The equation $f=0$ represents a loading surface in strain space. For the uniaxial tensile case, the equivalent uniaxial strain in Eq. (16) is straightforward. It is the axial strain. However, for general states of stress, damage evolution should be related to some scalar quantity, function of the state of strain. There are, to this regard, several proposals. For example, an appropriate definition for metals is rooted in the elastic stored energy (Peerlings et al. 1998):

$$ \tilde{\epsilon} = \sqrt{\frac{1}{E} \epsilon_{ij} C_{ijkl} \epsilon_{kl}} \qquad (17) $$

which is depicted in Figure 2. It is equivalent to the integrated version of the constitutive relation presented in section 2.1. For concrete, Mazars (1984) proposed the following form:

$$ \tilde{\epsilon} = \sqrt{\sum (<\epsilon_i>_+)^2} \qquad (18) $$

where $\epsilon_i$ are the principal strains. A third possibility, which is also mentioned by Peerlings, is the modified von Mises definition:

$$ \tilde{\epsilon} = \frac{k-1}{2k(1-2\nu)} I_1 + \frac{1}{2k} \sqrt{\frac{(k-1)^2}{(1-2\nu)^2} I_1^2 - \frac{12k}{(1+\nu)^2} J_2} \qquad (19) $$

$I_1$ and $J_2$ are the first invariant of the strain tensor and the second invariant of the deviatoric strain tensor respectively. Only the modified von Mises criterion leads to a new material parameter, namely the factor $k$. The parameter $k$ is the ratio between uniaxial compressive and uniaxial tensile strengths. These criteria are plotted in figure 2 (with $k = 10$). The loading surfaces (Eqs. 17-19) are closed contours around the origin. The dashed lines represent the constant uniaxial compression and uniaxial tension stress paths.

The evolution of damage has the same form as in the previous sections:

$$ \begin{cases} \text{if } f(\tilde{\epsilon}, \kappa) = 0 \text{ and } \dot{f}(\tilde{\epsilon}, \kappa) = 0 \text{ then } \begin{cases} d = h(\kappa) \\ \kappa = \tilde{\epsilon} \end{cases} \text{ where } \dot{d} \ge 0 \\ \text{otherwise } \begin{cases} \dot{d} = 0 \\ \dot{\kappa} = 0 \end{cases} \end{cases} \qquad (20) $$

The function $h(\kappa)$ is specific, depending on different models in the literature. For tension only, exponential softening can be used:

$$ h(\kappa) = 1 - \frac{\kappa_0}{\kappa} (1 - \alpha + \alpha e^{-\eta(\kappa - \kappa_0)}) \qquad (21) $$
where $κ_0, α, η$ are model parameters.

In order to capture the differences of mechanical responses of the material in tension and in compression, Mazars proposed to split the damage variable into two parts and used the equivalent strain defined in Eq. (18):

$$d = \alpha_t d_t + \alpha_c d_c \qquad (22)$$

where $d_t$ and $d_c$ are the damage variables in tension and compression, respectively.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Contour plots for $\tilde{\epsilon}$ for the elastic stored energy (top left), Mazars definition (top right), and the modified von Mises expression (bottom), after Peerlings et al. (1998).</figcaption></figure>

They are combined with the weight coefficients $α_t$ and $α_c$ defined as functions of the principal values of the strains $\epsilon_{ij}^t$ and $\epsilon_{ij}^c$, due to positive and negative stresses (see Mazars, 1984).

$$\epsilon_{ij}^t = (1-d) C_{ijkl}^{-1} \sigma_{kl}^t, \quad \epsilon_{ij}^c = (1-d) C_{ijkl}^{-1} \sigma_{kl}^c \qquad (23)$$

$$\alpha_t = \sum_{i=1}^{3} \left( \frac{\langle \epsilon_i^t \rangle \langle \epsilon_i \rangle_+}{\tilde{\epsilon}^2} \right)^\beta, \quad \alpha_c = \sum_{i=1}^{3} \left( \frac{\langle \epsilon_i^c \rangle \langle \epsilon_i \rangle_+}{\tilde{\epsilon}^2} \right)^\beta \qquad (24)$$
In uniaxial tension $\alpha_t=1$ and $\alpha_c=0$. In uniaxial compression $\alpha_c=1$ and $\alpha_t=0$. Hence, $d_t$ and $d_c$ can be obtained separately from uniaxial tests. The purpose of exponent $\beta$ is to reduce the effect of damage on the response of the material under shear compared to tension (Pijaudier-Cabot et al. 1991). The evolution of damage is provided in an integrated form, as a function of the variable $\kappa$ (Mazars, 1984):

$$ 
\begin{aligned}
d_t &= 1 - \frac{\kappa_0(1 - A_t)}{\kappa} - \frac{A_t}{\exp(B_t(\kappa - \kappa_0))} \\
d_c &= 1 - \frac{\kappa_0(1 - A_c)}{\kappa} - \frac{A_c}{\exp(B_c(\kappa - \kappa_0))}
\end{aligned}
\qquad (25) 
$$

<figure><img src="image_4.png" /><figcaption>Figure 3. Uniaxial response of the model by Mazars (1984).</figcaption></figure>

Figure 3 shows the uniaxial response of the model in tension and compression with the following parameters: $E_0=30000$ MPa, $\nu_0=0.2$, $\kappa_0=0.0001$, $A_t=1$, $B_t=15000$, $A_c=1.2$, $B_c=1500$, $\beta=1$.

### 3.2.3 Damage coupled to plasticity

Elastic damage models or elastic plastic laws are not totally sufficient to correctly capture the constitutive behavior of concrete. In some cases, the calculation of the damage variable is a key point. For instance, it can become essential when coupled effects are considered (coupling between damage and permeability, damage and porosity ...). Picandet et al. (2001) observed in experiments that there is a relationship between damage and the gas permeability of the material. Damage, measured according to the unloading slope during cyclic compressive loading is related to the variation of permeability, independently from the type of concrete tested (normal concrete, high strength concrete, fiber-reinforced concrete). Thus, the capability of the constitutive model to capture the unloading behavior is essential if a proper evaluation of the permeability is required.
An elastic damage model is not appropriate as irreversible strains cannot be captured: a zero stress corresponds to a zero strain and the value of the damage is thus overestimated (figure 4a). An elastic plastic relation is not adapted (even with softening, see for example Grassl et al, 2002) as the unloading curve always follows the same elastic slope (figure 4b). An appealing alternative consists in combining these two approaches into an elastic-plastic-damage law. The softening behavior and the decrease in the elastic modulus are so reproduced by the damage part while the plasticity effect accounts for the irreversible strains. With this formulation, experimental unloading can be correctly simulated (figure 4c).

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Unloading response of elastic damage, elastic plastic and elastic plastic damage models.</figcaption></figure>

The plastic model is governed by the following set of equations:

$$ 
\begin{aligned}
d\varepsilon_{ij} &= d\varepsilon_{ij}^e + d\varepsilon_{ij}^p \\
\sigma_{ij}^t &= C_{ijkl}^0 \varepsilon_{kl}^e \\
d\varepsilon_{ij}^p &= d\lambda \cdot m_{ij}(\sigma^t, k_h) \\
dk_h &= d\lambda \cdot h(\sigma^t, k_h)
\end{aligned} 
\qquad (26) $$

with the effective stress defined as

$$ \sigma_{ij}^t = \frac{\sigma_{ij}}{(1-d)} \qquad (27) $$
where $\epsilon$, $\epsilon^e$ and $\epsilon^p$ are respectively the total, elastic and plastic strains, $\lambda$ the plastic multiplier and $k_h$ the hardening parameter ($0 \le k_h \le 1$). $m$ and $h$ are the flow vectors defined by:

$$m_{ij} = \frac{\partial f}{\partial \sigma_{ij}^t} \\ h = \frac{\sqrt{\frac{2}{3} \frac{\partial f}{\partial \sigma_{ij}^t} \cdot \frac{\partial f}{\partial \sigma_{ij}^t}}}{\zeta(\sigma^t)} \qquad (28)$$

$f$ is the yield surface and $\zeta$ a function of the stress (Etse and Willam, 1994).

$$f = \bar{\rho}^2 (J_2^t) - \frac{k(k_h, I_1^t) \bar{\rho}_c^2 (I_1^t)}{r^2(\theta, I_1^t)} \qquad (29)$$

with $I_1^t$ and $J_2^t$ the effective stress invariant and $\theta$ the Lode angle function of the second and third stress invariant and ranging from [-$\pi$/6; $\pi$/6]. The evolution of the plastic multiplier is finally given by:

$$f(\sigma^t, k_h) \le 0, \quad d\lambda \ge 0, \quad f(\sigma^t, k_h) \cdot d\lambda = 0 \qquad (30)$$

This local problem is solved with an iterative procedure associated to a closest point projection algorithm (see Perez-Foguet et al, 2000). Figure 5 shows the evolution of the yield surface with an increasing hardening parameter $k_h$ for simple compression $\bar{\rho} = K_1 \sqrt{J_2^t}$, $\bar{\xi} = K_2 I_1^t$, where $K_1$ and $K_2$ are two positive constants).

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Plastic yield surface. Evolution with the hardening parameter for simple compression (Lode = $\pi$/6).</figcaption></figure>
As the hardening parameter has a limited value of 1, once it has reached this critical level, hardening is not allowed any more and the yield surface becomes a failure one (constant effective stress).

Figure 6 gives the axial stress – strain curve for a loading in tension. As the development of damage is predominant during simple tension tests, the damage and plastic-damage models are almost equivalent.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Stress strain curve for simple tension test. Elastic plastic damage formulation.</figcaption></figure>

Cyclic compression is a second elementary test. Experimental results are taken from (Sinha et al, 1964). Figure 7 illustrates the numerical response for simple damage law (without any plasticity). With this type of relation, a zero stress corresponds to a zero strain.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Cyclic compression test. Elastic damage simulation.</figcaption></figure>
The numerical response of the elastic plastic damage model is given in figure 8. Damage simulates the global softening behaviour of concrete while the plastic part reproduces quantitatively the evolution of the irreversible strains. Experimental and numerical unloading slopes are thus similar, contrary to the simple damage formulation response. If this difference could seem negligible, it is in fact essential if a correct value of the damage needs to be captured. The elastic damage model overestimates *d* whereas the full plastic damage constitutive law provides more acceptable results.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Cyclic compression test. Elastic plastic damage formulation.</figcaption></figure>

Figure 9 and 10 illustrates the differences between the two approaches in terms of volumetric behaviors.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Cyclic compression test. Volumetric behaviour for elastic damage model.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Cyclic compression test. Volumetric behaviour for the plastic damage model.</figcaption></figure>

While, with the simple damage law, the volumetric strains keep negative values, the introduction of plasticity exhibits a change in the volumetric response from contractant (negative volumetric strains) to dilatant, a phenomenon which is experimentally observed (see Sfer et al, 2002 for example).

The introduction of plasticity associated with the development of damage plays thus a key role in the numerical simulation of a cyclic compression test. Moreover, the volumetric response, that was totally misevaluated by the elastic damage model, is correctly simulated with the full formulation. In figure 11, numerical results are compared with experiment for different levels of confinement pressures in triaxial experiments (P = 0, 1.5, 4.5, 9, 30 and 60 MPa). Simulations and experiment give similar results. The peak position is quantitatively well reproduced (except for 1.5 MPa).

The global evolution is also correct: the maximum of the axial stress increases with the pressure and the softening part is less and less significant. When the initial hydrostatic pressure takes higher values, the damage part of the model plays a minor role and plasticity effect becomes predominant.

## 3.3 Damage Induced Anisotropy

Microcracking is usually geometrically oriented as a result of the loading history on the material. In tension, microcracks are perpendicular to the tensile stress direction, in compression microcracks open parallel to the compressive stress direction. Although a scalar damage model, which does not account for directionality of damage, might be a sufficient approximation in usual applications, i.e. when tensile failure is expected with a quasi-radial loading path, damage induced anisotropy is required for more complex loading histories. The influence of crack closure is needed in the case of alternated loads: microcracks may close and
the effect of damage on the material stiffness disappears. Finally, plastic strains need to be incorporated in the formulation when the material unloads in compression. The following section describes a constitutive relation based on elasto-plastic damage which addresses these issues (Fichant et al. 1999).

<figure><img src="image_12.png" /><figcaption>Figure 11. Triaxial test with increasing confinement. Axial stress - strain curves for low hydrostatic pressures. Straight lines (black markers) correspond to simulation, dotted lines (white markers) to experiment (on top (a) 0 MPa, (b) 1.5 MPa, (c) 4.5 MPa, (d) 9 MPa).</figcaption></figure>

### 3.3.1 Principle

The model is based on the approximation of the relationship between the overall stress (latter simply denoted as stress) and the effective stress in the material defined here by the equation:

$$ \sigma_{ij}^{t} = C_{ijkl}^{0} \varepsilon_{kl}^{e} \quad \text{or} \quad \sigma_{ij}^{t} = C_{ijkl}^{0} (C^{\text{damaged}})_{klmn}^{-1} \sigma_{mn} \qquad (31) $$

where $\sigma_{ij}^{t}$ is the effective stress component, $\varepsilon_{kl}^{e}$ is the elastic strain, and $C_{ijkl}^{\text{damaged}}$ is the stiffness of the damaged material. We define the relationship between the stress and the effective stress along a finite set of directions of unit vectors $\vec{n}$ at each material point:
$$ \sigma = (1-d(n))n_i \sigma_{ij}^t n_j, \quad \tau = (1-d(n)) \sqrt{\sum_{i=1}^{3} (\sigma_{ij}^t n_j - (n_k \sigma_{kl} n_l) n_i)^2} \quad (32) $$

$\sigma$ and $\tau$ are the normal and tangential components of the stress vector respectively. $d(n)$ is a scalar valued quantity which introduces the effect of damage in each direction $\vec{n}$.

The basis of the model is the numerical interpolation of $d(n)$ (called damage surface) which is approximated by its knowledge over a finite set of directions. The stress is solution of the virtual work equation:

$$ \begin{gathered} \text{find } \sigma_{ij} \text{ such that } \forall \varepsilon_{ij}^* \\ \frac{4\pi}{3} \sigma_{ij} \varepsilon_{ij}^* = \int_S \left( \left[ (1-d(n))n_k \sigma_{kl}^t n_l n_i + (1-d(n))(\sigma_{ij}^t n_j - n_k \sigma_{kl}^t n_l n_i) \right] \varepsilon_{ij}^* n_j \right) d\Omega \end{gathered} \quad (33) $$

Depending on the interpolation of the damage variable $d(n)$, several forms of damage induced anisotropy can be obtained. Here, it is defined by three scalars in three mutually orthogonal directions. It is the simplest approximation that yields anisotropy of the damaged stiffness of the material. The material is orthotropic with a possibility of rotation of the principal axes of orthotropy.

The stiffness degradation occurs mainly for tensile loads. Hence, the evolution of damage will be indexed on tensile strains. The evolution of damage is controlled by a loading surface $f$, which is similar to Eq. (6):

$$ f(n) = n_i \varepsilon_{ij}^e n_j - \varepsilon_d - \chi(n) \quad (34) $$

$\chi$ is an hardening softening variable which is interpolated in the same fashion as the damage surface. The initial threshold of damage is $\varepsilon_d$. The evolution of the damage surface is defined by an evolution equation inspired from that of an isotropic model:

$$ \begin{gathered} \text{if } f(n^*) = 0 \text{ and } n_i^* d\varepsilon_{ij}^e n_j^* > 0 \\ \text{then} \quad \left\{ \begin{array}{l} dd(n^*) = \left[ \frac{\varepsilon_d (1 + a(n_i^* \varepsilon_{ij}^e n_j^*))}{(n_i^* \varepsilon_{ij}^e n_j^*)^2} \exp(-a(n_i^* \varepsilon_{ij}^e n_j^* - \varepsilon_d)) \right] n_i^* d\varepsilon_{ij}^e n_j^* \\ d\chi(n^*) = n_i^* d\varepsilon_{ij}^e n_j^* \end{array} \right. \\ \text{else } dd(\vec{n}^*) = 0, d\chi(n^*) = 0 \end{gathered} \quad (35) $$

The model parameters are $\varepsilon_d$ and $a$. Note that the vectors $\vec{n}^*$ are the three principal directions of the incremental strains whenever damage grows. After an incremental growth of damage, the new damage surface is the sum of two ellipsoidal surfaces: the one corresponding to the initial damage surface, and the ellipsoid corresponding to the incremental growth of damage.
the new damage surface is the sum of two ellipsoidal surfaces: the one corresponding to the initial damage surface, and the ellipsoid corresponding to the incremental growth of damage.

### 3.3.2 Coupling with plasticity

Same as for the isotropic model, we decompose the strain increment in an elastic and plastic one. The evolution of the plastic strain is controlled by a yield function which is expressed in term of the effective stress in the undamaged material. We have implemented in this anisotropic damage model the yield function due to Nadai (1950). It is the combination of two Drucker-Prager functions $F_1$ and $F_2$ with the same hardening evolution:

$$F_i = \sqrt{\frac{2}{3} J_2^t + A_i \frac{I_1^t}{3} - B_i w} \quad (36)$$

where $J_2^t$ and $I_1^t$ are the second invariant of the deviatoric effective stress and the first invariant of the effective stress respectively. $w$ is the hardening variable and $(A_i, B_i)$ are four parameters ($i=1,2$) which were originally related to the ratios of the tensile strength to the compressive strength denoted $\gamma$ and of the biaxial compressive strength to the uniaxial strength denoted $\beta$:

$$A_1 = \sqrt{2} \frac{1-\gamma}{1+\gamma}, \quad A_2 = \sqrt{2} \frac{\beta-1}{2\beta-1}, \quad B_1 = 2\sqrt{2} \frac{\gamma}{1+\gamma}, \quad B_2 = \sqrt{2} \frac{\beta}{2\beta-1} \quad (37)$$

These two ratios are constant in the model: $\beta = 1.16$ and $\gamma = 0.4$. The evolution of the plastic strains is associated to these surfaces. The hardening rule is given by:

$$w = q p^r + w_0 \quad (38)$$

where $q$ and $r$ are model parameters, $w_0$ defines the initial reversible domain in the stress space, and $p$ is the effective plastic strain.

### 3.3.3 Crack closure effects

A decomposition of the stress tensor into a positive and negative part is introduced again: $\sigma = \langle\sigma\rangle_+ + \langle\sigma\rangle_-$, where $\langle\sigma\rangle_+$, and $\langle\sigma\rangle_-$ are the positive and negative parts of the stress tensor. The relationship between the stress and the effective stress defined in Eq. (31) of the model is modified:

$$\sigma_{ij} n_j = (1-d(n)) \langle \sigma \rangle_{+ij} n_j + (1-d_c(n)) \langle \sigma \rangle_{-ij} n_j. \quad (39)$$

$d_c(n)$ is a new damage surface which describes the influence of damage on the response of the material in compression. Since this new variable refers to the same physical state of
degradation as in tension, $d_c(n)$ is directly deduced from $d(n)$. It is defined by the same interpolation as $d(n)$ and along each principal direction $i$, we have the relation:

$$d_c^i = \left( \frac{d_j(1 - \delta_{ij})}{2} \right)^\alpha, \quad i \in [1, 3] \qquad (38)$$

$\alpha$ is a model parameter. The constitutive relations contain 6 parameters in addition to the Young's modulus of the material and the Poisson's ratio. Their determination benefits from the fact that in tension, plasticity is negligible. Once the evolution of damage in tension has been fitted, the remaining parameters are fitted from a compression test. Figure 12 shows a typical uniaxial compression-tension response of the model corresponding to a concrete with a tensile strength of 3 MPa and a compressive strength of 40 MPa.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Uniaxial tension-compression response of the anisotropic model (longitudinal (1), transverse (2) and volumetric (v) strains as functions of the compressive stress)</figcaption></figure>

## 3.4 Non Local Damage

We turn now attention to structural computations with the isotropic damage model. It is well known that when strain softening is encountered (see. Fig. 3 for instance), some difficulties arise.

Physically, after microcracking has taken place in the first part of the inelastic behavior of concrete, a macroscopic mechanism develops in the form of a damage band. The strains tend
to localize in a specific area of finite dimension (macrocrack) but the model predicts this
dimension to be zero, with a vanishing energy dissipation at failure. The stress response of a
material point cannot be described locally but has also to consider its neighbourhood in order
to take into account microcracks interaction (Askes, 2000). Mathematically (Benallal et al,
1993, Peerlings et al, 1996), when a material point enters the softening regime, the type of the
governing equations changes from elliptic to hyperbolic. The mathematical description
becomes ill-posed since the initial and boundary conditions that are meaningful for the elliptic
problem are usually not relevant for the hyperbolic part. Numerically, the local formulation of
the problem yields a response which is dependant on the smoothness and the orientation of the
mesh. Let us consider for example a bar with a variable section under tensile loading whose
boundary conditions are given in figure 13. In the damage distribution in figure 14 (a), damage
localises only in the most loaded element. If we consider the extreme case of an infinitely fine
mesh, this geometrical localisation in a single element yields a fracture without energy
dissipation. The local formulation thus exhibits a mesh dependency problem and a
regularisation technique has to be added in order to achieve a mesh independent result and
damage localisation over several finite elements (Fig. 14b).

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Boundary conditions</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Damage distribution (a) local simulation (b) non local simulation. Clear zones correspond to damaged one.</figcaption></figure>

It is now established that non locality, in a gradient or integral format, is mandatory for a
proper, consistent, modelling of fracture (see e.g. de Borst et al. 1993). It avoids the
difficulties encountered upon material softening and strain localisation. Within a single
approach, it encompasses both crack initiation (for which continuum models are very well
fitted) and crack propagation (for which discrete fracture approaches have been developed).
We shall describe in the following integral and gradient non local damage models.
### 3.4.1 Integral Model

Consider for instance the scalar damage model in which evolution of damage is controlled by the equivalent strain $\tilde{\varepsilon}$ introduced by Mazars. The principle of nonlocal continuum models with local strains is to replace the equivalent strain $\varepsilon$ with its average (Pijaudier-Cabot and Bazant, 1987):

$$ \bar{\varepsilon}(x) = \frac{1}{V_r(x)} \int_{\Omega} \psi(s) \tilde{\varepsilon}(s+x) ds \quad \text{with} \quad V_r(x) = \int_{\Omega} \psi(s) ds \qquad (39) $$

where $\Omega$ is the volume of the structure, $V_r(x)$ is the representative volume at point $x$, and $\psi(s)$ is the weight function, for instance:

$$ \psi(s) = \exp\left(-\frac{4\|s\|^2}{l_c^2}\right) \qquad (40) $$

$l_c$ is the internal length of the non local continuum. $\bar{\varepsilon}$ replaces the equivalent strain in the evolution of damage. In particular, the loading function becomes $f(\bar{\varepsilon}, \kappa) = \bar{\varepsilon} - \kappa$. As we will see in section 3.5, it should be noticed that this model is easy to implement in the context of explicit, total strain models. Its extension to plasticity and to implicit incremental relations is awkward. The local tangent stiffness operator relating incremental strains to incremental stresses becomes non symmetric, and more importantly its bandwidth can be very large due to non local interactions. This is one of the reasons why gradient damage models have become popular over the past few years.

### 3.4.2 Gradient damage model

A simple method to transform the above non local model to a gradient model is to expand the equivalent strain into Taylor series truncated for instance to the second order:

$$ \tilde{\varepsilon}(x + s) = \tilde{\varepsilon}(x) + \frac{\partial \tilde{\varepsilon}(x)}{\partial x} s + \frac{\partial^2 \tilde{\varepsilon}(x)}{\partial x^2} \frac{s^2}{2!} + \dots \qquad (41) $$

Substitution in Eq. (39) and integration with respect to variable $s$ yields:

$$ \bar{\varepsilon}(x) = \tilde{\varepsilon}(x) + c^2 \nabla^2 \tilde{\varepsilon}(x) \qquad (42) $$

where $c$ is a parameter which depends on the type of weight function in Eq. (39). Its dimension is $m^2$ and it can be regarded as the square of an internal length. Substitution of the new expression of the non local equivalent strain in the non local damage model presented above yields a gradient damage model. Computationally, this model is still delicate to implement because it requires higher continuity in the interpolation of the displacement field.
This difficulty can be solved if an implicit format of the gradient damage model is used. Eq. (42) is replaced with

$$ \bar{\varepsilon}(x) - c^2 \nabla^2 \bar{\varepsilon}(x) = \tilde{\varepsilon}(x) \qquad (43) $$

Here, the definition of the non local equivalent strain is implicit. It is the solution of a Fredholm equation. As shown by Peerlings et al. (1996), the implicit form is in fact an exact representation of the integral relation devised by Pijaudier-Cabot and Bazant (1987), provided a specific exponential weight function is used.

### 3.4.3 Extension of the Gradient Approach

In the anisotropic damage models described in section 3.3, the evolution of damage is directional. The evolution of damage is also directional in the microplane-based models and in the smeared crack models. The extension of the gradient damage model to anisotropy, or to any strain-based damage model, can be performed in a very systematic way (Godard, 2001, 2003). The Fredholm equation is written for each strain component $\varepsilon_{nn}$:

$$ \bar{\varepsilon}_{nn} - c^2 \nabla^2 (\bar{\varepsilon}_{nn}) = \varepsilon_{nn} \qquad (44) $$

Therefore, a non local (gradient type) strain tensor is computed within each iteration. It is from this tensor that an equivalent strain can be computed, or that directional damage growth can be devised depending on the type of damage model that is implemented. This extended form of a gradient model has been implemented in the general purpose finite element code “Code_Aster” at EDF.

<figure><img src="image_16.png" /><figcaption>Figure 15. Bending beam.</figcaption></figure>
## 3.5 Computer implementation of damage models

Since the major effect of damage is a secant stiffness reduction when microcracking progresses, a non-linear computational procedure based on the secant material stiffness seems quite appropriate. This type of algorithm is also widely used when an integral non local damage model is implemented. The major reason is that the tangent stiffness operator (in the case of loading) is far from being straightforward because of non locality.

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Comparison between two non local approaches for different values of the internal length parameter.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Distribution of damage. Dark zones correspond to damaged one (D = 1 for black zones, 0<D<1 for grey zones and D = 0 for white zones)</figcaption></figure>

The integral relation due to the non local term is discretised according to the finite element mesh used for the analysis and an usual quadrature rule is employed for its evaluation. In finite element calculations, the weight function is chopped off: the weights that are less than 0.001
are set to zero. The actual volume of integration does not span over the entire volume of the solid and the calculation of the integrals requires less computer time and memory as the number of neighbouring integration points is reduced. It should be kept in mind, however, that the secant stiffness algorithm presents serious difficulties, mainly there is no proof of convergence of such an algorithm. In large scale finite element computations, a standard Newton-Raphson algorithm is highly preferable, as convergence is more robust. Such an algorithm relies on the derivation of the consistent tangent stiffness of the material response.

The implementation of the scalar damage model, coupled to the gradient approach is different. The equilibrium equations and the Fredholm equations are solved as a coupled problem. The non local strain tensor can be discretised according to the same finite element grid as the displacements. For a better robustness, a standard Newton-Raphson algorithm is used with a consistent tangent operator that ensures quadratic convergence, if properly implemented.

In order to illustrate results with the integral and gradient damage models, we consider the classical case of a three point bending beam containing a notch. Two different meshes are used (a COARse 890 element and a MEDIum 1938 element). Figure 15 shows the beam geometry. Figure 16 presents the comparison between the non local integral and regularised approaches for the medium mesh ($s_c$ is equal to the square root of c). Choosing the appropriate non local parameters, the two techniques are almost equivalent. Figure 17 presents the damage distribution in the half beam for the gradient model. One can see that it is mesh independent.

## 3.6 Model calibration

The calibration of the model parameters in the non local (integral or gradient) damage model is facing the difficulty of getting at the same time the parameters involved in the evolution of damage $\kappa_0$, $A_t$, $B_t$, and the internal length $l_c$. The first set is a point-wise piece of information while the influence of the internal length appears in a boundary value problem with non homogeneous strains (by definition, $\tilde{\varepsilon} = \bar{\varepsilon}$). It follows that it is not possible to calibrate the model parameters in the damage evolution law, without considering the computation of a complete boundary value problem up to failure. Hence the model parameter calibration has to rely on an inverse analysis technique.

The model parameters are optimised so that computational results are as close as possible to the experimental results. The optimisation process used here is based on a Levenberg-Marquardt algorithm, which performs the minimisation of the following functional $\mathfrak{I}(\bar{P})$:

$$ \mathfrak{I}(\bar{P}) = \frac{1}{2} \sum_{i=1}^{3} \left( \frac{R^i(\bar{P}) - R^{i \text{exp}}}{R^{i \text{exp}}} \right)^2 \quad (45) $$

where $R_{\text{exp}}^i$ is the experimental response for test $i$ ($i \in [1,3]$), and $R^i(\bar{P})$ is the numerical response corresponding to the vector of model parameters $\bar{P}$. As we will see next, the
calibration cannot be performed on a single set of experimental data, e.g. three point bend experiments. One possibility is to perform the calibration from several experiments, and for instance size effect test data.

### 3.6.1 Size effect tests

Three point bend experiments have been carried out by Le Bellégo et al. (2000) on notched geometrically similar specimens of various height D = 80, 160, and 320 mm, of length L = 4D, and of thickness b = 40 mm kept constant. The length-to-height ratio is L/D = 4, the span-to-height ratio is l/D = 3, and the notch length is D/10 (Fig. 18). These are mortar specimens with a water-cement ratio of 0.4 and a cement-sand ratio 0.46 (all by weight). The maximal sand grain size is dₐ = 3mm. Experimental and computational data will be also interpreted with the help of Bazant's size effect law, in one of its simplest form (Bazant and Planas, 1998). The nominal strength is obtained with the formula:

$$ \sigma = \frac{3 Fl}{2 b(0.9D)^2} \qquad (46) $$

where F is the peak load. The size effect law reads:

$$ \sigma = \frac{Bf'_{t}}{\sqrt{1 + D / D_{0}}} \qquad (47) $$

<figure><img src="image_19.png" /><figcaption>Figure 18. Geometry of the three point bend size effect tests</figcaption></figure>

D₀ is a characteristic size, f'ₜ is the tensile strength of the material, and B is a geometry – related parameter.

Figure 19 shows the size effect law in a log-log plot. It can be seen that the nominal value of the material strength decreases as the size of the beam increases. This is typical of a damage process in which the process zone is size independent (controlled by the internal length) compared to the structure. For very small sizes, the process zone covers all the specimen and a material strength criterion is expected, while for infinitely large specimens the process zone size is negligible and Linear Elastic Fracture Mechanics should apply.
The calibration of the model parameters from inverse analysis relies first on an accurate finite element representation of each size of specimen. Experience shows that the accuracy of the finite element discretisation in the fracture process zone depends on the ratio of the size of the finite element to the internal length. Therefore, geometrically similar finite element meshes should not be used for the calibration. In order to avoid some possible mesh bias, the finite element size in the process zone should be kept constant and small enough compared to the internal length (at least one third of the internal length). Obviously, the internal length is not known a priori and this requirement is checked at the end of the calibration.

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Experimental results of size effect in three point bend tests.</figcaption></figure>

The optimisation technique has to start from an initial set of model parameters obtained from a manual (trial and error) calibration. This manual process depends on the exact form of the constitutive relation. Therefore, it is difficult to devise a general technique. Still it can be carried out in a rather systematic way as follow: first the model parameters are fitted such that the load deflection curve of one size of specimen is correctly described (usually the medium size). As an initial guess, the internal length lies between 3d<sub>a</sub> and 5d<sub>a</sub>, where d<sub>a</sub> is the maximum aggregate size of concrete according to Bazant and Pijaudier-Cabot (1989). The initial damage threshold and the parameters in the equation controlling the damage growth are varied only (in fact, the Young's modulus of the material could be obtained as well from the initial stiffness of the beam). Then, the internal length is adjusted, the other parameters being kept constant, so that the numerical responses for the three specimen sizes are as close as possible to the experimental results. This process requires many calculations, each trial iteration involves three finite element computations. Figure 20 shows such a fit obtained manually.
<figure><img src="image_21.png" /><figcaption>**Figure 20.** Manual fit of the damage parameters from size effect tests.</figcaption></figure>

Because there are experimental errors and because the above calibration process is based on a variation of the internal length which is disconnected to the variation of the other model parameters, there are several sets of model parameters which may provide an acceptable fit of the tests data. By acceptable, we mean that the numerical responses are as much as possible within the extreme curves due to experimental dispersion (or within an arbitrary distance from them). The difference between these sets ought to be magnified. For this purpose, the size effect method is quite useful (see Le Bellégo et al. 2003).

<figure><img src="image_22.png" /><figcaption>**Figure 21.** Interpretation of quasi-identical fits on the size effect plot.</figcaption></figure>

Figure 21 shows how several fits obtained manually with different values of the internal length can be distinguished on the size effect plot.
<figure><img src="image_23.png" /><figcaption>**Figure 22.** Automatic calibration of the model parameters on the middle size experiment for different - fixed - internal lengths.</figcaption></figure>

For each set of model parameters, the constants in the size effect law have been fitted and the plot on Fig. 21 has been constructed. It is clear on this figure that none of them are consistent with the experimental data. While the quality of the fit is quite hard to decide upon on load – deflexion curves, it is much easier to examine it on the size effect plot.

<figure><img src="image_24.png" /><figcaption>**Figure 22.** Automatic calibration of the model parameters. The dotted lines are the model responses and the plain lines are, for each size, the minimum and maximum curves obtained from the experimental data.</figcaption></figure>
In order to exemplify the lack of objectivity of the model calibration (i.e. the lack of uniqueness of the result) when a single size of beam is used as the experimental data set, Fig. 22 shows the results of the automatic calibration on the middle size specimen. Three fixed different values of the internal length have been considered. For each one, the optimisation has been performed yielding different values of the damage parameters. The computed responses are the same, in good agreement with the experimental data. Bending of notched beams is quite often used for parameter identification because it is relatively easy to perform compared to tensile tests that are very sensitive to boundary conditions (Carmeliet, 1999). This figure shows that the model cannot be calibrated from inverse analysis of a single load deflexion curve.

Figure 23 shows the result of the optimisation from the three sizes of specimens simultaneously. Compared to the manual fit in Fig. 20, the present set of model parameter provides computational results that are closer to the experiments. We tried to change the initial set of model parameter in order to investigate its stability. Overall, the final sets did not differ very much. The largest difference was found on the internal length, which lied between 40 mm and 50 mm.

## 3.7 Closure

Continuous damage mechanics is a theory which aims at describing the mechanical effect of cracking and void growth in an elastic material. Isotropic and anisotropic damage models have been presented in this chapter. There is a rather generic technique which permits a rational derivation of damage models, given a level of anisotropy which is defined a priori. In order to be more representative of the concrete mechanical response, damage must be coupled to plasticity. A simple, and classical, technique is to assume that the plastic part of the model is written in term of the effective stress. The same applies to the modelling of rock masses.

Strain-softening, which is typical of concrete and rocks, induces spurious strain localisation. Integral or gradient damage model are solutions to circumvent this problem. A central issue is the calibration of the model parameters, and in particular the experimental determination of the internal length. This can be achieved from size effect test data with the help of inverse analysis of structural responses. It is important, however, to underline that a proper calibration *cannot* be achieved from the fit of a single experiment. Size effect experiments on geometrically similar specimens, or bending and pure tension tests are required.

## Acknowledgements

Financial support from the European Commission through the MAECENAS project (FIS5-2001-00100), through the DIGA RTN network (HPRN-CT-2002-00220), and from the partnership between Electricité de France and the R&DO group through the MECEN project are gratefully acknowledged.
## 3.8 References

Askes H., 2000, Advanced spatial discretisation strategies for localised failure. mesh adaptivity and meshless methods, PhD Thesis, Delft University of Technology, Faculty of Civil Engineering and Geosciences.

Bazant, Z.P., and Pijaudier-Cabot, G., 1989, Measurement of the characteristic length of non local continuum, J. Engrg. Mech. ASCE, 115, 755-767.

Bazant, Z.P. and Planas, J., 1998, Fracture and size effect in concrete and other quasibrITTLE Materials, CRC Press.

Benallal A., Billardon R., Geymonat G., 1993, Bifurcation and rate-independent materials, Bifurcation and stability of dissipative systems, CISM Lecture Notes 327, Springer, 1-44.

de Borst R., Sluys L.J., Muhlhaus H.B., and Pamin J., 1993, Fundamental issues in finite element analyses of strain localisation, Engrg. Comput., 10, 99-121.

Carmeliet J., 1999, Optimal estimation of gradient damage parameters from localisation phenomena in quasi-brittle materials, Int.J. Mech. Cohesive Frict. Mats., 4, 1-16.

Etse, G. Willam, K., 1994, Fracture energy formulation for inelastic behavior of plain concrete. J. Engrg. Mech. ASCE, 120, 1983-2011.

Fichant S., La Borderie C., and Pijaudier-Cabot G., 1999, Isotropic and anisotropic descriptions of damage in concrete structures, Int.J. Mech. Cohesive Frict. Mats., 4, pp. 339-359.

Grassl, P. Lundgren, K. Gylltoft, K., 2002, Concrete in compression: a plasticity theory with a novel hardening law, International Journal of Solids and Structures 39:5205-5223.

Godard V., 2001, Modélisation mécanique non locale des matériaux fragiles: délocalisation de la déformation, Stage de DEA de Mécanique, Université Pierre et Marie Curie.

Godard V., 2003, Modélisation non locale à gradients de déformation, Documentation de référence Code_Aster R5.04.02, www.code-aster.org.

La Borderie C., 1991, Phénomènes unilatéraux dans un matériau endommageable, Thèse de Doctorat de l'Université Paris 6.

Le Bellego C., Dubé J-F., Pijaudier-Cabot G., B. Gérard, 2003, Calibration of nonlocal damage model from size effect tests, European Journal of Mechanics A/Solids, 22, 33-46.

Le Bellego, C., Gérard B., and Pijaudier-Cabot G., 2000, Chemomechanical effects in mortar beams subjected to water hydrolysis, J. Engrg. Mech. ASCE, 126, 266-272.

Lemaitre J., 1992, A course on damage mechanics, Springer-Verlag.

Mazars J., 1984, Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture du béton de structure, Thèse de Doctorat ès Sciences, Université Paris 6.

Nadai A., 1950, Theory of flow and fracture of solids, Vol. 1, second edition, Mc Graw Hill, New York, p. 572.

Peerlings R.H., de Borst R., Brekelmans W.A.M., de Vree J.H.P., 1996, Gradient enhanced damage for quasi-brittle materials, Int. J. Num. Meth. Engrg., 39, 3391-3403.

Peerlings R.H., de Borst R., Brekelmans W.A.M, Geers M.G.D., 1998, Gradient enhanced modelling of concrete fracture, Int.J. Mech. Cohesive Frict. Mats., 3, 323-343.

Picandet, V. Khelidj, A. Bastian, G., 2001, Effect of axial compressive damage on gas permeability of ordinary and high-performance concrete, Cement and Concrete Research, 31, 1525-1532.
Pijaudier-Cabot G. and Bazant Z.P., 1987, Nonlocal damage theory, J. of Engrg. Mech., ASCE, 113, 1512-1533.

Pijaudier-Cabot G., Mazars J., and Pulikowski J., 1991, Steel-concrete bond analysis with nonlocal continuous damage, J. of Structural Engineering, ASCE, 117, 862-882.

Sfer, D., Carol I. Gettu R., Etse, G., 2002, Study of the behavior of concrete under triaxial compression, J. Engrg. Mech. ASCE, 128, 156-163.

Sinha, B.P. Gerstle, K.H. Tulin, L.G., 1964, Stress strain relations for concrete under cyclic loading. Journal of the American concrete institute, 195-211.