on Time domain RTE

Overview

The time-dependent radiative transport equation (RTE) coupled with the material temperature equation describes the evolution of radiation and heat in a participating medium. We derive the decoupled integral formulation for the scalar flux and discuss the grey approximation used in numerical simulations.

🏷️ Mathematical Formulation

The coupled system for radiation intensity $u(x,t,\hat{s},\nu )$ and material temperature $T(x,t)$ is given by:

$$\begin{array}{rl}& \frac{1}{c}\frac{\partial u}{\partial t}+\hat{s}\cdot \nabla u+\sigma u=\sigma B(\nu ,T)\\ & C_v\frac{\partial T}{\partial t}={\int }_0^{\infty }{\int }_{S^{d-1}}\sigma (u-B(\nu ,T))d\hat{s}d\nu \end{array}$$

where $\sigma =\sigma (x,\nu ,T)$ is the absorption coefficient, $c$ is the speed of light, and $C_v$ is the material heat capacity. $B(\nu ,T)$ is the Planck function:

$$B(\nu ,T)=\frac{2h{\nu }^3}{c^2}\frac{1}{\exp (h\nu /kT)-1}$$

🏷️ Grey Approximation

In many applications, we assume the material properties are independent of frequency $\nu$ (the grey medium approximation). Integrating the Planck function over all frequencies yields:

$${\int }_0^{\infty }B(\nu ,T)d\nu =\frac{ac}{4\pi }{T}^4$$

where $a$ is the radiation constant. Defining the scalar flux $\varphi (x,t)={\int }_{S^{d-1}}ud\hat{s}$ and assuming isotropic emission, the equations simplify to:

$$\begin{array}{rl}& \frac{1}{c}\frac{\partial u}{\partial t}+\hat{s}\cdot \nabla u+\sigma u=\frac{\sigma ac{T}^4}{|S^{d-1}|}\\ & C_v\frac{\partial T}{\partial t}=\sigma (\varphi -ac{T}^4)\end{array}$$

🏷️ Optical Depth and Heterogeneous Media

In realistic media, the absorption coefficient $\sigma$ is a function of position, frequency, and local temperature, $\sigma =\sigma (x,\nu ,T)$. This introduces strong nonlinearities into the coupled system. For the grey approximation, we define the optical depth $\tau$ along a ray between two points $z_1$ and $z_2$:

$$\tau (z_1,z_2)={\int }_0^{|z_1-z_2|}\sigma (z_1-\mu \hat{s})d\mu$$

The attenuation of radiation intensity $u$ is then governed by $\exp (-\tau )$, representing the probability of a photon traveling between these points without being absorbed.

Varying Opacity

When $\sigma$ varies spatially, the integral formulation for the scalar flux $\varphi$ becomes:

$$\varphi (x,t)={\int }_{\Omega }\frac{\exp (-\tau (x,x^{\prime }))}{|S^{d-1}|r^{d-1}}\left(\sigma (x^{\prime })ac{T}^4(x^{\prime },t-r/c)\right)d{x}^{\prime }$$

For numerical efficiency in the treecode, if $\sigma$ is relatively smooth, we can approximate the optical depth using the average opacity between the target and the source cluster: $\tau (x,x^{\prime })\approx \bar{\sigma }|x-x^{\prime }|$. Smoothness is critical here, as the treecode’s far-field approximation assumes the kernel (which depends on $\sigma$) varies slowly across the source clusters. High-frequency oscillations or sharp discontinuities in $\sigma$ would require much smaller clusters or higher-order multipole expansions to maintain accuracy.

🏷️ Integral Formulation

By treating the transport equation as a first-order linear PDE along the characteristic rays $x(p)=x_0+p\hat{s}$ and $t(p)=t_0+p/c$, we can write the solution using the method of characteristics. Let $z=(x,t)$ and $\tilde{s}=(\hat{s},1/c)$, the transport equation becomes:

$$\tilde{s}\cdot {\nabla }_{z}u+\sigma u=H(z)$$

where $H=\sigma ac{T}^4/|S^{d-1}|$. Integrating along the ray from the boundary or from $t=-\infty$:

$$u(z,\hat{s})={\int }_0^{\infty }\exp \left(-{\int }_0^p\sigma (z-\mu \tilde{s})d\mu \right)H(z-p\tilde{s})dp$$

Integrating over all directions $\hat{s}$ gives the scalar flux $\varphi (z)$:

$$\varphi (x,t)={\int }_{S^{d-1}}{\int }_0^{\infty }\exp (-\tau (x,x-p\hat{s}))H(x-p\hat{s},t-p/c)dpd\hat{s}$$

Changing variables $x^{\prime }=x-p\hat{s}$ (where $r=|x-x^{\prime }|=p$ and $dpd\hat{s}=r^{-(d-1)}d{x}^{\prime }$), we obtain the general volume integral:

$$\varphi (x,t)={\int }_{\Omega }\frac{\exp (-\tau (x,x^{\prime }))}{|S^{d-1}|r^{d-1}}\left(\frac{\sigma (x^{\prime })ac{T}^4(x^{\prime },t-r/c)}{|S^{d-1}|}\right)d{x}^{\prime }$$

This is an integral with a retarded Green’s function, where the attenuation factor is determined by the accumulated optical depth between the source and the target.

🏷️ Integro-differential Equation for Material Temperature

By substituting the integral expression for $\varphi (x,t)$ into the material temperature equation, we obtain a closed-form nonlinear integro-differential equation. It is often computationally convenient to solve for the radiation energy density $E(x,t)=ac{T}^4(x,t)$. Using the chain rule, the evolution of $E$ is governed by:

$$\frac{\partial E}{\partial t}=4ac{T}^3\frac{\partial T}{\partial t}=\frac{4\sigma (ac{)}^{1/4}{E}^{3/4}}{C_v}\left(\varphi [E](x,t)-E(x,t)\right)$$

where $\varphi [E]$ represents the retarded integral operator:

$$\varphi [E](x,t)={\int }_{\Omega }\frac{\exp (-\tau (x,x^{\prime }))}{|S^{d-1}|r^{d-1}}\sigma (x^{\prime })E(x^{\prime },t-|x-x^{\prime }|/c)d{x}^{\prime }$$

This formulation highlights the non-local and retarded nature of the coupling. In the limit of high opacity ($\sigma \to \infty$), the system relaxes towards local thermodynamic equilibrium where $\varphi \approx E$.

And solving this can apply some fast algorithm like treecode or FMM as we did before.

🏷️ Numerical Algorithm

The solver utilizes a time-marching scheme that alternates between evaluating the non-local scalar flux and updating the local energy density through a semi-implicit relaxation.

🏁 1. Initialization

Initialize the material temperature $T^0(x)$ and compute the equilibrium radiation energy density $E^0(x)=ac(T^0{)}^4$. The initial emission source is set to $S^0(x)=\sigma (x)E^0(x)$.

⏳ 2. Time-stepping Loop

For each time step $t_n\to t_{n+1}$, the following operations are performed:

A. Retarded Potential Integration Evaluate the scalar flux ${\varphi }^n(x)$ via the retarded volume integral over the domain $\Omega$:

$${\varphi }^n(x)=\sum _j{w}_j\frac{\exp (-\tau (x,x_j))}{|S^{d-1}|r_j^{d-1}}S(x_j,t_n-r_j/c)$$

• Treecode: Clusters distant sources to accelerate the summation from $O(N^2)$ to $O(N\log N)$.
• Composite Simpson: Computes the accumulated optical depth $\tau (x,x_j)$ along the ray path using a multi-interval quadrature of $\sigma$.
• History Interpolation: Retrieves the source term $S$ at retarded times using a temporal history buffer and linear interpolation.

B. Semi-implicit Energy Relaxation Evolve the radiation energy $E^{n+1}$ by solving the stiff material coupling equation. We employ a semi-implicit Euler scheme for numerical stability as $\sigma$ increases:

$$E^{n+1}=\frac{E^n+\Delta t{\kappa }^n{\varphi }^n}{1+\Delta t{\kappa }^n},\text{where }{\kappa }^n=\frac{4\sigma (ac{)}^{1/4}(E^n{)}^{3/4}}{C_v}$$

C. Source Synchronization Update the isotropic emission source $S^{n+1}=\sigma E^{n+1}$ and append the new values to the history buffer for use in future time steps.

🏷️ 2D Implementation and Treecode Acceleration

In the 2D spatial domain, the scalar flux $\varphi (x,t)={\int }_{S^1}u(x,t,\hat{s})d\hat{s}$ satisfies an integral equation involving a retarded Green’s function. By approximating the kernel, we can express $\varphi$ as:

$$\varphi (x,t)={\int }_{\Omega }\frac{\exp (-\sigma r/c)}{2\pi r}S(x^{\prime },t-r/c)d{x}^{\prime }$$

where $r=|x-x^{\prime }|$ and $S$ is the source term proportional to $\sigma ac{T}^4$. This integral can be accelerated using a Treecode algorithm. The core idea is to cluster distant sources and represent their contribution using a single far-field approximation at the cluster’s center of mass.

Treecode for Retarded Potentials

For time-dependent problems, the treecode must account for the retarded time $t_{ret}=t-r/c$. When evaluating the contribution of a distant node $N$ to a target point $x$, we use:

$${\varphi }_N(x,t)\approx \frac{\exp (-\sigma R/c)}{2\pi R}\sum _{j\in N}{w}_{j}S(x_j,t-R/c)$$

where $R$ is the distance from $x$ to the center of node $N$. This reduces the complexity from $O(N^2)$ to $O(N\log N)$.

📊 Numerical Verification

We implemented the time-domain RTE solver in 2D using an optimized quadtree-based treecode. The system was discretized on an 80x80 grid with a refined time step $\Delta t=0.05$ to capture the fast transient dynamics of radiation transport. The simulation features a complex initial thermal profile (multiple Gaussian peaks and a periodic background) and smooth heterogeneous fields for both opacity $\sigma (x,y)$ and heat capacity $C_v(x,y)$.

The animation illustrates the high-resolution dissipation and expansion of the thermal energy through radiation over time, influenced by the spatially varying material properties. The treecode maintains $O(N\log N)$ complexity, allowing for efficient simulation of large-scale problems while ensuring numerical stability via semi-implicit energy relaxation.

• Source Code: on_Time_domain_RTE.py in the /codes/2016 Fall/ directory.
• Results: The final temperature distribution shows the expected smoothing and expansion of the thermal front.

Links

• on Adjoint Framework --- The adjoint method can be applied to this transport system for parameter estimation of $\sigma$.
• on Transport Numerical --- Early notes on transport equation discretization.