Finite-Element Method for Radiation Diffusion in Nonisothermal and Nonhomogeneous Media
Following from: Diffusion approximation in multidimensional radiative transfer problems; Radiation of isothermal volumes of scattering medium: An error of the diffusion model
Leading to: Radiative transfer in multidimensional problems: A combined computational model; Radiation heat transfer in supersonic nozzle of solid-propellant rocket engine
Thermal radiation of a scattering medium is described in a diffusion approximation by the linear boundary-value problem with a simple second-order differential operator:
Therefore, it is not difficult to proceed to the variational formulation of the problem when the solution yields the minimum of the following functional:
The order of the differential operator in (2) is reduced by one compared with (1a). As a result, the transition to the variational formulation is convenient for use of the finite-element method (FEM) with a low local approximation order of the function to be found. With the unified program blocks of FEM, the arrangement of computer codes for calculation of the multidimensional radiation field is considerably simplified.
The extensive possibilities of FEM in the description of complex form regions with nonhomogeneous properties of the medium provide the universality of the algorithm. In accordance with the general FEM scheme, we divide a computational region into subregions, or finite elements, and suppose that the unknown function can be presented as a polynomial in every element. For the sake of simplicity, according to a paper by Dombrovsky and Barkova (1986), only triangular elements are considered with linear approximation of the function Iλ0 in the element. The total functional χλ is presented as a sum of functionals χλe for individual elements. In every element,
where Iλ,i0 are the nodal values of Iλ0 at the element corners and Ni is the linear form function that in the Cartesian coordinates is as follows:
where Δe is the element area. The condition of minimum for every functional χλe gives equations
having the following matrix form:
Matrices Aik, Bik, Cik and vectors Fi, Gi are calculated under the assumption of constant Dλ, αλ, and T in the element and linear γw and Bλ(Tw) at the element boundary. As a result, we have the following relations for the plane geometry:
If node i is not a boundary node, then Cik ≡ 0 and Gi ≡ 0. In the opposite case, for the boundary side (i,l) we obtain
where hil is the length of the corresponding side of the element, and δil is the Kronecker delta.
For problems with axial or spherical symmetry, the expressions (6-7) have a slightly different form. For example, in axisymmetrical cases matrices Aik, Bik and vector Fi should be multiplied by the average radius of the element, whereas matrix Cil and vector Gi, by the average radius of the boundary side.
By transfer to the whole computational region and taking into account the additivity of the functional, we obtain the system of linear algebraic equations for the unknown function at Np nodes of the finite-element mesh:
where matrices and vectors contain the components of the corresponding matrices and vectors of individual elements placed in accordance with the general nodes numbering. The formation of coefficients and subsequent band-matrix solution for Eq. (8) can be accomplished by using the codes given in a monograph by Zienkiewicz et al. (2005).
Some computational results for the emissivity of a nonscattering medium in square channels are given in Tables 1 and 2. The same variant is calculated by use of the analytical solution (Table 1) and numerically by use of the FEM (Table 2). The symmetrical triangulation scheme with equal numbers of the intervals Nint of the rectangular mesh were used for the computational region (see Fig. 1).
Table 1. Analytical solution for emissivity of a nonscattering medium in a square channel
|τ0 = 0.2||0.6||1.0||3.0||5.0|
Table 2. Numerical solution for emissivity of a nonscattering medium in a square channel
|x/x0||ε : (Nint = 10)/ (Nint = 20)|
|τ0 = 0.2||0.6||1.0||3.0||5.0|
Figure 1. The schematic of a finite-element pattern for computational region 0 ≤ x,y ≤ 1 at Nint = 10.
A comparison of the data presented in Tables 1 and 2 confirms the high accuracy of the numerical calculations, even at Nint = 10, at least in the case of lower optical thicknesses. It is interesting that some oscillations of the calculated values of emissivity at high optical thickness are observed. The oscillations are similar to those in the analytical solution with a limited number of series terms Nser. The increase of Nser in the analytical solution influences the solution convergence in the same way as that of Nint by employing the FEM. Therefore Eq. (6), from the article Radiation of isothermal volumes of scattering medium: An error of the diffusion model, can also be used to estimate the required number of mesh intervals in the FEM solution. The conclusion of the abovementioned article on the maximum number of characteristic functions for nondegenerated two-dimensional problems is also applicable to the FEM numerical solution, i.e., one can rely on the sufficiently detailed triangulation with Nint ≤ 50 for the majority of practical problems.
Dombrovsky, L. A. and Barkova, L. G., Solving the two-dimensional problem of thermal-radiation transfer in an anisotropically scattering medium using the finite element method, High Temp., vol. 24, no. 4, pp. 585-592, 1986.
Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z., The Finite Element: Its Basis and Fundamentals, 6th ed., Oxford: Elsevier, 2005 (see also the 1st ed.: Zienkiewicz, O. C., The Finite Element Method in Engineering Science, London: McGraw Hill, 1967).