Numerical simulation on antenna temperature field of complex structure satellite in solar simulator

Acta Astronautica 65 (2009) 1098 – 1106 www.elsevier.com/locate/actaastro

Numerical simulation on antenna temperature field of complex structure satellite in solar simulator Yang Liua, b,∗ , Guo-hui Lic , Li-xiang Jiangd a Marine Engineering College, Dalian Maritime University, Dalian 116023, China b Department of Engineering Mechanics, Tsinghua University, Beijing 10084, China c School of Electronic and Information Engineering, Dalian Jiaotong University, Dalian 116026, China d Beijing Institute of Spacecraft Engineering Environment, Beijing 10094, China

Received 22 October 2007; accepted 5 March 2009 Available online 17 April 2009

Abstract A thermal network model is developed for studying the temperature variation of complicated structure satellite surfaces. The solar incident areas, the infrared and solar radiation transfer coefficients among surfaces are numerically simulated by means of Monte Carlo ray tracing (MCRT) method in model. The non-uniformity and the instability of solar radiation, which plays the important roles in simulating outer-space heat flux designation parameters by solar simulator, are studied for analysis variation of antenna temperature fields in detail. Results showed non-uniformity irradiation effects are greater than those of instability for this kind of geometry sheltering structure. © 2009 Elsevier Ltd. All rights reserved. Keywords: Thermal network; Monte Carlo ray tracing; Temperature; Numerical simulation; Solar simulator

1. Introduction Complicated structure satellite is a kind of spacecraft which is composed of a series of surfaces with different physical properties and sheltering each other. The traditional simulation method of space outer heat flux makes use of infrared lamp as heater. Infrared irradiation has an unparallel and non-uniformity character and leads to difficult to carry out thermal balance test for those of sheltering surfaces due to estimating electric voltage and current un-precisely [1,2]. For high attitude spacecraft with complicated sheltering surfaces, solar simulator is superior to infrared simulator for having a ∗ Corresponding author at: Marine Engineering College, Dalian Maritime University, Dalian 116023, China. Tel.: +86 13717916118. E-mail address: [email protected] (Y. Liu).

0094-5765/$ - see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.actaastro.2009.03.017

better parallelism, uniformity and spectrum irradiation characters from xenon lamp heater [3]. But, solar simulator would be very costly and time-taking for keeping heat sink and vacuum. Numerical simulation technology is the effective way to predict test results by proper physical and mathematical model and saving testing costs and time to a great extent. For solving surface temperature field of thermal radiation, the thermal network method and the finite element method are used in general [4–10]. But, diffused-grey character assumption of surface and view-factor solving method among surfaces confined their application and accuracy. For multifunction complicated spacecraft, their surface character should be expressed as the sum of a specula reflection and a diffuse reflection component [11]. Furthermore, view-factor considered surfaces as ideal black body and absolute geometry relationship,

Y. Liu et al. / Acta Astronautica 65 (2009) 1098 – 1106

it can lead to ignore the physical behavior of surfaces, i.e. reflecting behaviors. So, these defects are limited their applications in the field of engineering. Up to now, thermal network model, which is applied to solving surface temperature field in solar simulator, has not been reported. In order to overcome shortcomings as above, the Monte Carol ray tracing (MCRT) method is used for solving radiation transfer coefficient (RTC). It considered not only space geometry position but also all kinds of reality physical and reflecting property. It is superior to view-factor method for RTC. Meanwhile, non-uniformity and instability uniformity character, which act as the most important designation parameters for solar simulator, is studied on affecting variation of surface temperature field.

1099

6

9

10

1

7

3

11

8 Fig. 2. Calculation mapping.

2. Calculation objects and conditions 2.1. Calculation objects The calculation object is showed in Fig. 1. The device box surfaces are treated as studying objects emphatically, which are located in the rear position of antenna surface. Their temperature distribution and the variation play a very important role to keep instruments inside the device box in a good condition. The convex and concave surfaces take a sheltering function to all surfaces. For global system as shown in Fig. 2, plane surfaces consist of the 1–8th, 11th and 12th surfaces, convex and concave cone surfaces consist of the 9th and 10th surfaces.

instability are equal to ± 4% and ± 1%, respectively. The rotational axis of the spacecraft with respect to the sun irradiation is set as the contact point between 1st and rotates one cycle by itself. Concave cone surface is radiated directly at an angle of 0◦ . The ratio of absorption of polyimide membrane is equal to 0.35 and the ratio of emission is 0.65 in the surface (1, 2, 4–9, 11, 12). The optica1 solar reflector (OSR) membrane in 3rd surface is set as 0.13 and 0.79 and the white paint is 0.20 and 0.87 in 10th surface for them. 3. Mathematical model

2.2. Calculation conditions 3.1. Thermal network model The designation parameters of solar simulator are used, where temperature of heat sink is equal to 100 K, vacuum degree is equal to 10−5 Pa, solar constant is equal to 1353.0 W/m2 , the non-uniformity and the

Surface is divided into micro-surfaces by grid line, and it is represented by nodes i and j following topology sequences. Each node represents the thermal performance and optical characters, as shown in Fig. 2. Neglecting inner heat source, thermal network equation is given in Eq. (1) for node i in enclosure. ND  j=1

Fig. 1. Antenna and prototype of a satellite.

D ji (t j − ti ) +

NR 

R ji (t 4j − ti4 ) = Q i = 0

(1)

j=1

where Dji is the conduct network coefficient, Rji is the radiation network coefficient between node i and node j, NR is the number of radiation network nodes, ND is the number of conduct network nodes, Qi is the total thermal flux received by node i. Heat conduction is classified into general heat conduction and contact heat conduction. Here, Dji is

1100

Y. Liu et al. / Acta Astronautica 65 (2009) 1098 – 1106

defined by Eq. (2) for ignoring general heat conduction. D ji = K ji

A ji L ji

energy bundles Ni is emitted from i surface and is absorbed by j surface Nij , RTC is given in

(2) RT Ci j =

where Kji is the heat conduction coefficient, Aji is the cross-sectional area and Lji is the distance between nodes i and j. The nodes, which would be take into action with (i, j), are (i−1, j), (i+1, j), (i−1, j−1), (i+1, j−1), (i, j−1), (i−1, j+1), (i, j+1) and (i+1, j+1). Rji is given in Eq. (3) for opaque surfaces. in f ra

R ji = 0 · RT C j,i

· j · Aj

(3)

where 0 is Stenfen–Boltzman’s constants (0 = 5.67× 10−8 w/m2 K4 ), RTCij infra is the infrared radiation transfer coefficient, j is the emissivity of j surface, Aj is the micro-area of j surface. Qi is the radiation energy of solar incidence for node i, which includes the solar direct incidence and reflection radiation energy. It is defined as Q i = i · Aisun · Sc + × Asun j · Sc ·

N Rsun

(1 −  j )

j=1 RT C sun ji

Ni j Ni

(5)

When bundles from i surface is making collision with j surface, it may be absorbed or reflected. As far as reflection is concerned, how to judge it, will take effects on calculating accuracy greatly. The novel judgment condition is used. For infrared irradiation process, if generating random number is less than , bundle can be absorbed, the tracing course is stopped. Otherwise, it will be reflected. For solar irradiation process, if generating random number is less than , bundle can be absorbed and tracing is ceased. It can be seen that ratio of absorption and emission acts as the different terms for infrared and solar irradiation in solar simulator. Furis processed by this judgment solving thermore, Asun j with MCRT method. This way has not been reported in spacecraft temperature field literatures. 3.2. Convergence and error of MCRT

(4)

where the first item is considered as the solar direct incident radiation energy for node i, the second term is reflecting radiation energy by other surfaces to node i. i is absorptivity of node i, Sc is the solar constant (Sc = 1353 w/m2 ), NRsun is the number of reflecting nodes, RT Cisun is the solar radiation transfer coeffij is the solar irradiation areas. Exchanged cient, Asun j radiation energy is classified into two categories: one exists among surfaces and the other lies in solar radiation directly. As for infrared RTC, it can be solved by integral method of geometry view-factor or energy balance method of net radiation, i.e. Gebhart’s factor. But, these solving methods are limited to the hypothesis of diffuse and grey condition. For this kind of complicated spacecraft in solar simulator, which exists in the mirror reflection surface (3rd) and solar irradiation, they can not be used any more. According to Kirhhoff laws, ratio of absorption is equal to emission and not with respective to wavelength in infrared irradiation ( = ). But for solar irradiation, ratio of absorption is not equal to emission and has a close relationship with wavelength. So, as mentioned above, in order to comprehend consider optical characters and mirror reflection characters, MCRT is utilized to process these difficult problems [12–14]. MCRT is based on the mathematical statistics and belongs to the stochastic method. Lots of discrete

The absorption or reflection effects will take place between i surface emitting rays and hitting with j surface and has a repeatability and independent behavior [15]. So, the Bernoulli probability distribution follows, as given below: P(k) = Cnk p k (1 − p)n−k , 0  k n

(6)

where Cnk = N !/k!(N − k)! The interval of RTC confidence degree is expressed in Eq. (7) and its estimated value is p − Wc s  RT Ci j  p + Wc s  RT Ciej − Wc

(7)

RT Ciej (1 − RT Ciej )

 RT Cit j  RT Ciej

N + Wc



RT Ciej (1 − RT Ciej ) N

(8)

where p is the estimated value of general sample, s is the mean square deviation of p, w is the quantile fractile. Error expression of RTC is given in Eq. (9). It can be deduced in Eq. (10) for enclosure system.  RT Ciej (1−RT Ciej ) Err Di j ≡|RT Cit j −RT Ciej |  Wc N (9)

Y. Liu et al. / Acta Astronautica 65 (2009) 1098 – 1106

1101

Table 1 Surface describing equations. Surface number

Describing equations

Origin point of surface coordination system

Geometry field of each surface in local coordination system X axis

z=0 z = −80 i = −100 x = +100 y = −100 y = −100 y = −100 z = +30 x2 +y2 = −6.25z2 x2 +y2 = 6.25z2 y = −100 z = +30

1 2 3 4 5 6 7 8 9 10 11 12



 Err Di j  = Err D  Wc

(0,0,0) (0,0, −80) (−100, 0, −40) (100, 0, −40) (0, −100, −40) (0, 100, −40) (0, −150, −95) (0, −225, +30) (0,0,0) (0,0,0) (0, −150, −95) (0, −225, +30)

Diej (1 − Diej )

(−100, (−100, (0,0) (0,0) (−100, (−100, (−250, (−250, (0,0) (0,0) (−250, (−250,



N

Err D  Wc

RT Ci j  =

N n n n 1  1  1 RT C = 1= i j 2 2 n n n i=1 j=1

  1 Err D    Wc

n

1− N

(11)

(12)



=Wc

n−1 N n2

+250) +250)

(−100, +100) (−100, +100) (−100, +100) (−100, +100) (0,0) (0,0) (0,0) (−75, +75) (−200, +200) (−200, +200) (0,0) (−75,+75)

(0, 0) (0, 0) (−40, +40) (−40, +40) (−40, +40) (−40, +40) (−125, +125) (0, 0) (0, +80) (0, +80) (−125, +125) (0, 0)

F(x, y, z) = c1 x 2 + c2 y 2 + c3 z 2 + c4 x y + c5 x z + c6 yz + c7 x + c8 y + c9 z + c10

(14)

where ci (i = 1,2, . . . ,10) is the coefficient. Global and local surface coordination system is employed. The origin point of system coordination lies in the cone convex peak point, and the positive direction is defined by outer-normal direction of surface. The mathematical equations of surfaces in Cartesian coordinates are given in Table 1. The transition from system coordination to local coordination is written as Eq. (15) and from local to system coordination as Eq. (16): ⎡x ⎤

⎡ cos  cos  cos  ⎤ ⎡ x ∗ ⎤ ⎡ x ⎤ 1 0 1 1 p ⎥⎢ ∗⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎣ y p ⎦ = ⎣ cos 2 cos 2 cos 2 ⎦ ⎣ y p ⎦ + ⎣ y0 ⎦ (15) p

i=1

1 n

+100) +100) +250) +250)

Z axis

may be obtained by different equation coefficients. (10)

where RT Cit j is the true value of unknown RTC, RT Ciej is evaluation value of RTC, based on the number of N bundles for sample space. By terms of Jensen inequality, Eq. (10) can be derived as 

  RT C e  1 − RT C e   ij ij

+100) +100)

Y axis

(13)

where n is the number of surfaces. It can be seen from Eq. (13), RTC may be considered as a function of N and n. It provides the minimum number of tracing bundle rays of surface on the conditions that mean error and confidence degree level had been determined. Given that confidence degree level is equal to 0.95 and n is equal to 12, the minimum tracing rays are valued at 300 million rays and error is less than 2.5%.

z ∗p z0 ⎡ x ∗ ⎤ ⎡ cos  cos  cos  ⎤ ⎡ x − x ⎤ 1 2 3 p 0 p ⎢ ∗⎥ ⎢ ⎥⎢ ⎥ (16) ⎣ y p ⎦ = ⎣ cos 1 cos 2 cos 3 ⎦ ⎣ y p − y0 ⎦ zp

cos 3 cos 3 cos 3

z ∗p

cos 1 cos 2 cos 3

where (x0 , y0 , z0 ) is the origin point of local system coordination and (cos i , cos i , cos i , i = 1, 2, 3) is the direction cosine. The solar incidence direction vector is given below: m  x

3.3. Surface describing equation Surface space position can be expressed as the binomial equation in Eq. (14). Different surface positions

z p − z0

my mz

 sin  cos   = sin  sin  cos 

(17)

where  is the zenith angle and  is the round circular angle.

1102

Y. Liu et al. / Acta Astronautica 65 (2009) 1098 – 1106

4. Results and discussion

Fig. 3. It can be seen that temperature distribution took on the quite non-uniformity tendency, and the maximum and minimum values are varied greatly at different angles. This may be explained by the sheltering effects of surfaces, infrared and solar absorption or reflection characteristics. As mentioned above, infrared and solar RTC will play the important factors and affect directly the calculation results.

4.1. Distribution of temperature field The temperature field distribution of antenna and other satellite surfaces is calculated at the 0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300 and 330 irradiation angular, with respect to ideal instability and uniform solar irradiation in solar simulator, as shown in

t (K)

t (K)

z

-0.1 -0.2 -0.3

-0.2

-0.1 y 0

0.1

0.2

0.2

0.1

-0.1

0

-0.2

x

0 -0.1 z

0

290 280 270 260 250 240 230 220 210 200 190 180 170 160 150 140 130 120 110

-0.2 -0.3

-0.2

-0.2

-0.1

-0.1 0 y

0

x

0.1

0.1 0.2

280 270 260 250 240 230 220 210 200 190 180 170 160 150 140 130 120 110

0.2 t (K)

t (K)

-0.1

0 -0.1 -0.2

-0.2

-0.3

-0.3

-0.2

-0.2 -0.1 0 y

0 0.1 0.2

-0.2

-0.2

-0.1

-0.1 0 y

x

0.1

340 320 300 280 260 240 220 200 180 160 140 120

z

z

0

320 300 280 260 240 220 200 180 160 140 120

0.2

-0.1 0 0.1

0.1 0.2

0.2

x

t (K) t (K)

0.2

0 -0.2

0.2

x

0

0 -0.2

-0.2

-0.2

0.2

y

0 x

Fig. 3. Calculation results of temperature field.

-0.2

y

-0.2

0.2

0 z

z

0

320 300 280 260 240 220 200 180 160 140 120

320 300 280 260 240 220 200 180 160 140 120

Y. Liu et al. / Acta Astronautica 65 (2009) 1098 – 1106

0.2

0 z -0.2

t (K)

t (K)

320 300 280 260 240 220 200 180 160 140 120

320 300 280 260 240 220 200 180 160 140 120

0.2

y

y

0

1103

0

-0.2

-0.2

0 -0.1

0.2 0 x

0.2

0

-0.2

-0.2

-0.2

x

z

t (K)

t (K)

0 y

z

-0.1

-0.2

-0.2 0.2

0 x

0

0

y

0

0.2

-0.1 -0.2

-0.2

-0.2

0.2 0.1

-0.1 -0.2

0

x

t (K)

t (K) 300 280 260 240 220 200 180 160 140 120

0.2

0.2 0

x

x

0

0.2 -0.2

-0.2 0

y

-0.1

0 -0.2

-0.2

270 260 250 240 230 220 210 200 190 180 170 160 150 140 130 120 110

y

0

0.2

250 240 230 220 210 200 190 180 170 160 150 140 130 120 110

z

0.2

280 270 260 250 240 230 220 210 200 190 180 170 160 150 140 130 120 110

0 -0.2

z

-0.1 -0.2

z

Fig. 3. (Continued.)

4.2. Instability effects Irradiation instability value is designed at ± 1% deviation for irradiation intensity, that is E = E solar ± 1%E solar . It can be seen in Fig. 4, the worst instability

will increase the deviation of maximum average temperature, especially in 3 and 4 surfaces. Maximum deviation value is larger than normal condition about 1.6 K. A better instability can decrease the temperature deviation for all surfaces.

Y. Liu et al. / Acta Astronautica 65 (2009) 1098 – 1106

2.6 2.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0

Surf_2 normal worst

3.5

Surf_1 normal worst

3.0 2.5 T (K)

T (K)

1104

2.0 1.5 1.0 0.5 0.0

1

2

3

4

5

6

7

8

9

1

10 11 12 13

2

3

4

6

7

8

9

10 11 12 13

3.5

Surf_3 normal worst

2.5

5

Solar incident angle (/)

Solar incident angle (/)

Surf_4 normal worst

3.0

2.0 1.5

T (K)

T (K)

2.5 2.0 1.5

1.0

1.0

0.5 0.5

0.0

0.0

1

2

3

4

5

6

7

8

9

10 11 12 13

1

2

3

4

5

6

7

8

9

10 11 12 13

Solar incident angle (/)

Solar incident angle (/)

Fig. 4. Maximum temperature deviation of the normal and the worst irradiation instability.

4.3. Non-uniformity effects Irradiation non-uniformity value is designed at ± 4% deviation for irradiation intensity, that is U =U0 ±4%U0 . The non-uniformity parameter is defined as U = ±

E max − E min 100% E max + E min

(18)

where U is non-uniformity of irradiation, Emax and Emin are maximum and minimum values of the solar irradiation intensity, respectively. It is difficult to make a relationship between irradiation intensity and tracing rays. The new parameter solar incidence area Asolar is established for the first time for studying non-uniformity irradiation according to MCRT methods. Solar incident area is the sum of the solar direct incident area from solar incidence directly and the solar un-direct area from the other reflecting surface. In MCRT method, the direction and space position of ray

are determined stochastically. The course of calculating the direct incident area is the emitted ray from surface to outer environment along the inverse-direction of solar incidence. If rays can be sheltered by other surfaces, it is considered as an un-shined point. Otherwise, it is treated as a shined point. The number of shined points is numerically stated and is divided by the total number of simulation rays, and then multiplied by the projection area. So, the solar direct incident area would be obtained as follows: Asolar = f (x0 , y0 , z 0 , 0 , 0 , x, y, z)

(19)

Uniformity can alternate the emitting direction along solar incidence direction and lead to alter irradiation intensity. It can be written as Asolar = f (x0 , y0 , z 0 , 0 , 0 , x, y, z) − f (x0 , y0 , z 0 , 1 , 1 , x, y, z)

(20)

Y. Liu et al. / Acta Astronautica 65 (2009) 1098 – 1106 7

1105

8 Sruf_1 normal worst

6

Surf_2 normal worst

7 6

5 T (K)

T (K)

5 4 3

4 3

2

2

1

1 0

0 1

7

2

3

4

5 6 7 8 9 10 11 12 13 Solar incident angle (/)

1

Surf_3 normal worst

6

3

4

5 6 7 8 9 10 11 12 13 Solar incident angle (/)

8 Surf_4 normal worst

7 6

5 4

T (K)

T (K)

2

5 4

3 3 2

2

1

1 0

0 1

2

3

4

5

6

7

8

9

10 11 12 13

Solar incident angle (/)

1

2

3

4

5

6

7

8

9

10 11 12 13

Solar incident angle (/)

Fig. 5. Maximum temperature deviation of the normal and the worst irradiation non-uniformity.

Surface irradiation intensity explained is equal to the product of the solar constant Sc and the surface Asolar , as in Eq. (19). E = Sc Asolar

(21)

Eq. (20) builds relationship uniformity with maximum and minimum number of light bundles within radiation zones. So, non-uniformity is given in Eq. (21): U = ± = ±

E max − E min 100% E max + E min Asolar,max − Asolar ,min Asolar,max + Asolar ,min

100%

(22)

During the course of calculation, non-uniformity is replaced by the variation in Asolar , which can be obtained by MCRT method.

It can be seen in Fig. 5, the worst non-uniformity will increase the deviation of maximum average temperature, especially in 3 and 4 surfaces. Maximum deviation value is larger than normal condition about 2.4 K. A better instability can decrease the temperature deviation for all surfaces.

4.4. Comparison between instability and non-uniformity Surfaces 2 and 3 are considered as typical surfaces to study maximum deviation of temperature variation, effected by the worst instability and non-uniformity. In Fig. 6, it can be seen that variation amplitude, caused by the radiation non-uniformity, is larger than that of instability effects. Maximum value is reached up to 6.5 K.

1106

Y. Liu et al. / Acta Astronautica 65 (2009) 1098 – 1106 9

7 Surf_2

8

instability

instability

non-uniform

7

non-uniform

5

6 5

T (K)

T (K)

Surf_3

6

4 3

4 3 2

2 1

1 0

0 1

2

3

4

5

6

7

8

Solar incidence angle (/)

9

1

2

3

4

5

6

7

8

9

Solar incidence angle (/)

Fig. 6. Deviation comparison between non-uniformity and instability of irradiation.

5. Conclusions (1) Combining thermal network model with MCRT method it can be used to predict the temperature distribution of sheltering surfaces and provide reference data for thermal test in solar simulator. (2) The novel definition of infrared and solar reflection judge and solar incidence by MCRT is contributed comprehensively for considering the effects of optical characteristics and complicated geometry factors for temperature in success. (3) The non-uniformity of irradiation has a greater effect on temperature variation than that of instability. Acknowledgments This study was sponsored by the pre-research projects of China Academy of Space Technology (CAST) under the Grant no. 900311. Furthermore, we sincerely appreciate the direction of Professor Ben-cheng Huang of CAST. References [1] B.C. Huang, Y.L. Ma, Space Environment Test Technology of Spacecraft, National Defenses Industry Press, Beijing, China, 2002, pp. 137–154. [2] G.R. Min, S. Guo, Thermal Control of Spacecraft, second ed. Science Press, Beijing, China, 1998, pp. 456–473. [3] E.K. Latvala, Solar Simulation Methods, AIAA Journal (1970) 38282–38298. [4] G.R. Min, J.G. Hu, The research of space thermal simulation testing methods for future large spacecraft, Acta Astronutica 9 (10) (1981) 321–329.

[5] A.K. Oppenheim, Radiation analysis by the network method, ASME Journal of Heat Transfer. 78 (4) (1956) 725–735. [6] J.H. Weng, Z.F. Pan, G.R. Min, Thermal network analysis of an experimental spacecraft, China Space Science and Technology 10 (5) (1998) 56–60. [7] J.H. Weng, Z.F. Pan, A correction method for spacecraft thermal network and its coefficients, China Space Science and Technology 8 (4) (1995) 10–15. [8] J.H. Weng, Z.F. Pan, A method of calculating experimental heat flux on a arbitrary shaped convex surfaces of spacecraft, China Space Science and Technology 14 (2) (1994) 11–18. [9] H.P. Lee, Application of finite-element method in the computation of temperature with emphasis on radiative exchange, AIAA Journal (1972) 274–281. [10] H.P. Lee, C.E. Jackson, Finite-element solution for a combined radiative conductive analysis with mixed diffuse-specular surface characteristics, AIAA Journal (1975) 682–689. [11] E.M. Sparrow, E.R.G. Eckcrt, V.K.A. Jonsson, An enclosure theory for radiation exchange between specularly and diffusely reflecting surfaces, Journal of Heat Transfer. 84 (11) (1962) 294–300. [12] D. Polder, M.V. Hove, Theory of radiative heat transfer between closely space bodies, Physical Review B 4 (10) (1971) 303–314. [13] H.P. Tan, G.M. Cui, L.M. Ruan, Combining domain decomposition method with Monte Carlo method to determine the temperature profile and infrared radiance of tank, Journal of Engineering Thermophysics 19 (3) (1998) 340–344. [14] H.P. Tan, G.M. Cui, L.M. Ruan, The Monte Carlo method and parallel algorithm in infrared mapping model of terrain object, Journal of Infrared and Millimeter Waves 17 (6) (1998) 417–423. [15] J.D. Malt, P.J. Burns, Performance, accuracy, and convergence in a three dimensional Monte Carlo irradiative heat transfer simulation, Numerical Heat Transfer, Part B 19 (2) (1991) 191–209.