Theoretical background

Limace.jl uses a Galerkin projection method to discretize a linearized version of the non-linear rotating magnetohydrodynamic (MHD) equations, following works of Bullard and Gellman (1954), Ivers and Phillips (2008) and Gerick and Livermore (2024). For a very detailed account of all steps, the interested reader is referred to Gerick and Livermore (2024). Here, a quick and simplified overview of the mathematical and physical tools used in the model is given.

In non-linear form, the momentum equations of the incompressible fluid and the induction equation, are written as

\[\begin{align} \frac{\partial \mathbf{U}}{\partial t} + \left(\boldsymbol{\nabla}\times\mathbf{U}\right)\times\mathbf{U} + 2\boldsymbol{\Omega}\times\mathbf{U} &= -\frac{1}{\rho}\nabla P + \frac{1}{\rho\mu_0} \left(\boldsymbol{\nabla}\times\mathbf{B}\right)\times\mathbf{B} + \nu\boldsymbol{\nabla}^2\mathbf{U} + \mathbf{F}, \\ \frac{\partial \mathbf{B}}{\partial t} &= \boldsymbol{\nabla}\times\left(\mathbf{U}\times\mathbf{B}\right) + \eta \boldsymbol{\nabla}^2\mathbf{B}, \end{align}\]

with $\mathbf{U}$ the velocity, $\mathbf{B}$ the magnetic field, $\boldsymbol{\Omega}$ the rotation axis, $\rho$ the fluid density, $P$ the reduced hydrodynamic pressure, $\mu_0$ the magnetic permeability of free space, $\nu$ the kinematic viscosity, $\mathbf{F}$ some additional body force, and $\eta$ the magnetic diffusivitiy.

The velocity, magnetic and pressure fields a linearized, so that

\[\begin{align} \mathbf{U}(\mathbf{r},t) & = \mathbf{U}_0(\mathbf{r})+ \mathbf{u}(\mathbf{r}) e^{\lambda t}, \\ \mathbf{B}(\mathbf{r},t) & = \mathbf{B}_0(\mathbf{r})+ \mathbf{b}(\mathbf{r}) e^{\lambda t}, \\ P(\mathbf{r},t) & = P_0(\mathbf{r})+ p(\mathbf{r}) e^{\lambda t}, \end{align}\]

with $\lambda=-\sigma+\mathrm{i}\omega$, with $\sigma$ the damping rate and $\omega$ the frequency of the oscillatory perturbation to the steady background. Removing the steady part and neglecting higher order terms, the linearized MHD equations read

\[\begin{align} \lambda\mathbf{u} =& -\left(\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{U}_0 -\left(\mathbf{U}_0\cdot\boldsymbol{\nabla}\right)\mathbf{u} -2\boldsymbol{\Omega}\times\mathbf{u} - \frac{1}{\rho}\nabla p\\ &+ \frac{1}{\rho\mu_0}\left(\left(\boldsymbol{\nabla}\times\mathbf{b}\right)\times\mathbf{B}_0+\left(\boldsymbol{\nabla}\times\mathbf{B}_0\right)\times\mathbf{b}\right) + \nu \boldsymbol{\nabla}^2\mathbf{u},\nonumber\\ \lambda\mathbf{b} =& \boldsymbol{\nabla}\times\left(\mathbf{U}_0\times\mathbf{b}\right) + \boldsymbol{\nabla}\times\left(\mathbf{u}\times\mathbf{B}_0\right) + \eta \boldsymbol{\nabla}^2\mathbf{b}. \end{align}\]

Galerkin projection

The linearized equations are projected onto trial vectors $\boldsymbol{\xi}_i$, so that

\[f_{ij} = \int \boldsymbol{\xi}_i \cdot \mathbf{f}\left(\mathbf{u}_j,\mathbf{b}_j, \mathbf{U}_0, \mathbf{B}_0\right)\,\mathrm{d}V,\]

where $\mathbf{f}$ is any of the terms in the considered MHD equations and $\boldsymbol{\xi}_i = [\mathbf{u}_i, \mathbf{b}_i]$.

Poloidal and toroidal basis vectors

Due to the divergence free condition on the velocity and magnetic field, i.e. the flow is incompressible and no magnetic monopoles, it is convenient to decompose the fields into poloidal and toroidal components.

\[\begin{align} \mathbf{u} &= \sum_i \alpha_i\mathbf{u}_i = \sum_{l,m,n} \alpha^S_{lmn}\mathbf{S}^\mathbf{u}_{lmn} + \sum_{l,m,n} \alpha^T_{lmn}\mathbf{T}^\mathbf{u}_{lmn},\\ \mathbf{b} &= \sum_i \beta_i\mathbf{b}_i = \sum_{l,m,n} \beta^S_{lmn}\mathbf{S}^\mathbf{b}_{lmn} + \sum_{l,m,n} \beta^T_{lmn}\mathbf{T}^\mathbf{b}_{lmn}, \end{align}\]

with the respective poloidal $\mathbf{S}$ and toroidal $\mathbf{T}$ basis vectors

\[\begin{align} \mathbf{S}_{lmn} & = \boldsymbol{\nabla}\times\boldsymbol{\nabla}\times S_{ln}(r)Y_l^m(\theta,\phi)\mathbf{r}, \\ \mathbf{T}_{lmn} & = \boldsymbol{\nabla}\times T_{ln}(r) Y_l^m(\theta,\phi)\mathbf{r}. \end{align}\]

Here, $Y_l^m$ is the (fully normalized) spherical harmonic of degree $l$ and order $m$. The scalar functions $S_{ln}(r)$ and $T_{ln}(r)$ are defined for each basis (See also Bases) and are constructed to satisfy the desired boundary condition for each $l,n$.

We therefore need to consider all combinations of poloidal and toroidal vector combinations in the projection of the forces. This leads to several long coupling terms, especially for the Lorentz force and induction term. The integrals of these coupling terms over the spherical sourfaces can be reduced to Adam-Gaunt and Elsasser integrals, and it remains to calculate the radial integral. The detailed equations implemented in Limace.jl are given in Gerick and Livermore (2024).