Numerical analysis of application limits of Vyalov’s formula for an ice-soil wall thickness

Artificial ground freezing is an efficient technique which allows one to sink the mine shafts under complex hydrogeological conditions. The aim of artificial ground freezing is a creation of a temporary wall of frozen soil around the intended excavation. To estimate an optimal thickness of an ice-soil wall, the Vyalov’s formula is widely used by engineers. The article is devoted to analysis of Vyalov’s formula on the basis of the numerical simulation. The numerical simulation has been conducted by the finite element method. For the simulation of a stress-strain state of an ice-soil wall a new computational scheme has been proposed. The scheme is based on Vyalov’s design layout for a vertical shaft sinking and takes into account a soil layer beyond the excavation bottom. A mechanical behavior of frozen soil is described by Vyalov’s constitutive relations. As a result, it has been shown that values of the wall thickness given by Vyalov’s formula do not agree with the ones obtained by the numerical simulation. In order to conform results given by Vyalov’s formula and the numerical simulation, two modifications of the formula have been proposed.


INTRODUCTION
t the present time mineral deposits characterize of an increase of a depth of an occurrence and a complication of hydrogeological conditions of building of a vertical mine shaft for their elaboration. Artificial ground freezing (AGF) is an efficient technique which allows sinking of the mine shafts through weak, unstable, water saturation A soil stratums. The aim of AGF is a creation of a temporary wall of frozen soil around the intended excavation. Ice in soil pores changes its mechanical and hydrodynamical properties, so at certain temperature the ice-soil wall can withstand rock pressure from unfrozen rock mass and eliminate groundwater filtration. The freezing of ground is carried out by a series of wells that are drilled around the excavation. Due to circulation of liquid refrigerant inside the pipes, heat extracts from the near-wellbore domain and the ice-soil retaining structure is established. Strength and reliability of the ice-soil wall determine safety of a shaft sinking and geotechnical works. Mechanical behavior of the wall depends on properties of soil that its formed. Presence of ice and unfrozen water in frozen soil leads to occurrence of rheological properties that characterize development of slow deformations (creep) and losing strength during a prolonged load action [1]. As a result, the ice-soil wall could significantly deform under rock pressure [2]. A lot of experimental studies have been carried out to study time-dependent deformation of frozen soil under uniaxial and multiaxial compression tests [1][2][3][4][5][6][7][8][9][10]. Creep deformation of frozen soil depends on long-term strength limit [3]. If the stress does not exceed the long-term strength limit, deformation rate tends to zero and failure time is infinity. In the case evolution of deformation is characterized as attenuate creep. If the stress is larger than the long-term strength, nonattenuate creep develops. In general, the creep process can be divided on three stage: unsteady creep with a decreasing deformation rate, viscoplastic flow with an almost constant rate, and progressive flow with an increasing rate after that failure arises. In [1,3] it is noted that the long-term strength limit of frozen soil significantly depends on temperature of frozen soil. In [4] from uniaxial creep tests of warm ice-rich sand it has been established that creep characteristics are affected by small temperature variation. In [5,6] on the basis of triaxial tests it has been shown that a decrease of temperature leads to an improvement of strength properties. In [7] a study of an influence of thermal gradient on the frozen sand has been carried out. It has been determined that a decrease of thermal gradient leads to an increase of the long-term strength limit. In [8] it is noted that rheological properties of frozen soils also depend on ice and unfrozen water saturations, grain-size distribution, mineral content. In [9] it has been concluded that cryostructure has an influence on creep behavior of frozen soil. In [10] based on triaxial creep tests it has been shown that coarse-grained content in frozen silty clay causes a rise of it strength. For describing the phenomena of creep of frozen soils various constitutive relations have been developed. To reflect changes in soils that arises during creep process at the microscale into the macroscale, micromechanics [11], thermodynamics [12,13] and damage mechanics [10] are used. However, to determine material parameters for this models difficult and costly experimental studies have to be conducted. As a result, for engineering purposes macroscale phenomenological models are widely used [14]. In [3,14,15] extensive reviews of various models describing creep under complex stress state and various stages of deformation are presented. It could be divided in two types: models that take into account an influence of confining pressure on creep deformation and models in that the influence is supposed to be negligible. In [16,17] Nishihara and Burgers models have been improved to a new creep constitutive models that describe viscoplasticity deformation of frozen soils under high confining pressure. One of basic model for analysis of creep deformation of ice-soil retaining structure is Vyalov's constitutive relations [1,2]. The relations are based on the Norton-Bailey power law in that a modulus depends on time and temperature. An influence of confining pressure on creep of frozen soil is neglected. The advantage of the model is that material parameters for describing deformation in complex stress state can be estimated from simple uniaxial loading tests. Experimental verification of the Vyalov's relations has been performed in uniaxial and multiaxial creep tests in wide range of temperature and ice saturation of frozen soils [1,2,18]. In [19] experimental values of creep deformation of beams of frozen soil have been compared to values obtained by a finite element modeling. In [20] experimental and numerical studies of interaction of footings with frozen sand have been performed. In [21,22] an improvement of Vyalov's relations has been performed by consideration of an influence of volumetric ice content on creep deformation. In [1,23] on the basis of Vyalov's relations an analytical formula for an estimation of an optimal thickness of an ice-soil wall from a condition of maximal creep deformation has been derived. Nowadays Vyalov's formula is widely applied by engineers for design parameters of an ice-soil retaining structure and determination of AGF regimes [24]. However, analytical derivation of the formula is accompanied by a series of strong assumptions that could not agree with real conditions of a shaft sinking. In [25] deformation of an ice-soil wall of Callovian sandy loam has been estimated by numerical modelling with using Vyalov's relations for various depths of a shaft sinking. It has been shown that deformation of an ice-soil wall with thickness estimated by Vyalov's formula exceeds admissible values for depths of the sinking more than 100 meters. The present work is devoted to analysis of Vyalov's formula on the basis of numerical modeling. The simulation of a sinking of a vertical shaft under AGF is performed according to Vyalov's design layout of a shaft sinking [1] modified by consideration of a soil layer beyond the excavation bottom. It is considered three type of soils: chalk, sand and clay that are prevalent on potash deposits in the Petrikov Site in Belarus. Results of the numerical simulation is used for modifications of Vyalov's formula to reconcile values of the wall thickness obtained by the formula to the numerical ones.

MATHEMATICAL MODELS FOR DEFORMATION OF AN ICE-SOIL RETAINING STRUCTURE
n [1] a design layout for a vertical shaft sinking with using AGF has been proposed Fig. 1. The ice-soil retaining structure is a hollow cylinder of inner radius a and outer radius b that encloses unfrozen soil. As a shaft is sunk an inner surface of the cylinder is exposed. The shaft sinking is performed gradually with an installation of a shaft tubbing lining. Therefore, only a height h of the inner surface is unsupported. This part of the surface is deformed under rock pressure p during the time pr t that requests for building of the lining. At that the upper end of the cylinder is supposed to be fixed by the tubbing lining, as it prevents of deformation of the soil to the cylinder center. According to the scheme in [16] the following system of equation in axial-symmetric configuration for describing a stressstate state of an ice-soil wall has been proposed: , , , , where ( , , ) r z  -cylindrical coordinates, u -displacement vector, σ -Cauchy stress tensor, ε -strain tensor, a , binner radius and outer radius of the ice-soil cylinder. The geometry of the ice-soil cylinder is presented in Fig 1. As it can I be seen from (1)-(3) axial stress is supposed to be constant along the height of the cylinder and axial strain is supposed to be constant along the radius of the cylinder. The Vyalov's constitutive relations for frozen soil is generally accepted. The relationship between stress and strain are written as -deformation coefficient depending on temperature T of the frozen soil and time t of loading action, m ,  ,  -material constants. It is assumed that frozen soil is isotropic, homogenous material. The volumetric creep strain of frozen soil is taken to be zero: where tr -trace of a tensor. The hydrostatic stress trσ has no an influence on the strain. A contribution of the elastic strain is neglected. An additional Hencky's relations are taken into the account: From engineering practice and experiments on a physical modeling of an operation of an ice-soil wall it is well known that under an action of rock-pressure radial strain on the inner surface of the cylinder could attain significant values that threaten of failure of freezing pipes and deviation on designed characteristics of the mine shaft. Therefore, the problem of estimation of the ice-soil wall thickness E b a   under that the radial displacement r u on the inner surface do not exceed a limiting value  is stated: On the basis of analytical solution of the problem (1)-(11) taking into account the condition (12) in [16] the following estimation of an optimal thickness E has been obtained: The formula (13) is named Vyalov's formula. To analytically solve the problem (1)-(11) a series of assumptions have been purposed. The key assumption relates to expression for distribution of shear strain rz  on the upper end of the ice-soil cylinder. According to the assumption a distribution of the shear strain is approximated by the following expression: The assumption allows significantly simplifying mathematical treatment and obtaining formula convenient for engineering calculations. Also it should be noted that the problem (1)- (11) has been solved for small neighborhoods of the ends of the ice-soil cylinder.
In order to analysis Vyalov's formula a numerical simulation of stress-strain state of the ice-soil cylinder has been carried out. Vyalov's design layout has been modified by taking into account a soil layer beyond the bottom of the excavation as usual it performs in modelling of vertical shafts (Fig. 2). The layer also consists of the hollow ice-soil cylinder and unfrozen soil inside it. The system of equations considered in the simulation is written as: where grad -the gradient operator, div -the divergence operator. Boundary conditions are identical to (5), (6). The other boundaries of the cylinder are free. The constitutive relations for describing inelastic strain of the ice-soil cylinder under the rock pressure coincide with Vyalov's relations. An additional it is taken into account elastic part of strain. The unfrozen soil it is supposed to be isotropic, homogenous material that can undergo only elastic deformation. The simulation is performed in axial-symmetric configuration with using the finite element method. The geometrical domain is divided on mapped elements of the second order.
The rock pressure p acting on the ice-soil retaining structure is estimated by the standard engineering formula [26]: where l p -effective lateral pressure, w p -hydrostatic pressure, H -depth of bottom of the soil stratum within that the ice-soil cylinder is placed,  l -average density of soils in the rock mass where a vertical shaft is constructed,  -friction angle of unfrozen soil, C -cohesion coefficient of unfrozen soil,  w -density of groundwater, w H -groundwater hydraulic head.

RESULTS OF THE NUMERICAL SIMULATION
he numerical simulation has been conducted for three types of soils: clay, sand and chalk. Elastic    In Fig. 3 -5 distributions of radial displacement r u in an excavation obtained by the numerical simulation are presented.
The thickness S of the ice-soil wall is estimated by Vyalov's formula (13). For the depth up to 300 m Vyalov's formula gives an overvalued estimation of optimal thickness of the ice-soil wall independent on considered soil. In chalk stratum the obtained displacements do not exceed the admissible value for all considered range of the depths. At the depth 500 m the radial displacement of the ice-soil wall of frozen clay exceeds the admissible value by 1 cm. At the same time the displacement of the wall of frozen sand at this depth exceeds the admissible value more than two times. Qualitatively the distributions for soils under consideration coincide. A maximum value of the displacement attains on the inner surface of the ice-soil cylinder. Thus it can be concluded that Vyalov's formula gives an overvalued estimation of the optimal thickness of the wall of quasi-brittle carbonate rocks such as chalk. For more ductile soils such as clay and sand Vyalov's formula gives an overvalued estimation of the optimal thickness of the wall for the depths less than 300 m and an undervalued the one for the large depths. Distributions of shear strain rz  in depending on radius on the top end of the ice-soil cylinder are presented in Fig 6. It can be seen that for soils under consideration values of the strain obtained by the numerical simulation significantly differ from the ones given by the Eqn. (14). Normalization of the distribution on the maximum value it allows one to conclude that there is a qualitative discrepancy between the numerical and analytical distributions (Fig 7).

MODIFICATIONS OF VYALOV'S FORMULA
he results of the numerical simulation presented in the previous section show that for the thickness of the ice-soil wall given by Vyalov's formula the radial displacement either is less than the admissible value or exceeds it. Thus to reconcile Vyalov's formula to the results of the numerical simulation it has to be modified. From the plots of the shear strain obtained by the numerical simulation (Fig. 6, Fig. 7) it can be seen that the all curves converge to a limit curve under the increasing of the depth. As a consequence, it can be concluded that at the large depth H the shear strain rz  on the top end of the ice-soil cylinder do not reflect a relationship between a maximal value of the radial displacement r u of the inner surface of the cylinder and rock pressure p . Thus a modification of Eqn. (14) does not allow one to modify Vyalov's formula that it coincide with the result of the numerical modeling. Fig. 8 shows estimations of the thickness E of the ice-soil wall obtained by Vyalov's formula and the numerical simulation on the basis of the condition (12) for clay, sand and chalk stratums. The rock pressure p corresponds to the depth H within the range from 100 m to 500 m. By analysis of the presented data it has been noted that for a certain value of the rock pressure Vyalov's formula can be conformed to the results of the numerical simulation by a simple way: T where Vyalov E -a value of the thickness given by Vyalov's formula (13), 1 E -a modified value.
In Fig. 9 (a)  A value Vyalov E depends on a value of the rock pressure acting on the ice-soil wall. Therefore, it can be deduced that the Eqn. (20) is correct within a range of the rock pressure. In order to describe the dependency of the thickness of the icesoil wall obtained by the numerical simulation within the all range of rock pressure the following modification of Vyalov's formula has been proposed: Here the function  ( ) g g p is a quadratic function where coefficients i  , ( 1,2,3) i  , are depended on considered frozen soil. In Tab. 4 values of the coefficients for clay, sand and chalk are listed. Fig. 9 (b) shows estimations of the thickness of the ice-soil wall of frozen clay, sand and chalk given by the Eqn. (21). It can be seen that the proposed equation allows one to describe the results of the numerical simulation both qualitatively and quantitatively.

CONCLUSIONS
his work is devoted to analysis of Vyalov's formula for determination of an optimal thickness of an ice-soil retaining structure on the basis of a numerical simulation. For the numerical simulation of a stress-strain state of the ice-soil wall a new computational scheme has been proposed. The scheme is based on Vyalov's design layout for a vertical shaft sinking and takes into account a soil layer beyond the excavation bottom. The layer consists of a hollow ice-soil cylinder and unfrozen soil inside it. A mechanical behavior of frozen soil is described in according to Vyalov's constitutive relations. In comparison with an analytical approach used for derivation of the formula in the numerical simulation it has not to make additional assumptions concerning to distributions of stress and strain. As a result of the numerical simulation it has been shown that if a depth of soil stratum occurrence is less than 300 m, Vyalov's formula gives an overvalued estimation of the wall thickness for the all soils under consideration. For chalk, radial displacement of the wall is less than an admissible value within all range of the studied depths. For sand and clay, the displacement exceeds the admissible value if the depth is large than 300 m. Thus, in quasi-brittle carbonate rock stratum an ice-soil retaining structure with a wall thickness estimated by Vyalov's formula has excess margin of safety. In ductile soil stratums, that are occurrence at large depth, Vyalov's formula gives an understated estimation of the wall thickness. It has been studied of an analytical expression for distribution of shear strain on the top end of the cylinder. For soils and depths under consideration the distributions given by the expression do not coincide qualitatively and quantitatively with the ones obtained by the numerical simulation. Since Vyalov's formula is derived on the basis of the analytical expression, this is a reason due that estimates of the wall thickness obtained by Vyalov's formula do not agree to the numerical results. In order to conform estimations of the wall thickness given by the Vyalov's formula to the results of the numerical simulations, two modifications of the formula have been proposed. The first modification could be useful for engineering computations. According to the modification if a value of a thickness of the ice-soil wall given by Vyalov's formula is less than 11 m than it has to be divided on 2 to coincide with the value given by the numerical simulation. This rule is fulfilled for the all studied soils. The second modification intends including in Vyalov's formula a quadratic function on a rock pressure that acts on an ice-soil wall. Coefficients of the function depends on soil from that consists of the ice-soil wall. The modification allows one to describe the results of the numerical simulation qualitatively and quantitatively within all range of the rock pressure. However, to determine the coefficients of the function the numerical simulations has to be conducted for typical stratums of an elaborated deposit. Thus, it can be concluded that for soil stratums occurrence at small depths a use of Vyalov's formula in design computations causes overvalued costs on the AGF process. To prevent dangerous breakdown at large depths of a shaft sinking, estimates given by Vyalov's formula should be corrected by the numerical simulation.