Theoretical analysis of frozen wall dynamics during transition to ice holding stage

Series of calculations for the artificial freezing of the rock mass during construction of mineshafts for the conditions of a potash mine in development was carried out. Numerical solution was obtained through the finite element method using ANSYS software package. Numerical dependencies of frozen wall thickness on time in the ice growing stage and ice holding stage are obtained for two layers of the rock mass with different thermophysical properties. External and internal ice wall boundaries were calculated in two ways: by the actual freezing temperature of pore water and by the temperature of –8 °С, at which laboratory measurements of frozen rocks' strength were carried out. Normal operation mode of the freezing station, as well as the emergency mode, associated with the failure of one of the freezing columns, are considered. Dependence of a decrease in frozen wall thickness in the ice holding stage on the duration of the ice growing stage was studied. It was determined that in emergency operation mode of the freezing system, frozen wall thickness by the –8 °C isotherm can decrease by more than 1.5 m. In this case frozen wall thickness by the isotherm of actual freezing of water almost always maintains positive dynamics. It is shown that when analyzing frozen wall thickness using the isotherm of actual freezing of pore water, it is not possible to assess the danger of emergency situations associated with the failure of freezing columns.

pressure is made with the assumption of a uniform temperature distribution in the entire FW volume [1,12]. This temperature is usually lower than the actual freezing temperature of pore water in the rock mass.
The first method for selecting boundary isotherms assumes that the average volume temperature for the FW is approximately equal to the temperature, for which the rock strength and creep parameters was calculated. Moreover, near FW boundaries, temperature of frozen rocks is higher than the temperature in calculations of rocks' strength and creep parameters. The second method considers the temperature at FW boundaries to be equal to the temperature in calculations of rocks' strength and creep parameters. In this case, average volume temperature for the FW is lower than the temperature, at which the rocks' strength and creep parameters was calculated. In this regard, the first method can lead to an underestimation of the actual FW strength and creep parameters characteristics, and the second method -to overstatement. Deeper analysis of the advantages and disadvantages of the two indicated methods for selecting the boundary isotherms of the FW should be carried out by means of conjugate thermophysical and geomechanical calculations, which is beyond the scope of the thermophysical study described in this paper.
Based on the experience described in the literature on artificial freezing of rocks during mineshafts construction [2,10,11], it can be accepted that both methods of selecting boundary isotherms of the FW are appropriate in the ice growing stage of the rock mass when the temperature of the freezing brine (coolant) is minimal, and its flow rate is maximum. However, during transition to ice holding stage, which corresponds to higher brine temperatures and/or lower brine flow rates, relevance of these two methods is not obvious. This is primarily due to the fact that, upon transition to ice holding stage, temperature field near freezing columns can vary greatly [9]. In turn, this can lead to a local decrease in FW thickness to values lower than those required by the condition of strength and creep.
Emergency situations during the operation of artificial freezing systems raise even more questions about relevance of the mentioned above methods. Emergency mode can be associated, for example, with the failure of one or more freezing columns [4,8]. In this case, temperature field of the frozen rock mass in the vicinity of the damaged freezing columns can change strongly.
In the Russian and foreign literature on mining thermophysics, the issue of reducing FW thickness, calculated by different isotherms during the transition to ice holding stage in case of emergency situations during the operation of the artificial freezing system has not been studied. Moreover, this issue is relevant for performing the correct thermophysical calculations for freezing of the rock mass and FW formation.
Aim of this study is a theoretical analysis of the dynamics for FW thickness, calculated by various isotherms during the transition from to ice growing stage to ice holding stages. Particular attention is paid to the analysis of FW thickness in the event of an emergency situation associated with the failure of the freezing column.
Mathematical formulation of the problem. Horizontal layer of the rock mass was considered in the work. To carry out a theoretical analysis of heat transfer in a rock mass layer during its artificial freezing, a two-dimensional unsteady heat conduction case with a moving phase transition boundary was solved. In the model of the rock mass, following physical processes were taken into account that affect the heat distribution in it: "water-ice" phase transition, conductive heat transfer (heat conduction) in the rock mass, heat transfer between the rock mass and the circulating coolant in the columns.
When stating a mathematical problem, following assumptions are made: 1) rock mass has isotropic and homogeneous thermophysical properties in frozen and unfrozen zones; 2) vertical heat transfer is negligible compared to horizontal one; 3) temperature of the brine in the freezing columns is constant; 4) freezing columns are equidistant and oriented vertically; deviations of column positions from a perfectly vertical project direction are not considered.
The second assumption is applicable only if the considered rock mass layer has a sufficiently large thickness (more than 10 m), and the simulation time interval is not very large (less than 200 days).
Taking into account the third and fourth assumptions, rotational symmetry appears in the case. Therefore, it suffices to consider the sector of the rock mass, limited by the two main FW planes [11]. Figure 2 shows a simplified geometric model of the frozen rock mass layer, which was used for numerical calculations.
In the geometric model (Fig.2), boundary zone I represents the internal boundary of the calculated area. This zone was introduced to exclude an acute angle from the model and ensure the stability of the numerical solution. Zone I should be offset to the axis of rotational symmetry as much as possible.
Mathematical formulation of the heat conduction problem with a moving boundary of the phase transition (or the Stefan problem) in the enthalpy form is written as follows: where H -specific enthalpy of water-saturated rock mass, J/m 3 ; x, y -physical coordinates, m; tphysical time, s; λ , λ lq sd -thermoconductivity of the mass in the unfrozen and frozen zones, respectively, W/(m·°C); , lq sd c c -specific heat capacity of the mass in the unfrozen and frozen zones, respectively, J/(kg·°C); ρ , ρ lq sd -density of the mass in the unfrozen and frozen zones, respectively, kg/m 3 ; T lq -onset temperature of pore water crystallization (liquidus temperature), °С; T sd -melting temperature of pore ice (solidus temperature), °С; φ ice -volume fraction of ice in the pores of the rock mass, m 3 /m 3 ; L -specific heat of water crystallization, J/kg; n -porosity of the mass;  wdensity of pore water, kg/m 3 ; ( ) fb T t -temperature of the walls in freezing well, °С; T 0 -temperature of undisturbed rock mass at a distance from the freezing circuit, °С;  -heat transfer co- Phase transformations of pore water in the model are taken into account by setting the specific enthalpy function H(T) (3). Surge of this function in the vicinity of the phase transition temperatures is determined by the latent heat of the phase transition in a unit volume of the water-saturated rock mass. Typical piecelinear form of the specific enthalpy function H(T) used in this work for numerical calculations is shown in Fig.3. Calculation of coefficient for the heat transfer between the rock mass and the freezing brine in (4) is determined by the formula for the transition mode of the brine flow [7]: where d -internal diameter of the pipeline, m;  -characteristic cross-section of the flow region of the brine, which exchanges heat with the rock mass, m; λ c -brine thermal conductivity, W/(m•°С);  с -brine density, kg/m 3 ; V c -brine velocity, averaged by cross section, m/s; η с -dynamic viscosity of the brine, Pa•s. Structurally, freezing columns are made of two pipes -an external freezing and an internal supplying. Freezing brine first moves down the internal flow region, and then up the external flow region between the outer and inner tubes. Ratio of the cross-section areas for the external and internal flow regions to simplify the calculation can be taken equal to 1. In this case, characteristic crosssection of the flow region of the brine is  = 2 d . Volume fraction of ice in the pores of the rock mass is determined as follows: Emergency operation mode of the freezing system is modeled by turning off (or failure) one of the freezing columns at the moment of transition ice holding stage. Shutdown is modeled by setting  = 0 at the moment of transition ice holding stage.
Such a choice for the moment of system's emergency operation occurrence is due to the following considerations. Firstly, occurrence of malfunctions in the freezing system seems to be most likely in system's transient modes of operation. Secondly, in this case it will be possible to carry out the most informative and understandable assessment of the decrease in FW thickness, since the comparison will be made with the design value of FW thickness, achieved in the ice growing stage in the normal mode of operation of all freezing columns.
Numerical simulation. Numerical solution for the case of FW formation is carried out using the finite element method in the Transient Thermal module of the ANSYS software package.  As input data for numerical calculation, initial data for the project of the rock mass freezing in the shafts of a potash mine under construction in the Republic of Belarus were used.
Two rock layers were studied: the least heat conducting (chalk layer) and the most heat conducting (clay layer) in the considered range of rock freezing (0-275 m). Main thermophysical properties of the two considered rock layers are presented in Table 1: It was assumed that freezing brine with a temperature of -35 °С is supplied to the columns during ice growing stage. During transition to ice holding stage, temperature of the brine for 5 days evenly rises to -20 °C, brine flow rate is constant and taken equal to 240 m 3 /h. FW thickness was determined in two ways -by the actual temperature T sd of complete freezing of water in the pore space of the mass and by the isotherm T d = -8 °С, at which the strength and creep parameters of the rock was determined. Actual freezing temperature of pore water for the clay layer is T sd = -0.9 °С, while for the chalk layer T sd = -0.58 °C.
Radius of the external boundary for the calculated area was taken to be 40 m, whereas the radius of the internal boundary was 0.25 m. Radius of the freezing column contour was 8 m. Outer radius of the freezing column was 0.073 m and inner radius was 0.068 m. Distance between the centers of neighboring freezing columns was 1.2 m.
Transition to ice holding stage varies in the range of 12.5-100 days with an increment of 12.5 days. Figure 4 shows the dependencies of FW thickness E on time obtained as a result of numerical simulation. Numerical dependencies are presented for two layers (chalk and clay), two operating modes (normal and emergency), two methods for determining FW thickness (by isotherms T sd and T d ). Fig.4, a and c shows the transition to emergency mode to occur after 50 days, whereas Fig.4, b and d -after 100 days since the start of freezing.
Chalk freezing is much slower compared to clay (Fig.4). This is due to the lower thermal conductivity of the chalk and the high water content in it. Differences between the graphs of FW thickness E(t) in the ice growing and ice holding stages, as well as in the normal and emergency operation modes of the freezing station, are most significant for the case of calculating FW thickness by the isotherm T d = -8 С. This is due to the fact that this isotherm is closer to the contour of the freezing columns; therefore, a change in the operating mode of the columns has a stronger effect on it than on the more remote T sd isotherm. Also, when the operating mode of the columns changes, temperature change near the boundary of the actual phase transition of pore water (i.e. near T sd isotherm) occurs more slowly due to the additional inertia of the temperature field in this zone, associated with the duration of the phase transition process (which is determined by the latent heat of the phase transition and water content in the rock).
Inertia of the phase transition process can be easily demonstrated by writing the equation for motion of the phase transition front, which for one-dimensional radial Stefan problem [3,6] has the form: where R(t) -radial coordinate of the phase transition front, m; q lq -heat flow from the unfrozen zone, W/ m 2 ; q sd -heat flow to frozen zone, W/m 2 . Fig.4 shows that for all layers and all modes of freezing system operation, during transition to ice holding stage, a decrease in FW thickness to values lower than required by the project is observed. With an earlier transition to ice holding stage (50 days), a decrease in FW thickness is manifested more significantly than in the case of a later transition (100 days). This is explained by the fact that after 100 days of ice growing stage the rock mass manages to cool to lower temperatures and to a greater distance than for 50 days, as a result of which the heat influx, coming from the unfrozen zone to the phase transition front, decreases. In emergency mode of operation of the freezing system, decrease in FW thickness by the isotherm T sd = -8 С manifests itself more significantly.
FW thickness reduction analysis. An informative parameter that allows an analysis of the decrease in FW thickness in the ice holding stage and emergency operation of the system is the maximum decrease in FW thickness: where E A -design FW thickness achieved by the end of the ice growing stage, m; E P -minimal FW thickness in the ice holding stage, m. Table 2 presents the values of FW thicknesses and their maximum decrease for the two layers under consideration by the isotherm -8 °С.  Figure 5 shows the numerical dependence of E by the -8 °C isotherm on time of ice growing stage of the rock mass for emergency operation mode of the freezing system. Schedule for the normal operation mode of the freezing system is not given, since E is almost everywhere equal to zero: value is different E from zero only in a small time interval from 35 to 50 days and does not exceed 0.15 m.
Figure shows that in the emergency operation mode of the system, E is a complex nonmonotonic function of ice growing stage time. Moreover, for the two rock layers considered in this work, this function has a significantly different form.
Growth curve E in the interval from zero to the first local maximum corresponds to a situation in which failure of the freezing column leads to a decrease in FW thickness by the -8 °C isotherm to zero. With a further increase in time of ice growing stage, function ( ) E t  for chalk decreases monotonously, tending to a zero value. Decrease of ( ) E t  in this case is due to the fact that the longer the time of ice growing stage, the stronger the mass has time to freeze in the unfrozen zone and the smaller the decrease in FW thickness in ice holding stage will be.
The ( ) E t  curve for clay decreases only until time instant of 37.5 days, after which the decrease is replaced by an increase. This behavior of ( ) E t  curve for clay is associated with different dynamics of the internal and external fronts of the phase transition (Fig.6). In the interval of 12.5-37.5 days, the decrease in ( ) E t  for clay is due to the fact that the internal FW front moves deeper into the mass faster than the displacement of the external FW front to the contour of freezing wells. Further, displacement rate of the external FW front begins to exceed the advancement speed of the internal FW front, as a result of which, in time interval of 37.5-62.5 days the function ( ) E t  increases. After 62.5 days the displacement of the external FW front again becomes smaller than the advancement of the internal FW front, therefore ( ) E t  decreases and tends to zero. Analysis of Fig.5 allows to conclude that the decrease in FW thickness by the -8 °C isotherm in ice holding stage in emergency operation mode connected with the failure of one of the freezing columns is significantly higher than zero, i.e. over some period of time, actual FW thicknesses will be lower than the values required by the project, which is unacceptable. When analyzing FW thickness by isotherms of actual freezing of rocks, it is impossible to "catch" such a decrease in FW thickness in ice holding stage. Moreover, in practice, there may be a situation when FW thickness in the isotherm T sd will continue to increase, and FW thickness by the isotherm -8 °С will decrease to zero. In this case, it is advisable to develop methods and measures aimed at eliminating the situation when the actual FW thickness will be lower than the design one. One of such measures may be a temporary change in the parameters of the freezing brine in ice holding stage (decrease in temperature, increase in flow rate). In addition, a possible method to increase the reliability of the artificial freezing system in the event of emergency shutdown of freezing columns may be to regulate the distance between adjacent freezing columns at the stage of developing a project for freezing a rock e -chalk; f -clay, 87.5 days of ice growing stage