A direct collocation meshless approach with upwind scheme for radiative transfer in strongly inhomogeneous media

Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

A direct collocation meshless approach with upwind scheme for radiative transfer in strongly inhomogeneous media Kang Luo, Zhi-Hong Cao, Hong-Liang Yi n, He-Ping Tan School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, PR China

a r t i c l e in f o

abstract

Article history: Received 31 August 2013 Received in revised form 12 November 2013 Accepted 21 November 2013 Available online 4 December 2013

A direct collocation meshless (DCM) method with upwind scheme is employed for solving the radiative transfer equation (RTE) for strongly inhomogeneous media. The trial function is constructed by a moving least-squares (MLS) approximation. The upwind scheme is implemented by moving the support domain of MLS approximation to the opposite direction of each streamline. To test computational accuracy and efficiency of the upwind direct collocation meshless (named UPCM) method, various problems in 1-D and 2-D geometries are analyzed. For the comparison, we also present cases of both the DCM method for the first-order RTE (employed by Tan et al. [1]) and the DCM for the MSORTE (a new second-order form of radiative transfer equation proposed by Zhao et al. [2]). The results show that the proposed method is more accurate and stable than the DCM method (no upwinding) based on both the RTE and MSORTE. Computationally, it is also faster. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Radiative transfer Upwind collocation meshless method Strongly inhomogeneous media Numerical properties

1. Introduction The radiative transfer equation (RTE) is an integrodifferential equation widely used in different research areas to model wave motions such as light propagation through a turbulent atmosphere, electromagnetic waves propagating in plasmas, light propagation in biological tissue, or radiative heat transfer. Although it is a very well-known equation, analytical solutions are only available for the simplest problems. Therefore, numerical solutions to the RTE are necessary in practical applications. Many numerical methods have been developed to solve the RTE in semitransparent media in recent years. These methods can be mainly classified into two kinds. The first kind of method is based on ray tracing, such as the ray tracing method [3,4], the zonal method [5,6], the DRESOR (Distributions of Ratios of Energy Scattered or Reflected) method [7], the Monte Carlo method [8–10], and the

n

Corresponding author. Tel.: þ86 451 86412674. E-mail addresses: [email protected] (H.-L. Yi), [email protected] (H.-P. Tan). 0022-4073/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jqsrt.2013.11.015

discrete transfer method [11–13]. These kinds of methods do not explicitly rely on the differential form of RTE because the simulation process of these methods is more physically based and the ray path is analytically determined. However, these methods have some obvious drawbacks, such as often being time-consuming even for relatively simple problems, or being difficult to deal with anisotropic scattering, or being difficult to be implemented to complex geometries due to the geometrical shielding. The second kind of method is based on the discretization of partial differential equations, such as spherical harmonics approximations [14,15], discrete ordinate methods (DOM) [16,17], the finite volume method (FVM) [18,19], the finite element method (FEM) [20,21], the collocation spectral method [22–24], and the Meshless method [25–27]. The methods based on the discretization of partial differential equation have the advantages of high efficiency and excellent flexibility to deal with multidimensional complex geometries. However, these methods suffer much from numerical oscillations in some cases caused by the convection-dominated property of the RTE. The RTE is a first-order partial differential equation in each angular direction. Compared to the convection-diffusion

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

Nomenclature a a d I i, j, k G M n Nsol N P pj q r S s t T V w x

coefficients for MLS approximation vector of coefficient a diameter of the support domain, m radiation intensity, W/m2 sr general spatial indices incident radiation, W/m2 number of discrete directions unit normal vector total number of solution nodes nodal basis function vector of legendre polynomial legendre polynomial of jth order heat flux, W/m2 distance between points x and xi, m source term of the radiative transfer equation unit vector in a given direction computation time, s temperature, K solution domain weight function vector of optical location

Greek symbols β

67

γ upwind factor κa absorption coefficient, m  1 κs scattering coefficient, m  1 ε wall emissivity μm ; ηm ; ξm direction cosine in direction, m λ support domain amplifying factor s Stefan-Boltzmann constant, W/m2 K4 τ optical length Φ scattering phase function Ω vector of radiation direction Ω solid angle, sr ω scattering albedo Subscripts b i m j, k w

black body node index direction index legendre polynomial order index wall

Superscript m, m′ T

direction index transposition

extinction coefficient, m  1

equation, the RTE can be considered as a special kind of convection-dominated equation without the diffusion term. As mentioned before, the methods based on the discretization of partial differential equations may suffer from nonphysical oscillations in the numerical results, especially for the models with inhomogeneous media where some regions have very small/zero extinction coefficient or discontinuous extinction coefficient. To make the numerical discretization schemes correctly model the transfer process, two kinds of special stabilization techniques are often used: (1) transforming the RTE into a numerical stable equation, for example, the secondorder partial differential equation, and (2) taking various numerical stabilization schemes, such as upwind scheme or artificial diffusion which are often used in FDM (finite difference method), FVM, FEM and the meshless method. The second-order form of radiative transfer equation (SORTE) contains second-order derivative terms which are known to have the characteristic of diffusion and good numerical properties. Besides, the diffusion term introduced in the analytical transformation process is natural and consistent with the original first-order equation. In recent years, three different transformed equations have been proposed [2,28–32]. One is the even-parity (EPRTE) formulation of the RTE mentioned in Refs. [28–30] which is a second-order partial differential equation of the even parity of radiative intensity. Another is the second-order radiative transfer equation (SORTE) proposed recently [31,32] which uses radiative intensity as solution variable and is more convenient to be applied to the complex radiative transfer equation, but has some problems in

dealing with inhomogeneous media due to the existence of the reciprocal of extinction coefficient in the equation. The third is a new second-order form of radiative transfer equation (named MSORTE) proposed to overcome the singularity problem of the SORTE. Zhao et al. [2] applied a meshless method based on the moving least square approximation to verify the versatility and performance of the MSORTE in solving radiative transfer in very strongly inhomogeneous media, and stable and accurate results were obtained. Although the three types of SORTE are more stable than RTE, the iterative solutions of the SORTE are more time consuming. In addition, compared to the first-order RTE, SORTE is of higher differentiability of radiative intensity with respect to the ray trajectory coordinates, which will limit the application of SORTE. An alternative approach to reduce the nonphysical oscillations in the results is the introduction of RTE with the upwind scheme which can obtain the numerical stable solutions under less computational cost compared to SORTE. A number of upwinding schemes have been developed for FDM, FEM and FVM. It is apparent that the upwind effect, achieved by whatever means, is needed only in the direction of flow. However, it is not easy to design such methods for multidimensional cases. Hughes and Brooks [33] introduced the ‘Streamline Upwind (SU) method’ where the artificial diffusion operator was constructed to act only in the flow direction. However, the SU method cannot meet the consistency of the RTE, resulting in excessively diffuse solutions, and at the cost of accuracy. Hughes and Brooks [34] proposed an improved upwind

68

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

scheme ‘Streamline Upwind Petrov-Galerkin (SUPG) method’ which consistently introduced an additional stability term in the upwind direction, but the choice of the related parameters in the method is still not so straightforward. Both the SU and SUPG methods were applied to solve the first-order RTE using the Galerkin finite element method (GFEM) by Liu and Zhao [35,36], and stable results were obtained in some cases. It has been seen from extensive literature survey that no work, to the authors' knowledge, is available that apply the upwind scheme to solve RTE using the meshless method. In this article, we employ the direct collocation meshless method with upwind scheme (UPCM) to solve the first-order RTE. The trial functions are constructed by a moving least-squares (MLS) approximation. The upwind scheme is constructed by moving the support domain of MLS to the streamline direction, which is natural and convenient to conduct. Some cases are examined to test the stability, accuracy and efficiency of the UPCM, and the results are compared with those of the DCM method based on the first-order RTE and MSORTE. Besides, we analyze the relationship between the stability of numerical results and the upwind factors and examine how to obtain the optimal upwind factors for different diameters of support domain.

2. Mathematical formulation 2.1. Discrete-ordinates equation of radiative transfer

ð1Þ m

m

where Ω is the transport direction, and S term, given by Ωm ¼ iμm þ jηm þ kξm

Sm ¼ κ a I b þ

κs M m′ m′m m′ ∑ I Φ w 4π m′ ¼ 1

is the source

ð2aÞ

1  εw ∑ I m′ jnwdsm′ jwm′ π jnwdsm′ j o 0 w

ðΩmd∇Þ2 I m þ Ωmd∇ðβI m Þ ¼ Ωmd∇S

ð4Þ

With boundary conditions Im w ¼ εw I bw þ



1 εw m′
m′
∑jnwdsm′ j o 0 I m′ w ; w nwds π

m m Ωmd∇I m w þ βI w ¼ S ;

jnwdsm j o 0

jnwdsm j 40 ð5aÞ ð5bÞ

The boundary conditions of the MSORTE contain two parts: the inflow boundary and the outflow boundary. The inflow boundary condition Eq. (5a) is conventionally called the Dirichlet boundary or the essential boundary condition, which is the same as the boundary condition of the first-order RTE. The outflow boundary of MOSRTE (Eq. (5b)) is the Neumann boundary or natural boundary, and is constructed by applying the RTE to the outflow boundary points.

To calculate the divergence of radiative heat flux in the energy equation, the incident radiation G distribution needs to be known. It can be calculated as Z Z 2π Z π GðrÞ ¼ Iðr;sÞdΩ ¼ I sin θ dθ dφ ð6Þ 4π

φ¼0

θ¼0

The radiative heat flux in the x and y directions can be calculated as Z qR;x ðrÞ ¼

Z Iðr;sÞðsdex ÞdΩ ¼ 4π

2π

φ¼0

Z

π

θ¼0

I sin θ cos φ sin θ dθ dφ

ð7aÞ ð2bÞ

Z qR;y ðrÞ ¼

With the boundary condition for the opaque surface given by Im w ¼ εw I bw þ

As mentioned in Section 1, three types of second-order radiative transfer equation – EPRTE, SORTE and MSORTE – have been proposed in recent years. Owing to the best stability of the MSORTE in dealing with the inhomogeneous media we follow the argument of Zhao [2] to give the formulation of MSORTE. Directly applying the streaming operator dI=ds ¼ ðΩd∇ÞI to Eq. (1) yields

2.3. Incident radiation and radiative heat flux

Consider the radiative transfer in the enclosure filled with semitransparent media having opaque surfaces. The discrete-ordinate equation of radiative transfer can be written as Ωmd∇I m þ βI m ¼ Sm

2.2. Second-order discrete-ordinates equation of radiative transfer

Z Iðr;sÞðsdey ÞdΩ ¼ 4π

2π φ¼0

Z

π

θ¼0

I sin θ sin φ sin θ dθ dφ

ð7bÞ

ð3Þ 3. Direct collocation meshless discretization

where I m w is the radiation intensity at the boundary at the direction m; β ¼ κa þκ s is the extinction coefficient; κa and κ s are the absorption and scattering coefficients, respectively; μm , ηm , and ξm are the direction cosine; Φ is the scattering phase function; εw is the wall emissivity; nw is the unit normal vector of the boundary surface; sm′ is the unit vector in the direction m′; and wm′ is the weight corresponding to the direction m′.

For the application of DCM, the moving least-squares (MLS) approximation is employed to construct the trial functions using the collocation points, and the DCM is employed to construct the linear equations. Tan et al. [1] applied the DCM method to solve the first-order RTE, and we follow the argument of Lin [37] and Zhao [2] to give the formulation of MLS and DCM.

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

69

And the local polynomial approximation [Eq. (8)] can be formulated into the equivalent nodal basis expansion as

3.1. Moving least-squares approximation Consider a spatial sub-domain V x in the neighborhood of a point x denoted as the domain of the definition of the MLS approximation for the trial function at x. To approximate the distribution of a generic trial function I m ðxÞ in V x , over a number of local nodes {xi}, i ¼ 1; 2; …; n, the m MLS approximation I~ ðxÞ of I m ðxÞ, 8 x A V x , is given by p

m I~ ðxÞ ¼ ∑ pj ðxÞaj ðxÞ ¼ pT ðxÞaðxÞ

ð8Þ

j¼0

N

m I~ ðxÞ ¼ ∑ Ni ðxÞI m ðxi Þ

ð17Þ

i¼1

The first- and second-order partial derivatives of the radiative intensity are required when solving the RTE and MSORTE. Here, the evaluation of the partial derivatives is by the ‘diffuse derivative’ approach [38], and yields N

N

i¼1

i¼1

∇I~ ðxÞ ¼ ∑ ∇Ni ðxÞI m ðxi Þ ¼ ∑ ∇pT ðxÞΦ  1 ðxÞpðxi Þwðjxi  xjÞI m ðxi Þ m

where pj is the jth polynomial basis function and aj is the corresponding expansion coefficient, and pT can be chosen as linear

ð18Þ

ð9Þ

N

N

i¼1

i¼1

∇∇′I~ ðxÞ ¼ ∑ ∇∇′N i ðxÞI m ðxi Þ ¼ ∑ ∇∇′pT ðxÞΦ  1 ðxÞpðxi Þwðjxi  xjÞIm ðxi Þ m

pT ðxÞ ¼ ½1; x; y

ð19Þ

or quadratic: pT ðxÞ ¼ ½1; x; y; x2 ; xy; y2 

ð10Þ

In this article, the complete second-order polynomial basis in Eq. (10) is used for MLS approximation. The least square approach can then be formulated with the help of the local inner product as * p m minJ ðaÞ ¼ min J I~ ðxÞ  I m ðxÞ J 2 ¼ min ∑ p ðxÞaj  I m ðxÞ; aj

x

p

aj

aj

j¼0

j

ð20Þ

+ m

∑ pk ðxÞak  I ðxÞ

k¼0

where ∇ and ∇′ are considered as gradient operator or any first-order partial differential operator. To define the weight function wðrÞ, the fourth-order Wendland's function [39], given by Eq. (20), is chosen that is compactly supported and continuous up to the fourth-order derivatives. ( ð1  r=dÞ6 ð35ðr=dÞ2 þ 18ðr=dÞ þ 3Þ; ðr=dÞ r 1 wðr=dÞ ¼ 0 ðr=dÞ 4 1

ð11Þ x

A solution to the minimization problem (using the extreme value condition ∂J x ðaÞ=∂aj ¼ 0) yields p

∑ 〈pk ðxÞ; pj ðxÞ〉x aj ¼ 〈pk ðxÞ; I m ðxÞ〉x

where r is the distance between xi and x, and d is the diameter of the support domain given by d ¼ λðp þ 1ÞΔh

ð21Þ

ð12Þ

where p is the order of polynomial basis, Δh is the local discrete length, and λ is the support domain amplifying factor.

ð13Þ

3.2. Collocation-based meshless approach

j¼0

which can be written in the matrix form as Φa ¼ f where   Φ ¼ Φkj k ¼ 1…p;j ¼ 1…p ¼ ½〈pk ðxÞ; pj ðxÞ〉x j ¼ 1…p;

a ¼ ½aj j ¼ 1…p ¼ Φ  1 f

ð14bÞ

The basic idea of the DCM approach for solving the RTE is to force the residual of the equations to be zeros on each collocation point. Using the MLS approximation mentioned above and substituting the shape function [Eq. (17)] and its first-order partial derivatives [Eq. (18)] into Eq. (1), we get

f ¼ ½f k k ¼ 1…p ¼ ½〈pk ðxÞ; I m ðxÞ〉x k ¼ 1…p

ð14cÞ

βðxÞI m ðxÞ  ∑ Ni ðxÞβi I m i ;

k ¼ 1…p

ð14aÞ

Substituting Eq. (14) into Eq. (8), and the locally approximated intensity field can be expressed as m I~ ðxÞ ¼ pT ðxÞaðxÞ ¼ pT ðxÞΦ  1 ðxÞfðxÞ

N

i¼1

N

SðxÞ  ∑ Ni ðxÞSi ðxÞ

and the discretization form of the DCM approach can be obtained based on the RTE and its boundary condition as N

N

i¼1

i¼1

∑ fΩmd∇Ni ðxÞ þ βi N i ðxÞgI m i ¼ ∑ N i ðxÞSi ðxÞ

N

¼ ∑ pT ðxÞΦ  1 ðxÞpðxi Þwðjxi xjÞI m ðxi Þ

ð22Þ

i¼1

ð23aÞ

i¼1

ð15Þ A nodal basis function (at node i) from the MLS approximation can then be written as Ni ðxÞ ¼ pT ðxÞΦ  1 ðxÞpðxi Þwðjxi  xjÞ

ð16Þ

N

∑ Ni ðxw ÞI m i ¼ εw I bw þ

i¼1

1  εw m′ m′ ∑jnwdsm′ j o 0 I m′ w jnwds jw π ð23bÞ

We also obtained the discretization form of the DCM approach based on the MSORTE and its inflow and outflow

70

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

Trial functions

I (Ω )

Shifted local sub-domain Test functions

Test functions

Fig. 2. DCM upwinding scheme.

Fig. 1. DCM method without upwinding.

boundary conditions as N

N

i¼1

i¼1

m ∑ fðΩmd∇Þ2 Ni ðxÞ þ βi Ωmd∇N i ðxÞgI m i ¼ ∑ Ω d∇N i ðxÞSi ðxÞ

ð24aÞ N

∑ N i ðxw ÞI m i ¼ εw I bw þ

i¼1

1  εw m′ m′ ∑jnwdsm′ j o 0 I m′ w jnwds jw ; π

jnwdsm j 40

Fig. 3. DCM upwinding scheme: specification.

ð24bÞ

N

N

i¼1

i¼1

∑ fΩmd∇Ni ðxw Þ þ βi N i ðxw ÞgI m i ¼ ∑ N i ðxw ÞSi ðxw Þ;

jnwdsm j o 0

ð24cÞ By applying Eq. (23) or Eq. (24) to each collocation point, we can obtain an equation set in each direction, and by iterative computations, we can obtain the solutions of the DCM approach based on both the first-order RTE and MSORTE. The solution procedure of the UPCM is the same as that of the DCM approach based on the RTE, except for the support domain of the trial function. 4. DCM upwinding scheme The DCM method for radiative transfer is based on the discrete ordinate equation. The residual of the equation is forced to be zero on each collocation point, and the MLS approximation is employed for constructing the trial functions. Different spaces for the test and trial functions can be used, as shown in Fig. 1. Therefore one of the very natural ways to construct the upwinding scheme is to choose different trial and test functions. This can be done by a number of ways. For example, the local sub-domain can be shifted opposite to each ray direction, as shown in Fig. 2. Here, we use the same spaces for the trial and test functions, that is, we employ the same support and the same interpolation scheme (MLS) for the trial function. In this case, the local sub-domain at xi is no longer coincident with the support for the test functions at xi, but the size is the same.

In particular, the local sub-domain can be shifted from xi to xi  γr i Ω, as shown in Fig. 3, where Ω is the ray direction, r i is the size of the support domain for the test functions, which is equal to the size of the local domain at xi, and γ is the upwind factor which is related to the support domain amplifying factor λ [in Eq. (21)]. The optimal upwind factor is given by ( γ¼

1

1 rλ o 1:5

1:45

1:5 rλ

ð25Þ

Unlike the upwind scheme of the FDM, FVM and FEM, the upwind factor of the meshless method for RTE can be greater than one.

5. Numerical examples To investigate the accuracy and the computational cost of the UPCM method in solving the first-order RTE, extensive numerical calculations are performed and the results are compared with those of the DCM method based on the first-order RTE and MSORTE. Examples for 1D and 2D problems are given as follows.

5.1. One-dimensional test cases

Case 1. Gaussian-shaped radiative source term inside a plane-parallel slab with black surfaces

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

To show that the DCM based on RTE suffers from numerical oscillations, and the UPCM based on RTE can eliminate oscillations, we consider the radiative transfer problems in a non-scattering uniform refractive index plane-parallel medium having black surfaces. The radiative source term of the medium has a profile similar to the Gaussian function. A similar problem was also investigated in Refs. [2,40,41], and the governing equation of this problem is modeled as follows: μ

2 dI 2 þκ a I ¼ e  ðz  cÞ =α ; dx

z; c A ½0; 1

ð26Þ

with the boundary condition given by Ið0; μÞ ¼ κ a 1 e  c

2

=α2

;

μ40

ð27aÞ

2

;

ð27bÞ

Ið1; μÞ ¼ κa 1 e  ð1  cÞ

=α2

μo0

where c is a parameter denoting the location of the peak of the Gaussian profile and α is a parameter which can change the shape of the Gaussian profile. The exact solution to this problem in the case of μ 40 can be written as

  2     pffiffiffi κa z α π κa α κa z þc Iðz; μÞ ¼ Ið0; μÞexp   exp  μ 4μ μ 2μ      ακ a c z ακ a c  erf  erf þ þ ð28Þ α α 2μ 2μ This problem features large gradient of intensity distribution in the medium which was investigated by many researchers [32,40] to test the performance of the numerical method for RTE, such as the FEM and meshless method. The UPCM method based on RTE is applied to this case for c¼0.5, α ¼0.02, and the three different optical thicknesses, namely 0.1, 1 and 10. The upwind factor and support domain amplifying factor are set to 1, and about 80 nodes and 20 angular directions are considered for this case. The distribution of the radiative intensity in the direction μ ¼ 0:51 obtained by the UPCM method is shown in Fig. 4(a) and compared with the exact analytical solutions. The results of DCM based on the first-order RTE under the same condition are also plotted in Fig. 4(b). By comparison, we find that the results of the UPCM method are more stable and more accurate, and have a

good agreement with the analytical solutions, but the results obtained by the DCM method show enormous spurious numerical oscillations. This case demonstrates the good numerical properties of the UPCM scheme. Case 2. Infinite slab with inclusion layers In the above case, we examined the numerical properties of the UPCM method for radiative transfer in a homogeneous slab. In this case, for further investigation of the accuracy, stability and efficiency of the UPCM method, we explore radiative transfer in an inhomogeneous slab with two absorbing-scattering inclusion layers, and compare the results with those of DCM based on the first order RTE and MSORTE. The configuration of the slab is shown in Fig. 5. The thickness of the slab is L and the thickness of the layers is d¼0.2L. The slab is uniformly divided into 200 subdivisions and the angular discretization with 20 discrete directions is used for the Gauss integration. The extinction coefficient of the inclusion layers is βI ð ¼ τ=LÞ and the other portion of the slab is a void. Both the walls of the slab are black. The temperatures of the left wall and inclusion layers are T L ð ¼ 1000 KÞ and T g ð500 KÞ, respectively, and the right wall is kept cold (T R ¼ 0K). Fig. 6 presents the heat flux and incident radiation distributions for different values of extinction coefficient and scattering albedo. The numerical results of the UPCM method (λ ¼ 1; γ ¼ 1) based on the first-order RTE are compared with those of the DCM method based on MSORTE, and the results obtained by the discontinuous spectral element method [41] are taken as the reference (exact) results. As shown in Fig. 6(a), the dimensionless heat flux has a ladder-like distribution at low extinction coefficient, resulting from attenuation of radiant energy only in the inclusion layers. For the large extinction coefficient, the heat flux sharply declines in the first inclusion layer and has an arc distribution in the second inclusion layer. Similar effects can be found in the distribution of the incident radiation as shown in Fig. 7(a). In addition, the heat fluxes at the left wall decrease with the increase of scattering albedos as shown in Fig. 6(b). However, the opposite effect can be found in the distribution of the incident radiation as shown in Fig. 7(b), which

0.08

0.08

κα = 0.1

Exact UPCM

κα = 0.1

Exact RTE 0.06 I,W/m2sr

I,W/m2sr

0.06

71

κα = 1 0.04

κα = 1

0.04

0.02 0.02

κα = 10

κα = 10 0.00 0.0

0.2

0.4

0.6 x

0.8

0.00 1.0

0.0

0.2

0.4

0.6

0.8

1.0

x

Fig. 4. Radiative intensity distribution in the infinite slab by the meshless method with and without the upwind scheme: (a) with upwind scheme and (b) without upwind scheme.

72

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

gradient in the frequency domain will have high-frequency harmonic component, which will induce large numerical errors for the central difference discretization of the firstorder RTE as mentioned in ref. [2]. Furthermore, the spurious oscillations are more intense in the location with small/zero extinction coefficient. Fig. 8 presents the relative error of the UPCM and DCM (based on the first-order RTE and MSORTE) for β¼10 and 1. The relative error is defined as J Gsolved  Gexact J =J Gexact J . It is seen that larger error appears for a bigger value of the extinction coefficient. Besides, the UPCM based on RTE shows the lowest level of relative error and has the fastest convergence speed. To investigate the computational cost of the above numerical methods, three computer codes for these three numerical methods have been written and run on a personal computer with Intel Core 2 Quad processor. The computation times are shown in Table 1. From Table 1, we find that the computation time increases with the increase in scattering coefficient (κ s ). Besides, the UPCM based on RTE requires the least computation time. This is because the UPCM method based on the first-order RTE needs less iterative steps than the DCM method based on the firstorder RTE and needs less time in each iterative step than the DCM method based on MSORTE. Furthermore, the grid number required in the UPCM is lower than that in the DCM based on the first-order RTE and MOSRTE to obtain the same target accuracy, as seen in Fig. 8. The influence of upwind factor on the relative error is investigated for β ¼ 10 and ω ¼ 0. Two different values of support domain

means that the high scattering albedos can make the radiative energy spread to all angular directions. It is seen that the heat flux and incident radiation obtained by the upwind meshless method based on the first-order RTE agree very well with the reference results, and are more accurate than those of DCM based on MSORTE, especially at the place with discontinuous extinction coefficients as can be seen from Fig. 6. However, the results obtained by the meshless method based on first-order RTE are totally spoiled by the spurious oscillations as shown in Fig. 6(c)–(d) and Fig. 7(c)–(d), resulting from the steep gradient of intensity distribution. A function with steep

L

TL

β =0

β =β I

L /5

L /5

β =0

β =βI

β=

L /5

0

TR

L /5

z Fig. 5. Configuration of the inhomogeneous slab with two inclusion layers.

a

b

1.0

1.0

Exact UPCM MSORTE

0.8

q(z)/σ T4

q(z)/ σ T 4

0.8 0.6 Exact UPCM MSORTE

0.4

β = 0.1, 1, 10 ω=0

0.2

ω = 0.1, 0.5, 0.9

0.6

β = 10

0.4

0.2

0.0

0.0 0.0

0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

z/L

c

0.8

1.0

d

2.00

β = 10, 1, 0.1 ω=0

1.50

1.75

Exact RTE

1.75

Exact RTE

1.50 1.25 q(z)/σ T4

1.25 q(z)/ σ T4

0.6 z/L

1.00 0.75

1.00

0.50

0.25

0.25

0.00

ω = 0.1, 0.5, 0.9 β = 10

0.75

0.50

0.00 0.0

0.2

0.4

0.6 z/L

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

z/L

Fig. 6. Radiative heat flux distribution in the inhomogeneous slab by the UPCM method based on RTE and the DCM method based on RTE and MSORTE for different extinction coefficients and scattering albedos.

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

2.1

3.0

Exact UPCM MSORTE

2.7

1.8

2.4 Exact UPCM MSORTE

1.2 0.9 β = 0.1, 1, 10

0.6

ω = 0.1, 0.5, 0.9

2.1 G(z)/ σ T4

G(z)/ σ T

4

1.5

73

β = 10

1.8 1.5 1.2 0.9

ω=0

0.6

0.3

0.3 0.0 0.0

0.0 0.0

0.2

0.4

0.6

0.8

1.0

0.2

0.4

0.6

0.8

1.0

z/L

z/L

6 Exact RTE

5 β = 10, 1, 0.1

4

ω=0

G(z)/ σ T4

G(z)/ σT4

4

Exact RTE

5

3

3

2

2

1

1

0

ω = 0.1, 0.5, 0.9 β = 10

0 0.0

0.2

0.4

0.6

0.8

0.0

1.0

0.2

0.4

z/L

0.6

0.8

1.0

z/L

Fig. 7. Incident radiation distribution in the inhomogeneous slab by the UPCM method based on RTE and the DCM method based on the RTE and the MSORTE for different extinction coefficients and scattering albedos.

102

β = 10 β=1

101 RTE 100 Relative error

amplifying factor (λ ¼ 1; 1:5) are considered and the results are presented in Fig. 9. It is seen that the relative error reduces to a minimum as γ increases to its optimal value, and the numerical results will become instable with the further increase of γ. Besides, the optimal upwind factor γ depends on the support domain amplifying factor λ, which is given in Eq. (25). This case demonstrates good numerical stability, high accuracy and efficiency of the upwind collocation meshless method for solving radiative transfer in strongly inhomogeneous media with a discontinuous distribution of the extinction coefficient.

10-1 10-2

MSORTE

-3

10

10-4 UPCM 10-5 10-6

5.2. Two-dimensional test cases

100

200

300

400

500

600

700

800

900

1000

Node number

Case 1. Square enclosure with a Gaussian emissive field In this case, we consider radiative transfer in a square enclosure with a non-scattering medium. The emissive filed in the medium has a Gaussian profile. This case can be considered as an extension of the one-dimensional case studied in Section 5.1 to two dimensions, and is taken here to especially demonstrate the stability and accuracy of the UPCM methods in multi-dimensions. The governing equation of this problem is given as μ

dI dI þη þκ a I ¼ S; dx dy

x; y A ½0; 1

ð29aÞ

Fig. 8. Spatial convergence test of the UPCM method based on RTE and the DCM method based on the RTE and the MSORTE for solving the incident radiation for different optical thicknesses of inclusions.

S ¼ e  ½s~ ðx;yÞ  c

2

=α2

ð29bÞ

in which the incident direction is selected as μ ¼ 1 and pffiffiffi η ¼ 0; s~ ðx; yÞ ¼ 2x is a distance parameter in the incident direction; the emission profile parameters are selected as pffiffiffi c ¼ 2=2 and α ¼ 0:02. The unit square enclosure is depicted in Fig. 10(a) and the source term S is presented in Fig. 10(b). The boundary condition is prescribed on the

74

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

inflow boundary with null emission. Here, the UPCM method based on the first-order RTE is applied to this case for λ ¼ 1 and γ ¼ 1 and the four different extinction coefficients, namely 1, 5, 10 and 100 m  1. The square enclosure is divided uniformly into 1681 nodes. The meshless method based on the first-order RTE and MSORTE are also used to solve the intensity distribution in the medium, and the results are presented in Fig. 11(a) and (b), respectively. The UPCM results are shown in Fig. 11(c). As shown in Fig. 11(a), strong spurious wiggles appear in the solution obtained by the meshless method based on the RTE, and the spurious wiggles have a short wavelength about two times that of grid spacing. The spurious wiggles are caused by the steep gradient of intensity distribution. A function with steep gradient in frequency domain will have a high-frequency harmonic component, which will induce numerical errors for the center difference discretization of the first-order RTE, as indicated by the theoretical analysis presented by Zhao et al. [2]. Furthermore, the magnitude of the spurious wiggles tends to decrease with the increase in the extinction coefficient. It is seen from Fig. 11(b) that the meshless method based on MSORTE can obtain stable results, which is attributed to the second-order differential (diffusion) term. However, those results suffer from spurious diffusion due to the first-order continuity requirements of intensity for the MSORTE, which leads to the smooth intensity distribution in some place where there should have sharp corner as shown in Fig. 11(b). Fig. 11(c) presents the numerical solution of the UPCM method based on RTE, and it is seen that the major characteristics agree qualitatively well with the theoretically predicted trends given in Ref. [2], which demonstrate good numerical stability and high accuracy of the UPCM method in solving multidimensional problems in strongly inhomogeneous media. In this case, the numerical properties of the UPCM method based on the first-order RTE are shown to be better than those of the DCM method based on the MSORTE.

Table 1 Comparison of the computation times for DCM and UPCM. Equation Numerical method

ks

RTE

0 2 0.5 4 1 4 3 6 5 8 9 22

DCM

Iteration number

Computation time (s) 1.246 2.258 2.477 2.693 5.065 6.089

MSORTE DCM

0 0.5 1 3 5 9

2 4 4 5 6 7

2.033 3.528 3.628 4.873 8.312 10.42

RTE

0 2 0.5 4 1 4 3 5 5 6 9 12

1.248 2.262 2.482 2.545 3.684 4.672

UPCM

10-1 λ =1 λ = 1.5

Relative error

10-2

10-3

10-4

0.5

0.6

0.7

0.8

0.9

1.0

1.1

1.2

1.3

1.4

1.5

Case 2. Radiative equilibrium in a square enclosure

Upwind factor γ

In this case, the radiative equilibrium (∇dq ¼ 0) in an isotropically scattering media is considered. The temperature

Fig. 9. Relative error of the incident radiation with the change of upwind factor for different support domain amplifying factors.

1

1 Non-participating 0.8

0.8

0.6

0.6 y

y

β = βI

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6 x

0.8

1

0

0

0.2

0.4

x

0.6

0.8

Fig. 10. Configuration of (a) square enclosure with a Gaussian emissive field and (b) source term S.

1

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

75

Fig. 11. Radiative intensity distribution by (a) the DCM method based on the RTE, (b) the DCM method based on the MSORTE, and (c) the UPCM method based on the RTE for different extinction coefficients and scattering albedos.

of the bottom wall is kept T w1 ¼ 1000 K, and the other walls are kept cold (0 K). As depicted in Fig. 12, the ‘boundary induced ray effect’ exists because of the non-uniform and

discontinuity loading from the boundaries. As a result, fluctuations are visible in the heat flux distribution at the top wall while using the numerical methods based on the

76

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

Fig. 11. (continued)

Top wall (T = 0K) ∇

Ω

0K

0K

Bottom wall (T =1000K) Fig. 12. Schematic diagram to illustrate the ‘boundary induced ray effect’.

discretization form of RTE. This case was studied by several researchers [41–43] and serves as a good test case for verifying the performance of the numerical methods developed. Therefore, the stability and accuracy of the UPCM method for the radiative equilibrium problems need to be tested by this case, and the UPCM results are compared with the quasi-exact solution of Crosbie and Schrenker [42]. The results by other numerical methods were reported in Refs. [41,43]. The UPCM method based on the first-order RTE is applied to this case for λ ¼ 1 and γ ¼ 1. The spatial grid system used here is 15  15 and the radiative parameters

are set to be β¼1 m  1, ω ¼ 0 and ε ¼ 1. First, three different angular discretization schemes, namely, N θ  Nφ ¼ 8  16, Nθ  Nφ ¼ 12  24 and Nθ  Nφ ¼ 16  32 are adopted. The dimensionless net radiative heat flux distributions on the top wall are shown in Fig. 13(a) and compared with the quasi-exact solution by Crosbie and Schrenker [42]. It can be seen that obvious ‘wiggles’ exist in the result of the UPCM when N θ  Nφ ¼ 8  16, which results from the lesser accuracy of the angular discretization, while with the refining of the angular discretization, the ‘wiggles’ are effectively mitigated and smooth results are obtained. Fig. 13(b) plots the radiative heat flux distribution at the top wall of the square enclosure obtained by the UPCM method, the standard Galerkin spectral element method (GSEM) [43] and the discontinuous spectral element method (DSEM) [41] under the same spatial and angular discretization, respectively, and they are all compared to the quasi-exact solution of Crosbie and Schrenker [42]. It is seen that the UPCM results are free of ‘wiggles’ while the GSEM and DSEM still suffer from it at the same angular discretization of N θ  Nφ ¼ 16  32. It is well known that two types of errors exist in the numerical solution of the RTE, namely the ‘ray effect’ and ‘false scattering’. Similar to other numerical methods, such as the GSEM and the DSEM, simultaneous refinements of both the spatial and the angular discretizations are needed for the UPCM to obtain accurate results. Owing to its upwind scheme, the UPCM method suffers little from false scattering. Therefore, the ‘wiggles’ in the results obtained by the UPCM method can be effectively mitigated and

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

77

Fig. 13. Dimensionless net radiative heat flux distributions on the top wall of the square enclosure: (a) by the UPCM method under different angular discretizations, and (b) by the GSEM, DSEM and UPCM methods under the same angular discretization.

0.45

0.24

0.40 0.35

0.22

4

qw/ Tw

4

qw/ Tw

0.30 0.25 0.20

= 0.1 =1 = 10

0.15 0.10

0.18 0.16 0.14

0.05 0.00 0.0

=0 = 0.5 =1

0.20

0.2

0.4

0.6

0.8

0.0

1.0

0.2

0.4

0.6

0.8

1.0

x/L

x/L

0.25

4

qw/ Tw

0.20

= 0.1 = 0.5 =1

0.15

0.10

0.05

0.00 0.0

0.2

0.4

0.6

0.8

1.0

x/L Fig. 14. Dimensionless net radiative heat flux distributions on the top wall of the square enclosure for various radiative parameters: (a) effect of β, (b) effect of ω, and (c) effect of ε.

fewer angular directions are required to achieve the same numerical precision. Besides, the ray effect can also be bated effectively through angular refinement. Generally, the UPCM method is convergent, accurate and effective to solve the radiative transfer and can effectively bate the ‘boundary-induced ray effect.’ With the validation of the UPCM method for the radiative equilibrium problem, we also investigate the effect of various radiative parameters, namely, the extinction

coefficient β, scattering albedo ω and wall emissivity ε. The spatial grid system used here is 15  15 and the angular space is discretized into Nθ  N φ ¼ 16  32. Fig. 14 (a) illustrates the dimensionless heat flux distribution at the top wall of the square enclosure filled with non-scattering media (ω¼0) enclosed in black walls under different values of β, namely, β ¼ 0:1, 1 and 10 m  1. It is seen that the results increase with the decrease of extinction coefficients. The effect of scattering albedo ω is shown in Fig. 14(b) for

78

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

β¼1 m  1 and black walls. The dimensionless heat flux distribution at the top wall is demonstrated to be little affected by the scattering albedo. Fig. 14(c) shows the results for different wall emissivities ε, namely, ε ¼ 0:1, 0.5 and 1, and it can be seen that the dimensionless heat flux increases with the increase in wall emissivity. Case 3. Anisotropically scattering in semicircular enclosure with internal circle In this case, the UPCM method is applied to simulate the radiative transfer in an irregular geometry with an anisotropically scattering media. A semicircular enclosure with an inner circle is considered as shown in Fig. 15. Radii of the semicircle and inner circle are 1 m and 0.2 m, respectively. The medium enclosed by the semicircle and the inner circle is maintained at a constant temperature of 1000 K. The enclosure walls are assumed cold and black. The spatial grid system used here is 2895 and the angular space is discretized into Nθ  N φ ¼ 16  32. Many researchers have solved this problem [44–46]. The scattering phase function of the media is the F2 phase function given by Eq. (28) with the asymmetry factor of 0.66972: 8

Φ ¼ ∑ C j P j ðμÞ

ð30Þ

j¼0

where the P j are the Legendre polynomials; the C j are the expansion coefficients: C 0 ¼ 1:0, C 1 ¼ 2:00917, C 2 ¼ 1:56339, C 3 ¼ 0:67407, C 4 ¼ 0:22215, C 5 ¼ 0:04725, C 6 ¼ 0:00671, C 7 ¼ 0:00068, and C 8 ¼ 0:00005.

Tw= 0K ew = 1 Rc = 0.2

0.4

Tg = 1000 K

6. Conclusions In the present work, the direct collocation meshless method with the upwind scheme (UPCM) is applied to solve the radiative transfer equation for the first time, and the MLS approximation is employed for constructing the trial functions. The upwind scheme is constructed by shifting the local sub-domain of MLS approximation to the opposite direction of each streamline for the firstorder RTE. The relationship between the upwind factors and support domain amplifying factors is also provided. Several cases for 1D and 2D radiative transfer problems have been considered to investigate the accuracy, stability and computational cost of the UPCM method. For the sake of comparison, these cases have also been solved by the

y

Rs = 1

The upwind scheme is constructed by shifting the local sub-domain from xi to xi  γr i Ω in each streamline direction. Although this step will increase the computational cost, it is outside the loop and is conducted only once for a case. However, the DCM method based on MSORTE increases the computational time in each iterative step. Besides, to obtain the same target accuracy, the meshless method with the upwind scheme needs a smaller support domain and fewer iterative steps than the direct meshless method without the upwind scheme. Therefore, the UPCM method is more efficient than the DCM method based on MSORTE. A comparison of the non-dimensional radiative heat flux on the bottom wall is presented in Fig. 16(a) and (b). Two different absorption coefficients of 1 and 0.1 m  1 are considered using the UPCM and DCM methods based on the first-order RTE. The UPCM solutions are in good agreement with the exact solutions presented by Byun et al. [44], and the maximum deviation is 1.8%. However, the results obtained by the DCM method based on RTE are unstable, especially in the left and right corners of the semicircular enclosure as shown in Fig. 16. In addition, the computational time of the UPCM is almost the same as that of DCM based on the RTE. This case demonstrates that the UPCM based on the RTE can stably and accurately solve radiative transfer in irregular enclosure.

x

Fig. 15. Schematic of semicircular enclosure with the inner circle.

0.6

0.10 0.09

0.5

qw / σ T

4 g

qw / σ T

4 g

0.08

0.4

0.3

0.06

DCM based on RTE UPCM based on RTE Ref. [44]

0.2 -1.0 -0.8 -0.6 -0.4 -0.2

0.0 x

0.2

0.4

0.07

0.05

0.6

0.8

1.0

DCM based on RTE UPCM based on RTE Ref. [44]

0.04 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 x

0.2

0.4

0.6

0.8

1.0

Fig. 16. Non-dimensional radiative heat flux distribution over the bottom wall for DCM and UPCM method, (a) κ a ¼ 1 m  1 , and (b) κ a ¼ 0:1 m  1 .

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

collocation meshless method based on the first-order RTE and MSORTE. The UPCM based on RTE is demonstrated to be numerically stable and accurate in solving radiative transfer in strongly inhomogeneous media, media with void region and even with discontinuous extinction coefficient. The UPCM approach based on the RTE shows a higher accuracy and stability than DCM based on RTE. In some special cases, the numerical performance of the UPCM method based on the first-order RTE is shown to be even better than that of the DCM method based on MSORTE. Under the same node system, the computation time for UPCM is very close to that for the DCM approach based on the first-order RTE and far less than that for the DCM approach based on MSORTE. Besides, the UPCM method can be easily used to treat complex geometry.

Acknowledgments This work was supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant no. 51121004) and the National Natural Science Foundation of China (Grant no. 51176040). A very special acknowledgment is made to the editors and referees who made important comments to improve this paper. References [1] Tan JY, Liu LH, Li BX. Least-squares collocation meshless approach for coupled radiative and conductive heat transfer. Numer Heat Transfer B Fundam 2006;49:179–95. [2] Zhao JM, Tan JY, Liu LH. A second order radiative transfer equation and its solution by meshless method with application to strongly inhomogeneous media. J Comput Phys 2013;232:431–55. [3] Wallace JR, Elmquist KA, Haines EA. A ray tracing algorithm for progressive radiosity. Comput Gr 1989;23:315–24. [4] Modest MF. Radiative heat transfer. New York: McGraw-Hill; 1993. [5] Siegel R, Howell JR. Thermal radiation heat transfer. 4th ed. London: Taylor & Francis; 2001. [6] Yuen WW. The multiple absorption coefficient zonal method (MACZM), an efficient computational approach for the analysis of radiative heat transfer in multidimensional inhomogeneous nongray media. Numer Heat Transfer B Fundam 2006;49:89–103. [7] Zhou HC, Cheng Q, Huang ZF, He C. The influence of anisotropic scattering on the radiative intensity in a gray, plane-parallel medium calculated by the DRESOR method. J Quant Spectrosc Radiat Transfer 2007;104(1):99–115. [8] Howell JR, Perlmutter M. Monte Carlo solution of thermal transfer through radiant media between gray walls. ASME J Heat Transfer 1964;86:116–22. [9] Guo ZX, Aber J, Garetz BA, Kumar S. Monte Carlo simulation and experiments of pulsed radiative transfer. J Quant Spectrosc Radiat Transfer 2002;73:159–68. [10] Wang AQ, Modest MF, Haworth DC, Wang L. Monte Carlo simulation of radiative heat transfer and turbulence interactions in methane/air jet flames. J Quant Spectrosc Radiat Transfer 2008;109:269–79. [11] Krishna NA, Mishra SC. Discrete transfer method applied to radiative transfer in a variable refractive index semitransparent medium. J Quant Spectrosc Radiat Transfer 2006;102:432–40. [12] Rath P, Mishra SC, Mahanta P, Saha UK, Mitra K. Discrete transfer method applied to transient radiative transfer problems in participating medium. Numer Heat Transfer A Appl 2003;44:183–97. [13] Mishra SC, Lankadasu A. Transient conduction-radiation heat transfer in participating media using the lattice Boltzmann method and the discrete transfer method. Numer Heat Transfer A Appl 2005;47: 935–54.

79

[14] Mengüç MP, Viskanta R. Radiative transfer in three-dimensional rectangular enclosures containing inhomogeneous, anisotropically scattering media. J Quant Spectrosc Radiat Transfer 1985;33:533–49. [15] Mengüç MP, Iyer RK. Modeling of radiative transfer using multiple spherical harmonics approximations. J Quant Spectrosc Radiat Transfer 1988;39:445–61. [16] Sakami M, Charette A, Le Dez V. Application of the discrete ordinates method to combined conductive and radiative heat transfer in a two-dimensional complex geometry. J Quant Spectrosc Radiat Transfer 1996;56:517–33. [17] Ismail KAR, Salinas CS. Application of multidimensional scheme and the discrete ordinate method to radiative heat transfer in a twodimensional enclosure with diffusely emitting and reflecting boundary walls. J Quant Spectrosc Radiat Transfer 2004;88:407–22. [18] Kim C, Kim MY, Yu MJ, Mishra SC. Unstructured polygonal finitevolume solutions of radiative heat transfer in a complex axisymmetric enclosure. Numer Heat Transfer B Fundam 2010;57:227–39. [19] Lygidakis GN, Nikolos IK. Using the finite-volume method and hybrid unstructured meshes to compute radiative heat transfer in 3-D geometries. Numer Heat Transfer B Fundam 2012;62:289–314. [20] Liu LH, Zhang L, Tan HP. Finite element method for radiation heat transfer in multi-dimensional graded index medium. J Quant Spectrosc Radiat Transfer 2006;97:436–45. [21] Grissa H, Askri F, Ben Salah M, Ben Nasrallah S. Prediction of radiative heat transfer in 3D complex geometries using the unstructured control volume finite element method. J Quant Spectrosc Radiat Transfer 2010;111:144–54. [22] Li BW, Tian S, Sun YS, Hu ZM. Schur-decomposition for 3D matrix equations and its application in solving radiative discrete ordinates equations discretized by Chebyshev collocation spectral method. J Comput Phys 2010;229:1198–212. [23] Sun YS, Ma J, Li BW. Chebyshev collocation spectral method for three-dimensional transient coupled radiative-conductive heat transfer. ASME J Heat Transfer 2012;134:092701. [24] Sun YS, Li BW. Prediction of radiative heat transfer in 2D irregular geometries using the collocation spectral method based on bodyfitted coordinates. J Quant Spectrosc Radiat Transfer 2012;113: 2205–12. [25] Liu LH, Tan JY, Li BX. Meshless approach for coupled radiative and conductive heat transfer in one-dimensional graded index medium. J Quant Spectrosc Radiat Transfer 2006;101:237–48. [26] Tan JY, Zhao JM, Liu LH, Wang YY. Comparative study on accuracy and solution cost of the first/second-order radiative transfer equations using the meshless method. Numer Heat Transfer B Fundam 2009;55:324–37. [27] Wang CA, Sadat H, Ledez V, Lemonnier D. Meshless method for solving radiative transfer problems in complex two-dimensional and three-dimensional geometries. Int J Therm Sci 2010;49: 2282–8. [28] Fiveland WA, Jessee JP. Comparison of discrete ordinates formulations for radiative heat transfer in multidimensional geometries. J Thermophys Heat Transfer 1995;9:47–54. [29] Liu J, Chen YS. Examination of conventional and even-parity formulations of discrete ordinates method in a body-fitted coordinate system. J Quant Spectrosc Radiat Transfer 1999;61:417–31. [30] Sadat H. On the use of a meshless method for solving radiative transfer with the discrete ordinates formulations. J Quant Spectrosc Radiat Transfer 2006;101:263–8. [31] Hassanzadeh P, Raithby GD. Finite-volume solution of the secondorder radiative transfer equation: accuracy and solution cost. Numer Heat Transfer B Fundam 2008;53:374–82. [32] Zhao JM, Liu LH. Second-order radiative transfer equation and its properties of numerical solution using the finite-element method. Numer Heat Transfer B Fundam 2007;51:391–409. [33] Hughes TJR, Brooks A. A multidimensional upwind scheme with no crosswind diffusion. Finite element methods for convection dominated flows, AMD 1979;34:19–35. [34] Brooks AN, Hughes TJR. Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput Meth Appl Mech Eng 1982;32:199–259. [35] Liu LH. Finite Element simulation of radiative heat transfer in absorbing and scattering media. J Thermophys Heat Transfer 2004;18:555–7. [36] Zhao JM, Tan JY, Liu LH. A deficiency problem of the least squares finite element method for solving radiative transfer in strongly inhomogeneous media. J Quant Spectrosc Radiat Transfer 2012;113: 1488–502.

80

K. Luo et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 135 (2014) 66–80

[37] Lin H, Atluri SN. Meshless local Petrov–Galerkin (MLPG) method for convection diffusion problems. Comput Modell Eng Sci (CMES) 2000;1(2):45–60. [38] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992;10:307–18. [39] Wendland H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv Comput Math 1995;4:389–96. [40] Liu LH. Least-squares finite element method for radiation heat transfer in graded index medium. J Quant Spectrosc Radiat Transfer 2007;103:536–44. [41] Zhao JM, Liu LH. Discontinuous spectral element method for solving radiative heat transfer in multidimensional semitransparent media. J Quant Spectrosc Radiat Transfer 2007;107:1–16. [42] Crosbie AL, Schrenker RG. Radiative transfer in a two-dimensional rectangular medium exposed to diffuse radiation. J Quant Spectrosc Radiat Transfer 1984;31(4):339–72.

[43] Zhao JM, Liu LH. Spectral element method with adaptive artificial diffusion for solving the radiative transfer equation. Numer Heat Transfer B Fundam 2008;53(6):536–54. [44] Byun DY, Baek SW, Kim MY. Investigation of radiative heat transfer in complex geometries using blocked-off, multiblock, and embedded boundary treatments. Numer Heat Transfer A Appl 2003;43:807–25. [45] Amiri H, Mansouri SH, Safavinejad A. Combined conductive and radiative heat transfer in an anisotropic scattering participating medium with irregular geometries. Int J Therm Sci 2010;49: 492–503. [46] Zhang Y, Yi HL, Tan HP. Natural element method for radiative heat transfer in a semitransparent medium with irregular geometries. J Comput Phys 2013;241:18–34.