A simplified formulation of wire-plate corona discharge in air: application to the ion wind simulation

1 A simplified formulation of wire-plate corona discharge in air: application to the ion wind simulation M R Bouazza1, K Yanallah1, F Pontiga2 and J H Chen3 1Laboratoire de Génie Electrique et des Plasmas, Université de Tiaret, Algeria 2Departamento de Física Aplicada II, Universidad de Sevilla, Spain 3Department of Mechanical Engineering, University of Wisconsin-Milwaukee, Milwaukee, WI 53211, USA Abstract. The spatial distribution of charged particles (electrons, negative ions and positive ions) and electric field have been evaluated using a semi-analytical approach of the positive and negative corona discharge for a wire-to-plate electrode system. Thus, approximate formulas useful for the characterization and control of corona discharge devices are provided, which helps to significantly reduce computational costs. Based on the obtained results, the electro-hydrodynamic (EHD) force generated by the corona discharge has been determined, and it has been used in the Navier-Stokes equations to compute the spatial distribution of the gas velocity. As a result, the influence of the corona plasma region in the flow pattern, particularly in the vicinity of the corona electrode, has been brought to light, which helps to understand the different flow velocities observed in positive and negative coronas. Moreover, the influence of voltage, wire radius, and inter-electrode separation on the electric wind velocity has been investigated. Keywords: Corona discharge; electric wind; semi-analytical approach; wire-to-plane electrode geometry; fluid model. 1. Introduction Corona discharge occurs when high voltage is applied between two electrodes, and at least one of them has a sharp curvature. A strong electric field is then created around the sharp electrode, causing a local breakdown and the flow of electric current between the two electrodes. When charged carriers (electrons and ions) traversing the inter-electrode space collide with neutral gas particles, they transmit a part of their momentum to them, which sets the gas in motion. The flow thus generated is generally referred to as electric wind, ion wind or EHD flow. This phenomenon has been known for a long time [1], but it still remains under study due to its complexity. The electric wind can be generated in corona discharges as well as in other types of electric discharges [2], and its mechanism depends on the nature of the discharge itself. Electric wind has many growing industrial and research applications, such as heat transfer enhancement [3, 4], EHD pumps and micro-pumps [5, 6], EHD flow control [7] and ionic loudspeakers [8]. Moreover, practical devices exploiting electric wind usually have additional advantages over other technologies, including a simple design, small weight and size, no need for moving parts, low cost, and a long lifetime. In the last decades many theoretical and experimental efforts have been conducted to increase our understanding on the EHD flow produced by corona discharges. For instance, the pioneer model of Robinson [9] related the velocity of the ionic wind to the current intensity of the electric discharge, © 2018. This manuscript version is made available under the CC-BY-NC-ND 4.0 license By No no https://creativecommons.org/licenses/by-nc-nd/4.0/ You can view the Version of Record for this article at: https://doi.org/10.1016/j.elstat.2018.02.001 2 which was confirmed in his experiments. Béquin et al. [10] investigated the velocity flow distribution in negative point-to-plane corona in air. They developed a one-dimensional model of neutral particle velocity along the corona discharge axis, and they measured the velocity profiles using Laser Doppler Anemometry (LDA). The discrepancy between the modeling and the measurements were attributed to the omission of the ring vortex formation in the model, and to the perturbation that tracing particles may have on the corona discharge. For the case of positive DC corona discharge, Lacoste et al. [11] were able to estimate the gas velocity using a simplified model, and they compared the model predictions with the measurements obtained using LDA. The measured velocity profile was reasonably fitted by the model, but the radius of the corona plasma region was an adjustable parameter in their model. More recently, a wide variety of numerical techniques and modeling approaches have been applied to investigate and characterize the generation of ion wind in electric coronas. For example, Zhao and Adamiak [12, 13] used a hybrid numerical algorithm to solve both the Gauss and charge transport equations. Once the EHD force was determined, Navier-Stokes equations were integrated using the finite-volume method in order to evaluate the gas flow. The authors followed an iterative procedure to ensure the consistency between the corona simulation and the fluid flow programs. Ahmedou and Havet [14] have simulated the airflow in a flat duct, in which corona discharge is generated by a single or multiple fine-wire electrodes. The electrical problem (Gauss and charge conservation equations) was coupled to the fluid dynamics problem (mass, momentum and energy conservation equations) and they were solved in the drift region using the finite element method. Since the ionization layer at the vicinity of the wire was not modeled, the authors have to formulate a proper boundary condition for the space charge density on the corona wire. In order to reduce the computational costs, Seimandi et al. [15] have proposed an asymptotic model for steady wire-to-wire corona discharges, in which the inter-electrode space is divided in three regions: two thin ionization regions around the electrodes and one larger ion-drift region. A simplified kinetic model was considered for each region, which allowed the establishment of quasi-analytical solutions for the electric field and the particle fluxes, and to estimate the velocity of the ionic wind from Euler equations. Cagnoni et al. [16] have carried out a numerical study with the aim of modeling the EHD cooling of a condensation radiator. They have applied a staggered solution algorithm to solve the strongly coupled set of partial differential equations that governs the problem. Very recently, Chen et al. [17] have studied numerically and experimentally the ionic wind generation in a point-to-cylinder negative corona discharge operating in the Trichel regime. They observed that flow velocity and the EHD body force have opposite directions in the ionization region, close to the tip, and in the ion drift region, further away from the tip. In the present work, the generation of EHD wind will be modeled for a parallel wire-to-plate electrode geometry, both for the positive and negative polarities. Corona discharges are time-dependent phenomena, and a complete description of the problem would require solving a fully dynamic numerical model, which demands very long computational times. Therefore, in this study, a stationary or steady-state model will be used to formulate the governing equations of the corona discharge for electrons, positive ions and negative ions. This type of model allows the evaluation of the spatial distribution of time-averaged corona parameters with less numerical effort, and they have been widely used in the scientific literature [15, 18, 19]. The governing equations will be solved semi-analytically, so that simple relations of the time-averaged spatial distributions of ions, electrons and their fluxes are obtained. These approximated solutions are useful in modeling studies and applications of coronas and, in particular, to evaluate the electrical force that originates the EHD motion. Thus, in a second step, the electrical force will be used in the Navier-Stokes equation to compute the 2D spatial distribution of the gas velocity. Special attention will be paid in this study to the ion wind structure in the vicinity of the corona wire for each polarity, as well as to the influence of the geometrical 3 parameters (wire radius and inter-electrode separation) and electrical parameters (applied voltage) on the ion wind intensity. The present study is based in previous models of corona discharge, which have been applied to different electrode geometries [18, 20-22]. These models have been validated by comparing their results with other numerical techniques [20, 22] and, indirectly, with the experimental measurements of the ozone produced by the corona discharge [21]. This work also follows the same strategy as that adopted in [23], where the electric wind generated in positive corona discharge was evaluated. However, the emphasis of that work was placed on the effect of the EHD motion on the spatial distribution of the chemical species generated by the corona discharge. In contrast, the goal of the present study is to investigate the differences in the ion wind produced by positive and negative coronas, and how it is affected by geometrical and electrical parameters. 2. Governing equations The correct assessment of the charged particles densities and of the electric field is essential to have a good approximation of the EHD force acting on the fluid. This force is an input parameter in the Navier-Stokes equations, which ultimately will determine the gas flow profile. Therefore, in this section, the mathematical model of the positive and negative corona discharge will be firstly presented in section 2.1. Then, the expression for the EHD force will be defined in section 2.2. Finally, in section 2.3, the mechanical equations governing the fluid flow will be introduced. 2.1. The corona discharge model Corona discharge is assumed to take place in air at atmospheric pressure and room temperature using a parallel wire-to-plate electrode configuration. DC high voltage is applied to the wire, while the plate is grounded. The radius of the wire and the distance between the wire and the plate will be designated as r0 and d, respectively. The electrical discharge is supposed to have translational symmetry in the direction of the wire, so that it can be modeled as a plane 2D problem. Furthermore, it is assumed that the electrical discharge is symmetric with respect to the symmetry axis of the electrode configuration. These assumptions are quite reasonable when modeling positive coronas. However, negative coronas may appear either as a glow, with a very uniform distribution of current, or concentrated into small 24, 25]. The first mode, termed Vann Bush and Snyder, is favored by the used of smoothly polished wire. Therefore, the previous assumptions fit well into this first mode of negative corona. Corona discharge can be successfully simulated using a fluid approximation, which consists in solving a set of continuity equations for the charged particles coupled with the Gauss equation for the electric field. Usually, only three generic types of charged particles (electrons, positive ions, and negative ions) need to be considered for the physical modelling [26, 27]. In the case of a positive corona, negative ions and electrons spread over a short distance (~10 m) from the wire [22], but the negative ion density is much weaker than the electron density. Therefore, for the purpose of this work, negative ions can be ignored in the positive corona. Assuming a stationary discharge, the governing equations can be written as , e e J J (1) , p e J J (2) 4 , n e J J (3) 0 0 , p e n eN N NE (4) where subscripts e, p and n correspond to electrons, positive ions and negative ions; Ji and Ni denotes the flux and the number density of particles of type i (i = e, p and n), respectively; E is the electric field, and are the ionization and attachment coefficients, respectively; 0 is the air permittivity; and e0 is the elementary charge. The flux of each type is given as Ji = i Ni E, where i is the electrical mobility of particle i. As shown in previous works, bipolar coordinates [28, 29] are particularly adequate to formulate electrostatic problems in wire-plate electrode geometry, since the coordinate curves = const. coincides with the Laplacian field lines (figure 1). Therefore, in order to solve semi-analytically (1)- (4), the same strategy as that adopted in [21] for the positive corona will be used here, namely, the electric field lines in the corona discharge will be approximated by the Laplacian field lines. The justification of this approximation lies in fact that field lines computed with more complex models do not differ significantly from Laplacian field lines [30]. The corona discharge model will also be extended to include the case of the negative polarity, and some additional approximations will be introduced later in the vicinity of the wire, which will greatly simplify the solution. Taking into account the above considerations, the set of equations (1)-(4) can be written using the bipolar coordinates system ( , ) as: 1 e e dSJ J S dl , (5) 1 p e dSJ J S dl , (6) 1 n e dSJ J S dl , (7) 0 0 1 p e n d e SE N N N S dl , (8) Figure 1. Bipolar coordinates system ( , ) used to describe the corona discharge in the wire-to-plate electrode geometry. Coordinate lines are shown as dashed lines. (Not to scale). 5 where the upper (lower) sign is for positive (negative) corona, E is the magnitude of the electric field; Ji = i Ni E is the magnitude of the flux of i-th particles; dl Sd is the infinitesimal displacement along the electric field line ( = const.) in the direction of the plate, and S is the scale factor of coordinates and . The equations defining the scale factor and the correspondence between Cartesian and bipolar coordinates are sin sinh , , , cosh cos cosh cos cosh cos a a a S x y (9) where distance a is related to the wire radius, r0, and the electrode separation, d, as 2 2 0 0 ( )a d r r . Note that a d for a thin wire ( 0 r d ). The ionization and attachment coefficients to be used in (5)-(7) have been fitted from the data provided by Bolsig [31, 32] for air as follows, 4 5 1 1 exp , 1.1 10 2.0 10 V/cm B A E E (10) 2 3 4 5 1 2 3 4 , 1.1 10 2.0 10 V/cmE E E E (11) where A1 = 6.4×103 cm 1, B1 = 1.9×105 Vcm 1, 1 = 13.60 cm 1, 2 = 8.20×10 5 V 1, 3 = 7.78×10 12 V2 cm and 4 = 1.042 ×10 15 V 3 cm2. From the continuity equations (5)-(7), it can be shown that the product JS, where e p n J J J J , is a constant along the field lines, const.JS (12) Therefore, the current density measured at the wire, j0 = e0 J0, and at the plate, jT = e0 JT, can be linked to each other as 0 0 0 cosh cos 1 cos T T T S j j j S , (13) where subscripts 0 and T denote that the corresponding quantities are evaluated on the wire and the plate, respectively. The corona current distribution over the ground plate, jT, is formulated from law [33, 34] as 0 2 2 2 0 1 cos (0)cos (0) sin 1 cos m m T d r j j j a d r , (14) where m = 4.5 (m = 4.65) for positive (negative) corona, and j(0) is the current density at the center of the plate (on the axis of symmetry). The value j(0) can be determined from the total current I as follows, /2 2 0 0 0 2 2( ) (0) cos (0) m T I j wdx d r w j d d wj K , (15) where w is the length of the wire and K = 1.44 (K = 1.40) for positive (negative) corona. Following the approach of Chen and Davison [18], the space between the electrodes can be split into two different zones: the corona plasma region and the drift region. Ionization, electron attachment, and electron-induced chemical reactions take place within the corona plasma region, in the vicinity of the wire. The drift region extends from the outer boundary of the plasma region up to the ground 6 electrode. No further plasma activity occurs in drift region, since the number density of electrons is small. Therefore, in this zone, ions injected from the corona plasma region drift towards the collecting plate under the effect of the electric field. The boundary between the plasma region and the drift region must be situated at a distance from the wire where the density of electrons is sufficiently weak. Since the wire diameter is very small compared with the electrode separation, bipolar coordinates can be approximated by polar coordinates inside the corona plasma region. Therefore, the scale factor S becomes S r(S0/r0), with r the radial distance to the centre of the wire, and the infinitesimal displacement along the electric field line satisfies dl dr. This approximation will greatly simplify the solution of (5)-(8) in the plasma region. 2.1.1 Electric field As shown in previous studies [22], most of the space charge accumulates in the drift region. Therefore, the time-averaged electric field distribution in the plasma region can be obtained by integrating (8) without space charge, which gives 0 0 0 0 E S E r E S r , (16) where subscript 0 refers to the wire. In the drift region, the space charge is mainly constituted by ions with the same polarity as that of the corona wire, that is, /( ) /( ) i i i T T i N J E J S ES , where i = p in the positive corona, and i = n and in the negative corona. Therefore, integration of (8) gives [21], 0 1 2 3 2 21 0 0 1 /sin 2 1 cos sinh , tan tanh ( , ) 1 cos tan sin 2 cosh cos i a E S E c S , (17) where ci = 2jT /( 0 i) and E0 is the electric field on the wire surface, which can be estimated using 35]. This expression of the electric field can actually be used in whole inter-electrode space, since the contribution of the second addend inside the square root becomes negligible in the corona plasma region, where 0, and (17) becomes (16). 2.1.2 Electron flux Since electrons are confined to the plasma region, bipolar coordinates can be approximated by polar coordinates, and the ionization and attachment coefficients can be written as 1 2 expA B r , (18) 2 3 4 1 2 3 a a a ar r r , (19) where 2 1 0 0 /( )B B E r , a1 = 1, 2 2 0 0 ( ) a E r , 2 3 3 0 0 ( ) a E r and 3 4 4 0 0 ( ) a E r . Then, the flux of electrons can be obtained by direct integration of (5) 0 0 0 0 0 0 0 exp exp ( ) l r e e e e r S J S J J dl dr J f r S S , (20) where 2 1 0 1 4 2 2 0 1 0 3 2 2 2 0 0 1 1 1 1 exp exp( ) exp( ) 2 a r A a f r B r B r a r r a r B r r r r (21) 7 and Je0 is the electron flux on the wire surface. As usual, the upper (lower) sign corresponds to positive (negative) corona. 2.1.3 Positive ion flux Within the plasma region, the flux of positive ions can be obtained from (6) as follows 0 0 0 0 0 1 2 0 0 1exp ( ) l r p p e p e r S r J J S J dl J A r B r J f r dr S S r , (22) where Jp0 is the flux of positive ions on the wire. In the positive corona, it must be Jp0 = 0, owing to the polarity of the wire. In the negative corona, the current on the wire is contributed by positive ions and electrons, but 0 0 p e J J . Therefore J0 Jp0, and flux of positive ions on the wire can be linked to the current on the plate using (13) 0 0 0 1 T p T S J j e S . (23) Equation (22) is particularly useful in the case of a negative corona, since positive ions are confined to the plasma region. However, in a positive corona, positive ions accumulate in the drift region. In such a case, a more straightforward evaluation of the positive ion flux can be obtained using (12) T T p e S J J J S , (24) which further simplifies in the drift region to / p T T J S J S , since p e J J . 2.1.4 Negative ion flux As stated previously, in the positive corona, the number density of negative ions is much weaker than that of electrons, and they are confined to a short distance from the wire. Therefore, they will not contribute significantly to the electric force and will be omitted in the present study. In contrast, in the negative corona, negative ions are present both in the corona plasma region and in the drift region. In the corona plasma region, the flux of negative ions can be deduced from (7) as 0 0 2 3 4 0 1 0 2 3 0 1 1 ( ) l r n n e e r S a a a J J S J dl r a J f r dr S S r r r r , (25) where it has been taken into account that Jn0 = 0, due to the polarity of the wire. However, once the electron flux and the positive ion flux have been previously determined, a more convenient relation can be obtained using (12), T T n p e J S J J J S , (26) which is valid in the whole discharge gap. In the drift region, where n e J J and n p J J , the above equation can be approximated as / n T T J S J S . 2.1.5 Electron flux on the corona wire The value of the electron flux on the corona wire, Je0, is required in (20), (22) and (25) in order to evaluate Je, Jp and Jn, respectively. 8 In the positive corona, the flux of positive ions on the wire must be null (Jp0 = 0) and, as explained before, the flux of negative ions can be ignored compared to the flux of electrons. Therefore, the electron flux on the wire surface can be linked to the total current density on the plate using (13) 0 0 0 1 T e T S J j e S . (27) However, in the negative corona, positive ions constitute the major contributors to the current density on the wire. Thus, a different strategy must be followed in order to determine the electron flux on the corona electrode. The electron current reaches its maximum value in the corona plasma region when the ionization coefficient equals the attachment coefficient, = , which corresponds to an electric field strength EL 30 kV/cm and to a radial distance 0 0 / L L r E r E . At this location, the positive and negative ion currents are very similar, so that ( ) ( ) n L p L J r J r , and the total current density can be linked to that on the plate using (12) 0 0 0 ( ) ( ) ( ) ( ) 2 ( ) T T e L p L n L L e L p L L e J S e J r J r J r S e J r J r S . (28) This equation can be used as condition to determine the electron density on the wire surface. Substitution of (20) and (22) in (28) leads to 0 0 0 0 1 2 / 2 exp ( ) ( ) L T T er L L r r J S S J A r B r f r dr r f r . (29) 2.2. The EHD force expression Charged particles generated in the corona discharge are accelerated by the electric field and gain kinetic energy. Moreover, they transfer their momentum to the surrounding neutral gas molecules through collisions, which results in a 36], and assuming that the air motion has a negligible influence on the corona discharge, the EHD force density acting on the fluid can be expressed as 0 p n e e N N N F E . (30) Since Ji = i Ni E (i = e, p and n), the electric force can also be written as a function of the flux of charged species, 0 p n e p n e eJ J J F . (31) 2.3. Gas flow model The fluid flow can be described by Navier-Stokes equation and the continuity equation, which represent the conservation of momentum and mass, respectively. Assuming that the air motion can be approximated as an incompressible Newtonian flow, these equations can be written as ,p t V F V V V (32) 0, V (33) 9 where V is the gas velocity, is the kinematic air viscosity, p is the gas pressure, and is the air density. In this work, the flow is assumed to be confined between two parallel plates, with the bottom plate (the ground electrode) located at y = 0 and the top plate at y = 2d + 2r0 (see figure 2). The channel width is taken as 10 cm, the entrance being situated at x = 5 cm and the exit at x = 5 cm. The length of the channel, along the z-direction, is assumed to be infinite. The corona wire runs parallel to the z-axis, with its center at the coordinates x = 0 and y = d + r0. Boundary conditions for the airflow are straightforward. Firstly, no-slip conditions are imposed at the plates and at the wire surface. At the entrance of the channel (x = 5 cm), the gas flow is assumed to be parallel to the x-direction, with a constant velocity V0. At the exit (x = 5 cm), the gas pressure is imposed to be equal to the atmospheric gas pressure. The electrohydrodynamic motion induced by the corona discharge will be, in general, turbulent [12, 14]. Therefore, the usual approach of time-averaging the conservation equations of momentum and mass, which leads to the so called Reynolds-Averaged Navier Stokes (RANS) equations, has been followed here. The resulting system of partial differential equations has then been closed using the standard ktwo-equation model of turbulence [37]. Thus, the turbulent viscosity, t, the turbulent kinetic energy, k, and its rate of dissipation, , are linked through the relation t = Cµ2Cµ = 0.09, Prandtl numbers k and have been given the values k = 1.00 and = 1.30, and constants C1 and C2 of the model are set as C1 = 1.44 and C2 = 1.92 [37]. The resulting system of equations has been solved for the stationary regime with the help of the computational fluid dynamics software FLUENT [38], which uses the finite volume method. For that purpose, the pressure-based solver SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) has been applied. 3. Results 3.1 EHD force and ion wind velocity The numerical results presented in this section have been obtained for a corona wire of radius r0 = 100 µm and an inter-electrode separation d = 0.5 cm. The positive and negative high voltage applied to the corona wire have been chosen as = ±8 kV, to which corresponds a current intensity per unit of wire length of 30.11 µA/cm and 37.6 µA/cm, respectively [39, 30]. As we are mainly interested in studying the effect of the EHD contribution to the air motion, the gas velocity at the entrance of the channel was assigned a small value, V0 = 2.5 cm/s. In any case, the fluid flow within the channel will be turbulent, because the ion wind velocity generated by the discharge of the corona is substantially higher than V0. Figure 2 . Schematic picture of the wire plate corona reactor . 16 is favored, and they have a much larger size and velocity than in the positive corona. The maximum velocity, V cm/s, is reached on the symmetry axis. Although this velocity is smaller than the ones observed in the drift region, the corresponding velocity gradients are larger. Figure 10 highlights the differences in the EHD induced velocity along the symmetry axis between the positive and a negative corona discharge. The air velocity is zero at the corona wire and at the ground plate, due to the non-slip condition. However, the gas velocity increases faster and reaches a higher speed (Vmax 3.8 m/s) in the positive corona than in the negative corona (Vmax 2.3 m/s). Also, in the positive corona, the point where the maximum speed is attaint is closer to the wire. The prediction of the model of a higher gas velocity in the positive corona is in agreement with the experimental observations of Kozlov and Solovyov [41] and other investigators. As explained by these authors, the ionic wind velocity depends on the ion mobility (the larger the ion mobility the lower the ionic wind velocity) and the mobility of positive ions is lower than that of negative ones ( p = 2×10 m2V s and n = 2.5×10 m2V s in this work). However, according to the results of the present study, the structure of the corona plasma region is another important factor that needs to be taken into account: in the negative corona, the EHD force changes its direction from the ionization to the drift region and, therefore, the fluid is accelerated towards the plate less efficiently than in the case of a positive corona. Comparison with experiments In this subsection, modeling results are compared with the experimental measurements reported by Elagin et al. [4] for a parallel wire-plate electrode configuration. They used the particle image velocimetry (PIV) method to visualize the flow, and the purpose of their study was to investigate the cooling enhancement produced by the ionic wind. Therefore, the plate was set horizontally below the wire, and it was heated to a steady temperature of about 65 K above the ambient temperature. However, since the gas velocity due to natural convection was very small compared to that produced by the ion wind, the velocity pattern and velocity values measured by the PIV technique are essentially due to EHD gas motion. Furthermore, in their experiments, negative corona was characterized by the is not compatible with the assumptions of the model. As a result, ionic wind jets emitted from these localized spots resulted in an unstable and chaotic flow of difficult study. Therefore, comparison with experiments will be here limited to the case of positive corona. In order to match the experimental set-up used by Elagin et al., the top boundary (upper plate) in our modeling configuration (figure 2) has been suppressed. Therefore, excluding the lower plate and the wire, where a non-slipping condition is imposed, the rest of boundaries are treated as open boundaries, where pressure is assumed to be equal to the atmospheric pressure. However, these modifications do not affect the validity of the electrohydrodynamic force density previously obtained. Figure 11 shows the numerical simulation of the flow pattern corresponding to the same experimental conditions as in one of the two experiments reported by Elagin et al., using a 11 cm long corona wire of radius 50 µm, and for an applied voltage of 13.5 kV. The general structure of the ionic wind velocity distribution matches closely the observations of Elagin et al., with an acceleration zone near the wire, a uniform flow zone in the central jet, and the formation of a stagnation zone when the jet collides with the plane. Then, the ionic wind jet flows sideways along the plate in opposite directions with a high velocity. A quantitative comparison between modeling and simulation is shown in Table 1, where the maximum the ion wind velocity is given for two different sets of electrical and geometrical parameters. 17 Table 1. Maximum ion wind velocity magnitude in wire-plate positive corona discharge for a wire of radius r0 = 50 µm and length w = 11 cm. Comparison between numerical modeling and experiments. Vmax (m/s) d (cm) (kV) I (µA) numerical experimental [4] 1.5 13.5 230 3.65 3.5 4 26 110 2.74 2.75 3.2 Parametric Study In this section, the influence of different parameters on the electric wind velocity induced by the corona discharge will be investigated, both in the positive and in the negative polarity. For simplicity, the gas velocities will be only shown along the axis of symmetry. Three different parameters will be here considered: the applied voltage, the wire radius and the inter-electrode separation. Effect of the applied voltage Figure 12 shows the electric wind velocity along the y-axis corresponding to a positive (figure 12a) and a negative (figure 12b) corona discharge for three different values of the applied voltage: ±7 kV, ±9 kV, and ±11kV. The current intensities per unit of length corresponding to these voltages are, in the positive corona, 15.8 µA/cm, 47.4 µA/cm and 91.1 µA/cm, respectively, and, in the negative corona, µ µ µA/cm, respectively. As stated by (34), the raise of the applied voltage and, therefore, of the corona current intensity, implies an augmentation of the EHD force. Thus, an increase of the fluid velocity may be expected, as confirmed by the results presented in figure 12. Figure 11. 2D- spatial distribution of the gas velocity magnitude corresponding to a positive corona discharge between a wire and a plate in an open channel. Streamlines are drawn as white solid lines. Wire radius: r0 = 50 µm, wire length: w = 11 cm, wire-plate distance: d = 1.5 cm, applied voltage: = +13.5 kV, current intensity: I = 230 µA. 18 In the positive corona (figure 12a), the flow velocity increases very rapidly beyond 100 µm from the wire: the higher the applied voltage, the faster the velocity gradient. However, in the vicinity of the wire, the velocity of the recirculating vortices also increases with the applied voltage. This augmentation is due to the shear with the adjacent mainstream flow, whose velocity has grown more rapidly with voltage. As previously said, in the positive corona, these vortices have mainly a hydrodynamic origin. This is confirmed by the fact that the length of the recirculating region increases as Reynolds number, Re, increases [40] (Re = Ud/ , where U is the main stream flow velocity, d is the wire diameter and is the kinematic viscosity). At lower voltages, the recirculating flow tends to disappear. In contrast, in the negative polarity, the air velocity is directed towards the wire along the first 300 µm, and towards the plate in the rest of the space. The velocity in the corona plasma region increases faster with the applied voltage than the drift region: the ratio of the maximum velocities in the plasma region and in the drift region is about 11% at 7 kV, and increases up to 21% at 9 kV. Finally, contrary to the positive corona, the coordinate at which the fluid velocity changes its sign is rather insensitive to the applied voltage, due to the active role played by the electric force in sustaining the vortices. Different published data suggest an almost linear dependence of the averaged wind velocity with the square root of the discharge current [9, 41]. These observations are coincident with the results obtained from the present model since, as shown in the insets of figure 12, the maximum velocity scales as the square root of the current intensity, both in the positive and in the negative corona. Effect of the wire radius The effect of varying the wire radius on the EHD velocity is illustrated in figure 13 for an applied voltage of ±7 kV and an inter-electrode separation of 0.5 cm. Three values of the wire radius have been considered: 50 µm, 100 µm, and 150 µm. 35], the radius of the corona wire affects both the electric field on the cathode, 1/2 0 0 29.8 (1 0.301/ ) E r kV/cm, where r0 is in cm, and the critical voltage for the onset of corona discharge, 0 0 0 ln(2 / ) c E r d r . Decreasing the radius of the corona wire lowers the corona Figure 12. Gas velocity magnitude along the axis of symmetry (x = 0) in (a) positive corona and (b) negative corona, for three different values of the applied voltage: = ±7, ±9 and ±11 kV. Wire radius: r0 = 100 µm, inter - electrode separation: d = 5 mm. 19 onset voltage and more current will flow for the same applied voltage. In particular, the corona currents per unit of length corresponding to the radii mentioned above are, for the positive polarity, 26.1 µA/cm, 15.8 µA/cm and 5.89 µA/cm, respectively, and, for the negative polarity, µA/cm, respectively. Consequently, the smaller the wire radius, the higher the current intensity, the stronger the EHD force, and the faster the gas velocity. This correlation is readily observed on the results presented in figure 13: the air velocity in the drift region more than doubles when the radius is reduced from 150 µm to 50 µm, both in the positive and in the negative corona. Moreover, as in the previous case, the peak values of the flow velocity increases linearly with the square root of the current intensity (see insets in figure 13). The variation of the wire radius also affects the flow velocity in the close vicinity of the wire. Since Reynolds number is proportional to the diameter of the wire, increasing the wire radius favors a longer recirculation region below the wire. In the positive corona (figure 13a), for the chosen applied voltage, the recirculation region is very weak, but this effect can be appreciate in the fact that the gas velocity directed towards the plate starts growing at a longer distance from the wire surface when the radius of the wire is increased. In the negative corona, an additional factor comes into play. According to (19), the outer boundary of the ionization layer is located around 0 0 / L L r r E E , where EL is the electric field for which the ionization and attachment coefficients become equal. Thus, the extent of the ionization region, 0 0 ( / 1) L r E E , is enlarged when the radius of the corona wire becomes thicker and, correspondingly, there is an increase of the area where the electric force and the air motion is directed towards the wire (see figure 13b). The maximum fluid velocities in this region are similar for 50 µm and 100 µm, but larger than for 150 µm. Effect of the inter-electrode separation Figure 14 shows the effect of increasing the inter-electrode gap on the air velocity along the y-axis. Three values of the inter-electrode separations (5 mm, 7.5 mm and 1 cm) have been considered, for an applied voltage of ±7 kV and wire radius of 100 µm. Figure 13. Gas velocity magnitude along the axis of symmetry (x = 0) in (a) positive corona, = 7 kV, and (b) negative corona, = 7 kV, for three different values of the wire radius: r0 = 50 m, 100 m and 150 m. Inter - electrode separation: d = 5 mm. 20 As shown in the previous section, the critical voltage for the onset of corona discharge, 0 0 0 ln(2 / ) c E r d r , is also affected by the electrode separation. As the distance between the wire and the plate is augmented, the corona onset voltage increases and, therefore, less current will flow for the same applied voltage. In the present case, the corona currents per unit of length corresponding to the inter-electrode gaps mentioned above are, for the positive corona, 15.8 µA/cm, 4.36 µA/cm and 1.53 µ A/cm, respectively, and, for the negative corona, .7 µA/cm and µA/cm, respectively. The reduction of current intensity means a lower EHD force and a decrease of the flow intensity in the drift region, which can be easily appreciated in figure 14. In the positive corona, the fluid velocity in the drift region decreases by nearly a factor 3 when the interelectrode gap is doubled. In the negative corona, the corresponding reduction is of about a factor of 2. In both polarities, the point where the gas velocity reaches its maximum value moves away from the wire as the electrode separation is increased. As in the previous cases, the maximum flow velocity scales as the square root of the current intensity (see insets in figure 14). In the negative corona, the secondary vortices beneath the wire are importantly affected by the separation of the electrodes, though [35], the electric field on the wire surface is unaffected by the distance between the electrodes. However, since the current intensity decreases when the inter-electrode separation is raised, the density of electrons and ions in the ionization region must lower, and correspondingly the EHD force directed towards the wire. This fact, added to the lower velocity of the main stream flow around the wire, which favors a shorter recirculation region, leads to an important reduction of the extension and intensity of the velocity in this region. Therefore, as shown in figure 14b, the recirculation zone is progressively reduced and becomes negligible for an inter-electrode distance of 1cm. 4. Conclusions A simplified formulation of the positive and negative corona discharge in air has been proposed to compute the gas flow generated by the EHD force. The approach assumes that the discharge gap is divided in two zones: the corona plasma region and the ion drift region. The first zone is a thin layer close to the corona electrode where ionization and electron attachment take place. The second zone corresponds to the rest of inter-electrode space, where ions with the same polarity as that of the corona Figure 14. Gas velocity magnitude along the axis of symmetry (x = 0) in (a) positive corona, = 7 kV, and (b) negative corona, = 7 kV, for three different values of the inter-electrode separation: d = 0.5, 0.75 and 1 cm . Wire radius: r0 = 100 m. 21 electrode carry the electric current. Semi-analytical expressions for the electric field and the current densities of electrons, positive ions and negative ions are obtained coupled with the continuity equations of charged particles. These approximated solutions can be evaluated with a minimal computational cost, which makes them interesting for the analysis and optimization of devices that employ corona discharge in different technological applications. Here, they have been used to obtain the expression of the EHD force responsible for the gas flow, which has been inserted in Navier-Stokes equations in order to compute the flow distribution. The results of the simulation have shown that the structure of the plasma region of the corona discharge has a significant effect on the EHD gas flow. In the positive corona, the negative space charge is confined to a layer of negligible extent around the corona wire. Therefore, the electric force is directed towards the plate in the whole inter-electrode gap. In spite of this, two tiny vortices are observed beneath the wire, with recirculation velocities opposite to those in the adjacent large roll. These vortices have a hydrodynamic origin, since they are produced by the downwards air flow coming from above the wire. However, the electric force tends to damps down these vortices. In contrast, in the negative corona, the EHD force is directed towards the wire in the inner part of the plasma region, where ionization prevails over attachment and positive ions accumulate. Therefore, the electrical force helps to sustain the vortices generated underneath the wire. As a consequence, the effective length along which the fluid is accelerated towards the plate is reduced, precisely in the region where the electric force is more intense. Thus, the electric wind velocity observed in the negative corona is smaller than in the positive corona, even though the negative corona current is higher that the positive corona current at a given voltage. The simulation results have been compared with experimental measurements for the case of a positive corona, and a good agreement has been found. The intensity of the electric wind is affected by the applied voltage, the wire radius and the interelectrode separation but, in all cases, the maximum wind speed in the drift region scales as the square root of the current intensity, as reported in experimental observations. In particular, at a fixed voltage, reducing the wire radius and the inter-electrode distance leads to a faster EHD flow velocity. The secondary vortices located beneath the wire are also affected by these parameters, especially by the inter-electrode separation: they progressively die down when the gap between electrodes is increased. Acknowledgments This research has been funded by the Spanish Government Agency de Economía y under contract no. FIS2014-54539-P. References [1] Cobine J D (1958) Gaseous Conductors: Theory and Engineering Applications (New York: Dover) [2] control by non-thermal plasma actuators J. Phys. D: Appl. Phys. 40 605 636 [3] Molki M and Damronglerd P (2006) Electrohydrodynamic enhancement of heat transfer for developing air flow in square ducts Heat Transfer Eng. 27 35 45 [4] Elagin I A, Yakovlev V V, Ashikhmin I A, and Stishkov Yu K (2016) Experimental Investigation of Cooling of a Plate by Ionic Wind from a Corona-Forming Wire Electrode Technical Physics 61 1214 1219 [5] Zhao L and Adamiak K (2016) EHD flow produced by electric corona discharge in gases: From fundamental studies to applications (a review) Part. Sci. Technol. 34 63-71 22 [6] Birhane Y T, Lin S-C and Lai F C (2017) Flow characteristics of a two-stage EHD gas pump in a circular pipe IEEE Trans. Ind. Appl. 53 2461-2470 [7] Léger L, Moreau E and Touchard G (2002) Effect of a DC corona electrical discharge on the airflow along a flat plate IEEE Trans. Ind. Appl. 38 1478- 1485 [8] Bastien F (1987) Acoustics and gas discharges: application to loudspeakers J. Phys. D: Appl. Phys. 20 1547 1557 [9] Robinson M (1961) Movement of air in the electric wind of the corona discharge AIEE Transactions 80 143-150 [10] Béquin Ph, Castor K and Scholten J (2003) Electric wind characterisation in negative point-to- plane corona discharges in air Eur. Phys. J. AP 22 41-49 [11] Lacoste D A, Pai D and Laux C (2004) Ion wind effect in a positive dc corona discharge in atmospheric pressure air 42nd AIAA Aerospace Sciences Meeting and Exhibit (Reno, Nevada, 5-8 January 2004) [12] Zhao L and Adamiak K (2009) Effects of EHD and external airflows on electric corona discharge in point-plane/mesh configurations IEEE Trans. Ind. Appl. 45 16-21 [13] Zhao L and Adamiak K (2013) Numerical simulation of the effect of EHD flow on corona discharge in compressed air IEEE Trans. Ind. Appl.49 298-304 [14] Ahmedou S O and Havet M (2009) Effect of process parameters on the EHD airflow J. Electrostat. 67 222 227 [15] Seimandi P, Dufour G and Rogier F (2009) An asymptotic model for steady wire-to-wire corona discharges Math Comput Model 50 574-583 [16] Cagnoni D, Agostini F, Christen T, Parolini N, Stevanovic I and de Falco C (2013) Multiphysics simulation of corona discharge induced ionic wind J. Appl. Phys. 114 233301 [17] Chen S, Nobelen J C P Y and Nijdam S (2017) A self-consistent model of ionic wind generation by negative corona discharges in air with experimental validation Plasma Sources Sci. Technol. 26 095005 [18] Chen J H and Davidson J H (2003) Model of the negative DC corona plasma: comparison to the positive dc corona plasma Plasma Chem. Plasma Process. 23 83-102 [19] Zhang J and Adamiak K (2007) A multi-species DC stationary model for negative corona discharge in oxygen; point-plane configuration J. Electrostat. 65 459 464 [20] Yanallah K and Pontiga F (2012) A semi-analytical stationary model of a point-to-plane corona discharge Plasma Sources Sci. Technol. 21 045007 [21] Yanallah K, Pontiga F and Chen J H (2013) A semi-analytical study of positive corona discharge in wire plane electrode configuration J. Phys. D: Appl. Phys.46 345202 [22] Yanallah K, Pontiga F, Meslem Y and Castellanos A (2012) An analytical approach to wire-to- cylinder corona discharge J. Electrostat. 70 374-383 [23] Yanallah K, Pontiga F, Bouazza M R and Chen J H (2017) The effect of the electric wind on the spatial distribution of chemical species in the positive corona discharge J. Phys. D: Appl. Phys 50 335203 [24] Vann Bush P and Snyder T R (1987) Laboratory analyses of corona discharges J. Electrostat. 19 83 99 [25] Chang J S, Lawless P A and Yamamoto T (1991) Corona discharge processes IEEE Trans. Plasma Sci. 19 1152 1166 [26] Morrow R (1985) Theory of negative corona in oxygen Phys. Rev. A 32 1799 1809 23 [27] Kulikovsky A A (1997) The mechanism of positive streamer acceleration and expansion in air in J. Phys. D: Appl. Phys. 30 1515 1522 [28] Morse P M and Feshbach H (1953) Methods of Theoretical Physics (New York: McGraw-Hill) [29] Spiegel M R, Lipschutz S and Liu J (1998) Mathematical Handbook of Formulas and Tables (New York: McGraw-Hill). [30] Sigmond R S (1986) The unipolar corona space charge problem, J. Electrostat. 18 249 272 [31] Hagelaar G J M BOLSIG+: Electron Boltzman Equation Solver [online] http://www.bolsig.laplace.univ-tlse.fr (accessed July 2015) [32] Hagelaar G J M and Pitchford L C (2005) Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models Plasma Sources Sci. Technol. 14 722- 33 [33] Warburg E (1927) Characteristik des spitzenstromes, Handbuch der Physik (Berlin: Springer). [34] Zebboudj Y 2000 Measurements of current and electric field distributions beneath a positive DC wire-to-plane corona using a linear biased probe IEE Proc.-Sci. Meas. Technol. 147 74 80 [35] Peek F W (1929) Dielectric Phenomena in High Voltage Engineering (New York: McGraw-Hill). [36] Boeuf J P and Pitchford L C (2005) Electrohydrodynamic force and aerodynamic flow acceleration in surface dielectric barrier discharge J. Appl. Phys. 97 103307 [37] Versteeg H K and Malalasekera W 2007 An Introduction to Computational Fluid Dynamics. The Finite Volume Method (Harlow: Pearson Education) [38] ANSYS® Academic Research, Release 15.0 http://www.ansys.com/ [39] Deutsch W (1933) Über die Dichteverteilung unipolarer Ionenströme Ann. Phys. 408 588 612. [40] Guyon E, Hulin J-P, Petit L and Mitescu C D (2015) Physical Hydrodynamics (New York: Oxford University Press). [41] Kozlov B A and Solovyov V I (2007) Electric wind in electrode systems with corona points Technical Physics 52 892 897