COMSOL 6.3 - Frozen Inclusion

This benchmark example simulates phase change in porous media. The model studies the melting process of an ice inclusion within a porous medium, and it demonstrates how to couple Darcy’s Law and the Heat Transfer in Porous Media interfaces including phase change. It is one of the test cases from the INTERFROST project — a comparison project for coupled thermo-hydraulic systems with respect to climate change in permafrost regions.

Model Definition

The model setup is shown in Figure 1. A flow channel 3 meters in length and 1 meter in width contains a squared ice inclusion of 33.3 cm in length, centered at the coordinates x = 1 m and y = 0.5 m.

Figure 1: Model geometry, initial temperature distribution, and boundary conditions.

The initial temperature of the ice inclusion is Tini = −5°C. The temperature in the porous channel is Tinw = +5°C. The fluid flow is driven by a hydraulic head gradient of ΔH = 3% over the length of the channel.

The initial pressure and velocity field are computed from a steady state flow simulation for a uniformly distributed temperature and permeability (as if the ice inclusion was not present). Considering the symmetry of the benchmark test model, only the lower half of the domain is modeled.

The following simplifications have been made in the INTERFROST benchmark and therefore are also assumed in this example: No thermal dispersion is included in the heat transfer equation, and water density and dynamic viscosity are considered constant with respect to the temperature (see Table 1).

Model equations

The mass conservation for water and ice in a porous medium with porosity εp is given by:

where ρi and ρw are the ice and liquid water densities, and Si and Sw are the saturations of ice and water, respectively. Assuming constant porosity and ice density, and using the relationship Si = 1-Sw the left hand side becomes:

(1)

To account for the density change of water, the compressibility is introduced:

Using Equation 1 together with this definition the mass conservation becomes:

The feature in the Darcy’s Law interface allows the definition of a user-defined storage coefficient Sp (SI unit: 1/Pa) which is written as:

(2)

to obtain the first term on the left hand side. The right hand side of the equation can be incorporated using a mass source Qm according to:

The liquid water saturation Sw depends on the phase change as:

(3)

Here, Sw, res it the residual liquid water saturation and fphtr is the phase transition function which is defined as follows:

(4)

In the INTERFROST benchmark it is assumed that the mushy ice zone extends from 0°C to −1°C. Therefore, W, which corresponds to the slope of the phase transition function, is set to 0.5 K.

The heat transfer in the porous medium is described by the following equation:

(5)

here (ρC)eq is the equivalent value of density ρ (kg/m3) and heat capacity C at constant pressure (J/kg·K), keq is the effective thermal conductivity (W/m·K), T is the temperature (K), and Q is a heat source (W/m3). The variable Cw is the effective fluid’s heat capacity at constant pressure.

Model data

The benchmark example describes the ice-to-water phase change of an ice inclusion within a porous channel. The values for thermal properties of water, ice, and solid matrix are presented in Table 1.

Table 1: Material properties of ice, water, and solid matrix.
Material property	Ice	Water	Solid Matrix
Density	920 kg/m3	1000 kg/m3	2650 kg/m3
Heat capacity at constant pressure	2060 J/(kg·K)	4182 J(kg·K)	835 J(kg·K)
Thermal conductivity	2.14 W/(m·K)	0.6 W/(m·K)	9 W/(m·K)

A latent heat of melting, L= 334 kJ/kg is used. Additional physical parameters used in the example are summarized in Table 2.

Table 2: Additional physical parameters.
Parameter	Description	Value
μw	Dynamic viscosity of water	1.793·10-3 Pa·s
ε	Porosity	0.37
β	Effective compressibility	10-8 Pa-1
Swres	Residual water saturation	0.05
Ω	Impedance factor	50
kint	Intrinsic permeability	1.3·10-10 m2

Results and Discussion

The simulation is run for 56 hours which is the time until the frozen inclusion is completely melted and the cooler temperatures have been convected out of the channel. However, it is also interesting to display the effect of the molten inclusion on the pressure and temperature distribution; Figure 2 shows the pressure distribution after 9 hours of simulated time when there is still some ice present. As gravity has been neglected, only the pressure effects of the ice inclusion are shown.]

Figure 2: Pressure distribution after 9 hours when there is still some ice present in the model domain.

The temperature field is shown at the same time (after 9 hours) in Figure 3 and the liquid water saturation in Figure 4. Both figures indicate the ice inclusion and its influence on the temperature field. The low temperatures are convected in the downstream part of the domain.

Figure 3: Surface plot showing the temperature distribution after 9 hours.

Figure 4: Surface plot showing the liquid water saturation after 9 hours. Effectively the areas with liquid water appear in blue whereas the areas of ice appear in white.

The plots below have been used as performance measures in the INTERFROST benchmark, and therefore are also included here:

The evolution of the minimum of the simulated temperature field as a function of time is shown in Figure 5.

Figure 5: The minimum temperature in the model domain as a function of time.

Figure 6 shows the time evolution of the total heat flux exiting the system at the downstream boundary.

Figure 6: Total heat flux exiting the system via the outflow boundary as a function of time.

Figure 7 displays the evolution of the total liquid water volume in the domain as a function of time.

Figure 7: Total liquid water volume in the domain as a function of time.

Notes About the COMSOL Implementation

A residual liquid phase is present in soils even for temperatures significantly below 0°C. This residual water is not subject to phase change or convection. The definition as well as the node within the interface allow to define such . According to Ref. 1, the residual saturation Sw, res is set to 0.05.

Reference

1. C. Grenier and others, “Groundwater flow and heat transport for systems undergoing freeze-thaw: Intercomparison of numerical simulators for 2D test cases,” Adv. Water Resour., vol. 114, pp. 196–218, 2018.

Subsurface_Flow_Module/Heat_Transfer/frozen_inclusion

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Right-click and choose .

In the tree, select > > .

Right-click and choose .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

The parameters defining the material properties of the solid matrix, water, and ice have to be declared. They are available in a file and you can import them as follows:

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file frozen_inclusion_parameters.txt.

Geometry 1

Now, the model domain is drawn as two rectangles, one representing the porous domain, the other the frozen inclusion within the porous domain. Create the geometry as follows and note that due to symmetry only half of the model domain is drawn.

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 3.

In the text field, type 0.5.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.333.

In the text field, type 0.333/2.

Locate the section. From the list, choose .

In the text field, type 1.

In the text field, type 0.5-0.333/4.

In the toolbar, click

Materials

Define the materials in the next step. Introduce them as empty material nodes, first. As soon as the physics has been defined, the material node menu will show you which properties are needed for the simulation. You can then just fill in the values defined in the list.

Solid Matrix

In the window, under right-click and choose .

In the window for , type Solid Matrix in the text field.

Water

Right-click and choose .

In the window for , type Water in the text field.

Ice

Right-click and choose .

In the window for , type Ice in the text field.

Porous Material 1 (pmat1)

Right-click and choose > .

In the window for , locate the section.

From the list, choose .

Locate the section. Click

Fluid 1 (pmat1.fluid1)

In the window, click .

In the window for , locate the section.

From the list, choose .

Solid 1 (pmat1.solid1)

In the window, click .

In the window for , locate the section.

From the list, choose .

In the θs text field, type 1-epsilon_p.

With the following steps you define the residual liquid water as immobile fluid and its volume as a product of porosity and predefined residual water saturation.

Immobile Fluid 1 (pmat1.imfluid1)

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

In the θimf text field, type epsilon_p*Sw_res.

Heat Transfer in Porous Media (ht)

Continue with setting up the physics. Since some heat transfer variables are needed to define the saturation variables as explained in the Model Definition section, start with Heat Transfer in Porous Media (ht).

In the window, under click .

In the window for , locate the section.

In the Tref text field, type T_initial_w.

Porous Matrix 1

In the window, under > > click .

In the window for , locate the section.

From the list, choose , because the properties of a pure solid material are given.

Now add the immobile fluid to the node. It will automatically refer to the settings you have made earlier.

Porous Medium 1

In the window, click .

Immobile Fluids 1

In the toolbar, click

and choose .

Fluid 1

In the window, click .

In the window for , locate the section.

From the pA list, choose .

Locate the section. From the u list, choose .

Phase Change Material 1

In the toolbar, click

and choose .

In the window for , click to expand the section.

Inflate the sketch to visualize which material you have to choose for which temperatures. In this case, material 1 should be ice, material 2 water. Enter a phase transition function that has yet to be defined.

Locate the section. From the list, choose . In the α1 → 2 text field, type f_phtr(T).

In the L1 → 2 text field, type L.

Locate the section. From the list, choose .

Now define the phase transition function.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
W	0.5[K]	0.5 K	Parameter for phase transition function

Definitions (comp1)

Analytic 1 (an1)

In the window, expand the > node.

Right-click > and choose > .

In the window for , type f_phtr in the text field.

Locate the section. In the text field, type exp(-((T-273.15)/W)^2)*(T<273.15)+1*(T>=273.15).

In the text field, type T.

Locate the section. In the table, enter the following settings:


Argument	Unit
T	K

To ensure that you have entered everything correctly, plot the function as follows.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	T	270	275	0	K

Click

Now continue setting up the physics for the interface.

Heat Transfer in Porous Media (ht)

Initial Values 1

In the window, under > click .

In the window for , locate the section.

In the T text field, type T_initial_w.

Initial Values 2

In the toolbar, click

and choose .

Select Domain 2 only.

In the window for , locate the section.

In the T text field, type T_initial_i.

Temperature 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

In the T0 text field, type T_initial_w.

Outflow 1

In the toolbar, click

and choose .

Select Boundary 9 only.

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 3, 6, and 8 only.

Definitions (comp1)

For Darcy’s Law the water saturation has to be defined as temperature-dependent variable. More precisely it depends on the phase transition function which determines the fraction of water as a function of temperature during phase change. The permeability is then defined as a function of water saturation.

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Sw	Sw_res+(1-Sw_res)*f_phtr(T)		Water saturation
kr	max(10^(-Omegaepsilon_p(1-Sw)),1e-6)		Relative permeability
kappa_w	k_int*kr	m²	Permeability

Darcy’s Law (dl)

In the window, under click .

In the window for , click to expand the section.

From the list, choose .

Change the shape function to linear to reduce the computational time and increase numerical stability.

Porous Medium 1

In the window, under > click .

In the window for , locate the section.

From the list, choose . In the Sp text field, type Sw*epsilon_p*beta.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the εp list, choose . In the associated text field, type epsilon_p.

Mass Source 1

In the toolbar, click

and choose .

The mass source accounts for the addition of water during melting as described in the Model Definition section.

In the window for , locate the section.

From the list, choose .

Locate the section. In the Qm text field, type epsilon_p*(rho_i-rho_w)*d(Sw,t).

The flow is driven by an imposed gradient of the hydraulic head. This is defined as a parameter and could be varied in a parametric sweep in future simulations (as it was done in the INTERFROST project, see Reference).

Hydraulic Head 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

In the H0 text field, type delta_H_dLx*3[m].

Hydraulic Head 2

In the toolbar, click

and choose .

Select Boundary 9 only.

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 3, 6, and 8 only.

Materials

Having defined the physics, now fill in the empty expressions in the node.

Porous Material 1 (pmat1)

In the window, under > click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Permeability	kappa_iso ; kappaii = kappa_iso, kappaij = 0	kappa_w	m²	Basic

Solid Matrix (mat1)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rho_s	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	lambda_s	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	Cs	J/(kg·K)	Basic

Water (mat2)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rho_w	kg/m³	Basic
Dynamic viscosity	mu	mu_w	Pa·s	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	lambda_w	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	Cw	J/(kg·K)	Basic

Ice (mat3)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	lambda_i	W/(m·K)	Basic
Density	rho	rho_i	kg/m³	Basic
Heat capacity at constant pressure	Cp	Ci	J/(kg·K)	Basic