PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINE return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

moo chlMpSS—p.“

 

NUMERICAL SIMULATIONS OF A DEVELOPING TURBULENT WALL JET
ALONG A THERMALLY ACTIVE SURFACE

By

Qingtian Wang

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Mechanical Engineering

1 999

ABSTRACT

NUMERICAL SIMULATIONS OF A DEVELOPING TURBULENT WALL JET
ALONG A THERMALLY ACTIVE SURFACE

By

Qingtian Wang

In support of defogging studies of automotive vehicles, a 20 steady state
numerical simulation and phase change modeling for a developing wall jet
experimental facility were performed with a commercially available CFD code,
FLUENT.

The velocity, temperature, and turbulence fluctuations were predicted with
the standard k-e, the realizable k-s, and the Reynolds stress turbulence models.
The simulation results were compared to experimental measurements after they
were outer and inner scaled. The numerical predictions are in good agreement
with the experimental data in outer scaling for isothermal and non-isothermal
simulations without phase change.

Two phase change models were developed to perform the mass transfer
simulation. The surface phase model was implemented, and gives promising
water vapor condensation predictions in comparison with experimental drip-off
tests. The volumetric phase change model is analytically established. and needs

future work on implementation.

Copyright by
Qingtian Wang
1 999

To Xiaohui and my family

iv

ACKNOWLEDGMENTS

I wish to express my heartfelt thanks to my advisor, Dr. John J. McGrath,
for his constructive supervision and constant support during the entire program.

I would like to thank my committee members, Dr. Tom Shih and Dr.
Abraham Engeda, for their advice and assistance.

The experimental data in the present study were collected and processed
by Paul Hoke and Rene Neust. I must thank them as their generous contribution
makes all the experimental confirmations possible. Meanwhile, I would like to
thank Dr. Bashar AdbdulNour and Ford Motor Company for initializing the
geometry for the 20 model and the continuous technical support during the
course of the study.

I want to thank some of my fellow graduate students for their help: Bin
Lian, Gloria Elliot, AI Aksan, and Shi Zheng.

I also want to thank Dr. Jayanta Sanyal and Dr. Aniruddha Mukhopadhyay
at FLUENT, INC., for their help in developing the user defined functions.

TABLE OF CONTENTS

TABLE OF CONTENT ........................................................................................ VI
LIST OF TABLES ............................................................................................... IX
LIST OF FIGURES .............................................................................................. X
LIST OF SYMBOLS .......................................................................................... XII
CHAPTER 1 INTRODUCTION AND BACKGROUND INVESTIGATION ......... 1
1.1 PROBLEM STATEMENT ................................................................. 1
1.2 BACKGROUND INVESTIGATION AND PREVIOUS WORK ........... 3
1.2.1 Wall Jet: Basic Illustration ................................................... 3

1.2.2 Previous Work: Velocity Distribution and
Turbulence Quantities ....................................................... 5

1.3

CHAPTER2
2. 1

2.2
2.3

2.4

1.2.3 Previous Work: Temperature Distribution and

Phase Change ................................................................... 6
SUMMARY .................................................................................... 10
NUMERICAL METHODOLOGY ................................................. 12
BASIC MODELS AND EQUATIONS ............................................. 12
2.1.1 Continuity and Momentum Equations ............................... 12
2.1.2 Energy Equation ............................................................... 13
2.1.3 Species Transport Equation (Concentration Equation) ..... 15
DISCRETIZATION AND MESH GENERATION ............................ 16
TURBULENCE MODELS .............................................................. 19
2.3.1 Standard k-e Model ........................................................... 22
2.3.2 Realizable k-e Model ........................................................ 27
2.3.3 Reynolds Stress Model (RMS) ......................................... 30
NEAR WALL TREATMENT ........................................................... 31

2.4.1 The Turbulence Boundary Layer and

vi

the Law of the Wall ........................................................... 32

2.4.2 Near Wall Modeling in FLUENT ........................................ 36
2.5 SUMMARY .................................................................................... 39
CHAPTER 3 SIMULATION WITHOUT PHASE CHANGE ............................. 40
3.1 SOLUTION PROCEDURE ............................................................ 40
3.1 .1 Pre-processing ................................................................. 41
3.1.2 Model Control ................................................................... 42
3.1.3 Near Wall Treatment ........................................................ 43
3.1.4 Unit System ...................................................................... 43
3.1.5 Material Properties ........................................................... 44
3.1.6 Operating Conditions ........................................................ 45
3.1.7 Boundary Conditions ........................................................ 45
3.1.8 Solution Control ................................................................ 47
3.1.9 Convergence Judgement ................................................. 48

3.2 RESULTS AND DISCUSSIONS OF THE
ISOTHERMAL SIMULATION ........................................... 49
3.2.1 Velocity Profiles ................................................................ 50
3.2.2 Turbulence Fluctuations ................................................... 54

3.3 RESULTS AND DISCUSSIONS OF THE
NON-ISOTHERMAL SIMULATION .................................. 57
3.4 SUMMARY ..................................................................................... 61
CHAPTER 4 SIMULATIONS WITH PHASE CHANGE MODELS .................. 87
4.1 SOME BASICS ABOUT MASS TRANSFER .................................. 88
4.1 .1 General Mixture ................................................................ 88
4.1.2 Water Vapor and Dry Air .................................................. 90
4.2 SURFACE PHASE CHANGE MODEL ........................................... 91
4.2.1 Model Analysis ................................................................. 92
4.2.2 Model Implementation ...................................................... 94
4.2.3 Results and Comparisons ................................................ 97
4.3 VOLUMETRIC PHASE CHANGE MODEL ................................... 108

vii

4.3.1 Model Analysis ............................................................... 108

4.3.2 Model Implementation .................................................... 1 13
4.3.3 Results and Comparisons .............................................. 116
4.4 SUMMARY AND DISCUSSION ................................................... 116

CHAPTER 5 COCLUSIONS AND RECOMMENDATIONS

 

FOR FUTURE WORK .................................................... 1 13
5.1 SIMULATIONS WITHOUT PHASE CHANGE .............................. 113
5.2 SIMULATIONS WITH PHASE CHANGE ..................................... 119
APPENDICES ............... .................................................... 122
A. SIMULATIONS OF OUTER-SCALED SNS ................................... 123
B. SIMULATIONS OF INNER-SCALED
FLUCTUATIONS WITH RSM ....................................................... 124
C. SIMULATION OF TEMPERATURE PROFILS
(OUTER SCALING) ....................................................................... 126
D. SIMULATION OF THE lNNER-SCALED TEMPERATURE
PROFILE AND THE HEAT TRANSFER COEFFICIENT ............... 127
E. UDF TO IMPOSE SATURATION MASS FRACTION AS
BOUNDARY CONDITION AND CALCULATE THE VAPOR MASS
FLUX (SURFACE PHASE CHANGE MODEL) ............................ 123
F. COMPARISON OF NUMERICAL PREDITION WITH THE
SIMILARITY SOLUTION (LAMINAR FLOW, FLAT PLATE) ........ 130
G. PREDICTION OF MOLE FRACTION PROFILES OF WATER
VAPOR BY THE SURFACE PHASE CHANGE MODEL .............. 131
H. USER DEFINED FUNCTION FOR THE VOLUMETRIC PHASE
CHANGE MODEL (METHOD (1)) ................................................. 134
I. USER DEFINED FUNCTION FOR THE VOLUMETRIC PHASE
CHANGE MODEL (METHOD (2)) ................................................. 133
LIST OF REFERENCES ................................................................................... 142

viii

 

LIST OF TABLES

TABLE 1 DEFAULT VALUES OF PHYSICAL PROPERTIES
OF AIR IN FLUENT ......................................................................... 44

LIST OF FIGURES

FIGURE PAGE
1. ILLUSTRATION OF EXPERIMENTAL SETUP ............................................ 2
2. CONFIGURATION OF A 2-D WALL JET ..................................................... 4
3. THERMAL CONDITIONS OF A WALL JET ................................................. 5
4. SCHEMATIC OF A HEAT/MASS TRANSFER PROBLEM
OVER A FLAT PLATE ........................................................................ 7
5. CONTROL VOLUME USED FOR DISCRETIZATION ............................... 17
6 ZONES IN THE TURBULENCE BOUNDARY LAYER
(REPRODUCED FROM [13]) ............................................................ 33
7. GEOMETRY AND BOUNDARY TYPES OF THE SIMULATION .............. 42
3 DEVELOPMENT OF VELOCITY PROFILE .............................................. 63
. PLOT OF VELOCITY MAGNITUDE ......................................................... 63
1o. CONTOURS OF VELOCITY MAGNITUDE .............................................. 64
11. CONTOURS OF STATIC PRESSURE ..................................................... 64
12. VISUALIzATION OF VELOCITY VECTORS ............................................ 65
13. VISCOSITY-AFFECTED SUBLAYER (REY <200) .................................... 66
14. COMPARISON OF VELOCITY PROFILES
(OUTER SCALING, ske, A) ............................................................ 67
15. COMPARISON OF VELOCITY PROFILES
(OUTER SCALING, ske, B) ............................................................ 63
16. COMPARISON OF VELOCITY PROFILES
(OUTER SCALING, rke, A) ............................................................. 69
17. COMPARISON OF VELOCITY PROFILES
(OUTER SCALING, rke, B) ............................................................. 7o
13. COMPARISON OF VELOCITY PROFILES
(OUTER SCALING, RSM, A) .......................................................... 71
19. COMPARISON OF VELOCITY PROFILES
(OUTER SCALING, RSM, B) .......................................................... 72

20.
21.
22.
23.
24.

25.

26.
27.
28.
29.

30.

31.
32.

33.

35.
36.

37.

38.
39.
40.
41.
42.

COMPARISON OF TURBULENT VISCOSITY RATIO (A) ....................... 73

COMPARISON OF TURBULENT VISCOSITY RATIO (3) ....................... 74
COMPARSION OF VELOCITY PROFILES (INNER SCALING) ............... 75
COMPARISON OF 0* DISTRIBUTION ..................................................... 76
OVERALL INNER-SCALED PREDICTIONS OF

0* DISTRIBUTION (ske) ................................................................. 77
OVERALL INNER-SCALED PREDICTIONS OF

0" DISTRIBUTION (rke AND RSM) ................................................. 73
COMPARISON OF SNS (OUTER SCALING, A) ...................................... 79
COMPARISON OF SNS (OUTER SCALING, B) ...................................... 30
COMPARISON OF u" (INNER SCALING, A) .......................................... 31
COMPARISON OF II'+ (INNER SCALING, B) .......................................... 32
COMPARISON OF 7* (INNER SCALING) .............................................. 33
COMPARISON OF W (INNER SCALING) ........................................... 34
COMPARISON OF TEMPERATURE PROFILES

IN OUTER SCALING ...................................................................... 35
COMPARISON OF lNNER-SCALED TEMPERATURE

PROFILES AT W = 10.3 .............................................................. 36
20 FLAT PLATE GEOMETRY .................................................................. 97
ILLUSTRATION OF FACE AND CELL VALUES .................................... 100
COMPARISON OF MASS TRANSFER

COEFFICIENT (LAMINAR, FLAT PLATE) ................................... 101
RATIO OF ANALYTICAL To NUMERICAL

SOLUTIONS FOR h... (LAMINAR, FLAT PLATE) ........................ 101
PREDICTIONS OF VAPOR MASS FLUX ON THE WALL ..................... 103
COMPARISON OF RELATIVE HUMIDITY AT XNV =10.3 ...................... 104
COMPARISON OF CONDENSATION RATE OF THE PLATE ............... 105
COMPARISON OF SHERWOOD NUMBER ........................................... 103
CONTROL VOLUME FOR PHASE CHANGE ANALYSIS ...................... 111

xi

lDifim
DI

Psat

Pr

qllapp

qflo

LIST OF SYMBOLS

specific heat at constant pressure (J/kg-K)

mass diffusivity (mzls)

diffusion coefficient of species i' in the mixture (mzls)
effective diffusion coefficient due to turbulence (m2/s)
total internal energy (J/m3)

heat transfer coefficient (W/mZ-K)

local mass transfer coefficient (m/s)

mass diffusion flux of species i' in i direction (kg/mz-s)
mass flux at the wall (kg/mz-s)

kinetic energy and its dissipation rate

effective thermal conductivity (W/m-K)

Lewis number Le = L33 = 9.
PT D

mass of species i' (kg)
volumetric rate of the creation of the water vapor (kg/ma-s)

turbulent Mach number
saturation pressure (Pa)

partial pressure of vapor in the mixture (Pa)
Prandtl number Pr = X—
a

apparent heat flux (W/mz)

heat flux at y = 0 (wall surface) (W/mz)

xii

Rey
Raw
3'.
Sm

SC

SCI

Sh

SNS

Umax
U1/2

mass source (or sink) of species i' (kg/m3-s)
turbulent Reynolds number

Reynolds number based on the jet nozzle width
mass generation rate of species i' (kg/ma-s)

source term in the mass continuity equation (kg/m3-s)

Schmidt number Sc = %

Sc = fi— is the turbulent Schmidt number (default value 0.7)

I
I

 

Sherwood number based on the width of jet nozzle w, D

—'2

u

2
max

streamwise normal stress, SNS = 100

 

normalized temperature

dew point temperature (K)

temperature at the jet nozzle (K)

temperature of the plate surface (K or °C)
reference temperature (K)

velocity components in x, y and z direction (m/s)
x and y velocity components of water vapor (m/s)
friction velocity (m/s)

velocity magnitude (m/s)

maximum velocity magnitude (m/s)

one half of the maximum velocity magnitude (m/s)

xiii

 

U191 velocity at the jet nozzle (Ms)
W width of the jet nozzle (slot width) (m)

x‘, y", u*, v“ wall coordinates

x) mole fraction of species i'

ym y location where the velocity takes the maximum value in the profile
(M)

y', u' wall units

Greek Letters

 

a thermal diffusivity a = k (mzls)
,.

B coefficient of thermal expansion

62 y location where velocity magnitude takes the value of U/2 (m)

6ij kronecker delta

3M momentum eddy diffusivity (m2/S)

3H thermal eddy diffusivity (mZ/s)

4) relative humidity

d).- mass fraction of species i'

(I)... main stream mass fraction

)3 dynamic viscosity (kg/s-m)

pea effective viscosity (kgIs-m)

pt turbulent viscosity (kgls-m)

v kinetic viscosity (mzls)

xiv

p...
to

I);

density of the mixture (kg/m3)

mass concentration of species i' (kg/m3)

mass concentration at the wall (kg/m3)

mass concentration in the main stream (kg/m3)
shear stress at the surface of the wall (Pa)

stress tensor (Pa)

Chapter 1

INTRODUCTION AND BACKGROUND INVESTIGATION

1.1 Problem Statement

The motivation for the present study was to provide a Computational Fluid
Dynamics (CFD) simulation of experimental results related to basic research of
the defogger in motor vehicles. Defogging effects on the interior windshield
surface is an important safety factor for automobiles. It is generally investigated
by experimental tests. Meanwhile, with the development of CFD techniques and
its obvious benefits compared to conventional experimental methods, it is
desirable to perform satisfactory numerical simulations and predictions to
improve the defogger design. The current study is serves as an initial step
towards this goal.

An experimental facility was constructed at the Heat Transfer Lab in
Department of Mechanical Engineering, Michigan State University [1]. It defines
the geometry and boundary conditions for the CFD computation. Therefore, the
data Obtained from this facility are utilized as experimental verification of the CFD
predictions. A commercially available CFD code FLUENT (unstructured version
5.0) was adopted as the tool to implement the numerical Simulations. The code
runs on a Windows NT 4.0 platform. The hardware is a DELL PC with a 400MHz
Pentium II CPU and 256MB RAM.

This thesis is comprised of two foci:

(1) Prediction of the steady state wall jet velocity and temperature
distributions with currently available numerical models that come with the
commercial solver; and

(2) Theoretical analysis and modeling of the fogging phase change
phenomenon and its numerical implementation with FLUENT.

The experimental setup produced by Hoke [1] is depicted in Figure 1.

Flow

System Inlet

Steam Supply
Armstrong Model 90
Steam Humidifier
Condensate Return
Duct

Flow Conditioner

9’91“ 9N2"

(Honeycomb and Screens)

7. 2—Dimensional Contraction
8. Splitter Plate
9. Thermally Active Test Section
10. Sub Atmospheric Test Chamber
11. Probe Traverse System
12. Calibration Nozzle

‘0 13. Wind Tunnel Base

14. Fluid Collection

15. Prime Mover

16. Exit Throttle Valve

16 z

  
 
         
   

11

Wall jet
Flow
13

‘4 Flow

Figure 1. Illustration of Experimental Setup
The region of interest for both the numerical and experimental
investigations is the thermally active test section 9 (there are actually two
identical regions due to the symmetry of the facility). The design of the

experimental facility makes it possible to provide an adjustable flow pattern in this

region. The driving force for the whole system is the pressure differential
between the system inlet 1 and wind tunnel base 13 created by the prime mover
15. The pressure driven flow forms two identical wall jets along the thermally
active plate surface in 9 after passing through the contraction 7 and being
separated by the splitter plate 8. The pressure differential, the temperature of the
plate surface and the relative humidity of the flow can all be adjusted according
to the experimental needs [1].

The fluid studied is a mixture of water and non-condensable air. The
numerical simulation includes investigation of the velocity field, the turbulence
quantities, and the temperature field. Meanwhile, since the plate surface
temperature T3 is adjustable, a phase Change (condensation or evaporation) of
water will occur under certain thermal circumstances. When T3 is below the dew
point of water, for instance, there will be liquid water condensing onto the wall
surface, accompanied with latent heat release. This is of most concern in
windshield defogger designing. In turn, modeling the phase Change phenomenon

becomes very necessary.

1.2 Background Investigation and Previous Work

1.2.1 Wall Jet: Basic Illustration

A wall jet is defined by Launder and Rodi [3] as a shear flow directed
along a wall where, by virtue of the initially supplied momentum, at any station,
the streamwise velocity over some region within the shear flow exceeds that in

the external stream. Some other researchers define a wall jet as a jet of fluid

impinging onto a wall at an angle from 0 to 90 degrees [3] [4]. As a first step of
defogger research, the impinging angle of the currently studied wall jet is zero.
Figure 2 defines some basic velocity field configurations of a two
dimensional wall jet flowing tangentially to the wall surface (zero impinging
angle). In the figure, w is the width of the jet exit nozzle, U)... is the velocity
(magnitude) at the Slot exit, Umax is the maximum velocity of the velocity profile,
Um = Umax/ 2, ylm is the location where the velocity U = Um... and 82 is the y

location where U = Um».

Q

T- .. U10- Unun

E\\\\\\\\\\\\\\ r- (q

 

 

 

 

5:
////° //// .//// f/// 7‘

Figure 2. Configuration of a 2-D Wall Jet
If the temperature of the plate surface is different from that of the main
stream of the wall jet, the temperature distribution in the wall jet will be of interest.
Figure 3 shows a schematic view of the thermal conditions of a typical wall jet

with Ts < Tjet, where Ts is the temperature of the isothermal wall surface, Tjet is

the temperature at the jet nozzle, T00 is the far field temperature, and Tmax is the

maximum temperature.

Q

 

E\\\\\\\\\\\X\ ._ .4
ll
"—1
It

 

L

 

 

 

Them ally Active Wall Swface

////o //// //// //// *"

Figure 3. Thermal Conditions of a Wall Jet
Note that the fluid of the wall jet is a mixture of condensable water and
non-condensable air. The concentration distribution of water vapor is of
importance from the standpoint of phase change and mass transfer. The mass
flux on the wall surface due to condensation/evaporation is especially of interest

for the study of a defogger.

1.2.2 Previous Work: Velocity Distribution and Turbulence Quantities
The nature of a wall jet has been explored by many former researchers.
Launder and Rodi [2] gave the most comprehensive review of the previous

work on the flow field up to 1980. Their attention was mainly focused on the

velocity field and turbulence quantities. Heat and mass transfer were not
included.
In 1992, Wygnanski, Katz and Horev [5] published their experimental data

in the region where the downstream distance is greater than 20 slot width

(x Z 20W). They claimed that the bulk of the flow is self-Similar by outer scaling
(U/Umax as a function of y/62), where the term self-similar (or self-preserving)
means that a dimensionless quantity depends on only one lateral or transverse
dimensionless variable. Hence, the distribution preserves its normalized shape
as the flow proceeds downstream [6]. It is believed that the wall jet was well
developed for the case of their research, which is beyond the region of interest
for the current study.

Eriksson, Karlson and Persson [7] published velocity and turbulence data

for a wider region in 1998. In the lateral direction, their data ranged between 0 to

150 slot widths (0 S x S 150W ). In the transverse direction, they paid special
attention to the near wall region and collected data below y+ = 4 (y+, u+ are inner
scaling coordinates) by using the Laser-Doppler (LDV or LDA) technique. The x
and y components Of the velocity vector can be distinguished with this technique.
They found that the data in the range 40 S x/ W S 150 was reasonably consistent
with similarity, claiming that the mean velocity profiles were self-similar in inner
scaling up to y+ = 100-200, and that the turbulence quantities showed similarity

in a much shorter range in terms of y+.

1.2.3 Previous Work: Temperature Distribution and Phase Change

As mentioned previously, if the temperature of the wall surface is different
from that of the free stream, the temperature distribution is of interest. Moreover,
note that the fluid studied is a mixture of condensable water and non-
condensable air. The concentration of water vapor is influenced by the
temperature distribution and phase change (condensation/evaporation). For the
study Of the windshield defogger, the phase change and mass transfer
phenomena, especially that on the wall surface and in the near-wall region, are
issues of utmost importance.

A simple situation is the so-called ”forced convection condensation on a
flat plate in the presence of a non-condensable gas”, where, in the current study
the non-condensable gas is dry air. This is shown in Figure 4, which also shows

the dimensional nomenclature and coordinates.

 

 

 

 

U... T... m.
——-—ev-
Y
Water- Air Mixture
3 5!. Liqu'd Film
x
o \ P ]

 

Ts
Figure 4. Schematic of a heat/mass transfer problem over a flat plate

U”, T00 and mo are free stream velocity, temperature and mass fraction of
water vapor, respectively. BL is the thickness of the liquid water film on the plate
surface due to condensation.

Sparrow, Minkowycz and Saddy [8] worked on the case of laminar flow in

the 19608. The liquid film on the plate surface was included in their analysis by

both numerical and integral methods. It was assumed that the vapor was
saturated in the main stream and at the interface between the gaseous phase
and liquid film. Furthermore, they assumed the phase change only occurred at
the interface. The possible volumetric condensation of water vapor to water
droplets in the boundary layer of the water-air mixture was not considered.
Therefore, no alteration was needed in the continuity, energy and concentration
(species diffusion) equations for the binary gaseous mixture. They claimed that
the resistance of the liquid—gas interface to the forced condensation is negligible.
Another investigation of the laminar case was by Hijikata and Mari [9] in
1973. The possibility of volumetric condensation of vapor (i.e. the fog formation)
was considered to satisfy the local thermodynamic equilibrium of the water-air
mixture. Therefore, the mixture became a two-species (air and water), two-phase
(gaseous air-vapor and liquid water droplets) system in their model.
Corresponding alterations in the governing equations were made to capture the
volumetric phase change (condensation of vapor to liquid droplets) phenomenon.
The boundary layer equations were solved by the integral method, assuming the
liquid film is ”very thin" so that the temperature of the liquid film and that of the
interface were assumed equal to the temperature of the plate surface. Their
analysis showed that, for a vapor-air mixture with a small temperature difference
between the free stream and the plate surface (less than 50°C), the effect of
volumetric condensation on the heat and mass flux is dependant upon the
relationship between the temperature and concentration difference. If the

difference of temperature is dominant compared to that of concentration, the heat

flux to the liquid film will not be altered very much by the volumetric
condensation, and the mass flux is smaller for the case with condensation than
that without. In other words, for the case with thermodynamic equilibrium
considered, the temperature distribution does not vary drastically, while the slope
of the concentration profile becomes smaller away from the plate surface.
Othenlvise, if the concentration difference (i.e. the difference of water vapor mass
fraction between the main stream and plate surface) plays a more important role,
the heat flux to the liquid film is mainly produced by the latent heat released due
to the phase Change of mass flux.

The same “thin film” assumption was adopted by Legay-Desesquelles and
Prunet-Foch [10]. By using the finite difference method, they numerically solved
both laminar and turbulent cases of the problem with volumetric phase change
considered. The numerical calculation was compared to experimental data. Their
calculation showed that in the case when phase change occurred in the
boundary layer, larger heat and mass transfer at the wall were predicted
compared to the case of dry air. These results appear to be consistent with the
case in Hijikata and Mori's study, in which the concentration difference is
dominant compared to the difference of temperature. This seems reasonable
since the temperature difference in Legay—Desesquelles and Prunet-Foch's study
is less than 20 °C.

The most comprehensive analysis of the laminar case was carried out by
Matuszkiewicz and Vernier [11] in 1991. Instead of adopting the " thin film"

assumption, the effect of liquid film was taken into account analytically.

Moreover, the liquid phase in the main stream was modeled by analyzing the
mass fraction of liquid water (droplets). The concept of "source terms" in the
species diffusion and energy equations was explicitly introduced to take care of
the overall phase Change phenomenon. Their results showed that for fluid with
Lewis number Le < 1 (which is the case for vapor-air mixture), the saturated
condition does not hold all the way through the boundary layer unless the droplet
(water liquid) mass fraction in the main stream has a certain minimum value. In
other words, if the droplet mass fraction is less than the minimum value, there
exits regions in the boundary layer where the water vapor is superheated. For Le
> 1, the saturated condition holds throughout the boundary layer whatever value
the main stream droplet mass fraction takes. They found that the heat flux to the
wall surface is only slightly influenced by the presence of droplets in the
condensation boundary layer. Also, the "thin film” assumption is only valid when
the vapor mass fraction is very low. That is, the assumption is valid for the air-
vapor system since the mass fraction of water vapor is usually less than 0.1.

It is also worth mentioning that in the analysis of all the researchers, the
assumption that the interface is impermeable to non-condensable air was made.
Therefore, the transverse velocity component of non-condensable air was set to

be zero for the boundary (or jump) condition at the interface.

1.3 Summary

As can be seen from the previous analysis, the task of simulating a wall jet

is complicated, especially when it is coupled with heat transfer and phase change

10

phenomenon. Extensive modeling work is needed for the coupling of heat and
mass transfer. The goal of the current study is to tackle this task numerically with
the commercially available tool FLUENT. A better understanding of this problem
through the CFD study can serve as the first step towards the optimization of

windshield defogger design.

11

Chapter 2

NUMERICAL METHODOLOGY

The commercial software FLUENT is utilized as a workbench to carry
numerical computation and implement necessary physical models developed by
the user. It should be pointed out that it is not the focus of the current study to
improve the numerical methodology (such as the turbulence models) of the code.
Nevertheless, the comparison between the simulation and experimental data
may provide suggestions on choosing the appropriate numerical methods
(turbulence models, for instance).

It is beneficial, however, to know what is available in commercial CFD
package. Following is a brief introduction to the relative numerical methodology

employed by FLUENT, and some necessary background physical theories.

2.1 Basic Models and Equations

2.1.1 Continuity and Momentum Equations

FLUENT solves conservation equations for mass and momentum for all
flows. For flows with heat transfer, the energy conservation equation is solved.
For flows involving species which are mixing or reacting, FLUENT solves the
species conservation equation. Additional transport equations are also solved
when the flow is turbulent.

The continuity equation, or conservation equation of mass, solved by

FLUENT takes the form [12]:

12

28+3 .I
a: 6x,

s (2.1)

"I

This is a general form of continuity equation and is valid for both
incompressible and compressible flows. Sm (kg/ma-s) includes the source term of
mass added to the continuous phase from the dispersed second phase (eg. due
to vaporization of liquid droplets) and any user-defined sources.

Later in the current study, a user-defined source term will be introduced to
take into account the condensation/evaporation of water.

In an inertial (non-accelerating) system, the momentum equation (Navier—

Stokes Equation) in the idirection solved by FLUENT is:

a a 07,--
5W1)+§Wiui)=’%+§+%i+fi (22)

where p is the static pressure, pg; is the gravitational body force, F. (Nlm3) are the
external body forces in the idirection and other model-dependent source terms
such as porous media user-defined sources, and r.) is the stress tensor given by

[13]:

Du. Du]. 2 Oak
..= __‘+__ __ _§.. 2.3
TU ,u[ j i] 3;: k U ( I

where u is the dynamic viscosity and 8;) is the Kronecker delta function (6., =1 if i =

j and 6., = 0 if i at: j). The second term on the right-hand side is the effect of volume

dilation.

2.1.2 Energy Equation

13

FLUENT can solve a heat transfer problem within both fluid and solid
regions in the computational domain. The energy equation solved in FLUENT
takes the following form:

gt—(pE)+axii[u,(pE + p)]= 36,-[k‘4’ gx—TI—ghflj. +uj(r,j)fl] + s, (2.4)
S), (W/m3) is the source term for Chemical reaction, and any other user-defined
volumetric source. The first three terms on the right-hand side of equation (2.4)
represent heat transfer due to conduction, species diffusion and viscous
dissipation, respectively. ken is the effective thermal conductivity; and (Tij)eff is the
deviatoric stress tensor. They will be discussed more in turbulence models. J) is
the diffusion flux of species j'.

In equation (2.4), E is the total energy

2
E=h-£+Ei"— (2.5)
p

where h is the sensible enthalpy. If (D) is the mass fraction of species j', h is
defined as

h = 20 fh 1.. (2.6)
for an ideal gas. And for incompressible flow,

h = 20,72) + £ (2.7)
1' .0

where h) is given by

T
h, = [ 7:”.de (2.3)

T",

The reference temperature TM in FLUENT is 298.15 K.

14

2.1.3 Species Transport Equation (Concentration Equation)

FLUENT models the mixing and transport of chemical species by solving
the conservation equation of convection, diffusion and reaction sources for each
individual species. Therefore, it can be used to model Species transport with or
without chemical reaction. The four reaction modeling approaches available are:
- Generalized finite formulation
I Mixture fraction/PDF formulation
- Non-equilibrium (flamelet) formulation
- Premixed combustion formulation
None of these models will be used for the present problem.

FLUENT predicts the local mass fraction for each species, (Dr, by solving
the convectionodiffusion equation for the i'th species. The general form of this

equafionis

—a-(p¢. )+—a—(pu.<l> )=—axiJ,., +R +S (2.9)

1

Equation (2.9) is also known as the "concentration equation” [14], in which
R;- is the mass source (or sink) due to chemical reaction, and Sr (kg/m3-S) is the
mass generation rate by addition from the dispersed phase and other user-
defined sources.

In laminar flows, J.-,. is the diffusion flux of species i' due to the
concentration gradient in the i direction. FLUENT uses the dilute approximation,

under which Jr) can be written as

15

J..- = -pD ——'— (2.10)

where Dl',m (mzls) is the diffusion coefficient of species i' in the mixture, and is
also known as the mass diffusivity or diffusion constant [14]. Equation (2.10)
takes the same form as Fick's Law for a binary mixture.

In turbulent flows, the mass diffusion flux is computed by

J... = (an... JEEP (2.11)

where Sc, = 113—- is the turbulent Schmidt number with a default value of 0.7, m is

the turbulent eddy viscosity, and D. is the effective diffusion coefficient due to
turbulence.

The transport of enthalpy due to species diffusion

i'=l

is not negligible for many mixture flows, especially when the Lewis number

 

is far from unity (If air (Pr = i = 0.7 ) is the main fluid, Sc = ‘3 of water is 0.6, and
at

Le is 0.86 [14]). Here a = k

 

is the thermal diffusivity, v is the kinetic viscosity,
P

and D is the mass diffusivity.

2.2 Discretizatlon and Mesh Generation

16

The finite volume method (FVM) is used in FLUENT to convert the
governing equations to algebraic equations that can be solved numerically.

In the finite volume approach, the integral form of the governing equations
is solved. The discrete nature of the computational model is recognized at the
outset. This feature is shared in common with the finite element method (FEM).

Instead of solving the equations at points of the computational domain, as
is the case for the finite difference method (FDM), in FVM the physical domain is
divided into small volumes (or areas for 2-D cases) known as control volumes.
The physical law is satisfied over a finite region rather than only at a point as the
mesh is shrunk to zero. The governing equations are integrated for each control
volumes, which yields discrete equations that convert each quantity to a control-
volume basis. A distinctive character of FVM is that a ”balance" of a physical
quantity is always maintained for the control volume in the neighborhood of each

grid point.

Figure 5. Control Volume used for Discretization

17

FLUENT stores the discrete values of the physical quantity at the cell
centers (CO and C1 in Figure 5) of the control volumes as shown in Figure 5.
The convection terms in the governing equations require the face values, which
are interpolated from the center values. This is obtained by an upwind scheme.
"Upwinding' means that the face value is derived from the quantities in the
upstream cell relative to the normal velocity direction.

The first order upwind accuracy was employed in the present study as no
Significant difference was observed compared to the second order scheme in
predicting the outer scaled velocity profiles (outer scaling will described in section
3.2.1).

An advantage of the finite volume method (FVM) over the finite difference
method (FDM) is that the use of an unstructured mesh is allowed, since arbitrary
volumes can be utilized to subdivide the physical domain. This is very important
when the physical domain is highly irregular and complicated, where a structured
mesh would be difficult to construct. Also, since the integral equations are solved
directly in the physical domain, no coordinate transformation from physical to
computational domain is needed [13].

FLUENT 5.0 Is an unstructured solver. It uses internal data structures to
assign an order to the cells, faces and grid points in the mesh, and to maintain
contact between adjacent cells. Therefore, it does not require (i, j, k) indexing to
locate neighboring cells. The solver does not force an overall structure or
topology on the grid [12]. This provides a flexibility to use the grid topology that is

best for the specific problem of interest. In addition, for a complicated physical

18

domain, an unstructured mesh is much easier to generate compared to a
structured one. When an unstructured approach is employed, the unstructured
mesh generator (Gambit is the so-called pre-processor used by FLUENT to
generate the mesh.) is used to create the grids automatically, and the user can
concentrate his/her effort on defining the configuration of interest.

As a trade-off, unstructured grids require the connectivity information to be
stored cell-by-cell as opposed to block-by-block. Additional storage, therefore, is
necessary compared to the structured approach. From the standpoint of solver
efficiency, unstructured solvers are usually not as computationally efficient as
their structured counterparts [13]. Also when doing mesh generation, one needs
to keep in mind that numerical diffusion can be minimized when the flow is
aligned with the mesh [12].

In the present study, the mesh was refined from y" E 3 to y+ E 1 for the
wall adjacent cells of the thermally active surface without significant difference
observed in the prediction for the outer-scaled mean velocity (more about y+ and

outer scaling can be found in section 2.4.1 and 3.2.1, respectively).

2.3 Turbulence Models
The Reynolds number based on the jet nozzle width, Rew, for the current

study is

RC = pUjerw

w z1.3x104
,u

The flow is believed to be turbulent based on this number.

19

Turbulent fluid motion is an irregular condition of flow in which the various
quantities show a random variation with time and Space coordinates so that
statistically distinct average values can be discerned [15]. In the conventional

Reynolds decomposition, the randomly Changing flow variables are defined as
f = f+ f’ (2.12)
where f is the symbolic flow variable, f is the time-averaged quantity, and f '

is the fluctuating component. The time-averaged quantity, f , is defined as

f = Alt-[lwfdt (2.13)

where At should be large compared to the period of fluctuations of the
turbulence, but small compared to the time constant for any slow variations
associated with ordinary unsteady flow.

By definition, for symbolic flow variables f and g

7'20. fg'=0. E42: 71E=i+§ (2.14)

Note that 77" a: 0. In fact, the root mean square of the velocity fluctuations
is defined as the turbulence intensity.

Compared to laminar flows, turbulent flows are usually associated with
higher values of friction drag and pressure drop, and the diffusion rates of scalar
quantities are usually larger.

It is known that the unsteady Navier—Stokes equations are the governing
equations for turbulent flows in the continuum regime. However, solving those
equations directly, which is known as Direct Numerical Simulation (DNS) of

turbulent flows, is limited by current computational resources. Although

20

techniques like parallel processing are showing promising progress [16], it is still
computationally too expensive to Simulate directly in practical engineering
problems. Instead, the exact governing equations are manipulated to remove the
small scales. This results in a modified set of equations that are less expensive
to solve. The modification, however, brings in additional unknown variables.
Turbulence models, therefore, are needed to connect these variables to the
known quantities and achieve “closure” to the modified equation set.

With the above-mentioned Reynolds decomposition, the instantaneous
(exact) continuity and momentum equations can be time (Reynolds, or
ensemble) averaged. Dropping the over bar on the mean velocity, t7 , the so-
called Reynolds-Averaged Navier-Stokes (RANS) equations can be written in the

Cartesian tensor form as:

_+_ i)=0 (2.15)

 

. , 611.7 I
6x 6x, 6x, 6x, 3 6x, ax.

1

Compared to the instantaneous Navier-Stokes equations, now the velocities and

other variables represent time-averaged (mean) values.

 

The additional terms, - pui'u 1., , are called Reynolds stresses, which

represent the turbulence effects. These terms are the ones that need to be
modeled in order to close equation (2.16). The turbulence models employed in
the current study are briefly introduced. More details for the standard k-e model

are given as an example of how the turbulence models work.

21

It is worth noting that there does not exist single turbulence model that is
superior for all kinds of problems. The choice of turbulence model generally
depends on the physics of the specific problem, practical experience for a
specific class of problems, the available computational resources, and the
amount of time required for the simulation. Therefore, a basic understanding of

various models is needed for an appropriate choice.

2.3.1 Standard k-e Model

Two—equation models are the Simplest ”complete models” Of turbulence
[12], in which the turbulent velocity and length scales are independently
determined by solving two separate transport equations. The standard k-3 model
in FLUENT belongs to this turbulence model class. It is a semi-empirical model
based on the calculation of the turbulence kinetic energy (k) and its dissipation

rate (3), where the turbulent kinetic energy (k) is defined as

 

k=—u.u. (2.17)

l l

and 3 can be written as

3/2
a = C, 5r- (2.13)

where Cd is a drag coefficient approximately equal to 1, and L is the diameter of
an imaginary fluid mass packet [17].
Note that In FLUENT, a related quantity output, the turbulent intensity, is given by

(2.19)

22

where k is the turbulent kinetic energy, and V is a reference mean velocity
magnitude specified by the user in ”reference values” with a default value of
1m/s. This definition is inconsistent with the usual definition of turbulence
intensity where the numerator in (2.19) is the local mean fluctuation, and V is the
local mean velocity. When imposed as boundary conditions, however, the
turbulence intensity takes the definition in the usual sense.

The transport equation of k is derived from the complete momentum
equations (detailed procedure and be found in reference [17]). The transport
equation for 3, however, is obtained by physical reasoning and bears little
resemblance to its mathematically exact counterpart [12].

The assumption in the derivation of the model is that the flow is fully
turbulent, and that the effects of molecular viscosity are negligible. Therefore, the

standard k-3 model is valid only for fully turbulent flows.

Boussinesq Hypothesis and Modeling of Turbulent Viscosity

 

To model the Reynolds stresses, — pui'uj' , the Boussinesq hypothesis

 

I

6a

_ 'u _ $4.91!. _3. p“. __i.§ (220)
p“: j I“: 6x}. 6x, 3 ”tax. 7)- .

is adopted in k-3 models, which suggests that the apparent turbulent shear
stresses is related to the rate of mean strain through an apparent scalar, namely,
the turbulent or ”eddy“ viscosity. In equation (2.20), M is the eddy viscosity, k is
the turbulent kinetic energy, and Si, is the Kronecker delta.

The turbulent viscosity, m, is modeled in terms of k and 3:

23

k2
u. =pC.—g— (2.21)

where C,I is a constant. It is Clear that two additional transport equations for k and

3 are needed to close the system.

Transport Equations for the Standard k-c Model
The two transport equations for the turbulent kinetic energy. k, and its

dissipation rate, 3, are

Bk 6 , 6k
pfi=g[[p+g_]?ax_]+0k +Gb -p£—YM (2.22)
r‘ k r‘
DE 6 , 63 3 32
p5 = EKA +%]5x7]+0.. 7‘16. wean—Cap; (2.23)

In equation (2.22) and (2.23), G). represents the generation of turbulent kinetic
energy due to the mean velocity gradients. Gt, represents the generation of
turbulent kinetic energy due to buoyancy. YM represents the contribution of the
fluctuating dilation in compressible turbulence due to the overall dissipation rate.
C18, C26 and C3. are constants. ck and c. are the turbulent Prandtl numbers for k
and 3, respectively.

Furthermore, Gk, Ga and YM are modeled in terms of the known flow
variables so that the entire system is brought to a "closure":

Gk represents the production of turbulent kinetic enerQY. and is defined as

 

r Tau},
Gk =-puiuj 6x

 

(2.24)

Or, to be consistent with the Boussinesq hypothesis,

24

G, = 71,52 (2.25)
where S is the modulus of the mean rate-of-strain tensor, defined as

SE ZSS

(2.26)
with the mean strain rate S.) given by

6..
s, =1 95—4—1 (2.27)
, 2 3x, 6x,

G, is written as

126T

G = —
. 43.1,” ax!

(2.23)

where Pr. is the turbulent Prandtl number for energy with a default value 0.85 in

FLUENT. B is the coefficient of thermal expansion, defined as

p = —%[%)p (2.29)

YM takes into account the compressibility effects for high Mach number
flows:

Y” = 2p£M,2 (2.30)

where M( is the turbulent Mach number, defined as

M, = J: (2.31)
a

where a is the speed of sound, and k is the turbulent kinetic energy.

Model Constants

25

The constants In equations (2.21 ), (2.22) and (2.23) take the following
default values in FLUENT:

C,,=1.44, C2€=l.92, Cfl=0.09, a,=1.0, a=1.3

These values come from experiments with air and water for fundamental
turbulent shear flows, and have been found to work well for a wide range of wall-
bounded and free shear flows. Users are able to adjust these default values in

FLUENT.

Heat and Mass Transfer Modeling
Turbulent heat transfer is modeled using the concept Of Reynolds' analogy
to turbulent momentum transfer. The "modeled" energy equation has already

been given by equation (2.4).
a a 6 6T
abEhgrluJ/fi +p)]= EI-{kcfl E-ghflj. +uj(tij)efl]+sh

The term involving (1.1).“ represents the viscous heating. The deviatoric stress

tensor, (13).... is defined as

Bu,- au, 2 an,
(7.]. )4, = #47 [5; + a] - Eflefl 5:5.)- (2.32)

1'
Note the difference between equation (2.32) and (2.3). The effective thermal

conductivity k.“ and M are given by
k4, = k + k, (2.33)
and

#4)” = I” + I“: (234)

26

where k. is the turbulent thermal conductivity

_ cpfl!

kl
Pr,

(2.35)

with the default value of the turbulent Prandtl number, Pr. set to 0.85.
Similar treatment is applied to turbulent mass transfer. Note that in

equation (2.11), the turbulent Schmidt number takes the default value 0.7.

2.3.2 Realizable k-3 Model
There are two main deficiencies in traditional k-3 models.
First, by combining the Boussinesq hypothesis (2.19) and the definition of

eddy viscosity (2.20) and keeping 0,, as a constant, a negative value may be

obtained for the normal stress, “7, which is supposed to be positive. Meanwhile,
the so-called Schwarz inequality for shear stresses can be violated. This "non-
realizable" case happens when the strain is large enough [12]. The way to
ensure "realizable” is to make 0,, variable according to the mean flow (mean
deformation) and the turbulence (k, 3).

Second, the modeled transport equation of 3 in traditional k-3 models is
believed to be responsible for the unexpectedly poor prediction of the spreading
rate for axisymetric jets (the well-known round-jet anomaly).

The realizable k-3 model is intended to address the above mentioned
weaknesses by adopting:

- A new eddy-viscosity formula with a variable C.

27

- A new transport equation for 3 based on the dynamic equation of the

mean-square vorticity fluctuation

Modeling the Turbulent Viscosity
In the realizable k-3 model, the turbulent viscosity takes the same form as

that in the standard k-3 model (2.20).

k2
#1 =flp—
8

However, C,l is no longer a constant:

1

 

 

 

C = 2.36
,. U'k ( l
A. + A.
8
where
u’ a szs, +6562.) (2.37)

and

0,]. = Q”. — 23mm,

where '5; is the mean rate-of-rotation tensor viewed in the rotating reference

frame with the angular velocity (0].. The model constants A0 and A. are given by
A0 =4.04, A, =J6cos¢

where

S.S.S.
¢=—;-arccos(J6W), W: 03:1...

Cat
II
to
CI.)

 

28

Therefore, Cu is a function of the mean strain and rotation rates, the

angular velocity of the system rotation, and the turbulence fields (k, 3).

Transport Equations for the Realizable k-e Model

In the realizable k-3 model, the transport equations for k and 3 are

Bk 6 ,3: 6k

__=_ +——’——+G +G— -Y 2.38

pm 639“” 6,169.] I. ,, p3 M ( )
and

D3 6 ,u 63‘ 32 3

_=_ +_z__+ Ss— -————+C -CG 2.39

'01): and-Ky akjaxj] pC' pC2k+M ”k 3‘ b ( )
where

CI = max[0.43,—"—]

n+5

and

77 = Sk / 3‘
Note that the transport equation for k is the same as that for the standard k-3

model, while the one for 3 is quite different.

Model Constants

C,,=1.44, C,=1.9, e,=1.0, a'=1.2

Heat and Mass Transfer Modeling

29

The modeling of heat and mass transfer for the realizable k—3 model is the
same as that of the standard k-3 model.

It has been found that the realizable k-3 model is substantially better in
predicting flows like rotating homogeneous Shear flows, free flows including jets
and mixing layers. channel and boundary layer flows and separate flows. In

particular, it resolves the round-jet anomaly which occurs in standard k-3 model.

2.3.3 The Reynolds Stress Model (RMS)

The Reynolds Stress Models differ from the k-3 models in that the isotropic
eddy-viscosity assumption is abandoned, and turbulent shearing stress is not
assumed to be proportional to the rate of mean strain. That is, in a 2-D

incompressible flow,

_‘V,. 2.3:
pu Maya]:-

Reynolds Stress Models are more general than those based on the
Boussinesq assumption and are expected to give better predictions for flows with
sudden Change in the mean strain rate or with effects such as streamline
curvature or gradients in the Reynolds normal stresses [13]. Meanwhile, it is also
more computationally expensive compared to the two-equation models.

The Reynolds Stress Model (RSM) provided in FLUENT closes the
Reynolds-averaged Navier-Strokes equations by solving the transport equations
for the Reynolds stresses and an equation for the dissipation rate. The transport

equations for the Reynolds stresses are derived by taking the moments of the

30

exact momentum equations. For instance, in 2-D incompressible flow, the Ith
exact momentum equation is multiplied by the fluctuating velocity component u,’
and then the product of u.’ and the jth momentum equation are added. The

results are then Reynolds-averaged:

 

u,.'Nj +uj'N, = 0 (2.40)
where N. and N) are the ith and 1th components of the Navier-Stokes equation,
respectively. This kind of manipulation results in four additional equations in 20
flows, and seven additional equations in 3D flows.

Several terms in the exact transport equation (2.40) are unknown, and
turbulence models are required to close the equations. Details of modeling for
these terms can be found in [12].

In the Reynolds stress model, FLUENT solves the equation for turbulent
kinetic energy (k) that is essentially identical to that of the standard k-3 model in
the entire computational domain, but the k Obtained is used only for boundary
conditions. For cases other than the boundary condition, k is obtained by taking
the trace of the Reynolds stress tensor (2.17):

k=§mm

The dissipation tensor, 3.], is modeled as

2
3,] =35g.(p£+YM) (2.41)

where the scalar dissipation rate, 3, is computed the same way as it is for the

standard k-3 model.

31

In the Reynolds Stress Model, the turbulent viscosity and heat/mass

transfer are modeled similarly to the standard k-3 model.

2.4 Near Wall Treatment

The presence of walls Significantly affects turbulent flows, and is obviously
of interest for the study of wall jets. First of all, as with laminar flows, the velocity
field is affected by the no-slip condition at the wall surface. Meanwhile, from the
standpoint of turbulence, the flow is also affected in non-trivial ways due to the
existence of the wall. In the region very close to the wall, Viscous damping
reduces tangential velocity fluctuations, and the kinematic blocking reduces the
normal fluctuations. Farther away from the wall surface in the near-wall region,
however, the turbulence is rapidly augmented by the production of turbulent
kinetic energy due to the large gradients in mean velocity [12].

It is in the near wall region that flow variables undergo drastic Changes,
and the transport of momentum and other scalars is most vigorous. Walls are
also the main source of turbulence. Therefore, the near-wall modeling is critical
for the numerical simulations. Good resolution in the near wall region determines
the success of the wall-bounded turbulent flow predictions.

The turbulence models introduced in Chapter 2 are primarily valid for
turbulent core flows (the flow in the regions somewhat far from walls). Therefore,
near wall treatment is necessary to make these models work for wall-bounded

flows as well.

32

2.4.1 The Turbulence Boundary Layer and the Law of the Wall
Experiments show that the turbulence boundary layer can be subdivided

into two zones, namely, the inner region and the outer region as shown in Figure

 

  
 

 

6.
30
""'""—"' TYPICAL VELOCITY PROFILE
———— $61,, 1.. +515
25F- ,
_-—— u C y I
INNER REGTW
LAN OF THE HALL ZONE I
20>-
OUTER REGION ———>l
+ 151-—

   
  
   

I
L_ FULLY TURBULENLL:

VISCOUS
"’sueunrr:n"’l I (05mm zone

lor-

I I t l I l I l I I l I
I 2 5 IO 20 50 100 200 500 1000 2000 5000 10,000
I.

.7
Figure 6. Zones in the turbulence boundary layer (reproduced from [13])

 

 

rsurrtn 203911 I
l

 

The boundary layer equation can be written as [17]
fl '65 —§—5+i3[ 6‘7 “7“] (2.42)

Since the size of u' fluctuation corresponding to v' appears to depend on

the steepness of the average 6 profile [17], the notation of eddy shear stress

makes sense:

33

7 6a
- pu v = pa” 6? (2.43)

where the empirical function 3M is called momentum eddy diffusivity. Or, to keep
the consistency with the turbulent viscosity, it also has the meaning of kinetic

viscosity due to turbulence, and can be written as v.. Thus, in equation (2.42), the

molecular diffusion shear stress [1% is augmented by the time-averaged eddy

shear stress — 0W. The apparent Shear stress ram, therefore, can be given as

5&7 -,-, _
Tripp =#—_puv =p(v+€M)5y_

5y (2.44)
Note that 3.9. is a flow parameter, not a fluid property.

Substituting (2.43) into the boundary layer equation (2.42) yields
(2.45)

There exits a region close enough to the wall that the left-hand side of
equation (2.45) is sufficiently small, while the longitudinal pressure gradient g—xl:

is zero as in the case of uniform flow parallel to a flat plate. This region is also

characterized by an apparent shear stress that does not vary with y:

637 2'
M ———=—°— 2.46
(we )6y p ( )

in which To is the actual shear stress at the wall surface, that is, the apparent

shear stress, Tam, at y = 0 where the Reynolds stress — p 'v' vanishes. The

constant apparent shear stress assumption in equation (2.46) defines the inner
region in Figure 6.
To nondimensionalize equation (2.45) in the inner region, the so-called
friction velocity was introduced in the turbulence literature:
l/2
u, = (L0) (2.47)
p

With the following nondimensionalization know as the wall coordinates

 

 

. t7 . F
u = —— , v = —
“1' u!

x+ = "ur , y+ -_— yu, (2.43)
V V

 

(HEM-)5“ =1 (2.49)

which determines the velocity distribution in the inner region.
With equation (2.49), the inner region can be further subdivided into three
sublayers according to the relationship between the momentum eddy diffusivity

3M (V() and v as shown in Figure 6:
(1) The viscous sublayer, where v >> 3M.
(2) The buffer zone, where v and 3... are of the same order of magnitude.
(3) The fully turbulent sublayer (or the log-law zone), where 3.4 >> v.
In the viscous sublayer, the flow is almost laminar-like. Dropping the term

v/eM in equation (2.49), we learn that velocity has a linear profile:

u: = y: (2.50)

35

In the fully turbulent sublayer, the effects of turbulence are dominant. The

 

equation (2.49) reduces to
523:,“ :1 (2.51)
V

and the velocity profile takes the form [17]
u+ =Alny+ +B (2.52)
where A and B are two empirical constants. Equation (2.52) is known as the law

of the wall, which is valid in the log-law zone in Figure 6.

2.4.2 Near Wall Modeling In FLUENT
There are two approaches available in FLUENT to model the near wall

region: the wall functions and the two-layer zonal model.

Wall Functions

In this approach, the viscous-affected inner region (the viscous sublayer
and buffer zone in Figure 6) is not resolved. Instead, the semi-empirical formulas
called ”wall functions“ are adopted to ”bridge” the inner region between the wall
and the fully turbulent sublayer.

First, the wall units u' and y' were introduced:

Chg/5‘7 . pC,:/‘k’!'5y
= , y :—
To/p fl

u

 

(2.53)

where k is the turbulent kinetic energy.
When y' < 11.225, FLUENT applies the laminar stress-strain relationship

for the velocity profile:

36

u' = y' (2.54)

When y' > 11.225, the law of the wall is employed as
u’ = £1110: ‘) (2.55)

where 1: is the von Karman's constant (= 0.42), E (= 9.81 )is an empirical
constant. (The logarithmic law for mean velocity (2.55) is known to be valid for y'
> 30 ~ 60.) It is worth noting that the law of the wall employed in FLUENT (2.55)
is based on the wall units y' and u’ rather than the wall coordinates y+ and u‘.

The Reynolds' analogy between momentum and energy transport gives a
similar logarithmic law for mean temperature [12].

If the wall function approach is chosen in k-3 models or RSM, the k
equation is solved in the entire domain with the boundary condition

1%.

6n 0 (2.56)

imposed at the wall surface, where n is the local coordinate normal to the wall.

In the wall-adjacent cells, the production of kinetic energy. Gk, is computed

 

from [12]
G, s to Q = r, ———f°l— (2.57)
Kpcildk/Zy
And instead of solving the 3 equation, 3 is computed from
C%k%
3 = ,. (2.58)
xy

in the wall-adjacent cells [12].

37

The wall function approach gives acceptable predictions for most wall-
bounded flows with high Reynolds number. However, it is inadequate for low
Reynolds number flows or when near wall effects are pervasive, in which the
underlying hypotheses for the wall functions cease to be valid (e.g., flow though a

small gap as is the case for the present study).

Two-Layer Zonal Model

In this approach, the turbulence models are modified in order to resolve
the viscosity-affected region all the way to the viscous sublayer. The entire
domain is subdivided into two regions: the viscosity-affected region and the fully-
turbulent region. The demarcation of the two regions is determined by the so-
called turbulent Reynolds number, Rey, based on the wall-distance. The

definition of the turbulent Reynolds number is

E lax/Ty
.V ’1

Re (2.59)

where y is the normal distance from the wall at the cell centers.
In the viscosity-affected region where Rey < 200, a one-equation model is
employed [12]. The same momentum and k equations described earlier are

solved. However, the 3 field is computed from

 

3/2
3 = k1 (2.60)
and the turbulent viscosity, m, is computed by
A. = 726‘. J31. (2.61)

38

where the length scale formulas are given as

1 Rev

. =0Iy l-exp - A. (2.62)
l — 1 Re-" 263
p "Cly "CX _ A” ( ' )

The constants in the above equations are

 

 

c, = xC‘3/‘, A, =2c,, A = 70 (2.64)

y a
In the fully turbulent region where Rey > 200, the k-3 models and RSM are
turned on.
Since the near wall effects are considered crucial for the present study,
the two-layer-zonal model is believed to be more suitable for the near treatment

compared to the wall function approach.

2.5 Summary

FLUENT provides comprehensive models in single-phase, single-species
simulations. No existing multi-phase model, however, was found suitable to
capture the condensation/evaporation phenomenon in the present study.
Nevertheless, in FLUENT certain portions of the solver are open to the users by
User-Defined source terms in the governing equations. It is one of the efforts of
the current study to model the phase change effect by implementation of

appropriate User Defined Functions (UDF'S).

39

Chapter 3

SIMULATION WITHOUT PHASE CHANGE

This chapter is focused on the velocity and temperature simulations
without phase change, the so-called "dry air" condition, which happens when the
temperature is higher than the dew point of water vapor everywhere in the
computational domain. Simulation results are compared with experimental and

published data.

3.1 Solution Procedure
The basic steps to take in order to solve a problem with FLUENT are

given as the following:

- Create model geometry and grid with the pre-processor Gambit

. Choose appropriate solver for 2D or 3D modeling when starting the solver
FLUENT 5

- Import the grid into the solver

- Select the solver formulation: steady or unsteady, explicit or implicit, etc.

. Choose the basic equations to be solved, and additional models need: i.e.
multi-species

- Specify material properties

- Impose boundary conditions

- Adjust the solution control parameters

40

- Initialize the flow field
. Calculate a solution
- Save and examine the results
I If necessary, refine the grid or consider revisions to the numerical or physical
model

The geometry and mesh generation can be completed by the pre-
processor Gambit. The mesh file is then imported into the solver (FLUENT 5 is
actually the solver and post-processor.) The solver also accepts l-DEAS
universal files, NASTRAN files, PATRAN neutral files, and ANSYS prep7 files.
FLUENT also has a fairly strong post-processing function to analyze and

examine the computation results.

3.1.1 Pre-processing

Strictly Speaking, the experimental facility is not really a 2D platform. The
inlet and outlet are not of the same length in the z direction of Figure 1, and
neither of them goes across all the way through the z direction of the facility.
However, since the measurements are taken mainly at the center portion of the z
direction of the facility, it is believed that the third direction effect is not significant,
and the whole geometry can be treated as a ZD model in the x, y plane.

As mentioned in Chapter 1, the curvature of the contraction is a special
design with the aim that the velocity profile distributes evenly along the y
direction at the jet nozzle after the flow passes through the contraction [1]. This

curve is difficult to create directly with the pre-processor. Therefore, the geometry

41

of the computation was imported into Gambit as an IGES file created by
AutoCAD as shown in Figure 7.

PressureNelocity Inlet

— Symmetry Plane

Contraction H

/JetNozzte (w-2an)

 

‘ —- ThennaIy Active Plate
(with Constant Temperature and a total Length of 25 am)

 

i

Receiver 80:

——-— Symmetry Plane

 

 

 

PressureOmtet

Figure 7. Geometry and Boundary Types of the Simulation

Note that the symmetry of the geometry is taken into account to reduce

the computational effort.

The mesh was generated with Gambit with quadrilateral cells. It is well
structured and highly clustered in the thermally active test section (see Figure 1),
especially in the region near the thermally active plate (See details in the near

wall treatment section). The mesh file was then imported into the solver.

3.1.2 Model Control

- Solver: Segregated (Governing equations are solved sequentially. i.e.,

segregated from one another. Because the governing equations are non-

42

linear (and coupled), several iterations of the solution loop must be performed
before a converged solution is obtained [12].)

- Formulation: Implicit

- Space: 2D

I Time: Steady

- Energy: Enabled for non-isothermal simulation

- Viscous (Turbulence Model): k-3 or RSM

3.1.3 Near Wall Treatment
AS was discussed in Chapter 2, given the limitations of the wall function
approach and specific needs of the near wall resolution, the two-layer-zonal

model is adopted as the near wall treatment for the current study.

3.1.4 Unlt System

lntemally, FLUENT stores all the parameters and performs all the
calculations in SI units. Nevertheless, one can work with different unit systems,
even inconsistent units because a correct set of conversion factors is built in the
solver, which converts the units one wants to use into the standard SI unit
system. The restriction is that boundary profiles, source terms, custom field
functions, and user-defined functions must be in the SI units.

The length unit of the AutoCAD design was in inches, and was converted
into meters before the calculation. This is for the sake of consistency since all the

other parameters work in the SI system.

43

3.1.5 Material Properties

Because no phase change effect is taken into account at the current
stage, the fluid is considered as ”dry air”. The default physical properties of air
stored in the FLUENT database (listed in Table 1) were employed.

Table 1. Default Values of Physical Properties of Air in FLUENT

 

 

 

 

 

 

Density p 1.225 (kg/m3)
Specific Heath 1006.43 (J/kg-K)
Thermal Conductivity k 0.0242 (W/m-K)
Viscosity l1 1.7894e-05 (kg/m-K)
Molecular Weight 28.966 (kg/kmol)

 

 

Also, since the experimental data were measured mostly in the range Of
common room temperature, it is assumed that physical properties of air do not
vary significantly with temperature. It is assumed that all the physical properties
are constants and the air is incompressible. This assumption makes the energy
equation uncoupled from the continuity and momentum equations. Therefore, it is
not necessary to solve the flow equations (continuity, momentum and turbulence)
and the energy equation simultaneously. Instead, the temperature field can be
calculated based upon simulation of the velocity field (the solution of the flow

equafions)

 

3.1.6 Operating Conditions
Operating pressure, pop, is important for incompressible ideal gas flows
because it directly determines the density: the incompressible ideal gas law

computes density as

=ffl

pRT

As for compressible flows,

= pop +P
RT

 

where p is the gauge pressure, R is the universal gas constant, and T is the
temperature in Kelvin.

The operating pressure for the current simulation is set to be 101325 Pa
(1atm). The reference pressure location (where the gauge pressure is always
zero) can also be specified. But when pressure boundaries are involved, the
reference pressure location is ignored since it is no longer needed [12].

FLUENT uses gauge pressure in calculation. When absolute pressure,

pans, in needed, it is obtained by

pubs =pop +p

3.1.7 Boundary Conditions
As mentioned earlier, the driving force for the system in Figure 6 is the
pressure differential between the inlet and outlet. Pressure inlet and outlet

boundary conditions are supposed to be imposed for them. However, it was

45

found that the simulation results were almost the same when the boundary
condition at the inlet was specified as a velocity condition with a constant value.

The thermal boundary conditions that need to be specified are the
temperature at the inlet, backflow temperature at the pressure outlet (backflow is
predicted numerically at the pressure outlet as seen in Figure 12), and the
constant temperature on the surface of the thermally active plate.

FLUENT defaults the boundary condition to be adiabatic wall, where the
heat flux equals zero and all the velocity components are zeros. So apart from
the two symmetry planes in Figure 6, the only boundary conditions needed for
the simulation are:

Inlet:
I Total gauge pressure: 0 Pa; or constant velocity: 1.6896 m/s

(according to the experimental measurement)

1
. . p, = P. + -lV|’
Note that the Bernoulli's equatIon 2

holds at the pressure
inlet boundary condition for incompressible flows, where p. is the total pressure;
p. is the static pressure; and M is the velocity magnitude.

- Turbulence intensity: 2.5% (as was measured in experiment);

Hydraulic diameter: 0.284 m (calculated from the real geometry).

- Total temperature: ~ 298.15 K (room temperature)

Total temperature is the temperature at the thermodynamic state that
would exist if the fluid were brought to zero velocity. Static temperature is the

temperature that is measured moving with the fluid. For incompressible flows, the

total temperature is equal to the static temperature.

46

Outlet:
- Gauge pressure: -60 Pa
- Backflow turbulence intensity: ~30% (estimated)
Backflow hydraulic diameter: 0.42 m

The value of backflow turbulence intensity is not known from experimental
measurement. Nevertheless, it was found from calculations that the simulation
results in the thermally active section are not sensitive to this value since this
region is far away from the pressure outlet.

- Backflow total temperature: same as the temperature of the inlet
Thermally Active Plate:

- Wall boundary with a constant temperature: ~ 278.15 K (adjustable)

It is worth noting that, from the standpoint of comparison, it is the types of
boundary conditions that are of importance, but not the magnitudes. Since the
simulation results are usually normalized before they compared to other data,
many input values for the boundary conditions do not effect the results very much

after normalization.

3.1.8 Solution Control

The segregated solver makes it possible to solve the governing equations
separately. Since the physical properties are assumed to be constants as was
discussed in the section of material properties, the following strategy was taken

in order to save some running time.

47

- Perform velocity simulation by selecting flow equations and while keeping the
energy equation deselected
- Perform temperature simulation by turning on the energy equation and
deselecting the flow equations after the velocity solution converged
The default under-relaxation factors are kept for the k-3 model. However, it
was found difficult to get converged solutions with the Reynolds stress model
(RSM). Therefore, more conservative under-relaxation factors are employed for
RSM.
Most of the simulations are performed with the first order upwind
discretization since no major improvement was achieved with the second order
discretization when the velocity data in the wall jet region are compared to the

experimental data.

3.1.9 Convergence Judgement

There Is no universal method for judging convergence. The residual
definitions are employed in FLUENT. The residual of a general variable
computed by FLUENT‘S segregated solver is the imbalance of the variable
summed over all the computational cells. It is scaled using a scaling factor
representative of the flow rate of the variable through the domain since it is
difficult to judge convergence by examining the non-scaled residuals. The
solution is considered converged when the scaled residuals go below a certain

set of default values that can be changed by the user. The computation will be

48

stopped when the residuals of the governing equations are less than the
corresponding default values.

It was found that the default values of FLUENT for the convergence
criterion are not satisfactory in terms of experimental verification. The
convergence judgement employed here was not by the value of the residuals
themselves, but by their behavior. It was ensured that the residuals continue to
decrease (or remain low) for several iterations (say 50 or more) before

concluding that the solution has converged.

3.2 Results and Discussions of the Isothermal Simulation

The simulation without phase change is first performed for an isothermal
condition, which means that the temperature on the plate temperature equals
that of the jet (inlet room temperature). The velocity field and turbulence variables
tested from the experimental facility and some other published data are utilized
for comparison with the numerical simulation.

The velocity distribution of the flow field is of great interest in the wall jet
study. The data in the downstream self-similar region appear to be more
interesting to most researchers as introduced in Chapter 1. For the study of the
defogger flow, however, more attention needs to be paid in the developing region
since on the interior windshield surface of a vehicle it is in the developing region
that the defogging effects need to be applied.

The experimental facility (depicted in Figure 7) to be simulated has a jet

nozzle with a slot width of 2 cm, and the thermally active plate spans 25 cm

49

along the streamwise direction of the jet. Six streamwise locations (xlw = 0, 1.59,
3.18, 4.45, 7.62 and 10.8) were chosen to plot velocity profiles from the
simulation in order to compare with the experimental measurements. Figure 8
Shows the development of the velocity profile simulated by the realizable k-3
model. A plot of over-lapped velocity profiles at the six tested locations is given in
Figure 9. The abscissa of the plot represents the y distance from the plate
(surface of the wall). The ordinate is the magnitude of velocity.

Some visualization of the simulation results is provided to give some basic
impressions of the flow field. Figure 10 and Figure 11 are the contour of velocity
magnitude and the static pressure, respectively. Figure 12 shows the
visualization Of velocity vectors in different regions of the flow. Figure 13 depicts
the viscosity-affected sublayer at various locations for the two-layer-zonal model
for the near wall treatment, where the turbulent Reynolds number, Rey, is less

than 200.

3.2.1 Velocity Profiles

Data In the Outer Scaling
As in most wall jet literature, the velocity profiles were first normalized with
outer scaling. That is, the velocity magnitude is normalized by the maximum

value of the profile, and the y distance is normalized by 62 described in Chapter
1. The numerical simulation was performed with the standard k-3 (ske in the

Figures) models, the realizable k-3 (rke in the Figures), and the Reynolds stress

50

model (rsrn in the Figures). The comparison between the CFD simulation and the
experimental data is given by Figure 14 through 19.

It can be observed that the velocity profile predictions of realizable k-3
model and RSM are very similar. According to the test data, they both give a
better velocity prediction at the jet nozzle (xlw = 0) and the region close to the
nozzle (xlw = 1.59) compared to the standard k-3 model. This means that the
realizable k-3 model and RSM better simulate the effect of the upstream
contraction. As mentioned earlier, the design purpose of the contraction is to
provide an evenly distributed velocity profile at the jet nozzle. This is proved by
the experimental measurement at xlw = 0. The standard k-3 model appears to
over-predict the friction effect of the curved surface of the contraction so that the
velocity magnitude close to the side of the curved surface is smaller than that
measured.

On the other hand, however, the standard k-3 model gives better
predictions as the jet develops along the wall. Especially, for the location where
the velocity profile takes the maximum value (y... in Figure 2), the standard k-3
model gives much better predictions than the other two models. This can be
explained by the fact that the standard k-3 model was originally developed from
experiments with air and water for fundamental turbulent Shear flows [12]. The
empirical model constants of the standard k-3 model have been found to be
appropriate for plane jets and plane shear layers [17] and a wide range of wall-

bounded and free shear flows [12].

51

The comparison of the development of velocity profiles suggests that the
realizable k-3 model and RSM appear to under-predict the diffusion of the
turbulence dominated core (free Shear layer) as the jet develops along the wall
surface. This can be clearly seen from the development of the turbulent viscosity
ratio, (1.41, predicted by the two different models, as shown in Figure 20 and
Figure 21. The profile of the ratio of turbulent viscosity to molecular viscosity
suggests the development of turbulence. From Figure 20 and 21, it is observed
that the standard k—3 model predicts a higher value for the turbulent viscosity at
turbulence-dominated core compared to the other two models. This may result in

different predictions of the mass transfer according to equation (2.11).

Data In Inner Scaling

The wall coordinates (2.48) scheme introduced in Chapter 2 is also
referred to as the inner scaling. Figure 22 shows the velocity comparison of the
CFD simulation to experimental data at the initial stage (xlw = 1.59) and the
relatively well developed one (xlw = 10.8) with inner scaling.

Wygnanski et al [5] argued that, in the well-developed region far away

from the jet nozzle, the constant A in the log-law
u+ =Alog(y+)+B (3.1)
does not vary with the jet Reynolds number Rew, while B does. Their

measurements showed that A always equals 5.5, while B varies between 4.9 to

9.2 as Rew goes from 3700 to 19300. Given the jet Reynolds number of the

52

present study (13000 as given in section 2.3), B can be linearly interpolated as a
very rough approximation. Thus, equation (3.1) can be written as [18]

u+ = 5.510g(y+) + 8.83 (3.2)
The plot of this equation is also given in Figure 22 as a reference.

Worth mentioning is that the simulation data of u“ in Figure 22 is strictly

according to the wall coordinates definition (2.48). But the u" in the experimental

data is actually the velocity magnitude (U = W ) divided by the friction
velocity (u.) since the technique employed in the measurements in not capable Of
detecting the individual velocity components (u, V). Although v is considered
negligible compared to u (which makes U E u), the difference between U and u
may still account for the discrepancy between the numerical and experimental
data in Figure 22.

Figure 23 compares the inner-scaled distribution of u” predicted by the
standard k-3 model to the published data by Eriksson et al [7]. For reference, the
overall inner-scaled predictions by the standard k-3 model are given in Figure 24.

The experimental data of the u+ distributions of Eriksson et al [7] and Hoke
et al [18] are quite different in the region close to the jet nozzle. Nevertheless, in
the region farther away from the nozzle, the data by Hoke et al [18] appear to
agree with the equation (3.1) by Wygnanski et al [5] (xlw = 10.8 in Figure 22).
The CFD simulations appear to match well with the experimental data by
Eriksson et al [7] (xlw =10 in Figure 23).

The above comparisons suggest that the wall jet can be considered

relatively well-developed farther away from the nozzle when the inner-scaled

53

data tend to become self-similar. In the near-nozzle region, however, the wall jet
is still in the developing state, the inner-scaled data are dependent on the
upstream geometry situations for different experimental conditions.

Since the same near wall approach (two-layer zonal model) was adopted
for all the three different turbulence models, it is not surprising to see that the
predictions of three different turbulence models do not differ from each other very
much in inner scaling (Figure 22). All the CFD predictions give lower values of u”
compared to the experimental data by Hoke et al.

In the relatively well-developed region (xlw =10 in Figure 23), the
predictions match with the data of Eriksson et al fairly well, especially when y“ <
100. That is, the two-layer zonal model gives better predictions in the well-
developed region of the wall jet.

In fact, the numerical simulations give similar predictions of u+ in both the
near-nozzle developing region and the relatively well-developed region for y“ less
than 100 as shown in Figure 24. This suggests that the two-layer zonal model
tends to predict the distribution of u+ in both regions as if the wall jet is well-
developed. The Similar predictions can be found in the predictions of the other

two turbulence models. The compiled Realizable k-3 (rke) and Reynolds Stress

Model (RSM) predictions of the inner-scaled velocity profiles are given in Figure
25.

3.2.2 Turbulence Fluctuations

Data in Outer Scaling
In experimental studies, the so-called Streamwise Normal Stress (SNS) is

used to account for the velocity fluctuation due to turbulence. It is defined as

—,2

SNS =100l‘j‘2

max

(3.3)

 

where 3 is the fluctuation of the velocity, and Um...x is the local maximum velocity
of the velocity profile at the location x/w.

In the numerical Simulations, the fluctuations are assumed to be isotropic
in the k-3 models. The magnitude of 57 can be obtained by using equation (2.19)
and setting the reference velocity to be 1 (m/s). For comparison, I? is obtained
the same way in the Reynolds stress model. Comparisons at the developing
region (xlw =1.59) and developed region (xlw =10.8) are plotted in Figure 26.

AS can be observed, by using equation (2.19) to get the fluctuations, none
Of the models provides satisfactory agreement with the experimental data. The

standard k-3 model gives worse predictions in the developing region than the
other two models.

It is worth noting, nevertheless, the individual stress components ( W,
v "v ,u "v and w’w 'for 2D flows) can be obtained by the Reynolds stress model
(RSM). So the SNS can be obtained directly by the u'u 'stress in RSM instead of
using equation (2.19) for 3.

Figure 27 shows that the SNS calculated directly from the u'u 'stress from

the RSM matches the experimental data much better than the ones predicted by

the isotropic k-3 models (see Figure 26).

55

For reference, the compiled simulations of outer-scaled SNS are included

in Appendix A.

Data in Inner Scaling

In the inner scaling, the velocity fluctuation is non-dimensionalized as

72
u” = “ (3.4)

u

f

 

where ut is the friction velocity defined in equation (2.47).

Figure 28 shows the inner-scaled comparison of predicted u" distributions
and experimental data (Exp in figure) in the developing region (xlw = 1.59) and
the developed region (xlw =10.8). The fluctuation boundary conditions at the
pressure inlet are the same for all three turbulence models: turbulence intensity

equals 2.5%. The u'+ distribution by the Reynolds stress model (uu_rsm in figure)

was calculated directly from the 'u' stress.

None of the models produce predictions that are in very good agreement
with the experimental data. The agreement between the standard k-3 model (ske
in figure 26) predictions and the experimental data is especially poor in the
developing region. In contrast to the inner scaled comparison of mean velocity,
the experimental data are not consistently larger than the predictions at xlw =
1.59. Thus, It is difficult to explain this solely on the basis of the experimental

difficulty in calculating the friction velocity ut (or the surface Shear stress To).

Eriksson et al [7] obtained experimental data for 07+, 7 and 07+. Their

data showed that the fluctuations appear to become self-similar when xlw is

56

larger than 40, which suggests that the fluctuations of the current study are
developing in the entire length of the plate (x/w s 12.5). Also, their upstream
conditions are different from the present study. Since the fluctuations in the
developing region are believed to be highly dependent upon the upstream
conditions (especially the turbulence intensity as the boundary condition at the
inlet), the value of the direct comparison between the current simulations and the

published data is somewhat questionable. Figure 29, 30 and 31 are the
comparison of the RSM simulation and Eriksson et al's measurement of 07+,
17+ and 337, respectively. The available streamwise locations for the published
data are xlw = 0, 5 and 10 for E“ (Figure 29), and xlw = 5 and 10 for 7+ and

077+ (Figure 30 and 31).
For reference, the Simulations of inner-scaled fluctuations with RSM are

included in Appendix B.

3.3 Results and Discussions of the Non-Isothermal Simulation

The study of the heat transfer and temperature distribution in a wall jet is
rarely found in the literature. It is very important, nevertheless, in terms of the
defogger research because the temperature distribution directly affects the
concentration distribution of water vapor, and in turn, influences the mass
exchange of water on the interior surface of the windshield (the wall surface).

With the assumption of constant physical properties, the velocity
simulation is independent of temperature field. However, the success of the

temperature field simulation depends on how well the velocity field is predicted.

57

Data in Outer Scaling

In the outer scaling, the temperature is normalized as

7': T4: (3.5)
T T

inlet — s

 

where T is the local temperature; T,. is the temperature held constant on the wall
surface; and TI"... is the temperature at the inlet of the whole system, which is
believed to equal the temperature in the main stream of the jet. For the present
simulation, Tm... = 297.15 K (24 °C), and T.3 = 277.15 K (4 °C).

The three turbulence models, namely, the standard k—3 model (ske), the
realizable k-3 model (rke), and the Reynolds Stress Model (rsm), predict virtually
identical temperature profiles in the outer scaling, as can be seen in the xlw =
10.8 profile in Figure 32. Actually, the predictions of the three turbulence models
differ the most at xlw =10.8. Even so, the difference is still difficult to detect. Also
included in Figure 32 are the comparison Of the experimental data and the
standard k—3 prediction at x/w =1.59, and xlw = 4.45.

The temperature predictions have better agreement with the experimental
data in the relatively well-developed region (xlw =10.8) compared to that in the
developing region (xlw =4.45). This appears to be consistent with the velocity
predictions shown in Figure 15: At xlw = 4.45, the experimental velocity profile
has a larger gradient than the predicted one. And so is the case for the

temperature profile at the same streamwise location in Figure 32.

58

For reference, the simulations of outer-scaled temperature profiles with

the standard k-3 model are included in Appendix C.

Data in Inner Scaling
The energy equation of a flat plate turbulent boundary layer is given as
[17]

air-Ni = 3[(a + 3,, 3i] (3.6)
0x 6y 0y 6y
where 0 is the thermal diffusivity; and the empirical function an is known as the
thermal eddy diffusivity.
Analogous to the constant rapp assumption made earlier in Chapter 2, it is
assumed that the left hand side of equation (3.6) is negligible sufficiently close to

the wall, and the apparent heat flux q"app does not depend on y. Thus, equation

(3.6) becomes

 

(a+3',,)%::=[(a+3,,)%] (3.7)
That is,
(a + 3,, )6;- = [:5 (3.3)

Introducing the wall coordinates in (2.48), equation (3.8) becomes

pCPur_a_—_ =—1__
-q: W“ T") a/v+a../v (3'9)

Defining the temperature in wall coordinates as

59

 

T+(x+,y+)=(To—T)£-CC—u' (3.10)
40
equation (3.9) becomes
. dy“
T = 3.
Jo l/Pr+(l/Pr, )(3,, /v) ( 11)

where Pn is the turbulence Prandtl number, defined as

Pr, = f!— (3.12)

8}!
Equation (3.12) says that the temperature profile in the constant q"app region is
governed by the Prandtl number, the turbulent Prandtl number, and, via the
turbulence viscosity ratio 31)./V, the velocity distribution in the same region.
With similar assumptions and derivation for velocity in Chapter 2, the

temperature profile can be obtained from (3.11) [17]:

Pry+a(}’+ < yz‘SL)

 

T+ = + Pr + + + (3.13)
PryCSL +41“ }: 2(y >yCSL)
yCSL

where flu is the dimensionless thickness of a conduction sublayer in which the

molecular mechanism outweighs the eddy transport of heat [17]. Good
agreement with temperature measurements can be achieved when the empirical

constants take the following values [17]:
Pn E 0.9, it E 0.41, yESLE 13.2 (3_14)

Note that equation (3.13) is valid only in the region where ram and q'o are both

constants.

60

As can be seen in equation (3.14), in the fully turbulent region where q"app
and raw are still constant, the temperature profile "I" is analytically the same as
the law of the wall for the velocity profile.

The log-law form of the temperature profile in equation (3.13) for air (Pr =
0.7) is

T“ =2.1951ny+ +3.53 (3.15)
if the empirical constants take the values in equation (3.14).

Figure 33 is the comparison of inner-scaled temperature profiles. The plot
of equation (3.15) is also included in Figure 33 for reference. Note that the
default value of Pr( in FLUENT is 0.85 instead of 0.9 in equation (3.15).

Since the numerical prediction matches with the log-law form
quantitatively, it is believed that the relatively poor agreement between the
experimental and the numerical data is due to the difficulty in measuring the
velocity and temperature close enough to the wall so that the shear stress and
heat flux at the wall surface (needed in the wall coordinates calculation) can be
calculated more accurately.

For reference, the simulation of the inner-scaled temperature profile is
included in Appendix D. Also include is the heat transfer coefficient on the wall

surface, defined in FLUENT as

 

h = ‘10 (3.16)

where To is the temperature on the wall surface; Tm, is the reference temperature

specified by the user as the inlet temperature (main stream temperature in the

61

jet). Note that q'o is positive when the direction of the heat flux is into the

computational domain (fluid).

3.4 Summary
For the velocity simulation, the three turbulence models employed
(standard k-3, realizable k-3, and Reynolds Stress model) give fair predictions for

the velocity field compared to the experimental measurements with the two-layer
zonal approach for the near wall treatment. The standard k-3 model gives better
outer-scaled predictions in the relatively well-developed region far from the jet
nozzle. The Reynolds stress model, on the other hand, gives better predictions
for the turbulent fluctuations.

For the temperature simulation, the three turbulence models give virtually
identical predictions. The simulation results are in fair agreement with
experimental measurements in the outer scaling form, and appear to be
reasonable compared to the inner-scaled analytical temperature distribution in
the log-law form. The reason for the discrepancy between the simulation and the

inner-scaled experimental measurements is yet to be found out.

62

 

 

 

 

 

E41104“ ”HI"

3...... ’ - - - - ' - ‘ - ' - "””””“'"‘““"’lllllll

:"°°"°° ' ”""""""lllllllllll

#39404“

3...... ' ' ' ‘ ' ‘ ' ‘ ' ' ‘ """‘“""“""‘””“lllllllll

:6190-06 b—l

Velocity V0110" cm BY V°'°°"V mm (M) FLUENT 5.0 (2d. $03.93

 

 

 

Figure 8. Development of Velocity Profile

 

 

 

 

 

 

 

 

1.20.+01 -
1.00.401 - f/ """"
i Til,
3.00.400 '
g .
Velocity coo-+00 ~ “1‘;
Magnitude . l ‘,
(m/s) L i".
4.00.+00 ( (.
2.00.400 1L
o.ooo+00 - -, --’”'.-- .--, -
0 0.01 0.02 0.03 0.04 0.05 0.03
Position (m)
Velocity Magnitude AugOQ, 1999
FLUENT 5.0 (2d. segregated, I'ke)

 

 

 

Figure 9. Plot of Velocity Magnitude

63

1.04a001
‘71-'33 9.40””
0.363000

7.31.900

0.00000

‘ Coat s 0‘ Valoc' HWa [ah I
0" 6’ ] HUI" 5.0 (ad. segregated. rh]

Figure 10. Contours of Velocity Magnitude

 

 

6010-02

5 -7 .50“
I -1 5204-01
4290-0-01
6.050001

4.020001

 

. 4.50am

-5.35o+01
4.11.061
T‘

-6 3804-01

 

-7 540-001

Canton of Static Preasue (pascal) Dec1998
BC‘s: P(lrlet_tolial) = OPa. l(lnlet)- - 10%, Wallet _statlc)= -60Pa FLUENT 5 0 (26. segege1ted. rite)

Figure 11. Contours of Static Pressure

 

 

1.010101

 

 

IIIJIIIIIIIIIIIIII

 

l_ 9.110100
6.100100
1.090100
6.070200
5.oso+oo
4.050100
3.04uoo
2.02am '
101-+00

6.790-05

 

 

 

 

Velocity Vectors Colored By Velocity Magitude (m/s) 20, 1999

A09
FLUENT 5.0 (2d, seg'egated, ke)

 

 

1.010401

H 9.110400

_ 3.100100

7.090“

6.070400

5.060100

4.0504“

3.040100

2.020%

1.01000

 

LJIIIIIITIIIIITJ

6190-05

/

/

/,///

\

'\
\
'\
\
\
\

 

 

 

Velocity Vectors Colored By Velocity Magnitude (m/s)

Aug 20, 1999
FLUENT 5.0 (2d, sewegated, ke)

 

 

Figure 12. Visulization of Velocity Vectors

65

 

 

1.400402

1.200402

 

3),; 3.00am

 

 

 

Viscosity-Affected Stblayer In the Wel Jet Area Aug 09, 1999
FLUENT 5.0 (2d, segegated, rke)

 

Figure 13. Viscosity-Affected Sublayer (Rey <200)

66

 

ComperleonofVeloclty Profileetxlw=o I H I
(Outer Scaling)

 

 

 

1.5 2

 

Comparison of Velocity Profile at xlw - 1.59
(Outer Scaling)

 

  
   
  
 
   

A" lew‘;1h.59—.Ske.
‘ CI xlw=1.59,test

 

 

 

Com perleon of Velocity Profile at xlw-3.18
(Outer Scaling)

 

. 421672313. ‘ske ~

 

 

Figure 14. Comparison of Velocity Profiles (Outer Scaling, ske, A)

67

Comparison otVeIocIty Profileatxlw I 4.45
(Outer Scaling)

     

 

r——"xlw=4:45, ske
El xlw=4.45. test

 

 

i 0.6
5 0.4
0.2
o - .
o 0.5 1 1.5 2
Y’s:

c.5223... of 17.1.31. ProfiIO at tutu-7.02

 

 

 

(Ouhr Scaling)
1
0 8 ‘_xM =7.62. ske
. ' D xlw =7.62. test
j 0.6 '
i
0.2
0 .
0 0.5 1 1.5 2
We,
Comparieon ofVeIoclty Profile at xlw I 10.8
(Outer Scaling)
1

 

    
  

x/w=10.8,ske
- D x/w=10.8.test

 

 

 

0.6
l 30.4 5
0.2
b
0 -
0 0.5 1 1.5 2

V“:

Figure 15. Comparison of Velocity Profiles (Outer Scaling, ske, B)

68

63am 31...... 4.3111. .3... 0 ”
(Outer Sealing)

 

 

 

 

 

Comparison ofVeIoclty Profile at xlw I 1.59 A
(Outer Scaling)

 

 

 

 

7 3an 311921661.) 5.31.3.3...
(mung)

    
  

 

‘_"3</iv‘=3.13.rke i
i D XIW=3.18.LBSTI

 

 

o 0.5 1 1.5 2
via.

#2; w“ -HE... _z_z ___I

Figure 16. CompaTISOn Of VelocityPTOfilesfi(OuterScali-ng,_rke, A)

69

' 66W" Qrvaloclty PrOfilO .. .1. 4 4.45
(Otter Scaling)

    
  

 

 

 

    
  

 

 

 

    
  

 

 

 

1
_* we 145.11... .
0'8 a xlw=4.45.testl
I 0.6
3 0.4
0.2
0 - .
0 0.5 1 1.5 2
visa
W Comparison of Velocity Profile at xlw-7.62
(OW Scaling)
1
,_w=7.62.rke
0.8 u xlw=7.62.tast‘
0.6
3'
0.2
0 v P
0 0.5 1 1.5 2
W5:
Comparison ofVeIocIty Profile at xlw - 10.8
(01hr Scaling)
1
0 8 —— ‘ W168. rite"a
' c) xlw=10.8. test.
} 0.6 y "” ” ’
f '5’ 0.4 .
0.2
o v b
0 0.5 1 1.5 2
via:

Figure 17. comparison art/sass Profiles (Outer scaling, rke, B)

70

Comparison ofVeIocIty Profile at xlw e 0 7
(Outer Scaling)

 

 

 

1.5 2

 

Comparison of Velocity Profile at xlw I 1.59

 

 

 

 

(OuhrScaIIng)
1
0.8 276:1.592‘rsm
g x/w=1.59.test
g 0.6 . —— 4.....2 ,
D 0.4
l 0.2
0 ..
0 0 5 1 1 5 2
Via:

Comparison of Velocity Profile at xlw-3.18
(Outer Scaling)

 

    
  

_ ' "SE/$3.16, Isa?
g D xlw=3.18,test;

 

 

0 0.5 1 1.5 2

Figure 18. Comparison of Velocity Profiles (Outer Scaling, rsm, A)

71

Cornperleon of Velocity Pnollle et xlw - 4.45
(Oubr Scaling)

 

  
    
    
  
  
   
  

>;MA=;.4S.fim-
' u ma.4s.mr

 

 

 

0 0.5 1 1.5 2

Cornperieon of Velocity Prollle et xlw-7.62
(Oubr Scaling)

 

‘-—4x/w;7.82. r‘sm‘.
u xlw=7.62. test

 

 

0 0.5 1 1.5 2
VB:

Cornperleon of Velocity Prollle et xlw - 10.8
(Outer Scaling)

 

xIJ=iOE mi,
0 x/w=10.8. test.

 

 

 

; 0.6
: s'...
0.2

o - I s

o 0.5 1 1.5 2

Y’aa

Figure 19. Comparison of Velocity Profiles (Outer Scaling, rsm, B)

72

ComperleonotTurbuientWecoeltyRetloetx/wao

'—o’” : xlw-=0" Tim?

 

 

 
 

 

 

 

200 , 0'“: '-
,< u xlw-0.!“ i
1505 HAW in W“
Hill-100 4
5° ,
0 """""""
0 0 005 0.01 0.015 0.02
Hm)

Comp-bondTubulentVbcoolty Momma-1.50

 

  
 
  

 

W1 .59. eke 1

.74”: $44.55;.»
i
i _D “"951592'.“

L

 

 

E 100 1
50 J
o v--.- ‘ .
0 0.005 0.01 0.015 0.02
y (In)

Comp-hon of‘l’ubulentVleooelty Ratio et xlw I 3.18
g—o TmeFtiEsm—g
MAB. eke ;
‘ _o Wyn-3.18 j

150

 

    

 

 

 

0 0.005 0.01 0.015 0.02

Figure 20. Comparison of Turbulent Viscosity Ratio (A)

73

Comperleon ot Turbulent Vlecoelty Retlon et xlw I 4.45

 

200
'—°- xlw=4.45, rms

x/w=4.45.ske"
‘ Di xlw:4.45.rke '

 

150

 
  
    

 

 

 

 

W" 100
50
O

O 0.005 0.01 0.015 0.02

Comperleon o! Turbulent Vlecoelty Retlo et xlw I 7.52

 

 

  

 

 

 

250
_"'°' xlw=7.62. rsm
20° ‘ xlw=7.62. ske
WP
100
w .......
0
0 0.005 0.01 0.015 0.02
y (mi
Comperleon of Turbulent Viecoelty Retlo et xlw I10.8
300

 

   
   
 

—°" xlw-10.8.rsm.
25° . xMI10.8,ske1
zoo , Piggyyvarosye i

 

  

 

 

 

up 150
100
50
o
o 0.005 0.01 0.015 0.02
y (In)

Figure 21. Comparison of Turbulent Viscosity Ratio (B)

74

 

 

Comparison of Velocity Profile at xlw =1 .59 (Inner Scaling)

—-——y+=5.5log(u+)+8.83 —x/w=1.59,ske D xlw=1.59,rke
29 y A xlw=1.59,rsm + Exp1.59

 

1 10 100 1000 10000

 

 

 

Comparison of Velocity Profile at xlw =10.8

(Inner Scaling)
--—-u+ = 5.5iog(y+) + 8.83—W108. ske g x/w=10.8. rke
29 j A adw=10.8, rsm 4. Exp 10.8

 

1 10 100 1000 10000

 

Figure 22. Comparison of Velocity Profiles (Inner Scaling)

75

 

 

 

 

Comparison of u’ at xlw I 0

 

__§/w=6, at; ” '. _... .. .
I xlw=0.ErIkeson I
r“ wnm,_u W ,, .

 

 

 

1 10 100

1000

 

Comerleon of u’ at xlw = 5

 

 

 

Companion oi u’ Distribution at xlw I10

 

 

 

i— 11111310316“

: . xlvw=10.Enkseonl

 

 

 

 

Figure 23. Comparison of u“ Distribution

76

10000

 

u'I'

iu'IerScdedWodtyusblbuim
(Wk-eMJm-uyum liTredrmrt)

----- xM=0,ske o xlw=1.59,sle —xlw=8.18,ske
19 -1 a xlw=4.45,sie —)dw=7.62,ske ~~1dw=10.aske
17
15 «j
13
11

dwUINCD

 

Figure 24. Overall Inner-Scaled Predictions of u+ Distribution (ske)

77

inner Scaled Velocity Distribution
(Realizable k-e Model, Two-Layer Zonal Wall Treatment)

    
   

 

 
  
 
  
  

 

19 ~ ,. ,,,,, - k ,
. . 1 a xlw=1.59.rkem
, ‘ ‘ i xlw=3.18,rke* '
14 ‘ - xlw=4.45,rkel
i f ' .—xlw=7.62,rke
'5 9 I’ll. L W199.
4 /I’
I .
-1 '
1 10 100 1000 10000
y-l-
Inner-Scaled Velocity Distribution
(RSM,Two-LayerZonal) 7 7 * We 7_
19—1 - '- - xlw=0.rsm .
' D W1.59, ism J ‘
} xlw=3.18,rsml i
14 1 ‘ - xlwd.45.rsmt ‘
' E—xlw=7.62, rsm '
-. . -~ -x/w=10.8. ' -
g 9 ,1 L -- w-.. .3... 1
, /
4 1.....- «NJ
-1 T ,__-_
1 10 100 1000 10000

y'l-

Figure 25. Overall Inner-Scaled Predictions of 0+ Distribution (rke and RSM)

78

mmmwusmxlw 1.59

 

(OtterSchng)
3.5 -_-e W - u -.
, if _ + Exp1.59 ;
. 3 + ' o xlw=1.59,rke‘
2.51 .. __ ' n )dw=1.59,ske=

 

 

 

 

ShoamuiseNomlStressfiNS)atx/w=10.8

 

 

   

(Outer Scaling)
3'5 ++ - —-)dw=10.8,rsm
+++ ; + Exp10.8
3 4+++ +++++ ; o x/w=10.8, rke
1 g xlw=10. 83113 -0

 

.' ‘\
.' \‘
\

 

 

Figure 26. Comparison of SNS (Outer Scaling, A)

79

CormaflsonofSNSatx/w=1.59(0uter$mling)
4

. f "I €8.50 7 z
3.5 + ‘

 

3 . - ""_"""”=.1-59' new

2.5-
2 -
1.5 -

1

0.5
0 .

   

 

 

 

 

1
Y“:

 

i '4I-OMExp1o.8' j
1 ”F103 '19-'31“?

   
  

 

 

 

 

0 0.5 1 1.5
we

Figure 27. Comparison of SNS (Outer Scaling, B)

80

u'+

u'+

Cornpa'leon of u'+ Detrlbution a w I 1.59 (Inner Scding)

 

 

 

 

 

 

 

 

 

 

5
A xlw=1.59, rle g xlw=1.59, sie
4 + Exp 1.59 .._x/w=1.59, uu_rsm
3 _
2
1
0 I ,
1 10 w 1000 10000 '
Compulsion of u'+ Detrlbutlon a xlw = 10.8 (inner Scdlng)
5
A xlw=10.8, die n xlw=10.8, sle
+ 5010.8 _x/wfi0.8, w_rsm
4
3 .
2 a.
1
0 - .
1 10 (g) 1000 10000

Figure 28. Comparison of u" (Inner Scaling, A)

81

Compulsonofu'+atxlw=0

 

 

 

 

 

 

 

 

 

 

 

 

 

6 —xlw=0.rsm
5 . 9. nyO-Ffiksson
4 1
' 4»
a 3
2
1
0 O .1 ... .00”. f ‘ ,
1 10 100 1000
y+
Comperieonofu'+etxlw=5
6
, w~ . v e
.—-xlw=5.rem
. e xlw=5.Eriksson 9 e
41 T * ” ” “
-; 3
2
0
1 10000
yd-
Comparieon of u'+ at xlw-10
6 . 7
5 + F—xthtonm
é.-. #831“?in
4
4-
.: 3.
2
1
0
1 10000

 

Figure 29. Comparison of u'+ (inner Scaling, B)

82

 

 

Comparison of v'+ at xlw = 5

 

 

 

 

 

 

 

 

 

 

 

 

 

10000
Vt
Comparison of v'+ at xlw = 10
3 . «W,
'_.xlw=10,rsm ‘
2-5 «g . x/w=10,Eriksson; .00
4.
.> 1.5
1
0.5
0 1 .
1 10 100 1000 10000

 

Figure 30. Comparison of 7' (Inner Scaling)

83

 

 

 

 

u'v'+

Comparison of u'v'+

 

_ mama
. xlw=5,Bi<sson e

 

 

1 0000

 

 

 

 

u'v'+

Comparison of u'v'+

 

;__“ 1:710:10, rsrn

A ‘v'..- #3191???" : I'-

 

 

 

 

 

Figure 31. Comparison of E.” (Inner Scaling)

 

 

Comparison of Temperature Profiles at xlw I1.50

 

 

(Outer Scaling)

‘ 1 . t

1 4" 1+ 1+ +1- +1» +1-
0.8 <
0.8 + Exporunental

g. ‘ W159i“?
0.4
0.2
0

 

 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
W8:

Comparison or Temperature Profiles at xlw I4.45

 

 

 

 

 

 

 

(Outer Scaling)
1 + + + + 47
+ 1
0.8
i 0.5 3 + Experimental ‘
‘ {—XIW34.45. SRO
‘- _. __ _. . _
‘ 0.4
0.2
0 -
0.2 0.4 0.6 0.8 1
W52
Wotfwfiofiiedmflu
(Ohm
1 Eu A A ..
“+ ’+ ‘+ '+ H'
0.8
i 0.0 A + W
1 l- + i- - - W108.an
;_ -xfw=10.8.rlo 1
0.4 + ‘—r W- W10.8.d0“
0.2
0

 

 

Figure 32. Comparison of Temperature Profiles in Outer Scaling

85

Comptlson of Temperature Profiles at xlw I10.8 (Inner Scaling)

 

 

 

 

 

 

 

1 1 0 1 00 1000
y-l-

Figure 33. Comparison of Inner-Scaled Temperature Profiles at xlw = 10.8

86

Chapter 4

SIMULATIONS WITH PHASE CHANGE MODELS

The mass transfer of water, especially that which occurs on the interior
surface of the windshield, is the most important factor that affects the visibility of
the driver in an automotive vehicle. Therefore, a successful simulation of the
mass transfer is a crucial step towards computer-aided design for the windshield
defogger optimization.

Analogous to the relationship between temperature distribution and heat
transfer, the concentration distribution of water needs to be predicted in order to
simulate mass transfer. Analogous to the heat transfer, the simulation of the
velocity field forms an important basis of mass transfer predictions as can be
seen in the concentration equation (2.9). On the other hand, the phase change
phenomenon involved in defogging will affect the temperature distribution. With
thermodynamic equilibrium taken into account, the concentration distribution
(mass transfer) is coupled with the temperature distribution (heat transfer) of
water, even if the velocity distribution is assumed to be unaffected by the phase
change.

In this chapter, the impact of phase change (condensation/evaporation) on
heat and mass transfer will be investigated, and the task of modeling the phase
change phenomenon will be considered with respect to the commercial CF 0

code FLUENT 5.0.

87

4.1 Some Basics about Mass Transfer
Before the phase change models are introduced, it beneficial to introduce

a brief review about some basic ideas in mass transfer.

4.1.1 General Mixture
Suppose a single-phased mixture has n species. The mass concentration

(also known as the partial density) of species i', pr, is defined as

,,=_L 4.1
A V ( )

where m.- is the mass of species i' in the mixture sample; and V is the volume of

the mixture sample. The density of the mixture, p, is the sum of the partial

densities:

pip.» (4.2)

("=1

The mass fraction of species i', (Dr, is defined as

<0, = E'- (4.3)
p

which is the main variable involved in the concentration equation

6 a 6
EWr')+'a;(pur¢r ) = ”EJN + R." + Sr" (2'9)

1

and Fick's law is

60>.
J.,.=— D. —'— 2.10
(J p 1.!!! ax ( )

l

88

The following equations hold when the sum of the volumetric rates of

creation for all the individual species is zero:

1 n
u - P H pi'ui'
l I!
= — " " 4.4
v pgnm ()

where u and v are the bulk or mass-averaged velocity components of the
mixture; and Ur and v.- are the x and y velocity components of species i'.
The mole fraction of species i' in a sample mixture containing N moles is

defined as
x. = —" (4.5)
where N.- is the number of moles of i' in that sample. Clearly,
N=:M. Mm
.~=1
and

Za=1 an

m..
N" = ._l 4e8
' Mr ( )
m
N=—- 40

e'_l

.

The relationship between the mass and mole fractions of species i' is

89

xi. =y—d), (4.11)
Mr

It is worth noting that, for an ideal gas mixture,

P.
.. =-—'— 4.12
x, P ( )

where P,- is the partial pressure of species i'; and

n

10:212. (4.13)

r
i'II

4.1.2 Water Vapor and Dry Air
With Fick's law, the velocity component of air in the binary vapor-air

mixture can be written as [10]

ua=u_mpi[_&_), ., =._&m03[_m_] (4.14)
p. 6x p..+p. p.+p.

where u, and Val are the x and y velocity components of air; u and v are the bulk
or mass-averaged velocity components of the mixture.

For a water vapor and dry air mixture, the relative humidity is defined as

 

 

¢= . = . (4.15)

where P. is the partial pressure of vapor in the mixture; and P33, is the saturation
pressure of the vapor. For an ideal gas mixture under a constant total pressure,
P, (say 1 atm),

P P/P x
= V = V = v 4.16
4’ P. P./P 0.)... ‘ ’

 

 

where x" is the mole fraction of vapor; and (xv)sa( is the saturation mole fraction.

90

At a given pressure, there existsa temperature dependent "moisture
capacity" of air, which is the maximum amount of moisture the air can hold. At
this point. the air is saturated, and the relative humidity is 100 percent. Any
further drop in temperature or addition of moisture results in condensation of
water vapor into liquid water (dew) in order to keep the thermodynamic
equilibrium. The dew point temperature, up, is defined as the temperature at
which condensation begins if the air is cooled at constant pressure [19]. Note that
as a mixture cools at constant pressure, the partial pressure of vapor Pv remains
constant until the temperature drops below the dew point.

There exists a definite relation between the saturation pressure, Psat, and
the saturation temperature, Tm, during a phase change process. By treating the
vapor as an ideal gas and the latent heat of condensation/evaporation as a

constant, this relation can be described by the Clapeyron-Clausius equation [19]

5?— ;5 i—J— (4.17)
P, w R T, T2 s...

where L is the assumed constant latent heat at some average value; and R is the
gas constant of water vapor. L, R, T1, and T2 are all expressed in the SI unit

system.

4.2 Surface Phase Change Model

For the time being, the situation of forced convection condensation on a
flat plate in the presence of a non-condensable gas will be taken as an example
to illustrate the ideas of phase change modeling. The analysis of evaporation has

the same nature.

91

4.2.1 Model Analysis

Due to the limitation of the current version of FLUENT, it is difficult to
include the liquid condensate film (see Figure 4) in the simulation. The “thin film"
assumption employed by many researchers mentioned in Chapter 1 is adopted in
the present model.

As in Figure 4 (reproduced below), the liquid film which remains on the

surface is assumed to be infinitesimally thin so that it has no influence on the

 

 

 

 

 

velocity field and heat transfer.
U... T... 1p,
——e-
y
# Water-Mr Mixture
I: t 5; Liqu'd Film
0 '— J

 

T.
The key approximation of the surface phase change model is that the

condensation occurs ONLY on the surface of the plate.

The fluid is assumed to be an ideal gas mixture of water vapor and dry air.
In other words, the flow is in single (gaseous) phase with binary species (vapor
and air) in the computational domain. The boundary of the computational domain
actually starts from the interface between the gaseous phase and liquid film.
Since the phase change (condensation or evaporation) is assumed to occur
ONLY on the surface, there are no droplets of liquid water formed volumetrically
in the computational domain. No thermodynamic equilibrium is required, the

vapor remains in the gaseous phase even when the local temperature is below

92

the dew point. Thus, volumetrically, neither latent heat release nor mass
consumption needs to be accounted for in the gaseous phase.

The latent heat released at the interface does not affect the temperature
distribution because the temperature is always held constant on the surface. The
flow in the computational domain can not "see” the latent heat released when
condensation occurs on the surface, but can see only the constant temperature
condition. By keeping the surface temperature constant, the latent heat released
due to condensation is absorbed into the plate surface instead of the flow.

Again, analogous to convection heat transfer, the driving force of the mass
transfer onto the surface of the plate is the concentration differential of vapor
between the interface and the main stream of the flow.

With the above assumptions, the boundary conditions at the interface can
be defined as:

I The temperature is the same as the one on the wall surface, Ts;

. The mass fraction of water vapor is the one corresponding to the

saturated concentration of vapor at T8;

- The no-slip (u=0), impermeable surface model (v=0) is adopted for the

velocity.

Strictly speaking, the impermeable surface condition imposed at the
interface is another approximation [14] (called impermeable surface model).
Actually, the interface is impermeable to air only, and is ”sucking” (condensation)
or "blowing" (evaporation) for the species of water vapor. The exact velocity

boundary condition at the interface should be

93

 

=0 (4.18)

where V, is the y velocity component of air. With equation (4.14), the equivalent

statement of equation (4.17) is

VI _ =p_a-i-_,0LD_6_ __€2_ (4.19)
'v‘” p. 6y p. +p.

However, the v = 0 approximation is justified only when the concentration of the

species of interest is "low", namely, lower than a critical level [14]:

pw _ poo < SCI—n (4.20)
p

where So is the Schmidt number (0.6 for air-vapor mixture); n is 1/3 when So is
greater than 0.5, and 1/2 when So is less than 0.5. The air-vapor mixture meets
this critical condition in the range of normal room temperature, which makes the

impermeable surface model (v = 0) valid for the current study.

4.2.2 Model Implementation

The ideal gas mixture of non-condensable air and condensable water
vapor was first set up in the ”Define. .. Materials. . ." panel in FLUENT [12], with
air as the main species. The density of the mixture is computed with the ideal gas

law for an incompressible flow [12]

Pop
RT __i'
iM'II

where P0,, is the operating pressure (1 atm). Theoretically, this makes all of the

governing equations mentioned in Chapter 2 coupled together even if the

physical properties other than the density are assumed constants. Therefore, all
the governing equations should be solved together, which makes the simulation
more time-consuming.

As mentioned earlier, there exists a strong analogy between heat transfer
and mass transfer. The temperature distribution can be obtained by solving the
energy equation (2.4), from which FLUENT calculates the wall heat flux output
and make it available to the user. Ironically, the current version of FLUENT does
not provide the species mass flux output on the wall although the species
concentration field has been obtained by solving the species concentration
equation (2.9).

In order to get the vapor mass flux on the wall with equation (2.10), the
derivative of the vapor is needed. A User Defined Function (UDF) was developed
in a C program (attached as Appendix E), compiled, linked and run with the
solver to get the information needed for the mass flux. In this UDF, the mass
fraction of vapor is stored in the first User Defined Scalar, UDS(O) (C_UDSI(c,t,O)
and F_UDSl(f,t,0) in the code), which makes the derivatives of the vapor mass
faction available to the user (C_UDSI_G(c,t,0) in the code). For reference, the
saturated mole fraction corresponding to the local temperature is calculated with
equation (4.23) and stored in the second User Defined Scalar, UDS(1). The
mass flux was calculated with equation (2.10) and stored in the third User

Defined Scalar, UDS(2).

95

All the above operations were realized via the Define_Adjust function (or
Macro [12]) named 'spec_grad" in the code. The solver calls the Define_Adjust
function at the beginning of every iteration [12].

If P1 in equation (4.17) is set to be 1atm (equal to the mixture pressure of
the current study), the saturation pressure of the vapor at temperature T can be

obtained as

P)

(19,)”, = e” T' 7 [atm] (4.22)

where T1 is the temperature at which the saturation pressure of the vapor takes

the value of 1 atm. That is, T1 is the boiling point of water (373.15 K) at 1 atm.
The pressure of the mixture for the current study is about 1 atm. For an

ideal gas, the ratio of vapor partial pressure to the total pressure of the mixture

equals the mole fraction as described in equation (4.12). Therefore,

fibiffl L _'-_l
(3%)”, = (P, )‘m’ = e [atm] .__ efln I]

P l[atm]

 

 

(4.23)

The saturated mass fraction is then obtained via the following equation [11]

_ r(x.)..,
’1-(x,)m,(1—r)

 

(4%)... (4.24)

where r is the ratio of molecular weights of vapor (18 kg/kmol) and the air (29
kg/kmol). This calculation is consistent with equation (4.11), and is valid for non-
saturated mixture as well.

In the UDF of Appendix E, the saturated mass fraction of vapor computed
with equations (4.23) and (4.24) is imposed as the boundary condition for the
plate of interest (via the DEF INE_PROFILE function [12]).

96

According to equations (4.16) and (4.23), the relative humidity, 41, can be

obtained in FLUENT via manipulation in the "Define...Custom Field Functions...”

menu [12]:

 

4.2.3 Results and Comparisons

The Flat Plate Case
To begin with, the 2D geometry in Figure 34 was created for the case of

forced convection condensation on a flat plate in the presence of a non-

 

condensable gas.
Inlet T Outlet
__. .8 _.
°1

 

 

fl 0.8m ’5

Figure 34. 20 Flat Plate Geometry

 

 

 

The distance between the two parallel plates is 0.2 m. The lengths of
plates are both 0.8 m. The vapor-air mixture enters the system from the inlet with
a constant velocity of 3 m/s, and goes out of the system at the outlet where a
constant pressure is assumed. The lower plate is a thermally active one that is of
interest, with a constant temperature and the corresponding saturated vapor

mass fraction. The top plate is adiabatic with a zero gradient vapor mass fraction.

97

As seen in the comparisons in Appendix F , the numerical solutions for
velocity and temperature of laminar flow are in good agreement with the similarity
solutions [17] by Blasius and Pohlhausen without considering mass transfer (In
FLUENT, set the fluid to be air only.). The good agreement suggests the above
20 geometry is acceptable to perform the simulation for a flat plate case
illustrated in Figure 4.

For the mass transfer simulation, an ideal gas mixture of vapor and air is
set up with the ”Define. .. Materials. . ." menu [12]. The mass diffusivity, D, of the
vapor in the mixture is assumed to be a constant. In order to compare the
simulations with analytical solutions, the inlet temperature was first set equal to

the temperature on the lower plate, say 303.15 K (30 °C). And the inlet (or main
stream) mass fraction was set to be corresponding to 50% relative humidity (41).
With equation (4.16), the mole fraction, x... corresponding to the relative humidity
0 is obtained by

x... =¢(x.)... (4.25)
Since FLUENT requires the mass fraction as the input for the concentration
boundary condition, the corresponding mass fraction, (1).... is then obtained with
equation (4.24).

Alternatively, the inlet main stream mass fraction (I).o corresponding to the
relative humidity 0., can be calculated via the so called absolute or specific

humidity, 00, which is defined as

 

a) = m" (kg water vapor/kg dry air) (4.26)
m

98

The relationship between the relative humidity and absolute humidity is

a) = ——°'622¢R°' (4.27)
P - P

sat

where P“, can be computed from equation (4.22); and P is 1 atm.

The inlet mass fraction can then be calculated as

on = m... (4.28)
1+0)“,

 

The local mass transfer coefficient of vapor along the lower plate is of

interest in terms of mass transfer, which is defined as

h J” (4.29)

“In-p...

 

where JW is the mass flux at the wall; pW is the mass concentration at the wall;

and p... is the mass concentration in the main stream. With equation (2.10),

an
J =_ _ 4.30
. po(ay)_0 < 1

where p is the mixture density at the wall; D is the mass transfer coefficient
(mass diffusivity); and (I) is the mass fraction of vapor.

The analytical solution of h,n for the case of laminar flow over a flat plate is

[14]

 

1/3 1/2
hm=o.332[—‘i) [(4‘) 2 (4.31)
D x

V

when v/D 2 0.5 (valid for the air-vapor mixture), where v is the kinetic viscosity of

the fluid (assumed to be the same as dry air, 0.00001568867 m2/s).

99

The numerical solution of hm is obtained by equation (4.29), where JW is
replaced by UDS(2) at the wall. Note that the value of UDS(2) at the wall is the
value stored at the center of the wall-adjacent cell (as seen in Figure 35) since
the face value of the derivative of UDS(O) is not available. This approximation
requires the near wall mesh to be fine enough so that the mass fraction

distribution, UDS(O), can treated as linear.

 

 

Wall-Adjacent Cell _________. Cell Values

/

 

 

 

 

 

e
6/

/
.A’
i

/,. .

Face Values

Figure 35. Illustration of Face and Cell Values
The comparison of the numerical and analytical solutions for hm along the
plate is shown in Figure 36. The ratio of the analytical solution to the numerical
one is given in Figure 37.
The agreement between the analytical and numerical solutions suggests
that F LUENT is an acceptable tool to perform mass transfer problem without

volumetric phase change, at least for the case of laminar flow over a flat plate.

100

Comparison of Predicted Local Mass Transfer Coefficient
Profile with Analytical Solution (Laminar)

 

0.16 I
0.14 .
0.12

f ’ analytical
,_ ._ nurnericai

 

J A

 

 

 

0 0.2 0.4 0.6 0.8

x (m) ‘

Figure 36. Comparison of Mass Transfer Coefficient (Laminar, Flat Plate)

Ratio of Man Transfer Coefficients
(analytical reunite over numerical prediction)

 

 

 

ratio

 

 

 

9 0 0.2 0.4 0.6 0.8
X (m)

A -H -.._._.v.._.- _.__. -._.. h... __.,_..- -ww. .E_._.-__.__.__ .ufi .5

Figure '37. Ratio of Analytical to Numerical -Soiutionfisi'for- hm (Laminar, Flat Plate)

101

The Experimental Facility Case

The same surface phase change model was applied to the experimental
facility.

Similar to the flat plate case, the inlet mass fraction boundary condition is
obtained by equation (4.28). The mass fraction on the thermally active surface is
imposed as saturated according to the surface temperature via the same UDF in
Appendix E.

An important difference from the flat plate case is the treatment of the
turbulence effect on the mass transfer prediction. Instead of using equation
(2.10) for the laminar flow, equation (2.11) is adopted to calculate the mass flux,
Jr], for the turbulent flow. The impact of turbulence on the mass transfer is also
considered in the UOF in Appendix E. (When applied for laminar flows, the
turbulence viscosity pt (C_MU_T(c,t) in the code) is equal to zero.) As seen in
equation (2.11), the prediction of the turbulent viscosity pt is important for the
mass transfer simulation in turbulence flows since it affects the effective mass
diffusivity.

Only three cases of simulation were performed, as it is rather
computationally expensive with the current hardware resources:

- Case 1: The inlet temperature (approximately equal to the jet

temperature) of the air-vapor mixture is 25°C with a relative humidity of
60% (the corresponding mass fraction was calculated with equations
(4.27) and (4.28); the plate temperature is 5°C with saturation vapor

mass fraction.

102

- Case 2: The inlet/jet temperature is 24°C with a relative humidity of
75%; the plate temperature is set to be 4°C with saturation vapor mass
fraction.

- Case 3: The inlet/jet temperature is 25 °C with a relative humidity of

85%; the wall temperature is 5 °C with saturation vapor mass fraction.

The velocity of the jet nozzle is 10 m/s for all three cases.

Figure 38 represents the predictions of local vapor mass flux, Jw, at the
wall. Jw is also the flux of water condensed on the surface of the wall. Note that
the original value of Jw is negative from the CF D code since the direction of it is
out of the computational domain. The absolute value is taken in the figure as it

makes more physical sense in terms of water condensed on the wall.

Predictions of Vapor Mass Flux J. onto the Wall

 

n .+Case1‘

205-03 ' A Case 2 .
Poi???)

 

 

 

 

0 0.05 0.1 0.15 0.2 0.25
x (ml

Figure 38 Predictions of Vapor-Mass Flux on theWalI I
The predictions of mole fraction profiles of water vapor for the three cases

are included as Appendix G.

103

Note that with the current surface phase change model, the mole/mass
fraction is predicted by the concentration equation without considering local
thermodynamic equilibrium. it is possible, especially for higher inlet/jet humidity
(like in Case 3), that the predicted vapor mole/mass fraction is higher than the
local saturation mole/mass fraction. That is, the vapor is super cooled. The super
cooling tends to happen in the near wall region where the temperature is low
(possibly lower than the dew point).

To investigate the super cooling that can happen using predictions of the
current model, the relative humidity profiles at the location x/w =10.8 are obtained
by dividing the predicted mole fraction profile with the saturation mole fraction
profile corresponding to the predicted local temperature. As seen in Figure 39,
the surface phase change model predicts super cooling in the near wail region (0

s y s 3.0e—4 m for Case 3) when the inlet/jet humidity is relatively higher (Case

3).

Comparison of Relative Humidity Profiles at xlw I10.8
i (Surface Phase Change Model)
; 110% .- PP -4 -4 -. -_

l i
l D (:1
i 100% ~ D D --
1 la 0 D u .
a 1 + O D D
1 90% 7 . + O .. O ..
3 . O ,
f 1 + + a
70% 4‘ 7_ , + . 4. 1
l . +Case 1 OCase 2 DCase 3 1
60% ¥ PPP - ~~P - 2...“. Pt»,- xv , . ..
0 0.0002 0.0004 0.0006 0.0008 0.001
11'“)

Figure 39. Comparison of Relative Humidity at xlw =10.8

104

However, whether the super cooling really happens in the real
experiments is questionable. It is generally believed [9] [10] [11] that as the local
temperature drops below the dew point, there will be water droplets formed to
keep the local thermodynamic equilibrium, as will be considered later in the
volumetric phase change model. The current model assumes the vapor stays in
the gaseous phase even when it is super cooled.

Multiplying the averaged values of the mass flux predictions in Figure 38

(K) by the plate area (A), the global condensation rate, mp, for the plate can be

obtained: "'1’, = 7:.4 . To compare with experimental drip-off tests, the unit of mp

is defined in terms of grams per minute (g/min).

 

 

 

 

Corrpar‘lsonofCondensatlonMes
1 2
: mm 1) = 112,“ - 0.0815 ‘ Make) = 141 .31x + 0.013 X
i 1.75 R230.“ _ +
. 15 “test 2) . 130.58: - 0.4315 M“) ‘ 131'“ + ”209
I g 15* 122209735
i v managers-02281 . x
. 5 125“ R'=0.9895 mm
1 1‘ ¥ r-W
3 E 0.75. MIME”
5 0.5 Q I
i e w 'I'
‘ 0.25 let
0 -
0.000 0.002 0.004 0.008 0.008 0.010 0.012 0.014
e -A e 912,886,- _ g e _. _w
e test1 A test2 o test3 x ske
+ rka ~ Linearitest1) Linear(test2) Linear(test3) .,
~ Hmlskel ‘1,-U"°ar-rl<.°)--, ‘

Figure 40. Comparison of CondenHSation Rate of the Plate

105

Figure 40 is the comparison of the predicted condensation rate on the
plate with experimental drip-off tests. The abscissa (Apv, kg/m3) in the plot is the
difference of vapor mass concentration between the inlet/jet and the surface of
the wall. As seen in the figure, the predicted condensation rate is very linear. The
intercept of the curve is very close to_ the origin, which suggests that as the
concentration differential is zero there is no drip-off of water.

The slope of the curves have the meaning of the averaged mass transfer

coefficient, Hm, if divided by the area of the plate. The linearity of the predictions

suggests that Hm does not vary with the concentration differential so long as the
jet velocity stays the same (10 m/s). As seen in the figure, the predicted slope is
close to the experimental data, while the absolute values of the condensation
rate are less than 20% off from the tested data. The realizable k-e model gives
lower condensation rates because it predicts lower turbulent viscosity. As seen in
equation (2.11), a lower 0. will result in lower predictions for the mass flux. The
author offers the following possible reasons for the discrepancies in the above
comparison:

- The mass flux of the numerical calculation accounts for the water
VAPOR that goes onto the interface between the gaseous mixture and
the liquid film. It is not clear whether the total amount of vapor going
through the interface will COMPLETELY condense into water liquid. If
not, it may be a reason why the experimental data have lower values

than the predicted ones.

106

I it is suspected that, experimentally, there may exist a flow transition
from laminar to turbulence on the plate surface. if the flow is laminar-
like at the developing region close to the jet nozzle, the mass transfer
of vapor onto the plate will be reduced. However, this numerical
simulation does not include the flow transition.

I Experimentally, the temperature in the area close to the edges of the
plate may be higher than it is supposed to, under the influence of the
ambient flow with higher temperature. This will result in higher
saturation mass fraction in the corresponding area, and lower
concentration difference between the main stream and the plate, which
will reduce the mass flux of vapor onto the plate.

I Another possible reason for the disagreement may be the three-
dimensional effects at the end of the thermally active plate. The
numerical simulation is performed with a 20 model, which is an
appropriate approximation close to the jet nozzle according to the
experimental measurements in the z direction. But at the end of the
plate, the 2D approximation was not experimentally validated. If the
three-dimensional effect is significant in this region, it may be a factor
that causes the disagreement.

Sh = hmw

The Sherwood number based on the slot width, D , in a study

 

similar to the current one was measured by Mabuchi et al [20], where hm is the
local mass transfer coefficient, w is the width of the jet nozzle, and D is the mass

diffusivity. The comparison of the numerical calculations and the published data

107

for the turbulent case (tripped data in [20]) is shown in Figure 41. In this case, the
flow is turbulent from the beginning of the wall jet, and there is no flow transition

from laminar to turbulence.

Comparison of Sherwood Number

 

 

140 . 77 77 7 7 7 7

7 I Re=1.78E+04, Mabuchietal .
120 A. Re=1.30E+04,ske
100 . 7 A Re=9.14E+03.Mabuchietal

 
 
 
 

A I
“A“ .I'IIII....-

 

 

 

 

A
40 AAAAAAA“ ‘d‘q
201
o ,
0 2 4 6 8 10 12 14
xlw

’ "Page '41. 60.1.3430; ofSiierwood number-) 1‘

The possible explanation for the above comparisons is that, in the
developing region of the wall jet, the Sherwood number may depend more on the
velocity of the jet, UM, than on the width of the jet nozzle, w. The higher the jet
velocity is, the higher the Sherwood numbers are. Although the slot width based
Reynolds number for the current study (Re 5 1.3e+4) lies between the ones in

the published data, the jet velocity for all the published data ( 2 14.4 m/s) is

greater than that for the current study (10 m/s).

4.3 Volumetric Phase Change Model

4.3.1 Model Analysis

108

The simple case of the forced convection condensation on a flat plate in
the presence of a non-condensable gas situation is again employed as the
example for the model illustration.

in the volumetric phase change model, the same "thin film” assumption is
adopted as in the surface phase change model, which means that the liquid film
on the wall has no influence on the velocity field and the heat transfer. Moreover,
the impermeable surface model is still assumed valid at the interface between
the liquid film and the gaseous vapor-air mixture. The computation starts from the
interface where the mass fraction of vapor is saturated.

The key point of the volumetric phase change model is that the possible
formation of liquid droplets (fog) in the volume of the flow is considered, which is
the major difference compared to the surface phase change model. In the
previous surface phase change model, the volumetric condensation is not
considered, which means that the thermodynamic equilibrium is only maintained
on the wall surface. Thus, there will be no phase change occurring volumetrically
even if the vapor is super-cooled. However, as mentioned earlier, when the local
temperature is below the dew point, it is possible that the water vapor will
condense into liquid water droplets (fog) to keep local thermodynamic equilibrium
in the volume of the flow. When the condensation occurs, there will be latent heat
released, which will affect the temperature distribution. Meanwhile, there will be a
loss (or sink) of mass in the gaseous phase as the water vapor condenses out
into the liquid phase. Analytically, the condensation phenomenon will cause an

additional source term (sink term for evaporation) in the energy equation, and a

109

sink term (source term for evaporation) in the vapor concentration equation and
the overall mass continuity equation, assuming that the phase change does not
affect the form of the momentum (Navier—Stokes) equations.

Note that in FLUENT, the user can not manipulate the basic forms of the
governing equations. The only possibility for to the user is to use the user defined
source terms in these equations.

In order to account for the impact of the volumetric phase change, the
source terms in the energy, concentration, and continuity equations need to be

figured out, and expressed in terms of the flow variables being solved in the

 

governing equations.

Besides the same assumptions as in the surface phase change model, the
additional assumptions in the volumetric phase change model include:

I The volume of the condensed droplets in the gaseous phase is

negligible so that classical Navier-Stokes equations are still valid.

I The velocity components of gas and dispersed phases are equal.

I The temperature of the gas and the dispersed liquid phase are equal.

I There is no interaction between the droplets.

With the above assumptions, the only effects of the volumetric phase
change are the additional source terms in the energy, the vapor concentration
and the overall continuity equations. The fluid is still modeled as a single
(gaseous) phase, dual species mixture (air and vapor) in the computation. The

phase change effects are accounted for with the source terms.

110

To develop the source terms needed, the control volume shown in Figure

42 is considered in the gaseous phase.

H a
[p.v. +5 (m. >444

a _
pvuv +_— pvuv Ax]Ay
p.u.Ay [ 6y( )

——>l

a
33
B
3

,.

n'rfoAy

P
J:
I?

 

 

L—Ax—4 ’*

Figure 42. Control Volume for Phase Change Analysis
First, the mass conservation of water vapor is considered in the AxAy
control volume. The conservation principle is the statement that the net flow of
the water vapor into the control volume must equal the rate of vapor

accumulation inside the control volume [14]:

%MAy = p..u.Ay -[p.u. + Bax-(p.u.)Ax +~~]Ay
(4.32)

+ pvvax — [pvvv + gy—(pvvv )Ay + - - :le + mfoAy

111

 

The velocity components uV and v, refer strictly to the motion of water vapor

relative to the control volume. They should be distinguished from the mass-

averaged velocity components (u and v) of the bulk flow in equation (4.4).
Worth noting is that the product pvuv represents the mass flux of water

vapor in the x direction, per unit area normal to the x direction. The term 01': on

the right hand side of equation (4.32) is the volumetric rate of the creation of the

. ‘

water vapor, which is also the source term, 8,... needed in the vapor

concentration and overall continuity equations to account for the phase change.

The vapor generation rate 111’: is negative (sink) when the water vapor is

 

consumed during the condensation process.

Also note that since the fluid in the computation accounts for the gaseous
phase only, the relationship of velocity components of the individual species and
that of the bulk flow represented in equation (4.4) does not hold any more. The

reason is that the sum of the volumetric rates of creation for the two species (air
and vapor) is not zero, but tint.

The next step is to write the source term m: (or S...) in terms of the flow
variables solved in the governing equations mentioned in Chapter 2.

Dividing by AxAy and invoking the limit (Ax, Ay) —> 0, the vapor

conservation statement (4.32) can be reduced to

mr=%+%(p.u.)+%(p.v.) (4.33)

Since the present study is for the steady state, the source term Sm can be written

as

112

 

S... = mr= §@,u,)+%(p,v,) (4.34)

Similar to the analysis for equation (4.14), the velocity components of vapor can

be written as

,7 =,_&:&D_‘Z[_£v_], =,_mDi[_£.__] (4735,
I). 6x p.+p. p. 6y p.+p.

According to the statement of the volumetric phase change model, the air
and vapor are the only species in the single gaseous phase. Therefore, p = (3V

+pa. By the definition of mass fraction in equation (4.3), equation (4.35) can be

 

 

further written as

u, =u— {Di—((1),), v, = v- (I: 075(5) (4.36)
and

P. = pq)‘, (4.37)

Substituting equations (4.36) and (4.37) into equation (4.34), the source term Sm

can be written as

s..=—a—[p¢.(u- 1 Biting]+—a—[p¢.[v——‘—D—a—(<I>.)]] (4.38)

ax 0,81: 8y (Dvfiy

 

Now, all the variables in the expression of Sm are the ones being solved by the
governing equations in FLUENT. This source should go into both the vapor
concentration equation (2.9) and the overall continuity equation (2.1). Term R.-- in

equation (2.9) drops out since there is no chemical reaction.

113

 

Accompanying the mass source term of vapor due to phase change, there
is an energy source term 8., for the energy equation (2.4) to account for the latent
heat released during the phase change. Apparently,

S, = —LS,,, (4.39)

where L is the latent heat of water vapor, which is assumed to be a constant in

the normal temperature range inside a vehicle.

4.3.2 Model Implementation

Another user defined function (UDF) was developed to implement the
volumetric phase change model, which is included as Appendix H.

introducing the condensation into the computation is accomplished by
calculating the saturation mass fraction of water vapor according to the local
temperature (0.)”, via equation (4.24), and then comparing it to (by, the vapor
mass fraction at the same location returned by the solver via the concentration
equation (2.9) with source term Sm included. If the (D, is greater than (chasm, it
means that the vapor is super-cooled at this location. In order to keep the local
thermodynamic equilibrium, (D, is adjusted down to the local (0),)...» Physically,
the adjustment introduces a vapor mass sink and the release of latent heat in the
volume of the flow where condensation occurs. Mathematically, the adjustment
produces finite source terms Sm and 8., in the governing equations. On the other
hand, if (by is smaller than ((003,... no adjustment is made for (1)., and the source

terms Sm and 8.. are both zeros. The solver keeps iterating to solve the governing

114

equations with the above-mentioned adjustment until a converged solution is
reached.

in the UDF included as Appendix H, the adjustment of vapor mass fraction
is made in the DEFINE_ADJUST function (macro) named as spec _grad. For
comparison, the local saturation mass fraction of water vapor (0).)“. is stored in
the user defined scalar UDS(3). The adjusted vapor mass fraction is stored in
UDS(O). The vapor mass flux components p,,uV and pvvv are stored in UDS(1)
and UDS(2), respectively. Again, the purpose of storing the flow variables into a
user defined scalars is to get the relative derivatives of the variables that can not
be returned directly by the solver to the UDF but are necessary to specify the
source terms needed in the governing equations. The cell values of the
derivatives of the UDS are available to the user when developing the UDF code.
None of the transport equations of the user defined scalars (UDS's) is solved.

One of the difficulties about using the UDS is that assigning the cell values
for a UDS (as seen in Figure 35) does not automatically assign the face values to
it. if no special assignment is made and the transport equation of the UDS is not
solved, the face values will be defaulted as zeros. However, the face values of
the above mentioned UDSs are very important because they actually serve as
the boundary conditions for the flow variables the UDSs stand for. Therefore, the
face values of the UDSs should be assigned to be consistent with the
corresponding cell values.

In the current UDF, UDS(1) and UDS(2) are functions of the derivatives of

UDF(0) as shown in equation (4.36). However, unlike the cell values, the face

115

 

 

values of the derivative of a UDS are not available to the user, which makes it

difficult to specify the face values of UDS(1) and UDS(2). Although currently

there is no good solution to this problem according to the technical support

personnel of FLUENT, lNC., the author tried two methods to work around this

problem.

(1) One is to specify the face values of UDS( 1) and UDS(2) to be the same as
the cell values at the wall-adjacent cells (as in the UDF of Appendix H)

(2) the other is to extrapolate face values UDS(1) and UDS(2) with their cell
value and derivatives at wall-adjacent cell (as in the UDF included as

Appendix I).

4.3.3 Results and Comparisons

Unfortunately, no valuable results were obtained by the above mentioned
implementation, even for the flat plate case.

Through observation of the behaviors of the residuals, the iteration is
highly unstable. The residuals jump back and forth as the iteration goes on. The
computational results vary from one iteration to another, and never come to a
stable and converged solution.

The unstableness in method (2) of the implementation is much worse than
that in method (1 ), which suggests the iteration is sensitive to the face values
(boundary conditions) of UDS(1) and UDS(2) (p.,uv and pvvv). Although it gives
better stability compared to method (2), method (1) is theoretically an incorrect

way to get the face values of UDS(1) and USD(2) since the vapor mass fraction,

116

 

(1),, on the wall surface is definitely different from that at the center of the wall-

adjacent cell. Besides, as mentioned earlier, none of the ways gives stable and

converged results.

4.4 Summary and Discussion

Despite the questionable assumption that no water droplets (fog) form
volumetrically to keep the thermodynamic equilibrium, the practical value of the
surface phase change model for defogger designs is yet to be further
investigated by comparison with experimental measurements.

It is believed that the volumetric phase change model is in the concept-
ready stage with the theoretical analysis of the model if source terms
manipulation is the way to handle the phase change problem. The author offers
two possible reasons for the incompleteness of the model implementation.

I From the standpoint of model implementation, apparently imposing the
boundary values of pvuv and p,,vv by method (2) is still not appropriate
or needs to be refined.

I Form the standpoint of the the formulation source terms, Sm and Sh are
highly non-linear with respect to the flow variables solved in the
governing equations, which makes the equations more difficult to
solve. Since the algorithm to solve the equations is not open to the
users, it is uncertain whether the algorithm is compatible with the non-

linearity of the source terms.

117

 

Chapter 5

COCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK

5.1 Simulations without Phase Change

A 20 steady state simulation of the wall jet experimental facility was
performed with three different turbulence models, namely, the standard k-e
model, the realizable k-e model, and the Reynolds stress model (RSM). The two-
layer zonal model was employed as the near wail treatment for all three
turbulence models.

in outer scaling, the numerical results of velocity field from all three
turbulence models are in good agreement with the experimental data. The
standard k-e model gives even better mean velocity profile predictions, especially
in the relatively well-developed region of the jet. in addition, it gives better match
of the ym/Sz prediction to the experimental measurements compared to the other
two models. The Reynolds stress model (RSM) provides better predictions for
the turbulence fluctuations. it is not generally expected that the standard k-e
model (or the any version of k-s model) will out perform the RSM in predicting the
mean velocity profiles. The reason for the unexpected model performance in the
current study is not clear.

In inner scaling, the numerical results of the velocity field by the different

models are virtually identical in the near wail region (y“ < 100), but all the

118

 

 

predictions for the mean velocity are quantitatively smaller than the experimental
data.

The temperature simulation results of different turbulence models are
virtually identical, and in good agreement with the experimental test in outer
scaling. in inner scaling, the predictions for the mean temperature profiles are
lower than the experimental data.

Compared to the measurements, the velocity and temperature predictions
in the developing region (xlw s 7.62) are not as good as that in the well
developed region (xlw 2 10) since the turbulence models were originally
developed for fully developed turbulent flows. For a simulation with a commercial
CFD code, there is little room for improvement of predictions unless more

suitable models for developing turbulent flows appear.

5.2 Simulations with Phase Change

Two phase change models (surface and volumetric) were developed to
simulate the mass transfer of water vapor to the cooled wail of a wall jet.

The surface phase change model is realizable from the viewpoint of
implementation, and is computationally stable as no user defined source terms
were introduced in the governing equations. The shortcoming of this model is the
questionable assumption that the phase change ONLY occurs on the wall
surface, which may be one of the reasons for the discrepancy between the
numerical prediction and experimental drip-off test data. Numerically, the

validation of this model needs to be further investigated in the following range:

119

 

I The validation of the impermeable surface model mentioned in Chapter
4, which makes the mass transfer completely analogous to the heat
transfer.

I The limits on velocity, temperature, and vapor concentration conditions

in order for the key assumption of the model to hold.

The volumetric phase change model is believed to be in a concept-ready
state. Several possible reasons may account for the incompleteness of the model
implementation with the commercial code. Further work is needed:

I It is difficult to implement the volumetric model even with the aid of the user
defined function. if feasible at all, the user defined function for the model
implementation needs to be further refined and improved in a much deeper
level with the technical support from FLUENT, INC..

I The validation of the impermeable surface model may play an even more
important role in the volumetric phase change model compared to the surface
phase change model, noting it is not exactly the correct boundary condition
for velocity on the interface.

I The possible investigation of the compatibility between the algorithm of the
code and the highly nonlinear source terms.

For both the surface and volumetric models, the impact of the turbulence
viscosity prediction on the vapor mass transfer, and the possible flow transition
from laminar to turbulent on the wall may be interesting for future studies.

Although there may exist plenty of room for development from the users'

side, as a commercially available software, the current version of FLUENT is not

120

 

developed enough to handle the mass transfer simulation with phase change. if
possible for the algorithm of the code, it would be beneficial to the users if the

code

has the velocity components of each individual species available so
that the user can both specify them as boundary conditions (vv at the
wall equals zero for the current case) and use them in user defined
source terms;
I has the gradients of the all the flow variables solved in the governing

equations available in the user defined functions;

 

I allows the user to specify the movement of the wall with user defined
functions;

I automatically assigns the face value of a user defined scalar
consistent with its cell value when the transport equation of the scalar

is not being solved.

121

 

APPENDICES

122

 

APPENDIX A

SIMULATIONS OF OUTER-SCALED SNS

SNS Simulation (ske, Outer Scaled)

 

2.5

 

 
  
   
 
 
     

 

 

 

 

 

 

 

T xlw-0. we
2 D xlw-1.59. site;
. f A inn-4’ .45, sire.
1.5 . . .\ X xlw-7.82
1 ‘ I W‘tla’ m
0.5 .
0
0 0 5 1 1.5 2 2 5
v15;
SNS Simulation (rite, Outer Scaled)
2.5 . . .
i xlw-0, rite
‘ D W159. I‘kO
2 0 xlw-3.18, rite ‘
~ .. xlw-4.45, rue;
1-5 ' x W762. rite?
W108.“
1 1
0.5
0
0 0 5 1 Via: 1 5 2 2 5
SNS Simulation (RSM, Outer Scaled)
2.5 . , -

 

xlw-0, ram

0 xlw=1.59. ram

0 xlw-3.18, ram
' xlw-4.45, rsm .

+ xlw=7.82, rsrn 7‘

xlw-107.8, ram ‘

  
 
 
 
   

 

 

 

123

SIMULATIONS OF INNER-SCALED FLUCTUATIONS WITH RSM

APPENDIX B

RSM Simulation of u'+ (inner Scaled)

 

 

 

  

 

 

 

 

 

 

   
  

 

 

3.5
. I“! xflv=0TrsmA Cl 3t/w=1.59,rsm‘
3 ‘ 1 o xlw=3.18,rsm xlw=4.45. rsm ,
‘ + xlw=7.62,rsm -xlw=10.8,rsm
i + 2 4 n
' I
a 1.5
1 -
i
0.5
0
1 10 100 1000 10000
y+
1‘ RSM Simulation of v'+ (Inner Scaled)
; 2.5 _,- , 7 ,, ..
' . xlw=0,rsm D xlw=1.59,rsm*
I 2 I 0 xlw=3.18.rsm xlw=4.45,rsm.
L + Nix/£762,7rsm xm=10.8,rsm,
1.5
.t
; >
1 ‘
0.5
i . "
I 0
1 10 100 1000 10000
17"

124

 

w'+

 

RSM Simulation of u'v'+ (Inner Scaled)

 

 

  

 

 

 

 

 

 

 

 

4 __ ,5. 2m , .52
I xlw=0, rsm E] xlw=1.59. rsm
3 o xlw=3.18, rsm xlw=4.45. rsm
.1 +1 xlw=7.62,rsm xlw=10.8,rsm:
:1 1
3
0 P 1.— . ~
‘I 1 10900
-1
-2
y+
RSM Simulation of w'+
3 -- , -, _
xlw=0, rsm u x/w=1.59, rsm
2.5 . 0 XM=3.18. rsm xlw=4.45, rsm
.1. xlw=7.62.rsm xlw=10.8,rsm
L P , 7‘: "1»
2 ‘ In“; 3::
2' :2”
1-5 -._ 5P
'D”: +83 I123
U . 3 z,
1 0°82: I P 5):}:
“it” + “‘b :25:
0.5 - 006538 .«m 2128.
o J. . ,,
1 10 100 1000 10000
y-l-

125

 

44 T*

APPENDIX C

SIMULATION OF TEMPERATURE PROFILS (OUTER SCALING)

 

 

 

1 . 9H"
09 I
. W
. «a?
0.8 3:-
0.7 ’ if
x/w=0,ske
0'6 1' u xlw=1.59,ske‘
0.5? I o )dw=3.18,ske
0.4 g xlw=4.45, skefi
0.3 + x/w=7.62. eke»
0.2 . xlw=10.8,sie'
0.1 '
0 i
0 0.2 0.4 3,152 0.6 0.8 1

126

 

APPENDIX D

SIMULATION OF THE lNNER—SCALED TEMPERATURE PROFILE AND THE
HEAT TRANSFER COEFFICIENT

Simulation of Temperature Profiles (Inner Scaling)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

20 we am a L
9—log-law form 1
| xlw=0. aka 1 q
o xlw=1.59.ska , r
15 .. o xlw=3.18, ske ; 5' -_
A xlw=4.45, ske ‘
+ xlw=7.62,ske 1 1
' 11 " :"._’£’_W=19:.8v €595
,1 10 _,..
,... 0
5 + ’5‘.
w'w,"
0 9’”
1 10 y+ 100 1000
Prediction of Surface Heat Transfer Coefficient for
the Wall Surface
200 ,
_wall surface, ske
150 " f
.c 100 4
50 a J
O
0 4 8 12
xlw

 

 

 

127

APPENDIX E

UDF TO IMPOSE SATURATION MASS FRACTION AS BOUNDARY
CONDITION AND CALCULATE THE VAPOR MASS FLUX

(SURFACE PHASE CHANGE MODEL)

#include "udf.h"

#define L 2400.0e3 /* latent heat of water [J/kg] */
#define b L*18 / 8314.4

#define a b / (273.15 + 100.)

DEFINE_ADJUST(spec_grad, domain)

/* This function assigns the saturated mole fraction of
vapor to UDSI(1), and calculated mass fraction of
vapor to UDSI(O) */

Thread *t;
cell_t c;
face_t f;

thread_loop_c (t,domain)
{
begin_c_loop_all (c,t)
{
float MLFS = exp(a - b / C_T(c,t));

if (NULL != THREAD_STORAGE(t,SV_UDS_I(1)) &&
NULL != THREAD-STORAGE(t,SV_UDS_I(O)))
{
C_UDSI(c,t,l) = MLFS;
C_UDSI(c,t,O) = C_YI(c,t,O);
I
I
end_c_loop_all (c,t)

}

thread_loop_f (t,domain)
{ if (NULL != THREAD_STORAGE(t,SV_UDS_I(1)) &&
NULL != THREAD_STORAGE(t,SV_UDS_I(O)))
{ begin_f_loop (f,t)
{ float FMLFS = exp(a - b / F_T(f,t));

F_UDSI(f,t,l)
F_UDSI(f,t,O)

FMLFS;
F_YI(f,t,O);

}
end_f_loop (f,t)

128

 

}

thread_loop_c (t,domain)
/* Stores the y direction mass flux
of vapor in UDS(2) */
{
begin_c_loop_all (c,t)
{
if (NULL != THREAD_STORAGE(t,SV_UDS_I(O)) &&
( ))

NULL != T_STORAGE_R_NV(t,SV_UDSI_G o )

{
float diff_eff = 2.88e-5 +
C_MU_T(c,t)/O.7);
C_UDSI(c,t,2) = - C_R(c,t)* diff_eff *
C_UDSI_G(c,t,O)[l];
}

}
end_c_loop_all (c,t)

}

 

DEFINE_PROFILE(plate_mf, t, position)

/* This function specifies the concentration profile of
water vapor as saturated on boundary wall surface
according to the local surface temperature. */

face_t f;

begin_f_loop (f,t)

{
float FMLFS = exp(a — b / F_T(f,t));
F_PROFILE(f,t,pOSition) = (FMLFS * 18./29.) /
(l. ‘ FMLFS * (l. - 18./29.));

}
end_f_loop (f,t)

129

 

APPENDIX F

COMPARISON OF NUMERICAL PREDITION WITH THE SIMILARITY

SOLUTION (LAMINAR FLOW, FLAT PLATE)

Comparleon wlth Bladua Solution

 

 

 

 

 

 

 

 

 

 

1.2
1 I
0.8
g 0.6 +
0.4 I ~ . .
__Blasuus Solution
0.2 a . MmericaISolution;
0
0 2 4 6 8 1O 12 14
1]
Compardon with Pohlhausen Solution
1.2
1 4 .
0.8 1
¢ 0.6
0.4 ,1 . _ ‘
l—thlhausen SOIUIIOD
0.2 . NJmarical Solution 1
0
0 5 10 15

130

 

APPENDIX G

PREDICTION OF MOLE FRACTION PROFILES OF WATER VAPOR BY THE
SURFACE PHASE CHANGE MODEL

Case 1:

MoleFracflon ProfilesowaorforCase1

 

IT'S—"m“,‘éié "W" ‘2 raw-3.59: ske * Axlw=3.18, éké

0.024 «
Xx/w=4.45, ske +xlw=7.62, ske Ox/w=10.8, ske

 

gggggmm 0 r em

0.016 +0

Mole Fraction

0.012

 

 

0.008 . -
0 0.002 0.004 x(m) 0.006 0.008 0.01 ;

131

 

WQJJQREISE e I
. .0 p .

p
8
on

  

Mole Fraction of Vapor for Case 2

 

; * game, Ska * D xlw=1.59, ske A xlw=3.18, Ska”

' x W445, ska + xlw=7.62, ske o xlw=10.8, ske ‘
ma; gamma. 2. a a ama. WW
0 [[2ng

 

r

 

 

x (m)

132

0 0.002 0.004 0.006 0.008 0.01 ‘

 

Case 3:

 

 

 

Mole Fraction ofVaporfor Case 3
0.024 ".an'.;3 a :55 , .(.R.V/;x‘.71..7
° onafiiibbbbb
. ° :1 5
1 5 0.02 1 C9 6 oX/w=0,sle .
i=3 3% ’gx/w=1.59,sle
i 8 1
ill- 0.016 . f Axlw=3.18,sle
i 3 xlw=4.45, skef
I o 1
3 E +X/W=7.62, ske:
i 0.012 No xlw=10.8, ske;
0.008 , T
0 0.002 0.004 0.006 0.008 0.01

x(m) V 5

133

 

#include "udf.h"
#include "sg.h"

APPENDIX H

USER DEFINED FUNCTION FOR THE VOLUMETRIC PHASE CHANGE

MODEL (METHOD (1 ))

#define L 2400.0e3 /* latent heat of water [J/kg] */
#define b L*18 / 8314.4
#define a b / (273.15 + 100.)

DEFINE_ADJUST(spec_grad, domain)

{

Thread *t;

cell_t c;
face_t f;

float MLFS; /* saturation mole fraction */
float MSFS; /* saturation mass fraction */
float UVPR; /* x velocity of vapor */
float VVPR; /* y velocity of vapor */

thread_loop_c (t,domain)

{

begin_c_loop_all (c,t)

{

}

MLFS = exp(a - b / C_T(c,t));
/* equation (4.23) */

MSFS = (MLFS * 18./29.) / (1. - MLFS *
(1. - 18./29.));
/* equation (4.24) */

if (C_YI(c,t,0) > MSFS) C_YI(c,t,0) = MSFS;
if (C_YI(C,t,O) <= 0.0)

C_YI(c,t,O) = 1.0e-12;
/* adjust the mass fraction of vapor if it is
higher than the saturation value. */

if (NULL 1= THREAD_STORAGE(t,SV;UDS_I(3)) &&
NULL 1= THREAD_STORAGE(t,SV_UDS_I(0)))
{

C_UDSI(c,t,3)
C_UDSI(c,t,O)

MSFS;
C_YI(c,t,0);

end_c_loop_all (c,t)

134

 

thread_loop_f (t,domain)
{
if (NULL l- THREAD_STORAGE(t,SV_UDS_I(3)) &&
NULL 1: THREAD_STORAGE(t,SV;UDS_I(O)))
{

begin_f_loop (f,t)

{
float FMLFS = exp(a - b / F_T(f,t));
float FMSFS (FMLFS * 18./29.) /
(1. - FMLFS * (1. - 18./29.));

F_UDSI(f,t,3) = FMSFS;
F_UDSI(f,t,O) = F_YI(f,t,O);

}
end_f_loop (f,t)

}

thread_loop_c (t,domain)
{
if (NULL la THREAD_STORAGE(t,SV;UDS_I(O)) &&
NULL 1: T_STORAGE_R_NV(t,SV;UDSI_G(O)))
I

begin_c_loop_all (c,t)
{
float diff_eff = 2.88e-5 +
(C_MU_T(c,t)/0.7);
UVPR = C_U(c,t) ~ diff_eff *
C_UDSI_G(c,t,0)[OI/C_UDSI(c,t,O);
VVPR = C_V(c,t) - diff_eff *
C_UDSI_G(c,t,0)[1]/C_UDSI(c,t,O):

C_UDSI(c}t,1)

C_R(c,t) * C_UDSI(c,t,0)
* UVPR;

C_R(c,t) * C_UDSI(c,t,O)
* vva;

C_UDSI(c,t,2)

}

end_c_loop_all (c,t)

}

thread_loop_f (t,domain)
/* assign face value with way (1) */
{
if (NULL 1: THREAD_STORAGE(t,SV_UDS_I(0)) &&
NULL 1: THREAD_STORAGE(t,SV;UDS_I(l)) &&
NULL s= THREAD_STORAGE(t,SV;UDS_I(2))I

begin_f_loop (f,t)

{

cell_t cell = F_C0(f,t);
Thread *c_thread = THREAD_TO(t);

135

 

if (NULL I:
T_STORAGE_R_NV(c_thread,SV_UDSI_G(0)I)
{
F_UDSI(f,t,1) =
C_UDSI(cell,c_thread,1);
F_UDSI(f,t,2) =
C_UDSI(cell,c_thread,2);
}
}

end_f_loop (f,t)

}

DEFINE_PROFILE(p1ate_mf, t, position)
/* specify saturation mass fraction at the boundary */

face_t f;
begin_f_loop (f,t)
{
float FMLFS = exp(a - b / F_T(f,t));
F_PROFILE(f,t,position) = (FMLFS * 18./29.) /
(1. - FMLFS * (1. - 18./29.));

}

end_f_loop (f,t)

}

DEFINE_SOURCE(mass_src, c, t, dS, eqn)

/* source term of continuity and concentration equation */

{

float source;
if (NULL 1= T_STORAGE_R_NV(t,SV_UDSI_G(1)) &&

NULL 1: T_STORAGE_R_NV(t,SV;UDSI_G(2)) &&
NULL l= THREAD_STORAGE(t,SV;UDS_I(0)I &&
NULL 1= THREAD_STORAGE(t,SV;UDS_I(3)I)

{
if (C_UDSI(c,t,0) < C_UDSI(c,t,3))
source = 0.;
else
source = C_UDSI_G(c,t,1)[0] +
C_UDSI_G(c,t,2)[1];
/* equation (4.34) */
}
dS[eqn]=0;

return source;

}

DEFINE_SOURCE(energy_src, c, t, dS, eqn)
/* source term for energy equation */

{

float source;
if (NULL I= T_STORAGE_R;NV(t,SV_UDSI_G(1)) &&

136

 

NULL 1: T_STORAGE_R_NV(t,SV;UDSI_G(2)) &&
NULL 1= THREAD_STORAGE(t,SV_UDS_I(O)) &&
NULL (a THREAD_STORAGE(t,SV;UDS_I(3)II

if (C_UDSI(c,t,0) < C_UDSI(c,t,3))
source = 0.;
else
. source = -L * (C_UDSI_G(c,t,1)[O] +
C_UDSI_G(c,t,2)[1]);
}
dS[eqn]=0;
return source;

}

 

137

 

APPENDIX l

USER DEFINED FUNCTION FOR THE VOLUMETRIC PHASE CHANGE
MODEL (METHOD (2))

#include "udf.h"
#include 'sg.h"
#define L 2500.0e3 /* latent heat of water [J/kg] */

#define b L*18 / 8314.4
#define a b / (273.15 + 100.)

DEFINE_ADJUST(spec_grad, domain)
{

Thread *t;

cell_t c;

face_t f;

float MLFS;
float MSFS;
float UVPR;
float VVPR;

thread_loop_c (t,domain)

{
begin_c_loop_all (c,t)
{
MLFS = exp(a - b / C_T(c,t));
MSFS = (MLFS * 18./29.) / (1. - MLFS *
(l. - 18./29.));
if (C_YI(C,t,0) > MSFS) C_YI(C,t,0) = MSFS;
if (C_YI(c,t,0) <= 0.0)
C_YI(c,t,O) = 1.0e-12;
if (NULL 1: THREAD_STORAGE(t,SV4UDS_I(3)I &&
NULL I: THREAD_STORAGE(t,SV;UDS_I(0)))
{
C_UDSI(c,t,0) = C_YI(c,t,0);
C_UDSI(c,t,3) = MSFS;
}
}
end_c_loop_all (c,t)
}
thread_loop_f (t,domain)
{

if (NULL 1: THREAD_STORAGE(t,SV;UDS_I(3)) &&
NULL 1: THREAD_STORAGE(t,SV;UDS_I(0)))
{

138

 

‘11....

begin_f_loop (f,t)

float FMLFS = exp(a - b / F_T(f,t));
float FMSFS a (FMLFS * 18./29.) /
(1. - FMLFS * (1. - 13./29.));

F_UDSI(f,t,0)
F_UDSI(f,t,3)

F_YI(f,t,0);
FMSFS;

end_f_loop (f,t)

threadwloop_c (t,domain)

if (NULL 1= THREAD_STORAGE(t,SV_UDS_I(0)) &&
NULL 1= T_STORAGE_R_NV(t,SV;UDSI_G(0)))

begin_c_loop_all (c,t)
{
float diff_eff = 2.88e-5 +
(C_MU_T(c,t)/0.7);
UVPR a C_U(c,t) - diff_eff *
C_UDSI_G(c,t,O)[0]/C_UDSI(c,t,O);
VVPR a C;V(c,t) - diff_eff *
C;UDSI_G(c,t,0)[1]/C_UDSI(c,t,O);
C_UDSI(c,t,1) = C_R(c,t) * C_UDSI(c,t,O)
* UVPR;
C_UDSI(c,t,2) = C_R(c,t) * C_UDSI(c,t,0)
* VVPR;

end_c_loop_a11 (c,t)

{
}
}
}
{
{
}
}

}

thread_loop_f (t , domain)
/* assign the face value with way (2) */

{

if (NULL 1= THREAD_STORAGE(t,SV;UDS_I(1)) &&
NULL s= THREAD_STORAGE(t,SV_UDS_I(2)))

{

begin_f_loop (f,t)

{

cell_t cell a F_C0(f,t);
Thread *c_thread = THREAD_T0(t);

if (NULL 1:
T_STORAGE_R_NV(c_thread,SV_UDSI_G(1))
&& NULL I=
T_STORAGE_R_NV(c_thread,SV;UDSI_G(2)))
{

if (THREAD_ID(t)==2)

{

139

float xf[ND_ND],xc[ND_ND],
ds[ND_ND].dn;

F_CENTROID(xf,f,t):
C_CENTROID(xc,cell,c_thread);
NV_VV(ds, =, xf,-,xc);

dn = NV;MAG(ds);

F_UDSI(f,t,1) =
C_UDSI(cell,c_thread,1) - dn * C_UDSI_G(cell,c_thread,1)[1];

F_UDSI(f.t,2) =
C_UDSI(cell,c_thread,2) - dn * C_UDSI_G(cell,c_thread,2)[1];

}
else
{
F_UDSI(f,t,1) — 0.,
F_UDSI(f,t,2) = 0.,
}

}
}

end_f_loop (f,t)

}

thread_loop_c (t,domain)
/* store y mass flux in UDS(4) */

{
begin_c_loop_all (c,t)
{
if (NULL 1= THREAD_STORAGE(t,SV;UDS_I(O)) &&
NULL 12 T_STORAGE_R_NV(t,SV;UDSI_G(O)))
{
float diff_eff = 2.88e-5 +
(C_MU;T(C,t)/0.7);
C_UDSI(c,t,4) = — C_R(c,t)* diff_eff *
C_UDSI_G(c,t,0)[1];
}
end_c_loop_a11 (c,t)
}

DEFINE_PROFILE(p1ate_mf, t, position)

{

face_t f;

begin_f_loop (f,t)

{
float FMLPS = exp(a - b / F_T(f,t));
F_PROFILE(f,t,position) = (FMLFS * 18./29.) /
(1. - FMLFS * (1. — 18./29.));
}

140

 

end_f_loop (f,t)

DEFINE_SOURCE(mass_src, c, t, dS, eqn)
{
float source;
if (NULL 1: T_STORAGE_R_NV(t,SV;UDSI_G(1)) &&
NULL 1: T_STORAGE_R_NV(t,SV;UDSI_G(2)I &&
NULL 1: THREAD_STORAGE(t,SV;UDS_I(0)) &&
NULL I= THREAD_STORAGE(t,SV;UDS_I(3)I)

 

{
if (C_UDSI(c,t,0) < C_UDSI(c,t,3))
source = 0.;
else
source = C_UDSI_G(c,t,1)[0] +
C_UDSI_G(c,t,2) [1];
}
dSleqn]=0;
return source;
}
DEFINE_SOURCE(energy_src, c, t, dS, eqn)
{
float source;
if (NULL 1: T_STORAGE_R_NV(t,SV;UDSI_G(1)) &&
NULL != T_STORAGE_R_NV(t,SV;UDSI_G(2)) &&
NULL 1: THREAD_STORAGE(t,SV;UDS_I(0)) &&
NULL != THREAD_STORAGE(t,SV;UDS_I(3)I)
{
if (C_UDSI(c,t,O) < C_UDSI(c,t,3))
source = 0.;
else
source = -L * (C_UDSI_G(c,t,1)[0] +
C_UDSI_G(C,t,2) [1] );
}
d3[eqn]=0;
return source;
}

141

 

 

 

LIST OF REFERENCES

142

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

LIST OF REFERENCES

Hoke, P. B. (1998). EXPERIMENTAL MEASUREMENT OF SLIT
RESPONSE FUNCTION AND CORRECTION APPLIED TO INFRARED
THERMOGRAPHIC MEASUREMENTS, Thesis for MS degree,
Department of Mechanical Engineering, Michigan State University.

Launder, B. E. and Rodi, W. (1981 ). THE TURBULENT WALL JET, Prog.
Aerospace Sci., Vol.19, pp.81-128.

Schwarz, W. H. and Cosart, W. P. (1961). THE TWO-DIMENSIONAL
TURBULENT WALL JET, Journal of Fluid Mechanics, Vol. 10, pp.481-
495.

 

Kruke, V. and Eskinazi, S. (1964). THE WALL-JET IN A MOVING
STREAM, Journal of Fluid Mechanics Vol. 20, pp. 555-579.

Wygnanski, L, Katz, Y. and Horev, E. (1992). ON THE APPLICABILITY
OF VARIOUS SCALING LAWS TO THE TURBULENT WALL JET, J.
Fluid Mech., vol. 234, pp. 669-690.

Potter, M. and Foss, J. (1982). FLUID MECHANICS, pp. 422, Great Lakes
Press, Inc.

Eriksson, J. G., Karlsson, R. I. and Persson, J. (1998). AN
EXPERIMENTAL STUDY OF A TWO-DIMENSIONAL PLANE
TURBULENT WALL JET, Experiments in Fluids 25. pp. 50—60.

Sparrow, E. M., Minkowycz, W. J. and Saddy, M. (1967). FORCED
CONVECTION CONDENSATION IN THE PRESENCE OF
NONCONDENSABLES AND INTERFACIAL REISISTANCE, Int. J. Heat

Mass Transfer Vol. 10, pp. 1829-1845

Hijikata, K. and Mori, Y (1973). FORCED CONVECTIVE HATE
TRANSFER OF A GAS WITH CONDENSING VOPOUR AROUND A
FLAT PLATE, Heat Transfer-Jap. Res. 2, pp. 81-101.

143

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

Legay-Desesquelles, F. and Prunet-Foch, B (1986). HEAT AND MASS
TRANSFERWITH CONDENSATEION IN LAMINAR AND TURBULENT
BOUNDARY LAYERS ALONG A FLAT PLATE, Int. J. Heat Mass
Transfer, Vol. 29, NO. 1, pp.95-105.

Matuszkiewicz, A. and Vernier, Ph. (1991 ). TWO-PHASE STRUCTURE
OF THE CONDENSATION BOUNDARY LAYER WITH A NON-
CONDENSING GAS AND LIQUID DROPLETS, Int. J. Multiphase Flow,
Vol. 17, No. 2, pp. 213-225.

FLUENT INCORPORATED (1998). FLUENT 5 User’s Guide. I

Tannehill, J. C., Anderson, D. A. and Fletcher, R. H. (1997).
COMPUTATIONAL FLUID MECHANICS AND HEAT TRANSFER,
Second Edition. Taylor & Francis.

 

Bejan, A. (1993). HEAT TRANSFER. John Wiley & Sons, Inc.
Hinze, J. O. (1975). TURBULENCE, 2 ed., McGraw-Hill, New York.

Posey, S., Kremenetsky, M. and Jassim, A. (1999). PARALLEL
PROCESSING IS BRINGING CFD CLOSER TO THE FLUIDS OF THE
REAL WORLD. Mechanical Engineering, March 1999, pp. 57-59.

Bejan, A. (1995). CONCTION HEAT TRANSFER, SECOND EDITION.
John Wiley 8. Sons, Inc.

Hoke, P. B., Wang, Q., McGrath, J. J. and AbduINour B. S. (1999).
EXPERIMENTAL AND NUMERICAL STUDY OF A CONDENSING FLOW
IN A DEVELOPING WALL JET. Turbulence and Shear Flow Phenomena,
First lntemational Symposium, DoubIeTree Resort, Santa Barbara, CA.

September 12-15. (Publishing Pending)

Cengel, Y. A. and Boles, M. A. (1989). THERMODYNAMICS - AN
ENGINEERING APPROACH. McGraw-Hill, Inc.

Mabuchi, I. and Kumada, M. (1972). STUDY ON HEAT TSRTANSFER TO
TURBULENT JETS WITH ADJACENT BOUNDARIES (1 REPORT,
FLOW DEVELOPMENT AND MASS TRANSFER IN PLANE

144

TURBULENT WALL JET), Bulletin of JSME, Vol 15, No. 88, pp. 1236-
1245.

145