A modified finite particle method: multi-dimensional statics and dynamics

691 III International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2013 M. Bischoff, E. O˜nate, D.R.J. Owen, E. Ramm & P. Wriggers (Eds) A MODIFIED FINITE PARTICLE METHOD: MULTI-DIMENSIONAL STATICS AND DYNAMICS DOMENICO ASPRONE1, FERDINANDO AURICCHIO2,3, ANDREA MONTANINO2,3AND ALESSANDRO REALI2,3 1Dipartimento di Ingegneria Strutturale (DIST) Universit`a degli Studi di Napoli, via Claudio 21, 80125 Napoli (Italy) 2Dipartimento di Ingegneria Civile ed Architettura (DICAr) Universit`a degli Studi di Pavia, via Ferrata 3, 27100 Pavia (Italy) 3Centro di Simulazione Numerica Avanzata (CeSNA) Istituto Universitario di Studi Superiori, piazza della Vittoria 15, 27100 Pavia (Italy) Key words: Particle methods, collocation methods, projection methods, elasto-dynamics Abstract. We present the so-called Modified Finite Particle Method (MFPM), that is a recent methodology of approximation of differential operators, based on the projection of the Taylor series of a function u(x) on a set of projection functions. In particular, we discuss the generalization of MFPM formulation to the multi-dimensional case, extending the methodological procedure adopted for the one-dimensional case. Moreover, we address the extension to dynamics and solve problems with an explicit time integration scheme. Finally, we apply the MFPM to an elasto-statics (a perforated plate under tension) and two elasto-dynamics (a two-dimensional bar under a quasi-impulsive load and a quarter of an annulus under a sinusoidal body load) benchmarks. When an analytical solution is available we calculate the corresponding convergence orders of the error, always obtaining the expected second-order accuracy 1 INTRODUCTION In recent years, meshless numerical methods have become increasingly important due to their characteristic of being totally free of grids or meshes. This is particularly useful when dealing with the numerical simulation of problems implying large deformations or high velocity impacts, where classical element-based methods (e.g., the Finite Element Method) suffer from pathologies like excessive element distortion, spurious numerical errors, and mesh sensitivity. On the contrary, meshless methods overcome these difficulties since the nodes are not “rigidly” connected between themselves. As a consequence, the approximation is carried out taking into account the current distribution of the particles. 1 A modified finite particle method: Multi-dimensional statics and dynamics 692 Domenico Asprone, Ferdinando Auricchio, Andrea Montanino and Alessandro Reali Among meshless methods, we focus in particular on Particle Methods. The first proposed particle method has been the Smoothed Particle Hydrodynamics (SPH), introduced by Lucy [1] and Gingold and Monaghan [2], for the study of astrophysical problems. In this method the whole domain Ω is discretized in a finite number of particles, each one characterized by a given mass, velocity, and energy. Therefore the state of a continuum is represented by the properties of the particles, and as a consequence SPH is particularly suitable for the description of the behaviour of fluids and ideal gases. The starting point, in 1D, is the identity f(xi)=Ω f(x)δ(x−xi)dx (1) where δ(x−xi) is the Dirac Delta function. The first approximation consists of replacing the Dirac Delta function with a bell-shaped smooth function W(x−xi,h), called smoothing function or kernel function: f(xi)≃Ω f(x)W(x−xi,h)dx (2) then the derivatives are computed according to f(n)(xi)≃(−1)nΩ f(x)W(n)(x−xi)dx (3) Equation (2) is referred as the kernel evaluation of f(x). The amount his the smoothing length, that is the distance where the smoothing function is significantly different from zero. In order to numerically evaluate the integrals of Equations (2) and (3), the whole domain is divided into a number of subdomains, each associated to a particle xj; then the integrals are replaced by summations over j. In proximity of the boundary, Equations (2) and (3) are inaccurate because the smoothing function is not completely developed: as a consequence, many other method have been developed to overcome this deficiency. One of the first has been introduced by Liu et al. [3], that introduced a polynomial corrective function in order to restore the effectiveness of Equation (2) also at the boundary, but leaving unsolved the inaccuracy of Equation (3). Chen et al. [4], in the Corrective Smoothed Particle Method (CSPM) achieve accurate formulas for the kernel evaluation by manipulating Taylor series expansion up to the zero-th order, then the following approximations are carried out by considering the following terms of the Taylor series, and using the previous approximations. In the Modified Smoothed Particle Hydrodynamics (MSPH) [5], the unknown function and its derivatives are evaluated at the same time. This procedure prevents the method from error propagation, but increases the computational cost, since at each particle a system of equations has to be solved. In Asprone et al. [6] a Modified Finite Particle Method (MFPM) has been developed starting from the MSPH [5]. The kernel evaluation of the unknown function is not computed, since it is assumed that the value of the function is not a real unknown of the 2 693 Domenico Asprone, Ferdinando Auricchio, Andrea Montanino and Alessandro Reali problem. Moreover a 1D elasto-statics test is solved using this method as well as others taken from the literature (the original SPH formulation [1], the RKPM [3], the CSPM [4]) and the results are compared. In Asprone et al. [7] MFPM is extended to a 2D scalar problem, while in Asprone et al. [8] the solution of some 2D elasto-statics and elasto-plastics problems are presented. In this work we recall the 1D formulation of MFPM. We introduce the extension of the MFPM to the multi-dimensional case, and give an analytical expression of the discrete form of the differential operators, then we recall the equations of structural elastodynamics and show their discrete form by MFPM, then we show the treatment of the boundary conditions for explicit time integration. We apply the method to some elastostatics and elasto-dynamics tests. Finally we draw some conclusions and illustrate the further developments of the method. 2 MODIFIED FINITE PARTICLE METHOD: MULTI-DIMENSIONAL FORMULATION In the following we show the extension of the Modified Finite Particle Method to the three-dimensional case, recalling the same steps of the 1D procedure presented in Asprone et al. [6]. We consider the Taylor series expansion of an unknown function u(x x x) up to the second order and multiply it by 9 projection functions Wk i=Wk(x x x−x x xi), k=1, ..., 9; then, integrating over the domain, we obtain a set of equations of the type Dxu(x x xi)Ω (x−xi)Wk idV +Dyu(x x xi)Ω (y−yi)Wk idV +Dzu(x x xi)Ω (z−zi)Wk idV + 1 2D2 xxu(x x xi)Ω (x−xi)2Wk idV +1 2D2 yyu(x x xi)Ω (y−yi)2Wk idV + 1 2D2 zzu(x x xi)Ω (z−zi)2Wk idV +D2 xyu(x x xi)Ω (x−xi)(y−yi)Wk idV + D2 yzu(x x xi)Ω (y−yi)(z−zi)Wk idV +D2 xzu(x x xi)Ω (x−xi)(z−zi)Wk idV = Ω (u(x x x)−u(x x xi))Wk idV (4) 3 694 Domenico Asprone, Ferdinando Auricchio, Andrea Montanino and Alessandro Reali that can be rewritten in matrix form as: A A Ai               Dxu(x x xi) Dyu(x x xi) Dzu(x x xi) D2 xxu(x x xi) D2 yyu(x x xi) D2 zzu(x x xi) D2 xyu(x x xi) D2 yzu(x x xi) D2 xzu(x x xi)               =               �Ω(u(x x x)−u(x x xi)) W1 idV �Ω(u(x x x)−u(x x xi)) W2 idV �Ω(u(x x x)−u(x x xi)) W3 idV �Ω(u(x x x)−u(x x xi)) W4 idV �Ω(u(x x x)−u(x x xi)) W5 idV �Ω(u(x x x)−u(x x xi)) W6 idV �Ω(u(x x x)−u(x x xi)) W7 idV �Ω(u(x x x)−u(x x xi)) W8 idV �Ω(u(x x x)−u(x x xi)) W9 idV               (5) The choice of the projection functions is performed with the only requirement that, at each particle, the matrix A A Aiis non singular. In our tests we choose      W1 i=x−xi W2 i=y−yi W3 i=z−zi      W4 i=(x−xi)2 W5 i=(y−yi)2 W6 i=(z−zi)2     W7 i=(x−xi)(y−yi) W8 i=(y−yi)(z−zi) W9 i=(x−xi)(z−zi) The domain is divided into finite subdomains ∆Vj, one for each particle x x xj, according to the Voronoi tessellation procedure; for each particle an influence region Ωiis also defined, depending, in SPH-based methods, on the smoothing length. In MFPM we do not define a fixed value of the smoothing length, but we prefer to set the number of particles to be considered for the approximation of derivatives. For all the particles such that xj/∈Ωiwe pose that Wk i(x x x=x x xj) = 0 for all k=1, ..., 9. Then, the integrals are approximated with summations, and the approximation schemes for the first and second spatial derivatives at each particle xiare obtained by inverting (5). Therefore, we are able to write the analytical approximation schemes of the spatial derivatives. In particular, the approximation schemes of the first derivatives are:                        Dx,ij = 9 � k=1 Ei 1k�Wk ij∆Vj−δij � h Wk ih∆Vh� Dy,ij = 9 � k=1 Ei 2k�Wk ij∆Vj−δij � h Wk ih∆Vh� Dz,ij = 9 � k=1 Ei 3k�Wk ij∆Vj−δij � h Wk ih∆Vh� (6) 4 695 Domenico Asprone, Ferdinando Auricchio, Andrea Montanino and Alessandro Reali and the schemes for the second derivatives are:                                                          D2 xx,ij = 9 � k=1 Ei 4k�Wk ij∆Vj−δij � h Wk ih∆Vh� D2 yy,ij = 9 � k=1 Ei 5k�Wk ij∆Vj−δij � h Wk ih∆Vh� D2 zz,ij = 9 � k=1 Ei 6k�Wk ij∆Vj−δij � h Wk ih∆Vh� D2 xy,ij = 9 � k=1 Ei 7k�Wk ij∆Vj−δij � h Wk ih∆Vh� D2 yz,ij = 9 � k=1 Ei 8k�Wk ij∆Vj−δij � h Wk ih∆Vh� D2 xz,ij = 9 � k=1 Ei 9k�Wk ij∆Vj−δij � h Wk ih∆Vh� (7) where E E Eiis the inverse of A A Ai A 2D formulation of the method is simply achieved by considering only the derivatives in the xand ydirections, and k=1,2,4,5,7 from equation (4) to equation (7). The three-dimensional subdomains ∆Vjare replaced by planar subdomains ∆Aj. 3 ELASTICITY In the following we introduce the linear elasto-dynamics problem in the three-dimensional space and show how it can be approximated with the Modified Finite Particle Method. We consider an elastic body on a domain Ω, subjected to internal forces b b b, constrained displacements ¯ s s son the Dirichlet boundary ΓDand the traction ¯ t t ton the Neumann boundary ΓN. The equations that govern the problem are                      ρ∂2s s s ∂t2=∇·σ σ σ+b b bx x x∈Ω σ σ σn n n=¯ t t t(t)x x x∈ΓN s s s=¯ s s s(t)x x x∈ΓD s s s(x x x, 0) = s s s0(x x x)x x x∈Ω ∂s s s ∂t(x x x, 0) = ˙ s s s0(x x x)x x x∈Ω (8) where ρis the mass density of the material, n n nis the outward normal vector at the boundary, s s s=s s s(x x x, t) is the vectorial displacement field, whose components are u=u(x x x, t), 5 696 Domenico Asprone, Ferdinando Auricchio, Andrea Montanino and Alessandro Reali v=v(x x x,t) and w=w(x x x, t), σ σ σ=C(∇s s s)Sis the symmetric Cauchy stress tensor, Cis the fourth order linear elastic isotropic tensor, with components Cijkl =λδijδkl +µ(δikδjl +δilδjk)(9) where λand µare the Lam´e constants, which can be expressed in terms of the Young modulus Eand the Poisson ratio ν. By making explicit Equation (8) with respect to the components of the displacement u,vand w, we obtain      ρu,tt =(λ+2µ)u,xx +µ(u,yy +u,zz )+(λ+µ)(v,xy +w,xz )+bx ρv,tt =(λ+2µ)v,yy +µ(v,xx +v,zz )+(λ+µ)(u,xy +w,yz )+by ρw,tt =(λ+2µ)w,zz +µ(w,xx +w,yy )+(λ+µ)(u,xz +v,yz )+bz (10) where the subscripts preceded by a comma stand for the operation of partial derivative. In the spirit of collocation methods, we enforce the discrete form of Equation (10) for internal particles, and the discrete form of the boundary conditions on the boundary particles, according to the approximation schemes (6) and (7). We obtain a linear system in the form K K Kˆ s s s=f f f(11) where the components of f f fare ρ¨ ˆ s s s−b b bfor the rows associated to internal particles, and ¯ s s s or ¯ t t tin case of Dirichlet or Neumann boundary particles respectively. For elasto-statics applications, the form (11) is sufficient to solve the problem: in this case the time derivative is zero, and the system can be inverted; therefore both the internal and external particle values are found at once. In case of elasto-dynamics, we have to discretize also the time derivative. We choose an explicit second order scheme ¨ ˆ s s sn=ˆ s s sn+1 −2ˆ s s sn+ˆ s s sn−1 ∆t2(12) where ∆tis the time step. The equations of the system (11), collocated on internal particles, become � j Kij ˆsn j=ρˆsn+1 i−2ˆsn i+ˆsn−1 i ∆t2−bn i(13) while the equations collocated on boundary particles, that do not undergo the time derivative, are in the form � j Kij ˆsn+1 j=¯un+1 i(14) Equations (14) cannot be solved by explicit time integration, since the values of ˆsn+1 j may depend, in case of Neumann boundary conditions, on the values of the internal 6 697 Domenico Asprone, Ferdinando Auricchio, Andrea Montanino and Alessandro Reali particles at the same time step n+1. To overcome this difficulty, we perform a static condensation of K K K, and separate the equations collocated on internal particles from the ones collocated on the boundary. Even the degrees of freedom are separated, and so the final form of Equation (11) is �K K KII K K KIB K K KBI K K KBB ��ˆ s s sI ˆ s s sB�=�ρ¨ ˆ s s sI−b b bI ¯ u u u�(15) where ¯ u u uis the vector of the imposed displacements or the imposed stresses at the boundary, and K K KII,K K KIB,K K KBI,K K KBB are the minors of the matrix K K Krelated to the internal and boundary particles. From the second set of equations of (15) we compute u u uB=K K K−1 BB(¯ u u u−K K KBIˆ s s sI)(16) then we replace it into the first set of (15), and obtain (K K KII −K K KIBK K K−1 BBK K KBI)ˆ s s sI=ρ¨ ˆ s s sI−K K KIBK K K−1 BB¯ u u u−b b bI(17) where the amount K K KII −K K KIBK K K−1 BBK K KBI is the modified stiffness matrix, namely ˜ K K K. Equations (17) form an unconstrained ordinary differential equation system which can be solved with a suitable time integration scheme, even explicit, (e.g. Equation (12)), taking care to respect the eventual limitations of the time step. The unknowns of this system are the values of the unknown functions at the internal particles. Equation (16) can be used to retrieve the values of the functions at the boundary particles. 4 APPLICATIONS In the following we apply the mmethod presented in the earlier sections. First we introduce a 2D statics problem, that is, the classical test of an infinitely extended plate with a central hole under a uniform remote stress. Regarding dynamics problems, we show the wave propagation in a two-dimensional bar, and a quarter of annulus under a sinusoidal body load. 4.1 A quarter of perforated plate The geometry of this problem is depicted in Figure 1. The radius of the internal hole is a=0.2 The equations that govern the problem are the 2D statics version of (10). The boundary conditions are      σ σ σn n n=0 0 0on Γ1and Γ4 σ σ σn n n·t t t= 0 and s s s·n n n= 0 on Γ2and Γ5 σ σ σn n n=[σ00]Ton Γ3 (18) 7 698 Domenico Asprone, Ferdinando Auricchio, Andrea Montanino and Alessandro Reali Ox yaa σ0 Γ1 Γ2 Γ3 Γ4 Γ5 Figure 1: The model of a quarter of square plate with a central hole where n n nis the outward normal, t t tis the unit vector tangent to the boundary, and σ0is the remote stress. We solve the problem considering a reduced domain, on which we impose boundary conditions according to the exact solution, that is, in terms of stresses, σxx =σ01−a2 r23 2cos 2θ+ cos 4θ+3a4 2r4cos 4θ(19a) τxy =σ0−a2 r21 2sin 2θ+ sin 4θ+3a4 2r4sin 4θ(19b) σyy =σ0−a2 r21 2cos 2θ−cos 4θ−3a4 2r4cos 4θ(19c) where (r, θ) are the polar coordinates, θbeing measured from the positive x-axis counterclockwise. We introduce the Stress Intensity Factor (SIF), that is the ratio between the maximum σxx and the remote stress σ0. For the test under consideration, the analytical solution provides SIF = 3. We solve the problem with the following values of the data E= 100000 ,ν=0.33 ,σ 0= 100 (20) and then we compare the analytical value of the SIF with the numerical results. The distribution of σxx stress is shown in Figure 2(a). The error is computed as er=|SIFan −SIFnum| |SIFan|=|3−SIFnum| 3(21) The convergence diagram of the error is shown in Figure 2(b), where Nis the total number of particles used for the numerical solution. 8 699 Domenico Asprone, Ferdinando Auricchio, Andrea Montanino and Alessandro Reali 0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1 0 50 100 150 200 250 300 x y (a) σxx in a square with a central hole 10 20 40 80 160 320 0.01 0.1 1 10 Original MFPM 2nd order slope √N er (b) Logaritmic convergence diagram of the error of the SIF in a square plate with a central hole Figure 2: The quarter of perforated plate: numerical solution and convergence diagram 4.2 Dynamics of a bar under quasi-impulsive load We consider a two-dimensional bar under a quasi-impulsive load. The geometry is depicted in figure 3, where L= 1 and H=0.2. σ(t) L Hx y O Γ1Γ2 Γ3 Γ4 Figure 3: Geometry of the bar under quasi-impulsive load The equations that govern the problem are the 2D plane strain version of (10); the boundary conditions are      u= 0 and v= 0 on Γ1 σyy = 0 and τxy = 0 on Γ2and Γ4 σxx =σ(t) and τxy = 0 on Γ3 (22) where σ(t)=σ0exp �−b(t−t0)2�is the quasi-impulsive load on the right side of the bar; the test has been performed considering a Poisson ratio equal to zero, so to reproduce a one-dimensional test. We also set E= 100 and ρ= 100. For this test an analytical solution is available for σxx(x, y), since the analytical propagation speed c=�E/ρ is known. The missing data of this problem are: σ0= − 100, b= 100, t0=0.3(23) 9