Analysis of Elastic–Plastic Problems Using the Improved Interpolating Complex Variable Element Free Galerkin Method

: A numerical model for the two-dimensional nonlinear elastic–plastic problem is proposed based on the improved interpolating complex variable element free Galerkin (IICVEFG) method and the incremental tangent stiffness matrix method. The viability of the proposed model is veriﬁed through three elastic–plastic examples. The numerical analyses show that the IICVEFG method has good convergence. The solutions using the IICVEFG method are consistent with the solutions obtained from the ﬁnite element method using the ABAQUS program. Moreover, the IICVEFG method shows greater computing precision and efﬁciency than the non-interpolating meshless methods.


Introduction
In the science and engineering fields, nonlinear elastic-plastic problems are very common, and are complicated due to the nonlinearity of the plastic state. In addition, the corresponding analytical solutions are difficult to be obtain [1][2][3]. Thus, it is important to study efficient and accurate computational models for solving nonlinear elastic-plastic problems, particularly those problems characterised by geometric [4][5][6] and material nonlinearities [7][8][9]. In recent years, the meshless method has played a great role in numerical simulation. Unlike conventional mesh-based numerical methods [10,11], the meshless method is built on discrete nodes, which does not cause mesh distortion. Therefore, the meshless method is more suitable for nonlinear large deformation and fracture problems, as well as complicated porous structures [12][13][14][15][16][17].
As a popular numerical computing method, the meshless method has developed various branches, such as the smoothed particle hydrodynamics method [18], the famous moving least squares (MLS) approximation [11], the reproducing kernel particle method [19], and so on. In this paper, we mainly study methods that improve upon the original MLS approximation in order to overcome shortcomings which lead to low efficiency and low accuracy [20][21][22]. The shortcomings include relying on numerous nodes, possessing a non-interpolating shape function, and possibly producing an ill-conditioned final discrete system.
To overcome inefficiency, complex variables were introduced into the MLS approximation, and the complex variable moving least squares (CVMLS) approximation was proposed [20]. Since the CVMLS approximation can use a complex variable to express the two directional variables, the CVMLS approximation has a higher computing efficiency [21,22]. However, this approximation lacks specific mathematical implications. Thus, the improved complex variable moving least squares (ICVMLS) approximation was proposed, which brings in a conjugated basis function and a new functional [23][24][25].

Brief Descriptions of the Two-Dimensional Elastic-Plastic Problem
Within the problem domain Ω, the classic equilibrium equation of a two-dimensional elastic-plastic problem is [34] L T .
For an arbitrary point z = x 1 + ix 2 ∈ Ω, b 2 ) is the body force rate field, and L is the differential operator matrix, The relationship between the strain rate field u 2 ) is described as and the stress-strain relationship is ε. (4) Considering that this problem is a plane stress problem, the constitutive matrices D are different in the elastic and plastic state, which are shown in the following, respectively [35].
where E is the elastic modulus, ν is the Poisson's ratio, and σ is the equivalent stress Assume that the stress-strain relationship in the plastic state satisfies the linear hardening model, the plastic modulus for hardening H is expressed as where E is the tangent modulus. In addition, σ ij is the stress deviation which can be written as The velocity and natural boundary conditions are where u 2 ) is the prescribed velocity vector on the velocity boundary Γ u , t 2 ) is the given traction rate, and n is the unit vector of the outer normal direction on the natural boundary Γ t .

The IICVEFG Method for the Elastic-Plastic Problem
In this paper, the IICVMLS method is applied to disperse the integral weak form of the elastic-plastic problem, and then the corresponding IICVEFG method is presented [30]. According to Equations (1) and (14), the integral weak form is expressed as In the IICVMLS method, the trial function of velocity where .
Then the velocity matrix . u(z) in Equation (15) can be written as where . .
According to Equation (27), the strain rate of point z in Equation (3) can be rewritten and simplified as . where Then the stress rate of point z in Equation (4) is expressed as Substitute Equations (27), (32) and (35) into Equation (15), then consider the arbitrary of δ . U(z) T . Finally, we obtain the discrete matrix equation Mathematics 2021, 9, 1967 5 of 16 where The shape function of Equation (17) in the IICVMLS method has the interpolating feature, which results in the fact that Equation (13) can be enforced into Equation (36) directly. Therefore, the final matrix Equation (36) is more succinct than that in the non-interpolating complex variable meshless method. Thus, the IICVEFG method can, theoretically, solve the problem with greater precision and efficiency.

Incremental Tangent Stiffness Matrix Method
In this part, the incremental tangent stiffness matrix method is applied to solve the nonlinear Equation (36) in the plastic state, due to the nonlinearity of the constitutive relationship in Equation (6) [36,37]. In this method, the total load is divided into sufficiently small loads, which are applied step by step. Then, the nonlinear process is decomposed into a sequence of approximately linear processes.
In each approximately linear process, the relation between the stress increment tensor ∆σ and strain increment tensor ∆ε is which is linear, and D ep is only connected with the known stress and strain.
To solve this problem, first, we must determine whether the structure produces plastic deformation when applying the total external force. Under this circumstance, assume that the structure is completely in the elastic state; then the displacement, strain, stress, equivalent stress, and the maximum equivalent stress σ max of all nodes and Gauss points are obtained by using the elastic constitutive relationship in Equation (5). Denote the material yield stress as σ s .
If σ max ≤ σ s , the structure is still in the elastic state, and the above solutions are the final solutions of this problem. If σ max > σ s , the structure has produced plastic deformation, and the external force should be applied using the incremental method. Define the elastic limit load F 0 = σ s σ max F; then the displacement, stress, and strain in the elastic state can be obtained. Behind the elastic state, we use the incremental method to apply the rest of the external force, where N is the total number of load steps. In the i − th load step, ∆F i is the incremental force. Then, the solving equation is obtained, where the matrix K is related to the known stress σ i−1 in the previous load step, and ∆U i is the displacement increment in the i − th load step. Then, depending on ∆U i , we can obtain the strain increment ∆ε i and the stress increment ∆σ i in the current load step. Then the stress in the i − th load step can be obtained as, In each load step, it is still necessary to choose the elastic or plastic constitutive relationship according to the relation between the stresses σ max and σ s . Repeat the above process until the whole external force is applied. Eventually, the obtained displacement, stress, and strain are the final solutions of the elastic-plastic problem. The whole computing procedure is described in Figure 1.

Numerical Examples
Since there are no analytical solutions for nonlinear elastic-plastic problems, the ABAQUS program is used as a reference in this paper. In this section, to validate the numerical advantages of the IICVEFG method, three numerical examples under different constraints and forces are simulated and analysed. The solutions obtained from the IICVEFG method are compared with the solutions obtained from the ABAQUS program and the non-interpolating ICVEFG method [38]. To guarantee that the computed solutions are more accurate, in each load step the Equation (41) is modified to where the second term on the right side is the equivalent node load, and F i is the total equivalent node load in the i − th load step, Mathematics 2021, 9,1967 7 of 16 Based on Equation (40), the above expression can be rewritten as

Numerical Examples
Since there are no analytical solutions for nonlinear elastic-plastic problems, the ABAQUS program is used as a reference in this paper. In this section, to validate the numerical advantages of the IICVEFG method, three numerical examples under different constraints and forces are simulated and analysed. The solutions obtained from the IICVEFG method are compared with the solutions obtained from the ABAQUS program and the non-interpolating ICVEFG method [38].
For the following examples, the linear basis function and the 4 × 4 Gauss points are used in the IICVEFG and ICVEFG methods. The singular weight function and cubic spline weight function are used in the IICVEFG and ICVEFG methods, respectively [23,30]. In the ABAQUS program, the C3D8R element is selected. Additionally, the external load is applied with N = 100 steps. In this paper, two-dimensional elastic-plastic problems are regarded as plane stress problems, and the tangent modulus E in Equation (11) is given as For error analysis, define the error between the ABAQUS program and other numerical methods, and the traditional relative error at a specific node, respectively, where n is the number of the nodes, and i represents the IICVEFG and ICVEFG methods in this paper. For convergence analysis, define the following variance between the different node distributions, where u ij I represents the average displacement value of node I under two different kinds of node distributions i and j.

A cantilever Beam Constrained with a Concentrated Force
The first example is a cantilever beam; the free end is constrained with a concentrated force, as displayed in weight function are used in the IICVEFG and ICVEFG methods, respectively [23,30]. In the ABAQUS program, the C3D8R element is selected. Additionally, the external load is applied with 100 = N steps. In this paper, two-dimensional elastic-plastic problems are regarded as plane stress problems, and the tangent modulus E′ in Equation (11) is given For error analysis, define the error between the ABAQUS program and other numerical methods, and the traditional relative error at a specific node, respectively, where n is the number of the nodes, and i represents the IICVEFG and ICVEFG methods in this paper. For convergence analysis, define the following variance between the different node distributions, where ij I u represents the average displacement value of node I under two different kinds of node distributions i and j .

A cantilever Beam Constrained with a Concentrated Force
The first example is a cantilever beam; the free end is constrained with a concentrated force, as displayed in Figure 2. The geometric parameters are: and the depth is the unit length. The material parameters are:  For this example, we mainly discuss the influence of the parameter scale d max and different node distributions (as described in Figure 3) on the presented IICVEFG method. In the ABAQUS program, 64 × 8 square four-node elements are. From Figure 3, it is clear that, with the decrease in the d max value and the increase in the total number of nodes, the numerical solutions obtained from the IICVEFG method at point (8,0) are closer to the solutions of the ABAQUS program. From Figure 4 and the error analysis of Figure 5, we can also see that under d max = 2 and a 65 × 9 node distribution, the numerical deflection u 2 of the nodes on the x 1 axis of the IICVEFG method is consistent with the ABAQUS solutions.          Figure 6 shows the convergence of the IICVEFG method between different node distributions. We can conclude that the greater the number of total nodes is, the smaller the variance is. Therefore, the IICVEFG method has good numerical convergence. We should note that the IICVEFG method may spend more CPU computing time as the number of total nodes increases.   Figure 6 shows the convergence of the IICVEFG method between different node distributions. We can conclude that the greater the number of total nodes is, the smaller the variance is. Therefore, the IICVEFG method has good numerical convergence. We should note that the IICVEFG method may spend more CPU computing time as the number of total nodes increases.  Figure 6 shows the convergence of the IICVEFG method between different node distributions. We can conclude that the greater the number of total nodes is, the smaller the variance is. Therefore, the IICVEFG method has good numerical convergence. We should note that the IICVEFG method may spend more CPU computing time as the number of total nodes increases.

A Cantilever Beam Constrained with a Distributed Load
The second example is a cantilever beam where the upper boundary is constrained with a distributed load, as described in Figure 8. The geometric and material parameters are the same as in the first example. The distributed load is For this example, we mainly discuss the merits of the IICVEFG method in terms of computing accuracy and efficiency, as compared to the non-interpolating ICVEFG method. In the IICVEFG method, 4 17 × nodes are used. The grid distribution in the ABAQUS program is same as the first example.
Under different parameter scalings max d , Figure 9 shows the relative errors of the

A Cantilever Beam Constrained with a Distributed Load
The second example is a cantilever beam where the upper boundary is constrained with a distributed load, as described in Figure 8. The geometric and material parameters are the same as in the first example. The distributed load is p = −1 N/m.

A Cantilever Beam Constrained with a Distributed Load
The second example is a cantilever beam where the upper boundary is constrained with a distributed load, as described in Figure 8. The geometric and material parameters are the same as in the first example. The distributed load is For this example, we mainly discuss the merits of the IICVEFG method in terms of computing accuracy and efficiency, as compared to the non-interpolating ICVEFG method. In the IICVEFG method, 4 17 × nodes are used. The grid distribution in the ABAQUS program is same as the first example.
Under different parameter scalings max d , Figure 9 shows the relative errors of the  For this example, we mainly discuss the merits of the IICVEFG method in terms of computing accuracy and efficiency, as compared to the non-interpolating ICVEFG method. In the IICVEFG method, 17 × 4 nodes are used. The grid distribution in the ABAQUS program is same as the first example.
Under different parameter scalings d max , Figure 9 shows the relative errors of the deflection u 2 at point (8, 0) between the IICVEFG method and ABAQUS program. It is clear that when d max = 2, the relative error is close to zero. Thus, we choose these data for the following discussion. It widely known that the greater the value of d max is, the more CPU time the IICVEFG method requires. Mathematics 2021, 9, x FOR PEER REVIEW 13 of 18 With the increase in the penalty factor α the relative error is gradually reduced, and when 12 10 0 . 1 × = α the relative error is sufficiently small and close to the error of the IICVEFG method. It is worth mentioning that this selection process is time consuming. Therefore, the IICVEFG method has a greater computing accuracy and efficiency than the ICVEFG method. Taking the ICVEFG method as a comparison method, the suitable node distribution, d max , and penalty factor α should be chosen. For this example, 21 × 11 nodes and d max = 3.9 are used. The penalty factor α is chosen according to Figure 10. From the relative error analysis, there are two observations: (1) The solutions of the IICVEFG method have no relation with the penalty factor α, because in the IICVEFG method special techniques are unnecessary when applying the essential boundary conditions; (2) With the increase in the penalty factor α the relative error is gradually reduced, and when α = 1.0 × 10 12 the relative error is sufficiently small and close to the error of the IICVEFG method. It is worth mentioning that this selection process is time consuming. Therefore, the IICVEFG method has a greater computing accuracy and efficiency than the ICVEFG method. With the increase in the penalty factor α the relative error is gradually reduced, and when 12 10 0 . 1 × = α the relative error is sufficiently small and close to the error of the IICVEFG method. It is worth mentioning that this selection process is time consuming. Therefore, the IICVEFG method has a greater computing accuracy and efficiency than the ICVEFG method. After deciding the values of the relative parameters, the deflection u 2 of the nodes on the x 1 axis are acquired from the ABAQUS program, as well as the IICVEFG and ICVEFG methods. The results are shown in Figure 11. The solutions are in good agreement with each other. Moreover, for a single calculation, the CPU time spent in the IICVEFG method is 30.1 s, which is much less than the 538.0 s spent in the ICVEFG method under similar accuracy. After deciding the values of the relative parameters, the deflection 2 u of the nodes on the 1 x axis are acquired from the ABAQUS program, as well as the IICVEFG and ICVEFG methods. The results are shown in Figure 11. The solutions are in good agreement with each other. Moreover, for a single calculation, the CPU time spent in the IICVEFG method is s 1 . 30 , which is much less than the s 0 . 538 spent in the ICVEFG method under similar accuracy.   Identically to the first example, the relation between the load and the deflection u 2 at the point (8, 0) is shown in Figure 12. For this example, the structure begins to produce the plastic deformation when p = −0.2 N/m, and up to p = −0.4 N/m the structure is completely in the plastic stage.
After deciding the values of the relative parameters, the deflection 2 u of the nodes on the 1 x axis are acquired from the ABAQUS program, as well as the IICVEFG and ICVEFG methods. The results are shown in Figure 11. The solutions are in good agreement with each other. Moreover, for a single calculation, the CPU time spent in the IICVEFG method is s 1 . 30 , which is much less than the s 0 . 538 spent in the ICVEFG method under similar accuracy.

A Rectangular Plate with a Central Hole under a Distributed Traction
The final example is a rectangular plate with a circular hole in the centre. This structure bears a distributed traction at the right boundary, as described in Figure 13. The geometric parameters are: L = 10 m, h = 4 m, r = 1 m, and the depth is the unit length. The material parameters are: E = 1.0 × 10 5 Pa, ν = 0.25, and σ s = 250 Pa. The traction load is p = 1000 N/m.

A rectangular Plate with a Central Hole under a Distributed Traction
The final example is a rectangular plate with a circular hole in the centre. This structure bears a distributed traction at the right boundary, as described in Figure 13. The geometric parameters are: , and the depth is the unit length.
The material parameters are:  , respectively. Furthermore, under similar computing accuracy, the IICVEFG method requires s 84 . 408 to compute, which is less than the s 78 . 474 for the ICVEFG method. For this example, we mainly discuss the applicability of the IICVEFG model for complex structures. First, for the IICVEFG and ICVEFG methods, the numerical solutions were obtained under the same node distribution 23 × 15. In the IICVEFG and ICVEFG methods d max = 2. Moreover, in the ICVEFG method α = 1.0 × 10 9 . As a reference, 40 × 8 grids are distributed in the ABAQUS program. Figure 14 shows the displacement u 1 of the nodes on x 2 = −2 m. We can see that the solutions obtained from the IICVEFG and ICVEFG methods are in good agreement with the solutions of the ABAQUS program. The relative errors at point (5, −2) acquired from the IICVEFG and ICVEFG methods are 0.144% and 0.1611%, respectively. Furthermore, under similar computing accuracy, the IICVEFG method requires 408.84 s to compute, which is less than the 474.78 s for the ICVEFG method.

A rectangular Plate with a Central Hole under a Distributed Traction
The final example is a rectangular plate with a circular hole in the centre. This structure bears a distributed traction at the right boundary, as described in Figure 13. The geometric parameters are:    to compute, which is less than the s 78 . 474 for the ICVEFG method.  Figure 15 describes the grid distribution in the ABAQUS program, node distribution in the IICVEFG method, and the corresponding distributions of displacement u 1 . We can see that the displacement distributions obtained from the ABAQUS program and IICVEFG method are similar. see that the displacement distributions obtained from the ABAQUS program and IICVEFG method are similar.
For this example, the relation of the traction and the displacement 1 u at point ) 0 , 5 ( is also discussed in Figure 16. Around the value 100 N/m p = the plastic deformation occurs; and around 300 N/m p = the deformation of the whole structure is plastic deformation.  For this example, the relation of the traction and the displacement u 1 at point (5, 0) is also discussed in Figure 16. Around the value p = 100 N/m the plastic deformation occurs; and around p = 300 N/m the deformation of the whole structure is plastic deformation. Figure 15 describes the grid distribution in the ABAQUS program, node distribution in the IICVEFG method, and the corresponding distributions of displacement 1 u . We can see that the displacement distributions obtained from the ABAQUS program and IICVEFG method are similar.

Conclusions
In this paper, a more efficient and accurate numerical model for the two-dimensional nonlinear elastic-plastic problems is established, based on the IICVEFG method and the incremental tangent stiffness matrix method. Through numerical analyses we found that, in the IICVEFG method, the variance is smaller as the number of total nodes increases. Thus, the IICVEFG method has good convergence. The solutions of the IICVEFG method are consistent with the solutions of the ABAQUS program. In addition, the error between the IICVEFG method and ABAQUS program is less than that between the ICVEFG method and ABAQUS program. Furthermore, under similar accuracy constraints, the CPU time spent on the IICVEFG method is 30.1 s, which is much less than 538.0 s spent on the ICVEFG method. Thus, the IICVEFG method has greater numerical accuracy and efficiency than the non-interpolating ICVEFG method for different elastic-plastic problems. This validates that the IICVEFG method is sufficiently adaptable.