# Continuum damage modelling and some computational issues

Gilles Pijaudier-Cabot, Ludovic Jason

► To cite this version:

Gilles Pijaudier-Cabot, Ludovic Jason. Continuum damage modelling and some computational issues. Revue Française de Génie Civil, 2002, 6 (6), pp.991 - 1017. 10.1080/12795119.2002.9692728. hal-01005996

HAL Id: hal-01005996

https://hal.science/hal-01005996

Submitted on 28 Jan 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 and some computational issues

Gilles Pijaudier-Cabot* — Ludovic Jason**

* R&DO

Laboratoire de Génie Civil de Nantes – Saint Nazaire

Ecole Centrale de Nantes and Institut Universitaire de France

1 rue de la Noe

BP 92101

F-44321 Nantes cedex 3

** R&DO, Laboratoire de Génie Civil de Nantes – Saint Nazaire

Ecole Centrale de Nantes and EDF – R&D

1 rue de la Noe

BP 92101

F-44321 Nantes cedex 3, France

**ABSTRACT.** Continuum damage mechanics is a framework for describing the variations of the elastic properties of a material due to microstructural degradations. This paper 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. In the second part of this contribution, computational issues in damage mechanics related to iterative schemes and solution control in non linear computations are considered. The paper concludes with an example of 3D finite element computation of a reinforced concrete beam, as part of a benchmark initiated by Electricité de France.

**KEYWORDS:** Damage, failure analysis, concrete, non linear finite element computations.
# 1. Introduction

This paper is concerned with the presentation of continuum damage models and their implementation in non linear finite element analyses. Continuum damage means that the mechanical effects of progressive micro cracking, 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 several version 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 with more specific applications: cases where the evolution of damage is specified in an integrated form, where damage induced inelastic strains are introduced, and finally where crack closure effects are modelled.

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 is discussed. The relation with smeared crack model is also presented. Non local damage is discussed 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 and convergence properties of the tangent implementation. Solution control techniques whenever bifurcation or snap-back occurs are recalled. The paper concludes on an example of 3D computation of a reinforced concrete bending beam.

# 2. Scalar Damage Model

## 2.1. Incremental 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:

$$ \epsilon_{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)\epsilon_{ij}C_{ijkl}\epsilon_{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} \epsilon_{ij} C^0_{ijkl} \epsilon_{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} \quad [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} \quad [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. For the integration of this constitutive relation, an Euler
forward integration technique may be used. This method consists in an elastic predictor/ damage corrector technique. The elastic predictor is:

$$ \delta\sigma_{ij}^{e} = (1 - d)C_{ijkl}^{0}\delta\epsilon_{kl} \quad [9] $$

and the damage corrector is:

if $\frac{1}{2}(\epsilon + \delta\epsilon) : C^0 : (\epsilon + \delta\epsilon) - Z > 0$ then

$$ \sigma^{\text{new}} = \sigma + \delta\sigma^{e} - \delta\sigma^{\text{dam}}, \quad \delta\sigma^{\text{dam}} = \{\langle H(Z) \rangle(C^0 : \epsilon) \otimes (C^0 : \epsilon)\} \delta\epsilon \quad [10] $$

else $\sigma^{\text{new}} = \sigma + \delta\sigma^{e}$

As opposed to plasticity, there is no risk for divergence when the loading point is mapped back onto the loading surface. The reason for this is that the loading surface is in fact a one-dimensional segment in the Y space. In fact, the finite element implementation of such a model is exactly the same as the finite element implementation of one dimensional plasticity because damage is a scalar, not a second order tensor like the plastic strains in the 3D case.

## 2.2. 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, 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 [11] $$

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^+$.

Anelastic 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 anelastic strains due to damage in tension when the material is subjected to compression:

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

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

$$ \begin{aligned} \epsilon_{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 [13] $$

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 [14] $$

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 [15] $$
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}
\epsilon_{ij} = \frac{1}{E_0(1-d_t)} \langle \sigma \rangle_{ij}^{+} + \frac{1}{E_0(1-d_c)} \langle \sigma \rangle_{ij}^{-} \\
\qquad + \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} \\
Y_t - Z_t = 0 \\
Y_c - Z_c = 0 \\
\frac{1}{(1-d_t)} = F_t(Y_t, Z_t) \\
\frac{1}{(1-d_c)} = F_c(Y_c, Z_c)
\end{array}
\right.
\qquad [16]
$$

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]\bar{x} = 0 \quad \text{with } \bar{x}^T = \begin{Bmatrix} \sigma & Y_t & Y_c & \frac{1}{1-d_t} & \frac{1}{1-d_c} \end{Bmatrix} \qquad [17]
$$

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 nonlinearities: 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.

**2.3. Integrated Damage Model**

In continuum damage mechanics, 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{\varepsilon}, \kappa) = \tilde{\varepsilon} - \kappa \quad [18]$$

where $\tilde{\varepsilon}$ 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 uni-axial tensile case, the equivalent uniaxial strain in Eq. [18] is straightforward. It is the axial strain if the lateral strains are neglected. 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{\varepsilon} = \sqrt{\frac{1}{E} \varepsilon_{ij} C_{ijkl} \varepsilon_{kl}} \quad [19]$$

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

$$\tilde{\varepsilon} = \sqrt{\sum (<\varepsilon_i>)^2} \quad [20]$$

where $\varepsilon_i$ are the principal strains. A third possibility, which is also mentioned by Peerlings, is the modified von Mises definition. It is written as follows:

$$\tilde{\varepsilon} = \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} \quad [21]$$

$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 uni-axial compressive and uni-axial tensile strength. These criteria are plotted in figure 2 (with $k = 10$). The loading surfaces (Eqs. 19-21) 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:

$$\text{if } f(\tilde{\varepsilon}, \kappa) = 0 \text{ and } \dot{f}(\tilde{\varepsilon}, \kappa) = 0 \text{ then } \begin{cases} \dot{d} = h(\kappa) \\ \kappa = \tilde{\varepsilon} \end{cases} \text{ where } \dot{d} \ge 0 \quad [22] \\ \text{otherwise } \begin{cases} \dot{d} = 0 \\ \dot{\kappa} = 0 \end{cases}$$

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

$$h(\kappa) = 1 - \frac{\kappa_0}{\kappa} (1 - \alpha + \alpha e^{-\eta(\kappa - \kappa_0)}) \quad [23]$$
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. [20]:

$$d = α_t d_t + α_c d_c \quad [24]$$

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 function of the principal values of the strains $ε_{ij}^t$ and $ε_{ij}^c$, due to positive and negative stresses (see Mazars, 1984).

$$ε_{ij}^t = (1 - d)C_{ijkl}^{-1}σ_{kl}^t, \quad ε_{ij}^c = (1 - d)C_{ijkl}^{-1}σ_{kl}^c \quad [25]$$

$$α_t = \sum_{i=1}^{3} \left( \frac{\langle ε_i^t \rangle \langle ε_i \rangle_+}{\tilde{\epsilon}^2} \right)^{\beta}, \quad α_c = \sum_{i=1}^{3} \left( \frac{\langle ε_i^c \rangle \langle ε_i \rangle_+}{\tilde{\epsilon}^2} \right)^{\beta} \quad [26]$$
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 evolution of damage is provided in an integrated form, as a function of the variable $\kappa$ (Mazars, 1984):

$$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))} \qquad [27]$$

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

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). Figure 3 shows the uniaxial response of the model in tension and compression with the following parameters: $E_0=30000$ MPa, $v_0=0.2$, $\kappa_0=0.0001$, $A_t=1$, $B_t=15000$, $A_c=1.2$, $B_c=1500$, $\beta=1$.

### 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 are observed
when the material unloads in compression. The following section describes a constitutive relation based of elasto-plastic damage which addresses these issues (Fichant et al. 1999).

### 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 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 [28] $$

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 $\bar{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} \left( \sigma_{ij}^t n_j - (n_k \sigma_{kl} n_l) n_i \right)^2} \qquad [29] $$

$\sigma$ and $\tau$ are the normal and tangential components of the stress vector respectively. $d(n)$ is a scalar valued quantity which introduce the effect of damage in each direction $\bar{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:

$$ \text{find } \sigma_{ij} \text{ such that } \forall \varepsilon_{ij}^* \qquad \frac{4\pi}{3} \sigma_{ij} \varepsilon_{ij}^* = \int_S \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) \cdot \varepsilon_{ij}^* n_j \, d\Omega \qquad [30] $$

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 \epsilon_{ij}^e n_j - \epsilon_d - \chi(n) \quad [31]$$

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

if $f(n^*) = 0$ and $n_i^* d\epsilon_{ij}^e n_j^* > 0$

$$\text{then} \quad \left\{ \begin{array}{l} dd(n^*) = \left[ \frac{\epsilon_d(1 + a(n_i^* \epsilon_{ij}^e n_j^*))}{(n_i^* \epsilon_{ij}^e n_j^*)^2} \exp(-a(n_i^* \epsilon_{ij}^e n_j^* - \epsilon_d)) \right] n_i^* d\epsilon_{ij}^e n_j^* \\ d\chi(n^*) = n_i^* d\epsilon_{ij}^e n_j^* \end{array} \right. \quad [32]$$

else $dd(\bar{n}^*) = 0$, $d\chi(n^*) = 0$

The model parameters are $\epsilon_d$ and $a$. Note that the vectors $\bar{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.

## 3.2. Coupling with plasticity

We decompose the strain increment in an elastic and plastic one:

$$d\epsilon_{ij} = d\epsilon_{ij}^e + d\epsilon_{ij}^p \quad [33]$$

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 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 [34]$$

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 [35]$$

These two ratios will be kept 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 [36]$$

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. 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. [29] of the model is modified:

$$\sigma_{ij} n_j = (1 - d(n)) \langle\sigma\rangle_{+ij}^t n_j + (1 - d_c(n)) \langle\sigma\rangle_{-ij}^t n_j. \quad [37]$$

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Uniaxial tension-compression response of the anisotropic model (longitudinal (1), transverse (2) and volumetric (v) strains as functions of the compressive stress)</figcaption></figure>
$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 4 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.

### 3.4. Damage and smeared crack models

Historically, smeared crack models have been developed for modelling concrete fracture fifteen years before continuous damage models started to become popular. In smeared crack models the directionality of the material decohesion is a fundamental characteristic, which was not included into the first damage models. Therefore, it is natural to relate damage models to smeared crack models (de Borst and Gutierrez, 1999). In the latter model, the material response is defined in a local coordinate system $(n, s)$. Direction $n$ is the normal to the plane in which the greatest positive normal strain denoted as $\epsilon_{nn}$ is found. In the $(n, s)$ system, the secant stress-strain relation is:

$$\sigma_{n,s} = D_{ns}^s \epsilon_{n,s} \qquad [39]$$

with $\sigma_{n,s} = [\sigma_{nn} \ \sigma_{ss} \ \sigma_{ns}]^T$, $\epsilon_{n,s} = [\epsilon_{nn} \ \epsilon_{ss} \ \epsilon_{ns}]^T$, and the secant stiffness $D_{ns}^s$ is

$$D_{ns}^s = \begin{bmatrix} \frac{(1 - d_1)E}{1 - (1 - d_1)\nu^2} & \frac{(1 - d_1)\nu E}{1 - (1 - d_1)\nu^2} & 0 \\ \frac{(1 - d_1)\nu E}{1 - (1 - d_1)\nu^2} & \frac{E}{1 - (1 - d_1)\nu^2} & 0 \\ 0 & 0 & (1 - d_2)G \end{bmatrix} \qquad [40]$$

where $d_1$ and $d_2$ are two damage parameters. $(1 - d_2)$ is the degradation of the shear stiffness $G$ and can be related to the shear retention factor in traditional
smeared crack models. The evolution of these two damage parameters is defined with the help of a loading function in the (n,s) coordinate system: $f(\epsilon_{nn}, \kappa) = \epsilon_{nn} - \kappa$ and the appropriate Kuhn-Tucker conditions. $d_1$ and $d_2$ are functions of the history variable $\kappa$, same as in the other damage models. If we introduce $\phi$ as the angle between the (x,y) coordinate system and the (n,s) coordinate system, the strain and stress components in the two systems can be related as follow:

$$ \epsilon_{n,s} = T_{\epsilon}(\phi)\epsilon_{x,y} \quad \text{and} \quad \sigma_{n,s} = T_{\sigma}(\phi)\sigma_{x,y} \qquad [41] $$

where $T_{\epsilon}(\phi)$ and $T_{\sigma}(\phi)$ are the appropriate transformation matrices. The secant relation in Eq. [42] becomes:

$$ \sigma_{x,y} = T_{\sigma}^{-1}(\phi) D_{ns}^{s} T_{\epsilon}(\phi) \epsilon_{x,y} \qquad [42] $$

This equation, and the loading function $f(\epsilon_{nn}, \kappa)$, incorporate the fixed smeared crack model and the rotating crack model as well. The difference is that in the fixed crack model the angle $\phi$ is constant whereas in the rotating crack model, the coaxiality requirement enforces that the n-direction is always the major principal strain direction.

As we can see, smeared crack models are indeed damage models, with two damage coefficients and some rule upon which the principal directions of damage may change. In the fixed crack model, the principal directions of damage are fixed; in the rotating crack model, they rotate according to the condition of coaxiality. There is a similarity between the above formulation of the smeared crack model and the microplane-based damage model described in section 3. The integral which defines the overall stress in the microplane model, is the same as in the model by Fichant [Eq. 30]. This integral can be transformed in order to arrive to a format that is very similar to the multiple fixed crack model:

$$ \sigma_{ij} = \sum_{n} \omega^{n} \left( \left[ \begin{array}{c} (1 - d(\bar{n})) E_{N} n_{k} \epsilon_{kl} n_{l} n_{i} \\ + (1 - D(\bar{n})) E_{T} (\epsilon_{ij} n_{j} - n_{k} \epsilon_{kl} n_{l} n_{i}) \end{array} \right] \cdot n_{j} \right) \qquad [43] $$

where $\omega^n$ is a weighting factor.

## 4. Non Local Damage

We turn now attention to the non local generalisation of damage models. 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, or the contribution by M. Jirasek in this volume). 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).

## 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 $\tilde{\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 [44] $$

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 [45] $$

$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 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.

## 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 [46] $$

Substitution in Eq. [44] and integration with respect to variable $s$ yields:

$$ \bar{\varepsilon}(x) = \tilde{\varepsilon}(x) + c^2 \nabla^2 \tilde{\varepsilon}(x) \qquad [47] $$
where c is a parameter which depends on the type of weight function in Eq. [47]. 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. [47] is replaced with

$$ \bar{\epsilon}(x) - c^2 \nabla^2 \bar{\epsilon}(x) = \tilde{\epsilon}(x) \quad [48] $$

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 an exponential weight function is used.

### 4.3. Extension of the Gradient Approach

In the anisotropic damage models described in section 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). The Fredholm equation is written for each strain component $\epsilon_{nn}$:

$$ \bar{\epsilon}_{nn} - c^2 \nabla^2 (\bar{\epsilon}_{nn}) = \epsilon_{nn} \quad [49] $$

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.

## 5. Some Computational Issues

### 5.1. Secant Stiffness Algorithm

Since the major effect of damage is a secant stiffness reduction when micro-cracking progresses, a non-linear computational procedure based on the secant material stiffness seems quite appropriate. Total displacement v.s. total force relations can be used at each load increment:

$$ K^{\sec} \bar{u} = \bar{f} \quad [50] $$
$K^{\text{sec}}$ is the secant stiffness matrix of the discrete solid which is analysed:

$$K^{\text{sec}} = \int_V (1-d) B^T C^0 B dv \qquad [51]$$

$$\text{with } \varepsilon = B\bar{u}$$

where $\bar{u}$ is the vector of nodal displacements. During each iteration, a new displacement field is computed. Convergence is obtained when the residual forces are less than a maximum tolerance. These residual forces at iteration n are:

$$\delta \bar{f}_n = K_n^{\text{sec}} \bar{u}_n - \bar{f}, \quad \text{with} \quad K_n^{\text{sec}} = \int_V (1-d_n) B^T C^0 B dv \qquad [52]$$

In most applications, the norm of the local residual forces at each node and the norm of the residual forces vector over the entire structure should be less than the maximum tolerance. These residuals are non dimensionalised (e.g. by dividing by the applied load vector).

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. It is possible to obtain an expression of the rate constitutive relations where we can see that the local tangent operator is non local, i.e. does not depend on the behaviour of a single material point:

$$\dot{\sigma}(x) = \left\{ (1 - d_0) C^0 \dot{\varepsilon}(x) - \left( \frac{\partial h(\kappa)}{\partial \kappa} C^0 : \varepsilon_0(x) \right) \int_V \psi(s) \dot{\varepsilon}(x + s) dv \right\} \qquad [53]$$

One central feature of this model is that the rate of damage and the rate of stress are explicit functions of the rate of strain. Compared to non local plasticity where the incremental strain is split into a plastic and an elastic part which are unknown in advance, the consistency condition is here an integral relation whose kernel is constant and known once for all in an integration step. This property makes the time integration of the constitutive relation more simple than in non local plasticity. It is also a reason why gradient plasticity (de Borst et al. 1993) is generally preferred to non local (integral) plasticity because the consistency condition in non local plasticity is an implicit integral relation which is complex to solve iteratively.

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. Nevertheless, the bandwidth of the tangent stiffness matrix is much larger with the non local model than with the local model. Averaging expands the number of non zero terms simply because non local terms correspond to long range interactions. The tangent stiffness matrix of the finite element model is, however, non symmetric for two reasons sketched on figure 5:

1. Non locality is an interaction that exists only when damage grows. Consider two points A and B sufficiently close to each other so that the non local interaction is strong. Damage grows at point A and remains constant at point B. The growth of damage in A is non local and depends on the incremental strain at point B. The behaviour of point B is local and does not depend on the incremental strain at point A. The influence functions between points A and B are not mutually equal. Therefore the corresponding terms placed symmetrically in the stiffness matrix are not equal.

2. The volume $V_r(x)$ in Eq. [53] is not constant. In particular, this quantity which normalises the averaging decreases near free boundaries. The influence of a point A located near a free surface on a point B located farther from the free surface is not the same as the influence of B on A because the averages are different.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** *Sources of non symmetry of the tangent stiffness matrix*</figcaption></figure>

With the secant algorithm, the non local model is relatively easy to implement since the equilibrium equations remain standard. Their weak form is the same as for the classical, local, damage model.

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.

## 5.2. Consistent Tangent Stiffness and Convergence

In this section, we are going to focus on the scalar damage model, coupled to the gradient approach. The implementation of damage induced anisotropy following the same procedure remains to be performed still.

Consider again the discrete equations of equilibrium in a total displacement format:

$$ \int_V (1-d) B^T C^0 B dv \bar{u} = \bar{F} \quad [54] $$

where $\bar{u}$ and $\bar{F}$ are the unknown nodal displacement and external forces respectively. The tangent stiffness is formally an unknown in this equation as it depends on the displacements (through the unknown distribution of damage). The fundamental principle of the Newton-Raphson method is the construction of a series of approximation of the residual $R(\bar{u})$

$$ R(\bar{u}) = \int_V (1-d) B^T C^0 B dv \bar{u} - \bar{F} \quad [55] $$

which converges to zero for a finite number of iterations, and for a fixed value of the external loads. For this purpose, the residual at iteration $i+1$ is computed from that at iteration $i$ according to a first order Taylor expansion:

$$ R(\bar{u}_{i+1}) = R(\bar{u}_i) + \left( \frac{\partial R}{\partial \bar{u}} \right)_i \delta \bar{u}_i \quad [56] $$

with

$$ \bar{u}_{i+1} = \bar{u}_i + \delta \bar{u}_i \quad [57] $$

The corrective term in the displacement field $\delta \bar{u}_i$ is evaluated as the solution of the system:

$$ \int_V B^T \left( \frac{\partial \sigma}{\partial \epsilon} \right)_i B dv \delta \bar{u}_i = -R(\bar{u}_i) \quad [58] $$

In this equation, the matrix that is on the left hand-side term is the consistent tangent stiffness matrix. Its evaluation requires an exact derivation of the constitutive relation. The advantage is that with this matrix, quadratic convergence of the algorithm is observed.
In the modified Newton-Raphson scheme, the tangent matrix is substituted with the elastic matrix:

$$ \int_V B^T C^0 B dv \delta \bar{u}_i = -R(\bar{u}_i) \quad [59] $$

The price to pay to this simplification is that convergence is not quadratic. It is linear. On the other side, the evaluation of this elastic matrix is less time consuming. Recall that the tangent stiffness matrix needs to be updated at each time step.

In the gradient damage model, 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, with the same interpolation function N. Eq. [58] becomes:

$$ \begin{bmatrix} K_{aa} & K_{ae} \\ K_{ea} & K_{ee} \end{bmatrix} \begin{pmatrix} \delta \bar{u}_i \\ \delta \bar{\epsilon}_i \end{pmatrix} = - \begin{pmatrix} R(\bar{u}_i) \\ K_{ee} \bar{\epsilon}_i - \int_V N^T \epsilon_i dv \end{pmatrix} \quad [60] $$

with the matrices K defined as

$$ \begin{align} K_{aa} &= \int_V (1-d) B^T C^0 B dv, & K_{ae} &= \int_V B^T \frac{\partial \sigma}{\partial \bar{\epsilon}} N dv \nonumber \\ K_{ea} &= \int_V -N^T B dv, & K_{ee} &= \int_V (c^2 B^T B + N^T N) dv \tag{61} \end{align} $$

In order to avoid stress oscillations, the non local strain tensor can be also discretised with interpolation functions that are linear while the interpolation of displacements is quadratic. This is what has been implemented in the FE software Code ASTER. The iterations are carried out until the relative residual forces are less than a given tolerance:

$$ \left\| \frac{R(\bar{u}_i)}{\bar{F}} \right\| \le \text{Tolerance} \quad [62] $$

The difficulty in this algorithm is the proper calculation of the tangent stiffness at the material level. Very often, small errors, or approximations, in the derivation result in a loss of the quadratic convergence of the Newton-Raphson scheme. In particular cases such as the damage model proposed by Mazars and presented Eqs. [20, 22-27], the tangent operator is very difficult to derive properly. The reason is that the differentiation of the factors $\alpha_t, \alpha_c$ is very uneasy. In the calculations presented in this paper, the variation of these coefficients is neglected. Hence quadratic convergence is expected for radial loading only (constant $\alpha_t, \alpha_c$).
Figure 6 shows the evolution of the relative residual during the iterations for several initial states of damage. The computation corresponds to uniaxial compression. The Newton-Raphson scheme [Eq. 58] and the modified scheme Eq. [59] have been compared on this plot. Quadratic convergence yields an evolution of the residual which follows a parabola approximately whereas linear convergence results into a straight line on this plot. It is clear that the modified scheme yields a slowest convergence, resulting into a larger number of iteration for the same tolerance.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution of the relative residuals for the Newton-Raphson and modified Newton-Raphson schemes</figcaption></figure>

### 5.3. Solution Control

Upon localisation of damage and strain, i.e. in the course of the failure process, uniqueness of the equilibrium equations can be lost. Bifurcation can occur in the continuum solution and should be expected for the discrete approximation as well, and stability might be lost too. Therefore, finite element codes must be equipped with solution checkings capable of detecting loss of uniqueness and loss of stability. If bifurcation occurs, the solution should be controlled so that it follows the stable path. If stability is lost, there should be the possibility for switching to indirect displacement control so that the finite element solution can be traced during instability. In the following we shall consider the rate equation of equilibrium of the discrete problem written generically as:

$$ \int_V B^T \dot{\boldsymbol{\sigma}} d\boldsymbol{v} = K \dot{\bar{\boldsymbol{u}}} - \bar{\boldsymbol{F}} \qquad [63] $$
where **K** is the tangent stiffness matrix of the discrete solid. Uniqueness of the solution can be assessed by considering the equations of equilibrium, Eq. [63]. For the sake of simplicity, attention is restricted to the case where the loading is not stationary as bifurcation points and limit points rarely coincide. If the solution of the rate boundary value problem is non unique, the tangent stiffness matrix should be singular :

$$ \det(K) = 0 \qquad [64] $$

However, **K** is not a single valued matrix and depends on the loading conditions. Rigorously, all the possible (loading-unloading) combinations should be investigated. In most applications, **K** is calculated for the incrementally linear comparison solid only. It means that, locally, the tangent modulus computed at each material point corresponds to loading if the consistency condition is met at the considered time step of the calculation. Our experience shows that the loss of uniqueness for the comparison solid occurs before it can be detected for any other loading-unloading configurations.

At a state of equilibrium under dead load, the equilibrium of a discrete system is critical if there exists a kinematically admissible field $\dot{\vec{u}}$ such that:

$$ \dot{\vec{u}}K\dot{\vec{u}} = 0 \qquad [65] $$

Again **K** is assumed to be single valued. This condition is equivalent to:

$$ \begin{cases} \det(\vec{K}) = 0 \\ \vec{u} \text{ orthogonal to } \vec{K\vec{u}} \end{cases} \qquad [66] $$

The loss of stability does not necessarily occur at a bifurcation point because the tangent stiffness operator is not symmetric, and theoretically it may be possible to find a velocity field such that the second condition in Eq. [66] is satisfied. However it is easy to show (Pijaudier-Cabot and Huerta, 1991) that Eq. [66] is equivalent to:

$$ \det(S) = 0 \qquad [67] $$

where **S** is the symmetric part of **K**.

During the calculation and once the conditions for stability or uniqueness are not satisfied, it is necessary to investigate whether there exists a stable response of the structure and what are the different possible paths in the post-bifurcation regime. As the loading progresses, the loss of stability is encountered first. If the criterion for uniqueness is not met, i.e. if **K** is non singular, then the solution is unique and the equilibrium of the structure is critical if the vanishing eigenvector of **S** is collinear to the solution of the problem Eq. [63].

If this is the case, load control or displacement control techniques fail to trace the snap-back portion of the response curve of the structure, and other algorithms
should be implemented to circumvent the singularities beyond critical points. In fact, most of the continuation techniques are based on the introduction of another variable which governs the load factor and induces an extra-equation. Some of the available methods are indirect displacement control (de Borst, 1986) or the well known arc-length control method (Riks, 1979; Crisfield, 1981).

<figure><img src="image_8.png" /><figcaption>**Figure 7. Finite element mesh for three point bending beams. The load is applied at the top-right corner (axis of symmetry) and vertical displacements are fixed at the bottom left corner**</figcaption></figure>

As an example, let us consider a notched bending beam. Figure 7 shows the finite element mesh (half of it due to the symmetry) that has been implemented. Figure 8 shows a comparison of the responses computed for the non local gradient damage model and for the non local integral damage model. In both cases, an indirect control displacement procedure is needed. Indeed, a displacement controlled computation would yield an unstable response due to snap-back.

<figure><img src="image_9.png" /><figcaption>**Figure 8. Non Local computations of notched three point bending beams with the gradient (lcgrad) and integral (lcint) models for several internal lengths**</figcaption></figure>
Finally, the case where the criterion for uniqueness is not met must be investigated. Again we know that there is at least one solution denoted as $\dot{\vec{u}}^*$ to the rate equilibrium problem. Let us call this solution the fundamental solution. The objective now is to perturb the fundamental solution adequately in order to obtain the solution over the stable path. This perturbation is realised by means of the right hand eigenvector $\vec{v}$ associated to the vanishing eigenvalue of the tangent stiffness matrix $K$ because it can be demonstrated (de Borst, 1986 and 1988) that for any arbitrary scalar $\beta$, the velocity field $\dot{\vec{u}}^*+\beta\vec{v}$ is also a solution of the problem:

$$K(\dot{\vec{u}}^*+\beta\vec{v}) = \dot{\vec{F}} \quad [67]$$

For more details on this specific subject, see Pijaudier-Cabot and Huerta (1991) and de Borst (1988).

## 6. Conclusion - Example of Computation

As a conclusion, the 3D computation of a reinforced concrete bending beam is presented in this section. This test case belongs to the benchmark study that has been carried out by EDF (Ghavamian, 1999).

<figure><img src="image_10.png" /><figcaption>**Figure 9.** *Three point bending beam*</figcaption></figure>
The beam and the load system are shown in Fig. 9. The computation has been carried out with a scalar damage model which derives from the damage induced anisotropic model presented in section 3. The equations are the same, except that the interpolation of damage is carried out over a single direction instead of three mutually orhtogonal ones. The finite element mesh is made of brick elements for concrete and elasto-plastic bar elements for the longitudinal reinforcements and for the stirrups. Figure 10 shows the load deflection curve.

<figure><img src="image_11.png" /><figcaption>**Figure 10. Load deflection curve**</figcaption></figure>

Figures 11 and 12 show the distribution of damage in two cross sections of the beam. The first one is longitudinal and the second one is a transverse cross section at mid-span.

<figure><img src="image_12.png" /><figcaption>**Figure 11. Distribution of damage in a longitudinal cross section.** The load is applied at the top-right corner (axis of symmetry) and vertical displacements are fixed at the bottom left corner</figcaption></figure>
On these figures, the dark grey colour corresponds to zero damage and the white colour corresponds to complete damage. Overall, the distribution of damage is very much in agreement with experiments, which show regularly spaced cracking.

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Distribution of damage in a transverse cross section at mid-span</figcaption></figure>

As seen on these figures, 3D effects are captured in addition to the map of regularly spaced vertical cracks and diagonal cracks. The computation has been stopped when the applied vertical displacement reached –1 cm.

It is on this type of 3D finite element computations, which correspond to major industrial applications, that several issues involving computational damage are faced. One can mention the pertinence of the constitutive relation, of course, but also with respect to “non mechanical” issues involved in safety analyses (occurrence of leaks, ageing, ...). The robustness of the finite element implementation is another one, in which the simple question “is a solution expected from the finite element calculation?” is far from being trivial in large scale analyses. Finally, the quality “at large” of finite element simulation of failure still remains an outstanding issue which should be addressed if this type of modelling is intended to be more widely used in engineering practice.

### Acknowledgements

Partnership of Electricité de France with R&DO is gratefully acknowledged. The authors would like to thank Prof. A. Huerta (UPC Barcelona, Spain) for helpful discussions about the implementation of the gradient non local damage model.

## 7. References

* Crisfield M.A. (1981), "A Fast Incremental / Iterative Solution Procedure that Handles Snap-Through", *Comp. and Struct.*, Vol. 13, pp. 55-62.
de Borst R., Sluys L.J., Muhlhaus H.B., and Pamin J. (1993), "Fundamental Issues in Finite Element Analyses of Straain Localisation", *Engrg. Comput.*, Vol. 10, pp. 99-121.

de Borst R. (1988), "Bifurcation in Finite Element Models with Non-Associated Flow Law", *Int. J. for Num. and Analytical Meth. in Geomechanics*, Vol. 12, pp. 99-116.

de Borst R. (1986), Non Linear Analysis of Frictional Materials, Dissertation, TU-Delft.

de Borst R. and Gutierrez M. (1999), "A unified framework for concrete damage and fracture models including size effects", *Int. J. Fracture*, Vol. 95, pp. 261-277.

Fichant S., La Borderie C., and Pijaudier-Cabot G. (1999), "Isotropic and Anisotropic Descriptions of Damage in Concrete Structures", *Int. J. Mechanics of Cohesive Frictional Materials*, Vol. 4, pp. 339-359.

Ghavamian S. (1999), Evaluation tests on models of non-linear behaviour of cracking concrete, using three dimensional modelling, Benchmark EDF/Division R&D, CR - MMN 99/232.

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.

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

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.*, Vol. 39, pp. 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. Mechanics of Cohesive Frictional Materials*, Vol. 3, pp. 323-343.

Pijaudier-Cabot G. and Bazant Z.P. (1987), "Nonlocal Damage Theory", *J. of Engrg. Mech.*, ASCE, Vol. 113, pp.1512-1533.

Pijaudier-Cabot G., Mazars J., and Pulikowski J. (1991), "Steel-Concrete Bond Analysis with Nonlocal Continuous Damage", *J. of Structural Engineering*, ASCE, Vol. 117, No.3, pp. 862-882.

Pijaudier-Cabot G., and Huerta A. (1991), "Finite Element Analysis of Bifurcation in Nonlocal Strain Softening Solids", *Comp. Meth. in Applied Mech. and Engrg.*, Vol. 90, pp. 905-919.

Riks E. (1979), "An Incremental Approach to the Solution of Snapping and Buckling Problems", *Int. J. Solids and Struct.*, Vol. 15, pp. 529-551.