Order Reduction of the Radiative Heat Transfer Model for the Simulation of Plasma Arcs

# Order Reduction of the Radiative Heat Transfer Model for the Simulation of Plasma Arcs

L. Fagiano, R. Gati L. Fagiano and R. Gati are with ABB Switzerland Ltd., Corporate Research, Segelhofstrasse 1K, Baden-Daettwil, CH-5405. Corresponding author: L. Fagiano (lorenzo.fagiano@ch.abb.com)
###### Abstract

An approach to derive low-complexity models describing thermal radiation for the sake of simulating the behavior of electric arcs in switchgear systems is presented. The idea is to approximate the (high dimensional) full-order equations, modeling the propagation of the radiated intensity in space, with a model of much lower dimension, whose parameters are identified by means of nonlinear system identification techniques. The low-order model preserves the main structural aspects of the full-order one, and its parameters can be straightforwardly used in arc simulation tools based on computational fluid dynamics. In particular, the model parameters can be used together with the common approaches to resolve radiation in magnetohydrodynamic simulations, including the discrete-ordinate method, the P-N methods and photohydrodynamics. The proposed order reduction approach is able to systematically compute the partitioning of the electromagnetic spectrum in frequency bands, and the related absorption coefficients, that yield the best matching with respect to the finely resolved absorption spectrum of the considered gaseous medium. It is shown how the problem’s structure can be exploited to improve the computational efficiency when solving the resulting nonlinear optimization problem. In addition to the order reduction approach and the related computational aspects, an analysis by means of Laplace transform is presented, providing a justification to the use of very low orders in the reduction procedure as compared with the full-order model. Finally, comparisons between the full-order model and the reduced-order ones are presented.

Keywords: Arc simulations, Radiative heat transfer, Model order reduction, Nonlinear estimation, Nonlinear model identification

## 1 Introduction

The switching performance of circuit breakers depends strongly on the behavior of the electric arc that originates when the contacts are opened in presence of relatively large electric current values [27, 7]. In turn, the arc dynamics are influenced by multiple interacting physical phenomena which, together with the short timescale of the arcing event and the large values of temperature and pressure, increase the complexity and difficulty of understanding, carrying out experiments, and deriving numerical models of the switching behavior. Computational fluid dynamic (CFD) approaches are being used in both public and private research efforts to simulate the time evolution of the plasma that carries the current during the interruption process, see e.g. [11, 6, 2, 13]. The CFD simulations are often coupled with solvers for the electro-magnetic (EM) phenomena, resulting in sophisticated multi-physics simulation tools (see [6, 13]) that allow one to gather an insight of what is actually happening during the current interruption process - aspects that are very difficult to quantify with direct measurements for the above-mentioned reasons. Such simulation tools provide a significant added value to explain the results of experimental tests and to support the development of switchgear devices, however they also bring forth an important issue in addition to the inherent difficulty of plasma modeling: the need to find a good balance between the accuracy of the employed physical models and their computational complexity. The modeling of radiative heat transfer during the arcing process is an illuminating example of such an issue.

Radiation is one of the most important cooling mechanisms during switching, as it redistributes the heat produced by the current flowing through the plasma. Hence, accurate models of radiation are of fundamental importance to simulate the arc behavior, which is, due to the physics of radiation at the temperatures present in the plasma, a formidable task. Typically, the core of the arc is heated up to 25,000K, meaning that the electromagnetic radiation emitted by ions, atoms, and molecules of several different species (e.g. nitrogen, oxygen, or copper) have to be taken into account. The relevant window of the electromagnetic spectrum ranges from Hz-Hz, corresponding to wavelengths between m and m, i.e. from infrared to ultraviolet.
The main difficulty in simulating the electromagnetic radiation emitted by an arc derives from the complexity of the emission spectrum, where the relevant property, the absorption coefficient, changes by many orders of magnitude at spectral lines of which several 10,000 exist in the range under consideration. The propagation of radiative heat in space for each frequency is modeled, under assumptions that are reasonable for the arcing phenomena encountered in switchgear devices, by a first-order differential equation taking into account the absorption and the emission of radiation along the direction of propagation. The energy removed from the arc is with this defined by the temperature, pressure, and composition distribution within. Due to the complexity of the emission spectrum, a simple discretization according to frequencies leads to hundred thousands of very thin frequency bands; within each one of such bands the absorption coefficient can be assumed to be constant for fixed temperature, pressure and composition of the gaseous medium. This however, leads to the same number of three dimensional field equations which needed to be solved. Given the large number of finite volumes that have to be considered in CFD simulations of realistic geometries (see e.g. [22]), the use of such a large-scale radiation model is not feasible. Hence, there is the need to derive models for the radiative heat transfer with much lower complexity, possibly without compromising too much the accuracy as compared with the large-scale model. This issue has been tackled by several contributions in the literature [14, 23, 12, 20]. Most of the existing approaches consist in discretizing the fraction of the EM spectrum of interest into few bands, and assuming for each one of them some averaged absorption properties. The mentioned approaches have the advantage of being quite simple to implement, however in principle one should optimally choose ad-hoc different bands and averaged absorption coefficients as a function of pressure, temperature and chemical composition of the considered medium, since the radiation parameters are affected by all these aspects. If the bands and/or the averaged absorption coefficients are not chosen in an appropriate way, the model accuracy is worse and more bands are needed to improve it, resulting in a relatively large number of bands (6-10) in order to achieve accurate results with respect to the original, large-scale model. Hence, this procedure can be time consuming, not trivial to carry out in a systematic way, and ultimately suboptimal in terms of complexity/accuracy compromise.

In this paper, we present a new approach to derive small-scale, band-averaged models of the radiative heat transfer. We first describe the problem of radiation modeling from a novel perspective, where the aim is to approximate the input-output behavior of a large scale, linear-parameter-varying (LPV) dynamical system with that of a low-order one. The large scale system has one input (black-body intensity), one output (radiated intensity), three scheduling parameters (temperature, pressure and composition), and a large number of internal states (one for each considered frequency of the EM spectrum), while the low-order LPV system has the same input, output and scheduling parameters, but just a handful of internal states. From this point of view, the problem can be classified as a model-order reduction one [17]. Then, using classical tools for the analysis of signals and dynamical systems, we provide evidence that indeed models with quite low order (typically 2-3 bands) can be already good enough to capture the main behavior of the full-order model. Finally, we tackle the order reduction problem by using nonlinear system identification techniques (see e.g. [18]), where we collect input-output data from the large-scale system and use it to identify the parameters of the reduced-order model. The approach results in a nonlinear optimization problem (nonlinear program - NLP) with a smooth non-convex cost function and convex constraints, which are needed to preserve the physical consistency of the reduced-order model. We show through examples that the obtained reduced-order models enjoy a high accuracy with respect to the full-order one, while greatly reducing the computational times. As compared with the existing approaches, the method proposed here has the significant advantage of being systematic, i.e. there is no need to tailor or tune it for each different composition of the absorbing/emitting medium. The main user-defined parameter is the desired number of frequency bands in the reduced-order model, which can be then increased gradually until the desired tradeoff in terms of model quality vs. complexity is reached.

The paper is organized as follows. Section 1.1 provides a description of the problem we want to solve. In Section 2, such a problem is analyzed from a system’s perspective and connections are made to the order reduction of a large-scale LPV dynamical system. The proposed computational approach is described in Section 3, finally results are presented in Section 4 and conclusions and future developments are discussed in Section 5.

### 1.1 Problem formulation

Let us consider a region in space containing a hot gaseous mixture of different components. Each component is present in a fraction (e.g. of mass or of mole) . We indicate with the composition vector of the medium. We consider a line of propagation along which we want to compute the radiated heat, and denote with the position of a point lying on such a line. As the line crosses the hot matter, the temperature , pressure and composition change with (see Figure 1 for a graphical visualization).

The temperature, pressure and composition of the gas lie in some sets of interest, and respectively, defined as follows:

 T=[Tmin,Tmax]⊂R+P=[pmin,pmax]⊂R+Y={y∈[0,1]r:r∑i=1yi=1} (1)

Typical values defining the sets and are K, K, Pa and Pa.

Assuming without loss of generality that at the radiated heat intensity, for a given frequency of the EM spectrum, is equal to a given value , the distribution of the intensity as a function of can be computed through the following ordinary differential equation:

 dI(x,ν)dx=α(T(x),p(x),y(x),ν)(Ibb(T(x),ν)−I(x,ν)),I(0,ν)=Iν,0, (2)

where is the absorption coefficient and is the black body intensity, both pertaining to the frequency . For a fixed value of , depends on temperature, pressure and composition, while is a function of temperature only.

In equation (2), it is implicitly assumed that the radiation distribution has already converged to a steady state. Since the temporal dynamics of the radiated heat are much faster than the timescale of the phenomena of interest in arc simulations, such an assumption is reasonable. Similarly, scattering effects have been also neglected as they represent a negligible term for the application considered here. For a more complete form of equation (2) the interested reader is referred to [26].

Thus, for each frequency the corresponding ODE (2) is characterized by two parameters, namely the black body intensity and the absorption coefficient . The dependency of the latter on renders the equation nonlinear. From a physical perspective, the black body intensity is a source of radiation, while the current intensity is a sink: the change of radiated intensity in an infinitesimal space interval is the difference between these two contributions, scaled by the absorption coefficient. The black-body intensity as a function of the frequency and of temperature is given by Planck’s law [26]:

 Ibb(T,ν)=2hν3c21ehνkBT−1, (3)

where Js is the Planck constant, ms is the speed of light in vacuum and JK is the Boltzmann constant. The behavior of as a function of for some temperature values is shown in Figure 2.

For the purpose of this study, the total black-body intensity per unit of surface and solid angle (black-body radiance) is also needed, given by Stefan-Boltzmann’s law:

 ¯¯¯Ibb(T)≐∞∫0Ibb(T,ν)dν=σSBπT4 (4)

where WmK is the Stefan-Boltzmann constant. From equation (3) and Figure 2 it can be noted that, for the temperature values of interest for arc simulations in circuit breakers (3,000K-25,000K), most of the black-body radiance (i.e. the area below the solid curves in Figure 2) is contributed by the frequencies in the range of approximately Hz-Hz.

While the dependence of the black-body intensity on and is well-known from the above physical laws, the absorption coefficient is a much more uncertain quantity. Gaseous media have absorption spectra (i.e. the function relating to ) characterized by sharp lines at specific frequencies, whose values depend on the composition of the mixture and whose number is typically very large. Moreover, the absorption coefficient at each of such frequencies depends strongly on temperature and (less markedly) on pressure. Broadening effects due to pressure are also important as they contribute to a spread of the absorption lines over the nearby frequencies. Physical models to capture such complex absorption spectra have been proposed in the literature, [4, 16, 3, 8]. Moreover, databases of experimentally measured absorption data for several mixtures in specific frequency bands are available, typically for relatively low temperatures [21]. In the following, we will assume that is a known function of , in the domain of interest; we call such information the “base data” :

 D≐{α(T,p,y,ν);∀T∈T,∀p∈P,∀y∈Y,∀ν∈R+} (5)

As an example of information contained in , Figure 3 shows the absorption spectrum for a mixture of 50% silver, 25% air and 25% hydrogen at 16,300K and 10Pa. For a fixed frequency of such a spectrum, Figure 4 shows the corresponding absorption coefficient as a function of temperature and pressure (note the logarithmic scale for the absorption coefficient in this plot).

Given the base data, the total radiated heat intensity is computed by integrating the intensity contributed by each frequency of the EM spectrum:

 Itot(x)≐∞∫0I(x,ν)dν. (6)

In order to obtain a tractable problem, a first step is to adopt a fine discretization of the frequency domain such that (6) becomes a sum over a finite number of terms. Since the radiated intensity is relevant only in a specific region of the spectrum (compare Figure 2), the frequency discretization can be coarser outside such a region and finer inside, in particular around the peaks of the absorption spectrum. In this way, a finite number of frequency values are considered, each one being the middle point of an interval of the spectrum. Then, equation (6) can be re-written as:

 I(x)≃N∑i=0I(x,νi)Δνi. (7)

Equation (7), together with (2)-(3) and with the base data (5) evaluated at , form a high-dimensional model of the radiated heat intensity, named the Full-Order Model (FOM). The high dimensionality of the FOM comes from the fact that a large number of frequencies is taken into account in the discretization, so that the approximation error is small. Typical values of for the conditions of our interest are in the order of 1-2. Due to the large number of considered frequencies, the FOM can be used effectively only in very simple cases, for example in one-dimensional problems. In fact, in order to solve the radiative distribution in two- and three-dimensional cases, as it is needed for example in CFD simulations of real devices, one would have to solve many equations (whose nature depends on the chosen method, like the so-called P-N methods, the discrete ordinate method or the photohydrodynamic approach, see [26], [10]) for each one of the considered frequencies, leading quickly to an intractable problem. In this paper, we tackle this issue by deriving a method to compute low complexity models for the radiative heat transfer. More specifically, we consider the following problem:

###### Problem 1

: Given the full-order model defined by (2)-(3), (7) and the base data , derive a model with the same structure, i.e. where the total intensity is the sum of a finite number of contributions obtained by partitioning the frequency domain, but where the number of such partitions is very small, while still capturing accurately the total radiated heat intensity.

We call such a simplified model the Reduced Order Model (ROM). The ROM can then be effectively used to model the propagation of heat via radiation in many applications, including full three-dimensional simulations of plasma arcs encountered in switchgear devices.

We remark that, in light of Problem 1, in this work we will consider the radiated intensity computed with the FOM as “exact”. We will thus evaluate the quality of a given ROM by assessing the discrepancy between the radiative heat intensity given by the latter and the one given by the FOM. In other words, we don’t consider here the accuracy of the FOM with respect to the real-world behavior of gaseous media. Indeed, the topic of modeling accurately and/or measuring the absorption spectrum of a given gas as a function of temperature and pressure (i.e. to compute the base data ) is by itself an important and active research area [4, 16, 3, 8], however it is outside the scope of this work, which is focused on the simplification of the FOM into a computationally tractable model. On the other hand, the method presented in this paper does not depend on the specific absorption spectrum, i.e. it can be applied systematically to any base data, and it yields quite small discrepancies between the derived ROM and the employed FOM, such that the error between the ROM and the real behavior of the considered medium depends ultimately only on the quality of the FOM.

## 2 A system’s perspective of radiative heat transfer

### 2.1 Equivalent input-output models of the radiative heat transfer

As a preliminary step to address Problem 1, we re-write the FOM in a slightly different form, which is convenient to show that this model can be seen as a single-input, single-output (SISO) Linear Parameter-Varying dynamical system. Let us define the spectral emissivity as:

 e(T,ν)≐Ibb(T,ν)¯¯¯Ibb(T), (8)

i.e. the ratio between the total black body intensity and the one pertaining to each frequency of the EM spectrum. By construction we have and . The course of is illustrated in Figure 5.

By inserting the spectral emissivity (8) in equation (2), the radiated heat intensity for each frequency value of the FOM can be equivalently computed as:

 dI(x,νi)dx=α(T(x),p(x),y(x),νi)(e(T(x),νi)¯¯¯Ibb(T(x))−I(x,νi)),I(0,νi)=Iνi,0. (9)

Equation (9) is just a re-writing of (2) but, together with equation (7), it highlights the fact that the FOM can be seen as a dynamical system where is the independent variable, is an exogenous scalar input, are internal states, a scalar output, and are space-dependent parameters. In virtue of the fact that the function relating the temperature to the total black-body intensity is known, one can consider the latter as input to the system hence putting into evidence the Linear Parameter Varying (LPV) structure of the model:

 d¯¯¯I(x)dx=A(T(x),p(x),y(x))¯¯¯I(x)+B(T(x),p(x),y(x))¯¯¯Ibb(T(x))Itot=C¯¯¯I(x), (10)

where

 ¯¯¯I(x)≐[I(x,ν1),…,I(x,νN)]T∈RN×1A(T,p,y)≐diag([−α(T,p,y,ν1),…,−α(T,p,y,νN)]T)∈RN×NB(T,p,y)≐[α(T,p,y,ν1)e(T,ν1),…,α(T,p,y,νN)e(T,νN)]T∈RN×1C≐[Δν1,…,ΔνN]∈R1×N. (11)

A block-diagram of the FOM from this new point of view is shown in Figure 6(a).

Such a system’s perspective of the radiative heat transfer equations is the fundamental step at the basis of all the developments described in the remainder of this paper. We can now introduce more precisely the structure of a candidate reduced-order model (ROM) meant to approximate the FOM. Let us denote the bands of the ROM with , with (e.g. ), and the corresponding vector of inner states with . Then, we can write the equations describing the ROM as:

 d^¯¯¯I(x)dx=^A(T(x),p(x),y(x))^¯¯¯I(x)+^B(T(x),p(x),y(x))¯¯¯Ibb(T(x))^Itot(x)=^C^¯¯¯I(x), (12)

where

 ^A(T,p,y)≐diag([−^α(T,p,y,μ1),…,−^α(T,p,y,μN)]T)∈RM×M^B(T,p,y)≐[^α(T,p,y,μ1)^e(T,p,y,μ1),…,^α(T,p,y,μN)^e(T,p,y,νN)]T∈RM×1^C≐[1,…,1]∈R1×M. (13)

Figure 6(b) gives a graphical representation of the ROM. Practically speaking, from (12)-(13) one can see that each component of vector accounts for a certain portion of the total radiated heat intensity, in complete analogy with the FOM (10), the only difference being the number of frequency bands, which in the ROM is much smaller than in the FOM. For each band , the parameters and have thus the meaning of “equivalent” absorption coefficient and spectral emissivity, respectively. The task of deriving a ROM for the radiated heat transfer from the FOM is referred to as model order reduction. Computing a ROM is equivalent to assigning suitable values to and as a function of the underlying parameters and . In particular the collection of all bands has to form a non-overlapping partition covering the whole EM spectrum. Then, the value of has to correspond to the integral of the spectral black body intensity (3) over the frequency band pertaining to . Hence, one can equivalently state that computing a ROM amounts to choose, for each pair of pressure and composition values, a partition of the EM spectrum (which defines the equivalent emissivity as a function of temperature) and the courses of the corresponding equivalent absorption coefficients as a function of temperature. In previous contributions in the literature, e.g. [20, 23], the task of defining the ROM has been carried out by picking a finite number of bands covering the EM spectrum and then computing the equivalent absorption coefficients through some averaging procedure on the portion of the absorption spectrum contained in each band. This approach is simple to implement but it has the drawback of not being systematic, since both the choice of the band cuts and the averaging of the absorption coefficients have to be made by the user, without an immediate link to the accuracy of the resulting ROM. In the next sections, we will present a new order reduction approach to compute the partitioning of the EM spectrum and the values of in a systematic way, that yields quite accurate results as compared with the FOM.

About this last point, i.e. the accuracy of the ROM, we note that the input of the FOM and of the ROM is exactly the same, corresponding to the total black-body radiation for the considered temperature profile, . Therefore, it is quite intuitive that the discrepancy between the total intensity given by the FOM, , and the one predicted by the ROM, , for the same temperature profile (i.e. the same distribution of ), represents a reasonable indicator of the accuracy of the reduced order model. In other words, the error signal will be considered to evaluate the goodness of a given ROM. This choice is motivated by the fact that the total radiated heat intensity is the main quantity of interest predicted by the ROM when it is embedded in multi-physics simulations of arc plasma, since it is used to compute the heat transferred from the plasma volume to the walls, and also (through its divergence) the heat redistributed within the plasma volume. Hence, the ROM should reproduce this quantity as accurately as possible with respect to the FOM, given the same spatial distribution of temperature, pressure and chemical composition.

Before going to the details of the proposed method to derive the ROM, a sensible question to be addressed is whether the approximation problem we are dealing with has a reasonably good solution or not. More specifically, recall that we aim to approximate the input-output behavior of a large scale system, with hundred of thousands of internal states, with that of a small-scale one, with at most a handful of internal states. It is not immediately clear if there exist such a low-order ROM still capable of delivering high approximation accuracy, since this aspect depends on the characteristics of the FOM, i.e. on the underlying physics of the radiative heat transfer. In the next section, we exploit the system’s perspective described above to provide an intuition that indeed a ROM with a handful of bands can capture most of the input-output behavior of the FOM.

### 2.2 Frequency-domain analysis

Analyzing the complexity of the FOM in its form (10) in the domain of position is not straightforward due to the very large number of internal states, each one following its own dynamic evolution. Besides noticing that the FOM is given by the sum of a large number of non-interacting, asymptotically stable first-order systems, all driven by the same input and whose (position dependent) poles are given by , little less can be said.

However, if we consider fixed values of temperature, pressure and composition, and respectively, and we assume that only infinitesimal perturbations of temperature take place in space, such that the absorption properties of the medium can be assumed constant, we immediately notice that the FOM becomes a linear-parameter-invariant (LPI) system, to which well-assessed tools in systems theory and signal processing can be applied. In particular, after establishing the analogy between the position in the FOM with the continuous time variable in dynamical systems, we can study the input-output response of the FOM to such infinitesimal temperature variations by applying the Laplace transform [25] to its equations and deriving the transfer function from its input to its output , where is the Laplace variable:

 G(s)≐Itot(s)¯¯¯Ibb(s)=C(sI−A(¯¯¯¯T,¯¯¯p,¯¯¯y))−1B(¯¯¯¯T,¯¯¯p,¯¯¯y), (14)

In (14), are the matrices given in equation (11) and evaluated at the chosen temperature, pressure and composition values, I is the identity matrix of suitable order and denotes a matrix inverse operation. A tool commonly used to analyze the transfer function of a dynamical system is the Bode diagram of the corresponding frequency response, obtained by evaluating the magnitude and phase of , where assumes here the physical meaning of the frequency of purely sinusoidal oscillations (in space) of the black-body intensity (i.e. of temperature), with infinitesimal amplitude, with the period measured in (e.g. rad/m corresponds to a period of oscillation of the input of roughly m). As an example, the Bode diagram of the FOM frequency response obtained by fixing K, Pa and containing 50% silver, 25% air and 25% hydrogen (whose absorption spectrum is shown in Figure 3) is shown in Figure 7.

It can be noted that the overall behavior of the FOM is that of a low-pass filter, whose frequency response is shaped by the contributions of the large number of internal states present in the system. As a quantitative example regarding this case, a sinusoidal oscillation of the black-body intensity (due to an oscillation in temperature) with a period of 2m would give rise, in the gaseous medium, to a sinusoidal distribution of radiated heat intensity with unchanged amplitude and phase with respect to what would happen in vacuum (compare Figure 7 with rad/m), i.e. the amplitude of the oscillations would be equal to that of the input black-body intensity, and the spatial distribution would show almost zero phase shift. On the other hand, the spatial distribution of radiated heat induced by a sinusoidal oscillation of the black-body intensity with a period of 2m (compare again Figure 7, with rad/m) would have an amplitude equal to only about with respect to that of the input, with a phase lag of about deg.

If we consider now the order reduction of such a dynamical system, we see that this is a standard problem of model order reduction of LPI systems, for which a well-established literature exist, see e.g. [17, 1]. Thus, we can use one of the existing approaches to derive a reduced order model that approximates the FOM for the chosen values of and . After deriving the ROM, the related transfer function can be computed as (compare equation (13)):

 ^G(s)≐^Itot(s)¯¯¯Ibb(s)=^C(sI−^A(¯¯¯¯T,¯¯¯p,¯¯¯y))−1^B(¯¯¯¯T,¯¯¯p,¯¯¯y). (15)

A comparison between the Bode diagram of the latter and that of the FOM reveals that up to very small gain values of about , such that the corresponding radiated heat intensity is negligible, the frequency responses of the two models are practically super-imposed, hence showing a very good agreement between them. Moreover, applying this procedure for many values of and , chosen by gridding their respective domains and , shows that such a good agreement is obtained always with no more than four-five bands in the ROM. A good agreement up to a gain of about is obtained with just two bands in the ROM, for oscillation periods of fractions of millimeters. Indeed, this level of accuracy would be enough for the sake of arc plasma CFD simulations, where the resolution of the spatial discretization of the considered volume is of the order of m. An example of the obtained results is depicted in Figure 7, too. In particular, the ROM whose frequency response is shown in the figure has five bands, , and the corresponding absorption coefficients and emissivities (and frequency boundaries) are equal to

 ^α(¯¯¯¯T,¯¯¯p,¯¯¯y,μ1)=1.110−1m−1^e(¯¯¯¯T,¯¯¯p,¯¯¯y,μ1)=7.210−1[0,1.61015]Hz^α(¯¯¯¯T,¯¯¯p,¯¯¯y,μ2)=9.210−1m−1^e(¯¯¯¯T,¯¯¯p,¯¯¯y,μ2)=2.110−1(1.61015,2.41015]Hz^α(¯¯¯¯T,¯¯¯p,¯¯¯y,μ3)=1.3101m−1^e(¯¯¯¯T,¯¯¯p,¯¯¯y,μ3)=5.410−2(2.41015,3.11015]Hz^α(¯¯¯¯T,¯¯¯p,¯¯¯y,μ4)=1.2102m−1^e(¯¯¯¯T,¯¯¯p,¯¯¯y,μ4)=1.210−2(3.11015,3.81015]Hz^α(¯¯¯¯T,¯¯¯p,¯¯¯y,μ5)=1.7103m−1^e(¯¯¯¯T,¯¯¯p,¯¯¯y,μ5)=3.910−3(3.81015,+∞)Hz

A comparison between these values and Figure 7 shows that the values of the ROM absorption coefficients correspond to the dominant poles of the FOM, as the intuition would suggest.

Overall, the analysis reported so far provides if not a rigorous proof at least an indication that the problem we are dealing with has a reasonable solution, using ROMs of quite low order. In section 3, we describe in details the solution approach that we propose to deal with the linear-parameter-varying case.

### 2.3 Position discretization

Before proceeding further, it is convenient to introduce the discretized versions of the FOM and of the ROM, where the position variable is taken at nodes , that are equally spaced by an interval . The latter has to be chosen according to the features of the problem at hand, trading off computational speed with a sufficiently fine discretization, which can capture well the fastest transients of the model’s input and output. As a rule of thumb, can be chosen, where is the smallest space-scale of interest in the problem (typically in our case m, as mentioned above). Hence, all the space-dependent variables (i.e. ) are now evaluated at discrete position values. The discretization of the radiation models is carried out by assuming that such variables are constant between two subsequent position nodes, and , and then computing the explicit integration of (10) (for the FOM) and (12) (for the ROM). In particular, for the full-order model we have:

where the matrices and are computed as:

and

 a(T,p,y,νi)=e−α(T,p,y,νi)Δx,i=1,…,M.

The matrix in (16) is the same as in (10), since the output equation is static and thus it is not changed by the discretization of the position . Similarly, for the reduced-order model we have:

where the matrices and are computed as:

and

 ^a(T,p,y,μi)=e−^α(T,p,y,μi)Δx,i=1,…,M. (20)

Also for the ROM the matrix in (18) is the same as in (12).

The reason why we employ such a space discretization is twofold: on the one hand, it is needed to obtain a finite-dimensional computational problem, on the other hand it improves the computational efficiency of the method (in particular by using a constant discretization step ).

## 3 Model order reduction of the radiative heat transfer equation

### 3.1 Solution approach

Considering the analysis of Section 2.2, one can be tempted to derive the ROM, i.e. the functions and , by gridding the domains and and for each triplet compute the equivalent absorption coefficient and emissivity of the corresponding LPI model, using well-assessed and efficient model order reduction techniques. Then, the ROM can be obtained by interpolating among the computed reduced-order LPI models. This approach could work well if only pressure and composition dependence were considered, since the sensitivity of the base data on these values is mild and hence one can be confident that interpolating among the bands computed at different triplets yields correct results. In other words, for a given temperature , the absorption spectrum of the FOM does not change dramatically between two neighboring pairs and within the pressure and composition intervals of interest for the application considered here, so that for each band of the ROM it is safe to interpolate between the coefficients and , computed independently by means of LPI model order reduction, to obtain the ROM coefficients for generic values of with . However, this procedure would not achieve good results when temperature dependence is taken into account as well. In fact, the dependence of the absorption spectrum on temperature is very strong (compare Figure 4), so that the ROM coefficients pertaining to the same band (e.g. ) but computed at two different temperature values, even with the same pressure and composition values, might be completely unrelated to each other, and interpolating between them can give highly inaccurate results. Driven by these considerations, we adopt a hybrid strategy, where we grid the domains and and for each pair we derive the partition of the EM spectrum in bands and the related functions that define the ROM. Since the corresponding FOM is now parameter varying (because we let the temperature distribution change while fixing only pressure and composition), this problem falls in the class of model-order reduction of LPV systems, for which, differently from the LPI case, few results exist in the literature [28, 24], and their practical applicability to systems with states, like the FOM in our problem, is not straightforward. For the above reasoning, in the remainder of this section it is assumed, unless otherwise stated, that a fixed pair of pressure and composition values has been chosen, and that the only space-varying variable is the temperature . The complete ROM can be then obtained by repeating the LPV order-reduction for all the pairs chosen by gridding the respective domains, and then interpolating among the obtained values to compute the equivalent absorption coefficients and frequency bands for a generic triplet .

###### Remark 1

The simplification introduced by fixing pressure and composition is made possible thanks to the above-discussed particular properties of the absorption spectra of the considered gaseous media. Indeed, such a simplification by itself is not required for the order-reduction technique described in the following, which may straightforwardly be applied also with varying pressure and composition, but at the price of higher computational requirements. In our experience, using constant pressure and composition in the order reduction computation yields accurate enough results for the applications of interest.

In the following, we propose to address the LPV order reduction problem with a nonlinear system identification approach (see e.g. [18]), in which we search, within a given set of possible ROMs, the one which is closest to the FOM according to a pre-defined optimality criterion. In the next sections, this task is brought in the form of a tractable optimization program. It has to be noted at this point that, in the literature on system identification, there exist several contributions devoted to the problem of identification of LPV systems, see e.g. [5] and the references therein. However, such results are not applicable in our case, due to the additional constraints that are present on the ROM, namely the need to preserve a specific structure where the total radiated intensity is the sum of the contributions given by the frequency bands. Differently from such previous approaches, the one proposed here is able to take into account these constraints, since it is specifically tailored for the considered application.

#### 3.1.1 Cost function and model set

In order to have an optimization problem that can be solved with common numerical techniques, two main ingredients need to be defined: the set of reduced-order models of the form (18) where we carry out our search, denoted with (model set), and a cost function giving a measure of how much a given ROM is close to the FOM. The model set should represent the limitations that we want to impose on the ROM in order to account for the physics of the problem. A ROM is fully characterized by its parameters , which in the considered settings are, for each band , functions of temperature only. To be consistent with the underlying physical phenomena, the equivalent absorption coefficients should be positive (to retain an asymptotically stable model):

 ^α(T,¯¯¯p,¯¯¯y,μi)>0,∀T∈T,i=1,…,M (21)

and the equivalent emissivities should lie in the interval and sum to one over all the considered temperature range, in analogy with (8), so that the black-body limit is not violated:

 0≤^e(T,¯¯¯p,¯¯¯y,μi)≤1,∀T∈T,i=1,…,MM∑i=1^e(T,¯¯¯p,¯¯¯y,μi)=1,∀T∈T (22)

Since equation (20) is invertible, we can equivalently consider the functions to define the model set. We select this alternative for the sake of computational efficiency, as we will discuss more in details in section 3.2.3. Then, the constraints (21) can be re-written as:

 0≤^a(T,¯¯¯p,¯¯¯y,μi)<1,∀T∈T,i=1,…,M (23)

As a final step to define , we choose a finite parametrization of functions and , in order to obtain a finite dimensional model set (hence also a finite dimensional optimization problem):

 ^a(T,¯¯¯p,¯¯¯y,μ)=f^a(T,θ^a(¯¯¯p,¯¯¯y,μ)),θ^a(¯¯¯p,¯¯¯y,μ)∈Rnθ^a^e(T,¯¯¯p,¯¯¯y,μ)=f^e(T,θ^e(¯¯¯p,¯¯¯y,μ)),θ^e(¯¯¯p,¯¯¯y,μ)∈Rnθ^e, (24)

where are chosen, once again, to tradeoff the flexibility of the parametrization with computational efficiency. In particular, a convenient choice for is the class of piecewise affine functions of temperature, while is taken such that the corresponding parameters define the partitioning of the EM spectrum in a finite number of bands. These choices have the advantage of yielding a convex model set (see section 3.2.2 for details), hence improving the computational efficiency and stability of the procedure. For a given pair , let us collect the parameters and for all the ROM bands in a single vector :

 θ≐⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣θ^a(¯¯¯p,¯¯¯y,μ1)⋮θ^a(¯¯¯p,¯¯¯y,μM)θ^e(¯¯¯p,¯¯¯y,μ1)⋮θ^e(¯¯¯p,¯¯¯y,μM)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦∈Rnθ, (25)

where is the total number of optimization variables. Then, considering (22)-(25), we can define the model set as

 H≐⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩θ∈Rnθ:0≤f^a(T,θ^a(¯¯¯p,¯¯¯y,μi))<1,i=1,…,M0≤f^e(T,θ^e(¯¯¯p,¯¯¯y,μi))≤1,∀T∈T,i=1,…,MM∑i=1f^e(T,θ^e(¯¯¯p,¯¯¯y,μi))=1,∀T∈T⎫⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪⎭ (26)

As to the cost function , this should represent a criterion by which the accuracy of the ROM is evaluated. As discussed in section 2, an indicator of the accuracy of a ROM, suitable for the considered application, is related to the error between its total output intensity profile and that of the FOM, given the same temperature profile. Motivated by this consideration, we select a series of temperature profiles, and then evaluate the discrepancy between the corresponding outputs generated by the FOM and those given by the ROM. More specifically, let us consider a space interval , chosen such that for some . Similarly to the choice of , the value of depends on the problem at hand: a typical choice is three-four times the largest space-scale of interest (see section 3.2.1 for details). Thus, in the interval we have position nodes , with and . Consider now a finite sequence of discretized temperature profiles . In face of each of such temperature profiles, the FOM provides a corresponding profile of the radiated intensity, , and similarly the ROM, for a given value of provides an approximated one, . Let us define the weighted error profile as

 ΔIj(θ)≐wj(~Itot,j(~Tj)−~^Itot,j(~Tj,θ)),j=1,…,LT, (27)

where is a scalar weight (the specific choice of is discussed in section 3.2.1). Then, we take the cost function as:

 J(θ)=LT∑j=1ΔITj(θ)ΔIj(θ), (28)

i.e. the sum, over all error profiles and over all position nodes of each profile, of the weighted intensity error squared.

#### 3.1.2 Nonlinear Program formulation

We can now write the model order reduction problem in a computationally tractable form:

 minθ∈HJ(θ) (29)

where and are given in (26) and (28), respectively. Namely, the aim of the optimization problem (29) is to search, within the model set , the value of the parameter vector that achieves the best fitting criterion . Problem (29) is a finite-dimensional nonlinear program (NLP) which, for a suitable choice of the model parametrization (24), can be tackled with state-of-the-art numerical methods [19]. As we show in the examples of section 4, a typical dimension of the optimization variable is , corresponding for example to a 3-order model with parameters defining the functions and two parameters defining the frequency cuts in the EM spectrum that separate the three bands . In general, the cost function (28) is non-convex with respect to the decision variable , so that only a local solution to problem (29) can be computed efficiently with such problem dimensions. As we show in the examples of section (4), the obtained solutions nevertheless yield ROMs that are very accurate with respect to the FOM.
In the next section, we discuss several aspects involved in the formulation and solution of (29), including the design of the input temperature profiles and the choice of the weights in the cost function (28), the choice of model parametrization, the numerical approach to solve the NLP, finally some ways to improve the efficiency and stability of the numerical optimization.

### 3.2 Computational aspects

#### 3.2.1 Input design, weighting coefficients and initial conditions

From an application’s perspective, the position interval , the resolution of the position discretization , and the temperature profiles should be representative of the typical temperature distributions in the plasma that is generated during the switching process. About the choice of , considering the discussion in section 2.3, a reasonable tradeoff between accuracy and computational speed is m, while for a good choice in our experience is m, for arcs whose width is between m and m, corresponding to current values in the range 1-40A. Finally, for the temperature profiles, we select a bell-shaped function, where we can adjust the maximum temperature reached as well as the steepness of the rising and falling slopes. Then, we generate a series of such profiles by cycling through different maximum temperatures and transient slopes. Furthermore, from a system identification’s perspective, given the nonlinearity of the system’s equations, different initial conditions for the internal states of the model should be considered, and the input should be designed in order to excite the system in a broad range of frequencies. For the former aspect, we replicate the temperature profile several times, so that the system is presented more than once with the same temperature profile, but each time starting from the internal states resulting from the previous profile. For the second aspect, white noise processes are a well-known choice [15] to excite the system’s dynamics over all frequencies, hence we super-impose to the computed profiles a uniformly distributed white noise temperature signal whose amplitude is a fraction (e.g. ) of the highest temperature in the original profile. With these choices, typical temperature profiles are shown in Figure 8(a)-(b), together with the resulting intensity distributions computed with the FOM for a mixture of pure air.

Due to the fourth-order dependence of the total black-body intensity on temperature (see (4)), the radiated intensity obtained with different temperature profiles can have different orders of magnitude, compare for example Figure 8(a) with Figure 8(b), where the maximum temperatures are about 10K and K, respectively. If no corrective measure is taken, such a difference among the various intensity profiles would lead to poor accuracy of the reduced order models, due to the biasing towards the higher intensity values in the cost function. In order to compensate this effect, we select the weights in (27) on the basis of the radiated heat intensity given by the FOM when the corresponding temperature profile, , is considered. In this way, all the error profiles (27) are normalized with respect to the corresponding intensity; typical specific choices for the weights include the maximum or the average radiated heat, i.e.:

 wj=max(~Itot,j),j=1,…,LTorwj=mean(~Itot,j),j=1,…,LT

Finally, it is worth mentioning the choice of the initial conditions for the radiated heat intensity, for both the FOM and the ROM, i.e. the values of vectors and needed to simulate the models (16) and (18), respectively, in the interval . To this end, we assume that for (i.e. outside such an interval) the incoming radiation in positive direction corresponds to the steady-state intensity at ambient temperature , i.e. for the FOM and

 ^Ij(x,μi)=^e(Ta,¯¯¯p,¯¯¯y,μi)¯¯¯Ibb(Ta),i=1,…,M (30)

for the ROM. The rationale behind this choice is that in the arc simulations related to switchgear devices one can assume that initially the boundaries of the computational domain are at ambient temperature and radiate the corresponding intensity into the volume filled with plasma. Additionally, the boundaries can not reach temperatures higher than about 3,000 K (depending on the material), such that the related radiated heat is anyway negligible with respect to the intensity emitted in the plasma volume where current is flowing.

#### 3.2.2 Model parametrization and model set

The choice of the model parametrization, i.e. of functions and in (24), is crucially important for the accuracy of the obtained results and for the solution of the numerical optimization. For the first aspect, one shall choose a rich enough family of functions, such that the data from the FOM can be reproduced with small errors. For the second aspect, the best choice would be a parametrization leading to a convex model set (26), so that sequential quadratic programming and line search algorithms [19] can be used efficiently. For the equivalent absorption coefficients , i.e. function , a choice that meets both requirements is a piecewise affine parametrization, where is, for each band , a vector containing the values of the coefficients at a finite number of pre-defined temperature nodes , such that , chosen by the user (e.g. equally spaced):

The temperature nodes must include the values at the boundaries (1) of the domain , i.e. and .
Then, for a given temperature value , the function in (24) is computed by interpolating linearly among the values of corresponding to the neighboring temperature nodes:

 f^a(T,θ^a(¯¯¯p,¯¯¯y,μi))=λ(T)Tθ^a(¯¯¯p,¯