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:
(1a) |
(1b) |
Therefore, it is not difficult to proceed to the variational formulation of the problem when the solution yields the minimum of the following functional:
(2) |
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,
(3) |
where I_{λ,i}^{0} are the nodal values of I_{λ}^{0} at the element corners and N_{i} is the linear form function that in the Cartesian coordinates is as follows:
(4) |
where Δ_{e} is the element area. The condition of minimum for every functional χ_{λ}^{e} gives equations
(5) |
having the following matrix form:
(6a) |
(6b) |
(6c) |
Matrices A_{ik}, B_{ik}, C_{ik} and vectors F_{i}, G_{i} are calculated under the assumption of constant D_{λ}, α_{λ}, and T in the element and linear γ_{w} and B_{λ}(T_{w}) at the element boundary. As a result, we have the following relations for the plane geometry:
(7a) |
If node i is not a boundary node, then C_{ik} ≡ 0 and G_{i} ≡ 0. In the opposite case, for the boundary side (i,l) we obtain
(7b) |
where h_{il} 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 A_{ik}, B_{ik} and vector F_{i} should be multiplied by the average radius of the element, whereas matrix C_{il} and vector G_{i}, 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 N_{p} nodes of the finite-element mesh:
(8) |
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 N_{int} 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
x/x_{0} | ε | ||||
τ_{0} = 0.2 | 0.6 | 1.0 | 3.0 | 5.0 | |
0 | 0.344 | 0.738 | 0.897 | 0.999 | 1.000 |
0.1 | 0.344 | 0.736 | 0.895 | 0.998 | 1.000 |
0.2 | 0.342 | 0.730 | 0.888 | 0.998 | 1.000 |
0.3 | 0.340 | 0.720 | 0.876 | 0.996 | 1.000 |
0.4 | 0.336 | 0.705 | 0.858 | 0.992 | 0.999 |
0.5 | 0.331 | 0.686 | 0.834 | 0.984 | 0.998 |
0.6 | 0.325 | 0.661 | 0.801 | 0.970 | 0.995 |
0.7 | 0.318 | 0.631 | 0.758 | 0.942 | 0.984 |
0.8 | 0.310 | 0.593 | 0.703 | 0.886 | 0.953 |
0.9 | 0.300 | 0.548 | 0.631 | 0.776 | 0.858 |
1 | 0.289 | 0.495 | 0.538 | 0.550 | 0.552 |
Table 2. Numerical solution for emissivity of a nonscattering medium in a square channel
x/x_{0} | ε : (N_{int} = 10)/ (N_{int} = 20) | ||||
τ_{0} = 0.2 | 0.6 | 1.0 | 3.0 | 5.0 | |
0 | 0.345/0.344 | 0.740/0.738 | 0.900/0.898 | 1.016/1.004 | 1.043/1.013 |
0.1 | 0.344/0.344 | 0.736/0.736 | 0.895/0.896 | 0.994/1.003 | 0.991/1.013 |
0.2 | 0.343/0.342 | 0.732/0.730 | 0.891/0.889 | 1.015/1.002 | 1.043/1.013 |
0.3 | 0.340/0.340 | 0.720/0.720 | 0.876/0.877 | 0.991/1.001 | 0.990/1.012 |
0.4 | 0.336/0.336 | 0.707/0.706 | 0.862/0.859 | 1.010/0.997 | 1.043/1.012 |
0.5 | 0.331/0.331 | 0.686/0.686 | 0.834/0.835 | 0.981/0.990 | 0.989/1.011 |
0.6 | 0.325/0.325 | 0.663/0.662 | 0.805/0.802 | 0.991/0.976 | 1.040/1.008 |
0.7 | 0.318/0.318 | 0.631/0.631 | 0.759/0.759 | 0.941/0.948 | 0.979/0.998 |
0.8 | 0.310/0.310 | 0.595/0.594 | 0.707/0.704 | 0.915/0.894 | 1.012/0.870 |
0.9 | 0.300/0.300 | 0.547/0.548 | 0.629/0.631 | 0.776/0.784 | 0.866/0.879 |
1 | 0.289/0.289 | 0.493/0.494 | 0.534/0.535 | 0.550/0.554 | 0.573/0.547 |
Figure 1. The schematic of a finite-element pattern for computational region 0 ≤ x,y ≤ 1 at N_{int} = 10.
A comparison of the data presented in Tables 1 and 2 confirms the high accuracy of the numerical calculations, even at N_{int} = 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 N_{ser}. The increase of N_{ser} in the analytical solution influences the solution convergence in the same way as that of N_{int} 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 N_{int} ≤ 50 for the majority of practical problems.
REFERENCES
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).