Reverse Monte Carlo Method for Transient Radiative Transfer in Participating Media

The Monte Carlo (MC) method has been widely used to solve radiative transfer problems due to its flexibility and simplicity in simulating the energy transport process in arbitrary geometries with complex boundary conditions. However, the major drawback of the conventional (or forward) Monte Carlo method is the long computational time for converged solution. Reverse or backward Monte Carlo (RMC) is considered as an alternative approach when solutions are only needed at certain locations and time. The reverse algorithm is similar to the conventional method, except that the energy bundle (photons ensemble) is tracked in a timereversal manner. Its migration is recorded from the detector into the participating medium, rather than from the source to the detector as in the conventional MC. There is no need to keep track of the bundles that do not reach a particular detector. Thus, RMC method takes up much less computation time than the conventional MC method. On the other hand, RMC will generate less information about the transport process as only the information at the specified locations, e.g., detectors, is obtained. In the situation where detailed information of radiative transport across the media is needed the RMC may not be appropriate. RMC algorithm is most suitable for diagnostic applications where inverse analysis is required, e.g., optical imaging and remote sensing. 1 In this study, the development of a reverse Monte Carlo method for transient radiative transfer is presented. The results of non-emitting, absorbing, and anisotropically scattering media subjected to an ultra short light pulse irradiation are compared with the forward Monte Carlo and discrete ordinates methods results. INTRODUCTION The recent research on the propagation of ultra-short light pulse inside the absorbing and scattering media has lead to some interesting applications in the area of material properties diagnostic, optical imaging, remote sensing, etc. The time scales of such processes are usually in the order of 10 to 10 seconds. In the case of remote sensing using short light pulse, the pulse width is in the order of 10 seconds. The corresponding spatial and temporal variations of radiation intensity in these processes are comparable. Therefore, the consideration of the transient term in the radiation transport equation is necessary. The simulation of transient radiation process is more complex than that in the steady state due to the hyperbolic wave equation coupled with the in-scattering integral term. Several numerical strategies have been developed, which include discrete ordinate method (Sakami, et al. 2000; 2002), finite volume method (Chai, 2003; Lu, et al. 2003), integral equation models (Tan and Copyright © 2003 by ASME Hsu, 2001; Wu, 1999; Wu and Wu, 2000), and Monte Carlo method (Brewster and Yamada, 1995; Hsu, 2001a). The Monte Carlo method is a numerical technique of solving various sciences and engineering problems by the simulation of random variables. Monte Carlo is one of the most versatile and widely used numerical methods, very suitable for solving multidimensional problems, especially when deterministic solutions are difficulty to obtain. Monte Carlo simulations have in the past, and continue to, consume a significant fraction of high performance computing time. With the advancement of the low cost Beowulf cluster (Sterling, et al. 1995) and inherently high parallel efficiency of Monte Carlo method (Siegel and Howell, 2002; Sawetprawichkul, et al. 2002), its usage will only grow over time. The theoretical foundation of the method had been known long before digital computers were available. Of course, the use of the Monte Carlo method became practical only with the advent of computers and high-quality pseudorandom number generators. In the early days of the digital computer, significant amount of computational work was carried out with Monte Carlo simulations of the neutron transport. Radiation heat transfer is one of the disciplines that take great advantage of the Monte Carlo method (Howell, 1968). However, the Monte Carlo method can be computationally very time consuming due to its error bound being limited by the number of sampling in the form of O(N), in which N is the sampling number. That is why much of the effort in Monte Carlo development has been in construction of variance reduction methods that speed up the computation. The other approach, which has shown promising results, is the utilization of deterministic number sequences instead of random numbers for integration (O’Brien, 1992; Hsu, 2001)– this belongs to the socalled quasi-Monte Carlo (QMC) method (Holton, 1960; Niederreiter, 1978). In theory, the error bound of QMC is proportional to O(N(logN)), where D is the dimension of integration or the number of deterministic number sequence. In practice, depending on the quality of the number sequences and the nature of the problem, the error bound is typically between O(N) to O(N) (Hsu, 2001b; Bratley, et al. 1992). In the case of detecting radiation signals at some selected positions and/or given time intervals, the reverse Monte Carlo (RMC) method becomes very advantageous. Since computational results in the whole spatial and temporal domains are not always necessary, then either MC or QMC 2 becomes unnecessarily inefficient. The RMC is based on the reciprocity principle in radiative transfer theory (Case, 1957). The reverse algorithm is similar to the conventional method, except that the energy bundle (photons) is tracked in a time-reversal manner. Its migration is recorded from the detector into the participating medium, rather than from the source to the detector as in the conventional MC. There is no need to keep track of the bundles that do not reach a particular detector. Thus, RMC method takes up much less computation time than the conventional MC method. On the other hand, one should always note that RMC will generate far less information than the corresponding MC method. Collins, et al. (1972) reported the earliest work on RMC relevant to radiation transfer calculations. Adam and Kattawar (1978) adopted the same concept of Collins, et al. in their study of spherical shell atmospheric radiation. They also provided justifications of RMC, which will be expanded in this study. In a series of rocket plume base heating calculations by Nelson (1992), RMC was also found to be a practical tool to include various effects, e.g., spectral, scattering. To provide a more rigorous theoretical foundation to RMC, Walters and Buckius (1992) utilized the reciprocity relations developed by Case (1957). Recently, Modest (2003) used their results and applied to singular light source problems. The only RMC treatments for transient radiation processes reported so far, to the authors' knowledge, are by Wu and Wu (2000) and Andersson-Engels, et al. (2000). However, the details of the RMC algorithm were not given in either paper and the RMC results of the latter work were inconsistent. This study will present a detailed and validated RMC method to simulate transient radiation process, particularly for light pulse propagation within scattering media, with validated solutions. The utilization of a Beowulf cluster for simulation demonstrates the RMC method has the potential for real time inverse analysis. Although only onedimensional geometry is considered in this study, the algorithm can be readily extended into multidimensional geometry, which is to be treated in a follow-up study. TRANSIENT RADIATION TRANSPORT Consider a one-dimensional slab containing an absorbing and scattering medium, the radiative transfer equation (RTE) in a given direction ŝ is