Acessibilidade / Reportar erro

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

Abstract

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.

Keywords
FE analysis -concrete- nonlocal damage; error indicator; least squares smoothing; error convergence rate; step mesh adaptation

1 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 (1958Kachanov, L. M. (1958). On the time to failure under creep conditions. Izv. AN SSSR, Otd. Tekhn. Nauk 8: 26-31.) and developed by Rabotnov (1969Rabotnov, YN. (1969) Creep problem in structural members. North Holland, Amsterdam.). 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 (1984Mazars, J. (1984). Application de la mécanique de l’endommagement au comportement non linéaire et à la rupture du béton de structure, Thèse de doctorat d’état de l’Université Paris VI.) and Mazars et al (2014)Mazars J., Hamond F., Grange S., (2014) A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings, Materials and Structures.. 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. (1995De Vree, J.H.P., Brekelmans, W.A.M., Gils, M.A.J., (1995). Comparison of nonlocal approaches in continuum damage mechanics. Computers and Structures 55: 581-588.) seems to be very used in the literature. For our simulations, the model named RICRAG that is developed by Richard et al. (2010Richard, B, Ragueneau, F., Cremona, C., Adelaide L. (2010) Isotropic continuum damage mechanics for concrete under cyclic loading: Stiffness recovery, inelastic strains and frictional sliding. Engineering Fracture Mechanics, Elsevier, 77 (8):1203-1223.) 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, 1986Bazant Z.-P., (1986).”Mechanics of distributed cracking” Appl Mech Rev, ASME, 39: 675-705.). 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, 1988Needleman, A. (1988) Material rate dependence and mesh sensitivity in localization problems Comput. Meth. Appl. Mech. Engrg. Vol. 67 p.69-85.). To eliminate this mesh dependence, different localization limiters have been developed. The method performed by Pijaudier-cabot and Bazant (1987Pijaudier-Cabot G, Bazant Z -P. (1987) “Nonlocal damage theory” Journal of Engineering Mechanics, ASCE, vol 113, n°10:1512-1533.) 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. (2011Giry C., Dufour F., Mazars J. (2011) Stress based nonlocal damage model, International Journal of Solids and Structures, 48: 3431-3443.). 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 (1978Babuska, I., Rheinboldt, W.C., (1978). Error estimates for adaptive finite element computations. SIAM J. Numer. Anal. vol. 15, n°4: 736-754.) and based on the residual of the equilibrium equations; the estimators linked to the smoothing of finite element stresses performed by Zienkiewicz and Zhu (1987Zienkiewicz, O.C., Zhu, J.Z., (1987). A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. for Num. Meth. in Engng. 24, 337-357. and 1992); and the estimator based on the errors in the constitutive relation developed by Ladeveze (1975Ladeveze, P., (1975). Comparaison de modèles de milieux continus, Ph.D. thesis, Paris VI, university.) and Ladeveze and Leguillon (1983)Ladeveze, P., Leguillon, D. (1983). Error estimate procedure in the finite element method and applications. SIAM J. Num. Anal. Vol. 20 N° 3: 483-509.. A review of different a posteriori estimators can be found in Zienkiewicz (2006)Zienkiewicz, O.C. (2006). The background of error estimation and adaptivity in finite element computations. Comput. Methods Appl. Mech. Engrg. 195: 207-213. and Verfürth (2013Verfürth R. (2013) A Posteriori Error Estimation Techniques for Finite Element Methods. Oxford University Press: Oxford.). 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 (2000Rodríguez-Ferran, A., Huerta, A. (2000) Error estimation and adaptivity for nonlocal damage models. Int. J. Solids Struct. 37:7501-7528.). Furthermore, this method has been improved by the recent contributions of Gerasimov et al. (2015Gerasimov, T; Stein E. and Wriggers P. (2015) Constant-free explicit error estimator with sharp upper error bound property for adaptive FE analysis in elasticity and fracture Int. J. Numer. Meth. Engng; 101:79-126.). On the other hand, Comi and Perego (2002Comi, C., Perego, U., (2002) Finite element strategies for damage assessment up to failure. In Proceedings of the 6th National Congress SIMAI, Chia Laguna, Italy.) adapted the error in the constitutive relations cited above to a nonlocal damage model they developed (Comi and Perego, 2001Comi, C., Perego, U., (2001) Numerical aspects of nonlocal damage analyses. Revue européenne des éléments finis, 10:227-242.). 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 (1974Hinton, E., Campbell, J.S. (1974). Local and global smoothing of discontinuous finite element functions using a least squares method. Int. J. Numer. Meth. Engrg. 8: 461-480.) and Hinton et al. (1975)Hinton E., Scott F.C., Ricketts R.E. (1975) - Local least squares stress smoothing for parabolic isoparametric elements - Int. J. for Num. Meth. in Eng. 9: 235 - 256.. 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 well-known applications of simple and double edge notched beams respectively called SENB and DENB tests.

2 NONLOCAL DAMAGE FORMULATION

In this paper, the RICRAG model developed by Richard et al. (2010Richard, B, Ragueneau, F., Cremona, C., Adelaide L. (2010) Isotropic continuum damage mechanics for concrete under cyclic loading: Stiffness recovery, inelastic strains and frictional sliding. Engineering Fracture Mechanics, Elsevier, 77 (8):1203-1223.) 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)Cast3m. (2018). www-Cast3m.cea.fr
www-Cast3m.cea.fr...
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:

ρΨ = 1 2 { κ 3 ( ( 1 d ) ε kk + 2 ε kk + 2 ) + 2 ( 1 d ) με ij D ε ij D + 2 ( ε ij D ε ij π ) ( ε ij D ε ij π ) + γα ij α ij } + Η ( z ) (1)

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. εij is the second order total strain tensor and εijD is the second order deviatoric total strain tensor such as εijD= εij13εkkδij where δij 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 (2dμ(εijDεijπ)(εijDεijπ)) represent the sliding mechanism by contact with friction where εijπ is the second order sliding tensor.

Finally, the last terms (γαijαij) and Η(z) 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, αij 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, 2010Richard, B, Ragueneau, F., Cremona, C., Adelaide L. (2010) Isotropic continuum damage mechanics for concrete under cyclic loading: Stiffness recovery, inelastic strains and frictional sliding. Engineering Fracture Mechanics, Elsevier, 77 (8):1203-1223.), the first state equations can be deduced such as:

σ ij = ρΨ ε ij = κ 3 ( ( 1 d ) ε kk + ε kk + ) δ ij + 2 ( 1 d ) με ij D + 2 ( ε ij D ε ij π ) (2-a)

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:

Y = ρ Ψ d = 1 2 { κ 3 ε k k + 2 + 2 μ ε i j D ε i j D + 2 μ ( ε i j D ε i j π ) ( ε i j D ε i j π ) } (2-b)

A threshold surface, denoted fd is introduced:

f d = Y ¯ ( Y 0 + Z ) such as d ˙ = λ ˙ d . f d Y ¯ (3)

where λ˙d is a Lagrange multiplier and Z a thermodynamic force associated to isotropic hardening:

Z = ρ Ψ z = Η ( z ) z (4)

Y¯ denotes the energy variable driving damage and Y0 an initial threshold. This variable Y¯ is defined through the decomposition between direct extension mechanisms (tension) and induced extension mechanisms (compression):

Y ¯ = β Y D i r + Y I n d (5)

with

Y D i r = 1 2 ε i j D i r E i j k l ε k l D i r (6)

Y I n d = 1 2 ε i j I n d E i j k l ε k l I n d (7)

where Eijkl 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:

ε ij Dir = ε ij + H ( ε ij + σ ij + ) (8)

ε i j I n d = ε i j ε i j D i r (9)

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. (2011Giry C., Dufour F., Mazars J. (2011) Stress based nonlocal damage model, International Journal of Solids and Structures, 48: 3431-3443.). 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 Y¯ is replaced by the nonlocal one Y¯nl such as the damage threshold surface becomes:

f d n l = Y ¯ n l ( Y 0 + Z ) a n d d ˙ = λ ˙ d . f d Y ¯ n l (10)

Y ¯ n l ( x ) = Ω ( x ) W ( x s ) Y ¯ ( s ) d s Ω ( x ) W ( x s ) d s (11)

W(xs) is the weight or also called Gauss distribution function defining the interaction between the considered point located with x and the neighboring points located with s inside the volume of the structure Ω(x) such as:

W ( x s ) = exp [ ( 2 x s l c . ρ ( x , σ prin ( s ) ) ) 2 ] (12)

lc is a characteristic length which limits the domain of neighboring points and ρ(x,σprin(s)) is a scalar which describes the interaction between points. It depends on the point location x and σprin(s) the principal stress state on the surrounding point s expressed in its principal stress reference system (s, U1(s),U2(s),U3(s)) such as:

σ p r i n ( s ) = i = 1 3 σ i ( s ) ( U i ( s ) U i ( s ) ) (13)

is the tensor product and σi(s) (i=1,3) are the principal stresses in this reference system. More details can be found in Giry et al. (2011Giry C., Dufour F., Mazars J. (2011) Stress based nonlocal damage model, International Journal of Solids and Structures, 48: 3431-3443.). Finally the expression of the damage variable can be written:

d = 1 1 1 + ( A Dir H ( ε ij + σ ij + ) + A Ind ( 1 H ( ε ij + σ ij + ) ) ) ( Y ¯ nl Y 0 ) (14)

ADir and AInd 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 (1980Terrien, M. (1980). Emission acoustique et comportement mécanique post-critique d’un béton sollicité en traction. BLPC, vol 105, p 65-72.) and Ramtani (1990Ramtani, S. (1990). Contribution à la modélisation du comportement multi-axial du béton endommagé avec description du caractère unilatéral. Ph.D. thesis, Paris VI, university.).

Figure 1
Ric-Rag uniaxial behaviors model compared with experimental results (Richard et al., 2010Richard, B, Ragueneau, F., Cremona, C., Adelaide L. (2010) Isotropic continuum damage mechanics for concrete under cyclic loading: Stiffness recovery, inelastic strains and frictional sliding. Engineering Fracture Mechanics, Elsevier, 77 (8):1203-1223.)

3 THE DEVELOPED DAMAGE ERROR INDICATOR

3.1 The reference problem

We consider a domain Ù with boundary Γ such as Γ= ΓD U ΓF where ΓD and ΓF are two open disjoints parts. The boundary displacement field is given on ΓD. To simplify, we assume that Ωis clamped on ΓD. On the boundary ΓF, a density of forces denoted F is applied. Ω is also submitted to volume forces of density f (figure 2). The notation n represents the unit outward normal vector on Γ.

Figure 2
Setting of the problem.

The problem consists in finding the scalar damage field d, the displacement vector field u and the stress tensor field σ in Ω satisfying the previous equations (1) to (14) and the kinematics (15) and equilibrium conditions (16) and (17) as follows:

u = 0 o n Γ D (15)

d i v σ + f = 0 i n Ω (16)

σ n = F o n Γ F (17)

The constitutive relation being written in the expressions (2-a) and (14) such as for all point of Ω the strain tensor ε= 12(grad u+gradtu) is assumed to be linearized. div and grad are respectively the divergence and gradient operators of tensor and vector valued functions.

3.2 The approximated problem

Finite element method in displacements consists in finding in Ω the displacement field denoted by uh and satisfying:

The kinematic conditions: uh of F.E type we have:

u h = 0 o n Γ D (18)

The equilibrium equations:

Ω [ T r ε ( u h ) σ ( v h ) ] d Ω + Ω f T v h d Ω + Γ F F T v h d Γ = 0 v h V w i t h V = { v d e f i n e d a n d r e g u l a r e n o u g h o n Ω s u c h a s v = 0 o n Γ D } (19)

the constitutive relations (2-a) and (14) which can be rewritten as:

σ h ij = κ 3 ( ( 1 d h ) ε h kk + ε h kk + ) δ ij + 2 ( 1 d h ) με h ij D + 2 d h μ ( ε h ij D ε h ij π ) (20-a)

d h = 1 1 1 + ( A Dir H ( ε h ij + σ h ij + ) + A Ind ( 1 H ( ε h ij + σ h ij + ) ) ) ( Y ¯ h nl Y h 0 ) (20-b)

We consider the triplet (σ, u, d) fields on Ω as the exact solution of the reference problem settled above in section (3.1), and the triplet (σh, uh, dh) 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 ΓD. However, the approximated stress field σh 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 σh, εh and (Y¯hnlYh0) the computed solution dh is also different from the exact solution d.

3.3 The exact definition of the error

By considering the triplets cited above (σ, u, d) and (σh, uh, dh) defined on Ω, let set the quantity:

e = [ σ σ h σ , Ω 2 + C d . Ω ( d d h ) 2 ] 1 / 2 (21)

with Cda non-zero positive and constant parameter that has the same dimension as the error energy norm σσhσ,Ω defined as:

σ σ h σ , Ω = ( Ω ( σ σ h ) T C 1 ( σ σ h ) ) 1 / 2 (22)

C is a four order symmetric tensor of elasticity . By definition, the quantity e represents the exact definition of the error associated to the obtained solution (σh, uh, dh) field. It is trivial that for any triplet of finite element solution the error e is equal to zero if and only if (σh, uh, dh) 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 e 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.

3.4 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 d fields in (21) by a stress σ˜h and damage d˜h 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 and Zhu (1987Zienkiewicz, O.C., Zhu, J.Z., (1987). A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. for Num. Meth. in Engng. 24, 337-357. and 1992). 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, 1974Hinton, E., Campbell, J.S. (1974). Local and global smoothing of discontinuous finite element functions using a least squares method. Int. J. Numer. Meth. Engrg. 8: 461-480.). 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) σσhσ,Ω by σ˜hσhσ,Ω such as:

σ ˜ h σ h σ , Ω = ( Ω ( σ ˜ h σ h ) T C 1 ( σ ˜ h σ h ) ) 1 / 2 (23)

In addition, the same principle can be applied to the exact damage field d which is replaced by a smoothed finite element damage d˜h so that the second term (Ω(ddh)2dΩ) of the equation (21) is replaced by Ω(d˜hdh)2dΩ with 0d˜h1. This approach is motivated by the fact that not only the internal damage values d 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 es can be computed by the following expression:

e s = [ σ ˜ h σ h σ , Ω 2 + C d . Ω ( d ˜ h d h ) 2 ] 1 / 2 (24)

In order to obtain the smoothed stress σ˜h and damage d˜h 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 ZZ1 and local ZZ2, these two methods developed by Zienkiewicz and Zhu (1987Zienkiewicz, O.C., Zhu, J.Z., (1987). A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. for Num. Meth. in Engng. 24, 337-357.) and Zienkiewicz and Zhu (1992)Zienkiewicz, O. C., Zhu, J. Z., (1992). The Superconvergent Patch Recovery and adaptive finite element refinement. Computer Methods in Applied Mechanics and Engineering, 101, Iss 1-3: 207-224. seem to be actually among the best known. The ZZ1 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 ZZ2 this method is not only less precise but also quite expensive in terms of CPU time calculations. For its part, the local ZZ2 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 (1974Hinton, E., Campbell, J.S. (1974). Local and global smoothing of discontinuous finite element functions using a least squares method. Int. J. Numer. Meth. Engrg. 8: 461-480.) 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 σ˜h,i(xs,ys) is a least squares fit to selected values of σhi(xs,ys) 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 (2011Delmas, J. (2011) https://www.code-aster.org/V2/doc/v10/fr/man_r/r3/r3.06.03.pdf
https://www.code-aster.org/V2/doc/v10/fr...
) 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 σh(x,y) 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 σ˜h,i(x,y) can be written as:

σ ˜ h , i ( x , y ) = P . a i = j = 1 p a i , j P j ( x , y ) i = 1 to 3 for 2 D problems (25)

where P contains appropriate polynomial terms of expansion, and ai the vector of unknown parameters. For each component i, the vector ai is determined by the resolution of system of equations, which is obtained by the minimization of the functional:

χ i = Ω ( σ h i ( x , y ) σ ˜ h i ( x , y ) ) 2 d x d y (26)

The values of σhi(x,y) being known at the Gauss points, the minimum will be reached if and only if:

χ i a i , j = 0 j = 1 t o p (27)

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, σ˜h,i(x,y) may be obtained by interpolation using the shape functions as:

σ ˜ h , i ( x , y ) = j = 1 n N j ( x , y ) σ ˜ h , i ( x j , y j ) (28)

Nj(x,y) is the shape function associated with the node j on the considered finite element,

σ˜h,i(xj,yj) is the unknown value of the smoothed stress component i at the node j

n is the number of nodes per element.

By adopting for Nj(x,y) the same interpolation polynomial function with the same order as in (25) the condition (27) can be rewritten as follows:

χ i σ ˜ h , i ( x j , y j ) = 0 j = 1 t o n (29)

The discrete local smoothing method consists to substitute the functional χi by the following summation written in the reference element:

χ ˜ i = s = 1 n p g ( σ h i ( ξ s , η s ) σ ˜ h i ( ξ s , η s ) ) 2 = s = 1 n p g ( σ h i ( ξ s , η s ) j = 1 n N j ( ξ s , η s ) σ ˜ h , i ( ξ j , η j ) ) 2 (30)

(ξs,ηs) are the parametric coordinates of the gauss point

σhi(ξs,ηs) is the known unsmoothed stress component value at the gauss point s

σ˜h,i(ξj,ηj) is the unknown smoothed stress component value at the node j

Nj(ξs,ηs) is the shape function value at the gauss point s.

npg: is the total number of gauss point in the element

n is the number of nodes per element.

The minimization condition leads to the following system:

s = 1 n p g j = 1 n N k ( ξ s , η s ) N j ( ξ s , η s ) σ ˜ h i , j ( ξ j , η j ) = s = 1 n p g N k ( ξ s , η s ) σ h i ( ξ s , η s ) k = 1, , n (31)

which can be written in a matrix form as:

M { σ ˜ h i , nodes } = P { σ h i , Gauss } (32)

Being respectively square (n x n) and rectangular (n x npg) the matrices M and P 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 {σ˜hi,nodes} 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 (1974Hinton, E., Campbell, J.S. (1974). Local and global smoothing of discontinuous finite element functions using a least squares method. Int. J. Numer. Meth. Engrg. 8: 461-480.) 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 d˜h(x,y) can be recovered from the gauss point values of dh(xs,ys) by minimizing γ˜i in the reference element:

γ ˜ i = s = 1 n p g ( d h ( ξ s , η s ) d ˜ h ( ξ s , η s ) ) 2 = s = 1 n p g ( d h ( ξ s , η s ) j = 1 n N j ( ξ s , η s ) d ˜ h ( ξ j , η j ) ) 2 (33)

dh(ξs,ηs) is the known unsmoothed stress component value at the gauss point s

d˜h(ξj,ηj) is the unknown smoothed stress component value at the node j

This minimization leads to the same equation system type with the same matrices M and P which have already been calculated in (31):

M { d ˜ h , n o d e s } = P { d h , G a u s s } (34)

However, it should be noted that some corrections must sometimes be made when the smoothed damage field d˜h presents some values slightly greater than 1. This correction consists in normalizing this field d˜h to 1. While this may be a disadvantage, nevertheless, we have noticed that this normalization of d˜h 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:

ε = [ σ ˜ h σ h σ , Ω 2 σ ˜ h + σ h σ , Ω 2 + Ω ( d ˜ h d h ) 2 Ω ( d ˜ h + d h ) 2 ] 1 2 (35)

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:

ε σ = σ ˜ h σ h σ , Ω σ ˜ h + σ h σ , Ω (36)

which will be called “the relative stress error” is the error carried by the stress field solution. The second term is noted εd, and called by “the relative damage error”. It consists on the error due to the solution represented by the scalar field of damage solution.

ε d = ( Ω ( d ˜ h d h ) 2 d Ω Ω ( d ˜ h + d h ) 2 d Ω ) 1 2 (37)

Finally, the relation between ε, εσ and εd can be written in the following form:

ε 2 = ε σ 2 + ε d 2 (38)

On the other hand, let us consider E a part of Ω . Then, we define the local contribution of E to the relative global error ε by the quantity εE, such as:

ε E = [ σ ˜ h σ h σ , E 2 σ ˜ h + σ h σ , Ω 2 + E ( d ˜ h d h ) 2 d E Ω ( d ˜ h + d h ) 2 d Ω ] 1 2 (39)

In practical situations, E 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:

ε 2 = E ε E 2 (40)

4 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 (1993Schlangen, (1993). Experimental and numerical analysis of fracture processes concrete. Ph.D. thesis, Delft University of Technology.), Carpinteri et al. (1993Carpinteri, A., Valente, S., Ferrara, G., Melchiorri, G., (1993). Is mode II fracture energy a real material property. Computers and Structures 48: 397-413.) and the numerical contributions of Geers et al. (2000Geers MGD, De Borst R, Peerlings RHJ. (2000) Damage and crack modeling in single-edge and double-edge notched concrete beams. Engineering Fracture Mechanics; 65:247-261.) and recently Gerasimov et al. (2015Gerasimov, T; Stein E. and Wriggers P. (2015) Constant-free explicit error estimator with sharp upper error bound property for adaptive FE analysis in elasticity and fracture Int. J. Numer. Meth. Engng; 101:79-126.). 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. (2002Huerta, A., Rodríguez-Ferran A., Díez P. (2002) Error estimation and adaptivity for nonlinear FE analysis. Int. J. Appl. Math. Comput. Sci.12:59-70.) and Rodríguez-Ferran and Huerta, (2000Rodríguez-Ferran, A., Huerta, A. (2000) Error estimation and adaptivity for nonlocal damage models. Int. J. Solids Struct. 37:7501-7528.). 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.

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, 2010Richard, B, Ragueneau, F., Cremona, C., Adelaide L. (2010) Isotropic continuum damage mechanics for concrete under cyclic loading: Stiffness recovery, inelastic strains and frictional sliding. Engineering Fracture Mechanics, Elsevier, 77 (8):1203-1223.) and by adopting the following material parameters of table 1. These parameters are also explained in the section 2 above.

Table 1
Material parameters

Figure 3
Single Edge Notched Beam SENB test

Figure 4
Mesh refinement examples of SENB test

After resolution, the damage distribution is shown in figure 5 below. The arched form and the numerical crack trajectories obtained of the damaged area under the notch is identical to those obtained by several authors such as Schlangen (1993Schlangen, (1993). Experimental and numerical analysis of fracture processes concrete. Ph.D. thesis, Delft University of Technology.) Rodríguez-Ferran and Huerta, (2000Rodríguez-Ferran, A., Huerta, A. (2000) Error estimation and adaptivity for nonlocal damage models. Int. J. Solids Struct. 37:7501-7528.), Geers et al. (2000Geers MGD, De Borst R, Peerlings RHJ. (2000) Damage and crack modeling in single-edge and double-edge notched concrete beams. Engineering Fracture Mechanics; 65:247-261.), Jirásek and Grassl (2008Jirásek, M. and Grassl, P. (2008) Evaluation of directional mesh bias in concrete fracture simulations using continuum damage models. Engineering Fracture Mechanics, 75:1921-1943.) and Gerasimov et al. (2015Gerasimov, T; Stein E. and Wriggers P. (2015) Constant-free explicit error estimator with sharp upper error bound property for adaptive FE analysis in elasticity and fracture Int. J. Numer. Meth. Engng; 101:79-126.) and others. In addition, due to the bending effect another damage triangular shape zone can be observed on the same figure 05. This bending effect was also reported in Rodríguez-Ferran and al, (2001)Rodríguez-Ferran A., Arbós I. and Huerta A. (2001): Adaptive analysis based on error estimation for nonlocal damage models.-Revue européenne des éléments finis (Special issue: Numerical Modelling in Damage Mechanics), Vol. 10, No. 2-4: 193-207. who used another damage model: the MVM or the Modified Von Mises Model.

Figure 5
Damage distribution for undeformed and deformed shape (amplified 100 times)

The advantage of the present indicator consists in its ability to dissociate between the error stress field εσ and the error damage field εd. Therefore, this allows the possibility of obtaining three different error representations 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 εd 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 (2000Rodríguez-Ferran, A., Huerta, A. (2000) Error estimation and adaptivity for nonlocal damage models. Int. J. Solids Struct. 37:7501-7528.) who used a residual type error estimator based on local computations.

Figure 6
Contribution to relative global error ε = 4.89% for mesh with 26816 elements and 27249 nodes (h= 2.5 mm)

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

Moreover, for the RICRAG model that we used, the influence of the mesh refinement on the areas of the damaged zones appears clearly by comparing the results in the three illustrations on the figure 8 below. Indeed, for 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 εd 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 quasi-constant 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.

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

Figure 9
Evolution of ε, εσ and εd in function of the loading intensity for a mesh with h= 2.5 mm

Although this diminution is important, however, in order to perform an automatic step-by-step mesh adaptation, the knowledge of the error convergence rate evolution remains essential. Indeed let us recall the well known a priori relation between the relative global error ε and the mesh finite element size h such as:

ε = C . h p (41)

where C is a non-zero positive constant and p the rate of convergence of the error. The exponent p is equal to 1 for regular finite element solutions and remains less than 1 with the presence of singularities or the presence of steep gradients areas. The evolutions of the relative errors ε, εσ and εd in function of the total number of degrees of freedom of the nodes are illustrated in figure 10. On the other hand, for selected values of the applied load intensities, the evolutions of these three errors in function of the size element h of the successive refined meshes are depicted in figure 11, figure 12 and figure 13. These curves illustrate the convergence properties of these different errors ε, εσ and εd.

However, figures 11, figure 12 and figure 13 show also the influence of the load intensity on the error convergence rates p. More precisely, it can be easily noted with the logarithmic scales representation of these figures, that one another most important result consists in the important variation with the loading intensities of the convergence rate p of εd and ε. This important result is also showed for each loading step in figure 14. Indeed, in the framework of nonlinear calculations and with the current requirement to develop an automatic step-by-step mesh adaptation, this variation would constitute a major difficulty because for each step loading, it becomes necessary to compute this error convergence rate p. Hence, this is also quite expensive in terms of CPU time calculations.

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

Figure 11
global error evolutions with the uniform refinement mesh for four different loading values

Figure 12
stress error evolutions with the uniform refinement mesh for different loading values

Figure 13
damage error evolutions with the uniform refinement mesh for different loading values

Figure 14
SENB test error convergence rates p for three errors ε, εσ and εd in function of the loading intensity

5 EFFECTIVITY INDEX

Once defined, the effectiveness properties of this indicator should be investigated. For this purpose, the effectivity index noted here by γ seems to be the most used tool. It is defined generally as follows:

γ = e s e (42)

with es the approximated error provided by the present indicator such as established in (24) and e the exact error computed with respect to the exact solution. However, since this later is not available, e is obtained by using a reference solution corresponding to a very refined mesh procedure.

Moreover, the calculation of the absolute error es requires the knowledge of the constant Cdintroduced in (21). However, this parameter remains unknown so that one can only compute the relative errors: å, å ó and å d (where this constant Cd is simplified). Thus, this does not allows computing directly this approximated error es. 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:

γ σ = σ ˜ h σ h σ , Ω σ ˜ h r e f σ h r e f σ , Ω (43)

where σ˜href and σhref are respectively the smoothed and non-smoothed reference solutions cited above and σ˜h and σh 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:

γ d = ( Ω ( d ˜ h d h ) 2 d Ω ) 1 2 ( Ω ( d ˜ h r e f d h r e f ) 2 d Ω ) 1 2 (44)

where d˜href and dhref are respectively the smoothed and non-smoothed damage reference solutions cited above and d˜h and dh 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.

Table 2
Effectivity indexes for different loading levels

6 CHARACTERISTIC LENGTH INFLUENCES

In the framework of the nonlocal approach, the characteristic length lc 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 lc 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 QUA4 elements of size h equal to 2.5 mm and with the same properties as in section 4 above but with different values of lc 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 lc (10 and 20 mm), the damaged arc zones are completely formed. Nevertheless, a larger bandwidth for the case lc 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 lc 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 lc 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 lc, 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 Y¯nl.This latter depends directly on the characteristic length parameter lc . It should be also noted that lcis generally determined by identification with experimental results.

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

Figure 16
Damage error distributions of SENB tests obtained under the same intensity of loading (138.F) and simulated with different values of lc (=10, 20, 30 and 40 mm)

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

Table 4
loading intensity values corresponding to the first appearance of damage for different values of lc.

For the second aspect considering the influence of lc on the error convergence rate, five cases corresponding to five different values of lc (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 lc, 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.

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

However, despite all these curves present either strong variation either weak values, the first two curves (lc equals 10 and 20 mm) are those with relatively high p-values. Hence, for the direct mesh adaptation performed manually for each increment loading (i.e. without automatic mesh adaptation for all the steps) the difficulty can be reduced by the choice of low values of lc. However, it is not often possible to make this choice because this parameter lc is obtained practically by identification with experimental results. Moreover, more investigations with other tests and other damage models should be performed with future researches.

6 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 h-generation 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. 1986Ladeveze, P., Coffignal, G., Pelle, J.P., (1986). Accuracy of elastoplastic and dynamic analysis, in Accuracy estimates and adaptative refinements in Finite Element computations. Chapter 11, 181-203, Babuska, Gago, Oliveira, Zienkiewicz Editors, J. Wiley.) for a relative global error å if:

{ f o r ε = ε 0 , t h e a c c u r a c y p r e s c r i b e d b y t h e u s e r N * i s m i n i m a l ( N * i s t h e t o t a l n u m b e r o f e l e m e n t s o f t h e m e s h T * ) (45)

In order to solve the optimization problem of (45), we adopt the following iterative procedure:

  • - Computation on a coarse mesh T

  • - Computation of both the global relative errors å,εσ , εdand the global, stress and damage error cards

  • - Determination of the optimized mesh T*provided by the calculations of a prescribed size card obtained by using the method developed by Coorevits et al. (1995Coorevits, P., Ladeveze, P., Pelle, J.P., (1995). An automatic procedure for finite element analysis in 2D elasticity, Computer Methods in Applied Mechanics and Engineering, 121: 91-120.)

  • - and a second computation of the relative errors on the new mesh T* is carried out.

If the obtained relative global error å is greater than ε0, this procedure is repeated as far as a precision close to ε0 is reached. It is worth noting that this prescribed size card is determined by the computation of size modification coefficients rE on each element E of the mesh T:

r E = h E * h E (46)

where hE denotes the size of E and hE* represents the size that must be imposed to the elements of T* in the region of E in order to ensure optimality. The computation of these coefficients rE uses the rate of convergence of the error p explained in section 4. Therefore, to compute these coefficients rE, a technique detailed in Coorevits et al. (1994Coorevits, P., Ladeveze, P., Pelle, J.P., (1994). Mesh optimization for problems with steep gradient areas. Engineering computations. 11: 129-144.) that automatically takes into account the steep gradient regions is used. The mesh T* 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 (in 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.

5.1 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 lc 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.

Figure 18
initial mesh with 765 elements and 789 nodes

Figure 19
deformed shape and damage distribution (amplified 63 times)

Figure 20
initial mesh stress error card with åó =11.37%

Figure 21
initial mesh damage error card with åd=4.43%

Figure 22
element contribution card to a relative global error of å= 12.20%

Figure 23
prescribed element size card for a prescribed error of å0 =5%.

Figure 24
optimized mesh of 5958 qua4 elements with a final accuracy of å = 4.86%

Figure 25
global error card of the optimized mesh (with å=4.86%)

Figure 26
optimized mesh stress error card with åó = 3.39%

Figure 27
optimized mesh damage error card with åd = 3.48%

5.2 Example 02: The DENB test (Double edge notched beam)

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 lc 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 å = 4.95%, 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.

Figure 28
DENB test (Geers et al. 2000Geers MGD, De Borst R, Peerlings RHJ. (2000) Damage and crack modeling in single-edge and double-edge notched concrete beams. Engineering Fracture Mechanics; 65:247-261.)

Figure 29
Initial mesh with 842 QUA4 elements and 892 nodes.

Figure 30
Non smoothed (on the left) and smoothed (on the right) damage field distributions

Figure 31
Initial mesh stress error card with åó= 9.50%

Figure 32
Initial mesh damage error card with åd=5.65%

Figure 33
Initial mesh global error card with a relative global error å= 11.05%

Figure 34
Prescribed element size card for a prescribed error of å0 =5%

Figure 35
optimized mesh of 6110 qua4elements with a final accuracy of å = 4.95%

Figure 36
global error card of the optimized mesh (with å=4.95%)

Figure 37
optimized mesh stress error card with åó = 2.72%

Figure 38
optimized mesh damage error card with åd = 4.14%

6 Discussion and concluding remarks

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.

References

  • Babuska, I., Rheinboldt, W.C., (1978). Error estimates for adaptive finite element computations. SIAM J. Numer. Anal. vol. 15, n°4: 736-754.
  • Bazant Z.-P., (1986).”Mechanics of distributed cracking” Appl Mech Rev, ASME, 39: 675-705.
  • Carpinteri, A., Valente, S., Ferrara, G., Melchiorri, G., (1993). Is mode II fracture energy a real material property. Computers and Structures 48: 397-413.
  • Cast3m. (2018). www-Cast3m.cea.fr
    » www-Cast3m.cea.fr
  • Comi, C., Perego, U., (2001) Numerical aspects of nonlocal damage analyses. Revue européenne des éléments finis, 10:227-242.
  • Comi, C., Perego, U., (2002) Finite element strategies for damage assessment up to failure. In Proceedings of the 6th National Congress SIMAI, Chia Laguna, Italy.
  • Coorevits, P., Ladeveze, P., Pelle, J.P., (1994). Mesh optimization for problems with steep gradient areas. Engineering computations. 11: 129-144.
  • Coorevits, P., Ladeveze, P., Pelle, J.P., (1995). An automatic procedure for finite element analysis in 2D elasticity, Computer Methods in Applied Mechanics and Engineering, 121: 91-120.
  • Delmas, J. (2011) https://www.code-aster.org/V2/doc/v10/fr/man_r/r3/r3.06.03.pdf
    » https://www.code-aster.org/V2/doc/v10/fr/man_r/r3/r3.06.03.pdf
  • De Vree, J.H.P., Brekelmans, W.A.M., Gils, M.A.J., (1995). Comparison of nonlocal approaches in continuum damage mechanics. Computers and Structures 55: 581-588.
  • Geers MGD, De Borst R, Peerlings RHJ. (2000) Damage and crack modeling in single-edge and double-edge notched concrete beams. Engineering Fracture Mechanics; 65:247-261.
  • Hinton, E., Campbell, J.S. (1974). Local and global smoothing of discontinuous finite element functions using a least squares method. Int. J. Numer. Meth. Engrg. 8: 461-480.
  • Hinton E., Scott F.C., Ricketts R.E. (1975) - Local least squares stress smoothing for parabolic isoparametric elements - Int. J. for Num. Meth. in Eng. 9: 235 - 256.
  • Gerasimov, T; Stein E. and Wriggers P. (2015) Constant-free explicit error estimator with sharp upper error bound property for adaptive FE analysis in elasticity and fracture Int. J. Numer. Meth. Engng; 101:79-126.
  • Giry C., Dufour F., Mazars J. (2011) Stress based nonlocal damage model, International Journal of Solids and Structures, 48: 3431-3443.
  • Huerta, A., Rodríguez-Ferran A., Díez P. (2002) Error estimation and adaptivity for nonlinear FE analysis. Int. J. Appl. Math. Comput. Sci.12:59-70.
  • Jirásek, M. and Grassl, P. (2008) Evaluation of directional mesh bias in concrete fracture simulations using continuum damage models. Engineering Fracture Mechanics, 75:1921-1943.
  • Kachanov, L. M. (1958). On the time to failure under creep conditions. Izv. AN SSSR, Otd. Tekhn. Nauk 8: 26-31.
  • Ladeveze, P., (1975). Comparaison de modèles de milieux continus, Ph.D. thesis, Paris VI, university.
  • Ladeveze, P., Leguillon, D. (1983). Error estimate procedure in the finite element method and applications. SIAM J. Num. Anal. Vol. 20 N° 3: 483-509.
  • Ladeveze, P., Coffignal, G., Pelle, J.P., (1986). Accuracy of elastoplastic and dynamic analysis, in Accuracy estimates and adaptative refinements in Finite Element computations. Chapter 11, 181-203, Babuska, Gago, Oliveira, Zienkiewicz Editors, J. Wiley.
  • Mazars, J. (1984). Application de la mécanique de l’endommagement au comportement non linéaire et à la rupture du béton de structure, Thèse de doctorat d’état de l’Université Paris VI.
  • Mazars J., Hamond F., Grange S., (2014) A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings, Materials and Structures.
  • Rabotnov, YN. (1969) Creep problem in structural members. North Holland, Amsterdam.
  • Needleman, A. (1988) Material rate dependence and mesh sensitivity in localization problems Comput. Meth. Appl. Mech. Engrg. Vol. 67 p.69-85.
  • Pijaudier-Cabot G, Bazant Z -P. (1987) “Nonlocal damage theory” Journal of Engineering Mechanics, ASCE, vol 113, n°10:1512-1533.
  • Ramtani, S. (1990). Contribution à la modélisation du comportement multi-axial du béton endommagé avec description du caractère unilatéral. Ph.D. thesis, Paris VI, university.
  • Rodríguez-Ferran, A., Huerta, A. (2000) Error estimation and adaptivity for nonlocal damage models. Int. J. Solids Struct. 37:7501-7528.
  • Rodríguez-Ferran A., Arbós I. and Huerta A. (2001): Adaptive analysis based on error estimation for nonlocal damage models.-Revue européenne des éléments finis (Special issue: Numerical Modelling in Damage Mechanics), Vol. 10, No. 2-4: 193-207.
  • Richard, B, Ragueneau, F., Cremona, C., Adelaide L. (2010) Isotropic continuum damage mechanics for concrete under cyclic loading: Stiffness recovery, inelastic strains and frictional sliding. Engineering Fracture Mechanics, Elsevier, 77 (8):1203-1223.
  • Schlangen, (1993). Experimental and numerical analysis of fracture processes concrete. Ph.D. thesis, Delft University of Technology.
  • Terrien, M. (1980). Emission acoustique et comportement mécanique post-critique d’un béton sollicité en traction. BLPC, vol 105, p 65-72.
  • Verfürth R. (2013) A Posteriori Error Estimation Techniques for Finite Element Methods. Oxford University Press: Oxford.
  • Zienkiewicz, O.C., Zhu, J.Z., (1987). A simple error estimator and adaptive procedure for practical engineering analysis. Int. J. for Num. Meth. in Engng. 24, 337-357.
  • Zienkiewicz, O. C., Zhu, J. Z., (1992). The Superconvergent Patch Recovery and adaptive finite element refinement. Computer Methods in Applied Mechanics and Engineering, 101, Iss 1-3: 207-224.
  • Zienkiewicz, O.C. (2006). The background of error estimation and adaptivity in finite element computations. Comput. Methods Appl. Mech. Engrg. 195: 207-213.
  • Available Online: June 26, 2018.

Publication Dates

  • Publication in this collection
    2018

History

  • Received
    02 Jan 2018
  • Reviewed
    29 May 2018
  • Accepted
    22 June 2018
Individual owner www.lajss.org - São Paulo - SP - Brazil
E-mail: lajsssecretary@gmsie.usp.br