Physics of Radiation Transfer
tau = l \kappa=\frac{l}{\lambda_p
|
beta = \frac{u}{c
|
tau << 1 | streaming limit
|
tau >> 1 \mbox{, } \beta \tau << | static diffusion limit
|
tau >> 1 \mbox{, } \beta \tau >> | dynamic diffusion limit
|
Equations of Radiation Hydrodynamics
frac{\partial \rho}{\partial t} + \nabla \cdot \left ( \rho \mathbf{v} \right ) =
|
frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) + \nabla \cdot \left ( \rho \mathbf{vv} \right ) = -\nabla P+\mathbf{G
|
frac{\partial e}{\partial t} + \nabla \cdot \left [ \left ( e + P \right ) \mathbf{v} \right ] = c G^
|
frac{\partial E}{\partial t} + \nabla \cdot \mathbf{F} = -c G_0
|
frac{1}{c^2} \frac{\partial \mathbf{F}}{\partial t} + \nabla \cdot \mathbf{P} = -\mathbf{G
|
where
E=\int_0^\infty d \nu \int d \Omega I(\mathbf{n}, \nu
|
mathbf{F}=\int_0^\infty d \nu \int d \Omega \mathbf{n} I(\mathbf{n}, \nu
|
\mathbf{P}=\int_0^\infty d \nu \int d \Omega \mathbf{nn} I(\mathbf{n}, \nu
|
and
G^0 = \int_0^\infty d \nu \int d \Omega \left [ \kappa (\mathbf{n}, v)I(\mathbf{n}, \nu)-\eta(\mathbf{n},v) \right
|
\mathbf{G} = \int_0^\infty d \nu \int d \Omega \left [ \kappa (\mathbf{n}, v)I(\mathbf{n}, \nu)-\eta(\mathbf{n},v) \right ] \mathbf{n
|
Simplifying assumptions
- If the flux spectrum of the radiation is direction-independent then we can write the radiation four-force density in terms of the moments of the radiation field
^0=\gamma [ \gamma^2 \kappa_{0E} + (1-\gamma^2)\kappa_{0F}]E-\gamma \kappa_{0P}a_RT_0^4-\gamma(\mathbf{v} \cdot \mathbf{F}/c^2)[\kappa_{0F}-2\gamma^2(\kappa_{0F}-\kappa_{0E})] - \gamma^3(\kappa_{0F}-\kappa+{0E})(\mathbf{vv}) : \mathbf{P}/c^
|
mathbf{G} = \gamma \kappa_{0F}(\mathbf{F}/c)-\gamma \kappa_{0P}a_RT_0^4 (\mathbf{v}/c) - \left [ \gamma^3( \kappa_{0F} - \kappa_{0E})( \mathbf{v}/c)E + \gamma \kappa_{0F}(\mathbf{v}/c) \cdot \mathbf{P} \right ] + \gamma^3 ( \kappa_{0F} - \kappa_{0E}) \left [ 2 \mathbf{v} \cdot \mathbf{F}/c^3 - ( \mathbf{vv}) : \mathbf{P}/c^3 \right ] \mathbf{v}
|
where
kappa_{0P | comoving-frame Planck function weighted opacity
|
kappa_{0E | comoving-frame radiation energy weighted opacity
|
kappa_{0F | comoving-frame radiation flux weighted opacity
|
- If the radiation has a blackbody spectrum then kappa_{0E}=\kappa_{0P
- If the radiation is optically thick, then
mathbf{F_0}(\nu_0) \propto -\nabla E_0(\nu_0)/\kappa_0(\nu_0) \propto -[\partial B(\nu_0, T_0) / \partial T_0](\nabla T_0)/\kappa_0(\nu_0
which implies that kappa_{0F}^{-1}=\kappa_{0R}^{-1}=\frac{\int_0^\infty d \nu_0 \kappa_0(\nu_0)^{-1}[\partial B(\nu_0,T_0)/\partial T_0]}{\int_0^\infty d \nu_0 [\partial B(\nu_0, T_0)/\partial T_0]
- In the optically thin regime, \mathbf{F}_0(\nu_0)| \rightarrow cE_0(\nu_0 so we would have kappa_{0F}=\kappa_{0E however assuming a blackbody temperature in the optically thin limit may be any more accurate than assuming that kappa_{0F}=\kappa_{0R
Flux limited diffusion
The flux limited diffusion approximation drops the radiation momentum equation in favor of
mathbf{F}_0=-\frac{c\lambda}{\kappa_{0R}}\nabla E_
where lambd is the flux-limiter
lambda = \frac{1}{R} \left ( \coth R - \frac{1}{R} \right )
=\frac{|\nabla E_0|}{\kappa_{0R}E_0
which corresponds to a pressure tensor
mathbf{P}_0=\frac{E_0}{2}[(1-R_2)\mathbf{I}+(3R_2-1)\mathbf{n}_0\mathbf{n}_0
_2=\lambda+\lambda^2R^
If we Lorentz boost the comoving terms into the lab frame and keep terms necessary to maintain accuracy we get:
^0=\kappa_{0P} \left ( E-\frac{4 \pi B}{c} \right ) + \left ( \frac{\lambda}{c} \right ) \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}} - 1 \right ) \mathbf{v} \cdot \nabla E - \frac{\kappa_{0P}}{c^2} E \left [ \frac{3-R_2}{2}v^2 + \frac{3R_2-1}{2}(\mathbf{v} \cdot \mathbf{n})^2 \right ] + \frac{1}{2} \left ( \frac{v}{c} \right ) ^2 \kappa_{0P} \left ( E - \frac{4 \pi B}{c} \right )
mathbf{G} = -\lambda \nabla E + \kappa_{0P} \frac{\mathbf{v}}{c} \left ( E - \frac{4 \pi B}{c} \right ) - \frac{1}{2} \left ( \frac{v}{c} \right ) ^2 \lambda \nabla E + 2 \lambda \left ( \frac{\kappa_{0P}}{\kappa_{0R}} - 1 \right ) \frac{(\mathbf{v} \cdot \nabla E )\mathbf{v}}{c^2
Numerics of Flux Limited Diffusion
If we plug the expressions for the radiation 4-momentum back into the gas equations and keep terms necessary to maintain accuracy we get:
frac{\partial }{\partial t} \left ( \rho \mathbf{v} \right ) + \nabla \cdot \left ( \rho \mathbf{vv} \right ) = -\nabla P\color{green}{-\lambda \nabla E
frac{\partial e}{\partial t} + \nabla \cdot \left [ \left ( e + P \right ) \mathbf{v} \right ] = \color{red}{-\kappa_{0P}(4 \pi B-cE)} \color{green}{+\lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1 \right ) \mathbf{v} \cdot \nabla E} \color{blue}{-\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E
frac{\partial E}{\partial t} \color{red}{ - \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E} = \color{red}{\kappa_{0P} (4 \pi B-cE)} \color{green}{-\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1 \right )\mathbf{v}\cdot \nabla E} \color{green}{-\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E \right )} \color{blue}{+\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E}
For static diffusion, the terms in blue with v2/c can be dropped and the system can be split into the usual hydro update (black), radiative source terms (green) using time centered radiation energy, and a coupled implicit solve (red) for the radiation energy density and thermal energy density (ie temperature). If the opacity is independent of temperature and radiation energy density, then the implicit solve only involves the radiation energy density. Otherwise some sort of sub-cycling would be required.
Operator splitting
Krumholz et al. perform Implicit Radiative, Explicit Hydro, Explicit Radiative
In AstroBEAR this would look like:
- Initialization
- Prolongate, d, p, e, E, Edot
- Step 1
- Overlap d, p, e, E and do physical BC's
- Do IR which updates e0, and E0 using d1, e1, E1
- Update E2*mbc using Edot2*mbc
- Update e2*mbc using E2*mbc, Edot2*mbc, and e2*mbc
- Update Edot0 using pre IR and post IR E0
- Ghost e2*mbc, Embc+1, Edotmbc+1
- Do first EHmbc
- Do ERmbc —- Terms with grad E can be done without ghosting since EH did not change E. The del dot vE term needs time centered face centered velocities which can be stored during the hydro update.
- Store Edot in child arrays to be prolongated
- Step 2
- Overlap d, p, e, E, and do physical BCs
- Do IR which updates e0, and E0 using d1, e1, E1
- Update Edot0 using pre IR and post IR E0
- Update E1 using Edot1
- Ghost embc, E1, Edot1
- Do second EH0
- Do ER0 —- Terms with grad E can be done without ghosting since EH did not change E. The del dot vE term needs time centered face centered velocities which can be stored during the hydro update.
Implicit Discretization
For now we will assume that kappa_{0P and kappa_{0R are constant over the implicit update and we will treat the energy as the total internal energy ignoring kinetic and magnetic contributions. In this case we can solve the radiation energy equations:
frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} (4 \pi B(T)-cE) = \nabla \cdot \mathbf{F} + \kappa_{0P} (4 \pi B(T)-cE
frac{\partial e}{\partial t} = - \kappa_{0P} (4 \pi B(T)-cE
where mathbf{F} = \frac{c\lambda}{\kappa_{0R}} \nabla
Expanding about e0
Of course even if the opacity is independent of energy and radiation energy, the above combined system of equations is still non-linear due to the dependence on Temperature of the Planck Function (T.
If we ignore the changes in the Temperature due to heating during the implicit step which would feed back into the source function. We can improve this by writing
(T) = B \left ( T_0+dT \right ) = B \left ( T_0 \right ) + \left . \frac{\partial B}{\partial T} \right | _{T_0} \frac{\partial T}{\partial e} de = \frac{c}{4 \pi} a_R \left ( T_0^4 + 4T_0^3\Gamma de \right ) = B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )
where
Gamma = \frac{\partial T}{\partial e} = \frac{(\gamma-1)}{n k_B
Then the system of equations becomes
frac{\partial E}{\partial t} = \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E + \kappa_{0P} \left [4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ]
frac{\partial e}{\partial t} = - \kappa_{0P} \left [ 4 \pi B_0 \left ( 1 + 4\Gamma \frac{e-e_0}{T_0} \right )-cE \right ]
1D
Which we can discretize for (1D) as
^{n+1}_i-E^{n}_i = \left [ \alpha^n_{i+1/2} \left ( E^{*}_{i+1}-E^{*}_{i} \right ) - \alpha^n_{i-1/2} \left ( E^{*}_{i}-E^{*}_{i-1} \right ) \right ] - \epsilon^n_i E^{*}_i + \phi^n_i e^{*}_i + \theta^n_i)
^{n+1}_i-e^{n}_i = \epsilon^n_i E^{*}_i - \phi^n_i e^{*}_i - \theta^n_i
where the diffusion coefficient is given by
alpha_{i+1/2}=\frac{\Delta t}{\Delta x^2} \frac{c \lambda_{i+1/2}}{\kappa_{i+1/2}} \mbox{ where } \kappa_{i+1/2} = \frac{\kappa^n_{i}+\kappa^n_{i+1}}{2} \mbox{ and } \lambda_{i+1/2} = \frac{1}{R_{i+1/2}} \left ( \coth R_{i+1/2} - \frac{1}{R_{i+1/2}} \right )
and where
_{i+1/2} = \frac{\left | E^n_{i+1}-E^n_{i} \right | }{2 \kappa_{i+1/2} \left ( E^n_i+E^n_{i+1} \right )
and
epsilon^n_i=c\Delta t \kappa^n_{0P,i
represents the number of absorption/emissions during the time step
and
phi = \epsilon^n_i \frac{4 \pi}{c} B \left ( T^n_i \right ) \left ( \frac{4\Gamma}{T^n_i} \right )
theta = \epsilon^n_i \frac{4 \pi}{c} B \left ( T^n_i \right ) \left ( 1 - 4\Gamma \frac{e^n_i}{T^n_i} \right )
and we can think of the radiative flux as
frac{\Delta t}{\Delta x}\mathbf{F}^n_{i+1/2} = \alpha^n_{i+1/2} \left ( E^{*}_{i+1} - E^{*}_i \right )
Now all the terms on the right hand side that are linear in E or e have been written as E* or e* because there are different ways to approximate E* (e*). For Backward Euler we have
^*_i = E^{n+1}_
and for Crank Nicholson we have
^*_i = \frac{1}{2} \left ( E^{n+1}_i + E^n_i \right )
or we can parameterize the solution
^*_i = \psi E^{n+1}_i + \bar{\psi}E^n_
where bar{\psi} = 1-\ps
Backward Euler has psi= and Crank Nicholson has psi=1/
Forward Euler has psi=
In any event in 1D we have the following matrix coefficients
left [ 1 + \psi \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^{n+1}_i - \left ( \psi \alpha^n_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha^n_{i-1/2} \right ) E^{n+1}_{i-1} - \left ( \psi \phi^n_i \right ) e^{n+1}_i=\left [ 1 - \bar{\psi} \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^n_i+\bar{\psi}\phi^n_i e^n_i + \theta^n_
left ( 1 +\psi \phi^n_i \right ) e^{n+1}_i - \left ( \psi \epsilon^n_i \right )E^{n+1}_i =\left ( 1 - \bar{\psi}\phi^n_i \right ) e^n_i + \left ( \bar{\psi} \epsilon^n_i \right ) E^n_i-\theta^i_n
If we ignore the change in the Planck function due to heating during the implicit solve, it is equivalent to replacing the term with psi \phi^n_i e^{n+1}_ with psi \phi^n_i e^n_ in which case the equations simplify to
left [ 1 + \psi \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^{n+1}_i - \left ( \psi \alpha^n_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha^n_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ 1 - \bar{\psi} \left( \alpha^n_{i+1/2} + \alpha^n_{i-1/2} + \epsilon^n_i \right ) \right ] E^n_i+\phi^n_i e^n_i + \theta^n_
^{n+1}_i = e^n_i + \epsilon^n_i \left [ \left ( \psi E^{n+1}_i + \bar{\psi} E^{n}_i \right ) - \frac{4 \pi}{c} B \left ( T^n_i \right ) \right ]
In this case the first equation decouples and can be solved on it's own, and then the solution plugged back into the second equation to solve for the new energy.
2D etc…
For 2D or 3D we have more connections to add to the matrix elements but it is very straight forward… There will be additional alpha terms for each dimension, but everything else stays the same.
Initial solution vector
For the initial solution vector, we can just use Edot from the parent update (or last time step if we are on the coarse grid) to guess E, and then we can solve for the new e given our guess at the new E using the same time stepping (Backward Euler, Crank Nicholson, etc…).
^{n+1}_i = E^{n}_i+\dot{E}^{n}_i \Delta
^{n+1}_i = \frac{1}{1+\psi \phi^n_i} \left \{ \left ( \psi \epsilon^n_i \right )E^{n+1}_i + \left ( 1 - \bar{\psi}\phi^n_i \right ) e^n_i + \left ( \bar{\psi} \epsilon^n_i \right ) E^n_i-\theta^i_n \right \}
Coarse-Fine Boundaries
Since we are doing our implicit solves first, we can use time interpolated solutions for the implicit solve for non-refined ghost zones. To do this we just need Edot. The opacities etc… in the ghost zones can be obtained from the hydro terms.
And the radiative implicit heating in coarse ghost cells can be done along with the initial solution vector so they are available for the hydro update.
Physical Boundary Conditions
Open boundaries
We would like the radiation to leave at the free streaming limit. So
frac{c \lambda}{\kappa_{0R}} \nabla E = \mathbf{F} = cE\mathbf{n} = \frac{c \lambda_g}{\kappa_g}\frac{ \left ( E-E_g \right ) }{\Delta x
Clearly if we set _g = and lambda_g =\kappa_g \Delta we should get the correct flux.
This corresponds to an
alpha = c \frac{\Delta t}{\Delta x
So we would just modify alph and zero out the matrix coefficient to the ghost zone
Thermally emitting boundary
Another possible boundary condition would be to have the edge of the grid be adjacent to some thermal emitter. This corresponds to setting the radiation energy density in the ghost zones to the Planck function.
_g = \frac{4 \pi}{c} B(T_g however, we still need an opacity which could be defined from the fluid properties of the ghost region… This would essentially be a thermally emitting boundary with the temperature, density, opacity, etc… derived from the hydro boundary type. If the boundary type was extrapolating, and the density and temperature uniform, one could have a constant thermal energy spectrum… However we would need to either specify the opacity and temperature of the ghost zone - or the various fluid properties needed to reconstruct the opacity and temperature of the ghost zone - or make an extra call to set physicalBC - since this may happen in between a hydro step and a physical boundary update.
Reflecting Boundary
Reflecting boundary should be fairly straightforward. This an be achieved by setting alpha_g = which zeros out any flux - and has the same effect as setting ^{n+1}_g=E^{n+1}_
Constant radiative flux
To have a constant radiative flux we must have
alpha_g \left ( E^{n+1}_i-E^{n+1}_g \right ) = F_0 \frac{\Delta t}{\Delta x
Which we can solve for
^{n+1}_g = E^{n+1}_i - \frac{F_0 \Delta t}{\alpha_g \Delta x
but when we plug this into the coefficient matrix the terms with alpha_ cancel and we just get _0 \frac{\Delta t}{\Delta x in the source vector
Summary
Boundary | alph |
|
Open | \frac{\Delta t}{\Delta x | 0
|
User-Defined opacity and Temperature | alpha_ | alpha_g \frac{4 \pi}{c} B(T_g
|
Reflecting | |
|
User-Defined Flux | | _0 \frac{\Delta t}{\Delta x
|