Gas Heating Mechanisms in Atmospheric Pressure Helium Dielectric-Barrier Discharges Driven by a kHz Power Source

The present work investigates contributions of different heating mechanisms and power efficiency of atmospheric-pressure helium dielectric-barrier discharges (APHeDBDs) containing a small amount of N2 for temperature measurements by developing the numerical methodology combining the one-dimensional (1D) plasma fluid model (PFM) and 3D gas flow model (GFM) with simulated results validated by measurements including the discharge power consumption and temperature distribution. The discharge dynamics are modeled by the 1D PFM for evaluating the average heating source considering elastic collision, ion Joule heating, and exothermic reactions as the source term of energy equation solved in the 3D GFM. The simulated current density reaches 29 A m−2 which is close to that measured as 35 A m−2. The simulated power consumption is 2.0 W which is in good agreement with the average measured power consumption as 2.1 W. The simulated average gas temperature in the reactive zone is around 346 K which is also close to the rotational temperature determined. The analysis shows that elastic collision and ion Joule heating are dominant heating mechanisms contributing 23.9% and 65.8% to the heating source, respectively. Among ion species, N2+ and N4+ are dominant species contributing 44.1% and 50.7% to the heating source of ion Joule heating, respectively. The simulated average total heating source is around 5.6 × 105 W m−3 with the maximum reaching 3.5 × 106 W m−3 in the sheath region due to the contribution of ion Joule heating.


Introduction
Atmospheric pressure helium dielectric-barrier discharges (APHeDBDs) have been developed intensively in different fields including surface treatments [1][2][3] and plasma medicine [4][5][6] due to the efficient generation of reactive species and stable operating conditions. Typically, APHeDBDs are driven by kHz power sources since they are simple and can be scaled up for applications easily. Therefore, APHeDBDs driven by kHz power sources can be considered as one of the important plasma sources for the aforementioned applications.
It is common to consider gas temperature as one of critical discharge parameters since gas temperature influences rate constants of chemical reactions and densities of background gases, determining the discharge chemistry. Gas temperature is also the key parameter affecting transport properties (i.e., diffusivities), which may lead to different transport behavior. Moreover, it is critical to control the gas temperature for applications with materials that are sensitive to heating load. Hence, it is essential to control the gas temperature of APHeDBDs.
The gas temperature of a helium discharge is typically obtained by determining the rotational temperature from the emission spectrum of N 2 second positive system (SPS, N 2 (C 3 Π u )→N 2 (B 3 Π g )) measured via an optical emission spectrometer (OES) with the addition of low N 2 concentration (e.g.,~0.1%) [7][8][9]. The helium discharge can be sustained at a wide range of gas temperatures ranging from room temperature (e.g., 300 K) to a few hundred degrees Kelvin [9][10][11], depending on the operating conditions. Therefore, it is necessary to control and determine the gas temperature properly. From the application perspective, it is important to realize the heating mechanisms for optimizing an appropriate reactor and controlling at proper operating conditions. However, the measurements of the rotational temperatures at various operating conditions provide no further information on heating mechanisms of complex discharge behavior. Fortunately, numerical simulations provide another alternative to reveal details of discharge dynamics.
APHeDBDs have been studied numerically for years due to their importance for applications [12][13][14][15][16]. In general, the simulated results were obtained efficiently with the one-dimensional (1D) plasma fluid model (PFM) by solving the species continuity equations, the Poisson equation, and the electron energy density equation to model the spatial and temporal changes of species densities, electric field, and electron temperature, respectively, for capturing the complex discharge dynamics. Although APHeDBDs can be operated as uniform glow-like discharges that are modeled appropriately by the 1D PFM, it was identified that APHeDBDs can also be operated as the patterned non-uniform discharges [17,18]. Hence, it is important to examine the discharge uniformity before employing the 1D PFM to model the APHeDBDs. Moreover, it was reported that a small amount of N 2 (e.g., a few ppm) may play the dominant role in APHeDBDs due to Penning ionization and associative processes among excited helium species, helium ions, and N 2 molecules [12,15,16]. Therefore, it is essential to model APHeDBDs with the proper concentration of N 2 added for temperature measurements.
The mechanisms of gas heating in atmospheric pressure helium discharges have been studied for years [9,[19][20][21][22]. The gas heating effects of atmospheric pressure helium microplasmas sustained by the dc power sources were investigated with either the 1D or 2D PFM in conjunction with the 2D gas flow model (GFM) considering ion Joule heating as the heating source for modeling the gas temperature in the discharge [9,19]. It was identified that most (~87%) of the total input power goes to gas heating and the gas temperature increases to~650 K due to the ion Joule heating. Kang et al. explored numerically the gas temperature effect of an APHeDBD excited by a kHz sinusoidal voltage with the 1D PFM accompanied by the energy equation of the background gas considering the ion Joule heating and exothermic reactions as heating mechanisms [20]. It was reported that the gas temperature plays a significant role in the discharge mode. Recently, the gas heating effects of atmospheric pressure helium discharges driven by radio-frequency power sources were studied by solving the 1D PFM with the energy equation of the background gas considering elastic collision, ion Joule heating, and exothermic reactions as heating mechanisms for obtaining the gas temperature distribution in the discharge [21,22]. It was shown that elastic collision is the dominant heating mechanism in the bulk region while the ion Joule heating and exothermic reactions are dominant in the sheath region.
Although the gas heating effects of helium discharges have been studied numerically for various configurations with different operating conditions, detailed thermal characterization of APHeDBDs generated in a planar reactor driven by the kHz power source remains rare and insufficient for gaining a full understanding. The present work investigates the gas heating behavior of APHeDBDs generated in the planar reactor driven by a kHz power source with both numerical simulations and experimental measurements. The APHeDBDs are modeled by the 1D PFM with the helium chemistry considering the addition of N 2 gas (used for the measurement of the rotational temperature) for evaluating heating sources of different mechanisms including elastic collision, ion Joule heating, and exothermic reactions. The discharge uniformity is confirmed by the intensified CCD camera (ICCD) to justify the assumption Appl. Sci. 2020, 10, 7583 3 of 20 of uniform discharge for using the 1D PFM to simulate the APHeDBDs. The average heating source calculated by the 1D PFM is provided for the energy equation solved in the 3D GFM to model the gas temperature of the reactor. The simulated cyclic average power consumption of the APHeDBD is validated with experimental measurements. The rotational temperature in the reactive zone is determined and the surface temperature of the reactor is detected by the infrared (IR) thermal imager to validate the 3D GFM. The power efficiency is studied and the contributions of different heating mechanisms are presented. Section 2 introduces the experimental configuration and the numerical method developed. The model validation is presented in the first part of Section 3 and the heating behavior is investigated as the second part of this section. Section 4 concludes based on the major findings of this work.

Experimental Configuration
The experimental configuration used in this work is shown in Figure 1a. A planar reactor is designed with parallel glass plates as transparent electrodes for generating APHeDBDs. One of the surfaces of each glass plate is coated with indium tin oxide as the conductive layer. The reactive area is a circle with a diameter of 30 mm and the discharge gap is fixed at 5 mm. The sinusoidal external power source is supplied with the voltage frequency and amplitude controlled at 16.5 ± 0.1 kHz and 3.5 ± 0.1 kV, respectively. The helium flowrate is controlled at 2 slm and a small amount of N 2 admixture (0.1%) is added for measuring the rotational temperature. The discharge voltages and currents are measured by the high-voltage probe (Tektronix P6015A) and current probe (Magnelab CT-C1), respectively, with data stored by the oscilloscope (Keysight Technologies DSOX3054T, 500 MHz and 5.0 GSa s −1 ). A monitor capacitor with the capacitance of 230 pF is installed at the location after the reactor to measure charge transferred for determining the power consumption by integrating the Lissajous curve [23]. To confirm the discharge uniformity, an ICCD camera (Princeton Instruments PI-MAX4) with an ultraviolet (UV) lens (Nikon PF10545MF-UV) is applied to take photographs from the reactor outlet during the breakdown phase. An OES (AvaSpec ULS2048LUSB2) with the optical fiber located at the lateral side of the reactor as shown in Figure 1b is employed for measuring the average emission spectrum to determine the rotational temperature. The infrared thermal imager (AVIO F50A) is applied to detect the surface temperature of the reactor from the lateral side of the reactor (similar to that shown in Figure 1b) for validating the temperature distribution simulated by the 3D GFM. The gas temperature is measured after the system reaches its thermal equilibrium as the steady-state gas temperature.

Numerical Method
The 1D PFM and 3D GFM are developed to model the discharge dynamics and steady-state reactor temperature for investigating the thermal behavior of APHeDBDs ignited within the planar reactor. The discharge is simplified as the 1D PFM with its uniformity confirmed by the ICCD photograph taken. The average volumetric heating source calculated by the 1D PFM is adopted as the source term of the energy equation considered in the 3D GFM to model the reactor temperature. The simulation will be performed until the gas temperature used in the 1D PFM is consistent with that obtained by the 3D GFM as described later.

One-Dimensional Plasma Fluid Model (1D PFM)
This section briefs the 1D PFM developed previously [24] with the framework of discharge chemistry including helium and nitrogen related species. Figure 2a shows the numerical domain of the 1D PFM. The discharge gap is 5.0 mm and the thickness of each glass plate is 0.7 mm as used in experimental measurements. The non-uniform mesh is employed with the mesh size covering Appl. Sci. 2020, 10, 7583 4 of 20 from 30 µm to 100 µm for capturing discharge dynamics appropriately in the sheath and bulk regions. The mesh independence test has been conducted to ensure the mesh quality.

Numerical Method
The 1D PFM and 3D GFM are developed to model the discharge dynamics and steady-state reactor temperature for investigating the thermal behavior of APHeDBDs ignited within the planar reactor. The discharge is simplified as the 1D PFM with its uniformity confirmed by the ICCD photograph taken. The average volumetric heating source calculated by the 1D PFM is adopted as the source term of the energy equation considered in the 3D GFM to model the reactor temperature. The simulation will be performed until the gas temperature used in the 1D PFM is consistent with that obtained by the 3D GFM as described later.

One-Dimensional Plasma Fluid Model (1D PFM)
This section briefs the 1D PFM developed previously [24] with the framework of discharge chemistry including helium and nitrogen related species. Figure 2a shows the numerical domain of the 1D PFM. The discharge gap is 5.0 mm and the thickness of each glass plate is 0.7 mm as used in experimental measurements. The non-uniform mesh is employed with the mesh size covering from The continuity equation of each species considered is solved as: where n is the species density, Г is the species flux, and the subscripts e, i, and uc represent electron, ion, and neutral species, respectively. The source term S of each continuity equation is evaluated from the discharge chemistry introduced next to account for the generation and destruction reactions of each species. The species flux terms calculated with the drift-diffusion approximation are written as: Γ e,i = sign(q e,i )µ e,i n e,i E − D e,i ∂n e,i ∂x (2) where q is the species charge, µ is the species mobility, E (= − ∂ϕ ∂x , ϕ is the potential) is the electric field, and D is the species diffusivity. from the experimental cross section data.
The discharge heating source is obtained from the integration of the power released by different heating mechanisms during the cycle as the cyclic average heating source. From the 1D PFM, the cyclic average heating source of the discharge is calculated as: where the subscript "CA" represents cyclic average quantity calculated, tcycle is the period of one cycle, Ni is the number of ion species, Rl is the reaction rate of the lth reaction, ∆ is the energy released from the reaction. The first term of the integrand considers the energy released from elastic collisions. The second term of the integrand considers the ion Joule heating under the assumption that energy acquired by ions from the electric field is transferred thoroughly to background gases due to efficient collisions among ions and molecules of background gases. The last term of the integrand considers the energy released from exothermic reactions listed in Table 1.  The Poisson equation is solved to calculate the electric field built inside the discharge gap as: where ε (= ε 0 ε r , ε 0 is the permittivity of free space and ε r is the relative permittivity of each region) is the permittivity of each region and N c represents the number of charged species. The electron energy density equation is solved to evaluate the electron temperature for determining rate constants of electron-related reactions and electron transport properties (i.e., D e and µ e ) as: 3 m e M k n e k B ν m,k T e − T g − n e N e j=1 ε j k j n j where n ε = 3 2 n e k B T e is the electron energy density, k B is the Boltzmann constant, T e is the electron temperature, m e is the electron mass, M is the mass of the background gas (i.e., He or N 2 ), ν m is the frequency of elastic collision, N B is the number of background gases (i.e., N B = 2 for He and N 2 ), T g is the gas temperature, ε j is the threshold energy of the jth electron-impact reaction, k is the rate constant, and N e is the number of electron-impact reactions considered in the chemical model. On the right-hand side of Equation (5), the first term is the electron Joule heating evaluating the energy acquired by electrons from the externally applied voltage. The second term calculates the energy loss due to elastic collisions among electrons and neutral background gases (i.e., He and N 2 ). The last term calculates the energy consumption due to electron-impact reactions considered in the discharge chemistry.
It is assumed that charged species accumulate and neutral species quench at the glass surface. The instantaneous voltage of the externally applied power source is assigned on the powered electrode and the potential of the grounded electrode is kept at zero as the boundary conditions for solving the Poisson equation.
The rate constants of electron-impact reactions and electron transport properties are calculated before the simulation by using the Boltzmann equation solver (i.e., Bolsig+ [25]) as look-up tables. The mobilities of ions are adopted from the literature [16,26] with diffusivities evaluated from the Einstein relation. The diffusivities of neutral species are evaluated from the binary diffusion theory [27] with relative parameters that can be found in [28].
The reactions considered in the discharge chemistry are listed in Table 1. The charged species considered include electrons, He + , He 2 + , N 2 + , and N 4 + , to model the conduction currents properly through ionization reactions, association reactions, and recombination reactions. The excited species of helium including He m and He 2 * play essential roles for introducing ionization pairs through Penning ionization. The excited species of nitrogen (i.e., are considered to account for the energy released from exothermic reactions. Interactions among helium and nitrogen related species are also included to mimic the realistic discharge chemistry.
The discharge heating source is obtained from the integration of the power released by different heating mechanisms during the cycle as the cyclic average heating source. From the 1D PFM, the cyclic average heating source of the discharge is calculated as: where the subscript "CA" represents cyclic average quantity calculated, t cycle is the period of one cycle, N i is the number of ion species, R l is the reaction rate of the lth reaction, ∆E is the energy released from the reaction. The first term of the integrand considers the energy released from elastic collisions. The second term of the integrand considers the ion Joule heating under the assumption that energy acquired by ions from the electric field is transferred thoroughly to background gases due to efficient collisions among ions and molecules of background gases. The last term of the integrand considers the energy released from exothermic reactions listed in Table 1.

Three-Dimensional Gas Flow Model (3D GFM)
To simulate the steady-stage gas temperature of the reactor, the commercial software CFD-ACE 2015 was used for building the 3D GFM considering the flow dynamics of the reactor. This section summarizes the governing equations modeled in the 3D GFM while detailed descriptions can be found in [33]. Figure 2b shows the numerical domain including the reactive zone, glass plate, Teflon block, and ambient air environment. Due to the reactor symmetry, the 3D GFM is simplified with the symmetric plane assigned to reduce the model size without losing the accuracy.
The steady-state flow dynamics of the background gases is modeled by solving the mass and momentum conservation equations as: where ρ is the density of the background gas (i.e., air or He), → V is the velocity vector, p represents the pressure of the background gas, τ is the viscous stress tensor, and → S M is the momentum source term to calculate the buoyance force resulting in the cooling of natural convection on the reactor surface. He m + He * 2 → e + 2He + He He m + He m → e + He + He m + He * 2 → e + He + He He He The steady-state gas temperature of the reactor is modeled by solving the energy equation (total enthalpy) as: x , h is the enthalpy and u represents the velocity component) is the total enthalpy, λ is the thermal conductivity of the background gas, and S GFM is the average volumetric heating source obtained from the 1D PFM. The S GFM is applied to the reactive zone as the average volumetric heating source evaluated as: where the subscript "SCA" represents the spatial cyclic average quantity calculated, N cell is the cell number used in the 1D PFM for the discharge gap, w is the cell width, d gap is the discharge gap as 5 mm, and the subscript m represents the mth cell.
In the discharge gap, the helium flow is provided as the inlet boundary condition from the top of the reactor. We test whether the small amount of N 2 addition had only a minor effect on the simulated gas temperature; therefore, the small amount of N 2 addition is neglected in the 3D GFM. The outlet boundary condition is assigned on boundaries of ambient air around the reactor to allow air flows in and out of the boundaries as the natural convection is built. The solid boundaries are assigned as zero-gradient condition as: The thermal conductivities of the Teflon block and glass plate are 0.25 W m −1 K −1 and 1.1 W m −1 K −1 , respectively. The temperature-dependent thermal properties of helium and air are evaluated by the kinetic theories which can be found in [33]. This is tested and it is concluded that the simulated temperature distribution is insensitive to the size of the Teflon block. The mesh independence test has been conducted for the 3D GFM to ensure the mesh quality.

Model Coupling
The 1D PFM and 3D GFM are developed to model the discharge dynamics and gas temperature of the reactor separately and combined with the weak coupling as shown in Figure 3 since it is efficient to obtain the individual solution separately for two models with different characteristic times [34]. The simulation is performed from the 1D PFM by adopting the initial gas temperature T g,PFM as 400 K until the quasi-steady solution is obtained. The cyclic average power consumption and average volumetric heating source are monitored as the criteria of the quasi-steady solution and are written as: where the superscripts k and k − 1 are the quantities obtained in the current and previous time steps, respectively. The cyclic average power consumption of the simulated discharge P CA−PFM is evaluated by integrating from the voltage and total current at each time step obtained in the 1D PFM over the whole cycle as: Appl. Sci. 2020, 10, 7583 9 of 20 where I total is the total current including the conduction currents of charged species and the displacement current. Then quasi-steady solution of the 1D PFM is used to calculate the average volumetric heating source which is adopted as the source term of the energy equation considered in the 3D GFM for modeling the temperature of the reactor. After obtaining the steady-state solution of the 3D GFM, the average gas temperature in the reactive zone is calculated as T g,GFM . If the gas temperature obtained from the 3D GFM is not consistent with the gas temperature T g,PFM used in the 1D PFM, the T g,PFM is replaced with the T g,GFM for running the 1D PFM again. The converged solution will be obtained until the T g,PFM used is consistent with the T g,GFM calculated with the following criterion: Generally, two or three iterations are required to obtain the converged solution. The converged solution is analyzed and presented in Section 3.

Results and Discussion
The first part of this section validates the method developed. The discharge uniformity is confirmed with the ICCD photograph to justify the simplification of 1D PFM. The simulated discharge characteristics including the current density and power consumption are compared with experimental measurements. The simulated average gas temperature of the reactive zone is compared with the rotational temperature determined, and the surface temperature of the reactor is compared with that detected by the IR thermal imager. The second part of this section analyzes the dominant heating mechanisms of APHeDBDs. The contribution of each heating mechanism is presented.

Results and Discussion
The first part of this section validates the method developed. The discharge uniformity is confirmed with the ICCD photograph to justify the simplification of 1D PFM. The simulated discharge characteristics including the current density and power consumption are compared with experimental measurements. The simulated average gas temperature of the reactive zone is compared with the rotational temperature determined, and the surface temperature of the reactor is compared with that detected by the IR thermal imager. The second part of this section analyzes the dominant heating mechanisms of APHeDBDs. The contribution of each heating mechanism is presented.

Discharge Uniformity
In this work, the APHeDBD is simplified with the 1D PFM under the assumption of uniform discharge. However, it is noted that APHeDBDs may not be uniform under various operating conditions [17,18]. The discharge uniformity under current operating conditions is confirmed by the ICCD photograph taken during the breakdown of a typical cycle with the 1 µs exposure time as shown in Figure 4. The uniform discharge irradiation occurs near the surface of the glass plate covering the grounded electrode due to excited species generated by energetic electrons in the cathode fall region. Uniform discharges generated in the planar reactors under similar operating conditions were reported by other groups with similar discharge structure presented [16,17], supporting the uniform discharges observed in this work. The confirmed discharge uniformity justifies the application of the 1D PFM for modeling the APHeDBDs.

Power Consumption
The power consumption measured is obtained by integrating the Lissajous curve measured in a typical cycle from the capacitor as shown in Figure 5. It is noted that only a few data points are plotted for presentation purposes and the solid line is regressed from measured data to calculate the power consumption. The average power consumption measured from measurements is around 2.1 W. Figure 6 shows the comparison between the measured and simulated current densities of a typical cycle. It is noted that the experimental current densities in the positive half period (HP) are similar to those in the negative HP, which shows that the discharge is symmetrical under current operating conditions. The simulated current density reaches around 29 A m −2 which agrees with that measured on average as 35 A m −2 . The simulated cyclic average power consumption can be calculated from the integration of the externally applied voltage and simulated current as 2.0 W which agrees with the measured average power consumption as 2.1 W.

Power Consumption
The power consumption measured is obtained by integrating the Lissajous curve measured in a typical cycle from the capacitor as shown in Figure 5. It is noted that only a few data points are plotted for presentation purposes and the solid line is regressed from measured data to calculate the power consumption. The average power consumption measured from measurements is around 2.1 W. Figure 6 shows the comparison between the measured and simulated current densities of a typical cycle. It is noted that the experimental current densities in the positive half period (HP) are similar to those in the negative HP, which shows that the discharge is symmetrical under current operating conditions. The simulated current density reaches around 29 A m −2 which agrees with that measured on average as 35 A m −2 . The simulated cyclic average power consumption can be calculated from the integration of the externally applied voltage and simulated current as 2.0 W which agrees with the measured average power consumption as 2.1 W. typical cycle. It is noted that the experimental current densities in the positive half period (HP) are similar to those in the negative HP, which shows that the discharge is symmetrical under current operating conditions. The simulated current density reaches around 29 A m −2 which agrees with that measured on average as 35 A m −2 . The simulated cyclic average power consumption can be calculated from the integration of the externally applied voltage and simulated current as 2.0 W which agrees with the measured average power consumption as 2.1 W.   The simulated and measured current densities of atmospheric-pressure helium dielectric-barrier discharges (APHeDBDs) generated in the planar reactor. Figure 7 shows the average emission spectrum of N2 SPS measured from the reactive zone for determining the rotational temperature. The red solid line labeled as "Trot" is the theoretical spectrum calculated from the commercial software Specair with theories of quantum mechanics for comparing with the average measured emission spectrum to determine the rotational temperature. The rotational temperature determined is 350 K. The simulated and measured current densities of atmospheric-pressure helium dielectric-barrier discharges (APHeDBDs) generated in the planar reactor. Figure 7 shows the average emission spectrum of N 2 SPS measured from the reactive zone for determining the rotational temperature. The red solid line labeled as "T rot " is the theoretical spectrum calculated from the commercial software Specair with theories of quantum mechanics for comparing with the average measured emission spectrum to determine the rotational temperature. The rotational temperature determined is 350 K. Figure 7 shows the average emission spectrum of N2 SPS measured from the reactive zone for determining the rotational temperature. The red solid line labeled as "Trot" is the theoretical spectrum calculated from the commercial software Specair with theories of quantum mechanics for comparing with the average measured emission spectrum to determine the rotational temperature. The rotational temperature determined is 350 K.  The simulated heating power calculated by the 1D PFM in the reactive zone is around 1.8 W.

Reactor Temperature
In other words, around 90% of the average power consumption contributes to gas heating under current operating conditions, which is similar to the simulated result obtained for a helium microplasma operated under atmospheric pressure [9]. The volumetric heating source calculated by Equation (10) is adopted as the source term of the energy equation considered in the 3D GFM for modeling the temperature distribution of the reactor. Figure 8 shows the gas temperature in the symmetric plane of the reactor. The helium gas flows from the inlet on the reactor top and exits from the outlet on the reactor bottom. The helium flow is supplied at room temperature (300 K) and is heated as the flow passes the reactive zone. It is observed that the high-temperature zone moves slightly downward from the central region of the reactive zone because of the convective effect with the highest gas temperature reaching around 373 K. The average simulated gas temperature of the reactive zone is around 346 K which is close to the rotational temperature determined.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 21 The simulated heating power calculated by the 1D PFM in the reactive zone is around 1.8 W. In other words, around 90% of the average power consumption contributes to gas heating under current operating conditions, which is similar to the simulated result obtained for a helium microplasma operated under atmospheric pressure [9]. The volumetric heating source calculated by Equation (10) is adopted as the source term of the energy equation considered in the 3D GFM for modeling the temperature distribution of the reactor. Figure 8 shows the gas temperature in the symmetric plane of the reactor. The helium gas flows from the inlet on the reactor top and exits from the outlet on the reactor bottom. The helium flow is supplied at room temperature (300 K) and is heated as the flow passes the reactive zone. It is observed that the high-temperature zone moves slightly downward from the central region of the reactive zone because of the convective effect with the highest gas temperature reaching around 373 K. The average simulated gas temperature of the reactive zone is around 346 K which is close to the rotational temperature determined. To further validate the 3D GFM developed, the simulated surface temperature of the reactor is compared with that detected by the IR thermal imager as presented in Figure 9a and b, respectively. Similarly, the high-temperature zone of the simulated surface temperature moves downward from the central region of the reactive zone as explained previously, which agrees with the temperature distribution detected by the IR thermal imager. The simulated highest surface temperature is 362 K, To further validate the 3D GFM developed, the simulated surface temperature of the reactor is compared with that detected by the IR thermal imager as presented in Figure 9a,b, respectively. Similarly, the high-temperature zone of the simulated surface temperature moves downward from the central region of the reactive zone as explained previously, which agrees with the temperature distribution detected by the IR thermal imager. The simulated highest surface temperature is 362 K, which is in good agreement with that detected by the IR thermal imager as 365 K with an uncertainty of 2 K.

Analysis of Atmospheric-Pressure Helium Dielectric-Barrier Discharges (APHeDBDs) Gas Heating
To understand the gas heating behavior of APHeDBDs, the simulated cyclic average power consumption is analyzed and the contributions of different heating mechanisms are presented. Generally, the simulated results capture the thermal behavior of the APHeDBDs.

Analysis of Atmospheric-Pressure Helium Dielectric-Barrier Discharges (APHeDBDs) Gas Heating
To understand the gas heating behavior of APHeDBDs, the simulated cyclic average power consumption is analyzed and the contributions of different heating mechanisms are presented. Figure 10 shows the details of the simulated spatial cyclic average power consumption of different mechanisms considered in the 1D PFM. The cyclic average power consumption of the discharge is equal to the power acquired by charged species as:

Analysis of Power Consumption
where S SCA−ElectronJoule is the spatial cyclic average electron Joule heating and S SCA−IonJoule is the spatial cyclic average ion Joule heating. Moreover, the S SCA−ElectronJoule can be further categorized as: where S SCA−Elastic is the heating source due to elastic collisions among electrons and background gases, S SCA−Exothermic is the heating power released from exothermic reactions considered in the discharge chemistry, and S SCA−Others is the power consumed for generating other reactive species in the discharge. Figure 10 shows that around 90% of the cyclic average power consumption P CA−PFM goes to gas heating through elastic collision, ion Joule heating, and exothermic reactions considered. It is observed that elastic collision and ion Joule heating are the dominant heating mechanisms contributing 23.9% and 65.8% of the power consumption, respectively, to the heating source, which is different from the simulated results obtained in [21,22] for helium discharges excited by the radio-frequency power sources showing that elastic collision is the dominant heating source. Although a few exothermic reactions release a large amount of energy in each reaction, the exothermic reactions considered contribute only a minor amount of energy to the heating source due to the low reaction rates.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 21 Figure 10 shows the details of the simulated spatial cyclic average power consumption of different mechanisms considered in the 1D PFM. The cyclic average power consumption of the discharge is equal to the power acquired by charged species as:

Analysis of Power Consumption
where is the spatial cyclic average electron Joule heating and is the spatial cyclic average ion Joule heating. Moreover, the can be further categorized as: where is the heating source due to elastic collisions among electrons and background gases, is the heating power released from exothermic reactions considered in the discharge chemistry, and is the power consumed for generating other reactive species in the discharge. Figure 10 shows that around 90% of the cyclic average power consumption goes to gas heating through elastic collision, ion Joule heating, and exothermic reactions considered. It is observed that elastic collision and ion Joule heating are the dominant heating mechanisms contributing 23.9% and 65.8% of the power consumption, respectively, to the heating source, which is different from the simulated results obtained in [21,22] for helium discharges excited by the radio-frequency power sources showing that elastic collision is the dominant heating source. Although a few exothermic reactions release a large amount of energy in each reaction, the exothermic reactions considered contribute only a minor amount of energy to the heating source due to the low reaction rates.  Figure 11 shows the contributions of elastic collisions among electrons and different background neutral gases. Although the energy exchange of each elastic collision process between an electron and a background molecule is inefficient due to large mass difference, the high collision  Figure 11 shows the contributions of elastic collisions among electrons and different background neutral gases. Although the energy exchange of each elastic collision process between an electron and a background molecule is inefficient due to large mass difference, the high collision frequency results in a remarkable contribution to the heating source as reported by other groups [9,21,22] in atmospheric pressure helium discharges. The collisions of electrons and helium atoms contribute to more than 99.9% of the heating source of elastic collision while the collisions of electrons and N 2 molecules contribute to less than 0.1% of the heating source of elastic collision since the concentration of helium atoms is much higher than that of the N 2 molecules.

Elastic Collision
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 21 frequency results in a remarkable contribution to the heating source as reported by other groups [9,21,22] in atmospheric pressure helium discharges. The collisions of electrons and helium atoms contribute to more than 99.9% of the heating source of elastic collision while the collisions of electrons and N2 molecules contribute to less than 0.1% of the heating source of elastic collision since the concentration of helium atoms is much higher than that of the N2 molecules. Figure 11. The contribution of each background gas to the heating source of elastic collision.

Ion Joule Heating
The contribution of each species to the ion Joule heating is shown in Figure 12. It is assumed that the power acquired by the ion species is transferred to background gases thoroughly due to the efficient collisions among ion species and background gases with comparable masses. Among ion species, the species N2 + and N4 + contribute 44.1% and 50.7%, respectively, to the heating source of ion Joule heating as the dominant species though the concentration of N2 is much lower than that of helium.
To further study the contribution of each ion species, Figure 13 shows the spatial average density changes of ion species during the breakdown in a typical positive HP. For the presentation purpose, the horizontal axis is rescaled. As the breakdown occurs, the number densities of He + and He2 + increase simultaneously due to: It is observed that the number density of N2 + increases faster than those of He + and He2 + because of the high rate constant of Penning ionization: As the density of N2 + increases, the density of N4 + increases due to the fast associative process: As a result, the species N4 + is the dominant ion species. Similar simulated results can be found in [13,15,16] showing that N4 + is the dominant ion species if the concentration of N2 admixture is higher than a certain level (e.g., ~tens of ppm) in APHeDBDs. Therefore, the species N4 + has the highest contribution to the ion Joule heating.

Ion Joule Heating
The contribution of each species to the ion Joule heating is shown in Figure 12. It is assumed that the power acquired by the ion species is transferred to background gases thoroughly due to the efficient collisions among ion species and background gases with comparable masses. Among ion species, the species N 2 + and N 4 + contribute 44.1% and 50.7%, respectively, to the heating source of ion Joule heating as the dominant species though the concentration of N 2 is much lower than that of helium.
To further study the contribution of each ion species, Figure 13 shows the spatial average density changes of ion species during the breakdown in a typical positive HP. For the presentation purpose, the horizontal axis is rescaled. As the breakdown occurs, the number densities of He + and He 2 + increase simultaneously due to: e + He → 2e + He + (R4) 2He + He + → He + He + 2 (R17) It is observed that the number density of N 2 + increases faster than those of He + and He 2 + because of the high rate constant of Penning ionization: As the density of N 2 + increases, the density of N 4 + increases due to the fast associative process: As a result, the species N 4 + is the dominant ion species. Similar simulated results can be found in [13,15,16] showing that N 4 + is the dominant ion species if the concentration of N 2 admixture is higher than a certain level (e.g.,~tens of ppm) in APHeDBDs. Therefore, the species N 4 + has the highest contribution to the ion Joule heating.

Spatial Distributions of Heating Sources
Before investigating the spatial distributions of heating sources, it is essential to understand the cyclic average spatial distributions of species densities of electrons and dominant ions as shown in Figure 14 with the distribution of electron temperature included. It is noted that the sheath regions are built near the glass surfaces (i.e., at 0.7 mm and 5.7 mm) since the intensive electric fields (not shown) are built due to the distribution of charged species. The intensive electric field that is built results in energetic electrons with high electron temperature, which leads to high reaction rates in the sheath region. Similar simulated results with sheath regions were observed for atmospheric pressure helium dielectric barrier discharges excited by kHz power sources [16,35]. It is observed that the electron temperature reaches around 3.9 eV as the maximum in the sheath region, while the

Spatial Distributions of Heating Sources
Before investigating the spatial distributions of heating sources, it is essential to understand the cyclic average spatial distributions of species densities of electrons and dominant ions as shown in Figure 14 with the distribution of electron temperature included. It is noted that the sheath regions are built near the glass surfaces (i.e., at 0.7 mm and 5.7 mm) since the intensive electric fields (not shown) are built due to the distribution of charged species. The intensive electric field that is built results in energetic electrons with high electron temperature, which leads to high reaction rates in the sheath region. Similar simulated results with sheath regions were observed for atmospheric pressure helium dielectric barrier discharges excited by kHz power sources [16,35]. It is observed that the electron temperature reaches around 3.9 eV as the maximum in the sheath region, while the

Spatial Distributions of Heating Sources
Before investigating the spatial distributions of heating sources, it is essential to understand the cyclic average spatial distributions of species densities of electrons and dominant ions as shown in Figure 14 with the distribution of electron temperature included. It is noted that the sheath regions are built near the glass surfaces (i.e., at 0.7 mm and 5.7 mm) since the intensive electric fields (not shown) are built due to the distribution of charged species. The intensive electric field that is built results in energetic electrons with high electron temperature, which leads to high reaction rates in the sheath region. Similar simulated results with sheath regions were observed for atmospheric pressure helium dielectric barrier discharges excited by kHz power sources [16,35]. It is observed that the electron temperature reaches around 3.9 eV as the maximum in the sheath region, while the density of N 4 + reaches around 5.2 × 10 16 m −3 . Both the species densities of electron and N 4 + reach the level of Figure 15 shows the cyclic average heating sources of different heating mechanisms. The simulated average total heating source is around 5.6 × 10 5 W m −3 with the maximum reaching 3.5 × 10 6 W m −3 in the sheath region due to the contribution of ion Joule heating. The total heating source reaches the level of 10 5 W m −3 in the bulk region due to the contribution of elastic collision. It is observed that the heating source of elastic collision increases to the maximum near the sheath region since the local electron density and electron temperature are higher than those in the central region. Similarly, the heating sources of ion Joule heating and exothermic reactions reach the maxima because of the high species densities generated and the electric field built in the sheath region.    Figure 15 shows the cyclic average heating sources of different heating mechanisms. The simulated average total heating source is around 5.6 × 10 5 W m −3 with the maximum reaching 3.5 × 10 6 W m −3 in the sheath region due to the contribution of ion Joule heating. The total heating source reaches the level of 10 5 W m −3 in the bulk region due to the contribution of elastic collision. It is observed that the heating source of elastic collision increases to the maximum near the sheath region since the local electron density and electron temperature are higher than those in the central region. Similarly, the heating sources of ion Joule heating and exothermic reactions reach the maxima because of the high species densities generated and the electric field built in the sheath region.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 21 density of N4 + reaches around 5.2 × 10 16 m −3 . Both the species densities of electron and N4 + reach the level of 10 16 m −3 in the central bulk region. Figure 15 shows the cyclic average heating sources of different heating mechanisms. The simulated average total heating source is around 5.6 × 10 5 W m −3 with the maximum reaching 3.5 × 10 6 W m −3 in the sheath region due to the contribution of ion Joule heating. The total heating source reaches the level of 10 5 W m −3 in the bulk region due to the contribution of elastic collision. It is observed that the heating source of elastic collision increases to the maximum near the sheath region since the local electron density and electron temperature are higher than those in the central region. Similarly, the heating sources of ion Joule heating and exothermic reactions reach the maxima because of the high species densities generated and the electric field built in the sheath region.

Conclusions
The present work studies the thermal behavior of APHeDBDs with a small amount of N 2 addition for temperature measurements by numerical simulations and experimental measurements.
The dynamics of APHeDBDs generated in a planar reactor is modeled by the 1D PFM with the discharge chemistry including He and N 2 related species to calculate the spatial cyclic average heating source considering elastic collision, ion Joule heating, and exothermic reactions. The average heating source obtained by the 1D PFM is adopted as the source term of the energy equation considered in the 3D GFM for modeling the temperature of the reactor. The 1D PFM and 3D GFM are coupled weakly to obtain the converged solution until the gas temperature used in the 1D PFM is consistent with that obtained in the 3D GFM.
The discharge uniformity is verified by the ICCD photograph to justify the application of the 1D PFM for modeling the uniform discharge. The simulated current density reaches 29 A m −2 which is close to that measured as 35 A m −2 . The simulated power consumption is around 2.0 W which is in good agreement with the average measured power consumption of 2.1 W. The simulated average gas temperature in the reactive zone is about 346 K which is also close to the rotational temperature determined from the emission spectrum of N 2 SPS. The high-temperature zone of the simulated surface temperature of the reactor moves slightly downward due to the convective effect, which is similar to the temperature distribution detected by the IR thermal imager. Moreover, the maximum simulated surface temperature is 362 K which is close to that measured as 365 K. In general, the numerical simulations capture the thermal behavior of the APHeDBDs.
The analysis shows that elastic collision and ion Joule heating are dominant heating mechanisms contributing 23.9% and 65.8% to the heating source, respectively. Detailed analysis shows that elastic collisions of electrons and helium atoms contribute more than 99.9% of the heating source of elastic collision since the concentration of He atoms is much higher than that of N 2 molecules. Among ion species, N 2 + and N 4 + are dominant species contribute 44.1% and 50.7% to the heating source of ion Joule heating, respectively. The simulated result shows that the Penning ionization of excited He species and N 2 molecules leads to the abundant generation of N 2 + while the fast associative process results in the dominant ion N 4 + . The simulated average total heating source is around 5.6 × 10 5 W m −3 with the maximum reaching 3.5 × 10 6 W m −3 in the sheath region due to the contribution of ion Joule heating.