L Papageorgiou

A C Metaxas

G E Georghiou

Three-dimensional numerical modelling of gas discharges at atmospheric pressure incorporating photoionisation phenomena

A three-dimensional (3D) numerical model for the characterization of gas discharges in air at atmospheric pressure incorporating photoionisation through the solution of the Helmholtz equation is presented. Initially, comparisons with a two-dimensional (2D) axi-symmetric model are performed in order to assess the validity of the model. Subsequently several discharge instabilities (plasma spots and low pressure inhomogeneities) are considered in order to study their effect on streamer branching and off-axis propagation. Depending on the magnitude and position of the plasma spot, deformations and off-axis propagation of the main discharge channel were obtained. No tendency for branching in small (of the order of 0.1cm) overvolted discharge gaps was observed.

Introduction

Gas discharges at atmospheric pressure have been attracting a lot of attention recently due to the inherent advantages they offer in producing plasma, namely low cost production without the need for costly enclosures. Plasma medicine [START_REF] Fridman | Applied plasma medicine[END_REF] and plasma actuators [START_REF] Boeuf | Electrohydrodynamic force in dielectric barrier discharge plasma actuators[END_REF] are two areas where these discharges are currently being used widely. Due to their nature, the experimental characterization of such discharges is very difficult: their temporal evolution is of the order of a few nanoseconds while their spatial resolution is in the range of a few micrometers. Additionally, during their formation and propagation strong non-linear phenomena take place making their characterization even more challenging.

A common phenomenon observed in atmospheric pressure gas discharges in air is streamers; these are ionisation waves that propagate with velocities much higher than the maximum drift velocity of electrons. Their spatial structure is inherently three-dimensional; their experimental investigation is very difficult rendering numerical modelling a very important tool in determining a proper understanding of streamer physics [START_REF] Georghiou | Numerical modelling of atmospheric pressure gas discharges leading to plasma production[END_REF][START_REF] Papageorghiou | Two-dimensional modelling of a nitrogen DBD at atmospheric pressure[END_REF]. Until now, numerous papers have been published involving one-dimensional (1D) and two-dimensional (2D) streamer modelling, elucidating several characteristics of the streamer nature [START_REF] Georghiou | Effect of photoemission on the streamer development and propagation in short uniform gaps[END_REF]. Two-dimensional axi-symmetric models are predominantly employed allowing the study of some phenomena in a three-dimensional setting as well. Work based on 1.5D and 2D models [START_REF] Babaeva | Effect of inhomogeneities on streamer propagation: I. Intersection with isolated bubbles and particles[END_REF][START_REF] Ebert | The multiscale nature of streamers[END_REF][START_REF] Naidis | On streamer interaction in a pulsed positive corona discharge[END_REF] has been conducted by making several assumptions in an attempt to understand several aspects of streamer branching and streamer interaction. However, one important drawback of the 2D axi-symmetric models is that they restrict gas discharge modelling in configurations and conditions that are inherently symmetric around the axis of symmetry.

Phenomena like (a) streamer branching, (b) electrode configurations with no cylindrical symmetry, (c) interaction between streamers and (d) non homogeneous surfaces, which have inherently a threedimensional character, cannot be dealt with in a 2D axi-symmetric setting. Complete numerical characterization has to be based on fully three-dimensional models such as the ones developed by Hallac et al [START_REF] Hallac | Three-Dimensional Characterization of Short-Gap Streamers in Air at Atmospheric Pressure[END_REF], Pancheshnyi [START_REF] Pancheshnyi | Numerical simulation of filamentary discharges with parallel adaptive mesh refinement[END_REF], Luque [START_REF] Luque | Interaction of streamer discharges in air and other oxygen-nitrogen mixtures[END_REF] and Kulikovsky [START_REF] Kulikovsky | Three-dimensional simulation of a positive streamer in air near curved anode[END_REF]. The 3D models of Luque [START_REF] Luque | Interaction of streamer discharges in air and other oxygen-nitrogen mixtures[END_REF] and Kulikovsky [START_REF] Kulikovsky | Three-dimensional simulation of a positive streamer in air near curved anode[END_REF] are axi-symmetric with an expansion of the coordinate around the axis of symmetry. Questions on resolution and subsequently accuracy in the direction arise in regions away from the axis of symmetry in the radial direction. In the work of Hallac et al [START_REF] Hallac | Three-Dimensional Characterization of Short-Gap Streamers in Air at Atmospheric Pressure[END_REF] photoionisation was introduced via a truncated implicit subgridding method. Pancheshnyi et al [START_REF] Pancheshnyi | Numerical simulation of filamentary discharges with parallel adaptive mesh refinement[END_REF] presented results from an efficient parallel fully 3D streamer model but without the incorporation of photoionisation. Here we make a contribution by fully incorporating photoionisation in order to provide the whole picture during the formation and propagation of the streamer. Photoionisation phenomena are taken into account by solving a set of Helmholtz equations as suggested by Luque et al [START_REF] Luque | Photoionization in negative streamers: Fast computations and two propagation modes[END_REF] and Bourdon et al [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF] (this is in contrast to the computationally demanding integral method developed by Zeleznyak et al [START_REF] Zheleznyak | Photoionisation of nitrogen and oxygen mixtures by radiation from gas discharges[END_REF], introduced in 2D models by Kulikovsky [16]). The approach proposed by Luque et al [START_REF] Luque | Photoionization in negative streamers: Fast computations and two propagation modes[END_REF] is mainly based on fitting the absorption function of photo-ionization radiation with a sum of exponential functions leading to a set of integrals, each of which can effectively be interpreted as an integral solution of a separate Helmholtz differential equation. After this equivalent representation is established the integral problem of estimating the photoionization can be obtained solved by solving the set of Helmholtz differential equations, instead of direct directly evaluating the of integrals [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF]. In our case, photoionization occurs throughout absorption by oxygen molecules, of the radiation emitted by nitrogen molecules in singlet excited states b 1 , b' 1 u + and c 4 ' 1 u + in the wavelength range 98-102.5 nm [START_REF] Naidis | On photoionization produced by discharges in air[END_REF].

In this work a preliminary examination on how several discharge instabilities (plasma spots and low pressure inhomogeneities) could affect streamer branching and off axis propagation is conducted. For such studies a fully 3D numerical scheme is required. First we validate our new 3D numerical scheme with a 2D axi-symmetric one by considering a situation were both schemes are equally valid. After the validation phase, the three dimensional model can be used with confidence to study the problems stated above. Plasma spots can arise by irradiation of the gap with UV laser pulses; for example, whereas low pressure inhomogeneities represent fluctuations in the pressure that could arise in the gap. The study of such parameters is very instructive in the further understanding of streamer branching and gas discharge development. The results demonstrate clearly that, depending on the magnitude and position of the plasma spot, streamer head deformation and off axis propagation of the discharge channel are observed. On the other hand no tendency for branching in small (of the order 0.1cm) overvolted discharge gaps has been observed.

Governing equations

The differential equations that describe the cold plasma are the Poisson equation for the electric field, the continuity equations for charged particles and the Helmholtz equation for the photoionisation phenomenon. These are shown below:

(

(3)

(4) (5) (6) (7)
(8) ions respectively and e D is the electron diffusion coefficient. The symbols , , ep and np denote the ionisation, attachment, electron-positive-ion recombination and negative-ion-positive-ion recombination coefficients respectively. All of these quantities, with the exception of recombination coefficients, need to be known as a function of the electric field. In the present model the negative ions are considered to have negligible influence on the timescales considered for this problem. Φ is the electric potential, j j n q the densities of the charged species and s r the surface charge. S ph is the source term due to photoionisation, 2 O p is the partial pressure of oxygen and I is the emission of photons in the discharge volume assumed to be proportional to the number of ionizing collisions which are given by the product e e a n W ⋅ ⋅ r . and are fitting parameters [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF]. The parameters for the two-exponential fit are presented in the table 1.

0 j j S j N q r -∇ ⋅ ∇Φ - + = ε ( ) ( ) ( )
(
S N W N W N N S S N W N N N N S S N W N N α η β α β β η β = - - + = - - + = - r r r r 2

Table 1: Parameters of the two exponential fit

A (cm -2 Torr -2 ) (cm -1 Torr 1 ) 1 0.0021 0.0974 2 0.1775 0.5877

The computational cost for 3D gas discharge modelling is very high. This means that special attention should be paid in choosing the mathematical relations that yield the transport coefficients. These relations have to be as mathematically simple as possible without loss of accuracy. For this reason, the transport coefficients by Kang et al [START_REF] Kang | Numerical study on influences of barrier arrangements on dielectric barrier discharge characteristics[END_REF] and the drift velocity by Morrow and Lowke [START_REF] Morrow | Streamer propagation in air[END_REF] have been used.

The boundary conditions of the 3D model are shown below (see Table 2): At the anode the electrons are convected freely through the metallic electrode and the ions show the same behaviour at the cathode. Secondary electron emission from the cathode due to ion bombardment is taken into account by setting the gamma coefficient equal to 0.01. For the Helmholtz equations ( 8) Dirichlet boundary conditions equal to zero are introduced at the boundaries of the computational region.

Results

In this paper a plane to plane electrode configuration with a gap size of 0.1 cm is considered. An initial Gaussian electron and positive ion density (plasma spot 3 (

) n cm -given by the expression below) with maximum value A= 10 10 cm -3 and width = 0.005 cm [START_REF] Hallac | Three-Dimensional Characterization of Short-Gap Streamers in Air at Atmospheric Pressure[END_REF] are introduced in the gap at the cathode and a voltage of 5000 V is applied in the inter-electrode space.

2 2 2 0 0 0 2 2 2 0 0 0 ( ) ( ) ( ) (2 ) (2 )
(2 )

x x y y z z n A e σ σ σ -- -- -- + + ⋅ ⋅ ⋅ = ⋅ (9)
An unstructured mesh, as shown in Figure 1, was employed and optimized in order to reduce as much as possible the computational cost. The mesh is gradually becoming finer following an onion-like pattern as approaching the axis of the discharge. For the Poisson and continuity equations Lagrangequadratic elements were used while for the Helmholtz equations Lagrange linear elements were employed. The solver of the model is based on the Finite Element Method and uses Newton's iterations for the solution of the non-linear system of equations. With that approach we demand the whole system of equations to converge simultaneously. The time stepping algorithm is based on a fifth order Backward Differentiation Formula (BDF) which is an implicit method for the numerical integration the differential equations. The adopted numerical techniques allow the use of relatively large time-steps not restricted by the Courant-Friedrichs-Lewy (CFL) condition. The CFL condition slows down most of the high resolution schemes which are used in streamer modelling for the solution of the hyperbolic continuity equations. The evolution of the discharge reveals two well known phases: the avalanche and the streamer propagation phase.

Streamer evolution

The initial electrons under the influence of the electric field accelerate and multiply via ionisation processes as presented in Figure 2. The deformation of the initial Gaussian shape is apparent close to the anode where space charge effects become important. In Figure 3, the electron density close to the anode is depicted. The velocity of the avalanche is about 2.3x10 7 cm/s. This value is in accordance with results from 2D models, revealing that no phase error is introduced in the results during this part of the discharge.

The streamer formation and propagation phase is much more computationally demanding than the avalanche one. The density gradients at the streamer head are very steep and the phenomenon advances with higher velocities. For this reason a moving window of high mesh resolution is set to follow the streamer head movement. The size of the mesh that gives accurate results is of the order of 5 µm. In Figure 4 the evolution of the electron density on the axis of symmetry of the discharge is shown. The results are monotonic and positive. In figure 5 the evolution of the electric field on the axis of symmetry is also presented. The electric field magnitude at the head of the streamer has a value of 200kV/cm whereas in the ionized channel behind the propagating front towards the cathode the value is approximately 25kV/cm. It is observed that the shielding in the quasi-neutral channel behind the streamer head is of the order of 50% compared to the geometric field. Similar values for the electric field in the plateau of the discharge have been also been observed in [START_REF] Marode | Physics and applications of atmospheric non-thermal air plasma with reference to environment[END_REF]. In Figure 6 the electron density when the streamer is in the middle of the gap is presented. The discharge channel that is formed has a radius of approximately 50 µm, a value that is closely dependant on the choice of photoionization radiation model parameters [START_REF] Marode | Physics and applications of atmospheric non-thermal air plasma with reference to environment[END_REF], and on the width ( 0 ) of the initial Gaussian plasma spot. For 0 = 25 µm a much narrower streamer was observed in accordance with the results presented in [START_REF] Papageorgiou | and Three dimensional (3D) modeling of a streamer in air at atmosperic pressure incorporating photoionization phenomena[END_REF].

The streamer velocity and diameter increase as the streamer approaches the cathode. This behavior agrees with already published 2D results [START_REF] Kulikovsky | Two-dimensional simulation of the positive streamer in N2 between parallel-plate electrodes[END_REF] but deviates from experimental observations that show streamer propagation with a constant channel radius in plane to plane gaps [START_REF] Veldhuizen | Pulsed positive corona streamer propagation and branching[END_REF]. This difference is attributed to the longer gap employed for the experimental observations allowing the streamer radius to reach a steady state value.

Comparison of 2D axi-symmetric and 3D models

The comparison between 2D axi-symmetric and 3D results is very instructive as it will first of all provide confidence in the validity of the newly developed 3D model. In Figure 7 the initial photoionisation radiation from both models is shown. The solution of the 2D axi-symmetric and the 3D model gives almost identical results. It is worth noting that the three-exponential Helmholtz photoionisation model [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF] was also implemented giving a very similar discharge evolution. In our work the two-exponential model was finally adopted in order to reduce the computational cost. In Figures 8 and9, the electron density profiles on the axis of symmetry of the discharge during the avalanche and streamer phases for the two different models (2D and 3D) are presented. In Figure 10 the current density calculated by Sato's formula [START_REF] Sato | Discharge current induced by the motion of charged particles[END_REF] for the 2D and 3D model is also provided yielding very good agreement. The comparison is performed for exactly the same external conditions and mesh resolution. During the avalanche phase the resolution is set to 15 µm whereas during the streamer phase a moving window consisting of elements of size 5 µm follows the streamer head. At the avalanche and streamer phases the comparison between the 2D axi-symmetric and fully 3D results yield very good agreement (see Figures 8 and9). This is in disagreement with the results in [START_REF] Pancheshnyi | Numerical simulation of filamentary discharges with parallel adaptive mesh refinement[END_REF] where it is claimed that the 3D model gives more diffusive results compared to the 2D ones. In Figure 8 there are no differences in the density profiles whereas in Figure 9, during the streamer phase, the density profile of the 3D model is slightly ahead of the corresponding 2D profile. This difference is due to interpolation errors that are introduced in the 3D solution every time the moving window is updated to follow the streamer head. The discharge current (Figure . 10) for both cases is almost identical. The hump observed at about 3 ns is due to the absorption of the initial electron distribution at the anode (see Morrow and Lowke Figure 1 [START_REF] Kang | Numerical study on influences of barrier arrangements on dielectric barrier discharge characteristics[END_REF]).

Cathode Anode

Figure. [START_REF] Naidis | On streamer interaction in a pulsed positive corona discharge[END_REF] Comparison of the electron density profiles on the axis of symmetry between the 2D (red line) and 3D (blue line) model during the avalanche phase at 3.0 ns.

Figure. 9

Comparison of the electron density profiles on the axis of symmetry between the 2D (red line) and 3D (blue line) model during the streamer phase at 7.6 ns.

Cathode Anode Cathode Anode

Another crucial parameter in gas discharge numerical modelling and especially in 3D problems is the computational cost (see Table 3). The simulation time for the discharge evolution in 2D was approximately 20 minutes, whereas in 3D the calculation takes one week. In 2D the degrees of freedom (DoFs) are of the order of 36,000 while in 3D, when the streamer head is at the middle of the gap the resolution has to be very high, of the order of 4,000,000. The need for RAM in 3D is also very high and approximately 15 GB in this case. This shows clearly that 2D models have to be adopted in cases where the cylindrical symmetry is preserved and especially in large gaps; the computational cost of the 3D models is very high. The above simulations were performed on a system with four Intel Xeon 5500 series cores.

The effect of plasma spot and pressure inhomogeneities on streamer propagation

Following the validation of the 3D model, we will now use it to investigate the intrinsic parameters which are responsible for the formation of discharge branching and off-axis propagation. According to the literature at least four causes are involved in streamer branching: Laplacian instabilities [START_REF] Ebert | The multiscale nature of streamers[END_REF], pressure instabilities [START_REF] Babaeva | Effect of inhomogeneities on streamer propagation: I. Intersection with isolated bubbles and particles[END_REF], streamer interactions [START_REF] Briels | Positive streamers in air and nitrogen of varying density: experiments on similarity laws[END_REF] and photoionisation fluctuations [START_REF] Nasser | Impulse streamer branching from Lichtenberg figure studies[END_REF]. Here we consider the effect of plasma spots and low pressure instabilities (bubbles in [START_REF] Babaeva | Effect of inhomogeneities on streamer propagation: I. Intersection with isolated bubbles and particles[END_REF]) on streamer propagation. Plasma spots usually have a spherical form with a radius of several micrometers and a charge density within them that resembles a Gaussian distribution. They represent laser pulses or UV radiation that could be applied in the discharge volume close to the propagating front of the streamer. Additionally, a neutral plasma spot could also describe, in a simplified manner, the interaction of the streamer head with an already formed discharge channel [START_REF] Briels | Positive streamers in air and nitrogen of varying density: experiments on similarity laws[END_REF]. Plasma spots are introduced in the model as an additional source term in the continuity equations for electrons and ions. Their Gaussian form has a maximum value in the range of 10 20 -10 22 cm -3 s -1 and the spots are applied for a short subnanosecond time interval. In this work the effect of (i) the distance between the plasma spot and the streamer, (ii) the size of the plasma spot and (iii) the charge density in the plasma spot on streamer propagation is considered. Low pressure instabilities, on the other hand, can appear in the discharge volume in a chaotic way [START_REF] Babaeva | Effect of inhomogeneities on streamer propagation: I. Intersection with isolated bubbles and particles[END_REF]. The magnitude of these fluctuations is dependent on the external conditions of the discharge and is expected to be more significant in the atmosphere (e.g. lightning). In the proposed model they have a spherical spatial form (bubbles) and they are introduced off-axis close to the propagating front. In their volume the ratio E/N is elevated compared to that of ambient air, influencing in this way the transport parameters of the charged species. In this study the relative position of the bubbles and the pressure difference have been studied.

Plasma spots

In figure 11 the formation and propagation of a streamer in the presence of a plasma spot is presented. The plasma spot has been placed 100 µm from the axis, 110 µm away from the streamer head and the maximum charge density production rate is 10 22 cm -3 s -1 . This instability is applied for 0.1 ns and corresponds to the application of a very fast laser pulse in the discharge volume. Electrons under the influence of the electric field around the streamer head are extracted from the spot, accelerated towards the head and multiplied via ionization collisions. When their multiplication becomes important, Figures 11(b) and 11(c), the propagating front is deformed. A small diagonal subchannel (symbol A in Figure 11(c)) is formed and bridges the gap between the head and the spot. After the formation of the sub-channel, the streamer continues the propagation but off axis, parallel to the initial axis of the discharge (figure 11(d)). The main reason for this pattern of propagation is the strong external geometric field that forces any discharge to eventually move in a direction parallel to the field lines. In order to have diagonal propagation it is necessary to produce seed electrons in the diagonal direction. In an overvolted plane to plane gap the field lines feed seed electrons in directions perpendicular to the electrodes making the off-axis diagonal propagation almost impossible.

It is anticipated that most of the instabilities which arise in such a discharge gap will not be able to produce branching but small changes in the initial direction of the discharge channel. This becomes clear if we compare the evolution of the discharge channel without the addition of the plasma spot and the off axis propagation after introducing the plasma spot, figures 12(a) and 12(b).

Even if the Laplacian instability [START_REF] Ebert | The multiscale nature of streamers[END_REF] case is considered, it is questionable that diagonal branching could arise. What could be expected in that case is an initial splitting of the main streamer head, the formation for example of two additional streamer channels, a diagonal propagation of these channels for a while due to the repulsion of the heads and after, as the repulsion reduces with distance, again propagation parallel to the electric field lines. In our opinion the most probable scenario for branching in overvolted parallel gaps is the streamer interactions: reconnection [START_REF] Briels | Positive streamers in air and nitrogen of varying density: experiments on similarity laws[END_REF] or attraction due to non local photoionisation [START_REF] Luque | Interaction of streamer discharges in air and other oxygen-nitrogen mixtures[END_REF]. Reconnection is caused by electrostatic attraction of a subsequent streamer with a conducting channel left by an early streamer that has already reached the charged electrode with opposite charge and has changed polarity [START_REF] Briels | Positive streamers in air and nitrogen of varying density: experiments on similarity laws[END_REF].

As Luque shows [START_REF] Luque | Interaction of streamer discharges in air and other oxygen-nitrogen mixtures[END_REF] two streamers of the same polarity do not always repel. At atmospheric pressure and at certain distance away the streamers attract each other: a cloud of electrons is created between the two streamers which makes them coalesce into a single wider one.

We have performed a parametric analysis on how the magnitude and the position of the several discharge instabilities influence the discharge channel. In Figure 13(a) and 13(b) the influence of a Gaussian plasma spot placed 150 µm and 100 µm from the axis of symmetry is depicted. In this case the maximum number density of the spot is chosen to be 10 22 cm -3 s -1 . The role of the distance is quite obvious: the interaction between the streamer head and plasma spot is mainly due to the electric field of the space charge. For distances larger than a critical value, the electron multiplication (in the vicinity of the spot) via ionization collisions is not enough to disturb the propagation path (Figure 13(a)). When the spot is closer, the streamer head deforms thus revealing a tendency for splitting (Figure 13(b)). Keeping in mind that in order for the streamer channel to be deformed, the distance between the streamer head and plasma spot has to be lower than a critical one, the effect of the magnitude of the plasma Gaussian disturbance placed at 100 µm off axis was examined. In the previous cases the maximum density was A = 10 22 cm -3 s -1 , and in figures 14(a) and 14(b) results for A = 10 20 cm -3 s -1 and A = 10 22 cm -3 s -1 are presented respectively. The deformation of the streamer head is again obvious. However, as the maximum charge density of the plasma spot is reduced so is the splitting of the discharge channel, leading to the absence of off -axis propagation when the value is reduced below a critical value.

Low pressure fluctuations

The last instability examined deals with low pressure fluctuations that may arise as streamers propagate. In this work a spherical bubble of radius 40 µm at a distance of 50 µm and 100 µm away from the axis of the main discharge was introduced. Initially the pressure in the bubble was 700 Torr (~8% decrease). The ratio E/N inside the bubble is hence larger than outside. As the streamer head approaches, there is a further increase in E/N due to the field enhancement ahead of the avalanche front. As photoionization supplies electrons in the bubble, the elevated E/N produces a more rapid rate of ionization than in the rest of the gap. This ionization produces a source of electrons that drift upstream into the gap between the bubble and the approaching avalanche front. The drift velocity of the electrons is also higher in the bubble than outside. As a result, the streamer accelerates into the bubble [START_REF] Babaeva | Effect of inhomogeneities on streamer propagation: I. Intersection with isolated bubbles and particles[END_REF] (see Figure 15). In Figures 15(a), 15(b) and 15(c) the effect on streamer propagation of the bubble pressure and the distance between the bubble and the streamer head is shown. No branching or splitting of the front can be observed in this case. Further calculations have shown that a decrease in the bubble pressure of ~30% (this is a very extreme case) could deform considerably the shape of the streamer head but cannot produce diagonal streamer branching as can be seen clearly in Figure 15(c).

No streamer branching could be reproduced under these conditions with the application of the 3D model. What were clear and reproducible were the streamer head splitting and the off-axis propagation under specific conditions. It is believed that the small gap length and the electrode configuration plays an important role in the development of streamer branching. Numerical investigations in point to plane geometries as well as longer gaps have to be performed in order to investigate further the development and evolution of branching and the role of Laplacian instabilities. Figure. [START_REF] Zheleznyak | Photoionisation of nitrogen and oxygen mixtures by radiation from gas discharges[END_REF] Influence of the low pressure bubble instability on streamer propagation.

CONCLUSIONS

Results have been presented from a full three dimensional gas discharge model incorporating photoionisation radiation via the solution of the Helmholtz equations. This has been compared with the two dimensional counterpart. The Helmholtz equations were solved by incorporating the twoexponential fit approach. The avalanche and streamer propagation phase were analysed and a comparison between the 2D axi-symmetric and the 3D model was performed. The two different approaches gave very similar results during the avalanche and streamer phases. The computational cost of the 3D models is very high showing clearly that 2D models have to be adopted in cases where the cylindrical symmetry is preserved and especially in large gaps. The role of two different discharge instabilities in streamer branching was examined. Plasma spots and low pressure fluctuations during c) Electron density (cm -3 ) at 7.6 ns. Low pressure bubble (500 torr) at 50 µm from the axis of symmetry the streamer propagation were simulated. Depending on the magnitude and position of the plasma spot, deformations and off-axis propagation of the main discharge channel was observed. Streamer splitting was less sensitive to low pressure discontinuities. Pressure instabilities start to be important when their relative difference to the ambient air is more than 30%. The small gap length and plane to plane electrode configuration are parameters that are believed to limit streamer branching.

  used for the 3D simulations b) Closer look of the mesh used for the 3D simulations

Figure 1 .

 1 Figure1. Three-dimensional mesh used for streamer propagation. The mesh is gradually becoming finer following an onion-like pattern as approaching the axis of the discharge.

Figure. 2 Figure. 3

 23 Figure. 2 Electron density profile evolution on the axis of symmetry during the avalanche phase Anode Cathode

Figure. 4 Figure. 5 Figure. 6

 456 Figure. 4 Evolution of the electron density profile on the axis of symmetry during the streamer phase

Figure. 7

 7 Figure. 7 Initial photoionisation radiation on the axis of symmetry for the 2D (red line) and 3D (blue line) model

Figure. 10

 10 Figure. 10 Comparison of current density profiles between the 2D (red line) and 3D (blue line) model during the discharge evolution.

  a) Electron density (cm -3 ) at 7.2 ns. b) Electron density (cm -3 ) at 7.4 ns. c) Electron density (cm -3 ) at 7.6 ns.d) Electron density (cm -3 ) at 7.9 ns.

Figure. 11 Figure. 12 Figure. 13

 111213 Figure. 11Electron density evolution during the streamer phase having introduced a plasma spot (red doted circle) offaxis close to the streamer head. The maximum charge density production rate of the plasma spot is 10 22 cm -3 s -1 and the spot is applied for 0.1 ns.

Figure. 14

 14 Figure. 14 Influence of the magnitude of the plasma spot instability on streamer propagation.

a)

  Electron density (cm -3 ) at 7.6 ns. Low pressure bubble (700 torr) at 50 µm from the axis of symmetry b) Electron density (cm -3 ) at 7.6 ns. Low pressure bubble (700 torr) at 100 µm from the axis of symmetry.

  )

	e i n ∂ + ∇ ⋅ N t N t N t ∂ ∂ + ∇⋅ ∂ ∂ + ∇⋅ ∂	e e i i n n N W N W S i N W S -∇⋅ ∇ e n D N r r = = r	e	=	S	e
	e	e	e	e	e	ep e p	ph
	i	e	e	ep e p	pn p n	ph
	n	e	e	pn p n

Table 2 :

 2 Boundary conditions used

	Electrons	Ions	Voltage
	Anode

Table 3 :

 3 Computational effort for 2D and 3D simulations

	Model	DoF	RAM	Time
	2D	36,000	300 Mb	20 min
	3D	4,000,000 15 Gb	1 week

ACKNOWLEDGEMENTS

The authors would like to thank Drs Richard and Vivienne Morrow for some very stimulating technical discussions and editing of the paper.


M Kotsonis
email: m.kotsonis@tudelft.nl

S Ghaemi

L Veldhuis

F Scarano

Measurement of the Body Force Field of Plasma Actuators

Keywords: numbers: 52.75.-d, 52.70.-m, 47.85.L-, 47.80.Jk

Measurement of the body force field of plasma actuators

Introduction

Plasma actuators have been studied extensively in recent years, as flow control devises. Their low power consumption, lack of moving parts and robustness render them ideal for flow manipulation. Plasma actuators have been used in several studies aiming at separation control [START_REF] Post | Separation control on high angle of attack airfoil using plasma actuators[END_REF], turbulent drag reduction [START_REF] Jukes | Turbulent drag reduction by surface plasma through spanwise flow oscillation[END_REF], boundary layer control [START_REF] Jacob | Boundary layer flow control using ac discharge plasma actuators[END_REF][START_REF] Huang | Plasma actuators for separation control of low-pressure turbine blades[END_REF][START_REF] Seraudie | Effect of plasma actuators on a flat plate laminar boundary layer in subsonic conditions[END_REF], noise reduction [START_REF] Li | The use of plasma actuators for bluff body broadband noise control[END_REF] and laminar-turbulent transition delay [START_REF] Grundmann | Active cancellation of artificially introduced Tollmien-Schlichting waves using plasma actuators[END_REF]. Comprehensive reviews on plasma actuators for aerodynamic flow control have been published recently [START_REF] Moreau | Airflow control by non-thermal plasma actuators[END_REF][START_REF] Corke | Dielectric barrier discharge plasma actuators for flow control[END_REF].

There are several implementations of the actuators the most popular being the Dielectric Barrier Discharge (DBD) (figure 1). These are based on the ionization of air using AC High Voltage (HV). In the most conventional configuration, two electrodes are employed. The electrodes are separated by a dielectric layer which prohibits arc forming and allows accumulation of the ionized gas in the vicinity of the exposed electrode. The exact mechanisms of interaction between the weakly ionized gas and neutral air are still debated. It is, nevertheless, commonly accepted that collisional processes between the heavy plasma species (mostly ions) and neutral air is responsible for the momentum transfer [START_REF] Likhanskii | Numerical modeling of dbd plasma actuators and the induced air flow[END_REF]. At a macroscopic scale the model of an exerted body force on the fluid seems to describe the effect with reasonable accuracy.

A large amount of simulation studies has been conducted in order to simulate the underlying physics of the ionization process [START_REF] Jayaraman | Modeling of dielectric barrier discharge plasma actuator[END_REF][START_REF] Lagmich | Model description of surface dielectric barrier discharges for flow control[END_REF][START_REF] Boeuf | Electrohydrodynamic force and aerodynamic flow acceleration in surface dielectric barrier discharge[END_REF]. These vary in model complexity, from simple phenomenological models to first principles fluid models [START_REF] Jayaraman | Modeling of dielectric barrier discharge-induced fluid dynamics and heat transfer[END_REF]. Extended simulations for multi-species fluids have also been investigated [START_REF] Singh | Modeling of dielectric barrier discharge plasma actuator with atmospheric air chemistry[END_REF]. In the majority of these modeling approaches the final goal is to determine the exerted body force on the fluid. A number of studies on flow solvers implementing plasma actuators have successfully coupled the effect of the actuator with the flow dynamics [START_REF] Orlov | Modeling and experiment of leading edge separation control using sdbd plasma actuators[END_REF][START_REF] He | Numerical and experimental analysis of plasma flow control over a hump model[END_REF]. The capabilities of the actuators suggest the feasibility of their implementation in several flow control scenarios. It is, therefore, desirable to have an efficient and accurate model of the actuator effect.

The experimental investigation on the actuator's operation has also received substantial attention. Several studies focus on the investigation of the induced velocity field [START_REF] Enloe | Time-correlated force production measurements of the dielectric barrier discharge plasma aerodynamic actuator[END_REF][START_REF] Forte | Optimization of a dielectric barrier discharge actuator by stationary and non-stationary measurements of the induced flow velocity: Application to airflow control[END_REF][START_REF] Opaits | Surface charge in dielectric barrier discharge plasma actuators[END_REF] and the produced thrust from the actuators [START_REF] Abe | Experimental study for momentum transfer in a dielectric barrier discharge plasma actuator[END_REF][START_REF] Dong | Experimental study of a DBD surface discharge for the active control of subsonic airflow[END_REF][START_REF] Hoskinson | Force measurements of single and double barrier dbd plasma actuators in quiescent air[END_REF]. These methods are robust and fast but they present some drawbacks which render the results difficult to be used as an input in flow solvers, an area that so far has been served only by numerical plasma modeling studies. An inherent disadvantage is that a part of the measured thrust does not correspond to the actuator effect but to the shear forces developed between the flow and the surface. Another shortcoming is also their inability to provide any information on the spatial distribution of the force field, information that is vital in cases of discrete, unsteady and localized actuation essential to the majority of active closed-loop flow control concepts. To the authors knowledge no experimental investigation on the determination of both the magnitude and the spatial distribution of the body force vector field has been published to this time.

This paper aims at the development of a novel experimental approach for determining the amplitude and spatial distribution of the body force exerted on the flow by a plasma actuator. The proposed method will be evaluated against well established thrust measurement techniques used as a benchmark. A description of the methodology for the new body force method as well as the benchmark thrust measurement techniques will be given in the following section. A description of the experimental setup and a discussion on the results is finally provided.

Methodology

The determination of the induced body force can be done experimentally in several ways. Due to the nature of the momentum transfer mechanisms direct measurements of the ionized species density in the plasma region can only be obtained using molecular level diagnostics. The species density can be used to derive the Coulomb forces exerted by the electric field. Taking into account the complexity of such approach, only indirect methods based on the induced flow field or reaction forces can be feasibly employed in the time span such an experiment allows.

In this study two techniques are applied to determine thrust produced by the operation of the actuator. These will be used as benchmark cases for validating and verifying the developed force estimation technique. More specifically the techniques employed are: direct thrust measurements using a highly sensitive load cell and thrust calculation via momentum balance of the induced flow field.

Body Force estimation

The proposed body force measurement technique determines the spatial distribution of the plasma induced body force from time-resolved velocity field. The required information is obtained using a high-speed PIV system as it is shown in a subsequent section of this paper. An analysis of the velocity field using the full Navier-Stokes (NS) equations can provide the body force. The 2D incompressible NS equations in the presence of body forces read as:

∂U ∂t + U • ∇U -ν ∇  U = -∇p + F ρ ( 1
)
where U is the 2D velocity field, p is the static pressure, ν is the kinematic viscosity of the fluid and ρ is the density. A first approach in deriving the body force (F ) is to apply (1) only for the first moments after the actuation. During this stage the flow is expected to accelerate only around the region of the exerted body force where the rest of the field is still relatively quiescent. For initial conditions of zero velocity the convective, viscous and pressure gradient terms are assumed relatively small and can be neglected. This reduces (1) to:

∂U ∂t = F ρ (2)
which is basically Newton's second law. This method will be referenced as the 'reduced method ' throughout this paper.

In order to use the full NS equations to derive the force it is necessary to calculate all the terms involved. The acceleration, convective and viscous terms are obtained from the available spatio-temporal data of the velocity field. In the absence of body force terms the only unknown is the pressure gradient which can be obtained using (1) from the measured velocity field. In the case of the plasma actuator, the body force term appears as one extra unknown. Physically, the pressure gradient term and the body force term are equivalent and one can be used to explain the phenomena related to the other. Nevertheless, the expression of the plasma effect through a volume distributed body force is a convenient way of imposing the effect in a flow solver. To be able to bypass the problem of one extra unknown, two major assumptions must be made: 1) the body force remains quasisteady over a large number of HV cycles and 2) the pressure gradient prior to the actuation is zero. The assumption for a quasisteady force is justified by taking into consideration the difference in time scales involved between one HV cycle and the acceleration of the flow. As mentioned by Enloe et al [START_REF] Enloe | Time-correlated force production measurements of the dielectric barrier discharge plasma aerodynamic actuator[END_REF], the body force production is not constant within one HV cycle. Nevertheless, the high frequency of operation (order of kHz) and the hydrodynamic 'inertia' of the flow reduce the effect to quasisteady in time. With these assumptions, (1) is differentiated in time:

∂ 2 U ∂t 2 + ∂(U • ∇U ) ∂t -ν ∂(∇ 2 U ) ∂t = - ∂(∇p) ∂t (3)
The first assumption states that the body force is steady in time and its time derivative is reduced to zero. Equation (3) now involves only one unknown which is the time derivative of the pressure gradient. Integrating (3) back in time the pressure gradient can be calculated from:

t 0 ( ∂ 2 U ∂t 2 + ∂(U • ∇U ) ∂t -ν ∂(∇ 2 U ) ∂t )dt = -(∇p) + A (4)
In the process of integration a constant A appears which has to be defined. At this point the separation between the pressure gradient due to the flow and the plasma body force is established. Based on the second assumption ((∇p)| t=0 = 0) the integration constant is set to zero. Equation ( 4) is plugged back into (1) and leaves the body force term as the only unknown. Since the time gradient of the NS is used, this method will be referenced as the 'gradient method ' throughout this paper.

In contrast to the two benchmark techniques the reduced and gradient methods have the advantage of providing spatial distribution of body force vectors instead of an integrated thrust value. The results from such analysis can be valuable for validation of numerical models of plasma actuators as well as an input into flow solvers investigating flow control concepts.

Direct thrust measurement

Direct thrust measurement is a straightforward technique which involves attaching the actuator to a load cell and measuring the force exerted by the induced flow on the actuator. By means of Newton's third law this force is equal and opposite to the total force the actuator exerts on the flow. The value of thrust can be valuable for parametric and optimization studies aiming at improving the actuator. Nevertheless, further clarifications must be made on what exactly is being measured. It has been reported repeatedly and also verified in this work that the induced flow field from the continuous operation of the actuator resembles a typical wall jet. As such, shear forces developing between the wall surface of the actuator and the accelerating fluid are also included in the load cell measurement. Another issue with this technique and in general with integral techniques is the lack of spatial information.

Thrust calculation from velocity measurements

The second benchmark technique is based on measurements of the velocity field in the vicinity of the plasma actuator, obtained by PIV. It involves the definition of a control volume in the flow where the momentum flux is calculated at the boundaries of the domain and equilibrium with the internal forces is assumed (figure 2) ‡. This approach is well established in experimental and numerical investigations for the calculation of aerodynamic drag of airfoils and wings [START_REF] Anderson | Fundamentals of Aerodynamics[END_REF]. In a recent experimental study [START_REF] Hoskinson | Force measurements of single and double barrier dbd plasma actuators in quiescent air[END_REF], it was used for estimating the thrust of plasma actuators. The calculated thrust includes both the body force from the actuator and also the friction force at the wall. This technique also suffers from lack of any information on the spatial distribution of the body force. More specifically applying the momentum balance equation on the control volume gives:

T = abcd (ρU • dS)U + abcd pdS ( 5
)
where T is the integrated force exerted on the fluid. Equation ( 5) is developed for the boundaries taking account the no-slip condition at the wall. This implies that the ‡ For all subsequent figures and references to spatial dimensions and distances, the displayed coordinate system is used. The origin lies at the downstream end of the exposed electrode and on the flat surface momentum flux across the wall is zero. Note that the boundaries where the unit vector perpendicular to the surface (dS) is negative take a minus sign:

• x -direction T x = ρ da u x 2 dy + ρ cd u x u y dx -ρ bc u x 2 dy + da pdy - bc pdy (6)
• y -direction

T y = ρ da u x u y dy + ρ cd u y 2 dx -ρ bc u x u y dy + cd pdx - ba pdx (7)
Due to the fact that pressure is not readily available from velocity data, a further assumption needs to be made. This states that if the control volume boundaries are far enough from the bulk of the plasma body force (near the inner electrode edge) pressure can be considered uniform and equal. Similar assumption is typically made in the case of airfoil drag estimation. It should also be noted here that for the x-direction the calculated force (T x ) contains not only the contribution of the plasma actuator (T p x ) but also the shear force between the flow and the surface (T s x ). These assumptions reduce ( 6) and ( 7) to:

• x -direction T p x + T s x = ρ da u x 2 dy +ρ cd u x u y dx -ρ bc u x 2 dy (8)
• y -direction

T y = ρ da u x u y dy +ρ cd u y 2 dx -ρ bc u x u y dy (9)

Experimental Setup

The actuator

In this investigation DBD actuators are fabricated using thin rectangular copper electrodes made out of self-adhesive copper tape (figure 1). The electrodes are 10 mm in length (x direction), with zero horizontal gap and 60 µm of thickness. Their effective spanwise length (z direction) is 300 mm for the load cell measurements and 200 mm for the PIV measurements. The electrodes are separated by two dielectric layers of polyimide Kapton tape. The thickness of each layer is 50.8 µm. The total thickness of the dielectric, including adhesive, is approx. 110 microns. The upper electrode is connected to the HV output cable of a TREK 20/20C HV amplifier (20 kV pp , 20 mA) while the lower electrode is grounded. The actuator operation is controlled remotely via a computer workstation were the driving signal is created in LabView software and is sent to the amplifier via a Digital/Analog (D/A) converter. The amplifier provides direct readings of the output voltage and current through internal measurement probes. While the internal voltage probe gives sufficiently accurate readings the internal current probe has been found to be too slow to resolve high frequency current fluctuations. This is typically the case with discharge currents occurring during plasma actuators operation. To resolve this, a resistance is placed between the lower electrode and the grounding cable and voltage is measured across it indicating the discharge current.

To reduce the effect of any external disturbances on the measurements, the tests are conducted in a closed box made from Plexiglass (PMMA) to provide optical access. This box is used in both the load cell and PIV experiments. Applied voltages of 8 to 16 kV pp and carrier frequencies of 1 to 4 kHz are tested. The waveform of the applied voltage is a sine. The thrust measurements are a part of an extensive work [START_REF] Kotsonis | Experimental study on the body force field of dielectric barrier discharge actuators[END_REF] which involves additional geometric variations of the actuator. An overview of the actuator and test parameters for all measurements is shown in table 1.

Load cell

For the direct thrust measurements an electronic ME-messsysteme KD-40S strain-gage load cell is used. The nominal force range is ± 2 N with an accuracy of 0.1% (2 mN ). A light PMMA plate of 5 mm thickness, which carries the electrodes and dielectric, is directly mounted on the load cell. The connection of the electrodes to the HV and ground connectors is via thin copper wires to minimize signal contamination from external vibrations. Tests are repeated three times to minimize any error due to load cell relaxation. The measurement is continued for ten seconds before actuation and thirty seconds after in order to ensure that the steady flow regime is achieved. The thrust is then calculated by comparing the time averaged readings from the two time periods.

PIV setup

Time-resolved PIV has been applied to characterize the flow field in the vicinity of the plasma actuator. This technique provides the required spatial and temporal resolution in order to characterize the transient behavior of the thin wall jet induced by the actuator.

A two component PIV configuration has been chosen since the large span of the actuator ensures minimal 3D effects. The two-dimensionality of the flow is also verified by hotwire measurements. The plasma actuator is placed flushed on the bottom of the Plexiglas box. A Photron Fastcam SA1 high speed CCD camera of 1024×1024 pixels (full sensor size) is used to image the field-of-view (FOV). Image acquisition has been conducted at 10 KHz rate in single-frame mode. This implies a time separation of 100 µs between successive images. A Nikkor 105 mm objective is set at f-stop 4 and is used along with extension tubes in order to achieve 0.8 magnification and a FOV of 15×6 mm. To achieve the high sample rate needed for these measurements, the sensor is cropped to a size of 1024×512 pixels. The air in the Plexiglas box is seeded with olive oil droplets of approximately 1 µm diameter generated by a TSI atomizer. The particles at the mid span of the actuator are illuminated by a light sheet of 2mm thickness generated by a Quantronix Darwin-Duo laser system with an average output of 80 W at 3 kHz. The images are analyzed using Davis 7.4 (Lavision GmbH) by cross-correlating successive images. Final interrogation window size of 12×12 pixels and overlap factor of 75% are used. The interrogation windows are elongated in the wall normal direction using a 4:1 aspect ratio in order to obtain higher spatial resolution. The velocity vectors are returned on a grid of 13 vectors per mm.

Results

Load cell thrust measurements

The load cell measurement offers a robust approach and has been applied to a large number of test cases. Detailed results can be found in [START_REF] Kotsonis | Experimental study on the body force field of dielectric barrier discharge actuators[END_REF]. In this paper the load cell results are only presented for the voltage (V ) and frequency (f ) values identical to the one used for the PIV measurement.

The test cases studied in this paper are shown in table 1. For voltage variation three cases with base frequency of 1, 1.5, and 2 kHz are done, while for frequency variation four cases with base voltage of 8, 9, 10 and 11 kV pp are performed. For the frequency variation the maximum tested frequency is limited by the slew rate of the HV amplifier (350V /µs) and is reduced with increase of voltage amplitude. For the entire study, a sinusoidal waveform has been used for the high voltage signal.

Thrust results for voltage and frequency variation are presented in figure 3. Both voltage and frequency seem to influence the magnitude of thrust albeit in different manners. Thrust is related to voltage with a power law. Similar behavior has been previously observed by Abe et al [START_REF] Abe | Experimental study for momentum transfer in a dielectric barrier discharge plasma actuator[END_REF]. In contrast, the relation with frequency is linear.

PIV measurements: estimation of thrust

The dataset from the high speed PIV measurement is used for the thrust estimation using the momentum balance technique. After the initial acceleration stage, the induced flow field develops into a wall jet. In the steady state regime the momentum balance equations do not require evaluation of the acceleration terms. Steady state is reached in approximately 30 ms after actuation for the case of 8 kV pp . In this respect it is safe to average for approximately 400 frames spanning from 43 ms after the actuation to 83 ms. Apart from the processing of the raw PIV data and time averaging no other data processing is required. The time averaged flow field with the selected control volume is shown in figure 4 for the case of 12 kV pp . The momentum increase mainly occurs in the x direction while a weak suction effect is observed upstream of the inner edge of the electrodes. The suction effect is an indication of a pressure gradient which is potentially strong in the vicinity of the actuator. As mentioned, any existing pressure gradient cannot be readily resolved using the time averaged velocity data. It is therefore necessary to choose the control volume such that the pressure effects are minimized. The chosen control volume is a rectangular area spanning form x = -3 to x = 5 mm and y = 0 to y = 2.5 mm.

The momentum balance equation in x and y directions calculated using the terms in ( 8) and ( 9) reveals the significance of each boundary of the control volume. As it is shown in table 2 and also apparent from the velocity field, the maximum contribution to thrust comes from the outflowing horizontal velocity from section bc. This implies that the very simple approach of calculating thrust from only one cross section of the jet would still deliver satisfying results. Nevertheless, since field data are available the method used herein takes into account all sections of the control volume.

The thrust as a function of both applied voltage and frequency are shown in figure 5. The same power law behavior as observed in the direct thrust measurements is observed here for thrust in x direction with different voltages. A slight increase of thrust magnitude in y direction is observed with increase of the applied voltage. For different frequencies, the linear dependence of thrust appears to hold. A more concise comparison of the results from different techniques is performed at a subsequent section.

PIV measurements: estimation of body force

The two described approaches, (reduced method ) and (gradient method ), are applied to the velocity field measured using PIV. Instantaneous snapshots of the evolving flowfield are presented in figure 6. The time values indicate elapsed time from the start of the actuation. The flow reacts immediately to the actuation with initial velocity components already starting to appear approximately 0.1 ms after the start of the actuation. The flow accelerates in a region starting from the end of the upper electrode and extending over half the covered electrode (x = 0 to 5 mm). The impulsive forcing on the flow creates a starting vortex which develops and moves out of the FOV in less than 40 ms after the start of the actuation. It is apparent that the general forcing of the actuator produces a wall jet which in turn creates a secondary weak suction effect upstream. At this point, some specific features of the induced flow field should be discussed due to their effect on the results of the proposed technique. In a previous study [START_REF] Benard | Capabilities of the dielectric barrier discharge plasma actuator for multi-frequency excitations[END_REF], high frequency velocity fluctuations are registered in the vicinity of the actuator. These fluctuations are a product of the carrier frequency at which the actuator is operating and reflect the underlying ion-neutral collision processes in the plasma region as reported in previous numerical studies [START_REF] Likhanskii | Modeling of dielectric barrier discharge plasma actuator in air[END_REF]. These fluctuations appear in the PIV measurements of this work as well, corresponding to the carrier frequency for each case. More in depth analysis of this phenomenon does not fall in the scope of this study. The proposed methods rely on the resolution of the accelerating flow and interpretation of that to a body force. In this process the first and second time derivatives of velocity are required. Due to the high frequency fluctuations, these terms cannot be explicitly resolved without large artifacts. Based on the large difference between the time scale governing the acceleration of the flow and the time scale of these fluctuations one further assumption must be made. This states that although the flow, and in extent the body force, fluctuate in time the latter can be considered as quasi steady. This assumption can be justified by the fact that the rise time of the flow (accelerating part) is considerably larger than the fluctuation period. This is also a major assumption made in a large majority of first-principles modeling studies on the plasma operation. Based on this, the PIV time-series of the flow field snapshots is filtered in time with a 6th order lowpass Chebychev filter as shown in figure 7. The cutoff frequency is chosen based on the carrier frequency of the HV signal corresponding to the test case under investigation. For the majority of cases the cutoff frequency is half the carrier frequency. By filtering the high frequency fluctuations, the development of the velocity field occurs smoothly and the time derivatives can be calculated with reasonable accuracy.

The reduced method involves the utilization of the very first moments (first 0.5 ms) after actuation when all other terms apart from acceleration are negligible. The acceleration, convective, viscous and pressure gradient components of the NS are shown in figure 8 as a function of time after actuation (t = 0). The terms are calculated in x direction from the time-filtered PIV data for a given point located well into the expected force field. As anticipated, the acceleration term dominates the very first moments of actuation with a high and almost constant value. Convection and pressure gradient are initially negligible and eventually dominate the event while viscous effects have a

∂ U ∂ t + U • ∇ U -ν ∇ 2 U = -∇p + F ρ reduced method
x=3 mm, y=0.5 mm x small influence during the entire actuation period. It should be mentioned here that the pressure gradient term in figure 8 is calculated using the gradient method as discussed in section 2.1. The initial acceleration-dominant time span is visually confirmed to be 0.5 ms and provides reasonable results. The spatial distribution of the body force, as calculated by the reduced method, is shown in figure 9 for the case of 12kV pp applied voltage and 2kHz carrier frequency. The PIV snapshots used for this case cover the first 0.5 ms after actuation and the force field is simply averaged over this span. The bulk of the force appears to be just downstream the inner electrode edge while the horizontal extent spans almost a quarter of the covered electrode. In the wall normal direction the force field extends no more than 1 mm which confirms the very thin wall jet reported by this and numerous other investigations. The vertical force seems to be very weak compared to the horizontal component. Its major direction is towards the wall especially near the edge of the exposed electrode. For this case the cutoff frequency of the Chebychev filter is 1000 Hz.

For the same case (12 kV pp , 2 kHz), the spatial distribution of the body force using the gradient method is shown in figure 10. The first 20 ms after actuation are used and the force field is averaged over this time span. The general shape of the field is identical to the results from the reduced method. Nevertheless, the magnitude of the force attained with the gradient method is slightly larger than the reduced method's result. This small discrepancy can be attributed to modeling errors inherent in the eliminated terms of the reduced method as well as the ambiguity in selecting the initial acceleration time span. These facts support the higher accuracy of the gradient method.

One additional advantage of the gradient method over the reduced method is its relative independence on the time span used for the force calculation. In contrast to the reduced method there is no term elimination involved and no consideration should be made in the dominance of one term or another. Applying this method for the first 20 ms after actuation gives an almost steady in time body force. Several snapshots of the instantaneous total force field are shown in figure 11. The temporal coherence of the force distribution is an extra indication of the validity of the method as this was an important assumption made during the formulation of the technique.

The force distributions for different voltage and frequencies are presented in figure 12 and figure 13, respectively. For voltage variation the trend of increasing force with increased voltage holds. Of interest is the spatial extent of the force field along the dielectric. As the voltage increases, the active phase of the negative half cycle, where the electron avalanche occurs, increases in duration as reported in the excellent paper by Likhanskii [START_REF] Likhanskii | Modeling of dielectric barrier discharge plasma actuator in air[END_REF]. This allows the plasma to propagate further downstream. Respectively, in the positive half cycle the corona regime is sustained for a longer period enhancing the movement of positive ions downstream. This is also verified by the force distributions for different frequencies. In this case, the spatial extent of the plasma cloud remains almost the same as the applied voltage is fixed. The magnitude of the force, nevertheless, is increased since the momentum exchange between the ionized gas and neutral air is performed at a higher rate. Furthermore, a jump in maximum force also occurs from the case of 3 kHz to the case of 4 kHz as well as a secondary vertical force component towards the surface. This is an indication that 4-5 kHz presents a threshold between two different regimes of operation. A low frequency regime exists (below 3-4 kHz for these plasma scales), where the frequency is slow enough to allow full electron-ion recombination between the negative and positive half cycles. In contrast, for higher frequencies (above 4-5 kHz) the electron cloud remaining on the dielectric from the negative cycle is not completely destroyed before the positive cycle begins. The residual electrons are then pulled back to the exposed electrode causing further avalanche ionization as they move. As reported by Likhanskii, in this case the generated positive ions move towards and along the dielectric surface imparting momentum to the air. This might explain the secondary vertical force components located around x = 3mm in the 3 and 4 kHz cases. This is also mentioned by Soloviev and Krivtsov [START_REF] Soloviev | Surface barrier discharge modelling for aerodynamic applications[END_REF] where the discharge formation appears to be affected by the background particle concentration fields left from the previous HV cycles. Soloviev and Krivtsov suggest that the momentum input is dependent on the temporal overlap between two subsequent discharge cycles.

Comparison

An overall comparison between the thrust values attained by the proposed body force estimation methods and the two benchmark techniques is presented in figure 14 for different voltages and frequencies. It should be noted here that the thrust value from the proposed PIV techniques is calculated by integrating the body force value over the

spatial domain.

For different voltages, a good agreement is observed between the different techniques. In x direction (horizontal body force) the momentum balance is located at the lower error margin of the load cell measurements. The thrust estimated by the body force methods is systematically larger than the two other techniques although still in the error margin of the load cell for the majority of cases. This is expected as both, the load cell measurements and the momentum balance calculations contain the shear stress at the wall. Additionally, the previously mentioned difference between the reduced and gradient methods is observed. This is attributed to the term elimination inherent in the reduced method.

For different frequencies, the absolute error stays in the same level while the relative error between the different methods is larger. This is partly due to the low applied voltage of 10 kV pp that is used for this test series as the slew rate of the currently used HV amplifier is the limiting factor in increasing both voltage and frequency. The low voltage results in lower body force and the results are more susceptible to interference from external sources and other measurement errors. It is expected that the relative error would decrease with higher applied voltage.

Conclusions

A novel technique for the estimation of the body force field of plasma actuators is developed and tested. This involves the use of time-resolved PIV data on the evolving flowfield during actuator operation. The force amplitude and spatial distribution can be calculated through the estimation of the individual terms of the 2D incompressible NS equations and is conducted in two distinct ways. First, a dimensional analysis is performed on the temporal and spatial evolution of the velocity field. An initial period dominated by pure acceleration is identified and the latter is expanded explicitly in body force terms. Second, the full NS is estimated. The decoupling between the force and pressure gradient terms is performed based on the temporal stability of the body force and the initial zero pressure gradient.

The results from the tested methods are compared against well-established thrust measurement techniques and demonstrate sufficient agreement. Based on the comparison several comments can be made on the proposed techniques. The reduced and gradient techniques produced consistently larger values of horizontal force than direct thrust measurements or momentum balance on the velocity data. This is mainly attributed to the inherent inclusion of the wall shear force, opposite to the plasma body force, in the total thrust value obtained by the two benchmark techniques. The proposed techniques negate this effect since the calculation is done in the fluid volume not including the wall effects. Furthermore, a comparison can be made between the reduced and gradient methods. Although attractive because of its simplicity, the reduced method presents some drawbacks in respect to the gradient method. These are based on term elimination as well as on the semi arbitrary choice of applied time length.

The general shape of the force field seems to be dictated by voltage and frequency for the current geometric configuration of the actuator. The shape resembles an elongated ellipse with the maximum force arising near the inner electrode edge. Increasing voltage has two major effects on the field being an increase in amplitude and an extension in the downstream direction over the covered electrode. This is attributed to the increase of the electric field and consequently the extension of the active phase for each of the two half cycles of the actuation signal. Compared to the voltage, the relation with frequency is slightly different since there is no visible change in shape but only in amplitude. A possible jump in operation regimes is visible for higher frequencies. A residual electron avalanche in the positive half cycle seems to occur for the cases of 3 and 4 kHz although no definite conclusion can be made at this point.

Figure 1 .

 1 Figure 1. Geometrical configuration and operation of the plasma actuator.

Figure 2 .

 2 Figure 2. The control volume for the momentum balance equations.

Figure 3 .

 3 Figure 3. Thrust (T x ) per length of actuator for different applied voltages and frequencies, as measured by the load cell.

Figure 4 .

 4 Figure 4. Time averaged velocity field and control volume for momentum balance (12kV pp , 2 kHz).

Figure 5 .

 5 Figure 5. Thrust per length of actuator for different applied voltages and frequencies, as calculated by momentum balance on the PIV data.

Figure 6 .

 6 Figure 6. Temporal evolution of the total velocity field (U) for the actuation case of 12 kV pp , 2 kHz.

Figure 7 .

 7 Figure 7. Raw and filtered signal of velocity in the developing flowfield.

Figure 8 .

 8 Figure 8. Temporal evolution of the NS terms governing the development of the flow in x direction ( 12kV pp , 2kHz).

Figure 9 .

 9 Figure 9. Spatial distribution and direction of the body force using the reduced method (12kV pp , 2kHz).

Figure 10 .Figure 11 .

 1011 Figure 10. Spatial distribution of the body force using the gradient method (12kV pp , 2kHz).

FFigure 12 .

 12 Figure 12. Body force distributions for different voltages as calculated by the gradient method (f = 2kHz).

Figure 13 .

 13 Figure 13. Instantaneous total force field obtained using the gradient method (10kV pp ).

Figure 14 .

 14 Figure 14. Thrust values from the tested methods as a function of applied voltage (a) and frequency (b).

Table 1 .

 1 Test cases for thrust and body force investigation.

	Parameter	Value	Base value
	upper electrode length (l u )	1 0mm	-
	lower electrode length (l l )	1 0mm	-
	dielectric thickness (t d )	1 1 0 µm	-
	horizontal gap (g)	0 mm	-
	Load cell measurements
	applied voltage (V app )	8t o1 4kV pp (steps of 0.5 kV )	8, 9, 10, 11 kV pp
	frequency (f )	0 . 5t o4kHz (steps of 0.25 kHz)	1, 1.5, 2 kHz
	PIV measurements
	applied voltage (V app )	8t o1 6kV pp (steps of 2 kV )	1 0kV pp
	frequency (f )	1 t o 4 kHz (steps of 1 kHz)	2kHz

Table 2 .

 2 Section contribution in thrust (mN/m act ) for case of 12 kV pp , 2 kHz.

	section T x	T y
	da	-0.131 0.0015
	cd	0.0160 -0.1344
	bc	3.9828 -0.0235
	total	3.8673 -0.1564

Acknowledgments

This work is done in the framework of the CleanEra project at the TU Delft (www.cleanera.tudelft.nl). The authors would like to acknowledge Leo Molenwijk and Stefan Bernardy for their technical assistance.