Fortran file condlim.f90

In this section we describe the subroutines and the functions of the file condlim.f90. They allow to impose initial conditions, boundaries conditions and the source terms of the equations considered. We remind that the code SFEMaNS can approximate the Navier-Stokes equation, the temperature equation and the Maxwell equations for various kind of problems (variable density, solid obstacle in the fluid domain, variable magnetic permeability). A template of the file condlim.f90 is available in the following directory: ($SFEMaNS_DIR)/TEMPLATE.

Settings for the Navier-Stokes equations

General case

The code SFEMaNS can approximate the solutions of the Navier-Stokes equations in an axisymmetric domain \(\Omega_{\bu}\). For a problem with constant density, constant viscosity and without solid obstacle, they are written as follows:

\begin{align*} \partial_t\bu+\left(\ROT\bu\right)\CROSS\bu - \frac{1}{\Re}\LAP \bu +\GRAD p &=\bef, \\ \DIV \bu &= 0, \\ \bu_{|\Gamma_1} &= \bu_{\text{bdy}} , \\ {\bu \cdot \textbf{n}}_{|\Gamma_2} &= 0 , \\ \bu_{|t=0} &= \bu_0, \\ p_{|t=0} &= p_0. \end{align*}

where \(\textbf{n}\) is the outward normals of \(\partial\Omega_{\bu}\). We also define \(\Gamma_1\) and \(\Gamma_2\) such that \(\Gamma_1 \cup \Gamma_2=\partial \Omega_{\bu}\) and \( \Gamma_1 \cap \Gamma_2 = \emptyset \). The parameter \(\Re\) is the kinetic Reynolds number.

Remark: for magnetohydrodynamic problem, the Lorentz force is added in the right hand-side of the Navier-Stokes equations. It is equal to \((\ROT\bH)\times \bB\).

To approximate such problems with the code SFEMaNS, the following subroutines/functions need to be edited. All required information to edit them are provided in the links below.

The subroutine init_velocity_pressure initializes the velocity field, the pressure and the pressure increment.

The function source_in_NS_momentum defines the source term \(\textbf{f}\).

The function vv_exact defines the Dirichlet boundary condition \(\bu_{\text{bdy}}\). It can also be used by the subroutine init_velocity_pressure to initialize the velocity field.

The function pp_exact can be used by the subroutine init_velocity_pressure to initialize the pressure. We note that Dirichlet boundary conditions on the pressure are not implemented.

The function extension_velocity defines the velocity field in a solid domain where the Navier-Stokes equations are not approximated. It is used when the domains of the temperature or the magnetic field are larger than the domain where the Navier-Stokes equations are approximated.

Problem with variable density and/or dynamical viscosity

The code SFEMaNS can approximate the solutions of the Navier-Stokes equations with variable density and viscosity in an axisymmetric domain \(\Omega_{\bu}\). We assume that the immiscible fluids are stratified. For each interface betweem two fluids, we introduce a level set function, denoted by \(\varphi\). It is equal to 0 in one fluid and 1 in the other. The interface is represented by \(\varphi^{-1}(\{1/2\})\). The level set is the solution of the following equation in \(\Omega_{\bu}\):

\begin{align*} \partial_t \varphi + \bu \cdot \GRAD \varphi = f_\varphi. \end{align*}

For a two phase field, we recontruct the density and dynamical viscosity as follows:

\begin{align*} \rho=\rho_1+(\rho_2-\rho_1) F(\varphi), \\ \eta=\eta_1+(\eta_2-\rho_1) F(\varphi), \end{align*}

with \(\eta_i\) and \(\rho_i\) are the dynamical viscosities and densities of the two fluids. The function \(F(\phi)\) is either the identity (called linear reconstruction) or a piece-wise polynomial function (called reg reconstruction).

We remind that the Navier-Stokes equations are approximated with the momentum \(\bm=\rho \bu\) as dependent variable. We refer to the section Extension to multiphase flow problem for more information on the mathematical formulation of such problems.

In addition to the previous functions described above, the following functions of the file condlim.f90 need to be edited. All required information to edit them are provided in the links below.

The subroutine init_level_set initializes the level set(s).

The function source_in_level_set defines the source term \(f_\varphi \).

The function level_set_exact can be used by the subroutine init_level_set to initialize the level set(s) (the boundary conditions on the level set follow from the ones on the velocity).

Problem with non axisymmetric solid obstacle

The code SFEMaNS can consider solid non axisymmetric obstacle \(\Omega_\text{obs}\) in the domain \(\Omega_{\bu}\) where the Navier-Stokes equations are approximated. The velocity in \(\Omega_\text{obs}\) can be non zero and the domain \(\Omega_\text{obs}\) can be time dependent. Such feature are made possible thanks to a penalty method. We refer to the section Extension to non axisymmetric geometry for more information on the mathematical formulation of such problems.

To consider such problems, the following functions needs to be set properly.

The function penal_in_real_space defines the penalty function \(\chi\) that is equal to 1 in the fluid domain \(\Omega_{\bu} \setminus \Omega_\text{obs} \) and 0 in the domain \(\Omega_\text{obs}\). We note that the penalty function can be continuous (have it go from the values 0 to 1 in a few mesh cells around the interface between the solid and the fluid domain).

The function imposed_velocity_by_penalty defines the velocity in the domain \(\Omega_\text{obs}\).

Problem with ferrofluid

The code SFEMaNS can consider ferrofluid. For such problems, the action of the magnetic field on the velocity field is reported by the Kelvin force. This force is equal to \(g(T) \GRAD \bH^2 \) where \(\bH\) is the magnetic field, \(T\) is the temperature and \(g(T)\) is a coefficient that depends of the temperature. We note that the definition of this coefficient can depend of the problem considered.

To consider such problems, the following function needs to be edited.

The function kelvin_force_coeff defines the coefficient \(g(T)\) depending of the temperature.

Settings for the temperature equation

The code SFEMaNS can approximate the solution of the temperature equation in an axisymmetric domain \(\Omega_{T}\). These equations are written as follows:

\begin{align*} \partial_t T+ \bu \cdot \GRAD T - \kappa \LAP T &= f_T, \\ T_{|\Gamma_1} &= T_\text{bdy} , \\ \partial_{\textbf{n}}T_{|\Gamma_2} &= 0, \\ T_{|t=0} &= T_0, \end{align*}

where \(\textbf{n}\) is the outward normals of \(\partial\Omega_T\). We also define \(\Gamma_1\) and \(\Gamma_2\) such that \(\Gamma_1 \cup \Gamma_2=\partial \Omega_T\) and \( \Gamma_1 \cap \Gamma_2 = \emptyset \). The parameter \(\kappa\) is the thermal diffusivity. We note that for thermohydrodynamic problems, we assume that \(\Omega_{\bu} \subset \Omega_T\).

To approximate such problems with the code SFEMaNS, the following subroutines/functions need to be edited. All required information to edit them are provided in the links below.

The subroutine init_temperature initializes the temperature.

The function source_in_temperature defines the source term \(f_T\).

The function temperature_exact defines the Dirichlet boundary condition \(T_{\text{bdy}}\). It can be used by the subroutine init_temperature to initialize the temperature.

Settings for the Maxwell equations

General case

The code SFEMaNS can approximate the solutions of the Maxwell equations in an axisymmetric domain \(\Omega_\text{mxw}\). The domain is the union of an axisymmetric conducting domain \(\Omega_{\bH}\) and an axisymmetric insulating domain \(\Omega_\phi\). The Maxwell equations are written as follows:

\begin{align} \begin{cases} \partial_t (\mu^c \mathbf{H}) + \nabla \times \left(\frac{1}{\Rm \sigma}\nabla \times \mathbf{H} \right) = \nabla\times (\bu^\text{given} \times \mu^c \mathbf{H}) + \nabla \times \left(\frac{1}{\Rm \sigma} \nabla\times \mathbf{j} \right) & \text{in } \Omega_{\text{H}}, \\ \text{div} (\mu^c \mathbf {H}) = 0 & \text{in } \Omega_{\text{H}} ,\\ - \partial_t \DIV (\mu^v \GRAD( \phi)) = 0 & \text{in } \Omega_\phi ,\\ \bH \times \bn^c + \nabla \phi \times \bn^v = 0 & \text{on } \Sigma, \\ \mu^c \bH \cdot \bn^c + \mu ^v \nabla \phi \cdot \bn^v = 0 & \text{on } \Sigma, \\ \bH \times \bn = \bH_{\text{bdy}} \times \bn & \text{on } \Gamma_{\bH}^1,\\ \left( \frac{1}{\Rm \sigma} \left( \ROT (\mathbf{H}) - \mathbf{j} \right) - \bu \times \mu \mathbf{H} \right) \times \bn = \textbf{a} \times \bn & \text{on } \Gamma_{\bH}^2,\\ \phi = \phi_{\text{bdy}} & \text{on } \Gamma_\phi,\\ \bH_{|t=0}= \bH_0, \\ \phi_{|t=0}= \phi _0, \end{cases} \end{align}

where \(\Sigma= \Omega_{\bH} \cap \Omega_\phi \), \(\Gamma_{\bH}^1 \cup \Gamma_{\bH}^2= \partial \Omega_{\bH}\setminus \Sigma \), \(\Gamma_{\bH}^1 \cap \Gamma_{\bH}^2=\emptyset\) and \(\Gamma_\phi= \partial \Omega_\phi\setminus \Sigma \). We denote by \(\textbf{n}\) the outward normal vector to \(\Gamma_{\bH}^1 \cup \Gamma_{\bH}^2\). The function \(\textbf{j}\) is a source term to define. The parameter \(\Rm\) is the magnetic Reynolds number. We also denoted by \(\sigma\) the conductivity in \(\Omega_{\bH}\) and by \(\mu^c\) the magnetic permeability in in \(\Omega_{\bH}\). They both are piecewise constant. The magnetic permeability in the insulating domain is denoted \(\mu^v\). The velocity \(\) comes from the Navier-Stokes equations or functions of the file condlim.f90. We remind that the scalar potential satisfies the relation \(\bH=\GRAD \phi\) in \(\Omega_\phi\). We note that for magnetohydrodynamic problems, we assume that \(\Omega_{\bu} \subset \Omega_T \subset \Omega_\text{mxw}\).

To approximate such problems with the code SFEMaNS, the following subroutines/functions need to be edited. All required information to edit them are provided in the links below.

The subroutine init_maxwell initializes the magnetic field \(\bH\) and the scalar potential \(\phi\).

The function Hexact defines the Dirichlet boundary condition \(\bH_{\text{bdy}}\). It can be used by the subroutine init_maxwell to initialize the magnetic field.

The function Phiexact defines the Dirichlet boundary condition \(\phi_{\text{bdy}}\). It can be used by the subroutine init_maxwell to initialize the scalar potential.

The function Vexact defines the velocity field \(\bu_\text{given}\) that is a time independent function. It is only used when the solutions of the Navier-stokes equations are not approximated.

The function Jexact_gauss defines the source term \(\textbf{j}\).

The function Eexact_gauss defines the boundary data \(\textbf{a}\).

Quasi-static approximation

The code SFEMaNS allows to approximate the Maxwell equation in a quasi-static regime with the magnetic field \(\bB=\mu\bH\) as dependent variable. The total magnetic field \(\bB\) is written \(\bB_b+\bb\) where \(\bB_b\) represent a basic state that is time independent. The Maxwell equations are approximated with \(\bb\) where the time derivative and the term \(\nabla\times (\bu^\text{given} \times \mathbf{b})\) are disregarded. The term \(\nabla\times (\bu^\text{given} \times \mathbf{B}_b)\) is still computed and that the source term \(\textbf{j}\) has to report the action of the term \(\ROT( \frac{1}{\sigma\Rm}\ROT \bB_b)\). The equations can be written as follows.

\begin{align*} \ROT \left( \frac{1}{\sigma\Rm} \ROT \textbf{b} \right) = \nabla\times (\bu^\text{given} \times \mathbf{B}_b) + \ROT\left( \frac{1}{\sigma\Rm}\textbf{j}\right). \end{align*}

Remark: the Lorentz force in the Navier-Stokes equations is taken equal to \((\ROT\bH_b)\times \bb + (\ROT \bh)\times \bB_b\) with \(\bh =\mu \bb\) and \(\bH_b=\mu \bB_b\).

To approximate a quasi-static regime of the Maxwell equation you need to set properly the following function:

The function H_B_quasi_static defines the base states \(\bB_b\) and \(\bH_b\).

All required information to edit them are provided in the above link.

Problem with a given magnetic permeability variable in \((r,\theta,z,t)\)

The code SFEMaNS can take into account magnetic permeability that are variable in \((r,\theta,z,t)\) with \(t\) being the time. It is done by approximating the Maxwell equations with the magnetic field \(\bB\) as dependent variable. Moreover, the diffusion term \(\ROT (\frac{1}{\sigma\Rm}\ROT(\frac{\bB}{\mu}))\) is made explicit. The numerical scheme is stabilized with a term that involves a function \(\overline{\mu}(r,z)\) that has to be smaller than \(\mu(r,\theta,z,t)\) for all \((r,\theta,z,t)\). We refer to the section Extension to magnetic permeability variable in time and azimuthal direction for more information on the mathematical formulation of such problems.

Remark:

1. If the magnetic permeability does not depend of the time and the azimuthal coordinate \(\theta\), the term \(\ROT( \ROT (\frac{1}{\mu}\bB))\) can be treated implicitly and no stabilization term is added.
2. If the magnetic permeability is a piece wise constant function with re (one value on each axisymmetric domains where the

To approximate such problems with the code SFEMaNS, the following subroutines/functions need to be edited. All required information to edit them are provided in the links below.

The function mu_bar_in_fourier_space. If the magnetic permeability does not depend of \((\theta,t)\), this function defines the magnetic permeability \(\mu\). If \(\mu\) depends of the azimuthal direction \(\theta\) and/or the time \(t\), this function defines a function \(\overline{\mu}(r,z)\). It satisfies the relation \(\overline{\mu}(r,z)\leq \mu(r,\theta,z,t)\) for all \((r,\theta,z,t)\) considered.

The function grad_mu_bar_in_fourier_space defines the gradient in the radial and vertical directions of the function mu_bar_in_fourier_space described above.

The function mu_in_real_space defines the magnetic permeability in \(\Omega_\text{mxw}\) when its depends of the time and/or the azimuthal direction.

Multiphase problem with variable electrical conductivity

The code SFEMaNS can approximate the Navier-Stokes equations for stratified immiscible fluids. The electrical conductivity of these fluids may be different. As a consequence, the electrical conductivity is variable in \((r,\theta,z,t)\) with t being the time. To approximate the Maxwell equations for such problems, the diffusion term \(\ROT (\frac{1}{\sigma\Rm}\ROT(\frac{\bB}{\mu}))\) is made explicit. The numerical scheme is stabilized with a term that involves a function \(\overline{\sigma}(r,z)\) that has to be smaller than \(\sigma(r,\theta,z,t)\) for all \((r,\theta,z,t)\). We refer to the section Extension to multiphase flow problem for more information on the mathematical formulation of such problems.

Remark:

1. The magnetic permeability \(\mu\) needs to be constant.
2. The Maxwell equations can either be approximated with \(\bB\) or \(\bH\) as dependent variable.

To approximate such problems with the code SFEMaNS, the following function needs to be edited.

The function sigma_bar_in_fourier_space. It defines the stabilization term \(\overline{\sigma}\) such that \(\overline{\sigma}(r,z)\leq \sigma\) for all \((r,\theta,z,t)\) considered.