Exact solutions of melting and solidification problems

From Thermal-FluidsPedia

1 Governing Equations of the Solidification Problem
2 Dimensionless Form of the Governing Equations
3 Exact Solution of the One-Region Problem
4 Exact Solution of the Two-Region Problem
5 References
6 Further Reading
7 External Links

Governing Equations of the Solidification Problem

The physical model of the solidification problem to be investigated in this subsection is shown in Fig. 1(b) from boundary conditions at solid-liquid interface, where a liquid PCM with a uniform initial temperature $T i,$ which exceeds the melting point $T m,$ is in a half-space $x > 0.$ At time $t = 0,$ the temperature at the boundary $x = 0$ is suddenly decreased to a temperature $T 0,$ which is below the melting point of the PCM. Solidification occurs from the time $t = 0.$ This is a two-region solidification problem as the temperatures of both the liquid and solid phases are unknown and must be determined. It is assumed that the densities of the PCM for both phases are the same. Natural convection in the liquid phase is neglected, and therefore the heat transfer mechanism in both phases is pure conduction. The temperature in the solid phase must satisfy

$\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 > 0 \qquad \qquad(1)$ ${T_1}(x,t) = {T_0}\quad \quad x = 0,\quad t > 0 \qquad \qquad(2)$

For the liquid phase, the governing equations 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 > 0 \qquad \qquad(3)$

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

${T_2}(x,t) = {T_i}\quad \quad x > 0,\quad t = 0 \qquad \qquad(5)$

The boundary conditions at the interface are

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

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

Before obtaining the solution of the above problem, a scale analysis of the energy balance eq. (7) is performed. At the solid-liquid interface, the scales of derivatives of solid and liquid temperature are

$\frac{{\partial {T_1}}}{{\partial x}} \sim \frac{{{T_m} - {T_0}}}{s}$ $\frac{{\partial {T_2}}}{{\partial x}} \sim \frac{{{T_i} - {T_m}}}{s}$

The scale of the interface velocity is

$\frac{{ds}}{{dt}} \sim \frac{s}{t}$

Substituting the above equations into eq. (7), one obtains

${k_1}\frac{{{T_m} - {T_0}}}{s} - {k_2}\frac{{{T_i} - {T_m}}}{s} \sim \rho {h_{s\ell }}\frac{s}{t}$

The scale of the location of the solid-liquid interface is obtained by rearranging the above equation, i.e.,

$\frac{{{s^2}}}{{{\alpha _1}t}} \sim {\rm{Ste}}\left( {1 - \frac{{{k_2}}}{{{k_1}}}\frac{{{T_i} - {T_m}}}{{{T_m} - {T_0}}}} \right) \qquad \qquad(8)$

where

${\rm{Ste}} = \frac{{{c_{p1}}({T_m} - {T_0})}}{{{h_{s\ell }}}} \qquad \qquad(9)$

is the Stefan number. Named after J. Stefan, a pioneer in discovery of the solid-liquid phase change phenomena, the Stefan number is a very important dimensionless variable in solid-liquid phase change phenomena, The Stefan number represents the ratio of sensible heat, ${c_{{p_1}}}({T_m} - {T_0}),$ to latent heat, ${h_{s\ell }}$ . For a latent heat thermal energy storage system, the Stefan number is usually very small because the temperature difference in such a system is small, while the latent heat ${h_{s\ell }}$ is very high. Therefore, the effect of the sensible heat transfer on the motion of the solid-liquid interface is very weak, and various approximate solutions to the phase change problem can be introduced without incurring significant error. It can be seen from eq. (8) that the effect of heat conduction in the liquid phase can be neglected if ${T_i} - {T_m} \ll {T_m} - {T_0}$ or ${k_2} \ll {k_1}$ . In that case, eq. (8) can be simplified as

$\frac{{{s^2}}}{{{\alpha _1}t}} \sim {\rm{Ste}} \qquad \qquad(10)$

which means that the interfacial velocity increases with increasing $\Delta T = \left| {{T_m} - {T_0}} \right|$ or decreasing ${h_{s\ell }}$ .

Dimensionless Form of the Governing Equations

The governing eqs. (1) – (7) can be nondimensionalized by introducing the following dimensionless variables:

$\left. \begin{array}{l} \theta = \frac{{{T_m} - T}}{{{T_m} - {T_0}}}\,\,\,\,{\theta _i} = \frac{{{T_m} - {T_i}}}{{{T_m} - {T_0}}}\,\,\,X = \frac{x}{L}\quad S = \frac{s}{L}\quad \tau = \frac{{{\alpha _1}t}}{{{L^2}}} \\ {N_\alpha } = \frac{{{\alpha _2}}}{{{\alpha _1}}}\,\,\,\,\,\,{N_k} = \frac{{{k_2}}}{{{k_1}}}\,\,\,\,{\rm{Ste}} = \frac{{{c_{{p_1}}}({T_m} - {T_0})}}{{{h_{s\ell }}}} \\ \end{array} \right\} \qquad \qquad(11)$

where $L$ is a characteristic length of the problem and can be determined by the nature of the problem or requirement of the solution procedure. The dimensionless governing equations are as follows:

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

${\theta _1}(X,\tau ) = 1\quad \quad X = 0,\quad \tau > 0 \qquad \qquad(13)$

$\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 > 0 \qquad \qquad(14)$

Figure 1: Dimensionless temperature distribution in the PCM

${\theta _2}(X,\tau ) \to {\theta _i}\quad \quad X \to \infty ,\quad \tau > 0 \qquad \qquad(15)$

${\theta _2}(X,\tau ) = {\theta _i}\quad \quad X > 0,\quad \tau = 0 \qquad \qquad(16)$

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

$- \frac{{\partial {\theta _1}}}{{\partial X}} + {N_k}\frac{{\partial {\theta _2}}}{{\partial X}} = \frac{1}{{Ste}}\frac{{dS}}{{d\tau }}\quad \quad X = S(\tau ),\quad \tau > 0 \qquad \qquad(18)$

Dimensionless temperature distribution in a PCM can be qualitively illustrated by Fig. 1. It can be seen that the dimensionless temperature distribution is similar to that of a melting problem. Equations (12) – (18) are also valid for melting problems, provided that the subscripts “1” and “2” represent liquid and solid, respectively. The following solutions of one-region and two-region problems will be valid for both melting and solidification problems.

Exact Solution of the One-Region Problem

Boundary condition of 1st kind

If the initial temperature of a PCM equals its melting point of the PCM ( $T i = T m)$ , the melting or solidification problem is a one-region problem referred to as a Stefan problem. The temperature distribution in one phase is unknown while the temperature in the other phase remains at the melting point. Therefore, the temperature in only one phase needs to be solved. The mathematical description of a Stefan problem can be obtained by simplifying eqs. (12) – (18), i.e.,

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

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

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

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

where the subscript “1” for the liquid phase has been dropped for ease of notation. The temperature distribution of this problem can be constructed on the basis of the heat conduction solution in a semi-infinite body (Ozisik, 1993), i.e.,

$\theta (X,\tau ) = 1 + B{\rm{erf}}(X/2{\tau ^{1/2}}) \qquad \qquad(23)$

where $B$ is an unspecified constant.

It can be verified that eq. (23) satisfies eqs. (19) and (20). The constant $B$ in eq. (23) can be determined by substituting eq. (23) into eq. (21), i.e.,

0 = 1 + B erf(S / 2τ 1 / 2)

Since $B$ is a constant, $S / 2τ 1 / 2$ must also be a constant in order for the above equation to be satisfied. This constant can be represented by $λ,$ so

$\lambda = S/2{\tau ^{1/2}} \qquad \qquad(24)$

Thus,

$B = - \frac{1}{{{\rm{erf}}(\lambda )}}$

By substituting the above expression of $B$ into eq. (23), the dimensionless temperature distribution in phase 1 is obtained:

$\theta (X,\tau ) = 1 - \frac{{{\rm{erf}}(X/2{\tau ^{1/2}})}}{{{\rm{erf}}(\lambda )}} \qquad \qquad(25)$

Substituting eq. (25) and eq. (24) into eq. (22), the following equation is obtained for the constant $λ$ :

$\lambda {e^{{\lambda ^2}}}{\rm{erf}}(\lambda ) = \frac{{{\rm{Ste}}}}{{\sqrt \pi }} \qquad \qquad(26)$

The constant $λ$ is a function of the Stefan number only; this constitutes the appeal of solving the problem by using the dimensionless form of governing equations. Once the constant $λ$ is obtained, the temperature distribution $θ(X,τ)$ and the location of the solid-liquid interface $S (τ)$ can be obtained by eqs. (25) and (24) respectively. The exact solution of the one-region phase change problem can also be obtained by using a similarity solution.

Boundary condition of 2nd kind

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

A solid PCM with a uniform initial temperature equal to the 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. Assume that the densities of the PCM for both phases are the same, and that natural convection in the liquid phase is negligible. We can find the transient location of the solid-liquid interface.

The melting problem under consideration is shown in Fig. 2. Since the temperature of the solid phase is equal to the melting point and only the temperature in the liquid phase is unknown, this is a one-region melting problem. The governing equation, and the corresponding initial and boundary conditions of the problem are as follows

$\frac{{{\partial ^2}T}}{{\partial {x^2}}} = \frac{1}{\alpha }\frac{{\partial T}}{{\partial t}}\quad \quad 0 < x < s(t),\quad t > 0 \qquad \qquad(27)$ $- k\frac{{\partial T(x,t)}}{{\partial x}} = {q''_0}\quad \quad x = 0,\quad t > 0 \qquad \qquad(28)$

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

$- k\frac{{\partial T}}{{\partial x}} = \rho {h_{s\ell }}\frac{{ds}}{{dt}}\quad \quad x = s(t),\quad t > 0 \qquad \qquad(30)$

where the subscript 1 for liquid phase has been dropped to simplify the notation. Introducing the following dimensionless variables:

$\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} \qquad \qquad(31)$

where $L$ is a characteristic length of the problem, eqs. (27) – (30) are nondimensionalized as

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

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

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

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

It is assumed that the temperature distribution in the liquid phase region can be constructed in a fashion similar to the exact solution of conduction in a semi-infinite body with boundary conditions of the second kind:

$\theta (X,\tau ) = A + 2\sqrt \tau {\rm{ierfc}}\left( {\frac{X}{{2\sqrt \tau }}} \right) \qquad \qquad(36)$

where $i e r f c$ is a differential complementary error function:

${\rm{ierfc}}(Z) = \frac{1}{{\sqrt \pi }}{e^{ - {z^2}}} - Z{\rm{erfc}}(Z) \qquad \qquad(37)$

The unknown variable $A$ in eq. (36) can be obtained by substituting eq. (36) into eq. (34), i.e.,

$A = - 2\sqrt \tau {\rm{ierfc}}\left( {\frac{S}{{2\sqrt \tau }}} \right) \qquad \qquad(38)$

Thus, the temperature distribution in the liquid phase is

$\theta (X,\tau ) = 2\sqrt \tau \left[ {{\rm{ierfc}}\left( {\frac{X}{{2\sqrt \tau }}} \right) - {\rm{ierfc}}\left( {\frac{S}{{2\sqrt \tau }}} \right)} \right] \qquad \qquad(39)$

Substituting eq. (39) into eq. (35) yields an ordinary differential equation about the location of the solid-liquid interface, i.e.,

$\frac{{dS}}{{d\tau }} = {\rm{Ste}} \cdot {\rm{erfc}}\left( {\frac{S}{{2\sqrt \tau }}} \right) \qquad \qquad(40)$

which is subject to the following initial condition:

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

A closed-form expression of $S (τ)$ cannot be obtained from eq. (40) because the separation of variables in eq. (40) is impossible. Equation (40) can be solved numerically to determine the transient location of the solid-liquid interface.

The above exact solution in dimensional form was first obtained by El-Genk and Cronenberg (1979). Cho and Sunderland (1981) pointed out that there is a contradiction in this solution, since eq. (36) does not satisfy the partial differential equation (32) unless $A$ in eq. (36) is a constant. However, $A$ in eq. (36) is not a constant because it is a function of dimensionless time $τ$ [see eq. (38)]. The interested reader can find the detailed evaluation of El-Genk and Cronenberg’s (1979) solution in Cho and Sunderland (1981). Nevertheless, this is a first attempt to obtain the exact solution of melting/solidification in a semi-infinite body with boundary conditions other than those of the first kind. In addition to the melting/solidification under boundary conditions of the first and second kinds discussed above, the boundary condition of the third kind (convective heating or cooling at surface) is a more generalized boundary condition for many applications, such as freezing of biological materials. Melting and solidification under the boundary condition of the third kind can be solved by a procedure similar to the above.

Freezing of biological tissue

Figure 3: Freezing of a slab under convective cooling.

A slab of biological material with a thickness of 2L is cooled by a fluid with temperature of ${T_\infty }$ , and the convective heat transfer coefficient is $h$ . It is assumed that the biological material can be treated as a single-component PCM with a well defined melting point, $T m$ , and its initial temperature is at $T m$ . The frozen process is so slow that heat transfer in a frozen layer can be regarded as a pseudo-steady-state process. The time it takes to freeze the entire slab can be estimated by a procedure below.

Since both sides of the slab are cooled, only half of the problem needs to be considered. The latent heat released by the biological material due to freezing is

$q'' = \rho {h_{s\ell }}\frac{{ds}}{{dt}}$

Assuming steady-state conduction in the frozen layer, the latent heat released by freezing must be conducted through the frozen layer, i.e.,

$q'' = ({T_m} - {T_\infty }){\left( {\frac{s}{k} + \frac{1}{h}} \right)^{ - 1}}$

Combining the above two equations, we have

$\frac{{{T_m} - {T_\infty }}}{{\rho {h_{s\ell }}}}dt = \left( {\frac{s}{k} + \frac{1}{h}} \right)ds \qquad \qquad(42)$

The time it takes to freeze the entire slab, $t f$ , can be obtained by integrating eq. (42)

$\int_0^{{t_f}} {\frac{{{T_m} - {T_\infty }}}{{\rho {h_{s\ell }}}}dt} = \int_0^L {\left( {\frac{s}{k} + \frac{1}{h}} \right)ds}$

i.e.,

${t_f} = \frac{{\rho {h_{s\ell }}}}{{{T_m} - {T_\infty }}}\left( {\frac{{{L^2}}}{{2k}} + \frac{L}{h}} \right) \qquad \qquad(43)$

which can be nondimensionalized as

${\tau _f} = \frac{{\rm{1}}}{{{\rm{Ste}}}}\left( {\frac{1}{2} + \frac{{\rm{1}}}{{{\rm{Bi}}}}} \right) \qquad \qquad(44)$

where $τ f = α t f / L 2$ , ${\rm{Ste}} = {c_p}({T_m} - {T_\infty })/{h_{s\ell }}$ , $Bi = h L / k$ are dimensionless time, Stefan number and Biot number, respectively. When ${\rm{Bi}} \to \infty$ , the problem becomes the boundary condition of the first kind [ $T(0,t) = {T_\infty }$ ] and eq. (44) becomes

${\tau _f} = \frac{{\rm{1}}}{{{\rm{2Ste}}}} \qquad \qquad(45)$

Exact Solution of the Two-Region Problem

If the initial temperature of the PCM is not equal to its melting point of the PCM ( ${T_i} \ne {T_m})$ , the melting or solidification problem becomes a two-region problem, referred to as a Neumann problem in the literature. Temperature distributions in both the liquid and solid phases are unknown and must be solved. Equations (12) – (18) provide the complete mathematical description of a Neumann problem. Based on the heat conduction solution of a semi-infinite body, the temperature distribution in the PCM can be constructed as follows:

${\theta _1}(X,\tau ) = 1 + A{\rm{erf(}}X/2{\tau ^{1/2}}) \qquad \qquad(46)$

${\theta _2}(X,\tau ) = {\theta _i} + B{\rm{erfc}}[X/2({N_\alpha }{\tau ^{1/2}})] \qquad \qquad(47)$

where $A$ and $B$ in eqs. (46) and (47) are unspecified constants, and erfc is the complementary error function, defined as

${\rm{erfc}}(z) = 1 - {\rm{erf}}(z) \qquad \qquad(48)$

It should be noted that eqs. (46) and (47) satisfy eqs. (12) – (16). The constants $A$ and $B$ can be determined by using boundary condition (17), i.e.,

1 + A erf(S / 2τ 1 / 2) = 0

θ i + B erfc[S / 2(N α τ 1 / 2)] = 0

Since $A$ and $B$ are constants, $S / 2τ 1 / 2$ must also be a constant in order for the above two equations to be satisfied. This constant can be represented by $λ,$ so

$\lambda = S/2{\tau ^{1/2}} \qquad \qquad(49)$

Thus, the constants $A$ and $B$ can be determined as

$A = - \frac{1}{{{\rm{erf}}(\lambda )}}$ $B = - \frac{{{\theta _i}}}{{{\rm{erfc}}(\lambda /N_\alpha ^{1/2})}}$

The temperature distributions in both phases can be obtained by substituting the above equations for $A$ and $B$ into eqs. (46) and (47), i.e.,

${\theta _1}(X,\tau ) = 1 - \frac{{{\rm{erf}}[X/2{\tau ^{1/2}}]}}{{{\rm{erf}}(\lambda )}} \qquad \qquad(50)$

${\theta _2}(X,\tau ) = {\theta _i}\left[ {1 - \frac{{{\rm{erf}}[X/2{{({N_\alpha }\tau )}^{1/2}}]}}{{{\rm{erf}}(\lambda /N_\alpha ^{1/2})}}} \right] \qquad \qquad(51)$

Substituting eqs. (50), (51), and (49) into eq. (18), the following equation is obtained for the constant $λ$ :

$\frac{{{e^{ - {\lambda ^2}}}}}{{{\rm{erf}}(\lambda )}} + \frac{{{N_k}{\theta _i}}}{{N_\alpha ^{1/2}}}\frac{{{e^{ - \lambda /{N_\alpha }}}}}{{{\rm{erfc}}(\lambda /N_\alpha ^{1/2})}} = \frac{{\lambda \sqrt \pi }}{{{\rm{Ste}}}} \qquad \qquad(52)$

Equation (52) can be solved by using an iterative method with under-relaxation because its left-hand side is a very complicated function of $λ$ . Once $λ$ is obtained, the temperature distributions – $θ 1 (X,τ)$ and $θ 2 (X,τ)$ – and the location of the solid-liquid interface $S (τ)$ can be obtained from eqs. (50), (51), and (49), respectively.

References

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.

Ozisik, M.N., 1993, Heat Conduction, 2^nd ed., Wiley-Interscience, New York.

External Links

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

Exact solutions of melting and solidification problems

From Thermal-FluidsPedia

Contents

Governing Equations of the Solidification Problem

Dimensionless Form of the Governing Equations

Exact Solution of the One-Region Problem

Boundary condition of 1st kind

Boundary condition of 2nd kind

Freezing of biological tissue

Exact Solution of the Two-Region Problem

References

Further Reading

External Links

Views

Personal tools

Navigation

Search