<%BANNER%>

Computational Modeling of Glow Discharge-Induced Fluid Dynamics


PAGE 1

COMPUTATIONAL MODELING OF GL OW DISCHARGE-INDUCED FLUID DYNAMICS By BALAJI JAYARAMAN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2006

PAGE 2

Copyright 2006 by Balaji Jayaraman

PAGE 3

This document is dedicated to my parents and all my teachers from whom I have learnt all my life.

PAGE 4

iv ACKNOWLEDGMENTS I would like to thank my advisor, Professo r Wei Shyy, for his advice, inspiration, patience and mentoring. He constantly prov ided insightful suggestions when I needed help and motivation. I would also like to express gratitude to my committee members: Dr. Siddharth Thakur, Professors Corin Segal, Toshi Nishida, Lou Cattafesta and William Lear. A special mention needs to be made of Dr.Tha kur for his suggestions regarding the thesis content. I have also gained significantly fr om the various technical discussions which I have had with him. Finally, I am forever i ndebted to my parents for their understanding and support when required. This acknowledg ement will not be complete without the mention of the role played by my rese arch group members for their support and friendship.

PAGE 5

v TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES............................................................................................................vii LIST OF FIGURES...........................................................................................................ix NOMENCLATURE.........................................................................................................xii ABSTRACT.....................................................................................................................xi v 1 INTRODUCTION............................................................................................................1 1.1 Overview of Plasma Generation and Di electric Barrier Discharges (DBD’s).......1 1.2 Plasma Actuator Operating Mechanism.................................................................3 1.3 Numerical Challenges.............................................................................................6 1.4 Overview of Present Approach...............................................................................7 1.5 Outline of the Dissertation......................................................................................9 2 LITERATURE REVIEW...............................................................................................12 2.1 Review of Experimental Discharge Plasma Studies.............................................12 2.1.1 Early Discharge Studies.............................................................................13 2.1.2 Mechanism of Plasma-Induced Flow Applications....................................15 2.1.3 Aerodynamic and Control Strategies Using Glow Discharges..................20 2.2 Review of Discha rge Modeling Studies...............................................................25 2.2.1 Phenomenological Modeling of the Plasma Actuator................................26 2.2.2 First Principles-Based Plasma Modeling Approaches...............................27 2.2.2.1 Review of dielectric barrie r discharge (DBD) modeling studies.....27 2.2.2.2 Review of plasma actuator modeling studies...................................33 3 NUMERICAL MODELING OF GLOW DISCHARGE PLASMA..............................36 3.1 Plasma Modeling Hierarchy.................................................................................36 3.1.1 Kinetic Modeling Approach.......................................................................38 3.1.2 Hybrid Modeling Approach.......................................................................42 3.1.3 Fluid Modeling Approach..........................................................................44 3.2 Numerical Modeling of Glow Discharge Plasmas...............................................47 3.2.1 Plasma Discharge Governing Equations....................................................48

PAGE 6

vi 3.2.2 Numerical Algorithm..................................................................................59 3.2.3 1-D Parallel Plate Disc harge Modeling in Helium.....................................74 3.3 Integration of Plasma Dynamics with Fluid Dynamics........................................78 4 MODELING FLUID FLOW WI TH PLASMA EFFECTS............................................81 4.1 Body Force Formulation.......................................................................................82 4.2 Navier-Stokes Equations with Body Force...........................................................84 4.3 Modeling a Plasma Actuator Usi ng a Phenomenological Force Model...............85 4.3.1 Plasma Generation in an Asy mmetric Electrode Configuration................85 4.3.2 Linearized Electric Body Force..................................................................88 4.3.3 Problem Description...................................................................................91 4.3.4 Effect on Flow Structure............................................................................93 4.3.5 Effect on Heat Transfer..............................................................................97 4.3.6. Thrust Estimation......................................................................................98 4.3.7 Effect of Plasma on Aerodynamics..........................................................102 4.4 Self-Consistent Fluid-Mode l-Based Plasma Dynamics......................................108 4.4.1 Discharge Actuator Setup.........................................................................109 4.4.1.1 Boundary condition........................................................................111 4.4.1.2 Plasma species initial condition.....................................................112 4.4.2 Plasma Structure.......................................................................................113 4.4.3 Geometric Effects.....................................................................................125 4.4.3.1 Impact of electrode spacing............................................................125 4.4.3.2 Impact of lower electrode size.......................................................127 4.4.4 Impact of Applied Voltage and Frequency..............................................128 4.4.4.1 Frequency sensitivity......................................................................129 4.4.4.2 Applied voltage sensitivity.............................................................130 4.5. Summary.........................................................................................................131 5 SUMMARY AND FUTURE WORK..........................................................................134 5.1 Summary and Conclusions.................................................................................134 5.2 Future Work........................................................................................................137 LIST OF REFERENCES.................................................................................................139 BIOGRAPHICAL SKETCH...........................................................................................149

PAGE 7

vii LIST OF TABLES Table page 1 Representative EHD operating mechanisms...................................................................18 2 Summary of experimental plasma actuator studies........................................................22 3 Summary of turbulent flow control studies....................................................................25 4 Summary of low frequency plasma modeling studies....................................................31 5 Summary of high frequency plasma modeling studies..................................................33 6 Summary of plasma-induced fl uid flow modeling studies.............................................35 7 Kinetic model studies of discharge plasmas...................................................................40 8 Summary of hybrid modeling studies.............................................................................43 9 Plasma chemistry considered in a helium discharge.......................................................46 10 Representative time scales of various processes in a 1-D helium discharge model operating at 10KHz, 1.5KV rms, calculate d based on values of parameters...........60 11 Summary of governing equations used in different studies..........................................61 12 Different types of boundary conditions for the ions and electrons used in modeling studies.......................................................................................................................7 6 13 Change in drag coefficient and drag force with Reynolds number............................104 14 Change in lift coefficient and lift force with Reynolds number.................................104 15 Effect of electrode position on aerodynamic and heat transfer properties..................107 16 Summary of property models employe d for the He discharge simulation.................111 17 Summary of boundary conditions for the different variables.....................................113 18 Domain averaged force over the volta ge cycle for different waveforms....................125 19 Domain averaged force over the voltage cycle for different configurations..............127

PAGE 8

viii 20 Force dependence on the lower electrode size............................................................128 21 Domain averaged force over the volta ge cycle for different frequencies...................130 22 Mean domain averaged force over the time cycle for different voltage.....................130

PAGE 9

ix LIST OF FIGURES Figure page 1 Experimental photograph illustrating the EHD effects on a laminar jet of titanium tetrachloride vapor....................................................................................................11 2 Generation of glow discharg e in various arrangements..................................................17 3 1-D Numerical modeling results showing diffe rent regions of the discharge gap at the instant when the current is maximum......................................................................18 4 A 3-D surface plot of photomultiplier light out put from the plasma as a function of time and chordwise distance, for a single AC cycle of the applied voltage.....................21 5 Separation control on a NACA66-18 wing at a Re=79000, 16 degrees angle of attack.23 6 Streamer front propagation at diffe rent instants in a DC voltage..................................35 7 Levels of plasma modeling.............................................................................................38 8 Parallel plate arrangement for modeling DBD discharge...............................................76 9 Evolution of discharge current over one period of the discharge cycle..........................76 10 Evolution of the gap voltage over one full cycle..........................................................77 11 Comparison of peak current species density.................................................................78 12 Schematic illustration of the plas ma-fluid flow modeling framework.........................83 13 Schematic of the actuator arrangement with approximate shape of electric field lines. The thickness of the electrodes d is 0.1 mm............................................................86 14 The line A-B constitutes the plasma fl uid boundary using linear approximation. The electric field strength outside this line is not strong enough to ionize the air..........89 15 The computational domain used for calcula tions shows the flat plate with the upper electrode...................................................................................................................93 16 Velocity profiles compared at the four different stations.............................................95 17 Effect of the various parameters as s een form the velocity profiles at ST4.................96

PAGE 10

x 18 Effect of the body force on the wall shear at ST4.........................................................96 19 Temperature profile along the wall for Re=1027 and Pr=0.7.......................................97 20 Comparison of normalized heat flux fo r different voltages and frequencies................98 21 Comparison of the experimental and calcu lated plasma force for different voltages.101 22 A schematic of the computational domain with the flat plate and the electrode. The electrode has negligible thickne ss in the above arrangement................................103 23 Drag variation with angle of attack for various Re. The plasma is operated at 4KV rms and 3KHz................................................................................................................103 24 Variation of lift with Re for different angl es of attack. The plasma is operated at 4KV rms and 3 KHz........................................................................................................106 Nusselt number variation with angle of attack for various Reynolds number. The plasma is operated at 4KV rms and 3KHz frequency........................................................106 26 Streamlines showing flow separation for Re=333 and an angle of attack twenty degrees for different positi ons of the electrode......................................................108 27 A representative 2-D asymme tric discharge arrangements........................................110 28 Sensitivity to initial conditi on (20KHz, 1KV applied voltage)..................................114 29 Photo intentisty measurements to depict the observed discharge evolution over a timecycle.......................................................................................................................115 30 Time dependent plasma behavior (5 KHz, 1KV applied voltage)..............................115 31 Ion density evolution in the domain ove r time (5 KHz, 1KV applied voltage)..........116 32 Force evolution over one time cycle (5 KHz, 1KV applied voltage).........................117 33 Asymmetric sawtooth voltage waveforms used in the experimental and numerical studies.....................................................................................................................122 34 Experimental results comparing the discharge observed from the two different sawtooth waveforms...............................................................................................123 35 Impact of waveform on the species number density...................................................123 36 Thrust measurements showing the more efficient......................................................125 37 Axial force integrated over the domain with time (Fx (dynes)) for different waveforms..............................................................................................................125

PAGE 11

xi 38 3-Different configurations to stud y the impact of electrode overlap..........................126 39 Time evolution of species number density over one time-cycle (20KHz, 1KV voltage signal).....................................................................................................................126 40 Impact of the geometric configuration on the domain integrated force field (20KHz, 1KV voltage signal)...............................................................................................127 41 Three different lower electrode lengths mode led to investigate the effect of increased surface area on the discharge..................................................................................128 42 Impact of lower electrode size....................................................................................128 43 Charge density and force variation over time for different frequencies.....................129 44 Charge density and force variatio n over time for different voltages..........................130 45 Fx power-law dependence on voltage is co mpared for two different numerical models131

PAGE 12

xii NOMENCLATURE jn,in ,en Species number density N Gas number density ju,iv ,ev Velocity Potential E Electric field T Temperature of various species ,, iej Species internal energy i Ratio of specific heats ,, iej Species mobilities Di,e,j Species diffusivities mi,e,j Species mass Rc, Rd Creation/destruction of species ,cD Creation/destruction frequencies KC,KD /,/cDNN Rij Number density gain/loss from reaction Sij Number density gain/loss from inelastic collisions ij Elastic collision frequencies qj Charge carried by jth species e Elementary electronic charge

PAGE 13

xiii 0 Permittivity of free space K Boltzmann constant Rm Momentum loss from collision (,) mie Momnetum loss collision frequency Km /mN Pe Electron pressure Ren Energy loss term Qe Electron heat conduction T Species thermal energy w Driving/operating time scale c Creation/ionization characteristic time dr Species drift characteristic time diel Dielectric relaxation characteristic time Ci Ion/species sound speed Ce Electron sound speed Cfluid Fluid sound speed Ftave Time-averaged body force Fcoll Force resulting from species collision losses

PAGE 14

xiv Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATIONAL MODELING OF GL OW DISCHARGE-INDUCED FLUID DYNAMICS By Balaji Jayaraman August, 2006 Chair: Wei Shyy Major Department: Mechanic al and Aerospace Engineering Glow discharge at atmospheric pressure using a dielectric ba rrier discharge can induce fluid flow and operate as an actuator for flow control. The largely isothermal surface plasma generation realized above can modify the near-wall flow structure by means of Lorentzian collisions between the ionized fluid and the neutral fluid. Such an actuator has advantages of no moving parts, performance at atmospheric conditions and devising complex control st rategies through the applied voltage. However, the mechanism of the momentum coupling between the plasma and the fluid flow is not yet adequately understood. In the present work a modeling framework is presented to simulate athermal, non-equilibrium plasma di scharges in conjunction with low Mach number fluid dynamics at atmospheric pressure The plasma and fluid species are treated as a two-fluid system exhibiting a few decades of length and time scales. The effect of the plasma dynamics on the fluid dynamics is devised via a body force treatment in the Navier-Stokes equations. Two different approaches of different degrees of fidelity are

PAGE 15

xv presented for modeling the plasma dynamics. The first approach, a phenomenological model, is based on a linearized force distribution approxima ting the discharge structure, and utilizing experimental guidance to deduce the empirical constants. A high fidelity approach is to model the plasma dynamics in a self-consistent manner using a first principle-based hydrod ynamic plasma model. The atmospheric pressure regime of interest here enables us to employ local equilibrium assumptions, signifying efficient collisional energy exchange as against thermal heating from inelastic collision processes. The time scale rati os between convection, diffusion, and reaction/ionization mechanisms are O(107), making the system computationally stiff. To handle the stiffness, a sequential finitevolume operator-splitting algorithm capable of conserving space charge is developed; the approach can handle time-step sizes in the range of the slowest species convection time-scale. The Navier-Stokes equations representing the fluid dynamics are solved using a well-establish ed pressure-based algorithm. A one-dimensional two-species plasma model was employed as a test case for validation purposes. The momentum coupling is primarily caused by the combination of factors which include discharge chemistry, individual species tr ansport properties, geometric construction and the nature of the insulator and electrode material. Overall, the paraelectric momentum coupling mechanism is due to the cumulative effect over time of the force field in the domain, as seen from our computations. Parametric studies conducted on the operati ng variables such as voltage. Frequency and geometric arrangements indicated strong agreement with the observed experimental work. The applied voltage indicated a pow er-law dependence on the voltage for the measured force in the domain.

PAGE 16

1 CHAPTER 1 INTRODUCTION This dissertation deals with the computat ional modeling of discharge processes and the accompanying plasma-fluid interaction with focus on the asymmetric electrode dielectric barrier plasma actuator and its appli cations. In particular we focus our attention on the operating mechanism and thermo-fluid effects of the discharge plasma actuator (see Figure 1) which generates efficient su rface plasma using a dielectric barrier arrangement capable of interacting with th e near-wall flow structure through momentum injection. In this ch apter we will present (i) a brie f overview and background of the discharge operating processes in the actuato r (ii) the mechanism of operation of the plasma actuator and momentum injection to the flow (iii) the state of the research of the discharge-induced fluid mechanics and their numerical modeling. Finally we will present an outline of the various issues addressed and the overall organi zation of the present dissertation. 1.1 Overview of Plasma Generation and Dielectric Barrier Discharges (DBD’s) The plasma is defined as a quasineutral particle system in the form of gaseous or fluid-like mixtures of free electrons and ions frequently also cont aining neutral particles [1]. Based on the extent of ionization, one obs erves either fully ionized or weakly ionized plasma. Typically, such ionization is realized in gases when a sufficiently large potential is applied under conditions of low pressures such as in vacuum discharges [1]. However, the advent of dielectric barrier discharges (or in short DBD’s) presents an efficient mode

PAGE 17

2 of volume plasma generation at atmospheric pr essure conditions at which the potential for applications are enormous. Dielectric barrier discharges have been studied since mid-nine teenth century when they were identified in the context of ozone generation (Siemens [2]). Since then, they have been the focus of extensive interest generated by their use in a wide variety of applications such as surface treatment, exci mer ultraviolet lamps [3], large area flat plasma display panels [3] etc. Recent applications have taken this a step further by using the concept of the DBD to generate surface electro-hydrodyna mic (EHD) flows [4] which have potential for flow control, therma l management among other applications. DBD’s are very attractive for industrial a pplications because th ey can provide nonequilibrium plasma conditions at atmospheri c pressure. Typically, discharges at high pressure tend to develop into filamentary or arc type discharges which are unsuitable for certain applications. However, in recent year s, homogeneous glow discharges have also been obtained at atmospheric pressure which can be used for large volume applications, surface treatment. The DBD provides a mechan ism for generating uniform discharge at atmospheric pressures without arc formation. Typical discharge configurations use planar or cylindrical arrangements with a dielectric layer (e.g. glass, quartz, polymers) placed between the electrodes. The discharge gaps us ually range from less than a millimeter to a few centimeters depending on the application [3]. The most interesting property of these discharges is that, at about atmospheric pressure, the breakdown is initiated in the form of a large number of independent current fila ments which are known as microdischarges. These discharges are characterized as weakly ionized plasma channels which resemble high pressure glow discharges. Due to the ch arge build-up at the dielectric surface in a

PAGE 18

3 short time after breakdown, the electric field in the discharge region is reduced, thus making it self-limiting. The short duration of the discharge current limits the energy dissipation and consequently results in litt le gas heating (hence, the plasma being athermal) and reduces power consumption. In recent years, homogeneous glow discharges have been generated at atmospheri c pressures in DBD configurations [4-7]. Such a discharge generation mechanism has been used by Roth [5] to devise a source of uniform and stable surface plasma. It has b een known that EHD is one of the efficient modes of flow control in an ionized gas [4]. The EHD induced flow is schematically illustrated below in Figure 1. The figure shows the effect of the plasma generated EHD flow on a surface using smoke flow visualiza tion as reproduced from [4]. The dominant effect on the flow can be seen from the bendi ng of the jet and re-e merging as a wall jet. This can be easily used to inject momentum in regions of adverse pressure gradients, boundary layers etc. Significant attention has been devoted to boundary layer flow control using methods such as synthetic jets [8]. The present method has the advantage of non-moving parts, small size, performance at atmospheric conditi ons and capacity for complex unsteady flow actuation. But, issu es such as power consumption, force generation mechanism and optimum operating c onditions have still no t been adequately addressed. In this dissertation, we will fo cus on developing numerical models (both first principle-based high-fidelit y as well as phenomenological approaches) to study the momentum generation and the thermo-f luid potential of these devices. 1.2 Plasma Actuator Operating Mechanism The exact mechanism of electrohydr odynamic (EHD) flow generation is not sufficiently understood although the concept behind the force generation is widely perceived to be that of Lorentzian collisi ons. The paraelectric effects arising from the

PAGE 19

4 electric field gradient acceler ate the ions which transfers mo mentum to the neutral fluid [4]. The effect of the plasma on the fluid is like a localized body force [9] on the neutral particles. Most of the investigations so fa r have been based on e xperimental observations and empirical arguments using approximate models [4,9] until recently. The above mentioned EHD effects involve the interac tion between the disc harge physics and the near-wall flow physics. In orde r to systematically explain th e plasma fluid interaction, detailed modeling studies are required. Such studies should attempt to address the complex plasma physics and its interaction w ith flow physics. The coupled plasma-fluid problem is inherently non-linear by nature, ju st because one is forced to treat the fluid and plasma as separate entities which interact It will be shown later in chapters 3 and 4 that the disparity in the fluid and plasma time scales can be used to decouple their solution methodologies. Moreover, the EHD flow can itself be of a different length scale as compared to the mean flow. Resolving these multiple scales is a challenging computational task. Also, the modeling of plasma dynamics in itself involves a convection-diffusion-reaction system operatin g over widely varying time and length scales. The interaction between th e aerodynamic motion of an electrically conducting medium and an electromagnetic field is of great interest in astro-phys ics; interstellar gas masses etc and is studied under magneto-aer odynamics. The potential of electromagnetic effects in enhancing aerodyna mics performance have been discussed and reviewed by Shang et al. [10]. Mechanism of flow control using plas ma studied in non-equilibrium hypersonic flows is fundamentally different from the EHD type mechanism described above. Here,

PAGE 20

5 the ionized flow modified in the presence of electromagnetic fields represents the interplay of aerodynamics, electromagnetism, chemical physics, and quantum physics. This is encountered in hypersonic flight s where, the air mixture bounded by the bow shock consists of highly energized internal energy modes and at high temperatures around 5000 K ionization is enc ountered. This ionized air mixt ure will interact with the electromagnetic field causing Lorentz force e ffects to affect the aerodynamic forces. This Lorentz force becomes significant at high altitudes and high speeds which are used to accelerate or decelerate a gas continuously at subsonic or supersonic speeds without choking even in constant cross section cha nnels. However, this is achieved by coupling the velocity and the temperature of the conduc ting media through Joule heating effects. In other words, the Joule heating from the electr omagnetic effects can si gnificantly alter the thermal equilibrium. Research on shock wave propagation in plasma has been studied by Rivir et al. [11] and high speed flow contro l using weakly ionized plasma by Leonov et al. [12]. The studies reveal that at high sp eed non-equilibrium flows the weakly ionized gases can substantially modify a traveling shock wave. A propagati ng shock wave in a plasma decreases in strength and the wave front disperses. Significant aerodynamic drag reduction via plasma injection in the stagnati on region has been reported in the studies of Rivir et al. The origin of this drag reduc tion is not completely understood, but the nonequilibrium thermodynamics and the electromagnetic forces play significant roles. The modification of the flow field struct ure due to non-equilibrium thermodynamic phenomena is caused in two ways. The energy from the weakly ionized plasma energy source modifies the local gas temperature a nd reduces the Mach num ber upstream of the bow shock. Secondly, the non-equilibrium ener gy distribution among internal degrees of

PAGE 21

6 excitation can also alter the thermodynamic pr operties within the shock layer. Also, in low supersonic flows, the high amounts of energy deposition upstream of the shock can make the flow locally subsonic. The modeli ng of the plasma effects in non-equilibrium effects also involves significan t complications as intrinsic plasma properties such as transport properties, conductivity etc. need to be calculated from kine tic modeling studies and not from the equilibrium Maxwellian dist ribution assumption used in fluid models. However, the focus of the present research is the low Mach number flow regime where non-equilibrium effects are minimal and, moreover, the plasma dynamics is mainly shaped by the electromagnetic effects and to a lesser extent by the background low-speed neutral flow. 1.3 Numerical Challenges Computational studies of the dielectric ba rrier discharges have been attempted for parallel plate arrangements with emphasis on modeling the plasma phys ics in detail [13]. The time and length scales operating in the pl asma extend over a vast range, warranting kinetic models in the extreme case (highly non-equilibrium, low pressure, high frequency (RF) discharge conditions) to relatively si mpler plasma-fluid models for equilibrium, high pressure plasmas. The latter is still an approximate model, but valid in certain situations such as at high pressures where equilibrium velocity di stribution functions and isotropy assumptions for pressure gradients can be used in spite of the large electrostatic effects. Even modeling the evolution of the pl asma discharge with reasonable complexity using a fluid model is computationall y challenging and expensive due to strong non-linearities between the plasma-f luid equations and the electrostatics widely varying spatial and temporal s cales, especially near the electrodes.

PAGE 22

7 Most of the computational work in the area of plasma modeling has been done for high frequency discharges and for parallel pl ate arrangements which can be extended to complex arrangements used in flow control studies. Keeping in mind the computational cost, 2-D or higher dimensional discharg e modeling studies are less common. More recently Roy et al.[14-18] have developed a finite-element method based 2-D modeling of such plasma-based actuators by solving for the discharge species a nd neutral transport. Although the model used only three species (ions electrons and ne utrals) in a 2-D framework, the computational time was signifi cant (a week of computational time for a typical actuator configuration) However, the plasma and the neutral fluid display threedimensional characteristics as reported from ex perimental observations. Also, the extent of plasma chemistry to be included is a compromise between exact physics and computational cost. Involving a large number of species will mean that we have to solve many transport equations and consequently a hugely coupled system However, we will see in later chapters that the impact of th e ion species (both posit ive and negative which are largely responsible for the collision dom inated momentum transfer) on the evolution of the discharge and electric force field is significant as deduced from the discharge structure. 1.4 Overview of Present Approach The key to modeling the DBD effects in fluid dynamics is to achieve realistic distributions of the species densities and their momentum in the domain which interacts with the neutral fluid by solving the plasma fl uid equations. As a first step to model these effects we developed a phenomenological mode l [9] based on a linear force distribution in the domain. This linearized body force mode l was later adopted by Gaitonde et al. [14] for modeling plasma-based separation control in a NACA 0015 wing section. This model

PAGE 23

8 was an attractive option because of the difficu lty in achieving efficient multi-dimensional and self-consistent plasma dynamics simula tions coupled with largely unsteady fluid dynamics. This was followed up by a first principles -based self-consistent modeling of the plasma-dynamics using a hydrodynamic model in helium gas. This multi-fluid formulation to model the radio-frequency disc harge in helium gas gives the spatial and temporal evolution of the charges species which is decoupled from the neutral fluid dynamics. The body force calculated from this data provides a more self-consistent way of modeling the plasma-fluid flow interaction. The wide range of operating length a nd time scales characterizing both the discharge and flow physics resulting in a highl y stiff system of equations and is a major limiting factor in achieving fully-coupled multi-dimensional plasma-fluid simulations. In this dissertation, we clarify some of the important computational issues pertaining to handling the multiple scales, while offering insight into improving our computational capability for modeling practical, multi-dimensional physics. Both sequential and fully implicit approach es [15] have been used for discharge modeling. In the present multiple scale problem, we have employed an operator-split sequential solution algorithm in which the solution procedure can be adapted to handle the individual processes efficiently and reali ze overall gain in computation. This method is integrated with a multi-block finite -volume algorithm capable of handling 3-D curvilinear grids [19]. The method is empl oyed to model the plasma dynamics in an asymmetric electrode configuration si milar to that shown in Figure 1.

PAGE 24

9 1.5 Outline of the Dissertation The remaining document is organized as follows. In chapter 2, we present a review of the various studies on discharge plasma to provide a reasonable understanding of their generation and opera ting mechanisms. The first section is devoted to experimental studies aimed at analyzing th e discharge operating characteristics and structure followed by some potential fluid fl ow control applications using the surface plasma arrangement. A brief glimpse is provi ded into the possible potential of these plasma actuators in a turbulent regime al ong with laminar separation control, drag reduction etc. The next section reviews th e various numerical modeling efforts of discharges operating in high and low frequency regimes. These give a view of the various versions of models employed to study the pl asma discharge in one and two dimensions. Further, we also look at some of the c ontemporary efforts at modeling the discharge plasma actuators in 2-D. In chapter 3, an overview of the various levels of plasma m odeling is presented starting from the studies which solve the Bo ltzmann equation in a kinetic model which are numerically intractable fo r large-scale problems. This is followed by a discussion on the hybrid methods which use both kinetic as well as fluid desc riptions with the goal of developing an efficient, feasible and accurate description of the processes. The last section presents in detail the various features of fluid mode ls including the various levels of complexity that can be achieved and the set of valid approximations. Motivated by our interest in efficiently handling the plasma equations, we present an overview of the techniques employed to handle such multiple scale processes. A solution framework for the coupled species transport equations and the electrostati cs is presented. The wide range of scales needs to be handled usi ng time-split methods and operator splitting

PAGE 25

10 techniques are employed for the source term treatment. With the large difference in the plasma and neutral fluid flow time scales the Navier-Stokes e quations are treated independently from the discharge model. Howe ver, the effect of th e glow discharge is included in the flow equations as a local body force term which can be evaluated from the quasi steady solution of the plasma-fluid equations. Having developed the modeling framework, a one-dimensional representative modeling study of a parallel plate glow discharge in helium is attempted as a valid ation excercise against prior modeling [20]. In chapter 4, a body force type coupling is developed in detail to facilitate interaction between the plasma dynamics a nd the fluid dynamics. A phenomenological body force model based on experimental obser vations, developed independent of the 1-D plasma model results, is used here to inves tigate qualitatively the induced flow structure. As a sequel to the above, the fully self-consis tent discharge dynamics developed from the hydrodynamic plasma model is presented to pr ovide insight into the discharge structure and the force generation mechanism. Chapter 5 summarizes the present research status and enunciates the future research goals to be attained.

PAGE 26

11 Figure 1 Experimental photograph illustrati ng the EHD effects on a laminar jet of titanium tetrachloride vapor. The top pi cture shows the unaffected jet when the plasma is non-existent and the botto m picture shows the effect on the jet which bends and forms a wall jet. Photograph reproduced from [4]. Glow discharge plasma Dielectric layer Electrodes Induced flow

PAGE 27

12 CHAPTER 2 LITERATURE REVIEW Flow control through various mechanisms has been extensively researched in literature. Recently, the application of surf ace glow discharge plasmas as flow control devices has attracted particular interest as re vealed in experimental studies undertaken in the last few years [5,21-27]. However, num erical studies to complement the above experimental efforts are not too common until recently [28-29,14-18]. The operating mechanism of the high frequency discharges be tween dielectric barriers has been studied using both experimental and numerical mode ling techniques, in the context of surface coating, sterilization, ozone production, plasma display panels (PDP’s) [3,6,30,31] etc. In the previous chapter we presented a brief overview of the dielectric barrier discharge plasma actuator, its operating mechanism, potential for thermo-fluid applications and some insight into the current status of the numerical aspect of the research and the associated challenges. Here, we review in de tail some of the rele vant experimental and numerical studies on discharge plasmas and fluid flows to help better understand the significance of this dissertati on. The first section 2.1 is de voted to experimental studies pertaining to discharge actuator mechanisms and applications. Th e second section 2.2 reviews the various plasma dynamics modeling efforts and concludes with the numerical studies attempted at modeling the discharge-induced fluid flows. 2.1 Review of Experimental Discharge Plasma Studies In this section we will review some of the experimental efforts to modify flow using the plasma based actuators and their operating mechanism. Initial studies of

PAGE 28

13 Dielectric barrier discharges, characterized the by the presence of numerous micro discharges which lead to filame ntation of atmospheric air. In these studies, the frequency of the applied voltage usually ranged from line frequency to a few hundred kilohertz. This was first investigated in the context of ozone production. A detailed review of various efforts undertaken is presented in [3]. 2.1.1 Early Discharge Studies It was Kanazawa et al.[5] in the late 1980s, whose work indicated that homogenous glow discharges as against filamentary ones can be obtain ed under certain conditions. Since then the potential of DBD’s in the contex t of plasma display pa nels for applications like flat television screens [3] etc. was reali zed, which led to spurt of research activity aimed at furthering the understanding of plasma dynamics applications. The experimental work on the plasma actuator development took off in the 1990s, when Roth et al.[32] formally devised the generation mechanism of uniform or homogeneous plasma discharges at atmospheric pressure which he called as the OAUGDP or One Atmosphere Uniform Glow Discharge Plasma. The OAUGDP is capable of operating at atmospheric pressure in air and other gases without significant generation of heat. It is nonequilibrium plasma, far from thermodynamic equilibrium. This means that the individual species such as the electron possess energies which are very far apart from the heavier ions as well as the neutral species. The generated plasma is more uniform in the breakdown region and attractive for surface and large volume app lications. The active species thus generated had potential for use in sterilization, surface decontamination, increasing surface wettability and surface energy of materials [9]. A schematic of the OAUGDP generator is given in Figure 2. This plasma generation devi ce can be operated

PAGE 29

14 in a wide range of geometrical configurations such as, between two parallel plates or as surface plasma which can cover a flat or curved surface such as an aerofoil. The asymmetric arrangement in Figure 2b is the source of the surface plasma over an aerofoil, with the direction of the ge nerated electro hydrodynami c (EHD) flow given by the placement of the submerged electrode, i.e. the flow is always seen to be generated from the upstream electrode to the downstream electrode. This E HD flow is caused by the ‘paraelectric’ effect [4] the term given to the flow (say 1-10 m/s) created because of the curved electric field lines as in an asymmetric arrangement shown above in Figure 2c. Only an asymmetric arrangement as shown in Figure 2b generates us able wall jet type flow. There are other symmetric arrangemen ts studied in [4] which by virtue of construction produce vortex flow as against wall jets. The mome ntum acquired by the ions in such an electric fiel d is coupled to the neutral fluid through Lo rentzian collisions. Roth et al.[33] identifies use of these EHD flows or its variants as control tools in aerodynamic boundary layers, pumping gases th rough tubes or ducts etc. A few key observations of the OAUGDP from th ese experimental studies [33] are Plasma is non-equilibrium and cold with time-dependent features. Discharge exhibits classical glow discharge characteristic s with identifiable regions as identified by computer modeling in [20,34] Ionization occurs at ‘Stoletow’ po int with minimum possible energy. By non-equilibrium, it is implied that the plasma does not have thermodynamic equilibrium and the kinetic temperatures of the species are far from equal to each other. This further implies that the background ga s need not be heated till thermodynamic equilibrium. Three prominent EHD based flow effects discussed in [4] are presented in Table 1.

PAGE 30

15 2.1.2 Mechanism of Plasma-I nduced Flow Applications In this section we will present efforts to address the paraelectric flow mechanism. Experimental and analytical studies were conducted by Roth et al.[4,7,33]. Here, EHD body force is modeled as the electrostatic fo rce acting on the charge d particles which acts on the neutral gas. The genera tion of paraelectric effects as a consequence of electric field gradients is illustrated in [32-33]. A brief description of the analytical argument of Roth et al.[33] is given below. If E, is the electric field (V/m) and c (C/cu.m) is the charge density, then the local electric body force FE (N) according to Roth et al.[7] is ECFE (N) (1.1) Using Maxwell’s equation, we get 0.CE (1.2) Hence, 2 00011 .. 22Ed FEEEEE dx (1.3) for a one-dimensional case. This formulation for the body force allows viewing it as an electrostatic pressure gradient where, 2 01 2EPE. (1.4) This simplified 1-D interpretation gives a feel for the nature of paraelectric body force effect, but cannot be ex tended to higher dimensions. Mo reover, evaluating the exact distribution of E, requires solving for the char ge particle distributi on. Nevertheless, the above electrostatic pressure argument give s a feel for the observed phenomenon. The ‘paraelectric’ effect can be interpreted as the pumping of charged particles towards

PAGE 31

16 regions of increasing electric field gradient. Also, in a quiescent medium, the electrostatic pressure acts as a driving pressure gradie nt causing the flow from a region of higher pressure to that of a lower pressure as repo rted from the smoke fl ow visualization study in [4]. More recently, Enloe et al.[21-22] studied, the plasma morphology and operating mechanism in an asymmetric electrode setup as shown in Figure 2b by means of optical measurements using a Photomultiplier (PMT) tube. The experimental study revealed characte ristic temporal and spatial structure. Computer modeling efforts for parallel plate arrangements in [20, 34] showed that the DBD generated in those conditions possessed ch aracteristic low pressure glow discharge qualities with regions such as the Faraday dark space, cathode glow, positive column etc. This is illustrated in the 1D modeling study by Lee at al [34] as shown in Figure 3. The figure shows the distribution of the electric fi eld and species densitie s at the instant of peak current which are characteristic of glow discharges. The temporal structure of the DBD generated in an asymmetric arrangement indicated noticeably different discharge patterns in the two half cycles. The light inte nsity used as a measure of the discharge was consistently regular when the exposed elect rode was negative-going and highly irregular when positive-going. Further, it was noticed that a consistently increasing voltage form was necessary to sustain the discharge and th e resulting light intensity depends heavily on the rate of increase. Also to study the asym metry effects further, asymmetric saw-tooth waveforms having vastly different positive-going and negative-going duty cycles were used to model the discharge. The resulting outcome clearly established the preferential phase for the actuator operati on during a whole time-cycle.

PAGE 32

17 (a) parallel plate setup (b) A single asymmetric electrode set up and its effect on fluid flow. (c) Illustration of distorted field lines in an asymmetric arrangement Figure 2 Generation of glow discharge in vari ous arrangements such as between parallel plates (a) or as a surface discharge (b) and the curved electric field lines in such a case [9]. outgoing wall jet AC voltage Source ( 5-10 KV, 1-20KHz [4] ) insulated electrode substrate insulatio n exposed electrode p lasma Incoming flow Flow bends Upper electrode Dielectric Lower electrode Electric field lines DIELECTRIC HIGH VOLTAGE POWER SOURCE ELECTRODES

PAGE 33

18 Table 1 Representative EHD operating mechanisms Type of EHD effect Max speed Mechanism Comments Paraelectric body force effect due to net charge density in the plasma 1-10 m/s Electric field gradients Very large collision frequencies (7GHz) ensure adequate momentum by Lorentzian collisions. EHD convection driven by DC ion mobility drift 300 m/s DC electric field over a generated OAUGDP causes ion drift Impractical due to high electric fields of the order of 104V/cm requirement. Peristaltic EHD plasma acceleration 100 m/s Series of paraelectric actuators operating at a phase difference Gives a higher induced flow velocity to work withmost useful and flexible of the three Note: See [33] for more information. Figure 3 1-D Numerical modeling results from L ee et al.[34] showing different regions of the discharge gap at the instant when the current is maximum for a 10KHz, 1.5KV voltage. acathode fall; bne gative glow; cfaraday dark space; dpositive column While the light intensity measurements of Enloe et al.[21-22], revealed the significant temporal and spatial structure of the discharge in the actuator, it still did not improve the understanding of the nature of momentum coupling. In a follow-up of the above study, Baird et al.[35] investigated the acoustic effects of the plasma actuator. The study revealed that the momentum transfer to the neutral fluid was different in the two halves and the acoustic pattern produced was a radiation pattern characteristic of a

PAGE 34

19 coherently-driven system indicating possible comp ressibility effects. Near electrode fluid density measurements over time by Enloe et al.[36], indicate density variation in space akin to a pumping-type mechan ism resulting in a wall-jet. However, hot-wire temperature and velocity measurements by Jukes et al.[37-38] around the exposed el ectrode show only a slight increase in ambient temperatur e of two degrees, thus reducing the possibility of buoyancy driven effects. In the present dissertation we proposed as part of our phenomenological model (Shyy et al.[9]) that, the asymmetries intr oduced through the electrode arrangement and the dielectrics along with the asymmetric disc harge structure had a si gnificant role to play in the generation of EHD based paraelectric fl ow or simply glow-discharge induced flow. It can be seen from Figure 4, that the distinct difference in the light intensity measurements between the two half cycles for asymmetric electrode arrangement. The first half cycle provides a more uniform disc harge of high intensity while the other half cycle provides an overall weak light intens ity. Overall, the studies on the discharge actuator summarily observed the following: Plasma structure is different in the tw o half cycles which was a fundamental assumption in the development of our phe nomenological model (Shyy et al.[9]) and was also suggested in [21] (see Figure 4) It is possible caused by the geometric asymmetry. Material asymmetry of the exposed and s ubmerged electrodes only influenced the net plasma intensity and a weak indu ced flow was present even with both electrodes insulated. Direction of flow is purely determined by electrode geometry for the actuator configuration discussed above. Mechanism of momentum coupling is yet to be adequately understood with possible reasons ranging from densit y-variation effects from acoustic measurements to contradicting temperature measurements indicating minimal buoyancy related effects.

PAGE 35

20 It was also perceived from the outcome of our phenomenological model [9] that the net induced flow was a product of the vari ous competing mechanisms listed above. Further, it also highlighted the glaring n ecessity for detailed multidimensional modeling of the actuator and plasma dynamics towards developing an accurate modeling strategy. Parametric studies and the impact of th e geometry on the actuator performance were studied by Enloe et al.[22] and Van Dyke n et al.[26]. Specifically, the impact of voltage frequency, waveforms and geometrical parameters were investigated. Larger thickness of the exposed electrode resulted in more thrust which indicated possible ion concentration near the exposed electrode as they are the primary momentum imparting collision species. Also, a larger thickness can mean a sharper electrode edge while transitioning from the conductor to the insulator surface which can lead to stronger electric fields in the region. The charge species concentration in the plasma will lead to significant Debye shielding effect s which can affect el ectric field distribu tion. In trying to provide an analytical insight, Enloe et al, discusses the eff ect of Debye shielding which further illustrates the inadequacy of Ro th’s electrostatic pressure argument. 2.1.3 Aerodynamic and Control Strategies Using Glow Discharges Now that we have looked into some of the efforts to understand the plasma actuator operating mechanism, we will look at it s various areas of application. The above discussed potential of the surface plasma to generate EHD flow effects whether wall jet type flows or re-circulating flows can be ha rnessed for a variety of applications. The injection of momentum into the boundary layer and delaying of separation has been primary motives in aerodynamic surface flow control.

PAGE 36

21 Stronger half cycle Weaker half cycle Stronger half cycle Weaker half cycle Figure 4 A 3-D surface plot of photomultiplier light output from the plasma as a function of time and chordwise distance, for a si ngle AC cycle of the applied voltage. (Reproduced from Enloe et al.[23]) However, most of the efforts [4,25,39-40] listed indicate effective performance of the plasma actuator from low (10,000) to moderate (100000) Reynolds numbers, essentially limited by energy transfer from th e plasma. Van Dyken et al.[26] use a single plasma actuator at 7% chord on a NACA0015 airf oil resulting in subs tantial increases in lift and reduction in drag beyond the baseline stall angle of attack. Santhanakrishnan et al.[41] present the pos sibility of using plasma-based flow control for Unmanned air vehicles (UAV’s) to control st all and enhance lift, especially in the range of low to moderate Reynolds number range(104-105). Goskel et al.[42] did an experimental investigation of steady and uns teady plasma wall jets for se paration control of Micro-airvehicles in the low Reynolds number regime (~ 20,000). It was observed that a gain in lift by a factor of two and the loss of low Re ynolds number hysteresis was achieved using discharge power inputs in the range of 1.2mW/ cm. Corke et al.[43], presents the idea of using an unsteady plasma actuator to achieve flow control performance on an airfoil comparable to that of a moving plane flap. Th is application is interesting for the reason that the unsteady operation was observed to re duce the actuator power consumption up to

PAGE 37

22 ninety percent. Related work includes stud ies by Corke et al.[24], Rivir et al.[11] and Leonov et al.[12]. Table 2 Summary of experimental plasma actuator studies Author and year Work description Comments Roth et al(1998) [4] Boundary layer flow control using plasma actuator arrays Increase in drag. EHD flow altered boundary layer profile. Corke et al(2000) [24] Phased plasma arrays for unsteady flow control Corke et al.(2002) [25] Flow control on a NACA 0009 wing Lift enhancement and also that of drag (higher skin friction) when typically using single actuators. Post et al.(2003) [27] Separation control on a NACA 0066 wing Increased lift to drag ratio and delays separation. Roth et al.(2003) [39] Flow re-attachment using Paraelectric and peristaltic EHD effects on a NACA0015 wing Peristaltic acceleration is by a traveling electrostatic wave List et al.(2003) [40] Laminar separation control on a cascade turbine blade Separation is quelled at moderate Re~50000 Enloe et al.(2004) [21-23] Experimental visualization of the generated plasma structure in space and time, effect of geometric construction Plasma structure helped gain insight into the operating mechanism Rivir et al, (2003) [45] Turbine flow control with plasma Rivir et al.(2004) [11] AC and Pulsed plasma flow control AC plasma controls flow through wall jet mechanism and pulsed DC Plasma through heating mechanism Leonov et al.(2003) [12] Effect of plasma induced separation Uses plasma as a local heat source Chan (2005) [44] Flow control using acoustic effects A symmetric plasma actuator generates vertical structures interacting with shear layer instabilities. Enloe et al.(2005) [36] Temporal density, velocity measurements near the plasma actuator Indicated density variation slope near the electrode edge as possible reason for momentum coupling. Baird et al.(2005)[35] Acoustic testing of plasma actuator. Acoustic patterns showed a characteristic coherently driven system indicating compressibility effects. Jukes et al.(2006) [37] Characterization of dischargeinduced wall fo w mechanism The discharge actuator characteristics are investigated using hot-wire measurements in the near wall region.

PAGE 38

23 In [24] a plasma array with phased inputs with unsteady flow control is investigated. Separation control on a NACA66 wi ng is studied by Post et al.[27] and is illustrated in Figure 5. The plasma operation significantly reduces the size of the separation bubble. In [12] th e plasma generation is used as a local heating mechanism thereby inducing boundary layer separation. Ri vir et al.[11] compare the paraelectric surface flow control mechanism with the heating mechanism of Leonov et al.[12]. A summary of experimental effort s is presented in Table 2. Figure 5 Separation control on a NACA66-18 wi ng at a Re=79000, 16 degrees angle of attack, reproduced from Post et al.[17] A novel approach to apply plasma based EHD flow control to subsonic flows is being investigated by Chan [44]. The acousti c properties of low rectangular cavities are being investigated for the potential use of plasma as flow-control devices. In the

PAGE 39

24 experiments, the electrodes of the plasma actuators were laid out streamwise to the flow. A symmetric electrode arrangement (where the lower electrode is wider compared to the exposed upper electrode and placed exactly below it to ensure symmetry) which can generate vortical structures is used upstream of the cavity to inhibit the feedback process in the shear layer which is key to the generation of large amplitude frequency tones. The plasma actuators were reported to attenuate co mpletely the discrete tones produced by the cavity. The acoustic effects can act on the time-scale smaller than the fluctuating frequency. There have also been efforts to use pl asma based actuators in the context of turbulent boundary layer flow control. EHD based turbulent flow control concepts has been presented in recent articles by Soldati et al.[45-47]. As in the present effort, the electrohydrodynamic effects are modeled as a variable body force term in the equations describing fluid motion. The complex interac tion between the EHD effects, wall structure of turbulence and the mean flow were studi ed. The EHD flows affect the turbulence field by increasing both dissipation as well as production thus maintaining the balance. However, overall drag reduction was observed with consistent decrease in the Reynolds stresses. The ionization in thes e cases is through elec trostatic precipitato rs and efforts are being made to use OAUGDP in this context. W ilkinson [48] investigates the use of an oscillating surface plasma wave for turbulent dr ag reduction. Jukes et al.[38] looks at the turbulent boundary layer control for skin fr iction drag reduction us ing the surface plasma actuator. A summary of turbulent flow control literature relevant to the present study is tabulated in Table 3.

PAGE 40

25 Table 3 Summary of turbulent flow control studies Author and year Work description Soldati et al, (2001) [47] Turbulent boundary layer control from EHD flows. Drag reduction is achieved by modifying the Reynolds stress contributions Soldati et al, (1998) [45] Turbulence control using large scale EHD flows Du et al, (2002) [49] Drag reduction in wall-bounded turbulence using a transverse traveling wave Dhanak et al, (1999) [50] Reducing turbulent wall friction through spanwise wall oscillations Wilkinson (2003) [48] Turbulent drag reduction using oscillating surface plasma 40 % drag reduction was achieved, but ineffective at frequencies < 100 Hz Jukes et al.(2004) [38] Turbulent boundary layer control for drag reduction. 2.2 Review of Discharge Modeling Studies In general, the numerical modeling effo rts are developed in conjunction with experimental studies to determine their fidelity. The complex physiochemical and spatiotemporal processes associ ated with the glow discharge plasma also make numerical modeling essential for their understanding and in modeling applications of plasma-based devices. Experimental studies as reviewed in section 2.1 were limited in their ability to reveal plasma structure and adequately ex plain the observed plasma dynamics and other mechanisms. Specifically, the apparently contradicting arguments for the momentum coupling based on compressibility effects prop osed by Enloe et al.[35] as against the observed lack of significant bulk heating due to the plasma reported by Jukes et al.[38], highlight the necessity for a high-fidelity first principles-based numerical modeling approach. Also, most of the preliminary mo deling attempts were empirical [5,22] and argumentative in nature. In order to establis h a competent mechanism for the paraelectric

PAGE 41

26 and other induced flow effects observed, de tailed simulation of the plasma dynamics in an asymmetric electrode arrangement is necessa ry. In this section we will review some of the plasma and plasma actuator modeling e fforts to establish the background for the present dissertation. Section 2.2.1 consists of the various phenomenological/simplified modeling approaches employed by various inve stigators to model the discharge-induced fluid dynamics while, in section 2.2.2 we wi ll present some of th e more high-fidelity discharge modeling studies. 2.2.1 Phenomenological Modeling of the Plasma Actuator Most of the preliminary plasma actuator modeling efforts to study the fluid flow coupling was either simplified phenomenologica l approaches or models with reduced complexity. In the present dissertation we develop our own phenomenological model to study the plasma-induced fluid flows. This model involved the development of an empirical approach to calcul ate the time-averaged force to study the momentum coupling with the neutral fluid. The various discha rge parameters were assumed spatial or temporal variations with the closure of the problem related constants requiring experimental data for validation. Our analytic al-empirical model (see Chapter 4) based on a linear force distribution in the domain wa s later adopted by Gaitonde et al. [14] for modeling plasma-based separati on control in a NACA 0015 wi ng section. This linearized body force model was an attractive option beca use of the difficulty in achieving efficient multi-dimensional and self-consistent plasma dynamics simulations coupled with the fluid dynamics. Hall et al.[51] use a poten tial flow approach to mode l the aerodynamic effects of the plasma actuator. In particular, the ac tuator mechanism is modeled as a doublet (source-sink) to mimic the pumping action of the discharge on the neutral fluid.

PAGE 42

27 Experimental validation was again required as a closure to determine the doublet strength. In a different approach, Orlov et al.[52], employ a lumped-element circuit model to account for the spatial and temporal variation of the discharge and mimic the plasma body force effect on the neutral fluid. While the above studies showed qualitative and reasonable qualitative corres pondence with the experimental observations, they were limited in their ability to confidently predict as well as explain the observed mechanisms and the plasma-fluid interaction. 2.2.2 First Principles-Based Plasma Modeling Approaches Most of these studies were undertaken eith er as part of studying dielectric barrier discharges (DBD’s), between parallel plates operating in the kilohertz regime (low frequency) or the capacitively coupled radio frequency (high frequency) discharges for plasma applications with special emphasis on the plasma actuator. The knowledge of parallel plate discharge modeling is esse ntial for gaining insight and developing a numerical model for the plasma actuator. This is divided in to two subdivisionsthe first one review the numerical challenges and approaches in the modeling of DBD’s in general, while the latter places sole emphasi s on the discharge-coupled fluid flows which are directly relevant to the present dissertation. 2.2.2.1 Review of dielectric barr ier discharge (DBD) modeling studies The modeling of low frequency (KHz) barri er discharges was attempted in the nineties, when the potential of uniform glow discharges was realized and experimental studies related to flow cont rol and other surface processe s gained momentum. In the various studies reviewed here, the focus is on parallel plate gl ow discharges controlled by a dielectric barrier in gases such helium, neon, argon, nitrogen and even air [53]. By low

PAGE 43

28 frequency we mean a range 50Hz to a few KHz which is much smaller than the radio frequency regimes (MHz). The modeling of glow discha rges involves the considerat ion of mechanisms such as the reaction kinetics and plasma dynami cs. With the high degree of complexity involved, the choice of mode l is usually simplified tailor ing to the application. An overview of plasma models of varying complex ity will be presented in the next chapter. The approximation chosen in different studies are mostly based on choice and convenience, resulting in each model differing sl ightly from the other while retaining the original qualitative behavior. One of the firs t numerical studies of dielectric barrier controlled discharges at atmospheric pressure was studied by Eliasson et al.[54] using a two-dimensional fluid model for the plasma. Almost all modeling studies employ a fluid model as against solving a Boltzmann equation b ecause of the relatively easier nature of solution technique for the former. The fluid m odel consists of the tr ansport equations for the charged and neutral species concentrati on, momentum and energy. The first three moments of the Boltzmann equation are considered for modeling purposes. The momentum and energy equations are also solved depending on the complexity of the model. The numerical solution consists of thr ee steps which all are interconnected. The first step is the solution of the Boltzma nn equation to obtain the electron energy distribution function in the gas under investigation. This is e ssential, for the properties of electron transport such as mobility; diffusivity mean energy etc vary as a function of E/N (E is electric field and N is the gas number density). The above transport properties are obtained solving a steady state Boltzmann e quation under a constant field condition. This

PAGE 44

29 assumption is reasonable because the characte ristic time for energy equilibration under a changing electric field at at mospheric pressure is 10-11s while the characteristic time for electric field fluctuation is 10-9s [55]. This makes it possible to assume energy equilibration to be an instantaneous proce ss. The solution of the Boltzmann equation for the electron energy distribution and other trans port properties has be en tabulated over a wide range of E/N and can be used during th e next phase which is the transport module. In this second step, the transport equations of species number density, momentum and energy are solved. The electric field E is solved from the Poisson equation 0. E (1.5) where is the net charge density. Also, the charge build up on the dielectric as the discharge evolves and the external circuit s hould be added to the model at this stage. Discussed above is brief plas ma modeling framework using fl uid models. Eliasson et al, study a two dimensional model of the dielectric barrier discharge in three gasesXenon, nitrogen and oxygen. Most mode ling studies use iner t gases in their st udies so as to obtain stable glow discharges. Golubovskii et al, [13] use a multi-spe cies 1-D fluid model to study the different modes of operation of a helium discharge at atmospheric pressure. The dielectric barrier discharg es were seen to operate in two modesa glow discharge mode, displaying characteristic low pressure discharge st ructure and a Townsend mode characterized by multiple current or discharge peaks. The formation of the various modes is governed by design parameters such as ga p width, dielectric barrier thickness and applied voltage. Specifically, thin barriers and wide gaps led to a glow discharge mode and a Townsend mode was realized otherwis e. Similar behavior was reported in the

PAGE 45

30 experimental studies by Mangolini et al, [56]. Two-dimensional model studies were performed by Golubovksii et al., [57] Massines et al., [58] a nd Xu et al., [59]. All these studies included substantial plasma chemistry to account for the exact species concentrations. However, the significant aspect of these studies was it revealed significant two-dimensional structur e for the discharge. Xu et al [59] state that the radial expansion for the discharges can significantly alter the built up charge distribution on the dielectric surface. This has been observed in a number of two-dimensional studies such as Xu et al.[59] and Golubovskii et al, [60]. Mangolini et al.[56] state that the multiple pulses as seen in a Townsend mode are manifestations of the radial inhomogeneities. The radial expansion of the current pulse cau ses a non-uniform surface charge distribution with the next breakdown occurring at th e discharge periphery Another important modeling aspect is the consideration plasma chemistry. The treatment of plasma chemistry varies from the most complex reaction mechanisms which give rise to conservation and handling of multiple ion species, metastable atoms or molecules. The Penning ionization effect also contributes to the complex plasma chemistry. The simplest treatment of plasma chemistry is the two sp ecies model as used by Ben Gadri et al.[20] Shin et al,[61] with more complex mu lti species treatments by Mangolini et al.[62] Massines et al.[6,58] Tochi kubo et al.[63] etc. Mangolini et al.[62] state that the metastable atoms in helium discharge contribute significantly to the pre-ionization through penning ionization mechanisms. This pre-ionization is important for the next breakdown. Handling more number of species im plies higher computational cost and is usual kept at a possible minimum. The choice of boundary condition for the transport equations differs in various studies which try to approximate the species behavior in the

PAGE 46

31 sheaths. Golubovskii et al.[60] discusses the nature of interaction between the charged particles and the dielectric boundary surface in the context of a homogeneous barrier discharge in nitrogen. The infl uence of different mechanisms of electron emission and the recombination are studied. A detailed discus sion of the different numerical modeling strategies, treatment of boundary conditions etc. are presented in later chapters. A summary of the different low frequency numer ical studies are tabulated in Table 4. Table 4 Summary of low frequency plasma modeling studies Author and year Work description Model Golubovskii et al.(2003) [13] Modeling of dielectric barrier discharges in helium at atmospheric pressure Multi-species 1-D fluid model. Investigated various modes of barrier discharge operation Golubovskii et al.(2004) [56] Modeling of nitrogen Multi-species 2-D model Dongsoo Lee at al, (2004) [34] Modeling of DBD’s using a helium-oxygen mixture 1-D multi-species fluid model Massines et al, (1998) [6] DBD modeling using He-Ar mixture 1-D two-species fluid model Massines et al, (2000) [58] DBD modeling using He-Ar mixture 1-D multi-species model Tochikubo et al, (1999) [63] DBD modeling of pure Helium 1-D multi-species model Mangolini et al, (2004) [62] Modeling of Heliumoxygen mixture 1-D multi-Species model Shin et al, (2003) [61] Modeling of He-Ar mixture 0-D, 2-species model Xu et al, (1998) [59] Modeling of nitrogen based mixture 2-D multi species model Eliasson et al, (1991) [54] Modeling in Xenon 2-D multi-species model High frequency modeling studies have been the focus of numerical model development for plasma simulations, in view of their extensive applications in material processing, plasma display panels (PDP’s) etc. While the nature and physics of these RF discharges are different from what we are interested, the study of these discharges is useful for their numerical modeling. Prior to the last decade, the simulations were basically one-dimensional either with a flui d model or a particle model of the plasma. Since then, two and three-dimensional studies have been successfully presented. A

PAGE 47

32 review of one-dimensional studies is pres ented by Govindan et al.[64]. The dominant type of plasma considered in these studies are low pressure, low-temperature and weakly ionized. High pressure glow disc harges have also been modele d lately by Yuan et al.[65]. The fluid model is based on the self-consistent solution of th e electron and ion continuity and momentum transfer equations coupled with the Poisson equation for solving the potential between the electrodes. Boeuf et al .[66] assume a local e quilibrium between the charge particle kinetic and the electric field. This local e quilibrium hypothesis serves as a closure relation in representi ng the actual Boltzmann equati on through the continuity and the momentum transfer equations. Although, nu merical modeling efforts to solve for the Boltzmann equation have been at tempted [67] they are less preferred to simpler fluid models. Comparison studies have been made between the particle and fluid simulations by Nitschke et al.[68] to study the range of validity of the fluid approximation for the plasma. At higher pressures, the fluid model and the local equilibrium turn out to be pretty good approximations for the RF plasma. More detailed discussion of the above will be presented in later chapters. However, fluid models of varying complexity have been summarized in Table 5. Solution of the transport equations in the fluid models is computationally intensive because of their very stiff nature in space and time. To illustrate this we will look at the streamer propagation modeling in a DC voltage by Dhali et al.[69]. The propagation of the streamer obtained by a 2-D modeling study is shown in Figure 6. The wave-like solution and sharp gradients are characteristic of streamer propagation and also low frequency discharges The system of equations is also highly non-linear with the modeling of multiple sp ecies transport equations along with the

PAGE 48

33 Poisson equation for the electric field. The va rious numerical techniques used to handle these issues will be reviewed in the following chapters. Table 5 Summary of high fre quency plasma modeling studies Author and year Work description Comments Boeuf et al, (1987) [66] 2-D fluid model 2 moment equations, two specie Nitschke et al, (1994) [68] 2-D fluid and particle model comparison 3 moment equations are used in fluid model and two species Lymberopoulos et al, (1995) [70] 2-D fluid model for Gaseous electronic conference (GEC) reactor Multi species plasma chemistry and three moment fluid formulation Campbell et al, (1995) [30] 2-D fluid model for a PDP Multi species plasma chemistry and two moment fluid equations Boeuf et al, (1995) [31] 2-D fluid model for Gaseous electronic conference (GEC) reactor comparison with experiments Two species plasma chemistry and three moment fluid equations Meunier et al, (1995) [71] Veerasingham et al, (1996) [72] 1-D fluid model for PDP’s Multi species plasma chemistry and two moment fluid equations Yuan et al, (2003) [65] 1-D fluid model for high pressure helium glow discharges Multi species plasma chemistry and three moment fluid equations 2.2.2.2 Review of plasma actuator modeling studies While numerical studies on discharge phys ics have been popular, computational modeling of the discharge actuator and its physics have been attempted only recently [1418,29,74] in parallel with our ongoing effort s. The modeling the DBD effects in fluid flows requires realistic distributions of the species densities and their momentum in the domain which interacts with the neutral flui d by solving the plasma fluid equations. This will in turn enable us to calculate the net electric force acting on th e charge species and acts as a body force on the neutral fluid. More recently, Roy et al.[15-16] proposed a

PAGE 49

34 self-consistent two-dimensional DBD fluid mo del for helium gas with application to separation control using finite element techni ques. This multi-fluid formulation to model the radio-frequency discharge in helium gas gives the spatial and temporal evolution of the charges species which is decoupled fr om the neutral fluid dynamics. The body force calculated from this data provides a more se lf-consistent way of modeling the plasmawall jet interaction. This study employed a finite-element method and solved the various species transport equations using a globally implicit procedure wh ere the system of charge continuity and momentum equations are assembled as part of a global matrix to solve for the solution vector. Singh et al.[17], in a related paper, present a parametric study of the different conditions in an asymme tric discharge configuration. Specifically, it was observed that the net body force producti on in the domain over a whole time-cycle produced a downward positive force for a typical configuration. Ku mar et al.[18] study the nature of the discharge and the resulting force field in the presence of a magnetic field. Also studied is the shape effect of th e electrodes in the event of a finite electrode thickness. All the above-men tioned modeling studies assume negligible electrode thickness. In the event of a finite-electrode thickness, the treatment of the dielectricelectrode edge can impact the near-wall force field. The common approach to be employed in the present study as well as in the previously mentioned modeling studies is that the plasma species and the neutral flui d are treated as a two-fluid system coupled through dynamic forces and pressure interac tions. More recently, the simulation of the asymmetric electrode discharge actuator has been extended to atmospheric air chemistry instead of the Helium gas [75]. Visbal et al.[ 76] explored the potenti al of these plasma actuators in the contro l of turbulent and transitional se parated flows by means of the

PAGE 50

35 discharged-induced near wall effects. Spec ifically, a pulsed act uator producing an unsteady discharge was employed to highlight the role of turbulen t and transition causing mechanisms in drag reduction as against a pure wall-jet momentum injection. A summary of the various modeling studies ha ve been listed in Table 6. Figure 6 Streamer front propa gation at different instants in a DC voltage (Dhali et al.[69]) Table 6 Summary of plasma-induced fluid flow modeling studies Author and year Work description Model Roy et al.(2005) [15-16] Self-consistent discharge modeling using hydrodynamic plasma model. Employed a 3-species modelions, electrons and neutrals for Helium discharge. Singh et al, (2005) [17] Parametric analysis of different discharge operating parameters and geometric effects Same as above. Kumar et al, (2005) [18] Discharge and fluid flow evolution in the present of a magnetic field and shape effects of the electrode edge Same model as above. Singh et al, (2006) [75] Discharge modeling extended to air chemistry and the corresponding flows effects were studied. 8-species model with atmospheric air chemistry Visbal et al, (2006) [76] Control of transitional and turbulent separated flows using phenomenological and first principles based models. Electron number density (per cu.cm) Axial distance (0-5mm) 0.1ns 1ns 2ns 2.5ns 3 ns

PAGE 51

36 CHAPTER 3 NUMERICAL MODELING OF GLOW DISCHARGE PLASMA In chapter 2 we reviewed some of the relevant experimental and numerical efforts in the modeling of the discharge actuator e ffects for fluid flow applications. We also highlighted some of the numerical issues arising from the multi-scale nature of the problem and how it can limit the computational capability. The presence of these large ranges of scales operating in the problem necessitates the use of modeling strategies which can range from kinetic level models to the more common hydrodynamic/fluid type models. In this chapter, we present a glim pse into some of these different modeling approaches and in the process, present an argument for the choice ‘fluid model’ used to describe the plasma dynamics in the present dissertation. In the second half of this chapter, we present in detail, the various num erical issues associated with fluid model and the presently developed solution algorith m which can efficiently handle the plasmafluid dynamics. Lastly, a one-dimensional pl asma dynamics simulation to model the discharge evolution between parallel plates is presented for validating the present capability. 3.1 Plasma Modeling Hierarchy Modeling of glow discharge plasmas has been researched considerably more in the context high frequency (radio frequency (M Hz) or KHz range) discharges than low frequency (usually smaller than KHz range) atmospheric pressure discharges. However, the studies of radio frequency (RF) discharg es do provide insight for handling the stiff, multi-scale and highly non-linear plasma equati ons coupled with the electrostatics. The

PAGE 52

37 fundamental numerical model describing th e high and low-freque ncy discharges is essentially the same with the lone difference coming from the operating frequency of the driving voltage. While this may impact the natu re of the discharge, it retains the same inherent physical characteristics in term s of length and time scales. The following discussion is equally applicable to bot h high and low frequency modeling studies although most of the work conducted in this dissertation focused on efficiently modeling discharge plasma actuators in the KHz regime. Computational modeling of discharge plas mas can be achieved using a variety of approaches. In general, we can have a fluid treatment a kinetic treatment a hybrid scheme which combines aspect s of the fluid and kinetic schemes. In this chapter we will deal with the differe nt types and levels of plasma modeling. The hierarchy of models based on complex ity will be in three areas namely, type of model (fluid or kinetic etc.) treatment of plasma chemistry space complexity (1D, 2-D). The various levels of plasma modeling are illustrated in Figure 7. The chapter is organized as follows. A brief overview of th e three broad classifications of plasma models will be presented with emphasis on the fluid model which is of interest in the present work. The modeling complexity aris ing from plasma chemistry and spatial treatment will be discussed in the context of the fluid model.

PAGE 53

38 Figure 7 Levels of plasma modeling 3.1.1 Kinetic Modeling Approach In this section we will give brief overview of the kinetic models and how they compare with the fluid model employed in th is work. Kinetic models typically involve the solution of the Boltzmann equation for the species velocity or energy distribution function in both space and time or the particle simulations using techniques such as Monte Carlo methods. Although the Boltzm ann equation can provide much more information than the continuum fluid models its solution is computationally intensive requiring an order of magnitude more computer time than the fluid calculation as kinetic models are inherently higher dimensional and the solutions are functions of the velocity, space and time, the additional degrees of freedom being described in the full phase space Glow discharge modeling Kinetic model Fluid model (valid at high pressures) Mass conservation Momentum conservation Energyconservation Boltzmann equation solving for velocity distribution function Particle models Involving classical equations of motion for the particles in a force field Monte-Carlo collision model (MCC) Based on stochastic or probabilistic modeling Particle-in-cell (PIC) technique Based on deterministic classical mechanics Hybrid model Combination of the fluid model along with one of the kinetic models.

PAGE 54

39 of the problem. For example, Riley et al.[ 77] use a 1D-2V formulation which implies solving for the particle velocity distribution function (VDF) in a single spatial direction and in two directions of the velocity space. Th is is in contrast to the fluid models where we assume an equilibrium distri bution function for the species. The coupling between the Boltzmann equation and the Maxwell’s equations for electrostatics is a difficult num erical problem as the partic le dynamics should include the strong body force effect of elec trostatics. The solution variable in the kinetic models is the species velocity distribution function. By assuming an equilibrium distribution function in the case of fluid models, the full Boltzmann equation can replaced by its first few moments corresponding to the conservatio n of mass, momentum and energy. In a continuum regime the plasma modeling st arts from the Boltzmann equation [55] ../ii iiicii collff cfFfm tt (2.1) The index ‘i’ refers to single species. Here f is the distribution function, c the charged particle velocity, F the force, the gradient operator in physical space, c the gradient operator in velocity space. The subscript ‘coll’ refers to the change of distribution due to collisions. F represents the force vector which can be the Lorentz force on the charge particles. A detailed discussion on the theory of Boltzmann equation and its relation to the fluid models for radio frequency glow disc harges is presented in [55]. The index ‘i’ implies that (2.1) will be written for every ki nd of species. To enable a self-consistent solution of the above, we need to solve for the electric field E by means of the Poisson equation. 0. E (2.2)

PAGE 55

40 Table 7 Kinetic model studies of discharge plasmas Author and year Work description Kushner [78], 1983 Monte-Carlo (MC) simulation in rf discharges. Boswell et al.[79], 1988 Particle-in-cell (PIC) simulation of parallel plate RF discharge Surendra et al.[80], 1991 Birdsall et al.[67], 1991 PIC-MC technique for rf parallel plate glow discharge in helium Hitchon et al.[81], 1991 Direct numerical integration of the Boltzmann equation using BGK model for the collision term DiCarlo et al.[82], 1989 Direct integration of the Boltzmann equation using Flux corrected transport. Riley et al.[77], 1994 Direct integration Boltzmann equation using method of characteristics with a finitedifference treatment of the collision integral. Solution is accelerated using a time-averaged equation along with the full equation. Another complexity in solving the Boltzmann equation (2.1) comes from the evaluation of the collision te rm which is the source of th e non-linearity. Typically, the collision term is modeled using a particle technique such as Monte Carlo methods (DSMCdirect simulation Mont e Carlo). Other methods su ch as the particle-in-cell (PIC) technique are also common. Kushner [ 78] uses a particle simulation of a radio frequency parallel plate discharge using a M onte Carlo model. Boswell et al.[79] used PIC techniques in a self-consistent simulation of parallel plate discharges. Birdsall [67] and Surendra et al.[80] study the PIC-M CC (Monte Carlo collis ion) technique for modeling discharges. As an improvement over the early PIC techniques to include binary collsions (due to electrostatic effects), the Monte Carlo collision model was incorporated. This is in a way, a hybrid of the PI C and the MCC techniques which existed independently. PIC/Monte Carlo models of ra dio frequency discharges provide kinetic information self-consistently by integrating the equations of motion of many particles

PAGE 56

41 representing ion and electrons with a simulta neous solution of the Poisson equation with a Monte Carlo treatment of the electron-neut ral, ion-neutral collisions. The charged particles are given an initial velocity and the position coordinates and the particle paths are followed through time using the equation of motion until a quasi-steady state condition is achieved. These simulations are th en run for an extended period of time to filter out statistical fluctuations and obtain smooth profiles. The Monte-Carlo algorithm is employed to model the charge-neutral collis ion events following which the new position and velocity distribution is upda ted for the next time instant. The PIC [67] by itself is a method based on deterministic classical mechanics of particles with prescribed force fields. In PIC, the simulation regi on is divided into a number of cells and the resultin g grid is used in the solution of a force field from which the force on each particle inside the cell ca n be determined. The equation of motion of each particle is integrated to get the new position and velocity. On the other hand, MonteCarlo simulation is a statistical technique whic h is basically probabili stic in nature and the particle velocities are obtai ned by directly integrating th e equations of motion. The key steps of a integrated particle-in-cel l and Monte-Carlo technique (PIC-MC) of Nitschke et al.[68] are the following: a. Particle position and velocities are upda ted each time step using the equation of motion. b. Particles which have crossed the system boundaries are accounted for. c. Monte-Carlo algorithm to determine probabl e collisions and the resulting states d. New particle positions are linearly interpolated to the nearest grid point to obtain ion and electron spatial density profile s. Poisson’s equation (Eq. (2.2)) is

PAGE 57

42 calculated again to determine the electric fiel d to be used in the force term in step a. Other solution methods for kinetic models include direct solution of Boltzmann equation as in Hitchon et al.[81] where the equa tion is directly integr ated numerically as a partial differential equation. Here the collis ion term is approximated by a BGK type model (apart from Monte Carlo collisio n models). Because of the strong convective nature of the equation, higher order oscillation preventing schemes such as Flux Corrected Transport (FCT) (DiCarlo et al.[82]) or the “convectiv e scheme” of Hitchon et al.[82] are employed. The direct solution of the Boltzmann equation is slowed down by the very small time resolution dictated by electron and collision relaxation times. This problem is overcome by Riley et al.[77] using a hybrid formulation which employs a time-averaged description. Here the Boltzmann equation is solved for onl y a few cycles at a time. The information obtained from the Boltzmann equation solution is used to solve a set of time averaged equations describing the heavy species. The aver aged equations are solved for thousands of RF cycles and then the Boltzmann equation is solved again and this process is repeated. This is in a way another type of hybrid approach. Another technique along similar lines where the faster varying scales are solved a fe w times and the information is used to model the slowly varying processes is used by Hi chton et al.[81]. This scale up technique accelerates the solution process significantly. Other hybrid models will be presented next. The various kinetic studies are summarized in Table 7. 3.1.2 Hybrid Modeling Approach The computationally intensive kinetic simulations (as compared to fluid models) are more accurate in the sense that they can used to model low-pressure (typically < 100mTorr [68]) or highly non-equilibrium situations for which fluid models are highly suspect. The

PAGE 58

43 hybrid models are typically employed when there is a need to preserve the accuracy and capability of the kinetic schemes and at the same time reducing the computational burden. Table 8 Summary of hybrid modeling studies Author and year Hybrid model Belenguer et al.[84] ‘Bulk-beam’ model where the bulk and high speed or beam electrons are treated as fluid and particle species respectively. Surendra et al.[80] 1990 ‘Bulk-beam’ model Fiala et al.[85] 1994 Fluid model with Monte-Carlo simulation for the ionization term. Hybrid models have been studied which combine aspects of both fluid as well kinetic models and thereby achieve a balance between accuracy and efficiency. The charge species in low-pressure discharges are generally not at local thermodynamic equilibrium and one cannot assume in general that these species attain a thermodynamic temperature. This lack of thermodynamic equilibrium compli cates the modeling of discharges, since the kinetic models determine the species velocity distribution function as against an assumed distribution function in the case of fluid m odels. From another perspective, under low pressures the thermal flux is comparable to the drift flux from the electric field and as a result the isotropic behavior for macroscopic prop erties such as pressure is lost. The kinetic models in such non-equilibrium cases give acc urate results as compared to fluid models. One type of hybrid model treats the fast or high energy electrons called as ‘beam’ electrons as particles and the slow electrons or ‘bulk’ and ions as fluid and achieves a self-consistent simulation of the entire discha rge while addressing the non-equi librium nature of the high energy electrons. Such a model has been first introduced by Belenguer et al.[84] and used by Surendra et al.[80]. Another type of hybrid procedure is used by Fiala et al.[85] where

PAGE 59

44 the ionization source term is re-calculated by a Monte Carlo simulation after a few hundred time-steps of the fluid equations. 3.1.3 Fluid Modeling Approach Fluid models generally consist of a few moments of the Boltzmann equation for the various species with an assumed distribution function. However, as discussed previously, the validity of the fluid type de scription of the plasma is determined by the existence of a local thermodynamic equilibrium. This issue is significant while modeling low-pressure discharges. Nitschke et al.[68] compare th e particle-in-cell model with a fluid model simulation of a radio frequency discharge in helium. The results from the two models were compared over a wide range of app lied voltages and pressures. Considerable agreement between the studies was observe d for pressures greater than 100 mTorr [68,55], which is much smaller than the atmo spheric pressure regime of 760 Torr. As a result most of the atmospheric pressure disc harge simulations use fluid models which are potentially attractive for computational studies if the limitations to their validity are adequately addressed. In this dissertation, we are interested in the modeling of discharge processes at atmospheric pressures. Fluid models usually consist of the firs t three moments of the Boltzmann equation which includes the species continuity, mo mentum and energy conservation equations. The fluid species usually consist of the ions equations and neutral species. The number of species is typically dependent on the extent of plasma chemistry included in the model. A typical fluid description of the model as used by Colella et al.[86] is presented below. Species continuity equation: (2.3) j ij i i iR u n t n .

PAGE 60

45 Species momentum equation: jR j ij i i R j ij i i j i ij i j i j i i i i i i i i i i i iij ijR u m R u m u u n m m m m E n q kT n u u n m t u n m0 | 0 |. (2.4) This is a generalized form of the moment um equation incorporating various phenomena. The third term on the left hand side represen ts the partial pressure gradient of the ith species. The right side includes all source term sthe Lorentz force, momentum loss from collision and momentum losses from productio n and destruction, in that order. Species energy equation: j ij ij R j ij j j j j R j ij i i i i j j i j i i j j j j i i i ij i j i j i i i i i i i i i i i i i i i i i iS R u u m R u u m u u m m u u m u u m n m m m m E u n q u kT n u u u m n u u m n tij ij 0 | 0 | 2. 2 1 2 1 . 2 1 2 . 2 1 . 2 1(2.5) where 1 i i ikT Here the summation notation denotes the summation over the different kind of species. In assuming a fluid model, each plasma component (index ‘j’) is assumed to have a near Maxwellian dist ribution. The infinite set of moment equations is truncated to include just the moments corresponding to the conservation of mass, momentum and energy for each component. The indices ‘i’ and ‘j’ represents the particular species and its complementary ones, respectively. The summation notation indicates summation over all particle species which includes el ectrons, ions and neutral species. These transport equations are solved in conjunction with the Maxwell’s equations for the electric potential which is written as

PAGE 61

46 i i in q01 and E. (2.6) The number of equations that are solved in the system is determined by the number of species considered in the plasma chemistr y and can result in a hugely coupled system. However in reality, the number of equations act ually solved is very manageable with the various approximations to the model and the system of equations. One of the common assumptions [20] is approximating the plasma chemistry to include just a single species of positive ion and negative ion or electron. Also, in many cases, the neutral species are not considered because of their lack of c ontribution to the net charge. However there exists non-charged species such as metastables whose dynamics can differ from those of the neutral species. Such species are usua lly considered when the complex plasma chemistry is important for the modeling study. This is illustrated below in the context of a helium atmospheric pressure glow discha rge where a model with complex set of equations, as given in Tochikubo et al.[63]. Table 9 Plasma chemistry considered in a helium discharge 1 3 2 2 1 3 2 3 332 2 2 2 2 22 22 22 2 22 eHeHee eHeHeSe eHeHeSe HeHeHeHe HeeHee HeeHee HeSHeHeh HeSHeHeHe HeSeHee HeSHeSHeHee Source : Tochikubo et al. [63]

PAGE 62

47 The plasma chemistry shown in Table 9 needs to be modele d taking into account the reaction rates to account for the species c oncentrations and transport. Predominantly, modeling studies are one-dimensional (see Ta ble 4), followed by two-dimensional studies which are less frequent, while three dimensional studies are rare. It s hould be noted that distinct higher dimensional effects have been observed in studies by Mangolini et al.[62,56]. The fluid modeling studies have already been summarized in Table 4 and Table 5 Comparison of the results of different numerical models for the same set of modeling parameters has been compiled by Surendra [87] for a 1-D RF discharge simulation. The results indicate that fluid models, although less reliable than kinetic models, give reasonably accurate predictions of discharg e properties for sufficiently high pressure. Also, the kinetic models, say PIC-MC, compar e well with the hybrid models. The fluid models are generally reported to show great er variation than kinetic/hybrid schemes. Since we are interested in the discharge ope rating at atmospheric pr essure, we will focus on the study of fluid models for this dissertation. 3.2 Numerical Modeling of Glow Discharge Plasmas In this section, we will present the de tailed numerical modeling framework for the atmospheric pressure glow discharge pl asma. The plasma-fluid-Maxwell equations presented in the previous chapter represent the generalized set of equations. In reality, however, the set of equations solved are mu ch different and less complex from the ones presented, after the various problem speci fic modeling approximations applied to the system. In this section, we will first discu ss the various approximations that can be used from which we can arrive at the final set of equations to be solved. Following this, the numerical nature of these equations will be di scussed to highlight the difficulties that will

PAGE 63

48 be encountered during its solution and the techniques to handle them. Finally a solution algorithm will be developed to handle the sy stem of intricately coupled equations. 3.2.1 Plasma Discharge Governing Equations The fluid model for the plasma can be described by the following equations which describe the first three moments of the Bolt zmann equation for the various species. They are the continuity, momentum and energy equations. The plasma is considered as a multicomponent fluid comprised of two types of pr imary species, namely, ions and electrons. A generalized description of th e plasma model was presented in section 3.1.3 (Eqs. (2.3)(2.6)). The various modeling studies use a variety of approximations to simplify the system of equations. In the following sec tion, we will present the commonly used versions and highlight the specific va riants used in different studies. We will describe the composition of the fluid by means of electrons and ions. By ions we imply all charged species which have a mass comparable with the neutral particles. As seen in the previous chapter, the uncharged metastable species also play a significant role in the plasma chemistry. But we shall neglect them in favor of the present focus towards developing an efficient modeling strategy for these plasma fluid equations. Now, we will briefly develop the form of the fluid equations to be used. The continuity equations for the positive ion and electron species are given by ,,.e eeceden nvRR t (2.7) ,,.i iicidin nvRR t (2.8) The subscripts e and i stand for electrons and positive ions, respectively. The parameter Rc is the ionization or creation term and Rd, the recombination or destruction term. With no loss of generality, we can have a similar equation for negative ions. The

PAGE 64

49 terms on the right hand sides of eq. (2.7)-(2.8 ) represents the production from ionization and the loss from destruction, in that or der. These two terms are obtained from the collision terms in the actual Boltzmann equati on(2.1). The exact forms of these terms are obtained from the type closure model one uses to represent the particle phenomena at the fluid model. This can be done in a variety of ways as mentioned in [55] such as assuming a velocity distribution f unction (Maxwellian) or by assuming phenomenological expressions for the collision terms. The gene ral form of the source terms is given as cc dd R n R n (2.9) where c and d are the creation and destruction fr equency which are function of known parameters such number density ni,e velocity ui,e etc. The species velocities are given by ui and ue for the ions and electrons respectively. These velocities are obtained from the momentum equations given below. The genera lized momentum equation is given in(2.4). It can be written simply as .eeeeeeeeeeemenmvnmvvpnmqER t (2.10) .iiiiiiiiiiiminmvnmvvpnmqER t (2.11) where Rmi,me represents the collisional momentum loss of the ith species or electron and Pi,e is the partial pressure which can otherwise be written as ni,ekTi,e. k is the Boltzmann constant, Ti,e is the species temperature and mi is the species mass. This form of expression for the pressure tensor assumes isot ropic behavior and is strictly valid only for a Maxwellian distribution. The momentum loss Rmi,me is usually written as meeeemeRmnv (2.12)

PAGE 65

50 miiiimi R mnv (2.13) where mi,me is the frequency of collisions which result in momentum transfer. The second term on the right hand side represents the body force from the electric field with qi,e being the species charge. Combini ng (2.7)-(2.8) and (2.9) we have ,,.e eeceedeen nvnn t (2.14) ,,.i iiciidiin nvnn t (2.15) and substituting (2.12)-(2.13) into (2.10)-(2.11) we have .eeeeeeeeeeeeeeenmvnmvvnkTnqEmnv t (2.16) .iiiiiiiiiiiiiiinmvnmvvnkTnqEmnv t (2.17) Now we will simplify the momentum equation (2 .16)-(2.17) using (2.14)-(2.15), to yield. ,,.e eeeeeeeeecede eeeeeeemev nmnmvvnmv t nkTnqEnmv (2.18) ,,.i iiiiiiiiicidi iiiiiiimiv nmnmvvnmv t nkTnqEnmv (2.19) or ,,1 .ee eeee ecedee mememeeemeemenkT vvvqE vv tnmm (2.20) ,,1 .ii iiii icidii mimimiiimiiminkT vvvqE vv tnmm (2.21)

PAGE 66

51 where time derivative is a substantial derivative. Now we try to simplify the momentum equation using a scaling argument. The first term in (2.20)-(2.21), which is the local derivative scales, as the velocity times the ra tio of the driving volta ge frequency and the collision frequency. This ratio is usually very small even for driving voltages in the megahertz range as the collision frequenc y is of the order of a few gigahertz. ,, ,,1ieie mimemimevv O t (2.22) Hence the local derivative can be neglecte d. For the convection term, momentum loss collision frequency mi,me scales as the ratio of the thermal velocity and the mean free path. ,, ,, ,.eiie eiie memiTvv vOv Lv (2.23) As long as the thermal velocity is comp arable to the specie s drift velocity ( ,1ie Tv v ) and we are in the continuum regime ( 1 Kn L ), the convection term will be very small compared to the right hand side (force and collision loss terms). However, for the heavier ions the thermal velocity is comparable to the drift velocity at very low pressures and hence the inertial terms become significant. Hence, at low pressures the distribution function departs significantly from the Maxwel lian or equilibrium distribution function. This makes the diffusivity and pressure tensor s anisotropic. But at higher pressures, the above simplification is not too costly and th e ion momentum equation has the same form as the electron momentum equation. The third te rm on the left hand side of (2.20)-(2.21) is the ratio of the creation or destructi on frequency and the momentum loss collision

PAGE 67

52 frequency. It is usually expected that th e collision processes are more frequent when compared to one which causes ionization and the ratio is negligible. The final simplified form of the momentum equation can be written as ee e e eemeemenkT qE v nmm (2.24) ii i i iimiiminkT qE v nmm (2.25) Now let us define the species mobility i,e and the diffusivity Di,e as , ,,ie ie iemimeq m (2.26) , ,,ie ie iemimekT D m (2.27) The final form of the momentum equation becomes ee e e enD v E n (2.28) ii ii inD Ev n (2.29) This form of the momentum equation is cal led as the drift-diffusion approximation. Further variants of this form have been used in many studies such as by Boeuf et al.[31] etc., and is shown below. e e e e en D v E n (2.30) ii ii iDn Ev n (2.31)

PAGE 68

53 Now we will briefly derive the simplified form of the energy equation. The general energy equation can be written as ..Ree eeeeeeeeeenn nvvPQnqEv t (2.32) where is the particle energy, Pe the electron pressure, Qe the heat conduction and Ren the energy loss from collisions etc. As before assuming a Maxwellian distribution for the velocity, we have ,2 3eeeTePnkTn (2.33) and ,2 3Te eeeeQKTK k (2.34) Here T is the thermal energy component. The thermal conductivity Ke can be written in terms of the diffusion coefficient De, from kinetic theory results as 3 2ee enDk K (2.35) The fact that all the parameters such as conduc tivity, diffusivity are written as scalars is because we assume isotropy. Substituting(2.33),(2.34) and (2.35) into (2.32) we get the complete energy equation as ,2 ..R 3ee eeeeeeeeTeeeeenn nvnvnDnqEv t (2.36) The form of the energy loss rate Ren as well as the parameters Rmi,me, Rc, Rd, the mobility and diffusivity etc, needs to be identified to close the system. These parameters are usually related to some measurable macr oscopic properties. Boeuf et al.[31] and Gogolides et al.[55] among othe rs assume that the above ener gy loss rates etc., depend on

PAGE 69

54 the mean electron energy as they do under equilibrium conditions when the electron energy gain is locally balanced by the losses. If we assume the following form for Ren [31, 55] /enLeLe R KENNnn (2.37) then the energy equation under equilibrium conditions yields ee LqvE K N (2.38) At equilibrium, all the mean electron properties such as e, De and e etc. can be assumed to be a function of E/N. This dependence of the mean electron properties on the local electric field as against the velocity distri bution function (in the Bo ltzmann equation) is called as the local field approximation. This approximation rules out any non-local effects for these properties. The values of th e electron properties are the same as that obtained in a uniform electric field condition, where the collision term balances the force term. Detailed discussions are presented in [66,64,72,55] Hence, from (2.38), Ren is also a function of E/N. A similar form is used for other parameters such as Rmi,me, Rc, Rd. The equilibrium dependencies of the various pa rameters on the ratio E/N are tabulated by solving a Boltzmann equation or a Monte Carl o code under the appropriate conditions. It should be noted that in some hybrid approach es (e.g. Kushner and co-workers [88-89]), the fluid energy equation is dropped and the kinetic energy equation is solved using a 2-D Monte Carlo method. The ion energy seems to have little eff ect on the ionization and attachment rates and hence is not solved unless required. To summarize, the electron inertial terms are insignificant considering their negligible mass. Hence, the momentum balance reduces to the drift diffusion form. For the ions, the iner tial term is significant at low pressures (<

PAGE 70

55 100 mTorr [55]) and is included in some studies. At high pressures, however, a driftdiffusion approximation is considered reason able [31,55]. A summary of conservation equations for the mass, momentum and ener gy used to model high frequency glow discharges is given below. A similar set of equations is presented by Gogolides [55]. Continuity equation: Electrons: ,,.//e eeceedeen nvKENNnKENNn t (2.39) Ions: ,,.//i iiciidiin nvKENNnKENNn t (2.40) Momentum equation: Electrons: eee eeennv EnD (2.41) Ions: iii iiinnv EnD (for high pressure discharges) (2.42) ,,.i iiiiiiiiicidi iiiiiiimiv nmnmvvnmvNKK t nkTnqEnmvNK (2.43) (for low pressure discharges) Energy equation: Electrons: ,2 3 ./ee eeeeeeeeTe eeeLen nvnvnD t nqEvKENNn (2.44) Ions: No energy equation is considered as i on energies do not affect ionization processes. To summarize, the various assumptions made in the development of the fluid equations are: 1. the pressure tensor is isotropic 2. the drift energy is comparable to the thermal energy

PAGE 71

56 3. heat flux is proportional to the electron temperature gradient 4. the mean electron-neutral collision ra tes depend only on the electron mean energy. Most of the numerical modeling studies of hi gh frequency glow discharges use whole or part of the above system of equations along with the Poisson equation for the electric field which is written as. 0.iee Enn (2.45) The values of the parameters K, mobility, di ffusivity etc. are either obtained by solving the spatially homogeneous Boltzmann equation, or assigned curve fitted expressions from experimental data or a combination of bot h. A standard model for the ionization coefficient is the Townsend expression pres ented in [61] which is shown below. 0.416(34)E p e iEp Ke N (2.46) Here p is the pressure. Some of the numeri cal issues associated with the system of equations discussed above are as follows. Th e set of moment equations described above are intricately coupled and s trongly non-linear with a wide ra nge of spatial and temporal scales. The solution profiles for the charge particle densities c ontain steep unsteady gradients especially in the sheaths. The resolution of these sharp gradients accurately requires well refined grids leading to time -intensive computations. Sharply varying production terms such as ionization processes re sult in stiff source terms which also limit the time step. For these reasons, 3-D and to a lesser extent 2-D modeling studies are less frequently attempted. The continuity equation for the ions and electrons are essentially inhomogeneous Euler equations (as diffusion is mostly small compared to drift) of gas

PAGE 72

57 dynamics augmented with the electrostatic force, collisional and source/sink terms. Solution of this system poses challenges due to its multi-time-scale nature which we will illustrate now. The critical time scales in the problem are listed below. 1. The discharge operating frequency 2. The ion drift 3. The non-linear source termionization, collision etc. 4. The electron drift 5. The dielectric relaxation time The discharge operating frequency time scale is given by 1 (2.47) where is the driving frequency. In order to illustrate the time scales we will use the 1D modeling result from Lee et al.[34] which is shown in Figure 3. Th e driving voltage is 1.5KV with a frequency of 10 KHz. It can be seen from the figure that the appropriate length scale for the problem would be the size of some physical length such as the positive column or the cathode fall (Figure 3) which is 0.05 cm in a 0.3 cm gap. The drift velocity from the drift-diffusion approximation isE The values of the mobility for the helium ions and electrons can be obt ained from literature and the electric field value is the peak value in the positive column (~ 20 KV/cm). Then the drift time scales for the ions and electrons can be calculated as the ratio of the width of the positive column and the drift velocities. dr drion ionl E (2.48)

PAGE 73

58 dr drelec elecl E (2.49) To evaluate the source term time scale, we will use the Townsend formula for calculating the ionization coefficient. From (2.9) the ionization term Rc can be written in terms of the ionization frequency c. From (2.37) and (2.46) we can see that c is of the form 0.4161 (34)E p ce cEpe (2.50) which gives the ionization time scale c as the reciprocal of the frequency. We use the above form of the ionization frequency in ca lculating the timescale in Table 10. Now we shall describe the origin of the dielectric relaxation time scale and evaluate it. This timescale represents the relaxation of the elec tric field in the medium to adjust to the change in the charge density and is usually the most stringent of all the time-step limitations. Taking the time derivative of the Poisson equation(2.45), we get 0.ienn e E ttt (2.51) Now substitute for right hand side from the continuity equations (2.39)-(2.40), which gives 0..eeiiEe nvnv t (2.52) Getting rid of the divergence operator, we get 0 eeiiEe nvnvC t (2.53)

PAGE 74

59 Substituting for the velocities using the d rift diffusion approximation (2.28) yields 0 eeiiiieeEe nEnEDnDnC t (2.54) Eq. (2.54) can be viewed as the evolution of the electric field due to the space charge and leads to a time scale as 0 diel eeiienn (2.55) This restriction is highly severe and is usually overcome by a semi-implicit solution of (2.54). In other words the above restricti on can be viewed as the CFL condition for the current continuity equation. Using the values of the various parameters from the 1-D model of Lee et al.[34], the magnitudes of the various time scales are calculated and tabulated in Table 10. The dielectric relaxation time becomes real small as the number density increases which is the case for high frequency (MHz) discharges and hence needs to be treated accordingly. As can be seen from Table 10, time scales involved in the problem vary over a wide range and this needs to be taken into account while developing a numerical algorithm. If we take into account adequate grid resolution in the sheath regions, this time scale will be further reduced. 3.2.2 Numerical Algorithm Modeling of plasma fluid equations has b een extensively attempted in the last decade or so predominantly for 1-D applications and occasionally 2-D applications. Emphasis was placed on plasma chemistry [ 89] and studying the limitations of the model by comparison with experiments or kinetic th eory based models [68,31,67,80]. Efforts to develop efficient solution algorithms for such systems such as Boeuf et al.[31] Hammond et al.[90] Golubovskii et al.[56] Colella et al.[86] etc. will be discussed in the following

PAGE 75

60 section. These are representative studies an d a more comprehensive review has been presented in the previous chapter 2. Except for Golubovskii et al. where an atmospheric pressure glow discharge between two parallel plates is studied, all other efforts are devoted to low pressure high frequency di scharges. The set of governing equations employed by the various efforts differ from each and are summarized in Table 11. While the equations have the same form as given in (2.39)-(2.44) the methods of evaluating parameters such as the mobility, ionization et c. may differ form each other. It may be recalled that the two forms of the momentum equations ((2.42) and (2.43)) which include a fully approximated drift-diffusion form and the full transport equation. Both these forms are used in the modeling studies, but the simpler drift-diffusion form is preferred and valid at high pressures for ions. It is generally always used for electrons because of its light mass causing negligible inertial effects. Table 10 Representative time scales of va rious processes in a 1-D helium discharge model operating at 10KHz, 1.5KV rms, calculated based on values of parameters. Time scale Value (s) Operating frequency timescale, w (2.47) 1e-4 Ion drift drion (2.48) 2.5e-7 Electron drift drelectron (2.49) 2.5e-9 Dielectric relaxation diel (2.55) 2e-9 Ionization timescale c (2.50) 1.4e-10 Note: Calculate based on parameters taken from Lee et al.[34] It can also be seen from Table 11 that the electron energy equation is not solved for high pressure simulations (Lee et al.[34] Massines et al.[6]) when the local field approximation is valid. Not solving for the en ergy equation, however, prevents us from understanding the energetics such as collisional heating or Joule heating. While, the choice of governing equations and their approximations vary from one study to the other,

PAGE 76

61 the numerical characteristics remain fundame ntally the same. The governing equations are essentially the inhomogeneous Euler eq uations of gas dynamics augmented with electrostatic forces, collisional and source/sink terms. The system of equations represents multiple time scales (Table 10) and widely se parated spatial scales. The spatial scales can be adequately resolved by using local grid refinement techniques [116] in the regions of sharp gradients coupled with higher order discretization methods which have been successfully applied for modeling neutral fluid flows. Table 11 Summary of governing equa tions used in different studies Governing equations used Continuity Momentum Energy Author, year Topic Ion electron ion electron ion Electron Boeuf et al.[31] 1995 2-D rf discharge in argon at 100mTorr (low pressure) (2.40)(2.39)(2.42)(2.41) (2.44) Colella et al.[86] 1999 Develops a conservative finite difference method for plasma fluid equations for high density low pressure plasmas (2.40)(2.39)(2.43)(2.41) (2.44) Hammond et al.[90] 2002 Numerical modeling of low-pressure, high density rf discharges in Helium (2.40)(2.39)(2.43)(2.41) (2.44) Lee et al.[34] 2004 Golubovskii et al.[56] 2004 1-D parallel plate discharge in helium at atmospheric pressure (2.40)(2.39)(2.42)(2.41) Roy et al. [15,115] 2005 2-D asymmetric discharge actuator modeling in air/He (300 torr) (2.40)(2.39)(2.42)(2.41)

PAGE 77

62 It should be kept in mind that the choice of convection scheme is also determined by its performance in the range of Peclet numbers encountered and also in the presence of a source term [91]. In the following discussion we will restrict ourselves to the set of governing equations as presented in (2.39)-( 2.44).The challenge here is to efficiently handle multiple time scales, and the treatment of the reaction source terms. Such complications are more common in the modeling of reacting flows where detailed modeling of chemical kinetic mechanisms can result in a stiff system of governing equations. The focus of the present dissertati on includes the development of an efficient solution algorithm to handle the stiff multi-scale problem as seen in the context of the plasma dynamics. The discharge equations as described in (2.39)-(2.44) involve wavelike ionization fronts which travel through the doma in before charging the dielectric or any other boundary surface. The solution procedure needs to be designed to handle the multiple scales efficiently. However, one distinct difference between chemically reacting flows which pose similar challenges is the electr ic field effect. It has been shown in [86] that the plasma sound speed is modified by th e electric field drift wi th the effective ion sound speed given by eii i ikTkT c m (2.56) Here Te refers to the electron temperature which is usually very high compared to the ion temperature (which is closer to the neutral particle temperat ure). Thus the modified sound speed is much higher than the neutral sound spee d (as ion temperature is ~ neutral). Note that the first term in the square root with Te comes from the approximation of the electric field in the electrostatic pressure. In the ab sence of a driving elect ric field, this term vanishes giving the neutral fluid sound speed (the second component in the square root).

PAGE 78

63 When Te >> Ti ~T, the electric field dominates the thermodynamic pressure. Similarly, the characteristic speed for electrons can be obtained as eee e ekTkT c m .(2.57) The above expression for the ions shows that the electric drift component is more significant than the thermal component. The ra tio of the two species sound speeds goes as the inverse of the square root of their resp ective mass. The electrons being lighter than ions are faster. In general, from the range of electric fields considered in this problem, the ion and electron sound speeds are much higher th an that of the neutral flow and hence the bulk of the acoustic effects are caused by the driving electric field and not the low speed neutral flow. This is different from normal chemically reacting flows where the species characteristic speeds are comparable to that of the neutral fluid. This distinction needs to be kept in mind while dealing with the a bove strongly convective equations where the reaction and convection time scales are comparable. Integrating such a system comprising ma ny time scales gets very demanding in terms of computational resources. Both implicit and explicit techniques have been used to march in time with the choice usually based on the extent of stiffness of the system, expected time accuracy etc. Implicit integra tion methods are attractiv e in the sense that they are unconditionally stable, but are not time-accurate and in some non-linear situations are even computationally burdensom e. This is necessitated by the need to find the root to a highly non-linear equation which might entail very small time-steps to achieve convergence. Also, there is a burde n of huge computational requirement. This method works reasonably well when the overall timescales are not very far apart. Several studies aimed at improving the efficiency of stiff integration schemes have been studied

PAGE 79

64 in the context of modeling reacting flow a nd atmospheric chemistry which involve taking advantage of the problem specific approximations to the Jacobian. In other examples, the stiff variation of the reacti on terms which is typically e xponential makes it advantageous to study the variation of the logarithmic quant ities as against orig inal variables [90]. Another alternative to the implicit procedur es used above is the use of matrix free Newton-Krylov methods, which have been reviewed in [92]. Both sequential and fully implicit approaches have been used for discharge modeling. Roy et al.[15-16] used a globally im plicit finite element procedure where the system of species continuity and momentum equations are assembled as part of a global matrix to solve for the solution vector. Sequential approaches have also been employed such as in the study of Hammond et al.[90]. Hammond et al.[90] uses a hybrid formulation to solve the plasma fluid equati ons, where the ions species being slower are treated by an explicit 4th order Runge-Kutta method followed by an implicit treatment of the fast electron equations. The implicit Eu ler (first order) and implicit Runge-Kutta (second order) methods were used with Ne wton-Raphson iterations to overcome the nonlinearity. Further, the resulting Jacobian is simplified by neglecting the weak off diagonal contributions reducing to a block tri-diagona l system. However, taking into account the computational overhead from the iterative procedure for the fully implicit treatment and the cost of inverting the matrix more than compensated for the gain over explicit procedures. Implicit procedures though are us ed in the context of obtaining steady state solution with stability and not time accuracy as the consideration. Acceleration to steady state is usually achieved using preconditioni ng techniques (for e.g. [93]) which decrease the condition number of the system with widely separated scales. These techniques have

PAGE 80

65 been extended to time accurate computations by Merkle [94]. Schwer et al.[95] uses a consistent splitting of operators with implic it transport computations of reacting flows using local preconditioning (Merkle et al.[94]). Recent modeling strategies use semi-implicit hybrid procedures which seem to be more attractive especially in the context of operator split methods. Here the stiff portions of the system are handled using efficient solvers and integrated with the non-stiff solution procedure. A common practice to handle the stiff reaction terms is to integrate the reaction ODE using a stiff ODE solver (CVODE, DVODE, LSODE) over a time step determined by a non-stiff process. These are essentially based on the backward differentiation formulae of Gear. This integration time-scale in most problems is the convection time-scale of species which is considerably larger. Such a procedure makes for a sequential approach and thereby enable s handling smaller time scales by allowing integration sub-cycles for the relevant portion. However, stability and coupling demands appropriate corrector steps and handling of splitting and conservation errors. Source term treatment in plasma simulations has been discussed by Hagelaar et al.[96] where a linearized implicit correction is used to alleviate the stiffness by means of a semi-implicit solution. This is similar to another common technique which is the pointimplicit treatment where the source terms ar e treated implicitly with a single Newton iteration. Here the source term is linearized and written in terms of a Jcobian which needs to be evaluated every time step. Also, if there are a large number of species with complicated chemistry, inverting the sparse Jacobian becomes expensive and it is common practice to resort to simplifications to accelerate convergence. Updating of the Jacobian is also done only when absolutely necessary and is not usually done for every

PAGE 81

66 step unless it seriously affects convergence. It is shown in [ 97] that this point implicit method is a special case of a generalized time scale preconditioner where the stiff processes are scaled in time to the speed of slower processes. Preconditioning techniques used to handle multiple scales to solve steady state problems and are reviewed by Turkel [93]. Implicit treatment of s tiff terms provides solutions which are accurate at the slow scales and stable at the fast scales. Overa ll, a fully implicit technique is rarely used, especially with the strong coupling between the various species equations and the resulting non-linearity. However, with the us e of efficient solution techniques such as implicit multigrid methods and efficient solv ers like LU-SGS of Yo on et al.[98] high speed reacting flows with strong convection ha ve been successfully handled implicitly. In spite of using inversion techniques as above the species chemistry Jacobian makes it a block inversion which can get expensive. Hence, it is common practice to simplify the species chemistry either adaptively or uniformly or to treat the reaction terms in a partially implicit fashion so that block inversion can be reduced to a scalar inversion. Explicit methods have been developed to deal with stiff problems, such as stabilized multi-stage Rung-kutta schemes, Runge-Kutta-Chebychev (RKC) [99] schemes which possess extended real stability intervals. However, with the stability range dependent on the number of stages employe d, the computation increases for handling very stiff problems and is almost comparable to implicit schemes. A review of explicit schemes for stiff problem can be found in Verw er [99]. Explicit time-scale splitting using computational singular perturbation (CSP) is presented by Valorani et al.[100]. Here, the fast and slow time scales are split such that the algorithm marches with the slow time

PAGE 82

67 scale with the faster time scales being take n into account at the e nd of the integration time-step as a correction. Operator splitting methods are used in reacting flow computations and are discussed in [101]. The treatment of source te rm using operator splitting has been popular in the computation of reacting flows and di scussed in Najm et al.[102] and Oran and Boris [101,102]. Here, the stiff reaction operat or is treated as a stiff ODE using the DVODE integrator [103] and is embedded as part of a projection method. Colella et al.[86] employs a time-split integration algorithm of the species equations in a sequential form to decouple the three possible stiff pr ocesses namely ionization, convection and dielectric relaxation. The different speci es equations are handled in a sequential formulation in a predictor-corrector approach wherein the evolution of the various species is appropriately coupled to the electric field. In the present multiple study, we have employed an operator-split sequential solution algorithm. Some examples of studies which employed such ideas are Najm et al. [102] in combustion problems, Verwer et al. [99] in atmospheric chemistry modeling and Tyson et al. [104] in chemotaxis models. He re, the solution procedure can be adapted to handle the individual processes efficiently and realize overall ga in in computation. However, such an approach requires careful attention to stability considerations and performance is highly dependent on the physics of the problem. In certain situations the presence of competing stiff processes can l ead to a system with much larger dynamical time-scale than dictated by the individual processes. Such dynamic equilibrium usually cannot be predetermined in general for a non-linear advection-diffusion-reaction system. The stiffness of the reaction part is typi cally overcome by using stiff integration

PAGE 83

68 procedures in ODE integration packages [ 105] such as the ones based on the backwarddifference formulae (BDF). In using time-sp lit algorithms for processes operating in a range of time-scales, the choice of time-step si ze is typically determined by the smallest time scale, but need not necessarily be chosen as such. To speed-up the solution procedure, an intermediate time scale is chosen to advance the overall system in time, while the faster processes are advanced by sub-cycling within the time-step. In the present study, the time-step dictated by the sl ower ion species convection is targeted to march the full discharge system while sub-cycling is used for the faster processes. Also, a predictor-corrector approach is employed to ensure sufficient coupling between the electric field and the species densities. A s trong coupling is essential for achieving stable time-accurate simulations while using a sufficiently large global time-step. This method is integrated with a multi-block finite -volume algorithm capable of handling 3-D curvilinear-grids [19]. The implemented solution algorithm conserves space charge and avoids dielectric relaxation time step restrictio ns. In the following section, an outline of the algorithm is presented. The split solution algorithm consists of the following steps: 1. Predictor step 2. Solving the Poisson equation 3. Corrector step 1. Predictor step: At the beginning of the n+1th time level and for the kth species, we haven kn,n kv,nE and n k as the known quantities. Here the species continuity equations (2.39)-(2.40), along with the drift-diffusi on momentum equations (2.41)-(2.42) are integrated using lagged values for the various coefficients (as they are a function of the

PAGE 84

69 electric field E). The source term is integrated using a higher-order (4th-order) BDF using the CVODE solver [103]. The convection and di ffusion operators can be treated either implicitly or explicitly. In this case, we will employ a second-order upwind for the convection term and second-order central difference for the diffusion term. The continuity equation can be written as ,.k kkckkn nvn t (2.58) As can be seen from Table 10, the presence of processes with disparate time scales can possibly be better handled using operator splitting. Th ree types of splitting are popular, namely, (i) First-order splitting (ii) Strang splitting (iii) Source splitting. (i) First-order splitting: The first-order splitting [101] can be written symbolically as kM n kknTtStn (2.59) where S is the reaction operato r integrated using the ODE solver CVODE [103] and T is the transport operator. k M is the number of substeps used for the transport term integration to march to the global timestept Therefore, we have kt t M (2.60) (ii) Strang splitting: The Strang’s time splitting method [106] solves for the conservation law without the source term and the ordi nary differential equation (ODE) with the unsteady term and the source term in an alternating fashion. This method is formally second-order in time which is achieved by splitting the operators symmetrically. In this

PAGE 85

70 case, the transport term integration is usua lly split into two halves to achieve the symmetry since the ODE solver used for the reaction part is computationally burdensome. /2/2kkMM n kknTtStTtn (2.61) It is worth noting that, in both the first-order and Strang splitting procedures the initial guess for the reaction part is not directly from the previous time-step, but after a half or full time-step of the transport term integrati on. This results in the introduction of stiff transients in the solution which are nothing but an artifact of the splitting errors arising from the lagged treatment of the coupled proc esses. These can have significant impact in the presence of strong non-linearities. Even though the Strang splitting is formally second-order, it is rarely achieved [107] in certain stiff problems, where it is known to deteriorate to lower order accuracy. To overcom e the solution discontinuities which gave rise to stiff transients in th e above two splitting methods, the source splitting treats the transport as a piecewise constant source. (iii) Source splitting: In this study, we limit the discu ssion to first-order source splitting [99]. For the source spli tting, we can write n kknStn (2.62) and n kk knn n t (2.63). Therefore, we have n kkknTtntn (2.64)

PAGE 86

71 where S is the reaction operato r integrated using the ODE solver CVODE [105] and T is the transport operator. k M is the number of substeps used for the transport term integration to march to the global timestep t Thus, we have kt t M (2.65) The ODE solver employs the following: A 5th-order BDF for time integration. Newton iteration for non-linearity. A direct method with a banded treatment of the Jacobian. Normal mode with subcycling within the time-step. Relative and absolute tolerances of 1e-12 and 1e-14 respectively. The above strict tolerances were chosen so that the ODE integration is almost exact. After the predictor-step, we have obtained the guessed values fo r the species number densities based on the lagged value for th e electric field. The strong super-linear dependence of the convection on the electric field dictates the overall stability and requires strong coupling between the variables. 2. Solving the Poisson equation: Here we solve a semi-implicit version of the Poisson equation to overcome the space charge stability constraint leading to the dielectric relaxation time-step constraint. This time step re striction is usually one of the most severe and is equivalent to the CFL type stability criterion for the current continuity equation. The Poisson equation is 0.iee Enn (2.66) The origin of the restriction has been illust rated earlier in this section in (2.51)-(2.55). This restriction is primarily caused by the nonimplicit treatment of the electric field in

PAGE 87

72 the species momentum equation. In order to overcome this, a linearized implicit treatment for the species number densities is used. 111 00.n nnnnn ie ieienn ee Ennnnt tt (2.67) which can be written using the sp ecies continuity equations as 1 0...n nnn ieiieee Enntnvnv (2.68) Note that the only term in the right hand side which contains the velocity and hence the electric field (through the drift-diffusion form of the momentum equation) is the one inside the brackets witht By treating this term implicitly we can overcome the dielectric relaxation time step restriction. The implicit version of the equation can be written as 1 1 0 11 0 11 0... .. ..n nnn ieiiee nnnn ieiiee nnnnnn ieiiiieeeee Enntnvnv e nntnvnv e nntnEDnnEDn (2.69) Here the ion density at the time level n+1 is approximated by the predicted density at the end of step 1. The species velocities are directly substituted using the drift-diffusion approximation so that the electric field can be treated implicitly. In the event, the full species momentum equations are solved in step 1, the predicted velocity thus obtained (with the bar) will be used in the Poisson equation. In that case, the non-linearity will be difficult to overcome and will need a Newton-Raphson type treatment. The equation (2.69) leads to a symmetric system for the el ectric field which needs to be inverted. This

PAGE 88

73 can be done using a point Gauss Seidel method with successive over relaxation (GSSOR) or a line SOR. The cost of this st ep is possibly one of the most demanding and hence needs to be performed as less frequen tly as possible. The choice of the global timestep t will be determined taking this into account. 3. Corrector step: At the end of the previous step, the predicted species densities knand the electric field at the new time level En+1 are available. Now we will obtain the corrected densities at the new time level. It is worth noting here that the corrector step is needed to ensure adequate coupling between the electric field and the species number densities. In other words, we require appropr iate coupling between the species continuity equation and the momentum equation (or the drift-diffusion equation). Hence, the corrector step is the same as the predic tor step, but performed with the updated coefficients using En+1. It is also worth noting that the if the diffusion term and the source terms are not strongly impacted by the electric field, then the corresponding operators can be removed from the corrector step as correct ion is required only for the convection term in such cases. Stability. The stability criterion for the above predictor-corrector approach with stiff substeps is not straightforward. Of the time s cales listed in Table 10, the ionization source term and the dielectric relaxation time step restrictions are not binding in the above algorithm. The implicit treatment for the r eaction source term makes it unconditionally stable and is limited only by the non-linearity. The predictor-corrector formulation with the semi-implicit solution of the Poisson equation ensures space charge stability overcoming the dielectric relaxation time re striction due to the non-linearity. With the split integration procedures for the convection and other processes, the stability of the

PAGE 89

74 system will be determined by the CFL conditi ons for the individual steps. The global integration time-stept is determined by the slower moving ion CFL condition. Here is the CFL number for the different species. max ,iCFLit v x (2.70) In the case of a number of different species considered, the choice of the global time-step is usually the convection time scales of th e slowest species, unless there is a huge variation. This is important if we are consid ering the neutral species in reaction chemistry whose dynamics are determined by the slower moving neutral fluid. Once this is fixed, the time-step for the faster electron species and the corresponding number of sub-steps M, can be obtained, so as to satisfy the electron CFL limit. However, studies [108] indicate that the CFL stability limits of the split algorithm deviates from that of unsplit algorithm depending on the number of sub-steps M, employed. Knio et al.[108] shows that the critical CFL number due to el ectron convection sub-stepping decreases monotonically as the number of sub-steps in creases in two-dimensional modeling studies of reacting flows. Also, one-dimensional stud ies indicate a limiting value of the critical CFL being achieved as the number of sub-steps is increased. In spite of the stability criterion becoming stringent the overall computational savings is generally substantial. 3.2.3 1-D Parallel Plate Discharge Modeling in Helium We study the following test problem to validate our computational capability. A 1-D parallel plate (5mm gap be tween electrodes) discharge at atmospheric pressure using just two-species chemistry, driven by a vo ltage of 1.5 KV at 10 KHz was modeled as shown in Figure 8. This model problem has been studied by Ben Gadri et al.[20] and Massines et al.[6]. The governing equations ar e the species continuity equations ((2.39)-

PAGE 90

75 (2.40)) and the momentum equations with the drift-diffusion approximation ((2.41)(2.42)). No energy equation is solved for this model. The gas used is helium with a small impurity of Argon (0.5 percent) to simula te the Penning effect arising from the impurities. The parallel plates used are cove red by a layer of dielect ric with permittivity 9.0 and 0.6 mm thick. A detailed description of the arrangement is pres ented in [20]. The various transport parameters, such as the species mobility, diffusivity and ionization coefficient are tabulated as a function of th e electric field using the ‘BOLSIG’ database. A simplified kinetic model with various collis ion cross sections is used to develop the database for the various types of species and these are locally interpolated in the code. The boundary conditions are such that the fluxes near the boundary charge the dielectric which is very significant in terms of maintaining a controlled discharge. However there are inherent difficulties in the fluid model of the plasma in representing the near wall sheath regions. Hence differe nt sets of boundary conditions have been employed in different studies [34,58,62]. For th e 1-D species continuity equation we need two boundary conditions in space for the number density. A few of the different choices of boundary conditions employed are shown below in Table 12. A detailed discussion of the issues in the choice of boundary conditions for fluid models of plasma is presented by Hagelaar et al. [109]. In th e present discussion we employ set of boundary conditions 1 from Table 12. A uniform time-step of 10-9 seconds is used for the time-integration and the domain is uniformly discretized into 201-points. A fine grid with double the number of points is also used and th e results are compared with th at of Massines et al.[6]. The discharge current over a period of one full cycle after a quasi-steady state has been attained during the discharge is self-sustain ed and is periodic as shown in Figure 9.

PAGE 91

76 Table 12 Different types of boundary conditions for the ions and electrons used in modeling studies. Ion BC’s Electron BC’s Serial no Cathode Anode Cathode Anode 1 0n x N=0 Sec emission Flux 0n x 2 0 n x 0 n x Sec emission Flux 0 n x 3 4thnv En 4thnv 4thnv+sec emission 4thnv En Note: Vth is the thermal flux at the boundary. Figure 8 Parallel plate arrangement for modeling DBD discharge Figure 9 Evolution of discharge current over one period of the discharge cycle u pp er electrode Lower electrode Vext=1.5KV @ 10KHz 0 V Dielectric with dielectric constant 9.0( / 0) and thickness 0.6 mm(Ll,Lu) Gap=3.8 mm(L)

PAGE 92

77 The curves compare the coarse and fine g rid results obtained in the present model with that of Massines et al.[6]. As can be seen from Figure 9, the current peak representing the instant after the breakdown slightly differs from that obtained in the reference study. Also, there is a slight difference in the magnitude of the peak current. The possible source of difference could be traced to the slight differences in the choice of material data for the various species employed. However, the overall physics is adequately represented although qua ntitative differences exist. Figure 10 Evolution of the gap voltage over one full cycle The gap voltage which is the net voltage between the two dielectric layers (or the air gap to be specific) after a ccounting for the charge build up on the dielectric surfaces is shown in Figure 10. The deviation from the reference results is again seen, but the coarse and fine grid solutions are not much differe nt, indicating possible grid convergence. The peak gap voltage is at the instant just before the breakdown and it drops sharply after the discharge peak, with large currents as seen in the previous figure reducing the voltage. The distribution of the two species especial ly, the ions and the electrons, during the

PAGE 93

78 discharge are of interest, esp ecially for the present study where the force acting on the fluid needs to be accurately calculated. Th e species number density distribution over the domain at the instant of the peak current density is presented in Figure 11. The comparison shows that the present calculations have a ten percent error in the peak value of the ion and electron densities. The fine gr id results show identical peak densities and hence not shown here. The very small time step required in the above study makes operator split treatment of the various stiff processes necessary. This split approach will be even more valuable when the plasma dynamics need to be solved coupled with the fluid dynamics. As the faster fluid dynamics e volves, they will interact significantly with the evolution of the discharge. Figure 11 Comparison of peak current species density distributions of the present effort with Massines et al.[6]. 3.3 Integration of Plasma Dynamics with Fluid Dynamics The above discussion pertains to the modeling of the high frequency glow discharge independently which needs to be in corporated in the Navier-Stokes calculations for modeling the coupled system. The first issu e is the effect of the neutral flow on the plasma dynamics which needs to be captured. The actual species momentum equation in

PAGE 94

79 a neutral flow would need to include the pr essure gradient contribution from the fluid dynamics. We can argue the total pressure as a sum of the partial pressures of the plasma species and the neutral fluid particles. The momentum equation (2.41)-(2.43) should also account for the neutral fluid pressure gradie nt. Hence the modified momentum equations in the drift-diffusion form are Electrons: eeeeeefluidnvnEDnp (2.71) Ions: iiiiiifluidnvnEDnp (2.72) where fluid p is the pressure effect of the neutral fluid. Since we are using the two-fluid model for the plasma-fluid mixture and with the plasma sound speed being much greater than the neutral flow sound speed, the effect of the neutral fluid transport on the species dynamics will be much smaller than the electric drift for the charged species which makes the two pressure gradients in (2.71)-(2.72) negligible. However, the neutral reactant species such as the metastables will still convect on the fluid dynamic scales and evolve at a different rate. Therefore the modeling strategy can use a coupled solution method where the above time-split integration of the plasma model is embedded into a fractional step Navier-Stokes solver such as PISO [110]. In such an arrangement, the corrector step just corrects for the fluid dyn amics part of the plas ma species equations. This is an approach which is used in the modeling of reacting flows and is computationally demanding. An approximate approach would be to neglect the effect of the neutral fluid flow effects in the plas ma discharge equations. This would entail a discharge solution independent of the fluid dynamics which can be solved just once and then be incorporated into the Navier-Stokes equations as shown in Chapter 4. Since the most time consuming part of the solution algorithm is the modeling of the plasma

PAGE 95

80 dynamics, the latter approach makes for a computationally efficient solution procedure. In chapter 4 a framework for representing the plasma effects as a time-averaged body force in the fluid flow is presented. An em pirical expression for the body force based on assumed discharge characteristics is used to model the plasma-fluid effects.

PAGE 96

81 CHAPTER 4 MODELING FLUID FLOW WITH PLASMA EFFECTS In this chapter a framework is presented to model the effects of the discharge on the fluid dynamics. The model for the morphology and time evolution of the glow discharge was presented in the previous chapter. In order to couple the plasma dynamics to the neutral flow an appropriate strategy ne eds to be used. This usually depends on the problem at hand and the operating conditions. Fo r example, the flow physics needs to be fully coupled and integrated with the discharg e plasma simulation if the time scales of the two processes are comparable. Specifically, the largest of the plasma timescales in Table 10 which is typically the inverse of the operating frequency, need to be comparable to the mean flow time scale for the above to be true. However, for the low Reynolds number in the incompressible limit we are interested in the mean flow responds much slower than the discharge evolution. As an example, the typical fluid flow time-scale for air over a 20 cm chord airfoil at Re=100 is of the order of a few seconds while the discharge operates at a few KHz frequency. Hence, a decoupled treatment for the glow discharge with the fluid flow solver would be a reasonable approximation here. In this chapter, a brief description of the body force treatment of the glow discharge plasma will be presented. This is followed by a detailed presentation of the phenomenological body force model and the resulting glow discharge-induced flow effects. Here, an empirical formulation of the body force based on an assumed plasma structure was used to study the resulting flow patterns and the aerodynamics which revealed qu alitatively similar characteristics to the experiments. The next section presents the results obtained from the fully self-consistent

PAGE 97

82 plasma dynamics calculations performed on th e asymmetric electrode actuator operating on Helium gas using a plasma-fluid model as discussed in the previous chapter. The observed results provide some interesting insight into the mechanism of plasma operation, structure and the resulting force ge neration. Also, the effect of the various input voltage waveforms, operating parame ters (the applied volta ge and frequency), electrode spacing and arrangements are studied. 4.1 Body Force Formulation It is worth recalling that the momentum im parted to the fluid by the discharge is through Lorentzian collisions of the heavier ion species with the neutral particles. We will now revisit the right hand side of the ion momentum equation (2.21) which represents the collision term. Although ther e will be different ki nd of collisions, the primary collision loss will be that between the heavier ions and the neutral particles. If we assume that the drift-diffusion approximation (2.25) for the momentum equation is valid, then the net loss due to collision becomes (,)iiiiiiimicollnqEnkTvnmFxt (3.1) Here, the first term is the contribution from the electric drift and the second is the gradient of the thermodynamic pressure which reduces to more of a diffusion flux if the thermal effects are minimal. If we neglect th e diffusion flux, then we get the loss due to collision balanced by the Lorentz force acti ng on the charged ions. The electrons being lighter are not significant pl ayers in terms of momentum transfer. Consequently the electrons reach mean kinetic energy levels higher than those of the heavy components. This loss due to collision Fcoll can be viewed as the instantaneous local body force acting on the neutral fluid particles and can be mode led as a source term in the Navier-Stokes

PAGE 98

83 equation. Once we have the instantaneous distribution of the ion number density ni and the electric field E, the instantaneous body force Fcoll can be calculated as in (3.1). This force acts on a much different time-scale than the characteristic time of the flow and it is important that it be appropriately averaged over time. If we represent the integration timescale by and the time-averaged collision force by Ftave, we have 0,,taveFxFxtdt (3.2) Figure 12 Schematic illustration of the plasma-fluid flow modeling framework. The choice of this time-scale is determined by the smallest of the flow scales which need to be resolved so that the corres ponding effect of the for ce in that time scale can be realized. It is worth noting that if is much larger than the time period of the operating voltage cycle, then the averaged force is independent of time. In other words, Dielectric barrier discharge modeling gives ni ( ion number density) and ne (electron number density) Calculate instantaneous force F=E(x,y,t)e(ni-ne) Average F over a flow time scale to get Ftave,x Solve flow field x tave xy xxF y uv x p u t u, 2 Substitute body force in Navier-Stokes, eg: x-momentum eqn.

PAGE 99

84 the fluid will see the discharge as a source of constant average body force in time. However, if the flow reacts faster than the time period of the driving voltage cycle (which happens to be the largest plasma time scal e) then the unsteady forcing effects of the plasma body force will be felt by the flow and hence plasma transients need to be modeled. This is important while modeli ng subsonic compressible flows or turbulent flows with low plasma frequency. A schematic illustration of the modeling approach is presented in Figure 12. 4.2 Navier-Stokes Equations with Body Force The numerical model consists of the c ontinuity equation and the two dimensional momentum equations for a steady incompre ssible viscous flow. The body force terms, which are added to the momentum equations, ca rry the effect of the plasma discharge on the fluid flow. The fluid is assumed to be in compressible in view of the plasma induced jet being a low Reynolds number and essentially isothermal phenomenon. In the following, we offer a complete set of governing equations, in the two-dimensional form. ABC D txy (3.3) Au v (3.4) 2 xx xyu Bup uv (3.5) 2 xy yyv Cuv vp (3.6)

PAGE 100

85 ,0tavex taveyDF F (3.7) where xy is the shear stress and Ftave,x ,Ftave,y are the body force components in a 2-D case, as given in (3.2). 4.3 Modeling a Plasma Actuator Using a Phenomenological Force Model Here we try to model the plasma actuator used by Roth et al.[4] and Corke et al.[25] with dimensions as shown in Figure 13 where Lu, Ll and d are the length of the upper, lower electrodes and gap distances respectively. Th e plasma generated in this instance is weakly ionized and non-thermal (~ room temp erature) where the electrons alone have a large mean kinetic energy. The ions are more or less in thermal equilibrium with the neutral species because of the extensive collisions with neutrals in a weakly ionized state. A phenomenological plasma model without a deta iled simulation is used to depict the glow discharge. The motive of this study wa s to develop a framework which would take into account detailed balances between momentum transfer mechanisms such as convection, hydrostatic pressure, viscous st resses and the body force arising from the Lorentzian collisions. Now we present a brief phenomenological discussion of the plasma generation in an asymmetric electrode conf iguration based on our understanding of the plasma structure in a high frequency parall el plate (Massines et al.[6]) arrangement. 4.3.1 Plasma Generation in an Asymmetric Electrode Configuration At these high frequencies the charged sp ecies (mainly electrons) attain high velocities, leaving a very short time for the species to interact and recombine. The ions, being less mobile are believed to be trappe d between the electrodes without actually reaching the electrode surfaces. This ion-tra pping mechanism is one of the keys for

PAGE 101

86 sustaining the discharge as they ensure adequate seed particles for further breakdown at the beginning of each cycle. This is one of the reasons why at low frequencies (typically <100Hz) when there is enough time to charge the electrodes, the discharge is no more periodic and starts resembling the glows in a DC case. The defining characteristic of the arrangement shown in Figure 13 is the asymmetry which causes the electric field lines to be curved with dominant fringe effects. This distorts the electric field distribution which is far from being uniform and can modify th e plasma species dynamics sufficiently. The electric field distribution is st rongest in the region closest to the inner edges of the two electrodes which will also be the region of dense ionization. This strength decreases in magnitude as one moves far away from the electrodes. Figure 13 Schematic of the actuator arrangement with approximate shape of electric field lines. The thickness of th e electrodes d is 0.1 mm It was proposed that asymmetry of the anode and cathode conf iguration as seen by the fluid is responsible for the uni-directi onal flow generation. Let us look at a single cycle of the discharge. When the voltage ris es, the electrons leave the negative electrode (cathode and upper electrode for this half cy cle) surface and move towards the anode and in the process are involved in collision, ionization etc. Further, the anode surface (lower electrode) here is not the metal but the dielectric surface which is in contact with the fluid. We can call this as the pseudo anode. Based on Figure 13, the surface area of this Lu LL d

PAGE 102

87 pseudo anode is much larger than that of th e cathode. Also the distribution of charges on the dielectric will be different than on an equi valent electrode conductor, as the dielectric prevents movement of charges on the surface. As a result, we have a concentration of charges on the cathode (upper electrode) and a less dense distribution on the dielectric surface, caused both by the asymmetry in the electrode structure and material. Thus we have a much stronger electric field near th e cathode and comparatively weak field near the anode. On voltage reversal the memory voltage due to the space charge separation aids recombination which appears as the resi dual current. Our view is that the weaker electric field near the pseudo cathode (pseudo anode during the previous half cycle) renders the initial accelerating force aiding recombination small when compared to the corresponding force during the previous half cycl e. This, in turn results in a lower energy level in the region of the plasma when compared with the first half cycle. Since, in the recombination phase, the lack of electron em ission will lead to species depopulation. The subsequent breakdown during this half cycl e will depend on the presence of a critical number of electrons with th e threshold energy. At low frequency of operation this is usually violated and hence there is little or nor breakdown at all. In any case, the discharge in this half cycle is weaker than the corresponding other half. This nonuniform plasma generation coupled with the high-frequency operation can hypothetically explain the thrust produced in one direction alone. It should be stressed that while the electrons are of much higher ve locity than ions, due to thei r small mass their ability to transfer momentum is insignificant. Instead, between the plasma and the ambient fluid, momentum transfer is realized via ions. As can be seen from(3.1), the transferred momentum is a function of both the number density distribution a nd the electric field

PAGE 103

88 strength. In the present modeling effort we will use the electric field as the key indicator to represent the body force and just cons ider a uniform species number density distribution. In this study we have incorporated the geometric effect of the electrodes and the resulting impact on the paraelectric force in the computational model. In the later part of this ch apter, we will look at a more detailed first principlesbased modeling study to further our understanding of the momentum generation mechanism. 4.3.2 Linearized Electric Body Force The objective of the following formulation is to sketch a qualitative expression for the body force term acted by the plasma and on the fluid. As one can see from Figure 13, the electric field lines, as show n are concentrated at the catho de and are almost uniformly distributed on the anode. It is reasonable to simply represent the field lines to be parallel in most of the region except the small space near the cathode. Thus as shown in Figure 14, we can linearize the field variation in space without computing the detailed electric field. Specifically, the field lines here are such that the strength of the field decreases as one move far away from the source. This vari ation of E can be mathematically written a y k x k E E2 1 0 (3.8) Where E0 is the electric field in the darkened region in Figure 13 In general, the constant E0 is large and, k1 and k2 are two positive constants which represent the gradient of electric field intensity along the two mutually perpendicular directions, namely x and y. The positive nature of these two constants ensures that the electric field intensity decreases as one move along the positive directions of the axes. E0 can be approximated as

PAGE 104

89 0U E d (3.9) where d (Figure 13) is the distance of sepa ration between the two electrodes in the x direction. The constants are evaluated by using the condition that the field strength is the breakdown value at the plasma fl uid boundary. Hence we can find k1 and k2. Figure 14 The line A-B constitutes the plasma fluid boundary using linear approximation. The electric field strength outside this line is not strong enough to ionize the air. The components of the electric field are given by 2 2 2 1 2k k Ek Ex (3.10) 2 2 2 1 1k k Ek Ey (3.11) The period of interest during a part icular cycle is the small time t during which the plasma discharge takes place. As it is here th at we have the high discharge currents and the resulting momentum transfer. The magnitude of this t is small compared to the fluid time scales, especially so for the frequenc y range we are dealing with. Thus the body force components along the xand y-direction, fx and fy, are calculated as y A a B x b

PAGE 105

90 xxcc f Ee (3.12) yycc f Ee (3.13) This force f, acting on the charged particles is imparted to the neutral molecules. Here, c is the charge species number density and ec, the electronic charge. The domain in which the body force acts is determined by the value of E. We use the function in the expression for the body force. 1 for crE E ; 0 for crE E ; The Ecr in this case is the break down electric field strength, Eb and represents the voltage boundary. This force acts only during the time t (during which the plasma is formed) and that corresponds to one half cycle, as already discussed. The force during the second half of the cycle can be neglected as it ha s little or no plasma formation. The high frequency of the discharge (around 5 KHz) ma kes it reasonable to consider the force acting on the fluid as a constant. Hence the force can be time averaged over a complete cycle although it exists only for a small time t per cycle. x tavexx aft F ft T (3.14) y taveyy aft F ft T (3.15) This time period of the cycle Ta is then period of the applied voltage and (reciprocal of the time period) is the frequency of the applied voltage. These two components of the body force can be added to the Navier-Stoke s equation. We can infer from the above expressions that, when the frequency is zero, the body force term vanishes which is

PAGE 106

91 consistent the expectations. The‘t’ used above can be e xpected to vary with But this is not the case, as the discharge cannot continue. The memory voltage bounds the t in a small range independent of the frequency. This explains the lack of induced flow in the case of DC voltage generated plasma. In the present model, no attempt is made to establish a criterion for the plasma is capable of inducing the fluid flow. Such a criterion requires a more detailed modeling of the plas ma dynamics which is presented in section 4.4. The force calculated in (3.14) and (3.15) can be substituted in the Navier-Stokes equations (3.3)-(3.7). 4.3.3 Problem Description In order to model the plasma induced flow effects a single actuator configuration on flat plate as shown in Figure 13 is included as one boundary of the computational domain as shown in Figure 15. We will only focus on the vicinity of the electrode to investigate the detailed flow field. The fl at plate is 20.5 mm long, with the upper electrode fitted 12mm from the leading edge of the plate. The electrode has the length 0.5 mm and the height 0.1 mm. The height of the domain is 10 mm. The flat plate constitutes the lower boundary of the computational domai n, with its leading edge 1mm from the vertical boundary of the grid. As shown in Figure 15, three grid blocks constitute the computational grid. The total number of grid points employed is 290x100. The chosen mesh distribution is clustered to focus th e computational effort on the flow phenomena occurring around and after the upper electrode wh ere the body force is applied. This grid size and distribution can pr oduce satisfactory solutions in the present study. The key dimensionless parameters in the present model are: (i) Reynolds (Reynolds) number:

PAGE 107

92 m m m mL V (ii) The ratio of the body force to the inertial force is the same: e e e e m m m mF L V F L V where L is the characteristic length; V is the free stream velocity; F is the body force; is the density; is the coefficient of viscosity; The subscripts m and e represent the model and experimental values respectively. For the CFD simulations, we have employed the pressure-based algorithm detailed in [112-113]. The convection scheme we have used is the second order upwind scheme [113]. Boundary condition: The no-slip condition when applied at the solid wall yields, 0 v u at 0 y (3.16) At the outlets, zero velocity gradient is assigned taking advantage of the negligible influence of the plasma, far downstream of the electrode. The plasma discharge boundary is specified by 0 c for bE E (3.17)

PAGE 108

93 Figure 15 The computational domain used for cal culations shows the flat plate with the upper electrode. The various dimensions, the boundary specifications and the stations used for analyzing results are indicated. 4.3.4 Effect on Flow Structure The various modeling parameters used were the same as that used in the experimental studies [4]. As against experiments which used multiple actuators, only a single actuator was used. Further, heat transf er computations were also performed in which the contact flat plate in Figure 15 is maintained at twice the ambient fluid temperature for the present discussion. A ll temperatures are non-dimensionalized. The various parameters used are summarized below. = Frequency of applied voltage = 3 kHz c = Charge density (electrons) = 1.01011 / cm3 Ua = Applied voltage = 4 kV rms Eb = Breakdown electric field strength = 30 kV / cm [1] t = Discharge time = 67 s (obtained from Massines [6]) a = Height of the plasma = 1.5 mm (par ameters used to determine the linearized electric field distribution) b = Width of the plasma = 3 mm (obtained from picture) d=distance between the plates=0.25mm (mm) OUTLET (Zero pressure gradient) INLET ST 1 (4.75) ST2(12.0) ST3(13.9) ST4(17.3) OUTLET SOLID (no slip) 1.0 12.0 0.5 8.0 9.9 0.1 electrode

PAGE 109

94 Lu=length of upper electrode is 0.5 mm Ll=length of upper electrode is 3.0 mm An issue of note here is the choice of c. The value used above is a rough estimate for air based on the data from literature and will vary depending on the gas employed. Specifically for the helium gas we use in sec tion 4.4, the number density is of the order 109/cm3. Consequently the force calculated should be expectedly weaker. The values of the various parameters for the different cases in the study are listed below: Free stream velocity, V: 2, 4, 5 and 10 m/s Applied voltage, Ua: 3, 4 and 5 kV rms Applied frequency, : 2, 3, 4 and 6 kHz The corresponding Reynolds numbers (laminar flow) based on the above combinations and the characteristic length, Ll = 3 mm (The length of the lower electrode.)are: 411, 822, 1027 and 2054. To study the evolution of the velocity profile and temperature profile four different stations (Figure 15) ha ve been chosen for discussion purposes. In all the results, the reference values for the operating parameters are 4KV rms for the voltage at 3 KHz with a mean flow of 5 m/s whic h are maintained unless otherwise stated. The effect of the plasma operation on the various regions of the flow can be seen from the velocity profiles at the various sta tions is shown in Figure 16. As expected, the maximum peak value of the velocity is obtai ned downstream of the electrode (i.e. in the high pressure region). The overshoot of velocity above the free-stream velocity indicates the performance index of the actuator. Overal l, the velocity structure in the downstream of the electrode resembles that of a wall je t. With the mass being conserved at each vertical plane, the free stream velocity upstre am of the electrode does not attain the farfield value, in trying to compensate for th e induced mass flux in the wall jet region. The

PAGE 110

95 free stream approaches the far-field value asym ptotically as one move far away from the electrode. Figure 16 Velocity profiles compared at the four different stations In the following studies, the parameters which are maintained constant are assumed the reference values of 4KV rms, 3 KHz and 5 m/s for the applied voltage, frequency and inlet velocity respectively. In Figure 17a comparative study of the induced jet velocity profiles for different operating parameters is presented. Figure 17a gives the normalized velocity profile at ST4 for the various free stream velocities. As expected, the plasma actuator is most effective at lower Reynolds numbers. Figure 17b and Figure 17c give the normalized profiles for the different applie d voltages and frequencies. The peak jet velocity and the body force are both direct func tions of the frequency and applied voltage in the linearized body force model and ther efore increase with the actuator operating parameters. It is worth mentioning that with t being held constant, the frequency variation becomes a load factor study. Ho wever, detailed plas ma modeling studies indicate significantly affected plasma struct ure at different frequencies and this might play a more serious role than just a load f actor. In Figure 18, the wall shear profile at ST4 is presented. The wall shear is normalized by the corresponding peak value for the no-

PAGE 111

96 plasma case. Closer to the wall, the wall jet flow reduces the shear stress to a large negative value and farther away the velocities gr adually reduce to the free stream value. Figure 17 Effect of the various parameters as seen form the velocity profiles at ST4 Figure 18 Effect of the body force on the wall shear at ST4 (b) Effect of the frequency of applied voltage. (a) Effect of free stream ve locity on velocity profiles.

PAGE 112

97 4.3.5 Effect on Heat Transfer In addition to the flow field, the effect of the induced jet on heat transfer from the surface was also studied. Figure 19 presents the temperature profile along the wall at the different stations. In Figure 19a, the temper ature gradient along the normal at ST2, ST3 and ST4 is smaller than that at ST1. This is due to the viscous layer being thinner at ST1. The Reynolds number considered in 1027 (bas ed on the lower electrode length) and Prandtl (a) Plasma off (b) plasma on Figure 19 Temperature profile along the wall for Re=1027 and Pr=0.7 number of 0.7. With the plasma on, the te mperature gradients at ST2-ST4 increase but there is no appreciable change at ST1. Figur e 20 presents the heat flux from the surface for different applied voltages and frequencie s. In Figure 20a and Figure 20b, the effect of the plasma is felt only from a small region upstream of the electrode which is the region near ST2. All heat flux values were normalized using the corresponding plasma off value. The maximum heat flux in the ra nge of frequencies and voltages studied was about four-five times the plasma off value, at a point immediately downstream of the electrode.

PAGE 113

98 (a) Comparison of Normalized heat flux for different applied voltages (b) Comparison of Normalized heat flux for di fferent frequencies of applied voltage Figure 20 Comparison of normalized heat flux for different voltages and frequencies 4.3.6. Thrust Estimation In order to assess the accuracy of the phenomenological body force model developed in the previous sections, we try to compare the thrust produced with the experimental studies [26]. A particular issue pertaining to the model is the procedure for employing the experimental data to fix up the modeling parameters such as the duty cycle, charge density etc. In the followi ng, we present the force prediction under the quiescent condition, in accordance with th e experiment conducted by Van Dyken et al.[26]. In these experiments, the force measur ed is the shear caused by the induced flow on the contact surface in a quiescent envir onment. After incorporating the duty cycle

PAGE 114

99 factor, the force acting on a fluid element from Eqs. (3.12)-(3.15) and acting on the volume is xccxt F eEAdx T (3.18) where A is the area of the face of the differe ntial volume perpendicular to the orientation of dx and Ex is the component of the electric field in the x-direction. Under the linear field approximation, the net force in the x-direction, Fx, acting on the contact surface is obtained by integrating along the discha rge operating region b is given by, 0 b xxcct FEneAdx T = 01()b xcc ot EkxneAdx T (3.19) where nx is the direction cosine, illustrated in Figure 14. Here, k2 is neglected as y=0 on the surface. From the Maxwell equation we get 1 0 ccdEe k dx (3.20) The absolute values are taken so that we can conveniently ignore the sign conventions while defining the charge on an electron. Usin g the condition that the electric field is equal to the breakdown field at the edge of the plasma along the surf ace of the plate we have from (3.20) 0 0()b c c E E b e (3.21) The equation (3.20) when substituted in the expression for the force in (3.19) gives, 0110 0()b xxt F EkxnkAdx T (3.22) Substituting for k1 and evaluating the integral we get

PAGE 115

100 22 0 02b xx E Et FnA T (3.23) Note that, by assuming the linear variation for E, we are fixing the direction cosine nx from the choice of the breakdown lengths a and b. These two equations (3.21) and (3.23) relate the three parameters, namely E0,c and b with the thrust. The duty cycle also vari es with the voltage but is assumed constant in the present discussion. All other parameters in the expression for the thrust in (3.23) are fixed within the modeling framework. For a start, let us assume the knowledge of the thrust F and the plasma breakdown length b fo r a particular value of the voltage. This helps us to calculate the charge number density c from (3.21) and (3.23). In order to extend the knowledge of two parameters, c and b, we assume a linear proportionality for above parameters with the applied voltage. The experimental results by Van Dyken et al. [26] present the force generated for various voltages with an electrode 0.625mm wide and Teflon dielectric. Lets us start with a volta ge V=4.66KV, for which, from the direct measurement, F=0.3969mN and from light inte nsity measurements given in [26] one deduces that b=2.8 mm for the present electro de configuration. The various parameters that are assumed invariant in the present framework are 12 0 620.77 8.85210/ 1.1510 0.89xt T Faradm Am n The above duty cycle for a sine waveform of the applied voltage has been measured experimentally by Van Dyken et al.[26]. By substituting the above data into (3.21) we get

PAGE 116

101 c=1.59 x 1011/cm3. As we can see, this calculated nu mber density is roughly of the same order as the value assumed earlier. Now, th e charge number density and the length b are fixed for the assumed voltage. The values of the same for any other voltage can be scaled as: 11 22 c cV V and 11 22bV bV These values can be substituted back (3. 21) and (3.23) to obtain the force. The comparison between the experimentally reco rded force and the calculated force is presented in Figure 21. Figure 21 Comparison of the experimental a nd calculated plasma force for different voltages Here the Factual represents experimentally measured force while the Fcalculated represents the predicted values. While there is a quantita tive deviation between the experimental and computational values, the trends between them are consistent. The different results likely arise from linearization of the electric fiel d as well as simplification made during the scaling process described above.

PAGE 117

102 4.3.7 Effect of Plasma on Aerodynamics The effect of the plasma operation on aerodynamic performance was investigated using the same model as above for the plasma modeling. Corke et al.[25,27] demonstrated that the plasma operation enhances lift as well as drag due to increased shear as well as possible downstream separation of the induced wall jet. In this study, computations were performed in consultati on with available experimental evidence [2122,25] to estimate the generated thrust as we ll as the aerodynamics. The computational model consists of a flat plate placed in a free stream of fl uid as shown in Figure 22. The electrode arrangement serving as an actuator is placed at the middle of the plate. The computational domain encloses the flat plate and electrode arrangement. The boundary conditions for the domain are as shown in Figure 22. To reduce the computational domain and to attain the desirable accuracy, th e plate is chosen to be 1cm long, so that the upper electrode, 0.5mm long and 0.25mm th ick, can be adequately modeled without excessive computational resources. The main interest here is to probe the effect of the glow discharge-induced fluid flow on the existing flow structures, for both attached and separated fluid flows. The number of grid points used for the simulation is 161x102. The grid distribution is focused towards capturing the modified flow field in the vicinity of the plasma. The grid points are more clustered near the plate and near the leading edge. The numerical simulation employs a mulitblock model solved by a second-order upwind scheme. Detailed discussions can be found in Sh yy et al.[112-113]. It is to be noted that, wherever not mentioned, the values for th e voltage, frequency and Reynolds number are 4KV (rms), 3KHz and 333 respectively. The values of the various parameters used are Free stream velocity, V: 0.2, 0.4, 0.5m/s Reynolds number based on plate length of 0.01m: 133, 266,333

PAGE 118

103 Angle of attack: 0, 8 and 20 degrees Figure 22 A schematic of the computational dom ain with the flat plate and the electrode. The electrode has negligible thickness in the above arrangement (a) Re=133 (b) Re=266 (c) Re=333 Figure 23 Drag variation with angle of attack for various Re. The plasma is operated at 4KV rms and 3KHz Fluid flow: The effect of the plasma on various flows at different Reynolds numbers and angles of attack are studied. In the presence of the plasma, the momentum loss will Electrode

PAGE 119

104 include the glow discharge-induced body force effect. Figure 23 presents the variation of the drag coefficient with angle of attack for a given Reynolds number in the presence of the plasma. Irrespective of Re, the plasma induced additional drag, which corroborates with the results reported in [25,27]. More importantly, the magnitude of increase in the drag coefficient CD appears insensitive to angles of at tack for a particular Re for the two lower angles of attack (0 and 8 degrees) when the flow is attached. For a 20-degree angle of attack where the flow separates, the plas ma effect modifies the separation region and the form drag. The decreased impact of the pl asma force for higher Re has been explicitly pointed out by Shyy et al.[9]. The body force term used in the momentum equation is independent of Re and hence becomes smaller when normalized by the dynamic head. We can calculate the increase in drag forceD defined by (1.24) in stead of the drag coefficient, for a given angle of attack. D = 21 2DCU (1.24) Table 13 Change in drag coefficient and drag force with Reynolds number Re 133266333 D C 2.80.710.46 21 () 2 D DNUC 6.726.8166.9 Table 14 Change in lift coefficient and lift force with Reynolds number Re(AoA=8) 133266333 lC 1.510.4350.297 21 () 2lLNuC 3.6244.1764.455 It is observed in Table 14 that the drag force is largely independent of the Reynolds number in the absence of flow sepa ration indicating that, the drag increase is

PAGE 120

105 mainly an artifact of the large plasma force. Similar trends were observed in the experimental results presented by Corke et al .[25]. Figure 24 presents the variation of lift coefficient with Re for a given angle of atta ck. In general, the plasma on case yields a higher lift. In Figure 24 and Tabl e 14, the increase in lift coefficient becomes smaller as Re increases. However, for separated flow s as shown in Figure 24b, there are more substantial increases in lift compared to th e attached flow (Figure 24a) with the same glow discharge parameters and Reynolds number. As indicated in Table 14, and consistent with the drag force, the increase in lift force appears to stabilize with larger Re. Heat transfer: The heat transfer enhancement result ing from glow discharge has been briefly discussed by Shyy et al.[9]. In the pr esent study the heat transfer is studied for various Re and angles of attack. Figure 25 presents the Nusselt number variation with the angle of attack for various Re. The Nusselt number here is given by H l Nu K T (3.25) where H is the heat flux, l is the characteristic length of the plate, K is the coefficient of thermal conductivity and T is the temperature difference between the ambient fluid and the hot surface, which is set to a non-dimensi onal value of unity. In all cases, the heat transfer is enhanced during plasma operation and especially so when the angle of attack is high and flow separates. Moreover, the gl ow discharge-enhanced Nusselt number only exhibits modest variations versus the angle of attack. The plasma helps reattach the flow and enhances heat transfer. The modest devia tion in Nu for the high angle of attack cases is due to the incomplete removal of the re-circulation zone. One way to further improve the heat transfer rate is to place multiple el ectrodes along the stream wise direction to induce stronger fluid flows and also eliminate the re-circulation zone.

PAGE 121

106 (a) Angle of attack= 8 degrees (b) Angle of attack=20 degrees Figure 24 Variation of lift with Re for different angles of attack. The plasma is operated at 4KV rms and 3 KHz. (a) Re=133 (b) Re=266 (c ) Re=333 Figure 25 Nusselt number variation with angle of attack for various Reynolds number. The plasma is operated at 4KV rms and 3KHz frequency Effect of electrode position on the flow: To investigate the effect of the electrode placement position on the flow, two positions (middl e of the plate and at a distance of 10

PAGE 122

107 percent length (of the plate) from the l eading edge) were considered. The study was conducted for an angle of attack of 20-degrees and Re=333. Shown in Table 15 is a comparison of the aerodynamic and heat transfer data. It appears that the effect of the electrode position on the drag is insignificant, as the observed data seem to be an artifact of the plasma body force. As summarized in Table 15, there are only minor changes in aer odynamics and heat transfer coefficients. Table 15 Effect of electrode position on aerodynamic and heat transfer properties Figure 26 compares the streamlines of the flow over a flat plate, with Re=333 and an angle of attack of 20 degrees. Figure 26a representing the no plasma case shows a large flow separation. In Figure 26b the electr ode is placed at the middle of the plate, while in Figure 26c, the electrode is closer to the leading edge. While the upstream placed electrode can more effectively elim inate the flow separation, both electrode locations produce very comparable outcome in aerodynamics and heat transfer. It seems that the overall impact of the glow discharge, namely, the induced wall jet is not sensitive to the details of the local flow structures fo r the Re considered, making the present device very attractive for active control. Similar, experimental observations have also been reported by Corke et al.[24-25]. In fact Corke et al, report efficient separation control at angle of attacks well past the stall angle. Also a typical pl asma actuator with capacities Electrode position Aerodynamic data Middle 10%from leading edge No plasma Drag coefficient 1.01.0 0.124 Lift coefficient 1.221.12 0.636 Nusselt number 28.3529.57 18.04

PAGE 123

108 described above leads to substantial pressure recovery and also manifold increase in the lift to drag ratio. (a) Plasma off (b) Plasma on with electrode at 0.5l. (c ) Plasma on with electrode at 0.1L from the leading edge Figure 26 Streamlines showing flow separation for Re=333 and an angle of attack twenty degrees for different positions of the electrode 4.4 Self-Consistent Fluid-Model-Based Plasma Dynamics The key to modeling DBD effects in fluid dynamics is to achieve realistic distributions of the species number densitie s and their momentum in the domain which interacts with the neutral fluid by solving the plasma-fluid equations along with the Poisson equation for the potential (Eqs. (2.39) -(2.45)). We will use a simplified form of the plasma model as present it below for convenience. Continuity equation:

PAGE 124

109 Electrons: .e eeeieien nvnSrnn t (3.26) Ions: .i iieieien nvnSrnn t (3.27) Here, the source terms on the right hand side represent the reaction/ionization processes which result in the creation (S) or de struction (r) of the species, as applicable. It is very common to write the ionization coefficient S as a function of E/N. Momentum equation: The momentum equation for the species at high pressures can be reduced to the drift diffusion form which ne glects the inertial and unsteady terms and balances the thermodynamic pressure gradient with the drift force and collision terms. Electrons: eeeeeenEnDnv (3.28) Ions: iiiiiinEnDnv (for high-pressure discharges) (3.29) Electric field equation: The electric field E is obtained us ing the solution of the Poisson equation, given by 0.ie denn E (3.30). The solution algorithm employed has been discussed ex tensively in the previous chapter. 4.4.1 Discharge Actuator Setup In this section, we will model the tw o-dimensional radio frequency dielectric barrier discharge (DBD) in helium gas at high pressure in an asymmetric configuration as shown in Figure 27 The electrodes are 2mm in length and the insulated bottom electrode is shifted downstream by 2mm. A similar case has been studied by Roy et al.[15]. The helium discharge is modeled at a pressure of 300 To rr, and a temperature of 300K and is driven by an AC voltage of 1 KV (peak voltage) operating at 5 KHz. The thickness of the

PAGE 125

110 electrodes is negligible and the dielectric thickness is 5 mm. The electron temperature is assumed to be 1eV (~11600K), while the ions a nd the neutrals are essentially in thermal equilibrium at 300K. The various transpor t properties and property relationships are available in literature and we will use the ones used by Roy et al.[15] for the present study. These are summarized in Table 16. X Y 0 0.002 0.004 0.006 0.008 0.01 0 0.002 0.004 0.006 0.008 0.01 sin( t) Insulateddomain d=4.5 Plasma+fluiddomain d=1.055I N L E TOUTLETO U T L E TINTERFACE X Y 0 0.002 0.004 0.006 0.008 0.0 1 0.006 0.008 0.01SatisfyCurent continuityd Ni e/ d n = 0dNi,e/dn=0d Ni e/ d n = 0PlasmaboundaryconditionsNe-specify thermalflux, dNi/dn=0 Figure 27 A representative 2-D asymmetric discharge arrangements (a) Complete domain with insulator indicating boundary conditions for fluid flow (b) Top half of the domain indicating boundar y conditions for plasma species and the region where the results are presented.

PAGE 126

111 The computational domain employed is a square domain of size 1cm x 1 cm x 1 cm (unit cm in the normal or spanwise direc tion). The grid consists of 127 x 61 points for the 2-D case as is shown in Figure 27b. The dielectric cons tant in the fluid/discharge domain is 1.0055; the permittivity of vacuum and that of the insulator is 4.5. For investigation purposes two different freque ncies are considered, namely, 5 KHz and 10 KHz, respectively, with all other conditions remaining constant. A global time-step of 108 s is used for the computations. The initial number density in the plasma is 1.0x1015 /m3 for all the different species. Table 16 Summary of property models employed for the He discharge simulation. Transport/reaction properties Models/values employed i ( ion mobility) = 3 3211810 1810/ EpcmVs p for 11/25 EpVcmtorr = 4 211 1.54.11027.44 1 / / cmVs pEp Ep for 11/25 EpVcmtorr e (electron mobility) = 211eene cmVs m where 1210/en s ieS (species ionization model) = 1 0.414 4.4exp /epEs Ep r(recombination coefficient) =209/231.0910/eTnms iD(ion diffusivity) =2500/ cms eD(electron diffusivity) = 2/e ekT cms e (viscosity of He gas) =522.010/ Nsm Note: Data presented is the same used in [15]. 4.4.1.1 Boundary condition The boundary conditions are schematically presented in Figure 27b. The homogeneous Neumann boundary condition is app lied for the electrostatic potential at the open boundaries while the Dirichlet boundary condition is used at the electrode:

PAGE 127

112 At the exposed electrode: 0sin2 f t, where 0 is the peak value. At the submerged electrode: 0 For the plasma species modeling, the domain boundary away from the insulator/electrode surface is assigned a zero gradient condition assuming insignificant impact far away from the fluid-actuator in terface. At the dielect ric surface, the drift current and the displacement current from the gas domain is balanced with the displacement current inside the insulator. At the electrode, the treatment is slightly different. The electrons are assumed to be isothermal (at 11600K or 1eV) at boundaries. At the exposed electrode, the thermal flux towards the wall is considered while it is neglected when the drift is away from the wall. For the ions or the heavier species, the drift effects are significant and he nce a zero gradient condition ( 0in n ) is applied. The boundary conditions are summarized in Table 17. 4.4.1.2 Plasma species initial condition The initial condition used to start the simulation requires specifying the starting number densities for the various species. Ideally the simulation should start from a very small charge concentration, but to accelerate th e evolution of the transient behavior to a periodic steady state solution we will use a weakly preionized neutral concentration of 109/ cm3. To look into the sensitivity of the so lution to the initial guess we look at two initial guesses of 1010, 2.5*1010 cm-3. We will us a high frequency of 20 KHz and 1 KV acting on the geometry in Figure 27a for the present discussion. Solutions are compared after 70 time-cycles at which time a periodic st eady state has been attained as shown in Figure 28. There is a significant differen ce in the domain averaged number-density values, but the net charge densities and the resu lting force field seem to be much closer to

PAGE 128

113 each other for the two initial guesses as the solution appro aches a periodic steady state solution. This shows that the dynamics of th e operating discharge plasma actuator is only a function of the difference in the values of Ni-Ne or more generally the net charge in the region. The dependence of Ni and Ne on the in itial guess is simply due to the non-linear nature of the recombination term in the transport equation (Eq (3.26),(3.27)). It was observed that the periodic steady state solution was attained much quicker with a starting guess of 109/cm3 and was used as the initial condition for the results presented in the rest of the discussion. Table 17 Summary of boundary conditi ons for the different variables Variable Open boundary El ectrode/dielectric surface potential Zero gradient At the submerged electrode: 0 At the exposed electrode: 0sin2 f t,01 K V The dielectric is not a bo undary for the potential. Electron,Ion species Zero gradient Electrode: Electrons: Away from the electrode: flux=0 If drift is toward the electrode: flux=NeVeth where Veth is the thermal flux Ions: Away from the electrode: flux=0 If drift is toward the electrode: Zero gradient On the dielectric: Allow surface charge accumulation. To do this the current continuity is established. i.e. At the dielectric-gas interface th e species flux is determined by 0 gasgas dielectricdielectric iiee gasE E e NVNV tt 4.4.2 Plasma Structure In the following discussion we will anal yze the discharge structure in space and time with reference to the experimental observations. The experimental results which we will present here are for qualitative comparison purposes and were run with air as the

PAGE 129

114 working gas using an applied voltage of 6KV and 2.5 KHz frequency. The numerical simulation shown is for 1KV and 5 KHz app lied signal for the configuration described above in Figure 27b and is based on Helium as the working gas. This discussion will shed some light into the operating mechanism of the actuator. (a) Number density (cm-3) comparison after 70 time cycles (b) Ni-Ne (cm-3) and Fx comparison after 70 time-cycles (c) Time history of Fx integrated over the domain Figure 28 Sensitivity to initial condition (20KHz, 1KV applied voltage)

PAGE 130

115 Figure 29 Photo intentisty measurements of En loe et al.[21](pictures reproduced for this discussion) in trying to depict the obs erved discharge evolution over a timecycle. Figure 30 Time dependent plasma behavior (5 KHz, 1KV applied voltage) (d) absolute value of discharge current calculated (c ) Fx (domain in tegrated) over one time-cycle. (a) Input voltage signal and Domainaveraged ion charge density (b) Domain-averaged species charge density and net charge density variation over a time period (b) Current and light measurements show discharge asymmetry (a ) Applied voltage and discharge light intensity measurements

PAGE 131

116 Figure 29 shows the experimentally observe d light intensity measurements of Enloe et al.[21] to depict the observed discharge evolution over time for an asymmetric plasma actuator. While the light intensity measurements cannot be directly compared to the numerically calculated quantities, we wi ll comment on some of the key qualitative features in the observed results. (a) Domain-avera ged ion number density (per cu.cm) over time. Instant A Instant B Instant C Instant D X Y 0 0.002 0.004 0.006 0.008 0.01 0.006 0.008 0.01 Instant E Instant F Figure 31 Ion density evolution in the domain over time (5 KHz, 1KV applied voltage)

PAGE 132

117 Instant G Instant H Instant I Instant J (b) Instantaneous snapshots of th e ion density at different instants over a cycle. Figure 31. Continued. t/T(non-dimensionaltime) Fx(N/cu.m) 0 0.2 0.4 0.6 0.8 1 -1 0 1 2C I B D EF G H A J (a) Domain averaged force (N/cu.m) variation in time Instant A Instant B Figure 32 Force evolution over one time cycle (5 KHz, 1KV applied voltage)

PAGE 133

118 Instant C Instant D Instant E Instant F Instant G Instant H Instant I Instant J (b) Instantaneous snapshots of the fo rce field at different in stants over a cycle. Figure 32. Continued. Self-limiting Discharge: It commonly believed that the breakdown is quenched and controlled by the accumulated charge since the discharge terminates on a dielectric surface and hence the self-limiting nature of the discharge. This can be observed from the experimental results in Figure 29a where a continuously varying applied voltage is necessary to sustain the discharge. Here, the lig ht intensity observed is zero as soon as the observed potential reaches either maxima or minima. This phenomenon can be observed

PAGE 134

119 in the numerical results as well. Figure 30a a nd b show the mean density variation in time for the different species and the net charge de nsity. It is known that in the positive halfcycle the ions move away from the electrode and tend to get deposited on the surface. Hence the decrease in the average ion density over the domain as the averaging process doesn’t account for the surface deposition. Sim ilar to the experimental observation this decrease in ion density is stopped when the in crease in voltage is zero near the maxima. This can be clearly seen in the evolution of the ion density in the domain as shown in Figure Specifically in the first quarter of the time-cycle as can be seen in instants A,B and C, the ion concentration really accumula tes near the dielectri c surface and the free concentration in the domain at its minimum. At instants D, E and F, the applied voltage starts decreasing causing the built up potential from the accumulation to dominate resulting in releasing the ion species into the bulk region. The ion density reaches its maximum value at instant F. At this point th e change in polarity e nhances the release of ions which travel towards the negative exposed electrode. However, there is no ion buildup during this part of the cycle as the elec trode being a conductor absorbs the incoming charge flux. At instant G(1.5 ) Figure a shows the ion numbe r density as stabilizing due to the equilibrium attained between the absorption at the electrode and the accumulated ion release from the insulator surface. This behavior continues for the rest of the cycle. Force variation: As seen from the current measurem ents (Figure 29b and [21]) much larger fluctuations are observed for the positive half-cycle. The corresponding domain averaged force variation in time presented in Fi gure a shows that the positive-half cycle is more efficient and has a predominantly positive force, while the less efficient negative half-cycle has both positive and negative com ponents. This is clearly seen from the

PAGE 135

120 evolution of the force field over time as shown in Figure b. It is well known that the bulk of the momentum transfer is achieved by th e heavier ion species ion species in the domain. As can be seen from Figure b, the positive half-cycle corresponding to instants A,B,C,D and E correspond to the accumulati on of ion species on the insulator surface from the bulk of the domain. Since the subm erged electrode is downstream, the bulk of the transport in this half is in the positive direction and to a lesser extent in the negative direction. Hence the strong positive force re gions near the insulator surface downstream of the electrode and a weaker component ups tream of the electrode for instants A,B,C and D as seen in Figure b. Starting from instant E at which time the charge build-up reverses the polarity of the potential in the region. This happens for half a cycle corresponding to instants E and F during which the weak negative force region near the downstream insulator surface exits due to the ion transport towards the upstream electrode. An interesting phenomenon happens at this instant in the cycle. It should be noted that during this half that the electrons are the species which get accumulated on the dielectric surface resulting in the charge bui ld-up. However, the electrons being lighter species travel faster and hence the faster charge build-up and a quicker reversal of polarity in the domain. There is the qualitative similarity between the instants G,H, I and B,C,D in Figure 32b, but the process is only more quicker. The above mechanism renders the force weaker in this negative half, as the discharge is quenched early. This force generation efficiency is of importance here as the slower response time of the fluid will only see the net overall force generation over a few time-cycles and is less sensitive to the fluctuation over time.

PAGE 136

121 Asymmetry in the discharge: The above discussion on the time evolution of the force field indicates two totally different mechanis ms in the two half-cycles. The first halfcycle shows a pure ionization and deposition of the slower ion species on the insulator surface while the second half-cycle shows the exactly similar behavior for the electrons but at a faster rate. This difference in the mobilites of the two species and the geometric asymmetry totally alters the discharge evoluti on in the second half-cycle. This asymmetry can be observed in the experimental measurements as well. The light intensity measurements Figure 29a and Figure 29b show that the intensity measured over the two half-cycles are distinctly different. Compari ng with the voltage waveform, it can be seen that the discharge light emission is irregular in the positive going half-cycle and is more uniform in the negative going half-cycle. It wo uld be useful to establish a correlation between the observed light intensity measurements from the experimental study and the net charge variation shown in Figure 30b. Specifi cally, the strength of the light intensity measurement over the time-cycle follows the same trend as the absolute value of the time rate of change of the net charge (Ni-Ne) in the domain, given by (()) deNiNe dt While this may sound arbitrary at first, further insi ght reveals that this time derivative measure is equivalent to the discharge current from the domain. As Figure 29b from the experimental study shows, the light intensity measure follow the same trend as the magnitude of the discharge current. While not relevant in the present study, the above discussion leads us on to anot her key aspect which will determine nature of the force field and the momentum coupling. This is th e presence of the ne gative ions and other charged species by including much more extensive chemistry as in the case of atmospheric air modeling. Since, the discharge current variation over time is heavily

PAGE 137

122 dependent on the mobility of the different charged species; the presence of negative ions will significantly alter the di scharge evolution in time. Impact of Voltage Waveform: The observed behavior for the harmonic waveform indicated that the positive going half (irregular) had a different discharge characteristic as compared to the negative going half (uniform discharge) leading to the asymmetry which causes the uni-directional momentum coupling with the fluid flow. Also, in the earlier discussion the positive sawtooth waveform has a much faster positive-going and a slower negative going cycle. The negative sawtooth has a faster negative going and slower positive going cycle. -1000 -800 -600 -400 -200 0 200 400 600 800 1000 00.20.40.60.81t/T (non-dimensional time)Applied voltage (V) Figure 33 Asymmetric sawtooth voltage waveforms used in the experimental and numerical studies For both waveforms faster voltage vari ation results in stronger discharge intensities. The positive sawtooth (Figure 34a) waveform has a stronger positive-going phase which results in an irregular discha rge for the stronger peak. Negative sawtooth (Figure 34b) has a weaker positive-going phase and hence the irregular weak discharge is seen for a longer time. (a) Input sawtooth voltage signal (7KV, 5.8 KHz) reprinted from Enloe etal.[ 21 ] (b) Different input voltage signals (1KV, 10 KHz) used in the numerical study

PAGE 138

123 As discussed in the earlier section, the light intensity can be directly correlated to the discharge current or (()) deNiNe dt The ion density (Figure 35a) and the net charge density variation (Figure 35b) are shown in Figure 35. As mentioned, the maximum slopes are observed for the two sawtooth waveforms in the early part of the cycle as compared to the harmonic waveforms are much larger. (a) Positive sawtooth (b) Negative sawtooth Figure 34 Experimental results of Enloe et al.[21] comparing the discharge observed from the two different sawtooth waveform s (pictures reproduced from [21] for this discussion) (a) Ion density for the different waveforms (b) Net charge density Figure 35 Impact of waveform on the species number density Efficiency of the Waveform: Enloe et al. [21] measures the thrust produced by the actuator for both the above sawtooth waveform s using a mass balance ([21]) to assess efficiency of the two half-cycles. The pos itive sawtooth has a faster varying positive

PAGE 139

124 going voltage and a slower varying negative goi ng voltage and vice versa for the negative sawtooth. The experimental thrust measurements in Figure 36 clearly shows that the positive sawtooth is more efficient compared to the negative sawtooth waveform. This is because unlike the light intensity measurements, the force orientation tends to depend heavily on the geometric asymmetry and henc e strongly dependent on the polarity of the exposed and the submerged electrodes. On the other hand the magnitude of the force is a function of the strength of the discharge and charge concentration available in the domain. The positive sawtooth has a stronger po sitive going part which results in a strong breakdown while the negative sawtooth has a stronger negative going part. One can infer from the force response to the harmonic waveform (Figure a) that the positive going part (exposed electrode is positive) is possibly more efficient than the negative going part. The domain-averaged force responses for the sawtooth waveforms are shown in Figure 37. The force for the positive and negative sawtooth waveforms starts the time-cycle in similar fashion with negative slope. For the negative sawtooth, the force continues on its negative slope and attains a negative peak du ring the fast varying early phase (negativegoing voltage) and then steadily increases to a positive value during the slow-varying positive-going phase. Table 18 shows the time-average value of the force in the domain for the three different waveforms. The harmoni c waveform seems to be the most efficient waveform. The negative sawtooth is the least efficient with the positive sawtooth falling in between. The force profile over time indicat es that the positive going phase is the most efficient part of the voltage cycle for the asymmetric actuator modeled in this study.

PAGE 140

125 Figure 36 Thrust measurements of Enloe et al. (reproduced with permission from [22]) showing the positive sawtooth waveform as being more efficient Figure 37 Axial force integrated over the domain with time (Fx (dynes)) for different waveforms (10Khz, 1KV applied voltage) The green vertic al line indicates the instant of change in slope in the sawtooth waveform. Table 18 Domain averaged force over the voltage cycle for different waveforms Waveform Average force over the cycle (dynes) Positive sawtooth 0.0068 Negative Sawtooth 0.0063 Sinusoidal 0.0081 4.4.3 Geometric Effects 4.4.3.1 Impact of electrode spacing The asymmetry observed in the earlier discussion seems to be heavily orchestrated by the asymmetry in the geometry. So to investigate this we looked at three

PAGE 141

126 different geometric constructions as shown below. The discharge was generated by a 20KHz, 1KV sinusoidal applied voltage waveform. Figure 38 3-Different configurations to study the impact of electrode overlap Configuration I has two el ectrode lengths (2*2mm) ho rizontal separation between exposed and submerged electrode centerlin es. Configuration II has a one electrode separation and configuration III has zero separation. Figure 39 Time evolution of species number density over one time-cycle (20KHz, 1KV voltage signal) Based on the observed average force values obtained in Table 19, it looks like configuration II is the most efficient while configuration III is the least efficient. Configuration I and II both display strong as ymmetry in the domain averaged force over the time-cycle. It can be observed that config uration I results in a weaker average force and that there might exist an optimum gap which gives the best performance. This

PAGE 142

127 asymmetry tends towards a more symmetric force variation in time as the electrode gap is reduced. Configuration III being the symmetric arrangement shows a symmetric variation over time. While the force in itself shows a smaller value as compared to the other cases, the symmetry enhances the overall time-averaged force in the domain. Figure 40 Impact of the geometric configuration on the domain integrated force field (20KHz, 1KV voltage signal) Table 19 Domain averaged force over the voltage cycle for different configurations Geometry Average force over the cycle (dynes) Configuration I 0.0071 Configuration II 0.0087 Configuration III 0.0065 4.4.3.2 Impact of lower electrode size Here we fix the length of the top electrode and vary the bottom electrode to gauge the impact of the area cove red the by the net surface discha rge formation. The actuator operates using a 20KHz, 400 V sinusoidal voltage waveform. Figure 42 and Table 20 indicate that the weakest for ce field is obtained for the 2mm electrode case while there is steady increase in the net force value as the size of the bottom electrode increases. This effect is because the larger lower electrode obviously resu lts in a much larger surface discharge generation and hence the stronger force field. However, once the lower electrode becomes sufficiently large, a balance between size and strength will be reached

PAGE 143

128 signifying the most efficient configurati on. Similar trends have been observed by experimental studies such as by Enloe et al.[22]. Figure 41 Three different lower electrode lengt hs modeled to investigate the effect of increased surface area on the discharge Figure 42 Impact of lower electrode size Table 20 Force dependence on the lower electrode size Lower electrode size Average force over the cycle (dynes) 2 mm0.00023 3 mm0.00065 4 mm0.0010 4.4.4 Impact of Applied Voltage and Frequency The frequencies used were 20KHz, 10 KHz and 5 KHz with a voltage of 1KV. The voltages used were 1KV, 500V, 400V and 250 V with a 20KHz frequency.

PAGE 144

129 4.4.4.1 Frequency sensitivity (a) Averaged net charge density (b) Averaged horizontal force component Figure 43 Charge density and force variation over time for different frequencies The 5 KHz frequency case seems to show a much larger fluctuation (Figure 43a) in the net charge (over the full time-cycle) as compared to the higher frequencies. This is because the smaller frequency enables the charge in the domain to completely settle down on the insulator surface before voltage reversal occurs. This surface charge accumulation might not be so complete for the higher frequencies. The incomplete accumulation of the charge on the insulator surface for higher frequencies in the force variation over time is shown in Figure 43b. The 5 KHz frequency case shows a sharp transition of Fx upon voltage reversal indicative of charge release from the insulator surface. On the other hand, the higher freque ncy cases display a slower transition from the negative value in the domain till it stab ilizes to a value near zero. The 5 KHz frequency apparently shows a more effici ent time-cycle as compared to the higher frequencies. But this can be misleading, since the corresponding duty cycle is much larger. In other words, one cycle of the 5 KHz frequency actuator is equivalent to 4 cycles of the 20 KHz actuator in time and hence will have a stronger impact when acting on a slower process.

PAGE 145

130 Table 21 Domain averaged force over the voltage cycle for different frequencies Frequency Average force over the cycle (dynes) 20 KHz0.0059 10 KHz0.0095 5 KHz0.0125 (a) Averaged Net charge density (b) Averaged horizontal force component Figure 44 Charge density and force varia tion over time for different voltages Table 22 Mean domain averaged force over the time cycle for different voltage Voltage Average force over the cycle (dynes) 250 V0.0008 400 V0.0017 500 V0.0026 1000 V0.006 4.4.4.2 Applied voltage sensitivity The force shows a power-law dependence on the applied voltage. i.e. 2.1 F xV In Figure 45 three different force m easurements are compared for their voltage. The experimental (work of Van Dyken et al.[26]) and the phenomenological model data are the same as shows in Figure 21. It can be seen from the figure that all three sets of data corresponding to the plasma -fluid model, empirical approach and the experimental measurements indicate compar able dependence on the applied voltage.

PAGE 146

131 Figure 45 Fx power-law dependence on voltage is compared for two different numerical modelsphenomenological and the first principle-based plasma fluid model along with experimental data of [26]. The arrows poi nt to the set of axes chosen for each of the curves. 4.5. Summary In this chapter, we developed and pr esented a modeling framework to study the operating mechanism and the thermo-fluid poten tial of a dielectric barrier discharge actuator. A phenomenological model was developed as suming a linearized electric field distribution and constant species number densit y in the region. This approach is basically equivalent to solving a lowest order repr esentation of the problem. Although it required experimental characterization to achieve closure, it was able to reasonable capture some of the wall-jet type flow features and othe r observed experimental outcomes. However, the need for extensive empirical data, such as the specie charge density, duty-cycle coeeficient restricts the scope of this model. In order to evolve the capability to the next levela high fidelity, first principles based approach using a plasma-fluid (or plasma-hydrodynamic) continuum formulation

PAGE 147

132 to model the discharge dynamics. The two-species plasma simulations revealed that the generation of the uni-directional momentum coupling is primarily caused by the combination of factors listed below: o The asymmetry of the geometric arrangement o Asymmetry of the voltage waveform o Mobility of the dominant species resulting from the ionization Overall, the paraelectric momentum coupli ng mechanism is due to the existence of a dominant force in one direction for a bulk of the time-period and for a lesser time in the other direction. The cumulative effect over time is primarily responsible for the unidirectional coupling. Efficiency of the plasma actuator was analyzed using asymmetric sawtooth waveforms which revealed that the pos itive going voltage is the most efficient force generating part of the half-cycle. The ge ometric effects analyzed revealed that the presence of a larger lower el ectrode will increase the discharge formation region, but will lead to diminishing returns sue to a weakening discharge. Parametric studies such as sensitivity to voltage and frequency were conducted. The effect of frequency was quite modest as compared to the sensitivity to the applied voltage which showed a power-law dependenc e on the voltage for the measure force in the domain. This power-law relationship showed substantial similarity to the experimental measurements and the phenomenological model. It is worth noting here, that for the helium gas the observed species number density (O(109)/cm3) is a couple of orders of magnitude smaller than that observed for air (O(1011)/cm3). Also, the force measurement shown is for a unit centimeter length of the actuator in the spanwise direction and explains the very small realized in the high-fidelity

PAGE 148

133 approach. The quadratic dependence of the for ce on the voltage is significant for realizing large momentum effects using relatively small increases in voltage.

PAGE 149

134 CHAPTER 5 SUMMARY AND FUTURE WORK 5.1 Summary and Conclusions A glow discharge-induced fluid dynami cs model has been proposed where the plasma and the neutral flow are modeled as a two-fluid system with the plasma dynamics and the fluid dynamics solved separately. As a first step we, present an argument to decouple the plasma dynamics and the fluid dynamics using the disparity in the timescales of the associated physics. The pres ent regime of interest, which is the low Reynolds number flow, operates at a much slow er time-scale as compared to the faster discharge breakdown and evolution processes, generated by the high frequency applied voltage. The success of this multi-scale handling of the plasma-fluid flow system will now be determined by an efficient treatment of the plasma model. To integrate the plasma dynamics and the fluid dynamics, a body force model was developed to represent the effects of the discharge. The instantaneous body force which is the local Lorentz force acting on the charged species operates on the plasma scales and need to be averaged over the smallest fluid dynamic time scale. Two modeling approaches have been developed. (i) A phenomenological modeling approach ba sed on linearized force distribution and experimentally observed discharge structure was developed. This model though requires the identification of empirical constants, showed reasonable qualitative agreement with the experimentally observed flow features This model also captured parametric

PAGE 150

135 dependencies such as the dependence of the force on the voltage comparable to that obtained by the first-principle based a pproaches and experimental studies. This approach obviates the need for computing accurate descriptions of the plasma species and hence, th e electric field distribution. This model, though required experimental characterization to achieve closure, was able to capture some of the reasonable flow features such as, the wall je t formation and other observed experimental outcomes. (ii) In order to evolve the capability to the next levela high fide lity, first principles based –approach using a plasma-fluid (or plasma-hydrodynamic) continuum formulation to model the discharge dynamics. First, a one-dimensional two-species plasma model with reduced chemistry served as a sample test problem for validating the input modeling data and in identifying the various time-s cales and modeling limitations. The model was in reasonable quantitative agreement with the experimental results. To extend to multiple dimensions, an operator splitting algorithm capable of handling multiple time-scales was employed. The resulting multi-scale algorithm enabled the integration of the overall system at the slowest species convection time-scale which is 100 times the dielectric relaxation time-scale. The asymmetric plasma discharge in helium has at 300 torr and 300 K modeled using the hydrodynamic-plasma model showed characteristic temporal and spatial structure which was consistent with the experimentally observed results of Enloe et al .[21] for atmospheric air. Furthermore, the following conclusions can be drawn. The two-species plasma simulations revealed that the generation of the uni-directional moment um coupling is primarily caused by the combination of factors listed below:

PAGE 151

136 The difference in the mobility or convecti on of the two plasma species consideredthe faster electron species and the slower ion species Difference in the behavior of the abso rbing conductor surface (electrode) and the self-limiting charge build-up of the electrode. Lastly, the geometric asymmetry due the downstream submerged electrode is closely coupled to the two above mentioned observations. Overall, the paraelectric momentum coup ling mechanism is due to the cumulative effect over time of the force field in the domain, as seen from our computations. The asymmetric sawtooth waveforms revealed that the positive going voltage is the most efficient force generating part of the half-cyc le. The geometric effects analyzed revealed that the presence of a larger lower electrode will increase the discharge formation region, but will lead to diminishing returns sue to a weakening discharge. Parametric studies such as sensitivity to voltage and frequency were conducted. The applied voltage revealed a power-law dependence on the voltage for the measured force in the domain. This power-law relationship showed bears resemblance to the experimental measurements and the phenomenological model. Finally, the neutral fluid on which the plasma acts is assumed incompressible because of the low speed regime of interest in this dissertation and the experimental observations indicating very little heati ng and hence minimal to no buoyancy related effects. The effects considered in this dissert ation are purely mechanical in nature as we assumed the discharge to be in local th ermal equilibrium by assuming an efficient collision process leading to more moment um exchange and less thermal heating. However, future efforts might attempt to investigate the potential of these actuators for

PAGE 152

137 converting the input electic power to gene rate gas heating effects. Based on the phenomenological model for air, the current actuator configuration presented in this dissertation displayed sufficient potential for flow control applications for Re in the range of O(1000). The plasma induced body force effect is similar in nature to the pressure gradient and stronger voltages can enhance the potential performance range of these actuators to O(10000-100000) taking advant age of the power-law dependence. 5.2 Future Work Based on the accomplishments to date, further efforts are proposed in the following to further advance the predictive capabilities, and to extend the approach to engineering applications including drag reduction, thermal management, and flow control. The multi-scale nature of the solution can be further aided by using local grid refinements techniques near the regions of sharp ionization wave propagation. The time and length scales in these regions are of O(~106), and consequently, can be better resolved with grid refinement. As evidenced in Chapter 4, the discharge is heavily characterized by breakdown phases followed by the convection phase where the charge species accumulate on th e material surface. Th ese processes act in quick succession resulting in the pres ence of ionization wave-type phenomena which is the reason for the stiffness in the present case of plasma modeling. Performance, stability and accuracy cons iderations of the sequential algorithm employed in the present study should be furth er assessed with different physical and operating parameters.

PAGE 153

138 The study of fluid–plasma interactions is common in the area of flow-control applications with a fluid such as air. In such cases the plasma chemistry needs to be modeled in a mixture of nitrogen and oxygen. For example, the role of chemistry, i.e., how the presence of both heavy negativ e and positive ions can affect resulting flow structure needs to be investigated. In the present effort, only Helium gas was considered due to the wealth of modeling l iterature available. However, as we have seen, the discharge evolution is critically dependent on the nature and mobility of the individual species concerned. Hence, in order to get a realistic representation of the performance of the actuator in air, we will need to expand the chemistry modeling.

PAGE 154

139 LIST OF REFERENCES 1. S. Hippler, M. Pfau, K.H. Schmidt, and Schoenbach, Low Temperature Plasma Physics: Fundamental Aspects and Applications, Wiley-VCH, Berlin (2001). 2. W.Siemens, “Plasma discharges,” Poggendroffs Ann. Phys. Chem. 102, pp 66 (1857). 3. U.Kogelschatz, B.Eliasson and W.Egli, “Fro m ozone generators to flat television screens: history and future potential of di electric-barrier discharges,” Pure. Appl. Chem., 71(10), pp 1819-1828 (1999). 4. J.R.Roth, D.M.Sherman and S.P. Wilkins on, “Boundary layer flow control with a one atmosphere uniform glow discharge,” AIAA paper 98-0328, 36th AIAA Aeropsace Sciences Meeting, Reno NV, January 12-15 (1998). 5. S.Kanazawa, M. Kogoma, T.Moriwaki and S.Okazaki, “Stable glow plasma at atmospheric pressure,” Phys. D. Appl. Phys. 21, pp 838 (1988). 6. F.Massines, R.Ben Gadri, A.Rabehi, Ph.Decomps, P.Segur and Ch.Mayoux, “Experimental and theoretical study of a gl ow discharge at atmospheric pressure controlled by dielectric barrier,” J. Appl. Phys, 83(6), pp 2950-2957 (1998). 7. J.R.Roth, “Electrohydrodynamically induced airflow in a one atmosphere uniform glow discharge surface plasma”, Paper 6P-67, Proc. 25th IEEE International Conference on Plasma Science, Raleigh, NC June 1-4, pp 291 (1998). 8. B.L.Smith and A.Glezer, “The formation and evolutionof synthetic jets,” Phys. of Fluids, 10(9), pp 2281–2297 (1998). 9. W.Shyy, B.Jayaraman and A.Anderson, “Mode ling of glow-discharge induced flow dynamics,” J. Appl. Phys., 92(11), pp 6434-6443 (2002). 10. J.S.Shang, “Recent research in magneto-aerodynamics,” Prog. Aero. Sci., 37, pp 120, (2001). 11. R.Rivir, A.White, C.Carter, B.Ganguly, J. Jacob, A.Forelines and J.Crafton, “AC and pulsed plasma flow control,” AIAA paper 2004-0847, 42nd Aerospace Sciences Meet & Exhibit, Reno, NV, Jan 5-8 (2004).

PAGE 155

140 12. S.Leonov, V.Bityurin and D.Yarantsev, “The effect of plasma induced separation,” AIAA paper 2003-3853, 34th AIAA Plasmadynamics and Lasers Conference, Orlando, FL, Jun 23-26 (2003). 13. Y.B.Golubovskii, V.A Maiorov, J.Behnke and J.F.Behnke, “Modeling of homogeneous barrier discharge in helium at atmospheric pressure,” J.Phys. D: Appl. Phys. 36, pp 39-49 (2003). 14. D.V.Gaitonde, M.R.Visbal and S.Roy, “Control of flow past a wing section with Plasma-based Body forces,” AIAA paper 2005-5302, 36th Plasma Dynamics and Lasers Conference, Toronto, 6-9 June (2005). 15. S.Roy and D.V.Gaitonde, “Modeling surface discharge effects of atmospheric RF on gas flow control,” AIAA paper 2005-0160, 43rd AIAA Aeropsace Sciences Meeting, Reno NV, January 10-13 (2005). 16. S.Roy and D.V.Gaitonde, “Multidimensional collisional dielectric barrier discharge for flow separation control at atmosp heric pressures,” AIAA paper 2005-4631, 35th AIAA Fluid dynamics conference, To ronto, Canada, 6-9 June (2005). 17. K.P.Singh and S.Roy, “Simulation of an asymmetric single dielectric barrier plasma actuator,” J. App. Phys. 98, pp 083303-1 (2005). 18. H.Kumar and S.Roy, “Multidimenisonal hydrodynamic plasma-wall model for collisional plasma discharges with and without magnetic-field effects,” Phys. Of Plasmas, 12, pp 093508 (2005). 19. S.Thakur, J.Wright and W.Shyy, “STREA M: A computational fluid dynamics and heat transfer Navier-Stokes solver: theory and applications,” Application manual, Gainesville, FL: Streamline Numerics, Inc. (2002). 20. R.Ben Gadri, “One atmospheric glow di scharge structure revealed by computer modeling,” IEEE Trans. Plasma Sci., 27(1) (1999). 21. C.L.Enloe, T.E.Mclaughlin, R.D.Van Dyken and K.D. Kachner, “Mechanisms and responses of a single dielectric barrier plasma actuator: plasma morphology,” AIAA J. 42(3) (2004). 22. C.L.Enloe, T.E.Mclaughlin, R.D.Van Dyken, K.D. Kachner, E.J.Jumper, T.C.Corke, M.Post and O.Haddad “Mechanisms and responses of a single dielectric barrier plasma actuator: geometric effects,” AIAA J. 42(3) (2004). 23. C.L.Enloe, T.E.Mclaughlin, R.D.Van Dyken and J.C.Fischer, “Plasma structure in the aerodynamic plasma actuator,” AIAA paper 2004-0844, 42nd Aerospace Sciences Meet & Exhibit, Reno, NV, Jan 5-8, 2004. 24. T.C.Corke and E.Matlis, “Phased plasma a rrays for unstaeady flow control,” AIAA paper 2000-2323, Fluids 2000, Denver, CO, June 19-22 (2000).

PAGE 156

141 25. T.C.Corke, E.J.Jumper, M.L.Post, D.O rlov and T.E.Mclaughlin, “Application of weakly-ionized plasmas as wing flow-control devices,” AIAA paper 2002-0350, 41st Aerospace Sciences Meeting & E xhibit, Reno, NV, Jan 6-9 (2003). 26. R.Van Dyken, T.E.Mclaughlin and C.L.Enloe, “Parametric investigations of a single dielectric barrier plasma actuator,” AIAA paper 2004-846, 42nd Aerospace Sciences Meet & Exhibit, Reno, NV, Jan 5-8 (2004). 27. M.L.Post and T.C.Corke, “Seperation contro l on high angle of attack airfoil using plasma actuators,” AIAA paper 2003-1024, 41st Aerospace Sciences Meet & Exhibit, Reno, NV, Jan 6-9 (2003). 28. B.Jayaraman, S.Thakur and W.Shyy, “Mode ling of dielectric barrier discharge and resulting fluid dynamics,” AIAA paper 2006-0686, 44th Aerospace Sciences Meeting & Exhibit, Reno, NV, January 9-12 (2006). 29. B.Jayaraman, W.Shyy and S.Thakur, “M odeling of fluid dynamics and heat transfer induced by dielectric barrier plasma actuator,” J. Heat Trans., Invited paper in the Special Edition in Honour of Dr S.V. Patankar, submitted for review. 30. R.B.Campbell, R.Veerasingham and R.T.Mcgrath, “A two-dimensional multispecies fluid model of the plasma in an AC plasma display panel,” IEEE Trans. on Plasma Sci., 23(4), pp 698-708 (1995). 31. J.P.Boeuf and L.C.Pitchford, “Two-dimensi onal model of a capacitively coupled rf discharge and comparisons with experiments in the gaseous electronics conference reference reactor,” Phys. Rev. E, 51(2), pp 1376-1390 (1995). 32. J.R.Roth, Industrial Plasma Engineering: Vol. I-Principles, IOP, Bristol (1995). 33. J.R.Roth, Industrial Plasma Engineering: Vol. IIApplications to Non-thermal Plasma Processing, IOP, Bristol (2001). 34. D.Lee, J.M.Park, S.H.Hong and Y.Kim, “N umerical simulation on mode transition of atmospheric dielectric barrier discha rge in helium-oxygen mixture,” IEEE Trans. on Plasma Sci., 33(2), pp 949-957 (2005). 35. C.Baird, C.L.Enloe, T.E.Mclaughlin and J.W.Baughn, “Acoustic testing of dielectric barrier discharge (DBD) pl asma actuator,” AIAA paper 2005-0565, 43rd Aerospace Sciences meeting and Exhibit, Reno, NV, Jan 10-13 (2005). 36. C.L.Enloe, T.E.Mclaughlin, G.I.Font and J.W.Baughn, “Parameterization of temporal structure in the single dielec tric barrier aerodynami c plasma actuator,” AIAA paper 2005-0564, 43rd AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan 10-13 (2005).

PAGE 157

142 37. T.N.Jukes, K.S.Choi, G.A.Johnson and S.J. Scott, “Characterization of surface plasma-induced wall flows through velocity and temperature measurements”, AIAA J., 44(4), pp 764-771 (2005). 38. T.N.Jukes, K.S.Choi, G.A.Johnson and S.J.Scott, “Turbulent boundary-layer control for drag reduction using surface plasma,” AIAA paper 2004-2216, 2nd AIAA Flow Control Conference, Portland, Oregon, Jun 28Jul 1 (2004). 39. J.R.Roth, H.Sin, R.C.M.Madhan and S.P.Wilkinson, “Flow re-attachment and acceleration by paraelectric and perist altic electrohydrodynamic effects,” AIAA paper 2003-0351, 41st Aerospace Sciences Meet & Exhibit, Reno, NV, Jan 6-9 (2003). 40. J.List, A.R.Byerly, T.E.Mclaughlin and R. D.Vandyken, “Using plasma actuator to control laminar separation on a linear cascade turbine blade,” AIAA paper 20031026, 41st Aerospace Sciences Meet & Exhi bit, Reno, NV, Jan 6-9 (2003). 41. A.Santhanakrishnan, N.J.Pern, K.Ramkumar A.Simpson and J.D.Jacob, “Enabling flow control technology for low speed UAV’s,” AIAA paper 2005-6960, Infotech at aerospace, Arlington, Virg inia, Sept. 26-29 (2005). 42. B.Goskel, D.Greenblatt, I Rechenberg, C.N.Nayeri and C.O.Paschereit, “Steady and unsteady plasma wall jets for separation and circulation control,” AIAA paper 2006-3686, 3rd AIAA Flow Control Conference, San Francisco, CA, Jun 5-8 (2006). 43. T.C.Corke, B.Mertz and M.P.Patel, “Plasma flow control optimized airfoil,” AIAA paper 2006-1208, 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan 9-12 (2006). 44. S.Chan, X.Zhang and S.Gabriel, “The attenuation of cavity tones using plasma actuators,” AIAA paper 2005-2802, 11th AI AA/CEAS Aeroacoustics Conference (26th AIAA Aeroacoustics Conference) 23 25 May 2005, Monterey, California (2005). 45. R.Rivir, A.White, C.Carter, B.Ganguly, A.Forelines and J.Crafton, ”Turbine flow control, plasma flows,” AIAA paper 2003-6055, 41st Aerospace Sciences Meet & Exhibit, Reno, NV, Jan 6-9 (2003). 46. A.Soldati and S.Banerjee, “Turbulence modification by large-scale organized electrohydrodynamic flows,” Phys. of Fluids, 10(7), pp 1742-1756 (1998). 47. A.Soldati and C.Marchioli, “Prospects fo r modulation of turbulent boundary layer by EHD flows,” Turbulence Structure and Modulation, Springer-Verlag, Wien, pp 119-160 (2001).

PAGE 158

143 48. S.P.Wilkinson, “Investigation of an osci llating surface plasma for turbulent drag reduction,” AIAA paper 2003-2199, 41st Aerospace Sciences Meet & Exhibit, Reno, NV, Jan 6-9 (2003). 49. Y.Du, V.Symeonidis and G.E.Karniadak is, “Drag reduction in wall-bounded turbulence via a transverse traveling wave ,” J. Fluid Mech., 457, pp 1-34 (2002). 50. M.R.Dhanak and C.Si, “On reduction of turbulent wall friction through spanwise wall oscillations,” J. Fluid Mech., 383, pp 175-195 (1999). 51. K.D.Hall, E.J.Jumper, T.C.Corke and T.E.Mclaughlin, “Potential flow model of a plasma actuator as a lift enhancement device,” AIAA paper 2005-0783, 43rd Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan 10-13 (2005). 52. D.M.Orlov and T.C.Corke, Numerical simulation of aerodynamic plasma actuators,” AIAA paper 2005-1083, 43rd Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan 10-13 (2005). 53. J. Tepper and M.Lindmayer, “Investig ations on two different kinds of homogeneous barrier discharges at atmos pheric pressure,” Proc. of HAKONE VII, 7th Int. Symp. On High Pressure Low Temperature PlasmaChemistry, Greifswald, vol. 1, Sep 10-13, pp 38 (2000). 54. B. Eliasson and U.Kogelschatz, “Modeling and applications of silent discharge plasmas,” IEEE Trans. Plasma Sci. 19(2), pp 309-323 (1991). 55. E.Gogolides and H.H.Sawin, “Continuum modeling of radio-frequency glow discharges. I. Theory and results for electropositive and electronegative gases,” J. Appl. Phys. 72(9), pp 3971-3987 (1992). 56. L.Mangolini, K.Orlov, U.Kortshagen, J. Heberlein and U.Kogelschatz, “Radial structure of a low-frequency atmospheric -pressure glow discharge in helium,” Applied physics letters, 80(10), pp 1722-1724 (2002). 57. Y.B.Golubovskii, V.A Maiorov, J.Behnke and J.F.Behnke, “Study of homogeneous glow-like discharge in nitrogen at atmosphe ric pressure,” J.Phys. D: Appl. Phys. 37, pp 1346-1356 (2004). 58. P.Segur and F.Massines, “The role of numerical modeling to understand the behavior and to predict the existence of an atmospheric pressure glow discharge controlled by a dielectric ba rrier,” Proc. Of the XIIIth Int. Conf. of Gas Discharges and their applications, Glasgow, Sep. 3-8, pp 15-24 (2000). 59. X.Xu and M.J.Kushner, “Multiple microd ischarges dynamics in dielectric barrier charges,” J. Appl. Phys. 84(8), pp 4153-4160 (1998).

PAGE 159

144 60. Y.B.Golubovskii, V.A Maiorov, J.Behnke and J.F.Behnke, “Influence of interaction between charged particles and dielectric surface over a homogeneous barrier discharge in nitrogen,” J.Phys. D: Appl. Phys. 35, pp 751-761 (2002). 61. J.Shin and L.L.Raja, “Dynamics of puls e phenomena in helium dielectric-barrier atmospheric-pressure glow discharges,” J. Appl. Phys. 94(12), pp 7408-7415 (2003). 62. L.Mangolini, C.Anderson, J.Heberlein and U.Kortshagen, “Effects of current limitation through the dielectric in atmospheric pressure glows in helium,” J.Phys D: Appl. Phys., 37, pp 1021-1030 (2004). 63. F.Tochikubo, T.Chiba and T.Watanabe, “S tructure of low-frequency helium glow discharge at atmospheric pressure between parallel plate dielectric electrodes,” Jpn. J. Appl. Phys. 38, pp 5244-5250 (1999). 64. T.R.Govindan and M.Meyappan, “One-dimensional modeling studies of the gaseous electronic conference RF refere nce cell,” J. Research, NIST, 100(4), pp 463-472 (1995). 65. X.Yuan and L.L.Raja, “Computational study of capacitively coupled high-pressure glow discharges in helium,” IEEE Tran s. Plasma Sci., 31(4), pp 495-503 (2003). 66. J.P.Boeuf, “Numerical modeling of rf gl ow discharges,” Phys. Rev. A, 36(6), pp 2782-2792 (1987). 67. C.K.Birdsall, “Particle-in-cell chargedparticle simulations, plus Monte Carlo collisions with neutral atoms,” IEEE Tr ans. Plasma Sci., 19(2), pp 65-85 (1991). 68. T.E.Nitschke and D.B.Graves, “A comparison of particle in cell and fluid model simulations of low-pressure radio freque ncy discharges,” J. Appl. Phys. 76(10), pp 5646-5660 (1994). 69. S.K.Dhali and P.F.Williams, “Two-dimensional studies of streamers in gases,” J. Appl. Phys. 62, pp 4696 (1987). 70. D.P Lymberopoulos and D.J.Economou, “T wo-dimensional self-consistent radio frequency plasma simulations relevant to the gaseous electronics conference RF refernce cell,” J. Research of the NIST, 100(4) (1995). 71. J.Meunier, Ph.Belenguer and J.P.Boeuf, “Num erical model of an ac plasma display panel cell neon-xenon mixtures,” J. Appl. Phys 78(2), pp 731-745 (1995). 72. R. Veerasingham, R.B.Campbell and R .T.Mcgrath, “One-dimensional single and multipulse simulations of the ON/OFF voltages and the bistable margin for He, Xe filled plasma display panels,” IEEE Trans. Plasma Sci., 24(6), pp 1399-1410 (1996).

PAGE 160

145 73. W.M.Hilbun and B.J.Case, “Preliminary de velopment of a com putational model of a dielectric barrier discharge,” AIAA paper 2005-1176, 43rd Aerospace Sciences Meeting and Exhibit, Reno, NV, Jan 10-13 (2005). 74. B.Jayaraman and W.Shyy, “Flow control a nd thermal management using dielectric glow discharge concepts,” AIAA paper 2003-3712, 33rd AIAA Fluid Dynamics Conference and Exhibit, Orlando, June 23-26 (2003). 75. K.P.Singh, S.Roy and D.V.Gaitonde, “Mode ling of dielectric barrier discharge plasma actuator with atmospheric air chemistry,” AIAA paper 2006-3381, 37th AIAA Plasmadynamics and Lasers conference, SanFrancisco, CA, 5-8 June (2006). 76. M.R.Visbal, D.V.Gaitonde and S.Roy, “Con trol of transitional and turbulen flows using plasma-based actuators”, AIAA paper 2006-3230, 36th AIAA Fluid Dynamics Conference and Exhibit, San Francisco, CA, 5-8 June, 2006. 77. M.E.Riley, K.E.Greenberg, G.A.Hebner and P.Drallos, “Theoretical and experimental study of low-temperature, capacitiviley coupled, radio-frequency helium plasmas,” J. Appl. Phys, 75(6), pp 2789-2798 (1994). 78. M.J.Kushner, “Mechanisms for power de position in Ar/SiH4 capacitively coupled RF discharges,” IEEE Trans. Pl asma Sci., 14(4), pp. 188 (1985). 79. R.W.Boswell and I.J.Morey, “Self-consistent simulation of a parallel-plate rf discharge,” Appl. Phys. Lett. 52(1), pp 21-23 (1988). 80. M.Surendra and D.B.Graves, “Particle simulations of radio-frequency glow discharges,” IEEE Trans. Plas ma Sci., 19(2), pp 144-157 (1991). 81. W.Nicholas G.Hitchon, T.J.Sommerer and J .E Lawler, “A self-consistent kinetic plasma model with rapid convergence,” IEEE Trans. Plasma Sci., 19(2), pp 113121 (1991). 82. J.V.DiCarlo and M.J.Kushner, “Solving the spatially dependent Boltzmann’s equation for the electron velo city distribution using flux corrected transport,” J. Appl. Phys. 66(12), pp 5763-5774 (1989). 83. W.N.G.Hitchon, D.J.Koch and J.B.Adams, “An Efficient scheme for convection dominated transport,” J. Co mput. Phys. 83, 79 (1989). 84. Ph.Belenguer and J.P.Boeuf, “Transition between different regimes of rf glow discharges,” Phys. Rev. A, 41, pp 4447 (1990). 85. A.Fiala, L.C.Pitchford and J.P.Boeuf, “Two-dimensional, hybrid model of lowpressure glow discharges,” Phys. Rev. E, 49(6), pp 5607 (1994).

PAGE 161

146 86. P.Colella, M.R.Dorr and D.D.Wake, “A c onservative finite difference method for the numerical solution of plasma fluid equations,” J. Comp. Phys. 149, pp 168-193 (1999). 87. M.Surendra, “Radio frequency discharge benchmark model comparison,” Plasma Sources Sci. Technol. 4, pp 56-73 (1995). 88. P.L.G. Ventzek, T.J.Sommerer, R.J.Hoekstra and M.J.Kushner, “Two-dimensional hybrid model of inductively coupled plasma sourced for etching,” Appl. Phys Lett. 62, pp 3247-3249 (1993). 89. T.J.Sommerer and M.J.Kushner, “Numeri cal investigation of the kinetics and chemistry of rf glow discharg e plasmas sustained in He, N2, O2, He/ N2/ O2, He/CF4/O2, and SiH4/NH3 using a Monte Carlo-fluid hybr id model,” J. Appl. Phys. 71(4), pp 1654-1673 (1992). 90. E.P.Hammond, K.Mahesh and P.Moin, “A Numerical method to simulate radiofrequency plasma discharges,” J.Comp. Phys., 176, pp 402-429 (2002). 91. Shyy, W. “A Study of finite difference appr oximations to steady-state, convection dominated flow problems,” J. Comp. Phys. 57, pp 415-438 (1985). 92. D.A.Knoll and D.E.Keyes, “Jacobian-fr ee Newton Krylov Methods: a Survey of Approaches and Applications,” J. Comp. Phys. 193, pp 357-397 (2004). 93. E. Turkel, “Review of preconditioning me thods for fluid dynamics,” NASA CR189712, ICASE Report no 92-47 (1992). 94. C.L.Merkle, J.Y.Sullivan, P.E.O.Buelow and S.Venkateswaran, “Computations of flows with arbitrary equations of state,” AIAA J. 36(4), pp. 515, 1998. 95. D.A.Schwer, P.Lu, W.H.Green Jr. and V.Se miao, “A Consistent-splitting approach to computing stiff steady-state reacting flows with adaptive chemistry,” Comb. Theory Modeling, 7, pp 383-399 (2003). 96. G.J.M. Hagelaar and G.M.W. Kroesen, “Speeding up fluid models for gas discharges by implicit treatment of the electron energy source term,” J. Comp. Phys. 159(1), pp 1-12 (2000). 97. T.R.A.Bussing and E.M.Murman, “Finitevolume method for the calculation of compressible chemically reacting flows,” AIAA J. 26(9), pp 1070-1078. 98. S.Yoon and A.Jameson. “Lower-upper symmetric Gauss-Seidel method for the Euler and Navier-Stokes Equations ,” AIAA J. 26(9), pp1025 (1988). 99. J.G.Verwer, “Explicit Runge-Kutta methods for parabolic partial differential equations,” Appl. Numer. Math. 22, pp 359-379 (1996).

PAGE 162

147 100. M.Valorani and D.A. Goussis, “Explicit time-scale splitting algorithms for stiff problems: auto-ignition of gaseous mixt ures behind a steady shock,” J. Comp. Phys. 169 pp 44-79 (2001). 101. E.S.Oran and J.P.Boris, Numerical Simulation of Reactive Flow, Cambridge university press, 2nd Edition (2001). 102. H.N.Najm, P.S.Wyckoff and O.M. Knio, “A semi-implicit numerical scheme for reacting flow,” J. Comp. Phys. 143, pp 381-402 (1998). 103. P.N.Brown, G.D.Byrne, and A.C.Hindmar sh, “VODE: A variable coefficient ODE solver,” SIAM J. Sci. Stat. Comput. 10, pp 1038 (1989). 104. R.Tyson, L.G.Stern and R.J.Leveque, “Fractional steps method applied to chemotaxis model,” J. Math. Biol. 41, pp 455-475 (2000). 105. S.D.Cohen, A.C.Hindmarsh, “CVODE: a stiff/nonstiff ODE solver in C,” Comp.Phys. 10(2), pp 138-143 (1996). 106. G.Strang, “On the construction and compa rison of difference schemes,” SIAM J. Num. Anal. 21, pp 994-1011 (1969). 107. B.Sportissee, “An analysis of operator sp litting techniques in the stiff case,” J. Comp. Phys. 161, pp 140-168 (2000). 108. O.M.Knio, H.N.Najm and P.S.Wyckoff, “A semi-implicit numerical scheme for reacting flow,” J.Comp. Phys. 154, pp 428-467 (1999). 109. G.J.M.Hagelaar, F.J. de Hoog and G.M.W. Kroesen, “Boundary conditions in fluid models of gas discharges,” Phs. Rev. E. 62(1), pp 1452-1454 (2000). 110. C.Enloe, T.Mclaughlin, R.Van Dyken and J.Fischer, “Plasma structure in the aerodynamic plasma actuator,” AIAA paper 2004-0846, 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, Jan. 5-8 (2004). 111. S.Thakur and J.Wright, “A multi-block operator splitting algorithm for unsteady flows at all speeds in complex geometries,” Int. J. Num. Meth. Fluids, 46, pp 383413 (2004). 112. W.Shyy, S.Thakur, H.Ouyang, J.Liu and E.Blosch, Computatitional Techniques for Complex Transport Phenomena, Cambridge University Press, Cambridge, United Kingdom (1997). 113. W.Shyy, Computational Modeling for Flui d Flow and Interfacial Transport, Elsevier, Amsterdam, Netherlands (1994).

PAGE 163

148 114. M.Surendra, D.B.Graves and G.M.Jellum, “Self-consistent model of a directcurrent glow discharge: Treatment of fast electrons,” Phys. Rev. A, 41(2), pp 11121126 (1990). 115. D.V.Gaitonde, M.R.Visbal and S.Roy, “A coupled approach for plasma-based flow control simulations of wing sections,” AIAA paper 2006-1205, 44th AIAA Aerospace Engineering Sciences Meeting and Exhibit, Reno, NV, 9-12 Jan (2006). 116. P.Colella, M.R.Dorr and D.D.Wake, “Numerical solution of plasma fluid equations using locally refined grids,” J. Comp. Phys. 152, pp 550-583 (1999).

PAGE 164

149 BIOGRAPHICAL SKETCH Balaji Jayaraman was born in the city of Madras in southern India in 1980. After completing his preliminary schooling, he was fortunate to pursue his Bachelor of Technology degree in naval architecture and ocea n engineering at the Indian Institute of Technology, Madras. On completion of his B.Tech, he began his graduate school to pursue his Ph.D. at UF under Dr. Wei Shyy in the fall of 2001.


Permanent Link: http://ufdc.ufl.edu/UFE0015702/00001

Material Information

Title: Computational Modeling of Glow Discharge-Induced Fluid Dynamics
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0015702:00001

Permanent Link: http://ufdc.ufl.edu/UFE0015702/00001

Material Information

Title: Computational Modeling of Glow Discharge-Induced Fluid Dynamics
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0015702:00001


This item has the following downloads:


Full Text











COMPUTATIONAL MODELING OF GLOW DISCHARGE-INDUCED FLUID
DYNAMICS
















By

BALAJI JAYARAMAN


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2006

































Copyright 2006

by

Balaji Jayaraman




























This document is dedicated to my parents and all my teachers from whom I have learnt
all my life.















ACKNOWLEDGMENTS

I would like to thank my advisor, Professor Wei Shyy, for his advice, inspiration,

patience and mentoring. He constantly provided insightful suggestions when I needed

help and motivation.

I would also like to express gratitude to my committee members: Dr. Siddharth

Thakur, Professors Corin Segal, Toshi Nishida, Lou Cattafesta and William Lear. A

special mention needs to be made of Dr.Thakur for his suggestions regarding the thesis

content. I have also gained significantly from the various technical discussions which I

have had with him. Finally, I am forever indebted to my parents for their understanding

and support when required. This acknowledgement will not be complete without the

mention of the role played by my research group members for their support and

friendship.
















TABLE OF CONTENTS



A C K N O W L E D G M E N T S ................................................................................................. iv

LIST OF TABLES ............... ......... ... .. ......... ............................ vii

LIST OF FIGU RE S ......... ........................... ............ ........... ix

N O M E N C L A T U R E ......................................................................................................... x ii

ABSTRACT ................................. ....... .............. xiv

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

1.1 Overview of Plasma Generation and Dielectric Barrier Discharges (DBD's) .......1
1.2 Plasma Actuator Operating M echanism ...................................... ............... 3
1.3 N um erical C hallenges.......... ..... .................................................. ................ .6
1.4 O verview of Present A pproach................................................................... ......7
1.5 O utline of the D issertation ........................................................... ............... 9

2 L ITE R A TU R E R E V IEW ....................................................................... ..................12

2.1 Review of Experimental Discharge Plasma Studies............................................12
2.1.1 E arly D ischarge Studies ....................... ...................... ............... ... 13
2.1.2 Mechanism of Plasma-Induced Flow Applications................................ 15
2.1.3 Aerodynamic and Control Strategies Using Glow Discharges ..................20
2.2 Review of Discharge Modeling Studies ............... .............................25
2.2.1 Phenomenological Modeling of the Plasma Actuator..............................26
2.2.2 First Principles-Based Plasma Modeling Approaches ............................27
2.2.2.1 Review of dielectric barrier discharge (DBD) modeling studies .....27
2.2.2.2 Review of plasma actuator modeling studies...............................33

3 NUMERICAL MODELING OF GLOW DISCHARGE PLASMA .............................36

3.1 P lasm a M odeling H ierarchy ...................................................... .....................36
3.1.1 Kinetic Modeling Approach ............. ............................................. 38
3.1.2 H ybrid M odeling Approach ............................................ ............... 42
3.1.3 Fluid M odeling A approach ........................................ ....... ............... 44
3.2 Numerical Modeling of Glow Discharge Plasmas ............................................47
3.2.1 Plasma Discharge Governing Equations ..................................................48









3.2.2 Numerical Algorithm ....................... ..... ....... ........................ 59
3.2.3 1-D Parallel Plate Discharge Modeling in Helium................................74
3.3 Integration of Plasma Dynamics with Fluid Dynamics.......................................78

4 MODELING FLUID FLOW WITH PLASMA EFFECTS ............................................81

4 .1 B ody F force F orm ulation ............................................................ .....................82
4.2 Navier-Stokes Equations with Body Force.........................................................84
4.3 Modeling a Plasma Actuator Using a Phenomenological Force Model ..............85
4.3.1 Plasma Generation in an Asymmetric Electrode Configuration ................85
4.3.2 Linearized Electric B ody Force.............................................................. 88
4.3.3 Problem D description ........................................................ ............... 91
4.3.4 Effect on Flow Structure ........................................ ........................ 93
4.3.5 Effect on H eat Transfer ........................................ ......................... 97
4.3.6. Thrust E stim ation .............................................. ............. ............... 98
4.3.7 Effect of Plasma on Aerodynamics ....................................................102
4.4 Self-Consistent Fluid-Model-Based Plasma Dynamics.................................... 108
4.4.1 Discharge Actuator Setup............................. .............. .............. 109
4.4.1.1 Boundary condition .......................................111
4.4.1.2 Plasma species initial condition ...............................112
4.4.2 Plasm a Structure .......................................... .......... .. ........... .... 113
4.4.3 G eom etric Effects .......................................................... ............. 125
4.4.3.1 Impact of electrode spacing............................. 125
4.4.3.2 Im pact of low er electrode size .............................................127
4.4.4 Impact of Applied Voltage and Frequency ...........................................128
4.4.4.1 Frequency sensitivity.................................. ........................ 129
4.4.4.2 A applied voltage sensitivity............................................................130
4.5. Sum m ary .............................................................................. 131

5 SUMMARY AND FUTURE WORK .............................................................. 134

5.1 Sum m ary and Conclusions ...........................................................................134
5.2 Future W ork .................................... ............................... ........ 137

LIST OF REFEREN CES ........................................................... .. ............... 139

BIOGRAPHICAL SKETCH ............................................................. ............... 149
















LIST OF TABLES


Table page

1 Representative EHD operating mechanisms........................ ....................... 18

2 Summary of experimental plasma actuator studies ................................ ............... 22

3 Summary of turbulent flow control studies ...................... ...............25

4 Summary of low frequency plasma modeling studies..................................................31

5 Summary of high frequency plasma modeling studies...............................................33

6 Summary of plasma-induced fluid flow modeling studies ..........................................35

7 Kinetic model studies of discharge plasmas ....................................... ............... 40

8 Summary of hybrid modeling studies................................ ........................ ......... 43

9 Plasma chemistry considered in a helium discharge.....................................................46

10 Representative time scales of various processes in a 1-D helium discharge model
operating at 10KHz, 1.5KV rms, calculated based on values of parameters. .........60

11 Summary of governing equations used in different studies............... ...................61

12 Different types of boundary conditions for the ions and electrons used in modeling
stu d ie s ...................................... .................................................... 7 6

13 Change in drag coefficient and drag force with Reynolds number ..........................104

14 Change in lift coefficient and lift force with Reynolds number ..............................104

15 Effect of electrode position on aerodynamic and heat transfer properties................07

16 Summary of property models employed for the He discharge simulation. ...............11

17 Summary of boundary conditions for the different variables...................................113

18 Domain averaged force over the voltage cycle for different waveforms.................. 125

19 Domain averaged force over the voltage cycle for different configurations ............127









20 Force dependence on the lower electrode size................................. ............... 128

21 Domain averaged force over the voltage cycle for different frequencies .................130

22 Mean domain averaged force over the time cycle for different voltage .....................130
















LIST OF FIGURES


Figure page

1 Experimental photograph illustrating the EHD effects on a laminar jet of titanium
tetrachloride vapor .................. ............................. .... ........ .. ........ .... 11

2 Generation of glow discharge in various arrangements...............................................17

3 1-D Numerical modeling results showing different regions of the discharge gap at the
instant when the current is m maximum ........................................... ............... 18

4 A 3-D surface plot of photomultiplier light output from the plasma as a function of time
and chordwise distance, for a single AC cycle of the applied voltage...................21

5 Separation control on a NACA66-18 wing at a Re=79000, 16 degrees angle of attack.23

6 Streamer front propagation at different instants in a DC voltage.............................35

7 Levels of plasm a m modeling .................................... .... ......... ................................... 38

8 Parallel plate arrangement for modeling DBD discharge ......................................76

9 Evolution of discharge current over one period of the discharge cycle.....................76

10 Evolution of the gap voltage over one full cycle ..................................... ..................77

11 Comparison of peak current species density.........................................................78

12 Schematic illustration of the plasma-fluid flow modeling framework.........................83

13 Schematic of the actuator arrangement with approximate shape of electric field lines.
The thickness of the electrodes d is 0.1 mm ..................... ............................... 86

14 The line A-B constitutes the plasma fluid boundary using linear approximation. The
electric field strength outside this line is not strong enough to ionize the air..........89

15 The computational domain used for calculations shows the flat plate with the upper
e le ctro d e .......................................................................... 9 3

16 Velocity profiles compared at the four different stations ..........................................95

17 Effect of the various parameters as seen form the velocity profiles at ST4 ................96









18 Effect of the body force on the wall shear at ST4 ....................................................96

19 Temperature profile along the wall for Re=1027 and Pr=0.7................... ..............97

20 Comparison of normalized heat flux for different voltages and frequencies................98

21 Comparison of the experimental and calculated plasma force for different voltages. 101

22 A schematic of the computational domain with the flat plate and the electrode. The
electrode has negligible thickness in the above arrangement .............................103

23 Drag variation with angle of attack for various Re. The plasma is operated at 4KV rms
an d 3K H z ......................................................................... 10 3

24 Variation of lift with Re for different angles of attack. The plasma is operated at 4KV
rm s an d 3 K H z .................................................. ................ 10 6

Nusselt number variation with angle of attack for various Reynolds number. The plasma
is operated at 4KV rms and 3KHz frequency ................................ ............... 106

26 Streamlines showing flow separation for Re=333 and an angle of attack twenty
degrees for different positions of the electrode ............................................... 108

27 A representative 2-D asymmetric discharge arrangements .................... ........ 110

28 Sensitivity to initial condition (20KHz, 1KV applied voltage) ................................114

29 Photo intentisty measurements to depict the observed discharge evolution over a time-
cy c le ...................................... .................................................. 1 1 5

30 Time dependent plasma behavior (5 KHz, 1KV applied voltage)...........................115

31 Ion density evolution in the domain over time (5 KHz, 1KV applied voltage) ..........116

32 Force evolution over one time cycle (5 KHz, 1KV applied voltage) .........................117

33 Asymmetric sawtooth voltage waveforms used in the experimental and numerical
stu d ie s ....................................................................... 12 2

34 Experimental results comparing the discharge observed from the two different
saw tooth w av reform s.................................................................... .. ................... 123

35 Impact of waveform on the species number density...................................................123

36 Thrust measurements showing the more efficient.................... .................. ................125

37 Axial force integrated over the domain with time (Fx (dynes)) for different
w aveform s ..................................... ................................. .......... 125









38 3-Different configurations to study the impact of electrode overlap..........................126

39 Time evolution of species number density over one time-cycle (20KHz, 1KV voltage
sig n a l) ................................................... ................... ................ 1 2 6

40 Impact of the geometric configuration on the domain integrated force field (20KHz,
1K V voltage signal) ......................... .... ................ ............................ 127

41 Three different lower electrode lengths modeled to investigate the effect of increased
surface area on the discharge....................................................... ............... 128

42 Impact of lower electrode size ......................................... ......... ........................... 128

43 Charge density and force variation over time for different frequencies.................129

44 Charge density and force variation over time for different voltages ..........................130

45 Fx power-law dependence on voltage is compared for two different numerical modelsl31















NOMENCLATURE


n j, n, Species number density

N Gas number density

u v,, Velocity

Potential

E Electric field

T Temperature of various species

S,e, Species internal energy

7, Ratio of specific heats

/p,,eg Species mobilities

Di,e,j Species diffusivities

mi,ej Species mass

R,, Rd Creation/destruction of species

vi, VD Creation/destruction frequencies

Kc,KD v, /N, vDlN

Rij Number density gain/loss from reaction

Number density gain/loss from inelastic
collisions

VIj Elastic collision frequencies

qj Charge carried by jth species

e Elementary electronic charge










o, Permittivity of free space

K Boltzmann constant

Rm Momentum loss from collision

Vm,,,~ Momnetum loss collision frequency

Km Vm/N

Pe Electron pressure

Ren Energy loss term

Qe Electron heat conduction

ET Species thermal energy

r, Driving/operating time scale

T, Creation/ionization characteristic time

Td, Species drift characteristic time

Tdel Dielectric relaxation characteristic time

Ci Ion/species sound speed

Ce Electron sound speed

Cfluid Fluid sound speed

Ftave Time-averaged body force

Force resulting from species collision
olllosses
losses















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

COMPUTATIONAL MODELING OF GLOW DISCHARGE-INDUCED FLUID
DYNAMICS

By

Balaji Jayaraman

August, 2006

Chair: Wei Shyy
Major Department: Mechanical and Aerospace Engineering

Glow discharge at atmospheric pressure using a dielectric barrier discharge can

induce fluid flow and operate as an actuator for flow control. The largely isothermal

surface plasma generation realized above can modify the near-wall flow structure by

means of Lorentzian collisions between the ionized fluid and the neutral fluid. Such an

actuator has advantages of no moving parts, performance at atmospheric conditions and

devising complex control strategies through the applied voltage. However, the

mechanism of the momentum coupling between the plasma and the fluid flow is not yet

adequately understood. In the present work, a modeling framework is presented to

simulate athermal, non-equilibrium plasma discharges in conjunction with low Mach

number fluid dynamics at atmospheric pressure. The plasma and fluid species are treated

as a two-fluid system exhibiting a few decades of length and time scales. The effect of

the plasma dynamics on the fluid dynamics is devised via a body force treatment in the

Navier-Stokes equations. Two different approaches of different degrees of fidelity are









presented for modeling the plasma dynamics. The first approach, a phenomenological

model, is based on a linearized force distribution approximating the discharge structure,

and utilizing experimental guidance to deduce the empirical constants.

A high fidelity approach is to model the plasma dynamics in a self-consistent

manner using a first principle-based hydrodynamic plasma model. The atmospheric

pressure regime of interest here enables us to employ local equilibrium assumptions,

signifying efficient collisional energy exchange as against thermal heating from inelastic

collision processes. The time scale ratios between convection, diffusion, and

reaction/ionization mechanisms are O(107), making the system computationally stiff. To

handle the stiffness, a sequential finite-volume operator-splitting algorithm capable of

conserving space charge is developed; the approach can handle time-step sizes in the

range of the slowest species convection time-scale. The Navier-Stokes equations

representing the fluid dynamics are solved using a well-established pressure-based

algorithm. A one-dimensional two-species plasma model was employed as a test case for

validation purposes. The momentum coupling is primarily caused by the combination of

factors which include discharge chemistry, individual species transport properties,

geometric construction and the nature of the insulator and electrode material. Overall, the

paraelectric momentum coupling mechanism is due to the cumulative effect over time of

the force field in the domain, as seen from our computations.

Parametric studies conducted on the operating variables such as voltage. Frequency

and geometric arrangements indicated strong agreement with the observed experimental

work. The applied voltage indicated a power-law dependence on the voltage for the

measured force in the domain.














CHAPTER 1
INTRODUCTION

This dissertation deals with the computational modeling of discharge processes and

the accompanying plasma-fluid interaction with focus on the asymmetric electrode

dielectric barrier plasma actuator and its applications. In particular we focus our attention

on the operating mechanism and thermo-fluid effects of the discharge plasma actuator

(see Figure 1) which generates efficient surface plasma using a dielectric barrier

arrangement capable of interacting with the near-wall flow structure through momentum

injection. In this chapter we will present (i) a brief overview and background of the

discharge operating processes in the actuator (ii) the mechanism of operation of the

plasma actuator and momentum injection to the flow (iii) the state of the research of the

discharge-induced fluid mechanics and their numerical modeling. Finally we will present

an outline of the various issues addressed and the overall organization of the present

dissertation.

1.1 Overview of Plasma Generation and Dielectric Barrier Discharges (DBD's)

The plasma is defined as a quasineutral particle system in the form of gaseous or

fluid-like mixtures of free electrons and ions, frequently also containing neutral particles

[1]. Based on the extent of ionization, one observes either fully ionized or weakly ionized

plasma. Typically, such ionization is realized in gases when a sufficiently large potential

is applied under conditions of low pressures such as in vacuum discharges [1]. However,

the advent of dielectric barrier discharges (or in short DBD's) presents an efficient mode









of volume plasma generation at atmospheric pressure conditions at which the potential

for applications are enormous.

Dielectric barrier discharges have been studied since mid-nineteenth century when

they were identified in the context of ozone generation (Siemens [2]). Since then, they

have been the focus of extensive interest generated by their use in a wide variety of

applications such as surface treatment, excimer ultraviolet lamps [3], large area flat

plasma display panels [3] etc. Recent applications have taken this a step further by using

the concept of the DBD to generate surface electro-hydrodynamic (EHD) flows [4] which

have potential for flow control, thermal management among other applications.

DBD's are very attractive for industrial applications because they can provide non-

equilibrium plasma conditions at atmospheric pressure. Typically, discharges at high

pressure tend to develop into filamentary or arc type discharges which are unsuitable for

certain applications. However, in recent years, homogeneous glow discharges have also

been obtained at atmospheric pressure which can be used for large volume applications,

surface treatment. The DBD provides a mechanism for generating uniform discharge at

atmospheric pressures without arc formation. Typical discharge configurations use planar

or cylindrical arrangements with a dielectric layer (e.g. glass, quartz, polymers) placed

between the electrodes. The discharge gaps usually range from less than a millimeter to a

few centimeters depending on the application [3]. The most interesting property of these

discharges is that, at about atmospheric pressure, the breakdown is initiated in the form of

a large number of independent current filaments which are known as microdischarges.

These discharges are characterized as weakly ionized plasma channels which resemble

high pressure glow discharges. Due to the charge build-up at the dielectric surface in a









short time after breakdown, the electric field in the discharge region is reduced, thus

making it self-limiting. The short duration of the discharge current limits the energy

dissipation and consequently results in little gas heating (hence, the plasma being

athermal) and reduces power consumption. In recent years, homogeneous glow

discharges have been generated at atmospheric pressures in DBD configurations [4-7].

Such a discharge generation mechanism has been used by Roth [5] to devise a source of

uniform and stable surface plasma. It has been known that EHD is one of the efficient

modes of flow control in an ionized gas [4]. The EHD induced flow is schematically

illustrated below in Figure 1. The figure shows the effect of the plasma generated EHD

flow on a surface using smoke flow visualization as reproduced from [4]. The dominant

effect on the flow can be seen from the bending of the jet and re-emerging as a wall jet.

This can be easily used to inject momentum in regions of adverse pressure gradients,

boundary layers etc. Significant attention has been devoted to boundary layer flow

control using methods such as synthetic jets [8]. The present method has the advantage of

non-moving parts, small size, performance at atmospheric conditions and capacity for

complex unsteady flow actuation. But, issues such as power consumption, force

generation mechanism and optimum operating conditions have still not been adequately

addressed. In this dissertation, we will focus on developing numerical models (both first

principle-based high-fidelity as well as phenomenological approaches) to study the

momentum generation and the thermo-fluid potential of these devices.

1.2 Plasma Actuator Operating Mechanism

The exact mechanism of electrohydrodynamic (EHD) flow generation is not

sufficiently understood although the concept behind the force generation is widely

perceived to be that of Lorentzian collisions. The paraelectric effects arising from the









electric field gradient accelerate the ions which transfers momentum to the neutral fluid

[4]. The effect of the plasma on the fluid is like a localized body force [9] on the neutral

particles. Most of the investigations so far have been based on experimental observations

and empirical arguments using approximate models [4,9] until recently. The above

mentioned EHD effects involve the interaction between the discharge physics and the

near-wall flow physics. In order to systematically explain the plasma fluid interaction,

detailed modeling studies are required. Such studies should attempt to address the

complex plasma physics and its interaction with flow physics. The coupled plasma-fluid

problem is inherently non-linear by nature, just because one is forced to treat the fluid

and plasma as separate entities which interact. It will be shown later in chapters 3 and 4

that the disparity in the fluid and plasma time scales can be used to decouple their

solution methodologies. Moreover, the EHD flow can itself be of a different length scale

as compared to the mean flow. Resolving these multiple scales is a challenging

computational task. Also, the modeling of plasma dynamics in itself involves a

convection-diffusion-reaction system operating over widely varying time and length

scales.

The interaction between the aerodynamic motion of an electrically conducting

medium and an electromagnetic field is of great interest in astro-physics; interstellar gas

masses etc and is studied under magneto-aerodynamics. The potential of electromagnetic

effects in enhancing aerodynamics performance have been discussed and reviewed by

Shang et al. [10].

Mechanism of flow control using plasma studied in non-equilibrium hypersonic

flows is fundamentally different from the EHD type mechanism described above. Here,









the ionized flow modified in the presence of electromagnetic fields represents the

interplay of aerodynamics, electromagnetism, chemical physics, and quantum physics.

This is encountered in hypersonic flights where, the air mixture bounded by the bow

shock consists of highly energized internal energy modes and at high temperatures

around 5000 K ionization is encountered. This ionized air mixture will interact with the

electromagnetic field causing Lorentz force effects to affect the aerodynamic forces. This

Lorentz force becomes significant at high altitudes and high speeds which are used to

accelerate or decelerate a gas continuously at subsonic or supersonic speeds without

choking even in constant cross section channels. However, this is achieved by coupling

the velocity and the temperature of the conducting media through Joule heating effects. In

other words, the Joule heating from the electromagnetic effects can significantly alter the

thermal equilibrium. Research on shock wave propagation in plasma has been studied by

Rivir et al. [11] and high speed flow control using weakly ionized plasma by Leonov et

al. [12]. The studies reveal that at high speed non-equilibrium flows the weakly ionized

gases can substantially modify a traveling shock wave. A propagating shock wave in a

plasma decreases in strength and the wavefront disperses. Significant aerodynamic drag

reduction via plasma injection in the stagnation region has been reported in the studies of

Rivir et al. The origin of this drag reduction is not completely understood, but the non-

equilibrium thermodynamics and the electromagnetic forces play significant roles. The

modification of the flow field structure due to non-equilibrium thermodynamic

phenomena is caused in two ways. The energy from the weakly ionized plasma energy

source modifies the local gas temperature and reduces the Mach number upstream of the

bow shock. Secondly, the non-equilibrium energy distribution among internal degrees of









excitation can also alter the thermodynamic properties within the shock layer. Also, in

low supersonic flows, the high amounts of energy deposition upstream of the shock can

make the flow locally subsonic. The modeling of the plasma effects in non-equilibrium

effects also involves significant complications as intrinsic plasma properties such as

transport properties, conductivity etc. need to be calculated from kinetic modeling studies

and not from the equilibrium Maxwellian distribution assumption used in fluid models.

However, the focus of the present research is the low Mach number flow regime where

non-equilibrium effects are minimal and, moreover, the plasma dynamics is mainly

shaped by the electromagnetic effects and to a lesser extent by the background low-speed

neutral flow.

1.3 Numerical Challenges

Computational studies of the dielectric barrier discharges have been attempted for

parallel plate arrangements with emphasis on modeling the plasma physics in detail [13].

The time and length scales operating in the plasma extend over a vast range, warranting

kinetic models in the extreme case (highly non-equilibrium, low pressure, high frequency

(RF) discharge conditions) to relatively simpler plasma-fluid models for equilibrium,

high pressure plasmas. The latter is still an approximate model, but valid in certain

situations such as at high pressures where equilibrium velocity distribution functions and

isotropy assumptions for pressure gradients can be used in spite of the large electrostatic

effects. Even modeling the evolution of the plasma discharge with reasonable complexity

using a fluid model is computationally challenging and expensive due to

* strong non-linearities between the plasma-fluid equations and the electrostatics

* widely varying spatial and temporal scales, especially near the electrodes.









Most of the computational work in the area of plasma modeling has been done for

high frequency discharges and for parallel plate arrangements which can be extended to

complex arrangements used in flow control studies. Keeping in mind the computational

cost, 2-D or higher dimensional discharge modeling studies are less common. More

recently Roy et al.[14-18] have developed a finite-element method based 2-D modeling

of such plasma-based actuators by solving for the discharge species and neutral transport.

Although the model used only three species (ions, electrons and neutrals) in a 2-D

framework, the computational time was significant (a week of computational time for a

typical actuator configuration). However, the plasma and the neutral fluid display three-

dimensional characteristics as reported from experimental observations. Also, the extent

of plasma chemistry to be included is a compromise between exact physics and

computational cost. Involving a large number of species will mean that we have to solve

many transport equations and consequently a hugely coupled system. However, we will

see in later chapters that the impact of the ion species (both positive and negative which

are largely responsible for the collision dominated momentum transfer) on the evolution

of the discharge and electric force field is significant as deduced from the discharge

structure.

1.4 Overview of Present Approach

The key to modeling the DBD effects in fluid dynamics is to achieve realistic

distributions of the species densities and their momentum in the domain which interacts

with the neutral fluid by solving the plasma fluid equations. As a first step to model these

effects we developed a phenomenological model [9] based on a linear force distribution

in the domain. This linearized body force model was later adopted by Gaitonde et al. [14]

for modeling plasma-based separation control in a NACA 0015 wing section. This model









was an attractive option because of the difficulty in achieving efficient multi-dimensional

and self-consistent plasma dynamics simulations coupled with largely unsteady fluid

dynamics.

This was followed up by a first principles-based self-consistent modeling of the

plasma-dynamics using a hydrodynamic model in helium gas. This multi-fluid

formulation to model the radio-frequency discharge in helium gas gives the spatial and

temporal evolution of the charges species which is decoupled from the neutral fluid

dynamics. The body force calculated from this data provides a more self-consistent way

of modeling the plasma-fluid flow interaction.

The wide range of operating length and time scales characterizing both the

discharge and flow physics resulting in a highly stiff system of equations and is a major

limiting factor in achieving fully-coupled multi-dimensional plasma-fluid simulations. In

this dissertation, we clarify some of the important computational issues pertaining to

handling the multiple scales, while offering insight into improving our computational

capability for modeling practical, multi-dimensional physics.

Both sequential and fully implicit approaches [15] have been used for discharge

modeling. In the present multiple scale problem, we have employed an operator-split

sequential solution algorithm in which the solution procedure can be adapted to handle

the individual processes efficiently and realize overall gain in computation. This method

is integrated with a multi-block finite-volume algorithm capable of handling 3-D

curvilinear grids [19]. The method is employed to model the plasma dynamics in an

asymmetric electrode configuration similar to that shown in Figure 1.









1.5 Outline of the Dissertation

The remaining document is organized as follows. In chapter 2, we present a

review of the various studies on discharge plasma to provide a reasonable understanding

of their generation and operating mechanisms. The first section is devoted to

experimental studies aimed at analyzing the discharge operating characteristics and

structure followed by some potential fluid flow control applications using the surface

plasma arrangement. A brief glimpse is provided into the possible potential of these

plasma actuators in a turbulent regime along with laminar separation control, drag

reduction etc. The next section reviews the various numerical modeling efforts of

discharges operating in high and low frequency regimes. These give a view of the various

versions of models employed to study the plasma discharge in one and two dimensions.

Further, we also look at some of the contemporary efforts at modeling the discharge

plasma actuators in 2-D.

In chapter 3, an overview of the various levels of plasma modeling is presented

starting from the studies which solve the Boltzmann equation in a kinetic model which

are numerically intractable for large-scale problems. This is followed by a discussion on

the hybrid methods which use both kinetic as well as fluid descriptions with the goal of

developing an efficient, feasible and accurate description of the processes. The last

section presents in detail the various features of fluid models including the various levels

of complexity that can be achieved and the set of valid approximations. Motivated by our

interest in efficiently handling the plasma equations, we present an overview of the

techniques employed to handle such multiple scale processes. A solution framework for

the coupled species transport equations and the electrostatics is presented. The wide

range of scales needs to be handled using time-split methods and operator splitting









techniques are employed for the source term treatment. With the large difference in the

plasma and neutral fluid flow time scales, the Navier-Stokes equations are treated

independently from the discharge model. However, the effect of the glow discharge is

included in the flow equations as a local body force term which can be evaluated from the

quasi steady solution of the plasma-fluid equations. Having developed the modeling

framework, a one-dimensional representative modeling study of a parallel plate glow

discharge in helium is attempted as a validation exercise against prior modeling [20].

In chapter 4, a body force type coupling is developed in detail to facilitate

interaction between the plasma dynamics and the fluid dynamics. A phenomenological

body force model based on experimental observations, developed independent of the 1-D

plasma model results, is used here to investigate qualitatively the induced flow structure.

As a sequel to the above, the fully self-consistent discharge dynamics developed from the

hydrodynamic plasma model is presented to provide insight into the discharge structure

and the force generation mechanism.

Chapter 5 summarizes the present research status and enunciates the future

research goals to be attained.























laminar jet of titanium
tetrachloride vapor .


-15 mm





(a) Plasma Off














(b) Plasma On, E ~ 4.5 kV rms, F = 3 kHz

Figure 1 Experimental photograph illustrating the EHD effects on a laminar jet of
titanium tetrachloride vapor. The top picture shows the unaffected jet when
the plasma is non-existent and the bottom picture shows the effect on the jet
which bends and forms a wall jet. Photograph reproduced from [4].














CHAPTER 2
LITERATURE REVIEW

Flow control through various mechanisms has been extensively researched in

literature. Recently, the application of surface glow discharge plasmas as flow control

devices has attracted particular interest as revealed in experimental studies undertaken in

the last few years [5,21-27]. However, numerical studies to complement the above

experimental efforts are not too common until recently [28-29,14-18]. The operating

mechanism of the high frequency discharges between dielectric barriers has been studied

using both experimental and numerical modeling techniques, in the context of surface

coating, sterilization, ozone production, plasma display panels (PDP's) [3,6,30,31] etc. In

the previous chapter we presented a brief overview of the dielectric barrier discharge

plasma actuator, its operating mechanism, potential for thermo-fluid applications and

some insight into the current status of the numerical aspect of the research and the

associated challenges. Here, we review in detail some of the relevant experimental and

numerical studies on discharge plasmas and fluid flows to help better understand the

significance of this dissertation. The first section 2.1 is devoted to experimental studies

pertaining to discharge actuator mechanisms and applications. The second section 2.2

reviews the various plasma dynamics modeling efforts and concludes with the numerical

studies attempted at modeling the discharge-induced fluid flows.

2.1 Review of Experimental Discharge Plasma Studies

In this section we will review some of the experimental efforts to modify flow

using the plasma based actuators and their operating mechanism. Initial studies of









Dielectric barrier discharges, characterized the by the presence of numerous micro

discharges which lead to filamentation of atmospheric air. In these studies, the frequency

of the applied voltage usually ranged from line frequency to a few hundred kilohertz.

This was first investigated in the context of ozone production. A detailed review of

various efforts undertaken is presented in [3].

2.1.1 Early Discharge Studies

It was Kanazawa et al.[5] in the late 1980s, whose work indicated that homogenous

glow discharges as against filamentary ones can be obtained under certain conditions.

Since then the potential of DBD's in the context of plasma display panels for applications

like flat television screens [3] etc. was realized, which led to spurt of research activity

aimed at furthering the understanding of plasma dynamics applications. The

experimental work on the plasma actuator development took off in the 1990s, when Roth

et al.[32] formally devised the generation mechanism of uniform or homogeneous plasma

discharges at atmospheric pressure which he called as the OAUGDP or One Atmosphere

Uniform Glow Discharge Plasma. The OAUGDP is capable of operating at atmospheric

pressure in air and other gases without significant generation of heat. It is non-

equilibrium plasma, far from thermodynamic equilibrium. This means that the individual

species such as the electron possess energies which are very far apart from the heavier

ions as well as the neutral species. The generated plasma is more uniform in the

breakdown region and attractive for surface and large volume applications. The active

species thus generated had potential for use in sterilization, surface decontamination,

increasing surface wettability and surface energy of materials [9]. A schematic of the

OAUGDP generator is given in Figure 2. This plasma generation device can be operated









in a wide range of geometrical configurations such as, between two parallel plates or as

surface plasma which can cover a flat or curved surface such as an aerofoil.

The asymmetric arrangement in Figure 2b is the source of the surface plasma over

an aerofoil, with the direction of the generated electro hydrodynamic (EHD) flow given

by the placement of the submerged electrode, i.e., the flow is always seen to be generated

from the upstream electrode to the downstream electrode. This EHD flow is caused by

the 'paraelectric' effect [4] the term given to the flow (say 1-10 m/s) created because of

the curved electric field lines as in an asymmetric arrangement shown above in Figure 2c.

Only an asymmetric arrangement as shown in Figure 2b generates usable wall jet type

flow. There are other symmetric arrangements studied in [4] which by virtue of

construction produce vortex flow as against wall jets. The momentum acquired by the

ions in such an electric field is coupled to the neutral fluid through Lorentzian collisions.

Roth et al.[33] identifies use of these EHD flows or its variants as control tools in

aerodynamic boundary layers, pumping gases through tubes or ducts etc. A few key

observations of the OAUGDP from these experimental studies [33] are

* Plasma is non-equilibrium and cold with time-dependent features.

* Discharge exhibits classical glow discharge characteristics with identifiable regions
as identified by computer modeling in [20,34]

* Ionization occurs at 'Stoletow' point with minimum possible energy.

By non-equilibrium, it is implied that the plasma does not have thermodynamic

equilibrium and the kinetic temperatures of the species are far from equal to each other.

This further implies that the background gas need not be heated till thermodynamic

equilibrium. Three prominent EHD based flow effects discussed in [4] are presented in

Table 1.









2.1.2 Mechanism of Plasma-Induced Flow Applications

In this section we will present efforts to address the paraelectric flow mechanism.

Experimental and analytical studies were conducted by Roth et al.[4,7,33]. Here, EHD

body force is modeled as the electrostatic force acting on the charged particles which acts

on the neutral gas. The generation of paraelectric effects as a consequence of electric

field gradients is illustrated in [32-33]. A brief description of the analytical argument of

Roth et al.[33] is given below. If E, is the electric field (V/m) and p, (C/cu.m) is the

charge density, then the local electric body force FE (N) according to Roth et al.[7] is

FE = Ep,. (N) (1.1)

Using Maxwell's equation, we get


V.E =c (1.2)
go

Hence,

1dF 1) (1.3)d
F, = sEVE = oV.(EE) = COE2 (1.3)
2 dx 2

for a one-dimensional case. This formulation for the body force allows viewing it as an

electrostatic pressure gradient where,

1 2
PE = E2 (1.4)
2

This simplified 1-D interpretation gives a feel for the nature of paraelectric body

force effect, but cannot be extended to higher dimensions. Moreover, evaluating the exact

distribution of E, requires solving for the charge particle distribution. Nevertheless, the

above electrostatic pressure argument gives a feel for the observed phenomenon. The

'paraelectric' effect can be interpreted as the pumping of charged particles towards









regions of increasing electric field gradient. Also, in a quiescent medium, the electrostatic

pressure acts as a driving pressure gradient causing the flow from a region of higher

pressure to that of a lower pressure as reported from the smoke flow visualization study

in [4]. More recently, Enloe et al.[21-22] studied, the plasma morphology and operating

mechanism in an asymmetric electrode setup as shown in Figure 2b by means of optical

measurements using a Photomultiplier (PMT) tube.

The experimental study revealed characteristic temporal and spatial structure.

Computer modeling efforts for parallel plate arrangements in [20,34] showed that the

DBD generated in those conditions possessed characteristic low pressure glow discharge

qualities with regions such as the Faraday dark space, cathode glow, positive column etc.

This is illustrated in the 1-D modeling study by Lee at al [34] as shown in Figure 3. The

figure shows the distribution of the electric field and species densities at the instant of

peak current which are characteristic of glow discharges. The temporal structure of the

DBD generated in an asymmetric arrangement indicated noticeably different discharge

patterns in the two half cycles. The light intensity used as a measure of the discharge was

consistently regular when the exposed electrode was negative-going and highly irregular

when positive-going. Further, it was noticed that a consistently increasing voltage form

was necessary to sustain the discharge and the resulting light intensity depends heavily on

the rate of increase. Also to study the asymmetry effects further, asymmetric saw-tooth

waveforms having vastly different positive-going and negative-going duty cycles were

used to model the discharge. The resulting outcome clearly established the preferential

phase for the actuator operation during a whole time-cycle.















DIELECTRIC
N


(a) parallel plate setup




Incoming
flow
I


.1.4.41+4.C


Flow bei



-I


exposed
electrode
substrate


ELECTRODES


HIGH
VOLTAGE
POWER
SOURCE


nds


nds
plasma



Outgoing wall jet



insulated
electrode


AC voltage
Source
(5-10 KV, 1-20KHz
[4])


insulatio


(b) A single asymmetric electrode setup and its effect on fluid flow.


(c) Illustration of distorted field lines in an asymmetric arrangement
Figure 2 Generation of glow discharge in various arrangements such as between parallel
plates (a) or as a surface discharge (b) and the curved electric field lines in
such a case [9].


Lower
electrode







18


Table 1 Representative EHD operating mechanisms
Type of EHD effect Max speed Mechanism Comments

Very large collision
Paraelectric body force frequencies (7GHz)
effect due to net charge 1-10 m/s Electric field gradients ensure adequate
density in the plasma momentum by
Lorentzian collisions.

Impractical due to high
DC electric field over a
EHD convection driven electric fields of the
by DC ion mobility dift 300 m/s generated OAUGDP order of 104V/
causes ion drift
requirement.
of par c Gives a higher induced
Series ofparaelectric
Peristaltic EHD plasma 100 m/s rs low velocity to work
100 m/s actuators operating at a .
acceleration p d rc with- most useful and
phase differenceof the three
flexible of the three


Note: See [33] for more information.










im


/
/'


cathoe


Figure 3 1-D Numerical modeling results from Lee et al.[34] showing different regions of
the discharge gap at the instant when the current is maximum for a 10KHz,
1.5KV voltage, a- cathode fall; b- negative glow; c- faraday dark space; d-
positive column

While the light intensity measurements of Enloe et al.[21-22], revealed the

significant temporal and spatial structure of the discharge in the actuator, it still did not

improve the understanding of the nature of momentum coupling. In a follow-up of the

above study, Baird et al.[35] investigated the acoustic effects of the plasma actuator. The

study revealed that the momentum transfer to the neutral fluid was different in the two

halves and the acoustic pattern produced was a radiation pattern characteristic of a


-' '' i '









coherently-driven system indicating possible compressibility effects. Near electrode fluid

density measurements over time by Enloe et al.[36], indicate density variation in space

akin to a pumping-type mechanism resulting in a wall-jet. However, hot-wire temperature

and velocity measurements by Jukes et al.[37-38] around the exposed electrode show

only a slight increase in ambient temperature of two degrees, thus reducing the possibility

of buoyancy driven effects.

In the present dissertation we proposed as part of our phenomenological model

(Shyy et al.[9]) that, the asymmetries introduced through the electrode arrangement and

the dielectrics along with the asymmetric discharge structure had a significant role to play

in the generation of EHD based paraelectric flow or simply glow-discharge induced flow.

It can be seen from Figure 4, that the distinct difference in the light intensity

measurements between the two half cycles for asymmetric electrode arrangement. The

first half cycle provides a more uniform discharge of high intensity while the other half

cycle provides an overall weak light intensity. Overall, the studies on the discharge

actuator summarily observed the following:

* Plasma structure is different in the two half cycles which was a fundamental
assumption in the development of our phenomenological model (Shyy et al.[9]) and
was also suggested in [21] (see Figure 4). It is possible caused by the geometric
asymmetry.

* Material asymmetry of the exposed and submerged electrodes only influenced the
net plasma intensity and a weak induced flow was present even with both
electrodes insulated.

* Direction of flow is purely determined by electrode geometry for the actuator
configuration discussed above.

* Mechanism of momentum coupling is yet to be adequately understood with
possible reasons ranging from density-variation effects from acoustic
measurements to contradicting temperature measurements indicating minimal
buoyancy related effects.









It was also perceived from the outcome of our phenomenological model [9] that the

net induced flow was a product of the various competing mechanisms listed above.

Further, it also highlighted the glaring necessity for detailed multidimensional modeling

of the actuator and plasma dynamics towards developing an accurate modeling strategy.

Parametric studies and the impact of the geometry on the actuator performance

were studied by Enloe et al.[22] and Van Dyken et al.[26]. Specifically, the impact of

voltage frequency, waveforms and geometrical parameters were investigated. Larger

thickness of the exposed electrode resulted in more thrust which indicated possible ion

concentration near the exposed electrode as they are the primary momentum imparting

collision species. Also, a larger thickness can mean a sharper electrode edge while

transitioning from the conductor to the insulator surface which can lead to stronger

electric fields in the region. The charge species concentration in the plasma will lead to

significant Debye shielding effects which can affect electric field distribution. In trying to

provide an analytical insight, Enloe et al, discusses the effect of Debye shielding which

further illustrates the inadequacy of Roth's electrostatic pressure argument.

2.1.3 Aerodynamic and Control Strategies Using Glow Discharges

Now that we have looked into some of the efforts to understand the plasma

actuator operating mechanism, we will look at its various areas of application. The above

discussed potential of the surface plasma to generate EHD flow effects whether wall jet

type flows or re-circulating flows can be harnessed for a variety of applications. The

injection of momentum into the boundary layer and delaying of separation has been

primary motives in aerodynamic surface flow control.









Weaker half
cycle









Stronger half
cycle


Figure 4 A 3-D surface plot of photomultiplier light output from the plasma as a function
of time and chordwise distance, for a single AC cycle of the applied voltage.
(Reproduced from Enloe et al.[23])

However, most of the efforts [4,25,39-40] listed indicate effective performance of

the plasma actuator from low (10,000) to moderate (100000) Reynolds numbers,

essentially limited by energy transfer from the plasma. Van Dyken et al.[26] use a single

plasma actuator at 7% chord on a NACA0015 airfoil resulting in substantial increases in

lift and reduction in drag beyond the baseline stall angle of attack. Santhanakrishnan et

al.[41] present the possibility of using plasma-based flow control for Unmanned air

vehicles (UAV's) to control stall and enhance lift, especially in the range of low to

moderate Reynolds number range(104-105). Goskel et al.[42] did an experimental

investigation of steady and unsteady plasma wall jets for separation control of Micro-air-

vehicles in the low Reynolds number regime (~ 20,000). It was observed that a gain in lift

by a factor of two and the loss of low Reynolds number hysteresis was achieved using

discharge power inputs in the range of 1.2mW/cm. Corke et al.[43], presents the idea of

using an unsteady plasma actuator to achieve flow control performance on an airfoil

comparable to that of a moving plane flap. This application is interesting for the reason

that the unsteady operation was observed to reduce the actuator power consumption up to










ninety percent. Related work includes studies by Corke et al.[24], Rivir et al.[11] and

Leonov et al.[12].

Table 2 Summary of experimental plasma actuator studies
Author and year Work description Comments


Roth et al(1998) [4]


Corke et al(2000) [24]


Corke et al I 1, '"2 [25]


Post et al 1,11' i) [27]


Roth et al.(2003) [39]


List et al 21' 211 ) [40]


Enloe et al 1 "4,1 [21-23]


Rivir et al, (2003) [45]


Rivir et al 121,i14 [11]


Leonov et al.(2003) [12]


Chan (2005) [44]



Enloe et al.(2005) [36]



Baird et al.(2005)[35]



Jukes et al.(2006) [37]


Boundary layer flow control
using plasma actuator arrays

Phased plasma arrays for
unsteady flow control

Flow control on a NACA 0009
wing

Separation control on a NACA
0066 wing
Flow re-attachment using
Paraelectric and peristaltic
EHD effects on a NACA0015
wing
Laminar separation control on a
cascade turbine blade
Experimental visualization of
the generated plasma structure
in space and time, effect of
geometric construction
Turbine flow control with
plasma

AC and Pulsed plasma flow
control

Effect of plasma induced
separation

Flow control using acoustic
effects


Temporal density, velocity
measurements near the plasma
actuator


Acoustic testing of plasma
actuator.


Characterization of discharge-
induced wall fow mechanism


Increase in drag. EHD flow
altered boundary layer profile.


Lift enhancement and also that
of drag (higher skin friction)
when typically using single
actuators.
Increased lift to drag ratio and
delays separation.

Peristaltic acceleration is by a
traveling electrostatic wave

Separation is quelled at
moderate Re-50000

Plasma structure helped gain
insight into the operating
mechanism



AC plasma controls flow
through wall jet mechanism
and pulsed DC Plasma through
heating mechanism
Uses plasma as a local heat
source
A symmetric plasma actuator
generates vertical structures
interacting with shear layer
instabilities.
Indicated density variation
slope near the electrode edge
as possible reason for
momentum coupling.
Acoustic patterns showed a
characteristic coherently
driven system indicating
compressibility effects.
The discharge actuator
characteristics are investigated
using hot-wire measurements
in the near wall region.









In [24] a plasma array with phased inputs with unsteady flow control is

investigated. Separation control on a NACA66 wing is studied by Post et al.[27] and is

illustrated in Figure 5. The plasma operation significantly reduces the size of the

separation bubble. In [12] the plasma generation is used as a local heating mechanism

thereby inducing boundary layer separation. Rivir et al.[ll] compare the paraelectric

surface flow control mechanism with the heating mechanism of Leonov et al.[12]. A

summary of experimental efforts is presented in Table 2.


Figure 5 Separation control on a NACA66-18 wing at a Re=79000, 16 degrees angle of
attack, reproduced from Post et al.[17]

A novel approach to apply plasma based EHD flow control to subsonic flows is

being investigated by Chan [44]. The acoustic properties of low rectangular cavities are

being investigated for the potential use of plasma as flow-control devices. In the









experiments, the electrodes of the plasma actuators were laid out streamwise to the flow.

A symmetric electrode arrangement (where the lower electrode is wider compared to the

exposed upper electrode and placed exactly below it to ensure symmetry) which can

generate vortical structures is used upstream of the cavity to inhibit the feedback process

in the shear layer which is key to the generation of large amplitude frequency tones. The

plasma actuators were reported to attenuate completely the discrete tones produced by the

cavity. The acoustic effects can act on the time-scale smaller than the fluctuating

frequency.

There have also been efforts to use plasma based actuators in the context of

turbulent boundary layer flow control. EHD based turbulent flow control concepts has

been presented in recent articles by Soldati et al.[45-47]. As in the present effort, the

electrohydrodynamic effects are modeled as a variable body force term in the equations

describing fluid motion. The complex interaction between the EHD effects, wall structure

of turbulence and the mean flow were studied. The EHD flows affect the turbulence field

by increasing both dissipation as well as production thus maintaining the balance.

However, overall drag reduction was observed with consistent decrease in the Reynolds

stresses. The ionization in these cases is through electrostatic precipitators and efforts are

being made to use OAUGDP in this context. Wilkinson [48] investigates the use of an

oscillating surface plasma wave for turbulent drag reduction. Jukes et al.[38] looks at the

turbulent boundary layer control for skin friction drag reduction using the surface plasma

actuator. A summary of turbulent flow control literature relevant to the present study is

tabulated in Table 3.









Table 3 Summary of turbulent flow control studies
Author and year Work description

Turbulent boundary layer control from EHD
Soldati et al, (2001) [47] flows. Drag reduction is achieved by
modifying the Reynolds stress contributions
Turbulence control using large scale EHD
Soldati et al, (1998) [45]flows
flows
u et a, ( ) [ Drag reduction in wall-bounded turbulence
using a transverse traveling wave
D k et ( ) [] Reducing turbulent wall friction through
Dhanak et al, (1999) [50] s e o
spanwise wall oscillations
Turbulent drag reduction using oscillating
.Wilkinon (2003surface plasma 40 % drag reduction was
Wilkmson (2003) [48] .
achieved, but ineffective at frequencies <
100 Hz
Juks et a 4) [] Turbulent boundary layer control for drag
Jukes et al.(2004) [38] rdto
reduction.


2.2 Review of Discharge Modeling Studies

In general, the numerical modeling efforts are developed in conjunction with

experimental studies to determine their fidelity. The complex physiochemical and

spatiotemporal processes associated with the glow discharge plasma also make numerical

modeling essential for their understanding and in modeling applications of plasma-based

devices. Experimental studies as reviewed in section 2.1 were limited in their ability to

reveal plasma structure and adequately explain the observed plasma dynamics and other

mechanisms. Specifically, the apparently contradicting arguments for the momentum

coupling based on compressibility effects proposed by Enloe et al.[35] as against the

observed lack of significant bulk heating due to the plasma reported by Jukes et al.[38],

highlight the necessity for a high-fidelity first principles-based numerical modeling

approach. Also, most of the preliminary modeling attempts were empirical [5,22] and

argumentative in nature. In order to establish a competent mechanism for the paraelectric









and other induced flow effects observed, detailed simulation of the plasma dynamics in

an asymmetric electrode arrangement is necessary. In this section we will review some of

the plasma and plasma actuator modeling efforts to establish the background for the

present dissertation. Section 2.2.1 consists of the various phenomenological/simplified

modeling approaches employed by various investigators to model the discharge-induced

fluid dynamics while, in section 2.2.2 we will present some of the more high-fidelity

discharge modeling studies.

2.2.1 Phenomenological Modeling of the Plasma Actuator

Most of the preliminary plasma actuator modeling efforts to study the fluid flow

coupling was either simplified phenomenological approaches or models with reduced

complexity. In the present dissertation we develop our own phenomenological model to

study the plasma-induced fluid flows. This model involved the development of an

empirical approach to calculate the time-averaged force to study the momentum coupling

with the neutral fluid. The various discharge parameters were assumed spatial or

temporal variations with the closure of the problem related constants requiring

experimental data for validation. Our analytical-empirical model (see Chapter 4) based on

a linear force distribution in the domain was later adopted by Gaitonde et al. [14] for

modeling plasma-based separation control in a NACA 0015 wing section. This linearized

body force model was an attractive option because of the difficulty in achieving efficient

multi-dimensional and self-consistent plasma dynamics simulations coupled with the

fluid dynamics.

Hall et al.[51] use a potential flow approach to model the aerodynamic effects of

the plasma actuator. In particular, the actuator mechanism is modeled as a doublet

(source-sink) to mimic the pumping action of the discharge on the neutral fluid.









Experimental validation was again required as a closure to determine the doublet

strength. In a different approach, Orlov et al.[52], employ a lumped-element circuit

model to account for the spatial and temporal variation of the discharge and mimic the

plasma body force effect on the neutral fluid. While the above studies showed qualitative

and reasonable qualitative correspondence with the experimental observations, they were

limited in their ability to confidently predict as well as explain the observed mechanisms

and the plasma-fluid interaction.

2.2.2 First Principles-Based Plasma Modeling Approaches

Most of these studies were undertaken either as part of studying dielectric barrier

discharges (DBD's), between parallel plates operating in the kilohertz regime (low

frequency) or the capacitively coupled radio frequency (high frequency) discharges for

plasma applications with special emphasis on the plasma actuator. The knowledge of

parallel plate discharge modeling is essential for gaining insight and developing a

numerical model for the plasma actuator. This is divided into two subdivisions- the first

one review the numerical challenges and approaches in the modeling of DBD's in

general, while the latter places sole emphasis on the discharge-coupled fluid flows which

are directly relevant to the present dissertation.

2.2.2.1 Review of dielectric barrier discharge (DBD) modeling studies

The modeling of low frequency (KHz) barrier discharges was attempted in the

nineties, when the potential of uniform glow discharges was realized and experimental

studies related to flow control and other surface processes gained momentum. In the

various studies reviewed here, the focus is on parallel plate glow discharges controlled by

a dielectric barrier in gases such helium, neon, argon, nitrogen and even air [53]. By low









frequency we mean a range 50Hz to a few KHz which is much smaller than the radio

frequency regimes (MHz).

The modeling of glow discharges involves the consideration of mechanisms such

as the reaction kinetics and plasma dynamics. With the high degree of complexity

involved, the choice of model is usually simplified tailoring to the application. An

overview of plasma models of varying complexity will be presented in the next chapter.

The approximation chosen in different studies are mostly based on choice and

convenience, resulting in each model differing slightly from the other while retaining the

original qualitative behavior. One of the first numerical studies of dielectric barrier

controlled discharges at atmospheric pressure was studied by Eliasson et al.[54] using a

two-dimensional fluid model for the plasma. Almost all modeling studies employ a fluid

model as against solving a Boltzmann equation because of the relatively easier nature of

solution technique for the former. The fluid model consists of the transport equations for

the charged and neutral species concentration, momentum and energy. The first three

moments of the Boltzmann equation are considered for modeling purposes. The

momentum and energy equations are also solved depending on the complexity of the

model.

The numerical solution consists of three steps which all are interconnected. The

first step is the solution of the Boltzmann equation to obtain the electron energy

distribution function in the gas under investigation. This is essential, for the properties of

electron transport such as mobility; diffusivity, mean energy etc vary as a function of E/N

(E is electric field and N is the gas number density). The above transport properties are

obtained solving a steady state Boltzmann equation under a constant field condition. This









assumption is reasonable because the characteristic time for energy equilibration under a

changing electric field at atmospheric pressure is 1011s while the characteristic time for

electric field fluctuation is 10-9s [55]. This makes it possible to assume energy

equilibration to be an instantaneous process. The solution of the Boltzmann equation for

the electron energy distribution and other transport properties has been tabulated over a

wide range of E/N and can be used during the next phase which is the transport module.

In this second step, the transport equations of species number density, momentum and

energy are solved.

The electric field E is solved from the Poisson equation


V.E= (1.5)
g0

where p is the net charge density. Also, the charge build up on the dielectric as the

discharge evolves and the external circuit should be added to the model at this stage.

Discussed above is brief plasma modeling framework using fluid models. Eliasson et al,

study a two dimensional model of the dielectric barrier discharge in three gases- Xenon,

nitrogen and oxygen. Most modeling studies use inert gases in their studies so as to

obtain stable glow discharges. Golubovskii et al, [13] use a multi-species 1-D fluid model

to study the different modes of operation of a helium discharge at atmospheric pressure.

The dielectric barrier discharges were seen to operate in two modes- a glow discharge

mode, displaying characteristic low pressure discharge structure and a Townsend mode

characterized by multiple current or discharge peaks. The formation of the various modes

is governed by design parameters such as gap width, dielectric barrier thickness and

applied voltage. Specifically, thin barriers and wide gaps led to a glow discharge mode

and a Townsend mode was realized otherwise. Similar behavior was reported in the









experimental studies by Mangolini et al, [56]. Two-dimensional model studies were

performed by Golubovksii et al., [57] Massines et al., [58] and Xu et al., [59]. All these

studies included substantial plasma chemistry to account for the exact species

concentrations. However, the significant aspect of these studies was it revealed

significant two-dimensional structure for the discharge. Xu et al, [59] state that the radial

expansion for the discharges can significantly alter the built up charge distribution on the

dielectric surface. This has been observed in a number of two-dimensional studies such as

Xu et al.[59] and Golubovskii et al, [60]. Mangolini et al.[56] state that the multiple

pulses as seen in a Townsend mode are manifestations of the radial inhomogeneities. The

radial expansion of the current pulse causes a non-uniform surface charge distribution

with the next breakdown occurring at the discharge periphery. Another important

modeling aspect is the consideration plasma chemistry. The treatment of plasma

chemistry varies from the most complex reaction mechanisms which give rise to

conservation and handling of multiple ion species, metastable atoms or molecules. The

Penning ionization effect also contributes to the complex plasma chemistry. The simplest

treatment of plasma chemistry is the two species model as used by Ben Gadri et al.[20]

Shin et al,[61] with more complex multi species treatments by Mangolini et al.[62]

Massines et al.[6,58] Tochikubo et al.[63] etc. Mangolini et al.[62] state that the

metastable atoms in helium discharge contribute significantly to the pre-ionization

through penning ionization mechanisms. This pre-ionization is important for the next

breakdown. Handling more number of species implies higher computational cost and is

usual kept at a possible minimum. The choice of boundary condition for the transport

equations differs in various studies which try to approximate the species behavior in the










sheaths. Golubovskii et al.[60] discusses the nature of interaction between the charged

particles and the dielectric boundary surface in the context of a homogeneous barrier

discharge in nitrogen. The influence of different mechanisms of electron emission and the

recombination are studied. A detailed discussion of the different numerical modeling

strategies, treatment of boundary conditions etc. are presented in later chapters. A

summary of the different low frequency numerical studies are tabulated in Table 4.

Table 4 Summary of low frequency plasma modeling studies
Author and year Work description Model
Modeling of dielectric Multi-species 1-D fluid model.
Golubovskii et al.(2003) [13] barrier discharges in helium Investigated various modes of
at atmospheric pressure barrier discharge operation
Golubovskii et al.(2004) [56] Modeling of nitrogen Multi-species 2-D model

Modeling ofDBD's using a .
Dongsoo Lee at al, (2004) [34] eling of D uing 1-D multi-species fluid model
helium-oxygen mixture
DBD modeling using He-Ar
Massines et al, (1998) [6] D modeling using He-Ar 1-D two-species fluid model
mixture
DBD modeling using He-Ar
Massines et al, (2000) [58] m e u 1-D multi-species model
mixture
DBD modeling of pure .
Tochikubo et al, (1999) [63] DBD modeling of pure 1-D multi-species model
Helium
.S \Modeling of Helium-
Mangolini et al, (2004) [62] Modeling of Helium- 1-D multi-Species model
oxygen mixture
Shin et al, (2003) [61] Modeling of He-Ar mixture 0-D, 2-species model
Modeling of nitrogen based .
Xu et al, (1998) [59] i i e 2-D multi species model
mixture
Eliasson et al, (1991) [54] Modeling in Xenon 2-D multi-species model


High frequency modeling studies have been the focus of numerical model

development for plasma simulations, in view of their extensive applications in material

processing, plasma display panels (PDP's) etc. While the nature and physics of these RF

discharges are different from what we are interested, the study of these discharges is

useful for their numerical modeling. Prior to the last decade, the simulations were

basically one-dimensional either with a fluid model or a particle model of the plasma.

Since then, two and three-dimensional studies have been successfully presented. A









review of one-dimensional studies is presented by Govindan et al.[64]. The dominant

type of plasma considered in these studies are low pressure, low-temperature and weakly

ionized. High pressure glow discharges have also been modeled lately by Yuan et al.[65].

The fluid model is based on the self-consistent solution of the electron and ion continuity

and momentum transfer equations coupled with the Poisson equation for solving the

potential between the electrodes. Boeuf et al.[66] assume a local equilibrium between the

charge particle kinetic and the electric field. This local equilibrium hypothesis serves as a

closure relation in representing the actual Boltzmann equation through the continuity and

the momentum transfer equations. Although, numerical modeling efforts to solve for the

Boltzmann equation have been attempted [67] they are less preferred to simpler fluid

models. Comparison studies have been made between the particle and fluid simulations

by Nitschke et al.[68] to study the range of validity of the fluid approximation for the

plasma. At higher pressures, the fluid model and the local equilibrium turn out to be

pretty good approximations for the RF plasma. More detailed discussion of the above

will be presented in later chapters. However, fluid models of varying complexity have

been summarized in Table 5. Solution of the transport equations in the fluid models is

computationally intensive because of their very stiff nature in space and time. To

illustrate this we will look at the streamer propagation modeling in a DC voltage by Dhali

et al.[69]. The propagation of the streamer obtained by a 2-D modeling study is shown in

Figure 6. The wave-like solution and sharp gradients are characteristic of streamer

propagation and also low frequency discharges. The system of equations is also highly

non-linear with the modeling of multiple species transport equations along with the










Poisson equation for the electric field. The various numerical techniques used to handle

these issues will be reviewed in the following chapters.


Table 5 Summary of high frequency plasma modeling studies
Author and year Work description
Boeufet al, (1987) [66] 2-D fluid model

2-D fluid and particle model
Nitschke et al, (1994) [68] comparison
comparison

2-D fluid model for Gaseous
Lymberopoulos et al, (1995) [70] electronic conference (GEC)
reactor


Campbell et al, (1995) [30]



Boeufet al, (1995) [31]



Meunier et al, (1995) [71]
Veerasingham et al, (1996) [72]



Yuan et al, (2003) [65]


2-D fluid model for a PDP


2-D fluid model for Gaseous
electronic conference (GEC)
reactor comparison with
experiments


1-D fluid model for PDP's


1-D fluid model for high
pressure helium glow
discharges


Comments
2 moment equations,
two specie
3 moment equations
are used in fluid model
and two species
Multi species plasma
chemistry and three
moment fluid
formulation
Multi species plasma
chemistry and two
moment fluid
equations
Two species plasma
chemistry and three
moment fluid
equations
Multi species plasma
chemistry and two
moment fluid
equations
Multi species plasma
chemistry and three
moment fluid
equations


2.2.2.2 Review of plasma actuator modeling studies

While numerical studies on discharge physics have been popular, computational

modeling of the discharge actuator and its physics have been attempted only recently [14-

18,29,74] in parallel with our ongoing efforts. The modeling the DBD effects in fluid

flows requires realistic distributions of the species densities and their momentum in the

domain which interacts with the neutral fluid by solving the plasma fluid equations. This

will in turn enable us to calculate the net electric force acting on the charge species and

acts as a body force on the neutral fluid. More recently, Roy et a1.[15-16] proposed a









self-consistent two-dimensional DBD fluid model for helium gas with application to

separation control using finite element techniques. This multi-fluid formulation to model

the radio-frequency discharge in helium gas gives the spatial and temporal evolution of

the charges species which is decoupled from the neutral fluid dynamics. The body force

calculated from this data provides a more self-consistent way of modeling the plasma-

wall jet interaction. This study employed a finite-element method and solved the various

species transport equations using a globally implicit procedure where the system of

charge continuity and momentum equations are assembled as part of a global matrix to

solve for the solution vector. Singh et al.[17], in a related paper, present a parametric

study of the different conditions in an asymmetric discharge configuration. Specifically,

it was observed that the net body force production in the domain over a whole time-cycle

produced a downward positive force for a typical configuration. Kumar et al.[18] study

the nature of the discharge and the resulting force field in the presence of a magnetic

field. Also studied is the shape effect of the electrodes in the event of a finite electrode

thickness. All the above-mentioned modeling studies assume negligible electrode

thickness. In the event of a finite-electrode thickness, the treatment of the dielectric-

electrode edge can impact the near-wall force field. The common approach to be

employed in the present study as well as in the previously mentioned modeling studies is

that the plasma species and the neutral fluid are treated as a two-fluid system coupled

through dynamic forces and pressure interactions. More recently, the simulation of the

asymmetric electrode discharge actuator has been extended to atmospheric air chemistry

instead of the Helium gas [75]. Visbal et al.[76] explored the potential of these plasma

actuators in the control of turbulent and transitional separated flows by means of the










discharged-induced near wall effects. Specifically, a pulsed actuator producing an

unsteady discharge was employed to highlight the role of turbulent and transition causing

mechanisms in drag reduction as against a pure wall-jet momentum injection. A summary

of the various modeling studies have been listed in Table 6.

I ---. ---I-- ,--


Electron
number density
(per cu.cm)


Axial distance (0-5mm)
Figure 6 Streamer front propagation at different instants in a DC voltage (Dhali et
al.[69])

Table 6 Summary of plasma-induced fluid flow modeling studies


Author and year


Roy et al.(2005) [15-16]


Singh et al, (2005) [17]


Kumar et al, (2005) [18]


Singh et al, (2006) [75]



Visbal et al, (2006) [76]


Work description


Self-consistent discharge modeling
using hydrodynamic plasma model.

Parametric analysis of different
discharge operating parameters and
geometric effects
Discharge and fluid flow evolution in
the present of a magnetic field and
shape effects of the electrode edge
Discharge modeling extended to air
chemistry and the corresponding flows
effects were studied.
Control of transitional and turbulent
separated flows using
phenomenological and first principles
based models.


Model
Employed a 3-species
model- ions, electrons
and neutrals for Helium
discharge.

Same as above.


Same model as above.

8-species model with
atmospheric air
chemistry














CHAPTER 3
NUMERICAL MODELING OF GLOW DISCHARGE PLASMA

In chapter 2 we reviewed some of the relevant experimental and numerical efforts

in the modeling of the discharge actuator effects for fluid flow applications. We also

highlighted some of the numerical issues arising from the multi-scale nature of the

problem and how it can limit the computational capability. The presence of these large

ranges of scales operating in the problem necessitates the use of modeling strategies

which can range from kinetic level models to the more common hydrodynamic/fluid type

models. In this chapter, we present a glimpse into some of these different modeling

approaches and in the process, present an argument for the choice 'fluid model' used to

describe the plasma dynamics in the present dissertation. In the second half of this

chapter, we present in detail, the various numerical issues associated with fluid model

and the presently developed solution algorithm which can efficiently handle the plasma-

fluid dynamics. Lastly, a one-dimensional plasma dynamics simulation to model the

discharge evolution between parallel plates is presented for validating the present

capability.

3.1 Plasma Modeling Hierarchy

Modeling of glow discharge plasmas has been researched considerably more in the

context high frequency (radio frequency (MHz) or KHz range) discharges than low

frequency (usually smaller than KHz range) atmospheric pressure discharges. However,

the studies of radio frequency (RF) discharges do provide insight for handling the stiff,

multi-scale and highly non-linear plasma equations coupled with the electrostatics. The









fundamental numerical model describing the high and low-frequency discharges is

essentially the same with the lone difference coming from the operating frequency of the

driving voltage. While this may impact the nature of the discharge, it retains the same

inherent physical characteristics in terms of length and time scales. The following

discussion is equally applicable to both high and low frequency modeling studies

although most of the work conducted in this dissertation focused on efficiently modeling

discharge plasma actuators in the KHz regime.

Computational modeling of discharge plasmas can be achieved using a variety of

approaches. In general, we can have

* a fluid treatment

* a kinetic treatment

* a hybrid scheme which combines aspects of the fluid and kinetic schemes.

In this chapter we will deal with the different types and levels of plasma modeling. The

hierarchy of models based on complexity will be in three areas namely,

* type of model (fluid or kinetic etc.)

* treatment of plasma chemistry

* space complexity (ID, 2-D).

The various levels of plasma modeling are illustrated in Figure 7. The chapter is

organized as follows. A brief overview of the three broad classifications of plasma

models will be presented with emphasis on the fluid model which is of interest in the

present work. The modeling complexity arising from plasma chemistry and spatial

treatment will be discussed in the context of the fluid model.











Glow discharge modeling



Kinetic model



Boltzmann equation Particle models
solving for velocity distribution Involving classical equations of
function motion for the particles in a
force field

Fluid model (valid at high
pressures) Monte-Carlo Particle-in-cell
Mass conservation collision model (PIC) technique
Momentum conservation (MCC) Based on
Energv conservation Based on stochastic deterministic
or probabilistic classical mechanics
modeling



Hybrid model
Combination of the fluid
model along with one of
the kinetic models.


Figure 7 Levels of plasma modeling

3.1.1 Kinetic Modeling Approach

In this section we will give brief overview of the kinetic models and how they

compare with the fluid model employed in this work. Kinetic models typically involve

the solution of the Boltzmann equation for the species velocity or energy distribution

function in both space and time or the particle simulations using techniques such as

Monte Carlo methods. Although the Boltzmann equation can provide much more

information than the continuum fluid models, its solution is computationally intensive

requiring an order of magnitude more computer time than the fluid calculation as kinetic

models are inherently higher dimensional and the solutions are functions of the velocity,

space and time, the additional degrees of freedom being described in the full phase space









of the problem. For example, Riley et al.[77] use a 1D-2V formulation which implies

solving for the particle velocity distribution function (VDF) in a single spatial direction

and in two directions of the velocity space. This is in contrast to the fluid models where

we assume an equilibrium distribution function for the species.

The coupling between the Boltzmann equation and the Maxwell's equations for

electrostatics is a difficult numerical problem as the particle dynamics should include the

strong body force effect of electrostatics. The solution variable in the kinetic models is

the species velocity distribution function. By assuming an equilibrium distribution

function in the case of fluid models, the full Boltzmann equation can replaced by its first

few moments corresponding to the conservation of mass, momentum and energy. In a

continuum regime the plasma modeling starts from the Boltzmann equation [55]


f+c.Vf /+.V,(f /m)=f (2.1)
at t' \ co

The index 'i' refers to single species. Here f is the distribution function, c the charged

particle velocity, F the force, V the gradient operator in physical space, Vc the gradient

operator in velocity space. The subscript 'coll' refers to the change of distribution due to

collisions. F represents the force vector which can be the Lorentz force on the charge

particles. A detailed discussion on the theory of Boltzmann equation and its relation to

the fluid models for radio frequency glow discharges is presented in [55]. The index 'i'

implies that (2.1) will be written for every kind of species. To enable a self-consistent

solution of the above, we need to solve for the electric field E by means of the Poisson

equation.

V.E = (2.2)
G0









Table 7 Kinetic model studies of discharge plasmas
Author and year Work description

Ku er [78, 193 Monte-Carlo (MC) simulation in rf
Kushner [78], 1983
discharges.
Particle-in-cell (PIC) simulation of parallel
Boswell et al.[79], 1988
plate RF discharge
Surendra et al.[80], 1991 PIC-MC technique for rf parallel plate glow
Birdsall et al.[67], 1991 discharge in helium
Direct numerical integration of the
Hitchon et al.[81], 1991 Boltzmann equation using BGK model for
the collision term
DiCo et al.[, 19 Direct integration of the Boltzmann equation
DiCarlo et al.[82], 1989
using Flux corrected transport.
Direct integration Boltzmann equation using
method of characteristics with a finite-
Riley et al.[77], 1994 difference treatment of the collision integral.
Solution is accelerated using a time-averaged
equation along with the full equation.


Another complexity in solving the Boltzmann equation (2.1) comes from the

evaluation of the collision term which is the source of the non-linearity. Typically, the

collision term is modeled using a particle technique such as Monte Carlo methods

(DSMC- direct simulation Monte Carlo). Other methods such as the particle-in-cell

(PIC) technique are also common. Kushner [78] uses a particle simulation of a radio

frequency parallel plate discharge using a Monte Carlo model. Boswell et al.[79] used

PIC techniques in a self-consistent simulation of parallel plate discharges. Birdsall [67]

and Surendra et al.[80] study the PIC-MCC (Monte Carlo collision) technique for

modeling discharges. As an improvement over the early PIC techniques to include binary

collsions (due to electrostatic effects), the Monte Carlo collision model was incorporated.

This is in a way, a hybrid of the PIC and the MCC techniques which existed

independently. PIC/Monte Carlo models of radio frequency discharges provide kinetic

information self-consistently by integrating the equations of motion of many particles









representing ion and electrons with a simultaneous solution of the Poisson equation with

a Monte Carlo treatment of the electron-neutral, ion-neutral collisions. The charged

particles are given an initial velocity and the position coordinates and the particle paths

are followed through time using the equation of motion until a quasi-steady state

condition is achieved. These simulations are then run for an extended period of time to

filter out statistical fluctuations and obtain smooth profiles. The Monte-Carlo algorithm is

employed to model the charge-neutral collision events following which the new position

and velocity distribution is updated for the next time instant.

The PIC [67] by itself is a method based on deterministic classical mechanics of

particles with prescribed force fields. In PIC, the simulation region is divided into a

number of cells and the resulting grid is used in the solution of a force field from which

the force on each particle inside the cell can be determined. The equation of motion of

each particle is integrated to get the new position and velocity. On the other hand, Monte-

Carlo simulation is a statistical technique which is basically probabilistic in nature and

the particle velocities are obtained by directly integrating the equations of motion. The

key steps of a integrated particle-in-cell and Monte-Carlo technique (PIC-MC) of

Nitschke et al.[68] are the following:

a. Particle position and velocities are updated each time step using the equation of

motion.

b. Particles which have crossed the system boundaries are accounted for.

c. Monte-Carlo algorithm to determine probable collisions and the resulting states

d. New particle positions are linearly interpolated to the nearest grid point to obtain

ion and electron spatial density profiles. Poisson's equation (Eq. (2.2)) is









calculated again to determine the electric field to be used in the force term in step

a.

Other solution methods for kinetic models include direct solution of Boltzmann

equation as in Hitchon et al.[81] where the equation is directly integrated numerically as a

partial differential equation. Here the collision term is approximated by a BGK type model

(apart from Monte Carlo collision models). Because of the strong convective nature of the

equation, higher order oscillation preventing schemes such as Flux Corrected Transport

(FCT) (DiCarlo et al.[82]) or the convectivee scheme" of Hitchon et al.[82] are employed.

The direct solution of the Boltzmann equation is slowed down by the very small time

resolution dictated by electron and collision relaxation times. This problem is overcome by

Riley et al.[77] using a hybrid formulation which employs a time-averaged description.

Here the Boltzmann equation is solved for only a few cycles at a time. The information

obtained from the Boltzmann equation solution is used to solve a set of time averaged

equations describing the heavy species. The averaged equations are solved for thousands of

RF cycles and then the Boltzmann equation is solved again and this process is repeated.

This is in a way another type of hybrid approach. Another technique along similar lines

where the faster varying scales are solved a few times and the information is used to model

the slowly varying processes is used by Hichton et al.[81]. This scale up technique

accelerates the solution process significantly. Other hybrid models will be presented next.

The various kinetic studies are summarized in Table 7.

3.1.2 Hybrid Modeling Approach

The computationally intensive kinetic simulations (as compared to fluid models) are

more accurate in the sense that they can used to model low-pressure (typically < 100mTorr

[68]) or highly non-equilibrium situations for which fluid models are highly suspect. The









hybrid models are typically employed when there is a need to preserve the accuracy and

capability of the kinetic schemes and at the same time reducing the computational burden.

Table 8 Summary of hybrid modeling studies

Author and year Hybrid model

'Bulk-beam' model where the bulk and high
Belenguer et al.[84] speed or beam electrons are treated as fluid and
particle species respectively.
Surendra et al.[80] 1990 'Bulk-beam' model

Fluid model with Monte-Carlo simulation for the
Fiala et al.[85] 1994 ion n
ionization term.


Hybrid models have been studied which combine aspects of both fluid as well

kinetic models and thereby achieve a balance between accuracy and efficiency. The charge

species in low-pressure discharges are generally not at local thermodynamic equilibrium

and one cannot assume in general that these species attain a thermodynamic temperature.

This lack of thermodynamic equilibrium complicates the modeling of discharges, since the

kinetic models determine the species velocity distribution function as against an assumed

distribution function in the case of fluid models. From another perspective, under low

pressures the thermal flux is comparable to the drift flux from the electric field and as a

result the isotropic behavior for macroscopic properties such as pressure is lost. The kinetic

models in such non-equilibrium cases give accurate results as compared to fluid models.

One type of hybrid model treats the fast or high energy electrons called as 'beam' electrons

as particles and the slow electrons or 'bulk' and ions as fluid and achieves a self-consistent

simulation of the entire discharge while addressing the non-equilibrium nature of the high

energy electrons. Such a model has been first introduced by Belenguer et al.[84] and used

by Surendra et al.[80]. Another type of hybrid procedure is used by Fiala et al.[85] where









the ionization source term is re-calculated by a Monte Carlo simulation after a few hundred

time-steps of the fluid equations.

3.1.3 Fluid Modeling Approach

Fluid models generally consist of a few moments of the Boltzmann equation for the

various species with an assumed distribution function. However, as discussed previously,

the validity of the fluid type description of the plasma is determined by the existence of a

local thermodynamic equilibrium. This issue is significant while modeling low-pressure

discharges. Nitschke et al.[68] compare the particle-in-cell model with a fluid model

simulation of a radio frequency discharge in helium. The results from the two models

were compared over a wide range of applied voltages and pressures. Considerable

agreement between the studies was observed for pressures greater than 100 mTorr

[68,55], which is much smaller than the atmospheric pressure regime of 760 Torr. As a

result most of the atmospheric pressure discharge simulations use fluid models which are

potentially attractive for computational studies if the limitations to their validity are

adequately addressed. In this dissertation, we are interested in the modeling of discharge

processes at atmospheric pressures.

Fluid models usually consist of the first three moments of the Boltzmann equation

which includes the species continuity, momentum and energy conservation equations.

The fluid species usually consist of the ions, equations and neutral species. The number

of species is typically dependent on the extent of plasma chemistry included in the model.

A typical fluid description of the model as used by Colella et al.[86] is presented below.

Species continuity equation:

)n (2.3)
n + V.(nu ) = R (2.3)









Species momentum equation:

(mnu) +V.(mznuu,) +V(n,kT)= q,n,EZ mm- nyu, -u,)+ m,u,R, + m,u,R
at Jn, +MJ |JRO <0
(2.4)

This is a generalized form of the momentum equation incorporating various phenomena.

The third term on the left hand side represents the partial pressure gradient of the ith

species. The right side includes all source terms- the Lorentz force, momentum loss from

collision and momentum losses from production and destruction, in that order.

Species energy equation:


+ ,, n -mu +e & + v. -mun i u, + su, + V.nku = q,n,.E -


I 2 n7a n-y, n u, .uI n *u^.u + n n M1 ,.u )+ 6, 6 (2 .5)
n w +mn 2- ~

+ imu,.u +, R,0 + I -nmJu.uJ + SJRj +z1S1J
i% ,2 J

kT
where s, Here the summation notation denotes the summation over the different
Y,-1

kind of species. In assuming a fluid model, each plasma component (index 'j') is assumed

to have a near Maxwellian distribution. The infinite set of moment equations is truncated

to include just the moments corresponding to the conservation of mass, momentum and

energy for each component. The indices 'i' and 'j' represents the particular species and

its complementary ones, respectively. The summation notation indicates summation over

all particle species which includes electrons, ions and neutral species.

These transport equations are solved in conjunction with the Maxwell's equations

for the electric potential which is written as










A== yZqn, andE=-Vq. (2.6)


The number of equations that are solved in the system is determined by the number of

species considered in the plasma chemistry and can result in a hugely coupled system.

However in reality, the number of equations actually solved is very manageable with the

various approximations to the model and the system of equations. One of the common

assumptions [20] is approximating the plasma chemistry to include just a single species

of positive ion and negative ion or electron. Also, in many cases, the neutral species are

not considered because of their lack of contribution to the net charge. However there

exists non-charged species such as metastables whose dynamics can differ from those of

the neutral species. Such species are usually considered when the complex plasma

chemistry is important for the modeling study. This is illustrated below in the context of

a helium atmospheric pressure glow discharge where a model with complex set of

equations, as given in Tochikubo et al.[63].

Table 9 Plasma chemistry considered in a helium discharge
e+He ->He+ +2e
e +He He(21S)+e
e +He He(23S)+e
He+ +2He He+ +He
He+ +2e -He* +e
He2 + 2e- 2He+e
He (2S)+ He 2He + hv
He (23 S) + 2He -* He + He
He (2S)+e He +e
He(23S)+He(23S)- He+ +He+e
Source : Tochikubo et al. [63]









The plasma chemistry shown in Table 9 needs to be modeled taking into account

the reaction rates to account for the species concentrations and transport. Predominantly,

modeling studies are one-dimensional (see Table 4), followed by two-dimensional studies

which are less frequent, while three dimensional studies are rare. It should be noted that

distinct higher dimensional effects have been observed in studies by Mangolini et

al.[62,56]. The fluid modeling studies have already been summarized in Table 4 and Table

5 Comparison of the results of different numerical models for the same set of modeling

parameters has been compiled by Surendra [87] for a 1-D RF discharge simulation. The

results indicate that fluid models, although less reliable than kinetic models, give

reasonably accurate predictions of discharge properties for sufficiently high pressure.

Also, the kinetic models, say PIC-MC, compare well with the hybrid models. The fluid

models are generally reported to show greater variation than kinetic/hybrid schemes.

Since we are interested in the discharge operating at atmospheric pressure, we will focus

on the study of fluid models for this dissertation.

3.2 Numerical Modeling of Glow Discharge Plasmas

In this section, we will present the detailed numerical modeling framework for the

atmospheric pressure glow discharge plasma. The plasma-fluid-Maxwell equations

presented in the previous chapter represent the generalized set of equations. In reality,

however, the set of equations solved are much different and less complex from the ones

presented, after the various problem specific modeling approximations applied to the

system. In this section, we will first discuss the various approximations that can be used

from which we can arrive at the final set of equations to be solved. Following this, the

numerical nature of these equations will be discussed to highlight the difficulties that will









be encountered during its solution and the techniques to handle them. Finally a solution

algorithm will be developed to handle the system of intricately coupled equations.

3.2.1 Plasma Discharge Governing Equations

The fluid model for the plasma can be described by the following equations which

describe the first three moments of the Boltzmann equation for the various species. They

are the continuity, momentum and energy equations. The plasma is considered as a multi-

component fluid comprised of two types of primary species, namely, ions and electrons.

A generalized description of the plasma model was presented in section 3.1.3 (Eqs. (2.3)-

(2.6)). The various modeling studies use a variety of approximations to simplify the

system of equations. In the following section, we will present the commonly used

versions and highlight the specific variants used in different studies.

We will describe the composition of the fluid by means of electrons and ions. By

ions we imply all charged species which have a mass comparable with the neutral

particles. As seen in the previous chapter, the uncharged metastable species also play a

significant role in the plasma chemistry. But we shall neglect them in favor of the present

focus towards developing an efficient modeling strategy for these plasma fluid equations.

Now, we will briefly develop the form of the fluid equations to be used. The continuity

equations for the positive ion and electron species are given by

ne +V.(nv ) = R, Rd (2.7)
St

,+ V.(nv) Rd, (2.8)
aSt
The subscripts e and i stand for electrons and positive ions, respectively. The

parameter Re is the ionization or creation term and Rd, the recombination or destruction

term. With no loss of generality, we can have a similar equation for negative ions. The









terms on the right hand sides of eq. (2.7)-(2.8) represents the production from ionization

and the loss from destruction, in that order. These two terms are obtained from the

collision terms in the actual Boltzmann equation(2.1). The exact forms of these terms are

obtained from the type closure model one uses to represent the particle phenomena at the

fluid model. This can be done in a variety of ways as mentioned in [55] such as assuming

a velocity distribution function (Maxwellian) or by assuming phenomenological

expressions for the collision terms. The general form of the source terms is given as

Re = nv
R c nv(2.9)
Rd = nd

where u, and Ud are the creation and destruction frequency which are function of known

parameters such number density ni,e, velocity ui,e etc. The species velocities are given by

ui and ue for the ions and electrons respectively. These velocities are obtained from the

momentum equations given below. The generalized momentum equation is given in(2.4).

It can be written simply as


- (nemeve) + V.(nemvev) = -VPe + nemqE R (2.10)


-(nmv) + V.(nmv ) = -Vp, + nmqE R (2.11)
at

where Rmi,me represents the collisional momentum loss of the ith species or electron and

Pi,e is the partial pressure which can otherwise be written as ni,ekTi,e. k is the Boltzmann

constant, Ti,e is the species temperature and mi is the species mass. This form of

expression for the pressure tensor assumes isotropic behavior and is strictly valid only for

a Maxwellian distribution. The momentum loss Rmi,me is usually written as

Rm = meneevWme (2.12)









Rm, = m nv.vm, (2.13)

where umi,me is the frequency of collisions which result in momentum transfer. The

second term on the right hand side represents the body force from the electric field with

qi,e being the species charge. Combining (2.7)-(2.8) and (2.9) we have

e + V.(nve) =V,ene d,e e (2.14)


an +V.(n,)= vc,,= -vd,z (2.15)
at

and substituting (2.12)-(2.13) into (2.10)-(2.11) we have

-(nmeve )+ V.(nemeveve) = -V(nekTe) + neE men eeVe (2.16)
at


a (nmv )+V.(nmvy,) = -V(nakT)+ n ,E- mnvv, (2.17)

Now we will simplify the momentum equation (2.16)-(2.17) using (2.14)-(2.15), to yield.

av
neme e +neme (v eV)ve + ne e (Ve d,e)+
at (2.18)
V (nekT -) nefe = -nemeveVme

av,
at (2.19)
V (nkT,) nqE = -n, mvym

or


1 vm v v e me (nek me
+ (Ve.V) + (Vce de= -ve (2.20)
me at Vme Vme tnemee meVme

1 av v v V(nkT) q E
v +(v,.V) + (v, -Vd,) =-v (2.21)
vm, t Vm, vm, n l m lV, m Vm,









where time derivative is a substantial derivative. Now we try to simplify the momentum

equation using a scaling argument. The first term in (2.20)-(2.21), which is the local

derivative scales, as the velocity times the ratio of the driving voltage frequency and the

collision frequency. This ratio is usually very small even for driving voltages in the

megahertz range as the collision frequency is of the order of a few gigahertz.


m1 ,e (2.22)
Vmn,me at Vm,me

Hence the local derivative can be neglected. For the convection term, momentum loss

collision frequency )mi,me scales as the ratio of the thermal velocity and the mean free

path.


(v, -V) V O v, (2.23)
Vme,m L vT


As long as the thermal velocity is comparable to the species drift velocity ( v, <1) and
VT


we are in the continuum regime (- Kn <1), the convection term will be very small
L

compared to the right hand side (force and collision loss terms). However, for the heavier

ions the thermal velocity is comparable to the drift velocity at very low pressures and

hence the inertial terms become significant. Hence, at low pressures the distribution

function departs significantly from the Maxwellian or equilibrium distribution function.

This makes the diffusivity and pressure tensors anisotropic. But at higher pressures, the

above simplification is not too costly and the ion momentum equation has the same form

as the electron momentum equation. The third term on the left hand side of (2.20)-(2.21)

is the ratio of the creation or destruction frequency and the momentum loss collision









frequency. It is usually expected that the collision processes are more frequent when

compared to one which causes ionization and the ratio is negligible. The final simplified

form of the momentum equation can be written as

V (nkT, ) qE
V k= -ve. (2.24)
emeVme meVme

V(tnkT) qE
= -v (2.25)
tnm Vm, m Vm1

Now let us define the species mobility pi,e and the diffusivity Di,e as


, q (2.26)
l,eVmi,me



m,e mi,me

The final form of the momentum equation becomes

V(nD) E =-ve (2.28)
e


u-D ,E = -v, (2.29)


This form of the momentum equation is called as the drift-diffusion approximation.

Further variants of this form have been used in many studies such as by Boeuf et al.[31]

etc., and is shown below.

DeV ) E = -ve (2.30)
en

V(n uE=-v (2.31)
n,









Now we will briefly derive the simplified form of the energy equation. The general

energy equation can be written as


n + V. (n, e + vP + Qe) nqE.v, = Re, (2.32)
at

where s is the particle energy, Pe the electron pressure, Qe the heat conduction and Ren the

energy loss from collisions etc. As before, assuming a Maxwellian distribution for the

velocity, we have

2
Pe = nkT = ncE, (2.33)

and


Qe = -KeVT =-KeV (2.34)
3 k K }

Here ST is the thermal energy component. The thermal conductivity Ke can be written in

terms of the diffusion coefficient De, from kinetic theory results as

3neDek
Ke (2.35)
2

The fact that all the parameters such as conductivity, diffusivity are written as scalars is

because we assume isotropy. Substituting(2.33),(2.34) and (2.35) into (2.32) we get the

complete energy equation as

+V+V. ne 1 +-n v+-neDV(Ere) -n .ve=Rn (2.36)
at 3

The form of the energy loss rate Ren as well as the parameters Rmi,me, R1, Rd, the mobility

and diffusivity etc, needs to be identified to close the system. These parameters are

usually related to some measurable macroscopic properties. Boeuf et al.[31] and

Gogolides et al.[55] among others assume that the above energy loss rates etc., depend on









the mean electron energy as they do under equilibrium conditions when the electron

energy gain is locally balanced by the losses. If we assume the following form for Ren

[31, 55]

Re = K (E / N) Nn, = VLn (2.37)

then the energy equation under equilibrium conditions yields


KL (2.38)
N

At equilibrium, all the mean electron properties such as se, De and -es etc. can be assumed

to be a function of E/N. This dependence of the mean electron properties on the local

electric field as against the velocity distribution function (in the Boltzmann equation) is

called as the local field approximation. This approximation rules out any non-local

effects for these properties. The values of the electron properties are the same as that

obtained in a uniform electric field condition, where the collision term balances the force

term. Detailed discussions are presented in [66,64,72,55] Hence, from (2.38), Ren is also

a function of E/N. A similar form is used for other parameters such as Rmi,me, Re, Rd. The

equilibrium dependencies of the various parameters on the ratio E/N are tabulated by

solving a Boltzmann equation or a Monte Carlo code under the appropriate conditions. It

should be noted that in some hybrid approaches (e.g. Kushner and co-workers [88-89]),

the fluid energy equation is dropped and the kinetic energy equation is solved using a 2-D

Monte Carlo method.

The ion energy seems to have little effect on the ionization and attachment rates

and hence is not solved unless required. To summarize, the electron inertial terms are

insignificant considering their negligible mass. Hence, the momentum balance reduces to

the drift diffusion form. For the ions, the inertial term is significant at low pressures (<








100 mTorr [55]) and is included in some studies. At high pressures, however, a drift-

diffusion approximation is considered reasonable [31,55]. A summary of conservation

equations for the mass, momentum and energy used to model high frequency glow

discharges is given below. A similar set of equations is presented by Gogolides [55].

Continuity equation:

Electrons: Lne + V.(nev) = K', (E N)Nn, Kd, (E / N) Nn (2.39)
at

Ions: nL + V.(nv) = K, (E/N)Nn -Kd,i (E/IN) Nn (2.40)
at

Momentum equation:

Electrons: ne/eE -V(nD)= nve (2.41)

Ions: n E V(nD) = nv (for high pressure discharges) (2.42)


n -n, + +n n, (v.V)v, +nm,v,N(K, -Kd,,)
a t (2.43)
+V (nkT) -nqE = -nmvNK,,

(for low pressure discharges)

Energy equation:

e) +v. Ve + nee eD V(s )V
Electrons: at 3 e e ) (2.44)
-neqE.v, = KL (E / N) N

Ions: No energy equation is considered as ion energies do not affect ionization processes.

To summarize, the various assumptions made in the development of the fluid equations

are:

1. the pressure tensor is isotropic

2. the drift energy is comparable to the thermal energy









3. heat flux is proportional to the electron temperature gradient

4. the mean electron-neutral collision rates depend only on the electron mean

energy.

Most of the numerical modeling studies of high frequency glow discharges use whole or

part of the above system of equations along with the Poisson equation for the electric

field which is written as.


V.E = (~4 ,- e) (2.45)


The values of the parameters K, mobility, diffusivity etc. are either obtained by solving

the spatially homogeneous Boltzmann equation, or assigned curve fitted expressions from

experimental data or a combination of both. A standard model for the ionization

coefficient is the Townsend expression presented in [61] which is shown below.

16
04
,uEp
K /= (34)e (2.46)


Here p is the pressure. Some of the numerical issues associated with the system of

equations discussed above are as follows. The set of moment equations described above

are intricately coupled and strongly non-linear with a wide range of spatial and temporal

scales. The solution profiles for the charge particle densities contain steep unsteady

gradients especially in the sheaths. The resolution of these sharp gradients accurately

requires well refined grids leading to time-intensive computations. Sharply varying

production terms such as ionization processes result in stiff source terms which also limit

the time step. For these reasons, 3-D and to a lesser extent 2-D modeling studies are less

frequently attempted. The continuity equation for the ions and electrons are essentially

inhomogeneous Euler equations (as diffusion is mostly small compared to drift) of gas









dynamics augmented with the electrostatic force, collisional and source/sink terms.

Solution of this system poses challenges due to its multi-time-scale nature which we will

illustrate now. The critical time scales in the problem are listed below.

1. The discharge operating frequency

2. The ion drift

3. The non-linear source term- ionization, collision etc.

4. The electron drift

5. The dielectric relaxation time

The discharge operating frequency time scale is given by


r =1 (2.47)


where co is the driving frequency. In order to illustrate the time scales we will use the 1-

D modeling result from Lee et al.[34] which is shown in Figure 3. The driving voltage is

1.5KV with a frequency of 10 KHz. It can be seen from the figure that the appropriate

length scale for the problem would be the size of some physical length such as the

positive column or the cathode fall (Figure 3) which is 0.05 cm in a 0.3 cm gap. The

drift velocity from the drift-diffusion approximation is yE The values of the mobility

for the helium ions and electrons can be obtained from literature and the electric field

value is the peak value in the positive column (- 20 KV/cm). Then the drift time scales

for the ions and electrons can be calculated as the ratio of the width of the positive

column and the drift velocities.


dT ldr (2.48)
,dr nE
pro










rdr,elec dr (2.49)
elecE

To evaluate the source term time scale, we will use the Townsend formula for

calculating the ionization coefficient. From (2.9) the ionization term Re can be written in

terms of the ionization frequency v,. From (2.37) and (2.46) we can see that v, is of the

form

16

vC = peEp(34)ep = (2.50)
c

which gives the ionization time scale Tc as the reciprocal of the frequency. We use the

above form of the ionization frequency in calculating the timescale in Table 10. Now we

shall describe the origin of the dielectric relaxation time scale and evaluate it. This

timescale represents the relaxation of the electric field in the medium to adjust to the

change in the charge density and is usually the most stringent of all the time-step

limitations.

Taking the time derivative of the Poisson equation(2.45), we get


(V.E) = ne. (2.51)
at so at at

Now substitute for right hand side from the continuity equations (2.39)-(2.40), which

gives


V.E -V.(nev -n) (2.52)
at so

Getting rid of the divergence operator, we get

= e (ne- ne )+C (2.53)
at so









Substituting for the velocities using the drift diffusion approximation (2.28) yields

OE e
- -(enlE -, n E +D Vn, -DVn, )+C (2.54)
at so

Eq. (2.54) can be viewed as the evolution of the electric field due to the space charge and

leads to a time scale as


=el (2.55)
e (l ne + u, n) 1)

This restriction is highly severe and is usually overcome by a semi-implicit solution

of (2.54). In other words the above restriction can be viewed as the CFL condition for the

current continuity equation. Using the values of the various parameters from the 1-D

model of Lee et al.[34], the magnitudes of the various time scales are calculated and

tabulated in Table 10. The dielectric relaxation time becomes real small as the number

density increases which is the case for high frequency (MHz) discharges and hence needs

to be treated accordingly. As can be seen from Table 10, time scales involved in the

problem vary over a wide range and this needs to be taken into account while developing

a numerical algorithm. If we take into account adequate grid resolution in the sheath

regions, this time scale will be further reduced.

3.2.2 Numerical Algorithm

Modeling of plasma fluid equations has been extensively attempted in the last

decade or so predominantly for 1-D applications and occasionally 2-D applications.

Emphasis was placed on plasma chemistry [89] and studying the limitations of the model

by comparison with experiments or kinetic theory based models [68,31,67,80]. Efforts to

develop efficient solution algorithms for such systems such as Boeuf et al.[31] Hammond

et al.[90] Golubovskii et al.[56] Colella et al.[86] etc. will be discussed in the following









section. These are representative studies and a more comprehensive review has been

presented in the previous chapter 2. Except for Golubovskii et al. where an atmospheric

pressure glow discharge between two parallel plates is studied, all other efforts are

devoted to low pressure high frequency discharges. The set of governing equations

employed by the various efforts differ from each and are summarized in Table 11. While

the equations have the same form as given in (2.39)-(2.44) the methods of evaluating

parameters such as the mobility, ionization etc. may differ form each other. It may be

recalled that the two forms of the momentum equations ((2.42) and (2.43)) which include

a fully approximated drift-diffusion form and the full transport equation. Both these

forms are used in the modeling studies, but the simpler drift-diffusion form is preferred

and valid at high pressures for ions. It is generally always used for electrons because of

its light mass causing negligible inertial effects.

Table 10 Representative time scales of various processes in a 1-D helium discharge
model operating at 10KHz, 1.5KV rms, calculated based on values of
parameters.
Time scale Value (s)
Operating frequency timescale, r, (2.47) le-4
Ion drift rdron (2.48) 2.5e-7
Electron drift drelectron (2.49) 2.5e-9
Dielectric relaxation Tel (2.55) 2e-9
Ionization timescale rc (2.50) 1.4e-10
Note: Calculate based on parameters taken from Lee et al.[34]

It can also be seen from Table 11 that the electron energy equation is not solved

for high pressure simulations (Lee et al.[34], Massines et al.[6]) when the local field

approximation is valid. Not solving for the energy equation, however, prevents us from

understanding the energetic such as collisional heating or Joule heating. While, the

choice of governing equations and their approximations vary from one study to the other,










the numerical characteristics remain fundamentally the same. The governing equations

are essentially the inhomogeneous Euler equations of gas dynamics augmented with

electrostatic forces, collisional and source/sink terms. The system of equations represents

multiple time scales (Table 10) and widely separated spatial scales. The spatial scales can

be adequately resolved by using local grid refinement techniques [116] in the regions of

sharp gradients coupled with higher order discretization methods which have been

successfully applied for modeling neutral fluid flows.

Table 11 Summary of governing equations used in different studies
Governing equations used
Author, year Topic Continuity Momentum Energy
Ion electron ion electron ion Electron


Boeuf et
al.[31] 1995




Colella et
al.[86] 1999





Hammond et
al.[90] 2002



Lee et al.[34]
2004
Golubovskii et
al.[56] 2004


Roy et al.
[15,115] 2005


2-D rf discharge
in argon at
100mTorr (low
pressure)
Develops a
conservative
finite difference
method for
plasma fluid
equations for
high density low
pressure plasmas
Numerical
modeling of
low-pressure,
high density rf
discharges in
Helium
1-D parallel
plate discharge
in helium at
atmospheric
pressure
2-D asymmetric
discharge
actuator
modeling in
air/He (300 torr)


(2.40) (2.39) (2.42) (2.41)


(2.40) (2.39) (2.43)






(2.40) (2.39) (2.43)





(2.40) (2.39) (2.42)




(2.40) (2.39) (2.42)


(2.44)


(2.41)


(2.44)


(2.41)


(2.44)


(2.41)




(2.41)









It should be kept in mind that the choice of convection scheme is also determined

by its performance in the range of Peclet numbers encountered and also in the presence of

a source term [91]. In the following discussion we will restrict ourselves to the set of

governing equations as presented in (2.39)-(2.44).The challenge here is to efficiently

handle multiple time scales, and the treatment of the reaction source terms. Such

complications are more common in the modeling of reacting flows where detailed

modeling of chemical kinetic mechanisms can result in a stiff system of governing

equations. The focus of the present dissertation includes the development of an efficient

solution algorithm to handle the stiff multi-scale problem as seen in the context of the

plasma dynamics. The discharge equations as described in (2.39)-(2.44) involve wavelike

ionization fronts which travel through the domain before charging the dielectric or any

other boundary surface. The solution procedure needs to be designed to handle the

multiple scales efficiently. However, one distinct difference between chemically reacting

flows which pose similar challenges is the electric field effect. It has been shown in [86]

that the plasma sound speed is modified by the electric field drift with the effective ion

sound speed given by


c = (2.56)
m

Here Te refers to the electron temperature which is usually very high compared to the ion

temperature (which is closer to the neutral particle temperature). Thus the modified sound

speed is much higher than the neutral sound speed (as ion temperature is neutral). Note

that the first term in the square root with Te comes from the approximation of the electric

field in the electrostatic pressure. In the absence of a driving electric field, this term

vanishes giving the neutral fluid sound speed (the second component in the square root).









When Te >> Ti -T, the electric field dominates the thermodynamic pressure. Similarly,

the characteristic speed for electrons can be obtained as


ce =k .(2.57)
me

The above expression for the ions shows that the electric drift component is more

significant than the thermal component. The ratio of the two species sound speeds goes as

the inverse of the square root of their respective mass. The electrons being lighter than

ions are faster. In general, from the range of electric fields considered in this problem, the

ion and electron sound speeds are much higher than that of the neutral flow and hence the

bulk of the acoustic effects are caused by the driving electric field and not the low speed

neutral flow. This is different from normal chemically reacting flows where the species

characteristic speeds are comparable to that of the neutral fluid. This distinction needs to

be kept in mind while dealing with the above strongly convective equations where the

reaction and convection time scales are comparable.

Integrating such a system comprising many time scales gets very demanding in

terms of computational resources. Both implicit and explicit techniques have been used to

march in time with the choice usually based on the extent of stiffness of the system,

expected time accuracy etc. Implicit integration methods are attractive in the sense that

they are unconditionally stable, but are not time-accurate and in some non-linear

situations are even computationally burdensome. This is necessitated by the need to find

the root to a highly non-linear equation which might entail very small time-steps to

achieve convergence. Also, there is a burden of huge computational requirement. This

method works reasonably well when the overall timescales are not very far apart. Several

studies aimed at improving the efficiency of stiff integration schemes have been studied









in the context of modeling reacting flow and atmospheric chemistry which involve taking

advantage of the problem specific approximations to the Jacobian. In other examples, the

stiff variation of the reaction terms which is typically exponential makes it advantageous

to study the variation of the logarithmic quantities as against original variables [90].

Another alternative to the implicit procedures used above is the use of matrix free

Newton-Krylov methods, which have been reviewed in [92].

Both sequential and fully implicit approaches have been used for discharge

modeling. Roy et al.[15-16] used a globally implicit finite element procedure where the

system of species continuity and momentum equations are assembled as part of a global

matrix to solve for the solution vector. Sequential approaches have also been employed

such as in the study of Hammond et al.[90]. Hammond et al.[90] uses a hybrid

formulation to solve the plasma fluid equations, where the ions species being slower are

treated by an explicit 4th order Runge-Kutta method followed by an implicit treatment of

the fast electron equations. The implicit Euler (first order) and implicit Runge-Kutta

(second order) methods were used with Newton-Raphson iterations to overcome the non-

linearity. Further, the resulting Jacobian is simplified by neglecting the weak off diagonal

contributions reducing to a block tri-diagonal system. However, taking into account the

computational overhead from the iterative procedure for the fully implicit treatment and

the cost of inverting the matrix more than compensated for the gain over explicit

procedures. Implicit procedures though are used in the context of obtaining steady state

solution with stability and not time accuracy as the consideration. Acceleration to steady

state is usually achieved using preconditioning techniques (for e.g. [93]) which decrease

the condition number of the system with widely separated scales. These techniques have









been extended to time accurate computations by Merkle [94]. Schwer et a1.[95] uses a

consistent splitting of operators with implicit transport computations of reacting flows

using local preconditioning (Merkle et al.[94]).

Recent modeling strategies use semi-implicit hybrid procedures which seem to be

more attractive especially in the context of operator split methods. Here the stiff portions

of the system are handled using efficient solvers and integrated with the non-stiff solution

procedure. A common practice to handle the stiff reaction terms is to integrate the

reaction ODE using a stiff ODE solver (CVODE, DVODE, LSODE) over a time step

determined by a non-stiff process. These are essentially based on the backward

differentiation formulae of Gear. This integration time-scale in most problems is the

convection time-scale of species which is considerably larger. Such a procedure makes

for a sequential approach and thereby enables handling smaller time scales by allowing

integration sub-cycles for the relevant portion. However, stability and coupling demands

appropriate corrector steps and handling of splitting and conservation errors.

Source term treatment in plasma simulations has been discussed by Hagelaar et

al.[96] where a linearized implicit correction is used to alleviate the stiffness by means of

a semi-implicit solution. This is similar to another common technique which is the point-

implicit treatment where the source terms are treated implicitly with a single Newton

iteration. Here the source term is linearized and written in terms of a Jcobian which needs

to be evaluated every time step. Also, if there are a large number of species with

complicated chemistry, inverting the sparse Jacobian becomes expensive and it is

common practice to resort to simplifications to accelerate convergence. Updating of the

Jacobian is also done only when absolutely necessary and is not usually done for every









step unless it seriously affects convergence. It is shown in [97] that this point implicit

method is a special case of a generalized time scale preconditioner where the stiff

processes are scaled in time to the speed of slower processes. Preconditioning techniques

used to handle multiple scales to solve steady state problems and are reviewed by Turkel

[93]. Implicit treatment of stiff terms provides solutions which are accurate at the slow

scales and stable at the fast scales. Overall, a fully implicit technique is rarely used,

especially with the strong coupling between the various species equations and the

resulting non-linearity. However, with the use of efficient solution techniques such as

implicit multigrid methods and efficient solvers like LU-SGS of Yoon et al.[98] high

speed reacting flows with strong convection have been successfully handled implicitly. In

spite of using inversion techniques as above, the species chemistry Jacobian makes it a

block inversion which can get expensive. Hence, it is common practice to simplify the

species chemistry either adaptively or uniformly or to treat the reaction terms in a

partially implicit fashion so that block inversion can be reduced to a scalar inversion.

Explicit methods have been developed to deal with stiff problems, such as

stabilized multi-stage Rung-kutta schemes, Runge-Kutta-Chebychev (RKC) [99] schemes

which possess extended real stability intervals. However, with the stability range

dependent on the number of stages employed, the computation increases for handling

very stiff problems and is almost comparable to implicit schemes. A review of explicit

schemes for stiff problem can be found in Verwer [99]. Explicit time-scale splitting using

computational singular perturbation (CSP) is presented by Valorani et al.[100]. Here, the

fast and slow time scales are split such that the algorithm marches with the slow time









scale with the faster time scales being taken into account at the end of the integration

time-step as a correction.

Operator splitting methods are used in reacting flow computations and are

discussed in [101]. The treatment of source term using operator splitting has been popular

in the computation of reacting flows and discussed in Najm et al.[102] and Oran and

Boris [101,102]. Here, the stiff reaction operator is treated as a stiff ODE using the

DVODE integrator [103] and is embedded as part of a projection method. Colella et

al.[86] employs a time-split integration algorithm of the species equations in a sequential

form to decouple the three possible stiff processes namely ionization, convection and

dielectric relaxation. The different species equations are handled in a sequential

formulation in a predictor-corrector approach wherein the evolution of the various

species is appropriately coupled to the electric field.

In the present multiple study, we have employed an operator-split sequential

solution algorithm. Some examples of studies which employed such ideas are Najm et al.

[102] in combustion problems, Verwer et al. [99] in atmospheric chemistry modeling and

Tyson et al. [104] in chemotaxis models. Here, the solution procedure can be adapted to

handle the individual processes efficiently and realize overall gain in computation.

However, such an approach requires careful attention to stability considerations and

performance is highly dependent on the physics of the problem. In certain situations the

presence of competing stiff processes can lead to a system with much larger dynamical

time-scale than dictated by the individual processes. Such dynamic equilibrium usually

cannot be predetermined in general for a non-linear advection-diffusion-reaction system.

The stiffness of the reaction part is typically overcome by using stiff integration









procedures in ODE integration packages [105] such as the ones based on the backward-

difference formulae (BDF). In using time-split algorithms for processes operating in a

range of time-scales, the choice of time-step size is typically determined by the smallest

time scale, but need not necessarily be chosen as such. To speed-up the solution

procedure, an intermediate time scale is chosen to advance the overall system in time,

while the faster processes are advanced by sub-cycling within the time-step. In the

present study, the time-step dictated by the slower ion species convection is targeted to

march the full discharge system while sub-cycling is used for the faster processes. Also, a

predictor-corrector approach is employed to ensure sufficient coupling between the

electric field and the species densities. A strong coupling is essential for achieving stable

time-accurate simulations while using a sufficiently large global time-step. This method

is integrated with a multi-block finite-volume algorithm capable of handling 3-D

curvilinear-grids [19]. The implemented solution algorithm conserves space charge and

avoids dielectric relaxation time step restrictions. In the following section, an outline of

the algorithm is presented.

The split solution algorithm consists of the following steps:

1. Predictor step

2. Solving the Poisson equation

3. Corrector step

1. Predictor step: At the beginning of the n+lth time level and for the kth species, we

have n, v,,E" and vt as the known quantities. Here the species continuity equations

(2.39)-(2.40), along with the drift-diffusion momentum equations (2.41)-(2.42) are

integrated using lagged values for the various coefficients (as they are a function of the









electric field E). The source term is integrated using a higher-order (4th-order) BDF using

the CVODE solver [103]. The convection and diffusion operators can be treated either

implicitly or explicitly. In this case, we will employ a second-order upwind for the

convection term and second-order central difference for the diffusion term. The

continuity equation can be written as


+ V.(n V )= c,k k (2.58)


As can be seen from Table 10, the presence of processes with disparate time scales

can possibly be better handled using operator splitting. Three types of splitting are

popular, namely,

(i) First-order splitting

(ii) Strang splitting

(iii) Source splitting.

(i) First-order splitting: The first-order splitting [101] can be written symbolically as

7k = TMk (At')S(At) n (2.59)

where S is the reaction operator integrated using the ODE solver CVODE [103] and T is

the transport operator. MIis the number of substeps used for the transport term

integration to march to the global timestep At. Therefore, we have

At
At' = (2.60)
M k

(ii) Strang splitting: The Strang's time splitting method [106] solves for the conservation

law without the source term and the ordinary differential equation (ODE) with the

unsteady term and the source term in an alternating fashion. This method is formally

second-order in time which is achieved by splitting the operators symmetrically. In this









case, the transport term integration is usually split into two halves to achieve the

symmetry since the ODE solver used for the reaction part is computationally

burdensome.

nk TAk (At') S(At) T^ (At')nn (2.61)

It is worth noting that, in both the first-order and Strang splitting procedures the initial

guess for the reaction part is not directly from the previous time-step, but after a half or

full time-step of the transport term integration. This results in the introduction of stiff

transients in the solution which are nothing but an artifact of the splitting errors arising

from the lagged treatment of the coupled processes. These can have significant impact in

the presence of strong non-linearities. Even though the Strang splitting is formally

second-order, it is rarely achieved [107] in certain stiff problems, where it is known to

deteriorate to lower order accuracy. To overcome the solution discontinuities which gave

rise to stiff transients in the above two splitting methods, the source splitting treats the

transport as a piecewise constant source.

(iii) Source splitting: In this study, we limit the discussion to first-order source splitting

[99]. For the source splitting, we can write

k = S(At)n"n (2.62)

and

nk nn
Ihk k (2.63).
At

Therefore, we have

Sk= (T(At')+ k At')) n" (2.64)









where S is the reaction operator integrated using the ODE solver CVODE [105] and T is

the transport operator. MIis the number of substeps used for the transport term

integration to march to the global timestep At. Thus, we have

At
At' = (2.65)


The ODE solver employs the following:

A 5th-order BDF for time integration.
Newton iteration for non-linearity.
A direct method with a banded treatment of the Jacobian.
Normal mode with subcycling within the time-step.
Relative and absolute tolerances of le-12 and le-14 respectively.

The above strict tolerances were chosen so that the ODE integration is almost exact.

After the predictor-step, we have obtained the guessed values for the species number

densities based on the lagged value for the electric field. The strong super-linear

dependence of the convection on the electric field dictates the overall stability and

requires strong coupling between the variables.

2. Solving the Poisson equation: Here we solve a semi-implicit version of the Poisson

equation to overcome the space charge stability constraint leading to the dielectric

relaxation time-step constraint. This time step restriction is usually one of the most severe

and is equivalent to the CFL type stability criterion for the current continuity equation.

The Poisson equation is


V.E = (n ) (2.66)
so

The origin of the restriction has been illustrated earlier in this section in (2.51)-(2.55).

This restriction is primarily caused by the non-implicit treatment of the electric field in








the species momentum equation. In order to overcome this, a linearized implicit treatment

for the species number densities is used.


V.E" =e(n -n) -n At 8 on' (2.67)
so So at at)

which can be written using the species continuity equations as

V.E-1 (e< -At(V.(n,,) V.(n e))" (2.68)


Note that the only term in the right hand side which contains the velocity and hence the

electric field (through the drift-diffusion form of the momentum equation) is the one

inside the brackets withAt. By treating this term implicitly we can overcome the

dielectric relaxation time step restriction. The implicit version of the equation can be

written as

V.E" e (" At(V. (n,) V.(nev ))

e n -n-At V n v .') -nv (2.69)
so
= (- n" At (V. ,E" V(Dq )) V.( eE.1 V(Dee)m))


Here the ion density at the time level n+1 is approximated by the predicted density at the

end of step 1. The species velocities are directly substituted using the drift-diffusion

approximation so that the electric field can be treated implicitly. In the event, the full

species momentum equations are solved in step 1, the predicted velocity thus obtained

(with the bar) will be used in the Poisson equation. In that case, the non-linearity will be

difficult to overcome and will need a Newton-Raphson type treatment. The equation

(2.69) leads to a symmetric system for the electric field which needs to be inverted. This









can be done using a point Gauss Seidel method with successive over relaxation (GS-

SOR) or a line SOR. The cost of this step is possibly one of the most demanding and

hence needs to be performed as less frequently as possible. The choice of the global time-

step At will be determined taking this into account.

3. Corrector step: At the end of the previous step, the predicted species densities 7 and

the electric field at the new time level En+1 are available. Now we will obtain the

corrected densities at the new time level. It is worth noting here that the corrector step is

needed to ensure adequate coupling between the electric field and the species number

densities. In other words, we require appropriate coupling between the species continuity

equation and the momentum equation (or the drift-diffusion equation). Hence, the

corrector step is the same as the predictor step, but performed with the updated

coefficients using En+. It is also worth noting that the if the diffusion term and the source

terms are not strongly impacted by the electric field, then the corresponding operators can

be removed from the corrector step as correction is required only for the convection term

in such cases.

Stability. The stability criterion for the above predictor-corrector approach with stiff sub-

steps is not straightforward. Of the time scales listed in Table 10, the ionization source

term and the dielectric relaxation time step restrictions are not binding in the above

algorithm. The implicit treatment for the reaction source term makes it unconditionally

stable and is limited only by the non-linearity. The predictor-corrector formulation with

the semi-implicit solution of the Poisson equation ensures space charge stability

overcoming the dielectric relaxation time restriction due to the non-linearity. With the

split integration procedures for the convection and other processes, the stability of the









system will be determined by the CFL conditions for the individual steps. The global

integration time-step At, is determined by the slower moving ion CFL condition. Here X

is the CFL number for the different species.

AF = vma (2.70)
Ax

In the case of a number of different species considered, the choice of the global time-step

is usually the convection time scales of the slowest species, unless there is a huge

variation. This is important if we are considering the neutral species in reaction chemistry

whose dynamics are determined by the slower moving neutral fluid. Once this is fixed,

the time-step for the faster electron species and the corresponding number of sub-steps

M, can be obtained, so as to satisfy the electron CFL limit. However, studies [108]

indicate that the CFL stability limits of the split algorithm deviates from that of unsplit

algorithm depending on the number of sub-steps M, employed. Knio et al.[108] shows

that the critical CFL number due to electron convection sub-stepping decreases

monotonically as the number of sub-steps increases in two-dimensional modeling studies

of reacting flows. Also, one-dimensional studies indicate a limiting value of the critical

CFL being achieved as the number of sub-steps is increased. In spite of the stability

criterion becoming stringent the overall computational savings is generally substantial.

3.2.3 1-D Parallel Plate Discharge Modeling in Helium

We study the following test problem to validate our computational capability. A

1-D parallel plate (5mm gap between electrodes) discharge at atmospheric pressure using

just two-species chemistry, driven by a voltage of 1.5 KV at 10 KHz was modeled as

shown in Figure 8. This model problem has been studied by Ben Gadri et al.[20] and

Massines et al.[6]. The governing equations are the species continuity equations ((2.39)-









(2.40)) and the momentum equations with the drift-diffusion approximation ((2.41)-

(2.42)). No energy equation is solved for this model. The gas used is helium with a small

impurity of Argon (0.5 percent) to simulate the Penning effect arising from the

impurities. The parallel plates used are covered by a layer of dielectric with permittivity

9.0 and 0.6 mm thick. A detailed description of the arrangement is presented in [20]. The

various transport parameters, such as the species mobility, diffusivity and ionization

coefficient are tabulated as a function of the electric field using the 'BOLSIG' database.

A simplified kinetic model with various collision cross sections is used to develop the

database for the various types of species and these are locally interpolated in the code.

The boundary conditions are such that the fluxes near the boundary charge the

dielectric which is very significant in terms of maintaining a controlled discharge.

However there are inherent difficulties in the fluid model of the plasma in representing

the near wall sheath regions. Hence different sets of boundary conditions have been

employed in different studies [34,58,62]. For the 1-D species continuity equation we need

two boundary conditions in space for the number density. A few of the different choices

of boundary conditions employed are shown below in Table 12. A detailed discussion of

the issues in the choice of boundary conditions for fluid models of plasma is presented by

Hagelaar et al. [109]. In the present discussion we employ set of boundary conditions 1

from Table 12. A uniform time-step of 10-9 seconds is used for the time-integration and

the domain is uniformly discretized into 201-points. A fine grid with double the number

of points is also used and the results are compared with that of Massines et al.[6]. The

discharge current over a period of one full cycle after a quasi-steady state has been

attained during the discharge is self-sustained and is periodic as shown in Figure 9.







76


Table 12 Different types of boundary conditions for the ions and electrons used in
modeling studies.
Ion BC's Electron BC's
Serial no
Cathode Anode Cathode Anode

On Sec On
1 = 0 N=0 emission -= 0
Ox Flux Ox
On On Sec On
2 0 0 emission -= 0
Ox Ox Flux ax
nv
nvt nvh th +SCe nvh
3 uEn + th4 ec En
4 4 4
emission
Note: Vth is the thermal flux at the boundary.


V,,e=1.5KV ( 10KHz


Dielectric with
dielectric constant
9.0(s/ so) and
thickness 0.6
mm(Li,Lu)


U up c Ic-r d


Figure 8 Parallel plate arrangement for modeling DBD discharge

150
present-201
100
S- Massines et al
.1 ...... present-401


50
i' o? -" -- V ^ --
i


Gap=3.8 mm(L)


time microsecondss)
Figure 9 Evolution of discharge current over one period of the discharge cycle










The curves compare the coarse and fine grid results obtained in the present model

with that of Massines et al.[6]. As can be seen from Figure 9, the current peak

representing the instant after the breakdown slightly differs from that obtained in the

reference study. Also, there is a slight difference in the magnitude of the peak current.

The possible source of difference could be traced to the slight differences in the choice of

material data for the various species employed. However, the overall physics is

adequately represented although quantitative differences exist.




*/

_^ /










U 0 U U 4U U t 9U It
time ( microsecond3)

Figure 10 Evolution of the gap voltage over one full cycle

The gap voltage which is the net voltage between the two dielectric layers (or the

air gap to be specific) after accounting for the charge build up on the dielectric surfaces is

shown in Figure 10. The deviation from the reference results is again seen, but the coarse

and fine grid solutions are not much different, indicating possible grid convergence. The

peak gap voltage is at the instant just before the breakdown and it drops sharply after the

discharge peak, with large currents as seen in the previous figure reducing the voltage.

The distribution of the two species especially, the ions and the electrons, during the
The distribution of the two species especially, the ions and the electrons, during the










discharge are of interest, especially for the present study where the force acting on the

fluid needs to be accurately calculated. The species number density distribution over the

domain at the instant of the peak current density is presented in Figure 11. The

comparison shows that the present calculations have a ten percent error in the peak value

of the ion and electron densities. The fine grid results show identical peak densities and

hence not shown here. The very small time step required in the above study makes

operator split treatment of the various stiff processes necessary. This split approach will

be even more valuable when the plasma dynamics need to be solved coupled with the

fluid dynamics. As the faster fluid dynamics evolves, they will interact significantly with

the evolution of the discharge.

x 10"
--- electron -present
ion- present
Sion-Massines
electron- Massines
B 4

3

S2





Axial direction (0-0.5 cm)


Figure 11 Comparison of peak current species density distributions of the present effort
with Massines et al.[6].

3.3 Integration of Plasma Dynamics with Fluid Dynamics

The above discussion pertains to the modeling of the high frequency glow

discharge independently which needs to be incorporated in the Navier-Stokes calculations

for modeling the coupled system. The first issue is the effect of the neutral flow on the

plasma dynamics which needs to be captured. The actual species momentum equation in









a neutral flow would need to include the pressure gradient contribution from the fluid

dynamics. We can argue the total pressure as a sum of the partial pressures of the plasma

species and the neutral fluid particles. The momentum equation (2.41)-(2.43) should also

account for the neutral fluid pressure gradient. Hence the modified momentum equations

in the drift-diffusion form are

Electrons: neve = neE- V (Den)- Vpfl ,d (2.71)

Ions: nz, = n,,E V(Dn,)- VpflUd (2.72)

where Pfld is the pressure effect of the neutral fluid. Since we are using the two-fluid

model for the plasma-fluid mixture and with the plasma sound speed being much greater

than the neutral flow sound speed, the effect of the neutral fluid transport on the species

dynamics will be much smaller than the electric drift for the charged species which

makes the two pressure gradients in (2.71)-(2.72) negligible. However, the neutral

reactant species such as the metastables will still convect on the fluid dynamic scales and

evolve at a different rate. Therefore the modeling strategy can use a coupled solution

method where the above time-split integration of the plasma model is embedded into a

fractional step Navier-Stokes solver such as PISO [110]. In such an arrangement, the

corrector step just corrects for the fluid dynamics part of the plasma species equations.

This is an approach which is used in the modeling of reacting flows and is

computationally demanding. An approximate approach would be to neglect the effect of

the neutral fluid flow effects in the plasma discharge equations. This would entail a

discharge solution independent of the fluid dynamics which can be solved just once and

then be incorporated into the Navier-Stokes equations as shown in Chapter 4. Since the

most time consuming part of the solution algorithm is the modeling of the plasma






80


dynamics, the latter approach makes for a computationally efficient solution procedure.

In chapter 4 a framework for representing the plasma effects as a time-averaged body

force in the fluid flow is presented. An empirical expression for the body force based on

assumed discharge characteristics is used to model the plasma-fluid effects.














CHAPTER 4
MODELING FLUID FLOW WITH PLASMA EFFECTS

In this chapter a framework is presented to model the effects of the discharge on

the fluid dynamics. The model for the morphology and time evolution of the glow

discharge was presented in the previous chapter. In order to couple the plasma dynamics

to the neutral flow an appropriate strategy needs to be used. This usually depends on the

problem at hand and the operating conditions. For example, the flow physics needs to be

fully coupled and integrated with the discharge plasma simulation if the time scales of the

two processes are comparable. Specifically, the largest of the plasma timescales in Table

10 which is typically the inverse of the operating frequency, need to be comparable to the

mean flow time scale for the above to be true. However, for the low Reynolds number in

the incompressible limit we are interested in, the mean flow responds much slower than

the discharge evolution. As an example, the typical fluid flow time-scale for air over a 20

cm chord airfoil at Re=100 is of the order of a few seconds while the discharge operates

at a few KHz frequency. Hence, a decoupled treatment for the glow discharge with the

fluid flow solver would be a reasonable approximation here. In this chapter, a brief

description of the body force treatment of the glow discharge plasma will be presented.

This is followed by a detailed presentation of the phenomenological body force model

and the resulting glow discharge-induced flow effects. Here, an empirical formulation of

the body force based on an assumed plasma structure was used to study the resulting flow

patterns and the aerodynamics which revealed qualitatively similar characteristics to the

experiments. The next section presents the results obtained from the fully self-consistent









plasma dynamics calculations performed on the asymmetric electrode actuator operating

on Helium gas using a plasma-fluid model as discussed in the previous chapter. The

observed results provide some interesting insight into the mechanism of plasma

operation, structure and the resulting force generation. Also, the effect of the various

input voltage waveforms, operating parameters (the applied voltage and frequency),

electrode spacing and arrangements are studied.

4.1 Body Force Formulation

It is worth recalling that the momentum imparted to the fluid by the discharge is

through Lorentzian collisions of the heavier ion species with the neutral particles. We

will now revisit the right hand side of the ion momentum equation (2.21) which

represents the collision term. Although there will be different kind of collisions, the

primary collision loss will be that between the heavier ions and the neutral particles. If we

assume that the drift-diffusion approximation (2.25) for the momentum equation is valid,

then the net loss due to collision becomes

nqE- V (nkT) = vnmm, = FNozz(x,t) (3.1)

Here, the first term is the contribution from the electric drift and the second is the

gradient of the thermodynamic pressure which reduces to more of a diffusion flux if the

thermal effects are minimal. If we neglect the diffusion flux, then we get the loss due to

collision balanced by the Lorentz force acting on the charged ions. The electrons being

lighter are not significant players in terms of momentum transfer. Consequently the

electrons reach mean kinetic energy levels higher than those of the heavy components.

This loss due to collision F,,oo can be viewed as the instantaneous local body force acting

on the neutral fluid particles and can be modeled as a source term in the Navier-Stokes









equation. Once we have the instantaneous distribution of the ion number density ni and

the electric field E, the instantaneous body force Fcoll can be calculated as in (3.1). This

force acts on a much different time-scale than the characteristic time of the flow and it is

important that it be appropriately averaged over time. If we represent the integration

timescale by c and the time-averaged collision force by Ftave, we have


F (x, r) = F(x,t)dt
0


(3.2)


Dielectric barrier discharge modeling gives ni (
ion number density) and ne (electron number
density)


Calculate instantaneous force
F=E(x,y,t)e(ni-ne)


Average F over a flow time scale to get Ftave,x

I J
Substitute body force i Navier-Stokes, eg: x-momentum
eqn.
8pu pu2+p- Txx (pu v- y
+ + =F
8t 8x y tave ,


Solve flow field

Figure 12 Schematic illustration of the plasma-fluid flow modeling framework.

The choice of this time-scale r is determined by the smallest of the flow scales

which need to be resolved so that the corresponding effect of the force in that time scale

can be realized. It is worth noting that if c is much larger than the time period of the

operating voltage cycle, then the averaged force is independent of time. In other words,









the fluid will see the discharge as a source of constant average body force in time.

However, if the flow reacts faster than the time period of the driving voltage cycle (which

happens to be the largest plasma time scale) then the unsteady forcing effects of the

plasma body force will be felt by the flow and hence plasma transients need to be

modeled. This is important while modeling subsonic compressible flows or turbulent

flows with low plasma frequency. A schematic illustration of the modeling approach is

presented in Figure 12.

4.2 Navier-Stokes Equations with Body Force

The numerical model consists of the continuity equation and the two dimensional

momentum equations for a steady incompressible viscous flow. The body force terms,

which are added to the momentum equations, carry the effect of the plasma discharge on

the fluid flow. The fluid is assumed to be incompressible in view of the plasma induced

jet being a low Reynolds number and essentially isothermal phenomenon. In the

following, we offer a complete set of governing equations, in the two-dimensional form.

0A 0B C -
-+ = D (3.3)
8t 8x ay








B = pu2 + p-x (3.5)
puv r


pv
C = puv- r (3.6)
pV' + P ryy









0
D= F, .. (3.7)
Ltave,y

where cxy is the shear stress and Ftave,x ,Ftave,y are the body force components in a 2-D

case, as given in (3.2).

4.3 Modeling a Plasma Actuator Using a Phenomenological Force Model

Here we try to model the plasma actuator used by Roth et al.[4] and Corke et al.[25] with

dimensions as shown in Figure 13 where Lu, L1 and d are the length of the upper, lower

electrodes and gap distances respectively. The plasma generated in this instance is

weakly ionized and non-thermal (- room temperature) where the electrons alone have a

large mean kinetic energy. The ions are more or less in thermal equilibrium with the

neutral species because of the extensive collisions with neutrals in a weakly ionized state.

A phenomenological plasma model without a detailed simulation is used to depict the

glow discharge. The motive of this study was to develop a framework which would take

into account detailed balances between momentum transfer mechanisms such as

convection, hydrostatic pressure, viscous stresses and the body force arising from the

Lorentzian collisions. Now we present a brief phenomenological discussion of the plasma

generation in an asymmetric electrode configuration based on our understanding of the

plasma structure in a high frequency parallel plate (Massines et al.[6]) arrangement.

4.3.1 Plasma Generation in an Asymmetric Electrode Configuration

At these high frequencies the charged species (mainly electrons) attain high

velocities, leaving a very short time for the species to interact and recombine. The ions,

being less mobile are believed to be trapped between the electrodes without actually

reaching the electrode surfaces. This ion-trapping mechanism is one of the keys for