Basic Equations of Ideal One-Fluid Magnetohydrodynamics (Part III & IV)

FEATURED IMAGE CREDIT: U.R. Christensen from the Nature article: “Earth Science: A Sheet-Metal Dynamo”

The image shows the overall distortion of magnetic field lines inside the core and its’ effect on the magnetic field outside. 

This post shall continue to derive the principal equations of ideal one-fluid magnetohydrodynamics. Here I shall derive the continuity equation and the vorticity equation. Furthermore, I shall show that the Boussinesq approximation results in a zero-valued divergence. I have consulted the following works while researching this topic:

Davidson, P.A., 2001. An Introduction to Magnetohydrodynamics. 3-4,6. 

Murphy, N., 2016. Ideal Magnetohydrodynamics. Lecture presentation. Harvard-Smithsonian Center for Astrophysics. 

Consider a fluid element through which fluid passes. The mass of this element can be represented as the volume integral of the material density \rho:

\displaystyle m=\iiint \rho dV. (1)

Take the first order time derivative of the mass, we arrive at

\displaystyle \dot{m}=\frac{d}{dt}\iiint\rho dV = \iiint \frac{\partial \rho}{\partial t}dV. (2).

We now make a slight notation change; let the triple integral be represented as

\displaystyle \dot{m} = \int_{V}\frac{\partial \rho}{\partial t}dV. (3)

Now, the mass flux through a surface element dV = \hat{n}dV is \rho \textbf{v} \cdot dV. Thus, the integral of the mass flux is

\displaystyle \int_{V}(\rho\cdot \textbf{v})dV = \dot{m}. (4)

We may write this as

-\displaystyle \oint_{\partial V}(\rho\cdot \textbf{v})dV= \int_{V}\frac{\partial \rho}{\partial t}dV. (5)

Now, we invoke Gauss’s theorem (also known as the Divergence Theorem) of the form

\displaystyle \int_{V}(\nabla \cdot \rho \textbf{v})dV=-\oint_{\partial V}(\rho \cdot \textbf{v})dV, (6)

we get

\displaystyle \int_{V}\bigg\{\frac{\partial \rho}{\partial t}+(\nabla \cdot \rho \textbf{v})\bigg\}dV=0. (7)

Since the integral cannot be zero, the integrand must be. Therefore, we get the continuity equation:

\displaystyle \frac{\partial \rho}{\partial t}+ (\nabla \cdot \rho \textbf{v})=0. (8)

Recall the following equation from a previous post (specifically Eq. (6) of that post)

\displaystyle \frac{\partial \textbf{v}}{\partial t}+(\nabla\cdot \textbf{v})\textbf{v}=-\frac{1}{\rho}\nabla\bigg\{P+\frac{B^{2}}{2\mu_{0}}\bigg\}+\frac{(\nabla\cdot \textbf{B})\textbf{B}}{\mu_{0}\rho}, (8)

Now we define the concept of vorticity. Conceptually, this refers to the rotation of the fluid within its velocity field. Mathematically, we define the vorticity \Omega to be the curl of the fluid velocity:

\Omega \equiv \nabla \times \textbf{v}. (9)

Now recall, the vector identity

\displaystyle \nabla \frac{\textbf{v}^{2}}{2}=(\nabla \cdot \textbf{v})\textbf{v}+\textbf{v}\times \Omega, (10.1)

which upon rearrangement is

(\nabla \cdot \textbf{v})\textbf{v}=\nabla \frac{\textbf{v}^{2}}{2}-\textbf{v}\times\Omega. (10.2)

Using Eq.(10.2) with Eq.(8) we get

\displaystyle \frac{\partial \textbf{v}}{\partial t}=\textbf{v} \times \Omega -\nabla \bigg\{\frac{P}{\rho}-\frac{\textbf{v}^{2}}{2} \bigg\}+\nu \nabla^{2} \textbf{v}. (11)

Recall our definition of vorticity. Upon taking the curl of Eq.(11) we arrive at a variation of the induction equation (see post)

\displaystyle \frac{\partial \Omega}{\partial t}=\nabla \times (\textbf{v} \times \Omega)+\nu \nabla^{2}\Omega. (12.1)

Next, we invoke another vector identity

\displaystyle \nabla \times (\textbf{v} \times \Omega) = (\Omega \cdot \nabla)\textbf{v}-(\textbf{v} \cdot \nabla)\Omega (12.2)

Using Eq. (12.2) in Eq. (12.1) yields the vorticity equation

\displaystyle \frac{\partial \Omega}{\partial t}+(\nabla \cdot \textbf{v})\Omega = (\nabla \cdot \Omega)\textbf{v}+\nu \nabla^{2}\Omega. (13)

Returning to the continuity equation:

\displaystyle \frac{\partial \rho}{\partial t}+(\nabla \cdot \rho \textbf{v})=0.

I will now show that the flow does not diverge. In other words, there are no sources in the fluid velocity field. The Boussinesq approximation’s main assertion is that of isopycnal flow (i.e. flow of constant density). Therefore, let \rho = \rho_{0}>0. Substitution into the continuity equation yields the following

\displaystyle 0+ \rho_{0}(\nabla \cdot \textbf{v})=0. (14)

The temporal derivative is a derivative of a non-zero constant which is itself zero. This just leaves \rho_{0}(\nabla \cdot \textbf{v})=0. Now, since \rho_{0}\neq 0, this then means that

\displaystyle (\nabla \cdot \textbf{v})=0. (15)

Hence, the divergence of the fluid’s velocity field is zero.

Solution to Laplace’s Equation

This post deals with the familiar (to the physics student) Laplace’s equation. I am solving this equation in the context of physics, instead of a pure mathematical perspective. This problem is considered most extensively in the context of electrostatics. This equation is usually considered in the spherical polar coordinate system. A lot of finer details are also considered in a mathematical physics course where the topic of spherical harmonics is discussed. Assuming spherical-polar coordinates, Laplace’s equation is

\displaystyle \frac{1}{r^{2}}\frac{\partial}{\partial r}\bigg\{r^{2}\frac{\partial \psi}{\partial r}\bigg\}+\frac{1}{r^{2}\sin{\theta}}\frac{\partial}{\partial \theta}\bigg\{\sin{\theta}\frac{\partial \psi}{\partial \theta}\bigg\}+\frac{1}{r^{2}\sin{\theta}}\frac{\partial^{2}\psi}{\partial \phi^{2}}=0, (1)

Suppose that the function \psi\rightarrow \psi(r,\theta,\phi). Furthermore, by the method of separation of variables we suppose that the solution is a product of eigenfunctions of the form R(r)Y(\theta,\phi). Hence, Laplace’s equation becomes

\displaystyle \frac{Y}{r^{2}}\frac{d}{dr}\bigg\{r^{2}\frac{d^{2}R}{dr^{2}}\bigg\}+\frac{R}{r^{2}\sin{\theta}}\frac{\partial}{\partial \theta}\bigg\{\sin{\theta}\frac{\partial Y}{\partial \theta}\bigg\}+\frac{R}{r^{2}\sin^{\theta}}\bigg\{\frac{\partial^{2}Y}{\partial \phi^{2}}\bigg\}=0. (2)

Furthermore, we can separate further the term Y(\theta, \phi) into \Theta(\theta)\Phi(\phi). Rewriting (2) and multiplying by r^{2}\sin^{2}{\theta}R^{-1}Y^{-1}, we get

\displaystyle \frac{1}{R}\frac{d}{dr}\bigg\{r^{2}\frac{dR}{dr}\bigg\}+\frac{\sin{\theta}}{\Theta}\frac{d}{d\theta}\bigg\{\sin{\theta}\frac{d\Theta}{d\theta}\bigg\}+\frac{1}{\Phi}\frac{d^{2}\Phi}{d\phi^{2}}=0.(3)

Bringing the radial and angular component to the other side of the equation and setting the azimuthal component equal to a separation constant -m^{2}, yielding

\displaystyle -\frac{1}{R}\frac{d}{dr}\bigg\{r^{2}\frac{dR}{dr}\bigg\}-\frac{\sin{\theta}}{\Theta}\frac{d}{d\theta}\bigg\{\sin{\theta}\frac{d\Theta}{d\theta}\bigg\}=\frac{1}{\Phi}\frac{d^{2}\Phi}{d\phi^{2}}=-m^{2}.

Solving the right-hand side of the equation we get

\displaystyle \Phi(\phi)=A\exp({+im\phi}). (4)

Now, we set the azimuthal component equal to m^{2} and carry out the derivative in the angular component to get the Associated Legendre equation:

\displaystyle \frac{d}{d\mu}\bigg\{(1-\mu^{2})\frac{d\Theta}{d\mu}\bigg\}+\bigg\{l(l+1)-\frac{m^{2}}{1-\mu^{2}}\bigg\}\Theta=0, (5)

where I have let \mu=\cos{\theta}. The solutions to this equation are known as the associated Legendre functions. However, instead of solving this difficult equation by brute force methods (i.e. power series method), we consider the case for which m^{2}=0. In this case, Eq.(5) simplifies to Legendre’s differential equation discussed previously. Instead of quoting the explicit form of the Legendre polynomials, we equivalently state the Rodrigues formula for these polynomials

\displaystyle P_{l}(\mu)=\frac{1}{2^{l}l!}\frac{d^{l}}{d\mu^{l}}(\mu^{2}-1)^{l}. (6)

However, the solutions to Eq.(5) are the associated Legendre functions, not polynomials. Therefore, we use the following to determine the associated Legendre functions from the Legendre polynomials:

\displaystyle {\Theta_{l}}^{m}(\mu) = (1-\mu^{2})^{\frac{|m|}{2}}\frac{d^{|m|}}{d\mu^{|m|}}P_{l}(\mu). (7)

For real solutions |m|\leq l and for complex solutions -l \leq m \leq l, thus we can write the solutions as series whose indices run from 0 to l for the real-values and from -l to l for the complex-valued solutions.

Now we turn to the radial part which upon substitution of -l(l+1) in place of the angular component, we get

\displaystyle \frac{d}{dr}\bigg\{r^{2}\frac{dR}{dr}\bigg\}-l(l+1)R=0,

and if we evaluate the derivatives in the first term we get an Euler equation of the form

\displaystyle r^{2}\frac{d^{2}R}{dr^{2}}+2r\frac{dR}{dr}-l(l+1)R=0. (8)

Let us assume that the solution can be represented as

\displaystyle R(r)=\sum_{j=0}^{\infty}a_{j}r^{j}.

Taking the necessary derivatives and substituting into Eq.(8) gives

\displaystyle r^{2}\sum_{j=0}^{\infty}j(j-1)a_{j}r^{j-2}+\sum_{j=0}^{\infty}2rja_{j}r^{j-1}-l(l+1)\sum_{j=0}^{\infty}a_{j}r^{j}=0. (9)

Simplifying the powers of r we find that

\displaystyle \sum_{j=0}^{\infty}[j(j-1)+2j-l(l+1)]a_{j}r^{j}=0.

Now, since \displaystyle \sum_{j=0}^{\infty}a_{j}r^{j}\neq 0, this means that

\displaystyle j(j-1)+2j-l(l+1)=0. (10)

We can simplify this further to get

\displaystyle j^{2}-l^{2}+j-l=0.

Factoring out (j-l), we arrive at

\displaystyle (j-l)((j+l)+1)=0 \implies j=l ; j=-l-1.

Since this must be true for all values of l, the solution therefore becomes

\displaystyle R(r)= \sum_{l=0}^{\infty}Br^{l}+Cr^{-l-1}. (11)

Thus, we can write the solution to Laplace’s equation as

\displaystyle \psi(r,\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=0}^{\infty}[Br^{l}+Cr^{-l-1}]\Theta_{l}^{m}(\mu)A\exp({im\phi}),  (12)

for real solutions, and

\displaystyle \psi(r,\theta,\phi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}[Br^{l}+Cr^{-l-1}]\Theta_{l}^{m}(\mu)A\exp({im\phi}), (13)

for complex solutions.