Document Type : Research Paper


1 Assistant Professor, Department of Chemistry, Payame Noor University, Tehran, Iran

2 professor ,Photonics Research Institute, Institute of Science and High Technology and Environmental Sciences, Graduate University of Advanced Technology, Kerman, Iran

3 Associate professor,Department of Mechanical Engineering, Graduate University of Advanced Technology, Kerman, Iran

4 2 Associate professor ,Photonics Research Institute, Institute of Science and High Technology and Environmental Sciences, Graduate University of Advanced Technology, Kerman, Iran 3 Associate professor,Department of Mechanical Engineering, Graduate University of Advanced Technology


In this paper, using a one-dimensional simulation model, the reforming process of sour gas, i.e. CH4, CO2, and H2S, to the various charged particles and syngas in a dielectric barrier discharge (DBD) reactor is studied. An electric field is applied across the reactor radius, and thus a non-thermal plasma discharge is formed within the reactor. Based on the space-time coupled finite element method, the governing equations are solved, and the temporal and spatial profiles of different formed charged species from sour gas inside the plasma reactor are verified. It is observed that the electric field increases radially towards the cathode electrode. Moreover, the electron density growth rate at the radial positions closer to the cathode surface is smaller than the one in the anode electrode region. Furthermore, as time progresses, the positive ions density near the anode electrode is higher. In addition, the produced syngas density is mainly concentrated in the proximity of anode dielectric electrode.


Natural gas or methane (CH4) with carbon dioxide (CO2), hydrogen sulfide (H2S), and other sulfides is called sour gas (acidic gas). Major natural gas pollutants such as H2S and CO2 cause the corrosion and damage to the natural gas processing equipment and pipelines. It has serious consequences on the health of the ecosystems and living organisms and major financial losses. H2S density in a gas reservoir can vary from zero to 98%. In fact, its toxicity level is comparable with the hydrogen cyanide (HCN) gas. The sour gas is undesirable as it is seriously harmful, and its poisonous sulfur components are highly corrosive (Marsland et al., 1989(.

The natural gas can be reformed by means of many methods. So far, depending on conditions and the industrial demands, many different technologies such as steam methane reforming, steam methane and oxygen reforming, partial oxidation, catalytic partial oxidation, auto-thermal reforming, carbon dioxide reforming, combined auto-thermal reforming, sulfur passivated reforming, and gas heated reforming are applied. However, the partial oxidation method has taken much more industrial attentions (Khoshnoodi, and Mohammad, 2005).

Abdel-Aal et al. studied the conversion of sour natural gas to syngas (H2 and CO) based on the non-catalytic partial oxidation of the sour natural gas (Abdel-Aal et al., 1999). Their method was integrated with a process with ammonia or methanol. Moreover, the optimization of the syngas production via non-catalytic auto-thermal partial oxidation of methane was studied by Khoshnoodi et al. (Khoshnoodi, and Mohammad, 2005). Their software was based on the minimization of the total Gibbs energy.

Over the last two decades, significant efforts have been made on the industrialization of plasma discharges. These efforts are mainly including the cleaning of exhaust gas and ores. Usually, the plasma discharge reactors are considered as gliding arc plasma and cylindrical dielectric barrier discharge (DBD). In comparison with other methods such as selective catalytic reduction method, they have demonstrated higher efficiency and fewer side problems (Tao et al., 2008).

Plasma technologies for gas reforming have taken significant attention recently as many problems with chemical techniques can now be eliminated. The main role of plasma discharge in sour gas reforming is to provide enough energy to create the required free radicals by the process (Lieberman and Lichtenberg, 1994; Becker and Loffhagen, 2009). Catalytic effects of non-thermal plasma on the production of hydrogen and the separation of hydrogen sulfide as well as the fracture of natural gas have been investigated by other researchers )Hu et al., 2005).

Czernichowski (2001) used plasma discharge as an activation medium for the partial oxidation of natural or associated gases along with the syngas production. They showed that plasma discharges play the role of catalysts. Moreover, they are very active and provide the necessary energies for such reactions.

The production of syngas based on the reforming of CH4-CO2 mixture using the cold and thermal plasma discharge processes has been studied by Tao et al. (2011). They particularly focused on attaining higher conversions at high feed-gas flow rate and reducing the energy consumption to meet industrial requirements of particular applications. Hence, the emphasis was on the electron density, plasma temperature, and reactor configuration. They found that the energy conversion efficiency and the treatment capacity of the process can be improved via the optimization of plasma discharge operating conditions and designing the DBD reactor.

Using a gliding arc plasma reactor, Thanompongchart et al. (2014) experimentally studied the production of the syngas through the biogas reforming process. They studied the effects of biogas composition (CH4/CO2), input power, biogas flow rate, conversion of CH4 and CO2, and energy consumption in the presence of air.

The plasma assisted conversion of pyrolysis gas (Pyrogas) fuel to syngas was studied by Odeyemi et al., 2012). The effectiveness of the plasma reactor based on gliding arc plasma assisted to remove the light and heavy hydrocarbons contained in Pyrogas was shown. Moreover, the conversion of hydrocarbons and carbon dioxide to syngas using both the steam reforming reactions and dry CO2 reactions at atmospheric pressure conditions was demonstrated (Odeyemi et al., 2012). 

So far, using various methods, many theoretical and experimental works have been performed on the reforming of sour gas at different operation conditions. Generally, the noble gases such as argon, helium etc. are added to the sour gas in DBD reactors to produce more ionized carriers inside the reactor. Mostly, the main concentration has been on the gas discharge mechanism and the structure of the reactor. Thus, a comprehensive study on the conversion of only the sour natural gas without assuming any additive and in a real industrial geometry is still to appear. Hence, in this paper, using a one-dimensional fluid computational model based on the space-time coupled finite element method (FEM), the temporal and spatial behavior of only the sour gas reforming and the formation of different species from the sour gas is studied. A special attention has been paid to the production of syngas from the sour gas in the DBD reactor.  

2. Simulation model

Generally, in a suitable and complete mathematical model to study the conversion of sour gas molecules in a DBD reactor, the following phenomena should be taken into accounted: (a) the spatial and temporal variations of every formed species and their corresponding reaction rates, and (b) the ionization process inside the reactor. Because of the significant role of these phenomena inside the DBD reactor, the complete simulation of different formed species inside the plasma reactor, encompassing these phenomena together, could be interesting.

As shown in Figure 1, the DBD reactor is considered to be cylindrical with the inner and outer electrodes and the density of background sour gas is assumed to be constant. The inner electrode is connected to earth (cathode) and the outer electrode (anode) is covered with quartz as a dielectric and connected to a high voltage direct current (HVDC) power source. The quartz dielectric is chosen as the dielectric material because of its high melting point and transparency. As depicted in Figure 1, this tube has a thickness of 0.001 m. The internal and external electrodes are both assumed to be steel. The diameters of the inner and outer electrode are 0.04 m and 0.05 m respectively.






Figure 1

A schematic representation of the DBD plasma reactor.

As gas moves along the reactor axis, the application of electric potential between the two electrodes leads to the formation of a radial electric field. Subsequent electron acceleration leads to the ionization of the background sour gas. During the ionization process, free radicals, ions, and neutral atoms are created inside the DBD reactor.

Commonly, the spatial and temporal macroscopic description of the gas discharge inside the reactor is via solving the fluid continuity equations for different species together with Poisson’s equation to derive the electric potential. Here, without loss of generality, it is assumed that the axial and azimuthal behavior of the formed plasma discharge inside the plasma reactor is uniform. Hence, the spatial description of the problem is mathematically one dimensional (radial direction only).

The continuity equations for all the formed species inside the DBD reactor is expressed as follows:


where, ni is the number density, and Γi expresses the flux for species i. Ri,m is the reaction rate between species i and species m. Generally, the reaction between species A and B producing species C and D inside the DBD reactor is written as follows:


Here, reaction rates are obtained through those relations that are based on the density of each species (Lieberman and Lichtenberg, 1994; Becker and Loffhagen, 2009). In fact, each every reaction rate is directly proportional to the product of the density of each reactant raised to the power of that reactant coefficient from Equation 2, i.e.:


where, k is the performance reaction constant, which is assumed to remain unchanged in the reactor environment (Lieberman and Lichtenberg, 1994; Becker and Loffhagen, 2009). In this work, to find the reaction constant, two different approaches are considered. Firstly, for some reactions, the experimental data for their reaction rates are simply available in the literature (Westley, 1982; Dorai, 2000). Secondly, for reactions with no available experimental data, the reaction rate constants are calculated using the total collision cross sections in terms of the electrons energy. To this end, the following relationship between collision cross sections and the constant reaction rate is used (Janev and Murakami, 2001; Prigogine, 2003):


where, IP is a parameter close (but not always equal) to the ionization or appearance potential for a given ionization channel (expressed in electron-Volts (eV)). E is the collision energy, and Bi(i=1…N) are fitting coefficients; and N is determined from those conditions for achieving an standard deviation of the fit from the data smaller than 3-5% (Janev and Murakami, 2001). The performance reaction rate constant (k) and the collision energy (E) are written as given below (Prigogine, 2003):


Furthermore, it must be noted that in the calculation of the constant reaction rate, the most important reactions occurring in the DBD reactor are taken into the account. Thus, as presented in Table 1, the constant rates for all the reactions inside the DBD reactor are calculated at room temperature (about 300 K) and at atmospheric pressure.

Table 1

All reaction rates for different species formed in the non-thermal DBD reactor.

Reaction rate constants ( )







The flux term in the continuity fluid equation for all species (Equation 1) is based on the momentum conservation of each species. The corresponding flux term for each species depends on the electrical mobility and spatial diffusivity of that species, i.e.:


where, μi and Di are the mobility and diffusion coefficients for species i respectively. The plus or minus sign in Equation 7 accounts for the sign of the charged particle (Lieberman and Lichtenberg, 1994; Becker and Loffhagen, 2009).

Moreover, to derive electric potential inside the DBD reactor, the Poisson’s equation must be solved. Thus, its general form is defined by (Lieberman and Lichtenberg, 1994):


where, E is the electric field, and j is the electric potential. Moreover, writing down the charge density ρ in terms of the number density of the charged species; the Poisson’s equation can be rewritten as follows:


where, e0 is the vacuum permittivity coefficient; ni and qi are the number density and electrical charge of the ith charged species respectively (Lieberman and Lichtenberg, 1994).

To solve the coupled Equations 1 and 8, a stabilized FEM has been employed. To achieve a stabilized spatial discretization of the continuity equations for the charged particles, an upwind Petrov–Galerkin FEM is used. Moreover, based on the standard Galerkin FEM, a discretization process is performed on the continuity equation of the uncharged particles as well. In all the calculations, the time discretization is carried out through an implicit single step scheme; thus, there is no need of extreme small time steps.

In order to obtain a unique solution for the continuity fluid equations together with Poisson’s equation based on the geometry presented in Figure 1, necessary boundary conditions must be imposed (Dirichlet and Neumann boundary conditions). However, the applied boundary conditions for the DBD reactor are similar to those which can be found in the existing literature (Hagelaar et al., 2000). Generally, the following boundary condition is used to account for the particles flux which may be included as the secondary-electron emission:


where, n is the normal vector pointing toward the reactor wall, and vth is the thermal velocity inside the reactor, which can be written as: vth= (8kBTm)1/2. Moreover, the number a is defined as follows (Hagelaar et al., 2000):


For electrons, as a special case, the particles flux due to secondary emission is added to the system and is defined as follows:


where, subscript e indicates the electron charge and the summation term is over the ion species impinging the reactor wall. Moreover, the secondary emission coefficient is the average number of electrons emitted per impinging ion. Additionally, the boundary condition of 3000 VDC is applied to the cylindrical wall of the DBD reactor.

Moreover, only the mobility coefficients for ions and electrons are included (Raju, 2005). The mobility of ions is calculated according to the Langevin equation (Bleecker and Bogaerts, 2006):


where, α is the polarization of background gas per units of cubic angstroms and its value for various gases is presented in the existing literature on gaseous discharges (Bleecker and Bogaerts, 2006). In this work, they have been calculated as 0.0012, 0.001, 0,0014, and 0.0009 m2/Vs for CO2+, O2+, CH4+, and O- species respectively. The ion diffusion coefficient is calculated from Einstein relation. Moreover, the diffusion coefficients for neutral species are calculated using the distribution coefficients of Lennard-Jones (Bleecker and Bogaerts, 2006; Marrero and Mason, 1972; Jasper and Miller, 2014; Laurendeau, 2005) as follows:


In this work, the formation of 23 different species is taken into the account. Thus, 23 continuity equations together with the Poison’s equation are solved. Then, the continuity equations are transformed into the Galerkin weak formulation. To this end, firstly, the continuity equations are converted into the first order partial differential equations (Becker et al., 2009). Later, the obtained equations are multiplied by the test functions expressed in Table 2. In Appendix A, the derivation of the temporal and spatial integration of Equations 1 and 8 is presented. Finally, the boundary conditions are applied, and the obtained stiffness matrix for the whole problem is solved.

Table 2

Specific test functions corresponding to all the formed species inside the DBD reactor.

Index No.


Test function

Index No.


Test function









































































3. Simulation results

In this section, using the space-time coupled FEM based on the Galerkin scheme, Equations 1 and 8 are solved numerically, and the temporal and spatial behavior of different formed species inside the DBD reactor are described. Moreover, t0 and tf are the initial and final moments respectively and h represents the time step. The used simulation parameters in this work are presented in the Tables 3 and 4.

Generally, when an electron comes out from the cathode electrode via the cathode bombardment, it will be accelerated by the applied electric field and passes through the gaseous medium. Its movement causes the ionization and excitation of the background gas molecules. The production of newly born charged particles inside the DBD reactor results in the distortion of the applied electric field distribution inside the reactor and may affect its performance. The employed model describes the injection of electrons from cathode surface and their transport along the radial direction inside the reactor (Lieberman and Lichtenberg, 1994).

Table 3

Simulation data.



Mobility, μe, m2/Vs (for electrons)


Temperature, T, K (for electrons)


Potential amplitude, KV


Diffusion coefficient (for electrons)


Initial density, m-3 (for all species)


Initial density, m-3 (for electrons)


Initial density, m-3 (for H2S,CO2,CH4)


Total Time, μs


Time-steps, s


Table 4

Diffusion coefficients for the various formed species inside the DBD reactor.


Diffusion coefficient, D, m2/s











































Figure 2 shows the spatial and temporal distribution of the electric field across the DBD reactor. As seen, the electric field decreases in the proximity of the anode dielectric. This is due to the formation of negatively charged particles and their accumulation in the region closer to the anode dielectric which is known as ionization region. On the other hand, the positive ions slowly drift towards the cathode. Thus, the electric field will increase at the positions in the proximity of the cathode electrode, as seen in Figure 2. The observed behavior of the electric field across the reactor is in good agreement with the finding by Khodja et al. (2012). Additionally, via a numerical study, a similar result was reported by Kang et al. (2003). Moreover, our result corroborates the numerical findings by Becker et al. (2009) and the analysis of the microdischarges in an asymmetric DBD system with argon as the background gas (Becker et al., 2013).


Figure 2 

The spatial and temporal variation of the electric field inside the DBD reactor.

To complement the discussions on the electric field across the DBD reactor, the temporal and spatial variations of electric potential (E=gradV) are presented in Figure 3. It almost shows similar variations at different radial positions. Due to the increased ionization, a drastic temporal increase in the electric potential is observable and follows the temporal variations of the electric field. On the other hand, the electric potential varies almost linearly versus radius (Lieberman and Lichtenberg, 1994).


Figure 3 

The spatial and temporal variation of the electric potential inside the DBD reactor.

The considered continuity model for the DBD reactor helps us to understand the dynamical behavior of the charged particles inside the reactor. Generally, electrons play a significant role in the plasma growth phenomena inside the reactor. The temporal and spatial (radial) variations of electron density inside the plasma reactor are shown in Figure 4. As can be seen, after applying the electric field (Ea) along the radial direction, the electron density increases toward the reactor walls. However, owing to the space charge formation and their induced electric field (Es), the total electric field (Ea+Es) decreases. This results in the saturation of electron density as time progresses. Owing to the charge accumulation on the dielectric anode, the growth of electron density on the cathode surface is smaller than that of the anode electrode. This is in good agreement with the finding by Nikonov et al. (1999). Moreover, Braun et al. (1992) observed similar spatial behavior for electron density in a DBD reactor. In addition, our result corroborates the finding by Khodja et al. (2012) for DBD in a gas mixture of Ne–Xe.


Figure 4

The spatial and temporal demonstration of the number density of electrons inside the DBD reactor.

Contrary to the electrons, the electric field accelerates the produced positive ions (CO2+) in the direction of the applied electric field. If a CO2+ ion gets enough acceleration from the electric field, then it collides with the background gas molecules and ionizes them. However, this phenomenon rarely happens in the non-thermal plasmas and CO2+ ions are less effective than electrons in the impacting ionization processes (Becker et al., 2013). As shown in Figure 5(a), the density of CO2+ ions grows in the proximity of the anode dielectric. It is clear that an increase in the density of CO2+ ions is higher near the anode electrode than in the proximity of cathode dielectric electrode. Generally, in the DBD reactor, positive ions are distributed evenly over the inter-electrode radial distance. Figure 5(a) shows that as time progresses, the number of CO2+ ions near the anode electrode (ionization region) increases. Furthermore, the behavior of positive ions corroborates the findings by Braun et al. (2013). They used a 2D modeling for the micro-discharges in a device in which one electrode is covered with a dielectric. On the other hand, as seen in Figure 5(b), the number density of negative oxygen ions (O-) increases due to the electron attachment in the radial positions closer to the anode dielectric. However, they are relatively fewer than the formed positive ions inside the reactor. These findings are in agreement with the reported results by Fridman (2008) and Raizer et al. (1997). Furthermore, similar spatial behavior for the density of ions in a DBD system was observed in the computational work performed by Khodja et al.





Figure 5

The spatial and temporal variations of the number density of a) CO2+ ions and b) O- ions inside the DBD reactor.

The spatial and temporal variations of syngas (H2 and CO) production from the sour gas in the DBD reactor are presented in Figure 6. Similar to the positive and negative charged species, the density distribution of the molecular hydrogen and carbon monoxide are mainly concentrated in the ionization region (in the proximity of anode dielectric). However, in comparison with the other formed species inside the reactor, the production rate of syngas is quite slow. In comparison with the other conventional methods such as chemical schemes for the syngas production from the sour gas, the existing experimental reports with the DBD reactor show that the production rate for CO is significantly higher. The DBD reactor works at much lower temperatures (room temperature) and the energy costs are very low. Moreover, it does not cause corrosion in the gas reforming system, and there is no need for catalyzer. Additionally, contrary to the chemical methods, the production of sulfide dioxide is not reported experimentally (Abdel-Aal et al., 1999; Lim et al., 2009).




Figure 6

The spatial and temporal variations of the number density of a) H2 molecules and b) CO molecules.

4. Conclusions

In this paper, a one-dimensional computational fluid model was used to study the reforming of the sour gas to the different charged species and syngas production inside a DBD reactor. A decrease and an increase in the electric field were observed towards the anode dielectric and the cathode electrode respectively. It was seen that, owing to the greater ionization of the neutrals in the ionization region of the reactor, the effect of electric field on the different space charge densities at the radial positions closer to the anode electrode is higher. The number density of O- was seen to increase at the radial positions closer to the anode dielectric, and it is higher than the carbon dioxide positive ions density. Furthermore, the production of syngas mainly occurs at the radial positions near the anode dielectric. Finally, in this work, the input and output flow of sour gas was kept constant. However, if the temporal and spatial variations of the input and output gas velocity would be considered, a complete description could be achieved for the generation of all the formed species in the DBD reactor. Moreover, the performed simulation model was one dimensional; however, a three dimensional fluid simulation will certainly give more detailed results although a higher processing power is definitely necessary.


We would like to acknowledge the Institute of Science and High Technology and Environmental Sciences for financial support within the project No. 7/3/94/5147-2/12/94.



: Mobility for electrons, [m2/Vs]


: Temperature, [K]


: Number density for i species, [1/m3]


: Diffusion coefficient for i species, [m2/s]


: Test functions corresponding to the formed i species


: Electric Potential, [kV]


: Number density for electron, [1/m3]


: Reaction rate constants for i species, [m3/s]

Appendix A

The temporal and spatial integration of Equations 1 and 8 are derived as given below:























































Abdel-Aal, H, Shalabi, M., D. Al-Harbi, and Hakeem, T. ,Simulation of the Direct Production of Synthesis Gas from Sour Natural Gas by Noncatalytic Partial Oxidation (NCPO), Thermodynamics and Stoichiometry, Industrial & Engineering Chemistry Research, Vol. 38, No. 3, p. 1069-1074, 1999.
Becker, M., Loffhagen, D., and Schmidt, W., A Stabilized Finite Element Method for Modeling of Gas Discharges, Computer Physics Communications, Vol. 180, p. 1230-124, 2009.
Becker, M., Loffhagen, D., and Schmidt, W., A Stabilized Finite Element Method for Modeling of Gas Discharges, Computer Physics Communications, Vol. 180, p. 1230-1241, 2009.
Becker, M., Hoder, T., and Brandenburg, R., Analysis of Microdischarges in Asymmetric Dielectric Barrier Discharges in Argon, Journal of Physics D: Applied Physics, Vol. 46, No. 35, p. 355-203, 2013.
Bleecker, K. De, Bogaerts, A., and Goedheer, W., Detailed Modeling of Hydrocarbon Nanoparticle Nucleation in Acetylene Discharges, Physical Review E, Vol. 73, p.026405-16, 2006.
Braun, D., Gibalov, V., and Pietsch, G., Two-dimensional Modeling of the Dielectric Barrier Discharge in Air, Plasma Sources Science and Technology, Vol. 1, p. 166, 1992.
Czernichowski, A., Glid, Arc Assisted Preparation of the Synthesis Gas from Natural and Waste Hydrocarbons Gases, Oil & Gas Science and Technology, Vol. 56, No. 2, p. 181-198, 2001.
Dorai, R., Modeling of Plasma Remediation of NOx Using Global Kinetic Models Accounting for Hydrocarbons, University of Illinois at Urbana-Champaign, 2000.
Fridman, A., Plasma Chemistry: Cambridge University Press, 2008.
Hagelaar, G., De Hoog, F., and Kroesen, G., Boundary Conditions in Fluid Models of Gas Discharges, Physical Review E, Vol. 62, p. 1452, 2000.
Hu, Y., Li, G., Yang ,Y., Gao, X., and Lu, Z., Hydrogen Generation from Hydro-ethanol Reforming by DBD-plasma, International Journal of Hydrogen Energy, Vol. 37, p. 1044-1047, 2012.
Janev, R. K., Murakami, I., Kato, T., and Wang, J., Cross Sections and Rate Coefficients for Electron-Impact Ionization of Hydrocarbon Molecules, National Institute for Fusion Science, Toki, Gifu (Japan), Vol. 77, p. 161-213, 2001.
Jasper, A. W. and Miller, J. A., Lennard–Jones Parameters for Combustion and Chemical Kinetics Modeling from Full-dimensional Intermolecular Potentials, Combustion and Flame, Vol. 161, p. 101-110, 2014.
Kang ,W. S., Park, J. M., Kim, Y., and Hong, S. H. , Numerical Study on Influences of Barrier Arrangements on Dielectric Barrier Discharge Characteristics, Plasma Science, p. 504-510 IEEE,2003.
Khodja, K. and Belasri, A., Sheath Formation Study in Ne–Xe DBD Discharge, Radiation Effects and Defects in Solids, Vol. 167, p. 734-742, 2012.
Khoshnoodi, M. and Mohammad, M. D., Optimization of Syngas Production in Noncatalytic Auto Thermal Partial Oxidation of Methane Process, the First Iranian Conference on Combustion, 2005.
Laurendeau, N. M., Statistical Thermodynamics: Fundamentals and Applications, Cambridge University Press, 2005.
Lieberman, M. A. and Lichtenberg, A. J., Principles of Plasma Discharges and Materials Processing, MRS Bulletin, Vol. 30, p. 899-901, 1994.
Lim, M. S., Hong, M. S., and Chun, Y. N., Production of Synthesis Gas from Methane Using Compression Ignition Reformer, Korean journal of Chemical Engineering, Vol. 26, No. 4, p. 1022-1027, 2009.
Marrero, and Mason, E. A., Gaseous Diffusion Coefficients, Journal of Physical and Chemical Reference Data, Vol. 1, p. 3-118, 1972.