Application of Computational Methods

From Thermal-FluidsPedia

(Difference between revisions)
Jump to: navigation, search
 
(10 intermediate revisions not shown)
Line 1: Line 1:
-
The computational methods presented in the previous sections 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.
+
{{Comp Method for Forced Convection Category}}
-
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, Figure 4.27, including the conjugate wall effect for both heating with and without evaporation from the free surface.   
+
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<ref name="R2005">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.</ref>.
-
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.  Figure 4.28 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.  
+
[[Image:Fig4.27.png|thumb|400 px|alt=Controlled liquid impinging jet onto a rotating disk|Controlled liquid impinging jet onto a rotating disk]]
 +
[[Image:Fig4.28.png|thumb|400 px|alt=Computational domain for a 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:
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:
Line 18: Line 21:
<math>\rho \frac{DV}{Dt}=-\nabla p+\nabla \cdot \left( \mu \nabla V \right)+\rho \mathbf{g}+\mathbf{F}</math>
<math>\rho \frac{DV}{Dt}=-\nabla p+\nabla \cdot \left( \mu \nabla V \right)+\rho \mathbf{g}+\mathbf{F}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(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.   
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.   
Line 27: Line 30:
<math>\nabla \cdot \left( V\rho h \right)=\nabla \cdot \left( k\nabla T \right)</math>
<math>\nabla \cdot \left( V\rho h \right)=\nabla \cdot \left( k\nabla T \right)</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(3)}}
|}
|}
In the solid region, only conduction is solved.
In the solid region, only conduction is solved.
Line 36: Line 39:
<math>\nabla ^{2}T=0</math>
<math>\nabla ^{2}T=0</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(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:
+
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:
{| class="wikitable" border="0"
{| class="wikitable" border="0"
Line 45: Line 48:
<math>\frac{\partial \varepsilon }{\partial t}+V\cdot \nabla \varepsilon =0</math>
<math>\frac{\partial \varepsilon }{\partial t}+V\cdot \nabla \varepsilon =0</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(5)}}
|}
|}
-
The interface between fluids is represented by a piecewise linear approach, similar to the work of Youngs (1982), 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. (4.327).  A continuum surface force method proposed by Brackbill et al. (1992) is used to model surface tension.
+
The interface between fluids is represented by a piecewise linear approach, similar to the work of Youngs <ref name="Y1982">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.</ref>, 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. <ref name="B1992">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.</ref> is used to model surface tension.
{| class="wikitable" border="0"
{| class="wikitable" border="0"
Line 54: Line 57:
<math>F=\sigma \frac{\rho K\nabla \varepsilon }{0.5\left( \rho _{\ell }+\rho _{v} \right)}</math>
<math>F=\sigma \frac{\rho K\nabla \varepsilon }{0.5\left( \rho _{\ell }+\rho _{v} \right)}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(6)}}
|}
|}
-
The curvature, K, is defined as:
+
The curvature, ''K'', is defined as:
{| class="wikitable" border="0"
{| class="wikitable" border="0"
Line 63: Line 66:
<math>K=-\frac{\nabla ^{2}\varepsilon }{\left| \nabla \varepsilon  \right|}</math>
<math>K=-\frac{\nabla ^{2}\varepsilon }{\left| \nabla \varepsilon  \right|}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(7)}}
|}
|}
The fluid properties are calculated by the volume weighted average:
The fluid properties are calculated by the volume weighted average:
Line 72: Line 75:
<math>\varphi =\varepsilon \varphi _{\ell }+\left( 1-\varepsilon  \right)\varphi _{v}</math>
<math>\varphi =\varepsilon \varphi _{\ell }+\left( 1-\varepsilon  \right)\varphi _{v}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(8)}}
|}
|}
-
The general fluid property,  
+
The general fluid property, <math>\varphi </math>, represents either density, viscosity or thermal conductivity.  The enthalpy is calculated using a mass-weighted average instead of a volume-weighted average.
-
<math>\varphi </math>
+
-
, represents either density, viscosity or thermal conductivity.  The enthalpy is calculated using a mass-weighted average instead of a volume-weighted average.
+
{| class="wikitable" border="0"
{| class="wikitable" border="0"
Line 83: Line 84:
<math>h=\frac{\left[ \varepsilon \rho _{\ell }c_{\ell }+\left( 1-\varepsilon  \right)\rho _{v}c_{v} \right]\left( T-T_{ref} \right)}{\rho }</math>
<math>h=\frac{\left[ \varepsilon \rho _{\ell }c_{\ell }+\left( 1-\varepsilon  \right)\rho _{v}c_{v} \right]\left( T-T_{ref} \right)}{\rho }</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(9)}}
|}
|}
The boundary conditions are as follows:
The boundary conditions are as follows:
-
At the inlet  
+
At the inlet <math>\left( r=r_{in} \right)</math>
-
<math>\left( r=r_{in} \right)</math>
+
-
 
+
{| class="wikitable" border="0"
{| class="wikitable" border="0"
Line 95: Line 94:
<math>\left( 0<x<\delta _{in} \right),\text{ }v=v_{in},\text{ }T=T_{in}</math>
<math>\left( 0<x<\delta _{in} \right),\text{ }v=v_{in},\text{ }T=T_{in}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(10)}}
|}
|}
Line 103: Line 102:
<math>\left( -d_{d}<x<0 \right),\text{ }\frac{\partial T}{\partial r}=0</math>
<math>\left( -d_{d}<x<0 \right),\text{ }\frac{\partial T}{\partial r}=0</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(11)}}
|}
|}
-
At the collar  
+
At the collar <math>\left( x=\delta _{in},\text{ }r_{in}<r<r_{hin} \right)</math>
-
<math>\left( x=\delta _{in},\text{ }r_{in}<r<r_{hin} \right)</math>
+
Line 114: Line 112:
<math>u=v=0,\text{ }w=\omega r,\text{ }\frac{\partial T}{\partial x}=0</math>
<math>u=v=0,\text{ }w=\omega r,\text{ }\frac{\partial T}{\partial x}=0</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(12)}}
|}
|}
-
At the disk surface  
+
At the disk surface <math>\left( x=0,\text{ }r_{in}<r<r_{out} \right)</math>
-
<math>\left( x=0,\text{ }r_{in}<r<r_{out} \right)</math>
+
-
 
+
-
 
+
{| class="wikitable" border="0"
{| class="wikitable" border="0"
Line 126: Line 121:
<math>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}</math>
<math>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}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(13)}}
|}
|}
-
At the disk bottom  
+
At the disk bottom <math>\left( x=-d_{d} \right)</math>
-
<math>\left( x=-d_{d} \right)</math>
+
-
 
+
{| class="wikitable" border="0"
{| class="wikitable" border="0"
Line 137: Line 130:
<math>\left( r_{in}<r<r_{hin} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}=0</math>
<math>\left( r_{in}<r<r_{hin} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}=0</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(14)}}
|}
|}
Line 145: Line 138:
<math>\left( r_{hin}<r<r_{hout} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}={q}''_{w}</math>
<math>\left( r_{hin}<r<r_{hout} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}={q}''_{w}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(15)}}
|}
|}
Line 153: Line 146:
<math>\left( r_{hout}<r<r_{out} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}=0</math>
<math>\left( r_{hout}<r<r_{out} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}=0</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(16)}}
|}
|}
-
At the outer boundaries
+
At the outer boundaries <math>\left( r=r_{hin},\delta _{in}<x<100\delta _{in} \right)</math>,<math>\left( r_{hin}<r<r_{out},x=100\delta _{in} \right)</math>,<math>\left( r=r_{out},\delta _{in}<x<100\delta _{in} \right)</math>
-
<math>\left( r=r_{hin},\delta _{in}<x<100\delta _{in} \right)</math>,<math>\left( r_{hin}<r<r_{out},x=100\delta _{in} \right)</math>,<math>\left( r=r_{out},\delta _{in}<x<100\delta _{in} \right)</math>
+
{| class="wikitable" border="0"
{| class="wikitable" border="0"
Line 163: Line 155:
<math>p=p_{ref}</math>
<math>p=p_{ref}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(17)}}
|}
|}
-
 
{| class="wikitable" border="0"
{| class="wikitable" border="0"
|-
|-
| width="100%" |<center>
| width="100%" |<center>
-
<math>\text{if }\left( V\cdot n>0 \right)\text{ then }\frac{\partial T}{\partial \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\rightharpoonup}$}} {n}}=0\text{, else }T=T_{ref}</math>
+
if <math>\left( V\cdot n>0 \right)</math> then <math>\frac{\partial T}{\partial n}=0</math>, else <math>T=T_{ref}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(18)}}
|}
|}
At the liquid-vapor interface
At the liquid-vapor interface
Line 180: Line 171:
<math>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</math>
<math>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</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(19)}}
|}
|}
Line 188: Line 179:
<math>\left. \mu _{\ell }\frac{\partial V\cdot t}{\partial n} \right|_{\ell }=\left. \mu _{v}\frac{\partial V\cdot t}{\partial n} \right|_{v}</math>
<math>\left. \mu _{\ell }\frac{\partial V\cdot t}{\partial n} \right|_{\ell }=\left. \mu _{v}\frac{\partial V\cdot t}{\partial n} \right|_{v}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(20)}}
|}
|}
heating:
heating:
Line 197: Line 188:
<math>\left. k_{\ell }\frac{\partial T}{\partial n} \right|_{\ell }=\left. k_{v}\frac{\partial T}{\partial n} \right|_{v}</math>
<math>\left. k_{\ell }\frac{\partial T}{\partial n} \right|_{\ell }=\left. k_{v}\frac{\partial T}{\partial n} \right|_{v}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(21)}}
|}
|}
evaporation:
evaporation:
Line 204: Line 195:
|-
|-
| width="100%" |<center>
| width="100%" |<center>
-
<math>T=T_{sat}</math>
+
<math>\begin{matrix}{}\\\end{matrix}T=T_{sat}</math>
</center>
</center>
-
|{{EquationRef|(1)}}
+
|{{EquationRef|(22)}}
|}
|}
-
Rice et al. (2005) solved this problem numerically with the result presented below.
+
Rice et al. <ref name="R2005"/> 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.
-
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:<BR>
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:<BR>
-
1. Solve momentum equations<BR>
+
1. Solve momentum equations<BR>
-
2. Solve continuity (pressure correction) equation, and update pressure and face mass flow rate<BR>
+
2. Solve continuity (pressure correction) equation, and update pressure and face mass flow rate<BR>
-
3. Repeat steps 1 and 2 until converged, for each time step<BR>
+
3. Repeat steps 1 and 2 until converged, for each time step<BR>
-
A co-located finite volume computational scheme as described in Section 4.8, 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.
+
A co-located finite volume computational scheme as described in [[Computational methodologies for forced convection |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 ''δ<sub>in</sub>'', there were 25 cell rows, with a growth rate of 1.01.  The second layer, spanning from ''δ<sub>in</sub>'' to 3''δ<sub>in</sub>'', consisted of 30 cell rows with a growth rate of 1.02.  The third layer, spanning from 3''δ<sub>in</sub>'' to 10''δ<sub>in</sub>'', consisted of 25 cell rows, and had a growth rate of 1.05.  In the final layer, spanning from 10''δ<sub>in</sub>'' to 100''δ<sub>in</sub>'', 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 ''rh<sub>in</sub>'' and ''r<sub>out</sub>''.  The grid spaned such a great distance in the axial direction (100''δ<sub>in</sub>''), because convergence is greatly increased with the added cells, so the computations ran more rapidly.
 +
 
 +
Typical results for local film thickness ''δ'' and Nussult number, <math>\text{Nu}={q}''r_{in}/[k_{\ell }(T_{w}-T_{in})]</math>, 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<ref>Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.</ref>.
 +
 
 +
==References==
 +
{{Reflist}}

Current revision as of 02:32, 27 July 2010

 Related Topics Catalog
Computational methodologies for forced convection
  1. One-Dimensional Steady-State Convection and Diffusion
    1. Central Difference Scheme
    2. Upwind Scheme
    3. Hybrid Scheme
    4. Exponential and Power Law Schemes
    5. A Generalized Expression of Discretization Schemes
  2. Multidimensional Convection and Diffusion Problems
  3. Numerical Solution of Flow Field
    1. Special Difficulties
    2. Staggered grid
    3. Pressure Correction Equation
    4. The SIMPLE Algorithm
  4. Numerical Simulation of Interfaces and Free Surfaces
  5. 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
Controlled liquid impinging jet onto a rotating disk
Computational domain for a 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<x<\delta _{in} \right),\text{ }v=v_{in},\text{ }T=T_{in}

(10)

\left( -d_{d}<x<0 \right),\text{ }\frac{\partial T}{\partial r}=0

(11)

At the collar \left( x=\delta _{in},\text{ }r_{in}<r<r_{hin} \right)


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}<r<r_{out} \right)

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}<r<r_{hin} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}=0

(14)

\left( r_{hin}<r<r_{hout} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}={q}''_{w}

(15)

\left( r_{hout}<r<r_{out} \right):\text{ }\left. k\frac{\partial T}{\partial x} \right|_{S}=0

(16)

At the outer boundaries \left( r=r_{hin},\delta _{in}<x<100\delta _{in} \right),\left( r_{hin}<r<r_{out},x=100\delta _{in} \right),\left( r=r_{out},\delta _{in}<x<100\delta _{in} \right)

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.