# Integral approximate method for solidification of binary solution

Solidification of a binary solution

The mushy zone formation in binary solutions on a horizontal cold surface will be studied using an integral approximate method. The physical model of the solidification problem is shown in Fig. 1 (Zhang and Faghri, 1997). An ammonium chloride water solution with initial temperature Ti and mass fraction ωi fills a half-space x > 0. At time t=0, the wall temperature at x=0 is suddenly reduced to a temperature Tw, which is higher than the eutectic temperature of the ammonium chloride water solution. The solidification process starts from the cold wall and the mushy zone grows upward. There is no solid phase because the cold wall temperature, Tw, is above the eutectic temperature. Therefore, this is a two-region problem with temperatures in the mushy and liquid zones, as well as the location of interface between these two zones, as unknowns. During the solidification process, a denser and colder solution appears in the mushy zone because the solution near the ice surface rejects the solute (Braga and Viskanta, 1990). Thus, the liquid phase is hydrodynamically stable and no natural convection occurs in the liquid phase. Furthermore, the following assumptions are made in order to simplify the analysis:

1. The growth of the mushy zone is controlled by heat conduction because the thermal diffusivity of the salt solution is 100 times greater than the mass diffusivity (Fang et al., 1984).

2. The properties of the solid and liquid are constant within the liquid and solid, but different from each other, and the properties of the mushy zone are weighted according to the local solid fraction.

3. The densities of the solid and liquid phases are the same, i.e., the density change during the solidification process is neglected.

Based on the above assumptions, the governing energy equations of the solidification problem can be given as follows:

Mushy zone:

${{\rho }_{mu}}\frac{\partial }{\partial t}({{c}_{mu}}{{T}_{mu}})=\frac{\partial }{\partial x}\left( {{k}_{mu}}\frac{\partial {{T}_{mu}}}{\partial x} \right)+{{\rho }_{mu}}{{h}_{s\ell }}\frac{\partial f}{\partial t}\quad \quad 0\le x\le s\qquad\qquad(1)$
$T={{T}_{w}}\quad \quad x=0\qquad\qquad(2)$

Liquid zone:

${{\rho }_{\ell }}{{c}_{p\ell }}\frac{\partial {{T}_{\ell }}}{\partial t}={{k}_{\ell }}\frac{{{\partial }^{2}}{{T}_{\ell }}}{\partial {{x}^{2}}}\quad \quad x\ge s\qquad\qquad(3)$
${{T}_{\ell }}={{T}_{i}}\quad \quad x\to \infty \qquad\qquad(4)$
${{T}_{\ell }}={{T}_{i}}\quad \quad t=0\qquad\qquad(5)$

Mushy zone and liquid zone interface:

${{T}_{mu}}={{T}_{\ell }}={{T}_{s}}\quad \quad x=s\qquad\qquad(6)$
${{k}_{mu}}\frac{\partial {{T}_{mu}}}{\partial x}={{k}_{\ell }}\frac{\partial {{T}_{\ell }}}{\partial x}\quad \quad x=s\qquad\qquad(7)$
$s=0\quad \quad t=0\qquad\qquad(8)$

where the last term in eq. (1) represents the latent heat released during solidification. Since phase change occurs within the entire mushy zone, the latent heat appears as a source term in the energy equation (1), instead of appearing in a boundary condition at interface, as was the case for solid-liquid phase change of single-component PCMs. f is the local solid mass fraction in the mushy zone, which is the same as the volume fraction of solid if the density is not changed during the phase change process. It can be determined from the phase diagram by

$f=\frac{\omega ({{T}_{mu}})-\omega }{\omega ({{T}_{mu}})}\qquad\qquad(9)$

where ω(Tmu) is the mass fraction obtained by liquidus 1 in the phase diagram, Fig. 2.11. Ts in eq. (6) is a liquidus line temperature distinguishing the mushy zone and the liquid zone. The liquidus line temperature can be determined from the liquidus line equation of the phase diagram of ammonium chloride water solution based on the initial concentration. The properties in the mushy zone are mass weighted according to the local solid fraction as follows:

${{\rho }_{mu}}={{\rho }_{\ell }}={{\rho }_{s}}\qquad\qquad(10)$
${{c}_{pmu}}=({{c}_{ps}}-{{c}_{p\ell }})f+{{c}_{p\ell }}\qquad\qquad(11)$
${{k}_{mu}}=({{k}_{s}}-{{k}_{\ell }})f+{{k}_{\ell }}\qquad\qquad(12)$

By defining the following dimensionless variables,

\left. \begin{align} & {{\theta }_{mu}}=\frac{{{T}_{mu}}-{{T}_{s}}}{{{T}_{s}}-{{T}_{w}}}\,\,\,\,\,\,\,\,\begin{matrix} {{\theta }_{\ell }}=\frac{{{T}_{\ell }}-{{T}_{s}}}{{{T}_{s}}-{{T}_{w}}} & {{\theta }_{i}}=\frac{{{T}_{i}}-{{T}_{s}}}{{{T}_{s}}-{{T}_{w}}}\,\,\,\,\,X=\frac{x}{L} \\ \end{matrix}\,\,\, \\ & S=\frac{s}{L}\,\,\,\,\tau =\frac{{{\alpha }_{\ell }}t}{{{L}^{2}}}\,\,\,\,\text{Ste}=\frac{{{c}_{p\ell }}({{T}_{s}}-{{T}_{w}})}{{{h}_{s\ell }}}\,\,\,\,\,{{R}_{c}}=\frac{{{c}_{ps}}}{{{c}_{p\ell }}}\,\,\,\,\,{{R}_{k}}=\frac{{{k}_{s}}}{{{k}_{\ell }}} \\ \end{align} \right\}\qquad\qquad(13)

the governing equations and boundary conditions become

\begin{align} & \frac{\partial }{\partial \tau }\left[ (({{R}_{c}}-1)f+1){{\theta }_{mu}} \right] \\ & =\frac{\partial }{\partial X}\left[ (({{R}_{k}}-1)f+1)\frac{\partial {{\theta }_{mu}}}{\partial X} \right]+\frac{1}{Ste}\frac{\partial f}{\partial \tau }\begin{matrix} {} & 0\le X\le S \\ \end{matrix} \\ \end{align}\qquad\qquad(14)
${{\theta }_{mu}}=-1\quad \quad X=0\qquad\qquad(15)$
$\begin{matrix} \frac{\partial {{\theta }_{\ell }}}{\partial \tau }=\frac{{{\partial }^{2}}{{\theta }_{\ell }}}{\partial {{X}^{2}}} & X>S \\ \end{matrix}\qquad\qquad(16)$
${{\theta }_{\ell }}={{\theta }_{i}}\quad \quad X\to \infty \qquad\qquad(17)$
${{\theta }_{\ell }}={{\theta }_{i}}\quad \quad \tau =0\qquad\qquad(18)$
${{\theta }_{mu}}={{\theta }_{\ell }}=0\quad \quad X=S\qquad\qquad(19)$
$\frac{\partial {{\theta }_{mu}}}{\partial X}=\frac{\partial {{\theta }_{\ell }}}{\partial X}\quad \quad X=S\qquad\qquad(20)$
$\begin{matrix} S=0 & \tau =0 \\ \end{matrix}\qquad\qquad(21)$

The exact analytical solution of eq. (16) can be obtained by using the exact solution of heat conduction in a semi-infinite body, i.e.,

${{\theta }_{\ell }}={{\theta }_{i}}\left[ 1-\frac{\text{erfc}(X/\sqrt{4\tau })}{\text{erfc}(S/\sqrt{4\tau })} \right]\qquad\qquad(22)$

which exactly satisfies eqs. (16) – (19).

For solidification on a cold isothermal surface, the thickness of the mushy zone, S, can be expressed as (Braga and Viskanta, 1990; Ozisik, 1993; Tien and Geiger, 1967)

$\text{S=2}\lambda \sqrt{\tau }\text{ }\qquad\qquad(23)$

where λ is a constant. Integrating eq. (14) over the interval of (0,S) yields the integral equation of the mushy zone:

\begin{align} & \frac{d}{d\tau }\int_{0}^{S}{\left[ (({{R}_{c}}-1)f+1){{\theta }_{mu}} \right]}dX \\ & ={{\left. \frac{\partial {{\theta }_{mu}}}{\partial X} \right|}_{X=S}}-[({{R}_{k}}-1)f+1]{{\left. \frac{\partial {{\theta }_{mu}}}{\partial X} \right|}_{X=0}}+\frac{1}{Ste}\frac{d}{d\tau }\int_{0}^{S}{fdX} \\ \end{align}\qquad\qquad(24)

where ${{\left. \left( \partial {{\theta }_{mu}}/\partial X \right) \right|}_{X=S}}$ can be found from eq. (21) via eq. (17).

Assuming that the temperature distribution in the mushy zone is a second-order polynomial function, and determining the coefficients, one obtains the temperature distribution in the mushy zone:

${{\theta }_{mu}}=\left[ \frac{2\lambda {{e}^{-{{\lambda }^{2}}}}{{\theta }_{i}}}{\text{erfc}(\lambda )} \right]\left( \frac{X-S}{S} \right)+\left[ \frac{2\lambda {{e}^{-{{\lambda }^{2}}}}{{\theta }_{i}}}{\text{erfc}(\lambda )}-1 \right]{{\left( \frac{X-S}{S} \right)}^{2}}\qquad\qquad(25)$

To solve the solidification problem using the integral approximate method, it is necessary to know the distribution of the solid mass fraction f in the mushy zone. The expression of the solid fraction as a single value of temperature in the mushy zone – as in Braga and Viskanta (1990) – is impossible to use here because the integral term in eq. (25) will be very difficult to obtain. In their integral solution of the solidification in a semi-infinite region without liquid superheat, Tien and Geiger (1967) assume a linear distribution of the solid fraction in the mushy zone. In that instance, the solid fraction distribution is found to have no significant effect. Cao and Poulikakos (1991) assumed that the solid fraction, and its derivative with respect to temperature, are constants when solving freezing problems of a binary alloy saturating a packed bed of spheres. For the sake of simplicity, it is also assumed here that the solid fraction is a linear function in the mushy zone, i.e.,

$f={{f}_{w}}\left( 1-\frac{X}{S} \right)\qquad\qquad(26)$

where fw in eq. (26) is the solid mass fraction at the cold isothermal surface, which depends on the cold wall temperature Tw and the initial concentration of the ammonium chloride water solution ωi. It can be determined from the phase diagram by the lever rule (see Example 2.4):

${{f}_{w}}=1-\frac{{{\omega }_{i}}}{\omega ({{T}_{w}})}\qquad\qquad(27)$

where ω(Tw) can be determined by the liquidus line equation:

$\omega ({{T}_{w}})=1.678\times {{10}^{-3}}-1.602\times {{10}^{-2}}{{T}_{w}}-2.857\times {{10}^{-4}}T_{w}^{2}-4.491\times {{10}^{-6}}T_{w}^{3}\qquad\qquad(28)$

Substituting eqs. (25) and (26) into eq. (24) and considering eq. (23), an algebraic equation of λ is obtained:

$\lambda ={{\left\{ \frac{2[({{R}_{k}}-1){{f}_{w}}+1]-\frac{2\lambda {{e}^{-{{\lambda }^{2}}}}{{\theta }_{0}}}{\text{erfc}(\lambda )}[({{R}_{k}}-1){{f}_{w}}+2]}{\frac{2\lambda {{e}^{-{{\lambda }^{2}}}}{{\theta }_{0}}}{\text{erfc}(\lambda )}\frac{{{f}_{w}}+2}{12}+\frac{3({{R}_{c}}-1){{f}_{w}}+4}{12}+\frac{{{f}_{w}}}{2Ste}} \right\}}^{\frac{1}{2}}}\qquad\qquad(29)$

which is simply an algebraic equation that can be solved iteratively. After the value of λ is obtained, the dimensionless thickness of the mushy zone can be obtained from eq. (23), and the temperature distribution in the liquid and mushy zones can be obtained by eqs. (22) and (25). In order to compare these results based on the integral method with Braga and Viskanta’s (1990) experimental results, the following dimensionless variable should be introduced:

$\eta =\frac{X}{2\sqrt{\tau }}\qquad\qquad(30)$

The temperature distributions in the liquid and mushy zones then become

${{\theta }_{m}}=\left[ \frac{2\lambda {{e}^{-{{\lambda }^{2}}}}{{\theta }_{i}}}{erfc(\lambda )} \right]\left( \frac{\eta -\lambda }{\lambda } \right)+\left[ \frac{2\lambda {{e}^{-{{\lambda }^{2}}}}{{\theta }_{i}}}{erfc(\lambda )}-1 \right]{{\left( \frac{\eta -\lambda }{\lambda } \right)}^{2}}\qquad\qquad(32)$

The thermal properties of the ammonium chloride water solution at specified concentrations of 5%, 10%, and 15% can be found in Zeng and Faghri (1994a) and Cao and Poulikakos (1991). In Braga and Viskanta’s (1990) experimental investigation, the solidification of NH4Cl-H2O was performed for three different cases. Zhang and Faghri (1997) compared the mushy zone thickness obtained using the above integral method and by Braga and Viskanta (1990). It was found that the result obtained by the integral approximate solution agreed well with Braga and Viskanta’s similarity solution and the experimental results.

## References

Braga, S. L., and Viskanta, R., 1990, “Solidification of a Binary Solution on a Cold Isothermal Surface,” International Journal of Heat and Mass Transfer, Vol. 33, pp. 745-754.

Cao, W.Z., and Poulikakos, D., 1991, “Freezing of a Binary Alloy Saturating a Packed Bed of Spheres,” AIAA Journal of Thermophysics and Heat Transfer, Vol. 5, pp. 46-53.

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

Fang, L.J., Cheung, F.B., Linehan, J.H., and Pedersen, D.R., 1984, “Selective Freezing of a Dilute Salt Solution on a Cold Ice Surface,” ASME Journal of Heat Transfer, Vol. 106, pp. 385-393.

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

Tien, R.H., Geiger, G.E., 1967, “A Heat-Transfer Analysis of the Solidification of A Binary Eutectic System,” ASME Journal of Heat Transfer, Vol. 91, pp. 230-234.

Zeng, X., and Faghri, A, 1994a, “Temperature-Transforming Model for Binary Solid-Liquid Phase-Change Problems Part I: Physical and Numerical Scheme,” Numerical Heat Transfer, Part B, Vol. 25, pp. 467-480.

Zhang, Y., and Faghri, A., 1997, “Analysis of Freezing in an Eccentric Annulus,” ASME Journal of Solar Energy Engineering, Vol. 119, No. 3, pp. 237-241.