Lagrangian meshfree finite difference particle method with variable smoothing length for solving wave equations

Abstract

In a Lagrangian meshfree particle-based method, the smoothing length determines the size of the support domain for each particle. Since the particle distribution is irregular and uneven in most cases, a fixed smoothing length sometime brings too much or insufficient neighbor particles for the weight function which reduces the numerical accuracy. In this work, a Lagrangian meshfree finite difference particle method with variable smoothing length is proposed for solving different wave equations. This pure Lagrangian method combines the generalized finite difference scheme for spatial resolution and the time integration scheme for time resolution. The new method is tested via the simple wave equation and the Burgers’ equation in Lagrangian form. These wave equations are widely used in analyzing acoustic and hydrodynamic waves. In addition, comparison with a modified smoothed particle hydrodynamics method named the corrective smoothed particle method is also presented. Numerical experiments show that two kinds of Lagrangian wave equations can be solved well. The variable smoothing length updates the support domain size appropriately and allows the finite difference particle method results to be more accurate than the constant smoothing length. To obtain the same level of accuracy, the corrective smoothed particle method needs more particles in the computation which requires more computational time than the finite difference particle method.

Keywords

Lagrangian approach meshfree method wave equation Burgers’ equation finite difference particle method smoothed particle hydrodynamics

Introduction

Numerical methods have been widely implemented to model acoustic and hydrodynamic waves. Compared with mesh-based methods,^1–4 meshfree methods receive much concern because field points used in this method are arbitrarily distributed and material properties are no more assigned to elements. The element-free Galerkin method,⁵ the method of radial basis function method,⁶ the singular boundary method,⁷ and other meshfree methods^8,9 had been used to address certain wave problems. In this work, we proposed a pure Lagrangian meshfree particle-based approach for solving different time-dependent wave equations.

In recent years, Lagrangian meshfree particle-based methods are fast developing and widely used, since it is easily adapted to modeling fluid dynamic problems with complex or time-evolving boundaries for single or multiple phases like the numerical simulation of dam break flow,¹⁰ flooded city,¹¹ and bubble collapse.¹² Introducing such a method to solve the wave equation can bring these advantages to this field. Among all Lagrangian meshfree methods, the smoothed particle hydrodynamics (SPH) method is one of the earliest and widely applied methods in different fields. The SPH method was first pioneered independently by Lucy¹³ and Gingold and Monaghan¹⁴ to solve astrophysical problems in 1977. Details about the SPH computation can be found in the literature.^15–19 The SPH and its modified methods were also used to solve different wave equations.^20–22 However, it is still a challenge to find a particle-based computational method for solving wave equations with high-order accuracy.

Considering the generalized finite difference (GFD) scheme for spatial derivatives with arbitrarily distributed points, we proposed a Lagrangian meshfree finite difference particle method (FDPM) for solving wave equations. Meshfree GFD approximation was first discussed for fully arbitrary meshes by Jensen²³ in 1972. Perrone and Kao²⁴ also contributed to the development of this method at that time. Subsequently, a variation using the moving least squares method was proposed by Liszka and Orkisz.²⁵ The meshless finite difference approximation or the GFD approximation that we used came from Benito et al.,²⁶ who gave a discussion about the influence of several factors in the GFD approximation, and a comparison between the GFD scheme and the element-free Galerkin method in solving Laplace equation is also presented in Gavete et al.²⁷ and Benito et al.²⁸ The GFD scheme is more accurate than the Galerkin method. So the FDPM is supposed to be a Lagrangian particle-based approach developed from the high-accuracy GFD scheme. It should be noticed that there are some other particle-based methods improving the accuracy with the Taylor series expansion. Chen et al.^29,30 proposed the corrective smoothed particle method (CSPM) by introducing the Taylor series expansion to normalize the kernel function. After that, Zhang and Batra^31,32 have proposed the modified smoothed particle method (MSPH) and the symmetrical smoothed particle hydrodynamics (SSPH) method with a higher order accuracy. Although these methods used the Taylor series expansion to modify the kernel function, the FDPM still have its advantages.

The FDPM has several advantages compared with the conventional SPH. First, the method is kernel gradient free.³³ Only the kernel function itself and the positions of each particle are used to compute the spatial derivative. Second, the method can be easily extended to higher orders. We propose the second- and fourth-order schemes for the FDPM, and describe how to obtain higher order schemes. It can be seen that only few lines of code need to be changed to obtain a high-order FDPM, which is simple for users. Third, the second derivative can be obtained without special limitations on the kernel function. And so the viscous term in the Burgers’ equation can be obtained directly without introducing any artificial parameters.

Although the FDPM based on the Taylor series expansion shows more accuracy than the traditional SPH method, the motion of the particles still affects the accuracy a lot. As the time goes on, the particles change from regular distribution to irregular, and the support domain for some particles becomes too small to have sufficient neighbor particles or too large to cost unnecessary computing resources. The same problem can also be found in the application of the SPH method and some related work can be found in the literature.^34–36 Recently, Price and Monaghan³⁷ considered a smoothing length in the application of the smoothed particle magnetohydrodynamics. Qiang and Gao³⁸ implicitly coupled the density evolution equation with variable smoothing length to deal with the large density gradient and large smoothing length gradient problems. All these ways to change the smoothing length are proposed with their own pros and cons, and they need to use the physical density to update the length, which is not suitable for this work. So we show another way to update the smoothing length only based on particle positions and propose different approaches to preserve the symmetry of particle interactions.

This article is organized as follows. In section “FDPM,” the FDPM is proposed. Section “Review of SPH method” provides a review of the SPH method. Sections “Computation of the simple wave equation in Lagrangian form” and “Computation of the nonlinear Burgers’ equation in Lagrangian form” evaluate the accuracy and efficiency of the proposed method in solving the Lagrangian simple wave equation and the nonlinear Burgers’ equation. Section “Computation of the acoustic wave equation in Lagrangian form” summarizes the results of this work.

FDPM

GFD scheme for solving spatial derivative

In the FDPM, the GFD scheme is utilized to solve the spatial derivative in governing equations. Consider a particle i surrounded by particles j = 1, 2, …, N, with all N + 1 of them being in a compact support domain. The values of an infinitely differentiable function F at the positions of the particles i and j are defined as F_i and F_j. Let us expand in the Taylor series of F at the particle i

F_{j} = F_{i} + h_{j} \frac{\partial F_{i}}{\partial x} + \frac{h_{j}^{2}}{2} \frac{\partial^{2} F_{i}}{\partial x^{2}} + \frac{h_{j}^{3}}{6} \frac{\partial^{3} F_{i}}{\partial x^{3}} + \frac{h_{j}^{4}}{24} \frac{\partial^{4} F_{i}}{\partial x^{4}} + \cdot \cdot \cdot

(1)

\begin{matrix} F_{j} = F_{i} + h_{j} \frac{\partial F_{i}}{\partial x} + k_{j} \frac{\partial F_{i}}{\partial y} + \frac{h_{j}^{2}}{2} \frac{\partial^{2} F_{i}}{\partial x^{2}} + \frac{k_{j}^{2}}{2} \frac{\partial^{2} F_{i}}{\partial y^{2}} + h_{j} k_{j} \frac{\partial^{2} F_{i}}{\partial x \partial y} \\ + \frac{h_{j}^{3}}{6} \frac{\partial^{3} F_{i}}{\partial x^{3}} + \frac{k_{j}^{3}}{6} \frac{\partial^{3} F_{i}}{\partial y^{3}} + \frac{h_{j}^{2} k_{j}}{2} \frac{\partial^{3} F_{i}}{\partial x^{2} \partial y} + \frac{h_{j} k_{j}^{2}}{2} \frac{\partial^{3} F_{i}}{\partial x \partial y^{2}} + \frac{h_{j}^{4}}{24} \frac{\partial^{4} F_{i}}{\partial x^{4}} \\ + \frac{k_{j}^{4}}{24} \frac{\partial^{4} F_{i}}{\partial y^{4}} + \frac{h_{j}^{3} k_{j}}{6} \frac{\partial^{4} F_{i}}{\partial x^{3} \partial y} + \frac{h_{j}^{2} k_{j}^{2}}{4} \frac{\partial^{4} F_{i}}{\partial x^{2} \partial y^{2}} + \frac{h_{j} k_{j}^{3}}{6} \frac{\partial^{4} F_{i}}{\partial x \partial y^{3}} + \cdot \cdot \cdot \end{matrix}

(2)

where equations (1) and (2) are for one and two dimensions, respectively. x and y are the space coordinates with h_j = x_j – x_i and k_j = y_j – y_i.

Ignoring high-order terms, the approximation of F is denoted by f

f_{j} = f_{i} + h_{j} \frac{\partial f_{i}}{\partial x} + \frac{h_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial x^{2}} + \frac{h_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial x^{3}} + \frac{h_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial x^{4}}

(3)

\begin{matrix} f_{j} = f_{i} + h_{j} \frac{\partial f_{i}}{\partial x} + k_{j} \frac{\partial f_{i}}{\partial y} + \frac{h_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial x^{2}} + \frac{k_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial y^{2}} + h_{j} k_{j} \frac{\partial^{2} f_{i}}{\partial x \partial y} \\ + \frac{h_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial x^{3}} + \frac{k_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial y^{3}} + \frac{h_{j}^{2} k_{j}}{2} \frac{\partial^{3} f_{i}}{\partial x^{2} \partial y} + \frac{h_{j} k_{j}^{2}}{2} \frac{\partial^{3} f_{i}}{\partial x \partial y^{2}} + \frac{h_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial x^{4}} \\ + \frac{k_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial y^{4}} + \frac{h_{j}^{3} k_{j}}{6} \frac{\partial^{4} f_{i}}{\partial x^{3} \partial y} + \frac{h_{j}^{2} k_{j}^{2}}{4} \frac{\partial^{4} f_{i}}{\partial x^{2} \partial y^{2}} + \frac{h_{j} k_{j}^{3}}{6} \frac{\partial^{4} f_{i}}{\partial x \partial y^{3}} \end{matrix}

(4)

After rearranging these equations, the sum of these expressions for all particles j with multiplying the weight functions is obtained

\begin{matrix} \sum_{j = 1}^{N} (f_{i} - f_{j} + h_{j} \frac{\partial f_{i}}{\partial x} + \frac{h_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial x^{2}} + \frac{h_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial x^{3}} + \frac{h_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial x^{4}}) \\ W (h_{j}, r_{s}) = 0 \end{matrix}

(5)

\begin{matrix} \sum_{j = 1}^{N} (\begin{matrix} f_{i} - f_{j} + h_{j} \frac{\partial f_{i}}{\partial x} + k_{j} \frac{\partial f_{i}}{\partial y} + \frac{h_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial x^{2}} + \frac{k_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial y^{2}} + h_{j} k_{j} \frac{\partial^{2} f_{i}}{\partial x \partial y} \\ + \frac{h_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial x^{3}} + \frac{k_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial y^{3}} + \frac{h_{j}^{2} k_{j}}{2} \frac{\partial^{3} f_{i}}{\partial x^{2} \partial y} + \frac{h_{j} k_{j}^{2}}{2} \frac{\partial^{3} f_{i}}{\partial x \partial y^{2}} + \frac{h_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial x^{4}} \\ + \frac{k_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial y^{4}} + \frac{h_{j}^{3} k_{j}}{6} \frac{\partial^{4} f_{i}}{\partial x^{3} \partial y} + \frac{h_{j}^{2} k_{j}^{2}}{4} \frac{\partial^{4} f_{i}}{\partial x^{2} \partial y^{2}} + \frac{h_{j} k_{j}^{3}}{6} \frac{\partial^{4} f_{i}}{\partial x \partial y^{3}} \end{matrix}) \\ W (h_{j}, k_{j}, r_{s}) \end{matrix} = 0

(6)

where W(h_j,r_s) and W(h_j,k_j,r_s) are the weight functions in one and two dimensions, respectively, and r_s represents the size of the support domain, which is always called the smoothing length in the SPH method. For a circular shape domain, r_s stands for the radius of the support domain.

The closest nodes to the particle i are selected as particles j, and these particles should be in the support domain at the same time. In addition, quintic splines are selected in both CSPM and FDPM as the weight function in this work. More details about the choice of weight functions can be found in the literarture.^19,26,33,39 In the following equations, we use W to denote the weight function.

Define the functions G₁ in one dimension and G₂ in two dimensions as

\begin{matrix} G_{1} (f) = \sum_{j = 1}^{N} \\ {[(f_{i} - f_{j} + h_{j} \frac{\partial f_{i}}{\partial x} + \frac{h_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial x^{2}} + \frac{h_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial x^{3}} + \frac{h_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial x^{4}}) W]}^{2} \end{matrix}

(7)

\begin{matrix} G_{2} (f) = \sum_{j = 1}^{N} \\ {[(\begin{matrix} f_{i} - f_{j} + h_{j} \frac{\partial f_{i}}{\partial x} + k_{j} \frac{\partial f_{i}}{\partial y} + \frac{h_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial x^{2}} + \frac{k_{j}^{2}}{2} \frac{\partial^{2} f_{i}}{\partial y^{2}} + h_{j} k_{j} \frac{\partial^{2} f_{i}}{\partial x \partial y} \\ + \frac{h_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial x^{3}} + \frac{k_{j}^{3}}{6} \frac{\partial^{3} f_{i}}{\partial y^{3}} + \frac{h_{j}^{2} k_{j}}{2} \frac{\partial^{3} f_{i}}{\partial x^{2} \partial y} + \frac{h_{j} k_{j}^{2}}{2} \frac{\partial^{3} f_{i}}{\partial x \partial y^{2}} + \frac{h_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial x^{4}} \\ + \frac{k_{j}^{4}}{24} \frac{\partial^{4} f_{i}}{\partial y^{4}} + \frac{h_{j}^{3} k_{j}}{6} \frac{\partial^{4} f_{i}}{\partial x^{3} \partial y} + \frac{h_{j}^{2} k_{j}^{2}}{4} \frac{\partial^{4} f_{i}}{\partial x^{2} \partial y^{2}} + \frac{h_{j} k_{j}^{3}}{6} \frac{\partial^{4} f_{i}}{\partial x \partial y^{3}} \end{matrix}) W]}^{2} \end{matrix}

(8)

According to equations (5) and (6), the norms G₁ and G₂ are equal to 0, so we obtain

\frac{\partial G_{1} (f)}{\partial (\frac{\partial f_{i}}{\partial x})} = 0, \frac{\partial G_{1} (f)}{\partial (\frac{\partial^{2} f_{i}}{\partial x^{2}})} = 0, \frac{\partial G_{1} (f)}{\partial (\frac{\partial^{3} f_{i}}{\partial x^{3}})} = 0, \frac{\partial G_{1} (f)}{\partial (\frac{\partial^{4} f_{i}}{\partial x^{4}})} = 0

(9)

\begin{matrix} \frac{\partial G_{2} (f)}{\partial (\frac{\partial f_{i}}{\partial x})} = 0, \frac{\partial G_{2} (f)}{\partial (\frac{\partial f_{i}}{\partial y})} = 0, \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{2} f_{i}}{\partial x^{2}})} = 0, \\ \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{2} f_{i}}{\partial y^{2}})} = 0, \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{2} f_{i}}{\partial x \partial y})} = 0, \\ \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{3} f_{i}}{\partial x^{3}})} = 0, \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{3} f_{i}}{\partial y^{3}})} = 0, \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{3} f_{i}}{\partial x^{2} \partial y})} = 0, \\ \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{3} f_{i}}{\partial x \partial y^{2}})} = 0, \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{4} f_{i}}{\partial x^{4}})} = 0, \\ \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{4} f_{i}}{\partial y^{4}})} = 0, \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{4} f_{i}}{\partial x^{3} \partial y})} = 0, \\ \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{4} f_{i}}{\partial x^{2} \partial y^{2}})} = 0, \frac{\partial G_{2} (f)}{\partial (\frac{\partial^{4} f_{i}}{\partial x \partial y^{3}})} = 0 \end{matrix}

(10)

From equations (9) and (10), the following equations can be obtained

A_{1} D_{1} = B_{1}

(11)

A_{2} D_{2} = B_{2}

(12)

where

A_{1} = [\begin{matrix} \sum_{j = 1}^{N} h_{j}^{2} W^{2} & \sum_{j = 1}^{N} \frac{1}{2} h_{j}^{3} W^{2} & \sum_{j = 1}^{N} \frac{1}{6} h_{j}^{4} W^{2} & \sum_{j = 1}^{N} \frac{1}{24} h_{j}^{5} W^{2} \\ \sum_{j = 1}^{N} \frac{1}{4} h_{j}^{4} W^{2} & \sum_{j = 1}^{N} \frac{1}{12} h_{j}^{5} W^{2} & \sum_{j = 1}^{N} \frac{1}{48} h_{j}^{6} W^{2} \\ \sum_{j = 1}^{N} \frac{1}{36} h_{j}^{6} W^{2} & \sum_{j = 1}^{N} \frac{1}{144} h_{j}^{7} W^{2} \\ SYM & \sum_{j = 1}^{N} \frac{1}{576} h_{j}^{8} W^{2} \end{matrix}]

(13)

A_{2} = [\begin{matrix} \sum_{j = 1}^{N} h_{j}^{2} W^{2} & \sum_{j = 1}^{N} h_{j} k_{j} W^{2} & \sum_{j = 1}^{N} \frac{1}{2} h_{j}^{3} W^{2} & \sum_{j = 1}^{N} \frac{1}{2} h_{j} k_{j}^{2} W^{2} & \dots & \sum_{j = 1}^{N} \frac{1}{6} h_{j}^{2} k_{j}^{3} W^{2} \\ \sum_{j = 1}^{N} k_{j}^{2} W^{2} & \sum_{j = 1}^{N} \frac{1}{2} h_{j}^{2} k_{j} W^{2} & \sum_{j = 1}^{N} \frac{1}{2} k_{j}^{3} W^{2} & \dots & \sum_{j = 1}^{N} \frac{1}{6} h_{j} k_{j}^{4} W^{2} \\ \sum_{j = 1}^{N} \frac{1}{4} h_{j}^{4} W^{2} & \sum_{j = 1}^{N} \frac{1}{4} h_{j}^{2} k_{j}^{2} W^{2} & \dots & \sum_{j = 1}^{N} \frac{1}{12} h_{j}^{3} k_{j}^{3} W^{2} \\ \sum_{j = 1}^{N} \frac{1}{4} k_{j}^{4} W^{2} & \dots & \sum_{j = 1}^{N} \frac{1}{12} h_{j} k_{j}^{5} W^{2} \\ ⋱ & ⋮ \\ SYM & \sum_{j = 1}^{N} \frac{1}{36} h_{j}^{2} k_{j}^{6} W^{2} \end{matrix}]

(14)

D_{1} = {\frac{\partial f_{i}}{\partial x} \frac{\partial^{2} f_{i}}{\partial x^{2}} \frac{\partial^{3} f_{i}}{\partial x^{3}} \frac{\partial^{4} f_{i}}{\partial x^{4}}}^{T}

(15)

D_{2} = {\frac{\partial f_{i}}{\partial x} \frac{\partial f_{i}}{\partial y} \frac{\partial^{2} f_{i}}{\partial x^{2}} \frac{\partial^{2} f_{i}}{\partial y^{2}} \dots \frac{\partial^{4} f_{i}}{\partial x \partial y^{3}}}^{T}

(16)

B_{1} = {\begin{matrix} \sum_{j = 1}^{N} (f_{j} - f_{i}) h_{j} W^{2} \\ \sum_{j = 1}^{N} (f_{j} - f_{i}) \frac{h_{j}^{2}}{2} W^{2} \\ \sum_{j = 1}^{N} (f_{j} - f_{i}) \frac{h_{j}^{3}}{6} W^{2} \\ \sum_{j = 1}^{N} (f_{j} - f_{i}) \frac{h_{j}^{4}}{24} W^{2} \end{matrix}}

(17)

B_{2} = {\begin{matrix} \sum_{j = 1}^{N} (f_{j} - f_{i}) h_{j} W^{2} \\ \sum_{j = 1}^{N} (f_{j} - f_{i}) k_{j} W^{2} \\ \sum_{j = 1}^{N} (f_{j} - f_{i}) \frac{h_{j}^{2}}{2} W^{2} \\ \sum_{j = 1}^{N} (f_{j} - f_{i}) \frac{k_{j}^{2}}{2} W^{2} \\ ⋮ \\ \sum_{j = 1}^{N} (f_{j} - f_{i}) \frac{h_{j}^{2} k_{j}^{3}}{6} W^{2} \end{matrix}}

(18)

Since the matrices A₁ and A₂ are symmetrical, equations (11) and (12) can be solved using the Cholesky factorization. In this way, the spatial derivative can be obtained.

Time integration scheme

In our Lagrangian approach, calculation of the particle motion and time integration are performed with second-order leapfrog integration. In this scheme, the equations for updating particle positions and velocity are

v_{i} (t + \frac{1}{2} Δ t) = v_{i} (t - \frac{1}{2} Δ t) + Δ t \frac{D_{0} v_{i} (t)}{Dt}

(19)

r_{i} (t + Δ t) = r_{i} (t) + Δ t v_{i} (t + \frac{1}{2} Δ t)

(20)

where v _i is the velocity of the fluid particle i, t is the time of the previous time step, and Δt is the time step.

Equation (19) starts with the initial velocity offset given by

v_{i} (\frac{1}{2} Δ t) = v_{i} (0) + \frac{1}{2} Δ t \frac{D_{0} v_{i} (0)}{Dt}

(21)

Variable smoothing length

In the FDPM, the computation of the weight function depends on the support domain. The support domain around a particle decides its neighbor particles, so the shape and size of the support domain are important. The most used domain is the circular (two-dimensional (2D)) or spherical (three-dimensional (3D)) domain because of its feasible and stable property. A circular support domain in the computational domain Ω is shown in Figure 1. It is obvious that the circular radius decides the size of the support domain, so this radius of the support domain (r_s) is used in the computation and named as the smoothing length or the kernel length. An initial smoothing length (r_s0) is needed before the computation.

Figure 1.

Support domain (dot dash circle) and the smoothing length (r_s) for particle i in the computational domain (Ω).

Here we proposed a simple approach to dynamically update r_s so that the number of the neighbor particles remains relatively constant as

r_{s} (i) = r_{s 0} (i) {(\frac{\sum_{j = 1}^{N} | [x (j) - x (i)] W |}{\sum_{j = 1}^{N} | [x_{0} (j) - x_{0} (i)] W |})}^{\frac{1}{d}}

(22)

where the particle j is the neighbor particle of the particle i, N is the total number of neighbor particles of the particle i, x is the position of the particle, 0 means the value at the initial time, and d is the number of dimensions.

Since r_s is set to vary in both time and space, each particle has its own r_s. If r_s(i) is not equal to r_s(j), the support domain of the particle i may cover the particle j but the particle j does not cover the particle i. In order to overcome this problem, several methods can be used to obtain the symmetric r_s(i,j) as follows

r_{s} (i, j) = min (r_{s} (i), r_{s} (j))

(23)

r_{s} (i, j) = max (r_{s} (i), r_{s} (j))

(24)

r_{s} (i, j) = \frac{r_{s} (i) + r_{s} (j)}{2}

(25)

r_{s} (i, j) = \frac{2 r_{s} (i) r_{s} (j)}{r_{s} (i) + r_{s} (j)}

(26)

In this work, we use the third one to get the symmetric r_s(i,j).

Review of SPH method

The integral representation for f(r) used in the SPH method can be written as

< f (r) > = \int_{Ω} f (r') W (r - r', r_{s}) d r'

(27)

where f is a function of the position vector r , W is the smoothing kernel function, and r_s denotes the influence area of the smoothing function.

The kernel approximation of f(r) and its derivative used in the SPH method can be described as the summation of neighboring particles as

< f (r) > = \sum_{j = 1}^{N} \frac{m_{j}}{ρ_{j}} f (r_{j}) W_{ij}

(28)

< \nabla \cdot f (r) > = - \sum_{j = 1}^{N} \frac{m_{j}}{ρ_{j}} f (r_{j}) \cdot \nabla_{i} W_{ij}

(29)

where N indicates the number of particles in the support domain, m_j is the mass of the particle j, and ρ_j is the density of the particle j. In the SPH convention, the kernel approximation operator is marked by the angle brackets <>.

Computation of the simple wave equation in Lagrangian form

Simple wave equation in Lagrangian form

The simple wave equations⁴⁰ with only the first-order partial differential terms are always used as fundamental tests in acoustic and hydrodynamic wave problems. This equation can also be seen as the advection equation. To evaluate the FDPM algorithm, we consider the numerical solution of the one-dimensional (1D) simple wave equation in Lagrangian form

\frac{d u}{d t} = (u - 1) \frac{\partial u}{\partial x}

(30)

where u is the velocity, t represents the time, x is the coordinate, and the Lagrangian derivative is

\frac{d}{d t} = \frac{\partial}{\partial t} + u \frac{\partial}{\partial x}

(31)

with the initial condition defined over −20 m ≤ x ≤ 60 m as

u (x, 0) = 0.5 \exp [- \ln 2 {(\frac{x}{3})}^{2}]

(32)

The boundary condition at both ends is set as u = 0 as a Dirichlet boundary condition.

The exact solution to the problem is available as

u (x, t) = 0.5 \exp [- \ln 2 {(\frac{x - t}{3})}^{2}]

(33)

Nondimensional error definition

The numerical accuracy is evaluated by the root mean square error and the maximum error, which are defined as

ε_{RMS} (S) = \frac{\sqrt{\frac{1}{NT} \sum_{k = 1}^{NT} {| S_{num} (k) - S_{ana} (k) |}^{2}}}{\sqrt{\frac{1}{NT} \sum_{k = 1}^{NT} {| S_{ana} (k) |}^{2}}}

(34)

ε_{MAX} (S) = \frac{max_{1 \leq k \leq NT} | S_{num} (k) - S_{ana} (k) |}{max_{1 \leq k \leq NT} | S_{ana} (k) |}

(35)

where NT is the number of all particles in the computational domain.

In this section, the FDPM is based on the second-order Taylor series expansion.

Validation of the computation

Numerical results for FDPM with variable smoothing length in solving the Lagrangian simple wave equation are shown in Figure 2. The initial particle spacing is set as Δx = 0.15 m, the time step Δt = 0.03 s, and r_s = 3.5Δx. At the start, the particles are distributed regularly along the computational domain.

Figure 2.

Comparison between FDPM (second order with variable r_s) and exact solutions in solving the simple wave equation in Lagrangian form: (a) t = 5 s and (b) t = 40 s.

The wave propagates from the left to right in the figure, and a peak appears during the propagation. Comparing with the exact solutions, it can be seen that the FDPM with variable smoothing length accurately simulates the waves, and the peak value is obtained precisely.

Accuracy and efficiency evaluation

Different wave propagation models with time changing from 2 to 40 s are simulated with CSPM and FDPM with or without variable smoothing length. The root mean square and maximum error of the velocity are shown in Figure 3.

Figure 3.

Comparison between CSPM and FDPM (second order with constant or variable r_s) solutions of the simple wave equation in Lagrangian form at different times: (a) ε_RMS(u) and (b) ε_MAX(u).

At the beginning, the CSPM error is obviously larger than the FDPM error, and little difference can be found between the FDPM with constant and variable r_s. After that, the errors of three methods increase along with the increase of time. Numerical results with the root mean square and maximum error from maximum to minimum are the CSPM, FDPM with constant r_s, and FDPM with variable r_s, respectively.

The effect of the initial smoothing length r_s0 is shown in Figure 4. r_s₀/Δx for CSPM and FDPM with constant r_s ranges from 2.5 to 8 which ensures at least four neighbor particles at the beginning. When r_s0/Δx is smaller than 2.5, both methods cannot obtain the final results. For the FDPM with variable smoothing length, r_s0/Δx should be larger than 2.5.

Figure 4.

Comparison between CSPM and FDPM (second order with constant or variable r_s) solutions of the simple wave equation in Lagrangian form with different initial smoothing length r_s0: (a): ε_RMS(u) and (b) ε_MAX(u).

As shown in Figure 4, both root mean square and maximum numerical error increase along with the increase of r_s0/Δx. We can see from the line graph that CSPM solution error is much larger than the error of FDPM. The FDPM with using variable r_s is the most accurate one.

Different computational cases corresponding to the CSPM and FDPM results with similar numerical error are shown in Table 1. The purpose of the computation is to obtain the CPU time cost of the CSPM and FDPM with variable r_s when the numerical error is in the same level. The performance of computation is measured using the Intel Core i3-6100 processor (3.7 GHz) with 4 GB RAM. It can be seen from the data that the CSPM needs smaller particle spacing to reduce numerical error compared with the FDPM with variable smoothing length.

Table 1.

Efficiency comparison between CSPM and FDPM with variable r_s in solving the simple wave equation.

	r _s	x	Δx	ε _RMS	ε _MAX	CPU time (s)
CSPM	Constant	[−20,200]	0.07	0.0197	0.0228	1515
FDPM	Variable	[−20,200]	0.15	0.0201	0.0189	1445

CSPM: corrective smoothed particle method; FDPM: finite difference particle method.

Computation of the nonlinear Burgers’ equation in Lagrangian form

Burgers’ equation in Lagrangian form

Because of having the nonlinear convection terms and the occurrence of viscosity, Burgers’ equations are a simple form of the nonlinear partial differential equations. Hence, they have often been used as a test model for ensuring the accuracy and stability of numerical methods. For this reason, Burgers’ equations are considered in this work for testing the applicability of the FDPM to nonlinear dynamic problems.

The Lagrangian form of 1D Burgers’ equation is

\frac{d u}{d t} = υ \frac{\partial^{2} u}{\partial x^{2}}

(36)

where u is the velocity, t represents the time, x is the coordinate, and ν is the kinematic viscosity.

Considering the time-dependent Burgers’ equation, the initial condition for solving the Burgers’ equation is defined over −0.5 m ≤ x ≤ 1.5 m as

u (x) = \sin (π x)

(37)

The exact Fourier solution⁴¹ to the problem is available as

u (x, t) = 2 π υ \frac{\sum_{n = 1}^{\infty} c_{n} \exp (- n^{2} π^{2} υ t) n \sin (n π x)}{c_{0} + \sum_{n = 1}^{\infty} c_{n} \exp (- n^{2} π^{2} υ t) \cos (n π x)}

(38)

where the coefficients are defined as

c_{0} = \int_{0}^{1} \exp (- \frac{1 - \cos (π x)}{2 π υ}) d x

(39)

c_{n} = 2 \int_{0}^{1} \exp (- \frac{1 - \cos (π x)}{2 π υ}) \cos (n π x) d x

(40)

The number n is chosen such that the coefficient c_n is less than 1.0e–15.

Validation of the computation

Considering the nonlinear Burgers’ equation, Figure 5 shows the comparison of the FDPM results with exact solutions. The initial particle spacing is set as Δx = 0.025 m, the time step Δt = 0.0001 s, r_s = 5.8Δx, and υ = 100 m²/s. At the beginning, the particles are distributed regularly along the computational domain. The FDPM with second-order Taylor series expansion is used.

Figure 5.

Comparison between FDPM (second order with variable r_s) and exact solutions in solving the Burgers’ equation in Lagrangian form.

The figure shows that the FDPM with variable r_s solves the Burgers’ equation with good results compared with the exact solutions. It should be noticed that, since the particle moves with time, the particle’s position changes from regular to irregular along with the computation, which can also be seen from the figure. With irregular particle distribution, the variable smoothing length can work well.

Accuracy and efficiency evaluation

The effect of the time changing from 0.15 to 0.30 s is shown in Figure 6.

Figure 6.

Comparison between CSPM and FDPM (second order with constant or variable r_s) solutions of the Burgers’ equation in Lagrangian form at different times: (a) ε_RMS(u) and (b) ε_MAX(u).

As shown in Figure 6, both root mean square and maximum numerical error increase along with the increase of t. We can see from the line graph that the CSPM solution error is much larger than the error of FDPM. The FDPM with using variable r_s is the most accurate one. When t = 0.3 s, the ε_RMS values from the CSPM, FDPM with constant r_s, and FDPM with variable r_s solutions reach about 0.04, 0.012, and 0.001, respectively.

Different wave propagation models with r_s0/Δx changing from 2 to 8 are simulated with CSPM and FDPM with or without variable smoothing length. The root mean square and maximum error of the velocity are shown in Figure 7.

Figure 7.

Comparison between CSPM and FDPM (second order with constant or variable r_s) solutions of the Burgers’ equation in Lagrangian form with different initial smoothing lengths r_s0: (a) ε_RMS(u) and (b) ε_MAX(u).

When r_s0/Δx is small, the CSPM error is obviously larger than the FDPM, and difference can be found between the FDPM with constant and the FDPM with variable r_s. When r_s0/Δx becomes larger, the errors of the three methods increase along with the r_s0/Δx increase, both the root mean square and maximum error from maximum to minimum correspond to the CSPM, FDPM with constant r_s, and FDPM with variable r_s, respectively.

For the nonlinear Burgers’ equation, different computational cases corresponding to the CSPM and FDPM results with similar numerical error are shown in Table 2. It can be seen from the data that the CSPM needs smaller particle spacing to reduce numerical error compared with the FDPM with variable smoothing length. With numerical error in the same level, the FDPM with variable r_s consumes 196 s to compute, while the CSPM costs 208 s. The FDPM with variable r_s is more efficient than the CSPM for this case.

Table 2.

Efficiency comparison between CSPM and FDPM with variable r_s in solving the Burgers’ equation.

	r _s	x	Δx	ε _RMS	ε _MAX	CPU time (s)
CSPM	Constant	[−0.5,21.5]	0.01	0.0047	0.0223	208
FDPM	Variable	[−0.5,21.5]	0.05	0.0045	0.0166	196

CSPM: corrective smoothed particle method; FDPM: finite difference particle method.

Computation of the acoustic wave equation in Lagrangian form

Acoustic wave equation in Lagrangian form

In order to test the performance of the proposed method in solving the 2D problem, sound propagation is simulated using the FDPM. Governing equations for acoustic waves in Lagrangian form can be found in Zhang et al.,²¹ which are defined as

\frac{d δ ρ}{d t} = - ρ_{0} \nabla \cdot v

(41)

\frac{d v}{d t} = - \frac{1}{ρ_{0}} \nabla δ ρ

(42)

\frac{d δ p}{d t} = c_{0}^{2} \frac{d ρ'}{d t}

(43)

where ρ′ is the density change, t is the time, ρ₀ is the density, v is the velocity, p is the sound pressure, and c₀ is the sound speed.

A Gaussian pulse model of sound propagation is built. The center of the computation region is located at the origin of the coordinates, the computation domain is a square, and the length of the square is 60 m. At the time t = 0, the distribution of the sound pressure of the Gaussian pulse model in the computation domain can be written as

δ p (t = 0, x, y) = α_{1} ρ_{0} c_{0}^{2} \exp [- α_{2} (x^{2} + y^{2})]

(44)

where α₁ = 15/c₀², α₂ = ln2/16, and c₀ = 340 m/s.

Acoustic wave

The simulation results of the Gaussian pulse of sound propagation are shown in Figure 8. Figure 8(a) and (b) shows the contour of sound pressure at t = 0.01 and 0.03 s, respectively. It can be seen that the acoustic wave propagates all around in the form of a circle, and the center of the circle is located at the origin of the coordinate. As the time goes on, the radius of the circle increases, while the amplitude of the sound waves decreases.

Figure 8.

Sound pressure contour at different times obtained by FDPM: (a) t = 0.01 s and (b): t = 0.03 s.

For providing a better comparison, the sound pressure of the particles at y = 0 m are obtained and compared with the theoretical solution, as shown in Figure 9. In the figure, the red solid line represents the theoretical solution and the black scattered points represent the simulation results.

Figure 9.

Comparison between the FDPM (second order with variable r_s) and theoretical results.

From the figure, it can be seen that, at y = 0 m, the FDPM simulation results possess similar smoothing trends compared to the theoretical solution. With the time increasing, the amplitude decreases. Generally, the FDPM algorithms and codes can simulate the Gaussian pulse of sound propagation well.

Conclusion

A Lagrangian meshfree FDPM is proposed to numerically solve different wave equations. A variable smoothing length approach for updating the support domain size is suggested. The numerical method is tested with the simple wave equation and the Burgers’ equation in Lagrangian form. Comparison with the CSPM, FDPM with constant smoothing length, and FDPM with variable smoothing length is made. In summary, we have shown the following:

The proposed FDPM solves the simple wave equation and the Burgers’ equation with good accuracy;

An updating smoothing length approach based on particle distributions is given. It is found to be able to reduce numerical error when the particles are distributed irregularly;

Compared with the CSPM, the FDPM with variable smoothing length is more accurate with different initial particle spacing, and so it needs fewer particles to maintain the same level of numerical error which costs less CPU time in solving wave equations.

Footnotes

Handling Editor: Kai Bao

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was supported by the National Natural Science Foundation of China (No. 51741910), and the Fundamental Research Funds for the Central Universities (WUT: 2017IVB081).

ORCID iD

Yong Ou Zhang

References

Chai

Gong

et al . Hybrid smoothed finite element method for two-dimensional underwater acoustic scattering problems. Ocean Eng 2016; 116: 129–141.

Sanchez-Orgaz

Denia

Martinez-Casas

et al . 3D acoustic modelling of dissipative silencers with nonhomogeneous properties and mean flow. Adv Mech Eng 2014; 6: 537935.

Chai

Lei

et al . Analysis of coupled structural-acoustic problems based on the smoothed finite element method (S-FEM). Eng Anal Bound Elem 2014; 42: 84–91.

Durran

. Numerical methods for wave equations in geophysical fluid dynamics. New York: Springer Press, 2013.

Zhang

Cheng

et al . The improved element-free Galerkin method for three-dimensional wave equation. Acta Mech Sinica 2012; 28: 808–818.

Dehghan

Shokri

. A meshless method for numerical solution of the one-dimensional wave equation with an integral condition using radial basis functions. Numer Algorithms 2009; 52: 461–477.

Chen

et al . Singular boundary method: three regularization approaches and exterior wave applications. Comp Model Eng 2014; 99: 417–443.

Shivanian

. Meshless local Petrov-Galerkin (MLPG) method for three-dimensional nonlinear wave equations via moving least squares approximation. Eng Anal Bound Elem 2015; 50: 249–257.

Hon

Schaback

. Direct meshless kernel techniques for time-dependent equations. Appl Math Comput 2015; 258: 220–226.

10.

Yang

Liu

Peng

et al . Numerical modeling of dam-break flow impacting on flexible structures using an improved SPH–EBG method. Coastal Eng 2016; 108: 56–64.

11.

Ihmsen

Cornelis

Solenthaler

et al . Implicit incompressible SPH. IEEE T Vis Comput Gr 2014; 20: 426–435.

12.

Tian

Qiu

et al . Numerical solution on spherical vacuum bubble collapse using MPS method. J Eng Gas Turbines Power 2010; 132: 102920.

13.

Lucy

. A numerical approach to the testing of the fission hypothesis. Astron J 1977; 82: 1013–1024.

14.

Gingold

Monaghan

. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Month Not R Astron Soc 1977; 181: 375–389.

15.

Liu

. Smoothed particle hydrodynamics (SPH): an overview and recent developments. Arch Comput Method E 2010; 17: 25–76.

16.

Liu

Zong

. An overview on smoothed particle hydrodynamics. Int J Comput Methods 2008; 5: 135–188.

17.

Monaghan

. Smoothed particle hydrodynamics and its diverse applications. Annu Rev Fluid Mech 2012; 44: 323–346.

18.

Springel

. Smoothed particle hydrodynamics in astrophysics. Annu Rev Astron Astr 2010; 48: 391–430.

19.

Liu

Smoothed particle hydrodynamics: a meshfree particle method. Singapore: World Scientific, 2003.

20.

Zhang

Ouyang

et al . SPH simulation of sound propagation and interference. In: Proceedings of the 5th international conference on computational methods, Cambridge, 28–30 July 2014. UK.

21.

Zhang

Ouyang

et al . Efficient SPH simulation of time-domain acoustic propagation. Eng Anal Bound Elem 2016; 62: 112–122.

22.

Zhang

Ouyang

et al . SPH simulation of acoustic waves: effects of frequency, sound pressure, and particle spacing. Math Probl Eng 2015; 2015: 348314.

23.

Jensen

. Finite difference techniques for variable grids. Comput Struct 1972; 2: 17–29.

24.

Perrone

Kao

. A general finite difference method for arbitrary meshes. Comput Struct 1975; 5: 45–57.

25.

Liszka

Orkisz

. The finite difference method at arbitrary irregular grids and its application in applied mechanics. Comput Struct 1980; 11: 83–95.

26.

Benito

Urena

Gavete

. Influence of several factors in the generalized finite difference method. Appl Math Model 2001; 25: 1039–1053.

27.

Gavete

Benito

. Improvements of generalized finite difference method and comparison with other meshless method. Appl Math Model 2003; 27: 831–847.

28.

Benito

Urena

Gavete

. Solving parabolic and hyperbolic equations by the generalized finite difference method. J Comput Appl Math 2007; 209: 208–233.

29.

Chen

Beraun

Carney

. A corrective smoothed particle method for boundary value problems in heat conduction. Int J Numer Meth Eng 1999; 46: 231–252.

30.

Chen

Beraun

Jih

. Completeness of corrective smoothed particle method for linear elastodynamics. Comput Mech 1999; 24: 273–285.

31.

Zhang

Batra

. Modified smoothed particle hydrodynamics method and its application to transient problems. Comput Mech 2004; 34: 137–146.

32.

Batra

Zhang

. SSPH basis functions for meshless methods, and comparison of solutions with strong and weak formulations. Comput Mech 2008; 41: 527–545.

33.

Huang

Lei

Liu

et al . A kernel gradient free (KGF) SPH method. Int J Numer Meth Fl 2015; 78: 691–707.

34.

Hernquist

Katz

. TREESPH: a unification of SPH with the hierarchical tree method. Astron J Suppl Ser 1989; 70: 419–446.

35.

Springel

Hernquist

. Cosmological smoothed particle hydrodynamics simulations: the entropy equation. Month Not R Astron Soc 2002; 333: 649–664.

36.

Monaghan

. SPH compressible turbulence. Month Not R Astron Soc 2002; 335: 843–869.

37.

Price

Monaghan

. Smoothed particle magnetohydrodynamics II. Variational principles and variable smoothing-length terms. Month Not R Astron Soc 2004; 348: 139–152.

38.

Qiang

Gao

. A new SPH equation including variable smoothing lengths aspects and its implementation. In: Proceedings of the international symposium on computational mechanics, Beijing, China, 30 July–1 August 2007.

39.

David

Dennis

. An analysis of 1-D smoothed particle hydrodynamics kernels. J Comput Phys 1996; 126: 165–180.

40.

Gurbatov

Rudenko

Saichev

. Waves and structures in nonlinear nondispersive media. China: Higher Education Press, 2011.

41.

Liao

. An implicit fourth-order compact finite difference scheme for one-dimensional Burgers’ equation. Appl Math Comput 2008; 206: 755–764.