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 looks like

\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}

\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 [ n \hat{b} \left ( \kappa_\parallel T^{\lambda_\parallel} - \frac{n \kappa_\perp}{B^2} T^{\lambda_\perp} \right ) \left ( \hat{b} \cdot \nabla T \right ) + \frac{n \kappa_\perp}{B^2} T^{\lambda_\perp} \nabla T \right ]

or

\rho c_v \frac{\partial T}{\partial t} = \nabla \cdot \left [ n \hat{b} \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( \hat{b} \cdot \nabla T^{\lambda_\parallel+1} \right ) - \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( \hat{b} \cdot \nabla T^{\lambda_\perp + 1} \right ) \right ) + \frac{n \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 [ n b_i \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( b_j \partial_j T^{\lambda_\parallel+1} \right ) - \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( b_j \partial_j T^{\lambda_\perp + 1} \right ) \right ) + \frac{n \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 [ n b_i \left ( \frac{\kappa_\parallel}{\lambda_\parallel+1} \left ( b_j \partial_j T^{\lambda_\parallel+1} \right ) - \frac{n \kappa_\perp}{B^2 \left ( \lambda_\perp+1 \right )} \left ( b_j \partial_j T^{\lambda_\perp + 1} \right ) \right ) + \frac{n \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 n b_i \frac{\kappa_\parallel}{\lambda_\parallel+1} b_j \partial_j T^{\lambda_\parallel+1} + \partial_i n \left ( \delta_{ij} - b_i b_j \right ) \frac{n \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 n 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} = n 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 - \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 + \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 t}{\Delta x^2}

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

\alpha_{\pm i, \pm j} = \pm \pm \frac{F_{ij}\Delta t}{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}	n 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 t}{\Delta x^2}
\alpha_{\pm j}	\pm \frac{E_j \Delta t}{2 \Delta x} + \frac{F_{jj}\Delta t}{\Delta x^2}
\alpha_{\pm i, \pm j}	\pm \pm \frac{F_{ij}\Delta t}{4 \Delta x^2}

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 t}{\Delta x^2}
\alpha_{\pm j}	\pm \frac{E_j \Delta t}{2 \Delta x} + \frac{F\Delta t}{\Delta x^2}
\alpha_{\pm i, \pm j}	0

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

Explicit Case

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

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 ] = T'_0