# Nonlocal Damage Model with Evolving Internal Length

Gilles Pijaudier-Cabot, Khalil Haidar, Jean-François Dubé

► To cite this version:

Gilles Pijaudier-Cabot, Khalil Haidar, Jean-François Dubé. Nonlocal Damage Model with Evolving Internal Length. International Journal for Numerical and Analytical Methods in Geomechanics, 2004, Computational Mechanics of Concrete and Concrete Structures, 28 (7-8), pp.633-652. 10.1002/nag.367. hal-01007119

HAL Id: hal-01007119

https://hal.science/hal-01007119

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
# Non-local damage model with evolving internal length

Gilles Pijaudier-Cabot¹, Khalil Haidar¹ and Jean-François Dubé²

¹R&DO, Institut de Recherches en Génie Civil et Mécanique – GeM/UMR6183, Ecole Centrale de Nantes,
BP 92101, F-44321 Nantes Cedex, France

²LMGC/UMR5508, Université de Montpellier II, Place Eugène Bataillon, F-34095 Montpellier Cedex, France

A modified non-local damage model with evolving internal length, inspired from micromechanics, is developed. It is shown in particular that the non-local influence between two points in the damaged material depends on the value of damage at each of these points. The resulting weight function is non-symmetric and truncated. Finite element results and strain localization analysis on a one-dimensional problem are presented and compared to those of the original non-local damage model. It is shown that in the course of damage localization, the incremental strain profiles expand according to the modified non-local model, instead of shrinking according to the original constitutive relation. Comparisons with experimental data on model materials with controlled porosity are also discussed. Acoustic emission analyses provide results with which the theoretical model is consistent qualitatively. This model also opens the path for durability mechanics analyses, where it has been demonstrated that the internal length in the non-local model should evolve with environmentally induced damage.

KEY WORDS: concrete; damage; cracking; internal length

## 1. INTRODUCTION

Non-local continuum damage, in an integral or gradient format, is a consistent general concept for macroscopic modelling of failure in quasi-brittle materials. The non-local aspect is a requisite for a realistic description of fracture, including crack inception, crack propagation and also structural size effect which is a consequence of the existence of a finite size fracture process zone (FPZ).

Among the governing parameters in the non-local damage model, the internal length plays a pivotal role as it controls the size of the fracture process zone. In most existing proposals, the internal length is a constant parameter. In fact, one could state that it is already so difficult to
measure this parameter with a relatively good accuracy that considering a variation of it may
hardly be envisioned without solid experimental evidences.

There are some theoretical indications, however, which suggest that the internal length should change in the course of the fracture process. With the help of micromechanics, Bazant [1, 2] has proposed a description of the interactions between cracks and voids in the course of failure. The resulting weight function has been derived from statistical analysis by Bazant and Jirasek [3] and it seems quite natural from these studies that it should evolve during failure. Geers has also observed that, in a gradient damage model, the internal length should change so that the continuum-based constitutive relations can describe fracture with a displacement discontinuity across the crack faces when damage is equal to one [4]. If the introduction of an internal length in the constitutive relations results from microcracks and their interactions, there is also a simple reasoning which tells that without microcracks the internal length should vanish. In the course of the degradation process, the internal length should evolve.

Recent experiments on model cementitious materials with controlled porosity, i.e. initial
damage, have attempted to capture the evolution of the size of the fracture process zone with
growing damage and growing initial porosity [5, 6]. As we will see in this paper, acoustic
emission analyses show clearly that the FPZ is widening in the course of its development. An
increase of the width of the FPZ with initial porosity/damage is also observed.

This paper is an attempt to integrate these observations in a modified non-local damage model. In the first part, an analysis based on micromechanics is presented. It is based on several simplifying assumptions and provides qualitative arguments about non-locality of damage in the case where the microcrack density in the material is small. It is shown that the additional strain induced in the material due to the voids depends on the applied remote stress, on the sizes of the voids, and on their distance to the point at which this strain is computed. Experimental results on model materials with initial porosity are recalled in the second part. Acoustic emission analyses provide quantitative results on the development of the FPZ, consistent with the micromechanical analysis and over values of damage that cover a wider range, from diffuse microcracking to macrocrack propagation. An enhanced non-local damage model is developed in the third part. The non-local influence of one point on another depends on the state of strain at the first one; it does not exist without damage. The consequences of this enhancement on the inception of strain localization are considered from theoretical and numerical points of views in the last section. This non-local damage model should be also capable of capturing the evolution of the internal length due to the variation of the initial porosity in a heterogeneous material.

## 2. SIMPLIFIED MICROMECHANICAL MODEL

The objective of the foregoing simplified analysis is to obtain qualitative expressions of non-local effects induced by the presence of voids in an elastic, homogeneous material. These expressions deal with the local strains in the material—in between the voids essentially. In the non-local damage model developed at the end of this paper, we shall consider that it is these local strains that control damage growth and the non-local effects will be inspired from the present derivation.

Consider an infinite isotropic two-dimensional solid subjected to a remote uniform stress field $\sigma_\infty$. Our purpose is qualitative and we shall consider for more simplicity $\sigma_\infty$ as isotropic. Under this simplification, the considered problem is essentially one dimensional, the strain and stress
being volumetric. This solid contains microcracks caused by the loading history that are distributed in the material arbitrarily. Each crack induces a modification of the local stress and strain fields. In order to evaluate this perturbation, we assume that the microcracks are far enough from each other. Under this assumption and as we will see in the next section, the induced remote stress and strain fields due to a crack are similar to those generated by a circular void. Therefore, cracks are going to be replaced by circular voids in the foregoing derivations. The calculation of the perturbation stresses and strains between the voids uses the superposition theorem. It follows the method proposed by Kachanov [7], although we will deal here with a much simpler, not to say simplistic, problem. Note that cracks or voids are far from each others and that the perturbation stresses are computed far from them. We look at a material that has been slightly damaged only. The following derivation cannot be applied to cases where the material contains a large number of microcracks, or close to material failure when damage exhibits a localized pattern.

## 2.1. Interaction between two voids in an elastic material

Consider the medium described in Figure 1. The diameters of the two circular inclusions, denoted as $S_1$ and $S_2$, are $a_1$ and $a_2$, respectively. We may decompose this problem into two

<figure><img src="image_2.png" /><figcaption>Figure 1. Superposition scheme for an isotropic solid with a two circular voids.</figcaption></figure>
sub-problems:

* *Sub-problem I:* The solid is considered without any inclusion. It is subjected to the remote tensile volumetric stress $\sigma_{\infty}$.

* *Sub-problem II:* Inside each inclusion, a distribution of normal pressures $P_1 = P_2 = -\sigma_{\infty}$ is applied.

Superposition of these two sub-problems provides the distribution of stress inside the medium containing voids.

Consider now sub-problem II. Again, we can apply the principle of superposition (Figure 2) in order to compute the interaction stress fields due to the presence of each inclusion, and subsequently the distribution of internal normal pressure inside the voids: in sub-problem II-1 (Figure 2), the inclusion $S_1$ is considered alone. It is loaded by an unknown normal pressure $\bar{P}_1$, which causes on the imaginary location of $S_2$ a normal pressure $p_{21}$. We have a similar situation for sub-problem II-2. In order to evaluate $\bar{P}_1$ and $\bar{P}_2$, we may again use the superposition principle and write

$$ 
\begin{aligned}
P_1 &= \bar{P}_1 + p_{12} \\
P_2 &= p_{21} + \bar{P}_2
\end{aligned} 
\qquad (1) $$

Note that in an exact derivation, $\bar{P}_1, p_{21}, \bar{P}_2, p_{21}$ should not be normal pressures only. Tangential distributed forces should also be present. These terms are neglected in our equations, but a complete derivation, that includes tangential terms, would not change the conclusions at which we will arrive as far as the mathematical form of the leading terms entering in the expression of the pressures are concerned. It would be, however, slightly more complex.

The interaction pressures have the following expressions:

$$ p_{12} \propto -\bar{P}_2 \times \frac{a_2^2}{r^2} \qquad (2a) $$

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

<figure><img src="image_4.png" /><figcaption>Figure 2. Superposition scheme in subproblem II.</figcaption></figure>
and

$$p_{21} \propto -\overline{P_1} \times \frac{a_1^2}{r^2} \qquad (2b)$$

where r is the distance between the centres of the two voids. It is the perturbation stress field due to a circular inclusion loaded by an internal pressure obtained from the well-known Eshelby solution (see e.g. Reference [8]). It is also similar to the long-range asymptotic crack influence function derived by Bazant [2]. According to our assumptions, the first Equation (2a) holds far from void 2 and the second holds far from void 1. These expressions are singular r goes to zero, but they are not correct in such a case since we deal with a diffuse void distribution only.

We can recast Equation (1) into an algebraic system format:

$$\begin{bmatrix} P_1 \\ P_2 \end{bmatrix} = \begin{bmatrix} 1 & -l_{12}^2 \\ -l_{21}^2 & 1 \end{bmatrix} \begin{bmatrix} \overline{P_1} \\ \overline{P_2} \end{bmatrix} \qquad (3)$$

with the notations $l_{12} \propto a_2/r$ and $l_{21} \propto a_1/r$. These interaction terms, also called influence coefficients, are assumed to be constant over each void location. Since $P_1 = P_2 = -\sigma_\infty$ the solution of this system of two equations is

$$\overline{P}_1 = -\sigma_\infty \left( \frac{1 + l_{12}^2}{1 - l_{21}^2 l_{12}^2} \right) \quad \text{and} \quad \overline{P}_2 = -\sigma_\infty \left( \frac{1 + l_{21}^2}{1 - l_{12}^2 l_{21}^2} \right) \qquad (4)$$

Since the voids (or cracks) are far away from each other, $r \gg a_1, a_2$ and they are assumed to be constant over the void inner surfaces. Equation (4) may be further simplified:

$$\overline{P}_1 \approx -\sigma_\infty(1 + l_{12}^2) \quad \text{and} \quad \overline{P}_2 \approx -\sigma_\infty(1 + l_{21}^2) \qquad (5)$$

If the tangential forces had been considered, they would be proportional to $l_{12}^2$ and $l_{21}^2$ on each void surfaces, respectively, which indicates that the normal and tangential terms that reflect the influence of interaction are of a similar format, with different pre-factors (angular functions) in the influence coefficients.

## 2.2. Stress and strain fields in an elastic material with voids

The extension to the problem of n circular interacting voids $S_i$ ($i = 1 - n$) of diameter $a_i$, arbitrarily distributed in an elastic matrix follows exactly the same technique, with the same simplifying assumptions (Figure 3). The remote traction $\sigma_\infty$ is transformed into distributed normal forces $P_i = -(\sigma_\infty \cdot n_i)$ acting inside each void $S_i$, where $n_i$ is the outward normal vector to the inner contour $\Gamma_i$ of each void $S_i$. The effect of the pressure inside the voids $S_j$ on void $S_i$ is written as the sum of the $p_{ij}$ ($i \neq j$), where $p_{ij}$ represents the normal pressure at the imaginary location of $S_i$ produced by $S_j$:

$$p_{ij} \propto -\overline{P}_j \times \frac{a_j^2}{r_{ij}^2} \qquad (6)$$

$r_{ij}$ being the distance between the centres of voids i and j. Again, we replace the problem of n interacting voids by the superposition of n problems in which each void is considered alone in
<figure><img src="image_5.png" /><figcaption>Figure 3. Superposition scheme for an isotropic solid with n circular voids.</figcaption></figure>

the infinite solid and subjected to an unknown pressure $\bar{P}_j$. Superposition yields

$$P_i = \bar{P}_i + \sum_{j=1(j \neq i)}^{n} p_{ij} = \bar{P}_i - \sum_{j=1(j \neq i)}^{n} (l_{ij}^2 \times \bar{P}_j) \quad (7)$$

where $l_{ij} \propto a_j/r_{ij}$ are the influence coefficients representing the average normal pressure at the boundary of the imaginary void $S_i$ caused by a unit uniform normal pressure inside void $S_j$ when all other voids do not exist. Since the voids are sufficiently far from each others, Equation (5) may easily be generalized to

$$\bar{P}_i \approx P_i \left( 1 + \sum_{j=1(j \neq i)}^{n} l_{ij}^2 \right) \quad (8)$$

On each contour $\Gamma_i$ of each voids we define

$$\sigma_i = -\bar{P}_i \quad (9)$$

where $\sigma_i$ is the volumetric stress that equilibrates the normal pressure applied inside void $i$, to which the effect of interaction is added in order to equilibrate the remote stress $\sigma_\infty$ on the void surface. For simplicity, we define also $l_{ij} = 0$. With these notations Equation (8) becomes

$$\sigma_i = \sigma_\infty \left( 1 + \sum_{j=1}^{n} l_{ij}^2 \right) \quad (10)$$
Let us now evaluate the stress acting in the material, between the voids. According to the superposition scheme, the resulting local stress field at each point $x_k$ of the isotropic medium is the sum of the contributions due to each internal pressure inside each void plus the remote field:

$$ \sigma(x_k) = \sigma_{\infty}(x_k) + \sum_{i} \sigma_{i}(x_i) l_{ik}^{2} \qquad (11) $$

where $l_{ik} \propto a_i/r_{ik}$ and $r_{ik}$ is the distance between the centre of void $i$ and point $x_k$. Substitution of Equation (10) into the previous equation yields

$$ \begin{aligned} \sigma(x_k) &= \sigma_{\infty}(x_k) + \sum_{i} \sigma_{\infty}(x_i) \left( 1 + \sum_{j=1}^{n} l_{ij}^{2} \right) l_{ik}^{2} \\ &= \sigma_{\infty}(x_k) + \sum_{i} \sigma_{\infty}(x_i) l_{ik}^{2} + \sum_{i} \sum_{j=1}^{n} l_{ij}^{2} l_{ik}^{2} \sigma_{\infty}(x_i) \end{aligned} \qquad (12) $$

where we have also considered that the remote stress is not necessarily constant over the body.

The stress acting in the material in between the voids is increased compared to the remote stress $\sigma_{\infty}$. This is due to the presence of voids which reduce the net resisting cross-section of the material and it is expected that the local stresses are higher than the remote stress which is also the average stress applied to the material containing voids. In the right-hand side of this equation, the second term is the perturbation induced by the voids in the material. The third term results from the interaction between the voids. It is a product of two small quantities $l_{ik}^2 l_{ij}^2$. Therefore, we shall neglect this term compared to the others, which means that interactions between voids are assumed to induce a second-order effect compared to the perturbation stress due to the voids. Interactions ought to be included if the voids become close to each other's or if one wants to evaluate the perturbation stress nearby a void.

The stress acting at point $x_k$ of the material due to the remote stress reduces to

$$ \sigma(x_k) \approx \sigma_{\infty}(x_k) + \sum_{i} l_{ik}^2 \sigma_{\infty}(x_i) \qquad (13) $$

This equation is typically non-local, the voids produce a stress amplification. At each material point, this amplification depends on the distance to the void considered, the remote stress, and the square of the radius of the void $a_i$. One should not pay attention to the exponent 2 in Equation (13) as it is related to the fact that the problem is two dimensional. In 3D, this exponent is 3 in all the above equations (at the numerator and at the denominators). Equation (13) is very similar to the one derived by Bazant [2] with a more simple influence function here.

We are going now to derive the local strain field with a view to the expression that controls the growth of damage locally. Since the material in between the voids is elastic, we have

$$ \varepsilon(x_k) = \frac{\sigma(x_k)}{K} = \frac{\sigma_{\infty}(x_k)}{K} + \sum_{i} l_{ik}^2 \frac{\sigma_{\infty}(x_i)}{K} \qquad (14) $$

where $K$ is the bulk stiffness of the isotropic elastic medium. We can see here that the strain in the material is modified due to the voids. Same as for the stress, there is a non-local amplification. Our goal is to recover an expression of the strain that is similar to the quantity that controls damage in continuum relations. For this, we need first to define such a constitutive relation. We will use the isotropic damage model:

$$ \sigma_{\infty}(x_k) = K(1 - d(x_k))\varepsilon_{\infty}(x_k) \qquad (15) $$
Substitution in Equation (14) yields

$$ \varepsilon(x_k) = (1 - d(x_k))\varepsilon_\infty(x_k) + \sum_i l_{ik}^2 (1 - d(x_i))\varepsilon_\infty(x_i) \quad (16) $$

Equations (14), (16) provide the actual, local, strain in the material in between the voids. If void nucleation is strain controlled, it is this strain which should control the growth of material damage in a continuum sense.

It may be interesting to modify Equation (16) in order to exhibit the correction to the remote strain $\varepsilon_\infty$ due to the presence of the voids. For this, we divide Equation (16) by $(1 - d)$, which is possible because damage is small. We obtain

$$ \bar{\varepsilon}(x_k) = \frac{\varepsilon(x_k)}{(1 - d(x_k))} = \varepsilon_\infty(x_k) + \sum_i l_{ik}^2 \frac{(1 - d(x_i))}{(1 - d(x_k))} \varepsilon_\infty(x_i) \quad (17) $$

$\bar{\varepsilon}(x_k)$ is derived from the local strain. It is the sum of two quantities: the remote strain and the corrections due to the presence of surrounding voids. Same as for Equations (14), (16) it is non-local, but it can also be seen as a weighted average of the remote strain distribution, a quantity which is very similar to those used in integral non-local models in order to control damage growth (see e.g. Reference [9]). It is this similarity which will inspire us the modified non-local model discussed in Section 4. Obviously, the definition of $\bar{\varepsilon}(x_k)$ is rather arbitrary and any expression derived from Equation (14) may serve the purpose of deriving a continuum model.

The non-local component in Equation (17), due to the presence of the voids, is a function of the remote strain field and damage. In particular, the terms $l_{ik}^2$ are functions of the square of the void radius $a_i^2$ which can be related to damage, according to many homogenization schemes (see e.g. Reference [10]). These theories deal, however, with uniformly distributed damage and in this case, the non-local value of damage equals the local one. Considering Equation (6), we shall replace the interaction terms $l_{ik}^2$ by a function of the local strain, denoted $f(\varepsilon_\infty(x_i))$ generically, divided by the square of the distance from the void to the considered point of co-ordinate $x_k$. Since damage is small with small gradients too $(1 - d(x_i))/(1 - d(x_k)) \approx 1$, and Equation (17) is further simplified:

$$ \bar{\varepsilon}(x_k) = \varepsilon_\infty(x_k) + \sum_i \frac{f(\varepsilon_\infty(x_i))}{r_{ik}^2} \varepsilon_\infty(x_i) \quad (18) $$

Non-local effects appears only when voids are present. They should grow when the number of voids increases (influence of function $f(\varepsilon_\infty(x_i))$) and as the spacing between the void decreases (when $r_{ik}$ decreases). We shall try to derive a constitutive model which has the same characteristics in the fourth section. Before that, we are going to discuss experimental results which exhibit the variation of the width of the process zone as damage develops. Since the width of the process zone is related to non-local effects, and more precisely to the internal length in non-local damage models, we are going to obtain additional data on how non-locality evolve in the course of damage growth, without considering the restrictive assumptions devised in this theoretical analysis.
### 3. EXPERIMENTS WITH ACOUSTIC EMISSION ANALYSIS

Acoustic emission (AE) is an experimental tool well suited for monitoring microcracking processes. Several works have focused on relating acoustic emission characteristics to the properties of the fracture process zone [11] and using AE source location analysis to evaluate damage localization [12, 13]. In the present study, attention is focused on the evolution of the FPZ in the course of damage propagation.

The experimental programme which is briefly recalled here (for details, see Reference [5] or Reference [6]), aimed at determining the mechanical characteristics (compressive strength, elastic modulus and fracture energy) of mortar specimens with a controlled microstructure. This control was achieved by adding inclusions of weak mechanical characteristics (polystyrene beads) in a matrix having a reduced porosity. Specimens were made with a mix which consisted of ordinary Portland cement CPA-CEMI 52.5, polystyrene beads, normal—density fine sand with a maximum size of 2 mm, a superplasticizing agent and water, with a cement/sand ratio of 0.46 and a water/cement ratio of 0.4. Expanded polystyrene spheres of 3–7 mm mean diameter were used as aggregate in the mix design. Expanded polystyrene beads consist essentially of air 98% and are considered as voids.

Four different mixes of densities 2.0, 1.8, 1.6 and 1.4, having polystyrene content $g$ of 13–22, 31 and 39%, respectively, were tested in addition to the reference material of density 2.3 (mortar without inclusions). Three point bend tests were performed on notched concrete specimens of depths were $D = 160$ mm and thickness $b = 40$ mm. The length to depth ratio was $L/D = 8:3$ and the span to depth ratio was $I/D = 2.5$. A notch of depth $D/6$ and thickness 1.5 mm was placed in each specimen using a steel plate in the mould before casting.

Transducers were placed around the expected location of the process zone, in a triangular array (Figure 4). The AE localization measurements were performed for three different material densities (2.3, 1.8) in order to compare the effect of the material density on the crack band width. The cumulative locations of acoustic events throughout the test is shown in Figure 5 for the material with density 2.3. The accuracy of the localization of AE events was ±5 mm. In order to process these results and to visualize the evolution of the FPZ, the specimen was divided into a

<figure><img src="image_6.png" /><figcaption>Figure 4. Specimen geometry and AE transducer locations.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5. Cumulative location of AE events: for the experiment on material density 2.3.</figcaption></figure>

square grid (1 cm × 1 cm) and the AE events located within each element of the grid were counted. This technique allowed the determination of the shape of fracture process zone in its final stage, but incremental counts (in between two loading stages) provided also the evolution of the FPZ as the failure process developed. Figure 6 shows the profiles of the number of events as a function of the horizontal co-ordinate, perpendicular to the direction of the crack propagation, for material densities 2.3 and 1.8 in the course of loading. The macrocrack is located at $x = 67$ mm. It is clear that the width of the FPZ enlarges in the course of damage progression. Furthermore, the width is larger for the porous material compared to the reference one (without inclusions). Figure 7 shows the profiles where the energy of the acoustic events has been considered instead of their number (same data interpretation technique). We may observe that the trends are very similar to those in Figure 6.

The evolution of the width of the FPZ is shown in Figure 8. On the horizontal axis is the recorded load applied to the beam. It increases and then decreases in the post peak regime. The width of the FPZ is defined as the intersection between the profiles given in Figure 6 and a horizontal line placed at 20% of the maximum number of event. The value of the threshold of 20% is chosen rather arbitrarily; different values would provide different results for the width of the FPZ but the evolution of the width in the course of crack propagation would be the same still [5]. We can observe several different characteristics: first, the width is increasing; second the slope of the two curves is approximately the same in this co-ordinate system; third, there are different values of the width of the FPZ at the origin depending on the initial porosity. It is worth noting that acoustic events are still recorded for very low values of the load in the post
<figure><img src="image_8.png" /><figcaption>Figure 6. Evolution of the cumulative events in a cross-section perpendicular to the crack propagation at different loading state: (a) material density 2.3; and (b) material density 1.8. The maximum number of events increases when the crack forms.</figcaption></figure>

peak regime. A possible reason for this is that even when the macrocrack is formed, there are still bridges of material in between the crack faces. Since some load is transmitted, microcracking may still occur. Another possibility is that some events may be due to friction across the crack faces. These events have not been filtered off from the records of acoustic emission. One would expect, however, that events due to friction be concentrated nearby the macrocrack and would not contribute to the acoustic activity away from the macrocrack.

Finally, let us remark that the scale on the horizontal axis is not a linear function of time (or displacement since these are displacement-controlled experiments). Very low values of the load in the post peak regime are reached for very large values of time. In a plot of the width of the FPZ versus time, the experimental data would not follow a linear curve but rather a curved one reaching a maximum width asymptotically.

The model derived in Section 2 can be adapted in order to account for the initial porosity of the material. Equation (13) becomes

$$ \sigma(x_k) \approx \sigma_\infty(x_k) + \sum_i l_{ik}^2 \sigma_\infty(x_i) + \sum_j l_{jk}^2 \sigma_\infty(x_j) \quad (19) $$
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7. Evolution of the cumulative energy due to acoustic events in a cross-section perpendicular to the crack propagation at different loading state: (a) material density 2.3; and (b) material density 1.8. The maximum energy increases when the crack forms.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8. Evolution of the width of the fracture process zone vs the applied load for the two material densities. The load on the horizontal axis is increasing and then decreasing in the post peak regime.</figcaption></figure>
where $l_{jk}^2 \sigma_\infty(x_j)$ is the perturbation stress due to each initial void $j$ of diameter $a_j$. The derivation of the continuum model follows the same steps as in Section 2, under the same assumptions and considering that the initial voids in the model material are sufficiently far from each others. The continuum constitutive relation which accounts for the initial porosity is chosen as $\sigma = (1-d)(1-V)K_\varepsilon$ where $V$ is damage due to the initial porosity. This model is similar to the one devised by Gérard et al. [14] among others, for chemo-mechanical damage. The strain governing damage growth becomes

$$ \tilde{\varepsilon}(x_k) = \varepsilon_\infty(x_k) + \sum_i \frac{f(\varepsilon_\infty(x_i))}{r_{ik}^2} \varepsilon_\infty(x_i) + \sum_j \frac{V(x_j)}{r_{jk}^2} \varepsilon_\infty(x_j) \quad (20) $$

## 4. MODIFIED NON-LOCAL DAMAGE MODEL

In the original non-local damage model, the width of the FPZ is related to the internal length. According to bifurcation analyses [15] it is in fact proportional to the internal length. Therefore, the experimental data presented in the previous section show that the internal length should increase in the course of damage. Furthermore, micromechanics shows that non-local effects increase with damage. Again that can be captured considering that the internal length should increase in a non-local model, increasing thus the relative weight of neighbouring points in the computation of the non-local averages. In this section, we are going to enhance the original model accordingly. Before that, let us first recall the non-local damage model from which we will derive the modified one. The stress strain relation is

$$ \sigma_{ij} = (1 - d)\Lambda_{ijkl}\varepsilon_{kl} \quad (21) $$

where $\sigma_{ij}$ and $\varepsilon_{ij}$ are the components of the stress and strain tensors, respectively, $(i,j,k,l \in [1, 3])$, $\Lambda_{ijkl}$ are the initial stiffness moduli, and $d$ is the scalar damage variable. The material is isotropic, with $E$ and $\nu$ the initial Young's modulus and Poisson's ratio, respectively. For the purpose of defining damage growth, the equivalent strain is introduced:

$$ \tilde{\varepsilon} = \sqrt{\sum_{i=1}^{3} \langle \varepsilon_i \rangle_{+}^2}. \quad (22) $$

where $\langle \cdot \rangle_+$ is the Macauley bracket and $\varepsilon_i$ are the principal strains. In the non-local damage model, the variable $\tilde{\varepsilon}$ which enters in the equations governing damage growth is

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

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

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

$l_c$ is the internal length of the non-local continuum. The rest of the constitutive model, and in particular the equations governing damage growth, can be found for instance in Reference [9].

We look now for an enhancement of this model which accounts for the main characteristics of Equation (18) of the micromechanical model. According to Equation (18) and at a given point
x, the non-local influence of a void located in its neighbourhood at co-ordinate z depends on the value of the scalar function $f(\varepsilon(z))$. It vanishes if there is no void, when damage is zero and $f(\varepsilon(z)) = 0$, and it increases afterwards. In a continuum model, the non-local contribution to the equivalent strain $\bar{\varepsilon}(x)$ of the local equivalent strain $\tilde{\varepsilon}(z)$ at point z in Eq. (23) should be weighted following the same principle. It is zero if damage is equal to zero at point z and it increases with damage. Whether the interactions are neglected of not in the micromechanical model Eq. (12), one will arrive to the same conclusions, with a slightly more complex form of the non-local influence functions.

There are several possibilities for capturing such an effect in the above non-local model. The first one can be to replace Eq. (21) with

$$ \bar{\varepsilon}(x) = \frac{1}{V_r(x)} \int_{\Omega} \psi(x-s)f(\tilde{\varepsilon}(s)).\tilde{\varepsilon}(s) \, ds \quad \text{with} \quad V_r(x) = \int_{\Omega} \psi(x-s)f(\tilde{\varepsilon}(s)) \, ds \quad (25) $$

The weight of point s, $\psi(x-s)f(\tilde{\varepsilon}(s))$, is increasing as $f(\tilde{\varepsilon}(s))$ increases, when damage increases at this point. A second possibility is to modify the value of the internal length in the weight function directly in order to capture the increase of the weight:

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

In this expression which will be further used in this paper, the internal length is increasing. A linear function

$$ l_c(\tilde{\varepsilon}) = \alpha f(\tilde{\varepsilon}) + l_{c0} \quad (27) $$

is a simple mathematical expression that can be considered. We noted in Section 2 that the function $f(\tilde{\varepsilon})$ depends on the value of local damage. This is the reason why it is expressed in term of the local equivalent strain. In the absence of specific information on this function, it shall be considered in the applications to have exactly the same format as the function which controls the growth of damage in the constitutive relation, with the same model parameters, although it has a local argument instead of a non-local one. Hence, it will vary between 0 and 1. $l_{c0}$ is related to initial defects in the material. It will be equal to zero in the following example. $\alpha$ is a material constant; in the example shown in this paper it is equal to the maximum value of the internal length $l_c$.

This modified non-local damage model can be looked at as the integral equivalent of the model devised by Geers et al. [4] where the internal length increases in the course of damage. It can be implemented similarly. Our formulation has, however, an important feature that is different from Geers' proposal: a material point where damage is equal to zero does not influence its neighbour where damage is non-zero. Conversely a material point at which damage is non-zero has an influence on a point where damage is zero.

Another possible enhancement of the model deals with the exact form of the weight function. Equation (18) suggests that it should include a factor $r^2$ or $\|x-z\|^2$ (in 2D and with an exponent $-3$ in 3D) but the micromechanical model in Section 2 is based on assumptions that are too restrictive for that. Cracks and voids are assumed to be far from each others and the perturbation stress is evaluated far from them. A more accurate form of the weight function has been derived following a more comprehensive analysis of the interaction between microcracks or voids [2, 3]. $\|x-z\|$ enters in the weight function to a power which is not $(-2)$ only. There are higher-order exponents which come from crack interaction too (see also Equation (12)). The
difficulty that is inherent to such an analysis is that the interaction depends on the topology of the defects in the elastic material (distribution, orientation with respect to the loads, etc.). The expression of the weight function becomes much more complex, and stress dependent too. This is the reason why the model still relies on a Gauss distribution of the weights (modified by damage). Other distributions could be used as well.

An enhanced non-local model can again be devised for materials with an initial distribution of voids such as those used in the experiments in Section 3. One possibility is a subsequent modification of the weight function

$$ \psi(x - s, \tilde{e}(s), V(s)) = \exp\left(-\frac{4\|x - s\|^2}{l_c^2(\tilde{e}(s), V(s))}\right) \quad \text{with } l_c(\tilde{e}, V) = \alpha f(\tilde{e}) + \beta V + l_{c0} \quad (28) $$

where $\beta$ is an additional constant which accounts for the variation of internal length due to the initial material porosity. Of course such a proposal has to be calibrated and compared with experimental data.

## 5. ONE-DIMENSIONAL STRAIN LOCALIZATION ANALYSIS

The aim of this section is to investigate failure due to strain localization with the proposed formulation and to compare with results from the original non-local damage model. For that, a one-dimensional setting is considered and we look first at computational results.

The example is a bar loaded in tension with an initial imperfection (a decrease of 20% of the Young's modulus) as presented in Figure 9. Owing to this defect, localization is triggered and it is possible to visualize the development of the damage zone. The geometry of the bar and the material parameters in the damage model are provided in Table I. The bar length is 100 mm. The finite element mesh is made of 32 constant size elements (smaller that the internal length). A variation of the weight function according to Equations (26), (27) has been used (see the parameters $\alpha$ and $l_{c0}$ in Table I). For the original non-local model, the internal length is $l_c = 20$ mm. It is also the maximum value in the modified model with evolving internal length.

<figure><img src="image_12.png" /><figcaption>Figure 9. One-dimensional bar loaded in tension.</figcaption></figure>
<table><caption>Table I. Characteristics of the bar and the material parameters.</caption><thead><tr><th>Test Configuration and model parameters</th><th>Quantity</th><th>Value</th></tr></thead><tbody><tr><td rowspan="2">Geometry</td><td>Length of the bar</td><td>100 mm</td></tr><tr><td>Cross-section</td><td>1 mm<sup>2</sup></td></tr><tr><td rowspan="2">Elastic parameters</td><td>Young's modulus</td><td>33 000</td></tr><tr><td>Young's modulus of the defect element</td><td>24 750</td></tr><tr><td rowspan="5">Non-local parameters</td><td>α</td><td>20 mm</td></tr><tr><td>ε<sub>D0</sub></td><td>1.5 × 10<sup>-4</sup></td></tr><tr><td>A<sub>t</sub></td><td>1</td></tr><tr><td>B<sub>t</sub></td><td>6000</td></tr><tr><td>l<sub>c0</sub></td><td>0</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Figure 10. Damage evolution in the one-dimensional tension test in the original (a) and modified (b) non-local damage model.</figcaption></figure>

We compare in Figure 10 the evolutions of damage in the bar obtained with the original non-local damage model with fixed internal length (a) and with the present one with an evolving internal length (b). This figure shows that the modified non-local model yields a damage profile with a sharp boundary. The kinetics of damage growth is different. The width of the damage zone expands according to the modified model, instead of having a relatively constant width
according to the original constitutive relation with a fixed internal length. This is in qualitative agreement with AE data.

These evolutions of the damage localization profiles can be further explained with the help of a one-dimensional bifurcation analysis on an infinite bar in an initial state of homogeneous tensile strain and damage ($\epsilon_0$, $d_0$). The analysis follows the steps described, for instance, in Reference [15]. The incremental constitutive relations reads

$$d\sigma = (1 - d_0)E \, d\epsilon - E\epsilon_0 \frac{\partial g(\bar{\epsilon})}{\partial \bar{\epsilon}} \, d\bar{\epsilon} \quad (29)$$

where $\partial g(\bar{\epsilon})/\partial \bar{\epsilon}$ stands for the derivative of the evolution law of damage (provided in an integrated form usually), assuming that incremental damage is strictly positive everywhere. The increment of the non-local strain can be computed without difficulties:

$$d\bar{\epsilon}(x) = \frac{\int_{\Omega} \psi(x-s, \epsilon_0) \, d\epsilon(s) \, ds}{\int_{\Omega} \psi(x-s, \epsilon_0) \, ds} \quad (30)$$

The weight function used here is given in Equation (26) where the initial strain has been substituted to the equivalent strain, which is equal to the uniaxial tensile strain here.

Equilibrium requires that the incremental stress be homogeneous $d\sigma_{,x} = 0$. Substitution of Equations (29), (30) yields the following governing equation:

$$(1 - d_0)E \, d\epsilon_{,x} - E\epsilon_0 \frac{\partial g(\bar{\epsilon})}{\partial \bar{\epsilon}} \frac{\int_{\Omega} \psi(u, \epsilon_0) \, d\epsilon(x+u)_{,x} \, du}{\int_{\Omega} \psi(u, \epsilon_0) \, du} = 0 \quad (31)$$

in which $u = x - s$. Compared to the equation obtained with the original damage model, this expression is quite similar, except that the initial state of strain enters in the expression of the weight function. The solution to this equation is obtained by considering that admissible incremental displacements are harmonic: $du(x) = A \exp(i\xi x)$ where $\xi$ is the angular frequency and $A$ is the amplitude of the bifurcation mode. Substitution of these possible solutions in the equation of equilibrium yields

$$A \exp(i\xi x) \xi^2 \left[ (1 - d_0)E - E\epsilon_0 \frac{\partial g(\bar{\epsilon})}{\partial \bar{\epsilon}} \bar{\Psi}^* \right] = 0 \quad (32)$$

$$\bar{\Psi}^* = \frac{\int_{\Omega} \psi(u, \epsilon_0) \exp(i\xi u) \, du}{\int_{\Omega} \psi(u, \epsilon_0) \, du} \quad (33)$$

$\bar{\Psi}^*$ is the Fourier transform of the modified weight function. The admissible angular frequency of the localized modes is

$$\xi = \sqrt{-\frac{2}{\alpha f(\epsilon_0) + l_{c0}}} \ln\left(\frac{(1 - d_0)E}{(E\epsilon_0) \frac{\partial g(\bar{\epsilon})}{\partial \bar{\epsilon}}}\right) \quad (34)$$
<figure><img src="image_14.png" /><figcaption>Figure 11. Variation of the admissible wavelength (mm) of the localized modes according to the original and enhanced damage models.</figcaption></figure>

We have substituted in the above Equation (27) which provides here the variation of the internal length with the initial state of homogeneous strain. Figure 11 shows the variation of the wavelength $2\pi/\xi$ of the localized incremental strain as a function of the initial state of strain $\epsilon_0$. For this plot, we have used the model parameters in Table I. The comparison with the original damage model in which the internal length is constant $l_c = \alpha + l_{c0}$ is also provided. We can see that the evolution of the wavelengths of the localized mode are different according to the original and modified damage models. With the modified one, the wavelength increases in the course of damage and with the original one, it shrinks in the course of damage growth. These results are quite consistent with the numerical ones in Figure 10. For relatively high values of the initial strain (initial damage), the two models are quite equivalent and Figure 11 shows that the wavelength tends to decrease slightly instead of increasing. This is not representative because it corresponds to initial states located far beyond the peak strain in the softening regime. It is quite seldom that localization is so much delayed that it is initiated when damage is quite high, far in the post peak regime of the stress strain relation.

## 6. CONCLUSIONS

Micromechanics indicates that non-local effects increase in the course of damage. Experiments show also that the width of the FPZ is increasing with damage and with the initial porosity of the material.

In the simplest proposal of a non-local model, a constant value of the internal length is selected. The above features, however, cannot be captured and it is more adequate to implement a non-local damage model in which the internal length is increasing. Numerical analyses and localization analysis show that, in the course of failure, the width of the damaged zone grows, and that the wavelength of localized mode increases too. These results are qualitatively consistent with experimental data deduced from acoustic emission analysis.

In the enhanced non-local model, the influence of point *z* located in the neighbourhood of point *x* at which an non-local average is computed depends on the local value of damage at
point *z* that enters in the expression of the internal length at point *z*. This is a difference with similar models with an evolving internal length where the internal length is a function of damage or strain at the point *x* where the non-local quantity is computed. This difference is motivated by the simplified micromechanical analysis discussed in this contribution.

It should be underlined that the enhanced non-local model with a variation of the internal length remains to be implemented in 2D and 3D analyses. The model implementation requires also the calibration of additional parameters compared to the standard non-local model (see Equation (28)). Already, and under the simplification that the internal length is function of initial damage only, the evolution of the internal length with initial damage has been obtained by Haidar *et al*. [5] from inverse analysis based on 2D non-linear finite element computations. It is a linear function of the porosity, as suggested in Equation (28)). The same technique (see also Reference [16]) could be used for calibrating the two constants in Equation (27) with the help of the experimental data presented in Section 3. The experimental input in the inverse analysis would be the load deflection curves and the evolution of the width of the FPZ with applied displacement. The load deflection curve would be computed from 2D FE analysis and the width of the damage zone could be obtained from damage maps, following the same technique as for the acoustic emission analysis. The functional used in the minimization process leading to the model calibration would have to be modified in order to account for the adjustment of the evolution of the width of the FPZ. Such a calibration procedure is the subject of future works.

Finally, there are some results which point out that when damage is not only due to
mechanical loads, but also to chemical attack such as calcium leaching, there is a variation of
the internal length [17, 18]. One may view calcium leaching as inducing a variation of the
material porosity (as the material is leached) prior to applying any mechanical load on the
sample. The proposed constitutive model provides at least some tracks for explaining and
capturing these observed variations of the internal length.

ACKNOWLEDGEMENTS

This study has been performed with the support of the ACI ‘Jeunes Chercheurs’ from the French ministry
of research. The authors would also like to thank the reviewers for their constructive comments.

REFERENCES

1. Bazant ZP. Why continuum damage is nonlocal: micromechanics arguments. *Journal of Engineering Mechanics (ASCE)* 1990; **117**:1070–1087.

2. Bazant ZP. Nonlocal damage theory based on micromechanics of crack interactions. *Journal of Engineering Mechanics (ASCE)* 1994; **120**:593–617.

3. Bazant ZP, Jirasek M. Damage nonlocality due to microcrack interactions: statistical determination of crack influence function. In *Fracture and Damage in Quasibrittle Structures*, Bazant ZP *et al*. (eds). E&F Spon Pubs: London UK, 1994; 3–17.

4. Geers MGD, de Borst R, Brekelmans WAM, Peerlings RHJ. Strain-based transient gradient damage model for failure analysis. *Computer Methods in Applied Mechanics and Engineering* 1998; **160**:133–154.

5. Haidar K, Pijaudier-Cabot G, Dubé JF, Loukili A. Correlation between the internal length, the fracture process zone and size effect in mortar and model materials. *Materials and Structures*, 2004; in press.

6. Haidar K. Damage modelling of concrete structures—numerical approaches and microstructure effects on the rupture properties. *Ph.D. Thesis*, Ecole Centrale de Nantes, France, 2002 (in French).

7. Kachanov M. Elastic solids with many cracks—a simple method of analysis. *International Journal of Solids Structures* 1987; **23**:23–43.
8. Mura T. *Micromechanics of Defects in Solids* (2nd edn). Martinus Mijhoff Pubs.: The Hague, The Netherlands, 1987.

9. Pijaudier-Cabot G, Mazars J. Damage models for concrete. In *Handbook of Materials Behaviour Models*, Lemaitre J (ed.). vol. 2. Chapter 6—Damage Models. Academic Press: New York, 2001; 500–513.

10. Lemaitre J. *A Course on Damage Mechanics* (2nd edn). Springer: Berlin, 1996.

11. Maji AK, Shah SP. Process zone and acoustic emission in concrete. *Experimental Mechanics* 1988; **28**:27–33.

12. Berthelot JM, Robert JL. Modeling concrete damage by acoustic emission. *Journal of Acoustic Emission* 1987; **7**: 43–60.

13. Otsuka K, Date H, Kurita T. Fracture process zone in concrete tension specimens by X-ray and AE techniques. In *Fracture Mechanics of Concrete Structures*, Mihashi *et al*. (eds). *Proceedings of FRAMCOS-3*, AEDIFICATIO Publishers: Freiburg, Germany, 1998; 3–16.

14. Gérard B, Pijaudier-Cabot G, La Borderie C. Coupled diffusion-damage modelling and the implications on failure due to strain-softening. *International Journal of Solids and Structures* 1998; **35**:4107–4120.

15. Pijaudier-Cabot G, Benallal A. Strain localisation and bifurcation in a Nonlocal continuum. *International Journal of Solids and Structures* 1993; **30**:1761–1775.

16. Le Bellégo C, Dubé JF, Pijaudier-Cabot G, Gérard B. Calibration of non local damage model from size effect tests. *European Journal of Mechanics A/Solids* 2003; **22**:33–46.

17. Le Bellégo C, Pijaudier-Cabot G, Gerard B, Dubé JF, Molez L. Coupled chemical and mechanical damage in calcium leached cementitious structures. *Journal of Engineering Mechanics (ASCE)* 2003; **129**:333–341.

18. Ulm FJ, Heukamp FH, Germaine JT. Residual design strength of cement based materials for nuclear waste storage systems. *Nuclear Engineering and Design* 2002; **211**:51–60.