Numerical studies of the voltage amplitude effect on plasma characteristics in atmospheric pressure driven by dual LF-RF frequency discharge

1.   Introduction

Atmospheric pressure dielectric barrier discharges (DBDs) refer to the generation of plasma in an open atmospheric environment, characterized by lower macro gas temperatures, higher electron densities, and a wider variety of active particle species [1]. These DBDs possess favorable parameter characteristics and unique technological advantages, demonstrating significant potential for research and practical applications in enhancing material surface properties [2, 3], film deposition [4], vapor deposition, and environmental science and engineering [5]. Notably, their application in high-tech fields, particularly in biomedicine, has gained considerable attention in recent years [6–10]. Recently, the dual-frequency driven atmospheric pressure DBD technology was introduced, building upon traditional single-source driving methods by incorporating an additional driving source. This approach leverages the nonlinear synergy and electrical asymmetry between the dual-frequency sources, significantly altering electron power absorption during the discharge process and consequently modifying the electron distribution [11].

Dual-frequency excitation has been proposed to regulate plasma generation and the ion flux at the surface, particularly in terms of ion acceleration across the sheath [12, 13]. Numerous studies have explored the application of high frequencies in the MHz range, with recent investigations linking kHz and MHz frequencies to atmospheric pressure discharge plasmas [14, 15]. Liu et al demonstrated that superimposing radio frequency (RF) voltage onto low frequency (LF) improves the uniformity of the discharge [16]. Bazinette et al conducted diagnostic studies on dual-frequency driven DBDs using a 50 kHz/9 MHz combination [17]. Research indicates that dual-frequency application can effectively alter the morphology and energy transfer processes within the discharge plasma sheath, impacting the electron microscopic dynamics and further influencing plasma chemical dynamics, thereby enhancing the understanding of plasma characteristics [18]. Magnan et al investigated mode conversion in atmospheric pressure dielectric barrier LF-RF dual-frequency discharge, revealing that varying the amplitude of the dual-frequency voltage and gas composition allows for effective control of discharge modes, including α mode and α-γ mode conversions [19]. In our previous work, we explored various aspects of optimizing and controlling key parameters in DBD plasma. Specifically, we examined the effects of different external conditions such as frequency matching, harmonic number, fundamental frequency, and phase matching on discharge plasma parameters [20, 21]. From an energy transfer perspective, we analyzed how various physical effects influence the micro-dynamic behavior of charged particles during discharge, including the plasma sheath, electron collision ionization, and the spatio-temporal distribution of electron absorption power [22–25].

In the current study, we employ dual medium frequencies of 50 kHz and 5 MHz to investigate the amplitude effects of dual-frequency on electron dynamic behavior and the spatio-temporal distribution of particles. The paper is structured as follows. Section 2 introduces the simplified simulation model and relevant methodology. Section 3 presents and discusses the simulation results, and section 4 provides a conclusion.

2.   Model description

In present research, the fluid model was applied to simulate, which used the finite element calculation software COMSOL Multiphysics, the modeling and processing of the plasma discharge structure explored in this work include selecting boundary conditions, governing equations, setting initial conditions, establishing chemical reaction processes, and related constants, etc [24, 26].

The density of electrons, mean electron energy can be determined based on calculating the charged particles under the drift-diffusion approximation. The density continuity equation was described as follows:

The energy continuity equation is:

In the above formula, n_e refer to particle density of electrons, $n_{\mathit\varepsilon}$ means energy density, ${\boldsymbol{E}}$ represents the electric field strength, ${\boldsymbol{\varGamma}}_{\text{e}}$ denotes flux vector of electrons [27].

${\boldsymbol{\varGamma}}_{\text{e}}$ and ${\boldsymbol{\varGamma}}_{\varepsilon}$ can be expressed as:

${\mu _{\text{e}}}$ and ${\mu _{\varepsilon} }$ refer to mobility coefficients. ${\mu _{\varepsilon} }$ was expressed as:

${D_{\text{e}}}$ and ${D_{\varepsilon} }$ are the diffusion coefficients. ${D_{\text{e}}}$ and ${D_{\varepsilon} }$ were expressed as:

According to the formula, e is the elementary charge, $k_{\mathrm{B}}$ is the Boltzmann constant. R_e is the source term of the particle, S_en means the energy loss:

In the above formula, x_j is the molar fraction, N_n refers to neutral number density, and $\Delta {\varepsilon _j}$ is reaction energy loss. Assuming that a plasma chemical reaction system is composed of k = 1, …, Q particles and j = 1, …, N reactions, the heavy particle transport equation that all particles except Q satisfy is:

Here, $\rho$ is mixture density, $_{ }R_k$ is mass change fraction, $w_k$ is mass fraction, $\boldsymbol{u}$ represents speed. The definition of the diffusion flux vector j_k is as follows:

The research object in this work is the diffusion model of mixture averaging. For this model, the definition of ${\boldsymbol{V}}_k$ is as follows:

Here, $D_{k,\mathrm{m}}$ is average diffusion coefficient of the mixture, $_{ }^{ }\mathit{\mathrm{\mathit{D}}_k^{\mathrm{T}}}$ is thermal diffusivity coefficient. Calculating the molar mass M_n is based on the molar mass and mass fraction of all particle molecules:

The x_k appearing in equation (15) can be calculated from the mass fraction and average molar mass of the particle:

The electric field of the discharge gap can be determined via the following formula:

According to the formula, $\varepsilon_0$ refers to vacuum permittivity. The dielectric material is MgO with $\varepsilon _{\text{r}}$ = 9.

Argon is chosen as the working gas for simplicity. Four species, electrons, Ar^*, Ar⁺ and ${\mathrm{Ar}}^+_ 2$, are taken into considerations of argon chemical kinetics with detail information being listed. The research involves the chemical reaction in discharge, the related parameters can be seen from table 1.

Table  1.  The related parameters in discharge.

Index	Species	$\varepsilon_{\mathrm{th}}$ (eV)	μpN	DpN (m−1 s−1)
No. 1	$\mathrm{e}$	${\text{0}}$	${\mu _{\text{p}}}N({{\varepsilon }})$	Einstein
No. 2	${\text{A}}{{\text{r}}^*}$	${\text{11}}{\text{.55}}$	No charge	${\text{6}}{\text{.45}} \times {\text{1}}{{\text{0}}^{20}}$
No. 3	${\text{A}}{{\text{r}}^{\text{ + }}}$	${\text{15}}{\text{.76}}$	${\mu _{{\text{A}}{{\text{r}}^ + }}}N(E/N)$	Einstein
No. 4	${\text{Ar}}_{\text{2}}^{\text{ + }}$	${\text{14}}{\text{.5}}$	${\mu _{{\text{Ar}}_2^ + }}N(E/N)$	Einstein
Index	Reaction		Rate coefficient
R1	${\mathrm{e }}+ {\text{Ar}} \to {\text{A}}{{\text{r}}^*} + {\mathrm{e}}$		${k_1}(\varepsilon )$
R2	${\mathrm{e }}+ {\text{A}}{{\text{r}}^*} \to {\text{Ar}} + {\mathrm{e}}$		${k_{\text{2}}}(\varepsilon )$
R3	$\text{e}+\text{A}\text{r}^*\to\text{A}\text{r}^{\text{+}}+\text{e + e}$		${k_{\text{3}}}(\varepsilon )$
R4	$\text{e}+\text{Ar}\to\text{A}\text{r}^{\text{+}}+\text{e + e}$		${k_{\text{4}}}(\varepsilon )$
R5	$\text{A}\text{r}^*+\text{A}\text{r}^*\to\text{A}\text{r}^{\text{+}}+\text{e + Ar}$		${k_{\text{5}}}(\varepsilon ){\text{ = 5}} \times {\text{1}}{{\text{0}}^{{{ - 16}}}}$
R6	${\text{e}} + {\text{A}}{{\text{r}}^ + } \to {\text{A}}{{\text{r}}^{\text{*}}}$		$k_{\text{6}}(\varepsilon)\text{ = 4}\times\text{1}\text{0}^{-19}T_{\mathrm{e}}^{-0.5}$
R7	${\text{e}} + {\text{e}} + {\text{A}}{{\text{r}}^ + } \to {\text{A}}{{\text{r}}^{\text{*}}} + {\text{e}}$		$k_{\text{7}}(\varepsilon)\text{ = 5}\times\text{1}\text{0}^{-33}T_{\mathrm{e}}^{-\text{4}\text{.5}}$
R8	${\text{e}} + {\text{Ar}}_2^ + \to {\text{A}}{{\text{r}}^{\text{*}}}{\text{ + Ar}}$		$k_{\text{8}}(\varepsilon)\text{ = 5}\text{.38}\times\text{1}\text{0}^{-14}T_{\mathrm{e}}^{-\text{0}\text{.66}}$
R9	${\text{Ar + Ar}} + {\text{Ar}}_{}^ + \to {\text{A}}{{\text{r}}^{}}{\text{ + Ar}}_2^ +$		${k_{\text{9}}}(\varepsilon ){\text{ = 2}}{\text{.5}} \times {\text{1}}{{\text{0}}^{{{ - 37}}}}$

 | Show Table

DownLoad: CSV

In this work, the diameter of the dielectric layer is set to 0.01 m, the initial electron density is 10¹³ m⁻³ and the dual frequency of the driving source is set to f_L = 50 kHz and f_H = 5 MHz, respectively. The gas gap in two dielectric layer is 2 mm.

The calculation time step is 5×10⁻¹². The narrow area grid resolution is 1 and the grid is divided into 200 units, which satisfies the Debye length. The number of the grid divided in 200 units is of high quality. It could solve the sharp gradient of electron and ion density in the gas gap well.

3.   Results and discussions

The present research aimed to reveal the amplitude effect on electron dynamic behavior and particles’ generation during the discharge process, Thereafter, we analyze the electron dynamics and the method of controlling different modes.

Figure 1 presents the different amplitude components of the LF and RF component voltage amplitude, which determine the different waveforms of applied voltage, resulting in significant differences in the current waveform of air gap discharge. In figure 1 the evolution of the applied voltage shows a trend of LF voltage components dominating the voltage amplitude characteristics when V_H = 150 V as constant and V_L = 300 V and 450 V, respectively. The discharge current waveform maintains an obvious periodic pulse peak appearance in figure 1. This is caused by the dominant LF voltage component and the rapid oscillations generated by the RF component while the RF voltage is 150 V. Comparing figures 1(a)–(f), it can also be observed that the edge morphology of the sheath layer in the spatio-temporal evolution structure is similar to the spatio-temporal evolution structure of the air gap current and aligns with the electric field distribution.

Figure  1.  The distributions of the applied voltage (blue lines) and discharge current density (red lines). (a) Under the V_H = 150 V, V_L = 150 V condition; (b) under the V_H = 300 V, V_L = 150 V condition; (c) under the V_H = 450 V, V_L = 150 V condition; (d) under the V_H = 150 V, V_L = 900 V condition; (e) under the V_H = 300 V, V_L = 900 V condition; (f) under the V_L = 450 V, V_H = 900 V condition.

DownLoad: Full-Size Img PowerPoint

When RF component or LF component dominates individually, the gas gap discharge current is more similar to the current characteristics during a single low- or high- frequency voltage discharge. When the amplitude of LF voltage dominates, the discharge current waveform maintains the characteristics of the gas gap current when discharging under a single LF voltage, which still has an obvious periodic pulse peak appearance, while also having a clear difference, namely the rapid oscillation caused by the presence of RF components.

3.1   The effect of LF component

In the present section, we analyzed the effect of LF amplitude component on both electron dynamic behavior and the distribution of particle density. From figure 2 we can see that by controlling the dual-frequency power supply frequency to f_L = 50 kHz, changing the voltage of the low-frequency component while maintaining the RF component voltage at 150 V, there are many phenomena observed corresponding to this configuration. From figures 2(a)–(d), it can be seen that in the discharge gas gap, the argon ion density in the bulk region is obviously higher compared to that near the boundary, which decreases to form a sheath, and the change of the particle density is accompanied by RF oscillation characteristics, which are caused by RF modulation effects. The argon ion density shows a peak in figures 2(e)–(h), and when the condition changed to V_H = 150 V, V_L = 900 V, a double-peak structure appears. From figures 2(i)–(l) we can see that the Ar^* density distribution at the edge of the sheath gradually approaches a continuous state from a discontinuous state, with collapse throughout the entire sheath region. This indicates that the distribution of electron density remains consistent with the distribution of Ar⁺ density. However, there are significant differences in the density distributions in the sheath zone, sheath edge, and discharge plasma bulk zone.

Figure  2.  The distributions of (a)–(d) the computed electron density, (e)–(h) the computed argon ion density and (i)–(l) the computed metastable argon density. f_L = 50 kHz, f_H = 5 MHz. (a), (e) and (i) under the V_H = 150 V, V_L = 100 V condition; (b), (f) and (j) under the V_H = 150 V, V_L = 300 V condition; (c), (g) and (k) under the V_H = 150 V, V_L = 500 V condition; (d), (h) and (l) under the V_H = 150 V, V_L = 900 V condition.

DownLoad: Full-Size Img PowerPoint

Moreover, the sheath region is significantly compressed, and the spatio-temporal evolution distribution structure of Ar^* exhibits continuous distribution structure at the sheath boundary, which is closely related to the strengthening of RF modulation effects. Observing the evolution of the electric field intensity in figure 3, it can be seen that this index at the sheath boundary is very high, while that in the discharge plasma bulk zone is relatively low and exhibits a pulse structure. As LF voltage component increases, electric field strength gradually increases. Due to the small RF voltage amplitude component at this moment, there is no significant change in this discharge mode.

Figure  3.  The distributions of the computed electric field. f_L = 50 kHz, f_H = 5 MHz. The RF component is set to V_H = 150 V as a constant, the LF component is set to V_L = 100 V, 300 V, 500 V and 900 V. (a) Under the V_H = 150 V, V_L = 100 V condition; (b) under the V_H = 150 V, V_L = 300 V condition; (c) under the V_H = 150 V, V_L = 500 V condition; (d) under the V_H = 150 V, V_L = 900 V condition.

DownLoad: Full-Size Img PowerPoint

Figure 4 shows the electron impact ionization rate and absorption power profiles changed with the increasing LF voltage. Under the condition of V_H = 150 V as a constant, changing the LF component from V_L = 100 V to 900 V, the electron absorption power near the transient anode is higher compared to that at the cathode. Meanwhile, the distribution changes to a multi-pulse structure in figure 4, and in the case of low RF voltage components, the discharge state transforms to the glow-like discharge.

Figure  4.  The distributions of (a)–(d) the computed electron impact ionization rate and (e)–(h) the absorption power profiles. f_L = 50 kHz, f_H = 5 MHz. The RF component is set to V_H = 150 V as a constant, the LF component is set to V_L = 100 V, 300 V, 500 V and 900 V. (a) and (e) under the V_H = 150 V, V_L = 100 V condition; (b) and (f) under the V_H = 150 V, V_L = 300 V condition; (c) and (g) under the V_H = 150 V, V_L = 500 V condition; (d) and (h) under the V_H = 150 V, V_L = 900 V condition.

DownLoad: Full-Size Img PowerPoint

Increasing the RF component to 300 V, the distributions of particle density are shown in figure 5. According to figure 5, the maximum electron density is located near the sheath region. The electron density in the bulk plasma zone remains at high level. At the same time, argon ions have a similar spatio-temporal distribution trend.

Figure  5.  The distributions of (a)–(d) the computed electron density, (e)–(h) the computed argon ion density and (i)–(l) the computed metastable argon density. f_L = 50 kHz, f_H = 5 MHz. (a), (e) and (i) under the V_H = 300 V, V_L = 100 V condition; (b), (f) and (j) under the V_H = 300 V, V_L = 300 V condition; (c), (g) and (k) under the V_H = 300 V, V_L = 500 V condition; (d), (h) and (l) under the V_H = 300 V, V_L = 900 V condition.

DownLoad: Full-Size Img PowerPoint

From the observed graph, it can be seen that when the LF voltage component rises to the dominant position, electrons and Ar⁺ are mainly distributed at the edge of sheath and discharge plasma bulk region. For the expansion of the discharge plasma bulk zone, the electron density presents a symmetrical saddle shaped distribution in that region. The maximum density distribution area is located at both ends of that region, and there is a periodic strength variation with the applied parameters in this area. However, due to resolution reasons, it is not particularly obvious. It can be considered that when the amplitude of the dual frequency voltage is comparable, although bulk heating dominates, it is not uniformly heated. Compared to figure 2 of V_H = 150 V, under the condition that the RF voltage component enhances to 300 V, metastable argon density nearly remains in a continuous distribution state, mainly in the sheath region.

The electric field decreases from the cathode to the anode, as shown in figure 6. When the RF voltage is set to 300 V, it can be seen that when the LF voltage is below 500 V, the computed electron field profiles present a RF dominant distribution characteristic. When the LF voltage increases to 500 V or above, the electric field distribution near the sheath exhibits a pulse peak structure, similar to the LF dominant discharge characteristics. Comparing the distribution when the RF voltage component is 150 V, it is observed that RF voltage component increases.

Figure  6.  The distributions of (a)–(d) the computed electron field profiles. f_L = 50 kHz, f_H = 5 MHz. The RF component is set to V_H = 300 V as a constant, the LF component is set to V_L = 100 V, 300 V, 500 V and 900 V. (a) Under the V_H = 300 V, V_L = 100 V condition; (b) under the V_H = 300 V, V_L = 300 V condition; (c) under V_H = 300 V, V_L = 500 V condition; (d) under the V_H = 300 V, V_L = 900 V condition.

DownLoad: Full-Size Img PowerPoint

The distributions of electron absorption power are shown in figure 7. From figure 7 it can be indicated that the ionization rate changes with increasing the LF voltage amplitude. From absorption power, the ionization rate and electron absorption power profiles occur near the plasma bulk region, which presents the transition dominated by RF component characteristics resembling LF characteristics. Figure 7 shows that the discharge plasma bulk region expands and further compresses the plasma sheath region, and the LF modulation effect becomes more obvious. At the edge of the sheath region, there is a higher electron density in every half LF cycle, and RF oscillations are observed, gradually becoming apparent with the increase of LF voltage. This corresponds to the synergistic effect of RF and LF components in the electric field distribution.

Figure  7.  The distributions of (a)–(d) the computed electron impact ionization rate and (e)–(h) the computed electron absorption power profiles. f_L = 50 kHz, f_H = 5 MHz. The RF component is set to V_H = 300 V as a constant, the LF component is set to V_L = 100 V, 300 V, 500 V and 900 V. (a) and (e) under the V_H = 300 V, V_L = 100 V condition; (b) and (f) under the V_H = 300 V, V_L = 300 V condition; (c) and (g) under the V_H = 300 V, V_L = 500 V condition; (d) and (h) under the V_H = 300 V, V_L = 900 V condition.

DownLoad: Full-Size Img PowerPoint

3.2   The effect of HF component

In figure 8 the computed metastable argon density has a tendency towards consistency. The spatio-temporal evolution distribution structure presents a continuous distribution structure at the sheath boundary, which is closely related to the strengthening of RF modulation effects when RF voltage is set at 450 V. Compared with comparable dual frequency, when the RF component dominates, the uniformity of particle density in the discharge plasma bulk zone is better, but the overall particle density presents a downward trend, which is related to the stronger RF modulation effect. A large number of electrons are heated not only in the sheath but also in the plasma bulk region. There is a high density of electrons and ions in the plasma bulk region.

Figure  8.  The distributions of (a)–(d) the computed electron density, (e)–(h) the computed argon ion density and (i)–(l) the computed metastable argon density. (a), (e) and (i) under the V_L = 300 V, V_H = 150 V condition; (b), (f) and (j) under the V_L = 300 V, V_H = 450 V condition; (c), (g) and (k) under the V_L = 900 V, V_H = 150 V condition; (d), (h) and (l) under the V_L = 900 V, V_H = 450 V condition.

DownLoad: Full-Size Img PowerPoint

Under the condition that the LF component is respectively set to V_L = 300 V and V_L = 900 V, the evolution distributions of particle density varying with RF component are shown in figure 8. From figure 8 we can observe that with the increase of RF components, the discharge tends to be stable when V_H = 450 V. When the LF components is set to 150 V, the electron density is concentrated on the bulk region, while the argon ion density spreads to the sheath region. This indicates that the metastable argon density shows collapse due to electron heating being compressed. Additionally, as the voltage increases, the additional energy is applied to heat the metastable particles, not the ions.

In figure 9, the electron field is similar to the condition when the low-voltage is set to 300 V, but the pulse structure seems to be much larger with the increase of voltage components. The electric field indicates that when changing the dominant power of the applied voltage, it can present different discharge state. It is in line with the mode change proposed in the reference [28], and in figure 10 the relevant impact ionization rate and electron absorption power profiles present significant differences when V_H = 150 V. Also, when the RF component is set to 150 V, the impact ionization and electron absorption power are all manifested in the form of pulses in the plasma bulk region. When RF component increases to 450 V, the impact ionization rate and the absorption power results are dense near the sheath layer. When the electron density enhances to the corresponding threshold, the metastable argon density is mostly in sheath.

Figure  9.  The distributions of the computed electron field. f_L = 50 kHz, f_H = 5 MHz. The RF component is set to V_H = 150 V and 450 V, the LF component is set to V_L = 300 V and 900 V. (a) Under the V_L = 300 V, V_H = 150 V condition; (b) under the V_L = 300 V, V_H = 450 V condition; (c) under the V_L = 300 V, V_H = 500 V condition; (d) under the V_L = 300 V, V_H = 900 V condition.

DownLoad: Full-Size Img PowerPoint

Figure  10.  The distributions of (a)–(d) the computed electron impact ionization rate and (e)–(h) the computed electron absorption power profiles. f_L = 50 kHz, f_H = 5 MHz. The RF component is set to V_H = 150 V and 450 V, the LF component is set to V_L = 300 V and 900 V. (a) and (e) under the V_L = 300 V, V_H = 150 V condition; (b) and (f) under the V_L = 300 V, V_H = 450 V condition; (c) and (g) under the V_L = 900 V, V_H = 150 V condition; (d) and (h) under the V_L = 900 V, V_H = 450 V condition.

DownLoad: Full-Size Img PowerPoint

The distributions of electron absorption power under different LF/RF components are shown in figure 10. From figure 10 we can observe that when high frequency dominates, the electrons density distribution exhibits the same spatio-temporal evolution characteristics. When RF dominates, the electrons have no time to respond to the rapid oscillations of with the RF electric field and relevant process. Electrons were limited within the bulk region and oscillate back and forth with the changes of the electric field, which can provide a larger initial electron density for the next discharge period and form a positive feedback effect under certain conditions. On the contrary, when LF component dominates, especially when V_L = 900 V, the periodic enhancement of electric field intensity distribution is very obvious, showing discharge characteristics of low-frequency discharge. The electric field strength has significantly increased, and the spatio-temporal evolution structure at the edge of the sheath has become more uniform and smoother. Electrons are unable to keep pace with the rapid changes of the electric field and stay in the discharge plasma bulk region, resulting in a higher electron density distribution in that region. Within the sheath region, the particle density is lower and sheath is wider. With the increase of the applied voltage, more electrons are bound in the plasma bulk region, which is conducive to the accumulation of ions in the sheath. Because the electric field of the sheath enhances, the electron absorption power and the electron impact ionization rate increase.

4.   Conclusions

In summary, a DBD system driven by dual LF/RF excitation was studied in present research. Our findings reveal that the distributions of particle densities—specifically electrons, argon ions, and metastable argon ions—vary with different LF/RF applied voltage amplitudes. Electron heating is a crucial factor influencing the excitation mode. When the RF component is kept constant while the LF component varies from 100 to 900 V, the electron density distribution transitions from a pulse-like structure to a continuous structure, with electron heating primarily occurring in the plasma bulk region. As the RF component increases, electron heating begins to compress toward the sheath regions. When the RF component predominates, the peak values of electron heating and electron density are observed in the sheath region. Additionally, the electric field structure shifts from a multi-pulse pattern to a continuous one. This change is closely linked to the inability of electrons to follow the electric field due to strong collisions. These results enhance our understanding of dual LF/RF modulated atmospheric pressure DBDs and have the potential to expand their applications in various fields.