A boundary value problem of heat transfer within DBD-based plasma jet setups

1.   Introduction

Cold atmospheric pressure plasmas (CAPPs) refer to the ionized gaseous phase, being out of thermodynamic equilibrium and thus rich in reactive species, which is sustained close to the ambient pressure. On the other hand, dielectric barrier discharges (DBDs) prevent transition to arcs and thus to electrical discharge regimes associated with raised temperatures. Hence, unique physicochemical features may be combined by means of DBD-based CAPPs, and this concept has led to cutting-edge technologies in the fields of material processing, aeronautics ^[1], environmental remediation ^[2], etc. Of special importance is the plasma biomedicine field, where DBD-based CAPPs are explored as tools for wound healing, sterilization ^[3], cancer treatment ^[4], tooth bleaching, etc.

Towards such biomedical applications, DBD-based CAPPs are usually engineered in the form of the so-called "plasma jets" ^[5], where a gas is channeled into the atmospheric air through a dielectric capillary tube, while it is simultaneously subjected to a time-varying electric field. The latter is developed by driving a pair of electrodes of appropriate geometry (e.g., coaxial) with high voltages of various waveforms (e.g., ac-continuous or ns-pulsed in the kHz range). Thus, a complicated interplay between charged species, excited neutrals, chemical radicals, photons, electric fields ^[6], fluid fields ^[7,8,9,10], etc. and the bio-specimen itself, takes place. However, the benefits of the interaction between this plasma reactive phase and the specimen may be compensated by the rise of the gas temperature with respect to the specimen's thermal tolerance. In other words, dominant design and operation parameters that are considered during the optimization of the plasma jet setups must be correlated with the thermal fields of various heating mechanisms ^[11].

Several works have been devoted to the thermal properties of CAPPs, whereas they mostly refer to experimental and numerical techniques ^{[12,13,14,15,16,17,18,19,20]}. Herein, the interest in investigating the thermal properties of plasma jets, as a special case of CAPPs, is combined with the challenge of producing an initial, boundary value problem of heat conduction within such setups. The ultimate motivation for this research is to obtain analytical formulae for the temperature field in the system that can be employed in plasma simulations, instead of time-consuming numerical solutions. These solutions will provide an easy means of estimating the spatial distribution of the temperature for different geometric and operation characteristics, which are important factors for the properties of the plasma. Furthermore, novel physical results can be derived from analytical techniques and they can be utilized for the verification of numerical results. Under these aspects, the present article has endeavored to show that an analytical approach to problems involving CAPP jet setups is conceivable.

Namely, the theoretical analysis considers a two-region problem, which corresponds to the dielectric barrier tube of the DBD and the external grounded electrode, wherein our goal is to calculate the temperature distribution in the axisymmetric configuration of the setup. In order to construct the related initial and boundary value problem and proceed to the solution, our analytical methodology is primarily based on the selection of a suitably located cylindrical coordinate system ^[21] for modeling purposes with respect to the current investigation. Therein, the domain of field activity is divided into the two subsectors of the dielectric and the electrode-metal, in which classical heat diffusion partial differential equations ^[22] are considered. These are supplemented by the appropriate conditions, i.e., the Dirichlet and Neumann continuity conditions on the interface of the two solids, an imposed given heat conduction rate assumed on the boundary of the dielectric with the working gas, and the Robin-type heat convection condition on the boundary of the electrode with the outer environment. Otherwise, an initial room temperature condition is imposed for both the areas of thermal activity, completing the well-pose of the problem. Henceforth, the well-known method of separation of variables is applied to both the heat equations, while a standard technique, based on the method of asymptotic kernels ^[23], is necessary for handling the non-homogeneous boundary conditions. Following these steps, the time-dependent temperature fields in each region are expanded in terms of cylindrical eigenfunctions ^[23] and decaying time exponents, in which rotational symmetry is implied, featured by the conditions. Then, the aforementioned problem is solved straightforwardly, resulting in either analytical or semi-analytical expressions for the incorporated unknown constant coefficients, while a pair of dual relations of Fourier-Bessel series ^[24] for a single unknown also appear due to the different boundaries. The latter are merged to one and only relationship, which is handled appropriately so as to be reduced to an infinite algebraic linear system of equations for the implicated set of constants, by taking its projection to Bessel and Neumann functions. Therein, the obtained handy system is manipulated via well-known cut-off techniques, while the temperature fields are provided in a compact fashion via closed-type convergent series expansions and they are implemented so as to show their response graphically. The results reveal the thermal behavior both in the bulk and on the surfaces of the solid materials, which is of interest for future elaboration either using alternative approaches, e.g., by solving in the Laplace domain and transforming back to the time domain ^[25,26] or proceeding to possible experimental processing. In order to demonstrate the efficiency of the present method, the analytical solution is juxtaposed with numerical results, being obtained using the advanced Finite Integration Technique (FIT) ^[27,28,29], wherein an excellent agreement is obtained.

2.   Heat Transfer Configuration and Analytical Solution

Let us consider two very long coaxial cylinders, the dielectric tube (alumina) of thermal diffusivity ${a_d}$ with inner radius ${R_d}$ and the electrode (brass) of thermal diffusivity ${a_e}$ with inner radius ${R_e}$. Furthermore, ${k_d}$ and ${k_e}$ are the respective thermal conductivities (see Figure 1). The system is initially set at the same temperature with the free-environment, whose convective heat transfer coefficient is $h$, contains air at room temperature ${T_\infty }$, while plasma kinetics produces a known conduction heat $q$. Hence, in terms of a properly adjusted cylindrical coordinate system $(\rho, \varphi, z)$, in which the $z$ -variable is eliminated as a fair approximation of the big length of the plasma jet setup, we consider two different temperature distributions, the first one within the dielectric ${T_d}(\rho, \varphi, t)$ and the second one within the electrode ${T_e}(\rho, \varphi, t)$. Then the boundary and initial value problem that we have to solve admits two separate heat partial differential equations in cylindrical geometry, which assume

and

both defined for every $\varphi \in [0, 2\pi)$ and $t > 0$, accompanied by the appropriate boundary conditions

and the initial condition of common temperature for the dielectric tube and the electrode

completing thus the set of a well-posed boundary value problem that describes the heat transfer within our system.

Our first task in order to solve (1)-(7), is to transform (3) and (6) into homogeneous conditions. To this end, we use the method of asymptotic kernels, according to which the solution under consideration is rewritten as the summation of a general solution and a particular time-independent solution (resembling the asymptotic behavior as $t \to + \infty$) via

and

According to our argumentation, functions ${K_d}(\rho, \varphi)$ and ${K_e}(\rho, \varphi)$ are chosen, such as to solve the problems

with boundary condition

while

with boundary condition

On the other hand, the main fields ${\tilde T_d}(\rho, \varphi, t)$ and ${\tilde T_e}(\rho, \varphi, t)$, implying (10) and (12) into relations (1)-(7) and according to decompositions (8) and (9), satisfy the boundary value problems for every $\varphi \in [0, 2\pi)$ and $t > 0$

and

with boundary conditions

and initial condition of the system

whereas the common temperature is set by the final solution.

In the sequel, we work as follows. Due to conditions (11) and (13), functions ${K_d}(\rho, \varphi)$ and ${K_e}(\rho, \varphi)$ depend only on the $\rho$ -variable, i.e., ${K_d}(\rho, \varphi) \equiv {K_d}(\rho)$ and ${K_e}(\rho, \varphi) \equiv {K_e}(\rho)$, respectively, henceforth (10) and (12) become eventually ordinary differential equations, since $\partial {K_d}(\rho)/\partial \varphi = 0$ and $\partial {K_e}(\rho)/\partial \varphi = 0$, whereas applying (11) and (13) independently of $\varphi$ -variable, it is readily obtained

and

where ${c_d}, {c_e} \in \mathbb{R}$ are arbitrary constants, appropriately chosen within the forthcoming steps. On the other hand, in order to manipulate the fields ${\tilde T_d}(\rho, \varphi, t)$ and ${\tilde T_e}(\rho, \varphi, t)$ into the partial differential equations (14) and (15), respectively, we apply the well-known method of separation of variables, according to which we suppose a solution of the general form

for $\varphi \in [0, 2\pi)$ and $t > 0$, wherein, since (14) and (15) retain the same structure, it is doable to perform the procedure once. Hence incorporating two separation constants $\kappa$ and $\lambda$, we substitute (23) into either (14) or (15) to get

accordingly to the domain of interest and

whereas either $\rho \in ({R_d}, {R_e})$ ($j = d$) or $\rho \in ({R_e}, {R_s})$ ($j = e$). Solving the second equality of (25) and demanding the proper periodic azimuthal behavior, we are led to

whilst by combination of (26) with (24) for $\kappa = - \mu _j^2$ with $j = d, e$, we obtain the Bessel-type ordinary differential equation

whereas, again, either $\rho \in ({R_d}, {R_e})$ ($j = d$) or $\rho \in ({R_e}, {R_s})$ ($j = e$), whose solutions are the Bessel functions ${J_n}$ and the Neumann functions ${N_n}$, that is

Finally, the time-dependence for $\kappa = - \mu _j^2$ with $j = d, e$ within (24) yields

Gathering all the formulae obtained above through (26), (28) and (29), the general solutions are provided via

for every $\varphi \in [0, 2\pi)$ and $t > 0$, where $j = d$ if $\rho \in ({R_d}, {R_e})$ and $j = e$ if $\rho \in ({R_e}, {R_s})$. The symbolic feature " $\int {\sum\limits_{{\mu _j}} \cdots }$ " corresponds either to series definition in the case where ${\mu _j}$ take discrete values or to integral notation if ${\mu _j}$ are continuous parameters for $j = d, e$. Due to the nature of conditions (17) and (18), the azimuthal dependence must be omitted by imposing ${c_n} = {d_n} = 0$ for $n \geqslant 1$ and keeping the terms for $n = 0$. Consequently, by the new definition of $A_0^{{\mu _j}} \equiv a_0^{{\mu _j}}{c_0}$ and $B_0^{{\mu _j}} \equiv b_0^{{\mu _j}}{c_0}$, the general solution (30), which corresponds to the solution of (14) and (15), reduces to

and

whereas the set of $A_0^{{\mu _d}}$, $B_0^{{\mu _d}}$ with parameter ${\mu _d}$ and $A_0^{{\mu _e}}$, $B_0^{{\mu _e}}$ with parameter ${\mu _e}$ must be evaluated by (16)-(20) with respect to (21) and (22).

Bearing in mind that $\partial /\partial \rho = {\mu _j}\partial /\partial ({\mu _j}\rho)$ and the derivative property ${J'_0}({\mu _j}\rho) = - {J_1}({\mu _j}\rho)$ and ${N'_0}({\mu _j}\rho) = - {N_1}({\mu _j}\rho)$ over the argument for $j = d, e$, putting (16) into (31), we reach

while applying (19) into (32), we get

for every value of the parameters ${\mu _d}$ and ${\mu _e}$, respectively. Substituting (33) into (31) and (34) into (32), we obtain the homogeneous part of the solutions via

for $\rho \in ({R_d}, {R_e})$ and $t > 0$, while

for $\rho \in ({R_e}, {R_s})$ and $t > 0$, wherein $C_0^{{\mu _d}} \equiv A_0^{{\mu _d}}{J_1}({\mu _d}{R_d})$ and $C_0^{{\mu _e}} \equiv A_0^{{\mu _e}}\left[{h{J_0}({\mu _e}{R_s}) - {k_e}{\mu _e}{J_1}({\mu _e}{R_s})} \right]$ are the newly defined constants, respectively. We precede with the interface conditions (17) and (18), which, by virtue of (35) and (36), as well as (21) and (22), yield

and

with $t \geqslant 0$, respectively. The utility of relationships (37) and (38) is twofold, since the proper matching for any $t \geqslant 0$ acquires two things. Firstly, we are forced to consider that

which immediately offer us

providing the two arbitrary constants in the fields (21) and (22). Once done, implying (39) into (37) and (38), we proceed to the second step, i.e., in order to retain the same exponents in the time variable, we imply

and therein, for notational convenience, we define as $\mu \equiv {\mu _e}$, therefore ${\mu _d} = \mu \sqrt {{\alpha _e}/{\alpha _d}}$. Hence, following the new parameter notation, we interchange $F_0^\mu \equiv C_0^{{\mu _d}}$ and $G_0^\mu \equiv C_0^{{\mu _e}}$ for the same reason and we readily define convenient functions that enter into the fields (35) and (36), as well as into conditions (37) and (38), those being

and

Within the same aspect, we introduce the ratio $C_0^\mu$ of the newly-defined defined constants $C_0^{{\mu _d}}$ and $C_0^{{\mu _e}}$ as

so from (37) and (38) we readily obtain

where the prime denotes differentiation with respect to the argument for each case. The latter are two transcendental conditions, being solved numerically for the evaluation of the ratio $C_0^\mu$ from (44) and the parameter $\mu$. Once done, parameter ${\mu _d}$ is readily recovered via (41), while we remain with only one coefficient to calculate, either $F_0^\mu$ or $G_0^\mu$, since the ratio (44) is a known quantity. Incorporating the above relationships (37)-(45) to the fields (35) and (36), we may directly rewrite them as

and

Substituting (46) and (47) into (8) and (9) (taking into account the azimuthal independence), while attaching (21) and (22) with (40), we get

for every $\rho \in [{R_d}, {R_e}]$ and $t > 0$, while

for every $\rho \in [{R_e}, {R_s}]$ and $t > 0$.

Thus, the initial conditions (7) or equivalently (20), by virtue of (48) and (49), render the expressions

and

which comprise two relations that, beyond the fact that they contain the $\rho$ -variable for $\rho \in [{R_d}, {R_s}]$, they also depend upon the physical and geometrical characteristics of the system under consideration. The last tricky step is to calculate the remaining coefficient $G_0^\mu$ from the dual series conditions (50) and (51), where despite their simplicity, special treatment is required according to the following. Since we attain two relationships with the same unknown, our first task includes the merging of (50) and (51) into a single expression, which means that they can be written as

where

and

that conveniently reduces our problem to that of solving one condition (52) with the single unknown coefficient $G_0^\mu$. Our second task is now to eliminate the $\rho$ -variable from (52) and towards this direction we multiply relationship (52) by ${J_0}(\mu '\rho \, )\rho d\rho$ and we integrate overall from ${R_d}$ to ${R_s}$ in an appropriate manner, based on integration either from ${R_d}$ to ${R_e}$ or from ${R_e}$ to ${R_s}$, according to (53) and (54). During this procedure, we observe that functions (42) and (43), which appear on the left-hand side of (52), due to definition (53), involve the zeroth-order Bessel functions ${J_0}$ and Neumann functions ${N_0}$, which are defined both for the same and for different arguments. Hence, it is apparent that we need to analytically calculate the definite integrals

for every value of $\mu$, $\mu '$ and

for $\mu \ne \mu '$, while

when $\mu = \mu '$, being given in a convenient fashion. Otherwise, the right-hand side of (52) includes simple logarithmic and constant functions, due to definition (54), consequently, we need to evaluate the trivial formulae

and, similarly,

as well as

where the latter is defined for $\rho \in [{R_d}, {R_s}]$. In the sequel, the outcomes (55)-(63), whose evaluation is mainly based on the use of standard recurrence relations of Bessel functions ^[23], are readily implemented into (52) (after being multiplied by ${J_0}(\mu '\rho \, )\rho d\rho$ and integrated with respect to the interval $[{R_d}, {R_s}]$), taking into account (53) and (54) as well, hence, we obtain straightforwardly

wherein

and

with $\mu, \mu ' \in \mathbb{R}$. Here, we observe that in order to manipulate condition (64), we are obliged to set the parameter $\mu$ so as to obtain discrete values, that is $\mu \equiv {\mu _{0, m}}$ for every $m \geqslant 1$ (similarly $\mu ' \equiv {\mu _{0, m'}}$ with $m' \geqslant 1$), whereas the zeroth order of the Bessel functions is readily implied to the notation. Then, the aforementioned series-integral symbolism becomes now a series, due to the eigenvalue property of parameter $\mu$, thus $\int {\sum\limits_\mu \cdots } \, \, \, \to \, \, \, \sum\limits_{m = 1}^{ + \infty } \cdots$, and therein our analysis is restricted to a classical infinite series semi-analytical solution.

Under this aim, we recapitulate for convenience our findings in the present analysis, hence the temperature distributions within the dielectric tube and the electrode that come from expressions (48) and (49), admit

for every $\rho \in [{R_d}, {R_e}]$ and $t \geqslant 0$, while

for every $\rho \in [{R_e}, {R_s}]$ and $t \geqslant 0$, where formulae (42) and (43) become

and

for $m \geqslant 1$, while $C_0^m$ and ${\mu _{0, m}}$ are the solutions of the two transcendental relationships (45), which provide us with

wherein $C_0^m$ for $m \geqslant 1$ is obtained analytically from the left-hand side of (73), while the evaluation of ${\mu _{0, m}}$ for $m \geqslant 1$ is subject to further numerical handling (see in the sequel for a more detailed and accurate explanation on that matter), as it is revealed from the right-hand side of the set of relationships (73). Otherwise, $G_0^m$ for $m \geqslant 1$ can be evaluated as the outcome of the infinite system of linear algebraic equations

resulting from relation (64) with formulae (65)-(68) and (55)-(63), in which the involved mathematical quantities assume the following complicated, yet handy, expressions in terms of the geometrical and physical parameters of the system, i.e.,

and

while

with

and

for $m \ne m'$, whilst

for $m = m'$, all above provided for every $m, m' \geqslant 1$.

This sequence of systems of equations (74) with (75)-(84) appears often in problems of mathematical physics and our effort is then limited in developing a solution, derived from a partial differential equation problem with mixed boundaries. To this end, we introduce an idea, based on usual cut-off techniques in order to solve infinite linear systems. According to this methodology, we truncate the series (74) by the imposition of an indispensable upper limit $M \in {\mathbb{N}^ * }$, so as the series become finite, due to the fact that we may write $\sum\limits_{m = 1}^{ + \infty } \cdots \, \, \, \to \, \, \, \sum\limits_{m = 1}^M \cdots$, and facilitate the procedure. Thereafter, (74) reduces to a finite algebraic linear system of equations if we force the number of the incorporated equations to be finite as well, that is $m' = 1, 2, ..., M'$, wherein $M' \in {\mathbb{N}^ * }$. This system becomes quadrature only if $M = M' \equiv L \in {\mathbb{N}^ * }$ and once done, it is rewritten in matrix form as

in which $m, m' = 1, 2, ..., L$, while $\mathbb{A}$ is the squared-type invertible matrix of the coefficients of the unknowns, ${\bf{x}}$ is the vector of the unknown coefficients and the ${\bf{b}}$ is vector of the known constants. Thus, each time we repetitively solve quadrature systems numerically for increasing $L$ and the standard cut-off is applied in order to decide the truncation amplitude ${L_{\max }}$, for which the number of the evaluated constant coefficients are sufficient to obtain the expected accuracy in order for the series to converge. The solution of the system (85) yields the sought coefficients ${\bf{x}}$, i.e., the unknown constant coefficients $G_0^m$ for every $m = 1, 2, ..., {L_{\max }}$.

The final solution of the problem is eventually achieved via the temperature functions (69) and (70), followed by the conveniently chosen functions (71) and (72), in which the constant coefficients $C_0^m$ and ${\mu _{0, m}}$ for $m = 1, 2, ..., {L_{\max }}$ are the solutions of relationships (73), while the set of constants $G_0^m$ for $m = 1, 2, ..., {L_{\max }}$ that satisfies the easily amenable infinite set of linear relationships (74) with (75)-(84), are given by (85), i.e., ${\bf{x}} = {\mathbb{A}^{ - 1}}{\bf{b}}$, when the incorporated inverse matrix ${\mathbb{A}^{ - 1}}$ is given. Under the aim of presenting a compact expression for the temperature field, we utilize the definition of the Heaviside step function $H(x) = 1$ for $x \geqslant 0$, otherwise zero, and we incorporate fields (69) and (70) into

which provides the temperature in both the dielectric tube and the external electrode in a handy analytical fashion. Therefore, the value of the temperature on the outer cover of the electrode for $\rho = {R_s}$, being of our main concern, is provided theoretically according to (86) and (70) (or from (70) directly) as

where ${T_s}(t) \equiv {T_e}({R_s}, t)$ for any $t \geqslant 0$, while in practice the upper infinite limit of the series (87) must be truncated to ${L_{\max }}$ for engineering applications.

In order to demonstrate the above analysis, we devote the rest of the analysis to the numerical implementation and validation of the obtained analytical expressions for the temperature distribution in the domains of interest, i.e., within the dielectric tube and the external electrode. Towards this direction, we provide specific indicative values of the involved physical parameters and properties, appearing during our analysis ^[22]. To this end, we consider a representative value for the heat conduction $q = 18624\, {\text{W}}{{\text{m}}^{ - 2}}$ in plasma jets and $h = 13.2\, {\text{W}}{{\text{m}}^{ - 2}}{{\text{K}}^{ - 1}}$ for such systems, while we have ${k_d} = 30\, {\text{W}}{{\text{m}}^{ - 1}}{{\text{K}}^{ - 1}}$ and ${k_e} = 110\, {\text{W}}{{\text{m}}^{ - 1}}{{\text{K}}^{ - 1}}$, wherein we directl calculate ${a_d} = {k_d}/{\rho _d}{c_{p, d}} = 9.9 \times {10^{ - 6}}\, {{\text{m}}^2}{{\text{s}}^{ - 1}}$ and ${a_e} = {k_e}/{\rho _e}{c_{p, e}} = 3.3 \times {10^{ - 5}}\, {{\text{m}}^2}{{\text{s}}^{ - 1}}$, since it is ${\rho _d} = 3800\, {\text{kg}}{{\text{m}}^{ - 3}}$, ${\rho _e} = 8600\, {\text{kg}}{{\text{m}}^{ - 3}}$ and ${c_{p, d}} = 800\, {\text{J}}\, {\text{k}}{{\text{g}}^{ - 1}}{{\text{K}}^{ - 1}}$, ${c_{p, e}} = 388\, {\text{J}}\, {\text{k}}{{\text{g}}^{ - 1}}{{\text{K}}^{ - 1}}$, respectively for each case and the geometrical characteristics of the system yield ${R_d} = 5.7 \times {10^{ - 4}}{\text{m}}$, ${R_e} = 1.25 \times {10^{ - 3}}{\text{m}}$ and ${R_s} = 5 \times {10^{ - 3}}{\text{m}}$. Note that these values assume units in ${\text{K}}$ for the temperature, however our plots use the corresponding units in ${}^oC$. Given these representative values, we initially have to numerically solve the right-hand-side of (73) so as to find the proper values of ${\mu _{0, m}}$ for $m \geqslant 1$ that force the new-defined function

to become zero, i.e., $h(x) = 0$. To ensure that all the roots were computed, the sign of the equation was checked at every interval between two consecutive roots. It was found that the roots do not follow a specific pattern, so it was necessary to use very small steps in the control procedure. These roots are given in Figure 2, hence we have solved the equivalent relationship of the second equation in (73) for the evaluation of ${\mu _{0, m}}$, $m \geqslant 1$.

Having found the roots ${\mu _{0, m}}$ for $m \geqslant 1$, it is necessary to calculate the terms (75)-(84), which are necessary for computing the linear system of algebraic equations (85). For different values of $m$ and $m'$ the procedure was straightforward in (75)-(82), while for the diagonal elements, where $m = m'$, it is obvious that we have to refer to relations (83) and (84). The resulting coefficients for the system (85) are obtained via the application of the cut-off method. Thereafter, the temperature distribution $T$ in the dielectric tube and the external electrode follows from relationship (86), wherein its variation with respect to time $t$ at the middle thickness of the dielectric tube $({R_d} + {R_e})/2$ and the external electrode $({R_e} + {R_s})/2$, as well as with respect to the radius $\rho$ of the system at $3000\, s$, which is approximately steady state, is implemented and depicted in Figure 3 and Figure 4, respectively. Obviously, the derivation of the temperature refers to the transient heat transfer problem, as shown in Figure 3, even though the verification in Figure 4 was for steady temperature. The results of the finalized semi-analytical solution have been compared against numerical simulations carried out using the Finite Integration Technique (FIT). FIT is a grid-based method very similar to the Finite Difference in Time Domain (FDTD), originally introduced for the solution of the Maxwell equations, using a system of staggered orthogonal grids ^[27] and which has been extended to acoustics, elastodynamics (the extension is known as Elastodynamics Finite Integration Technique (EFIT)) ^[28], as well as to heat equation problems ^[29]. We indicate that for numerical convenience purposes the temperature $T$ has been re-gauged so as the initial temperature of the system is $0\, {}^oC$ instead of room temperature.

The analytical solution computations are performed in two steps. The first (and most tedious one) concerns the calculation of the eigenvalues. For the (non-optimized) root finding algorithm used in the present implementation, this time was $17\, s$. The second step comprises the inversion of the system matrix for the estimation of the development coefficients and the evaluation of the series. The computation time for this second part has been counted to $80\, ms$. These computational times have to be compared with the respective ones required for the numerical solution of the problem. Since we are dealing with a diffusion problem in the time domain, the considered time window has to be sampled and the discrete linear system of equation obtained by the application of the method's grid (either using FIT or FEM) must be inverted at each time instance (implicit formulation). For the specific example examined here, the total computational time reached is $4\, s$. Clearly, there is a computational overhead when we need to evaluate the eigenvalues of the analytical solution. This calculation however has to be carried out only once, independently from the excitation signal. Every other evaluation is of the order of $ms$, whereas the numerical solution has to be re-run for each new excitation. The gain in performance increases dramatically when moving at higher dimensions, the cost, however, that one has to pay is that the analytical solution becomes cumbersome.

At this stage, we provide the necessary discussion about the relative errors during our methodology. There are three sources of error that enter in the computation of the analytical solution, i.e., the error of the eigenvalue calculation, the numerical evaluation of the special functions and the truncation of the development series. The first one is controlled by the numerical solution of the transcendental equation (88), which is determined by the method threshold and the round-off error. The evaluation error of the special (Bessel) functions is determined by the specific numerical implementation of the computational platform (in our case Matlab) to which one usually has little control, unless dedicated implementations are used. With higher orders, overflows occur; therefore, in our implementation, the Bessel functions are always evaluated in form of ratios. The function evaluation error can be assumed having negligible impact to the results. The third one (number of series term taken into account) is the one with the strongest impact to the solution precision. There is a trade-off between precision and efficiency and one has to resort to numerical experimentation in order to pick-up the optimal number. In addition, care must be taken with higher-order terms, since the root finding algorithm becomes less accurate due to overflows (despite the normalization of the Bessel functions). Taking all these factors into consideration and after numerical experimentation, we have concluded that a number of 36 terms gives the most satisfactory results in terms of the computational time and precision. Concerning the issue of error variation as function of the position raised, one can expect that it may be higher at the interfaces, where the continuity is imposed by matching the series terms.

The analytical solution is compared with results provided by the FIT method for a given test configuration, whose parameters have been carefully chosen in order to represent a realistic case. Since the two calculations (analytical and numerical one) are independent and given that the FIT computations have been carried out using a generic code, validated for a number of different test problems (in fact the FIT solver is part of the CIVA commercial simulation platform), we accept that the agreement of the two data sets validates the correctness of the presented solution. It can be seen that the calculated initial condition displays a great variation, probably due to the over defined nature of the algebraic system. According to the problem formulation and expression (87), the temperature on the surface of the electrode as $t \to + \infty$ (or large time scaling for applications) converges to the temperature $160.84\, {}^oC$, which is the plateau in Figure 3. On the other hand, for $t = 0$ (initial condition), the room temperature is readily recovered from (87), as it is shown in Figure 3.

The analytical technique presented herein, deals with the 1D three-layer problem of the plasma cell, described in Figure 1. It can be easily extended to address finite multilayer configurations with a consequent increase, however, of the solution complexity. A basic limitation of the method is that it remains applicable to piecewise homogeneous geometries, i.e., material gradients cannot be tackled without drastic modifications, such as the introduction of hybrid development bases mixing the herein applied Fourier basis with a spatial mesh.

3.   Conclusions and discussion

The present work was devoted to the presentation of an analytical technique for the determination of the temperature in a representative plasma jet reactor. We studied the heat transfer process from the plasma area to the outer environment through the solid components of the reactor, i.e., the dielectric tube and the grounded electrode.

The physical problem itself was mathematically formulated with respect to the cylindrical geometry, which is employed, since it fits the plasma jet reactor. We separated the area of thermal activity into two distinct domains, matching the one on the dielectric tube and the other of the external electrode, considering them as coaxial cylinders. Under this aim, we constructed appropriately the corresponding initial and boundary value problems in each domain with either Dirichlet or Neumann conditions and the solution technique was based on the method of separation of variables in the cylindrical coordinate system and the imposition of a particular method of asymptotic kernels for the manipulation of the non-homogeneous boundary conditions. The time-dependent temperature fields within the tube and the electrode were constructed in terms of cylindrical harmonic eigenfunctions and the analytical solution was obtained in view of the restriction of rotational symmetry, due to the nature of the boundary conditions. During the analysis, all the implicated unknown constant coefficients were evaluated straightforwardly, except one. For this last unknown, we were obliged to solve a pair of dual Fourier-Bessel series relations, which were handled simultaneously as a single relationship and led to infinite linear systems for the unknowns via integration techniques, involving mostly Bessel and Neumann eigenfunctions among other simple functions. These systems were solved with the aim of a repetitive cut-off method until the expected accuracy was achieved.

The efficiency of the analysis was demonstrated by implementing numerically the expressions of the final formula, which satisfy the initial and boundary value problem under consideration, while a Finite Integration Technique (FIT) was applied for the validation of the methodology. It was shown that the distribution of the temperature behaves as expected, as far as the spatial dependence is concerned, while as the time evolves, it converges.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.