A SPH-GFDM Coupled Method for Elasticity Analysis

Abstract

SPH (smoothed particle hydrodynamics) is one of the oldest meshless methods used to simulate mechanics of continuum media. Despite its great advantage over the traditional grid-based method, implementing boundary conditions in SPH is not easy and the accuracy near the boundary is low. When SPH is applied to problems for elasticity, the displacement or stress boundary conditions should be suitably handled in order to achieve fast convergence and acceptable numerical accuracy. The GFDM (generalized finite difference method) can derive explicit formulae for required partial derivatives of field variables. Hence, a SPH–GFDM coupled method is developed to overcome the disadvantage in SPH. This coupled method is applied to 2-D elastic analysis in both symmetric and asymmetric computational domains. The accuracy of this method is demonstrated by the excellent agreement with the results obtained from FEM (finite element method) regardless of the symmetry of the computational domain. When the computational domain is multiply connected, this method needs to be further improved.

1. Introduction

Smoothed particle hydrodynamics (SPH) is a Lagrangian and mesh-free technique for computational modeling of continuum systems such as solids and fluids. The SPH method, which was originally introduced by Gingold and Monaghan [1] and Lucy [2] in astrophysics, has been extended to many research fields including fluid dynamics, heat transfer, and high velocity impact problems in computational domains of various kinds.

With the increasing application of SPH to more and more engineering problems, many challenges have emerged in dealing consistently with boundary conditions. These difficulties are mainly caused by two aspects. Firstly, the shape of the boundary of the computational domain is usually complicated in practical engineering problems. Secondly, for those particles near the boundary, their support domain will be truncated. Due to these difficulties, the accuracy of the SPH approximation is reduced dramatically near the boundary and then these will have an influence on the overall accuracy of the simulation. In view of this, proper boundary treatments have been an ongoing concern for the successful implementation of the SPH method.

In the last few years, many approaches have been proposed to improve the boundary condition implementation in SPH method for fluid problems. In the existing literature, there are mainly two kinds of methods to treat the boundary conditions (namely, virtual particle methods and semi-analytical methods). The virtual particle methods consist of boundary particle force methods, the image particle method, and the dynamic particle method. The main idea of the boundary particle force method [3] is that some specified boundary particles are placed on the boundary and these particles exert repulsive forces to the inner fluid particles to prevent the fluid particles from penetrating the wall. This method is easy to implement, but unfortunately, it contains some empirical parameters which cannot be determined in advance.

In the image particle method [4], image particles are placed outside of the fluid region. These particles have the same pressure, density, and temperature as the inner fluid particles. When there is an external force acting on the fluid particles, the pressure of the image particles must be adjusted. These image particles have the same perpendicular component of the corresponding fluid particle velocity but the opposite tangential component. The image particles are generated through reflecting or mirroring the corresponding fluid particles along the wall in each time step. The disadvantage of this model is that when the boundary changes with time, it is difficult to determine how the image particle can be quickly generated and placed.

The dynamic particle method was firstly introduced by Liu [5]. In this model, the dynamic particles are introduced and fixed outside the solid boundary to approximate the field variables of the fluid particle. Hence, the accuracy is improved to some extent. In the meantime, the field variables of the dynamic particle are acquired from the governing equation of the fluid particles, such as the state equation and continuity equation. It should be pointed out that though the numerical accuracy is improved, the boundary deficiency of the dynamic particles exists, which may lead to pressure oscillations.

The main idea of the semi-analytical method is based on a so-called renormalizing factor to reestablish the governing equations for the fluid particle near the boundary. This factor depends on the location of the fluid particle relative to the fixed wall and the shape of the boundary. Kulasegaram et al. [6] once proposed a new formula to determine the renormalizing factor of the missing area of the kernel support for the particles near the boundary. This model was later improved by Ferrand [7] through solving an extra dynamic equation for the renormalizing factor. This method is more accurate but too troublesome.

For the mixed finite element method, the displacement field and stress field inside the element are represented by a node displacement vector and internal force vector. Then, the mixed model can be obtained by applying the generalized variational principle. A mixed-type finite element approximation for radiation problems using fictitious domain method written in the form of pseudo-differential operator was proposed by Nasir et al. [8], which is also efficiently applied to compute the solution of the radiation problem outside the computational domain and to compute the far-field pattern. The mixed finite element method utilizes the relatively simple interpolation function, which are also widely applied in thermodynamics [9,10], but the coefficient matrix of the simultaneous equations it derived is not positive definite, which limits the wide application of this method to a certain extent.

When SPH is applied to heat transfer problems and solid mechanics, the support domain truncation problem become harder. This is due to the fact that in these situations, the Neumann and Robin boundary conditions are almost unavoidable. Consequently, the literature on SPH application in solid mechanics is sparse.

M.A. Esmaili Sikarudi et al. [11] proposed smoothing directional derivatives method and manipulated Taylor series method when SPH was applied to heat transfer. Wenxiao Pan et al. [12] proposed a continuous boundary force (CBF) method for SPH to solve the Navier–Stokes equations subjected to the Robin boundary condition. In the CBF method, a volumetric force term is added to the momentum equation and the Robin boundary condition is replaced by the homogeneous Neumann boundary condition. Härdi, Simon, et al. [13,14] enhanced SPH by using the finite particle method to numerically solve the shallow water equations and the results agreed well with the analytical solution.

For elasticity problems, the Dirichlet boundary condition and Neumann boundary condition are commonly encountered due to the displacement and stress boundary conditions. When SPH is applied to elasticity problems, the boundary conditions must be carefully treated.

To ensure the accuracy of the numerical solution, the derivatives of the field function must be approximated as accurate as possible in region near the boundary of the computational domain. In this paper, we developed a SPH-GFDM coupled method and apply it to the analysis of 2D linear elastic problems. The numerical results are compared with commercial FEM software COMSOL to verify its accuracy. This method takes in both advantages of SPH and GFDM, which makes it a relatively more robust meshless method for solving PDEs.

2. Methodology

2.1. Basic Ideas and Differential Operators

The main idea of SPH is based on the integral interpolation. The integral representation of a field function $f\left(\mathit{x}\right)$ can be written as:

$$f(\mathit{x})={\displaystyle \int f}\left({\mathit{x}}^{\prime }\right)\delta \left(\mathit{x}-{\mathit{x}}^{\prime }\right)\mathrm{d}{\mathit{x}}^{\prime }$$

(1)

where $\mathit{x}=\left(\mathit{x},y\right)$, $\mathit{x}’$ is an auxiliary integral coordinate and $\delta$ is the Dirac delta function:

$$\delta \left(\mathit{x}-{\mathit{x}}^{\prime }\right)=\left\{\begin{array}{ll}\infty ,\hfill & \mathit{x}={\mathit{x}}^{\prime }\hfill \\ 0,\hfill & \mathit{x}\ne {\mathit{x}}^{\prime }\hfill \end{array}\right.$$

(2)

$${\displaystyle {\int }_{-\infty }^{+\infty }\delta \left(\mathit{x}-{\mathit{x}}^{\prime }\right)d\mathit{x}}=1$$

(3)

In numerical computation, usually a smooth function $W$, is used instead of the $\delta$ function:

$$〈f(\mathit{x})〉\approx {\displaystyle \int f}\left({\mathit{x}}^{\prime }\right)W\left(\mathit{x}-{\mathit{x}}^{\prime },h\right)\mathrm{d}{\mathit{x}}^{\prime }$$

(4)

where < > represents the SPH approximation operator and h represents the smoothing length of the smoothing function. The smoothing function must be normalized over its support domain:

$${\displaystyle \underset{\Omega }{\int }W(\mathit{x}-\mathit{x}’)\mathrm{d}\mathit{x}=1}$$

and the commonly used smooth function are as follows:

(1) Gaussian kernel function:

$$W(R,h)={\alpha }_d{e}^{-R^2}$$

(5)

where $R=\left|\mathit{x}-\mathit{x}’\right|/h$ and ${\alpha }_d=1/{\pi }^{1/2}h,1/\pi h^2\mathrm{and}1/{\pi }^{3/2}{h}^3$ in one-, two-, and three-dimensional space, respectively.

(2) Cubic spline function:

$$W(R,h)={\alpha }_d\times \left\{\begin{array}{ll}\frac{2}{3}-R^2+\frac{1}{2}{R}^3\hfill & 0\le R<1\hfill \\ \frac{1}{6}{(2-R)}^3\hfill & 1\le R<2\hfill \\ 0\begin{array}{c}\\ \end{array}\hfill & R\ge 2\hfill \end{array}\right.$$

(6)

where ${\alpha }_d=1/h,15/7\pi h^2\mathrm{and}3/2\pi h^3$ in one-, two-, and three-dimensional space, respectively.

(3) Quintic spline function:

$$W(R,h)={\alpha }_d\times \left\{\begin{array}{ll}{(3-R)}^5-6{(2-R)}^5+15{(1-R)}^5\hfill & 0\le R<1\hfill \\ {(3-R)}^5-6{(2-R)}^5\hfill & 1\le R<2\hfill \\ {(3-R)}^5\hfill & 2\le R<3\hfill \\ 0\hfill & R\ge 3\hfill \end{array}\right.$$

(7)

where ${\alpha }_d=120/h$, $7/478\pi h^2$ and $3/359\pi h^3$ in one-, two-, and three-dimensional space, respectively.

In SPH, the computational domain is discretized by a set of particles. Then the above integral form can be replaced by the following summation:

$$\langle f({\mathit{x}}_i)\rangle ={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}}f\left({\mathit{x}}_j\right)W\left({\mathit{x}}_i-{\mathit{x}}_j,h\right)$$

(8)

$${\nabla }_i{W}_{ij}=\frac{{\mathit{x}}_i-{\mathit{x}}_j}{\left|{\mathit{x}}_i-{\mathit{x}}_j\right|}\frac{\partial W_{ij}}{\partial r_{ij}}=\frac{{\mathit{x}}_i-{\mathit{x}}_j}{r_{ij}}\frac{\partial W_{ij}}{\partial r_{ij}}=\frac{{\mathit{x}}_{ij}}{r_{ij}}\frac{\partial W_{ij}}{\partial r_{ij}}$$

(9)

where $m_j$ and ${\rho }_j$ are mass and density of particles, respectively. The summation is performed over all the particles inside the support domain.

Similar to Equation (8), the first partial derivatives including gradient and divergence operators can also be expanded by SPH rules:

$$〈\nabla \cdot f\left({\mathit{x}}_i\right)〉={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}}f\left({\mathit{x}}_j\right)\cdot {\nabla }_i{W}_{ij}$$

(10)

Unfortunately, the numerical discretization precision of above equation is poor. According to Monaghan [15], the following formula is most widely used in existing literature:

$$〈\nabla \cdot f\left({\mathit{x}}_i\right)〉=\frac{1}{{\rho }_i}{\displaystyle \sum _{j=1}^N{m}_j}\left[f\left({\mathit{x}}_j\right)-f\left({\mathit{x}}_i\right)\right]\cdot {\nabla }_i{W}_{ij}$$

(11)

Discretizing the second-order derivative is also important since many PDEs in physics are of second-order. One approach is to use the Formula (4) twice [16,17,18], which involves the summation over all the particles. The other one is to employ second-order derivative of smoothing function [19], but this may result in large error for some specified problems. In order to improve the numerical accuracy of the second-order derivative in SPH, Revenga [20] proposed a flexible approach to discretize the second-order derivative and the general approximate formulas are as follows:

$${\left(\frac{{\partial }^{2}f}{\partial {\mathit{x}}^2}\right)}_i={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}}\left(4\frac{{\mathit{x}}_{ij}^2}{{\left|r_{ij}\right|}^2}-1\right)\left(f_i-f_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}$$

(12)

$${\left(\frac{{\partial }^{2}f}{\partial y^2}\right)}_i={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}}\left(4\frac{y_{ij}^2}{{\left|r_{ij}\right|}^2}-1\right)\left(f_i-f_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}$$

(13)

$${\left(\frac{{\partial }^{2}f}{\partial \mathit{x}\partial y}\right)}_i={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}}\left(4\frac{{\mathit{x}}_{ij}{y}_{ij}}{{\left|r_{ij}\right|}^2}\right)\left(f_i-f_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}$$

(14)

where $\left({\mathit{x}}_{ij},y_{ij}\right)=\left({\mathit{x}}_i-{\mathit{x}}_j,y_i-y_j\right)$.

2.2. The Governing Equations for Two-Dimensional Elastic Problems and Their Discrete SPH Forms

For plane stress problems, the equilibrium equations are:

$$\frac{\partial {\sigma }_{\mathit{x}}}{\partial \mathit{x}}+\frac{\partial {\tau }_{y\mathit{x}}}{\partial y}+F_{\mathit{x}}=0,\frac{\partial {\tau }_{\mathit{x}y}}{\partial \mathit{x}}+\frac{\partial {\sigma }_y}{\partial y}+F_y=0$$

(15)

where $\sigma$ and $\tau$ are normal and shear stresses, respectively, $F_{\mathit{x}}$ and $F_y$ are body forces.

The constitutive equations are Hooke’s law:

$${\epsilon }_{\mathit{x}}=\frac{1}{E}\left({\sigma }_{\mathit{x}}-\mu {\sigma }_y\right),{\epsilon }_y=\frac{1}{E}\left({\sigma }_y-\mu {\sigma }_{\mathit{x}}\right),{\gamma }_{\mathit{x}y}=\frac{2(1+\mu )}{E}{\tau }_{\mathit{x}y}$$

(16)

and the strain-displacement relations for plane stress problem are:

$${\epsilon }_{\mathit{x}}=\frac{\partial u}{\partial \mathit{x}},{\epsilon }_y=\frac{\partial v}{\partial y},{\gamma }_{\mathit{x}y}=\frac{\partial v}{\partial \mathit{x}}+\frac{\partial u}{\partial y}$$

(17)

where $\epsilon$ is normal strain and $\gamma$ are shear strain, $\mu$ and $E$ are Poisson’s ratio and Young’s modulus, respectively.

Expressing the equilibrium equations by the displacements yields the Navier equations of plane stress problems:

$$\begin{array}{l}\frac{E}{1-{\mu }^2}\left(\frac{{\partial }^{2}u}{\partial {\mathit{x}}^2}+\frac{1-\mu }{2}\frac{{\partial }^{2}u}{\partial y^2}+\frac{1+\mu }{2}\frac{{\partial }^{2}v}{\partial \mathit{x}\partial y}\right)+F_{\mathit{x}}=0,\\ \frac{E}{1-{\mu }^2}\left(\frac{{\partial }^{2}v}{\partial y^2}+\frac{1-\mu }{2}\frac{{\partial }^{2}v}{\partial {\mathit{x}}^2}+\frac{1+\mu }{2}\frac{{\partial }^{2}u}{\partial \mathit{x}\partial y}\right)+F_y=0\end{array}$$

(18)

For a well-posed elasticity problem, appropriate boundary conditions must be specified to ensure a unique solution. There are typically two kinds of boundary conditions in elasticity problems.

(1) Displacement boundary conditions. For these, the displacements are given on the boundary:

$$u=\overline{u},v=\overline{v}$$

(19)

where $\overline{u}$ and $\overline{v}$ are the components of the prescribed displacement on the boundary.

(2) Stress boundary conditions. For these, the tractions are specified on the boundary:

$$t_{\mathit{x}}=\overline{\mathit{x}},t_y=\overline{Y}$$

(20)

where $\overline{\mathit{x}}$ and $\overline{Y}$ are the components of the surface tractions on the boundary.

According to Equations (12)–(14), the discrete SPH forms of Equation (18) are as follows:

$$\begin{array}{c}\hfill \frac{E}{1-{\mu }^2}\left({\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}\left(4\frac{{\left({\mathit{x}}_{ij}\right)}^2}{{\left(r_{ij}\right)}^2}-1\right)}\left(u_i-u_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}+\frac{1-\mu }{2}{\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}\left(4\frac{{\left(y_{ij}\right)}^2}{{\left(r_{ij}\right)}^2}-1\right)}\left(u_i-u_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}+\right.\\ \hfill \left.\frac{1+\mu }{2}{\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}\left(4\frac{{\mathit{x}}_{ij}{y}_{ij}}{{\left(r_{ij}\right)}^2}\right)}\left(v_i-v_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}\right)+F_{\mathit{x}}=0\\ \hfill \frac{E}{1-{\mu }^2}\left({\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}\left(4\frac{{\left(y_{ij}\right)}^2}{{\left(r_{ij}\right)}^2}-1\right)}\left(v_i-v_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}+\frac{1-\mu }{2}{\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}\left(4\frac{{\left({\mathit{x}}_{ij}\right)}^2}{{\left(r_{ij}\right)}^2}-1\right)}\left(v_i-v_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}+\right.\\ \hfill \left.\frac{1+\mu }{2}{\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}\left(4\frac{{\mathit{x}}_{ij}{y}_{ij}}{{\left(r_{ij}\right)}^2}\right)}\left(u_i-u_j\right)\frac{1}{r_{ij}}\frac{\partial W}{\partial r}\right)+F_y=0\end{array}$$

(21)

2.3. Treatment of the Boundary Conditions

In this paper, we propose a SPH-GFDM, namely the SPH-Generalized Finite Difference Method coupled approach to deal with the boundary condition implementation difficulties when SPH is applied to elasticity problems.

The generalized finite difference method (GFDM) is another meshless collocation method. It was first proposed by Liszka et al. [21,22], later improved and extended by many other authors [23,24,25,26]. Before this study, this method has been successfully applied to heat transfer, elastic analysis and other engineering problems [27,28,29,30,31,32,33,34]. The main idea of the method is to combine the moving-least squares (MLS) approximation and the Taylor series expansion to derive explicit formulae for the required partial derivatives of unknown field variables. Without loss of generality, let us consider the following second-order linear differential equation in 2D domain:

$$a_1\frac{\partial U}{\partial \mathit{x}}+a_2\frac{\partial U}{\partial y}+a_3\frac{{\partial }^{2}U}{\partial {\mathit{x}}^2}+a_4\frac{{\partial }^{2}U}{\partial y^2}+a_5\frac{{\partial }^{2}U}{\partial \mathit{x}\partial y}=f(\mathit{x},y)$$

(22)

where $f\left(\mathit{x},y\right)$ is a known function, $U\left(\mathit{x},y\right)$ is the field function, $a_i\left(i=1,2,3,4,5\right)$ is the coefficient, and $a_i$ is usually a given function or constant.

To obtain the explicit GFDM formulae for partial derivatives in a given differential equation and boundary conditions, a regular or irregular set of nodes is scattered in the computational domain and on the boundary. For a specified node $\left({\mathit{x}}_0,y_0\right)$, namely the central node, N nearest nodes $\left({\mathit{x}}_i,y_i\right)\left(i=1,2,\dots ,N\right)$, called supporting nodes or neighbors, will be identified within a prescribed region (usually a circle with a given radius) to form the support domain. It is similar to the concept of support domain in SPH. Every node inside the computational domain and on the boundary has an associated support domain.

Let $U_0$ be the value of the field function at the central node $\left({\mathit{x}}_0,y_0\right)$ and $U_i\left(i=1,2,\dots ,N\right)$ are field function values at the rest of the nodes inside the support domain. Expanding the field function value $U_i$ around the central node $\left({\mathit{x}}_0,y_0\right)$ using Taylor series expansion, we will get:

$$U_i=U_0+h_i\frac{\partial U_0}{\partial \mathit{x}}+k_i\frac{\partial U_0}{\partial y}+\frac{h_i^2}{2}\frac{{\partial }^2{U}_0}{\partial {\mathit{x}}^2}+\frac{k_i^2}{2}\frac{{\partial }^2{U}_0}{\partial y^2}+h_i{k}_i\frac{{\partial }^2{U}_0}{\partial \mathit{x}\partial y}+\dots ,i=1,2,\dots ,N$$

(23)

where $h_i={\mathit{x}}_i-{\mathit{x}}_0$, $k_i=y_i-y_0$.

By truncating the Taylor series after the second-order derivatives (this depends on the partial differential equation being solved. For example, if the PDE is of forth-order, the Taylor series should be truncated after the forth-order derivatives), we then define the following residual function B(U):

$$B\left(U\right)={\displaystyle \sum _{i=1}^N{\left[\left(U_0-U_i+h_i\frac{\partial U_0}{\partial \mathit{x}}+k_i\frac{\partial U_0}{\partial y}+\frac{h_i^2}{2}\frac{{\partial }^2{U}_0}{\partial {\mathit{x}}^2}+\frac{k_i^2}{2}\frac{{\partial }^2{U}_0}{\partial y^2}+h_i{k}_i\frac{{\partial }^2{U}_0}{\partial \mathit{x}\partial y}\right){\omega }_i\right]}^2}$$

(24)

where ${\omega }_i=\omega (h_i,k_i)$ are the weighting coefficients at point $\left({\mathit{x}}_i,y_i\right)$ with respect to the central node and the weighting coefficients are usually taken to be:

$${\omega }_i=1-6S^2+8S^3-3S^4,d_i\le d_m$$

(25)

where $d_i=\sqrt{h_i^2+k_i^2}$ is the distance between points $\left({\mathit{x}}_0,y_0\right)$ and $\left({\mathit{x}}_i,y_i\right)$, and $d_m$ is the radius of the support domain of the selected node, $S=d_i/d_m$. The weighting coefficients are smaller for the nodes far away from the central node within the support domain than those nearby. Other choices of weighting functions are also possible, and the reader can find them in Refs [7,9].

For brevity, we define:

$$\begin{array}{l}{D}_U={\left\{\frac{\partial U_0}{\partial \mathit{x}},\frac{\partial U_0}{\partial y},\frac{{\partial }^2{U}_0}{\partial {\mathit{x}}^2},\frac{{\partial }^2{U}_0}{\partial y^2},\frac{{\partial }^2{U}_0}{\partial \mathit{x}\partial y}\right\}}^T\\ p_i=\left\{h_i,k_i,h_i^2/2,k_i^2/2,h_i{k}_i\right\},i=1,2,\dots ,N\\ P={\left[p_1,p_2,\dots ,p_N\right]}^T\end{array}$$

(26)

Then the residual function $B\left(U\right)$ can be written in the following matrix form:

$$B(U)={\left(P{D}_U+U_0-U\right)}^{T}W\left(P{D}_U+U_0-U\right)$$

(27)

where

$$U={\left\{U_1,U_2,\cdots ,U_N\right\}}^T,U_0={\left\{U_0,U_0,\cdots ,U_0\right\}}^T,W=\mathrm{diag}\left({\omega }_1^2,{\omega }_2^2,\cdots ,{\omega }_N^2\right)$$

(28)

By minimizing $B\left(U\right)$ with respect to the partial derivatives $D_U$:

$$\partial B\left(U\right)/\partial D_U=0$$

(29)

the following linear equation system will be obtained for node $\left({\mathit{x}}_i,y_i\right)$ in the domain:

$$A{D}_U=b$$

(30)

where

$$A=P^{T}WP$$

(31)

and

$$b=P^{T}W\left(U-U_0\right)$$

(32)

According to above equations, the partial derivative vector $D_U$ can be expressed as:

$$D_U=A^{-1}b=A^{-1}{P}^{T}W\left(U-U_0\right)$$

(33)

For the central node $\left({\mathit{x}}_0,y_0\right)$ in the computational domain, the partial derivatives of the unknown field function can be approximated by a linear combination of field function values at its neighboring nodes within the support domain. Substituting Equation (33) Into the Governing Equation (22) leads to:

$$m_0{U}_0+{\displaystyle \sum _{i=1}^N{m}_i{U}_i}-f\left({\mathit{x}}_0,y_0\right)=0$$

(34)

where the coefficients $m_0$ and $m_i^{}\left(i=1,2,\dots ,N\right)$ can be acquired from Equation (22). For all nodes within the computational domain, we repeat the same steps and a system of linear equations can be established. By solving this equation system, we will get the numerical results of the field function $U$ at each node within the computational domain.

Taking advantage of both SPH and GFDM, we propose a SPH-GFDM coupled method to solve PDEs. The main steps of our SPH-GFDM coupled method are as follows:

(1) Discretize the entire computational domain by a set of arbitrarily distributed nodes inside the domain and on the boundaries. The nodes distribution should be as regular as possible.

(2) For the nodes far away from the boundaries whose support domain are complete, the SPH discretization form is accurate enough. SPH will thus be adopted. As for the nodes near and on the boundaries, their support domains are truncated by the boundary. Hence, the GFDM approximation will be used to discretize the governing equation and implement the boundary conditions. Then, every node inside the computational domain and on the boundaries will correspond to an algebraic equation. This results in a large-scale system of linear equations.

(3) Solve the system of linear algebraic equations and we will get the numerical solutions of the problem on all discrete nodes.

It is necessary to note that, however, the approach this paper proposed is based on the condition that the forces in the calculation domain reach an equilibrium state, and the equations are also deduced on this basis. For the transient problems, the smoothed particles in the calculation domain are motioning at an accelerated speed, which is another case for further study.

3. Results and Discussion

In order to verify the accuracy of GFDM-SPH coupled method proposed in the previous section, two benchmark test examples are examined for both symmetric and asymmetric computational domains. The numerical solutions are compared with those obtained by a commercial FEM code (COMSOL), in which a high-quality mesh (small skewness and proper aspect ratio) is adopted to ensure the accuracy of the numerical results. A mesh convergence study was conducted to ensure the results did not vary with mesh size.

3.1. A Rectangular Metal Plate with Two Ends Fixed under Prescribed Uniform Loads

As shown in Figure 1a, a 4 m × 1 m rectangular plate whose two ends are fixed is present, its upper and lower sides are under a uniformly distributed load q = 50 N/m. The material constants are ρ = 7800 kg/m³, E = 80 GPa and μ = 0.32, respectively. This is a mixed-boundary problem; the upper and lower sides of the plate are prescribed for Neumann boundary conditions and the left and right sides of it are assigned with Dirichlet boundary conditions. As shown in Figure 1d, the blue points in the inner region represent the nodes in which the governing equations are discretized by SPH, and the red points near and on the boundary are those in which the governing equations and boundary conditions are discretized by GFDM. Figure 1b shows the location of the sample point P(1.0, 0.5) and line AB and Figure 1c is the mesh generated by COMSOL.

Figure 2 shows the contours of stress and displacement obtained by FEM. Figure 3 shows those by SPH-GFDM. Figure 4 shows the y-direction displacement v and x-direction displacement u on line AB obtained by SPH-GFDM and FEM.

Table 1. The x and y-direction maximum displacements (Unit: m).

	SPH-GFDM	FEM
${\left|u\right|}_{ma\mathit{x}}$	4.08 × 10−9	4.12 × 10−9
${\left|v\right|}_{ma\mathit{x}}$	1.75 × 10−8	1.76 × 10−8

displays the absolute value of the maximum displacement in x and y-direction given by SPH-GFDM and FEM and

Table 2. The stress values in sample point P (Unit: Pa).

	σx	σy	τxy	σSig_eq
SPH-GFDM	2.28 × 10−6	7.30 × 10−7	−148.82	257.76
FEM	0.00	0.00	−149.45	263.76

shows the stress values in sample point P obtained by these two methods.

As show in these figures and tables, the numerical results of the displacement field and stress field obtained by the SPH-GFDM coupled method in this simple case are accurate compared to the results obtained by FEM of a high-quality mesh. The maximum relative errors are not more than 3%.

3.2. A L-Shape Plate with Two Ends Fixed under Prescribed Uniform Loads

As shown in Figure 5a, a L-shape plate whose left and upper ends are fixed is present, its right and lower sides are under a uniformly distributed load q = 50 N/m. The material constants are ρ = 7800 kg/m³, E = 80 GPa and μ = 0.32, respectively. As shown in Figure 5d, the blue points in the inner region represent the nodes in which the governing equations are discretized by SPH, and the red points near and on the boundary are those in which the governing equations and boundary conditions are discretized by GFDM. Figure 5b shows the location of the sample points A(0.0, 1.0), B(4.0, 1.0), P₁(2.0, 1.8), P₂(2.2, 1.8), P₃(2.2, 2.0) and line AB and Figure 5c is the mesh generated by COMSOL.

Figure 6 shows the contours of stress and displacement obtained by FEM. Figure 7 shows those by SPH-GFDM. Figure 8 shows the y-direction displacement u and x-direction displacement v on line AB obtained by SPH-GFDM and FEM.

Table 3. The x and y-direction maximum displacements (Unit: m).

	SPH-GFDM	FEM
${\left|u\right|}_{ma\mathit{x}}$	2.20 × 10−9	2.21 × 10−9
${\left|v\right|}_{ma\mathit{x}}$	2.20 × 10−9	2.21 × 10−9

displays the absolute value of the maximum displacement in x and y-direction given by SPH-GFDM and FEM and

Table 4. The stress values in sample point P₁ (Unit: Pa).

	σx	σy	τxy	σSig_eq
SPH-GFDM	40.28	37.41	16.23	48.28
FEM	42.77	37.08	18.25	51.16

,

Table 5. The stress values in sample point P₂ (Unit: Pa).

	σx	σy	τxy	σSig_eq
SPH-GFDM	43.25	43.25	9.93	45.17
FEM	42.46	42.44	11.83	47.14

and

Table 6. The stress values in sample point P₃ (Unit: Pa).

	σx	σy	τxy	σSig_eq
SPH-GFDM	37.41	40.28	16.23	48.28
FEM	37.07	42.74	18.17	51.06

show the stress values in sample points P₁, P_2, and P₃ obtained by these two methods.

As the figures and data in the tables illustrated, for L-shape plate, the y-direction displacement u and x-direction displacement v on line AB and the normal stress at sample points P₁, P_2, and P₃ obtained by SPH-GFDM match well with the results obtained by FEM, within the maximum relative error at 5.8%. However, the deviation of shear stresses and Von-Mises stresses obtained by these two methods can reach 16.1%, probably due to the asymmetry of this L-shape plate, which will be further analyzed in the later section.

3.3. A Square Plate with A Hole in the Center Whose Left Side Is Fixed and the Right Side Is under Uniform Normal Loads

As shown in Figure 9a, a 4 m × 4 m square plate with a 1 m × 1 m hole in the center is under a uniform loads q = 50 N/m on its right side, its left side is totally fixed and the surface of the inner hole is free. The material constants are ρ = 7800 kg/m³, E = 80 GPa and μ = 0.32, respectively. For simplicity, the influence of the gravity is neglected and the distribution of scattered nodes is shown in Figure 9d. The blue points in the inner region represent the nodes in which the governing equations are discretized by SPH. The red points near the boundary are those in which the governing equations. and boundary conditions are discretized by GFDM. Figure 9b shows the location of the sample points A(2, 3.25) and B(0.75, 2) and Figure 9c is the mesh generated by COMSOL.

Compared with the previous test examples, this problem is relatively more complex. The computational domain of this problem is a multiple connected region. The left side of the plate is fixed, which is a Dirichlet boundary condition. The other sides of the plate are assigned with Neumann boundary conditions, including the free boundaries on which the tractions are zero.

Figure 10a–f are the contours of stress and displacement obtained by FEM and Figure 11a–f are those given by SPH-GFDM.

Table 7. The x and y-direction maximum displacements (Unit: m).

	SPH-GFDM	FEM
${\left|u\right|}_{ma\mathit{x}}$	3.26 × 10−9	3.32 × 10−9
${\left|v\right|}_{ma\mathit{x}}$	6.33 × 10−10	6.34 × 10−10

shows the absolute value of the maximum displacement in x and y-direction given by SPH-GFDM and FEM.

Table 8. The stress values in sample point A (Unit: Pa).

	σx	σy	τxy	σSig_eq
SPH-GFDM	−62.57	−4.21	2.57	60.25
FEM	−63.29	−4.32	2.76	61.43

and

Table 9. The stress values in sample point B (Unit: Pa).

	σx	σy	τxy	σSig_eq
SPH-GFDM	−21.68	−9.80	0.00	18.41
FEM	−22.81	−10.00	0.00	19.80

shows the stress values in sample points A and B obtained by SPH-CFDM and FEM, respectively. Roughly speaking, the numerical results obtained by the SPH-GFDM method are similar to those given by FEM with a high-quality mesh, as shown in Figure 2, Figure 3, Figure 6, Figure 7, Figure 10, and Figure 11. But there still exist some differences: the maximum x-direction displacement given by SPH-GFDM is relatively smaller than that given by COMSOL, and the y-direction normal stress near the left side of the square hole obtained by SPH-GFDM is also smaller than that given by FEM. The possible reasons for this may result from two aspects: Firstly, there exist stress singularities in this example. Theoretically speaking, when the plastic behaviors of the material are not considered here, the stress values in the four corner points of the inner square hole are infinite. Secondly, the numerical results given by COMSOL (a FEM-based software) are ‘weak’ solutions and those obtained by SPH-GFDM are ‘strong’ solutions. As we know, FEM does not directly discretize the governing equation and the boundary conditions. Instead, it depends on the equivalent integral weak forms of the problem. However, both SPH and GFDM are meshless direct collocation methods and in SPH-GFDM coupled method. We discretize the partial derivatives in the governing equation and the boundary conditions by a linear combination of the field function values of the nodes within the support domain. For elastic problems, the ‘weak’ solutions are usually more stable and accurate than the ‘strong’ solutions. These two reasons may result in the differences.

4. Conclusions

In this study, we demonstrate the applicability of the SPH–GFDM coupled method for the analysis of 2-D elastic problems. The main idea of this method is to numerically solve PDEs with prescribed boundary conditions. The partial derivatives of the field functions at nodes within the interior region of the computational domain and near the boundary are discretized by the SPH and GFDM approximations, respectively. Results are obtained by solving the combined algebra equations associated with the two approximations. Due to the least square essence of GFDM, Dirichlet and Neumann boundary conditions in GFDM can be easily implemented with relatively high accuracy. This improves the overall numerical accuracy of the SPH-GFDM coupled method. Two benchmark examples of both symmetric and asymmetric computational domains are given to verify the stability and accuracy of the SPH-GFDM coupled method. It is demonstrated that, as the results and comparisons show, this method behaves as reliable as FEM with high-quality meshes in some simple cases in which the computational domains are simply connected, which means better computational efficiency and less time consuming. Yet, considering to the relatively high deviation of shear stresses and Von-Mises stresses obtained by SPH-GFDM and FEM in some cases, the proposed method still needs to be improved, especially when applied to the case where the computational domains are multiply connected and the boundary shapes are more complex.