Extended Peridynamics and Parameter Optimization Study Based on Moving Least Squares Method

Based on the theory of peridynamics, the least squares and the moving least squares method are proposed to ﬁt the physical information at nondiscrete points. It makes up for the shortcomings of the peridynamic method that only solves the discrete nodes and cannot obtain the physical information of other blank areas. The extended method is used to ﬁt the one-way vibration problem of the rod, and the curve of the displacement of a nondiscrete node in the rod is extracted with time. The ﬁtted displacement results are compared with the theoretical results to verify the feasibility of the ﬁtting method. At the same time, the parameters in the ﬁtting of the moving least squares method are optimized, and the eﬀects of diﬀerent tight weight functions and inﬂuence ranges on the results are analyzed. The results show that when the weight function is a power exponential function, the ﬁtting eﬀect increases with the decrease in the coeﬃcient. When the weight function is a cubic spline weight function, a better ﬁtting eﬀect is obtained. And in the case of ensuring the ﬁtting result, the aﬀected area should be reduced as much as possible, and the calculation eﬃciency and precision can be improved.


Introduction
e failure of solid materials and structures is a classical problem in mechanical research, and it is also the focus of attention in mechanical, civil, water conservancy, chemical, and other engineering fields. Under different conditions, the cracks in the engineering materials will expand and eventually form macrocracks, resulting in structural strength reduction or even damage, which brings serious safety problems. However, due to the existence of microcracks with various geometry sizes and heterogeneous spatial distributions, the theoretical solutions of crack growth are usually not available, and alternatively, the numerical simulation is commonly applied to replicate the fracturing scenario.
Although numerical simulation methods bring some convenience and can simulate the damage [1,2], the traditional numerical methods have some drawbacks in terms of simulating material damage; for example, the finite element method encounters the problem of a crack tip singularity during simulation, and it is necessary to introduce external criteria when using the finite element method to simulate the problem of material crack propagation. An extended finite element theory [3] is proposed for dealing with the above problems, but it is still necessary to introduce external criteria to deal with the discontinuous displacement; furthermore, it is impossible to simulate complex impact damage problems. Discrete element methods [4,5] based on discontinuous assumptions have achieved satisfactory results in simulating materials, but they cannot solve the characteristics of continuous change. Molecular dynamics that can describe microscopic characteristics are used to simulate material damage and crack propagation, but the methods usually require a long computation time [6].
Recently, a peridynamic method [7] based on nonlocal thoughts has been developed. In the method, it is not necessary to differentiate the displacement. By solving the integral equation that describes the constitutive behavior of the material, the problem of the tip singularity in the finite element method can be solved [8]. In addition, it has its own failure criterion in the constitutive relationship of peridynamics; thus, there is no need to introduce external laws to govern the evolution of structural damage. Peridynamics has the characteristics of a meshless method and molecular dynamics and solves problems for which the molecular dynamics calculation efficiency is too low. Huang et al. [9] proposed an extended peridynamic approach for quantitative analysis of elastic deformation and simulation of both quasistatic and dynamic crack propagation. Gerstle et al. [10] proposed the peridynamic model by adding pairwise peridynamic moments to simulate linear elastic materials with varying Poisson's ratios. Zhou et al. [11] introduced the theory of the peridynamics to the failure of rock materials. Ma and Li [12] used the peridynamics to determine energy absorption characteristics of ordinary glass under impact load.
When using the peridynamics to model and calculate, it is necessary to divide the model into a series of nodes with physical information. e resolution of discretization is usually limited by the computation resources, in which only finite nodes are generated; thus, the physical information at nonnodes cannot be obtained. e meshless method [13][14][15][16][17] does not need to generate mesh in numerical calculation but constructs the interpolation function discrete control equation according to some arbitrarily distributed coordinate points, which can simulate various complex physical fields (such as displacement field and temperature field) conveniently. e least square method has always been a common method for curve and surface fitting. e moving least square method can improve the calculation efficiency and accuracy by introducing the tight support function and the base function, and the meshless method can be used for fitting [18,19]. e rest of the paper is organized as follows: in Section 2, we introduce the basic theory of the peridynamics including the modeling and calculation process. In Section 3, the physical information at the discrete points is fitted by using the least square method and moving least square method. In Section 4, the fitting solution is compared with the theoretical solution to verify the effectiveness of the method. In Section 5, the influence of the weight function and the radius of the influence area on the results are discussed, respectively. Finally, the conclusions are drawn in Section 6.

Basic eory.
Peridynamics belongs to the category of nonlocal methods. Compared with other nonlocal methods, what is most characteristic of peridynamics is the use of spatial integral formula instead of differential formula. is can solve the problem that the differential formula cannot be calculated when the space is discontinuous. Compared with traditional continuous mechanics, the basic idea of peridynamics is to interact with other material points in a certain region. e space R is discretized into point elements, as shown in Figure 1, and the basic equation of a material point x at time t is where x′ is the point around the material point x, ρ is the density, f is the peridynamic force, and b is the body force.
Here, we define that a point within a distance from the material point x will interact with its center point x. us, points that interact with x at a certain time t can be represented by a set: e displacements of the material points x and x′ at time t are u(x, t) and u(x ′ , t), respectively. e relative position vector of the two material points at the initial moment is denoted as x ′ − x � ξ, and the relative displacement vector generated between the two points after the elapsed time t is erefore, the relative position vector after the deformation occurs is ξ + η. In addition, ρ(x) and € u(x, t) on the left side of the basic equation are the density and acceleration of the material point x at time t, respectively, and on the right side of the basic equation, b(x, t) is the external force received by the material point x, and it can be used to apply boundary conditions. e function f is the interaction force between the points x and x′, and the material point x is coordinated by other points in the region to maintain the motion or equilibrium state. We term the interaction between point x and x′ as "bond interaction" and idealized the behavior as the elastic force of the spring. Because, in traditional continuous mechanics, the interaction between units is based on the force of direct contact, the concept of a bond at a finite distance is a fundamental difference between the theory of peridynamics and the theory of traditional continuity.

Constitutive Model.
In microelastic materials, the local interaction force function f between discrete points can be derived from the micropotential energy: e micropotential energy here is a scalar function that characterizes the interaction relationship between material points, specifically expressed as the energy per square of the unit volume stored in a single bond. us, the strain energy density at one point can be expressed as e 1/2 term in equation (4) indicates that the strain energy stored in the substance point x has half of the energy of the interaction bond. According to the theory of bondbased peridynamics, the direction of the interaction force between the material points is consistent with the direction of deformation. So, the constitutive force f can be expressed as 2 Mathematical Problems in Engineering e deformation of the bond between the material points is represented by the elongation s. e elongation not only characterizes the mechanical relationship between the points but also controls the breakage of the bond. It is expressed as e elongation value can be positive or negative, representing the stretching or compression of the bond, respectively. When the bond tensile deformation exceeds the critical elongation s 0 , the bond will undergo an irreversible break.
e specific derivation of the linear microelastic model can be found in the work by Silling and Askari [20], in which the expression of the micropotential energy function is Failure happens when the deformation is beyond a critical value, s 0 , i.e., the so-called critical relative elongation, and the bonds break. s 0 is computed based on the fracture energy of a material. e energy per unit fracture length for complete separation of the two halves of the body is the fracture energy G 0 : In the above formula, c is the microelastic modulus, which is similar to the elastic modulus in the traditional continuous theory and characterizes the deformation characteristics of the material. According to equations (3) and (5), the expression of the constitutive force function f can be obtained as where μ is a scalar function that controls the breakage of the bond. Its specific form is as follows: A sudden failure occurs when the stress reaches a peak. is is similar to the damage characteristics of a brittle material. erefore, peridynamics does not need to introduce additional criteria to determine whether the material is damaged, and it defines the damage at a point by considering the breakage of the bond at a point: A D value of 0 indicates that the damage has not occurred at the point x, and a D value of 1 indicates that the point is destroyed all. e numerical implementation method for peridynamics can be referred to in [11].

Discretization.
In order to solve the peridynamic model, it is necessary to discretize the whole model area with uniformly distributed particles. e physical properties of each element are transferred through the geometric center point (or other points) of the element, as shown in Figure 2.
e other physical parameters such as external load also act on the node, so the information at the node is obtained when using the mobilization dynamics calculation: where n is the number of time steps and V is the volume of points in the neighborhood. For the material points on the neighborhood boundary that are not in the region, the corresponding volume needs to be corrected and calculated. It can be reduced in proportion according to the relationship between the material points on the neighborhood boundary and the neighborhood radius according to [20]. e time step discrete scheme uses the Verlet velocity difference scheme.

Least Square Fitting.
As illustrated in the previous section, the peridynamic method often employs uniformly distributed particles to disperse the material area of the whole model in numerical calculation, and the physical properties of each element are transmitted through the geometric center of the element. However, after the numerical calculation, we can only get the physical information (such as displacement and temperature) at the center of the discretization. e physical information at a certain point between the discrete points cannot be obtained. In order to obtain the physical information at a certain point, we can increase the discrete points, but we cannot discretize infinitely small in the discretization model, which will lead to the increase in the calculation cost. Instead, we can apply the physical information of the known points around to fit the information of the nondiscrete point area. Let there be a series of physical information (such as displacement and temperature) of the peridynamics discrete calculation point coordinate In order to make the polynomial most approximate to the known discrete point, we make the sum of the squares of the difference between the fitting solution of the polynomial and the known solution to the minimum value: If S obtains the minimum value, then the partial derivative of equation (13) is 0, so in order to obtain the value of each coefficient a k (k � 0, 1 . . . k) in the k-order polynomial, the partial derivative of right a k (k � 0, 1, . . . , k) of equation (11) is obtained, respectively, and then K + 1 equation is obtained: By combining equation (14), the solution can be read: e above formula can be written in the form of matrix: where erefore, we can get the global approximate function by knowing the physical information of the discrete points (such as displacement and temperature) and then fit the physical information at any point. Because the calculation efficiency of the peridynamics is lower than that of other methods, we can discrete fewer the peridynamics calculation points in space under the condition of ensuring the accuracy, and the rest areas can be obtained by fitting equations to improve the efficiency of calculation.

Moving Least Square
Fitting. When using the least square method to fit, if the degree of the trial function polynomial is high, the ill conditioned equations are easy to be generated in the process of fitting, which will affect the fitting results. At this time, we can divide the region of discrete points into several subregions. Although this division can avoid the occurrence of oscillation, it also affects the continuity and smoothness of curves or surfaces. In order to solve the problem of ill conditioned equations and the influence of smoothness, two new concepts are introduced in the MLS: tight support weight function and base function.
e concept of the tight support weight function is similar to that of the shadow region in the peridynamics, which means that the physical information at a nondiscrete point is only affected by the discrete points in the neighborhood, and it also reflects the interaction in a certain region with the peridynamics. e physical information (such as displacement and damage value) x i (i � 1, . . . , N) at the peridynamic nodes u i � u(x i ) in the region is known. Our purpose is to construct a global approximation function u h (x) of physical information in the whole solution domain, and then, the function to be solved can be locally approximated as follows: where x is the fitting point, x is the known discrete point in the neighborhood of the point to be fitted, e tight support weight function is consistent with the scalar control function μ in the peridynamics. e solution domain is discretized in the peridynamics nodes, and a weight function w i (x) � w(x − x i ) is defined at each node. e weight function w i is greater than 0 in the limited region around node x i and 0 outside the region. is region is the support region of weight function w i . It is similar to the effect of neighborhood radius δ in the peridynamics, and it only works in the region. e introduction of weight function w i can solve the problem of segment fitting discontinuity, which embodies the characteristics of fitting locality. e selection of coefficient a i (x) enables trial function u h (x) to obtain the best fitting result in its solution domain. e weighted square sum of the errors of u h (x) at all discrete points in the solution domain is e above formula can be expressed in the form of matrix: where u � u 1 , u 2 , . . . , u N , Take the extreme value for J to get the coefficient a(x): where A(x) and B(x) are written as Furthermore, we can get erefore, the expression of trial function u h (x) is where Φ(x) is a shape function: Compared with the least square method, the major difference between the moving least square method and the least square method is that the coefficient vector and the basis function are used to build the fitting function. e coefficient vector is a function of system coordinate, and the concept of tight support function is introduced. ese improvements can reduce the influence of continuity and smoothness in curve and surface fitting and only interact in neighborhood like the peridynamics.

Weight Function and Influence Area.
e introduction of weight function is of great significance. When fitting number of discrete points of the peridynamics, it can be localized. Several principles need to be followed when selecting the fitting function:

Mathematical Problems in Engineering
(1) e weight function needs to be smooth, continuous, and differentiable. (2) Unit decomposition.
(3) It is not 0 in the influence area of weight function but 0 outside the area. (4) e weight in the influence area of the weight function is inversely proportional to the distance from the center point.
e commonly used weight functions are power index weight function, multiple spline weight function, and so on. e power index weight function is as follows: where d � d/d max is the relative distance, d is the distance from the calculation point to the node, and d max is the influence radius of the node, that is, the range of each weight function. e three-power spline weight function is shown in the following formula: e four-power spline weight function is shown in the following formula: When the moving least square is used for fitting, the range of the weight function at each discrete point should be determined, which can be consistent with the selection of the peridynamic neighborhood. e optimization of action radius will be introduced in the following section.

Verification Scheme.
In this section, the transverse vibration of members is taken as an example. First, the peridynamics method is used to model and calculate the displacement of each node with time. en, the least square and the moving least square are used to fit the displacement at the nondiscrete point. e fitting solution is compared with the theoretical solution to verify the effectiveness of the method.
It is assumed that the member will produce strain under the initial tensile load, and the transverse vibration will occur with time after the load is removed, as shown in Figure 3. e theoretical solution of the displacement of any point on the member changing with time after load removal is as follows: e peridynamics is used for discretization and calculation, and the discretization method is in accordance with the abovementioned and [11]. 1000 particles are dispersed in the direction of vibration, and 3 particles are dispersed at the fixed position at the left end. e particle occupied volume is ΔV � 1 × 10 − 9 m 3 , the particle spacing is Δ � 0.001 m, the neighborhood radius is δ � 3.015Δ, and the time step is Δt � 1.95 × 10 − 7 s. e volume correction coefficient is as mentioned above, the total time steps are 25000, and the calculation process is as shown above. FORTRAN language is used to write and calculate. e displacement of each particle with time is extracted by using the results of the peridynamics. In order to improve the calculation speed and efficiency, the least square and moving least square are used to fit the displacement of nondiscrete points, and the fitting results are compared with the theoretical solution. e general higher-order polynomials are selected in least square fitting, while the threepower spline weight function is temporarily selected in compact support function fitting by the moving least square method. e selection of influence region is consistent with the selection of neighborhood radius of the peridynamics.

Results.
e nondiscrete points around the coordinate x � 0.4995 of the discrete points of the peridynamics are selected for fitting. In order to verify the generality, x � 0.4991 and x � 0.4997 are selected for fitting. e comparison between the fitting results and the theoretical results at x � 0.4997 is shown in Figure 4. e comparison shows that the fitting results of the least square method and the moving least square method are close to the theoretical solution at x � 0.4991, which shows that it is feasible to fit the physical information at nondiscrete points. e comparison of the ML and the MLS fitting results is shown in Figure 5. We find that the fitting results of the LS and the MLS are very close, and they can fit the displacement of nondiscrete points well. We can get the following conclusion: when the data are less, both methods can be used.  When the data are more and change greatly, the MLS is recommended, which can better avoid the emergence of ill conditioned equations.
In this simulation, the CPU running time of LS and MLS for one point is 1.1090659 s and 1.1363598 s, respectively, and the CPU model is Intel Core i7-7700hq. We found that the running time of LS method is less in this test.
In order to compare with other numerical methods, the results of ANSYS transient analysis using FEM are also shown in Figure 5. We find that the results of several different numerical calculations are close to the theoretical solution.

Further Verification.
In order to further verify the effectiveness of the method, a two-dimensional problem, that is, a rectangular isotropic plate subjected to uniaxial uniform tension, is tested in this section. e tension is applied in the boundary layer region by body force concept in the peridynamics, as shown in Figure 6. e length of the plate is 1 m with a width of 0.5 m, elastic modulus of E � 200 GPa, Poisson's ratio of ] � 1/3, and density of ρ � 7850kg/m 3 . e values of initial tensile stress are 200 MPa. e discretization and calculation methods are in accordance with the abovementioned. 5000 particles are dispersed, and the total time steps are 3000. e theoretical solution of the displacement of any point on the member changing with location is e nondiscrete point at the coordinate x � 0.253 and y � 0.123 is fitted by the discrete points in the same way. Because the results of the two fitting methods in the previous example are almost the same, we only use the MLS here.
e comparison between the fitting results and the theoretical results is shown in Figure 7. e comparison shows that the fitting results of the moving least square method are close to the theoretical solution, which further shows that it is feasible to fit the physical information at nondiscrete points.

Optimization of Tight Supported Weight Function.
Based on the discrete data of the peridynamics, the influence of power index weight function of different coefficients and multiple spline weight function on the fitting results is analyzed under the same basis function.
In order to facilitate the comparison of fitting results, a series of near-field dynamic discrete point displacements are set to satisfy y � [(x 2 − 1) 3 + 0.5]sin(2x). e series of discrete points are fitted by power index weight function and spline weight function, respectively. e fitting curve is compared with a series of points to evaluate the fitting effect. e comparison between the fitting curve of power index weight function with different coefficients and the original discrete point is shown in Figure 8. e comparison between the fitting curve with three-power and four-power spline weight function and the original discrete point is shown in Figure 7.
It can be seen from the comparison of Figure 8 that the power index weight functions with different coefficients have a great influence on the fitting results. When the coefficient is large, the fitting effect in the data boundary area will produce some fluctuations, and with the reduction in the coefficient, the fluctuations will also weaken, but the calculation cost will also increase.
Compared with Figure 9, the fitting results of threepower and four-power spline weight functions are close to the original discrete points. e fitting effect of three-power spline weight function is slightly higher than that of quartic spline weight function. e fitting effect of four-power spline weight function will deviate to a certain extent when the original discrete point data change greatly.
To sum up, the power index weight function with a smaller selected coefficient and the three-power spline weight function can achieve better fitting results under the condition of not significantly increasing the calculation cost.

Optimization of Influence Range.
e radius of neighborhood is generally 3 times the discrete distance of model particles, but there is no uniform method to select the influence range of weight function. erefore, it is necessary to analyze the influence of weight function range on the results. e selection of discrete data is the same as the previous section. e influence range of weight function is 0.5, 1, and 1.5 times of neighborhood radius, and the same base function and weight function are used for curve fitting. We compare the fitting curve with the discrete points to evaluate the fitting effect. e comparison figure is shown in Figure 10.
Compared with the data in Figure 10, different influence ranges have great influence on the fitting results. With the increase in the influence range, there is a certain degree of deviation in the original discrete point where the data change greatly, and the larger influence area will also increase the calculation cost, so the smaller influence area is selected when selecting the influence area.

Conclusions
In the paper, we use the least square method and moving least square method to fit the physical information in the peridynamics. e conclusions are as follows: (1) In the calculation of the peridynamics, the whole region cannot be discretized infinitely. In order to obtain the physical information of nondiscrete points, the least square method and the moving least square method can be used for fitting solution under the condition of ensuring the accuracy. e moving least square method can reduce the occurrence of ill conditioned equation and ensure the continuity and smoothness of the fitting results. (2) Choosing different weight functions has great influence on fitting results. With the increase in power index weight function coefficient, the fitting results will generate data fluctuation near the edge. e fitting effect of the four-power spline weight function will deviate to a certain extent when the data changes greatly, and the fitting effect of the three-power spline weight function is better.
(3) Selecting different weight function influence ranges will affect the fitting results. With the expansion of the influence ranges, the calculation cost will increase, and the fitting effect will be biased, so the influence ranges of the weight function will be reduced under the condition that the fitting result can be obtained.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.