MHD equations

Units in Astrobear

Astrobear uses something like rationalized electromagnetic units (extra factor of \sqrt{4 \pi} in the electric and magnetic fields) - or Lorentz-Heaviside but scaling the E field by c and the charge density \rho (and current J ) by \frac{1}{c}

AstroBEAR	Lorentz-Heaviside	EMU	Gaussian	Physical Units
E	c E^{LH}	\frac{1}{\sqrt{4\pi}} E^{EMU}	\frac{c }{\sqrt{4\pi}}E^{G}	\mbox{g}^{1/2} \mbox{cm}^{1/2} \mbox{s}^{-2}
\rho	\frac{1}{c}\rho^{LH}	\sqrt{4 \pi}\rho^{EMU}	\frac{\sqrt{4 \pi} }{c}\rho^{G}	\mbox{g}^{1/2} \mbox{cm}^{-5/2}
J	\frac{1}{c} J^{LH}	\sqrt{4 \pi} J^{EMU}	\frac{\sqrt{4 \pi}}{c} J^{G}	\mbox{g}^{1/2} \mbox{cm}^{-3/2} \mbox{s}^{-1}
B	B^{LH}	\frac{1}{\sqrt{4\pi}} B^{EMU}	\frac{1}{\sqrt{4\pi}} B^{G}	\mbox{g}^{1/2} \mbox{cm}^{-1/2} \mbox{s}^{-1}

Or - Using the approach in the appendix of Jackson, we have

k_1 = \frac{c^2}{4 \pi}

k_2 = \frac{1}{4\pi}

k_3 = 1

\alpha = 1

\mu_0 = 1

\epsilon_0 = \frac{1}{c^2}

This allows us to write Maxwell's equations as

\nabla \cdot \mathbf{E} = c^2 \rho

\nabla \times \mathbf{B} = \mathbf{J} + \frac{1}{c^2} \frac{\partial \mathbf{E}}{\partial t}

\nabla \times \mathbf{E} = - \frac{\partial \mathbf{B}}{\partial t}

\nabla \cdot \mathbf{B} = 0

as well as

Lorentz Force Law	\mathbf{F} = q \left ( \mathbf{E} + \mathbf{v} \times \mathbf{B} \right)
Coulomb's Law	\mathbf{F} = -\frac{c^2}{4 \pi} \frac{q_1 q_2}{r^2}\hat{\mathbf{r}}

Most of the time we don't care about E, \rho, or J, however for the Hall MHD terms, we need to determine the electron charge in our system of units.

For our system of equations, the elementary charge is

e = 4.80320425 \times 10^{-10} \mbox{statC} = 1.70269157 \times 10^{-9} \mbox{hsu} = 5.67956774 \times 10^{-20} \left [\mbox{g}^{1/2} \mbox{cm}^{1/2} \right]

Now for MHD we assume that the plasma is neutral (\rho = 0) and that \mathbf{J} >> \frac{1}{c^2}\frac{\partial \mathbf{E}}{\partial t} and instead of tracking the electric field \mathbf{E} through time, we will use Ohm's law to calculate the electric field.

This leaves us with the following set of equations

Induction equation	\frac{\partial \mathbf{B}}{\partial t} = -\nabla \times \mathbf{E}
Ampere's equation (without Maxwell's correction)	\mathbf{J} = \nabla \times \mathbf{B}
Lorentz Force Law	\frac{\partial \rho \mathbf{v}}{\partial t}=\mathbf{J} \times \mathbf{B}
Generalized Ohm's Law	\frac{m_e}{n e^2} \frac{\partial J}{\partial t} - \frac{1}{ne}\nabla \cdot \mathbf{P_e} = \mathbf{E} + \mathbf{v} \times \mathbf{B} - \eta\mathbf{J} - \frac{1}{ne} \mathbf{J} \times \mathbf{B}

In addition we have equations for the total electromagnetic energy involving the net poynting flux \left ( -\nabla \cdot \mathbf{S} \right) and work done on the gas \left (-\mathbf{J} \cdot \mathbf{E} \right) .

Poynting's theorem (electromagnetic energy)	\frac{\partial e}{\partial t} = - \nabla \cdot \mathbf{S} - \mathbf{J}\cdot \mathbf{E}
Poynting vector	\mathbf{S} = \mathbf{E} \times \mathbf{B}

Equations for momentum, magnetic energy, and magnetic field

Subsituting Ampere's equation and the definition of the Poynting vector into the equations for momentum, magnetic energy, and the magnetic field, we arrive at

\frac{\partial \rho \mathbf{v}}{\partial t} = \mathbf{J} \times \mathbf{B} = \left (\nabla \times \mathbf{B} \right) \times \mathbf{B} = \nabla \cdot \left [ \mathbf{BB}-\frac{B^2}{2} \mathbf{I} \right ]

\frac{\partial e}{\partial t} = -\nabla \cdot \left ( \mathbf{E} \times \mathbf{B} \right) - \left ( \nabla \times \mathbf{B} \right) \cdot \mathbf{E}

\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right)

Ohm's Law

Not the momentum equation only involves the magnetic field \mathbf{B} and is the same regardless of the electric field \mathbf{E}. Note it is also already in a form involving a total divergence making it straightforward to implement in a conservative finite volume scheme. The energy and induction equation, however, will depend on our choice of \mathbf{E}.

The \mathbf{E} field can be calculated using a Generalized Ohm's law

\frac{m_e}{n e^2} \frac{\partial J}{\partial t} - \frac{1}{ne}\nabla \cdot \mathbf{P}_e = \mathbf{E} + \mathbf{v} \times \mathbf{B} - \eta\mathbf{J} - \frac{1}{ne} \mathbf{J} \times \mathbf{B}

Under under certain conditions, we can ignore various terms resulting in the following approximations

Resistive + Hall MHD	\mathbf{E} = -\mathbf{v} \times \mathbf{B} + \eta \mathbf{J} + \frac{1}{ne} \mathbf{J} \times \mathbf{B}
Resistive MHD	\mathbf{E} = -\mathbf{v} \times \mathbf{B} + \eta \mathbf{J}
Ideal MHD	\mathbf{E} = - \mathbf{v}\times \mathbf{B}

Also - the Biermann Battery term can be included \left ( \frac{1}{ne} \nabla \cdot \mathbf{P}_e \right) which then allows for magnetic fields to be generated without the need for nascient fields.

Ideal MHD

For ideal MHD we substitute

\mathbf{E} = \mathbf{E_{EMF}} = -\mathbf{v} \times \mathbf{B}

into the equations above and after some vector calculus (including a very messy vector quadruple product with a del operator), manage to write the energy term as a total divergence (and the field as a curl)

\frac{\partial e}{\partial t} = \nabla \cdot \left ( \left (\mathbf{v} \times \mathbf{B} \right ) \times \mathbf{B} \right) + \left ( \nabla \times \mathbf{B} \right) \cdot \left ( \mathbf{v} \times \mathbf{B} \right) = \nabla \cdot \left [ \mathbf{B} \left ( \mathbf{B} \cdot \mathbf{v} \right) - \frac{B^2}{2} \mathbf{v} \right]

\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right)

Resistive MHD

For resistive MHD, it will be convenient to define the resistive component of the electric field \color{green}{\mathbf{E}_{\scriptscriptstyle R} = \eta \mathbf{J} = \eta \nabla \times \mathbf{B}}

Then our total electric field is given by

\mathbf{E} = -\mathbf{v} \times \mathbf{B} \color{green}{+ \mathbf{E}_{\scriptscriptstyle R}}

We then need to add this contribution to the induction equation \color{green}{-\nabla \times \mathbf{E}_{\scriptscriptstyle R}} and the energy equation where it enters both the Poynting flux -\color{green}{\nabla \cdot \left ( \mathbf{E}_{\scriptscriptstyle R} \times \mathbf{B} \right)} and the work done on the the charge distribution \color{green}{ - \mathbf{J} \cdot \mathbf{E}_{\scriptscriptstyle R}}.

However the work done on the charge distribution would need to also be added to the thermal energy cancelling any changes in the total energy due to the \color{green}{\mathbf{J} \cdot \mathbf{E}_{\scriptscriptstyle R}} term. This leaves us with.

\frac{\partial e}{\partial t} = \nabla \cdot \left [ \mathbf{B} \left ( \mathbf{B} \cdot \mathbf{v} \right) - \frac{B^2}{2} \mathbf{v} \color{green}{- \mathbf{E}_{\scriptscriptstyle R} \times \mathbf{B}} \right]

\frac{\partial \mathbf{B}}{\partial t} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right) \color{green}{ - \nabla \times \mathbf{E}_{\scriptscriptstyle R}} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right) \color{green}{- \nabla \times \left ( \eta \nabla \times \mathbf{B} \right )}

Note the resistive part of the induction equation

\frac{\partial \mathbf{B}}{\partial t} = \color{green}{- \nabla \times \left ( \eta \nabla \times \mathbf{B} \right )}

can be expanded as

\frac{\partial \mathbf{B}}{\partial t} = \color{green}{- \nabla \times \left ( \eta \nabla \times \mathbf{B} \right ) = \eta \nabla^2 \mathbf{B} - \vec{\nabla} \mathbf{B} \cdot \nabla \eta}

Which hints to the diffusive nature of the resistive term. It also involves a second derivative so the time stepping criteria will be \Delta t \propto \Delta x^2 and an implicit solve might be worth considering.

Hall+Resistive MHD

As in resistive MHD, it will be convenient to define the Hall component of the electric field \color{red}{\mathbf{E}_{\scriptscriptstyle H} =} \color{red}{ \frac{1}{ne} \mathbf{J} \times \mathbf{B} = \frac{1}{ne} \left ( \nabla \times \mathbf{B} \right ) \times \mathbf{B}}

Then our total electric field is given by

\mathbf{E} = -\mathbf{v} \times \mathbf{B} \color{green}{+ \mathbf{E}_{\scriptscriptstyle R}} \color{red}{+\mathbf{E}_{\scriptscriptstyle H}}

Again we need to add the additional electric field component to the energy equation (in both the Poynting Flux and the work on the charge distribution) as well as the induction equation. Also note in this case the electric field is perpendicular to the current - so it does no work on the charge distribution.

\mathbf{J} \cdot \color{red}{\mathbf{E}_{\scriptscriptstyle H}} = \mathbf{J} \cdot \color{red}{\frac{1}{ne} \mathbf{J} \times \mathbf{B}} = 0

So as in the Resistive case, only the Poynting Flux term needs to be added to the energy equation.

\frac{\partial \rho \mathbf{v}}{\partial t} = \left (\nabla \times \mathbf{B} \right) \times \mathbf{B} = \nabla \cdot \left [ \mathbf{BB}-\frac{B^2}{2} \mathbf{I} \right ]

\frac{\partial e}{\partial t} = \nabla \cdot \left [ \mathbf{B} \left ( \mathbf{B} \cdot \mathbf{v} \right) - \frac{B^2}{2} \mathbf{v} \color{green}{ - \mathbf{E}_{\scriptscriptstyle R} \times \mathbf{B}} \color{red}{ -\mathbf{E}_{\scriptscriptstyle H} \times \mathbf{B}} \right]

\frac{\partial B}{\partial t} = \nabla \times \left ( \mathbf{v} \times \mathbf{B} \right) \color{green}{- \nabla \times \mathbf{E}_{\scriptscriptstyle R}} \color{red}{- \nabla \times \mathbf{E}_{\scriptscriptstyle H}}

Comparing the advective and the Hall terms, we can calculate the length scale for magnetic gradients below which the Hall term becomes larger than the advective term.

L = \frac{B}{n e v}

For the lab experiments, B = 5 T, v = 60 \mbox{km/s}, n=1\times 10^{18} \mbox{cm}^{-3} which works out to give

L = \frac{5\times 10^4 G}{\sqrt{4\pi} \left ( 1\times 10^{18}\right) \left( 6 \times 10^6 \right) \left ( 5.67956774 \times 10^{-20} \right )} \frac{\mbox{g}^{1/2} \mbox{cm}^{-1/2} \mbox{s}^{-1}}{\mbox{cm}^{-3} \mbox{cm}^1 \mbox{s}^{-1} \mbox{g}^{1/2} \mbox{cm}^{1/2} } = .04 \mbox{cm}

If we assume v = v_{A} = \frac{B}{\sqrt{\rho}} then we get

L = \frac{\mu^{1/2}}{n^{1/2} e}

and given the definition of the plasma ion frequency

\omega_{pi}=\left (\frac{4 \pi n Z^2 e^2}{m_i} \right)^{1/2}

Discretization

For resistive MHD, we calculate \mathbf{E} = \eta \nabla \times \mathbf{B} at cell edges for which each component can be calculated by differencing the 4 cell faces that share the edge.

Getting the curl of B at a cell edge is straight forward as it involves difference the 4 adjacent cell faces. For example, the Z component would be given by

E^z_{i+1/2,j+1/2,k} = \eta_{i+1/2,j+1/2,k} \left ( \frac{B^y_{i+1,j+1/2,k} - B^y_{i,j+1/2,k}}{\Delta x} - \frac{B^x_{i+1/2,j+1,k} - B^x_{i+1/2,j,k}}{\Delta y} \right)

For Hall MHD we need to calculate \xi = \left (\nabla \times \mathbf{B} \right) \times \mathbf{B} along with n at cell edges.

\xi_i = B_k d_l B_m \epsilon_{ijk} \epsilon_{jlm} = B_k d_l B_m \left(\delta_{kl} \delta_{im} - \delta_{km} \delta_{il} \right )

So we have

\xi_z = B_x \left ( d_x B_z - d_z B_x \right) + B_y \left ( d_y B_z - d_z B_y \right)

So we need

• edge centered estimates of the perpendicular fields (which can be found by averaging the two adjacent faces)
• transverse gradients of the parallel field (which can be found by differencing adjacent face centered averages of the parallel field (found by averaging the adjacent cell centered fields)
• parallel gradients of the transverse field

Using b to represent cell centered values, B to represent face centered values, we have

\xi^z_{i+1/2,j+1/2,k} = \frac{B^x_{i+1/2,j,k} + B^x_{i+1/2,j+1,k}}{2} \left ( \frac{b^z_{i+1,j,k} + b^z_{i+1,j+1,k} - b^z_{i,j,k}-b^z_{i,j+1,k}}{2 \Delta x} - \frac{B^x_{i+1/2, j, k+1} + B^x_{i+1/2,j+1,k+1}-B^x_{i+1/2,j,k-1}-B^x_{i+1/2,j+1,k-1}}{4 \Delta z} \right )

+\frac{B^y_{i,j+1/2,k} + B^y_{i+1,j+1/2,k}}{2} \left ( \frac{b^z_{i,j+1,k} + b^z_{i+1,j+1,k} - b^z_{i,j,k}-b^z_{i+1,j,k}}{2 \Delta y} - \frac{B^y_{i, j+1/2, k+1} + B^y_{i+1,j+1/2,k+1}-B^y_{i,j+1/2,k-1}-B^y_{i+1,j+1/2,k-1}}{4 \Delta z} \right )

Now if we represent edge centered perpendicular B fields \beta^x_{i+1/2,j+1/2,k} = \frac{B^x_{i+1/2,j,k}+B^x_{i+1/2,j+1,k}}{2}

and face centered parallel fields

\alpha^z_{i,j+1/2,k} = \frac{b^z_{i,j,k}+b^z_{i,j+1,k}}{2}

this simplifies to

\xi^z_{i+1/2,j+1/2,k} = \beta^x_{i+1/2,j+1/2,k} \left ( \frac{\alpha^z_{i+1,j+1/2,k} - \alpha^z_{i,j+1/2,k}}{\Delta x} - \frac{\beta^x_{i+1/2, j+1/2, k+1} -\beta^x_{i+1/2,j+1/2,k-1}}{2 \Delta z} \right )

+\beta^y_{i+1/2,j+1/2,k} \left ( \frac{\alpha^z_{i+1/2,j+1,k} - \alpha^z_{i+1/2,j,k}}{\Delta x} - \frac{\beta^y_{i+1/2, j+1/2, k+1} -\beta^y_{i+1/2,j+1/2,k-1}}{2 \Delta z} \right )

So to calculate 'z' edges we need adjacent face centered z fields and edge centered transverse fields

These fields won't be used for any of the other directional updates - so can be discarded after each dimension