# Integral approximate solution of melting and solidification problems

### From Thermal-FluidsPedia

The melting and solidification problems that can be solved by exact solution are very limited, so it is necessary to introduce some approximate solution techniques. The integral approximate method is a very effective approach to solve melting/solidification problems.

## Contents |

## One-Region Problem

### Boundary condition of the 1st kind

The dimensionless governing equations for melting and solidification problems have the same form if the dimensionless variables are defined in exact solutions of melting and solidification problems were used. The procedure for solving one-region phase change problems will be demonstrated by solving a one-region melting problem. To solve the problem by integral approximate solution, the thermal penetration depth must be specified. For a one-region melting problem, only the temperature distribution in the liquid phase needs to be solved, because the temperature in the solid phase remains uniformly equal to the melting point of the PCM. Therefore, the thickness of the liquid phase is identical to the thermal penetration depth. Integrating the energy equation:

with respect to *X* in the interval of (0,*S*), one obtains

where

Substituting the following equation for energy balance at interface

into eq. (1) yields the integral equation of the one-region problem:

Assume now that the temperature distribution in the liquid phase is the following second-degree polynomial:

where *A*_{0}, *A*_{1}, and *A*_{2} are three unknown constants. The following two equations

can be used to determine the constants in eq. (4). Another appropriate boundary condition at the solid-liquid interface is needed.

Differentiating

one obtains

i.e.,

Substituting the following two equations

into the above equation yields

Equation (5) is an additional boundary condition at the solid-liquid interface; it can be used to determine the coefficients in eq. (4). After the constants in eq. (4) are determined, the temperature distribution in the liquid phase becomes

Substituting eq. (6) into the integral equation (3) leads to an ordinary differential equation for *S*(τ):

The initial condition of eq. (7) is

The right-hand side of eq. (7) is a constant and its solution is

where

Figure 1 shows the comparison of λ obtained by the exact solution and the integral approximate solution. The integral approximate solution agrees very well for a small Stefan number, but the difference increases along with increasing Stefan number. For a latent heat thermal energy storage system, the Stefan number is usually less than 0.2, so the integral approximate solution can provide sufficiently accurate results for that case.

### Boundary Condition of the 2nd kind

A solid PCM with a uniform initial temperature at its melting point, *T*_{m}, is in a half-space, *x* > 0. At time *t* = 0, a variable heat flux, *q*''_{0}(*t*), is suddenly applied to the surface of the semi-infinite body. Assume that the densities of the PCM for both phases are the same, and that natural convection in the liquid phase is negligible. We will find the transient location of the solid-liquid interface.

This problem is the same as the exact solutions of melting and solidification problems under the boundary condition of the 2nd kind except the heat flux at the surface is a function of time, i.e.,

where *q*''_{0} is a reference heat flux and *f*(*t*) is a given function. The dimensionless governing equation and the corresponding initial and boundary conditions of the problem are

where the nondimensional variables in the above eqs. (12) – (15) are defined by

Integrating eq. (12) with respect to *X* in the interval (0,*S*), and considering the boundary conditions at the surface, the integral equation of the problem is obtained:

where

Integrating eq. (16) with respect to τ in the interval (0, τ) yields

Assuming that the temperature distribution in the liquid phase is

where the constants, *A*_{0},*A*_{1}, and *A*_{2} are three unspecified constants that can be determined from eqs. (13), (14) and (5). After all unknown constants are determined, the temperature distribution in the liquid phase becomes

where

Substituting eq. (19) into eq. (18), one obtains

For the case of constant heat flux, i.e., *f*(τ) = 1, eq. (21) can be simplified to

## Two-Region Problem

The physical model of a melting problem is shown in Fig. 2: a solid PCM with a uniform initial temperature *T*_{i}, which is below its melting point *T*_{m}, is in a half-space, *x* > 0. At time *t* = 0, a constant heat flux, *q*''_{0}, is suddenly applied to the surface of the semi-infinite body. Because the initial temperature of the PCM is below its melting point, melting does not begin until after the wall temperature reaches the melting point. Therefore, the problem can be divided into two sub-problems: (1) heat conduction over the duration of preheating, and (2) the actual melting process. An integral approximate method will be employed to solve both sub-problems (Zhang et al., 1993).

### Duration of Preheating

At the beginning of heating, no melting occurs, and the problem is a pure conduction problem with a boundary condition of the second kind. Its mathematical description is as follows:

where *t*_{m} is the duration of preheating. This problem is solved by integral approximate method. Assuming the temperature profile is a second-degree polynomial, one obtains the temperature profile:

where δ is the thermal penetration depth. It can be obtained by substituting eq. (27) into the integral equation, which can in turn be obtained by integrating eq. (23) in the interval of (0, δ). The result is

The highest temperature of the semi-infinite body occurs at its surface (*x* = 0) and can be expressed as

Melting occurs when the surface temperature reaches the melting point, *T*_{m}, and the corresponding thermal penetration depth is

Then the duration of preheating can be obtained by combining eqs. (28) and (30),

The temperature distribution at time *t*_{m} is

### Governing Equations for the Melting Stage

After melting starts, the governing equations in the different phases must be specified separately. The temperature in the liquid phase satisfies

The governing equation and corresponding boundary and initial conditions for the solid phase are

At the solid-liquid interface, the following boundary conditions are necessary to link solutions in the liquid and solid phases:

By defining the dimensionless variables as follows:

where *S*_{c} is a subcooling parameter that signifies the ratio of the sensible preheat needs to raise the temperature of the body before melting takes place, eqs. (33) – (39) become

where Δ_{m} and τ_{m} can be obtained by substituting eq. (40) into eqs. (30) and (31), i.e.,

### Integral Approximate Solution

Figure 3 shows the physical model represented by the above dimensionless governing equations. After the melting process starts, the solid-liquid interface moves along the positive *x*-axis, and the thermal penetration depth continuously increases. The solid phase temperature over the range of *S*(τ) < *X* < Δ(τ) is affected by the boundary condition at the surface. The temperature in the region beyond the thermal penetration depth, Δ, is not affected and remains at –*S**c*.

Integrating eq. (43) over the interval (*S*,Δ), then applying the definition of the thermal penetration depth and eq. (46), yields the integral equation of the solid phase:

where

The temperature distribution in the solid phase is assumed to be a second-degree polynomial, i.e.,

_{2}=

*B*

_{0}+

*B*

_{1}(

*X*−

*S*) +

*B*

_{2}(

*X*−

*S*)

^{2}

where the constants *B*_{0},*B*_{1}, and *B*_{2} can be obtained from the boundary condition, eq. (46) and the definition of the thermal penetration depth ( and The final form of the temperature distribution in the solid phase is

Substituting eq. (52) into eq. (50) yields a relationship between the location of the solid-liquid interface and the thermal penetration depth:

Integrating the differential equation of liquid phase eq. (41) over the interval of (0,*S*) and applying the boundary conditions eqs. (42) and (46) yields the integral equation of the liquid phase:

where

Assuming that the temperature profile in the liquid phase is a second-degree polynomial function,

where *A*_{0},*A*_{1}, and *A*_{2} need to be determined using appropriate boundary conditions. Equations (42) and (46) provided two conditions; and the third condition needs to be obtained by differentiating eq. (46), i.e.,

i.e.,

Substituting eqs. (41) and (47) into eq. (57) yields

which is an additional boundary condition at the solid-liquid interface, and can be used to determine the coefficients in eq. (56). After the constants in eq. (56) are determined, the temperature distribution in the liquid phase becomes

where

Substituting the integral eqs. (54) and (50) into eq. (47) yields

Integrating both sides of eq. (61) with respect to τ within (τ_{m},τ), one obtains

Substituting eqs. (51), (52), (55), and (59) into eq. (62), the resulting interface equation is

Equation (53) is a first-order ordinary differential equation, while eq. (63) is an algebraic equation. They can be solved simultaneously with an implicit numerical method to determine the location of the solid-liquid interface and the thermal penetration depth. If there is no subcooling, Sc=0, the solid-liquid interface location can be obtained by the following expression:

## Ablation under Constant Heat Flux Heating

Ablation is defined as the removal of material from the surface of an object by vaporization, chipping, or other erosive process. One application is when tire manufacturers design tires for automotive race events, such as NASCAR races. They have to balance the expected tire wear on the track versus the amount of heat that will be removed via ablation. If they underestimate the wear, the tire overheats and fails as too little heat is removed. If they overestimate the wear, the tire temperature never reaches its optimal operating condition to maximize adhesion. Ablation is also an effective means of protecting the surfaces of missiles and space shuttles from high-rate aerodynamic heating during atmospheric re-entry. It is a sacrificial cooling method, because the protective layer is partially destroyed. The advantage of the ablative cooling process is its self-regulation: the rate of ablation is automatically adjusted in response to the heating rate. The most commonly used materials for ablative cooling are PCMs with higher melting points (such as glass, carbon, or polymer fiber) in combination with an organic binder.

From the heat transfer point of view, ablation is a special case of melting in which the liquid phase is completely removed immediately upon its production. Therefore, ablation can be considered as a one-region melting problem where heating occurs directly on a solid-liquid interface. If the spacecraft maintains a constant velocity during re-entry, one can assume that the entire process occurs under constant heat flux heating. Similar to the melting of a subcooled solid under constant heat flux, ablation does not start simultaneously with heating, because the initial temperature of the ablating material, *T*_{i}, is usually below its melting point, *T*_{m}. The melting will start only after a period of preheating that allows the surface temperature to reach the melting point of the ablating material, *T*_{m}. The preheating problem is a pure conduction problem and can be solved analytically. Since the surface of the ablating material is constantly moving and the liquid phase does not accumulate, ablation can be described by using a coordinate system, x, that is attached to and moves with the ablating velocity, Ua (see Fig. 4). The ablation material moves with a velocity –Ua in the moving coordinate system. The energy equation in a moving coordinate is (Eckert and Drake, 1987):

Since the ablating materials usually have very low thermal conductivity, it is reasonable to assume that the ablation occurs in a semi-infinite body. Therefore, the initial and boundary conditions of eq. (65) are

The energy balance at the surface is

Shortly after ablation begins, the process enters steady-state and the ablation velocity becomes a constant. Equation (65) can be simplified as

where the thermal properties are assumed to be independent of temperature. Equation (70) can be treated as a first-order differential equation of *d**T* / *d**x*, and its solution is

where *C*_{1} is an integrating constant. Equation (71) can be further integrated to yield

The constants *C*_{1} and *C*_{2} can be determined from eqs. (67) – (68), and the temperature distribution in the ablating material becomes

Substituting eq. (73) into eq. (69) yields

which can be rearranged to obtain the ablation velocity:

where *S*_{c} is the subcooling parameter defined as

The fraction of heat removed by ablation over the total heat flux is

Substituting the ablation velocity from eq. (75) into eq. (77), one obtains

If the subcooling parameter is *S*_{c} = 0.2, 83% of the total heat flux will be removed by ablation.

## Solidification/Melting in Cylindrical Coordinate Systems

### Boundary condition of the 1st kind

Application of the integral approximate solution to one-dimensional solid-liquid phase change problems – including ablation – in a Cartesian coordinate system has been discussed in the preceding sections. Since tubes are widely used in shell-and-tube thermal energy storage devices, it is necessary to study melting and solidification in cylindrical coordinate systems as well. The polynomial temperature distribution is a very good approximation of the one-dimensional problem in the Cartesian coordinate system, but it can result in very significant error if it is used to solve for the phase change heat transfer in a cylindrical coordinate system. This is because heat transfer area for a cylindrical coordinate system varies with the coordinate *r* instead of remaining constant. Therefore, the temperature distribution in the cylindrical coordinate system has to be modified by taking into account the variation of the heat transfer area.
Solidification around a cylinder with radius of *r*_{i}, as shown in Fig. 5, will be investigated in order to demonstrate application of the integral approximate solution in the cylindrical coordinate system. An infinite liquid PCM has a uniform initial temperature equal to the melting point of the PCM, *T*_{m}. At time *t* = 0, the inner surface of the cylinder suddenly decreases to a temperature *T*_{0}, which is below the melting point of the PCM. Since the temperature of the liquid PCM equals the melting point of the PCM and the temperature in the solid phase is unknown, this is a one-region solidification problem. Conduction controls the solidification process, because the temperature in the liquid phase is uniformly equal to the melting point of the PCM. The governing equation and the initial and boundary conditions of this problem are

where the subscript 1 for solid phase has been dropped for ease of notation.
Defining dimensionless variables,

eqs. (79) – (82) become

The above eqs. (84) – (87) are also valid for melting around a hollow cylinder when used with the appropriate dimensionless variables.
From elementary heat transfer, it is known that the logarithmic function is the exact solution of the steady-state heat conduction problem in a cylindrical wall. Therefore, we can assume that the temperature distribution has a second-order logarithmic function of the form (Zhang and Faghri, 1996a):

where is an unknown variable. Equation (88) automatically satisfies eqs. (85) and (86), and can be obtained by differentiating eq. (86)

which can be rearranged to get

Substituting eqs. (84) and (86) into eq. (89) yields the following expression:

Substituting eq. (88) into eq. (90) gives the following expression for :

Substituting eqs. (88) and (91) into eq. (87) leads to an ordinary differential equation for the locations of the solid-liquid interface

which is subjected to the following initial condition:

Integrating equation (92) over the time interval (0,τ) results in the following equation for the location of the solid-liquid interface:

which is valid for solidification around a cylinder with a constant inner temperature. However, in practical applications, cooling inside the tube is usually achieved by the flow of cooling fluid through the tube. Therefore, the boundary condition at the inner surface of the tube should be a convective boundary condition instead of an isothermal boundary condition. The following example demonstrates application of the integral approximate method to the solution of solidification around a cylindrical tube convectively cooled from inside.

### Boundary condition of the 3rd kind

An infinite liquid PCM has a uniform initial temperature equal to the melting point of the PCM, *T*_{m}. At time *t* = 0, a cooling fluid with temperature *T*_{i} flows inside the tube. The heat transfer coefficient between the cooling fluid and the inner surface of the tube is *h*_{i}. We can find the transient location of the solid-liquid interface.

The governing equations of the problem can also be represented by eq. (79) – (82) except that eq. (80) needs to be replaced by the following expression:

The dimensionless governing equations of the problem are the same as eqs. (84) – (87), but eq. (85) needs to be replaced by the dimensionless form of eq. (95), i.e.,

where *B**i* in eq. (96) is the Biot number defined as

and θ is defined as

It is also assumed that the temperature distribution has a second-order logarithmic function of the form

The two unknown variables *A* and *B* in eq. (99) can be determined by substituting eq. (99) into eqs. (96) and (90), with the result

Substituting eq. (100) into eq. (101), an equation for *A* is obtained as follows:

Substituting eq. (99) into eq. (87), the ordinary differential equation for the location of the solid-liquid interface is obtained:

which is subject to the initial condition specified by eq. (93).
The temperature of the inner surface of the tube is

The solid-liquid interface location can be obtained by numerical solution of eq. (103) and (104). It should be noted that A becomes 1 if the Biot number becomes infinite [see eq. (102)]. In that case the temperature of the inner surface becomes 1 and eq. (103) reduces to eq. (92).

## References

Eckert, E.R.G., and Drake, R.M., 1987, *Analysis of Heat and Mass Transfer*, Hemisphere, Washington, DC.

Faghri, A., and Zhang, Y., 2006, *Transport Phenomena in Multiphase Systems*, Elsevier, Burlington, MA.

Faghri, A., Zhang, Y., and Howell, J. R., 2010, *Advanced Heat and Mass Transfer*, Global Digital Press, Columbia, MO.

Zhang, Y., Chen, Z., and Wang Q., 1993, “Analytical Solution of Melting in a Subcooled Semi-Infinite Sold with Boundary Conditions of the Second Kind,” *Journal of Thermal Science*, Vol. 2, pp. 111-115, Science Press, Beijing.

Zhang, Y., and Faghri, A., 1996a, “Semi-Analytical Solution of Thermal Energy Storage System with Conjugate Laminar Forced Convection,” *International Journal of Heat and Mass Transfer*, Vol. 39, pp. 717-724.