Version 164 (modified by 12 years ago) ( diff ) | ,
---|
- Physics of Radiation Transfer
- Numerics of Flux Limited Diffusion
Physics of Radiation Transfer
Typically when we discuss the radiation field we use the spectral intensity
which is a function of frequency, position, and direction. This is very similar to the phase space density used in deriving the fluid equations except that- light always travels at 'c', so the velocity dependence is just a direction dependence.
- Furthermore, photons can have different frequencies, so there is an extra dimension to the phase space.
- Instead of storing the phase space density of photons, the spectral intensity is the phase space density of energy flux…
Going between photon number and energy just involves a factor of
and going from energy density to energy flux density just involves a factor of so we have
This can also be seen by considering the differential energy
where the number of photons traveling normal to the surface dA that cross the surface dA in time dt is just the number of photons in the volume dV = dA c dt (assuming the photons are headed normal to dA)…
so we also have
which gives
There are also parallels between the radiation transport equation and the Boltzmann transport equation except that photons don't experience body forces, and photon's don't collide with each other, but scatter off of particles.
If we solve the transport equation along a characteristic
we have
where
and then we can divide through by
we get
Now if we define
which gives
we can write the transport equation in the simplest form
although the RHS is now more difficult to evaluate as
Also if we include scattering then the source function can depend on the mean radiative flux and the transport equation becomes an integro-differential equation that must be solved iteratively…
There are also a few important dimensionless numbers to consider:
streaming limit | |
static diffusion limit | |
dynamic diffusion limit |
Equations of Radiation Hydrodynamics
Equations of Radiation Hydrodynamics
Some of what follows is taken from Krumholz et al. 2007
where the moments of the specific intensity are defined as
Radiation Energy moments | Corresponding fluid moments |
and the radiation 4-force density is given by
If we had a closure relation for the radiation pressure then we could solve this system. For gas particles, collisions tend to produce a Boltzmann Distribution which is isotropic and gives a pressure tensor that is a multiple of the identity tensor. Photons do not "collide" with each other and they all have the same velocity 'c' but in various directions. If the field were isotropic than
but in general where 'f' is the Eddington Tensor.
Simplifiying assumptions
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
where
comoving-frame Planck function weighted opacity | |
comoving-frame radiation energy weighted opacity | |
comoving-frame radiation flux weighted opacity |
- If the radiation has a blackbody spectrum then
- If the radiation is optically thick, then
which implies that
- In the optically thin regime, so we would have however assuming a blackbody temperature in the optically thin limit may be any more accurate than assuming that
Flux limited diffusion
Flux limited diffusion
The flux limited diffusion approximation drops the radiation momentum equation in favor of
where
is the flux-limiterwhich corresponds to a pressure tensor
If we Lorentz boost the comoving terms into the lab frame and keep terms necessary to maintain accuracy we get:
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:
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), and a coupled implicit solve (red) for the radiation energy density and thermal energy density (ie temperature).
Operator Splitting 1
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.
Explcit Update 1
Eplicit Update 1
The extra terms in the explicit update due to radiation energy are as follows:
These can be discretized as follows:
where
where
and
and
and
and
and
Implicit Update 1
Implicit Update 1
For now we will assume that
and 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:where
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
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
where
Then the system of equations becomes
which will be accurate as long as
orWe can calculate the time scale for this to be true using the evolution equation for the energy density
which gives
Implicit Discretization 1
Which we can discretize for (1D) as
where the diffusion coefficient is given by
and where
and
represents the number of absorption/emissions during the time step
and
and we can think of the radiative flux as
Time Discretization
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
and for Crank Nicholson we have or we can parameterize the solution whereBackward Euler has
and Crank Nicholson hasForward Euler has
In any event in 1D we have the following matrix coefficients
Now since the second equation has no spatial dependence, we can solve it for
and plug the result into the first equation to get a matrix equation involving only one variable.
If we ignore the change in the Planck function due to heating during the implicit solve, it is equivalent to setting
This gives the following equations: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…).
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 (Free streaming) boundaries
We would like the radiation to leave at the free streaming limit. So
Clearly if we set
and we should get the correct flux.This corresponds to an
So we would just modify
and zero out the matrix coefficient to the ghost zoneConstant Slope Boundary
Here we want the flux to be constant so energy does not pile up near the boundary. If we cancel all derivative terms on both sides of the cell, this will effectively match the incoming flux with the outgoing flux. This can also be done by setting
Periodic Boundary
This is the same as internal zones - it just maps the neighbor cell to be across the domain. Hypre has built in functionality for this under for the Struct Interface
User defined radiation field/Coarse Fine boundary
This will be the boundary at internal coarse-fine boundaries, but could also be used at the physical boundary if the radiation energy were specified.
Reflecting/ZeroSlope Boundary
Reflecting boundary should be fairly straightforward. This an be achieved by setting
which zeros out any flux - and has the same effect as setting orConstant radiative flux
To have a constant radiative flux we must zero out terms involving the gradient and just add
in the source vectorSummary
Numerical value | Boundary | ||||
0 | RAD_FREE_STREAMING | ||||
1 | RAD_EXTRAPOLATE_FLUX | ||||
2 | INTERNAL/PERIODIC | ||||
3 | RAD_REFLECTING | ||||
4 | RAD_USERDEFINED_BOUNDARY/AMR_BOUNDARY | ||||
5 | RAD_USERDEFINED_FLUX |
Alternative Splitting Method
It is possible that with included the terms in orange in a semi-implicit method, the dynamic diffusion regime may be stable… In any event, it costs very little to add all of the terms in orange to the implicit solve (using the old velocity). Then the momentum update can be done explicitly - though using a time centered radiation field.
Operator Splitting 2
Operator splitting 2
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
- Store EDot in child arrays to be prolongated
- Do first EHmbc
- Do ERmbc —- Only momentum terms.
- Step 2
- Overlap d, p, e, E, and do physical BCs
- Do IR which updates e0, and E0 using d1, e1, E1
- Update Embc using Edotmbc
- Update e2*mbc using E2*mbc, Edot2*mbc, and e2*mbc
- Update Edot0 using pre IR and post IR E0
- Ghost embc, E1, Edot1
- Store EDot in child arrays to be prolongated
- Do second EH0
- Do ER0 —- Only momentum terms.
Explcit Update 2
Eplicit Update 2
The extra terms in the explicit update due to radiation energy are as follows:
These can be discretized as follows:
and
Implicit Update 2
Implicit Update 2
For now we will assume that
and 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:which we can also write as
where
and
Now we can linearize f about e0
so that the first equation can be written as
and then discretized as
where
which can be solved for
Then if we take the semi-discretized equation for E
and then plugin the solution for en+1i we get
which simplifies to
Now we have 1 equation with 1 variable that we can solve implicitly using hypre, and then we can use En+1 and En to construct E* which we can plug into the equation for en+1
Expanding f about e0
Expanding
where
and we can identify
Then the equation for e becomes
which will be accurate as long as
orWe can calculate the time scale for this to be true using the evolution equation for the energy density
which gives
Implicit Discretization 2
Now we can discretize
as
which along with the other terms gives
where the diffusion coefficient is given by
and where
and
and
and
and
and
and
where
and
and
and
and
and
and
Which we can arrange into the following form
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 terms for each dimension - and the velocity components (vx) will change, 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
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
For each boundary type we need to specify
in terms of other quantities (ie including En+1_i). The vg and kappag terms should come from the hydro boundary conditions so we just need an equation for EngOpen (Free streaming) boundaries
We would like the radiation to leave at the free streaming limit. But we could have to modify the opacity in the ghost zone to keep the radiation energy positive
Clearly if we set
and we should get and a free streaming flux ofwhich also corresponds to an
Constant Slope Boundary
Periodic Boundary
This is the same as internal zones - it just maps the neighbor cell to be across the domain. Hypre has built in functionality for this under for the Struct Interface
User defined radiation field/Coarse Fine boundary
This will be the boundary at internal coarse-fine boundaries, but could also be used at the physical boundary if the radiation energy were specified.
Reflecting/ZeroSlope Boundary
Constant radiative flux
Attachments (1)
- fld.png (28.2 KB ) - added by 12 years ago.
Download all attachments as: .zip