Unified Gas-kinetic Wave-Particle Methods III: Multiscale Photon Transport

Unified Gas-kinetic Wave-Particle Methods III: Multiscale Photon Transport

Weiming Li liweiming@csrc.ac.cn Chang Liu cliuaa@connect.ust.hk Yajun Zhu zhuyajun@mail.nwpu.edu.cn Jiwei Zhang jwzhang@csrc.ac.cn Kun Xu makxu@ust.hk Applied and Computational Mathematics Division, Beijing Computational Science Research Center, Beijing 100193, China National Key Laboratory of Science and Technology on Aerodynamic Design and Research, Northwestern Polytechnical University, Xi’an, Shaanxi 710072, China Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China HKUST Shenzhen Research Institute, Shenzhen 518057, China
Abstract

In this paper, we extend the unified gas-kinetic wave-particle (UGKWP) method to the multiscale photon transport. In this method, the photon free streaming and scattering processes are treated in an un-splitting way. The duality descriptions, namely the simulation particle and distribution function, are utilized to describe the photon. By accurately recovering the governing equations of the unified gas-kinetic scheme (UGKS), the UGKWP preserves the multiscale dynamics of photon transport from optically thin to optically thick regime. In the optically thin regime, the UGKWP becomes a Monte Carlo type particle tracking method, while in the optically thick regime, the UGKWP becomes a diffusion equation solver. The local photon dynamics of the UGKWP, as well as the proportion of wave-described and particle-described photons are automatically adapted according to the numerical resolution and transport regime. Compared to the -type UGKS, the UGKWP requires less memory cost and does not suffer ray effect. Compared to the implicit Monte Carlo (IMC) method, the statistical noise of UGKWP is greatly reduced and computational efficiency is significantly improved in the optically thick regime. Several numerical examples covering all transport regimes from the optically thin to optically thick are computed to validate the accuracy and efficiency of the UGKWP method. In comparison to the -type UGKS and IMC method, the UGKWP method may have several-order-of-magnitude reduction in computational cost and memory requirement in solving some multsicale transport problems.

keywords:
Radiative transfer equations, Diffusion equation, Wave-particle formulation, Monte Carlo particle method, Unified gas-kinetic scheme
journal: Journal of LaTeX Templates

1 Introduction

The radiative transfer equation describes photon propagation in the background medium and has important applications in the fields of astrophysics davis2012 (), atmospheric physics marshak20053d (), optical imaging klose2002optical () and so on. In this paper, we focus on the gray radiative transfer equation with isotropic scattering, which reads

(1)

where is the specific intensity in time , space , and angle . Here is the speed of light, is the scattering coefficient, is the absorption coefficient, and is an internal source of photons.

There are typically two categories of numerical methods for solving the radiative transfer equations. The first category consists of the deterministic methods with different ways of discretizing and modeling, such as the discrete ordinate method (DOM) hunter2013comparison (); coelho2014advances (); chen2015 (); roos2016conservation () and the moment methods frank2006partial (); carrillo2008numerical (); vikas2013radiation (); alldredge2016approximating (). The second category is the stochastic methods, for example, the Monte Carlo method fleck1971 (); lucy1999computing (); hayakawa2007coupled (). The Monte Carlo method is a very popular method for solving the radiative transfer problems. In comparison with the deterministic methods, it is more efficient in optically thin regime especially for the multidimensional cases, and it does not suffer from the ray effect. However, as a particle method, it unavoidably has statistical noise. Also, it becomes inefficient when it comes to diffusive regime. In diffusive regime where the mean free path is small, photons may go through a huge number of scatterings during their lifetimes. Direct simulation of each scattering process for all particles makes the Monte Carlo method very expensive for capturing the diffusive solutions. On the other hand, in the diffusive regime the photon transport process could be well described by the diffusion equation, which could be solved efficiently. Based on this observation, many hybrid methods have been developed in order to improve the overall efficiency in different regimes fleck1984random (); giorla1987random (); densmore2007hybrid (); densmore2012hybrid (), where the Monte Carlo method is used in the optically thin regions and the diffusion equation is applied to the optically thick regions. However, as far as we know, there is still no unified principle for accurate domain decomposition for different regimes.

Another approach towards releasing the stiffness issue in the diffusive regime is to develop asymptotic-preserving (AP) DOM klar1998asymptotic (); naldi1998numerical (); jin1999efficient (); jin2000uniformly (); mieussens2013asymptotic (); sun2015asymptotic1 (); sun2015asymptotic (); sun2017multidimensional (); sun2017implicit (); sun2018asymptotic () . One of the examples is the unified gas-kinetic scheme (UGKS), which couples the particles’ transport and collision process using a multiscale flux function obtained from the integral solution of the kinetic model equation. The cell size and time step are not restricted by the mean free path and mean collision time. It was developed initially in the field of rarefied gas dynamics xu2010unified (); xu2015direct (); jiang () and has been applied to the field of radiative transfer mieussens2013asymptotic (); sun2015asymptotic1 (); sun2015asymptotic (); sun2017multidimensional (); sun2017implicit (); sun2018asymptotic (), plasma transport liu2017unified (), and disperse multi-phase flow LIU2019264 (). Since it is a discrete ordinate based numerical scheme, it has no statistical noise, but unavoidably suffers from the ray effect.

The UGKS provides a general methodology to construct multiscale simulation methods for the transport equations xu2015direct (); xu-liu-pof (). It consists of two governing equations for the microscopic distribution function and macroscopic flow variables on the mesh size and time step scales, and a multiscale evolution solution for the interface distribution function. The time-dependent evolution solution at a cell interface covers the dynamics from the particle free transport to the hydrodynamic wave interaction with the variation of the ratio of the time step to the local particle collision time. The original UGKS is constructed based on the discrete velocity method (DVM) or DOM formulations xu2010unified (). However, the corresponding purely particle version and wave-particle version of UGKS can be constructed as well liu-zhu-xu (); zhu-liu-zhong-xu (). In this work, we develop a novel unified gas-kinetic wave-particle (UGKWP) method for the multiscale photon transport, which combines the advantages of UGKS and the particle method. To facilitate understanding, we first describe a purely particle-based unified gas-kinetic particle (UGKP) method. In the UGKP method, the photons are described by the particle transport and collision, and this process is controlled by a multiscale integral formulation. More specifically, the Monte Carlo particle model is used to discretize the angular direction of the photon’s movement. Based on the particles’ transport nature in the discretized physical space, particles are categorized into two groups. Given a fixed time step, the free-stream particles are accurately tracked by following the trajectories of the simulation particles, while those particles that get scattered within the given time step are eliminated and re-sampled according to their macroscopic variables at the new time level. The fluxes across a cell interface from different type particles are taken into account to update cell averaged macroscopic variables. In such a way, the multiscale process is preserved by coupling particle free streaming and collision. Based on UGKP, a more efficient UGKWP method is proposed. Instead of representing all photons by simulation particles, a proportion of photon is represented by an analytical distribution function in UGKWP. Therefore, part of the evaluation of the fluxes, which were computed by particles in the UGKP method, could be done analytically with significant reduction in computational cost and statistical noise, especially in the diffusive regime. The multiscale flux function of the UGKS is precisely preserved in the UGKWP implementation. In the diffusive regime, the resulting algorithm would become a standard central difference scheme for the diffusion equation without any particle and DOM coordinate in the velocity space. In the optically thin regime, it gives a particle tracking method same as the Monte Carlo method. In the transition regime, the ratio of the time step over particle collision time determines the transport dynamics between the above two limits. The UGKWP is not a hybrid method of domain decomposition, but is a method to recover the multiscale modeling in UGKS xu2015direct () with the help of wave and particle in order to achieve a highest efficiency in the transport simulation in different regimes.

The rest of this paper is organized as follows. Section 2 briefly recalls the basic idea of the unified gas-kinetic scheme (UGKS) for the linear transport equation. Section 3 presents the UGKP method for linear photon transport. The UGKWP method is described in Section 4, while the extension to radiation-material coupling system is discussed in Section 5. In Section 6, numerical tests are presented to demonstrate the accuracy and efficiency of the UGKWP method. The final section is the conclusion.

2 Review of the UGKS for the linear transport equation

The unified gas-kinetic scheme (UGKS) was initially developed for the problems in the field of rarefied gas dynamics xu2010unified (); xu2015direct (), and has also been successfully applied to problems in radiative transfer under the finite volume framework mieussens2013asymptotic (); sun2015asymptotic1 (); sun2015asymptotic (); sun2017multidimensional (); sun2017implicit (); sun2018asymptotic (). In this section, we review the basic idea of the UGKS using the example of the multidimensional linear transport equation in a purely scattering medium.

Consider

(2)

which gives a non-dimensional equation

(3)

where . The non-dimensionalization is used same as that in mieussens2013asymptotic (), and is the ratio between the typical mean free path and the macroscopic length scale.

The UGKS is based on a finite volume framework. Consider a single cell . Denote as the area covered by this cell, as the interfaces of its boundary , as the outward normal of , and as the volume of cell . Define

to be the averaged specific intensity over the spatial cell, and

to be the averaged energy density function over a spatial cell. Under the finite volume framework, the discretizations of the two fundamental governing equations for microscopic and macroscopic variables in UGKS are,

(4)

and

(5)

where the microscopic and macroscopic flux terms are respectively

(6)

and

(7)

Based on the above microscopic and macroscopic governing equations, the key ingredient of the UGKS is the construction of the multiscale flux function by adopting the integral solution of the kinetic model equation (3). The integral solution of equation (3) along the characteristic line gives

(8)

which is used to construct the numerical fluxes in equation (4) and (5). The integral solution couples transport with particle collisions, and bridges the kinetic and the hydrodynamic scale dynamics. Note that the dynamics in the evolution equations of UGKS through equations (4), (5), and (8), is related to time step . The ratio of the time step to the local particle collision time determines the transport regime of the UGKS.

The numerical flux for microscopic and macroscopic variable updates is based on the piecewise linear initial reconstruction of and at the beginning of each time step. The discretization of is based on the DOM in the original UGKS. It has been proved in mieussens2013asymptotic () that when equals , the UGKS tends to a finite volume scheme which is consistent with free transport solution. More specifically, for 1D case with uniform grid, the scheme becomes

In the diffusion limit, with a uniform mesh the UGKS scheme becomes a standard central difference method for the limit diffusion equation as tends to . Specifically, for 1D case with uniform grid, in the diffusive limit the scheme is

More details on the asymptotic analysis of the UGKS for the radiative transfer equation can be found in mieussens2013asymptotic () and sun2015asymptotic1 ().

Following the methodology of the UGKS, we will construct two particle-based algorithms with multiscale transport property for recovering transport physics from the kinetic scale to the hydrodynamic scale, i.e., the unified gas-kinetic particle (UGKP) method and the unified gas-kinetic wave-particle (UGKWP) method. The UGKP is a purely particle based method whose computational cost and statistical noise remain the same in all flow regimes. The UGKWP method is an improvement of the UGKP method such that the distribution function is represented by the combination of simulation particles and an analytical distribution function, which makes it more efficient and less noisy than the UGKP in the diffusive regime. In both algorithms, for the kinetic scale particle free transport, the particle trajectories are tracked precisely; for those particles suffering collisions, their macroscopic variables will be updated and used to re-sample these particles subsequently. The multiscale particle methods for equation (8) are constructed through the tracking and re-sampling of particles with the help of updated macroscopic variables. We will first introduce the UGKP method in the next section.

3 The unified gas-kinetic particle (UGKP) method

In this section, we present a unified gas-kinetic particle (UGKP) method. Following the direct modeling methodology of UGKS, the evolution of microscopic simulation particle is coupled with the evolution of macroscopic energy in UGKP method. The method presented in this and the next section is based on a single linear transport equation,

(9)

In the following two subsections, we will introduce the multiscale evolution equations of particle and macroscopic energy. For simplicity, the method will be presented for the two-dimensional case on uniform Cartesian mesh. Its extension to non-uniform mesh and 3D is straightforward.

3.1 Multiscale Particle evolution

The evolution of the simulation particle follows the integral solution of the linear transport equation (9),

(10)

The local equilibrium can be expanded as

(11)

and the integral solution can be written as

(12)

The first order approximation gives

(13)

and the second order approximation gives

(14)

The reformulated integral solution Eq.(12) is the multiscale governing equation for particle evolution in UGKP method, which implies that the simulation particle has a probability of to free stream and has a probability of to get scattered. The scattered particle follows a velocity distribution , and the first order approximation Eq.(13) is used in this paper. Different from the traditional Monte Carlo method which tracks the trajectory of each particle as shown in Fig. 1, UGKP tracks particle trajectory until its first scattering time, and then evolves the velocity distribution function. The scattered particles will be temporarily removed and then get re-sampled at from their evolved distribution, as shown in Fig. 2. The free-streaming particles are named as the particles, and the scattered particles are named as the particles. The total energy of the particles is called energy. In UGKP we do not need to resolve each scattering process by dividing every time step into several sub time steps and therefore achieve a high efficiency compared to the traditional Monte Carlo method.

The first scattering time is defined as the time when the first scattering happens to the particle. According to Eq.(12), the time can be obtained for each particle as

where is a random number which follows a uniform distribution on . Simulation particles can be categorized into two groups based on . The first group of simulation particle satisfies

and does not suffer any collision during a time step , and this particle position is updated by

The second group of particle satisfies , those particles suffer one or multiple times of scattering during , with the assumption for simplification. Therefore, these particles will only be accurately tracked for and be re-sampled from at . is obtained from the evolved macroscopic energy. In the UGKP method, the evolution of particle and macroscopic energy are closely coupled. The multiscale evolution equation of macroscopic energy is given in the next subsection.


Figure 1: Particle trajectory in Monte Carlo method for with sub time steps , , and .

Figure 2: UGKP’s treatment of particle evolution for .

3.2 Multiscale macroscopic energy evolution

The update of the macroscopic variable follows the finite volume framework

(15)

where the flux is the diffusive flux and is the contribution of the free streaming flux. Based on (12), the diffusive flux is computed same as that in the UGKS,

We adopt the piecewise linear reconstruction for ,

(16)

and the explicit central difference discretization of ,

where , are the barycenters of cells and with the same interface . Then the diffusive flux is given by

(17)

where

The free transport flux is contributed by the particle free streaming, which is computed by

(18)

where and are respectively the weight and velocity of the particle , and is the set of all particles which are transported across interface during the free streaming process. It recovers the other part of UGKS flux

(19)

With the diffusive flux Eq. (17) and the free transport flux Eq. (18), the macroscopic energy can be updated by

(20)

With the updated , the energy of the scattered particles can be obtained and be used to re-generate particles,

(21)

where is the set of particles in cell at . In summary, the evolution of particle provides the free-stream flux for macroscopic energy evolution, and the updated energy provides the distribution for the corresponding particles. The evolution of microscopic and macroscopic quantities are closely coupled, and the algorithm of UGKP is shown in Fig. 4.


Figure 3: UGKP re-sampling process at .

Figure 4: Flowchart of UGKP.

4 Unified gas-kinetic wave-particle (UGKWP) method

4.1 The multiscale evolution of particles and macroscopic variable

In UGKP, the particles are re-sampled from an equilibrium distribution at . According to the integral solution (10), particles have a probability to get scattered in the next time step from to . We define the scattered particles in as particles and the non-scattering particles as particles. As the particles will be deleted after within , and their contribution to the free transport flux could be computed analytically, only the particles need to be re-sampled and they keep on free transport in the step from to . For the particles, we only need to store the distribution function instead of re-sampling particles from it. Based on this observation, we propose a more efficient unified gas-kinetic wave-particle (UGKWP) method.

The unified gas-kinetic wave-particle(UGKWP) method improves UGKP in two aspects:

  1. Only particles are sampled as shown in Fig. 5;

  2. The free transport flux contributed by can be calculated analytically.


Figure 5: UGKWP re-sampling process at .

The finite volume energy evolution follows

(22)

where the diffusive flux is computed using Eq. (17), and the free transport flux contributed by is calculated by

The free transport flux contributed from the particles in Eq. (22) is computed by particle tracking,

(23)

where is the set of all particles which are transported across the cell interface during the free streaming process. In the particle free transport process, the re-sampled particles will always have and keep free transport without collision in the whole time step. The other particles from will transport according to their individual generated . With the solution of , Eq.(21) is used to get , which is divided into

and

and only particles from will be generated and keep on free transport in the whole next time step.

The algorithm of UGKWP is shown in Fig. 6.


Figure 6: Flowchart of UGKWP.

4.2 Properties of the UGKWP algorithm

The unified gas-kinetic wave-particle method satisfies the following properties:

  1. The microscopic and macroscopic evolutions are consistently evaluated based on the wave-particle decomposition. The macroscopic energy density is the sum of the energy from the un-scattering particles and the scattering one . The evolution of particles can be evaluated analytically.

  2. In the diffusive limit, and , essentially no particle can be freely transported and survived within a time step. Therefore, no particle is re-sampled from the energy. The algorithm becomes a time-explicit central difference solver of the diffusion equation. For instance, in 1D case with uniform grid, the scheme tends to the following limiting equation

    Note that in the diffusion regime, different from DOM method there is no discrete velocity space discretization in UGKWP.

  3. In the free transport limit, and , each particle is traced exactly by free transport with probability 1. In this case, the UGKWP could recover the exact solution for individual particle.

5 Extension to the coupled system of gray radiative transfer and material energy evolution

This section extends UGKWP to solve the coupled system of gray radiative transfer and material temperature equation,

(24)

Define and , then the second equation can be re-written as

(25)

The implicit Monte Carlo method proposed by Fleck and Cummings in fleck1971 () has been shown to be an effective technique for solving non-linear, time-dependent, radiative transfer problems and is widely used in the radiative transfer community. Fleck’s implicit Monte Carlo method uses an effective scattering process to approximate the absorption and emission of radiation by the background medium. This treatment allows it to take larger time steps than that in a purely explicit method. Here the similar semi-implicit discretization for material temperature evolution will be employed. Specifically, Eq. (25) is discretized by

which gives

(26)

With the definition

substituting Eq. (26) into Eq. (24) yields

An operator splitting approach is used to solve the above system subsequently by the linear kinetic equation

(27)

the radiation energy exchange,

(28)

and the update of material energy. Here Eq. (27) is solved using the algorithm introduced in Section 4. The exact solution of Eq. (28) is

(29)

where is solved from Eq. (27) by UGKWP. The exponential decay term in Eq. (29) is implemented by modifying the weight of particles, while the source term in Eq. (29) is added to energy. Afterwards, the energy change of particles is summed up, and the material temperature and are updated by energy conservation. The above UGKWP can be further improved if the technique in sun2015asymptotic () is incorporated into the current scheme, where both macroscopic equations for the radiation energy and material energy are solved iteratively first before updating the radiation intensity .

6 Numerical Experiments

In this section, we present six numerical examples to validate the proposed UGKWP method. In all these examples, the UGKWP method obtains results which are in good agreement with that of the Monte Carlo method. At the same time, the UGKWP is more efficient and less noisy compared with the Monte Carlo method in the near diffusive regime. All computations are performed in sequential code on a computer with Intel i7-8700K CPU and GB memory. In all simulations, the time step is determined by , with for 1D examples and for 2D examples.

6.1 Inflow into purely scattering homogeneous medium

We first consider the behaviour of the UGKWP method for purely scattering homogeneous medium. Tests in this section are for non-dimensional linear equation

defined on the semi-infinite spatial domain with an isotropic inflow condition imposed on the left boundary. For the numerical simulation, the spatial domain is taken to be . Inflow boundary condition is imposed at with the incoming specific intensity . The initial value is for all .

The results of both the UGKWP method and the Monte Carlo method are obtained using grids in space. As we are targeting to develop a method that automatically bridges the optically thin and optically thick regimes, the parameters cover the rarefied (), the intermediate (), and the diffusive () regimes, as defined in jin1998diffusive (). When is small, i.e. in the diffusive regime, the current method can use a much larger cell size and time step than the particle mean free path and collision time. Also, due to the exponential factor in re-sampling, fewer particles are used in UGKWP than the Monte Carlo method when is small, and UGKWP becomes much more efficient and less noisy than the Monte Carlo method in those cases. Table 1 compares the computational time of the Monte Carlo method and the UGKWP method for different , while Table 2 compares the maximum number of particles used. Notice that no particle is generated in UGKWP method when and the method only solves a diffusion equation. For , the Monte Carlo method takes more than half an hour while the UGKWP method takes only seconds, showing that the UGKWP method is much more efficient than the Monte Carlo method in the diffusive regimes.

MC UGKWP
s s
s s
s s
Table 1: Comparison of computational time of the Monte Carlo method and UGKWP for inflow into homogeneous medium under different parameters.
MC UGKWP
Table 2: Comparison of the maximum number of particles used by the Monte Carlo method and UGKWP for inflow into homogeneous medium under different parameters.

In Figure 7 the numerical results of the UGKWP method are compared with the solutions of the Monte Carlo method for different .

(a) Comparison of between UGKWP and the Monte Carlo solution for .
(b) Comparison of between UGKWP and the Monte Carlo solution for .
(c) Comparison of between UGKWP and the Monte Carlo solution for .
Figure 7: Numerical results of homogeneous medium.

In this example, the UGKWP method gives solutions which are almost identical with the Monte Carlo solution for all flow regimes. Also, note that the UGKWP solution is smooth for , while the Monte Carlo solution has statistical noise.

6.2 Inflow into purely scattering heterogeneous medium

We next consider the UGKWP method for purely scattering heterogeneous medium. For this example, the set-up for the computational domain as well as that of the initial and boundary conditions are the same as in the previous section. Tests in this section are for the dimensional linear equation

Both the UGKWP method and the Monte Carlo method use grids in space.

We calculate two test cases. For the first test case, we take . In this case, radiation passes through an optically thick medium near the left boundary and gets into the region with a gradually reducing of optically thickness. In Figure 8(a) we compare the solutions of the UGKWP method and the Monte Carlo method for and they agree with each other closely. Also, in this case where the medium is relatively optically thick overall, the UGKWP produces solutions much smoother than those given by the Monte Carlo method. For this test case, the Monte Carlo method uses a maximum of particles and takes seconds. The UGKWP method uses a maxium of particles and takes seconds. Therefore, the UGKWP method is much more efficient than the Monte Carlo method for this test case.

For the second test case, we take , implying that the medium gets more and more optically thick from the left to right boundaries. In Figure 8(b) the numerical result of the UGKWP method is compared with the solution of the Monte Carlo method at and they are almost identical except for the statistical noise. The Monte Carlo method uses a maximum of particles and spends seconds while the UGKWP method uses a maximum of particles and spends seconds in computing this case. In this test case the medium is relatively optically thin overall, and the UGKWP method has similar time cost as the Monte Carlo method.

(a) Comparison of between UGKWP and the Monte Carlo solution in thick-to-thin heterogeneous medium.
(b) Comparison of between UGKWP and the Monte Carlo solution in thin-to-thick heterogeneous medium.
Figure 8: Numerical results of heterogeneous medium.

6.3 Marshak wave problem

This section studies the Marshak wave problem where radiation is coupled with material medium. Consider the system

for semi-infinite domain with a computational one . We take absorption coefficient to be , the speed of light to be , the parameter and the specific heat to be . The initial material temperature is set to be and initially material and radiation energy are at equilibrium. A constant isotropic inflow radiation intensity with an equivalent radiation temperature of is imposed on the left boundary. Both the UGKWP method and the Monte Carlo method use grids in space. For this example, the implicit Monte Carlo method uses a maximum of particles, while the UGKWP method uses a maximum of particles. It takes the implicit Monte Carlo method seconds to compute until and the UGKWP method seconds.

In Figure 9(a), the computed radiation wave front at time , and are given, while Figure 9(b) presents the computed material temperature. The solutions of the UGKWP method are shown to be consistent with those of the Monte Carlo method for both radiation and material temperature.

(a) Comparison of radiation temperature between UGKWP and the Monte Carlo solution for Marshak wave problem.
(b) Comparison of material temperature between UGKWP and the Monte Carlo solution for Marshak wave problem.
Figure 9: Numerical results of Marshak wave problem.

6.4 Line-source problem in purely scattering homogeneous medium

The next problem we look at is the line-source problem in purely scattering medium. In this test case, we consider equation (3) for and . The spatial 2D computational domain is . The initial density distribution is

which means that initially all particles are concentrated at and they spread out over the time. This is a particularly difficult problem for the based methods because they suffer from severe ray effect. Previous studies have used this problem for comparisons between methods brunner2002forms () and as a test case for schemes to mitigate the ray effect camminady2019ray ().

Both the UGKWP method and the Monte Carlo method use cells in space. The UGKWP method uses particles while the Monte Carlo method uses particles and the same spatial mesh as the UGKWP method. Because this example considers the kinetic regime, the number of particles used in the UGKWP method is almost the same as that used in the Monte Carlo method and the computation time is similar, with seconds for the UGKWP method and seconds for the Monte Carlo method. The UGKWP method is slightly more expensive due to its additional overheads.

Fig.10 compares the solution of between the UGKWP method, the Monte Carlo method and the method for time . For the computation, we use an equal-weight quadrature set, cells for the spatial mesh and the DOM-UGKS. The UGKWP solution is consistent with the Monte Carlo solution. Both solutions display a sharp spike near the wave front and has a smooth nonzero region behind the wave front, and neither suffers from the ray effect. On the other hand, the solution is qualitatively incorrect due to the ray effect. This example demonstrates that the UGKWP method preserves the advantage of the Monte Carlo method in the rarefied regime without suffering from the ray effect.

(a) Contour plot of the solution of by the Monte Carlo method as a function of the spatial coordinate in the line-source problem at .
(b) Contour plot of the solution of by the UGKWP method as a function of the spatial coordinate in the line-source problem at .
(c) Contour plot of the solution of by the method as a function of the spatial coordinate in the line-source problem at .
(d) Slice at of the solution of by the Monte Carlo method and the UGKWP method as a function of in the line-source problem at .
Figure 10: Numerical results of the line-source problem.

6.5 Emitting isotropic source in lattice medium

In this section we study a problem with multiple medium mcclarren2010robust () by considering the following equation

The diffusion term is solved by using a backward Euler method. The two-dimensional physical space consists of a set of squares belonging to a strongly absorbing medium and a background of weakly scattering medium. The specific layout of the problem is given in Fig.11, where the light green regions and the white region are purely scattering with and ; the red regions are pure absorbers with and . In the white region in the center there is an isotropic source . The initial material temperature is and initially the radiation and material temperature are at equilibrium. We take , and to simulate the case.

Figure 11: Layout of the lattice problem.

Both the UGKWP method and the Monte Carlo method use a mesh of in physical space. The results at time are summarized in Fig. 12 with the natural logarithm of density distribution shown in contour and slice, and they are in agreement with each other. This is a test case which includes material in both the kinetic and the diffusive regimes, but most of the time particles travel in the kinetic regime. The Monte Carlo method uses particles and took seconds, while the UGKWP method uses particles and took seconds, showing the UGKWP method to be more efficient than the Monte Carlo method for this case where there is mixed medium.

(a) Contour plot of the solution of by the Monte Carlo method as a function of the spatial coordinate in the lattice problem at .
(b) Contour plot of the solution of by the UGKWP method as a function of the spatial coordinate in the lattice problem at .
(c) Slice at of the solution of by the Monte Carlo method and UGKWP as a function of in the lattice problem at .
(d) Slice at of the solution of by the Monte Carlo method and UGKWP as a function of in the lattice problem at .
Figure 12: Numerical results of the lattice problem.

6.6 A hohlraum problem

This examples considers a hohlraum problem similar to the one studied in mcclarren2010robust (). We study equation (24) with . The problem layout is given in Fig.13. This is a problem of purely absorbing medium where the absorption coefficient relies on the material temperature. As the material temperature varies, the values of the absorption coefficient cover a wide range, presenting challenges to numerical methods.

Figure 13: Layout of the hohlraum problem. In the blue regions , and . The white region is vacuum, .

On the left boundary there is an isotropic source with the specific intensity . Initially radiation and material temperature are at equilibrium and the initial material temperature is .

The UGKWP method uses grids in space while the Monte Carlo method use grid in space. We compute for . The Monte Carlo method uses a maximum of particles and takes hours. The UGKWP method uses a maximum of particles and takes seconds, therefore the UGKWP method is times more efficient than the Monte Carlo method for this test case. For this test case, the initial absorption coefficient on the boundary is extremely large, and the Monte Carlo method took an especially long time over the first time step when the inflow radiation first heats up the boundary material. The UGKWP method, on the other hand, has no such problems.

The authors of mcclarren2010robust () pointed out solution to this problem should contain two characteristics: the central block should be heated non-uniformly; and the region behind the block with respect to the source should have less radiation energy then regions within the line of sight of the source. Their studies showed that the diffusion equation could not capture these two characteristics. In Fig. 14 and 15, we compare the solution of the UGKWP method and the Monte Carlo method for energy density and material temperature at and they are all consistent with each other. Both solutions agree with each other such that there is less radiation to the right of the block. This shows that the UGKWP method keeps the accuracy of the Monte Carlo method and is more accurate than the diffusion approximation in problems which are essentially multiscale in nature.

(a) Contour plot of the solution of by the Monte Carlo method as a function of the spatial coordinate in the hohlraum problem at .
(b) Contour plot of the solution of by the UGKWP method as a function of the spatial coordinate in the hohlraum problem with at .
(c) Contour plot of the solution of by the Monte Carlo method as a function of the spatial coordinate in the hohlraum problem at .
(d) Contour plot of the solution of by the UGKWP method as a function of the spatial coordinate in the hohlraum problem with at .
Figure 14: Numerical results of the hohlraum problem I.
(a) Comparison of sliced solution at of the solution of by the UGKWP method and the Monte Carlo method as a function of coordinate in the hohlraum problem at .
(b) Comparison of sliced solution at of the solution of by the UGKWP method and the Monte Carlo method as a function of coordinate in the hohlraum problem at .
Figure 15: Numerical results of the hohlraum problem II.

7 Conclusion

In this paper, for the first time a unified gas-kinetic wave-particle method is proposed to simulate photon transport. Due to the particle wave decomposition and their dynamic coupling, the UGKWP method becomes a multiscale method to simulate transport physics with the most efficient way in different regimes. For the linear transport equation, this method recovers the solution of the diffusion equation in the optically thick limit without constraint on the time step being less than the photon’s mean collision time. At the same time, it gives the exact solution in the free transport regime. In the transition regime, both macroscopic wave and kinetic particle contribute to the radiative transfer, and their weights and coupling depend on the ratio of the time step to the local particle’s mean collision time. The UGKWP method is also extended to the coupled radiation-material system. With the inclusion of energy exchange, the UGKWP method can give excellent simulation results in different regimes. A few benchmark problems are tested to show the performance of the current scheme. The accuracy and efficiency of the UGKWP method are fully confirmed. For the multiscale transport, in the cases with the co-existence of different regimes the UGKWP may improve the efficiency on several-order-of-magnitude in comparison with the purely particle methods. The UGKWP takes also advantages of both particle and macroscopic solver. The UGKWP has totally removed the ray effect and and has a much improved efficiency in comparison with DOM-type UGKS. Based on the direct modeling of the transport physics in the time step scale xu2015direct (), such as the relationship between and in the current study, we have three versions of UGKS for the multiscale transport simulations, i.e., the DOM or DVM type, the purely particle formulation, and the wave-particle one. The advantages and disadvantages of these different discretization under the same UGKS framework will be further investigated, and the extension of these schemes to complex systems, such as radiative-hydrodynamics, plasma, multi-component system, and reactive flow, will be constructed.

Acknowledgement

The current research is supported by Hong Kong research grant council (16206617), National Science Foundation of China (11472219,11772281,91852114,11771035, U1530401).

References

  • (1) S. W. Davis, J. M. Stone, Y.-F. Jiang, A radiation transfer solver for athena using short characteristics, The Astrophysical Journal Supplement Series 199 (1) (2012) 9.
  • (2) A. Marshak, A. Davis, 3D radiative transfer in cloudy atmospheres, Springer Science & Business Media, 2005.
  • (3) A. D. Klose, U. Netz, J. Beuthan, A. H. Hielscher, Optical tomography using the time-independent equation of radiative transfer—part 1: forward model, Journal of Quantitative Spectroscopy and Radiative Transfer 72 (5) (2002) 691–713.
  • (4) B. Hunter, Z. Guo, Comparison of quadrature schemes in dom for anisotropic scattering radiative transfer analysis, Numerical Heat Transfer, Part B: Fundamentals 63 (6) (2013) 485–507.
  • (5) P. J. Coelho, Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media, Journal of Quantitative Spectroscopy and Radiative Transfer 145 (2014) 121–146.
  • (6) S.-S. Chen, B.-W. Li, Y.-S. Sun, Chebyshev collocation spectral method for solving radiative transfer with the modified discrete ordinates formulations, International Journal of Heat and Mass Transfer 88 (2015) 388–397.
  • (7) T. Roos, T. Harms, C. du Toit, Conservation of scattered energy and asymmetry factor in the new rotationally symmetric spherical discretisation scheme, International Journal of Heat and Mass Transfer 101 (2016) 205–225.
  • (8) M. Frank, B. Dubroca, A. Klar, Partial moment entropy approximation to radiative heat transfer, Journal of Computational Physics 218 (1) (2006) 1–18.
  • (9) J. A. Carrillo, T. Goudon, P. Lafitte, F. Vecil, Numerical schemes of diffusion asymptotics and moment closures for kinetic equations, Journal of Scientific Computing 36 (1) (2008) 113–149.
  • (10) V. Vikas, C. Hauck, Z. Wang, R. O. Fox, Radiation transport modeling using extended quadrature method of moments, Journal of Computational Physics 246 (2013) 221–241.
  • (11) G. W. Alldredge, R. Li, W. Li, Approximating the method by the extended quadrature method of moments for radiative transfer in slab geometry, Kinetic & Related Models 9 (2).
  • (12) J. Fleck, J. Cummings, An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport, J. Comput. Phys. 8 (3) (1971) 313–342.
  • (13) L. Lucy, Computing radiative equilibria with Monte Carlo techniques, Astronomy and Astrophysics 344 (1999) 282–288.
  • (14) C. K. Hayakawa, J. Spanier, V. Venugopalan, Coupled forward-adjoint Monte Carlo simulations of radiative transport for the study of optical probe design in heterogeneous tissues, SIAM Journal on Applied Mathematics 68 (1) (2007) 253–270.
  • (15) J. Fleck Jr, E. Canfield, A random walk procedure for improving the computational efficiency of the implicit Monte Carlo method for nonlinear radiation transport, Journal of Computational Physics 54 (3) (1984) 508–523.
  • (16) J. Giorla, R. Sentis, A random walk method for solving radiative transfer equations, Journal of Computational Physics 70 (1) (1987) 145–165.
  • (17) J. D. Densmore, T. J. Urbatsch, T. M. Evans, M. W. Buksas, A hybrid transport-diffusion method for Monte Carlo radiative-transfer simulations, Journal of Computational Physics 222 (2) (2007) 485–503.
  • (18) J. D. Densmore, K. G. Thompson, T. J. Urbatsch, A hybrid transport-diffusion Monte Carlo method for frequency-dependent radiative-transfer simulations, Journal of Computational Physics 231 (20) (2012) 6924–6934.
  • (19) A. Klar, An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit, SIAM journal on numerical analysis 35 (3) (1998) 1073–1094.
  • (20) G. Naldi, L. Pareschi, Numerical schemes for kinetic equations in diffusive regimes, Applied mathematics letters 11 (2) (1998) 29–35.
  • (21) S. Jin, Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations, SIAM Journal on Scientific Computing 21 (2) (1999) 441–454.
  • (22) S. Jin, L. Pareschi, G. Toscani, Uniformly accurate diffusive relaxation schemes for multiscale transport equations, SIAM Journal on Numerical Analysis 38 (3) (2000) 913–936.
  • (23) L. Mieussens, On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models, Journal of Computational Physics 253 (2013) 138–156.
  • (24) W. Sun, S. Jiang, K. Xu, An asymptotic preserving unified gas kinetic scheme for gray radiative transfer equations, Journal of Computational Physics 285 (2015) 265–279.
  • (25) W. Sun, S. Jiang, K. Xu, S. Li, An asymptotic preserving unified gas kinetic scheme for frequency-dependent radiative transfer equations, Journal of Computational Physics 302 (2015) 222–238.
  • (26) W. Sun, S. Jiang, K. Xu, A multidimensional unified gas-kinetic scheme for radiative transfer equations on unstructured mesh, Journal of Computational Physics 351 (2017) 455–472.
  • (27) W. Sun, S. Jiang, K. Xu, An implicit unified gas kinetic scheme for radiative transfer with equilibrium and non-equilibrium diffusive limits, Communications in Computational Physics 22 (4) (2017) 889–912.
  • (28) W. Sun, S. Jiang, K. Xu, An asymptotic preserving implicit unified gas kinetic scheme for frequency-dependent radiative transfer equations, International Journal of Numerical Analysis & Modeling 15.
  • (29) K. Xu, J.-C. Huang, A unified gas-kinetic scheme for continuum and rarefied flows, Journal of Computational Physics 229 (20) (2010) 7747–7764.
  • (30) K. Xu, Direct modeling for computational fluid dynamics: construction and application of unified gas-kinetic schemes, World Scientific, 2015.
  • (31) D. Jiang, M. Mao, J. Li, X. Deng, An implicit parallel UGKS solver for flows covering various regimes, Advances in Aerodynamics, (2019) 1:8 https://doi.org/10.1186/s42774-019-0008-5.
  • (32) C. Liu, K. Xu, A unified gas kinetic scheme for continuum and rarefied flows v: multiscale and multi-component plasma transport, Communications in Computational Physics 22 (5) (2017) 1175–1223.
  • (33) C. Liu, Z. Wang, K. Xu, A unified gas-kinetic scheme for continuum and rarefied flows vi: Dilute disperse gas-particle multiphase system, Journal of Computational Physics 386 (2019) 264 – 295.
  • (34) K. Xu, C. Liu, A paradigm for modeling and computation of gas dynamics, Physics of Fluids 29 (026101).
  • (35) C. Liu, Y. Zhu, K. Xu, Unified gas-kinetic wave-particle methods i: continuum and rarefied gas flow, arXiv preprint arXiv:1811.07141.
  • (36) Y. Zhu, C. Liu, C. Zhong, K. Xu, Unified gas-kinetic wave-particle methods ii: multiscale simulation on unstructured mesh, arXiv preprint arXiv:1903.11861.
  • (37) S. Jin, L. Pareschi, G. Toscani, Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations, SIAM Journal on Numerical Analysis 35 (6) (1998) 2405–2439.
  • (38) T. A. Brunner, Forms of approximate radiation transport, SAND2002-1778, Sandia National Laboratory.
  • (39) T. Camminady, M. Frank, K. Küpper, J. Kusch, Ray effect mitigation for the discrete ordinates method through quadrature rotation, Journal of Computational Physics.
  • (40) R. G. McClarren, C. D. Hauck, Robust and accurate filtered spherical harmonics expansions for radiative transfer, Journal of Computational Physics 229 (16) (2010) 5597–5614.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
349776
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description