Introduction

Atmospheric-pressure non-thermal plasmas, providing formation and reaction of radicals in ambient condition, has recently received a growing attention in manufacturing processes and chemistry applications such as surface modification and material synthesis [1,2,3,4], conversion of waste and natural gases [5, 6], and N2 fixation into ammonia or NOx [7, 8]. However, the energy efficiency of these plasma processes remains a key bottleneck against its wider industrial acceptance. Therefore, improving the selectivity and energy efficiency related to them has been the main focus of research studies carried out in this area. Increasing the energy efficiency and selectivity of product species can be achieved by optimisation of plasma reactors and making synergy between plasma systems and other previously common process approaches.

Natural gas is an alternative to petroleum that can be considered as a bridge into other renewable and sustainable fuel or chemical sources. Largely composed of methane, natural gas production is steadily increasing due to the new, previously inaccessible unconventional sources, such as methane hydrates or shale gas. As a principal component of natural gas, methane (CH4) is the important hydrocarbon feedstock for the synthesis of fuels and chemicals. Non-thermal plasma has exhibited promising potential in the methane conversion [9]. Methane activation in discharge plasma is a rapidly expanding research topic in both science and engineering. This includes life cycle assessment and eco-efficiency analysis at industrial scale in process engineering field [10, 11], and making synergy between plasma and catalysts to control the selectivity of products, product yield and conversion rate in chemistry field [12,13,14,15]. In plasma physics field, the focus is on the design of new plasma reactors and application of different discharge types for the purpose of controlling and tuning the reaction mechanism, heat transfer rate, and, more important, increasing energy efficiency [16,17,18,19,20,21]. The methane conversion processes are classified in two main areas: (i) oxidative conversion to syngas or \(\text{C}_1\)-oxygenates and (ii) non-oxidative conversion to C2H4 and \(\text{C}_3\text{H}_6\) as the desired products.

Oxidative Coupling of Methane

Tu and Whitehead developed an alternating-current (AC) gliding arc (GA) reactor to the co-generation of syngas and value-added carbon nanomaterials by plasma dry reforming of methane [22]. Bie et al. studied the dominated pathways for methane conversion into oxygenates and syngas in atmospheric pressure DBD using a 1D model [23]. Park et al. explored conversions of \(\text{CH}_4\) and \(\text{CO}_2\), energy efficiency, and syngas generation costs in a \(\text{CH}_4/\text{CO}_2\) AC-pulsed tornado GA plasma [24]. Indarto et al. achieved high yield production of methanol by a single-stage plasma reactor from non-catalytic conversions of \(\text{CH}_4\) and \(\text{NO}_2\) [25]. Lefkowitz et al. used the schlieren imaging and optical emission spectra to reveal the evolution characteristics during and between the discharge pulses for the nanosecond repetitively pulsed discharge in a \(\text{CH}_4\)/air mixture [26]. They observed that the flame-development time for methane/air mixtures in the flow tube depends on the pulse repetition frequency and total ignition energy. Levko and Raja studied the active species production by the \(\text{CH}_4\)/air streamer discharge and concluded that steamers initiated by the positive and negative trigger voltage cannot be considered symmetrical [27]. More recently, Qian et al. presented a fluid modelling study of the radical species generation by a positive streamer in a dense \(\text{CH}_4\)/air DBD reactor [28].

Non-oxidative Coupling of Methane

Chiremba et al. [29] used \(\text{BaTiO}_3\) pellets as packing material in a DBD reactor for direct nonoxidative conversion of methane. The methane conversion increased with the specific energy density and it was dependent on plasma power and flow rate, but the \(\text{H}_2\) selectivity did not experience a significant change by plasma power and flow rate variations. Also, when the specific energy density value was lower than 100 kJ/L, ethylene was the major product among all \(\text{C}_2\) species. Bogaerts et al. investigated influences of \(\text{N}_2\) dilution on the \(\text{CH}_4\) conversion into \(\text{H}_2\) in a \(\text{CH}_4/\text{N}_2\) dielectric barrier discharge (DBD) reactor [30]. Fridman et al. applied nanosecond pulses both in a DBD reactor and in an atmospheric pressure glow discharge (APGD) reactors, and their results indicated that both discharges provided \(\text{CH}_4\) conversion [31]. Their research illustrated that the low temperature discharge plasma can provide new approaches to highly efficient \(\text{CH}_4\) conversion.

For methane activation, alternating current (AC) plasma generates \(\text{CH}_x\) and H radicals, for which their successive interaction leads to the formation of \(\text{C}_2\) hydrocarbons, hydrogen, and C-containing deposits onto the electrodes and reactor walls. The change of discharge mode to pulsed plasma operation shifts the distribution of products. In comparison with AC power sources, nanosecond pulsed high voltage supplies show a higher energy efficiency and better conversion performance in methane plasma conversion processes [12, 16, 32]. The plasmas created with these power sources are more stable, homogenous, and chemically active. They produce fast ionization waves which propagate with velocities of up to \(10^{10}\) cm/s throughout the discharge gap. This fast ionization waves makes the plasma more uniform. In the atmospheric pressure range, the nano-second pulsed discharges appear in a glow mode and/or powerful sparks regime. Also, they are completely non-thermal and their transition to arc regime is restricted by pulse off-time interval. The carbon deposition and coke formation in plasmas driven by nano-second power sources are low and heat generation which leads to overheating is mitigated. In addition, these pulsed plasmas have a higher quenching rate which suppresses the subsequent dehydrogenation reactions. Therefore, they can be considered as a strongly effective approach to achieving high selectivity to saturated hydrocarbons. For these plasmas energy delivering rate is high and the amplitude of the voltage pulse is bigger than the breakdown threshold. These features cause that deposited electrical energy, rather than going into gas heating, is directly channelled into the electron-induced chemical reactions, the excitation of high-threshold electronic states of molecules and on their dissociation and ionization, and is spent to increase vibrational temperatures which are beneficial when catalyst materials are coupled with plasma. Also, due to a high reduced electric field, they can generate high densities of charged particles at low values of specific deposited energy, which leads to higher energy efficiencies. All the above-mentioned features can be controlled by a wide range of adjustable parameters of pulsed power supply such as pulse repetition frequency, pulse shape, pulse rise time, pulse fall time, pulse width, and pulse voltage amplitude [16, 21, 32,33,34].

For the atmospheric pressure nanosecond pulsed methane plasma, the in-situ and real-time diagnosis methods for particle densities are still rare, so numerical simulation becomes an alternative way to study these types of plasmas. Although simulation studying spatial-temporal evolution characteristics of particle densities and reaction pathways of producing and losing particles in the non-equilibrium AC methane plasma have been attempted a lot, numerical studies of the conversion in the nanosecond pulsed discharge plasma have been less reported so far. The self-consistent 2-D or 3-D models on the nanosecond time scale should be developed to further research the reaction pathway contributions and mechanisms of producing and losing reaction species. Therefore, in this paper, the application of ns-pulses in a DBD reactor with one electrode containing charge injection points is discussed in terms of methane conversion, selectivity towards the key products and the formation rate of carbon deposits. It is discussed how the pulse shape and duration impacts the properties of discharges and affects the rate of species production. To this end, a multi-species fluid model that combines momentum, mass transport equations for each species as well as energy balance relations with the Poisson equation to anticipate the behaviour of the system will be used. In the next section, the structure of the reactor and the employed model will be described while section III presents and discusses the results of the simulation. Finally, section IV is devoted to the conclusion and summarizing the main findings.

Description of the Model and Structure

The structure of the plasma reactor considered in this work is shown in Fig. 1. It consists of two electrodes (grounded and powered) and a dielectric barrier layer with a dielectric constant of 3.0 that covers the upper electrode (grounded electrode). The bottom electrode (powered) is engraved to a periodic cone-shape sharp array, establishing charge injection phenomena [35]. The dielectric layer has a thickness of 0.5 mm and the distance between the tip of cones and surface of the dielectric layer is 0.5 mm. The radius of the circular cone base is 0.2 mm and the height of cones and separation distance between them are fixed to 0.5 mm and 1 mm, respectively. Also, since in reality, there is not an infinite sharp point the curvature radius of the cone tip was selected to be 25 μm. In a conventional DBD system, the gas breakdown voltage is governed by Paschens law which accounts for Townsend processes (electron impact ionization in the bulk of the plasma and secondary electron emission from the electrodes. The flux of the secondary electron emission is proportional to the velocity of ions, which is directly dependent on the electric field intensity. It is obvious that around sharp points where there is an accumulation of the electric field, the secondary emission process can be more intensive. These sharp points are defined as charge injection parts. In addition, the high reduced electric field from the surface irregularities creates a reactive atmosphere because electrons accelerate and collide energetically with other species, creating more reactive species. A high reduced electric field in a DBD generates high densities of charged particles at low values of specific deposited energy, which leads to higher energy efficiencies. However, the fabrication of electrodes with very sharp points is difficult and needs a higher financial cost. Therefore, in the design of the current structure, we tried to choose the smallest charge points that do not need more financial cost. The Knudsen number (Kn) for our structure at atmospheric pressure range is still sufficiently small, so a 3-D fluid approach would be able to replicate the real system and provide a correct thorough understanding of the discharge behaviour. The 3-D modelling of the whole system is very complex and requires a high run-time cost. However, with bearing in mind the periodicity of system, a 2D axisymmetric fluid modelling of a periodic unit will also be able to mimic the actual problem. This approach has widely been used in literature in simulation studies related to DBD devices used in gas conversion applications [36, 37]. The fluid approach used here is based on solving a set of coupled differential equations that express the conservation of mass, momentum, and energy, for the different plasma species. The mass conservation equation is given as:

$$\begin{aligned} {{\partial {n_j}} \over {\partial t}} + \overrightarrow{\nabla }.{\overrightarrow{\Gamma }_j} = {\omega _j}, \end{aligned}$$
(1)

where \(n_j\) is density of species \(j=i,e,n\) (ion, electron and neutral, respectively), \({\overrightarrow{\Gamma }_j}\) stands for flux and \({\omega _j}\) refers to species source term. The fluxes are calculated by the drift-diffusion approximation as:

$$\begin{aligned} \overrightarrow{{\Gamma _j}} = - {D_j}\overrightarrow{\nabla }{n_j}\pm {\mu _j}\overrightarrow{E} {n_j}, \end{aligned}$$
(2)

where the first term on the right-hand side describes the diffusion while the second term determines the electric drift which is zero for neutrals. Here, \(D_j\) and \(\mu _j\) are the diffusion and mobility coefficients of the species j, respectively, and \(\overrightarrow{E}\) is the electric field with the signs plus and minus for ions and electrons, respectively.

Fig. 1
figure 1

a Simplified structure of the DBD reactor, b 3D unit cell of the reactor

Full size image

To obtain the ion mobility coefficient \(\mu _{j}\) in \({\text{m}}^{2}/ {\text{Vs}}\), we first determine the binary mobility coefficients \(\mu _{ij}\) in the background gas i as follows:

$$\begin{aligned} \mu _{i,j}=0.514\frac{T_{\rm {gas}}}{P_{\rm {tot}} \sqrt{m_{ij} \alpha _{i}}}, \end{aligned}$$
(3)

where \(T_{\text {gas}}\) is the gas temperature in Kelvin, \(P_{\text {tot}}\) is the total gas pressure in Pascal, \(m_{ij}\) is the reduced molecular mass in amu units

$$\begin{aligned} m_{ij}=\frac{m_{i}m_{j}}{m_{i}+m_{j}}, \end{aligned}$$
(4)

and \(\alpha _{i}\) in Å\(^{3}\) represents the polarizability of the background gas i. The overall ion mobility coefficient \(\mu _{j}\) can then be calculated by Blanch\(^{,}\)s law,

$$\begin{aligned} \frac{\rm {p}_{tot}}{\mu _{j}}=\sum _{i=\rm {background}} \frac{P_{i}}{\mu _{i,j}}, \end{aligned}$$
(5)

where the summation is over each background gas i. The partial pressure \(\hbox {P}_{i}\) of species i is determined via concentration \(n_{i}\) and the equation of state. Finally, the ion diffusion coefficient can be derived from the Einstein relation

$$\begin{aligned} D_{j}=\frac{k_{\rm {B}}T_{\rm {ion}}}{e}\mu _{j}, \end{aligned}$$
(6)

where e is the electron charge, \(k_{\text {B}}\) is the Boltzmann constant and \(T_{\text {ion}}\) represents the ion temperature. The mobility coefficient, \(\mu _{j}\) is zero for the neutrals, but to determine diffusion coefficient \(D_{j}\) in \({\text{m}}^{2}{\text{s}}^{-1}\) of a neutral species j, first the binary diffusion coefficient \(D_{ij}\) is calculated as

$$\begin{aligned} D_{ij}=\frac{3}{16}\frac{k_{B}T_{\rm {gas}}}{p_{\rm {tot}}}\frac{(2\pi k_{\rm {B}}T_{\rm {gas}} /m_{ij})^{1/2}}{\pi \sigma ^{2}_{ij}\Omega _{D}(\Psi )}, \end{aligned}$$
(7)

where \(\sigma _{ij}\) is the binary collision diameter in Å,

$$\begin{aligned} \sigma _{ij}=\frac{\sigma _{i}+\sigma _{j}}{2}, \end{aligned}$$
(8)

and \(\Omega _{D}(\Psi )\) is the diffusion collision integral given by

$$\begin{aligned} \Omega _{D}=\frac{A}{\Psi ^{B}}+\frac{C}{e^{D\Psi }}+\frac{E}{e^{F\Psi }}+\frac{G}{e^{H\Psi }}, \end{aligned}$$
(9)

with \(\Psi =T_{\text {gas}}/\epsilon _{ij}\), \(\epsilon _{ij}=(\epsilon _{i}\times \epsilon _{j})^{0.5}\), and constants \(A=1.06\), \(B=0.16\), \(C=0.19\), \(D=0.48\), \(E=1.04\), \(F=1.53\), \(G=1.76\), and \(H=3.89\), where \(\epsilon _{j}\) and \(\sigma _{i}\) are Lennard-Jones parameters [38]. The overall neutral diffusion \(D_{j}\) can then be obtained from Blanch\(^{,}\)s law. Finally the electron transport coefficient are interpolated from look-up tables created by Boltzmann equation solver.

The source term \({\omega _j}\) in equation (1) is obtained by considering the volume reactions in which species are created or lost as bellow:

$$\begin{aligned} \omega _{j}=\sum \limits _k {[(a_{jk}^R - a_{jk}^L)} {k_k}\prod \limits _l {n_l^L]}, \end{aligned}$$
(10)

where \(a_{jk}^R\) and \(a_{jk}^L\) are the right-hand side and left-hand side stoichiometric coefficients of species j in reaction k , \(k_k\) is the reaction rate coefficient and \(n_l^L\) is the density of the l-th species in the left-hand side of reaction k.

In atmospheric and high-pressure non-thermal plasma with low ionization degree ions which have comparable mass with neutrals are usually in thermal balance with the background gas. Since in this work we assume a constant pressure and constant temperature for background gas the ions which are in thermal balance with neutral have also constant temperature. Therefore, the energy balance equations for ions and neutrals are disregarded in this model and it is assumed that ions have the same temperature of background gas and their temperature is constant throughout the reactors [36, 37]. The energy balance equation will be solved only for the electrons, which is expressed as:

$$\begin{aligned} {{\partial {n_\varepsilon }} \over {\partial t}} + \overrightarrow{\nabla }.{\overrightarrow{\Gamma }}_{\varepsilon }+ \overrightarrow{\Gamma }_e.\overrightarrow{E}= {S_\varepsilon }, \end{aligned}$$
(11)

where \(n_\varepsilon\) is the electron energy density, \(S_\varepsilon\) is the energy loss/gain due to collisions, and the flux vector for electron energy \({\mathbf {\Gamma }_\varepsilon }\) is given as

$$\begin{aligned} \overrightarrow{{\Gamma _\varepsilon }} = -{\textstyle {5 \over 3}}{\mu _e}\overrightarrow{E} {n_\varepsilon } - {\textstyle {5 \over 3}}{D_e}\overrightarrow{\nabla }{n_\varepsilon }, \end{aligned}$$
(12)

The set of fluid equations is completed with Poisson equation for electric potential:

$$\nabla .\left( {\varepsilon _{r} \varepsilon _{0} {\mathbf{E}}} \right) = - \nabla .\left( {\varepsilon _{r} \varepsilon _{0} \nabla V} \right) = e\left( {\sum\limits_{i} {n_{i} } - n_{e} } \right),{\text{ }}$$
(13)

\(\varepsilon _0\) is the permittivity of vacuum, \(\varepsilon _r\) the relative permittivity of the material or the gas, V is the potential and e is the elementary charge. Notice that, unlike in the discharge gap region, the transport and energy equations are not considered in the dielectric area and only the Laplace equation for electric potential is solved inside it. The plasma species considered in this work are listed in Table 1. Here, twenty species including electron, ions, radicals, and neutrals are assumed to participate in the plasma chemistry. Collision processes and reactions occurring between these species along with their reaction rate and threshold energies are presented in Table 2. A global model that considers more reactions and species has been developed and its results was used to select these main reactions and species, they play main role in underlying reaction pathways of DBD methane plasma. The global model used here is an atmospheric-pressure version of the model presented in reference [39], which deals with low-pressure plasma. In this model, the balance equations for heavy species and electron energy equation as well as gas temperature equation are solved with considering species loss on reactor surface. The electron density is calculated by plasma neutrality condition. Ions loss on surfaces is taken into account by thermal flux with different secondary emission coefficients on the dielectric layer and electrodes, while surface loss of neutral species is determined by sticking coefficients calculated by surface reactions. A consistent set of 367 gas-phase reactions involving the 36 defined species was built to describe the methane plasma chemistry. This model is applied for the same volume and same surface areas of the reactor presented in Fig. 1. A constant power between 1 W and 100 W is fixed to supply the reactor. After each run, species with higher densities are selected and the main reaction producing and consuming these species are determined. This process is repeated for different input powers and then results are compared with each other. Finally, we chose 20 species with 30 main reactions that produce and consume them. The reaction mechanism must be closed, which means that for each species there is at least a reaction consuming that species and a reaction producing it. During the numerical implementation, the rate constant of the interaction of the electron with background gases and neutrals is extrapolated from look-up tables constructed by the Boltzmann equation solver, in which the rate and transport coefficients as functions of electron mean energy are arranged. The rate constants of ion-neutral and neutral-neutral reactions have been exploited from reference [40]. The surface charge accumulation on the dielectric layer due to the difference in fluxes between the electrons and the ions is taken into account using the following boundary conditions:

$$\begin{aligned} \overrightarrow{n} .({\overrightarrow{D} _1} - {\overrightarrow{D} _2}) = \sigma , \end{aligned}$$
(14)

where \(\mathbf {n}\) is the unit normal, \(\overrightarrow{D} _1\) and \(\overrightarrow{D} _2\) are the electric displacement fields on both sides of the boundary and surface charge density, \(\sigma\), is calculated by

$$\begin{aligned} {{\partial \sigma } \over {\partial t}} = \overrightarrow{n} .{\overrightarrow{J} _i} - \overrightarrow{n} .{\overrightarrow{J} _e}. \end{aligned}$$
(15)

where, \(\overrightarrow{J} _i\) and \(\overrightarrow{J} _e\) are the total ion and electron current densities at dielectric surface, respectively. The upper electrode was grounded while the bottom electrode is powered by a voltage waveform combined of a repetitive negative nanosecond pulses and a positive dc bias. In our previous work we showed how this combination of negative pulse and positive dc can lead to higher efficiency and can be used to control the process in plasma gas conversion applications [35]. The pulses is designed with 3 index labels, corresponding to pulse rise time (\(t_{r}\)), pulse peak width (\(t_{p}\)) and pulse fall time (\(t_{f}\)) in nanoseconds, respectively. Pulses with different rise times, different widths, and different fall times are shown in Fig. 2. Where each pulse parameter changes independently and any change in a parameter does not affect other pulse parameters. The pulse shape is governed by the following relation

$$\begin{aligned} {\rm {V(t)}} = \left\{ \begin{array}{cc}{{V_0}} &{} {t \le T/2} \\ {{V_0}+{V_p}({\rm {t}} - {\rm {T/2}})/{t_r}} &{} {T/2 \le t \le T/2 + {t_r}} \\ {{V_0}-{V_p}} &{} {T/2 + {t_r} \le t \le T/2 + {t_r} + {t_p}} \\ {V_0}+ {V_p}(1.0 - {{t - T/2 - {t_p} - {t_r}} \over {{t_f}}}) &{} {T/2 + {t_r} + {t_p} \le t \le T/2 + {t_r} + {t_p} + {t_f}} \\ {{V_0}} &{} {T/2 + {t_r} + {t_p} + {t_f} \le t} \\ \end{array} \right. \end{aligned}$$
(16)

with positive dc bias \({V_0} = 400 V\) and negative pulse amplitude \({V_p} = -3300 V\). Where, as an example, single pulse 5.0–2.0–7.0 corresponds to \(t_r\)=5.0 ns, \(t_p\)=2.0 ns, and \(t_f\)=7.0 ns. The repetition frequency of pulsed voltage is assumed to be constant at 10 kHz [41,42,43,44].

Fig. 2
figure 2

Shape of nanosecond pulse of − 3300 V amplitude with repetition rate of 10 kHz and different rise times (a), different peak widths (b), and different fall times (c) superimposed on a dc bias of 400 V

Full size image

In the lateral side of the computational unit cell, the periodic condition must be applied. To this end, the fluxes of electron and its energy are set to zero. Also, zero charge boundary condition are assumed for Poisson equation i.e defining the normal electric displacement equal to zero. The combination of these boundary conditions implies periodicity, with a normal zero gradient of charged species across the boundary. The other boundary conditions, considered surface reactions and the parameters of secondary electron emission are as same as those of reference [40]. The field emission effects are not considered because these effects appear at size lengths smaller than \(10\,\mu \hbox {m}\) [45].

Table 1 Species considered in the reaction mechanism
Full size table
Table 2 Reactions Mechanism
Full size table

Simulation Results and Discussion

All of the coupled equations with their boundary conditions are solved in the COMSOL software, which is based on finite element method. The pressure is fixed to 1 atm and the background gas temperature is in room temperature. The gas temperature is assumed to be uniform throughout the reactor and constant during whole process. In DBD reactors having discharge gap less than 1 mm, gas temperature is nearly constant throughout the gap. To more generalized investigation, the background gas heating and heat transfer models can be connected with current approach, but it makes the problem more challenging and complex and need a higher run-time. Therefore, here we assume that temperature is 300 K, which in the experiment can be established by cooling systems connected to the reactor. A non-uniform mesh was applied to the computational domain as it was dense near the electrodes and the dielectric layers surfaces, specially near the cone tip, while it becomes coarser in other regions like the gap center and near the periodic boundary layer. The size of mesh elements near the cone tip was smaller than the cone tip curvature in order to capture charge injection phenomena and the effects associated to high electric field intensity. The mesh size near the walls was also smaller than the electron Debye length calculated for reaction conditions. Since the time scales for pulse on-time (in the ns range) and off-time (100 \(\mu\)s) were different by more than 3 orders of magnitude, the time steps for these intervals were not same. To capture any quick change of voltage in pulse-on interval, the time step was made much smaller. Also, the time steps were smaller than inverse of frequency of momentum transfer collisions. As mentioned before, the net voltage applied to the reactor is a composition of a positive DC bias voltage and a negative pulse voltage. All Figs. 3, 4, 5, 6, 7, 8, 9, 10, 11 and 12 are presented during pulse on-time, when both pulse voltage and dc bias voltage are simultaneously applied to system. Also, t = 0 in this figures is correspond to \(t=100 \mu s\) in Fig. 2. It means that initial values in these figure are the steady state values of the parameters in pulse off-time period (when only positive bias voltage is applied to system). In other word, before applying a single pulse, the simulation was carried out for the positive dc voltage. Once the plasma reached steady state the simulation of pulse on-time interval started with initial conditions calculated form pulse off-time interval simulation. The (instant) selectivity to the main products was calculated, which are spatiality averaged over the discharge gap domain. The instant selectivity for hydrogen and hydrocarbons are given by Eqs. (16 and 17), respectively, following the approach introduced in [40].

$$\begin{aligned} {S_{{H_2}}} = {{{n_{{H_2}}}} \over {2 \times {n_{C{H_4},\,converted}}}} \times 100\,\% \end{aligned}$$
(17)
$$\begin{aligned} {S_{{C_x}{H_y}}} = {{x \times {n_{{C_x}{H_y}}}} \over {{n_{C{H_4},\,converted}}}} \times 100\,\% \end{aligned}$$
(18)

where density is updated in each time step. The \(\text{CH}_4\) conversion is not reported here as it was very small for a single period of voltage profile, which this work is focused on. A significant \(\text{CH}_4\) conversion can be obtained for a total duration of a few seconds. The latter requires a very high computational time for the 2D model and it was beyond the scope of the present paper. It is noteworthy to mention that in a plasma reactor although the final amount of conversion is dependent on reactor structure and input power, time to reach a steady-state conversion depends mainly on the average time which the neutrals spend in plasma volume before leaving the reactor. In a steady state, rates of reactions that produce and consume a species are in balance with each other and the density of species does not change after that. Time-average densities of plasma species (electrons and ions) and radicals in all voltage periods are nearly the same, because input power for each voltage period is constant, while the density of product molecules and background gas change with time. The background gas initially is pure methane but with the conversion process the background gas change to a mixture of methane with other neutral product gases. In steady state, background mixture gas does not change and the mole fraction of each neutral species is constant. The time to reach this steady state depends on the rate of reactions between different species. The rate of conversion initially is fast because the density of methane gas is higher but it decreases by the reduction in density of methene and become zero in the steady state. If we consider reaction 22 and assume that conversion of \(50\,\%\) is achieved in steady state, the rate of this reaction gives:

$$\begin{aligned} {{d{n_{C{H_4}}}} \over {dt}} = - {k_{22}}{n_{C{H_4}}}{n_{CH}}. \end{aligned}$$
(19)

With taking account that average density of radical CH is \(1.0\times 10^{+15} m^{-3}\le n_{CH}\le 1.0\times 10^{+16} m^{-3}\) and \(k_{22}= 9.74\times 10^{-17}\) \(m^3s^{-1}\) we can integrate the equation from initial density \({{n_{0,\text{C}{\text{H}_4}}}}\) to \(50\,\%\) conversion of \(\text{CH}_4\), (\(={0.5{n_{0,\text{C}{\text{H}_4}}}}\)) as below:

$$\begin{aligned} \int _{{n_{0,C{H_4}}}}^{0.5{n_{0,C{H_4}}}} {{{d{n_{C{H_4}}}} \over {{n_{C{H_4}}}}}} = - \int _0^{{t_f}} {{k_{22}}{n_H}dt} \end{aligned}$$
(20)

This gives the value of time reaching the steady state as \(t_f\simeq 0.64\,s\,\,to\,\,6.4\,s\)

Fig. 3
figure 3

The spatial distributions of different time-averaged plasma parameters: the electron density, the ion density, the electron temperature, and the electric field during a single 5.5–2.0–5.5 pulse

Full size image

Figure 3 illustrates the spatial distributions of plasma parameters averaged over a single pulse. The maximum of electric field and electron energy distributions appear around the cone tip (sharp point), as expected. The mean electron density shows a broader peak than that of total ion density. With comparing the distribution of electron and ion, it is found out that ion distribution peak occurs close to cone tip while for the electron it appears far away from the electrode tip, nearly in the centre of discharge gap. In our previous work on argon plasma for this structure [35], we showed that, during each repetition period, both electron and ion density distribution peaks perform a reciprocating motion between the cone tip and the dielectric surface. The same type of reciprocating motion for peak densities of electron and ion has also been seen here for \(\text{CH}_4\) plasma. The electron and ion avalanches start to move away from the cone tip in the beginning of each pulse. Once they reach the dielectric layer surface during pulse width time \(t_{p}\) they are pushed back to cone tip region over pulse fall time and pule off-time. Periods of reciprocating motions for electron and ions are different. The later causes that the maximum of time-averaged densities of electron and ion do not occur near the cone tip or near dielectric surface and they occur in different positions.

Fig. 4
figure 4

a–c Main plasma parameters as a function of time in a single 5.5–2.0–5.5 pulse: a spatially averaged electron and ion densities, b overall production rate of different products, c selectivity to carbon containing species and \(\text{H}_2\) selectivity

Full size image

The time evolution of spatially averaged electron and ion densities, the production rate of different products, and selectivities are illustrated in Fig. 4a–c, respectively. It can be seen form Fig. 4a that before the pulse fall time, density figures of both electron and total ions are similar, experiencing a sharp rise. Then, during the fall time, the behaviours of electron and total ions are different. The density of electrons is dramatically decreased while for the ions the density profile still continues increasing in value. The difference in behaviour of electron and ions densities is related to mobility of these species. The electrons due to having higher mobility in comparison with the ions are able to response more quickly to the alternating electric field. Therefore, the electrons respond quickly as compared with the ions, while it takes more time till the ion start to respond to the new value of electric field. According to Fig. 4b, among six dominant neutral molecules produced in the reaction, hydrogen molecule has highest production rate at all times except in its local minimum which occurs when net applied voltage is zero in pulse fall time interval (\(t\simeq 12\,\hbox {ns}\)). Over pulse rise and width times the production rate of all species increases, expect for \(\text{C}_3\text{H}_8\) that experiences a decrease, more quickly in pulse width time. The behaviour of production rate of different species is not same during pulse fall time. From the start of pulse fall time \(t=7.5\,\hbox {ns}\), until \(t\simeq 12\,\hbox {ns}\), when net applied voltage is zero, the production the rate of \(\text{H}_2\) and \(\text{C}_2\text{H}_4\) decreases, while the rate of production of \(\text{C}_2\text{H}_6\) and \(\text{C}_3\text{H}_8\) increases. The rates of \(\text{C}_2\text{H}_2\) and \(\text{C}_3\text{H}_6\) production remain constant during in this time interval. At the critical time-point \(t\simeq 12\,\hbox {ns}\), the net applied voltage makes a transition from negative to positive value, so all plasma species experience a change in their production rate, although for some species this change remains relatively small. Figure 4c compares the selectivity towards the main reaction products after a single pulse. It can be seen that the selectivity of \(\text{C}_2\text{H}_4\), \(\text{C}_2\text{H}_6\), and \(\text{C}_3\text{H}_8\) decreases while selectivities of other products experience an increase in value. The next sections focus more on these.

Change in Pulse Rise Time

Figure 5 shows the effect of pulse rise times in the 5–7 ns range on the production rate of \(\text{H}_2\), ethylene and ethane. A significant increase in production rates can be observed at larger times, as the production shifts to later times with an increase of pulse rise time. The maximum production rate shows a growth, which is the highest for ethane and smallest for \(\text{H}_2\). Figure 5 also shows that an increase of pulse rise time from 5 to 7 ns reduces a sharp change in \(\text{H}_2\) concentration at critical time corresponding to a net applied voltage of zero. It also influences the production rates during pulse width time, pulse fall time and even pulse off-time.

Fig. 5
figure 5

Spatially averaged production rates of \(\text{H}_2\), \(\text{C}_2\text{H}_4\) and \(\text{C}_2\text{H}_6\) as functions of time for three pulses with different rise times: 5.0–2.0–5.0 ns (red), 6.0–2.0–5.0 ns (green), 7.0–2.0–5.0 (blue) (Color figure online)

Full size image
Fig. 6
figure 6

Selectivity to the main hydrocarbon species \(\text{C}_2\text{H}_4\) and \(\text{C}_2\text{H}_6\), and selectivity to \(\text{H}_2\) as functions of time for same pulses in Fig. 5: 5.0–2.0–5.0 ns (red), 6.0–2.0–5.0 ns (green), 7.0–2.0–5.0 (blue) (Color figure online)

Full size image

Figure 6 compares the selectivity towards three species for three different input pulse voltages that are different in pulse rise time. It can be seen from the figure the \(\text{C}_2\text{H}_4\) selectivity is considerably higher than that for ethane. An increase in pulse rise time changes the selectivity mainly over pulse off-time. During pulse on-time, change in the selectivities is only seen in pulse fall interval. An increase of pulse rise time from 5 to 7 ns decreases the \(\text{C}_2\text{H}_4\) selectivity from 29 to 20 \(\%\) while the \(\text{C}_2\text{H}_6\) selectivity drops from 22 to 17 \(\%\). At the same time, more \(\text{H}_2\) is produced and the \(\text{H}_2\) selectivity increases from 4 to 7 \(\%\) of hydrogen. For higher pulse rise times the difference between selectivities of ethylene and ethane are less than that of lower rise times.

Change in Pulse Width

In Fig. 7, the effect of changing pulse peak width on production rate of \(\text{H}_2\), ethylene and ethane is investigated. The figure indicates that with increasing peak width form 2 to 4 ns, the product distribution at the end of pulse (t = 14 ns) as well as during pulse-off time (t = 30 ns) demonstrates significant changes. The \(\text{C}_2\text{H}_6\) production rate shows the highest increase, followed by \(\text{C}_2\text{H}_4\) rate and then by \(\text{H}_2\) rate. At a total width of 14 ns, more time is also required to restore plasma to conditions observed during the initial positive bias voltage interval.

Fig. 7
figure 7

Spatially averaged production rates of \(\text{H}_2\), \(\text{C}_2\text{H}_4\) and \(\text{C}_2\text{H}_6\) as functions of time for three pulses with different pulse peak widths: 5.0–2.0–5.0 ns (red), 5.0–3.0–5.0 ns (green), 5.0–4.0–5.0 (blue) (Color figure online)

Full size image
Fig. 8
figure 8

Selectivity to the main hydrocarbon species \(\text{C}_2\text{H}_4\) and \(\text{C}_2\text{H}_6\), and selectivity to \(\text{H}_2\) as functions of time for same pulses in Fig. 7: 5.0–2.0–5.0 ns (red), 5.0–3.0–5.0 ns (green), 5.0–4.0–5.0 (blue) (Color figure online)

Full size image

The effect of pulse peak-width change on the \(\text{H}_2\) selectivity as well as \(\text{C}_2\text{H}_4\) and \(\text{C}_2\text{H}_6\) selectivities is shown in Fig. 8. The \(\text{C}_2\text{H}_4\) selectivity (measured during 30 ns) decreases from 27 to \(21.5\,\%\) as the pulse peak width increases from 2 to 4 ns. At the same time, the \(\text{C}_2\text{H}_6\) selectivity drops less significantly from 20 to \(17\,\%\). However, near a time of 9 ns, an inverse situation was observed when the ethane selectivity was just slightly higher at wider pulses. The \(\text{H}_2\) selectivity increases up to 1.5 times with an increase of total pulse width from 12.0 to 14.0 ns. It is found out that with changing total pulse width the selectivity of \(\text{C}_2\text{H}_4\) changes more significantly than selectivity of \(\text{C}_2\text{H}_6\).

Change in Pulse Fall Time

Figure 9 compares the production rate of \(\text{H}_2\), ethylene and ethane as a function of pulse fall time in the range 5–7 ns. This corresponds a total pulse duration between 12 and 14 ns. The increase of pulse fall time results in growth of production rate of all species after 12 ns. Similar to the effect of very short rise time described above, a sharp peak in the \(\text{H}_2\) production rate was observed near the critical time of 11 ns at the shortest pulse fall time. Increasing the fall time reduced the depth of this peak and increased the overall \(\text{H}_2\) production rate.

Fig. 9
figure 9

Spatially averaged production rates of \(\text{H}_2\), \(\text{C}_2\text{H}_4\) and \(\text{C}_2\text{H}_6\) as functions of time for three pulses with different fall times: 5.0–2.0–5.0 ns (red), 5.0–2.0–6.0 ns (green), 5.0–2.0–7.0 (blue) (Color figure online)

Full size image
Fig. 10
figure 10

Selectivity to the main hydrocarbon species \(\text{C}_2\text{H}_4\) and \(\text{C}_2\text{H}_6\), and selectivity to \(\text{H}_2\) as functions of time for same pulses in Fig. 9: 5.0–2.0–5.0 ns (red), 5.0–2.0–6 ns (green), 5.0–2.0–7.0 (blue) (Color figure online)

Full size image

Figure 10 highlights the effects of pulse fall time increase on selectivity of three species. Similar to effects of increasing pulse rise time and pulse width, increase in pulse fall time decreases the selectivity of \(\text{C}_2\text{H}_4\) and \(\text{C}_2\text{H}_6\) while it increases the selectivity of hydrogen molecule. However, in comparison with them, changes of selectivity due to pulse fall time increase is low. Change in selectivity of neutrals is seen from start time of pulse fall time, which implys that the rate of voltage change is a important key factor to control the selectivity.

Discussion of the Results

To understand the change in selectivities and overall production rate of products by variation of pulse parameters, the contribution of each reaction channel on the production and consumption of each neutral product should be investigated in detail. The instant rate coefficient of reactions that electron is one of its reactants (R1–R17 in Table 2) is determined by the instant reduced electric field. Since the reduced electric field is dependent on input voltage the rate coefficient profile of these reactions is governed by input voltage profile and any change in voltage profile must directly change the rate of these reactions. The calculated rate coefficients of these reactions during pulse on-time for pulse 5.5–2–5.5 are illustrated in Fig. 11a. It is clear that the rate coefficients of electron reactions are governed by the input net voltage profile (see Fig. 2). There are two dip values in each rate profile that correspond to time points that net voltage is zero, a result of the combination of the positive dc and negative pulse voltages. This figure implies that any change in voltage profile, which can be established by varying pulse parameters, can directly change rates of electron reactions that have the main roles in the plasma reaction mechanism. These electron reactions, which we here call elementary or primary reactions, in addition to ionization processes include dissociation processes that generate active radicals along with electrons and ions. The radicals and ions produced by primary reactions interact with each other, in reactions called here as secondary reactions, to produce new species. The rates of these secondary reactions (R18–R30 in Table 2) as functions of time for pulse 5.5–2–5.5 are presented in Fig. 11b. This figure reveals that rates of secondary reactions, like primary reactions, are affected by voltage change during pulse on-time, although for secondary reactions the rate profiles are not similar with pulse shape. The secondary reactions are important in the reaction mechanism not only in pulse on-time but also during pulse off-time when electron reactions have less contribution in the conversion process. Electron density according to Fig. 4 is too small during pulse off-time, so electron reactions must have smaller participation in the reaction mechanism during that time interval. Generally, any change in input pulse voltage profile and its characteristics change initially the rate coefficient of primary reactions (electron reactions), and then it affects the rate of secondary reactions in which radical species are their reactants. The rates of the reactions are not the same. The reactions having a high rate have the main role in the whole process and they receive the effect of voltage change very quickly than reactions having a small rate. In other words, all reactions consuming and producing a species do not have the same sufficiency for increase and decrease in the species density. To justify this, the next figure compares the participation of the reactions in producing and consuming background gas and neutral products \(\text{H}_2\), \(\text{C}_2\text{H}_4\), and \(\text{C}_2\text{H}_6\). Although it is presented for pulse condition 5.5–2–5.5 the chemistry pathways and reactions discussed for it are applicable for any input pulse condition. The chemistry pathways and reactions do not change by pulse shape variation and only their rates are affected by it, as discussed above. Figure 12 is prepared by calculation of relative contributions of each reaction in consuming and producing of species. The dominant reactions consuming the feedstock gas (\(\text{CH}_4\)) are

$$\begin{aligned} \rm {R6:} \,\,\,\,\,\, e^{-}+CH_{4}\rightarrow CH_{3}+H+e{-}, \end{aligned}$$
(21)
$$\begin{aligned} \rm {R7:} \,\,\,\,\,\, e^{-}+CH_{4}\rightarrow CH_{2}+H_{2}+e_{-}, \end{aligned}$$
(22)
$$\begin{aligned} \rm {R22:} \,\,\,\,\,\, CH_{4}+CH\rightarrow C_{2}H_{4}+H, \end{aligned}$$
(23)

that time evolutions of their contribution in overall consumption are presented in Fig. 12a. Reaction channels R6 and R7 have the same behavior. As both plots show a peak near the end of pulse width time and a dip at the critical point time when the net applied voltage is zero. After that, they firstly experience a quick increase and then a slight increase that continues until they reach a new steady state. The plot for reaction channel R22 shows a reverse behavior as compared to those for R6 and R7. Initially, R22 has a high rate, while after \(t=12\,\hbox {ns}\), R6 and R7 become the main channels consuming background gas. The main reason for this is that the electron density before pulse on-time, according to Fig. 4a, is very small, so rates of reactions 6 and 7 are too small and they have less participation in background gas consumption. When the negative single pulse with higher voltage amplitude is on, background electrons receive more energy to ionize the background gas which leads not only to consuming of background gas but also causes the increase of electron density. An increase in electron density again increases the rate of reactions 6 and 7 during the pulse on-time. Lower electron density during pulse-off time is related to electron higher mobility that plays the main factor in the loss of electrons to surfaces, while radical species without electric charge are not affected by the electric field and have a longer lifetime than electrons, so they play the main role in the consumption of background gas during the pulse-off time. Hydrogen molecule, according to Fig. 12b, is consumed by reaction R15 (\(e^{-}+\text{H}_{2}\rightarrow 2\text{H}+e{-}\)) while two main reactions

$$\begin{aligned} \rm {R7:}\,\,\,\,\,\,e^{-}+CH_{4}\rightarrow CH_{2}+H_{2}+e_{-}, \end{aligned}$$
(24)
$$\begin{aligned} \rm {R21:}\,\,\,\,\,\,CH_{2}+H\rightarrow CH+H_{2}, \end{aligned}$$
(25)

compete with each other in producing \(\text{H}_2\). Reaction R21 is dominant in producing \(\text{H}_2\) during pulse off-time (\(t=0\)) while during pulse on-time with the increase of electron density and electron energy the reaction R7 surpass reaction R21 and become a more important reaction. \(\text{C}_2\text{H}_4\) molecules is dominantly consumed by reaction R27: \(\text{C}_{2}\text{H}_{4}+\text{CH}_{4}+\text{H}\rightarrow \text{C}_{2}\text{H}_{5}+\text{CH}_{4}\) which its rate is nearly constant during pulse on-time (see Fig. 11b). Two bellow reactions produce \(\text{C}_2\text{H}_4\) molecules during the whole process:

$$\begin{aligned} \rm {R12:}\,\,\,\,e^{-}+C_{3}H_{8}\rightarrow C_{2}H_{4}+CH_{4}+e^{-} \end{aligned}$$
(26)
$$\begin{aligned} \rm {R22:}\,\,\,\,\,\,CH_{4}+CH\rightarrow C_{2}H_{4}+H \end{aligned}$$
(27)

The reaction R22 has the main contribution in the production, nearly \(90\%\) while reaction R12 participate in the production only during pulse on-time with the contribution of less than \(10\%\). Finally, Fig. 12a shows how contributions of reactions consuming and producing \(\text{C}_2\text{H}_6\) vary during pulse on-time. It reveals that the production of \(\text{C}_2\text{H}_6\) is only governed by reaction R23, which has a constant rate during all times. Consumption \(\text{C}_2\text{H}_6\) is established by reactions R10 and R20. The latter is dominant in the pulse off-time while reaction R10 becomes the dominant reaction by pulse voltage application. If the rate of a reaction presented in Fig. 12 changes the selectivity of species that is produced or consumed by that reaction will change. This implies that overall production rates and selectivities will change with pulse shape if reactions in Fig. 12 are dependent on voltage and their rate is governed by input pulse shape. The dependency of reactions to voltage and their change during pulse on-time was presented and discussed in Fig. 11 and it explained that rate of most of reaction reactions change by pulse shape.

Fig. 11
figure 11

Rate coefficients of electron reactions (R1–R17 in Table 2) (a) and rates of secondary reactions (R17–R30 in Table 2) (b) during pulse on-time for pulse 5.5–2.0–5.5

Full size image
Fig. 12
figure 12

a Contribution of each dominant channel reaction in consuming and producing the \(\text{C}_\text{H}4\) (a) \(\text{H}_2\) (b), \(\text{C}_2\text{H}_4\) (c), and \(\text{C}_2\text{H}_6\) (d) during pulse on-time for pulse 5.5–2.0–5.5

Full size image

Conclusions

The reactions pathway of methane plasma affected by input voltage pulse shape in a micro-DBD was investigated. To this end, a micro-DBD including an electrode with special shape, having sharp points aimed to establish charge injection effects, was assumed to be powered by a negatively pulse voltage superimposed on positive DC voltage. A multi component fluid model was used to simulate the system and to investigate the production rate of different neutral products, instant selectivity of neutral products and time evolution of plasma parameters for different shapes of input pulse voltage. The shape of input pulse voltage was changed by changing its characteristic times such as pulse rise time, pulse width, and pulse fall time.

Results of study showed that reaction rates in methane plasma formed inside the structure are strongly dependent on input pulse shape. Any small increase in one of characteristic times increased the production rate of neutral species not only during the pulse on-time but also its effects was seen over pulse off-time. The electrons respond very quickly to voltage change while heavy species like ion needed more time to update itself with new electric field. The production rates around the critical time point, when the net applied voltage becomes zero in pulse fall interval, were more affected by input pulse shape and their profiles changed remarkably with that. Time evolution of selectivities of produced hydrocarbons and hydrogen molecule were studied for different input pulses. The selectivity for hydrocarbons increased by increasing of pulse rise time and pulse fall time, more vigorously after pulse on-time. A important parameter that affected the selectivity of products was rate of voltage change, which changes by change in pulse rise and fall times. Generally, results explicitly proved that reaction pathways of methane plasma, even any other plasma gases, can be controlled by input pulse shape in order to increase the selectivity of a desired hydrocarbon and/or amount of a specific products.

In outlook for next work, we attend to include ferroelectric packing material in the structure and also add the electrohydrodynamics effect in our study. In this way, due to hysteresis loop of polarization of ferroelectric materials the effects of electric field of pulse can be transferred to pulse off-time, leading to increase of the plasma activity after pulse on-time. Also, due to high intensity of electric field around the sharp points, electrohydrodynamics flow can rotate the background gas to make plasma more uniform.