{'fraction_non_alphanumeric': 0.028114285714285713, 'fraction_numerical': 0.0096, 'mean_word_length': 4.8649463806970505, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlocal plastic softening algorithm and numerical investigation based on the Drucker–Prager criterion

Strain softening in geotechnical materials often leads to the ill condition of boundary value problems. That is, the regional governing equations or partial differential equations are no longer applicable under the original fixed solution conditions, and the over non-localization theory is used to solve the finite element equations to ensure the possibilities of the governing equations. Based on the Drucker–Prager (D–P) criterion, a nonlocal elastoplastic constitutive model and its numerical algorithm for geomaterials are established while considering the softening plasticity. Strain localization leads to high strain gradients and severe deformations in the zone. The introduction of characteristic lengths establishes the relationship between microstructure and macroscopic deformation. Accordingly, the surroundings will affect a certain point of a damaged rock material. When a Gaussian weight function is introduced, and the characterization is concerned on the damage point, the influence gradually decreases, and the characteristic length controls the influence range as the distance increases. When using the closest point projection method to calculate the stress regression of the D–P criterion, the plastic multiplier is weighted and averaged by the Gaussian weight function, and the local plastic multiplier is converted into a nonlocal plastic multiplier for numerical analysis. The results show that with the introduction of the weight function and the feature length, the localized zone width no longer depends on the grid size, that is, the grid sensitivity is overcome, and the result is more in line with the physical nature.


Introduction
For many materials, such as concrete, strain softening is a general stage for modeling the behavior in the post-peak region. As the softening goes on, a shear band will form and eventually reach failure. The material strength gradually decreases with the strain increase, leading to strain localization in a certain part of the material [1][2][3][4] . In a numerical calculation, the shear zone failure greatly differs with the number of meshes, that is, strain softening is dependent on the mesh size. To overcome this situation, regularization can be realized by introducing a material characteristic length into constitutive models, such as the micropolar [1] , strain gradient [5][6][7] , and nonlocal [8][9][10][11] models, which prevents strains from localizing into infinitely narrow bands when the mesh is refined. In a nonlocal algorithm, the plastic multiplier plays a vital role when the stress is out of the yield surface. Many models cannot handle this well [12,13] ; therefore, the plastic multiplier solution is generally considered in the algorithm presented in this paper. Meanwhile, the weight function (Gauss, triangular, and exponential) plays a vital role in the nonlocal theory. The Gauss function is smoother when the peak transitions to both ends, which is more in line with the actual situation. In this work, we select the Gauss-type function as the weight function and weight the plastic variable average. The closest point projection method has the characteristics of stability and convergence [14][15][16] . In this work, the evolution of nonlocal variables is added in the closest point projection method. All current strains are assumed to be elastic strains. Stress is calculated and brought into the yield function. If it is inside the yield surface, then 2 the variable is updated. If it falls outside the yield surface, it must be iterated. First, a set is established to store the plasticity multiplier, which is the most critical step in transforming the weighted average of the plastic variables in the yield function into nonlocal variables. The nonlocal plastic variables are then brought into the yield function. If it is in the yield surface, exit is performed; otherwise, the iteration continues until the condition is met. Second, you update the set each time as you loop over the Gauss points. The set of stored plastic multipliers is reset to zero at the end of all Gauss point cycles. The influence of the characteristic length on the value of the plastic variables and the range of the plastic zone is further analyzed on the basis of the algorithm. Finally, based on existing theoretical knowledge, the relationship between the characteristic length and the nonlocal plastic variables is assumed to be expressed by an exponential function brought into the algorithm to draw images and stress-strain curves. Its rationality is verified through an analysis.
The remainder of this paper is structured as follows: Section 2 presents an analysis of the weight function, explains the role of the weight function through data calculation and how it is applied to the yield function (i.e., how the plastic variable changes from a local quantity to a nonlocal quantity), and introduces the closest point projection method and how to incorporate nonlocal algorithms; Section 3 presents the Drucker-Prager model as an example for verifying the algorithm feasibility; and Section 4 introduces the theoretical knowledge of some scholars in the study of the characteristic length and verifies the rationality of the hypothetical exponential function through a numerical simulation. (1) is the center, and is the point around the center. Figure 1 depicts a function image of a Gauss distribution. Figure 1. Gauss function. The Gauss function is theoretically infinite, but actually close to zero beyond a certain distance from the center. The farther away from the Gauss point in the graph, the smaller the influence. Therefore, the influence within the given range, which is the external influence, is not considered zero. In Figure 2, a Gauss point was taken in each triangle. The center Gauss point is represented by p, with coordinates. The other arbitrary Gauss point is denoted by q, with coordinates. The weight is calculated as follows:

Make
, and assume that the characteristic length is . Outside this range, the Gauss graphic shows that the value is almost zero. First, we have L = 0.5 cm. The weighted average of the points within the influence range is homogenized. For example, the three weights 0.1, 0.6, and 0.9 become 0.0625, 0.3750, and 0.5625, respectively, after the homogenization treatment. That is, all weights obtained after the homogenization are added up to 1. We have 442 elements, and each element has a Gauss point. A set is assumed as . Table 1 lists the weight of each Gauss point in the model: In Table 1, represents all the columns in the first row of the set, that is, the first Gauss point and the weight of the Gauss points within the influence range, while represents all the data in the first column of the set. An analysis of the table data showed that the weight on the diagonal was the largest, as presented by the number marked in yellow in the table. This Gauss point is the center point; thus, the weight within its influence range decreases with the increase of the distance. For example, the second number in the first row is smaller than the four elements, indicating that the second Gauss point is closer to the center Gauss point. is a bridge effect in local and nonlocal, which is the weighted average of the local Gauss point values in all affected ranges. The next section will explain how to apply into the yield function.

Integral nonlocal variable
The nonlocal plasticity theory holds that it is only necessary to treat the plastic variables in the yield function as nonlocalized. Assuming that the plasticity variable is , we have (3) It is a standard non-localization theory if m = 1. It is a progressing non-localization theory if m > 1 [9,12,13] . In this work, m is 2. Therefore, is the weight set that is similar to the previous part. Assume that the yield function is: . is the Cauchy stress.
is the nonlocal plastic variable. In the iterative process, stress can only be generated by the elastic strain; therefore, The form of writing rate is (7) When judging whether or not the stress state is returned to the yield surface, the plasticity variable in the yield function is the weighted average. Hence, we must determine whether or not the yield function is less than or equal to 0. The consistency condition should be (8) In summary, we simply converted the plastic variable to a nonlocal variable .

Nonlocal numerical algorithms
As shown in the previous section, the yield function was assumed to be: . is the Cauchy stress. is the local plastic variable.
is the nonlocal plastic variable. Assume that it converged at the n step. , , and are known conditions given . Assuming that all current strains are elastic strains, we have (9) (10) The stress and the localized plasticity variable after the first iteration can be found after linearizing (1) and (4). The specific algorithm is given below. This time, must be transformed to a nonlocal plastic variable. This is the core of the algorithm. Assuming that the current Newton iteration is performed at the Gauss integration point, we have (13) is the local plastic multiplier for all Gauss integration points. We initialized zero before looping through all Gauss points and updated each Gauss point per cycle.

Closest point projection method (join weight function)
1. Assuming that the current Newton iteration is performed at the Gauss integration point , initialize , , , and given . The stress update at the Gaussian point belongs to the inner iteration and is also the core of the whole algorithm. The outer loop is given below. The Newton iteration is performed on the element as follows: 1. Input data (i.e., coordinates, element connections, and boundary conditions). 2. Update the stress, strain, plasticity variable, and consistent tangent stiffness matrix . 3. Update the imbalance force and the assembly stiffness matrix.
(29) 4. Calculate the displacement of each incremental step by modifying the stiffness matrix and the imbalance force.
(30) (31) 5. Call the algorithm in the above table . Calculate. 6. Update the imbalance force and the assembly stiffness matrix.
(32) Determine the converge according to the imbalance force. If it converges, return to 2 and continue to the next cycle; otherwise, return to 4 and use the Newton iteration until convergence.

Applications
(33) is the Cauchy stress.
is the first invariant. is the nonlocal plastic variant. H is the hardening or softening modulus. Table 2 lists the material parameters [17] .
(34)  Figure 2. Geometric dimensions and boundary conditions. In Figure 2, the geometric size is . The material strength in the black area is 2.25 MPa, while that in the other areas is 2.5 MPa. The goal was to make the black area yield first, that is, to produce a localized strain. The upper part was subjected to displacement loading. The center point of the lower part was subjected to constraints in two directions. The rest was subjected to constraints in the vertical direction only.   In Figure 3, the numbers of the elements in the grid are 442, 792, and 1232. The characteristic length in the calculation process is 6.25 cm. By comparison, the size of the plastic zone and the plastic value were approximate. Figure 4 shows the stress-displacement curves under different grid sizes. The curves almost coincided in the elastic and plastic stages. The abovementioned results prove well that the algorithm can overcome the grid sensitivity to some extent. The characteristic length in the above calculation was 6.25 cm. Accordingly, 2.25 cm, 4.25 cm, and 6.25 cm were taken for the calculation to illustrate the influence of the characteristic length.   Figure 6. Displacement curves at different characteristic lengths. In Figure 5, the following conclusions can be drawn through a comparative analysis: first, the influence range gradually increased with the characteristic length increase; and second, with the characteristic length increase came a more uniform plastic zone value and a greater influence of the surrounding area on the failure point. Figure 6 depicts the stress-displacement curves corresponding to the three characteristic lengths. The curve analysis showed that with the characteristic length decrease came a higher slope after the yield point and faster material destruction. In concrete materials, the material damage gradually increases, and the effective area gradually decreases with the progress of destruction. Therefore, the characteristic length should be a variable. The influence of the microvariables on the characteristic length will be introduced in the subsequent section. The microvariables and the characteristic length will be coupled and included in the algorithm.

Influence of microvariables on the feature scale
At present, research on the nonlocalization theory is quite extensive, while that on the characteristic length is relatively few, and there is no unified standard. However, in research on the material failure stage (e.g., for metal and concrete), the characteristic length has an important impact on its failure. In terms of theoretical research, Bazant puts the following points: first, as an important parameter of "localization restriction," the characteristic length has an important influence on the plastic or damage zone size ( Figure  5), but it should be ensured that characteristic size is always positive in the evolution process to avoid the "zero energy mode." Second, in the damage evolution process, the characteristic length will be greatly affected by the gradual expansion of microdefects and the gap reduction. This can give us the research idea of coupling the size of the characteristic length with the plastic variable (i.e., the size of the characteristic length is affected by the plastic variable). According to Bazant's theory, to avoid the characteristic length from being zero, the exponential function must be selected to affect the size of the characteristic lengths. The initial characteristic length is set as . The characteristic length of the variables is . The nonlocal plastic variable is ( Figure 6), which is the point at which the characteristic length is the same as the initial characteristic length from the elastic to plastic stage. The characteristic length gradually decreases as the plastic value increases, but is always greater than zero. The value of the nonlocal plastic variable gradually increases. Assume that the relationship of and is .
(38) Equation (38) indicates that the characteristic length is a variable in the numerical simulation. The characteristic length value and the influence range gradually decrease with the development of the microdefects. The weight set algorithm and the application of Eq. (38) in the algorithm are presented below: The biggest difference between the above table and the weight set formed by the algorithm presented above is that the characteristic length is a variable. That is, the nonlocal plastic variables are different in different calculation stages, as well as the characteristic length size. The characteristic length gradually decreases with the variable increase, which is quite similar to Bazant's point of view. However, in practical materials, the characteristic length should remain unchanged when the damage reaches a certain extent. In this algorithm, the initial characteristic length is assumed to be 6.25 cm. When the characteristic length is less than 2.25 cm, 2.25 cm remains unchanged.    Figure 7, when the characteristic length is a variable, the plastic zone range is small, and the plastic value is large. We obtain the opposite result when the characteristic length remains constant. In short, the characteristic length gradually decreases with the calculation progress. Meanwhile, the black curve in Figure  8 represents the curve when the characteristic length is constant, while the red curve represents the curve when the characteristic length is a variable. In the red curve, the slope gradually decreases in the first and second halves of the plastic phase similar to the curve with the characteristic length of 2.25 cm in Figure 6.

Conclusion
The following conclusions are drawn from this study: (1) The classical continuum theory is based on the assumption of material uniformity. Matter is believed to be continuously distributed in space. When strain localization occurs due to the high strain gradient and severe deformation in the belt, the assumption of material uniformity is no longer valid. It is precisely because of the reflection of the characteristic scale that the partial differential equation is no longer deformed during the simulation process of the strain localization, ensuring the well-posedness of the governing equation.
(2) The mesh sensitivity of plastic softening materials in the finite element calculation is overcome by combining the nonlocalized integral with the closest point projection algorithm. The result conforms to the physical essence and can handle different softening models (e.g., linear and curve types).
(3) Under the same model conditions and different characteristic scales, as the feature scale gradually increases, the shear band width gradually increases, while the plasticity value gradually decreases.
(4) By coupling the microscopic variables with the characteristic scale during the calculation process, the shear band width decreases as the characteristic scale decreases.