# Mixture model for solidification of binary solution

The purpose here is to develop and improve the temperature transforming model so as to deal more effectively with binary solid-liquid phase-change problems, with emphasis on the inclusion of the effect of the mushy zone in the energy equation. The essential feature of the model proposed by Zeng and Faghri (1994a, b) is its separation of the coupled effects of temperature and concentration on the latent heat evolution within the mushy zone.

The assumptions made by Zeng and Faghri (1994a) using the mixture model (mass-averaged) include: (1) the thermophysical properties of two phases are constant, although not necessarily equal, and (2) the solid and liquid in the mushy zone are assumed to be in thermal equilibrium. The governing equations will be given for solidification of a binary solution in a two-dimensional rectangular cavity cooled from the sidewall.

The mixture continuity equation is $\frac{D\tilde{\rho }}{Dt}+\tilde{\rho }\nabla \cdot \mathbf{\tilde{V}}=0\qquad\qquad(1)$

where the mixture density is defined as $\tilde{\rho }={{\varepsilon }_{s}}{{\rho }_{s}}+{{\varepsilon }_{\ell }}{{\rho }_{\ell }}\qquad\qquad(2)$

and the volume fractions of solid and liquid ${{\varepsilon }_{s}}\text{ and }{{\varepsilon }_{\ell }}$ satisfy ${{\varepsilon }_{s}}+{{\varepsilon }_{\ell }}=1\qquad\qquad(3)$

The mass-averaged velocity in eq. (1) is defined as $\mathbf{\tilde{V}}={{f}_{s}}{{\mathbf{V}}_{s}}+{{f}_{\ell }}{{\mathbf{V}}_{\ell }}\qquad\qquad(4)$

where ${{\mathbf{V}}_{s}}$ and ${{\mathbf{V}}_{\ell }}$ are velocities of solid and liquid phases, respectively; the mass fractions of solid and liquid satisfy ${{f}_{s}}+{{f}_{\ell }}=1\qquad\qquad(5)$

The volume fraction $\varepsilon$ and the mass fraction f are related by ${{f}_{s}}=\frac{{{\rho }_{s}}{{\varepsilon }_{s}}}{{{\varepsilon }_{s}}{{\rho }_{s}}+{{\varepsilon }_{\ell }}{{\rho }_{\ell }}}\qquad\qquad(6)$ ${{f}_{\ell }}=\frac{{{\rho }_{\ell }}{{\varepsilon }_{\ell }}}{{{\varepsilon }_{s}}{{\rho }_{s}}+{{\varepsilon }_{\ell }}{{\rho }_{\ell }}}\qquad\qquad(7)$

The momentum equations in x- and y-directions are $\frac{\partial (\tilde{\rho }\tilde{u})}{\partial t}+\nabla \cdot (\tilde{\rho }\mathbf{\tilde{V}}\tilde{u})=-\frac{\partial p}{\partial x}+\nabla \cdot \left( {{\mu }_{\ell }}\frac{{\tilde{\rho }}}{{{\rho }_{\ell }}}\nabla \tilde{u} \right)-\frac{{{\mu }_{\ell }}}{K}\frac{{\tilde{\rho }}}{{{\rho }_{\ell }}}(\tilde{u}-{{u}_{s}})+X\qquad\qquad(8)$ $\frac{\partial (\tilde{\rho }\tilde{v})}{\partial t}+\nabla \cdot (\tilde{\rho }\mathbf{\tilde{V}}\tilde{v})=-\frac{\partial p}{\partial y}+\nabla \cdot \left( {{\mu }_{\ell }}\frac{{\tilde{\rho }}}{{{\rho }_{\ell }}}\nabla \tilde{v} \right)-\frac{{{\mu }_{\ell }}}{K}\frac{{\tilde{\rho }}}{{{\rho }_{\ell }}}(\tilde{v}-{{v}_{s}})+Y\qquad\qquad(9)$

where us and vs are the components of the solid velocity, Vs, in the x- and y-directions, respectively. X and Y are components of body force in the x- and y- directions. The energy equation is $\frac{\partial }{\partial t}(\tilde{\rho }\tilde{h})+\nabla \cdot (\tilde{\rho }\mathbf{\tilde{V}}\tilde{h})=\nabla \cdot (\tilde{k}\nabla T)-\nabla \cdot [\tilde{\rho }({{h}_{\ell }}-\tilde{h})(\mathbf{\tilde{V}}-{{\mathbf{V}}_{s}})]\qquad\qquad(10)$

where the mixture enthalpy is defined as $h={{f}_{s}}{{h}_{s}}+{{f}_{\ell }}{{h}_{\ell }}\qquad\qquad(11)$

The mass species equation is \begin{align} & \frac{\partial }{\partial t}(\tilde{\rho }\omega )+\nabla \cdot (\tilde{\rho }\mathbf{\tilde{V}}\omega )=\nabla \cdot (\tilde{\rho }D\nabla \omega ) \\ & +\nabla \cdot [\tilde{\rho }D\nabla ({{\omega }_{\ell }}-\omega )]-\nabla \cdot [\tilde{\rho }({{\omega }_{\ell }}-\omega )(\mathbf{\tilde{V}}-{{\mathbf{V}}_{s}})] \\ \end{align}\qquad\qquad(12)

where the species mass fraction is defined as $\omega ={{f}_{s}}{{\omega }_{s}}+{{f}_{\ell }}{{\omega }_{\ell }}\qquad\qquad(13)$

The permeability of the two-phase mushy zone, which is modeled as a porous medium, is $K={{K}_{0}}\frac{\varepsilon _{\ell }^{3}}{{{(1-{{\varepsilon }_{\ell }})}^{2}}}\qquad\qquad(14)$

where K0 is the permeability constant. The thermal conductivity of the mixture is $\tilde{k}={{\varepsilon }_{s}}{{k}_{s}}+{{\varepsilon }_{\ell }}{{k}_{\ell }}\qquad\qquad(15)$

Under constant specific heat assumption, the enthalpy and temperature are related by the specific heats as ${{h}_{s}}=\int_{0}^{T}{{{c}_{ps}}dT}={{c}_{ps}}T\qquad\qquad(16)$ ${{h}_{\ell }}=\int_{0}^{T}{{{c}_{p\ell }}}dT+{{h}_{0}}={{c}_{p\ell }}T+{{h}_{0}}\qquad\qquad(17)$

where h0 is the reference enthalpy for the liquid phase. The mass diffusivity, in light of the assumption neglecting diffusion in the solid phase, is $D={{f}_{\ell }}{{D}_{\ell }}\qquad\qquad(18)$

Substituting eqs. (16) and (17) into eq. (11), the enthalpy becomes $\tilde{h}={{\tilde{c}}_{p}}T+{{f}_{\ell }}{{h}_{0}}\qquad\qquad(19)$

where ${{\tilde{c}}_{p}}={{f}_{s}}{{c}_{p}}_{s}+{{f}_{\ell }}{{c}_{p}}_{\ell }\qquad\qquad(20)$

Substituting eq. (19) into eq. (10), the energy equation becomes (Zeng and Faghri, 1994a) $A(T,\omega )\frac{\partial T}{\partial t}+\nabla \cdot (\tilde{\rho }\mathbf{V}{{c}_{p\ell }}T)=\nabla \cdot (\tilde{k}\nabla T)-B(T,\omega )\qquad\qquad(21)$

where A(T,ω) is an effective heat coefficient, $A(T,\omega )=\tilde{\rho }{{\tilde{c}}_{p}}+a(T,\omega )\frac{\partial {{f}_{\ell }}}{\partial T}\qquad\qquad(22)$

and B(T,ω) is a heat source term $B(T,\omega )=-\nabla \cdot (\tilde{\rho }{{h}_{0}}\mathbf{\tilde{V}})+\nabla \cdot \left\{ \tilde{\rho }{{f}_{s}}[{{h}_{0}}+({{c}_{p\ell }}-{{c}_{ps}})T]{{\mathbf{V}}_{s}} \right\}-a(T,\omega )\frac{\partial {{f}_{\ell }}}{\partial \omega }\frac{\partial \omega }{\partial t}\qquad\qquad(23)$

where $a(T,\omega )=\tilde{\rho }\left[ ({{c}_{p\ell }}-{{c}_{ps}})T+\frac{({{\rho }_{\ell }}-{{\rho }_{s}}){{{\tilde{c}}}_{p}}T+{{\rho }_{\ell }}{{h}_{0}}}{{{\rho }_{\ell }}-{{f}_{\ell }}({{\rho }_{\ell }}-{{\rho }_{s}})} \right]\qquad\qquad(24)$

The temperature-transforming model presented by Zeng and Faghri (1994a) can be considered as a combination of the apparent heat capacity method and the latent heat source term method designed to deal with the latent heat evolution in binary phase-change systems. The latent heat evolution with reference to temperature evolution is accounted for in the energy equation by defining the effective heat coefficient A(T,ω), which is, in turn, similar to the apparent heat capacity $\tilde{\rho }d\tilde{h}/dT$, while the latent heat evolution – with reference to concentration variation – is incorporated into the source term B(T,ω), which is similar to the source-based method. The model develops and expands the fixed-grid method to a binary phase-change system and links the source-based and the apparent heat capacity methods. For the case of $\partial {{f}_{\ell }}/\partial \omega =0$, that is, $A(T,\omega )=A(T)\qquad\qquad(25)$ $B(T,\omega )=B(T)\qquad\qquad(26)$

eq. (21) reduces to the apparent heat capacity formulation. For the case of $\partial {{f}_{\ell }}/\partial T=0$, $A(T,\omega )=\tilde{\rho }{{\tilde{c}}_{p}}\qquad\qquad(27)$ $B(T,\omega )=B(T,\omega )\qquad\qquad(28)$

eq. (21) represents the so-called source-based energy equation.

It should be pointed out that, in an iterative solution of the discrete equation, the extra term $a(T,\omega )(\partial {{f}_{\ell }}/\partial T)$ in the effective heat coefficient A(T,ω), which is positive, will play the role of a damping factor, in that the temperature of the nodes in which phase change is occurring will change more slowly from one iteration to the next. This stability may represent a main advantage of the apparent heat capacity method. A disadvantage, however, is that if the liquid fraction curve ${{f}_{\ell }}(T,\omega )$ is very steep or contains some discontinuities, instabilities may occur in the solution procedure. This tendency is reduced in the present formulation, because the binary phase-change system generally releases or absorbs latent heat over a relatively wide temperature range. Furthermore, the energy equation has been transformed into a nonlinear equation with a single dependent variable T. Over the iterative loop that solves the energy equation, the extra term $-a(T,\omega )(\partial {{f}_{\ell }}/\partial \omega )(\partial \omega /\partial T)$ in the source term is seen only as a function of temperature.

Zeng and Faghri (1994b) solved the solidification of an NH4Cl solution in a two-dimensional rectangular cavity of dimensions 36×144 mm2 (width × height). Initially, the cavity is insulated and charged with a superheated solution at temperature Ti and concentration ωi= 30% of NH4Cl by weight. Solidification is initiated at t = 0 by holding the left vertical boundary to a temperature Tc < Te, while maintaining the right vertical boundary at the initial temperature. Since the initial concentration of the solution was greater than the eutectic concentration ωe the phase change occurred within mushy zone 2, as shown in Figure 2.11. It is further assumed that: (1) the densities of the two phases are equal and constant, and (2) the velocity of the solid phase is zero. The components of body force in the x- and y- directions are: $X=0\qquad\qquad(29)$ $Y=\tilde{\rho }g[{{\beta }_{T}}(T-{{T}_{ref}})+{{\beta }_{s}}({{\omega }_{ref}}-\omega )]\qquad\qquad(30)$

where βT and βs are thermal and solutal expansion coefficients, respectively. With the additional assumptions, a(T,ω) and B(T,ω) are reduced to $a(T,\omega )=\tilde{\rho }\left[ ({{c}_{p\ell }}-{{c}_{ps}})T+{{h}_{s\ell }} \right]\qquad\qquad(31)$ $B(T,\omega )=-a(T,\omega )\frac{\partial {{f}_{\ell }}}{\partial \omega }\frac{\partial \omega }{\partial t}\qquad\qquad(32)$

where $\frac{\partial {{f}_{\ell }}}{\partial T}=\frac{{{T}_{m}}-{{T}_{\ell }}}{(1-{{k}_{p}}){{({{T}_{m}}-T)}^{2}}}\qquad\qquad(33)$ $\frac{\partial {{f}_{\ell }}}{\partial \omega }=\frac{{{T}_{m}}-{{T}_{e}}}{(1-{{k}_{p}})({{T}_{m}}-T){{\omega }_{e}}}\qquad\qquad(34)$

where Tm is the fusion temperature for ω = 1, and kp is the equilibrium partition ratio.

Zeng and Faghri (1994b) solved the above governing equations using a finite volume method. A 40×40 grid, in which the space step in the y-direction remained constant, while the variation in the x-direction was expressed as a symmetric arithmetic progression, was used to compute this problem with a time step of 3 s. Iterations performed within a time step were terminated when the residual mass within each control volume fell below 10-5, and when the following values fell below 0.001%: (1) changes in the minimum and maximum concentration within each control volume, (2) the average cold wall heat extraction, and (3) total liquid fraction.

## References

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

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.

Zeng, X., and Faghri, A., 1994b, “Temperature-Transforming Model for Binary Solid-Liquid Phase-Change Problems Part II: Numerical Simulation,” Numerical Heat Transfer, Part B, Vol. 25, pp. 481-500.