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

HAL
open science

# New continuous strain-based description of concrete's damage-permeability coupling

David Bouhjiti, Mohamad Ezzedine El Dandachy, Frédéric Dufour, Stefano Dal Pont, Matthieu Briffaut, Julien Baroth, Benoît Masson

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

## To cite this version:

David Bouhjiti, Mohamad Ezzedine El Dandachy, Frédéric Dufour, Stefano Dal Pont, Matthieu Briffaut, et al.. New continuous strain-based description of concrete's damage-permeability coupling. International Journal for Numerical and Analytical Methods in Geomechanics, 2018, 42 (14), pp.1671-1697. 10.1002/nag.2808. hal-02045996

HAL Id: hal-02045996

https://hal.univ-grenoble-alpes.fr/hal-02045996

Submitted on 14 Mar 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# New continuous strain-based description of concrete's damage-permeability coupling

D. E.-M. Bouhjiti<fnref n="1" />,<fnref n="2" /> ID | M. Ezzedine El Dandachy<fnref n="5" /> | F. Dufour<fnref n="1" />,<fnref n="3" /> | S. Dal Pont<fnref n="1" /> ID |
M. Briffaut<fnref n="1" /> | J. Baroth<fnref n="1" /> | B. Masson<fnref n="4" />

## Correspondence

D. E.-M. Bouhjiti, Grenoble INP Partnership Foundation—Industrial Chair PERENITI, Grenoble, France.

Email: david.bouhjiti@3sr-grenoble.fr

F. Dufour, Univ. Grenoble Alpes, CNRS, G INP, 3SR, Grenoble, France.

Email: frederic.dufour@3sr-grenoble.fr

## Funding information

France's National Research Agency, Grant/Award Number: ANR-11-RSNR-012; Industrial Chair PERENITI (EDF-SEPTEN/DTG/CIH)

## Summary

This contribution aims at improving the finite element modelling of concrete's damage-permeability coupling using continuous and elastic-based damage approaches. With this aim in view, a review of existing local and nonlocal damage formulations is performed showing how they differ in terms of damage estimation, crack propagation, and spatial discretization restrictions (when the energy-based Hillerborg criterion is applied). Another critical review of damage-based permeability laws is achieved to pinpoint their limitations. Indeed, the use of a bounded and strictly increasing damage variable does not allow (1) an accurate description of permeability in the presence of macrocracks that are physically unbounded and (2) a proper description of the permeability decrease as those macrocracks are closed. Accordingly, a new strain-based permeability law is suggested. Its analytical validation is based on the BIPEDE tensile test for one loading-unloading cycle. This new formulation intends to enhance the predictiveness of damage-permeability coupling models leading to a better assessment of large structure's durability and serviceability.

## KEYWORDS

concrete, cracking, finite elements, regularization, strain-based permeability model

# 1 | INTRODUCTION

The tightness of strategic and large structures such as dams and nuclear containment buildings is due to the low permeability of concrete and, mostly, to the continuous control and monitoring of their cracking state throughout their operational lifespan. With that regard, the ability to predict the evolution of damage and cracking states in time, in addition to the resultant permeability, are key issues for the durability and serviceability assessment process.

As it has been demonstrated by numerous experimental works, the permeability of concrete is strongly related to its damage state:

1. Darcy's mode: When damage is low and diffuse, the increase of concrete's permeability is due to the deployment of microcracks within its volume. This leads to a generalized increase of the porous network's connectivity. In the absence of strain localization, concrete is seen as a damage-dependent porous media that justifies the hypothesis of a transfer mode governed by Darcy's law.<fnref n="1" />⁻<fnref n="4" />

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, Grenoble INP, Institute of Engineering Univ. Grenoble Alpes, 3SR, F-38000 Grenoble, France</footnote>

<footnote n="2">²Grenoble INP Partnership Foundation—Industrial Chair PERENITI, Grenoble, France</footnote>

<footnote n="3">³Head of the Chair PERENITI (EDF SEPTEN/DTG/CIH), Grenoble, France</footnote>

<footnote n="4">⁴Electricité De France (EDF-SEPTEN), Lyon, France</footnote>

<footnote n="5">⁵Faculty of Engineering, Beirut Arab University, Beirut, Lebanon</footnote>
2. Poiseuille's mode: As the loading increases, those microcracks localize into macrocracks, which are descriptive of a high and spatially localized strain. The associated permeability in that case is due to the fluid's flow through the macrocracks' openings and is commonly described using the so-called Poiseuille's law.<sup>5-7</sup>

Whether the permeability increase is viewed as an increase in concrete's porous connectivity (Darcy's mode) or as a flow through an equivalent opening (Poiseuille's law), the evolution of permeability from a microcracked (and poorly damaged) to a macrocracked state (and highly damaged) seems to be lacking experimental investigation and characterization. On the one hand, there are the complex questions of (1) the equivalence between damage and its associated cracking state and (2) the threshold from which concrete state changes from a microcracked to a macrocracked one. On the other hand, there is the dependence of the obtained results, in terms of damage estimation and leakage rate measurement, on the experimental set-ups and considered numerical models. This makes the generalization of any test results to all loading cases and configurations hard to achieve. Still, attempts have been made to match the 2 extreme modes (low diffuse and highly localized damage states) using damage-pondered logarithmic or linear laws<sup>8,9</sup> or the linear superimposition principle.<sup>6,10,11</sup> As all previous permeability laws are damage based, the obtained permeability fields can largely differ depending on the used matching strategy and considered damage formulation (local or nonlocal). This issue is problematic and needs to be addressed to enhance the damage-permeability coupling.

One can already notice from the used terminology 2 distinct approaches when dealing with concrete's increase of permeability under loadings: a first one discrete and explicitly related to cracks development and a second continuous based on an elastic-damage concept:

1. **Fracture-based crack models:** These models are based on elastic (or plastic) fracture mechanics to explicitly describe discontinuities within the concrete volume. Amongst this class, one can enumerate cohesive crack models,<sup>12</sup> lattice and network models,<sup>13-15</sup> and LSM-based* inclusive models using either nodal enrichment XFEM<sup>†</sup> or elemental enrichment EFEM<sup>‡</sup>.<sup>16</sup> For such approaches, both, the location and opening of the cracks are known and quantifiable, which facilitate the use of Poiseuille's law. However, for the precracked state, a continuous analysis of the domain is required, which raises the question of the transition from continuous (noncracked) to discrete (cracked) modelling. In Grassl and Bolander<sup>15</sup> and Larsson et al,<sup>17</sup> it is suggested to use a continuous strategy until crack initiation before pursuing calculations in the cracked area using a discrete approach. The threshold to switch control from continuous to cohesive modelling can be set arbitrarily based on strain, stress, or energetic quantities<sup>15,18,19</sup> or defined according to an error analysis between strong discontinuity (SD) gradient-based approach (as a valid approximation of cohesive models) and weak discontinuity (WD) integral-based approach (continuous).<sup>20</sup> It is worth mentioning, however, that the SD approach is a particular case of the WD methods using a sharper and singular weighting function (Green's function).<sup>21</sup> This shows that the WD method and the continuous damage approach can be used to describe equivalently the discontinuous nature of cracks.

2. **Continuous damage-based models:** These models are based on damage mechanics, which regards cracking of concrete as a continuous loss of rigidity.<sup>22</sup> In the particular case of strain-based damage criterion, the damage variable (usually denoted *d*) evolves from 0 to 1 and is computed based on the gap between a certain strain damage threshold and the developed equivalent strain. The way this so-called equivalent strain is computed distinguishes local models from gradient-enhanced ones and integral nonlocal-based ones.<sup>20,21</sup> For the first class (local models), the equivalent strain is computed for each integration point and is mesh sensitive if not regularized. For the second type of models (implicit gradient and nonlocal approaches), mesh sensitivity is alleviated as an internal length is introduced to account for a constant and material-dependent microstructural scale (usually related to the aggregates' size<sup>23</sup>) over which energy dissipation occurs. Nonlocal strategies are generally more time consuming compared with local ones and require more refined mesh to preserve their nonlocal nature. This might lead to considerable computational memory and time cost when dealing with large concrete structural volumes. One can also notice the constant bridging with fracture mechanics even though the continuum media are used: Local damage models are usually regularized using the Hillerborg energy-based approach,<sup>24</sup> and the equivalent crack openings<sup>20,25</sup> and paths<sup>26</sup> are postprocessed using adapted algorithms to overcome the incompatibility between the crack's discrete nature and the used continuous approach. On that issue, the dissipated energy in a continuum domain is supposed to fully

<footnote>*Level-set method</footnote>

<footnote>†Extended finite element method</footnote>

<footnote>‡Energy finite element method</footnote>
contribute to crack propagation and is, consequently, converted into an equivalent crack length.²⁷ This hypothesis remains strong since, experimentally, only a part of the dissipated fracture energy contributes to crack propagation ahead of the crack tip; this amount being intrinsic to the fracture process zone (FPZ) development is estimated at 40% of the total fracture energy.²⁸⁻³⁰ The remaining 60% is associated with frictional energy dissipation. Eventually, despite particular shortcomings related to spurious boundary effects for nonlocal models and mesh dependency for local ones, the ability of continuous elasticity-based damage models to describe initiation and propagation of mode I cracks is not altered when proper regularization and identification techniques are used. Moreover, continuous finite element (FE) models remain, within the currently available computational capacities, the most efficient support for large concrete structures' (~10 m) modelling.

In this contribution, only the permeability associated with mode I cracking is of interest. Accordingly, mode II and mode III cracking modes should be addressed differently as the physics of their associated permeability differ. To serve our purpose, continuous damage-based models are retained and discussed according to their validity domains, their evaluation of damage around the FPZ, their evaluation of the crack opening in the localization point, and their effect on the resulting permeability. So this paper is twofold:

1. The first part of this paper recalls the mathematical framework of the used local, SD, and WD nonlocal modelling strategies (using a strain-based damage criterion) prior to permeability estimation. It also demonstrates the limits of existing damage-based permeability laws that lack a state variable and mathematical formwork able to describe accurately the permeability's evolution from a low damage state to a macrocracked one.

2. The second part introduces a new strain-based permeability law, addressing those drawbacks. It is more physically representative of transfer properties during the transition phase from lowly to highly damaged state. It also describes better the experimental results from the BIPEDE tensile test (one loading-unloading cycle) compared with existing damage-based permeability laws.

## 2 | DAMAGE-BASED DESCRIPTION OF DAMAGE-PERMEABILITY COUPLING

In the framework of a staggered approach, the mechanical calculation (ie, damage and crack opening estimation) represents the first step of the damage-permeability coupling. In this section, the damage formulation principles are recalled along with the postprocessing steps required to define an equivalent crack opening in the continuous domain. In terms of the dissipated fracture energy, energetic criteria are also addressed to ensure mesh independency of local models and softening laws shape independency of both local and nonlocal damage formulations.

### 2.1 | Preview of local and nonlocal damage models

#### 2.1.1 | Isotropic damage modelling

In the case of isotropic damage description, the general behaviour law relating the stress $\sigma$ to the strain $\varepsilon$ variable writes

$$ \sigma_{1D} = (1 - d)E: \varepsilon, \qquad (1) $$

where $E$ is the nondamaged stiffness and $d$ is the isotropic damage variable increasing from 0 to 1 depending on the following loading function:

$$ \begin{gathered} f(\overline{\varepsilon_{eq}}, \varepsilon_{d0}) = \overline{\varepsilon_{eq}} - \varepsilon_{d0} \\ \overline{\varepsilon_{eq}} = \|\|\varepsilon_{eq}\|\|_{t}, \end{gathered} \qquad (2) $$

where $\overline{\varepsilon_{eq}}$ is an equivalent strain (see, for instance, the $\mu$-Mazars damage criterion presented in Mazars et al³¹) that can only increase with respect to the second law of thermodynamics starting from an initial value of $\varepsilon_{d0}$, $\varepsilon_{d0}$ being the strain threshold starting from which the loss of stiffness is obtained. $\|\|f(t)\|\|_{t} = \max(f(\tau \le t))$ is a norm that ensures the monotonic increase over time of a time-dependent function $f$.
Several mathematical models exist for the postpeak behaviour $\sigma = f(\varepsilon \ge \varepsilon_{d0})$. Three types are retained herein for illustration: linear (Equation 3—Figure 1A),²¹ bilinear (Equation 4—Figure 1B),²⁸⁻³⁰ and exponential (Equation 5—Figure 1C)²²,³¹ laws.

$$ \text{Linear (LIN): } \sigma = \begin{cases} \frac{E \varepsilon}{\varepsilon - \varepsilon_u} f_t & \text{if:} \varepsilon \le \varepsilon_{d0} \\ \frac{\varepsilon - \varepsilon_u}{\varepsilon_{d0} - \varepsilon_u} & \text{if:} \varepsilon_{d0} \le \varepsilon \le \varepsilon_u \\ 0 & \text{if:} \varepsilon_u \le \varepsilon \end{cases} \quad d = \begin{cases} 1 - \frac{\varepsilon_{d0} \bar{\varepsilon}_{eq} - \varepsilon_u}{\varepsilon_{eq} \varepsilon_{d0} - \varepsilon_u} & \text{if:} \varepsilon_{d0} \le \bar{\varepsilon}_{eq} \le \varepsilon_u \\ 1 & \text{if:} \varepsilon_u \le \bar{\varepsilon}_{eq} \end{cases}, \quad (3) $$

where $\varepsilon_u$ is the failure strain ($d = 1$).

$$ \text{Bi-linear (BIL): } \sigma = \begin{cases} \frac{E \varepsilon}{\varepsilon - \varepsilon_u' - f_t} & \text{if:} \varepsilon \le \varepsilon_{d0} \\ \frac{\varepsilon - \varepsilon_u}{\varepsilon_{d0} - \varepsilon_u} & \text{if:} \varepsilon_{d0} \le \varepsilon \le \varepsilon_b \\ \frac{\varepsilon - \varepsilon_u}{\varepsilon_b - \varepsilon_u} \frac{f_t}{n} & \text{if:} \varepsilon_b \le \varepsilon \le \varepsilon_u \\ 0 & \text{if:} \varepsilon_u \le \varepsilon \end{cases} \quad d = \begin{cases} 1 - \frac{\varepsilon_{d0} \bar{\varepsilon}_{eq} - \varepsilon_u'}{\varepsilon_{eq} \varepsilon_{d0} - \varepsilon_u'} & \text{if:} \varepsilon_{d0} \le \bar{\varepsilon}_{eq} \le \varepsilon_b \\ 1 - \frac{1}{n} \frac{\varepsilon_{d0} \bar{\varepsilon}_{eq} - \varepsilon_u}{\varepsilon_{eq} \varepsilon_b - \varepsilon_u} & \text{if:} \varepsilon_b \le \bar{\varepsilon}_{eq} \le \varepsilon_u \\ 1 & \text{if:} \varepsilon_u \le \bar{\varepsilon}_{eq} \end{cases}, \quad (4) $$

where $\varepsilon_u$ and $\varepsilon_u'$ are the failure strains if only the first and second slopes are considered respectively. $n$ ($n = 3^{30}$) defines the tensile strength at intersection between the 2 slopes at $\varepsilon = \varepsilon_b = \frac{\varepsilon_{d0} + (n-1)\varepsilon_u'}{n}$.

$$ \text{Exponential (EXP): } \sigma = \begin{cases} \frac{E \varepsilon}{E e^{-B_t (\bar{\varepsilon}_{eq} - \varepsilon_{d0})} \varepsilon} & \text{if:} \varepsilon \le \varepsilon_{d0} \\ \frac{0}{1 - e^{-B_t (\bar{\varepsilon}_{eq} - \varepsilon_{d0})}} & \text{if:} \varepsilon \ge \varepsilon_{d0} \end{cases} \quad d = \begin{cases} 0 & \text{if:} \varepsilon \le \varepsilon_{d0} \\ 1 - e^{-B_t (\bar{\varepsilon}_{eq} - \varepsilon_{d0})} & \text{if:} \varepsilon \ge \varepsilon_{d0} \end{cases}, \quad (5) $$

where $B_t$ defines the shape of the postpeak behaviour law. One can notice that theoretically the value of 1 for damage is never reached since it requires an infinite equivalent strain $\bar{\varepsilon}_{eq}$.

## 2.1.2 | Local, SD, and WD damage formulations

Within a nonlocal formulation framework,³²,³³ the estimation of the equivalent strain $\bar{\varepsilon}_{eq}$ requires a general spatial weighting function $\phi$. As shown in Equation 6, $\phi$ depends on an internal length $l_c$ and the developed stress field $\sigma$ (in the case of stress-based models³³). $l_c$ and $\sigma$ define a finite influence domain $D_c$ that is a material-dependent constant³⁴ representative of the size of the FPZ (assumed to be equal to 3 times the maximum aggregates' size $D_c = l_{FPZ} \approx 3d_{max}$²³). $D_c$ is, accordingly and by definition, representative of the area over which energy dissipation occurs for a given crack propagation. And $\phi$ defines the spatial distribution of such dissipation around the point of interest (ie, the point of localization—Figure 2).

<figure><img src="image_4.png" /><figcaption>**FIGURE 1** Parameters of the (A) linear (LIN), (B) bilinear (BIL), and (C) exponential (EXP) softening laws [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**FIGURE 2** (A) Objective ($D_c = l_{FPZ} = 6$ cm—Table 2) vs (B) subjective ($l_c = 6$ cm—Table 1) usage of the weighting functions in nonlocal models [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

$$ \varepsilon_{eq}(x) = \int_{\Omega} \bar{\phi}(x-y, l_c, \sigma(y)) \varepsilon_{eq}^{loc}(x) dy \bar{\phi} = \frac{\phi(x-y, l_c, \sigma)}{\int_{\Omega} \phi(x-y, l_c, \sigma) dy} \quad (6) $$

where $\varepsilon_{eq}^{loc}$ estimation involves the first and second strain tensor's invariants.³¹

Depending on the used weighting function $\phi$ (Table 1), different models are obtained. By setting $\phi$ equal to the Dirac function $\delta$, local damage models are obtained, the Green's function $G$ leads to the second-order gradient-enhanced models, which are descriptive of an SD gradient-based approach,²⁰,²¹ and the classical Gaussian function $\rho$ leads to the WD formulation²⁰,²¹ and stress-based (SWD) integral models.³³ It is also possible to derive a stress-based second-order gradient-enhanced models (SSD) when the Green's function is used instead if the Gaussian one. In these 2 last cases, the effective internal length $l'_c$ depends on, both, the internal length $l_c$ and the stress distribution around the point of interest. It is worth underlining that all damage models, herein, depend on a strain-based damage criterion ($\varepsilon_{eq}^{loc}$ in Equation 6); the “stress-based” terms refer rather to the energetic regularization (preventing mesh dependency—$\varepsilon_{eq}$ in Equation 6). Compared with SD and WD formulations, stress-based approaches SSD and SWD are more representative of the stress concentration ahead of the crack tip and unloading behind it. Stress-based formulations lead to a more physical regularization close to stress-free edges (geometrical boundaries or fully developed cracks). Also, they can

<table><caption>**TABLE 1** Examples of subjective weighting functions and their associated types of damage models</caption><thead><tr><th>Formulation</th><th>Weighting Function φ</th><th></th><th></th></tr></thead><tbody><tr><td>Local (LOC)</td><td>δ(x - y) = <table><tr><td>1, x - y = 0</td></tr><tr><td>0, x - y ≠ 0</td></tr></table></td><td></td><td>-</td></tr><tr><td>Second-order gradient enhanced (SD)</td><td>G<sub>ε</sub>(x - y) = e<sup>-α|x - y|/l'<sub>c</sub></sup></td><td>α = 1<br>α = 1/√2</td><td>Dufour et al<sup>20</sup><br>Peerlings et al<sup>21</sup></td></tr><tr><td>Integral nonlocal (WD)</td><td>ρ<sub>ε</sub>(x - y) = e<sup>- (α|x - y|/l'<sub>c</sub>)<sup>2</sup></sup></td><td>α = 2<br>α = 1/√2</td><td>Dufour et al<sup>20</sup><br>Peerlings et al<sup>21</sup></td></tr><tr><td>Stress-based integral nonlocal (SWD)</td><td>ρ<sub>σ</sub>(x - y) = e<sup>- (α|x - y|/l'<sub>c</sub>)<sup>2</sup></sup><br>l'<sub>c</sub> = l<sub>c</sub> f<sub>t</sub> √<table><tr><td>(sin<sup>2</sup>φ cos<sup>2</sup>θ / σ<sub>I</sub>(y)<sup>2</sup>) + (sin<sup>2</sup>φ sin<sup>2</sup>θ / σ<sub>II</sub>(y)<sup>2</sup>) + (cos<sup>2</sup>φ / σ<sub>III</sub>(y)<sup>2</sup>)</td></tr></table><sup>-1</sup> = l<sub>c</sub>(σ/f<sub>t</sub>)<br>φ and θ are angles in the principal stress basis (σ<sub>I</sub>, σ<sub>II</sub>, σ<sub>III</sub>)</td><td>α = 2</td><td>Giry et al<sup>33</sup></td></tr><tr><td>Stress-based second-order gradient enhanced (SSD)</td><td>G<sub>σ</sub>(x - y) = e<sup>-α|x - y|/l'<sub>c</sub></sup><br>l'<sub>c</sub>=l<sub>c</sub> f<sub>t</sub> √<table><tr><td>(sin<sup>2</sup>φ cos<sup>2</sup>θ / σ<sub>I</sub>(y)<sup>2</sup>) + (sin<sup>2</sup>φ sin<sup>2</sup>θ / σ<sub>II</sub>(y)<sup>2</sup>) + (cos<sup>2</sup>φ / σ<sub>III</sub>(y)<sup>2</sup>)</td></tr></table><sup>-1</sup> = l<sub>c</sub>(σ/f<sub>t</sub>)<br>φ and θ are angles in the principal basis (σ<sub>I</sub>, σ<sub>II</sub>, σ<sub>III</sub>)</td><td></td><td>-</td></tr></tbody></table>
prevent damage attraction by boundaries due to the truncation of the interaction domain of integration points without requiring a preconditioning of the denominator in Equation 6. Consequently, as the stress state is reduced due to damage, the effective internal length tends to zero and the weighting function tends to become singular and equal to a Dirac impulsion for the considered position. In other words, and for a cracked zone, the evolution of SWD and SSD damage models is equivalent to the one using local models.

Generally speaking, these weighting functions (Table 1) define the shape of the FPZ over an influence domain sized $D_c = \int \phi(u, l_c)du$ (Equation 7).

$$ 
\begin{align} 
\text{Influence domain equivalent to a segment sized } D_c: \quad & D_c = \int_{-\infty}^{\infty} \phi(u, l_c) du \nonumber \\ 
\text{Influence domain equivalent to a circle with a radius } D_c: \quad & D_c = \frac{1}{2D\sqrt{\pi}} \int_{-\infty}^{\infty} \phi(u, l_c) du \tag{7} \\ 
\text{Influence domain equivalent to a sphere with a radius } D_c: \quad & D_c = \sqrt[3]{\frac{3}{4\pi}} \int_{-\infty}^{\infty} \phi(u, l_c) du \nonumber 
\end{align} 
$$

Naturally, this influence domain should be the same as the designated FPZ size ($D_c = l_{FPZ}$) over which the shape is defined. This means that, for an objective description of energy dissipation along the FPZ (Figure 2A), the influence domain $D_c = l_{FPZ}$ should be independent from the used weighting function to prevent a biased computation of the damage variable $d$ and the equivalent strain $\bar{\epsilon}_{eq}(x)$ (Figure 2B). However, the internal length $l_c$ should be seen as a fitting parameter and needs to be defined by solving Equation 7. $l_c$ can be considered as a subjective internal length of the non-local models whereas $D_c = l_{FPZ}$ is a physical parameter that is objective regarding the nonlocal formulations. For the weighting functions in Table 1, the objective formulations in Table 2 are obtained.

In the case of local formulations and in the absence of an internal length, the influence domain $D_c$ is brought to the unit of the spatial discretization $h$ (also called the mesh characteristic length). This requires an additional regularization step to prevent mesh dependency. Yet it remains problematic in terms of the maximal cracks' number estimation, which is more related to the FPZ size than the spatial discretization $h$.

### 2.1.3 | Crack opening estimation

As damage is described according to a continuous approach, a postprocessing step is required to compute the equivalent crack opening in the considered continuous domain. Herein, the 1D strategy defined in Mattalah et al²⁵ and accounting for elastic unloading is retained (Equation 8). The equivalent cracking strain $\epsilon_{ck}$ then writes

$$ \epsilon_{ck} = d \cdot \langle \epsilon_1 \rangle_+, \quad (8) $$

<table><caption>**TABLE 2** Examples of objective weighting functions and their associated types of damage models</caption><thead><tr><th>Formulation</th><th>Weighting Function φ</th></tr></thead><tbody><tr><td>SD</td><td>G<sub>ε</sub>(x-y) = e<sup>-<sup>|x-y|</sup>&frasl;<sub>α<sub>FPZ</sub></sub></sup></td></tr><tr><td>WD</td><td>ρ<sub>ε</sub>(x-y) = e<sup>-<sup>(|x-y|)/β<sub>FPZ</sub></sup>2</sup></td></tr><tr><td>SWD</td><td>ρ<sub>σ</sub>(x-y) = e<sup>-<sup>|x-y|</sup>&frasl;<sub>βD'<sub>c</sub></sub></sup></td></tr><tr><td>D'<sub>c</sub>=l<sub>FPZ</sub></td><td>f<sub>t</sub>√<sup>(sin<sup>2</sup>φ cos<sup>2</sup>θ + <sup>sin<sup>2</sup>φ sin<sup>2</sup>θ</sup>&frasl;<sub>σ<sub>II</sub>(y)<sup>2</sup></sub> + <sup>cos<sup>2</sup>φ</sup>&frasl;<sub>σ<sub>III</sub>(y)<sup>2</sup></sub>)</sup><sup>-1</sup> =l<sub>FPZ</sub>(σ/f<sub>t</sub>)</td></tr><tr><td>SSD</td><td>φ and θ are angles in the principal stress basis (σ<sub>I</sub>, σ<sub>II</sub>, σ<sub>III</sub>)<br>G<sub>σ</sub>(x-y) = e<sup>-<sup>|x-y|</sup>&frasl;<sub>aD<sub>c</sub></sub></sup><br>D'<sub>c</sub> = l<sub>FPZ</sub> f<sub>t</sub>√<sup>(sin<sup>2</sup>φ cos<sup>2</sup>θ + <sup>sin<sup>2</sup>φ sin<sup>2</sup>θ</sup>&frasl;<sub>σ<sub>II</sub>(y)<sup>2</sup></sub> + <sup>cos<sup>2</sup>φ</sup>&frasl;<sub>σ<sub>III</sub>(y)<sup>2</sup></sub>)</sup><sup>-1</sup> =l<sub>FPZ</sub>(σ/f<sub>t</sub>)<br>φ and θ are angles in the principal stress basis (σ<sub>I</sub>, σ<sub>II</sub>, σ<sub>III</sub>)</td></tr></tbody></table>
where $\epsilon_1$ is the elastic strain in the first principal direction (supposed to be the mode I cracking direction) and $\langle x \rangle_+ = \max(x; 0)$.

And for all weighting functions in Table 2 (in addition to local models), the equivalent crack opening $w_{ck}$ is computed according to Equation 9.²⁰

$$w_{ck} = \epsilon_{ck} \frac{\int_{\Omega} \phi \, dy}{\phi(0)} = \epsilon_{ck} D_c. \qquad (9)$$

Equation 9 remains also true for local modelling strategies considering an influence domain equal to the spatial discretization unit $D_c = h$. One can notice that, for nonlocal formulations, by using Equation 9 and Table 2, such definition does not depend on the used weighting function nor the subjective internal length $l_c$ ($D_c = l_{FPZ} = c^{ste}$). Also, it is worth underlining that the crack opening estimation is based on the integrated strain over the influence domain $D_c$. This leads to the implicit hypothesis of one crack per influence domain. Physically, such hypothesis is only convenient if the influence domain is related to the size of the FPZ. Consequently, for nonlocal models, it suffices to choose $D_c = l_{FPZ}$ to be physically representative. However, for local models, the influence domain is equal to the mesh characteristic size $h$ (user-defined), which can be bigger or smaller than $l_{FPZ}$. From a conceptual point of view and following a 1D tensile load, this can be problematic in the presence of nonnull residual stresses leading to multicracking: For $h \ll l_{FPZ}$, the number of cracks per $l_{FPZ}$ might be overestimated as more than one FE within the FPZ length can be damaged. On the contrary, for $h \gg l_{FPZ}$, and following an inverse reasoning, the number of cracks can be underestimated. Eventually, this means that local models can only be predictive regarding the number of cracks if the mesh discretization is equal to the FPZ size. Otherwise, an additional hypothesis has to be considered requiring a prior knowledge of the cracks' number $n_{fiss}$ per FE—in which case each crack opening writes $w_{ck} = \frac{\epsilon_{ck}}{n_{fiss}} l_{FPZ}$.

### 2.1.4 | Hillerborg energy-based criterion

Another issue related to the definition of the crack opening is its dependence on the used softening law (and damage evolution law) with respect to the considered fracture energy. In the case of nonlocal formulations, the weighting functions insure a mesh independency of the dissipated energy. However, for local models, mesh independency is not guaranteed by default and requires the application of regularization techniques such as the Hillerborg energy-based criterion.²⁴ It allows the implementation of a material descriptive behaviour law that energy is constant ($\sigma, \epsilon$) regardless of the used spatial discretization (case of local approaches) and regardless of the postpeak softening laws (case of local and nonlocal models). Practically, the criterion (Equation 10) leads to the identification of the postpeak parameters for a given influence domain $D_c$ and a given fracture energy $G_F$. In the case of local models, the Hillerborg criterion has a true regularization meaning since $D_c = h$ is mesh dependent. Whereas for nonlocal formulations, Equation 10 becomes mesh independent as $D_c = l_{FPZ}$ regardless of the used spatial discretization and should be looked at as an identification formula of the postpeak parameters.

This writes for all damage models in Table 2 and softening laws (Equations 3-5):

$$\sigma \le f_t \quad \forall \epsilon$$

$$\text{For non-local models: } \frac{G_F}{l_{FPZ}} = \int_{\Omega} \left( \int_{0}^{+\infty} \bar{\phi}(\mathbf{x}-\mathbf{x}_0, l_{FPZ}, \sigma(\mathbf{x}_0)) \sigma \, d\epsilon \right) d\mathbf{x} \quad \text{For local models: } \frac{G_F}{h} = \int_{0}^{+\infty} \sigma \, d\epsilon, \qquad (10)$$

where $G_F$, $f_t$, and $E$ are known parameters identified at the specimen scale. One can notice that the residual stress when damage is equal to 1 needs to be null for the dissipated energy to be finite and defined.

The application of Equation 10 leads to different validity domains of the postpeak softening laws depending on the damage formulation (local or nonlocal). For the sake of conciseness, the theoretical developments and their analyses are given in Appendix A. Eventually, one can refer to Tables 3 and 4 (summarizing the validity domain of the softening laws described above [Equations 3-5]) to select the convenient softening law and damage modelling strategy in Table 2. For other softening laws, the same analysis should be performed. An insight on the effects of regularization and damage formulation on the local and global mechanical responses is highlighted in Appendix B based on an FE analysis of a 1D beam tensile test. It shows how the energy-based regularization restrains spatial discretization and improves damage modelling objectivity. However, one should keep in mind that energy-based regularization does not ensure local
objectivity in terms of damage evaluation within the FPZ (better objectiveness is achieved for the crack opening evaluation at the point of localization) nor the global response objectivity in the presence of structural effects such as snap backs. This issue also affects the permeability evaluation if based on the damage variable as it will be detailed in the next section.

## 2.2 | Limitations of damage-based descriptions of permeability

As mentioned in the introductory part, 2 transfer modes have been extensively studied: a first one through a low and diffuse damaged medium and a second through a localized crack.

<table><caption>TABLE 3 Relationship between the postpeak parameters, influence domain, and fracture energy</caption><tbody><tr><td>LOC<br>SD<br>WD</td><td>$$ \begin{cases} \frac{G_F}{D_c E} = \frac{\varepsilon_{d0}^2}{2} + \int_{\varepsilon_{d0}}^{+\infty} (1-d)\varepsilon \, d\varepsilon \\ \sigma \le f_t \quad \forall \varepsilon \end{cases} $$</td><td>LIN</td><td>$$ D_c \in \left[ 0; 2 \frac{G_F E}{f_t^2} \right] \quad \begin{cases} \varepsilon_{d0} = \frac{f_t}{E} \\ \varepsilon_u = \frac{2G_F}{f_t D_c} \end{cases} $$</td></tr><tr><td></td><td>$$ \begin{cases} \frac{G_F}{D_c E} = \frac{\varepsilon_{d0}^2}{2} + \int_{\varepsilon_{d0}}^{+\infty} (1-d)\varepsilon \, d\varepsilon \\ \sigma \le f_t \quad \forall \varepsilon \end{cases} $$</td><td>BIL</td><td>$$ D_c \in \left[ 0; \frac{4G_F E}{5f_t^2} \right] \quad \begin{cases} \varepsilon_{d0} = \frac{f_t}{E} \\ \varepsilon_u = \frac{4G_F}{5f_t D_c} \\ \varepsilon_u = \frac{4 + 6n G_F}{5f_t D_c} \end{cases} $$</td></tr><tr><td></td><td>$$ \begin{cases} \frac{G_F}{D_c E} = \frac{\varepsilon_{d0}^2}{2} + \int_{\varepsilon_{d0}}^{+\infty} (1-d)\varepsilon \, d\varepsilon \\ \sigma \le f_t \quad \forall \varepsilon \end{cases} $$</td><td>EXP</td><td>$$ D_c \in \left[ D_{c, \min}; \frac{2G_F E}{5f_t^2} \right] \quad \begin{cases} B_t = \frac{\varepsilon_{d0} + \sqrt{\frac{4G_F}{E D_c} - \varepsilon_{d0}^2}}{2\left(\frac{G_F}{E D_c} - \frac{\varepsilon_{d0}^2}{2}\right)} \\ \varepsilon_{d0} = \frac{\ln\left(\frac{f_t}{E}B_t\right) + 1}{B_t} \end{cases} \quad \text{Requires numerical solving} $$</td></tr><tr><td></td><td>$$ \begin{cases} \frac{G_F}{D_c E} = \frac{\varepsilon_{d0}^2}{2} + \int_{\varepsilon_{d0}}^{+\infty} (1-d)\varepsilon \, d\varepsilon \\ \sigma \le f_t \quad \forall \varepsilon \end{cases} $$</td><td>EXP</td><td>$$ D_c \in \left[ \frac{2G_F E}{5f_t^2}; \frac{2G_F E}{f_t^2} \right] \quad \begin{cases} B_t = \frac{\varepsilon_{d0} + \sqrt{\frac{4G_F}{E D_c} - \varepsilon_{d0}^2}}{2\left(\frac{G_F}{E D_c} - \frac{\varepsilon_{d0}^2}{2}\right)} \\ \varepsilon_{d0} = \frac{f_t}{E} \end{cases} $$</td></tr><tr><td>SSD</td><td>$$ \begin{cases} \frac{G_F}{E} = l_{FPZ} \frac{\varepsilon_{d0}^2}{2} + \int_{\Omega} \left( \int_{\varepsilon_{d0}}^{+\infty} (1-d)e^{-\frac{2|\mathbf{x}-\mathbf{x}_0|}{E} (1-d)\varepsilon} \right)^{l_{FPZ}} \varepsilon \, d\varepsilon \\ \sigma \le f_t \quad \forall \varepsilon \end{cases} $$</td><td></td><td>Requires numerical solving</td></tr><tr><td>SWD</td><td>$$ \begin{cases} \frac{G_F}{E} = l_{FPZ} \frac{\varepsilon_{d0}^2}{2} + \int_{\Omega} \left( \int_{\varepsilon_{d0}}^{+\infty} (1-d)e^{-\pi \left( \frac{|\mathbf{x}-\mathbf{x}_0|}{E} (1-d)\varepsilon \right)^2} \varepsilon \, d\varepsilon \right) \\ \sigma \le f_t \quad \forall \varepsilon \end{cases} $$</td><td></td><td>Requires numerical solving</td></tr></tbody></table>
<table><caption>**TABLE 4** Validation domains of the softening laws and local/nonlocal damage modelling approaches</caption><thead><tr><th rowspan="2">Spatial discretization</th><th colspan="2">Possible modelling strategy</th><th colspan="4">Possible softening law</th></tr><tr><th>Non-local (D<sub>c</sub> = l<sub>FPZ</sub>)</th><th>Local (D<sub>c</sub> = h)</th><th>D<sub>c</sub> ∈ [0; D<sub>c, min</sub>]</th><th>D<sub>c, min</sub> / (2G<sub>F</sub>E / 5 f<sub>t</sub><sup>2</sup>)</th><th>D<sub>c, min</sub> / (2G<sub>F</sub>E / 5 f<sub>t</sub><sup>2</sup>) / (4G<sub>F</sub>E / 5 f<sub>t</sub><sup>2</sup>)</th><th>D<sub>c, min</sub> / (4G<sub>F</sub>E / 5 f<sub>t</sub><sup>2</sup>) / (2G<sub>F</sub>E / 5 f<sub>t</sub><sup>2</sup>)</th></tr></thead><tbody><tr><td>h &lt; l<sub>FPZ</sub></td><td>yes</td><td>yes</td><td>LIN</td><td>yes</td><td>yes</td><td>yes</td></tr><tr><td>h = l<sub>FPZ</sub></td><td>yes but not required</td><td>yes</td><td>BIL</td><td>yes</td><td>yes</td><td>no</td></tr><tr><td>h &gt; l<sub>FPZ</sub></td><td>no</td><td>yes</td><td>EXP</td><td>no</td><td>yes but not recommended*</td><td>yes</td></tr></tbody></table>

### 2.2.1 | Permeability through a low and diffuse damaged medium

For diffuse damage, several models are available (Table 5) based on the fitting of experimental results.<sup>3-5,9,10,36,37</sup> These models directly relate the damage variable *d* to the apparent intrinsic permeability denoted *k*<sub>D</sub> (D in *k*<sub>D</sub> stands for Darcy's mode). Even though they offer interesting insights on the evolution of permeability with damage, their general use remains limited given their strong dependency on the experimental set-ups and loading states. On the one hand, the homogeneity of the domain should not be jeopardized to define a proper behaviour law (hence the validity limitation to low damage values). On the other hand, the definition of damage based on experimental results remains— at best—in a mean sense (not a local one) and still requires the definition of a damage model that eventually affects the permeability estimation. For instance, models in Choinska et al<sup>3,4</sup> and Picandet et al<sup>5</sup> have been both developed based on compressive tests and for small and diffuse damage values. However, lower permeability values are measured in Choinska et al<sup>3,4</sup> given that the performed analysis covered the whole specimen instead of the central part where damage increases. In that regard, results from Picandet et al<sup>5</sup> are more reliable since they cover and characterize only the damaged zone. Assuming the equivalence between compressive and tensile damage at low damage values, tests in previous studies<sup>3-5</sup> led to the definition of a damage-based exponential permeability law for 0 ≤ *d* ≤ 0.15. As *d* increases, the exponential law shows a considerable overestimation of the apparent permeability. Consequently, the Taylor expansion of the exponential law at *d* = 0 is preferred (up to the third order) herein as suggested in Choinska et al.<sup>9</sup> Another issue raised by experimental results in Choinska et al<sup>3,4</sup> is the existence of irreversible *k*<sub>D</sub><sup>i</sup> and reversible *k*<sub>D</sub><sup>r</sup> parts of damage-induced permeability at low levels (Equation 11). Such irreversible tendencies have also been observed in Bary<sup>35</sup>

<table><caption>**TABLE 5** Examples of low damage-permeability coupling</caption><thead><tr><th>Model k<sub>D</sub>(d) With d: Damage Variable and k<sub>0</sub>:<br>Intrinsic Permeability</th><th>Comment/Ref.</th></tr></thead><tbody><tr><td>k<sub>D</sub> = k<sub>0</sub>e<sup>(&alpha;d)<sup>&beta;</sup></sup></td><td>Picandet et al<sup>5</sup>; based on a global analysis of damage in a compressive test (unloaded specimen)<br>Choinska et al<sup>3,4</sup>; based on a global analysis of damage in a compressive test (unloaded specimen)<br>Bary<sup>35</sup>; based on BIPEDE tensile test results<sup>36</sup></td></tr><tr><td>k<sub>D</sub> = k<sub>0</sub> + 1 + &sum;<sup>3</sup><sub>n=1</sub> ((&alpha;d)<sup>n&beta;</sup> / n!)</td><td>Mozayan Kharazi et al<sup>10</sup>; developed from Picandet et al<sup>5</sup> to include the tensile strains induced by axial compressive loads (Poisson effect)</td></tr><tr><td>k<sub>D</sub> = k<sub>0</sub>e<sup>(&alpha;d)<sup>&beta;</sup></sup> (1 + (&alpha;'d)<sup>&beta;</sup>)</td><td>Picandet et al<sup>5</sup>; Taylor expansion to prevent the overestimation of permeability at high damage values (unloaded specimen)<br>Choinska et al<sup>3,4</sup>; Taylor expansion to prevent the overestimation of permeability at high damage values (unloaded specimen)</td></tr><tr><td>k<sub>D</sub> = k<sub>0</sub>10<sup>&alpha;d</sup></td><td>Choinska et al<sup>3,4</sup>; based on a global analysis of damage in a compressive test (loaded specimen)<br>Gawin et al<sup>37</sup>; based on BIPEDE tensile test results in Bary<sup>35</sup> and Gérard<sup>36</sup></td></tr></tbody></table>

<footnote>*shows a hardening behaviour and underestimates the damage threshold ε<sub>d0</sub> ≤ f<sub>t</sub>/E.</footnote>
and Gérard³⁶ through a BIPEDE tensile test (concrete disc subjected to uniaxial tension—Figure 8) where advanced tensile loading has been applied until the development of macrocracks. This justifies the higher values of damage (up to 80%) that are considered in the models suggested in Gawin et al³⁷ and Dal Pont et al³⁸ based on the BIPEDE test.³⁶ Such models are debatable since at this high damage level, strain localization leads to the development of macrocracks that are governed rather by the Poiseuille's mode than Darcy's. As for the observed residual permeability and in the presence of macrocracks, it can be due to a given residual crack opening or, more likely, to a residual diffuse damage in the cracked domain. Indeed, for a complete unloading, the macrocrack can be considered fully closed whereas the network of microcracks around remains active. Such configuration favours the use of Darcy's law rather than Poiseuille's, which is valid for crack openings that are at least higher than 20 µm.⁶ Accordingly, the irreversibility of permeability is herein associated with concrete damage and not with a residual macrocrack opening.

Eventually, for whatever damage-permeability coupling model, one can notice that the variation of their parameter has more effect on the computed permeability compared with the effect of the damage softening law variation (Figure 3).

For the purpose of further developments, the reference permeability law retained herein is the following:

$$k_D = \begin{cases} k_0 \left( 1 + \sum_{n=1}^{3} \frac{(\alpha d)^{n\beta}}{n!} \right) \left( 1 + (\alpha' d)^{\beta'} \right); & \text{if loaded} \\ k_0 \left( 1 + \sum_{n=1}^{3} \frac{(\alpha d)^{n\beta}}{n!} \right) & \text{; if unloaded} \end{cases} \quad \text{with } 0 \le d \le 0.15 \quad (11)$$

with values issued from experimental results in Picandet et al⁵ ($\alpha$ = 11.3, $\beta$ = 1.6) and based on experimental observations in Choinska et al³,⁴ $\alpha'$ ≈ $\alpha$, $\beta'$ ≈ $\beta$.

It is important to underline that the previous permeability law has been derived in a local damage modelling sense despite the fact that only the mean response (mean damage of the influence domain) is accessible experimentally. Consequently, and from a numerical point of view, the computed permeability over the influence domain (Equation 12) is strongly dependent on the damage formulation (LOC/SW/SD/SSW/SSD in Figure B4). This has a nonnegligible effect on $k_D$ estimation (Figure 4): For nonlocal formulations and in the case of an unloaded concrete at $d$ = 15%, $k_D$ is 10 times lower than the reference value (LOC descriptive of the FPZ mean response). Gaps are 15 times higher for $d$ = 100%. This simply means that the fitting coefficients in Table 5 for Equation 11 are only valid for local models and need to be redefined in the case of nonlocal modelling of damage. Another alternative, more physically representative of the FPZ's mean response, is to associate Darcy's transfer mode with a mean quantity that is less damage dependent (at least from a local point of view), for instance, the mean damage value that would require an additional

<figure><img src="image_6.png" /><figcaption>**FIGURE 3** Comparison between different low damage-permeability coupling laws using LIN/BIL/EXP softening laws and LOC model with ($h = D_c = l_{FPZ} = 6$ cm) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 4 Effect of damage formulation on the computed damage-permeability coupling in the case of Darcy's transfer mode for an unloaded concrete $k_D = k_0 + \sum_{n=1}^{3} \frac{(11.3d)^{n-1.64}}{n!}$ using (A) LIN, (B) BIL, and (C) EXP softening laws</figcaption></figure>

postprocessing as performed for the crack opening in Equation 9 or directly the equivalent crack opening that would be representative (in a hydric sense) of the flow through the damaged volume.

Another drawback of Equation 11, also related to its damage dependence, is the absence of a path between loaded specimen and unloaded one that hinders the description of permeability evolution during a loading-unloading cycle.

$$k_D^{\text{FPZ}} = \frac{1}{D_c} \int_{-D_c/2}^{+D_c/2} k_D(d(x)) \, dx. \qquad (12)$$

### 2.2.2 | Permeability through a macrocracked medium

Poiseuille's law is usually used to describe permeability through a macrocrack$^{5,6,8,11,39}$ (Equation 13). It consists of writing the apparent permeability $k_P$ (P in $k_P$ stands for Poiseuille) as a function of the squared crack opening value$^8$ (which is unbounded in opposition with the damage variable):

$$k_P = \xi \frac{w_{ck}^2}{12} = \frac{w_{ck,eff}^2}{12}, \qquad (13) \\ w_{ck,eff} = \sqrt{\xi w_{ck}}$$

where $\xi \le 1$ is a factor to account for the roughness and shape effects of the through crack domain.$^{5,6}$ When $\xi = 1$, the domain is supposed to be perfectly smooth. The expression suggested in Rastiello et al$^6$ based on real-time water permeability evolution through a localized crack during a splitting test (Equation 14) seems more appropriate compared to others where $\xi$ is considered constant (between 0.01 and 0.1 in Picandet et al$^5$). Indeed, as the crack opening becomes larger, the roughness effects are supposed to be relatively less and $\xi$ should tend towards 1 (Figure 5). It is worth mentioning that experimental results cover only crack openings up to 160 µm whereas the condition of $\xi \le 1$ is verified until $w_{ck} = 268$ µm. So an extrapolation is performed for crack opening values between 160 and 268 µm.

$$\xi = \min((\alpha_w w_{ck})^{\beta_w}; 1) \qquad (14) \\ (\alpha_w, \beta_w) = (3726.71, 1.19) \text{ for } 20 \ \mu\text{m} \le w_{ck} \le 160 \ \mu\text{m}$$

Numerically, a macrocracked FE domain sized $D_c$ contains a macrocrack of an opening $w_{ck}$ and also uncracked domain around sized $D_c - w_{ck}$. The permeability $k_F$ (F stands for fracture in $k_F$) when considering the cracked domain only writes$^{40}$

<footnote>$^8$Hypothesis: crack = 2 parallel planes</footnote>
<figure><img src="image_8.png" /><figcaption>**FIGURE 5** Crack's roughness effect on (A) the effective crack opening and (B) the through crack permeability</figcaption></figure>

$$k_F = k_P \varepsilon_{ck} = k_P \frac{w_{ck}}{D_c} = \frac{1}{12 D_c} \min \left( \alpha_w \beta_w w_{ck}^{3+\beta_w}; w_{ck}^3 \right). \quad (15)$$

### 2.2.3 | Damage-based matching laws

Permeability matching laws describe the evolution of the permeability from a lowly diffuse damage state ($d < 15\%$) to a macrocracked state governed by Poiseuille's law. Even if experimental investigation of this intermediate domain lacks, several matching strategies have been proposed based on mathematical extrapolations between Darcy's and Poiseuille's modes (Table 6).

A comparison between the different matching laws is depicted in Figure 6A-C:

1. In literature,⁶,⁹,¹⁰,⁴² one can find a linear matching between the 2 modes where the increase of permeability, when damage exceeds 15%, is only due to the increase of the crack opening value. However, for values of damage between 15% and 90%, the effective crack opening remains quasinull, which leads to a constant permeability within that domain.

2. Other linear models suggest a damage-weighted addition of the 2 modes.³,⁹ This has the advantage of not having to limit the damage value to 15% since at high damage values, only the Poiseuille's mode is activated. Still, this solution seems to underestimate permeability at low damage values when the crack's roughness effect is considered (timorous lines in Figure 6A-C). On the other hand, when roughness is not considered (bold lines in Figure 6A-C), the damage-pondered linear law overestimates Darcy's permeability at high damage values leading to a non-physical decrease of permeability as strain localizes and Poiseuille's mode is activated.

3. Lastly, the damage-pondered logarithmic law³,⁴¹ (Equation 16) offers a smoother transition from Darcy's to Poiseuille's mode. However, when the crack's roughness effects are accounted for, it shows a nonphysical decrease of permeability as damage is initiated and also underestimates the resulting permeability for $0 \le d \le 0.15$. In the

<table><caption>**TABLE 6** Examples of damage-permeability matching laws</caption><thead><tr><th>Model k = f(k<sub>D</sub>, k<sub>F</sub>)</th><th>Ref.</th></tr></thead><tbody><tr><td>k = k<sub>D</sub> + k<sub>F</sub></td><td>k<sub>D</sub> = k<sub>0</sub>e<sup>(&alpha;<sub>w</sub>&radic;2&nu;<sub>d</sub>)<sup>&beta;</sup></sup>; d &le; 0.15<br>k<sub>F</sub> = &xi;<sup>w<sub>ck</sub><sup>3</sup></sup>&frasl;<sub>12 D<sub>c</sub></sub><br>k<sub>D</sub> = k<sub>0</sub><br>k<sub>F</sub> = &xi;<sup>w<sub>ck</sub><sup>3</sup></sup>&frasl;<sub>12 D<sub>c</sub></sub></td><td>Mozayan Kharazi et al<sup>10</sup><br>Rastiello et al<sup>6</sup> and Asali et al<sup>41</sup></td></tr><tr><td>k = (1 &minus; d)k<sub>D</sub> + d k<sub>F</sub></td><td>k<sub>D</sub> = k<sub>0</sub>e<sup>(&alpha;<sub>w</sub>d)<sup>&beta;</sup></sup><br>k<sub>F</sub> = &xi;<sup>w<sub>ck</sub><sup>3</sup></sup>&frasl;<sub>12 D<sub>c</sub></sub><br>w<sub>ck</sub> = f(d)</td><td>Choinska et al<sup>3,9</sup></td></tr><tr><td>k = k<sub>D</sub><sup>1 &minus; d</sup>k<sub>F</sub><sup>d</sup></td><td>k<sub>D</sub> = k<sub>0</sub>e<sup>(&alpha;<sub>w</sub>d)<sup>&beta;</sup></sup><br>k<sub>F</sub> = &xi;<sup>w<sub>ck</sub><sup>3</sup></sup>&frasl;<sub>12 D<sub>c</sub></sub></td><td>Ezzedine El Dandachy<sup>42</sup></td></tr></tbody></table>
case of the BIL softening law, the decrease also occurs at $d \sim 80\%$ as Poiseuille's mode has more weight even if its contribution remains less than Darcy's.

Herein, the log-type matching is retained for further developments (Equation 16). Two main issues need to be addressed: The first one is the weighting of Darcy's and Poiseuille's modes so as to prevent permeability decrease with damage. And the second is the damage dependence of the matching law preventing an accurate description of permeability as concrete is unloaded. Indeed, for a fully damaged element ($d = 1$), the only active mode is Poiseuille's. During an unloading cycle, the damage variable due to tensile loads remains equal to one (one is reminded that the damage variable can only increase) whereas the crack opening and Poiseuille's flow decrease. Such configuration leads to an unrealistic residual permeability that is either overestimated ($k = k_D(d = 1)$) in the case of linear matching³ or null in the case of damage-pondered laws ($k = k_F(d = 1) = 0$ in the absence of any residual crack opening).⁴²

$$k = k_F^d k_D^{1-d} \qquad (16)$$

## 3 | NEW STRAIN-BASED DESCRIPTION OF DAMAGE-PERMEABILITY COUPLING

The previous section has presented different damage-based matching laws and the related drawback that can be identified. Several enhancements can be suggested aiming at a more physical description of permeability evolution from a low diffuse damage to a highly localized damage state. Particularly, the foreseen law should

<figure><img src="image_9.png" /><figcaption>**FIGURE 6** Damage-based matching laws using (A,D) LIN, (B,E) BIL, and (C,F) EXP softening laws with $k_0 = 10^{-17}m^2$, $D_c = 6$ cm, $k_D = k_0 \left( 1 + \sum_{n=1}^{3} \frac{(11.3 d)^{n1.64}}{n!} \right) (1 + (11.3d)^{1.64})$, and $k_F = \xi \frac{w_{ck}^3}{12 D_c}$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
* limit damage model dependence on the permeability law (especially for Darcy's transfer mode),

* ensure a monotonously increasing permeability with damage,

* distinguish reversible and irreversible parts of permeability, and

* account for loading-unloading cycles.

## 3.1 | Matching law

1. Strain-based permeability for Darcy's transfer mode: Let us first express the permeability as the sum of a reversible and an irreversible part in line with experimental results in literature.3,4,35,36 The reversible permeability due to damage $k'_D$ would be associated with the transfer properties of the reversible closure and opening of microcracks or the elastic behaviour of the damaged porous network subjected to loading. The irreversible part $k'_D$ would be descriptive of the fluid's flow through the residual and irreversible openings of the mechanically microcracked domain. In opposition to the reversible part, the irreversible one can only increase in time depending on the loading level. Its maximum value is obtained when the macrocracked state is reached. At this stage, it is supposed that, in the vicinity of the crack lips, concrete is unloaded and the maximum damage level has been achieved. In other words, the microcracks network around the macrocrack is within its residual mechanical state. From that point on, any evolution of the transfer properties is directly associated with the variation of the macrocrack opening. As a consequence, and within the suggested framework, the residual opening of macrocracks is supposed null; irreversibility is rather associated with the microcracked network. In a hydraulic sense, it is possible to define an equivalent residual macrocrack opening leading to the same permeability evaluation. However, from a physical point of view, a residual opening of the macrocrack would be more representative of the crack lips' mismatch during the closure or the surface alteration (during a compressive cycle) and should be addressed separately. In Equation 11, the total permeability $k_D$ and irreversible one $k'_D$ have already been defined. So what lacks is the definition of the reversible part $k'_D = k_D - k'_D$, a loading path between the total and irreversible permeability values and a state variable to prevent damage dependence of $k_D$. This last issue is particularly of interest to limit the dependence of Darcy's transfer mode modelling on the damage formulation (local or nonlocal—Figure B4).

$$k_D = k_D^i + k_D^r \qquad (17)$$

with $k_D^i(d) = k_0 \left(1 + \sum_{n=1}^3 \frac{1}{n!} (\alpha d)^{n\beta}\right)$ for $0 \le d \le 15\%$ and $k_D^r = \left(1 + (\alpha' d)^{\beta'}\right) \left(1 + \sum_{n=1}^3 \frac{1}{n!} (\alpha d)^{n\beta}\right) - k_D^i$. For a macrocracked area, the residual permeability would be then equal to $k_0 \left(1 + \sum_{n=1}^3 \frac{1}{n!} (\alpha)^{n\beta}\right)$. With $k_0 = 10^{-17}$, $\alpha = 11.3$, and $\beta = 1.64$, the obtained residual permeability is $k_D^i(d) = 2.6 \times 10^4 k_0$, which is equivalent to a flow through an effective crack opening of 57 µm. This is obviously not physically admissible. Therefore, $k_D^i(d)$ needs to be bounded according to a certain maximum damage threshold $d_{lim}$ from which the residual permeability remains constant. Based on the validity domain of Poiseuille's law in Equation 15 ($k_D^i(d_{lim}) \le k_F(w_{ck} = 20 \mu m)$) and Darcy's law in Equation 11, $d_{lim}$ verifies: $0.2 \le d_{lim} \le 0.5$ and the irreversible part of permeability writes

$$ \begin{aligned} k_D^i &= k_0 \left( 1 + \sum_{n=1}^3 \frac{1}{n!} (\alpha \overline{d_n})^{n\beta} \right) \\ \overline{d_{lim}} &= \gamma \left| \overline{d^r} \right|_t + (1 - \gamma) \min \left( \left| \overline{d^r} \right|_t ; d_{lim} \right) \\ \overline{d^r} &= \begin{cases} 0; & \text{if } \varepsilon_I^+ = 0 \\ \frac{w_{ck}}{\varepsilon_I^+ D_c}; & \text{if } \varepsilon_I^+ > 0 \end{cases} \end{aligned} \qquad (18) $$

with $0 \le \overline{d^r} \le 1$ a “reversible” and strain-based variable to prevent damage dependence of Darcy's law and $d_{lim}$ a fitting parameter depending on the measured residual permeability. $\gamma$ is a loading path parameter varying between 0 and 1 chosen herein arbitrarily as $0 \le \gamma = \frac{\varepsilon_{eq}}{\varepsilon_{eq}} \le 1$: when $\varepsilon_{eq} = \overline{\varepsilon_{eq}}$, $\gamma = 1$, and $\overline{d_{lim}} = \left| \overline{d^r} \right|_t$—during the unloading phase: $\gamma = \frac{\varepsilon_{eq}}{\varepsilon_{eq}} \to 0$ and $\overline{d_{lim}} \to \min \left( \left| \overline{d^r} \right|_t ; d_{lim} \right)$.
Accordingly, the reversible part of Darcy's contribution writes

$$k_D^r = k_0 (\alpha' \overline{d^r})^{\beta'} \left(1 + \sum_{n=1}^{3} \frac{1}{n!} (\alpha \overline{d^r})^{n\beta}\right). \quad (19)$$

2. Strain-based permeability for Poiseuille's transfer mode: Poiseuille's mode in Equation 15 is already expressed in terms of the crack opening rather than damage. Such formulation does not require any additional modification. The challenge regarding Equation 15 remains related to the experimental quantification of the crack's roughness and tortuosity effects on the measured Poiseuille's flow.

3. Log-type matching: As demonstrated in Section 2.2.3, the use of damage as a weighting function does not allow a proper description of permeability at low damage values, as it does not ensure an increasing permeability with damage and, being monotonously increasing, does not allow an accurate description of the permeability of an unloaded concrete. Accordingly, the foreseen weighting function needs to be reversible and only active within the intermediate zone $0.15 \le d < 0.99$. Such pondering can be achieved using the previously defined “reversible” damage variable (using $\epsilon_{eq}$ instead of $\overline{\epsilon}_{eq}$ in Equations 3-5), and its associated weighting function has to be decreasing only when $0.15 \le d$. This leads to a new strain-based (and damage independent) permeability law (Equation 20):

$$k = (k_F)^{f^r} (k_D)^{1-f^r} \\ f^r = \left( \frac{1}{1-d_{\text{lim}}} <\overline{d^r} - d_{\text{lim}}>_+ \right)^{\delta} \quad (20)$$

with $\delta$ a fitting parameter defining the rate of pondering between the 2 modes. In the absence of further experimental results, $\delta$ should be adjusted to get an increasing permeability in the loading direction.

Finally, and as shown in Figure 7, the new matching law (Equation 20) is more accurate and physically acceptable for both the loading and unloading phases compared with the original logarithmic one (Figure 6): The permeability transition from Darcy's to Poiseuille's mode is achieved successfully with respect to the monotonously increasing permeability with damage—the permeability transition from the loaded to the residual (unloaded) configuration is properly defined so as to insure a nonnull residual permeability.

This new formulation remains valid whether unilateral effects and residual strains are present in the damage model or not.³¹ In a local sense (for the FPZ characterization), one should keep in mind that the fitting of Equation 20 depends on the used damage softening law. Particularly, the condition of increasing permeability leads to a smaller value of $\delta$ ($\delta = 0.5$) in the case of BIL law compared with LIN/EXP softening laws ($\delta = 0.8$). This is due to the slow evolution of damage after initiation when using the BIL law to meet with energetic criteria. However, still in a local sense, its dependence on the damage formulation remains limited given the strain-dependent nature of Equation 20. Indeed, from a theoretical point of view, the estimation of crack opening values is sensibly the same using local or nonlocal models

<figure><img src="image_10.png" /><figcaption>**FIGURE 7** New strain-based permeability law using (A) LIN $\delta = 0.8$, (B) BIL $\delta = 0.5$, and (C) EXP $\delta = 0.8$ softening laws with $k_0 = 10^{-17} m^2$, $D_c = 6$ cm, $\alpha' \approx \alpha$, $\beta' \approx \beta$, $d_{\text{lim}} \approx 0.15$, $\alpha = 11.3$, $\beta = 1.64$. Bold lines: $\alpha_w = 1 m^{-1}$, $\beta_w = 0$ and timorous lines: $\alpha_w = 3726.71 m^{-1}$, $\beta_w = 1.19$ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 8 Schematic description of the BIPEDE tensile test³⁶ [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

when the energetic criteria are respected (see Appendix B). The effect of damage formulation is expected to be more pronounced in the presence of structural effects (snap backs). In that case, different fittings would be required depending on the used damage formulation.

## 3.2 | Analytical validation based on the BIPEDE tensile test

Experimental results covering the evolution of concrete's permeability under loading-unloading cycles for low and high damage values are not numerous in the literature. In this contribution, the BIPEDE tensile tests developed in Gérard³⁶ is considered. It consists of applying uniaxial tensile loads to a disc glued to steel plates (Figure 8). Tension is applied to the plates and is transmitted to the disc until cracks develop.

The disc is 40 mm thick and has a diameter of 110 mm. The plates are holed so as to allow fluid (water) transfer through the concrete disc (hole's diameter is 55 mm). Test results are provided in terms of water flow through the concrete specimen $q$ (m/s).³⁶ The specimen's equivalent permeability can be estimated using Equation 21 supposing that the stationary state is constantly achieved.

$$k (\text{m}^2) = \frac{e\mu}{S \Delta P} q \qquad (21)$$

with $\Delta P = 0.25$ MPa, $e = 40$ mm, $\mu = \mu_{water} = 10^{-3}$ Pa.s, and $S = \pi \frac{\phi^2}{4} = 2375$ mm².

The purpose of this paragraph is not to show the validity of a given mechanical model but to demonstrate that the suggested matching law in Equation 20 can be fitted to realistic experimental results. It is worth mentioning that, experimentally, one or 2 cracks were obtained (not forcibly centred). Herein, one crack is supposed to develop and local strains at the point of localization are directly used to define—according to a given softening law—the resultant damage and crack opening values (explicit analytical scheme). Also, in the framework of this contribution, only homogeneous

<figure><img src="image_12.png" /><figcaption>FIGURE 9 Fitting of the suggested strain-based permeability law to experimental results for loading and unloading phases with $k_0 = 5 \times 10^{-22}$m², $D_c = 6$ cm, $\alpha_w = 10$. $m^{-1}$, $d_{lim} = 0.13$, $\beta_w = 1.19$, $\alpha' \approx \alpha$, $\beta' \approx \beta$, $\alpha = 11.3$, $\beta = 1.64$. (A) LIN, (B) BIL, and (C) EXP softening laws using strain-based approach (Equation 20) and damage-based one (Equation 16) [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
FE modelling strategies of concrete are mentioned. However, nothing prevents the use of the suggested strain-based matching law (Equation 20) in the case of mesoscopic modelling strategies where concrete is considered as a biphasic material (aggregates and cement paste).<sup>43</sup>

For the numerical application, the considered concrete properties are the following: $f_t = 4.5$ MPa, $G_f = 80$ N/m, and $E = 30$ GPa. The influence domain $D_c$ is estimated at $2 d_{max} = 6$ cm. In Figure 9, the adjustment of Equation 20 to experimental results is performed using the least squared method where only roughness effects ($\alpha_w$) and irreversibility parameter $d_{lim}$ are adjusted. In terms of crack opening evolution, analytical developments lead to the same results for LOC/SD/WD/SSD/SWD approaches when the energetic Hillerborg criterion is properly applied (Appendix B). By comparing the obtained strain-based curve (Equation 20) with damage-based one (Equation 16) in Figure 9, one can notice the accuracy of the newly suggested formula to describe permeability evolution during both the loading and unloading phases using LIN/BIL and EXP softening laws. For this test, the roughness effect on the computed Poiseuille's permeability is less pronounced compared with the one in Rastiello et al.<sup>6</sup> For that reason, the decrease of permeability with damage using Equation 16 is limited at low damage values. As shown in Figure 6, this drawback of damage-based approaches is more pronounced as the roughness and tortuosity effects increase. As for the unloading phase, it is clear that Equation 20 leads to a more physically admissible results as Equation 16 provides a null residual permeability. The strong scattering of experimental data as concrete is fully unloaded does not allow a clear understanding of whether the loading and unloading permeability paths are superimposed or not. With that regard, an in-depth experimental investigation of permeability evolution of an unloaded concrete is still lacking. For this contribution, the path is considered reversible and no hysteresis is included. This can be improved, if clearly demonstrated by experimental results, by modifying the shape of Equation 19 for the unloading phase accordingly. Nevertheless, our main interest is served since the suggested strain-based permeability law in Equation 20 allows appropriate fitting of realistic experimental results regardless of the used damage modelling strategy (in the absence of structural effects).

## 4 | CONCLUSIONS

In this contribution, a new strain-based permeability log-type matching law is suggested instead of damage-based ones. Being exclusively strain dependent, the new law does not depend on unilateral effects or damage formulation (as far as the local behaviour of the FPZ is concerned). As demonstrated through the BIPEDE tensile test, it allows better description of permeability decrease during the unloading phase and ensures a monotonously increasing permeability as the loading advances. As convincing experimental works of permeability evolution during the unloading phase are lacking, no hysteresis is considered herein for the reversible permeability part. Our future work shall investigate concrete hydric behaviour more thoroughly for more than one loading-unloading cycle to better characterize the permeability during unloading phases.

## ACKNOWLEDGEMENTS

The first author is grateful for the support provided by EDF-SEPTEN/DTG/CIH within the Chair PERENITI agreement with the Grenoble INP Partnership Foundation. The Chair PERENITI partners shall not in any circumstances be deemed liable for the content of this publication which is only binding its author.

The second author is grateful for the financial support provided by France's National Research Agency (ANR-11-RSNR-012).

The 3SR lab is part of the LabEx Tec 21 (Investissements d'Avenir—Grant agreement No. ANR-11-LABX-0030).

## ORCID

D. E.-M. Bouhjiti ID http://orcid.org/0000-0001-6357-5205

S. Dal Pont ID http://orcid.org/0000-0002-0348-5918

## REFERENCES

1. Gérard B, Breysse D, Ammouche A, Houdusse O, Didry O. Cracking and permeability of concrete under tension. *Mater Struct.* 1996;29(187):141-151. https://doi.org/10.1007/BF02486159
2. Picandet V, Khelidj A, Bastian G. Effect of axial compressive damage on gas permeability of ordinary and high performance concrete. *Cem Concr Res.* 2001;31(11):1525-1532. https://doi.org/10.1016/S0008-8846(01)00546-4

3. Choinska M. Effet de la température, du chargement mécanique et de leurs interactions sur la perméabilité du béton de structure. PhD thesis ECN Nantes France 2006

4. Choinska M, Khelidj A, Chatzigeorgiou G, Pijaudier-Cabot G. Effects and interactions of temperature and stress-level related damage on permeability of concrete. *Cem Concr Res.* 2007;37(1):79-88. https://doi.org/10.1016/j.cemconres.2006.09.015

5. Picandet V, Khelidj A, Hervé B. Cracks effects on gas and water permeability of concretes. *Cem Concr Res.* 2009;39(6):537-547. https://doi.org/10.1016/j.cemconres.2009.03.009

6. Rastiello G, Tailhan JL, Rossi P, Dal Pont S. Macroscopic probabilistic cracking approach for the numerical modeling of fluid leakage in concrete. *Ann Solid Struct Mech.* 2015;7(1):1-16. https://doi.org/10.1007/s12356-015-0038-6

7. Rastiello G, Desmettre C, Tailhan JL, Dal Pont S. Modeling of fluid leakage through multi-cracked RC structural elements using a numerical probabilistic cracking approach. *Mater Struct.* 2016;49(8):3095-3108.

8. Pijaudier-Cabot G, Dufour F, Choinska M. Permeability due to the increase of damage in concrete: from diffuse to localize damage distributions. *J Eng Mech.* 2009;135(9):1022-1028. https://doi.org/10.1061/(ASCE)EM.1943-7889.0000016

9. Choinska M, Dufour F, Pijaudier-Cabot G. Matching permeability law from diffuse damage to discontinuous rack opening. *Proc. Fracture Mechanics of Concrete and Concrete Structures.* 2007;1:541-547.

10. Mozayan Kharazi M, Davenne L, Ghavamian S. A modeling methodology in predicting the leakage tightness of a pressure containment vessel without liner. *Proceedings SMIRT'22 2013 California USA.*

11. Ezzedine El Dandachy M, Briffaut M, Dufour F, Dal Pont S. An original semi-discrete approach to assess gas conductivity of concrete structures. *Int J Numer Anal Methods Geomech.* 2017;41(6):940-956. https://doi.org/10.1002/nag.2655

12. Hillerborg A, Modeer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. *Cem Concr Res.* 1976;6(6):773-782. https://doi.org/10.1016/0008-8846(76)90007-7

13. Bolander JE, Saito S. Fracture analyses using spring networks with random geometry. *Eng Fract Mech.* 1998;61(5-6):569-591. https://doi.org/10.1016/S0013-7944(98)00069-1

14. Cusatis G, Bažant Z, Cedolin L. Confinement-shear lattice model for concrete damage in tension and compression: II. Computation and validation. *J Eng Mech.* 2003;129(12):1449-1458. https://doi.org/10.1061/(ASCE)0733-9399(2003)129:12(1449)

15. Grassl P, Bolander J. Three-dimensional network model for coupling of fracture and mass transport in quasi-brittle geomaterials. *Materials.* 2016;9(9):782. https://doi.org/10.3390/ma9090782

16. Olivier J, Huespe AE, Sánchez PJ. Comparative study on finite elements for capturing strong discontinuities: E-FEM vs X-FEM. *Comput Methods Appl Mech Eng.* 2006;195(37-40):4732-4752. https://doi.org/10.1016/j.cma.2005.09.020

17. Larsson R, Steinman P, Runesson K. Finite element embedded localization band for finite strain plasticity based regularized strong discontinuity. *Mech Cohesive-Frict Mater.* 1998;4(2):171-194. https://doi.org/10.1002/(SICI)1099-1484(199903)4:2<171::AID-CFM81>3.0.CO;2-X

18. Simone A, Wells GN, Sluys LJ. From continuous to discontinuous failure in a gradient enhanced continuum damage model. *Comput Methods Appl Mech Eng.* 2003;192(41-42):4581-4607. https://doi.org/10.1016/S0045-7825(03)00428-6

19. Comi C, Mariani S, Negri M, Perego U. A one dimensional variational formulation for quasibrittle fracture. *Mech Mater Struct.* 2006;1(8):1323-1342. https://doi.org/10.2140/jomms.2006.1.1323

20. Dufour F, Pijaudier-Cabot G, Choinska M, Huerta A. Extraction of a crack opening from a continuous approach using regularized damage models. *Comput Concr.* 2008;5(4):375-388. https://doi.org/10.12989/cac.2008.5.4.375

21. Peerlings RHJ, Geers MGD, de Brost R, Brekelmans WAM. A critical comparison of nonlocal and gradient enhanced softening continua. *J Solids Struct.* 2001;38(44-45):7723-7746. https://doi.org/10.1016/S0020-7683(01)00087-7

22. Mazars J, Pijaudier-Cabot G. Continuum damage theory: application to concrete. *J Eng Mech.* 1989;115(2):345-365. https://doi.org/10.1061/(ASCE)0733-9399(1989)115:2(345)

23. Bažant ZP, Pijaudier-Cabot G. Measurements of characteristic length of nonlocal continuum. *J Eng Mech.* 1989;115(4):755-767. https://doi.org/10.1061/(ASCE)0733-9399(1989)115:4(755)

24. Hillerborg A. Application of fracture mechanics to concrete: summary of a series of lectures. Report TVBM; Vol. 3030 Division of Building Materials, LTH, Lund University 1988 Link: http://portal.research.lu.se/portal/en/publications/application-of-fracture-mechanics-to-concrete--summary-of-a-series-of-lectures-1988(26e32ce6-6cc9-44a0-831d-0a32d0ede1ae).html

25. Mattalah M, La Borderie C, Maurel O. A practical method to estimate crack openings in concrete structures. *Int J Numer Anal Methods Geomech.* 2010;34(15):1615-1633. https://doi.org/10.1002/nag.876

26. Bottoni M, Dufour F, Giry C. Topological search for the crack pattern from a continuum mechanical computation. *J Eng Struct.* 2015;99(15):346-359. https://doi.org/10.1016/j.engstruct.2015.05.005
27. Mazars J, Pijaudier-Cabot G. From damage to fracture mechanics and conversely: a combined approach. *J Solids Struct.* 1996;33(20-22):3327-3342. https://doi.org/10.1016/0020-7683(96)00015-7

28. Planas J, Elices M, Guinea GV. Measurement of the fracture energy using three-point tests. Part 2: influence of bulk energy dissipation. *Mater Struct.* 1992;25(5):305-316. https://doi.org/10.1007/BF02472671

29. Bažant ZP. Concrete fracture models: testing and practice. *Eng Fract Mech.* 2002;69(2):165-205. https://doi.org/10.1016/S0013-7944(01)00084-4

30. Bažant ZP, Becq-Giraudon E. Statistical prediction of fracture parameters of concrete and comparison of testing methods. *Cem Concr Res.* 2002;32(4):529-556. https://doi.org/10.1016/S0008-8846(01)00723-2

31. Mazars J, Hamon F, Grange S. A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings. *Mater Struct.* 2015;48(11):3779-3793. https://doi.org/10.1617/s11527-014-0439-8

32. Pijaudier-Cabot G, Bažant ZP. Non-local damage theory. *J Eng Mech.* 1987;113(10):1512-1533. https://doi.org/10.1061/(ASCE)0733-9399(1987)113:10(1512)

33. Giry C, Dufour F, Mazars J. Stress-based nonlocal damage model. *J Solids Struct.* 2011;4(25-26):3431-3443. https://doi.org/10.1016/j.ijsolstr.2011.08.012

34. Bažant ZP, Pijaudier-Cabot G. Nonlocal continuum damage localization instability and convergence. *J Appl Mech.* 1988;55(2):287-293. https://doi.org/10.1115/1.3173674

35. Bary B. Etude du couplage hydraulique-mécanique dans béton endommagé. PhD LMT Cachan France 1996.

36. Gérard B. Contribution des couplages mécaniques-chimie-transfert dans la tenue à long-terme des ouvrages de stockage de déchets radioactifs. PHD ENSC France 1996.

37. Gawin D, Pesavento F, Schrefler B. A simulation of damage-permeability coupling in hydro-thermo-mechanical analysis of concrete at high temperature. *Commun Numer Methods Eng.* 2002;18(2):113-119. https://doi.org/10.1002/cnm.472

38. Dal Pont S, Durand S, Schrefler B. A multiphase thermo-hydro-mechanical model for concrete at high temperatures—finite element implementation and validation under LOCA load. *Nucl Eng Des.* 2007;237(2):2137-2150. https://doi.org/10.1016/j.nucengdes.2007.03.047

39. Jason L, Masson B. Comparison between continuous and localized methods to evaluate the flow rate through containment concrete structures. *Nucl Eng Des.* 2014;277:146-153. https://doi.org/10.1016/j.nucengdes.2014.06.010

40. Pijaudier-Cabot G, La Borderie C. Mechanical damage, chemical damage and permeability in quasi-brittle cementitious materials. *Eur J Environ Civil Eng.* 2009;13(7-8):963-982. https://doi.org/10.1080/19648189.2009.9693163

41. Asali M, Capra B, Mazars J, Colliat JB. Numerical strategy for forecasting the leakage rate of inner containments in double-wall nuclear reactor buildings. *J Adv Concrete Technol.* 2016;14(8):408-420. https://doi.org/10.3151/jact.14.408

42. Ezzedine El Dandachy M. Characterization and modelling of permeability of damaged concrete: application to reinforced concrete structures. PhD thesis UGA Grenoble France 2016.

43. Nguyen D, Lawrence C, La Borderie C, Mattalah M, Nahas G. A mesoscopic model for a better understanding of the transition from diffuse damage to localized damage. *Eur J Environ Civil Eng.* 2010;14(6-7):751-776. https://doi.org/10.1080/19648189.2010.9693261

44. Fichant S, La Borderie C, Pijaudier-Cabot G. Isotropic and anisotropic descriptions of damage in concrete structures. *Mech Cohesive-Frict Mater.* 1999;4(4):339-359. https://doi.org/10.1002/(SICI)1099-1484(199907)4:4<339::AID-CFM65>3.0.CO;2-J

45. Cast3m software on www-cast3m.cea.fr

**How to cite this article:** Bouhjiti DE-M, Ezzedine El Dandachy M, Dufour F, et al. New continuous strain-based description of concrete's damage-permeability coupling. *Int J Numer Anal Methods Geomech.* 2018;1-27. https://doi.org/10.1002/nag.2808

# APPENDIX A

## HILLERBORG REGULARIZATION OF LOC/WD/SD DAMAGE MODELS

The application of Equation 10 to different softening laws allows the identification of their remaining unknown parameters (Table 1).7

### A.1 | LOC/SD/WD models

For LOC/SD/WD models, the weighting function does not depend on the stress variable. This leads to a unified description of their associated Hillerborg regularization equation:
$$ \left\{ \begin{array}{l} \frac{G_F}{D_c E} = \frac{\varepsilon_{d0}^2}{2} + \int_{\varepsilon_{d0}}^{+\infty} (1-d)\varepsilon \, d\varepsilon \\ \sigma \le f_t \quad \forall \varepsilon \end{array} \right. \qquad (A1) $$

* Linear softening law: In Equation 3 (LIN), the unknowns are the damage threshold $\varepsilon_{d0}$ and the failure strain $\varepsilon_u \ge \varepsilon_{d0}$ (Figure 1A). After simple calculations (Equation 10), one can derive

$$ D_c \le 2 \frac{G_F E}{f_t^2} \left\{ \begin{array}{l} \varepsilon_{d0} = \frac{f_t}{E} \\ \varepsilon_u = 2 \frac{G_F}{f_t D_c} \end{array} \right. \qquad (A2) $$

For $D_c = 2 \frac{G_F E}{f_t^2}$, the computed behaviour is perfectly brittle: $\varepsilon_u = \varepsilon_{d0}$ and the fracture energy per an element sized $D_c$ is equal to the elastic one.

* Bi-linear softening law: In Equation 4 (BIL), the unknowns are the damage threshold $\varepsilon_{d0}$, the failure strains $\varepsilon_u'$ and $\varepsilon_u$, and the n parameter (Figure 1B). The number of parameters exceeds the number of available equations. For that reason, additional hypotheses need to be considered$^{29,30}$. The first slope is defined using the intrinsic fracture energy that represents 40% of the total fracture energy. Eventually, the second slope is deduced to verify the global Hillerborg criterion in Equation 10.

$$ D_c \le \frac{4}{5} \frac{G_F E}{f_t^2} \left\{ \begin{array}{l} \varepsilon_{d0} = \frac{f_t}{E} \\ \varepsilon_u' = \frac{4}{5} \frac{G_F}{f_t D_c} \\ \varepsilon_u = \frac{4 + 6n}{5} \frac{G_F}{f_t D_c} \end{array} \right. \qquad (A3) $$

Compared to the linear softening law, values of $D_c$ are more restrained and cannot exceed $\frac{4}{5} \frac{G_F E}{f_t^2}$ given the use of the intrinsic energy ($\frac{2}{5} G_F$) rather than the total one $G_F$.

* Exponential softening law: In Equation 5 (EXP), the unknowns are the damage threshold $\varepsilon_{d0}$ and the shape parameter $B_t$ (Figure 1C). Calculations are more complex than the 2 previous ones but eventually lead to the following system of equations:

$$ D_{c, \min} < D_c \le 2 \frac{G_F E}{f_t^2} \left\{ \begin{array}{l} B_t = \frac{\varepsilon_{d0} + \sqrt{\frac{4G_F}{E D_c} - \varepsilon_{d0}^2}}{2 \left( \frac{G_F}{E D_c} - \frac{\varepsilon_{d0}^2}{2} \right)} \\ D_{c, \min} < D_c \le \frac{2}{5} \frac{G_F E}{f_t^2} : \varepsilon_{d0} = \frac{\ln\left(\frac{f_t}{E}B_t\right) + 1}{B_t} \le \frac{f_t}{E} \\ \frac{2}{5} \frac{G_F E}{f_t^2} \le D_c \le 2: \quad \varepsilon_{d0} = \frac{f_t}{E} \end{array} \right. \qquad (A4) $$

To prevent the damage threshold $\frac{f_t}{E}$ from being inferior to $\frac{f_t}{E}$, another regularization strategy is suggested in Mazars et al.$^{31}$ It consists of considering the linearized postpeak behaviour law (EXPL—tangent line at $\varepsilon = \frac{f_t}{E}$—Equation A5) instead of the exponential one (Figure A1).
<figure><img src="image_13.png" /><figcaption>FIGURE A1 (A) Total (Equation A4) vs (B) linearized (Equation A6) energy-based regularization applied to the exponential softening law [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

This simplifies calculations (Equation A6) since by definition, the yield strength becomes equal to the ultimate tensile strength of concrete. However, it does not reduce the number of unknowns. Indeed, even if the damage threshold is known, the amount of the fracture energy needed for regularization remains to be identified. This amount, denoted $\alpha G_F$, should allow the area under the behaviour law to approaches the best $G_F/D_c$.

$$ \varepsilon \geq \varepsilon_{d0}:\sigma = E(1 - B_t\varepsilon_{d0}) (\varepsilon - \varepsilon_{d0}) + E\varepsilon_{d0}, \quad (A5) $$

$$ 0 \le D_c \le \min\left(\frac{1}{5}; \alpha\right) 2 \frac{G_F E}{f_t^2} \begin{cases} \varepsilon_{d0} = \frac{f_t}{E} \\ B_t = \frac{E}{f_t} \left( 1 + \frac{D_c f_t^2}{2\alpha G_F E - D_c f_t^2} \right) \\ G_F^{\text{EXPL}} = \frac{D_c f_t^2}{E} \left( \frac{1}{2} + \frac{1}{1 + \frac{D_c f_t^2}{2\alpha G_F E - D_c f_t^2}} + \frac{1}{1 + \frac{D_c f_t^2}{2\alpha G_F E - D_c f_t^2}} \right)^2 \le G_F \end{cases} \quad (A6) $$

where $G_F^{\text{EXPL}}$ stands for the dissipated energy using the linearized law regularization.

A simple functional analysis of the condition $G_F^{\text{EXPL}} \le G_F$ for $D_c = \frac{2 G_F E}{5 f_t^2}$ leads to the conclusion that $G_F^{\text{EXPL}} \to G_F$ only when $\alpha \to +\infty$. Hence, the computed fracture energy will always be underestimated compared to the desired one leading to the computation of a more brittle behaviour. Eventually, the choice of $\alpha = 1$ is the only one acceptable with respect to $D_c \le 2 \frac{G_F E}{f_t^2}$. The regularization of the exponential law is controversial for small mesh sizes ($D_c \le \frac{2G_F E}{5 f_t^2}$). However, more errors are induced when using the linearized regularization technique recommended in Mazars et al³¹ leading to the computation of a largely overestimated brittleness. So the total energy regularization should be retained instead, keeping in mind that the damage initiation might occur before reaching the threshold $\frac{f_t}{E}$ for small influence domains and that the crack opening is somewhat overestimated after initiation (in comparison with linear and bilinear softening laws). One might also refer to the postpeak softening exponential law in Fichant et al,⁴⁴ which has the same validity domain as the linear one and does not show any complications in terms of energy regularization.
<table><caption>TABLE A1 Identification of the LIN/BIL/EXP softening laws parameters using the Hillerborg energy based</caption><thead><tr><th rowspan="2">Formulation<br>Influence Domain D<sub>c</sub><br>Spatial Discretization h</th><th>LOC</th><th>LOC</th><th>LOC</th><th>LOC</th><th>LOC</th><th colspan="2">SD/WD</th><th colspan="2">SSD/SWD</th></tr><tr><th>h</th><th></th><th></th><th></th><th></th><th colspan="2">l<sub>FPZ</sub></th><th colspan="2"></th></tr><tr><th></th><th>D<sub>c, min</sub> = L/60</th><th>L/50 ∈ [D<sub>c, min</sub>; 2 G<sub>F</sub>E/5 f<sub>t</sub><sup>2</sup>]</th><th>2 G<sub>F</sub>E/5 f<sub>t</sub><sup>2</sup></th><th>l<sub>FPZ</sub></th><th>4 G<sub>F</sub>E/5 f<sub>t</sub><sup>2</sup></th><th>L/60</th><th>L/16</th><th>L/60</th><th>L/16</th></tr></thead><tbody><tr><td rowspan="2">LIN</td><td>ε<sub>d0</sub>, µm/m</td><td>150</td><td></td><td></td><td></td><td></td><td>592</td><td></td><td>907</td><td></td></tr><tr><td>ε<sub>u</sub>, µm/m</td><td>2200</td><td>1800</td><td>756</td><td>592</td><td>355</td><td>592</td><td></td><td></td><td></td></tr><tr><td rowspan="3">BIL</td><td>ε<sub>d0</sub>, µm/m</td><td>150</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ε'<sub>u</sub>, µm/m</td><td>888</td><td>711</td><td>302</td><td>237</td><td>142</td><td>237</td><td></td><td>593</td><td></td></tr><tr><td>ε<sub>u</sub>, µm/m</td><td>4900</td><td>3900</td><td>1700</td><td>1300</td><td>782</td><td>1300</td><td></td><td>3260</td><td></td></tr><tr><td rowspan="2">EXP</td><td>ε<sub>d0</sub>, µm/m</td><td>2.4</td><td>71</td><td>150</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>B<sub>t</sub></td><td>2467</td><td>3053</td><td>6619</td><td>8194</td><td>14274</td><td>8194</td><td></td><td>5781</td><td></td></tr></tbody></table>

## A.2 | SSD/SWD models

On the other hand, SSD and SWD models' regularization should in principle account for the constant decrease of their associated effective internal length due to damage localization ($D'_{c} = l_{FPZ} \frac{\sigma}{f_t}$ → 0 and $\phi \to \delta$ when $d \to 1$ and $\sigma \to 0$). As a consequence, the identified parameters are bounded by those of a local modal ($D_c = h$) and those of a strain-based nonlocal model ($D_c = l_{FPZ}$). As a result, the solving of Equation 10 is not obvious due to the presence of nonelementary functions. If SSD and SWD are regularized according to LOC models ($D_c = h < l_{FPZ}$), then the energy available for dissipation is overestimated whereas the regularization according to WD and SD models leads to its underestimation (Table A1). Eventually, in the case of SSD and SWD formulations, a numerical solving of Equation 10 is required.

## A.3 | Discussion

One can notice that the use of LIN/BIL and EXP softening laws is conditioned by the size of the influence domain to meet energetic criteria. Particularly, the use of LIN and EXP laws seem the least constraint by the influence domain $D_c$ compared with the BIL softening that requires additional hypothesis (estimation of the intrinsic fracture energy). One should pay attention to the use of the EXP law for refined mesh (in the case of local models) or small FPZ sizes (in the case of nonlocal modelling). For $D_c \ge \frac{2E G_F}{5 f_t^2}$, the computed behaviour is quasibrITTLE to brittle as $D_c$ increases. The maximum value for which Equation 10 remains verified is $D_c = 2 \frac{E G_F}{f_t^2}$ for a perfectly brittle law. However, for $D_c \le \frac{2E G_F}{5 f_t^2}$, the system of equations is strongly non-linear and requires numerical solving. Moreover, in order for the solution to be acceptable, the damage threshold should remain strictly positive, hence the existence of a minimal value $D_{c, min}$. In opposition to the LIN and BIL softening laws, the computed behaviour for EXP law shows a hardening phase for small $D_c$ values: To dissipate the same fracture energy and not overestimate the tensile strength, the damage threshold needs to be reduced, which induces a hardening behaviour until the tensile strength is reached (Figure B2A,B). The main drawback, though, is the underestimation of the damage threshold (at least from a numerical point of view) $\epsilon_{d0} \le \frac{f_t}{E}$ (Table A1) and overestimation of the failure strain at low damage values. This result is controversial, though, from a physical point of view since the damage strain threshold is physically constant and cannot be equal to zero nor depend on the mesh size (for local models). This restrains even more the use of the EXP law when both energetic and damage strain threshold are considered.

# APPENDIX B

## EFFECT OF LOC/SD/WD/SSD/SWD FORMULATIONS ON THE GLOBAL AND LOCAL BEHAVIOUR

To illustrate the utility of the regularization steps presented above, a 1D beam under pure tensile loads is considered. It has a length of $L = 1$ m, is fixed on one side, and undergoes a tensile load on the other (Figure B1). To prevent a diffuse
<figure><img src="image_14.png" /><figcaption>**FIGURE B1** 1D beam under tensile test with a predefined position of default for damage localization</figcaption></figure>

damage state and be able to characterize strain and damage distribution around a localized crack, a centred default is considered.

The selected weak point has a tensile strength that is reduced by 1% compared to the rest of the beam. The following mechanical properties are considered: $f_t = 4.5$ MPa, $G_f = 80$ N/m, and $E = 30$ GPa. The FPZ width is supposed to be $D_c = l_{FPZ} \approx 3 * d_{max} = 6$ cm ($d_{max} = 2$ cm).

As for the spatial discretization, values of $h$ are issued from Table 4 defining the validity domain of each postpeak law and local/nonlocal damage modelling: $h = \frac{L}{60} = 1.6$ cm $= D_{c, \min} < l_{FPZ}$, $h = \frac{L}{50} = 2$ cm $< l_{FPZ}$, $h = \frac{L}{21} = 4.7$ cm $\approx \frac{2}{5} \frac{G_f E}{f_t^2} < l_{FPZ}$, $h = \frac{L}{16} = 6$ cm $= l_{FPZ}$, and $h = \frac{L}{10} = 10$ cm $\approx \frac{4}{5} \frac{G_f E}{f_t^2} > l_{FPZ}$. The identified parameters for LIN, BIL, and EXP laws are summarized in Table A1.

The implemented behaviour laws are depicted in Figure B2. One can notice, in the case of local models, how the computed brittleness increases with the size of the considered FE. Figure B2A,B shows the hardening behaviour of the EXP softening law for small mesh sizes. In the case of stress-based nonlocal models, one can notice how the area under the $(\sigma, \varepsilon)$ curve is higher than the one of strain-based nonlocal models given the decrease of their internal length with stress concentration at the crack tip (Table A1 and Figure B2D,F). Nevertheless, the dissipated energy around the FPZ remains the same.

## B.1 | Local response

By local response, it is meant the damage distribution and energy dissipation over the influence domain $D_c$ (FPZ), and crack opening evolution at the weak point (WP) in Figure B1. Their estimation can be performed using explicit and analytical developments. In Figure B3, the maximum damage value at the weak point $d_{WP}$ is depicted in the case of local and nonlocal formulations, whereas in Figure B4, damage distribution over the influence domain $d_{FPZ}$ is shown for nonlocal configurations (the damage variable is homogeneous over the mesh size $h$ in the case of local models). In Figure B5, the crack strain evolutions at the weak point are provided. At this position, one can see the dependence of the damage and crack opening evolutions on both the used softening law (Figure 1) and damage modelling strategy (local Figure B5A-E and nonlocal Figure B5D,F/Figure B4). In terms of the maximum damage value (at the weakest point—Figure B3), more gaps between local and nonlocal approaches are observed when the mesh size $h$ recedes from the FPZ size $l_{FPZ}$. This result underlines again the previous remark mentioned in Section 2.1.3 for crack opening estimation that local damage models are only comparable with nonlocal ones when the considered spatial discretization is equal to the FPZ size (according to the hypothesis of one crack per FPZ width). For $D_c = h = l_{FPZ} = 6$ cm, absolute differences between the damage evolutions at the weakest point (Figures B3D,F and B4) due to the use of different softening laws and different weighting functions remain less than 0.2. As for the equivalent crack opening values, absolute differences due to the shape of the softening laws are less than 3 µm (Figure B5D,F). Within these margins, all damage modelling strategies and softening laws can be considered equivalent in terms of stress–strain laws (Figure B4A) and energy dissipation over the FPZ (Figure B6). But even for the same dissipated fracture energy $G_F^{FPZ}$, it is worth mentioning that nonlocal modelling of damage around the FPZ is more realistic (particularly stress-based approach), as it decreases away from the localization point (Figure B4A), compared with local models where damage is constant over the influence domain.

## B.2 | Global response

As one can see in Figure B2C-F, the LIN and EXP behaviour laws are about the same within their respective validity domains in Table 4. One should expect similar response at the structural scale as well. As for the BIL law, more gaps are expected given the different energy dissipation paths. Herein, the study is limited to the EXP law and LOC, WD, and SWD damage formulations (similar tendencies are expected for LOC, SSD, and SSD). The obtained structural responses are shown in Figure B7A-C. They are expressed in terms of the normalized global stress (force over the beam's section) vs the displacement over the beam's length. Figure B7D-F shows the postprocessed crack opening as a function of the imposed displacement. The numerical problem is solved using Cast3m FE code$^{45}$ in which an arc-length
<figure><img src="image_15.png" /><figcaption>**FIGURE B2** Behaviour law using LIN/BIL/EXP softening laws with respect to the Hillerborg energy criterion. For local models: (A) $h = L/60$, (B) $h = L/50$, (C) $h = L/21$, (D) $h = L/16$, and (E) $h = L/10$. For nonlocal models with $l_{FPZ} = 6$ cm: (D) SD/SW and (F) SSD/SSW</figcaption></figure>

technique is implemented to trace the equilibrium path when the global response exhibits eventual snap backs. From Figure B7A,D, one can observe that the stress-based nonlocal formulations SWD is more sensitive to the spatial discretization compared with strain-based ones WD; the more refined the mesh is, the more brittle the global structural behaviour gets. This result is expected given that the performed regularization (Table 3) is performed based on a continuous solving where the influence domain tends to zero as damage localizes. In FE analysis such configuration is never reached since, by definition, the influence domain cannot be less than the mesh characteristic length $h$. Consequently, the computed brittleness is underestimated as the mesh size gets bigger. In other words, the use of refined mesh in the case of stress-based models is recommended in line with the regularization in Table 3 to obtain realistic responses. Generally and for the same mesh size $h$, the SWD result remains bounded by LOC and WD responses (Figure B7B,E) as previously underlined. As for LOC models (Figure B7C,F), the hardening behaviour of the EXP law obtained for $h<\frac{2G_F E}{5 f_t^2}$ induces a mesh-dependent global response (such results are not obtained for LIN and BIL laws, which remain brittle to quasibrITTLE). However, for $h>\frac{2 G_F E}{5 f_t^2}$ and $\epsilon_{d0} = f_t/E$, the Hillerborg regularization ensures the objectiveness of the global response.

## B.3 | Discussion

As it has been explicitly shown in this section, the use of continuous damage-based models is able to describe concrete cracking whether local or nonlocal approaches are selected. However, the obtained results can largely differ in the absence of additional energetic and usage restraints. With that regard, one can define 3 criteria according to which a given modelling strategy can be selected or discarded:

* Local response objectivity: It refers to the strain, stress, damage, and crack opening estimation within the considered domain of influence. At this scale, and in terms of local quantities' (particularly strain and damage) distribution away from the localization point, nonlocal models are more physically representative compared with local ones. However, the stress-based damage models are the only ones able to reproduce the decrease of the influence domain ($D_c'_{1D} = l_{FPZ} (\sigma/f_t)$ in Table 2) due to the stress unloading in the vicinity of crack lips.³³ Such differences are expected to
<figure><img src="image_16.png" /><figcaption>**FIGURE B3** Damage variable at the weakest point using LIN/BIL/EXP softening laws with respect to the Hillerborg energy criterion. For local models: (A) $h = \frac{L}{60}$, (B) $h = \frac{L}{50}$, (C) $h = \frac{L}{21}$, (D) $h = \frac{L}{16}$, and (E) $h = \frac{L}{10}$. For nonlocal with $l_{FPZ} = 6$ cm: (D) SD/SW and (F) SSD/SSW</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**FIGURE B4** Effects of the weighting functions and softening laws on the computed damage profiles (A, LIN-SD; B, BIL-SD; C, EXP-SD; D, LIN-WD; E, BIL-WD; F, EXP-WD; G, LIN-SSD; H, BIL-SSD; I, EXP-SSD; J, LIN-SWD; K, BIL-SWD; L, EXP-SWD)</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIGURE B5 Hillerborg energy-based regularization and the computed crack opening using LIN/BIL/EXP softening laws. For local models: (A) $h = L/60$, (B) $h = L/50$, (C) $h = L/21$, (D) $h = L/16$, and (E) $h = L/10$. For nonlocal models with $l^{FPZ} = 6$ cm: (D) SD/SW and (F) SSD/SSW</figcaption></figure>

affect the modelling of permeability as well especially if it is damage based. When the mean mechanical behaviour over the influence domain is considered (represented by the averaged crack opening or mean damage and strains), and regardless of the used postpeak law, local models show approximately the same response as nonlocal ones when the mesh size is equal to the FPZ size. In other words, local and nonlocal models are locally comparable when they both (1) have the same influence domains and (2) verify the same energetic criterion. This is true also in terms of the cracks' number as only one crack can develop within a given influence domain; physically, the influence domain of a crack has the same size as the FPZ over which dissipation occurs.

* Semiglobal response objectivity: By semiglobal is meant the comparison of local and global quantities as usual experiments provide a partial description of concrete behaviour at the local scale (FPZ). For instance, in the case of splitting or flexion tests, the maximal crack opening (local quantity) is monitored along with the applied force (global

<figure><img src="image_19.png" /><figcaption>FIGURE B6 Effects of the weighting functions and softening laws on the computed fracture energy (A) LIN-SD/WD/SSD/SWD (B) BIL-SD/WD/SSD/SWD (C) EXP-SD/WD/SSD/SWD</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**FIGURE B7** Effect of spatial discretization on the computed global behaviour. Normalized global stress vs imposed displacement: (A) nonlocal $h \le l_{FPZ}$, (B) local and nonlocal $h = l_{FPZ}$, and (C) local. Crack opening at the weak point vs imposed displacement: (D) nonlocal $h \le l_{FPZ}$, (E) local and nonlocal $h = l_{FPZ}$, and (F) local</figcaption></figure>

quantity) whereas the stress state around the FPZ is (experimentally) uncharacterized. Since structural effects are considered, the objectiveness of the semiglobal response is conditioned by the global one and vice versa. Obviously, such aim can only be achieved for the mean mechanical response when local and nonlocal models have the same domain of influence (local objectiveness).

* Global response objectivity: In the case of local models, the objectiveness of the global (or structural) response is ensured when the Hillerborg energy-based and damage threshold criteria are applied. Within the validity domain of each postpeak law, the global response scarcely varies with the mesh size. For nonlocal models, the stress-based ones seem particularly mesh dependent as their influence domain tends to localize within a defined mesh size even if, locally at the FPZ level, the dissipated energy is constant. Their mesh dependency seems to decrease as the spatial discretization gets smaller, which can become heavily time consuming for large structures. Generally speaking, at the structural scale, the local, nonlocal, and stress-based nonlocal models can hardly provide the same response even though the energy criterion is applied. Nevertheless, when local objectiveness is ensured at least for the mean behaviour of the FPZ, LOC and WD/SD remain distant from the SWD/SSD configurations that show naturally more brittleness at the structural scale.