# Application of Computational Methods

The computational methods presented in the previous topics provide the solution techniques for solving the governing equations and boundary conditions. However, it is equally important to mathematically set up the physical problem before applying the numerical techniques. In many practical problems for external flow one needs to deal with conjugate effects, free surface, and non-conventional geometry including multiple domain regions with non-uniform boundary conditions. In order to show these effects, and the approach to deal with them, an application is presented below that includes all these non-conventional effects[1].

Controlled liquid impinging jet onto a rotating disk
Computational domain for a controlled liquid impinging jet onto a rotating disk.

The physical problem that will be presented in this section is heat and mass transfer from a controlled liquid impinging jet over a rotating disk, as shown in the first picture, including the conjugate wall effect for both heating with and without evaporation from the free surface. First we need to develop the conservative governing equations and boundary conditions for fluid mechanics, heat transfer, and volume of fluid (VOF) for a controlled liquid impinging jet on a rotating disk with a free surface. The second figure to the right shows a computational domain for the three regions: disk, liquid, and vapor, with an extended domain to include the edge effects. The basic assumptions are two dimensional, constant properties, incompressible, and laminar flow. For evaporation, the effect of mass transfer is neglected. It is also assumed that there is constant heat flux along the heated section.

Flow enters the disk between two circular plates, one being the collar, and the other being the disk. The space between the collar and the disk is δin. After the flow leaves the entrance region between the collar and the disk at rhin, the flow turns from an internal fully developed flow to a free surface flow. The heater provides a constant heat flux at the outside wall, from the bottom of the aluminum disk. The heat is conducted through the disk to the fluid. The basic assumptions of the problem at hand are that the flow field is incompressible, and the fluid is considered to be in the laminar flow regime while all of the fluid properties are constant. When evaporation is being modeled, no mass transfer is modeled and the effects are only considered in the energy equation, because the evaporation rate of the liquid is much less than the inlet mass flow rate of the liquid. The Navier-Stokes equations are solved to compute the fluid flow. The continuity and momentum equations are as follows:

 $\nabla \cdot V=0$ (1)
 $\rho \frac{DV}{Dt}=-\nabla p+\nabla \cdot \left( \mu \nabla V \right)+\rho \mathbf{g}+\mathbf{F}$ (2)

Since the flow field is assumed to be independent of the temperature field, a steady-state energy equation is solved in the fluid region once the flow field has been resolved.

 $\nabla \cdot \left( V\rho h \right)=\nabla \cdot \left( k\nabla T \right)$ (3)

In the solid region, only conduction is solved.

 $\nabla ^{2}T=0$ (4)

The free surface is tracked by the Volume of Fluid (VOF) method, where the volume fraction of the fluid, ε, is tracked through each computational cell. The VOF equation is:

 $\frac{\partial \varepsilon }{\partial t}+V\cdot \nabla \varepsilon =0$ (5)

The interface between fluids is represented by a piecewise linear approach, similar to the work of Youngs [2], in order to greatly limit numerical diffusion of the interface. Surface tension effects are modeled in the numerical simulations. The surface tension forces are represented by the F term in eq. (2). A continuum surface force method proposed by Brackbill et al. [3] is used to model surface tension.

 $F=\sigma \frac{\rho K\nabla \varepsilon }{0.5\left( \rho _{\ell }+\rho _{v} \right)}$ (6)

The curvature, K, is defined as:

 $K=-\frac{\nabla ^{2}\varepsilon }{\left| \nabla \varepsilon \right|}$ (7)

The fluid properties are calculated by the volume weighted average:

 $\varphi =\varepsilon \varphi _{\ell }+\left( 1-\varepsilon \right)\varphi _{v}$ (8)

The general fluid property, $\varphi$, represents either density, viscosity or thermal conductivity. The enthalpy is calculated using a mass-weighted average instead of a volume-weighted average.

 $h=\frac{\left[ \varepsilon \rho _{\ell }c_{\ell }+\left( 1-\varepsilon \right)\rho _{v}c_{v} \right]\left( T-T_{ref} \right)}{\rho }$ (9)

The boundary conditions are as follows: At the inlet $\left( r=r_{in} \right)$

 $\left( 0 (10)
 $\left( -d_{d} (11)

At the collar $\left( x=\delta _{in},\text{ }r_{in}

 $u=v=0,\text{ }w=\omega r,\text{ }\frac{\partial T}{\partial x}=0$ (12)

At the disk surface $\left( x=0,\text{ }r_{in}

 $u=v=0,\text{ }w=\omega r,\text{ }\left. k\frac{\partial T}{\partial x} \right|_{\ell }=\left. k\frac{\partial T}{\partial x} \right|_{S}$ (13)

At the disk bottom $\left( x=-d_{d} \right)$

 $\left( r_{in} (14)
 $\left( r_{hin} (15)
 $\left( r_{hout} (16)

At the outer boundaries $\left( r=r_{hin},\delta _{in},$\left( r_{hin},$\left( r=r_{out},\delta _{in}

 p = pref (17)
 if $\left( V\cdot n>0 \right)$ then $\frac{\partial T}{\partial n}=0$, else T = Tref (18)

At the liquid-vapor interface

 $p_{\ell }-p_{v}-\left. \mu _{\ell }\frac{\partial V\cdot n}{\partial n} \right|_{\ell }+\left. \mu _{v}\frac{\partial V\cdot n}{\partial n} \right|_{v}=\sigma K$ (19)
 $\left. \mu _{\ell }\frac{\partial V\cdot t}{\partial n} \right|_{\ell }=\left. \mu _{v}\frac{\partial V\cdot t}{\partial n} \right|_{v}$ (20)

heating:

 $\left. k_{\ell }\frac{\partial T}{\partial n} \right|_{\ell }=\left. k_{v}\frac{\partial T}{\partial n} \right|_{v}$ (21)

evaporation:

 $\begin{matrix}{}\\\end{matrix}T=T_{sat}$ (22)

Rice et al. [1] solved this problem numerically. Air was used as the vapor for all of the simulations. The viscosity and thermal conductivity of air were lessened by an order of magnitude for selective cases, and were found not to change the results. The boundary condition at the interface is therefore essentially a zero shear stress for all cases, an insulated boundary for heating cases, and Tsat for evaporation.

As was noted before, the flow field is considered to be independent of the thermal field; therefore the flow field is first solved as a time-dependent solution, until a steady state solution has been reached. Once a steady-state solution has been reached, the energy equation is solved for both the fluid and solid regions, as a steady-state solution. The criterion used to determine a steady-state solution was a mass flow rate of liquid leaving the computational domain within 0.05% of the mass flow rate of liquid entering the domain for greater than 0.05 seconds. The flow field was solved using the following procedure:
1. Solve momentum equations
2. Solve continuity (pressure correction) equation, and update pressure and face mass flow rate
3. Repeat steps 1 and 2 until converged, for each time step

A co-located finite volume computational scheme as described in computational method, where both the flow-field variables and the pressure are stored in the cell centers, was used to solve the governing equations. The pressure was discretized in a manner similar to a staggered-grid scheme, while the pressure and velocity were coupled using the SIMPLE algorithm described in the previous section. All of the convective terms in the governing equations were discretized using a second-order upwind scheme.

The grid used for the numerical simulations consisted of rectangular-shaped cells, which were produced in four axial layers in the liquid regions, and one axial layer in the solid region. In the first layer of the fluid region, spanning from the disk surface to δin, there were 25 cell rows, with a growth rate of 1.01. The second layer, spanning from δin to 3δin, consisted of 30 cell rows with a growth rate of 1.02. The third layer, spanning from 3δin to 10δin, consisted of 25 cell rows, and had a growth rate of 1.05. In the final layer, spanning from 10δin to 100δin, there were 20 cell rows, with a growth rate of 1.2. In the solid region, there were 15 evenly spaced cell rows. In the radial direction, there were 15 evenly spaced cell columns between rin and rhin, and 75 evenly spaced cell columns between rhin and rout. The grid spaned such a great distance in the axial direction (100δin), because convergence is greatly increased with the added cells, so the computations ran more rapidly.

Typical results for local film thickness δ and Nussult number, $\text{Nu}={q}''r_{in}/[k_{\ell }(T_{w}-T_{in})]$, as a function of radial distance are presented in the two pictures to the right respectively. The heating boundary condition was used for all cases with an inlet temperature of 20ºC to 40ºC. The evaporation boundary condition was used for all cases with an outlet temperature of 100ºC[4].

## References

1. 1.0 1.1 Rice, J., Faghri, A., and Cetegen, B.M., 2005, “Analysis of a Free Surface Film from a Controlled Liquid Impinging Jet over a Rotating Disk including Conjugate Effects, with and without Evaporation,” Int. Journal of Heat and Mass Transfer, Vol. 48, pp.5192-5204.
2. Youngs, D.L., 1982, “Time-Dependent Multimaterial Flow with large Fluid Distortion,” numerical Methods for Fluid Dynamics, K.W. Morton and M.J. Baines, eds. Academic Press, pp. 273-285.
3. Brackbill, J.U., Kothe, D.B., and Zemach, C., 1992, “A Continuum Method for Modeling Surface Tension,” Journal of Computational Physics, Vol. 100, pp. 335-354.
4. Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.