Simulation of Coal Ash Deposition on Modern Turbine Nozzle Guide Vanes

Thesis

Presented in Partial Fulfillment of the Requirements for the Degree of Master of Science

in the Graduate School of The Ohio State University

By:

Brett Jordan Barker, B.S.

Graduate Program in Aeronautical and Astronautical Engineering

The Ohio State University

2010

Master’s Examination Committee:

Dr. Jeffrey P. Bons, Advisor

Dr. James Gregory

Dr. Ali Ameri

Copyright by

Brett Barker

2010
 ABSTRACT

A numerical study of the physical mechanisms associated with the deposition of

coal ash particulate on a blade passage of the first stage of the GE-E3 turbine was

performed. The software package FLUENT was used to model the flow through the

turbine passage and predict the trajectory of particles injected into the flow using an

Euler-Lagrangian two phase approach. The standard k-ω turbulence model was used for

its accuracy at modeling flow very near to the wall. A critical impact velocity sticking

model and a critical viscosity sticking model were compared for their usability and

accuracy predicting deposition locations and quantities. The critical impact velocity

sticking model was chosen due to the extensive work that has been done at BYU to

calibrate the model to experimental data. The critical viscosity sticking model proved

simple to use but required additional calibration beyond the initial modeling. The

simulations were compared to experimental results obtained from the Turbine Reacting

Flow Rig (TuRFR) at The Ohio State University. Finally, modifications were made to the

hub end wall geometry to mitigate deposition growth. The hub end wall inlet angle to the

axial was decreased by 30 degrees and showed an over-all 18% decrease in deposition

mass. The area near the hub wall and pressure surface of the vane showed marked

improvement; however, there was a slight increase in deposition along the mid-span

region of the turbine vane.

ii
 ACKNOWLEDGEMENTS

I would like to thank Dr. Jeffrey Bons for his guidance and support through my Graduate

career. He provided me with opportunities to learn and encouraged me every step of the

way.

I want to give a special thanks to Dr. Ali Ameri for his patience and availability during

the long nights of computational work. Without him much of this work would not have

been possible. I would also like to thank Dr. James Gregory for being on my committee

and for being a wonderful teacher.

I also want to show my appreciation for all of my fellow Graduate students, Chris Smith,

Josh Webb, Prashanth Shankara, and Kyle Gompertz and the rest of the AARL research

members, for their support and help through this process.

iii
 VITA

May 30, 1986 ………………………………….Born in Sidney, Ohio

June 2008 ………………………………………B.S. Aerospace Engineering,

The Ohio State University

2008 – Present ………………………………Graduate Research Associate,

Department of Aerospace Engineering

The Ohio State University

PUBLICATIONS

1. Lewis, S., Barker, B., Bons, J.P., Ai, W., and Fletcher, T.H. "Film Cooling

Effectiveness and Heat Transfer Near Deposit-Laden Film Holes." GT2009-

59567. Presented at the IGTI 2009 conference and recommended for publication

in the Journal of Turbomachinery publication, 2009

2. Smith, C., Barker, B., Clum, C., Bons, J.P. “Deposition in a Turbine Cascade with

Combusting Flow.”Presented at the ASME Turbo Expo 2010: Power for Land,

Sea, and Air (2010)

iv
 AWARDS

IGTI/ASME Heat Transfer Committee 2009 Best Paper Award for “Film Cooling

Effectiveness and Heat Transfer Near Deposit-Laden Film Holes” S. Lewis, B. Barker,

J.P. Bons, W. Ai, T.H. Fletcher. GT2009-59567

FIELDS OF STUDY

Major Fields of Study: B.S. Aerospace Engineering

M.S. Aerospace Engineering

v
 TABLE OF CONTENTS

ABSTRACT........................................................................................................................ii

ACKOWLEDGEMENTS..................................................................................................iii

VITA...................................................................................................................................iv

FIELDS OF STUDY...........................................................................................................v

TABLE OF CONTENTS....................................................................................................vi

LIST OF FIGURES.............................................................................................................x

LIST OF TABLES...........................................................................................................xiii

LIST OF EQUATIONS....................................................................................................xiv

NOMENCLATURE........................................................................................................xvii

CHAPTER 1........................................................................................................................1

INTRODUCTION...............................................................................................................1

1.1 Background........................................................................................................1

1.2 Literature Review...............................................................................................2

1.3 Other Deposition Research................................................................................3

CHAPTER 2........................................................................................................................5

FLUID MODELING...........................................................................................................5

2.1 Fluid Modeling Package....................................................................................5

2.1.1 Grid Generation..............................................................................................5

2.2 Fluid Equations of Motion.................................................................................9

vi
 2.3 Turbulence Modeling.......................................................................................11

2.3.1 Near Wall Flow.............................................................................................12

2.4 Boundary Conditions.......................................................................................13

2.4.1 Inlet Boundary Condition.............................................................................14

2.4.2 Outlet Boundary Condition...........................................................................15

2.4.3 Wall Boundary Condition.............................................................................16

2.4.4 Periodic Boundary Condition.......................................................................16

2.5 Experimental Description................................................................................17

2.6 Operating Conditions.......................................................................................19

2.7 Flow Solution Results......................................................................................20

CHAPTER 3......................................................................................................................23

PARTICULATE MOTION...............................................................................................23

3.1 Particle Trajectory Modeling...........................................................................23

3.1.1 Eulerian-Lagrangian Model..........................................................................23

3.1.2 Particle Forces...............................................................................................26

3.2 Injection Properties..........................................................................................31

3.3 Random Walk Model.......................................................................................32

3.3.1 Perturbation Velocity Corrections................................................................34

3.4 Particle Material Properties..............................................................................35

3.4.1 Ash Composition..........................................................................................36

3.4.2 Size Distribution...........................................................................................37

3.5 Particle Trajectories.........................................................................................40

vii
CHAPTER 4......................................................................................................................42

STICKING MODELING...................................................................................................42

4.1 Impact, Sticking, Capture Efficiency Definitions............................................42

4.2 Wall-Particle Interaction..................................................................................43

4.3 Critical Impact Velocity Sticking Model.........................................................44

4.3.1 Particle Impact Properties.............................................................................45

4.4 Critical Viscosity Sticking Model....................................................................47

4.4.1 Ash Compensation Model.............................................................................48

4.4.2 Viscosity-Temperature Relationship............................................................49

4.4.3 Relationship Limits.......................................................................................51

4.4.4 Critical Viscosity Sticking Criteria...............................................................53

4.5 Pre-Existing Deposition Factor........................................................................55

4.6 Detachment Model...........................................................................................55

4.6.1 Individual Particle Detachment....................................................................56

4.6.2 Large Scale Detachment...............................................................................59

4.7 Sticking Model Comparisons...........................................................................61

4.8 Experimental Comparison...............................................................................67

4.9 Particle Size Deposition Distribution ..............................................................72

CHAPTER 5......................................................................................................................75

END WALL MODIFICATION........................................................................................75

5.1 Reasoning for end wall modification...............................................................75

5.2 Method of adjustment......................................................................................75

viii
 5.3 Exit and inlet profile changes..........................................................................78

5.4 Deposition augmentation.................................................................................80

5.5 Recommendations for improvements..............................................................82

CHAPTER 6......................................................................................................................84

CONCLUSION AND FUTURE WORK..........................................................................84

6.1 Conclusion.......................................................................................................84

6.2 Future Work.....................................................................................................85

REFERENCES..................................................................................................................87

APPENDIX........................................................................................................................93

A.1 Critical Impact Velocity Sticking Model UDF...............................................93

A.2 Critical Viscosity Sticking Model UDF..........................................................99

A.3 Random Walk Correction UDF....................................................................105

A.4 Local Flow Properties Injection UDF...........................................................115

ix
 LIST OF FIGURES

Figure 1 – VKI Grid............................................................................................................6

Figure 2 – E3 Annulus (Colored by Radial Position) .........................................................7

Figure 3 – E3 geometry.......................................................................................................7

Figure 4 – Mid-span of E3 Turbine vane grid....................................................................8

Figure 5 – Modified E3 geometry.......................................................................................9

Figure 6 – Universal laws of the wall, Schlichting and Gersten [35] ...............................12

Figure 7 – Fitted capture efficiencies obtained from 2-D CFD modeling versus measured

values [3] ...........................................................................................................................15

Figure 8 – Image of the TuRFR.........................................................................................17

Figure 9 – Image of vane holder........................................................................................18

Figure 10 – Mach distribution at mid-span........................................................................21

Figure 11 – Pressure distribution at mid-span...................................................................21

Figure 12 – Temperature distribution at mid-span............................................................22

Figure 13 – Total Pressure at Exit Plane............................................................................22

Figure 14 – Classification of flow regimes for gas-solid flows [17] ................................24

Figure 15 – Particle injection locations.............................................................................31

Figure 16 – RMS values of velocity in the boundary layer [13] ......................................34

Figure 17 – Particle size distribution results from Coulter counter for Bituminous coal fly

ash [38] ..............................................................................................................................38

Figure 18 – Size Distribution Modeling............................................................................39

x
Figure 19 – 1 μm diameter particle trajectory...................................................................40

Figure 20 – 20 μm diameter particle trajectory.................................................................41

Figure 21 – 50 μm diameter particle trajectory.................................................................41

Figure 22 – dimensionless sticking forces versus particle diameter [6] ...........................43

Figure 23 – Transition modified deposition model using critical viscosity approach

[41].....................................................................................................................................48

Figure 24 – Viscosity – Temperature relationship for various ash types [41] ..................52

Figure 25 – Sawtooth wave and random number generated..............................................54

Figure 26 – Geometric features of a spherical particle stuck to a smooth surface [16].....55

Figure 27 – Cone and cylinder modeling...........................................................................59

Figure 28 – Coefficient of drag for a cylinder for vs. Reynolds numbers.........................60

Figure 29 – Large Scale Deposition Removal Comparison...............................................61

Figure 30 – Number of Captured Particles along Axial Chord.........................................63

Figure 31 – Mass % of Captured Particles along Axial Chord..........................................64

Figure 32 – VKI Surface Temperature along Pressure Surface.........................................65

Figure 33 – Sticking Model Comparison of Deposition Mass % locations.......................67

Figure 34 – Low Temperature Comparison between results.............................................68

Figure 35 – Medium Temperature Comparison between results.......................................69

Figure 36 – Time lapse of deposition growth in TuRFR facility.......................................70

Figure 37 – Turbine vane after cool down.........................................................................71

Figure 38 – Turbine vane after cool down with suspended deposition.............................72

Figure 39 – Deposition locations for various particle diameters.......................................73

xi
Figure 40 – Impact and Deposit Mass % locations............................................................76

Figure 41 – Flow visualization of horseshoe vortex for E3 (Vortex magnified) ..............77

Figure 42 – Flow Visualization of horseshoe vortex for E3M (Vortex Magnified) ..........78

Figure 43 – Inlet Mach.......................................................................................................79

Figure 44 – Exit Mach.......................................................................................................79

Figure 45 – Exit flow angle...............................................................................................80

Figure 46 – Impact Mass % Locations..............................................................................81

Figure 47 – Deposition Mass % Locations........................................................................82

xii
 LIST OF TABLES

Table 1 – Experimental Operating Conditions..................................................................19

Table 2 – Simulation Operating Conditions......................................................................20

Table 3 – Particle Forces [37] ...........................................................................................30

Table 4 – Particle Material Properties BYU......................................................................36

Table 5 – Coal fly ash composition (%w.t.) [41] ..............................................................37

Table 6 – Particle size distribution for FLUENT simulation.............................................38

Table 7 – Fitted Values of Coefficients for Compositional Dependence of B[36] ..........50

Table 8 – Coal Ash Chemical Composition used at BYU.................................................51

Table 9 – Viscosity-Temperature relationship constants...................................................52

Table 10 – VKI Sticking Model Comparison Results.......................................................62

Table 11 – E3 Sticking Model Comparison Results..........................................................66

Table 12 – Particle Size Comparison Results....................................................................73

Table 13 – E3 Sticking Model Comparison Results..........................................................81

xiii
 LIST OF EQUATIONS

Equation 1 - Mass Conservation Equation..........................................................................9

Equation 2 - Momentum Conservation Equation................................................................9

Equation 3 - Stress Tensor.................................................................................................10

Equation 4 - Energy Conservation Equation.....................................................................10

Equation 5 - Total Energy..................................................................................................10

Equation 6,7 - Transport Equations for Standard k-ω Turbulence Model.........................11

Equation 8 - Non-Dimensional Wall Distance Definition.................................................13

Equation 9 - Non-Dimensional Wall Velocity Definition.................................................13

Equation 10 - Friction Velocity Definition........................................................................13

Equation 11 - Isentropic Pressure Relationship.................................................................14

Equation 12 - Turbulence Intensity....................................................................................15

Equation 13 - Particle Force Balance for Unit Mass.........................................................26

Equation 14 - Spherical Particle Drag Force.....................................................................26

Equation 15 - Summation of Particle Forces.....................................................................27

Equation 16 - Saffman Lift Force......................................................................................28

Equation 17 - General form of Saffman Lift Force...........................................................28

Equation 18 - Deformation Tensor....................................................................................28

Equation 19 - Knudsen Number........................................................................................29

Equation 20 - Mean Free Path...........................................................................................29

Equation 21 - Average Molecule Speed............................................................................29

xiv
Equation 22 - Local Flow Velocity....................................................................................32

Equation 23 - Instantaneous Velocity Definition...............................................................32

Equation 24 - RMS Values for Isotropic Flow..................................................................32

Equation 25 - Time Integral...............................................................................................33

Equation 26 - Lagrangian Integral Time............................................................................33

Equation 27 - Time Scale for Eddy....................................................................................33

Equation 28 - Time for Particle to Traverse Eddy.............................................................33

Equation 29 - Rosin-Rammler Distribution.......................................................................39

Equation 30 - Impact Efficiency Definition.......................................................................42

Equation 31 - Sticking Efficiency Definition....................................................................42

Equation 32 - Capture Efficiency Definition.....................................................................42

Equation 33 - Van Der Waals Force..................................................................................44

Equation 34 - Modified Van Der Waals Force..................................................................44

Equation 35 - Critical Impact Velocity..............................................................................45

Equation 36, 37, 38 - Composite Young’s Modulus.........................................................45

Equation 39 - Third-Order Particle Young’s Modulus Relation.......................................46

Equation 40 - First-Order Particle Young’s Modulus Relation.........................................46

Equation 41 - Exponential Particle Young’s Modulus Relation........................................46

Equation 42 - Modified Exponential Particle Young’s Modulus Relation........................46

Equation 43 - Average Temperature..................................................................................46

Equation 44 - Viscosity Sticking Probability.....................................................................47

Equation 45 - Non-Bridging Oxygens to Tetrahedral Oxygens Ratio (NBO/T) ..............49

xv
Equation 46 - Viscosity Temperature Relationship...........................................................49

Equation 47, 48, 49, 50, 51, 52 - Viscosity Temperature Relationship Coefficient

Equations............................................................................................................................49

Equation 53 - Sticking Probability Exponential Fit...........................................................52

Equation 54 - Sawtooth Function......................................................................................53

Equation 55 - Detachment Model Summation of Moments..............................................56

Equation 56 - Detachment Model Summation of Moments Reduced..............................57

Equation 57 - Attached Spherical Particle Drag................................................................57

Equation 58 - Cunningham Correction Factor...................................................................57

Equation 59 - Stokes Coefficient of Drag Equation..........................................................57

Equation 60 - Modified Attached Spherical Particle Drag................................................58

Equation 61 - Attached Spherical Particle Contact Radius...............................................58

Equation 62 - Composite Young’s Modulus.....................................................................58

Equation 63 - Critical Wall Shear Velocity.......................................................................58

Equation 64 - Large Scale Detachment Moment Balance.................................................60

Equation 65 - Large Scale Sticking Force.........................................................................60

xvi
 NOMENCLATURE

A Hamaker Constant

AL Low Temperature Coefficient

AH High Temperature Coefficient

Acp Cross Sectional Area of Particle

BL Low Temperature Coefficient

BH High Temperature Coefficient

CD Coefficient of Drag

CP Heat Capacity at Constant Pressure

Cu Cunningham Correction Factor

CFD Computational Fluid Dynamics

Dp Diameter of particle

DPM Discrete Particle Model

E Total Energy

E Composite Young’s Modulus

Ep Particle Young’s Modulus

Es Surface Young’s Modulus

External Force Vector

FD Drag Force of Particle

Van der Waals Sticking Force

FL Fluid Lift Force

xvii
Gk Generation of Turbulent Kinetic Energy

Gω Generation of ω

I Unit Tensor

I Turbulence Intensity

Diffusion Flow of the Specifies j

Kc Composite Young’s Modulus

Kr Local Velocity Gradient

Knudsen Number

Le Length Scale of Eddy

M Mach Number

P Static Pressure

Po Total Pressure

Ps(T) Probability of Sticking Function

R Specific Gas Constant

ReDH Reynolds Number based on Hydraulic Diameter

Rep Reynolds Number based on Particle Properties

RNG Random Number Generator

S Ratio of Particle Density to Fluid Density

Sh Heat Source

Sk Kinetic Energy Source

Sm Mass Source

Sω ω Source

xviii
T Static Temperature

TAVG Average Temperature

TE Time Scale for Eddy

Tg Temperature of Gas

Time Integral

TL Lagrangian Integral Time

Tp Temperature of Particle

Ts Critical Sticking Temperature

TADF Turbine Accelerated Deposition Facility

TuRFR Turbine Reaction Flow Rig

UDF User Defined Function

Vcr Critical Impact Velocity

WA Work of Sticking

Yk Dissipation of k due to Turbulence

Yω Dissipation of ω due to Turbulence

a Random Function Period

Average Molecule Speed

dij Deformation Tensor

dm Mean Particle Diameter

dp Diameter of Particle

f FD Correction Factor for Near Wall

Body Force Vector

xix
h Enthalpy

h Cone Height

hc Height of Cone Center

k Turbulent Kinetic Energy

keff Effective Conductivity

n Spread Parameter

t Time

tcross Time for Particle to Cross Eddy

Velocity Vector

u+ Non-Dimensional Wall Velocity

u* Friction Velocity

Instantaneous Velocity

uc Fluid Velocity at Particle Center

up Particle Velocity

uτc Critical Wall Shear Velocity

y+ Non-Dimensional Wall Distance

z Distance between Particle and Surface

α Relative Approach between Particle and Surface

α2 Volume Fraction of the Dispersed Phase

Ratio of Specific Heats

Normal Distribution Function

Mean Free Path

xx
μ Molecular Viscosity

μcrit Critical Sticking Viscosity

μTp Viscosity of Particle at Tp

ν Kinematic Viscosity

νp Poisson Ratio for the Particle

νs Poisson Ratio for the Surface

Fluid Density

p Particle Density

Stress Tensor

τx12 Particle Relaxation Time

τk Kolmogorov Time Scale

τt1 Lagrangian Integral Time Scale

τw Wall Shear Stress

ω Specific Dissipation Rate

xxi
 CHAPTER 1: INTRODUCTION

1.1 Background

Traditional turbine engine fuel is becoming less available due to the political and

environmental concerns of oil. To meet with the demand, modern engineers are using

fuels derived from new sources. The direct liquefaction and gasification of coal has

provided an alternative fuel for jet engines and power generators. Coal is widely available

and is inexpensive compared to the highly refined oil based fuels. There are some

drawbacks to using coal derived fuels, however. Coal has negative impacts on the

environment and produces more pollutants during combustion. The products of

combustion contain airborne particulates known as fly ash. This fly ash deposits on the

internal walls of the engine and degrades the system over time. The deposition negatively

affects the power generation of the engine and surface roughness resulting from such

deposition augments the heat transfer of the high temperature fluid to the surfaces.

Removal of the deposits is expensive and time consuming as the afflicted parts need to be

removed for repair. The parts require extensive cleaning to remove the deposition.

Understanding how the particulate adheres to the walls and under what conditions is

crucial to reducing the amount of deposition. Reducing the deposition will reduce the

time spent cleaning the engine and prolong the life of the hardware.

1
1.2 Literature Review

There have been many studies looking into the causes and effects of deposition on

turbine hardware. Small scale deposition can increase the surface roughness. Bons

concluded that increases in surface roughness can decrease turbine performance [8].

Abuaf supported the finding that losses are associated with increased roughness and

found that the heat transfer was increased with additional roughness [1]. Large scale

deposition was found to be more damaging to turbine hardware. Kim et. al. conducted

experiments exploring the effects of volcanic ash on turbine hardware [24]. They found

that deposition can clog film cooling holes and can lead to failure of the turbine. Dunn

et. al. also studied the effects of volcanic ash and found that damage due to deposition

was related to the turbine inlet temperature, concentration of particulate and the material

properties of the volcanic ash [14]. Dunn et. al. confirmed that deposition can clog film

cooling holes and that deposition is a major issue for modern engines with high

combustor temperatures. Dunn et. al. also noted that aircraft engines that ingest particle

laden flow have increased difficulty in restarting and that engines exposed to sufficient

levels of particulate can be damaged beyond repair. Sundaram et. al. studied the effects of

deposition on film cooling and found that film cooling effectiveness degrades as

deposition forms near and in film cooling holes [40]. Lewis et. al. found that the location

of deposition around film cooling holes can affect the heat transfer [28]. Deposition that

forms between and downstream of film cooling holes causes the high temperature free

stream to flow into the deposition valleys and increases heat transfer into the surface. In

the rare cases where deposition forms exclusively upstream of the film cooling holes, a

2
cooling benefit was observed. Lawsen et. al. conducted experiments in a low speed wind

tunnel using low melt wax to investigate the sticking mechanism and formation of

deposits on a flat plate [26]. Jensen et. al. and Crosby et. al. constructed the Turbine

Accelerated Deposition Facility (TADF) and observed the growth of deposits on one inch

coupons[11,21]. This facility was operated at actual engine temperatures and used coal

fly ash for particulate. Smith et. al. constructed the Turbine Reaction Flow Rig (TuRFR)

and tested actual turbine hardware at engine temperatures and velocities. Deposition from

coal fly ash was found to form on the pressure surface and leading edge of the turbine

vanes. The negative effects of deposition have demanded that the mechanisms of

deposition formation be better understood yet the difficulties in operating experiments at

engine operating conditions has supported the need for computer modeling of deposition.

1.3 Other Deposition Research

Numerical models of deposition have been the focus of much research due to the

freedoms of computational fluid dynamic (CFD) research. Hossain et. al. developed a

model to track particles through piping with bends [20]. They found that small diameter

particles were less likely to deposit on the walls because they had a higher tendency to

follow the flow. Larger particles with higher inertia were more likely to deposit

downstream of the bends. Longmire looked into the particle trajectories and developed

corrections for particle dispersion when predicting particle flow near surfaces [29]. El-

Batsh et. al. developed a sticking model to quantify the criteria for particle sticking and

detachment [15]. They tested a 2D turbine vane cascade to study the locations where

3
deposition was most likely. Ai et. al. modified the sticking model developed by El-Batsh

and calibrated the sticking model to match experimental results from the TADF [2]. Tafti

et. al. developed a sticking model using the coal ash composition to determine a sticking

probability based on the viscosity-temperature relationship [41]. The focus of the study

presented in this thesis is using the sticking models developed by El-Batsh et. al., Ai et.

al. and Tafti et. al. on a 3D turbine vane with end walls and determining the effects of end

walls on deposition. This study explores the 3D effects on particle trajectories, as well as

the capture efficiencies for particles of varying size. This study also develops

mechanisms for deposition accumulation and removal for transient studies.

4
 CHAPTER 2: FLUID MODELING

2.1 Fluid modeling package

The understanding of deposition that was found from experiments was

incorporated into a numerical model. The modeling process is done in two steps. The first

step is the carrier phase simulation step and the second is the discrete particle modeling

(DPM). The modeling was accomplished using the software package FLUENT 6.3.26.

FLUENT is a commercially available computational fluid dynamics (CFD) solver

program and was chosen for its prevalence in industry. The modeling in this study can be

applied using other CFD models as well. FLUENT solves flow equations for finite-

volumes and is used to study the effects of fluid dynamic forces and heat transfer on

aerodynamic bodies. FLUENT was used at the Ohio Super Computer center (OSC).

FLUENT can be modified for specific applications using User Defined Functions

(UDF). UDF’s can be written to adjust certain parameters of the modeling to meet the

requirements of the simulation. Three UDF’s were used in post processing for this study

and are described later.

2.1.1 Grid generation

Definition of the computational domain and production of the required grid for

use in FLUENT was achieved using a grid generator. The number and shape of these

finite volumes affects the quality and accuracy of the solution. There are several

5
programs used for grid generation with each having their strengths. The grid generator

used for this study was Grid Pro. Grid Pro is a commercially available grid generator that

generates the grid automatically from a user defined topology. Surfaces are provided

analytically or from a CAD model.

Three geometries were used in this study. The first is the VKI turbine vane. This

geometry is two dimensional and is an individual blade from a cascade. The grid is

shown below in Fig. 1.

Figure 1 - VKI Grid

This grid has 5168 cells with the near-wall cells having y+ values around 1.

The second geometry is from the first stage of the General Electric Energy

Efficient Engine (GE-E3) turbine. This geometry is three dimensional and is an individual

vane from the annular stage. The geometry includes the casing and hub end walls.

6
Figure 2 – E3 Annulus (Colored by Radial Position)

Figure 3 – E3 geometry

7
The grid for the individual passage has 568512 cells and is show in Fig. 3. A 2D cross

section of the mid-span is shown in Fig. 4.

Figure 4 – Mid-span of E3 Turbine vane grid

The third geometry in this study is the E3 vane, but with a modified hub end wall (E3M).

The end wall was modified to study the effects it has on the deposition growth. The end

wall surface is a surface of revolution formed by a hub contour (Generatrix) rotated

around the axis of the machine. This hub contour curve was altered as shown in Fig. 5.

8
 (a) (b)

Figure 5 – Modified and unmodified hub end wall contour of the E3 geometry, (a)

Analytical function of hub end wall, (b) resultant modified surface

The inlet angle of the hub was changed approximately 30 degrees and the total inlet area

decreased by about 10%.

2.2 Fluid equations of motion

FLUENT solves the equations of fluid motion to determine the flow properties.

The equations of motion are the conservation of mass and momentum and are known as

the Navier-Stokes equations. The differential forms of the mass and momentum

conservation equations are defined below [18]

[1]

[2]

9
where

[3]

where ρ is the fluid density, is the velocity vector, Sm is mass due to a dispersed second

phase (e.g. vaporization of liquid droplets) in the finite volume, is a body force, and

is the external body force term. also contains other source terms dependent on the

model and user-defined sources. The term is a stress tensor where µ is the molecular

viscosity and I is the unit tensor. The momentum equation is a vector equation and has an

equation for each spatial dimension of the case. These equations are solved for each finite

element and their information is passed on to the next element.

For compressible flow and flows with heat transfer, the energy equation is

included. The energy equation is defined below

[4]

where

[5]

where keff is the effective conductivity, is the diffusion flow of the species j, h is the

enthalpy, and Sh is a source of additional heat in the flow (e.g. exothermic chemical

reactions).

10
2.3 Turbulence Modeling

There are several turbulence models that are used in numerical studies; however

only a few have been used with particle transport and deposition studies. Common

turbulence models are of the k-ε family, where k is the turbulent kinetic energy and ε is

the dissipation rate of kinetic energy. A variation of this model is the k-ω turbulence

model which uses ω, the turbulent frequency defined as ω=ε/k. El-Batsh used the

Standard k-ε and the RNG k-ε with Standard wall functions and two-layer zonal model

for studying deposition [16]. This choice was made due to the limitations of turbulence

models available in the earlier versions of FLUENT. Ai et al. used the standard k-ω with

RANS to determine the flow field and heat transfer in their study of deposition on turbine

coupons [2]. k-ω was chosen to avoid using near-wall functions that could incorrectly

model particle trajectories. Ajersh et. al. found that the k-ω turbulence model gave better

predictions of the near-wall flow structures compared to the k-ε model [4]. Wilcox found

that k-ω also handles non-zero pressure gradients better than the k-ε model [47]. For

these reasons, k-ω was used as the turbulence model. The k-ω model is defined by the

equations:

[6]

[7]

11
where Gk represents the generation of turbulent kinetic energy due to mean velocity

gradients, Gω is the generation of ω, Yk and Yω are the dissipation of k and ω due to

turbulence and Sk and Sω are user-defined source terms.

2.3.1 Near wall flow

Flow near a solid surface will form a boundary layer to separate the free stream

flow from the no-slip velocity of the surface. The flow in the boundary layers has been

studied extensively over the years and there has been much development in predicting the

velocity fields. For turbulent boundary layers the boundary layer can be described as two

distinct regions; the viscous sub-layer and the fully turbulent layer. The limits of these

regions are shown in Fig. 6,

Figure 6 - Universal laws of the wall, Schlichting and Gersten [35]

12
The dimensionless terms, u+ and y+, are defined as:

[8]

[9]

where,

[10]

u* is the friction velocity, ν is kinematic viscosity, and τw is the wall shear stress. Several

numerical functions have been developed to describe the two regions of the turbulent

layer for flow over flat plates, channels and tubes.

2.4 Boundary Conditions

The boundary conditions are set to create the flow conditions in the study. The

inlet boundary condition provides the inlet flow. In the case of the turbine vanes, this will

be the exhaust from the combustor. The outlet boundary condition is the exhaust for the

flow through the domain. The outlet for the turbine vanes will normally be the inlet to the

next stage; however, this study has the flow exiting to atmospheric conditions. The wall

boundary conditions are the geometry of the turbine vane. To reduce the computational

load, periodic boundary conditions were used. The periodic boundaries are set in between

adjacent turbine vanes. The periodic boundary conditions as used enforces that the flow

in the annulus between vanes are identical.

13
2.4.1 Inlet boundary conditions

The inlet boundary was set as a “Pressure Inlet” for this study. The

pressure inlet defines the total and static pressure of the flow at the inlet boundary. The

inlet Mach number of the flow is determined by the relation

[11]

The static pressure at the inlet was set to achieve an inlet Mach number of around 0.1.

This value is typical for the exit Mach of a combustor [5]. The total pressure was set to

the operating conditions of the experiment as described below.

The inlet total temperature was set to match the temperatures experienced in

actual engines. The free stream temperature from experimental work done by Crosby et.

al. has been shown to have a significant effect on the accumulation of deposition [11]. Ai

et. al. calibrated their numerical model to match capture efficiencies for various free

stream temperatures as shown in Fig. 7 [3].

14
Figure 7 – Fitted capture efficiencies obtained from 2-D CFD modeling versus measured

experimental values [3]

The combustion process is highly turbulent and can affect the flow turbulence

level entering the turbine stage. The inlet turbulence intensity was set using the relation

for duct flow[18]:

[12]

where ReDH is the Reynolds numbers based on the hydraulic diameter which is set to the

airfoil chord.

2.4.2 Outlet boundary conditions

The outlet boundary was set as a “pressure outlet”. The pressure outlet specifies

the static pressure at the outlet boundary. The Mach number at the exit for the bulk flow

15
can be approximated using equation 11 with the inlet total pressure, outlet static pressure,

and the specific heat ratio, γ, as 1.36 for the high temperatures.

2.4.3 Wall boundary conditions

The wall boundaries guide the flow through the domain and represent the vane

geometry. The walls are smooth and use the “no-slip” condition for near-wall flow. The

wall thermal properties were tested under two conditions. The first condition was a

constant wall temperature. This temperature was set to 90% of the total inlet temperature.

This model is less computationally sensitive but does not describe the physical conditions

experienced with flows with varying heat transfer. The second condition was to set the

wall heat flux to zero. The adiabatic conditions would allow the wall temperature to vary

along the wall. This condition is computationally more difficult to converge but provides

a more accurate model of the vane surface temperatures. Both of these conditions are for

a non-cooled turbine vane. The adiabatic wall case was found to better approximate the

surface temperatures and is presented in this thesis.

2.4.4 Periodic boundary conditions

The periodic boundary condition was used to model the passage flow as a

repeating pattern of identical flow fields for each vane. This is valid when used for a

single row of vanes with circumferentially uniform (or periodic) upstream boundary

conditions. With this, a single turbine vane was modeled.

16
2.5 Experimental Description

This study modeled the Turbine Reacting Flow Rig (TuRFR) as described by

[10]. Fig. 8 shows a schematic of the TuRFR facility. The facility tests turbine hardware

at temperatures and velocities typical in engines.

Figure 8 – Image of the TuRFR

Air flow is brought in through the combustor where fuel and particulate are added. The

fuel is combusted and the hot gases accelerate through the cone. The flow is stabilized in

the equilibration tube. The transition piece changes from circular to rectangular. The flow

is then sent through a rectangular to annular transition. Two turbine vane doublets are

located at the end of this transition. The exit of the turbine vanes open to the atmosphere.
17
The vanes are held in place in the vane holder as shown in Fig. 9. The vanes can be

supplied with coolant for film cooling.

Figure 9 – Image of vane holder

The operation of the facility is described by Smith et. al.[38].

The turbine vanes used in this study are part of the first stage turbine in the

General Electric CFM-56 engine. The CFM-56 turbine vanes feature cooling holes

located near the leading edge and along the pressure surface. The CFM-56 geometry was

unavailable for CFD simulation so the GE-E3 turbine vane was used instead. Due to this

reason, this study will focus on the trends of deposition rather than the specific deposition

locations on the turbine vanes.

18
2.6 Operating conditions

The TuRFR was operated at three temperatures and one mass flow. The operating

conditions of these tests are shown in table 1.

Table 1 – Experimental Operating Conditions

Case Can Inlet Recovery Inlet Exit Static Exit Reynolds

Pressure Temperature Mach Pressure [Pa] Number

[Pa] [K]

CFM56 -1 206842 1284 0.092 98595 360,500

CFM56 -2 144789 1314 0.098 98595 351,900

CFM56-3 193053 1338 0.093 98595 345,500

The simulations utilized the same operating pressures and mass flow and used three

temperatures as shown in Table 2.

19
 Table 2 – Simulation Operating Conditions

Case Inlet Total Inlet Total Inlet Inlet Exit Static Exit Static

Pressure Temperature Mach Turbulence Pressure Temperature

[Pa] [K] Intensity % [Pa] [K]

VKI - 1 179263 1283 0.095 4.000 98595 300

VKI - 2 179263 1338 0.095 4.000 98595 300

VKI - 3 179263 1477 0.095 4.000 98595 300

E3 - 1 179263 1283 0.095 4.335 98595 300

E3 - 2 179263 1338 0.095 4.360 98595 300

E3 - 3 179263 1477 0.095 4.418 98595 300

E3M - 1 179263 1283 0.095 4.335 98595 300

E3M - 2 179263 1338 0.095 4.360 98595 300

E3M - 3 179263 1477 0.095 4.419 98595 300

The results of the flow solution and deposition locations will be explored below.

2.7 Flow solution results

Below are the flow properties for each geometry at mid-span. Fig. 10 shows the

Mach distribution. The mach distribution for E3 and E3 M are very similar despite the

difference in end wall geometry. The static pressure distribution is shown in Fig. 11 and

the static temperature distribution is shown in Fig. 12.

20
 (a) (b) (c)

Figure 10 – Mach distribution at mid-span, (a) VKI, (b) E3, (c) E3M

[Pa]

(a) (b) (c)

Figure 11 – Pressure distribution at mid-span, (a) VKI, (b) E3, (c) E3M

21
 [K]

(a) (b) (c)

Figure 12 – Temperature distribution at mid-span, (a) VKI, (b) E3, (c) E3M

Turbine blade rows are designed so that each row provides the following row the correct

flow angles and magnitudes. Another important parameter is the total pressure loss across

the blade row. The total pressure at the exit plane for the E3 and E3M cases are shown in

Fig. 13. The E3 and E3M have very similar exit total pressure maps and the modification

to the end wall does not have a major effect on the over-all performance of the stage. The

mass averaged total pressure at the exit was 166268 Pa and 167281 Pa for the E3 and

E3M respectively.

Pressure Side

[Pa]
C
H a
u s
b e

W W
a a
l l
l l
Suction Side
(a) (b)

Figure 13 – Total Pressure at Exit Plane, (a) E3, (b) E3M

22
 CHAPTER 3: PARTICULATE MOTION

3.1 Particle trajectory modeling

The particle trajectory modeling was done using the DPM module in FLUENT.

The DPM uses the Euler-Lagrange approach to track the particle trajectory. The DPM

module allows for flow induced forces on particles to be computed as they determine the

flow path. This method is described below.

3.1.1 Eulerian-Lagrangian model

The Eulerian-Lagrangian model splits the particulate modeling into two phases.

The first phase (Eulerian) is solving the flow solution absent of particulate. The flow

solution phase is treated as a continuum and the Navier-Stokes equations are solved

throughout the fluid domain. The second phase is tracking a large number of dispersed

particles traveling through the flow. The trajectory of each particle is predicted and stored

for analysis. Zhang & Chen found this method to provide accurate results that are easy to

interpret but are computationally expensive [48].

Goesbet et. al. showed in their studies that there can be three types of coupling for

solid phase particles in turbulent flows [19]. The first coupling is One-Way; the

turbulence in the flow affects the particle trajectories and dispersion. The second

coupling is Two-Way and has the particles affecting the turbulence. The third coupling is

23
Four-Way and includes the particles interacting with each other. Fig. 14 shows the

various gas-solid flows where each coupling method is appropriate.

Figure 14 – Classification of flow regimes for gas-solid flows [17]

α2 is the volume fraction of the dispersed phase, τk, is the Kolmogorov time scale in

seconds, τt1 is the Lagrangian Integral time scale in seconds and τx12 is the particle

relaxation time.

Typically, the volumetric fraction of particles in an engine is around 0.1 parts per

million by weight (ppmw) [46]. Coal fly ash particulate has a much higher density than

the fluid and the number of particles is low enough for inter-particle collisions to be

neglected. Kulick et. al. and Kaftori et. al. have shown that for volume fractions less than

10E-6 the turbulence modification to particles in the flow is negligible[25]. For these

reasons, One-Way coupling is used in this study and the number of particles injected into

the flow is kept low.

24
 The DPM tracking in FLUENT is used for modeling particle trajectories under

One-Way coupling. The particle is tracked from the injection point until it meets with one

of three fates. The first fate is escaped and means the particle exited the domain. For

accurate results, the particle should escape through the outlet boundary. Particles can

escape through the inlet boundary as well but this has been avoided by extending the inlet

boundary far enough upstream to prevent particles bouncing past the inlet.

The second fate is abort and means the particle trajectory has been canceled. This

signifies that the particle was trapped by the flow. The UDF for the sticking models

assigns the abort fate to particles that deposit on the vane surface. Particles that do not

deposit return to the flow and continue their trajectory. Particles that stick can be

removed by the flow as described in Chapter 4 Section 6. Their trajectories were not

computed after leaving the surface.

The third fate is incomplete and means the particle did not reach the other two

fates before reaching the specified number of time steps. This is encountered when

particles are stuck in circulatory flow and cannot escape the domain. FLUENT has a

maximum number of steps a particle can take and the incomplete fate is used to prevent

an infinite loop. For accuracy the number of incomplete particles should be kept less

than 1% of the injected particles.

25
3.1.2 Particle forces

FLUENT predicts the trajectory of a discrete phase particle (or droplet or bubble)

by integrating the force balance on the particle. The particle inertia is balanced with the

forces acting on the particle and can be written per unit mass in the x direction as:

[13]

where u is the fluid phase velocity, up is the particle velocity, μ is the molecular viscosity

of the fluid, ρ is the fluid density, ρp is the particle material density, FD is the drag force

per unit particle mass, and dp is the particle diameter. FD is defined as:

[14]

where CD is the coefficient of drag for the particle. The definition for FD is valid for

Reynolds numbers less than 1. The second term on the right hand side in equation 13 is

the acceleration of fluid around the particle. Fx is an additional acceleration term that

appears under specific situations. These acceleration terms are described in further detail

below.

In this study, the following assumptions were made regarding the interaction

between the fluid phase and the particles. The particles are assumed to be perfect spheres

that are considered as points located at the center of sphere. The various forces acting on

the particle were described by Rudinger as [33]:

26
 [15]

These forces have been discussed extensively by El-Batsh et. al. and many others so only

a brief description is provided here.

The drag force is the Stokes drag law and is dependent on the relative velocity of

the particle and the fluid. The Drag force is the most dominant force especially for

particle Reynolds numbers less than 100. The drag force uses the Stokes law for

Reynolds numbers less than 1, the modified Stokes law for 1<Reynolds numbers<500

and Newton’s law for 500<Reynolds numbers<2x105. The range of Reynolds numbers

seen in this study were less than 1for all particles.

The mass effect term is the force of accelerating the fluid around the particle. The

mass that is accelerated is equal to one half of the displaced mass of the fluid. The history

term is the force required to maintain the flow pattern. These effects with the Basset force

are negligible for particle densities much greater than the fluid density.

The Saffman’s lift force is caused by the shear of the surrounding fluid which

causes a non-uniform pressure distribution acting on the particle. If the particle is leading

the fluid motion then the lift force pushes the particle down the velocity gradient towards

the wall. If the particle lags the flow then the lift force moves up the velocity gradient

towards the free stream. This force has been shown to affect the deposition velocity and

Wang et. al. has shown that neglecting this force can reduce the deposition rate slightly

[44]. The Saffman Lift force was originally given by Saffman as [34]:

27
 [16]

where Kr is the local velocity gradient. Li and Ahmadi developed a more generalized

version of the Saffman Lift source so that it would apply to three dimensional shear fields

[27]:

[17]

where S is the ratio of particle density to fluid density and K=2.594 is a constant

coefficient. dij is the deformation tensor and is given as:

[18]

The buoyancy force is important for large particles with large changes of position.

Since the particles in this study are small and the domain is small in size, the gravitational

force is expected to be negligible. The buoyancy force is usually negligible for particle

deposition studies and is neglected for this study.

The Magnus force is the lift developed due to rotation of the particle. The lift is

caused by the pressure difference that occurs with the difference in flow on opposite sides

of the sphere. Kallio & Reeks found that this force is an order of magnitude less than the

Saffman lift forces in most cases and could be neglected [23].

The Brownian and Thermophoretic forces are important for sub-micron particles.

The Brownian force is caused by random interactions with agitated gas molecules in the

fluid. The thermophoretic force is caused by the unequal momentum exchange between

28
particle and fluid. Temperature gradients in the flow can cause the momentum exchange

to be higher in the direction of decreasing temperature. Both of these forces are neglected

since the particles used in the study are greater than 0.03 µm.

Another force that can be important is the rarefaction force. The rarefaction force

is the influence of gas particles colliding with the coal ash particle in a non-continuum

regime. The Knudsen number, KN, is the ratio of the mean free path of the gas molecule

and coal ash particle size:

[19]

where λ is the mean free path of gas particles which is defined as:

[20]

[21]

where R is the specific gas constant.

29
 Table 3 – Particle Forces [37]

FORCE Domain of Importance Included

1 Drag Dominant force for particle motion; Rep<100, spherical YES

particles assumed. Stokes' law: Rep < 1; Modified Stokes' 1<Rep < 500

Law: 1<Rep < 500; Newton's Law: 500< Rep< 2 x 105

2 Rarefaction Important for sub-micron particles (<1µm) for non-continuum NO

Effect regime; Continuum: Kn < 0.1; Transition: 0.1 < Kn < 10 0.1<Kn<10

Free-molecule: Kn > 10. Important for Kn>0.02

3 Virtual Mass Important for small values of particle material density to gas NO

Effect density; ρ p / ρ f << 1 ρ p / ρ f >> 1

4 Basset Important for small values of particle material density to gas NO

density; ρ p / ρ f << 1 ρ p / ρ f >> 1

5 Saffman Lift Non-trivial magnitudes only in the viscous sub-layer; slight YES

decrease in deposition rate if neglected. More accurate deposition rate

6 Magnus Lift force due to particle rotation; at least an order of NO

magnitude lower than Saffman force in most regions; No particle rotation in the flow

7 Gravity Important in Stokes Regime for large particles NO

Very small particles

8 Thermophoretic Important for sub-micron particles and Kn < 2 NO

Kn > 10

9 Brownian Important for dp < 0.03 µm NO

dp > 0.03 µm

10 Intercollision Important for high-volume fraction of particles NO

Low volume fraction of particles

30
3.2 Injection Properties

The particles in this study were injected at specified locations across the inlet

boundary as shown in Fig. 15. There are 900 injection points located at 30 radial

positions and 30 angular positions. The injection points are grouped into three groups.

Group A is located near the hub end wall, group C is located near the case end wall and

group B is located in between A and C.

Figure 15 – Particle injection locations

Particles were assumed to be in thermal equilibrium with the flow at the inlet and were

injected at the local flow velocity and direction. Each injection point injects one particle

at a time.

31
3.3 Random walk model

The dispersion of particles can be modeled using the stochastic tracking, also

known as random walk, in FLUENT. FLUENT predicts the trajectory of the particles

using the mean velocity parameter, [18]. For turbulent flows, the instantaneous random

velocity, , is used to model dispersion. The local flow velocity is determined by the

relation

[22]

The random effects of turbulence are based on the varying nature of .

The method for determining is based on the random walk model. The model

assumes small scale eddies exist in the flow. These eddies are random and have

instantaneous flow velocities determined by the following relation

[23]

where ζ is a normally distributed random number and the second term is the RMS

fluctuating component. Assuming isotropy the RMS values are defined as

[24]

The time scale of the eddy is used to determine when the instantaneous velocity is

calculated. The integral time, TI, is found using the relation

32
 [25]

Infinitely small particles with zero drift velocity will have the integral time become the

fluid Lagrangian integral time, TL. TL is approximated using the relation

[26]

where CL is an empirically determined constant and ω is the specific dissipation rate from

the k-ω turbulence model. CL is set to 0.30 for k-ε and k-ω turbulence models. The time

scale of the eddy, TE, is defined as

TE = 2 TL [27]

The time for the particle to cross the eddy is defined as

[28]

where τ is the particle relaxation time and LE is the eddy length scale. If tcross is shorter

than the lifetime of the eddy, then the particle has left the current eddy into a new one and

is recalculated. If TE passes before the particle has crossed the eddy, the eddy decays.

A new eddy is formed and the instantaneous velocity is recalculated.

33
3.3.1 Perturbation velocity corrections

The assumption of the RMS fluctuating components being isotropic is accurate

for free stream flow. However, near the wall, the RMS values become highly anisoptric.

FLUENT was found to over-predict the wall normal perturbation velocities, v’+, by more

than an order of magnitude [13]. Fig. 16 shows the magnitudes of the perturbation

velocities near the wall.

Figure 16 –RMS values of velocity in the boundary layer [13]

During preliminary simulations, particles were prone to being trapped in strong eddies

near the leading edge of the turbine vane. These particles would have incomplete

trajectories and would lead to inaccurate particle impact locations.

Attempts to correct this in FLUENT have been suggested with some success.

Dehbi replaces the calculation for u’, v’, and w’ with functions using Y+ values [13].

These functions are based on channel flow and lead to some improvement in the

34
accuracy. Longmire used corrections based on wall functions curve fitted to discrete

Navier-Stokes results [29].

For this study channel flow is too simplistic to model the complex geometries in a

turbine passage. To reduce the error associated with the invalid wall-normal

perturbations, the instantaneous velocities were set to zero for the regions where Y+

values are less than 200. This means the particles will use the mean flow velocities in

equation 13 to determine the particle trajectory in this region. Particles outside of this

limit will have the fluctuating velocity values calculated based on the isotropic

assumption. The Y+ value of 200 was chosen by Longmire as the limit where the error

would still be negligible [29]. A study was performed using a Y+ value of 100 as the limit

and there were no changes in the results. For this study, the interest is in the general

location of deposition and the deviations to the particle trajectories near the wall are

ignored for Y+ values less than 200. The deviations in the free stream are more important

for determining the directions the particle will travel.

3.4 Particle material properties

The material properties of the particle are important in determining how the

particle responds to changes in the flow. The density of the coal fly ash, ρ, the heat

capacity at constant pressure, CP, and the thermal conduction value, K, are used to

determine the thermal diffusivity which is a measure of how quickly the particle changes

temperature with the flow. The material properties were not known for the fly ash used in

35
the TuRFR experiment so the material properties with a similar coal ash particle from

Brigham Young were used [2].

Table 4 – Particle Material Properties BYU [2]

ρ (kg/m3) CP (J/kg-K) K(W/m-K)

990 984 0.5

3.4.1 Ash composition

The chemical composition for the coal fly ash is dependent on the type of coal

used in the combustion process. Coal is found in many geographical regions and has

slight variations in chemical composition. Table 5 shows the chemical properties from

different coal deposits.

36
 Table 5 – Coal fly ash composition (%w.t.) [41]

Lethabo HN115 KL1 HNP01 WV UF WY ILL Pitt ND SE Ohio

SiO2 60.33 42.3 47.1 50.1 60.43 46.86 35.59 46.62 50.37 23.68 32.9

Al2O3 29.73 34.5 35.3 32.9 30.69 23.24 15.15 14.41 21.04 7.94 20.3

Fe2O3 2.68 6.17 4.72 8.42 3.03 21.02 7.53 26.80 21.23 9.82 40.6

TiO2 1.37 2.24 1.96 1.51 2.02 1.00 1.40 0.73 1.14 0.47 -

P2O5 0.41 0.55 0.19 0.34 0.69 0.69 3.02 1.00 0.90 3.85 -

CaO 3.58 8.55 5.67 1.37 0.42 1.68 18.92 2.74 1.37 18.43 2.93

MgO 1.19 1.00 0.75 0.62 0.67 1.00 4.76 0.72 0.65 7.44 -

Na2O 0.14 0.21 0.26 0.45 0.27 0.27 2.10 0.88 0.53 10.20 -

K2O 0.45 0.76 1.33 2.17 0.97 2.65 1.01 3.15 2.00 1.35 2.48

S 0.10 0.04 0.03 0.01 0.81 1.60 10.53 2.94 0.78 16.82 -

MnO 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -

SO3 - - - - - - - - - - 0.83

Silicon, Aluminum, and Iron are the common elements in coal ash. The TuRFR

experiment uses the ash from south east (SE) Ohio. Notice that the ash is high in iron

content compared to all other ash types.

3.4.2 Size distribution

The size of the coal ash particle is important in determining the types of forces

that act on its trajectory in a flow. Advanced filters are commonly used in modern

turbines to remove particulate above 10 µm in diameter. The particulate size distribution

for the TuRFR facility is shown in Fig. 17.

37
 6.00%
Mean = 14.117 µm
Median = 12.038 µm
5.00% Size (µm) Mass %
0-5 23.16%
5-10 31.15%
4.00% 10-20 25.81%
20-40 14.82%
40-60 3.82%
Mass % 3.00% 60-100 1.23%

2.00%

1.00%

0.00%
0 20 40 60 80 100
Diameter (µm)

Figure 17 – Particle size distribution results from Coulter counter for Bituminous coal fly

ash [38]

Nearly half of the particles in the TuRFR facility would normally be filtered in modern

turbine engines. Since the TuRFR results are being simulated in this study, the study will

include the larger size particles as well. The distribution of particles used in the

simulation is shown in table 6.

Table 6 – Particle size distribution for FLUENT simulation

Particle Diameter (µm) 1 10 20 30 40 50

Mass % 15 30 25 20 15 10

The logarithmic form of the Rosin-Rammler distribution equation was used to model the

size distribution. The Rosin Rammler distribution equation is of the form [18]:

38
 [29]

where dp is the diameter of the particle, dm is the mean diameter and n is the spread

parameter. Yd is the percentage of particles with diameters greater than dp and the

percentage is assumed to be exponential. For the size distribution used in the TuRFR

facility the Rossin-Rammler function was fitted to the experimental data using dm = 12

and n= 2.78. The Rossin-Rammler fit is shown in Fig. 18 with the experimental size

distribution.

35
30
25
Mass %

20
15
10
5
0
0 20 40 60 80 100

Experiment Rosin Rammler

Particle Diameter (μm)

Figure 18 – Size Distribution Modeling

39
3.5 Particle Trajectories

Particle diameter has a major effect on the trajectory the particle will take.

Particles with larger diameters have more inertia due to their increased mass. Smaller

particles have less mass and tend to be more susceptible to forces acting on them from the

flow. Fig. 19-21 show the trajectories for particles with diameters of 1, 20, 50 μm.

Particles enter from the inlet on the left and will exit at the outlet to the right. Particles

that pass through the upper/lower bounds will reappear at the opposite side with the same

trajectory.

Figure 19 – 1 μm diameter particle trajectory

40
 Figure 20 – 20 μm diameter particle trajectory

Figure 21 – 50 μm diameter particle trajectory

1 μm particles follow the flow closely and tend to not contact the surface. Larger

diameter particles will likely strike the surface of the vane and depending on the

conditions will reflect multiple times before exiting the domain. Particles bouncing

between vanes will increase the number of impacts considerably.

41
 CHAPTER 4: STICKING MODELING

4.1 Impact, sticking, capture efficiency definitions

When modeling deposition, it is useful to define three parameters to measure the

performance of the simulations. The impact, sticking and capture efficiencies are defined

below:

[30]

[31]

[32]

The impact efficiency is useful for measuring how likely particles are to hit the surface.

Small particles are more likely to follow the flow and less likely to hit the surface. Larger

particles will have more inertia and will hit the surfaces more. The impact efficiency will

be higher for larger particles than smaller particles if the same number of both particles is

injected. The sticking efficiency shows how likely particles stick to the surface when the

particles contact the surface. The capture efficiency shows how much particulate mass

deposits on the surface. Based on these definitions, impact efficiency is equal to the

sticking efficiency and capture efficiency combined.

42
4.2 Wall-particle interaction

There are several forces that help a particle stick to a surface. The most common

sticking forces are the van der Waals force, electrostatic force and liquid bridging. The

electrostatic force is the attraction of charges that build up on the particle and the a

surface. The van der Waals force is the collection of forces that attract molecules together

that do not include the electrostatic force. The liquid bridge is the force of a semi-molten

particle sticking to a surface. Fig. 22 shows the relation of these forces for various

particle diameters.

Figure 22 – dimensionless sticking forces versus particle diameter [6]

It can be seen that smaller particles will have a stronger sticking force compared to larger

particles. For temperatures below the molten temperature of the particulate, the van der

43
Waals force is the dominant sticking force. The other two forces are neglected for this

study.

The van der Waals force can be expressed in the form:

[33]

where A is the Hamaker constant and z is the distance between the particle and surface.

For particles that contact the surface, the particle deforms and the equation for Fpo

becomes:

[34]

where ka is a constant equal to 3π/4 and WA is the work of sticking that is dependent on

the material [15]. For silicon, WA is 0.039 J and is the value that is used in this study.

4.3 Critical impact velocity model

The critical impact velocity is a model that uses the impact velocity component

normal to the surface to determine if the particle sticks to the surface. The critical impact

velocity is used as criterion for determining whether the particle sticks or reflects and is

affected by various properties of the flow and particle. The critical velocity in this study

is given as [9]:

44
 [35]

[36]

[37]

[38]

where Vcr is the critical velocity, E is the composite Young’s modulus, Es is the surface

Young’s modulus, EP is the Young’s modulus of the particle, νs is the surface poisson

ratio and νP is the particle poisson ratio. Particles that impact the surface with a normal

velocity below the critical velocity will deposit. Those particles that impact with

velocities above the critical impact velocity will reflect off the surface and continue their

trajectory until they leave the domain or deposit on the surface further downstream.

4.3.1 Particle impact properties

The Young’s modulus for the particle is difficult to measure experimentally due

to the high temperatures the particles are in. The Young’s modulus is sensitive to the

temperature of the particle and can vary drastically at high temperatures. The material

properties of the particle can also affect the Young’s modulus. Due to the complications

of measuring the Young’s modulus for the particle, it is approximated using empirically

derived functions. El-Batsh used the Young’s modulus as a value that can be determined

empirically from experimental deposition results. El-Batsh developed a function of the

45
Young’s modulus to match his predicted capture efficiencies to experimental results. The

function that El-Batsh used was [15]:

for Tp>1100 K [39]

for Tp>1300 K [40]

Ai developed a function for the particle Young’s modulus based on coupon studies[2].

The function Ai used was [3]:

[41]

where Tg is the free stream gas temperature above the surface. Ai & Fletcher made

modifications to the function to include the surface temperature effects. The new

correlation is[3]:

[42]

[43]

Another modification was made in this study to the function developed by Ai & Fletcher.

The temperature near the surface varies due to the thermal boundary layer. For flows

where the surface is at a lower temperature than the free stream, the flow right above the

surface will be lower temperature than the free stream. Depending on the thermal

properties of the particle, the particle temperature may be at a higher temperature than the

flow near the surface. For this study, TAVG was modified by replacing the Tg in equation

43 with the particle temperature Tp.

46
4.4 Critical viscosity model

One drawback of the critical velocity model is the difficulty in calibrating the

Young’s modulus. Changes in coal ash properties require separate calibrations. Tafti et.

al. developed a sticking model based on the particle viscosity [41]. The particle viscosity

changes with temperature and the relationship can be predicted based on the coal ash

properties.

Tafti et. al. based their model on the probability of the particle sticking. The ash

softening temperature is the critical sticking temperature, TS. Particles above this

temperature will have a sticking probability of one. Particles with temperature much less

than the critical sticking temperature will have a sticking probability of zero. For particles

with temperatures in-between, the probability function is found using [41]:

[44]

where µcrit is the viscosity of the particle at the critical sticking temperature and µTp is the

viscosity of the particle at the current particle temperature. The variation of the

probability function with temperature is shown in Fig. 23.

47
Figure 23 – Transition of modified deposition model using critical viscosity approach

[41]

4.4.1 Ash composition model

Senior et. al. developed a method to determine the coal ash particle viscosity

based on the coal ash chemical composition. Silicate melts are made up of the building

block SiO44- and can form networks that compose the particle glass structure. There are

three categories of cations and each interacts with the network in different ways. The first

category is known as the glass formers which form the basic anionic polymer unit. Si4+,

Ti4+, and P5+ are glass formers. The second category is the modifiers, which disrupt the

polymer chains bonding oxygen. Ca2+, Mg2+, Fe2+, K+ and Na+ are modifiers and are

effective at terminating chains in the structure. The third category is the amphoterics,

which act either as formers or modifiers. Al3+, Fe3+, B3+ are amphoterics. Modifier ions

disrupt the glass structure and tend to lower the viscosity. Amphoterics interact with the

modifiers and can form bonds to balance their charge. If there is an abundance of

amphorterics, they can act as modifier ions as well. The parameter that relates the balance

48
of amphoterics and modifiers is given by the ratio of non-bridging oxygens to tetrahedral

oxygens (NBO/T) given by Vargas[43]:

[45]

4.4.2 Viscosity versus temperature relationship

N’Dala et. al. have shown that the temperature dependence of viscosity of silicate

and aluminosilicate melts can be described by[30]:

[46]

where A and B are constants which depend on the chemical composition. Senior

conducted experimental tests to develop a curve fit using A, B and NBO/T. The model

developed by Senior is split into two temperature regions with AL being A for low

temperatures and AH being A for high temperatures. The functions for A and B developed

by Senior et. al. are [36]:

[47]

NBO/T≥1.3 [48]

0.2≤NBO/T<1.3 [49]

0.0≤NBO/T<0.2 [50]

NBO/T<0.0 [51]

[52]
49
where α is the ratio of mole fraction network modifiers to the sum of network modifiers

and amphoterics, and N is the mole fraction of SiO2 in the glass. The coefficient bi is

shown in table 7.

Table 7 – Fitted Values of Coefficients for Compositional Dependence of B[36]

High Temperature Low Temperature

b0 -224.98 -7563.46

b1 636.67 24431.69

b2 -418.70 -17685.4

b3 823.89 32644.26

b4 -2398.32 -103681.0

b5 1650.56 74541.33

b6 -957.94 -46484.8

b7 3366.61 146008.4

b8 -2551.71 -104306.0

b9 387.32 21904.63

b10 -1722.24 -68194.8

b11 1432.08 48429.31

The coefficients A and B are used in equation 46 to determine the viscosity at the particle

temperature. There will be a viscosity based on the low temperature coefficients and the

high temperature coefficients. The higher of these two viscosities is the value used.

50
 4.4.3 Relationship limits

One limitation of the method developed by Senior et. al. is that the model does

not work with ash containing high amounts of iron or no modifier ions. Nowok looked

into finding a correction that would describe the chemical nature of iron in the ash

composition and how to predict the viscosity from the iron content [31]. The limitation of

this correction is knowing the ratio between the Fe3+ and Fe2+ ions at temperature. Due to

the complexities associated with the high iron content, the viscosity-temperature

relationship could not be determined for the ash used in the TuRFR experiment. Instead,

the ash composition from the BYU facility was used due to its lower iron content and the

fact that the sticking model had been calibrated to use the ash [3].

Table 8 – Coal Ash Chemical Composition used at BYU[3]

Mass Element Weight %

Bulk Apparent
Mean
Density Density Na Mg Al Si P S K Ca Ti V Fe Ni
Diameter

13.4 µm 0.99 g/cc 1.98 g/cc 6.9 3.6 18 47 1.6 1.8 2.6 8.7 1.6 0 6.4 0

For this ash composition, the constants in equation 46 were determined using the method

described by Senior et. al. The viscosity-temperature relationship for the coal fly ash at

BYU is shown in Fig. 24 with fly ash used by Tafti et. al.

51
 Table 9 – Viscosity-Temperature relationship constants

A B

High Temperature -11.3059 15.9561

Low Temperature -44.2546 49.8791

1.E+09
1.E+08 BYU Ash
1.E+07 Lethabo
Viscosity (Pa-s)

1.E+06 Hn115
1.E+05 KL1
1.E+04 HNPO1
1.E+03
WV
1.E+02
UF
1.E+01
WY
1.E+00
ILL
1250 1300 1350 1400 1450 1500 1550 1600
Pitt
Temperature (K)

Figure 24 – Viscosity – Temperature relationship for various ash types [41]

If the ash composition is unavailable or if there is sufficient experimental data, the

viscosity model can be used by calibrating the constants A and B, and the critical sticking

temperature. Tafti et. al. used this method for studying PVC particles in a heated flow

[41]. The sticking probability used for ash can not be applied to PVC but PVC is assumed

to have the same exponential trend for the viscosity and temperature relationship.

[53]

52
The sticking probability at zero Kelvin, PS(T0) is set to have PS equal to 1.00E-6 and at

86% of the melting temperature of the PVC. k is set to have Ps(TMelt) equal to 1.0 at the

melting temperature. This calibration can be done with any particle material as long as

the melting temperature of the material is known.

4.4.4 Critical viscosity model sticking criteria

With the viscosity-temperature relationship the probability of sticking can be

determined by specifying the critical sticking temperature. The softening temperature for

most ash particles was found to be around 1478 K so 1500 K was used as the critical

sticking temperature in this study[42]. Particles with a temperature above this critical

temperature were assumed to always stick.

To determine if a particle sticks for temperatures below the critical sticking

temperature a random number generator (RNG) was used to provide a number between 0

and 1. There are several types of RNG’s used in literature. The randomness of these

RNG’s is dependent on a seed variable being different each time the function is used. A

modified sawtooth wave was used with a flow variable as a seed. The form of the

modified sawtooth wave is given as:

[54]

53
where X is the random number, t is the seed variable and a is the period of the function. a

can be adjusted so that small changes in the seed variable will lead to a nearly random set

of numbers between 0 and 1.

1
0.9
0.8
0.7
0.6
X(t)

0.5
0.4
0.3
0.2
0.1
0
0 0.5 1 1.5 2 2.5 3 3.5 4

Sawtooth wave Random number

Seed

Figure 25 – Sawtooth wave and random number generated

Fig. 25 shows the sawtooth wave with a linearly changing seed and a random number

based on a seed derived from the local flow properties.

When a particle impacts a surface, the viscosity of the particle is determined

based on the particle temperature and the probability of sticking is determined based on

the critical sticking viscosity. The RNG is used to provide a number between 0 and 1. If

the random number generated is greater than the probability of sticking, the particle

reflects off the surface and returns into the flow. If the random number is less than the

54
probability of sticking, the particle sticks to the surface. Higher particle temperatures will

have higher probabilities of sticking.

4.5 Pre-Existing-Deposition factor

Experiments have shown (see section 4.8) that deposits are slow to grow initially

but grow rapidly over time. Initial deposits act as collectors for future particles to stick.

The Pre-Existing-Deposition (PED) factor is a factor that increases the likelihood of the

particle sticking in the simulation. For the critical velocity model, the critical velocity

was increased by a factor of 10 on cells with pre-existing deposition and for the critical

viscosity model the viscosity of the particle was increased by a factor of 1.1. These

factors were chosen to increase the likelihood of a particle depositing but have not been

calibrated. The result of the PED factor is an increase of deposition near the leading edge

and a decrease near the trailing edge. Particles that would reflect off the leading edge

without the PED factor were more likely to be captured using the PED factor.

4.6 Detachment Model

Deposits can detach from the surface if the shear forces are sufficient to overcome

the sticking force. Many particles stick temporarily before being dislodged from the

surface. The amount of deposit that is on the blade is determined by the balance of

particles sticking and particles detaching.

55
4.6.1 Individual particle detachment

Soltani & Ahmadi studied different mechanisms for particle detachment and

found that for spherical particles they are detached predominantly by rolling [39]. Sliding

and lifting are less important. Das et. al. also found that rolling was the more dominant

removal mechanism [12]. El-Batsh used this reasoning to develop a model for particle

detachment [16]. El-Batsh goes into deep detail about this mechanism so a brief summary

is described below.

Since the primary removal mechanism is the particle rolling off the surface, the

model for determining if a particle detaches is the sum of moments acting on the

spherical particle.

Figure 26 – Geometric features of a spherical particle stuck to a smooth surface [16]

[55]

56
where α is the relative approach between the particle and the surface, FL is the fluid lift

force acting on the particle, and a is the radius of the contact surface. α is small for most

cases and is much smaller than Dp/2. FL is also small compared to the other forces.

Equation 55 becomes:

[56]

For spherical particles and Reynolds number less than 1 FD can be written as:

[57]

where uc is the fluid velocity at the center of the particle, Acp is the cross sectional area of

the particle, f is the correction factor for the particle being near the wall as provided by

Soltani & Ahmadi, and Cu is the Cunningham correction factor which is defined as [39]:

[58]

Using f=1.7 and Acp = (π/4)Dp2 and Stokes drag for CD, FD becomes:

[59]

[60]

57
The contact radius, a, can be found using the relation developed by Soltani & Ahmadi

and Rimai et. al. [32,39]:

[61]

where,

[62]

Kc is the composite Young’s modulus and combines the material properties of the particle

with those of the surface. Using the sticking force, Fpo, from equation 34, equation 56 can

be rewritten to form the critical wall shear velocity uτc:

[63]

If the wall friction velocity is greater than uτc the particle detaches from the surface.

The critical viscosity model proposed by Tafti et. al. does not include a

detachment mechanism. For this study, the detachment model from El-Batsh was

incorporated with the viscosity model to form a hybrid model. The detachment

mechanism is important for determining deposition growth and should be included for all

modeling.

58
4.6.2 Large scale deposit detachment

Experiments (cf. section 4.8) have shown that some large scale deposits are pulled

off the surface. It is believed these deposits are removed by the additional pressure drag

that exists by their large size (Reynolds number 1). A method was developed to

simulate this removal mechanism by modeling the deposits as right circular cones.

Referring to Fig. 27 the deposit accumulation was modeled as a cone for which h is the

total height of the cone, hC is the height to the center of the cone, and B is the area of the

cone base.

dh
h
FD
hC

Fpo

Figure 27 – Cone and cylinder modeling

The volume of the cone is based on the amount of particulate mass and particulate

density on a surface face. The base area of the cone is set to the surface area of the cell

face. The height of the cone is determined by the volume of the cone and the base area.

As more particulate mass is added to the surface face, the cone grows higher from the

surface. The drag on the cone is modeled by representing the cone as a series of cylinders

with variable radius. The coefficient of drag was determined for each cylinder based on

59
the Reynolds number for the cylinder. The relation between Reynolds number and CD is

shown in Fig. 28.

Figure 28 – Coefficient of drag for an infinite cylinder vs. Reynolds numbers [45]

The drag on each cylinder was then integrated for the cone to determine the pressure

drag. The same methodology for the detachment of a single particle was applied to the

entire cone. Equation 56 would become:

[64]

where Fpo is now,

[65]

If the sum of moments due to the pressure drag exceeds the moment due to the sticking

force, the cone will be removed. If the sticking force prevailed, then the cone would

remain. This model does not take into account deposits merging together from adjacent
60
cells. This model is a first step at predicting the complex removal mechanisms that occur

with large scale deposits and neglects the complex geometries of deposits that form in

real life. The effects of this model are shown below.

Pressure Side

C
H a
u s
b e

W W
a a
l l
l l

(a) (b)

Figure 29 – Large Scale Deposition Removal Comparison, (a) No Cone Removal, (b)

Cone Removal

Tall deposits with small bases are likely to be removed by the addition of pressure drag

acting on the deposit. The results from Fig. 29 show that the large scale deposits would

detach. Deposits that remain are likely to merge together and increase their likelihood of

remaining attached.

4.7 Comparison between sticking models

The major difference between the sticking models is that the critical impact

velocity model accounts for the particle velocity as well as the particle temperature. The

61
critical viscosity model is more simplistic by being dependent entirely on the particle

temperature. For flows with wide ranges in velocity magnitudes the sticking models are

expected to produce different results. The critical impact velocity and critical viscosity

sticking models were compared for the VKI and the E3 turbine vanes.

For the VKI case, particles were injected at the inlet boundary with the local

temperature and velocity. The sticking models were compared for the high temperature

VKI case as shown in table 10. 0.0001 kg of particulate was injected with 5 particle

diameters as shown in table 2. The particles either passed through the domain or

deposited on the vane. The injection results for both sticking models are shown in table

10.

Table 10 - VKI Sticking Model Comparison Results

Critical Velocity Critical Viscosity

Total Injected # 600 600

Aborted # 216 224

Incomplete # 0 0

Escaped # 384 376

Impact Efficiency % 77 70

Sticking Efficiency % 75 35

Capture Efficiency % 58 45

Table 10 shows that there are a large number of impacts. This is likely due to large

particles reflecting multiple times before depositing or escaping. The critical viscosity

62
model shows that more of the injected particles deposited on the surface; however the

capture efficiency is lower. The lower capture efficiency shows that smaller particles

were more likely to be captured. Fig. 30 shows the number of particles, regardless of size,

that were captured along the pressure surface of the vane at midspan. Notice that the

viscosity model predicts more captures near the leading edge while the velocity model

predicts more deposits near the trailing edge.

X vs Deposit #
30
25
20
15
#

10
5
0
0 0.2 0.4 0.6 0.8 1

Viscous model Velocity Model

Axial Distance/Axial Chord

Figure 30- Number of Captured Particles along Axial Chord

The low capture efficiency near the leading edge is due to the particles impacting the

leading edge with high normal velocities. Although the viscosity model shows more

captured particles near the leading edge, the percentage of mass depositing near the

leading edge is similar for both sticking models. Fig. 31 shows the percentage of injected

mass that deposited on the surface.

63
 X vs Deposit Mass %
8
7

Injected Mass %
6
5
4
3
2
1
0
0 0.2 0.4 0.6 0.8 1

Viscous model Velocity Model

Axial Distance/Axial Chord

Figure 31 – Mass % of Captured Particles along Axial Chord

Fig. 31 shows higher amounts of mass that deposited near the trailing edge which is due

to large particles favoring the trailing edge. The differences in the models can be

explained by the sticking criteria. The temperature near the trailing edge drops

significantly. The viscosity model is more sensitive to changes in temperature so a

decrease in captures is expected.

64
 1480

1475

Surface Temperature (K)

1470

1465

1460

1455

1450
0 0.2 0.4 0.6 0.8 1 1.2
Axial Distance/Axial Chord

Figure 32 – VKI Surface Temperature along Pressure Surface

Fig. 21 shows that larger particles tend to bounce between the adjacent vanes and there

are a high number of impacts near the trailing edge. The particles lose some momentum

due to the drag forces resisting cross-streamline flow. Reflecting particles that impact the

trailing edge have a lower normal velocity and are more likely to deposit according to the

velocity sticking model.

For the E3 case, particles were injected at the locations shown in Fig. 15 at the

local temperature and velocity. The particle mass flow was again 0.0001 and was

distributed among 5 particle diameters given by table 2. The E3 high temperature

adiabatic wall case was used for the comparison. The injection results are given in table

11.

65
 Table 11 – E3 Sticking Model Comparison Results

Critical Velocity Critical Viscosity

Total Injected # 22525 20115

Aborted # 9283 7197

Incomplete # 77 52

Escaped # 13165 12866

Impact Efficiency % 63 69

Sticking Efficiency % 92 49

Capture Efficiency % 58 34

The critical velocity and critical viscosity sticking models differed significantly for this

test capture and sticking efficiencies although the impacts were similar. The difference is

likely due to the calibration that was performed by Ai et. al. for the critical velocity

model. Ai et. al. calibrated the model for temperatures ranging from 1300 K to 1450 K.

The critical viscosity model was created using the known properties of coal ash

particulate and no calibration was performed. The simplicity of the viscosity model

requires that the constants A, B and the critical sticking temperature be calibrated using

experimental results. Using the ash composition model is a good initial start for the

values of A and B. The critical sticking temperature can also have a significant effect on

the amount of deposition at a given temperature. The critical sticking temperature of the

coal ash used in the simulation was not known and likely led to the discrepancy between

the models. Fig. 33 shows the deposition locations for the two models.

66
 Pressure Side

C
H a
u s
b e

W W
a a
l l
l l

(a) (b)

Figure 33 – Sticking Model Comparison of Deposition Mass % locations, (a) Critical

Velocity, (b) Critical Viscosity

The deposition locations for the two models show a majority of the deposition being on

the pressure surface near the trailing edge. The viscosity model predicts lower overall

deposition. Due to time constraints, the critical impact velocity sticking model will be

used for studying the end-wall modification in Chapter 5.

4.8 Experimental comparison

The simulations results were compared with the TuRFR facility operating at the

conditions listed in table 1. The amount of deposition on the turbine vanes in the

experiment was very difficult to quantify due to the frailty of the deposition.

Measurements were limited to optical scans of the topography. This work is ongoing and
67
as of this study is unavailable. The results with the simulations will be compared

qualitatively with the experiment to study the deposition locations and relative amounts.

The low temperature case is compared in Fig. 34. The low temperature was found

to be at the lower limit of deposition growth in the facility.

Pressure Side

C
a
H s
u e
b
W
W a
a l
l l
l

(a) (b)

Figure 34 – Low Temperature Comparison between results, (a) Simulation, (b) TuRFR

Experiment

The deposition in the simulation was less than 0.1% of the injected mass. The experiment

shows a light dusting of particulate but no major deposits.

The medium temperature case is compared in Fig. 35. This experimental test

shows a great amount of deposition along the leading edge and pressure side.

68
 Pressure Side

C
a
H
s
u
e
b

W
W
a
a
l
l
l
l

(a) (b)

Figure 35 – Medium Temperature Comparison between results, (a) Simulation, (b)

TuRFR Experiment

The simulation shows increased amounts of deposition compared to the low temperature

case. The simulation also predicts a significant amount of deposition near the trailing

edge which is in contrast to the experimental results. There are multiple explanations for

this discrepancy. One limitation of the current modeling is the deposition growth. The

simulation predicts the particulate trajectories at the time of injection. Transient results

are not modeled. The cone removal mechanism is the first step in exploring the time after

the initial injection of particulate. The E3 geometry in this study does not have film

cooling holes which serve as deposition sites in the experiment.

A time lapse video of the formation of deposition was taken with the TuRFR

facility for a test at 1310 K with no film cooling. The video shows the film cooling holes

69
 acting as the starting locations of deposits and the growth of the deposition. Fig. 36

shows the growth of the deposits on the pressure surface of the vane.

C
FC
a hole H
s u
e b
t=0 t=25 s t=29 s
W W
a a
l l
l l

Lost
Dep. Dep.
Detach near
T.E.

t=35 s t=79 s t=86 s

Figure 36 – Time lapse of deposition growth in TuRFR facility

At t= 25 seconds deposition has formed in the film cooling holes. By t=35 seconds the

deposition has filled in the regions between the film cooling holes. The majority of the

deposition is half way to the trailing edge. At t= 79 seconds, large deposits begin flaking

70
off. By t= 86 seconds large amounts of the initial deposition near the trailing edge is

gone. Fig. 37 shows the vane after the conclusion of the test.

Pressure Side

H C
u a
b s
e
W
a W
l a
l l
l

Figure 37 – Turbine vane after cool down

The deposition in the trailing edge is significantly less than the amount seen during the

test at operating conditions. This is likely due to the thermal expansion of the blades at

the temperature. During cool down the blades contract and the brittle deposits flake off

during this time. The amount of deposition lost during this time is currently unknown;

however, it is clear that the deposition near the trailing edge is most susceptible to being

removed. Inspection of the remaining deposition shows some deposits that are suspended

above the vane leading edge surface as seen in Fig. 38. The fragility of the deposits

requires that the deposits be compared while at temperature.

71
 Figure 38 - Turbine vane after cool down with suspended deposition

4.9 Particle Size Deposition Distribution

An important part of this study is the effect of the size distribution of particles.

Particles follow different trajectories based on their size and mass and the likelihood of

the particle depositing can be affected by its diameter. Filtration advancements have

made the filtration of large particles possible. Understanding the effects of particle size is

important for optimizing filtration techniques to eliminate the particles most susceptible

to depositing. A study was performed using the E3 high temperature adiabatic wall case

to determine the location and amount of deposition for a single particle size. A range of

particle sizes were injected and their deposition locations are shown in Fig. 39. Table 12

shows the capture efficiencies for various particle sizes.

72
 # of Particles

(a) (b) (c)

Figure 39 – Deposition locations for various particle diameters, (a) 1μm, (b) 5μm,

(c) 10μm

Table 12 – Particle Size Comparison Results

1μm 2μm 5μm 8μm 10μm 20μm

Total Injected # 900 900 900 900 900 900

Impact Efficiency % 0.77 2.33 18.9 43.7 68.7 119.9

Sticking Efficiency % 100 100 100 100 100 41.2

Capture Efficiency % 0.77 2.33 18.9 43.7 68.7 49.3

It can be seen from these results that increasing particle diameter above 1μm greatly

increases the amount of particulate that deposits. The particles above 2μm in diameter are

more likely to impact the surface and due to the high surface temperature are more likely

to deposit. There is a cut off above 10 μm where large particles tend to reflect off the

surface more, but there is still a significant amount of deposition. For the 20μm particles

the impact efficiency is above 100% which indicates that some of the particles reflect

73
multiple times before being captured. It can be seen from this study that particles above

2μm should be filtered to reduce the amount of deposition significantly.

74
 CHAPTER 5: END WALL MODIFICATION

5.1 Reasoning for end wall modification

A significant amount of time has gone into the design and manufacturing of

turbine vanes to achieve the power requirements of the turbine. The turbine hub and case

walls are far less complex than the turbine vane and modifying the end walls is seen as

the most ideal form of preventing deposition growth. The E3M was tested under similar

conditions as the TuRFR experiments to determine the benefits of end wall modification.

5.2 Method of adjustment

An observation of the simulation results indicates a distinct pattern of deposition

on the pressure surface of the vane close to the hub end wall. A line of deposits was

observed running parallel with the end wall as shown in Fig. 40. The number of impacts

in this region is also high, which indicates the increased deposition was related to the

flow field.

75
 Pressure Side

C
H
a
u
s
b
e

W
W
a
a
l
l
l
l

(a) (b)

Figure 40 – Impact and Deposit Mass % locations

Streamlines were injected into the flow to visualize the flow in this region. It became

clear from Fig. 41 that a powerful horseshoe vortex existed on the junction of the turbine

vane and hub wall. The horseshoe vortex can pull particles into the wall increasing the

likelihood that particles will stick.

76
 Streamwise
Vorticity

Figure 41 – Flow visualization of horseshoe vortex for E3 (Vortex magnified)

The deposition on the casing side was significantly less. The large inlet angle associated

with the hub wall causes stronger spanwise flow to develop in the pressure surface

boundary layer. This spanwise flow pulls more particles to the surface. The case wall has

a less significant inlet angle and the spanwise flow is far less powerful. This limits the

number of particles being pulled to the surface. The hub end wall was adjusted to have a

less severe inlet angle in the E3M case. The case end wall was left the same. The vortex

for this new configuration is shown in Fig. 42.

77
 Streamwise
Vorticity

Figure 42 – Flow Visualization of horseshoe vortex for E3M (Vortex Magnified)

5.3 Exit and Inlet profile changes

The low Mach flow exiting the combustor means that modification of the end

walls will affect the flow upstream and downstream of the modification. This will affect

the design of the combustor as well as the later stages of the turbine. Therefore,

modification of the end wall should be limited to reduce excessive changes to the

neighboring engine subsystems. The inlet Mach number is shown in Fig. 43. The changes

are limited to the region near the hub end wall. The majority of the fluid remains

unaffected.

78
 Pressure Side

H C
u a
b s
e
W
a W
l a
l l
l

(a) (b)

Figure 43 – Inlet Mach, (a) E3, (b) E3M

Downstream of vane, the Mach number, total pressure and exit flow angle are important

for the performance of later turbine stages. The exit total pressure was shown in Fig. 13

and showed little difference between the E3 and modified E3. The exit Mach number is

shown in Fig. 44.

Pressure Side

C
H a
u s
b e

W W
a a
l l
l l

(a) (b)

Figure 44 – Exit Mach, (a) E3, (b) E3M

The exit Mach values are similar indicating that the modified hub has little effect on the

Mach number downstream. Fig. 45 shows the flow angle relative to the axial axis. The
79
angles are very similar across the span. The modified hub end wall has little effect on the

flow properties. This satisfies the requirement that the end wall modification be minute so

that it can be applied to engines without extensive redesign of the engine components.

Pressure Side

C
H a
u s
b e

W W
a a
l l
l l

(a) (b)

Figure 45 – Exit flow angle, (a) E3, (b) E3M

5.4 Deposition augmentation

The deposition result for the E3M is shown in table 13 and Figs. 46 and 47. The

capture efficiency for the E3M was less than the standard E3. The impacts were

significantly less as well. The impacts in the region near the hub were much less,

although there was a slight increase in the mid-span of the vane.

80
 Table 13 – E3 Sticking Model Comparison Results

E3 E3M

Total Injected # 22525 18105

Aborted # 9283 8826

Incomplete # 77 260

Escaped # 13165 9019

Impact Efficiency % 63 54

Sticking Efficiency % 92 87

Capture Efficiency % 58 47

Pressure Side

C
H a
u s
b e

W W
a a
l l
l l

(a) (b)

Figure 46 – Impact Mass % Locations, (a) E3, (b) E3M

81
 Pressure Side

C
H
a
u
s
b
e

W
W
a
a
l
l
l
l

(a) (b)

Figure 47 – Deposition Mass % Locations, (a) E3, (b) E3M

The deposition decreased in the area of the vortex when the hub wall was modified. The

decrease in capture efficiency also shows that there was an improvement in the

deposition amount. There is a slight increase in deposition in the mid-span region.

Although there is an over-all decrease in deposition the new location of the deposition

could be the source of future problems.

5.5 Recommendations for improvements

The modification of the hub end wall was very simple and appeared to mitigate

the formation of deposits. The modification was the first step in exploring methods of

reducing deposition and there is plenty of room for optimization. Further decrease of the

inlet angle could yield additional results, although there are limitations to how much the

82
angle can be changed before major changes would need to be made to the combustor

outlet. Modification of the casing could also be explored but it is unclear if the benefit

would be as great as the hub wall modification.

83
 CHAPTER 6: CONCLUSIONS AND FUTURE WORK

6.1 Conclusion

A numerical study of the physical mechanisms associated with the deposition of

coal ash particulate on a blade passage of the first stage of the GE- E3 turbine was

performed. The flow was calculated using the software package FLUENT and particulate

was tracked using a Lagrangian approach. Two sticking models were developed and

compared with experimental data. The critical viscosity sticking model was found to

require additional calibration apart from using the viscosity-temperature relationship

derived from the ash composition. The viscosity model was simple to implement and

with proper calibration could be quickly applied to different particulate-laden flows. For

this study, the basic calibration described by Tafti et. al. was insufficient; however, the

experimentally specific calibration described by Tafti et. al. for unknown or complicated

viscosity-temperature relationships may yield better results. The critical impact velocity

sticking model was found to match the initial deposition observed from experimental

tests in the TuRFR facility. The initial deposition locations and quantities were similar

upon injection of particulate but diverged from the experiment as the test advanced.

Modification of the fluid flow due to deposition growth and experimental operation

limitations are likely the source of these errors. Finally, the hub end wall of the turbine

passage was modified to observe the effects of the end walls on deposition. Adjustment

of the hub end wall showed an 18% decrease in over-all deposition on the turbine vane

84
and reduced a significant amount of deposition near the hub and vane interface. The end

wall adjustment also caused a slight increase in deposition in the mid-span region of the

turbine vane. These findings support end wall modification as a viable means of reducing

deposition growth. This paper has shown that:

 The 3D flow features in turbine passages can greatly affect the particle trajectories

and deposition amount

 The balance of deposition removal and accumulation mechanisms is important for

determining how much and how quickly deposition forms on a surface

 The size of particles in the flow can significantly affect the amount of deposition,

with particles less than 2μm in diameter being less likely to deposit while

particles above are highly likely to deposit

6.2 Future Work

Future studies will further develop the transient effects of deposition modeling.

One feature recently made available in FLUENT is the real-time grid adaptation. This

feature can be used to adjust the geometry based on user-defined criteria. The fluid flow

is automatically updated as the geometry varies due to deposition growth and the

performance of the vane can be evaluated after the injection of particulate. The large

scale deposition removal mechanisms will also be developed further to more accurately

predict the locations of deposition formation. Additionally, film cooling will be modeled

on the turbine vane to study the effects of film cooling on deposition. The film cooling

benefit of reducing surface temperatures is offset by the tendency of film cooling holes to

85
serve as sites for deposition collection. The blockage of film cooling holes is also a major

issue with deposition growth. Finally, further modifications to the end walls will be

explored. This study was the first step in exploring the effects of end wall modification

and allow for optimization.

86
 REFERENCES

1) Abuaf, N., Bunker, R.S., Lee, C.P., “Effects of Surface Roughness on Heat Transfer

and Aerodynamic Performance of Turbine Airfoils,” J. Turbomachinery. Vol. 120,

July 1998, pp 522-529.

2) Ai, W., Fletcher, T.H., “Computational Analysis of Conjugate Heat Transfer and

Particulate Deposition on a High Pressure Turbine Vane,” presented at ASME Turbo

Expo 2009: Power for Land, Sea, and Air. June, 2009, Florida, United States. Paper

#: GT2009-59573.

3) Ai, Weiguo. "Deposition of Particulate from Coal-Derived Syngas on Gas Turbine

Blades Near Film Cooling Holes." Thesis. Brigham Young University, 2009. Print.

4) Ajersch, P., Zhou, J.-M., Ketler, S., Salcudean, M., and Gartshore, I.S., 1995,

"Multiple Jets in a Crossflow: Detailed Measurements and Numerical Simulations,"

International Gas Turbine and Aeroengine Congress and Exposition, Houston, TX,

USA.

5) Barringer, M. D., K. A. Thole, and M. D. Polanka. "Effects of Combustor Exit

Profiles on High Pressure Turbine Vane Aerodynamics and Heat Transfer." ASME

Turbo Expo 2006: Power for Land, Sea, and Air (2006). Print.

6) Berbner, S. and Loeffler, F. (1994), Influence of High Temperatures on Particle

Adhesion, Powder Technology, Vol. 78, No. 3, pp. 273-280.

87
7) Bons, J. P., S. T. McClain, Z. J. Wang, X. Chi, and T. I. Shih. "A Comparison of

Approximate versus Exact Geometrical Representations of Roughness for CFD

Calculations of Cf and St." Journal of Turbomachinery (2007). Print.

8) Bons, Jeffrey P. "A Review of Surface Roughness Effects in Gas Turbines." Journal

of Turbomachinery 132 (2010). Print.

9) Brach, R. and P. Dunn, "A Mathematical Model of the Impact and Adhesion of

Microspheres," Aerosol Science and Technology, 16(1), 51-64 (1992).

10) Cramer, K., “Design Construction and Preliminary Validation of the Turbine

Reacting Flow Rig,” 2009, The Ohio State University, Master’s Thesis, Aeronautical

and Astronautical Engineering.

11) Crosby, J.M., Lewis, S., Bons, J.P, Ai, W., Fletcher, T.H., “Effects of Particle Size,

Gas Temperature, and Metal Temperature on high Pressure Turbine Deposition in

Land Based Gas Turbines from Various Synfuels,” presented at ASME Turbo Expo

2007: Power for Land, Sea, and Air. May 14-17, 2007, Montreal, Canada. Paper #:

GT2007-25731.

12) Das, S. K.; Sharma, M. K. and Schechter, R. S. (1995), Adhesion and Hydrodynamic

Removal of Colloidal Particles from Surfaces, Particulate Science and Technology,

Vol.13, No. 3-4, pp. 227-247.

13) Dehbi. A, A CFD model for particle dispersion in turbulent boundary layer flows,

Nucl. Eng. Des. 238 (2008), pp. 707–715

88
14) Dunn, M. G., A. J. Baran, and J. Miatech. "Operation of Gas Turbine Engines in

Volcanic Ash Clouds." Journal of Engineering for Gas Turbines and Power (1996).

Print.

15) El-Batsh, H. and H. Haselbacher, "Numerical Investigation of the Effect of Ash

Particle Deposition on the Flow Field through Turbine Cascades," IGTI, Amsterdam,

Netherlands, GT-2002-30600 (2002).

16) El-Batsh, Hersham. "Modeling Particle Deposition on Compressor and Turbine Blade

Surfaces." Thesis. Vienna University of Technology, 2001. Print.

17) Elghobashi, S., “Particle-laden turbulent flows: direct simulation and closure

models”, Applied Scientific Research, 48: 301-314, 1991.

18) FLUENT 6.3.26 User Guide, 2006, FLUENT Inc

19) Gouesbet, G., A. Berlemont and P. Desjonquères, "Particle Lagrangian simulation in

turbulent flows," Int. J. Multiphase Flow 16, 19 (1990).

20) Hossain, A., and J. Naser. "CFD Investigation of Particle Deposition around Bend Ina

Turbulent Flow." 15th Australasian Fluid Mechanics Conference (2004). Print.

21) Jensen, J.W., Squire, S.W., Bons, J.P., Fletcher, T.H., “Simulated Land-Based

Turbine Deposits Generated in an Accelerated Deposition Facility,” J.

Turbomachinery. Vol. 127, July 2005, pp 462-470.

22) Kaftori, D., Hetsroni, G., and Banerjee, S., “Particle behavior in the turbulent

boundary layer. I. Motion, deposition, and entrainment”, Phys. Fluids 7, 1095-1106,

1995

89
23) Kallio GA, Reeks MW. 1989. A numerical simulation of particle deposition in

turbulent boundary layers. Int. J. Multiph. Flow 15:433–46

24) Kim, J., Dunn, M.G., and Baran, A.J. et al, “Deposition of Volcanic Materials in the

Hot Sections of Two Gas Turbine Engines,” J. Engr. Gas Turbines & Power. Vol 115,

Jul 1993, pp 641-651.

25) Kulick, J.D., Fessler, J.R., and Eaton, J.K., “Particle response and turbulence

modification in fully developed channel flow”, J. Fluid Mech., 277, 109-134, 1994

26) Lawson, S.A., Thole, K., “The Effects of Simulated Particle Deposition on Film

Cooling,” presented at ASME Turbo Expo: Power for Land, Sea, and Air. June 8-12,

2009, Orlando, Florida USA. Paper #: GT2009-59109.

27) Li, A. and Ahmadi, G., “Dispersion and Deposition of Spherical Particles from Point

Sources in a Turbulent Channel Flow”, Aerosol Science and Technology, 16:209-226,

1992.

28) Lewis, S., Barker, B., Bons, J.P., Ai, W., and Fletcher, T.H. "Film Cooling

Effectiveness and Heat Transfer Near Deposit-Laden Film Holes." GT2009-59567.

Presented at the IGTI 2009 conference and recommended for publication in the

Journal of Turbomachinery publication, 2009

29) Longmire, Pamela. "Computational Fluid Dynamics (CFD) Simulations of Aerosol in

a U-Shaped Steam Generator Tube." Thesis. Texas A&M University, 2007. Print.

30) N’Dala, I.; Cambier, F.; Anseau, M. R.; Urbain, G. Viscosity of Liquid Feldspars.

Part I: Viscosity Measurements.Br. Ceram. Trans. J. 1984,83, 108-112.

90
31) Nowok, Jan W. "Viscosity and Structural State of Iron in Coal Ash Slags under

Gasification Conditions." Energy & Fuels (1994). Print.

32) Rimai, D. S.; Demejo, L. P. and Bowen, R. C. (1994), Mechanics of Particle

Adhesion, J. Adhesion Science Technology, Vol. 8, No. 11, pp. 1333-1355.

33) Rudinger, G., “Fundamentals of gas-particle flow”, Elsevier Scientific Pub. Co.,

Amsterdam, New York, 1980

34) Saffman, P. G. (1965), The Lift on a Small Sphere in a Slow Shear Flow, Journal of

Fluid Mechanics, Vol. 22, Part 2, pp. 385-400.

35) Schlichting, H. and Gersten, K. (2000), Boundary Layer Theory, 8th Edition,

Springer.

36) Senior, C. L., and Srinivasachar, S., "Viscosity of Ash Particles in Combustion

Systems for Prediction of Particle Sticking." Energy and Fuels (1995). Print.

37) Shankara, Prashanth S. "CFD Simulation and Analysis of Molten Particulate

Deposition on Gas Turbine Vanes." Thesis. The Ohio State University, 2010. Print.

38) Smith, C., Barker, B., Clum, C., Bons, J.P. “Deposition in a Turbine Cascade with

Combusting Flow.”Presented at the ASME Turbo Expo 2010: Power for Land, Sea,

and Air (2010) Soltani, M., and Ahmadi, G., “On particle adhesion and removal

mechanisms in turbulent flows”, J. Adhesion Sci. Technology, Vol. 8, No. 7, pp-763-

785 (1994)

39) Sundaram, N., Thole, K., “Effects of Surface Deposition, Hole Blockage, and

Thermal Barrier Coating Spallation on Vane Endwall Film Cooling,” J.

Turbomachinery. Vol 129, July 2007, pp 599-607.

91
40) Tafti, Danesh K., and Sreedharan, S.S., "Composition Dependent Model for the

Prediction of Syngas Ash Deposition with the Application to a Leading Edge Turbine

Vane." ASME Turbo Expo 2010: Power for Land, Sea, and Air (2010). Print.

41) Thermal properties of coal ash, MONA LEDESMA AND L.L. ISAACS,1990 -

Materials Research Society

42) Vargas, Signe. "Straw and Coal Ash Rheology." Thesis. Technical University of

Denmark, 2001. Print.

43) Wang, J., Shirazi, S.A., “A CFD Based Correlation for Erosion Factor for Long-

Radius Elbows and Bends”, Journal of Energy Resources Technology 125, 26 (2003)

44) Weisstein, Eric. "Cylinder Drag -- from Eric Weisstein's World of Physics."

ScienceWorld. 2007. Web. 06 Aug. 2010.

<http://scienceworld.wolfram.com/physics/CylinderDrag.html>.

45) Wenglarz, R.A., Wright, I.G., “Alternative Fuels for Land-Based Turbines,”

published in proceedings of the “Workshop on Materials and Practices to Improve

Resistance to Fuel Derived Environmental Damage in Land-and Sea-Based

Turbines,” October 22-24, 2002, Colorado School of Mines, Golden Colorado.

46) Wilcox, D.C. Turbulence Modeling for CFD, La Cañada, CA DCW Industries (1993)

47) Zhang Z, Chen Q. Comparison of the Eulerian and Lagrangian methods for predicting

particle transport in enclosed spaces. Atmospheric Environment 2007; 41(25):5236–

48.

92
 APPENDIX

A.1 Critical Impact Velocity Sticking Model UDF

#include "udf.h"
#include "dpm.h"
#include "mem.h"
#include "sg.h"
#define nu_s 0.27
#define nu_p 0.27
#define YMIN 0.0
#define YMAX 0.4
#define UMEAN 1.0
#define B 1./7.
#define DELOVRH 0.5
#define VISC 1.7894e-05
/*#define RGAS (UNIVERSAL_GAS_CONSTANT/MW)*/
#define Tdatum 288.15
#define NUM_UDM 6

real ParticleTotalMass;
real P_Mass[6];
real P_Impact_Mass[6];
real P_Stick_Mass[6];

DEFINE_DPM_BC(udfdeposition8,p,t,f,f_normal,dim)
{

#if !RP_HOST
real crit_vel,alpha,Tavg,Ep,vcr=0.,Es,k1,k2,calc,E,utau1,dudz,wall_shear_stress;
real Cu,cbar,lamda,kn,val,kc,ucws,utc,h,del,ufr,ff,utau,area,wall_shear_force,utaunew;
real vn=0.;
real utau2, utc1,vpabs=0.,MassImpact;
real nor_coeff=0.8;
real tan_coeff=0.4;
real R=287.;
real tem_Mass=0.;
real tem_Particle_Dia=0.;
real A[ND_ND],ds,es[ND_ND],A_by_es,dr0[ND_ND],ivu=0.,jvv=0.,kvw=0.,du=0.;
real tauwall1=0.,tauwall2=0.,tauwall3=0.,wallfricv1=0.,wallfricv2=0.,wallfricv3=0.;
real yplus,wallfricv4=0.,wallfricv5=0.,uvel,vvel,wvel,dumag=0.;
FILE *fp;
Thread *tcell=P_CELL_THREAD(p);
cell_t c=P_CELL(p);
/*Message("boo1\n");*/
Domain *d;

real normal[3];
93
int i,idim=dim;
real Wa = 0.039;
d=Get_Domain(1);
real x[ND_ND];
for (i=0.; i<idim; i++)
normal[i] = f_normal[i];
/*Message("boo2\n");*/

C_UDMI(c,tcell,11) += 1.0;
/*Message("boo3\n");*/

if(p->type==DPM_TYPE_INERT)
{
/*||||||||||||||||||||||||||||||||computing normal velocity|||||||||||||||||||||||||||||*/
for(i=0.;i<idim;i++)
{
vn += p->state.V[i]*normal[i];
vpabs += pow(p->state.V[i],2.);
}
vpabs = pow(vpabs,.5);
/*Message("boo4\n");*/
/*|||||||||||||||||||||||||||||computing critical velocity||||||||||||||||||||||||||||||*/
Tavg = (F_T(f,t)+P_T(p))/2.;
/*Avg temp is Tavg)*/
/*particle temp is P_T(p))*/
/*face temp from tcell is F_T(f,tcell));*/
/*normal face temp is F_T(f,t));*/

/*NEW CORRELATION*/
Ep=(3.*(pow(10.,20.)))*exp(-0.02365*Tavg);
/*Young's Modulous is for particle is Ep*/
Es=Ep;
k1 = (1.-(nu_s*nu_s))/(3.14*Es);
k2 = (1.-(nu_p*nu_p))/(3.14*Ep);
/*particle density is P_RHO(p))*/
/*Ai particle density is p->init_state.rho);*/
calc = (5.*3.14*3.14*(k1+k2))/(4.*(pow(P_RHO(p),1.5)));
E = 0.51*(pow(calc,(2./5.)));
/*young's modulus is E*/
vcr = pow(((2.*E)/P_DIAM(p)),(10./7.));
/*Message("boo5\n");*/
/*capture velocity is vcr*/
/*particle diameter is P_DIAM(p))*/
/*Ai particle dia is p->init.state.diam*/
/*|||||||||||||||||||||||||||Calculating mass of particles impacting||||||||||||||||||||*/
/*particle mass is P_MASS(p)*/

MassImpact=P_FLOW_RATE(p)*pow(10,9);
C_UDMI(c,tcell,12)+=MassImpact;
/*Message("boo6\n");*/

/*Work of sticking - assumed Wa of silicon */

94
 cbar = sqrt(((8.*C_RGAS(c,tcell)*C_T(c,tcell))/3.14));
lamda = (2.*C_MU_L(c,tcell))/(C_R(c,tcell)*cbar);
/*Ai gas viscosity is VISC, C_MU_EFF(c,tcell)*/
kn = (2.*lamda)/P_DIAM(p);
Cu = 1.+kn*(1.2+(0.41*exp(-0.88/kn))); /*Cunningham Correction Factor*/

Tavg = (F_T(f,t)+C_T(c,tcell))/2.; /*Avg of particle & surface temperature*/

/*Young's Mod of part & surf equal*/

/*Calculating composite Young's Modulus*/

val = ((1.-(pow(nu_s,2.)))/Es)+((1.-(pow(nu_p,2.)))/Ep);
kc = (4./3.)/val;

/*Calculating critical wall shear velocity*/

ucws = ((Cu*Wa)/(C_R(c,tcell)*P_DIAM(p)))*(pow((Wa/(P_DIAM(p)*kc)),(1./3.)));
utc = sqrt(ucws);
/*wall shear vel is utc*/
utc1=pow(ucws,0.5);
/*Message("boo7\n");*/
/*critical wall shear vel is ,utc,utc1); */
/*Calculating wall friction velocity*/

/*new wall friction calculation*/

/*Message("Got to here...\n");*/
Thread *t0=THREAD_T0(t);
cell_t c0=F_C0(f,t);
BOUNDARY_FACE_GEOMETRY(f,t,A,ds,es,A_by_es,dr0);
/*wall distance is ds*/
tauwall1=C_MU_EFF(c,tcell)*(du/ds);
wallfricv1=sqrt(tauwall1/C_R(c,tcell));
/*wallfrictionvelocity 1 wallfricv1*/
/*Message("Got this far...\n");*/
/*NEW 2*/

/*Message("boo8\n");*/

utaunew = sqrt((2.*C_MU_L(c,tcell)*vpabs)/(C_R(c,tcell)*P_DIAM(p)));
/*Message("boo10\n");*/
/*wallfrictionvelocity of weiguo is utaunew*/

/*fp=fopen("Impact2.txt","a");*/
/*fprintf(fp,"%6.3f %6.2f %6.2f %f %f %f %f\n",
P_DIAM(p),vn,vcr,Ep,MassImpact,ucws,utaunew);*/
/*fclose(fp)*/

tem_Particle_Dia=P_DIAM(p)*pow(10,6);

if(tem_Particle_Dia<1)
P_Impact_Mass[0]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>1 && tem_Particle_Dia<3)
P_Impact_Mass[1]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>3 && tem_Particle_Dia<5)

95
 P_Impact_Mass[2]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>5 && tem_Particle_Dia<7)
P_Impact_Mass[3]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>7 && tem_Particle_Dia<10)
P_Impact_Mass[4]+=P_FLOW_RATE(p)*pow(10,9);
else
P_Impact_Mass[5]+=P_FLOW_RATE(p)*pow(10,9);

/*Message("boo11\n");*/

fp=fopen("Impact_Location.txt","a");
fprintf(fp,"%f %f %f %E\n",P_POS(p)[0],P_POS(p)[1],P_POS(p)[2],P_MASS(p));
fclose(fp);

/*Message("E is:%f\n",E);
Message("Ep is:%f\n",Ep);*/
Message("vn is:%f\n",vn);
Message("vcr is:%f\n",vcr);

if(C_UDMI(c,tcell,15)!=0.0)
{
/*Message("Vcr:%f\n",vcr);*/
vcr=10*vcr;
/*Message("Vcr:%f\n",vcr);*/
}

if(abs(vn)>vcr)
{
C_UDMI(c,tcell,13) += 1.0;
C_UDMI(c,tcell,14) += P_FLOW_RATE(p)*pow(10,9);

for(i=0;i<idim;i++)
p->state.V[i]-=vn*normal[i];

/*Apply tangential coefficient of restitution. */

for (i=0; i<idim; i++)
p->state.V[i]*=tan_coeff;

/*Add reflected normal velocity. */

for (i=0; i<idim; i++)
p->state.V[i]-=nor_coeff*vn*normal[i];

/* Store new velocity in state0 of particle*/

for (i=0; i<idim; i++)
p->state0.V[i]=p->state.V[i];

return PATH_ACTIVE;
}

else
/*particle deposits*/

96
 if(wallfricv1 < utc)
{

/*num of particles deposited or num of hits*/

C_UDMI(c,tcell,15) += 1.0;

/* mass of particles deposited*/

C_UDMI(c,tcell,16) += P_FLOW_RATE(p)*pow(10,9);
tem_Mass=P_FLOW_RATE(p)*pow(10,9);
/*Message("boo13\n");*/

if(tem_Particle_Dia<1.)
P_Stick_Mass[0]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>1 && tem_Particle_Dia<3)
P_Stick_Mass[1]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>3 && tem_Particle_Dia<5)
P_Stick_Mass[2]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>5 && tem_Particle_Dia<7)
P_Stick_Mass[3]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>7 && tem_Particle_Dia<10)
P_Stick_Mass[4]+=P_FLOW_RATE(p)*pow(10,9);
else
P_Stick_Mass[5]+=P_FLOW_RATE(p)*pow(10,9);

/*Message("boo14\n");*/

/*fp=fopen("Stick2.txt","a");*/
/*fprintf(fp,"%f %f %6.2f %f %6.2f %6.2f %f\n",
tem_Particle_Dia,tem_Mass,vn,vcr,Ep,ucws,utaunew);*/
/*fclose(fp);*/

fp=fopen("Deposition_Location.txt","a");
fprintf(fp,"%f %f %f %E\n",P_POS(p)[0],P_POS(p)[1],P_POS(p)[2],P_MASS(p));
fclose(fp);

/* real volume,Re,r,R,h,Cd,D,theta,BA,CA;

volume=C_UDMI(c,tcell,16)/(pow(10,9)*P_RHO(p));
BA=F_AREA(A,f,t);
h=3*volume/BA;
R=sqrt(BA/3.14);
theta=atan(h/R);
D=0;

Message("height:%E\n",h);
Message("Base Radius:%E\n",R);
Message("Angle:%f\n",theta*57.295);

for (i=1; i<=3; i++)

{
r=h/(tan(theta))*(1+(1-i)/3);

97
Re=C_R(c,tcell)*sqrt(pow(C_U(c,tcell),2)+pow(C_V(c,tcell),2)+pow(C_W(c,tcell),2))*2*r/C_MU_EFF(c,
tcell);
if(Re<100)
Cd=24/(pow(Re,0.7));
else
Cd=3/(pow(Re,0.25));

CA=(r+(r-h/(3*tan(theta))))*h/3;
D=D+0.5*Cd*(pow(C_U(c,tcell),2)+pow(C_V(c,tcell),2)+pow(C_W(c,tcell),2))*CA;
}

Message("Drag is:%E\n",D);
Message("Fpo is:%E\n",0.75*3.14*Wa*R*R*3/h);
if(D>=0.75*3.14*Wa*R*R*3/h)
{
C_UDMI(c,tcell,15)=0.0;
C_UDMI(c,tcell,16)=0.0;
Message("!!!!!Deposit Removed!!!!!\n");
}
*/

}
/*Message("boo15, did it stick?");*/
return PATH_ABORT;
#endif
}

/*Setting UDM locations and initializing them to 0*/

DEFINE_ON_DEMAND(reset_UDMsnew)
{
int i=0;
Thread *t;
Domain *d;
cell_t c;
face_t f;

d=Get_Domain(1);
Message("Setting UDMs \n");
thread_loop_c(t,d)
{
begin_c_loop(c,t)
{
for(i=11;i<17;i++)
C_UDMI(c,t,i)=0.0;
}
end_c_loop(c,t)
}

98
}

/*Macro to sum up all the UDM values and display*/

A.2 Critical Viscosity Sticking Model UDF

#include "udf.h"
#include "dpm.h"
#include "mem.h"
#include "sg.h"
#include "stdlib.h"
#include "stdio.h"
#include "math.h"
#define nu_s 0.27
#define nu_p 0.27
#define YMIN 0.0
#define YMAX 0.4
#define UMEAN 1.0
#define B 1./7.
#define DELOVRH 0.5
#define VISC 1.7894e-05
/*#define RGAS (UNIVERSAL_GAS_CONSTANT/MW)*/
#define Tdatum 288.15
#define NUM_UDM 6

real ParticleTotalMass;
real P_Mass[6];
real P_Impact_Mass[6];
real P_Stick_Mass[6];

DEFINE_DPM_BC(udfdeposition8,p,t,f,f_normal,dim)
{

#if !RP_HOST
/*Message("boo\n");*/
real crit_vel,alpha,Tavg,Ep,vcr=0.,Es,k1,k2,calc,E,utau1,dudz,wall_shear_stress;
real Cu,cbar,lamda,kn,val,kc,ucws,utc,h,del,ufr,ff,utau,area,wall_shear_force,utaunew;
real vn=0.;
real utau2, utc1,vpabs=0.,MassImpact;
real nor_coeff=0.8;
real tan_coeff=0.4;
real R=287.;
real tem_Mass=0.;
real tem_Particle_Dia=0.;
real A[ND_ND],ds,es[ND_ND],A_by_es,dr0[ND_ND],ivu=0.,jvv=0.,kvw=0.,du=0.;
real tauwall1=0.,tauwall2=0.,tauwall3=0.,wallfricv1=0.,wallfricv2=0.,wallfricv3=0.;
real yplus,wallfricv4=0.,wallfricv5=0.,uvel,vvel,wvel,dumag=0.;
real RANDOM;
real Ac,Bc,Sticky,StickyCr,Sticky_ratio;
FILE *fp;
Thread *tcell=P_CELL_THREAD(p);
cell_t c=P_CELL(p);
99
/*Message("boo1\n");*/
Domain *d;

real normal[3];
int i,idim=dim;
real Wa = 0.039;
d=Get_Domain(1);
real x[ND_ND];
for (i=0.; i<idim; i++)
normal[i] = f_normal[i];
/*Message("boo2\n");*/

C_UDMI(c,tcell,11) += 1.0;
/*Message("boo3\n");*/

if(p->type==DPM_TYPE_INERT)
{
/*||||||||||||||||||||||||||||||||computing normal velocity|||||||||||||||||||||||||||||*/
for(i=0.;i<idim;i++)
{
vn += p->state.V[i]*normal[i];
vpabs += pow(p->state.V[i],2.);
}
vpabs = pow(vpabs,.5);
/*Message("boo4\n");*/
/*|||||||||||||||||||||||||||||computing critical velocity||||||||||||||||||||||||||||||*/
Tavg = (F_T(f,t)+P_T(p))/2.;
/*Avg temp is Tavg)*/
/*particle temp is P_T(p))*/
/*face temp from tcell is F_T(f,tcell));*/
/*normal face temp is F_T(f,t));*/

/*NEW CORRELATION*/
Ep=(3.*(pow(10.,20.)))*exp(-0.02365*Tavg);
/*Young's Modulous is for particle is Ep*/
Es=Ep;
k1 = (1.-(nu_s*nu_s))/(3.14*Es);
k2 = (1.-(nu_p*nu_p))/(3.14*Ep);
/*particle density is P_RHO(p))*/
/*Ai particle density is p->init_state.rho);*/
calc = (5.*3.14*3.14*(k1+k2))/(4.*(pow(P_RHO(p),1.5)));
E = 0.51*(pow(calc,(2./5.)));
/*young's modulus is E*/
vcr = pow(((2.*E)/P_DIAM(p)),(10./7.));
/*Message("boo5\n");*/
/*capture velocity is vcr*/
/*particle diameter is P_DIAM(p))*/
/*Ai particle dia is p->init.state.diam*/
/*|||||||||||||||||||||||||||Calculating mass of particles impacting||||||||||||||||||||*/
/*particle mass is P_MASS(p)*/

MassImpact=P_FLOW_RATE(p)*pow(10,9);
C_UDMI(c,tcell,12)+=MassImpact;

100
 /*Message("boo6\n");*/

/*Work of sticking - assumed Wa of silicon */

cbar = sqrt(((8.*C_RGAS(c,tcell)*C_T(c,tcell))/3.14));
lamda = (2.*C_MU_L(c,tcell))/(C_R(c,tcell)*cbar);
/*Ai gas viscosity is VISC, C_MU_EFF(c,tcell)*/
kn = (2.*lamda)/P_DIAM(p);
Cu = 1.+kn*(1.2+(0.41*exp(-0.88/kn))); /*Cunningham Correction Factor*/

Tavg = (F_T(f,t)+C_T(c,tcell))/2.; /*Avg of particle & surface temperature*/

/*Young's Mod of part & surf equal*/

/*Calculating composite Young's Modulus*/

val = ((1.-(pow(nu_s,2.)))/Es)+((1.-(pow(nu_p,2.)))/Ep);
kc = (4./3.)/val;

/*Calculating critical wall shear velocity*/

ucws = ((Cu*Wa)/(C_R(c,tcell)*P_DIAM(p)))*(pow((Wa/(P_DIAM(p)*kc)),(1./3.)));
utc = sqrt(ucws);
/*wall shear vel is utc*/
utc1=pow(ucws,0.5);
/*Message("boo7\n");*/
/*critical wall shear vel is ,utc,utc1); */
/*Calculating wall friction velocity*/

/*new wall friction calculation*/

/*Message("Got to here...\n");*/
Thread *t0=THREAD_T0(t);
cell_t c0=F_C0(f,t);
BOUNDARY_FACE_GEOMETRY(f,t,A,ds,es,A_by_es,dr0);
/*wall distance is ds*/
tauwall1=C_MU_EFF(c,tcell)*(du/ds);
wallfricv1=sqrt(tauwall1/C_R(c,tcell));
/*wallfrictionvelocity 1 wallfricv1*/
/*Message("Got this far...\n");*/
/*NEW 2*/

/*Message("boo8\n");*/

utaunew = sqrt((2.*C_MU_L(c,tcell)*vpabs)/(C_R(c,tcell)*P_DIAM(p)));
/*Message("boo10\n");*/
/*wallfrictionvelocity of weiguo is utaunew*/

/*fp=fopen("Impact2.txt","a");*/
/*fprintf(fp,"%6.3f %6.2f %6.2f %f %f %f %f\n",
P_DIAM(p),vn,vcr,Ep,MassImpact,ucws,utaunew);*/
/*fclose(fp)*/

tem_Particle_Dia=P_DIAM(p)*pow(10,6);

if(tem_Particle_Dia<1)

101
 P_Impact_Mass[0]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>1 && tem_Particle_Dia<3)
P_Impact_Mass[1]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>3 && tem_Particle_Dia<5)
P_Impact_Mass[2]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>5 && tem_Particle_Dia<7)
P_Impact_Mass[3]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>7 && tem_Particle_Dia<10)
P_Impact_Mass[4]+=P_FLOW_RATE(p)*pow(10,9);
else
P_Impact_Mass[5]+=P_FLOW_RATE(p)*pow(10,9);

/*Message("boo11\n");*/

fp=fopen("Impact_Location.txt","a");
fprintf(fp,"%f %f %f %E\n",P_POS(p)[0],P_POS(p)[1],P_POS(p)[2],P_MASS(p));
fclose(fp);

/*||||||||||||||||||||||||Random number Generator|||||||||||||||||||||||||||||||||||||||||||||||||||||||*/

/*RANDOM=sqrt(sin(P_T(p)*100)*sin(P_T(p)*100));*/
RANDOM=(2*(pow(sin(C_T(c,tcell)*C_U(c,tcell)),2)-
floor(pow(sin(C_T(c,tcell)*C_U(c,tcell)),2)+0.5))+1)/2;
Message("RANDOM is:%f\n",RANDOM);
Message("p is:%f\n",p);
fp=fopen("Random.txt","a");
fprintf(fp,"%f\n",RANDOM);
fclose(fp);
/*||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||*/

/*|||||||||||||||||||||||||||||||||||Viscosity Model||||||||||||||||||||||||||||||||||||||||||||||||||||*/
/*Ash Compisition Model - Weiguo, Hot>1037, critical viscosity is in Hot range */

/*Critical Temp Coefficients*/

Ac=-11.3059025;
Bc=15.95606;
StickyCr=pow(10,Ac+pow(10,3)*Bc/1500);

if(P_T(p)>1500)
Sticky_ratio=1;
else if(P_T(p)<1037)
{
/*Cold Coefficients*/
Ac=-44.25456;
Bc=49.87914;
Sticky=pow(10,Ac+pow(10,3)*Bc/P_T(p));
Sticky_ratio=(StickyCr*1500)/(Sticky*P_T(p));
}
else
{
Ac=-11.3059025;
Bc=15.95606;

102
 Sticky=pow(10,Ac+pow(10,3)*Bc/P_T(p));
Sticky_ratio=(StickyCr*1500)/(Sticky*P_T(p));
}

if(C_UDMI(c,tcell,15)!=0.0)
{
Sticky_ratio=1.1*Sticky_ratio;
}

if(RANDOM>Sticky_ratio)
{
C_UDMI(c,tcell,13) += 1.0;
C_UDMI(c,tcell,14) += P_FLOW_RATE(p)*pow(10,9);

for(i=0;i<idim;i++)
p->state.V[i]-=vn*normal[i];

/*Apply tangential coefficient of restitution. */

for (i=0; i<idim; i++)
p->state.V[i]*=tan_coeff;

/*Add reflected normal velocity. */

for (i=0; i<idim; i++)
p->state.V[i]-=nor_coeff*vn*normal[i];

/* Store new velocity in state0 of particle*/

for (i=0; i<idim; i++)
p->state0.V[i]=p->state.V[i];

return PATH_ACTIVE;
}

else
/*particle deposits*/
if(wallfricv1 < utc)
{

/*num of particles deposited or num of hits*/

C_UDMI(c,tcell,15) += 1.0;

/* mass of particles deposited*/

C_UDMI(c,tcell,16) += P_FLOW_RATE(p)*pow(10,9);
tem_Mass=P_FLOW_RATE(p)*pow(10,9);
/*Message("boo13\n");*/

if(tem_Particle_Dia<1.)
P_Stick_Mass[0]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>1 && tem_Particle_Dia<3)
P_Stick_Mass[1]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>3 && tem_Particle_Dia<5)
P_Stick_Mass[2]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>5 && tem_Particle_Dia<7)

103
 P_Stick_Mass[3]+=P_FLOW_RATE(p)*pow(10,9);
else if(tem_Particle_Dia>7 && tem_Particle_Dia<10)
P_Stick_Mass[4]+=P_FLOW_RATE(p)*pow(10,9);
else
P_Stick_Mass[5]+=P_FLOW_RATE(p)*pow(10,9);

/*Message("boo14\n");*/

/*fp=fopen("Stick2.txt","a");*/
/*fprintf(fp,"%f %f %6.2f %f %6.2f %6.2f %f\n",
tem_Particle_Dia,tem_Mass,vn,vcr,Ep,ucws,utaunew);*/
/*fclose(fp);*/

fp=fopen("Deposition_Location.txt","a");
fprintf(fp,"%f %f %f %E\n",P_POS(p)[0],P_POS(p)[1],P_POS(p)[2],P_MASS(p));
fclose(fp);
}

}
return PATH_ABORT;
#endif
}

/*Setting UDM locations and initializing them to 0*/

DEFINE_ON_DEMAND(reset_UDMsnew)
{
int i=0;
Thread *t;
Domain *d;
cell_t c;
face_t f;

d=Get_Domain(1);
Message("Setting UDMs \n");
thread_loop_c(t,d)
{
begin_c_loop(c,t)
{
for(i=11;i<17;i++)
C_UDMI(c,t,i)=0.0;
}
end_c_loop(c,t)
}
}

/*Macro to sum up all the UDM values and display*/

104
A.3 Random Walk Correction UDF(excerpt from Longmire, [Longmire])

/* The function do_stochastics shows an implementation of the standard stochastic

* tracking in FLUENT and an additional function for boundary layer treatment.
* Limitations:
* - Uses only spherical drag
* - tested for steady state simulations
* - does not work with accuracy control of particle tracks
* - stochastic tracking is activated for ALL defined injections
*
* The function do_wall_stoch also takes into account the effect of fluid
* boundary layer damping of turbulent fluctuations.
* Eddy lifetime is stored in scalar p->user[1]. The residence time of a particle
* in a given eddy is stored in p->user[0]
*
* Usage:
* - compile and load udfs
* - activate drag law
* - activate time step routines
* - for walls with reflection activate DPM boundary condition
* - activate 11 particle scalars and hook SCALAR_UPDATE function
* - hook DEFINE_INIT function
* - activate 11 UDM locations
* - after calc. of flow field execute the set_udm function for setting
* Y_plus values for particles in wall nearest cell (coarse grid), or
* if the sublayer is resolved choose set_udm_finegrid (may take a long time to execute)
* - finally perform particle tracking. The DPM_DRAG function evaluates the particle Y_plus
* and chooses the appropriate time scale for eddy interaction.
*
*
*
*/
#include "udf.h"
#include "dpm.h"
#include "cxsurf.h" /* for screen output with CX_Message(...) */
#include "surf.h"
#include "sg.h"
#define C_MU_FACTOR(I) (I->time_scale_constant)
#define C_MU 0.09
#define MY_WALL_SHEAR_X(f,t) F_STORAGE_R_N3V(f,t,SV_WALL_SHEAR)[0] /* FORCES!!
divide by face area */
#define MY_WALL_SHEAR_Y(f,t) F_STORAGE_R_N3V(f,t,SV_WALL_SHEAR)[1]
#define MY_WALL_SHEAR_Z(f,t) F_STORAGE_R_N3V(f,t,SV_WALL_SHEAR)[2]
#define THERMO 0
#define INTERPOLATION 1
#define DEBUGGER 0
cxboolean do_stochastics = TRUE;
cxboolean do_wall_stoch = FALSE;
cxboolean random_eddy_lifetime = FALSE;
cxboolean do_thermophoretic = TRUE;
real y_T = 11.8;
int update_interval = 100001;
int update_interval_up = 1000000;
105
int update_interval_old = 5000000;
static real therm_diff_factor = 100.;
static real yplus_thresh = 200.0;
static real dt_min = 1.0e-06;
static real te_min = 1.0e-05;
static int msg_count = 0;
/* global function declarations */
static real get_yplus_particle(cell_t c, Thread *t, real *pp);
static real get_tl_plus(real);
static void get_rms_vel_fluc(real,real,real*,real*,real*,real*,cphase_state_t*,cell_t c0, Thread *t0);
static real half_zero(real,real);
enum udm_names{
Y_PLUS,
TAU_WALL,
N_VEC_X,
N_VEC_Y,
N_VEC_Z,
FACE_X,
FACE_Y,
FACE_Z,
TWALL,
QWALL,
WALL_HIT
};
enum p_names{
E_PLUS,
T_EDDY,
Y_PLUS_P,
LAMBDA1,
LAMBDA2,
LAMBDA3,
T_FORCE_X,
T_FORCE_Y,
T_FORCE_Z,
N_EDDY,
NEW_FORCE
};
/********************/
real get_tl_plus(real yplus)
{
real TL_plus = 1.0;
#if 0
if (yplus>5.0 && yplus<yplus_thresh)
{
TL_plus = 7.122 + 0.5731*yplus - 0.00129*yplus*yplus;
}
else if (yplus <= 5.0)
{
TL_plus = 10.0;
}
#endif
TL_plus = 1./(4.529 + 0.0116*pow(yplus,1.75) + 0.768*sqrt(yplus) );
return TL_plus;

106
}
/********************/
void get_rms_vel_fluc(real yplus, real vel_fric, real *Vprime_vec, real *I, real *J, real *K,cphase_state_t*
c,cell_t c0, Thread *t0)
{
real I_mag;
Vprime_vec[0] = 0.0;
Vprime_vec[1] = 0.0;
Vprime_vec[2] = 0.0;

/*Message("yplus is: %E\n",yplus);*/

/* calc. anisoptropic local fluctuation velocities ...*/

/* and assign the local coordinate system vectors ... */
I[0] = c->V[0];
I[1] = c->V[1];
#if RP_3D
I[2] = c->V[2];
#endif
I_mag = NV_MAG(I);
NV_S(I,/=,I_mag); /* normalize */
J[0] = C_UDMI(c0,t0,N_VEC_X);
J[1] = C_UDMI(c0,t0,N_VEC_Y);
#if RP_3D
J[2] = C_UDMI(c0,t0,N_VEC_Z);
#endif
NV_CROSS(K,I,J); /* K = I x J */
#if DEBUGGER
CX_Message("I0 = %f, I1 = %f, I2 = %f\n",I[0],I[1],I[2]);
CX_Message("J0 = %f, J1 = %f, J2 = %f\n",J[0],J[1],J[2]);
CX_Message("K0 = %f, K1 = %f, K2 = %f\n",K[0],K[1],K[2]);
#endif
/* ... finally transform everything back into cartesian coordinates */
NV_S(I,*=,Vprime_vec[0]); /* x-component stored in I*/
NV_S(J,*=,Vprime_vec[1]); /* y-component stored in J*/
NV_S(K,*=,Vprime_vec[2]); /* z-component stored in K*/
/*re-assign fluctuation velocities*/
Vprime_vec[0] = I[0] + J[0] + K[0];
Vprime_vec[1] = I[1] + J[1] + K[1];
Vprime_vec[2] = I[2] + J[2] + K[2];
#if 0
Vprime_vec[0] = 0.*vel_fric;
Vprime_vec[1] = 0.*vel_fric;
Vprime_vec[2] = 0.*vel_fric;
#endif
}
/* function to determine a particle's dimensionless wall distance */
real get_yplus_particle(cell_t c,Thread *t,real *pp)
{
real Y_plus_p,Y_p;
real p_dist[ND_ND];
real A[ND_ND];
real Xf[ND_ND];

107
 if(C_UDMI(c,t,Y_PLUS) < 0.0)
return -1.0;
A[0] = C_UDMI(c,t,N_VEC_X);
A[1] = C_UDMI(c,t,N_VEC_Y);
#if RP_3D
A[2] = C_UDMI(c,t,N_VEC_Z);
#endif
Xf[0] = C_UDMI(c,t,FACE_X);
Xf[1] = C_UDMI(c,t,FACE_Y);
#if RP_3D
Xf[2] = C_UDMI(c,t,FACE_Z);
#endif
NV_VV(p_dist,=,Xf,-,pp);
Y_p = fabs(NV_DOT(A,p_dist));
Y_plus_p = Y_p*sqrt(C_UDMI(c,t,TAU_WALL)*C_R(c,t))/C_MU_L(c,t);
return Y_plus_p;
}
/* distribute user RP variables */

DEFINE_DPM_DRAG(dpm_stochastic_drag,Re,p)
{
cphase_state_t *c = &(p->cphase); /* pointer to continuous phase properties, velocity etc. (dpm.h)
*/
real Vrel;
real drag_factor;
real k,eps,eplus,TL_plus = 1.0; /* dimensionless time scale */
real Y_plus_p = 200.0; /* dimensionless particle dist. */
cell_t c0 = RP_CELL(&(p->cCell)); /* index of cell in which currently tracked particle is in */
Thread *t0 = RP_THREAD(&(p->cCell)); /* cell thread to which above cell belongs */
real coeff;
cxboolean check_time = FALSE;
int dim = ND_ND;
real pv[3] = {0.,0.,0.}; /* particle velocity vector */
real pp[3] = {0.,0.,0.}; /* particle position vector */
real I[3] = {0.,0.,0.}; /* streamwise direction vector at particle location */
real J[3] = {0.,0.,0.}; /* wall normal vector */
real K[3] = {0.,0.,0.}; /* spanwise direction vector */
int i;
real Vprime_rms;
real Vprime_vec[3]; /* fluct. vel. vector */
real vel_fric; /* u_tauwall , frictional velocity */
real y_dist;
NV_V(pv, =, p->state.V);
NV_V(pp, =, p->state.pos); /* position of particle, x,y,z-coord. */
#if !INTERPOLATION
c->V[0] = C_U(c0,t0);
c->V[1] = C_V(c0,t0);
#if RP_3D
c->V[2] = C_W(c0,t0);
#endif
#endif /* INTERPOLATION */

108
#if !PARALLEL
update_interval_old = (int)RP_Get_Real("update_old");
update_interval_up = (int)RP_Get_Real("update_new");
yplus_thresh = RP_Get_Real("yplus_threshold");
dt_min = RP_Get_Real("delta_t_min");
te_min = RP_Get_Real("t_eddy_min");
#endif
Y_plus_p = get_yplus_particle(c0,t0,pp);
/* CX_Message("Y_plus_p = %f, px = %f, py = %f, pz = %f\n",Y_plus_p,pp[0],pp[1],pp[2]); */
if(Y_plus_p < yplus_thresh && Y_plus_p > -0.1 )
{
do_stochastics = FALSE;
do_wall_stoch = TRUE;
}
else
{
do_stochastics = TRUE;
do_wall_stoch = FALSE;
if(p->user[1]>0.)
{
p->eddy_time=DPM_SMALL;
p->user[T_EDDY] = DPM_SMALL;
}
}
y_dist = Y_plus_p/sqrt(C_UDMI(c0,t0,TAU_WALL)*C_R(c0,t0))*C_MU_L(c0,t0);
if(y_dist < 0.5*P_DIAM(p) && Y_plus_p > -0.1)
{
do_stochastics = FALSE;
do_wall_stoch = FALSE;
Trap_Particle(p);
}
#if RP_2D
if (rp_axi_swirl)
{
dim = 3;
pv[2] = p->state.V[2];
}
else
{
pv[2]=0;
pp[2]=0.;
}
#endif
/* toggle for comparison - remove later*/
/* do_stochastics = TRUE; */
/* do_wall_stoch = FALSE; */
if (do_wall_stoch)
{
if (p->eddy_time <= DPM_SMALL)
{
check_time = TRUE;
p->user[N_EDDY] ++;
p->user[E_PLUS]=0.; /* reset interaction time */

109
 vel_fric = sqrt(C_UDMI(c0,t0,TAU_WALL)/C_R(c0,t0));
get_rms_vel_fluc(Y_plus_p,vel_fric,Vprime_vec,I,J,K,c,c0,t0);
#if DEBUGGER
CX_Message("I0 = %f, I1 = %f, I2 = %f\n",I[0],I[1],I[2]);
CX_Message("J0 = %f, J1 = %f, J2 = %f\n",J[0],J[1],J[2]);
CX_Message("K0 = %f, K1 = %f, K2 = %f\n",K[0],K[1],K[2]);
CX_Message("U' = %f, V' = %f, W' =
%f\n",Vprime_vec[0],Vprime_vec[1],Vprime_vec[2]);
#endif
/*TL_plus = get_tl_plus(Y_plus_p);*/ /* calc. dimensionless time scale based on
particle y_plus */
eplus = get_tl_plus(Y_plus_p);
eps = eplus*vel_fric*vel_fric*vel_fric*vel_fric*C_R(c0,t0)/C_MU_L(c0,t0);
k=
0.5*(Vprime_vec[0]*Vprime_vec[0]+Vprime_vec[1]*Vprime_vec[1]+Vprime_vec[2]*Vprime_vec[2]);
/* calc. modified time scale */
p->eddy_time = TL_plus*C_MU_L(c0,t0)/C_UDMI(c0,t0,TAU_WALL); /* eq.
5 of proposal*/
p->eddy_time = 0.15*k/eps;
if (random_eddy_lifetime)
p->eddy_time *= -log(dpm_uniform_random(p)); /* Eq. 23.2-32 User's
Guide, characteristic eddy
lifetime*/
else
p->eddy_time *=2; /* Eq. 23.2-31 User's Guide, characteristic eddy
lifetime */
if (p->state.time == 0.0) p->eddy_time *= dpm_uniform_random(p);
if (p->eddy_time < te_min)
p->eddy_time = te_min;
/*store eddy time*/
p->user[E_PLUS] = eplus; /* time scale constant */
p->user[T_EDDY] = p->eddy_time; /* store eddy lifetime */
p->user[Y_PLUS_P] = Y_plus_p;
p->user[NEW_FORCE] = 1.;
/* velocity fluctuations are calculated using a Gaussian deviate random number
*/
for (i=0; i<3; i++) p->user[3+i] = dpm_gauss_random(p); /* random
fluctuations constant for t_eddy */
/*for (i=0; i<3; i++) p->user[LAMBDA1+i] = 1.0; /* random
fluctuations constant for t_eddy */
for (i=0; i<3; i++) p->V_prime[i] = p->user[LAMBDA1+i] *
Vprime_vec[i]; /******************/
if (msg_count%update_interval_up==0)
{
CX_Message("NEW: particle Y+ = %f , t_eddy = %g
, eplus = %f , dt = %g\n",Y_plus_p,p->user[T_EDDY],p->user[E_PLUS],P_DT(p));
CX_Message("NEW: u' = %f , v' = %f , w' = %f\n"
,p->V_prime[0],p->V_prime[1],p->V_prime[2]);
msg_count++;
}
}
else
{

110
 check_time = FALSE;
/* p->eddy_time = p->user[1] - p->user[0]; */
vel_fric = sqrt(C_UDMI(c0,t0,TAU_WALL)/C_R(c0,t0));
get_rms_vel_fluc(Y_plus_p,vel_fric,Vprime_vec,I,J,K,c,c0,t0);
for (i=0; i<3; i++) p->V_prime[i] = p->user[LAMBDA1+i] * Vprime_vec[i];
/*****************/
if (msg_count%update_interval_old==0)
{
CX_Message("OLD: particle Y+ = %f , remaining eddy-
lifetime = %g ( = %3.2f \% ), dt = %g\n",Y_plus_p,p->eddy_time,100*p->eddy_time/p-
>user[T_EDDY],P_DT(p));
CX_Message("OLD: u' = %f , v' = %f , w' = %f\n",p-
>V_prime[0],p->V_prime[1],p->V_prime[2]);
msg_count++;
}
}
msg_count++;
for (i=0; i<dim; i++) c->V[i] += p->V_prime[i];/* add fluctuating vel. to mean cont.
phase vel*/
#if 0
CX_Message("u = %f ; v = %f ; w = %f ; t = %f ; Y_plus_p = %f ; Y =
%f\n ",c->V[0], c->V[1],c->V[2],p-
>state.time,Y_plus_p,Y_plus_p*C_MU_L(c0,t0)/sqrt(C_UDMI(c0,t0,TAU_WALL)*C_R(c0,t0)));
#endif
}
/* standard stochastic tracking outside of boundary layer */
if (do_stochastics)
{
if (p->eddy_time <= DPM_SMALL)
{
check_time = TRUE;
p->user[9] ++;
p->eddy_time = C_MU_FACTOR(p->injection) * c->tke / c->ted; /* Eq. 23.2-
23 User's Guide, fluid
Lagrangian integral time */
if (random_eddy_lifetime)
p->eddy_time *= -log(dpm_uniform_random(p)); /* Eq. 23.2-32 User's
Guide, characteristic eddy lifetime */
else
p->eddy_time *=2; /* Eq. 23.2-31 User's Guide, characteristic eddy
lifetime */
if (p->state.time == 0.0) p->eddy_time *= dpm_uniform_random(p);
if (p->eddy_time < te_min)
p->eddy_time = te_min;
/* CX_Message("tke %e, ted %e, eddy_time %e\n", c->tke, c-
>ted, p->eddy_time); */
/* velocity fluctuations are calculated using a Gaussian deviate
random number */
if (sg_rsm)
{
real random[3];
real chc11=0., chc12=0., chc13=0.;
real chc22=0., chc23=0., chc33=0.;

111
 chc11 = sqrt(C_RUU(c0,t0));
chc12 = C_RUV(c0,t0)/chc11;
chc22 = sqrt( MAX(SMALL, C_RVV(c0,t0)-chc12*chc12) );
#if RP_2D
if (sg_swirl)
{
#endif
chc13 = C_RUW(c0,t0)/chc11;
chc23 = (C_RVW(c0,t0)-chc12*chc13)/chc22;
chc33 = sqrt( MAX(SMALL, C_RWW(c0,t0)-
chc13*chc13-chc23*chc23) );
#if RP_2D
}
else
chc33 = sqrt(C_RWW(c0,t0));
#endif
for (i=0; i<3; i++) random[i] = dpm_gauss_random(p);
p->V_prime[0] = random[0] * chc11;
p->V_prime[1] = random[0] * chc12 + random[1] * chc22;
#if RP_2D
if (rp_axi)
#endif
p->V_prime[2] = random[0] * chc13 + random[1] * chc23
+random[2] * chc33;
#if RP_2D
else
p->V_prime[2] = 0.;
if (rp_axi_swirl)
{
/* tracking done in 3-D, rotate fluctuation
velocity to particle */
/* position. No need for this if isotropic
fluctuations */
real R = sqrt(SQR(p->state.pos[1]) +
SQR(p->state.pos[2]));
real over_R =
R/SQR(MAX(R,DPM_SMALL));
real Vr=p->V_prime[1];
real Vt=p->V_prime[2];
p->V_prime[1] = (Vr*p->state.pos[1] -
Vt*p->state.pos[2])*over_R;
p->V_prime[2] = (Vr*p->state.pos[2] +
Vt*p->state.pos[1])*over_R;
}
#endif
}
else
{
Vprime_rms = sqrt(2./3. * c->tke); /* Eq. 23.2.27 User's
Guide, RMS fluctuating components */
for (i=0; i<3; i++) p->V_prime[i] = dpm_gauss_random(p) *
Vprime_rms; /* Eq. 23.2.28-30 User's Guide, velocity fluctuations */
}

112
 }
else
{
/* rescale the remaining portion of eddy life time and rms values */
p->eddy_time *= (p->old_ted/p->old_tke) * (c->tke/c->ted);
for (i=0; i<dim; i++) p->V_prime[i] *= sqrt(c->tke/p->old_tke);
check_time = FALSE;
}
p->old_tke = c->tke;
p->old_ted = c->ted;
for (i=0; i<dim; i++) c->V[i] += p->V_prime[i];/* add fluctuating vel. to mean
cont. phase vel*/
}
Vrel = 0.0;
for (i=0; i<dim; i++)
Vrel += SQR(c->V[i]-pv[i]);
Vrel = sqrt(Vrel);
if (rp_visc)
p->Re = c->rho * P_DIAM(p) * Vrel / c->mu;
else
p->Re = 0.;
coeff = SphereDragCoeff(p->Re);
if (P_DIAM(p) != 0.0)
drag_factor = coeff * c->mu / ( P_RHO(p) * P_DIAM(p) * P_DIAM(p));
else
drag_factor = 1.;
if (check_time && (do_stochastics||do_wall_stoch) ) /* !!! */
{
real drag_f = MAX(drag_factor, ACCURACY);
/* limit the eddy_time by the eddy crossing time too */
real Le = C_MU_FACTOR(p->injection) / 1.225 * pow(c->tke,1.5)/ c->ted;
if (Le < Vrel/drag_f)
p->eddy_time = MIN(p->eddy_time,-2*log(1. -
Le/(Vrel/drag_f))/drag_f); /* Eq. 23.2.33 User's Guide, particle eddy crossing time */
/*CX_Message("std.-tracking: TL = %g \n",p->eddy_time);*/
}
return coeff;
}
/* set particle trajectory integration time step size */
DEFINE_DPM_TIMESTEP(dpm_ts_by_eddy_time,p, dt)
{
#if !PARALLEL
dt_min = RP_Get_Real("delta_t_min");
te_min = RP_Get_Real("t_eddy_min");
#endif
if (do_stochastics)
{
if (0) /* (msg_count%update_interval==0) */
{
CX_Message("dt %e, eddy_time %e\n", dt, p->eddy_time);
msg_count++;
}
if (dt>p->eddy_time)

113
 return p->eddy_time;
}
if (do_wall_stoch)
{
if (0) /* (msg_count%update_interval==0) */
{
CX_Message("dt %e, eddy_time %e\n", dt, p->eddy_time);
msg_count++;
}
/*if (dt>p->eddy_time) return p->eddy_time;*/
dt = p->user[T_EDDY]/10.;
if (dt < dt_min) dt = dt_min;
}
return dt;
}
/* name particle scalars for post-processing */
DEFINE_DPM_SCALAR_UPDATE(dpm_sc_interact,cell,thread,initialize,p)
{
cell_t c = RP_CELL(&p->cCell); /* Get Cell and Thread from */
Thread *t = RP_THREAD(&p->cCell); /* Particle Structure using new macros*/
real pp[3] = {0.,0.,0.}; /* particle position vector */
NV_V(pp, =, p->state.pos); /* position of particle, x,y,z-coord. */
if (initialize)
{
/* this is the initialization call */
p->user[E_PLUS] = 0.; /* e_plus */
p->user[T_EDDY] = 0.; /* eddy lifetime */
p->user[Y_PLUS_P] = get_yplus_particle(c,t,pp); /* y_plus value at particle location */
p->user[LAMBDA1] = dpm_gauss_random(p);
p->user[LAMBDA2] = dpm_gauss_random(p);
p->user[LAMBDA3] = dpm_gauss_random(p);
p->user[T_FORCE_X] = 0.; /* thermophoretic body force x */
p->user[T_FORCE_Y] = 0.; /* thermophoretic body force y */
p->user[T_FORCE_Z] = 0.; /* thermophoretic body force z */
p->user[N_EDDY] = 0.; /* number of eddies */
p->user[NEW_FORCE] = 0.; /* flag for thermophoretic model */
strcpy(user_particle_vars[E_PLUS].name,"e_plus");
strcpy(user_particle_vars[E_PLUS].label,"e_Plus");
strcpy(user_particle_vars[T_EDDY].name,"t_eddy");
strcpy(user_particle_vars[T_EDDY].label,"t_Eddy");
strcpy(user_particle_vars[Y_PLUS_P].name,"y_plus");
strcpy(user_particle_vars[Y_PLUS_P].label,"Y_Plus");
strcpy(user_particle_vars[LAMBDA1].name,"lambda_1");
strcpy(user_particle_vars[LAMBDA1].label,"Lambda_1");
strcpy(user_particle_vars[LAMBDA2].name,"lambda_2");
strcpy(user_particle_vars[LAMBDA2].label,"Lambda_2");
strcpy(user_particle_vars[LAMBDA3].name,"lambda_3");
strcpy(user_particle_vars[LAMBDA3].label,"Lambda_3");
strcpy(user_particle_vars[T_FORCE_X].name,"t_force_x");
strcpy(user_particle_vars[T_FORCE_X].label,"t_Force_x");
strcpy(user_particle_vars[T_FORCE_Y].name,"t_force_y");
strcpy(user_particle_vars[T_FORCE_Y].label,"t_Force_y");
strcpy(user_particle_vars[T_FORCE_Z].name,"t_force_z");

114
 strcpy(user_particle_vars[T_FORCE_Z].label,"t_Force_z");
strcpy(user_particle_vars[N_EDDY].name,"n_eddy");
strcpy(user_particle_vars[N_EDDY].label,"N_Eddy");
strcpy(user_particle_vars[NEW_FORCE].name,"new_force");
strcpy(user_particle_vars[NEW_FORCE].label,"New_Force");
}
else
{
/* not needed */
}
}
/* UD boundary condition for trapping particles on a wall and mark position of impact */
…

A.4 Local Flow Property Injection UDF

/*********************************************************************/
/* UDF that initializes particle injection properties */
/*********************************************************************/

#include "udf.h"
#include "cxiface.h"

DEFINE_DPM_INJECTION_INIT(injection_velocity,I)
{
Particle *p;

loop(p,I->p) /* Standard Fluent Looping Macro to get

particle streams in an Injection */
{
cell_t c;
Thread *t;
CX_Cell_Id cx_cell;
cx_cell=*CX_Find_Cell_With_Point(P_POS(p));
c=RP_CELL(&cx_cell);
t=RP_THREAD(&cx_cell);

P_VEL(p)[0] = C_U(c,t);
P_VEL(p)[1] = C_V(c,t);

#if RP_3D
P_VEL(p)[2] = C_W(c,t);
#endif

P_T(p)=C_T(c,t);
}
}

115