Open Access
Add to BibSonomyAbstract
Based on the kinetic theory of photons, a new lattice Boltzmann method for the simulation of 3D radiation transport is presented. The method was successfully validated with Monte Carlo simulations of radiation transport in optical thick absorbing and non-absorbing turbid media containing either isotropic or anisotropic scatterers. Moreover, for the approximation of Mie-scattering, a new iterative algebraic approach for the discretization of the scattering phase function was developed, ensuring full conservation of energy and asymmetry after discretization. It was found that the main error sources of the method are caused by linearization and ray effects and suggestions for further improvement of the method are made.
© 2016 Optical Society of America
1. Introduction
Radiation transfer is an important issue in several fields of engineering. Typical examples are radiative heat transfer or photon transport in technical optics. In many of those applications, fluids are present and affect radiation transfer by their physical properties. In particular in turbid particle-containing media, spatial concentration gradients of particles do influence the local absorption, emission and scattering of radiation. Conversely, radiation may lead to locally varying fluid temperature, which affects temperature-dependent fluid properties, convective and diffusive transport phenomena or reaction rates within the fluid. Reaction rates might be also linked directly to irradiance, for instance in the case of photocatalytic reactions [1] or the cultivation of photosynthetic microorganisms [2,3]. Accordingly, for accurate prediction, design and optimization of such systems, the coupling between radiation transport, flow and reaction rates has to be taken into account. Owing to the complexity and non-linearity of the arising mathematical problem, analytic solutions can be found just for very few simplified cases, necessitating the application of numerical methods.
Most commonly used for numerical calculations of radiative transfer are Monte Carlo methods (MCM), Discrete Ordinate methods (DOM) and Finite Volume methods (FVM). In case of MCM, a sample of individual photons is traced from the point of emission to the point of absorption. Since no angular discretization is required, MCM are considered to be very accurate, even if complex problems of radiation transfer must be solved. However, because of their statistical nature, MCM may lead to noisy outcomes [6]. Moreover, the size of the traced photon sample must ensure statistical certainty which leads to slow convergence, particularly if optical thick media are simulated [7]. In contrast, both, DOM and FVM, are applied frequently as a compromise between accuracy and speed of convergence [8] and enable also transient computations [9–11]. These methods are in essence similar and based on a directional discretization of the Radiative Transfer Equation (RTE) and integration by means of specific quadrature rules. Besides well-known error sources such as the ray effect or numerical smearing (also referred to as false scattering) [12, 13], angular false scattering as a third type of error was considered [14, 15]. It is caused by inaccurate discretization of the scattering phase function and becomes significant for strong anisotropic scattering.
In the recent past, as a new approach, lattice Boltzmann methods (LBM) have been developed for the simulation of radiation transport. The idea behind seems to be straightforward since the RTE as the governing equation of radiation transport can be interpreted as a Boltzmann-type transport equation for photons [4, 5]. By combining both, accuracy and simplicity, lattice Boltzmann methods were established in the past as powerful tools for the simulation of several physical problems. Originating from the kinetic theory of gases [16], these methods were primarily developed to study fluid dynamical problems. Many research was done to extend the methodological framework to complex flow problems, such as turbulence, thermohydrodyamics, multiphase, reactive or non-newtonian flows and many others [17–26]. However, the existence of LBM for applications beyond fluid dynamics, such as quantum mechanics [27] or phonon transport [28], shows the principal applicability of the lattice Boltzmann framework to a wide range of systems which can be characterized by means of statistical physics. With regard to coupled flow and radiation transport, new lattice Boltzmann methods for radiation transfer potentially allow the fast and efficient simulation of these systems by using a consistent set of numerical methods on structured grids, even in complex geometries such as porous media. In addition, it is well known that LBM show great opportunities for parallel computing, which is of special importance for calculations of radiative transfer since the multidimensionality of the problem causes significant computational costs [7].
However, the development of lattice Boltzmann methods for the simulation of radiation transfer is still a new field of research. Pioneer work in this field was done by Asinari et al. [29], Ma et al. [30] and DiRienzo et al. [31] and applied to benchmark cases in one and two dimensions. While those approaches dealt with isotropic scattering media only, Bindra and Patil [32] extended the approach of [30] to anisotropic scattering media by approximating the in-scattering integral by summation over a Gaussian quadrature. Moreover, it was shown that interparticle collision of photons are negligible and modeling of absorption, emission and scattering by means of source terms is sufficient. Further work was less fundamental than rather focused on applications. In particular, the simulation of collimated radiation by means of lattice Boltzmann methods was objected. For that purpose, based on the work of Asinari et al. [29], a flux splitting approach was designed by Mishra and Vernekar [33] and applied for the transient simulation of collimated radiation in one dimension [34]. Here, the time stepping scheme was explicit. Zhang et al. [35] treated transient transfer of radiation in one dimension by using an implicit discretization for the time derivative in the RTE. The method was applied for transient simulations of radiative transfer in index graded media [36] as well as for the transient simulation of scattered collimated radiation in combination with a flux splitting scheme [37, 38]. McCulloch and Bindra [39] extended the approach given by [32] to higher order angular discretization schemes in two dimensions and applied their method for coupled radiation and convective heat transfer. Recently, further developments concerned the extension of LBM solvers for radiation transfer to three-dimensional, isotropically scattering media [40].
Although progress has been made in recent years, the development of a consistent methodological framework is still outstanding. For instance, directional weights are either constructed geometrically [29, 31] or derived from moment equations [39]. Moreover, most of the works mentioned previously did not consider radiation transfer in media containing anisotropic scatterers. To close this gap, the simulation of radiative transfer in three-dimensional domains with anisotropic scatterers by means of lattice Boltzmann methods is shown in this work for the first time. The method is derived from the kinetic theory of photons and takes fundamental differences to particle kinetics into account. For angular discretization, sets of discrete directions in three dimensions are developed ensuring the fulfillment of angular moment equations. In addition, to take scattering by Mie-scatterers into account, an innovative algebraic approach is applied to discretize the Henyey-Greenstein phase function under full conservation of energy and asymmetry. Finally, the method is applied in numerical examples to study its accuracy in comparison to Monte Carlo methods.
2. Mathematical formulation
2.1. Brief review of kinetic theory of photons
The origin of lattice Boltzmann methods lies in the kinetic theory of gases. Thus, for the calculation of radiative fields by means of LBM it is essential to consider the kinetic theory of photons and take differences to particle kinetics into account. In case of gas particles, the fundamental idea is that particle motion follows Newtonian mechanics and particles can interact via interparticle collisions [4]. On a mesoscopic scale, microkinetics of single particles are captured statistically by the Boltzmann equation, Eq. (1), which governs the kinetics of a whole particle ensemble.
Herein f is the velocity density distribution, which defines the probability of finding a gas particle in a phase space volume ΔxΔξ at time t, whereby x denotes the position vector and ξ the momentum vector. The velocity density distribution changes in time by convection, by influence of external forces F on gas particles of mass m or due to interparticle collisions as defined by the collision operator Γ. In gas kinetic theory, a famous formulation of Γ is given by the Bhatnagar-Gross-Krook (BGK) approximation [16, 19]. It is assumed that the velocity density distribution relaxes due to interparticle collisions towards local equilibria, which are defined by the Maxwell-Boltzmann velocity distribution for gas particles of equal mass. More general, Γ can be interpreted as the difference between losses and gains of gas particles per unit time due to interactions within the phase space volume under consideration. According to this interpretation, the collision operator can be rewritten in the following way.
A kinetic equation similar to Eq. (1) can be derived for photons by replacing f by a quantity ϕ, which is known as the photon density distribution. By definition, ϕ measures the probability of finding an unpolarized photon in a phase space volume ΔxΔp, where p is the momentum vector of photons, whose magnitude is p = hν/c [4]. For photons, an element within the three-dimensional momentum space Δp consists of an element dν in the one-dimensional frequency space and an element in the two-dimensional directional space, represented by a solid angle dΩ, which contains all vectors n = p/p for p ∈ Δp. By using this definitions, a phase space volume for photons becomes to
where c is the speed of light and h the Planck constant. As an important complement, according to Liouville’s theorem, the photon number density in phase space is conserved in absence of collisions [5].Radiative fields are commonly expressed in terms of the specific intensity Iν, which is defined as the amount of energy dE, carried by photons within a frequency band dν into a solid angle dΩ per unit time dt and unit area dA [6]. The specific intensity of photons is related to the photon density distribution by the expression
and hence, since c and h are constants, Iν/ν3 can be used to derive the kinetic equation of photons [4]. Although it is known that gravity and electromagnetic forces act on the momentum of photons, in the following it is assumed that momentum of photons changes due to interaction with matter only. By using the general formulation of the collision operator as given by Eq. (2) and multiplying with ν3, the kinetic equation for photons arises.The indexes G and L again mark gains and losses of photons within the phase space volume under consideration. Similar to Eq. (1), Eq. (5) governs the kinetics of a photon ensemble on a mesoscopic scale. However, with regard to the mechanism of collisions, fundamental differences between photons and particles exist. While particles exchange momentum due to interparticle collisions, changes in momentum of photons occur by interaction with surrounding matter. As a first mechanism, absorption and isotropic emission of photons by matter occurs. If matter is in local thermal equilibrium (LTE), the momentum distribution of photons also moves towards an equilibrium, whose frequency distribution is Plankian [41]. As a second mechanism, elastic scattering of photons occurs in turbid media [42], which affects only the directional component of the photon momentum vector while frequency stays constant [4]. By taking these mechanisms of interaction into account, Eq. (5) becomes to the well known Radiation Transfer Equation (RTE), which gives a description of photon kinetics in turbid media in terms of the specific intensity.
The right-hand side of Eq. (6) capture the extinction of photons due to absorption and scattering by matter, the re-emission of photons with a specific intensity Iν,e and the local in-scattering of photons into the direction of propagation n from directions n′. Here, the scattering phase function Φ represents the angular distribution of scattered photons according to the scattering properties of the surrounding matter. The extinction coefficient βν = κν + σν measures the interaction of photons with matter per unit volume due to absorption and scattering, both individually expressed by the absorption (κν) and scattering (σν) coefficients, respectively. Thus, βν is related to the mean free path (MFP) lf of photons by
Eq. (7) is valid for homogeneous media only, while for non-homogeneous media both, the local MFP lf (x) as well as the length scale of its spatial variability, should be evaluated to avoid underestimation of free path lengths [43]. If a medium is homogeneous, the probability for collision of photons with matter along a path of length s can be shown to be Pcoll(s) = 1 − exp(s/lf) [6].
The local radiative field can be characterized by moment equations of the specific intensity by integration over the momentum space. A general formulation for the radiative moment Mk of order k of the specific intensity reads as
where ψk(c) is a polynomial of c. Only the zeroth (mean intensity J), first (radiative flux) and second (radiation pressure) moments have physical meaning [6]. Nevertheless, higher order moments of the radiation field are computed for some angular discretization schemes in the DOM [44]. Since the discretization in LBM is also based on the satisfaction of several moment equations [45], it is important to take the different definitions of momentum space for particles and photons into account.2.2. Lattice Boltzmann method
The RTE as given by Eq. (6) has to be discretized both, in time and phase space. Moreover, the scattering phase function Φ needs to be discretized to solve in-scattering integral in the discretized momentum space. For this purpose an algebraic technique will be applied to the Henyey-Greenstein phase function.
2.2.1. Time discretization
By using the substantial derivative, Eq. (6) can be reformulated as a pseudo-linear ordinary differential equation.
Herein, Sν denotes the source term
where ων = σν/βν is the scattering albedo. Making usage of the fundamental theorem of calculus, an exact expression for the left-hand side of Eq. (9) can be reached by integration over a small time intervall Δt along a characteristic line c = cn [45]. The time integral of the right-hand side is usually approximated by a rectangular rule so that the time-discrete formulation of Eq. (9) reads as follows.2.2.2. Discretization of phase space
In LBM, space is discretized by a structured lattice, whereby the lattice nodes represent discrete positions between which photons are able to propagate, as shown in Fig. 1. Accordingly, as a characteristic of LBM, discretization of momentum is connected to the discretization of space. An assumption in LBM is that particle or photons stream freely between the nodes of the lattice and collisions occur on the nodes [16]. This means that the grid spacing should be in the order of the MFP, as defined by Eq. (7). This local description of photon kinetics holds in the continuum limit as characterized a low Knudsen number Kn = lf/L, which relates the MFP to a characteristic length of the system L. In other words, local fluctuations by collisions do not affect the mean radiation properties on a macroscopic scale, if many collisions occur.
Fig. 1 2D schematic representation of a discretized phase space in the LBM. Arrows symbolize photon propagation between lattice nodes.
As a simplification, monochromatic radiation is assumed from now on, and thus, the six-dimensional phase space reduces to five dimensions. For convenience the notation of the specific intensity is changed to I and moreover, all indexes referring to frequency-dependency are dropped. To solve the RTE on the lattice, I is discretized into discrete specific intensities Ii, which propagate with constant speed along characteristic lines ci. Thereby the choice of the ci is constrained by the arrangement of the lattice as discussed above. On a lattice node, the macroscopic quantities of the radiative field must be evaluated by the calculation of radiative moments. Thus, also in a discretized momentum space, the exact evaluation of moment equations, Eq. (8), is required, which again necessitates the approximation of the integrals by a quadrature scheme. Assuming isotropy in radiative equilibrium, the integral is of the type ∫aψ(x)dx and can be evaluated by a Gaussian quadrature, so that Eq. (8) becomes to
where wi are the weight coefficients, corresponding to the abscissas of the quadrature as given by a set of ci. It should be noted that Eq. (12) approximates the inner integral in Eq. (8), while the integration in frequency space is not considered. However, Eq. (12) represents the approximation of an integral on the surface of the unit sphere, which is in common to the calculation of radiative moments in the DOM. Since a structural similarity between both methods is known [46], design principles for the quadrature rule can be adopted from the DOM. Accordingly, the quadrature must ensure rotational invariance, symmetry and conservation of energy [44]. A common ansatz function for the polynomial ψk(ni) is where li, mi, ni are the direction cosines of a quadrature abscissa, which must lie on the unit sphere, so that . The moment order k is related to the exponents by k = k1 + k2 + k3. By choosing k1 = k2 = 0 and k3 = k, the requested condition of energy conservation is ensured [44]. From symmetry, it follows that all odd moments become zero. For even moments a general solution was derived for the absolute value of moments of order k [47]. For the specific case of k3 = k, the quadrature must be able to fulfill the following conditions exactly.For a set of given quadrature abscissas as defined by the lattice arrangement, Eqs. (12)–(14) can be reformulated in terms of a linear equation system and solved for the unknown quadrature weights.
In a three-dimensional structured lattice, the solution of Eq. (15) is similar to Gaussian-type quadrature schemes as constructed by Lebedev [48] for the integration on the unit sphere. The derivation of the quadrature stands in contrast to previous approaches [29, 31, 39], wherefore important differences to fluid LBM shall be stressed. First, photons propagate with constant speed so that the quadrature integrates on the unit sphere. This is an important difference to common LBM, where particle velocities differentiate in direction and speed. In consequence, the weighting factors do not incorporate different lengths of photon propagation, which may be present in a cubic lattice. To address these effects in the present model, the incorporation of different propagation lengths is realized by direction-specific integration as discussed below. Second, the different equilibrium functions for particles and photons lead to different structures of the integrals which have to be calculated to evaluate the moment equations. In contrast to photons, for particles the equilibrium is Maxwellian and the integral is approximated by Gauss-Hermite quadratures, which, unlike the Lebedev quadrature, require a “zero-velocity”. In a physical interpretation, this difference can be explained by the fact that photons travel with constant speed and, thus, no resting energy has to be taken into account.
As stated before, the choice of the discrete directions or quadrature abscissas is constrained by the three-dimensional lattice. Three possible arrangements are shown in Fig. 2. In accordance to the usual nomenclature, the sets of discrete directions are termed DnQm, where n is the dimension of the problem and m the number of discrete directions within one set.
Fig. 2 Sets of discrete directions in a three-dimensional lattice. Arrows symbolize the discrete directions of photon propagation between lattice nodes.
By coupling the discretization of time, Eq. (11), and phase space, a fully discretizated representation of Eq. (6) can be reached. The in-scattering integral is replaced by a quadrature rule. The resulting model equation for the LBM framework reads as follows.
Herein, the left-hand and right-hand sides are denoted as streaming and collision terms, respectively. With regard to Fig. 2, it becomes clear that, except from the D3Q6 set, the distance between two nodes in the lattice depends on the direction of propagation. Because the distance Δxi = ciΔt is included in the streaming term, Eq. (16) is valid for the D3Q6 set only. Since the speed of light is constant, streaming on a structured lattice requires direction-dependent time steps Δti. Because the time integration uses non-uniform time-steps, the amount of energy carried by the discrete specific intensities depends not only on scattering, emission and absorption but also on the direction dependent time-step. This numerical effect becomes important for the computation of in-scattering, since the in-scattering term couples different discrete specific intensities. As a correction, the discretized specific intensity is normalized per unit time and corrected to the specific time-step. In consequence, a limitation of the LBM is the computation of steady-state intensity fields only. With this modification, the resulting model equation becomes to the following expression.
Table 1 sums up the sets of discrete directions and resulting quadrature weights for different discretizations of momentum space for the LBM.
Table 1. Quadrature sets for different discretizations of momentum space. The quadrature abscissas result from all possible permutations and multiplication by −1 of the typical vector elements.
2.2.3. Discretization of the scattering phase function
The evaluation of the in-scattering term requires a discretized representation of the phase function Φ, which represents the angular distribution of scattered photons. For media containing Mie-scatterers, the Henyey-Greenstein phase function has often been used to approximate Mie-scattering and is given by
Herein, g is the asymmetry factor, which is equal to the mean cosine of the scattering angle θ. The asymmetry factor can take any value within the interval [−1 ≤ g ≤ 1] to represent anisotropic as well as isotropic (g = 0) scattering. Both, the continuous as well as the discrete representation of Φ must ensure the conservation of energy and asymmetry and fulfill the following conditions.
Violation of these conditions causes inaccurate computation of the scattering kernel since either energy sources arise or the shape of the scattering phase function is not conserved in the discrete phase space. The latter error source is known from the DOM and termed as “false angular scattering” [14, 15]. To reduce the effects of inaccurate discretization, normalization methods can be applied to ensure conservation of energy or asymmetry. Recently, a normalization method was developed which conserves both, energy and asymmetry of the discretized scattering phase function [49,50]. The idea behind this technique is that a set of discrete values Φi,j, which does not fulfill the conditions as given by Eqs. (19a) and (19b), is normalized by
so that the discrete values Φ̃i,j fulfill the conservation of energy or asymmetry. The computation of the normalization factors Ai,j leads to an under-determined linear equation system, which is solved by QR-decomposition and computation of the minimum norm of the linear equation system.In this work, the approach of [49, 50] is adopted to discretize the scattering phase function iteratively. The basic idea behind is that the computed set of discrete Φ̃i,j is more accurate after normalization, if the set of Φi,j was already a good initial solution. However, there is no solution known for Φi,j except from the isotropic case (g = 0), where Φi,j = 1 is true for all i, j ≤ m. Since this solution is exact, it is a good approximation for a slightly different value of g1 = giso + Δg and can be used as an initial solution for the normalization. The normalized solution ensures conservation of energy and asymmetry and is again a good approximation for Φ̃i,j at a new slightly different g2 = g1 + Δg. The procedure is repeated until a desired value of g is reached and a discretized set of Φi,j is computed. Compared to the explicit calculation of Φi,j by integration of Φ(θ) and averaging over discrete solid angles dΩi and dΩj, an advantage of the iterative technique is that the calculation is fully algebraic and no integration bounds have to be defined on the unit sphere. To specify the step size Δg, a discretization error εg can be computed as [51]
where ϕ1 is a m × m matrix containing the values of Φi,j for a desired value of g, reached with a certain Δg and ϕ2 is the solution for the same g, but computed with the next bigger step size. If the step size is changed by a factor r, the order of error reduction p can be computed as The discretization error at different step sizes is shown in Fig. 3 for the three quadrature sets as defined in table 1 and a final value of g = 0.9. It can be seen that the error decreases linearly with step size so that the accuracy of the iterative discretization is of first order. For a step size Δg = 1E–04, the magnitude of the discretization error becomes 1E–06 or lower. At this step size, the required computational time for the discretization was between 2 and 34.5 seconds on a 3.4 GHz machine, depending on the number of discrete directions. Regardless of the chosen step size, the discretized phase functions fulfill the conditions as given by Eqs. (19a) and (19b) exactly within the floating-point relative accuracy.Fig. 3 Left: Discretization error with respect to the step size Δg for a final value of g = 0.9 (see text for details). Right: Discretized (symbols) and continuous (solid line) Henyey-Greenstein phase function for g = 0.9 and Δg = 1E–04.
2.2.4. Nondimensionalization
The model equation, Eq. (17), still represents physical units. Although simulations can be carried out by using physical units, it is convenient to transform the model equations into a non-dimensional representation. With regard to the RTE, an adequate choice of reference quantities would be a characteristic length L0, a characteristic mean specific intensity J0 and a characteristic speed, which is the speed of light c. By using these reference quantities, a set of dimensionless quantities (index D) can be defined.
The dimensionless extinction coefficient is also known as the optical depth τ, which needs to be equal in both, the physical and the non-dimensional representation. Since β is the inverse of the MFP, the optical depth is the inverse of the Knudsen number Kn. As a restriction, LBM are applicable for simulations at low Knudsen numbers so that Kn ≈ 1E–02 [16], which corresponds to optical thick media (τ ≫ 1) [6]. With regard to Eq. (23) it turns out that the dimensionless speed equals unity, since photons propagate with constant speed of light. Furthermore, it holds that lD,0 = 1 and JD,0 = 1. In a discretized representation, the non-dimensional length becomes to
where Δx is the grid spacing and N the number of nodes in the lattice used to discretize the characteristic length. As mentioned above, the grid spacing should be in the order of the MFP. From Eqs. (23) and (24) follows that tD = lD/cD or ΔtD = ΔxD/cD. As discussed before, with the LBM the steady-state radiation field can be computed and, hence, the model becomes time-independent. In the dimensionless notation, Eq. (17) reads as where the dimensionless length in the lattice is . The dimensionless mean specific intensity can be computed as2.3. Algorithm
The implementation of the LBM follows the classical lattice Boltzmann algorithms. After initialization, the following procedure is repeated until a desired degree of accuracy is reached. First, the collision of photons with matter is calculated for all discrete directions according to the right-hand side of Eq. (25). In the second step, the boundary conditions are applied, leading to an intermediate state of the specific intensity.
The third step is the so-called streaming step, where the new specific intensity is calculated. The accuracy is measured in terms of the relative difference of the discretized specific intensity between two computational steps.2.4. Stability analysis
As a criterion for stability, the parameter settings must ensure that the specific intensity remains positive. A limiting case exists, if radiation is transported in one discrete direction and in-scattering results only from forward scattering. From Eq. (25) it follows for an absorbing, scattering and non-emitting turbid media that in that limiting case the condition
must be fulfilled to ensure positivity of specific intensity. Figure 4 shows valid parameter combinations for the D3Q26 discretization, which guarantee stable simulations. The limitation of stability was found to be caused by the direction vectors 7–14 (see table 1) for both, the D3Q14 and D3Q26 discretization. As stated before, the condition of low Knudsen numbers necessitates high optical thickness, so that in particular simulations of isotropically scattering or mainly absorbing turbid media requires fine grid spacing to ensure numerical stability. In contrast, for strong forward scattering, also simulations on coarse grids can run stable.Fig. 4 Stability map for the D3Q26 discretization. Left: Contour plot of the discretized in-scattering term ωw7Φ7,7 with respect to anisotropy g and scattering albedo ω. Right: Stable and unstable regions with respect to optical depth, grid spacing and in-scattering for the liming case (see text). The limit between the stable and unstable regions is given by Eq. (30).
3. Results and discussion
3.1. Simulation cases
To validate the LBM, simulations of radiative transfer in an absorbing, scattering and non-emitting turbid media were carried out and compared to Monte Carlo simulations. The simulation domain is a cubic enclosure with one emitting surface at which monochromatic radiation of mean intensity J0 enters the domain. The other five surfaces are assumed to be black and cold. The domain, as sketched in Fig. 5, was assumed to initially dark. As a fist case, diffuse emission of radiation from the emitting boundary is considered. In a second case, the radiation source is modified in order to emit collimated radiation of intensity Ic. For this case, the source is changed to a square-shaped area around the origin O with an edge length of 0.2 L0 from which the collimated radiation enters the domain normal to the emitting surface.
Fig. 5 Sketch of the simulated domain with evaluated lines and planes. The origin of the coordinate system is marked with O. Radiation is emitted from the wall at x = 0 (yellow). Lines Ly are located at y1 = 0, y2 = 0.2, y3 = 0.4 and y4 = 0.48. Lines Lx are located at x1 = 0.25, x2 = 0.50 and x3 = 0.75. Planes Pz are located at z1 = 0 and z2 = −0.25.
For the D3Q26 LBM, the spatial resolution was set to N = 101 nodes in each direction after a grid independence test. Steady-state was assumed to be reached if the maximum local error in the domain (see Eq. (29)) was below a threshold of 1E–06. The reference MC solutions were obtained with an Open source solver [52], which is documented elsewhere [53]. For the MC simulations, the cubic domain was splitted into 1003 cubic voxels. In total, 1E07 photons were tracked in each simulation. The number of photons per voxel Nvox was normalized to the incoming number of photons per surface area so that the dimensionless mean intensity becomes to
In total 20 different parameter combinations of optical depth, anisotropy factor and scattering albedo were simulated for both cases under consideration.
3.1.1. Diffuse emission
At the boundary, the relation between the mean intensity and the specific intensity is given by
In the discretized momentum space of the LBM, Eq. (32) is expressed by Eq. (26) with the restriction of Ii(θi) = 0 for θi ≤ 0. Moreover, for diffuse radiation, the specific intensity per unit solid angle is constant so that Iiwi = const. applies for all directions i with θi > 0. The discrete specific intensities at the emitting wall can be calculated by solving the equation system that arises from the conditions as mentioned. At the cold walls, the boundary condition is Ii(θi) = 0 for all directions. For the MCM the boundary condition at the emitting boundary is given by Eq. (32). At the cold walls, photons are fully absorbed.
3.1.2. Collimated radiation
Collimated radiation is added to the LBM model by using a flux splitting approach [33]. Here, it is distinguished between collimated and diffuse radiation, which both are connected by an one-way coupling since diffuse radiation can not be focused again by scattering. Hence, the transport of collimated radiation is affected by extinction only as given by Lambert’s Law for which an analytic solution is known.
The extinction of collimated radiation is accompanied with an increase of the diffuse specific intensity via in-scattering. This can be included by adding a source-term Sc,i to Eq. (25), which is
where ΔIc,D = Ic,D(xD + Δxc,D) − Ic,D(xD). The discrete scattering phase function is represented by a vector Φc,i, which is calculated as described in section 2.2.3. The boundary condition was set to JD = 1 for the LBM. For the MCM, the angle of emission was fixed to the surface normal.3.2. Comparison with Monte Carlo simulations in scattering turbid media
Simulation results from the LBM and the MCM for the mean intensity in a scattering turbid media (ω = 1) are exemplarily shown for g = 0.85 in Fig. 6. The predictions by the LBM are in good agreement to results obtained by the MCM, although the momentum space is discretized in the LBM with only 26 abscissas, while for the MCM no angular discretization is required. It should be noted that for the case of collimated radiation (right-hand side of Fig. 6) the intensity profiles on lines Ly3 and Ly4 result from scattered radiation only, which indicates that the computation of the discretized scattering phase function by the iterative method leads to accurate results. However, if the optical depth decreases, the deviation between the LBM prediction and the MCM increases, what can be particularly noticed on line Ly1. Results of comparable consistency to those shown in Fig. 6 were found for lower anisotropy (g = 0.7, g = 0), but, in contrast, for strong anisotropy (g = 0.95) the LBM was found to become less accurate.
Fig. 6 Profiles of the dimensionless mean intensity JD on plane Pz1 along lines Ly1 (blue), Ly3 (red) and Ly4 (green) for g = 0.85 and different optical depths in a pure scattering turbid medium (ω = 1). Solid lines: lattice Boltzmann, symbols: Monte Carlo. Left: diffuse emission. Right: collimated radiation (see text for details).
A more extensive comparison of results obtained by the LBM and the MCM for the cases of diffuse and collimated radiation, different anisotropy factors and optical depths is shown in Fig. 7. For anisotropy factors g = 0.7 and g = 0.85 the data from LBM and MCM match well, independent from τ. Also for isotropic scattering, good agreement of results was achieved. For strong anisotropy (g = 0.95), the predicted intensity by the LBM increasingly deviates from the MCM for decreasing optical depth.
Fig. 7 Simulation data of the mean intensity obtained by the MCM and LBM for different levels of anisotropy. Data is evaluated on lines Ly1 – Ly4 and Lx1 – Lx3 of planes Pz1 and Pz2 for both methods. Colors symbolize different optical depths: τ = 100 (red), τ = 50 (blue), τ = 33.333 (green). Solid line: bisector.
The results indicate that numerical errors depend on anisotropy and optical depth for which the following error sources are conceivable. A first error is given by the linearization of the RTE during the discretization (see section 2.2.1). It is known that intensity is decreased exponentially due to extinction. By linearization, the decay of intensity and consequently also the amount of energy being in-scattered into a certain solid angle are overestimated. In case of isotropy or weak anisotropy, the linearization error overestimates sideward scattering of radiation while for strong anisotropy, focussing of radiation due to forward scattering occurs. Hence, the extent of the linerization is basically determined by the optical depth and the grid spacing while its effect on the spatial distribution of the mean intensity is determined by anisotropy. If grid spacing or optical depth is reduced, the linearized approximation should match the exponential decay of intensity better, but, however, the results shown in Figs. 6 and 7 do not reflect the expected behavior.
To investigate this phenomenon, ray effects caused by the propagation of the specific intensity on the lattice have to be considered. For that, the assumption shall be made that specific intensity propagates into a solid angle dΩ and further that τΔxD < 1 holds between two nodes of the lattice (as it is the case here for τ = 50 and τ = 33.333). In that case only a portion of the specific intensity is scattered and the remaining energy propagates by free streaming. On the lattice, however, there is no difference between forward scattering and free streaming in a certain direction. In consequence, the solid angle in which the radiation propagates artificially decreases as sketched in Fig. 8, which again leads to a ray effect. From a physical point of view, the ray effect is tantamount to a violation of Liouville’s theorem, which states conservation of intensity is phase space [5]. Strong anisotropic scattering pronounces the ray effect, while it is counteracted by low and medium anisotropy, high optical depth and coarse grid spacing. Another possible way to reduce ray effects is a better angular resolution and usage of higher order angular schemes [12, 15]. Similar to LBM for fluids, higher order angular schemes would also potentially allow the applicability of the method to higher Knudsen numbers [19], but at the cost of increasing computational effort.
Fig. 8 Schematic representation of the ray effect. In case of free streaming of radiation, the solid angle dΩ is artificially reduced to dΩ* while no radiation streams into the surrounding space.
The total numerical error, thus, results from two error sources, namely linearization error and ray effect. The extent of the total error results from an interplay of anisotropy, optical depth and grid spacing. Table 2 shows the root mean square error (RMSE) between the LBM and the MCM after outlier detection. Data was obtained by evaluation of lines Ly1 – Ly4 and Lx1 – Lx3 of planes Pz1 and Pz2 for both simulated cases described in section 3.1 and each data set consisted 2400 data points before elimination of outliers. Outliers were defined to be outside the interval [Q1 − a(Q3 − Q1), Q3 + a(Q3 − Q1)], where Qi are the quartiles of each data set and the parameter a was set to a = 3.
Table 2. Root mean square error of the mean intensity obtained by the LBM compared to the MCM. Data was evaluated along lines Ly1 – Ly4 and Lx1 – Lx3 of planes Pz1 and Pz2 (see Fig. 5).
Except from the isotropic case, decreasing the optical depth and, thus, the linearization error, did not result into more accurate prediction of the LBM. Moreover, amplification of the ray effect by strong anisotropy lead to an excessive increase of the total error, which partially is compensated by high optical depths. This results indicate that the ray effect is the dominating error source for the simulation of a scattering turbid media. Hence, to improve the accuracy of the method, future work should deal primarily with this effect. A numerical correction could be reached by implementation of artificial scattering terms, which guarantee the conservation of the specific intensity in phase space. The iterative discretization method (section 2.2.3) potentially allows to design discrete artificial scattering kernels for this purpose.
3.3. Comparison with Monte Carlo simulations in scattering and absorbing turbid media
The analysis made so far shall be extended to radiation transfer in a scattering, absorbing and non-emitting turbid medium. Table 3 shows the RMSE of results obtained by the LBM relative to the MCM. Similar to the simulation of radiation transfer in a non-absorbing medium, the ray effect acts as the dominant error source, even if the computation of absorption is affected by the linearization error only. This result was consistent for both simulated cases (see section 3.1).
Table 3. Root mean square error of the simulated mean intensity in a scattering and absorbing turbid medium obtained by the LBM compared to the MCM. Data was evaluated along lines Ly1 – Ly4 of planes Pz1 and Pz2 (see Fig. 5).
The observation that the ray effect dominates the numerical error is plausible since for the present case, the scattering albedo was set to ω = 0.75 and thus, scattering dominated over absorption. It can be assumed that for a low scattering albedo the linearization error becomes more dominant, although this was not investigated in the present study and the dependency of the different numerical error on the scattering albedo is not finally clarified. By comparison of tables 2 and 3 it can be seen that the dependency of the numerical error on the optical depth is less pronounced for the non-absorbing case. Hence, it can be concluded that the error caused by the ray effect is increasingly damped by absorption the greater the optical depth becomes.
The accurate numerical calculation of absorption requires fine grid spacing or low optical depths. Alternatively, further improvement of the method could be reached by using a semi-analytic approach, similar to the modeling of collimated radiation as given by Eqs. (33) and (34). Since this approach was also used to calculate the in-scattering from collimated radiation, a semi-analytic hybrid potentially neutralizes the linearization error.
3.4. Computational effort of the lattice Boltzmann method
The impact of simulation parameters on the computational effort is depicted in Fig. 9. For pure scattering, the number of iterations to reach steady state increases almost linearly with increasing optical depth for all applied anisotropy factors. However, the choice of the anisotropy factor dominates the computational effort, which is the highest for isotropic scattering. This result is expectable since isotropy corresponds with the slowest propagation of intensity on macroscopic scales. In contrast, for ω = 0.75 the increase of the optical depth accelerates convergence independent from anisotropy. It can be expected that this observation changes drastically in case of absorbing and emitting participating media, since isotropic emission of radiation decreases the speed of radiative transfer on macroscopic length scales.
Fig. 9 Number of iterations at steady-state with respect to optical depth and anisotropy for ω = 1 (blue) and ω = 0.75 (red). Dashed lines: linear regression.
4. Conclusion
A new lattice Boltzmann method for the simulation of three-dimensional radiation transport was developed and validated with Monte Carlo simulations. The method was used for the simulation of radiation transport in turbid media containing either isotropic and anisotropic scatters, as well as absorbing particles. To the best of the authors knowledge, this is the first time that a lattice Boltzmann method was successfully applied for the simulation of anisotropic radiation transfer in three dimensions. Moreover, for the approximation of Mie-scattering, a new iterative algebraic technique for the discretization of the scattering phase function was introduced, which ensures full conservation of energy and the mean scattering angle after discretization.
It was found that the accuracy of the lattice Boltzmann method is affected by two sources of error, namely linearization error and ray effect. For the investigated simulation cases, the latter dominates the total numerical error, independently of optical properties of the turbid media. However, it can be assumed that the dominating error type depends strongly on the scattering albedo. Based on the analysis of numerical errors, the implementation of an artificial scattering kernel and semi-analytic computation of extinction was suggested to further improve the accuracy of the method.
The development of a lattice Boltzmann method for radiation transport is motivated by the interpretation of the RTE as a Boltzmann-type transport equation in the framework of statistical physics. Based on this interpretation, the RTE can be discretized analogous to the discretization of the Boltzmann equation in established lattice Boltzmann methods for fluid dynamics if characteristics of radiation, in particular the constant speed of light, are considered. The establishment of a consistent methodological set of computational methods seems to be promising for solving multiphysical problems with high accuracy in engineering applications.
Acknowledgments
The authors are grateful to Antonio Delgado and Giovanni Luzi (LSTM, FAU Erlangen-Nuremberg) for the helpful discussions and productive collaboration.
References and links
1. M. R. Hoffmann, S. T. Martin, W. Choi, and D. W. Bahnemann, “Environmental applications of semiconductor photocatalysis,” Chem. Rev. 95(1), 69–96 (1995). [CrossRef]
2. A. D. Jassby and T. Platt, “Mathematical formulation of the relationship between photosynthesis and light for phytoplankton,” Limnol. Oceanogr. 21(4), 540–547 (1976). [CrossRef]
3. P. H. C. Eilers and J. H. C. Peters, “A model for the relationship between light intensity and the rate of photosynthesis in phytoplankton,” Numer. Heat Transfer Part B 57(2), 126–146 (2010).
4. J. Oxenius, Kinetic Theory of Particles and Photons (Springer, 1986). [CrossRef]
5. R. L. Liboff, Kinetic Theory. Classical, Quantum and Relativistic Descriptions (Springer, 2003).
6. M. F. Modest, Radiative Heat Transfer (Academic, 2013).
7. E. Meinköhn, G. Kanschat, R. Rannacher, and R. Wehrse, “Numerical methods for multidimensional radiative transfer,” Reactive Flows, Diffusion and Transport, 485–526 (Springer2007). [CrossRef]
8. P. J. Coelho, “Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media,” J. Quant. Spectrosc. Radiat. Transf. 145(16), 121–146 (2014). [CrossRef]
9. K. Mitra and S. Kumar, “Development and comparison of models for light-pulse transport through scatteringabsorbing media,” Appl. Opt. 38, 188–196 (1999). [CrossRef]
10. Z. Guo and S. Kumar, “Discrete ordinates solution of short pulse laser transport in two-dimensional turbid media,” Appl. Opt. 40(19), 3156–3163 (2001). [CrossRef]
11. Z. Guo and S. Kumar, “Three-dimensional discrete ordinates method in transient radiative transfer,” J. Thermophys Heat Transfer 16(3), 289–296 (2002). [CrossRef]
12. J. C. Chai, H. S. Lee, and S. V. Patankar, “Ray effect and false scattering in the discrete ordinate method,” Numer. Heat Transfer Part B 24(4), 373–389 (1993). [CrossRef]
13. P. J. Coelho, “The role of ray effects and false scattering on the accuracy of the standard and modified discrete ordinates methods,” J. Quant. Spectrosc. Radiat. Transf. 73, 231–238 (2002). [CrossRef]
14. B. Hunter and Z. Guo, “Comparison of quadrature schemes in DOM for anisotropic scattering radiative transfer analysis,” Numer. Heat Transfer Part B 63, 485–507 (2013). [CrossRef]
15. B. Hunter and Z. Guo, “Numerical smearing, ray effect, and angular false scattering in radiation transfer computation,” Int. J. Heat Mass Transfer 81, 63–74 (2015). [CrossRef]
16. D. Hänel, Molekulare Gasdynamik (Springer, 2004).
17. S. Chen, S. P. Dawson, G. D. Doolen, D. R. Janecky, and A. Lawniczak, “Lattice methods and their applications to reacting systems,” Comput. Chem. Eng. 19(6–7), 617–649 (1995). [CrossRef]
18. S. Chen and G. D. Doolen, “Lattice Boltzmann method for fluid flows,” Annu. Rev. Fluid Mech. 30, 329–364 (1998). [CrossRef]
19. S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Oxford University, 2001).
20. A. J. C. Ladd and R. Verberg, “Lattice Boltzmann simulation for particle fluid suspensions,” J. Stat. Phys. 104 (5/6), 1191–1251 (2001). [CrossRef]
21. D. Yu, R. Mei, L. S. Luo, and W. Shyy, “Viscous flow computations with the method of lattice Boltzmann equation,” Prog. Aerosp. Sci. 39, 329–367 (2003). [CrossRef]
22. C. K. Aidun and J. R. Clausen, “Lattice Boltzmann method for complex flows,” Annu. Rev. Fluid Mech. 42, 439–472 (2010). [CrossRef]
23. D. Anderl, M. Bauer, C. Rauh, U. Rüde, and A. Delgado, “Numerical simulation of adsorption and bubble interaction in protein foams using a lattice Boltzmann method,” Food Funct. 5, 755–763 (2014). [CrossRef] [PubMed]
24. D. Anderl, S. Bogner, C. Rauh, U. Rüde, and A. Delgado, “Free surface lattice Boltzmann with enhanced bubble model,” Comput. Math. Appl. 67 (2), 331–339 (2014). [CrossRef]
25. D. Anderl, M. Bauer, C. Rauh, U. Rüde, and A. Delgado, “Numerical simulation of bubbles in shear flow,” Proc. Appl. Math. Mech. 14 (1), 667–668 (2014). [CrossRef]
26. D. Anderl, M. Bauer, C. Rauh, U. Rüde, and A. Delgado, “Foam flow through channels with the lattice Boltzmann method,” Comput. Math. Appl. (2014, submitted).
27. S. Palpacelli and S. Succi, “Quantum lattice Boltzmann simulation of expanding Bose-Einstein condensates in random potentials,” Phys. Rev. E 77, 066708 (2008). [CrossRef]
28. A. Nabovati, D. P. Sellan, and C. H. Amon, “On the lattice Boltzmann method for phonon transport,” J. Comput. Phys. 230, 5864–5876 (2011). [CrossRef]
29. P. Asinari, S. C. Mishra, and R. Borchiellini, “A lattice Boltzmann formulation to the analysis of radiative heat transfer problems in a participating medium,” Numer. Heat Transfer Part B 57(2), 126–146 (2010). [CrossRef]
30. Y. Ma, S. K. Dong, and H. P. Tan, “Lattice Boltzmann method for one-dimensional radiation transfer”, Phys. Rev. E 84, 016704 (2011). [CrossRef]
31. A. F. DiRienzo, P. Asinari, R. Borchiellini, and S. C. Mishra, “Improved angular discretization and error analysis of the lattice Boltzmann method for solving radiative heat transfer in a participating medium,” Int. J. Numer. Methods Heat Fluid Flow , 21(5), 640–662 (2011). [CrossRef]
32. H. Bindra and D.V. Patil, “Radiative or neutron transport modeling using a lattice Boltzmann equation framework,” Phys. Rev. E 86, 016706 (2012). [CrossRef]
33. S. C. Mishra and R. R. Vernekar, “Analysis of transport of collimated radiation in a participating media using the lattice Boltzmann method,” J. Quant. Spectrosc. Radiat. Transf. 113(16), 2088–2099 (2012). [CrossRef]
34. R. R. Vernekar and S. C. Mishra, “Analysis of transport of short-pulse radiation in a participating medium using lattice Boltzmann method,” Int. J. Heat Mass Transfer 77, 218–229 (2014). [CrossRef]
35. Y. Zhang, H. Yi, and H. Tan, “One-dimensional transient radiative transfer by lattice Boltzmann method,” Opt. Express 21(21), 24532 (2013). [CrossRef] [PubMed]
36. Y. Zhang, H.-L. Yi, and H.-P. Tan, “The lattice Boltzmann method for one-dimensional transient radiative transfer in graded index gray medium,” J. Quant. Spectrosc. Radiat. Transf. 137, 1–12 (2014). [CrossRef]
37. Y. Zhang, H.-L. Yi, and H.-P. Tan, “Short-pulsed laser propagation in a participating slab with Fresnel surfaces by lattice Boltzmann methods,” Int. J. Heat Mass Transfer 80, 717–726 (2015). [CrossRef]
38. Y. Zhang, H.-L. Yi, and H.-P. Tan, “Lattice Boltzmann method for short-pulsed laser transport in multi-layered medium,” J. Quant. Spectrosc. Radiat. Transf. 155, 75–89 (2015). [CrossRef]
39. R. McCulloch and H. Bindra, “Coupled radiative and conjugate heat tranfer in participating media using lattice Boltzmann methods,” Comput. Fluids 124, 261–268 (2015).
40. A. Mink, G. Thäter, H. Nirschl, and M. J. Krause, “A 3D lattice Boltzmann method for light simulation in participating media,” J. Comput. Sci. (2016, in press). [CrossRef]
41. G. C. Pomraning, Linear Kinetic Theory in Stochastic Mixtures (World Scientific, 1991).
42. M. Jonasz and G. R. Fournier, Light Scattering by Particles in Water (Academic, 2007).
43. A. B. Davis and A. Marshak, “Photon propagation in heterogeneous optical media with spatial correlations: enhanced mean-free-paths and wider-than-exponential free-path distributions,” J. Quant. Spectrosc. Radiat. Transf. 84, 3–34 (2004). [CrossRef]
44. R. Koch, W. Krebs, S. Wittig, and R. Viskanta, “Discrete ordinates quadrature schemes for multidimensional radiative transfer,” J. Quant. Spectrosc. Radiat. Transf. 53(4), 353–372 (1995). [CrossRef]
45. X. He and L. S. Luo, “Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation,” Phys. Rev. E 56, 6811–6818 (1997). [CrossRef]
46. T. Abe, “Derivation of the lattice Boltzmann method by means of the discrete ordinate method for the Boltzmann equation,” J. Comput. Phys. 131, 241–246 (1997). [CrossRef]
47. W. Unno and E. A. Spiegel, “The Eddington approximation in the radiative heat equation,” Publ. Astron. Soc. Jap. 18(2), 85–95 (1966).
48. V. I. Lebedev and D. N. Laikov, “A quadrature formula for the sphere of the 131st algebraic order of accuracy,” Dokl. Math. 59(3), 477–481 (1999).
49. B. Hunter and Z. Guo, “Conservation of asymmetry factor in phase function discretization for radiative transfer analysis in anisotropic scattering media,” Int. J. Heat Mass Transfer 55(5–6), 1544–1552(2012). [CrossRef]
50. B. Hunter and Z. Guo, “Phase-function normalization in the 3-D discrete-ordinates solution of radiative transfer -Part 1: conservation of scattered energy and asymmetry factor,” Numer. Heat Transfer Part B 62, 203–222 (2012). [CrossRef]
51. J. H. Ferziger and M. Peric, Computational Methods for Fluid Dynamics (Springer, 2002). [CrossRef]
52. S. L. Jacques, T. Li, and S. Prahl, “mcxyz.c, a 3D Monte Carlo simulation of heterogeneous tissues,” http://omlc.org/software/mc/mcxyz/index.html.
53. L.-H. Wang, S. L. Jacques, and L.-Q. Zheng, “MCML - Monte Carlo modeling of photon transport in multilayered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]
