===================================

Viscosity and Resistivity in AstroBEAR

===================================

Viscosity

The viscosity term in the Navier-Stokes Equation can be separated as:

dv/dt = \nabla \cdot T
For the applications we are considering, the fluid is isotropic. Under this assumption, the tensor can be expressed as the following combination of two scalar dynamic viscosities \mu and \mu':

T = 2 \mu \nabla v + \mu' \nabla \cdot v I

where I is the unit tensor.
It can be proven that the stress tensor should have zero trace, which leads to the following relation between \mu and \mu':

\mu' = -2/3\mu

Combining these relations, we can write down the final form of the separated viscous equation:

dv/dt=\nabla \cdot (2 \mu \nabla v) - 2/3 \nabla (\mu \nabla \cdot v)

Notice that this equation already counts the effect of all the elements in a viscosity tensor, assuming the fluid itself is isotropic.

Resistivity

The resistivity equation can be written as:
\frac{\partial \textbf{B}}{\partial t}=\nabla \times (\eta \nabla \times \textbf{B})

Here, the magnetic field are located at face centers. Using the face centered magnetic field, we can calculate the edge centered currents, which then give a source term.
The curl of the edge centered source term gives a face centered \delta B, which can then be applied to the face centered magnetic field.

===================================

Implementation

===================================

Viscosity

The viscous source term is calculated rigorously by first finding the tensor elements at the surrounding face centers of a cell center considered: The \mu sources 2\mu \nabla v are naturally located at the cell centers, the \mu' sources however needs an averaging to obtain the cell centered velocity divergence.
Therefore the scheme requires a 20 point stencil in 3-D configuration, which makes the implicit scheme a bit complicated. The following graph shows such a stencil in 2-D setting:



Using this stencil setting, one can get face centered velocity divergence for the cell considered, and then calculate the gradient components located at the cell center.
In the viscous plasma context, the physical viscosity can be calculated rigorously. The plasma viscosity is proportional to the product of ion mean free path and the ion thermal energy:

\mu \sim \lambda_{ifp}v_{thermal}

In our code, we offer the option to calculate the realistic viscosity using the following equation:

\mu = \frac{m_{i}}{3 \sqrt{2} \pi d_{i}^2} (\frac{8kT}{\pi m_{i}})^{\frac{1}{2}}

The ion mass and local temperature should be readily available, but a good estimate of ion diameter d_{i} usually requires look up a ionic radius table.
As an estimate, one can use that of hydrogen, about 0.1~0.2 nm.
One can see that the viscosity itself is anisotropic and depends on the local temperature.
That is why we have both the \mu and \mu' terms in our viscosity equation, and have to calculate face centered source terms before applying the final divergence.

The viscosity interface is as follows:
The viscosity source term can be turned on in astrobear by putting the viscosity flag to true:
lViscous = .true.
The second parameter VisType determines which type of viscosity one requires: 1 for constant viscosity; 2 for real physical viscosity.
The third parameter viscosity is the user defined constant viscosity.
The viscosity module is tested to be working compatibly with implicit thermal conduction module and uniform gravity module.

Resistivity

For each cell, the resistive source term is calculated by taking the curl of face centered magnetic field. This results in a 12 dimension array that stores the current source terms running over the 12 edges of a cell in 3-D. This array is then used to calculate the new face centered magnetic field. The stencil for this explicit solver is a 3 by 3 cube surrounding the cell we want to update. See the following diagram:



In cgs units, the real magnetic diffusivity can be written as:

\eta = c^2/4\pi \sigma

where \sigma is the material conductivity which depends on the electron-ion collision rate \nu_{ei} which in term depends on the electron temperature T^{-3/2}.
The cross field diffusivity can thus be expressed as function of electron temperature

\eta = 0.14 m_{e}^{1/2} c^2 Z \ln \Lambda (k_B T_{e})^{-3/2}

where \ln \Lambda is the Coulomb logarithm:

\ln \Lambda = 16.3+1.5 \ln T - 0.5\ln n$, (for $T < 4.2 \times 10^{5}$ K)\\$\ln \Lambda =$$22.8+1.5 \ln T - 0.5\ln n$, (for $T > 4.2 \times 10^{5}$ K)

THe Coulomb logarithm for interested number density is plotted below:

The parallel diffusivity however, is dependent on

F(Z) = \frac{1+1.198Z+0.222Z^2}{1+2.996Z+0.753Z^2}

For Z = 1 case, parallel diffusivity is about half of cross field diffusivity.

Currently in the code, the realistic resistivity gets calculated at each cell edges by averaging the temperatures surrounding that edge. So the temperature effect is taken into full consideration. But the Coulomb log is assigned as a constant, although it may change depending on the temperature and particle density.

The resistivity interface is as follows:
The resistive source term can be turned on in astrobear by putting the resistive flag to true:
lresistive = .true.
The second parameter ResType determines which type of resistivity one requires: 1 for constant computational resistivity; 2 for real physical resistivity.
The third parameter resistivity is the user defined constant resistivity.

===================================

Multiphysics Integration

===================================

Taking viscosity as an example. The goal for the multiphysics modules is to calculate a change of velocisty inside a grid, using the known hydro vars.
This integration can be done before any hyperbolic updates. We denote the hyperbolic update, integrated source term calculation and explicit source term calculation as H, S and E.
The operators applied are as follows:

E S^{1/2} H S^(1/2)

But this operator does not commute with its own transpose, meaning that

E S^{1/2} H S^(1/2) \neq S^{1/2} H S^(1/2) E

This is obvious with for example, viscosity, in which the E operator raises temperature and lowers the velocity profile. The hyperbolic solver will have a different response to the data modified by E comparing to that of unmodified.
A symmetric operator splitting can be achieved if we do

E^{1/2} S^{1/2} H S^(1/2)E^{1/2}

Since the commutator of this operator with its transpose will be zero.
For any given grid, the external explicit solver requires two steps of updating. The tricky part of the solver is to figure out what to use at the grid boundaries. The three different boundary situations are summarized below:
(1)For any physical boundary, use the physical boundary condition.(outflow: impose a fixed velocity; and so on)
(2)For any internal boundary to direct neighbors, get the needed data directly by communication.
(3)For any internal boundary to parent level, we need the following ghosting scheme:



When advancing dt, we have 10 ghost zones, we also need to explicitly store \dot{v}_a dt/2, the time derivative of velocity advection to be used for adjacent child grids.
The steps are:
1) apply E(dt/2) to internal + ghost 1~8 using ghost 9 and 10.
2) apply H(dt/2) to internal + ghost 1~4 using ghost 5~8.
3) apply another H(dt/2) to internal using ghost 1~4.
4) manually update ghost 1~2 using \dot{v}_a dt/2.
5) apply E(dt/2) to internal again using the updated ghost 1~2.

From here, we can see two problems of putting the explicit operator into the source term integration:
1) each complete update requires two more ghost zones on parent-child boundaries which does not fit into the source term subcyclilng. (N subcycling requires 2N ghost zones)
2) for internal boundary to direct neighbors, the two grids need to communicate in order to . This again does not fit into the source term subcycling. (there is no communications happening during subcycling in the current code)