SFEMaNS
version 5.3
Reference documentation for SFEMaNS
|
This section describes the approximation techniques that are used in SFEMaNS. Everything is based on weak formulations.
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.
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.
\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}\).\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*}
\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\).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:
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.
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:
\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.\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\).\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.
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:
\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}
\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.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})}\).
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.
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.
\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\})\).\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.\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.\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:
\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}
\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\).\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}
\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}\).Remark: The last two steps of the above algorithm (known as projection method) can be replaced by an artificial compression method that reads:
\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}
\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:
\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.The heat equations is approximated as follows.
\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.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.
\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]\).\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.
The Maxwell equations are approximated as follows:
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\).
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.
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\).