SFEMaNS  version 5.3
Reference documentation for SFEMaNS
Numerical approximation

This section describes the approximation techniques that are used in SFEMaNS. Everything is based on weak formulations.

Fourier/Finite element representation

SFEMaNS uses a hybrid Fourier/Finite element formulation. The numerical approximation of any 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}(r,z,t) \cos(m\theta) + f_h^{m,\sin}(r,z,t) \sin(m\theta), \end{align*}

with \((r,\theta,z)\) the cylindrical coordinates, \(t\) the time and \(M\) the number of Fourier modes considered. The functions \(f_h^{m,\cos}\) and \(f_h^{m,\sin}\) belongs to a finite element space. The main advantage of this decomposition is that the Fourier modes \(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 basis functions \((\phi_j)_{j \in J}\) of the finite element space of the meridian section, the following set \(\{\phi_j \cos_m\}_{j\in J, m \in [|0,M|]} \cup \{\phi_j \sin_m\}_{j\in J, m \in [|1,M|]}\) is a basis of the approximation space we work with.

Approximation spaces

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 equations, that heat equations and Maxwell's equations either involve \(\mathbb{P}_1\) or \(\mathbb{P}_2\) Lagrange finite elements. The following defines the approximation spaces used for each dependent variable.

  1. 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}^{\bu,2D}\) of \(\textbf{V}_{h}^{\bu,2D}\) composed of the vector fields that are zero at the boundary of \(\Omega_{c,f}\).
  2. 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*}

  3. 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 stepping

To present the time stepping, we introduce a time step \(\tau\) and denote by \(f^n\) the approximation of \(f\) at the time \(t_n=n \tau\). When approximating the Navier-Stokes, the heat, and the Maxwell equations, the time stepping can be summarized formulated as follows:

  1. 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\).
  2. 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}\).
  3. 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}\).
  4. 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:

  1. Initialization of the velocity field, the pressure and the pressure increments.
  2. 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.
  3. 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\).
  4. 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:

  1. 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}

  2. 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.
  3. 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.

  1. 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\})\).
  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.
  3. 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.
  4. 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{\rho^n-\rho^{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.
  5. 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:

  1. 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}

  2. 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\).
  3. 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}

  4. 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)\).

Remark: The last two steps of the above algorithm (known as projection method) can be replaced by an artificial compression method that reads:

  1. 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{\alpha}{\overline{\rho}} \GRAD (\DIV(\bm^{n+1}-\bm^n))) = \frac{2}{\Re}\DIV( \eta^n\epsilon(\bu^n)) - \alpha \GRAD (\DIV \bu^n) - \DIV(\bm^n\times\bu^n) - \GRAD(p^n) +\textbf{f}^{n+1}. \end{align}

  2. Update the pressure as follows:

    \begin{align} p^{n+1}= p^n - \alpha \DIV \bu^{n+1} \end{align}

    where $$ is a tunable constant and \(\overline{rho} = \min(\rho) \).

Remarks:

  1. 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.
  2. 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.

  1. Initialization of the temperature.
  2. 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.

  1. 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 \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]\).
  2. 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\). We refer to Guermond et al. (2007) for more details.

Approximation of the Maxwell equations with H

The Maxwell equations are approximated as follows:

  1. Initialization of the magnetic field \(\bH^c\), the scalar potential \(\phi\) and the magnetic pressure \(p_\text{m}^c\).
  2. 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 +\int_{\Omega_v} \muv\GRAD\phi^{n+1}\SCAL \GRAD\varphi - \int_{\partial\Omega_v} \muv\varphi \bn\SCAL \GRAD \phi^{n+1}\\ & + \beta_1\left(\int_{\Omega_c} \mu^c\GRAD p_\text{m}^{c,n+1}\SCAL\bb - \int_{\Omega_c} \mu^c\bH^{c,n+1}\SCAL\GRAD q + \int_{\Omega_c} h^{2(1-\alpha)}\GRAD p_\text{m}^{c,n+1}\SCAL \GRAD q + \int_{\Omega_c} h^{2\alpha}\DIV (\mu^c \bH^{c,n+1} )\DIV (\mu^c \bb)\right)\\ & +\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 )\\ & +\beta_3 \int_{\Sigma_{\mu}} h^{-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 )\\ & +\beta_1 \int_{\Sigma_{\mu}} h^{-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) + \beta_2 \int_\Sigma h^{-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)\\ & + \beta_1 \int_\Sigma h^{-1} \left( { \mu^c\bH ^{c,n+1}}\cdot \bn_1^c + {\GRAD \phi^{n+1}}\cdot \bn_2^c\right ) \SCAL ({\mu^c}\bb\cdot \bnc + \GRAD\varphi \cdot \bnv)\\ & + \int _{\Gamma_c} \frac{1}{\sigma R_m} \ROT \bH ^{c,n+1} \cdot ( \bb \CROSS \bnc) + \beta _3\left( \int_{\Gamma_c} h^{-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) +\beta_3 \int_{\Gamma_c} h^{-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^{c,n+1}=\dfrac{3\phi^{c,n+1}-4\phi^{c,n}+\phi^{c,n-1}}{2}\), \(\bH^{*,n+1}=2\bH^{c,n}-\bH^{c,n-1}\). 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\).

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.

  1. Initialization of the magnetic field \(\bB^c\), the scalar potential \(\phi\) and the magnetic pressure \(p_\text{m}^c\).
  2. 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 + \int_{\Omega_v} \mu^v\GRAD\phi^{n+1}\SCAL \GRAD\varphi - \int_{\partial\Omega_v} \mu^v\varphi \bn\SCAL \GRAD \phi^{n+1} \\ & + \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 + {\GRAD \phi ^{n+1}} \cdot \bn_2^c\right) \SCAL \left(\overline{{\mu^c}}\bb\cdot \bn^c + \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_N}(\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\).

    Remark: 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 conducting domain \(\Omega_c\).