UFDC Home  myUFDC Home  Help 



Full Text  
COMPUTATIONAL MODELING OF GLOW DISCHARGEINDUCED 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 PlasmaInduced 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 PrinciplesBased 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 1D 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 NavierStokes 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 SelfConsistent FluidModelBased 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 plasmainduced 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 1D 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 1D Numerical modeling results showing different regions of the discharge gap at the instant when the current is m maximum ........................................... ............... 18 4 A 3D 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 NACA6618 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 plasmafluid 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 AB 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 2D 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 3Different configurations to study the impact of electrode overlap..........................126 39 Time evolution of species number density over one timecycle (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 powerlaw 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 Timeaveraged 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 DISCHARGEINDUCED 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 nearwall 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, nonequilibrium plasma discharges in conjunction with low Mach number fluid dynamics at atmospheric pressure. The plasma and fluid species are treated as a twofluid 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 NavierStokes 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 selfconsistent manner using a first principlebased 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 finitevolume operatorsplitting algorithm capable of conserving space charge is developed; the approach can handle timestep sizes in the range of the slowest species convection timescale. The NavierStokes equations representing the fluid dynamics are solved using a wellestablished pressurebased algorithm. A onedimensional twospecies 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 powerlaw 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 plasmafluid 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 thermofluid effects of the discharge plasma actuator (see Figure 1) which generates efficient surface plasma using a dielectric barrier arrangement capable of interacting with the nearwall 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 dischargeinduced 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 fluidlike 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 midnineteenth 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 electrohydrodynamic (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 buildup at the dielectric surface in a short time after breakdown, the electric field in the discharge region is reduced, thus making it selflimiting. 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 [47]. 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 reemerging 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 nonmoving 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 principlebased highfidelity as well as phenomenological approaches) to study the momentum generation and the thermofluid 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 nearwall 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 plasmafluid problem is inherently nonlinear 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 convectiondiffusionreaction 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 astrophysics; interstellar gas masses etc and is studied under magnetoaerodynamics. 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 nonequilibrium 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 nonequilibrium 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 nonequilibrium 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 nonequilibrium 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 nonequilibrium 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 nonequilibrium effects are minimal and, moreover, the plasma dynamics is mainly shaped by the electromagnetic effects and to a lesser extent by the background lowspeed 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 nonequilibrium, low pressure, high frequency (RF) discharge conditions) to relatively simpler plasmafluid 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 nonlinearities between the plasmafluid 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, 2D or higher dimensional discharge modeling studies are less common. More recently Roy et al.[1418] have developed a finiteelement method based 2D modeling of such plasmabased actuators by solving for the discharge species and neutral transport. Although the model used only three species (ions, electrons and neutrals) in a 2D 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 plasmabased separation control in a NACA 0015 wing section. This model was an attractive option because of the difficulty in achieving efficient multidimensional and selfconsistent plasma dynamics simulations coupled with largely unsteady fluid dynamics. This was followed up by a first principlesbased selfconsistent modeling of the plasmadynamics using a hydrodynamic model in helium gas. This multifluid formulation to model the radiofrequency 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 selfconsistent way of modeling the plasmafluid 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 fullycoupled multidimensional plasmafluid 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, multidimensional physics. Both sequential and fully implicit approaches [15] have been used for discharge modeling. In the present multiple scale problem, we have employed an operatorsplit 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 multiblock finitevolume algorithm capable of handling 3D 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 2D. 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 largescale 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 timesplit 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 NavierStokes 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 plasmafluid equations. Having developed the modeling framework, a onedimensional 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 1D plasma model results, is used here to investigate qualitatively the induced flow structure. As a sequel to the above, the fully selfconsistent 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 nonexistent 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,2127]. However, numerical studies to complement the above experimental efforts are not too common until recently [2829,1418]. 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 thermofluid 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 dischargeinduced 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 110 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 nonequilibrium and cold with timedependent 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 nonequilibrium, 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 PlasmaInduced 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 [3233]. 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 onedimensional 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 1D 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.[2122] 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 1D 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 negativegoing and highly irregular when positivegoing. 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 sawtooth waveforms having vastly different positivegoing and negativegoing duty cycles were used to model the discharge. The resulting outcome clearly established the preferential phase for the actuator operation during a whole timecycle. 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 (510 KV, 120KHz [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 110 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 1D 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.[2122], 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 followup 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 ' coherentlydriven 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 pumpingtype mechanism resulting in a walljet. However, hotwire temperature and velocity measurements by Jukes et al.[3738] 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 glowdischarge 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 densityvariation 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 recirculating 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 3D 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,3940] 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 plasmabased 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(104105). Goskel et al.[42] did an experimental investigation of steady and unsteady plasma wall jets for separation control of Microair 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 [2123] 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 reattachment 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 Re50000 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 hotwire 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 NACA6618 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 flowcontrol 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 timescale 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.[4547]. 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 wallbounded 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 plasmabased 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 highfidelity first principlesbased 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 dischargeinduced fluid dynamics while, in section 2.2.2 we will present some of the more highfidelity 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 plasmainduced fluid flows. This model involved the development of an empirical approach to calculate the timeaveraged 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 analyticalempirical model (see Chapter 4) based on a linear force distribution in the domain was later adopted by Gaitonde et al. [14] for modeling plasmabased separation control in a NACA 0015 wing section. This linearized body force model was an attractive option because of the difficulty in achieving efficient multidimensional and selfconsistent 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 (sourcesink) 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 lumpedelement 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 plasmafluid interaction. 2.2.2 First PrinciplesBased 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 dischargecoupled 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 twodimensional 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 109s [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 multispecies 1D 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]. Twodimensional 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 twodimensional 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 twodimensional 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 nonuniform 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 preionization through penning ionization mechanisms. This preionization 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 Multispecies 1D 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 Multispecies 2D model Modeling ofDBD's using a . Dongsoo Lee at al, (2004) [34] eling of D uing 1D multispecies fluid model heliumoxygen mixture DBD modeling using HeAr Massines et al, (1998) [6] D modeling using HeAr 1D twospecies fluid model mixture DBD modeling using HeAr Massines et al, (2000) [58] m e u 1D multispecies model mixture DBD modeling of pure . Tochikubo et al, (1999) [63] DBD modeling of pure 1D multispecies model Helium .S \Modeling of Helium Mangolini et al, (2004) [62] Modeling of Helium 1D multiSpecies model oxygen mixture Shin et al, (2003) [61] Modeling of HeAr mixture 0D, 2species model Modeling of nitrogen based . Xu et al, (1998) [59] i i e 2D multi species model mixture Eliasson et al, (1991) [54] Modeling in Xenon 2D multispecies 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 onedimensional either with a fluid model or a particle model of the plasma. Since then, two and threedimensional studies have been successfully presented. A review of onedimensional studies is presented by Govindan et al.[64]. The dominant type of plasma considered in these studies are low pressure, lowtemperature and weakly ionized. High pressure glow discharges have also been modeled lately by Yuan et al.[65]. The fluid model is based on the selfconsistent 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 2D modeling study is shown in Figure 6. The wavelike solution and sharp gradients are characteristic of streamer propagation and also low frequency discharges. The system of equations is also highly nonlinear 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] 2D fluid model 2D fluid and particle model Nitschke et al, (1994) [68] comparison comparison 2D 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] 2D fluid model for a PDP 2D fluid model for Gaseous electronic conference (GEC) reactor comparison with experiments 1D fluid model for PDP's 1D 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.[1516] proposed a selfconsistent twodimensional DBD fluid model for helium gas with application to separation control using finite element techniques. This multifluid formulation to model the radiofrequency 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 selfconsistent way of modeling the plasma wall jet interaction. This study employed a finiteelement 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 timecycle 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 abovementioned modeling studies assume negligible electrode thickness. In the event of a finiteelectrode thickness, the treatment of the dielectric electrode edge can impact the nearwall 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 twofluid 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 dischargedinduced 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 walljet 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 (05mm) Figure 6 Streamer front propagation at different instants in a DC voltage (Dhali et al.[69]) Table 6 Summary of plasmainduced fluid flow modeling studies Author and year Roy et al.(2005) [1516] Singh et al, (2005) [17] Kumar et al, (2005) [18] Singh et al, (2006) [75] Visbal et al, (2006) [76] Work description Selfconsistent 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 3species model ions, electrons and neutrals for Helium discharge. Same as above. Same model as above. 8species 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 multiscale 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 onedimensional 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, multiscale and highly nonlinear plasma equations coupled with the electrostatics. The fundamental numerical model describing the high and lowfrequency 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, 2D). 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) MonteCarlo Particleincell 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 1D2V 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 selfconsistent 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 MonteCarlo (MC) simulation in rf Kushner [78], 1983 discharges. Particleincell (PIC) simulation of parallel Boswell et al.[79], 1988 plate RF discharge Surendra et al.[80], 1991 PICMC 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 timeaveraged 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 nonlinearity. 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 particleincell (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 selfconsistent simulation of parallel plate discharges. Birdsall [67] and Surendra et al.[80] study the PICMCC (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 selfconsistently 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 electronneutral, ionneutral 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 quasisteady 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 MonteCarlo algorithm is employed to model the chargeneutral 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 particleincell and MonteCarlo technique (PICMC) 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. MonteCarlo 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 timeaveraged 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 lowpressure (typically < 100mTorr [68]) or highly nonequilibrium 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 'Bulkbeam' 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 'Bulkbeam' model Fluid model with MonteCarlo 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 lowpressure 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 nonequilibrium 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 selfconsistent simulation of the entire discharge while addressing the nonequilibrium 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 recalculated by a Monte Carlo simulation after a few hundred timesteps 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 lowpressure discharges. Nitschke et al.[68] compare the particleincell 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 ny, 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% , 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 noncharged 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 onedimensional (see Table 4), followed by twodimensional 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 1D 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 PICMC, 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 plasmafluidMaxwell 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 uD ,E = v, (2.29) This form of the momentum equation is called as the driftdiffusion 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 nonlocal 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 coworkers [8889]), the fluid energy equation is dropped and the kinetic energy equation is solved using a 2D 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 electronneutral 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 nonlinear 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 timeintensive computations. Sharply varying production terms such as ionization processes result in stiff source terms which also limit the time step. For these reasons, 3D and to a lesser extent 2D 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 multitimescale 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 nonlinear 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 driftdiffusion 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 timestep 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 semiimplicit 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 1D 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 1D applications and occasionally 2D 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 driftdiffusion form and the full transport equation. Both these forms are used in the modeling studies, but the simpler driftdiffusion 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 1D 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) le4 Ion drift rdron (2.48) 2.5e7 Electron drift drelectron (2.49) 2.5e9 Dielectric relaxation Tel (2.55) 2e9 Ionization timescale rc (2.50) 1.4e10 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 2D 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 lowpressure, high density rf discharges in Helium 1D parallel plate discharge in helium at atmospheric pressure 2D 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 multiscale 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 timeaccurate and in some nonlinear situations are even computationally burdensome. This is necessitated by the need to find the root to a highly nonlinear equation which might entail very small timesteps 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 NewtonKrylov methods, which have been reviewed in [92]. Both sequential and fully implicit approaches have been used for discharge modeling. Roy et al.[1516] 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 RungeKutta method followed by an implicit treatment of the fast electron equations. The implicit Euler (first order) and implicit RungeKutta (second order) methods were used with NewtonRaphson iterations to overcome the non linearity. Further, the resulting Jacobian is simplified by neglecting the weak off diagonal contributions reducing to a block tridiagonal 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 semiimplicit 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 nonstiff 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 nonstiff process. These are essentially based on the backward differentiation formulae of Gear. This integration timescale in most problems is the convection timescale of species which is considerably larger. Such a procedure makes for a sequential approach and thereby enables handling smaller time scales by allowing integration subcycles 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 semiimplicit 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 nonlinearity. However, with the use of efficient solution techniques such as implicit multigrid methods and efficient solvers like LUSGS 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 multistage Rungkutta schemes, RungeKuttaChebychev (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 timescale 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 timestep 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 timesplit 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 predictorcorrector approach wherein the evolution of the various species is appropriately coupled to the electric field. In the present multiple study, we have employed an operatorsplit 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 timescale than dictated by the individual processes. Such dynamic equilibrium usually cannot be predetermined in general for a nonlinear advectiondiffusionreaction 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 timesplit algorithms for processes operating in a range of timescales, the choice of timestep size is typically determined by the smallest time scale, but need not necessarily be chosen as such. To speedup the solution procedure, an intermediate time scale is chosen to advance the overall system in time, while the faster processes are advanced by subcycling within the timestep. In the present study, the timestep dictated by the slower ion species convection is targeted to march the full discharge system while subcycling is used for the faster processes. Also, a predictorcorrector approach is employed to ensure sufficient coupling between the electric field and the species densities. A strong coupling is essential for achieving stable timeaccurate simulations while using a sufficiently large global timestep. This method is integrated with a multiblock finitevolume algorithm capable of handling 3D curvilineargrids [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 driftdiffusion 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 higherorder (4thorder) 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 secondorder upwind for the convection term and secondorder 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) Firstorder splitting (ii) Strang splitting (iii) Source splitting. (i) Firstorder splitting: The firstorder 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 secondorder 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 firstorder and Strang splitting procedures the initial guess for the reaction part is not directly from the previous timestep, but after a half or full timestep 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 nonlinearities. Even though the Strang splitting is formally secondorder, 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 firstorder 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 5thorder BDF for time integration. Newton iteration for nonlinearity. A direct method with a banded treatment of the Jacobian. Normal mode with subcycling within the timestep. Relative and absolute tolerances of le12 and le14 respectively. The above strict tolerances were chosen so that the ODE integration is almost exact. After the predictorstep, we have obtained the guessed values for the species number densities based on the lagged value for the electric field. The strong superlinear 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 semiimplicit version of the Poisson equation to overcome the space charge stability constraint leading to the dielectric relaxation timestep 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 nonimplicit 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.E1 (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 driftdiffusion 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 nAt 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 driftdiffusion 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 nonlinearity will be difficult to overcome and will need a NewtonRaphson 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 driftdiffusion 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 predictorcorrector 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 nonlinearity. The predictorcorrector formulation with the semiimplicit solution of the Poisson equation ensures space charge stability overcoming the dielectric relaxation time restriction due to the nonlinearity. 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 timestep 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 timestep 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 timestep for the faster electron species and the corresponding number of substeps 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 substeps M, employed. Knio et al.[108] shows that the critical CFL number due to electron convection substepping decreases monotonically as the number of substeps increases in twodimensional modeling studies of reacting flows. Also, onedimensional studies indicate a limiting value of the critical CFL being achieved as the number of substeps is increased. In spite of the stability criterion becoming stringent the overall computational savings is generally substantial. 3.2.3 1D Parallel Plate Discharge Modeling in Helium We study the following test problem to validate our computational capability. A 1D parallel plate (5mm gap between electrodes) discharge at atmospheric pressure using just twospecies 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 driftdiffusion 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 1D 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 timestep of 109 seconds is used for the timeintegration and the domain is uniformly discretized into 201points. 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 quasisteady state has been attained during the discharge is selfsustained 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 Icr d Figure 8 Parallel plate arrangement for modeling DBD discharge 150 present201 100 S Massines et al .1 ...... present401 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 SionMassines electron Massines B 4 3 S2 Axial direction (00.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 NavierStokes 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 driftdiffusion 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 twofluid model for the plasmafluid 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 timesplit integration of the plasma model is embedded into a fractional step NavierStokes 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 NavierStokes 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 timeaveraged 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 plasmafluid 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 timescale 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 dischargeinduced 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 selfconsistent plasma dynamics calculations performed on the asymmetric electrode actuator operating on Helium gas using a plasmafluid 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 driftdiffusion 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 NavierStokes 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 timescale 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 timeaveraged 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(nine) Average F over a flow time scale to get Ftave,x I J Substitute body force i NavierStokes, eg: xmomentum eqn. 8pu pu2+p Txx (pu v y + + =F 8t 8x y tave , Solve flow field Figure 12 Schematic illustration of the plasmafluid flow modeling framework. The choice of this timescale 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 NavierStokes 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 twodimensional form. 0A 0B C  + = D (3.3) 8t 8x ay B = pu2 + px (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 2D 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 nonthermal ( 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 iontrapping mechanism is one of the keys for 