An improvement for singularity in 2-D corrective smoothed particle method with application to Poisson problem

An improvement for singularity in 2-D corrective smoothed

particle method with application to Poisson problem

Citation for published version (APA):

Hou, Q., Tijsseling, A. S., & Mattheij, R. M. M. (2010). An improvement for singularity in 2-D corrective smoothed particle method with application to Poisson problem. (CASA-report; Vol. 1055). Technische Universiteit

Eindhoven.

Document status and date: Published: 01/01/2010

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

providing details and we will investigate your claim.

EINDHOVEN UNIVERSITY OF TECHNOLOGY

Department of Mathematics and Computer Science

CASA-Report 10-55

September 2010

An improvement for singularity in 2-D corrective

smoothed particle method with application

to Poisson problem

by

Q.Z. Hou, A.S. Tijsseling, R.M.M. Mattheij

Centre for Analysis, Scientific computing and Applications

Department of Mathematics and Computer Science

Eindhoven University of Technology

P.O. Box 513

5600 MB Eindhoven, The Netherlands

ISSN: 0926-4507

An improvement for singularity in 2-D

corrective smoothed particle method with

application to Poisson problem

Q.Z. Hou, A.S. Tijsseling, R.M.M. Mattheij

AbstractThe corrective smoothed particle method (CSPM) is an improvement to the conventional smoothed particle hydrodynamics (SPH) method. It overcomes the problems of tensile instability and boundary deficiency in SPH. As a specially designed numerical method for time-dependent problems, it has been applied to transient heat conduction and structural dynamics problems. In this paper the 2-D CSPM is used to solve time-independent problems, where an extra numerical in-version needs special attention. To permit an invertible matrix for small smoothing length and hence enhance the computational efficiency, an improvement is proposed to the original CSPM. The accuracy and efficiency of the CSPM and the proposed modification is demonstrated by the numerical example.

1 Introduction

Smoothed particle hydrodynamics (SPH) is a meshfree, transient, Lagrangian par-ticle approach proposed for solving 3D astrophysical problems in 1977. It was ex-tended to solid mechanics at the beginning of the 1990s and from then on a number of researchers were attracted to this new field and consequently many new mesh-free methods were developed for an increasing number of applications. The history, development and applications can be found in the book by Liu and Liu [7] and the most recent reviews papers [8, 9].

As is well known, the SPH method does not have zero-order consistency in the boundary region and it suffers from tensile instability problem encountered in elas-todynamics. To address these two shortcomings, Chen et al. [2] proposed a correc-tive smoothed particle method (CSPM). Like SPH, the CSPM is a specially designed method for time-dependent problems. However, as an approximation technique for

Q.Z. Hou, A.S. Tijsseling, R.M.M. Mattheij

Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands e-mail: q.hou@tue.nl, a.s.tijsseling@tue.nl, r.m.m.mattheij@tue.nl

2 Q.Z. Hou, A.S. Tijsseling, R.M.M. Mattheij

spatial operators, it is suitable for time-independent problems as well. The purpose of this paper is to introduce the meshfree CSPM to solve time-independent Poisson problems.

The rest of this paper is organized as follows. As CSPM is developed from SPH, Section 2 briefly reviews the SPH method. The CSPM in two space dimensions is presented in Section 3. An improvement to the original CSPM is proposed to release the restriction on the smoothing length in relation to non-singularity of the interpo-lation matrix. In Section 4 the original and modified CSPM are applied to solve the Poisson equation, and the accuracy and efficiency are compared and discussed. Due to using small smoothing length, the modified CSPM may not only improve the numerical accuracy, but also increase the computational efficiency. Section 5 gives some conclusions.

2 Smoothed particle hydrodynamics

In the SPH method, a function f (x) is approximated by the integral ˆ

f (x) :=

∫

Ωf (ξ)W dξ, (1)

where ˆf (x) is the approximated value of f (x), W := W (x−ξ, h) is the smoothing or kernel function, and h is a dilation parameter named smoothing length. It determines the influence domain or support of the kernel function W with respect toξ.

The gradient of f (x) is approximated by ∇ f (x) :=∫

Ωf (ξ)∇Wdξ, (2)

where∇W := ∇W(x −ξ, h). Definitions (1) and (2) are generally called kernel es-timates. The approximation accuracy of (1) and (2) largely depends on the ker-nel function W and the smoothing length h. The kerker-nel function W has to fulfill a number of conditions such as normalization, compactness, positivity, symmetry, and etc [8]. In SPH, the most widely used kernel function in a two-dimensional domain Ω is the cubic spline function

W (q) := 15 7πh2      2/3− q2+ q3/2, 06 q < 1, (2− q)3/6, 16 q < 2, 0, q> 2, (3)

where q :=|x −ξ|/h is the relative distance between two points.

For numerical analysis, the continuous integrals in (1) and (2) are approximated by proper numerical quadrature rules. The value of ˆf (x) at xi, i.e. ˆfi:= ˆf (xi), is

An improvement for singularity in 2-D corrective smoothed particle method 3 ˆ fi=. N

∑

where fj:= f (xj), Wi j:= W (xj− xi, h). The particle mass is computed from its

volume Vjand densityρjas mj=ρjVj, and N is the number of particles in domain

Ω. Similarly, kernel estimates of the derivatives (2) are approximated as fi,α=. N

∑

∂xα , andα denotes the direction index, the range of which equals the dimension of the domain. Note that (4) and (5) can be interpreted as interpolations based on nodes j weighted by Wi j

mj

ρj and Wi j,α

mj

ρj, re-spectively. Consequently, SPH is often referred to as a weighted interpolation tech-nique.

3 Corrective smoothed particle method

Based on Taylor series expansions and kernel estimates in SPH, a corrected kernel estimates have been developed by Chen et al. [2]. It is designed typically for time-dependent problems and has been applied to transient heat conduction [2], elasto-dynamics [4], Burgers’ equation [1] and elastoplastic elasto-dynamics [5]. The CSPM in two space dimensions is described below. To release limitations on choosing small smoothing length h and also reduce the bandwidth of the system matrix, an improve-ment is proposed herein.

3.1 Original CSPM

The following equation is the backbone in CSPM for generating the derivative ap-proximations of a function at any point xi= (x1i, x2i) in a 2-D domainΩ,

∫ Ωf (x)Wi,γdΩ = fi ∫ ΩWi,γdΩ + fi,α ∫ Ω(x α_{− x}α i )Wi,γdΩ + fi,αβ 2! ∫ Ω(x α_{− x}α i)(xβ− xβi)Wi,γdΩ + ··· , (6) where Wi,γ:=∂W(x−x_∂xγi,h) and fi,αβ := ∂ 2_{f (x} i)

∂xα_∂xβ; the superscriptsα,β,γ= 1, 2 denote spacial components. Repetition of a spatial index implies summation over its range in (6) only, which is derived by expanding the Taylor series for f around xi,

4 Q.Z. Hou, A.S. Tijsseling, R.M.M. Mattheij

Neglecting the second-order and higher derivative terms on the RHS of (6) and approximating the integrals by Riemann sums, the following equations for the two first derivative approximations are obtained:

A_αβif_i,β = F_αi, (7) in which A_αβi:= N

∑

∑

Notations xji:= xj− xiand fji:= fj− fiare used. Substituting Wi,γδ =∂ 2_{W (x−x}

i,h)

∂xγxδ

for Wi,γin (6) and neglecting the third-order and higher derivative terms, the

follow-ing three equations for the second derivative approximations are obtained:

B_ηθig_θi=ϕ_ηi, (8)

in which B_ηθi:=ζ N

∑

∑

∑

The symbols B_ηθi, g_θiandϕ_ηirepresent B_γδαβi, f_i,αβandϕ_γδi, respectively. Integers

η,θ = 1 to 3, and α,β,γ,δ = 1 to 2. The index η corresponds to index γδ as follows: 1↔ 11, 2 ↔ 22 and 3 ↔ 12, as doesθtoαβ. The constantζ= 1−δ_αβ/2, whereδ_αβis the Kronecker delta. Elements A_αβiand B_ηθiconstitute a 2×2 matrix Aiand a 3×3 matrix Bi, respectively. For the reason mentioned at the end of Section

2, we call Aiand Biinterpolation matrices.

The corrected kernel estimates of the first- and second-order derivatives are solved from (7) and (8), respectively. Apparently, matrices Ai and Bi are not

al-lowed to be singular. To ensure nonsingularity, every particle must at least have three neighbours in its support. This requirement is particularly true for particles in the boundary region. The reason can be explained by using the evenly spaced particles in a rectangular domain as shown in Fig. 1.

For the particle i at the corner, three quarters of its circular support fall outside the domain. Consequently, to ensure at least three other particles in its support, the constant radius R := 2h must be larger than√2∆ (i.e. h >√2∆/2), where ∆ is the distance between two adjacent particles in one direction (see Fig. 1). Under this condition at least 5 and 8 other particles are in the support of the edge particle j and the interior particle k, respectively, as shown in Fig. 1.

An improvement for singularity in 2-D corrective smoothed particle method 5

Fig. 1 Uniformly distributed particles in a 2-D domain and their supports: i – corner particle, j – edge particle, k – interior particle.

3.2 Modified CSPM

It is reported that the smaller smoothing length is used, the better results are promised [2, 3]. It is also obvious that the less particles are involved in the ker-nel approximation, the smaller bandwidth of the system matrix will be generated and accordingly less computational effort is needed, which is especially remarkable in 2- and 3-D problems. Moreover, If central differences are employed in FDM to approximate the Laplacian, a 5-point stencil is obtained [6]. Accordingly, one may ask whether it is possible to weaken the requirement of the inversibility of matrices Bi such that only 4 neighbouring particles are needed in the support of particle k

(green circle in Fig. 1). Much effort has been put into obtaining a positive answer. We found that for h∈ (∆/2,√2∆/2] the singularity arises from the mixed derivative g3iin (8). By excluding it from gθiand changing the other two terms accordingly,

two equations are generated as e

B_ηθieg_θi= eϕ_ηi. (9)

Nowη,θ= 1 to 2, and the other terms are the same as defined before. This is the modified 2-D CSPM for the second derivative approximations proposed herein.

3.3 Numerical implementation

The implementation of the CSPM for a 2-D Poisson problem will be described. First, for a given particle i we establish the relationship between the second deriva-tive approximations fi,αβand the nodal unknowns fjby solving the set of equations

(8) or (9), where fi,α are known from (7). Second, the relationship between the

6 Q.Z. Hou, A.S. Tijsseling, R.M.M. Mattheij

and fi,yy. Third, the system matrix K is assembled according to a given numbering

of the particles. The RHS term p consists of the nodal values of a given function. Fourth, the boundary conditions are enforced by appropriately modifying rows of K and p. Then the system equation with imposed boundary conditions is completed as

Kf = p. (10)

Finally, the nodal unknowns f are solved from the sparse linear system (10) by any suitable numerical method. Gauss elimination is used herein.

4 Results and discussion

In order to demonstrate the efficiency and accuracy of the CSPM and the proposed improvement, a standard Poisson equation is considered and the results are com-pared with the analytical solutions.

For a classic Poisson problem with a unit square domain, the RHS function is

5

4ex+y/2. The following Dirichlet boundary conditions: u(x, 0) = ex, u(x, 1) = ex+1/2,

u(0, y) = ey/2and u(1, y) = ey/2are applied to the boundaries. The analytical solu-tion of this problem is u(x, y) = ex+y/2. In the simulation, 40 particles are evenly spaced on the boundary and 81 particles are uniformly distributed in the interior. The numerical simulations are depicted in Fig. 2, where the exact solutions and the FDM simulations are plotted as well. It is clear as shown in Fig. 2 that all the re-sults agree very well with each other. The differences can not be distinguished from Fig. 2.

The numerical results at different points are shown in Table 1. For comparison, the exact values and relative error are also presented. It is noted that the CSPM results have excellent agreements with the exact solutions. The relative errors for the concerned positions are less than 0.02%.

Table 1 Comparison of CSPM results with the exact solution for a 2-D Poisson problem. Method\ (x, y) (0.1, 0.4) (0.2, 0.6) (0.5, 0.5) (0.7, 0.8) CSPM 1.3500 1.6490 2.1173 3.0043 Exact 1.3499 1.6487 2.1170 3.0042 Errora_{(%) 0.0134 0.0156 0.0151 0.0052} a_Error=|uei−uni| |ue i| × 100

The effects of smoothing length h and number of particle N on the numerical accuracy and convergence rate for CSPM (Section 3.1) are listed in Table 2. As mentioned in Section 3, numerical results do not exist for r := h/∆ <√2/2 because of singularity of Bi. It is clear that the max-norm error is decreasing for increasing

An improvement for singularity in 2-D corrective smoothed particle method 7 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1 2 3 4 5 x (a) y u 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1 2 3 4 5 x (b) y u 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1 2 3 4 5 x (c) y u 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1 2 3 4 5 x (d) y u

Fig. 2 Comparison of numerical simulations by the (a) FDM and CSPM with (b) h = 0.8∆ and (c)

h = 1.1∆ with (d) the exact solution.

more or less the same. For a particle number N = 11× 11, the CPU time is 8.2 ms for r = 0.9 and 9.1 ms for r = 1.1. The reason is that the bandwidth of system matrix Kis changed from 9 (r = 0.9) to 13 (r = 1.1) as the smoothing length increases; see Fig. 1.

Table 2 The effects of h and N on max-norm error and convergence rate of the original CSPM.

r\ N 11×11 21×21 41×41 rate

0.7 – – – –

0.9 2.63× 10−4 6.58× 10−5 1.65× 10−5 2.00 1.1 3.27× 10−4 8.52× 10−5 2.20× 10−5 2.01

The results of the modified CSPM (Section 3.2) are listed in Table 3. Three more points should be noted. First, the smoothing length h can be smaller than√2∆/2, and the corresponding errors for all three Ns are also smaller. Second, comparing

8 Q.Z. Hou, A.S. Tijsseling, R.M.M. Mattheij

with Table 2 for r >√2/2 = 0.707, the numerical errors and convergence rates are the same as the original CSPM. The reason is that very close system matrix K is obtained for the original and modified CSPM. Last, with N = 11× 11 particles, the CPU time for the three cases are 7.4 ms, 7.9 ms and 8.6 ms, respectively, less than that in CSPM. This is because for r = 0.7 the bandwidth of K is 5 in the modified CSPM; for r = 0.9 or 1.1, although the bandwidth is the same as in the original CSPM, the dimension of eBiis 2 instead of 3 of Bi. Therefore, the modified CSPM

not only allows smaller smoothing length h, but is more efficient.

Table 3 The effects of h and N on max-norm error and convergence rate of the modified CSPM.

r\ N 11×11 21×21 41×41 rate

0.7 1.42× 10−4 3.60× 10−5 9.04× 10−5 2.00 0.9 2.63× 10−4 6.58× 10−5 1.65× 10−5 2.00 1.1 3.27× 10−4 8.52× 10−5 2.20× 10−5 2.01

5 Conclusions

In this article, the corrective smoothed particle method (CSPM) is used to solve 2-D Poisson problem. The accuracy and efficiency of the CSPM are demonstrated by one numerical example. It was found that the smaller the smoothing length, the better results are obtained. Moreover, to reduce the restriction on choosing the smoothing length, an improvement is proposed, which increases the numerical accuracy and reduces the bandwidth of the system matrix and hence improves the computational efficiency.

Acknowledgements The authors are grateful to the Chinese Ministry of Education for financially supporting the PhD studies of the first author.

References

1. Chen, J., Beraun, J.: A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems. Comput. Meth. Appl. Mech. Eng. 190, 225–239 (2000)

2. Chen, J., Beraun, J., Carney, T.: A corrective smoothed particle method for boundary value problems in heat conduction. Int. J. Numer. Methods Eng. 46, 231–252 (1999)

3. Chen, J., Beraun, J., Jih, C.: Completeness of corrective smoothed particle method for linear elastodynamics. Comput. Mech. 24, 273–285 (1999)

4. Chen, J., Beraun, J., Jih, C.: An improvement for tensile instability in smoothed particle hydro-dynamics. Comput. Mech. 23, 279–287 (1999)

5. Chen, J., Beraun, J., Jih, C.: A corrective smoothed particle method for transient elastoplastic dynamics. Comput. Mech. 27, 177–187 (2001)

An improvement for singularity in 2-D corrective smoothed particle method 9

6. Leveque, R.: Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. SIAM, Philadelphia (2007)

7. Liu, G., Liu, M.: Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientific, Singapore (2003)

8. Liu, M., Liu, G.: Smoothed particle hydrodynamics (SPH): An overview and recent develop-ments. Arch. Comput. Meth. Eng. 17, 25–76 (2010)

PREVIOUS PUBLICATIONS IN THIS SERIES:

Number Author(s)

Title

Month

10-51

10-52

10-53

10-54

10-55

T.I. Seidman

A. Muntean

R.V. Polyuga

R.V. Polyuga

A. van der Schaft

R.V. Polyuga

A. van der Schaft

Q.Z. Hou

A.S. Tijsseling

R.M.M. Mattheij

Fast-reaction asymptotics

for a time-dependent

reaction-diffusion system

with a growing nonlinear

source term

Discussion on: “Passitivity

and structure preserving

order reduction of linear

port-Hamiltonian systems

using Krylov subspaces”

Structure preserving

moment matching for

port-Hamiltonian systems:

Arnoldi and Lanczos

Model reduction of

port-Hamiltonian systems as

structured systems

An improvement for

singularity in 2-D

corrective smoothed

particle method with

application to Poisson

problem

Sept. ‘10

Sept. ‘10

Sept. ‘10

Sept. ‘10

Sept. ‘10