# Damage Models for Concrete

Gilles Pijaudier-Cabot, Jacky Mazars

► To cite this version:

Gilles Pijaudier-Cabot, Jacky Mazars. Damage Models for Concrete. Jean Lemaitre Handbook of Materials Behavior Models, 2, Elsevier, pp.500-512, 2001, Failures of materials, 978-0-12-443341-0. 10.1016/B978-012443341-0/50056-9. hal-01572309

HAL Id: hal-01572309

https://hal.science/hal-01572309

Submitted on 6 Aug 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
Damage Models for Concrete

GILLES PIJAUDIER-CABOT¹ and JACKY MAZARS²

¹Laboratoire de Génie Civil de Nantes Saint-Nazaire, Ecole Centrale de Nantes, BP 92101, 44321 Nantes Cedex 03, France

²LMT-Cachan, ENS de Cachan, Université Paris 6, 61 avenue du Président Wilson, 94235, Cachan Cedex, France
# 1 ISOTROPIC DAMAGE MODEL

## 1.1 VALIDITY

This constitutive relation is valid for standard concrete with a compression strength of 30–40 MPa. Its aim is to capture the response of the material subjected to loading paths in which extension of the material exists (uniaxial tension, uniaxial compression, bending of structural members) [4]. It should not be employed (i) when the material is confined ( triaxial compression) because the damage loading function relies on extension of the material only, (ii) when the loading path is severely nonradial (not yet tested), and (iii) when the material is subjected to alternated loading. In this last case, an enhancement of the relation which takes into account the effect of crack closure is possible. It will be considered in the anisotropic damage model presented in Section 3. Finally, the model provides a mathematically consistent prediction of the response of structures up to the inception of failure due to strain localization. After this point is reached, the nonlocal enhancement of the model presented in Section 2 is required.

## 1.2 BACKGROUND

The influence of microcracking due to external loads is introduced via a single scalar damage variable $d$ ranging from 0 for the undamaged material to 1 for completely damaged material. The stress-strain relation reads:

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

$E_0$ and $v_0$ are the Young's modulus and the Poisson's ratio of the undamaged material; $\varepsilon_{ij}$ and $\sigma_{ij}$ are the strain and stress components, and $\delta_{ij}$ is the Kronecker symbol. 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_{ijkl}^{0}$ is the stiffness of the undamaged material. This energy is assumed to be the state potential. The damage energy release rate is

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

with the rate of dissipated energy:

$$ \dot{\phi} = -\frac{\partial\rho\psi}{\partial d} $$
Since the dissipation of energy ought to be positive or zero, the damage rate is constrained to the same inequality because the damage energy release rate is always positive.

## 1.3 EVOLUTION OF DAMAGE

The evolution of damage is based on the amount of extension that the material is experiencing during the mechanical loading. An equivalent strain is defined as

$$ \tilde{\epsilon} = \sqrt{\sum_{i=1}^{3} (\langle \epsilon_i \rangle_{+})^2} \quad (3) $$

where $\langle \cdot \rangle_{+}$ is the Macauley bracket and $\epsilon_i$ are the principal strains. The loading function of damage is

$$ f(\tilde{\epsilon}, \kappa) = \tilde{\epsilon} - \kappa \quad (4) $$

where $\kappa$ is the threshold of damage growth. Initially, its value is $\kappa_0$, which can be related to the peak stress $f_t$ of the material in uniaxial tension:

$$ \kappa_0 = \frac{f_t}{E_0} \quad (5) $$

In the course of loading $\kappa$ assumes the maximum value of the equivalent strain ever reached during the loading history.

If $f(\tilde{\epsilon}, \kappa) = 0$ and $\dot{f}(\tilde{\epsilon}, \kappa) = 0$, then

$$ \begin{cases} d = h(\kappa) \\ \kappa = \tilde{\epsilon} \end{cases} \quad \text{with } \dot{d} \ge 0, \text{ else} \begin{cases} \dot{d} = 0 \\ \dot{\kappa} = 0 \end{cases} \quad (6) $$

The function $h(\kappa)$ is detailed as follows: in order to capture the differences of mechanical responses of the material in tension and in compression, the damage variable is split into two parts:

$$ d = \alpha_t d_t + \alpha_c d_c \quad (7) $$

where $d_t$ and $d_c$ are the damage variables in tension and compression, respectively. They are combined with the weighting coefficients $\alpha_t$ and $\alpha_c$, defined as functions of the principal values of the strains $\epsilon_{ij}^t$ and $\epsilon_{ij}^c$ due to positive and negative stresses:

$$ \epsilon_{ij}^t = (1 - d)C_{ijkl}^{-1}\sigma_{kl}^t, \quad \epsilon_{ij}^c = (1 - d)C_{ijkl}^{-1}\sigma_{kl}^c \quad (8) $$
$$ \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 \quad (9) $$

Note that in these expressions, strains labeled with a single indicia are principal strains. 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$:

$$ \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} \quad (10) $$

## 1.4 IDENTIFICATION OF PARAMETERS

There are eight model parameters. The Young's modulus and Poisson's ratio are measured from a uniaxial compression test. A direct tensile test or three-point bend test can provide the parameters which are related to damage in tension ($\kappa_0$, $A_t$, $B_t$). Note that Eq. 5 provides a first approximation of the initial threshold of damage, and the tensile strength of the material can be deduced from the compressive strength according to standard code formulas. The parameters ($A_c$, $B_c$) are fitted from the response of the material to uniaxial compression. Finally, $\beta$ should be fitted from the response of the material to shear. This type of test is difficult to implement. The usual value is $\beta = 1$, which underestimates the shear strength of the material [7]. Table 1 presents the standard intervals for the model parameters in the case of concrete with a moderate strength.

<table><caption>TABLE 1 STANDARD Model Parameters</caption><tr><td>E<sub>0</sub> ≈ 30,000–40,000 MPa</td></tr><tr><td>v<sub>0</sub> ≈ 0.2</td></tr><tr><td>κ<sub>0</sub> ≈ 1 × 10<sup>-4</sup></td></tr><tr><td>0.7 ≤ A<sub>t</sub> ≤ 1.2</td></tr><tr><td>10<sup>4</sup> ≤ B<sub>t</sub> ≤ 5 × 10<sup>4</sup></td></tr><tr><td>1 ≤ A<sub>c</sub> ≤ 1.5</td></tr><tr><td>10<sup>3</sup> ≤ B<sub>c</sub> ≤ 2 × 10<sup>3</sup></td></tr><tr><td>1.0 ≤ β ≤ 1.05</td></tr></table>
<figure><img src="image_2.png" /><figcaption>FIGURE 1 Uniaxial response of the model.</figcaption></figure>

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

## 2 NONLOCAL DAMAGE

The purpose of this section is to describe the nonlocal enhancement of the previously mentioned damage model. This modification of the model is necessary in order to achieve consistent computations in the presence of strain localization due to the softening response of the material [8].

### 2.1 VALIDITY

As far as the type of loading is concerned, the range of validity of the nonlocal model is exactly the same as the one of the initial, local model. This model, however, enables a proper description of failure that includes damage initiation, damage growth, and its concentration into a completely damaged zone, which is equivalent to a macrocrack.
## 2.2 PRINCIPLE

Whenever strain softening is encountered, it may yield localization of strains and damage. This localization corresponds to the occurrence of bifurcation, and a surface (in three dimension) of discontinuity of the strain rate appears and develops. When such a solution is possible, strains and damage concentrate into a zone of zero volume, and the energy dissipation, which is finite for a finite volume of material, tends to zero. It follows that failure occurs without energy dissipation, which is physically incorrect [1].

Various remedies to this problem can be found (e.g., [5]). The basic idea is to incorporate a length, the so-called internal length, into the constitutive relation to avoid localization in a region of zero volume. The internal length controls the size of the region in which damage may localize. In the nonlocal (integral) damage model, this length is incorporated in a modification of the variable which controls damage growth (i.e., the source of strain softening): a spatial average of the local equivalent strain.

## 2.3 DESCRIPTION OF THE MODEL

The equivalent strain defined in Eq. 3 is replaced by its average $\bar{\varepsilon}$:

$$ \bar{\varepsilon}(x) = \frac{1}{V_r(x)} \int_{\Omega} \psi(x-s) \bar{\varepsilon}(s) ds \quad \text{with} \quad V_r(x) = \int_{\Omega} \psi(x-s) ds \quad (11) $$

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

$$ \psi(x-s) = \exp \left( \frac{4 \|x-s\|^2}{l_c^2} \right) \quad (12) $$

where $l_c$ is the internal length of the nonlocal continuum. The loading function (Eq. 4) becomes $f(\varepsilon^{-}, \chi) = \varepsilon^{-} - \chi$. The rest of the model is similar to the description provided in Section 1.

## 2.4 IDENTIFICATION OF THE INTERNAL LENGTH

The internal length is an additional parameter which is difficult to obtain directly by experiments. In fact, whenever the strains in specimen are homogeneous, the local damage model and the nonlocal damage model are, by definition, strictly equivalent ($\bar{\varepsilon} = \tilde{\varepsilon}$). This can be viewed also as a simplification, since all the model parameters (the internal length excepted)
are not affected by the nonlocal enhancement of the model if they are
obtained from experiments in which strains are homogeneous over the
specimen.

The most robust way of calibrating the internal length is by a semi-inverse
technique which is based on computations of size effect tests. These tests are
carried out on geometrically similar specimens of three different sizes. Since
their failure involves the ratio of the size of the zone in which damage can
localize versus the size of the structure, a size effect is expected because the
former is constant while the later changes in size effect tests. It should be
stressed that such an identification procedure requires many computations,
and, as of today, no automatic optimization technique has been devised for it.
It is still based on a manual trial-and-error technique and requires some
experience. An approximation of the internal length was obtained by Bazant
and Pijaudier-Cabot [2]. Comparisons of the energy dissipated in two tensile
tests, one in which multiple cracking occurs and a second one in which failure
is due to the propagation of a single crack, provided a reasonable
approximation of the internal length that is compared to the maximum
aggregate size $d_a$ of concrete. For standard concrete, the internal length lies
between $3d_a$ and $5d_a$.

## 2.5 HOW TO USE THE MODEL

The local and nonlocal damage models are easily implemented in finite
element codes which uses the initial stiffness or secant stiffness algorithm.
The reason is that the constitutive relations are provided in a total strain
format. Compared to the local damage model, the nonlocal model requires
some additional programming to compute spatial averages. These quantities
are computed according to the same mesh discretization and quadrature as for
solving the equilibrium equations. To speed the computation, a table in
which, for each gauss point, its neighbors and their weight are stored can be
constructed at the time of mesh generation. This table will be used for any
subsequent computation, provided the mesh is not changed. Attention should
also be paid to axes of symmetry: as opposed to structural boundaries where
the averaging region lying outside the structure is chopped, a special
averaging procedure is needed to account for material points that are not
represented in the finite element model.

The implementation of the nonlocal model in an incremental format is
awkward. The local tangent stiffness operator relating incremental strains to
incremental stresses becomes nonsymmetric, and, more importantly, its
bandwidth can be very large because of nonlocal interactions.
# 3 ANISOTROPIC DAMAGE MODEL

## 3.1 VALIDITY

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 on elastoplastic damage which addresses these issues. This anisotropic damage model has been compared to experimental data in tension, compression, compression-shear, and nonradial tension-shear. It provides a reasonable agreement with such experiments [3].

## 3.2 PRINCIPLE

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

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

where $\sigma_{ij}^t$ is the effective stress component, $\epsilon_{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 $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_{nk} n_l) n_i]^2} \quad (14) $$

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

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

find $\sigma_{ij}$ such that $\forall \epsilon_{ij}^*$:

$$ \frac{4\pi}{3} \sigma_{ij} \epsilon_{ij}^* = \int_S \left( \left[ (1-d(\mathbf{n}))n_k \sigma_{kl}^t n_l n_i + (1-d(\mathbf{n})) (\sigma_{ij}^t n_j - n_k \sigma_{kl}^t n_l n_i) \right] \cdot \epsilon_{ij}^* n_j \right) d\Omega \quad (15) $$

Depending on the interpolation of the damage variable $d(\mathbf{n})$, several forms of damage-induced anisotropy can be obtained.

## 3.3 DESCRIPTION OF THE MODEL

The variable $d(\mathbf{n})$ is now defined by three scalars in three mutually orthogonal directions. It is the simplest approximation which 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. In compression or tension-shear problems, plastic strains are also of importance and will be added in the model. When the loading history is not monotonic, damage deactivation occurs because of microcrack closure. The model also incorporates this feature.

### 3.3.1 Evolution of Damage

The evolution of damage is controlled by a loading surface $f$, which is similar to Eq. 4:

$$ f(\mathbf{n}) = n_i \epsilon_{ij}^e n_j - \epsilon_d - \chi(\mathbf{n}) \quad (16) $$

where $\chi$ is a 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:

$$ \begin{array}{l} \text{If } f(n^*) = 0 \text{ and } n_i^* d\epsilon_{ij}^e n_j^* > 0 \\ \text{then} \\ \left\{ \begin{array}{l} dd(\mathbf{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^* \quad (17) \\ \\ d\chi(\mathbf{n}) = n_i^* d\epsilon_{ij}^e n_j^* \\ \\ \text{else } dd(n^*) = 0, \quad d\chi(n) = 0 \end{array} \right. \end{array} $$

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

### 6.13.3.3.2 Coupling with Plasticity

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

$$d\epsilon_{ij} = d\epsilon_{ij}^e + d\epsilon_{ij}^p \quad (18)$$

The evolution of the plastic strain is controlled by a yield function which is expressed in terms of the effective stress in the undamaged material. We have implemented the yield function due to Nadai [6]. 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 (19)$$

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 (20)$$

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 with these surfaces. The hardening rule is given by

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

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.

### 6.13.3.3.3 Crack Closure Effects

Crack closure effects are of importance when the material is subjected to alternated loads. During load cycles, microcracks close progressively and the tangent stiffness of the material should increase while damage is kept constant. A decomposition of the stress tensor into a positive and negative part is introduced: $\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. 14 of the model is modified:

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

where $d_c(\mathbf{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(\mathbf{n})$ is directly deduced from $d(\mathbf{n})$. It is defined by the same interpolation as $d(\mathbf{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] \quad (23) $$

where $\alpha$ is a model parameter.

### 6.13.3.4 IDENTIFICATION OF PARAMETERS

The constitutive relations contain six parameters in addition to the Young's modulus of the material and the Poisson's ratio. The first series of three parameters ($\epsilon_d$, $a$, $\alpha$) deals with the evolution of damage. Their determination benefits from the fact that, in tension, plasticity is negligible, and hence $\epsilon_d$ is directly deduced from the fit of a uniaxial tension test. If we assume that in uniaxial tension damage starts once the peak stress is reached, $\epsilon_d$ is the uniaxial tensile strain at the peak stress (Eq. 5). Parameter $a$ is more difficult to obtain because the model exhibits strain softening. To circumvent the difficulties involved with softening in the computations without introducing any nonlocality (as in Section 2), the energy dissipation due to damage in uniaxial tension is kept constant whatever the finite element size. Therefore, $a$ becomes an element-related parameter, and it is computed from the fracture energy. For a linear displacement interpolation, $a$ is the solution of the following equality where the states of strain and stresses correspond to uniaxial tension:

$$ h\phi = G_f, \quad \text{with} \quad \phi = \int_0^\infty \int_\Omega [\dot{d}(\bar{\mathbf{n}})n_k\sigma_{kl}^t n_l n_i] n_j d\Omega d\epsilon_{ij} \quad (24) $$

where $\phi$ is the energy dissipation per unit volume, $G_f$ is the fracture energy, and $h$ is related to the element size (square root of the element surface in a two-dimensional analysis with a linear interpolation of the displacements). The third model parameter $\alpha$ enters into the influence of damage created in tension on the compressive response of the material. Once the evolution of damage in tension has been fitted, this parameter is determined by plotting the decrease of the uniaxial unloading modulus in a compression test versus
the growth of damage in tension according to the model. In a log-log coordinate system, a linear regression yields the parameter $\alpha$.

The second series of three parameters involved in the plastic part of the constitutive relation is $(q, r, w_0)$. They are obtained from a fit of the uniaxial compression response of concrete once the parameters involved in the damage part of the constitutive relations have been obtained.

Figure 2 shows a typical uniaxial compression-tension response of the model corresponding to concrete with a tensile strength of 3 MPa and a compressive strength of 40 MPa. The set of model parameters is: $E = 35,000$ MPa, $v = 0.15$, $f_t = 2.8$ MPa (which yields $\epsilon_d = 0.76 \times 10^{-4}$); fracture energy: $G_f = 0.07$ N/mm; other model parameters: $\alpha = 12$, $r = 0.5$, $q = 7000$ MPa, $\omega_0 = 26.4$ MPa.

## 3.5 HOW TO USE THE MODEL

The implementation of this constitutive relation in a finite element code follows the classical techniques used for plasticity. An initial stiffness algorithm should be preferred because it is quite difficult to derive a consistent material tangent stiffness from this model. Again, the evolution of

<figure><img src="image_3.png" /><figcaption>FIGURE 2 Uniaxial tension-compression response of the anisotropic model (longitudinal [1], transverse [2], and volumetric [v] strains as functions of the compressive stress).</figcaption></figure>
damage is provided in a total strain format. It is computed after incremental
plastic strains have been obtained. Since the plastic yield function depends on
the effective stress, damage and plasticity can be considered separately (plastic
strains are not affected by damage growth). The difficulty is the numerical
integration involved in Eq. 15, which is carried out according to Simpson's
rule or to some more sophisticated scheme.

REFERENCES

1. Bazant, Z.P. (1985). Mechanics of distributed cracking. *Applied Mech. Review* 39: 675-705.

2. Bazant, Z.P., and Pijaudier-Cabot, G. (1989). Measurement of the characteristic length of nonlocal continuum. *J. Engrg. Mech. ASCE* 115: 755-767.

3. 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* 4: 339-359.

4. 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, France.

5. Muhlhaus, H. B., ed. (1995). *Continuum Models for Material with Microstructure*, John Wiley.

6. Nadai, A. (1950). *Theory of Flow and Fracture of Solids*, p. 572, vol. 1, 2nd ed., New York: McGraw-Hill.

7. Pijaudier-Cabot, G., Mazars, J., and Pulikowski, J. (1991). Steel-concrete bond analysis with nonlocal continuous damage. *J. Structural Engrg. ASCE* 117: 862-882.

8. Pijaudier-Cabot, G., and Bazant, Z. P. (1987). Nonlocal damage theory. *J. Engrg. Mech. ASCE* 113: 1512-1533.