This section introduces the spaces of approximation and the time marching considered. The weak formulations of each equations (Navier-Stokes, temperature and Maxwell) are also described. Additionnal information are given for specific feature of the code such as the entropy viscosity method, multiphase flow problem, etc.

Fourier/Finite element representation

The SFEMaNS code uses a hybrid Fourier/Finite element formulation. The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode independently, modulo nonlinear terms that are made explicit. The variables are then approximated on a meridian section of the domain with a finite element method.

The numerical approximation of a function $f$ is written in the following generic form:

\begin{align*} f(r,\theta,z,t)=f_h^{0,\cos}(r,z,t) + \sum_{m=1}^M f_h^{m,\cos} \cos(m\theta) + f_h^{m,\sin} \sin(m\theta), \end{align*}

with $(r,\theta,z)$ the cylindrical coordinates, $t$ the time and $M$ the number of Fourier modes considered. The unknown $f_h^{m,\cos}$ and $f_h^{m,\sin}$ can be approximated independtly modulo the computation of nonlinear terms. Introducing the functions $\cos_m = \cos(m\theta)$, $\sin_m = \sin(m\theta)$ and a basis functions $(\phi_j)_{j \in J}$ of the finite element space of the meridian section results in $(\phi_j \cos_m)_{j\in J, m \in [|0,M|]} \cup (\phi_j \sin_m)_{j\in J, m \in [|1,M|]}$ being a basis of the space of approximation.

Space of approximation

We set the number of Fourier mode to $M$. We define the meridian sections $\Omega_{c,f}^{2D}$, $\Omega_{T}^{2D}$, $\Omega_{c}^{2D}$, $\Omega_v^{2D}$ and $\Omega^{2D}$ of the domains $\Omega_{c,f}$, $\Omega_{T}$, $\Omega_{c}$, $\Omega_v$ and $\Omega$ introduced in the section SFEMaNS presentation. We consider $\left\{ \mathcal{T}_h \right\}_{h > 0}$ a family of meshes of the meridian plane $\Omega^{2D}$ composed of disjoint triangular cells $K$ of diameters at most $h$. For a given $h>0$, we assume that the mesh $\mathcal{T}_h$ can be restricted to each sub domain of interests. These sub-meshes are denoted $\mathcal{T}_h^{c,f}$, $\mathcal{T}_h^{T}$, $\mathcal{T}_h^{c}$ and $\mathcal{T}_h^{v}$. The approximation of the solutions of the Navier-Stokes, heat and Maxwell equations either involve $\mathbb{P}_1$ or $\mathbb{P}_2$ Lagrange finite elements. The following defines the space of approximation used for each dependent variable.

The velocity field $\bu$ and the pressure $p$ are respectively approximated in the following spaces:
\begin{align*} \textbf{V}_{h}^\bu := \left\{ \textbf{v} =\sum\limits_{k=-M}^M \textbf{v}_h^k (r,z) e^{ik \theta} ; \textbf{v}_h^k \in \textbf{V}_{h}^{\bu,2D} \text{, } \overline{\textbf{v}_h^k}=\textbf{v}_h^{-k} \text{, } -M \leq k \leq M \right\} ,\\ S_{h}^p := \left\{ q_h= \sum\limits_{k=-M}^M q_h^k (r,z) e^{i k \theta} ; q_h^k \in S_{h}^{p,2D} \text{, } \overline{q_h^k}=q_h^{-k} \text{, } -M \leq k \leq M \right\}, \end{align*}
where we introduce the following finite element space:
\begin{align*} \textbf{V}_{h}^{\bu,2D} : = \left\{ \textbf{v}_h \in C^0(\overline{\Omega_{c,f}^{2D}}) ; \textbf{v}_h|_K \in \mathbb{P}_2^6 \text{ } \forall K \in \mathcal{T}_h^{c,f} \right\} ,\\ S_{h}^{p,2D} : = \left\{ q_h \in C^0(\overline{\Omega_{c,f}^{2D}}) ; q_h|_K \in \mathbb{P}_1^2 \text{ } \forall K \in \mathcal{T}_h^{c,f} \right\} . \end{align*}
We also introduce the subspace $\textbf{V}_{h,0}^{2D}$ of $\textbf{V}_{h}^\bu$ composed of vector field that are zero on the boundaries of $\Omega_{c,f}$.
The temperature $T$ is approximated in the following space:
\begin{align*} S_h^T := \left\{ q_h= \sum\limits_{k=-M}^M q_h^k (r,z) e^{i k \theta} ; q_h^k \in S_{h}^{T,2D} \text{, } \overline{q_h^k}=q_h^{-k} \text{, } -M \leq k \leq M \right\}, \end{align*}
where we introduce the following finite element space:
\begin{align*} S_{h}^{T,2D} : = \left\{ q_h \in C^0(\overline{\Omega_{T}^{2D}}) ; q_h|_K \in \mathbb{P}_2^2 \text{ } \forall K \in \mathcal{T}_h^{T} \right\} . \end{align*}
The magnetic field $\bH^c$ and the scalar potentiel $\phi$ are respectively approximated in the following spaces:
\begin{align*} \textbf{V}_h^{\bH^c} := \left\{ \textbf{b}= \sum\limits_{k=-M}^M \textbf{b}_h^k (r,z) e^{ik \theta} ; \textbf{b}_h^k \in \textbf{V}_h^{\bH^c,2D}, \; -M \leq k \leq M \right\}, \\ S_h^{\phi} := \left\{ \varphi= \sum\limits_{k=-M}^M \varphi_h^k (r,z) e^{ik \theta} ; \varphi_h^k \in S_{h}^{\phi,2D}, \; -M \leq k \leq M \right\}, \end{align*}
where we introduce the following finite element spaces:
\begin{align*} \textbf{V}_{h}^{\bH^c, 2D} : = \left\{ \textbf{b}_h \in C^0(\overline{\Omega_{c}^{2D}}); \textbf{b}_h|_K \in \mathbb{P}_{l_\bH}^6 \text{ } \forall K \in \mathcal{T}_h^{c} \right\} ,\\ S_{h}^{\phi, 2D} : = \left\{ \varphi_h \in C^0(\overline{\Omega_{v}^{2D}}) ; \varphi_h|_K \in \mathbb{P}_{l_\phi}^2 \text{ } \forall K \in \mathcal{T}_h^v , \right\} \end{align*}
with $l_\bH\in\{1,2\}$ and $l_\phi\in\{1,2\}$ such that $l_\phi \geq l_\bH$.

Time marching

To present the time marching, we introduce a time step $\tau$ and denote by $f^n$ the approximation of a function at the time $t_n=n \tau$. When approximating the Navier-Stokes, heat and Maxwell equations, the time marching can be summarized by the four following steps:

Initialization of the temperature $T^0, T^1$, velocity fields $\bu^0, \bu^1$ , the dynamical pressure $p^0, p^1$ , the magnetic field $\bH^0, \bH^1$ and the scalar potential $\phi^0, \phi^1$.
Approximation of $T^{n+1}$ after the nonlinear terms are computed with extroplation involving $T^n, T^{n-1}, \bu^n$ and $\bu^{n-1}$.
Approximation of $\bu^{n+1}$ and $p^{n+1}$ after the nonlinear terms are computed with extroplation involving $\bu^n, \bu^{n-1}, \bH^n$ and $\bH^{n-1}$.
Approximation of $\bH^{n+1}$ and $\phi^{n+1}$ after the nonlinear terms are computed with extroplation involving $\bu^{n+1}, \bH^n$ and $\bH^{n-1}$.

Weak formulation and extensions

This section introduces the weak formulations implemented in SFEMaNS and additional features/extensions of the code. The notations introduced previously, such as the domain of approximation for each equations or the time step $\tau$, are unchanged.

Hydrodynamic setting

Approximation of the Navier-Stokes equations

The approximation of the Navier-Stokes equations is based on a rotational form of the prediction-correction projection method detailed in Guermond and Shen (2004). As the code SFEMaNS approximates the predicted velocity, a penalty method of the divergence of the velocity field is also implemented.

The method proceeds as follows:

Initialization of the velocity field, the pressure and the pressure increments.
For $n\geq0$ let $\bu^{n+1}$, that matches the Dirichlet boundary conditions of the problem, be the solutions of the following formulation for all $\textbf{v}$ in $\textbf{V}_{h,0}^\bu$:
\begin{equation} \label{eq:SFEMaNS_weak_from_NS_1} \int_{\Omega_{c,f}} \frac{3}{2 \tau} \textbf{u}^{n+1} \cdot \textbf{v} + \frac{2}{\Re} \varepsilon(\textbf{u}^{n+1}) : \GRAD \textbf{v} + \frac{\text{c}_\text{div}}{\Re} \DIV\bu^{n+1} \DIV \bv= - \int_{\Omega_{c,f}} ( \frac{-4 \textbf{u}^n + \textbf{u}^{n-1}}{2 \tau} + \GRAD ( p^n +\frac{4\psi^n - \psi^{n-1}}{3} ) ) \cdot \textbf{v} \\ + \int_{\Omega_{c,f}} ( \textbf{f}^{n+1} - (\ROT \textbf{u}^{*,n+1} ) \times \textbf{u}^{*,n+1} ) \cdot \textbf{v} , \end{equation}
where $\text{c}_\text{div}$ is a penalty coefficient, $\textbf{u}^{*,n+1}=2\textbf{u}^n-\textbf{u}^{n-1}$ and $\varepsilon(\textbf{u})=\GRAD^s \textbf{u} = \frac{1}{2} (\GRAD \textbf{u} +(\GRAD \textbf{u})^{\sf T})$. We remind that $\Re$ is the kinetic Reynolds number and $\textbf{f}$ a source term.
Computation of $\psi^{n+1}$ and $\delta^{n+1}$ solutions in $S_h^p$ of:
\begin{equation} \label{eq:SFEMaNS_weak_from_NS_2} \int_{\Omega_{c,f}} \GRAD \psi^{n+1} \cdot \GRAD q = \frac{3}{2 \tau} \int_{\Omega_{c,f}} \textbf{u}^{n+1} \cdot \GRAD q, \end{equation}

\begin{equation} \label{eq:SFEMaNS_weak_from_NS_3} \int_{\Omega_{c,f}} q \delta^{n+1} = \int_{\Omega_{c,f}} q \DIV \textbf{u} ^{n+1}, \end{equation}
for all $q$ in $S_h^p$.
The pressure is updated as follows:
\begin{equation} \label{eq:SFEMaNS_weak_from_NS_4} p^{n+1} = p^n + \psi^{n+1} - \frac{2}{\Re} \delta^{n+1} - \frac{\text{c}_\text{div}}{\Re}\delta^{n+1}. \end{equation}

Remark: The update of the pressure slighlty differs from Guermond and Shen (2004). It is due to the use of the strain rate tensor $\varepsilon(\textbf{u}) $ in the diffusion operator and of a penalty method of the velocity field divergence.

Entropy viscosity for under resolved computation

Hydrodynamic problems with large kinetic Reynolds number introduce extremely complex flows. Approximating all of the dynamics's scales of such problems is not always possible with present computational ressources. To address this problem, a nonlinear stabilization method called entropy viscosity is implemented in SFEMaNS. This method has been introduced by Guermond et al. (2011). It consists in introducing an artifical viscosity, denoted $\nu_{E}$, that is taken proportional to the default of equilibrium of an energy equation.

This implementation of this method in SFEMaNS can be summarized in the three following steps:

Define the residual of the Navier-Stokes at the time $t_n$ as follows:
\begin{equation} \textbf{Res}_\text{NS}^n:= \frac{\bu^n-\bu^{n-2}}{ 2 \tau} -\frac{2}{\Re} \DIV (\varepsilon(\textbf{u}^{n-1})) + \ROT (\bu^{*,n-1}) \times \bu^{*,n-1} + \GRAD p^{n-1} -\textbf{f}^{n-1} , \end{equation}
Compute the entropy viscosity on each mesh cell K as follows:
\begin{equation} \label{eq:SFEMaNS_NS_entropy_viscosity} \nu_{E|K}^{n}:=\min\left(c_\text{max} h \|\bu^{n-1}\|_{\bL^\infty(K)}, c_\text{e} h^2 \frac{\|\textbf{Res}_\text{NS}^n \cdot \bu^{n-1}\|_{\bL^\infty(K)}}{\|\bu^{n-1}\|_{\bL^\infty(K)}^2}\right), \end{equation}
with $h$ the local mesh size of the cell K, $c_\text{max}=\frac{1}{2}$ for ${\mathbb P}_1$ finite elements and $c_{\max}=\frac{1}{8}$ for ${\mathbb P}_2$ finite elements. The coefficient $c_\text{e}$ is a tunable constant in the interval $(0,1)$. It is generally set to one.
When approximating $\bu^{n+1}$, the term $-\DIV(\nu_{E}^n \GRAD \bu^n)$ is added in the left handside of the Navier-Stokes equations.

Thus defined, the entropy viscosity is expected to be smaller than the consistency error in the smooth regions. In regions with large gradients, the entropy viscosity switches to the first order viscosity $\nu_{\max|K}^n:=c_\text{max} h_K \|\bu^{n-1}\|_{\bL^\infty(K)}$. Note that $\nu_\max^n$ corresponds to the artifical viscosity induces by first order up-wind scheme in the finite difference and finite volume litterature.

Remark: To facilitate the explicit treatment of the entropy viscosity, the following term can be added in the left handside of the Navier-Stokes equations:

\begin{equation} \label{eq:SFEMaNS_NS_LES_c1} - \DIV( c_1 h \GRAD (\bu^{n+1}-\bu^{*,n+1})). \end{equation}

with $h$ the local mesh size and $c_1$ is a tunable constant. The coefficient $c_1$ should be of the same order of $c_\text{max} \|\bu\|_{\bL^\infty(\Omega_{c,f})}$.

Extension to non axisymmetric geometry

A penalty method of Pasquetti et al. (2008) is implemented so the code SFEMaNS can report of the presence of non axisymmetric solid domain in $\Omega_{c,f}$. Such solid domains can either be driving the fluid or represents an obtacle to the fluid motion when their velocity is zero. The domain $\Omega_{c,f}$, where the Navier-Stokes equations are approximated, is splitted into a fluid domain $\Omega_\text{fluid}$ and a solid domain $\Omega_\text{obs}$. These sub domains can be non axisymmetric and time dependent. The penalty method introduces a penalty function $\chi$. It is used to force the velocity field approximated by the Navier-Stokes equations to match the given velocity field of the solid in $\Omega_\text{obs}$. This penalty function is defined as follows:

\begin{equation} \label{eq:SFEMaNS_NS_penal_1} \chi = \left\{ \begin{array}{c} 1 \text{ in } \Omega_\text{fluid}, \\ 0 \text{ in } \Omega_\text{obs}. \end{array} \right. \end{equation}

The velocity field is updated as follows:

\begin{equation} \label{eq:SFEMaNS_NS_penal_2} \frac{3\bu^{n+1}}{2\tau} - \frac{2}{\Re} \DIV (\varepsilon(\textbf{u}^{n+1})) - \frac{\text{c}_\text{div}}{\Re} \GRAD (\DIV\bu^{n+1}) = -\GRAD p^{n} + \chi^{n+1} \left(\frac{4\bu^n - \bu^{n-1}}{2\tau} - \GRAD( \frac{4\psi^n-\psi^{n-1}}{3}) \right) \\ + \chi^{n+1} \left( - ( \ROT \bu^{*,n+1} ) \times\bu^{*,n+1} + \textbf{f}^{n+1} \right) + (1 - \chi^{n+1}) \frac{3\bu^{n+1}_\text{obst}}{2\tau}, \end{equation}

with $\bu_\text{obs}$ the given velocity of the solid obstacle. As the subdomains $\Omega_\text{fluid}$ and $\Omega_\text{obs}$ can be time dependent so is the penalty function $\chi$. Note that the original scheme is recovered where $\chi=1$.

Remark: The update of the pressure is not modified.

Extension to multiphase flow problem

The code SFEMaNS can approximate two phase flow problems. The governing equations can be written as follows:

\begin{equation} \label{eq:SFEMaNS_NS_multiphase_1} \partial_t \rho + \DIV( \textbf{m}) = 0, \end{equation}

\begin{equation} \label{eq:SFEMaNS_NS_multiphase_2} \partial_t(\textbf{m}) + \DIV(\textbf{m}{\otimes}\bu) - \frac{2}{\Re} \DIV(\eta \varepsilon(\bu)) = -\GRAD p + \textbf{f}, \end{equation}

\begin{equation} \label{eq:SFEMaNS_NS_multiphase_3} \DIV \bu = 0, \end{equation}

where $\rho$ is the density, $\eta$ the dynamical viscosity, $\bm=\rho\bu$ the momentum, $\textbf{f}$ a forcing term and $\varepsilon(\bu)=\GRAD^s \bu = \frac12 (\GRAD \bu +(\GRAD \bu)^\sf{T})$. The densities, respestively dynamical viscosities, of the two fluids are denoted $\rho_0$ and $\rho_1$, respectively $\eta_0$ and $\eta_1$.

The approximation method is based on the following ideas.

Use of a level set method to follow the interface evolution. The method consists of approximating $\varphi$ that takes value in $[0,1]$ solution of:
\begin{align*} \partial_t \varphi + \bu \cdot \GRAD \varphi=0. \end{align*}
The level set is equal to 0 in a fluid and 1 in the other fluid. The interface between the fluid is represented by $\varphi^{-1}(\{1/2\})$.
Use the momentum as dependent variable for the Navier-Stokes equations. The mass matrix becomes time independent and can be treated with pseudo-spectral method.
Rewritte the diffusive term $- \frac{2}{\Re} \DIV(\eta \varepsilon(\bu)) $ as follows:
\begin{align*} - \frac{2}{\Re} \DIV(\eta \varepsilon(\bu)) = - \frac{2}{\Re} \DIV(\overline{\nu} \varepsilon(\bm)) - \left( \frac{2}{\Re} \DIV(\eta \varepsilon(\bu)) - \frac{2}{\Re} \DIV(\overline{\nu} \varepsilon(\bm)) \right) \end{align*}
with $\overline{\nu}$ a constant satisfying $\overline{\nu}\geq \frac{\eta}{\rho}$. The first term is made implicit while the second is treated explicitly. The resulting stiffness matrix is time independent and does not involve nonlinearity.
The level set and Navier-Stokes equations are stabilized with the same entropy viscosity. For each mesh cell $K$ and each time iteration $n$, the entropy viscosity $\nu_E$ is defined as follows:
\begin{align*} \nu_{E|K}^{n}:=\min\left(c_\text{max} h \|\bu^{n-1}\|_{\bL^\infty(K)}, c_\text{e} h^2 \frac{\ \max( |\textbf{Res}_\text{NS}^n \cdot \bu^{n-1}\|_{\bL^\infty(K)}, |\text{Res}_\rho^n \|\bu{n-1}\|^2| }{\|\bu^{n-1}\|_{\bL^\infty(K)}\|\bm^{n-1}\|_{\bL^\infty(K)}} \right), \end{align*}
where
\begin{align*} \textbf{Res}_\text{NS}^n= \frac{\bm^n-\bm^{n-2}}{ 2 \tau} -\frac{1}{\Re} \DIV (\eta^{n-1}\epsilon(\bu^{n-1})) + \DIV(\bm^n{\otimes}\bu^n) + \GRAD p^{n-1} -\textbf{f}^{n-1} , \end{align*}
and
\begin{align*} \text{Res}_\rho^n= \frac{\bm^n-\bm^{n-2}}{ 2 \tau} + \DIV (\bm^{m-1}). \end{align*}
To facilitate the explicit treatment of the entropy viscosity, the term $- \DIV( c_1 h \GRAD (\bu^{n+1}-\bu^{n}))$, respectively $-\DIV( c_1 h \GRAD (\varphi^{n+1}-\varphi^n))$, can be added in the left handside of the Navier-Stokes, respectively of level set equation.
A compression term that allows the level set to not get flatten over time iteration is added. It consists of adding the following term in the right handside of the level set equation:
\begin{align*} - \DIV \left(c_\text{comp}\nu_E h^{-1} \varphi(1-\varphi)\frac{\GRAD\varphi}{\|\varphi\|}\right). \end{align*}
The coefficient $c_\text{comp}$ a tunable constant in $[0,1]$. We generally set $c_\text{comp}=1$.

After initialization, the first time order algorithm (BDF1) proceeds as follows:

Compute $\varphi^{n+1}$ solution of
\begin{align} \frac{\varphi^{n+1}-\varphi^n}{\tau} = - \bu^n \cdot \GRAD \varphi^n + \DIV \left( \nu_E^n\GRAD \varphi^n - c_\text{comp} \nu_E^n h^{-1} \varphi^n(1-\varphi^n)\frac{\GRAD\varphi^n}{\|\varphi^n\|} \right). \end{align}
Reconstruct $\rho^{n+1}$ and $\eta^{n+1}$ as follows:
\begin{align*} \rho^{n+1} = \rho_0 + (\rho_1 - \rho_0) F(\varphi^{n+1}), \qquad \eta = \eta_0 + (\eta_1 - \eta_0) F(\varphi^{n+1}), \end{align*}
where $F$ is either equal to the identity, $F(\varphi)=\varphi$, or a piecewise ponylomial function defined by:
\begin{align*} F(\varphi) = \begin{cases} 0 & \text{if $\varphi - 0.5\le -c_{\text{reg}}$}, \\ \frac12 \left(1+\frac{(\varphi-0.5)((\varphi-0.5)^2 - 3 c_{\text{reg}}^2)}{-2c^3_{\text{reg}}}\right) & \text{if $|\varphi - 0.5| \le c_{\text{reg}}$}, \\ 1 & \text{if $c_{\text{reg}} \le \varphi - 0.5$}. \end{cases} \end{align*}
The tunable coefficient $c_\text{reg}$ lives in $[0,0.5]$. We generally set $c_\text{reg}=0.5$.
Compute $\bm^{n+1}$ solution of:
\begin{align} \frac{\bm^{n+1}-\bm^n}{\tau} - \frac{2\overline{\nu}}{\Re}\DIV(\epsilon(\bm^{n+1})-\epsilon(\bm^n)) = \frac{2}{\Re}\DIV( \eta^n\epsilon(\bu^n)) - \DIV(\bm^n\times\bu^n) - \GRAD(p^n+\phi^n) +\textbf{f}^{n+1}. \end{align}
Update the pressure as follows:
1. Solve $\phi^{n+1}$ solution of
  \begin{align*} - \LAP \phi^{n+1} = \frac{\rho_\text{min}}{\tau} \DIV \bu^{n+1}, \end{align*}
  with $\rho_\text{min}=\min(\rho_0,\rho_1)$ and $\bu^{n+1}=\frac{1}{\rho^{n+1}}\bm^{n+1}$.
2. Set $p^{n+1}=p^n + \phi^{n+1} - \frac{\eta_\text{min}}{\Re} \DIV \bu^{n+1}$ with $\eta_\text{min}=\min(\eta_0,\eta_1)$.

Remarks:

This method can be used to approximate problems with a stratification or an inclusion of $n\geq 3$ fluids. One level set is approximated per interface between two fluids. The fluids properties are reconstructed with recursive convex combinations.
MHD multiphase problems with variable electrical conductivity between the fluids can also be considered. The electrical conductivity in the fluid is reconstructed with the level set the same way the density and the dynamical viscosity are. The magnetic field $\bH^{n+1}$ is updated as follows:
\begin{equation*} \frac{3\bH^{n+1}-4\bH^n+\bH^{n-1}}{2\tau} + \ROT \left( \frac{1}{\overline{\sigma}\Rm} \ROT ( \bH^{n+1}-\bH^{*,n+1}) \right) = - \ROT\left( \frac{1}{\sigma\Rm} \ROT \bH^{*,n+1} \right) \\ + \ROT (\bu^{n+1}\times \mu^c \bH^{*,n+1}) + \ROT \left( \frac{1}{\sigma\Rm} \textbf{j}^{n+1} \right) \end{equation*}
with $\bH^{*,n+1}=2\bH^{n+1}-\bH^n$ and $\overline{\sigma}$ a function depending of the radial and vertical coordinates $(r,z)$ such that $\overline{\sigma}(r,z)\leq \sigma(r,\theta,z,t)$ for all $(r,\theta,z,t)$ considered.

Heat equation's weak formulation

The heat equations is approximated as follows.

Initialization of the temperature.
For all $n\geq0$ let $T^{n+1}$, that matches the Dirichlet boundary conditions of the problem, be the solution of the following formulation for all $v\in S_h^T$:
\begin{equation} \label{eq:SFEMaNS_weak_form_temp} \int_{\Omega_T} \frac{3 C }{2 \tau}T^{n+1} v + \lambda \GRAD T^{n+1} \cdot \GRAD v = - \int_{\Omega_T} \left( C \frac{4 T^n -T^{n-1}}{2 \tau} - \DIV (T^{*,n+1} \bu^{*,n+1}) + f_T^{n+1}\right) v, \end{equation}
where $T^{*,n+1}=2 T^n - T^{n-1}$. We remind that $C$ is the volumetric heat capacity, $\lambda$ the thermal conductivty and $f_T$ a source term.

Magnetic setting

The code SFEMaNS uses $\bH^1$ conforming Lagrange finite element to approximate the magnetic field. As a consequence, the zero divergence condition on the magnetic field cannot be enforced by standard penalty technique for non-smooth and non-convex domains. To overcome this obstacle, a method inspired of Bonito and Guermond (2011) has been implemented. This method consists of introducing a magnetic pressure denoted $ p_\text{m}$ and proceeds as follows.

Add the term $-\mu^c \GRAD p_\text{m}^c$ in the right handside of the magnetic field $\bH^c$ equation with $p_\text{m}^c$ the solution in $\Omega^c$ of:
\begin{align*} - \DIV( h_\text{loc}^{2(1-\alpha)} \GRAD p_m^{c,n+1} ) &= - \DIV( \mu^c \GRAD \bH^{c,n+1}) , \\ p_m^{c,n+1}|_{\partial \Omega_c} &= 0, \end{align*}
where $h$ is the local mesh size and $\alpha$ a constant parameter in $[0.6,0.8]$. It is set to 0.6 in the following.
Add the term $ -\DIV(\mu^v \GRAD p_\text{m}^v)$ in the right handside of the scalar potential $\phi$ equation with $p_\text{m}^v$ the solution in $\Omega^v$ of:
\begin{align*} \LAP p_m^{v,n+1} = \LAP \phi^{n+1}, \\ \GRAD p_m^{v,n+1} \cdot \textbf{n} |_{\partial \Omega_v} = 0. \end{align*}

We note that the magnetic pressure can be eliminated from the equation of the scalar potential $\phi$. We refer to Guermond et al. (2011) for more details. The approximation space used to approximate $ p_\text{m}^c$ is the following:

\begin{align*} S_h^{p_\text{m}^c} := \left\{ \varphi= \sum\limits_{k=-M}^M \varphi_h^k (r,z) e^{ik \theta} ; \varphi_h^k \in S_{h}^{p_\text{m}^c,2D}, \; -M \leq k \leq M \right\}, \end{align*}

where we introduce the following finite element space:

\begin{align*} S_{h}^{p_\text{m}^c, 2D} : = \left\{ \varphi_h \in C^0(\overline{\Omega_{c}^{2D}}) ; \varphi_h|_K \in \mathbb{P}_1^2 \text{ } \forall K \in \mathcal{T}_h^c , \right\}. \end{align*}

In addition, an interior penalty method is used to enforce the continuity conditions across the interfaces $\Sigma_\mu$ and $\Sigma$ introduced in the section SFEMaNS presentation. We refer to Guermond et al. (2007) for more details.

In the following, we assume that the potential scalar $\phi$ matches the Dirichlet boundary conditions of the problem on the set $\Gamma_{v,D}$ that is included in $\partial \Omega_v \setminus \Sigma$.

Approximation of the Maxwell equations with H

The Maxwell equations are approximated as follows:

Initialization of the magnetic field $\bH^c$, the scalar potential $\phi$ and the magnetic pressure $p_\text{m}^c$.
For all $n\geq 1$, computation of $(\bH^{c,n+1},\phi^{n+1},p_\text{m}^{c,n+1})$ solutions of the following formulation for all $b\in \bV_h^{\bH^c} $, $\varphi\in S_h^{\phi}$ and $q\in S_h^{p_\text{m}^c} $.

\begin{align*} & \int_{\Omega_c}\mu^c \frac{D\bH^{c,n+1}}{\Delta t}\SCAL \bb +\int_{\Omega_c} \frac{1}{\sigma R_m} \ROT \bH ^{c,n+1}\cdot \ROT \bb +\int_{\Omega_v} \muv\frac{\GRAD D\phi^{n+1}}{\Delta t}\SCAL \GRAD\varphi \\ & + \frac{\beta_1}{R_m} \int_{\Omega_c} \mu^c\GRAD p_\text{m}^{c,n+1}\SCAL\bb + \frac{\beta_1}{R_m}\int_{\Omega_c} \frac{1}{\sigma_\text{min} \mu_\text{min}^2} (\frac{h}{D_{\Omega_c}})^{2\alpha} \DIV (\mu^c \bH^{c,n+1} )\DIV (\mu^c \bb) \\ & - \frac{\beta_1}{R_m} \int_{\Omega_c} \mu^c\bH^{c,n+1}\SCAL\GRAD q + \frac{\beta_1}{R_m} \int_{\Omega_c} \sigma_\text{min} \mu_\text{min}^2 D_{\Omega_c}^2 (\frac{h}{D_{\Omega_c}})^{2(1-\alpha)} \GRAD p_\text{m}^{c,n+1}\SCAL \GRAD q \\ & +\int_{\Sigma_{\mu}} \left \{ \frac{1}{\sigma R_m} \ROT {\bH ^{c,n+1}} \right \} \cdot \left ( { \bb_1}\times \bn_1^c + { \bb_2}\times \bn_2^c\right )\\ & +\frac{\beta_3}{R_m} \int_{\Sigma_{\mu}} \frac{1}{\sigma_\text{min}D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{-1} \left( { \bH_1^{c,n+1}}\times \bn_1^c + {\bH_2^{c,n+1}}\times \bn_2^c\right ) \SCAL \left ( { \bb_1}\times \bn_1^c + { \bb_2}\times \bn_2^c\right )\\ & +\frac{\beta_1}{R_m} \int_{\Sigma_{\mu}} \frac{1}{\sigma_\text{min} \mu_\text{min}^2 D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{2\alpha-1} \left({ \mu^c_1\bH_1^{c,n+1}}\cdot \bn_1^c + {\mu^c_2 \bH_2^{c,n+1}}\cdot \bn_2^c\right ) \SCAL \left ( {\mu^c_1}{ \bb_1}\cdot \bn_1^c + {\mu^c_2}{ \bb_2}\cdot \bn_2^c\right )\\ & +\int_{\Sigma} \frac{1}{\sigma R_m} \ROT {\bH ^{c,n+1}} \cdot \left( { \bb }\times \bn^c + \nabla \varphi ^{n+1}\times \bn^v\right) \\ & + \frac{\beta_2}{R_m} \int_\Sigma \frac{1}{\sigma_\text{min} D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{-1} \left( {\bH^{c,n+1}}\CROSS \bn_1^c + {\GRAD \phi^{n+1}}\CROSS \bn_2^c\right ) \SCAL (\bb\CROSS \bnc + \GRAD\varphi\CROSS \bnv)\\ & + \frac{\beta_1}{R_m} \int_\Sigma \frac{1}{\sigma_\text{min} \mu_\text{min}^2 D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{2\alpha-1} \left( { \mu^c\bH ^{c,n+1}}\cdot \bn_1^c + {{\mu^v}\GRAD \phi^{n+1}}\cdot \bn_2^c\right ) \SCAL ({\mu^c}\bb\cdot \bnc + {\mu^v}\ \GRAD\varphi \cdot \bnv)\\ & + \int _{\Gamma_c} \frac{1}{\sigma R_m} \ROT \bH ^{c,n+1} \cdot ( \bb \CROSS \bnc) + \frac{\beta_3}{R_m} \left( \int_{\Gamma_c} \frac{1}{\sigma_\text{min} D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{-1} \left( { \bH^{c,n+1}}\CROSS \bn^c \right ) \SCAL (\bb\CROSS \bnc) \right )\\ & =\\ & \int_{\Omega_c} \left( \frac{1}{\sigma R_m}\bj^s + \bu^{n+1} \times \mu^c \bH^{*,n+1} \right ) \cdot \ROT \bb + \int _{\Sigma_{\mu}} \left \{ \frac{1}{\sigma R_m}\bj^s + \bu^{n+1} \times \mu^c \bH^{*,n+1} \right \} \cdot \left( { \bb_1}\times \bn_1^c + { \bb_2}\times \bn_2^c\right )\\ & +\int_{\Sigma}\left ( \frac{1}{\sigma R_m} \bj^s + \bu^{n+1} \times \mu^c \bH^{*,n+1} \right)\cdot \left ( { \bb }\times \bn^c + \nabla \varphi \times \bn^v\right) +\int_{\Gamma_c}(\ba \times \bn) \cdot \left ({\bb} \times \bn \right) + \int_{\Gamma_v} (\ba \times \bn) \cdot (\nabla \varphi \times \bn)\\ & + \int_{\Gamma_c} \left ( \frac{1}{\sigma R_m}\bj^s + \bu^{n+1} \times \mu^c \bH^{*,n+1} \right )\cdot ( \bb \CROSS \bnc) +\frac{\beta_3}{R_m} \int_{\Gamma_c} \frac{1}{\sigma_\text{min} D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{-1} \left( {\bH}_\text{bdy}^{c,n+1}\CROSS \bn^c \right) \SCAL (\bb\CROSS \bnc) , \end{align*}

where we set $D\bH^{c,n+1}=\dfrac{3\bH^{c,n+1}-4\bH^{c,n}+\bH^{c,n-1}}{2}$, $D\phi^{n+1}=\dfrac{3\phi^{n+1}-4\phi^{n}+\phi^{n-1}}{2}$, $\bH^{*,n+1}=2\bH^{c,n}-\bH^{c,n-1}$, $\sigma_\text{min} = \min(\sigma)$ and $\mu_\text{min} = \min(\mu)$. $D_{\Omega_c}$ is the diameter of $\Omega_c$. We use the operator $\{\cdot\}$ defined by $\{f\}=\frac{f_1+f_2}{2}$ on the interface $\Sigma_\mu$. The constants $\beta_1, \beta_2$ and $\beta_3$ are penalty coefficients. They are normalized by $(\sigma\Rm)^{-1}$ so their value can be set to one in the data file. The function $\bH_\text{bdy}^{c}$ is a user function used to impose Dirichlet boundary conditions on the surface $\Gamma_c=\partial\Omega_c \setminus \Sigma$. We also introduced $\Gamma_v=\partial \Omega_v \setminus ( \Sigma \cup \Gamma_{\phi,D} )$ and the function $\textbf{a}$ that satisfies:

\begin{align*} \textbf{a} = \textbf{E} \times \textbf{n} \text{ on } \Gamma_c \cup \Gamma_v, \end{align*}

with $\textbf{E}$ the electric field and $\textbf{n}$ the outward normal to the surface considered.

Remarks:

In the code, the term $ \frac{D(\phi^{n+1})}{\Delta t}$ is replaced by $\frac{3\phi^{n+1}}{2 \Delta t} + \frac{\phi^{n+1}}{10}$. It is consistent with the algorithm used. Indeed the initial conditions $\phi^0$ and $\phi^1$ are such that the initial magnetic fields $\GRAD \phi^0$ and $\GRAD \phi^1$ have a zero divergence. The term $ \frac{\phi^{n+1}}{10}$ is added to reduce the error on $\phi^{n+1}$ when using large time steps.
In the conducting domain $\Omega_c$, the Ohm's law states that: $\textbf{E}= \frac{1}{R_m\sigma} (\ROT \bH - \textbf{j}) - \bu \times (\mu \bH) $.
One can show that $\textbf{a}\times \textbf{n} = - \textbf{E}$.

Extension to magnetic permeability variable in time and azimuthal direction

The use of a Fourier decomposition in the azimuthal direction leads us to use the magnetic field $\bB^c=\mu\bH^c$ as dependent variable of the Maxwell equations in the conducting domain. The mass matrix becomes time independent and can be computed with pseudo-spectral methods. To get a time independent stiffness matrix that does not involve nonlinearity, the diffusive term $\ROT \left(\frac{1}{\sigma\Rm} \ROT\frac{\bB^c}{\mu^c} \right)$ is rewritten as follows:

\begin{align*} \ROT \left( \frac{1}{\sigma\Rm} \ROT \frac{\bB^c}{\mu^c} \right) = \ROT \left( \frac{1}{\sigma\Rm \overline{\mu^c}} \ROT\frac{\bB^c}{\mu^c} \right) + \ROT \left( \frac{1}{\sigma\Rm} \ROT ((\frac{1}{\mu^c}-\frac{1}{\overline{\mu^c}})\bB^c) \right) \end{align*}

with $\overline{\mu^c}$ a function depending of the radial and vertical coordinates $(r,z)$ such that $\overline{\mu^c}(r,z)\leq \mu^c(r,\theta,z,t)$ for all $(r,\theta,z,t)$ considered. The first term is then made implicit while the term involving $\frac{1}{\mu^c}-\frac{1}{\overline{\mu^c}}$ is treated explicitly.

Under the previous notations and assuming,

\begin{align*} \overline{\mu^c} &= \mu^c \quad \text{on} \quad \Sigma, \\ \overline{\mu^c} &= \mu^c \quad \text{on} \quad \Sigma_{\mu},\\ \overline{\mu^c} &= \mu^c \quad \text{on} \quad \Gamma_c, \end{align*}

the Maxwell equations are approximated as follows.

Initialization of the magnetic field $\bB^c$, the scalar potential $\phi$ and the magnetic pressure $p_\text{m}^c$.
For all $n\geq 1$, computation of $(\bB^{c,n+1},\phi^{n+1},p_\text{m}^{c,n+1})$ solutions of the following formulation for all $b\in \bV_h^{\bH^c} $, $\varphi\in S_h^{\phi}$ and $q\in S_h^{p_\text{m}^c} $.

\begin{align*} & \int_{\Omega_c}\frac{D\bB^{c,n+1}}{\Delta t}\SCAL \bb + \int _{\Omega_c} \frac{1}{\sigma R_m} \ROT \frac{\bB ^{c,n+1}}{\overline{\mu^c}}\cdot \ROT \bb + \int_{\Omega_v} \mu^v\frac{\GRAD D\phi^{n+1}}{\Delta t}\SCAL \GRAD\varphi \\ & + \frac{\beta_1}{R_m} \int_{\Omega_c} \overline{\mu^c}\GRAD p_\text{m}^{c,n+1}\SCAL\bb + \frac{\beta_1}{R_m} \int_{\Omega_c} \frac{1}{\sigma_\text{min} \mu_\text{min}^2} (\frac{h}{D_{\Omega_c}})^{2\alpha} \overline{\mu^c} \DIV \bB^{c,n+1} \DIV \bb \\ & - \frac{\beta_1}{R_m} \int_{\Omega_c} \bB^{c,n+1}\SCAL\GRAD q + \frac{\beta_1}{R_m} \int_{\Omega_c} \sigma_\text{min} \mu_\text{min}^2 D_{\Omega_c}^2 (\frac{h}{D_{\Omega_c}})^{2(1-\alpha)} \GRAD p_\text{m}^{c,n+1}\SCAL \GRAD q \\ & +\int _{\Sigma_{\mu}} \left\{ \frac{1}{\sigma R_m} \ROT \frac{\bB ^{c,n+1}}{\overline{\mu^c}} \right \} \cdot \left ( { \bb_1}\times \bn_1^c + { \bb_2}\times \bn_2^c\right ) \\ & +\frac{\beta_3}{R_m} \int_{\Sigma_{\mu}} \frac{1}{\sigma_\text{min}D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{-1} \left( \frac{\bB_1^{c,n+1}}{\overline{\mu^c}_1}\times \bn_1^c + \frac{\bB_2^{c,n+1}}{\overline{\mu^c_2}}\times \bn_2^c \right) \SCAL \left ( { \bb_1}\times \bn_1^c + { \bb_2}\times \bn_2^c\right )\\ & +\frac{\beta_1}{R_m} \int_{\Sigma_{\mu}} \frac{1}{\sigma_\text{min} \mu_\text{min}^2 D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{2\alpha-1} \left( {\bB_1^{c,n+1}}\cdot \bn_1^c + {\bB_2^{c,n+1}}\cdot \bn_2^c\right) \SCAL \left( \overline{\mu^c_1}{ \bb_1}\cdot \bn_1^c + \overline{\mu^c_2}{ \bb_2}\cdot \bn_2^c\right )\\ & +\int _{\Sigma} \frac{1}{\sigma R_m} \ROT \frac{\bB ^{c,n+1}}{\overline{\mu^c}} \cdot \left( {\bb }\times \bn^c + \nabla \varphi \times \bn^v\right) \\ & + \frac{\beta_2}{R_m} \int_{\Sigma} \frac{1}{\sigma_\text{min} D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{-1} \left( \frac{\bB^{c,n+1}}{\overline{\mu^c}}\CROSS \bn_1^c + {\GRAD \phi ^{n+1}}\CROSS \bn_2^c\right) \SCAL (\bb\CROSS \bn^c + \GRAD\varphi\CROSS \bn^v)\\ & + \frac{\beta_1}{R_m} \int_{\Sigma} \frac{1}{\sigma_\text{min} \mu_\text{min}^2 D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{2\alpha-1} \left( {\bB ^{c,n+1}}\cdot \bn_1^c + {{\mu^v}\GRAD \phi ^{n+1}} \cdot \bn_2^c\right) \SCAL \left(\overline{{\mu^c}}\bb\cdot \bn^c + {\mu^v} \GRAD\varphi \cdot \bn^v \right ) \\ & + \int_{\Gamma_c} \frac{1}{\sigma R_m} \ROT \frac{\bB ^{c,n+1}}{\overline{\mu^c}} \cdot ( \bb \CROSS \bn^c) + \frac{\beta_3}{R_m} \left( \int_{\Gamma_c} \frac{1}{\sigma_\text{min} D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{-1} \left( \frac{ \bB^{c,n+1}}{\overline{\mu^c}}\CROSS \bn^c \right ) \SCAL (\bb\CROSS \bn^c) \right)\\ & =\\ & \int_{\Omega_c} \frac{1}{\sigma R_m} \ROT (\langle \overline{\mu^c},{\mu^c} \rangle {\bB^{*,n+1}})\cdot \ROT \bb \\ & \int_{\Omega_c} \left( \frac{1}{\sigma R_m}\bj^s + \bu^{n+1} \times \bB^{*,n+1} \right )\cdot \ROT \bb + \int_{\Sigma_{\mu}} \left \{ \frac{1}{\sigma R_m}\bj^s + \bu^{n+1} \times \bB^{*,n+1} \right \} \cdot \left( { \bb_1}\times \bn_1^c + { \bb_2}\times \bn_2^c\right )\\ & +\int_{\Sigma}\left( \frac{1}{\sigma R_m} \bj^s + \bu^{n+1} \times \bB^{*,n+1} \right ) \cdot \left ( { \bb }\times \bn^c + \nabla \varphi \times \bn^v\right) +\int_{\Gamma_c}(\ba \times \bn) \cdot \left ({\bb} \times \bn \right) + \int_{\Gamma_v}(\ba \times \bn) \cdot (\nabla \varphi \times \bn)\\ & + \int_{\Gamma_c} \left ( \frac{1}{\sigma R_m}\bj^s + \bu^{n+1} \times \bB^{*,n+1} \right )\cdot ( \bb \CROSS \bn^c) +\frac{\beta_3}{R_m} \int_{\Gamma_c} \frac{1}{\sigma_\text{min} D_{\Omega_c}} (\frac{h}{D_{\Omega_c}})^{-1} \left({\bH}_\text{bdy}^{c,n+1}\CROSS \bn^c \right ) \SCAL (\bb\CROSS \bn^c), \end{align*}

where we set $\bB^{*,n+1}=2\bB^n-\bB^{n-1}$, $\langle \overline{\mu^c},{\mu^c}\rangle=\frac{1}{\overline{\mu^c}}- \frac{1}{\mu^c}$, $\sigma_\text{min} = \min(\sigma)$ and $\mu_\text{min} = \min(\mu)$. $D_{\Omega_c}$ is the diameter of $\Omega_c$. The function $\bH_\text{bdy}^{c}$ is a user function used to impose Dirichlet boundary conditions on the surface $\Gamma_c=\partial\Omega_c \setminus \Sigma$.

Remarks:
1. The above formulation can also be used to approximate problem where the magnetic permeability depends of the radial and vertical coordinate $(r,z)$ only. In that case, we set $\overline{\mu^c}=\mu^c$ in the whole conductiving domain $\Omega_c$.
2. Same than for the algorithm with $\bH$, the term $ \frac{D(\phi^{n+1})}{\Delta t}$ is replaced by $\frac{3\phi^{n+1}}{2 \Delta t} + \frac{\phi^{n+1}}{10}$.