COMSOL 6.4 - Tin Melting Front

This example demonstrates how to model phase transition by a moving boundary interface according to the Stefan problem. It is adapted from the benchmark study in Ref. 1.

A square cavity containing both solid and liquid tin is submitted to a temperature difference between left and right boundaries. Fluid and solid parts are solved in separate domains sharing a moving melting front (see Figure 1). The position of this boundary through time is calculated according to the Stefan energy balance condition.

Figure 1: Square cavity with moving phase interface.

In the molten tin domain, motion generated by natural convection is expected due to the temperature gradient. This motion, in turn, influences the front displacement.

Model Definition

The geometry presented in Figure 2 shows a square of side length 10 cm filled with pure tin. The left and right boundaries are maintained at 508 K and 503 K, respectively. Because the fusion temperature of pure tin is 505 K, both liquid and solid phases co-exist in the square.

Figure 2: Geometry and boundary conditions at initial time.

The initial temperature distribution is assumed to vary linearly in the horizontal direction as shown in Figure 3.

Figure 3: Initial temperature profile.

The melting front is the vertical line initially located at x = 6 cm where the temperature is 505 K.

The fluid motion in the liquid part on the left is modeled by the incompressible Navier–Stokes equations with Boussinesq approximation, as described in Gravity and The Boussinesq Approximation sections in the CFD Module User’s Guide, with Tref = Tf the fusion temperature of tin. The reduced pressure formulation is used to enhance the accuracy of the buoyancy forces since they are relatively small compared to the other terms in the momentum balance.

As the metal melts, the solid-liquid interface moves toward the solid side. The energy balance at this front is expressed by the Stefan condition:

(1)

where ΔH is the latent heat of fusion, equal to 60 kJ/kg, v (SI unit: m/s) is the front velocity vector, n is the normal vector at the front, and Φl and Φs (SI unit: W/m2) are the heat fluxes coming from the liquid and solid sides, respectively (see Figure 4).

Figure 4: Heat fluxes at the melting front.

Table 1 reviews the material properties of tin (Ref. 2) used in this model.

Table 1: Material properties of tin.
Parameter	Description	Value
ρ0	Density	7500 kg/m3
Cp	Heat capacity	200 J/(kg·K)
k	Thermal conductivity	60 W/(m·K)
αp	Coefficient of thermal expansion	2.67·10-4 K-1
μ	Dynamic viscosity	6.0·10-3 Pa·s
Tf	Fusion temperature	505 K
ΔH	Latent heat of fusion	60 kJ/kg

Results and Discussion

Figure 5 shows the velocity profile in the fluid domain. The convective cell due to buoyancy increases the melting speed at the upper part of the cavity.

Figure 5: Velocity profile in the fluid at the end of the simulation.

At the end of the simulation, the melting front does not move anymore because balance between left and right adjacent fluxes has been reached.

In Figure 6, the temperature profile is represented jointly by a surface plot for temperature distribution and a velocity arrow plot.

Figure 6: Temperature profile at the end of the simulation.

Notes About the COMSOL Implementation

The quantities Φl and Φs, illustrated in Figure 4, can be computed on each side of the melting front using the up and down operators. The components of Φl − Φs are then given by up(ht.tfluxx)-down(ht.tfluxx) and up(ht.tfluxy)-down(ht.tfluxy). However, this method evaluates the temperature gradient which may lead to imprecisions due to the mesh discretization. Instead, the quantity (Φl − Φs) ⋅ n, involved in Equation 1, is more precisely evaluated through the Lagrange multiplier for temperature, T_lm. This variable is available when weak constraints are enabled in the region of interest, as it is the case here with the fixed temperature constraint at the melting front. For more information about weak constraints, refer to the Weak Constraint section in the COMSOL Multiphysics Reference Manual.

To handle the melting front movement, a mesh deformation is necessary. During such a transformation, matter from solid tin is removed while the same amount of liquid tin is added to the fluid. The appropriate tool for deforming the mesh without reflecting any expansion or contraction effects to the material properties is the node under .

References

1. F. Wolff and R. Viskanta, “Solidification of a Pure Metal at a Vertical Wall in the Presence of Liquid Superheat,” Int.J. Heat and Mass Transfer, vol. 31, no. 8, pp. 1735–1744, 1988.

2. V. Alexiades, N. Hannoun, and T.Z. Mai, “Tin Melting: Effect of Grid Size and Scheme on the Numerical Solution,” Proc. 5th Mississippi State Conf. Differential Equations and Computational Simulations, pp. 55–69, 2003.

Heat_Transfer_Module/Phase_Change/tin_melting_front

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

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
k_Sn	60[W/(m*K)]	60 W/(m·K)	Thermal conductivity
Cp_Sn	200[J/(kg*K)]	200 J/(kg·K)	Specific heat capacity
alpha_Sn	2.67e-4[1/K]	2.67E-4 1/K	Coefficient of thermal expansion
mu_Sn	6e-3[Pa*s]	0.006 Pa·s	Kinematic viscosity
rho_Sn	7500[kg/m^3]	7500 kg/m³	Density
DelH	60[kJ/kg]	60000 J/kg	Latent heat of fusion
Tf	505[K]	505 K	Melting point
Th	508[K]	508 K	Hot wall temperature
Tc	503[K]	503 K	Cold wall temperature

Geometry 1

Square 1 (sq1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.1.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	0.06

Select the checkbox.

Clear the checkbox.

In the toolbar, click

Click the

button in the toolbar.

Materials

Tin (Solid)

In the toolbar, click

In the window for , type Tin (Solid) in the text field.

Select Domain 2 only.

Before defining the material properties, set the solid and fluid domains in the physics interfaces to let COMSOL Multiphysics flag what properties you need to specify.

Laminar Flow (spf)

In the window, under click .

Select Domain 1 only.

Heat Transfer in Solids and Fluids (ht)

In the window, under click .

In the window for , locate the section.

In the Tref text field, type Tf.

Fluid 1

In the window, under > click .

Select Domain 1 only.

Materials

Tin (Solid) (mat1)

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
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	k_Sn	W/(m·K)	Basic
Density	rho	rho_Sn	kg/m³	Basic
Heat capacity at constant pressure	Cp	Cp_Sn	J/(kg·K)	Basic