Global model simulation of OH production in pulsed-DC atmospheric pressure helium-air plasma jets

1 Introduction

The fast growing application of complex water-containing plasma systems in various fields, such as material processing [1, 2], bio-medical applications [3], surface modification [4], puts new requirements on the level of their understanding.

Special attention in the last decade was put on production of atmospheric pressure plasmas rich with OH radicals, as OH is one of the most reactive species in such systems [5]. Density and the underpinning production pathways of OH radicals can be determined both by the direct measurements [6, 7, 8] and by means of numerical modelling. Two groups of models have been mainly applied, global (volume averaged or zero-dimensional) models [9, 10, 11] and multidimensional fluid (1D or 2D) models [12, 13, 14, 15]. Despite neglecting the spatial derivatives, global model offers an opportunity of a fast assessment of chemical composition of complex plasmas, and determination of the main production pathways for all plasma species. On the other hand, modelling of plasmas formed in gas mixtures which contain humid air and/or water imposes restrictions on multidimensional models, mainly due to extremely large computational time and hardware resources.

In a recent paper [14] the result of a numerical study of the effect of water content on OH production in a pulsed-dc plasma jet in helium was presented. The main production and loss mechanisms of OH were analysed for 1%, 0.1% and 0.01% H₂O contents in air. Spatial and temporal distributions of air mole fraction and streamer discharge parameters (the electron density and temperature, electric field and total ionization rates) were calculated by means of fluid modelling. The plasma chemistry set consists of 21 different species and 71 elementary reactions. This study is further extended to examination of the effects of pulse voltage rising time on discharge characteristic and plasma composition of a same gas mixture, especially on O and OH densities distributions [15].

On the other hand, global models have been used for an estimation of the chemical composition of atmospheric pressure plasmas created in mixtures containing helium as dominant gas with water [9], with real air [10], or with oxygen and water or humid-air [11]. Papers [16, 17] present an analysis of the influence of the humidity level on chemistry of reactive species in RF-driven atmospheric-pressure helium-oxygen mixture (0.5% of O₂) plasmas.

The aim of this paper is to test the applicability of global model in a simulation of a pulsed-dc plasma jet in helium, described in [14]. The obtained results, particularly density and the main chemical processes responsible for creation and destruction of OH radicals, will be compared with the results of a fluid model [14]. With this idea, we have created a zero-dimensional code by using a comprehensive list of the chemical processes with rate coefficients from Appendix and supplementary data files of papers [16, 17]. Model comprises a reaction scheme with 1425 reactions for kinetics of 68 species. The results of the fluid simulation in [14] have been used as input parameters, and as a basis for testing of our model results. A similar comparison between fluid and global model has been made for the modelling of a streamer dielectric barrier discharge [18] and for atmospheric pressure radio-frequency plasmas [19].

The details of global model and the way how input parameters were adopted from Ref. [14] are summarized in Section 2. The analysis of our simulation’s results and their comparison with results obtained by fluid modelling follows in Section 3, and finally concluding remarks are given in Section 4.

2 Model description

The global model is taken in a reduced form and only the equation that describes the balance of particle number density is included. In addition, we assume that analysed plasma system is spatially homogeneous leading to the governing equations which include only the external terms of generation and loss in particle collisions [20], in the following form:

In the previous equation n_α denotes concentration of species α, and the first and the second terms denote sum of the generation and loss terms of species α through specific volumetric chemical reactions and include two-body and three body processes. The rate coefficients for electron-impact processes are given as a function of the electron temperature T_e in [eV] (mean electron energy) [16, 17].

Numerous global models [20], including the present work, are made for a uniform infinite system by using the mean energy as a parameter connecting different sources of input data. Determination of the mean electron energy is a very important issue in global modelling and implies solving of the electron power balance equation, simultaneously with particle balance equation for each plasma species. The mean electron energy is then assumed to be ϵ = 3/2 n_eT_e (with T_e in [eV]) based on the assumption of Maxwell-Boltzmann electron energy distribution function (EEDF). Another approach to obtain T_e is to use results of a fluid model applied to the same plasma system, but with simplified chemistry set, as an input to the global model [12, 16]. We chose the former approach for our work, with idea to compare results of the plasma chemical part of the global and fluid models, with the same discharge parameters (i.e. electron density and electron temperature).

The simulation domain is fully described in papers [14, 15] and Figure 1 shows rough sketch of schematics of the plasma jet. When the positive voltage pulse is applied, starting from a needle tip area (anode), the discharge zone expands 2.3mm radially to the dielectric wall, and 3mm axially to the quartz plate. The grounded metal plate is assumed to be behind the quartz plate. The working gas is pure helium flowing into domain filled with humid air.

Figure 1

(Colour online) Schematic of simulated plasma jet from [14] with marked position and values of: a) electron temperature (red contours), b) maxima of OH density (magenta contours), c) higher electron density (brown contour), and d) lower electron densities (aqua blue contours), estimated from 2D model results [14]

The estimation of input parameter from 2D simulation results is very important for application of global model, which basically deals with volume average quantities. Since the results in [14] were present in form from which volume average electron density and electron temperature cannot be correctly extracted, we decided to cover the whole interval of n_e and T_e. The spatial distributions of n_e and T_e presented in Figures 3a and 3b [14] are marked in Figure 1, together with position of maximal OH production. It is obvious that figure represents final stage of plasma plume (or streamer head) growth, with n_e and T_e distribution in that moment. According to figure 1, the region of maximal OH production overlaps with region where electron density is around 1.5·10¹³cm⁻³ and electron temperature between 6eV and 9eV. Bearing in mind fast streamer development and non-homogeneity of discharge process, [19, 21, 22, 23, 24], the early phase of streamer growth, with n_e and T_e distribution which differ from those in Figure 1, needs to be taken in consideration. So, we chose to run our code with the electron density’s values between 1.0·10¹³cm⁻³ and 6.0·10¹³cm⁻³, and with the range of electron temperatures from 1.5eV to 9.0eV, estimated from the spatial distributions of n_e and T_e presented in Figures 3a and 3b [14] and in Figure 1. Another important input parameter for modelling of atmospheric pressure plasma jet expanding in free air is the helium-air mixing ratio in the jet zone. In our calculation we estimated that, according to spatial overlapping of the position of OH density maxima and region with n_e ≥ 1.5·10¹¹cm⁻³ in Figure 1 with the air molar fraction distribution in Figure 1b in [14], the value around 5% can be taken as the maximum of the air amount in zone of plume propagation. On the other hand, the results of numerical simulation in recent paper [24] showed that the air amount in helium higher than 0.5% slows down streamer development in He/air mixture, through lowering of the electron energy peak value. Based on these results, we limited the air mole fraction to 0.01 (1% of air in helium). Finally, the input parameters in our simulation are as follows: air mole fraction in plasma zone is 1%, H₂O contents in air are 1%, 0.1% and 0.01%, and peak values of N₂ and O₂ are (0.78, 0.21), (0.788, 0.211), (0.7888, 0.2111), respectively.

The system of rate-equations is solved by MATLAB ODE15s solver, with relative and absolute tolerances equal to 10⁻¹² and 10⁻⁶, respectively. The pulse duration was 30ns and the time-step of 0.3ns was chosen. Similar modelling procedure was applied in case of a surface microdischarge in humid air at atmospheric pressure, up to 1000s [11]. The numerical procedure was tested by applying ODE45 solver for non-stiff system of differential equation, but the final concentrations were similar as in the calculation with ODE15s solver.

3 Results and discussion

We have started analysis of global model results with testing the sensitivity of global model results on the air content in helium. We made calculation with 5% and 1% of air in helium, in both case with 1% of the H₂O in air. All results presented in this section were obtained at 30ns after the pulse application. The results of calculation are presented in Figure 2 together with OH density interval evaluated from minimal and maximal value of spatial distribution in Figure 5a from [14]. It is obvious that the global model in the case with 5% of air (upper part of Figure 2) gives OH densities which are up to two orders of magnitude higher than the peak value in 2D simulation [14]. Only at the lowest values of n_e and T_e, our results lie in the interval of OH densities predicted by the fluid model in [14]. In the case with 1% of air and with the same amount of water, calculated curves intersect OH interval from [14] for all three values of n_e, but still are over targeted interval at the higher electron temperatures. In addition to above presented, this calculation additionally supports choice of 1% of air as input parameter.

Figure 2

(Colour online) OH density as a function of T_e determined by global model with: 1% of air (full lines with symbols); and 5% of air (dashed line with symbol) compared on interval of [OH] from [14] (parallel dashed lines) for 1% of H₂O in air

Observed disagreement between results of global and fluid modelling is to some extent expected for the case of nanoseconds pulsed discharge, where streamer breakdown occurs. The growth of streamer is associated with rapidly development of ionization front, caused by the large gradient of electric field which produces hot electrons with T_e up to 9eV [14, 15]. Indeed, this rapid development imposes some restrictions even on the fluid models, as have been discussed in papers [21, 22]. On the other hand, the data point in Figure 2 was calculated with assumption of the constant electron density at given electron temperature, in whole modelled domain. So, our global simulation does not include temporal development of n_e and T_e, and comparison with results of fluid simulation [14] is noted straightforward and quantitative, but more qualitative.

As it is noted above, in fluid simulation [14] the water fractions in air were 1%, 0.1% and 0.01%, and the calculated peak OH densities, from the Figure 5d of paper [14], are 6.3·10¹²cm⁻³, 7.0·10¹¹cm⁻³ and 7.0·10¹⁰cm⁻³, respectively. As a test, we have run our code with these H₂O fractions, but only for n_e = 3.0·10¹³cm⁻³. The Figure 3 shows that in all three cases the global model predictions of OH density are over interval calculated by the fluid model [14], for the electron temperatures higher than 3 eV. Also, the results of our simulation, for given H₂O amounts gives similar dependence of [OH] on T_e, with values which decrease by an order or magnitude with decreasing H₂O content.

Figure 3

(Colour online) OH density as a function of T_e obtained by global model with 1% of air for different H₂O fractions (full lines with symbols) compared to results of fluid model[14] (parallel dashed lines). Vertical arrows connect compared data sets

Calculated results with n_e = 1·10¹³cm⁻³ and n_e = 6·10¹³cm⁻³ are symmetrically positioned below and above that for n_e = 3·10¹³cm⁻³ in Figure 3, as is expected from dependencies in Figure 2, and they are not shown for clarity of presentation. So, one can conclude that overlapping of our results and results of fluid modelling [14], with 0.1% and 0.01% H₂O will be as in the case of 1% of H₂O shown in Figure 2, and we will proceed further analysis only for 1% of H₂O in air.

Another important task in modelling of atmospheric pressure plasmas generated in mixtures which contain water is revealing the main production and losses of OH radicals, as is done in [5, 8, 13, 14]. In order to reveal the main production and loss mechanisms, we have calculated the percentage contribution of all processes which determined OH radical’s kinetics, according the reactions list in [16, 17], with n_e = 6·10¹³cm⁻³ and for the T_e interval 1.5eV to 9eV. The results are presented in Table 2, and in Figure 4 and Figure 5. The percentage contributions in the case of n_e = 1·10¹³cm⁻³ are not presented, but they are discussed in text. The Table 2 lists the main OH production and loss processes, which are labelled with R1-R5 and R6-R12, respectively.

Figure 4

Contribution of different processes from Table 2 on OH production as a function of the electron temperature for n_e = 6·10¹³cm⁻³

Figure 5

Contribution of different processes from Table 2 to OH destruction as a function of the electron temperature for n_e = 6·10¹³cm⁻³

Compared to results of fluid model [14], processes R2 and R4 remain dominant in OH production in case of the global model, unlike the electron-impact dissociation of H₂O with contribution that is now several orders of magnitude lower and is omitted from Table 2. Similar results are obtained in [4], where contribution of electron induced water dissociation to overall OH production is about 2%.

Electron impact dissociation of water is recognized as one of dominant mechanisms of OH production in plasmas with water content around several thousand ppm (particles per million) [5]. Furthermore, global model reveals additional processes of OH production (R1, R3 and R5) which were not included in the list of processes considered in fluid models [13, 14]. Dissociation of H₂O by O(¹S) is of particular interest, since it carries more than half of the OH production at low electron temperature, and around 7% at 9 eV, as it is shown in Figure 4.

On the other hand, the contribution of process R4 on OH production significantly increases at higher T_e, mostly caused by increased production of H₂O⁺ ions, since the rate coefficient for H₂O⁺ recombination with electrons (Table 2) does not depend on the electron temperature [15]. According to our results, the production of H₂O⁺ ions dominantly goes through charge transfer reactions with nitrogen ions (N⁺, N2+ and N4+ ), for the entire T_e interval from 1.5eV-9eV. Further analysis reveals that nitrogen ions production pathways are based on processes of Penning ionization with He metastable and on He⁺ ions conversion. So, as a final conclusion, an increase of contribution of process R4 in the OH production (Figure 4), and in the total OH density (Figures 2, 3) with T_e, is related to the increase of electron-impact rate coefficients, given as function of T_e [16, 17] (in Arrhenius parametric form). The calculations results with low electron density, n_e = 1·10¹³cm⁻³ show that contributions of processes R4 and R5 are approximately reduced for factor 2.5 in whole T_e interval, and now processes R1 and R4 (booth including H₂O⁺ ions) carry around 67% of OH production at 9eV, while R2 and R3 dominate at low T_e.

Figure 5 shows an analysis of the percentage contribution of processes R6-R10 on destruction of OH radicals. Our simulation reveals that the electron dissociation of OH (R8) dominates at low T_e, mainly caused by the fact that the electron density is kept constant during 30 ns of simulation period, despite temporal growth during streamer development [14, 21, 22]. On the other hand, the processes involving positive ions He⁺ and N⁺ start to dominate in OH destruction at higher T_e, above 7eV [8]. This is a direct consequence of an increased production of positive ions through electron impact and Penning ionization. The quenching of OH by N and O atoms (R9 and R10) was stated in [5] as an important loss mechanism, and has around of 20% contribution according to our results.

The percentage contributions of destruction processes shown in Figure 5 are slightly affected by decreasing of electron density. The calculation with n_e = 1·10¹³cm⁻³ reveals that relative contributions of R8 in this case are approximately equal at low T_e and differ within 2% at higher T_e, compared to values in Figure 5. The net contribution of processes R6 and R7 is now few percent lower at higher T_e, with decreasing of R6 contribution and increasing of R7 contribution within 5%. The three-body recombination reaction H + OH + M → M + H₂O was stated in [5, 13, 14] as the main OH loss mechanism. According to results of our simulation, its contribution is approximately two orders of magnitude lower then contribution of R6-R8, in the whole modelled range of electron temperatures. Except for differences in kinetic scheme and used rate coefficients for this process [14, 16], observed disagreement may be explained by higher density of He⁺ and N⁺ ions according to our simulation, since global model calculation, generally, could give higher ion densities.

The question that is imposing after presented analysis is: whether the results of our global model may be improved? Having in mind that low temperature atmospheric pressure plasmas are strongly non-equilibrium, it is obvious that the choice of electron energy distribution functions (EEDF) would affects the kinetics of species [25, 26, 27]. Particularly, the non-equilibrium EEDF for higher energies will be quite different from equilibrium Maxwell-Boltzmann (MB) EEDF, due to considerable inelastic losses in processes of electronic excitation and ionization of helium gas. The rate coefficients for electronic collision calculated with assumption of MB EEDF and with non-equilibrium one, by solving of Boltzmann equation (BE), may be different by several orders of magnitude and lead to quite different chemical composition of modelled plasmas [25, 26, 27]. Some of rate coefficients in our calculation, taken from [16, 17] were expressed as a parametric function of electron temperature T_e in [eV], with assumption of MB EEDF. The results presented in Figures 2-5, were calculated in this manner. On the other hand, with idea to test the influence of non-equilibrium EEDF on global model results, the rate coefficients for electron-molecules processes were obtained by solving Boltzmann equation with the two-term approximation solver BOLSIG+ [28]. In BOLSIG+ calculation, we have used cross section data from the MORGAN database [29] for He, O₂, N₂ and H₂O as an input, for 1% of air in mixture and 1% of H₂O in air (78% of N₂ and 21% of O₂). Figure 6 shows the test results, through the impact of non-equilibrium rate coefficients on OH densities. It is obvious that BE curves are significantly shifted to the lower OH densities, and they now approach to OH interval from [14] at higher T_e, especially for n_e = 6·10¹³cm⁻³. Decreasing of OH density in BE case is a direct consequence of a drop of H₂O⁺ density, since its production is now reduced by using of lower rate coefficient.

Figure 6

(Colour online) Influence of non-equilibrium EEDF on global model results. The OH densities drops are marked by a black arrow. Dotted curves are calculated with the non-equilibrium rate coefficients (BE EEDF)

Figure 7

The plasma composition calculated by means of global model with 1% of air content, with the water content is given as percentage of the air: 1% (black bars) and 0.01% (white bars)

Another fact that would be stressed is concerned on choice of electron-impact rate coefficients. The fundamental characteristic of a streamer discharge is formation of a streamer head, an area with increased electric field strength, in front of which electron temperature reaches very high values (around 9 eV) and ionization frequency is maximized, as can be seen in Figures 3b, c, d in ref. [14] and in papers [21, 22]. So, adopting of the rate coefficients from [16, 17], where they were used for T_e up to 3eV, in the modelling of streamer-like discharge with T_e up to 9eV, is questionable, as was stated in literature [5, 20]. And finally, in Eq. (1), which describes time evolution of particle density, diffusion loss and loss by the flux of particles directed to the wall are omitted. Taking into consideration very short pulse duration of 30ns and atmospheric pressure, diffusion can be negligible [18]. From Figure 1a in [14] characteristic diffusion length for cylindrical geometry is estimated to be 1/λ² = 2.235·10⁶m⁻² (for r = 2.25mm and z = 3mm), and with OH diffusion coefficient in He, D_OH= 0.87·10⁻⁴ m²/s [30] we obtain diffusion loss frequency ≈ 200 1/s. Compared to dominant destruction mechanisms R6-R10, OH loss due to diffusion is four orders of magnitude lower.

Why use a global model? Despite a large number of different species and comprehensive set of processes that determine kinetics of different species, numerical solving of the system of Eq. (1) is not time-consuming and the calculation time for a single data point in Figures 2-6 is shorter than 10 seconds. The second advantage of global modelling is reflected through fast assessment of the plasma’s chemical composition. In each run of code, we obtained densities of all 68 species included in Table 1, and we are able to investigate the changes in chemical composition induced by different input parameters. So, global model allows us to assess the densities of many reactive oxygen and nitrogen species, as shown in Figure 7. The comparison of plasma’s chemical compositions for two different amounts of water in air (1% vs 0.01%) is also presented in Figure 7. The main difference arises in density of water-based species, with concentrations that are several orders of magnitude lower in case with 0.01% of H₂O in air. The overall amount of oxygen and nitrogen reactive species is not affected by changes in the water content.

4 Conclusion

According to presented results, the intrinsic approximation of spatial homogeneity in global model imposes a significant limitation on its application in simulation of dynamic systems, such is a pulsed plasma jet. Rapid development of streamer (in form of plasma plumes or bullets) in space and time, ultimately requests using of higher order models, able to correctly describe streamer parameters [21, 22]. On the other hand, global model’s advantages, such as short computational time, comprehensive list of chemical reactions, fast assessment of chemical composition, fast determination of underlying chemistry pathways for each plasma species, make global model still very useful for modelling of atmospheric pressure plasmas. For those reasons such models are helpful in industry.

In making comparison of our results, we have employed data from the literature for a fluid model and assumed that it is essentially correct, in an attempt to test the applicability of the global model. In doing so one has to be aware that fluid models suffer from their own limitations, including the treatment of non-locality, spatially inhomogeneous plasma production, assumptions about the energy distribution function and using of the swarm data and the need to simultaneously solve a large number of equations thereby necessitating a simpler chemical model [21, 22]. All of these reasons could also be regarded as sources of discrepancy. We have however, taken information from the fluid model as initial conditions for our model in order to put comparisons on equal footing. Making comparison between a streamer model and a global model is particularly difficult due to small radial dimension of plasma while spatial diffusion would extend profiles of long lived species over a larger volume.