A three-dimensional model of streamer discharges in unsteady airflow

A 3D fluid model has been developed to simulate streamer discharges in unsteady airflow. The model couples the drift–diffusion equation for charged particles, the Navier–Stokes equations for air, the Poisson’s equation for the electric field, and the Helmholtz equation for photoionization. It allows us to study electrical discharges at different timescales defined by light and heavy particles and to investigate the effects of unsteady airflow. The model treats the time integration in an implicit manner to allow longer time steps, which makes the simulation of long-duration discharges feasible. Moreover, the model uses an unstructured mesh allowing the calculation around solid bodies with complex geometries, and uses adaptive mesh refinement to lower the computation time. The validity and accuracy of the model has been verified by comparing its results with published results, which compares simulations in steady air from six different streamer codes. Our results are consistent and among the most accurate in terms of charge conservation. In order to investigate the influence of wind on streamer discharges, we present results from simulation of a long-duration discharge, in which two successive positive streamers are initiated from a positive polarity electrode in presence of a transverse airflow. This simulation shows that the impact of airflow on positive streamers is driven by the ions, and therefore the airflow effects are seen in ions timescale. Interestingly, we observe that the positive streamer channel, while tilting in the direction of the wind, remains attached to the surface of the electrode. The subsequent positive streamer emerges from the charges remaining from the initial streamer, which have been moved over the electrode surface toward the trailing edge. This mechanism shows and explains the clear tilting of the successive positive streamers in the direction of the wind.


Introduction
Streamers are thin ionized channels that initiate large-scale electric discharges such as sparks and lightning, and a deeper * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. understanding of their nature will shed light on the fascinating dynamics of long sparks. Streamers have been extensively studied both by numerical and experimental means [1][2][3][4][5]. Streamer dynamics are highly nonlinear and small changes in the discharge conditions can lead to large variations in streamer dynamics. The fast propagation and small dimensions of streamers do not allow experiments to study the microscopic characteristics of streamers with high resolution. Thus, numerical investigation of streamer discharges have been used to understand streamer dynamics, and numerous models of varying complexity have been developed. We can for instance mention particle, hybrid particle/fluid and fluid models [3,[6][7][8][9].
In this paper, we use a fluid description for cold atmospheric plasma. We are interested in unsteady fluctuations in air density and also the influence of the air velocity. A few streamer models have tried to study the streamer propagation in non-uniform air [10,11]. There are two main physical mechanisms from which unsteady airflow might affect the electric discharges. Firstly, the breakdown electric field is modified due to variations in air pressure, which consequently affects the streamer dynamics [10,12]. Secondly, the airflow affects the motion of charged particles by transferring momentum. The role of this mechanism in electric discharges can be significant, yet only a few research papers have investigated it [13,14]. Therefore, we have included Navier-Stokes and turbulence equations to model the airflow.
The model can study these effects in an accurate and realistic manner by coupling electric discharge set of equations with Navier Stokes and turbulence equations. It is fully three dimensional, allowing calculation of realistic flows around complex structures, and of streamers branching and interactions with each other [15]. The model also allows the incorporation of structures with complicated geometries and is built on an unstructured mesh to capture the curves and complexities of the structure. Unstructured meshes have previously been used in streamer simulations [16].
The aim of the model is to capture the effects of unsteady airflow on electric discharges. Part of these effects involve the motion of ions, and therefore occur on longer timescales than the streamer timescales. Our model integrates implicitly over time, removing the Courant-Friedrichs-Lewy (CFL) restriction on the time step size, which allows us to model longduration electric discharges of hundreds of microseconds.
The model is built on the commercial software Ansys Fluent. Ansys Fluent is based on the finite-volume method and allows the implementation and modification of equations and numerical schemes associated to electrical discharges. Furthermore, it supports adaptive refinement on unstructured mesh which is considered necessary for our model.
In the following sections, we describe the physics of the model, its implementation, and the numerical schemes used. Later, we verify the model accuracy by comparing to a test case from [17]. In the last section, we present a simulation of longduration discharge which shows the ability of the code to couple the airflow and electric discharge dynamics, and explains the influence of the wind on the discharge.

Physical model
Charged particles dynamics are modeled by drift-diffusion equations [18][19][20]. + n e ν i − S r e−p − S r n−p , (1.2) ∂n n ∂t = ∇ · (D n ∇n n ) + ∇ · (v air + μ n En n ) + n e ν a − S r n−p , (1.3) where n e , n p , and n n are the densities of electrons, positive ions and negative ions and E is the electric field vector. The loss terms S r e−p and S r n−p describe the charged particle losses due to electron-ion and ion-ion recombination processes [21]. The scalar parameters D and μ are the diffusion and mobility coefficients, and ν i , and ν a are the ionization and attachment rate coefficients. These swarm parameters are approximated based on the local field and they are calculated using the Bolsig + software or determined from experiments [22][23][24][25][26][27][28].
We assume that the superposition law holds and the velocities of particles can be approximated by adding the air velocity v air to the drift velocity. The drift-diffusion assumption for the velocity is valid for electrons as the time constant of momentum transfer is short compared to timescales of electric field variations in the streamer head [7,29]. Although the ions have a higher inertia than electrons, we included the drift-diffusion model for ions in order to provide accurate results at microsecond timescales. The model is justified since the ion velocities are smaller than that of electrons and the ion displacement in streamer timescales is negligible. The effect of eventual delay to relax to the drift-diffusion velocity does not have significance in long timescales.
Poisson's equation for electrostatics (1.4) is used for the calculation of the electric potential: where φ is the electric potential, 0 is the vacuum permittivity, and ρ c is the charge density. The electric field is the gradient of the electrical potential and is calculated through: Photoionization effects are modeled based on the approximation proposed by Luque et al [30] of the integral model of photoionization from Zhelezniak et al [31]: where I(r 2 ) is the emission of photons from differential volume at r 2 , and S ph (r 1 ) is the photoionization source term at r 1 . The method proposed by Luque et al [30] approximates the absorption function, f(|r 1 − r 2 |), by a series of exponential functions: where p O 2 is the partial pressure of molecular oxygen. The integral equation, (1.6), can then be calculated by solving a set of Helmholtz differential equations. The fitting parameters are taken from [32]. The set of Helmholtz equations are [32]: (1.8) with I(r) defined as: where p q /(p + p q ) is a quenching factor, ξ is the photionization efficiency, ν u is the electron impact excitation frequency of level u, and S i (r) is the electron impact ionization source term [32]. Finally, the photoionization source term is:

Numerical model
The general framework for implementing equations comes from Ansys Fluent which allows us to implement user-defined scalar (UDS) transport equations [33]. We used this capability to implement the drift-diffusion equations for charged particles, Poisson's equation for the electric potential and the Helmholtz equations for photoionization. The general form of the transport equation provided by Ansys Fluent is [33]: where, v is the velocity of the transport of the variable φ, D is the diffusion coefficient, and S φ is the source term. The drift-diffusion equations of charged particles have been implemented using the general format of the transport equation (2.1).
Poisson's equation and the Helmholtz equations have been implemented using the transport equation (2.1). We removed the unsteady and advection terms from (2.1) and replaced the source term with the source terms from (1.4) and (1.8). The diffusion term in (2.1) is discretized by a central difference scheme.

Spatial discretization
Ansys Fluent allows the use of several different schemes to calculate the advection flux. In our simulations, we used a second order upwind scheme and a slope limiter to bound the face value to the values of the cells straddling the face.
where n c 0 and n f are respectively the value of the variable at the upwind cell center and face center, and vector δr c 0 − f , connects the cell center to the face center. ζ is a scalar between 0 and 1 and is calculated as follows: where ζ is chosen such that the value of the n f would always be bounded between n c0 and n c1 . When ζ = 1 we have a 2nd order upwind scheme, and when ζ = 0 we have a 1st order upwind scheme.
The discretization of the diffusion flux is done using the central difference, a second order scheme. This scheme assumes a linear profile for the variable between the center of the two cells straddling the face.
The gradient calculation is carried out using the least square cell based method [34]. This method provides comparable accuracy with lower computation time than Green Gauss node based method [35,36]. It preserves a 1st order spatial accuracy in an irregular mesh and a 2nd order spatial accuracy in a regular mesh [35].

Time discretization
Ansys Fluent solves the UDS transport equations in an implicit way, and allows the implementation of user defined time integration schemes. In this paper, we present results from two different integration schemes. The first is the classical 2nd order backward difference scheme that provides accurate results in small time steps. The second is an implicit ELP scheme designed to allow the use of long time steps.
The 2nd order backward difference scheme is not accurate in high electric field regions especially if the time step size is large. To remedy this shortcoming, we implemented an implicit ELP integration scheme.
We have implemented an integration scheme for electrons and ions following [37], which uses an exact integration for the linear part describing charged particles growth and decay (ELP). The motivation behind this is the error in high electric field regions around streamer tip [1,4,38]. If the time step is not small enough, as compared to ionization timescale, the second order scheme introduces a relatively large error. The implicit ELP integration scheme integrates the source term due to direct electron impact ionization in an exact manner, and therefore does not suffer the inaccuracies in high electric field regions at large time steps.
To derive the ELP scheme, we write the equation describing the electron dynamics: n e = (ν i − ν a )n e + N e , (2.4) where n e is the electron density and N represents the terms due to photoionization, recombination and transport fluxes. We integrate the linear growth and decay term, (ν i − ν a )n e , in an exact manner and the other term, N, is integrated using an implicit 2nd order scheme. The derivation of the ELP integration scheme for electron density, n e , is similar to the one described in [37], which we modified to incorporate adaptive time stepping. The derivation of the scheme and the terms that are required to be passed to Fluent are summarized in appendix A. ELP schemes are not directly applicable to ions. Ion growth and decay dynamics include terms that are not linear w.r.t. ion densities, but linear to electron density.
(2.5) Figure 1. Contours of R p (2.7) variable near the streamer tip is shown on the left; the contours of the R p opt (2.8) variable is plotted on the right. The variable R p opt allows us to identify the streamer propagation direction and refine the mesh accordingly.
The ion dynamics has a large term due to electron impact ionization, and is solved by taking into account the ionization from the linear part of (2.4). The linear part of the growth and decay dynamics is integrated in an exact manner, and the other parts of the dynamics are integrated implicitly with second order accuracy. A detailed derivation of the integration scheme for ions is also provided in appendix A. The advantages of using the implicit ELP scheme over the second order backward difference scheme are presented in section 3. The spatial and temporal discretization presented above lead to a pseudo-linear system of equations: where the matrix A includes the coefficients of variables and vector b includes the source terms and the explicit contributions from previous time steps. Both A and b are dependant on the solution variables at the current time step. Fluent solves this system of equations in an iterative way, allowing A and b to update after each iteration. The algebraic multigrid method (AMG) has been used to solve the set of linearized equations at each iteration [39]. The solution convergency at the end of each time step is checked based on the scaled residuals defined by Ansys Fluent [33].

Adaptive mesh refinement
Due to the multiscale nature of streamer discharges [40], it is crucial for 3D streamer simulations to use adaptive mesh refinement. A typical length of charge separation in a streamer head is 4 μm which can be several orders of magnitude smaller than the size of the domain, and therefore without adaptive mesh refinement the computation time would be extremely high. In this work, we use the adaptive mesh refinement from Ansys Fluent that is based on the hanging node method [33].
An important aspect of adaptive mesh refinement is the criterion to choose the refinement region. Our criterion for refinement is based on the parameter: where S ph 2 is the 2nd term of the photoionization source (1.10) and |E| is the electric field magnitude. The contour plots of R p is shown in figure 1 (left). Including both photoionization and electric field has two advantages. Firstly, the refiner does not refine the streamer body due to low electric field inside the streamer body. Secondly, it refines the high electric field region around the electrode only when there is high electron density near the electrode in the form of a patch or a streamer, otherwise it skips the regions around the electrode. Consequently, the high density mesh is focused in the streamer head and the region in front of it. The algorithm works as follows: • We evaluate the value of max(S ph 2 |E|) in the domain, and then we calculate the value of R p from (2.7) for every cell. • If the value of R p in a cell is larger than some pre-specified values {R p 1 , R p 2 , R p 3 }, then the cell's volume should be smaller than a pre-specified value {V p 1 , V p 2 , V p 3 }, otherwise the cell is considered large, and is marked for refinement. The mesh refinement layers in front of the streamer are shown by {p 1 , p 2 , p 3 }, and {V p 1 , V p 2 , V p 3 } are the corresponding maximum cell volume allowed in that refinement level.
that the maximum dimension of the cell is smaller than the ionization length (1/α), where α is the townsend coefficient of ionization.
The above algorithm leads to different refinement levels covering the streamer tip, resulting in a mesh density that gradually increases toward the streamer tip. Figure 2 shows the mesh around the streamer tip and the different mesh density levels. Figure 2 (left) shows that the refinement region fully covers the streamer tip. However, the refinement region also spreads Mesh around the streamer tip before (left) and after (right) using the optimizer algorithm based on R p opt (2.8) to limit mesh refinement to the region in front of the streamer. Contours of the charge density ρ c is also plotted to show the full coverage of the streamer tip. to the side of the streamer tip which is considered unnecessary, as the streamer is not propagating in that direction. To improve the algorithm, we have added a second criterion that prevents refinement of cells in these regions, which is based on the optimization parameter R p opt defined as: and is in fact the cosine of the angle between the gradient of the photoionization source term and electrons velocity vector. Figure 1 (right) shows the contours of R p opt around the streamer, and that the unnecessary mesh refinement on streamer sides can be avoided if we prevent refinement of cells where |R p opt | < 0.75. Figure 2 compares the mesh around the streamer tip before (left) and after (right) applying this improvement. We note that the above criteria do not depend on the direction of the electron velocity, the refinement algorithm is therefore applicable for both positive and negative streamers.
We also have an extra criterion that refines the streamer body up to 10 μm. The criterion is based on the positive ion density such that only cells having a density above a given threshold are chosen for refinement.
Finally, a criterion has been added to moderately refine the mesh at the vicinity of the electrode; it is simply based on the distance between the cell and the electrode.

Streamer model comparison
The accuracy of the model is verified by comparing the model results to other streamer codes. Bagheri et al [17] have performed a comparison of six streamer codes from different research groups. They compared the accuracy and performance of the codes with respect to three test cases. Here, we have simulated the third test case, which is an axisymmetric simulation of a streamer discharge in air. The configuration of the simulation domain is shown in figure 3. In order to enhance the electric field, an spherical Gaussian patch of positive ions with a max density of 5 × 10 18 m −3 and an e-folding radius of 0.4 mm is placed in the domain. The top plane electrode is at potential φ = 18.75 kV and the lower electrode is grounded. The transport and rate coefficients, and photoionization parameters are calculated in the same way as in [17]. It is to be noted, that we performed the simulation using a 3D model instead of an axisymmetric one. Furthermore, we adopt the same The results from our model is shown by the DTU-X ps-Y μm tag, where X and Y are respectively the time step size and the minimum mesh size in the streamer tip. The other results are from [17]. Table 1. The settings and the run times of the codes compared in [17]. The data in this table are taken from [17]. All the simulations were performed in a machine using Intel(R) Xeon(R) E3-1271 v3 CPU running at 3.6 GHz [17].  [17] to refer to different streamer codes either by the name of the institute or the country of origin. We performed three simulations with different time step and minimum mesh sizes. The first one uses a time step size of 2.5 ps and a minimum mesh size of 2.5 μm; it is tagged 'DTU-2 ps-2.5 μm' in figure 4. In the second simulation ('DTU-2 ps-5 μm'), we used a time step size of 2 ps, and a minimum mesh size of 5 μm. In the third simulation ('DTU-10 ps-5 μm'), the time step and the minimum mesh size were respectively 10 ps and 5 μm. Figure 3 (2nd and 3rd panels) shows the contour plots of electron density and electric field at 15 ns from 'DTU-2 ps-5 μm' simulation. In figure 4, we compare some important characteristics of streamer simulations with the study by Bagheri et al [17]. Figure 4 (top left) compares the maximum electric field at the tip of the streamer at different lengths. The maximum electric field from our model is in good agreement with the results from CWI and DE. The results from 'DTU-2 ps-2.5 μm' is almost identical to the results from the DE model. Interestingly, the result from 'DTU-10 ps-5 μm" simulation is closer to the results from 'DTU-2 ps-2.5 μm' than the result from 'DTU-2 ps-5 μm' simulation. The reason is that the 10 ps time step size corresponds to a CFL value in the streamer tip region that is close to 1, which leads to lower numerical diffusion and dispersion, and hence increases the accuracy of the solution. The comparison of the models showed that the results from our model are in agreement with other streamer codes, specifically with the codes from CWI and DE. Our code did exceptionally well in conserving the charges when using long time steps. This comparison demonstrate the accuracy of the model to simulate streamer simulation in steady air. In the next section, we simulate a long-duration electric discharge in unsteady air.
The settings and the run times of the DTU model and of the models from [17] are reported in tables 1 and 2. For DTU model, we performed the same simulation in an axisymmetric configuration with both fixed mesh and adaptive mesh refinement. It should be mentioned that implicit schemes require more computational time per time step; however, they are not limited by the CFL requirement or the dielectric relaxation time as opposed to explicit models, and can use longer time steps [41]. Specifically, as we will see in the next section, the implicit models can be very effective in pulsed streamers simulations by using long time steps between the streamer pulses.
In order to compare the performance of the time integration schemes, we repeated the previous simulations using the 2nd order backward difference scheme. Figure 5 compares the results of the implicit integration scheme to the results obtained using the implicit ELP scheme along with the results from CWI and DE models. Figure 5 (left) compares the streamer tip electric field vs streamer length. Figure 5 (right) shows the variation of L(t) − vt in time. The implicit and implicit ELP schemes perform similarly at the small time step size of 2 ps. However, increasing the time step size to 10 ps leads to an inaccurate solution for the 2nd order implicit integration scheme, while the implicit ELP scheme remains accurate.
To conclude the comparison, both schemes are stable when taking large time steps. However, inaccuracies develop when using the 2nd order implicit integration scheme. The implicit ELP integration scheme remains stable and accurate while using long time steps, which is essential to model longduration discharges. Furthermore, the implicit ELP scheme allows the solution to converge faster in each time step, leading to a lower computation cost.

Positive streamers exposed to a transverse wind
We modeled the behavior of positive streamers in the presence of lateral wind, based on an experiment performed in [14]. The experiment shows that the positive streamer corona when subjected to a transverse wind, tilts in the direction of the wind. Here, we have devised a simulation to capture the effect of the airflow, and investigate the main mechanisms that affect positive streamers in low speed airflow. This simulation demonstrates the ability of our model to simulate long discharges coupled with unsteady airflow, and also helps us understand the physical mechanism explaining the tilting.
First, we need to point out some peculiarities of this simulation: (a) The effect of the airflow on electric discharges reveals itself at the timescales corresponding to heavy particles, and therefore the duration of this simulation needs to be long (around 0.1 ms), which renders the simulation very challenging. However, the implicit ELP integration scheme allows us to take very large time steps and makes this simulation feasible. (c) In this simulation, we have streamers that reach the stagnation point where the electric field in their tip fades away, the propagation stops and the streamer vanishes. Previously, Pancheshnyi [43] and Starikovskiy [44], using a first order fluid model and local field and density approximation, simulated positive streamers that eventually came into stagnation. In their simulations, the streamer head became narrower while the electric field at the streamer tip grew indefinitely and finally the simulation could not continue [43,44]. We performed an axisymmetric simulation of a positive streamer using the local field and density approximation for ionization and experienced the same behavior, see figure 6. The electric field in the streamer tip grows unboundedly and eventually the simulation stops, similar to the behavior observed in [43,44]. However, the local field and density approximation of the ionization source term is not accurate at high fields, as a large contribution to ionization comes from fast, high energy electrons which will spread the ionization over a significant length. Therefore, we smooth the ionization source term in high electric field regions to mimic the non-local aspect of electron ionization. We solved the following equation to smooth the ionization source term: where L is the smoothing length scale, and S ni is the smoothed ionization source term, used in the drift-diffusion equations for electrons and positive ions (1.1) and (1.2). We have set the boundary condition at electrode surface for the smoothing equation to Neumann zero (∂S ni /∂n = 0). We can note that the total number of ionization events in the domain is not altered by this smoothing process. By integrating both sides of (4.1) over the entire domain, we have: Figure 6. This plot shows the evolution of the streamer tip electric field in an axisymmetric simulation without the airflow and using the same domain and parameters as the 3D simulation presented in figure 8. As the streamer propagates, it narrows down and the field in the streamer tip rises indefinitely and eventually the simulation stops, similar to the behavior described in [43,44]. and, since the integral of the diffusion term over the volume cancels. Figure 7 compares the ionization source term S i using the local field approximation in the streamer tip with the smoothed ionization source term S ni . The smoothing length scale, L, is limited to 4 μm in high field regions, and as shown in figure 7 (top), the smoothing process slightly changes the ionization source term near the tip, and as we get away from the streamer tip, the smoothing effects disappear. The smoothed ionization source term solves the issue of the unbounded electric field for a non-propagating streamer.
In figure 8, we show the simulation domain in a pointplane electrode configuration with the point anode at 8.5 kV electric potential. In order to experiment with wind speed magnitude coherent with the ones found near the tip of wind turbine blades, we set a lateral wind of velocity v air = 50 m s −1 , blowing transversely. To initiate the first streamer, we place a neutral Gaussian plasma patch of electrons and positive ions on the tip of the electrode to initiate the streamer. The maximum density of the patch is 10 14 m −3 and the e-folding radius is 0.14 mm. The charged particles boundary condition at the point electrode surface is set to absorb all the incoming charge fluxes, and no out flux.
The air streamlines around the electrode are shown in figure 9. Compressible Navier Stokes equations for an ideal gas along with shear-stress transport k-ω turbulent model [45] is used to calculate the airflow around the electrode.
The 8.5 kV electric potential applied to the point electrode leads to a maximum electric field of 110 kV cm −1 on the electrode tip. The streamer emerges from the plasma patch and propagates in the inhomogeneous field created by the anode while a uniform lateral wind is blowing to the streamer. During the initial phase after the streamer initiation, the streamer  tip electric field rises to around 190 kV cm −1 and the streamer velocity V s reaches 0.45 mm ns −1 (see figure 10). As the streamer length increases, the streamer thickness shrinks and the electric potential in the streamer tip decreases, such that there is a steady decline in the electric field and streamer velocity. The streamer propagation slows down till it reaches stagnation, at this point as seen in figure 10 (right) the electric field in the streamer tip fades out.
The plots of the electron density in figure 11, show the propagation of the streamer over time. The last panel in figure 11 shows the electron density at 50 μs, which is long after the streamer has stopped propagating and clearly shows the effect of the airflow and electron decay. The streamer reaches a maximum channel length of 3.15 cm. In order to compare the channel length to the one measured in experiments, we try to choose an experiment with a single streamer to avoid the streamer interaction effects. However, it seems most experiments in atmospheric pressure lead to a collection of streamers rather than single streamers [5]. Here, we have chosen an experiment from [46] that produces relatively few streamers, and hence the streamer interactions are low. In one of their experiments, they captured positive streamer discharges in a gas composition of 20% oxygen and 80% nitrogen at a pressure of 200 mbar. In the experiment, the electrode with an 11 kV potential initiated a few streamers that crossed a 16 cm gap. If we scale down the streamer length using the Townsend similarity law to account for the low pressure, we will have an equivalent streamer length of 3.2 cm in atmospheric pressure. This result in an electrode potential to streamer length ratio of 3.4 kV cm −1 , the corresponding ratio in our simulation is around 2.7 kV cm −1 , which is in good agreement. Our results are also in good agreement with the results from [44], where they report an average electric field inside the streamer of 3.25 kV cm −1 .
After the streamer stops propagating and the electric field in the streamer tip vanished, the simulation enters the quiescent phase. In this phase, the density of the charged particles created mostly by the 1st streamer are continuously decreasing. The main mechanisms removing charged particles are electron absorption by the anode, attachment, and ion-ion and electron-ion recombination. The total number of charged particles in the discharge gap is plotted in figure 12 (left), and we can notice a large decrease in the total number of particles after the 1st streamer. The electron density inside the   domain is reduced to a greater extent than the positive ions, mostly due to absorption by the anode. A similar behavior has been observed in [18]. The positive charge density remainder from the 1st streamer dampens the field around the electrode, inhibiting the emergence of a new streamer. Figure 12 (right), shows that the electric field near the electrode rises over time as the positive ions density is reduced and pushed further away from the anode. Note that the implicit ELP integration scheme was essential in performing this simulation, as we had to use very long time steps.
During the quiescent phase, we observe a very interesting effect. The plasma channel remaining from the initial streamer, remains attached to the electrode surface while being moved by the air, as shown in figure 13. The streamer remains attached  to the surface, notwithstanding the separation of the airflow from the electrode surface. The reason for this behavior is that the electrons inside the streamer channel are continuously moving toward the electrode, and they ionize the air in their path and maintain the channel (see figure 13). Therefore, the streamer channel remains constantly connected to the electrode surface. It should be noted that the electric field inside the streamer channel near the electrode, due to charge removal, is high enough to provide sufficient ionization, which is shown in the ionization source plots in figure 13 (bottom row). The significance of this sustained plasma channel connection will be understood in the next section, where we discuss the initiation of a new streamer from the remaining charges in this plasma channel.
During the quiescent phase the electric field keeps rising as shown in figure 12 (right), and eventually it will reach the Laplacian field in the absence of screening charges allowing the possible initiation of a 2nd streamer. However, the rise in the electric field was slow and in order to speed up the simulation, after 50 μs, we increased artificially the electrode potential to 10.5 kV with a rise time of 100 ns. Then the charge density remaining from the 1st streamer and transported by the wind over the surface of the electrode, developed into a new streamer. The plots of electron density in figure 14 clearly Figure 15. Contours plots of electron density n e and electric field |E| for the 2nd streamer at different time instants. Streamer is initiated from the trailing side of the electrode, and starts following the pre-ionized channel left by the 1st streamer. Later, it changes its direction and starts to follow the field lines as the pre-ionized density in the channel lowers. The streamer reaches stagnation at a length of 3.75 cm, shown in the 3rd plot, and during the next 50 μs it continuously tilts in the direction of the wind. The last plot shows the tilted streamer channel which is still attached to the electrode at 104 μs. show the initiation of a streamer from the residual charges of the 1st streamer.
Initially, the streamer starts to follow the pre-ionized path remaining from the 1st streamer, as positive streamers tend the follow the region with a higher free electron density [47]. Later, as the streamer move away from the electrode, the electron density in the pre-ionized path decreases and the streamer starts to change direction and follows the field lines instead. In figure 15, the plots of electron density and electric field clearly demonstrate the curved path of the 2nd streamer. In the final 1 cm of the 2nd streamer, it shifts slightly toward left due to the field created by the residual charges of the 1st streamer. The new streamer propagates for 3.75 cm, which is longer than the 1st streamer as the potential of the electrode has been increased from 8.5 kV to 10.5 kV. Figure 16, shows the evolution of the electric field and streamer velocity for the 2nd streamer. The maximum length and the propagation time of the 2nd streamer is longer than the 1st one due to a higher electrode potential. Similarly to the 1st streamer, the 2nd streamer stops propagating, and another quiescent phase starts with a sharp reduction of the charged particle densities and a steady increase in electric field near the electrode. We modeled this phase until 104 μs, which show an even more pronounced tilting of the 2nd streamer plasma channel in the direction of the wind. The channel also remains attached to the electrode.

Conclusion
A 3D model is presented to couple the dynamics of streamer discharges and unsteady airflow. The model allows us to investigate the effects of unsteady airflow on electric discharge characteristics at different timescales. The effects of unsteady airflow on streamer discharges through charge redistribution manifest in longer timescales than the typical timescales of a streamer discharge. We therefore implemented a model that allows large time step sizes by using an implicit ELP time integration scheme. The model also uses an unstructured mesh allowing the incorporation of solid bodies with complex geometries. Comparing the results of the model with a test case from literature, shows a good agreement. The ability of the model to couple streamer discharges and unsteady airflow was demonstrated by a simulation of the initiation and termination of two successive positive streamers in transverse airflow.
The results of the simulation can be summarized by several points: • During its propagation, the 1st positive streamer had only a small tilting in the direction of the wind due to its short duration. It shows that the effect of the low speed airflow on the electric discharges manifests itself on longer timescales than the ones of the streamer propagation driven by electron timescales. • After the stagnation of the 1st streamer, the streamer channel is moved by the air. However, it remains attached to the surface of the electrode, even though the air streamlines suggest that the streamer channel should separate from the electrode surface. This is an intrinsic characteristic of positive streamers and is due to the flow of electrons from the streamer channel toward the electrode and subsequent ionization inside the plasma channel near the electrode surface. It is uncertain that negative streamers remain similarly attached since electrons move away from the electrode. This simulation shows that the effects of the low speed airflow on positive streamers mainly result from the charged particle redistribution in the domain which will affect the subsequent streamer discharges.

Acknowledgments
This project has received funding from the European Unions Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement 722337. The computational resources were provided by DTU Computing Center (DCC).

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.

Appendix A
This appendix gives a derivation of the integration scheme for charged particles and describe how it is implemented in Ansys Fluent. The derivation of the integration scheme for electrons is similar to [37]. In the following, we show the modifications we brought to the scheme to allow a varying time step size and its implementation in Ansys Fluent. The derivation of the integration scheme for ions also follows [37] with a higher level of complexity coming from the fact that the linear term in ion dynamics is proportional to electron density and not the ion density.

A.1. 2nd order implicit ELP integration scheme for electrons
Starting from the Taylor series expansion of the electron density n e n+1 at the time step n + 1 given by: where λ n+ 1 2 = (ν i n + ν i n+1 − ν a n − ν a n+1 )/2 is the effective ionization frequency averaged between steps n and n + 1, and N e n represents the other source terms, which include the variations due to fluxes, photoionization and recombination. The first term is the familiar exponential growth term and the second term can be simplified by changing the order of summation such that we have: with Q j defined as [37] by: In order to have a 2nd order implicit ELP integration scheme the integration scheme needs to equalize the Taylor series up to the 2nd order terms: Assuming the following form for the integration scheme: n e n+1 = n e n e λ n+ 1 2 Δt n + γΔt n N e n+1 + βΔt n N e n (A. 8) and expanding N e n+1 in Taylor series, we have: Δt n − γ (A.10) The integration scheme to calculate n e n+1 now depends on N n , which includes the flux terms and the source terms from the time step n. In order to have an implicit solver, we derive a similar expression for n e n−1 and combine the two expressions to remove the dependency on N n . We proceed as follows: The Taylor series expansion of the electron density at previous time step, n e n−1 : after following a similar derivation path as that of n e n+1 , we have: where λ n− 1 2 = (ν i n + ν i n−1 − ν a n − ν a n−1 )/2 is the effective ionization frequency averaged between steps n and n − 1, and Δt n−1 , with Q j defined as: We combine the approximations of n e n+1 and n e n−1 in such a way to remove the explicit dependence on N e n . The expressions of the n e n+1 and n e n−1 are: where the ionization and attachment source term are removed from the equation, and the unsteady term, ∂n e ∂t , is approximated by N e n+1 .
At the initial time step we do not have the value of n e n−1 . Therefore, we have to integrate the terms represented by N e n+1 using a first order scheme. The scheme to calculate the n e n+1 would be: Ansys Fluent expects the unsteady term in the UDS transport equations to be decomposed into implicit and explicit terms, they are respectively called 'apu' and 'su' in the Fluent manual [48], and are given by: .
We have limited the use of the implicit ELP integration scheme to the regions of the domain where the electric field magnitude is strictly above the breakdown field, i.e. where the electron density grows exponentially. In the other regions, the integration is done by a classical 2nd order backward difference scheme. The two schemes will provide the same level of accuracy near the breakdown field, therefore there is no need for a special treatment of the interface between the cells using the exponential and the 2nd order backward difference schemes.

A.2. 2nd order implicit ELP integration scheme for ions
We now describe the derivation of the integration scheme for the ions. Starting from the Taylor series expansion of the ions density at the next time step, n i n+1 : where the jth time derivative of the ions density is given by: There are four different terms in the rhs The first term is the ion density in the current time step. The second term is a closed form analytical expression and needs no further simplification. The fourth term also needs no further simplification. The third term is simplified by changing the order of summation: with Q j (x) defined as [37] by: Then n i n+1 reads: In order to have a 2nd order implicit ELP integration scheme the integration scheme needs to equalize the Taylor series up to the 2nd order terms. Assuming the following form of the integration scheme: n i n+1 = n i n + n e n ν i Δt n − 1 + Δt 2 n (ζN e n+1 + ηN e n ) + Δt n (γ i N i n+1 + β i N i n ).
(A. 34) Expanding the N i n+1 and N e n+1 in Taylor series we have: N i n+1 = ∇ · (D p ∇n p ) + ∇ · (v air + μ p En p ) − S r e−p − S r n−p , (A.55) where the ionization source term has been removed from the equation, and the unsteady term, ∂n p /∂t, is approximated by N i n+1 (A.51).
In the initial time step there is no information about the previous time step so we have to use the first order method: Δt n−1 (γ i Δt n + A 1 γ i Δt n−1 ) −1 .
(A.59) Similar to electrons, the implicit ELP scheme for ions is applied to the regions of the domain where the electric field magnitude is strictly above the breakdown field. In the other regions, the integration is done by a classical 2nd order backward difference scheme.