Course Hero has millions of student submitted documents similar to the one

below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

Course Hero has millions of student submitted documents similar to the one below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

and Development application of a dispersed two-phase flow capability in a general multi-block Navier Stokes solver
Anant P. Shah
Thesis submitted to the faculty of the Virginia Polytechnic Institute and State University in partial fulfillment of the requirements for the degree of
Master of Science In Mechanical Engineering
Danesh K. Tafti - Chair Ishwar K. Puri Pavlos P. Vlachos
December 6, 2005 Blacksburg, VA
Keywords: Dispersed Phase, Lagrangian Particle Tracking, Sand, Internal Cooling Ribbed Duct, Channel Flow, Large Eddy Simulation
Copyright 2005, Anant Shah
Development and application of a dispersed two-phase flow capability in a general multi-block Navier Stokes solver Anant P. Shah
ABSTRACT Gas turbines for military applications, when operating in harsh environments like deserts often encounter unexpected operation faults. Such performance deterioration of the gas turbine decreases the mission readiness of the Air Force and simultaneously increases the maintenance costs. Some of the major factors responsible for the reduced performance are ingestion of debris during take off and landing, distorted intake flows during low altitude maneuvers, and hot gas ingestion during artillery firing. The focus of this thesis is to study ingestion of debris; specifically sand. The region of interest being the internal cooling ribbed duct of the turbine blade. The presence of serpentine passages and strong localized cross flow components makes this region prone to deposition, erosion, and corrosion (DEC) by sand particles. A Lagrangian particle tracking technique was implemented in a generalized coordinate multi-block Navier-Stokes solver in a distributed parallel framework. The developed algorithm was validated by comparing the computed particle statistics for 28 microns lycopodium, 50 microns glass, and 70 microns copper with available data [2] for a turbulent channel flow at Re=180. Computations were performed for a particle-laden turbulent flow through a stationary ribbed square duct (rib pitch / rib height = 10, rib height / hydraulic diameter = 0.1) using an Eulerian-Lagrangian framework. Particle sizes of 10, 50, and 100 microns with response times (normalized by friction velocity and hydraulic diameter) of 0.06875, 1.71875, and 6.875 respectively are considered. The calculations are performed for a nominal bulk Reynolds number of 20,000 under fully developed conditions. The carrier phase was solved using Large Eddy Simulation (LES) with Dynamic Smagorinsky Model [1]. Due to low volume fraction of the particles, one-way fluid-particle coupling was assumed. It is found that at any given instant in time about 40% of the total number of 10 micron particles are concentrated in the vicinity (within 0.05 Dh) of the duct surfaces, compared to 26% of the 50 and 100 micron particles. The 10 micron particles are more sensitive to the flow features and are prone to preferential concentration more so than the larger particles. At the side walls of the duct, the 10 micron particles exhibit a high potential to erode the region in the vicinity of the rib due to secondary flow impingement. The larger particles are more prone to eroding the area between the ribs and towards the center of the duct. At the ribbed walls, while the 10 micron particles exhibit a fairly uniform propensity for erosion, the 100 micron particles show a much higher tendency to erode the surface in the vicinity of the reattachment region. The rib face facing the flow is by far the most susceptible to erosion and deposition for all particle sizes. While the top of the rib does not exhibit a large propensity to be eroded, the back of the rib is as susceptible as the other duct surfaces because of particles which are entrained into the recirculation zone behind the rib.
iii
Acknowledgement This thesis has been the most rewarding experience of my academic career. I always wondered how in this materialistic world, can an individual sacrifice many of the worldly-pleasures and work for something that he or she is not sure of being useful to the world. This experience has helped me realize that research is more than solving equations and running simulations. It is the hope of an individual that his or her work will help some Einstein of our century to better understand Mother Nature. Even though all of us don't have the patience of Edison, our work can count for one of the 10,000 trials. I take this opportunity to show respect for all the Ph.D.'s and the professors who are working day and night to make our planet a cheaper and safer place to live in. The peaceful nights of many people are the result of your sleepless nights in the research labs. I realized that one of the reasons for the quality life in USA is its encouragement to science and innovation. This realization has come through my interaction with lot of people who have supported and encouraged me to give my best. I take this opportunity to thank all the people who have helped me directly or indirectly. First, I would like to thank my advisor Dr. Danesh Tafti. His dedication to achieve excellence and his belief to be a constructive part of the whole have inspired me to put that extra effort and motivated me to keep going in spite of failures. I owe him lots of gratitude for having me shown this way of research. He could not even realize how much I have learned from him. Also, I am indebted to him for his financial help. Although this project was not funded, he made sure that I did not have to bear the burden of paying tuition fees. In addition, during summers he paid stipend from his personal funds to make sure I had the best stay. Lastly, his appreciation of my efforts by gifting me a car is an honor by itself. Dr. Tafti, I have been very fortunate to have an advisor like you. Thanks for everything. I would like to thank my committee members for reviewing my work and providing valuable suggestions. I appreciate your help and time. My sincere thanks to the HPCFD group: Aroon Viswanathan, Ali Rozati, Evan Sewall, Shi-Ming Li, Mohammad Elyyan, Pradeep Gopalkrishnan, and Keegan Delaney. You have been very helpful and cooperative. I would like to specially thank Aroon (for answering all the weird questions and accommodating me in summer (Thanks Ashvin and Vikram)), Evan (for his tips on presentations and writing codes), and Ali (for his help on System X and his free lessons on the "Important" things in life). Thanks you guys. I have been privileged to stay with the coolest people at Virginia Tech: Harshit Shah, Abhijeet Jhala, Hardik Patel, Devdutta Bhosale, and Aditya Ketkar. Thanks guys. I'll miss the tight mocking schedules followed by us. I wish you the very best of life. A special thanks to Harshit (who helped me with my application for the Masters program), Abhijeet, and Hardik for going out of their way to help me. This thesis would not have been possible without the encouragement and support of Hetal (my girlfriend -your passion and zeal for life have taught me how simple it is to be happy), my parents (my career is a result of your sacrifices and belief in me), and Aakash (brother). It is hard to put in words how much I love you all. Finally, I am very grateful to aunt Palak for supporting me in the most challenging phase of my life. Her love and encouragement are invaluable to the completion of this thesis.
iv
TABLE OF CONTENTS Chapter 1: Introduction...............................................................................1 Chapter 2: Numerical Procedure....................................................................7 Chapter 3: Particle Tracking Algorithm.........................................................15 Chapter 4: Turbulent Channel Flow at Re = 180...............................................29 Chapter 5: Internal Cooling Ribbed Duct........................................................46 Chapter 6: Summary and Conclusions............................................................72 References Appendix A: Modified drag coefficient near solid walls.......................................80 Appendix B: Newton's iterative procedure for solving the trilinear function...............81 Appendix C: A pseudo FORTRAN program combining the particle location and cell search technique.....................................................................84 Appendix D: Flowchart Lagrangian Particle Tracking Algorithm.........................86 Appendix E: Wall collision models...............................................................90 Vita
v
LIST OF FIGURES Chapter 3 Figure 3.1: Particle locations at times t and t+t in a general structured 3-D mesh........17 Figure 3.2: Search criteria for hexahedral cells.................................................20 Figure 3.3: Demonstration of the cell search scheme..........................................21 Figure 3.4: Domain decomposition and data structure in GenIDLEST.....................22 Figure 3.5: An example of unstructured block topology. East face (+) of one block exchanges information with the north face (+) of an adjoining block. The axes are also arbitrarily oriented with each other...............................................................23 Figure 3.6: Decomposition of a physical domain in a multi-block framework. The numbers represent the block and processor numbers. It is assumed that each processor has one block and that the processor number is same as the block number.....................24 Chapter 4 Figure 4.1: Computational domain non-dimensionalized by channel half width ()......39 Figure 4.2: Time evolution of the mean flow velocity on the top and bottom wall........39 Figure 4.3: RMS Fluctuations of normal turbulent stresses (b) turbulent shear stress.....39 Figure 4.4: Comparison of mean streamwise velocity predicted by current model with particle statistics by Wang and Squires, 1996.Particles compared: 70m copper, 50m glass, and 28m lycopodium......................................................................40 Figure 4.5: Mean stream-wise velocity (a). Lycopodium (b). Glass (c). Copper ..........41 Figure 4.6: Comparison of velocity fluctuations for (a) Lycopodium (b). Glass (c). Copper...........................................................................................42 Figure 4.7: Number density profile (a). Lycopodium (b). Glass (c). Copper...............43 Figure 4.8: Number of particles colliding with the wall as a % of the total number of particles (a). Lycopodium (b). Glass (c). Copper...............................................44
vi
Figure 4.9: Instantaneous particle at t+ = 6, 0 < y+ <5, Re = 180 (a). Lycopodium (b).Glass(c).Copper.................................................................................45 Figure 4.10: Instantaneous particle at t+ = 12, 170 < y+ < 180, Re = 180 (a). Lycopodium (b). Glass(c). Copper...........................................................45 Chapter 5 Figure 5.1: Schematic of the cooling of a modern gas turbine blade........................47 Figure 5.2: Computational domain non-dimensionalized by hydraulic diameter...........61 Figure 5.3: (a) Mean streamline distribution in the z-symmetry plane. Reattachment occurs at 4.1e downstream of rib. The leading edge eddy extends between 0.8-0.9e upstream of rib. (b)Mean lateral or spanwise flow velocity (w b) in the vicinity of the smooth wall..........................................................................................62 Figure 5.4: Instantaneous distribution of sand at t = 27 (Dh / u) a) 10m, b) 50m, c) 100m.............................................................................................63 Figure 5.5: Instantaneous streamwise velocity at a plane z / D h = 0.05 .....................63 Figure 5.6: Number of particles impinging the side wall in 1.3 milliseconds. (a) 10m, (b) 50m, (c) 100m...............................................................................64 Figure 5.7: Impact velocity on the side wall. (a). 10m (b). 50m (c). 100m............64 Figure 5.8: Impact angle on the side wall. (a). 10m (b). 50m (c). 100m...............64 Figure 5.9: Measure of potential for erosion and deposition at side wall ( 106) (a) 10m, (b) 50m, (c)100m...................................................................65 Figure 5.10: Number of particles impinging the ribbed (bottom) wall in 1.3 milliseconds (a). 10m (b). 50m (c). 100m..................................................................65 Figure 5.11: Impact velocity on the ribbed wall. (a). 10m (b). 50m (c). 100m.......66 Figure 5.12: Impact angle on the ribbed wall. (a). 10m (b). 50m (c). 100m..........66 Figure 5.13: Measure of potential for erosion and deposition at bottom wall ( 106). (a) 10m, (b) 50m, (c)100m.....................................................67
vii
Figure 5.14: Number of particles impinging the front surface of the rib in 1.3 milliseconds. (a) 10m, (b) 50m, (c) 100m.................................................67 Figure 5.15: Impact velocity on the front surface of the rib. (a). 10m (b). 50m (c). 100m............................................................................................68 Figure 5.16: Impact angle on the front surface of the rib (a). 10m (b). 50m (c). 100m............................................................................................68 Figure 5.17: Measure of potential for erosion and deposition on front surface of the rib ( 106). (a) 10m, (b) 50m, (c)100m.....................................................69 Figure 5.18: Measure of potential for erosion and deposition on top surface of the rib ( 106). (a) 10m, (b) 50m, (c)100m.....................................................69 Figure 5.19: Number of particles impinging the back surface of the rib in 1.3 milliseconds. (a) 10m, (b) 50m, (c) 100m.................................................70 Figure 5.20: Impact velocity on the back surface of the rib (a). 10m (b). 50m (c). 100m............................................................................................70 Figure 5.21: Impact angle on the back surface of the rib (a). 10m (b). 50m (c). 100m............................................................................................71 Figure 5.22: Measure of potential for erosion and deposition on back surface of the rib ( 106). (a) 10m, (b) 50m, (c)100m.....................................................71
viii
LIST OF TABLES
Table 4-1: Channel properties used in the calculations........................................31 Table 4-2: Particle properties for a turbulent channel flow at Re=180.....................31 Table 5-1: Flow and geometrical parameters used in the calculation........................49 Table 5-2: Properties of the carrier phase........................................................49 Table 5-3: Particle properties used in the calculation..........................................50
ix
CHAPTER 1. INTRODUCTION "Particle laden flows are characterized by interaction between two different kinds of matter. The difference between the matters can be their thermodynamic state, called the phase (e.g., gas, liquid, or solid) and/or their chemical composition" [3]. Such flows occur in wide ranging industrial applications. Examples include spray combustion, turbomachinery operating in polluted environments, drug delivery, material processing, filtration, fluidized beds, fluid-particle transport through pipes and ducts, nuclear reactor cooling, and pollution control. In addition, global issues like bio-terrorism, CO2 emissions, and disposal of nuclear waste involve particle-laden flows. In spite of such varied applications these flows are poorly understood. A better understanding of these flows is required to develop devices that can counter the current threats and can help in achieving a cleaner and safer environment [4]. In order to enhance our current knowledge level in particle laden flows, efforts are needed in modeling of fundamental processes like particle dispersion, mixing, concentration, phase changes (evaporation, condensation, boiling, melting), and inter-phase heat, mass, and momentum transfer. The goal of this thesis was to develop a dispersed phase model to understand particle dispersion and concentration in complex turbulent flows, and to apply the model to study particle transport in a stationary internal cooling ribbed duct. Such flows are encountered in gas turbines operating in polluted environments. The primary issues of concern being deposition, erosion, and corrosion (DEC) due to ingestion of high speed particles in the intake air of the compressor. These particles can be ingested during low altitude maneuvers when the inlet of the gas turbine acts like a big vacuum and ingests foreign
1
objects with the intake air. In addition, previous experiments show that a vortex system can develop between the ground and the air inlet which imparts momentum to these particles such that they reach the air intake. Also, particles can get ingested due to a sand storm. In addition, the use of alternate fuels based on the gasification of coal, biomass, and fuel oils introduce corrosive vapors and molten ash, which combined with high temperatures and high gas stream velocities, can cause DEC of turbine components. The extent of damage done by DEC is directly dependent on two factors: (a) impaction of solid or molten particulates on surfaces, which lead to erosion and deposition. (b). Condensation of corrosive vapors on surfaces, which initiates thermal corrosion and deposition. Jensen et al. [5] in their study on deposition in a land based turbine cite the following references on the effect of erosion and deposition on the engine efficiency: A 6-10% loss in adiabatic efficiency was reported by Ghenaiet et al. [6] for an axial fan subjected to 6 hours of sand ingestion. Kim et al. [7] reported a strong decrease in the film cooling effectiveness due to blockage of the film cooling holes by volcanic ash. An interesting behavior was observed by Bons [8]. He found that the increase in surface roughness due to erosion led to an increase in heat transfer by 50% and that in skin friction by 300%.
The objective of this thesis is to provide data on the number of particles impacting at a given location, their average impingement velocity and angles. These results combined with empirical relations correlating material properties and particle impingement data can give the amount of material removed from a certain location in a given operation time of the gas turbine. DEC not only decreases the durability of
2
components but also results in unexpected operation faults. Previous studies [9]-[14] on erosion of turbomachinery are concentrated on studying the compressor sections and the turbine blades. While, experiments by [15] have shown that particles do reach the internal cooling air system of the gas turbine, there is no available data in the literature that can help in understanding this phenomenon. This dictates the need for more detailed particle transport data to analyze the effect of foreign matter on the internal geometries of a gas turbine. Further, this data will be valuable to researchers in developing improved models to study dynamics of particulate matter in the post combustor hot gases and their interaction with the hot gas flow path leading to DEC. In addition to the useful practical applications, particle transport through an internal cooling ribbed duct represents an interesting case for studying particle behavior in anisotropic flows with wall-bounded turbulence.
Previous studies [16]-[21] on particle-laden turbulent flows are focused on isotropic turbulence, plane channels, and pipes. From these studies the following general conclusions are reached:
An inertial bias in the particle trajectories results in a preferential concentration of particles in regions of low vorticity.
Particles with Stokes number closer to Kolmogorov scales exhibit stronger effects of preferential concentration.
Particles tend to accumulate in low speed streaks near the wall.
In addition, studies have been done to understand particle transfer mechanisms in turbulent boundary layers. Marchioli and Soldati [22] examined the reasons for non-
3
uniform concentrations of heavy particles in a turbulent boundary layer by performing a Direct Numerical Simulation (DNS) in a vertical channel flow at Re=150. They inferred that particle transfer was affected by the coherent structures which control the turbulence regeneration cycle near the wall and that sweeps (characterized by inward motion of high speed fluid, i.e., instantaneous streamwise velocity ( u ') > 0, and instantaneous crossstream velocity ( v ' < 0) and ejections (characterized by outward motion of low speed ) fluid, i.e., instantaneous streamwise velocity ( u ') < 0, and instantaneous cross-stream ) velocity ( v ' > 0) were responsible for preferential concentration of particles in the near wall region. However, only few studies have been reported on particle transport in ducts with inhomogeneous turbulence. Madabhushi and Vanka [23] performed a LES study of unladen turbulent flow in a square duct at Re=360. They found that secondary flows were able to convect fluid from the centre of the duct to the near wall region. Winkler et. al. [24] used LES to study the preferential concentration of particles in a fully developed turbulent square duct flow at Re=360. They observed that particle dispersion was affected due to the presence of secondary flows and that particles were seen to accumulate in regions of high strain rate, and low swirling strength. They also noted that particle accumulation in regions of high vorticity increased with an increase in the particle Stokes number. A Direct Numerical Simulation (DNS) was carried out by Sharma [4] to study particle concentration in a straight square duct at Re=300. He concluded that particles with Stokes number closer to Kolmogorov time scales showed no deposition but were seen to accumulate within the turbulent boundary layer. Also,
particle dispersion showed similar behavior as predicted in flows with isotropic turbulence like channel flows. Considering the previous studies [25] and the flow features
4
of an internal cooling duct, like, the presence of strong localized cross-flow components near the side walls with mean spanwise velocities reaching 30% of the mean streamwise velocity, and the presence of eddies in front and behind of the rib, it is expected that particle dispersion may be significantly enhanced as a result of strong secondary flows.
Finally, particle-laden flows are best studied in the Lagrangian frame of reference where the dispersed phase is described by recording the detailed histories of individual particles. The motion of these particles is highly sensitive to the unsteady nature of the flow and the relevant turbulence length scales. The most accurate information about the instantaneous flow structure is obtained from DNS. In DNS, the exact Navier-Stokes equations are solved without any modeling or empiricism. Although highly accurate, DNS is restricted to low Reynolds number flows due to the very fine spatial and temporal grid requirements. In comparison, a less computationally intensive method is LES. LES resolves the large energetic scales and models the smaller scales, which are thought to be universal. Within the last two decades, numerous DNS studies have appeared in the literature; but data spanning much wider ranges of flow Reynolds number are needed for applying the results in industrial settings. The current study provides data for a nominal bulk Reynolds number of 20,000, which is often the case in industrial flow applications.
To the authors knowledge this thesis is the first study on particle transport in an internal cooling ribbed duct. Three challenges were faced in obtaining results for the present work. First, developing the model for a general curvilinear system; second, implementing the model in an existing flow solver Generalized Incompressible Direct
5
and Large-Eddy Simulations of Turbulence (GenIDLEST [26]) with parallel capabilities, and third, applying the particle physics to a stationary internal cooling ribbed duct. The approach to these challenges is discussed in detail as follows. Section 2 describes the numerical procedure for solving the carrier phase and the dispersed phase. Section 3 describes in detail the development of a Lagrangian particle tracking algorithm to solve the dispersed phase and its implementation in the flow solver. Section 4 discusses the validation of the algorithm by comparing the computed particle statistics of 28 micron lycopodium, 50 micron glass, and 70 micron copper with available data [2] for a turbulent channel flow at Re=180. Section 5 presents results on particle impingement (velocity, angle, and number of particles) on the surfaces of a stationary duct. Finally, Section 6 presents a summary and conclusions from the present simulations.
6
CHAPTER 2. NUMERICAL PROCEDURE
This section describes the numerical procedure for solving the carrier phase and the dispersed phase. The carrier phase is solved using Large Eddy Simulation coupled with Dynamic Smagorinsky Model [1] and the dispersed phase is solved using a Lagrangian approach. The chapter proceeds as follows. It starts with a brief description of the governing equations for the carrier phase followed by the approach used for solving the dispersed phase in a Lagrangian framework.
2.1 Carrier Phase
The governing non-dimensional incompressible fluid equations written in general coordinates are as follows: Continuity:
j ( gU j )=0,
(2.1)
Momentum:
t
(
g ui +
)
j
j gU u i = -
j
j g (a )i p + j
1 1 + Re Ret
gg
jk ui k
,
(2.2)
where a i are the contravariant basis vectors (the notation (a j )k is used to denote the k-th component of vector a j .
j
a
= j / x k ,
k
g
is the Jacobian of the transformation, gij
gU j = g a j k uk
are the elements of the contravariant metric tensor,
( )
is the contravariant
7
flux vector, ui is the Cartesian velocity vector. The Reynolds number is defined as
Re = u c Lc
, where uc and Lc are a characteristic velocity and length scale of the flow,
respectively. The overbar denotes grid filtered quantities with an implicit top-hat filter,
G
. Re t , the inverse of the non-dimensional eddy-viscosity, is modeled as:
1 = C s 2 ( g )2/3 S Ret
,
(2.3)
where S is the magnitude of the strain rate tensor given by S = 2 S ik S ik . The strain rate tensor is given by:
S ik = u i u k 1 (a m )k + (a m )i , 2 m m
(2.4)
and the Smagorinsky constant, C s 2 , is obtained via the Dynamic subgrid stress model [1]. The fluid field equations are discretized with a conservative finite-volume formulation using a second-order central difference scheme on a non-staggered grid topology The Cartesian velocities, and pressure are calculated and stored at the cell center, whereas contravariant fluxes are stored and calculated at the cell faces. A projection method is used for the time integration of the discretized continuity and momentum equations. The temporal advancement is performed in two steps, a predictor step, which calculates an intermediate velocity field, and a corrector step, which calculates the updated velocity at the new time step by satisfying discrete continuity. Details about the method can be found in [26], and application to ribbed duct flow in [25].
8
2.2 Dispersed Phase
The dispersed phase is solved using the Lagrangian approach which is described below.
2.2.1 Lagrangian approach
The basic objective of a Lagrangian Particle Model is to integrate the particle equation of motion, described in the Lagrangian frame. The equation of motion for a particle can be expressed by the following set of ODE's:
mp dU ip = dt Fi p
dX ip = U ip dt
(2.5)
where, m p is the mass of the particle and U ip , X ip represent the velocity & position of the particle in the i direction. Fi p is the force acting on the particle in the i direction. A number of forces can act on a finite-inertia particle. The physical significance of various forces is explained below. For a detailed description of the terms the reader should refer to [27].
Fip = Drag force + Added Mass force + History effect + Gravitational force + Buoyancy force + Lift force + Intercollision force + Brownian force + Thermophoresis + Magnus force
The first term on the right hand side of the above equation represents the Stokes drag, which is the force exerted on the particle due to the relative velocity between the fluid and the particle and acts in a direction opposing the relative flow. In most industrial
9
applications, particle motion is dominated by the drag force. The second term describes the acceleration of the fluid near the particle surface from fluid velocity to the particle velocity. Mass of fluid that undergoes this acceleration is called "carried mass" which is equal to one half of the displaced mass of the fluid. The acceleration of the fluid near the particle causes the flow around the particle to differ from that in the steady motion. The force required to maintain the flow pattern was approximated by Basset and is represented by the third term. This force depends on the history of the particle trajectory and hence is called the "history effect". The effect of the added mass and Basset forces is negligible for particles with density substantially large than the fluid density. Saffman's lift force [28] is caused by the shear of the surrounding fluid which results in a nonuniform pressure distribution around the particle. This force assumes non-trivial magnitudes only in the viscous sublayer. McLaughlin [29] showed that even in the viscous sublayer it is an order of magnitude smaller than the normal component of the Stokes drag force. Force exerted due to inter particle collisions is represented by the inter-collision force. It assumes importance for high volume fraction of the dispersed phase. Brownian and thermophoresis forces are important in the study of submicron particles, and are a result of random molecular motion and forces induced due to temperature gradients, respectively. Magnus force is a result of the rotation of the particle which causes differential pressures perpendicular to the flow. The particle equation of motion used for the present study is based on the following assumptions:
The particles are rigid spheres and that they are considered as points located at the center of the sphere. This assumption does induce some inaccuracy for particles sizes
10
which are greater than the near-wall grid spacing. But, considering previous studies [30], [2] we found that it is the best approximation for calculating particle drag.
The particle density is substantially larger than the fluid density. Elgobashi and Truesdell [31] showed that for particles with ( p / f >> 1) , the only significant forces are the Stokes drag, the buoyancy, and the Basset forces and that the Basset force was always an order of magnitude smaller than the drag and buoyancy forces.
The contribution of lift force to particle motion is small and can be neglected. Wang
et al. [32] have shown that neglecting the lift force results in a slight decrease in the
deposition rate.
Due to low volume fraction of the particles inter-particle collisions are negligible. Particles do not affect the fluid turbulence. Experiments by Kulick et al.[33], and Kaftori et al. [34] have shown that for low volume fractions the turbulence modifications are negligible. Also, in the near-wall region where the particle concentration may be locally large, the turbulence intensities are modified by a very small amount and can be neglected [22].
For the particle sizes considered in the present study, subgrid-scales have a negligible effect on particle trajectories. Wang and Squires [35] investigated the effect of subgrid scale fluctuations on the trajectory of particles in a vertical turbulent channel flow for bulk Reynolds number up to 79400. They found that the difference in results obtained with and without the added fluctuating velocities was negligible and that they compared well with the DNS results of McLaughlin [29] Yeh and Lei [36] performed LES calculations using Smagorinsky subgrid-scale model to study particle motion in homogeneous isotropic turbulence. They found that the particle motion is
11
governed by the large scales and obtained good agreement with experimental data. Armenio et al. [37] in their study on the effect of the subgrid scales on the particle motion concluded that particles with finite inertia are less sensitive to the subgird scale fluctuations and that for small filter sizes the difference in statistics obtained with a filtered field and with a DNS field are very small. Considering previous studies the effect of subgrid scale fluctuations on the particle motion is neglected. The filter size chosen for the current study resolves the turbulent scales well below the inertial range as shown in [25].
With the above assumptions, the force balance equation takes the form,
f 3 CD du ip =- u - v u ip - vif dt p 4 dp
(
)
( )
(2.6)
where C D = Drag coefficient, u ip = Particle velocity in the i direction, and vif = Fluid velocity in the i direction. The drag coefficient C D =
24 1 + 0.15 Re0.687 ; given by Clift p Re p
u-vdp
et al. [38], which is valid for a particle Reynolds number
up to 700 was used
for all the simulations. The drag coefficient can be modified near solid walls. This effect can be incorporated in the particle equation by multiplying the drag term by a wallcoefficient C wall . The expressions for the wall-coefficient are provided in Appendix A. It is worth mentioning that the expressions available in the literature are valid for particle Reynolds number within the Stokes limit. No relevant expressions are available for particle Reynolds number higher than the Stokes limit. Hence for the present simulations the effect of modified drag coefficients
12
is
not
taken
into
consideration.
Non-dimensionalizing equation (2.6) by characteristic velocity and length scale (uc, Lc ) and using "+" for non-dimensional variables the particle equation of motion becomes:
du ip + dt
+
= -
1
+ p
(1 + 0 . 15 Re
0 . 687 p
)(u
i
p+
- v if +
)
(2.7)
where p is the particle response time given by:
2 d pS
p =
18
, S=
p f
p+ , ui =
p p u ip f + v if + , vi = , p = = = St p uc uc Lc / u c f
f
The ratio p /
is the Stokes number for a given particle and it gives an
estimate as to how quickly a particle responds to the flow structure. For Stokes numbers much less than unity, the particle responds almost instantaneously to the carrier phase and hence at a given instant in time, the particle and fluid velocities are nearly identical. The converse is true for Stokes number much greater than unity. In this case, the particle carries its own inertia and takes a long time to adjust to external conditions. Between these two extremes are particles that respond to some of the flow structures. Hence for the same flow, particles with different Stokes number exhibit different concentrations. Therefore the right characterization of the Stokes number is critical as it largely determines the particle behavior in a specified flow field.
The governing equations of particle motion are integrated using the standard third order Adams-Bashforth scheme described below:
n +1 = n + t
23 n 16 n -1 5 n - 2 - + 12 12 12
13
A fixed time step is used to advance the particles. It is equal to the time step used to advance the Eulerian flow field. In order to numerically integrate the particle equation of motion, fluid velocities are needed at the particle location. Second order Lagrange polynomials are used to interpolate fluid velocities at the particle location.
14
CHAPTER 3. PARTICLE TRACKING ALGORITHM
The previous chapter described the governing equations for the carrier phase and the dispersed phase. In order to integrate the particle equation of motion, it is necessary to know the current location of the particle. This demands an efficient algorithm for searching and locating particles in complex geometries. There are two types of particle tracking algorithms available in the literature, i.e., C-space algorithms, and P-space algorithms. C-space (computational space) algorithms transform the curvilinear grid in physical domain to a cartesian grid in the computational space using Jacobian transformations. In comparison, P-space (physical space) algorithms avoid grid transformation but require complex point location and interpolation techniques. A study by [39] shows that the P-space algorithms are more efficient and accurate than the Cspace algorithms. Hence, a P-space algorithm was developed to search and locate the particles. This chapter proceeds as follows. The first section discusses the objective of a particle search and locating technique followed by a detailed description of the particle tracking algorithm. Finally, the implementation of the algorithm in a parallel multi-block framework is discussed.
3.1 Objective
The objective of a particle tracking algorithm is to describe the dispersed phase by recording the detailed histories of individual particles and it does so by following their trajectories. Figure 3.1depicts the motion of one such particle along an arbitrary trajectory
15
as it traverses through a hexahedral mesh within a tracking time step t. There are two possibilities for this particle. It can either remain in the same cell or can traverse to any of the adjoining cells. The probability of the particle leaving its previous cell is governed by the local flow features and the forces acting on the particle, and is not known a priori. Considering a 3-D grid, the leaving particle can move to any of the 27 surrounding cells at the next time-step, assuming a given particle traverses no more than one cell in time t. But, near the boundaries where the mesh is extremely fine, particles may traverse more than one cell in a given time step. This increases the number of cells a particle can move to in a single time step. A number of particle tracking algorithms try to avoid this situation by using an adaptive technique where the time-step is chosen such that the particle spends at least 5-10 time steps in a given cell. Although these techniques may limit the number of possibilities, they may incur additional computational costs due to small time-steps. Our objective is to develop the model such that it locates the particle in minimum number of searches and uses an efficient time-step for solving the carrier and the dispersed phase. Finally, it is important to note that in reality the physical domain consists of arbitrary shaped hexagonal cells and that the cells shown in Figure 3.1 are for simplicity of explanation. From the above discussion we can conclude that the particlelocating algorithm must answer two questions:
First, is the particle inside or outside the cell it was in at the previous time step? Second, if it has moved out of its previous cell within a time step t, in which control volume can the particle possibly reside? This helps in reducing the number of searches for a given particle.
16
Considering, the tracking of hundreds of thousands of particles for hundreds of thousands of time-steps, the particle tracking algorithm has a direct impact on the size and complexity of two-phase problems that can be modeled with the available resources [40]. West North
t t+t
East
Low
High Figure 3.1: Particle locations at times t and t+t in a general structured 3-D mesh.
South
3.2 Point location in a hexahedral cell
We now turn our attention to the first step of our particle tracking algorithm, i.e., how to know if the particle is in or out of a eulerian cell? Consider a domain R 3 decomposed into N* hexahedral elements such that for
l m , the following exists:
l m =
( null - set ) no intersecti on [41] one common surface
For a given particle x p = ( x p , y p , z p ) in we wish to know whether x p is in l . Let, p1 , p 2 , p 3 , p 4 , p 5 , p 6 , p 7 , p 8 be the eight vertices of a hexahedral cell, where,
*
N = ni x n j x nk , where ni = number of cells in i direction, n j = number of cells in j direction, nk = number of cells in k direction
17
p1 = x(i, j , k ) p 2 = x(i + 1, j, k ) p 3 = x(i, j + 1, k ) p 4 = x(i + 1, j + 1, k ) p 5 = x(i, j, k + 1) p 6 = x(i + 1, j, k + 1) p 7 = x(i, j + 1, k + 1) p 8 = x(i + 1, j + 1, k + 1).
We use the widely known trilinear interpolation function which determines the location of a point x p = (x p , y p , z p ) in physical space from a given natural (non-dimensional) coordinate ( , , ) . The function is defined as follows:
x p ( , , ) =
(p1 (1 - ) + p 2 (1 + )) (p 3 (1 - ) + p 4 (1 + )) (p 5 (1 - ) + p 6 (1 + )) (p 7 (1 - ) + p 8 (1 + ))
(1 - ) + (1 - ) + (1 + ) (1 - ) + (1 + ) (1 + )
/8
(3.1)
Since the above function consists of non-linear products, a direct inversion is not possible. It needs to be solved numerically using an iterative scheme such as NewtonRaphson method; the details of which are explained in Appendix B. To locate a particle at (xp, yp, zp) in a given cell; ( , , ) [- 1,1] should be satisfied. For example,
x p ( , , )
x p ( , , )
0 ,0 , 0
= x centroid of l .
1,1,1
= p8 ; = p1 , and so on.
x p ( , , )
-1, -1, -1
18
Once the particle is located in hexahedron l , the value of any scalar ( ) known at the vertices of l can be obtained at the particle location by,
x
p , y p ,z p
=
( (1 - ) + (1 + )) ( (1 - ) + (1 + )) ( (1 - ) + (1 + )) ( (1 - ) + (1 + ))
(1 - ) + (1 - ) + (1 + ) (1 - )+ (1 + ) (1 + )
/8
(3.2)
3.3 Cell search scheme
The previous section discussed the first aspect of our particle tracking algorithm, i.e., to find out if the particle is inside or outside the hexahedral cell. It remains to establish a criterion to choose the control volume to apply the above technique. It is comparatively easy to identify the correct cell when the particles are injected into the domain as the particles are usually introduced at the boundaries. The problem arises when particles cross cell or block boundaries. As explained in section 3.1 a particle with an arbitrary (i,j,k) location at time t can move to any one of the surrounding 27 cells at time t+t. Hence an efficient cell-search scheme is required which identifies the correct eulerian cell in a minimum number of searches. The point location technique described in section 3.2 gives a good indication for the next search cell if the particle being tracked has moved out of the previous cell location. It is worth mentioning that the particle search starts from the cell in which the particle resided at the previous time step. Then, based on the values of the unknowns , , and , the search criterion advances to the next cell. The search criterion is described below.
19
-1 1
-1 1 -1 1 1 1 1
Go to cell i-1, j+1, k
-1 -1 1
Go to cell i, j+1, k
-1 1 -1 1
Go to cell i+1, j+1, k
1 -1 1
Go to cell i-1, j, k
-1 -1 -1
i, j, k Go to cell i, j-1, k
-1 1 -1 -1
Go to cell i+1, j, k
1 -1 -1
Go to cell i-1, j-1, k
Go to cell i+1, j-1, k
Figure 3.2: Search criteria for hexahedral cells Figure 3.2 shows the search criteria used to advance to the next cell when a particle crosses a cell boundary. As may be observed, the decision for the next cell depends on the magnitude and sign of the unknowns. It is important to clarify that for simplicity of explanation it has been assumed that the particle does not cross the high (+) or low (-) faces. In reality, similar rules apply for a particle traversing through these faces. For example, for a given particle, if -11, -11, and 1 is obtained then it indicates that the particle has crossed the high face of the cell and hence the next cell searched will be (i, j, k+1). Similarly, if 1,1,1 then (i+1, j+1, k+1) cell will be chosen.
20
9
t +t
10
11
12 t+t
5
t
6 t 2
7 3
8 4
1
Figure 3.3: Demonstration of the cell search scheme.
The cell search scheme can be further explained by an example as shown in Figure 3.3. The motions of the particles (black and blue) have been exaggerated to demonstrate the scheme. The numbers in green are the cell numbers. Suppose the black particle moves from cell 1 to cell 10 in time t. The algorithm will start searching for the black particle in cell 1 because its previous location is known to be cell 1. But, the point location scheme will give the values of the unknowns' , and to be greater than one. According to Figure 3.2, cell 6 will be chosen for searching the particle. Using the point location technique on this cell will indicate that -1 1 , and 1. The cell search scheme will advance to cell-10, where the particle will be located. Similarly, the search for the blue particle will follow the path marked by dotted blue arrows, i.e., 5-10-11-12. Most of the particles get located in one or two searches. The situation gets complicated when a particle crosses a block boundary. This is discussed in the next section. A pseudo FORTRAN program combining the particle location and cell-search technique is given in Appendix C.
21
3.4 Implementation of the particle tracking algorithm in a multi-block framework This section discusses the implementation of the above techniques in an existing multi-block Eulerian solver; in our case GenIDLEST [26]. The section starts by giving a brief introduction on the multi-block framework adopted by GenIDLEST, followed by a detailed discussion of the primary issues in the integration of the particle tracking algorithm within the adopted framework.
3.4.1 Multi-block framework of GenIDLEST
Processor domain Physical domain Mesh blocks Virtual Cache blocks
Figure 3.4: Domain decomposition and data structure in GenIDLEST [26].
GenIDLEST is a finite volume code which solves the coupled momentum-energy equations using an overlapping; multi-block, structured mesh topology. Figure 3.4 shows the multi-block framework used by GenIDLEST. Following are its features:
The physical complexity of the geometry is handled by decomposing the domain into smaller units called blocks. In each block, the equations are mapped from physical (x ) to logical/computational space ( ) by a boundary conforming transformation x = x ( ) , where x = ( x, y , z ) and = ( , , ) . Inter-block connectivity can be structured or unstructured. In a structured inter-block topology an East (+) face adjoins to a West (-) face, whilst in an unstructured inter-block topology a 22
face boundary can adjoin a or face with arbitrary axes orientations. Figure 3.5 illustrates the unstructured topology, in which the + (east) face of one block interfaces with a + (north) face of the adjoining block. Boundary information between adjoining faces undergoes appropriate axis translations and rotations to align the data consistently with the local coordinate system when passed from one block to the other.
Axis 2
Axis 1
BLOCK 1
Axis 1
BLOCK 2
Axis 2
East Face
North Face
Figure 3.5: An example of unstructured block topology. East face (+) of one block exchanges information with the north face (+) of an adjoining block. The axes are also arbitrarily oriented with each other [26].
The adopted framework provides an opportunity to extract parallelism at multiple levels. Depending on the total number of mesh blocks and the degree of parallelism sought, each processor can have multiple blocks residing on it as shown in the Figure 3.4. Given the number of blocks on a processor, arrays are dimensioned (ni,nj,nk,mb), where ni, nj, nk are the number of nodes in , , and directions, and
mb is the number of blocks on that processor. It should be noted that the current
code uses a distributed memory network where each processor carries its own data and has no access to the data on other processors. Also, a given processor can send
23
data only to its neighboring processor, i.e., processor on its immediate east, west, north, south, high, and low. Similar restriction holds true for data transfer between blocks. Finally, the block numbers on a particular processor are numbered locally, i.e., if two processors have two blocks each, then the local block numbers on each processor will be 1 and 2, while the global block numbers will be 1, 2, 3, and 4. A mapping is provided to convert a local block number to a global number. Such a description results in efficient use of memory.
Open boundary
7 8 2 1 9 3 10 4 6 5
Wall boundary Particle injection
Figure 3.6: Decomposition of a physical domain in a multi-block framework. The numbers represent the block and processor numbers. It is assumed that each processor has one block and that the processor number is same as the block number.
The above framework is further explained by an example in Figure 3.6. Figure 3.6 shows the division of a 2D domain into 10 blocks. For simplicity of explanation we will consider a 2D case. The numbers represent the block numbers. Depending on the number of processors to be used for the calculation, a corresponding number of blocks are assigned to a processor. Assume that 10 processors are used for the calculation and that
24
each processor has one block. Further, for the sake of discussion assume that the processor numbers are the same as the block numbers. The distributed memory implementation requires that processor 1 can send data only to processor 2; viz., its neighbor on the north. On the contrary, processor 9 can send data to any one of the three neighboring processors; viz., processor 8 on the west, 10 on the east, and 3 on the south. Now, given the current framework following were the primary issues considered while developing the particle tracking algorithm: 1. 2. How to distribute the particles at the beginning of the calculation? What happens when a particle crosses a block boundary? Where should the particle data be sent in order to locate the particle in a minimum number of searches?
3.4.2 Initial distribution of particles
Suppose N particles need to be distributed in the domain. These particles are distributed by inputting the number of particles in a given block (N
block).
A sum of the
number of particles on each block gives the total number of particles N. This gives a flexibility of injecting the particles only at the inlet or distributing them in the entire domain. For a homogenous distribution of particles a random number generator is used to distribute the N block particles within the number of cells in a given block. Mathematically, the particles are distributed as follows: N block = N particles / cell in a block * (ni-1) * (nj-1) * (nk-1) N = N block * total blocks
25
3.4.3 Particle crossing a block boundary
During the particle evolution, situations arise where a particle crosses the block boundary and enters into one of its neighboring or diagonal blocks. When this happens, the particle search algorithm will not find the particle in its current block. At this point the destination block is also not known. Therefore it remains to decide to which block the particle information should be sent so that it can be located in a minimum number of searches. When a particle crosses a block boundary the following steps are performed: First, it is checked if the particle was located in a corner cell at the previous time step. If it was not located at any of the corners then its information is passed to the block adjacent to the face through which the particle is assumed to cross. For example, when the blue particle is not found in block 5 its information is first sent to block 6, as it has the largest possibility of having the particle. Now, suppose that the particle is located at one of the corners. For example, the orange particle in Figure 3.6. This particle can go to any of the blocks shown by the arrows. A judicious choice of where to check first is made by considering the values of , (2D case), and the component of particle velocities in the respective directions. The following questions help in making an efficient decision: 1. Depending on the values of , and , is it possible to exclude any of the blocks linked to + (East), - (West), + (North), - (South) faces? First, if the value of
, or [-1,1] then it cannot pass through that face. Second, check if any of the
violated conditions gives a wall boundary. For example, consider the orange particle trajectory (dotted) from block 3 to block 1. A check on the previous location of the orange particle in block 3 would give values of the unknowns as,
26
-1, and -1. The condition -1 indicates that the particle has crossed the
South face. But, it's known that the south face of block 3 is a wall boundary and hence particles cannot cross that face. This helps in taking an efficient decision by reducing the number of blocks available to a particle. 2. After the above checks, the following terms are calculated: X1 = Abs( * up), X2 = Abs( * vp) The maximum violator of X1 and X2 is chosen to decide on the block to send the particle information. Hence, if X1 > X2, then the particle information is sent to the block linked to the sign () face, i.e., for > 1 - East face, and for < -1 - West face is chosen. Once the block number is known, the corresponding processor number is
obtained through the mapping between global and local block numbers on a processor. However, before sending the particle data, a check is performed to find out whether the particle left through a periodic face. If it did leave through a periodic face then the following operation is performed:
Let the particle coordinates be x p , y p , z p before sending the data. Considering the
periodic lengths in the x, y, and z directions, as Lx , L y , and Lz respectively, the particle coordinates are readjusted as follows.
x p = x p Lx y p = y p Ly z p = z p Lz
The periodic length is added when a particle leaves from the west, south, or low faces while, it is subtracted when the particle leaves from the east, north, or high faces.
27
Finally, it is important to know, when a particle moves out of the domain. The answer is simple. If the particle is not located in the block it was at its previous time step and if the values of the , suggest that there is no block linked to that face then the particle is considered to be out of the domain. Its information is recorded in a lost file and then removed from the calculation. For example, the green particle in block 7. Consider that it moved out of the domain in a tracking time step t. A check on its previous location would give a value of 1. This implies that the particle left through the north face. But, it is known that there is no block linked to the north face of block 7. In such a situation the particle is considered out of the calculation domain. A flowchart of the entire algorithm is provided in Appendix D.
28
CHAPTER 4. TURBULENT CHANNEL FLOW AT Re = 180
The objective of this chapter is to validate the developed Lagrangian particle model in a turbulent shear flow for which computational results exist. The accuracy of the model is assessed by comparing the computed particle statistics with available data [2] for 70 m copper, 50 m glass and 28 m lycopodium particles for a turbulent channel flow at Re = 180 based on a characteristic length (Lc = ) and velocity scale (uc = u). Here is the channel half-width and u is the friction velocity.
4.1 Computational details
The computational domain is shown in Figure 4.1. The 3-D channel extends 2 in X-direction, 2 in Y-direction, and in Z-direction, where is the half-channel width. No slip boundary conditions are applied at the channel walls with periodic boundary condition in the homogeneous X and Z directions. The flow field was calculated using 64 x 64 x 64 grid points in the x, y, and z directions respectively. The grid spacing used in wall coordinates was x+ = 17.67 and z+ = 8.84 in the
homogeneous wall-parallel directions, and y+ = 0.9 in the inhomogeneous wall-normal direction. The maximum spacing at the center of the channel was y+ = 6. The non-
dimensional time step used was 1 x 10-3. After specifying a fixed mean pressure gradient in the flow direction, the calculation was initialized with a perturbed mean turbulent channel flow profile, which was allowed to evolve in time till the solution reached a stationary state. The time evolution of the stationary mean flow velocity at the top and
29
bottom walls is shown in Figure 4.2. Suffice to say that, for Re = 180, 64 x 64 x 64 grid points are enough to capture the particle dynamics near the wall. This can be corroborated from Figure 4.3 which shows a comparison of the RMS fluctuations of normal turbulent stresses and turbulent shear stresses (normalized by friction velocity) calculated by the current code with previously published Direct Numerical Simulation (DNS) results of Kim et al. [42]. As may be observed, very good agreement is obtained with the DNS data in wall coordinates. The slight under-prediction of the peak values is attributed to the fact that Kim et al. [42] use spectral functions with a much higher resolution than the current study.
4.2 Properties
Table 4-1and Table 4-2 respectively; describe the channel properties and the particle properties used in the simulation. These properties are identical to those used by Wang and Squires [2]. The particle Stokes number used for the current simulations is nondimensionalized by the outer scales, i.e., u and . The particle Stokes number can also be defined using the inner scales, i.e., u and . Such scaling becomes important near the walls where the viscous effects dominate. Here, Td is the time required for the particles to be independent of their initial locations. It should be noted that no explicit value of the development time for glass particles was given by Wang and Squires and that the current value was obtained from Fukagata et al. [43] who performed a LES calculation for a particle-laden turbulent channel flow at Re=180. They used a square root relation between the particle response time and the development time (Td).
30
Table 4-1: Channel properties used in the calculations
u
0.294 m/s 0.009 m 1.47 X 10-5 m2/s 1.02145
g+
Table 4-2: Particle properties for a turbulent channel flow at Re=180 Lycopodium
d p in m
p in kg/m
Td in sec
3
Glass
50 2500 2.5 / u 0.650 3.611 x 10-3
Copper
70 8800 6 / u 4.500 2.5 x 10-2
28 610 0.5 / u 0.048 2.67 x 10-4
( Stp ) , u ( Stp ), u
31
4.3 Simulation overview
The Eulerian flow field was time advanced to reach a statistically stationary state. At this point 250,000 particles were then randomly distributed throughout the channel with an initial velocity same as the fluid velocity at the particle location. Wang and Squires [2] showed that this sample size gave a smooth representation of the particle statistics. The particle position was updated by integrating the following equation:
du ip + 1 = - + 1 + 0.15 Re 0.687 (u ip + - vif + ) + g + p + dt p
(
)
(4.1)
+ where g =
g u2
, and p = particle response time.
The first term on the right hand side of equation (4.1) is the Stokes drag and the second term is the gravitational force acting in the streamwise direction. A detailed description of equation (4.1) is provided in chapter 2.
Elastic collisions were assumed for the particles contacting the wall. A particle was assumed to contact the wall when its center was one radius away from the wall. If a particle contacts the wall a with normal velocity U 2p , its normal velocity after impact is given by U 2p = -U 2p (elastic collision). The other components of the velocity remain unchanged. A description of other wall-collision models used in particle-laden flows is provided in Appendix D. Periodic conditions were applied on the particles leaving from the stream-wise and span-wise directions. These particles were re-introduced in the domain at the corresponding location with the same exit velocity. Before calculating any statistics, particles were time advanced till they became independent of their initial
32
locations. Following Wang and Squires [2], copper particles were integrated for 6/u (6 non-dimensional time units) in order to allow the particles to adjust to the flow. In comparison, glass particles were integrated for 2.5 /u (2.5 non-dimensional time units), while lycopodium particles were integrated for 0.5 /u (0.5 non-dimensional time unit). After the development time, particle statistics were accumulated at every one dimensionless time unit for a total time interval of 6 non-dimensional units, i.e., 6 < tu/ < 12 for copper, 2.5 < tu/ < 8.5 for glass, and 0.5 < tu/ < 6.5 for lycopodium. Statistics were calculated by dividing the channel into slabs parallel to the wall and then averaging over both, the channel halves and time.
4.4 RESULTS AND DISCUSSION
4.4.1 Mean streamwise velocity
Figure 4.4 shows a comparison of the mean streamwise velocity predicted by the current model with earlier simulation by Wang and Squires [2]. As may be observed, good agreement was found for all the particle sizes. It can be seen that the inertial copper particles which respond slowly to the turbulent structures have a higher velocity than the fluid throughout the channel. Wang and Squires [2] noted that the "small plateau" in the near wall region is due to the fact that the copper particles transported from the outer regions of the channel carry high velocity and retain a significant fraction of their streamwise momentum after colliding with the walls, which was also observed by Kulick
et al. [16]. Lycopodium particles which respond well to the changes in flow field, act like
33
tracer particles and closely follow the fluid. Between these two particle sizes are the glass particles that respond to some turbulent structures. It can be observed that they exhibit a high velocity near the wall, but closely follow the fluid after a y+ > 10. Figure 4.5 (a)-(c) shows the mean streamwise velocity, respectively, for lycopodium, glass, and copper particles at different time intervals. A constant decrease in the streamwise velocity is observed in the near wall region. It can be seen that the decrease is more for the heavier particles compared to the lighter particles.
4.4.2 RMS velocity fluctuation
Figure 4.6 (a)-(c) compares the velocity fluctuations predicted by the current model with those calculated by Wang and Squires [2]. As may be observed a good agreement is obtained for all the particle sizes. It is clear from the figures that with the increase in Stokes number the streamwise velocity fluctuations increase while the normal and spanwise velocity fluctuations decrease. Compared to the fluid velocity fluctuations in Figure 4.3 it was observed that near the wall the streamwise velocity fluctuations of copper particles are stronger than those of the fluid while, near the centerline; values are similar to the fluid.
34
4.4.3 Number density profiles
Figure 4.7 (a)-(c) show the time development of the number density profiles, respectively, for lycopodium, glass, and copper particles. calculated as follows: The number density is
NDENSITY =
number of particles in a computational cell / cell volume total number of particles / total volume of the domain
Lycopodium particles: It can be seen that in the near wall region (0<y+<4), lycopodium
particles exhibit a steep increase in the number density profile, i.e., from 5 at 2 < tu/ < 3 to 12 at 11 < tu/ < 12. This is a result of the accumulation of these particles near the walls. It should be noted that the increase in the number density profile with time decreases with the passage of time. Comparing the number density profiles accumulated for each non-dimensional unit from 2 < tu/ < 3 to 11 < tu/ < 12; in the region y+ > 5, one finds that the number density remains almost constant.
Glass particles: Figure 4.7 (b) depicts the evolution of particle number density profiles
for glass particles. Similar behavior as the case of Lycopodium particles was observed. The only difference being the rate increase in the number density profiles with time. As compared to a decreasing rate for the lycopodium particles, the rate of increase remains almost steady for the glass particles. As noted for the lycopodium particles, the increase in the number density profiles is restricted to a thin region near the wall. A slight decrease is observed in the region (4 < y+ < 20) followed by a constant profile for the region y+ > 25.
35
Copper particles: In comparison to the glass and lycopodium particles, a mixed
behavior is observed for the heavy copper particles. An investigation of Figure 4.7 (c) reveals that a constant increase in the number density profile is observed for the region near the wall, i.e., 0<y+<4, for the entire period over which the statistics were accumulated. Comparing the number density profile (for the near wall region, i.e. 0<y+<4) accumulated at time 6 < tu/ < 7 with that at 11 < tu/ < 12, it is found that lycopodium and glass exhibit a high accumulation (number density increasing from 5 to 15 for lycopodium and number density increasing from 7 to 16 for glass) while the copper particles exhibit a relatively low accumulation (number density increasing from 2.2 to 2.6). Considering the region away from the wall, it can be seen that the number
density decreases in the region 10 < y+ < 90 followed by an increase in the region y+ > 90 as compared to a constant number density profile for the lycopodium and glass particles.
4.4.4 Number of particles colliding with the wall
Figure 4.8 (a) (c) illustrates the time development of the percentage of the total number of particles colliding with the wall during the period over which statistics are accumulated. The percentage is calculated as follows:
NHITS = number of particles colliding with the wall X 100 total number of particles (= 250000)
A comparison of the plots reveals that initially the fraction of particles colliding with the wall is more for the particles with higher Stokes number. This indicates that the particles are still adjusting to the flow. Comparing the profiles after the development time of
36
copper, i.e., t+ = 6, it can be observed that the copper particles achieve a steady state as compared to an increasing percentage of the lycopodium and glass particles. Close analyses of the plots reveal that for the lycopodium and glass particles, this fraction is almost steady after time, t+ = 11. Further, it can be seen that after the development time, the fraction of particles colliding with the wall is higher for glass particles (8-10%) as compared to copper (4-6%) and lycopodium (2-5%) particles.
4.4.5 Preferential concentration
Preferential concentrations are known to exist in the near wall region of a turbulent channel. Following conclusions are reached from past studies on turbulent channel flows:
Eaton and Fessler [18] have shown that inertial bias in the particle trajectories causes preferential concentration of particles in regions of low vorticity or high strain rate.
Wang and Maxey [19] have shown that particles with Stokes number closer to Kolmogorov scales exhibit stronger effects of preferential particles.
Rouson and Eaton [20], and Pedinotti et al. [17] have shown that particles tend to accumulate in low speed streaks near the wall.
In order to observe the above findings, snapshots of instantaneous particle distributions were taken at t+ = 6 for all the particle sizes as shown in Figure 4.9 (a)-(c). All the particles within a y+ < 5 are plotted. As expected, the lycopodium particles exhibit
37
stronger preferential concentration than the glass and copper particles. Comparing the streaky structures of the particles, it can be concluded that the lighter lycopodium particles accumulate in the low speed streaks as compared to less organized distribution of the glass particles and a random (or uniform) distribution of the heavy copper particles. It is evident from the distributions that with the increase in Stokes number, the homogeneity in particle distribution increases. Recent investigations by Fessler et al. [21] have shown that preferential concentration exists throughout the channel. This can be corroborated from Figure 4.10 (a)-(c) which depicts the instantaneous distribution of the particles at the channel centerline, i.e., 170 < y+ < 180. Similar to the near wall region, the copper particles are randomly distributed whereas the lycopodium and glass particles exhibit varying degrees of preferential concentration. Similar behavior was observed by Wang and Squires [2] for a turbulent channel flow at Re= 644.
38
Figure 4.1 Computational domain non-dimensionalized by channel half width ()
16
15.8
mean velocity
15.6
15.4
15.2
15 140
160
non-dimensional time
180
200
220
240
260
Figure 4.2: Time evolution of the mean flow velocity on the top and bottom wall
2.6 2.4 2.2
symbols DNS (KMM) present urms
0.8
urms.vrms,wrms
1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 10
0
turbulent shear stress
2
0.6
0.4
wrms
symbols
DNS (KMM) present
vrms
0.2
y+
10
1
10
2
0
10
-1
10
0
y+
10
1
10
2
(a).
(b).
Figure 4.3: (a) RMS Fluctuations of normal turbulent stresses; (b) turbulent shear stress.
39
Figure 4.4: Comparison of mean streamwise velocity predicted by current model with particle statistics by [2]. Particles compared: 70m copper, 50m glass, and 28m lycopodium
40
(a).
(b).
(c). Figure 4.5: Mean stream-wise velocity (a). Lycopodium (b). Glass (c). Copper
41
(a).
(b).
(c). Figure 4.6: Comparison of velocity fluctuations for (a) Lycopodium (b). Glass (c). Copper
42
(a).
(b).
(c). Figure 4.7: Number density profile (a). Lycopodium (b). Glass (c). Copper
43
(a).
(b).
(c). Figure 4.8: Number of particles colliding with the wall as a % of the total number of particles (a). Lycopodium (b). Glass (c). Copper
44
(a).
(a).
(b).
(b).
(c).
(c).
Figure 4.9: Instantaneous particle distribution at t+ = 6, 0 < y+ <5, Re = 180 (a). Lycopodium (b). Glass (c). Copper
Figure 4.10: Instantaneous particle distribution at t+ = 12, 170 < y+ <180, Re = 180 (a). Lycopodium (b). Glass (c). Copper
45
CHAPTER 5. INTERNAL COOLING RIBBED DUCT
5.1 Cooling of turbine blades:
Gas turbines used for propulsion and power generation are subjected to high temperatures (1200-1500 oC) to improve the thermal efficiency and power output. With the increase in turbine inlet temperature the heat transferred to the turbine blade increases. Such high operating temperatures are far above the permissible limit for the blade material. A variation in the temperature within the blade material produces thermal stresses which directly impacts the blade life and hence the durability of the gas turbine. Hence for safe operation of the turbine it is necessary to cool the turbine blade. The cooling of a turbine blade can be divided in three zones as shown in Figure 5.1, i.e., leading edge cooling, cooling of the pressure and suction surfaces, and trailing edge cooling. In modern gas turbines, the required cooling air is extracted from individual compressor stages. This air enters from the base of the blade and flows through a multipass duct or serpentine passage to provide internal convective cooling. In order to enhance heat transfer, the surface of these serpentine passages is roughened by turbulence promoters. Turbulence promoters can take varying shapes, sizes and configurations, ranging from ribs to pin fins, the former being commonly used. The cooling air after flowing through the internal passages is directed to the film cooling holes to form a protective layer on the surface of the turbine blade. The protective layer prevents direct contact between the hot gases and the turbine blade material. The leading edge is cooled by jet impingement with film cooling while; the trailing edge is cooled by pin fins as shown in Figure 5.1.
46
Figure 5.1: Schematic of the cooling of a modern gas turbine blade [51]
47
Gas turbines operating in deserts ingest sand with the intake air. Previous studies have shown that this air which is responsible for providing internal convective cooling is capable of eroding the compressor blades to a limit that the gas turbine cannot provide adequate power. The objective of this chapter is to apply the developed model to understand the particle transport in an internal cooling duct which is downstream of the compressor. Of particular interest are the following: 1. How do particles of different sizes (10m, 50m, and 100m) impinge on the duct walls and ribs? What is the average impingement velocity and angle of the particles hitting the walls? 2. Which are the areas prone to erosion due to prolonged operation of the gas turbine in deserts and other dusty environments?
5.2 Computational details
Large-Eddy simulations at Re=6660 (nominal Reb =
u b Dh
=20,000), based on
friction velocity (uc = u ) and hydraulic diameter (Lc = Dh) are carried out on a grid consisting of 128 x 128 x 128 cells in the streamwise, wall-normal, and spanwise directions, respectively [25]. It is assumed that the flow is fully-developed and only one rib pitch is included in the computational domain. The computational domain normalized by the hydraulic diameter is shown in Figure 5.2. A square rib (e/Dh=0.1 and P/e=10) is placed on two sides of the duct, where e is the rib height and P the streamwise pitch. High mesh density was used in the vicinity of the rib and the duct surfaces to resolve the
48
turbulent boundary/shear layers, which is crucial to the accurate prediction of turbulence and particle dynamics. The grid spacing in the vicinity of walls was such that the first mesh point was within + < 1 with 4 5 mesh points within 10 wall units. In the streamwise and spanwise directions the grid spacing varied between 5 and 30 wall units. LES by Tafti [25], established that the turbulent flow and heat transfer are predicted within experimental uncertainty. The non-dimensional time step used was 5 x 10-5. No slip boundary condition is applied on the walls in the normal and spanwise directions. Periodic boundary condition is assumed in the streamwise direction. Table 5-1 - Table 5-3 describe the flow and geometrical parameters and the particle properties used in the calculations. Under these conditions, the three particle sizes exhibit Stokes numbers ranging from 0.07 to 6.9.
Table 5-1: Flow and geometrical parameters used in the calculation:
u ub Dh
56.05 m/s 168.15 m/s 0.01 m 8.49 X 10-5 m2/s
Table 5-2: Properties of the carrier phase CARRIER PHASE - AIR Temperature 800 K Pressure 20 atm Density 8.8 kg/m3
49
Table 5-3: Particle properties used in the calculation: DISPERSED PHASE - SAND 10 50 1650 1650 p in kg / m 3
d p in m
100 1650 6.87500 20.62500
p u / Dh p u b / Dh
0.06875 0.20625
1.71875 5.15625
5.3 Choice of particle sizes
Van der Walt and Nurick [44] conducted experiments on a dust filtered helicopter turbine engine and concluded that only particles with diameter less than 100m are ingested into the engine while the others were filtered out. In addition, they observed that all particles below 100m can be considered sparse and that inter-particle collisions can be neglected. This again reconfirms our previous assumption of one-way coupling. Also, Jensen et al. [5] in their study of deposition in land based turbine observed that particle sizes between 0-80 microns are ingested into the engine. Although their objective was to study deposition in a land based turbines we believe that it resembles the situation of take off and landing. and 100m. Based on this data, three particle sizes are studied: 10m, 50m,
50
5.4 Simulation Overview
Once a statistical stationary stage for the carrier phase is established, 268,800 particles are distributed uniformly in the computational domain. The particles are allowed to evolve for 23.5 non-dimensional time units which corresponds to between 340 particle time constants for the 10 m particles, 13.7 time constants for 50 m particles, and 3.42 time constants for 100 m particles. During integration, the number of particles in the calculation domain remains constant. Particles leaving the domain at the periodic boundary are reintroduced into the domain at the exit velocity at the equivalent location.
A particle is said to impinge at a given location when its center is one radius away from the wall. Elastic collisions are assumed, although in reality some particle energy will be imparted to the wall.
When a particle impinges at a given location, its impingement velocity, which is the magnitude of its velocity vector and the angle it makes with the surface is calculated. Particle statistics are accumulated for 7 (23.5 to 30.5) non-dimensional time units.
At the end of the simulation, an average impact velocity and angle is calculated at all impingement locations on the surfaces. Although the whole domain is considered for the calculation, owing to the symmetry along the y and z directions, the average particle statistics for a quadrant are presented. This quadruples the statistical sample size. All calculations utilize 32 processors of the 1100 compute nodes (dual 2.3 Ghz Xserve G5) available at the Virginia Tech supercomputer (System X). Each time step with the dispersed phase consumes about 12-16 hours of wall clock time.
51
5.5 Results and Discussion
Before describing the results, it is important to clarify that the number of particles impinging on the walls of the duct are reported for a finite time, which in this case is 7 non-dimensional time units or 1.3 milliseconds. To put this time in perspective, a typical gas turbine service cycle is 10,000-25,000 hours.
5.5.1Mean Flow Features
Figure 5.3(a) shows the mean streamline pattern at the center of the duct (z=0.5). The mean flow is characterized by a leading edge eddy at the rib-wall junction, the counter-rotating eddy in the rib wake, the main recirculation region and the recirculation region on top of the rib. In the vicinity of the smooth walls (z=0.0) the flow field becomes strongly three-dimensional with mean cross flow velocities (wb) approaching 30% of the mean streamwise velocity (ub). Figure 5.3 (b) shows contours of wb in a plane z=0.05 in the vicinity of the smooth wall. Strong localized cross flow components are found to occur. Of particular interest is the high lateral velocity moving towards and impinging on the smooth wall within the confines of the shear layer at the leading edge of the rib.
5.5.2 Instantaneous Particle Distributions
Figure 5.4(a)-(c) shows the instantaneous distribution of sand particles, respectively, for 10, 50, and 100 micron sizes. These snapshots were taken after
52
advancing the particles for 27 dimensionless time units from their initial random locations. Any particle within a distance of 0.05Dh from the wall or rib is plotted. Flow direction is shown by an arrow in the plot. The following observations can be made:
Calculating the fraction of particles within a distance of 0.05Dh from the walls and ribs, it was found that 40% of the 10 micron particles fall in that region as compared to 26% and 27%, respectively, for 50 and 100 micron particles.
The smaller 10 micron particles which closely follow the fluid are seen to accumulate at the center of the side-walls. In comparison, the 50 and 100 micron particles are homogeneously distributed. This can be attributed to the fact that particles with Stokes number closer to Kolmogorov scales exhibit stronger effects of preferential concentration [45]. In addition, previous studies [46]-[47] in wall-bounded flows show that particles tend to accumulate in regions of low streamwise velocity. This can be corroborated by considering the instantaneous streamwise velocity at a plane 0.05 from the side-wall, shown in Figure 5.5. It can be observed that in the center region of the side wall where high accumulation (relatively darker regions) is seen, instantaneous stream wise velocity is between 0-10 % of the mean bulk velocity.
Another clear trend with particle size is the distribution in the vicinity of the ribbed wall. The number of 10 m particles near the ribbed wall is much less than that present for the larger particle sizes. This can be attributed to the presence of the separated shear layer on the rib and the resulting large separation zone behind the rib (Figure 5.3 a). For the 10 m particles, very few particles are entrained into the separation zone because they are prone to follow the fluid streamlines. On the other
53
hand, the larger particles because of their inertia infiltrate into the separation zone and exhibit a larger concentration near the ribbed wall.
5.5.3 Characteristics of Time-Averaged Particle Impingement
In order to calculate the material removed from a surface, empirical correlations relating the material properties of the surface with the impact velocity and angles are needed. The purpose of our study is to map the impact velocity and angles for the given particle sizes, which when combined with the empirical correlations can give a measure of the amount of material removed from the surface during a given operation time of the gas turbine. Three quantities are presented for each surface:
The total number of particles impinging in 1.3 milliseconds out of a total population of Ntot = 268,800 particles.
The average magnitude of impingement velocity. The average impingement angle with respect to the surface. Finally, to estimate the potential for erosion or deposition the three quantities are
combined to define the fraction of energy of incoming particulates which is potentially transferred to the surface of the duct. In order to estimate this quantity, an equivalent particle mass flow rate is calculated and is given by:
m p = N tot m p / (P / ub )
(
)
(5.1)
which assumes that the mean particle velocity is the same as the bulk flow velocity. Hence the total flow of particulate energy per unit time is given by:
2 E ptot = 1 / 2m pub
(5.2)
54
In a given time period, a measure of the rate of energy which can be potentially transferred to the surface due to particle impingement is given by:
E pimpinge = 1 / 2 Nimpinge m p / t u pimpinge sin 2
(
)(
)
(5.3)
where Nimpinge is the number of particles that impinge at a given location, upimpinge and is the impingement velocity magnitude and angle, respectively. The potential for erosion and deposition is defined by the fraction
= E pimpinge / E ptot
(5.4)
It is noted that is a fraction of incoming energy and although may be smaller for 50 and 100 microns, the total incoming energy is larger for the larger particles.
5.5.4 Smooth Side wall
Flow near the smooth wall is characterized by a streamwise velocity and a spanwise cross-flow which varies in sign and magnitude as shown in Figure 5.3-b. Figure 5.6 (a)-(c) shows the number of particles colliding with the smooth wall in 7 dimensionless time units or 1.3 milliseconds for half the duct. It can be observed that for each particle size, impingement is highest above and on either side of the rib1. A careful investigation of the plots reveal that in the region above the rib within the confines of the shear layer where the flow accelerates, the number of particles colliding with the smooth wall decrease with an increase in the particle Stokes number. Also, in the region behind the rib in its wake, where outward fluid motion occurs (Figure 5.3 b), the number of
1
The extended horizontal and vertical regions of low impingement density at the rib edges is a result of the fine grid density in that region and the lower probability of a particle being located in those cells.
55
particles impinging on the smooth wall increases with an increase in the particle response time. A possible explanation for this phenomenon can be attributed to the ability or inability of a particle to respond to the fluid motion. The small 10 micron particles which closely follow the fluid respond well to this outward motions and hence cause little impingement as compared to the inertial 50 and 100 micron particles. These particles which are governed by their inertia and respond slowly may continue on their trajectories and impinge on the smooth wall.
Figure 5.7 (a)-(c) shows the impingement velocity, respectively, for 10, 50, and 100 micron particles. Following observations are made:
On the smooth wall region at the center of the duct, the 10 micron particles impinge at varying velocities as compared to a uniform impingement velocity of the 50, and 100 micron particles. The higher impact velocities of the 50 and 100 micron particles indicate that they carry more momentum than the smaller particles.
In the shear layer above the rib, the velocity distribution of the 10 micron particles correlates with the secondary flow in Figure 5.3 b, whereas the larger particles have no preferential impingement velocities in this region. In the rib wake, consistent with the observation in Figure 5.6 (a-c), impingement velocities increase steadily as the particle size increases because the smaller particles are forced away from the wall by the secondary flow velocities.
Particles impinge with a uniform impact velocity (20-30 % of the mean bulk velocity) on the smooth wall in front of the rib (0 < y/Dh < 0.1).
56
The average angles at which these particles impact the smooth wall are shown in Figure 5.8 (a)-(c). It can be observed that in the high impingement region particles have a very shallow angle indicating that streamwise motion of the fluid dominates the particle trajectories and that they move nearly parallel to the smooth wall. In comparison, particles impinge at an angle between 10-25o on the smooth-wall in front of the rib (0<y/Dh<0.1). This indicates that particles get trapped in the instantaneous secondary vortices produced in front of the rib and at some instant impinge on the smooth wall. Figure 5.9 (a-c) plots the fraction 106 which consolidates the data in Figure 5.6 Figure 5.8 gives a measure of the potential for erosion and deposition at the side wall of the duct. For the 10 micron particles, has the largest values in the region of secondary flow impingement in the vicinity of the rib. As the particles get bigger, has a much broader coverage and spreads to the center of the duct. Interestingly, the 50 micron particles have the broadest coverage and the largest potential for erosion.
5.4.5 Ribbed Wall
Figure 5.10 (a)-(c) shows the number of particles colliding with the ribbed wall. It can be seen that upstream of the rib, the 50 and 100 micron particles exhibit uniform impingement as compared to selective impingement by the 10 micron particles, which exhibit the highest impingement at the base of the rib, indicating that the smaller particles are entrained into the instantaneous vortical structures which form in this region. It is observed that the 50 micron particles cause more impingement than the 100 micron particles upstream of the rib. In comparison, downstream of the rib, particle impingement
57
increases with an increase in particle Stokes number. Considering the region in the immediate vicinity of the rib, it can be seen that high impingement is caused by the 100 micron particles. Figure 5.11(a)-(c) shows the average impingement velocity of the particles on the ribbed wall. It can be observed that 10 and 50 micron particles impinge with low velocities downstream of the rib as compared to upstream of the rib. In comparison, the 100 micron particles impinge at higher velocities both in the front and back of the rib. Figure 5.12 (a)-(c) shows the angles at which these particles impinge on the wall. The low impingement angles indicate that particles impact tangentially to the wall. Considering the high impingement angles in the immediate vicinity of the rib, it appears that particles get entrained in the unsteady secondary vortices (Figure 5.3-a) and impinge on the bottom wall. The potential for erosion occurring at the ribbed wall is of the same order of magnitude as the side wall. Figure 5.13 (a-c) plots for the three particle sizes. There is a gradual increase in the magnitude of as the particle size increases. While 50 micron particles have a uniform broad coverage, the 100 micron particles show much higher values downstream of the rib in the reattachment region.
5.4.6 Ribs
Figure 5.14 (a)-(c) shows the number of particles impinging on the front surface of the rib. From all the other surfaces, the front of the rib has by far the largest amount of impingement. The impingement density increases with an increase in particle size. It is interesting to note that the impingement angle steadily increases with particle size in Figure 5.16 (a)-(c) as well. This indicates that the 10 micron particles react to the change
58
in flow direction much quicker than the larger particles, which due to their inertia impinge more directly with larger angles. However, even for the 100 micron particles, the impingement angle does not exceed 50 degrees. It is also observed that the impingement angles exhibit alternating high-low values. Interestingly this pattern is also found at the back of the rib and only in the impingement angles. This may be a consequence of localized vortical structures which exist in these regions and the finite sampling time. As with the angles, there is also a steady increase in impingement velocity with particle size as shown in Figure 5.15 (a)-(c). Overall, the front of the rib is most prone to erosion as indicated by the distribution of shown in Figure 5.17 (a-c). The magnitude of is approximately one to two orders of magnitude larger than that on the ribbed and side wall. It is observed that the 100 micron particles are the most damaging.
In contrast the top of the rib undergoes the least amount of damage as shown in Figure 5.18 (a-c) is an order of magnitude smaller than that experienced on the walls of the duct. Surprisingly, the back of the rib is not exempt from damage as one would expect. Particle impingement on the back of the rib is shown in Figure 5.19 (a)-(c). The flow in this region is characterized by a primary separation zone as shown in Figure 5.3 (a). Hence only particles which are entrained into the recirculation zone are positioned to impinge on the back of the rib. It is observed that particle impingement increases with the increase in particle Stokes number. As pointed out earlier, the 10 micron particles for the most part follow the flow streamlines (see Figure 5.10a-c) over the rib and do not get entrained into the separation region behind the rib. Considering the impingement angles in Figure 5.21(a)-(c), it can be seen that the 50 micron particles which are somewhat
59
more "agile" than the 100 micron particles and are able to change direction in the recirculating flow impinge at higher angles. In comparison, the larger 100 micron particles after infiltrating in the recirculation zone follow their trajectories and impinge on the back surface of the rib at lower angles. While for all particle sizes, impingement angles are in the range from 15 to 50 degrees, there is a substantial difference in the impingement velocity between the 10-50 micron particles and the 100 micron particles. Figure 5.20 (a-c) shows that the 10 and 50 micron particles impinge at a low velocity (2030% of the mean bulk velocity) as compared to the 100 microns (40-50% of the mean bulk velocity). Figure 5.22 (a-c) shows the potential for damage at the back of the rib. It is observed that the fraction of kinetic energy potentially transferred to the back of the rib is highest for the 100 micron particles and is comparable to that at the front of the rib.
60
1
Y
0.5
0 0 0.2 0.4 0.8 0.6 0.4 0.8 1 0 0.2
1
Flow
Z
0.6
X
Figure 5.2: Computational domain non-dimensionalized by hydraulic diameter [25].
61
0.6
128 , DSM
3
center plane, z=0.5
0.5
0.4
y
0.3
0.2
0.1
0
0
0.1
0.2
0.3
0.4
0.5
x
0.6
0.7
0.8
0.9
1
(a)
0.6
3
128 , DSM plane, z=0.05
0.05
0.5
0.4
y
-0.01
0.3
0.2
7
0.1
-0.19 -0.31
-0. 0
-0.1
-0.04
0
0
0.1
0.2
0.3
0.4
0.5
x
0.6
0.7
0.11
0.0 5
0.8 0.9
0.02
1
(b) Figure 5.3:
(a) Mean
streamline distribution in the z-symmetry plane. Reattachment
occurs at 4.1e downstream of rib. The leading edge eddy extends between 0.8-0.9e upstream of rib. (b)Mean lateral or spanwise flow velocity (wb) in the vicinity of the smooth wall [25].
62
(a).
(b).
(c). Figure 5.4: Instantaneous distribution of sand at t = 27 Dh / u a) 10m, b) 50m, c) 100m.
Figure 5.5: Instantaneous streamwise velocity at a plane z / Dh = 0.05
63
(a).
(b).
(c).
Figure 5.6: Number of particles impinging the side wall in 1.3 milliseconds. (a) 10m, (b) 50m, (c) 100m.
(a)
(b)
(c)
Figure 5.7: Impact velocity on the side wall. (a). 10m (b). 50m (c). 100m
(a)
(b)
(c)
Figure 5.8: Impact angle on the side wall. (a). 10m (b). 50m (c). 100m
64
(a)
(b)
(c)
Figure 5.9: Measure of potential for erosion and deposition at side wall ( 106). (a) 10m, (b) 50m, (c)100m
(a)
(b)
(c)
Figure 5.10: Number of particles impinging the ribbed (bottom) wall in 1.3 milliseconds. (a) 10m, (b) 50m, (c) 100m.
65
(a)
(b)
(c)
Figure 5.11: Impact velocity on the ribbed wall. (a). 10m (b). 50m (c). 100m
(a)
(b)
(c)
Figure 5.12: Impact angle on the ribbed wall. (a). 10m (b). 50m (c). 100m
66
(a).
(b).
(c).
Figure 5.13: Measure of potential for erosion and deposition at bottom wall ( 106). (a) 10m, (b) 50m, (c)100m
(a)
(b)
(c) Figure 5.14: Number of particles impinging the front surface of the rib in 1.3 milliseconds. (a) 10m, (b) 50m, (c) 100m.
67
(a)
(b)
(c) Figure 5.15: Impact velocity on the front surface of the rib. (a). 10m (b). 50m (c). 100m
(a)
(b)
(c) Figure 5.16: Impact angle on the front surface of the rib (a). 10m (b). 50m (c). 100m
68
(a)
(b)
(c) Figure 5.17: Measure of potential for erosion and deposition on front surface of the rib ( 106). (a) 10m, (b) 50m, (c)100m
(a)
(b)
(c) Figure 5.18: Measure of potential for erosion and deposition on top surface of the rib ( 106). (a) 10m, (b) 50m, (c)100m
69
(a)
(b)
(c) Figure 5.19: Number of particles impinging the back surface of the rib in 1.3 milliseconds. (a) 10m, (b) 50m, (c) 100m.
(a)
(b)
(c) Figure 5.20: Impact velocity on the back surface of the rib (a). 10m (b). 50m (c). 100m
70
(a)
(b)
(c) Figure 5.21: Impact angle on the back surface of the rib (a). 10m (b). 50m (c). 100m
(a)
(b)
(c) Figure 5.22: Measure of potential for erosion and deposition on back surface of the rib ( 106). (a) 10m, (b) 50m, (c)100m
71
CHAPTER 6. SUMMARY AND CONCLUSIONS
A two-phase dispersed flow capability was developed for studying particle dynamics in complex turbulent flow. A parallelized 3-D particle tracking code was developed to solve the dispersed phase in a Lagrangian framework. The code was integrated in to an existing multi-block Navier Stokes solver (GenIDLEST [26]) to carry out large scale turbulent particle-laden flows on parallel computers. The developed algorithm was validated by comparing the computed particle statistics for 28 microns lycopodium, 50 microns glass, and 70 microns copper with available data [2] for a turbulent channel flow at Re=180. The particle statistics comprised of particle velocity fluctuations, instantaneous distribution near the wall and at the centre of the channel, and time development of streamwise velocity, number density profiles, and percentage of the total number of particles colliding with the walls. Statistics were calculated by dividing the channel into slabs parallel to the wall and then averaging over both, the channel halves and time. Good agreement was found for all the computed statistics.
LES calculations were performed for a bulk Reynolds number of 20,000 to study particle transport in an internal cooling duct with normal ribs. Particle sizes of 10, 50, and 100 microns with non-dimensional response times, respectively, 0.06875, 1.71875, and 6.875 were considered. The particles were time advanced for 23.5 dimensionless time units to become independent of their initial locations. Particle statistics were gathered for 7 dimensionless time units or 1.3 milliseconds. Although the whole domain was considered for the calculation, owing to the symmetry along the y and z directions, the average particle statistics for a quadrant were presented. Particle statistics, comprising of
72
number of particles impinging at a given location, and their average impact velocity and angle are presented for each surface. It is found that at any given instant in time about 40% of the total number of 10 micron particles are concentrated in the vicinity (within 0.05 Dh) of the duct surfaces, compared to 26% of the 50 and 100 micron particles. The 10 micron particles are more prone to preferential concentration than the larger particles. At the side walls, while the 10 micron particles exhibit a high potential to erode the region in the vicinity of the rib, the larger particles are more prone to eroding the area between the ribs and towards the center of the duct. At the ribbed walls, while the 10 micron particles exhibit a fairly uniform but low propensity for erosion, the 100 micron particles show a much higher tendency to erode the surface in the vicinity of the reattachment region. The rib face facing the flow is by far the most susceptible to erosion and deposition for all particle sizes. While the top of the rib does not exhibit a large propensity to be eroded, the back of the rib is as susceptible as the front surface of the rib because of particles which are entrained into the recirculation zone behind the rib.
We have generated useful new data on particle transport in a square ribbed duct using large-scale time-accurate numerical simulations. The results reported in this thesis are valuable information to study deposition and erosion in a ribbed duct. The computed data when combined with empirical relations correlating material properties and particle impingement velocity and angle can give a measure of the amount of material removed from the surface during a given operation time of the gas turbine. To the author's knowledge, these results are the first that use Large Eddy Simulation coupled with Lagrangian particle tracking for particle-laden flow in an internal cooling ribbed duct.
73
References:
1. Germano, M., Piomelli, U., Moin, P. and Cabot, W. H., "A dynamic subgrid-scale eddy viscosity model", Phys. Fluids A, vol. 3, pp. 1760-1765, 1991. 2. Wang, Q., Squires, K.D., "Large eddy simulation of particle-laden turbulent channel flow", Phys. Fluids 8 (5), May 1996. 3. Kleinstreuer, C., "Two-Phase Flow: Theory and Applications." Taylor and Francis, New York, NY, 2003. 4. Sharma, G., "Direct Numerical Simulation of particle-laden turbulence in a straight square duct", Masters Thesis, Mechanical Engineering, Texas A&M University, 2004. 5. Jensen, J.W., Squire, S.W., Bons, J.P., and Fletcher, T.H., "Simulated land-based turbine deposits generated in an accelerated deposition facility", International Gas Turbine and Aeroengine Congress and Exhibition, Vienna, Austria, June 13-17, 2004, GT2004-53324. 6. Ghenaiet, A., Elder, R.L., and Tan, S.C., "Particles and Trajectories through an axial fan and performance degradation due to sand ingestion," ASME, 2001-G-497, 2001. 7. 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. Eng. Gas Turbines Power, 115, 641-651. 8. Bons, J.P., "St and Cf augmentation for real turbine roughness with elevated freestream turbulence," ASME J. Turbomachinery, 124, 632-644. 9. Scala, S.M., Konrad, M., and Mason, R.B., "Predicting the performance of a gas turbine undergoing compressor blade erosion", 39th AIAA/ASME/SAE/ASEE Joint
74
Propulsion Conference and Exhibit, Huntsville, Alabama, July 20-23, 2003, AIAA 2003-5259. 10. Hathaway, M.D., Chen, J., and Webster, R., "Time accurate unsteady simulation of the stall inception process in the compression system of a US Army helicopter gas turbine engine", DoD HPCMP User's Group Conference (UGC), Bellevuew, Washington, June 9-13, 2003. 11. Tabakoff, W., and Hamed, A., "Installed engine performance in dust-laden atmosphere", AIAA/AHS/ASEE Aircraft Design Systems and Operations Meeting, San Diego, CA, October 31 November 2, 1984, AIAA-84-2488. 12. Tabakoff, W., Lakshminarasimha, A. N., and Pasin, M., "Simulation of compressor performance deterioration due to erosion " , 34th International Gas Turbine and Aeroengine Congress and Exhibition, Toronto, Canada, June 4-8, 1989. 13. Tabakoff, W., "Compressor erosion and performance deterioration", AIAA/ASME 4th Joint Fluid Mechanics, Plasma Dynamics and Laser Conference, Atlanta, GA, May 12-14, 1986. 14. Balan, C., Tabakoff, W., "A method of predicting the performance deterioration of a compressor cascade due to sand erosion", AIAA 21st Aerospace Sciences Meeting, Reno, Nevada, January 10-13, 1983, AIAA-83-0178. 15. Schneider, O., Dohmen, H.J., Benra, F.-K., Brillert, D., "Investigations of dust separation in the internal cooling air system of gas turbines", Proceedings of ASME Turbo Expo 2003, Atlanta, GA, June 16-19, 2003, GT2003-38293. 16. Kulick, J.D., Fessler, J.R., and Eaton, J.K., ``Particle response and turbulence modification in fully turbulent channel flow'', J. Fluid Mech. 277,109, 1994.
75
17. Pedinotti, S., Mariotti, G., and Banerjee, S., "Direct Numerical Simulation of particle behaviour in the wall region of turbulent flows in horizontal channels", Int. J. Multiphase Flow, 18, 6. 927-941, 1992. 18. Eaton, J.K., and Fessler, J.R., ``Preferential concentration of particles by turbulence", Int. J. Multiphase Flow, 20, 169, 1994. 19. Wang, L.-P., and Maxey, M.R., ``Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence,'' J. Fluid Mech., 256, 27, 1993. 20. Rouson, D.W.I, and Eaton, J.K., ``Direct numerical simulation of particles interacting with a turbulent channel flow,'' Proceedings of the 7th Workshop on Two-Phase Flow Predictions, edited by M. Sommerfeld, Erlangen, Germany, 1994. 21. Fessler, J.R., Kulick, J.D., and Eaton, J.K., ``Preferential concentration of heavy particles in a turbulent channel flow,'' Phys. Fluids, 6, 3742, 1994. 22. Marchioli, C., and Soldati, A., "Mechanisms for particle transfer and segregation in a turbulent boundary layer", J. Fluid Mech, 468, 283-315, 2002. 23. Madabhushi, R.K., and Vanka, S.P., "Large Eddy Simulation of turbulence-driven secondary flow in a square duct," Phys. Fluids, 3, 2734-2745. 24. Winkler, C.M., Sarma, L.R., and Vanka, S.P., "Preferential concentration of particles in a fully developed turbulent square duct flow," Int. J. of Multiphase flow, 30, 27-50, 2004. 25. Reprinted from Int. Journal of Heat and Fluid Flow, Vol. 26, Tafti, D.K., "Evaluating the role of subgrid stress modeling in a ribbed duct for the internal cooling of turbine blades," 92-104, 2005, with permission from Elsevier.
76
26. Tafti, D.K., "GenIDLEST A scalable parallel computational tool for simulating complex turbulent flows", Proc. ASME Fluids Engineering Division, FED-Vol. 256, ASME-IMECE, New York, 2001. 27. Rudinger, G., "Fundamentals of gas-particle flow", Elsevier Scientific Pub. Co., Amsterdam, New York, 1980. 28. Saffman, P. G., "The lift on a small sphere in a slow shear flow.", J. Fluid Mech. 22, 385-400, 1965. 29. McLaughlin, J. B., "Aerosol particle deposition in numerically simulated channel flow." Phys. Fluids A 1, 1211-1224, 1989. 30. Pan, Y., and Banerjee, S., "Numerical simulation of particle interactions with wall turbulence", Phys. Fluids 8 (10), October 1996. 31. Elgobashi, S., and Truesdell, G.C., "Direct simulation of particle dispersion in a decaying isotropic turbulence", J. Fluid Mech. 242, 655, 1992. 32. Wang, Q., Squires, K.D., Chen, M., and McLaughlin, J.B., "On the role of the lift force in turbulence simulations of particle deposition", Int. J. Multiphase flow 23, 749, 1997. 33. 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. 34. 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 35. Wang, Q., and Squires, K.D., "Large eddy simulation of particle deposition in a vertical turbulent channel flow", Int. J. Multiphase Flow 22, 667, 1996.
77
36. Yeh, F., and Lei, U., "On the motion of small particles in a homogeneous isotropic turbulent flow", Phys. Fluids A 3, 2571, 1991. 37. Armenio. V., Piomelli, U., Fiorotto, V., "Effect of the subgrid scales on particle motion", Phys. of Fluids, Vol. 11, No. 10, 1999. 38. Clift, R., Grace, J.R., and Weber, M.E., Bubbles, Drops and Particles, Academic Press, New York, 1978 39. Sadarjoen A, van Walsum Th, Hin A, Post FH, "Particle tracing algorithms for 3D curvilinear grids". Fifth Eurographics Workshop on Visualization in Scientific Computing, 1994. 40. Zhou, Q., and Leschziner, M.A., "An improved particle-locating algorithm for Eulerian-Lagrangian computations of two-phase flows in general coordinates", International Journal of Multiphase Flow, 25, 813-825, 1999. 41. Allievi, A., and Bermejo, R., "A generalized particle search-locate algorithm for arbitrary grids", Journal of Computational Physics, 132, 157-166, 1997. 42. Kim, J., Moin, P. & Moser, R. D., "Turbulence statistics in fully developed channel flow at low Reynolds number", J. Fluid Mech., 177, 133, 1987. 43. Fukagata, K., Zahrai S., and Bark, F.H., "Large Eddy Simulation of particle motion in a turbulent channel flow," ASME Fluids Engineering Division Summer Meeting, June 22-26, 1997, FEDSM'97. 44. Van der Walt, J.P., Nurick, A., "Erosion of dust-filtered helicopter turbine engines Part I: Basic theoretical considerations", J. Aircraft, Vol. 32, No 1, 1995. 45. Wang, L.-P., and Maxey, M.R., "Settling velocity and concentration distribution of heavy particles in homogeneous isotropic turbulence", J. Fluid Mech. 256, 27, 1993.
78
46. Pedinotti, S., Mariotti, G., and Banerjee, S., "Direct numerical simulation of particle behaviour in the wall region of turbulent flows in horizontal channels", Int. J. Multiphase Flow 18, 927, 1992. 47. Rouson, D. W. I., and Eaton, J. K., "Direct numerical simulation of particles interacting with a turbulent channel flow", Proc. of the 7th Workshop on Two-Phase Flow Predictions, edited by M. Sommerfeld (Erlangen, Germany, 1994). 48. Faxn, H., "Die Bewgung einer starren Kugel Lngs der Achse eines mit zaher Flsigkeit fefllten Rohres", Arkiv Mat. Astron. Fys., 17 (27), 1-28, 1923. 49. Wakiya, S.J., Research Report 9, Faculty of Engineering, Niigata Univ., Japan, 1960. 50. Sommerfeld, M., "Modeling of particle-wall collisions in confined gas-particle flows", Int. J. Multiphase flow, Vol. 18, No. 6, 905-926, 1992. 51. Han, J.C., "Recent studies in Turbine Blade Cooling," Int. J. of Rotating Machinery,
10(6), 443-457-2004.
79
Appendix A Modified drag coefficient near solid walls
It has been observed that particles near the wall experience an increase in the drag force. This effect can be incorporated in the particle equation (equation 2.6) by multiplying the drag term by a wall-coefficient C wall which is defined as:
C wall = 1 -
9 1 45 4 1 5 + 3 - - 16 8 256 16
-1
([48])
for wall-parallel directions
9 1 1- + 3 8 2
-1
([49])
for wall-normal directions
where, =
dp 2y
and y is the distance of the centre of the particle from the wall. As
opposed to wall bounded channel flow, the internal cooling ribbed duct has walls in the normal and spanwise direction which introduces additional complexity and increased computational time to incorporate the above effects. Hence for the present study the effect of increased drag coefficient near the wall is neglected.
80
Appendix B Newton's iterative procedure for solving the trilinear function.
Consider x p coordinate of the particle. Similar procedure holds for y p , and z p coordinates of the particle. Regrouping the trilinear function (equation 3.1) for the x coordinate,
x p = f0'
+ f1 *
+ f 2 *
+ f3 *
+
f 4 * * + f 5 * * + f 6 * * + f 7 * * * where,
f 0 '= (x8 + x7 + x6 + x5 + x4 + x3 + x2 + x1 ) / 8 f1 = (x8 - x7 + x6 - x5 + x4 - x3 + x2 - x1 ) / 8 f 2 = (x8 + x7 - x6 - x5 + x4 + x3 - x2 - x1 ) / 8 f 3 = (x8 + x7 + x6 + x5 - x4 - x3 - x2 - x1 ) / 8 f 4 = (x8 - x7 - x6 + x5 + x4 - x3 - x2 + x1 ) / 8 f 5 = (x8 - x7 + x6 - x5 - x4 + x3 - x2 + x1 ) / 8
f 6 = (x8 + x7 - x6 - x5 - x4 - x3 + x2 + x1 ) / 8
f 7 = (x8 - x7 - x6 + x5 - x4 + x3 + x2 - x1 ) / 8
Let,
f x = f 0 '- x p + f1 *
(
)
+ f 2 *
+ f3 *
+
f 4 * * + f 5 * * + f 6 * * + f 7 * * *
f x = 0 = f0 + f1 * + f 2 * +
f 3 * + f 4 * * + f 5 * * + f 6 * * + f 7 * * *
where, f 0 = f 0'- x p Similar equations can be written for y p , and z p .
g y = 0 = g0
+ g1 * + g 2 * +
g 3 * + g 4 * * + g 5 * * + g 6 * * + g 7 * * *
81
hz = 0 =
h0
+ h1 * + h2 * +
h3 * + h4 * * + h5 * * + h6 * * + h7 * * *
where,
g 0 = ( y8 + y7 + y6 + y5 + y 4 + y3 + y2 + y1 ) / 8 - y p h0 = (z8 + z 7 + z 6 + z5 + z 4 + z3 + z 2 + z1 ) / 8 - z p
Similarly g1 , g 2 , g 3 , g 4 , g 5 , g 6 , g 7 , and h1 , h2 , h3 , h4 , h5 , h6 , h7 can be written analogous to
f 1 to f 8 .
Equations (2)-(4) can be solved by Newton's iterative method by letting
df x = - f x , dg y = - g y , dh z = -h z
From equation (2),
df x = f1 * d f 4 * * d f 5 * * d f 6 * * d + f 2 * d + + f 3 * d +
+ f 4 * * d +
f 5 * * d +
+ f 6 * * d +
f 7 * * * d + f 7 * * * d + f 7 * * * d
Regrouping,
df x = ( f1 + f 4 * + f 5 * + f 7 * * ) * d +
( f 2 + f 4 * + f 6 * + f 7 * * ) * d + ( f 3 + f 5 * + f 6 * + f 7 * * ) * d
=-f x =0
Similarly,
dg y = (g1 + g 4 * + g 5 * + g 7 * * ) * d +
( g 2 + g 4 * + g 6 * + g 7 * * ) * d + ( g 3 + g 5 * + g 6 * + g 7 * * ) * d
= -g y = 0
82
dh z = (h1 + h4 * + h5 * + h7 * * ) * d +
(h2 + h4 * + h6 * + h7 * * ) * d + (h3 + h5 * + h6 * + h7 * * ) * d
z
= -h
= 0
Equations (5)-(7) can be solved simultaneously for d , d , d . An initial guess of
= 0, = 0, = 0 is used to start the iteration.
After each iteration,
new = old + d new = old + d new = old + d
If d , d , and d 10 -6 , then the solution is considered to be converged. If ( , , ) [- 1,1] then the point is inside the hexahedron l .
83
Appendix C
A pseudo FORTRAN program combining the particle location and cell search techniques:
Consider that n particles are at cell locations (in, jn, kn) and block mn at time t
Do n=1, total number of particles
Starting cell for ntb particle = i, j, k of the particle at time t
100
search_number = search_number + 1 If (search_number max_searches) then Go to Step 1 Else Go to 500 End if
Step 1:
Calculate f 0 , f1, f 2 , f 3 , f 4 , f 5 , f 6 , f 7
g 0 , g1 , g 2 , g 3 , g 4 , g 5 , g 6 , g 7 h0 , h1 , h2 , h3 , h4 , h5 , h6 , h7
Step 2:
Do k = 1, total number of iterations (max 20) Set , , = 0 Solve the following equations for d, d, d:
( f1 + f 4 * + f 5 * + f 7 * * ) * d + ( f 2 + f 4 * + f 6 * + f 7 * * ) * d + ( f 3 + f 5 * + f 6 * + f 7 * * ) * d = 0
84
(g1 + g 4 * + g5 * + g 7 * * ) * d + ( g 2 + g 4 * + g 6 * + g 7 * * ) * d + ( g 3 + g 5 * + g 6 * + g 7 * * ) * d = 0 (h1 + h4 * + h5 * + h7 * * ) * d + (h2 + h4 * + h6 * + h7 * * ) * d + (h3 + h5 * + h6 * + h7 * * ) * d = 0
new = old + d new = old + d new = old + d
if d , d , and d 10 -6 then solution has converged. Step 3:
Check if absolute (, , ) 1 then, particle is located.
If any of the condition is violated then
Get new cell using the following decisions: if( +1)new_i = min(i + 1, max_i) if( -1)new_i = max(1, i-1) if( +1)new_j = min(j + 1, max_j) if( -1)new_j = max(1, j - 1) if( +1)new_k = min(k + 1, max_k) if( -1)new_k = max(1, k - 1) Go to100
End if
Repeat the above procedure till the particle is located.
500
Particle is considered out of the current block. Its information has to be sent to the corresponding block. This is discussed in the implementation section. Go to the next particle.
End do (total number of particles)
85
Appendix D Lagrangian Particle Tracking Algorithm
86
Discretize the physical domain into blocks. Decide the number of processors to be used for the calculation. Calculate the number of blocks/processor, and distribute the corresponding blocks to the processors.
Solve the Eulerian flow field No injection of particles Check if the flow field has reached a statistically steady state
NO
YES
Inject particles in the domain. Give them an initial i, j, k location and a block number
Time loop with both phases begins
Start the particle loop Set search = 0 Apply the relevant wall collision scheme
Check if the particle is near the wall
YES
NO
87
100 search = search + 1 NO
search Max_search
YES
Search the particles starting from its previous cell location on the current block
Use the particle search technique discussed in section 3
YES
Set `found = 1'
Is the particle located?
NO
Set `found = 0'
Get the local block & neighboring processor number where the information of this particle needs to be sent.
NO
Did the particle leave through an open boundary?
YES
If the particle left the current block through a periodic face then adjust its coordinates according to the periodic lengths. Particle is considered out of the domain and is written in a lost file
Store the particle information
Store the particle information in a send buffer
Go to the next particle End of particle loop
88
Reorder the particle numbers such that particles with `found=1' are together. Nptime=n = Nptime=n-1 number of particles not found or lost from the current processor.
Send all the particle information stored in the send buffer to the respective processors. Receive particle information from the neighboring processors.
YES
Were there any particles received from the neighboring processors?
Repeat the above process. Go to 100
NO
Check if the particles on the remaining processors are located
NO
YES
Interpolate the fluid velocities at the particle location
Update the particle position by integrating the following ODE's:
d xp dt
= up
mp
d up dt
=
Fp
Go to the next time step End time loop
89
APPENDIX E Wall collision model
During the time advancement of particles, situations arise where a particle may contact the boundary. The basic idea of a wall collision model is to take into consideration the fact that when a particle collides with a boundary it loses a fraction of its momentum before being re-introduced in the flow. The term which relates this loss of momentum of an impinging particle is called the restitution coefficient (eR), which can be defined for the normal and tangential direction. The normal restitution coefficient is defined as the amount of momentum retained by the particle in a direction normal to the wall after colliding with the boundary. Mathematically, it is the ratio of the normal component of the particle velocity after collision to the normal velocity of the particle velocity before collision. For example, in a wall-bounded channel flow with walls in the normal direction the normal restitution ratio is given by, eR = Vafter collision / Vbefore collision. Similar description holds true for the restitution coefficient in the tangential direction. It should be noted that a particle is said to contact a boundary when it's within a specified distance from the wall. For most of the studies involving particle laden flows, a particle is considered to contact the wall when its center is one radius from the wall. A number of schemes may be applied for the above situation. A general description of different models used in particle transport studies is given below
90
A general wall collision model can be defined as follows:
Consider the velocity of a particle after impact to be superscripted by a "^". The following conditions need to be satisfied.
U n = - e R U n
U n t = U n t - (1 + e R ) U n
^ ^ ^
^
[17]
Here n is the unit vector normal to the boundary and n t is the unit vector tangential to the boundary. is the coefficient of solid friction. For = 0 , the collision is frictionless. Depending on the values of eR, n, nt, and different wall-collision models can be defined. Some of the models are defined below. For all the models described it is assumed that = 0 , and the tangential component of velocity vector remains unchanged.
Model 1: Deposition model
This model assumes that when the particle is within the specified distance it deposits on the wall. The particle is removed from the calculation or a new particle is introduced in the flow to maintain the particle flow rate.
Model 2: Irregular bouncing model
An irregular bouncing model can be classified into two general categories: In the first category, the model assumes a constant restitution ratio. It implies that all the particles retain the same fraction of their normal momentum after colliding with the wall.
91
In the second category, the model uses a restitution ratio which is a function of the material properties, material roughness, and the particle impingement velocity and the impingement angle. Details of such a model can be found in Sommerfeld [50].
Model 3: Elastic collision model
This model assumes that particles reflect elastically from the wall, which implies they retain all of their normal momentum. Hence the restitution ratio, for this model is 1 as opposed to 0 for a deposition model.
All the above models can be summarized based on the restitution ratio:
0 - Deposition Model 1 - Elastic Model
eR =
Normal f impact , U impact , Material - Irregular Model
(
)
Wall collision model for the present study
The treatment of boundary-particle interactions for the present simulations is described below: If a particle impacts the wall with a normal velocity U 2p , its normal velocity after impact is given by U 2p = - eR U 2p = -U 2p (elastic collision). The other components of the velocity remain unchanged.
92
VITA Anant Shah
Anant Shah is a proud Indian born in Ahmedabad on March 7th, 1982. He received his B.E. in Mechanical Engineering from L.D. College of Engineering in June of 2003 and then joined the Masters program at Virginia Tech in Jan of 2004. He is thankful to the Almighty for giving him an opportunity to pursue higher education.

**Find millions of documents on Course Hero - Study Guides, Lecture Notes, Reference Materials, Practice Exams and more. Course Hero has millions of course specific materials providing students with the best way to expand their education.**

Below is a small sample set of documents:

Virginia Tech - ETD - 03272001

Appendix I(VT Project Box Construction Manual)70Construction of a VT Project BoxTable of Contents List of Parts for the VT Project Box .73 Construction of the VT84 board..74 VT84 attachment and jumpering.77 VT84 wire connections and connectors

Virginia Tech - ETD - 01042001

MATERIAL CHARACTERIZATION AND LIFE PREDICTION OF A CARBON FIBER/THERMOPLASTIC MATRIX COMPOSITE FOR USE IN NON-BONDED FLEXIBLE RISERSBlair Edward RussellA thesis submitted to the Faculty of the Virginia Polytechnic Institute and State University i

University of Texas - AFFILIATES - 2007

Equation of State (EOS) Compositional ModelingS. G. Thomas M. Delshad M. F. WheelerCenter for Subsurface Modeling Institute of Computational Engineering and Sciences The University of Texas at AustinIndustrial Affiliates Meeting, Oct. '07Outli

University of Texas - AFFILIATES - 2007

A MULTISCALE MORTAR METHOD FOR EFFICIENT SOLUTION OF THE HETEROGENEOUS PRESSURE EQUATION Todd Arbogast Department of Mathematics and Center for Subsurface Modeling, Institute for Computational Engineering and Sciences (ICES) The University of Texas a

University of Texas - AFFILIATES - 2005

A Geomechanics Model Coupling Multipoint Flux and Galerkin FEMXiuli Gai Mary F. WheelerCenter for Subsurface Modeling ICES, UT AustinIvan YotovUniversity of PittsburghCSM, ICES UT AUSTINOutline Motivation Coupling MFE MPFA and Galerkin FEM

University of Texas - AFFILIATES - 2008

Domain Decomposition for Problems in ElasticityGergina Pencheva1 Vivette Girault2 Tim Wildey1 Mary F. Wheeler11Institute for Computational Engineering and Sciences, The University of Texas at Austin 2 Universite Pierre et Marie Curie, Paris VI,

University of Texas - AFFILIATES - 2006

Research Directions of Solvers for Flow in Porous MediaHector KlieCenter of Subsurface Modeling Institute of Computational Engineering and Sciences The University of Texas at Austin2006 Industrial Affiliates MeetingAustin, October 16-17, 2006O

University of Texas - AFFILIATES - 2005

Modeling Flow in Fracture Reservoirs: From Discrete to ContinuumAdolfo A. Rodriguez, Hector Klie, Mary Wheeler, Xiuli Gai CSM and Horacio Florez, USB.ObjectiveObtaining a continuum representation of a discrete fractured system represented by a te

University of Texas - AFFILIATES - 2008

Domain Decomposition for Poro-ElasticityV. Girault Laboratoire Jacques-Louis Lions, UPMC, Paris 6 G. Pencheva, M. Wheeler, T. Wildey Center for Subsurface Modelling, ICES, UT Austin Model and motivation Primal variational formulation Domain Decompos

University of Texas - AFFILIATES - 2005

A Numerical Technique for Simulating Non-planar Evolution of Hydraulic FracturesJaroon Rungamornrat and Mary F. WheelerCenter for Subsurface Modeling Institute for Computational Engineering and Sciences The University of Texas at AustinMark E. Me

University of Texas - USNCCM - 8

Computational Electromagnetics: Boundary Integral Equation Solutions from the Very Large to the Very SmallW.C. Chew Center for Computational Electromagnetics and Electromagnetics Laboratory Department of Electrical and Computer Engineering Universit

University of Texas - USNCCM - 8

rsph|ri~qqrvhEYqr6|YwW(~fqsti|slk(hYPYhfEY|%Y|slhfsdqrihqE(Tr9rYfi~s{9(kh(~qhg qe d qe t lg t j ln j d x t g l j p{ q e { e xg l t j d { j qg p q{ e l j g { q y j d e y q g j { e l { p { e j l g d x {e xg l t j d {

University of Texas - USNCCM - 8

LEVEL SET METHOD SIMULATION OF SURFACE & DEFECT EVOLUTION IN MATERIALS David J. Srolovitz Department of Mechanical and Aerospace Engineering Princeton University Princeton, NJ 08544-5263 USA srol@princeton.edu The level set method provides a robust a

University of Texas - M - 427

M427K Handout: Separable and Exact EquationsSalman Butt January 26, 2006This handout describes the procedure used to solve two specific types of first order differential equations: separable and exact equations Separable Equations Observe that we

University of Texas - M - 427

M427K Handout: PDE's: The Heat and Wave EquationsSalman Butt May 4, 20061. Next week, I will continue to hold office hours. Regular discussion will also be office hours. The Heat Equation Recall the setup and solution for the heat equation:Let's

Virginia Tech - ETD - 04192000

CHAPTER 5A PARAMETRIC STUDY OF THE FOUR-BOLT WIDE MOMENT END-PLATE CONNECTION5.1INTRODUCTIONThis chapter presents the results of a parametric study of the four-bolt widemoment end-plate connection. The behavior of this connection is quite co

Virginia Tech - ETD - 91497

Chapter 1IntroductionSince the invention of the column flotation process and its successful introduction in mineral processing plants, there has been a lot of interest in understanding the various mechanisms that take place during operation. Column

Virginia Tech - ETD - 062599

"The center must combine the maximum of congestion with the maximum of light and space." - Koolhaas, p.178COMMERCIALNatural ventilation patterns must be considered at every scale from the urban to the domestic as an alternative to artificial clim

Virginia Tech - ETD - 11012000

Chapter 1 Introduction1.1Statement of the ProblemSite characterization for geotechnical purposes is generally defined as "the identification and description of the subsurface strata within the areas of influence of a project" (Mitchell 1978). I

Virginia Tech - CS - 5504

Tomasulo Loop ExampleLoop: LD MULTD SD SUBI BNEZ F0 F4 F4 R1 R1 0 F0 0 R1 Loop R1 F2 R1 #8 This time assume Multiply takes 4 clocks Assume 1st load takes 8 clocks, 2nd load takes 1 clock Clocks for SUBI, BNEZ Show 2 iterationsLoad=8,1Add=3

University of Texas - CH - 301

Worksheet 11: The thermodynamics of chemical and physical processes Chemical Reaction H (kJ) BE (kJ) nga w (kJ) nsystem sign CH4g + 2O2g CO2g + 2H2Og Explain: Combustion reaction, would predict large calc. heat and spontaneous reaction. n = 0 means

University of Texas - CH - 301

Laude's CH301 Worksheet 1: EMR and Waves (Sections 1.1-1.9)(The textbook referenced is Jones & Atkins's Chemical Principle, 3rd edition) 1. Classical physics predicts that light can be transferred from object to object independent of frequency. This

University of Texas - CH - 301

University of Texas - CH - 302

CH 302 Spring 2008 Worksheet 14(Note:Tthere are specific rules about which groups get priority in numbering, but you aren't responsible for knowing them. So if your name has different numbers from mine, try flipping your numbering system around on t

University of Texas - CH - 302

CH302 Spring 2008 Worksheet 5 14 questions involving simple acid base equilbria and the approximations that make them simple. 1. The only water equilibrium for which we make no approximations is the case of pure water (amazing how simple something is

University of Texas - CH - 302

CH 302 Spring 2008 Worksheet 6 Key1. You have a 750 mL solution of 0.1 M methylamine. You can't find the Kb for methylamine but notice that the Ka for its conjugate acid is 1 x 10-9. What is the pH of the methylamine solution? Answer: Kw = KaKb = 1

University of Texas - CH - 302

CH 302 Spring 2008 Worksheet 61. You have a 750 mL solution of 0.1 M ammonia (Kb = 10-5). What is the pH of this solution? 2. You decide to titrate it against 1 M hydrochloric acid. When you've added 25 mL of the HCl to the solution, what is the pH?

University of Texas - CH - 302

CH 302 Spring 2008 Worksheet 81. NaH2PO4 (conc. = CNaH2PO4) is dissolved in water. Write the mass balance equation for this system.2. Write the charge balance equation for the solution in question 1.3. Write the charge balance equation for a sol

University of Texas - CH - 302

All of this is intended to be done without the aid of a calculator. All of the calculations are designed such that approximating should be straight-forward and produce a correct result. 1. Based on the physical constants involved, which colligative p

University of Texas - NQTHM - 1992

#| Copyright (C) 1994 by Alex Bronstein and Carolyn Talcott. Reserved. All RightsYou may copy and distribute verbatim copies of this Nqthm-1992 event script as you receive it, in any medium, including embedding it verbatim in derivative works, prov

University of Texas - NQTHM - 1992

#| Copyright (C) 1994 by Alex Bronstein and Carolyn Talcott. Reserved. All RightsYou may copy and distribute verbatim copies of this Nqthm-1992 event script as you receive it, in any medium, including embedding it verbatim in derivative works, prov