Comparison of Two Models for Radiative Heat Transfer in High Temperature Thermal Plasmas

Numerical simulation of the arc-flow interaction in high-voltage circuit breakers requires a radiation model capable of handling high-temperature participating thermal plasmas. The modeling of the radiative transfer plays a critical role in the overall accuracy of such CFD simulations. As a result of the increase of computational power, CPU intensive methods based on the radiative transfer equation, leading to more accurate results, are now becoming attractive alternatives to current approximate models. In this paper, the predictive capabilities of the finite volume method (RTE-FVM) and the P1 model are investigated. A systematic comparison between these two models and analytical solutions are presented for a variety of relevant test cases. Two implementations of each approach are compared, and a critical evaluation is presented.


Introduction
The design of modern high-voltage circuit breakers involves the simulation of high temperature thermal plasmas interacting with fluid flows in complex configurations.In such devices, an intense radiation from an electric arc is at the origin of coupled energy exchange mechanisms with a participating medium and confining boundaries.Among the critical issues in modeling such coupled phenomena is the accurate calculation of the radiative energy transfer within the gaseous medium and the fraction reaching the boundaries, causing them to ablate.In industrial CFD simulations, this results in a compromise between computational resources and modeling fidelity.Considering that radiation is the dominant mode of energy transfer, the predictive capabilities of two different radiation heat transfer approaches are investigated.These are the radiative transfer equation (RTE) which models the energy transfer as an integrodifferential equation, and the P1 model, an approximation resulting in a simpler Helmholtz type of equation.
The comparison of the methods being investigated is carried out by a series of test cases which are representative of the difficulties encountered in circuit breaker configurations.In increasing order of complexity, these difficulties are (1) the disparity of scale between the radiating arc and the enclosure, (2) the radiative transfer over long distances in a participating medium, (3) the discretisation of directional space in the presence of preferred directions in the arc geometry, and (4) the discontinuous optical properties between the arc and the surrounding medium.These are embodied in three cases computed with two implementations, the present code, MC 3 , and that of the commercial package, fluent for both models.With the availability of analytical solutions [1,2], a quantitative comparison of the methods is possible.
The goal is to illustrate concretely through these systematic comparisons some weaknesses of the P1 model and the possible improvement obtained with the RTE-finite volume method.

Radiative Heat Transfer
2.1.Governing Equations.Following [3], the conservation of the radiative energy along a ray path of direction s can be written as: where I ν is the spectral radiative intensity in a direction s at a frequency ν. κ ν is the coefficient of absorption and σ sν is the coefficient of scattering.β ν = κ ν + σ sν is known as the extinction coefficient.Finally, Φ( s i , s) is the phase function.
The difficulty in numerical solutions of (1) is the large number of variables for the radiative intensity I ν : three for space r, two for direction s, and one for the frequency ν.In the present context, certain simplifications are possible.Scattering can be neglected because the size of the particles are much smaller than the wave lengths present.Following [4,5], this is equivalent to let σ sν = 0, hence the last term of (1) vanishes and β ν is replaced by κ ν , giving The frequency dependence is handled by splitting the radiative spectrum into gray bands (typically five) as proposed by [6].
While inexpensive approaches such as the net emission coefficients [7] can be applied for the arc energy balance, these simplified 1D models do not provide the mechanism for the interaction between the arc and the surrounding medium and boundaries.In order to obtain a self-consistent CFD-type modeling, a more elaborate radiation model is needed.

The P1 Model.
Because of the high computing costs associated with the solution of the radiative transfer equation (RTE), affordable implementations resulting in a substantial simplification have been devised ( [3,8]).In the P1 model, it is assumed that the radiative intensity can be expressed as a Fourier series which is variable separated into a function depending on coordinates and another on the direction as follows: where I m l (r) corresponds to the position coefficients and Y m l ( s) to the spherical harmonics.
Taking the first term of this series, [3] shows that the radiative intensity can be related to G ν , the incident radiation, which leads to the P1 model for the incident radiation G ν : A major concern in practical applications of this model is the lack of a physical realistic set of boundary conditions.Following established practice in the literature, a zero flux can be imposed on the symmetry axis, and for the nozzle walls, the most appropriate physical condition is the Marshak condition which is expressed as In the present work, w = 1 is assumed as the walls can be treated as nonreflecting black bodies.The Helmholtz type governing (4) is solved by a finite volume discretization [8], where the incident radiation G is solved at the nodes of a triangular mesh, whereas the values of κ and the black body emission I b are cell centered.Assuming a linear variation of G inside a triangle, and using the Gaussian integration, this variable storage ensures the continuity of G, whereas I b and κ can be discontinuous.

The RTE Finite Volume Method (RTE-FVM).
Another approach is to solve (2) directly using a formal finite volume discretization of both the spatial domain and the directions.Initially proposed by Briggs et al. [9] from the field of neutron transport, this idea was adapted to the radiative heat transfer by pioneers like Raithby and Chui [10] and Chai et al. [11].
To reduce the number of unknowns, an axisymmetric scheme is devised as described in [12].The discretisation is obtained by sweeping the control volumes from an unstructured two-dimensional mesh through an angle ΔΦ, and the 2π solid angle is divided into N θ × N ψ control solid angles equally distributed by Δθ and Δψ.

Results for Analytical Test Cases
Three test cases for which analytical solutions are available will be presented.For each case, both the P1 model and the FVM model are applied, and for each model, the present in house implementation, MC 3 , as well as the fluent implementation will be reported.The main purpose of this comparison of two implementations of the same radiation model is to validate the present implementation and the conclusions.The small difference between MC 3 and fluent is attributable to variations in the grids and possible differences in accounting for boundary conditions.

Test Case I.
One important characteristic of circuit breaker configurations is the difference between the arc and the enclosure dimensions.Thus, to assess the capability of a model to predict long distance propagation, a test case is proposed in Figure 1 which approximates a circuit breaker configuration with an electric arc on the axis, radiating towards an outer wall through a layer of absorbing cold gas which emits no radiation.The enclosure is a cylinder of radius R of 1 m and a length of 2 m filled with a gray gas at uniform temperature T g with an absorption coefficient κ, from which the nondimensional optical length is defined  as τ = κL = κR.The computational grid consists of around 1500 control volumes.Three different values of the coefficient τ are used to represent low (τ = 0.1), medium (τ = 1.0), and high (τ = 5.0) optical lengths.A symmetry condition is imposed on the axis and the temperature is set at T = 0 on the other three walls.The high temperature region at T g is confined to the region inside the radius of r = 0.1.
Figure 2 shows the comparison of the normalized incident radiation, Q * , at the upper wall for the various methods compared with the analytical solution for three values (0.1, 1.0, 5.0) of the optical length τ.A first comment concerns the absolute value of the radiative heat flux arriving on the upper boundary.This value does not vary monotonically with respect to the absorption coefficient, and the maximum incident flux corresponds to the medium optical length of τ = 1.This behavior results from two opposing effects: the stronger emission of the gas inside the arc region which is proportional to κ and the stronger absorption of the cold region which has an exponential behavior.
As seen from Figure 2(a) which presents the results for τ = 0.1, the P1 model misses the axial dependency of the flux and predicts an almost constant value.For τ = 1 and τ = 5, the values predicted by P1 are, respectively, 20% higher and 400% lower than the analytical results.The RTE-FVM consistently shows a very good agreement with the analytical case and represents a significative improvement over the P1 model for all values of τ.Both RTE-FVM methods tested show a stepwise behavior for all values of τ, which can be attributed to a form of "ray effect".This is a well-known characteristic of the FVM [13,14], and results in a progressive alignment of the radiation intensities with the discrete directions used when a long distance propagation is present.These two implementations use a similar first order discretization and slightly diffuse but never produce nonphysical results.

Test Case II.
The purpose of the following test case is to analyse the behavior of the directional discretization in the presence of a privileged direction in the geometry of the arc.This is achieved by devising a conical temperature distribution, as shown in Figure 3.
The optical length has been arbitrarily set at τ = 1; the hot region is at 1000 K and the cold at 0 K.The angle of the cone is 26.5 • and the normal direction (116.5 • ) is that along which the most energy is transferred.The numerical results are carried out on a uniform triangular grid composed of 4072 elements and compared with the analytical solution of [15].Figure 4 shows the normalized incident radiation at the upper wall for two sets of directional discretizations, 2 × 2 and 10 × 10.The large difference between these two sets of results illustrates the critical importance of this parameter on the solution.Table 1 summarizes the results of a parametric study carried out to investigate the dependency of the solution on this aspect of the discretization.The first two columns give the number of directions used and the corresponding values of Δθ for the test.The column with the heading θ i gives the closest direction present in the set of discrete directions, relative to the direction of 116.5 • .The column headed |θ i − 116.5| gives the difference between the closest discrete direction and that of 116.5 • .
The accuracy of the results is evaluated based on the average error defined as the average relative distance between equally spaced points on the analytical solution and the FVM solution and is given in the last column of Table 1.First, one can see that FVM solutions converge to the analytical solution with the number of directions.Second, it can be noted that even if the direction at 116.5 • almost coincides with one of the discrete directions in the sets of 7×7 and 10× 10 directions, the average error varies almost monotonically with the number of directions and is not dependent on the value of |θ i − 116.5|.Thus, basically, the more directions  are used in the discretization, the more accurate the results, despite the presence of a privileged direction.

Test Case III.
This last analytical test case, shown in Figure 5, has been devised to reflect typical dimensions as well as the order of magnitude of the coefficient of absorption κ encountered in a circuit breaker chamber.The radiation takes place inside a cylinder, 20 mm long with a 10 mm radius, divided into two zones: a hot internal core from radius 0 to 5 mm, representing a plasma at 10,000 K, surrounded by a cold layer at 0 K, from radius 5 mm to 10 mm.The walls are cold (0 K) black bodies and two sets of values of κ, one for each zone, are defined to mimic radiation transport in circuit breakers.Indeed, the spectral behavior of SF 6 (used in the circuit breakers) can be roughly split into  two types of bands: (i) The "opaque band" (UV): typically κ = 10 m −1 in the plasma and κ = 10 5 m −1 in the cold zone.
The comparison of the two radiation models is carried out by considering two different results: the divergence of the radiative heat flux in the midplane, which provides the source terms in the energy conservation equations needed to achieve the CFD coupling, and the incident radiative heat flux at the radial wall, which will cause nozzle ablation.It has been shown [3] that, for all points in the domain, the divergence of the radiative heat flux can be evaluated as ∇ • Q rad = κ(G − 4πI b ).Again, the analytical results are computed by following the method of [2] with a slight modification to account for a piecewise constant coefficient of absorption κ.
The results for the divergence of the radiative heat flux (the source term) are presented in Figure 6 for the opaque band.Because the value of κ is very high in the cold region for this band, all the radiation emitted in the hot region is absorbed within a tiny strip next to the interface between the two zones.The two implementations of the RTE-FVM exhibit a diffuse behavior at the interface, namely that the radiation is absorbed over a larger region than expected which is very dependent on the mesh used.Thus, the magnitude of the peak at the interface is greatly influenced by the mesh density close to the interface.However, since the method is conservative, the same quantity of energy is absorbed regardless of the size of the region of absorption.The P1 model shows a similar behavior but requires a very fine mesh (up to 100,000 triangles whereas 8000 and 10 × 10 directions are sufficient for the FVM) to capture correctly the order of magnitude of G, the incident radiation.Nevertheless, the magnitude of I b largely dominates G and, as a consequence, the P1 model manages to predict the source term quite accurately even if G is ill predicted.In a practical configuration there would not be such a discontinuity between the two zones, and this test case represents the worst scenario that a simulation code can encounter.The RTE-FVM presented a stable and trustworthy physical behavior, while the P1 model predicted nonsymmetrical distribution of the source term for mesh sizes coarser than 50,000 triangles.One possible explanation is proposed: at the interface, the P1 model assumption concerning the isotropy of the radiation is clearly wrong and leads to ill predictions.Finally, because none of the radiation reaches the radial wall for the opaque band, there is no point in presenting the incident radiative heat flux at the wall.
The second set of results for the transparent band, are presented in Figure 7 for the source term and Figure 8 for the incident flux at the radial wall.The prediction of the source term does not seem to be a challenge for this case, because both the RTE-FVM and P1 models accurately compute its magnitude whatever the method or the implementation used.As for the opaque band, the black body emission I b largely dominates G, thus even if G is wrong, it does not significantly affect the overall result.Moreover, for this case the P1 model does not suffer from a strong dependency on the mesh size and requires fewer triangles to obtain a symmetrical solution.A mesh of 8000 triangles has been used for both methods (and 10 × 10 directions for the FVM).The results are different for the radiative flux at the wall, as shown in Figure 8.The P1 model fails to predict correctly the 2D effects of the radiative heat flux reaching the radial wall, whereas the RTE-FVM succeeds equally well as the previous tests.
This test case is particularly enlightening as it explains why the P1 model has been successfully used in the field of circuit breakers: it predicts with a sufficient accuracy, at a low CPU cost, the source term which allows the CFD coupling in the energy equation and gives the correct order of magnitude for the flux on the walls.However, the present analysis demonstrates that the mesh has to be extremely fine in the limiting case of opaque bands with strong gradients presented by the different physical quantities.Moreover, the  P1 model is not suitable for predicting the 2D effects on the radiative flux reaching the walls, which will hamper an accurate prediction of the ablation process.As a consequence, ablation will be predicted at wrong locations, affecting the CFD flow field predictions.On the other hand, the RTE-FVM model, although much more costly in CPU time, captures the trends of the actual physical behavior for a wide range of circumstances.

Conclusion
This paper has presented a comparison of radiative transfer approaches with an emphasis on the application to circuit breaker modelling arc simulation.The test cases have demonstrated the inaccuracy of the P1 model for small values of the absorption coefficients.The RTE-FVM has consistently provided better results when compared to analytical solutions.Moreover, for realistic circuit breaker chamber dimensions and physical parameters, results show a highly nonphysical behavior for the P1 model together with numerical instabilities inherent to the model itself (as an approximation of the governing equation is actually solved).The RTE-FVM on the other hand performs quite well and with a trustworthy behavior.
The P1 model is traditionally used for this type of application for one decade because of its simplicity and efficiency.However, this method leads to inaccurate and even non-physical solutions in some situations involving long propagation distances in 2D geometries.In addition, a major concern in practical applications of this model is the lack of a physical realistic set of boundary conditions.
Not surprisingly, the number of directions for the RTE-FVM discretization is crucial.The minimum number of directions required to obtain realistic results is approximately 6 × 6 on crude meshes, and around 10 × 10 is better for the kind of meshes used in this study and in circuit breaker modeling.Unit normal vector to the wall s:

Nomenclature
Unit vector in the s-direction (m) r: Radial position (m).

Figure 2 :
Figure 2: Results for test case I.

Figure 5 :
Figure 5: Geometry for test case III.

Figure 6 :
Figure 6: Test case III.Source term for an opaque band.

3 )Figure 7 :
Figure 7: Test case III.Source term for a transparent band.

Figure 8 :
Figure 8: Test case III.Incident radiative flux for a transparent band.

Table 1 :
Effect of the number of directions on the accuracy of the simulation.