# How to extract a crack opening from a continuous damage finite-element computation?

Frédéric Dufour, Gilles Pijaudier-Cabot, Marta Choinska, Antonio Huerta

► To cite this version:

Frédéric Dufour, Gilles Pijaudier-Cabot, Marta Choinska, Antonio Huerta. How to extract a crack opening from a continuous damage finite-element computation?. International Conference on Computational Fracture and Failure of materials and Structures, Jun 2007, Nantes, France. pp.1-8. hal-04703248

HAL Id: hal-04703248

https://hal.science/hal-04703248v1

Submitted on 19 Sep 2024

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
# How to extract the crack opening from a continuous damage finite-element computation ?

F. Dufour, G. Pijaudier-Cabot & M. Choinska

R&DO - Institut GeM, Ecole Centrale de Nantes, Nantes, France

A. Huerta

Lab. de Càlcul Numèric, Univ. Politècnica de Catalunya, Barcelona, Spain

Doubtlessly, for assessing the integrity and the strength of concrete structures to environmental attacks it would be ideal to model the transition between diffuse damage and a discontinuous description. We propose, in a first approach, to extract an equivalent crack opening using a non local Mazars’ damage model. We consider a discontinuous displacement field with a given jump $[U]$ from which we compute the equivalent non local strain $\tilde{\varepsilon}([U])$ governing the damage variable in Mazars’ model. Besides, from a classical FE computation, we obtain a distribution of the state variable $\tilde{\varepsilon}^{FE}$ with the same shape due to the non local weight function. Finally the displacement jump is computed by setting $\tilde{\varepsilon}([U]) = \tilde{\varepsilon}^{FE}$ at the discontinuity. The proposed approach is applied on a tension rod with an error measure between both solutions to evaluate how the non local solution approaches asymptotically the representation of a discontinuity. Results show good agreement for large damage.

# 1 INTRODUCTION

The assessment of the diffuse damage and the crack opening is a key parameter to estimate the integrity and the strength of concrete structures to environmental attacks ruled by diffusion and/or Darcy’s transfers. Transfer properties are even the main serviceability criteria of structures such as confinement vessels of nuclear power plants.

Continuum enhanced and integral damage models are capable to represent diffuse damage, crack initiation and propagation (Pijaudier-Cabot and Bažant 1987). They regard damage process as an ultimate consequence of a gradual loss of material integrity. These models, however, do not predict the crack opening.

Another approach to model failure consists in representing the discontinuity within the material explicitly (e.g. cohesive crack model (Hillerborg et al. 1976)). Those models relate the crack opening to the stress level and are based on the linear elastic (or plastic) fracture mechanics. They are usually used when the Fracture Process Zone size is negligible compared to the structure size. Furthermore one needs to know the crack location beforehand and set the discontinuous model at this place.

Trying to overcome disadvantages of both approaches (continuum damage and discrete crack)

many authors have tried to bridge both theories. The key parameter of such an approach is the threshold upon which the discontinuity is activated and the crack starts to open. Usually, it is considered that the discontinuity appears when damage, stresses or energy reach a certain threshold fixed beforehand. However, this transition remains arbitrary (Comi et al. (2007), Simone et al. (2003)) and the equivalence between state variables of both approaches is not properly defined.

Instead of trying to bridge these two classes of models as done by Mazars and Pijaudier-Cabot (1996), Planas et al. (1993) and many others, we use the nonlocal damage model to compute the initiation and propagation of the diffuse damage and extract the discrete key parameter, i.e. the crack opening, using the concept of strong discontinuity introduced by Simo et al. (1993) and widely used over the last decade (e.g. Oliver et al. (2002), Larsson et al. (1998)).

We first recall briefly the strong discontinuity approach, then the non local damage model and finally we apply our procedure to a tension bar where the crack location is fixed at the center.
## 2 NON LOCAL STRONG DISCONTINUITY

Let imagine the failure of an unidimensional elastic brittle bar in elongation, one thinks that at a late stage of damaging the bar will be splitted into two parts with a crack somewhere depending on material heterogeneities. We assumed that at failure the material is fully unloaded (no stress) according to damage models (see Equation 5) it also means unstrained except at the crack. The displacement field will be constant piecewise with a step at the crack position and its continuous part vanishes. This concept can be seen as an asymptotic development of the strong discontinuity approach since our displacement field is purely discontinuous whereas in the latter approach it is the sum of a continuous and a discontinouous field.

We choose this pure kinematic mode of loading (elongation rather than traction) since damage model are driven by positive strain (see Equations 7 and 8). Moreover computations need to be driven by strain to model the post-peak regime in finite element code. Let's take a bar of length L oriented by $\vec{x}$ clamped at $x = 0$ and with a constant velocity $v$ at the other end ($x = L$). Let's assume that the crack will occur at a specific abscissa $x = x_0$ known before hand (see Figure 1.(a)).

If we assume that the bar is elastic perfectly brittle, after failure the displacement profile (see Figure 1.(b)) will be:

$$u_{sd}(x, x_0) = [U] H_{\Gamma}(x - x_0) \quad (1)$$

where underscore sd stands for strong discontinuity, $H_{\Gamma}$ is the heaviside function and $[U]$ is the displacement jump across the crack. Once the crack location is known, the displacement jump is the unique unknown of our problem. From Equation 1, we derive the strain field (see Figure 1.(c)) with the classical symmetric gradient operator $\nabla^s$:

$$\varepsilon_{sd}(x, x_0) = \nabla^s u_{sd} = [U] \delta(x - x_0) \quad (2)$$

where $\delta(x)$ is the Dirac function. To this strain field is associated a non local strain (see Figure 1.(d)) according to the same procedure (see Equation 9) used in the non local damage model.

$$\bar{Y}_{sd}([U], x, x_0) = \frac{\int_{\Omega} \phi(x - s) \varepsilon_{sd}(s, x_0) ds}{\int_{\Omega} \phi(x - s) ds} \quad (3)$$

where $\phi(x)$ is the weight function. Using Equations 2 and 3, and the properties of the Dirac function we get:

$$\bar{Y}_{sd}([U], x, x_0) = \frac{[U] \phi(x - x_0)}{\int_{\Omega} \phi(x - s) ds} \quad (4)$$

which is a non local measure of the local displacement jump.

<figure><img src="image_2.png" /><figcaption>Figure 1: Conceptual model of an unidimensional bar in tension. (a) Geometry of the bar with boundary conditions, (b) Displacement profile, (c) Strain profile and (d) Regularized strain profile.</figcaption></figure>

## 3 NON LOCAL DAMAGE APPROACH

### 3.1 Damage model

The Mazars' damage model (Mazars 1984) is used in the finite element computations for simulating the progressive failure of the bar. In this model the stress is expressed as follows:

$$\sigma = (1 - D)\varepsilon \quad (5)$$

where D is the damage scalar variable generally written as:

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

where $\alpha_t$ and $\alpha_c$ depend on the strain tensor. In the case of an uniaxial loading in elongation ($\alpha_t = 1$ and $\alpha_c = 0$) Equation 6 is reduced to:

$$D = 1 - \frac{Y_{D0}(1 - A_t)}{Y} - \frac{A_t}{e^{[B_t(Y - Y_{D0})]}} \quad (7)$$

where $A_t$, $B_t$ and $Y_{D0}$ are the parameters of the model and Y is defined by:

$$Y = \max(Y, \varepsilon_{eq}) \quad (8)$$

with initially $Y = Y_{D0}$ and the equivalent strain defined as:

$$\varepsilon_{eq} = \sqrt{\sum_{i=1}^{3} \langle \varepsilon_i \rangle_{+}^2} \quad (9)$$
where $\langle \epsilon_i \rangle_+$ denotes the positive part of the principal strains. Because of softening it is well established that this model is submitted to mesh dependency due to damage localization. A regularization technique needs to be employed to keep the objectivity of the numerical problem.

## 3.2 Integral damage model

In order to overcome the loss of ellipticity of the solution, the integral version of the Mazars' damage model (Pijaudier-Cabot and Bažant 1987) can be used. In this approach the equivalent strain $\bar{\epsilon}_{eq}$ is a function of the strain surrounding the computational point.

$$ \bar{\epsilon}_{eq}(x) = \frac{\int_{\Omega} \phi(x-s) \epsilon_{eq}(s) ds}{\int_{\Omega} \phi(x-s) ds} \quad (10) $$

This is mechanically justified by the interactions (stress redistribution when micro-cracks occur) due to heterogeneities of the material.

Since $Y$ is the state variable associated to the internal variable $D$ it cannot decrease in time (Equation 8), however the strain and consequently the regularized strain decrease in the neighborhood of the crack (close but not at the crack) due to elastic unloading. Therefore $Y$ and $\bar{Y}_{sd}$ cannot be directly compared as the second one decreases when elastic unloading occurs. Therefore, we construct a new variable named $\bar{Z}_{eq}$ defined by:

$$ \bar{Z}_{eq} = \begin{cases} \bar{\epsilon}_{eq} & \text{if } \bar{\epsilon}_{eq} \le Y_{D0} \\ \bar{Y} & \text{otherwise} \end{cases} \quad (11) $$

## 3.3 Computation of the crack opening

The idea is to compare both regularized strain fields $\bar{Z}_{eq}$ and $\bar{Y}_{sd}$ as the first one represents the state variable of the damage model and the latter one is easily accessible by the strong discontinuity description of the failure. The only unknown $[U]$ is computed by setting that both regularized strain are equal at the place of the crack supposed to be known beforehand:

$$ [U] = \frac{\bar{Z}_{eq}(x_0) \int_{\Omega} \phi(x_0 - s) ds}{\phi(0)} \quad (12) $$

Both entire field are then determined and can be compared by defining an absolute error field $\Delta$ such as:

$$ \Delta(x, x_0, [U]) = \bar{Y}_{sd}(x, x_0, [U]) - \bar{Z}_{eq}(x) \quad (13) $$

and normalized it by the numerical field, one gets a relative error field:

$$ \Delta^r(x_0, [U]) = \frac{\int_{\Omega_d} \|\Delta(s, x_0, [U])\| ds}{\int_{\Omega_d} \bar{Z}_{eq}(s) ds} \quad (14) $$

As explained further in the numerical application, the integration domain is reduced to elements where damage has occured ($\Omega_d$).

Our analysis does not imply any modification of a finite element code as we extract the regularized strain field and compare it to an analytical one in a post-treatment phase. In order to use our approach for a wider range of application, the crack location will have to be found from the damage and/or strain fields out of the finite element computations.

# 4 NUMERICAL COMPARISON

We first present the numerical results with the Gaussian weight function and, at the end, we look at the influence of the shape of the weigth function.

## 4.1 Global response

We numerically solve the problem defined in Figure 1.(a) for a unit length bar and we force the crack position by setting a weak element at the center with a smaller Young modulus $E_{we}$. In a first step, the most often used Gaussian weight function is chosen:

$$ \phi_g(x-s) = \exp\left(\left(\frac{2\|x-s\|}{l_c}\right)^2\right) \quad (15) $$

where $l_c$ is the internal length of the model. The material parameters are presented at Table 1.

$$ \begin{align*} E &= 37,7 \text{ GPa} & Y_{D0} &= 10^{-4} \\ E_{we} &= 37,6 \text{ GPa} & A_t &= 1 \\ \nu &= 0,2 & B_t &= 14000 \\ & & l_c &= 0.32 \text{ m} \end{align*} $$

<caption>Table 1: Material parameters for finite element computations.</caption>

$B_t$ and $l_c$ are chosen such that the process zone does not reach the bar ends and is large enough to avoid snap-back since our computations are driven by the end displacement. Figure 2 shows the material behaviour with such a set of parameters.

<figure><img src="image_3.png" /><figcaption>Figure 2: Material response for the chosen set of parameters.</figcaption></figure>
An odd number of finite elements must be taken to always set the weak element (crack position) at the same place independently of the mesh density. We found out that for this given problem 61 elements is enough to reach the finite element convergence and thus the mesh independency domain on the numerical response. The global response is shown in Figure 3. In order to highlight the evolution, 6 loading steps are indicated by letters from A to F. We will refer to them throughout the analysis.

<figure><img src="image_4.png" /><figcaption>Figure 3: Global response of the bar in traction.</figcaption></figure>

## 4.2 Damage and strain profiles

Figure 4 shows the evolution of damage profile during loading process. As soon as damage starts to increase the width of the damaged zone is roughly twice the internal length and remains constant throughout the fracture process. Subsequently the damage reaches at the center of the bar a value very close to 1, smaller though, according to Equation 7.

<figure><img src="image_5.png" /><figcaption>Figure 4: Damage profile for particular loading steps.</figcaption></figure>

Figure 5 shows profiles along the bar for the particular loading steps from A to F of the non local equivalent strain $Z_{eq}$ issued from finite element computations and the regularized strain $\bar{Y}_{sd}$ of our conceptual approach.

Both profiles coincide at the center of the bar by definition of $[U]$ and have the shape of the weight function. One can have chosen another definition of $[U]$ by minimizing the difference between both profiles. Although our error measurement is not optimized as a

<figure><img src="image_6.png" /><figcaption>Figure 5: Numerical and analytical regularized strain profiles for particular loading steps.</figcaption></figure>

smaller error could be obtained, we kept our definition of $[U]$ as it sounds more physical. The numerical profile is the widest since the engineer strains for non local finite element approach do not localize in a single element but in a small band included in the process zone (see Figure 14). This is one aspect of the relative incapacity of non local damage models to represent exactly the formation of a macro-crack as non local interactions are still modelled even for large values of damage. Besides, if the engineer strains do localize in a single finite element the numerical solution would be mesh dependent.

## 4.3 Error measurement

From Figure 5, we decided to reduce the integration domain of our error to the damaged zone $\Omega_d$ in order to get a measurement independent of the structure size. Indeed, at early stages (A to E), the difference between $Z_{eq}$ and $\bar{Y}_{sd}$ profiles is mainly outside the fracture process zone since our asymptotic strong discontinuity concept is defined for fully unloaded material. Therefore, if integrated over the structural domain, the error would be directly related to the relative size of the FPZ and the structure size and would be much larger for a longer bar than a smaller one.

<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of the error vs loading.</figcaption></figure>

As expected using the asymptotic concept of the strong discontinuity approach the error (see Figure 6)
is quite large just after the peak load since the macro-crack is not yet clearly formed, in other words the material is not yet fully unloaded and unstrained outside the macro-crack. Afterwards the error decreases rapidly to reach a limit value of 4% when the macro-crack is widely opened. This value represents the level of ability of the integral non local damage model using the Gaussian weight function to simulate a failure in terms of kinematic variables.

Within the process zone the error is the largest where the slope of the weight function is the steepest (see Figure 7). As the error is due to a wider profile of the non local strain, the error, i.e. the vertical distance between profiles, is the largest where the non local strain profile is the steepest.

<figure><img src="image_8.png" /><figcaption>Figure 7: Error profiles for particular loading steps.</figcaption></figure>

Once the crack is fully opened the increase of its opening is equal to the displacement applied at the active end of the bar as no further strain are stored in the material. Therefore the plot of the displacement jump versus the imposed displacement tends to the bisecting line (see Figure 8).

<figure><img src="image_9.png" /><figcaption>Figure 8: Evolution of the displacement jump vs loading.</figcaption></figure>

## 4.4 Influence of the weight function

One of the numerical parameters of an integral damage model is the shape of the weight function. Whatever the shape, the function must have an integral

over the domain of interest equal to 1 in order to retrieve the local strain if it is an homogeneous field. We choose, in addition of the Gaussian function already used, two other weight functions, i.e. the bell function:

$$ \phi_b(x-s) = \begin{cases} \left(1 - \left(\frac{\|x-s\|}{0.9l_c}\right)^2\right)^2 & \text{if } \|x-s\| < 0.9l_c \\ 0 & \text{otherwise} \end{cases} \qquad (16) $$

and an exponential function:

$$ \phi_e(x-s) = \exp\left(-\left(\frac{2\|x-s\|}{0.6l_c}\right)\right) \qquad (17) $$

<figure><img src="image_10.png" /><figcaption>Figure 9: Damage profiles given by the three different weight functions.</figcaption></figure>

The bell and the Gaussian functions are quite similar with a zero continuous tangent at the center (see Figure 10). The exponential function is directly inspired by the function (18) proposed by Peerlings 1999 in order to get the analytical equivalence between integral and gradient approaches of non local damage models.

$$ \phi_p(x-s) = \frac{1}{4\pi l_c^2 \|x-s\|} \times \exp\left(-\left(\frac{\|x-s\|}{l_c}\right)\right) \qquad (18) $$

This original Peerlings' function has an infinite weight at the center ($x = s$) which is a major drawback for a numerical use.

Constants (i.e. 0.6 and 0.9) multiplying the internal length $l_c$ for the bell and the exponential weight functions are found such as to obtain the same width for the damage band, with a different profile though (see Figure 9). We choose this common features as the damage band width is directly related to the material properties and is the one usually measured experimentally by acoustic emission technique for instance (Granger et al. 2007).
<figure><img src="image_11.png" /><figcaption>Figure 10: Profile of the three weight functions compared.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Global force displacement response for the three weight functions.</figcaption></figure>

At the end, the three normalized weight functions are plotted at Figure 10 for comparison. As expected the steepest function (exponential) is also the narrowest with the largest weight at the center but all of them are still quite similars and yield roughly the same damage profile (see Figure 9).

Figure 12 and 13 show the comparison of the analytical profile with the numerical profile of the non local strain for the bell function and the exponential function, respectively. From those plots the exponential weight function seems more adapted to represent a strong discontinuity.

<figure><img src="image_13.png" /><figcaption>Figure 12: Profiles of both analytical and numerical regularized strains for the bell weight function.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Profiles of both analytical and numerical regularized strains for the exponential weight function.</figcaption></figure>

Indeed engineer strain profiles for the three weight functions clearly show in Figure 14 that the one given by the exponential is much narrower than others due to the specific shape of the weight function.

<figure><img src="image_15.png" /><figcaption>Figure 14: Comparison of engineer strain profiles for the three weight functions.</figcaption></figure>

At the end, as expected from previous profiles the error given by the bell function is similar to the Gaussian function about 4% whereas the exponential function gives a smaller error of about 2%.

<figure><img src="image_16.png" /><figcaption>Figure 15: Comparison of the error evolution during loading process for the three weight functions.</figcaption></figure>

# 5 CONCLUSIONS

We have derived an analytical model based on the strong discontinuity approach for the failure of a
bar loaded in traction and compare it from a kinematic viewpoint to finite element computational results given by a continuous non local damage model. Once the macro-crack is formed, our approach is capable to extract a crack opening displacement from a finite element computation using non local damage model with a good accuracy. This computation is done as a post-treatment and is, therefore, costless and easy to add in an existing finite element code.

Unfortunately, this method still requires some user's expertise since the error evolution vs. number of element is not monotonously decreasing. The computation of the displacement jump can be nearly null if the mesh is not fine enough to properly describe the strain profile. At the limit, one can use a mesh fine enough to reach mesh independency of the global response and still get the engineer strain localized in a single finite element and consequently get a perfect matching with the strong discontinuity approach. The error asymptotically increases with the number of elements to reach the “model error” of few percents. Therefore, this procedure requires a very fine mesh in the process zone to get the asymptotic error between non local models and the strong discontinuity description of the failure.

This error slightly depends on the shape of the weight function. The steeper at the peak is the weight function the smaller is the error. At the end, our approach is capable to compute the crack opening from a non local damage finite element calculation with an error around 2% for the exponential weight function.

As the numerical engineer strain profile across the crack depends on the shape of the weight function, some input should be given by new field measurement techniques (e.g. correlation) in order to select which one of the weight functions gives the closer strain response.

On one hand this approach needs to be extrapolated to real reinforced structures with multi crack pattern. The key issue will be to locate in space several cracks from the damage and/or the strain field(s) given by the finite element computation. The crack will have to be discretised somehow in order to apply our 1D approach along the crack to get the variation of the opening along the crack. We expect our approach to be valid as long as cracks do not bifurcate even for multi crack pattern as long as we are able to locate and treat individually each crack.

On the other hand, it needs to be validated using experimental tests. It will be performed at the specimen scale with no reinforcement in order to get a single crack whose opening can be measured precisely. We will use specimen made of mortar for having a fracture process zone small enough compared to the specimen size in order to highlight structural effect. We choose to perform those tests on concrete discs

loaded in compression in the diameter plane (splitting test). Measurement of the displacement field by age technique carried out during the tests will permit to extract the crack opening which will be, at the end, compared to the one obtained by our approach from a continuous non local damage modelling.
REFERENCES

Comi, C., S. Mariani, and U. Perego (2007). An extended fe strategy for transition from continuum damage to mode i cohesive crack propagation. *International Journal for Numerical and Analytical Methods in Geomechanics*. In press.

Granger, S., A. Loukili, G. Pijaudier-Cabot, and G. Chanvillard (2007). Experimental characterization of the self healing of cracks in an ultra high performance cementitious material: mechanical tests and acoustic emission analysis. *Cement and Concrete Research*. In Press.

Hillerborg, A., M. Modeer, and P. E. Pertersson (1976). Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. *Cement and Concrete Research* 6, 773-782.

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

Mazars, J. (1984). *Application de la mécanique de l'endommagement au comportement non linéaire et à la rupture de béton de structure*. Ph. D. thesis, Université Pierre et Marie Curie, In French.

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

Oliver, J., A. E. Huespe, M. D. G. Pulido, and E. W. V. Chaves (2002). From continuum mechanics to fracture mechanics: the strong discontinuity approach. *Engineering Fracture Mechanics* 69, 113-136.

Peerlings, R. H. J. (1999). *Enhanced damage mod-
elling for fracture and fatigue*. Ph. D. thesis,
Technische Universiteit Eindhoven.

Pijaudier-Cabot, G. and Z. Bažant (1987). Nonlocal damage theory. *Journal of Engineering Mechanics* 113, 1512–1533.

Planas, J., M. Elices, and G. V. Guinea (1993). Cohesive cracks versus nonlocal models: Closing the gap. *International Journal of Fracture* 63, 173–187.

Simo, J. C., J. Oliver, and F. Armero (1993). An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. *Computational Mechanics* 12, 277–296.

Simone, A., G. N. Wells, and L. J. Sluys (2003).
From continuous to discontinuous failure in a

gradient-enhanced continuum damage model.
*Computer Methods in Applied Mechanics and Engineering* 192(41-42), 4581–4607.