Governing Equations #
Continuity
\[\pdv{\rho}{t} + v \cdot \grad \rho = - \rho \grad \cdot v\]Momentum:
\[\rho \left( \pdv{v}{t} + v \cdot \grad v\right) + \grad p - \frac{1}{\mu_0} (\curl \vec B) \cross \vec B = \div \vec T \] \[\vec T = \text{ viscous stress tensor }\] \[\frac{1}{\mu_0} \curl \vec B = \vec j = \text{ current density }\]Magnetic Field
\[\pdv{B}{t} - \curl ( v \cross B) = - \curl \left[ \frac{\eta}{\mu_0} (\curl \vec B) \right]\] \[\eta = \text{ electrical resistivity }\]Energy:
\[\begin{aligned} \rho \left( \pdv{\mathscr{E}}{t} + v \cdot \grad \mathscr{E} \right) + p \div v & = \div \left( v \cdot \vec T \right) & \\ & + \div ( k \div T) \qquad \text{ (thermal conduction) }& \\ & + \frac{\eta}{\mu_0 ^2} (\curl B )^2 \qquad \text{ (resistive heating) } & \end{aligned}\] \[\mathscr{E} = \text{ specific internal energy}\] \[\rho \mathscr{E} = \frac{p}{\gamma - 1}\]These are the governing equations for magnetohydrodynamics (MHD). It’s important to understand where they come from. Starting from the Boltzmann equation, which describes how a distribution evolves over time
\[\pdv{f}{t} + \vec v \cdot \pdv{f}{\vec x} + \vec a \pdv{f}{\vec v} = \left. \pdv{f}{t} \right|_{\text{collisions}}\]We take moments of \( \vec v \) of the Boltzmann equation to get all of the above:
\[\vec u (\vec x , t) = \frac{\int \vec v f( \vec x , \vec v, t) \dd \vec v}{\int f(\vec x, \vec v, t) \dd \vec v}\]The 0th moment gives the conservation of mass relation, since . The Knudsen number \( Kn \) is the mean free path between collisions
\[Kn = \frac{ \lambda_{mfp}}{L}\]When the Knudsen number is sufficiently small, these equations work well. When Kn becomes larger than the simulation volume, then the equations of fluid dynamics break down and we get non-Maxwellian phenomena. Also, if \( L \) is very small (for example, a powerful shock), we’ll probably miss something.
When taking moments, notice that each moment will introduce the next moment
\[\pdv{}{t} \int \vec v ^{n} f \dd \vec v + \int \vec v ^{n+1} f \dd \vec v\]We get to stop when we define a closure relation that relates moments to higher moments. For example, the ideal gas law is a closure relation, since it relates temperature (\( \mathscr{E} \) to density. Implicit in these equations are the closure relations that we need to keep in mind.
If we drop all terms with B and the magnetic field equation, then we get the Navier-Stokes equations. These equations are a system and must be solved self-consistently. The RHS of each equation represents non-ideal contributions due to self collisions. Those include viscosity, resistivity, and conduction.
Using vector identities and the fact that \( \div B = 0 \), we can re-write the equations in conservative form.
\[\pdv{\rho}{t} + \div (\rho v) = 0\] \[\pdv{(\rho v)}{t} + \div \left[\rho v v - \frac{BB}{\mu_0} + \vec 1 \left( \frac{B^2}{2 \mu_0} + p \right) \right] = \div \vec T\] \[\pdv{B}{t} + \div [v B - B v] = \div \left( \frac{\eta}{\mu_0} \grad B \right)\] \[\pdv{e}{t} + \div \left[ \left( e + p + \frac{B^2}{2 \mu_0} \right) v - \frac{(B \cdot v) B}{\mu_0} \right] = \div \left[ v \cdot \vec T + k \grad T - \frac{\eta}{\mu_0} ( \curl B) \cross B \right]\]where the total energy \( e = \frac{p}{\gamma - 1} + \frac{rho v^2}{2} + \frac{B^2}{2 \mu_0 } \).
As a quick refresher on dyad notation, we say \( \vec A \vec B \) is notation for
\[\vec A \vec B = \begin{bmatrix} A_x B_x & A_x B_y & A_x B_z \\ A_y B_x & A_y B_y & \ldots \\ \ldots & \ldots & \ldots & \end{bmatrix}\]For that reason, sometimes it’ll be written as a matrix multiplication \( \vec A \vec B^T \)
\[\vec A \vec B^T = \begin{bmatrix} A_x \\ A_y \\ A_z \end{bmatrix} \begin{bmatrix} B_x & B_y & B_z \end{bmatrix}\]In compact form,
\[\pdv{Q}{t} + \div \vec F = \div \vec F_p \quad \text{(parabolic)}\]or
\[\pdv{Q}{t} + \pdv{\vec F}{x} + \pdv{\vec G}{y} + \pdv{\vec H}{z} = \pdv{\vec F_p}{x} + \pdv{\vec G_p}{y} + \pdv{\vec H_p}{z}\]Now, we can re-write the governing equations in a form that’s easier to work with, called conservation form. This form involves relating the time derivative of a quantity we are interested in to the divergence of a flux:
\[\pdv{Q}{t} = - \div \vec F\]Our governing relations are not in conservation form, but they can be. Take the conservation of mass relation
\[\pdv{\rho}{t} = - \div ( \rho \vec v)\]Writing our expressions in conservation form helps to ensure that numerical codes replicate real physics. If our simulation does not conserve something like mass or momentum due to some numerical errors, we’re not likely to get a physical result from the simulation.
On to the conservation of momentum, combining with our mass conservation relation
\[\pdv{}{t}(\rho \vec v) + \div \left[ \rho \vec v \vec v - \frac{\vec B \vec B}{\mu_0} + (\frac{B^2}{2 \mu_0} + p ) \vec 1 \right] = \div \vec T\]and the magnetic field equation
\[\pdv{\vec B}{t} + \div \left[ \vec v \vec B - \vec B \vec v \right] = \div \left[ \frac{\eta}{\mu_0} \grad \vec B \right]\]and the energy relation
\[\pdv{e}{t} + \div \left[ \left(e + p + \frac{B^2}{2 \mu_0}\right) \vec v - \frac{( \vec B \cdot \vec v) \vec B}{\mu_0} \right] = \div \left[ \ldots \right]\] \[\text{ total energy } \quad e = \frac{p}{\gamma - 1} + \frac{\rho v^2}{2} + \frac{B^2}{2\mu_0}\]So all of our governing equations can be written in conservation form. We can get to a more compact form which looks like
\[\pdv{\overline{Q}}{t} + \div \overline{\overline{F}} = \div \overline{\overline{F_p}}\]where \( \overline{Q} \) is just the vector containing all of our scalar conserved quantities, and \( \overline{\overline{F} }\) is a vector of all of our flux vectors.
\[\overline{Q} = \begin{bmatrix} \rho \\ \rho \mu \\ \rho v \\ \rho w \\ B_x \\ B_y \\ B_z \\ e \end{bmatrix}\] \[\overline{\overline{F}} = \vec F + \vec G + \vec H = \overline{\overline{F_x}} + \overline{\overline{F_y}} + \overline{\overline{F_z}}\] \[\overline{F} = \begin{bmatrix} \rho u \\ \rho u^2 - B^2 _x / 2 \mu + p \\ \rho u v - B_x B_y / \mu_0 \\ \rho u w - B_x B_z / \mu \\ 0 \\ u B_y - B_x v \\ u B_z - B_x w \\ (e + p + B^2 / 2 \mu_0 )u - \frac{\vec B \cdot \vec v}{\mu_0} B_x \end{bmatrix} = \text{ flux vector (x-component) }\]We can similarly calculate the other factors \( \overline{G} \) and \( \overline{H} \). Now let’s move on to \( \overline{F_p} \), the parabolic flux vector (as opposed to the hyperbolic fluxes \( \overline{F} \)). These fluxes typically come from “non-ideal” effects.
\[\overline{F_p} = \begin{bmatrix} 0 \\ 2 \mu \left(\pdv{u}{x} - \frac{ \div \vec v}{3} \right) \\ \mu \left( \pdv{\mu}{y} + \pdv{v}{x} \right) \\ \mu \left( \pdv{u}{z} + \pdv{w}{z} \right) \\ \frac{\eta}{\mu_0} \pdv{B_x}{x} \\ \frac{\eta}{\mu_0} \pdv{B_y}{x} \\ \frac{\eta}{\mu_0} \pdv{B_z}{x} \\ u \vec T _{xx} + v T_{yx} + w T_{zx} + k \pdv{T}{x} - \frac{\eta}{\mu_0} \left[ ( \curl B)_y B_z - ( \curl B)_z B_y \right] \end{bmatrix}\]20.2.2 The Euler Equations #
If we have a fluid which is unmagnetized, inviscid, and thermal conductivity is zero, we are left with the Euler equations
\[\pdv{Q}{t} + \pdv{F}{x} + \pdv{G}{y} + \pdv{H}{z} = 0\]or, in 2 dimensions
\[\pdv{Q}{t} + \pdv{F}{x} + \pdv{G}{y} = 0\] \[Q = \begin{bmatrix} \rho \\ \rho u \\ \rho v \\ e \end{bmatrix}\] \[F = \begin{bmatrix} \rho u \\ \rho u ^2 + p \\ \rho u v \\ u (e + p) \end{bmatrix}\] \[G = \begin{bmatrix} \rho v \\ \rho u v \\ \rho v^2 + p \\ v(e + p) \end{bmatrix}\]The Euler equations are the inviscid Navier-Stokes equations, and they exhibit many of the nonlinearities that the full N-S and other more complicated flows do, but with a much simpler form. Solution of the Euler equations represents an important advance for any algorithm. Some such important nonlinearities are shockwaves, contact discontinuities, rarefications. It turns out that most of those parabolic fluxes make the situation easier from a CFD standpoint, because they are dissipating effects which actually serve to suppress instabilities.