Phase Change

This example demonstrates how to model a phase change and predicts its impact on a heat transfer analysis. When a material changes phase, for instance from solid to liquid, energy is added to the solid. Instead of creating a temperature rise, the energy alters the material’s molecular structure. Heat consumed or released by a phase change affects fluid flow, magma movement and production, chemical reactions, mineral stability, and many other earth-science applications.

Figure 1: Material properties as functions of temperature.

This 1D example uses the Phase Change Material subfeature from the Heat Transfer Module to examine transient temperature transfer in a rod of ice that heats up and changes to water. In particular, the model demonstrates how to handle material properties that vary as a function of temperature.

This model proceeds as follows. First, estimate the ice-to-water phase change using the transient conduction equation with the latent heat of fusion. Next, compare the first solution to estimates that neglect latent heat. Finally, run additional simulations to evaluate impacts of the temperature interval over which the phase change occurs.

Model Definition

This example describes the ice-to-water phase change along a 1-cm rod of ice. At its left end the rod is insulated, and at the other temperature is maintained at 80 °C. Values for thermal properties depend on the phase. They are presented in Table 1, at 265 K for ice and 300 K for water.

Table 1: Material properties of ice and water
Material property	Ice (at 265 K)	Water (at 300 K)
Density	918 kg/m3	997 kg/m3
Heat capacity at constant pressure	2052 J/(kg·K)	4179 J(kg·K)
Thermal conductivity	2.31 W/(m·K)	0.613 W/(m·K)

The latent heat of fusion, lm, is 333.5 kJ/kg and the rod is initially at –20 °C.

During the ice-to-water phase change, the density is modified, resulting in a volume compression. The material coordinates express all transformations in the initial coordinate system, when ice occupies all the domain. Assuming that there is no mixing in the liquid phase, the conduction equation in material coordinates can be used. It simplifies the model since you do not need to calculate the velocity field resulting from density variations during phase change. The conduction equation in material coordinates reads

(1)

where ρ (kg/m3) is the density, Ceq (J/kg·K) is the effective heat capacity at constant pressure, keq is the effective thermal conductivity (W/m·K), T is temperature (K), and Q is a heat source (W/m3).

The material properties ρ and keq of water must be in material coordinates. Because values given in Table 1 come from measurements, they correspond to spatial coordinates. Hence, conversion into material coordinates is necessary. In 1D models, you just have to multiply by the ratio of densities, ρratio:

With this transformation, the density of water, ρ, in material coordinates is ρ = ρwaterρratio = ρice. This is consistent with conservation of mass because the integral of ρ over the geometry domain remains constant in time.

The boundary conditions for this model are

•

thermal insulation at x = 0;

•

fixed temperature at x = 0.01; the fixed temperature creates a temperature discontinuity at the start time. You can thus replace Thot by a smoothed step function Tright that increases the temperature from T0 to Thot in 0.1 s.

Results and Discussion

Figure 2 shows images of the temperature distribution in time, predicted with latent heat. The system is solid ice at t = 0, and water content increases with time.

Figure 2: Temperature estimates with latent heat at t = 0 s, 15 s, 30 s, 45 s, 60 s, 2 min, 3 min, 4 min, …, 20 min.

The distributions all level out around the 0 °C temperature point because not all of the energy is going toward a temperature rise; some is being absorbed to change the molecular structure and change the phase.

The solution in Figure 3 shows temperature estimates for the simulation without latent heat.

Figure 3: Temperature estimates without latent heat at t = 0 s, 15 s, 30 s, 45 s, 60 s, 2 min, 3 min, 4 min, …, 20 min.

A change of profile also occurs at 0 °C but is less visible. Because latent heat is not accounted for, this change is here due to the different thermophysical properties of water below and above 0 °C.

Figure 4 shows results for different solid-to-liquid intervals at three times. The smaller the interval, the sharper the bend in the temperature profile at zero temperature, T. In the simulations, narrowing the temperature interval to a step change, for example, comes at a large computational cost. In the figure, the results for the wide and narrow pulses compare closely.

Figure 4: Temperature estimates for different temperature intervals for latent heat consumption. Estimates are for dT intervals of 0.5 K (blue), 1 K (green), and 2 K (red) at t = 30 s (three curves at bottom), 5 min (three curves at middle), and 10 min (three curves on top).

References

1. S.E. Ingebritsen and W.E. Sanford, Groundwater in Geologic Processes, Cambridge University Press, 1998.

2. N.H. Sleep and K. Fujita, Principles of Geophysics, Blackwell Science, 1997.

3. D.L. Turcotte and G. Schubert, Geodynamics: Applications of Continuum Physics to Geological Problems, 2nd ed., Cambridge University Press, 2002.

Subsurface_Flow_Module/Heat_Transfer/phase_change

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

The interface with its feature together with the subfeature solves for the temperature and automatically calculates the equivalent conductivity and the equivalent specific heat capacity.

Click .

Click

In the tree, select .

Click

Geometry 1

Interval 1 (i1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Coordinates (m)
0
0.01

Click

Form Union (fin)

In the window, click .

In the window for , click

Click the

button in the toolbar.

Global Definitions

The following steps describe how the model parameters are defined.

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
T_trans	0[degC]	273.15 K	Transition temperature
dT	1[K]	1 K	Transition interval
lm	333.5[kJ/kg]	3.335E5 J/kg	Latent heat of fusion
T_0	-20[degC]	253.15 K	Initial temperature of the rod
T_hot	80[degC]	353.15 K	Temperature of hot water
rho_ice	918[kg/m^3]	918 kg/m³	Density of ice
rho_water	997[kg/m^3]	997 kg/m³	Density of water
rho_ratio	rho_ice/rho_water	0.92076	Ratio of densities

Step 1 (step1)

In the toolbar, click

and choose .

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

Locate the section. In the text field, type 0.05.

In the text field, type T_0.

In the text field, type T_hot.

Click to expand the section. Click

Materials

Ice

In the toolbar, click

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

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	2.31	W/(m·K)	Basic
Density	rho	rho_ice	kg/m³	Basic
Heat capacity at constant pressure	Cp	2052	J/(kg·K)	Basic
Ratio of specific heats	gamma	1	1	Basic

Water

In the toolbar, click

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

Select Domain 1 only.

Because the model is solved in material coordinates, the water density and thermal conductivity are converted.

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	0.613[W/(mK)]rho_ratio	W/(m·K)	Basic
Density	rho	rho_water*rho_ratio	kg/m³	Basic
Heat capacity at constant pressure	Cp	4179	J/(kg·K)	Basic
Ratio of specific heats	gamma	1	1	Basic

Heat Transfer in Fluids (ht)

Fluid 1

In the window, under click .

Phase Change Material 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the ΔT1 → 2 text field, type dT.

In the L1 → 2 text field, type lm.

Locate the section. From the list, choose .

Initial Values 1

In the window, click .

In the window for , locate the section.

In the T text field, type T_0.

Temperature 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

In the T0 text field, type T_right(t[1/s]).

Mesh 1

Follow the steps below to generate a relatively fine mesh of 120 elements.

Edge 1

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

In the text field, type 120.

Click

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

Click

In the dialog box, type 15 in the text field.

In the text field, type 60.

Click .

In the window for , locate the section.

Click

In the dialog box, type 120 in the text field.

In the text field, type 60.

In the text field, type 1200.

Click .

For more robust convergence, tighten the relative tolerance, which controls the size of the time steps taken by the solver.

In the window for , locate the section.

From the list, choose .

In the text field, type 1e-3.

In the toolbar, click

Results

Temperature (ht)

All the parameter values in this model have a time unit of seconds, so the output time you enter here gives a total simulation time of 20 minutes. Different output intervals can be generated by adding other range commands as it is done above. Within the first minute, solution data is stored every 15 seconds, whereas for the remaining simulation period, the data is only stored every 60 seconds.

A line plot of the temperature distribution along the rod for all times is automatically produced. To generate Figure 2, you only need to change the temperature unit.

Line Graph

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Phase Change Without Latent Heat

To analyze the impact of the latent heat terms on the phase change model, a parametric sweep with a single value is used to set the latent heat to 0 in a new study.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Study 2

Step 1: Time Dependent

In the window for , locate the section.

Click

In the dialog box, type 60 in the text field.

In the text field, type 15.

Click .

In the window for , locate the section.

Click

In the dialog box, type 120 in the text field.

In the text field, type 1200.

In the text field, type 60.

Click .

In the window for , locate the section.

From the list, choose .

In the text field, type 1e-3.

In the window, click .

Click to expand the section. Select the check box.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
lm (Latent heat of fusion)	0	J/kg

In the toolbar, click

Results

Temperature, No Latent Heat

In the window for , type Temperature, No Latent Heat in the text field.

To generate Figure 3, you only need to change the units in the automatically generated temperature plot.

Line Graph

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

To be able to keep track of the different studies, rename the datasets containing the solutions of and .

Solution 1, lm Included

In the window, expand the node, then click .

In the window for , type Solution 1, lm Included in the text field.

Solution 2, lm Excluded

In the window, under click .

In the window for , type Solution 2, lm Excluded in the text field.

Phase Change for Varying Transition Intervals

Solutions to the phase change problem vary with the range in temperatures dT over which you assume that the phase transition occurs. To visualize the impact of different transition widths sample results from the original simulation and compare those estimates to results from simulations with varying dT values.

Add Study

In the toolbar, click