Navier-Stokes model is used to model incompressible Newtonian fluid. It can be described by these conservative laws :

we complete this set of equations with the fluid constitutive law

with strain tensor \(D(\boldsymbol{u}_{f})\) defined by :

An alternative model is the Stokes model. It is valid in the case of small Reynolds number. It corresponds to the same formulation than Navier-Stokes equations but without the convective term \(\left( \boldsymbol{u}_{f} \cdot \nabla_{\mathrm{x}} \right) \boldsymbol{u}_{f}\) .

Let’s define a bounded domain \(\Omega \subset \mathbb{R}^p\) (\(p=2,3\)) decomposed into two subdomains \(\Omega_1\) and \(\Omega_2\). We denote \(\Gamma\) the interface between the two partitions. The goal of the level set method is to track implicitly the interface \(\Gamma(t)\) moving at a velocity \(\mathbf{u}\). The level set method has been described in [osher] and its main ingredient is a continuous scalar function \(\phi\) (the /level set/ function) defined on the whole domain. This function is chosen to be positive in \(\Omega_1\), negative in \(\Omega_2\) and zero on \(\Gamma\). The motion of the interface is then described by the advection of the level set function with a divergence free velocity field \(\mathbf{u}\):

A convenient choice for \(\phi\) is a signed distance function to the interface. Indeed, the property \(|\nabla \phi| = 1\) of distance functions eases the numerical solution and gives a convenient support for delta and Heaviside functions).

Nevertheless, it is known that the advection equation \eqref{eq:advection} does not conserve the property \(|\nabla \phi|=1\). Thus, when \(|\nabla \phi|\) is far from \(1\) we have to reset \(\phi\) as a distance function without moving the interface. To do so we can either use an Hamilton-Jacobi method or the fast marching method (see \cite{Winkelmann2007} for details about the fast marching method).

In two-fluid flow simulations, we need to define some quantities related to the interface such as the density, the viscosity, or some interface forces. To this end, we introduce the smoothed Heaviside and delta functions:

where \(\varepsilon\) is a parameter defining a ``numerical thickness'' of the interface. A typical value of \(\varepsilon\) is \(1.5 h\) where \(h\) is the mesh size of elements crossed by the iso-value \(0\) of the level set function.

The Heaviside function is used to define parameters having different values on each subdomains. For example, we define the density of two-fluid flow as \(\rho = \rho_2 + (\rho_1-\rho_2) H_{\varepsilon}(\phi)\) (we use a similar expression for the viscosity \(\nu\)). Regarding the delta function, it is used to define quantities on the interface. In particular, in the variational formulations, we replace integrals over the interface \(\Gamma\) by integrals over the entire domain \(\Omega\) using the smoothed delta function. If \(\phi\) is a signed distance function, we have : \(\int_{\Gamma} 1 \simeq \int_{\Omega} \delta_{\varepsilon}(\phi)\). If \(\phi\) is not close enough to a distance function, then \(\int_{\Gamma} 1 \simeq \int_{\Omega} |\nabla \phi| \delta_{\varepsilon}(\phi)\) which still tends to the measure of \(\Gamma\) as \(\varepsilon\) vanishes. However, if \(\phi\) is not a distance function, the support of \(\delta_{\varepsilon}\) can have a different size on each side of the interface. More precisely, the support of \(\delta_{\varepsilon}\) is narrowed on the side where \(|\nabla \phi|>1\) and enlarged on regions where \(|\nabla \phi|<1\). It has been shown in [cottet] that replacing \(\phi\) by \(\frac{\phi}{|\nabla \phi|}\) has the property that \(|\nabla \frac{\phi}{|\nabla \phi|}| \simeq 1\) near the interface and has the same iso-value \(0\) as \(\phi\). Thus, replacing \(\phi\) by \(\frac{\phi}{|\nabla \phi|}\) as support of the delta function does not move the interface. Moreover, the spread interface has the same size on each part of the level-set \(\phi=0\). It reads \(\int_{\Gamma} 1 \simeq \int_{\Omega} \delta_{\varepsilon}(\frac{\phi}{| \nabla \phi|})\). The same technique is used for the Heaviside function.

We use Feel++ to discretize and solve the problem. Equation \eqref{eq:advection} is solved using a stabilized finite element method. We have implemented several stabilization methods such as Streamline Upwind Diffusion (SUPG), Galerkin Least Square (GLS) and Subgrid Scale (SGS). A general review of these methods is available in [Franca92]. Other available methods include the Continuous Interior Penalty method (CIP) are described in \cite{Burman2006}. The variational formulation at the semi-discrete level for the stabilized equation \eqref{eq:advection} reads, find \(\phi_h \in {\mathbb R}_h^k\) such that \(\forall \psi_h \in {\mathbb R}_h^k\) :

where \(S(\cdot, \cdot)\) stands for the stabilization bilinear form (see section \ref{sec:membr-inext} for description of \({\mathbb R}_h^k\) and \(\mathbf{u}_h\)). In our case, we use a BDF2 scheme which needs the solution at the two previous time step to compute the one at present time. For the first time step computation or after a reinitialization we use an Euler scheme.

The two-fluid flow problem can be expressed by

Where \(\boldsymbol{f}_\phi\) is the force obtain by projection of the density of interfacial forces on the domain \(\Omega\)

The CFD Toolbox supports both the Stokes and the incompressible Navier-Stokes equations.

The fluid mechanics model (Navier-Stokes or Stokes) can be selected in json file:

The next step is to define the fluid material by setting its properties namely the density \(\rho_f\) and viscosity \(\mu_f\). In next table, we find the correspondance between the mathematical names and the json names.

A Materials section is introduced in json file in order to configure the fluid properties. For each mesh marker, we can define the material properties associated.

We start by a listing of boundary conditions supported in fluid mechanics model.

Body forces are also defined in BoundaryConditions category in json file.

The marker corresponds to mesh elements marked with this tag. If the marker is an empty string, it corresponds to all elements of the mesh.

programme avalaible :