## Abstract

Progressive fracture in quasi-brittle materials such as concrete, rocks, soils, is often treated as strain softening in continuum damage mechanics. Such constitutive relations favour spurious strain localization and ill-posedness of boundary value problems, and call for some kind of regularization. In the present work, two different approaches are presented: a partially regularized local damage model that adjusts the softening part of a stress-strain law depending on the size of the element, and a fully regularized non-local damage model that introduces the characteristic length as an additional material parameter controlling the size of the fracture process zone.

In addition, the strain softening of such models usually results in highly complex structural responses, including the snap-back type, and thus in this work we will be discussing non-linearity associated to damage modelling, and a global arc-length method for tracing the equilibrium path will be exposed.

Furthermore, in the context of non-local damage models it is crucial to work with fine spatial discretizations at the damage progress zone, so that elements are smaller than the characteristic length. In this regard, a mesh-adaptive technique has been implemented with the purpose of enhancing the efficiency of the numerical analysis.

Finally, two classical examples, the three-point bending test, and the single-edge notched beam test, are performed in order to analyse the mesh objectivity of the implemented integral-type non-local damage model, and assess the strengths and limitations of the mesh-adaptive procedure.

# Resum

L'evolució de la fractura en materials quasi-fragils com el formigó, les roques, o els sols, sovint és tractada amb estovament de deformacions dins el marc de la mecanica del dany continu. Aquestes relacions constitutives afavoreixen la localització fictícia de deformacions i el mal plantejament del les equacions i condicions de contorn del problema i, per tant, necessiten algun tipus de regularització. En aquest treball es plantegen dues metodologies diferents: un model de dany local parcialment regularitzat que ajusta la part d'estovament de la llei tensió-deformació en funció de la mida de l'element, i un model de dany no local plenament regularitzat que introdueix la longitud característica com un parametre del material addicional, que controla l'extensió de la zona de fractura.

D'altra banda, l'estovament de deformació d'aquests models normalment porta associades respostes estructurals complexes, incloent les de tipus snap-back i, per aixo, en aquest treball també exposarem les conseqüencies de la no-linealitat associada als problemes de dany, i presentarem un metode d'arc global per traçar el camí d'equilibri de la solució.

A més, en el marc dels models de dany no locals, és crucial treballar amb discretitzacions de l'espai prou fines a les zones de progrés del dany, per tal que els elements siguin menors que la longitud característica. En aquest sentit, s'ha implementat un tecnica de refinament de malla adaptatiu amb l'objectiu de millorar l'eficiencia dels analisis numerics.

Finalment, s'han dut a terme dos exemples classics, l'assaig d'una biga flectada per tres punts, i l'assaig d'una biga amb flexió no simetrica per quatre punts, amb l'objectiu de provar l'objectivitat de malla del model de dany no local, i analitzar els punts forts i limitacions del metode de refinament de malla adaptatiu.

# 1 Introduction and Objectives

Failure of quasi-brittle materials, particularly in geo-materials such as concrete, soils and rocks, is preceded by a gradual development of a non-linear fracture process zone and a localization of strain. Realistic failure analysis of such quasi-brittle structures requires the consideration of progressive damage due to micro-cracking, which is usually modelled by a constitutive law with strain softening. This typically results in highly non-linear structural responses and so efficient non-linear solvers based on arc-length control are needed for the numerical simulations.

If the damage parameter depends only on the strain state at the point under consideration, the boundary value problem becomes ill-posed and, as a result, numerical simulations exhibit a pathological mesh dependence and the energy consumed by the fracture process tends to zero as the mesh is refined.

Simple remedies based on the adjustment of the softening part of the stress-strain law, depending on discretization parameters like the size and type of the finite elements, try to enforce the correct energy dissipation by keeping the product of the width of the process zone and the local dissipation energy constant.

Advanced regularization methods introduce an additional material parameter, the characteristic length, which is related to the size and spacing of inhomogeneities and controls the width of the fracture process zone. These techniques, including non-local averaging methods, need sufficiently fine spatial discretizations in the process zone to resolve narrow bands of highly localized strains. In this regard, mesh-adaptive techniques can increase the efficiency of the analysis, automatizing the process.

In this context, the main objective of this work is to implement a robust isotropic damage model for the numerical modelling of geo-materials failure, that guarantees mesh objectivity of the solution. Furthermore, we will need to implement an incremental-iterative strategy based on the arc-length control, in order to deal with the strong non-linearities derived from the damage progression. Finally, an adaptive-mesh refinement technique would be a very attractive tool as a means to enhance the efficiency of the damage simulations.

This work is organized as follows. First of all, damage mechanics theory is introduced in the second chapter, stressing the strain localization phenomenon of classical local damage models and discussing two alternatives for regularizing the problem: the partially regularized local damage model, and the integral-type non-local damage model.

In the third chapter we will be reviewing all those tools that the damage model requires for its numerical implementation. This includes, an overview of the finite element method, a discussion on the non-linearity associated to damage problems, presenting a global arc-length strategy, and an introduction to mesh-adaptive techniques.

Finally, in the fourth chapter, two classical numerical examples in the damage mechanics field are performed in order to determine whether the non-local damage model is a robust method to model failure of quasi-brittle materials, and to test the implemented mesh-adaptive technique.

# 2 Damage Mechanics

## 2.1 Introduction

Many engineering materials subjected to unfavourable mechanical and environmental conditions undergo micro-structural changes that decrease their strength. Since these changes impair the mechanical properties of these materials, the term damage is generally used.

This effect is particularly relevant and difficult to predict in brittle or quasi-brittle materials such as concrete, rocks, mortar or other geo-materials. For instance, concrete is a highly heterogeneous, anisotropic, brittle material with a very complex non-linear mechanical behaviour due to the occurrence of localization of deformation. This localization of deformation can appear as cracks, if cohesive properties are dominant, or as shear zones, if frictional properties prevail. As a result of strain localization, material softening takes place and a significant reduction of the material stiffness during cyclic loading occurs. A good understanding of the mechanism of the formation of localization is of crucial importance since it acts as a precursor of fracture and failure.

In order to describe the behaviour of quasi-brittle materials, different approaches have been developed in the the last decades: continuum models within fracture mechanics [1], continuum damage mechanics [2,3], plasticity [4,5], coupled damage and plasticity [6,7], micro-plane theory [8,9], discrete models using a lattice approach [10,11] and the discrete element method (DEM) [12,13], to name but a few.

In the present work we have focused our efforts on developing an isotropic continuum damage model for geo-materials that works in a small deformation regime.

Particularly, this chapter will be devoted to present some basic notions on continuum damage mechanics, stressing the most relevant aspects of the implemented damage model; and then expose the strain localization phenomenon that arises in these kind of simulations, introducing two possible approaches to avoid the problems of localization.

## 2.2 Basic Concepts on Continuum Damage Mechanics

Continuum damage mechanics is a branch of continuum mechanics that describes the progressive loss of material integrity due to the propagation and coalescence of micro-cracks, micro-voids, and similar defects. These changes in the micro-structure lead to an irreversible material degradation, characterized by a loss of stiffness that can be observed on the macro-scale.

The term “continuum damage mechanics” was first used by Hult [14], but the concept of damage had already been introduced by Kachanov in 1958 in the context of creep rupture [15], and further developed by Rabotnov [16], Hayhurst [17], and Leckie and Hayhurst [18].

In the pioneering work of Kachanov [15] the concept of effective stress was introduced, and by using continuum damage he solved problems related to creep in metals. Rabotnov [16] gave the problem physical meaning by suggesting that the reduction of the sectional area was measured by means of the damage parameter. Thermodynamic formalism involved in the irreversible process of damage was developed by Lemaitre and Chaboche [19], and other important contributions to our knowledge about damage mechanics include: Mazars [20], Mazars and Pijaudier-Cabot [21], Chaboche [22], Simo and Ju [23], Ju [24,25], Oliver, Oller and Oñate [26], Oliver and Oller [27,28], Cervera and Chiumenti [29,30], etc.

The main objective of standard Continuum Damage Mechanics is to propose a continuum-mechanics based framework that allows to characterize, represent and model at the macroscopic scale the effects of distributed defects and their growth on the material behaviour. In order to fully understand these concepts, first it is important to make clear some basic notions of damage mechanics.

### 2.2.1 Description of classical uniaxial damage theory

Damage models work with certain internal variables that characterize the density and orientation of micro-defects. To introduce its basic concepts and understand the basis of damage mechanics, it is easier to begin with a uniaxial stress case.

 Figure 1: Idealized material for the description of the uniaxial damage theory.

For the present purpose, the material is idealized as a bundle of fibers parallel to the loading direction (Figure 1). Initially, all the fibers respond elastically, and the stress is carried by the total cross section of all fibers ${\textstyle S}$ (Figure 2.1). As the applied strain is increased some fibers start breaking (Figure 2.2). Each fiber is assumed to be perfectly brittle, meaning that the stress in the fiber drops down to zero immediately after a critical strain level is reached. However, since the critical strain can differ from one fiber to another, the effective area ${\textstyle {\bar {S}}}$ (the area of unbroken fibers that can still carry stress) will decrease gradually from ${\textstyle {\bar {S}}=S}$ to ${\textstyle {\bar {S}}=0}$. Of course, when the applied force diminishes (Figure 3.2), the affected fibers remain broken and so the damaged area of the specimen is irrecoverable.

 Figure 2: Scheme of a uniaxial damage model through a monotonic loading process.
 Figure 3: Scheme of a uniaxial damage model through a non-monotonic loading process.

It is important to make a distinction between the nominal stress ${\textstyle \sigma }$, defined as the force per unit initial area of the cross section, and the effective stress ${\textstyle {\bar {\sigma }}}$, defined as the force per unit effective area. The nominal stress enters the Cauchy equations of equilibrium on the macroscopic level, while the effective stress is the "actual" stress acting in the material micro-structure. From the condition of equivalence, ${\textstyle \sigma S={\bar {\sigma }}{\bar {S}}}$, we obtain:

 ${\displaystyle \sigma ={\frac {\bar {S}}{S}}{\bar {\sigma }}}$
(1)

The ratio of the effective area to the total area, ${\textstyle {\bar {S}}/S}$, is a scalar characterizing the integrity of the material. Within the classical approach, a very simple measure of the damage amplitude in a given plane is obtained by measuring the area of the intersection of all defects with that plane. Thereby, we can define the damage parameter as:

 ${\displaystyle d=1-{\frac {\bar {S}}{S}}={\frac {S-{\bar {S}}}{S}}={\frac {S_{d}}{S}}}$
(2)

where ${\textstyle S_{d}=S-{\bar {S}}}$ is the damaged part of the area. An undamaged material is characterized then by ${\textstyle {\bar {S}}=S}$, i.e., by ${\textstyle d=0}$. Due to propagation of micro-defects and their coalescence, the damage variable grows and at the late stages of degradation process it approaches asymptotically the limit value ${\textstyle d=1}$, corresponding to a complete damaged material with effective area reduced to zero.

In the simplest version of the presented scheme, each fiber is supposed to remain linear elastic up to the strain level at which it breaks. Consequently, the 1D effective stress ${\textstyle {\bar {\sigma }}}$ is related to the elastic strain of the material ${\textstyle \varepsilon }$ by the uniaxial Hooke's law:

 ${\displaystyle {\bar {\sigma }}=E\varepsilon }$
(3)

where E is the elastic modulus of the undamaged material. Combining (1), (2) and (3), it is easily seen that the constitutive law for the nominal stress ${\textstyle \sigma }$ takes the form:

 ${\displaystyle \sigma =(1-d)E\varepsilon }$
(4)

For the uniaxial model formulation, equation (4) must be complemented with the damage evolution law, which can be characterized by the dependence between the damage variable and the applied strain:

 ${\displaystyle {\begin{array}{cc}d=g(\varepsilon )&\quad 0\leq d\leq 1\end{array}}}$
(5)

Function ${\textstyle g}$ affects the shape of the stress-strain diagram and can be directly identified from a uniaxial test.

The evolution of the effective stress, damage variable, and nominal stress in a material that remains elastic up to the peak stress is shown in Figure 2. This description is valid only for monotonic loading by an increasing applied strain ${\textstyle \varepsilon }$. When the material is first stretched up to a certain strain level ${\textstyle \varepsilon _{1}}$ that induces damage ${\textstyle d_{1}=g(\varepsilon _{1})}$ and then the strain decreases (Figure 3), the damaged area remains constant and the material responds as an elastic material with a reduced Young's modulus ${\textstyle E_{1}=(1-d_{1})E}$. This means that, during unloading and partial reloading, the damage variable in (4) must be evaluated from the largest previously reached strain and not from the current strain ${\textstyle \varepsilon }$. This means that it is crucial to introduce an internal state variable ${\textstyle r}$ characterizing the maximum strain level reached in the previous history of the material up to a given time ${\textstyle t}$:

 ${\displaystyle r(t)=\max \left\{r_{0},\max _{\tau \leq t}\varepsilon (\tau )\right\}}$
(6)

The expression implies that ${\textstyle r(t)\geq r_{0}}$, where ${\textstyle r_{0}}$ is the so-called damage threshold, a material parameter that represents the value of strain at which damage starts. In this formula, ${\textstyle t}$ is not necessarily the physical time (it can be any monotonically increasing parameter controlling the loading process). The damage evolution law (5) is then rewritten in the form:

 ${\displaystyle d=g(r){\hbox{ with }}\left\{{\begin{array}{cc}g(r)=0&{\hbox{if }}r=r_{0}{\hbox{ }}\\0r_{0}{\hbox{ }}\end{array}}\right.}$
(7)

which remains valid not only during monotonic loading but also during unloading and reloading. The evolution of the effective stress, damage variable, and nominal stress in a non-monotonic test is shown in Figure 3. Note that, upon a complete removal of the applied force, the strain returns to zero (due to elasticity of the unbroken fibers), i.e., the pure damage model does not take into account any permanent strains. Nevertheless, the material state will be different from the initial virgin state because, according to (6) and (7), when the state variable ${\textstyle r}$ becomes greater than ${\textstyle r_{0}}$, the damage variable will not be zero again, thus the stiffness and strength mobilized in a new tensile loading process will be smaller than their initial values. Therefore, we can say that the loading history is always reflected by the value of the internal state variable ${\textstyle r}$.

To gain further insight, we can rewrite the constitutive law (4) in the form ${\textstyle \sigma =E_{sec}\varepsilon }$ where ${\textstyle E_{sec}=(1-d)E}$ is the apparent or damaged modulus of elasticity. Furthermore, instead of defining the variable ${\textstyle r}$ like (6), we can introduce a loading function ${\textstyle f(\varepsilon ,r)=\varepsilon {-}r}$ and postulate the loading-unloading conditions in the Kuhn-Tucker form:

 ${\displaystyle {\begin{array}{ccc}f\leq {0};&{\dot {r}}\geq {0};&{\dot {r}}f=0\end{array}}}$
(8)

The first condition indicates that ${\textstyle r}$ can never be smaller than ${\textstyle \varepsilon }$, while the second condition means that ${\textstyle r}$ cannot decrease. Finally, according to the third condition, ${\textstyle r}$ can grow only if the current values of ${\textstyle \varepsilon }$ and ${\textstyle r}$ are equal.

At this point, we can already summarize the basic components of the uniaxial damage theory:

• The stress-strain law in the secant format:

 ${\displaystyle \sigma =E_{sec}\varepsilon }$
(9)

with:

 ${\displaystyle E_{sec}=(1-d)E}$
(10)
• The law governing the evolution of the damage variable (7)
• The damage criterion, comprising:
 ${\displaystyle f(\varepsilon ,r)=\varepsilon {-}r}$
(11)

specifying the elastic domain ${\textstyle {\mathcal {E}}_{r}=\{\varepsilon |f(\varepsilon ,r)<0\}}$, i.e., the set of states for which damage does not grow.

### 2.2.2 A simple isotropic damage theory

The simplest extension of the uniaxial damage theory to general multiaxial stress states is achieved by the isotropic damage model with a unique scalar variable. Isotropic damage models are based on the simplifying assumption that the stiffness degradation is isotropic, i.e., the stiffness moduli corresponding to different directions decrease proportionally, independently of the direction of loading. Since an isotropic elastic material is characterized by two independent elastic constrains, a general isotropic damage model should deal with two damage variables. The model with a single variable makes use of the additional assumption that the relative reduction of all the stiffness coefficients is the same, in other words, that the Poisson's ratio is not affected by damage. In this regard, the damaged constitutive tensor is expressed as:

 ${\displaystyle {\boldsymbol {E}}_{sec}=(1-d){\boldsymbol {E}}}$
(12)

where ${\textstyle {\boldsymbol {E}}}$ is the elastic constitutive tensor of the intact material, and ${\textstyle d}$ is the damage variable. In the present context, ${\textstyle {\boldsymbol {E}}_{sec}}$ is the secant constitutive tensor that relates the total strain ${\textstyle {\boldsymbol {\varepsilon }}}$ to the total stress ${\textstyle {\boldsymbol {\sigma }}}$, according to:

 ${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {E}}_{sec}:{\boldsymbol {\varepsilon }}=(1-d){\boldsymbol {E}}:{\boldsymbol {\varepsilon }}}$
(13)

One can clearly notice that (12) is a generalization of (10), and (13) is a generalization of (9) and (4). Furthermore, equation (13) can alternatively be written as:

 ${\displaystyle {\boldsymbol {\sigma }}=(1-d){\boldsymbol {\bar {\sigma }}}}$
(14)

which is the multidimensional generalization of (1), and where ${\textstyle {\boldsymbol {\bar {\sigma }}}}$ is the effective stress tensor defined as:

 ${\displaystyle {\boldsymbol {\bar {\sigma }}}={\boldsymbol {E}}:{\boldsymbol {\varepsilon }}}$
(15)

Similar to the uniaxial case, we introduce a loading function ${\textstyle f}$ specifying the elastic domain and the states at which damage grows. The loading function now depends on the strain tensor ${\textstyle {\boldsymbol {\varepsilon }}}$, and on a variable ${\textstyle r}$ that controls the evolution of the elastic domain. Physically. ${\textstyle r}$ represents a scalar measure of the largest strain level ever reached in the history of the material. States for which ${\textstyle f({\boldsymbol {\varepsilon }},r)<0}$ are supposed to be below the current damage threshold. Damage can grow only if the current state reaches the boundary of the elastic domain, i.e., when ${\textstyle f({\boldsymbol {\varepsilon }},r)=0}$. Essentially, we can postulate the damage criterion for a multiaxial isotropic damage model with:

 ${\displaystyle f({\boldsymbol {\varepsilon }},r)=\varepsilon _{eq}({\boldsymbol {\varepsilon }})-r}$
(16)

where ${\textstyle \varepsilon _{eq}}$ is the equivalent strain, i.e., a scalar measure of the strain level, and ${\textstyle r}$ is the largest value of equivalent strain calculated in the previous deformation history of the material up to its current state. In this regard, (6) can now be generalized to:

 ${\displaystyle r(t)=\max \left\{r_{0},\max _{\tau \leq t}\varepsilon _{eq}(\tau )\right\}}$
(17)

An important advantage of an isotropic damage model is that the stress evaluation algorithm is usually explicit and there is no need to use an iterative solution for non-linear equations.

Thereby, for a particular strain increment, the corresponding stress is obtained by computing the current value of equivalent strain, updating the maximum previously reached equivalent strain and the damage variable, and reducing the effective stress according to (13). In essence, one must follow the scheme of Table 1:

 Input data for time ${\displaystyle t+1}$: ${\displaystyle {\boldsymbol {E}}}$, ${\displaystyle {\boldsymbol {\varepsilon }}^{t}}$, ${\displaystyle \Delta {\boldsymbol {\varepsilon }}}$, ${\displaystyle r^{t}}$ 1. Determine new strain vector: ${\textstyle {\boldsymbol {\varepsilon }}^{t+1}={\boldsymbol {\varepsilon }}^{t}+\Delta {\boldsymbol {\varepsilon }}}$ 2. Evaluate effective stresses: ${\textstyle {\boldsymbol {\bar {\sigma }}}^{t+1}={\boldsymbol {E}}:{\boldsymbol {\varepsilon }}^{t+1}}$ 3. Compute equivalent strain from the new strain vector: ${\textstyle \varepsilon _{eq}^{t+1}}$ 4. Update ${\textstyle r}$ with (17): If ${\textstyle (\varepsilon _{eq}^{t+1}>r^{t})\Rightarrow r^{t+1}=\varepsilon _{eq}^{t+1}}$ 5. Update damage variable: ${\textstyle d^{t+1}=g(r^{t+1})}$ 6. Compute stresses: ${\textstyle {\boldsymbol {\sigma }}^{t+1}=(1-d^{t+1}){\boldsymbol {\bar {\sigma }}}^{t+1}}$

#### 2.2.2.1 Equivalent strain forms

To some extent, the expression defining the equivalent strain plays a role similar to the yield function in plasticity, because it directly affects the shape of the elastic domain (Figure 4).

 Figure 4: 3D elastic domain for a generic equivalent strain.

There are numerous forms of equivalent strain in the literature, but the simplest choice is to define it as the Euclidean norm of the strain tensor:

 ${\displaystyle \varepsilon _{eq}=\|{\boldsymbol {\varepsilon }}\|={\sqrt {{\boldsymbol {\varepsilon }}:{\boldsymbol {\varepsilon }}}}={\sqrt {\varepsilon _{ij}\varepsilon _{ij}}}}$
(18)

or as the energy norm:

 ${\displaystyle \varepsilon _{eq}={\sqrt {\frac {{\boldsymbol {\varepsilon }}:{\boldsymbol {E}}:{\boldsymbol {\varepsilon }}}{E}}}={\sqrt {{\frac {1}{E}}E_{ijkl}\varepsilon _{ij}\varepsilon _{kl}}}}$
(19)

where ${\textstyle E_{ijkl}}$ are the components of the elastic constitutive tensor ${\textstyle {\boldsymbol {E}}}$ and normalization by Young's modulus ${\textstyle E}$ is introduced to obtain a strain-like quantity.

The two norms of ${\textstyle {\boldsymbol {\varepsilon }}}$ introduced above, lead to symmetric elastic domain in tension and compression. Nevertheless, several materials (rocks, concrete, ceramics, etc.) often show a non symmetric damage surface, i.e., the yield value in compression can be several times the value in tension. In order to overcome this limitation, it is necessary to adjust the definition of the equivalent strain.

Micro-cracks in concrete grow mainly if the material is stretched, and so it is natural to take into account only normal strains that are positive (tensile strains) and neglect those that are negative (compressive strains). This leads to the so-called Mazars definition of equivalent strains [31]:

 ${\displaystyle \varepsilon _{eq}=\|\langle {\boldsymbol {\varepsilon }}\rangle \|={\sqrt {\langle {\boldsymbol {\varepsilon }}\rangle :\langle {\boldsymbol {\varepsilon }}\rangle }}}$
(20)

or to its energetic counterpart:

 ${\displaystyle \varepsilon _{eq}={\sqrt {\frac {\langle {\boldsymbol {\varepsilon }}\rangle :{\boldsymbol {E}}:\langle {\boldsymbol {\varepsilon }}\rangle }{E}}}}$
(21)

where McAuley brackets ${\textstyle \langle .\rangle }$ denote the "positive part" operator. For scalars ${\textstyle \langle x\rangle =\max(0,x)}$, i.e.,

 ${\displaystyle \langle x\rangle =\left\{{\begin{array}{ll}x&{\hbox{if }}x>0{\hbox{ }}\\0&{\hbox{if }}x\leq {\hbox{ 0 }}\end{array}}\right.}$
(22)

For symmetric tensors, such as the strain tensor ${\textstyle {\boldsymbol {\varepsilon }}}$, the positive part is a tensor with the same principal directions ${\textstyle {\boldsymbol {n_{i}}}}$ as the original one, but replacing its principal values ${\textstyle \varepsilon _{i}}$ by their positive parts ${\textstyle \langle \varepsilon _{i}\rangle }$. The subscript ${\textstyle i}$ ranges from 1 to 3 (the number of spatial dimensions). If we consider the following diagonalized strain tensor:

 ${\displaystyle {\boldsymbol {\varepsilon }}={\begin{bmatrix}\varepsilon _{1}&0&0\\0&\varepsilon _{2}&0\\0&0&\varepsilon _{3}\end{bmatrix}}}$
(23)

the positive part of ${\textstyle {\boldsymbol {\varepsilon }}}$ is expressed as:

 ${\displaystyle \langle {\boldsymbol {\varepsilon }}\rangle ={\begin{bmatrix}\langle \varepsilon _{1}\rangle &0&0\\0&\langle \varepsilon _{2}\rangle &0\\0&0&\langle \varepsilon _{3}\rangle \end{bmatrix}}}$
(24)

Thereby, knowing that ${\textstyle \langle {\boldsymbol {\varepsilon }}\rangle :\langle {\boldsymbol {\varepsilon }}\rangle =tr(\langle {\boldsymbol {\varepsilon }}\rangle ^{T}\langle {\boldsymbol {\varepsilon }}\rangle )}$, definition (20) can be rewritten as:

 ${\displaystyle \varepsilon _{eq}={\sqrt {\sum _{i=1}^{3}\langle \varepsilon _{i}\rangle ^{2}}}}$
(25)

An alternative formula, called the modified von Mises definition [32], reads:

 ${\displaystyle \varepsilon _{eq}={\frac {\kappa {-1}}{2\kappa (1-2\nu )}}I_{1}+{\frac {1}{2\kappa }}{\sqrt {\left({\frac {\kappa {-1}}{1-2\nu }}I_{1}\right)^{2}+{\frac {12\kappa }{(1+\nu )^{2}}}J_{2}}}}$
(26)

where ${\textstyle \kappa }$ is a model parameter that sets the ratio between the uniaxial compressive strength and the uniaxial tensile strength, ${\textstyle \nu }$ is the Poisson's ratio of the material, ${\textstyle I_{1}}$ is the first invariant of the strain tensor, and ${\textstyle J_{2}}$ is the second invariant of the deviatoric strain tensor. Given a generic symmetric strain tensor ${\textstyle {\boldsymbol {\varepsilon }}}$:

 ${\displaystyle {\boldsymbol {\varepsilon }}={\begin{bmatrix}\varepsilon _{xx}&\varepsilon _{xy}&\varepsilon _{xz}\\\varepsilon _{xy}&\varepsilon _{yy}&\varepsilon _{yz}\\\varepsilon _{xz}&\varepsilon _{yz}&\varepsilon _{zz}\end{bmatrix}}}$
(27)

The first invariant ${\textstyle I_{1}}$ is the trace of the strain tensor:

 ${\displaystyle I_{1}=tr({\boldsymbol {\varepsilon }})=\varepsilon _{xx}+\varepsilon _{yy}+\varepsilon _{zz}}$
(28)

One can always decompose the strain tensor into its volumetric and deviatoric parts ${\textstyle {\boldsymbol {\varepsilon }}={\boldsymbol {\varepsilon }}_{v}+{\boldsymbol {\varepsilon }}_{d}}$:

 ${\displaystyle {\boldsymbol {\varepsilon }}_{v}={\begin{bmatrix}{\frac {I_{1}}{3}}&0&0\\0&{\frac {I_{1}}{3}}&0\\0&0&{\frac {I_{1}}{3}}\end{bmatrix}}}$
(29)
 ${\displaystyle {\boldsymbol {\varepsilon }}_{d}={\begin{bmatrix}\varepsilon _{xx}-{\frac {I_{1}}{3}}&\varepsilon _{xy}&\varepsilon _{xz}\\\varepsilon _{xy}&\varepsilon _{yy}-{\frac {I_{1}}{3}}&\varepsilon _{yz}\\\varepsilon _{xz}&\varepsilon _{yz}&\varepsilon _{zz}-{\frac {I_{1}}{3}}\end{bmatrix}}}$
(30)

From the deviatoric strain tensor ${\textstyle {\boldsymbol {\varepsilon }}_{d}}$ we can calculate ${\textstyle J_{1}}$ and ${\textstyle J_{2}}$:

 ${\displaystyle J_{1}=tr({\boldsymbol {\varepsilon }}_{d})=\varepsilon _{xx}-{\frac {I_{1}}{3}}+\varepsilon _{yy}-{\frac {I_{1}}{3}}+\varepsilon _{zz}-{\frac {I_{1}}{3}}=I_{1}-3{\frac {I_{1}}{3}}=0}$
(31)
 ${\displaystyle J_{2}}$ ${\displaystyle ={\frac {1}{2}}({\boldsymbol {\varepsilon }}_{d}:{\boldsymbol {\varepsilon }}_{d}-J_{1}^{2})={\frac {1}{2}}({\boldsymbol {\varepsilon }}_{d}:{\boldsymbol {\varepsilon }}_{d})}$ ${\displaystyle ={\frac {1}{3}}\left[\varepsilon _{xx}^{2}+\varepsilon _{yy}^{2}+\varepsilon _{zz}^{2}-\varepsilon _{xx}\varepsilon _{yy}-\varepsilon _{xx}\varepsilon _{zz}-\varepsilon _{yy}\varepsilon _{zz}\right]+}$ ${\displaystyle \varepsilon _{xy}^{2}+\varepsilon _{xz}^{2}+\varepsilon _{yz}^{2}}$
(32)

Another possibility is to work with an evolution of the model proposed by Simo and Ju [23], using the energy norm of the strain tensor and modifying it to take into account the different degradation in tension and compression of concrete-like materials. In this case, the equivalent strain takes the form:

 ${\displaystyle \varepsilon _{eq}=\left(\theta +{\frac {1-\theta }{\kappa }}\right){\sqrt {{\boldsymbol {\varepsilon }}:{\boldsymbol {E}}:{\boldsymbol {\varepsilon }}}}}$
(33)

where the parameter ${\textstyle \theta }$ is a weighting factor depending on the effective stresses ${\textstyle {\bar {\boldsymbol {\sigma }}}}$, given by:

 ${\displaystyle \theta ={\frac {\sum _{i=1}^{3}\langle {\bar {\sigma }}_{i}\rangle }{\sum _{i=1}^{3}|{\bar {\sigma }}_{i}|}}}$
(34)

As mentioned before, the parameter ${\textstyle \kappa }$ is defined by means of the ratio between the compression elastic limit ${\textstyle \sigma _{y}^{c}}$ and the tension elastic limit ${\textstyle \sigma _{y}^{t}}$, i.e.:

 ${\displaystyle \kappa ={\frac {\sigma _{y}^{c}}{\sigma _{y}^{t}}}}$
(35)

In the case of concrete ${\textstyle \kappa \approx {10}}$.

The equivalent strains presented in this subsection lead to different damage surfaces but all three share a common trait, the elastic limit in tension ${\textstyle \sigma _{y}^{t}}$ is much lower than the elastic limit in compression ${\textstyle \sigma _{y}^{c}}$ (Figure 5.b, 5.c, 5.d). In contrast, the equivalent strain obtained from the original energy norm definition (19) leads to a symmetric damage surface in which there is no difference between the elastic limit in tension and compression (Figure 5.a).

 Figure 5: Damage surfaces for different equivalent strains in the 2D principal stress space.

### 2.2.3 Damage evolution laws

As we stated for the uniaxial case, one of the basic components of a damage model is the law governing the evolution of the damage variable (7). There are various damage governing laws ${\textstyle g(r)}$ that can be effectively used to model damage growth in geo-materials, and in the following lines we will present some of the most widely used options.

Two typical choices to describe the evolution of damage above the damage threshold ${\textstyle r_{0}}$ are the exponential law [20]:

 ${\displaystyle g(r)=1-{\frac {r_{0}(1-A)}{r}}-A\exp\{-B(r-r_{0})\}}$
(36)

and the polynomial law [33]:

 ${\displaystyle g(r)=1-{\frac {1}{1+B(r-r_{0})+A(r-r_{0})^{2}}}}$
(37)

In (36) and (37) parameter ${\textstyle A}$ is associated to residual strength and parameter ${\textstyle B}$ controls the slope of the softening branch at the peak of the stress-strain curve.

Another option, especially suited for simplified analyses, is the linear softening law. Limiting the state variable ${\textstyle r}$ between the damage threshold ${\textstyle r_{0}}$ and a maximum admissible value ${\textstyle r_{max}}$, damage evolves according to:

 ${\displaystyle g(r)={\frac {r_{max}}{r_{max}-r_{0}}}\left(1-{\frac {r_{0}}{r}}\right)}$
(38)

which leads to a linear softening branch in the stress-strain curve.

The effect of using one form of damage evolution law or another is clearly seen in a uniaxial case. Figure 6 shows the aspect of the stress-strain curve in a unidimensional case using the three softening laws presented above. For all three cases we have used a damage threshold of ${\textstyle r_{0}=2\cdot {10}^{-4}}$ and a Young's modulus of ${\textstyle E=28000}$ MPa. The exponential law is represented with ${\textstyle A=0.8}$ and ${\textstyle B=15000}$, the polynomial law with ${\textstyle A=0.8}$ and ${\textstyle B=30000}$, and for the linear law we have defined ${\textstyle r_{max}=5\cdot {10}^{-4}}$. One can see that when the strain exceeds the damage threshold after the elastic branch ${\textstyle \varepsilon {>}r_{0}}$, a softening post-peak branch starts as a result of the damage growing.

 Figure 6: Unidimensional stress-strain curves for different damage evolution laws.

An alternative exponential softening model was proposed in [27]:

 ${\displaystyle g(r)=1-{\frac {r_{0}}{r}}\exp \left\{A_{f}\left(1-{\frac {r}{r_{0}}}\right)\right\}}$
(39)

The parameter ${\textstyle A_{f}}$ is obtained from the following expression:

 ${\displaystyle A_{f}=\left({\frac {G_{f}}{l_{f}(r_{0})^{2}}}-{\frac {1}{2}}\right)^{-1}\geq 0}$
(40)

where ${\textstyle G_{f}}$ is the specific fracture energy per unit area, ${\textstyle l_{f}}$ is the characteristic length for the fractured domain, usually taken as the characteristic length of the finite elements, and ${\textstyle r_{0}}$ is the damage threshold which, for the Simo and Ju model, can be computed from:

 ${\displaystyle r_{0}={\frac {\sigma _{y}^{t}}{\sqrt {E}}}}$
(41)

being ${\textstyle \sigma _{y}^{t}}$ the tensile strength of the material, and ${\textstyle E}$ the Young's modulus.

Another popular damage evolution law specifically designed for concrete was proposed by Mazars [31,20]. He introduced two damage variables, ${\textstyle d_{c}}$ and ${\textstyle d_{t}}$, that are computed from the same equivalent strain (25) using two different damage functions, ${\textstyle g_{c}}$ and ${\textstyle g_{t}}$. Function ${\textstyle g_{c}}$ is identified from the uniaxial compressive test while ${\textstyle g_{t}}$ corresponds to the tensile test. The evolution law governing damage ${\textstyle g(r)}$ results from the combination of the two parts:

 ${\displaystyle g(r)=\alpha _{t}^{\beta }g_{t}(r)+(1-\alpha _{t})^{\beta }g_{c}(r)}$
(42)

Functions characterizing the evolution of damage in compression ${\textstyle g_{c}}$ and in tension ${\textstyle g_{t}}$ were proposed in the exponential form:

 ${\displaystyle g_{c}(r)=1-{\frac {r_{0}(1-A_{c})}{r}}-A_{c}\exp\{-B_{c}(r-r_{0})\}}$
(43)
 ${\displaystyle g_{t}(r)=1-{\frac {r_{0}(1-A_{t})}{r}}-A_{t}\exp\{-B_{t}(r-r_{0})\}}$
(44)

where ${\textstyle A_{c}}$, ${\textstyle B_{c}}$, ${\textstyle A_{t}}$, ${\textstyle B_{t}}$ are material parameters related to the shape of uniaxial stress-strain diagrams.

The coefficient ${\textstyle \alpha _{t}}$ in (42) ranges from 0 to 1 and takes into account the character of the stress state. It is evaluated from:

 ${\displaystyle \alpha _{t}=\sum _{i=1}^{3}{\frac {\varepsilon _{ti}\langle \varepsilon _{i}\rangle }{\varepsilon _{eq}^{2}}}}$
(45)

where ${\textstyle \varepsilon _{ti}}$ are the principal strains due to positive effective stresses, i.e., the principal values of ${\textstyle {\boldsymbol {\varepsilon }}_{t}={\boldsymbol {E}}^{-1}:\langle {\boldsymbol {E}}:{\boldsymbol {\varepsilon }}\rangle }$:

 ${\displaystyle {\boldsymbol {\varepsilon }}_{t}=-{\frac {\nu }{E}}tr(\langle {\boldsymbol {\bar {\sigma }}}\rangle ){\boldsymbol {I}}+{\frac {1+\nu }{E}}\langle {\boldsymbol {\bar {\sigma }}}\rangle }$
(46)

being ${\textstyle {\boldsymbol {I}}}$ the identity tensor.

The parameter ${\textstyle \beta }$ in (42) was equal to 1 in the original version of the model [31]. When it is higher than 1, ${\textstyle \beta }$ allows to slow down the evolution of damage under shear loading (when principal stresses have different sign).

The definition of ${\textstyle \varepsilon _{ti}}$ tells us that if all principal stresses are negative then ${\textstyle \alpha _{t}=0}$ and ${\textstyle d=d_{c}=g_{c}(r)}$, and if all principal stresses are positive we have ${\textstyle \alpha _{t}=1}$ and ${\textstyle d=d_{t}=g_{t}(r)}$. These are the "purely compressive" and "purely tensile" stress states. For intermediate stress states the value of ${\textstyle d}$ is between ${\textstyle d_{c}}$ and ${\textstyle d_{t}}$, depending on the relative magnitudes of tensile and compressive stresses.

### 2.2.4 Elastic-damage tangent constitutive tensor

In this section we will present the tangent constitutive tensor, which gives us an important advantage, from a computational point of view, when dealing with incremental iterative solution procedures. Basically, with the tangent constitutive tensor one can achieve a quadratic convergence rate, whereas with the secant constitutive tensor, the solution is limited to a linear convergence.

The damaged or secant constitutive tensor ${\textstyle {\boldsymbol {E}}_{sec}}$ introduced in (12) and (13) links the total stress to total strain and plays the role of the tangent stiffness only for processes with elastic loading or unloading with constant damage (${\textstyle f<0}$). To obtain the tangent stiffness for a loading process with growing damage (${\textstyle f=0}$ and ${\textstyle {\dot {f}}=0}$), it is necessary to find the link between stress and strain increments. This concepts are clearly seen in a unidimensional case (see Figure 7). In this regard, knowing the difference between secant and tangent moduli (see Figure 8.a and 8.b) is useful to understand the generalization to secant and tangent constitutive tensors.

 Figure 8: Secant and tangent moduli.

In essence, we define the elastic-damage tangent constitutive tensor ${\textstyle {\boldsymbol {E}}_{tan}}$ as the one that satisfies the following relation:

 ${\displaystyle {\dot {\boldsymbol {\sigma }}}={\boldsymbol {E}}_{tan}:{\dot {\boldsymbol {\varepsilon }}}}$
(47)

Now, by considering the equation (13), we can obtain the rate of change of stress as follows:

 ${\displaystyle {\dot {\boldsymbol {\sigma }}}({\boldsymbol {\varepsilon }},d)}$ ${\displaystyle ={\frac {\partial {\boldsymbol {\sigma }}}{\partial {\boldsymbol {\varepsilon }}}}:{\dot {\boldsymbol {\varepsilon }}}+{\frac {\partial {\boldsymbol {\sigma }}}{\partial d}}{\dot {d}}=(1-d){\boldsymbol {E}}:{\dot {\boldsymbol {\varepsilon }}}-{\boldsymbol {E}}:{\boldsymbol {\varepsilon }}{\dot {d}}}$ ${\displaystyle =(1-d){\boldsymbol {E}}:{\dot {\boldsymbol {\varepsilon }}}-{\bar {\boldsymbol {\sigma }}}\otimes {\dot {d}}}$
(48)

from which we should distinguish two possible situations:

1. A process with elastic loading or unloading (${\textstyle {\dot {r}}=0\Rightarrow {\dot {d}}=0}$)
2. With ${\textstyle {\dot {d}}=0}$ equation (48) becomes ${\textstyle {\dot {\boldsymbol {\sigma }}}=(1-d){\boldsymbol {E}}:{\dot {\boldsymbol {\varepsilon }}}}$, from which we obtain that the elastic-damage tangent constitutive tensor coincides with the secant constitutive tensor:

 ${\displaystyle {\boldsymbol {E}}_{tan}={\boldsymbol {E}}_{sec}=(1-d){\boldsymbol {E}}}$
(49)
3. A loading process with growing damage (${\textstyle {\dot {r}}>0\Rightarrow {\dot {d}}>0}$)
4. In this case, in order to obtain an explicit expression of the tangent constitutive tensor, the damage rate ${\textstyle {\dot {d}}}$ should be expressed in terms of the strain rate ${\textstyle {\dot {\boldsymbol {\varepsilon }}}}$. This can be achieved by imposing the consistency condition ${\textstyle {\dot {f}}=0}$ in equation (16) and combining it with the rate of equation (7):

 ${\displaystyle {\dot {d}}={\frac {dg}{dr}}{\dot {r}}={\frac {dg}{dr}}{\dot {\varepsilon }}_{eq}={\frac {dg}{dr}}{\frac {\partial \varepsilon _{eq}}{\partial {\boldsymbol {\varepsilon }}}}:{\dot {\boldsymbol {\varepsilon }}}}$
(50)

For convenience, we introduce symbols ${\textstyle g'}$ for the derivative of the damage function ${\textstyle dg/dr}$, and ${\textstyle {\boldsymbol {\eta }}}$ for the second order tensor obtained from the partial derivative of the equivalent strain with respect to the strain tensor ${\textstyle \partial \varepsilon _{eq}/\partial {\boldsymbol {\varepsilon }}}$. Thereby, substituting ${\textstyle {\dot {d}}=g'{\boldsymbol {\eta }}:{\dot {\boldsymbol {\varepsilon }}}}$ into the rate of change of stress (48), we get:

 ${\displaystyle {\dot {\boldsymbol {\sigma }}}=(1-d){\boldsymbol {E}}:{\dot {\boldsymbol {\varepsilon }}}-{\bar {\boldsymbol {\sigma }}}\otimes (g'{\boldsymbol {\eta }}:{\dot {\boldsymbol {\varepsilon }}})=[(1-d){\boldsymbol {E}}-g'{\bar {\boldsymbol {\sigma }}}\otimes {\boldsymbol {\eta }}]:{\dot {\boldsymbol {\varepsilon }}}}$
(51)

and, consequently, the elastic-damage constitutive tensor results:

 ${\displaystyle {\boldsymbol {E}}_{tan}=(1-d){\boldsymbol {E}}-g'{\bar {\boldsymbol {\sigma }}}\otimes {\boldsymbol {\eta }}}$
(52)

where the form of ${\textstyle g'}$ and ${\textstyle {\boldsymbol {\eta }}}$ depend on the damage evolution law, and on the equivalent strain expression chosen for the damage model.

For instance, for a model with an equivalent strain ${\textstyle \varepsilon _{eq,1}}$ based on the energy norm of equation (19), tensor ${\textstyle {\boldsymbol {\eta }}_{1}}$ is given by:

 ${\displaystyle {\boldsymbol {\eta }}_{1}={\frac {\partial \varepsilon _{eq,1}}{\partial {\boldsymbol {\varepsilon }}}}={\frac {1}{2{\sqrt {\frac {{\boldsymbol {\varepsilon }}:{\boldsymbol {E}}:{\boldsymbol {\varepsilon }}}{E}}}}}{\frac {1}{E}}2{\boldsymbol {E}}:{\boldsymbol {\varepsilon }}={\frac {\bar {\boldsymbol {\sigma }}}{E\varepsilon _{eq,1}}}}$
(53)

and the resulting tangent constitutive tensor for this particular case is:

 ${\displaystyle {\boldsymbol {E}}_{tan,1}=(1-d){\boldsymbol {E}}-{\frac {g'}{E\varepsilon _{eq,1}}}{\bar {\boldsymbol {\sigma }}}\otimes {\bar {\boldsymbol {\sigma }}}}$
(54)

which preserves symmetry (${\textstyle {\boldsymbol {E}}_{tan,1}^{T}={\boldsymbol {E}}_{tan,1}}$). It must be noticed though that, for other definitions of equivalent strain, this kind of symmetry is lost.

The typical scheme for the numerical implementation of the elastic-damage tangent constitutive tensor is shown in Table (2).

 Input data for time ${\displaystyle t+1}$: ${\displaystyle {\boldsymbol {E}}}$, ${\displaystyle {\bar {\boldsymbol {\sigma }}}^{t+1}}$, ${\displaystyle r^{t}}$, ${\displaystyle r^{t+1}}$, ${\displaystyle d^{t+1}}$ If (${\textstyle r^{t+1}>r^{t}}$) Loading with growing damage (${\textstyle {\dot {r}}>0\Rightarrow {\dot {d}}>0}$): 1. Compute damage function derivative: ${\displaystyle g'=dg/dr}$ 2. Calculate equivalent strain partial derivative: ${\displaystyle {\boldsymbol {\eta }}=\partial \varepsilon _{eq}/\partial {\boldsymbol {\varepsilon }}}$ 3. Compute tangent constitutive tensor with (52): ${\displaystyle {\boldsymbol {E}}_{tan}^{t+1}=(1-d^{t+1}){\boldsymbol {E}}-g'{\bar {\boldsymbol {\sigma }}}^{t+1}\otimes {\boldsymbol {\eta }}}$ Else Elastic loading or unloading (${\textstyle {\dot {r}}=0\Rightarrow {\dot {d}}=0}$): ${\displaystyle {\boldsymbol {E}}_{tan}^{t+1}={\boldsymbol {E}}_{sec}^{t+1}=(1-d^{t+1}){\boldsymbol {E}}}$

## 2.3 Local and Non-local Damage Models

In previous sections we have presented the basic ideas of isotropic damage models with a unique scalar variable ${\textstyle d}$. Although these models are quite simple, they are often used to model the failure of concrete and other quasi-brittle materials that show important strain localization. If the damage parameter depends only on the strain state at the point under consideration, the numerical simulations exhibit a pathological mesh dependence, and physically unrealistic results are obtained.

This is the typical behaviour of the so-called local damage models, which are not able to properly describe both the thickness of localization and distance between them. They suffer from mesh sensitivity (for size and alignment) and produce unreliable results. The strains concentrate in one element wide zones and the computed force-displacement curves are mesh-dependent. The reason is that differential equations of motion change their type (from elliptic to hyperbolic in static problems) and the rate boundary value problem becomes ill-posed.

Thus, classical constitutive models require an extension in the form of a characteristic length to properly model the thickness of localized zones. Such extension can be done with micro-polar [34,35], strain gradient [36,37,38], viscous [39,40,41] and non-local terms [42,43,44,45]. In this work we have developed the latter approach.

Similarly to previous sections, in order to make it more understandable, we will present the inconveniences arisen due to strain localization with a uniaxial example. After that, the basic concepts of the implemented non-local damage model will be stated.

### 2.3.1 Strain localization phenomenon

The idea of modelling damaged concrete and other quasi-brittle materials as strain-softening continua, was not immediately accepted by all the scientific community. Actually, most of the early analyses were not truly objective and, upon mesh refinement, their results would not converge to a robust solution. Let us explain the nature of the problem by means of a unidimensional example.

Consider a straight bar with a constant cross section ${\textstyle S}$ and a total length ${\textstyle L_{0}}$ under uniaxial tension (Figure 9). The material is assumed to obey a simple stress-strain law with linear elasticity up to the peak stress ${\textstyle \sigma _{y}}$, followed by linear softening (Figure 10). If the bar is loaded in tension by an applied displacement ${\textstyle u}$ at its right extreme, the response remains linear elastic up to ${\textstyle u_{y}=L_{0}\varepsilon _{y}}$, instant at which the force transmitted by the bar (reaction at the support) attains its maximum value ${\textstyle F_{y}=\sigma _{y}S}$. After that, the resistance of the bar starts decreasing until the strain reaches ${\textstyle \varepsilon _{f}}$ and the transmitted stress completely disappears.

 Figure 9: Bar under uniaxial tension.
 Figure 10: Stress-strain diagram with linear softening.

Equilibrium equations imply that axial forces are constant along the bar and so the stress profile must remain also uniform (Figure 11). However, at a given stress level ${\textstyle \sigma _{1}}$ between zero and ${\textstyle \sigma _{y}}$, there are actually two values of strain, ${\textstyle \varepsilon _{1,u}}$ and ${\textstyle \varepsilon _{1,s}}$, that satisfy the constitutive equation (Figure 12). This is quite straightforward if one considers that each cross section can be either damaged, or undamaged. Indeed, an undamaged section will be on the elastic unloading branch with ${\textstyle \sigma _{1}=E\varepsilon _{1,u}}$, whereas a damaged one will fall in the softening branch with ${\textstyle \sigma _{1}=(1-d_{1})E\varepsilon _{1,s}}$.

 Figure 11: Axial force acting along the bar.
 Figure 12: Possible strain values corresponding to the same stress level.

Thereby, the strain profile along the bar does not have to be necessarily uniform. In fact, any piecewise constant strain distribution that jumps between the two possible strain values represents a valid solution (Figure 13). In Figure 13 we have denoted by ${\textstyle L_{s}}$ the cumulative length of the softening regions and by ${\textstyle L_{u}=L_{0}-L_{s}}$ the cumulative length of the unloading regions. When stress is completely relaxed, the strain in the unloading region is ${\textstyle \varepsilon _{u}=0}$ and the strain in the softening region is ${\textstyle \varepsilon _{s}=\varepsilon _{f}}$. Thus the total elongation of the bar in this case is ${\textstyle u_{f}=L_{u}\varepsilon _{u}+L_{s}\varepsilon _{s}=L_{s}\varepsilon _{f}}$. At this point, although ${\textstyle \varepsilon _{f}}$ is perfectly known from the linear softening law, ${\textstyle L_{s}}$ is totally undetermined. This means that the problem has infinite possible solutions with its corresponding post-peak branches of the load-displacement diagram (Figure 14). This fan of post-peak branches is bounded on one side by the solution with uniform softening (${\textstyle u_{f}=L_{0}\varepsilon _{f}}$) and on the other side by the solution with uniform unloading (${\textstyle u_{f}=0}$). The first limit corresponds to a totally damaged bar while the latter represents the case in which the bar is unloaded before any damage takes place. All the other solutions describe processes in which a part of the bar is damaged. However, it is not that easy to know which of all these solutions is the one that reflects better the actual failure process.

 Figure 13: Piecewise constant strain profile along the bar.
 Figure 14: Fan of possible post-peak branches of the load-displacement diagram.

The ambiguity is removed if imperfections are taken into account. The material properties and sectional dimensions of a real bar are not perfectly uniform. Thereby, supposing that there is a small region where the strength is lower than in the remaining portion of the bar, when the applied stress reaches the peak of that weaker region, softening starts and the stress decreases. Consequently, the material outside the damaged region must unload elastically because its strength has not been exhausted. This leads to the conclusion that the size of the softening region is related to the size of the region with minimum strength. The problem is that such a region can be arbitrarily small and so the corresponding softening branch is arbitrarily close to the elastic unload branch. Therefore, the standard strain-softening continuum formulation leads to a solution with several pathological features:

• The softening region is infinitely small.
• The load-displacement diagram always shows snap-back, regardless of the size of the structure and the ductility of the material.
• The total amount of energy dissipated during the failure process tends to zero.

From the mathematical point of view, these problematic features are related to the loss of ellipticity of the governing differential equation. The boundary value problem becomes ill-posed as a result, i.e., it no longer has a unique solution with continuous dependence on the given data.

From the numerical point of view, these inconveniences are manifested by a pathological sensitivity of the results to the size of the finite elements. For instance, let us suppose that the bar is discretized uniformly by ${\textstyle n}$ two-node elements with linear displacement interpolation and that the weakest region is located at the middle of the bar. The numerical algorithm will capture a very localized solution with a softening region extending over one only element. In consequence, the softening length will decrease as the number of elements increases (${\textstyle L_{s}=L_{0}/n}$) and thus the softening post-peak branch will depend completely on the number of elements of the mesh. Indeed, as it is shown in Figure 15, for ${\textstyle n=1}$ all the bar is damaged and the softening length is the total length of the bar ${\textstyle L_{s}=L_{0}}$, whereas for ${\textstyle n>1}$ the softening region is more localized with strains becoming especially important at the damaged element. In the limit case of ${\textstyle n\rightarrow \infty }$ the softening branch would coincide with the elastic branch.

 Figure 15: Mesh dependence of the numerical results. (a) Load-displacement diagrams for different number of elements. (b) Strain profiles for a prescribed imposed displacement.

In this section we have seen a problem that commonly arises in the simulation of damage processes involving quasi-brittle materials: the strain localization phenomenon. Although only a uniaxial case has been discussed, this problem is also present in multidimensional problems with similar effects on the numerical results.

The simple one-dimensional example presented in this section illustrates the essence of the problem with localization of inelastic strain into a zone of an arbitrarily small width. In uniaxial cases, localization occurs when the peak of the stress-strain diagram is reached, independently of the specific constitutive model used. In multiple dimensions, the analysis of the localization process is more complicated and the derivation of a criteria for the start of localization represents a challenging problem.

Mesh refinement in multiple dimensions leads to a reduction of the total dissipated energy (area under the load-displacement curve) with a lower peak load and a more brittle response. Apart from this, upon further refinement, one can also face convergence difficulties due to the abrupt change of strain distribution, from a smoothly distributed to a highly localized one. These effects will be shown more clearly with the simulation of a bi-dimensional case later on (see Chapter 4).

### 2.3.2 Regularization of the problem

In the present work, two different damage models were implemented in the research of a robust code that avoided pathological sensitivity of the finite element results to the mesh size.

#### 2.3.2.1 Partially regularized local damage model

As a first attempt, we tried a simple partially regularized local damage model, based on the crack band models [46].

This model is, essentially, an isotropic damage model following the classical local damage theory, in which an energy based adjustment of the stress-strain diagram, depending on the size of the element, is introduced in the definition of the damage parameter.

In this regard, the model makes use of the damage evolution law in (39), which depends on the characteristic length of the fractured domain ${\textstyle l_{f}}$ included in (40). This characteristic length of the fractured domain is actually what allows to partially regularize the model.

In essence, for elements larger than a prescribed limit length ${\textstyle l_{lim}}$, the fracture length takes the value of the characteristic length of the element ${\textstyle l_{e}}$, whereas for elements smaller than the limit length, the fracture length is taken as such limit length.

 ${\displaystyle l_{f}=\left\{{\begin{array}{ll}l_{e}&{\hbox{if }}l_{e}>l_{lim}{\hbox{ }}\\l_{lim}&{\hbox{if }}l_{e}\leq l_{lim}{\hbox{ }}\end{array}}\right.}$
(55)

Therefore, before calculating the damage parameter, one will have to compute the characteristic length of the element and then compare it with the limit length, a material parameter usually related to the maximum aggregate size of the composite material. The characteristic length of the element, on the other hand, can be computed from the dimensions of the element. In two-dimensional analysis, for instance, the characteristic length of the element can be defined as the diameter of the circle that contains the area of the element.

 ${\displaystyle l_{e}={\sqrt {\frac {4A_{e}}{\pi }}}}$
(56)

where ${\textstyle A_{e}}$ is the area of the considered element.

This approach is endowed with some, but not all, of the properties of fully regularized damage models. It can ensure a correct energy dissipation in a localized damage band, but the width of the numerically resolved fracture process zone depends on the element size and tends to zero as the mesh is refined. This is why this approach cannot be considered as a true localization limiter. It provides only a partial regularization of the problem in the sense that global response characteristics do not exhibit spurious mesh sensitivity, but the mesh-induced directional bias is still present.

Scaling of the stress-strain diagram is straightforward only for models that explicitly control the evolution of inelastic strain, e.g., for softening plasticity [47] or smeared crack models [27]. In those cases, the desired scaling effect is achieved by a modification of the hardening modulus (derivative of stress with respect to inelastic strain). Nevertheless, in continuum damage mechanics, non-linearity and softening are controlled by the damage evolution law, and the reduction factor ${\textstyle 1-d}$ multiplies the total strain. For this reason, it is not easy to scale only the post-localization part of strain while keeping the unloading part unaffected.

Furthermore, in some cases, diffuse softening damage patterns in certain parts of the structure can coexist with localized cracks in other parts, and they may even change during the loading process. In such cases it is virtually impossible to define a reasonable rule for the adjustment of the stress-strain diagram according to the element size.

#### 2.3.2.2 Fully regularized non-local damage model

The introduction of a characteristic length into the constitutive model, and the formulation of a non-local strain-softening model, have been shown to prevent the spurious localization of strain-softening damage, to regularize the boundary value problem, and to ensure numerical convergence to physically meaningful solutions.

In this regard, fully regularized description of localized inelastic deformation is achieved by a proper generalization of the underlying continuum theory, which can be done through two different approaches: generalization of the kinematic relations, i.e., continua with micro-structure (Cosserat-type continua or strain gradient theories), and continua with non-local strain (non-local elasticity); and generalization of constitutive equations, i.e., material models with gradients of internal variables, and materials models with weighted spatial averages of internal variables.

In this work we have worked with the second kind of generalization, because kinematic and equilibrium equations remain standard, and the notions of stress and strain keep their usual meaning. Moreover, in the generalization of constitutive equations through non-local models, we have focused on the integral-type methods.

Integral-type non-local models abandon the classical assumption of locality and admit that stress a certain point depends, not only on the state variables at that point, but also on the distribution of the state variables over the whole body, or over a finite neighbourhood of the point under consideration.

The first models of this type, proposed in the 1960s, aimed at improving the description of elastic wave dispersions in crystals. Non-local elasticity was further developed by Eringen, who later extended it to non-local elasto-plasticity [48,49]. Subsequently, it was found that certain non-local formulations could act as efficient localization limiters with a regularizing effect on problems with strain localization [43].

In a general manner, the non-local integral approach consists in replacing a certain variable by its non-local counterpart, obtained by weighted averaging over a spatial neighbourhood of each point under consideration.

Let ${\textstyle f({\boldsymbol {x}})}$ be some local field in a domain ${\textstyle V}$, the corresponding non-local field is defined as:

 ${\displaystyle {\tilde {f}}({\boldsymbol {x}})=\int _{V}\alpha ({\boldsymbol {x}},{\boldsymbol {\xi }})f({\boldsymbol {\xi }})d{\boldsymbol {\xi }}}$
(57)

where ${\textstyle \alpha ({\boldsymbol {x}},{\boldsymbol {\xi }})}$ is the chosen non-local weighting function.

In an infinite, isotropic and homogeneous medium, the weighting function depends only on the distance ${\textstyle D=\|{\boldsymbol {x}}-{\boldsymbol {\xi }}\|}$ between the source point ${\textstyle {\boldsymbol {\xi }}}$, and the receiver point ${\textstyle {\boldsymbol {x}}}$. Thereby, we usually write ${\textstyle \alpha ({\boldsymbol {x}},{\boldsymbol {\xi }})=\alpha _{0}(\|{\boldsymbol {x}}-{\boldsymbol {\xi }}\|)}$, where ${\textstyle \alpha _{0}(D)}$ is usually chosen as a non-negative function monotonically decreasing for ${\textstyle D\geq 0}$.

One possible ${\textstyle \alpha _{0}(D)}$ is the Gauss distribution function:

 ${\displaystyle \alpha _{0}(D)=\exp \left[-\left({\frac {2D}{l_{c}}}\right)^{2}\right]}$
(58)

where the characteristic length ${\textstyle l_{c}}$ is a material parameter reflecting the internal length of the non-local continuum.

If a bounded support is preferred, one can also truncate the previous function as follows:

 ${\displaystyle \alpha _{0}(D)=\left\{{\begin{array}{ll}\exp \left[-\left({\frac {2D}{l_{c}}}\right)^{2}\right]&{\hbox{if }}|D|\leq R{\hbox{ }}\\0&{\hbox{if }}|D|>R{\hbox{ }}\end{array}}\right.}$
(59)

where the interaction radius ${\textstyle R}$ is a parameter related to the characteristic length ${\textstyle l_{c}}$. In the present work, we have considered ${\textstyle R=l_{c}}$.

Another typical choice for the weighting function is the following truncated quartic polynomial function:

 ${\displaystyle \alpha _{0}(D)=\left\{{\begin{array}{ll}\left(1-{\frac {D^{2}}{R^{2}}}\right)^{2}&{\hbox{if }}|D|\leq R{\hbox{ }}\\0&{\hbox{if }}|D|>R{\hbox{ }}\end{array}}\right.}$
(60)

In essence, the interaction radius ${\textstyle R}$ corresponds to the smallest distance between points ${\textstyle {\boldsymbol {x}}}$ and ${\textstyle {\boldsymbol {\xi }}}$ at which the interaction weight ${\textstyle \alpha _{0}(\|{\boldsymbol {x}}-{\boldsymbol {\xi }}\|)}$ vanishes (for weighting functions with a bounded support) or becomes negligible (for weighting functions with an unbounded support). It represents a very important parameter because it controls the size of the softening region.

The interval, circle, or sphere of radius ${\textstyle R}$, centered at ${\textstyle {\boldsymbol {x}}}$, is called the domain of influence of point ${\textstyle {\boldsymbol {x}}}$. In the vicinity of the boundary of a finite body, it is simply assumed that the averaging is performed only on the part of the domain of influence that lies within the solid (Figure 16).

 Figure 16: Averaging zone when the domain of influence protrudes through the boundary of a solid.

Thereby, if a weighting function with bounded support is chosen, the non-local average at (57) will be calculated as a weighted sum over the values at all the finite element integration points ${\textstyle {\boldsymbol {\xi }}}$ lying within the non-local interaction radius ${\textstyle R}$.

In the application to softening materials, it is often required that the non-local operator do not alter a uniform field (consistency of order 0), which means that the weighting function must satisfy the normalizing condition:

 ${\displaystyle {\begin{array}{cc}\displaystyle \int _{V}\alpha ({\boldsymbol {x}},{\boldsymbol {\xi }})d{\boldsymbol {\xi }}=1&\forall {\boldsymbol {x}}\in V\end{array}}}$
(61)

In order to satisfy it, the weighting function is usually rescaled as:

 ${\displaystyle \alpha ({\boldsymbol {x}},{\boldsymbol {\xi }})={\frac {\alpha _{0}(\|{\boldsymbol {x}}-{\boldsymbol {\xi }}\|)}{\displaystyle \int _{V}\alpha _{0}(\|{\boldsymbol {x}}-{\boldsymbol {\zeta }}\|)d{\boldsymbol {\zeta }}}}}$
(62)

A suitable non-local damage formulation that restores well-posedness of the boundary value problem is obtained if damage is computed from the non-local equivalent strain [50].

In the loading function (16), the local equivalent strain ${\textstyle \varepsilon _{eq}}$ is replaced by its weighted spatial average:

 ${\displaystyle {\tilde {\varepsilon }}_{eq}({\boldsymbol {x}})=\int _{V}\alpha ({\boldsymbol {x}},{\boldsymbol {\xi }})\varepsilon _{eq}({\boldsymbol {\xi }})d{\boldsymbol {\xi }}}$
(63)

The internal state variable ${\textstyle r}$ is then the largest previously reached value of the non-local equivalent strain:

 ${\displaystyle r(t)=\max \left\{r_{0},\max _{\tau \leq t}{\tilde {\varepsilon }}_{eq}(\tau )\right\}}$
(64)

It is important to note that the damage variable is evaluated from the non-local equivalent strain ${\textstyle {\tilde {\varepsilon }}_{eq}}$, whereas the strains ${\textstyle {\boldsymbol {\varepsilon }}}$ used in (13) to compute the stresses are considered as local. This way, during the elastic range, when the damage variable remains equal to zero, the stress-strain relation is completely local. The process for the evaluation of stresses with this non-local damage model is schematically shown in Table 3.

 Input data for time ${\displaystyle t+1}$: ${\displaystyle {\boldsymbol {E}}}$, ${\displaystyle {\boldsymbol {\varepsilon }}^{t}}$, ${\displaystyle \Delta {\boldsymbol {\varepsilon }}}$, ${\displaystyle r^{t}}$ 1. Determine new strain vector: ${\textstyle {\boldsymbol {\varepsilon }}^{t+1}={\boldsymbol {\varepsilon }}^{t}+\Delta {\boldsymbol {\varepsilon }}}$ 2. Evaluate effective stresses: ${\textstyle {\boldsymbol {\bar {\sigma }}}^{t+1}={\boldsymbol {E}}:{\boldsymbol {\varepsilon }}^{t+1}}$ 3. Compute equivalent strain from the new strain vector: ${\textstyle \varepsilon _{eq}^{t+1}}$ 4. Calculate non-local equivalent strain from (63): ${\textstyle {\tilde {\varepsilon }}_{eq}^{t+1}}$ 5. Update ${\textstyle r}$ with (64): If ${\textstyle ({\tilde {\varepsilon }}_{eq}^{t+1}>r^{t})\Rightarrow r^{t+1}={\tilde {\varepsilon }}_{eq}^{t+1}}$ 6. Update damage variable: ${\textstyle d^{t+1}=g(r^{t+1})}$ 7. Compute stresses: ${\textstyle {\boldsymbol {\sigma }}^{t+1}=(1-d^{t+1}){\boldsymbol {\bar {\sigma }}}^{t+1}}$

# 3 Numerical Implementation

## 3.1 Introduction

So far, we have presented the theoretical basis of the classical continuum damage mechanics, including an important limitation regarding mesh dependency of the solution, and we have seen one regularization technique to avoid the cited inconvenience.

At this point, before going into example tests of continuum damage models, we must actually understand how the explained theory can be applied in a practical way that allows the performance of such experiments. Essentially, this chapter will present a set of methods and tools, directly or indirectly related to the implementation of a damage model, that are necessary to provide a proper framework for the simulation of damage mechanics problems.

In this regard, I would like to start this section by shortly presenting KRATOS, the framework in which the damage model has been implemented.

 Figure 17: KRATOS' logo.

KRATOS is an Open-Source framework for the implementation of numerical methods for the solution of engineering problems. All of its code is written in C++, and Python scripting language is used to define the main procedure of the problems, which significantly improves the flexibility of the framework in time of use. It features a "core and applications" approach, where standard tools (databases, linear algebra, search structures, etc.) come as a part of the core and are available as building blocks in the development of new applications which focus on the solution of particular problems of interest. KRATOS is designed to allow collaborative development by large teams of researchers focusing on modularity as well as on performance. Its ultimate goal is to simplify the development of new multi-disciplinary finite element programs.

Thereby, the implemented damage model must be understood as a "Damage Mechanics Application" inside that global framework, with all the consequences it implies. Perhaps the most important one is that the starting point of the "Damage Mechanics Application" was the already implemented linear-elastic Finite Element code from the "Solid Mechanics Application" of KRATOS. This means that, from the beginning of this work, I have been able to primarily focus on the implementation of the new damage model, knowing that the FEM core was properly working. The only drawback of this is that one needs to continuously adapt the new implementations to the previously set structure. That aside, it is actually advantageous to work inside KRATOS, being able to use the existing functions if convenient.

Knowing the framework of the monograph, in the following sections we are going to briefly review some basic notions on the classical Finite Element Method, focusing on the bi-dimensional elasticity theory (the one behind the examples performed). Then, we will discuss the non-linearity associated to the damage problem and its influence in the strategy chosen to solve it. Furthermore, we will review the basic aspects in the implementation of the damage model, pointing out the most relevant differences between a local and a non-local approach. Finally, we will introduce an adaptive mesh refinement technique that, although still under development, has already shown some interesting features that are worth mentioning.

## 3.2 Overview of the Finite Element Method

Most engineering structures are continuum, and so their behaviour cannot be properly represented in terms of a few number of discrete variables. A rigorous analysis of such structures requires the integration of the differential equations that govern the equilibrium of a generic differential element belonging to them.

The Finite Element Method (FEM) is a numerical technique for finding approximate solutions to boundary value problems for partial differential equations. It uses subdivisions of a whole problem domain, and variational methods to solve the problem by minimizing an associated error function. In the end, FEM connects many simple element equations over many small domains, to approximate a more complex equation over a larger domain.

The integration of partial differential equations helps the engineer to study uni, bi and three-dimensional structural problems, including its time evolution. In deed, although continuum structures are always three-dimensional, if the proper simplification hypothesis are fitted, one can accurately describe the behaviour of a structure by means of uni or bi-dimensional mathematical models, like in earth dams, tunnels, plates, etc.

At this point, in order to properly understand the concepts explained later on, it is important to state the generic stages for the analysis of a structure through the Finite Element Method. They can be summarized as follows:

• First of all, knowing the physics of a structure, the boundary conditions and applied loads, one must select a mathematical model that appropriately describes its behaviour. Besides, one needs to define the mechanical properties of the materials and the regime of deformations of the structure (small or large displacements, static, quasi-static or dynamic analysis, etc.). In this work, we are going to work with the elastic bi-dimensional FEM theory in a small deformation regime and quasi-static analysis. Furthermore, for the formulation of the equilibrium equations, we are going to use the Principle of Virtual Work (PVW).
• Once the mathematical model is chosen, the structure must be discretized in non-intersecting portions called finite elements. These elements are interconnected through nodes, the entities that store the main variables of the problem. In order to obtain the value of these variables at any point of the geometry, one just needs to interpolate the values stored at the nodes. The set of finite elements that conform all the model is the so called mesh.
• From the expression of the PVW one can obtain the stiffness matrix ${\textstyle {\boldsymbol {K}}^{(e)}}$ and the force vector ${\textstyle {\boldsymbol {f}}^{(e)}}$ for each element.
• Afterwards, the global stiffness matrix ${\textstyle {\boldsymbol {K}}}$ and the nodal force vector ${\textstyle {\boldsymbol {f}}}$ are obtained through the assembly of the previous elemental stiffness matrix and force vector.
• The resulting system of equations ${\textstyle {\boldsymbol {K}}{\boldsymbol {a}}={\boldsymbol {f}}}$ must be solved in order to obtain the nodal displacement vector ${\textstyle {\boldsymbol {a}}}$, i.e., the unknown variable of the FEM problem.
• From the obtained nodal displacement vector ${\textstyle {\boldsymbol {a}}}$, it is quite straightforward to calculate the strains ${\textstyle {\boldsymbol {\varepsilon }}}$ and stresses ${\textstyle {\boldsymbol {\sigma }}}$ at each element, as well as the reactions at the nodes with prescribed displacements.
• Finally, one must present and interpret the numerical results, assessing whether the solution is accurately enough captured, and introducing modifications in the previous stages if necessary. It is important to notice that the Finite Element Method is an approximate method, and so results are expected to have some solution errors, discretization errors and/or modeling errors.

The previous stages are schematically represented in Figure 18.

 Figure 18: Global flowchart of the analysis of a system through the Finite Element Method.

Since the basis of the damage simulations performed in this work is a two-dimensional elastic finite element code, from now on we will be focusing on this kind of FEM problems.

### 3.2.1 Two-dimensional elasticity theory

A large variety of structures with a great interest for engineers can be studied with the hypothesis of two-dimensional elasticity. All these structures have in common an approximate right prism shape. Nonetheless, depending on the proportion of such prism dimensions and on the load configuration, the problems involving these structures can be classified in one of the following groups:

• Plane Stress Problems: one of the dimensions of the structure (thickness) is much smaller than the other two, and loads are contained in the mid-surface of that dimension (Figure 19).
• Plane Strain Problems: one of the dimensions of the structure (length) is much larger than the other two, and loads act uniformly distributed along this length, and are contained in planes orthogonal to the axis passing through the center of gravity of all cross sections (Figure 20).
 Figure 19: Example of a plane stress problem.
 Figure 20: Example of a plane strain problem.

#### 3.2.1.1 Displacements, strains and stresses fields

In two-dimensional elasticity theory, a series of hypothesis are used in order to simplify the expressions defining displacements, strains and stresses. In the following lines we will be reviewing them.

The geometrical characteristics and load properties of a structure working in a plane state (plane stress or plane strain) permits stating the hypothesis that all sections perpendicular to the ${\textstyle z}$ axis suffer identical deformations contained in its own plane. This way, considering a generic cross section contained in the ${\textstyle x-y}$ plane of any of the structures in Figures 19 and 20, the displacement field of the structure can be perfectly defined if we know the displacement in the directions ${\textstyle x}$ and ${\textstyle y}$ for all points of that section. Thereby, the vector of displacements at any point is defined as:

 ${\displaystyle {\boldsymbol {u}}(x,y)=\left\{{\begin{array}{c}u(x,y)\\v(x,y)\end{array}}\right\}}$
(65)

where ${\textstyle u(x,y)}$ and ${\textstyle v(x,y)}$ are the displacements of the point in the directions ${\textstyle x}$ and ${\textstyle y}$, respectively.

From the displacements field (65) one can easily obtain the strains field following the general elasticity theory [51]:

 ${\displaystyle \varepsilon _{x}}$ ${\displaystyle =\displaystyle {\frac {\partial u}{\partial x}}}$ ${\displaystyle \varepsilon _{y}}$ ${\displaystyle =\displaystyle {\frac {\partial v}{\partial y}}}$ (66) ${\displaystyle \gamma _{xy}}$ ${\displaystyle =\displaystyle {\frac {\partial u}{\partial y}}+{\frac {\partial v}{\partial x}}}$ ${\displaystyle \gamma _{xz}}$ ${\displaystyle =\gamma _{yz}=0}$

In plane strain problems, there is the additional hypothesis that the strain ${\textstyle \varepsilon _{z}}$ is null, whereas in plane stress problems it is the stress ${\textstyle \sigma _{z}}$ that is zero. This makes that, in any case, it is no necessary to consider strain ${\textstyle \varepsilon _{z}}$ because it does not take part in the work of deformations since ${\textstyle \sigma _{z}\varepsilon _{z}=0}$. Thereby, at any point in a plane stress or plane strain state, we can define the vector of significant strains like:

 ${\displaystyle {\boldsymbol {\varepsilon }}=[\varepsilon _{x},\varepsilon _{y},\gamma _{xy}]^{T}}$
(67)

From equations in (66) we can deduce that tangential stresses ${\textstyle \tau _{xz}}$ and ${\textstyle \tau _{yz}}$ are zero. Moreover, for the same reasons of the previous paragraph, stress ${\textstyle \sigma _{z}}$ does not work, and so the significant stress vector can be defined as follows:

 ${\displaystyle {\boldsymbol {\sigma }}=[\sigma _{x},\sigma _{y},\tau _{xy}]^{T}}$
(68)

Once the strains and stresses have been presented, it is important to stablish the relation between them. This relation can be deduced by applying the previously stated hypothesis into the three-dimensional constitutive equation of elasticity [51]. The elastic stress-strain relation can be written like:

 ${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {E}}{\boldsymbol {\varepsilon }}}$
(69)

where ${\textstyle {\boldsymbol {E}}}$ is the elastic constitutive matrix:

 ${\displaystyle {\boldsymbol {E}}={\begin{bmatrix}E_{11}&E_{12}&0\\E_{21}&E_{22}&0\\0&0&E_{33}\end{bmatrix}}}$
(70)

For elastic isotropic materials the values in matrix ${\textstyle {\boldsymbol {E}}}$ are given by:

 ${\displaystyle {\hbox{Plane Stress }}}$ ${\displaystyle \qquad }$ ${\displaystyle {\hbox{Plane Strain }}}$ ${\displaystyle E_{11}=E_{22}={\frac {E}{1-\nu ^{2}}}}$ ${\displaystyle E_{11}=E_{22}={\frac {E(1-\nu )}{(1+\nu )(1-2\nu )}}}$ ${\displaystyle E_{12}=E_{21}=E_{11}\nu }$ ${\displaystyle E_{12}=E_{21}=E_{11}{\frac {\nu }{1-\nu }}}$ (71) ${\displaystyle E_{33}={\frac {E}{2(1+\nu )}}=G}$ ${\displaystyle E_{33}={\frac {E}{2(1+\nu )}}=G}$

with ${\textstyle E}$ being the Young's modulus, and ${\textstyle \nu }$ the Poisson's ratio.

#### 3.2.1.2 Natural coordinates and shape functions

One of the most important concepts of the Finite Element Method are the so called, shape functions. These functions ${\textstyle N}$ are defined at each node of a finite element and allow us to obtain the value of a nodal variable at any point of the element through interpolation.

We will present the expressions of the shape functions for the quadrilateral element of four nodes, and for the triangular element of three nodes, since these are the elements that have been used in this work. To do so, we will first introduce the concept of natural coordinates, that will be used in the expressions of the shape functions.

The natural coordinates ${\textstyle (\xi ,\eta )}$ are normalized coordinates, meaning that they range from -1 to +1, like it is shown in Figure 21. From this same figure it can be deduced that:

 ${\displaystyle \xi ={\frac {x-x_{c}}{a}}\qquad ;\qquad \eta ={\frac {y-y_{c}}{b}}}$
(72)

where ${\textstyle x_{c}}$ and ${\textstyle y_{c}}$ are the coordinates of the center of the element. Moreover:

 ${\displaystyle {\frac {d\xi }{dx}}={\frac {1}{a}}\qquad ;\qquad {\frac {d\eta }{dy}}={\frac {1}{b}}}$
(73)

and so a differential of area is obtained like:

 ${\displaystyle dxdy=abd\xi d\eta }$
(74)

and, as it will appear later on, the integral of a function ${\textstyle f(x,y)}$ over a rectangular element can be perform by means of the following transformation to the natural coordinates:

 ${\displaystyle \iint _{A^{(e)}}f(x,y)dxdy=\int _{-1}^{+1}\int _{-1}^{+1}g(\xi ,\eta )abd\xi d\eta }$
(75)
 Figure 21: Geometry of a generic rectangular element. Natural and Cartesian coordinates.

In order to properly understand the expressions of the shape functions, it is interesting to introduce the complete polynomials in two-dimensions.

The shape functions can reproduce with no error polynomial variations of order less than or equal to the order of the complete polynomial contained in those shape functions. It can be deduced then, that the higher is the order of that complete polynomial, the more accurate can be the finite element solution.

In a 2D natural coordinates system, a complete polynomial of order ${\textstyle n}$ can be written like:

 ${\displaystyle f(\xi ,\eta )=\sum _{i=1}^{p}\alpha _{i}\xi ^{j}\eta ^{k}\qquad ;\qquad j+k\leq n}$
(76)

where the number of terms of the polynomial is:

 ${\displaystyle p={\frac {(n+1)(n+2)}{2}}}$
(77)

This way, a complete linear polynomial (${\textstyle n=1\rightarrow p=3}$) would be:

 ${\displaystyle f(\xi ,\eta )=\alpha _{1}+\alpha _{2}\xi {+\alpha }_{3}\eta }$
(78)

An easy way to identify the terms of a 2D complete polynomial is by using the Pascal's triangle (Figure 22).

 Figure 22: Pascal's triangle in two dimensions.

The shape functions of many elements (e.g. the 4-node quadrilateral) have terms of incomplete polynomials. Those terms generate nodal variables that barely contribute to the improvement of the element accuracy. Therefore, between two elements with shape functions containing complete polynomials of the same order, it is advisable to use that with less nodal variables.

That aside, all shape functions of an element must satisfy the following two conditions:

• Nodal compatibility condition
 ${\displaystyle N_{i}(\xi _{j},\eta _{j})=\left\{{\begin{array}{l}1\quad {\hbox{if }}i=j{\hbox{ }}\\0\quad {\hbox{if }}i\neq j{\hbox{ }}\end{array}}\right.\quad {\hbox{for }}\quad i,j=\{1,2,\dots ,n\}}$
(79)
• Rigid body condition
 ${\displaystyle \sum _{i=1}^{n}N_{i}(\xi ,\eta )=1}$
(80)

where ${\textstyle i}$ and ${\textstyle j}$ represent node indexes, and ${\textstyle n}$ is the number of nodes of the finite element.

At this point, we can already present the expressions of the shape functions, starting with the rectangular Lagrangian element of four nodes.

The shape functions of this element are based on polynomial interpolations of Lagrange in two dimensions. This permits defining the shape function of a generic node like the product between two uni-dimensional Lagrange polynomials in each direction of the two natural coordinates ${\textstyle \xi }$ and ${\textstyle \eta }$ at that node. Thereby, if we consider ${\textstyle l_{I}^{i}(\xi )}$ like the Lagrange polynomial of order ${\textstyle I}$ in direction ${\textstyle \xi }$ at node ${\textstyle i}$, and ${\textstyle l_{J}^{i}(\eta )}$ the one of order ${\textstyle J}$ in direction ${\textstyle \eta }$, the shape function at that node is:

 ${\displaystyle N_{i}(\xi ,\eta )=l_{I}^{i}(\xi )l_{J}^{i}(\eta )}$
(81)

In particular, for the 4-node rectangular Lagrangian element, the uni-dimensional Lagrange polynomials in each direction ${\textstyle \xi }$ and ${\textstyle \eta }$ coincide with the shape functions of the 2-node bar element. Thereby, at a node ${\textstyle i}$ of the element we have:

 ${\displaystyle l_{1}^{i}(\xi )={\frac {1}{2}}(1+\xi \xi _{i})\qquad ;\qquad l_{1}^{i}(\eta )={\frac {1}{2}}(1+\eta \eta _{i})}$
(82)

And consequently, the shape function at node ${\textstyle i}$ results:

 ${\displaystyle N_{i}(\xi ,\eta )=l_{1}^{i}(\xi )l_{1}^{i}(\eta )={\frac {1}{4}}(1+\xi \xi _{i})(1+\eta \eta _{i})}$
(83)

The values of ${\textstyle \xi _{i}}$ and ${\textstyle \eta _{i}}$ can be seen in Figure 21 but are summarized in Table 4.

 Node ${\displaystyle \xi _{i}}$ ${\displaystyle \eta _{i}}$ 1 -1 -1 2 1 -1 3 1 1 4 -1 1

With regard to the triangular element of three nodes, we will present its shape functions also in terms of the natural coordinates, since it is how they are used in KRATOS.

If we define a natural coordinates system over the normalized geometry of a triangular element, like in Figure 23, we can formulate the shape functions for each of its three nodes like:

 ${\displaystyle N_{1}=1-\xi {-\eta }\quad ;\quad N_{2}=\xi \quad ;\quad N_{3}=\eta }$
(84)
 Figure 23: Natural coordinates in a triangular element.

In this form, the triangular formulation becomes very similar to the quadrilateral one since, as it can be seen from figures 21 and 23, the normalized triangle is like a portion of a normalized square, which is very interesting when defining isoparametric elements.

#### 3.2.1.3 Discretization of the displacements, strains and stresses fields

In summary, we have seen that thanks to the shape functions, one can obtain a variable at any point of the geometry through an interpolation from a limited number of nodal results.

Therefore, for a generic two-dimensional element of ${\textstyle n}$ nodes, we can write the displacement field like:

 ${\displaystyle u=\displaystyle \sum _{i=1}^{n}N_{i}u_{i}\quad ;\quad v=\displaystyle \sum _{i=1}^{n}N_{i}v_{i}}$
(85)

where ${\textstyle (u_{i},v_{i})}$ and ${\textstyle N_{i}}$ are the horizontal and vertical displacements and the shape function of node ${\textstyle i}$. Equations in (85) can be rewritten in matrix format like:

 ${\displaystyle {\boldsymbol {u}}={\begin{Bmatrix}u\\v\end{Bmatrix}}={\begin{bmatrix}N_{1}&0&\dots &N_{n}&0\\0&N_{1}&\dots &0&N_{n}\end{bmatrix}}{\begin{Bmatrix}u_{1}\\v_{1}\\\vdots \\u_{n}\\v_{n}\end{Bmatrix}}}$
(86)

or equivalently:

 ${\displaystyle {\boldsymbol {u}}={\boldsymbol {N}}{\boldsymbol {a}}^{(e)}}$
(87)

where ${\textstyle {\boldsymbol {u}}}$ is the vector of displacements at a point of the element,

 ${\displaystyle {\boldsymbol {N}}=[{\boldsymbol {N}}_{1},{\boldsymbol {N}}_{2},\dots ,{\boldsymbol {N}}_{n}]\qquad ;\qquad {\boldsymbol {N}}_{i}={\begin{bmatrix}N_{i}&0\\0&N_{i}\end{bmatrix}}}$
(88)

are the matrix of shape functions of the element, and the matrix of shape functions of node ${\textstyle i}$, respectively, and

 ${\displaystyle {\boldsymbol {a}}^{(e)}={\begin{Bmatrix}{\boldsymbol {a}}_{1}^{(e)}\\{\boldsymbol {a}}_{2}^{(e)}\\\vdots \\{\boldsymbol {a}}_{n}^{(e)}\end{Bmatrix}}\qquad ;\qquad {\boldsymbol {a}}_{i}^{(e)}={\begin{Bmatrix}u_{i}\\v_{i}\end{Bmatrix}}}$
(89)

are the vector of nodal displacements of the element, and the vector of nodal displacements of node ${\textstyle i}$.

In the case of the strain field, we can write:

 ${\displaystyle \varepsilon _{x}}$ ${\displaystyle ={\frac {\partial u}{\partial x}}=\sum _{i=1}^{n}{\frac {\partial N_{i}}{\partial x}}u_{i}}$ ${\displaystyle \varepsilon _{y}}$ ${\displaystyle ={\frac {\partial v}{\partial y}}=\sum _{i=1}^{n}{\frac {\partial N_{i}}{\partial y}}v_{i}}$ (90) ${\displaystyle \gamma _{xy}}$ ${\displaystyle ={\frac {\partial u}{\partial y}}+{\frac {\partial v}{\partial x}}=\sum _{i=1}^{n}\left({\frac {\partial N_{i}}{\partial y}}u_{i}+{\frac {\partial N_{i}}{\partial x}}v_{i}\right)}$

and in matrix format:

 ${\displaystyle {\boldsymbol {\varepsilon }}={\begin{Bmatrix}{\frac {\partial u}{\partial x}}\\{\frac {\partial v}{\partial y}}\\{\frac {\partial u}{\partial y}}+{\frac {\partial v}{\partial x}}\end{Bmatrix}}={\begin{bmatrix}{\frac {\partial N_{1}}{\partial x}}&0&\dots &{\frac {\partial N_{n}}{\partial x}}&0\\0&{\frac {\partial N_{1}}{\partial y}}&\dots &0&{\frac {\partial N_{n}}{\partial y}}\\{\frac {\partial N_{1}}{\partial y}}&{\frac {\partial N_{1}}{\partial x}}&\dots &{\frac {\partial N_{n}}{\partial y}}&{\frac {\partial N_{n}}{\partial x}}\end{bmatrix}}{\begin{Bmatrix}u_{1}\\v_{1}\\\vdots \\u_{n}\\v_{n}\end{Bmatrix}}}$
(91)

or

 ${\displaystyle {\boldsymbol {\varepsilon }}={\boldsymbol {B}}{\boldsymbol {a}}^{(e)}}$
(92)

where

 ${\displaystyle {\boldsymbol {B}}=[{\boldsymbol {B}}_{1},{\boldsymbol {B}}_{2},\dots ,{\boldsymbol {B}}_{n}]}$
(93)

is the deformation matrix of the element, and

 ${\displaystyle {\boldsymbol {B}}_{i}={\begin{bmatrix}\displaystyle {\frac {\partial N_{i}}{\partial x}}&0\\0&\displaystyle {\frac {\partial N_{i}}{\partial y}}\\\displaystyle {\frac {\partial N_{i}}{\partial y}}&\displaystyle {\frac {\partial N_{i}}{\partial x}}\end{bmatrix}}}$
(94)

is the deformation matrix of node ${\textstyle i}$.

Now, substituting equation (92) in the relation (69), we can express the stress vector like:

 ${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {E}}{\boldsymbol {\varepsilon }}={\boldsymbol {E}}{\boldsymbol {B}}{\boldsymbol {a}}^{(e)}}$
(95)

If there were initial strains or stresses, they could be taken into account by using the following expression:

 ${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {E}}({\boldsymbol {\varepsilon }}-{\boldsymbol {\varepsilon }}^{0})+{\boldsymbol {\sigma }}^{0}={\boldsymbol {E}}{\boldsymbol {B}}{\boldsymbol {a}}^{(e)}-{\boldsymbol {E}}{\boldsymbol {\varepsilon }}^{0}+{\boldsymbol {\sigma }}^{0}}$
(96)

#### 3.2.1.4 Equilibrium equations of the discretization

The expressions of the stiffness matrix ${\textstyle {\boldsymbol {K}}}$ and the force vector ${\textstyle {\boldsymbol {f}}}$ can be obtained form of the equilibrium equations of the discretization. To do so, we are going to start from the expression of the Principle of Virtual Works applied at the equilibrium of an element.

Let us suppose that uniformly distributed forces per unit area act over the body of the element (mass forces ${\textstyle {\boldsymbol {b}}}$), and uniformly distributed forces per unit length act over one of its sides (surface forces ${\textstyle {\boldsymbol {t}}}$). Moreover, supposing that the equilibrium of the element is achieved at the nodes, we define punctual forces acting at the nodes (nodal forces of equilibrium ${\textstyle {\boldsymbol {q}}}$) that must balance the forces that appear due to the element deformation as well as the rest of acting forces.

Thereby, we can write the Principle of Virtual Works applied at the element like:

 ${\displaystyle \iint _{A^{(e)}}\delta {\boldsymbol {\varepsilon }}^{T}{\boldsymbol {\sigma }}\;tdA=\iint _{A^{(e)}}\delta {\boldsymbol {u}}^{T}{\boldsymbol {b}}\;tdA+\oint _{l^{(e)}}\delta {\boldsymbol {u}}^{T}{\boldsymbol {t}}\;tds+[\delta {\boldsymbol {a}}^{(e)}]^{T}{\boldsymbol {q}}^{(e)}}$
(97)

The first term of the expression represents the work of the stresses over the virtual strains ${\textstyle \delta {\boldsymbol {\varepsilon }}}$, whereas the terms on the right represent the work of the mass forces ${\textstyle {\boldsymbol {b}}}$, the surface forces ${\textstyle {\boldsymbol {t}}}$ and the punctual forces ${\textstyle {\boldsymbol {q}}^{(e)}}$ over the virtual displacements ${\textstyle \delta {\boldsymbol {u}}}$ and ${\textstyle \delta {\boldsymbol {a}}^{(e)}}$. ${\textstyle A^{(e)}}$ and ${\textstyle l^{(e)}}$ are the area and the contour of the element, and ${\textstyle t}$ is the thickness. In plane stress problems, ${\textstyle t}$ is the real thickness of the structure, while in plane strain ${\textstyle t}$ is usually taken as 1.

From (87) and (92), we can write:

 ${\displaystyle \delta {\boldsymbol {u}}^{T}=[\delta {\boldsymbol {a}}^{(e)}]^{T}{\boldsymbol {N}}^{T}\quad ;\quad \delta {\boldsymbol {\varepsilon }}^{T}=[\delta {\boldsymbol {a}}^{(e)}]^{T}{\boldsymbol {B}}^{T}}$
(98)

Substituting (98) in (97), and rearranging terms, we obtain:

 ${\displaystyle [\delta {\boldsymbol {a}}^{(e)}]^{T}\left[\iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}\;tdA-\iint _{A^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {b}}\;tdA-\oint _{l^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {t}}\;tds\right]=[\delta {\boldsymbol {a}}^{(e)}]^{T}{\boldsymbol {q}}^{(e)}}$
(99)

Since the virtual displacements ${\textstyle \delta {\boldsymbol {a}}^{(e)}}$ are arbitrary, we can write:

 ${\displaystyle \iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}\;tdA-\iint _{A^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {b}}\;tdA-\oint _{l^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {t}}\;tds={\boldsymbol {q}}^{(e)}}$
(100)

Equation (100) represents the equilibrium between the nodal forces of equilibrium ${\textstyle {\boldsymbol {q}}^{(e)}}$ and the forces from the deformation of the element (first integral), the mass forces ${\textstyle {\boldsymbol {b}}}$ and the surface forces ${\textstyle {\boldsymbol {t}}}$. If we now substitute the stresses from expression (96), we get:

 ${\displaystyle \iint _{A^{(e)}}{\boldsymbol {B}}^{T}({\boldsymbol {E}}{\boldsymbol {B}}{\boldsymbol {a}}^{(e)}-{\boldsymbol {E}}{\boldsymbol {\varepsilon }}^{0}+{\boldsymbol {\sigma }}^{0})\;tdA-\iint _{A^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {b}}\;tdA-\oint _{l^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {t}}\;tds={\boldsymbol {q}}^{(e)}}$
(101)

and rearranging terms:

 ${\displaystyle \left[\iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {B}}\;tdA\right]{\boldsymbol {a}}^{(e)}-\iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {\varepsilon }}^{0}\;tdA+}$ ${\displaystyle +\iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}^{0}\;tdA-\iint _{A^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {b}}\;tdA-\oint _{l^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {t}}\;tds={\boldsymbol {q}}^{(e)}}$
(102)

which can be expressed like:

 ${\displaystyle {\boldsymbol {K}}^{(e)}{\boldsymbol {a}}^{(e)}-{\boldsymbol {f}}^{(e)}={\boldsymbol {q}}^{(e)}}$
(103)

where:

 ${\displaystyle {\boldsymbol {K}}^{(e)}=\iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {B}}\;tdA}$
(104)

is the elastic stiffness matrix of the element, and:

 ${\displaystyle {\boldsymbol {f}}^{(e)}={\boldsymbol {f}}_{\varepsilon }^{(e)}+{\boldsymbol {f}}_{\sigma }^{(e)}+{\boldsymbol {f}}_{b}^{(e)}+{\boldsymbol {f}}_{t}^{(e)}}$
(105)

is the vector of equivalent nodal forces of the element, with:

 ${\displaystyle {\boldsymbol {f}}_{\varepsilon }^{(e)}=\iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {\varepsilon }}^{0}\;tdA\quad ;\quad {\boldsymbol {f}}_{\sigma }^{(e)}=-\iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {\sigma }}^{0}\;tdA}$
(106)

being the vectors of equivalent nodal forces due to initial deformations and initial stresses, respectively;

 ${\displaystyle {\boldsymbol {f}}_{b}^{(e)}=\iint _{A^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {b}}\;tdA}$
(107)

is the vector of equivalent nodal forces due to uniformly distributed forces per unit area; and

 ${\displaystyle {\boldsymbol {f}}_{t}^{(e)}=\oint _{l^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {t}}\;tds}$
(108)

is the vector of equivalent nodal forces due to uniformly distributed forces per unit length.

The global equilibrium equation of the mesh is obtained by imposing that the sum of nodal forces of equilibrium at each node must equal the external nodal force:

 ${\displaystyle \sum _{e}{\boldsymbol {q}}_{i}^{(e)}={\boldsymbol {p}}_{j}}$
(109)

where the left part of the expression represents the sum of contributions of the vectors of nodal forces of equilibrium of the different elements that share the global node ${\textstyle j}$, and ${\textstyle {\boldsymbol {p}}_{j}}$ represents the vector of external punctual forces acting on such node.

Therefore, the global equilibrium equation of the mesh can be obtained by assembling the contributions of stiffness matrices and vectors of equivalent nodal forces of the different elements. The global matrix equation can be written like:

 ${\displaystyle {\boldsymbol {K}}{\boldsymbol {a}}={\boldsymbol {f}}}$
(110)

where ${\textstyle {\boldsymbol {K}}}$, ${\textstyle {\boldsymbol {a}}}$ and ${\textstyle {\boldsymbol {f}}}$ are, respectively, the stiffness matrix, the vector of nodal displacements, and the vector of equivalent nodal forces of the whole mesh.

#### 3.2.1.5 Isoparametric elements and numerical integration

At this point, we have already obtained the equilibrium equations of the discretization, but we have not presented yet a convenient numerical technique to integrate the expressions of the stiffness matrix and the vectors of equivalent nodal forces of the elements.

The integration of the expressions of the stiffness matrix and the vectors of equivalent nodal forces, will be performed by means of Gauss-Legendre quadratures. This technique allow us to integrate any function over a normalized domain, using tabulated Gauss point coordinates and weights. However, we need to transform first the integrals over the element domain into integrals over the normalized space of the natural coordinates.

In order to do so, we must first introduce a very important concept in the Finite Element Method: the isoparametric interpolation. In this regard, we will present the isoparametric quadrilateral elements, and the isoparametric triangular elements.

The concept of isoparametric interpolation comes from the usage of the same shape functions to interpolate the geometry and the displacements. Thereby, just as we saw in (85), we can express the geometry of a two-dimensional isoparametric element from the coordinates ${\textstyle x_{i}}$ and ${\textstyle y_{i}}$ of its nodes like:

 ${\displaystyle x=\sum _{i=1}^{n}N_{i}(\xi ,\eta )x_{i}\quad ;\quad y=\sum _{i=1}^{n}N_{i}(\xi ,\eta )y_{i}}$
(111)

where ${\textstyle N_{i}(\xi ,\eta )}$ are the element shape functions. Equations in (111) relate the cartesian coordinates of a point with the natural coordinates ${\textstyle \xi }$ and ${\textstyle \eta }$, which is very convenient for the transformation of the integration domain.

In this regard, in order to formulate the change of variables from cartesian coordinates to natural coordinates, we need to compute the determinant of the Jacobian matrix. This can be obtained by deriving the shape functions with respect to ${\textstyle \xi }$ and ${\textstyle \eta }$. If we consider ${\textstyle N_{i}(\xi ,\eta )=N_{i}(x(\xi ,\eta ),y(\xi ,\eta ))}$ then, applying the chain rule, we have:

 ${\displaystyle {\begin{array}{l}\displaystyle {\frac {\partial N_{i}}{\partial \xi }}=\displaystyle {\frac {\partial N_{i}}{\partial x}}{\frac {\partial x}{\partial \xi }}+{\frac {\partial N_{i}}{\partial y}}{\frac {\partial y}{\partial \xi }}\\\\\displaystyle {\frac {\partial N_{i}}{\partial \eta }}=\displaystyle {\frac {\partial N_{i}}{\partial x}}{\frac {\partial x}{\partial \eta }}+{\frac {\partial N_{i}}{\partial y}}{\frac {\partial y}{\partial \eta }}\end{array}}}$
(112)

or in matrix format:

 ${\displaystyle {\begin{Bmatrix}{\frac {\partial N_{i}}{\partial \xi }}\\{\frac {\partial N_{i}}{\partial \eta }}\end{Bmatrix}}={\begin{bmatrix}{\frac {\partial x}{\partial \xi }}&{\frac {\partial y}{\partial \xi }}\\{\frac {\partial x}{\partial \eta }}&{\frac {\partial y}{\partial \eta }}\end{bmatrix}}{\begin{Bmatrix}{\frac {\partial N_{i}}{\partial x}}\\{\frac {\partial N_{i}}{\partial y}}\end{Bmatrix}}={\boldsymbol {J}}^{(e)}{\begin{Bmatrix}{\frac {\partial N_{i}}{\partial x}}\\{\frac {\partial N_{i}}{\partial y}}\end{Bmatrix}}}$
(113)

where ${\textstyle {\boldsymbol {J}}^{(e)}}$ is the Jacobian matrix of the transformation of natural coordinates to cartesian coordinates. From (113) we can derive:

 ${\displaystyle {\begin{Bmatrix}{\frac {\partial N_{i}}{\partial x}}\\{\frac {\partial N_{i}}{\partial y}}\end{Bmatrix}}=\left[{\boldsymbol {J^{(e)}}}\right]^{-1}{\begin{Bmatrix}{\frac {\partial N_{i}}{\partial \xi }}\\{\frac {\partial N_{i}}{\partial \eta }}\end{Bmatrix}}={\frac {1}{\left|{\boldsymbol {J}}^{(e)}\right|}}{\begin{bmatrix}{\frac {\partial y}{\partial \eta }}&-{\frac {\partial y}{\partial \xi }}\\-{\frac {\partial x}{\partial \eta }}&{\frac {\partial x}{\partial \xi }}\end{bmatrix}}{\begin{Bmatrix}{\frac {\partial N_{i}}{\partial \xi }}\\{\frac {\partial N_{i}}{\partial \eta }}\end{Bmatrix}}}$
(114)

where ${\textstyle \left|{\boldsymbol {J}}^{(e)}\right|}$ is the determinant of the Jacobian.

The determinant of the Jacobian permits expressing the differential of area in natural coordinates [52] like:

 ${\displaystyle dx\;dy=\left|{\boldsymbol {J}}^{(e)}\right|d\xi \;d\eta }$
(115)

Now, using the isoparametric transformation (111) we can obtain the terms of the Jacobian:

 ${\displaystyle {\frac {\partial x}{\partial \xi }}=\sum _{i=1}^{n}{\frac {\partial N_{i}}{\partial \xi }}x_{i}\quad ;\quad {\frac {\partial x}{\partial \eta }}=\sum _{i=1}^{n}{\frac {\partial N_{i}}{\partial \eta }}x_{i}\quad ;\quad {\hbox{etc. }}}$

and so:

 ${\displaystyle {\boldsymbol {J}}^{(e)}=\sum _{i=1}^{n}{\begin{bmatrix}{\frac {\partial N_{i}}{\partial \xi }}x_{i}&{\frac {\partial N_{i}}{\partial \xi }}y_{i}\\{\frac {\partial N_{i}}{\partial \eta }}x_{i}&{\frac {\partial N_{i}}{\partial \eta }}y_{i}\end{bmatrix}}}$
(116)

Furthermore, substituting the expressions of equation (114) in (94) we obtain the deformation matrix of an isoparametric element in terms of the natural coordinates:

 ${\displaystyle {\boldsymbol {B}}_{i}(\xi ,\eta )={\frac {1}{\left|{\boldsymbol {J}}^{(e)}\right|}}{\begin{bmatrix}b_{i}&0\\0&c_{i}\\c_{i}&b_{i}\end{bmatrix}}}$
(117)

where

 ${\displaystyle b_{i}={\frac {\partial y}{\partial \eta }}{\frac {N_{i}}{\partial \xi }}-{\frac {\partial y}{\partial \xi }}{\frac {\partial N_{i}}{\partial \eta }}\quad ;\quad c_{i}={\frac {\partial x}{\partial \xi }}{\frac {N_{i}}{\partial \eta }}-{\frac {\partial x}{\partial \eta }}{\frac {\partial N_{i}}{\partial \xi }}}$
(118)

Finally, by applying the previous expressions, the elastic stiffness matrix (104) can be written like an integral over the normalized domain of the natural coordinates.

 ${\displaystyle {\boldsymbol {K}}^{(e)}=\iint _{A^{(e)}}{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {B}}\;tdxdy=\int _{-1}^{+1}\int _{-1}^{+1}{\boldsymbol {B}}^{T}(\xi ,\eta ){\boldsymbol {E}}{\boldsymbol {B}}(\xi ,\eta )\left|{\boldsymbol {J}}^{(e)}\right|\;td\xi d\eta }$
(119)

And for triangular elements we must consider the different domain of integration (see Figure 23):

 ${\displaystyle {\boldsymbol {K}}^{(e)}=\int _{0}^{1}\int _{0}^{1-\eta }{\boldsymbol {B}}^{T}(\xi ,\eta ){\boldsymbol {E}}{\boldsymbol {B}}(\xi ,\eta )\left|{\boldsymbol {J}}^{(e)}\right|\;td\xi d\eta }$
(120)

At this point we can already use the Gauss quadrature to solve the integrals of the previous expressions. First, though, let us remind the fundamentals of the Gauss-Legendre quadrature for rectangular and triangular domains.

The integral of a generic function ${\textstyle g(\xi ,\eta )}$ over the domain of natural coordinates of a quadrilateral element can be evaluated by means of a Gauss quadrature like:

 ${\displaystyle \int _{-1}^{+1}\int _{-1}^{+1}g(\xi ,\eta )d\xi d\eta =\int _{-1}^{+1}d\xi \left[\sum _{q=1}^{n_{q}}g(\xi ,\eta _{q})W_{q}\right]=\sum _{p=1}^{n_{p}}\sum _{q=1}^{n_{q}}g(\xi _{p},\eta _{q})W_{p}W_{q}}$
(121)

where ${\textstyle n_{p}}$ and ${\textstyle n_{q}}$ are the number of integration points in each direction ${\textstyle \xi }$ and ${\textstyle \eta }$; ${\textstyle \xi _{p}}$ and ${\textstyle \eta _{q}}$ are the natural coordinates of the integration points ${\textstyle p}$, ${\textstyle q}$ and ${\textstyle W_{p}}$, ${\textstyle W_{q}}$ the weights of that point corresponding to each direction.

On the other hand, the Gauss quadrature for a triangular element can be written like:

 ${\displaystyle \int _{0}^{1}\int _{0}^{1-\eta }g(\xi ,\eta )d\xi d\eta =\sum _{p=1}^{n_{p}}g(\xi _{p},\eta _{p})W_{p}}$
(122)

The coordinates and weights of the integration points are tabulated. And the number of integration points must be properly chosen, taking into account that a Gauss quadrature of order ${\textstyle n}$ in each natural direction permits integrating with no error a polynomial of order ${\textstyle \leq 2n-1}$ in the corresponding natural coordinate. In our case, we have been using 2x2 integration points for the quadrilateral elements of four nodes, and 1 integration point for the triangular elements of three nodes.

Now, substituting (121) in (119), we can obtain the expression of the elastic stiffness matrix for an isoparametric quadrilateral element, evaluated through numerical integration:

 ${\displaystyle {\boldsymbol {K}}^{(e)}=\int _{-1}^{+1}\int _{-1}^{+1}{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {B}}\left|{\boldsymbol {J}}^{(e)}\right|\;td\xi d\eta =\sum _{p=1}^{n_{p}}\sum _{q=1}^{n_{q}}\left[{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {B}}\left|{\boldsymbol {J}}^{(e)}\right|t\right]_{p,q}W_{p}W_{q}}$
(123)

And from (120) and (122) we obtain the expression corresponding to a triangular element:

 ${\displaystyle {\boldsymbol {K}}^{(e)}=\int _{0}^{1}\int _{0}^{1-\eta }{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {B}}\left|{\boldsymbol {J}}^{(e)}\right|\;td\xi d\eta =\sum _{p=1}^{n_{p}}\left[{\boldsymbol {B}}^{T}{\boldsymbol {E}}{\boldsymbol {B}}\left|{\boldsymbol {J}}^{(e)}\right|t\right]_{p}W_{p}}$
(124)

Thereby, we can see that the numerical integration of the stiffness matrix requires the evaluation of the Jacobian ${\textstyle {\boldsymbol {J}}^{(e)}}$, its determinant ${\textstyle \left|{\boldsymbol {J}}^{(e)}\right|}$, the deformation matrix ${\textstyle {\boldsymbol {B}}}$, the constitutive matrix ${\textstyle {\boldsymbol {E}}}$, and the thickness ${\textstyle t}$ at each integration point.

Similarly, the computation of any of the vectors of equivalent nodal forces that imply integrals over a quadrilateral element, for instance (107), can be performed like:

 ${\displaystyle {\boldsymbol {f}}_{b}^{(e)}}$ ${\displaystyle =\iint _{A^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {b}}tdxdy=\int _{-1}^{+1}\int _{-1}^{+1}{\boldsymbol {N}}^{T}{\boldsymbol {b}}\left|{\boldsymbol {J}}^{(e)}\right|td\xi d\eta }$ ${\displaystyle =\sum _{p=1}^{n_{p}}\sum _{q=1}^{n_{q}}\left[{\boldsymbol {N}}^{T}{\boldsymbol {b}}\left|{\boldsymbol {J}}^{(e)}\right|t\right]_{p,q}W_{p}W_{q}}$
(125)

For triangular elements the double sum is replaced by the simple sum of (122).

Finally, the computation of the vector of surface forces (108) is a bit different because the integral is performed over the contour of the element ${\textstyle l^{(e)}}$. In general, such contour represents a straight line ${\textstyle \xi =}$const or ${\textstyle \eta =}$const in the space of natural coordinates. Therefore, for example, for a contour of an isoparametric quadrilateral element corresponding to the straight line ${\textstyle \eta =1}$, the differential of length ${\textstyle ds}$ is computed like:

 ${\displaystyle [ds]_{\eta =1}}$ ${\displaystyle =\left[{\sqrt {dx^{2}+dy^{2}}}\right]_{\eta =1}=\left[{\sqrt {\left({\frac {dx}{d\xi }}\right)_{\eta =1}^{2}+\left({\frac {dy}{d\xi }}\right)_{\eta =1}^{2}}}\right]d\xi }$ ${\displaystyle =\left[{\sqrt {\left(\sum _{i=1}^{n}{\frac {dN_{i}}{d\xi }}x_{i}\right)_{\eta =1}^{2}+\left(\sum _{i=1}^{n}{\frac {dN_{i}}{d\xi }}y_{i}\right)_{\eta =1}^{2}}}\right]d\xi =c(\xi )d\xi }$
(126)

Substituting (126) in (108) we obtain a line integral that depends only on the natural coordinate ${\textstyle \xi }$ and that can be computed through a uni-dimensional quadrature like:

 ${\displaystyle {\boldsymbol {f}}_{t}^{(e)}=\oint _{l^{(e)}}\left[{\boldsymbol {N}}^{T}\right]_{\eta =1}{\boldsymbol {t}}\;t\;c(\xi )d\xi =\int _{-1}^{+1}{\boldsymbol {g}}(\xi )d\xi =\sum _{p=1}^{n_{p}}{\boldsymbol {g}}(\xi _{p})W_{p}}$
(127)

A more extended overview on the Finite Element Method can be found in [53,54].

## 3.3 Non-linearity Associated to Damage

In the previous section we reviewed fundamental concepts on the Finite Element Method centering our attention on the two-dimensional elasticity theory, which is the basis of the implemented damage code.

In this regard, when deriving the equilibrium equations of the discretization from the Principle of Virtual Works (97), we used an elastic stress-strain relation (69). Nonetheless, in the previous chapter we explained that in order to take into account the lost of stiffness due to damage progression, a different stress-strain relation had to be used.

If we now consider the stress-strain law (13) in voigt notation:

 ${\displaystyle {\boldsymbol {\sigma }}={\boldsymbol {E}}_{sec}{\boldsymbol {\varepsilon }}=(1-d){\boldsymbol {E}}{\boldsymbol {\varepsilon }}}$
(128)

and we introduce it into a general expression of the PVW (100), we obtain the following relation:

 ${\displaystyle \left[\iint _{A^{(e)}}{\boldsymbol {B}}^{T}(1-d){\boldsymbol {E}}{\boldsymbol {B}}tdA\right]{\boldsymbol {a}}^{(e)}-\iint _{A^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {b}}\;tdA-\oint _{l^{(e)}}{\boldsymbol {N}}^{T}{\boldsymbol {t}}\;tds={\boldsymbol {q}}^{(e)}}$
(129)

where ${\textstyle d}$ is the damage parameter.

Expression (129) can be rewritten like:

 ${\displaystyle {\boldsymbol {K}}^{(e)}(d)\,{\boldsymbol {a}}^{(e)}-{\boldsymbol {f}}^{(e)}={\boldsymbol {q}}^{(e)}}$
(130)

where

 ${\displaystyle {\boldsymbol {K}}^{(e)}(d)=\iint _{A^{(e)}}{\boldsymbol {B}}^{T}(1-d){\boldsymbol {E}}{\boldsymbol {B}}tdA}$
(131)

is the damaged or secant stiffness matrix of the element, and

 ${\displaystyle {\boldsymbol {f}}^{(e)}={\boldsymbol {f}}_{b}^{(e)}+{\boldsymbol {f}}_{t}^{(e)}}$
(132)

is the vector of equivalent nodal forces of the element, in which ${\textstyle {\boldsymbol {f}}_{b}^{(e)}}$ and ${\textstyle {\boldsymbol {f}}_{t}^{(e)}}$ coincide with (107) and (108), respectively.

Like we stated before, by assembling the contributions of the different elements of the mesh we can obtain the global matrix equation:

 ${\displaystyle {\boldsymbol {K}}(d)\,{\boldsymbol {a}}={\boldsymbol {f}}}$
(133)

where ${\textstyle {\boldsymbol {K}}(d)}$ is the global secant stiffness matrix ${\textstyle {\boldsymbol {K}}_{sec}}$.

By the way, we saw in Chapter 2 that ${\textstyle d=g(r)}$, where ${\textstyle r}$ is the internal state variable depending on the equivalent strains ${\textstyle \varepsilon _{eq}}$. Therefore, since the strains are directly related to the displacements through (92), we can say that damage also depends on the displacements, i.e. ${\textstyle d=f({\boldsymbol {a}})}$. Taking into account this last consideration, expression (133) results:

 ${\displaystyle {\boldsymbol {K}}_{sec}({\boldsymbol {a}})\,{\boldsymbol {a}}={\boldsymbol {f}}}$
(134)

which is a non-linear system of equations.

Thereby, the inclusion of damage mechanics in the classical elasticity theory introduces a non-linearity that must be taken into account when solving this kind of problems. In this section we are going to review how we solve such non-linear system of equations, but before let us introduce the fundamentals of non-linear problems that will help us understand the results of the simulations performed.

### 3.3.1 Fundamentals of non-linear problems

Let us begin by introducing the term response as a pictorial characterization of non-linearity of a structural system. Response is a graphical representation of the fundamental concept of equilibrium path. Through this representation many concepts can be illustrated and interpreted in physical, mathematical or computational terms.

In this regard, a solid mechanics problem is said to be non-linear when the response diagram is not proportional, i.e. when it does not follow Hooke's law.

The most widely used form of these pictures is the load-deflection response diagram. Given an action parameter (e.g. the applied load) and a response parameter (e.g. the displacement), one can plot a response diagram for the equilibrium path of any system. In fact, real problems have multidimensional response diagrams but, in order to make it understandable, they are usually represented from two representative variables. Figure 24 shows two possible non-linear response diagrams.

 Figure 24: Response diagrams: (a)Load-deflection diagram showing equilibrium path. (b)Diagram distinguishing primary from secondary equilibrium path.

#### 3.3.1.1 Terminology

Before going into the different kinds of non-linear behaviour, let us briefly introduce some basic terminology concerning response diagrams.

The continuous curve shown in a load-deflection diagram is called a path. Paths are piecewise smooth, that is, they have a continuous tangent except at some exceptional points.

Each point in the path represents a possible configuration state of the structure. If the path represents configurations in static equilibrium it is called equilibrium path.

The origin of the response diagram (zero load and zero deflection) is called the reference state because it is the configuration from which loads and deflections are measured.

For ideal cases, the reference state is unstressed and undeformed, and it is also an equilibrium state. This means that an equilibrium path passes through the reference state like in Figure 24.a.

The path that crosses the reference state is called the fundamental equilibrium path, fundamental path or primary path. The fundamental path extends from the reference state up to special states called critical points. Any path that is not a fundamental path but connects with it at a critical point is called a secondary equilibrium path or secondary path (see Figure 24.b).

Qualifiers "primary" and "secondary" are linked with the relative importance of these equilibrium paths in design. Most structures are designed to operate in the fundamental path, with some sort of safety factor against reaching a critical point. However, knowledge of secondary paths may be important in some aspects of the design process, for example in the assessment of structural behaviour under emergency scenarios.

Apart from the stated terminology, there are a pair of additional concepts involving the response diagram that must be pointed out.

The tangent to an equilibrium path may be informally viewed as the limit of the ratio between a force increment and a displacement increment. This is by definition the tangent stiffness associated to the representative force-displacement diagram.

The sign of the tangent stiffness is closely related to the question of stability of an equilibrium state. A negative stiffness is necessarily associated to unstable equilibrium, whereas a positive stiffness is necessary but not sufficient for stability.

If the load and deflection quantities are conjugate in the virtual work sense, the area under a load-deflection diagram may be interpreted as work performed by the system or energy spent in the deformation process.

#### 3.3.1.2 Special equilibrium points

As we have seen, certain points of an equilibrium path have special significance in the applications and thus receive special names. Of particular interest are critical, turning, and failure points.

Critical points can be classified in two groups: limit points, at which the tangent to the equilibrium path is horizontal; and bifurcation points, at which two or more equilibrium paths cross. At critical points the relation between the characteristic load and the associated deflection is not unique. At these points, the structure becomes physically uncontrollable, and so they have engineering significance from a design perspective.

Turning points are states at which the tangent to the equilibrium path is vertical. These are not critical points and have less physical significance, although they can be of interest in connection with the so-called "snap-back" phenomenon. However, turning points may have computational significance because they can affect the performance of certain "path following" solution methods.

Finally, points at which a path suddenly stops or breaks because of physical failure are called failure points. The phenomenon of failure may be local or global in nature. In the first case the structure may regain functional equilibrium after dynamically jumping to another equilibrium path. In the latter case the failure is catastrophic and so the structure cannot regain functional equilibrium.

#### 3.3.1.3 Linear response

A few lines above, we stated that a system is considered non-linear when the response is not proportional. In order to properly understand what a non-linear problem implies, let us first present what a linear response is.

A linear system is a mathematical model characterized by a linear fundamental equilibrium path for all possible choices of load and deflection variables. The consequences of such behaviour are the following:

• A linear structure can sustain any load level and undergo any displacement magnitude.
• There are no critical, turning or failure points.
• Response to different load systems can be obtained by superposition.
• Removing all loads returns the structure to the reference position

The requirements for such a model to be applicable are:

• Perfect linear elasticity for any deformation
• Infinitesimal deformations
• Infinite strength

These assumptions are not only physically unrealistic but mutually contradictory, and thus there are necessarily limits placed on the validity of a linear model. Nevertheless, the linear model can be a good approximation of portions of the non-linear response, for instance, the fundamental path in the vicinity of the reference state. Thereby, since for many structures this segment represents the operational range, the linear model is widely used in design calculations.

#### 3.3.1.4 Response behaviours

Once the linear regime is abandoned, the response diagram may present different basic types of behaviour. Let us present some interesting cases.

Figure 25 illustrates three monotonic types of response: linear, hardening, and softening. Symbols ${\textstyle R}$, ${\textstyle F}$ and ${\textstyle L}$ identify reference, failure and limit points, respectively.

 Figure 25: Basic behaviours for a non-linear response. (a)Linear until brittle failure. (b)Hardening. (c)Softening.

The response shown in Figure 25.a: linear until fracture, is characteristic of pure crystals, glassy, as well as certain high strength composite materials that contain such materials as fibers.

The response illustrated by Figure 25.b is typical of cable, netted and pneumatic (inflatable) structures. The stiffening effect comes from geometry adaptation to the applied loads.

A response such as in Figure 25.b is more common for structure materials, like concrete or steel, than the previous two. A linear branch is followed by a softening regime that may occur suddenly or gradually.

The diagrams of Figure 26 show a combination of basic behaviours that can complicate the response of the structure. Here ${\textstyle B}$ and ${\textstyle T}$ denote bifurcation and turning points.

 Figure 26: Complex non-linear responses. (a)Snap-through. (b)Snap-back. (c)Bifurcation. (d)Bifurcation with snap-back.

The snap-through response in Figure 26.a combines softening with hardening following the second limit point. The response branch between the two limit points has a negative stiffness and is therefore unstable. A response of this type is typical of slightly curved structures such as shallow arches.

The snap-back in Figure 26.b is an exaggerated snap-through, in which the response curve turns back in itself at the turning points. The equilibrium between the two turning points may be stable and physically realizable. This type of response is exhibited by folded and thin shell structures in which moving arch effects occur following the first limit point.

In the previous diagrams the response was a unique curve. Nevertheless, at bifurcation points of Figures 26.b and 26.c more than one response path is possible. The structure takes the path that is dynamically preferred over the others, i.e., the path that implies spending less energy. Bifurcation points can appear in any sufficient thin structure that experiences compressive stresses.

#### 3.3.1.5 Sources of non-linearities

A response diagram characterizes only the gross behaviour of a structure, as it could be observed by performing an experiment on a mechanical testing machine. Further insight into the source of non-linearity is required to capture such physical behaviour with mathematical and computational models.

For structural analysis we can encounter four sources of non-linear behaviour: geometric, material, force boundary conditions, and displacement boundary conditions.

Geometric non-linearities appear when the difference between the deformed and undeformed configurations is taken into account when setting up the strain-displacement and equilibrium equations. Examples of geometric non-linearities include: slender structures in aerospace, civil and mechanical engineering applications; tensile structures such as cables and inflatable membranes; metal and plastic forming, etc.

Material non-linearity is the case of systems whose constitutive behaviour depends on the current deformation state and possibly past history of the deformation. Other constitutive variables may be involved. We can find material non-linearity in structures undergoing non-linear elasticity, plasticity, visco-elasticity, damage, creep, etc.

Non-linear force boundary conditions are applied forces that depend on the deformation. The most important engineering applications concerns pressure loads of fluids. These include hydrostatic loads on submerged or container structures; aerodynamic and hydrodynamic loads (wind loads, wave loads, drag forces). Also gyroscopic and non-conservative followers forces are of mathematical interest.

Lastly, displacement boundary conditions non-linearities come from displacement boundary conditions that depend on the deformation of the structure. The most well-known application is the contact problem, in which no-interpenetration conditions are enforced on flexible bodies while the extent of the contact area is unknown. Non-structural applications of this problem include: ice melting, phase changes, flow in porous media, etc.

Our damage mechanics problem falls into the material non-linearity category, and we will be working with conservative loads and constant boundary conditions.

### 3.3.2 Strategies for the solution of non-linear problems

After reviewing the different response behaviours of non-linearity as well as the sources behind it, one can clearly see that it is not easy to deal with non-linear problems. Several difficulties may arise due to a number of reasons.

The proper definition of the non-linear problem is one of these reasons. Indeed, most of the non-linear models are build "neglecting" many of the non-linear effects presents in the problem (coupling between various phenomena, geometric non-linearity, material non-linearities, contact, etc.).

Another difficulty comes from the lack of uniqueness. In non-linear problems existence and uniqueness of the solution is not guaranteed. Neither it is possible to estimate the computational cost of finding one or more solutions.

Besides, even when one finds a possible solution, the reliability of that solution depends on the robustness of the algorithm used. Usually, parametric studies are necessary to verify the solution (type and number of finite elements used, size of the time steps, boundary conditions, etc.). Furthermore, numerical problems appear when the original problem is not sufficiently well posed, that is, it is not correctly formulated.

Taking into account these difficulties, let us present one of the most widely used strategies to solve non-linear problems like the one in this work.

#### 3.3.2.1 Incremental-iterative methods

Because of the lack of uniqueness of the solution, in many cases the correct solution is path dependent and so it depends on the path followed to reach a given equilibrium state.

In this regard, in many cases it is advisable to follow the physics of the problem and consider the solution, not only as the response to a given action, but as a full historical sequence of the successive states of equilibrium that go from the reference state to another one.

Apart from the stated above, it is interesting to follow the full history of a non-linear process because it gives more information on the mechanical behaviour of the system (engineering reason); and it also helps tracing the equilibrium path near critical points and facilitates convergence (mathematical reason).

Thereby, the method for solving non-linear solid mechanic problems consists on following the equilibrium path by using incrementation or continuity strategies. Assuming the action and the response is known at a given equilibrium state, a new equilibrium state is sought, located on the same branch of equilibrium path and a certain distance from the previous one.

Three different strategies for advance control are typically used:

1. Control of the action or force control: The action (force) is incremented and the corresponding incremental change in the response (displacement) is computed. This is the more commonly used method because of its simplicity, but it is unable to overcome limit points. (Figure 27.a)
2. Control of the response or displacement control: The response (displacement) is incremented and the corresponding incremental change in the action (force) is computed. It allows to overcome limit points, but it cannot deal with turning points. (Figure 27.b)
3. Mixed control or arc-length control: A nearby point on the equilibrium path is sought at a given distance from the previous equilibrium point, measured with a certain norm in the mixed space of action-response (load-displacement). This method is very general and allows overcoming limits and turning points, although its effectiveness depends on the selection of the appropriate distance norm. (Figure 27.c)

Figure 27 displays the geometric representation of the three strategies for advance control presented. In this figure, S represents the last solution point, P is the predicted point, and C is the converged solution after a corrective process.

 Figure 27: Geometric representations of different strategies for advance control: (a)Force control. (b)Displacement control. (c)Arc-length control

In order to clearly explain the incremental-iterative concept, let us formulate our case using the simple force control strategy as example.

The non-linear system of equations obtained in a damage mechanics problem (134) can be expressed like the residual between the internal and the external forces. Thereby, at step ${\textstyle n}$ of a force control strategy we can write:

 ${\displaystyle {\boldsymbol {R}}(^{n}{\boldsymbol {a}})={\boldsymbol {f}}_{int}(^{n}{\boldsymbol {a}})-{^{n}{\boldsymbol {f}}}_{ext}={\boldsymbol {0}}}$
(135)

where

 ${\displaystyle {\boldsymbol {f}}_{int}(^{n}{\boldsymbol {a}})={\boldsymbol {K}}_{sec}(^{n}{\boldsymbol {a}})\,^{n}{\boldsymbol {a}}}$
(136)

and ${\textstyle ^{n}{\boldsymbol {f}}_{ext}}$ is the vector of external forces at the load step ${\textstyle n}$, i.e. ${\textstyle ^{n}{\boldsymbol {f}}_{ext}\equiv {^{n}{\boldsymbol {f}}}}$.

Supposing that step ${\textstyle n}$ is in equilibrium, we must first increase the external load by a prescribed ${\textstyle {^{n+1}\Delta }{\boldsymbol {f}}_{ext}}$. Consequently, the new external force is:

 ${\displaystyle ^{n+1}{\boldsymbol {f}}_{ext}={^{n}{\boldsymbol {f}}}_{ext}+{^{n+1}\Delta }{\boldsymbol {f}}_{ext}}$
(137)

The unknown is the incremental displacement ${\textstyle ^{n+1}\Delta {\boldsymbol {a}}}$ that results from the previous force increment. Thereby:

 ${\displaystyle ^{n+1}{\boldsymbol {a}}={^{n}{\boldsymbol {a}}}+{^{n+1}\Delta }{\boldsymbol {a}}}$
(138)

Imposing equilibrium at the sought solution we obtain the equation that must be solved:

 ${\displaystyle {\boldsymbol {0}}}$ ${\displaystyle ={\boldsymbol {R}}(^{n+1}{\boldsymbol {a}})={\boldsymbol {f}}_{int}(^{n+1}{\boldsymbol {a}})-{^{n+1}{\boldsymbol {f}}}_{ext}}$ ${\displaystyle ={\boldsymbol {f}}_{int}({^{n}{\boldsymbol {a}}}+{^{n+1}\Delta }{\boldsymbol {a}})-{^{n}{\boldsymbol {f}}}_{ext}-{^{n+1}\Delta }{\boldsymbol {f}}_{ext}}$
(139)

There are many methods that can be applied to compute ${\textstyle ^{n+1}\Delta {\boldsymbol {a}}}$, but perhaps the most used are the secant matrix method and the tangent matrix method.

With the first method, the internal forces are described by the secant stiffness matrix at the previous equilibrium point:

 ${\displaystyle {\boldsymbol {f}}_{int}({^{n}{\boldsymbol {a}}}+{^{n+1}\Delta }{\boldsymbol {a}})\approx {\boldsymbol {K}}_{sec}(^{n}{\boldsymbol {a}})({^{n}{\boldsymbol {a}}}+{^{n+1}\Delta }{\boldsymbol {a}})}$
(140)

Then one just needs to solve the resulting linear system of equations:

 ${\displaystyle {\boldsymbol {K}}_{sec}(^{n}{\boldsymbol {a}})({^{n}{\boldsymbol {a}}}+{^{n+1}\Delta }{\boldsymbol {a}})={^{n}{\boldsymbol {f}}}_{ext}+{^{n+1}\Delta }{\boldsymbol {f}}_{ext}}$
(141)

or equivalently

 ${\displaystyle {\boldsymbol {K}}_{sec}(^{n}{\boldsymbol {a}})\,={^{n+1}{\boldsymbol {f}}}_{ext}}$
(142)

On the other hand, the tangent matrix method is based on a linearization of the problem. This means that the internal forces are represented through a Taylor series expanded to first order, neglecting the higher order terms:

 ${\displaystyle {\boldsymbol {f}}_{int}({\boldsymbol {a}}+\Delta {\boldsymbol {a}})={\boldsymbol {f}}_{int}({\boldsymbol {a}})+{\frac {\partial {\boldsymbol {f}}_{int}}{\partial {\boldsymbol {a}}}}\Delta {\boldsymbol {a}}+\dots }$
(143)

where the derivative of the internal forces with respect to the vector of nodal displacements is the global tangent stiffness matrix:

 ${\displaystyle {\boldsymbol {K}}_{tan}={\frac {\partial {\boldsymbol {f}}_{int}}{\partial {\boldsymbol {a}}}}}$
(144)

Substituting (143) in (139) we obtain:

 ${\displaystyle {\boldsymbol {0}}={\boldsymbol {f}}_{int}(^{n}{\boldsymbol {a}})+{\boldsymbol {K}}_{tan}(^{n}{\boldsymbol {a}})\,{\boldsymbol {a}}-{^{n}{\boldsymbol {f}}}_{ext}-{^{n+1}\Delta }{\boldsymbol {f}}_{ext}}$
(145)

and by applying (135), we finally get the linear system of equations to be solved:

 ${\displaystyle {\boldsymbol {K}}_{tan}(^{n}{\boldsymbol {a}})\,{\boldsymbol {a}}={^{n+1}\Delta }{\boldsymbol {f}}_{ext}}$
(146)

So far, we have presented a purely incremental (or prediction) method. It is not an iterative method, and for each increment of force, equilibrium is assumed. The approximations introduced induce small errors at each step, and so the accuracy of the technique strongly depends on the non-linearity of the problem and on the size of the load increments. As a result, progressive drift should be expected (see Figure 28).

 Figure 28: Drifting of the solution computed with a purely incremental method.

Thereby, in order to avoid the cited drifting, it is advisable to use an incremental-iterative (or prediction-correction) method. In essence, one must compute a prediction of the solution, and then correct it by means of iterations (see Figure 29). Each load level must be seen as a new non-linear problem that requires iterations to converge.

 Figure 29: Example of an incremental-iterative method. Drifting is avoided.

The most commonly used iterative methods are the Picard's method and the Newton-Raphson's method. The first one is based on the secant matrix method, whereas the latter relies on the tangent matrix method.

Picard's method is very simple. At the beginning of a step ${\textstyle n+1}$ we must compute the prediction like:

 ${\displaystyle {^{n+1}{\boldsymbol {a}}}^{0}=[{\boldsymbol {K}}_{sec}(^{n}{\boldsymbol {a}})]^{-1}\,[{^{n+1}}{\boldsymbol {f}}_{ext}]}$
(147)

After that, we must correct that prediction by iterating. At a generic iteration ${\textstyle i+1}$ we have:

 ${\displaystyle {^{n+1}{\boldsymbol {a}}}^{i+1}=[{\boldsymbol {K}}_{sec}(^{n+1}{\boldsymbol {a}}^{i})]^{-1}\,[{^{n+1}}{\boldsymbol {f}}_{ext}]}$
(148)

Note that the superscript on the right of a variable represents the iteration number, whereas the superscript on the left shows the load step.

Thereby, at the end of each iteration we need to check the convergence criterion. A good option is to check a ratio of displacements and a ratio of residuals:

 ${\displaystyle {\frac {\|{^{n+1}{\boldsymbol {a}}}^{i+1}-{^{n+1}{\boldsymbol {a}}}^{i}\|}{\|{^{n+1}{\boldsymbol {a}}}^{i+1}\|}}\leq tolerance_{a}\quad ;\quad {\frac {\|{\boldsymbol {R}}(^{n+1}{\boldsymbol {a}}^{i+1})\|}{\|{\boldsymbol {R}}(^{n+1}{\boldsymbol {a}}^{0})\|}}\leq tolerance_{R}}$
(149)

With Newton-Raphson's method, we start each step computing the prediction just like in Picard's method. In this case, though, it is calculated from an incremental update:

 ${\displaystyle ^{n+1}{\boldsymbol {a}}^{0}={^{n}{\boldsymbol {a}}}+[{\boldsymbol {K}}_{tan}(^{n}{\boldsymbol {a}})]^{-1}\,[{^{n+1}\Delta }{\boldsymbol {f}}_{ext}]}$
(150)

Since in Newton's method it is more convenient to work with increments of displacements, let us define the sought solution at step ${\textstyle n+1}$ and iteration ${\textstyle i+1}$ like:

 ${\displaystyle {^{n+1}{\boldsymbol {a}}}^{i+1}={^{n+1}{\boldsymbol {a}}}^{i}+{^{n+1}\delta }{\boldsymbol {a}}^{i+1}}$
(151)

The unknown ${\textstyle {^{n+1}\delta }{\boldsymbol {a}}^{i+1}}$ is obtained by linearizing the residual around ${\textstyle {^{n+1}{\boldsymbol {a}}}^{i}}$ and imposing equilibrium:

 ${\displaystyle {\boldsymbol {R}}({\boldsymbol {a}}+\delta {\boldsymbol {a}})={\boldsymbol {R}}({\boldsymbol {a}})+{\frac {\partial {\boldsymbol {R}}}{\partial {\boldsymbol {a}}}}\delta {\boldsymbol {a}}+\dots ={\boldsymbol {0}}}$
(152)

 ${\displaystyle {\frac {\partial {\boldsymbol {R}}}{\partial {\boldsymbol {a}}}}={\frac {\partial {\boldsymbol {f}}_{int}}{\partial {\boldsymbol {a}}}}={\boldsymbol {K}}_{tan}}$
 ${\displaystyle ^{n+1}\delta {\boldsymbol {a}}^{i+1}=-[{\boldsymbol {K}}_{tan}(^{n+1}{\boldsymbol {a}}^{i})]^{-1}{\boldsymbol {R}}(^{n+1}{\boldsymbol {a}}^{i})}$
 ${}^{}$