Study on Freezing-Thawing Failure Process of Natural Damage Sandstone Based on Cellular Automata

Natural damage such as internal cracks, pores, and discontinuous surfaces of rocks induces the formation, expansion, and interaction of mesofractures until macroscopic failure under the action of freeze-thaw and load. The failure process is the evolution of self-organization and has localization characteristics. The CT scanning test of sandstone with different freezing and thawing cycles was carried out. According to the theory of cellular automata analysis, combined with CT image processing technology, the numerical model of the mesostructure of the sandstone constructed with natural damage, and the numerical simulation test of the Brazilian splitting failure of natural damage sandstone was completed under freezing and thawing. The results show that the combination of CT image processing technology and CASRock numerical simulation software can dynamically describe the whole process of initial damage propagation, new crack generation, and penetration of rock containing natural damage under tension in a freeze-thaw environment. Due to random and disordered natural damage inside the rock, gradually changing to an orderly state under freezing and load, the rock shows discontinuous deformation after the plastic deformation reaches a certain degree and gathers at the initial damage expansion. The fracture process of natural damaged rock under freezing and load has continuous-discontinuous deformation characteristics. The existence of initial rock damage causes the rock’s localized deformation and failure characteristics. The form and size of the initial damage determine the direction of development of the surface of the sandstone fracture and the coalescence of the crack. The initial damage area of sandstone is expanded due to the freeze-thaw cycle. After loading, the local deformation and failure characteristics caused by initial damage have delay characteristics for penetration of the main cracks in the sandstone, which determines the formation and development direction of secondary cracks in the rock failure process.


Introduction
Under long geological action and various physical and chemical actions, the internal formation of initial damage structures such as joints, cracks, pores, and discontinuous surfaces causes the heterogeneity of rocks. Heterogeneity has an important influence on the failure process, the failure mode, and the rock stress-strain relationship [1][2][3][4]. Especially in a low-temperature environment and load, the rock undergoes the elastic-plastic stage and presents discontinuous deformation after the plastic deformation reaches a certain degree. This random and disordered initial damage gradually changes to an orderly state under the action of external forces or internal forces, gathers at the maximum fracture surface, and causes rock failure, which has a signif-icant impact on the safety and stability of engineering in cold regions.
With the progress of theory and technology, many scholars have carried out research on the heterogeneous characteristics, mechanical properties, and failure mechanisms of freeze-thaw rock using theoretical research, laboratory tests, and numerical simulation and have gained some research achievements [5][6][7][8][9][10]. The constitutive model of rock is closely related to the mechanical properties of rock. To perform a more accurate quantitative analysis of the mechanical properties of rocks after freeze-thaw cycles, Lei et al. [11] studied the influence of joint shear degradation caused by freeze-thaw cycles considering joint durability on joint cohesion and friction angle through a large number of shear tests. The sheer damage evolution model of joints under freeze-thaw cycles was established, and the feasibility and ability of the model were further verified experimentally and theoretically. Lin et al. [12] established the statistical damage constitutive model of rocks under freeze-thaw cycles based on the coupled damage hypothesis of rocks and combined with the commonly used logarithmic normal distribution and strain strength theory in engineering reliability analysis. According to Feng et al. [13], to assess the degree of deterioration of rock freeze-thaw damage in the cold region, the law of energy evolution of yellow sandstone was analyzed under different freeze-thaw cycles, and the freeze-thaw damage model was established according to the relative change in the dissipation energy ratio before and after freeze-thaw. By establishing the freeze-thaw damage variable and the load damage variable, Zhang et al. obtained the equation of the evolution of total damage and the constitutive model of the freeze-thaw rock under load [14].
In laboratory test research, to study the fracture mechanism of the joined rock mass, Wu et al. [15] made rock samples with different cross cracks with wire cutting equipment and studied the failure process of rock samples under uniaxial compression using an acoustic emission system and digital image correlation technology. Zhang et al. [16] processed intact granite into prefabricated fractured granite samples with fracture angles of 0, 30, and 45, respectively. Then, the chemical corrosion test of prefabricated broken granite was carried out. The effects of chemical corrosion and natural joints on the damage characteristics and strength of rock mass were studied. Shen et al. [17] prepared single fractured sandstone samples with different angles of decomposition using similar materials and observed differences in the local damage effects of different angles of decomposition and freeze-thaw cycles in the end region of the fracture. Yang et al. [18] used a hydraulic servo rigid test system to perform uniaxial compression tests on specimens with different rock bridge widths and rock bridge inclination and analyzed the influence of rock bridge width and inclination on the initiation, propagation, and ultimate failure mode of fractured fine sandstone. Lin et al. [19] think that cumulative freezethaw action is the leading cause of rock mass degradation in a severe cold area. The cumulative effect of cyclic freezing and thawing on rock mass joints was studied by membrane pressure sensors. The rule of change in the frost heaving pressure during the freeze-thaw cycle was continuously monitored by sensors.
The theoretical research results and laboratory research results deepen our understanding of the internal strength and failure mode of freeze-thaw rock under load. However, due to the opacity of the rock, the sealing of the test process, and the rapid development of damage, it is difficult to observe the detailed crack propagation and penetration process. Therefore, many numerical methods have been developed to simulate the initiation and propagation of cracks under loading. In order to study the failure mechanism of fractured rock, Zhang et al. [20] used the bond particle model in the particle flow code (PFC) to simulate the granite biaxial compression test. Yuan et al. [21] used a uniaxial compression test and a particle flow code (PFC2D) to analyze the influence of the angle of fracture drop on the propagation of the mechanical properties of the fracture and the failure mode of red sandstone given the phenomenon that fractured rock columns are prone to collapse under underground pressure in underground engineering. Yahaghi et al. [22] performed a three-dimensional numerical simulation using a self-developed three-dimensional hybrid finite discrete element method (HFDEM). Three-dimensional numerical simulation reveals the degradation and failure mechanism of sandstone under different freeze-thaw cycles. Chen et al. [23] took Brazil tensile test as an example, the microstructure model of porous sandstone was established by the FDEM-GBM model, and the microscale failure mechanism of sandstone during the Brazil tensile test was discussed.
In recent decades, the cellular automata (CA) method, a numerical simulation program with local consideration and parallel characteristics, has been proposed to simulate the failure process of rock under load. Cellular automata describe the strong nonlinear interaction between the internal elements of discrete dynamical systems by constructing random interaction rules between the internal elements (or cells) of rock, which is an effective method to analyze complex nonlinear systems. Pan et al. [24][25][26][27][28], combined with cellular automata, elastic-plastic theory, and rock mechanics, developed the cellular automata software (CASRock) for the fracture process of engineering rock mass, which can well simulate the failure process of heterogeneous rock and the fracture process of the engineering rock mass. Pan et al. used the cellular automata software to simulate the uniaxial compression failure process of strain-softening rock materials at the mesoscale and analyze the mechanism of type I and II curves generated in the uniaxial compression failure process of rock samples. Chen et al. [29] studied the propagation and coalescence of prefabricated cracks under true triaxial compression using a true triaxial compression test and cellular automata software for the rock mass failure process (CAS-Rock). The above research uses cellular automata to analyze the damage process of rock, which provides an important reference for the mesonumerical simulation of the rock damage process.
It can be seen from the literature that the crack initiation, propagation, and mechanical properties of freeze-thaw rock have been extensively studied. However, as a heterogeneous body, there are many micropores and microcracks inside the rock, and these natural damages will have an important impact on the failure mode and mechanical properties of the rock. At present, the research on the crack evolution 2 Geofluids process of rock containing natural damage under freezethaw cycles is limited. In this study, the theory of cellular automata is introduced into the study of the freeze-thaw rock failure process. Firstly, combined with CT image processing technology, considering the natural damage structure inside the rock, a two-dimensional numerical model is constructed to characterize the initial damage to red sandstone. Then, the tensile failure process of rocks with natural damage is studied in different freeze-thaw cycles. The local deformation and failure characteristics of rock in a freezethaw environment and the mechanism of stress degradation are discussed. The influence of natural damage dominated by pore damage on the macroscopic failure mechanism and mechanical properties of rock is analyzed. It provides a new theoretical and numerical analysis method for further understanding the heterogeneous mechanical properties of freeze-thaw rock.

Materials and Methods
2.1. Rock Sample Preparation. The sandstone used for the test came from Weinan, Shaanxi Province, China. Sandstone is cemented with argillaceous sandy gravel. The surface of the sandstone is brown-red, and the texture is uniform. The internal pores are small, and the compactness is good. The particle size is within 0.05~0.25 mm. The main minerals are quartz, mica, calcite, plagioclase, and clay. According to the International Society of Rock Mechanics procedure, the rock sample is processed into a standard cylinder with 50 × 50 mm. The rock samples are shown in Figure 1.
The NN-4B nonmetallic ultrasonic test analyzer was used to measure the longitudinal wave velocity of rock samples, and rock samples with similar wave velocities were selected. There were 15 rock samples selected. All rock samples were placed in a 105°C oven for 24 h and then placed in a vacuum saturation machine for a saturation test. The basic physical parameters of red sandstone are shown in Table 1.

Test Scheme
2.2.1. Freeze-Thaw Cycle Test. A programmable TMS9018-R20 cryogenic thermostat was used for the freeze-thaw cycle test. Before the freeze-thaw cycle, the rock samples were saturated with a vacuum pressure saturation instrument for 24 hours. According to the characteristics of diurnal temperature difference in cold regions and related test procedures, combined with the research results of relevant scholars [30][31][32], each freeze-thaw cycle lasted 16 h, the cycle temperature decreased from +20°C to −20°C in 4 h, and the temperature gradient was 0.167°C/min. Freezing was at −20°C for 4 h. Then, it was increased to +20°C in 4 h, and the temperature gradient was 0.167°C/min. Then, it was melted at +20°C for 4 h in water. The conditions of the freeze-thaw cycle are shown in Figure 2. At the end of each freeze-thaw cycle, the rock sample was saturated for 12 h to replenish the water on time. All rock samples were subjected to 0, 5, 10, 20, and 30 freeze-thaw cycles.

CT Test.
In this test, the GE LightSpeed 64 row VCT spiral CT machine of Xijing Hospital was used to scan the rock samples after freeze-thaw cycles, and the scanning layer was 20 layers. Since the same scanning position is essential to trace the evolution of the sandstone microstructure after different freeze-thaw cycles, a scanner bracket is designed to place the rock sample on the bracket, and the bracket is placed in the same position on the scanner at each scan.

Test Results.
Using the VC++ program, the DICOM format image obtained by the CT scanning test is converted to the BMP format, and the CT images of different scanning layers of rock samples are obtained. Figure 3 is the CT image of the 10th layer of rock sample A1 under different freezethaw cycles.

Failure Process Analysis of Freeze-Thaw
Rock Based on Cellular Automata Theory 3.1. Mechanical Analysis Model of Cellular Automata. It is assumed that the rock that contains natural damage is a discrete dynamic system in the time dimension under the action of freeze-thaw and load. The system is made up of cells with discrete and finite states according to certain local rules. If the rock unit of unit volume is regarded as a cell, the state of the rock cell can be defined as [33]: Wave speed v p (km·s -1 ) 3.832 3 Geofluids where u and v are the displacements of cell nodes in the x -direction and y-direction, respectively; t, E, and μ are the thickness, elastic modulus, and Poisson's ratio of a cell unit, respectively; f x and f y are the node forces in the x-direction and y-direction of the cell nodes, respectively; and ε p is the equivalent plastic strain of a cell unit.
The state of the rock cells is related to the state of adjacent rock cells, so the mechanical behaviour of the rock complex system in a freeze-thaw environment can be described by the interaction between the rock cells, as shown in Figure 4.
Under the action of the freeze-thaw environment and external load, the stress of the rock mineral element in the natural damaged area does not change the movement of  4 Geofluids the area outside the neighborhood of the element; that is, the "perception" of the damage of the mineral element is local, and only the stress and movement of the surrounding element can be determined according to the movement of the surrounding element, and the information outside the neighborhood is "numb." At the same time, the rock reaches an equilibrium state through the mutual transmission of displacement and stress between its particles. Under the action of freeze-thaw and external load, the distribution of damage stress is nonuniform. Due to the complexity of the rock fracture process, the established updating rules must be able to reflect the permanent deformation of rock cell units in the stress process. In this study, the equilibrium conditions of cells are used to establish the update rule for elastoplastic cellular automata [24]. The incremental stress-strain relationship is defined as follows: where ½D e and ½D p are the elastic and plastic matrices, respectively, by linear equations: Therefore, the equilibrium equation of the rock cell node can be expressed as follows: where ΔF i ′ is the sum of the equivalent nodal forces generated by the yielding element at the central node of the cell in all cell elements, and it is assumed that Among them, fΔF i ′g j is the equivalent nodal force array of the jth element of the cell node i. For the plane, the four-node isoparametric element fΔF i ′ g j is the 8 × 1 array; ½B e j represents the strain matrix for element j.
Assuming ΔF i ′ ′ is the equivalent nodal force component corresponding to node i in fΔF i ′g j , then In this way, the nonlinear term ΔF i ′ is regarded as a part of the external load, and then, Equation (6) becomes the equilibrium equation of elasticity. Each cell in the system updates its status according to the same rules. Through force displacement → force iteration, the increments Δu i ⟶ 0 and ΔF i ′ ⟶ 0 achieve self-organizing balance.

Basic Algorithm of a Two-Dimensional Physical Cellular Automata Model for Freeze-Thaw Sandstone
(1) The study area is divided into a block matrix of N × N; each block being a cell whose position is represented by ðx, yÞ. Using PLYZ ðx i , y i , tÞ represents the breakdown threshold of cells ðx i , y i Þ when the time step is t (equivalent to the strength of the material microunit). When the generalized energy stored in cells ðx i , y i Þ exceeds that stored in PLYZ ðx i , y i , tÞ, the cells will be destroyed. XBNL ðx i , y i , tÞ is used to represent the "generalized energy" stored in cells ðx i , y i Þ at time t (equivalent to the elastic strain energy of the material). PLNL ðtÞ denotes the total number of cell ruptures in the whole system at time step t (2) The breakdown threshold PLYZ ðx i , y i , 0Þ of each cell in the initial state was set according to the Weibull probability distribution. At the same time, the initial energy XBNL ðx i , y i , 0Þ of each cell, the probability P ðiÞ of energy transfer to different directions after cell rupture, the local energy attraction domain, the generalized energy input mode, and the energy dissipation coefficient NH ðjÞ were established, where the initial rupture threshold PLYZ ðx i , y i , 0Þ is an integer of 1~4, indicating the initial strength 5 Geofluids of the material microcell (cell). The initial cell energy XBNL ðx i , y i , 0Þ is an integer of 1~4, indicating that the original stress exists in the material element. The probability P ðiÞ of energy transfer to different directions after cell rupture (i is an integer of 1~4) represents the anisotropy of the material. The local energy attraction region is equivalent to the influence region of local macroweak regions, such as cracks in the material. The generalized energy input mode is the generalized energy particle number input to the system at each time step, indicating the different stress environments of the material system. The energy It is thought that the cell reaches destructive strength and transmits its energy to adjacent cells according to the rules. At this time, the energy of cells ðx i , y i Þ changes to 0; that is, XBNL ðx i , y i , tÞ = 0. When the energy of all cells is less than their rupture threshold, it is the end of a time step.
(5) According to the rules above, the number of cell ruptures can be obtained in each time step. Since the microfracture process of heterogeneous materials such as rocks has a corresponding relationship with their acoustic emission, the number sequence of cell rupture in the model can be regarded as the acoustic emission rate sequence of materials [34] (6) The total energy dissipation of the system is represented by the sum of the energy escaped from the boundary of the system and the energy dissipated after cell destruction. The whole system satisfies the law of energy conservation in the evolution process

Establishment of the Digital Model for Mesodamage Evolution of Freeze-Thaw Sandstone
The rock contains various pores and microcracks. CT images obtained by CT scanning include CT numbers for each point in the scanning plane. CT numbers can reflect the changes in different material densities, which provides a basis for the segmentation of matrix and pores in rock CT images. In this study, the genetic algorithm is combined with the maximum interclass variance method to complete the CT image segmentation of frozen-thawed rock [35] ( Figure 5). The black area in the figure represents the initial damage area of rock. Figure 4 is converted into the DXF vector diagram format, and the digital model of the mesostructure of rock with natural damage is finally established ( Figure 6). Based on the theory of cellular automata, the numerical simulation of the Brazilian splitting experiment with natural damaged sandstone was completed by using the CASRock software. The established model ( Figure 6) is divided into two parts: pores and matrix. The red part represents the damaged area, and the blue part represents the rock matrix. The model is 50 mm in diameter and is divided into 4353 units and 4423 nodes. According to the calculation theory of the CASRock numerical software, the characteristics of rock failure are reflected in the law of inhomogeneity distribution. Table 2 shows the numerical simulation parameters. Table 3 shows the physical parameters of sandstone under different freeze-thaw cycles.
The numerical model parameters in this study are obtained by iterative correction until the stress-strain curve and tensile strength obtained by simulation are consistent with those obtained by laboratory tests. To facilitate model calibration, model parameters are divided into the following categories: (1) macroscopic parameters, including elastic modulus (E) and Poisson's ratio (v), and (2) microscopic parameters, including initial cohesion, initial internal friction angle, initial tensile strength, residual cohesion, residual internal friction angle, residual tensile strength, homogeneity, and random seed number.
For (1), it can be obtained directly from the test. The most important step in modelling using the CASRock software is to correct the microscopic parameters in (2). The initial cohesion and initial internal friction angle of sandstone were obtained through a triaxial compression test, and the microscopic parameters of sandstone were selected according to the characteristics of freeze-thaw sandstone and the research results of relevant scholars [24][25][26][27][28]33]. Subsequently, the microparameters were fine-tuned until the obtained stress-strain curves and tensile strength were consistent with the test results. In the CASRock software, the distribution of material strength parameters is controlled by homogeneity m. The greater the value, the more uniform the rock, and the higher the strength. Therefore, the homogeneity value m is also related to the test results. After several simulation tests, the homogeneity was 8, and the random seed number is 10.
Displacement loading is used in the numerical calculation; the loading rate is 1 × 10 -6 m/s. At each loading step, the deformation and stress field of the sample were solved 7 Geofluids according to the cellular automata updating rule, and the yield state of the cell unit was judged according to the Mohr-Coulomb strength criterion. When the rock cell unit reaches the yield state, the stress drops, and plastic deformation occurs. The unbalanced force generated in this process is transmitted to the adjacent cell nodes; thus, the stress of the whole system is adjusted. In the adjustment process, if the rock cell unit breaks, the stress adjustment will continue until no cell unit breaks in the calculation model.

Analysis of Mechanical Characteristics of Sandstone
Failure. The Brazilian splitting test was carried out on sandstone samples to obtain their tensile strength and stressstrain curves [36], as shown in Table 4 and Figure 7. Through calculation, the stress-strain curves of sandstone under different freeze-thaw cycles (Figure 7) are obtained and compared with the experimental stress-strain curve.
It can be seen from Figure 7 that the stress-strain curve and tensile strength obtained by numerical simulation are consistent with the results of laboratory tests. The stress of the numerical simulation curve gradually decreases after the peak point, which is different from the experimental curve. The reasons are as follows: (1) In this study, CT images of the middle part of the sandstone are intercepted for a numerical simulation of the two-dimensional Brazilian splitting. The twodimensional model can only represent some characteristics of sandstone samples. This numerical model contains a large number of pores and microcracks and has typical localized deformation characteristics during loading. These localized deformation characteristics will have an important impact on the mechanical properties of sandstone (2) Numerical calculation principle of the CASRock software. The CASRock software calculates according to the local update rule of cellular automata; that is, the "perception" of rock cells is only local, it can only determine its stress and motion according to the motion of surrounding cells, and the information outside the neighborhood is "numb." The cellular automata model believes that the stress point process from the peak yield surface to the residual yield surface is gradual. The stress drop is completed in several steps, and the strength of the element gradually deteriorates. In the process of loading, due to the uneven stress distribution, the crack occurs first in the initial damage area, and the initial damage area will affect the stress state of surrounding rock cells, which has a strain-softening phenomenon in the macroscopic view. Therefore, the curve gap between the numerical simulation and the experimental postpeak stage is large. However, this does not affect the authenticity and accuracy of the numerical simulation. This study is limited to the study of the failure mode of sandstone with natural damage under the action of the freeze-thaw cycle and does not discuss the postpeak stage of the stress-strain curve The tensile strength curve and the fitting curve of the sandstone test results under different freeze-thaw cycles were obtained by laboratory tests and numerical simulation (Figure 8). The tensile strength variation of sandstone under different freeze-thaw cycles is as follows.
In the equation, σ t is tensile strength, N is freeze-thaw cycles, and R 2 = 0:98. According to Figure 8, the tensile strength of the sandstone decreases with an increase in the freeze-thaw cycle, indicating that the freeze-thaw cycle is one of the factors causing the decline in the mechanical properties of the sandstone. However, after 20 freeze-thaw     9 Geofluids cycles, the slope of the curve gradually slowed down, indicating that the degradation of sandstone by freeze-thaw cycles was limited.
The stress-strain curve (Figure 9) of sandstone under different freeze-thaw cycles can be divided into the compaction stage, elastic deformation stage, inelastic deformation stage, postpeak stage, and residual stress stage. With the increase of freeze-thaw cycles, the initial mesodamage of sandstone is further aggravated. Under load, freeze-thaw damage induces a decrease in the macroscopic mechanical properties of the rock, and the evolution of freeze-thaw damage is an important factor influencing the gradual weakening of the mechanical properties of the rock.

Analysis of the Sandstone Failure Process with Natural
Damage. In this study, the degree of fracture of the rock (RFD) of local fracture development is used to reflect the extent of damage expansion under load under the freezethaw cycle [33]. RFD is defined as follows: where gðθÞ is the shape function on the partial plane, θ is the lode angle, p is the mean stress, q is the equivalent shear stress, ε p v is the plastic body strain, ðε p v Þj limit is the limit value, γ is the equivalent plastic shear strain, ð γÞj limit is the limit value, F represents the peak stress, F < 0 represents the prepeak stage, and F > 0 represents the postpeak stage.
If the evolution of microcracks before the peak is consistent with the Wiebols-Cook criterion, A, B, and C in Equation (11) are the correlation coefficients of the Wiebols-Cook criterion. By the analogy of cohesion and internal fric- According to the numerical calculation, the crack evolution process of sandstone containing natural damage can be obtained under load (Figure 10). The stress concentration phenomenon first occurs near the loading point, indicating that it is easy to induce cracks earlier near the loading point. With the stress loading, the crack extends along the diameter direction, and microcracks appear near the centre of the diameter. At this point, the horizontal tensile stress reaches the maximum, and the tiny elements are destroyed (yellow and red spots). Subsequently, the crack was subjected to intense tension and developed to both ends of the disk, resulting in the continuous rupture of rock cell units under stress, and the cracks were connected to form macroscopic penetrating cracks.
After microcracks appear in the centre of the sandstone, an obvious local stress concentration occurs around the initial damage, and a secondary crack is generated with an angle between the direction of the maximum principal stress and the direction of similar original pores, and the two are connected. This connection mechanism results from the local tensile stress in the rock bridge being higher than the tensile strength of the rock material due to the superposition of stress fields in two adjacent natural damage regions. Since the initial damage area on the left side of the sandstone is more serious than that on the right side, the secondary cracks are mostly distributed on the left side of the sandstone. It shows that secondary cracks will be formed in areas with more natural damage. The evolution of the crack of sandstone with natural damage under freeze-thaw and load is mainly related to the degree of damage and the distribution of the damage. The initial damage determines the direction of crack coalescence development.

Effect of Freeze-Thaw on the Failure Mode of Sandstone.
The RFD evolution of sandstone under different freeze-thaw cycles is shown in Figure 11. When the stress of sandstone reaches the peak value, the strain energy represented by the green area expands to the edge of the sample instantly, the main crack expands rapidly along the diameter direction, and the secondary crack increases, indicating that the sandstone presents an avalanche failure.
The important feature of rock from the evolution of initial damage to the failure of the final fracture is the formation of penetrating main cracks. When t = 6 s, the spacing between the upper and lower main cracks of sandstone increases with the increase of the freeze-thaw cycle (the black arrow in Figure 11 represents the main crack spacing). The main crack of the unfrozen-thawed sandstone is penetrated when t = 4 s, the main crack of the sandstone after 5 freeze-thaw cycles is  10 Geofluids penetrated when t = 6 s, and the main crack of 30 freezethawed sandstone is penetrated when t = 10 s. It can be seen that with the increase in the freeze-thaw cycle, the main crack penetration time of the sandstone is delayed. The main reason for the delay in the penetration of the main crack in the freeze-thaw cycle is that the internal stress of the sandstone under load is transferred along the path of minimum resistance. The cementation force of the mineral particles in the initial damage area of the sandstone is weak, and part of the stress is transferred along the initial damage area, and the other part of the stress is transferred along the main crack direction. With an increase in freeze-thaw cycles, the initial damage to the sandstone extends continuously outward. Sandstone stress is transmitted mainly along the initial damage area, and the stress transmitted along the main crack direction is gradually reduced, resulting in the delay in penetration of the main crack.
The degree of morphology, distribution, and development of secondary cracks in freeze-thaw sandstone under load is the reflection of microdamage and its evolution. The expansion of secondary cracks is limited by the main cracks due to energy constraints. The evolution of the crack of the split sandstone is dominated by the main crack and the secondary crack. The main crack dominates the crack  11 Geofluids evolution process of unfrozen-thawed sandstone. With an increasing freeze-thaw cycle, the initial damage area of sandstone expands, and the secondary cracks gradually increase. When the energy released by the rock sample is greater than the ability of the rock itself to resist external damage, the local position of the rock is destroyed, which is considered a new damage area of the rock. Based on the statistics of the numerical simulation results, the area-time curve of the new sandstone-damaged area (the red and yellow areas in the RFD cloud picture) can be obtained ( Figure 12). With increasing time, the area of the new damage area of the sandstone increases. Combined with the stress-strain curve of the sandstone (Figure 9), it can be seen that there are few new damaged areas of the sandstone in 0~3 s, and the sandstone is in the compaction and elastic stage. 4~10 s, the area of the new damage zone increases sharply, a large number of microcracks initiate and propagate, and the main cracks penetrate. The sandstone is in the plastic stage and the postpeak stage. After 10 s, the area of the new damage zone increases slowly, and the sandstone is in the residual stress stage. The new damage zone is mainly located near the main crack.

Effect of Freeze-Thaw on the Cumulative Energy of
Acoustic Emission in Sandstone. The cumulative energy of the acoustic emission is an important indicator that reflects the failure characteristics of the sandstone. Since the number of acoustic emissions of rock is directly related to the generation of internal damage, it is assumed that the energy release rate of acoustic emission is proportional to the release rate of elastic strain energy of a unit when a Gaussian point yields, and the release of strain energy of unit yield corresponds to the release of acoustic emission energy. The plane elastic-plastic cell automata model was used to study the law of acoustic emission of frozen-thawed sandstone during loading failure.
In this numerical calculation, the four-node isoparametric element is used, and the number of Gaussian points that are obtained is the corresponding acoustic emission number. The cumulative acoustic emission energy curve versus time is obtained by numerical calculation (Figure 13). During the loading process, the initial damage inside the sandstone continues to expand and penetrate to form a large number of new fractures, and the dislocation friction between these new fractures will produce energy of acoustic emission. It can be seen from Figure 13 that the cumulative energy of acoustic emission of sandstone is related to the number of freeze-thaw cycles. With the increase in the freeze-thaw cycle, the cumulative energy of acoustic emission of sandstone decreases, and the slope of the curve becomes slow. It indicates that with increasing freeze-thaw cycles, the initial damage to the sandstone expands, the medium inside the rock gradually loosens, and the sandstone changes from brittle failure to ductile failure. The variation of energy in the freeze-thaw sandstone failure process helps us to understand the essential characteristics of tensile failure.

Conclusions
In this article, through the CT scan test of sandstone with different freezing and thawing cycles, based on the theory of cellular automata analysis, and combined with the CAS-Rock numerical analysis software, the numerical model of the mesostructure of sandstone with natural damage was constructed, and the numerical simulation of the Brazilian fracture failure of sandstone with natural damage under freezing and thawing was completed. Through the analysis of sandstone failure mode, local fracture, stress-strain curve, and acoustic emission energy curve, the influence of the freeze-thaw environment and initial damage structure on the macroscopic failure process and mechanical properties of rock is revealed.
(1) Based on a CT scan test, digital image processing technology, and CASRock numerical simulation software, a numerical model of the mesostructure of sandstone with natural damage is constructed. The model is highly consistent with the mesostructure of real sandstone, and it can truly reflect the damage characteristics of internal pores and microcracks of natural sandstone, which is of great significance for studying the stress failure process and crack evolution law of rock (2) Based on the theory of cellular automata analysis, the basic algorithm of the two-dimensional physical cellular automata model of freeze-thaw sandstone is proposed, which can realize the numerical analysis of the failure process under the condition of considering the heterogeneity of natural damaged sandstone and can dynamically describe the whole process of initial damage propagation, new crack generation, and penetration of natural damaged rock under freeze-thaw and load (3) Due to random and disordered natural damage within the rock that gradually changes to an orderly state under freeze-thaw and load, the rock shows discontinuous deformation after the plastic deformation reaches a certain degree and converges at the initial expansion of damage. Therefore, the fracture process of natural damaged rock under freeze-thaw and load has continuous-discontinuous deformation characteristics (4) The existence of initial damage aggravates the localized deformation of the rock. The crack evolution of rock with natural damage under freeze-thaw and load is mainly related to the initial damage degree and damage distribution. The initial damage determines the direction of crack coalescence development. Secondary cracks are formed in areas with more natural damage (5) The initial damage to sandstone continues to expand due to the freeze-thaw cycle. Under the load, the localized deformation of damage increases, and the penetration of the main crack delays with the number of freeze-thaw cycles

Data Availability
The data used to support the study is available within the article.

Conflicts of Interest
The authors declare that they have no conflict of interest.