Diffusion Solvers in AstroBEAR

Astrobear supports both implicit and explicit thermal diffusion that is both isotropic and anistropic. We first derive the equations for the implicit-anisotropic case, and then discuss simplifications that arise for isotropic and explicit versions.

Basic Equations

The anisotropic conduction equation is

\rho c_v\frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \left ( \chi_\parallel - \chi_\perp \right ) \left ( \hat{b} \cdot \nabla T \right ) + n \chi_\perp \nabla T \right ]

however, the conduction coefficients have a Temperature dependence.

\chi_\parallel = \kappa_\parallel T^{5/2} n^{-1}

\chi_\perp = \kappa_\perp \frac{n}{B^2 T^{1/2}}

Collecting power's of T

So we can rewrite the equations as

\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \left [ \hat{b} \left ( \kappa_\parallel T^{\lambda_\parallel} - \frac{n^2 \kappa_\perp}{B^2} T^{\lambda_\perp} \right ) \left ( \hat{b} \cdot \nabla T \right ) + \frac{n^2 \kappa_\perp}{B^2} T^{\lambda_\perp} \nabla T \right ]

or

\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \left [ \hat{b} \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( \hat{b} \cdot \nabla T^{\lambda_\parallel+1} \right ) - \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( \hat{b} \cdot \nabla T^{\lambda_\perp + 1} \right ) \right ) + \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \nabla T^{\lambda_\perp+1} \right ]

Einstein simplification

or in Einstein notation

\rho c_v \partial_t T = \partial_i \left [ b_i \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( b_j \partial_j T^{\lambda_\parallel+1} \right ) - \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( b_j \partial_j T^{\lambda_\perp + 1} \right ) \right ) + \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \partial_i T^{\lambda_\perp+1} \right ]

or

\rho c_v \partial_t T = \partial_i \left [ b_i \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( b_j \partial_j T^{\lambda_\parallel+1} \right ) - \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( b_j \partial_j T^{\lambda_\perp + 1} \right ) \right ) + \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp +1 \right )} \delta_{ij} \partial_j T^{\lambda_\perp+1} \right ]

or

\rho c_v \partial_t T = \partial_i b_i \frac{\kappa_\parallel}{\lambda_\parallel+1} b_j \partial_j T^{\lambda_\parallel+1} + \partial_i \left ( \delta_{ij} - b_i b_j \right ) \frac{n^2 \kappa_\perp}{B^2 \left ( \lambda_\perp + 1\right )} \partial_j T^{\lambda_\perp+1}

or

\rho c_v \partial_t T = \frac{\kappa_\parallel}{\lambda_\parallel+1} \partial_i b_i b_j \partial_j T^{\lambda_\parallel+1} + \frac{\kappa_\perp}{\left ( \lambda_\perp + 1\right )} \partial_i n^2 B^{-2} \left ( \delta_{ij} - b_i b_j \right ) \partial_j T^{\lambda_\perp +1}

Implicitization

Now we can rewrite the equation

\partial_t T = \displaystyle \sum_{\parallel, \perp} A \partial_i B_{ij} \partial_j T^{\lambda+1}

where the \parallel or \perp subscript on A, B_{ij}, and \lambda is implied.

A_\parallel = \frac{\kappa_\parallel}{\rho c_v \left ( \lambda_\parallel + 1\right ) } and A_\perp = \frac{\kappa_\perp}{\rho c_v \left (\lambda_\perp + 1 \right )}

and B_{\parallel, ij} = b_i b_j and B_{\perp, ij} = n^2 B^{-2} b_i b_j \left ( 1 - \delta_{ij} \right )

Now to solve this implicitly, we need to replace T with T_* where

T_* = T + \phi ( T' - T ) = (1-\phi) T + \phi T'

Note for Backward Euler, \phi = 1 and for Crank Nicholson, \phi = 1/2

\partial_t T = \displaystyle \sum_{\parallel, \perp} A \partial_i B_{ij} \partial_j T_*^{\lambda+1}

It is simpler to linearize the equation in terms of T_* and then subsitute then vice-versa - though both give the same answer

T_*^{\lambda + 1} \approx T^{\lambda + 1} + \left ( \lambda + 1 \right ) T^{\lambda} \left ( T_{*} - T \right ) = T^{\lambda + 1} + \left ( \lambda + 1 \right ) T^{\lambda} \phi \left ( T' - T \right ) = \left(1 - \phi \left ( \lambda + 1\right ) \right ) T^{\lambda + 1}+ \phi \left ( \lambda + 1 \right ) T^{\lambda} T'

So we have

\partial_t T = \displaystyle \sum_{\parallel, \perp}{A \partial_i B_{ij} \partial_j \left [ \left ( 1 - \phi \left ( \lambda + 1 \right ) \right )T^{\lambda + 1} + \phi \left ( \lambda + 1 \right ) T^{\lambda} T' \right ]}

\partial_t T = \displaystyle \sum_{\parallel, \perp}{A \partial_i B_{ij} \left [ \left ( 1 - \phi \left ( \lambda + 1 \right ) \right ) \partial_j T^{\lambda + 1} + \phi \left ( \lambda + 1 \right ) \partial_j T^{\lambda} T' \right ]}

\partial_t T = \displaystyle \sum_{\parallel, \perp}{A \partial_i B_{ij} \left [ C \partial_j T^{\lambda + 1} + D \partial_j T^{\lambda} T' \right ]}

where C= \left ( 1 - \phi \left ( \lambda + 1 \right ) \right ) and D = \phi \left ( \lambda + 1 \right )

Now we can expand the derivatives and get

Expand Derivatives

\partial_t T = \displaystyle \sum_{\parallel, \perp}{A \left ( \partial_i B_{ij} \right ) \left ( C \partial_j T^{\lambda + 1} + D \partial_j T^{\lambda} T' \right ) + A B_{ij} \left ( C \partial_i \partial_j T^{\lambda + 1} + D \partial_i \partial_j T^{\lambda} T' \right )}

We can also write this as

\partial_t T = \displaystyle \sum_{\parallel,\perp}{CE_j \partial_j T^{\lambda + 1} + DE_j \partial_j T^{\lambda} T' + CF_{ij} \partial_i \partial_j T^{\lambda + 1} + DF_{ij} \partial_i \partial_j T^{\lambda} T'}

where

E_j = A \left ( \partial_i B_{ij} \right )

F_{ij} = A B_{ij}

Discretization

We then need expressions for

\partial_t T = \frac{ T'_0 - T_0}{\Delta t}

\partial_j T^{\lambda+1} = \frac{T^{\lambda+1}_{\hat{j}} - T^{\lambda+1}_{-\hat{j}}}{2 \Delta x}

\partial_j T^{\lambda}T' = \frac{T^{\lambda}_{\hat{j}}T'_{\hat{j}} - T^{\lambda}_{-\hat{j}}T'_{-\hat{j}}}{2 \Delta x}

\partial_i\partial_j T^{\lambda+1} = \frac{T^{\lambda+1}_{\hat{i}+\hat{j}} - T^{\lambda+1}_{\hat{i}-\hat{j}}-T^{\lambda+1}_{-\hat{i}+\hat{j}} + T^{\lambda+1}_{-\hat{i}-\hat{j}}}{4 \Delta x^2}\left(1-\delta_{ij}\right ) + \frac{T^{\lambda+1}_{\hat{i}} - 2T^{\lambda+1}_{0} + T^{\lambda+1}_{-\hat{i}}}{\Delta x^2}\delta_{ij}

\partial_i\partial_j T^{\lambda}T' = \frac{T^{\lambda}_{\hat{i}+\hat{j}}T'_{\hat{i}+\hat{j}} - T^{\lambda}_{\hat{i}-\hat{j}}T'_{\hat{i}-\hat{j}}-T^{\lambda}_{-\hat{i}+\hat{j}}T'_{-\hat{i}+\hat{j}} + T^{\lambda}_{-\hat{i}-\hat{j}}T'_{-\hat{i}-\hat{j}}}{4 \Delta x^2}\left(1-\delta_{ij}\right ) + \frac{T^{\lambda}_{\hat{i}}T'_{\hat{i}} - 2T^{\lambda}_{0}T'_{0} + T^{\lambda}_{-\hat{i}}T'_{-\hat{i}}}{\Delta x^2}\delta_{ij}

Using the above definitions, we can write the discretized equation as

T_0 + \Delta t \displaystyle \sum_{\parallel, \perp} C \left [ \alpha_0 T^{\lambda+1}_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda+1}_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} \right ] = T'_0 - \Delta t \displaystyle \sum_{\parallel, \perp} D \left [ \alpha_0 T^{\lambda}_0T'_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda}_{\pm \hat{j}} T'_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda}_{\pm \hat{i} \pm \hat{j}} T'_{\pm \hat{i} \pm \hat{j}} \right ]

where

\alpha_0 = -\frac{2F_{ij}\delta_{ij}}{\Delta x^2}

\alpha_{\pm j} = \pm \frac{E_j }{2 \Delta x} + \frac{F_{jj}}{\Delta x^2}

\alpha_{\pm i, \pm j} = \pm \pm \frac{F_{ij}}{4 \Delta x^2}

Summary

	\parallel	\perp
A	\frac{\kappa_\parallel}{\rho c_v \left ( \lambda_\parallel + 1 \right ) }	\frac{\kappa_\perp}{\rho c_v \left (\lambda_\perp + 1\right )}
B_{ij}	b_i b_j	n^2 B^{-2} b_i b_j \left ( 1 - \delta_{ij} \right )
C	\left ( 1 - \phi \left ( \lambda_\parallel + 1 \right ) \right )	\left ( 1 - \phi \left ( \lambda_\perp + 1 \right ) \right )
D	\phi \left ( \lambda_\parallel + 1 \right )	\phi \left ( \lambda_\perp + 1 \right )

And then using those we can calculate

E_j	A \left ( \partial_i B_{ij} \right )
F_{ij}	A B_{ij}
\alpha_0	-\frac{2F_{ij}\delta_{ij} }{\Delta x^2}
\alpha_{\pm j}	\pm \frac{E_j }{2 \Delta x} + \frac{F_{jj}\Delta t}{\Delta x^2}
\alpha_{\pm i, \pm j}	\pm \pm \frac{F_{ij}}{4 \Delta x^2}

Effective diffusion coeffiient ==

The effecitve diffusion coefficient is \frac{n \chi}{\rho c_v} which determines the diffusion time scale etc…

Isotropic Case

For the Isotropic Case, the diffusion equation looks like

\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \kappa T^\lambda \nabla T

or

\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \frac{\kappa}{\lambda+1} \nabla T^{\lambda+1}

or in Einstein notation

\rho c_v \partial_t T = \partial_i \frac{\kappa}{\lambda+1} \partial_i T^{\lambda+1}

which we can rewrite as

\partial_t T = A \partial_i B \partial i T^{\lambda+1}

\partial_t T = A \partial_i B \delta_{ij} \partial j T^{\lambda+1}

This is the exact same form as before, but now

B_{ij}=\kappa \delta{ij}

is diagonal so the cross terms vanish

and \alpha_{\pm i, \pm j} = 0

and it is independent of direction, so it is just a scalar.

E_j	A \left ( \partial_j \kappa \right )
F	A \kappa
\alpha_0	-\frac{2F\delta_{ii} }{\Delta x^2}
\alpha_{\pm j}	\pm \frac{E_j }{2 \Delta x} + \frac{F}{\Delta x^2}
\alpha_{\pm i, \pm j}	0

Also if \kappa is a constant, the E_j terms drop out.

Time stepping

The implicit method requires approximating

T*^{\lambda+1} = T^{\lambda + 1} + \left ( \lambda + 1 \right ) T^{\lambda} \left [ \left ( T_{*} - T \right ) + \frac{\lambda}{2} \frac{T_{*} - T}{T} \right ]

but solving this accurately with a linear method requires

\frac{T_{*}-T}{T} << 1

but this restricts our time step. Using the instantaneous derivatives, we can calculate the maximum time step as

T_0'-T_0 \approx \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 T^{\lambda+1}_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda+1}_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} \right ]

so

\Delta t < \epsilon \frac{T_0}{\displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 T^{\lambda+1}_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda+1}_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} \right ]}

Explicit Case

For the Explicit case, we just set \phi = 0 which sets D = 0 and C=1 and we have

T_0 + \Delta t \displaystyle \sum_{\parallel, \perp} \left [ \alpha_0 T^{\lambda+1}_0 + \displaystyle \sum_{\pm j} \alpha_{\pm j} T^{\lambda+1}_{\pm \hat{j}} + \sum_{\pm i, \pm j,i \ne j} \alpha_{\pm i, \pm j} T^{\lambda+1}_{\pm \hat{i} \pm \hat{j}} \right ] = T'_0

Explicit time stepping

For the explicit scheme to be stable, we must have

\Delta t < \frac{1}{2\displaystyle \max_{\parallel, \perp, i} \left [ | \alpha_i T_i^{\lambda} | \right ]}

Boundary Terms

Explicit

Here, the boundary terms are quite straightforward. We can just use the regular boundary conditions (or user specified boundary conditions). The only difficulty arises when trying to set a heat flux since this involves solving a non-linear equation. One solution is to set the heat flux explicitly - instead of indirectly through the temprature. Then use extrapolating (pr reflecting) BC's for temperature so there is no heat flux initially calculated along the boundary - and then add in the flux in a separate step.

Implicit

For reflecting - or extrapolating boundary terms, we can just transfer the coefficients for the boundary terms and apply them to their mirror zone. The only thing to be careful about - is that for the anisotropic case, there are corner zones who's coefficients may need to be transferred twice. This just requires transferring corner zone coefficients, to edge zones, before transferring edge zones, to the central cell in the stencil.