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 starting 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.

Heat_Transfer_Module/Phase_Change/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.

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

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

Parameter name	Parameter value list	Parameter unit
dT (Transition interval)	0.5 1 2	K