A new a posteriori error indicator and its mesh optimization for nonlocal concrete damage models

Damage models are often used in finite element simulations of concrete structures. Evaluating the numerical quality of these computations consists in an important task which can be performed by using an error estimator. However, in the framework of nonlocal approaches these estimators seem to be much less abundant in the literature and their implementation requires considerable efforts. In this paper, a new error indicator with relatively less difficulties for its implementation is proposed. Based on local and least squares smoothing techniques, the convergence analyses of the developed indicator is carried out and discussed. The obtained results seem to be in good agreement with those obtained by the use of residual method. In addition the major difficulty to obtain an automatic step-by-step mesh adaptation is highlighted. Nevertheless for a chosen loading step, optimized calculations can be performed and example tests are presented.


INTRODUCTION
Damages in concrete structures are usually analyzed using the continuum theory of damage mechanics.This approach describes the influence of the progressive loss of material integrity due to the propagation and coalescence of microcracks and microvoids.Indeed, these changes in the microstructure lead to a degradation of material stiffness observed on the macro scale level.To characterize the density and orientation of microdefects a local representative internal variable of damage was introduced and a first damage model has been introduced by Kachanov (1958) and developed by Rabotnov (1969).Since many other isotropic and anisotropic models have been developed.For the isotropic models that we will limit ourselves, we can cite the model of Mazars (1984) and Mazars et al (2014).This model is very used for the monotonic loading cases.For its part, the modified Von Mises model also called MVM model by De Vree et al. (1995) seems to be very used in the literature.For our simulations, the model named RICRAG that is developed by Richard et al. (2010) and briefly exposed in the following section, will be used.This model seems to be not only relatively recent but also efficient.
The resolution of these nonlinear problems is often carried out by finite element approximation.However, it has been proved that these local models were no suitable for materials having a softening behavior (Bazant, 1986).Indeed, it was found that as soon as the mesh is refined, the localization tends to smaller areas.It follows that the dissipated energy during the rupture tends to vanish.This leads to a physical unacceptable result (Needleman, 1988).To eliminate this mesh dependence, different localization limiters have been developed.The method performed by Pijaudier-cabot and Bazant (1987) seemed to be for a long time among the most used.An improvement of this nonlocal regularization technique named stress based method is proposed by Giry et al. (2011).This will also be briefly exposed in the following section and used in our simulations.
Moreover, the simulated behavior is strongly influenced by discretization errors, which must be quantified and controlled.Otherwise, the effect of these errors on the results could be erroneously interpreted.To ensure the quality of the FE solution, an important task consists in evaluating the numerical quality of the solution by using a post-analyze or a posteriori error estimator.For linear problems, several different approaches leading to various estimators have been developed.In particular, we cite the estimator introduced by Babuska and Rheinboldt (1978) and based on the residual of the equilibrium equations; the estimators linked to the smoothing of finite element stresses performed by Zienkiewicz andZhu (1987 and1992); and the estimator based on the errors in the constitutive relation developed by Ladeveze (1975) and Ladeveze and Leguillon (1983).A review of different a posteriori estimators can be found in Zienkiewicz (2006) and Verfürth (2013).For nonlinear problems and especially for the non-linearity of damage mechanics, the works available are much less abundant.However, an estimator based on residuals and local computations has been developed by Rodríguez-Ferran and Huerta (2000).Furthermore, this method has been improved by the recent contributions of Gerasimov et al. (2015).On the other hand, Comi and Perego (2002) adapted the error in the constitutive relations cited above to a nonlocal damage model they developed (Comi and Perego, 2001).However, the implementation of these estimators seems requiring considerable efforts.
In this paper, a new error indicator using local and least squares smoothing techniques is presented.Indeed both the stress and damage field finite element solutions are smoothed and considered as reference solutions.Though ancient, the smoothed technique that we use seems to be not only simple but also efficient.It was developed by Hinton and Campbell (1974) and Hinton et al. (1975).In addition, the present error indicator is built by the explicit combination of two methods.The first is based on the error obtained by the energy standard difference between the smoothed stress solution and the finite element one.By the same manner, the second considers the error obtained by the difference between the smoothed damage solution and that obtained numerically.
On the other hand in order to obtain the mesh optimization with respect to a user prescribed error and thus to reduce the computations CPU costs, the knowledge of the properties of the present indicator is required.For this purpose, an analysis of the convergence properties, and also the efficiency of this indicator must be carried out and discussed.This optimization is generally obtained by a mesh adaptation operation performed by using a "size card" or a contribution card to the global error provided by the indicator.Thereafter by using both this indicator and procedures of mesh adaptivities, examples of optimized computations will be presented.It consists on the wellknown applications of simple and double edge notched beams respectively called SENB and DENB tests.

NONLOCAL DAMAGE FORMULATION
In this paper, the RICRAG model developed by Richard et al. (2010) is one we used during the implementation of the present error indicator.Known for its efficiency and robustness, this model is already implanted in the finite element code Cast3m (2018) developed by the CEA in France and used in our simulations.This section consists in a brief presentation of its basic principles and essential formulation.The exhaustive presentation of this model is beyond the scope of this paper and can be found in Richard et al. (2010).Formulated within the framework of the thermodynamics for irreversible processes a state potential was introduced such as it takes into account for elasticity, damage, sliding between cracked surfaces and hardenings as follows: By coupling damage with elasticity, the first and second terms (into accolades) of the equation (1) above introduce a single scalar damage variable d in order to represent the stiffness reduction.ñ is the material density,  and ì are the bulk and shear coefficients respectively. is the second order total strain tensor and  is the second order deviatoric total strain tensor such as  =  −   where δ is the second order Kronecker tensor.Let us recall that the scalar damage variable d evolves from 0 to 1 ie: from virgin to failed state material respectively.It should be also noted that the unilateral effect is taken into account by the separation of the total strain tensor into positive and negative parts as shown in the first term.The symbol < >+ denotes the positive part of the tensor.The third term 2  −   −  represent the sliding mechanism by contact with friction where  is the second order sliding tensor.
Finally, the last terms   and () are introduced to take into account the part of the strain energy retained by friction and hardening energy rates.Both kinematic and isotropic hardenings are considered;  is a material parameter,  the second order tensor associated to the kinematic hardening, z the internal variable related to the isotropic hardening and  the consolidation function.
After development (see Richard and al, 2010), the first state equations can be deduced such as: Latin American Journal of Solids and Structures, 2018, 15(9), e104 3/25 The second state equation in (2-b) allows obtaining the damage energy released rate Y.It can be obtained by differentiating the state potential with respect to the scalar damage variable: A threshold surface, denoted  is introduced: where  ̇ is a Lagrange multiplier and  a thermodynamic force associated to isotropic hardening: denotes the energy variable driving damage and  an initial threshold.This variable  is defined through the decomposition between direct extension mechanisms (tension) and induced extension mechanisms (compression): where  is the fourth order Hooke tensor computed from the elastic parameters  and ì. â is a parameter driving the dissymmetry of the threshold surface between tension and compression.The direct and induced extensions tensors are obtained through the following decompositions: H is the Heaviside function.Moreover, in order to avoid the mesh dependence of the strain localization phenomena, a nonlocal approach was used as a regularization technique by using the stress based method developed by Giry et al. (2011).In the following, we present this technique in an adapted version to the RICRAG model.Indeed, it consists in averaging the damage threshold surface in the vicinity of the current Gauss point.For this purpose, the local damage energy released rate  is replaced by the nonlocal one  such as the damage threshold surface becomes:  =  − ( + )   ̇=  ̇ . (10) ( − ) is the weight or also called Gauss distribution function defining the interaction between the considered point located with  and the neighboring points located with  inside the volume of the structure () such as: is a characteristic length which limits the domain of neighboring points and (,  ()) is a scalar which describes the interaction between points.It depends on the point location  and  () the principal stress state on the surrounding point s expressed in its principal stress reference system (,   (),   (),   ()) such as: () = ∑   () ⨂ is the tensor product and  () ( = 1,3) are the principal stresses in this reference system.More details can be found in Giry et al. (2011).Finally the expression of the damage variable can be written: A and A are the material parameters which can be interpreted as brittleness terms, respectively, in tension and in compression.
The uniaxial responses in traction and compression obtained with this model are presented in figure 1.There is a good reproduction of the experimental behavior obtained by Terrien (1980) and Ramtani (1990).The boundary displacement field is given on Γ .To simplify, we assume that  is clamped on Γ .On the boundary Γ , a density of forces denoted ℱ is applied.Ω is also submitted to volume forces of density  (figure 2).The notation  represents the unit outward normal vector on Γ.The constitutive relation being written in the expressions (2-a) and ( 14) such as for all point of Ω the strain tensor  =   (  +   ) is assumed to be linearized. and  are respectively the divergence and gradient operators of tensor and vector valued functions.

The approximated problem
Finite element method in displacements consists in finding in Ω the displacement field denoted by  and satisfying: The kinematic conditions: ∀  of F.E type we have: The equilibrium equations: the constitutive relations (2-a) and ( 14) which can be rewritten as: We consider the triplet (, , ) fields on Ω as the exact solution of the reference problem settled above in section (3.1), and the triplet ( ,  ,  ) fields as the finite element solution of the present approximated problem (18) to (20).If we compare these references and approximated problems, we can observe that both the finite element and exact solutions verify the kinematic conditions on Γ .However, the approximated stress field  obtained in (20-a) does not verify the equilibrium equation ( 16): in the displacement finite element method, the approximation is taken over the equilibrium equations.On the other hand, depending on tensor variables  ,   and  −  the computed solution  is also different from the exact solution d.

The exact definition of the error
By considering the triplets cited above (, , ) and ( ,  ,  ) defined on Ω, let set the quantity: with  a non-zero positive and constant parameter that has the same dimension as the error energy norm ‖ −   ‖ , defined as: is a four order symmetric tensor of elasticity .By definition, the quantity  represents the exact definition of the error associated to the obtained solution (  ,   ,   ) field.It is trivial that for any triplet of finite element solution the error  is equal to zero if and only if (  ,   ,   ) is the exact solution of the reference problem.In addition to the positive property of the error energy norm defined with the first term, the two terms of the error being squared, it is easy to observe that for any triplet of finite element solution the defined error  is either positive or equal to zero.
By definition, this error is composed of two terms that can be named respectively by stress and damage errors and computed separately.Indeed, for the first term, in order to compute the error energy norm, we consider the elasticity tensor as constant even the presence of damage.This dissociation between ó and d can also be justified by the two different state equations of the theory of the thermodynamics of irreversible processes written in (2-a) and (2-b) of the model description in section 2.

The developed error indicator
The exact solution is generally unknown.The idea consists to calculate an indicator of the error over the finite element solution by substituting the unknown exact stress  and damage  fields in ( 21) by a stress   and damage  fields with best rather properties of regularity and continuity between elements.Indeed, it is well known that the finite element solution is not smooth and regular enough compared to the exact solution.On the other hand, the choice of piecewise linear interpolation of displacement gives discontinuous solutions stresses (or other gradient).In addition it was found in the residual estimator works that for linear interpolation the term describing the discontinuity (i.e.jump) stresses predominates over the term of the residual equilibrium.Therefore, under these conditions, the error can be only estimated from the term of the jump.This idea is at the origin of the approach developed by Zienkiewicz andZhu (1987 and1992).This indicator consists in building from the finite element solution, a new stress field of high degree of smoothing and which is assumed to be as near as possible to the exact solution.This idea is completed by using some super convergence properties of the integration points in which the stresses are approximated with high accuracy (Hinton and Campbell, 1974).Thereafter, the error can be defined by the energy standard difference between this high degree solution and the finite element one and estimated by substituting the first term of (21) ‖ −   ‖ , by ‖  −   ‖ , such as: In addition, the same principle can be applied to the exact damage field  which is replaced by a smoothed finite element damage  so that the second term ∫ ( −  )  of the equation ( 21) is replaced by ∫  −   ℎ 0 ≤  ≤ 1.This approach is motivated by the fact that not only the internal damage values  are natively calculated at the integration points, but also linked to stress tensor  by the constitutive relations such as the equation (2-a) of the present damage model of RICRAG.
Finally an approximation of the error which will be denoted by   can be computed by the following expression: In order to obtain the smoothed stress   and damage  fields on Ω , several smoothing methods exist in the literature.It should be noted that all these techniques are based on least squares minimizations.Respectively called global ZZ 1 and local ZZ 2 , these two methods developed by Zienkiewicz and Zhu (1987) and Zienkiewicz and Zhu (1992) seem to be actually among the best known.The ZZ 1 method consists to resolve for each stress component a global system of equations for the entire total degree of freedom nodes number of the mesh.However compared to ZZ 2 this method is not only less precise but also quite expensive in terms of CPU time calculations.For its part, the local ZZ 2 method also called super convergent patch recovery or SPR technique produced a significant improvement of performance and accuracy.Indeed the smoothed stress fields are obtained by the resolution for each component of linear equation systems written on patches established by elements fully connected at corner nodes.Although equation systems are also resolved locally by patches this method remains expensive in terms of CPU time costs.In particular, this technique would be more expensive in the cases of meshes using a large number of elements and especially for the three-dimensional computations.This can be a major drawback.
On the other hand, the classical smoothing procedures developed by Hinton and Campbell (1974) consists in three different smoothing techniques: one global and two locals.These two locals are described as function smoothing and discrete smoothing.For the first, it is assumed that a smoothed function is least squares fit to its corresponding unsmoothed function.With the second local discrete smoothing method, it is assumed that  , ( ,  ) is a least squares fit to selected values of  ( ,  ) at the Gaussian integration points s in which the stresses are approximated with high accuracy.Although it is an oldest technique, this method seems not only giving good results in terms of accuracy but also it especially offers a less expensive CPU time of calculations.This remark is demonstrated by the formulation below and confirmed by Delmas (2011) in the framework of tests carried out with the ASTER finite element code developed by EDF in France.These are the reasons that motivated our choice for the implementation of this discrete smoothing local technique with the present error indicator.
Thus within the meaning of least squares and for the case of 2D space distribution that we limit ourselves for the present explanations, the tensors  (, ) can be smoothed by a polynomial function of order p such as for each component denoted by the index i, the smoothed stress field component  , (, ) can be written as: where  contains appropriate polynomial terms of expansion, and   the vector of unknown parameters.For each component , the vector   is determined by the resolution of system of equations, which is obtained by the minimization of the functional: The values of  (, ) being known at the Gauss points, the minimum will be reached if and only if: Moreover the smoothed stresses being the unknowns taken at the nodes j (j =1,n), thus for any point M(x, y) within any element,  , (, ) may be obtained by interpolation using the shape functions as: (, ) is the shape function associated with the node j on the considered finite element,  , ( ,  ) is the unknown value of the smoothed stress component  at the node   is the number of nodes per element.By adopting for  (, ) the same interpolation polynomial function with the same order as in (25) the condition ( 27) can be rewritten as follows: The discrete local smoothing method consists to substitute the functional χ by the following summation written in the reference element: ( ,  ) are the parametric coordinates of the gauss point  ( ,  ) is the known unsmoothed stress component value at the gauss point s  ,  ,  is the unknown smoothed stress component value at the node j  ( ,  ) is the shape function value at the gauss point s.
: is the total number of gauss point in the element  is the number of nodes per element.
The minimization condition leads to the following system: which can be written in a matrix form as: Being respectively square ( x ) and rectangular ( x npg) the matrices  and  are independent of the current element E of the mesh.Therefore they can be calculated only once in the reference element.This is the reason that allows the calculations to be less inexpensive in CPU time calculation.This is also the reason of our choice for this method.
However, it should be noted that the tensor components values σ , for ( i = 1,3) are computed with the resolution of system (32) over each element.Thus unlike the other smoothing approaches, the present method does not produce unique values of stresses at nodal points and therefore nodal averages must be calculated.Despite the calculation of nodal averages seems causing a priori loss of precision.However the numerical tests performed by Hinton and Campbell (1974) showed that this method of local stress smoothing with nodal averaging produces results in good agreement with theoretical ones.
Once the smoothed nodal component stress values obtained, smoothed stress field can be obtained over the entire mesh by the standard FE interpolation such as written in (28).
By the same manner  (, ) can be recovered from the gauss point values of  ( ,  ) by minimizing  in the reference element: However, it should be noted that some corrections must sometimes be made when the smoothed damage field  presents some values slightly greater than 1.This correction consists in normalizing this field  to 1.While this may be a disadvantage, nevertheless, we have noticed that this normalization of  does not affect the quality of the smoothed solution and its distribution.
Moreover, we associate to the error measure in (24) the global relative error defined by: Thus  consists of a global accuracy that allows evaluating the global quality of the finite element solution.Consisting on the sum of two independent terms, the particularity of this indicator is that  can be decomposed into two types of error.This constitute an advantage.Indeed, the first term: which will be called "the relative stress error" is the error carried by the stress field solution.The second term is noted  , and called by "the relative damage error".It consists on the error due to the solution represented by the scalar field of damage solution.In practical situations,  is an element of the mesh descretization associated to Ω.The local contributions allows obtaining information concerning the errors located on the structure.By construction, one has:

CONVERGENCE ANALYSIS
After having developed the present error indicator, the knowledge of its convergence properties being necessary, a benchmark numerical test will be carried out.It consists in the single edge notched beam called SENB concrete test (figure 3).Consisting in an anti-symmetrical four-point loading beam, this test is widely used in the literature.To name a few, we can cite, for example, the experimental and computational works performed by Schlangen (1993), Carpinteri et al. (1993) and the numerical contributions of Geers et al. (2000) and recently Gerasimov et al. (2015).In order to perform some comparisons, it should be noted that the adopted geometry, loads and supports are the same as that used by Huerta and al. (2002) and Rodríguez-Ferran and Huerta, (2000).The applied loads F1= 10F/11 and F2= F/11 are illustrated in this figure 3 such as F increases incrementally from (0) zero to (138F) and F= 0.7 KN.This convergence study will be essentially performed with the h-refinement method, which consists in uniform and successive subdivisions meshes as shown in figure 4. For this purpose the QUA4 four nodes quadrilateral element type will be used.
Latin American Journal of Solids and Structures, 2018, 15(9), e104 9/25 Moreover it should be noted that the implementation of this error indicator is achieved using the Cast3m finite element code developed by the CEA in France.Indeed a plane stress mode analysis is performed under the assumptions of small deformations with quasi-static incremental loading.As mentioned in the introduction, the tests are carried out with the use of the RICRAG nonlocal damage model (Richard and al, 2010) and by adopting the following material parameters of table 1.These parameters are also explained in the section 2 above.Poisson's ratio í

Brittleness in compression AInd
Tensile strength Ft 28 0.2 1.0 x 10 -2 1.0 x 10 -4 3.6 x 10 +6  or error cards.The first concerns the representation of the element contributions to the relative global error ε or global error card (see figure 6 for example).The second card represents the elements contributions to the stress error ε or stress error card.Finally the third shows the contribution of the elements to the damage error  or damage error card (see figure 7-a, figure 7-b, and figure 7-c) By observing the following figures (7-a) to (7-c) it is easy to note that the global errors and damage errors contributions are particularly marked in the large gradient damage zone located in the vicinity of the boundaries of the damaged zones.This result is in good agreement with the error distribution card obtained by Rodríguez-Ferran and Huerta (2000) who used a residual type error estimator based on local computations.the same loading intensity, the damaged area of the most refined mesh proves to be most extended and developed (figure 8-c) than the least refined one (figure 8-a).It should be noted that this observation was also relieved with the simulation of the two others nonlocal damage models of Mazars and MVM or modified Von Mises model.The evolution curves of the three errors illustrated in figure 09 below show that the damage error ε is relatively higher at the first steps following the appearance of damage.However, this error tends to decrease with the increase of loading.Therefore, according to (38), the relative global error  curve is also quite higher with the first steps and decreases substantially by getting closer to the relative stress error  .The latter presents a quasiconstant allure.This damage error diminution can be interpreted by the fact that numerical damage distribution field obtained at first is very irregular and tends to become less irregular after its propagation.
On the other hand, the constant level of the stress error in Figure 9, can be interpreted essentially by the concept of regularity of the stress field solutions and more precisely by the obtained stress gradient fields solutions that remains unchanged for all these increments.with  the approximated error provided by the present indicator such as established in ( 24) and  the exact error computed with respect to the exact solution.However, since this later is not available,  is obtained by using a reference solution corresponding to a very refined mesh procedure.Moreover, the calculation of the absolute error  requires the knowledge of the constant  introduced in (21).
However, this parameter remains unknown so that one can only compute the relative errors: å, åó and åd (where this constant  is simplified).Thus, this does not allows computing directly this approximated error  .To overcome this difficulty, the proposed solution consists to define simultaneously and separately two effectivity indexes.The first index is noted γó .It corresponds to the stresses error in energy norm considering the recovered stresses and it can be defined by: where  and  are respectively the smoothed and non-smoothed reference solutions cited above and    are the smoothed and non-smoothed obtained solutions.
On the other hand, the second index noted γd corresponds to the effectivity index of the error due to the damage field solution compared to the reference one.It can also be defined by: where  and  are respectively the smoothed and non-smoothed damage reference solutions cited above and    are the smoothed and non-smoothed damage obtained solutions.
Therefore, the efficiency of the present error indicator would be sufficiently accepted if both these two effectivity indexes are simultaneously close to unity or at least belonging to intervals whose lower bounds noted respectively C1s and C1d and upper bound noted C2s and C2d are generally close to unity and such as: C1s ≤ γó ≤ C2s and C1d ≤ γd ≤ C2d.
By applying the previous relations ( 43) and ( 44) to a reference solution obtained with a mesh of 60336 four nodes quadrilateral elements and 60985 nodes, several values of γó and γd are depicted in table 2 below under different loading levels.These results are derived from the SENB simulation test previously treated in section 4. Indeed, it should be recalled that the two applied loads F1 and F2 increased from zero to a maximal intensity value of 138.F (with F=0.7KN).For the chosen loading levels showed in table 2, the effectivity index γó values are constant and equal to 1.26.This shows that the stress errors in energy norm considering recovered stresses are overestimated by about 26%.On the other hand, for the γd index values, except the last loading steps, where γd equals 1.83, the other γd values are between 0.89 and 1.24 corresponding to a loading interval with lower and upper bounds respectively equals to 40F and 130.F (i.e.covering 92.1% of the total damage loading steps).Indeed, damage mechanism being started from a loading intensity of 36.8F until 138F, this interval corresponds to a percentage loading ratio of 92.1% equal to (130-36.8)/(138-36.8)%).In this interval, the γd values are closest to unity than the γó ones.However, beyond the value of 135F (ie for the last remaining 8%), the efficiency of this error indicator seems to be limited.

CHARACTERISTIC LENGTH INFLUENCES
In the framework of the nonlocal approach, the characteristic length  is considered among the most important parameters.Indeed, its influence on the results can be analyzed under two main aspects.The first consists in the influence on both damage and error distributions.The second aspect considers the  influence on the error convergence rates.
For the first aspect, we consider four cases of simple entailed beams simulated with a uniformly refined mesh of elements of size h equal to 2.5 mm and with the same properties as in section 4 above but with different values of   respectively equals to 10, 20, 30 and 40 mm.The obtained damage distributions for each beam under the same level of loading equal to 138.F, are illustrated in cases a) to case d) in figure 15 below.By comparing these figures a) to d), significant differences can be observed.Indeed, for the two lower values of   (10 and 20 mm), the damaged arc zones are completely formed.Nevertheless, a larger bandwidth for the case  equal to 20mm is observed.For the other values (lc equal to 30 and 40 mm), the arched damaged zones are not completely formed yet and significant differences in their shapes can be remarked.By comparing these figures, one can conclude that for this damage model of Ricrag, as  increases, as both the apparition and development of damage zones are delayed.
In addition, since the damage error card is closely linked to damage field distribution (see section 4 above), the influence of  on the distribution of this damage error distribution can be established.This is illustrated in Figure 16 and Table 3 below.Furthermore, for each value of  , the loading intensities corresponding to the first appearance of damage is shown in table 4.These results are in accordance with the equations ( 10) and ( 11) of the damage model established in section 2. Indeed, the damage threshold surface is written in function of the nonlocal damage energy released rate  .This latter depends directly on the characteristic length parameter  .It should be also noted that  is generally determined by identification with experimental results.For the second aspect considering the influence of  on the error convergence rate, five cases corresponding to five different values of  (10, 20, 30, 40 and 50 mm), have been studied.The results are presented in figure 17 below containing five curves and representing, for each value of  , the evolution of this convergence rate p as a function of the loading intensity.It can be observed that the strong variation of p previously shown in section 4 before is also confirmed for the three first values (10, 20 and 30 mm).However, for the two others, although this variation become less strong, weak values of p are obtained.These weak values constitute also another major difficulty for the mesh adaptation phase.Indeed, it is well known that, below a certain value of p approximately equal to 0.6, mesh adaptation procedure becomes rude and requires software mesh generator with high performances.

MESH ADAPTATION
The mesh adaptation procedure consists to both guarantee a certain level of precision to the finite element user and minimizing the computational time CPU costs.In this framework, the most frequently used procedure hgeneration will be used: the size and topologies of the elements are modified while the type of finite element shape functions on the different meshes remains the same.Indeed, a mesh T* is said to be optimal (Ladeveze et al. 1986 where ℎ denotes the size of E and ℎ * represents the size that must be imposed to the elements of  * in the region of E in order to ensure optimality.The computation of these coefficients  uses the rate of convergence of the error p explained in section 4. Therefore, to compute these coefficients  , a technique detailed in Coorevits et al. (1994) that automatically takes into account the steep gradient regions is used.The mesh  * is generated by an automatic mesher that would be able to accurately respect the obtained card sizes.For this purpose, a high level of performance and precision are required.For the present computations, the mesher implemented in Cast3m has been used.Moreover, by taking into account the difficulty showed in the previous section 4, we will limit ourselves to a mesh adaptation corresponding to one arbitrary chosen loading step.On the other hand, for the frequent cases of non-uniform mesh size elements (with h variable in the mesh) such as in our case (figure 13), the determination of the global error convergence rate p is more complex and difficult.A partial solution consists, in a first phase, in adapting the mesh with respect to the stress error card by adopting the predetermined constant value of p our case p = 0.59 for the SENB test, see figure 12).The second phase is most difficult.It consists in adapting the mesh with respect to the damage error card whose corresponding value p is (in this case) not only available but also non precise.Indeed, in the absence of available methods for determining p for the case of non-uniform meshes, we limited ourselves to manually determine p by successive remeshing iterations.This method is not only less accurate but also very expensive in terms of CPU time.Its improvement must be the subject of future researches.
The SENB first example of optimized calculations will be performed and presented in the following section 5.1 (figure 18).On the other hand, the section 5.2 consists on a mesh adaptation performed on double edge notched beam DENB test as treated in example 02.

Example 01: The SENB test (Simple edge notched beam)
This first example consists in the beam illustrated with its boundary conditions in figure 2. The same material properties as showed in table 1 and the same loads as explained in section 4 above are adopted.The characteristic length  is taken equal to 20mm.By applying the present indicator to an initial and coarse mesh made of 765 four nodes elements, and 789 nodes (figure 18), a global accuracy å of 12.20% is firstly obtained.The deformed shape and a damage field distribution are also represented in figure 19.The elements contributions to stress, damage and global errors of this coarse mesh can respectively be seen on figures 20, 21 and 22.For a prescribed error of å0 =5%, the resolution of the problem (45) gives a prescribed element size card as showed in figure 23.This card size is prescribed to the mesher used by the finite element code Cast3m.An optimized mesh made of 5958 four nodes elements and 5858 nodes (figure 24) is provided.The developed indicator provides also a global accuracy å = 4.86% with a relative stress error åó = 3.39% and a relative damage error of åd 3.48%.The elements contributions to global, stress and damage errors of this adapted mesh can respectively be seen on figures 25, 26 and 27.The second example is also illustrated with its boundary conditions and applied loads in figure 28.The same material properties as showed in table 1 are adopted.For this test, the characteristic length  is taken equal to 10mm.The initial mesh is made of 842 elements (figure 29).The smoothed (reference solution) and non smoothed (finite element solution) damage field distributions are also represented in figure 30.The elements contributions to stress, damage and global errors of this coarse mesh can respectively be seen on figures 31, 32 and 33.For a prescribed error of å0 =5%, the resolution of the problem (45) gives the element size card as shown in figure 34.This prescribed card exploited by the Cast3m mesher provides an optimized mesh made of 6110 four nodes elements and 5903 nodes (figure 35) with a global accuracy of å = a relative stress error of åó = 2.72% and a relative damage error of åd = 4.14%.The global, stress and damage error cards of this adapted mesh can respectively be seen on figures 36, 37 and 38.Based on local and least squares smoothing techniques, a new error indicator with relatively less difficulties for its implementation and with a significant reduced CPU time computations is developed.For the errors distribution fields the obtained results are in good agreement with those obtained by the use of residual method.However, for the global relative error, the obtained values are significantly larger than those provided by residuals.This can be explained by the fact that this global error is composed of two others (damage and stress errors) which have the same orders of magnitude.Therefore, this explain also the precision requirement of 5% instead of 2% for mesh adaptation tests.On the other hand, the convergence analyses highlighted a strong variation of the damage error convergence rate with step loading.This difficulty is most accentuated in the case of non uniform mesh element sizes.This constitutes the major difficulty to obtain an automatic step-by-step mesh adaptation.Nevertheless, for a chosen step and by limiting ourselves to determine manually the damage error convergence rate p, optimized calculations have been performed.However, it should be noted that in such case, adaptation procedure is not only less accurate but also expensive in terms of calculation time.Its improvement must be the subject of future researches.

Figure 2 :
Figure 2: Setting of the problem.
, (, ) = .  = ∑  ,  (, )  = 1 to 3 for 2D problems (25)Latin American Journal of Solids and Structures, 2018, 15(9), e104 7/25 relation between ,  and  can be written in the following form: =  +  (38)On the other hand, let us consider  a part of Ω .Then, we define the local contribution of  to the relative global error  by the quantity  , such as:

Figure 7 :
Figure 7: From left to right respectively: elements contributions to global error , to stress error  and to damage error  cards in the vicinity of the notch.

Figure 8 :
Figure 8: damage distributions in the vicinity of the notch for different meshes

Figure 10 :
Figure 10: Evolution of ε, εσ and εd in function of the total number of degree of freedom

A
Figure 15: Damage distribution of SENB test submitted to the same intensity of loading (138.F)and simulated with different values of  (=10, 20, 30 and 40 mm)

Figure 17 :
Figure 17: damage error convergence rate p in function of the loading intensities for different values of the characteristic length

Figure 18 :
Figure 18: initial mesh with 765 elements and 789 nodes

Figure 30 : 25 Figure 31 : 25 Figure 34 :
Figure 30: Non smoothed (on the left) and smoothed (on the right) damage field distributions ,  ) is the known unsmoothed stress component value at the gauss point s   ,  is the unknown smoothed stress component value at the node j This minimization leads to the same equation system type with the same matrices  and  which have already been calculated in (31): Latin American Journal of Solids and Structures, 2018, 15(9), e104 8/25  (

Table 1 :
Material parameters

Table 2 :
Effectivity indexes for different loading levels

Table 3 :
relative error values of SENB tests submitted to the same loading level of 138F and for different values of

Table 4 :
loading intensity values corresponding to the first appearance of damage for different values of  .
,  and the global, stress and damage error cards -Determination of the optimized mesh  * provided by the calculations of a prescribed size card obtained by using the method developed by Coorevits et al. (1995) -and a second computation of the relative errors on the new mesh  * is carried out.If the obtained relative global error å is greater than  , this procedure is repeated as far as a precision close to  is reached.It is worth noting that this prescribed size card is determined by the computation of size modification coefficients  on each element E of the mesh T: