Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer

Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

Contents lists available at ScienceDirect

3 5

Journal of Quantitative Spectroscopy & Radiative Transfer

7

journal homepage: www.elsevier.com/locate/jqsrt

9 11

Q2

13

Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer

15 Q1

Yujia Sun a,b, Xiaobing Zhang a, John R. Howell b

17 Q3

a

Q4

b

School of Energy and Power Engineering, Nanjing University of Science and Technology, PR China Department of Mechanical Engineering, The University of Texas at Austin, United States

19 21 23 25 Q5

a r t i c l e i n f o

abstract

Article history: Received 2 May 2016 Received in revised form 28 June 2016 Accepted 31 July 2016

This work investigates the performance of P1 method, FVM and SP3 method for 2D combined conduction and radiation heat transfer. Results based on the Monte Carlo method coupled with the energy equation are used as the benchmark solution. Effects of the conduction-radiation parameter and optical thickness are considered. Performance analyses in term of the accuracy of heat flux and temperature predictions and of computing time are presented and analyzed. & 2016 Elsevier Ltd. All rights reserved.

27 29 31

63

33

1. Introduction

35

Radiation heat transfer is important in many industrial devices, especially those with high temperature. It is important to both the evolution of the temperature field (through the volumetric radiative heat source) and the wall total heat flux. An overall thermal analysis often needs to take conduction, convection and radiation into consideration. Unlike conduction and convection, the radiative transfer equation (RTE) is an integro-differential equation, which depends on two angular dimensions, and this makes it difficult to solve. Even in the absence of scattering, combined heat transfer modes with radiation are sometimes difficult to converge due to the high nonlinearity that radiation brings. Given this condition, only numerical methods are applicable for real problems. The purpose to this work is to provide insight into the factors that are necessary to choose an appropriate RTE solver for a combined model problem. Three RTE solvers will be used in a simple conduction problem to point out the important things to consider. Many methods have been developed to solve the RTE, such as the Monte Carlo Method (MCM) [1–3], the zonal method [4], the diffusion approximation [5], the spherical harmonics method (PN) [6,7], the discrete transfer method (DTM) [8,9], the discrete ordinates method (DOM, or SN)

37 39 41 43 45 47 49 51 53 55 57 59

[10–13], the finite volume method (FVM) [14,15], the finite element method (FEM) [16,17], the meshless method [18], and so on. Among them, PN (e.g. P1), DOM and FVM have gained popularity and are implemented into general purpose CFD software, such as OpenFOAM and ANSYS Fluent, while they leave the choice of method up to the user. Tencer and Howell [19] performed a parametric study of the accuracy of M1 method, P1 method and SN method for several 2D benchmark problems. It is shown that SN method suffers from ray effects for Case 1 and provides very accurate results for Case 3 and the P1 method error is less than the M1 method while larger than the SN method for these two cases for parameters considered in that work. Numerical simulations of combined problems often start from zero or uniform internal temperature with different wall temperature (discontinuous or continuous). Thus different RTE solvers coupled with the energy equation iteration may lead to quite different convergence behavior, that is, combined problem convergence may rely largely on the performance of RTE solvers. Yuen and Takara [20] used the generalized exponential integral function to solve two dimensional combined conductive–radiative heat transfer. They concluded that the additive solution is an acceptable approach in estimating heat transfer and the diffusion approximation yields significant errors in both the temperature and heat flux

http://dx.doi.org/10.1016/j.jqsrt.2016.07.024 0022-4073/& 2016 Elsevier Ltd. All rights reserved.

61 Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

65 67 69 71 73 75 77 79 81 83 85 87 89 91

2

1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 53 55 57 59 61

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

predictions. Kim and Baek [21] studied the same problem using the central difference scheme for the heat diffusion equation and the discrete ordinates method for the RTE. They also included anisotropic scattering. Wu and Ou [22] took scattering into consideration for transient twodimensional combined problems and investigated the influence of aspect ratio, scattering albedo and conductionradiation parameters and concluded that the modified differential approximation is superior to the P1 method. The largest optical thickness considered was 5. Mishra et al. [23] investigated the computational efficiency of the collapsed dimension method and discrete transfer method for several radiation and combined mode heat transfer examples. These two methods were found to give the same results while the former was faster. Mishra [24] compared DTM, DOM and FVM methods for transient combined conduction and radiation heat transfer and found that DTM is the most time consuming while DOM is most efficient. Only an optical thickness of unity was considered. Wang et al. [25] developed a numerical iteration technique based on the flux conservation equation to solve 1D coupled conductionradiation problem and showed that it converges as rapidly as Thomas's method. Naraghi and Saltiel [26] compared the performance of the conventional successive substitution method, accelerated fixed points, and Newton–Raphson methods, and they considered the effect of variable conductivity, while the optical thickness was unity. Recently, there are papers dealing with combinations of different energy solvers and RTE solvers for these combined problems, such as the spectral collocation method [27,28], natural element method [29,30], meshless method [31,32], lattice Boltzmann method (LBM) coupled with DOM and FVM [33– 36], and LBM for both the energy equation and RTE [37]. These works mainly focus on developing new solvers for combined problems and how to choose an appropriate method for combined conduction and radiation heat transfer remains to be investigated. Although comparisons of different RTE solvers exist for a wide range of parameters [19], Tencer [38] analyzed the advantages and disadvantages of implementing different RTE solvers into conjugate heat transfer codes, while there is little work reporting their performance for combined heat transfer [38,39]. Dombrovsky [40,41] analyzed the choice of method for RTE when considering inhomogeneous properties and scattering. The choice of method may depend on accuracy and computing time, which are functions of many parameters, e.g., optical thickness, conduction-radiation number, scattering albedo, spectral properties. This kind of problem has not been widely investigated and this motivates us to do this work. This work aims to investigate the performance of P1 method, FVM and SP3 method for 2D combined conduction and radiation heat transfer. The P1 method and FVM are widely used for practical applications. The SP3 method is chosen in this work as it shows a distinct improvement over the P1 method with reasonable computing time [42]. Results obtained by MCM coupled with the energy equation serve as the benchmark solution, although MCM may introduce convergence instabilities considering its stochastic nature. It produced very smooth results for most of the cases studied. The results of accuracy in terms of heat

flux, temperature, and computing time will be presented and analyzed in detail.

63 65

2. Mathematical formulations

67

2.1. Energy equation

69

For steady combined conduction and radiation heat transfer problems, the energy equation is

71

∇ðk∇TÞ ¼ ∇ U qr

ð1Þ

73

where k is the thermal conductivity, the left-hand side accounts for conduction, and the right-hand side accounts for radiation, which makes the combined heat transfer complicated. The local divergence of radiative heat flux can be calculated as:

75

∇ U qr ¼ κð4σT 4 GÞ

81

ð2Þ

where κ is the absorption coefficient of the participating medium, σ is the Stefan-Boltzmann constant, and G is the incident radiation, which implicitly depends on the temperature field. This divergence of radiative heat flux (or the radiative heat source), specifically, G is obtained by MCM, P1 method, FVM, and SP3 method in this work. Existence of fourth power in the radiation term introduces high nonlinearity to the original simple diffusion equation. Incorporating the radiative heat source is usually done by loose coupling, in which ∇ U qr is computed based on initial conditions (basically temperature) by one of the RTE solvers. The radiative heat source is then inserted into the energy equation and a new temperature field is obtained, and the process is repeated until convergence. This separate solution procedure (or successive substitution) makes both the RTE and energy equation linear or slightly nonlinear and are commonly used in combined convection, combustion problems, in which the energy equation is more complicated.

77 79

83 85 87 89 91 93 95 97 99 101 103

2.2. P1 formulation

105

The P1 method is the simplest and most widely used RTE solver, and it consists of only one elliptic equation in terms of incident radiation, G
 1 1 ∇ U ∇G  G ¼  4πI b ð3Þ 3κ κ subject to the Marshak boundary condition: 

2ε 2 n U ∇G þ G ¼ 4πI bw ε 3κ 1 ∇G 3κ

109 111 113

ð4Þ

and the radiative heat flux is calculated by: qr ¼ 

107

115 117

ð5Þ

The P1 method is known to be efficient for optical thick conditions. For optical thin conditions, it not only converges very slowly but also yield results at poor accuracy. Despite this, it is still widely used due to its simplicity.

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

119 121 123

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

2.3. SP3 formulation

3

Modest's [42] formulation of the simplified spherical harmonics methods for RTE is chosen and only the SP3 method is considered here. The SP3 model has two simultaneous equations and two boundary conditions: Governing equations:

 1 1 2 ∇τ U ∇τ J 0 ¼ α0 J 0  J 2  I b 3 α1 3

 ð6Þ 3 1 5 4 ∇τ U α2 þ α0 J 2 2α0 ð J 0 I b Þ ∇τ J 2 ¼ 7 α3 3 3

5 7 9 11 13 15 17

Boundary conditions:   Jw 1 1 1 3α1 n U ∇τ J 0 ¼ 2 J 0  π  8 J 2   Jw 1 1 7 7α3 n U ∇τ J 2 ¼  8 J 0  π þ 24 J 2

ð7Þ

25

So the SP3 method has two elliptic partial differential equations, which are connected by the source term and boundary conditions. Once the J k have been determined, incident radiation and radiative flux are obtained from
 2 G ¼ 4π J 0  J 2 ð8Þ 3

27

qr ¼ 

19 21 23

4π ∇τ J 0 3α1

ð9Þ

2.4. FVM formulation 31 Followed Chai [43], the radiative transfer equation (RTE) in any direction b s is: 4

35

dI σT ¼  κI þ κ ds π

37

where κ is the absorption coefficient. It can be rewritten as:

39

dI ¼  κI þ S ds

41 43

ð10Þ

ð11Þ

where S is the source term given by σT 4 S¼κ π

ð12Þ

Integrating Eq. (11) over the solid angle, 45 47 49 51 53

are average face areas. For 2D cases, when marching from the first quadrant, the west and south intensities of the boundary are known values. Once I m P has been calculated, m Im E and I N are readily determined from the relation Eq. (15), and are set to be the west and south intensities for the next control volumes. The procedure is then repeated three times, starting from the other corners of the enclosure. The scattering is not considered in this work, but the scattering term is kept here for completeness of FVM. The radiative heat source ∇ U qr in the energy equation is: ! σT 4 G ð17Þ ∇ U qr ¼ κ 4π π where G is the incident radiation: Z 4π G¼ IdΩ 0

∂I m m ∂I m m ∂I m m D þ D þ D ¼  βI m ΔΩm þSm ΔΩm ∂x x ∂y y ∂z z where Dm is given by Z m b Ub Dm ¼ ðn s ÞdΩ ΔΩm

ð13Þ

ð14Þ

After some mathematical manipulations, the following is obtained Dm Dm A Dm m m x AEW m x AFB m I W þ xγ NS I m S þ γ x I B þ ΔVΔΩ SP γx x m m Dm A D A D A m x E x N x F γ x þ γ x þ γ x þ βΔVΔΩ

55

Im P ¼

57

where γ are the spatial weights and

59

AEW ¼ ð1  γ x ÞAE þ γ x AW ANS ¼ ð1  γ y ÞAN þ γ y AS

61

AFB ¼ ð1  γ z ÞAF þ γ z AB

ð15Þ

ð16Þ

63 65 67 69 71 73 75 77

ð18Þ

Once the intensity distribution is known, the heat flux can also be obtained: Z 4π qr ¼ I cos θ dΩ ð19Þ 0

2.5. Monte Carlo method

29

33

3

The forward ray-tracing Monte Carlo method is used in this work to generate benchmark solutions for comparisons with other RTE solvers [44]. Given a uniform absorption coefficient, the READ method is adopted to reduce the computing time. In the READ method, calculation of the distribution of radiation energy absorption is outside of the energy iteration loop, thus reducing the entire computing time. A detailed formulation is available in [44] and not repeated here.

79 81 83 85 87 89 91 93 95 97 99

3. Numerical treatments

101

3.1. Problem description

103

The commonly used benchmark conduction-radiation 105 problem [16,20,21,26] is chosen in this work. It is a rectangular enclosure with black boundaries, absorbing, 107 emitting and non-scattering media. Temperatures at the four boundaries are prescribed, and the optical thickness 109 and the conduction-radiation parameter have a wide range of values. The dimensionless temperature is obtained by 111 dividing by the bottom wall temperature (Fig. 1). Q6 113 3.2. Finite volume discretization 115 The RTE, energy equation, P1 equation and SP3 equation are all discretized by finite volume method, which is commonly used for CFD codes, thus making this work valuable for applications with CFD solvers. MCM particle tracking is also based on the finite volume discretization of the computing domain. Detailed discretization procedures are not shown here. They are available in many CFD textbooks [45].

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

117 119 121 123

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

4

1 3 5 7 9 11 13 Fig. 1. Problem geometry.

15 17

3.3. Linearization of the radiative heat source

19

27

As stated in Section 2.1, loose coupling is used for the radiative heat source. Due to the high non-linearity of the radiative heat source, the radiative heat source needs to be linearized to make the energy iteration stable, otherwise it will diverge after several iterations, especially for large optical thickness. Following Patankar [46], the radiative heat source is split into two parts thanks to the explicit term of fourth power of temperature 4κσT 4

29

4κσT 4 ¼ 4κσT 4 þ16κσT 3 ðT T  Þ

21 23 25

31 33 35

¼ 16κσT 3 T  12κσT 4

ð20Þ

The first term is absorbed into the coefficients of nodal energy equation and the second term is used in the calculation of the radiative heat source. This technique is used in all the RTE solvers when coupled with the energy equation.

This numerical set is used for all cases considered in this work, from optically thin to optically thick conditions. On the one hand, there is no great improvement of refining the spatial and angular discretization, and computing time dramatically increases with refined mesh. There is always a trade-off between accuracy and computing time, especially for FVM for optically thick cases. On the other hand, for the energy equation, there is no need to use a very fine mesh and it is thus desirable to evaluate different RTE's performances under the same grid size. Additionally, this may also be valuable for cases with non-gray gas radiation. For example, in WSGG/FSK/SLW model, different gray gases have different orders of magnitude of absorption coefficient, which ranges from 0 m  1 to several hundreds of m  1. For the clear gas with absorption coefficient 0 m  1, although it does not contribute to the total radiative heat source, it may have prominent effect on the wall radiative heat flux for some problems, like the one considered in this work, where the heat is mainly transferred from wall. One can skip the clear gas when solving the RTE during the energy iteration and solve it along with other gases for the wall heat flux when the iteration is converged. All these gray gases are modeled by an RTE solver with the same numerical and thermophysical parameters.

39

45 47 49 51 53 55 57 59 61

67 69 71 73 75 77 79 81 83 85

4. Results

89

Temperature distributions for different cases and are analyzed first and then detailed error patterns in terms of temperature and radiative heat flux for three methods are presented. After these, performances comparisons among these methods are discussed.

91

4.1. Temperature distributions for five optical thickness and four conduction-radiation numbers

93 95 97 99

3.4. Validations

43

65

87

37

41

63

To validate the current codes, MCM results are compared with available literature [20], which are commonly used as benchmark solutions. Fig. 2 shows the temperature profile along the vertical center line for τ ¼ 5 for four different conduction-radiation parameters. Agreement with the reference validates our code. The results for small optical thickness are not shown here as they are also almost identical. For other RTE solvers, namely, the P1 method, SP3 method and FVM, grid independence tests are performed to determine the optimal spatial discretization. Fig. 3 shows the effects of grid size and angular discretization on the temperature profile. The first number in the legend means grid size, the second the angular number in the polar direction, and the third in the azimuthal direction. With a large enough angular number, the effect of grid size is small and 43 is chosen here because the temperature profile of case 23 has relatively large derivations in regions near the hot wall and cold wall. The right-hand figure represents the effect of angular discretization on the temperature. The 2–4 case has large discrepancy with the two other cases and the difference between 4–8 and 4–12 is quite small, the 4–8 is chosen for further study.

Fig. 4 shows the temperature distributions for different optical thickness and four conduction-radiation numbers, ranging from optically thin to optically thick, and radiation dominated to conduction dominated. Figures in the same row have the same optical thickness and columns have the same conduction-radiation number. Numerical values for the graphs will be posted on-line as appendices for this paper. For small optical thickness, τ ¼0.1, radiation dominates for small Nc¼0.001. Differences between Nc¼0.01, Nc¼0.1, and Nc¼1.0 are relatively small, because conduction is the main heat transfer mode. The radiative heat source is small due to the optically thin condition, and this leads to small contribution of radiation exchange to the energy equations. When the optical thickness increases to 1.0, the temperatures become higher in most of the domain compared to case of τ ¼ 0.1, and the contour line is more uniform. Large temperature gradients occur at the four walls, indicating high non-linearity for this case. The effect of radiation on the temperature distribution decreases with the increase of Nc, which is similar to τ ¼0.1. And more importantly, the temperature become more linear with larger Nc. With further increasing optical thickness, τ ¼5.0 and τ ¼10.0, the same trend occurs for τ ¼1.0, although the temperature is more linear along the three cold walls. The

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

101 103 105 107 109 111 113 115 117 119 121 123

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1 3 5 7 9 11 13 15 17 19 21 23 25

effects of radiation are more prominent for larger Nc, which can be clearly seen for Nc¼0.1, in which the temperature contour is more spread. For very large optical thickness, τ ¼50, even for small Nc, the temperature shows a more uniform pattern. For example, the isotherm of 950 K covers a relatively larger region near the hot wall compared to other small optical thickness cases. The reason for this phenomenon is due to the diffusion-like behavior of radiation for large optical thickness. In other words, radiation is mainly exchanged only between neighboring regions. To sum up, radiation plays an important role for small Nc for all ranges of optical thickness considered, and shows highly non-linear behavior for small optical thickness and diffusion behavior for large optical thickness. Large optical thickness has a significant effect even for large Nc.

4.2. Mean relative error of temperature along the vertical center line for different methods In this section, the temperature errors for the three different methods are compared and analyzed. The results are based on the mean relative error of temperature along the vertical center line.

27 29 31 33

5

Fig. 5 shows relative error as a function of optical thickness and conduction-radiation number. Axes have a log scale. Fig. 5(a) shows the result for the P1 method. The maximum error occurs for case of τ ¼1.0, Nc ¼0.001, and is about 1.55%. The P1 method is less accurate for cases with small optical thickness and conduction-radiation number. The reason for its high accuracy is due to small contribution of radiation to energy transfer. For optical thickness larger than 10, the P1 method is less sensitive to the conduction-radiation number. For range of 1–10, its accuracy increases with the decrease of conduction-radiation number. The SP3 method shows superior accuracy to the P1 method. It has the largest error for the case τ ¼0.1, Nc¼0.001, and performs quite well for other cases. Its error is also small and not sensitive to N number for large optical thickness. Fig. 5(c) shows the error of the FVM. The angular discretization is 4 in polar angle and 8 in azimuthal angle. Further refining does not improve the results significantly and greatly increases the computing time. The FVM shows different error contours for the range of parameters considered. It has the largest errors for case with τ ¼1.0, Nc¼0.001, as for the P1 method, while keeping high accuracy for smaller optical thickness. For τ ¼5.0 and 10.0, the FVM has good performance and its accuracy changes within a small range for different Nc. For τ ¼ 50.0, it shows an opposite error dependence on Nc. Its error increases with the increase of Nc, and is largest for Nc ¼1.0. All the simulations for the three methods are performed with the same grid increments to allow direct comparison.

63

4.3. Mean relative error of different methods in term of radiative heat flux along the bottom wall

95

65 67 69 71 73 75 77 79 81 83 85 87 89 91 93

97

35 37 39 41 Fig. 2. Validation of MCM results.

In this section, the mean error in bottom wall radiative heat flux of the three methods are compared. Only radiative heat flux is considered here because the conductive heat flux is proportional to the temperature difference, and hence should share the same error pattern. Fig. 6(a) and (b) show the results for the P1 and SP3 methods. As seen, both these methods have large errors

99 101 103

43

105

45

107

47

109

49

111

51

113

53

115

55

117

57

119

59

121

61

Fig. 3. Grid Independence Tests for FVM for τ ¼ 5; Nc ¼ 0:01.

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

123

6

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

63

3

65

5

67

7

69

9

71

11

73

13

75

15

77

17

79

19

81

21

83

23

85

25

87

27

89

29

91

31

93

33

95

35

97

37

99

39

101

41

103

43

105

45

107

47

109

49

111

51

113

53

115 117

55 Fig. 4. Temperature contours for all cases.

57

119

59

121

61

123 Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

7

1

63

3

65

5

67

7

69

9

71

11

73

13

75

15

77

17

79

19

81

21

83

23

85

25

87

27

89

29

91

31

93

33

Fig. 5. Temperature error variations with τ and N.

97

35 37 39 41 43 45 47 49

for optical thickness smaller than τ ¼10.0, reaching 50%, and have acceptable errors for optically thick conditions, less than 10%. Their accuracies increase with the increase of optical thickness, and show weak dependence on conduction-radiation number. Contrary to the P1 and SP3 methods, the relative error of radiative heat flux for the FVM increases with the increase of optical thickness. The FVM has high accuracy for small and intermediate optical thickness (o¼10), while it performs poorly for τ ¼50. The radiative heat flux error reaches 60% irrespective of Nc. From the above analyses, the mean relative errors in terms of radiative heat flux of these three methods mainly depends on optical thickness regardless of the value of Nc. 4.4. Computing time patterns for three methods Besides the relative errors of temperature and heat flux, computing time is also important in practical applications and it is desirable that the RTE solver does not consume too much time, otherwise it will slow down the overall energy iterations. The above figures show the computing time contours of the three methods. For the P1 method and FVM, both Gauss-Seidel iteration (GS) method and ADI method are adopted. As seen from Fig. 7(a) and (b), the P1

As seen from Fig. 8, for small optical thickness, τ ¼1, the FVM is shown to be the fastest method. The SP3 method is much slower than the other two methods. Although there are differences in CPU time, their errors in temperature are both small, with SP3 the most accurate. Both P1 method and SP3 method failed to accurately predict the radiative heat flux along the bottom wall. So, in this condition, FVM is the best method.

53

57 59 61

method needs a very long time to converge for small optical thickness, no matter whether using GS or ADI, although the latter is faster than the former. The SP3 method has a similar pattern to the P1 method. For these two solvers, their worst performance occurs for the cases with large Nc (4 ¼0.01). Different from the P1 and SP3 methods, the FVM converges much faster in small and intermediate optical thickness, less than 8.45 s, while its computing time rises to more than 30 s for the case of τ ¼50.0. The mesh size for different optical thickness and methods are the same to allow for direct comparisons. 4.5. Comparisons of three methods for CPU time and relative error for different optical thickness

51

55

95

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

99 101 103 105 107 109 111 113 115 117 119 121 123

8

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

63

3

65

5

67

7

69

9

71

11

73

13

75

15

77

17

79

19

81

21

83

23

85

25

87

27

89

29

91

31

93 95

33 Fig. 6. Radiative heat flux errors for different τ and N.

97

35

99

37 39 41 43 45 47 49 51 53 55 57 59 61

For optical thickness¼ 5.0, the SP3 method is still the most time-consuming method, nearly 4–8 times slower than the other two methods. The mean error of temperature for all three methods is under 1.0%, and the SP3 method is slightly more accurate. In terms of wall radiative heat flux, the FVM remains very accurate in this condition. Although errors of the P1 method and SP3 methods are larger than FVM, their values decreased from 25% to 10%, compared to τ ¼ 1.0. In this scenario, FVM is the most appropriate method. The P1 method may also be used if wall radiative heat flux is not desired. Although this case can be considered as optical thick, the P1 method is not expected to give good results. For an optically thicker condition, τ ¼10.0, a different performance is observed. The P1 method becomes the fastest method, with the SP3 method still the slowest. And both the P1 and SP3 methods converge faster compared to τ ¼5.0. All three methods perform very well in terms of temperature and the SP3 method is the most accurate. For radiative heat flux, the P1 method has error of about 5.5%, the SP3 method about 4.3%, and FVM 6.5%. Although SP3 method is very accurate for both temperature and heat flux, its CPU time is relatively longer than the two other methods. Due to the small errors, both P1 method and

FVM can be used in this condition, and the P1 method is the best. For very large optical thickness, the FVM does not perform as well as P1 and SP3 methods. The P1 method is faster than FVM and SP3 method, which are similar. They are about 5 times slower than the P1 method. For temperature, the error of FVM increases with increases in Nc, and the P1 and SP3 methods are very accurate. Unexpectedly, P1 and SP3 methods perform very well in terms of radiative heat flux, whose errors are both about 5%, while FVM fails to get correct radiative heat flux values. Its error reaches 50% irrespective of Nc. For this condition, the P1 method is the ideal method. Based on the above analyses, a table is made to help researchers to choose among different solvers for the range of parameters considered in this paper. Although smaller optical thicknesses are not shown in this paper, FVM seems to be the appropriate one. For example, at very small optical thickness (such as τ ¼0.01), results show that radiation has very small effect on the temperature distribution, so the QoI in this scenario must be radiative heat flux, which the P1 method fails to predict correctly.

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

101 103 105 107 109 111 113 115 117 119 121 123

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

9

1

63

3

65

5

67

7

69

9

71

11

73

13

75

15

77

17

79

19

81

21

83

23

85

25

87

27

89

29

91

31

93

33

95

35

97

37

99

39

101

41

103

43

105

45

107

47

109

49

111

51 53

Fig. 7. Computing time of different methods for different τ and N.

4.6. Discussions

55 57 59 61

1. Table 1 is based on current simulations for a 2D enclosure. When applied to 3D problems, their accuracies are expected to be similar to 2D problems, while their computing time may be different, especially for optically intermediate cases. For optical thick (τ 4 ¼10) conditions, the P1 method may be still better than FVM.

Due to poor accuracy of P1 method for optically thin conditions, FVM is better in this range. For optically intermediate cases, FVM is expected to perform better than P1. Extensions to problems with scattering will be not that simple. Firstly, the scattering albedo is another important parameter that controls the conduction and radiation interaction. Secondly, convergence is slower with scattering, especially for the FVM. In this case, the algebraic

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

113 115 117 119 121 123

10

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

63

3

65

5

67

7

69

9

71

11

73

13

75

15

77

17

79

19

81

21

83

23

85

25

87

27

89

29

91

31

93

33

95

35

97

37

99

39

101

41

103

43

105

45

107

47

109

49

111

51

113

53

115 117

55 Fig. 8. Comparisons of computing time of different methods.

119

57 59 61

equation for each direction is coupled and hence needs more iteration to converge. This aspect needs more investigation.

2. Table 1 shows that FVM is superior to both the P1 and SP3 methods for optical thickness smaller than 5, while the latter are better for τ 4 ¼10. This may also apply to

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

121 123

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎

1

11

63

Table 1 Choosing table of RTE solvers for different parameters.

3

Nc

T

Q

Both

65

5

τ ¼ 0.1 τ¼1

0.001, 0.01, 0.1, 1 0.001, 0.01, 0.1, 1

FVM 4P1 4SP3 FVM 4P1 4SP3

FVM4 P1 4SP3 FVM4 P1 4SP3

FVM 4P1 4SP3 FVM 4P1 4SP3

67

7

τ¼5 τ ¼ 10 τ ¼ 50

FVM 4P1 4SP3 P1 4FVM 4SP3 P1 4FVM 4SP3 P1 4SP3 4FVM

FVM4 P1 4SP3 FVM4 P1 4SP3 P1 4FVM 4SP3 P1 4SP3 4FVM

FVM 4P1 4SP3 FVM 4P1 4SP3 P1 4FVM 4SP3 P1 4SP3 4FVM

69

9

0.001 0.01, 0.1, 0.1 0.001, 0.01, 0.1, 1 0.001, 0.01, 0.1, 1

73

11 13

different gray gases when using SLW/FSK. That is, the FVM method is expected to perform well for the gray gases with small and intermediate optical thickness while P1 method performs well for those with large optical thickness. It may be beneficial to use a hybrid method of the FVM with the P1 for different gray gases. Their performances for the overall heat source and radiative heat flux may not be easily determined from the conclusion of this work because the contribution of each gray gas is different and also problem dependent.

15 17 19 21 23 25

5. Conclusions

27

This work performs parametric studies of the P1 method, SP3 method and FVM for combined conduction and radiation heat transfer problems. Their accuracies and computing time are analyzed by comparing with MCM benchmark results. The P1 method converges slowly for small optical thickness and yields large errors in wall radiative heat flux. The SP3 method is the most timeconsuming method, while only improving the accuracy of heat flux slightly. The FVM performs very well for τ o ¼5.0 in both accuracy and computing time, but fails to accurately predict the radiative heat flux for large optical thickness, especially τ ¼50. For the cases considered in this work, FVM is recommended for τ o ¼5.0, and P1 for τ ¼10, 50.

29 31 33 35 37 39 41 43 45

Acknowledgments

This work is supported by the Natural Science Foundation of Jiangsu Province (Grant no. BK20131348) and Key 47 Laboratory Foundation of the People's Republic of China 49 Q8 (Grant no. 9140C300206120C30110) and the China Scholarship Council. Q7

51 53 55 57 59 61

71

References [1] Howell JR. The Monte Carlo method in radiative heat transfer. J Heat Transf 1998;120:547. http://dx.doi.org/10.1115/1.2824310. [2] Modest MF. Backward Monte Carlo simulations in radiative heat transfer. J Heat Transf 2003;125:57. http://dx.doi.org/10.1115/1.1518491. [3] Farmer JT, Howell JR. Comparison of Monte Carlo strategies for radiative transfer in participating media. Adv Heat Transf 1998;31: 333–429. http://dx.doi.org/10.1016/S0065-2717(08)70243-0. [4] M.F. Modest, The zonal method, in: Radiat. Heat Transf. Second Ed., 2003: p. 539–64. 〈http://dx.doi.org/10.1016/B978-012503163-9/50018-7〉.

[5] Deissler RG. Diffusion approximation for thermal radiation in gases with jump boundary condition. J Heat Transf 1964;86:240–5. [6] Ratzel AC, Howell JR. Two-dimensional radiation in absorbingemitting media using the p-n approximation. J Heat Transf 1983;105:333–40. http://dx.doi.org/10.1115/1.3245583. [7] Yang J, Modest MF. Elliptic PDE formulation of general, three-dimensional high-order PN-approximations for radiative transfer. J Quant Spectrosc Radiat Transf 2007;104: 217–27. http://dx.doi.org/10.1016/j.jqsrt.2006.07.017. [8] Cumber PS. Improvements to the discrete transfer method of calculating radiative heat transfer. Int J Heat Mass Transf 1995;38: 2251–8. http://dx.doi.org/10.1016/0017-9310(94)00354-X. [9] Coelho PJ, Carvalho MG. A conservative formulation of the discrete transfer method. J Heat Transf 1997;119:118–28. http://dx.doi.org/ 10.1115/1.2824076. [10] Fiveland WA. Discrete-ordinates solutions of the radiative transport equation for rectangular enclosures. J Heat Transf 1984;106: 699. http://dx.doi.org/10.1115/1.3246741. [11] Coelho PJ, High-Order Resolution Bounded Skew. Schemes for the discrete ordinates method. J Comput Phys 2002;175: 412–37. http://dx.doi.org/10.1006/jcph.2001.6899. [12] Chai JC, Lee HS, Patankar SV. Ray effect and false scattering in the discrete ordinates method. Numer Heat Transf Part B Fundam 1993;24:373–89. http://dx.doi.org/10.1080/10407799308955899. [13] Koch R, Becker R. Evaluation of quadrature schemes for the discrete ordinates method. J Quant Spectrosc Radiat Transf 2004;84: 423–35. http://dx.doi.org/10.1016/S0022-4073(03)00260-7. [14] Chai JC, Lee HS, V Patankar S. Finite volume method for radiation heat transfer. J Thermophys Heat Transf 1994;8:419–25. http://dx.doi.org/ 10.2514/3.559. [15] Raithby GD. Discussion of the finite-volume method for radiation, and its application using 3d unstructured meshes. Numer Heat Transf Part B Fundam 1999;35:389–405. http://dx.doi.org/10.1080/ 104077999275802. [16] Razzaque MM, Howell JR, Klein DE. Coupled radiative and conductive heat transfer in a two-dimensional rectangular enclosure with gray participating media using finite elements. J Heat Transf 1984;106:613. http://dx.doi.org/10.1115/1.3246723. [17] Liu LH, Zhang L, Tan HP. Finite element method for radiation heat transfer in multi-dimensional graded index medium. J Quant Spectrosc Radiat Transf 2006;97:436–45. http://dx.doi.org/10.1016/ j.jqsrt.2005.05.067. [18] Sadat H. On the use of a meshless method for solving radiative transfer with the discrete ordinates formulations. J Quant Spectrosc Radiat Transf 2006;101:263–8. http://dx.doi.org/10.1016/j.jqsrt. 2005.11.019. [19] Tencer, J, Howell, JR. A parametric study of the accuracy of several radiative transport solution methods for a set of 2-D benchmark problems. In: Proc. ASME 2013 Heat Transf. Summer Conf., Minneapolis, 2014: p. 1–10. [20] Yuen WW, Takara EE. Analysis of combined conductive-radiative heat transfer in a two-dimensional rectangular enclosure with a gray medium. J Heat Transf 1988;110:468–74. [21] Kim Taik Young, Baek Seung Wook. Analysis of combined conductive and radiative heat transfer in a two-dimensional rectangular enclosure using the discrete ordinates method. Int J Heat Mass Transf 1991;34:2265–73. http://dx.doi.org/10.1016/0017-9310(91) 90052-G. [22] Wu CY, Ou NR. Transient two-dimensional radiative and conductive heat transfer in a scattering medium. Int J Heat Mass Transf 1994;37: 2675–86. http://dx.doi.org/10.1016/0017-9310(94)90384-0. [23] Mishra SC, Talukdar P, Trimis D, Durst F. Computational efficiency improvements of the radiative transfer problems with or without conduction–a comparison of the collapsed dimension method and

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

75 77 79 81 83 85 87 89 91 93 95 97 99 101 103 105 107 109 111 113 115 117 119 121 123

12

1 3

[24]

5 7 9

[25]

[26]

11 [27]

13 [28]

15 17 [29]

19 21

[30]

23 [31]

25 27 29 31

[32]

[33]

[34]

Y. Sun et al. / Journal of Quantitative Spectroscopy & Radiative Transfer ∎ (∎∎∎∎) ∎∎∎–∎∎∎ the discrete transfer method. Int J Heat Mass Transf 2003;46: 3083–95. http://dx.doi.org/10.1016/S0017-9310(03)00075-9. Mishra SC, Kim MY, Maruyama S. Performance evaluation of four radiative transfer methods in solving multi-dimensional radiation and/or conduction heat transfer problems. Int J Heat Mass Transf 2012;55: 5819–35. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.05.078. Wang JY, Gao ZX, Lee CH. An iterative technique for coupled conduction–radiation heat transfer in semitransparent media. Numer Heat Transf Part A Appl 2015;67:1208–31. http://dx.doi.org/ 10.1080/10407782.2014.955370. Naraghi MHN, Saltiel CJ. Combined-Mode Heat Transfer in Radiatively Participating Media: Computational Considerations. Numer Heat Transf Part B Fundam 1992;21:253–67. Sun Y-S, Ma J, Li B-W. Chebyshev collocation spectral method for three-dimensional transient coupled radiative–conductive heat transfer. J Heat Transf 2012;134:92701. Sun YS, Li BW. Chebyshev collocation spectral approach for combined radiation and conduction heat transfer in one-dimensional semitransparent medium with graded index. Int J Heat Mass Transf 2010;53: 1491–7. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2009.11.047. Zhang Y, Ma Y, Yi H-L, Tan H-P. Natural element method for solving radiative transfer with or without conduction in three-dimensional complex geometries. J Quant Spectrosc Radiat Transf 2013;129: 118–30. Zhang Y, Yi HL, Tan HP. Natural element method analysis for coupled radiative and conductive heat transfer in semitransparent medium with irregular geometries. Int J Therm Sci 2014;76: 30–42. http://dx.doi.org/10.1016/j.ijthermalsci.2013.08.013. 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 Transf 2006;101:237–48. http://dx.doi.org/ 10.1016/j.jqsrt.2005.11.017. Liu LH, Tan JY. Meshless local Petrov–Galerkin approach for coupled radiative and conductive heat transfer. Int J Therm Sci 2007;46: 672–81. http://dx.doi.org/10.1016/j.ijthermalsci.2006.09.005. Sun Y, Zhang X. Analysis of transient conduction and radiation problems using lattice Boltzmann and finite volume methods. Int J Heat Mass Transf 2016;97:611–7. Sun Y, Zhang X. Analysis of transient conduction and radiation problems using the lattice Boltzmann and discrete ordinates

[35]

[36]

[37]

[38]

[39] [40]

[41] [42]

[43]

[44]

[45] [46]

methods. Numer. Heat Transf Part A Appl Int J Comput Methodol 2015;68:619–37. http://dx.doi.org/10.1080/10407782.2014.994406. Mondal B, Mishra SC. Lattice Boltzmann method applied to the solution of the energy equations of the transient conduction and radiation problems on non-uniform lattices. Int J Heat Mass Transf 2008;51: 68–82. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2007.04.030. Mishra SC, Roy HK. Solving transient conduction and radiation heat transfer problems using the lattice Boltzmann method and the finite volume method. J Comput Phys 2007;223:89–107. http://dx.doi.org/10. 1016/j.jcp.2006.08.021. Mishra SC, Poonia H, Das AK, Asinari P, Borchiellini R. Analysis of conduction-radiation heat transfer in a 2D enclosure using the lattice Boltzmann method. Numer Heat Transf Part A Appl 2014;66: 669–88. Tencer J, Howell JR. Coupling radiative heat transfer in participating media with other heat transfer modes. J Braz Soc Mech Sci Eng 2015. http://dx.doi.org/10.1007/s40430-015-0434-1. Howell JR, Personal View A. of 50 years of thermal radiation heat transfer research. Adv Heat Transf 2015;47:309–40. Dombrovsky LA. The use of transport approximation and diffusionbased models in radiative transfer calculations. Comput Therm Sci Int J 2012;4:297–315. Dombrovsky LA, Baillis D. Thermal radiation in disperse systems: an engineering approach title. New York: Begell House; 2010. Modest MF, Cai J, Ge W, Lee E. Elliptic formulation of the Simplified Spherical Harmonics Method in radiative heat transfer. Int J Heat Mass Transf 2014;76:459–66. http://dx.doi.org/10.1016/j. ijheatmasstransfer.2014.04.038. Chai JC, Lee HS, Patankar SV. Finite volume method for radiation heat transfer. J Thermophys Heat Transf 1994;8:419–25. http://dx.doi.org/ 10.2514/3.559. Hartnett JP, Irvine TF, Cho YI, Greene GA, Taniguchi H, Yang W-J, et al. In: Radiative heat transfer by the Monte Carlo method.Aca- Q9 demic Press; 1995. Ferziger, JH, Peric, M. Computational methods for fluid dynamics; 2002. 〈http://dx.doi.org/10.1016/S0898-1221(03)90046-0〉. Patankar S. In: Numerical heat transfer and fluid flow.CRC Press; 1–197.

33

Please cite this article as: Sun Y, et al. Evaluation of three different radiative transfer equation solvers for combined conduction and radiation heat transfer. J Quant Spectrosc Radiat Transfer (2016), http://dx.doi.org/10.1016/j. jqsrt.2016.07.024i

35 37 39 41 43 45 47 49 51 53 55 57 59 61 63 65 67