Stability of incompressible formulations enriched with X-FEM

Stability of incompressible formulations enriched with X-FEM G. Legrain aN. Mo¨es a,∗A. Huerta b aGeM Institute - ´ Ecole Centrale de Nantes / Universit´e de Nantes / CNRS 1 Rue de la No¨e, 44321 Nantes, France. bLaboratori de C`alcul Num`eric (LaC`aN) - Edifici C2, Campus Nord, Universitat Polit`ecnica de Catalunya - E-08034 Barcelona, Spain. Abstract The treatment of (near-)incompressibility is a major concern for applications involving rubber-like materials, or when important plastic flows occurs as in forming processes. The use of mixed finite element methods is known to prevent the locking of the finite element approximation in the incompressible limit. However, it also introduces a critical condition for the stability of the formulation, called the infsup or LBB condition. Recently, the finite element method has evolved with the introduction of the partition of unity. The eXtended Finite Element Method (XFEM) uses the partition of unity to remove the need to mesh physical surfaces or to remesh them as they evolve. The enrichment of the displacement field makes it possible to treat surfaces of discontinuity inside finite elements. In this paper, some strategies are proposed for the enrichment of mixed finite element approximations in the incompressible setting. The case of holes, material interfaces and cracks are considered. Numerical examples show that for well chosen enrichment strategies, the finite element convergence rate is preserved and the inf-sup condition is passed. Key words: Mixed formulation, X-FEM, Partition of Unity, inf-sup condition, incompressibility, holes, inclusions, fracture mechanics ∗Corresponding author Email addresses: [email protected] (G. Legrain), [email protected] (N. Mo¨es), [email protected] (A. Huerta). Preprint submitted to Elsevier Science 13 September 2005 Legrain, N. Moes and A. Huerta,Stability of incompressible formulations enriched with X-FEM, G., Computer Methods in Applied Mechanics and Engineering, Vol. 197, Issues 21-24, pp. 1835-1849, 2008 1 Introduction Displacement based finite element methods are nowadays abundantly used in engineer analysis. Indeed, they can solve a wide variety of problems, and have now been deeply mathematically investigated. However there still exists two main drawbacks for these methods. First, the treatment of incompressible or nearly incompressible problems necessitates the use of adapted formulations. If not, incompressibility constraint locks the approximation, leading for instance to non-physical displacement fields. Second, the generation and especially the update of the mesh in complex 3D settings for evolving boundaries such as cracks, material interfaces and voids still lacks robustness, and involves important human effort. Several techniques have been developed to respond to the locking issue. For instance, the selective-reduced-integration procedures [1–3] or the Bbar approach of Hugues [4] in which the volumetric part of the strain tensor is evaluated at the center of the element. Another way to avoid locking is to enhance the strain tensor in order to enlarge the space on which the minimization is performed, and meet the divergence-free condition (enhanced assumed strain methods, see [5–9]). Here, we will focus on two-field mixed finite element methods. The incompressibility constraint is weakened by the introduction of the pressure field. This alleviates locking at the price of additional pressure unknowns. However, mixed finite element methods are not stable in all cases, some of them showing spurious pressure oscillations if displacement and pressure spaces are not chosen carefully. To be stable, a mixed formulation must verify consistency, ellipticity and the so called inf-sup (or LBB) condition. The later is a severe condition which depends on the connection between the displacement and pressure approximation spaces. Stable mixed formulations can be obtained by stabilizing non-stable formulations with the use of parameters whose values may depend on the problem at hand. Otherwise, one has to work with approximation spaces which passe the inf-sup condition. To prove that a displacement-pressure pair satisfies the inf-sup condition is not a trivial task. However, a numerical test has been proposed by Brezzi and Fortin [10] then by Chapelle et al. [11], in order to draw a prediction on the fulfilment of the inf-sup. This test proved to be useful for the study of the stability of various mixed elements [12]. The second drawback of classical finite element methods (evolving boundaries) has been overcome by the development of alternative methods such as meshless methods in which the connectivity between the nodes is no longer obtained by the mesh, but by domains of influence which can be split by the boundaries. Moreover, the approximation basis can be enriched with functions coming from the physical knowledge of the problem. Note that in the incompressible limit, meshless methods are now known to lock [13,14] as classical finite elements. Thus, some strategies have been developed to circumvent this 2 issue [13,15]. Apart from mixing meshless methods and finite elements [16] another alternative to overcome the re-meshing issue in finite elements is to use the eXtended Finite Element Method (X-FEM) based on the partition of unity framework introduced by Babuˇska and Melenk [17]. Proper enrichment of the finite element basis makes it possible to model crack, material inclusions and holes with non-conforming meshes. The X-FEM method has been used for the simulation of a wide variety of problems such as fracture mechanics problems (2D[18–20], 3D[21–23], plates [24,25], cohesive zone modeling [26,27], dynamic fracture [28], nonlinear fracture mechanics[29–31]), holes [32,33], but also material inclusions [33,34] or multiple phase flows [35]. Here, we focus on the application of this method to mixed formulations for the treatment of holes, material inclusions and cracks in the incompressible limit. Bbar or selective-reduced formulations are not considered, because they do not seem to be generalized easily to enriched displacement fields. The main contribution of this paper is the design of enrichment strategies for the pressure and displacement fields, so that it leads to a stable formulation. The enrichment of mixed finite element approximations has already been used by Dolbow et al. [25] and Areias et al.[31] for fracture mechanics in plates and shells, and by Wagner et al. [36] for rigid particles in Stokes flow. However, the stability and the convergence of these approaches was not studied. The latest work concerning volumetric incompressibility was proposed by Dolbow and Devan [29]. In this papers, the authors focus on the application of the enhanced assumed strain method to X-FEM in large strain. This approach seems to lead to a stable low order formulation in the case of nearly incompressible nonlinear fracture mechanics. However, the stability of the method was not shown, and the influence of the near-tip enrichment was not studied. More precisely, it is not clear whether the near-tip enrichment could make this approach unstable, as the construction of an orthogonal enhanced strain field becomes difficult with non polynomial functions. The paper is organized as follows: first, the governing equations of incompressible linear elasticity are recalled. The conditions for the stability of mixed formulations are also reviewed. Next, some strategies are proposed to keep the stability of enriched finite elements. The case of holes, material interfaces and 2D cracks are presented. Finally, in a last section the stability of these strategies is investigated. 2 Governing Equations In this section, we focus on the design of stable mixed formulations for the treatment of incompressible elasticity. First, the equations governing incompressible linear elasticity are recalled. Then the inf-sup condition is presented together with a numerical test. 3 and εDis the deviatoric strain operator: εD=ε−εV 3I(5) When the material tends to incompressibility, the bulk modulus tends to infinity. This means that εVmust tend to zero in order to meet condition (2) (the displacement field must be divergence free in the incompressible limit). εV= div(u)−→ 0 as ν→0.5 (6) The strong form (1) is equivalent to the stationarity of a displacement potential Π: Π(u) = 1 2ZΩ ε(u) : C:ε(u)dΩ−ZΩ u·bdΩ−Z∂Ωt u·TddΓ (7) In order to model incompressible or almost incompressible problems, a two field principle is considered by introducing a second variable (the hydrostatic pressure p) in the potential (7): p=−κεV(u) = −1 3Tr(σ) (8) When κincreases, the volumetric strain εVdecreases and becomes very small. For total incompressibility, the bulk modulus is infinite, the volumetric strain is zero, and the pressure remains finite (of the order of the applied boundary tractions). The stress tensor is then expressed as: σ=−p I + 2µ εDin Ω (9) The solution of the governing differential equations (1) now involves two variables: the displacement field and the pressure field. Writing the two fields variational principle, the total potential for the u−pformulation is expressed as: χ(u, p) = 1 2ZΩεD(u) : C:εD(u)dΩ−ZΩ u·bdΩ−Z∂Ωt u·TddΓ (10) −1 2ZΩ p2 κdΩ−ZΩ p εV(u)dΩ Invoking the stationarity of χ(u, p) with respect to the two independent variables uand p, we obtain: ZΩδεD:C:εDdΩ−ZΩp δεVdΩ = R(δv) (11) −ZΩp κ+εVδp dΩ = 0 (12) where R(δv) represents the virtual work of the external loads. In the case of 5 expresses as: The existence of a stable finite element approximate solution (uh, ph)depends on choosing a pair of spaces Vhand Qhsuch that the following condition holds: inf qh∈Qhsup vh∈VhRΩqhdiv vhdΩ kvhk1kqhk0 ⩾β > 0 (17) where k · k1and k · k0indicates H1and L2norms respectively and βis independent of the mesh size h. If the inf-sup compatibility condition is satisfied, then there exists a unique uh∈ Vhand a ph∈ Qh(determined up to an arbitrary constant in the case of purely Dirichlet boundary conditions). 2.2.2 Numerical assessment of the inf-sup condition As seen before, the prediction of the stability of a mixed formulation involves the fulfillment of the inf-sup criterion. This criterion is however impossible to prove for practical situations. This is why the numerical evaluation of the inf-sup condition has received considerable attention [11,10]. This numerical evaluation, although not equivalent to the analytical inf-sup, gives indications on whether (17) is fulfilled or not for a given set of finite element discretizations. The numerical inf-sup test is based on the following theorem. Proposition 1 Let Muu and Mpp be the mass matrices associated to the scalar products of Vhand Qhrespectively and let µmin be the smallest non zero eigenvalue defined by the following eigenproblem: KupTMuu−1Kup q=µ2Mpp q(18) then the value of βis simply µmin. The proof can be found in [37] or [10]. The numerical test proposed in [11] consists in testing a particular formulation by calculating βusing meshes of increasing refinement. On the basis of three or four results it can be predicted whether the inf-sup value is probably bounded from underneath or, on the contrary, goes down to zero when the mesh is refined. The reliability of this test is demonstrated on several examples of elements for incompressible elasticity problems in [11]. In the following section this test is used to check the behavior of proposed enrichment strategies. However, we follow [11] and use only Suu = RΩ∇u:∇udΩ instead of Muu in (18). In order to perform the numerical infsup test, a sequence of successive refined meshes is considered. The objective is to monitor the inf-sup values, β, when hdecreases. If a steady decrease in log(β) is observed when h goes to zero, the element is predicted to violate the inf-sup condition and said to fail the numerical test. But, if the log(β) value is stable as the number of elements increases, the test is numerically passed. 7 3 X-FEM Discretization 3.1 Displacement field With classical finite elements, the approximation of a vector field uon an element Ωeis written as: u(x)|Ωe= nu X α=1 uαNα u(x) (19) where nuis the number of coefficients describing the approximation of the displacement over the element, uαis the αth coefficient of this approximation and Nα uis the vectorial shape function associated to the coefficient uα. Within the partition of unity, the approximation is enriched as: u(x)|Ωe= nu X α=1 Nα u uα+ nenr X β=1 aα βφu β(x) (20) where nenr is the number of enrichment modes, aα βis the additional dof associated to dof αand φu βstands for the βth scalar enrichment function. The number and the expression of the enrichment functions vary with the problem to model. The expression of this enrichment function will be recalled for holes, inclusions and fracture mechanics in the next sections. 3.2 Pressure field Using the same scheme, the pressure approximation is written as: p(x)|Ωe= np X α=1 Nα p(pα+aαφp(x)) (21) where npis the number of coefficients describing the approximation over the element, pαis the αth coefficient of this approximation and Nα pis the scalar shape function associated to the coefficient pα,aαis the additional dof associated to dof αand φpstands for the scalar pressure enrichment function. The key issue is the combined choice of enrichment functions φuand φpsuch that the whole enriched approximation (displacement and pressure) passes the inf-sup condition. 8 function are the most effective since they pass the inf-sup test. On the contrary, Heaviside based strategies do not pass this test and should not be considered. The formulations where only the linear part of the displacement is enriched (N˚4 and 8) with the ridge are interesting, since less degrees of freedom are involved in the approximation. This should be important in the context of an extension to three dimensional studies. 4.3 Incompressible fracture mechanics The resolution of compressible fracture mechanics problems has been extensively studied in the context of the X-FEM for both 2D [39,18,40,32,24] and 3D fracture mechanics [21,22]. The most common enrichment strategy consists in using the asymptotic displacement field as an enrichment for the displacement finite element approximation. In the context of incompressible media, the analytical asymptotic displacement field (Westergaard solution) is shown to be identical to the limit of the compressible one. The asymptotic evolution of the pressure field can be obtained also using the Westergaard solution. p(r, θ) = 2KI 3√2π r cos θ 2!+2KII 3√2π r sin θ 2!(25) φu=(√rsin θ 2!,√rcos θ 2!,√rsin θ 2!sin(θ), √rcos θ 2!sin(θ))(26) We use these expressions as an enrichment for the pressure field in the near-tip region. Thus, the enrichment basis for the pressure is expressed as:      φp I=1 √rcos θ 2 φp II =1 √rsin θ 2(27) Note that a classical Heaviside enrichment is considered for both pressure and displacement for nodes whose support is fully cut by the crack, and that only the Mini element is considered hereunder. 4.3.1 Convergence study Consider a domain Ω = [−1,1] ×[−1,1] under tension (see Figure 23). The tensions applied on the boundary of the domain are related to the exact ten20 lation. Moreover, the use of the geometrical enrichment leads to an improved convergence case, similar to the compressible rate. Finally, degrees of freedom can be saved by only enriching the linear part of the approximation. 5 Conclusion Some strategies for enriching existing mixed finite element methods have been presented. These strategies are natural extensions of the displacement-based X-FEM, and are shown to preserve the classical finite element convergence rate. The stability of these strategies has been shown through the numerical inf-sup test. However, quadratic-based elements could not be tested completely, as in some cases the linear interpolation of the level-set leads to a degraded rate of convergence. The construction and the validation of isoparametric quadratic elements will be the subject of a forthcoming paper. The method should also be applied to finite strain mechanics, as the fulfilment of the inf-sup condition seems to be a prerequisite to build efficient large strain formulations [42]. 24 References [1] D.J. Naylor. Stress in nearly incompressible materials for finite elements with application to the calculation of excess pore pressure. international journal for numerical methods in engineering, 8:443–460, 1974. [2] T.J.R. Hughes, R.L. Taylor, and J.F. Levy. High Reynolds number, steady, incompressible flows by a finite element method. John Wiley & Sons, 1978. [3] S.F. Pawsey and R.W. Clough. Improved numerical integration of thick slab finite elements. International journal for numerical methods in engineering, 3:275–290, 1971. [4] TJR Hughes. Generalization of selective reduced integration procedures to anisotropic and nonlinear media. International Journal for Numerical Methods in Engineering, 15:1413–1418, 1980. [5] RL. Taylor, PJ. Beresford, and EL. Wilson. A nonconforming element for stress analysis. International Journal for Numerical Methods in Engineering, 10(6):1211–1219, 1976. [6] EL. Wilson, RL. Taylor, WP. Doherty, and J. Ghaboussi. Incompatible displacement models In Numerical and Computer Models in Structural Mechanics. 1973. [7] JC. Simo and MS. Rifai. A class of mixed assumed strain methods and the method of incompatible modes. International Journal for Numerical Methods in Engineering, 29:1595–1638, 1990. [8] JC. Simo and F. Armero. Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes. International Journal for Numerical Methods in Engineering, 33:1413–1449, 1992. [9] B.D. Reddy and J.C. Simo. stability and convergence of a class of enhanced strain methods. SIAM J. Numer. Anal., 32(6):1705–1728, 1994. [10] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. Springer, New York, 1991. [11] D. Chapelle and K.J. Bathe. The inf-sup test. Computers & structures, 47(45):537–545, 1993. [12] K.J. Bathe. The inf-sup condition and its evaluation for mixed finite element methods. Computers & Structures, 79:243–252, 2001. [13] A. Huerta and S. Fern´andez-M´endez. Locking in the incompressible limit for the element free galerkin method. international journal for numerical methods in engineering, 50, 2001. [14] J. Dolbow and T. Belytschko. Volumetric locking in the element free galerkin method. Int. J. Numer. Methods Eng., 46(6):925–942, 1999. 25