An isotropic damage model for geomaterials in the KRATOS framework



=Agraïments=

Els mesos d'investigació i elaboració d'aquest treball han sigut molt profitosos, pero ha estat necessari estudiar molt i aprendre rapid. Després de tota la feina realitzada, crec que s'han assolit en bona mesura els objectius plantejats i aixo no hagués estat possible si no fos per tota la gent que tenia al voltant i que m'ha estat ajudant d'alguna manera.

Per aixo voldria agrair primer als meus pares, Lluis i Marta, que m'han donat suport en tot moment i s'han preocupat que mengés sempre bé.  També dono les gracies a el meu germa Oriol per haver-me acompanyat en les hores de redacció d'aquest document  quan ell estava d'examens finals, i a la meva germana MªAngels i el seu marit Gabriel per haver-me animat a sortir de tant en tant i passar una bona estona.

Estic especialment agraït pel suport del meu tutor, l'Eugenio Oñate, que ha sigut la persona que em va introduir en el tema tractat a la monografia i que m'ha permes aprendre tot el que sé avui.

Vull agrair també la col<math display="inline">\cdot </math>laboració del Miquel Portabella, sense qui no haguéssim pogut implementar el metode de refinament de malla adaptatiu d'aquest treball.

Per últim, també em vull recordar dels companys de CIMNE que m'han ajudat amb els dubtes que anaven sorgint durant l'elaboració d'aquest treball: el Miguel Ángel Celigueta, el Salva Latorre, el Ferran Arrufat, el Pablo Becker, el Miquel Santasusana, el Jordi Cotela i el Josep Maria Carbonell.

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

''Keywords'': quasi-brittle materials, continuum damage mechanics, material non-linearity, non-local damage models, mesh-adaptivity

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

''Paraules clau'': materials quasi-fragils, mecanica del dany continu, no-linealitat material, models de dany no locals, adaptabilitat de malla

=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 <span id='citeF-1'></span>[[#cite-1|[1]]], continuum damage mechanics <span id='citeF-2'></span><span id='citeF-3'></span>[[#cite-2|[2,3]]], plasticity <span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-4|[4,5]]], coupled damage and plasticity <span id='citeF-6'></span><span id='citeF-7'></span>[[#cite-6|[6,7]]], micro-plane theory <span id='citeF-8'></span><span id='citeF-9'></span>[[#cite-8|[8,9]]], discrete models using a lattice approach <span id='citeF-10'></span><span id='citeF-11'></span>[[#cite-10|[10,11]]] and the discrete element method (DEM) <span id='citeF-12'></span><span id='citeF-13'></span>[[#cite-12|[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 <span id='citeF-14'></span>[[#cite-14|[14]]], but the concept of damage had already been introduced by Kachanov in 1958 in the context of creep rupture <span id='citeF-15'></span>[[#cite-15|[15]]], and further developed by Rabotnov <span id='citeF-16'></span>[[#cite-16|[16]]], Hayhurst <span id='citeF-17'></span>[[#cite-17|[17]]], and Leckie and Hayhurst <span id='citeF-18'></span>[[#cite-18|[18]]].

In the pioneering work of Kachanov <span id='citeF-15'></span>[[#cite-15|[15]]] the concept of ''effective stress'' was introduced, and by using continuum damage he solved problems related to creep in metals. Rabotnov <span id='citeF-16'></span>[[#cite-16|[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 <span id='citeF-19'></span>[[#cite-19|[19]]], and other important contributions to our knowledge about damage mechanics include: Mazars <span id='citeF-20'></span>[[#cite-20|[20]]], Mazars and Pijaudier-Cabot <span id='citeF-21'></span>[[#cite-21|[21]]], Chaboche <span id='citeF-22'></span>[[#cite-22|[22]]], Simo and Ju <span id='citeF-23'></span>[[#cite-23|[23]]], Ju <span id='citeF-24'></span><span id='citeF-25'></span>[[#cite-24|[24,25]]], Oliver, Oller and Oñate <span id='citeF-26'></span>[[#cite-26|[26]]], Oliver and Oller <span id='citeF-27'></span><span id='citeF-28'></span>[[#cite-27|[27,28]]], Cervera and Chiumenti <span id='citeF-29'></span><span id='citeF-30'></span>[[#cite-29|[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.

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-barra_fleixos.png|600px|Idealized material for the description of the uniaxial damage theory.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-1|1]]). Initially, all the fibers respond elastically, and the stress is carried by the total cross section of all fibers <math display="inline">S</math> (Figure [[#img-2|2]].1). As the applied strain is increased some fibers start breaking (Figure [[#img-2|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 <math display="inline">\bar{S}</math> (the area of unbroken fibers that can still carry stress) will decrease gradually from <math display="inline">\bar{S}=S</math> to <math display="inline">\bar{S}=0</math>. Of course, when the applied force diminishes (Figure [[#img-3|3]].2), the affected fibers remain broken and so the damaged area of the specimen is irrecoverable.

<div id='img-2'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-monotonic_load_1d.png|600px|Scheme of a uniaxial damage model through a monotonic loading process.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Scheme of a uniaxial damage model through a monotonic loading process.
|}

<div id='img-3'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-non_monotonic_load_1d.png|600px|Scheme of a uniaxial damage model through a non-monotonic loading process.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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'' <math display="inline">\sigma </math>, defined as the force per unit initial area of the cross section, and the ''effective stress'' <math display="inline">\bar{\sigma }</math>, 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, <math display="inline">\sigma S=\bar{\sigma }\bar{S}</math>, we obtain:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma =\frac{\bar{S}}{S}\bar{\sigma } </math>
|}
| style="width: 5px;text-align: right;" | (1)
|}

The ratio of the effective area to the total area, <math display="inline">\bar{S}/S</math>, 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:

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>d = 1-\frac{\bar{S}}{S}=\frac{S-\bar{S}}{S}=\frac{S_d}{S} </math>
|}
| style="width: 5px;text-align: right;" | (2)
|}

where <math display="inline">S_d=S-\bar{S}</math> is the damaged part of the area. An undamaged material is characterized then by <math display="inline">\bar{S}=S</math>, i.e., by <math display="inline">d=0</math>. 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 <math display="inline">d=1</math>, 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 <math display="inline">\bar{\sigma }</math> is related to the elastic strain of the material <math display="inline">\varepsilon </math> by the uniaxial Hooke's law:

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\bar{\sigma }=E\varepsilon  </math>
|}
| style="width: 5px;text-align: right;" | (3)
|}

where E is the elastic modulus of the undamaged material. Combining ([[#eq-1|1]]), ([[#eq-2|2]]) and ([[#eq-3|3]]), it is easily seen that the constitutive law for the nominal stress <math display="inline">\sigma </math> takes the form:

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma=(1-d)E\varepsilon  </math>
|}
| style="width: 5px;text-align: right;" | (4)
|}

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

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\begin{array}{cc}d=g(\varepsilon ) & \quad 0 \le d \le 1 \end{array} </math>
|}
| style="width: 5px;text-align: right;" | (5)
|}

Function <math display="inline">g</math> 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 [[#img-2|2]]. This description is valid only for monotonic loading by an increasing applied strain <math display="inline">\varepsilon </math>. When the material is first stretched up to a certain strain level <math display="inline">\varepsilon _1</math> that induces damage <math display="inline">d_1=g(\varepsilon _1)</math> and then the strain decreases (Figure [[#img-3|3]]), the damaged area remains constant and the material responds as an elastic material with a reduced Young's modulus <math display="inline">E_1=(1-d_1)E</math>. This means that, during unloading and partial reloading, the damage variable in ([[#eq-4|4]]) must be evaluated from the largest previously reached strain and not from the current strain <math display="inline">\varepsilon </math>. This means that it is crucial to introduce an internal state variable <math display="inline">r</math> characterizing the maximum strain level reached in the previous history of the material up to a given time <math display="inline">t</math>:

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>r(t)=\max \left\{r_0,\max _{\tau \le t}\varepsilon (\tau )\right\} </math>
|}
| style="width: 5px;text-align: right;" | (6)
|}

The expression implies that <math display="inline">r(t)\ge r_0</math>, where <math display="inline">r_0</math> is the so-called ''damage threshold'', a material parameter that represents the value of strain at which damage starts. In this formula, <math display="inline">t</math> is not necessarily the physical time (it can be any monotonically increasing parameter controlling the loading process). The damage evolution law ([[#eq-5|5]]) is then rewritten in the form:

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>d=g(r)\hbox{ with }\left\{\begin{array}{cc}g(r)=0 & \hbox{if }r=r_0\hbox{ }\\ 0<g(r)\le{1} & \hbox{if }r>r_0\hbox{ } \end{array}\right. </math>
|}
| style="width: 5px;text-align: 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 [[#img-3|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 ([[#eq-6|6]]) and ([[#eq-7|7]]), when the state variable <math display="inline">r</math> becomes greater than <math display="inline">r_0</math>, 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 <math display="inline">r</math>.

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

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\begin{array}{ccc}f\le{0}; & \dot{r}\ge{0}; & \dot{r}f=0 \end{array} </math>
|}
| style="width: 5px;text-align: right;" | (8)
|}

The first condition indicates that <math display="inline">r</math> can never be smaller than <math display="inline">\varepsilon </math>, while the second condition means that <math display="inline">r</math> cannot decrease. Finally, according to the third condition, <math display="inline">r</math> can grow only if the current values of <math display="inline">\varepsilon </math> and <math display="inline">r</math> 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:

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sigma =E_{sec}\varepsilon  </math>
|}
| style="width: 5px;text-align: right;" | (9)
|}

with:

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>E_{sec}=(1-d)E </math>
|}
| style="width: 5px;text-align: right;" | (10)
|}
* The law governing the evolution of the damage variable ([[#eq-7|7]])
* The damage criterion, comprising:

* The loading function:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>f(\varepsilon ,r)=\varepsilon{-}r </math>
|}
| style="width: 5px;text-align: right;" | (11)
|}

specifying the elastic domain <math display="inline">\mathcal{E}_r=\{\varepsilon |f(\varepsilon ,r)<0\}</math>, i.e., the set of states for which damage does not grow.
* The loading-unloading conditions ([[#eq-8|8]])

===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:

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{E}_{sec}=(1-d)\boldsymbol{E} </math>
|}
| style="width: 5px;text-align: right;" | (12)
|}

where <math display="inline">\boldsymbol{E}</math> is the elastic constitutive tensor of the intact material, and <math display="inline">d</math> is the damage variable. In the present context, <math display="inline">\boldsymbol{E}_{sec}</math> is the secant constitutive tensor that relates the total strain <math display="inline">\boldsymbol{\varepsilon }</math> to the total stress <math display="inline">\boldsymbol{\sigma }</math>, according to:

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }=\boldsymbol{E}_{sec}:\boldsymbol{\varepsilon }=(1-d)\boldsymbol{E}:\boldsymbol{\varepsilon } </math>
|}
| style="width: 5px;text-align: right;" | (13)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }=(1-d)\boldsymbol{\bar{\sigma }} </math>
|}
| style="width: 5px;text-align: right;" | (14)
|}

which is the multidimensional generalization of ([[#eq-1|1]]), and where <math display="inline">\boldsymbol{\bar{\sigma }}</math> is the ''effective stress tensor'' defined as:

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\bar{\sigma }}=\boldsymbol{E}:\boldsymbol{\varepsilon } </math>
|}
| style="width: 5px;text-align: right;" | (15)
|}

Similar to the uniaxial case, we introduce a loading function <math display="inline">f</math> specifying the elastic domain and the states at which damage grows. The loading function now depends on the strain tensor <math display="inline">\boldsymbol{\varepsilon }</math>, and on a variable <math display="inline">r</math> that controls the evolution of the elastic domain. Physically. <math display="inline">r</math> represents a scalar measure of the largest strain level ever reached in the history of the material. States for which <math display="inline">f(\boldsymbol{\varepsilon },r)<0</math> 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 <math display="inline">f(\boldsymbol{\varepsilon },r)=0</math>. Essentially, we can postulate the damage criterion for a multiaxial isotropic damage model with:

* The loading function:

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>f(\boldsymbol{\varepsilon },r)=\varepsilon _{eq}(\boldsymbol{\varepsilon })-r </math>
|}
| style="width: 5px;text-align: right;" | (16)
|}
* The loading-unloading conditions ([[#eq-8|8]])

where <math display="inline">\varepsilon _{eq}</math> is the ''equivalent strain'', i.e., a scalar measure of the strain level, and <math display="inline">r</math> is the largest value of equivalent strain calculated in the previous deformation history of the material up to its current state. In this regard, ([[#eq-6|6]]) can now be generalized to:

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>r(t)=\max \left\{r_0,\max _{\tau \le t}\varepsilon _{eq}(\tau )\right\} </math>
|}
| style="width: 5px;text-align: 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 ([[#eq-13|13]]). In essence, one must follow the scheme of Table [[#table-1|1]]:


{| class="wikitable" style="text-align: left; margin: 1em auto;"
|+ <span id='table-1'></span>Table. 1 Isotropic damage model scheme.

|-
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Input data for time <math>t+1</math>: <math>\boldsymbol{E}</math>, <math>\boldsymbol{\varepsilon }^t</math>, <math>\Delta \boldsymbol{\varepsilon }</math>, <math>r^t</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 1. Determine new strain vector: <math display="inline">\boldsymbol{\varepsilon }^{t+1}=\boldsymbol{\varepsilon }^t+\Delta \boldsymbol{\varepsilon }</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 2. Evaluate effective stresses: <math display="inline">\boldsymbol{\bar{\sigma }}^{t+1}=\boldsymbol{E}:\boldsymbol{\varepsilon }^{t+1}</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 3. Compute equivalent strain from the new strain vector: <math display="inline">\varepsilon _{eq}^{t+1}</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 4. Update <math display="inline">r</math> with ([[#eq-17|17]]): If <math display="inline">(\varepsilon _{eq}^{t+1}>r^t) \Rightarrow r^{t+1}=\varepsilon _{eq}^{t+1}</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 5. Update damage variable: <math display="inline">d^{t+1}=g(r^{t+1})</math>

|- style="border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | 6. Compute stresses: <math display="inline">\boldsymbol{\sigma }^{t+1}=(1-d^{t+1})\boldsymbol{\bar{\sigma }}^{t+1}</math>

|}

====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 [[#img-4|4]]).

<div id='img-4'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-elastic_domain.png|600px|3D elastic domain for a generic equivalent strain.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\varepsilon _{eq}=\| \boldsymbol{\varepsilon }\| =\sqrt{\boldsymbol{\varepsilon }:\boldsymbol{\varepsilon }}=\sqrt{\varepsilon _{ij}\varepsilon _{ij}} </math>
|}
| style="width: 5px;text-align: right;" | (18)
|}

or as the energy norm:

<span id="eq-19"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\varepsilon _{eq}=\sqrt{\frac{\boldsymbol{\varepsilon }:\boldsymbol{E}:\boldsymbol{\varepsilon }}{E}}=\sqrt{\frac{1}{E}E_{ijkl}\varepsilon _{ij}\varepsilon _{kl}} </math>
|}
| style="width: 5px;text-align: right;" | (19)
|}

where <math display="inline">E_{ijkl}</math> are the components of the elastic constitutive tensor <math display="inline">\boldsymbol{E}</math> and normalization by Young's modulus <math display="inline">E</math> is introduced to obtain a strain-like quantity.

The two norms of <math display="inline">\boldsymbol{\varepsilon }</math> 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 <span id='citeF-31'></span>[[#cite-31|[31]]]:

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\varepsilon _{eq}=\| \langle \boldsymbol{\varepsilon }\rangle \| =\sqrt{\langle \boldsymbol{\varepsilon }\rangle :\langle \boldsymbol{\varepsilon }\rangle } </math>
|}
| style="width: 5px;text-align: right;" | (20)
|}

or to its energetic counterpart:

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\varepsilon _{eq}=\sqrt{\frac{\langle \boldsymbol{\varepsilon }\rangle :\boldsymbol{E}:\langle \boldsymbol{\varepsilon }\rangle }{E}} </math>
|}
| style="width: 5px;text-align: right;" | (21)
|}

where McAuley brackets <math display="inline">\langle . \rangle </math> denote the "positive part" operator. For scalars <math display="inline">\langle x \rangle =\max (0,x)</math>, i.e.,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\langle x \rangle = \left\{ \begin{array}{ll}x & \hbox{if }x>0\hbox{ } \\ 0 & \hbox{if }x \le \hbox{ 0 } \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;" | (22)
|}

For symmetric tensors, such as the strain tensor <math display="inline">\boldsymbol{\varepsilon }</math>, the positive part is a tensor with the same principal directions <math display="inline">\boldsymbol{n_i}</math> as the original one, but replacing its principal values <math display="inline">\varepsilon _i</math> by their positive parts <math display="inline">\langle \varepsilon _i\rangle </math>. The subscript <math display="inline">i</math> ranges from 1 to 3 (the number of spatial dimensions). If we consider the following diagonalized strain tensor:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\varepsilon }= \begin{bmatrix}\varepsilon _1 & 0 & 0 \\ 0 & \varepsilon _2 & 0 \\ 0 & 0 & \varepsilon _3 \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;" | (23)
|}

the positive part of <math display="inline">\boldsymbol{\varepsilon }</math> is expressed as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (24)
|}

Thereby, knowing that <math display="inline">\langle \boldsymbol{\varepsilon }\rangle :\langle \boldsymbol{\varepsilon }\rangle =tr(\langle \boldsymbol{\varepsilon }\rangle ^{T}\langle \boldsymbol{\varepsilon }\rangle )</math>, definition ([[#eq-20|20]]) can be rewritten as:

<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\varepsilon _{eq}=\sqrt{\sum _{i=1}^{3} \langle \varepsilon _i\rangle ^{2}} </math>
|}
| style="width: 5px;text-align: right;" | (25)
|}

An alternative formula, called the modified von Mises definition <span id='citeF-32'></span>[[#cite-32|[32]]], reads:

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (26)
|}

where <math display="inline">\kappa </math> is a model parameter that sets the ratio between the uniaxial compressive strength and the uniaxial tensile strength, <math display="inline">\nu </math> is the Poisson's ratio of the material, <math display="inline">I_1</math> is the first invariant of the strain tensor, and <math display="inline">J_2</math> is the second invariant of the deviatoric strain tensor. Given a generic symmetric strain tensor <math display="inline">\boldsymbol{\varepsilon }</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\varepsilon }= \begin{bmatrix}\varepsilon _{xx} & \varepsilon _{xy} & \varepsilon _{xz} \\ \varepsilon _{xy} & \varepsilon _{yy} & \varepsilon _{yz} \\ \varepsilon _{xz} & \varepsilon _{yz} & \varepsilon _{zz} \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;" | (27)
|}

The first invariant <math display="inline">I_1</math> is the trace of the strain tensor:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>I_1=tr(\boldsymbol{\varepsilon }) = \varepsilon _{xx}+\varepsilon _{yy}+\varepsilon _{zz} </math>
|}
| style="width: 5px;text-align: right;" | (28)
|}

One can always decompose the strain tensor into its volumetric and deviatoric parts <math display="inline">\boldsymbol{\varepsilon }=\boldsymbol{\varepsilon }_v+\boldsymbol{\varepsilon }_d</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (29)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (30)
|}

From the deviatoric strain tensor <math display="inline">\boldsymbol{\varepsilon }_d</math> we can calculate <math display="inline">J_1</math> and <math display="inline">J_2</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>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 </math>
|}
| style="width: 5px;text-align: right;" | (31)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>J_2 </math>
| <math>  =  \frac{1}{2}(\boldsymbol{\varepsilon }_d:\boldsymbol{\varepsilon }_d-J_1^2)=\frac{1}{2}(\boldsymbol{\varepsilon }_d:\boldsymbol{\varepsilon }_d)  </math>
|-
| style="text-align: right;" | 
| <math> = \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]+ </math>
|-
| style="text-align: right;" | <math> \varepsilon _{xy}^2+\varepsilon _{xz}^2+\varepsilon _{yz}^2 </math>
|}
| style="width: 5px;text-align: right;" | (32)
|}

Another possibility is to work with an evolution of the model proposed by Simo and Ju <span id='citeF-23'></span>[[#cite-23|[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:

<span id="eq-33"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\varepsilon _{eq}=\left(\theta +\frac{1-\theta }{\kappa }\right)\sqrt{\boldsymbol{\varepsilon }:\boldsymbol{E}:\boldsymbol{\varepsilon }} </math>
|}
| style="width: 5px;text-align: right;" | (33)
|}

where the parameter <math display="inline">\theta </math> is a weighting factor depending on the effective stresses <math display="inline">\bar{\boldsymbol{\sigma }}</math>, given by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\theta =\frac{\sum _{i=1}^3\langle \bar{\sigma }_i\rangle }{\sum _{i=1}^3|\bar{\sigma }_i|} </math>
|}
| style="width: 5px;text-align: right;" | (34)
|}

As mentioned before, the parameter <math display="inline">\kappa </math> is defined by means of the ratio between the compression elastic limit <math display="inline">\sigma _{y}^c</math> and the tension elastic limit <math display="inline">\sigma _{y}^t</math>, i.e.:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\kappa =\frac{\sigma _{y}^c}{\sigma _{y}^t} </math>
|}
| style="width: 5px;text-align: right;" | (35)
|}

In the case of concrete <math display="inline">\kappa \approx{10}</math>.

The equivalent strains presented in this subsection lead to different damage surfaces but all three share a common trait, the elastic limit in tension <math display="inline">\sigma _{y}^t</math> is much lower than the elastic limit in compression <math display="inline">\sigma _{y}^c</math> (Figure [[#img-5|5]].b, [[#img-5|5]].c, [[#img-5|5]].d). In contrast, the equivalent strain obtained from the original energy norm definition ([[#eq-19|19]]) leads to a symmetric damage surface in which there is no difference between the elastic limit in tension and compression (Figure [[#img-5|5]].a).

<div id='img-5'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-damage_surfaces.png|600px|Damage surfaces for different equivalent strains in the 2D principal stress space.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 ([[#eq-7|7]]). There are various damage governing laws <math display="inline">g(r)</math> 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 <math display="inline">r_0</math> are the exponential law <span id='citeF-20'></span>[[#cite-20|[20]]]:

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>g(r)=1-\frac{r_0(1-A)}{r}-A\exp \{-B(r-r_0)\} </math>
|}
| style="width: 5px;text-align: right;" | (36)
|}

and the polynomial law <span id='citeF-33'></span>[[#cite-33|[33]]]:

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>g(r)=1-\frac{1}{1+B(r-r_0)+A(r-r_0)^2} </math>
|}
| style="width: 5px;text-align: right;" | (37)
|}

In ([[#eq-36|36]]) and ([[#eq-37|37]]) parameter <math display="inline">A</math> is associated to residual strength and parameter <math display="inline">B</math> 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 <math display="inline">r</math> between the damage threshold <math display="inline">r_0</math> and a maximum admissible value <math display="inline">r_{max}</math>, damage evolves according to:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>g(r)=\frac{r_{max}}{r_{max}-r_0}\left(1-\frac{r_0}{r}\right) </math>
|}
| style="width: 5px;text-align: 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 [[#img-6|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 <math display="inline">r_0=2\cdot{10}^{-4}</math> and a Young's modulus of <math display="inline">E=28000</math> MPa. The exponential law is represented with <math display="inline">A=0.8</math> and <math display="inline">B=15000</math>, the polynomial law with <math display="inline">A=0.8</math> and <math display="inline">B=30000</math>, and for the linear law we have defined <math display="inline">r_{max}=5\cdot{10}^{-4}</math>. One can see that when the strain exceeds the damage threshold after the elastic branch <math display="inline">\varepsilon{>}r_0</math>, a softening post-peak branch starts as a result of the damage growing.

<div id='img-6'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-softening_laws.png|600px|Unidimensional stress-strain curves for different damage evolution laws.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Unidimensional stress-strain curves for different damage evolution laws.
|}

An alternative exponential softening model was proposed in <span id='citeF-27'></span>[[#cite-27|[27]]]:

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>g(r)=1-\frac{r_0}{r}\exp \left\{A_f\left(1-\frac{r}{r_0}\right)\right\} </math>
|}
| style="width: 5px;text-align: right;" | (39)
|}

The parameter <math display="inline">A_f</math> is obtained from the following expression:

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>A_f=\left(\frac{G_f}{l_f(r_0)^2}-\frac{1}{2}\right)^{-1}\ge 0 </math>
|}
| style="width: 5px;text-align: right;" | (40)
|}

where <math display="inline">G_f</math> is the specific fracture energy per unit area, <math display="inline">l_f</math> is the characteristic length for the fractured domain, usually taken as the characteristic length of the finite elements, and <math display="inline">r_0</math> is the damage threshold which, for the Simo and Ju model, can be computed from:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>r_0=\frac{\sigma _{y}^t}{\sqrt{E}} </math>
|}
| style="width: 5px;text-align: right;" | (41)
|}

being <math display="inline">\sigma _{y}^t</math> the tensile strength of the material, and <math display="inline">E</math> the Young's modulus.

Another popular damage evolution law specifically designed for concrete was proposed by Mazars <span id='citeF-31'></span><span id='citeF-20'></span>[[#cite-31|[31,20]]]. He introduced two damage variables, <math display="inline">d_c</math> and <math display="inline">d_t</math>, that are computed from the same equivalent strain ([[#eq-25|25]]) using two different damage functions, <math display="inline">g_c</math> and <math display="inline">g_t</math>. Function <math display="inline">g_c</math> is identified from the uniaxial compressive test while <math display="inline">g_t</math> corresponds to the tensile test. The evolution law governing damage <math display="inline">g(r)</math> results from the combination of the two parts:

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>g(r)=\alpha _t^\beta g_t(r)+(1-\alpha _t)^\beta g_c(r) </math>
|}
| style="width: 5px;text-align: right;" | (42)
|}

Functions characterizing the evolution of damage in compression <math display="inline">g_c</math> and in tension <math display="inline">g_t</math> were proposed in the exponential form:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>g_c(r)=1-\frac{r_0(1-A_c)}{r}-A_c\exp \{-B_c(r-r_0)\} </math>
|}
| style="width: 5px;text-align: right;" | (43)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>g_t(r)=1-\frac{r_0(1-A_t)}{r}-A_t\exp \{-B_t(r-r_0)\} </math>
|}
| style="width: 5px;text-align: right;" | (44)
|}

where <math display="inline">A_c</math>, <math display="inline">B_c</math>, <math display="inline">A_t</math>, <math display="inline">B_t</math> are material parameters related to the shape of uniaxial stress-strain diagrams.

The coefficient <math display="inline">\alpha _t</math> in ([[#eq-42|42]]) ranges from 0 to 1 and takes into account the character of the stress state. It is evaluated from:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\alpha _t=\sum _{i=1}^3\frac{\varepsilon _{ti}\langle \varepsilon _i\rangle }{\varepsilon _{eq}^2} </math>
|}
| style="width: 5px;text-align: right;" | (45)
|}

where <math display="inline">\varepsilon _{ti}</math> are the principal strains due to positive effective stresses, i.e., the principal values of <math display="inline">\boldsymbol{\varepsilon }_t=\boldsymbol{E}^{-1}:\langle \boldsymbol{E}:\boldsymbol{\varepsilon }\rangle </math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\varepsilon }_t=-\frac{\nu }{E}tr(\langle \boldsymbol{\bar{\sigma }}\rangle )\boldsymbol{I}+\frac{1+\nu }{E}\langle \boldsymbol{\bar{\sigma }}\rangle  </math>
|}
| style="width: 5px;text-align: right;" | (46)
|}

being <math display="inline">\boldsymbol{I}</math> the identity tensor.

The parameter <math display="inline">\beta </math> in ([[#eq-42|42]]) was equal to 1 in the original version of the model <span id='citeF-31'></span>[[#cite-31|[31]]]. When it is higher than 1, <math display="inline">\beta </math> allows to slow down the evolution of damage under shear loading (when principal stresses have different sign).

The definition of <math display="inline">\varepsilon _{ti}</math> tells us that if all principal stresses are negative then <math display="inline">\alpha _t=0</math> and <math display="inline">d=d_c=g_c(r)</math>, and if all principal stresses are positive we have <math display="inline">\alpha _t=1</math> and <math display="inline">d=d_t=g_t(r)</math>. These are the "purely compressive" and "purely tensile" stress states. For intermediate stress states the value of <math display="inline">d</math> is between <math display="inline">d_c</math> and <math display="inline">d_t</math>, 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 <math display="inline">\boldsymbol{E}_{sec}</math> introduced in ([[#eq-12|12]]) and ([[#eq-13|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 (<math display="inline">f<0</math>). To obtain the tangent stiffness for a loading process with growing damage (<math display="inline">f=0</math> and <math display="inline">\dot{f}=0</math>), it is necessary to find the link between stress and strain increments. This concepts are clearly seen in a unidimensional case (see Figure [[#img-7|7]]). In this regard, knowing the difference between secant and tangent moduli (see Figure [[#img-8|8]].a and [[#img-8|8]].b) is useful to understand the generalization to secant and tangent constitutive tensors.

<div id='img-7'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-loading_unloading_damage.png|600px|Uniaxial stress-strain diagram of a loading-unloading process.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Uniaxial stress-strain diagram of a loading-unloading process.
|}

<div id='img-8'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-E_sec_tan.png|600px|Secant and tangent moduli.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Secant and tangent moduli.
|}

In essence, we define the elastic-damage tangent constitutive tensor <math display="inline">\boldsymbol{E}_{tan}</math> as the one that satisfies the following relation:

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\dot{\boldsymbol{\sigma }}=\boldsymbol{E}_{tan}:\dot{\boldsymbol{\varepsilon }} </math>
|}
| style="width: 5px;text-align: right;" | (47)
|}

Now, by considering the equation ([[#eq-13|13]]), we can obtain the rate of change of stress as follows:

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\dot{\boldsymbol{\sigma }}(\boldsymbol{\varepsilon },d) </math>
| <math> =\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} </math>
|-
| style="text-align: right;" | 
| <math> =(1-d)\boldsymbol{E}:\dot{\boldsymbol{\varepsilon }}-\bar{\boldsymbol{\sigma }}\otimes \dot{d} </math>
|}
| style="width: 5px;text-align: right;" | (48)
|}

from which we should distinguish two possible situations:

<ol>

<li>A process with elastic loading or unloading (<math display="inline">\dot{r}=0 \Rightarrow \dot{d}=0</math>) </li>

With <math display="inline">\dot{d}=0</math> equation ([[#eq-48|48]]) becomes <math display="inline">\dot{\boldsymbol{\sigma }}=(1-d)\boldsymbol{E}:\dot{\boldsymbol{\varepsilon }}</math>, from which we obtain that the elastic-damage tangent constitutive tensor coincides with the secant constitutive tensor:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{E}_{tan}=\boldsymbol{E}_{sec}=(1-d)\boldsymbol{E} </math>
|}
| style="width: 5px;text-align: right;" | (49)
|}
<li>A loading process with growing damage (<math display="inline">\dot{r}>0 \Rightarrow \dot{d}>0</math>) </li>

In this case, in order to obtain an explicit expression of the tangent constitutive tensor, the damage rate <math display="inline">\dot{d}</math> should be expressed in terms of the strain rate <math display="inline">\dot{\boldsymbol{\varepsilon }}</math>. This can be achieved by imposing the consistency condition <math display="inline">\dot{f}=0</math> in equation ([[#eq-16|16]]) and combining it with the rate of equation ([[#eq-7|7]]):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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 }} </math>
|}
| style="width: 5px;text-align: right;" | (50)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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 }} </math>
|}
| style="width: 5px;text-align: right;" | (51)
|}

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

<span id="eq-52"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{E}_{tan}=(1-d)\boldsymbol{E}-g'\bar{\boldsymbol{\sigma }}\otimes \boldsymbol{\eta } </math>
|}
| style="width: 5px;text-align: right;" | (52)
|}

where the form of <math display="inline">g'</math> and <math display="inline">\boldsymbol{\eta }</math> 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 <math display="inline">\varepsilon _{eq,1}</math> based on the energy norm of equation ([[#eq-19|19]]), tensor <math display="inline">\boldsymbol{\eta }_1</math> is given by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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}} </math>
|}
| style="width: 5px;text-align: right;" | (53)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{E}_{tan,1}=(1-d)\boldsymbol{E}-\frac{g'}{E\varepsilon _{eq,1}} \bar{\boldsymbol{\sigma }}\otimes \bar{\boldsymbol{\sigma }} </math>
|}
| style="width: 5px;text-align: right;" | (54)
|}

which preserves symmetry (<math display="inline">\boldsymbol{E}_{tan,1}^{T}=\boldsymbol{E}_{tan,1}</math>). It must be noticed though that, for other definitions of equivalent strain, this kind of symmetry is lost. 

</ol>

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


{| class="wikitable" style="text-align: left; margin: 1em auto;"
|+ <span id='table-2'></span>Table. 2 Flowchart for the numerical implementation of the elastic-damage tangent constitutive tensor.

|-
| colspan='2' style="text-align: center;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Input data for time <math>t+1</math>: <math>\boldsymbol{E}</math>, <math>\bar{\boldsymbol{\sigma }}^{t+1}</math>, <math>r^t</math>, <math>r^{t+1}</math>, <math>d^{t+1}</math>

|-
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | If (<math display="inline">r^{t+1}>r^t</math>) 
| style="border-left: 2px solid;border-right: 2px solid;" | Loading with growing damage (<math display="inline">\dot{r}>0 \Rightarrow \dot{d}>0</math>):

|-
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 
| style="border-left: 2px solid;border-right: 2px solid;" | 1. Compute damage function derivative: 
|-
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 
| colspan='1' style="text-align: center;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>g'=dg/dr</math>

|-
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 
| style="border-left: 2px solid;border-right: 2px solid;" | 2. Calculate equivalent strain partial derivative:
|-
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 
| colspan='1' style="text-align: center;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\boldsymbol{\eta }=\partial \varepsilon _{eq}/ \partial \boldsymbol{\varepsilon }</math>

|-
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 
| style="border-left: 2px solid;border-right: 2px solid;" | 3. Compute tangent constitutive tensor with ([[#eq-52|52]]):
|-
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 
| colspan='1' style="text-align: center;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\boldsymbol{E}_{tan}^{t+1}=(1-d^{t+1})\boldsymbol{E}-g'\bar{\boldsymbol{\sigma }}^{t+1}\otimes \boldsymbol{\eta }</math>

|-
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | Else 
| style="border-left: 2px solid;border-right: 2px solid;" | Elastic loading or unloading (<math display="inline">\dot{r}=0 \Rightarrow \dot{d}=0</math>):

|- style="border-bottom: 2px solid;"
| style="text-align: center;border-left: 2px solid;border-right: 2px solid;" | 
| colspan='1' style="text-align: center;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\boldsymbol{E}_{tan}^{t+1}=\boldsymbol{E}_{sec}^{t+1}=(1-d^{t+1})\boldsymbol{E}</math>

|}

==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 <math display="inline">d</math>. 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 <span id='citeF-34'></span><span id='citeF-35'></span>[[#cite-34|[34,35]]], strain gradient <span id='citeF-36'></span><span id='citeF-37'></span><span id='citeF-38'></span>[[#cite-36|[36,37,38]]], viscous <span id='citeF-39'></span><span id='citeF-40'></span><span id='citeF-41'></span>[[#cite-39|[39,40,41]]] and non-local terms <span id='citeF-42'></span><span id='citeF-43'></span><span id='citeF-44'></span><span id='citeF-45'></span>[[#cite-42|[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 <math display="inline">S</math> and a total length <math display="inline">L_0</math> under uniaxial tension (Figure [[#img-9|9]]). The material is assumed to obey a simple stress-strain law with linear elasticity up to the peak stress <math display="inline">\sigma _y</math>, followed by linear softening (Figure [[#img-10|10]]). If the bar is loaded in tension by an applied displacement <math display="inline">u</math> at its right extreme, the response remains linear elastic up to <math display="inline">u_y=L_0\varepsilon _y</math>, instant at which the force transmitted by the bar (reaction at the support) attains its maximum value <math display="inline">F_y=\sigma _yS</math>. After that, the resistance of the bar starts decreasing until the strain reaches <math display="inline">\varepsilon _f</math> and the transmitted stress completely disappears.

<div id='img-9'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-bar_uniaxial.png|600px|Bar under uniaxial tension.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 9:''' Bar under uniaxial tension.
|}

<div id='img-10'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-stress-strain-lin-soft.png|600px|Stress-strain diagram with linear softening.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-11|11]]). However, at a given stress level <math display="inline">\sigma _1</math> between zero and <math display="inline">\sigma _y</math>, there are actually two values of strain, <math display="inline">\varepsilon _{1,u}</math> and <math display="inline">\varepsilon _{1,s}</math>, that satisfy the constitutive equation (Figure [[#img-12|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 <math display="inline">\sigma _1=E\varepsilon _{1,u}</math>, whereas a damaged one will fall in the softening branch with <math display="inline">\sigma _1=(1-d_1)E\varepsilon _{1,s}</math>.

<div id='img-11'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-axial_bar_unidim.png|600px|Axial force acting along the bar.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Axial force acting along the bar.
|}

<div id='img-12'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-unload_soft_strains.png|600px|Possible strain values corresponding to the same stress level.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-13|13]]). In Figure [[#img-13|13]] we have denoted by <math display="inline">L_s</math> the cumulative length of the softening regions and by <math display="inline">L_u=L_0-L_s</math> the cumulative length of the unloading regions. When stress is completely relaxed, the strain in the unloading region is <math display="inline">\varepsilon _u=0</math> and the strain in the softening region is <math display="inline">\varepsilon _s=\varepsilon _f</math>. Thus the total elongation of the bar in this case is <math display="inline">u_f=L_u\varepsilon _u+L_s\varepsilon _s=L_s\varepsilon _f</math>. At this point, although <math display="inline">\varepsilon _f</math> is perfectly known from the linear softening law, <math display="inline">L_s</math> is totally undetermined. This means that the problem has infinite possible solutions with its corresponding post-peak branches of the load-displacement diagram (Figure [[#img-14|14]]). This fan of post-peak branches is bounded on one side by the solution with uniform softening (<math display="inline">u_f=L_0\varepsilon _f</math>) and on the other side by the solution with uniform unloading (<math display="inline">u_f=0</math>). 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.

<div id='img-13'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-softening_length.png|600px|Piecewise constant strain profile along the bar.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Piecewise constant strain profile along the bar.
|}

<div id='img-14'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-fan_force_displ.png|600px|Fan of possible post-peak branches of the load-displacement diagram.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 <math display="inline">n</math> 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 (<math display="inline">L_s=L_0/n</math>) 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 [[#img-15|15]], for <math display="inline">n=1</math> all the bar is damaged and the softening length is the total length of the bar <math display="inline">L_s=L_0</math>, whereas for <math display="inline">n>1</math> the softening region is more localized with strains becoming especially important at the damaged element. In the limit case of <math display="inline">n\rightarrow \infty </math> the softening branch would coincide with the elastic branch.

<div id='img-15'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-mesh_dependence_bar.png|600px|Mesh dependence of the numerical results. (a) Load-displacement diagrams for different number of elements. (b) Strain profiles for a prescribed imposed displacement.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 Examples|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 <span id='citeF-46'></span>[[#cite-46|[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 ([[#eq-39|39]]), which depends on the characteristic length of the fractured domain <math display="inline">l_f</math> included in ([[#eq-40|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 <math display="inline">l_{lim}</math>, the fracture length takes the value of the characteristic length of the element <math display="inline">l_e</math>, whereas for elements smaller than the limit length, the fracture length is taken as such limit length.

<span id="eq-55"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>l_f=\left\{\begin{array}{ll}l_e & \hbox{if }l_e>l_{lim}\hbox{ } \\ l_{lim} & \hbox{if }l_e\le l_{lim}\hbox{ } \end{array}\right. </math>
|}
| style="width: 5px;text-align: 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.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>l_e=\sqrt{\frac{4A_e}{\pi }} </math>
|}
| style="width: 5px;text-align: right;" | (56)
|}

where <math display="inline">A_e</math> 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 <span id='citeF-47'></span>[[#cite-47|[47]]] or smeared crack models <span id='citeF-27'></span>[[#cite-27|[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 <math display="inline">1-d</math> 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 <span id='citeF-48'></span><span id='citeF-49'></span>[[#cite-48|[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 <span id='citeF-43'></span>[[#cite-43|[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 <math display="inline">f(\boldsymbol{x})</math> be some local field in a domain <math display="inline">V</math>, the corresponding non-local field is defined as:

<span id="eq-57"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\tilde{f}(\boldsymbol{x})=\int _V \alpha (\boldsymbol{x},\boldsymbol{\xi })f(\boldsymbol{\xi })d\boldsymbol{\xi } </math>
|}
| style="width: 5px;text-align: right;" | (57)
|}

where <math display="inline">\alpha (\boldsymbol{x},\boldsymbol{\xi })</math> is the chosen non-local weighting function.

In an infinite, isotropic and homogeneous medium, the weighting function depends only on the distance <math display="inline">D=\| \boldsymbol{x}-\boldsymbol{\xi }\| </math> between the ''source point'' <math display="inline">\boldsymbol{\xi }</math>, and the ''receiver point'' <math display="inline">\boldsymbol{x}</math>. Thereby, we usually write <math display="inline">\alpha (\boldsymbol{x},\boldsymbol{\xi })=\alpha _0(\| \boldsymbol{x}-\boldsymbol{\xi }\| )</math>, where <math display="inline">\alpha _0(D)</math> is usually chosen as a non-negative function monotonically decreasing for <math display="inline">D\ge 0</math>.

One possible <math display="inline">\alpha _0(D)</math> is the Gauss distribution function:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\alpha _0(D)=\exp \left[-\left(\frac{2D}{l_c}\right)^2\right] </math>
|}
| style="width: 5px;text-align: right;" | (58)
|}

where the characteristic length <math display="inline">l_c</math> 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:

<span id="eq-59"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\alpha _0(D)=\left\{\begin{array}{ll}\exp \left[-\left(\frac{2D}{l_c}\right)^2\right]& \hbox{if }|D|\le R\hbox{ } \\ 0 & \hbox{if }|D|>R\hbox{ } \end{array} \right.  </math>
|}
| style="width: 5px;text-align: right;" | (59)
|}

where the interaction radius <math display="inline">R</math> is a parameter related to the characteristic length <math display="inline">l_c</math>. In the present work, we have considered <math display="inline">R=l_c</math>.

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\alpha _0(D)=\left\{\begin{array}{ll}\left(1-\frac{D^2}{R^2}\right)^2 & \hbox{if }|D|\le R\hbox{ } \\ 0 & \hbox{if }|D|>R\hbox{ } \end{array} \right.  </math>
|}
| style="width: 5px;text-align: right;" | (60)
|}

In essence, the interaction radius <math display="inline">R</math> corresponds to the smallest distance between points <math display="inline">\boldsymbol{x}</math> and <math display="inline">\boldsymbol{\xi }</math> at which the interaction weight <math display="inline">\alpha _0(\| \boldsymbol{x}-\boldsymbol{\xi }\| )</math> 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 <math display="inline">R</math>, centered at <math display="inline">\boldsymbol{x}</math>, is called the domain of influence of point <math display="inline">\boldsymbol{x}</math>. 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 [[#img-16|16]]).

<div id='img-16'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-nonlocal_averaging_boundary.png|600px|Averaging zone when the domain of influence protrudes through the boundary of a solid.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 ([[#eq-57|57]]) will be calculated as a weighted sum over the values at all the finite element integration points <math display="inline">\boldsymbol{\xi }</math> lying within the non-local interaction radius <math display="inline">R</math>.

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:

<span id="eq-61"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\begin{array}{cc}\displaystyle \int _V \alpha (\boldsymbol{x},\boldsymbol{\xi })d\boldsymbol{\xi }=1 & \forall \boldsymbol{x}\in V \end{array} </math>
|}
| style="width: 5px;text-align: right;" | (61)
|}

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

<span id="eq-62"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\alpha (\boldsymbol{x},\boldsymbol{\xi })=\frac{\alpha _0(\| \boldsymbol{x}-\boldsymbol{\xi }\| )}{\displaystyle \int _V\alpha _0(\| \boldsymbol{x}-\boldsymbol{\zeta }\| )d\boldsymbol{\zeta }} </math>
|}
| style="width: 5px;text-align: right;" | (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 <span id='citeF-50'></span>[[#cite-50|[50]]].

In the loading function ([[#eq-16|16]]), the local equivalent strain <math display="inline">\varepsilon _{eq}</math> is replaced by its weighted spatial average:

<span id="eq-63"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\tilde{\varepsilon }_{eq}(\boldsymbol{x})=\int _V\alpha (\boldsymbol{x},\boldsymbol{\xi })\varepsilon _{eq}(\boldsymbol{\xi })d\boldsymbol{\xi } </math>
|}
| style="width: 5px;text-align: right;" | (63)
|}

The internal state variable <math display="inline">r</math> is then the largest previously reached value of the non-local equivalent strain:

<span id="eq-64"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>r(t)=\max \left\{r_0,\max _{\tau \le t}\tilde{\varepsilon }_{eq}(\tau )\right\} </math>
|}
| style="width: 5px;text-align: right;" | (64)
|}

It is important to note that the damage variable is evaluated from the non-local equivalent strain <math display="inline">\tilde{\varepsilon }_{eq}</math>, whereas the strains <math display="inline">\boldsymbol{\varepsilon }</math> used in ([[#eq-13|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 [[#table-3|3]].


{| class="wikitable" style="text-align: left; margin: 1em auto;"
|+ <span id='table-3'></span>Table. 3 Non-local damage model scheme.

|-
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | Input data for time <math>t+1</math>: <math>\boldsymbol{E}</math>, <math>\boldsymbol{\varepsilon }^t</math>, <math>\Delta \boldsymbol{\varepsilon }</math>, <math>r^t</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 1. Determine new strain vector: <math display="inline">\boldsymbol{\varepsilon }^{t+1}=\boldsymbol{\varepsilon }^t+\Delta \boldsymbol{\varepsilon }</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 2. Evaluate effective stresses: <math display="inline">\boldsymbol{\bar{\sigma }}^{t+1}=\boldsymbol{E}:\boldsymbol{\varepsilon }^{t+1}</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 3. Compute equivalent strain from the new strain vector: <math display="inline">\varepsilon _{eq}^{t+1}</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 4. Calculate non-local equivalent strain from ([[#eq-63|63]]): <math display="inline">\tilde{\varepsilon }_{eq}^{t+1}</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 5. Update <math display="inline">r</math> with ([[#eq-64|64]]): If <math display="inline">(\tilde{\varepsilon }_{eq}^{t+1}>r^t) \Rightarrow r^{t+1}=\tilde{\varepsilon }_{eq}^{t+1}</math>

|-
| style="border-left: 2px solid;border-right: 2px solid;" | 6. Update damage variable: <math display="inline">d^{t+1}=g(r^{t+1})</math>

|- style="border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | 7. Compute stresses: <math display="inline">\boldsymbol{\sigma }^{t+1}=(1-d^{t+1})\boldsymbol{\bar{\sigma }}^{t+1}</math>

|}

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

<div id='img-17'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-kratos_1.png|600px|KRATOS' logo.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 <math display="inline">\boldsymbol{K}^{(e)}</math> and the force vector <math display="inline">\boldsymbol{f}^{(e)}</math> for each element.

* Afterwards, the global stiffness matrix <math display="inline">\boldsymbol{K}</math> and the nodal force vector <math display="inline">\boldsymbol{f}</math> are obtained through the assembly of the previous elemental stiffness matrix and force vector.

* The resulting system of equations <math display="inline">\boldsymbol{K}\boldsymbol{a}=\boldsymbol{f}</math> must be solved in order to obtain the nodal displacement vector <math display="inline">\boldsymbol{a}</math>, i.e., the unknown variable of the FEM problem.

* From the obtained nodal displacement vector <math display="inline">\boldsymbol{a}</math>, it is quite straightforward to calculate the strains <math display="inline">\boldsymbol{\varepsilon }</math> and stresses <math display="inline">\boldsymbol{\sigma }</math> 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 [[#img-18|18]].

<div id='img-18'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-fem_scheme.png|600px|Global flowchart of the analysis of a system through the Finite Element Method.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-19|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 [[#img-20|20]]).

<div id='img-19'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-plane_stress.png|600px|Example of a plane stress problem.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' Example of a plane stress problem.
|}

<div id='img-20'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-plane_strain.png|600px|Example of a plane strain problem.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 <math display="inline">z</math> axis suffer identical deformations contained in its own plane. This way, considering a generic cross section contained in the <math display="inline">x-y</math> plane of any of the structures in Figures [[#img-19|19]] and [[#img-20|20]], the displacement field of the structure can be perfectly defined if we know the displacement in the directions <math display="inline">x</math> and <math display="inline">y</math> for all points of that section. Thereby, the vector of displacements at any point is defined as:

<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{u}(x,y)=\left\{\begin{array}{c}u(x,y)\\ v(x,y) \end{array}\right\} </math>
|}
| style="width: 5px;text-align: right;" | (65)
|}

where <math display="inline">u(x,y)</math> and <math display="inline">v(x,y)</math> are the displacements of the point in the directions <math display="inline">x</math> and <math display="inline">y</math>, respectively.

From the displacements field ([[#eq-65|65]]) one can easily obtain the strains field following the general elasticity theory <span id='citeF-51'></span>[[#cite-51|[51]]]:

<span id="eq-66"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\varepsilon _x </math>
| <math> = \displaystyle \frac{\partial u}{\partial x} </math>
|-
| style="text-align: right;" | <math> \varepsilon _y </math>
| <math> = \displaystyle \frac{\partial v}{\partial y} </math>
| style="width: 5px;text-align: right;" | (66)
|-
| style="text-align: right;" | <math> \gamma _{xy} </math>
| <math> = \displaystyle \frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} </math>
|-
| style="text-align: right;" | <math> \gamma _{xz} </math>
| <math> = \gamma _{yz}=0  </math>
|}
|}

In plane strain problems, there is the additional hypothesis that the strain <math display="inline">\varepsilon _z</math> is null, whereas in plane stress problems it is the stress <math display="inline">\sigma _z</math> that is zero. This makes that, in any case, it is no necessary to consider strain <math display="inline">\varepsilon _z</math> because it does not take part in the work of deformations since <math display="inline">\sigma _z\varepsilon _z=0</math>. Thereby, at any point in a plane stress or plane strain state, we can define the vector of significant strains like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{ \varepsilon }=[\varepsilon _x,\varepsilon _y,\gamma _{xy}]^T </math>
|}
| style="width: 5px;text-align: right;" | (67)
|}

From equations in ([[#eq-66|66]]) we can deduce that tangential stresses <math display="inline">\tau _{xz}</math> and <math display="inline">\tau _{yz}</math> are zero. Moreover, for the same reasons of the previous paragraph, stress <math display="inline">\sigma _z</math> does not work, and so the significant stress vector can be defined as follows:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{ \sigma }=[\sigma _x,\sigma _y,\tau _{xy}]^T </math>
|}
| style="width: 5px;text-align: right;" | (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 <span id='citeF-51'></span>[[#cite-51|[51]]]. The elastic stress-strain relation can be written like:

<span id="eq-69"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }=\boldsymbol{E}\boldsymbol{\varepsilon } </math>
|}
| style="width: 5px;text-align: right;" | (69)
|}

where <math display="inline">\boldsymbol{E}</math> is the elastic constitutive matrix:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{E}= \begin{bmatrix}E_{11} & E_{12} & 0 \\ E_{21} & E_{22} & 0 \\ 0 & 0 & E_{33} \end{bmatrix} </math>
|}
| style="width: 5px;text-align: right;" | (70)
|}

For elastic isotropic materials the values in matrix <math display="inline">\boldsymbol{E}</math> are given by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | 
| <math>\hbox{Plane Stress }</math>
| style="text-align: right;" | <math> \qquad </math>
| <math>\hbox{Plane Strain } </math>
|-
| style="text-align: right;" | 
| <math>E_{11}=E_{22}=\frac{E}{1-\nu ^2}</math>
| style="text-align: right;" | 
| <math>E_{11}=E_{22}=\frac{E(1-\nu )}{(1+\nu )(1-2\nu )} </math>
|-
| style="text-align: right;" | 
| <math>E_{12}=E_{21}=E_{11}\nu </math>
| style="text-align: right;" | 
| <math>E_{12}=E_{21}=E_{11}\frac{\nu }{1-\nu } </math>
| style="width: 5px;text-align: right;" | (71)
|-
| style="text-align: right;" | 
| <math>E_{33}=\frac{E}{2(1+\nu )}=G</math>
| style="text-align: right;" | 
| <math>E_{33}=\frac{E}{2(1+\nu )}=G  </math>
|}
|}

with <math display="inline">E</math> being the Young's modulus, and <math display="inline">\nu </math> 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 <math display="inline">N</math> 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 <math display="inline">(\xi ,\eta )</math> are normalized coordinates, meaning that they range from -1 to +1, like it is shown in Figure [[#img-21|21]]. From this same figure it can be deduced that:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\xi =\frac{x-x_c}{a} \qquad ;\qquad \eta =\frac{y-y_c}{b} </math>
|}
| style="width: 5px;text-align: right;" | (72)
|}

where <math display="inline">x_c</math> and <math display="inline">y_c</math> are the coordinates of the center of the element. Moreover:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{d\xi }{dx}=\frac{1}{a} \qquad ;\qquad \frac{d\eta }{dy}=\frac{1}{b} </math>
|}
| style="width: 5px;text-align: right;" | (73)
|}

and so a differential of area is obtained like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>dxdy=abd\xi d\eta  </math>
|}
| style="width: 5px;text-align: right;" | (74)
|}

and, as it will appear later on, the integral of a function <math display="inline">f(x,y)</math> over a rectangular element can be perform by means of the following transformation to the natural coordinates:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\iint _{A^{(e)}}f(x,y)dxdy=\int _{-1}^{+1}\int _{-1}^{+1}g(\xi ,\eta )abd\xi d\eta  </math>
|}
| style="width: 5px;text-align: right;" | (75)
|}

<div id='img-21'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-natural_coord.png|600px|Geometry of a generic rectangular element. Natural and Cartesian coordinates.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 <math display="inline">n</math> can be written like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>f(\xi ,\eta )=\sum _{i=1}^p\alpha _i\xi ^j\eta ^k \qquad ;\qquad j+k\le n </math>
|}
| style="width: 5px;text-align: right;" | (76)
|}

where the number of terms of the polynomial is:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>p=\frac{(n+1)(n+2)}{2} </math>
|}
| style="width: 5px;text-align: right;" | (77)
|}

This way, a complete linear polynomial (<math display="inline">n=1 \rightarrow p=3</math>) would be:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>f(\xi ,\eta )=\alpha _1+\alpha _2\xi{+\alpha}_3\eta  </math>
|}
| style="width: 5px;text-align: right;" | (78)
|}

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

<div id='img-22'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-pascal.png|600px|Pascal's triangle in two dimensions.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>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\} </math>
|}
| style="width: 5px;text-align: right;" | (79)
|}

* ''Rigid body condition''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sum _{i=1}^nN_i(\xi ,\eta )=1 </math>
|}
| style="width: 5px;text-align: right;" | (80)
|}

where <math display="inline">i</math> and <math display="inline">j</math> represent node indexes, and <math display="inline">n</math> 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 <math display="inline">\xi </math> and <math display="inline">\eta </math> at that node. Thereby, if we consider <math display="inline">l_I^i(\xi )</math> like the Lagrange polynomial of order <math display="inline">I</math> in direction <math display="inline">\xi </math> at node <math display="inline">i</math>, and <math display="inline">l_J^i(\eta )</math> the one of order <math display="inline">J</math> in direction <math display="inline">\eta </math>, the shape function at that node is:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>N_i(\xi ,\eta )=l_I^i(\xi )l_J^i(\eta ) </math>
|}
| style="width: 5px;text-align: right;" | (81)
|}

In particular, for the 4-node rectangular Lagrangian element, the uni-dimensional Lagrange polynomials in each direction <math display="inline">\xi </math> and <math display="inline">\eta </math> coincide with the shape functions of the 2-node bar element. Thereby, at a node <math display="inline">i</math> of the element we have:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>l_1^i(\xi )=\frac{1}{2}(1+\xi \xi _i) \qquad ;\qquad l_1^i(\eta )=\frac{1}{2}(1+\eta \eta _i) </math>
|}
| style="width: 5px;text-align: right;" | (82)
|}

And consequently, the shape function at node <math display="inline">i</math> results:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>N_i(\xi ,\eta )=l_1^i(\xi )l_1^i(\eta )=\frac{1}{4}(1+\xi \xi _i)(1+\eta \eta _i) </math>
|}
| style="width: 5px;text-align: right;" | (83)
|}

The values of <math display="inline">\xi _i</math> and <math display="inline">\eta _i</math> can be seen in Figure [[#img-21|21]] but are summarized in Table [[#table-4|4]].


{| class="wikitable" style="text-align: center; margin: 1em auto;"
|+ <span id='table-4'></span>Table. 4 Natural coordinates for the rectangular Lagrangian element of four nodes.
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  '''Node''' 
| style="border-left: 2px solid;" | <math>\xi _i</math>
| style="border-right: 2px solid;" | <math>\eta _i</math>
|- style="border-top: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" |  1 
| style="border-left: 2px solid;" | -1 
| style="border-right: 2px solid;" | -1
|-
| style="border-left: 2px solid;border-right: 2px solid;" | 2 
| style="border-left: 2px solid;" | 1 
| style="border-right: 2px solid;" | -1
|-
| style="border-left: 2px solid;border-right: 2px solid;" | 3 
| style="border-left: 2px solid;" | 1 
| style="border-right: 2px solid;" | 1
|- style="border-bottom: 2px solid;"
| style="border-left: 2px solid;border-right: 2px solid;" | 4 
| style="border-left: 2px solid;" | -1 
| style="border-right: 2px solid;" | 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 [[#img-23|23]], we can formulate the shape functions for each of its three nodes like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>N_1=1-\xi{-\eta}\quad ;\quad N_2=\xi \quad ;\quad N_3=\eta  </math>
|}
| style="width: 5px;text-align: right;" | (84)
|}

<div id='img-23'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-triangle_nat_coord.png|600px|Natural coordinates in a triangular element.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-21|21]] and [[#img-23|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 <math display="inline">n</math> nodes, we can write the displacement field like:

<span id="eq-85"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>u=\displaystyle \sum _{i=1}^nN_iu_i \quad ;\quad v=\displaystyle \sum _{i=1}^nN_iv_i </math>
|}
| style="width: 5px;text-align: right;" | (85)
|}

where <math display="inline">(u_i,v_i)</math> and <math display="inline">N_i</math> are the horizontal and vertical displacements and the shape function of node <math display="inline">i</math>. Equations in ([[#eq-85|85]]) can be rewritten in matrix format like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (86)
|}

or equivalently:

<span id="eq-87"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{u}=\boldsymbol{N}\boldsymbol{a}^{(e)} </math>
|}
| style="width: 5px;text-align: right;" | (87)
|}

where <math display="inline">\boldsymbol{u}</math> is the vector of displacements at a point of the element,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (88)
|}

are the matrix of shape functions of the element, and the matrix of shape functions of node <math display="inline">i</math>, respectively, and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (89)
|}

are the vector of nodal displacements of the element, and the vector of nodal displacements of node <math display="inline">i</math>.

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\varepsilon _x </math>
| <math> =\frac{\partial u}{\partial x}=\sum _{i=1}^n\frac{\partial N_i}{\partial x}u_i </math>
|-
| style="text-align: right;" | <math> \varepsilon _y </math>
| <math> =\frac{\partial v}{\partial y}=\sum _{i=1}^n\frac{\partial N_i}{\partial y}v_i </math>
| style="width: 5px;text-align: right;" | (90)
|-
| style="text-align: right;" | <math> \gamma _{xy} </math>
| <math> =\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) </math>
|}
|}

and in matrix format:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (91)
|}

or

<span id="eq-92"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\varepsilon }=\boldsymbol{B}\boldsymbol{a}^{(e)} </math>
|}
| style="width: 5px;text-align: right;" | (92)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{B}=[\boldsymbol{B}_1,\boldsymbol{B}_2,\dots ,\boldsymbol{B}_n]  </math>
|}
| style="width: 5px;text-align: right;" | (93)
|}

is the deformation matrix of the element, and

<span id="eq-94"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (94)
|}

is the deformation matrix of node <math display="inline">i</math>.

Now, substituting equation ([[#eq-92|92]]) in the relation ([[#eq-69|69]]), we can express the stress vector like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }=\boldsymbol{E}\boldsymbol{\varepsilon }=\boldsymbol{E}\boldsymbol{B}\boldsymbol{a}^{(e)} </math>
|}
| style="width: 5px;text-align: right;" | (95)
|}

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

<span id="eq-96"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;" | (96)
|}

====3.2.1.4 Equilibrium equations of the discretization====

The expressions of the stiffness matrix <math display="inline">\boldsymbol{K}</math> and the force vector <math display="inline">\boldsymbol{f}</math> 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 <math display="inline">\boldsymbol{b}</math>), and uniformly distributed forces per unit length act over one of its sides (surface forces <math display="inline">\boldsymbol{t}</math>). 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 <math display="inline">\boldsymbol{q}</math>) 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:

<span id="eq-97"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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)} </math>
|}
| style="width: 5px;text-align: right;" | (97)
|}

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

From ([[#eq-87|87]]) and ([[#eq-92|92]]), we can write:

<span id="eq-98"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;" | (98)
|}

Substituting ([[#eq-98|98]]) in ([[#eq-97|97]]), and rearranging terms, we obtain:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>[\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)} </math>
|}
| style="width: 5px;text-align: right;" | (99)
|}

Since the virtual displacements <math display="inline">\delta \boldsymbol{a}^{(e)}</math> are arbitrary, we can write:

<span id="eq-100"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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)} </math>
|}
| style="width: 5px;text-align: right;" | (100)
|}

Equation ([[#eq-100|100]]) represents the equilibrium between the nodal forces of equilibrium <math display="inline">\boldsymbol{q}^{(e)}</math> and the forces from the deformation of the element (first integral), the mass forces <math display="inline">\boldsymbol{b}</math> and the surface forces <math display="inline">\boldsymbol{t}</math>. If we now substitute the stresses from expression ([[#eq-96|96]]), we get:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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)} </math>
|}
| style="width: 5px;text-align: right;" | (101)
|}

and rearranging terms:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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+ </math>
|-
| style="text-align: right;" | <math> +\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)} </math>
|}
| style="width: 5px;text-align: right;" | (102)
|}

which can be expressed like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}^{(e)}\boldsymbol{a}^{(e)}-\boldsymbol{f}^{(e)}=\boldsymbol{q}^{(e)} </math>
|}
| style="width: 5px;text-align: right;" | (103)
|}

where:

<span id="eq-104"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}^{(e)}=\iint _{A^{(e)}}\boldsymbol{B}^T\boldsymbol{E}\boldsymbol{B}\;tdA </math>
|}
| style="width: 5px;text-align: right;" | (104)
|}

is the elastic stiffness matrix of the element, and:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{f}^{(e)}=\boldsymbol{f}_\varepsilon ^{(e)}+\boldsymbol{f}_\sigma ^{(e)}+\boldsymbol{f}_b^{(e)}+\boldsymbol{f}_t^{(e)} </math>
|}
| style="width: 5px;text-align: right;" | (105)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;" | (106)
|}

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

<span id="eq-107"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{f}_b^{(e)}=\iint _{A^{(e)}}\boldsymbol{N}^T\boldsymbol{b}\;tdA  </math>
|}
| style="width: 5px;text-align: right;" | (107)
|}

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

<span id="eq-108"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{f}_t^{(e)}=\oint _{l^{(e)}}\boldsymbol{N}^T\boldsymbol{t}\;tds </math>
|}
| style="width: 5px;text-align: right;" | (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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\sum _e\boldsymbol{q}_i^{(e)}=\boldsymbol{p}_j </math>
|}
| style="width: 5px;text-align: right;" | (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 <math display="inline">j</math>, and <math display="inline">\boldsymbol{p}_j</math> 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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}\boldsymbol{a}=\boldsymbol{f} </math>
|}
| style="width: 5px;text-align: right;" | (110)
|}

where <math display="inline">\boldsymbol{K}</math>, <math display="inline">\boldsymbol{a}</math> and <math display="inline">\boldsymbol{f}</math> 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 ([[#eq-85|85]]), we can express the geometry of a two-dimensional isoparametric element from the coordinates <math display="inline">x_i</math> and <math display="inline">y_i</math> of its nodes like:

<span id="eq-111"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>x=\sum _{i=1}^nN_i(\xi ,\eta )x_i \quad ;\quad y=\sum _{i=1}^nN_i(\xi ,\eta )y_i </math>
|}
| style="width: 5px;text-align: right;" | (111)
|}

where <math display="inline">N_i(\xi ,\eta )</math> are the element shape functions. Equations in ([[#eq-111|111]]) relate the cartesian coordinates of a point with the natural coordinates <math display="inline">\xi </math> and <math display="inline">\eta </math>, 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 <math display="inline">\xi </math> and <math display="inline">\eta </math>. If we consider <math display="inline">N_i(\xi ,\eta )=N_i(x(\xi ,\eta ),y(\xi ,\eta ))</math> then, applying the chain rule, we have:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (112)
|}

or in matrix format:

<span id="eq-113"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (113)
|}

where <math display="inline">\boldsymbol{J}^{(e)}</math> is the Jacobian matrix of the transformation of natural coordinates to cartesian coordinates. From ([[#eq-113|113]]) we can derive:

<span id="eq-114"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (114)
|}

where <math display="inline">\left|\boldsymbol{J}^{(e)}\right|</math> is the determinant of the Jacobian.

The determinant of the Jacobian permits expressing the differential of area in natural coordinates <span id='citeF-52'></span>[[#cite-52|[52]]] like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>dx\;dy=\left|\boldsymbol{J}^{(e)}\right|d\xi \;d\eta  </math>
|}
| style="width: 5px;text-align: right;" | (115)
|}

Now, using the isoparametric transformation ([[#eq-111|111]]) we can obtain the terms of the Jacobian:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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. } </math>
|}
|}

and so:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (116)
|}

Furthermore, substituting the expressions of equation ([[#eq-114|114]]) in ([[#eq-94|94]]) we obtain the deformation matrix of an isoparametric element in terms of the natural coordinates:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (117)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>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 } </math>
|}
| style="width: 5px;text-align: right;" | (118)
|}

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

For quadrilateral elements we have:

<span id="eq-119"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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  </math>
|}
| style="width: 5px;text-align: right;" | (119)
|}

And for triangular elements we must consider the different domain of integration (see Figure [[#img-23|23]]):

<span id="eq-120"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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  </math>
|}
| style="width: 5px;text-align: right;" | (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 <math display="inline">g(\xi ,\eta )</math> over the domain of natural coordinates of a quadrilateral element can be evaluated by means of a Gauss quadrature like:

<span id="eq-121"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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_pW_q </math>
|}
| style="width: 5px;text-align: right;" | (121)
|}

where <math display="inline">n_p</math> and <math display="inline">n_q</math> are the number of integration points in each direction <math display="inline">\xi </math> and <math display="inline">\eta </math>; <math display="inline">\xi _p</math> and <math display="inline">\eta _q</math> are the natural coordinates of the integration points <math display="inline">p</math>, <math display="inline">q</math> and <math display="inline">W_p</math>, <math display="inline">W_q</math> the weights of that point corresponding to each direction.

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

<span id="eq-122"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;" | (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 <math display="inline">n</math> in each natural direction permits integrating with no error a polynomial of order <math display="inline">\le 2n-1</math> 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 ([[#eq-121|121]]) in ([[#eq-119|119]]), we can obtain the expression of the elastic stiffness matrix for an isoparametric quadrilateral element, evaluated through numerical integration:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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_pW_q </math>
|}
| style="width: 5px;text-align: right;" | (123)
|}

And from ([[#eq-120|120]]) and ([[#eq-122|122]]) we obtain the expression corresponding to a triangular element:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;" | (124)
|}

Thereby, we can see that the numerical integration of the stiffness matrix requires the evaluation of the Jacobian <math display="inline">\boldsymbol{J}^{(e)}</math>, its determinant <math display="inline">\left|\boldsymbol{J}^{(e)}\right|</math>, the deformation matrix <math display="inline">\boldsymbol{B}</math>, the constitutive matrix <math display="inline">\boldsymbol{E}</math>, and the thickness <math display="inline">t</math> 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 ([[#eq-107|107]]), can be performed like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\boldsymbol{f}_b^{(e)} </math>
| <math> =\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 </math>
|-
| style="text-align: right;" | 
| <math> =\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_pW_q </math>
|}
| style="width: 5px;text-align: right;" | (125)
|}

For triangular elements the double sum is replaced by the simple sum of ([[#eq-122|122]]).

Finally, the computation of the vector of surface forces ([[#eq-108|108]]) is a bit different because the integral is performed over the contour of the element <math display="inline">l^{(e)}</math>. In general, such contour represents a straight line <math display="inline">\xi =</math>const or <math display="inline">\eta =</math>const in the space of natural coordinates. Therefore, for example, for a contour of an isoparametric quadrilateral element corresponding to the straight line <math display="inline">\eta=1</math>, the differential of length <math display="inline">ds</math> is computed like:

<span id="eq-126"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>[ds]_{\eta=1} </math>
| <math> =\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 </math>
|-
| style="text-align: right;" | 
| <math> =\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  </math>
|}
| style="width: 5px;text-align: right;" | (126)
|}

Substituting ([[#eq-126|126]]) in ([[#eq-108|108]]) we obtain a line integral that depends only on the natural coordinate <math display="inline">\xi </math> and that can be computed through a uni-dimensional quadrature like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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 </math>
|}
| style="width: 5px;text-align: right;" | (127)
|}

A more extended overview on the Finite Element Method can be found in <span id='citeF-53'></span><span id='citeF-54'></span>[[#cite-53|[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 ([[#eq-97|97]]), we used an elastic stress-strain relation ([[#eq-69|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 ([[#eq-13|13]]) in voigt notation:

<span id="eq-128"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\sigma }=\boldsymbol{E}_{sec}\boldsymbol{\varepsilon }=(1-d)\boldsymbol{E}\boldsymbol{\varepsilon } </math>
|}
| style="width: 5px;text-align: right;" | (128)
|}

and we introduce it into a general expression of the PVW ([[#eq-100|100]]), we obtain the following relation:

<span id="eq-129"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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)} </math>
|}
| style="width: 5px;text-align: right;" | (129)
|}

where <math display="inline">d</math> is the damage parameter.

Expression ([[#eq-129|129]]) can be rewritten like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}^{(e)}(d)\,\boldsymbol{a}^{(e)}-\boldsymbol{f}^{(e)}=\boldsymbol{q}^{(e)} </math>
|}
| style="width: 5px;text-align: right;" | (130)
|}

where

<span id="eq-131"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}^{(e)}(d)=\iint _{A^{(e)}}\boldsymbol{B}^T(1-d)\boldsymbol{E}\boldsymbol{B}tdA </math>
|}
| style="width: 5px;text-align: right;" | (131)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{f}^{(e)}=\boldsymbol{f}_b^{(e)}+\boldsymbol{f}_t^{(e)} </math>
|}
| style="width: 5px;text-align: right;" | (132)
|}

is the vector of equivalent nodal forces of the element, in which <math display="inline">\boldsymbol{f}_b^{(e)}</math> and <math display="inline">\boldsymbol{f}_t^{(e)}</math> coincide with ([[#eq-107|107]]) and ([[#eq-108|108]]), respectively.

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

<span id="eq-133"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}(d)\,\boldsymbol{a}=\boldsymbol{f} </math>
|}
| style="width: 5px;text-align: right;" | (133)
|}

where <math display="inline">\boldsymbol{K}(d)</math> is the global secant stiffness matrix <math display="inline">\boldsymbol{K}_{sec}</math>.

By the way, we saw in Chapter [[#2 Damage Mechanics|2]] that <math display="inline">d=g(r)</math>, where <math display="inline">r</math> is the internal state variable depending on the equivalent strains <math display="inline">\varepsilon _{eq}</math>. Therefore, since the strains are directly related to the displacements through ([[#eq-92|92]]), we can say that damage also depends on the displacements, i.e. <math display="inline">d=f(\boldsymbol{a})</math>. Taking into account this last consideration, expression ([[#eq-133|133]]) results:

<span id="eq-134"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{sec}(\boldsymbol{a})\,\boldsymbol{a}=\boldsymbol{f} </math>
|}
| style="width: 5px;text-align: right;" | (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 [[#img-24|24]] shows two possible non-linear response diagrams.

<div id='img-24'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-action-response.png|600px|Response diagrams: (a)Load-deflection diagram showing equilibrium path. (b)Diagram distinguishing primary from secondary equilibrium path.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-24|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 [[#img-24|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 [[#img-25|25]] illustrates three monotonic types of response: linear, hardening, and softening. Symbols <math display="inline">R</math>, <math display="inline">F</math> and <math display="inline">L</math> identify reference, failure and limit points, respectively.

<div id='img-25'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-non-linear_responses.png|600px|Basic behaviours for a non-linear response. (a)Linear until brittle failure. (b)Hardening. (c)Softening.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 25:''' Basic behaviours for a non-linear response. (a)Linear until brittle failure. (b)Hardening. (c)Softening.
|}

The response shown in Figure [[#img-25|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 [[#img-25|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 [[#img-25|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 [[#img-26|26]] show a combination of basic behaviours that can complicate the response of the structure. Here <math display="inline">B</math> and <math display="inline">T</math> denote bifurcation and turning points.

<div id='img-26'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-non-linear_responses_complex.png|600px|Complex non-linear responses. (a)Snap-through. (b)Snap-back. (c)Bifurcation. (d)Bifurcation with snap-back.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-26|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 [[#img-26|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 [[#img-26|26]].b and [[#img-26|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:

<ol>

<li>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 [[#img-27|27]].a)  </li>

<li>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 [[#img-27|27]].b) </li>

<li>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 [[#img-27|27]].c) </li>

</ol>

Figure [[#img-27|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.

<div id='img-27'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-incremental_solution.png|600px|Geometric representations of different strategies for advance control: (a)Force control. (b)Displacement control. (c)Arc-length control]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 ([[#eq-134|134]]) can be expressed like the residual between the internal and the external forces. Thereby, at step <math display="inline">n</math> of a force control strategy we can write:

<span id="eq-135"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{R}(^n\boldsymbol{a})=\boldsymbol{f}_{int}(^n\boldsymbol{a})-{^n\boldsymbol{f}}_{ext}=\boldsymbol{0} </math>
|}
| style="width: 5px;text-align: right;" | (135)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{f}_{int}(^n\boldsymbol{a})=\boldsymbol{K}_{sec}(^n\boldsymbol{a})\,^n\boldsymbol{a} </math>
|}
| style="width: 5px;text-align: right;" | (136)
|}

and <math display="inline">^n\boldsymbol{f}_{ext}</math> is the vector of external forces at the load step <math display="inline">n</math>, i.e. <math display="inline">^n\boldsymbol{f}_{ext}\equiv{^n\boldsymbol{f}}</math>.

Supposing that step <math display="inline">n</math> is in equilibrium, we must first increase the external load by a prescribed <math display="inline">{^{n+1}\Delta }\boldsymbol{f}_{ext}</math>. Consequently, the new external force is:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\boldsymbol{f}_{ext}={^n\boldsymbol{f}}_{ext}+{^{n+1}\Delta }\boldsymbol{f}_{ext} </math>
|}
| style="width: 5px;text-align: right;" | (137)
|}

The unknown is the incremental displacement <math display="inline">^{n+1}\Delta \boldsymbol{a}</math> that results from the previous force increment. Thereby:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\boldsymbol{a}={^n\boldsymbol{a}}+{^{n+1}\Delta }\boldsymbol{a} </math>
|}
| style="width: 5px;text-align: right;" | (138)
|}

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

<span id="eq-139"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\boldsymbol{0} </math>
| <math> =\boldsymbol{R}(^{n+1}\boldsymbol{a})=\boldsymbol{f}_{int}(^{n+1}\boldsymbol{a})-{^{n+1}\boldsymbol{f}}_{ext} </math>
|-
| style="text-align: right;" | 
| <math> =\boldsymbol{f}_{int}({^n\boldsymbol{a}}+{^{n+1}\Delta }\boldsymbol{a})-{^n\boldsymbol{f}}_{ext}-{^{n+1}\Delta }\boldsymbol{f}_{ext} </math>
|}
| style="width: 5px;text-align: right;" | (139)
|}

There are many methods that can be applied to compute <math display="inline">^{n+1}\Delta \boldsymbol{a}</math>, 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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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}) </math>
|}
| style="width: 5px;text-align: right;" | (140)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{sec}(^n\boldsymbol{a})({^n\boldsymbol{a}}+{^{n+1}\Delta }\boldsymbol{a})={^n\boldsymbol{f}}_{ext}+{^{n+1}\Delta }\boldsymbol{f}_{ext} </math>
|}
| style="width: 5px;text-align: right;" | (141)
|}

or equivalently

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{sec}(^n\boldsymbol{a})\,={^{n+1}\boldsymbol{f}}_{ext} </math>
|}
| style="width: 5px;text-align: right;" | (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:

<span id="eq-143"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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  </math>
|}
| style="width: 5px;text-align: right;" | (143)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{tan}=\frac{\partial \boldsymbol{f}_{int}}{\partial \boldsymbol{a}} </math>
|}
| style="width: 5px;text-align: right;" | (144)
|}

Substituting ([[#eq-143|143]]) in ([[#eq-139|139]]) we obtain:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;" | (145)
|}

and by applying ([[#eq-135|135]]), we finally get the linear system of equations to be solved:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{tan}(^n\boldsymbol{a})\,\boldsymbol{a}={^{n+1}\Delta }\boldsymbol{f}_{ext} </math>
|}
| style="width: 5px;text-align: right;" | (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 [[#img-28|28]]).

<div id='img-28'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-drifting.png|600px|Drifting of the solution computed with a purely incremental method.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-29|29]]). Each load level must be seen as a new non-linear problem that requires iterations to converge.

<div id='img-29'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-incremental-iterative.png|600px|Example of an incremental-iterative method. Drifting is avoided.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 <math display="inline">n+1</math> we must compute the prediction like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\boldsymbol{a}}^{0}=[\boldsymbol{K}_{sec}(^n\boldsymbol{a})]^{-1}\,[{^{n+1}}\boldsymbol{f}_{ext}] </math>
|}
| style="width: 5px;text-align: right;" | (147)
|}

After that, we must correct that prediction by iterating. At a generic iteration <math display="inline">i+1</math> we have:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\boldsymbol{a}}^{i+1}=[\boldsymbol{K}_{sec}(^{n+1}\boldsymbol{a}^i)]^{-1}\,[{^{n+1}}\boldsymbol{f}_{ext}] </math>
|}
| style="width: 5px;text-align: right;" | (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:

<span id="eq-149"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\|{^{n+1}\boldsymbol{a}}^{i+1}-{^{n+1}\boldsymbol{a}}^{i}\| }{\|{^{n+1}\boldsymbol{a}}^{i+1}\| } \le tolerance_a \quad ;\quad \frac{\| \boldsymbol{R}(^{n+1}\boldsymbol{a}^{i+1})\| }{\| \boldsymbol{R}(^{n+1}\boldsymbol{a}^{0})\| } \le tolerance_R </math>
|}
| style="width: 5px;text-align: right;" | (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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\boldsymbol{a}^0={^n\boldsymbol{a}}+[\boldsymbol{K}_{tan}(^n\boldsymbol{a})]^{-1}\,[{^{n+1}\Delta }\boldsymbol{f}_{ext}] </math>
|}
| style="width: 5px;text-align: right;" | (150)
|}

Since in Newton's method it is more convenient to work with increments of displacements, let us define the sought solution at step <math display="inline">n+1</math> and iteration <math display="inline">i+1</math> like:

<span id="eq-151"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\boldsymbol{a}}^{i+1}={^{n+1}\boldsymbol{a}}^i+{^{n+1}\delta }\boldsymbol{a}^{i+1} </math>
|}
| style="width: 5px;text-align: right;" | (151)
|}

The unknown <math display="inline">{^{n+1}\delta }\boldsymbol{a}^{i+1}</math> is obtained by linearizing the residual around <math display="inline">{^{n+1}\boldsymbol{a}}^i</math> and imposing equilibrium:

<span id="eq-152"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{R}(\boldsymbol{a}+\delta \boldsymbol{a})=\boldsymbol{R}(\boldsymbol{a})+\frac{\partial \boldsymbol{R}}{\partial \boldsymbol{a}}\delta \boldsymbol{a}+\dots =\boldsymbol{0} </math>
|}
| style="width: 5px;text-align: right;" | (152)
|}

Considering conservative loads we have:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\partial \boldsymbol{R}}{\partial \boldsymbol{a}}=\frac{\partial \boldsymbol{f}_{int}}{\partial \boldsymbol{a}}=\boldsymbol{K}_{tan} </math>
|}
|}

and so from the relation in ([[#eq-152|152]]) we can write:

<span id="eq-153"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\delta \boldsymbol{a}^{i+1}=-[\boldsymbol{K}_{tan}(^{n+1}\boldsymbol{a}^i)]^{-1}\boldsymbol{R}(^{n+1}\boldsymbol{a}^i) </math>
|}
| style="width: 5px;text-align: right;" | (153)
|}

Substituting ([[#eq-153|153]]) in ([[#eq-151|151]]) we get the final expression for the sought solution at any iteration:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\boldsymbol{a}}^{i+1}={^{n+1}\boldsymbol{a}}^i-[\boldsymbol{K}_{tan}(^{n+1}\boldsymbol{a}^i)]^{-1}\boldsymbol{R}(^{n+1}\boldsymbol{a}^i) </math>
|}
| style="width: 5px;text-align: right;" | (154)
|}

Finally, at the end of each iteration we just need to check the convergence criterion with ([[#eq-149|149]]).

Comparing Picard's and Newton-Raphson's methods, we must say that the first is cheap and robust, but it shows linear convergence rate. The latter, on the other hand, is much more expensive and less robust, but shows quadratic convergence.

To finish this section, we will present the arc-length strategy implemented in this work.

====3.3.2.2 Modified Riks-Wempner arc-length method====

The arc-length method is a solution strategy in which the path through a converged solution, at any step, follows a direction orthogonal to the tangent of the solution curve. In this procedure, both the load vector and the displacement field vary.

Thereby, the main differential aspect of the arc-length method with respect to other incremental-iterative techniques is that the former introduces an additional scalar variable in the residual equation, in order to account for the variation of the load vector. In this regard, we should rewrite expression ([[#eq-135|135]]) like:

<span id="eq-155"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{R}(\boldsymbol{a},\lambda )=\boldsymbol{f}_{int}(\boldsymbol{a})-\boldsymbol{f}_{ext}(\lambda )=\boldsymbol{K}_{sec}(\boldsymbol{a})\boldsymbol{a}-\lambda \boldsymbol{f}=\boldsymbol{0} </math>
|}
| style="width: 5px;text-align: right;" | (155)
|}

where <math display="inline">\lambda </math> is the load-level parameter that will be considered as a new unknown of the problem.

As a result, we have <math display="inline">N+1</math> degrees of freedom from the vector of nodal displacements <math display="inline">\boldsymbol{a}</math> and the load-level parameter <math display="inline">\lambda </math>, but we only have <math display="inline">N</math> equations in ([[#eq-155|155]]). Therefore, we need an additional constraint equation to complete the definition of the problem:

<span id="eq-156"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>h(\boldsymbol{a},\lambda )=0 </math>
|}
| style="width: 5px;text-align: right;" | (156)
|}

The form of the generic function <math display="inline">h((\boldsymbol{a},\lambda )</math> will depend on the chosen arc-length strategy. In our case, we chose the modified Riks-Wempner method (or Ramm's method).

In <span id='citeF-55'></span>[[#cite-55|[55]]] Ramm proposed the following relation for ([[#eq-156|156]]):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>h(\boldsymbol{a},\lambda )=[{^{n+1}\Delta }\boldsymbol{a}]^T\,[{^{n+1}\Delta }\boldsymbol{a}]+({^{n+1}\Delta }\lambda )^2-(\Delta S)^2=0 </math>
|}
| style="width: 5px;text-align: right;" | (157)
|}

where <math display="inline">\Delta S</math> is the arc-length increment, defined as the distance along the response curve <math display="inline">(\boldsymbol{a},\lambda )</math>; and <math display="inline">{^{n+1}\Delta }\boldsymbol{a}</math> and <math display="inline">{^{n+1}\Delta }\lambda </math> are, respectively, the increments in the displacement field and load parameter from step <math display="inline">n</math> to step <math display="inline">n+1</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\Delta }\boldsymbol{a}={^{n+1}\boldsymbol{a}}-{^{n}\boldsymbol{a}} \quad ;\quad {^{n+1}\Delta }\lambda ={^{n+1}\lambda }-{^{n}\lambda } </math>
|}
| style="width: 5px;text-align: right;" | (158)
|}

In these equations <math display="inline">({^n\boldsymbol{a}},{^n\lambda })</math> is the last converged equilibrium point on the curve <math display="inline">(\boldsymbol{a},\lambda )</math>, whereas <math display="inline">({^{n+1}\boldsymbol{a}},{^{n+1}\lambda })</math> is the converged equilibrium point after an increment of the arc-length <math display="inline">\Delta S</math> on the same curve.

Furthermore, Ramm limited the length of the displacement increment for the prediction of each step with:

<span id="eq-159"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>[^{n+1}\delta \boldsymbol{a}^0]^T[^{n+1}\delta \boldsymbol{a}^0]=\Delta S^2 </math>
|}
| style="width: 5px;text-align: right;" | (159)
|}

where <math display="inline">^{n+1}\delta \boldsymbol{a}^0</math> is the displacement increment between the predicted point <math display="inline">{^{n+1}\boldsymbol{a}^0}</math> and the last converged equilibrium point <math display="inline">{^n\boldsymbol{a}}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\delta \boldsymbol{a}^0={^{n+1}\boldsymbol{a}^0}-{^n\boldsymbol{a}} </math>
|}
| style="width: 5px;text-align: right;" | (160)
|}

Equation ([[#eq-159|159]]) represents that the incremental displacement of the prediction must fall on the hypersphere of radius <math display="inline">\Delta S</math> (see Figure [[#img-30|30]]).

<div id='img-30'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-ramms_arc_length.png|600px|Ramm's method for a two degree of freedom system.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 30:''' Ramm's method for a two degree of freedom system.
|}

Thereby, with the modified Riks-Wempner method, the solution along the curve is found by following a line (or hyperplane) orthogonal to the prediction, which can be expressed by:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>[^{n+1}\delta \boldsymbol{a}^0]^T[^{n+1}\delta \boldsymbol{a}^{i+1}]=0 </math>
|}
| style="width: 5px;text-align: right;" | (161)
|}

After this short geometrical interpretation of the method, let us present the main steps of the algorithm.

Like we saw for the Newton-Raphson incremental-iterative technique, we start each step computing the prediction. However, since here we have the additional unknown <math display="inline">\lambda </math>, we need to obtain a prediction for both the load factor and the displacement field.

Given a prescribed radius <math display="inline">^{n+1}\Delta S</math>, the prediction increment for the load factor is:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\delta \lambda ^0=\frac{{^{n+1}\Delta }S}{\sqrt{[^{n+1}\delta \boldsymbol{a}_f^0]^T[^{n+1}\delta \boldsymbol{a}_f^0]}} </math>
|}
| style="width: 5px;text-align: right;" | (162)
|}

and the corresponding to the displacement field is:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\delta \boldsymbol{a}^0=[{^{n+1}\delta }\lambda ^0][{^{n+1}\delta }\boldsymbol{a}_f^0] </math>
|}
| style="width: 5px;text-align: right;" | (163)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\delta \boldsymbol{a}_f^0=[\boldsymbol{K}_{tan}(^n\boldsymbol{a})]^{-1}\boldsymbol{f} </math>
|}
| style="width: 5px;text-align: right;" | (164)
|}

Again, in the arc-length method we have two unknowns, and so we define the sought solution for the displacement field like in Newton-Raphson's method ([[#eq-151|151]]), and the one for the load factor like:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\lambda }^{i+1}={^{n+1}\lambda }^i+{^{n+1}\delta }\lambda ^{i+1} </math>
|}
| style="width: 5px;text-align: right;" | (165)
|}

Thereby, as proposed by <span id='citeF-56'></span>[[#cite-56|[56]]], we obtain <math display="inline">{^{n+1}\delta }\boldsymbol{a}^{i+1}</math> from two distinct parts:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\delta }\boldsymbol{a}^{i+1}={^{n+1}\delta }\boldsymbol{a}_R^{i+1}+[{^{n+1}\delta }\lambda ^{i+1}][{^{n+1}\delta }\boldsymbol{a}_f^{i+1}] </math>
|}
| style="width: 5px;text-align: right;" | (166)
|}

and we calculate <math display="inline">{^{n+1}\delta }\lambda ^{i+1}</math> from:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\delta }\lambda ^{i+1}=-\frac{[^{n+1}\delta \boldsymbol{a}^0]^T[{^{n+1}\delta }\boldsymbol{a}_R^{i+1}]}{[^{n+1}\delta \boldsymbol{a}^0]^T[{^{n+1}\delta }\boldsymbol{a}_f^{i+1}]} </math>
|}
| style="width: 5px;text-align: right;" | (167)
|}

where

<span id="eq-168"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\delta }\boldsymbol{a}_R^{i+1}=[\boldsymbol{K}_{tan}(^{n+1}\boldsymbol{a}^i)]^{-1}\boldsymbol{R}(^{n+1}\boldsymbol{a}^i) </math>
|}
| style="width: 5px;text-align: right;" | (168)
|}

and

<span id="eq-169"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>{^{n+1}\delta }\boldsymbol{a}_f^{i+1}=[\boldsymbol{K}_{tan}(^{n+1}\boldsymbol{a}^i)]^{-1}\boldsymbol{f} </math>
|}
| style="width: 5px;text-align: right;" | (169)
|}

Finally, at the end of each iteration, we check the convergence criterion with ([[#eq-149|149]]).

As a concluding remark, we must say that Ramm's arc-length method is not the most accurate arc-length strategy, but offers the advantage of simplicity, and guarantees a solution as long as it physically exists.

We must also stress the importance of the arc-length radius <math display="inline">\Delta S</math>, which controls the advancing velocity in the incremental-iterative process. The larger this radius, the faster is followed the equilibrium path. Nevertheless, near limit points with abrupt change in the response, an excessively large radius may lead to difficulties in tracing the load-displacement response. On the contrary, if the radius is chosen too small, the computational cost to obtain the solution path is very expensive. For this reason, in this work we have been estimating the radius with the relation proposed by Crisfield <span id='citeF-57'></span>[[#cite-57|[57]]].

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^{n+1}\Delta S={^n\Delta }S\left(\frac{I_d}{{^nI}}\right)^{1/2} </math>
|}
| style="width: 5px;text-align: right;" | (170)
|}

where <math display="inline">I_d</math> is an input parameter of the desired number of iterations per step, and <math display="inline">^nI</math> is the number of iterations at the previous step <math display="inline">n</math>.

The above radius control leads to small increments when the response is most non-linear, and large increments when the response is most linear.

Moreover, in the first step we estimate the radius from:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^0\Delta S=\sqrt{[^0\delta \boldsymbol{a}]^T[^0\delta \boldsymbol{a}]} </math>
|}
| style="width: 5px;text-align: right;" | (171)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>^0\delta \boldsymbol{a}=[\boldsymbol{K}_{tan}(^0\boldsymbol{a})]^{-1}\boldsymbol{R}(^0\boldsymbol{a}) </math>
|}
| style="width: 5px;text-align: right;" | (172)
|}

with <math display="inline">^0\boldsymbol{a}</math> being the initial guess of the displacement field.

==3.4 Relevant Aspects on the Implementation of the Damage Models==

In the following lines we will be reviewing how the main "ingredients" of the implemented arc-length strategy are obtained for the two implemented damage models: the partially regularized local damage model, and the fully regularized non-local damage model.

The main purpose of the section is to clarify the differences in the practical implementation of both models, pointing out those aspects that have not been mentioned in previous sections.

Looking at the expressions of the variables ([[#eq-168|168]]) and ([[#eq-169|169]]) defining the unknowns of the problem, one can see that at each iteration we need to know the tangent stiffness matrix <math display="inline">\boldsymbol{K}_{tan}(^{n+1}\boldsymbol{a}^i)</math>, the residual vector <math display="inline">\boldsymbol{R}(^{n+1}\boldsymbol{a}^i)</math> and the vector of external forces <math display="inline">\boldsymbol{f}</math>. Since we are working with conservative loads, this last vector is calculated only at the first step of the process. The tangent matrix and the residual vector, though, must be computed at each iteration.

The residual vector was presented in ([[#eq-155|155]]), with the vector of internal forces defined in terms of the secant stiffness matrix. In KRATOS, however, it is not computed using that form. Instead, the vector of internal forces is obtained by assembling the elementary contributions of an equivalent relation extracted from the Principle of Virtual Works ([[#eq-100|100]]):

<span id="eq-173"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{f}_{int}^{(e)}(\boldsymbol{a})=\iint _{A^{(e)}}\boldsymbol{B}^T\boldsymbol{\sigma }\,tdA </math>
|}
| style="width: 5px;text-align: right;" | (173)
|}

where <math display="inline">\boldsymbol{\sigma }</math> is the stress vector.

Thereby, in order to clearly understand the solution procedure followed in the implemented damage models, in the next paragraphs we will be assessing how to compute the stress vector and the tangent stiffness matrix at each iteration, for the partially regularized local damage model, and for the non-local damage approach.

===3.4.1 Partially regularized local damage model===

When computing the elementary contributions of the vector of internal forces ([[#eq-173|173]]), it is necessary to calculate the stress vector at each integration point.

In the partially regularized damage model, this procedure is very straightforward because it is performed like in any classical isotropic damage model. Thereby, we can obtain the stresses at any integration point from the locally evaluated strains and damage parameter (see Table [[#table-1|1]]).

As mentioned in the previous section, the expression for the elemental tangent stiffness matrix can be derived from the same definition of the internal forces ([[#eq-173|173]]):

<span id="eq-174"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\boldsymbol{K}_{tan}^{(e)} </math>
| <math> =\frac{\partial \boldsymbol{f}_{int}^{(e)}}{\partial \boldsymbol{a}}=\frac{\partial }{\partial \boldsymbol{a}}\left[\iint _{A^{(e)}}\boldsymbol{B}^T\boldsymbol{\sigma }\,tdA\right]</math>
|-
| style="text-align: right;" | 
| <math> =\iint _{A^{(e)}}\frac{\partial \boldsymbol{B}^T}{\partial \boldsymbol{a}}\boldsymbol{\sigma }\,tdA+\iint _{A^{(e)}}\boldsymbol{B}^T\frac{\partial \boldsymbol{\sigma }}{\partial \boldsymbol{a}}\,tdA </math>
|}
| style="width: 5px;text-align: right;" | (174)
|}

Thereby, since the problem we are solving falls into the material non-linearity category and we are working in a small deformation regime we have:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\partial \boldsymbol{B}^T}{\partial \boldsymbol{a}}=\boldsymbol{0} \quad ;\quad \frac{\partial \boldsymbol{\sigma }}{\partial \boldsymbol{a}}=\frac{\partial \boldsymbol{\sigma }}{\partial \boldsymbol{\varepsilon }}\frac{\partial \boldsymbol{\varepsilon }}{\partial \boldsymbol{a}} </math>
|}
| style="width: 5px;text-align: right;" | (175)
|}

Now, taking into account the relations ([[#eq-47|47]]) and ([[#eq-92|92]]) we can write:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\partial \boldsymbol{\sigma }}{\partial \boldsymbol{\varepsilon }}=\boldsymbol{E}_{tan} \quad ;\quad \frac{\partial \boldsymbol{\varepsilon }}{\partial \boldsymbol{a}}=\boldsymbol{B} </math>
|}
| style="width: 5px;text-align: right;" | (176)
|}

And so ([[#eq-174|174]]) can be rewritten as:

<span id="eq-177"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{tan}^{(e)}=\iint _{A^{(e)}}\boldsymbol{B}^T\boldsymbol{E}_{tan}\boldsymbol{B}\,tdA </math>
|}
| style="width: 5px;text-align: right;" | (177)
|}

where <math display="inline">\boldsymbol{E}_{tan}</math> is the tangent constitutive tensor that we presented in section [[#2.2.4 Elastic-damage tangent constitutive tensor|2.2.4]].

It is interesting to notice that from expression ([[#eq-177|177]]) one obtains the the tangent stiffness matrix, but only by substituting the tangent constitutive tensor by the secant one <math display="inline">\boldsymbol{E}_{sec}</math>, one obtains the secant stiffness matrix (see [[#eq-131|131]]).

Figure [[#img-31|31]] shows the general scheme followed in the local damage model for obtaining the stress vector and the tangent stiffness matrix at any iteration.

<div id='img-31'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-local_damage_scheme.png|600px|General scheme for the evaluation of the stress vector and the tangent stiffness matrix at any iteration of the local damage model.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 31:''' General scheme for the evaluation of the stress vector and the tangent stiffness matrix at any iteration of the local damage model.
|}

===3.4.2 Integral-type non-local damage model===

Numerical implementation of the non-local damage model based on averaging of equivalent strain is relatively simple. The evaluation of stresses remains explicit, and no internal iteration is needed. All that one needs is to implement the algorithm of weighted spatial averaging and, before damage is evaluated, replace the local equivalent strain by its non-local counterpart.

The values of non-local equivalent strain must be traced at individual Gauss integration points of the finite element model, because these are the points at which stresses need to be evaluated.

Thereby, the averaging integral in ([[#eq-63|63]]), defining the non-local equivalent strain at a given integration point, is evaluated numerically. Thus, considering the weighting function defined in ([[#eq-62|62]]), we have:

<span id="eq-178"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\tilde{\varepsilon }_{eq}(\boldsymbol{x}_p)=\sum _qw_q\alpha _{pq}\varepsilon _{eq}(\boldsymbol{x}_q) </math>
|}
| style="width: 5px;text-align: right;" | (178)
|}

where <math display="inline">w_q</math> is a coefficient containing the product of the determinant of the Jacobian, the thickness, and the integration weights of Gauss point <math display="inline">q</math>; <math display="inline">\boldsymbol{x}_q</math> are the coordinates of the integration point; and <math display="inline">\alpha _{pq}</math> is the weight of non-local interaction between points <math display="inline">p</math> and <math display="inline">q</math>, defined as:

<span id="eq-179"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\alpha _{pq}=\frac{\alpha _0(\| \boldsymbol{x}_p-\boldsymbol{x}_q\| )}{\sum _rw_r\alpha _0(\| \boldsymbol{x}_p-\boldsymbol{x}_r\| )} </math>
|}
| style="width: 5px;text-align: right;" | (179)
|}

In the previous two equations, subscript <math display="inline">p</math> represents the receiver point under consideration, that can be any integration point in the solid domain, whereas indexes <math display="inline">q</math> and <math display="inline">r</math> correspond to source points. Furthermore, since the chosen weighting function <math display="inline">\alpha _0</math> has bounded support ([[#eq-59|59]]), <math display="inline">\alpha _{pq}</math> vanishes if the distance <math display="inline">\| \boldsymbol{x}_p-\boldsymbol{x}_q\| </math> between points <math display="inline">p</math> and <math display="inline">q</math> is larger than the interaction radius <math display="inline">R</math>. Therefore, the sums in ([[#eq-178|178]]) and ([[#eq-179|179]]) do not need to be taken over all Gauss points, but only over those that are located inside the domain of influence of point <math display="inline">p</math>. In this regard we must notice that, in order to account for the non-local interaction, at the damage process zone one must always use an element size smaller than the interaction radius. Otherwise the damage model would become local.

Therefore, each Gauss point must have a non-local interaction table that gives access to its neighbours. This table must be constructed, at the beginning of the problem, from a search of non-local neighbours.

In this work, the search of neighbours is performed by means of a grid-based algorithm. A general rectangular grid is defined in the entire domain and all the integration points are positioned in the cells. This way, the neighbour search that must be performed for each Gauss point is restricted to a limited number of cells, i.e., the ones that fall inside the interaction radius of the considered point (see Figure [[#img-32|32]]).

<div id='img-32'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-neighbour_search.png|600px|Grid-based non-local search.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 32:''' Grid-based non-local search.
|}

The stress evaluation procedure, repeatedly called during the incremental-iterative strategy, makes use of the non-local interaction tables when the non-local equivalent strain is computed. To obtain <math display="inline">\tilde{\varepsilon }_{eq}</math>, we first compute the local equivalent strains at all Gauss points, and then we calculate the non-local counterpart using ([[#eq-178|178]]).

With regard to the stiffness matrix, if one aims to obtain an iterative solver with quadratic convergence when working with a non-local damage model, it is necessary to construct the tangent stiffness matrix in a consistent manner.

Let us start rewriting the expression for the vector of internal forces ([[#eq-173|173]]) as a numerical integration over Gauss points:

<span id="eq-180"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{f}_{int}=\sum _pw_p\boldsymbol{B}^T(\boldsymbol{x}_p)\boldsymbol{\sigma }(\boldsymbol{x}_p) </math>
|}
| style="width: 5px;text-align: right;" | (180)
|}

Using the stress-strain law ([[#eq-128|128]]) and the standard strain-displacement relation ([[#eq-92|92]]), we can expand ([[#eq-180|180]]) as follows:

<span id="eq-181"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\boldsymbol{f}_{int} </math>
| <math> =\sum _pw_p(1-d_p)\boldsymbol{B}_p^T\boldsymbol{E}\boldsymbol{B}_p\boldsymbol{a} </math>
|-
| style="text-align: right;" | 
| <math> =\sum _pw_p(1-d_p)\boldsymbol{K}_p\boldsymbol{a} </math>
|}
| style="width: 5px;text-align: right;" | (181)
|}

where <math display="inline">\boldsymbol{K}_p=\boldsymbol{B}_p^T\boldsymbol{E}\boldsymbol{B}_p</math>. In the virgin state <math display="inline">d_p=0</math> at all integration points and so ([[#eq-181|181]]) can be rewritten as <math display="inline">\boldsymbol{f}_{int}=\boldsymbol{K}\boldsymbol{a}</math>, where <math display="inline">\boldsymbol{K}=\sum _pw_p\boldsymbol{K}_p</math> is the elastic stiffness matrix.

As we have already said, the tangent stiffness matrix is obtained by differentiating the internal forces with respect to the nodal displacements. Since the damage variable depends on the nodal displacements through the equivalent strain, we will need to compute first the derivative of the damage variable with respect to the displacement vector <math display="inline">\boldsymbol{a}</math>. Taking into account that <math display="inline">d=g(r)</math> ([[#eq-7|7]]), <math display="inline">r</math> depends on <math display="inline">\tilde{\varepsilon }_{eq}</math> ([[#eq-64|64]]), and <math display="inline">\tilde{\varepsilon }_{eq}</math> depends on <math display="inline">\boldsymbol{a}</math> through the interpolated strains, we can use the chain rule to write the derivative of the damage variable for an integration point <math display="inline">p</math> as follows:

<span id="eq-182"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\partial d_p}{\partial \boldsymbol{a}}=\frac{dg}{dr_p}\frac{dr_p}{d\tilde{\varepsilon }_{eq,p}}\frac{\partial \tilde{\varepsilon }_{eq,p}}{\partial \boldsymbol{a}}=g'_pt_p\frac{\partial \tilde{\varepsilon }_{eq,p}}{\partial \boldsymbol{a}} </math>
|}
| style="width: 5px;text-align: right;" | (182)
|}

where <math display="inline">g'_p</math> is the derivative of the damage evolution law with respect to the internal state variable <math display="inline">r</math>, and <math display="inline">t_p</math> is the loading-unloading factor that is 0 in an elastic loading or in an unloading regime, and 1 in a loading regime with growing damage:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>t=\frac{dr}{d\tilde{\varepsilon }_{eq}}=\left\{\begin{array}{ll}0 & \hbox{if }\tilde{\varepsilon }_{eq}<r\hbox{ } \\ 1 & \hbox{if }\tilde{\varepsilon }_{eq}=r\hbox{ and }\dot{\tilde{\varepsilon }}_{eq}\ge 0\hbox{ } \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;" | (183)
|}

Now, from expression ([[#eq-178|178]]), we can differentiate the non-local equivalent strain of a Gauss point <math display="inline">p</math> with respect to the nodal displacements:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\partial \tilde{\varepsilon }_{eq,p}}{\partial \boldsymbol{a}}=\sum _qw_q\alpha _{pq}\frac{\partial \varepsilon _{eq,q}}{\partial \boldsymbol{a}} </math>
|}
| style="width: 5px;text-align: right;" | (184)
|}

And by using the chain rule and the relation in ([[#eq-92|92]]) we obtain:

<span id="eq-185"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\frac{\partial \tilde{\varepsilon }_{eq,p}}{\partial \boldsymbol{a}}=\sum _qw_q\alpha _{pq}\frac{\partial \varepsilon _{eq,q}}{\partial \boldsymbol{\varepsilon }_q}\frac{\partial \boldsymbol{\varepsilon }_q}{\partial \boldsymbol{a}}=\sum _qw_q\alpha _{pq}\frac{\partial \varepsilon _{eq,q}}{\partial \boldsymbol{\varepsilon }_q}\boldsymbol{B}_q </math>
|}
| style="width: 5px;text-align: right;" | (185)
|}

the derivative of the equivalent strain with respect to the strain components <math display="inline">\partial \varepsilon _{eq,q}/\partial \boldsymbol{\varepsilon }_q</math> is a row matrix that will depend on the chosen form of equivalent strain.

Let us now rewrite the expression of the internal forces in ([[#eq-181|181]]) as follows:

<span id="eq-186"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{f}_{int}=\sum _pw_p\boldsymbol{K}_p\boldsymbol{a}-\sum _pw_p\boldsymbol{K}_p\boldsymbol{a}d_p </math>
|}
| style="width: 5px;text-align: right;" | (186)
|}

At this point, we can already differentiate ([[#eq-186|186]]) and substitute ([[#eq-182|182]]) to get a first expression for the tangent stiffness matrix.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\boldsymbol{K}_{tan} </math>
| <math> =\frac{\partial \boldsymbol{f}_{int}}{\partial \boldsymbol{a}}=\sum _pw_p\boldsymbol{K}_p-\left(\sum _pw_p\boldsymbol{K}_pd_p+\sum _pw_p\boldsymbol{K}_p\boldsymbol{a}\frac{\partial d_p}{\partial \boldsymbol{a}}\right)</math>
|-
| style="text-align: right;" | 
| <math> =\sum _pw_p(1-d_p)\boldsymbol{K}_p-\sum _pw_p\boldsymbol{K}_p\boldsymbol{a}\frac{\partial d_p}{\partial \boldsymbol{a}} </math>
|-
| style="text-align: right;" | 
| <math> =\sum _pw_p\boldsymbol{B}_p^T(1-d_p)\boldsymbol{E}\boldsymbol{B}_p-\sum _pw_p\boldsymbol{B}_p^T\boldsymbol{E}\boldsymbol{B}_p\boldsymbol{a}g'_pt_p\frac{\partial \tilde{\varepsilon }_{eq,p}}{\partial \boldsymbol{a}} </math>
|}
| style="width: 5px;text-align: right;" | (187)
|}

Note that the first term is the secant stiffness matrix, that coincides with the tangent matrix in an elastic loading or in an unloading regime (<math display="inline">t_p=0</math> <math display="inline">\forall p</math>). The term on the right is the non-local part of the tangent stiffness matrix which, substituting ([[#eq-185|185]]), results:

<span id="eq-188"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\boldsymbol{K}_{tan} </math>
| <math> =\boldsymbol{K}_{sec}-\sum _pw_p\boldsymbol{B}_p^T\boldsymbol{E}\boldsymbol{\varepsilon }_pg'_pt_p\sum _qw_q\alpha _{pq}\frac{\partial \varepsilon _{eq,q}}{\partial \boldsymbol{\varepsilon }_q}\boldsymbol{B}_q </math>
|-
| style="text-align: right;" | 
| <math> =\boldsymbol{K}_{sec}-\sum _{p,q}w_p\boldsymbol{B}_p^T\boldsymbol{\bar{\sigma }}_pg'_pt_pw_q\alpha _{pq}\frac{\partial \varepsilon _{eq,q}}{\partial \boldsymbol{\varepsilon }_q}\boldsymbol{B}_q </math>
|}
| style="width: 5px;text-align: right;" | (188)
|}

Defining for convenience the column matrix <math display="inline">\boldsymbol{b}^c_p=\boldsymbol{B}^T_p\boldsymbol{\bar{\sigma }}_p</math>, the row matrix <math display="inline">\boldsymbol{b}^r_q=\displaystyle \frac{\partial \varepsilon _{eq,q}}{\partial \boldsymbol{\varepsilon }_q}\boldsymbol{B}_q</math>, and the coefficient <math display="inline">w_{pq}=w_pw_q\alpha _{pq}</math>, ([[#eq-188|188]]) can be rewritten like:

<span id="eq-189"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{tan}=\boldsymbol{K}_{sec}-\sum _{p,q}w_{pq}g'_pt_p\boldsymbol{b}^c_p\boldsymbol{b}^r_q </math>
|}
| style="width: 5px;text-align: right;" | (189)
|}

The double index of the sum, caused by the non-local interaction, implies that the term on the right part of the expression ([[#eq-189|189]]) cannot be assembled from elementary contributions solely. Essentially, each pair of Gauss points <math display="inline">p</math> and <math display="inline">q</math> contributes to the global stiffness matrix with a block of the same size as that of the classical element stiffness matrix. The difference is that the assembling routine differs from the usual one because in this case one needs to take into account the elements of both points <math display="inline">p</math> and <math display="inline">q</math> (see Figure [[#img-33|33]]). In consequence, the global stiffness matrix is always non-symmetric and its bandwidth increases due to the non-local interaction.

<div id='img-33'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-nonlocal_assembly.png|600px|Assembly of non-local contributions to the global tangent stiffness matrix.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 33:''' Assembly of non-local contributions to the global tangent stiffness matrix.
|}

In order to avoid the additional non-zero entries that the non-local interaction introduces into the global stiffness matrix, some authors neglect the non-local terms by using <math display="inline">\alpha _{pq}=\delta _{pq}/w_p</math>, where <math display="inline">\delta _{pq}</math> is the Kronecker delta, and so <math display="inline">w_{pq}=\delta _{pq}w_q</math>. This way, equation ([[#eq-189|189]]) reduces to:

<span id="eq-190"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{tan}^{local}=\boldsymbol{K}_{sec}-\sum _{p}w_{p}g'_pt_p\boldsymbol{b}^c_p\boldsymbol{b}^r_p </math>
|}
| style="width: 5px;text-align: right;" | (190)
|}

in essence, the sum is performed over one index only. However, it must be noticed that the resulting local tangent matrix <math display="inline">\boldsymbol{K}_{tan}^{local}</math> is no longer consistent, and quadratic convergence is lost. Furthermore, although the secant matrix is always symmetric, the corrective term on the right part of expression ([[#eq-190|190]]) is symmetric only if the column matrix <math display="inline">\boldsymbol{b}^c</math> and the row matrix <math display="inline">\boldsymbol{b}^r</math> are linearly dependent, which is the case if <math display="inline">\partial \varepsilon _{eq}/\partial \boldsymbol{\varepsilon }</math> is a scalar multiple of the effective stresses <math display="inline">\boldsymbol{\bar{\sigma }}</math>. This is only true for the definition of the equivalent strain based on the damage energy release rate ([[#eq-19|19]]).

Therefore, the most important issue caused by non-locality is the evolutionary character of the profile of the stiffness matrix. For the simulation of quasi-brittle materials like concrete, the consistent stiffness matrix remains local through the elastic branch, and so the initial distribution of non-zero entries is the same as in the local case. Nevertheless, when the damage threshold is exceeded and the damage zone starts propagating, new non-zero entries appear due to the non-local interaction between Gauss points belonging to different elements, and the profile of the stiffness matrix must be dynamically adapted. The number of additional non-zero entries will depend on each particular case, but if a finer mesh is used in the expected softening zones, this number can be relatively high.

In conclusion, it must be noticed that non-local averaging certainly increases the computational cost with respect to the corresponding local model, but since the non-local model completely removes the pathological sensitivity to the mesh size, and partially alleviates the mesh-induced directional bias, this extra effort is indeed worthwhile.

Figure [[#img-34|34]] shows the general scheme followed in the non-local damage model for obtaining the stress vector and the tangent stiffness matrix at any iteration.

<div id='img-34'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-nonlocal_damage_scheme.png|600px|General scheme for the evaluation of the stress vector and the tangent stiffness matrix at any iteration of the non-local damage model.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 34:''' General scheme for the evaluation of the stress vector and the tangent stiffness matrix at any iteration of the non-local damage model.
|}

==3.5 Mesh-Adaptive Technique==

Non-local models lead to smooth solutions with a continuous variation of strain. However, to resolve narrow bands of highly localized strains, it is necessary to use sufficiently fine computational grids. Fortunately, the mesh must be fine only in the process zone, while the remaining part of the structure can be reasonably well represented by a coarser mesh. In general, the localization pattern is not known in advance, and it is actually tedious to construct suitably refined meshes by hand. Thereby, efficiency of the analysis can be greatly increased by means of an adaptive technique, which automates the whole process.

In this regard, one of the last implementations of this work has been an adaptive mesh refinement technique. Although still under development, it has already shown some interesting results, and so in this section we wanted to review the most important features of the procedure.

The starting point for the implementation of this adaptive technique was the work done by Miquel Portabella in its minor thesis <span id='citeF-58'></span>[[#cite-58|[58]]]. The adaptive mesh refinement technique of Portabella was aimed at static-linear-elastic two-dimensional cases, and so we had to adapt it to incremental-iterative damage mechanics problems. Due to the importance of achieving robust results regardless of the mesh, we have designed this adaptive technique for the non-local damage method.

In the following lines, we will present the general structure of the adaptive procedure, and then its basic stages will be briefly discussed.

The general algorithm of non-linear adaptive analysis can be described as follows. After reaching the equilibrium state and updating the solution, an error estimation is performed in order to evaluate the error distribution over the mesh. Then, a remeshing criterion uses the information about error distribution and determines the required mesh density. From this analysis, we can obtain a new spatial discretization using a mesh generator interface.

In a truly adaptive approach, after generating a new discretization, the data structures corresponding to the newly generated mesh are created, and the transfer of displacements and internal variables from the old mesh to the new one is performed. After the mapping, the internal variables are used together with the strain computed from the mapped displacements to update the internal state of each integration point on the new mesh (to achieve local consistency). Once the transfer has finished, the old discretization is deleted and the mapped configuration is brought into global equilibrium through iteration. Afterwards, the solution continues with the next incremental-iterative step.

Another possibility is to restart the analysis from the initial state after the new discretization is generated. This approach does not require the transfer of the current state from the old discretization to the new one, but from the computational point of view is less efficient than the truly adaptive approach, especially if the remeshing is done frequently.

Let us now present the main stages of the implemented adaptive procedure, pointing out the most relevant aspects of each one:

* '''Error Estimation'''

The implemented error estimator is based on the stress evaluation. Thus, let us define the error as the difference between an exact value of effective stresses <math display="inline">\bar{\boldsymbol{\sigma }}_e</math> and an approximate one <math display="inline">\bar{\boldsymbol{\sigma }}_a</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\epsilon }_{\bar{\sigma }}=\bar{\boldsymbol{\sigma }}_e-\bar{\boldsymbol{\sigma }}_a </math>
|}
| style="width: 5px;text-align: right;" | (191)
|}

Since the Finite Element Method is an approximate technique, it is not possible to know the exact value of stresses. Thereby, we must estimate the error and substitute the exact value of effective stresses by a "reasonably good" value of stresses. This reasonably good value is obtained from the extrapolation of the stresses to the nodes, and posterior smoothing <span id='citeF-59'></span>[[#cite-59|[59]]]. The approximate value of stresses is simply the default evaluated stresses of the FEM solution. Therefore, the estimated error can be defined as the difference between the smoothed effective stresses <math display="inline">\bar{\boldsymbol{\sigma }}_s</math>, and the calculated effective stresses <math display="inline">\bar{\boldsymbol{\sigma }}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\epsilon }_{\bar{\sigma }}\approx \bar{\boldsymbol{\sigma }}_s-\bar{\boldsymbol{\sigma }} </math>
|}
| style="width: 5px;text-align: right;" | (192)
|}

However, we will be working with integral measures of the error that take into account the error over all the element, and so we define the energetic norm of the error over an element as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: right;" | <math>\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| ^{(e)} </math>
| <math>=\left[\int _{\Omega ^{(e)}}\boldsymbol{\epsilon }_{\bar{\sigma }}^T\boldsymbol{E}^{-1}\boldsymbol{\epsilon }_{\bar{\sigma }}d\Omega \right]^{1/2} </math>
|-
| style="text-align: right;" | 
| <math> =\left[\int _{\Omega ^{(e)}}[\bar{\boldsymbol{\sigma }}_s-\bar{\boldsymbol{\sigma }}]^T\boldsymbol{E}^{-1}[\bar{\boldsymbol{\sigma }}_s-\bar{\boldsymbol{\sigma }}]d\Omega \right]^{1/2} </math>
|}
| style="width: 5px;text-align: right;" | (193)
|}

And, the square of the global error can be obtained from the sum of the squares of all the elemental errors:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| ^{2}=\sum _e(\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| ^{(e)})^2 </math>
|}
| style="width: 5px;text-align: right;" | (194)
|}

The smaller is the distance between the nodes of the mesh, the smaller will be the difference between the smoothed and non-smoothed stresses. Thereby, the presented energetic norm tends to zero as the size of the element diminishes, i.e., <math display="inline">\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| =\mathcal{O}(h^m)</math> where <math display="inline">h</math> is the element size and <math display="inline">m</math> is the order of the polynomial defining the shape functions.

* '''Remeshing Criterion'''

Determining whether a mesh must be refined or not requires to previously define some quality conditions based on the estimated error.

First, the energetic norm of the global error should be smaller than a certain percentage of the deformation energy:

<span id="eq-195"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| \le \eta \,U </math>
|}
| style="width: 5px;text-align: right;" | (195)
|}

where <math display="inline">\eta </math> is the permissible percentage of global error, and the deformation energy is obtained from:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>U^{2}=\sum _e(U^{(e)})^2 </math>
|}
| style="width: 5px;text-align: right;" | (196)
|}

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>U^{(e)}=\left[\int _{\Omega ^{(e)}}\bar{\boldsymbol{\sigma }}_s^T\boldsymbol{E}^{-1}\bar{\boldsymbol{\sigma }}_sd\Omega \right]^{1/2} </math>
|}
| style="width: 5px;text-align: right;" | (197)
|}

To determine whether condition ([[#eq-195|195]]) is fulfilled, it is convenient to work with a global error parameter defined as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\xi _g=\frac{\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| }{\eta \,U} </math>
|}
| style="width: 5px;text-align: right;" | (198)
|}

Thereby, when <math display="inline">\xi _g=1</math> the global error condition is perfectly fulfilled, for <math display="inline">\xi _g>1</math> mesh should be refined, and for <math display="inline">\xi _g<1</math> the mesh size could be larger.

Furthermore, taking into account that <math display="inline">\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| =\mathcal{O}(h^m)</math>, the new element size <math display="inline">\hat{h}^{(e)}</math> can be obtained with:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\hat{h}^{(e)}=\frac{h^{(e)}}{(\xi _g)^{\frac{1}{m}}} </math>
|}
| style="width: 5px;text-align: right;" | (199)
|}

If one aims at obtaining a selective refinement method, apart from ([[#eq-195|195]]), another condition concerning the error of each element must be simultaneously imposed. Depending on this local error condition, one can define different approaches <span id='citeF-60'></span><span id='citeF-61'></span>[[#cite-60|[60,61]]]. In this work, we have been working with a ''remeshing criterion based on global error equidistribution''.

This criterion distributes the global error uniformly among all the elements of the mesh, and so the elemental error should accomplish the following condition:

<span id="eq-200"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| ^{(e)}\le \frac{\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| }{\sqrt{n}} </math>
|}
| style="width: 5px;text-align: right;" | (200)
|}

where <math display="inline">n</math> is the number of elements of the mesh.

Like in the case of the global error, we can work with a local error parameter:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\xi _l^{(e)}=\frac{\sqrt{n}\,\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| ^{(e)}}{\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| } </math>
|}
| style="width: 5px;text-align: right;" | (201)
|}

with the same meaning that in the global case: <math display="inline">\xi _l^{(e)}=1</math> indicates an optimal element size, whereas <math display="inline">\xi _l^{(e)}>1</math> and <math display="inline">\xi _l^{(e)}<1</math> imply that the element is too large and too small, respectively.

This local parameter allows to define the new element size that accomplish ([[#eq-200|200]]):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\hat{h}^{(e)}=\frac{h^{(e)}}{(\xi _l^{(e)})^{\frac{2}{2m+d}}} </math>
|}
| style="width: 5px;text-align: right;" | (202)
|}

where <math display="inline">d</math> is the dimension of the problem.

In the end, the remeshing criterion results from the combination of the global error condition ([[#eq-195|195]]) and the local one ([[#eq-200|200]]). In consequence, the final refinement parameter of the element can be defined as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\xi ^{(e)}=\xi _g\,\xi _l^{(e)}=\frac{\sqrt{n}\,\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| ^{(e)}}{\eta \,U} </math>
|}
| style="width: 5px;text-align: right;" | (203)
|}

And the new element size is obtained with:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\hat{h}^{(e)}=\frac{h^{(e)}}{\gamma ^{(e)}} </math>
|}
| style="width: 5px;text-align: right;" | (204)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;" 
|-
| style="text-align: center;" | <math>\gamma ^{(e)}=(\xi _g)^{\frac{1}{m}}\,(\xi _l^{(e)})^{\frac{2}{2m+d}} </math>
|}
| style="width: 5px;text-align: right;" | (205)
|}

* '''Mesh Generator Interface'''

The mesh generator interface is the one from GiD, a pre-processing and post-processing software. Therefore, not only GiD meshes the geometry in the beginning of a problem, but it is also GiD which allows to obtain the new spatial discretization every time we need to adapt the mesh.

Thereby, after the error estimator and the remeshing criterion are applied, a "background mesh" file (*.bgm) is generated with the information of the new element sizes. Then, GiD allows loading this background mesh file along with the original mesh file so as to generate a new mesh according to the refinement parameter of each old element.

Furthermore, before starting any problem, one can always change GiD meshing preferences in order to modify the velocity of mesh transition and other parameters that will affect the original mesh as well as the subsequent adapted meshes.

* '''Mapping of Variables'''

Before one can continue with the solution of a problem, there is a last stage that must be performed in a truly adaptive approach: the mapping of primary unknowns and internal variables.

In essence, if one aims at continuing the analysis from the current state, instead of restarting it from scratch after every mesh refinement, it is necessary to apply some transfer algorithms for the displacements and internal history variables (in the present case, the state variable <math display="inline">r</math> governing the damage evolution).

Mapping of the primary unknowns (nodal displacements) is achieved by using the shape function projections. To do so, we must first place each new node inside an old element, by means of a grid-based search, and then we interpolate the displacements of the old nodes to the new one (see Figure [[#img-35|35]].a).

On the other hand, mapping of the internal state variables is done through a weighted spatial averaging, similar to the one used for the computation of the non-local equivalent strain. The difference is that, in this case, the source points are the integration points of the old mesh, and the receiver points are the integration points of the new mesh (see Figure [[#img-35|35]].b). Like before, another grid-based search is performed in order to determine the Gauss points of the old mesh that fall inside the interaction radius <math display="inline">R</math> of each new Gauss point.

<div id='img-35'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-mapping.png|600px|Mapping of variables: (a)Nodal displacements mapping. (b)Internal state variables mapping.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 35:''' Mapping of variables: (a)Nodal displacements mapping. (b)Internal state variables mapping.
|}

=4 Examples=

==4.1 Introduction==

At this point, we have already presented the fundamental concepts on damage mechanics theory, as well as the the most relevant aspects concerning the implementation of a damage model within the Finite Element Method.

The present chapter will be devoted to test and validate the implemented code, by solving two classical examples in the damage mechanics field: the three-point bending test, and the single-edge notched beam test.

For each one of the tests, two experiments will be performed. First, we will solve the problem with the two implemented damage models: the partially regularized local damage model, and the integral-type non-local damage model. The objective of this experiment is to analyse and compare both models, pointing out the strengths and limitations of the codes, and try to confirm that the non-local procedure is a valid and robust approach.

After that, we will solve the problem again, applying the implemented mesh-adaptive procedure to the non-local damage model. As we have already commented, this technique is still under development, but the examples performed will help us understand better the adaptive procedure and establish the next steps in its future evolution.

Both examples will be carried out in 2D, which is actually advantageous when testing a new code, given its reduced computational cost as compared to the corresponding three-dimensional case. Thereby, as we stated before, we will be using three-node triangular elements and four-node quadrilateral elements to represent the discretized domain. Regarding the solving strategy, we will be working with the arc-length method presented in the previous chapter, and thus we will be following the equilibrium path of the problems for a set of steps. Apart from that, self weight will not be considered in any of the examples.

Now, in order to properly understand the scheme followed for the performance of the examples, let us briefly present GiD.

GiD is a software for the pre-processing and post-processing of simulations. It has CAD tools for drawing the geometry, as well as a mesh generator for obtaining the spatial discretization of the model. Furthermore, GiD allows embedding customized modules specifically designed for a problem of interest, called "Problem Types", which enhance GiD pre-processing capacities so as to define all the relevant components of the model, i.e., boundary conditions, loads, material properties and computation parameters. In addition to this, the Problem Types are usually linked to the compiled KRATOS application that will solve the problem. Thereby, by using GiD along with an specific Problem Type, one can prepare a model and run the simulation from GiD itself. Lastly, GiD also offers post-processing tools, and so when the computation finishes, one can immediately display the obtained results in GiD. In essence, the solution of the problems follows the scheme in Figure 36.

<div id='img-36'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-examples_scheme.png|600px|Global structure for the solution of a problem.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 36:''' Global structure for the solution of a problem.
|}

==4.2 Three-Point Bending Test==

This test is performed with a notched beam subjected to three-point bending (TPB). The beam has a square cross section of 40 <math display="inline">\times </math> 320 mm, a span of 1280 mm, and the notch is 3 mm thick and extends over one tenth of the beam depth (see Figure [[#img-37|37]]).

<div id='img-37'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-three-point_bending_test.png|480px|TPB test. Problem statement. Distances in mm]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 37:''' TPB test. Problem statement. Distances in mm
|}

In order to reduce fictitious stress concentrations at the supports and at the loading point, both the displacement constraints and the load are imposed over a line of 20 mm, as shown in Figure [[#img-38|38]].

<div id='img-38a'></div>
<div id='img-38b'></div>
<div id='img-38'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-displacement_bc_tpbt.png|480px|Displacement constraints.]]
|- style="text-align: center; font-size: 75%;"
| (a) Displacement constraints.
|-
|[[Image:draft_Content_981184740-loads_tpbt.png|480px|Loads.]]
|- style="text-align: center; font-size: 75%;"
| (b) Loads.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 38:''' TPB test. Boundary Conditions.
|}

In this case, plane stress conditions have been assumed, and the geometry will be meshed by means of standard four-node quadrilaterals with 2 <math display="inline">\times </math> 2 integration points.

===4.2.1 Damage models testing===

In this first experiment we will solve the problem using the two implemented damage models. On the one hand, the partially regularized local damage approach will make use of the modified version of the equivalent strain proposed by Simo and Ju ([[#eq-33|33]]), and the damage evolution law presented in ([[#eq-39|39]]). On the other hand, the non-local damage approach, will be defined with the Mazars model, regarding the equivalent strain ([[#eq-25|25]]) and the damage evolution law ([[#eq-42|42]]). The weighting function for the non-local interaction will depend on a Gauss distribution function of bounded support ([[#eq-59|59]]).

The material parameters for the local damage model are summarized in Table [[#table-5|5]] and have been obtained after various attempts, trying to fit the experimental results in <span id='citeF-62'></span>[[#cite-62|[62]]]. For the non-local approach we are going to use the material parameters obtained in <span id='citeF-62'></span>[[#cite-62|[62]]] (see Table [[#table-6|6]]).


{| class="wikitable" style="text-align: left; margin: 1em auto;"
|+ <span id='table-5'></span>Table. 5 TPB test. Material parameters for the Simo-Ju local model.
|-
|  '''Parameter''' 
| '''Value''' 
|-
|  Young's modulus (<math display="inline">E</math>) 
| 38500 MPa 
|-
| Poisson's ratio (<math display="inline">\nu </math>) 
| 0.24 
|-
| Compressive strength (<math display="inline">\sigma _y^c</math>) 
| 45 MPa 
|-
| Tensile strength (<math display="inline">\sigma _y^t</math>) 
| 3.8 MPa 
|-
| Fracture energy (<math display="inline">G_f</math>) 
| 100 <math display="inline">J/m^2</math> 
|-
| Limit fracture length (<math display="inline">l_{lim}</math>) 
| 5 mm 

|}


{| class="wikitable" style="text-align: left; margin: 1em auto;"
|+ <span id='table-6'></span>Table. 6 TPB test. Material parameters for the Mazars non-local model.
|-
|  '''Parameter''' 
| '''Value''' 
|-
|  Young's modulus (<math display="inline">E</math>) 
| 38500 MPa 
|-
| Poisson's ratio (<math display="inline">\nu </math>) 
| 0.24 
|-
| Damage threshold (<math display="inline">r_0</math>) 
| <math display="inline">3\cdot{10}^{-5}</math> 
|-
| Parameter <math display="inline">A</math> in compression (<math display="inline">A_c</math>) 
| 1.25 
|-
| Parameter <math display="inline">B</math> in compression (<math display="inline">B_c</math>) 
| 1000 
|-
| Parameter <math display="inline">A</math> in tension (<math display="inline">A_t</math>) 
| 0.95 
|-
| Parameter <math display="inline">B</math> in tension (<math display="inline">B_t</math>) 
| 9000 
|-
| Characteristic length (<math display="inline">l_c</math>) 
| 40 mm 

|}

In order to assess the robustness of the models in terms of mesh sensitivity, we are going to solve the problem for different spatial discretizations. Thereby, here we will be using three unstructured meshes of quadrilaterals with a minimum size of 15 mm, 7 mm, and 3 mm, respectively (see Figure [[#img-39|39]]).

<div id='img-39a'></div>
<div id='img-39b'></div>
<div id='img-39c'></div>
<div id='img-39'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_m1.png|400px|Mesh 1: 2024 elements, 2191 nodes.]]
|[[Image:draft_Content_981184740-tpbt_m2.png|400px|Mesh 2: 2679 elements, 2859 nodes.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1: 2024 elements, 2191 nodes.
| (b) Mesh 2: 2679 elements, 2859 nodes.
|-
| colspan="2"|[[Image:draft_Content_981184740-tpbt_m3.png|300px|Mesh 3: 6543 elements, 6772 nodes.]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (c) Mesh 3: 6543 elements, 6772 nodes.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 39:''' TPB test. Spatial discretizations.
|}

Figure [[#img-40|40]] shows the relation between the applied load and the vertical deflection of the beam. As one can see from the discontinued equilibrium curves, we had serious difficulties in tracing the response of a full test. Once the limit point at the peak was surpassed, the ratio of residual forces would start oscillating without converging to the prescribed tolerance, even when reducing the arc-length radius. A stability study could reveal the presence of bifurcation points in the post-peak zone that should be carefully analysed in order to distinguish the stable branch (physically possible) from the unstable ones. However, another reason behind the convergence problems could be the use of a too global arc-length method, like the one implemented. Indeed, to account for the localized nature of quasi-brittle failure, a more specific control parameter, like the Crack Mouth Opening Displacement (CMOD), could help improving the convergence near snap-back zones.

That aside, if we look at the curves in Figure [[#img-40a|40a]] we can see that, although there is no relevant difference between the response obtained with meshes 1 and 2, the peak actually decreases with the finest mesh 3, and so the total dissipated energy. On the other hand, the response diagram in Figure [[#img-40b|40b]] shows practically the same peak load for the three meshes.

Thereby, only looking at the depicted curves, it seems that the partially regularized local damage model is more sensitive to changes in the spatial discretization than the non-local model.

<div id='img-40a'></div>
<div id='img-40b'></div>
<div id='img-40'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 75%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_local.png|450px|Partially regularized local damage model.]]
|- style="text-align: center; font-size: 75%;"
| (a) Partially regularized local damage model.
|-
|[[Image:draft_Content_981184740-tpbt_nonlocal.png|450px|Non-local damage model.]]
|- style="text-align: center; font-size: 75%;"
| (b) Non-local damage model.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 40:''' TPB test. Force-vertical deflection diagrams.
|}

However, in order to understand better the response of both models, let us present some snapshots with the damage distribution.

Figures [[#img-41|41]] and [[#img-43|43]] show an initial stage of damage progression in the local and non-local model, respectively. Just by looking at the size and shape of the damaged zone, we can clearly state the most differential trait of each model: localization on the one hand, and diffusion on the other. This concept is crucial to understand the behaviour seen in the load-deflection diagrams of Figure [[#img-40|40]].

Essentially, in the local approach, damage starts at the most stressed element and then "jumps" to the next one when the first is totally damaged. As a result, shape and direction of progression of damage strongly depends on the size and distribution of the elements of the mesh (see Figure [[#img-42|42]]).

<div id='img-41a'></div>
<div id='img-41b'></div>
<div id='img-41c'></div>
<div id='img-41'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 90%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_loc_m1_d.png|540px|Mesh 1.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1.
|-
|[[Image:draft_Content_981184740-tpbt_loc_m2_d.png|540px|Mesh 2.]]
|- style="text-align: center; font-size: 75%;"
| (b) Mesh 2.
|-
|[[Image:draft_Content_981184740-tpbt_loc_m3_d.png|540px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 3.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 41:''' TPB test. Damage initiation in the local model.
|}

<div id='img-42a'></div>
<div id='img-42b'></div>
<div id='img-42'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_loc_d_m1_z.png|400px|Mesh 1.]]
|[[Image:draft_Content_981184740-tpbt_loc_d_m3_z.png|400px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1.
| (b) Mesh 3.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 42:''' TPB test. Zoom of damage growing in the local model.
|}

On the other side, in Figures [[#img-43a|43a]], [[#img-43b|43b]] and [[#img-43c|43c]] we see a very similar diffusive damaged area. Indeed, in the non-local approach, damage size is controlled by the interaction radius <math display="inline">R</math> and thus, even when reducing the size of the elements, the damaged area is practically unaltered.

<div id='img-43a'></div>
<div id='img-43b'></div>
<div id='img-43c'></div>
<div id='img-43'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_nloc_m1_d.png|480px|Mesh 1.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1.
|-
|[[Image:draft_Content_981184740-tpbt_nloc_m2_d.png|480px|Mesh 2.]]
|- style="text-align: center; font-size: 75%;"
| (b) Mesh 2.
|-
|[[Image:draft_Content_981184740-tpbt_nloc_m3_d.png|480px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 3.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 43:''' TPB test. Damage initiation in the non-local model.
|}

This concept of damage localization and diffusion can also be seen in the evolution of damage obtained from the results of mesh 3 (Figure [[#img-44|44]]).

<div id='img-44a'></div>
<div id='img-44b'></div>
<div id='img-44c'></div>
<div id='img-44d'></div>
<div id='img-44e'></div>
<div id='img-44f'></div>
<div id='img-44'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_loc_d_1.png|400px|Local model. Initial stage.]]
|[[Image:draft_Content_981184740-tpbt_nloc_d_1.png|400px|Non-local model. Initial stage.]]
|- style="text-align: center; font-size: 75%;"
| (a) Local model. Initial stage.
| (b) Non-local model. Initial stage.
|-
|[[Image:draft_Content_981184740-tpbt_loc_d_2.png|400px|Local model. Intermediate stage.]]
|[[Image:draft_Content_981184740-tpbt_nloc_d_2.png|400px|Non-local model. Intermediate stage.]]
|- style="text-align: center; font-size: 75%;"
| (c) Local model. Intermediate stage.
| (d) Non-local model. Intermediate stage.
|-
|[[Image:draft_Content_981184740-tpbt_loc_d_3.png|400px|Local model. Advanced stage.]]
|[[Image:draft_Content_981184740-tpbt_nloc_d_3.png|400px|Non-local model. Advanced stage.]]
|- style="text-align: center; font-size: 75%;"
| (e) Local model. Advanced stage.
| (f) Non-local model. Advanced stage.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 44:''' TPB test. Evolution of damage propagation.
|}

The mesh-induced directional bias and mesh sensitivity observed in the local damage model are practically absent in the non-local model, and so the latter seems to be a more robust approach. However, apart from the mesh objectivity, it is important to assess whether the implemented model can properly reproduce the behaviour of the real beam. In this regard, we have compared the results obtained in the non-local model with experiments performed in <span id='citeF-62'></span>[[#cite-62|[62]]].

Thereby, Figure [[#img-45|45]] shows the force-vertical deflection diagram of the reference solution and the computed solution with mesh 3.

<div id='img-45'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 75%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_ref.png|450px|TPB test. Non-local model validation.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 45:''' TPB test. Non-local model validation.
|}

As we can see, the peak load is properly captured, but the post peak branch of the numerical solution falls faster than in the reference solution, and even shows a certain amount of snap-back behaviour. Moreover, looking at the elastic branch of the responses, it seems that the stiffness degradation starts before in the computed solution and, in consequence, the peak is slightly displaced to the right. Therefore, we can say that the behaviour of the numerical model seems more brittle than the observed in experiments.

===4.2.2 Mesh-adaptive example===

Finally, we are going to solve the same problem as before with the non-local damage model, but this time we will apply the implemented mesh-adaptive technique.

As we have commented, this adaptive procedure is a recent implementation still under development. As such, there are certain tools that are not currently available, but that need to be considered in the future. One of them is the possibility of measuring forces and displacements at points of interest in an adaptive mesh. Indeed, every time mesh is refined, nodes and elements numbering is altered, and so we cannot measure the forces and displacement like we did in the previous cases with a fixed mesh. Another limitation is that the current adaptive procedure has not any global indicator to initiate an error estimation and subsequent mesh refinement. This means that, before the problem starts, we must impose a certain number of refinements that will be performed uniformly distributed throughout all the incremental-iterative process.

Taking this into account, the main purpose of this subsection is simply to review an example of an adaptive procedure so as to understand the advantages and disadvantages of using this technique instead of a fixed mesh like before.

As commented in section [[#3.5 Mesh-Adaptive Technique|3.5]], we will be using a remeshing criterion based on global error equidistribution. Moreover, we have considered a permissible percentage of global error of <math display="inline">\eta=10%</math>, we have imposed three refinements, and the initial mesh is a uniform unstructured mesh of quadrilaterals with an average size of 10 mm. The resulting spatial discretizations throughout the process are represented in Figure [[#img-46|46]].

<div id='img-46a'></div>
<div id='img-46b'></div>
<div id='img-46c'></div>
<div id='img-46d'></div>
<div id='img-46'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_amr_m1.png|400px|Mesh 0: 4453 elements, 4618 nodes.]]
|[[Image:draft_Content_981184740-tpbt_amr_m2.png|400px|Mesh 1: 2977 elements, 3146 nodes.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 0: 4453 elements, 4618 nodes.
| (b) Mesh 1: 2977 elements, 3146 nodes.
|-
|[[Image:draft_Content_981184740-tpbt_amr_m3.png|400px|Mesh 2: 1598 elements, 1731 nodes.]]
|[[Image:draft_Content_981184740-tpbt_amr_m4.png|400px|Mesh 3: 1627 elements, 1759 nodes.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 2: 1598 elements, 1731 nodes.
| (d) Mesh 3: 1627 elements, 1759 nodes.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 46:''' TPB test. Evolution of the adaptive mesh.
|}

In order to understand the adaptive mesh process, it is interesting to see the evolution of damage for the different spatial discretizations (Figure [[#img-47|47]]), as well as the evolution of the refinement parameter of the elements <math display="inline">\xi ^{(e)}</math> (Figure [[#img-48|48]]).

<div id='img-47a'></div>
<div id='img-47b'></div>
<div id='img-47c'></div>
<div id='img-47d'></div>
<div id='img-47'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_amr_d_1.png|480px|Mesh 0.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 0.
|-
|[[Image:draft_Content_981184740-tpbt_amr_d_2.png|480px|Mesh 1.]]
|- style="text-align: center; font-size: 75%;"
| (b) Mesh 1.
|-
|[[Image:draft_Content_981184740-tpbt_amr_d_3.png|480px|Mesh 2.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 2.
|-
|[[Image:draft_Content_981184740-tpbt_amr_d_4.png|480px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (d) Mesh 3.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 47:''' TPB test. Progression of damage in the adaptive mesh.
|}

Let us remember that the adaptive procedure estimates the error of the mesh from the effective stresses in the structure. Thereby, the parts of the structure with higher stress variations will tend to be refined, whereas those zones with smoother stresses will tend to have larger elements. This explains the concentration of elements at the zone of damage progression (Figure [[#img-47|47]]) and at the points of support and load application. Furthermore, the refinement parameter shown in Figure [[#img-48|48]] represents the ratio between the actual element error and the desired error. Thereby, those elements depicted in blue (<math display="inline">\xi ^{(e)}<1</math>) are elements that could be larger, while the elements in red (<math display="inline">\xi ^{(e)}>1</math>) should be refined in order to reduce its error. Elements in green (<math display="inline">\xi ^{(e)}=1</math>) would be of optimal size.

<div id='img-48a'></div>
<div id='img-48b'></div>
<div id='img-48c'></div>
<div id='img-48d'></div>
<div id='img-48'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-tpbt_amr_ref_1.png|480px|Mesh 0.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 0.
|-
|[[Image:draft_Content_981184740-tpbt_amr_ref_2.png|480px|Mesh 1.]]
|- style="text-align: center; font-size: 75%;"
| (b) Mesh 1.
|-
|[[Image:draft_Content_981184740-tpbt_amr_ref_3.png|480px|Mesh 2.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 2.
|-
|[[Image:draft_Content_981184740-tpbt_amr_ref_4.png|480px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (d) Mesh 3.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 48:''' TPB test. Progression of the refinement parameter in the adaptive mesh.
|}

Looking at the number of elements in Figure [[#img-46|46]], we can say that the main advantage of an adaptive mesh technique is that it automatically optimizes the mesh for the problem we are solving, so that the number of elements can be drastically reduced, and so the computational cost of the problem.

Nevertheless, this procedure has also shown some inconveniences when solving damage problems with a non-local approach. For instance, if damage starts propagating abruptly, the refinement process can lead to a considerable large amount of elements in the the damage zone. This concentration of elements can noticeably increase the cost of the non-local interactions and so the computation can become very intensive. However, the good point is that once damage propagation stabilizes, the adaptive procedure readjusts automatically the mesh, and so one ends with an actually efficient spatial discretization.

An additional limitation is that the procedure does not differentiate compressive stresses from tensile stresses when estimating the error. For materials like concrete this can lead to excessively small elements in compressed zones, and so it could be convenient to give a heavier weight to the tensile stresses.

Finally, one should also notice that the adaptive mesh process will always have some error associated to the transfer of internal variables between meshes. Thereby, once we can measure forces and displacements in adaptive mesh, we should check whether the response obtained with a fixed mesh is reproduced accurately in an adaptive process.

==4.3 Single-Edge Notched Beam Test==

In this second example a single-edge notched beam (SENB) is subjected to non-symmetrical four-point bending. The analysed beam has a square cross section of 100 <math display="inline">\times </math> 200 mm, a span of 840 mm, and the notch is 10 mm thick and 40 mm depth (see Figure [[#img-49|49]]).

<div id='img-49'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-single-edge_notched_beam_test.png|480px|SENB test. Problem statement. Distances in mm]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 49:''' SENB test. Problem statement. Distances in mm
|}

Like in the previous case, in order to reduce fictitious stress concentrations at the supports and at the loading point, both the displacement constraints and the loads are acting over a line of 20 mm (see Figure [[#img-50|50]]).

<div id='img-50a'></div>
<div id='img-50b'></div>
<div id='img-50'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-displacement_bc_senbt.png|480px|Displacement constraints.]]
|- style="text-align: center; font-size: 75%;"
| (a) Displacement constraints.
|-
|[[Image:draft_Content_981184740-loads_senbt.png|480px|Loads.]]
|- style="text-align: center; font-size: 75%;"
| (b) Loads.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 50:''' SENB test. Boundary Conditions.
|}

Once again, we have assumed plane stress conditions, but in this case the geometry will be meshed by means of linear three-node triangular elements with 1 integration point.

===4.3.1 Damage models testing===

As we did for the three-point bending test, we will first solve the problem using the two implemented damage models. The local approach will be modelled with the same modified Simo and Ju model of the previous example, using the material parameters shown in Table [[#table-7|7]]. On the other hand, in this example the non-local damage model will be defined with the equivalent strain form of the modified von Mises model ([[#eq-26|26]]), and with the exponential damage evolution law that we presented in ([[#eq-36|36]]). The material parameters for the non-local approach  have been obtained from <span id='citeF-63'></span>[[#cite-63|[63]]], and are summarized in Table [[#table-8|8]].


{| class="wikitable" style="text-align: left; margin: 1em auto;"
|+ <span id='table-7'></span>Table. 7 SENB test. Material parameters for the Simo-Ju local model.
|-
|  '''Parameter''' 
| '''Value''' 
|-
|  Young's modulus (<math display="inline">E</math>) 
| 28000 MPa 
|-
| Poisson's ratio (<math display="inline">\nu </math>) 
| 0.1 
|-
| Compressive strength (<math display="inline">\sigma _y^c</math>) 
| 35 MPa 
|-
| Tensile strength (<math display="inline">\sigma _y^t</math>) 
| 3.2 MPa 
|-
| Fracture energy (<math display="inline">G_f</math>) 
| 140 <math display="inline">J/m^2</math> 
|-
| Limit fracture length (<math display="inline">l_{lim}</math>) 
| 5 mm 

|}


{| class="wikitable" style="text-align: left; margin: 1em auto;"
|+ <span id='table-8'></span>Table. 8 SENB test. Material parameters for the modified von Mises non-local model.
|-
|  '''Parameter''' 
| '''Value''' 
|-
|  Young's modulus (<math display="inline">E</math>) 
| 28000 MPa 
|-
| Poisson's ratio (<math display="inline">\nu </math>) 
| 0.1 
|-
| Damage threshold (<math display="inline">r_0</math>) 
| <math display="inline">1.5\cdot{10}^{-4}</math> 
|-
| Compressive to tensile strength ratio (<math display="inline">\kappa </math>) 
| 10 
|-
| Parameter <math display="inline">A</math> (<math display="inline">A</math>) 
| 0.8 
|-
| Parameter <math display="inline">B</math> (<math display="inline">B</math>) 
| 9000 
|-
| Characteristic length (<math display="inline">l_c</math>) 
| 10 mm 

|}

Just like before, we are going to solve the problem for different spatial discretizations. Thereby, we will be using three unstructured meshes of triangles with a minimum size of 8 mm, 5 mm, and 3 mm, respectively (see Figure [[#img-51|51]]).

<div id='img-51a'></div>
<div id='img-51b'></div>
<div id='img-51c'></div>
<div id='img-51'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_m1.png|400px|Mesh 1: 2216 elements, 1264 nodes.]]
|[[Image:draft_Content_981184740-senbt_m2.png|400px|Mesh 2: 3502 elements, 1923 nodes.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1: 2216 elements, 1264 nodes.
| (b) Mesh 2: 3502 elements, 1923 nodes.
|-
| colspan="2"|[[Image:draft_Content_981184740-senbt_m3.png|300px|Mesh 3: 7183 elements, 3796 nodes.]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (c) Mesh 3: 7183 elements, 3796 nodes.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 51:''' SENB test. Spatial discretizations.
|}

<div id='img-52a'></div>
<div id='img-52b'></div>
<div id='img-52'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 75%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_local.png|450px|Partially regularized local damage model.]]
|- style="text-align: center; font-size: 75%;"
| (a) Partially regularized local damage model.
|-
|[[Image:draft_Content_981184740-senbt_nonlocal.png|450px|Non-local damage model.]]
|- style="text-align: center; font-size: 75%;"
| (b) Non-local damage model.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 52:''' SENB test. Force-Crack Mouth Sliding Displacement curves.
|}

Figure [[#img-52|52]] shows the relation between the applied load and the Crack Mouth Sliding Displacement (CMSD). If we look at the curves in Figure [[#img-52a|52a]] we can see that in this case the peak and the dissipated energy also decrease as the mesh is refinement. However, this reduction is not so clear as in the previous example. On the other hand, the response diagram in Figure [[#img-52b|52b]] shows a peak load practically identical for the three meshes, but there are some differences in the response near the post-peak region.

Thereby, the response diagrams of the single-edge notched beam test seem to be in agreement with the ones obtained for the three-point bending test, although in this case the difference between the local and non-local approaches is not so clear. Now, in order to properly analyse the behaviour of both models, let us present some snapshots with the damage distribution.

Figures [[#img-53|53]] and [[#img-54|54]] show the damage progression in the peak region of the response, for the local and non-local model, respectively. Like before, from the size and shape of the damaged zone, we can clearly see a more localized response in the local damage model, and a more diffusive response in the non-local case.

Looking at the different damage patterns of the local model (Figure [[#img-53|53]]), we notice an additional vertical damage line that appears in the coarser meshes, stressing the idea that the local model suffers from mesh sensitivity.

<div id='img-53a'></div>
<div id='img-53b'></div>
<div id='img-53c'></div>
<div id='img-53'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 90%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_loc_m1_d.png|540px|Mesh 1.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1.
|-
|[[Image:draft_Content_981184740-senbt_loc_m2_d.png|540px|Mesh 2.]]
|- style="text-align: center; font-size: 75%;"
| (b) Mesh 2.
|-
|[[Image:draft_Content_981184740-senbt_loc_m3_d.png|540px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 3.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 53:''' SENB test. Damage progression in the local model.
|}

Regarding the non-local model (Figure [[#img-54|54]]), if we compare this example with the three-point bending test, we can clearly see that the progression zone here is restricted to a narrower region than before. The reason is that now the characteristic length defining the interaction radius is quite smaller than before.

<div id='img-54a'></div>
<div id='img-54b'></div>
<div id='img-54c'></div>
<div id='img-54'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_nloc_m1_d.png|480px|Mesh 1.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1.
|-
|[[Image:draft_Content_981184740-senbt_nloc_m2_d.png|480px|Mesh 2.]]
|- style="text-align: center; font-size: 75%;"
| (b) Mesh 2.
|-
|[[Image:draft_Content_981184740-senbt_nloc_m3_d.png|480px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 3.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 54:''' SENB test. Damage progression in the non-local model.
|}

Focusing on the results obtained with mesh 3, the comparative evolution of damage for both models is represented in Figure [[#img-55|55]].

<div id='img-55a'></div>
<div id='img-55b'></div>
<div id='img-55c'></div>
<div id='img-55d'></div>
<div id='img-55e'></div>
<div id='img-55f'></div>
<div id='img-55'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_loc_d_1.png|400px|Local model. Initial stage.]]
|[[Image:draft_Content_981184740-senbt_nloc_d_1.png|400px|Non-local model. Initial stage.]]
|- style="text-align: center; font-size: 75%;"
| (a) Local model. Initial stage.
| (b) Non-local model. Initial stage.
|-
|[[Image:draft_Content_981184740-senbt_loc_d_2.png|400px|Local model. Intermediate stage.]]
|[[Image:draft_Content_981184740-senbt_nloc_d_2.png|400px|Non-local model. Intermediate stage.]]
|- style="text-align: center; font-size: 75%;"
| (c) Local model. Intermediate stage.
| (d) Non-local model. Intermediate stage.
|-
|[[Image:draft_Content_981184740-senbt_loc_d_3.png|400px|Local model. Advanced stage.]]
|[[Image:draft_Content_981184740-senbt_nloc_d_3.png|400px|Non-local model. Advanced stage.]]
|- style="text-align: center; font-size: 75%;"
| (e) Local model. Advanced stage.
| (f) Non-local model. Advanced stage.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 55:''' SENB test. Evolution of damage propagation.
|}

We can see that, at the beginning of damage propagation (Figures [[#img-55a|55a]] and [[#img-55b|55b]]) both approaches show a very similar damage zone. This more localized damage zone in the non-local model is mainly due to the reduced interaction radius of this example. After that, the difference between both models is noticeable.

Like in the previous case, we have compared the results obtained in the non-local model with an experimental solution <span id='citeF-64'></span>[[#cite-64|[64]]].

Figure [[#img-56|56]] shows the relation between the applied load and the Crack Mouth Sliding Displacement of the reference solution and the computed solution with mesh 3.

<div id='img-56'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 75%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_ref.png|450px|SENB test. Non-local model validation.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 56:''' SENB test. Non-local model validation.
|}

Thereby, we can see that the obtained response is very similar to the reference one. What we should notice, though, is that the failure obtained in the computed solution seems more brittle as compared to the reference one. This is possibly consequence of an imprecise tracing of the equilibrium path, caused by the global arc-length strategy used. As commented before, there are specific control parameters for the arc length method that are probably more suited for this kind of quasi-brittle behaviour.

===4.3.2 Mesh-adaptive example===

Finally, we are going to solve the problem again, applying the implemented mesh-adaptive technique on the non-local damage model.

Like in the three-point bending test, we will be using the remeshing criterion based on global error equidistribution, with a permissible percentage of global error of <math display="inline">\eta=10%</math>. In this case we have imposed four refinements, and the initial mesh is a uniform unstructured mesh of triangles with an average size of 7 mm. The resulting spatial discretizations throughout the process are represented in Figure [[#img-57|57]].

<div id='img-57a'></div>
<div id='img-57b'></div>
<div id='img-57c'></div>
<div id='img-57d'></div>
<div id='img-57e'></div>
<div id='img-57'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_amr_m0.png|400px|Mesh 0: 7784 elements, 4048 nodes.]]
|[[Image:draft_Content_981184740-senbt_amr_m1.png|400px|Mesh 1: 6800 elements, 3624 nodes.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 0: 7784 elements, 4048 nodes.
| (b) Mesh 1: 6800 elements, 3624 nodes.
|-
|[[Image:draft_Content_981184740-senbt_amr_m2.png|400px|Mesh 2: 6590 elements, 3521 nodes.]]
|[[Image:draft_Content_981184740-senbt_amr_m3.png|400px|Mesh 3: 10626 elements, 5563 nodes.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 2: 6590 elements, 3521 nodes.
| (d) Mesh 3: 10626 elements, 5563 nodes.
|-
| colspan="2"|[[Image:draft_Content_981184740-senbt_amr_m4.png|300px|Mesh 4: 2542 elements, 1388 nodes.]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (e) Mesh 4: 2542 elements, 1388 nodes.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 57:''' SENB test. Evolution of the adaptive mesh.
|}

Just like before, we are going to review the evolution of damage for the different spatial discretizations (Figure [[#img-58|58]]), as well as the evolution of the refinement parameter of the elements (Figure [[#img-59|59]]).

Looking at Figure [[#img-58a|58a]], we see that there is no damage in the original mesh. Moreover, if we look also at Figure [[#img-59a|59a]], it can be deduced that the subsequent spatial discretization is mainly caused by the high compressive stresses at the loading points and at the supports. Thereby, as we stated in the previous example, it could be very interesting to modify the evaluated stresses during the error estimation, in order to account for the different resistance in tension and compression of some geo-materials like concrete.

However, from Figures [[#img-57c|57c]], [[#img-57d|57d]] and [[#img-57e|57e]] it seems clear that the adaptive procedure properly captures the damage process zone and adapts the mesh accordingly.

<div id='img-58a'></div>
<div id='img-58b'></div>
<div id='img-58c'></div>
<div id='img-58d'></div>
<div id='img-58e'></div>
<div id='img-58'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_amr_d_0.png|480px|Mesh 0.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 0.
|-
|[[Image:draft_Content_981184740-senbt_amr_d_1.png|480px|Mesh 1.]]
|- style="text-align: center; font-size: 75%;"
| (b) Mesh 1.
|-
|[[Image:draft_Content_981184740-senbt_amr_d_2.png|480px|Mesh 2.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 2.
|-
|[[Image:draft_Content_981184740-senbt_amr_d_3.png|480px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (d) Mesh 3.
|-
|[[Image:draft_Content_981184740-senbt_amr_d_4.png|480px|Mesh 4.]]
|- style="text-align: center; font-size: 75%;"
| (e) Mesh 4.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 58:''' SENB test. Progression of damage in the adaptive mesh.
|}

Another aspect to consider here is directly related to the abrupt failure of the computed solution seen in Figure [[#img-56|56]]. Indeed, as explained in the previous example, when damage grows abruptly, the adaptive technique refines excessively all the affected zone, as it can be seen from the predominant blue color in Figure [[#img-59d|59d]]. Nonetheless, once damage stabilises the adaptive procedure readjusts the spatial discretization and the resulting mesh [[#img-57e|57e]] is really efficient.

Finally, we can notice that the final damage pattern in Figure [[#img-58e|58e]] is a bit wider than the one shown in [[#img-55f|55f]]. The reason could be that   some extra damage diffusion is introduced due to the errors appearing during the transfer of internal variables between meshes. Therefore, it would be interesting to obtain the force-CMSD curve when using this adaptive technique, in order to assess the effect of these errors in the traced equilibrium path.

<div id='img-59a'></div>
<div id='img-59b'></div>
<div id='img-59c'></div>
<div id='img-59d'></div>
<div id='img-59e'></div>
<div id='img-59'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"
|-
|[[Image:draft_Content_981184740-senbt_amr_ref_0.png|480px|Mesh 0.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 0.
|-
|[[Image:draft_Content_981184740-senbt_amr_ref_1.png|480px|Mesh 1.]]
|- style="text-align: center; font-size: 75%;"
| (b) Mesh 1.
|-
|[[Image:draft_Content_981184740-senbt_amr_ref_2.png|480px|Mesh 2.]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 2.
|-
|[[Image:draft_Content_981184740-senbt_amr_ref_3.png|480px|Mesh 3.]]
|- style="text-align: center; font-size: 75%;"
| (d) Mesh 3.
|-
|[[Image:draft_Content_981184740-senbt_amr_ref_4.png|480px|Mesh 4.]]
|- style="text-align: center; font-size: 75%;"
| (e) Mesh 4.
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 59:''' SENB test. Progression of the refinement parameter in the adaptive mesh.
|}

=5 Conclusions and Future Research Lines=

We have been working in damage mechanics for almost half a year in the seek of a robust numerical method for the modelling of quasi-brittle materials failure, and we can point out the following conclusions.

First, we have seen that the partially regularization of the local damage model still shows mesh sensitiveness and presents severe mesh-induced directional bias, and so this model is not the best approach if one aims at achieving robust results. However, if the model is properly calibrated, it can become an interesting tool to obtain a first estimation of the damage propagation with a low computational cost.

On the other hand, the fully regularized non-local damage approach removes pathological mesh sensitivity and has proved to properly capture the peak load of the response diagrams in the three-point bending test, and in the single-edge notched beam test. In the former test, however, it has shown a more brittle post-peak branch as compared to the experimental results. Furthermore, an important aspect to note about this integral-type non-local approach is that the averaging performed to account for the non-local interaction, considerably increases the computational cost of the solution, as compared to the classical local approach, and also modifies the traditional assembly process of the global tangent stiffness matrix.

Concerning the non-linearity associated to damage, it must be stated that the implemented arc-length strategy has shown an irregular success when tracing the non-linear solution, and so it is probably not the best approach for the kind of problems we have been solving.

The implemented mesh-adaptive technique allows capturing the process damage zone and adjust the spatial discretization accordingly. This procedure generates large elements in zones with smooth stress distribution, and small elements in zones with strong variations of stresses, and so efficient mesh distributions can be obtained. However, the adaptive technique is still a recent implementation and there are some features that could be improved. For instance, the procedure does not differentiate compressive stresses from tensile stresses when estimating the error. This can lead to excessively fine meshes at compressed zones that are not actually transcendent for the damage progression of materials like concrete. Furthermore, errors in the transfer operations of internal variables introduce some extra damage diffusion that should be analysed in order to assess its influence in the traced equilibrium path.

Thereby, in the near future there are various aspects on the current code that should be enhanced. First, regarding the convergence of the non-linear solution, we could implement an advanced arc-length method with a specific control parameter (CMOD or CMSD) that accounts for the localized nature of quasi-brittle materials.

The measurement of forces and displacements in an adaptive mesh is also another important tool that would enable analysing the effect of the errors introduced into the response diagram of a structure, when transferring internal variables between different meshes.

Still in the mesh-adaptive technique, we could introduce some weights affecting the stresses evaluated in the error estimation process, in order to account for the different behaviour in tension and compression of many geo-materials. This way, even more optimized spatial discretizations could be generated, and the exaggerated concentrations of elements seen at the loading points and supports could be drastically reduced.

One last additional feature that should be implemented in the mesh-adaptive process, is a global error parameter that could be used to determine, at any load step, whether the error estimation process must be initialized or not.

Finally, as other research lines of interest, we should mention the extension of the implemented damage model and mesh-adaptive procedure to three-dimensional analysis; the coupling of damage and plasticity, necessary if one wants to simulate the interaction between concrete and steel; and the coupling of damage and fracture mechanics, which could make use of the information of the damage propagation to obtain the position and direction of an explicit fracture.

&nbsp;

===References===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' Bazant, Zdenek P and Cedolin, Luigi. ''Blunt crack band propagation in finite element analysis'', Volume 105. ASCE. Journal of the Engineering Mechanics Division 2 297&#8211;315, 1979

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' Dragon, A and Mroz, Z. ''A continuum model for plastic-brittle behaviour of rock and concrete'', Volume 17. Elsevier. International Journal of Engineering Science 2 121&#8211;137, 1979

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]''' Di Prisco, M and Mazars, J. ''Crush-crack': a non-local damage model for concrete'', Volume 1. Mechanics of Cohesive-frictional Materials 4 321&#8211;347, 1996

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]''' Willam, KJ and Warnke, EP. ''Constitutive model for the triaxial behavior of concrete'', Volume 19. ISMES, Bergamo, Italy 174, 1975

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]''' De Borst, René. ''Non-linear analysis of frictional materials'', 1986

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' Lemaitre, Jean. ''Coupled elasto-plasticity and damage constitutive equations'', Volume 51. Elsevier. Computer Methods in Applied Mechanics and Engineering 1 31&#8211;49, 1985

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' Klisiski, M and Mroz, Z. ''Description of inelastic deformation and degradation of concrete'', Volume 24. Elsevier. International journal of solids and structures 4 391&#8211;416, 1988

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]''' Bazant, Zdenek P and Ozbolt, Josko. ''Nonlocal microplane model for fracture, damage, and size effect in structures'', Volume 116. American Society of Civil Engineers. Journal of Engineering Mechanics 11 2485&#8211;2505, 1990

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]''' Jirásek, Milan. ''Numerical modeling of deformation and failure of materials''. Lecture notes, 1999

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]''' Herrmann, HJ and Kertész, J and De Arcangelis, L. ''Fractal shapes of deterministic cracks'', Volume 10. IOP Publishing. EPL (Europhysics Letters) 2 147, 1989

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]''' Vervuurt, Adri and Van Mier, Jan GM and Schlangen, Erik. ''Analyses of anchor pull-out in concrete'', Volume 27. Springer. Materials and Structures 5 251&#8211;259, 1994

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]''' Sakaguchi, H and Mühlhaus, HB. ''Mesh free modelling of failure and localization in brittle materials''. Pergamon: Oxford. Deformation and Progressive Failure in Geomechanics 15&#8211;21, 1997

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]''' Recarey Morfa, Carlos A and Oñate, E and Rojek, Jerzy and Canet, Juan Miquel and Zarate, Francisco and Cadoce, Gilles and Labra, Carlos A. ''Modelación del ensayo Brasileño empleando modelación discreta''. Ingeniería Civil 139, 2005

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]''' Hult, J. ''Creep in continua and structures''. Springer 137&#8211;155, 1974

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]''' Kachanov, LM. ''On the time of fracture under conditions of creep''. Izv. AN SSSR, Otd. Tekh. Nauk 8 26&#8211;35, 1958

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]''' Rabotnov, YN. ''Damage from creep'', Volume 2. Zhurn. Prikl. Mekh. Tekhn. Phys 113&#8211;123, 1963

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]''' Hayhurst, DR. ''Creep rupture under multi-axial states of stress'', Volume 20. Elsevier. Journal of the Mechanics and Physics of Solids 6 381&#8211;382, 1972

<div id="cite-18"></div>
'''[[#citeF-18|[18]]]''' Leckie, FA and Hayhurst, DR. ''Creep rupture of structures'', Volume 340. The Royal Society. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 1622 323&#8211;347, 1974

<div id="cite-19"></div>
'''[[#citeF-19|[19]]]''' Lemaitre, J and Chaboche, JL. ''Mécanique des matériaux solides, 1985''. Dunod, Paris

<div id="cite-20"></div>
'''[[#citeF-20|[20]]]''' Mazars, Jacky. ''A description of micro-and macroscale damage of concrete structures'', Volume 25. Elsevier. Engineering Fracture Mechanics 5 729&#8211;737, 1986

<div id="cite-21"></div>
'''[[#citeF-21|[21]]]''' Mazars, Jacky and Pijaudier-Cabot, Gilles. ''From damage to fracture mechanics and conversely: a combined approach'', Volume 33. Elsevier. International Journal of Solids and Structures 20 3327&#8211;3342, 1996

<div id="cite-22"></div>
'''[[#citeF-22|[22]]]''' Chaboche, J-L. ''The concept of effective strain applied to elasticity and viscoelasticity in the presence of anisotropic damage'' 1979-77, 1979

<div id="cite-23"></div>
'''[[#citeF-23|[23]]]''' Simo, JC and Ju, JW. ''Strain-and stress-based continuum damage models-I. Formulation'', Volume 23. Elsevier. International journal of solids and structures 7 821&#8211;840, 1987

<div id="cite-24"></div>
'''[[#citeF-24|[24]]]''' Ju, JW. ''Energy-based coupled elastoplastic damage models at finite strains'', Volume 115. American Society of Civil Engineers. Journal of Engineering Mechanics 11 2507&#8211;2525, 1989

<div id="cite-25"></div>
'''[[#citeF-25|[25]]]''' Ju, JW. ''On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects'', Volume 25. Elsevier. International Journal of Solids and structures 7 803&#8211;833, 1989

<div id="cite-26"></div>
'''[[#citeF-26|[26]]]''' Oller, S and Oñate, E and Oliver, J and Lubliner, J. ''Finite element nonlinear analysis of concrete structures using a “plastic-damage model”'', Volume 35. Elsevier. Engineering Fracture Mechanics 1 219&#8211;231, 1990

<div id="cite-27"></div>
'''[[#citeF-27|[27]]]''' Oliver, J and Cervera, M and Oller, S and Lubliner, J. ''Isotropic damage models and smeared crack analysis of concrete'', Volume 2 945&#8211;958, 1990

<div id="cite-28"></div>
'''[[#citeF-28|[28]]]''' Oliver, J and Cervera, M and Oller, S and Lubliner, J. ''A simple damage model for concrete, including long term effects'', Volume 2 945&#8211;958, 1990

<div id="cite-29"></div>
'''[[#citeF-29|[29]]]''' Cervera, M and Chiumenti, M and de Saracibar, C Agelet. ''Shear band localization via local J 2 continuum damage mechanics'', Volume 193. Elsevier. Computer Methods in Applied Mechanics and Engineering 9 849&#8211;880, 2004

<div id="cite-30"></div>
'''[[#citeF-30|[30]]]''' Cervera, M and Chiumenti, M. ''Mesh objective tensile cracking via a local continuum damage model and a crack tracking technique'', Volume 196. Elsevier. Computer Methods in Applied Mechanics and Engineering 1 304&#8211;320, 2006

<div id="cite-31"></div>
'''[[#citeF-31|[31]]]''' Mazars, Jacky. ''Application de la mécanique de l'endommagement au comportement non linéaire et a la rupture du béton de structure'', 1984

<div id="cite-32"></div>
'''[[#citeF-32|[32]]]''' De Vree, JHP and Brekelmans, WAM and Van Gils, MAJ. ''Comparison of nonlocal approaches in continuum damage mechanics'', Volume 55. Elsevier. Computers & Structures 4 581&#8211;588, 1995

<div id="cite-33"></div>
'''[[#citeF-33|[33]]]''' Pijaudier-Cabot, Gilles and Huerta, Antonio. ''Finite element analysis of bifurcation in nonlocal strain softening solids'', Volume 90. Elsevier. Computer methods in applied mechanics and engineering 1 905&#8211;919, 1991

<div id="cite-34"></div>
'''[[#citeF-34|[34]]]''' Mühlhaus, HB. ''SHEARBAND ANALYSIS IN ANTIGRANULOCYTES MATERIAL BY COSSERAT-THEORY'', Volume 56. SPRINGER VERLAG 175 FIFTH AVE, NEW YORK, NY 10010. Ingenieur Archiv 5 389&#8211;399, 1986

<div id="cite-35"></div>
'''[[#citeF-35|[35]]]''' Tejchman, J and Wu, W. ''Numerical study on patterning of shear bands in a Cosserat continuum'', Volume 99. Springer. Acta Mechanica 1-4 61&#8211;74, 1993

<div id="cite-36"></div>
'''[[#citeF-36|[36]]]''' Zbib, HM and Aifantis, EC. ''A gradient-dependent flow theory of plasticity: application to metal and soil instabilities'', Volume 42. American Society of Mechanical Engineers. Applied Mechanics Reviews 11S S295&#8211;S304, 1989

<div id="cite-37"></div>
'''[[#citeF-37|[37]]]''' Mühlhaus, HB and Alfantis, EC. ''A variational principle for gradient plasticity'', Volume 28. Elsevier. International Journal of Solids and Structures 7 845&#8211;857, 1991

<div id="cite-38"></div>
'''[[#citeF-38|[38]]]''' Meftah, F and Reynouard, JM. ''A multilayered mixed beam element in gradient plasticity for the analysis of localized failure modes'', Volume 3. Wiley Online Library. Mechanics of Cohesive-frictional Materials 4 305&#8211;322, 1998

<div id="cite-39"></div>
'''[[#citeF-39|[39]]]''' Loret, Benjamin and Prevost, Jean H. ''Dynamic strain localization in elasto-(visco-) plastic solids, Part 1. General formulation and one-dimensional examples'', Volume 83. Elsevier. Computer Methods in Applied Mechanics and Engineering 3 247&#8211;273, 1990

<div id="cite-40"></div>
'''[[#citeF-40|[40]]]''' Prevost, Jean H and Loret, Benjamin. ''Dynamic strain localization in elasto-(visco-) plastic solids, part 2. Plane strain examples'', Volume 83. Elsevier. Computer Methods in Applied Mechanics and Engineering 3 275&#8211;294, 1990

<div id="cite-41"></div>
'''[[#citeF-41|[41]]]''' Ehlers, W and Graf, T. ''Adaptive computation of localization phenomena in geotechnical applications''. Bifurcations and Instabilities in Geomechanics, Swets, Zeitlinger 247&#8211;262, 2003

<div id="cite-42"></div>
'''[[#citeF-42|[42]]]''' Bazant, Zdeněk P. ''Mechanics of distributed cracking'', Volume 39. American Society of Mechanical Engineers. Applied Mechanics Reviews 5 675&#8211;705, 1986

<div id="cite-43"></div>
'''[[#citeF-43|[43]]]''' Pijaudier-Cabot, Gilles and Bazant, Zdenek P. ''Nonlocal damage theory'', Volume 113. American Society of Civil Engineers. Journal of Engineering Mechanics 10 1512&#8211;1533, 1987

<div id="cite-44"></div>
'''[[#citeF-44|[44]]]''' Jirásek, Milan and Rolshoven, Simon. ''Comparison of integral-type nonlocal plasticity models for strain-softening materials'', Volume 41. Elsevier. International Journal of Engineering Science 13 1553&#8211;1602, 2003

<div id="cite-45"></div>
'''[[#citeF-45|[45]]]''' Rodríguez-Ferran, Antonio and Morata, Irene and Huerta, Antonio. ''A new damage model based on non-local displacements'', Volume 29. Wiley Online Library. International Journal for Numerical and Analytical Methods in Geomechanics 5 473&#8211;493, 2005

<div id="cite-46"></div>
'''[[#citeF-46|[46]]]''' Bazant, Zdenek P and Oh, Byung H. ''Crack band theory for fracture of concrete'', Volume 16. Springer. Matériaux et construction 3 155&#8211;177, 1983

<div id="cite-47"></div>
'''[[#citeF-47|[47]]]''' Pietruszczak, ST and Mroz, Z. ''Finite element analysis of deformation of strain-softening materials'', Volume 17. Wiley Online Library. International Journal for Numerical Methods in Engineering 3 327&#8211;334, 1981

<div id="cite-48"></div>
'''[[#citeF-48|[48]]]''' Eringen, A Cemal. ''On nonlocal plasticity'', Volume 19. Elsevier. International Journal of Engineering Science 12 1461&#8211;1474, 1981

<div id="cite-49"></div>
'''[[#citeF-49|[49]]]''' Eringen, A Cemal. ''On differential equations of nonlocal elasticity and solutions of screw dislocation and surface waves'', Volume 54. AIP Publishing. Journal of Applied Physics 9 4703&#8211;4710, 1983

<div id="cite-50"></div>
'''[[#citeF-50|[50]]]''' Bazant, Zdenek P and Jirásek, Milan. ''Nonlocal integral formulations of plasticity and damage: survey of progress'', Volume 128. American Society of Civil Engineers. Journal of Engineering Mechanics 11 1119&#8211;1149, 2002

<div id="cite-51"></div>
'''[[#citeF-51|[51]]]''' Timoshenko, S and Goodier, JN. ''Theory of elasticity, 1951'', Volume 412. New York

<div id="cite-52"></div>
'''[[#citeF-52|[52]]]''' Crisfield, Michael A. ''Finite Elements and Solution Procedures for Structural Analysis V. 1: Linear Analysis'', Volume 1. Pineridge Press Limited, 1986

<div id="cite-53"></div>
'''[[#citeF-53|[53]]]''' Oñate, Eugenio. ''Structural Analysis with the Finite Element Method. Linear Statics: Volume 1: Basis and Solids'', Volume 1. Springer Science & Business Media, 2010

<div id="cite-54"></div>
'''[[#citeF-54|[54]]]''' Oñate, Eugenio. ''Structural Analysis with the Finite Element Method. Linear Statics: Volume 2: Beams, Plates and Shells'', Volume 2. Springer Science & Business Media, 2013

<div id="cite-55"></div>
'''[[#citeF-55|[55]]]''' Ramm, Ekkehard. ''Strategies for tracing the nonlinear response near limit points''. Springer, 1981

<div id="cite-56"></div>
'''[[#citeF-56|[56]]]''' Batoz, Jean-Louis and Dhatt, Gouri. ''Incremental displacement algorithms for nonlinear problems'', Volume 14. Wiley Online Library. International Journal for Numerical Methods in Engineering 8 1262&#8211;1267, 1979

<div id="cite-57"></div>
'''[[#citeF-57|[57]]]''' Crisfield, Michael A. ''Non-linear finite element analysis of solids and structures: Advanced topics''. John Wiley & Sons, Inc., 1997

<div id="cite-58"></div>
'''[[#citeF-58|[58]]]''' Portabella Castany, Miquel. ''Un procedimiento para cálculo de estructuras por el método de elementos finitos con error prefijado utilizando refinamiento de malla adaptativa''. Universitat Politecnica de Catalunya, 2014

<div id="cite-59"></div>
'''[[#citeF-59|[59]]]''' Oñate, Eugenio. ''Cálculo de estructuras por el método de elementos finitos: análisis elástico lineal''. Centro Internacional de Métodos Numéricos en Ingeniería, 1992

<div id="cite-60"></div>
'''[[#citeF-60|[60]]]''' Bugeda, G and Oñate, E. ''Adaptive mesh refinement techniques for aerodynamic problems.'' 513&#8211;522, 1992

<div id="cite-61"></div>
'''[[#citeF-61|[61]]]''' Oñate, Eugenio and Bugeda, Gabriel. ''A study of mesh optimality criteria in adaptive finite element analysis'', Volume 10. MCB UP Ltd. Engineering computations 4 307&#8211;321, 1993

<div id="cite-62"></div>
'''[[#citeF-62|[62]]]''' Le Bellego, Caroline and normale supérieure de Cachan, École. ''Couplages chimie-mécanique dans les structures en béton attaquées par l'eau: Etude expérimentale et analyse numérique''. LMT-ENS Cachan, 2001

<div id="cite-63"></div>
'''[[#citeF-63|[63]]]''' Rodríguez Ferran, Antonio and Arbós, Ivan and Huerta, Antonio and others. ''Adaptive analysis based on error estimation for nonlocal damage models'', 2001

<div id="cite-64"></div>
'''[[#citeF-64|[64]]]''' Carpinteri, Alberto and Valente, Silvio and Ferrara, G and Melchiorrl, G. ''Is mode II fracture energy a real material property?'', Volume 48. Elsevier. Computers & structures 3 397&#8211;413, 1993