Continuous blending of SPH with finite elements

Continuous blending of SPH with finite elements ? Sonia Fern´andez-M´endez a, Javier Bonet b and Antonio Huerta a aLaboratori de C`alcul Num`eric (LaC`aN), Departament de Matem`atica Aplicada III, Universitat Polit`ecnica de Catalunya, M`odul C2, Jordi Girona 1-3, E-08034, Barcelona, Spain bCivil and Computational Engineering Center (C2EC), University of Wales Swansea, Singleton Park, Swansea, SA2 8PP, United Kingdom Abstract This paper proposes a methodology for the continuous blending of the finite element method and Smooth Particle Hydrodynamics. The coupled approximation with finite elements and particles, and the discretization of the boundary value problem with a coupled integration, are described. An integration correction is also proposed to stabilize the solution. Some numerical examples demonstrate the applicability of the method. Key words: SPH, Finite Elements, coupling, meshless, integration correction, continuous blending 1Introduction The special issues [1,2] and the present issue demonstrate the growing interest focused in meshless methods in the last decade. Meshless methods, such as the Element-Free Galerkin method (EFG) [3], the Reproducing Kernel Particle ?Sponsored by Ministerio de Ciencia y Tecnolog´ıa (grants: DPI2001-2204 and REN2001-0925-C03-01) and the Generalitat de Catalunya (grant: 2001SGR00257) Email addresses: [email protected] (Sonia Fern´andez-M´endez), [email protected] (Javier Bonet), [email protected] (Antonio Huerta). URL: www-lacan.upc.es (Antonio Huerta). Preprint submitted to Elsevier Preprint 27 October 2004 Fernández-Méndez, S.; Bonet, J. and Huerta, A., Continuous blending of SPH with finite elements, Computers and Structures, Vol. 83, Issues 17-18, pp. 1448-1458, 2005 Method (RKPM) [4] or the Smooth Particle Hydrodynamics method (SPH) [5– 9], among others, are competitive for some applications compared to classical mesh-based methods. In particular, SPH methods are widely used for fast-transient dynamic simulations, such as explosions or impact problems, because of their low computational cost and its ability to handle severe distortions [6,7]. Other meshless methods, such as EFG or RKPM, can also deal with large distortions and go beyond Finite Element (FE) computations [10], but with a higher computational cost (due to the use of Gauss quadratures or specific techniques to accurately integrate the weak form). Nevertheless, when remeshing is not an important issue (for instance in the absence of large distortions) FE computations are preferred by practitioners for two main reasons. First it should be noted that most of the users are familiar with FE methods, and mesh-free methods are still seen as research techniques, but this, of course, will change with time. Second, FE computations in the absence of remeshing can be very efficient from a computational point of view, and thus less costly than mesh-free computations. However, these advantages disappear when intensive remeshing is required. This is the case of problems with large distortions, impact, blast, etc. In these cases SPH goes beyond FE methods. Therefore, coupling FE and meshless methods seems a reasonable approach to take advantage of the best properties of both methods. FE are employed in almost everywhere and FE nodes are replaced by particles only in the regions where the limitations of FE become apparent. An important effort has been dedicated to the coupling of FE and meshless methods in the last years. Most of the literature is devoted to the coupling of FE with EFG or RKPM, see for instance [11–16]. However, the attention paid to the development of formulations linking FE with SPH methods is more modest [17,18]. Nowadays, there are several FE commercial or government-sponsored codes that also include SPH formulations, such as LS-Dyna, Pronto-3D, Europlexus and Sophia, among others. Nevertheless, in these codes the combination of SPH with FE is usually done via contact or sliding interfaces. For example, in impact problems is it usual to consider FE for the discretization of the projectile and an SPH formulation for the target, which suffers large distortions. The objective here is to discretize also the target using standard FE, and replace FE nodes by SPH particles only in regions where highly distorted elements appear during the computation. The existing contributions to link a FE region and an SPH region (with the same material) consider a transition layer of particles which are simultaneously treated as FE nodes and Lagrangian SPH particles. The influence between the different regions is taken into account with algorithmic and practical considerations. For instance, in [18] the computation of the stresses at an interface particle (node as seen from the 2 FE region) is done with contributions from the interface particles and from other standard particles, but there is no contribution from the neighboring FE nodes. Although, as it is commented in [18], this approach may be acceptable for many applications, the mathematical basis is not clear. Consequently, the standard properties (convergence, continuity, etc), one expects form FE and SPH formulations, are not guaranteed with this approach. The aim of this paper is to develop a continuous blending of SPH and FE, maintaining the reproducibility, convergence and stability properties of both methods. Several alternatives are proposed for the interpolation of a function and its derivatives, with the desired order of reproducibility (consistency) and, therefore, ensuring the convergence of the interpolation. The discretization of the problem is done with a Galerkin weak form, which is numerically integrated using Gauss quadrature in the FE region and particle integration in the rest of the domain. Section 2 recalls some basic concepts on the SPH method. Several possibilities for the approximation of the unknown and its derivatives, coupling FE and SPH, are proposed in Section 3. In particular, Section 3.1 recalls the continuous blending method, initially proposed in [13] to blend FE and Moving Least-Squares (MLS) particle approximations. The discretization of boundary value problems with a coupled integration is described in Section 4. Section 5 generalizes the integration correction, initially proposed in [19] in the SPH framework, to the coupling of SPH with FE. Finally some numerical examples are presented in order to demonstrate the applicability of the method and the proposed corrections. 2 Basic equations of SPH The basis for SPH, and other particle methods, is the approximation of an arbitrary function u(x) by the convolution b u(x) = ZΩ u(s)W(x−s, ρ)ds(1) where Wis a positive and compact supported function usually called kernel function or window function, see [15,19]. The support of Wis scaled by the dilation parameter ρ. The approximation b uwill converge to the exact function as W(x, ρ) approaches the Dirac delta function δ(x). Moreover, the kernel function Wcan be easily designed so that the approximation b uis exact for polynomials up to order m, that is b u=ufor any upolynomial of degree less or equal to m, see [20]. In that case, the approximation is said to have reproducibility of order m. Equivalently, some authors refer to this property as consistency of order m. 3 Given a set of points in the domain Ω, {xb}b∈P , usually called particles, the SPH approximation is obtained from (1) with a discrete evaluation of the integral. That is, as initially proposed by Lucy [5], the SPH approximation is b u(x) = X b∈P u(xb)Wb(x)ωb(2) where ωbis a volume associated to every particle, see [5], and Wb(x) := W(x−xb, ρ) is the weight or shape function associated to particle xb. Note that usually, even if approximation (1) reproduces polynomials up to order m, the reproducibility in the SPH approximation (2) is no longer satisfied (due to the numerical quadrature). This deficiency has important consequences for the resolution of boundary value problems in terms of accuracy, stability and convergence of the approximation. Several possibilities have been proposed in the last decade to overcome this problem. Some of the most popular particle approximations are the Corrected Smooth Particle Hydrodynamics (CSPH) [9,19], the Renormalized Meshless Derivative [8] or the mesh-free approximations based on MLS, see [15] for a general review. Remark 1 In the numerical resolution of a boundary value problem with SPH, the usual techniques are point collocation, or the discretization of a weak form with particle integration. Remark 2 In many applications, such as fast-transient dynamics, the corrections of the SPH method are formulated to verify reproducibility in the approximation of the derivatives, ∇u, [19]. The approximation of scalar quantities, such as the density, is done with a trivial 0-order correction. 3 Coupled approximation This section is devoted to the definition of coupled approximations using finite elements and particles. As shown in Figure 1, the domain Ω is discretized with a set of particles {xb}b∈P , marked with crosses, and a set of active nodes {xi}i∈N , marked with black dots. Other non-active nodes are considered to define the support of the FE shape functions (thus only associated to the geometrical interpolation). In the region where only FE are present a standard FE approximation is used, in the region where only particles have an influence a corrected particle approximation is considered. In the transition 4 the next section. 3.2 Coupled approximation of derivatives As it is commented in Remark 4, the continuous blending method provides an approximation for the derivatives with the desired reproducibility. Nevertheless, this section describes two alternative possibilities for the imposition of first-order reproducibility to the coupled approximation of ∇u. These approximations are not derived from an approximation of the function uand, therefore, they are closer to the classical SPH philosophy. In both cases the situation described in Proposition 1 is assumed. That is, (1) the FE interpolation is of order greater or equal to 1 and (2) the domain of influence of particles includes the region where FE does not have a complete basis. The first approach, similarly to [25], considers an approximation of the form ∇u(xa)≃X i∈N ∇Nh i(xa)ui+X b∈P gb(xa)ub(8) where the gradient shape functions associated to the particles are defined as gb(xa) = ωb[La∇Wb(xa) + εaδab].(9) The matrix Laand the vector εaare obtained imposing first-order reproducibility. That is, the approximation must be exact for constant functions, namely X i∈N ∇Nh i(xa) + X b∈P gb(xa) = 0(10) and for linear functions, namely X i∈N ∇Nh i(xa)⊗xi+X b∈P gb(xa)⊗xb=I.(11) Thus, the expressions for Laand εaare obtained after substitution of (9) in (10) and (11), La="I−X i∈N ∇Nh i(xa)⊗(xi−xa)# X b∈P ωb∇Wb(xa)⊗(xb−xa)  −1 (12) and εa=1 ωa −X i∈N ∇Nh i(xa)−LaX b∈P ωb∇Wb(xa) .(13) The second approach, similarly to [19], considers an approximation of the form ∇u(xa)≃X i∈N ∇Nh i(xa) [ui−ua] + X b∈P gb(xa) [ub−ua],(14) 8 where gb(xa) = Laωb∇Wb(xa).(15) Note that approximation (14) is exact for constant u, thus it is only necessary to impose the reproducibility condition for linear functions, that is X i∈N ∇Nh i(xa)⊗(xi−xa) + X b∈P gb(xa)⊗(xb−xa) = I.(16) After substitution of (15) in (16) the expression for Lais obtained, which exactly coincides with (12). However, this approach does not require the computation of vector εa. Remark 5 In the region where a complete FE base is present X i∈N ∇Nh i(xa) = 0 and X i∈N ∇Nh i(xa)⊗(xi−xa) = I. Thus, La=0,εa=0and both approximations coincide with a standard FE approximation. In the region where only particles are present, Nh i(xa) = 0 and both approximations coincide with standard SPH corrected approximations. 4 Discretization of the boundary value problem: coupled integration Two important decisions have to be taken for the numerical resolution of a boundary value problem: the choice of the functional approximation, which is commented in Section 3, and the formulation to be used for the verification of the Partial Differential Equations (PDEs) and the imposition of boundary conditions. In SPH methods it is usual to combine a weak form of the problem with particle integration, or to consider point collocation techniques for the strong imposition of the PDEs and the boundary conditions. In both cases all the computations are done at the particles, and all the quantities of interest are stored exclusively at the particles. This paper considers the use of a weak form with particle integration for several reasons. First of all, high order derivatives are avoided and, more important, the mathematical base of the weak form is well known, facilitating the analysis of the properties of the scheme. Furthermore, FE methods are based on the discretization of a weak form, thus, the most natural option for the coupling is considering the same weak form in the SPH region. To show the methodology for the discretization of the problem, and for the sake of clarity, a simple model problem is considered, the Laplace equation 9 solution. However, the solution is degraded due to the presence of too distorted elements in the punch region, see the zooms on the right. Figure 9 shows the solution obtained using the coupling of SPH and FE. The continuous blending method is used for the approximation of the displacements and its derivatives. Figure 7 shows the discretization with particles and FE, with active and non-active nodes. As usual in mesh-free methods [10], the use of an SPH approximation near the punch region allows to go beyond the FE computations. The SPH discretization avoids the problems due to the distortion of the FE. It is important to note that with an SPH discretization the refinement of the discretization near the punch region is easily performed, with no remeshing cost. Remark 10 Although it is not exploited in this example, the definition of the neighborhood between particles in the SPH region can be updated during the computations, leading to an updated Lagrangian formulation. The update of the “connectivity” between particles does not imply the usual projection errors of the FE computations (due to the projection of quantities between meshes). This is an important issue in problems with severe distortions, brittle fracture, etc. Remark 11 Obviously, the continuous blending of SPH and FE inherits the stability deficiencies of the SPH method. In problems with extremely large distortions, the SPH Lagrangian formulation will need updates of the reference configuration. Then, zero energy modes are more likely to be excited. When a reduced number of updates are necessary the error is small and can go unnoticed, but when frequent updates are performed the solution can be completely spoilt due to zero energy modes in the SPH approximation. To avoid numerical fracture in high-tension zones, standard techniques can be implemented [28–30]. Other stabilization techniques are currently under investigation and will be discussed in a forthcoming publication [31]. On the other hand, as it is seen in the example, a layer of FE can be used along the punch boundary in order to facilitate the enforcement of essential boundary conditions, or contact boundary conditions in a more general case. 7 Conclusions The examples demonstrate the applicability of the proposed methodology for the coupling of SPH and FE. 17 Apart from the obvious advantages of SPH for adaptivity processes, the use of an SPH discretization in the regions where severe distortions are present allows to go beyond FE computations. If needed, the “connectivity” between particles in the SPH region can be updated during the computations, without projection errors. In this situation, the continuous blending of SPH and FE inherits the stability properties of the SPH considered formulation. On the other hand, the use of a layer of FE along the boundary, coupled with the SPH interpolation, facilitates the enforcement of essential boundary conditions, or contact boundary conditions. The coupling of FE and mesh-free methods maintains the low computational cost of both methods, and takes advantage of their best properties. This method seems to be very promising for the solution of fast-transient dynamics problems, impact simulations, brittle fracture or explosions. Finally, the integration correction avoids, or at least attenuates, the presence of the possible instabilities in the SPH numerical solution (due to the particle integration). References [1] W. K. Liu, T. Belytschko, J. T. Oden, editors, Meshless methods, Comput. Methods Appl. Mech. Eng. 139 (1–4) (1996) 1–440. [2] J. S. Chen, W. Liu, editors, Meshless methods: Recent advances and new applications, Comput. Methods Appl. Mech. Eng. 193 (12–14) (2004) 933–1321. [3] T. Belytschko, Y. Y. Lu, L. Gu, Element free galerkin methods, Int. J. Numer. Methods Eng. 37 (2) (1994) 229–256. [4] W. K. Liu, S. Jun, S. Li, J. Adee, T. Belytschko, Reproducing kernel particle methods for structural dynamics, Int. J. Numer. Methods Eng. 38 (10) (1995) 1655–1679. [5] L. Lucy, A numerical approach to the testing of the fission hypothesis, Astron. J. 82 (1977) 1013–1024. [6] J. W. Swengle, D. L. Hicks, S. W. Attaway, Smoothed particle hydrodynamics stability analysis, J. Comput. Phys. 116 (1995) 123–134. [7] G. R. Johnson, S. R. Beissel, Normalized smoothing functions for SPH impact computations, Comput. Methods Appl. Mech. Eng. 39 (16) (1996) 2725–2741. [8] J. P. Vila, On particle weighted methods and smooth particle hydrodynamics, Math. Models Methods Appl. Sci. 9 (2) (1999) 161–209. 18 [9] J. Bonet, T. Lok, Variational and momentum preservation aspects of smooth particle hydrodynamic formulations, Comput. Methods Appl. Mech. Eng. 180 (1-2) (1999) 97–115. [10] S. Jun, W. K. Liu, T. Belytschko, Explicit reproducing kernel particle methods for large deformation problems, Int. J. Numer. Methods Eng. 41 (1998) 137–166. [11] T. Belytschko, D. Organ, Y. Krongauz, A coupled finite element–element-free Galerkin method, Comput. Mech. 17 (3) (1995) 186–195. [12] D. Hegen, Element free galerkin methods in combination with finite element approaches, Comput. Methods Appl. Mech. Eng. 135 (1-2) (1996) 143–166. [13] A. Huerta, S. Fern´andez-M´endez, Enrichment and coupling of the finite element and meshless methods, Int. J. Numer. Methods Eng. 48 (11) (2000) 1615–1636. [14] G. J. Wagner, W. K. Liu, Hierarchical enrichment for bridging scales and meshfree boundary conditions, Int. J. Numer. Methods Eng. 50 (3) (2001) 507–524. [15] T. Belytschko, A. Huerta, S. Fern´andez-M´endez, T. Rabczuk, Encyclopedia of Computational Mechanics, Vol. 1, John Wiley & Sons, 2004, Ch. 10 (Meshless methods). [16] A. Huerta, S. Fern´andez-M´endez, W. K. Liu, A comparison of two formulations to blend finite elements and mesh-free methods, Comput. Methods Appl. Mech. Eng. 193 (12-14) (2004) 1105–1117. [17] S. Attaway, M. Heinstein, J. Swegle, Coupling of smoothed particle hydrodynamics with the finite element method, Nuclear Eng. Design 150 (1994) 199–205. [18] G. R. Johnson, R. A. Stryk, S. R. Beissel, SPH for high velocity impact computations, Comput. Methods Appl. Mech. Eng. 39 (16) (1996) 2725–2741. [19] J. Bonet, S. Kulasegaram, Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations, Int. J. Numer. Methods Eng. 47 (6) (2000) 1189–1214. [20] W. K. Liu, Y. Chen, S. Jun, J. S. Chen, T. Belytschko, C. Pan, R. A. Uras, C. T. Chang, Overview and applications of the reproducing kernel particle methods, Arch. Comput. Methods Engrg. 3 (1) (1996) 3–80. [21] G. J. Wagner, W. K. Liu, Application of essential boundary conditions in meshfree methods: a corrected collocation method, Int. J. Numer. Methods Eng. 47 (8) (2000) 1367–1379. [22] A. Huerta, S. Fern´andez-M´endez, P. D´ıez, Enrichissement des interpolations d´ ´el´ements finis en utilisant des m´ethodes de particules, ESAIM-Math. Model. Numer. Anal. 36 (6) (2002) 1027–1042. [23] S. Fern´andez-M´endez, P. D´ıez, A. Huerta, Convergence of finite elements enriched with meshless methods, Numer. Math. 96 (1) (2003) 43–59. 19 [24] S. Fern´andez-M´endez, A. Huerta, Imposing essential boundary conditions in mesh-free methods, Comput. Methods Appl. Mech. Eng. 193 (12-14) (2004) 1257–1275. [25] J. Bonet, S. Kulasegaram, M. Rodr´ıguez-Paz, Variational formulation for the smooth particle hydrodynamics (SPH) simulation of fluid and solid problems, Comput. Methods Appl. Mech. Eng. 193 (12-14) (2004) 1245–1256. [26] S. Fern´andez-M´endez, J. Bonet, On the regularity and convergence of the integration correction in SPH. In preparation. [27] J. Bonet, R. D. Wood, Nonlinear continuum mechanics for finite element analysis, Cambridge University Press, Cambridge, 1997. [28] P. Randles, L. Libersky, Normalized SPH with stress points., Int. J. Numer. Methods Eng. 48 (10) (2000) 1445–1462. [29] P. W. Randles, A. G. Petschek, L. D. Libersky, C. T. Dyka, Stability of DPD and SPH., in: M. Griebel, M. A. Schweitzer (Eds.), Meshfree methods for partial differential equations, Vol. 26 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin, 2002, pp. 117–129, papers from the International workshop, Universit¨at Bonn, Germany, September 11-14, 2001. [30] T. Belytschko, S. Xiao, Stability analysis of particle methods with corrected derivatives., Comput. Math. Appl. 43 (3-5) (2002) 329–350. [31] Y. Vidal, Mesh-free methods for dynamic problems: incompressibility and large strains, Ph.D. thesis, Universitat Polit`ecnica de Catalunya, Barcelona, Spain (2005). 20