# Numerical simulation of melting and solidification

### From Thermal-FluidsPedia

Due to the strong nonlinearity of solid-liquid phase change phenomena and the moving boundary, the problems that can be solved via the analytical method are very limited. Exact solutions and some approximate solutions for melting and solidification have been introduced in preceding sections. The limitation of these methods is that they apply to conduction-controlled one-dimensional problems. Although some investigators have attempted to solve two-dimensional phase change problems by using analytical methods (Poots, 1962; Rathjen and Jiji, 1971; Budhia and Kreith, 1973), the cases investigated represented very simple and special geometries, such as a semi-infinite corner or a semi-infinite wedge. Generally, the analytical method will not work for two- or three-dimensional problems; this is particularly true for convection-controlled phase change processes. Most real applications of phase change occur in complex two- or three-dimensional geometries. In addition, natural convection in the liquid phase often plays a significant role in the phase-change processes. For such complex solid-liquid phase-change processes, analytical solutions will not work. Therefore, it is necessary to employ numerical methods.
A large number of numerical methods have been developed and they can be divided into two groups. The first group is called *strong numerical solution*, or *classical solution*. For this group, transformed coordinate systems or deformed grids are employed to deal with the location of the solid-liquid interface. In this methodology, complex geometric regions of solid or liquid phases are transformed into fixed, simple geometric regions through the coordinate transformation technique. At the same time, the governing equations and the boundary conditions become more complicated. The governing equations of the problem can be solved using the finite difference method. Such methods also deal successfully with multidimensional one-region or two-region problems with or without natural convection in the liquid phase. The disadvantage of these methods is that they are often difficult to program and thus require a significant amount of computer time.
The second group is called the *weak solution*, or *the fixed-grid solution*. In these methods, the governing equations are written for the entire phase change region, including liquid and solid phases. The location of the solid-liquid interface is determined after the converged temperature distribution is obtained. The three principal important methods in this group are the enthalpy method (Shamsunder and Sparrow, 1975; Crank, 1984), the equivalent heat capacity method (Bonacina et al., 1973), and the temperature-transforming model (Cao and Faghri, 1990). The temperature transforming model has the advantages of both the enthalpy method and the equivalent heat capacity method and can also account for the effect of natural convection in the liquid phase. The most significant advantage of the weak solution, as compared with the strong numerical solution, is its simplicity. In this section, various fixed-grid numerical solutions of the solid-liquid phase problem will be introduced.

## Contents |

### Enthalpy Method

In this methodology, the governing energy equation is written for the entire region of the PCM, including solid and liquid phases and the interface. The enthalpy method is introduced by analyzing a conduction-controlled, two-region melting problem in a finite slab, as shown in Fig. 1. It is assumed that the densities of the liquid and solid phase are identical (). The energy equation can be written as

The enthalpy, *h*, is a function of temperature, *T*:

The variation of enthalpy *h* with temperature *T* can be plotted as shown in Fig. 2. At the melting point of the PCM, the enthalpies of solid and liquid phases at the melting point are 0 and , respectively; there is a jump of enthalpy equal to the latent heat of the PCM. The thermal conductivity of the PCM, *k*, can be expressed as

**Failed to parse (syntax error): k(T) = \left\{ \begin{array}{l} {k_s}\quad T < {T_m} \\ {k_\ell }\,\,\,\,\,T > {T_m} \\ \end{array} \right \qquad \qquad(3)**

It can be verified that the energy equation for liquid and solid phases can be obtained simply by substituting eqs. (2) and (3) into eq. (1). The energy balance at the solid-liquid interface can be obtained by analyzing the energy balance for a control volume, which includes the solid-liquid interface.
Solving for temperature, *T*, from eq. (2) yields

**Failed to parse (syntax error): T = \left\{ \begin{array}{l} {T_m} + h/{c_{ps}} & & h \le 0 \\ {T_m} & & & 0 < h < {h_{s\ell }} \\ {T_m} + (h - {h_{s\ell }})/{c_{p\ell }} & h \ge {h_{s\ell }} \\ \end{array} \right \qquad \qquad(4)**

The initial and boundary conditions of the melting problem are

The discretization of space and time is shown in Fig. 3, in which location and time are represented by *j* and *n*. The temperature at *x* = (*j* − 1)Δ*x* and time *t* = *n*Δ*t* can be represented by the symbol . Equation (1) is discretized by using an explicit scheme which uses forward difference in time and central difference in space:

The thermal conductivities at the half-grid, and , are calculated using the harmonic mean method:

Equation (8) can be rewritten as

After the enthalpy distribution of the (*n* + 1)^{th} time step is obtained, the temperature distribution of the (*n* + 1)^{th} time zone can be obtained from eq. (4), i.e.,

The initial and boundary conditions are discretized as

The location of the solid-liquid interface at the (*n* + 1)^{th} time zone can easily be determined according to the temperature distribution. For example, if the enthalpy at a certain point *j* = *M* satisfies the solid-liquid interface will be

The solution procedure can be summarized as follows:

1. Determine the initial enthalpy at every node by using eq. (2).

2. Calculate the enthalpy after the first time step at nodes *j* = 2,...,*N* − 1 by using eq. (11).

3. Determine the temperature after the first time step at node *j* = 1,...,*N* by using eqs. (14), (12), and (15).

4. Find a control volume in which the enthalpy falls between 0 and and determine the location of the solid-liquid interface by using eq. (16).

5. Solve the phase-change problem at the next time step with the same procedure.

The solution procedure using the explicit scheme is attractive in its simplicity. However, there is a major drawback of this scheme – the limitation of stability. The limitation of stability should be familiar from prior experience in undergraduate heat transfer dealing with the numerical solution of unsteady state heat conduction. To satisfy the stability condition, the following Neumann stability criterion has to be satisfied:

where are thermal diffusivity of solid and liquid, respectively. Equation (17) is valid for one-dimensional problems only.

In addition to the above explicit scheme, one can use an implicit scheme – which is unconditionally stable – to overcome the drawback of potential instability of the explicit scheme. However, the solution process for an implicit scheme will be more complex than the solution for an explicit scheme because two unknown variables (enthalpy and temperature) are involved. The detailed procedure for the implicit scheme can be found in Crank (1984) or Alexiades and Solomon (1993).
The above enthalpy model treated enthalpy as a dependent variable in addition to the temperature, and discretized the energy equa¬tion into a set of equations which contain both *h* and *T*. For the implicit schemes, they actually treated all of the terms containing *T* = *T*(*h*) as a constant heat source term in the energy equation during iterations at each time step. This may cause some problems with convergence when *T* = *T*(*h*) is complicated and physical properties change significantly, as is the case for frozen heat pipe start-up, or when the boundary conditions are severe. Furthermore, when the energy equation contains a convective term, the previous methods have difficulties in handling the relationship between the convective and diffusive terms because of the two-dependent-variable nature of the equation.
Cao and Faghri (1989) proposed a simple strategy for trans¬forming the energy equation into a nonlinear equation with a single dependent variable *h*. Thus, solving a phase-change problem is equivalent to solving a non¬linear enthalpy equation, and existing algorithms are readily applicable with some modifications. The energy equation governing three-dimensional incompressible laminar flow with no viscous dissipation and with incorporation of the continuity equation in the Cartesian coordinate system is

with the state equation

In the cases with constant specific heats for each phase, in which phase change occurs at a single temperature, the temperature and enthalpy are related by eq. (4). Substituting eq. (4) into eq. (18), the energy equation can be transformed into the following form (Cao and Faghri, 1989):

where

**Failed to parse (syntax error): \Gamma (h) = \left\{ \begin{array}{l} {k_s}/{c_{ps}} & & h \le 0 \\ 0 & & 0 < h < {h_{s\ell }} \\ {k_\ell }/{c_{p\ell }} & & h > {h_{s\ell }} \\ \end{array} \right \qquad \qquad(21)**

**Failed to parse (syntax error): S(h) = \left\{ \begin{array}{l} 0 & & h \le 0 \\ 0 & & 0 < h < {h_{s\ell }} \\ h{k_\ell }/{c_{p\ell }} & h > {h_{s\ell }} \\ \end{array} \right \qquad \qquad(22)**

Cao and Faghri (1989) demonstrated that the above enthalpy transforming model is capable of handling complicated phase-change problems occurring both at a single temperature and over a temperature range with fixed grids. Due to the one dependent variable nature of the transformed equation, the convection and diffusion situations can be handled with appropriate algorithms. Comparisons have been made with the numerical results existing in the literature with a good agreement, showing that the model can prop¬erly predict the phase-change processes. The advan¬tages of this model based on enthalpy are that it allows us to avoid paying explicit attention to the nature of the phase-change front, and that it can be extended to accommodate complicated multidimensional problems with con¬vective terms without involving cumbersome math¬ematical schemes.

### Equivalent Heat Capacity Method

During the solid-liquid phase-change process, the PCM can absorb or release heat at a constant temperature *T*_{m}. This means that the temperature of the PCM does not change while it absorbs or releases heat, implying that the heat capacity of the PCM at temperature *T*_{m} is infinite. In the equivalent heat capacity method, it is assumed that the melting or solidification processes occur over a temperature range (*T*_{m} − Δ*T*, *T*_{m} + Δ*T*) instead of at a single temperature *T*_{m}. For a multicomponent system, Δ*T* can be chosen based on the range of phase change

temperature. For a single-component with well-defined melting point, Δ*T* should be as small as possible. Also, the latent heat is converted to an equivalent heat capacity of the PCM in the assumed temperature range. Thus, the specific heat of the PCM can be expressed as

**Failed to parse (syntax error): {c_p}(T) = \left\{ \begin{array}{l} {c_{ps}} & & & T < {T_m} - \Delta T \\ \frac{{{h_{s\ell }}}}{{2\Delta T}} + \frac{{{c_{ps}} + {c_{p\ell }}}}{2} & {T_m} - \Delta T < T < {T_m} + \Delta T \\ {c_{p\ell }} & & & T > {T_m} + \Delta T \\ \end{array} \right \qquad \qquad(23)**

which assumes that the temperature of the PCM is changed from *T*_{m} − Δ*T* to *T*_{m} + Δ*T* when latent heat, , is absorbed by the PCM during melting. During the solidification process, the PCM releases the latent heat and its temperature decreases from *T*_{m} + Δ*T* to *T*_{m} − Δ*T*. The equivalent specific heat in the mushy zone (*T*_{m} − Δ*T* < *T* < *T*_{m} + Δ*T*) includes the effect of both latent heat (the first term) and sensible heat (the second term). The relationship between specific heat and temperature in the equivalent heat capacity method is plotted in Fig. 4. For a three-dimensional conduction-controlled melting/solidification problem in the Cartesian coordinate system, the energy equation for the entire region of the PCM can be expressed as

where the thermal conductivity in eq. (24) is a function of temperature, *T*. Assuming that the thermal conductivity of the PCM in the two-phase region is a linear function of temperature, one obtains

**Failed to parse (syntax error): k(T) = \left\{ \begin{array}{l} {k_s} & & & {\rm{ }}T < {T_m} - \Delta T \\ {k_s} + \frac{{{k_\ell } - {k_s}}}{{2\Delta T}}(T - {T_m} + \Delta T){\rm{ }}{T_m} - \Delta T < T < {T_m} + \Delta T \\ {k_\ell } & & & {\rm{ }}T > {T_m} + \Delta T \\ \end{array} \right \qquad \qquad(25)**

The advantage of the equivalent heat capacity model is its simplicity. Equation (24) is simply the nonlinear heat conduction equation, and it appears that a conventional computational methodology for conduction problems is adequate for solving a solid-liquid phase change problem. However, many studies have revealed difficulties in the selection of time step, Δt, grid size, (Δx, Δy, Δz) and the phase change temperature range, ΔT. If these variables cannot be properly selected, the predicted location of the solid-liquid interface and the temperature may include some unrealistic oscillation. Therefore, although the equivalent heat capacity model leads to simple code development, it is not used as widely as the enthalpy model.

### Temperature-Transforming Model

The temperature-transforming model proposed by Cao and Faghri (1990) combines the advantages of the enthalpy and equivalent heat capacity models. For a three-dimensional conduction-controlled phase change problem, the governing equation in enthalpy form is

For a phase change occurring over a temperature range (*T*_{m} − Δ*T*, *T*_{m} + Δ*T*) with the specific heats assumed to be constant for each phase, the relationship between enthalpy and temperature can be plotted as in Fig. 5. This relationship can be analytically expressed as (see Problem 3.73)

By defining specific heat in the mushy zone as

eq. (27) becomes

**Failed to parse (syntax error): h(T) = \left\{ \begin{array}{l} {c_{ps}}(T - {T_m}) + {c_{ps}}\Delta T & & & T < {T_m} - \Delta T \\ \left( {{c_m} + \frac{{{h_{s\ell }}}}{{2\Delta T}}} \right)(T - {T_m}) + {c_m}\Delta T + \frac{{{h_{s\ell }}}}{2} & {T_m} - \Delta T < T < {T_m} + \Delta T \\ {c_{p\ell }}(T - {T_m}) + {c_{ps}}\Delta T + {h_{s\ell }} & & {\rm{ }}T > {T_m} + \Delta T \\ \end{array} \right \qquad \qquad(29)**

Equation (29) can be rewritten as

where *c*_{p}(*T*) and *b*(*T*) can be determined from eq. (29):

**Failed to parse (syntax error): {c_p}(T) = \left\{ \begin{array}{l} {c_{ps}} & & & T < {T_m} - \Delta T \\ {c_m} + \frac{{{h_{s\ell }}}}{{2\Delta T}} & & {T_m} - \Delta T < T < {T_m} + \Delta T \\ {c_{p\ell }} & & & T > {T_m} + \Delta T \\ \end{array} \right \qquad \qquad(31)**

**Failed to parse (syntax error): b(T) = \left\{ \begin{array}{l} {c_{ps}}\Delta T & & T < {T_m} - \Delta T \\ {c_m}\Delta T + \frac{{{h_{s\ell }}}}{2} & {T_m} - \Delta T < T < {T_m} + \Delta T \\ {c_{ps}}\Delta T + {h_{s\ell }} & T > {T_m} + \Delta T \\ \end{array} \right \qquad \qquad(32)**

Substituting eq. (30) into eq. (26) yields

where the thermal conductivity, *k*, is a function of temperature and can be obtained by eq. (25). The energy equation has been transformed into a nonlinear equation with a single dependent variable *T*. The nonlinearity of the phase change problem results from the large amount of heat released or absorbed at the solid-liquid interface. Also, in either the liquid or solid phases at some distance away from the interface, eq. (33) reduces to a linear equation.
The temperature-transforming model and the equivalent heat capacity model differ in significant ways. Comparing eq. (33) and eq. (24), one can conclude that the equivalent heat capacity model is a special case of eq. (33), with and *c*_{p} independent of the spatial variables *x*,*y*,*z* and time. This is the underlying reason why many studies using the equivalent heat capacity model have encountered difficulty in selecting the grid size and time step and have often produced physically unrealistic oscillatory results. In order to satisfy and the time step has to be small enough to assure that *c*_{p} and *b* are independent of both time and space variables – a difficult criterion to satisfy.
Equation (33) can be solved by many numerical methods, including the finite volume approach by Patankar (1980). To verify the accuracy of the temperature-transforming model, a numerical solution of the Stefan problem is presented. The physical model of the Stefan problem has been given in the previous section. Figure 6 shows a comparison between the exact solution and the numerical solution using the temperature transforming model for the dimensionless location of the melting front as a function of the square root of dimensionless time. The definition of dimensionless melting front location and dimensionless time is eq. from Exact Solutions of Melting and Solidification Problems.

The temperature-transforming model eliminates the time-step and grid-size limitations and is insensitive to the phase-change temperature range. Therefore, the model can properly handle phase change occurring over a temperature range (as in the case of multicomponent PCM phase change) or at a single temperature (as in the case of single-component PCM phase change). The temperature-transforming model has been used to solve many different solid-liquid phase change problems (Cao and Faghri 1990a, 1991, 1992; Zeng and Faghri 1994a, 1994b; Zhang and Faghri 1996b).

## References

Bonacina, C., Comini, G., Fasano, A., and Primicerio, A., 1973, “Numerical Solution of Phase Change Problems,” *International Journal of Heat and Mass Transfer*, Vol. 16, pp. 1285-1832.

Budhia, H., and Kreith, F., 1973, “Heat Transfer with Melting or Freezing in a Wedge,” *International Journal of Heat and Mass Transfer*, Vol. 16, pp. 195-211.

Cao, Y., and Faghri, A., 1990, “A Numerical Analysis of Phase Change Problem including Natural Convection,” ASME *Journal of Heat Transfer*, Vol. 112, pp. 812-815.

Cao, Y., and Faghri, A., 1989, “A Numerical Analysis of Stefan Problem of Generalized Multi-Dimensional Phase-Change Structures Using the Enthalpy Transforming Model,” *International Journal of Heat and Mass Transfer*, Vol. 32, pp. 1289-1298.

Crank, J., 1984, *Free and Moving Boundary Problems*, Clarendon Press, Oxford.

Poots, G., 1962, “An Approximate Treatment of Heat Conduction Problem Involving a Two-Dimensional Solidification Front,” *International Journal of Heat and Mass Transfer*, Vol. 5, pp. 339-348.

Rathjen, K.A., and Jiji, L.M., 1971, “Heat Conduction with Melting or Freezing in a Corner,” ASME *Journal of Heat Transfer*, Vol. 93, pp. 101-109.

Shamsunder, N., and Sparrow, E.M., 1975, “Analysis of Multidimensional Conduction Phase Change via the Enthalpy Model,” ASME *Journal of Heat Transfer*, Vol. 97, pp. 333-340.

Zhang, Y., and Faghri, A., 1996b, “Heat Transfer Enhancement in Latent Heat Thermal Energy Storage System by Using an External Radial Finned Tube,” *Journal of Enhanced Heat Transfer*, Vol. 3, pp. 119-127.