NUMERICAL MODELING OF PLASMA MULTI-ACTUATOR SYSTEM

The spatially periodic system of dielectric barrier discharge actuators intended for laminar flow control on a swept wing is simulated numerically. Novel scheme of the actuator system ensuring the discharge formation only near one edge of every exposed electrode is proposed. The mathematical model of dielectric barrier discharge in air is formulated in the drift-diffusion approximation without convective transfer with taking into consideration the following volumetric reactions: the ionization of nitrogen and oxygen by electron impact, the electron attachment to oxygen, the electron detachment from negative ions of oxygen, the ionion recombination, the electron-ion recombination. Two types of boundary conditions on dielectric surface are considered for comparison: the model of instantaneous recombination and the model of finite rates of recombination and electron desorption. Numerical simulation of both conventional and novel schemes of plasma actuators is carried out for one set of problem parameters. Enhanced energy efficiency of the proposed scheme is demonstrated.


INTRODUCTION
The concept of laminar flow control (LFC) method on a swept wing based on an attenuation of the cross-flow-type instability due to electrogasdynamic (EGD) force impact on three-dimensional boundary layer in the vicinity of a wing leading edge was proposed at TsAGI some years ago [1].This method is illustrated on Fig. 1.Practical realization of the proposed method seems to become realistic owing to the usage of plasma actuators operating on the base of dielectric barrier discharge (DBD) [2].Volumetric force F || generated in DBD and directed along a wing leading edge in xdirection will induce gas velocity in this direction thereby attenuate cross-flow velocity component in boundary layer.In turn, an attenuation of the cross-flow velocity results in a decrease of increments of spatial growth of steady-state modes of the cross-flow-type instability in transonic three-dimensional boundary layer [3].
Moreover one can expect that a "filling" of the velocity profile in viscous flow along a critical line of a swept wing due to flow acceleration in this direction will result in some rise of the stability of this flow and, hence, an increase in the critical Poll Reynolds number.
Practical realization of this LFC method demands development and optimization of multi-actuator system spatially periodic along a wing leading edge and ensuring EGD impact on large part of a wing surface.The optimal system must create a volumetric force F || large enough for prevention the laminar-turbulent transition caused by the cross-flow-type instability or the leading edge contamination at moderate electric power consumption.Preliminary theoretical estimations demonstrate possible damaging influence of volumetric heat release in DBD on boundary layer stability when the coefficient of energy efficiency of DBD-actuators is too small [4].In turn the coefficient of energy efficiency diminishes with a decrease of gas pressure because this coefficient is inversely proportional to the average drift velocity of charge carriers [5,6].Besides a miniaturization of DBD-actuators is necessary for practical realization of the mentioned LFC method because the boundary layer thickness in the vicinity of a swept wing leading edge of civil transonic airplane during cruise flight is very small (about 0.5 mm).
Experimental study of LFC on a swept wing in wind tunnels is very difficult through a necessity of simulation cruise flight conditions, namely, transonic Mach number, high Reynolds number, low turbulence level, low static pressure.In addition DBD-actuator system consists of a lot of geometrical and physical parameters influencing on its performance characteristics, first of all, volumetric force and heat release.Therefore parametric experimental study of plasma actuator system is very labor-consuming.All these reasons make actual preliminary numerical simulation of mentioned systems.
Last years several tens works devoted to both theoretical and experimental study of DBD-actuators are published annually [7,8].Numerous numerical modeling, as a rule, deal with a single actuator [9][10][11].The most detailed experimental investigations are executed with single actuator too [12,13].Experimental optimization of multi-actuator systems is more complex as compared to single actuator because of a presence of additional parameters, for example, the phase shift of the applied voltage on adjacent actuators, the distance between them, the widths of the exposed and insulated electrodes [14,15].Moreover the optimal system has to ensure a discharge formation only near one edge of every exposed electrode and exclude it at opposite edge.Main characteristics of various variants of the considered system of plasma actuators may be evaluated and compared on the base of numerical simulation with the use of proper mathematical models and computer codes.

STATEMENT OF BOUNDARY VALUE PROBLEM
The scheme of the proposed spatially periodic system of DBD-actuators is shown on Fig. 2. The electroconductive wing structure A and the main buried electrodes B are under constant zero electric potential, and alternating electric potential is applied to the exposed electrodes C and the additional buried electrodes D. Electrodes are separated by dielectric layers 1 and 2. If the alternating potential is of a large enough amplitude, DBD turns on in the vicinity of the right edges of the every exposed electrode.The oscillating volumetric electrostatic force arises in these regions due to charge separation in the discharge gap.The time-averaged horizontal component F || of this force is directed from left to right and accelerates boundary layer flow in this direction.The additional electrodes D permit to attenuate drastically the electric field strength near the left edges of the exposed electrodes thereby to avoid the discharge ignition in these regions.The usage of the electrodes D permits to reduce significantly the width of the exposed electrodes, the distance between them, and the thickness of the internal dielectric layer 1. Note, that the width of the exposed electrode is taken large enough or one its edge is covered with insulating film to avoid the formation of discharge at that edge in usual schemes of DBD-actuator with one buried electrode B [15].Numerical modeling of one spatial period of the described system is fulfilled in the calculation domain shown on Fig. 3 and consisting of subregions 1-4.Here C r is the right half of some exposed electrode and C l is the left half of the next exposed electrode.The thickness of the insulated electrodes B and D is taken to be zero.
Here ε 0 is the permittivity of vacuum, e is the elementary charge, are the volumetric concentrations of positive ions, negative ions, and electrons, respectively.Relative dielectric constant of air is taken to be equal 1.The external subregion 4 is introduced to satisfy the condition of zero electric potential at infinity.The following transformation of ycoordinate is used here to move the external boundary y=y ∞ at infinity: y / =y e (y ∞ -y e )/(y ∞ -y).
The periodicity conditions are stated on the left and the right boundaries of the calculation domain.The conditions of continuity of the electric potential and the induction are stated on the boundaries between the dielectric layers 1 and 2 and between the subregions 3 and 4. The potential continuity condition and the relation between the discontinuity of the electric induction and the superficial charge are stated on the upper surface of the external dielectric layer 2. Taking into account zero potential of the electrodes A and B and specified potential of the electrodes C and D the boundary conditions for the electric potential take the following form: Here Ω 12 denotes the boundary between the dielectric layers 1 and 2 except the electrodes B and D, Ω 23 denotes the open surface of the dielectric layer 2, Ω C and Ω D denote the surfaces of the exposed electrodes C r , C l and the insulated electrode D, respectively, ε 1 and ε 2 are the dielectric constants of the layers 1 and 2, σ is the superficial charge on the surface Ω 23 , V 0 and f are the amplitude and the frequency of the applied harmonic voltage.
The dynamics of charged particles is governed by three continuity equations for positive ions, negative ions, and electrons without taking into consideration a convective transfer [ Here Γ i (i=p, n, e) are the vectors of the drift-diffusion fluxes of charged particles,  and D are the coefficients of mobility and diffusion, N (cm -3 ) is the volumetric concentration of air molecules, k i is the ionization rate coefficient taking into account an origin of positive ions of both nitrogen and oxygen, k iN and k iO are the individual coefficients of nitrogen and oxygen ionization, k dr is the rate coefficient of the electron-positive ion recombination, k r is the rate coefficient of the ion-ion recombination, k at is the rate coefficient of attachment of electrons to oxygen molecules, k dt is the coefficient of detachment of electrons from negative ions, γ=10 16 E/N (V•cm 2 ) is the function of the reduced field, E (V/cm) is the absolute value of the electric field strength, T e and T p are the temperatures of electrons and ions.The reactions rate coefficients k iN , k iO , k dr , k r , k at , k dt are measured in cm 3 /s.The coefficient of the electron mobility μ e and the electron temperature T e depending on the reduced field function γ are calculated by the method described in [16].The mobility coefficients of ions are defined as μ p,n =C p,n (N 0 /N)(300/T p,n ) 0.3 (cm 2 V -1 s -1 ), C p =2.1, C n =3.2, N 0 =2.69•10 19 cm -3 .A distinction of ions T p,n and neutrals T temperatures is neglected.The diffusion coefficients D i (i=p, n, e) are calculated on the base of the Einstein relation using the mobility coefficients and characteristic temperature for electrons [16].The numbers 0.8 and 0.2 correspond to the content of the nitrogen and oxygen molecules in the total concentration of neutrals N.
The equation system (3) is solved in the subregion 3 only.The initial conditions for these equations are taken as follows: The periodicity conditions are stated on vertical boundaries of the subregion 3 The approximate conditions used on upper boundary imply negligible diffusive fluxes The boundary conditions for normal fluxes of ions and electrons Γ wi (i=p,n,e) on both dielectric and metallic surfaces are deduced in [16] with the use of the relaxation approximation of Boltzmann equation for the velocity distribution function of charged particles and the approximated Maxwell approach instead of more strict consideration of the Knudsen layer.The physical model of interaction between charged particles and dielectric surface in nitrogen [17] was completed for air and used for a statement the equations governing the accumulation of ions and electrons on dielectric surface.The following physical processes have been taken into consideration: full adsorption of all incident charged particles due to polarization forces, neutralization of positive ions with the origin of the superficial positive charge and emission of secondary electrons, desorption of electrons and negative ions with a finite rate, electron-ion and ion-ion recombination with a finite rate, detachment of electrons from negative ions.Boundary conditions on dielectric surface depend on the sign of the normal component of the electric field strength, but this influence is weak.Therefore the fluxes of charged particles and their superficial densities are governed at by the following equations: Here ν i (i=n, e) is the frequency of desorption of negative ions or electrons, ζ d is the coefficient of secondary ion-electron emission on dielectric surface, k dt is the same coefficient of detachment rate as in (3), c i (i=p,n,e) is the thermal velocity of corresponding charged particles, α re and α rn are the constants of the recombination rate of electrons and negative ions, respectively, with positive superficial charge.It is supposed that these constants are proportional to the mobility of proper negative particles.That is the relation α rn =α re μ n /μ e is used.The superficial charge in condition ( 2) is determined as The more simple boundary conditions on dielectric surface are stated in the framework of the model of instantaneous recombination [9].In this case the expressions for Γ wi and the equation of superficial charge variation take the form .

The boundary conditions for fluxes Γ wi (i=p,n,e) on open surfaces Ω C
/ of the external electrodes C r and C l remarkably depend on the electric field strength on these surfaces E w because the field reaches here its extreme values.Positive E w decelerates incident positive ions and secondary electrons, negative E w decelerates incident negative ions and electrons.These conditions take the following form [16]: Here ζ m is the coefficient of the secondary ion-electron emission on the electrode surface, λ i , c i , m i (i=p, n, e) are the mean free path, the thermal velocity, and the mass of corresponding particles, k B is the Boltzmann constant.The functions A and B reflect the decelerating effect of the electric field, A i and B i (i=p, n, e) relate to incident particles, A eζ relates to secondary electrons.It is supposed that the secondary electrons have the mean energy of 0.5 eV.The mean free paths of ions and electrons presenting in the argument of A and B are estimated according to formulas It must be emphasized, that the formulated mathematical model is not capable to describe streamer phenomena always taking place in DBD because these processes are strongly three-dimensional.However some justification for this model using consists in assumption that streamers influence significantly on total electric current and Joule dissipation but weakly on volumetric force because ionized gas is quasi-neutral in streamer channels.The stage of DBD saturation is explained by this feature of streamers [15].Therefore it seems evident that calculated Joule dissipation will be underestimated and, hence, the energy efficiency coefficient will be overestimated.

RESULTS OF CALCULATIONS
The formulated boundary value problem was solved by the finite element method with the use of the first order Lagrangians and strongly nonuniform triangular mesh.The total number of finite elements for considered geometry was about 2•10 5 .Results of numerical simulation of both conventional (without additional buried electrode) and novel systems of DBD-actuators carried out for one set of the problem parameters are presented and discussed below.
The geometric parameters indicated on Fig. 3 are the following: x 1 =0.1, x 2 =0.7, x 3 =0.8,x 4 =0.9, x e =1, y 1 =-0.3, y 2 =-0.1, y e =1, y ∞ =1.3 cm.The height of the exposed electrodes equals 55 μm and they are deepened in the dielectric layer 2 up to 5 μm.Both upper and lower verges of these electrodes are rounded with the radius of curvature 5 μm to eliminate electrostatic singularity arising on sharp corners.
The following physical parameters were used in calculations: the pressure and the temperature of air p=760 tor and T=300 K, the amplitude and the frequency of the applied voltage V 0 =7 kV and f=5 kHz, the dielectric constants of the insulation layers ε 1 =2.1 (Teflon) and ε 3 =3.5 (Kapton), the coefficients of the secondary ion-electron emission on the dielectric and the electrode surfaces ζ d =ζ m =0.05, the frequencies of electron and negative ion desorption from dielectric 5   n e   kHz, the electron-ion recombination coefficient on dielectric surface 7 10   re  cm 2 /s [17].The discharge behavior strongly depends on superficial charge distribution on the dielectric surface.This distribution becomes quasi-periodic in time only after several cycles of the applied voltage [12,16].Calculation of every cycle of DBD is very time-consuming.Therefore conclusions on main characteristics of DBD-actuators, namely, total volumetric force, Joule dissipation, and energy efficiency coefficient are made on the base of results of second cycle calculations.
Dependences of the conduction current are shown on Fig. 4 for conventional and novel schemes of plasma actuator system.Dashed curves here and then represent the electric potential of the exposed electrodes.The conduction current is calculated as integral over the open surfaces of the exposed electrodes C r and C l on fluxes of ions and electrons with taking into account the sign of charge carriers.The sign of the horizontal volumetric force in the considered geometry during both strokes of DBD is determined by positive ions.It seems to be in contradiction with known experimental data concerning critical role of negative ions in DBD in air [13,18].One possible explanation of this contradiction may be related with a scaling factor.Sizes of single plasma actuator studied in mentioned experiments were essentially greater than in the present modeling.Other reason may be connected with the overestimated value of the electron detachment rate coefficient dt k [19] used in the equations (3).According to the present calculations, the detachment of electrons and corresponding vanishing of negative ions play a dominant role in intervals between micro discharges during the forward stroke.In turn, additional preliminary calculations executed show that F || becomes positive during the forward stroke of discharge in novel system, in contrast to Fig. 5, b.These problems related with both scaling factor and detachment rate coefficient demand additional study and will be investigated in near future.During positive-going electric potential (the backward stroke of DBD) the discharge at first arises near the right edge C r in both systems.But then with the potential increase the electric field strength near the left edge C l in conventional system becomes sufficient for discharge ignition near this edge.The horizontal volumetric force generated near the left edge C l during the backward stroke is negative.Therefore the total horizontal force is remarkably less in conventional system as compared to novel one, as it is seen on Fig. 5.
Time dependences of Joule dissipation defined by the second expression (10) are shown on Fig. 7.It seems evident that Joule dissipation in convectional system is greater during both strokes of discharge because of two active edges on every exposed electrode.Maximal values of Q•during the backward stroke are much less than in forward stroke and equal about 75 W/m in conventional system and 50 W/m in novel one.About 40% of total power is spent during the forward stroke in both systems.
One of the main parameter of DBD-actuator characterizing its working capacity is the coefficient of energy efficiency E=<F || >/<Q> [20], where average volumetric force and dissipation is calculated according to The results discussed above and presented in second and third rows of Table 1 were obtained with using more complex boundary conditions on dielectric surface (7).For comparison calculations were executed with simpler model of instantaneous recombination (8) and appropriate results are presented in last row of the Table .Note, that according to experiments with a single DBD-actuator the coefficient of energy efficiency varies in the range (1÷7)•10 -4 s/m depending on geometric and physical parameters of actuator [21].Hence the calculated value of E for novel system may be considered as high enough.

CONCLUSION
Presented results on numerical modeling of DBD multi-actuator system demonstrate that the discharge formation at upstream edges of the exposed electrodes can be excluded.It can be achieved by installation of additional buried electrodes beneath these edges with the same electric potential as the exposed electrodes.The numerical estimations show that fairly large coefficient of energy efficiency can be obtained in the considered scheme of the multiactuator system.

Figure 1 .
Figure 1.The sketch of LFC method on a swept wing.

Figure 2 .
Figure 2. The scheme of DBD multi-actuator system.

Figure 3 .
Figure 3. Calculation domain.The electric potential φ and the electric field strength

Figure 4 .Figure 5 .
Figure 4. Conduction current in conventional (a) and novel (b) systems of plasma actuators.During the forward stroke of DBD when the potential of the exposed electrodes falls the discharge consists of a series of shot pulses (micro discharges) of duration about ( 5 2  )•10 -8 s separated by intervals of duration ( 5 2  )•10 -6 s.The current peaks in micro discharges are of the order of 1 A/m.During the backward stroke of DBD when the potential of the exposed electrodes grows the current is smoothed.Maximal values of the current during this stage achieve about 0.02 A/m in conventional scheme and 0.015 A/m in novel one.The principal distinction between two schemes of plasma actuator becomes clear in time dependences of spatial integrated horizontal force.The total horizontal force and Joule dissipation generated in one actuator of the considered systems are calculated as integrals over the area S 3 of the subregion 3 according to    dxdy e Q dxdy n n n eE F S e n p S e n p x         3 3 ,

Figure 6 .
Horizontal volumetric force in conventional (a) and novel (b) systems during the first quarter of the second cycle.

Figure 6 ,Figure 7 .
Figure6, a demonstrates distinctly that 1-2 micro discharges initiated near the right edge C r of the left exposed electrode is followed by one micro discharge arising near the left edge C l of the right exposed electrode in conventional system.Namely positive values of F || indicate the discharge ignition near C l .Negative superficial charge arising on dielectric surface near the right edge C r due to electrons deposition during micro discharge near C r attenuates significantly the electric field strength in this region.As a result, after 1-2 micro discharges taking place near C r the field strength becomes greater near the left edge C l after appropriate decrease of applied electric potential.The additional buried electrode D screens the left edge C l and prevents the discharge ignition here, as it is obvious from Fig.6, b. 10]

Table 1 .
Less horizontal volumetric force and greater Joule dissipation in conventional system stipulate essentially less energy efficiency coefficient, as one can see in Table1.Main characteristics of plasma actuator systems.