# Multi-fluid model

If spatial averaging is performed for each individual phase within a multiphase control volume, the multi-fluid model is obtained. Additional source terms are needed in these equations to account for the interaction between phases.

## Continuity Equation

The volume average of the continuity equation for the kth phase is obtained by taking extrinsic phase averaging from the local instantaneous continuity equation: $\left\langle \frac{\partial {{\rho }_{k}}}{\partial t} \right\rangle +\left\langle \nabla \cdot {{\rho }_{k}}{{\mathbf{V}}_{k}} \right\rangle =0\qquad \qquad(1)$

where the two terms on the left-hand side can be obtained by using $\left\langle \frac{\partial {{\Omega }_{k}}}{\partial t} \right\rangle =\frac{\partial \left\langle {{\Omega }_{k}} \right\rangle }{\partial t}-\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\Omega }_{k}}{{\mathbf{V}}_{I}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}$

and $\left\langle \nabla \cdot {{\Omega }_{k}} \right\rangle =\nabla \cdot \left\langle {{\Omega }_{k}} \right\rangle +\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\Omega }_{k}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}$

from averaging approaches, i.e., $\left\langle \frac{\partial {{\rho }_{k}}}{\partial t} \right\rangle =\frac{\partial \left\langle {{\rho }_{k}} \right\rangle }{\partial t}-\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\rho }_{k}}{{\mathbf{V}}_{I}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}$ $\left\langle \nabla \cdot {{\rho }_{k}}{{\mathbf{V}}_{k}} \right\rangle =\nabla \cdot \left\langle {{\rho }_{k}}{{\mathbf{V}}_{k}} \right\rangle +\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\rho }_{k}}{{\mathbf{V}}_{k}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}$

Substituting the above expressions into eq. (1), the volume-averaged continuity equation becomes $\frac{\partial \left\langle {{\rho }_{k}} \right\rangle }{\partial t}+\nabla \cdot \left\langle {{\rho }_{k}}{{\mathbf{V}}_{k}} \right\rangle =-\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\rho }_{k}}({{\mathbf{V}}_{k}}-{{\mathbf{V}}_{I}})\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}\qquad \qquad(2)$

The right-hand side of eq. (2) represents mass transfer per unit volume from all other phases to the kth phase due to phase change; it can be rewritten as $-\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\rho }_{k}}({{\mathbf{V}}_{k}}-{{\mathbf{V}}_{I}})\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}=\sum\limits_{j=1(j\ne k)}^{\Pi }{{{{{\dot{m}}'''}}_{jk}}}\qquad \qquad(3)$

where ${{{\dot{m}}'''}_{jk}}$ represents mass transfer per unit volume from the jth to the kth phase due to phase change. The value of ${{{\dot{m}}'''}_{jk}}$ depends on the phase change process that takes place in the multiphase system, and the conservation of mass requires that ${{{\dot{m}}'''}_{jk}}=-{{{\dot{m}}'''}_{kj}}$.

The extrinsic phase-averaged density, $\left\langle {{\rho }_{k}} \right\rangle$, is related to the intrinsic phase-averaged density, ${{\left\langle {{\rho }_{k}} \right\rangle }^{k}}$, by $\left\langle {{\rho }_{k}} \right\rangle ={{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}\qquad \qquad(4)$

Furthermore, the intrinsic phase-averaged density is equal to the density ρk.

Substituting eqs. (3) and (4) into eq. (2), and considering eq. $\left\langle {{\Phi }_{k}}{{\Psi }_{k}} \right\rangle ={{\varepsilon }_{k}}{{\left\langle {{\Phi }_{k}} \right\rangle }^{k}}{{\left\langle {{\Psi }_{k}} \right\rangle }^{k}}+\left\langle {{{\hat{\Phi }}}_{k}}{{{\hat{\Psi }}}_{k}} \right\rangle$ from Averaging approaches, the continuity equation for the kth phase becomes $\frac{\partial }{\partial t}\left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}} \right)+\nabla \cdot \left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}{{\left\langle {{\mathbf{V}}_{k}} \right\rangle }^{k}}+\left\langle {{{\hat{\rho }}}_{k}}{{{\mathbf{\hat{V}}}}_{k}} \right\rangle \right)=\sum\limits_{j=1(j\ne k)}^{\Pi }{{{{{\dot{m}}'''}}_{jk}}}\qquad \qquad(5)$

The dispersive term in eq. (5), $\left\langle {{{\hat{\rho }}}_{k}}{{{\mathbf{\hat{V}}}}_{k}} \right\rangle$, is generally small compared with ${{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}{{\left\langle {{\mathbf{V}}_{k}} \right\rangle }^{k}}$; it is assumed that it can be neglected. The continuity equation for the kth phase becomes $\frac{\partial }{\partial t}\left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}} \right)+\nabla \cdot \left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}{{\left\langle {{\mathbf{V}}_{k}} \right\rangle }^{k}} \right)=\sum\limits_{j=1(j\ne k)}^{\Pi }{{{{{\dot{m}}'''}}_{jk}}}\qquad \qquad(6)$

## Momentum Equation

The extrinsic phase-averaged momentum equation for the kth phase can be obtained by performing extrinsic phase-averaging on the momentum equation: $\left\langle \frac{\partial ({{\rho }_{k}}{{\mathbf{V}}_{k}})}{\partial t} \right\rangle +\left\langle \nabla \cdot ({{\rho }_{k}}{{\mathbf{V}}_{k}}{{\mathbf{V}}_{k}}) \right\rangle =\left\langle \nabla \cdot {{{\mathbf{{\tau }'}}}_{k}} \right\rangle +\left\langle {{\rho }_{k}}{{\mathbf{X}}_{k}} \right\rangle \qquad \qquad(7)$

where the body force per unit mass is assumed to be the same for different species for sake of simplicity. After evaluating each term in eq. (7), the multi-fluid volume-averaged momentum equation becomes(Faghri and Zhang, 2006) \begin{align} & \frac{\partial }{\partial t}\left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}{{\left\langle {{\mathbf{V}}_{k}} \right\rangle }^{k}} \right)+\nabla \cdot \left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}{{\left\langle {{\mathbf{V}}_{k}}{{\mathbf{V}}_{k}} \right\rangle }^{k}} \right) \\ & =\nabla \cdot \left( {{\varepsilon }_{k}}{{\left\langle {{{\mathbf{{\tau }'}}}_{k}} \right\rangle }^{k}} \right)+{{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}{{\mathbf{X}}_{k}}+\sum\limits_{j=1(j\ne k)}^{\Pi }{\left( \left\langle {{\mathbf{F}}_{jk}} \right\rangle +\left\langle {{{{\dot{m}}'''}}_{jk}} \right\rangle {{\left\langle {{\mathbf{V}}_{k,I}} \right\rangle }^{k}} \right)} \\ \end{align}\qquad \qquad(8)

where ${{\left\langle {{\mathbf{V}}_{k,I}} \right\rangle }^{k}}$ is intrinsic phase-averaged velocity of the kth phase at the interface. The difference between two adjacent phases results solely from the density difference between the two phases. $\left\langle {{\mathbf{F}}_{jk}} \right\rangle$ is an interactive force between the jth and the kth phase, and depends on the friction, pressure, and cohesion between different phases. Newton’s third law requires that the interactive forces satisfy $\left\langle {{\mathbf{F}}_{jk}} \right\rangle =-\left\langle {{\mathbf{F}}_{kj}} \right\rangle \qquad \qquad(9)$

The interactive force can be determined by $\left\langle {{\mathbf{F}}_{jk}} \right\rangle ={{K}_{jk}}\left( {{\left\langle {{\mathbf{V}}_{j}} \right\rangle }^{j}}-{{\left\langle {{\mathbf{V}}_{k}} \right\rangle }^{k}} \right)\qquad \qquad(10)$

where Kjk is the momentum exchange coefficient between phases j and k. Determining the momentum exchange coefficient is a very challenging task because interphase momentum exchange depends on the structure of the interfaces. If a secondary phase j is dispersed in the primary phase k, as is the case with the dispersed phase system summarized in Table 1.8, one can assume that the secondary phase is spherical in shape and an appropriate empirical correlation can be used to obtain the momentum exchange coefficient.

Since liquid-vapor flow is widely used in various applications, we will use liquid-vapor flow as an example to explain the determination of the momentum exchange coefficient. If liquid is the primary phase and vapor is the secondary phase, the vapor phase is dispersed in the liquid as vapor bubbles. If vapor is the primary phase and liquid is the secondary phase, the liquid phase is dispersed in the vapor as liquid droplets. Boysan (1990) suggested that the momentum exchange coefficient could be estimated by ${{K}_{jk}}=\frac{3}{4}{{C}_{D}}\frac{{{\varepsilon }_{j}}\left\langle {{\rho }_{k}} \right\rangle }{{{d}_{j}}}\left| {{\left\langle {{\mathbf{V}}_{j}} \right\rangle }^{j}}-{{\left\langle {{\mathbf{V}}_{k}} \right\rangle }^{k}} \right|\qquad \qquad(11)$

where phase k is the primary phase, and phase j is the secondary phase, and dj is the diameter of vapor bubbles or liquid droplets of the secondary phase j. CD is the drag coefficient based on the relative Reynolds number, which obtained by the following empirical correlations: {{C}_{D}}=\left\{ \begin{align} & \frac{24}{\operatorname{Re}}(1+0.15{{\operatorname{Re}}^{0.687}}) \\ & 0.44 \\ \end{align} \right.\begin{matrix} {} & \begin{align} & \operatorname{Re}\le 1000 \\ & \operatorname{Re}>1000 \\ \end{align} \\ \end{matrix}\qquad \qquad(12)

where $\operatorname{Re}=\frac{\left\langle {{\rho }_{k}} \right\rangle \left| {{\left\langle {{\mathbf{V}}_{j}} \right\rangle }^{j}}-{{\left\langle {{\mathbf{V}}_{k}} \right\rangle }^{k}} \right|{{d}_{j}}}{{{\mu }_{k}}}\qquad \qquad(13)$

## Energy Equation

The extrinsic phase-average of the energy equation is \begin{align} & \left\langle \frac{\partial ({{\rho }_{k}}{{h}_{k}})}{\partial t} \right\rangle +\left\langle \nabla \cdot {{\rho }_{k}}{{\mathbf{V}}_{k}}{{h}_{k}} \right\rangle \\ & =-\left\langle \nabla \cdot {{{\mathbf{{q}''}}}_{k}}_{k} \right\rangle +\left\langle {{{{q}'''}}_{k}} \right\rangle +\left\langle \frac{\partial {{p}_{k}}}{\partial t} \right\rangle +\left\langle {{\mathbf{V}}_{k}}\cdot \nabla {{p}_{k}} \right\rangle +\left\langle \nabla {{\mathbf{V}}_{k}}:\mathbf{\tau }{}_{k} \right\rangle \\ \end{align}\qquad \qquad(14)

which can be used to obtain the volume-averaged energy equation and the result is (Faghri and Zhang, 2006) \begin{align} & \frac{\partial }{\partial t}\left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}{{\left\langle {{h}_{k}} \right\rangle }^{k}} \right)+\nabla \cdot \left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k}} \right\rangle }^{k}}{{\left\langle {{\mathbf{V}}_{k}}{{h}_{k}} \right\rangle }^{k}} \right)=-\nabla \cdot \left\langle {{{\mathbf{{q}''}}}_{k}} \right\rangle +\left\langle {{{{q}'''}}_{k}} \right\rangle \\ & +{{\varepsilon }_{k}}\frac{D{{\left\langle {{p}_{k}} \right\rangle }^{k}}}{Dt}+\nabla \left\langle {{\mathbf{V}}_{k}} \right\rangle :\left\langle {{\mathbf{\tau }}_{k}} \right\rangle +\sum\limits_{j=1(j\ne k)}^{\Pi }{\left[ \left\langle {{{{q}'''}}_{jk}} \right\rangle +{{{{\dot{m}}'''}}_{jk}}{{\left\langle {{h}_{k,I}} \right\rangle }^{k}} \right]} \\ \end{align}\qquad \qquad(15)

where $\left\langle {{{{q}'''}}_{jk}} \right\rangle$ is the intensity of heat exchange between phase j and k. It can be obtained by using Newton’s law of cooling: $\left\langle {{{{q}'''}}_{jk}} \right\rangle =\frac{{{h}_{c}}\Delta {{A}_{j}}\left( {{\left\langle {{T}_{j}} \right\rangle }^{j}}-{{\left\langle {{T}_{k}} \right\rangle }^{k}} \right)}{\Delta {{V}_{j}}}\qquad \qquad(16)$

where hc is the convective heat transfer coefficient, ΔAj is the area of the interface between phases j and k, and ΔVj is the volume of the secondary phase in the volume element. Like the momentum exchange coefficient, the interphase heat transfer also depends on the structure of the interfaces. If a secondary phase j is dispersed in the primary phase k – as in the dispersed phase summarized in Table 1.8– the following empirical correlation recommended is widely used: $Nu=2+(0.4{{\operatorname{Re}}^{1/2}}+0.06{{\operatorname{Re}}^{2/3}})\Pr _{k}^{0.4}{{\left( \frac{{{\mu }_{k}}}{{{\mu }_{k,s}}} \right)}^{1/4}}\qquad \qquad(17)$

where the Reynolds number, Re, is obtained by eq. (13), the Nusselt number is defined as $Nu=\frac{{{h}_{c}}{{d}_{j}}}{{{k}_{k}}}\qquad \qquad(18)$

and all thermal properties of the primary phase are evaluated at Tk except μk,s, which is evaluated at Tj. Equation (17) is valid for $3.5<\operatorname{Re}<7.6\times {{10}^{4}}$ and $0.71<{{\Pr }_{k}}<380$, which covers a wide variety of problems.

If the secondary phase is liquid and the primary phase is vapor (gas), eq. (17) can be simplified to $Nu=2+0.6{{\operatorname{Re}}^{1/2}}\Pr _{k}^{1/3}\qquad \qquad(19)$

## Species

If the fluid undergoing phase change involves multiple components, it is also necessary to write the equation for conservation of the species mass in the kth phase. The extrinsic phase-average of conservation of species mass can be obtained by $\left\langle \frac{\partial {{\rho }_{k,i}}}{\partial t} \right\rangle +\left\langle \nabla \cdot {{\rho }_{k,i}}{{\mathbf{V}}_{k}} \right\rangle =-\left\langle \nabla \cdot {{\mathbf{J}}_{k,i}} \right\rangle +\left\langle {{{{\dot{m}}'''}}_{k,i}} \right\rangle \qquad \qquad(20)$

where each term can be obtained by $\left\langle \frac{\partial {{\rho }_{k,i}}}{\partial t} \right\rangle =\frac{\partial \left\langle {{\rho }_{k,i}} \right\rangle }{\partial t}-\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\rho }_{k,i}}{{\mathbf{V}}_{I}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}$ $\left\langle \nabla \cdot {{\rho }_{k,i}}{{\mathbf{V}}_{k}} \right\rangle =\nabla \cdot \left\langle {{\rho }_{k,i}}{{\mathbf{V}}_{k}} \right\rangle +\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\rho }_{k,i}}{{\mathbf{V}}_{k}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}$ $\left\langle \nabla \cdot {{\mathbf{J}}_{k,i}} \right\rangle =\nabla \cdot \left\langle {{\mathbf{J}}_{k,i}} \right\rangle +\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\mathbf{J}}_{k,i}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}$

Substituting the above expression into eq. (20), one obtains \begin{align} & \frac{\partial \left\langle {{\rho }_{k,i}} \right\rangle }{\partial t}+\nabla \cdot \left\langle {{\rho }_{k,i}}{{\mathbf{V}}_{k}} \right\rangle =-\nabla \cdot \left\langle {{\mathbf{J}}_{k,i}} \right\rangle +\left\langle {{{{\dot{m}}'''}}_{k,i}} \right\rangle \\ & -\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\rho }_{k,i}}({{\mathbf{V}}_{k}}-{{\mathbf{V}}_{I}})\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}+\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\mathbf{J}}_{k,i}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}} \\ \end{align}\qquad \qquad(21)

where the third and fourth terms in the right-hand side of eq. (21) represent mass source (or sink) of the ith component in the kth phase due to phase change from other phases to the kth phase, as well as mass transfer at the interface due to diffusion, respectively, i.e., $-\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\rho }_{k,i}}({{\mathbf{V}}_{k}}-{{\mathbf{V}}_{I}})\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}+\frac{1}{\Delta V}\int_{{{A}_{k}}}{{{\mathbf{J}}_{k,i}}\cdot {{\mathbf{n}}_{k}}d{{A}_{k}}}=\sum\limits_{j=1(j\ne k)}^{\Pi }{{{{{\dot{m}}'''}}_{jk,i}}}\qquad \qquad(22)$

where ${{{\dot{m}}'''}_{jk,i}}$ represents the mass source (or sink) of the ith component in phase k due to phase change from phase j to phase k, as well as diffusive mass transfer at the interface between phases j and k.

Substituting eq. (22) into eq. (21) and considering eq. $\left\langle {{\Phi }_{k}}{{\Psi }_{k}} \right\rangle ={{\varepsilon }_{k}}{{\left\langle {{\Phi }_{k}} \right\rangle }^{k}}{{\left\langle {{\Psi }_{k}} \right\rangle }^{k}}+\left\langle {{{\hat{\Phi }}}_{k}}{{{\hat{\Psi }}}_{k}} \right\rangle$ from Averaging approaches, the volume-averaged species equation becomes $\frac{\partial \left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k,i}} \right\rangle }^{k}} \right)}{\partial t}+\nabla \cdot \left( {{\varepsilon }_{k}}{{\left\langle {{\rho }_{k,i}} \right\rangle }^{k}}{{\left\langle {{\mathbf{V}}_{k}} \right\rangle }^{k}} \right)=-\nabla \cdot \left\langle {{\mathbf{J}}_{k,i}} \right\rangle +\left\langle {{{{\dot{m}}'''}}_{k,i}} \right\rangle +\sum\limits_{j=1(j\ne k)}^{\Pi }{{{{{\dot{m}}'''}}_{jk,i}}}\qquad \qquad(23)$

where the three terms on the right-hand side represent the effects of mass diffusion, species source/sink due to chemical reaction, and phase change at the interfaces.

## References

Boysan, F., 1990, A Two-Fluid Model for Fluent, Flow Simulation Consultants, Ltd., Sheffield, England.

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.