This page details the numerical methods implemented in TrioCFD for the incompressible model defined by the governing equations (incompressible, variable viscosity). Two spatial discretization methods are available: the Finite-Element Volume (VEF) method and the Finite-Difference Volume (VDF) method; only the VEF part is described here. The computational domain is denoted \(\Omega\).

Time scheme (valid for VEF and VDF)

After discretization, the matrix-vector system solved reads:

\[\left\{\begin{array}{rcl} \delta t^{-1}\mathbf{M}U^{n+1}+\mathbf{A}U^{n+1}+\mathbf{L}(U^{n})U^{n+1}+\mathbf{B}^{T}P^{n+1} & = & F^{n}+\delta t^{-1}\mathbf{M}U^{n},\\ \mathbf{B}U^{n+1} & = & 0, \end{array}\right. \]

where \(U^{n+1}\in\mathbb{R}^{N_{\mathbf{u}}}\) is the discretized velocity at time \((n+1)\delta t\) ( \(N_{\mathbf{u}}\) velocity degrees of freedom), \(P^{n+1}\in\mathbb{R}^{N_p}\) the discretized pressure ( \(N_p\) pressure degrees of freedom), and the bold matrices depend on the spatial discretization. To decouple velocity and pressure, the resolution is performed in three steps [21], [170] :

Prediction step: compute \(U^{*}\) solving \(\delta t^{-1}\mathbf{M}U^{*}+\mathbf{A}U^{*}+\mathbf{L}(U^{n})U^{*}+\mathbf{B}^{T}P^{n}=F^{n}+\delta t^{-1}\mathbf{M}U^{n}\) (at this step \(\mathbf{B}U^{*}\neq0\));
Pressure computation: compute \(P'\) solving \(\mathbf{B}\mathbf{M}^{-1}\mathbf{B}^{T}P'=\delta t^{-1}\mathbf{B}U^{*}\), then \(P^{n+1}=P'+P^{n}\);
Correction step: compute \(U^{n+1}\) solving \(\mathbf{M}U^{n+1}=\mathbf{M}U^{*}-\delta t\,\mathbf{B}^{T}P'\).

Note: Several solvers were tested (SIMPLE, SIMPLER, PISO) which, depending on the problems, showed relatively weak convergence. The solver used to date is inspired by [70].

VEF space scheme

The numerical method is based on the non-conforming Crouzeix-Raviart finite-element method [25], detailed in [48], [72], [57], [6]. For \(d=2\) (resp. \(d=3\)), a regular mesh \(\mathcal{T}_{h}=\cup_{\ell=1}^{N_T}T_{\ell}\) of the domain \(\Omega\) is composed of simplices (triangles in 2D, tetrahedra in 3D), with vertices \((S_{i})_{i=1}^{N_S}\); the boundary of \(T_{\ell}\) is made of edges (2D) or faces (3D), called "faces" in all cases. Let \(\overline{\mathcal{F}}_{h}=\cup_{k=1}^{\overline{N}_F}F_{k}\) be the set of all faces and \(\mathcal{F}_{h}=\cup_{k=1}^{N_F}F_{k}\) the set of interior faces; \(M_k\) denotes the barycentre of face \(F_k\) and \(\mathbf{n}_k\) its outward unit normal.

Let \(P_{1}(T)\) be the set of order-1 polynomials on \(T\). The velocity discretization space is

\[X_{h}:=\{\mathbf{v}_{h}\,|\,\forall T\in\mathcal{T}_{h},\,\mathbf{v}_{h}\in P_{1}(T)^{d}\,\text{ and }\,\forall F\in\mathcal{F}_{h}:\,[\mathbf{v}_{h}](\mathbf{x}_{F})=\mathbf{0}\} \]

where \(\mathbf{x}_{F}\) is the barycentre of face \(F\) and \([\mathbf{v}_{h}](\mathbf{x}_{F})\) is the jump of \(\mathbf{v}_{h}\) across \(F\). The space \(X_{h}\) is equipped with the semi-norm \(\|\mathbf{v}_{h}\|_{h}=(\sum_{\ell}|\mathbf{v}_{h|T_{\ell}}|_{1,T_{\ell}}^{2})^{1/2}\), and \(X_{0,h}:=\{\mathbf{v}_{h}\in X_{h}\,|\,\mathbf{u}_{h|\partial\Omega}=0\}\). With \(\lambda_i|_T\) the barycentric coordinate associated with vertex \(S_i|_T\) and \((\mathbf{e}^{\beta})_{\beta=1}^{d}\) the canonical basis vectors, the basis functions associated with \(X_{h}\) are the vectors \(\boldsymbol{\varphi}_{i}^{\beta}|_{T}=(1-d\lambda_{i}|_{T})\mathbf{e}^{\beta}\); \(\psi_j\) denotes the characteristic function associated with triangle \(T_j\).

Volume-element matrices for the P1NC/P0 approximation

Let \((\mathbf{u}_{h}^{n},p_{h}^{n})\in X_{h}\times L_{h}\) be the spatial approximation of \((\mathbf{u}^{n},p^{n})\), with \(\mathbf{u}_{h}^{n}=\sum_{\beta=1}^{d}\sum_{i=1}^{N_F}(U_{i}^{\beta})^{n}\boldsymbol{\varphi}_{i}^{\beta}\) and \(p_{h}^{n}=\sum_{\ell=1}^{N_T}P_{\ell}^{n}\psi_{\ell}\), and \(N_{\mathbf{u}}:=d\,N_F\).

Mass matrix \(\mathbf{M}\in\mathbb{R}^{N_{\mathbf{u}}\times N_{\mathbf{u}}}\): \(d\) diagonal blocks equal to \(\mathbf{M}_F\) with \((\mathbf{M}_F)_{i,j}=(\phi_i,\phi_j)_0\), where \(\phi_i|_T=(1-d\lambda_i|_T)\). In 2D the mass matrix is diagonal, \((\mathbf{M}_F)_{i,j}=\delta_{ij}\sum_{\ell\,|\,M_i\in T_{\ell}}\frac{|T_{\ell}|}{3}\), and equals the finite-volume mass matrix of Emonot's thesis [48]. This is no longer the case in 3D, but TrioCFD implements the finite-volume mass matrix, for which \((\mathbf{M}_F)_{i,j}=\delta_{ij}\sum_{\ell\,|\,M_i\in T_{\ell}}\frac{|T_{\ell}|}{4}\).
Stiffness matrix \(\mathbf{A}\in\mathbb{R}^{N_{\mathbf{u}}\times N_{\mathbf{u}}}\): \(d\) diagonal blocks equal to \(\mathbf{A}_F\) with \((\mathbf{A}_F)_{i,j}=(\nabla\phi_i,\nabla\phi_j)_0=\sum_{\ell\,|\,M_i,M_j\in T_{\ell}}|T_{\ell}|^{-1}\mathbf{S}_{i,\ell}\cdot\mathbf{S}_{j,\ell}\), where \(\mathbf{S}_{i,\ell}\) is the "face normal" vector of the face opposite vertex \(S_i\) in triangle \(T_{\ell}\). It also equals the finite-volume stiffness matrix of Emonot's thesis.
Coupling matrix \(\mathbf{B}\in\mathbb{R}^{N_T\times N_{\mathbf{u}}}\): composed of \(d\) blocks \(\mathbf{B}=(\mathbf{B}^{\beta})_{\beta=1}^{d}\), \(\mathbf{B}^{\beta}\in\mathbb{R}^{N_T\times N_F}\), with \((\mathbf{B}^{\beta})_{\ell,j}=-\mathbf{S}_{j,\ell}\cdot\mathbf{e}^{\beta}\).
** \(P_1\) pressure case:** when the pressure is \(P_1\), \(p_{h}^{n}=(P_{i}^{n})_{i}\in\mathbb{R}^{N_S}\) and \(\mathbf{B}\in\mathbb{R}^{N_S\times N_{\mathbf{u}}}\) with \((\mathbf{B}^{\beta})_{i,k}=-\sum_{\ell\,|\,M_k,S_i\in T_{\ell}}((d+1)d)^{-1}\mathbf{S}_{i,\ell}\cdot\mathbf{e}^{\beta}\).