# Estimation of crack opening from a two-dimensional continuum-based finite element computation

Frédéric Dufour, Grégory Legrain, Gilles Pijaudier-Cabot, Antonio Huerta

► To cite this version:

Frédéric Dufour, Grégory Legrain, Gilles Pijaudier-Cabot, Antonio Huerta. Estimation of crack opening from a two-dimensional continuum-based finite element computation. International Journal for Numerical and Analytical Methods in Geomechanics, 2012, 36 (16), pp.1813-1830. 10.1002/nag.1097 . hal-00854615

HAL Id: hal-00854615

https://hal.science/hal-00854615

Submitted on 3 Jul 2017

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Estimation of crack opening from a two-dimensional continuum-based
finite element computation

Frédéric Dufour¹, Grégory Legrain², Gilles Pijaudier-Cabot³ and Antonio Huerta⁴

¹Laboratoire 3S-R Grenoble-INP/Université Joseph Fourier/CNRS UMR 5521, France

²GeM-Institut de Recherche en génie civil et mécanique, Ecole Centrale de Nantes, France

³Laboratoire des Fluides Complexes, Université de Pau et des Pays de l'Adour, France

⁴Laboratori de Càlcul Numèric, Universitat Politècnica de Catalunya, Barcelona, Spain

Damage models are capable of representing crack initiation and mimicking crack propagation within a
continuum framework. Thus, in principle, they do not describe crack openings. In durability analyses
of concrete structures however, transfer properties are a key issue controlled by crack propagation and
crack opening. We extend here a one-dimensional approach for estimating a crack opening from a
continuum-based finite element calculation to two-dimensional cases. The technique operates in the case
of mode I cracking described in a continuum setting by a nonlocal isotropic damage model. We used the
global tracking method to compute the idealized crack location as a post-treatment procedure. The orig-
inal one-dimensional problem devised in Dufour et al. [4] is recovered as profiles of deformation orthog-
onal to the idealized crack direction are computed. An estimate of the crack opening and an error indicator
are computed by comparing finite element deformation profiles and theoretical profiles corresponding to a
displacement discontinuity. Two estimates have been considered: In the strong approach, the maxima of
the profiles are assumed to be equal; in the weak approach, the integrals of each profile are set equal.
Two-dimensional numerical calculations show that the weak estimates perform better than do the strong
ones. Error indicators, defined as the distance between the numerical and theoretical profiles, are less than
a few percentages. In the case of a three-point bending, test results are in good agreement with experimental
data, with an error lower than 10% for widely opened crack (> 40µm).

KEY WORDS: Crack opening; Crack location; Damage mechanics

# 1. INTRODUCTION

For many concrete structures, crack propagation and crack opening are key parameters needed to estimate durability. Cracks are preferential paths along which fluids or corrosive chemical species may penetrate inside structural elements exposed to aggressive environment. For concrete structures such as confinement vessels, reservoirs, or nuclear waste disposals, tightness to gas or liquids is also a major serviceability criterion governed by mass transfer characteristic quantities such as the permeability in a continuum description. Permeability is directly related to damage in concrete (see e.g., [2]), and it jumps several orders of magnitude upon macrocracking [29,7]. According to Poiseuille's law, the permeability of a cracked structure (with a single crack) is proportional to the square of the crack opening. Hence, quantifying crack openings in concrete structure upon failure analyses is an important issue.
Fictitious crack models are among the most popular approaches for this purpose. They are based on an explicit description of the discontinuity within the material (e.g., the cohesive crack model [8]). Crack opening is related to the stress level according to a relationship based on fracture mechanics. Cohesive crack models need proper algorithms for crack propagation, but more importantly, they are not capable of describing crack initiation. Enhanced continuum models for failure analyses (e.g., based on damage or on plasticity) are capable of representing diffuse cracking, macrocrack initiation, and possibly crack propagation [21,19]. Damage models for instance regard cracking as an ultimate consequence of a gradual loss of material integrity. These models, however, do not predict crack opening as they rely on a continuum approach to fracture. They need to be complemented by descriptions of displacement discontinuities and some bridges between continuous and discontinuous descriptions of failure.

Ideally, failure analyses that involve inception of cracking, crack propagation, and calculation of crack opening would require to merge the continuum damage approach and the discrete crack approach into a single, consistent, computational model. Bridges between damage and fracture have been devised in the literature (see e.g., [13,24]). They rely on the equivalence between the dissipation of energy due to damage and the energy dissipated to propagate a crack. The strong discontinuity approach initiated by [26] and widely used over the last decade (e.g., [18,10]) offers also the possibility of merging within the same formulation a continuous damage model for the bulk response and a cohesive model for the discontinuous part of the kinematics. It is a combination of continuum-discrete modeling that is appealing from the point of view of the physics of fracture. The pending issue in combining the continuum-based model for crack initiation and then a discrete crack model for propagation is, however, the threshold upon which the later prevails over the former description. Usually, it is considered that the discontinuity appears when damage, stresses, or strain energy reach a certain threshold fixed beforehand [3,27]. Consequently, estimates of crack opening may strongly depend upon this threshold, which is, in most cases, fixed rather arbitrarily.

The present article presents a technique aimed at computing the crack path and the crack opening directly from the continuum-based approach. It is an extension of the method proposed by Dufour and coworkers [4] to two-dimensional problems. Such an extension requires specific procedures for the detection of cracks and crack path, which were not considered in the prototype one-dimensional version of Dufour et al. [4]. These techniques are first discussed in the present article. Once the crack path is located, the displacement jump is derived by comparing the profiles of the regularized projection of the finite element (FE) deformation, perpendicular to the crack path, against theoretical profiles corresponding to a displacement discontinuity. These profiles are compared along straight lines perpendicular to the crack path. The proposed explicit post-processing calculation provides on top some indication on the accuracy of the calculation of the crack opening. In principle, this additional information could be used for switching from a continuum-based to a discontinuous-based formulation (e.g., according to the strong discontinuity approach).

Only mode I crack propagation is considered as it is the most common source of cracking in concrete (this encompasses also the cases of curved cracks if the principal strain direction rotate). The constitutive equations are based on the nonlocal version of the scalar damage model of Mazars and Pijaudier-Cabot [12], which describes fracture and damage growth due to extension. Basically, the detection of the crack path relies on the capability of the constitutive equations to capture the propagation of the same type of cracks. Once a crack has been detected however (and provided it does not propagate), the present estimate of the crack opening may still be used as it is computed from local extensions perpendicular to the crack path. There are no specific requirements on the crack path in this case.

First, we recall the continuum approach that will be considered, the (integral) nonlocal damage model. Then, the detection of the crack location in the computational domain and the estimate of its opening are discussed. Finally, numerical examples and comparisons with experimental data are discussed.

## 2. CONSTITUTIVE MODELING: NONLOCAL DAMAGE MODEL

The scalar isotropic damage model of Mazars and Pijaudier-Cabot [12] will be used in the FE computations for representing the progressive failure. This constitutive relation in its local version exhibits
strain softening, and a regularization technique shall be considered to avoid mesh dependency and ill-posedness of the governing equations of equilibrium. Nonlocal models are known to possess shortcomings such as spurious boundary effects on fracture propagation [9,22] or incorrect initiation of damage at a crack tip [28]. Still, these drawbacks in the model formulation do not alter their ability to capture a fully localized, mode I crack. Note also that the proposed estimate of the crack opening is independent from the details of the material model used for failure analysis. Only the variables needed to capture the crack location are directly related to the constitutive model.

In this model, the tensorial stress $\sigma$–strain $\varepsilon$ relationship is expressed as follows:

$$ \sigma = (1 - D)C : \varepsilon \quad (1) $$

where D is the damage scalar variable and C is the elastic stiffness tensor of the sound material. Damage is a combination of two components: $D_t$, which is damage due to tension based loads, and $D_c$ which is damage due to compression:

$$ D = \alpha_t D_t + \alpha_c D_c \quad (2) $$

$\alpha_t$ and $\alpha_c$ depend on both strain and stress tensors. Damage evolution laws for both traction and compression components read

$$ D_{t,c} = 1 - \frac{Y_{D0}(1 - A_{t,c})}{\bar{Y}} - \frac{A_{t,c}}{e^{(B_{t,c}(\bar{Y} - Y_{D0}))}} \quad (3) $$

where $A_t$, $A_c$, $B_t$, $B_c$, and $Y_{D0}$ are model parameters, and $\bar{Y}$ is defined by

$$ \bar{Y} = \max(\bar{Y}, \bar{\varepsilon}_{eq}) \quad (4) $$

with $\bar{Y} = Y_{D0}$ initially. The nonlocal equivalent strain [21] is defined as a weighted average of the local equivalent strain $\varepsilon_{eq}$:

$$ \bar{\varepsilon}_{eq}(x) = \frac{\int_{\Omega} \phi(\mathbf{x} - \mathbf{s}) \varepsilon_{eq}(\mathbf{s}) ds}{\int_{\Omega} \phi(\mathbf{x} - \mathbf{s}) ds} \quad (5) $$

The local equivalent strain follows Mazars’ definition:

$$ \varepsilon_{eq} = \sqrt{\sum_{i=1}^{3} < \varepsilon_i >_+^2} \quad (6) $$

where $< >_+$ denotes the positive part of the principal strain $\varepsilon_i$. Several weight functions $\phi(\mathbf{x} - \mathbf{s})$ exist in the literature; we chose the most often used, i.e. the Gaussian function

$$ \phi(\mathbf{x} - \mathbf{s}) = \exp\left(-\left(\frac{2\|\mathbf{x} - \mathbf{s}\|}{l_c}\right)^2\right) \quad (7) $$

where $l_c$ is the internal length of the model.

### 3. ESTIMATION OF CRACK OPENING

In the study of Dufour et al. [4], a technique was presented to estimate an equivalent crack opening from a nonlocal damage computation. The crack location was known a priori, and the computational domain was reduced to one dimensional. To extend this approach to more general cases, the location of the crack needs to be determined first. In the present article, we will deal with an “idealized crack”, that
is, a line (or a surface in the three-dimensional case) of maximum damage. On this line, damage is equal, or very close, to 1 and decreases toward 0 within the fracture process zone. Due to the assumptions in the constitutive model, the idealized crack is the result of damage growth due to extensions. Hence, the detection of the idealized crack path should be considered as a restriction to mode I crack propagation, which is not, however, a very important restriction for concrete structures. Equivalently, we will consider that the idealized crack is a line perpendicular to the direction of the maximum principal strain.

### 3.1. Localization of an idealized crack

Several techniques have been proposed to locate cracks and to update their length during their propagation. In Comi et al. [3], the authors proposed to fit a fourth-order polynomial on the damage field and then to propagate the crack in a direction that is perpendicular to the maximum curvature of the polynomial at the crack tip. The main drawback of this approach is that the accuracy may be limited if the curvature of the damage field is large. Moreover, the accuracy of the fitting may not be sufficient at the crack tip, which should lead to extra difficulties for the estimation of the crack direction. In the study of Mariani and Perego [11], this procedure was enhanced. It was applied on the stress field in a half disc centered at the crack tip to avoid the previously mentioned drawbacks. Alternatively, averaging approaches have been proposed in Refs. [14] and [30]. Later on, the directions pointing from the tip to the most damaged points at various radii are averaged on a V-shaped window.

Another technique called global tracking algorithm has been proposed by Oliver and Huespe [16,17]. It was used first within the strong discontinuity approach to evaluate the direction of crack propagation: The resolution of a heat-conduction-like problem leads to a scalar function whose iso-values represent all the possible directions of propagation. The selection of the iso-value emanating from the crack tip provides the direction of propagation of the crack. The technique has been modified by Feist and Hofstetter [6] to restrict the heat conduction problem on a subset of elements already or potentially crossed by the crack. It was also applied in the context of the eXtended Finite Element Method [15] by Dumstorff and Meschke [5], who compared various crack branching criteria. Here, we will use this technique to locate the crack path.

#### 3.1.1. Global tracking algorithm.
The determination of the crack location is a separate problem solved as a post-processing phase, after the mechanical calculation. Let us assume that there exists a vector $\vec{T}(x)$ (with unit norm), which is tangent to the crack defined as a line (or a surface in the three-dimensional case). The construction of the envelopes of $\vec{T}(x)$ supplies all the possible surfaces or lines at time $t$ (see Figure 1).

The envelopes of $\vec{T}(x)$ are described by function $\theta(x)$ whose level contours ($\theta(x) = \text{constant}$) define all the possible surfaces/lines where the equivalent strain is maximum, as described in Figure 1. The gradient of this function must be normal to $\vec{T}(x)$ in each point:

$$ \vec{T}(x) \cdot \text{Grad}\theta = \frac{\partial\theta}{\partial\vec{T}} = 0 \qquad (8) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. Global tracking algorithm: envelopes of the vector field $\vec{T}(x)$, possible crack path, and real crack.</figcaption></figure>
This condition can be reformulated as the following linear boundary value problem [17]:

$$ \mathrm{div}(\mathbf{K} \mathrm{Grad} \theta) = 0 \text{ in } \Omega \qquad (9) $$

where $\Omega$ is the domain occupied by the solid (or any subdomain containing the crack) and $\mathbf{K}$ is the second-order tensor defined as

$$ \mathbf{K}(\mathbf{x}) = \vec{T}(\mathbf{x}) \otimes \vec{T}(\mathbf{x}) \qquad (10) $$

Equation (9) is subjected to appropriate boundary conditions [17]:

$$ (-\mathbf{K} \mathrm{Grad} \theta) \cdot \vec{n} = 0 \quad \text{on } \partial_q \Omega \qquad (11) $$

$$ \theta = \theta_d \quad \text{on} \quad \partial_0 \Omega \qquad (12) $$

where $\vec{n}$ is the unit vector normal to $\partial_q \Omega$ and $\theta_d$ is the prescribed value for the Dirichlet boundary condition. The $\theta$ field can be assimilated as a temperature field, $-\mathbf{K} \mathrm{Grad} \theta$ as a heat flux, and $\mathbf{K}$ as an anisotropic conductivity tensor. If the Dirichlet boundary conditions are compatible with Equation (8), then a solution satisfying

$$ \theta(\mathbf{x}) \neq \text{constant}; \quad \frac{\partial \theta}{\partial \vec{T}} = 0 \qquad (13) $$

is the solution of the boundary value problem presented earlier.

To overcome the singularity of the problem ($\mathbf{K}$ is rank 1), the conductivity tensor is modified as [17]

$$ \mathbf{K}(\mathbf{x}) = \vec{T}(\mathbf{x}) \otimes \vec{T}(\mathbf{x}) + \varepsilon I \qquad (14) $$

where $\varepsilon$ is a small isotropic algorithmic conductivity and $I$ is the second-order identity tensor. Once the problem is solved, the crack path can be located as the iso-value of $\theta$, which passes at the Gauss point, with the largest principal strain. This point is expected to be the closest to the crack mouth, or the crack mouth itself.

### 3.1.2. Application to the location of the idealized crack.
Two main ingredients are defined: (i) the $\vec{T}$ field and (ii) the selection among iso-values on the basis of the state of deformation in the structure. Instead of dealing with the location of maximum damage or the location of maximum nonlocal equivalent strain, we assume here that the line or surface of maximum damage is perpendicular to the direction associated with the maximum principal strain $\varepsilon_{max}$. This consideration restricts the analysis to mode I crack opening where damage grows due to extension perpendicular to the crack path. The $\vec{T}$ field is thus perpendicular to the direction associated to $\varepsilon_{max}$. The boundary value problem in Equations (9), (11), and (12) is solved, and the $\theta$ field is obtained. The second operation consists in selecting the right iso-value. We assume that it is the iso-value that passes at the Gauss point where $\varepsilon_{max}$ is maximal. The algorithm that summarizes the process is presented in Figure 2. In practice, the thermal-like problem is not solved on the full structure, but only on the damaged zone. At the end, the iso-value representing the idealized crack is discretized into linear segments.

## 3.2. Reduction to a one-dimensional problem

The second step consists in evaluating the opening across this idealized crack. We extend the one-dimensional approach presented in Ref. [4] along the linear profiles $\Gamma$ perpendicular to the idealized crack. A set of lines is generated from the middle of each segment of the crack (see Figure 3). Along these lines of direction $\vec{N}$, the normal component of the FE local strain field defined as $\varepsilon_N = \vec{N} \cdot \varepsilon \cdot \vec{N}$ is computed, and then FE nonlocal strain $\bar{\varepsilon}_N$ is also computed (see Figure 4). It serves as an input for the estimate of crack opening. Note that we isolate here the normal component of the strain vector computed perpendicularly to the crack path because we intend to compute a crack
<figure><img src="image_3.png" /><figcaption>Figure 2. Algorithm for the location of the idealized crack.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. One-dimensional profiles (blue) generated from the mesh of the crack.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. One-dimensional profiles generated from the mesh of the crack.</figcaption></figure>

opening. This strain component is regularized prior to comparison to a regularized theoretical strain profile. At this stage, it is also important to remark that the nonlocal strain $\bar{\epsilon}_N$ is not related to the constitutive relation. As we will see next, it is computed to compare regularized strain profiles, and therefore any characteristic length, different from the internal length used in the nonlocal model, may be used. Of course, the characteristic length ought to be such that the weighted average can be computed accurately with the FE discretization at hand. Note also that the weight function itself may be different from that used in the constitutive relation too. In the foregoing calculations, we will use in most cases
the same nonlocal average for the constitutive relation and for the estimate of the crack opening. This is just for the sake of simplicity of the computational implementation.

The estimate of the crack opening follows the idea proposed by Dufour et al. [4]. Upon failure, the displacement field across a crack is a step function with a jump {$U$} at the crack location $x_0$. The strain field is a Dirac delta function, and the nonlocal strain has the same shape as the averaging function $\phi$ taken in the regularization technique defined in Equation (7) with an amplitude of {$U$}. The principle of the calculation of the crack opening estimate is to compare this theoretical nonlocal strain profile with the profile obtained from the FE calculation. The theoretical profile is a function of the crack opening {$U$}. The comparison with the FE profile yields this unknown.

In the original proposal, only the strong formulation was considered, i.e. the crack opening was computed so that both profiles are equal at their maximum $x=x_0$:

$$ \bar{\varepsilon}_{sd}(x_0) = \bar{\varepsilon}_{eq}(x_0) \Rightarrow \frac{\{U\}\phi(0)}{\int_{\Gamma} \phi(x_0 - s)ds} = \bar{\varepsilon}_{eq}(x_0) \Rightarrow \{U\} = \frac{\bar{\varepsilon}_{eq}(x_0)\int_{\Gamma} \phi(x_0 - s)ds}{\phi(0)} \quad (15) $$

Another possibility is to compute the crack opening from the condition that the integral of the two nonlocal strain profiles be the same:

$$ \int_{\Gamma} \bar{\varepsilon}_{eq}(s)ds = \int_{\Gamma} \bar{\varepsilon}_{sd}(s)ds \quad (16) $$

The resulting estimate of the crack opening will be further denoted as *weak estimate* as opposed to the original one denoted as *strong estimate* because maximum nonlocal strains are enforced to be the same in the latter.

Note that the estimate of crack opening is computed from the distribution of the nonlocal equivalent strain compared with a theoretical distribution. The estimate is a function of the current state of strain, and along the idealized crack, closure or opening may occur depending on the applied loads. Estimating the relative sliding between the crack lips would require a proper measure of shear strain that is taken into account neither in the nonlocal equivalent strain implemented in this article ($\bar{\varepsilon}_{eq}(x)$), which is related to local extensions, nor in the derivation of the theoretical distribution corresponding to a displacement discontinuity ($\bar{\varepsilon}_{sd}(x)$) because only mode I opening is considered in the (scalar) expression of the displacement jump {$U$}. Therefore, combined crack opening–crack lips sliding estimate should require an adaptation of the present technique.

In addition to these definitions, the distance between the theoretical and FE profiles gives an indicator of the quality of the solution:

$$ \varepsilon_I(x) = \frac{\int_{\Gamma} |\bar{\varepsilon}_{sd}(s) - \bar{\varepsilon}_{eq}(s)|ds}{\int_{\Gamma} \bar{\varepsilon}_{eq}(s)ds} \quad (17) $$

This indicator quantifies the distance between the actual shape of the nonlocal strain profile compared with a theoretical one describing a displacement discontinuity across the idealized crack. In the case of a damage zone, which does not correspond to a fully developed crack, this indicator provides some additional information: Along the crack path, if $\varepsilon_I(x)$ is above a given value, one may consider that the crack has not formed and that a crack opening displacement cannot be computed with some confidence. If the indicator is below this fixed value, it can be considered that the crack has formed and that a discontinuity has appeared because the FE nonlocal profile is close enough to the theoretical profile. Hence, the value of this indicator along the crack path defines the position of the actual macrocrack tip in a rather general context. Of course, the threshold value below which a discontinuity is considered to occur is somewhat arbitrary. As opposed to existing techniques in which the discontinuity is considered to have appeared if damage is above a fixed value [3,27], the present proposal provides some more consistent information, namely, the comparison between the FE strain distribution around the crack and a theoretical distribution corresponding to a discontinuity.
## 4. EXAMPLES

In this section, we are going first to illustrate the proposed post-processing technique on a three-point bending test. Comparisons with experimental data will be provided. A second example, illustrating the case of curved crack propagation, will also be discussed. In the foregoing examples, we shall use the same internal length in the calculation of the crack opening estimate Equations (15)–(16) and in the condition of damage growth.

### 4.1. Three-point bending test

A single-edge notched beam under three-point bending is considered (see Figure 5). The material parameters used in the nonlocal damage model are given in Table I.

The beam is computed by means of an arc-length procedure using the crack mouth opening displacement (CMOD) as a control parameter. We have chosen an unstructured mesh to avoid potential simplifications due to the symmetry of the problem and to keep a rather general solution processing. The final damage field is given in Figure 6. The load-deflection curve is shown in Figure 7.

The crack is located on the subset defined by the elements whose damage is greater than 0.8: This subset of the mesh and the idealized crack are given in Figure 8. The idealized crack is almost straight, even if the tracking area is nonsymmetric with respect to the center of the beam. One can also notice that the crack bend outward near the top of the beam: This stems from the influence of the loading platen that modifies the strain distribution in this area.

Now, we compute the strong and weak estimates of the crack opening corresponding to the final damage distribution. These estimates are going to be compared with the axial displacement jumps estimated along one-dimensional profiles as shown in Figure 9. In the figures, this jump denoted as *J* will be called *direct calculation*, whereas the weak or strong estimates will be denoted as *estimated* (after post-processing). The jump computed according to the direct calculation is a fit of the axial displacement with a step function. In principle, such a fit is rather easy to perform, probably as easy as the computation of the weak and strong estimates. The direct calculation of the jump is an upper bound of the crack opening displacement. If the actual axial displacement curve departs from a step function, it means that on each side of the crack some deformation occurs, which should not be taken

<figure><img src="image_6.png" /><figcaption>Figure 5. Three-point bending test (dimensions in millimeters).</figcaption></figure>

<table><caption>Table I. Material parameters (three-point bending).</caption><thead><tr><td>Meaning</td><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>Young's modulus</td><td><i>E</i></td><td>38.5 GPa</td></tr><tr><td>Poisson's coefficient</td><td><i>ν</i></td><td>0.24</td></tr><tr><td>Characteristic length</td><td><i>l</i><sub>c</sub></td><td>40 mm</td></tr><tr><td>Damage threshold</td><td><i>Y</i><sub>DO</sub></td><td>310<sup>-5</sup></td></tr><tr><td>Parameter A in compression</td><td><i>A</i><sub>c</sub></td><td>1.25</td></tr><tr><td>Parameter B in compression</td><td><i>B</i><sub>c</sub></td><td>1000</td></tr><tr><td>Parameter A in tension</td><td><i>A</i><sub>t</sub></td><td>1.0</td></tr><tr><td>Parameter B in tension</td><td><i>B</i><sub>t</sub></td><td>9000</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6. Final damage field.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Load-Deflection curve (three-point bending).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Mesh and idealized crack (shaded zone = tracking zone).</figcaption></figure>

into account in the computation of the crack opening. Besides the fact that the present technique is rather systematic and can be applied to very general cases with arbitrary meshes, there is another reason it is more effective to use nonlocal averaging in the proposed estimate of the crack opening: A standard FE mesh provides a poor description of potentially very discontinuous distributions of displacements. The actual value of the jump in the direct calculation may be strongly influenced by the FE sizes. Without any embedded discontinuity, a very slow convergence to the actual value of the displacement jump with the FE size is to be expected. The regularization technique in the nonlocal average provides continuous strain profiles, which may be easily be compared, even in the presence of displacement discontinuities. On top, relatively coarse meshes, with element sizes of one third of the internal length entering in the nonlocal damage model approximately
<figure><img src="image_10.png" /><figcaption>Figure 9. Typical axial displacement along one-dimensional profiles (dimensions in millimeters).</figcaption></figure>

are sufficient to reach a good accuracy in the case of simple structures such as tension bars or beams [1].

The evolution of the opening along the crack is plotted in Figure 10. It can be seen that the weak estimate provides a rather accurate prediction of the displacement jump. In both cases, the errors are moderate.

It is interesting, however, to examine how the estimate behaves upon mesh refinement. A 1994 quadrangular element unstructured mesh with a fine density along the symmetry axis is considered for this purpose (see Figure 11). The estimated opening is plotted in Figure 12. The accuracy of the

<figure><img src="image_11.png" /><figcaption>Figure 10. Estimation of the crack opening (the notch is on the left of the curve): (a) strong approach and (b) weak approach (dimensions in millimeters).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Fine mesh for the estimation of the crack opening.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12. Estimation of the crack opening (fine mesh, the notch is on the left of the curve): (a) strong approach and (b) weak approach (dimensions in millimeters).</figcaption></figure>

strong approach seems to be degraded upon mesh refinement. The error decreases as one approaches toward the tip of the crack. A close scrutiny of the profiles of FE nonlocal strain shows that they do not really exhibit a bell shape, but rather the distribution of nonlocal strain is close to a rectangular one. The weak approach is less sensitive to this difference.

Figure 13 shows the distribution of the error indicator along the idealized crack for the fine mesh. We may observe that the error indicator is rather small in the weak approach, less than 10% except at the bottom of the beam where boundary effect may appear and at the top of the idealized crack where full separation of the surfaces has not occurred yet (fracture process zone). The error indicator is much higher according to the strong approach. This is consistent with the comparisons between the direct calculation and the estimation of the crack opening in Figure 12.

## 4.2. Comparison with experimental results on the three-point bending test

To compare with experimental data, a three-point bending test has been performed on a notched beam. The dimensions of the beam are span = 400 mm, height = 100 mm, and thickness = 50 mm, and the notch is 20 mm high. The test is CMOD controlled.

To measure the crack length and the crack openings, a digital image correlation technique has been implemented (see Figure 14). The crack opening was estimated as the horizontal displacement jump in

<figure><img src="image_14.png" /><figcaption>Figure 13. Error indicator on the crack opening (fine mesh, the notch is on the left of the curve): (a) strong approach and (b) weak approach. $\epsilon_I$ stands for definition (17) of the error, and $\epsilon_D$ stands for the error of the estimated opening with respect to the direct calculation of the opening.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14. (a) Initial picture in undeformed state, (b) picture during the crack propagation, (c) horizontal displacement field, and (d) its three-dimensional view.</figcaption></figure>

the horizontal direction along the vertical crack. The picture frame size was restricted to 55 mm from the notch tip. The picture frame was discretized in 30 horizontal profiles to get the evolution of the crack opening along the height of the beam. For the values of the CMOD of 20 (corresponding to peak load), 30, 40, 50, 60, 80, 100, 150, and 200 µm, a linear curve describing the crack opening as a function of the distance to the top of the picture frame (55 mm from the tip of the notch) has been fitted from the 30 measurements (see Figure 15). The corresponding solid lines have been extrapolated as dashed lines up to the CMOD measurement device located at the bottom of the beam on the lips of the notch.

Parameters in the nonlocal damage model have been fitted from the experimental global response, force versus CMOD curve, with an internal length set equal to 3 times the maximum aggregate size. The fit is shown in Figure 16. The parameters are summarized in Table II.

For each value of the CMOD, the distribution of crack opening obtained in the experiments is compared with the weak and strong estimates. The two numerical estimates of the COD provide similar results (see Figure 17a and c), whereas for large values of the CMOD (200 µm in Figure 17b and d), the strong approach yields a larger error than does the weak approach. The weak approach gives an estimate of the crack opening with an error of 1% for large opening (160 µm) up to 30% for an opening of 5 µm at the crack tip for the peak load. The larger the crack opening, the better the estimation.

The numerical approach always underestimates the experimental crack opening, at least for the present three-point bending test. Although it is not on the safe side for an engineering use, it can be explained from crack propagation considerations: First, the stress state is closer to a plane stress condition at the surface of the beam and closer to a plain strain condition in the core of the beam,
<figure><img src="image_16.png" /><figcaption>Figure 15. Crack opening at different stages of the loading process.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16. Experimental and numerical force versus CMOD responses.</figcaption></figure>

<table><caption>Table II. Parameter fitting using nonlocal damage model.</caption><thead><tr><td>E (GPa)</td><td>A<sub>t</sub></td><td>B<sub>t</sub></td><td>Y<sub>D0</sub></td><td>l<sub>c</sub> (mm)</td></tr></thead><tbody><tr><td>30</td><td>0.9</td><td>4000</td><td>410<sup>-5</sup></td><td>8</td></tr></tbody></table>

which reduces the crack propagation due to confinement; second, due to casting, the material contains less aggregate close to the boundaries and is thus weaker in the sense that aggregates are obstacles for cracking. For these two reasons, the crack propagation is longer on the surface of the specimen than inside the beam. This is clearly illustrated if one looks at the experimental measurements of the crack opening for CMOD of 200 µm (see Figure 17b). The extension of the plot gives a zero opening above the top of the beam, i.e. the neutral axis is out of the beam. It should not be the case except if the applied load is 0. However, in Figure 16, one can see that for CMOD of 200 µm, the bearing capacity is not yet 0. Hence, the crack length ought to be smaller inside the beam compared with the surface of the beam. The two-dimensional plane stress computational model is fitted from the evolution of global quantities: force and CMOD. Thus, this numerical model provides an average crack length between the surface and the core of the beam, whereas the experimental measurements give the maximum values. This is the reason the computed crack length should always be a little smaller than the surface crack length, and the corresponding opening should follow the same trend.
<figure><img src="image_18.png" /><figcaption>Figure 17. Comparison between strong and weak approaches versus experimental crack opening for CMOD = (a) 50 µm and (b) 200 µm. Corresponding errors between numerical and experimental results for CMOD = (c) 50 µm and (d) 200 µm.</figcaption></figure>

### 4.3. Multiple curved crack propagation

In this example, we consider the case of the propagation of two curved cracks in a double-edge notched plate of 1 cm thick under tension (see Fig 18b). Among several specimens presented in Ref. [25] with different distances between the two notches, we have chosen the case of a distance of 10 mm (see Fig 18a)

<figure><img src="image_19.png" /><figcaption>Figure 18. (a) Geometry and (b) boundary conditions of the mechanical test.</figcaption></figure>
for which two curved cracks have been observed experimentally. The loading is carried out by means of an arc-length procedure [20] to follow the global response highlighting a snapback until complete failure for a top displacement of 8.2 µm.

The following set of parameters for the nonlocal damage model presented in Section 2 has been used for the numerical simulation (Table III).

At the final computational step, the crack location procedure is run using the principal strain direction field as a marker of the crack path opened in mode I. Because the procedure selects the isoline passing through the most deformed Gaussian point, only the main (most opened) crack is found in a first run. To find out the second crack path, one has to manually reduce the domain of analysis to the bottom half of the specimen. By running the same procedure a second time, the second crack path is found. Figure 19 shows the two crack paths passing in the central part of the two zones with the largest values of the damage variable *D*.

Figure 20 shows the estimate of the crack opening using the weak approach along the two crack paths. The crack opening in zones where only one crack is found and the sum of the two openings where the two cracks interact are both between 7 and 8 µm. Although these values are very close to the top displacement (8.2 µm), they stand below because the local crack orientation is not perpendicular to the loading direction.

Figure 21 In this example, which is primarily aimed at illustrating crack localization in the very complex case of curved crack propagation and interacting cracks, a rather coarse mesh has been used regarding the description of the notch tip, and a good quality of the estimate should not be expected. However, the model error stands below 15% except close to crack tips where some peak values around 30% are observed.

## 5. CONCLUSION

An extension of the technique proposed to estimate crack openings from the FE calculation by Dufour *et al.* [4] to two-dimensional computations of mode I fracture has been proposed. The first step is the

<table><caption>Table III. Model parameters for the numerical simulation using nonlocal damage model.</caption><thead><tr><td><i>E</i> (GPa)</td><td><i>A</i><sub>t</sub></td><td><i>B</i><sub>t</sub></td><td><i>Y</i><sub>DO</sub></td><td><i>l</i><sub>c</sub> (mm)</td></tr></thead><tbody><tr><td>31</td><td>0.95</td><td>9 000</td><td>9.710<sup>-5</sup></td><td>4</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>Figure 19. Comparison between the two crack paths and the damage field for a displacement of 8.2 µm.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 20. Crack opening using the weak approach along the two crack paths for a displacement of 8.2 µm.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21. Model error along the two crack paths for a displacement of 8.2 µm.</figcaption></figure>

calculation of the location of the idealized crack (line perpendicular to the direction of maximum principal strain in the fracture process zone). For this purpose, we have applied the global tracking method to compute the idealized crack location as a post-treatment procedure applicable to any mode I failure model. Then, the original one-dimensional problem devised by Dufour et al. [4] is recovered as profiles of normal deformation orthogonal to the idealized crack direction are considered. For each profile, an estimate of the crack opening and an error indicator are computed by comparing FE deformation profiles and theoretical profiles corresponding to a displacement discontinuity. Two estimates have been considered: In the strong approach, the maxima of the profiles are assumed to be equal; in the weak approach, the integrals of each profile are equal.

The two estimates are restricted to mode I crack opening, and due to the constitutive model implemented, only mode I crack propagation is considered too. In addition to the crack opening estimate, the technique yields also an error indicator, which quantifies whether a crack has fully developed or not. Along the idealized crack, if it is above a given value, one may consider that the crack has not formed and that a crack opening displacement cannot be computed with some confidence. If the indicator is below this fixed value, it can be considered that the crack has formed and that a discontinuity has appeared because the FE nonlocal profile is close enough to the theoretical profile. Hence, the value of this indicator along the crack path defines the position of the actual macrocrack tip. Of course, the threshold value below which a discontinuity is considered to occur needs to be decided *a priori*. It depends on the desired accuracy in the definition of a real crack, which depends on the constitutive relations. If they exhibit a zero stress for a finite strain, the threshold could be determined accordingly.
Two-dimensional numerical calculations show that the weak estimate performs better than does the strong one. Error indicators, defined as the distance between the numerical and theoretical profiles, of magnitude less than a few percentages may be achieved. Results have been compared with experimental measurements performed using the digital image correlation technique. Errors stand from less than 1% for large values of the opening (160 µm) to 30% for very small opening (5 µm). This should be expected because most probably a full crack has not formed yet when the opening is so small. The example of two interacting cracks in a tensile concrete plate shows that the crack detection technique is capable of handling curved and interacting cracks.

A direct extension of this work could be used to perform fluid flow calculations from the computed
crack openings (e.g., using Poiseuille's model). The values of crack openings may also be folded into
continuum-based hydromechanical models in which the permeability is a function of damage and also
of crack openings [23].

ACKNOWLEDGEMENTS

The financial support by the Agence Nationale de la Recherche under project Contifiss is gratefully acknowledged. The third author also acknowledges partial support from ERC Advanced Grant 27769 Failflow. Cédric Giry (PhD student at 3SR) is thanked for his help in running successfully some of the numerical simulations.

REFERENCES

1. Bažant ZP, Pijaudier-Cabot G. Nonlocal continuum damage, localization instability and convergence. *Journal of Applied Mechanics ASME* 1988; **55**:287–294.

2. Choinska M, Khelidj A, Chatzigeorgiou G, Pijaudier-Cabot G. Effects and interaction of temperature and stress level related damage on permeability of concrete. *Cement and Concrete Research* 2007; **37**:79–88.

3. Comi C, Mariani S, Perego U. An extended fe strategy for transition from continuum damage to mode i cohesive crack propagation. *International Journal for Numerical and Analytical Methods in Geomechanics* 2007; **310**(2):213–238.

4. Dufour F, Pijaudier-Cabot G, Choinska M, Huerta A. Extraction of a crack opening from a continuous approach using regularized damage models. *Computers and Concrete* 2008; **50**(4):375–388.

5. Dumstorff P, Meschke G. Crack propagation criteria in the framework of x-fem structural analyses. *International Journal for Numerical and Analytical Methods in Geomechanics* 2007; **310**(2):239–259.

6. Feist C, Hofstetter G. Three-dimensional fracture simulations based on the sda. *International Journal for Numerical and Analytical Methods in Geomechanics* 2007; **310**(2):189–212.

7. Hearn N, Lok G. Measurement of permeability under uniaxial compression-a test method. *ACI Materials Journal* 1998; **95**:691–695.

8. Hillerborg A, Modeer M, Pertersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. *Cement and Concrete Research* 1976; **6**:773–782.

9. Jirasek M, Rolshoven S, Grassl P. Size effect on fracture energy induced by non-locality. *International Journal for Numerical and Analytical Methods in Geomechanics* 2004; **280**(7-8):653–670.

10. Larsson R, Steinman P, Runesson K. Finite element embedded localization band for finite strain plasticity based on a regularized strong discontinuity. *Mechanics of Cohesive-Frictional Materials* 1998; **4**:171–194.

11. Mariani S, Perego U. Extended finite element method for quasi-brittle fracture. *International Journal for Numerical Methods in Engineering* 2003; **580**(1):103–128.

12. Mazars J, Pijaudier-Cabot G. Continuum damage theory - application to concrete. *Journal of Engineering Mechanics* 1989; **1150**(2):345–365.

13. Mazars J, Pijaudier-Cabot G. From damage to fracture mechanics and conversely: a combined approach. *International Journal of Solids and Structures* 1996; **33**:3327–3342.

14. Mediavilla J, Peerlings RHJ, Geers MGD. A robust and consistent remeshing-transfer operator for ductile fracture simulations. *Computers and Structures* 2006; **840**(8-9):604–623.

15. Moes N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**:131–150.

16. Oliver J, Huespe AE. Theoretical and computational issues in modelling material failure in strong discontinuity scenarios. *Computer Methods in Applied Mechanics and Engineering* 2004a; **1930**(27-29):2987–3014.

17. Oliver J, Huespe AE. Continuum approach to material failure in strong discontinuity settings. *Computer Methods in Applied Mechanics and Engineering* 2004b; **1930**(30-32):3195–3220.

18. Oliver J, Huespe AE, Pulido MDG, Chaves EWV. From continuum mechanics to fracture mechanics: the strong discontinuity approach. *Engineering Fracture Mechanics* 2002; **69**:113–136.

19. Peerlings RHJ, de Borst R, Brekelmans WAM, de Vree JHP. Gradient enhanced damage for quasi-brittle materials. *International Journal for Numerical Methods in Engineering* 1996; **39**:937–953.
20. Pegon P, Anthoine A. Numerical strategies for solving continuum damage problems with softening: Application to the homogenization of masonry. *Computers and Structures* 1997; **640**(1-4):623-642.

21. Pijaudier-Cabot G, Bažant Z. Nonlocal damage theory. *Journal of Engineering Mechanics* 1987; **113**:1512-1533.

22. Pijaudier-Cabot G, Dufour F. Nonlocal damage model: boundary and evolving boundary effects. *European Journal of Environmental and Civil Engrg.*, in press, 2010.

23. Pijaudier-Cabot G, Dufour F, Choinska M. Permeability due to the increase of damage in concrete: from diffuse to localised damage distributions. *Journal of Engineering Mechanics* 2009; **1350**(9):1022-1028. http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000016.

24. Planas J, Elices M, Guinea GV. Cohesive cracks versus nonlocal models: Closing the gap. *International Journal of Fracture* 1993; **63**:173-187.

25. Shi C, van Dam AG, van Mier JGM, Sluys LJ. Crack interaction in concrete. *Conference Proceedings of EUROMAT99*, 1999; **6**:27-30.

26. Simo JC, Oliver J, Armero F. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. *Computational Mechanics* 1993; **12**:277-296.

27. Simone A, Wells GN, Sluys LJ. From continuous to discontinuous failure in a gradient-enhanced continuum damage model. *Computer Methods in Applied Mechanics and Engineering* 2003; **1920**(41-42):4581-4607.

28. Simone A, Askes H, Sluys LJ. Incorrect initiation and propagation of failure in non-local and gradient-enhanced media. *International Journal of Solids and Structures* 2004; **41**:351-363.

29. Sugiyama T, Bremmer TW, Holm TA. Effect of stress on gas permeability in concrete. *ACI Materials Journal* 1996; **93**:443-450.

30. Wells G, Sluys LJ, de Borst R. Simulating the propagation of displacement discontinuities in a regularized strain-softening medium. *International Journal for Numerical Methods in Engineering* 2002; **530**(5):1235-1256.