Meshless Generalized Finite Difference Method Based on Nonlocal Differential Operators for Numerical Simulation of Elastostatics

Abstract

This study proposes an innovative meshless approach that merges the peridynamic differential operator (PDDO) with the generalized finite difference method (GFDM). Based on the PDDO theory, this method introduces a new nonlocal differential operator that aims to reduce the pre-assumption required for the PDDO method and simplify the calculation process. By discretizing through the particle approximation method, this technique proficiently preserves the PDDO’s nonlocal features, enhancing the numerical simulation’s flexibility and usability. Through the numerical simulation of classical elastic static problems, this article focuses on the evaluation of the calculation accuracy, calculation efficiency, robustness, and convergence of the method. This method is significantly stronger than the finite element method in many performance indicators. In fact, this study demonstrates the practicability and superiority of the proposed method in the field of elastic statics and provides a new approach to more complex problems.

1. Introduction

The numerical analysis of elastic statics plays a vital role in civil and mechanical engineering. Engineers can use the knowledge of elastic statics to predict the stress, strain, and displacement of structural members through numerical simulation before actual construction to evaluate the feasibility and safety of a project [1].

The finite element method (FEM) [2,3,4,5] is widely used in numerical analysis. However, when faced with problems involving large deformations and high-speed impacts, the grids may deform, leading to non-convergence. At the same time, when dealing with some complex problems, the computational cost of the FEM is usually higher compared to that of mesh-free methods [6]. The boundary element method (BEM) [7,8,9,10] focuses more on the processing and analysis of the boundary, reducing the discretization requirements of the internal points and improving the computational efficiency. The disadvantage is its limited ability to deal with nonlinear problems and the high computational cost of finding the internal solution. As one of the most common meshless methods, smoothed particle hydrodynamics (SPH) [11,12] has great advantages in dealing with complex geometric boundaries and large deformation problems. However, special treatment is usually required for boundary conditions, such as the arrangement of virtual particles [13]. Although this approach improves the accuracy and stability of the method, it significantly increases its computational cost. The Petrov–Galerkin method [14] also faces similar problems to those of the SPH method.

Based on the limitations of the aforementioned numerical simulation methods, this study introduces a new type of mesh-free method. This method refers to the theoretical basis of the two advanced calculation methods of the peridynamic differential operator (PDDO) and generalized finite difference method (GFDM) and combines their advantages. The primary objective of this study is to integrate the concept of Taylor expansion and the formulation of nonlocal integrals within the PDDO framework, leveraging the particle approximation method [15] frequently employed in SPH methods for discretizing these integrals. Subsequently, we aim to compute the necessary differential operators formulated through a weighted summation approach. Using differential operators to obtain the required partial derivatives, this step is crucial for the entire numerical simulation process because the accuracy of the differential operator directly determines the accuracy and reliability of the problem solving. By substituting these exact partial derivatives into the discrete equations constructed using the GFDM and solving the discrete equations, we can obtain physical quantities such as displacement and stress.

Peridynamics (PD) is a theory proposed by Silling et al. [16,17,18] to simulate the phenomenon of continuum mechanics. Unlike the traditional continuum mechanics method, PD describes the relationship between materials through the long-range interaction between material points, thus overcoming the difficulties in dealing with material fracture. Therefore, it has great advantages in dealing with cracks and fracture phenomena occurring inside the material. The peridynamic differential operator is a meshless method based on the nonlocal idea of PD proposed by Madenci [19,20] in 2016. Using the characteristics of Taylor series expansion and PD function, this method transforms the local differential into a widely applicable nonlocal integral form, realizes the numerical differential, and effectively avoids the discontinuity or singularity problems that may occur in the simulation process. The PDDO provides a nonlocal integral framework for this study, which makes approximations of the derivatives of physical quantities possible using the integral form without direct differentiation. When dealing with discontinuous problems, traditional differential methods encounter situations where the entire solution domain is not continuously derivable. In contrast, using the integral form to approximate the derivative is less susceptible to local discontinuities and can handle such discontinuities more precisely. Compared to the differential method, the integral form technique is more stable. This method has been widely used to simulate crack propagation [21], coupled stress [22], heat conduction [23], and other phenomena.

The generalized finite difference method, first proposed by Liszka [24] in 1980, is a meshless method with local control domain characteristics. It realizes the numerical solution of partial differential equations with the difference approximations of nodes in local regions. Since the GFDM is a meshless method that can arrange irregular nodes, it has unique advantages in dealing with problems with complex geometric boundaries. Many researchers have improved and perfected the GFDM through in-depth research over time [25,26,27], making its theoretical basis more mature and its application scope more extensive. The research of Hidayat [28] further demonstrates the effectiveness of the GFDM in solving two-dimensional elastic static problems. The GFDM is also widely used in the solution of partial differential equations [29] and the analysis of complex problems such as geotechnical engineering [30,31]. In these applications, the GFDM shows its unique advantages in dealing with irregular boundaries, multi-physics coupling, and nonlinear problems, thus playing a vital role in scientific research and engineering practice.

To evaluate the accuracy and efficiency of the method proposed in this article, we selected three classical examples for numerical simulation. In addition, the robustness and convergence of this method are discussed by changing the corresponding physical parameters [32]. Through the study of the above characteristics, the superiority of this method and its potential in simulating more complex problems are further illustrated.

2. Theoretical Knowledge

2.1. Deriving Partial Derivatives Using Nonlocal Differential Operators

In this section, the newly proposed nonlocal differential operator presented in this article is used to solve for partial derivatives. It mainly introduces how to use the Taylor series expansion and nonlocal properties of the PDDO, utilize the particle approximation method to discretize the corresponding differential operators, and incorporate the constructed partial derivatives into the GFDM of elastic constitutive relations to obtain variables such as displacement and stress.

Consider a field $f =f\left(X\right)$, where node $X$ serves as the point source, as illustrated in Figure 1. The influence range $kh$ is determined by the internal nodes $X^{\prime }$ within the vicinity of $X$. In this context, $h$ represents the smoothing length, and $k$ is the kernel coefficient. The distance between these points can be defined as $\mathit{\xi }=({\mathit{\xi }}_x,{\mathit{\xi }}_y) =X{-X}^{\prime }$ with $\left(\mathit{\xi } \le h\right)$. In meshless methods, the distances between these nodes and the size of the support domain are usually considered infinitesimal.

For a two-dimensional problem, where $X=\left(x,y\right)$, the field function $\mathit{f}$(${\mathit{X}}^{\mathbf{\prime }}$) = $\mathit{f}\left(\mathit{X}+\mathit{\xi }\right)$ is expanded using a second-order Taylor series:

$$\begin{array}{c}f\left(X+\mathit{\xi }\right)={\displaystyle \sum _{n_x=0}^2}{\displaystyle \sum _{n_y=0}^{2-n_x}}\frac{1}{n_x!n_y!}{\mathit{\xi }}_x^{n_x}{\mathit{\xi }}_y^{n_y}\frac{{\mathit{\partial }}^{n_x+n_y}f\left(X\right)}{\mathit{\partial }{x}^{n_x}\mathit{\partial }{y}^{n_y}}+T\left(X\right)\end{array}$$

(1)

where $n_x$ and $n_y$ are 0, 1, or 2 and $T\left(X\right)$ represents the higher-order remainder term.

Taking the second-order Taylor expansion at a point $X_j=X +{ \mathit{\xi }}^j$ within the influence domain of $X$ as an example,

$$\begin{array}{c}f\left(X+{\mathit{\xi }}^j\right)=f\left(X\right)+{\mathit{\xi }}_x^j\frac{\mathit{\partial }f\left(X\right)}{\mathit{\partial }x}+{\mathit{\xi }}_y^j\frac{\mathit{\partial }f\left(X\right)}{\mathit{\partial }y}+\frac{1}{2}{\mathit{\xi }}_x^j{\mathit{\xi }}_x^j\frac{{\mathit{\partial }}^{2}f\left(X\right)}{{\mathit{\partial }x}^2}+{ \mathit{\xi }}_x^j{\mathit{\xi }}_y^j\frac{{\mathit{\partial }}^{2}f\left(X\right)}{\mathit{\partial }x\mathit{\partial }y}+\frac{1}{2}{\mathit{\xi }}_y^j{\mathit{\xi }}_y^j\frac{{\mathit{\partial }}^{2}f\left(X\right)}{{\mathit{\partial }y}^2}+T\left(X\right) \end{array}$$

(2)

we define the following:

$$\begin{array}{c}{\mathit{d}}_j=[1,{\mathit{\xi }}_x^j,{\mathit{\xi }}_y^j,{\mathit{\xi }}_x^j{\mathit{\xi }}_x^j,{\mathit{\xi }}_x^j{\mathit{\xi }}_y^j,{\mathit{\xi }}_y^j{\mathit{\xi }}_y^j]\end{array}$$

(3)

$$\begin{array}{c}\mathit{C}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& \frac{1}{2}& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& \frac{1}{2}\end{array}\right]\end{array}$$

(4)

$$\begin{array}{c}𝛛\mathit{f}={\left[f \left(X\right),\frac{\mathit{\partial }f \left(X\right)}{\mathit{\partial }x},\frac{\mathit{\partial }f \left(X\right)}{\mathit{\partial }y},\frac{{\mathit{\partial }}^{2}f \left(X\right)}{{\mathit{\partial }x}^2},\frac{{\mathit{\partial }}^{2}f \left(X\right)}{\mathit{\partial }x\mathit{\partial }y},\frac{{\mathit{\partial }}^{2}f \left(X\right)}{{\mathit{\partial }y}^2}\right]}^T\end{array}$$

(5)

Then, Equation (2) can be simplified to:

$$\begin{array}{c}f\left(X_j\right)={\mathit{d}}_j\mathit{C}\times 𝛛\mathit{f}+T\left(X\right)\end{array}$$

(6)

Ignoring the higher-order terms, we apply weighted integration to both sides of Equation (6):

$$\begin{array}{c}\mathit{\int }{w}_j {\mathit{d}}_{\mathit{j}}^{\mathit{T}}f\left(X_j\right)d\Omega =\mathit{\int }{w}_j{\mathit{d}}_{\mathit{j}}^{\mathit{T}}{\mathit{d}}_j\mathit{C}d\Omega \times 𝛛\mathit{f}\end{array}$$

(7)

where $d\Omega$ in a two-dimensional problem refers to the area of the influence domain of the point source $X$ and $w_j=w (X-X_j, h)$ represents the weight function between the point source $X$ and a point $X_j$ within the computational domain. In the SPH method, this weight function can typically be a quintic spline kernel function, cubic spline kernel function, or Gaussian kernel function [15], with the quintic spline kernel function being presented as follows:

$$\begin{array}{c}w\left(R,h\right)=a_d\left\{\begin{array}{c}{\left(3-R\right)}^5-6({2-R)}^5+15({1-R)}^5 R\le 1\hfill \\ ({3-R)}^5-6({2-R)}^5\hspace{1em}\hspace{1em}\hspace{1em} 1<R\le 2\hfill \\ ({3-R)}^5\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}2<R\le 3\hfill \\ \hspace{1em} 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}3<R\hfill \end{array}\right.\end{array}$$

(8)

$$\begin{array}{c}R=\frac{\left|X-{ X}_j\right|}{h}\end{array}$$

(9)

where $a_d$ in one-, two-, or three-dimensional spaces is 1/120$h$, 7/478π$h^2$, or 3/359π$h^3$, respectively.

If we perform a second-order Taylor expansion for all particles within the entire computational domain with particle $X$ as the point source (assuming $M$ such points) and ignore the higher-order terms $T(X),$ we can obtain a system of equations with $𝛛\mathit{f}$ as the unknown variable:

$$\begin{array}{c}\mathit{\int }{\mathit{D}}^T\mathit{W}\mathit{F}d\Omega =\mathit{\int }{\mathit{D}}^T\mathit{W}\mathit{D}\mathit{C}d\Omega 𝛛\mathit{f}\end{array}$$

(10)

Inverting Equation (10) yields:

$$\begin{array}{c}𝛛\mathit{f}={\left(\mathit{\int }{\mathit{D}}^T\mathit{W}\mathit{D}\mathit{C}d\Omega \right)}^{-1}\mathit{\int }{\mathit{D}}^T\mathit{W}\mathit{F}d\Omega \end{array}$$

(11)

where:

$$\begin{array}{c}\mathit{D}=\left[\begin{array}{c}{d}_1\\ d_2\\ d_3\\ ⋮\\ d_M\end{array}\right]\end{array}$$

(12)

$$\begin{array}{c}\mathit{F}=\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\\ ⋮\\ f_M\end{array}\right]\end{array}$$

(13)

$$\begin{array}{c}\mathit{W}=\left[\begin{array}{ccccc}{w}_1& 0& 0& \dots & 0\\ 0& w_2& 0& \dots & 0\\ 0& 0& w_3& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ 0& 0& 0& 0& w_M\end{array}\right]\end{array}$$

(14)

In the process of converting integration into summation for discretization, the idea of the particle approximation method from the SPH method is referenced:

$$\begin{array}{c}f\left(X\right)={\displaystyle \mathit{\int }}f \left(X_j\right)w\left(X-X_j,h\right)d \Omega ={\displaystyle \sum _{j=1}^M}f \left(X_j\right)w\left(X-X_j,h\right)d{V}_j\end{array}$$

(15)

$$\begin{array}{c}{\displaystyle \sum _{j=1}^M}f \left(X_j\right)w\left(X-X_j,h\right)d{V}_j={\displaystyle \sum _{j=1}^M}f \left(X_j\right)w\left(X-X_j,h\right)\frac{m_j}{{\rho }_j}\end{array}$$

(16)

where $m_j$ and ${\rho }_j$ represent the mass and density, respectively, of a point $X_j$ within the computational domain with $X$ as the point source.

Combining Equations (11), (15), and (16) yields:

$$\begin{array}{c}𝛛\mathit{f}={ \left({\displaystyle \sum _{j=1}^M}{\mathit{D}}^T\mathit{W}\mathit{V}\mathit{D}\mathit{C}\right)}^{-1}{\displaystyle \sum _{j=1}^M}{\mathit{D}}^T\mathit{WVF}\end{array}$$

(17)

where $\mathit{V}$ represents:

$$\begin{array}{c}\mathit{V}=\left[\begin{array}{ccccc}\frac{m_1}{{\rho }_1}& 0& 0& \dots & 0\\ 0& \frac{m_2}{{\rho }_2}& 0& \dots & 0\\ 0& 0& \frac{m_3}{{\rho }_3}& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & 0\\ 0& 0& 0& 0& \frac{m_M}{{\rho }_M}\end{array}\right]\end{array}$$

(18)

From Equation (17), the numerical differentiation operator $\mathit{G}$ with $X$ as the point source can be obtained:

$$\begin{array}{c}\mathit{G}={ \left({\displaystyle \sum _{j=1}^M}{\mathit{D}}^T\mathit{WVDC}\right)}^{-1}{\displaystyle \sum _{j=1}^M}{\mathit{D}}^T\mathit{WV}={[1,\frac{𝛛}{𝛛x},\frac{𝛛}{𝛛y},\frac{{𝛛}^2}{{𝛛x}^2},\frac{{𝛛}^2}{𝛛x𝛛y},\frac{{𝛛}^2}{{𝛛y}^2}]}^{\mathit{T}}\end{array}$$

(19)

2.2. Elastic Constitutive of the Generalized Finite Difference Method

The governing equations for the elastic statics of a two-dimensional stress field, considering a constitutive model, where $\Omega$, ${\Gamma }_u$, and ${\Gamma }_t$ represent the computational domain, displacement boundary, and traction boundary respectively, can be expressed as:

$$\begin{array}{c} {\sigma }_{\mathit{ij},j }=-{\overline{f}}_i\mathit{in}\Omega \end{array}$$

(20)

$$\begin{array}{c}{\sigma }_{\mathit{ij}}{n}_j={\overline{t}}_i\mathit{on}{\Gamma }_t\end{array}$$

(21)

$$\begin{array}{c} u_i={\overline{u}}_i\mathit{on}{\Gamma }_u\end{array}$$

(22)

where ${\sigma }_{\mathit{ij}}$ is the stress tensor, ${\overline{f}}_i$ represents the known components of body force, $n_j$ and ${\overline{t}}_i$, respectively, denote the cosine of the exterior normal and known components of traction on the boundary ${\Gamma }_t$, $u_i$ represents the displacement components, and ${\overline{u}}_i$ represents the known displacement components on the boundary ${\Gamma }_u$, where $i=x,y$ and $j=x,y.$

Taking a two-dimensional plane stress problem as an example, Equations (20)–(22) can be expressed in terms of displacement components:

$$\begin{array}{c}\begin{array}{c}\frac{E}{1-v^2}\left[\begin{array}{cc}\frac{{\mathit{\partial }}^2}{{\mathit{\partial }x}^2}+b\frac{{\mathit{\partial }}^2}{{\mathit{\partial }y}^2}& a\frac{{\mathit{\partial }}^2}{\mathit{\partial }x\mathit{\partial }y}\\ a\frac{{\mathit{\partial }}^2}{\mathit{\partial }y\mathit{\partial }x}& b\frac{{\mathit{\partial }}^2}{{\mathit{\partial }x}^2}+\frac{{\mathit{\partial }}^2}{{\mathit{\partial }y}^2}\end{array}\right]\left[\begin{array}{c}{u}_x\\ u_y\end{array}\right]=-\left[\begin{array}{c}{\overline{f}}_x\\ {\overline{f}}_y\end{array}\right]\end{array}\end{array}$$

(23)

$$\begin{array}{c}\frac{E}{1-v^2}\left[\begin{array}{cc}{n}_x\frac{\mathit{\partial }}{\mathit{\partial }x}+n_{y}b\frac{\mathit{\partial }}{\mathit{\partial }y}& n_{y}b\frac{\mathit{\partial }}{\mathit{\partial }x}+n_{x}v\frac{\mathit{\partial }}{\mathit{\partial }y}\\ n_{y}v\frac{\mathit{\partial }}{\mathit{\partial }x}+n_{x}b\frac{\mathit{\partial }}{\mathit{\partial }y}& n_{x}b\frac{\mathit{\partial }}{\mathit{\partial }x}+n_y\frac{\mathit{\partial }}{\mathit{\partial }y}\end{array}\right]\left[\begin{array}{c}{u}_x\\ u_y\end{array}\right]=\left[\begin{array}{c}{\overline{t}}_x\\ {\overline{t}}_y\end{array}\right]\end{array}$$

(24)

$$\begin{array}{c}\left[\begin{array}{c}{u}_x\\ u_y\end{array}\right]=\left[\begin{array}{c}{\overline{u}}_x\\ {\overline{u}}_y\end{array}\right]\end{array}$$

(25)

where $E$ is the elastic modulus, $v$ represents the Poisson’s ratio, and $a$ and $b$ are defined as:

$$\begin{array}{c}a=\frac{1+v}{2}\end{array}$$

(26)

$$\begin{array}{c}b=\frac{1-v}{2}\end{array}$$

(27)

The fundamental goal of the GFDM is to approximate partial derivatives at each node using the function values of surrounding nodes. This approximation is achieved by solving for weighted coefficients, making the partial derivatives at a specific node approximate to a linear combination of the function values of the surrounding nodes. The main formula can be expressed as:

$$\begin{array}{c}{\left.\frac{{\mathit{\partial }}^{n}f}{\mathit{\partial }{x}^m\mathit{\partial }{y}^{n-m}}\right|}_{\left(x,y\right)=\left(x_i,y_i\right)}\approx {\displaystyle \sum _{j=1}^M}{a}_{\mathrm{ij} }f\left(x_j,y_j\right)\end{array}$$

(28)

By combining Equations (19), (23), and (28), we can obtain:

$$\begin{array}{c}\frac{E}{1-v^2}\left[\begin{array}{cc}{\mathit{G}}_1^4+b{\mathit{G}}_1^6& a{\mathit{G}}_1^5\\ ⋮& ⋮\\ {\mathit{G}}_{\mathrm{M}}^4+b{\mathit{G}}_{\mathrm{M}}^6& a{\mathit{G}}_M^5\\ a{\mathit{G}}_1^5& b{\mathit{G}}_1^4+{\mathit{G}}_1^6\\ ⋮& ⋮\\ a{\mathit{G}}_M^5& b{\mathit{G}}_M^4+{\mathit{G}}_M^6\end{array}\right]\left[\begin{array}{c}{u}_{x_1}\\ ⋮\\ u_{x_M}\\ u_{y_1}\\ ⋮\\ u_{y_M}\end{array}\right]=-\left[\begin{array}{c}{\overline{f}}_{x_1}\\ ⋮\\ {\overline{f}}_{x_{\mathit{M}}}\\ {\overline{f}}_{y_1}\\ ⋮\\ {\overline{f}}_{y_{\mathit{M}}}\end{array}\right]\end{array}$$

(29)

where ${\mathit{G}}_{\mathit{i}}^n$ denotes the row vector from the $\mathrm{n}$th row of the numerical differentiation operator based within the computational domain, with $X_i$, acting as the point source. In this study, the penalty function method is used to deal with the given displacement and stress boundary conditions. Assuming that a point $X_i$ is a node with a given boundary condition on the boundary, the displacement and stress boundaries are:

$$\begin{array}{c}e\left[\begin{array}{cc}{\mathit{G}}_i^1& 0\end{array}\right]\left[\begin{array}{c}{u}_{x_1}\\ ⋮\\ u_{x_M}\\ u_{y_1}\\ ⋮\\ u_{y_M}\end{array}\right]=e\left[{\overline{u}}_{x_i}\right]\end{array}$$

(30)

$$\begin{array}{c}e\left[\begin{array}{cc}0& {\mathit{G}}_i^1\end{array}\right]\left[\begin{array}{c}{u}_{x_1}\\ ⋮\\ u_{x_M}\\ u_{y_1}\\ ⋮\\ u_{y_M}\end{array}\right]=e\left[{\overline{u}}_{y_i}\right]\end{array}$$

(31)

$$\begin{array}{c}e\frac{E}{1-v^2}\left[n_x\begin{array}{cc}{\mathit{G}}_i^2+n_{y}b{\mathit{G}}_i^3& {n_{y}b\mathit{G}}_i^2+n_{x}v{\mathit{G}}_i^3\end{array}\right]\left[\begin{array}{c}{u}_{x_1}\\ ⋮\\ u_{x_M}\\ u_{y_1}\\ ⋮\\ u_{y_M}\end{array}\right]=e\left[{\overline{t}}_{x_i}\right]\end{array}$$

(32)

$$\begin{array}{c}e\frac{E}{1-v^2}\left[n_{y}v\begin{array}{cc}{\mathit{G}}_i^2+n_{x}b{\mathit{G}}_i^3& {n_{x}b\mathit{G}}_i^2+n_y{\mathit{G}}_i^3\end{array}\right]\left[\begin{array}{c}{u}_{x_1}\\ ⋮\\ u_{x_M}\\ u_{y_1}\\ ⋮\\ u_{y_M}\end{array}\right]=e\left[{\overline{t}}_{y_i}\right]\end{array}$$

(33)

where e is an extremely large penalty coefficient. If the boundary node $X_i$ has given displacements and stress boundary components, adding the boundary conditions determined by Equations (30)–(33) to the corresponding rows of Equation (29) yields:

$$\begin{array}{c}KU=\overline{f}\end{array}$$

(34)

The coefficient matrix $K$ is composed of the control equations for all nodes within the computational domain, together with the boundary equations for the boundary points, multiplied by the penalty coefficients e.

Inverting the system of equations enables the acquisition of the displacement components for each node:

$$\begin{array}{c}U=\left[\begin{array}{c}{u}_{x_1}\\ ⋮\\ u_{x_M}\\ u_{y_1}\\ ⋮\\ u_{y_M}\end{array}\right]=K^{-1}\overline{f}\end{array}$$

(35)

3. Numerical Simulation

In this section, the effectiveness and applicability of the proposed method are verified with numerical simulation, and key characteristics, such as accuracy, computational efficiency, robustness, and convergence, are compared. To accurately evaluate the performance of this method, this study uses two indicators: normalized root mean square error and absolute error.

$$\begin{array}{c}\left\{\begin{array}{c}{\mathit{error}}_u=\sqrt{{\displaystyle \sum _{k=1}^M}{\left(u_{k.\mathit{num}}-u_{k.\mathit{anal}}\right)}^2}/\sqrt{{\displaystyle \sum _{k=1}^M}{\left(u_{k.\mathit{anal}}\right)}^2}\\ {\mathit{error}}_{\sigma }=\sqrt{{\displaystyle \sum _{k=1}^M}{\left({\sigma }_{k.\mathit{num}}-{\sigma }_{k.\mathit{anal}}\right)}^2}/\sqrt{{\displaystyle \sum _{k=1}^M}{\left({\sigma }_{k.\mathit{anal}}\right)}^2}\end{array}\right.\end{array}$$

(36)

$$\begin{array}{c}Absoluteerror=\left|\mathit{Num}-\mathit{Anal}\right|,\end{array}$$

(37)

where ${\mathit{error}}_u$ and ${\mathit{error}}_{\sigma }$ represent the normalized root mean square errors of displacement and stress, respectively. $u_{k.\mathit{num}}$ and ${\sigma }_{k.\mathit{num}}$ are the numerical solutions for displacement and stress, respectively, $u_{k.\mathit{anal}}$ and ${\sigma }_{k.\mathit{anal}}$ are the analytical solutions for displacement and stress, respectively. $\mathit{Num}$ denotes the numerical solution obtained using this method, and $\mathit{Anal}$ represents the analytical solution.

Of note, the numerical simulation results of the following examples are based on the MATLAB R2023b program and were completed using an AMD Ryzen 9 5900HX eight-core processor.

3.1. Cantilever Beam Subjected to Shear Force

In the first example, as illustrated in Figure 2, we consider a cantilever beam with length $L = 48 \mathrm{m}$, width $D = 12 \mathrm{m}$, Poisson’s ratio $v= 0.3$, and elastic modulus $E= 30 \mathrm{MPa}.$ In addition, $I$ represents the moment of inertia, and for a rectangular beam of unit thickness, the moment of inertia $I ={ D}^3/12$. The integral of the distributed force, which is the shear load $P =1000 \mathrm{N}$, is also considered. For plane stress problems, the analytical solutions [33] for displacement and strain are:

$$\begin{array}{c}{u}_x=-\frac{\mathit{Py}}{6\mathit{EI}}\left[\left(6-3x\right)x+\left(2+v\right)\left(y^2-\frac{D^2}{4}\right)\right]\end{array}$$

(38)

$$\begin{array}{c}{u}_y=\frac{P}{6\mathit{EI}}\left[3v{y}^2\left(L-x\right)+\left(4+5v\right)\frac{D^{2}x}{4}+\left(3L-x\right)x^2\right]\end{array}$$

(39)

$$\begin{array}{c}{\sigma }_{\mathit{xx}}=-\frac{P\left(L-x\right)y}{I}\end{array}$$

(40)

$$\begin{array}{c}{\sigma }_{\mathit{yy}}=0\end{array}$$

(41)

$$\begin{array}{c}{\sigma }_{\mathit{xy}}=\frac{P}{2I}\left[\frac{D^2}{4}-y^2\right]\end{array}$$

(42)

3.1.1. Error Analysis

To evaluate the precision and computational speed of this method following adjustments to specific parameters, this study arranged nodes in configurations of 5 × 17, 7 × 25, 13 × 49, 16 × 61, 21 × 81, and 31 × 121 along the x- and y-axes, adopting a third-order Taylor expansion. This study employed computational approaches based on the cubic spline kernel function, quintic spline kernel function, and Gaussian kernel function, comparing these with the computational outcomes derived from the FEM.

As shown in

Table 1. Example 1: Comparison of displacement and stress errors and computational efficiency of each method at different numbers of nodes.

Method		N = 85	N = 175	N = 637	N = 976	N = 1701	N = 3751
Cubic spline	${\mathrm{error}}_{\mathrm{u}}$	1.1559 × 10−6	1.0396 × 10−6	7.4807 × 10−7	6.7813 × 10−7	5.5154 × 10−7	4.1265 × 10−7
	${\mathrm{error}}_{\mathrm{\sigma }}$	1.3159 × 10−6	1.1633 × 10−6	8.9742 × 10−7	7.5859 × 10−7	6.6085 × 10−7	5.4193 × 10−7
	Execution time (s)	0.0221	0.0436	0.1532	0.2605	0.6728	1.6735
Quintic spline	${\mathrm{error}}_{\mathrm{u}}$	1.1043 × 10−6	9.6119 × 10−7	7.0886 × 10−7	5.8304 × 10−7	4.8986 × 10−7	3.7597 × 10−7
	${\mathrm{error}}_{\mathrm{\sigma }}$	1.2215 × 10−6	1.1276 × 10−6	8.1298 × 10−7	6.8245 × 10−7	5.9064 × 10−7	4.8986 × 10−7
	Execution time (s)	0.0223	0.0483	0.1656	0.2906	0.6842	1.6972
Gaussian	${\mathrm{error}}_{\mathrm{u}}$	1.8407 × 10−6	1.6557 × 10−6	1.3924 × 10−6	1.1752 × 10−6	1.0110 × 10−6	9.0412 × 10−7
	${\mathrm{error}}_{\mathrm{\sigma }}$	1.9360 × 10−6	1.7710 × 10−6	1.5148 × 10−6	1.2638 × 10−6	1.0894 × 10−6	1.0310 × 10−6
	Execution time (s)	0.0234	0.0443	0.1352	0.2557	0.6335	1.6421
FEM	${\mathrm{error}}_{\mathrm{u}}$	1.5878 × 10−3	7.1529 × 10−4	1.8235 × 10−4	1.1561 × 10−4	6.3913 × 10−5	1.8732 × 10−5
	${\mathrm{error}}_{\mathrm{\sigma }}$	7.3210 × 10−2	4.3649 × 10−2	1.7055 × 10−2	1.2465 × 10−2	8.2747 × 10−3	2.4296 × 10−3
	Execution time (s)	0.0740	0.1238	0.3882	0.7605	1.4390	3.1675

, compared to the FEM, the proposed method demonstrates higher accuracy and computational efficiency, where the execution time was obtained using MATLAB’s tic-toc function. Moreover, it achieves greater accuracy with the quintic spline and cubic spline kernel functions. Although the scheme using the Gaussian kernel function has slightly lower accuracy, its computation time is marginally less than that of the aforementioned schemes. From a computational cost perspective, this method requires only a small number of nodes to achieve high-accuracy results. Because of the method’s provision of highly accurate differential operators, high-precision spatial derivatives are obtained in the process of deriving stress from displacement, resulting in stress errors very close to displacement errors.

Figure 3 and Figure 4 illustrate a comparison between the numerical solution for the stress components of a cantilever beam obtained with this method and the analytical solution, employing the quintic spline kernel function, with a node quantity N = 5 × 17 = 85 and a Taylor expansion order K = 3. The numerical solution closely matches the analytical solution in both magnitude and distribution.

3.1.2. Robustness Analysis

To further investigate the ability of the proposed method to maintain its predictive accuracy and operational stability when facing various challenges, changes, or external disturbances, this study introduces the perturbation coefficient μ, which is defined using the following formula:

$$\begin{array}{c}\left\{\begin{array}{c}{X}_{\mathit{new}}=X+\mu \times U_{\mathit{random}}\\ \left|{\mathrm{U}}_{\mathit{random}}\right|\le h\end{array}\right.\end{array}$$

(43)

where $X$ represents the coordinates of internal nodes (excluding boundary nodes), $h$ is the smoothing length, $U_{\mathit{random}}$ is a random displacement within the smoothing length, satisfying the condition $\left|U_{\mathit{random}}\right|\le h$, and $X_{\mathit{new}}$ represents the updated node coordinates. In this way, as the perturbation coefficient μ increases, the initial arrangement of nodes becomes more chaotic. Figure 5a–c show the initial arrangement of nodes when μ is set to 0, 0.3, and 0.6, respectively. As the perturbation coefficient μ increases, the initial distribution of nodes becomes more chaotic.

Figure 6 illustrates the relationship between the perturbation coefficient μ and displacement and stress errors. The displacement and stress errors associated with the quintic spline kernel function are most significantly impacted by the increase in node disorder. However, the range of error variation remains within an order of magnitude. In contrast, the impact on accuracy using the other two methods is very limited and almost negligible despite the increase in the perturbation coefficient μ.

From the analysis of the error variation relationship between the perturbation coefficient μ and each kernel function, the corresponding conclusion can be drawn: although the quintic spline kernel function is more sensitive to the influence of the disturbance amplitude, the method still has good robustness as a whole.

3.1.3. Convergence Analysis

This section explores the convergence of the methods used in this study, especially the specific effects of the number of nodes and the Taylor expansion order on the accuracy.

In this study, the initial node layout schemes of 85, 175, 637, 976, and 1701 nodes were adopted. Figure 7 shows a decreasing trend in the displacement error and stress error with the increase in the number of nodes N, and a monotonous negative correlation exists between the two, indicating the method has good convergence and stability. At the same time, as the Taylor expansion order increases from K = 3 to K = 4, the displacement error and stress error obtained with the method of different kernel functions decrease, which also confirms that the method has good convergence.

3.2. Analysis of the Surrounding Rock Compression for a Circular Tunnel

As shown in Figure 8, a circular tunnel surrounded by rock is subjected to uniform external pressure. The inner radius $R_1=1 \mathrm{m}$, the outer radius $R_2=11 \mathrm{m}$, with an elastic modulus $E=2.5 \mathrm{GPa}$ and a Poisson’s ratio of $v=0.3$. The external pressure $P_2=2.5 \mathrm{MPa},$ and the internal pressure $P_1=0 \mathrm{MPa}$. The analytical solutions [34] for displacement and stress are as follows:

$$\begin{array}{c}{u}_r=\left(\frac{1+v}{E}\right)\frac{r}{R_2^2-R_1^2}\left[\left(R_1^2{P}_1-R_2^2{P}_2\right)\left(1-2v\right)-\frac{R_1^2{R}_2^2}{r^2}\left(P_2-P_1\right)\right]\end{array}$$

(44)

$$\begin{array}{c}{u}_{\theta }=0\end{array}$$

(45)

$$\begin{array}{c}{\sigma }_r=\frac{1}{R_2^2-R_1^2}\left[\left(R_1^2{P}_1-R_2^2{P}_2\right)+\frac{R_1^2{R}_2^2}{r^2}\left(P_2-P_1\right)\right]\end{array}$$

(46)

$$\begin{array}{c}{\sigma }_{\theta }=\frac{1}{R_2^2-R_1^2}\left[\left(R_1^2{P}_1-R_2^2{P}_2\right)-\frac{R_1^2{R}_2^2}{r^2}\left(P_2-P_1\right)\right]\end{array}$$

(47)

3.2.1. Error Analysis

In this section, nodes were arranged tangentially and radially in configurations of 11 × 19, 19 × 29, 25 × 39, and 33 × 49. Figure 9 shows the comparison between the numerical solution and the analytical solution when the Taylor expansion order K = 4, using the quintic kernel function. According to

Table 2. Example 2: Comparison of displacement and stress errors and computational efficiency of each method at different node numbers.

Method		N = 209	N = 551	N = 975	N = 1617
Cubic spline	${error}_u$	7.2163 × 10−3	5.7165 × 10−3	4.5403 × 10−3	3.4731 × 10−3
	${error}_{\sigma }$	8.4958 × 10−3	6.6441 × 10−3	5.3892 × 10−3	4.2186 × 10−3
	Execution time (s)	0.0912	0.1732	0.3984	0.8547
Quintic spline	${error}_u$	6.1713 × 10−3	4.5487 × 10−3	3.1723 × 10−3	1.9475 × 10−3
	${error}_{\sigma }$	7.2931 × 10−3	5.5299 × 10−3	3.7452 × 10−3	2.4782 × 10−3
	Execution time (s)	0.1031	0.1973	0.4134	0.8835
Gaussian	${error}_u$	9.9101 × 10−3	8.6321 × 10−3	7.3814 × 10−3	6.1756 × 10−3
	${error}_{\sigma }$	1. 0927 × 10−2	9.2074 × 10−3	7.6133 × 10−3	6.2943 × 10−3
	Execution time (s)	0.1074	0.1867	0.4007	0.8379
FEM	${error}_u$	4.8610 × 10−4	2.5399 × 10−4	1.6901 × 10−4	1.1947 × 10−4
	${error}_{\sigma }$	8.0192 × 10−3	5.7098 × 10−3	4.4219 × 10−3	3.7964 × 10−3
	Execution time (s)	0.3944	0.7821	1.6741	3.5971

, the cubic spline kernel function, quintic spline function, and Gaussian kernel function all exhibit high accuracy, with displacement and stress errors on the order of 10⁻³. Although the accuracy is slightly lower than the FEM, the computation time is significantly less than that of the FEM. Additionally, compared to the FEM, the difference between displacement and stress errors in this method is minimal. This is attributed to the method’s differential operator providing high-accuracy partial derivatives, resulting in smaller errors in the stress components derived from displacement components.

3.2.2. Robustness Analysis

This section utilizes both regular and irregular nodes (with a perturbation coefficient μ $=$ 0.2), as shown in Figure 10. Given the symmetry of this example, the robustness analysis in this case only requires comparing the displacement and stress components along the x-axis. Figure 11 and Figure 12 display the numerical solutions for displacement and stress under regular and irregular initial node distributions, respectively, along with the absolute errors compared to the analytical solutions.

As shown in Figure 11, when employing a regular arrangement of initial nodes, the maximum absolute error of displacement along the x-axis only reaches the level of 2 × 10⁻⁵, and it is confined to a small portion within the inner radius range. When using an irregular node arrangement, the displacement along the x-axis increases progressively, with a maximum of 8 × 10⁻⁵. Although this is an increase compared to the regular node distribution, it still significantly differs from numerical solutions at the 10⁻³ level. Figure 12 shows that although the area of the absolute error of the stress becomes larger when the nodes are arranged irregularly, the maximum value decreases. From the perspective of displacement and stress errors, in this example, this method shows a low sensitivity to the confusion of the initial node arrangement.

3.2.3. Convergence Analysis

This section analyzes the displacement and stress errors when using different kernel functions and orders of Taylor expansions with node quantities of 209, 551, 975, and 1617. As shown in Figure 13, with an increase in the number of nodes, the precision of the calculations improves. In conjunction with Example 1 in Figure 7, the methods employing the three types of kernel functions all exhibit convergence, and the rate of decrease in errors gradually slows as the number of nodes increases. Additionally, using a higher order of Taylor expansion significantly enhances the accuracy of the method, bringing it closer to the analytical solution.

3.3. Homogeneous Slope Subjected to Self-Weight

Next, we consider a homogeneous slope solely subjected to self-weight [35], and Figure 14 displays its geometric model. The material properties are as follows: elastic modulus E = 80 MPa, unit weight γ = 19.62 kN/m³, and Poisson’s ratio v = 0.43. The model’s base is fully fixed, with normal constraints applied to its sides.

This research deployed the FEM to set up a dense network of nodes and elements (86,372 nodes and 85,650 quadrilateral elements) to derive a reliable reference solution for comparison with the results from the GFDM.

Figure 15 shows the computational outcomes of this approach versus the FEM, while

Table 3. Maximum values of displacement and stress components for the GFDM and FEM.

	Method	Max
Ux (m)	GFDM	0.425927
	FEM	0.425957
Uy (m)	GFDM	6.30321
	FEM	6.30359
σx (Pa)	GFDM	2.10600 × 106
	FEM	2.09875 × 106
σy (Pa)	GFDM	4.89859 × 106
	FEM	4.88162 × 106
τxy (Pa)	GFDM	2.07472 × 105
	FEM	2.07541 × 105

presents the maximum values for displacement and stress components from this method alongside the reference solution. Figure 15 and Table 3 demonstrate that the magnitude and distribution of results from this method closely match the reference solution obtained with the FEM.

4. Conclusions

This article introduces a meshless method that leverages the strengths of the PDDO and GFDM. This approach retains the nonlocal characteristics of the PDDO while reducing many preliminary assumptions. Unlike traditional meshless methods, it does not require special treatment of boundary conditions, making it easier to implement and more computationally efficient. Numerical simulations of several classic examples compared with traditional numerical methods (FEM) verify its superior performance in terms of accuracy, efficiency, robustness, and convergence. In particular, it demonstrates high adaptability and robustness in scenarios of irregular node arrangements, increased perturbation coefficients, and more nodes, confirming its potential in solving complex engineering problems.

In conclusion, the meshless numerical simulation method combining the PDDO and GFDM proposed in this study provides a new solution in the field of numerical analysis of elastic statics. Future work will focus on further expanding the application scope of this method, including exploring its effectiveness in dealing with more complex multi-physics problems and improving numerical algorithms to achieve higher accuracy and computational efficiency.