Local adaptive refinement method applied to solid mechanics

A good spatial discretization is of prime interest in the accuracy of the Finite Element Method. This paper presents a new refinement criterion dedicated to an h-type refinement method called Conforming Hierarchical Adaptive Refinement MethodS (CHARMS) and applied to solid mechanics. This method produces conformally refined meshes and deals with refinement from a basis function point of view. The proposed refinement criterion allow adaptive refinement where the mesh is still too coarse and where a strain or a stress field has a large value or a large gradient. The sensitivity of the criterion to the value or to the gradient ca be adjusted. The method and the criteria are validated through 2-D test cases. One limitation of the h-adaptive refinement method is highlighted: the discretization of boundary curves.


Introduction
Higher accuracy in Finite Element Method (FEM) can be obtained by adjusting spatial resolution of domain discretization. A solution is to perform a uniform mesh refinement until finest singularities are captured. This approach can lead to a cumbersome calculation if the studied structure exhibits large dimensions compared to the size of some local heterogeneity that drives the overall toughness, e.g. a crack across the thickness of a thin rod or a plate. An alternative is to perform local adaptive refinement. For applications we have in mind, we are looking to divide by about 10 the CPU time with a local adaptive refinement compared to a uniform mesh refinement at same finite element accuracy. Four main refinement methods can be listed: • the r-adaptive method consists in moving nodes of an initial mesh in the region of interest [12], • the s-adaptive method, where the initial mesh is overlaid with finer local meshes [11], • the p-adaptive method consists in increasing basis function degree where it is needed [4,24], • the h-adaptive method, where elements are split in some regions to decrease their size and locally increase the mesh density [4].
The major disadvantage of the p-method is its difficult use in case of geometrical singularities. Splitting elements with a h-type method generally leads to non-conformities. These methods have been also combined to increase their efficiency: hp-method [2,4] or hr-method [14,15] among others. A convenient mesh refinement method has both to keep good element shape quality in order to avoid inaccurate numerical predictions due to elongated/flat/degenerated elements and to lead to conformal refined meshes as required by FEM. Several methods have been developed to overcome this issue of non-conforming hand s-methods. For example, to restore continuity of finite element space, hanging (constrained) nodes, which are equipped with the Degrees Of Freedom (DOFs) expressed with DOFs of the active nodes, are added in suitable region. Another way is to use Lagrangian multipliers or penalty methods. An alternative is also to split extra elements in order to conveniently connect refined zone to the coarse one. A method based on this idea has been developed and extended to hexahedron in [19] to keep the number of degrees of freedom constant (r-method).
In this paper, a h-type refinement method called Conforming Hierarchical Adaptive Refinement MethodS (CHARMS) is studied. This method, initially developed in the framework of fluid mechanics [6,16] and more recently applied to solid mechanics [10,13], allows both (i) to preserve element shape quality while producing a conforming mesh and (ii) to perform local adaptive refinement and unrefinement. This method deals with refinement from a basis function point of view and not from a geometrical point of view. The general idea is to construct a hierarchy of nested meshes, where at a given level of refinement each basis functions can be expressed as a linear combination of basis functions associated to the immediate finer or coarser mesh. This refinement method is here used with new refinement criterion and is validated through 2-D test cases. Notice that CHARMS method is also valid in three dimensions and preserves conformal finite elements at all stages of refinement even with non-conformal geometric nodes present.
In Section 2, notations, local adaptive refinement procedure and refinement pattern are introduced. Although unrefinement can be easily performed, only mesh refinement is explained and validated in this study. A new refinement criterion is proposed in Section 3. Contrary to usual refinement criteria based on error estimations, this criterion is based on a coupling between intensity of a mechanical field, its gradient and local characteristic size of the mesh. Validation of the method applied to FEM in solid mechanics is exposed through 2-D test cases and compared to fine uniform meshes. In the last section, a limitation of the method is pointed out: the discretization of boundary curves.

Conforming, Hierarchical, Adaptive Refinement MethodS
CHARMS is here applied to conforming Lagrangian elements and based on the refinement/unrefinement of basis functions. This approach relies on the nested relation between approximation spaces X 0 ⊂ . . . ⊂ X J , J ≥ 1, generated by a set of basis functions where the coarsest space is X 0 and the finer spaces have increasing indexes up to X J . In this setting, adding/removing basis functions in B j allows both to refine/unrefine the approximation space X j and to ensure a linear independency between the set of basis functions in B j [6,13,16,17].
Given that X j ⊂ X j+1 , all basis functions in B j can be expressed as linear combination of basis functions in B j+1 . These linear combinations define children-parents relationships between two consecutive refinement levels. To unrefine a child basis function (level j + 1), this function is replaced by a linear combination of parents basis functions (level j).
In order to perform refinement using CHARMS, a refinement pattern satisfying compatibility requirements has to be defined to generate conformal meshes from a finite element point of view. This pattern must be constituted of the same type of elements than the coarse one. The edges of element have thus to be split in several equal edges, here in two equal edges: • P 1 -triangular elements are split in 4 identical triangles, Fig. 12, • Q 1 -quadrangle elements are split in 4 identical quadrangles, Fig. 13, • Q 1 -hexahedron elements are split in 8 identical hexahedra, Fig. 14, • P 1 -tetrahedron elements are based on the strategy given in [10] which ensure that the shape quality of refined tetrahedra does not degenerate.
The key element of the CHARMS method is the refinement equation that links the coarse basis functions to the refined basis functions. This equation depends on the refinement pattern. It can be generalized and written as follows for a basis function ϕ i related to a reference element K: where superscripts (0) and (1) denote the reference element and the refinement pattern respectively, N is the number of nodes, ϕ • are the basis functions and a (l) k is the node linked to the basis function ϕ (l) k .

Refinement criteria
In order to perform local adaptive refinement, refinement criteria have to be implemented in order to select basis functions that need to be refined. Classically refinement criteria used in the literature are based on error estimation [3,5,18,26,27]. In the case of error-based adaptivity, two approaches are possible: (i) one uses two meshes [5,9,17] for error estimation; (ii) error is obtained via solution of a global problem [17,18] or some local problems in elements [1,24]. Additionally, the error-controlled adaptivity can be driven either by the interpolation error [7][8][9] or by the approximation, modeling and total errors [1,21,24,25]. The estimators can either be recovery [16,18] or residual ones [1,21,24,25]. The former estimators were developed by Zienkiewicz [26]. In this case, the basic idea is to build a new stress field σ * from the stress field σ h obtained with the FEM, in interpolating σ * either by the same basis function used in the FEM [26] or by a polynomial in the so-called 'super-convergent points' in element patches of the domain [27]. Then, an error estimator |σ * q − σ h q | is calculated for each element q. Whereas the second construction of σ * is more time consuming, it shows better results [27] and does not underestimate the error. The interpolation error control is usually associated with the iterative adaptation, while the approximation error control may lead to adaptation with fixed number of steps (three-or four-step for hp-approximation, for example) [20,21,25]. The latter approach needs availability of the convergence theory in order to relate the discretization parameters to the error level. The estimators can either be recovery [17,18] or residual ones [1,21,24,25]. There also exists another approach which is called the goal-oriented adaptivity [22], where the error estimation for the quantity of interest [1,23] can be performed on the chosen part of the mesh.
All the above-mentioned criteria need additional calculations. In order to avoid these additional numerical costs, a new criterion, which belong to the iterative adaptation approach, have been developed. It is based both on the intensity of an ad hoc field F and on its gradient taking also into account the refinement level. The general form of the criterion reads where C crit ∈ [0, 1] is a threshold (C crit = 0 means that all basis functions are refined), F i and ∇F i represent the value of a field F and of its gradient respectively for a given node i, F max and ∇F max represent the largest value of F and of its gradient respectively over the structure. N dof is the number of degrees of freedom, || · || 2 represents the L 2 -norm and the power r stands for the refinement level. For applications we have in mind, the term |F i /F max | allows to refine areas submitted to high stress or strain where crack can be initiated, while the term ||∇F i || 2 /||∇F || max 2 allows to refine areas where stress or strain fields evolve rapidly in space, i.e. around a crack tip where already existing crack can propagate. The term (1/n) r allows to take into account the mesh size by favoring refinement on non refined elements where n is the subdivision number of an edge, namely 2 in our study. The power α ∈ [0, 1] allows to weight the influence of the field intensity or its gradient. The studied mechanical problem induces the choice of the field F and the power α. In this work, we focus on the von Mises stress field F = σ V M and the von Mises strain field F = ε V M : where σ i and ε i with i ∈ {I, II, III} are, respectively, principal stress and principal strains.

Role of power α
Although in the parametric studies presented, the mechanical behavior considered is independent of time, a fictitious time step is introduced in all numerical studies. This fictitious time step allows an arbitrary evolution of the value of the threshold C crit while the applied load remains constant. By this way, successive refinements can be observed and convergence analysis can be performed. The first test case deals with a square hole in a square plate, Fig. 2. The material properties and numerical parameters are shown in Table 1 Fig. 3 shows the refinements obtained for different values of α. These refinements are observed in an area of interest depicted by a blue rectangle in Fig. 2 and located in the vicinity of the corner of the square hole. As expected, when α is close to one, the refinement criterion (3) gives more weight to the intensity of the von Mises stress field than to its gradient and the mesh is refined where large stress values are observed, i.e. within the ligament of the structural problem with respect to the load direction. In Fig. 3(f) mesh is refined everywhere but above the square hole where the stresses are screened. On the contrary, when α is close to zero the refinement criterion only concerns the gradient of the von Mises stress field and the refinement takes place where some stress concentration occurs, here at the corner of the square hole, Fig. 3(a). Extrapolating this situation to crack problems, the role of the power α could be summarized as follows: when α is close to zero, the criterion (3) induces refinement where new cracks can occur, when α is close to one, the mesh is refined at the crack tip of existing cracks.

CHARMS refinement vs homogeneous refinement
At a given number of elements, a mesh refined in convenient zones can give more precise results than a homogeneous mesh. To quantify this enhancement induced by the refinement criterion (3), a normalized error on the overall strain energy is calculated at each node of a mesh where E refi and E ref are the global strain energy on a refined mesh and on a reference mesh respectively. The studied case is the square plate with a square hole above detailed. The reference mesh is a homogeneous fine mesh of the structure with about 120 000 quadrangles whose mesh size h ref is given in Table 2.
As previously mentioned, the adaptive mesh refinement is obtained by deceasing the value of the threshold C crit (by step of 0.02) whilst the load remains constant. Fig. 4 shows the evolution of this energy error (5) with respect to the number of elements for different values of α. The thick black curve shows the energy errors for homogeneous refinement while the colored curves correspond to CHARMS refinement. Both the CHARMS   Fig. 4. CHARMS refinements (top row), homogeneous refinements (bottom row). Stages 1 and 5 correspond to the same state refinement and the homogeneous refinement start from the same coarse mesh with about 7 000 elements (stage 1 in Fig. 4 and 5) and finish at the same fine mesh with about 120 000 elements (stage 5 in Fig. 4 and 5). At these stages 1 and 5, the energy error e is thus the same for the CHARMS approach and the homogeneous one, with a larger error at stage 1 since it corresponds to the coarser mesh. During stages from 2 to 4, the meshes for the CHARMS refinement do not correspond to the meshes of the homogeneous refinement: the refinement takes place where the criterion (3) is reached. At the first, second and third refinements (stages 2, 3 and 4) the enhancement proposed by the CHARMS approach is significant and the energy error is about two times smaller than for the homogeneous refinement. However, no theoretical estimate of this energy error has been carried out and the estimates in Fig. 4 is pragmatic. In this example, only two CHARMS refinement steps are allowed. Then, when the mesh size is 4 times smaller than in the initial coarse mesh, the refinement process takes place somewhere else and the interest of the CHARMS refinement disappears.
Moreover, we can notice in Fig. 4 that the energy error decreases with α. In this example, the high local stress concentration due to the square shape of the hole induces a predominant role of the gradient field in the criteria (α < 0.5).
Two main conclusions can be drawn on this example: (i) since successive refinements are allowed, the proposed method can significantly reduce the finite element error (here in an energy sense) compared to homogeneous refinements; (ii) the power α has to be adapted to each problem and in practical calculations, when the parametric studies on α are not available, the value of α can be set to the mid value 0.5.

Role of the field F
This test case deals with bimaterial square plate (Fig. 6) constituted with stiff inclusions in softer matrix. Material properties and numerical parameters are shown in Table 3 and 4. A constant vertical displacement is applied on the top of the specimen while the bottom remains fixed along the vertical direction. In this study again two refinement steps are applied, the power α is fixed to 0.4 and two fields of interest F are tested namely the von Mises stress Fig. 7(a) and the von Mises strain Fig. 7(b). As previously, the normalized error e (5) is calculated at each node of the refined mesh to evaluate the accuracy of the method. Again, the value of the threshold C crit is also reduced by 0.02 at each fictitious time step.   Fig. 7 shows that the obtained refined mesh clearly depends on the choice of the field of interest used in the criterion (3). When the field of interest is set to the von Mises strain, the refinement is located around the inclusion where maximum equivalent strain is reached. When von Mises strain is used, the refinement is located inside the inclusions where maximum   Fig. 8 shows the evolutions of energy normalized errors for the two criteria. For the two cases, the error rapidly decreases with respect to the number of elements and since the calculations are performed in linear elasticity the field of interest (stress or strain) has no influence on the FEM quality of the refined meshes.

Limitation: discretization of curved edge
In this section, a limitation of the CHARMS method is pointed out. Indeed, this method is based on the refinement of the basis functions (h-method) and not on remeshing. Geometric approximation of boundary curves are thus not improved during the refinement steps. We will highlight this limitation on the classical test case of an infinitely long cylinder loaded with an internal imposed displacement, with internal radius a and external radius b. The closed-form solution of this problem reads where the imposed internal displacement is u r (a) = u 0 = 0 and the external displacement u r (b) = 0. For this test case, a and b have been chosen so that ξ = 2.
The field F is the von Mises stress and the power α = 1 because there is no stress concentration in this example. The refinement pattern used is the P 1 -triangle element (Appendix A) and the maximum refinement level is 2. For the simulation, we consider a quarter of disc in plane strain assumption (Fig. 9). The material and numerical properties are presented in Table 5 and 6. The value of the threshold C crit is decreased by 0.1 at each time step so that at final time all basis functions are refined.
Knowing the closed-form solution, an exact relative error is calculated at the red dot exhibited in Fig. 9:ẽ • a b Fig. 9. Geometry of the internal displacement imposed cylinder: a = 100 mm is the internal radius and b = 200 mm is the external radius. The red point is the point of interest where u sim r is the displacement calculated at the red dot with adaptive mesh refinement and u exact r is the analytical displacement obtained with (6). In this particular case of an initial coarse mesh of a curved boundary, the curves in Fig. 10 exhibits that the error is more important for a refined mesh using CHARMS method than for a mesh geometrically refined. Effectively, unlike what is classically exhibited using geometric Fig. 10. Hollow cylinder: evolution of the exact relative error (7) calculated at the red node obtained with mesh refinement (blue curve) and with homogeneous meshes (black thick curve) with respect to the number of elements Fig. 11. Mesh obtained with the spatial discretization of a curved geometry. The real curvature is in black dots, the discretization obtained with CHARMS is in blue dots (the coarse mesh is the black square) and the discretization using a homogeneous mesh is in green dots. e d represents the spatial discretization error between the homogeneous mesh and the refined mesh mesh refinement and previous results ( Fig. 4 and 8), we do not observe a monotonous decrease in error during successive refinement steps using the CHARMS method. Fig. 11 illustrates the fact that geometric remeshing naturally induces a better edge interpolation. In fact using the CHARMS method coarser the initial mesh is, higher the geometrical error is while refining. The refinement method CHARMS cannot increase geometrical interpolation of an initial coarse mesh. Just as it is done in fluid mechanics for meshes to capture boundary layer behavior, it is important in solid mechanics for structures with large radii of curvature to use an initial mesh with good geometry interpolation before making refinement steps.

Conclusions
A local adaptive method using a basis function point of view is presented and validated in solid mechanics with elastic behavior. This method enables to refine elements without degenerating their shape quality and to implicitly handle the non-conformities. During a step of refinement the current basis functions are replaced by a linear combination of these functions based on a refinement pattern. A refinement criterion based on the value of an ad hoc field and its gradient has been developed. Compared to classical criteria, this approach is not time-consuming since it does not require any additional calculation but interpolated errors cannot be quantified a priori. The presented test cases show good results and put in evidence that the choice of the field of interest influences the region of refinement. At the end, one limitation of the h-adaptive refinement method is pointed out. If the specimen has a severe curvature, the initial mesh has to be fine enough because the CHARMS method, based on the interpolation of the shape functions, does not allow an improvement of the geometrical interpolation.
In future work, CHARMS will be applied in case of crack propagation where refinement and unrefinement can be performed adaptively. The refinement pattern is detailed in Fig. 13 and the corresponding equations are given in (2).   Fig. 14. Q 1 -cuboid refinement pattern, the refined basis functions coordinates and refinement equations