Abstract
Based on the theory of peridynamics, the least squares and the moving least squares method are proposed to fit 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 fit 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 fitted displacement results are compared with the theoretical results to verify the feasibility of the fitting method. At the same time, the parameters in the fitting of the moving least squares method are optimized, and the effects of different tight weight functions and influence ranges on the results are analyzed. The results show that when the weight function is a power exponential function, the fitting effect increases with the decrease in the coefficient. When the weight function is a cubic spline weight function, a better fitting effect is obtained. And in the case of ensuring the fitting result, the affected area should be reduced as much as possible, and the calculation efficiency and precision can be improved.
1. Introduction
The 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. The 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. The meshless method [13–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. The least square method has always been a common method for curve and surface fitting. The 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].
The 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.
2. Theory of the Peridynamics
2.1. Basic Theory
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. This 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. The 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. Thus, points that interact with x at a certain time t can be represented by a set:
The displacements of the material points x and x′ at time t are and , respectively. The relative position vector of the two material points at the initial moment is denoted as , and the relative displacement vector generated between the two points after the elapsed time t is . Therefore, the relative position vector after the deformation occurs is . In addition, and 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, is the external force received by the material point x, and it can be used to apply boundary conditions. The 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.
2.2. Constitutive Model
In microelastic materials, the local interaction force function f between discrete points can be derived from the micropotential energy:
The 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. Thus, the strain energy density at one point can be expressed as
The 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 bond-based 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
The deformation of the bond between the material points is represented by the elongation s. The elongation not only characterizes the mechanical relationship between the points but also controls the breakage of the bond. It is expressed as
The 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 , the bond will undergo an irreversible break. The 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, , i.e., the so-called critical relative elongation, and the bonds break. is computed based on the fracture energy of a material. The energy per unit fracture length for complete separation of the two halves of the body is the fracture energy :
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. This is similar to the damage characteristics of a brittle material. Therefore, 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. The numerical implementation method for peridynamics can be referred to in [11].
2.3. Discretization
In order to solve the peridynamic model, it is necessary to discretize the whole model area with uniformly distributed particles. The physical properties of each element are transferred through the geometric center point (or other points) of the element, as shown in Figure 2. The 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]. The time step discrete scheme uses the Verlet velocity difference scheme.
3. Moving Least Square Fitting
3.1. 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. The 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 the space as . We can construct the global approximation function about x, that is, the k-order polynomial about x. 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 in the k-order polynomial, the partial derivative of right of equation (11) is obtained, respectively, and then equation is obtained:
By combining equation (14), the solution can be read:
The above formula can be written in the form of matrix: where
Therefore, 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.
3.2. 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. The 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.
The physical information (such as displacement and damage value) at the peridynamic nodes in the region is known. Our purpose is to construct a global approximation function of physical information in the whole solution domain, and then, the function to be solved can be locally approximated as follows: where is the fitting point, is the known discrete point in the neighborhood of the point to be fitted, , is the basic function, m is the number of basic functions, , and is the undetermined coefficient. In the paper, the one-dimensional problem is taken as an example, as follows:(i)Linear basis: (ii)Square basis: (iii)Cubic basis:
The tight support weight function is consistent with the scalar control function in the peridynamics. The solution domain is discretized in the peridynamics nodes, and a weight function is defined at each node. The weight function is greater than 0 in the limited region around node and 0 outside the region. This region is the support region of weight function . It is similar to the effect of neighborhood radius in the peridynamics, and it only works in the region. The introduction of weight function can solve the problem of segment fitting discontinuity, which embodies the characteristics of fitting locality.
The selection of coefficient enables trial function to obtain the best fitting result in its solution domain. The weighted square sum of the errors of at all discrete points in the solution domain is
The above formula can be expressed in the form of matrix: where
Take the extreme value for J to get the coefficient :where and are written as
Furthermore, we can get
Therefore, the expression of trial function iswhere 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. The coefficient vector is a function of system coordinate, and the concept of tight support function is introduced. These improvements can reduce the influence of continuity and smoothness in curve and surface fitting and only interact in neighborhood like the peridynamics.
3.3. Weight Function and Influence Area
The 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: (1)The 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)The weight in the influence area of the weight function is inversely proportional to the distance from the center point.
The commonly used weight functions are power index weight function, multiple spline weight function, and so on. The power index weight function is as follows:where is the relative distance, is the distance from the calculation point to the node, and is the influence radius of the node, that is, the range of each weight function.
The three-power spline weight function is shown in the following formula:
The 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. The optimization of action radius will be introduced in the following section.
4. Numerical Verification
4.1. 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. Then, the least square and the moving least square are used to fit the displacement at the nondiscrete point. The 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. The left end of the member is a fixed end, with a length of 1 m and a cross section of 0.001 m rectangle, elastic modulus of , Poisson’s ratio of , and density of . The values of initial tensile strain and initial velocity are both 0. The theoretical solution of the displacement of any point on the member changing with time after load removal is as follows:
The 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. The particle occupied volume is , the particle spacing is , the neighborhood radius is , and the time step is . The 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.
The 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. The general higher-order polynomials are selected in least square fitting, while the three-power spline weight function is temporarily selected in compact support function fitting by the moving least square method. The selection of influence region is consistent with the selection of neighborhood radius of the peridynamics.
4.2. Results
The nondiscrete points around the coordinate of the discrete points of the peridynamics are selected for fitting. In order to verify the generality, and are selected for fitting. The comparison between the fitting results and the theoretical results at is shown in Figure 4.
(a)
(b)
The comparison shows that the fitting results of the least square method and the moving least square method are close to the theoretical solution at , which shows that it is feasible to fit the physical information at nondiscrete points.
The 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.
4.3. 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. The tension is applied in the boundary layer region by body force concept in the peridynamics, as shown in Figure 6. The length of the plate is 1 m with a width of 0.5 m, elastic modulus of , Poisson’s ratio of , and density of . The values of initial tensile stress are .
The discretization and calculation methods are in accordance with the abovementioned. 5000 particles are dispersed, and the total time steps are 3000. The theoretical solution of the displacement of any point on the member changing with location is
The nondiscrete point at the coordinate 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. The comparison between the fitting results and the theoretical results is shown in Figure 7.
The 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.
5. Parameter Optimization
5.1. 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 . The series of discrete points are fitted by power index weight function and spline weight function, respectively. The fitting curve is compared with a series of points to evaluate the fitting effect. The comparison between the fitting curve of power index weight function with different coefficients and the original discrete point is shown in Figure 8. The comparison between the fitting curve with three-power and four-power spline weight function and the original discrete point is shown in Figure 7.
(a)
(b)
(c)
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 three-power and four-power spline weight functions are close to the original discrete points. The fitting effect of three-power spline weight function is slightly higher than that of quartic spline weight function. The fitting effect of four-power spline weight function will deviate to a certain extent when the original discrete point data change greatly.
(a)
(b)
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.
5.2. Optimization of Influence Range
The 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. Therefore, it is necessary to analyze the influence of weight function range on the results.
The selection of discrete data is the same as the previous section. The 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. The comparison figure is shown in Figure 10.
(a)
(b)
(c)
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.
6. Conclusions
In the paper, we use the least square method and moving least square method to fit the physical information in the peridynamics. The 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. The 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. The 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
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported partially by the National Natural Science Foundation of China (nos. 51879150 and 41831278), Science and Technology Project of Qilu Transportation Development Group (no. 2016B20), and Major Scientific and Technological Innovation Projects in Shandong Province (no. 2019JZZY010428).