A novel simple practical thermal-hydraulic-mechanical (THM) coupling model with water-ice phase change

Computers and Geotechnics 118 (2020) 103357

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

A novel simple practical thermal-hydraulic-mechanical (THM) coupling model with water-ice phase change

T

Guofeng Lia, Ning Lia, , Yue Baib, Naifei Liuc, Mingming Hea, Min Yanga ⁎

a

Institute of Geotechnical Engineering, Xi'an University of Technology, Xi'an 710048, China Top International Engineering Corporation, Xi'an 710054, China c School of Civil Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China b

ARTICLE INFO

ABSTRACT

Keywords: Rock mechanics Water-ice phase change THM coupling

To reduce the threat from freeze-thaw weathering and to improve the safety of projects in permafrost areas, it is necessary to research freeze-thaw failure mechanisms and engineer safeguarding measures. Therefore, it is essential and effective to develop a reasonable, simple, practical and accurate thermal-hydraulic-mechanical (THM) coupling model with water-ice phase change. A simple model is embedded in the FLAC3D software package using FISH language. The model is developed by combining an ice-water phase change model and the THM coupling model. It is verified and extended via the uniform freeze-thaw test and the unidirectional freezethaw test, respectively. The results of the tests are in good agreement with laboratory tests, showing that the simplified algorithm is feasible. The favourable practical ability of the developed THM model makes it well suited for use in exploring multi-field change laws. It is an effective tool for analysing freeze-thaw failure mechanisms and engineering safeguard measures in permafrost regions.

1. Introduction Permafrost areas account for approximately 25% of the Earth’s land area and approximately 22.4% of China’s land [1]. Freeze-thaw weathering of rock masses poses serious threats, such as settlement and deformation on embankments [2–4], displacement of tower foundations [5], freezing injury in tunnels [6], and freeze-thaw collapse on slopes [7], and threatens the stability of geotechnical engineering in permafrost regions [8]. To reduce the occurrence of disasters, to reduce the loss of life and property, and to improve the safety and stability of projects, freeze-thaw failure mechanisms and engineering safeguard measures should be studied. They involve multi-field coupling problems under low-temperature environments, in which the water-ice phase change plays a dominant role. In addition to the relevant freeze-thaw tests, proposing a reasonable, simple and accurate numerical algorithm model is of primary importance. Water and temperature directly affect the effects of freeze-thaw weathering. The thermal field, hydraulic field and mechanical field in the rock mass interact with each other [8]. Under the freezing and thawing conditions, the interrelation relationship of the rock-water-ice system and the coupling mechanism of the thermal-hydraulic-mechanical (THM) multi-field are demonstrated in Fig. 1. When the temperature stays at a normal level, the thermal field is ⁎

influenced by convection heat transfer and moisture heat migration. The hydraulic field influences the mechanical field through pore pressure and is affected by the bulk strain. The mechanical field affects the thermal field via mechanical work and is affected by the thermal strain from the temperature change. However, a portion of the water will change into ice when the temperature drops below the freezing point. This phase change is the kernel of the relationship of the rock-water-ice system and the THM coupling mechanism. It is also the difference between the freeze-thaw process of a low-temperature rock mass and the freeze-thaw process of a general rock mass. Additionally, the phase change of the water in the pores or cracks is a unique process, and its freezing point changes with the confining pressure [8]. The phase change directly affects the temperature field. Releasing or storing latent heat allows the system to maintain a stable temperature when the temperature of the rock-water-ice system is lower or higher than the freezing point. The phase change indirectly impacts the seepage field and the stress field by the permeability coefficient and the frost heaving force. The expansion of the water volume produces additional pressure on the mechanical field during the freezing process. At the same time, the reduction in the permeability coefficient and the change in the pore pressure affect the seepage field. Once the ice pressure is greater than the tensile strength of the rock skeleton, the crack will propagate and extend. Under the condition of freezing and thawing, the compositions

Corresponding author. E-mail address: [email protected] (N. Li).

https://doi.org/10.1016/j.compgeo.2019.103357 Received 19 July 2019; Received in revised form 18 November 2019; Accepted 19 November 2019 Available online 27 November 2019 0266-352X/ © 2019 Elsevier Ltd. All rights reserved.

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

condition into a three-dimensional finite-difference program, and researched the coupled THM process with phase change by simulating the cooling test of an underground cavern for gas storage at low temperature. Huang et al. [17] deduced the governing equations for THM coupling of freezing rock under low temperature based on the energy conservation law, the mass conservation law and the principle of static equilibrium considering the water-ice phase transition. Tan et al. [18] presented a simple effective method in terms of the ‘‘three-zone’’ theory to obtain the temperature field distribution with ice-water phase change occurring in porous media. To date, studies that consider THM coupling are still limited, especially in water-ice phase change under a low-temperature environment. In addition, existing THM coupling models with phase change are too complex to be applied to actual frozen soil rock mass engineering construction for general engineers and scholars. Therefore, the authors feel it is highly necessary to conduct the current study, which aims to propose to a novel simple practical THM coupling model with water-ice phase change. This paper is structured as follows: first, the water-ice phase transition simplified algorithm model is proposed based on the energy conservation equation with water-ice phase change and the hypothesis of dividing energy stage of freezing-thawing process (Section 2). Second, the THM coupling cycle model is established by connecting the temperature, seepage and stress in series (Section 3.1). Third, a simple practical model of the THM coupling with phase change is built by combining the phase transition model and the THM coupling model (Section 3.2). Fourth, the validity of the combining model is verified via a uniform freeze-thaw test and unidirectional freeze-thaw test (Section 4), and the model is used for some applications (Section 5). Finally, the conclusions (Section 6) are indicated and presented.

Fig. 1. Interrelations of THM and the rock-water-ice system.

of the rock or soil materials include the skeleton, pores or fracture water, ice and air. To consider the main problems, the effects of the air are ignored [9]. The primary distinction between the THM coupling processes at low temperature and normal temperature is the participation of the waterice phase change. Earlier research on the THM coupling process in porous media mainly focused on the normal temperature, spurious coupling, full coupling, and a few fully coupled THM models with phase change. Wu et al. [10] applied a THM coupled model implemented in the COMSOL code to investigate the ground subsidence induced by the tunnel excavation. Chen et al. [11] proposed a THM model to reproduce the intricate coupling behaviours of compacted Gaomiaozi (GMZ) bentonite by realizing the hina-mock-up test with the help of the LAGAMINE software package. Zhao et al. [12] regarded the solid deformation, fluid seepage, and temperature field equations as independent subsystems for the THM coupling mathematical model. The scatter control equation was solved with the finite element discrete method in the geometry domain. Then, the discrete equation was obtained by using the difference method in the time domain. Bekele et al. [13] presented an isogeometric analysis (IGA)-based numerical model for simulation of the THM coupled processes in ground freezing. Tan et al. [9] established the governing equations for the coupling model depending on the basic law of water flow and heat transfer in porous media under the freeze-thaw condition and based on the theories of continuum mechanics, thermodynamics and segregation potential. Wu et al. [14] developed a numerical model of full coupling of the THM processes incorporating the equilibrium, motion, constitutive and compatibility equations, which provided good descriptions of the THM behaviour of subway tunnels in a tunnel geomaterial system (TGS) and their evolution. Cao et al. [15] designed a transient, three-dimensional numerical model in a novel and consistent formulation describing the fully coupled THM processes to simulate and analyse the enhanced geothermal system (EGS) heat extraction process. Teng et al. [16] developed a fully coupled THM model validated to have great applicability for coal deformation, gas flow and heat transport based on a series of thermal coal-gas interactions, such as thermal expansion, nonisothermal gas sorption and thermal fracturing. Kang et al. [8] proposed the frozen ratio function considering time-dependent behaviour, implemented the THM coupling equations under the freeze-thaw

2. Ice-water phase change model The successive freezing processes of water droplets on an ice surface [19], the detailed dynamic motions of a water droplet impacting on an ice surface [20], and the formation time of initial ice nuclei and freezing propagation velocity on super-hydrophobic surfaces [21] have been thoroughly researched at a macroscopic level. However, the processes of water flow and heat transport including ice-water phase change in a rock-water-ice system or in porous media are fairly complicated and significant in various fields of science and engineering. In nature, heat energy always transfers from high temperature to low temperature, and it finally reaches a relatively stable state or trends to balance. Nevertheless, the temperature reduction or increase in the rock-water-ice system, in addition to being affected by the ambient temperature, is also subjected to the latent heat of the water-ice phase change in the effective freezing temperature range. 2.1. Stages of the freeze-thaw process During the water-ice phase transition process, the system temperature remains near the freezing point, and the temperature of the rock matrix is not higher or lower than the temperature of the rock-water-ice system. In other words, in the effective freezing temperature range, releasing or storing the latent heat triggered by the water-ice phase change is offset by decreasing or increasing the system temperature caused by cooling or warming of the external environment temperature. Therefore, when the outside temperature is decreasing and the system temperature is about to become less than the freezing point, the latent heat releases, thus keeping the temperature of the system near the freezing point. The system temperature decreases again until the latent heat is completely released. In contrast, when the outside temperature is increasing and the system temperature is about to become higher than the freezing point, the latent heat is absorbed, thus keeping the temperature of the system near the freezing point. The system temperature rises again until the latent heat is completely stored. 2

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Therefore, the energy changes of each stage are assumed as follows:

expressed as follows [8,10,16,17,23]:

(I) Freezing stage A. For the cooling stage before the water-to-ice phase transition, the corresponding heat differential caused by the system temperature decrease, ΔQ1, is equal to the heat that the external environment absorbs, ΔQ2. B. For the super-cooling stage, the latent heat is stored; meanwhile, the corresponding heat differential caused by the system temperature decrease, ΔQ1, is equal to the heat that the external environment absorbs, ΔQ2, plus the latent heat, Lf. C. For the stage of the water-to-ice phase transition, the latent heat releases; in the meantime, the system temperature remains constant. That is, the system heat differential, ΔQ1, is zero. In other words, the latent heat released, Lf, offsets the heat absorbed by the external environment, ΔQ2, such that the system temperature is not reduced. D. For the cooling stage after the water-to-ice phase transition, the corresponding heat differential caused by the system temperature decrease, ΔQ1, is equal to the heat absorbed by the external environment, ΔQ2. (2). Thawing stage A. For the warming stage before the ice-to-water phase transition, the corresponding heat differential created by the system temperature increase, ΔQ1, is equal to the heat released from the external environment, ΔQ2. B. For the stage before the ice-to-water phase transition, the latent heat releases; meanwhile, the corresponding heat differential caused by the system temperature increasing, ΔQ1, is equal to the heat released from the external environment, ΔQ2, plus the latent heat, Lf. C. For the stage of the ice-to-water phase transition, the latent heat is stored; in the meantime, the system temperature remains constant. That is, the system heat differential, ΔQ1, is zero. In other words, the latent heat stored, Lf, offsets the heat released from the external environment, ΔQ2, such that the system temperature is not elevated. D. For the warming stage after the ice-to-water phase transition, the corresponding heat differential caused by the system temperature increase, ΔQ1, is equal to the heat released by the external environment, ΔQ2.

CT

T + t

·qT +

w Cw q w ·

T + L i nSr

ui + qvT = 0 t

(1) T

where T is the temperature of the rock [°C], t is time [s], C is the effective bulk specific heat of the rock mass [kJ/(m3·°C)], ρw is the density of water [kg/m3], Cw is the specific heat of water [J/kg. °C], qw is the fluid specific discharge, qT is the thermal flux,qvT is the volumetric heat source intensity, L is the latent heat of the water-ice phase transition [kJ/kg], ρi is the density of ice [kg/m3], n is the porosity, Sr is the saturation, and ui is the frozen water ratio. Among them, the first item is the system temperature difference caused by the heat change; the second item is the divergence of the temperature gradient or the heat conduction, which is in line with the Fourier law caused by the temperature gradient; the third item is the change in water and heat generated by seepage, the fourth item is the latent heat; and the fifth item is the heat source. 2.3. Developing a water-ice phase transition simplified algorithm model The latent heat term in the energy conservation equation with the water-ice phase transition is not built-in in the THM coupling energy balance equation of FLAC3D. Although scholars have used some methods to equivalently treat this oversight, the calculation theory is relatively complex for beginners or engineers. Therefore, in this paper, the previous simplified algorithm [24] is improved and upgraded. The latent heat term is replaced by the energy change due to the system temperature difference. The flow chart of the phase transition that considered the cumulative energy continuity of the unit bi-directional phase transition is shown in Fig. 2, and the default order is from freezing to thawing. The ideas are as follows: When the environment temperature, T(t), decreases gradually, the temperature of each unit is traversed. Until the element temperature, T(2), is lower than the freezing point (T0), the element enters the stage of the water-to-ice phase change. If the cumulative energy, ΣQi = ΣQiT 1 + ΔT• C , from the temperature difference, ΔT = |T2-T0|, is less than the cell latent heat, Q, it continues to reset the cell temperature to the freezing point, T2 = T0, and accumulate energy. When the cumulative energy is greater than the cell latent heat, ΣQi > Q, the water-to-ice phase transition of the element is complete. When the environment temperature, T(t), increasingly rises, the temperature of each unit is traversed. When the element temperature, T(2), is higher than the freezing point (T0), the element enters the stage of the ice-to-water phase change. If the cumulative energy, ΣQi = ΣQi-1ΔT• CT, from the temperature difference, ΔT=|T2-T0|, is greater than zero, it continues to reset the cell temperature to the freezing point, T2 = T0, and it continues to accumulate energy. When the cumulative energy is less than zero, ΣQi < 0, and the ice-to-water phase transition of the element is complete.

The overheating phenomenon is not found in the thawing stage of the samples; at the same time, the super-cooling phenomenon of a partial sample is also barely observed in the freezing stage [22]. In cases of high salt content, low water content, and little free water in the soil, the freezing temperature may not be reached. If the moisture content is very low, in a relatively short period of time, the free water in the soil will be completely frozen quickly such that the super-cooling stage and the constant-temperature stage may disappear. Furthermore, the times of latent heat release and storage are not certain, and the mechanism and the process stage of freezing and thawing are not clear, such that the stage before the phase transition, B, could be merged with the stage of phase transition, C, for simplicity in a simulation calculation.

3. THM coupling model with water-ice phase transition 3.1. THM coupling model The mechanisms of convective heat transfer in porous media considered in the FLAC3D implementation are forced convection, in which the heat is carried by the fluid motion, and free convection, which accounts for fluid motion caused by density differences due to temperature variations. The working assumptions are those of local thermal equilibrium and small density variations [23]. Considering that the multi-field coupling process involves unsaturated water migration, heat transport, and elasto-plastic deformation, the following assumptions are made to facilitate the research work: the rock mass is a multiphase system containing rock grain, water and ice. The rock material is homogeneous, isotropic and continuous. The fluid flow satisfies Darcy’s law, and the heat flux satisfies Fourier's law [8]. The coupling between

2.2. Energy conservation equation with the water-ice phase transition Phase transitions occur everywhere in nature and industry. Phase transitions involving water are the most important and common. The FLAC3D software package is proficient at simulations involving THM coupling geotechnical engineering problems. However, the built-in coupled THM calculation method does not include the phase change process [8]. THM coupling equations under the freeze-thaw condition and the freezing pressure could be applied using FISH programs in a calculating procedure. The energy conservation equation of unsaturated rock mass considering both convection and conduction can be 3

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Fig. 2. Ice water phase change model.

the mechanical and hydraulic calculations due to deformation can be implemented using the Biot coefficient. The temperature variation can lead to heat-induced expansion and contraction stress in the geo-materials. Hence, the mechanical and thermal processes can be coupled with each other through the thermal expansion coefficient [10]. The energy balance equation has been given in Section 2.2, and other governing equations are given below.

where P is the pore pressure, T is the temperature, M is the Biot modulus, is the volumetric strain, is the Biot coefficient, β is the volumetric thermal expansion of the porous matrix, and qw is the fluid specific discharge. (3). Fluid transport equation Assuming that the water fluid satisfies Darcy’s law, qw can be expressed as follows [10,13,16,23]:

(1). Momentum balance equation

qw =

The freezing rock is assumed to be a continuous isotropic medium; in the case of static equilibrium of the medium, the equilibrium equation or the general linear momentum balance equation can be expressed as follows [8,13,17]:

w

where σ is the total stress tensor, g is the vector of gravitational acceleration, and ρ is the average apparent density of the porous multiphase material [8,13,23].

= nSr [(1

ui )

w

+ ui i ] + (1

n)

s

(3)

+

T t

f (T

(6)

T0)]

Energy balance for convective diffusive heat transport is expressed by Eq. (1), where CT is the effective specific heat, defined as [8,23]

CT = nSr [(1

The mass conservation equation of water can be expressed as follows [23]:

t

0 [1

(4). Energy balance and heat transport equation

(2). Fluid mass balance equation

·qw

=

where f is the volumetric thermal expansion of the fluid and T0 is the reference temperature. In these equations, the fluid density variations due to temperature changes are significant only in their generation of buoyancy forces.

where n is the porosity, Sr is the saturation, and s , i , and w are the densities of the rock grain, ice and water, respectively. ui is the frozen ratio that represents the percentage of water changed into ice.

P =M t

(5)

w g)

where k is the fluid mobility coefficient (intrinsic permeability, κ, over dynamic viscosity, μw) and ρw is the fluid water density. In this equation, the fluid density is related to temperature changes by the linear equation [23]

(2)

· + g=0

k (P

ui )

w Cw

+ ui i Ci] + (1

n) s Cs

(7)

where n is the porosity; ui is the frozen ratio; Sr is the saturation; and s , i , and w are the densities of the rock matrix, ice and water, respectively. Cs ,Ci , and C w are the specific heats of the rock matrix, ice and water, respectively.

(4) 4

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

The conductive heat flux, assuming isotropic thermal conductivity, can be expressed using Fourier’s law as [8,13,23]

qT =

off successively. It is only necessary to ensure simultaneity for the calculation time of each field. Of course, parallel connecting or mixed connecting of each field embedded algorithm to simulate THM coupling is also feasible. However, these procedures will be more complicated, so they will take more time [24].

(8)

kT T

where kT is the effective thermal conductivity, which is isotropic in the advection formulation. The effective thermal conductivity is defined in terms of the solid and fluid thermal conductivities, k sT and k Tw , by the equation [8,23]

kT

= (1

n)k sT

+

3.2. THM coupling model with phase transition

(9)

nSr k Tw

The THM coupling model with phase transition is obtained by nesting the phase transition model and the THM coupling model. The flow chart is shown in Fig. 4. To keep the synchronicity of the THM time, the sub-step is set, and the step loop is also adopted. In a calculation cycle, detecting the environment temperature and element temperature, cumulating the heat quantity by applying the phase change model, updating the element temperature and parameters, and carrying out the THM coupling model are performed successively. The THM coupling model with phase transition is built-in the FLAC3D software using the FISH language. The basic grid model, initial conditions, changing boundary conditions, basic parameters changing relationships with time and temperature, etc., should be added to the calculation and analysis process.

(5). THM coupling Coupling to the mechanical constitutive equations is done via effective normal strain rates and pore pressure changes, which affect the effective stresses. The advective logic may be used in combination with any of the mechanical constitutive models available for FLAC3D. In the example of an elastic material in small strain, the mechanical constitutive relations are [23] ij

t

+

P t

ij

= 2G

ij

t

t

T t

ij

+ K

2 G · 3

kk

t

3

t

T t

ij

(10) where ij and ij are total stresses and strains; is the Biot coefficient; K and G are the bulk and shear moduli, respectively; t is linear thermal expansion coefficient of the solid matrix; and ij is the Kronecker delta. At present, there are many algorithm models of multi-field coupling, including the full coupling of governing equation relations and the simplified coupling of algorithm programs. Among them, the study of the full coupling of governing equation is more complicated, as shown above. Furthermore, the rock mass is usually rich in micro-cracks or flaws. If a more-complex inner structure is considered, the governing equations of THM coupling should be much more complicated [8]. To facilitate the work and engineering applications, some simplifications and assumptions need to be made to study the THM coupling of rock mass under freeze/thaw conditions. Regarding the latter, the simple THM coupling model established by connecting the basis embedded algorithm of the temperature field, seepage field and stress field in a series is more practical. The diagram is shown in Fig. 3. In the calculation process, the field switches are turned on or turned

4. Model validation This section will prove the validity and feasibility of the THM coupling model with the water-ice phase transition by using the uniform freeze-thaw test, case 1. It will also explore the deformation and temperature change law by extending the unidirectional freeze-thaw test, case 2. The temperature changes of all the boundaries on a cylindrical sample are uniform and simultaneous for the uniform freezethaw test, as shown in Fig. 5. The sample is kept frozen for 6 h under −15 °C or −30 °C conditions and then allowed to thaw for 6 h under 16 °C conditions [25]. However, the temperature changes of all boundaries on the cylindrical sample are non-uniform and simultaneous for the unidirectional freezing and thawing test, as shown in Fig. 6. The sample is kept frozen for 6 h under −15 °C conditions on the top boundary and then allowed to thaw for 6 h under 15 °C conditions on the top boundary, while the bottom boundary is be kept at 2 °C for the entire 12 h [26]. 4.1. Experimental setup The simulation test uses the FLAC3D software package. The specimen is cylindrical [25], with a diameter and height equal to 50 mm and 100 mm, respectively. The profiles, positions of all key monitoring points, boundary conditions, and initial conditions for the two cases are shown in Fig. 7. The M-C model, isotropic conduction convection model and isotropic seepage model are adopted when calculating the coupling of the stress field, temperature field and seepage field, respectively. The initial stress field is generated by the dead load of the rock sample and the initial temperature is room temperature, i.e., 16. At the same time, saturation of the rock is assumed. A normal stress boundary is applied only at the bottom of the specimen. It is assumed that there is no fluid exchanged with the outside for the boundary nodes, so an impermeable boundary is adopted. Because an air bath is used for the specimen surface of the constant-temperature box, the convection temperature boundary is employed to all the sample boundaries for case 1. However, because a liquid bath is used on the top and bottom of the temperature control board, convection temperature boundaries are only employed at the top and bottom sample boundaries, and adiabatic temperature boundaries are used on the side boundaries for case 2.

Fig. 3. THM coupling model. 5

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Fig. 4. THM coupling model with phase transition.

4.2. Parameter discussion As the changes in the Poisson’s ratio and internal friction angle in the freeze-thaw cycles are so small, they may be considered constant and equal to the initial values. It is assumed that only the elastic modulus and cohesion will change with the temperature and time of the freeze-thaw cycles. The initial parameters of the thermal field and hydrological field as well as stress field of the rock sample are reported in Table 1. The relationships of the compressive strength, Rb, the elastic modulus, E, and the internal cohesion, C, with the surrounding temperature, T, and the times of the freeze-thaw cycles, n, are as shown below [24–28]:

TR = 0.0004t 2 TE = 0.0005t 2 TC = 0.0001t 2 KR = KE = Kc =

Fig. 5. Model of uniform freeze-thaw test.

0.0127t + 0.6026 0.0147t + 0.5329 0.0078t + 0.8058

0.00001n3 + 0.0010n2 0.00002n3 + 0.0018n2 0.00002n3 + 0.0018n2

0.0310n + 1.0 0.0560n + 1.0 0.0583n + 1.0

(11)

(12)

where n is the number if freeze-thaw cycles, t is the surrounding temperature, and TR, TE, and TC are the corresponding reduction factors for the compressive strength, elastic modulus and internal cohesion, respectively. Kr, KE, and KC are the corresponding coefficients of freezing and thawing for the compressive strength, elastic modulus and internal cohesion, respectively. The specific heat (CT) and thermal conductivity coefficient (K) of rock can be considered to be unchanged in the relatively small temperature range. When the temperature is higher or lower than the freezing point, the water or ice parameters, respectively, are used. When the water-to-ice phase change or the ice-to-water phase change occurs, the water or ice parameters, respectively, are used as well. Based on the results [24–28], the thermal expansion coefficient (β) of the rock matrix is 0.3 × 10–5°C−1, but in the water-ice phase change, it is −1.5 × 10–5°C−1; if the ice-water phase change is almost restorable, it is −1.2 × 10–5°C−1. The convective heat transfer coefficient (h) is 5–25w/m2 °C and 20–100w/m2 °C under natural convection and forced convection, respectively. Thus, it is 5w/m2 °C and 20w/m2 °C for the stationary air and unidirectional freezing and thawing ends, respectively.

Fig. 6. Model of unidirectional freeze-thaw test.

4.3. Test results 4.3.1. Case 1: Uniform freezing and thawing test The sample undergoes freezing and thawing under (−15 °C/ 6

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Fig. 7. Profiles, key monitoring points, boundary conditions, initial conditions.

−30 °C ~ 16 °C) conditions. The temperature duration curve and freeze-thaw phases of the centre A3# node are shown in Fig. 8(a), and the axial radial strain duration curve is given in Fig. 8(b). The state of freezing and thawing with time is provided in Fig. 9. The net deformation with time is presented in Fig. 10. The temperature, axial strain and radial strain of the centre node are reported in Table 2. Fig. 8 and Table 2 indicate that during the time that the environment temperature stays at −15 °C for 6 h, the temperature of the sample entering into the freezing stage decreases gradually from room temperature, i.e., 16 °C to −15 °C, and experiences three phases in turn. In the cooling stage before the water-to-ice phase transition, C-B-WTI,

the temperature decreases approximately linearly for 57 min, and the cooling rate is approximately 0.3 °C/min. The strain decreases linearly and attains the maximum value at 57 min. The radial and axial strain rates are −0.4 × 10−6/min. In the stage of the water-to-ice phase transition, WTI, the temperature remains in the vicinity of the freezing point for approximately 63 min. Compared with the temperature, the strain delays approximately 18 min. The radial and axial strain rates are 2.7 × 10−6/min and 2.4 × 10−6/min, respectively. In the cooling stage after the water-to-ice phase transition, C-A-WTI, the temperature decreases approximately exponentially for 240 min. The radial and axial strains decrease by 13.8 × 10−6 and 13.5 × 10−6, respectively.

Table 1 Thermal hydraulic and mechanical properties of rock water and ice. Parameter Thermal properties Thermal conductivity coefficient of rock, Ks Thermal conductivity coefficient of water, Kw Thermal conductivity coefficient of ice, Ki Specific heat capacity of rock, Cs Specific heat capacity of water, Cw Specific heat capacity of ice, Ci Thermal expansion coefficient of rock, βs Thermal expansion coefficient of water, βw Thermal expansion coefficient of ice, βi Latent heat of fusion, Lf Hydraulic properties Hydraulic conductivity of rock, ks Saturability of rock, Sr Initial porosity of rock, n Mechanical properties Density of rock, ρs Density of water, ρw Density of ice, ρi Initial bulk modulus, K Initial shear modulus, G Initial internal cohesion, C Initial internal friction angle, Φ Initial compressive strength, Rb

Value 1(Above the ice point)

Value 2(Below the ice point)

Unit

2.7 0.52 / 816 4190 / 3.0 × 10−6 2.1 × 10−4 / 334

2.7 / 2.21 816 / 1880 −1.5 × 10−5 / −1.5 × 10−3 334

w/m·°C w/m·°C w/m·°C J/kg·°C J/kg·°C J/kg·°C °C−1 °C−1 °C−1 kJ/kg

1.0 × 10−7 1.0 0.1

1.0 × 10−10 1.0 0.1

cm/s – –

2100 1000 / 6.63 3.06 6.53 41.0 27.27

2100 / 910

kg/m3 kg/m3 kg/m3 GPa GPa MPa ° MPa

7

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Fig. 8. Temperature, axial radial strain, and stage division of the centre node under one freeze-thaw cycle. Note, ① is the cooling stage before the water to ice phase transition, C-B-WTI; ② is the stage of the water to ice phase transition, WTI; ③ is the cooling stage after the water to ice phase transition, C-A-WTI; ④ is the warming stage before the ice to water phase transition, W-B-ITW; ⑤ is the stage of the ice to water phase transition, ITW; and ⑥ is the warming stage after the ice to water phase transition, W-A-ITW.

Fig. 9. State of freezing and thawing with time.

Fig. 10. Net deformation with time (m). Note, the red net represents the original grid and the black net represents the freeze-thaw grid with a relative amplification coefficient of 3000. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

During the period when the environment temperature reaches room temperature and is maintained at 16 °C for 6 h, the temperature of the sample entering into thawing stages rises gradually from −15 °C to room temperature, i.e., 16 °C, and similarly experiences three phases in turn. In the warming stage before the ice-to-water phase transition, WB-ITW, the temperature approximately linearly increases for 45 min and the heating rate is approximately 0.3 °C/min. The strain decreases linearly and attains the maximum value at 45 min. The radial and axial strain rates are 0.4 × 10−6/min. In the stage of the ice-to-water phase transition, ITW, the temperature remains near the freezing point for approximately 48 min. The strain decreases rapidly. The radial and axial strain rates are −2.5 × 10−6/min and −2.2 × 10−6/min, respectively. The strain reaches the maximum value approximately 24 min after with the temperature. In the warming stage after the iceto-water phase transition, W-A-ITW, the temperature increases

approximately exponentially for 267 min, but the strain increases slowly to achieve a stable state. The radial and axial strains increase by 15.0 × 10−6. As can be observed from Figs. 9 and 10, when the environmental temperature remains at −15 °C, the sample temperature gradually decreases, and the sample is simultaneously narrowing in a nearly cylindrical fashion. Beginning with WTI, the sample freezes inwards, starting from the four corners, and then it freezes through the two ends and the side boundaries. The sample expands in the shape of slender waist dumbbells, towards the waist from both ends, until freezing completely. After it is placed into the room-temperature environment of 16 °C, the sample temperature rises gradually. Meanwhile, the sample thaws inwards, starting from the four corners, then through the two ends and the side boundaries also. The sample looks similar to a wooden beer barrel, which shrinks towards the centre at both ends. 8

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Table 2 Analysis of the temperature and axial radial strain of the centre node. Surrounding temperature °C

−15

16

−30

Stage

①

②

③

④

⑤

⑥

①

②

③

④

⑤

⑥

Centre temperature Interval time (min) Variation (°C) Rate of change (°C/min)

57.0 −16.4 −0.3

63.0 −0.4

240.0 −13.9

45.0 15.1 0.3

48.0 0.2

26.70 14.9

36.0 −16.9 −0.5

30.0 −0.3

294.0 −28.5

69.0 29.9 0.4

60.0 0.7

231.0 14.4

Axial strain Interval time (min) Variation (×10−6) Rate of change (×10−6/min)

57.0 −20.5 −0.4

81.0 193.9 2.4

222.0 −13.5

45.0 18.1 0.4

72.0 −158.9 −2.2

243.0 15.0

33.0 −28.0 −0.9

45.0 255.0 5.7

282.0 −33.8

69.0 44.7 0.7

75.0 −203.0 −2.7

216.0 7.2

Radial strain Interval time (min) Variation (×10−6) Rate of change (×10−6/min)

57.0 −20.5 −0.4

81.0 215.0 2.7

222.0 −13.8

45.0 18.4 0.4

72.0 −178.6 −2.5

243.0 15.0

36.0 −27.4 −0.8

45.0 283.2 6.3

279.0 −36.2

72.0 46.5 0.7

75.0 −226.0 −3.0

213.0 17.7

Then, it expands into a cylindrical shape until the thermal expansion completes after melting. The four corners are in a bi-directional heat transfer state everywhere such that the heat transfer and the freezethaw rate are faster. For C-B-WTI and W-B-ITW, the changes in temperature and strain are basically synchronous. At the same time, the rates of temperature change are equal, and the rates of strain change are also equal. The reason may be that it only relates to the thermal expansion and the cold contraction from the temperature conduction of the rock sample. For WTI, the water turns into ice such that the specific heat decreases and the thermal conductivity is larger, which delays the temperature conduction and leads to the reduction of the rate of the temperature change and strain change. For WTI and ITW, because of the latent heat participation and the changes in specific heat and thermal conductivity under the condition of keeping the system temperature constant, the rates of temperature change and strain change increase. There is one positive strain peak value in the freezing stage and the thawing stage. The strain peak value of the thawing stage is slightly larger, which is very similar to the phenomenon that temperature warming causes the frozen rock mass to melt, the slope to slide and the road to collapse in practical engineering. When the environmental temperature is −30 °C, the temperature sustaining times are 30 min and 60 min for WTI and ITW, respectively. However, the strain increasing time is 15 min greater for each of them. The radial strain rate is greater than the axial strain rate, and the strain rate for WTI is greater than the strain rate for ITW. Compared with −15 °C, the freezing environmental temperature decreases such that the C-B-WTI and WTI times are shortened by 36 min and 30 min, respectively. Moreover, the temperature change rates of C-B-WTI and WB-ITW are accelerated by −0.5 × 10−6/min and 0.4 × 10−6/min, respectively. The strain and strain change rate increase overall, but the strain time of WTI decreases. That is, the lower the freezing temperature is, the greater the frost heave strain rate is, and the shorter the time of reaching the strain peak value is. The WTI time is shorter than that of the laboratory for the two freezing temperatures, which may be related to the water content or porosity of the rock sample.

16

Fig. 11. Radial temperature changes with time.

Fig. 12. Radial strain changes with time.

16.

Figs. 11–14 and Table 3 indicate that during the whole freeze-thaw process, under the temperature gradient from the cold top end to the warm bottom end, the rock sample also experiences 6 stages (cold contraction, frost heaving, cold contraction, heat expansion, melt contraction and heat expansion) in turn. The 6 stages of the strain are relatively obvious, but the 6 stages of the temperature are not clear. The phase transition temperature platforms display only the 1# and 2# monitoring points, except for ITW. For C-B-WTI, the continuous radial strain times of each of the monitoring points are all less than the axial direction. The sustaining time of the radial strain and axial strain gradually decreases from the bottom to the top of the specimen and are 9 min and 12 min, respectively, at the top end. For W-B-ITW, the

4.3.2. Case 2: Unidirectional freezing and thawing test The top cooling end and the bottom warming end stay at −15 °C and 2 °C for 6 h, respectively, then the top end begins to warm to 15 °C, but the bottom end still remains at 2 °C for six hours. The axial temperature, axial strain and radial duration curves for each of the monitoring points are shown in Figs. 11–13. The temperature and strain comparison curves of the key points A3 of the uniform freezing and thawing test and the strategic point 3# of the unidirectional freezing and thawing test are shown in Fig. 14. The state of freezing and thawing with time and the net deformation with time are shown in Figs. 15 and 9

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

and the bottom of the specimen. However, the axial strain rates of WTI are all greater than that of ITW only at the waist of the specimen. For CA-WTI, the radial strain at the waist is greater than that at the ends. It gradually increases from the bottom to the top for W-A-ITW. However, for C-A-WTI and W-A-ITW, the axial strain at the waist is less than at the ends. As seen from Figs. 15 and 16, when the environmental temperature is maintained at −15 °C at the top and 2 °C at the bottom, the sample temperature gradually decreases. At the same time, the sample is narrowing in a nearly cylindrical fashion. Beginning with WTI, the sample is freezing heaving in the shape of a goblet, starting from the top cooling end through the waist to the bottom end. After the environment temperature is changed to 15 °C at the top and 2 °C at the bottom, the sample temperature gradually rises. Meanwhile, the sample thaws in a shape that is similar to a wooden beer barrel from both ends to the waist and gradually thins down to a bottle shape. Then, the sample expands in a cylindrical shape until the thermal expanding completes after melting. The rock sample is basically horizontally layered and freezes or thaws from the top end to the bottom end. During the freezing process, there are three states of water, water-ice and ice in order from the bottom to the top, and the ice state gradually increases. Eventually, the ice formation depth is 92.5 mm just past the 1# point. In the initial stage of melting, there are six separate layers of water, water-ice, icewater, ice, ice-water, and water in turn from the bottom to the top. Then, it gradually develops into three layers of water, ice-water, and water. The ice completely melted into water at 447 min. The strain duration time of the uniform freeze-thaw test is greater than that of the unidirectional freeze-thaw test before and during the phase change stage for the A3 point and 3# point, but the strain rate is smaller. For C-A-WTI, the radial strain and axial strain of the uniform freeze-thaw test are 180.7 × 10−6 and 159.9 × 10−6, respectively, and are less than those of the unidirectional freeze-thaw test. For W-A-ITW, the radial strain and axial strain of the uniform freeze-thaw test are 35.6 × 10−6 and 34.1 × 10−6, respectively, and are greater than those of the unidirectional freeze-thaw test. Although the freeze-thaw process is similar to the indoor one-way freeze-thaw test of Xia et al. [26], the overall radial strain is greater than the axial strain in the laboratory test. The reason may be related to the temperature propagation path, the height diameter ratio, the size effect of the specimen, and so on. Under the temperature gradient, the strains before and after the phase transition only relate to the expansion and contraction caused by the temperature conduction of the rock sample. The tendency of axial strain to decrease is more prominent of the unidirectional frozen-thaw test, and the maximum difference value is 67.3 × 10−6. Throughout the freeze-thaw cycle, the closer to the top, the larger the temperature gradient is, the faster the temperature change is.

Fig. 13. Axial strain changes with time.

Fig. 14. Radial axial strain and temperature of 3# change with time.

continuous radial strain times of each of the monitoring points are slightly less than the axial direction. The sustaining times of radial strain and axial strain are large at the waist and small at both ends and are 21 min and 24 min, respectively, at the middle of the waist. The continuous strain time of C-B-WTI is greater than the continuous strain time of W-B-ITW. The closer the bottom is, the greater the difference is for each monitoring point. For WTI and ITW, the continuous radial strain times of each of the monitoring points are all greater than that of the axial direction. The continuous time gradually decreases from the bottom to the top of the specimen. The strain continuous time of WTI is not less than that of ITW. The continuous times of radial and axial strain are 39 min and 21 min, respectively, at the top end for WTI, but they are 27 min and 15 min, respectively, for ITW. For WTI and ITW, the radial strain rates of each of the monitoring points are all less than that of the axial direction. The radial and axial strain rates gradually increase from the bottom to the top of the specimen. The radial and axial strain rates are 5.8 × 10−6/min and 17.4 × 10−6/min, respectively, at the top end for WTI but are 7.3 × 10−6/min and 19.8 × 10−6/min, respectively, for ITW. The radial strain rates of WTI are all less than ITW. The axial strain rates of WTI are all less than ITW only at the top

5. Discussion and application For a uniform freezing and thawing test, the change law of temperature is very close to the result of the once freeze-thaw laboratory test by Zhou et al. [22]. The strain duration curve, the strain peak value,

Fig. 15. State of freezing and thawing with time. 10

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Fig. 16. Net deformation with time (m). Note, the red net represents the original grid and the black net represents the freeze-thaw grid with a relative amplification coefficient of 3000. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

and the law that the overall radial strain is greater than the axial strain are all similar to the constant-temperature freezing and thawing indoor test of Kang et al. [25]. Only the strain change rate is greater than the indoor test for WTI and ITW, namely, the numerical simulation is quicker than the indoor test in the phase transformation stage. The strain result of the numerical simulation increases, but an increase is not obvious in the indoor test for W-B-ITW. The simulation test can realize the instantaneous temperature change of the boundary, but the actual indoor test still needs a certain time delay process. Therefore, the whole melting stage completes earlier and more quickly compared with the indoor test. In the whole freeze-thaw cycle, the dry rock sample exhibits the characteristics of heat expansion and cold contraction, and the deformation is elastic without residual deformation [26]. However, the saturated rock sample may experience six stages of cold contraction, frost heaving, cold contraction, heat expansion, melt contraction and heat expansion in turn. Additionally, the deformation is not elastic with obvious residual deformation. Therefore, for the saturated rock sample, the rock matrix still abides by the law of heat expansion and cold contraction. Due to the participation of the phase transition, the properties of cold contraction and cold expansion are presented when the temperature is higher and lower than the freezing point, respectively. The deformation caused by the extending of the pore space volume in the process of water freezing to ice is larger than the deformation caused by the cold contraction of the rock matrix. Thus, the overall deformation performance is a frost heave for WTI. The deformation caused by the closing of some micro-cracks in the process of ice thawing to water is larger than the deformation caused by the heat expansion of the rock matrix. Thus, the overall deformation is a melt condensation process for ITW. The residual deformation is composed of the deformation of thermal expansion and the deformation of unrecoverable shagginess. However, frozen geotechnical engineering, water loses and the rock mass strength reductions after thawing will lead to the intuitive sinking

representation. The simple algorithm model is continuously verified, optimized, promoted and applied from the beginning of the model presented. The change laws of temperature, deformation, stress and moisture on the cylindrical specimen model; the ideal slope model; and the MULI coal mine slope model are researched and described as follows under freezing and thawing conditions. 5.1. Case 1: cylindrical sample Based on this simple model, Li et al. [24] and Li [29] conducted a constant-temperature freeze-thaw test on a cylindrical rock sample that is similar the model presented in to this article. Assuming that the rock sample is saturated, the initial stress field is generated by the dead load of the rock sample, and the initial temperature is 15 °C. The impermeable boundary and the convective temperature boundary are employed. The Z-direction displacement constraint is applied to the bottom of the model only. Li [24] compared the three THM coupling algorithms and analysed the sensitivity of each influencing factor using one freeze-thaw (−15 °C ~ 15 °C) test on a rock sample. The results show that the temperature diachronic curves of the central key node of rock sample nearly are coincident and go through 6 stages (the cooling stage before freezing, the water-to-ice phase transition stage, the cooling stage after freezing, the warming stage before thawing, the iceto-water phase transition stage, and the warming stage after thawing) in turn. The strain duration curves of the side key points of the rock sample are also basically coincident and exhibit 5 stages (cooling contraction, frost expansion, back temperature hysteresis, thaw contraction, and warming expansion) in turn. Furthermore, the greater the convective heat transfer coefficient is, the higher the freezing temperature is, the higher the ambient temperature is, the smaller the porosity is, the shorter the duration time of the water-ice phase transition is and the greater the radial maximum freezing strain is. Additionally, the greater the thermal expansion coefficient and the time of

Table 3 Analysing the temperature and strain of each of the monitoring points. Monitoring point

Radial strain

Axial strain

1#

2#

3#

4#

5#

A3⊥

1#

2#

3#

4#

5#

A3//

Interval time (min) (−15 °C ~ 2 °C) (−15 °C ~ 2 °C) ④ (+15 °C ~ 2 °C) (+15 °C ~ 2 °C)

90.0 114.0 15.0 75.0

51.0 93.0 21.0 63.0

36.0 81.0 21.0 54.0

21.0 57.0 12.0 42.0

9.0 39.0 6.0 27.0

57.0 81.0 45.0 72.0

123.0 78.0 15.0 63.0

78.0 54.0 24.0 60.0

45.0 42.0 24.0 42.0

30.0 27.0 18.0 24.0

12.0 21.0 9.0 15.0

57.0 81.0 45.0 72.0

Rate of change (×10−6/min) (−15 °C ~ 2 °C) (+15 °C ~ 2 °C)

1.2 −1.9

2.3 −3.2

2.9 −3.9

4.1 −5.0

5.8 −7.3

2.7 −2.5

2.9 −3.8

5.8 −3.3

7.8 −5.7

12.7 −11.5

17.4 −19.8

2.4 −2.2

Final stable value (×10−6) ③ (−15 °C ~ 2 °C) ⑥ (+15 °C ~ 2 °C)

141.9 4.8

199.0 7.0

211.7 14.5

207.6 17.2

198.7 28.9

180.7 35.6

255.4 27.9

188.4 5.1

192.9 9.2

187.0 6.6

275.0 59.7

159.9 34.1

Note, the bold in the table indicates the comparative monitoring points for the uniform freeze-thaw test and the unidirectional freeze-thaw test. 11

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Fig. 17. Strata status of open-pit slope in permafrost regions.

freezing throwing cycles are, the greater the freezing strain of radial direction is, etc. Li [29] compared the three field laws of different coupling algorithms using a low-temperature uniaxial compression test and simulated the freeze-thaw cycles by introducing the freeze-thaw coefficient into the coupling program. The results show that the coupling algorithm of the THM series model is feasible and the subsequent temperature change is the same as that of the first freeze-thaw cycle; in addition, the subsequent strain is only accumulated on the basis of the previous strain. These results have some reference significance for designing related experiments about rocks under freezing throwing.

and its three field variation laws, the authors often use the ideal slope model, as shown in Fig. 18 below. Li [29] adopted a simple THM coupling model with water-ice phase change and used the ideal slope model to compare the THM variation laws of slope under different boundary temperature conditions. The authors proposed the relative safety factor based on the yield strain to analyse the factors affecting the stability of the slope. The factors that affect the slope stability under freeze-thaw cycles include the lithology, temperature, moisture, number of freeze-thaw cycles, stress state and geometry of the slope, etc. Among them, the number of freeze-thaw cycles and water migration are the two most important factors affecting the stability of the slope in the rock mass freeze-thaw zone. The height and angle of the slope have the second strongest effects on the slope stability, if regarding the month of the safety factor sudden change as the most dangerous month in a freeze-thaw cycle under the V-type temperature boundary conditions. Due to the effects of global warming, the temperature of permafrost is increasing year by year, and the slope undergoes freezing-thawing temperature alternation. The change in temperature affects the stability and safety of the project. The simplified algorithm model and the yield approach index (YAI) are used, based on FLAC3D and the slope ideal model, the change law on the region of damage or failure in a freezing and thawing cycle and the stability [31]. The volume percentage of YAI on the 12 months of the fifth freeze-thaw cycle is shown in Fig. 19. The result shows that the volume percentage of each month is basically the same. The part with YAI ≤ 0.6 actually is the seasonal frozen soil layer affected by the external environment. YAI ≥ 0.2 represents the structural safety, and YAI ≤ 0.6 represents the effective part for the freezingthawing rock slope. Since the YAI is just the degree of yield or damage on one unit, the local volume (YAI ≤ 0.6) weighted average may be

5.2. Case 2: ideal slope Open-pit mine slopes in permafrost regions have a special rock soil layer distribution [30], as shown in Fig. 17 below. The upper part of the permafrost layer is a seasonal frozen soil layer affected by seasonal temperature changes. After the slope excavation, due to the influence of the periodic external environment temperature, the permanent frozen soil layer within a certain depth of the slope surface will become a seasonal frozen soil layer. Once the excavation depth exceeds the maximum depth limit of the permanent frozen soil layer, the not-frozen layer under the permanent frozen soil layer will also become a seasonal frozen soil layer at a certain depth below the slope surface. Meanwhile, the phenomenon of freezing and thawing alternation will appear. To explore the stability of the open-pit mine slope in the permafrost region

Fig. 18. Ideal slope model.

Fig. 19. Volume percentage changes with YAI distribution. 12

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

conditions. The applicability of the local stability analysis method based on the YAI is discussed. The evolution process of freezing and thawing depth, the tension stress zone, shear strain increment zone and effective yield zone of each month are analysed, as shown in Fig. 21. The seasonal frozen layer of the slope shoulder begins to thaw at the end of March, and the maximum thawing depth reaches approximately 5–6 m from July to September. At the end of September, the temperature drops into negative temperatures. The seasonal frozen layer quickly back freezes from top to bottom. By the beginning of December, the seasonal frozen layer is completely frozen. The depth of freezing and thawing is very similar to the actual investigation. During the seasonal frozen layer thawing, due to the influence of surface thawing sliding, a certain tension stress zone will be generated at the shoulder and the foot of the slope. After the seasonal frozen layer back freezing, the tension stress region will gradually expand and penetrate the permafrost layer. However, the maximum tension stress part moves outward from the slope interior to the shoulder. The deformation zone distributes mainly in the shape of a triangle in the shoulder. The results also show that the new tensile plastic zone will move inward with the refreezing of the seasonal frozen layer. The shallow layer of the slope will have the possibility of sliding along the bottom ice surface of the seasonal frozen layer. This process is similar to the erosion of the slope roof. The interior of the slope is the key to the overall slope safety. If it is reinforced within a certain range of the lower part of the shoulder, the safety of the slope may be greatly improved. The local effective yield zone of the frozen rock soil slope is YYAI ≤ 0.6 and 0.2 < YYAI ≤ 0.6 is the damage zone; 0 < YYAI ≤ 0.2 is the adjacent damage zone; and YYAI ≤ 0.0 is the failure zone. The local effective yield zone combines the properties of the tension-shear failure and its triangular distribution. The failure zones concentrate on the shoulder and the damage zones distribute along the frozen-thaw surface. The effects of the freezing temperature, slope height, number of freeze-thaw cycles, thermal expansion coefficient and porosity on the local effective

Fig.20. YI changes with month.

used to represent the degree of local security of the freeze-thaw slope. Thus, the local yielding approach index (YI) on each month from 4th to 5th times freeze-thaw cycles is analysed, as shown in Fig. 20. The YI decreases rapidly in March, and the lowest month is from April to July. Then, the YI rises with water freezing. That is, the dangerous month is from April to July every year, and in this month, the seasonal frozen soil is most likely to exhibit local collapse or slip failure. The effect degree and the change trend of the influence factors for the slope failure and damage region are different, especially for the failure region on the water content and when the width of the road is extreme. The freezing temperature reduction and the slope height are still the main influencing factors. The porosity, the number of freeze-thaw cycles and the boundary temperature increment still seriously affect the slope stability. The slope surface is susceptible to the periodic impact of the external environment to produce the local slump or the erosion damage [30]. The simple algorithm model and ideal slope model are used to study the local damage failure evolution process and the influence of various factors for the slope stability under the freezing and thawing

Fig. 21. Evolution process on freezing-thawing. 13

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Fig. 22. Contour of shear increment.

yield zone are weakened in turn. Among them, the slope height, freezing temperature, and freezing and thawing times mainly lead to the accumulation of damage inside the slope. The slope angle is more likely to cause damage to the slope surface. The thermal expansion coefficient and porosity will simultaneously increase the damage zone and the failure zone of the slope. Xu [32] used the simplified algorithm model to study the five basic failure modes (bedding layer, cutting layer, reverse layer, fractured stratum and loose formation rock slope) of frozen rock slope under freezing and thawing conditions. Among them, bedding layer slope instability is more common. The main stratification and the slope angle of the bedding slope are between 0 and 45°. The rock layer develops more secondary fractures perpendicular to the main stratification. The fracture gradually penetrates the net under the action of the external force. The participation of water and repeated freezing and thawing intensified the process of crack development, causing the layered rock mass to block and flake off and collapse. When the inclination angles of structural stratification plane are 30° and 60°, the shear strain increment distribution is shown in Fig. 22. When the angle is 30°, the shear strain incremental zone is layered, and the slope has a tendency to undergo slippage failure. When the angle is 60°, there is a significant shear strain increment zone on the shoulder, and the slope may slide collapse at the shoulder. Based on the YAI analysis method, the factors affecting the slope stability under low-temperature freeze-thaw cycles are studied and analysed. The results show that the slope angle, thermal expansion coefficient, initial water content, slope height, number of freeze-thaw cycles, etc., are more likely to cause slope damage and gradually decrease in turn. While the number of freeze-thaw cycles, slope height, boundary temperature increase, thermal expansion coefficient and initial water content are more likely to cause slope failure, the influence degree is gradually reduced in turn. The best method to prevent the slope from being affected by the temperature field and weathering is to isolate the surface. The best method to reduce the crack propagation caused by the frost heave effect of the ice is to reduce the water content in the surrounding rock of the slope. Therefore, slope stability protection measures of the frozen rock edge covering the vegetation layer and the slope drainage were studied and analysed. The results show that both measures can effectively weaken the freeze-thaw damage of the slope. Li et al. [33] used a simplified algorithm model and a similar ideal slope model, based on the YAI analysis method, combined with excavation deformation, to systematically study the excavation key parameters of open-pit slopes in permafrost regions under freeze-thaw cycles. The results show that the slope of the frozen rock soil is mainly local damage failure, and the dangerous month is from May to July every year. Although the local damage index is smaller under cold season excavation, the low-temperature environment is not conducive to construction; therefore, warm season excavation may be used. To some extent, extending the excavation time of the slope could promote the external low temperature fully propagating into the slope, which will enhance the safety of the slope excavation. The excavation depth

equalling the permafrost lower limit could be used as the turning point of the damage mode and deformation mode. When the slope height passes through the permafrost lower limit, the failure mode will be similar to the ordinary slope. Although the deformation increases with the increase in the slope angle, the local damage degree decreases with the increase in the slope angle. Under the premise of ensuring the safety of the slope, the slope angle could be appropriately increased to reduce the extraction ratio. The ideal slope models and calculation parameters are obtained under certain assumptions and simplification conditions. Therefore, the research results have good applicability to homogeneous saturated rock soil slopes, but for the slope of heterogeneous, multi-fracture and thin layers, further research must be performed. 5.3. Case 3: MULI coal mine slope The MULI Coal Mining Area is located in the permafrost zone of the Qinghai-Tibet Plateau, with an altitude of 4000–4200 m. It is one of the largest coal resource zones in Qinghai Province. The seasonal frozen soil layer is approximately 5–6 m, and the permafrost layer is approximately 58–60 m. Due to the seasonal warm and cold changes in the climate, the surface layer of the frozen soil often melts in the warm season to form a seasonal melting layer. The annual minimum temperature is −34 °C, the highest temperature is 19.8 °C, the annual average temperature is −4.2 °C ~ −5.1 °C, and the annual average temperature is −1.0 °C ~ −3.5 °C. Every year in late April, the warm season begins the melting, reaching the maximum depth from the end of September to mid-October. The melting rate varies slightly from month to month, with an average of 0.01 m/d. At the end of September, the temperature begins to fall, entering the negative temperature season. The tops and bottoms quickly returned to freezing, and the melting layers were all frozen in the first quarter of December. The rocky soil of the slope in the permafrost region, especially the Quaternary coverings, are soft rock layers and fractured rock masses with weak interlayers, and they have experienced multiple creeping and hot melt slumps under repeated frost heaving and sedimentation. Coal mining has intensified the shrinkage of the permafrost area, and the scope of the melting area has been expanding. The upper limit of the permafrost has been increasing year by year, the surface water system has been changing year by year, the stability of the slope has been disturbed and destroyed, and the slope collapse is serious. Xu [32] used the simplified algorithm model, based on the evaluation method of YAI, to carry out stability analysis on the frozen rock slope of MULI Open-pit Coal Mine. The typical slope model of Qinghai MULI Open-pit Coal Mine is shown in Fig. 23. The temperature conduction ranges, the actual excavation range and the rock stratum distribution are considered in the model. For the sake of simplicity, the model considers only the influence of the weak structure of the shallow mudstone on the stability of the slope while ignoring the deeper mudstone interlayer. The YAI and the X-direction displacement curve of the slope are shown in Figs. 24 and 25, after the slope is subjected to 14

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

Fig. 26. Temperature changes curve with time on monitoring points of slope shoulder.

Fig. 23. Model and monitoring points of MULI coal mine slope.

Fig. 24. YI changes curve with the number of freezing-thawing cycles.

Fig. 27. Net deformation (m). Note, the black net represents the original grid and the red net represents the deformation grid with a relative amplification coefficient of 1000. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 25. X-direction displacement changes curve with the number of freezingthawing cycles.

multiple freeze-thaw cycles. The results show that June is also the most dangerous month. The YAI decreases with the increase in the number of freeze-thaw cycles, and YAI is nearly steady after 15 freeze-thaw cycles. The displacement of the key point of the lower slope has undergone a large increment; that is, the shallow slope may have a destructive tendency. The temperature diachronic curve of each horizontal key point in the permafrost layer is shown in Fig. 26. After excavation, freeze-thaw areas appear in a certain range of the slope shoulder affected by the ambient temperature. At a depth less than 5 m, the

Fig. 28. YAI distribution.

temperature basically changes with the change in the boundary temperature. At a depth greater than 5 m, the temperature has only a small delayed fluctuation. The temperature does not change with greater depth. After 15 freeze-thaw cycles, the mesh deformation and YAI distribution are shown in Figs. 27 and 28. The low elevation of the slope has a large horizontal deformation to the free surface, especially at the slope foot near the weak structure stratum. The smaller YAI zone 15

Computers and Geotechnics 118 (2020) 103357

G. Li, et al.

distributes mainly in the freezing and thawing range of the shallow slope and slope foot in front of the weak structure layer. There is a possibility of collapse at the top of the slope, and there is a tendency of slip failure at low elevation. Since the height of the slope excavation has penetrated the permafrost layer and is affected by the weak structural plane inside the slope, the whole slope body is the overall shear-slip failure along the structural plane to the foot of the slope. Within a certain range of shoulders and affected by the freeze-thaw cycle, it shows partial slip and collapse damage; in the middle and low elevation slopes and the original unfrozen layer of the slope foot, it is affected by the low-temperature freeze-thaw cycle, and the damage is continuously expanded. Then, surface slip damage may occur. In summary, the algorithm can well express the effect of THM coupling, and the results of temperature and strain are in good agreement with the laboratory test results, which indicates that the simple algorithm with phase change is feasible and has good practicability and applicability. However, due to a lack of space, the sample freeze-thaw cycle characteristics, the sample uniform and non-uniform freeze-thaw characteristics with different shapes and different height diameter ratio, will be discussed in the next paper in detail.

nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled, “A novel simple practical thermal–hydraulic–mechanical (THM) coupling model with water-ice phase change”. Acknowledgements This study is sponsored by the National Natural Science Foundation of China (Grants No. 51779207) and the National Natural Science Foundation of China (Grants No. 11572246). The financial support provided by these sponsors is greatly appreciated. References [1] Zhang M, Lai Y, Wu Q, Yu Q, Zhao T, Pei W, et al. A full-scale field experiment to evaluate the cooling performance of a novel composite embankment in permafrost regions. Int J Heat Mass Transf 2016;95:1047–56. https://doi.org/10.1016/j. ijheatmasstransfer.2015.12.067. [2] Yu F, Qi J, Lai Y, Sivasithamparam N, Yao X, Zhang M, et al. Typical embankment settlement/heave patterns of the Qinghai-Tibet highway in permafrost regions: formation and evolution. Eng Geol 2016;214:147–56. https://doi.org/10.1016/j. enggeo.2016.10.013. [3] Mu Y, Ma W, Wu Q, Sun Z, Liu Y, Qu G. Thermal regime of conventional embankments along the Qinghai-Tibet railway in permafrost regions. Cold Reg Sci Technol 2012;70:123–31. https://doi.org/10.1016/j.coldregions.2011.08.005. [4] Yuan C, Yu Q, You Y, Guo L. Deformation mechanism of an expressway embankment in warm and high ice content permafrost regions. Appl Therm Eng 2017;121:1032–9. https://doi.org/10.1016/j.applthermaleng.2017.04.128. [5] Guo L, Xie Y, Yu Q, You Y, Wang X, Li X. Displacements of tower foundations in permafrost regions along the Qinghai-Tibet power transmission line. Cold Reg Sci Technol 2016;121:187–95. https://doi.org/10.1016/j.coldregions.2015.07.012. [6] Yu W, Lu Y, Han F, Liu Y, Zhang X. Dynamic process of the thermal regime of a permafrost tunnel on Tibetan Plateau. Tunn UndergR Space Technol 2018;71:159–65. https://doi.org/10.1016/j.tust.2017.08.021. [7] Darrow MM, Daanen RP, Gong W. Predicting movement using internal deformation dynamics of a landslide in permafrost. Cold Reg Sci Technol 2017;143:93–104. https://doi.org/10.1016/j.coldregions.2017.09.002. [8] Kang Y, Liu Q, Huang S. A fully coupled thermo-hydro-mechanical model for rock mass under freezing/thawing condition. Cold Reg Sci Technol 2013;95:19–26. https://doi.org/10.1016/j.coldregions.2013.08.002. [9] Tan X, Chen W, Tian H, Cao J. Water flow and heat transport including ice/water phase change in porous media: numerical simulation and application. Cold Reg Sci Technol 2011;68:74–84. https://doi.org/10.1016/j.coldregions.2011.04.004. [10] Wu D, Deng T, Zhao R, Wang Y. THM modeling of ground subsidence induced by excavation of subway tunnel. Comput Geotech 2018;94:1–11. https://doi.org/10. 1016/j.compgeo.2017.08.013. [11] Chen L, Wang J, Liu Y, Collin F, Xie J. Numerical thermo-hydro-mechanical modeling of compacted bentonite in China-mock-up test for deep geological disposal. J Rock Mech Geotech Eng 2012;4:183–92. https://doi.org/10.3724/SP.J.1235.2012. 00183. [12] Zhao Y, Feng Z, Feng Z, Yang D, Liang W. THM (Thermo-hydro-mechanical) coupled mathematical model of fractured media and numerical simulation of a 3D enhanced geothermal system at 573 K and buried depth 6000e7000 M. Energy 2015;82:193–205. https://doi.org/10.1016/j.energy.2015.01.030. [13] Bekele YW, Kyokawa H, Kvarving AM, Kvamsdal T, Nordal S. Isogeometric analysis of THM coupled processes in ground freezing. Comput Geotech 2017;88:129–45. https://doi.org/10.1016/j.compgeo.2017.02.020. [14] Wu D, Zhang Y, Zhao R, Deng T, Zheng Z. A coupled thermal-hydraulic-mechanical application for subway tunnel. Comput Geotech 2017;84:174–82. https://doi.org/ 10.1016/j.compgeo.2016.12.006. [15] Cao W, Huang W, Jiang F. A novel thermal–hydraulic–mechanical model for the enhanced geothermal system heat extraction. Int J Heat Mass Transf 2016;100:661–71. https://doi.org/10.1016/j.ijheatmasstransfer.2016.04.078. [16] Teng T, Zhao Y, Gao F, Wang JG, Wang W. A fully coupled thermo-hydro-mechanical model for heat and gas transfer in thermal stimulation enhanced coal seam gas recovery. Int J Heat Mass Transf 2018;125:866–75. https://doi.org/10.1016/j. ijheatmasstransfer.2018.04.112. [17] Huang S, Liu Q, Cheng A, Liu Y, Liu G. A fully coupled thermo-hydro-mechanical model including the determination of coupling parameters for freezing rock. Int J Rock Mech Min Sci 2018;103:205–14. https://doi.org/10.1016/j.ijrmms.2018.01. 029. [18] Tan X, Chen W, Wu G, Yang J. Numerical simulations of heat transfer with ice–water phase change occurring in porous media and application to a cold-region tunnel. Tunn UndergR Space Technol 2013;38:170–9. https://doi.org/10.1016/j. tust.2013.07.008. [19] Jin Z, Cheng X, Yang Z. Experimental investigation of the successive freezing processes of water droplets on an ice surface. Int J Heat Mass Transf 2017;107:906–15. https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.001. [20] Jin Z, Zhang H, Yang Z. Experimental investigation of the impact and freezing processes of a water droplet on an ice surface. Int J Heat Mass Transf

6. Conclusions A simple practical coupled THM model with water-ice phase change was developed based on the water-ice phase change model and the THM model. The model shows excellent applicability when compared with the historic laboratory test. The main conclusions are as follows: (1) The water-ice phase change, which plays a significant role in the low-temperature multi-field coupling, is proposed based on the energy conservation equation with water-ice phase change and the hypothesis of dividing the energy stages of the freeze-thaw process. (2) The THM coupling model is established by connecting the temperature, seepage and stress in a series. A simple practical model of the THM coupling with phase transition is built by combining the phase change model and the THM coupling model. (3) The combining model is verified by the uniform freeze-thaw test and unidirectional freeze-thaw test. In addition, the results are in good agreement with the laboratory test. For C-B-WTI and W-BITW, the change in temperature and strain are basically synchronous. For WTI and ITW, the continuous time of strain delays by 18 min and 24 min, respectively, on the uniform freezing and thawing test under −15 °C ~ 16 °C conditions. For C-B-WTI and WB-ITW, the continuous times of radial strain of each monitoring points are less than the axial direction. The sustaining time of radial strain and axial strain gradually decreases from the bottom to the top for the former, while the sustaining time is larger at the waist and smaller at both ends for the latter. The favourable practical ability of the developed THM model indicates that it can be applied to explore the multi-field change laws under low-temperature freezing and thawing conditions. Furthermore, it can also be extended as an effective tool for analysing freeze-thaw failure mechanisms and engineering safeguard measures in permafrost regions. Author Contributions Section Li G. put forward and verified the THM coupling model as well as wrote this paper. Li N. provided overall guidance on the article and acted as a correspondent. Bai Y. modified the article language. The remaining authors contributed to data analysis and case studies. Declaration of Competing Interest We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any 16

Computers and Geotechnics 118 (2020) 103357

G. Li, et al. 2017;109:716–24. https://doi.org/10.1016/j.ijheatmasstransfer.2017.02.055. [21] Kim J, Jeon J, Kim DR, Lee K-S. Quantitative analysis of anti-freezing characteristics of superhydrophobic surfaces according to initial ice nuclei formation time and freezing propagation velocity. Int J Heat Mass Transf 2018;126:109–17. https:// doi.org/10.1016/j.ijheatmasstransfer.2018.06.023. [22] Zhou J, Tan L, Wei C, Wei H. Experimental research on freezing temperature and super-cooling temperature of soil. Rock Soil Mech 2015;36:777–85. [23] ITASCA. FLAC 3D (Version 3.0) User’s Manual. Minneapolis, Minnesota, USA: Itasca Consulting Group; 1997. [24] Li G, Li N, Liu N, Zhu C. Practical algorithm of THM coupling process with ice-water phase change based on FLAC3D. Chin J Rock Mech Eng 2017;36:3841–51. https:// doi.org/10.13722/j.cnki.jrme.2016.1643. [25] Kang Y, Liu Q, Zhao J, Zhang F. Research on frost deformation characteristics of rock and simulation of tunnel frost deformatlon in cold region. Chin J Rock Mech Eng 2012;31:2518–26. [26] Xia C, Li Q, Lv Z, Wang Y, Huang M. Comparative experimental study on frost deformation characteristics of saturated rock under uniform freezing and uni-

directional freezing conditions. Chin J Rock Mech Eng 2018;37:274–81. [27] Zhang H-M, Yang G. On the influence of moisture and freeze-thaw cycle on physical and mechanical properties of red sandstone. J Exp Mech 2013;28:635–41. [28] Fang Y, Qiao L, Chen X, Yang SJ. Experimental study of freezing-thawing cycles on sandstone in Yungang grottos. Rock Soil Mech 2014;35:2433–42. [29] Li G. Research on slope stability of open-pit coal mine under freezing-thawing cycles. Shaanxi: Xi'an University of Technology; 2016. [30] Li G, Li N, Liu N, Zhu C. Local stability of open-pit coal mine slope in permafrost regions. J Xi’an Univ Technol 2019;35:53–61. [31] Li G, Li N, Liu N. The stability of open-pit slope in permafrost regions under global warming. In: Zhang L, da Silva BG, Zhao C, editors. Proceedings of GeoShanghai 2018 international conference: rock mechanics and rock engineering, Singapore; 2018. p. 375–87. [32] Xu S. Stability evolution law of open-pit coal mine slope in high-altitude permafrost area. Shaanxi: Xi’an University of Technology; 2017. [33] Li G, Li N, Liu N, Yang M, Zhu C. Study on key parameters for excavation of open-pit slope in permafrost region. Water Resour Hydropower Eng 2018;49:8–17.

17