Research on the Elastoplastic Theory and Evolution Law of Plastic Zone Contours of Horizontal Frozen Walls under Nonuniform Loads

: The study of the changes in stress and deformation of frozen walls during excavation has always been a hot topic in underground freezing engineering, and the size of the plastic zone is an important basis for evaluating the stability of frozen walls. In response to the shortcomings in the current design of horizontal frozen walls, based on the internal excavation unloading conditions of the frozen wall in artificial ground freezing, an elastoplastic mechanical model for the interaction between a circular horizontal freezing wall and unfrozen soil mass is established under nonuniform loads to obtain the corresponding solutions for stress and displacement in the system. In this study, considering the shear stress of the plastic zone, the method for solving the traditional plastic zone contour equation is modified; consequently, the modified solution of the contour equation of the plastic zone for the frozen wall is obtained. Using this theoretical solution, the influence of the ex-ternal load 0 p and the lateral pressure coefficient  on the contour line of plastic zone and tensile stress zone are analyzed by combining the project case, the calculation results show that: the 485 . 0 =  is the critical point where the inner edge of the frozen wall just happens to have tensile stress. When 485 . 0 


Abstract:
The study of the changes in stress and deformation of frozen walls during excavation has always been a hot topic in underground freezing engineering, and the size of the plastic zone is an important basis for evaluating the stability of frozen walls. In response to the shortcomings in the current design of horizontal frozen walls, based on the internal excavation unloading conditions of the frozen wall in artificial ground freezing, an elastoplastic mechanical model for the interaction between a circular horizontal freezing wall and unfrozen soil mass is established under nonuniform loads to obtain the corresponding solutions for stress and displacement in the system. In this study, considering the shear stress of the plastic zone, the method for solving the traditional plastic zone contour equation is modified; consequently, the modified solution of the contour equation of the plastic zone for the frozen wall is obtained. Using this theoretical solution, the influence of the ex- , there will be also no tensile stress zone, with the increase of 0 p , the compression plastic zone contour evolves from the crescent shaped in the horizontal direction to the elliptical shaped, and there is no ear-shaped plastic zone in the whole evolution process. Based on our results, we show that our method can be used to provide a theoretical basis for the stability evaluation and parameter (thickness) design calculation of horizontal frozen walls under nonuniform loads.
Keywords: horizontal freezing method construction; circular frozen wall; nonuniform load; unload model; elastoplastic analysis

Introduction
Artificial ground freezing [1][2][3][4] is a special construction technology often used in the construction of underground structures in weak, water-bearing strata, such as shaft engineering, tunnel construction, and foundation pit excavation. Research on the stress and deformation of the frozen wall created during the artificial ground freezing process while 1 =  -i.e., on the elastic or elasticplastic characteristics of frozen wall under uniform load-is most common [5][6][7][8][9][10][11]. However, numerous measured data indicate that the in situ stress includes a certain degree of inhomogeneity, i.e., the ratio of horizontal stress to vertical stress (which is referred to as the lateral pressure coefficient  ) is not equal to 1 [12,13]. This inhomogeneous pressure plays an important role in the failure of underground structures.
Considering this, worldwide, research has been conducted considering nonuniform formation pressure on the surrounding rock of roadways [14][15][16][17][18][19]. Indeed, analytical solutions of models of excavation of circular tunnels in elastic homogeneous materials subject to nonuniform initial stresses (i.e., the solution known as Kirsch's solution) [20,21] are discussed in textbooks of rock mechanics and excavation analysis due to their importance as basic teaching/learning tools. Nevertheless, the research on the elastoplastic mechanical properties of the frozen wall formed by the artificial ground freezing method under nonuniform load conditions are still lacking and not explicitly developed [22,23]. In particular, under nonuniform load, the initial stress field of the frozen stratum is considerably different from that of the surrounding rock of roadways (as shown in Figures 1 and 2). In particular, as is clear from Figure 2, it is important to consider the interaction between the frozen wall and unfrozen soil in the frozen stratum while studying the mechanical characteristics of the frozen wall before the installation of a lining. However, as shown in Figure 1, before this installation, the surrounding rock of the roadway is just a single medium in the stratum; therefore, the latter case only involves an inhomogeneous stress analysis of a single medium. Therefore, the existing research results of the surrounding rock of roadways under nonuniform load are not suitable for analyzing the mechanical characteristics of horizontal frozen walls.  In the design theory of frozen walls, the vertical frozen wall theory has become relatively mature and has formed a relatively complete set of design theories [24]. However, the design of the horizontal frozen wall is not yet mature, mainly relying on the calculation formula for the vertical frozen wall of the neutral shaft in mining construction and combining the on-site experience method for design [25,26]. In recent years, scholars have focused on the design and research of horizontal frozen walls, which are still in the fields of numerical simulation, model testing, and on-site measurement. They mostly focus on the analysis of temperature and displacement fields [27][28][29][30]. Theoretical analysis only includes partial elastic design research, while elastic-plastic analysis has not yet seen relevant theoretical analysis research [6,22]. However, due to the different ground stresses in the vertical and horizontal directions, the external loads on the horizontal frozen wall are nonuniform, which fundamentally differs from the axisymmetric loads on the vertical frozen wall. Even the optimized design formula for the vertical frozen wall is theoretically not applicable to the horizontal frozen wall [23]. In order to gain a detailed understanding of the deformation and plastic zone contour characteristics of the horizontal frozen wall, and to design the horizontal frozen wall more reasonably, theoretical research on the elastic-plastic design of the horizontal frozen wall urgently needs to be addressed.
Traditionally, while studying the mechanical properties of frozen walls, the considered mechanical model was often assumed to be an infinitely long elastoplastic thickwalled cylinder under confining pressure loading conditions [6,31]. However, because the interaction between the frozen walls and surrounding unfrozen soil, as well as the radial unloading of the frozen wall, were often neglected, this led to an inaccurate analysis of the frozen walls, which consequently led to structural failures. However, in recent years, several scholars have used an unloading model to optimize their study of the mechanical properties of the frozen wall [8][9][10]32,33], which has led to suitable results. Nevertheless, analyses based on the elastoplastic theory for frozen walls using the excavation unloading model and considering the interaction between the frozen wall and unfrozen soil under nonuniform load conditions has not been conducted [34].
Therefore, in this study, considering the unloading effect of excavation, an elastoplastic mechanical model for the interaction between the horizontal frozen wall and surrounding unfrozen soil under nonuniform load conditions is established. In particular, the stress, deformation, and plastic zone contours in the elastic zone of the frozen wall are theoretically analyzed; in this case, the Mohr-Coulomb yield criterion is used as the plastic yielding criterion for the frozen wall. Based on the Mohr-Coulomb yield criterion and considering the shear stress in the plastic zone, the traditional method for calculating the plastic zone of the frozen wall is modified, and thus, a modified analytical solution of the contour equation for the plastic zone of the frozen wall is obtained. Finally, the distribution law of stress and displacement is analyzed for construction wherein the horizontal Frozen wall elastic zone Frozen wall plastic zone 1 r  r  r 0 r freezing method is used; based on this analysis, factors influencing the evolution of the elastoplastic zone are studied.

Basic Assumptions
1. The horizontal frozen wall and surrounding unfrozen soil are regarded as infinite thick-walled cylinders; in addition, the plane strain model is used to analyze the mechanical characteristics of the horizontal frozen wall. 2. The frozen wall is homogeneous and composed of ideal elastoplastic material, while the unfrozen soil is also homogeneous, but composed of linear elastic material. 3. The influence of the force of the unfrozen soil and frozen wall's own weight is ignored in the model under study. 4. All the frozen soil inside the frozen wall is excavated in one instant, and there is no support on the inner edge of the frozen wall after the excavation; i.e., the radial load on the inner edge is 0. 5. There is complete contact between the frozen wall and unfrozen soil; i.e., the radial stress and shear stress on the contact surface are equal; in addition, the radial displacement and circumferential displacement are equal as well. 6. The stress field before and after freezing remains unchanged and is the same as the initial stress field. The initial stress field in the Cartesian coordinate system is expressed as follows:  ; h is the buried depth;  is the soil weight; and  is the horizontal lateral pressure coefficient.

Mechanical Models
The mechanical model considered in our study is shown in Figures 2 and 3. In particular, Figure 2 shows the mechanical model for the initial stress field of the frozen stratum, while Figure 3 shows the mechanical model for the interaction between the frozen wall and surrounding rock after excavation and unloading. Although the contour shape of the elastoplastic zone of the frozen wall does not remain a regular circle under the action of a nonuniform stress field, for stress estimation, each contour can be approximately considered as a circle; then, the exact contour can be determined by solving for elastic stress of the frozen wall. Based on the abovementioned assumptions, the proposed mechanical model is divided into a frozen wall plastic zone, frozen wall elastic zone, and unfrozen soil elastic zone from the inside to the outside. In Figure 3, the stress boundary of the surrounding unfrozen soil at infinity is the initial stress field 0  , load at the interface between the elastic zone of the unfrozen soil and outer surface of the elastic zone of the frozen wall is 1 p , and load at the interface between the inner surface of the elastic zone of the frozen wall and outer surface of the plastic zone of the frozen wall is  p .

Basic Mechanical Equations
The plane strain mechanical model shown in Figure 3 is used to calculate the stress and deformation of the frozen wall under nonuniform loading. The basic equations of elastic mechanics used in solving this model are well known. In particular, the basic mechanical equations in polar coordinates given in the classical literature [14,35] are discussed here.
First, the equilibrium differential equation for the plane strain mechanical model can be expressed as follows: where r  and   represent the normal stress components in the radial and hoop directions, respectively, while   r represents the shear stress component. Second, the geometric equation for the plane strain mechanical model can be expressed as follows: where u and v represent the displacement components in the radial and hoop directions, respectively. In addition, r  and   represent the normal strain components in the radial and hoop directions, respectively, while   r represents the shear strain component.
Third, the physical equation for the plane strain mechanical model can be expressed as follows: where E and  are the elastic modulus and Poisson's ratio of the stratigraphic material, respectively.
The above basic mechanical equations are applicable to both unfrozen soil as well as frozen wall.

Solution of the Elastic Unloading Mechanical Model for the Frozen Wall
The elastic unloading mechanical model problem of the frozen wall can be obtained from the superposition of the equivalent initial stress field before excavation ( Figure 4) and virtual excavation disturbance stress field ( Figure 5) after the formation of the frozen wall.   surface between the frozen wall and surrounding unfrozen soil are also related trigonometrically.
Depending on the characteristics of the model boundary conditions and plane strain problem, when the semi-inverse solution method is used to solve the abovementioned problems, the plane problem using the stress method can be generalized as solving a stress function that satisfies a compatible equation. Let  denote the stress function, then, the compatible equation in polar coordinates is given as follows: If the stress function is known, the stress components can be obtained as follows: The stress functions of the unfrozen soil and frozen wall are represented using separable forms (including ) (r f i and  2 con ) (r g i ); these stress functions are similar and can be expressed as follows [36,37]: Substituting the stress functions of Equation (8) into the compatible Equation (6) yields the following result: Equation (9) is applicable to any arbitrary  . Hence, It is known that the constant term in the abovementioned stress function has no effect on the stresses and displacements of the unfrozen soil and frozen wall; therefore, we set the constant term in the stress function given by Equation (8) to 0. Hence, by solving the above differential equations, the appropriate stress function expression is obtained as follows: the undetermined parameters in the stress function satisfying the semi-inverse solution of the frozen wall are expressed; in contrast, when 2 = i , the undetermined parameters in the stress function satisfying the semi-inverse solution of the surrounding unfrozen soil are expressed. Substituting Equation (11) into the stress expression of Equation (7), the distribution function of the stress field for the frozen wall and surrounding unfrozen soil can be obtained as follows: (12) Because the stresses of the surrounding unfrozen soil and frozen wall are solved for in the same manner, we only describe the solution for the stress field of the frozen wall below. The excavation of frozen soil inside the frozen wall can be considered equivalent to a stress relief field acting on the inner edge of the frozen wall, which induces stress redistribution and deformation of the frozen wall. After the excavation of frozen soil inside the frozen wall, the relief stress inside the frozen wall can be expressed as follows: Furthermore, substituting Equation (13) into the physical as well as geometric equations and simplifying them, the elastic displacement of the frozen wall caused by the relief stress can be obtained as follows:

Determination of Undetermined Parameters
Based on the single value and continuous displacement conditions, we can say that According to the Saint-Venant principle, the original rock stress field at a sufficiently far distance will be correspondingly less affected by the excavation and unloading effect inside the frozen wall.
In particular, when  → r , the stress of the surrounding unfrozen soil after excavation and unloading of the frozen soil at the inner edge of the frozen wall can still be expressed by Equation (2). Comparing Equation (2) with Equation (12), we obtain: Furthermore, as previously mentioned, because the radial stress and shear stress at the inner edge of the frozen wall are equal to zero after excavation and unloading, the following relationship can be obtained: In a similar vein, because the radial stress, shear stress, radial displacement, and circumferential displacement at the contact surface between the frozen wall and surrounding unfrozen soil are equal, we obtain the following: Equations (8)- (11) are simultaneous equations that can be solved using MAPLE, which is a numerical software, to obtain the values of  (5) and (7). The expressions for stress and displacement of unfrozen soil can be obtained in a similar manner.

Preliminary Determination of the Plastic Zone Radius
Yang et al. [9] derived an iterative equation to solve for the radius of the plastic zone of a frozen wall based on different yield criteria and excavation unloading effects in the case when the in situ geostress field is uniform. However, in most cases, the horizontal and vertical stresses are not equal; i.e., 1   . Researchers worldwide typically study the plastic zone of surrounding rocks under nonuniform stress field conditions based on the assumption that the surrounding rock after excavation is in an elastic state, subsequently calculating the stress in this elastic state according to the elasticity theory; finally, the stress of the elastic zone is substituted into the plastic yield condition to ascertain whether the surrounding rock has yielded [38]. At this time, the contour shape of the plastic zone of the surrounding rock is no longer a regular circle; instead it has sickle, cross, or other shape [17]. Although this abovementioned method is only an approximation, it is helpful for stress estimation. Furthermore, a similar method can be used for the elastoplastic analysis of a horizontal circular frozen wall under the action of a nonuniform stress field.
In the plane strain state, the relationship between the principal stresses and three polar stress components can be expressed as follows: The Mohr-Coulomb yield criterion expressed in terms of principal stresses is as follows: where   sin 1 . Substituting Equation (19) into Equation (20), the expression of the yield condition can be obtained as follows: Then, substituting Equation (12) into Equation (21), and simplifying it, we obtain

Solution of Plastic Zone Stress
In the plastic zone of the frozen wall, the radial stress p r  is assumed to be the same as the stress distribution in its elastic zone [38]; in addition, the direction of the principal stress remains unchanged in the plastic and elastic zones. Based on these assumptions and considering Equation (21), the following relationship can be obtained: where  is the direction of the principal stress of the frozen wall at any angle  relative to the radius vector orientation of the stress ( p r  ， p   ， p r  ) at the calculation point. As indicated in Figure 6, the stress relationships can be represented in the form of a Mohr's stress circle; thus, the principal stress direction, radial stress, hoop stress, and shear stress distribution of the plastic zone of the frozen wall can be obtained as follows:

Modified Solution for the Plastic Zone Radius
The plastic zone contour equation given by Equation (22) does not consider the redistribution of stress in the plastic zone. In fact, when the frozen wall enters the plastic zone, its stress will be constantly adjusted and redistributed, resulting in continuous expansion of the plastic zone. In order to obtain a more accurate contour equation for the plastic zone than the previous one, a second approximation correction of the contour line of the plastic zone is performed considering stress redistribution of the plastic zone. In particular, for this second approximation, we use the Castner calculation method twice, wherein the initial plastic zone radius is calculated using the Castner method first, and then the stress values p r  and p r  of the calculated points are obtained using the stress equation of the plastic zone; in this case, the contour line of the plastic zone is assumed to be a circle (the circle is made per the radius of the plastic zone for the calculated point).

Engineering Examples and Analysis of the Evolution Law of the Elastoplastic Zone of the Frozen Wall
In the freezing design of the horizontal freezing method, the size of the plastic zone is an important basis for evaluating the stability of the frozen wall. In order to study the evolution law of the elastoplastic zone of the frozen wall under nonuniform load, we taking a horizontal frozen tunnel project as the background; in particular, the effects of the vertical external load 0 p and lateral pressure coefficient  on the plastic zone contour and range of the tensile stress zone of the horizontal frozen wall are analyzed.
The tunnel is constructed using artificial freezing method, and the designed average temperature of the freezing wall is −10 °C. The buried depth of the tunnel is 187.0~248 m, with a longitudinal slope of , and the initial pressure of the surrounding soil is between 0.3~0.9 MPa. The rock surrounding the tunnel is composed of N2L 2 argillaceous siltstone, sandy mudstone, sandy gravel layer, and N2L 2s loose sandstone. The elastic modulus and Poisson's ratio of frozen soil and unfrozen soil are respectively: , and the excavation radius of frozen wall is 3.3 m. The horizontal freezing curtain and freezing hole layout of the tunnel are shown in Figures 7  and 8, respectively, and the relevant geological parameters are listed in Table 1.   Artificial ground freezing is an unstable thermal conductivity problem with complex boundary conditions, and the solution to its freezing temperature field requires consideration of factors such as phase transition, moving boundaries, and internal heat sources. Based on the basic theories of soil moisture migration, heat change, thermodynamics, and poromechanics during freezing, as well as Harlan's coupled model of water and heat, COMSOL multi-physics finite element software and MATLAB were used to build the numerical calculation models of the freezing temperature field of A-A and B-B sections of the horizontal freezing curtain of the tunnel [2,39,40], as shown in Figure 9. Via finite element calculation and analysis, it was found that the thickness of the frozen wall of the tunnel was 2.76~3.25 m; for the convenience of calculating the thickness of the frozen wall, it is taken as 3 m.

Safety Analysis of Freezing Engineering
By substituting the relevant engineering and geological parameters into Equations  The stress values of the frozen wall under the elastic state can be obtained by substituting the undetermined parameters into Equation (12); the corresponding values for principal stresses, principal stress direction, and yield state for 12  k (k = 0, 1, 2, 3, 4, 5, 6) directions are listed in Table 3. Verification using Equation (21) indicates that plastic yield does not occur at the inner edge of the frozen wall of the tunnel; furthermore, the maximum circumferential stress at the inner edge of the frozen wall is 2.76 MPa, which meets the safety factor requirement of the compressive strength being twice the maximum circumferential stress. Under the actual working conditions for the tunnel considered in this study, for a frozen wall thickness of 3 m, it was calculated that, when the vertical external load increases to 1.83 MPa, the frozen wall reaches the elastic limit at the horizontal inner edge, and the vertical external load of the elastic limit is 2.03 times the design load, indicating that the tunnel freezing methodology leads to safe and reliable results. By substituting the parameters listed in Tables 1 and 2 in Equations (12) and (14), the stress and displacement distribution of the frozen wall in the No. 7 tunnel using the freezing method construction can be obtained as shown in Figures 10 and 11. (a) Radial stress (b) Circumferential stress (c) Shear stress As shown in Figure 10, because the frozen wall is subject to a nonuniform stress field, its radial stress increases at different gradients with the radius in each direction. In particular, when 1   , the radial stress in the horizontal direction increases faster than that in the vertical direction near the inner edge of the frozen wall; in addition, the magnitude of horizontal radial stress is larger than vertical radial stress. In contrast, near the outer edge, the radial stress in the vertical direction increases faster than that in the horizontal direction, and its value is also greater than that in the horizontal direction. However, for m r 31 . 5 = , the radial stresses in all directions are equal to 0.64 MPa. Furthermore, the circumferential stress in the horizontal direction is significantly greater than that in the vertical direction. Moving along the vector radius from the inner edge to the outer edge, the horizontal circumferential stress decreases, while the vertical circumferential stress remains almost unchanged; in addition, the horizontal circumferential stress at the outer edge of the frozen wall is less than the vertical circumferential stress at the outer edge. As with radial stress distribution, the circumferential stresses in all directions are equal to 1.38 MPa at m r 5.81 = . It should be noted that the shear stress is positive for angles in the range of 0-90° and 180-270°, whereas negative for angles in the range of 90-180° and 270-360°; the absolute value of shear stress is strictly symmetrical with respect to the horizontal, vertical, 45°, and 135° directions.
From Figure 11a, it can be seen that the radial displacement distribution of the elastic zone of the frozen wall shows that the vertical radial displacement of the frozen wall is greater than the horizontal radial displacement, and all of these displacements are positive values, i.e., indicating deformation to the excavation face; this can be attributed to the vertical outward load being is greater than the horizontal outward load. Furthermore, moving along a vector radius from the inner edge to outer edge, the radial displacement gradually decreases. From the circumferential displacement distribution shown in Figure 11b, we can see that the circumferential displacement is positive in the direction of 0-90° and 180-270°, negative in the direction of 90-180° and 270-360°, and the absolute value of circumferential displacement is strictly symmetrical with respect to the horizontal, vertical, 45°, and 135° directions. ), the tensile stress zone will clearly appear in the vertical direction at the inner edge of the frozen wall. In particular, when the value of 0 p is low, the frozen wall only appears as a crescent-shaped plastic zone within a certain range of the inner edge in the horizontal direction (e.g., contours of 2 MPa and 3 MPa in Figure 12a), and the radius of plastic zone in the direction of 0° and 180° is the largest. As 0 p increases, the plastic zone expands along the horizontal direction towards the outer edge of the freezing wall; in addition, it simultaneously expands along the inner edge of the freezing wall in the vertical direction. The contour of the plastic zone of the frozen wall gradually evolves from the crescent shape under low load to that of an ear-like shape (e.g., the contours of 5 MPa and 12 MPa in Figure 12a), and the plastic zone in the direction of 30-60° (the other three quadrants are similar in terms of the radius of the plastic zone) has the largest radius, which is the weakest position for plastic failure on the frozen wall. Furthermore, the horizontal inner edge of the freezing wall is in the plastic zone, while the outer edge is still in the elastic state; in contrast, for the vertical direction of the frozen wall, the inner edge lies in the tensile stress zone, whereas the outer edge lies in the elastic zone. If the nonuniformity of the two external loads is neither strong nor weak (as shown in Figure 12c with 0.5 =  ), there is no tensile stress zone in the entire frozen wall. In particular, when the value of 0 p is low, the shape of the plastic zone of the frozen wall is the same as that of

Influence of External Load on the Plastic Zone Contour
With an increase in 0 p , the contour of the plastic zone of the frozen wall gradually evolves from a crescent shape to ear-like shape (e.g., contours of 5 MPa and 12 MPa in Figure 12c), and the plastic zone in the direction of 30-60° (the other three quadrants are similar in terms of radius of the plastic zone) has the largest radius, which is the weakest position of plastic failure of the frozen wall, but its evolution process is slower than the case of 0.3 =  . Furthermore, the vertical direction of the frozen wall is entirely in the elastic zone, whereas inner edge and outer edge of the horizontal frozen wall lie in the plastic and elastic zones, respectively.
If the nonuniformity of the two external loads is weak (as shown in Figure 12e with 0.8 =  ), because the stress condition is similar to the stress state under uniform load ( 1 =  ), the contour of the plastic zone evolves from a crescent shape at the horizontal inner edge to an ellipse with an increase in the external load 0 p on the frozen wall (e.g., contour of 5 MPa in Figure 12e). Thus, in this case, no ear-shaped plastic zone is observed during the entire evolution process, and the radius of plastic zone in the direction of 0° and 180° is always the largest.

Influence of the Lateral Pressure Coefficient on the Plastic Zone Contour
The effect of the lateral pressure coefficient  on the extent and shape of the plastic zone of frozen wall for different cases is shown in Figure 13. In the case of 3MPa 0 = p , when the value of  is low, the contour of the plastic zone of the frozen wall is that of a short, thick crescent shape (e.g., contour of 0.3 =  in Figure 13a), and the range of the horizontal plastic zone is the largest. In contrast, if  is increased to 1, the crescent-shaped plastic zone gradually expands in the vertical direction along the inner edge of the frozen wall, and the horizontal plastic zone continuously reduces to form a long, thin, crescent-shaped contour, after which this contour becomes elliptical in shape (

= 
). When 1   and continues to increase, the horizontal plastic zone decreases gradually, while the vertical plastic zone increases and gradually separates in the horizontal direction, i.e., at 0°, finally forming two crescent shaped plastic zones in the vertical direction.
In the case of 5MPa 0 = p or 10 MPa, for the same values of  , the range of the tensile stress zone is the same as that of 3MPa 0 = p . In particular, when  =0.3, the plastic zone contour of the frozen wall is ear-shaped, and the plastic zone in the direction of 30-60° (the other three quadrants are similar in terms of the radius of the plastic zone) has the largest radius, which is the weakest position for plastic failure on the frozen wall.
When  increases to 1, the plastic zone of the frozen wall gradually closes in the vertical direction along the inner edge, and the horizontal plastic zone continuously reduces until an ellipse is formed (  =0.8).

Discussion on Stability Evaluation of Horizontal Frozen Wall under Uneven Load
Based on the abovementioned analysis of the influence of the lateral pressure coefficient  and vertical external load 0 p on the plastic zone contour and range of the tensile stress zone of the horizontal frozen wall, the following observations can be made: (1) When 0.485 ＜  , the tensile stress zone is clearly observable at the inner edge of the frozen wall in the vertical direction and its range is only related to  ; in particular, it increases with a decrease in  , but is independent of the magnitude of 0 p . Furthermore, 0 p only affects the magnitude of the tensile stress in the stress zone, but does not affect its range. Moreover, a reasonable vertical direction freezing reinforcement should be performed to increase the radius of the frozen wall and avoid tensile failure of the inner edge of the frozen wall. , there is no tensile stress zone as in the previous case, and the plastic zone in the 0° and 180° directions has the largest radius, indicating the weakest positions of the plastic wall of the frozen wall. Thus, the compressive plastic failure of the inner edge of the frozen wall in the horizontal direction should be avoided.
For circular horizontal freezing engineering, in practice,  should first be determined using an in situ stress test to judge whether tensile stress is present in the frozen wall, and the corresponding tensile failure of the tensile stress zone should be determined according to the tensile strength of the frozen soil obtained via the frozen soil mechanics test. By substituting the external load into the plastic contour equation of the frozen wall, the plastic zone can be determined; consequently, the maximum position of the plastic zone can be obtained. Finally, the freezing scheme can be improved based on the determined weak areas.

Conclusions
In this study, after deriving the analytical solution of the stress and displacement of the elastic-plastic zone of the frozen wall under nonuniform load, a specific engineering case is analyzed; in particular, the distribution law of stress and displacement of the frozen wall and corresponding influencing factors of the contour of its plastic zone are analyzed. Considering our observations, the following conclusions can be drawn: The point where

485
. 0 =  is the critical point at which the inner edge of the frozen wall just develops tensile stress. When 485 . 0   , the tensile stress zone is clearly observed at the inner edge of the frozen wall in the vertical direction, and its range is only related to  ; in particular, it increases with a decrease in  . Furthermore, 0 p only affects the magnitude of the tensile stress in the region, but does not affect its range. When 0.485   , no tensile stress zone is observed. When 61 . 0   , the frozen wall plastic zone contour evolves from a crescent shape to an ear shape as 0 p increases. In contrast, for 61 . 0   , the plastic zone contour evolves from a crescent shape in the horizontal direction to an elliptical shape with an increase in 0 p , and there is no ear-shaped plastic zone observed during the entire evolution process.
When 3MPa 0 = p , with an increase in  tending towards 1, the contour of the plastic zone of the frozen wall evolves from a short, thick crescent shape to long, thin crescent shape until the intersection finally becomes an elliptical plastic zone; in contrast, when 5MPa 0 = p or 10 MPa, for the same increase in  , the contour of the plastic zone of the frozen wall gradually closes along the inner edge in the vertical direction evolving from an ear shape to ellipse. Thus, this study provides a theoretical basis for the design and calculation of stress in horizontal frozen walls under nonuniform load conditions. Author Contributions: Methodology, S.P.; Validation, Y.X.; Formal analysis, S.P. and G.C.; Resources, L.P.; Writing-original draft, S.P. All authors have read and agreed to the published version of the manuscript.