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

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 vdot inside a grid, using the known hydro vars. It will need the two ghost zone's as illustrated below:

This integration can be done before any hyperbolic updates.