# Modelling Crack Propagation in Concrete Structures with a Two Scale Approach

Khalil Haidar, Jean-François Dubé, Gilles Pijaudier-Cabot

► To cite this version:

Khalil Haidar, Jean-François Dubé, Gilles Pijaudier-Cabot. Modelling Crack Propagation in Concrete Structures with a Two Scale Approach. International Journal for Numerical and Analytical Methods in Geomechanics, 2003, 27 (13), pp.1187-1205. 10.1002/nag.318 . hal-01007080

HAL Id: hal-01007080

https://hal.science/hal-01007080

Submitted on 28 Jan 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
# Modelling crack propagation in concrete structures with a two scale approach

Khalil Haidar¹, Jean François Dube² and Gilles Pijaudier-Cabot¹

¹R&DO, Laboratoire de Génie Civil de Nantes Saint Nazaire, Ecole Centrale de Nantes, 1 rue de la Noe, BP 92101, F-44321 Nantes Cedex 3, France

²Laboratoire de Mécanique et Génie Civil, Université Montpellier II/CNRS, Place Eugène Bataillon, cc048, F-34095 Montpellier Cedex 5, France

A simplified computational technique based on a refined global local method is applied to the failure analysis of concrete structures. The technique distinguishes the scale of the structure, modelled with large size finite elements, from the scale at which material non-linearity occurs due to progressive cracking and macro-crack propagation. The finite element solution is split into two parts: a linear elastic analysis on a coarse mesh over the entire structure and a non-linear analysis over a small part of the structure where a dense finite element grid is employed. In the non-linear calculation, a non-local damage model is implemented. These two computations are coupled with the help of an iterative scheme. The size and location of the region where a non-linear analysis is performed, is adapted to follow the development of the damage zone. Numerical examples of mode I fracture of concrete specimens with straight and curved cracks are presented.

**KEY WORDS:** concrete; damage; multiscale approaches

## 1. INTRODUCTION

The prediction of cracking and the evaluation of the long term safety of large concrete structures require a sound and robust computational method of assessment. On top of the difficulties involved in the description of the inception and propagation of cracks at the material level and the implementation of adequate numerical algorithms for solving non-linear problems, another issue is to be able to deal with finite element models which contain a large number of degrees of freedom. In the case of large size structures, such as cooling towers or containment vessels, there is a balance to achieve between two extremes: On one side practical numerical models are made of finite elements with a relatively large size (e.g. elements of the order of 1 m in the industrial example presented recently by Badel and et al. [1]). On the other side, non-linear constitutive relations which describe progressive cracking, e.g. based on plasticity or damage mechanics,
require finite elements which are rather small in order to achieve a good accuracy. One needs to describe material non-linearity which remains concentrated in very small regions of the structure, where cracking occurs. As an example, the accuracy of the discretization of damage zones (or fracture process zones) in the case of non-local damage models is directly related to the internal length entering in the constitutive relations. With the integral non-local damage model and with a linear interpolation of the displacement field, experience shows that a reasonable accuracy can be achieved when the finite element size is lower than one third of the internal length [2, 3], which means that the width of the fracture process zone is of the order of six finite elements roughly. The magnitude of the internal length is a few centimetres, say 50 mm typically. Hence, the maximum size of the finite elements in the fracture process zone should be about 17 mm if linear elements are used. One can easily figure that for very large structures, this requirement calls for mesh optimization between large regions of the structure, which stay in the elastic regime with small strain gradients, and small zones where stress concentration and/or progressive failure occur.

There are several ways to optimize the interpolation of the displacement field. *h* or *p*-mesh adaptations are a possibility and there are many of such proposals in the literature (see e.g. References [4–6]). It requires a robust error indicator, a mesh generator capable of creating an optimized grid, and efficient projection techniques of the discrete variables from one mesh to another, if meshes are to be changed in the course of the computation.

For a very large structure (such as a dam for instance), one may also consider that the fracture process occurs in such a small region that it can be viewed as a displacement or strain discontinuity. Many proposals followed this idea, including the strong discontinuity approach initiated by Simo *et al.* [7] and finite elements with embedded discontinuities [8–10]. These techniques require crack (or discontinuity) inception and propagation criteria. In some cases, adaptive meshing may be required too because accurate descriptions of crack initiation and crack propagation can be quite demanding in terms of finite element fineness. One drawback of such discrete models is the lack of direct relationship between local constitutive relations and the mechanical response of embedded discontinuities. With local constitutive relations, inception of damage, its localization into a narrow zone and the propagation of a macroscopic crack are the result of the structural analysis. A relationship between continuum damage modelling and fracture energy can be derived [11] but it remains difficult to derive the direction of propagation of a crack from such a continuum analysis. As a consequence embedded discontinuity approaches, in their various forms, leave aside some benefits of continuum descriptions.

The extended finite element approach is another possibility developed initially in the context of linear elastic fracture mechanics. Application for cohesive crack growth was performed recently by Moes *et al.* [12]. Further developments of this method are necessary before it can correctly predict crack growth paths in cohesive material (e.g. in concrete). The difficulty in X-FEM is that one has to decide *a priori* on the proper displacement enhancement in the region where damage occurs. In some cases, the shape of the added displacement field may be derived from strain localization analysis [13] and it is strain and damage dependent. Still, the enhancement is an approximation which does not consider distributed cracking, and thus requires a reasonably fine mesh for capturing distributed damage without the enhanced displacement fields.

It might be more attractive to combine continuum structural descriptions at different scales into a single analysis. The larger (global) scale aims at capturing the overall structural response with the help of a linear analysis, while the small (local) scale is enriched with a non-linear
material response that captures progressive damage and localization consistently. This
contribution investigates this last possibility. Among the various multiple scale analyses
techniques, a classical global-local method is implemented for this purpose. Attention is given
to the adaptation of the region where the small scale computation is performed as the damage
zone propagates. This two level method of analysis is illustrated in the case of a straight crack
propagating in a concrete plate. A second example is provided in order to show that the method
can handle curved crack propagation, which is very demanding in terms of mesh fineness even
when laboratory size specimens are considered.

## 2. MULTI-SCALE ANALYSES OF LARGE STRUCTURES

Over the years, a series of developments has been made for improving the efficiency of finite
element applications to large engineering problems. The simplest possible technique is the
uncoupled method which is widely used in engineering applications.

Global-local methods originated with Mote in 1971 [14] where finite element interpolants
were enriched with special functions. Early global-local structural schemes [15–18], focused on
the Rayleigh-Ritz method combined to the conventional finite element method over the entire
region of a structure, or to a local enhancement of the finite element description. Similar
contributions included the zooming method proposed by Hirai et al. [19, 20] and other authors
[21]. These methods usually involve refined meshes for the local regions containing stress
concentrations. Multiple scale finite element analyses have become more and more popular since
then. They can be classified into three different approaches: the superposition techniques,
multiple scale expansion methods, and domain decomposition techniques.

In the superposition technique, independently modelled structures (global and local meshes)
are considered. The method hinges on the hierarchical decomposition of the solution space into
global and local effects, and on the enforcement of solutions compatibility by prescribing
homogeneous boundary conditions at the global local interface. Among these is the composite
grid method (see e.g. References [22–24]) which provides a robust technique whose convergence
properties have been evaluated. The critical issues in the superposition method is the selection of
this global local interface. A rigorous mathematical analysis aimed at quantifying the ‘pollution’
effects of the localized phenomena on the global behaviour, and subsequently identifying the
optimal location of the interface has been carried out by Babuska et al. [25], and Fish and
Markolefas [26]. In Reference [26], an adaptive strategy has been devised to construct an
optimal discretization of the local and global solutions so that the local and global phenomena
of interest are resolved within the user specified accuracy.

Multiple scale expansion techniques are homogenization methods were field variables are
developed into asymptotic series and solved at each scale. They have been mainly implemented
in the field of composite materials where damage originates at a very small scale level, which is
not tractable within a single structure discretization. Some of these methods result into a non-
local formulation of damage that is similar to what will be used in the examples presented in this
paper [27] or may combine asymptotic expansions with global-local techniques [28]. A third,
smaller, scale level in which the composite microstructure is explicitly discretized, may also be
implemented in order to provide a more realistic description of local failure processes [29].

In domain decomposition techniques, the structure is divided into several sub-domains
connected to each other with interface elements (see for instance References [30–32]). Hence, the
discretization in each sub-domain is independent and can be decided depending on the type of material response expected (linear or non-linear) or on the state of stress. Such a technique has been also applied to hierarchical modelling of heterogeneous solids with a control of the quality of the numerical homogenization [33].

The above multi-scale methods provide quite robust computational tools. For the applications considered in this paper (large size concrete structures subjected to damage and cracking), the composite grid or domain decomposition methods seem for instance well suited. As the damage zones evolve, however, the domain in which the non-linear analysis is performed needs to be adapted. A step by step definition of the interface elements between local and global grids or sub-domains must be implemented.

In this contribution, we would like to focus on such a strategy for capturing a damage and cracked zone in a structure in the course of loading, i.e. on the adaptation of the size and location of the region where the non-linear computation is performed. We are going to select a simple zoom technique, the classical refined global-local method, due to Noor [16] and Mao and Sun [34] which avoids interface elements between the local and global region and the inherent adaptation of such a discretization. It is a first step enhancement of the uncoupled zoom method often used in industrial applications. This simplification is made at the price of a less robust multiple scale approach, in term of convergence properties especially. It should be emphasized that the foregoing strategy could be implemented similarly in more sophisticated multi-scale schemes.

### 3. ZOOM TECHNIQUE

#### 3.1. *Principle and damage zone tracking*

For the sake of illustration, consider the case of mode I crack propagation in plate (Figure 1(a)). The analysis and the decomposition of the problem into smaller sub-problems are performed as follows:

* In the first step, a linear analysis is performed on the coarse mesh (large scale). Approximate solutions of the displacements and stress fields are calculated. This is the initial solution.

* In the second step, the large finite elements, which should contain the fracture process zone and the crack tip, are taken for consideration. They are replaced with a refined mesh (Figure 1(b)), where the material follows a non-linear constitutive relation which captures the material degradation. Along the boundary of this mesh, the displacements obtained from the first step (or the global analysis) are enforced as boundary conditions. This computation forms the local solution.

* Because of the non-linear response of the local solution, the resulting forces acting at the external boundaries in the local solution are different from the forces acting along the edges ABCD of the large element in the initial solution. The corresponding reactions on the boundaries ABCD are applied as external forces in the global analysis, and a new global displacement field is obtained. At this stage, which is the third step, the large element ABCD has been removed from the mesh in the global, elastic, solution.

* Steps 2 and 3 are iterated until the forces at the external boundary of the local solution have stabilised from one iteration to the next one.
<figure><img src="image_2.png" /></figure>

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

<figure><img src="image_4.png" /><figcaption>Figure 1. (a) Global coarse mesh; (b) local refined region; and (c) zoom as a patch of finite elements.</figcaption></figure>

Several remarks should be made prior to dealing with crack propagation tracking:

* The use of the displacement field of the global analysis as a boundary condition for the local region is the simplest possibility. The enforcement of nodal forces (or distributed forces along the edge of element ABCD) might also have been chosen.

* The efficiency of this method relies on: (1) the detection of the inception of non-linear material response in the global solution; (2) a simple rule for tracking the evolution of the
zones where a local solution with a fine mesh is needed; and (3) the capacity to describe the
response of totally cracked elements in the case of friction or crack closure.

As far as the first item is considered, the location of the zones where a non-linear response
may occur is often controlled by the geometry of the structure and by the boundary conditions.
A reasonably accurate elastic analysis is often sufficient for the purpose of detecting where
cracking may initiate for instance. The inception of localized cracking in the local solution may
be triggered via a small notch if necessary (Figure 1(b)). It avoids unrealistic results where
damage is uniformly spread over all the elements in the local solution or cases where damage is
concentrated in the corners of the zoom due to spurious boundary effects.

For totally cracked elements, the issue of crack closure and friction between crack faces is
outside the scope of the present paper. It corresponds to the final stage of material degradation
locally. One may think to replace the crack with an interface element. Another possibility is to
keep the cracked element in the zoom instead of removing it and to rely on the local constitutive
response inside this region in order to capture crack closure or sliding. Consequently, the size of
the zoom region is bound to increase in the course of the development of the damage region.
In the examples of Section 4, the stiffness of elements containing a through crack will be set close
to zero.

Attention is focused on item 2, i.e. a simple method of tracking the evolution of damage (or
plastic) zones including crack tips in the fracture process zone. Let us assume that convergence is
achieved at a given load step. We check first for broken elements and look for the elements
which contain the tips of the damage zones at the scale of the global analysis. A broken element
contains a through crack, a band made of local elements where the material strength has been
exhausted completely. It corresponds to the cases where the material is entirely damaged,
according to a damage model, or has reached a critical strain in elasto-plasticity.
The new patch of global elements (or zoom) for the next load step is constructed as shown in
Figure 1(c). The local solution is preferably made of a patch of several large finite elements in
order to avoid boundary effects which may appear when such a zone is very close to the
boundary of the element in the global solution considered for zooming. For a straight crack
running parallel to the patch boundaries, the patch is made of the large finite element which
contains the damage zone and the crack tip, and the element sitting next to it in the propagation
direction observed during the previous loading step. If the tip of the fracture process zone is
located in element 2 on Figure 1(a), the local solution is formed with a patch made of elements 2
and 3. When element 2 is entirely cracked the patch is changed to element 3 and 4 and the
stiffness of element 2 is set close to zero. Note that according to this scheme, the patch should
never be entirely cracked after convergence. The size of the load steps should be controlled in
order to avoid such cases.

When the crack or damage zones are not straight or propagate in an arbitrary direction, the
patch is made of three large finite elements as shown in Figure 1(c). It is made of the element
which contains the damage zone and the crack tip and the surrounding elements in the direction
of propagation of the damage zone.

This process assumes that damage is located into small regions compared to the element size
in the global solution. In particular, this method fails when the damage zone extends over
several large finite elements without any crack. The patch of elements in the local solution ought
to be enlarged (or the mesh in the global solution might be modified) so as to encompass in the
local regions of the structure al the zones where material non-linearity occurs.
## 3.2. Formulation of the zoom technique

Consider a two-dimensional solid denoted $\Omega$. This solid is split into two parts: one has a linear elastic behaviour and is denoted $\Omega_G$; it is discretized with a relatively coarse mesh. The second part gathers the regions of the structure where material non-linearity occurs. It is denoted as $\Omega_L$ and it is discretized with a fine mesh. These regions may be plastic zones or in the present application damaged and cracked zones. The two zones are separated by the boundary $\Gamma_L$, and $\Omega = \Omega_G \cup \Omega_L$; $\Omega_G \cap \Omega_L = \Gamma_L$.

The equilibrium equations of the complete system $\Omega$ are expressed as

$$ \delta u^T \int_{\Omega} B^T \sigma \, d\Omega = \delta u^T F \quad (1) $$

where $\delta u$ is the variation of nodal displacement vector and $F$ is the nodal force vector applied at the boundary of the solid $\Omega$. The stress and displacement fields for this problem can be solved by superposition of the solutions of two problems (Figure 2). The first part of Equation (1) may be written as

$$ \delta u_G^T \int_{\Omega} B^T \sigma \, d\Omega = \delta u_G^T \int_{\Omega_G} B^T D_e B \cdot u_G \, d\Omega_G + \delta u_G^T \int_{\Omega_L} B^T \sigma_L \, d\Omega_L \quad (2) $$

where $\delta u_G$ refers now to the variation of the nodal displacements $u_G$ discretized with the coarse mesh in $\Omega_G$, $D_e$ is the elastic stiffness of the material in this region, and $\sigma_L$ is stress vector defined in the region $\Omega_L$. It is the solution of the following problem:

$$ \int_{\Omega_L} B_L^T \sigma_L \, d\Omega = \int_{\Gamma_L} N_L^T p_L \, dS \quad \text{with} \quad \sigma_L = D_{nl} B_L u_L \quad (3) $$

$u_L$ is the resulting displacement vector from the local fine mesh analysis, $N_L$ and $B_L$ are the shape function and its derivative according to the fine mesh respectively, and $D_{nl}$ is the constitutive stiffness matrix of the material in this region. In Equation (3), $p_L$ denotes the surface traction along $\Gamma_L$ (Figure 2).

We need now to define the boundary conditions in the two subproblems, the global and the local ones (Equations (2) and (3)). For the local one, pertaining to region of non-linear behaviour (Equation (3)), displacements are applied along the external boundary $\Gamma_L$. These

<figure><img src="image_5.png" /><figcaption>Figure 2. Superposition scheme.</figcaption></figure>
displacements denoted as $u_L$ are the projections of the displacements computed in the first subproblem (Equation (2)):

$$u_L(x, y) = N\delta_G \qquad (4)$$

where $N$ is the interpolation function in the global solution and $\delta_G$ is the vector of the generalized displacements at element nodes at the large scale. For the global problem (Equation (2)) the region $\Omega_L$ is substituted with equivalent applied forces. The forces denoted as $P_G$ in Figure 2 equilibrate the projection of the surface traction $p_L$ according to the interpolation of the coarse mesh:

$$P_G + \int_{\Gamma_L} N^T p_L \, dS = 0 \qquad (5)$$

Substitution of (5) and (3) into (2) yields

$$\int_{\Omega_G} B^T \sigma_G \, d\Omega_G = [F + P_G] \qquad (6)$$

The computation, and more specifically the iterative scheme between the global and local solutions, is carried out as follows:

1. An initial solution is obtained on the entire structure $\Omega$ assuming that the material is elastic and where $\Omega_L$ is replaced with large finite elements.

2. A local solution in $\Omega_L$ is obtained according to Equation (3) with the boundary conditions Equation (4).

3. The new pressure $P_G$ derived from the local solution and Equation (5) is applied as external forces in the global analysis performed now on $\Omega_G$.

4. Steps 2 and 3 are iterated until the variation of the interface traction between two iterations becomes small enough. The convergence criterion is

$$\frac{\|P_G^{i+1} - P_G^i\|}{\|P_G^i\|} \le e \qquad (7)$$

in which $i$ is the iteration number and $e$ is a given small tolerance; $e = 0.0005$ is used in computations and usually few iterations were needed.

It should be underlined that it cannot be proved that this simple iterative procedure converges in all instances and some more appropriate schemes such as those recalled in Section 2 might be preferred for a better robustness. Still, it possesses the advantage of being very simple to implement and quite versatile, since it does not depend on the constitutive relations implemented or on the local and global displacement field discretizations. Convergence will be checked numerically in one of the two examples carried out next.

## 4. EXAMPLES

Before presenting the illustrative examples, it is necessary to recall briefly the type of constitutive relations which have been implemented for the local FE analysis. We have chosen to use a scalar non-local damage model [3]. It should be emphasized that any non-linear constitutive relation could be implemented as well, provided that it is capable of describing failure and strain softening consistently. The zoom technique presented in Section 3 has been implemented in the
finite element package CAST3M. Two examples are presented. The first one deals with the wedge splitting test in the case of mode I cracking. The second example is the single-edge notched concrete specimen. Results of this method are compared with those performed with a global fine mesh, considered to be a reference solution.

## 4.1. Non-local damage model

In the reversible (elastic) domain, the stress-strain relation reads:

$$ \sigma = (1 - d)D_e : \varepsilon \qquad (8) $$

where $\varepsilon$ is the strain tensor, and $D_e$ is the elastic stiffness of the material, in the absence of damage. It will be assumed in the following that the material is initially isotropic. $d$ is the damage scalar variable which varies between 0 and 1. In the non-local version of this model, the growth of damage is controlled by the average $\bar{\xi}$ of an equivalent strain $\xi$ defined as

$$ \xi = \sqrt{\sum_{i=1}^{3} (\langle \varepsilon_i \rangle_{+})^2} \qquad (9) $$

where $\varepsilon_i$ are the principal strains and $\langle \rangle_+$ is the Macauley bracket. The average (non-local) equivalent strain is

$$ \bar{\xi} = \frac{1}{V_r(x)} \int_{\Omega} \psi(x-s) \xi(s) \, ds \quad \text{with} \quad V_r(x) = \int_{\Omega} \psi(x-s) \, ds \quad (10) $$

where $\Omega$ is the volume of the structure, $V_r(x)$ is the representative volume at point $x$, and $\psi(x-s)$ is the weight function [35]. The loading function is $f(\bar{\xi}, \kappa) = \bar{\xi} - \kappa$, where $\kappa$ is the hardening-softening variable. The contour $f(\bar{\xi}, \kappa) = 0$ defines the domain of reversible behaviour and the growth of damage is defined with classical Kuhn-Tucker conditions. The damage evolution laws are those defined by Mazars [36]. In particular, the tensile damage growth is

$$ d_t = 1 - \frac{\kappa_0(1 - A_t)}{\kappa} - \frac{A_t}{\exp(B_t(\kappa - \kappa_0))} \qquad (11) $$

Constants $A_t$, $B_t$ are model parameters. Overall, the constitutive relation has seven model parameters (in addition to the two elastic constants) whose range of variation is provided in Reference [3], for standard concrete.

## 4.2. Wedge splitting test

A plane stress model of the wedge splitting test is considered first (Figure 3). The horizontal displacements at nodes 1 and 2 are controlled, they are opposite and their absolute value is equal. First a reference fine mesh analysis is performed (Figure 3). A very dense mesh is generated in the fracture process zone, where crack propagation occurs. The interpolation of the displacements in the finite elements is linear. The model parameters are given in Table I. The rest of structure, outside the fracture process zone which is discretized with small rectangular elements, behaves elastically. In order to avoid some possible mesh bias, the finite element size in the process zone is small enough compared to the internal length given in Table I (at least one third of the internal length).

The zoom technique uses two finite element models. The global analysis (first step) employs a global coarse mesh as shown in Figure 4(a). The size of the finite element compared to the internal length is $h = 5l_c$ (larger than the FPZ width). The local analysis (second step) is
<figure><img src="image_6.png" /><figcaption>Figure 3. Wedge splitting test: finite element mesh for the reference analysis.</figcaption></figure>

<table><caption>Table I. Material properties in the wedge splitting test.</caption><tr><td colspan="4"><em>Concrete:</em></td></tr><tr><td>Young's modulus <i>E</i></td><td>32 000</td><td>MPa</td><td></td></tr><tr><td>Poisson's ratio <i>ν</i></td><td>0.2</td><td></td><td></td></tr><tr><td>Tensile strength <i>f</i><sub>t</sub></td><td>3.2</td><td>MPa</td><td></td></tr><tr><td>Compressive strength <i>f</i><sub>c</sub></td><td>35</td><td>MPa</td><td></td></tr><tr><td colspan="4"><em>Other model parameters:</em></td></tr><tr><td>κ<sub>0</sub></td><td>1 × 10<sup>4</sup></td><td></td><td></td></tr><tr><td><i>A</i><sub>c</sub></td><td>1.4</td><td><i>B</i><sub>c</sub></td><td>1500</td></tr><tr><td><i>A</i><sub>t</sub></td><td>1.0</td><td><i>B</i><sub>t</sub></td><td>10 000</td></tr><tr><td>β</td><td>1.06</td><td><i>l</i><sub>c</sub></td><td>20 cm</td></tr></table>

<figure><img src="image_7.png" /><figcaption>Figure 4. Zoom technique on the wedge splitting test: (a) global coarse mesh; and (b) local fine mesh.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5. Zoom technique on the wedge splitting test: load vs CMOD response.</figcaption></figure>

performed on a local region covered by a patch of 2 coarse elements since it is expected that the
crack is aligned with the mesh (Figure 1(c)). The refined mesh consists of 165 quadrilaterals in
every coarse element (Figure 4(b)). The element response in the local analysis follows the non-
local damage model (with the same parameters as in the reference solution). It should be noted
that the mesh density in the local analysis has been kept the same as in the fracture process zone
of the reference solution in order to assess the accuracy of the zoom technique itself,
independently from possible mesh effects. The global solution should provide a solution that is
as close as possible to the elastic response of the plate. This is required in order to approximate
the initial stiffness of the structure with a good accuracy for instance. This is the reason why we
refined the mesh near the points where the displacements are applied.

Figure 5 shows the plot of the applied load versus the crack mouth opening displacement (CMOD). We have plotted in the same figure the results obtained from the zoom technique (dashed curve) and the reference solution. A good agreement is observed between them. The reference solution looks slightly more brittle in the lower half of the softening regime.

The distribution of damage at different steps of calculation is plotted in Figure 6. Only the
region made of the patch of two elements numbered in Figure 4 is shown. Again, a good
agreement between the two calculations is observed. The damage zone propagates slightly faster
in the reference solution than in the zoom technique, which is consistent with a more brittle
response on Figure 5. Concerning the CPU time, it is divided by a factor higher than two with
the zoom technique compared to the reference solution.

4.3. *Numerical study of convergence*

Although the convergence of the present two level approach cannot be proved in a general
context, it is still feasible to examine convergence properties numerically. For this purpose, the
previous example is considered and we are looking for two characteristics: what should be the
minimum size of the zoom and what is the rate of convergence of the local-global iterative
approach. Note that separate convergence of the global and local iterative schemes alone is
already established: for the global calculation, it is rather obvious since it falls in the category of
elastic analyses with small deformations. Convergence of the local (non-linear) calculation has
been proved numerically on many instances (see e.g. References [3, 35]). The constitutive
<figure><img src="image_9.png" /><figcaption>Figure 6. Damage evolution in the wedge splitting test: (a) reference analysis; and (b) zoom technique.</figcaption></figure>

properties are the same as those in Table I, except for the internal length which is 22 cm (which will provide slightly different results compared to those in the previous section).

In order to investigate the influence of the size of the zoom, we consider a regular mesh for the wedge splitting test with an element size that is of the order of the internal length ($1.5l_c$ as shown in Figure 7(a)). The sizes of the zoom are patches of elements of the same height but with four gradually increasing width: 1, 3 and 5 times the size of coarse elements. This width is measured in a direction that is perpendicular to the direction of propagation of the fracture process zone. In the zoom, each element of the global mesh is subdivided into a regular grid of 5 × 5 elements, whose size is much smaller than the internal length in order to avoid discretization errors in the zoom as much as possible since the size of the finite element is less than one third the internal length. Figure 7(b) shows the patches of 1, 3 and 5 elements respectively. Figure 8 shows the load vs CMOD responses obtained for the different patches of elements and the corresponding damage zones. It is observed in Figure 8(a) that the peak load in particular, increases slowly as the size of the zoom increases. Recall that the material at the global scale is linear elastic. Material non-linearity is confined to occur in the zoom. If the zoom is not wide enough, the fracture process zone will remain confined within a region that is too small. As a consequence, the damage zone for the smallest patch is very different from the ones obtained with wider patches that are quasi-identical (Figure 8(b)). This constraint produces a concentration of the non-linear response in a small region and results into a smaller peak load, as if the internal length was decreasing. As the fracture process zone propagates and completely cracked large elements are removed, the influence of this constraint becomes smaller and the computed curves become similar. Convergence on the load vs CMOD curve can be considered to be achieved when the zoom encompasses 5 elements at the macro-scale. Note that for a proper description of the damage zone, a sufficient width of the zoom is three elements. A ratio of this order has been used in the previous section.
<figure><img src="image_10.png" /><figcaption>Figure 7. Wedge splitting test: convergence study.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 8. Convergence on the load vs CMOD curves (a) and of the damage zones (b).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 9. Evolution of the relative error with the number of iterations for the different patch sizes.</figcaption></figure>

Figure 9 shows the evolution of the residual error defined in Equation (7) as a function of the iteration number. This is the relative norm of the applied forces in the global analyses from one iteration to the next one. Patches of width 3 and 5 elements have been considered. The point at which this error is computed corresponds to a CMOD = 6.5 mm. It is observed that the error decreases rapidly. The rate of convergence (opposite of the slope of the curve) is decreasing however with the residual error. Nevertheless an acceptable value (less that 0.05%) is reached within a limited number of iterations. Oscillations are less pronounced as the size of the patch of elements in the local calculation increases and the average rate of convergence increases. This result shows that with a sufficiently large local zone, the method can be expected to provide a reasonably accurate approximate solution, at least better that the uncoupled technique.

4.4. Single-edge notched specimen

The second example is a single-edge notched concrete beam subjected to an anti-symmetric four
point loading. The anti-symmetric loading results in a curved crack path in the experiments
performed on a smaller specimen by Schlangen [37]. The crack starts from the right corner of the
notch and ends to the right of the lower right load platen.

The geometry and loading conditions are given in Figure 10(a). The notch depth in the centre of the beam is 40 cm and the notch width is 5 cm. Note that the size of the specimen has been very much enlarged compared to the original experiments. The load P is applied on the loading apparatus so that the point load close to the notch is P₁ = 10/11P and the point load near the beam end is P₂ = 1/11P. This loading apparatus will be considered as a rigid body in the calculation, including the load platens.

A plane stress calculation has been performed with a linear interpolation of the displacements. Figure 10(b) shows the finite element mesh used for the reference solution. The model parameters are those of the wedge splitting test. An indirect displacement control procedure has been used. The control parameter $\delta_2$ defined in Figure 10(a) is the relative vertical displacement of the neutral axis between the left support and the middle of the beam. At the middle of the beam, the relative vertical displacement is computed as the average vertical displacements of the two extreme (top and bottom) faces of the specimen.

For the zoom technique, the global coarse mesh (first step) is presented Figure 11(a). The size of the finite element is the same as in the previous example. The local analysis (second step) is
<figure><img src="image_13.png" /><figcaption>Figure 10. Single edge notched beam: (a) geometry and loads (dimension in cm); and (b) finite element mesh for the reference analysis.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11. Zoom technique: (a) global coarse mesh; and (b) local fine mesh.</figcaption></figure>

performed on a local region covered by four coarse elements, the element which contains the tip of the damaged zone and those around it in the direction of propagation of the damaged zone (Figure 11(b)). Again, the mesh density in the local analysis is the same as in the reference solution in the centre of the beam.
<figure><img src="image_15.png" /><figcaption>Figure 12. Single edge notched beam: load vs $\delta_2$ response.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13. Damage evolution in the single edge notched beam: (a) reference analysis; and (b) zoom technique.</figcaption></figure>

Figure 12 shows the plot of the applied load vs the deformation parameter $\delta_2$. The error is slightly larger than in the previous example, but still acceptable and the gain in CPU time is similar to that in the wedge splitting test. The evolution of damage at different steps of the calculations (indicated in Figure 12) is shown in Figure 13. In the first stage of the fracture process, damage is initiated at the right corner of the notch. Then, the damage zone at the notch continues to grow and follows a path which ends to the right of the lower right load platen eventually. We can see that the evolution of damage is slightly faster in the reference analysis than in the zoom technique. Note also that the width of the damage zone in the reference analysis is also slightly larger than in the zoom technique.
<figure><img src="image_17.png" /><figcaption>Figure 14. Final damage distribution in the single edge notched beam: (a) reference analysis; and (b) zoom technique.</figcaption></figure>

The final damage distribution is given in Figure 14. A curved damage zone is found in both cases, which is consistent again with the experiment on a smaller specimen. The damage distribution in Figure 14(b) results from those obtained in the successive patches of elements composing the local solution in the zoom technique. This is the reason why the top part of the distribution (represented in light grey) does not connect with the distribution in the bottom part. The top patch has been removed after the crack has crossed it. This mismatch indicates also one of the reasons why the load displacement curve according to the zoom technique softens more severely. It is often observed in finite element solutions with non-local models that the damage zone widens in the course of the fracture process, even after a crack has formed in its center. Figure 14(a) shows that it occurs in the reference solution. In the zoom technique however, this cannot occur because elements which are crossed by the crack are removed and the growth of the damage zone inside them is stopped after their removal. Hence damage is more distributed in the reference solution compared to the zoom technique. As a consequence the load deflection curve softens less because the energy dissipation is spread over a larger zone in the reference solution compared to the zoom solution.

## 5. CONCLUDING REMARKS

A simple strategy for technique for finite element analysis of large concrete structures with crack and damage propagation based on an adaptive zoom has been presented. This method has been applied to the failure analysis of concrete structures with the help of a non-local damage model. A global solution is obtained using a coarse grid, then the detailed stress and damage distribution near the crack (local region) are obtained iteratively by zooming on that area, refining the model and using the displacements from the coarser model as an input for refined mesh. The zoom comprises the element in which the tip of the damage zone is located and the
surrounding elements. It follows the propagation of the damage zone, and thus avoids any additional criterion for its propagation as opposed to fracture analyses.

In the presented applications, completely cracked elements are removed from the global and local meshes but they could be kept in the local analysis in order to capture cyclic (loading–unloading) effect or a subsequent different state of loads involving shear for instance. This requires of course that the constitutive relations are capable of capturing such stress histories.

Compared to the internal length in the non-local model, the width of the zoom should be at least 3–5 times the internal length. The discretization inside the zoom ought to meet the standard requirements for non-local analyses.

The damage propagation tracking is rather general in the sense that it can be implemented with various continuum models without changing the form of the discrete fields. It should be emphasized that the zoom technique implemented in this contribution, however, is simple and might present some computational weaknesses. More sophisticated multi-scale schemes might be preferred for a better robustness as far as convergence is concerned especially.

## REFERENCES

1. Badel P, Chavant C, Ghavamian S, Voldoire F. An industrial RC structure—application using damage mechanics, convergence difficulties and calculation technique. *CD-ROM Proceedings of ECCOMAS 2000*, Barcelona, 2000.

2. Carmeliet J. Optimal estimation of gradient damage parameters from localisation phenomena in quasi-brittle materials. *International Journal of Mechanics and Cohesive Frictional Materials* 1999; **4**:1 16.

3. Pijaudier-Cabot G, Mazars J, Pulikowski J. Steel concrete bond analysis with non-local damage. *Journal of Structural Engineering ASCE* 1991; **117**:862 882.

4. Belytschko T, Tabbara M. H-adaptive finite element methods for dynamic problems. *Computer and Structures* 1991; **55**:237 247.

5. Diez P, Arroyo M, Huerta A. Adaptivity based on error estimation for viscoplastic softening materials. *International Journal of Mechanics of Cohesive Frictional Materials* 2000; **5**:87 112.

6. Ortiz M, Quigley JJ. Adaptative mesh refinement in strain localisation problems. *Computer Methods in Applied Mechanics and Engineering* 1991; **90**:781 804.

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

8. Belytschko T, Fish J, Engelmann BE. A finite element with embedded localisation zones. *Computer Methods in Applied Mechanics and Engineering* 1988; **70**:59 89.

9. Jirasek M. Embedded crack models for concrete fracture. *In Computational Modelling of Concrete Fracture (EUROC)*, de Borst R, Bicanic N, Mang H, Meschke G (eds). Balkema: Rotterdam, 1998; 291 300.

10. Wells GN, Sluys LJ. Discrete analysis of localisation in three dimensional solids. *CD-ROM Proceedings of ECCOMAS 2000*, Barcelona, 11-14 Sept.

11. 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.

12. Moes N, Belytschko T. Extended finite element method for cohesive crack growth. *Engineering in Fracture Mechanics* 2002; **69**:813 833.

13. Patzak B, Jirasek M. Process zone resolution with extended finite elements. *Engineering Fracture Mechanics* 2003, in press.

14. Mote CD. Global local finite element. *International Journal for Numerical Methods in Engineering* 1971; **3**:565 574.

15. Dong SB. Global local finite element methods. In *State of the Art Surveys on Finite Element Technology*, Noor AK, Pilkey WD (eds). American Society of Mechanical Engineers, New York, 1983; 451 474.

16. Noor AK. Global local methodologies and their application to non-linear analysis. *Finite Element Analysis and Design* 1986; **2**:333 346.

17. Noor AK. 1986. Reduction method for the nonlinear analysis of symmetric anisotropic panels. In *Finite Element Methods for Nonlinear Problems*, Bergan PG, Bathe KJ, Wunderlich W (eds). Springer: Berlin, 389 407.

18. Noor AK. On making large non-linear problems small. *Computer Methods in Applied Mechanics and Engineering* 1982; **34**:955 985.

19. Hirai I, Wang BP, Pilkey WD. An efficient zooming method for finite element analysis. *International Journal for Numerical Methods in Engineering* 1984; **20**:1671 1683.

20. Hirai I, Uchiyama Y, Mizuta Y, Pilkey WD. An exact zooming method. *Finite Element Analysis and Design* 1985; **1**:61 69.

21. Cook RD. *Concepts and Applications of Finite Element Analysis*. Wiley: New York, 1974.
22. Fish J, Belsky V, Pandheeradi M. Composite grid method for hybrid systems. *Computer Methods in Applied Mechanics and Engineering* 1996; **135**:305 327.

23. Fish J, Suvorov A, Belsky V. Hierarchical composite grid method for global local analysis of laminated composite shells. *Applied Numerical Mathematics* 1997; **23**:231 248.

24. McCormick SF, Thomas JW. The fast adaptive composite grid (FAC) method for elliptic equations. *Mathematics with Computers* 1986; **46**:439 456.

25. Babuska I, Strouboulis S, Upadhay CS, Gangaraj SK. A posteriori estimation and adaptive control of the pollution-error in the h-version of the finite element method. Technical Note BN-1175, Institute for Physical Science and Technology, University of Maryland, College Park, MD, 1995.

26. Fish J, Markolefas S. Adaptive global local refinement strategy based on the interior error estimates of the h-method. *International Journal for Numerical Methods in Engineering* 1994; **37**:827 1031.

27. Fish J, Yu Q. Two-scale damage modelling of brittle composites. *Composite Science and Technology* 2001; **61**: 2215 2222.

28. Fish J, Shek KL. Multiscale analysis for composite materials and structures. *Composites Science and Technology* 2000; **60**:2547 2556.

29. Lee K, Moorthy S, Ghosh S. Multiple scale computational model for damage in composite materials. *Computer Methods in Applied Mechanics and Engineering* 1999; **172**:175 201.

30. Ladevze P, Dureisseix D. Une nouvelle stratgie de calcul micro/macro en mcanique des structures. *Comptes Rendus des Sciences de l' Academie des Sciences Paris* t.327, Srie IIb, 1999; 1237 1244.

31. Ladevze P, Nouy A. Une stratgie de calcul multichelle avec homognisation en espace et en temps. *Comptes Rendus Mecanique* 2002; **330**:683 689.

32. Le Tallec P. *Domain Decomposition Methods in Computational Mechanics*. North-Holland: Amsterdam, 1994.

33. Oden JT, Vermaganti K, Moes N. Hierarchical modelling of heterogeneous solids. *Computer Methods in Applied Mechanics and Engineering* 1999; **172**:3 25.

34. Mao KM, Sun CT. A refined global local finite element analysis method. *International Journal for Numerical Methods in Engineering* 1991; **32**:29 43.

35. Pijaudier-Cabot G, Bazant ZP. Nonlocal damage theory. *Journal of Engineering Mechanics ASCE* 1987; **113**: 1512 1533.

36. Mazars J. Application de la mcanique de l'endommagement au comportement non linaire et la rupture du bton de structure. *PhD Thesis*, 1984, Universit Paris VI, Paris, France (in French).

37. Schlangen E. Experimental and numerical analysis of fracture process in concrete. Dissertation, Delft University of Technology, Delft, the Netherlands, 1993.