# Comptes Rendus
Physique

Gilles Pijaudier-Cabot

**Fracture and permeability of concrete and rocks**

Volume 21, issue 6 (2020), p. 507-525

Published online: 14 January 2021

Issue date: 14 January 2021

https://doi.org/10.5802/crphys.38

**Part of Special Issue:** Prizes of the French Academy of Sciences 2019 (continued)

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/

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

Les Comptes Rendus. Physique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org

e-ISSN : 1878-1535
<figure><img src="image_3.png" /></figure>

Prizes of the French Academy of Sciences 2019 (continued) / *Prix 2019 de l'Académie des sciences (suite)*

# Fracture and permeability of concrete and rocks

## *Rupture et perméabilité des bétons et des roches*

Gilles Pijaudier-Cabot<sup>a, b</sup>

**Abstract.** Continuum Damage Mechanics provides a framework for the description of the mechanical response of concrete and rocks which encompasses distributed micro-cracking, macro-crack initiation, and then its propagation. In order to achieve a consistent setting, an internal length needs to be introduced to circumvent the difficulties inherent to strain softening and to avoid failure without dissipation of energy. Upon inserting this internal length, structural size effect is captured too. This paper reviews some the progresses achieved by the author since the introduction of the nonlocal damage model in 1987. Among them, the early proposals exhibited a proper description of the inception of failure but a poor one for complete failure since it is not straightforward to model a discrete cracking with a continuum approach. Candidate solutions, e.g. by considering a variable internal length are outlined. Then, the coupled effects between material damage and material permeability are considered. Is is recalled that the permeability of the material should be indexed on the damage growth in the regime of distributed cracking. Upon macro-cracking, there is a change of regime and it is the crack opening that controls the fluid flow in the cracked material. Both regimes may be captured with a continuum damage approach, however.

**Résumé.** La mécanique de l'endommagement fournit un cadre qui permet de décrire l'ensemble du processus de rupture d'un matériau quasi-fragile sollicité par un chargement mécanique, à savoir une micro fissuration distribuée tout d'abord, puis l'amorçage et la propagation d'une macro-fissure. Une longueur interne doit être introduite afin d'obtenir une énergie dissipée non nulle à la rupture. Cette longueur interne induit un effet de taille cohérent lui aussi avec les données expérimentales. Cet article passe en revue quelques des progrès réalisés par l'auteur depuis l'introduction du modèle d'endommagement non local en 1987. Parmi ceux-ci, les premiers modèles permettaient de bien décrire l'amorçage de la rupture mais moins bien la rupture complète, qui est une chose peu naturelle dans le contexte d'une description continue d'un solide. Des solutions possibles à ce problème, par exemple en faisant varier la longueur interne, sont évoquées. Puis, les effets couplés entre l'endommagement et la perméabilité d'un matériau sont abordés. La perméabilité du matériau doit être indexée sur la croissance des dommages dans le régime de fissuration distribuée. Lors de la macro-fissuration, il y a un changement de régime et c'est l'ouverture de fissure qui contrôle l'écoulement du fluide dans le matériau fissuré. Cependant, ces deux régimes peuvent être décrits avec un formalisme unique basée, au plan mécanique, sur la mécanique de l'endommagement.

<footnote><sup>a</sup> Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, Total, LFCR, Anglet, France</footnote>

<footnote><sup>b</sup> Institut Universitaire de France, France</footnote>

<footnote>*E-mail:* Gilles.Pijaudier-Cabot@univ-pau.fr</footnote>
**Keywords.** Damage, Cracking, Permeability, Size effect, Strain softening, Strain localisation, Internal length.

**Mots-clés.** Endommagement, Fissuration, Perméabilité, Effet d'échelle, Adoucissement, Localisation, Longueur interne.

# 1. Introduction

Fracture Mechanics—the theory that governs the propagation of cracks in solids—is over a century old. Pioneering motivations originated from the observation that, for crystals, the material strength was an order of magnitude less experimentally than that predicted by strength of material theories. The presence of flaws and cracks in the material, prior to any load, could explain such a mismatch. The playground of fracture mechanics has been initially that of metal alloys. Analytical and later on numerical models have been developed extensively, with applications to crack propagation under static, dynamic and fatigue loads. These developments were based upon the description of a crack in a solid as a displacement discontinuity resulting into an eigen state of stress in elasticity [1]. Because calculations show that the stress is infinite at the tip of the discontinuity, the crack propagation conditions could not be based on a stress-based criterion, but rather on a quantity called fracture toughness (see e.g. [2]), or on a fracture energy defined as the energy consumed in order to propagate a crack over an area of unit surface [3]. This theory relies on the assumption that the material is elastic and perfectly brittle and defines the conditions of propagation of the crack. It can be also easily extended to the case where the material is not elastic, typically if it has entered in the plastic regime in an infinitely small region about the crack tip (under the so-called small-scale yielding assumption).

In the 1970's, linear elastic fracture mechanics started to be applied to geomaterials, concrete and rocks, at first with limited success because these materials are not brittle, but quasi-brittle. The degradation processes involved during the crack propagation occur in such materials over a zone called the Fracture Process Zone (FPZ). Due to the size of material heterogeneities (grains which maybe up to several centimeters as opposed to the micron scale in alloys), this zone is much larger in concrete and rocks than in alloys and its size can no longer be considered as infinitely small. Consequently, the energy consumed upon crack propagation is not solely due to the propagation of the tip of macro-crack, it is also due to the propagation of the FPZ that surrounds the macro-crack.

In order to illustrate this important feature, let us consider the response of notched beams made of mortar under three-point bending loads (Figure 1). The mechanical set-up consists of a three-point-bend loading system controlled by a Crack Mouth Opening Displacement (CMOD) sensor. This sensor is a clip gage attached on both sides of the notch. The load is adjusted so that the CMOD is monotonically increasing with time. At the same time, the deflection at mid-span is measured with a laser extensometer. The deflection is the vertical displacement of a point located at mid-span and mid-height on the lateral surface of the beam. Typical mechanical responses will be illustrated in Section 3, let us focus here on data obtained from acoustic emission during the loading.

Acoustic Emission (AE) transducers are placed on the surface of the specimen in order to localise the acoustic emissions that are generated by the occurrence of micro-cracks. Indeed, upon their opening, micro-cracks generate a stress wave (like tiny earthquakes) that can be recorded by piezoelectric sensors. If three sensors at least are used and if the wave speed in the
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Three-point bending experiment on mortar beams. Mechanical set-up (left) and acoustic emission set-up (right).</figcaption></figure>

material is known, the records of the time of arrival of the stress wave on the sensors provide the location of the acoustic event by simple triangularisation.

In the present experiments [4], the accuracy of the localization is ±4 mm. Figure 2 shows the cumulative location of acoustic emission events throughout the test; the plotted points indicate the detected AE sources, we have also plotted the observed crack path that appeared after the test at complete failure. It is clear that the FPZ located around the crack (dark line) is large. Its width is about 7 cm here. The resulting fracture energy is no longer related to the propagation of a single discontinuity, but rather to the propagation of the FPZ.

The cohesive crack model was developed for this purpose by Barenblatt [5] and Dugdale [6] for ductile materials, and applied by Hillerborg et al. [7] to concrete. The theory assumes that the width of the FPZ (measured orthogonally to the crack surface) collapses onto the crack surface. Ahead of the crack tip, the cohesive zone follows a stress-relative displacement (or crack opening equivalently) softening relationship.

Fracture mechanics describes the conditions of crack propagation, which must be complemented with a criterion for crack initiation and crack orientation. Continuum Damage Mechanics (CDM) provides a broader framework. It covers distributed cracking, macro-crack initiation, and then its propagation within a single continuum setting. As we will see in Section 2, progressive micro-cracking is captured via the degradation of the elastic constants of the material. Its application to quasi-brittle materials dates back to the 1980's [8], but in order to achieve a consistent description of the fracture process, it is necessary to introduce, aside from the pointwise stress-strain constitutive relationship, an internal length that defines indirectly the size of the FPZ. This yielded the so-called nonlocal damage model [9] which opened the path to a wide variety of constitutive models, denoted as nonlocal continua with local strains [10]. Section 3 presents the non-local damage models and discusses the variety of possible localisation limiters that arose since then.

The issue of fracture of geomaterials should not be regarded as completely solved however. A typical example is encountered when considering coupled effect between fracture and permeability (e.g. those involved in hydraulic fracturing of rocks). Section 4 discusses this problem.

## 2. Continuum modelling of progressive cracking

The issue addressed in this section is the following: in a continuum setting, i.e. in a framework where the response of a material to a mechanical load is described by a local (pointwise) constitutive relation that relates the strain (derivative of displacements) to the stress (force
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Cumulative location of AE events in a three-point bending test (after [4]). The large circular dots represent the sensors and the black line the macro-crack.</figcaption></figure>

per unit surface of material), how should the progressive micro-cracking of the material be described? Typically, this question is related to two issues:

* First, it questions the equivalence between a solid that contains micro-cracks and a homogeneous solid with elastic properties that depend on the amount of micro-cracking. This equivalence can be performed with the help of homogenisation.

* Second, it is required to formalise how the response is, or not, irreversible due to the growth of micro-cracking: given a state of stress (or strain) and micro-cracking, it means that for a loading increment where micro-cracking progresses, the incremental mechanical response should be different compared to the case where the loading increment is such that micro-cracking does not change. This will be performed with a framework that is similar to elasto-plasticity (see e.g. [11]).

Let us first consider an approach based on elastic-brittle lattices where the second issue is captured with a simplified material response: lattice elements are elastic, up to a threshold strain where breakage occurs. Before the threshold, the element is elastic, after, it can no longer carry a load.

## 2.1. Failure in elastic-brittle disordered lattices

The lattice described in the following was used in the past for the study of the failure mechanism of quasi-brittle materials. The mechanical problem is substituted by an electrical analogue, simplifying the analysis within each lattice element from vectorial quantities (forces, displacements)
<figure><img src="image_6.png" /><figcaption>Figure 3. Lattice model.</figcaption></figure>

to scalar quantities (current, voltage difference). From the mechanical point of view, it is as if the lattice would be made of bars instead of beams. In the case of a lattice of infinite size and prior bond breakage, it is a discrete approximation of an elastic continuum with a fixed Poisson's ratio (see e.g. [12]).

The model, depicted in Figure 3, is a regular two-dimensional lattice. The lattice size is $L \times L$ where $L$ is related to the total number of bonds $n = 2L^2$. The strain $\epsilon$ is substituted with the voltage difference $v$, the stress $\sigma$ with the current $i$, and the Young's Modulus $E$ (material stiffness) with the conductance $G$. So in elasticity, we substitute the mechanical constitutive equation

$$ \sigma = E \cdot \epsilon \qquad (1) $$

with:

$$ i = G \cdot v. \qquad (2) $$

Periodicity is imposed along the boundaries and a constant unit jump of voltage is applied between the horizontal boundaries, corresponding to an imposed relative displacement in the vertical direction.

Every bond of the lattice behaves in the electric problem as an analogy with an elastic-brittle material. The conductance is equal to 1 and when it reaches a threshold current $i_c$, it falls to zero. $i_c$ differs from bond to bond following a uniform random distribution between 0 and 1. Figure 4 shows a typical current vs. voltage response.

The lattice does not represent a real material. In fact, we are interested in its scaling properties: as the lattice size tends to infinity, its response tends to a thermodynamic limit which is the response of a single point in a continuum approach. Hence, we shall look at the variables that describe the evolution of the lattice response with bond breakage, independently from its size. These variables are those which appear in a continuum formulation. For this, we introduce the moments calculated according to the formula:

$$ M_m = \int i^m N(i) di \qquad (3) $$

where $i$ is the current of each bond, $N(i)$ is the number of bonds whose current is $i$, and $m$ is the order of moment. The zero-order moment is the number of unbroken bonds, the first order moment is proportional to the average value of the current (stress). The second order moment is proportional to the overall conductance of the lattice (stiffness in the mechanical problem).

Figure 5 represents the average moments for different lattice sizes as a function of the moment of order 2 until the peak current. We observe that until the peak, the curves of $M_0, M_1, M_2$ and $M_3$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Typical intensity (load) versus voltage (relative displacement) response of a lattice of size 64 (after [13]).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Opposite of the logarithm of the moments of order 0 to 4 as a function of the opposite of the logarithm of the moment of order 2. For each order of moments, results for lattice sizes ranging from 16 to 128 are superimposed (after [13]).</figcaption></figure>

versus $M_2$ are almost independent from the size of the lattice. Hence, the average conductance is a parameter that may describe the distribution of current whatever the size of the lattice. It was also shown that the number of broken bonds could not characterise the state of damage in a size independent way [13].

## 2.2. Continuum damage

Because the average conductance (or lattice stiffness in the mechanical problem) is the parameter which represents local bond breakage in the continuum approach, the effect of micro-
cracking on the mechanical response of a material viewed as homogeneous may be properly captured with a degradation of its elastic stiffness. This is exactly the purpose of continuum damage mechanics. For the sake of simplicity (and for the sake of avoiding complex mathematical developments), we present in this section a one dimensional continuum damage model that can easily be extended to a 3D setting [8]. The constitutive relation given in (1) in the case of elasticity becomes:

$$ \sigma = (1 - D) \cdot E \cdot \varepsilon \quad (4) $$

where D is the damage variable ranging from 0 for a material without micro-cracks to 1 for a material that is totally cracked. Here the damage variable is a scalar, in 3D and if it is still a scalar, damage is isotropic, meaning that the degradation of stiffness is the same whatever the direction considered (random micro-crack orientation).

The growth of damage is defined following a framework that is common to many non-linear (rate-independent) constitutive relations. A distinction has to be made between an increment of strain during which micro-cracking grows and an increment of strain where it does not. The growth of damage is controlled by the tensile strain:

$$ \tilde{\varepsilon} = \frac{\varepsilon + |\varepsilon|}{2}. \quad (5) $$

The limit of the elastic domain is defined by:

$$ f(\varepsilon) = \tilde{\varepsilon} - Y = 0 \quad \text{with } Y = \max(Y_{D0}, \tilde{\varepsilon}) \quad \text{over the loading history} \quad (6) $$

Y is the history variable that keep tracks of the past applied loads that possibly produced some damage and $Y_{D0}$ is a damage threshold. The growth of damage is defined according to the following conditions:

$$ \text{If } f(\varepsilon) < 0 \text{ or if } f(\varepsilon) = 0 \text{ and } \dot{\varepsilon} \le 0 \text{ then } \dot{D} = 0 \quad (7a) $$

$$ \text{If } f(\varepsilon) = 0 \text{ and } \dot{\varepsilon} > 0 \text{ then } D = 1 - \frac{1}{\exp(B \cdot (\tilde{\varepsilon} - Y_{D0}))}. \quad (7b) $$

Equation (7a) corresponds to conditions where damage should not grow and (7b) corresponds to the case of growing damage due to an excess of tensile strain. Time derivatives in these equations represent the sequence of loading because the material response is rate-independent. Equations (5)–(7) set a framework for the non-linear response of the material due to micro-cracking that is very similar to elasto-plasticity. The difference is that the non-linear response is controlled here by strain-based quantities, whereas it is controlled by stress-based quantities in classical elasto-plasticity.

Figure 6 shows a typical tensile response obtained with the model. The material parameters (E, $Y_{D0}$, and B) are typically those of concrete. Under monotonically increasing strain and following the elastic regime and the peak stress which corresponds to the onset of damage, the response exhibits strain softening, that is a decreasing stress with increasing strain.

It is important here to consider the case where the strain is not monotonically increasing. Say that the strain has first increased so as to reach point A on the material response in Figure 6. Then, the strain rate is reversed and unloading occurs. Because of the conditions in (7a), damage stays constant. According to (4), the stress decreases linearly towards the origin of the curve. The slope is $E \cdot (1 - D)$ instead of the slope observed at the onset of loading which is E. The comparison of these two slopes is the evidence of damage in experimental results. At point A, a positive incremental strain was referred to as “loading” and (7b) applies in order to calculate the increment of damage. At point A, a negative incremental strain is referred to as “unloading” and damage is constant according to (7a). As we will see next, both incremental strains may yield the same incremental stress. This is the source of the localisation of strain and damage.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Stress–strain relation in tension ($E$ = 37.7 GPa, $Y_{D0}$ = $10^{-4}$, $B$ = 14,000).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Tensile bar.</figcaption></figure>

**3. Strain and damage localisation**

For many years, strain softening has been considered as something impossible in continuum mechanics, although multiscale analyses of progressive cracking were clearly pointing out that it was an expected property. The tangent modulus of the material in the softening regime $E_t$ (the ratio of the increment of stress divided by the increment of strain) is negative. Infinitesimal waves cannot propagate because their celerity, which involves the square root of the tangent modulus, is not real. This violates Hadamard's condition [14]. The consequence is that the incremental equations of equilibrium, along with the classical boundary conditions in statics form an ill-posed problem from the mathematical point of view (the governing differential equations are hyperbolic instead of being elliptic). This can be easily illustrated in the example of a softening bar subjected to tension [15].

**3.1. Localisation in a tensile bar**

Consider a tensile bar, clamped at one end and to which a monotonically increasing displacement $U$ is imposed at the other end (Figure 7).

The bar is initially at equilibrium, in a state of homogeneous strain $\epsilon_0$ and the corresponding displacement is $U_0 = L \cdot \epsilon_0$. The stress is derived directly from the constitutive relation, we assume that it belongs to the softening branch e.g. at point A in Figure 6. Since the strain is constant, the stress is also constant and therefore, the equilibrium conditions are verified.

Consider now that an incremental displacement $\delta U$ is applied from this initial state of equilibrium. A possible solution is that the strain increases and remains constant over the bar. There are, however, other possible solutions: assume that the profile of the incremental strain is set according to Figure 8: a part of the bar of length $h$ undergoes loading with a positive incremental
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Incremental perturbation due to an incremental displacement $\delta U$: strain profile (left) and constitutive law (right). The initial state corresponds to a displacement $U_0 = L \cdot \varepsilon_0$ where the strain $\varepsilon_0$ is uniform over the bar and lies in the strain softening regime.</figcaption></figure>

strain $\delta\varepsilon_t$, and the rest undergoes unloading with a negative incremental strain $\delta\varepsilon_u$. We say that the strain localises over the bar length $h$. The unloading and tangent moduli are fixed by the constitutive relation (considered here bilinear for simplicity, Figure 8). Equilibrium requires that the incremental stress is the same in both part of the bar:

$$E_t \cdot \delta\varepsilon_t = E_u \cdot \delta\varepsilon_u. \qquad (8)$$

Meanwhile, the compatibility between the imposed incremental displacement $\delta U$ and the incremental strain requires that:

$$\delta U = (L - h)\delta\varepsilon_u + h\delta\varepsilon_t. \qquad (9)$$

For each possible value of $h$, and for a given incremental displacement $\delta U$, it is possible to find a suitable set of incremental strains ($\delta\varepsilon_t, \delta\varepsilon_u$) that satisfies the two equations of the problem (Equations (8), (9)). In addition, the relation between the incremental force applied to the bar $\delta F$ and the incremental displacement $\delta U$ is:

$$\delta F = k \cdot \delta U \quad \text{with} \quad k = \frac{A \cdot E_u}{\left[ (L - h) + \frac{h E_u}{E_t} \right]} \qquad (10)$$

where $A$ is the area of the cross section of the bar. When $h = L$, the entire bar undergoes a positive incremental strain and the tangent stiffness $k = AE_t/L$ is negative. When $h$ tends to zero, the tangent stiffness becomes $k = AE_u/L$ and it is positive. For any value of $h$ in between, the tangent stiffness is given by (10). This simple problem shows that, there is an infinite number of solutions in response to an applied incremental displacement at the bar end. This is the reason why the mechanical problem is ill-posed. It is due to strain softening otherwise equilibrium cannot be satisfied, and it may occur right after the peak load is reached.

But then, which solution the bar is going to follow? The solution is not provided by the equilibrium equation, but rather by stability considerations. The path followed by the structure is the one that maximises dissipation of energy locally. It corresponds to $h = 0$, which means that positive incremental strain exists in a segment of the bar of vanishing length only, the rest of the bar is unloading. Accordingly, the response of the bar is the one shown in Figure 9.

Under increasing displacement, the bar is first elastic until it reaches the threshold of damage growth (Equation (6) is satisfied). According to Figure 6, this threshold is the peak load. It is at the peak load that strain has the first opportunity to localise, meaning that because the tangent modulus upon loading becomes negative, a distribution of the incremental strain illustrated in Figure 8 becomes admissible. The solution that maximises the local dissipation at this stage of the loading history corresponds to $h = 0$ because it yields the largest possible positive increment of
<figure><img src="image_12.png" /><figcaption>Figure 9. Response of a tensile softening bar.</figcaption></figure>

strain locally (infinite when $h = 0$) and the corresponding stiffness of the bar is provided by (10): $k = AE_u/L$. It is *exactly* the same stiffness as the one observed during loading, but this time, the increment of displacement is negative (9) and *the response of the bar will be exactly the same as during loading, followed in a reversed way* as the strain at the material point where damage grows increases towards infinity (and the force decreases to zero). When, complete failure is reached and the displacement may increase again, but with a vanishing applied force on the bar.

It follows that the energy dissipated during failure, which is the area under the response of the bar in this figure is zero—failure occurs without energy dissipation which is physically unrealistic. It should be stressed here that under monotonically increasing displacement, the force will suddenly jump from the peak load to zero, following the vertical dashed line in Figure 9. Then, the area under the force displacement curve is not seen as being zero. One should be cautious, however: at peak load, the bar is no longer stable meaning that a small perturbation of the control parameter (positive perturbation of the applied displacement) does not yield a small perturbation of the response (rather a large, finite size, jump). Upon the application of a positive incremental displacement, the bar jumps into a dynamic regime and the mechanical work is converted into kinematic energy. It is not converted into energy dissipated during fracture as the response of the bar snaps back during this dynamic process and still follows the response illustrated by the arrows in Figure 9.

## 3.2. Nonlocal damage model

In the above problem, the important feature is that upon strain localisation, the size of the localisation band ($h$ in the problem) is arbitrary. Therefore, a direct way to avoid the occurrence of an infinite number of solutions is to set the size of the segment of the bar where strain localises. This is achieved indirectly by the nonlocal damage model. The tensile strain (5) is substituted by its spatial average $\bar{\epsilon}$ in the equations governing the growth of damage [9]:

$$ \bar{\epsilon} = \frac{\int \phi(x-s)\tilde{\epsilon}(s) \, ds}{\int \phi(x-s) \, ds} \qquad (11) $$

where $\phi$ is a weight function, for instance the Gaussian function:

$$ \phi(x-s) = \exp\left[-\left(\frac{2|x-s|}{l}\right)^2\right] \qquad (12) $$

$l$ is the internal length of the model, a new model parameter. It is this internal length which sets the size of the zone in which strain localises, or more generally, the width of the FPZ observed in Figure 1.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Response of the bar in tension with the nonlocal damage model: force displacement response (left), damage profiles (right). The internal length is set large enough in order to better illustrate the growth of damage.</figcaption></figure>

This model is nonlocal (in the sense that a pointwise quantity is a spatial average) with local strains. Only the variable that controls damage is nonlocal, which is easier to handle compared to a fully nonlocal constitutive relation. With the model parameters given in Figure 6, for a bar of length 1 m, and an internal length $l = 0.32$ m, the response of the bar is shown in Figure 10 along with the evolution of the profiles of damage. Of course, inserting an internal length in a constitutive relation raises several questions, among them the issue of the experimental determination of this parameter. This length may be approached by considering a tensile test in which strain localisation occurs and comparing the response with an experiment in which strains are constrained to remain uniform (e.g. by gluing elastic fibers on the surface of the specimen). This comparison has been performed by Bažant and Pijaudier-Cabot [16] on concrete and provided an estimate of the internal length ($l \approx 3d$) where $d$ is the maximum size of the heterogeneity.

Another technique consists in considering size effect on fracture (for a detailed discussion of structural size effect see [17]). Upon strain localisation and fracture, the elastic energy stored in the specimen flows in the FPZ where it is dissipated by micro-cracking. Because the size of the fracture process zone is fixed, its ratio with the size of the specimen is changing when the later changes. The amount of flowing energy depends on the size of the structure and it yields size effect. Macroscopically, this induces a dependence of the fracture energy on the crack length [18] and also a dependence of the “nominal strength” of the structure on its size that can be illustrated by considering experiments on geometrically similar specimens.

Figure 11 shows experiments on three-point bending notched specimens of four different sizes similar to those in Figure 1. They are geometrically similar. Going from one size to another, dimensions are multiplied by two, except the thickness which is kept constant.

Because the beams are geometrically similar, they all have the same elastic stiffness. When the size of the beam is multiplied by two, and assuming that the peak load is reached when the nominal tensile strength is reached at the tip of the notch, the maximum load should, according to elasticity and strength of materials, be multiplied by two. This is not the case experimentally, as observed in Figure 11. Going from one size to the next one, the peak load is changed by a factor which is always less than 2. Therefore, the nominal tensile strength of the material decreases as the size of the beam increases. The nonlocal damage model captures such a size effect, and the inverse analysis of these experiments (e.g. using the finite element method) yields the internal length [20].

This original version of the nonlocal model suffers from several shortcomings. Among these, the model cannot capture complete failure with a constant width of the FPZ. The FPZ enlarges indefinitely upon failure and this is due to the averaging across the FPZ: two material points
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Geometrically similar bending specimens (left) and their mechanical response (right), after [19]. *n* = 4 is the smallest size and average of experimental response are shown.</figcaption></figure>

located on both sides of a crack surface do interact; e.g. enter in the spatial averaging centered on one side of the crack, which is not realistic. Enhancements of the nonlocal model aimed at solving this problem have been proposed in the literature, either by indexing the weight function on the stress state [21], or by defining the weight function from interactions between defects [22]. Also, the treatment of averaging near the boundary of a solid is somewhat arbitrary as the weight function is chopped off. This induces an improper description of crack initiation near the boundary of the solid [23] and also a bias in the description of size effect when unnotched beams are considered [24]. Again, the initial nonlocal model can be enhanced with a variation of the internal length nearby boundaries of the structure in order to solve this issue. These enhancements may be, however, regarded as second order ingredients in most practical cases, the important feature being the proper description of progressive damage, its localisation in order to form a FPZ, which in many instances corresponds to the maximum load that the structure can carry and is the important quantity for design.

### 3.3. Other localisation limiters

The nonlocal damage model belongs to a class of constitutive model called “localisation limiters”. An internal length is inserted in order to set a non-zero size of the FPZ, and thus to “limit” localisation. Since the early nonlocal damage proposal, there has been many other proposals inspired from the same point of departure. Here, we shall restrict ourselves to damage-based constitutive models, as many others are based on plasticity (see [25]), following the same principles.

The most popular one is the gradient damage model [26]. Instead of (11), the nonlocal strain is the solution of the following partial differential equation:

$$ \bar{\epsilon} - c\nabla^2\bar{\epsilon} = \bar{\epsilon} \qquad (13) $$

where $\nabla^2 x$ is the Laplacian of $x$, $c$ is a quantity with dimension $m^2$, the square of an internal length. In fact, this equation is exactly the same as (11), if a specific exponential weight function is used. Later on, Pijaudier-Cabot and Burlion [27] derived a more general gradient damage model, extending in the nonlinear regime the theory of materials with voids due to [28]. This theory, initially developed for porous materials such as bones, embeds a general variational framework [29] which is very similar to the gradient damage model proposed by Fremond and Nedjar [30] derived from the principle of virtual power. In the work of Cowin and Nunziatto,
however, the equation governing the growth of damage is not introduced a priori, it is intended to capture the effects of dilation centers, i.e. forces that produce a growth of the voids or cracks in the material. The equation governing the growth of damage reads [27]:

$$D - k_l \nabla^2 D - \mathcal{F} = 0 \qquad (14)$$

where $k_l$ is the square of an internal length and $\mathcal{F}$ is a function that defines the growth of damage, typically a function of the local strains. It is one of the governing equation of the boundary value problem, along with the equilibrium equation. Boundary conditions on damage must be enforced, raising the issue of damage initiation near a free boundary, same as in the integral model. These are set arbitrarily in the absence of any theoretical or physical motivation.

About two decades ago, a new kind of failure model became popular: the phase-field models. The phase-field theory was designed originally to model solidification, often a first order phase transition, and proved being very effective for this purpose [31]. The liquid and a solid state of the material are treated as two materials with the appropriate conditions that transform the former into the latter (and vice versa). The extension to failure is appealing, each state of damage describing the material as a specific “phase” (of micro-cracking), the equation controlling the growth of damage providing the conditions for “phase changes”.

Phase-field model have been first implemented in dynamic brittle fracture problems [32]. Later on Miehe and co-workers were among those who developed the theory for failure analyses (see e.g. [33]), and benefited from the regularized setting of the variational theory of linear elastic fracture (see e.g. [34]). The extension to cohesive fracture, or progressive cracking at the crack tip, is not trivial however, mainly because the crack opening needs to be computed. Verhoosel and de Borst [35] illustrated this difficulty. They also came out with the conclusion that for cohesive fracture, phase-field modelling of fracture and gradient damage model are very similar. The phase-field approach provides a better description of complete failure that can be achieved with the later by letting the internal length tend to zero upon complete failure. Recently also, the phase-field and continuum damage approaches have been unified into a single model which is capable to capture size and boundary effects for concrete [36].

# 4. Coupled damage and permeability evolutions

Up to now, we have been interested in the description of the mechanical response of a material undergoing progressive and localised micro-cracking due to mechanical loads. In practice, most geomaterials are porous and contain one or several fluids. These fluids are transported in the pore network of the material and there is a strong interest at looking at the interaction between the state of damage and cracking and the consecutive evolution of these transport properties. Typically, one may think to the stimulation of geological reservoirs that contain hydrocarbon [37], to the evolution of the tightness of nuclear vessels [38] or to the durability of concrete structures as the fluid present in the pore may carry aggressive species.

## 4.1. Evolving permeability of a damaged material

At first, let us consider again the lattice discussed in Section 2. In order to represent the influence of material damage on the permeability, it is reasonable to assume that when a bond fails in the mechanical lattice, it opens a larger path for fluid flow in the perpendicular direction. This basic scheme is depicted in Figure 12 [39]. A hydraulic lattice is attached to the mechanical one following this principle, with bonds that are perpendicular to the mechanical lattice (Figure 13).

The fluid flow at the micro-scale level (bond level) is described by Darcy's equation:

$$\vec{q} = -k \overrightarrow{\text{grad}(p)} \quad \text{and} \quad k = \frac{\kappa}{\mu} \qquad (15)$$
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Basic scheme for the coupled hydro-mechanical lattice analysis. When a mechanical bond fails, the permeability of a perpendicular bond increases suddenly.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 13.** The mechanical and the hydraulic lattices (after [39]).</figcaption></figure>

where $\bar{q}$ is the flow rate, $p$ is the pressure, and $k$ is the permeability of the material. This permeability is a function of the dynamic viscosity of the fluid $\mu$, and an intrinsic permeability denoted as $\kappa$ which depends on the pore network only (some tortuosity factor is also added in this formula traditionally). We are going to assume that when a bond fails, the permeability in the perpendicular direction is increased by an amplification factor which is very large, typically of the order of $10^6$.

Periodic boundary conditions are also applied to the hydraulic lattice. A constant drop of pressure equal to 1 is applied at the vertical boundaries of the hydraulic lattice. At each step of damage in the mechanical lattice, the flow rate is computed in the dual lattice and the various moments of its distribution are computed too.

It is worthwhile to point out that we analyse here the permeability of the unloaded material. There is no effect of the applied stress on the permeability in the reversible regime and it cannot be expected to be observed in the present computations. Following the same technique as in the mechanical problem, we look at the moments of the distribution of the hydraulic flow rate, and we try to find for which moment of the distribution of the local stress the plots collapse on the same curve for different sizes.

The second order moment in the mechanical problem is a natural candidate since it is this one which describes the evolution of mechanical damage. Figure 14 shows the evolution of the lattice permeability (second order moment of the flow rate distribution) as a function of the evolution of stiffness in the mechanical lattice. For all the sizes considered, the plots collapse onto the same curve for stiffness variations ranging from 0 to 0.2 approximately, before the peak load. For
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Permeability versus Young Modulus reduction according to the lattice analysis (after [39]).</figcaption></figure>

large variations of the stiffness, macro-cracking has occurred and the lattice size starts to play an important role because it controls the spacing between the macro-cracks. Results cannot be interpreted in this case.

It follows that for a correct description of the coupled effects between the growth of micro-cracking and the growth of permeability in a porous material, the permeability should depend on the damage variable. Typically, a power law is used in the literature [38]:

$$ \kappa \approx \kappa_0 \cdot 10^{A_D D} \qquad (16) $$

where $\kappa_0$ is the initial intrinsic permeability for the undamaged material and $A_D$ is a model parameter.

## 4.2. Permeability of a cracked material

At complete rupture, one or several macro-cracks are expected and fluid flow will be governed by these cracks. Poiseuille's law should be considered instead of Darcy's law. For a fluid flowing between two parallel plates, we have:

$$ \vec{q} = -k_p \overrightarrow{\text{grad}(p)} \quad \text{where } k_p = \frac{[u]^2}{12\mu} \qquad (17) $$

$k_p$ is the permeability, $[u]$ is the crack opening, i.e. the distance between the plates, and $\mu$ is the dynamic viscosity of the fluid. Some factor taking into account the crack roughness may be added too and provide a decrease of the permeability compared to the case of a perfectly flat crack. As opposed to (16), this equation does not define the permeability at a material point but rather it defines the fluid flow inside an opened crack. Typically, this quantity is several orders of magnitude larger than the fluid flow in the porous material governed by Darcy's law for tight materials such as concrete and rocks.

On one hand and prior to the localisation of damage, the fluid flow is described within a continuum description. On the other hand, after damage has localised, fluid flow is governed by the crack opening, and it is typically a “discrete” description in the sense that the crack opening is a displacement discontinuity in the material. The transition between the two regimes could be seen as the percolation of micro-cracks which connect themselves and form a macro-crack.

In order to capture the two regimes of the hydraulic problem in a single framework, similarly to the mechanical model which captures both distributed and localised micro-cracking, it is
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Maps of damage and fluid pressure in a hydraulic fracturing process in the presence of a vertical pre-existing joint (after [41]).</figcaption></figure>

necessary to calculate crack opening from the mechanical continuum model. This is not a trivial issue because at first, we assumed that the quantities describing the deformation in the material should be continuous and therefore, a displacement discontinuity is not allowed. Now there is a need to transform back from a continuum into a *discontinuum*. There are several possibilities for calculating this crack opening. The simplest one is to define the law that governs the growth of damage as a function of a crack opening displacement which is directly calculated as the product of the strain by a fixed length $h_c$. This length is related to the internal length in the nonlocal continuum, in fact it is the thickness of the FPZ which is proportional to the internal length in the nonlocal damage model. What is being done here is to collapse the FPZ onto a line, assuming that the opening of the crack is the relative displacement at opposite boundaries of the FPZ. This model is very close to the crack band model proposed by Bažant and Oh [40], and further quite similar to the cohesive crack model for concrete [7]. It is a very effective approach which has been used e.g. for modelling hydraulic fracturing [41].

Figure 15 shows the calculation of the propagation of a crack in a hydraulic fracture process and its interaction with an existing joint in the rock mass perpendicular to the crack propagation. The maps of the fluid pressure and damage have been superimposed. We can see that the joint arrests the initial crack and that new cracks develop at the tips of the joint in this calculation.

This simple approach assumes that the width of the FPZ is known in advance, that it is constant, and therefore that it does not depend on the applied loads or on the load history that yielded damage. If one wants to avoid this assumption, a possibility is to extract a crack opening displacement from the results of computations with the nonlocal damage model. The method, proposed originally by Dufour et al. [42] consists in comparing the nonlocal strain field $\bar{\epsilon}$ across the FPZ, calculated according to (11), to the analytical nonlocal strain field that would be obtained from a displacement discontinuity across the crack path. The displacement discontinuity is a heavyside function and therefore the nonlocal strain field reduces to the displacement jump times the (normalised) weight function. The comparison yields (1) the displacement jump that fits best the analytical field with the nonlocal strain field, and (2) an error which characterizes the quality of the fit, that is how the calculated profile is close to the analytical one. From this error, one may decide whether or not a macro-crack that connects the micro-cracks has formed (by comparing this error to a given threshold defined a priori), and in this case what is the crack opening displacement. Then, the crack permeability (17) can be calculated and within a contin-
uum setting, it is smeared over the FPZ so that the fluid flow in the FPZ is the same than inside the macro-crack [43].

A side result to this comparison addresses the following question: does the description of nonlocal (continuum) damage model converge to the description of a discrete crack when the material is totally damaged? This is not the case if the internal length in the nonlocal model is set constant [44, 45]. It should tend to zero upon full damage, or the weight function in the nonlocal averaging process should better capture the interactions between micro-cracks during their coalescence towards a macro-crack [22].

**5. Concluding remarks**

Failure due to progressive micro-cracking occurs in geomaterials such as rocks or concrete subjected to tensile loads. The material response is quasi-brittle. Micro-cracks propagate and are arrested or their path deviated due to the material heterogeneity, and this process cannot be captured by linear elastic fracture mechanics. Within a continuum mechanics framework, the material response may be described by continuum damage mechanics. Due to progressive micro-cracking, the material stiffness decreases. It results in strain softening.

In order to avoid spurious strain and damage localisation, it is necessary to introduce an internal length in the constitutive model so that failure occurs with a non-zero energy dissipation. This is the purpose on nonlocal continua with local strains, the nonlocal damage model being the first of this kind historically. A feature of such a model is size effect on the nominal strength, which is observed in laboratory tests and may have a great importance in civil engineering and geomechanics because experiments are always conducted on small specimens compared to the real structures. Since the apparent strength decreases when the size of the specimen increases, strength-based design (which is the usual engineering practice) ought to account for this size effect.

There are today a wide variety of failure models based on damage mechanics, phase field approaches or cohesive cracking that incorporate such an internal length. It has become a standard, although there is still a debate on whether this internal length should be viewed as a constant or not.

This consistent continuum approach to failure opened also the path to solving coupled problems. Here, we have considered hydro-mechanical coupled effects. Micromechanics shows that the permeability of the material ought to be indexed to the damage growth in the regime of distributed cracking. Upon macro-cracking, there is a change of regime and it is the crack opening that controls the fluid flow in the cracked material. Both regimes may be captured with a continuum damage approach, however.

To conclude, there is a wide variety of coupled problems involving fracture and damage coupled with other chemo-physical mechanisms. To name a few, let us mention calcium leaching yielding to the dissolution of the material (important for the long-term safety of waste containments), and salt or ice crystallisation in the pore structure of geomaterials and alcali-silica reactions in concrete which induces internal cracking. These problems, gathered today under the name of “durability mechanics” are intrinsically multi-disciplinary, and are not solely concerned with geomaterials. Biological materials and materials for energy storage are also of concern because they are porous materials which exchange fluids with their environment.

**Acknowledgements**

The author would like to thank Z. P. Bažant, N. Burlion, G. Chatzigeorgiou, M. Choinska, A. Delaplace, F. Dufour, D. Grégoire, A. Huerta, A. Khelidj, A. Krayani, V. Lefort, L. Rojas Solano, and
S. Roux for their great help and inspiration which lead to the results summarised in this paper. This work has been partially supported by the Investissement d'Avenir French programme (ANR-16-IDEX-0002) under the framework of the E2S UPPA hub Newpores.

References

[1] G. R. Irwin, "Analysis of stresses and strains near the end of a crack traversing a plate", *Trans. ASME, J. Appl. Mech.* **24** (1957), p. 361-364.

[2] J. F. Knott, *Fundamentals of Fracture Mechanics*, Butterworth and Co., Delft, the Netherlands, 1973.

[3] A. A. Griffith, "The theory of rupture", in *Proc. 1st Int. Conf. of Applied Mech.*, 1924, p. 55-63.

[4] K. Haidar, G. Pijaudier-Cabot, J. F. Dube, A. Loukili, "Correlation between internal length, fracture process zone and size effect in mortar and model materials", *Mater. Struct.* **38** (2005), p. 201-210.

[5] G. I. Barenblatt, "The mathematical theory of equilibrium cracks in brittle fracture", *Adv. Appl. Mech.* **7** (1962), p. 55-129.

[6] D. Dugdale, "Yielding of steel sheets containing slits", *J. Mech. Phys. Solids* **8** (1960), p. 100-108.

[7] A. Hillerborg, M. Modeer, P. E. Petersson, "Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements", *Cement Concr. Res.* **6** (1976), p. 773-782.

[8] J. Mazars, G. Pijaudier-Cabot, "Continuum damage theory—application to concrete", *J. Eng. Mech ASCE 115* (1989), p. 345-365.

[9] G. Pijaudier-Cabot, Z. P. Bažant, "Nonlocal damage theory", *J. Eng. Mech. ASCE 113* (1987), p. 1512-1533.

[10] Z. P. Bažant, G. Pijaudier-Cabot, "Nonlocal continuum damage, localization instability and convergence", *Trans. ASME, J. Appl. Mech.* **55** (1988), p. 287-294.

[11] J. Lemaitre, J. L. Chaboche, *Mécanique des matériaux solides*, Dunod, Paris, France, 1985.

[12] D. Krajcinovic, J. G. M. Van Mier, *Damage and Fracture of Disordered Materials*, CISM Courses and Lectures No. 410, Springer Verlag, Wien, Austria, 2000.

[13] A. Delaplace, G. Pijaudier-Cabot, S. Roux, "Progressive damage in discrete models and consequences on continuum modeling", *J. Mech. Phys. Solids* **44** (1996), p. 99-136.

[14] J. Hadamard, *Leçons sur la propagation des ondes et les équations de l'hydrodynamique*, Hermann, Paris, France, 1903.

[15] Z. P. Bažant, "Instability, ductility and size effect in strain-softening concrete", *J. Eng. Mech. ASCE 102* (1976), p. 331-344.

[16] Z. P. Bažant, G. Pijaudier-Cabot, "Measurement of the characteristic length of nonlocal continuum", *J. Eng. Mech. ASCE 115* (1989), p. 755-767.

[17] Z. P. Bažant, J. Planas, *Fracture and Size Effect in Concrete and Other Quasi-brittle Materials*, CRC Press, London, UK, 1998.

[18] Z. P. Bažant, P. A. Pfeiffer, "Determination of fracture energy from size effect and brittleness number", *ACI Mater. J.* (1987), p. 463-480.

[19] D. Grégoire, L. Rojas-Solano, G. Pijaudier-Cabot, "Failure and size effect for notched and unnotched concrete beams", *Int. J. Numer. Anal. Methods Geomech.* **37** (2013), p. 1434-1452.

[20] C. Le Bellego, J. F. Dube, G. Pijaudier-Cabot, B. Gérard, "Calibration of nonlocal damage model from size effect tests", *Eur. J. Mech. A 22* (2003), p. 33-46.

[21] C. Giry, F. Dufour, J. Mazars, "Stress-based nonlocal damage model", *Int. J. Solids Struct.* **48** (2011), p. 3431-3443.

[22] L. Rojas Solano, D. Grégoire, G. Pijaudier-Cabot, "Interaction based nonlocal damage model for failure in quasi-brittle materials", *Mech. Res. Commun.* **54** (2013), p. 56-62.

[23] A. Simone, H. Askes, L. J. Sluys, "Incorrect initiation and propagation of failure in non-local and gradient-enhanced media", *Int. J. Solids Struct.* **41** (2004), p. 351-363.

[24] A. Krayani, G. Pijaudier-Cabot, F. Dufour, "Boundary effect on weight function in nonlocal damage model", *Eng. Fract. Mech.* **76** (2009), p. 2217-2231.

[25] Z. P. Bažant, M. Jirasek, "Nonlocal integral formulations of plasticity and damage: survey of recent progress", *J. Eng. Mech. ASCE 128* (2002), p. 1119-1149.

[26] R. H. J. Peerlings, R. de Borst, W. A. M. Brekelmans, J. H. P. de Vree, "Gradient enhanced damage for quasibrITTLE materials", *Int. J. Numer. Methods Eng.* **39** (1996), p. 3391-3403.

[27] G. Pijaudier-Cabot, N. Burlion, "Damage and localisation in elastic materials with voids", *Mech. Cohesive Frict. Mater.* **1** (1996), p. 129-144.

[28] S. C. Cowin, J. W. Nunziato, "Linear elastic materials with voids", *J. Elast.* **13** (1983), p. 125-147.

[29] S. C. Cowin, M. A. Goodman, "A variational principle for granular materials", *Z. Angew. Math. Mech.* **56** (1976), p. 281-286.
[30] M. Fremond, B. Nedjar, "Endommagement et principe des puissances virtuelles", C. R. Acad. Sci., Paris II (1993), p. 857-864.

[31] N. Provatas, K. Elder, *Phase-field Methods in Material Science and Engineering*, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, 2010.

[32] A. Karma, D. A. Kessler, H. Levine, "Phase-field model for mode III dynamic fracture", Phys. Rev. Lett. 87 (2001), article no. 045501.

[33] C. Miehe, M. Hofacker, F. Welshinger, "A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator split", Comput. Methods Appl. Mech. Eng. 199 (2010), p. 2765-2778.

[34] B. Bourdin, G. Francfort, J. J. Marigo, *The Variational Approach to Fracture*, Springer, New York, USA, 2008.

[35] C. V. Verhoosel, R. de Borst, "A phase-field model for cohesive fracture", Int. J. Numer. Methods Eng. 96 (2013), p. 43-62.

[36] D. C. Feng, J. Y. Wu, "Phase-field regularised cohesive zone model and size effect of concrete", Eng. Fract. Mech. 197 (2018), p. 66-79.

[37] G. Pijaudier-Cabot, C. La Borderie, T. Rees, W. Chen, O. Maurel, F. Rey-Betbeder, A. de Ferron, *Electrohydraulic Fracturing of Rocks*, ISTE-Wiley, London, UK, 2016.

[38] L. Jason, G. Pijaudier-Cabot, S. Ghavamian, A. Huerta, "Hydraulic behaviour of a representative structural volume for containment buildings", Nucl. Eng. Des. 237 (2007), p. 1259-1274.

[39] G. Chatzigeorgiou, V. Picandet, A. Khelidj, G. Pijaudier-Cabot, "Coupling between progressive damage and permeability of concrete: analysis with a discrete model", Int. J. Numer. Anal. Methods Geomech. 29 (2005), p. 1005-1018.

[40] Z. P. Bažant, B. H. Oh, "Crack band theory for fracture of concrete", Mater. Struct. 16 (1983), p. 155-177.

[41] V. Lefort, O. Nouailletas, D. Grégoire, G. Pijaudier-Cabot, "Lattice modelling of hydraulic fracture: theoretical validation and interactions with cohesive joints", Eng. Fract. Mech. 235 (2020), article no. 107178.

[42] F. Dufour, G. Pijaudier-Cabot, M. Choinska, A. Huerta, "Extraction of crack opening from a continuous approach using regularised damage models", Comput. Concr. 5 (2008), p. 375-388.

[43] G. Pijaudier-Cabot, F. Dufour, M. Choinska, "Permeability due to the increase of damage in concrete: from diffuse to localised damage distributions", J. Eng. Mech. ASCE 135 (2009), p. 1022-1028.

[44] D. Grégoire, L. Rojas Solano, G. Pijaudier-Cabot, "Continuum to discrete transition in nonlocal damage models", Int. J. Multiscale Comp. Eng. 10 (2012), p. 567-580.

[45] G. Pijaudier-Cabot, D. Grégoire, "A review of nonlocal continuum damage: modelling of failure?", Netw. Heterog. Media 9 (2014), p. 575-597.