PDF
Artificial Ground Freezing for
Tunnel Excavation
Introduction
Artificial ground freezing is used in tunnel excavation when the soil is water-bearing and conventional stabilization methods are insufficient. Because it is an expensive technique, it is typically applied only to short, critical sections where it provides reliable ground support and ensures watertight conditions.
Numerical modeling helps in predicting the freezing process, optimizing the design, and reducing both risks and costs. For the most accurate predictions, a full Thermo-Hydraulic-Mechanical (THM) model should be used. This approach captures the coupled thermal, hydraulic, and mechanical processes — such as heat transfer, phase change, water migration, and changes in soil properties.
Model Definition
THM modeling of the freezing process captures the following processes:
•
Heat conduction in the frozen soil
•
Heat convection due to groundwater flow
•
Phase change during the freezing process
•
Water migration (cryosuction)
•
Poroelasticity and thermal expansion
Heat Conduction and Convection
The time-dependent heat transfer equation for porous media is
The effective volumetric heat capacity (ρCp)eff and the effective thermal conductivity keff of the porous medium is a combination of the fluid (f) and solid-matrix (s) properties:
where εp represents the porosity.
The fluid is transported by the Darcy velocity u and undergoes phase change between water and ice, which alters its properties. This effect is accounted for by using the apparent heat capacity method, described below.
Phase Change
The phase change is a key effect that leads to a strongly coupled THM system. The apparent heat capacity method describes the phase transition within a temperature interval, incorporating the latent heat as an additional term in the heat capacity.
The volume fractions of ice, θ1, and water, θ2, are computed using a phase transition function α12 = θ2. This phase transition function corresponds to the soil freezing characteristic curve (SFC), which describes the relationship between temperature and the unfrozen water content in the soil.
A common relation is the van Genuchten-type relationship (Ref. 1):
(1)
The van Genuchten parameters α, n, and m = 1 − 1/n are soil dependent and need to be determined by curve fitting to measurement data. The cryogenic suction head ψ(T) describes the relation between ice pressure, water pressure, and temperature
Water migration
The cryogenic suction head can be derived from the generalized Clapeyron equation
(2)
with the latent heat λ and the freezing temperature Tf. This describes the strong suction of water toward the freezing front.
Poroelasticity and Thermal Expansion
The flow velocity is determined by Darcy’s law:
In this example, the permeability κ depends on the water saturation as described by the van Genuchten relationship:
where κ0 is the intrinsic permeability, and n is the same van Genuchten parameter used for the phase change relationship.
The fluid pressure in the pores acts as a load on the soil. Biot’s poroelasticity model reads
where S is the stress tensor, αB the Biot–Willis coefficient, pA and pref are the absolute and reference pressure, and fV stands for volumetric loads like gravity. Deformation also leads to changes in porosity, which are taken into account by the poroelasticity coupling.
In addition to the pore pressure, ice formation leads to thermal strains due to the volumetric expansion of ice which is about 9%:
.
Model Setup
The construction site spans an area measuring 20 m in width, 40 m in length, and approximately 20 m in depth, with a variable surface topography. Eight freezing pipes, each with a diameter of 230 mm, are arranged in a semicircular pattern with a slight incline following the ground profile. The tunnel has a diameter of 5.6 m, and the freezing pipes are positioned 3.8 m from the tunnel center. This configuration allows for the formation of a sufficiently thick ice wall to support the excavation. The underlying material consists of an impermeable rock, where artificial ground freezing is not required. All other model parameters are listed in Table 1.
Table 1: Model Parameter
Name
Expression
Description
por
0.4
Porosity
kappa0
1e-10[m^2]
Permeability
lambda
334[kJ/kg]
Latent heat
T0
8[degC]
Ground temperature
Tc
-30[degC]
Cooling temperature
eta
0.05
Suction potential
n
0.46
Van-Genuchten parameter
The material properties for ice, water and the solid matrix are listed in Table 2.
Table 2: Material Properties
Material Property
Water
Ice
Porous Matrix
Density (kg/m3)
1000
920
2800
Dynamic Viscosity (Pa·s)
1e-3
-
-
Heat capacity at constant pressure (J/(kg·K)
4182
2060
800
Thermal conductivity (W/(m·K))
0.6
2.14
2.4
Young’s Modulus (MPa)
-
-
50
Poisson’s ratio (1)
-
-
0.3
Biot-Willis coefficient (1)
-
-
0.95
The Heat Transfer in Porous Media interface is used to describe the heat transport by conduction and convection, as well as accounting for phase change using Equation 1. The surface temperature and the groundwater temperature are set to T0. Heat flux boundary conditions are applied at the top and pipes surfaces, where the inner temperature is set to Tc.
For Darcy’s law, constant head boundary conditions are applied to generate a slow groundwater flow. To account for the amount of frozen water, a mass source is added according to
The Solid Mechanics interface assumes a linear elastic material under gravity loading. The bottom surface is fixed, while the vertical surfaces use roller conditions. The pipes are defined using a Rigid Connector boundary condition, which allows the pipes to move while keeping their shape.
First, a stationary study is computed without cooling to reach the initial state under gravity, an initial Darcy velocity contribution, and an initial (constant) temperature distribution. A subsequent time dependent study is used to model the freezing process over 30 days.
Results and Discussion
Figure 1 shows the deformation after 30 days, together with the flow field and the isosurface for ice with θ1 = 0.5 around the pipes.
Figure 1: Cut plane showing the deformation in cm and the flow field using streamlines. The isosurface indicates the ice front.
The plot for the ice volume fraction at different times is pictured in Figure 2.
Figure 2: Ice volume fraction at different times.
The ice thickness in the upstream direction and near the bottom surface is the most critical for stability. Therefore, the ice thickness is computed at different locations based on different criteria.
One approach is to use ice saturation as the criterion. In this case, a solid ice wall is assumed to form when the ice saturation exceeds a certain threshold (here, Si > 0.9).
Alternatively, a temperature-based criterion can be used, where the frozen zone is defined as the region where the temperature drops below a certain value (here, T<-5 °C).
The choice of the criterion and threshold values depends on the specific soil properties. The results are shown in Figure 3.
Figure 3: Ice thickness at different locations based on the saturation (solid) or temperature isotherm (dotted).
Reference
1. L.Sun, X. Tang, K. Ramzy Aboayanah, Q. Zhao, Q. Liu, and G. Grasselli “A coupled cryogenic thermo-hydro-mechanical model for frozen medium: Theory and implementation in FDEM,” J. Rock Mech. Geotech. Eng., vol. 16, no. 11,pp. 4335–4353, 2024; doi.org/10.1016/j.jrmge.2023.09.007.
Application Library path: Subsurface_Flow_Module/Poromechanics/artificial_ground_freezing
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
Focus on heat transfer first, as it is the key component in THM modeling and has impact on all other coupled processes. Once the Heat Transfer in Porous Media interface is properly set up, add the Poroelasticity multiphysics interface.
2
In the Select Physics tree, select Heat Transfer > Porous Media > Heat Transfer in Porous Media (ht).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies > Stationary.
6
Click  Done.
Geometry 1
1
In the Geometry toolbar, click Insert Sequence and choose Insert Sequence.
2
Browse to the model’s Application Libraries folder and double-click the file artificial_ground_freezing_geom_sequence.mph.
3
In the Geometry toolbar, click  Build All.
4
Click the  Wireframe Rendering button in the Graphics toolbar.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file artificial_ground_freezing_parameters.txt.
Next, define the functions for the cryogenic suction head (Equation 2) and the soil freezing characteristic curve (Equation 1). These are essential for capturing the phase change.
Definitions
Cryogenic Suction Head
1
In the Definitions toolbar, click  Analytic.
2
In the Settings window for Analytic, type Cryogenic Suction Head in the Label text field.
3
In the Function name text field, type hcs.
4
Locate the Definition section. In the Expression text field, type lambda*log(T/273.15[K])/g_const.
5
In the Arguments text field, type T.
6
Locate the Units section. In the Function text field, type m.
7
In the table, enter the following settings:
Argument
Unit
T
K
8
Locate the Plot Parameters section. In the table, enter the following settings:
Plot
Argument
Lower limit
Upper limit
Fixed value
Unit
√
T
270
275
0
K
9
Click  Plot.
Soil Freezing Characteristic Curve
1
In the Definitions toolbar, click  Analytic.
2
In the Settings window for Analytic, type Soil Freezing Characteristic Curve in the Label text field.
3
In the Function name text field, type sfc.
4
Locate the Definition section. In the Expression text field, type if(T<273.15[K], (1+abs(eta*hcs(T))^(1/(1-n)))^(-n), 1).
5
In the Arguments text field, type T.
6
Locate the Units section. In the table, enter the following settings:
Argument
Unit
T
K
7
Click to expand the Local Parameters section. In the table, enter the following settings:
Name
Expression
eta
0.05
n
0.46
8
Locate the Plot Parameters section. In the table, enter the following settings:
Plot
Argument
Lower limit
Upper limit
Fixed value
Unit
√
T
265
275
0
K
9
Click  Plot  and compare with the plot below.
Note the steep gradient near the freezing point. This can be numerically challenging as similar steep gradients will appear in multiple parts of the model.
Define variables for the water and ice saturation.
Variables 1
1
In the Model Builder window, right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Variables section.
3
In the table, enter the following settings:
Name
Expression
Unit
Description
Sw
nojac(sfc(T))
 
Water volume fraction
The nojac operator helps improve convergence by excluding certain terms from the Jacobian matrix. When a function has very steep gradients, it can make the solver unstable. By using nojac, the steep part of the function is temporarily ignored in this prediction, making it easier for the solver to handle.
4
In the table, enter the following settings:
Name
Expression
Unit
Description
Si
1-Sw
 
Ice volume fraction
Qm
por*(mat2.def.rho-mat1.def.rho)*d(sfc(T),TIME)
 
Mass source for freezing
To account for the mass transfer between liquid water and ice during phase change, the time derivative of the freezing curve must be explicitly evaluated.
Materials
Continue with adding materials for water, ice, and the porous material. The material properties are populated later, after the physics is set up. Then COMSOL automatically detects which material properties are required.
Water
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, type Water in the Label text field.
Ice
1
Right-click Materials and choose Blank Material.
2
In the Settings window for Material, type Ice in the Label text field.
Porous Material 1 (pmat1)
1
Right-click Materials and choose More Materials > Porous Material.
2
Select Domain 1 only.
Heat Transfer in Porous Media (ht)
Fluid 1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) > Porous Medium 1 click Fluid 1.
Phase Change Material 1
1
In the Physics toolbar, click  Attributes and choose Phase Change Material.
2
In the Settings window for Phase Change Material, locate the Phase Change section.
3
From the Phase transition function list, choose User defined. In the α1 → 2 text field, type Sw.
4
In the L1 → 2 text field, type lambda.
5
Locate the Phase 1 section. From the Material, phase 1 list, choose Ice (mat2).
6
Locate the Phase 2 section. From the Material, phase 2 list, choose Water (mat1).
Porous Matrix 1
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) > Porous Medium 1 click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the Define list, choose Solid phase properties.
The initial temperature is 8°C. This is also the ambient temperature and groundwater temperature.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the T text field, type T0.
Temperature 1
1
In the Physics toolbar, click  Boundaries and choose Temperature.
2
Select Boundary 1 only.
3
In the Settings window for Temperature, locate the Temperature section.
4
In the T0 text field, type T0.
Heat Flux: Surface
1
In the Physics toolbar, click  Boundaries and choose Heat Flux.
2
In the Settings window for Heat Flux, type Heat Flux: Surface in the Label text field.
3
Locate the Heat Flux section. From the Flux type list, choose Convective heat flux.
4
In the h text field, type 5.
5
In the Text text field, type T0.
6
Select Boundary 4 only.
For the cooling pipes, apply a heat flux boundary condition with a high heat transfer coefficient and a cooling temperature of -30°C. To ensure a smooth transition from the initial temperature to the cooling temperature, which improves solver stability, use a step function to decrease the pipe temperature from T0 to Tc.
Definitions
Step 1 (step1)
1
In the Definitions toolbar, click  More Functions and choose Step.
2
In the Settings window for Step, locate the Parameters section.
3
In the Location text field, type 4.
4
In the From text field, type T0.
5
In the To text field, type Tc.
6
Click to expand the Smoothing section. In the Size of transition zone text field, type 6.
Heat Transfer in Porous Media (ht)
Heat Flux: Pipes
1
In the Physics toolbar, click  Boundaries and choose Heat Flux.
2
In the Settings window for Heat Flux, type Heat Flux: Pipes in the Label text field.
3
Locate the Boundary Selection section. From the Selection list, choose Pipes.
4
Locate the Heat Flux section. From the Flux type list, choose Convective heat flux.
5
In the h text field, type 180.
6
In the Text text field, type step1(t[1/h]).
The heat transfer setup has been completed. Proceed with the poroelasticity multiphysics coupling.
Add Physics
1
In the Physics toolbar, click  Add Physics to open the Add Physics window.
2
Go to the Add Physics window.
3
In the tree, select Structural Mechanics > Poroelasticity > Poroelasticity, Solid.
4
Click the Add to Component 1 button in the window toolbar.
5
In the Physics toolbar, click  Add Physics to close the Add Physics window.
Solid Mechanics (solid)
Gravity 1
In the Physics toolbar, click  Global and choose Gravity.
Fixed Constraint 1
1
In the Physics toolbar, click  Boundaries and choose Fixed Constraint.
2
Select Boundary 3 only.
Roller 1
1
In the Physics toolbar, click  Boundaries and choose Roller.
2
Select Boundaries 1, 2, 5, and 38 only.
Rigid Connector 1
1
In the Physics toolbar, click  Boundaries and choose Rigid Connector.
2
In the Settings window for Rigid Connector, locate the Boundary Selection section.
3
From the Selection list, choose Pipes.
Rigid means the pipes can move as a whole, but its shape stays fixed.
Darcy’s Law (dl)
1
In the Model Builder window, under Component 1 (comp1) click Darcy’s Law (dl).
2
In the Settings window for Darcy’s Law, locate the Gravity Effects section.
3
Select the Include gravity checkbox.
4
Click to expand the Discretization section. From the Pressure list, choose Linear.
Darcy’s law only applies to the liquid water in the pores. If a pore is filled with ice, no flow can occur there. This is ensured with the following step.
Porous Matrix 1
1
In the Model Builder window, under Component 1 (comp1) > Darcy’s Law (dl) > Porous Medium 1 click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the εp list, choose User defined. In the associated text field, type por*Sw.
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Darcy’s Law (dl) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
Click the Hydraulic head button.
Add the mass source term to account for liquid water turning into ice.
Mass Source 1
1
In the Physics toolbar, click  Domains and choose Mass Source.
2
Select Domain 1 only.
3
In the Settings window for Mass Source, locate the Mass Source section.
4
In the Qm text field, type Qm.
Hydraulic Head 1
1
In the Physics toolbar, click  Boundaries and choose Hydraulic Head.
2
Select Boundary 1 only.
3
In the Settings window for Hydraulic Head, locate the Hydraulic Head section.
4
In the H0 text field, type 10[cm].
Hydraulic Head 2
1
In the Physics toolbar, click  Boundaries and choose Hydraulic Head.
2
Select Boundary 38 only.
Add Multiphysics
Now set up the remaining couplings by linking the Darcy velocity to the heat transfer interface to capture convection, and by adding the thermal expansion of ice.
1
In the Physics toolbar, click  Add Multiphysics to open the Add Multiphysics window.
2
Go to the Add Multiphysics window.
3
In the tree, select No Predefined Multiphysics Available for the Selected Physics Interfaces.
4
Find the Select the physics interfaces you want to couple subsection. In the table, clear the Couple checkbox for Darcy’s Law (dl).
5
In the tree, select Structural Mechanics > Thermal–Structure Interaction > Thermal Stress, Solid.
6
Click the Add to Component button in the window toolbar.
7
In the Physics toolbar, click  Add Multiphysics to close the Add Multiphysics window.
Multiphysics
Thermal Expansion 1 (te1)
1
In the Settings window for Thermal Expansion, locate the Thermal Expansion Properties section.
2
From the Input type list, choose Thermal strain.
3
From the dL list, choose User defined. In the associated text field, type 0.09*Si*por.
Heat Transfer in Porous Media (ht)
Fluid 1
1
In the Model Builder window, under Component 1 (comp1) > Heat Transfer in Porous Media (ht) > Porous Medium 1 click Fluid 1.
2
In the Settings window for Fluid, locate the Heat Convection section.
3
From the u list, choose Total Darcy velocity field (dl/porous1).
Define the permeability using a van Genuchten like relationship for the dependence of the permeability on the water saturation.
Definitions
Permeability
1
In the Definitions toolbar, click  Analytic.
2
In the Settings window for Analytic, type Permeability in the Label text field.
3
In the Function name text field, type kappa.
4
Locate the Definition section. In the Expression text field, type kappa0*Sw*(1-(1-Sw^(1/n))^n)^2.
5
In the Arguments text field, type Sw.
6
Locate the Local Parameters section. In the table, enter the following settings:
Name
Expression
n
0.46
7
Click  Plot.
The physics setup is complete, and the material properties can now be specified.
Materials
Porous Material 1 (pmat1)
1
In the Model Builder window, expand the Component 1 (comp1) > Materials > Porous Material 1 (pmat1) node, then click Porous Material 1 (pmat1).
2
In the Settings window for Porous Material, locate the Phase-Specific Properties section.
3
Click  Add Required Phase Nodes.
Fluid 1 (pmat1.fluid1)
1
In the Model Builder window, click Fluid 1 (pmat1.fluid1).
2
In the Settings window for Fluid, locate the Fluid Properties section.
3
From the Material list, choose Water (mat1).
Solid 1 (pmat1.solid1)
1
In the Model Builder window, click Solid 1 (pmat1.solid1).
2
In the Settings window for Solid, locate the Solid Properties section.
3
In the θs text field, type 1-por.
4
Locate the Material Contents section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
2800
kg/m³
Basic
Thermal conductivity
k_iso ; kii = k_iso, kij = 0
2.4
W/(m·K)
Basic
Heat capacity at constant pressure
Cp
800
J/(kg·K)
Basic
Porous Material 1 (pmat1)
1
In the Model Builder window, click Porous Material 1 (pmat1).
2
In the Settings window for Porous Material, locate the Homogenized Properties section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Density
rho
2800
kg/m³
Basic
Young’s modulus
E
50[MPa]
Pa
Young’s modulus and Poisson’s ratio
Poisson’s ratio
nu
0.3
1
Young’s modulus and Poisson’s ratio
Biot-Willis coefficient
alphaB
0.95
1
Poroelastic material
Permeability
kappa_iso ; kappaii = kappa_iso, kappaij = 0
kappa(Sw)
m²
Basic
Water (mat1)
1
In the Model Builder window, click Water (mat1).
2
In the Settings window for Material, locate the Material Contents section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Thermal conductivity
k_iso ; kii = k_iso, kij = 0
0.6
W/(m·K)
Basic
Density
rho
1000
kg/m³
Basic
Heat capacity at constant pressure
Cp
4182
J/(kg·K)
Basic
Dynamic viscosity
mu
1e-3
Pa·s
Basic
Ice (mat2)
1
In the Model Builder window, click Ice (mat2).
2
In the Settings window for Material, locate the Material Contents section.
3
In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Thermal conductivity
k_iso ; kii = k_iso, kij = 0
2.14
W/(m·K)
Basic
Density
rho
920
kg/m³
Basic
Heat capacity at constant pressure
Cp
2060
J/(kg·K)
Basic
Mesh 1
Free Triangular 1
1
In the Mesh toolbar, click  More Generators and choose Free Triangular.
2
Select Boundary 2 only.
Distribution 1
1
Right-click Free Triangular 1 and choose Distribution.
2
Click the  Select Box button in the Graphics toolbar. Draw a box around all pipe edges on the surface that is about to be meshed. This corresponds to the following.
3
Select Edges 10, 11, 15, 17, 22, 23, 27, 30, 34, 35, 39, 42, 46, 47, 51, 54, 58, 59, 63, 66, 70, 71, 75, 78, 82, 83, 87, 90, 94, 95, 99, and 101 only.
4
In the Settings window for Distribution, click  Build Selected.
Size
1
In the Model Builder window, under Component 1 (comp1) > Mesh 1 click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Extra fine.
4
Click the Custom button.
5
Locate the Element Size Parameters section. In the Maximum element size text field, type 2.
6
In the Maximum element growth rate text field, type 1.3.
Swept 1
In the Mesh toolbar, click  Swept.
Size 1
1
Right-click Swept 1 and choose Size.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
4
Locate the Element Size Parameters section.
5
Select the Maximum element size checkbox. In the associated text field, type 2.5.
6
Click  Build All.
Study 1
Use a stationary study step to compute the system with the pipe cooling turned off, so that gravity loads, groundwater flow, and heat balance are established. The results from this stationary study then serve as the initial conditions for the subsequent time-dependent simulation of the freezing process.
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, locate the Study Settings section.
3
Clear the Generate default plots checkbox.
Step 2: Time Dependent
1
In the Study toolbar, click  Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
From the Time unit list, choose d.
4
In the Output times text field, type range(0,1,30).
Step 1: Stationary
1
In the Model Builder window, click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
Select the Modify model configuration for study step checkbox.
4
In the tree, select Component 1 (comp1) > Heat Transfer in Porous Media (ht) > Heat Flux: Pipes.
5
Click  Disable.
6
In the Study toolbar, click  Get Initial Value.
Before computing the study, which takes about 1.5 hours to solve, prepare the solver and create a plot to monitor the freezing process during computation.
Results
In the Model Builder window, expand the Results node.
Cut Plane 1
1
In the Model Builder window, expand the Results > Datasets node.
2
Right-click Results > Datasets and choose Cut Plane.
3
In the Settings window for Cut Plane, locate the Plane Data section.
4
From the Plane list, choose XZ-planes.
5
In the Y-coordinate text field, type 20.
Ice Saturation and Stress
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Ice Saturation and Stress in the Label text field.
3
Click to expand the Plot Array section. From the Array type list, choose Square.
Surface 1
1
Right-click Ice Saturation and Stress and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type Si.
4
Locate the Coloring and Style section. From the Color table list, choose Pelagic.
Streamline 1
1
In the Model Builder window, right-click Ice Saturation and Stress and choose Streamline.
2
In the Settings window for Streamline, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Darcy’s Law > Velocity and pressure > dl.u,...,dl.w - Total Darcy velocity field (spatial frame).
3
Locate the Streamline Positioning section. From the Entry method list, choose Coordinates.
4
In the x text field, type 0.
5
In the y text field, type range(0,2,20).
6
Locate the Coloring and Style section. Find the Point style subsection. From the Type list, choose Arrow.
7
Click to expand the Plot Array section. Select the Manual indexing checkbox.
Surface 2
1
Right-click Ice Saturation and Stress and choose Surface.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Solid Mechanics > Stress > solid.misesGp - von Mises stress - N/m².
3
Locate the Coloring and Style section. From the Color table list, choose Prism.
Ice Saturation and Stress
Click the  Zoom Extents button in the Graphics toolbar.
Study 1
Step 2: Time Dependent
1
In the Model Builder window, under Study 1 click Step 2: Time Dependent.
2
In the Settings window for Time Dependent, click to expand the Results While Solving section.
3
Select the Plot checkbox.
4
From the Update at list, choose Time steps taken by solver.
5
In the table, enter the following settings:
Plot group
Plot window
Ice Saturation and Stress
Graphics
Solver Configurations
In the Model Builder window, expand the Study 1 > Solver Configurations node.
Solution 1 (sol1)
1
In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) node, then click Time-Dependent Solver 1.
2
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
3
From the Steps taken by solver list, choose Intermediate.
4
Select the Initial step checkbox. In the associated text field, type 1[h]  to capture the initial temperature drop in and around the pipes.
5
In the Study toolbar, click  Compute.
Results
Freezing Front and Deformation
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Freezing Front and Deformation in the Label text field.
Isosurface 1
1
Right-click Freezing Front and Deformation and choose Isosurface.
2
In the Settings window for Isosurface, locate the Expression section.
3
In the Expression text field, type Si.
4
Locate the Levels section. From the Entry method list, choose Levels.
5
In the Levels text field, type 0.5.
6
Locate the Coloring and Style section. From the Coloring list, choose Uniform.
7
From the Color list, choose White.
8
Clear the Color legend checkbox.
Transparency 1
Right-click Isosurface 1 and choose Transparency.
Surface 1
1
In the Model Builder window, right-click Freezing Front and Deformation and choose Surface.
2
In the Settings window for Surface, locate the Data section.
3
From the Dataset list, choose Cut Plane 1.
4
From the Solution parameters list, choose From parent.
5
Click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Solid Mechanics > Displacement > solid.disp - Displacement magnitude - m.
6
Locate the Expression section. From the Unit list, choose cm.
7
Locate the Coloring and Style section. From the Color table list, choose Prism.
Freezing Front and Deformation
In the Model Builder window, click Freezing Front and Deformation.
Streamline Surface 1
1
In the Freezing Front and Deformation toolbar, click  More Plots and choose Streamline Surface.
2
In the Settings window for Streamline Surface, locate the Data section.
3
From the Dataset list, choose Cut Plane 1.
4
From the Solution parameters list, choose From parent.
5
Click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Darcy’s Law > Velocity and pressure > dl.u,...,dl.w - Total Darcy velocity field (spatial frame).
6
Locate the Streamline Positioning section. From the Entry method list, choose Coordinates.
7
In the x text field, type 0.
8
In the y text field, type 20.
9
In the z text field, type range(0,1,18).
10
Locate the Coloring and Style section. Find the Point style subsection. From the Type list, choose Arrow.
11
From the Color list, choose Custom.
12
On Windows, click the colored bar underneath, or — if you are running the cross-platform desktop — the Color button  and select a dark gray tone.
Surface 2
1
Right-click Freezing Front and Deformation and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type z.
4
Locate the Coloring and Style section. Clear the Color legend checkbox.
5
From the Coloring list, choose Gradient.
6
From the Bottom color list, choose Custom.
7
Click Define custom colors.
8
Set the RGB values to 128, 128, and 128, respectively.
9
Click Add to custom colors.
10
Click Show color palette only or OK on the cross-platform desktop.
Selection 1
1
Right-click Surface 2 and choose Selection.
2
Select Boundaries 3, 5, and 38 only.
Surface 3
In the Model Builder window, right-click Freezing Front and Deformation and choose Surface.
Selection 1
1
In the Model Builder window, right-click Surface 3 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Pipes.
Material Appearance 1
1
In the Model Builder window, right-click Surface 3 and choose Material Appearance.
2
In the Settings window for Material Appearance, locate the Appearance section.
3
From the Appearance list, choose Custom.
4
From the Material type list, choose Steel.
Freezing Front and Deformation
1
In the Model Builder window, under Results click Freezing Front and Deformation.
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
Clear the Plot dataset edges checkbox.
4
In the Freezing Front and Deformation toolbar, click  Plot.
5
Click to expand the Title section. From the Title type list, choose None.
6
Click to expand the Quality section. From the Recover list, choose Within domains.
7
In the Freezing Front and Deformation toolbar, click  Plot.
Plot the ice thickness at different times.
Ice Saturation at Different Times
1
In the Results toolbar, click  2D Plot Group.
2
In the Settings window for 2D Plot Group, type Ice Saturation at Different Times in the Label text field.
3
Locate the Plot Array section. From the Array type list, choose Square.
4
From the Displacement list, choose Absolute.
5
In the Row displacement text field, type -25.
6
In the Column displacement text field, type 25.
Surface 1
1
Right-click Ice Saturation at Different Times and choose Surface.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type Si.
4
Locate the Coloring and Style section. From the Color table list, choose Pelagic.
Solution Array 1
1
Right-click Surface 1 and choose Solution Array.
2
In the Settings window for Solution Array, locate the Data section.
3
From the Time selection list, choose Interpolated.
4
In the Times (d) text field, type range(5,5,30).
5
Locate the Plot Array section. From the Array shape list, choose Square.
Add the time stamps in the next step.
Ice Saturation at Different Times
In the Model Builder window, under Results click Ice Saturation at Different Times.
Table Annotation 1
1
In the Ice Saturation at Different Times toolbar, click  More Plots and choose Table Annotation.
2
In the Settings window for Table Annotation, locate the Data section.
3
From the Source list, choose Local table.
4
In the table, enter the following settings:
x-coordinate
y-coordinate
Annotation
10
0
5d
35
0
10d
60
0
15d
10
-25
20d
35
-25
25d
60
-25
30d
5
Locate the Coloring and Style section. Clear the Show point checkbox.
6
From the Anchor point list, choose Upper middle.
7
In the Ice Saturation at Different Times toolbar, click  Plot.
8
Click the  Zoom Extents button in the Graphics toolbar.
Compare with Figure 2. Note that the thinnest part of the ice wall is in upstream direction close to the surface. To compute the ice thickness in this area at different times, follow the steps below.
Surface 1
1
In the Results toolbar, click  More Datasets and choose Surface.
2
Select Boundary 3 only.
Cut Line 2D 1
1
In the Results toolbar, click  Cut Line 2D.
2
In the Settings window for Cut Line 2D, locate the Data section.
3
From the Dataset list, choose Surface 1.
4
Locate the Line Data section. In row Point 1, set x to -20.
5
In row Point 2, set x to -10.
6
In row Point 1, set y to 5.
7
In row Point 2, set y to 5.
8
Select the Additional parallel lines checkbox.
9
In the Distances text field, type 10 20 30.
10
Click  PlotThe lines are at the bottom surface in upstream direction.
11
Click to expand the Advanced section.
Line Integration 1
1
In the Results toolbar, click  More Derived Values and choose Integration > Line Integration.
2
In the Settings window for Line Integration, locate the Data section.
3
From the Dataset list, choose Cut Line 2D 1.
4
From the Time selection list, choose Interpolated.
5
In the Times (d) text field, type range(10,5,30).
6
Locate the Multiple Curves section. Select the Evaluate each line separately checkbox.
7
Locate the Expressions section. Click  Clear Table.
Because of the soil freezing characteristic curve, a completely frozen state is never reached, and a small fraction of liquid water always remains. To evaluate when the frozen wall can be considered stable, different criteria can be used depending on soil type, flow and requirements. For example, an ice saturation of about 0.9 or, alternatively, the –5°C isotherm can be applied as indicators of a stable ice wall.
8
In the table, enter the following settings:
Expression
Unit
Description
1*(Si>0.9)
m
Ice thickness based on ice saturation
1*(T<-5[degC])
m
Ice thickness based on temperature isotherm)
9
Click  Evaluate.
Table 1
1
Go to the Table 1 window.
2
Click the Table Graph button in the window toolbar.
Results
Table Graph 1
1
In the Settings window for Table Graph, locate the Data section.
2
From the Plot columns list, choose Manual.
3
In the Columns list, choose Ice thickness based on ice saturation (m), Ice thickness based on ice saturation (m), Distance=10, Ice thickness based on ice saturation (m), Distance=20, and Ice thickness based on ice saturation (m), Distance=30.
4
Locate the Coloring and Style section. From the Width list, choose 2.
5
Click to expand the Legends section. Select the Show legends checkbox.
6
From the Legends list, choose Manual.
7
In the table, enter the following settings:
Legends
y=5m
y=15m
y=25m
y=35m
Table Graph 2
1
Right-click Results > 1D Plot Group 4 > Table Graph 1 and choose Duplicate.
2
In the Settings window for Table Graph, locate the Data section.
3
In the Columns list, choose Ice thickness based on temperature isotherm) (m), Ice thickness based on temperature isotherm) (m), Distance=10, Ice thickness based on temperature isotherm) (m), Distance=20, and Ice thickness based on temperature isotherm) (m), Distance=30.
4
Locate the Legends section. Clear the Show legends checkbox.
5
Locate the Coloring and Style section. Find the Line style subsection. From the Line list, choose Dotted.
6
From the Color list, choose Cycle (reset).
Ice Thickness
1
In the Model Builder window, under Results click 1D Plot Group 4.
2
In the Settings window for 1D Plot Group, type Ice Thickness in the Label text field.
3
Click to expand the Title section. From the Title type list, choose Manual.
4
In the Title text area, type Ice thickness based on saturation (solid) and temperature isotherm (dotted).
5
Locate the Legend section. From the Position list, choose Upper left.
6
In the Ice Thickness toolbar, click  Plot.