Skip to content

Governing Equations

Continuity

∂ρ∂t+v⋅∇ρ=−ρ∇⋅v \pdv{\rho}{t} + v \cdot \grad \rho = - \rho \grad \cdot v

Momentum:

ρ(∂v∂t+v⋅∇v)+∇p−1μ0(∇×B)×B=∇⋅T \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

T= viscous stress tensor  \vec T = \text{ viscous stress tensor }

1μ0∇×B=j= current density  \frac{1}{\mu_0} \curl \vec B = \vec j = \text{ current density }

Magnetic Field

∂B∂t−∇×(v×B)=−∇×[ημ0(∇×B)] \pdv{B}{t} - \curl ( v \cross B) = - \curl \left[ \frac{\eta}{\mu_0} (\curl \vec B) \right]

η= electrical resistivity  \eta = \text{ electrical resistivity }

Energy:

ρ(∂E∂t+v⋅∇E)+p∇⋅v=∇⋅(v⋅T)+∇⋅(k∇⋅T) (thermal conduction) +ημ02(∇×B)2 (resistive heating)  \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}

E= specific internal energy \mathscr{E} = \text{ specific internal energy}

ρE=pγ−1 \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

∂f∂t+v⋅∂f∂x+a∂f∂v=∂f∂t∣collisions \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 v \vec v of the Boltzmann equation to get all of the above:

u(x,t)=∫vf(x,v,t)dv∫f(x,v,t)dv \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 Kn is the mean free path between collisions

Kn=λmfpL 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 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

∂∂t∫vnfdv+∫vn+1fdv \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 (E \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 ∇⋅B=0 \div B = 0 , we can re-write the equations in conservative form.

∂ρ∂t+∇⋅(ρv)=0 \pdv{\rho}{t} + \div (\rho v) = 0

∂(ρv)∂t+∇⋅[ρvv−BBμ0+1(B22μ0+p)]=∇⋅T \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

∂B∂t+∇⋅[vB−Bv]=∇⋅(ημ0∇B) \pdv{B}{t} + \div [v B - B v] = \div \left( \frac{\eta}{\mu_0} \grad B \right)

∂e∂t+∇⋅[(e+p+B22μ0)v−(B⋅v)Bμ0]=∇⋅[v⋅T+k∇T−ημ0(∇×B)×B] \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=pγ−1+rhov22+B22μ0 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 AB \vec A \vec B is notation for

AB=[AxBxAxByAxBzAyBxAyBy…………] \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 ABT \vec A \vec B^T

ABT=[AxAyAz][BxByBz] \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,

∂Q∂t+∇⋅F=∇⋅Fp(parabolic) \pdv{Q}{t} + \div \vec F = \div \vec F_p \quad \text{(parabolic)}

or

∂Q∂t+∂F∂x+∂G∂y+∂H∂z=∂Fp∂x+∂Gp∂y+∂Hp∂z \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:

∂Q∂t=−∇⋅F \pdv{Q}{t} = - \div \vec F

Our governing relations are not in conservation form, but they can be. Take the conservation of mass relation

∂ρ∂t=−∇⋅(ρv) \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

∂∂t(ρv)+∇⋅[ρvv−BBμ0+(B22μ0+p)1]=∇⋅T \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

∂B∂t+∇⋅[vB−Bv]=∇⋅[ημ0∇B] \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

∂e∂t+∇⋅[(e+p+B22μ0)v−(B⋅v)Bμ0]=∇⋅[…] \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]

 total energy e=pγ−1+ρv22+B22μ0 \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

∂Q‾∂t+∇⋅F‾‾=∇⋅Fp‾‾ \pdv{\overline{Q}}{t} + \div \overline{\overline{F}} = \div \overline{\overline{F_p}}

where Q‾ \overline{Q} is just the vector containing all of our scalar conserved quantities, and F‾‾ \overline{\overline{F} } is a vector of all of our flux vectors.

Q‾=[ρρμρvρwBxByBze] \overline{Q} = \begin{bmatrix} \rho \\ \rho \mu \\ \rho v \\ \rho w \\ B_x \\ B_y \\ B_z \\ e \end{bmatrix}

F‾‾=F+G+H=Fx‾‾+Fy‾‾+Fz‾‾ \overline{\overline{F}} = \vec F + \vec G + \vec H = \overline{\overline{F_x}} + \overline{\overline{F_y}} + \overline{\overline{F_z}}

F‾=[ρuρu2−Bx2/2μ+pρuv−BxBy/μ0ρuw−BxBz/μ0uBy−BxvuBz−Bxw(e+p+B2/2μ0)u−B⋅vμ0Bx]= flux vector (x-component)  \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 G‾ \overline{G} and H‾ \overline{H} . Now let's move on to Fp‾ \overline{F_p} , the parabolic flux vector (as opposed to the hyperbolic fluxes F‾ \overline{F} ). These fluxes typically come from "non-ideal" effects.

Fp‾=[02μ(∂u∂x−∇⋅v3)μ(∂μ∂y+∂v∂x)μ(∂u∂z+∂w∂z)ημ0∂Bx∂xημ0∂By∂xημ0∂Bz∂xuTxx+vTyx+wTzx+k∂T∂x−ημ0[(∇×B)yBz−(∇×B)zBy]] \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

∂Q∂t+∂F∂x+∂G∂y+∂H∂z=0 \pdv{Q}{t} + \pdv{F}{x} + \pdv{G}{y} + \pdv{H}{z} = 0

or, in 2 dimensions

∂Q∂t+∂F∂x+∂G∂y=0 \pdv{Q}{t} + \pdv{F}{x} + \pdv{G}{y} = 0

Q=[ρρuρve] Q = \begin{bmatrix} \rho \\ \rho u \\ \rho v \\ e \end{bmatrix}

F=[ρuρu2+pρuvu(e+p)] F = \begin{bmatrix} \rho u \\ \rho u ^2 + p \\ \rho u v \\ u (e + p) \end{bmatrix}

G=[ρvρuvρv2+pv(e+p)] 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.