Validation of the thermo-chemical approach to modelling of CO2 conversion in sub-atmospheric pressure microwave gas discharges

The CO 2→ CO + 12 O2 conversion experiment [D’Isa et al 2020 Plasma Sources Sci. Technol. 29 105009] has been compared with thermo-chemical calculations. The experiment is a 2.45 GHz plasma torch with a straight channel in the effluent. The 1.5D model of the CO2/CO/O2/O/C mixture without turbulent transport has been applied with plasma acting only as a prescribed heat source. The parameter range covered is the specific energy input (SEI) 0.3–5 eV/molecule at pressure p = 0.9 bar, and SEI = 0.6–2 eV/molecule at p = 0.5, 0.2 bar. The calculated conversion χ is always close to the experimental values. At the same time, the calculated temperatures T deviate significantly from the experiment, especially for p = 0.2 bar. The calculated T were also found to be sensitive with respect to the uncertain model parameters, but χ is not sensitive. The net conversion in the model is driven to a large extent by the radial diffusion of CO and O from the hot core towards the wall and steep radial temperature gradients. The main factor that reduces the energy efficiency is the re-oxidation of CO at the edge of the hot plasma region and downstream. The 1.5D approximation applied has the principle limitation that the impact of a realistic bulk flow field on the chemical process cannot be studied. Hence, the results must be considered preliminary and have to be confirmed with a more elaborate and accurate model of the vortex-stabilized flow inside the reactor.


Introduction
Plasma chemical conversion of CO 2 into CO has the potential to become an important part of the future Carbon Capture and Utilization process chains [1].Remarkable progress has been achieved in this field in recent years, in particular in microwave (MW) flowing plasma reactors where chemical energy efficiencies (efficiencies related to the energy absorbed in plasma) of the conversion process up to 40%-50% were achieved [2][3][4][5].However, further increase in efficiency and productivity-especially when operating at atmospheric pressure-is required to make this technology mature for practical applications.
It has been demonstrated experimentally that in this type of reactors installing nozzles and orifices, as well as using different gas supply schemes has a significant impact on the process performance [2,5,6].Thus, it may be possible to further optimize the conversion process by influencing the reacting flow characteristics.A purely experimental optimization would require building numerous prototypes and executing a large number of experiments, which may be prohibitive in terms of time and resource consumption.A common practice for reducing the number of prototypes and accelerating process optimization is extensive use of numerical modelling for planning and evaluation of experiments.
A full numerical description of the CO 2 conversion in an MW plasma reactor must couple self consistently gas flow with chemical reactions, plasma processes and electromagnetic fields.At least in 2D, ideally in 3D.Apparently, such a fully self-consistent modelling poses serious difficulties even for modern simulation tools running on modern computers.For MW discharges, 3D models have so far been reported only for Ar plasma [7,8].2D modelling studies of plasmas sustained by radio-frequency waves including chemical reactions have been done for N 2 [9,10].Recently, a fully selfconsistent 2D modelling of the CO 2 plasma chemical conversion in MW reactors has been published [11], however, only for CO 2 diluted in Ar in the proportion 1:7.Most relevant for practical applications is the conversion at atmospheric pressure which does not require vacuum equipment.Fortunately, it turns out that for atmospheric and near-atmospheric pressures practically relevant results can be achieved already by much simpler non-self-consistent models.
Particularly extensive experimental studies of CO 2 conversion in MW plasmas at pressures ranging from 1 mbar to 1 bar have been carried out in devices with 2.45 GHz TE 10 wavessurfaguides, plasma torches and similar configurations.Along with the experimental activities substantial effort has been invested in the past decade by the community in the development of detailed state-to-state models of the CO 2 kinetics and plasma chemistry [12].The availability of these models enabled thorough analysis of the plasma-chemical mechanisms in weakly ionized plasma of MW plasma reactors on microscopic level.This analysis led to the conclusion that in the medium to atmospheric pressure range (⩾200 mbar) reactions with charged particles play a minor role in the net chemical process [13][14][15].The observed conversion can be largely explained by chemical non-equilibrium where the plasma only plays the role of the heat source [2,3,16,17].Moreover, in the pressure range in question the plasma is found to be contracted near the axis and the localization of this heat source in experiment can be reconstructed by applying the impedance matching technique [18] and from the measured radiation intensity [3,4,[19][20][21][22].
The observations summarized above open the possibility to use a purely thermo-chemical approach as a first stage of optimization of the nozzle and gas management configurations.The term 'thermo-chemical' is applied here to models where only the gas flow with chemical reactions is simulated.The plasma processes are not modelled explicitly and the plasma only enters as a prescribed heat source.In such models equations for the mass, energy and momentum (Navier-Stokes equation) conservation of the bulk flow are solved together with drift-diffusion equations with finite rate chemistry for each chemical component.The thermo-chemical models are well known and widely used in chemical engineering and in studies of combustion processes [23].The task of describing the chemical reactions is further simplified by the fact that in MW plasmas with relevant parameters no significant nonequilibrium between vibrational and translational-rotational temperatures has been detected [14,24].Thus, one can rely on readily available data on the reaction rate coefficients determined for conditions of temperature equilibrium.
The main goal of the present paper is to validate how the thermo-chemical approach works in principle by benchmarking it against the plasma torch experiment of D'Isa et al [4,22] with a straight channel in the effluent (without a nozzle).Emphasis is made on reproducing trends: dependencies of the target quantities on input power and flow rate as well as on pressure, rather than trying to match as exactly as possible each individual experiment.One method of the thermo-chemical analysis which can be found in the literature is modelling of the 2D or 3D single fluid flow without chemical reactions, and then applying 1D chemical kinetics models along flux tubes [25,26].In this case only the calculated temperature can be compared with experiment, and the local 1D analysis gives only qualitative insight into the interplay of the bulk flow and the chemical process.The focus of the present study is to calculate the final net conversion and to compare it with experiment.Therefore, another concept is used which goes back to Wolf et al [17].It is based on applying a simplified model for the bulk flow in a cylindrical discharge tube, but solved self-consistently with a set of transport equations for the individual components of the CO 2 /CO/O 2 /O/C mixture including chemical reactions.In this method the flow field is described only schematically, but it allows to calculate the conversion and energy efficiency with moderate computational overhead and compare them directly with experiment.
The transport model applied here can be called a '1.5D approximation'.Equations for heat balance and drift-diffusion equations are solved in 2D assuming axial symmetry.To calculate the bulk flow velocity a local 1D approximation is used instead of the full Navier-Stokes equation.This approximation uses the assumptions of constant static pressure and prescribes local mass conservation in axial flux tubes neglecting radial mass transfer.The tangential swirl which is always applied in the discharges in question for stabilization, is also not accounted for explicitly in this model.Whereas the constant pressure assumption is well justified for Mach numbers ≪1, the two other assumptions have no formal justification, but they allow great simplification of the calculations.At the same time, these are the principal limitations and shortcomings of the 1.5D approximation.Fluid dynamics simulations [6,27,28] have shown that in vortex-stabilized MW reactors complicated bulk flow patterns can build with reversed flows and convective cells.The influence of these flow patterns on the net chemical process cannot be studied within the 1.5D approximation where the axial flow is only a schematic 'pusher' which moves the reaction products and heat forwards.It should be stressed that especially because of this drawback of the approximation applied the results of the present study must be considered preliminary and have to be confirmed with more physically accurate models, which take into account realistic flow patterns.In the present study only a preliminary evaluation of the effect of the flow profile on the solution is done by performing numerical tests with zero axial bulk velocity near the axis-hollow flow profile.These tests show that in the straight cylindrical channel case considered here the variation of the flow field shape does modify the calculated convergence, but this modification is not large.
In the present paper only molecular diffusion and heat conduction coefficients are applied without taking into account the turbulent transport.There are two rationales for that.First, at high temperatures considered (up to 6000 K) the purely molecular transport is already strong.Second, the turbulence effects can significantly enhance the transport, but they are always added on top of the molecular transport that is always present in the system irrespective of the flow pattern and the level of fluctuation.Thus, the decision was taken to start with pure molecular transport as a first step to identify to what level of disagreement with experiment this will lead before addressing the complicated topic of turbulent transport theory.This is the main difference of the method of the present study compared to [17] where the turbulent viscosity is used as an adjustable parameter to match the measured peak temperature.Here it was not tried to enforce the solution to match the experimental temperature.Instead, the sensitivity with respect to the assumed distribution of the heat source is evaluated.This sensitivity evaluation is mandatory because the procedure for determining the localization of the plasma zone in experiment is known to be inaccurate.
Since in the simulations of the present study no specific value of the peak temperature is enforced by the model this parameter can also be compared with the measurements.In contrast to the conversion the agreement for the temperature was found to be largely bad, especially at the lowest investigated pressure 0.2 bar.Moreover, whereas the calculated conversion was found to be insensitive with respect to the parameters of the spatial distribution of the heat source, the calculated temperature was found to be very sensitive to these uncertain model parameters.The apparent decoupling of the peak temperature and the conversion is a reflection of the conversion mechanism as it appears in the model.The net transformation of CO 2 into CO in the model is to a large extent not driven by the thermal quenching [29,30] along the stream lines.An important driving process appears to be the radial diffusion at the edge of the hot plasma core along the steep temperature gradients that form there.The latter is the reason for the chemical non-equilibrium and partial mitigation of the backward recombination of CO.This conversion mechanism resembles the concept of the 'hot core' and 'cold shell' proposed for CO 2 conversion in [16,17] and adopted in [31] to CH 4 plasma reforming.
The rest of the paper is organized as follows.In section 2, the model is described, which includes the set of equations solved, the chemical mechanism and the transport coefficients.In section 3, a comparison of the model with the experiment [4,22] is discussed.Section 4.1 presents the assessment of sensitivity with respect to the uncertain model parameters.Section 4.2 gives some insight into the near-atmospheric pressure conversion process provided by the model emphasizing the role of radial transport.Section 5 summarizes the main conclusions.

Description of the model
It is convenient to divide the model into two blocks: the first is the transport equations and the second is the calculation of the coefficients of those equations-the 'physicochemical model'.This division reflects the software-technical implementation.

Transport equations
The set of transport equations applied here consists of the continuity equations for each chemical component, an energy balance equation and a schematic local 1D model for the axial velocity.Equations are solved assuming axial symmetry in cylindrical coordinates (z, r), where z is the coordinate along the axis and r is the distance from the axis.Only steady-state solutions are considered.
The continuity equations-convection-diffusion equations-read as follows: where ⃗ v is the average mass flow velocity, n i is the number density of the component i, S i is the specific volumetric particle source of the component i, n = ∑ i n i .Here and below ∑ i always stands for the sum taken over all chemical components.In (1) the diffusive flux is approximated by Fick's law with the effective mixture-averaged diffusion coefficient D im , rather than calculated by solving the full set of Stefan-Maxwell equations (see [23], chapter 3.5 there and [32]).The thermal diffusion and the pressure gradient (see [23], section 3.5.2) are neglected.The latter would contain the centrifugal effect of the swirl flow on the diffusion.
The following equation is solved for the gas temperature T-which is the common temperature of translational, rotational and vibrational modes, which are assumed to be in equilibrium: This equation is a reduced form of the total energy conservation equation for multi-species flow ( [33], section 18.3).In (2), λ is the effective heat conductivity of the mixture, Q is the specific volumetric heat source due to the plasma heating and h i is the variable part of the specific enthalpy per particle of the component i.The h i is defined so that h i (T ref ) = 0, where T ref is some reference temperature chosen to be the same for all components and the full enthalpy is defined as The constants H i (T ref ), which include the enthalpy of formation, are formally moved to the 'chemical' heat source S h calculated with the help of (1) as, To solve (2) numerically it is translated into the standard conduction-convection form by expressing h i as h i (T) = c i h (T) • T. Radiation heat losses are not taken into account in (2).Since only slow flows with Mach numbers M ≪1 are considered, kinetic energy of the bulk flow 1  2 ρv 2 is omitted in (3).The work of viscous stress is also omitted.Here, ρ = ∑ i m i n i , m i is the molecular weight of the component i.
The schematic model applied to find the bulk flow velocity ⃗ v is based on two assumptions.The first is the assumption of constant pressure p = const, which is justified, again, by small M. In all calculations that will be discussed below the peak local Mach numbers of the flow are always <0.1.The second is the assumption of zero bulk radial velocity, v r = 0, that is, the mass transfer in the radial direction is neglected.The latter assumption has no formal theoretical justification, but it allows to greatly simplify the calculations.In a channel with constant cross-section this condition implies local mass conservation.Thus, the axial velocity v z is readily calculated as, where Γ ρ (r) is determined by the boundary condition at the inlet.The total number density n is calculated from the equation of state, section 2.2.3.Equation (4) yields increased velocity at the locations with increased T, which is qualitatively correct.Sensitivity of the solution with respect to the flow details which go beyond this simple property can be evaluated to some extent by varying Γ ρ (r).
The set of equations ( 1), ( 2) and ( 4) is solved numerically using a self-written finite-volume code.Some more details on the implementation and numerics are given in appendix A.

Physico-chemical model
The physico-chemical model as it is defined here comprises: (a) Chemical mechanism: stoichiometric equations and reaction rate coefficients required to calculate the source terms S i in (1) (b) Transport coefficients D im and λ, which enter ( 1) and ( 2) (c) Thermodynamic data Five components are taken into account: CO 2 , CO, O 2 , O, C. The assumed chemical mechanism is summarized in table 1.The mechanism consisting of the processes N1-N3 goes back to [34]; this same mechanism was applied in [3,16].The processes N4 and N5 were added following the suggestion of [17] and because our own experience indicated that inclusion of CO dissociation is mandatory to obtain solutions with realistic temperatures.
The full set of chemical reactions included in the model is shown in table 2. The rate coefficients are calculated by applying the generalized Arrhenius equation: The coefficients A, β and E a are given in table 2.An evaluation of the different sets of rate coefficients in a microwave CO 2 plasma conversion experiment was performed in [14] on the basis of spatially resolved Raman measurements of concentrations and comparison of the reconstructed and experimental conversion and input power.The conclusion of [14] was that the best results are obtained with the rate coefficients taken from the chemical mechanism GRI-Mech 3.0 [35].This is a comprehensive chemical mechanism developed for modelling methane combustion.Therefore, most of the rate coefficients used in the present study are taken from that data set.For the process N1 only the low pressure limit is taken.The rate coefficients of the process N3 with M = O, O 2 and of the process N4, which are missing in [35], are taken from the data evaluations [36,37].
The equilibrium factors K eq , which appear in table 2, are calculated using the formula known from thermodynamics (see e.g.[23], equation (9.43)): where ν i are the stoichiometric factors of each chemical component in the corresponding process, by definition ν j > 0 for

Stoichiometric equation
A, β and Ea are the parameters of (5); Ea is in Kelvin; the units of A are chosen so that when T is in Kelvin the R (T) in ( 5) is in m 3 s −1 for second-order reactions, and in m 6 s −1 for third-order reactions.
For calculation of the equilibrium constants Keq see (6).
products and ν j < 0 for reactants, the sums are taken over all components that take part in the given reaction; G ref j (T) are Gibbs energies of the components at the reference pressure p ref , k B is the Boltzmann constant.

Transport coefficients.
The mixture-averaged diffusion coefficient D im in ( 1) is calculated by means of the expression derived in [32]: Here, D ij (T, p) is the binary diffusion coefficient calculated for the total pressure of the mixture p.According to [32], (7) gives correct approximation when i is a trace impurity whose concentration is small.One can also easily show that this formula gives exact results for binary mixtures.Most of the coefficients D ij used in the present study are taken from [38], equation (4.3-2), table 13.For the pairs missing in [38] Fuller scaling is used: [39], section 11-4, equation 11-4.4 and table 11-1.
The thermal conductivity of the mixture is calculated by applying the Wassiljewa equation (see [39], equation (10-6.1)): where λ i (T) are the thermal conductivities of individual components and ϕ ij are defined as ('Mason and Saxena Modification', equation (10-6.4) in [39]): where µ i is the dynamic viscosity.In the calculations, a simplified version of this formula is used obtained on the assumption that for the dependency of µ on the molecular weight m the scaling µ ∼ √ m holds.This scaling follows from the known theoretical expression for the viscosity derived for spherical molecules without internal degrees of freedom ( [40], equation (8.2-10)).Then, ( 9) is reduced to: It is worth noting that there was no particular reason to use specifically equations ( 8) and (10) for the calculation of λ.It is not unlikely that other simpler expressions which can be calculated faster, such as the formula suggested in [23], equation (12.119) there, will yield equally good results.To calculate the coefficients λ i for CO 2 , CO and O 2 the correlation formulas are used obtained by combining experimental data for low temperatures with theoretical extrapolations to higher temperatures.The λ i for CO 2 is taken from [41], equation (3), table 3. CO and N 2 are known to have very close thermal conductivities [42].Therefore, as λ CO the thermal conductivity of N 2 is used taken from [43], equations ( 5) and ( 2), tables I, II and IV.The λ O2 is taken from the same reference [43].Thermal conductivity of O atoms is a fit to the results of calculations reported in [44], table VI.For C atoms λ i is calculated from the diffusion coefficient defined by Fuller scaling by applying the theoretical formula that connects λ and the coefficient of self diffusion, see appendix B.
More detailed information on the evaluation and crosscheck of the literature data on D ij and λ i is given in supplemental material.6) are all taken from [45]; T ref = 0 K, p ref = 10 5 Pa.The equation of state is that of ideal gas p = nk B T.

Description of the experiment and its computational model
The experimental set-up will be only very briefly described here since its detailed description can be found in [22] and in the open access publication D'Isa et al 2020 [4].The discharge is ignited in a quartz tube that crosses a resonator where a 2.45 GHz wave is excited by a magnetron generator, see figure 1(a).The resonator design of Stuttgart University is used which consists of the main cylindrical resonator with E 010 wave (E-field is directed along the quartz tube axis) for continuous operation and the additional coaxial resonator serving for ignition of plasma [46,47].The working gas is injected through four tangential gas inlets.The discharge tube is followed by a heat exchanger connected to a vacuum pump.The conversion is obtained by measuring the mixture composition with a mass spectrometer [48].Optical access is available from the window on the top and through the side slit of the cylindrical resonator.Optical emission spectroscopy is applied to measure the spatial distribution of the plasma radiation as well as for the reconstruction of the gas temperature [24].
The computational model simplifies the real geometry by a straight cylindrical channel, figure 1(b).In particular, the tip of the coaxial resonator at the bottom is neglected.The origin of the coordinate system (z, r) corresponds to the centre of the heat source where it assumes its peak value, equation (11) below.This point is located on the discharge tube axis close to the middle of the cylindrical resonator.The internal diameter of the channel d = 26 mm, the length l from the gas inlet to the origin is set to l = 2l pl , where l pl is the estimated plasma length (see below).The length L from the origin to the end of the computational domain is set to L =1 m (in some cases to L =1.5 m).This choice ensures that at that location in all the simulation runs considered the gas temperature is guaranteed to be below 1500 K, the molar fraction of oxygen atoms is always <2•10 −5 , and to this point all the chemical transformations have been completed.
The boundary conditions of equations ( 1) and ( 2) are described in table 3. The influx Γ n in the reference model is set at the total flow rate divided by π d 2 /4.That is, it prescribes a uniform spatial distribution of the CO 2 flux through the inlet.Γ ρ in ( 4) is Γ n multiplied by the molecular weight of CO 2 .Single-fluid CFD calculations made for the experimental set-up of figure 1(a) and similar experiments [6,28,49] indicate that due to swirl imposed around the resonator tip hollow profiles of axial velocity may form or even convective cells with reversed flow.That is, a profile where the axial velocity near the axis is smaller than near the walls or is even in the opposite direction.The impact of convective cells on the net chemical effect cannot be directly tested within the applied 1.5D approximation.To roughly estimate the sensitivity of the solution with respect to the axial flow pattern a case was tested with Γ n (r) set to zero for r from 0-8 mm, which is the radius of the resonator tip.The corresponding numerical experiments are discussed in section 4.2 below.
In the experiment the quartz tube is actively cooled by pressurized air to keep its temperature low.There are signs that the temperature of the quartz tube wall T w is higher than the room temperature, but its exact values were not measured.Therefore, in the reference model T w = 300 K is specified everywhere on the wall on the implicit assumption that the solution must not be sensitive with respect to that boundary condition.This assumption was checked in a number of tests with elevated T w , see section 4.1 below (table 5).
The spatial shape of the power source is specified on the basis of the measured distribution of the plasma radiation intensity.The basic idea behind this method is that since the power is introduced into the gas via free electrons the specific volumetric input power should be approximately proportional to the electron density n e .At the same time, the line radiation intensity is proportional to n e as well.Therefore, the profile of the radiation intensity should roughly reflect the profile of the specific input power.A more thorough analysis [19,21] has confirmed that the profile of the 777 nm O-atom line radiation is indeed representative-to a certain degree of accuracy-of at least the radial power deposition profile.In D'Isa et al [4,22] this particular line was not used, but the total radiation, which is in contracted discharge mode is shown to be mainly due to the C 2 bands.Nevertheless, extra measurements done with the Table 3. Boundary conditions for the computational domain of figure 1(b).

Boundary
equation ( 1) equation ( 2) Γn is the inlet flux density, see section 3.1 filter confirmed that the distribution of the total radiation is the same as that of the oxygen only [4].
The measured radial radiation profile has approximately Gaussian shape.Therefore, in the calculations the spatial distribution of the specific input power is approximated by the following function: where the constant c is adjusted so that the integral of Q (z, r) over the computational volume is equal to the prescribed total input power.Whereas the use of a Gaussian function for the radial profile is well justified, this is less certain for the axial profile.Its use in the present paper goes back to the work [18] where Q (z, r) was reconstructed by applying the impedance matching technique and both z and r profiles were fitted by Gaussian functions.However, there are indications that for the experiments of D'Isa et al [4,22] a shape less peaked in the middle could be more appropriate.Therefore, while in the reference calculations (11) is always applied, several cases were repeated with flat z-profile to demonstrate that the final solution is not strongly affected by this assumption, see section 4.1.Parameters α r , α z in (11) are reconstructed from the 'plasma diameter' d pl and 'plasma length' l pl reported in [4], figures 9 and 11 there respectively.The plasma radius d pl /2 is defined in [4,22] as the radius where the radiation intensity is 15% of its maximum.Hence, α r is calculated by applying the following rule: exp The same rule is used to calculate α z with d pl replaced by l pl .It is known that the procedure of determining the characteristic dimension of the power deposition profile from the plasma optical size is not accurate.A specific example of the error that may be introduced by this method can be found in [21].There it is shown that for pressures <150 mbar the power deposition radius is by a factor 1.6 larger than the optical plasma radius.Therefore, evaluating the effect of increasing or reducing d pl and l pl is always part of the computational analysis when the experimental Q (z, r) is used.The impact of this variation will be discussed in section 4.1 below.

Comparing the reference model with the experiment
The results of a comparison of the calculations performed with the reference model described in section 3.1 above and the experiment of D'Isa et al [4,22] at quasi-atmospheric pressure of 0.9 bar are presented in figure 2. The peak temperature  [4,22] for pressure 0.9 bar and different flow rates.SEI is defined by (13), 'conversion' χ is defined by ( 14) and 'energy efficiency' η by (15).Dashed lines are obtained with the nominal value of d pl increased by a factor 1.5 and the dotted lines with d pl decreased by a factor 1.5, see section 4.1.
T max , conversion χ and energy efficiency η are plotted as functions of the specific energy input per molecule (SEI).The latter is defined as, where 'input power' is the power absorbed in the plasma.Measurements of the rotational temperatures [24] have shown that in the contracted mode, in a wide range of experimental parameters, T max always acquires similar values.Therefore, for all experiments, this value is set here at T max = 6000±500 K [4,24].The conversion is defined as, where 'CO outflux' in the model is the total flux of CO molecules through the outlet (see figure 1(b)).The experimental χ was determined by the mass spectrometer measurements after the heat exchange.The measured molar flows are translated into χ taking into account the increase in the number of particles after chemical transformation, as described in section VI.B in [48].The error bar of the measurements is set to ±1.6% according to [48].The energy efficiency is connected to χ by the equation: where ∆H f = 2.93 eV/molecule is the net enthalpy change in the chemical transformation CO 2 →CO+ 1 2 O 2 at T = 298.15K [45].
The temperature T max , figure 2(a), is not very well reproduced by the model.In most cases the (reference) model tends to overestimate T max by 20%-30%.However, in the case with the highest flow rate 40 slm the temperature may even be underestimated.Despite that, as one can see in figures 2(b) and (c), the χ and η are in good agreement with the experiment.The conversion is mostly close to the experimental error bars, and the trend observed in the experiment is reproduced, except for the largest SEI (lowest flow rate).The reduction of χ observed when the SEI is increased above 3 eV/molecule is not reproduced by the model which yields saturation of χ in this SEI range.In contrast, for η the largest deviation can be seen at SEI<1 eV/molecule where the experimental data indicate near saturation with decreased SEI, whereas in the model η still increases, although the values are formally within experimental error bars.[4,22] for pressure 0.5 bar and flow rate 20 slm.SEI is defined by ( 13), 'conversion' χ is defined by ( 14) and 'energy efficiency' η by (15).Dashed lines are obtained with the nominal value of d pl increased by a factor of 1.5 and the dotted lines with d pl decreased by a factor of 1.5, see section 4.1.

4.
Comparison of the model with the experiments [4,22] for pressure 0.2 bar and flow rate 20 slm.SEI is defined by ( 13), 'conversion' χ is defined by ( 14) and 'energy efficiency' η by (15).Dashed lines are obtained with the nominal value of d pl increased by a factor 1.5 and the dotted lines with d pl decreased by a factor 1.5, see section 4.1.
The results for pressure 0.5 bar are similar to those for 0.9 bar, figure 3. The situation is different for the lowest investigated pressure 0.2 bar, figure 4. In this case the model overestimates the temperature by up to a factor two, figure 4(a), to the extent that the modelling results are physically inconsistent.At temperatures above 8000 K, even in the state of local thermodynamic equilibrium, the plasma cannot be weakly ionized anymore [50], which contradicts the model assumptions.Surprisingly, despite this drastic mismatch the conversion and energy efficiency are still very well reproduced by the model, figures 4(b) and (c).The main disagreement is that in the experiment η is decreased with increased SEI, and in the model this trend goes in the opposite direction, although in both cases the trend is weak.
The calculated profiles of the gas temperature are compared with measurements in figure 5 which corresponds to the exact same experimental case as figure 6 in [4].The error bars of the measured temperature are set to ±500 K and the error bars of the spatial coordinates in the measurements are set to ±0.5 mm according to [4].The radial profile corresponds to the location '58 mm above the resonator bottom' [4].To translate this position into the model coordinates (shown in figure 1(b)) it was assumed that the centre of the heat source is located 30 mm above the resonator bottom, as can be estimated approximately from figure 10(d) in [4].One can see in figure 5 that not only the absolute values of temperatures are different, but also that both radial and axial profiles in the model are peaked whereas in the experiment they are nearly flat.

Sensitivity assessment
The sensitivity of the solution with respect to the shape of the heat source Q (z, r) has been investigated for six selected cases.They are listed in table 4 which presents the variation of the maximum temperature T max .In this table, the column 'reference' is the solution obtained with the reference model of section 3.1.The next two columns show the results of using for Q (z, r) other types of spatial shape than function (11).The column 'flat Q(z)' stands for the heat source described by the function: otherwise.
The column 'rect.Q' stands for the flat profile in both axial and radial directions.That is, for the 'rectangular shaped' source: Spatial distribution of the gas temperature in the case of 0.9 bar, 0.9 kW and 10 slm (same case as in figure 6 in [4]).(a) Radial profile at z =28 mm, which corresponds to '58 mm above the resonator bottom' in [4];  Other columns show the effect of increasing and decreasing the assumed plasma diameter d pl and plasma length l pl by a factor 1.5 compared to their values taken from [4,22].Whereas the effect of switching to the flat Q (z)-profile or rectangular Q (z, r) on T max is visible but limited, the effect of the variation of d pl , l pl is large.Its magnitude is further illustrated in figures 2(a)-4(a) where the solutions obtained with increased and reduced d pl are plotted as dashed and dotted lines, respectively.At lower pressures, the effect of the variation of d pl on T max is notably stronger than at the quasiatmospheric pressure.Figure 5 presents an example of the modification of the temperature profiles calculated with different model assumptions.Variations of d pl and l pl do not greatly affect the peaking of the radial and, subsequently, axial temperature profiles.At the same time, replacing the Gaussian Q (z)-profile with the flat profile leads to flattening of the calculated axial temperature profile as well, bringing its shape closer to that in the experiment.This result is an indication that a relatively flat Q (z) distribution could indeed be a more realistic assumption than the Gaussian profile applied in the reference model.
As one can see in figures 2(a)-4(a), the assumption that the real d pl is larger than its optical appearance can bring the calculated T max closer to the experimental values.This assumption would be in line with the observation made in [21].However, the 0.9 bar case with the lowest SEI (flow rate 40 slm) where the temperature calculated with increased d pl is too low would not fit into this hypothesis.Another factor that can strongly affect the calculated T max is the ionization, which is missing in the model.This factor could be, in particular, responsible for too high temperatures at the lowest pressure 0.2 bar.Because of the known limitations of the applied model and because from the beginning that was not the intention of the present study, it was not tried to fit the calculated temperature exactly into the measurements by adjusting the heat source shape.
Despite the large variation of T max in all cases listed in table 4 the conversions were found to vary only within ±0.3% at maximum, staying in most cases unchanged within three decimal digits.This result can also be seen in Figures 2(b)-4(b) and 2(c)-4(c) where the effect of increasing and decreasing d pl by a factor 1.5 is displayed by dashed and dotted lines, respectively.
The very weak influence of the large variation of the maximal temperature in the plasma core on the final net conversion χ has the following qualitative explanation.At 5000 K most of the CO 2 is dissociated and the temperatures of that level and above in the hot core only affect the degree of dissociation of CO into C and O there.This dissociation itself cannot change the net CO production effect because eventually the C always recombines back into CO.As will be discussed in section 4.2 below the net conversion here is driven mainly by the radial temperature gradient at the edge of the hot core where the temperature drops down to 4000 K and below.This gradient apparently does not change much with different Q (z, r) shapes; an example can be seen in figure 5(a).Therefore, the fact that χ is not sensitive with respect to the heat source distribution is within expectations.Nevertheless, the degree of this insensitivity seen in the simulations is surprising and so far has not found good theoretical explanation.
For the moment this result has been accepted solely as the outcome of numerical experiments.The observed insensitivity of the calculated conversion is in general a positive result which offers the possibility of relatively accurate calculation of χ (and η) despite large uncertainty in the assumed heat source shape.
Another model parameter that can be considered as uncertain is the wall temperature T w set as the boundary condition.In the reference model, T w = 300 K is assumed as the first 'educated guess'.Since no exact measurements of this temperature are available, extra tests were performed to check the sensitivity of the solution with respect to the increase in T w .The same cases as in table 4 were repeated with T w = 600 K and T w = 1000 K.At temperatures above 700 • C quartz glass starts to glow visibly.In the experiments this glowing was not observed, which means that the temperature was in any case ⩽700 • C. Therefore, the assumption T w = 1000 K represents the absolute upper boundary.Especially taking into account that as an extreme case this boundary condition is applied over the whole length of the channel.The results of the tests are presented in table 5.One can see that T w = 600 K does change the calculated conversion χ only in the third decimal digit.T w = 1000 K yields visible but small modification for 0.2 bar and for 0.9 bar with low SEI.Only in the two cases with SEI⩾2.5 eV/molecule is the modification substantial, although still quantitatively small.The effect on the maximum temperature T max is negligible in all cases.

Numerical insights into the conversion process
Since the agreement between the calculations and the experiment for 0.9 bar is relatively good, the model be used to gain some insight into the mechanism of the conversion process.In particular, insight into the spatial and transport effects.In the first approximation, the process can be separated in space into two regions.Into the 'plasma region' where the gas is heated and into the 'effluent region' downstream where the heating stops and the gas mixture cools down.
The composition of the hot mixture is illustrated in figure 6, which shows the calculated radial profiles of the molar fractions in the middle of the heat source (z = 0) and near the boundary between the 'plasma' and the 'effluent' (z = 28 mm).This latter is the same cross-section as that in figure 5(a).The example is given which corresponds to the same case as in figure 5 and to the case with twice as large input power.As expected, the gas in the centre is fully dissociated and consists only of CO/O/C.The profiles in the higherpower case are broader because of the larger assumed plasma diameter, d pl = 3.6 mm for 0.9 kW and d pl = 5.5 mm for 1.8 kW.One can see that, compared to the z = 0, in the cross-section at the end of the plasma region the molar fractions of CO 2 and O 2 at the centre and of CO and O at the edge are increased.This increase is the result of radial diffusion.
The significance of the diffusion process is revealed quantitatively by examining the particle balance in each individual flux tube.This examination showed that for all chemical components the contribution of the effective source due to the radial part of div( ⃗ Γ i ), see (1), is always comparable to that of the volumetric source S i .The importance of the radial transport is further illustrated in figures 6(c) and (f), which show the total radial fluxes of each component-fluxes integrated over cylindrical surfaces that start at the inlet and end at the outlet.Since in the present model there is no bulk radial convection, all these fluxes are purely diffusive.To make the quantitative role of the radial diffusion apparent, the fluxes are divided (normalized) by the flow rate.The normalized fluxes can be directly compared with the conversion: χ = 10% in 0.9 kW case and χ = 12.4% in 1.8 kW cases.The values of χ are similar in magnitude not only to the total fluxes shown in figures 6(c) and (f) by solid lines, but also to the partial fluxes that represent the plasma region only shown by dashed lines.Positive values in the flux diagrams correspond to the fluxes from the axis.As the the concentration gradients prescribe (figures 6(a), (b), (d) and (e)) CO and O diffuse from the axis towards the edge, and CO 2 diffuses to the centre.Note that the non-zero flux of CO 2 towards the centre at r < 2 mm, despite nearly zero molar fraction of CO 2 there, as seen in figures 6(a) Radial profiles of the molar fractions.(a), (d) is the middle of the heat source; (b), (e) is the same cross-section as in figure 5(a).(c), (f) are the total radial fluxes of the individual components-integrated from inlet to outlet, dashed lines-integrated only from inlet to z =28 mm.Case 0.9 bar, 10 slm, two different values of input power.and (d), reflects that the CO 2 content in the hot core upstream (z < 0) is yet higher than it becomes downstream (z > 0).
The CO produced in the hot core partly recombines back into CO 2 at the edge of the plasma region and in the effluent region due to the reverse reactions of the processes N1 and N2 (see table 1).This backward transformation determines the final net effect of the conversion process.The impact of the effluent on the overall energy efficiency can be illustrated by introducing the 'cumulative local efficiency' η (z) defined by modifying ( 14) and ( 15) as follows: Examples of the η (z) profiles for selected cases with different SEI are shown in figure 7. The maximum of η (z) is reached near the end of the plasma region, and its absolute value is almost the same in all three cases-around 40%.Note that the assumed plasma length in the last case l pl = 81 mm is larger than that in the two previous cases (l pl = 55 mm).Therefore, the position of the maximum is shifted to the right.The reduction of η (z) after reaching the peak value gets stronger with increased SEI, which produces the net effect seen in figure 2(c).The increase in the total axial CO 2 flux between the end of the plasma region and the outlet obtained in the model agrees qualitatively with the experimental results reported in [51].In [51], a similar microwave discharge was investigated, and the conversion was measured at two locations: near the plasma (downstream) and at the outlet.It was shown that the conversion determined near the plasma region is indeed much larger than the final value at the outlet.
The 1D representation of figure 7 can be further refined by considering two zones.The 'dissociation zone' where the net local volumetric source of CO 2 is negative, S CO2 (z, r) <0, and the 're-oxidation' zone where S CO2 (z, r) >0 3 .An example of the spatial distribution of S CO2 in both zones is shown in figure 8. Individual contributions of each zone can be added to the diagrams of figure 7 by introducing the following factors: It can readily be seen that due to particle conservation η (z) = η diss (z) − η reox (z).Figures 8(b) and (d) suggest that substantial re-oxidation can take place already at the edge of the plasma region, and not only in the effluent.The profiles of η reox (z) plotted as dashed lines in figure 7 confirm that the CO losses there lead to a substantial reduction of the peak η (z).This latter, if determined solely by η diss (z), could have reached almost 60%  13), and ηreox (z): equations ( 17) and ( 18), respectively.Vertical dashed line is the cross-section z =28 mm, the same as in figure 5(a).at maximum.Nevertheless, one can also clearly see that the major negative impact of the re-oxidation which is responsible for the drop in η with increased SEI is localized downstream from the plasma region in the effluent.
Essential process that eventually leads to non-zero net conversion here is the radial diffusion.This proposition can first be illustrated by the following simple 'back of envelope' estimate.The radial temperature gradient at the boundary of the hot core starting from the location where the temperature drops down to 3000-4000 K has the order of magnitude ∆T ∆r ∼ 1000 K 1 mm = 10 6 K m −1 , see e.g.figure 5(a).The diffusion coefficient at atmospheric pressure and at T ⩾3000 K has the magnitude D ≳10 −3 m 2 s −1 , see figure 1 in supplemental material, and the characteristic radial diffusion length of O and CO is 6).The estimated D and ∆r d translate into effective radial velocity of O atoms and CO molecules which diffuse from the hot core outwards u r ≈ D ∆r d ∼1 m s −1 .In turn, combining u r and ∆T ∆r yields the estimate of the effective cooling rate in the moving frame: ∆T ∆r • u r ∼10 6 K s −1 .A more elaborate example is given in figure 9, which shows the effective cooling rate in the frame moving with O-atoms extracted from a simulation run.One can see that u r is indeed around 1 m s −1 , and the cooling rate may locally even exceed 10 6 K s −1 .
From the basic 0D theories of the CO 2 thermal quenching [29,30] it is known that sufficiently fast cooling of a hot CO/O mixture leads to preferential recombination of O atoms into O 2 rather than recombination via the process CO + O + M → CO 2 + M. According to [29,30], for the  quenching process to reach the maximum efficiency cooling rates of 10 8 K s −1 and larger are required.However, a significant reduction in the CO re-oxidation and the net CO 2 conversion can be achieved starting already from the cooling rate of 10 6 K s −1 : see [30], section 5.7.2.The 'cooling rate' in figure 9(c) does not have exactly the same meaning as that in the basic 0D theories.In the diffusion process considered here different components have different velocities-u r of CO is smaller than that of O, and the interaction is complicated by the 'fresh' CO 2 moving in the opposite direction.Nevertheless, one can see that the results of [29,30] and the present calculations roughly coincide.
The consideration above provides the following interpretation of the conversion mechanism in the system in question.
The key non-equilibrium process which leads to the net ical effect is the radial diffusion.The combination of relatively large diffusion velocity and steep temperature gradient leads to effective cooling rates which favour recombination of O atoms into O 2 and in this way partially mitigate the backward re-oxidation of CO.This mitigation is far from being 100% efficient, but is sufficient to obtain a non-zero probability for CO from the hot core to penetrate into the cold surrounding gas.The 'non-ideality' of the process manifests itself in the formation of the re-oxidation zone as visualized in figure 8.
The idea that the CO 2 conversion in microwave discharges is pushed mainly by the transport processes in the radial direction rather than by thermal quenching in the flow direction goes back to the 'hot core/cold shell' model put forward in [16] and further developed in [17].The fact that radial transport alone can explain the conversion in the cases considered here is illustrated in the following numerical experiment.The axial velocity of the flow which passes through the hot core is set to zero by setting the influx Γ n (r) = 0 for r from 0 to 8 mm (see section 3.1).The absolute value of Γ n (r) is scaled to meet the specified total flow rate.The 8 mm radius is the radius of the resonator tip, figure 1(a), and the velocity profile specified in this way would correspond to the hollow profile forming due to the swirl imposed around the tip, see [49].Tests have been performed for the four 0.9 bar cases listed in table 4 and the results are presented in table 6.
One can see that the calculated conversion χ is not reduced but is even larger than in the reference model with uniform Γ n (r).In contrast to the effect of modifying the heat source distribution Q (z, r) discussed in section 4.1, the impact on the magnitude of χ is significant, but quantitatively is limited to an increase by 30% at maximum from the reference value.Comparison with figure 2(b) indicates that a more realistic flow model could improve agreement in the medium SEI range.Table 6 also shows that T max is sensitive with respect to Γ n (r) as well as to Q (z, r).With the same Q (z, r) shape the hollow Γ n (r) leads to increased T max .

Conclusions
The main idea of the thermo-chemical approach for modelling the CO 2 conversion in microwave flow reactors is that the plasma only enters the model as a prescribed fixed heat source whose spatial localization is taken from experiment.The mathematical problem is then reduced to modelling of a nonisothermal flow with chemical reactions.In the present paper this approach has been benchmarked against experimental data of D'Isa et al [4,22] by comparing them with 1.5D numerical transport calculations.The 1.5D transport model solves the 2D equations for the heat balance and for the particle balance of each chemical component and applies a schematic local 1D approximation for the axial flow.A basic physicochemical model of the CO 2 /CO/O 2 /O/C mixture has been proposed which comprises the minimal reaction mechanism and molecular transport coefficients without turbulent transport.A comparison has been performed for pressures 0.9, 0.5 and 0.2 bar.The SEI per molecule was varied for 0.9 bar between 0.3 and 5 eV/molecule, and for the two other pressures between 0.6 and 2 eV/molecule.
The calculated conversions (and energy efficiencies) are found to be always in good agreement with experiment, close to the experimental error bars.At the same time, deviations between the calculated and measured gas temperatures are always large.For 0.9 and 0.5 bar the maximum temperature T max is in most cases overestimated by the reference model by 20%-30%.Comparison of the radial and axial temperature profiles shows that besides the deviation of absolute values the calculated profiles appear more peaked than the experimental ones.For 0.2 bar the calculated T max can exceed the experimental value by almost a factor 2, getting significantly larger than 8000 K, which makes the results nonphysical.Nevertheless, even in these cases the calculated conversion agrees well with the measurements.
The thermo-chemical model requires the characteristic diameter and length of the heating zone d pl and l pl as input parameters.These are estimated in experiments from the intensity of plasma radiation and have relatively large uncertainty.Therefore, the sensitivity of the numerical solutions with respect to increasing and reducing d pl and l pl by a factor 1.5 has been investigated.It has been found that the calculated T max and the temperature profiles are sensitive with respect to the assumed d pl and l pl , but the conversion is not sensitive.
The results of the present study confirm that the conversion of CO 2 observed in contracted microwave plasmas can be explained by a purely thermo-chemical mechanism, which coincides with the analysis presented previously in [2-4, 16, 17].The results emphasize the important role played by radial diffusion in the net conversion process.That is, the diffusion of CO and O outwards from the hot plasma core and diffusion of CO 2 towards the hot core.A steep temperature gradient of the order of 10 6 Km −1 along the diffusion path that forms at the boundary of the hot core leads to-partial-mitigation of the backward recombination of CO, which produces the net chemical effect.The essential role of radial transport in the present configuration with a straight cylindrical channel is demonstrated explicitly by showing that the calculated conversion does not reduce even under the assumption of zero flow in the hot core near the axis.
In summary, the present study provides initial evidence that the thermo-chemical approach to modelling of the CO 2 to CO conversion at (near-) atmospheric pressure works.Since a large discrepancy between the model and the experiment was found for the gas temperature, this result is yet to be corroborated by a complete transport model that solves the full Navier-Stokes equation.In particular, the impact of the convective cells which form in the vortex stabilized flows [6,27,28] on the chemical process has to be clarified.Also, the benchmark performed here is restricted only to the flow configuration with straight channel in the effluent.The next mandatory step that has to be completed before thermo-chemical modelling can be used for practical process optimization is validation against experiments with nozzles in the effluent.In the present paper the turbulent transport is left out of the scope completely because the experimental conversion could be reproduced with purely molecular transport.In other works [17,31] vortex-induced turbulent mixing was identified as a potentially important driving force for the non-equilibrium process, whose enhancement in the nozzle flows could be the reason for their effect on the net conversion.module with separate unit-tests.The classes that implement spatial discretization were unit-tested by the Method of Manufactured Solutions.Heat transfer problems with the analytic solutions 1D in z or r direction were used for integrated tests of the solver of the whole system (1), ( 2) and (4).With a single component as well as several identical components (which must produce exactly the same solution).The convergence is monitored by controlling residuals of (1) and ( 2) and global balances.In all calculations discussed in the paper the errors in the global particle balance of any chemical component are always smaller than 2 • 10 −3 of the total flow rate, and the errors in the global energy balance are always smaller than 10 −3 of the total input power.The calculations were performed on relatively coarse grids with 120×40 cells.Grid size independence of the solution was checked by repeating the six cases listed in table 4 on larger 240×80 cells grids.
Flows with chemical reactions are mathematically stiff problems.To simulate them special operator splitting schemes are usually applied [54].This was not done in the present study-with negative consequences for numerical performance.To reach the level of convergence of the global balances reported above the time-step 10 −5 sec or smaller was required.At the same time, to achieve the steady state for the problems in question a simulation had to run over 0.1-10 sec of physical time.As a result, one simulation even on a 120×40 grid takes several hours or more.This experience emphasizes that applying appropriate operator splitting, such as the popular Wu-scheme [55], is highly desirable for practical simulations of thermo-chemical plasma conversion even if only a relatively small number of components, five in the present case, are taken into account.

Appendix B. Extract from kinetic theory
According to [40], equations (8.2-9) and (8.2-11), the following expressions can be derived for the coefficient of selfdiffusion D and thermal conductivity λ of a single component gas of spherically symmetric molecules without internal degrees of freedom: In these expressions, σ is the effective rigid sphere collision cross-section, c v is the molar heat capacity at constant volume, Ω (1,1) * and Ω (2,2) * are the normalized collision integrals.The latter are the total collision integrals Ω (1,1) , Ω (2,2) defined in [40], section 7.4d, equation (7.4-34), divided by the collision integrals calculated for the rigid sphere type of interaction for which the integration can be performed analytically: [40], equation (8.where A * = Ω (2,2) * Ω (1,1) * is shown to be nearly independent of T. The results of the calculation of A * for Lennard-Jones potential can be found in [40], table I-N.They demonstrate that in a very wide temperature interval A * only varies between 1.05 and 1.15.Therefore, in the calculations A * = 1.1 is always taken.The accuracy provided by (B.3) with D calculated by Fuller scaling ( [39], equation ) was tested by applying it to calculate λ for O atoms for which the accurate quantummechanical calculations [44] are available.This test has shown that in the temperature range 2000-10 000 K this method overestimates the results of [44] by 40% at maximum, which is a reasonable agreement.

Figure 2 .
Figure 2. Comparison of the model with the experiments[4,22] for pressure 0.9 bar and different flow rates.SEI is defined by (13), 'conversion' χ is defined by (14) and 'energy efficiency' η by(15).Dashed lines are obtained with the nominal value of d pl increased by a factor 1.5 and the dotted lines with d pl decreased by a factor 1.5, see section 4.1.

Figure 3 .
Figure 3.Comparison of the model with the experiments[4,22] for pressure 0.5 bar and flow rate 20 slm.SEI is defined by (13), 'conversion' χ is defined by (14) and 'energy efficiency' η by(15).Dashed lines are obtained with the nominal value of d pl increased by a factor of 1.5 and the dotted lines with d pl decreased by a factor of 1.5, see section 4.1.
Figure 5. Spatial distribution of the gas temperature in the case of 0.9 bar, 0.9 kW and 10 slm (same case as in figure6in[4]).(a) Radial profile at z =28 mm, which corresponds to '58 mm above the resonator bottom' in[4]; (b) profile along the axis, r = 0. Solid lines are the reference model of section 3.1, while dashed lines correspond to different cases listed in table 4. Dotted vertical line in sub-figure (b) marks the position of the cross-section shown in (a).

Figure 7 .
Figure 7. Examples of cumulative local energy efficiency for different values of SEI for pressure 0.9 bar; η (z) is defined by (13), and ηreox (z): equations (17) and (18), respectively.Vertical dashed line is the cross-section z =28 mm, the same as in figure5(a).

Figure 9 .
Figure 9. Example of the cooling rate of O atoms at three selected cross-sections, case 0.9 bar, 10 slm, 0.9 kW.(a) is the temperature profile; (b) is the effective diffusive velocity of O atoms; (c) is the cooling rate 'experienced by' the O atoms.'z = 0 mm' is the middle of the heat source and 'z = 28 mm' is the cross-section shown in figure 5(a).

Table 1 .
List of chemical processes.

Table 2 .
List of reactions.

Table 4 .
Maximum temperature Tmax, in kK, calculated on different assumptions for selected model cases.

Table 5 .
Comparison of the simulations with different values of the wall temperature Tw.