Introduction
Taber (1929; 1930) first recognized that frost heave requires the flow of unfrozen water towards freezing front and is different from the pore water expansion due to freezing. The suction force that moves the unfrozen water is termed as cryogenic suction. Radd and Oerstel (1973) proved that the cryogenic suction is governed by the generalized Clausius-Clapeyron equation. A lot of frost heave theories and models were developed applying the generalized Clausius-Clapeyron equation. For instance, the secondary frost heave theory (Miller 1978) explained the mechanics of discrete ice lens formation based on well-established soil physics, and further was modeled as the rigid ice lens model by O'Neill and Miller (1985).
A complete analysis of the multi-dimensional frost heave problems must be able to deal with the coupling of heat transfer, moisture transfer, and mechanical analysis. It is the rigid ice lens model, which has received the most remarkable credit in microscopic approach such as simulating an individual ice lens. However, to acquire quantitative results, especially for multidimensional frost heave problems in frozen ground engineering, the microscopic approach of the rigid ice lens model may not be the most reliable answer. This is often true in field conditions, in which empirical models based on laboratory frost heave test results are far more useful for predicting frost heave and induced load in multi-dimension than the microscopic models, which require a lot more of the complicated input parameters. A good example of the empirical models is the Segregation Potential (SP) concept.
The SP concept prevails because of its considerable simplification. The concept simplifies the frost heave theory with the notion that the heat and moisture transfer is considered a coupled process, on which the mechanical analysis is independent. Consequently, the pressure acting on the frozen fringe is simply equal to the overburden pressure using a non-deforming mesh (e.g. Kim et al. 2008; Konrad and Morgenstern 1984). However, the SP frost heave models using the non-deforming mesh can not respond to an induced stress due to frost heave displacement. Therefore, a consecutive frost heave law involving thermal, mass transfer, and stress coupling is necessary to understand the multi-dimensional pressure effect on SP using deforming mesh.
This chapter has two objectives. The first is to develop a consecutive frost heave model applying the SP concept in finite element application with deforming mesh. Also, it was found that in-situ Fairbanks silt had dependency upon the cooling rate and applied pressure as presented in Chapter 3. The other objective is to verify the SP values, which were obtained by the proposed fitting method using the developed frost heave model.
In the following sections, the fundamental function used in the developed frost heave model is described first. A description of the heat transfer model and the basic assumptions and implementation procedure of the developed frost heave model follow. Finally, the verification of the SP values with the effect of cooling rate and applied pressure and the developed frost heave model against a laboratory frost heave test is presented.
Porosity growth function applying the segregation potential concept
As far as the author knows, Blanchard and Fredmen (1985) were the first individuals to propose a multi-dimensional frost heave model for coupling heat transfer, moisture transfer, and mechanical analysis. Simulating individual ice lens was not considered in their model; instead, the ice lens growth was modeled as volume increase of the soil-ice-water mixture.
The idea that has been accepted by many frost heave studies includes control volume finite element method (Coutts 1991), porosity rate function (Michalowski 1993), and modified porosity rate function (Michalowski and Zhu 2006). These earlier work is expanded in this study.
First, volumetric expansion zone in this model is considered. As discussed in the previous chapters, frost heave can be described as a moisture transfer problem to a growing ice lens past the layered frozen fringe and the unfrozen soil. From an engineering point of view, the zone that is colder than the coldest growing ice lens is defined as a passive system. Moisture migrations were observed under temperature gradients in the passive system (e.g. Hoeskstra and Chamberlain 1964; Ishizaki et al. 1985). However, moisture transfer was much reduced in the passive system due in part to very low hydraulic conductivity of the frozen soil, and consequently the contribution to the total heave is negligible in both laboratory conditions (Mageau and Morgenstern 1980) and field condition (Slusarchuk et al. 1978). Therefore, volume expansion is considered only in the active system in this study. The active system is described in Figure 4.1. Schematically, the active system is composed of growing ice lens and frozen fringe. As the active system is assumed fully saturated, the active system is modeled as a mixture of soil-ice-water.
Next, the active system has two volumetric expansion zones: namely the in-situ freezing zone and segregation freezing zone (see Figure 4.1). The expansion due to in-situ heave occurs in the in-situ freezing zone, which is in between freezing temperature (T0) and in-situ freezing temperature (Tin). The in-situ heave is to account for the in-situ pore water expansion due to freezing. Therefore, it is directly related to the fraction of unfrozen water in the in-situ freezing zone. With accurate water intake rate measurements in a series of step-freezing tests shown in Chapter 2, the portion of unfrozen water is directly determined from a series of step-freezing test results as follows:
[4.0]
where ï¤ (t) = fraction of unfrozen water which changes in in-situ freezing at time t; hin (t) = in-situ heave at time t; nini = initial porosity; and X0(t) = depth of freezing front penetration at time t.
The maximum value of ï¤ was 0.9 in a series of step-freezing tests for the in-situ Fairbanks silt. It means that 90% of pore water freezing is the uttermost amount during freezing process. Assuming all in-situ heave occurred within the in-situ freezing zone, therefore, the unfrozen water contents (w) - temperature (T) relation is described as:
[4.0]
where w0 = unfrozen water content at freezing point; wu = residual unfrozen water content; and Tin = temperature at the end of in-situ freezing.
Figure 4.2 shows the comparison between the calculated unfrozen water content curve and the measured results of Fairbanks silt at the UAF frost heave experiment site (Huang et al 2004).
The segregation freezing zone is between a location in frozen fringe (Tsp) and the warmest growing ice lens temperature (Ts) as shown in Figure 4.1. Excess pore ice will occur in the segregation freezing zone due to the migrating water. Since it is assumed that all in-situ heave has occurred within the in-situ freezing zone, the mass of soil and unfrozen water are constant in the segregation freezing zone. Therefore, the incremental change of porosity in the segregation freezing zone is identical to the incremental change in volumetric ice content in response to the migrating water.
The incremental change in total porosity (ï„nt) is composed of porosity growth due to in-situ heave and segregation heave. Defining the initial water content at time t is w(t), the porosity growth due to in-situ heave (ï„nin) at time t+ï„t can be expressed as:
[4.0]
where ï±w(t) = volumetric fraction of water at time t; and ï¤max = maximum value of fraction of unfrozen water (e.g. 0.9 for in-situ Fairbanks silt).
To predict segregation heave, the SP concept is applied to the developed frost heave model. At each time step, the water migrating rate (v) is calculated as:
[4.0]
where SP = segregation potential; ï³ov = the stress acting on the segregation freezing zone; Tsp = cooling rate of the segregation freezing zone; P0 = pore water pressure at the freezing front; and gradTsp = temperature gradient in the segregation freezing zone.
The amount of migrating water is converted to volume strain. Thermal solution gives the effective area of segregation freezing zone (Asp) and determines the effective width of the segregation freezing zone in the element (lsp). Thus, the porosity growth due to segregation heave (ï„nsp) is calculated as:
[4.0]
whereï€ ï€ ï„nsp is named as SP porosity growth function here.
The incremental change in total porosity (ï„nt) is identical to the incremental change in volumetric strain. Assuming the initial porosity at time step t is n(t), the new porosity at time n(t+ï„t) is updated as:
[4.0]
The time-dependent volumetric fractions (ï±ï€ ) of the three-phase relations were defined as shown in Figure 4.3 and calculated as below:
[4.0]
where ï² = density; and subscripts, s, w, and i = soil, water, and ice, respectively.
Basic equation for heat transfer
The three basic modes of heat transfer are conduction, convection, and radiation. Although heat convection was considered in some models (e.g. Harlan 1973), it is usually negligible in soil heat transfer studies (Farouki 1981). Nixon (1975) reported that the convection component is two or three orders of magnitude smaller than the conduction component. Radiation components hardly contribute to the amount of heat transfer at all. Radiation accounts for less than 1% of the total heat transfer in sands and even less in fine-grained soils (Farouki 1981).
For the above mentioned reason, the developed frost heave model with isotropic thermal properties is based on conduction only and is governed by the following equation:
[4.0]
where C = volumetric heat capacity; L = volumetric latent heat; T = temperature; and ï¬ = thermal conductivity.
The thermal conductivity of soil mixture is determined by Johansen's method (Johansen 1977). The volumetric heat capacity of soil mixture is calculated as the sum of the volumetric heat capacities of the four phases multiplied by their volumetric fractions (Kay et al. 1977) as below:
[4.0]
where; ï² = density; c = mass heat capacity; ï± = volumetric fraction; and subscripts, s, w, i, a = soil, water, ice, and air, respectively.
Alternatively, eq. [4.0] can be expressed as:
[4.0]
where Cap = apparent volumetric heat capacity.
In this study, unfrozen water effect is directly taken into account in the phase change (e.g. Nixon 1983; Shen 1988) as:
[4.0]
where Lw= latent heat of fusion per unit mass of water.
Basic assumptions for the developed frost heave model
In order to solve this complex four-phase thermal transfer problem, several basic assumptions are made in this modeling study. They are as follows:
The soil is fully saturated (air does not exist);
Pore water pressure at the freezing front is zero;
Moisture migration occurs only in the active system and unfrozen soil. The moisture migration in the frozen zone on the cold side of the active ice lens is negligible;
Moisture migration occurs only by the liquid water form. The air phase and vapor transportation are not considered;
Effect of salt exclusion is negligible;
Effect of consolidation in unfrozen soil due to frost heave is negligible;
Volume of soil particles remains constant in freezing process;
Effect of freeze-thaw cycle is negligible; and
Segregation-freezing point depression under loading is negligible.
Implementation of the frost heave model
The developed frost heave model applying the SP porosity growth function is implemented by ABAQUS ver.6.8.1 finite element code. ABAQUS/CAE and ABAQUS/Standard are used for modeling. The ABAQUS/CAE is the user interface and the ABAQUS/Standard is the simulator. The ABAQUS/CAE creates an input file. The input file is composed of the following tasks: generation of finite element mesh, assignment of material properties, and definition of boundary condition. The created input file is submitted to the ABAQUS/Standard. The results are visualized in ABAQUS/CAE.
Also, ABAQUS/Standard has a user-defined subroutine function to extend its capabilities of the applications. The user-defined subroutine is called at all integration points of elements at each time step, and then is linked to the ABAQUS/Standard implementation. The user-defined subroutine is composed of four components - SDVINI, HETVAL, USDFLD, and UXPAN. The user-defined subroutines are written in FORTRAN language.
The flow chart of the user-defined subroutine is shown in Figure 4.4. First, SDVINI is used to define the initial material properties. For any particular time step, temperature distribution is first calculated using the user-defined subroutine HETVAL with an iterative procedure. The thermal results are used to determine the location of the active system. Once the location of the active system is determined, the change of unfrozen water content and gradTsp are calculated. Next, stress state is called from the previous time step using USDFLD. The user-defined function UEXPAN is then used to calculate the volumetric increase of ice due to both in-situ and segregation heave. The volumetric increase is passed to the next time step to carry out the thermal analysis.
Validation of the developed frost heave model
The developed frost heave model is verified against a series of step-freezing tests using in-situ Fairbanks silt. The soil sample is 100mm diameter and 115mm initial height. The boundary conditions are shown in Figure 4.5. Also, the initial and boundary conditions for the tests are summarized in Table 4.1.
A two-dimensional finite element geometry is generated to simulate a series of step-freezing tests. The geometry is divided by 1mm square elements. The mesh is only one element wide. The vertical boundaries are adiabatic and no displacement is allowed in the horizontal direction. Therefore, this two-dimensional finite element analysis is treated as a one-dimensional thermal-displacement analysis.
The soil sample is assumed fully saturated. The initial gravimetric water content is 0.325. The parameters for unfrozen water characteristics curve (eq. [4.0]) are defined as: w0 = 0.325 and wu = 0.0325 (= 0.1 x w0), respectively. The parameters for the frozen fringe are defined as: T0 = 0 oC, Tin = -0.1 oC, Tsp = -0.15 oC, and Ts = -0.35oC, respectively. Volumetric heat capacities are calculated by eq. [4.0]. The values of thermal properties are listed in Table 4.2. Thermal conductivity of the soil sample is 1.19W/(m x oC) for unfrozen state and 2.11W/(m x oC) for frozen state. The thermal conductivities were obtained by the thermal conductivity needle-probe method (Kim 2003), then calibrated using the correlations of Johansen (1977). The initial dry density and initial porosity of soil are calculated as 1451kg/m3 and 0.472, respectively.
The SP values, which are obtained by the proposed fitting method in Chapter 3, are utilized for the simulations. The SP values are calculated by:
[4.0]
where SPm = the maximum SP value at the peak; Tm = the cooling rate at which the maximum SP value occurs, T0 = the cooling rate to satisfy the SP at the formation of ice lens; and CRC = the critical rate of cooling, at which water intake starts. Figure 3.7 shows the relationship between the SP values and cooling rate at different applied pressure. The parameters for eq. ] are summarized in Table 4.3.
In these simulations, it is assumed that all porosity grows predominantly in the direction of the one-dimensional heat flow and no-overburden is applied instead of using the SP value at certain applied pressure. Therefore, the total strain increment in this simulation consists of the porosity growth only. However, examining multi-dimensional frost heave problem, the constraining of the deformations in orthogonal directions will influence the extent of frost heave development in the predominant direction. Additionally, the constraint will produce stress. Detail aspects of the anisotropy and mechanical analysis coupling will be presented in the next chapter.
The step-freezing test result at a 50kPa overburden pressure, STEP-4, is shown as an example of simulation. The initial temperature has a constant temperature gradient, which is 1oC at the bottom and 2.2oC at the top.
Figure 4.6a, b, and c show the variation of temperature, ice content, and porosity distributions in the soil sample at different time during the step-freezing process, respectively. The results are given for the following instances (in hours): 0, 1, 2, 4, 16, and 28.3 (at the formation of final ice lens).
Despite increasing ice content, porosity hardly grows due to a high cooling rate in the first 1hr. With delaying freezing front penetration due to the latent heat release, porosity starts to increase in response to in-situ heave. Once the cooling rate becomes smaller than the critical rate of cooling, porosity increases rapidly. Porosity increase curves up in a bulb shape after 16hr. The increasing trend in porosity distribution indicates the rhythmic ice lens banding. The coordinate system is updated each step taking into account the heave by deforming mesh.
The stress acting on the segregation heave zone is simply equal to the applied pressure in this one-dimensional analysis. The effect of pore water pressure at the freezing front is negligible for the in-situ Fairbanks silt because of its high hydraulic conductivity of the unfrozen soil. Therefore, calculating the cooling rate is most important to determine the SP value for in-situ Fairbanks silt according to eq. ].
The detailed analyses of segregation heave are shown in Figure 4.7. The variation of the simulated cooling rate is compared with the observation. The simulated values agreed well with the observed rates as shown in Figure 4.7b. In response to the cooling rate, the SP values are in agreement with the observed measurements during the growth of ice lenses (see Figure 4.7c).
The SP value starts to increase at approximately 10hr in response to reaching the critical rate of cooling. The simulated SP value gradually increases as the cooling rate decreases. The maximum SP value is obtained at approximately 20hr. After the peak, the SP value starts to decrease. The SP concept is not available to predict the final ice lens growth as discussed in Chapter 3. After the formation of final ice lens, the simulated segregation heave is overestimated (see Figure 4.7a).
The segregation heave was successfully simulated at 50kPa applied pressure considering the effect of cooling rate during transient condition. The rest of step-freezing tests at different applied pressure are simulated as well. The variation of simulated total heaves appears to fit the observed ones in each step-freezing test (see Figure 4.8a). The amounts of simulated total heaves at the formation of ice lens agreed well with those of observed ones as well (see Figure 4.8b). The simulated total heave decreased with increasing applied pressure. The proposed fitting method was also valid for the effect of the applied pressure.
The significance of this simulation is that a trial-and-error approach was not needed to adjust input parameters in order to obtain a good match with the observed measurements. The verification results of a series of step-freezing tests confirm the validity of the developed frost heave model and the SP values obtained by the proposed fitting method.
Summary
A frost heave model was developed applying the SP concept. Significant findings from this study are:
The volume expansion is considered in the active system, which simulates the growing ice lens and frozen fringe. Furthermore, the active system has two volumetric expansion zones: in-situ freezing zone and segregation freezing zone, respectively. The proposed model could reproduce the cooling rate of segregation freezing zone.
SP porosity growth function was created considering the effects of cooling rate and applied pressure. The proposed fitting method was applied to obtain the input SP values. The simulated SP values agreed well with the observed values in a series of step-freezing tests during transient freezing condition.
The developed frost heave model applying the SP porosity growth function simulated a series of step-freezing tests with remarkably accurate results. The model was implemented applying deforming mesh in a finite element code, ABAQUS. The developed frost heave model will be used for multi-dimensional frost heave simulation in following chapters.
Figures:
Figure 4.1 Active system of the proposed frost heave model.
Figure 4.2 Unfrozen water characteristic curve for in-situ Fairbanks silt.
Figure 4.3 Three-phase relations.
Figure 4.4 Flow chart of the user-defined subroutine.
Figure 4.5 Boundary conditions for simulated processes.
Figure 4.6 Freezing process for a step-freezing test at 50kPa applied pressure: (a) distribution of temperature; (b) ice content; and (c) porosity.
Figure 4.7 Freezing process for a step-freezing test at 50kPa applied pressure: (a) variation of segregation heave; (b) cooling rate; and (c) SP.
Figure 4.8 Comparison of total heave between simulated and observed for a series of step-freezing tests at different applied pressure: (a) variation; and (b) total heave amount at the formation of the final ice lens.
Tables:
Table 4.1 Boundary and initial conditions of simulations for a series of step-freezing tests.
Table 4.2 Thermal properties (Johnston 1981).
Table 4.3 Input SP parameters for simulated step-freezing tests.
