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.

1 One-Region Problem
- 1.1 Boundary condition of the 1st kind
- 1.2 Boundary Condition of the 2nd kind
2 Two-Region Problem
3 Ablation under Constant Heat Flux Heating
4 Solidification/Melting in Cylindrical Coordinate Systems
- 4.1 Boundary condition of the 1st kind
- 4.2 Boundary condition of the 3rd kind
5 References
6 Further Reading
7 External Links

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:

$\frac{{{\partial ^2}\theta }}{{\partial {X^2}}} = \frac{{\partial \theta }}{{\partial \tau }}\quad \quad 0 < X < S(\tau ),\quad \tau > 0$

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

$\mathop {\left. {\frac{{\partial \theta (X,\tau )}}{{\partial X}}} \right|}_{X = S(\tau )} - \mathop {\left. {\frac{{\partial \theta (X,\tau )}}{{\partial X}}} \right|}_{X = 0} = \frac{{d\Theta (\tau )}}{{d\tau }} \qquad \qquad(1)$

where

$\Theta (\tau ) = \int_0^{S(\tau )} \,\theta (X,\tau )dX \qquad \qquad(2)$

Substituting the following equation for energy balance at interface

$- \frac{{\partial \theta }}{{\partial X}} = \frac{1}{{{\rm{Ste}}}}\frac{{dS}}{{d\tau }}\quad \quad X = S(\tau ),\quad \tau > 0$

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

$- \frac{1}{{Ste}}\frac{{dS}}{{d\tau }} - \mathop {\left. {\frac{{\partial \theta (X,\tau )}}{{\partial X}}} \right|}_{X = 0} = \frac{{d\Theta (\tau )}}{{d\tau }} \qquad \qquad(3)$

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

$\theta (X,\tau ) = {A_0} + {A_1}\left( {\frac{{X - S}}{S}} \right) + {A_2}\mathop {\left( {\frac{{X - S}}{S}} \right)}^2 \qquad \qquad(4)$

where $A 0$ , $A 1$ , and $A 2$ are three unknown constants. The following two equations

$\theta (X,\tau ) = 1\quad \quad X = 0(\tau ),\quad \tau > 0$ $\theta (X,\tau ) = 0\quad \quad X = S(\tau ),\quad \tau > 0$

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

Differentiating

$\theta (X,\tau ) = 0\quad \quad X = S(\tau ),\quad \tau > 0$

one obtains

$d\theta = \frac{{\partial \theta }}{{\partial X}}dX + \frac{{\partial \theta }}{{\partial \tau }}d\tau = 0\quad \quad X = S(\tau )$

i.e.,

$\frac{{\partial \theta }}{{\partial X}}\frac{{dS}}{{d\tau }} + \frac{{\partial \theta }}{{\partial \tau }} = 0\quad \quad X = S(\tau )$

Substituting the following two equations

$\frac{{{\partial ^2}\theta }}{{\partial {X^2}}} = \frac{{\partial \theta }}{{\partial \tau }}\quad \quad 0 < X < S(\tau ),\quad \tau > 0$ $- \frac{{\partial \theta }}{{\partial X}} = \frac{1}{{{\rm{Ste}}}}\frac{{dS}}{{d\tau }}\quad \quad X = S(\tau ),\quad \tau > 0$

into the above equation yields

$Ste\mathop {\left( {\frac{{\partial \theta }}{{\partial X}}} \right)}^2 = \frac{{{\partial ^2}\theta }}{{\partial {X^2}}} \qquad \qquad(5)$

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

$\theta (X,\tau ) = \frac{{1 - \sqrt {1 + 2Ste} }}{{Ste}}\left( {\frac{{X - S}}{S}} \right) + \left( {\frac{{1 - \sqrt {1 + 2Ste} }}{{Ste}} + 1} \right)\mathop {\left( {\frac{{X - S}}{S}} \right)}^2 \qquad \qquad(6)$

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

Figure 1: Comparison between integral and exact solutions of the one-region conduction controlled melting problem

$S\frac{{dS}}{{d\tau }} = 6\frac{{1 - \sqrt {1 + 2Ste} + 2Ste}}{{5 + \sqrt {1 + 2Ste} + 2Ste}} \qquad \qquad(7)$

The initial condition of eq. (7) is

$S(\tau ) = 0\quad \quad \tau = 0 \qquad \qquad(8)$

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

$S = 2\lambda {\tau ^{{\textstyle{1 \over 2}}}} \qquad \qquad(9)$

where

$\lambda = \mathop {\left[ {3\frac{{1 - \sqrt {1 + 2Ste} + 2Ste}}{{5 + \sqrt {1 + 2Ste} + 2Ste}}} \right]}^{{\textstyle{1 \over 2}}} \qquad \qquad(10)$

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

$q''(t) = f(t){q''_0} \qquad \qquad(11)$

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

$\frac{{{\partial ^2}\theta }}{{\partial {X^2}}} = \frac{{\partial \theta }}{{\partial \tau }}\quad \quad 0 < X < S(\tau ),\quad \tau > 0 \qquad \qquad(12)$

$\frac{{\partial \theta }}{{\partial X}} = - f(\tau )\quad \quad X = 0,\quad \tau > 0 \qquad \qquad(13)$

$\theta (X,\tau ) = 0\quad \quad X = S(\tau ),\quad \tau > 0 \qquad \qquad(14)$

$- \frac{{\partial \theta }}{{\partial X}} = \frac{1}{{Ste}}\frac{{dS}}{{d\tau }}\quad \quad X = S(\tau ),\quad \tau > 0 \qquad \qquad(15)$

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

$\begin{array}{*{20}{c}} {\theta = \frac{{T - {T_m}}}{{{{q''}_0}L/k}},} & {X = \frac{x}{L}} \\ \end{array}\begin{array}{*{20}{c}} , & {S = \frac{s}{L}\begin{array}{*{20}{c}} , & {\tau = \frac{{\alpha t}}{{{L^2}}}\begin{array}{*{20}{c}} , & {{\rm{Ste}} = \frac{{{c_p}({{q''}_0}L/k)}}{{{h_{s\ell }}}}} \\ \end{array}} \\ \end{array}} \\ \end{array}$

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:

$- \frac{1}{{Ste}}\frac{{dS(\tau )}}{{d\tau }} + f(\tau ) = \frac{{d\Theta (\tau )}}{{d\tau }} \qquad \qquad(16)$

where

$\Theta (\tau ) = \int_0^{S(\tau )} \,\theta (X,\tau )dX \qquad \qquad(17)$

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

$\frac{{S(\tau )}}{{Ste}} + \Theta (\tau ) = \int_0^\tau \,f(\tau )d\tau \qquad \qquad(18)$

Assuming that the temperature distribution in the liquid phase is

$\theta (X,\tau ) = {A_0} + {A_1}\left( {\frac{{X - S}}{S}} \right) + {A_2}\mathop {\left( {\frac{{X - S}}{S}} \right)}^2$

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

$\theta (X,\tau ) = \frac{1}{2}\frac{{\left[ {1 - {{(1 + 4\mu )}^{1/2}}} \right]}}{{Ste}}\left( {\frac{{X - S}}{S}} \right) + \frac{1}{8}\frac{{\mathop {\left[ {1 - {{(1 + 4\mu )}^{1/2}}} \right]}^2 }}{{Ste}}\mathop {\left( {\frac{{X - S}}{S}} \right)}^2 \qquad \qquad(19)$

where

$\mu = f(\tau )S(\tau )Ste \qquad \qquad(20)$

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

$St{e^2}f(\tau )\int_0^\tau \,f(\tau )d\tau = \frac{\mu }{6}\left[ {\mu + 5 + {{(1 + 4\mu )}^{1/2}}} \right] \qquad \qquad(21)$

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

$S\left[ {SteS + 5 + {{(1 + 4SteS)}^{1/2}}} \right] = 6Ste\tau \qquad \qquad(22)$

Two-Region Problem

Figure 2: Melting in a subcooled semi-infinite body under constant heat flux

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:

$\frac{{{\partial ^2}{T_2}}}{{\partial {x^2}}} = \frac{1}{{{\alpha _2}}}\frac{{\partial {T_2}}}{{\partial t}}\quad 0 < x < \infty ,\quad 0 < t < {t_m} \qquad \qquad(23)$

$\frac{{\partial {T_2}}}{{\partial x}} = - \frac{1}{{{k_2}}}{q''_0}\quad x = 0,\quad 0 < t < {t_m} \qquad \qquad(24)$

${T_2}(x,t) \to {T_i}\quad x \to \infty ,\quad 0 < t < {t_m} \qquad \qquad(25)$

${T_2}(x,t) = {T_i},\quad 0 < x < \infty ,\quad t = 0 \qquad \qquad(26)$

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:

${T_2}(x,t) = {T_i} + \frac{{{{q''}_0}\delta }}{{2{k_2}}}\mathop {\left( {1 - \frac{x}{\delta }} \right)}^2 \qquad \qquad(27)$

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

$\delta = \sqrt {6{\alpha _2}t} \qquad \qquad(28)$

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

${T_s}(t) = {T_i} + \frac{{{{q''}_0}\delta }}{{2{k_2}}} \qquad \qquad(29)$

Melting occurs when the surface temperature reaches the melting point, $T m,$ and the corresponding thermal penetration depth is

${\delta _m} = \frac{{2{k_2}({T_m} - {T_i})}}{{{{q''}_0}}} \qquad \qquad(30)$

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

${t_m} = \frac{{2k_2^2{{({T_m} - {T_i})}^2}}}{{3{\alpha _2}{{q''}_0}^2}} \qquad \qquad(31)$

The temperature distribution at time $t m$ is

${T_2}(x,{t_m}) = {T_i} + ({T_m} - {T_i})\mathop {\left( {1 - \frac{x}{{{\delta _m}}}} \right)}^2 \qquad \qquad(32)$

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

$\frac{{{\partial ^2}{T_1}}}{{\partial {x^2}}} = \frac{1}{{{\alpha _1}}}\frac{{\partial {T_1}}}{{\partial t}}\quad \quad 0 < x < s(t),\quad t > {t_m} \qquad \qquad(33)$

$\frac{{\partial {T_1}(x,t)}}{{\partial x}} = - \frac{1}{{{k_1}}}{q''_0}\quad \quad x = 0,\quad t > {t_m} \qquad \qquad(34)$

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

$\frac{{{\partial ^2}{T_2}}}{{\partial {x^2}}} = \frac{1}{{{\alpha _2}}}\frac{{\partial {T_2}}}{{\partial t}}\quad \quad s(t) < x < \infty ,\quad t > {t_m} \qquad \qquad(35)$

${T_2}(x,t) \to {T_i}\quad \quad x \to \infty ,\quad t > {t_m} \qquad \qquad(36)$

${T_2}(x,t) = {T_i} + ({T_m} - {T_i})\mathop {\left( {1 - \frac{x}{{{\delta _m}}}} \right)}^2 \quad \quad x > 0,\quad t = {t_m} \qquad \qquad(37)$

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

${T_1}(x,t) = {T_2}(x,t) = {T_m}\quad \quad x = s(t),\quad t > {t_m} \qquad \qquad(38)$

${k_2}\frac{{\partial {T_2}}}{{\partial x}} - {k_1}\frac{{\partial {T_1}}}{{\partial x}} = \rho {h_{s\ell }}\frac{{ds}}{{dt}}\quad \quad x = s(t),\quad t > {t_m} \qquad \qquad(39)$

By defining the dimensionless variables as follows:

$\left. \begin{array}{l} {\theta _1} = \frac{{{c_{p1}}({T_1} - {T_m})}}{{{h_{s\ell }}}}\,\quad \,{\theta _2} = \frac{{{c_{p2}}({T_2} - {T_m})}}{{{h_{s\ell }}}}\,\quad \,Sc = \frac{{{c_{p2}}({T_m} - {T_i})}}{{{h_{s\ell }}}} \\ X = \frac{x}{{{\alpha _1}{\rho _1}{h_{s\ell }}/{{q''}_0}}}\quad S = \frac{s}{{{\alpha _1}{\rho _1}{h_{s\ell }}/{{q''}_0}}}\,\quad \,\Delta = \frac{\delta }{{{\alpha _1}{\rho _1}{h_{s\ell }}/{{q''}_0}}} \\ {\Delta _m} = \frac{{{\delta _m}}}{{{\alpha _1}{\rho _1}{h_{s\ell }}/{{q''}_0}}}\,\,\quad \tau = \frac{{{\alpha _1}t}}{{{{({\alpha _1}{\rho _1}{h_{s\ell }}/{{q''}_0})}^2}}}\quad {N_\alpha } = \frac{{{\alpha _2}}}{{{\alpha _1}}} \\ \end{array} \right\} \qquad \qquad(40)$

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

$\frac{{{\partial ^2}{\theta _1}}}{{\partial {X^2}}} = \frac{{\partial {\theta _1}}}{{\partial \tau }}\quad \quad 0 < X < S(\tau ),\quad \tau > {\tau _m} \qquad \qquad(41)$

$\frac{{\partial {\theta _1}(X,\tau )}}{{\partial X}} = - 1\quad \quad X = 0,\quad \tau > {\tau _m} \qquad \qquad(42)$

$\frac{{{\partial ^2}{\theta _2}}}{{\partial {X^2}}} = \frac{1}{{{N_\alpha }}}\frac{{\partial {\theta _2}}}{{\partial \tau }}\quad \quad S(\tau ) < X < \infty ,\quad \tau > {\tau _m} \qquad \qquad(43)$

${\theta _2}(X,\tau ) \to - Sc\quad \quad X \to \infty ,\quad \tau > {\tau _m} \qquad \qquad(44)$

${\theta _2}(X,\tau ) = Sc\left[ {\mathop {\left( {1 - \frac{X}{{{\Delta _m}}}} \right)}^2 - 1} \right]\quad \quad X > 0,\quad \tau = {\tau _m} \qquad \qquad(45)$

${\theta _1}(X,\tau ) = {\theta _2}(X,\tau ) = 0\quad \quad X = S(\tau ),\quad \tau > {\tau _m} \qquad \qquad(46)$

$- \frac{{\partial {\theta _1}}}{{\partial X}} + \frac{1}{{{N_\alpha }}}\frac{{\partial {\theta _2}}}{{\partial X}} = \frac{{dS}}{{d\tau }}\quad \quad X = S(\tau ),\quad \tau > {\tau _m} \qquad \qquad(47)$

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

${\Delta _m} = 2{N_\alpha }Sc \qquad \qquad(48)$

${\tau _m} = \frac{2}{3}{N_\alpha }S{c^2} \qquad \qquad(49)$

Integral Approximate Solution

Figure 3: Dimensionless temperature distribution for melting in a semi-infinite body under constant heat flux

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:

$\mathop {\left. { - \frac{{\partial {\theta _2}}}{{\partial X}}} \right|}_{X = S} = \frac{1}{{{N_\alpha }}}\frac{d}{{d\tau }}\left( {{\Theta _2} + Sc\Delta } \right) \qquad \qquad(50)$

where

${\Theta _2} = \int_S^\Delta \,{\theta _2}dX \qquad \qquad(51)$

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 ( $\mathop {\left. {{\theta _2}} \right|}_{X = \Delta } = - Sc$ and $\mathop {\left. {\partial {\theta _2}/\partial X} \right|}_{X = \Delta } = 0).$ The final form of the temperature distribution in the solid phase is

${\theta _2}(X,\tau ) = Sc\left[ {\mathop {\left( {\frac{{\Delta - X}}{{\Delta - S}}} \right)}^2 - 1} \right] \qquad \qquad(52)$

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

$\frac{{6{N_\alpha }}}{{\Delta - S}} = 2\frac{{dS}}{{d\tau }} + \frac{{d\Delta }}{{d\tau }} \qquad \qquad(53)$

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:

$\mathop {\left. {\frac{{\partial {\theta _1}}}{{\partial X}}} \right|}_{X = S} + 1 = \frac{{d{\Theta _1}}}{{d\tau }} \qquad \qquad(54)$

where

${\Theta _1} = \int_0^S \,{\theta _1}dX \qquad \qquad(55)$

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

${\theta _1} = {A_0} + {A_1}\left( {\frac{{X - S}}{S}} \right) + {A_2}{\left( {\frac{{X - S}}{S}} \right)^2} \qquad \qquad(56)$

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

$d{\theta _1} = \frac{{\partial {\theta _1}}}{{\partial X}}dX + \frac{{\partial {\theta _1}}}{{\partial \tau }}d\tau = 0\quad \quad X = S(\tau )$

i.e.,

$\frac{{\partial {\theta _1}}}{{\partial X}}\frac{{dS}}{{d\tau }} + \frac{{\partial {\theta _1}}}{{\partial \tau }} = 0\quad \quad X = S(\tau ) \qquad \qquad(57)$

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

$\mathop {\left( {\frac{{\partial {\theta _1}}}{{\partial X}}} \right)}^2 - \frac{1}{{{N_\alpha }}}\frac{{\partial {\theta _1}}}{{\partial X}}\frac{{\partial {\theta _2}}}{{\partial X}} = \frac{{{\partial ^2}{\theta _1}}}{{\partial {X^2}}} \qquad \qquad(58)$

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

${\theta _1}(X,\tau ) = \frac{S}{2}\mathop {\left( {\frac{{X - S}}{S}} \right)}^2 - \frac{1}{2}p\frac{{{X^2} - {S^2}}}{{{S^2}}} \qquad \qquad(59)$

where

$p = \left( {\frac{{{N_\alpha }ScS}}{{\Delta - S}} - \frac{1}{2}} \right) + \sqrt {\mathop {\left( {\frac{{{N_\alpha }ScS}}{{\Delta - S}} - \frac{1}{2}} \right)}^2 + S} \qquad \qquad(60)$

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

$\frac{{dS}}{{d\tau }} + \frac{d}{{d\tau }}\left( {{\Theta _2} + Sc\Delta } \right) + \frac{{d{\Theta _1}}}{{d\tau }} = 1 \qquad \qquad(61)$

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

$S + [{\Theta _2}(\tau ) - {\Theta _2}({\tau _m})] + Sc(\Delta - {\Delta _m}) + [{\Theta _1}(\tau ) - {\Theta _1}({\tau _m})] = (\tau - {\tau _m}) \qquad \qquad(62)$

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

$\frac{1}{2}{S^2} + (p + 3 + 2Sc)S + Sc(\Delta - {\Delta _m}) = 3(\tau - {\tau _m}) \qquad \qquad(63)$

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:

$S\left( {S + 5 + \sqrt {1 + 4S} } \right) = 6\tau \qquad \qquad(64)$

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.

Figure 4: Physical model for ablation.

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):

$\frac{\partial }{{\partial x}}\left( {k\frac{{\partial T}}{{\partial x}}} \right) = - {U_a}\rho {c_p}\frac{{\partial T}}{{\partial x}} + \rho {c_p}\frac{{\partial T}}{{\partial t}} \qquad \qquad(65)$

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

$\begin{array}{*{20}{c}} {T = {T_i},} & {t = 0} \\ \end{array} \qquad \qquad(66)$ $T = {T_m}\begin{array}{*{20}{c}} , & {x = 0} \\ \end{array} \qquad \qquad(67)$ $T = {T_i}\begin{array}{*{20}{c}} , & {\begin{array}{*{20}{c}} {\frac{{\partial T}}{{\partial x}} = 0,} & {} \\ \end{array}x \to \infty } \\ \end{array} \qquad \qquad(68)$

The energy balance at the surface is

$q'' - \rho {U_a}{h_{s\ell }} = - k{\left. {\frac{{\partial T}}{{\partial x}}} \right|_{x = 0}}\begin{array}{*{20}{c}} , & {x = 0} \\ \end{array} \qquad \qquad(69)$

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

$\frac{{{d^2}T}}{{d{x^2}}} = - \frac{{{U_a}}}{\alpha }\frac{{dT}}{{dx}} \qquad \qquad(70)$

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

$\frac{{dT}}{{dx}} = {C_1}{e^{ - \frac{{{U_a}}}{\alpha }x}} \qquad \qquad(71)$

where $C 1$ is an integrating constant. Equation (71) can be further integrated to yield

$T = - \frac{{{C_1}\alpha }}{{{U_a}}}{e^{ - \frac{{{U_a}}}{\alpha }x}} + {C_2} \qquad \qquad(72)$

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

$T = {T_i} + ({T_m} - {T_i}){e^{ - \frac{{{U_a}}}{\alpha }x}} \qquad \qquad(73)$

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

$q'' - \rho {U_a}{h_{s\ell }} = \rho {c_p}{U_a}({T_m} - {T_i}) \qquad \qquad(74)$

which can be rearranged to obtain the ablation velocity:

${U_a} = \frac{{q''}}{{\rho {h_{s\ell }}(1 + Sc)}} \qquad \qquad(75)$

where $S c$ is the subcooling parameter defined as

$Sc = \frac{{{c_p}({T_m} - {T_i})}}{{{h_{s\ell }}}} \qquad \qquad(76)$

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

$\frac{{{{q''}_a}}}{{q''}} = \frac{{\rho {U_a}{h_{s\ell }}}}{{q''}} \qquad \qquad(77)$

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

$\frac{{{{q''}_a}}}{{q''}} = \frac{1}{{1 + Sc}} \qquad \qquad(78)$

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

Figure 5: Solidification of an infinite liquid PCM around an ID-cooled cylinder.

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

$\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\frac{{\partial T}}{{\partial r}}} \right) = \frac{1}{\alpha }\frac{{\partial T}}{{\partial t}}\quad {r_i} < r < s(t)\quad t > 0 \qquad \qquad(79)$

$T(r,t) = {T_0}\quad r = {r_i}\quad t > 0 \qquad \qquad(80)$

$T(r,t) = {T_m}\quad r = s\quad t > 0 \qquad \qquad(81)$

$k\frac{{\partial T}}{{\partial r}} = \rho {h_{s\ell }}\frac{{ds}}{{dt}}\quad r = s\quad t > 0 \qquad \qquad(82)$

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

$\theta = \frac{{{T_m} - T}}{{{T_m} - {T_0}}}\quad R = \frac{r}{{{r_i}}}\quad S = {\frac{s}{r}_i}\quad \tau = \frac{{\alpha t}}{{r_i^2}}\quad Ste = \frac{{{c_p}({T_m} - {T_0})}}{{{h_{s\ell }}}} \qquad \qquad(83)$

eqs. (79) – (82) become

$\frac{1}{R}\frac{\partial }{{\partial R}}\left( {R\frac{{\partial \theta }}{{\partial R}}} \right) = \frac{{\partial \theta }}{{\partial \tau }}\quad 1 < R < S(\tau )\quad \tau > 0 \qquad \qquad(84)$

$\theta (R,\tau ) = 1\quad R = 1\quad \tau > 0 \qquad \qquad(85)$

$\theta (R,\tau ) = 0\quad R = S(\tau )\quad \tau > 0 \qquad \qquad(86)$

$- \frac{{\partial \theta }}{{\partial R}} = \frac{1}{{Ste}}\frac{{dS}}{{d\tau }}\quad R = S\quad \tau > 0 \qquad \qquad(87)$

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):

$\theta = 1 + \varphi \left( {\frac{{lnR}}{{lnS}}} \right) - (1 + \varphi )\mathop {\left( {\frac{{lnR}}{{lnS}}} \right)}^2 \qquad \qquad(88)$

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

$d\theta = \frac{{\partial \theta }}{{\partial R}}dR + \frac{{\partial \theta }}{{\partial \tau }}d\tau = 0,{\rm{ }}R = S$

which can be rearranged to get

$\frac{{\partial \theta }}{{\partial R}}\frac{{dS}}{{d\tau }} + \frac{{\partial \theta }}{{\partial \tau }} = 0 \qquad \qquad(89)$

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

$- Ste\mathop {\left( {\frac{{\partial \theta }}{{\partial R}}} \right)}^2 + \frac{1}{R}\frac{\partial }{{\partial R}}\left( {R\frac{{\partial \theta }}{{\partial R}}} \right) = 0\quad R = S\quad \tau > 0 \qquad \qquad(90)$

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

$2 + \varphi = \frac{{\sqrt {1 + 2Ste} - 1}}{{Ste}} \qquad \qquad(91)$

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

$\frac{{dS}}{{d\tau }} = \frac{{\sqrt {1 + 2Ste} - 1}}{{SlnS}} \qquad \qquad(92)$

which is subjected to the following initial condition:

$S(\tau ) = 1,\quad \quad \tau = 0 \qquad \qquad(93)$

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

$\frac{1}{2}{S^2}lnS - \frac{1}{4}({S^2} - 1) = \left( {\sqrt {1 + 2Ste} - 1} \right)\tau \qquad \qquad(94)$

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:

$k\frac{{\partial T}}{{\partial r}} = {h_i}(T - {T_i})\quad r = {r_i}\quad t > 0 \qquad \qquad(95)$

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

$\frac{{\partial \theta }}{{\partial R}} = Bi(\theta - 1)\begin{array}{*{20}{c}} , & {\begin{array}{*{20}{c}} {R = 1} & {\tau = 0} \\ \end{array}} \\ \end{array} \qquad \qquad(96)$

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

$Bi = \frac{{h{r_i}}}{k} \qquad \qquad(97)$

and $θ$ is defined as

$\theta = \frac{{{T_m} - T}}{{{T_m} - {T_i}}} \qquad \qquad(98)$

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

$\theta = A + B\left( {\frac{{lnR}}{{lnS}}} \right) - (A + B)\mathop {\left( {\frac{{lnR}}{{lnS}}} \right)}^2 \qquad \qquad(99)$

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

$B = - Bi(1 - A)lnS \qquad \qquad(100)$

$2A + B = \frac{{\sqrt {1 + 2ASte} - 1}}{{Ste}} \qquad \qquad(101)$

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

$2A - Bi(1 - A)lnS = \frac{{\sqrt {1 + 2ASte} - 1}}{{Ste}} \qquad \qquad(102)$

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

$\frac{{dS}}{{d\tau }} = \frac{{\sqrt {1 + 2ASte} - 1}}{{SlnS}} \qquad \qquad(103)$

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

$\theta (1,\tau ) = A \qquad \qquad(104)$

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.

External Links

Retrieved from "http://www.thermalfluidscentral.org/encyclopedia/index.php/Integral_approximate_solution_of_melting_and_solidification_problems"

Integral approximate solution of melting and solidification problems

From Thermal-FluidsPedia

Contents

One-Region Problem

Boundary condition of the 1st kind

Boundary Condition of the 2nd kind

Two-Region Problem

Duration of Preheating

Governing Equations for the Melting Stage

Integral Approximate Solution

Ablation under Constant Heat Flux Heating

Solidification/Melting in Cylindrical Coordinate Systems

Boundary condition of the 1st kind

Boundary condition of the 3rd kind

References

Further Reading

External Links

Views

Personal tools

Navigation

Search