Combination of a non-local damage model for quasi-brittle materials with a mesh-adaptive finite element technique


==Abstract==

Progressive fracture in quasi-brittle materials is often treated via strain softening models in continuum damage mechanics. Such constitutive relations favour spurious strain localization and ill-posedness of boundary value problems. The introduction of non-local damage models together with a characteristic length parameter controlling the size of the fracture process zone is known to regularize the problem. In order to account for the non-locality of these models, it is crucial to work with fine spatial discretizations at the damage progress zone. In this paper we present a non-local damage model in combination with a mesh-adaptive finite element technique that can help automatize the analysis of progressive fracture problems in an efficient manner. Classical two-dimensional examples are given to illustrate the presented approach.

'''Keywords''': Quasi-brittle materials, continuum damage mechanics, non-local damage models, finite element method, mesh-adaptivity

==1 Introduction==

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 <span id='citeF-14'></span>[[#cite-14|[14]]], but the concept of damage was introduced by Kachanov in 1958 in the context of creep rupture <span id='citeF-16'></span>[[#cite-16|[16]]]. In that work Kachanov introduced the concept of effective stress, and by using continuum damage he solved problems related to creep in metals. <span id='citeF-31'></span>[[#cite-31|[31]]] gave the problem physical meaning by suggesting that the reduction of the sectional area was measured by means of the damage parameter. The thermodynamic formalism involved in the irreversible process of damage was developed by <span id='citeF-18'></span>[[#cite-18|[18]]]. Other important contributions on damage mechanics include: <span id='citeF-22'></span> [[#cite-22|[22]]], <span id='citeF-33'></span> [[#cite-33|[33]]], <span id='citeF-26'></span> [[#cite-26|[26]]], <span id='citeF-24'></span> <span id='citeF-25'></span> [[#cite-24|[24, 25]]], <span id='citeF-8'></span> <span id='citeF-9'></span> [[#cite-8|[8, 9]]], to name but a few.

The behaviour of brittle or quasi-brittle materials such as concrete, rocks, mortar or other geo-materials is particularly difficult to predict. In those cases failure 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, 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 <span id='citeF-13'></span>[[#cite-13|[13]]].

If the damage parameter depends only on the strain state at the point under consideration, numerical simulations exhibit a pathological mesh dependence and the energy consumed by the fracture process tends to zero as the mesh is refined. 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 the distance between damaged zones. They suffer from mesh sensitivity (for size and alignment) and produce unreliable results. Strains concentrate in one element wide zones and the computed force-displacement curves are mesh-dependent. The reason behind these misbehaviours is that the differential equations of motion change their type (from elliptic to hyperbolic in static problems) and the boundary value problem becomes ill-posed <span id='citeF-2'></span>[[#cite-2|[2]]].

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 by micro-polar <span id='citeF-23'></span>[[#cite-23|[23]]], strain gradient <span id='citeF-34'></span>[[#cite-34|[34]]], viscous <span id='citeF-19'></span>[[#cite-19|[19]]] and non-local terms <span id='citeF-15'></span>[[#cite-15|[15]]]. In our model we have worked with the latter approach using a weighted spatial averaging of the internal variables. In this manner the kinematic and equilibrium equations remain standard, and the notions of stress and strain keep their usual meaning.

The first non-local 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 <span id='citeF-12'></span>[[#cite-12|[12]]] who later extended it to non-local elasto-plasticity <span id='citeF-11'></span>[[#cite-11|[11]]]. 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-30'></span>[[#cite-30|[30]]].

Non-local models lead to smooth solutions with a continuous variation of strain. However, to resolve narrow bands of highly localized strains using the finite element method it is necessary to use sufficiently fine computational grids. Fortunately, the mesh must be fine only in the damage progression 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 suitably construct refined meshes by hand. Thereby, the efficiency of the analysis can be greatly increased by means of an adaptive mesh refinement technique, which automates the whole process <span id='citeF-5'></span><span id='citeF-29'></span>[[#cite-5|[5, 29]]].

In the present work we present a robust non-local isotropic damage model for quasi-brittle materials that works in a small deformation regime, along with an adaptive-mesh finite element technique that permits adapting the spatial discretization in an optimal manner.

The paper is organized as follows. First, the basic concepts on continuum damage mechanics are introduced. Details are given on the basic components of the isotropic damage theory, and on the equivalent strain forms and damage evolution laws that have been implemented in this work.

Next, we review the regularization technique that has been used to overcome the problems associated to strain localization. The fundamental aspects of the integral-type non-local damage model derived are presented, pointing out the most relevant aspects of its numerical implementation. The method for estimating the error of the numerical solution and the mesh-adaptive scheme are explained in some detail.

Finally, two examples are presented showing that the combination of the non-local damage model and the mesh-adaptive technique is a robust method to model the failure of quasi-brittle materials.

==2 A simple isotropic damage model==

The simplest damage model for multiaxial stress states is the isotropic damage model with a simple scalar variable. This model is based on the 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. The stress-strain law is postulated as:

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

where <math display="inline">\boldsymbol{\sigma }</math> is the total stress tensor, <math display="inline">\boldsymbol{\varepsilon }</math> is the total strain tensor, <math display="inline">\bar{\boldsymbol{\sigma }}</math> is the effective stress tensor, <math display="inline">\mathbb{E}</math> is the elastic constitutive tensor of the intact material, and <math display="inline">d</math> is the scalar damage variable.

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 variable at a generic section of a material as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (2)
|}

where <math display="inline">S</math> and <math display="inline">\bar{S}</math> are respectively the total and the effective area of the section, and <math display="inline">S_d=S-\bar{S}</math> is the damaged part of the area. An undamaged material is characterized by <math display="inline">d=0</math>. Due to propagation and coalescence of micro-defects, 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 order to properly determine the evolution of the damage variable regardless of the loading case we must introduce a loading function <math display="inline">f</math> specifying the elastic domain and the states at which damage grows. The loading function 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. A typical definition for function <math display="inline">f</math> is

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>f(\boldsymbol{\varepsilon },r)=\varepsilon _{eq}(\boldsymbol{\varepsilon })-r </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

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> represents a scalar measure of the largest strain level ever reached in the previous deformation history of the material up to its current state, i.e.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (4)
|}

The above expression implies that <math display="inline">r(t)\ge r_0</math>, where <math display="inline">r_0</math> is the damage threshold, a material parameter that represents the value of equivalent 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).

We also postulate the loading-unloading conditions in the Kuhn-Tucker form:

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

The first condition indicates that <math display="inline">r</math> can never be smaller than <math display="inline">\varepsilon _{eq}</math>, while the second one 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 _{eq}</math> and <math display="inline">r</math> are equal.

The damage evolution law is defined as:

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>d=g(r)\hbox{  with }\left\{\begin{array}{cc}g(r)=0 & \hbox{ if }r=r_0\\ 0<g(r)\le{1} & \hbox{ if }r>r_0 \end{array}\right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

which holds not only during monotonic loading but also during unloading and reloading.

There are various damage governing laws <math display="inline">g(r)</math> that can be effectively used to model damage growth in quasi-brittle materials. In this work we adopt the exponential law proposed in [mazars1986description], which separates the damage in compression and tension as:

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

with:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (8)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (9)
|}

where parameters <math display="inline">A_t</math> and <math display="inline">A_c</math> are associated to residual strength, and parameters <math display="inline">B_t</math> and <math display="inline">B_c</math> control the slope of the softening branch at the peak of the stress-strain curve. The coefficient <math display="inline">\alpha _t</math> in ([[#eq-7|7]]) 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;width: 100%;" 
|-
| 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;white-space: nowrap;" | (10)
|}

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=\mathbb{E}^{-1}:\langle \mathbb{E}:\boldsymbol{\varepsilon }\rangle </math>.

In the second example of this work, however, we use the simplified version of the previous exponential law as:

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (11)
|}

Concerning the definition of the loading function <math display="inline">f</math> in ([[#eq-3|3]]), we must notice that the expression of the equivalent strain plays a role similar to the yield function in plasticity, because it directly affects the shape of the elastic domain.

Numerous forms of equivalent strain can be found in the literature. In this work, we have used two of them: the so-called <span id='citeF-20'></span>[[#cite-20|[20]]]:

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

and the modified version of <span id='citeF-10'></span>[[#cite-10|[10]]]:

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (13)
|}

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

<div id='img-1a'></div>
<div id='img-1b'></div>
<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG1a.png|250px|Mazars.]]
|[[Image:Draft_Samper_843551501-FIG1b.png|250px|Modified von Mises.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mazars.
| (b) Modified von Mises.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 1:''' Qualitative damage surfaces in the 2D principal stress space.
|}

As one can see in Figure [[#img-1|1]], the introduced forms of equivalent strain lead to non-symmetric damage surfaces in which the yield value in compression can be several times the value in tension. This is essential to account for the different strength in tension and compression of geo-materials such as concrete and rocks.

We summarize below the basic ingredients of the isotropic damage model:

* The stress-strain relation ([[#eq-1|1]])
* The damage criterion, comprised by:
**The loading function in ([[#eq-3|3]])
**The loading-unloading conditions in ([[#eq-5|5]])
*The law governing the evolution of the damage variable ([[#eq-6|6]])
*The expression defining the equivalent strain

==3 Non-local damage model==

Although isotropic damage 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.

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 <span id='citeF-4'></span>[[#cite-4|[4]]].

Integral-type non-local models abandon the classical assumption of locality and admit that stress at 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.

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{p})</math> be some local field in a domain <math display="inline">V</math>, the corresponding non-local field is defined as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{f}(\boldsymbol{p})=\int _V \alpha (\boldsymbol{p},\boldsymbol{q})f(\boldsymbol{q})d\boldsymbol{q} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">\alpha (\boldsymbol{p},\boldsymbol{q})</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{p}-\boldsymbol{q}\| </math> between the source point <math display="inline">\boldsymbol{q}</math>, and the receiver point <math display="inline">\boldsymbol{p}</math>. Thereby, we usually write <math display="inline">\alpha (\boldsymbol{p},\boldsymbol{q})=\alpha _0(\| \boldsymbol{p}-\boldsymbol{q}\| )</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 choice for <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;width: 100%;" 
|-
| 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;white-space: nowrap;" | (15)
|}

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-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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 \\ 0 & \hbox{ if }|D|>R \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

where the interaction radius <math display="inline">R</math> is a parameter usually 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>.

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

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG2.png|250px|Averaging zone near the boundary of a solid.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Averaging zone near the boundary of a solid.
|}

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 accomplish the normalizing condition:

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

In order to satisfy ([[#eq-17|17]]), the weighting function is rescaled as:

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha (\boldsymbol{p},\boldsymbol{q})=\frac{\alpha _0(\| \boldsymbol{p}-\boldsymbol{q}\| )}{\displaystyle \int _V\alpha _0(\| \boldsymbol{p}-\boldsymbol{x}\| )d\boldsymbol{x}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

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-3'></span>[[#cite-3|[3]]].

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

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

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-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (20)
|}

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-1|1]]) 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 fully local.

Numerical implementation of the non-local damage model based on averaging of the equivalent strain is relatively simple. The evaluation of the 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 the non-local equivalent strain must be traced at the individual Gauss integration points of the finite element mesh, because these are the points at which stresses are computed.

Thereby, the averaging integral in ([[#eq-19|19]]) is evaluated numerically. By considering the weighting function defined in ([[#eq-18|18]]) we can write:

<span id="eq-21"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\tilde{\varepsilon }_{eq,p}=\sum _qw_q\alpha _{pq}\varepsilon _{eq,q} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

where <math display="inline">w_q</math> is a coefficient containing the product of the determinant of the Jacobian and the integration weights of Gauss point <math display="inline">q</math>, 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-22"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\alpha _{pq}=\frac{\alpha _0(\| \boldsymbol{p}-\boldsymbol{q}\| )}{\sum _xw_x\alpha _0(\| \boldsymbol{p}-\boldsymbol{x}\| )} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

In the previous two equations, subscript <math display="inline">p</math> represents the receiver point under consideration, whereas indexes <math display="inline">q</math> and <math display="inline">x</math> correspond to source points. Since the chosen weighting function <math display="inline">\alpha _0</math> has bounded support ([[#eq-16|16]]), <math display="inline">\alpha _{pq}</math> vanishes if the distance 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-21|21]]) and ([[#eq-22|22]]) 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, note that at the damage process zone one must always use an element size smaller than the interaction radius, in order to account for the non-local interaction. Otherwise the damage model would become local.

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 domain of influence of the considered point (Figure [[#img-3|3]]).

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG3.png|300px|Grid-based non-local search.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' 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-21|21]]).

==4 Computation of the tangent 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 expressing the vector of internal forces as a numerical integration over the Gauss points. In order to make the explanation more understandable, from now on we will be using Voigt notation, i.e.

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

Using the stress-strain law ([[#eq-1|1]]) and the classic strain-displacement relation <math display="inline">\boldsymbol{\varepsilon }=\boldsymbol{B}\boldsymbol{a}</math>, we expand ([[#eq-23|23]]) as follows:

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

The tangent stiffness matrix is obtained by differentiating the internal forces with respect to the vector of nodal displacements <math display="inline">\boldsymbol{a}</math>. Since the damage variable depends on the nodal displacements through the equivalent strain, we will 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-6|6]]), <math display="inline">r</math> depends on <math display="inline">\tilde{\varepsilon }_{eq}</math> ([[#eq-20|20]]), and <math display="inline">\tilde{\varepsilon }_{eq}</math> depends on <math display="inline">\boldsymbol{a}</math> through the interpolated strains, we 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-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (25)
|}

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. Thus,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>t_p=\left\{\begin{array}{ll}0 & \hbox{ if }\tilde{\varepsilon }_{eq,p}<r_p \\ 1 & \hbox{ if }\tilde{\varepsilon }_{eq,p}=r_p\hbox{ and }\dot{\tilde{\varepsilon }}_{eq,p}\ge 0 \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

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

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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}}=\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;white-space: nowrap;" | (27)
|}

The derivative of the equivalent strain with respect to the strain vector <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-24|24]]) as follows:

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

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

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{tan}  =\frac{\partial \boldsymbol{f}_{int}}{\partial \boldsymbol{a}} </math>
|-
| style="text-align: center;" | <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;white-space: nowrap;" | (29)
|}

Note that the first term in the right hand side (r.h.s.) of ([[#eq-29|29]]) 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>). The second term in the r.h.s. is the non-local part of the tangent stiffness matrix. Substituting ([[#eq-27|27]]) into ([[#eq-29|29]]) yields:

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{tan}  =\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: center;" | <math>   =\boldsymbol{K}_{sec}-\sum _{p,q}w_p\boldsymbol{B}_p^T\boldsymbol{\bar{\sigma }}{_pg'}pw_qpw_q\alpha _{pq}\frac{\partial \varepsilon _{eq,q}}{\partial \boldsymbol{\varepsilon }_q}\boldsymbol{B}_q </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Defining for convenience the column matrix <math display="inline">\boldsymbol{\beta }^c_p=\boldsymbol{B}^T_p\boldsymbol{\bar{\sigma }}_p</math>, the row matrix <math display="inline">\boldsymbol{\beta }^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. ([[#eq-30|30]]) can be rewritten like:

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

The double index of the sum, caused by the non-local interaction, implies that the term on the right part of Eq. ([[#eq-31|31]]) cannot be assembled from element contributions only. 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> (Figure [[#img-4|4]]). Consequently, the global stiffness matrix is always non-symmetric and its bandwidth increases due to the non-local interaction.

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG4.png|200px|Non-local assembly process.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Non-local assembly process.
|}

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{K}_{tan}^{local}=\boldsymbol{K}_{sec}-\sum _{p}{g'}_pt_pw_{p}\boldsymbol{\beta }^c_p\boldsymbol{\beta }^r_p </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

where the sum is performed over one index only. Note that the resulting local tangent matrix <math display="inline">\boldsymbol{K}_{tan}^{local}</math> is no longer consistent, and quadratic convergence is lost.

Probably 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. However, 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.

The non-local averaging certainly increases the computational cost with respect to the corresponding local model. However, the non-local model completely removes the pathological sensitivity to the mesh size and partially alleviates the mesh-induced directional bias. Consequently, this extra computational effort is indeed worthwhile.

==5 Mesh-adaptive technique==

In this work we have used the error estimation and adaptive mesh refinement strategy presented by Oñate and co-workers [bugeda1992adaptive,bugeda1994methodology,onate1991adaptive,onate1993study,eugenio2009structural].

The general algorithm of non-linear adaptive FEM analysis implemented in this work is 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 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.

We present below the main stages of the implemented adaptive procedure:

===5.1 Error estimation===

The error estimator is based on the stress evaluation. We define the error as the difference between an exact value of the 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;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\epsilon }_{\bar{\sigma }}=\bar{\boldsymbol{\sigma }}_e-\bar{\boldsymbol{\sigma }}_a </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

Note that we compute the error with the effective stresses <math display="inline">\bar{\boldsymbol{\sigma }}</math> and not with the total stresses <math display="inline">\boldsymbol{\sigma }</math> because the latter tend to zero as damage grows.

A reasonably good value of the exact stresses is obtained from their extrapolation from the Gauss points to the nodes and a posterior smoothing [eugenio2009structural]. 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;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\epsilon }_{\bar{\sigma }}\approx \bar{\boldsymbol{\sigma }}_s-\bar{\boldsymbol{\sigma }} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

We will work with integral measures of the error. The energy norm of the error over an element is defined as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| ^{(e)} =\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: center;" | <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;white-space: nowrap;" | (35)
|}

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;width: 100%;" 
|-
| 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;white-space: nowrap;" | (36)
|}

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

===5.2 Remeshing===

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

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

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

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

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

with

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (39)
|}

To determine whether condition ([[#eq-37|37]]) 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;width: 100%;" 
|-
| style="text-align: center;" | <math>\xi _g=\frac{\| \boldsymbol{\epsilon }_{\bar{\sigma }}\| }{\eta \,U} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

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

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;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat{h}^{(e)}=\frac{h^{(e)}}{(\xi _g)^{\frac{1}{m}}} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

If one aims at obtaining a selective refinement method, apart from ([[#eq-37|37]]), another condition concerning the error of each element must be simultaneously imposed. In this paper, we have chosen a remeshing criterion based on a global error equidistribution [onate1993study].

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

<span id="eq-42"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (42)
|}

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

Like for 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;width: 100%;" 
|-
| 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;white-space: nowrap;" | (43)
|}

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 us to define the new element size that accomplishes ([[#eq-42|42]]):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (44)
|}

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

In the end, the remeshing criterion results from the combination of the global error condition ([[#eq-37|37]]) and the local one ([[#eq-42|42]]). 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;width: 100%;" 
|-
| 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;white-space: nowrap;" | (45)
|}

And the new element size is obtained with:

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

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| 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;white-space: nowrap;" | (47)
|}

===5.3 Mesh generator interface===

The mesh generator interface is the pre-processing and post-processing GiD software used to run the simulations [GiD_web]. Therefore, not only GiD meshes the geometry at the beginning of the numerical solution, 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 is generated with the information of the new element sizes. Then, GiD allows us to load 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.

===5.4 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 ones (Figure [[#img-5a|5a]]).

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 (Figure [[#img-5|5]]). Once again, 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-5a'></div>
<div id='img-5b'></div>
<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG5a.png|504px|Nodal displacements mapping.]]
|[[Image:Draft_Samper_843551501-FIG5b.png|548px|Internal state variables mapping.]]
|- style="text-align: center; font-size: 75%;"
| (a) Nodal displacements mapping.
| (b) Internal state variables mapping.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 5:''' Mapping of variables.
|}

==6 Examples==

In this section we will present two examples of application of the non-local damage model and the mesh adaptive technique proposed. The examples are the three-point bending test and the four-point shear test. For each case we will first assess the robustness of the non-local damage model when working with different spatial discretizations, and then we will solve the problem again applying the mesh-adaptive technique, in order to point out the strengths and limitations of the code in its current state.

Regarding the solving strategy, we will be following the equilibrium path of the problem for a series of steps through a global arc-length method in an incremental-iterative process. Self weight will not be considered.

===6.1 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. The notch is 3 mm thick and extends over one tenth of the beam depth (Figure [[#img-6|6]]).

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

Plane stress conditions have been assumed. The geometry was meshed by 4-noded quadrilaterals with 2 <math display="inline">\times </math> 2 integration points.

====Non-local damage model solution====

As it has been stated before, the non-local damage model was defined with the Mazars model, regarding the equivalent strain ([[#eq-12|12]]) and the damage evolution law ([[#eq-7|7]]). The material parameters were obtained from [le2001couplages] and are summarized in Table [[#table-1|1]].


{|  class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 TPB test. Material parameters for Mazars non-local damage model.
|-
|  '''Parameter''' 
| '''Value''' 
|-
|  Young 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 non-local damage model, we have also solved the problem with a local damage model based on the damage evolution law of [oliver1990isotropic] and a modified definition of the equivalent strain presented by Simo and Ju in [simo1987strain]. The parameters for this model are shown in Table [[#table-2|2]].


{|  class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 TPB test. Material parameters for Simo-Ju local damage model.
|-
|  '''Parameter''' 
| '''Value''' 
|-
|  Young 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 

|}

The problem was solved for different spatial discretizations. In this case we used three unstructured meshes of 4-noded quadrilaterals with a minimum size of 15 mm, 7 mm, and 3 mm, respectively (Figure [[#img-7|7]]).

<div id='img-7a'></div>
<div id='img-7b'></div>
<div id='img-7c'></div>
<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG7a.png|600px|Mesh 1: 2024 4-noded elements.]]
|[[Image:Draft_Samper_843551501-FIG7b.png|600px|Mesh 2: 2679 4-noded elements.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1: 2024 4-noded elements.
| (b) Mesh 2: 2679 4-noded elements.
|-
| colspan="2"|[[Image:Draft_Samper_843551501-FIG7c.png|600px|Mesh 3: 6543 4-noded elements.]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (c) Mesh 3: 6543 4-noded elements.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 7:''' TPB test. Spatial discretizations.
|}

Figure [[#img-8|8]] shows the relation between the applied load and the vertical deflection of the beam for each mesh. As one can see from the discontinued equilibrium curves, we had serious difficulties in tracing the response of a full test. The reason behind the convergence problems could be the use of a too global arc-length method. 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-8a|8a]] we can see that, although there is no relevant difference between the response obtained with meshes 1 and 2, the peak load actually decreases with the finest mesh 3, and so the total dissipated energy. On the other hand, the response diagram in Figure [[#img-8b|8b]] shows practically the same peak load for the three meshes.

<div id='img-8a'></div>
<div id='img-8b'></div>
<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG8a.png|600px|Local damage model.]]
|[[Image:Draft_Samper_843551501-FIG8b.png|600px|Non-local damage model.]]
|- style="text-align: center; font-size: 75%;"
| (a) Local damage model.
| (b) Non-local damage model.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 8:''' TPB test. Force-vertical deflection diagrams.
|}

Figure [[#img-9|9]] shows clearly the different behaviour of the local and non-local approaches. While in the local model the thickness of the damage pattern is of the order of the element size, in the non-local case the damage pattern is controlled by the interaction radius <math display="inline">R</math> and so the response remains virtually unaltered regardless of the mesh.

<div id='img-9a'></div>
<div id='img-9b'></div>
<div id='img-9c'></div>
<div id='img-9d'></div>
<div id='img-9e'></div>
<div id='img-9f'></div>
<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG9a.png|600px|Local model. Initial stage.]]
|[[Image:Draft_Samper_843551501-FIG9b.png|600px|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_Samper_843551501-FIG9c.png|600px|Local model. Middle stage.]]
|[[Image:Draft_Samper_843551501-FIG9d.png|600px|Non-local model. Middle stage.]]
|- style="text-align: center; font-size: 75%;"
| (c) Local model. Middle stage.
| (d) Non-local model. Middle stage.
|-
|[[Image:Draft_Samper_843551501-FIG9e.png|600px|Local model. Advanced stage.]]
|[[Image:Draft_Samper_843551501-FIG9f.png|600px|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 9:''' TPB test. Evolution of damage propagation.
|}

Furthermore, in order to check whether the implemented model can properly reproduce the behaviour of the real test, we have compared the results obtained in the non-local model with experiments performed in [le2001couplages].

Figure [[#img-10|10]] shows the force-vertical deflection diagram of the experimental solution and the computed solution with mesh 3.

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG10.png|420px|TPB test. Non-local model validation.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' 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 experimental solution, and even shows a certain amount of snap-back behaviour. 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 that observed in the experiment. This difference could be due to the approach used in the laboratory test to control the post-peak behaviour.

====Mesh-adaptive solution====

The same problem was solved again, using the mesh-adaptive technique described along with the non-local damage model.

As mentioned in Section [[#5.2 Remeshing|5.2]], we use a remeshing criterion based on the global error equidistribution. The objective of the example was to obtain efficient spatial discretizations for the changing damage states, limiting the global error to a 10% of the deformation energy <math display="inline">(\eta=10%)</math>. The initial mesh is an unstructured mesh of 4-noded quadrilaterals with an average size of 10 mm.

In this case three refinements have been performed. In order to understand the mesh-adaptive process, it is interesting to see the evolution of damage for the different meshes (Figure [[#img-11|11]]).

<div id='img-11a'></div>
<div id='img-11b'></div>
<div id='img-11c'></div>
<div id='img-11d'></div>
<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG11a.png|600px|Mesh 0: 4453 4-noded elements. ξ<sub>g</sub>=1.33]]
|[[Image:Draft_Samper_843551501-FIG11b.png|600px|Mesh 1: 2977 4-noded elements. ξ<sub>g</sub>=0.67]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 0: 4453 4-noded elements. <math>\xi _g=1.33</math>
| (b) Mesh 1: 2977 4-noded elements. <math>\xi _g=0.67</math>
|-
|[[Image:Draft_Samper_843551501-FIG11c.png|600px|Mesh 2: 1598 4-noded elements. ξ<sub>g</sub>=0.92]]
|[[Image:Draft_Samper_843551501-FIG11d.png|600px|Mesh 3: 1627 4-noded elements. ξ<sub>g</sub>=0.93]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 2: 1598 4-noded elements. <math>\xi _g=0.92</math>
| (d) Mesh 3: 1627 4-noded elements. <math>\xi _g=0.93</math>
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 11:''' TPB test. Progression of damage in the adaptive mesh.
|}

A quantitative measure of the suitability of the mesh is given by the global refinement parameter <math display="inline">\xi _g</math>. A value close to 1 means that the global error is similar to the goal percentage of error <math display="inline">\eta </math>. The value of 1.33 in Figure [[#img-11a|11a]] denotes that the mesh must be refined in those areas of the structure with higher stress variations, i.e., the zone near the notch, the point of application of the load, and the supports. On the other hand, since this procedure also optimizes the mesh enlarging the elements on those zones with smoother stresses, the number of elements can be reduced from 4453 to 2977.

In the subsequent meshes (Figures [[#img-11b|11b]], [[#img-11c|11c]] and [[#img-11d|11d]]) the global refinement parameters are smaller than 1, showing that the mesh-adaptive technique properly keeps the error below the prescribed threshold.

We can say that the main advantage of the 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. However, this procedure still does not distinguish compressive stresses from tensile stresses when estimating the error. In materials like concrete, this can lead to excessively small elements in compressed zones, and so it could be convenient to give more weight to the tensile stresses.

===6.2 Four-point shear test===

The test will be performed with a single-edge notched beam subjected to four-point shear (FPS). 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 (Figure [[#img-12|12]]).

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG12.png|450px|FPS test. Problem statement. Distances in mm]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' FPS test. Problem statement. Distances in mm
|}

Plane stress conditions have been assumed. The geometry was meshed using standard linear 3-noded triangles with one integration point.

====Non-local damage model solution====

The non-local damage model is defined with the equivalent strain form of the modified von Mises model ([[#eq-13|13]]), and the simplified exponential damage evolution law ([[#eq-11|11]]). The material parameters for the non-local approach  have been obtained from [rodriguez2001adaptive] and are summarized in Table [[#table-3|3]].


{|  class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 FPS test. Material parameters for modified von Mises non-local damage model.
|-
|  '''Parameter''' 
| '''Value''' 
|-
|  Young 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 

|}

Like in the first example, we have compared the non-local damage model with a local damage model (Table [[#table-4|4]]).


{|  class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-4'></span>Table. 4 FPS test. Material parameters for Simo-Ju local damage model.
|-
|  '''Parameter''' 
| '''Value''' 
|-
|  Young 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 

|}

In this case, we used three unstructured meshes of 3-noded triangles with a minimum size of 8 mm, 5 mm, and 3 mm, respectively (Figure [[#img-13|13]]).

<div id='img-13a'></div>
<div id='img-13b'></div>
<div id='img-13c'></div>
<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG13a.png|600px|Mesh 1: 2216 3-noded elements.]]
|[[Image:Draft_Samper_843551501-FIG13b.png|600px|Mesh 2: 3502 3-noded elements.]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 1: 2216 3-noded elements.
| (b) Mesh 2: 3502 3-noded elements.
|-
| colspan="2"|[[Image:Draft_Samper_843551501-FIG13c.png|600px|Mesh 3: 7183 3-noded elements.]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (c) Mesh 3: 7183 3-noded elements.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 13:''' FPS test. Spatial discretizations.
|}

In Figure [[#img-14|14]] we represent the relation between the applied load and the Crack Mouth Sliding Displacement (CMSD) for each mesh. The curves in Figure [[#img-14a|14a]] show that in this case the peak and the dissipated energy also decrease as the mesh is refined. However, this reduction is not so clear as in the previous example. On the other hand, Figure [[#img-14b|14b]] shows virtually the same peak load for the three meshes. Also, although we can notice some oscillations in the post-peak regions of meshes 1 and 2, we can say that the residual force at the right part of the graph seems to be pretty similar for all cases.

<div id='img-14a'></div>
<div id='img-14b'></div>
<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG14a.png|600px|Local damage model.]]
|[[Image:Draft_Samper_843551501-FIG14b.png|600px|Non-local damage model.]]
|- style="text-align: center; font-size: 75%;"
| (a) Local damage model.
| (b) Non-local damage model.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 14:''' FPS test. Force-Crack Mouth Sliding Displacement diagrams.
|}

Figure [[#img-15|15]] shows the damage progression for both models stressing the more localized nature of the local damage model. Regarding the non-local model, if we compare this example with the three-point bending test, we can clearly see that the damage 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-15a'></div>
<div id='img-15b'></div>
<div id='img-15c'></div>
<div id='img-15d'></div>
<div id='img-15e'></div>
<div id='img-15f'></div>
<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG15a.png|600px|Local model. Initial stage.]]
|[[Image:Draft_Samper_843551501-FIG15b.png|600px|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_Samper_843551501-FIG15c.png|600px|Local model. Middle stage.]]
|[[Image:Draft_Samper_843551501-FIG15d.png|600px|Non-local model. Middle stage.]]
|- style="text-align: center; font-size: 75%;"
| (c) Local model. Middle stage.
| (d) Non-local model. Middle stage.
|-
|[[Image:Draft_Samper_843551501-FIG15e.png|600px|Local model. Advanced stage.]]
|[[Image:Draft_Samper_843551501-FIG15f.png|600px|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 15:''' FPS test. Evolution of damage propagation.
|}

Figure [[#img-16|16]] shows the relation between the applied load and the Crack Mouth Sliding Displacement of a reference experimental solution [carpinteri1993mode] and the computed solution with the finest mesh.

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG16.png|420px|FPS test. Non-local model validation.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' FPS test. Non-local model validation.
|}

We can see that the obtained response is very similar to the experimental solution. What should be noticed, though, is that the failure obtained in the computed solution is slightly more brittle as compared to the experimental one. This is possibly consequence of an imprecise tracing of the equilibrium path, caused by the global arc-length strategy used.

====Mesh-adaptive solution====

This example is performed with the same permissible percentage of global error as in the previous one <math display="inline">(\eta=10%)</math>, but the initial mesh is an unstructured mesh of triangles with an average size of 7 mm.

In this case the number of refinements is four. The evolution of damage for the different spatial discretizations is represented in Figure [[#img-17|17]].

<div id='img-17a'></div>
<div id='img-17b'></div>
<div id='img-17c'></div>
<div id='img-17d'></div>
<div id='img-17e'></div>
<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Samper_843551501-FIG17a.png|600px|Mesh 0: 7784 3-noded elements. ξ<sub>g</sub>=1.84]]
|[[Image:Draft_Samper_843551501-FIG17b.png|600px|Mesh 1: 6800 3-noded elements. ξ<sub>g</sub>=1.05]]
|- style="text-align: center; font-size: 75%;"
| (a) Mesh 0: 7784 3-noded elements. <math>\xi _g=1.84</math>
| (b) Mesh 1: 6800 3-noded elements. <math>\xi _g=1.05</math>
|-
|[[Image:Draft_Samper_843551501-FIG17c.png|600px|Mesh 2: 6590 3-noded elements. ξ<sub>g</sub>=2.20]]
|[[Image:Draft_Samper_843551501-FIG17d.png|600px|Mesh 3: 10626 3-noded elements. ξ<sub>g</sub>=0.44]]
|- style="text-align: center; font-size: 75%;"
| (c) Mesh 2: 6590 3-noded elements. <math>\xi _g=2.20</math>
| (d) Mesh 3: 10626 3-noded elements. <math>\xi _g=0.44</math>
|-
| colspan="2"|[[Image:Draft_Samper_843551501-FIG17e.png|600px|Mesh 4: 2542 3-noded elements. ξ<sub>g</sub>=0.94]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (e) Mesh 4: 2542 3-noded elements. <math>\xi _g=0.94</math>
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 17:''' FPS test. Progression of damage in the adaptive mesh.
|}

Figure [[#img-17a|17a]] shows no damage before the first refinement. It can be deduced then that the subsequent spatial discretization is mainly caused by the high compressive stresses at the loading points and at the supports.

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

Another aspect to consider here is directly related to the abrupt failure of the computed solution seen in Figure [[#img-16|16]]. Indeed, when damage grows abruptly the global refinement parameter <math display="inline">\xi _g</math> increases up to a value of 2.20 (Figure [[#img-17c|17c]]). Consequently, the adaptive technique refines excessively all the affected zone, leading to a considerably large number of elements and a much smaller global refinement parameter (Figure [[#img-17d|17d]]). Nonetheless, once damage progress stabilises, the adaptive procedure readjusts the spatial discretization and the resulting mesh is actually efficient (Figure [[#img-17e|17e]]).

Finally, we can notice that the damage pattern in Figure [[#img-17e|17e]] is a bit wider than the one shown in Figure [[#img-15f|15f]]. 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, in the near future it would be important to assess the effect of these errors in the traced equilibrium path, and improve the transfer algorithms if necessary.

==7 Conclusions==

The implemented non-local damage model has shown a robust behaviour for changing spatial discretizations, avoiding pathological mesh sensitivity and mesh-induced directional bias. It has also proved to properly capture the experimental peak load of the response diagrams in the two examples analysed: the three-point bending test and the four-point shear test. However, the numerical solution has shown a more brittle post-peak branch as compared to the experimental results. 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.

That aside, it must be stated that the global arc-length strategy used in the simulations 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. Indeed, the implementation of an advanced arc-length method with a specific control parameter (CMOD or CMSD) that accounts for the localized nature of quasi-brittle materials could help improving the accuracy of the converged solution.

The implemented mesh-adaptive technique allows one to capture 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. Hence, efficient mesh distributions can be obtained. Nevertheless, the mesh adaptive technique is still a recent implementation and there are some features that will be improved in subsequent work. For instance, the procedure does not distinguish compressive stresses from tensile stresses when estimating the error, which can result in 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.

==Acknowledgements==

The authors are thankful to Mr. Miquel Portabella for his collaboration in the implementation of the mesh-adaptive technique.

This research was partially funded by the Advanced Grant project SAFECON of the European Research Council and project NUMEXAS of the FP7 of the European Commission.

==REFERENCES==

<div id="cite-1"></div>
[[#citeF-1|[1]]]  GiD the personal pre and post processor. http://www.gidhome.com/,
April 2015.

<div id="cite-2"></div>
[[#citeF-2|[2]]]  Z. P. Bazant, T. B. Belytschko, and T.-P. Chang. Continuum theory
for strain-softening. Journal of Engineering Mechanics, 110(12):1666-
1692, 1984.

<div id="cite-3"></div>
[[#citeF-3|[3]]]  Z. P. Bazant and M. Jirásek. Nonlocal integral formulations of plasticity
and damage: survey of progress. Journal of Engineering Mechanics,
128(11):1119-1149, 2002.

<div id="cite-4"></div>
[[#citeF-4|[4]]]  Z. P. Bazant and G. Pijaudier-Cabot. Nonlocal continuum damage,
localization instability and convergence. Journal of applied mechanics,
55(2):287-293, 1988

<div id="cite-5"></div>
[[#citeF-5|[5]]] G. Bugeda and E. Oñate. Adaptive mesh refinement techniques for
aerodynamic problems. In Proceedings of the International Conference
on Métodos Numéricos en Ingeniería y Ciencias Aplicadas. Universidad
de Concepción, Chile. 16-20 November 1992, pages 513-522. Published
by CIMNE, Barcelona, 1992.

<div id="cite-6"></div>
[[#citeF-6|[6]]] G. Bugeda and E. Oñate. A methodology for adaptive mesh refinement
in optimum shape design problems. Computing Systems in Engineering,
5(1):91-102, 1994.

<div id="cite-7"></div>
[[#citeF-7|[7]]]  A. Carpinteri, S. Valente, G. Ferrara, and G. Melchiorrl. Is mode ii
fracture energy a real material property? Computers & structures,
48(3):397-413, 1993.

<div id="cite-8"></div>
[[#citeF-8|[8]]]  M. Cervera and M. Chiumenti. Mesh objective tensile cracking via a
local continuum damage model and a crack tracking technique. Com-
puter Methods in Applied Mechanics and Engineering, 196(1):304-320,
2006.

<div id="cite-9"></div>
[[#citeF-9|[9]]]  M. Cervera, M. Chiumenti, and C. A. de Saracibar. Shear band localization
via local j 2 continuum damage mechanics. Computer Methods
in Applied Mechanics and Engineering, 193(9):849-880, 2004.

<div id="cite-10"></div>
[[#citeF-10|[10]]]  J. De Vree, W. Brekelmans, and M. Van Gils. Comparison of nonlocal
approaches in continuum damage mechanics. Computers & Structures,
55(4):581-588, 1995.

<div id="cite-11"></div>
[[#citeF-11|[11]]]  A. C. Eringen. On nonlocal plasticity. International Journal of Engi-
neering Science, 19(12):1461-1474, 1981.

<div id="cite-12"></div>
[[#citeF-12|[12]]] A. C. Eringen and D. Edelen. On nonlocal elasticity. International
Journal of Engineering Science, 10(3):233-248, 1972

<div id="cite-13"></div>
[[#citeF-13|[13]]] M. Fafard and B. Massicotte. Geometrical interpretation of the arc-length
method. Computers & structures, 46(4):603-615, 1993

<div id="cite-14"></div>
[[#citeF-14|[14]]] J. Hult. Creep in continua and structures. In Topics in applied contin-
uum mechanics, pages 137-155. Springer, 1974.

<div id="cite-15"></div>
[[#citeF-15|[15]]] M. Jir_asek and S. Rolshoven. Comparison of integral-type nonlocal
plasticity models for strain-softening materials. International Journal
of Engineering Science, 41(13):1553-1602, 2003.

<div id="cite-16"></div>
[[#citeF-16|[16]]] L. Kachanov. On the time of fracture under conditions of creep. Izv.
AN SSSR, Otd. Tekh. Nauk, (8):26-35, 1958.

<div id="cite-17"></div>
[[#citeF-17|[17]]] C. Le Bellego and É. normale supérieure de Cachan. 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-18"></div>
[[#citeF-18|[18]]] J. Lemaitre and J. Chaboche. Mécanique des matériaux solides, 1985.
Dunod, Paris.

<div id="cite-19"></div>
[[#citeF-19|[19]]]  B. Loret and J. H. Prevost. Dynamic strain localization in elasto-
(visco-) plastic solids, part 1. general formulation and one-dimensional
examples. Computer Methods in Applied Mechanics and Engineering,
83(3):247-273, 1990

<div id="cite-20"></div>
[[#citeF-20|[20]]] J. Mazars. Application de la mécanique de l'endommagement au com-
portement non linéaire et à la rupture du béton de structure. PhD
thesis, 1984.

<div id="cite-21"></div>
[[#citeF-21|[21]]] J. Mazars. A description of micro-and macroscale damage of concrete
structures. Engineering Fracture Mechanics, 25(5):729-737, 1986.

<div id="cite-22"></div>
[[#citeF-22|[22]]]  J. Mazars and G. Pijaudier-Cabot. From damage to fracture mechanics
and conversely: a combined approach. International Journal of Solids
and Structures, 33(20):3327-3342, 1996.

<div id="cite-23"></div>
[[#citeF-23|[23]]] H. Mühlhaus. Shearband analysis in antigranulocytes material by
cosserat-theory. Ingenieur Archiv, 56(5):389-399, 1986.

<div id="cite-24"></div>
[[#citeF-24|[24]]]  J. Oliver, M. Cervera, S. Oller, and J. Lubliner. Isotropic damage
models and smeared crack analysis of concrete. In Second international
conference on computer aided analysis and design of concrete structures,
volume 2, pages 945-958, 1990.

<div id="cite-25"></div>
[[#citeF-25|[25]]]  J. Oliver, M. Cervera, S. Oller, and J. Lubliner. A simple damage
model for concrete, including long term effects. In Second International
Conference on Computer Aided Analysis And Design of Concrete Struc-
tures, volume 2, pages 945-958, 1990

<div id="cite-26"></div>
[[#citeF-26|[26]]]  S. Oller, E. Oñate, J. Oliver, and J. Lubliner. Finite element nonlinear
analysis of concrete structures using a ‘’plastic-damage model’’.
Engineering Fracture Mechanics, 35(1):219-231, 1990.

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

<div id="cite-28"></div>
[[#citeF-28|[28]]]  E. Oñate and G. Bugeda. A study of mesh optimality criteria in adaptive
finite element analysis. Engineering Computations, 10(4):307-321,
1993.

<div id="cite-29"></div>
[[#citeF-29|[29]]]  E. Oñate and J. Castro. Adaptive mesh refinement techniques for structural
problems. In The finite element method in the 1990's, pages 133-
145. Springer, 1991.
<div id="cite-30"></div>
[[#citeF-30|[30]]]  G. Pijaudier-Cabot and Z. P. Bazant. Nonlocal damage theory. Journal
of Engineering Mechanics, 113(10):1512-1533, 1987

<div id="cite-31"></div>
[[#citeF-31|[31]]]  Y. Rabotnov. Damage from creep. Zhurn. Prikl. Mekh. Tekhn. Phys,
2:113-123, 1963

<div id="cite-32"></div>
[[#citeF-32|[32]]]  A. Rodríguez Ferran, I. Arbós, A. Huerta, et al. Adaptive analysis
based on error estimation for nonlocal damage models. 2001.

<div id="cite-33"></div>
[[#citeF-33|[33]]]  J. Simo and J. Ju. Strain-and stress-based continuum damage models-i.
formulation. International journal of solids and structures, 23(7):821-
840, 1987.

<div id="cite-34"></div>
[[#citeF-34|[34]]]  H. Zbib and E. Aifantis. A gradient-dependent ow theory of plasticity:
application to metal and soil instabilities. Applied Mechanics Reviews,
42(11S):S295-S304, 1989.