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

HAL
open science

# Boundary effect on weight function in nonlocal damage model

Abbas Krayani, Gilles Pijaudier-Cabot, Frédéric Dufour

► To cite this version:

Abbas Krayani, Gilles Pijaudier-Cabot, Frédéric Dufour. Boundary effect on weight function in nonlocal damage model. Engineering Fracture Mechanics, 2009, 76 (14), pp.2217-2231. 10.1016/j.engfracmech.2009.07.007. hal-00705454

HAL Id: hal-00705454

https://hal.science/hal-00705454

Submitted on 3 Jul 2017

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

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

Public Domain
# Boundary effect on weight function in nonlocal damage model

Abbas Krayani<sup>a</sup>, Gilles Pijaudier-Cabot<sup>b</sup>, Frédéric Dufour<sup>a</sup>

<sup>a</sup> R&DO, GeM – Institut de Recherche en Génie Civil et Mécanique, Ecole Centrale de Nantes/Université de Nantes/CNRS, 1 rue de la Noë, BP 92101, F-44321 Nantes Cedex 3, France

<sup>b</sup> Laboratoire des Fluides Complexes, UMR CNRS 5150, Université de Pau et des Pays de l'Adour, Allée du parc Montaury, F-64600 Anglet, France

Some insights on boundary effects in nonlocal damage modelling are addressed. Interaction stresses that are at the origin of nonlocality are expected to vanish at the boundary of a solid, in the normal direction to this boundary. Existing models do not account for such an effect. We introduce tentative modifications of the classical nonlocal damage model aimed at accounting for this boundary layer effect in a continuum modelling setting. Computations show that some nonnegligible differences may be observed between the classical and modified formulations. In a one dimensional spalling test, only the modified formulation provides a spall of finite nonzero thickness, whereas spalls smaller than the internal length cannot be obtained according to the original formulation. For the same set of model parameters, including the internal length, the fracture energy derived from the size effect test method is also very different according to both approaches. Parameters in the size effect laws for notched and unnotched specimens, obtained from computation of geometrically similar bending beams, are more consistent with the modified nonlocal model compared to the original nonlocal formulation.

Keywords:

Damage

Nonlocality

Boundary effects

Interactions

Size effect

## 1. Introduction

Most failure models for strain softening materials involve nonlocality. Whether nonlocality is introduced in an integral or in a gradient format, an internal length is added to the material description. Such constitutive relations provide consistent continuum failure models for progressive cracking in quasi-brittle materials (see e.g. [3]) or ductile failure in alloys (see e.g. [16]). In quasi-brittle materials at least, nonlocality finds its origin in the interaction between growing defects in the course of failure. When a microcrack opens, stresses are released and the stress field in the neighbourhood of the crack is modified accordingly which may induce some further cracking elsewhere. These interactions may be approximated following the superposition scheme due to Kachanov [14] for instance and folded into micromechanical damage based models (see e.g. Refs. [2,21,22]). There are at least two outcomes from these approaches: first, the weight function that is introduced in the nonlocal averaging, along with the internal length, is recovered; second this weight function depends on the state of damage and it is direction dependent with respect to the state of stress. Cracks may shield each other or amplify the interaction stresses acting in their neighbourhood.

Nearby the boundary of the solid, interactions between defects are expected to be different compared to those observed in the bulk material. Such boundary effects are among the pending issues in nonlocal modelling for which very little is known from an experimental or a theoretical point of view. In nonlocal models, boundaries are usually dealt with arbitrarily: in integral models the weight function involved in the nonlocal average is chopped off and normalized [20]. It follows that
Nomenclature

<table>
  <tr>
    <td>σ<sub>∞</sub></td>
    <td>remote stress field in finite solid</td>
  </tr>
  <tr>
    <td>σ<sub>ij</sub></td>
    <td>stress components</td>
  </tr>
  <tr>
    <td>ε<sub>kl</sub></td>
    <td>strain components</td>
  </tr>
  <tr>
    <td>C<sub>ijkl</sub></td>
    <td>components of the fourth-order elastic stiffness tensor</td>
  </tr>
  <tr>
    <td>D</td>
    <td>damage variable</td>
  </tr>
  <tr>
    <td>E</td>
    <td>Young's modulus</td>
  </tr>
  <tr>
    <td>ν</td>
    <td>Poisson's ratio</td>
  </tr>
  <tr>
    <td>ε<sub>eq</sub></td>
    <td>equivalent strain</td>
  </tr>
  <tr>
    <td>⟨ε<sub>i</sub>⟩<sub>+</sub></td>
    <td>positive principal strain</td>
  </tr>
  <tr>
    <td>D<sub>t</sub></td>
    <td>tensile damage</td>
  </tr>
  <tr>
    <td>D<sub>c</sub></td>
    <td>compressive damage</td>
  </tr>
  <tr>
    <td>Ψ(x - ξ)</td>
    <td>Ψ<sub>0</sub>(x - ξ), and Ψ<sup>*</sup>(x, ξ) weight functions</td>
  </tr>
  <tr>
    <td>ξ, x</td>
    <td>coordinate system</td>
  </tr>
  <tr>
    <td>ε̄<sub>eq</sub></td>
    <td>nonlocal equivalent strain</td>
  </tr>
  <tr>
    <td>Ω</td>
    <td>volume of the structure</td>
  </tr>
  <tr>
    <td>l<sub>c</sub></td>
    <td>internal length of the nonlocal continuum</td>
  </tr>
  <tr>
    <td>δ(x - ξ)</td>
    <td>Dirac delta function</td>
  </tr>
  <tr>
    <td>a</td>
    <td>minimum between the internal length and the distance from the point to the closest boundary</td>
  </tr>
  <tr>
    <td>b</td>
    <td>minimum between the internal length and the distance to the boundary of the solid in the orthogonal direction</td>
  </tr>
  <tr>
    <td>A<sub>t</sub>, A<sub>c</sub>, B<sub>t</sub>, B<sub>c</sub> and ε<sub>D<sub>0</sub></sub></td>
    <td>parameters in the evolution law for damage</td>
  </tr>
  <tr>
    <td>P<sub>u</sub></td>
    <td>ultimate load in the size effect tests</td>
  </tr>
  <tr>
    <td>σ<sub>N</sub></td>
    <td>nominal strength in size effect tests</td>
  </tr>
  <tr>
    <td>D<sub>b</sub></td>
    <td>parameter in size effect law for unnotched specimen</td>
  </tr>
  <tr>
    <td>f<sub>r∞</sub></td>
    <td>modulus of rupture for specimen of infinite size</td>
  </tr>
  <tr>
    <td>B and D<sub>0</sub></td>
    <td>parameters in size effect law for notched specimens. B depends on the geometry, D<sub>0</sub> is a characteristic size</td>
  </tr>
  <tr>
    <td>G<sub>f</sub></td>
    <td>fracture energy computed according to size effect law</td>
  </tr>
</table>

the influence of a point A located nearby a boundary on a point B located in the bulk of the solid is not the same as the influence of B on A. Due to the truncature of the interaction domain and to the renormalization, the weight function centered at point A and entering in the nonlocal averaging at point A is not the same as the weight function centered at point B and entering in the nonlocal averaging at point B. In the context of continuum damage modelling, this peculiarity of integral non-local models was pointed out many times (see e.g. [19,4]). It is at the origin of the loss of symmetry of the tangent operator in the nonlocal integral formulations. Because the nonlocal interactions are changing nearby the boundary of the solid, the constitutive formulation [20] does not derive from a thermodynamic free energy potential. The modified symmetric nonlocal damage theory due to Borino and coworkers [6] derives from a potential and fulfill thermodynamic principles. To this end, the weight function is modified near the boundary. The background for such a modification is the symmetry of the non-local interactions and energy considerations, it is not related to some specific boundary effect that would arise from the interaction between microcracking and the boundary of the solid.

In gradient enhanced models, the normal component of the gradient of the nonlocal variable is constrained to be zero on the boundary. In fact, the free boundary condition on the nonlocal variable is the same as the condition that would be induced by an axis of symmetry. It means that the nonlocal interactions nearby a boundary of the solid are the same as the nonlocal interactions that would be observed nearby an axis of symmetry. It can be hardly admitted, however, that the interaction between defects and a boundary surface is the same as the interaction between defects distributed symmetrically in a bulk material. Note that in displacement based gradient models [25,11], the displacements ought to be equal to the nonlocal displacements on the boundary of the solid. This is again a different boundary condition.

In any case, there is very little theoretical motivation for such boundary conditions in nonlocal models, at least some justification is lacking. It may be argued that boundary conditions are not very important. Generally, cracks propagate inside the structure and the fracture process zone is located in the bulk material. Initiation of cracking, however, very often occurs from the boundary of a solid. The simplest situation is that of bending beams. It is expected that the boundary effect may have some influence on the initiation condition of cracking and this is the primary motivation for the present study. Once a crack has propagated, it forms a new, evolutive, boundary of the solid and the nonlocal formulation should account for this additional boundary effect. This case is outside the scope of this paper, along with the issue of nonlocal effects nearby interfaces.

Our purpose is to provide some insight on the boundary effects induced by nonlocality and therefore to investigate the effect of subsequent modifications of nonlocal averaging nearby boundaries in integral damage models. Let us start with some intuitive argument about nonlocality nearby the boundary of a solid [23]. Consider a finite body that contains a population of microcracks or microvoids in an elastic matrix. Given a set of boundary conditions, the mechanical response of this body may be described following two techniques: in the first one, the elastic material containing the defects (cracks and voids) is homogenised. The result is a constitutive relation at each material point that depends on the defect density and
geometry. In the next section, we shall follow this track and consider that the influence of microcracking is described according to a continuum damage model.

In the second technique, each defect is described explicitly and the mechanical response of the solid is the result of the deformation of the elastic matrix and the deformation of the defects (e.g. microcrack openings) due to the applied loads. In the course of this calculation, the interaction between the defects is computed (e.g. using the superposition scheme due to Kachanov [14]). This interaction is very often understood as being at the origin of the nonlocality of the constitutive response of the homogenised material described e.g. according to a nonlocal damage model (see for instance Refs. [2,21,22]).

According to Kachanov's superposition scheme, the state of stress and the displacements in this body containing defects is computed as the sum of two sub-problems:

* Sub-problem I: the solid is considered without any defect and it is subjected to the applied boundary conditions. The result is a state of stress denoted as $\sigma_{\infty}$

* Sub-problem II: The solid contains the defects. It is free from any external load/applied displacement. Inside each defect, surface forces are applied so that the stress vectors due to $\sigma_{\infty}$ computed at the imaginary location of the defects in sub-problem I are exactly equilibrated. This is required in order to recover, after superposition, the fact that internal defect surfaces are free of any load.

Usually, this superposition scheme is applied to an infinite body and interaction between defects is considered only, but it can be applied in principle to the case of a finite body too. We shall focus next on sub-problem II, which means that the foregoing discussion is independent from the applied loads or displacements at the boundary of the solid. It is also important to note that in sub-problem II, and for each defect, the surface forces applied inside each defect do not directly equilibrate the stress vector computed at the imaginary location of this defect. The sum of the interaction forces due to the other defects, the interaction forces due to the boundary of the solid, and the applied surface forces inside this defect equilibrate the stress vectors due to $\sigma_{\infty}$ computed at the imaginary location of the defect.

Overall, the boundary of the solid remains always free from any load in sub-problem II. It means the stress vector due to the interactions, computed normal to the boundary, always cancels whatever the distribution of defects. If nonlocality of the constitutive response of the homogenised solid is understood as the consequence of interactions in a continuum setting, nonlocality should vanish on the free surface, in the normal direction to this surface. Close to the boundary, and in the normal direction to the boundary, there should be a layer in which nonlocality increases going farther from the boundary in order to reach the nonlocality expected in an infinite solid.

The purpose of this paper is to analyse the potential structural consequences of this boundary effect on nonlocality. We shall devise a modified phenomenological nonlocal model that satisfies the condition of local material response on the boundary of the solid only, in the normal direction to the boundary. We will then compare the original and modified formulation in structural analyses. A simplified one dimensional setting, a spalling test, will be considered first. Two dimensional examples will also be discussed and attention will be focused on the simulation of size effect bending tests on geometrically similar specimens that are notched and unnotched. We will show that the modified nonlocal model provides a more consistent fit of the size effect laws on both type of specimens than the original nonlocal model.

## 2. Modified nonlocal damage model

We are going to examine now how such results can be incorporated in a continuum formulation. We recall first the main equations involved in the considered damage model.

### 2.1. Local damage model

The classical stress-strain relation for this type of model reads:

$$ \sigma_{ij} = (1 - D)C_{ijkl}\varepsilon_{kl} \quad (1) $$

where $\sigma_{ij}$ and $\varepsilon_{kl}$ are the components of the stress and strain tensors, respectively $(i, j, k, l \in [1, 3])$ and $C_{ijkl}$ are the components of the fourth-order elastic stiffness tensor. The damage variable $D$ represents a measure of material degradation which grows from zero (undamage material with the virgin stiffness) to one (at complete loss of integrity). 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, a scalar equivalent strain $\varepsilon_{eq}$ is introduced, which quantifies the local deformation state in the material in terms of its effect on damage. In this contribution, Mazars' definition of the equivalent strain is used [17]:

$$ \varepsilon_{eq} = \sqrt{\sum_{i=1}^{3} (\langle \varepsilon_i \rangle_{+})^2} \quad (2) $$

where $\langle \varepsilon_i \rangle_{+}$ are the positive principal strains. Damage growth is governed by the loading function:

$$ g(\boldsymbol{\varepsilon}, k) = \varepsilon_{eq}(\boldsymbol{\varepsilon}) - k \quad (3) $$
$k$ equals the damage threshold $\varepsilon_{D_0}$ initially, and during the damage process it is the largest ever reached value of $\varepsilon_{eq}$. The evolution of damage is governed by the Kuhn-Tucker loading-unloading condition:

$$g(\varepsilon, k) \le 0, \quad \dot{k} \ge 0, \quad \dot{kg}(\varepsilon, k) = 0 \qquad (4)$$

The damage variable $D$ is determined as a linear combination of two damage variables $D_t$ and $D_c$, that represent tensile damage and compressive damage, respectively, by the help of two coefficients $\alpha_t$ and $\alpha_c$ which depend on the type of stress state [17]:

$$D = \alpha_t D_t + \alpha_c D_c \qquad (5)$$

$$D_{t,c} = 1 - \frac{1 - A_{t,c}}{\varepsilon_{eq}} + \frac{A_{t,c}}{\exp(B_{t,c}(\varepsilon_{eq} - \varepsilon_{D_0}))} \qquad (6)$$

Standard values of the model parameters in the damage have been given in Ref. [17].

## 2.2. Nonlocal formulation

In the integral-type nonlocal damage models, the local equivalent strain is replaced by its weighted average:

$$\bar{\varepsilon}_{eq}(\mathbf{x}) = \int_{\Omega} \Psi(\mathbf{x}, \xi) \varepsilon_{eq}(\xi) d\xi \qquad (7)$$

with $\Omega$ the volume of the structure and $\Psi(\mathbf{x}, \xi)$ the weight function. It is required that the nonlocal operator does not alter the uniform field, which means that the weight function must satisfy the condition:

$$\int_{\Omega} \Psi(\mathbf{x}, \xi) d\xi = 1 \quad \forall \mathbf{x} \in \Omega \qquad (8)$$

For this reason, the weight function is recast in the following form [20]:

$$\Psi(\mathbf{x}, \xi) = \frac{\Psi_0(\mathbf{x} - \xi)}{\Omega_r(\mathbf{x})} \quad \text{with} \quad \Omega_r(\mathbf{x}) = \int_{\Omega} \Psi_0(\mathbf{x} - \xi) d\xi \qquad (9)$$

where $\Omega_r(\mathbf{x})$ is a representative volume and $\Psi_0(\mathbf{x} - \xi)$ is the basic nonlocal weight function which is often taken as the polynomial bell-shaped function [3], or here as the Gauss distribution function:

$$\Psi_0(\mathbf{x} - \xi) = \exp \left( - \frac{4 \|\mathbf{x} - \xi\|^2}{l_c^2} \right) \qquad (10)$$

$l_c$ is the internal length of the nonlocal continuum. Preserving the uniform field in the vicinity of the boundary makes the averaging in Eq. (9) not symmetric with respect to its arguments $\mathbf{x}$ and $\xi$. This lack of symmetry leads to the nonsymmetry of the tangent operator [4,19,12]. A symmetric nonlocal formulation ensuring enhancement rule has been recently proposed [6]:

$$\Psi(\mathbf{x}, \xi) = \left(1 - \frac{\Omega_r(\mathbf{x})}{\Omega_\infty}\right) \delta(\mathbf{x} - \xi) + \frac{\Psi_0(\mathbf{x} - \xi)}{\Omega_\infty} \qquad (11)$$

where $\delta(\mathbf{x} - \xi)$ is the Dirac function and $\Omega_\infty$ is the representative volume in the infinite solid where it has a constant value. The first term is local, it vanishes for points far from the boundary and the original weight function in Eq. (9) is recovered. According to this modified formulation, the computation of the nonlocal equivalent strain becomes:

$$\bar{\varepsilon}_{eq}(\mathbf{x}) = \varepsilon_{eq}(\mathbf{x}) + \frac{1}{\Omega_{\infty}} \int_{\Omega} \Psi_{0}(\mathbf{x} - \xi)(\varepsilon_{eq}(\xi) - \varepsilon_{eq}(\mathbf{x})) d\xi \qquad (12)$$

The integrand in the right-hand-side term of this equation involves the difference between the local equivalent strain at the considered point and the local equivalent strain in its neighborhood. It is clearly a nonlocal contribution to the quantity (nonlocal equivalent strain) that controls damage.

## 2.3. Modified one-dimensional model

Let us consider first a one dimensional semi-infinite bar. The coordinate system is such that the origin $x = 0$ is at the extremity of this bar. We look at the nonlocal average at point $x$ located close to this extremity and we try to incorporate in the formulation the boundary effect.

A first possibility is inspired from consideration about micromechanics of interacting defects [22,23]. The weight function in Eq. (11) is modified in order to account for the boundary effect. The idea is to extend the solid outside the boundary $x = 0$ with a fictitious bar that has a state of strain symmetric with respect to this boundary. At each point $x$ of the bar, the nonlocal contribution is the classical expression minus a contribution due to $x^*$, where $x^*$ is the symmetric of $x$ with respect to the extremity (boundary) of the bar. The contribution due to $x^*$ is the opposite of the local equivalent strain times the weight of this point. On the boundary, it cancels with that of $x$ and the material response is always local. The resulting new weight function now denoted as $\Psi^*(x, \xi)$ is:
$$ \Psi^*(x, \xi) = \left(1 - \frac{\Omega_r(x)}{\Omega_\infty} + \frac{\Omega_r(x^*)}{\Omega_\infty}\right) \delta(x - \xi) + \frac{\Psi_0(x - \xi)}{\Omega_\infty} - \frac{\Psi_0(x^* - \xi)}{\Omega_\infty} \quad (13) $$

If $x$ is located exactly at the extremity of the bar, the material response is local since $x = x^*$. For a point $x$ sufficiently far from the extremity, the contribution due to $x^*$ vanishes and then this formulation is completely equivalent to the classical nonlocal average. Fig. 1 shows the contribution of the different terms of Eq. (13) for one-dimensional bar in uniform state of damage.

A second possibility, also investigated in this paper, uses directly the original nonlocal formulation of [20] (Eq. (9)), in which the following transformation of the coordinate system defined in the weight function is applied:

$$ \|x - \xi\| \rightarrow \|x - \xi\| \times \frac{l_c}{a} \quad (14) $$

where $a$ decreases as point $x$ is located closer to the boundary. $a$ is the minimum between the internal length and the distance from point $x$ to the extremity of the bar. When $a$ is equal to the internal length, the original formulation is recovered.

This alternative formulation is proposed because its extension to 2D and 3D geometries is simple to handle. The first proposal defined in Eq. (13) requires the definition of the symmetric counterpart of any point inside the structure (and sufficiently close to the boundary) with respect to the boundary of the solid. In the case of a 2D or 3D structures and nearby corners for instance, there are several overlapping symmetric counterparts of a given material point. Thus, this formulation requires tracking all these counterparts depending on the geometry of the structure and it is less easy to implement than the simple modification of the weight function for Gauss points close to the boundary of the solid defined in Eq. (14).

## 2.4. One-dimensional example: spalling test

In order to compare the initial and modified nonlocal formulations, let us consider a one dimensional dynamic tension test (split Hopkinson bar) as shown in Fig. 2. This example has been designed so that the location of localised failure may occur very close to one extremity of the bar, involving modified nonlocal interactions due to the extremity of the bar.

A square compression signal is generated in the bar. Upon reflection at the end of the bar, the compression signal turns into a tensile one. This signal is added to the incoming compression. If the absolute amplitude of the compression signal is greater than the tensile strength, failure is initiated at a distance from the boundary equal to half the signal length. Depending on the duration of the compression signal, it is possible to initiate failure in the material at any location, near the boundary or far from it.

The bar length is taken equal to 25 cm. The parameters used in this example are: the volumic mass $\rho = 1$ kg/m$^3$, the Young's modulus $E = 1$ MPa and the velocity boundary condition $c = 2.5$ mm/s applied at the left bar end. The other model parameters are $A_t = 1$, $B_t = 2$, $\epsilon_{D_0} = 1$ and the internal length $l_c$ is 3 cm (there is no damage in compression). A fixed mesh of 250 constant strain elements is used. Time integration is performed according to an explicit, central difference scheme. The signal length is calculated as $l = t_0\nu$ and its amplitude is $c/\nu$ where $\nu$ is $\sqrt{\frac{E}{\rho}}$.

Fig. 3 shows the case where damage is initiated in the middle of the bar, far away from the boundary. According to the discretisation, the strain inside each element is constant. Integration is performed according to a single Gauss point quadrature, as if the values of damage and nonlocal strain were constant, equal to the value computed at the center of each element. It is these values, constant over each element, that are plotted on the following strain and damage profiles. The results

<figure><img src="image_2.png" /><figcaption>Fig. 1. Diagram showing the contribution of the different terms of the modified weight function (Eq. (13)) for one-dimensional bar in a uniform state of damage. Results are obtained using a Gauss weight function (Eq. (10)) with internal length $l_c = 0.3L$.</figcaption></figure>
obtained with the original and modified formulations are exactly the same. If we consider that complete failure occurs where
damage is maximum, according to Fig. 3, the bar breaks in two pieces at the center.

Let us now study the case where damage initiates close to the boundary at a distance equal to 1.25 cm (the signal length is
2.5 cm). The results with the various weight functions are now quite different. According to the original formulation denoted
here as RWF (Fig. 4), damage is maximum at the extremity of the bar, once it has sufficiently developed. It is in this compu-
tation only that the kinetics of damage evolve during the spall formation as the damage peak in the profiles slowly drifts
towards the extremity of the bar. In the subsequent calculations performed with modified weight functions, the shape of
the damage profile during spall formation does not evolve and the final profiles of damage and nonlocal strain will be plotted
only. With the weight function modified according to [6] denoted as SLC, some improvement is observed as maximum dam-
age is not located at the extremity of the bar (Fig. 5). This is due to the modification of the weight function nearby the
extremity of the bar. The nonlocal component decreases and therefore, the local one becomes more predominant in the for-

<figure><img src="image_3.png" /><figcaption>Fig. 2. Principle of the dynamic tension test.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Damage (top) and strain (bottom) profiles with the original and modified models falling all in one plot when failure occurs at the center of the bar.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Damage (top) and strain (bottom) profiles for RWF when damage initiates nearby the extremity of the bar.</figcaption></figure>

mulation. Still, this case is not consistent with the considerations from crack interaction in Section 2. At the extremity of the bar, some nonlocality remains and the material response is not local. With the modified formulations denoted as MCLS and ERF (Figs. 6 and 7) for Eqs. (13) and (14), respectively, damage is maximum inside the bar, where tension initiates and there is no damage at the extremity.

Let us now define the thickness of the spall as the distance between the right-hand-side extremity of the bar and the closest point at which damage is equal to one (complete failure). The thickness of the spall is infinitesimally small according to the original RWF approach because maximum damage occurs at the extremity of the bar. As a matter of fact, if one wants to obtain a spall that is of the order of the internal length of the material or less, it is not possible with the RWF formulation. Damage develops initially with a maximum inside the bar but as it grows, the maximum of damage is attracted towards the extremity of the bar. As shown in Figs. 6 and 7, this is not the case with the modified nonlocal formulations mainly because interactions decrease nearby the extremity of the bar. According to the ERF and MCLS formulations of the weight function, the spall is approximately 1 cm thick. The two modified formulations provide very similar results. The thickness of the spall is approximately equal to one third of the internal length.

## 2.5. Modified 2D formulation

The ERF modified weight function is quite easy to extend to 2D calculations. The modified weight function reads:

$$ \|\mathbf{x} - \xi\| = l_c \sqrt{(x_1 - \xi_1)^2/a^2 + (x_2 - \xi_2)^2/b^2} \qquad (15) $$

in a coordinate system where subscript 1 refers to a vector that is normal to the closest boundary of the solid and subscript 2 refers to the orthogonal direction. a is the minimum between the internal length and the distance from the point to the closest boundary and b is assumed to be the minimum between the internal length and the distance to the boundary of the solid in the orthogonal direction. The weight function in the three dimensional case can also be extended following the same line,
<figure><img src="image_6.png" /><figcaption>Fig. 5. Damage (top) and strain (bottom) profiles for SLC when damage initiates nearby the extremity of the bar.</figcaption></figure>

without difficulties. From now on, we shall focus on this formulation and compare results with the original RWF nonlocal formulation.

## 3. Two dimensional examples

### 3.1. 2D plate in tension

Let us first discuss the simulation of a direct tension test of a plate (1 m × 1 m) that contains defects in order to trigger strain localisation. We look here at a double edge notched specimen subjected to a remote applied load (Fig. 8).

The numerical analysis has been carried out under the plane stress assumption. A regular mesh of 25 × 25 quadratic elements is used. The following parameters are used in this analysis: $E = 3.5 \times 10^{10}$ Pa, $v = 0.2$, $A_t = 0.95$, $B_t = 9000$, $A_c = 1.25$, $B_c = 1000$, $\varepsilon_{D_0} = 1.0 \times 10^{-4}$. The internal length $l_c$ is taken equal to 0.2 m. The two square notches are small compared to the internal length (four times smaller than the internal length). They are modelled as completely damaged finite elements.

The maps of damage with the two approaches are shown in the Fig. 9. With the two models (RWF and ERF), damage initiates nearby the tip of the notches, and then propagates towards the centre until failure occurs. The differences between the damage maps are nearby the notches only. With the classical nonlocal formulation RWF, the material located between the tip of the notch and the boundary is damaged completely. This is not the case with the modified model (ERF) because the interactions decrease close to the boundary.

There are rather small differences between the two damage maps and consequently it yields small differences on the corresponding load-displacement responses. There are, however, at least two cases where we expect a significative difference between the RWF and the ERF approaches. The first one is the case of a crack propagating towards a boundary. The fracture process zone (FPZ) is constrained by the boundary of specimen [8] and the energy required to propagate the crack should
<figure><img src="image_7.png" /><figcaption>Fig. 6. Damage (top) and strain (bottom) profiles for MLCS when damage initiates nearby the extremity of the bar.</figcaption></figure>

change. The specimen boundary will limit the development of the FPZ, and therefore, leads to a reduced fracture energy and strength [8, 7, 9, 10].

The second case is typically found in comparing specimens where a crack initiates from a surface with different geometries. The condition of onset of propagation should depend on the geometry of this boundary. We will focus in the following on this case, and more specifically on the prediction of size effect of geometrically similar notched and unnotched bending beams.

## 3.2. Three point bending size effect tests

A salient characteristics of nonlocal modelling is structural size effect (understood here as the dependence of nominal strength on the structural size). The response of geometrically similar specimens is not geometrically similar. The size of the FPZ is controlled by the internal length [18]. When the size of the structure changes, the ratio between the size of the FPZ, which is constant, and the size of the structure is changing too. This produces differences in energy release during the propagation of a crack and size effect on the nominal strength. Detailed explanations can be found in the textbook by Bažant and Planas [5], along with comparisons with experimental data on quasi-brittle materials such as concrete.

Our aim is here to compare the original (RWF) and modified (ERF) nonlocal damage formulations. We will start first with unnotched specimens and consider three point bending tests with three geometrically similar sizes (see Fig. 10). The specimens with various height $D = 80, 160, 320$ mm are referred to as small, medium and large beam, respectively. Simulations are 2D plane stress calculations. The element size was kept constant in the FPZ, and small enough compared to the internal length (at least three times smaller). The model parameters used for these simulations are: $E = 3.85 \times 10^4$ MPa, $v = 0.24$, $A_t = 0.95$, $A_c = 1.25$, $B_t = 9200$, $B_c = 1000$ and $\varepsilon_{D_0} = 3.0 \times 10^{-5}$. The internal length is equal to 10 mm.

### 3.2.1. Size effect on unnotched specimens

The nominal strength now is computed according to the formula:
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Damage (top) and strain (bottom) profiles for ERF when damage initiates nearby the extremity of the bar.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. 2D plate: geometry, loading setup and location of the notches in grey.</figcaption></figure>

$$ \sigma_N = \frac{9 P_u}{2 bD} \qquad (16) $$

It is the maximum tensile stress computed at peak load according to the elastic beam theory. The ultimate loads $P_u$ for the three sizes with the two approaches and the corresponding nominal strength are listed in Table 1.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Damage profiles with (a) RWF and (b) ERF.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Three point bending test: geometry and loading.</figcaption></figure>

For this set of specimens, we have also run the calculation with an internal length of 40 mm. The results are shown in Table 2.

We notice that these results are quite close to each other. As pointed out by Bazant, there is first a layer of distributed damage that forms at the bottom of the beam in the tensile part. At peak load, a crack forms in the center of the beam. The crack is perpendicular to the distributed damage layer, it starts from this layer, inside the beam and not directly from the boundary. Therefore, the effect of the boundary on the inception of the crack is small and peak loads are similar.

Let us now use Bazant's size effect law [5] for the case of unnotched beams:

$$ \sigma_N = f_{r\infty} \left( 1 + \frac{D_b}{D} \right) \quad (D \gg D_b) \qquad (17) $$

where $D_b$ and $f_{r\infty}$ are constants. The first constant is the thickness of the boundary layer in which damage is distributed prior to localised crack inception and the second one is the strength for a specimen of infinite size. This is the modulus of rupture

<table><caption>**Table 1**

Numerical results for three different sizes with RWF and ERF ($l_c$ = 10 mm).</caption><thead><tr><th>Size</th><th>D (mm)</th><th>P<sub>u</sub> (kN) (RWF)</th><th>σ<sub>N</sub> (MPA) (RWF)</th><th>P<sub>u</sub> (kN) (ERF)</th><th>σ<sub>N</sub> (MPA) (ERF)</th></tr></thead><tbody><tr><td>Small</td><td>80</td><td>65.732</td><td>3.70</td><td>64.360</td><td>3.62</td></tr><tr><td>Medium</td><td>160</td><td>123.792</td><td>3.48</td><td>122.588</td><td>3.45</td></tr><tr><td>Large</td><td>320</td><td>240.440</td><td>3.38</td><td>240.460</td><td>3.38</td></tr></tbody></table>
of an infinitely large specimen. $D_b$ and $f_{r\infty}$ are obtained by fitting our computed values of the nominal strengths with the size effect law (Table 3).

One can observe in this table that $f_{r\infty}$ is almost constant and does not depend on the internal length. The tensile strength of the material is, according to the constitutive relations, equal to 2 MPa. Here, the modulus of rupture $f_{r\infty}$ computed from the FE analyses of specimens of different sizes is approximately 1.65 times the tensile strength. This is close to the usual ratio of 1.5 times the tensile strength. According to dimensionnal analysis, $D_b$ should proportional to the internal length [5]. This proportionality is not exactly recovered in Table 3 for the two considered weight functions. This is certainly due to approximation errors that are inherent to the finite element discretisation. In the original RWF formulation, $D_b$ is approximately equal to the internal length. Due to the boundary effect, it decreases in the ERF formulation. This is again consistent with the fact that nonlocal effects are decreased nearby the bottom face of the beam.

### 3.2.2. Size effect on notched specimens

The issue is now failure at crack initiation from a notch. The notch located at mid-span is 0.2D high. In a local computation, the stress amplification due to the notch is much greater than in a nonlocal damage calculation (due to the averaging process in the later). Hence, the stress amplification should be greater near the notch for the modified (ERF) nonlocal model than for the original (RWF) one and a decrease of the load carrying capacity is expected. We shall focus here on computations with an internal length of 10 mm but trends are the same for other values of the internal length.

The ultimate loads $P_u$ for the three sizes with the two approaches are listed in the Table 4. The peak loads for the three sizes are quite different according to the two nonlocal formulations. At the peak, damage occurs nearby the tip of the notch only and this is where the modification of the nonlocal model is important. Such results, showing a great sensitivity of the mechanical response of specimens to the weight function, have already been pointed out, e.g. by Jirasek and coworkers [13].

For each size, the nominal strength, as defined by Bažant and Planas [5], is computed. It is the maximum tensile stress in the beam computed at the tip of the notch according to the elastic beam theory:

$$ \sigma_N = \frac{9}{2} \frac{P_u}{Db(0.8^2)} \qquad (18) $$

The numerical results are now interpreted with Bazant's size effect law for notched beams [1]:

<table><caption>**Table 2**
Numerical results for three different sizes with RWF and ERF ($l_c = 40$ mm).</caption><thead><tr><th>Size</th><th>D (mm)</th><th>P<sub>u</sub> (kN) (RWF)</th><th>σ<sub>N</sub> (MPa) (RWF)</th><th>P<sub>u</sub> (kN) (ERF)</th><th>σ<sub>N</sub> (MPa) (ERF)</th></tr></thead><tbody><tr><td>Small</td><td>80</td><td>87.296</td><td>4.91</td><td>79.242</td><td>4.46</td></tr><tr><td>Medium</td><td>160</td><td>144.488</td><td>4.06</td><td>140.486</td><td>3.95</td></tr><tr><td>Large</td><td>320</td><td>262.260</td><td>3.69</td><td>256.520</td><td>3.61</td></tr></tbody></table>

<table><caption>**Table 3**
$f_{r\infty}$ and $D_b$ for two internal lengths, original and modified weight functions.</caption><thead><tr><th></th><th>RWF<br>l<sub>c</sub> = 10 mm</th><th>ERF<br>l<sub>c</sub> = 10 mm</th><th>RWF<br>l<sub>c</sub> = 40 mm</th><th>ERF<br>l<sub>c</sub> = 40 mm</th></tr></thead><tbody><tr><th>f<sub>r\infty</sub></th><td>3.27</td><td>3.30</td><td>3.26</td><td>3.35</td></tr><tr><th>D<sub>b</sub></th><td>10.34</td><td>7.82</td><td>40.16</td><td>26.62</td></tr></tbody></table>

<table><caption>**Table 4**
Numerical results for three different sizes of the notched beam with RWF and ERF ($l_c = 10$ mm).</caption><thead><tr><th>Size</th><th>D (mm)</th><th>P<sub>u</sub> (kN) (RWF)</th><th>σ<sub>N</sub> (MPa) (RWF)</th><th>P<sub>u</sub> (kN) (ERF)</th><th>σ<sub>N</sub> (MPa) (ERF)</th></tr></thead><tbody><tr><td>Small</td><td>80</td><td>42.316</td><td>3.72</td><td>36.910</td><td>3.24</td></tr><tr><td>Medium</td><td>160</td><td>64.426</td><td>2.83</td><td>57.786</td><td>2.54</td></tr><tr><td>Large</td><td>320</td><td>97.454</td><td>2.14</td><td>90.586</td><td>1.99</td></tr></tbody></table>

<table><caption>**Table 5**
$B_{f_{r\infty}}$ and $D_0$ of the notched beam of different sizes for $l_c = 10$ mm.</caption><thead><tr><th></th><th>RWF</th><th>ERF</th></tr></thead><tbody><tr><th>B<sub>f<sub>r\infty</sub></sub></th><td>6.25</td><td>4.64</td></tr><tr><th>D<sub>0</sub></th><td>42.42</td><td>71.42</td></tr></tbody></table>
$$ \sigma_N = B f_{r\infty} (1 + D/D_0)^{-1/2} \qquad (19) $$

where B is a dimensionless geometry-dependent parameter and $D_0$ is a characteristic size. For each formulation, $D_0$ and $Bf_{r\infty}$, identified from a linear regression as explained in [5] are reported in Table 5.

There is a decrease of the maximum carrying capacity with the ERF formulation compared to the original nonlocal formulation. This decrease results in the size effect law into a decrease of $Bf_{r\infty}$, that is a decrease of $f_{r\infty}$ because B is constant and related to the geometry of the specimen only.

We may now compare the values of the modulus of rupture $f_{r\infty}$ obtained from size effect on unnotched specimens to those obtained with notched specimens. For that B is calculated first according to Rilem recommandations [24]. We obtain $B = 1.11$. Then, we take $f_{r\infty}$ obtained for the unnotched specimens, multiply by B, and compare the result with the value of the fit obtained from size effect on notched specimens. Ideally, the two values should be the same. Table 6 shows that the relative error on the prediction of $Bf_{r\infty}$ is three times smaller with the ERF formulation than with the RWF one. In view of the high sensitivity of the size effect parameters to the values of the peak loads [15], an error of the order of 30% is quite reasonable and the modified ERF formulation provides consistent results for both type of size effects whereas it is not the case for the original nonlocal formulation.

We may also observe in Table 5 that $D_0$ is increasing in the ERF calculations compared to the RWF calculations. It should be pointed out that the size effect formula in Eq. (19) holds for sizes that are large compared to $D_0$. D is close to $D_0$ and we may consider that it is more appropriate to implement here the universal size effect formula proposed by Bazant instead of two separate formulae for notched and unnotched specimens:

$$ \sigma_N = B f_{r\infty} (1 + D/D_0)^{-1/2} \cdot 1 + \left[ \left( \eta + \frac{D}{D_b} \right) \cdot (1 + D/D_0) \right]^{-1} \qquad (20) $$

where $\eta$ is taken here equal to 1. In order to obtain the parameters in this formula, we consider first the case of unnotched specimens. Eq. (20) reduces exactly to Eq. (17) and the corresponding constants can be fitted. We use then the data computed for notched specimens and the corresponding value of B computed above and we look for the value of $D_0$ that provides the best fit with the data. Fig. 11 shows these fits obtained for both nonlocal formulations.

These fits are obtained for $D_0 = 210$ mm with the RWF formulation and $D_0 = 140$ mm with the ERF formulation. Note that the agreement is again better with the modified nonlocal formulation than with the original nonlocal one. Finally, it is possible to compute the fracture energy from the size effect law.

$$ G_f = B f_{r\infty}^2 \cdot D_0 \cdot g/E \qquad (21) $$

<table><caption>**Table 6**
$Bf_{r\infty}$ computed from notched and unnotched specimens.</caption><thead><tr><th>Bf<sub>r∞</sub></th><th>RWF</th><th>ERF</th></tr></thead><tbody><tr><td>Computed from unnotched</td><td>3.63</td><td>3.66</td></tr><tr><td>Fit from notched</td><td>6.25</td><td>4.64</td></tr><tr><td>Relative error</td><td>72%</td><td>27%</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Fits of the universal size effect law on notched specimens according to the RWF and ERF formulations.</figcaption></figure>
where *E* is the Young's modulus and *g* is related to the geometry of the specimen (here *g* = 1.1116 according to [24]). For the *RWF* formulation, we obtain *G*<sub>f</sub> = 80 N/mm<sup>2</sup> and for the *ERF* formulation, we obtain *G*<sub>f</sub> = 54 N/mm<sup>2</sup>. The fracture energies computed according to the original and modified nonlocal formulation, and defined according to the size effect method, are very different.

## 4. Conclusions

We have presented two tentative modifications of the original nonlocal damage model that are consistent with the intuitive argument that nonlocality should vanish at the boundary of a solid (in the direction normal to the boundary) and should grow getting inside the solid. These local models satisfy the condition of local material response on the boundary of the solid, and in the direction normal to the boundary. The first one is a modification of the weight function nearby the boundary and the second one remaps the coordinate system defined in the weight function when the average is computed close to the boundary. The two formulations provide, in a 1D example at least, results that are very similar.

In the one-dimensionnal spalling test, we show that the modified formulations provide more realistic results with a spall of finite, nonzero, thickness. This cannot be achieved with the original nonlocal damage formulation. Even though we devised slight modifications of the nonlocal model, two dimensional computations show that some nonnegligible differences may be observed between the original and modified ones. It corresponds to cases where damage initiates close or at the boundary of the solid. This is in fact rather typical of many failure problems.

The simulation of notched and unnotched specimens of different sizes have been interpreted with the help of Bazant's size effect laws. With the modified *ERF* formulation, it is possible to fit data for notched and unnotched bending beams with the same set of model parameters including the modulus of rupture. It is not possible to achieve the same results with the original *RWF* formulation. The modulus of rupture needs to be changed in a range beyond possible calibration errors. The fracture energies are also quite different with the modified and original nonlocal damage models.

A problem is now to be able to discriminate the two nonlocal formulations on the basis of experimental evidences. This is beyond the scope of the present paper, nevertheless we may provide some preliminary ideas. A first possibility is to reproduce experimentally the spalling test modelled in the one-dimensionnal study in this paper. The calibration of test data from notched and unnotched specimens, and the comparisons of the parameters in the size effect laws for the two types of specimens, same as in the simulations discussed in this paper, is an indirect technique that could be followed, due to the simplicity of the experiments to be performed.

Further, several theoretical issues might be considered: a first one is the derivation of the same type of model within the gradient approximation and how should existing models, e.g. second order gradient damage models, be modified in order to account for this boundary effect. Dirichlet conditions might be applied at the boundary, better on the component that is normal to the boundary, but consequences should be assessed. In models based on nonlocal displacements [25,11], nonlocal displacement and local displacements are equal on the boundary of the solid. Proper descriptions of size effects for notched and unnotched should be checked. Finally, issues related to the shape of the stress distributions at the tip of notches and crack initiation could be also investigated in order to see whether the modified model provides more consistent results than the original one [26].

## References

[1] Bažant ZP. Size effect in blunt fracture: concrete, rock, metal. J Engng Mech 1984;110:518-35.

[2] Bažant ZP. Nonlocal damage theory based on micromechanics of crack interactions. J Engng Mech 1994;120:593-617.

[3] Bažant ZP, Jirasek M. Nonlocal integral formulations for plasticity and damage: survey of progress. J Engng Mech 2002;128:1119-49.

[4] Bažant ZP, Pijaudier-Cabot G. Nonlocal continuum damage, localization instability and convergence. J Appl Mech ASME 1988;55:287-94.

[5] Bažant ZP, Planas J. Fracture and size effect in concrete and other quasi-brittle materials, Boca Raton and London; 1998.

[6] Borino G, Failla B, Parinello F. A symmetric nonlocal damage theory. Int J Solids Struct 2003;40:3621-45.

[7] Duan K, Hu X. Specimen boundary induced size effect on quasi-brittle fracture. Strength Fract Complexity 2004;2(2):47-68.

[8] Duan K, Hu X, Wittmann H. Boundary effect on concrete fracture and non-constant fracture energy distribution. Engng Fract Mech 2003;70:2257-68.

[9] Duan K, Hu X, Wittmann H. Scaling of quasi-brittle fracture: boundary and size effect. Mech Mater 2006;38:128-41.

[10] Hu X, Duan K. Influence of fracture process zone height on fracture energy of concrete. Cem Concr Res 2004;34:1321-30.

[11] Jirasek M, Marfia S. Non-local damage model based on displacement averaging. Int J Numer Methods Engng 2005;63(1):77-102.

[12] Jirasek M, Patzàk B. Consistent tangent stiffness for nonlocal damage models. Comput Struct 2002;80:1279-93.

[13] Jirasek M, Rolshoven S, Grassl P. Size effect on fracture energy induced by non-locality. Int J Numer Anal Methods Geomech 2003;28:635-70.

[14] Kachanov M. Elastic solids with many cracks – a simple method of analysis. Int J Solids Struct 1987;23:23-43.

[15] Le Bellégo C, Dubé J-F, Pijaudier-Cabot G, Gérard B. Calibration of nonlocal damage model from size effect tests. Eur J Mech A/Solids 2003;22:33-46.

[16] LeBlond J, Perrin G, Devaux X. Bifurcation effects in ductile metals with damage delocalisation. J Appl Mech 1994;61:236-42.

[17] Mazars J. Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure. Ph.D. thesis, Université Pierre et Marie Curie; 1984.

[18] Mazars J, Pijaudier-Cabot G. From damage to fracture mechanics and conversely: a combined approach. Int J Solids Struct 1996;33:3327-42.

[19] Pijaudier-Cabot G. Non local damage. In: Muhlhaus HB, editor. Continuum models for materials with microstructure. John Wiley PubS.; 1995.

[20] Pijaudier-Cabot G, Bažant Z. Nonlocal damage theory. J Engng Mech 1987;113:1512-33.

[21] Pijaudier-Cabot G, Berthaud Y. Effets des interactions dans l'endommagement d'un milieu fragile. Formulation non locale. Comptes rendus de l'AcadTmie des sciences 1990;310:1577-82.

[22] Pijaudier-Cabot G, Haidar K, Dubé J-F. Non-local damage model with evolving internal length. Int J Numer Anal Methods Geomech 2004;28:633-52.
[23] Pijaudier-Cabot G, Krayani A, Dufour F. Comments on boundary effects in nonlocal damage based models. In: Yuan H, Wittmann FH, editors. Nonlocal modelling of failure of materials. Aedificatio Pubs.; 2007.

[24] Rilem. Recommendation. Mater Struct 1990;23:461-5.

[25] Rodriguez-Ferran A, Morata I, Huerta A. A new damage model based on non-local displacements. Int J Numer Anal Methods Geomech 2005;29:473-93.

[26] Simone A, Askes H, Sluys L. Incorrect initiation and propagation of failure in non-local and gradient-enhanced media. Int J Solids Struct 2004;41:351-63.