A second-order face-centred finite volume method for elliptic problems

TL;DR

This paper introduces a second-order face-centred finite volume method (FCFV) that defines solutions on mesh faces, offering robustness against mesh distortion and eliminating the need for gradient reconstruction, applicable to complex 3D geometries.

Contribution

It presents a novel second-order FCFV method based on a mixed formulation that computes solutions on faces, improving robustness and efficiency over traditional cell-centred and vertex-centred FV methods.

Findings

Method achieves second-order accuracy.

Insensitive to mesh distortion and stretching.

Effective in complex 3D geometries.

Abstract

A second-order face-centred finite volume method (FCFV) is proposed. Contrary to the more popular cell-centred and vertex-centred finite volume (FV) techniques, the proposed method defines the solution on the faces of the mesh (edges in two dimensions). The method is based on a mixed formulation and therefore considers the solution and its gradient as independent unknowns. They are computed solving an element-by-element problem after the solution at the faces is determined. The proposed approach avoids the need of reconstructing the solution gradient, as required by cell-centred and vertex-centred FV methods. This strategy leads to a method that is insensitive to mesh distortion and stretching. The current method is second-order and requires the solution of a global system of equations of identical size and identical number of non-zero elements when compared to the recently proposed first-order FCFV. The formulation is presented for Poisson and Stokes problems. Numerical examples are used to illustrate the approximation properties of the method as well as to demonstrate its potential in three dimensional problems with complex geometries. The integration of a mesh adaptive procedure in the FCFV solution algorithm is also presented.

Equations

$$\{\begin{array}{rlrl}-\text{boldmath∇unboldmath}\cdot \text{boldmath∇unboldmath}u& =s& & \text{in Ω,}\\ u& =u_D& & \text{on ΓD,}\\ \text{boldmathnunboldmath}\cdot \text{boldmath∇unboldmath}u& =t& & \text{on ΓN,}\end{array}$$

$$\Omega =\bigcup _{e=1}^{{\text{n}}_{\text{el}}}{\Omega }_e,{\Omega }_e\cap {\Omega }_l=\varnothing \text{ for }e\ne l,$$

$$\Gamma:=\Big{[}\bigcup_{e=1}^{\texttt{n}_{\texttt{el}}}\partial\Omega_{e}\Big{]}\setminus\partial\Omega$$

$$\partial {\Omega }_e:=\bigcup _{j=1}^{{\text{n}}_{\text{fa}}^e}{\Gamma }_{e,j},$$

$$\{\begin{array}{rlrl}\text{boldmathqunboldmath}+\text{boldmath∇unboldmath}u& =\text{boldmath0unboldmath}& & \text{in Ωe, and for e=1,…,nel,}\\ \text{boldmath∇unboldmath}\cdot \text{boldmathqunboldmath}& =s& & \text{in Ωe, and for e=1,…,nel,}\\ u& =u_D& & \text{on ΓD,}\\ \text{boldmathnunboldmath}\cdot \text{boldmathqunboldmath}& =-t& & \text{on ΓN,}\\ [[u\text{boldmathnunboldmath}]]& =\text{boldmath0unboldmath}& & \text{on Γ,}\\ [[\text{boldmathnunboldmath}\cdot \text{boldmathqunboldmath}]]& =0& & \text{on Γ,}\end{array}$$

$$\{\begin{array}{rlrl}{\text{boldmathqunboldmath}}_e+\text{boldmath∇unboldmath}{u}_e& =\text{boldmath0unboldmath}& & \text{in Ωe, }\\ \text{boldmath∇unboldmath}\cdot {\text{boldmathqunboldmath}}_e& =s& & \text{in Ωe,}\\ u_e& =u_D& & \text{on ∂Ωe∩ΓD,}\\ u_e& =\hat{u}& & \text{on ∂Ωe∖ΓD,}\end{array}$$

$$\{\begin{array}{rlrl}[[u\text{boldmathnunboldmath}]]& =\text{boldmath0unboldmath}& & \text{on Γ,}\\ [[\text{boldmathnunboldmath}\cdot \text{boldmathqunboldmath}]]& =0& & \text{on Γ,}\\ \text{boldmathnunboldmath}\cdot \text{boldmathqunboldmath}& =-t& & \text{on ΓN.}\end{array}$$

$$-{\int }_{{\Omega }_e}{\text{boldmathqunboldmath}}_e^{h}d\Omega ={\int }_{\partial {\Omega }_e\cap {\Gamma }_D}{u}_D{\text{boldmathnunboldmath}}_{e}d\Gamma +{\int }_{\partial {\Omega }_e\setminus {\Gamma }_D}{\hat{u}}^h{\text{boldmathnunboldmath}}_{e}d\Gamma ,$$

$$-{\int }_{{\Omega }_e}\text{boldmath∇unboldmath}v\cdot {\text{boldmathqunboldmath}}_e^{h}d\Omega +{\int }_{\partial {\Omega }_e}v({\text{boldmathnunboldmath}}_e\cdot {\widehat{\text{boldmathqunboldmath}}}_e^h)d\Gamma ={\int }_{{\Omega }_e}vsd\Omega$$

$${\text{boldmathnunboldmath}}_e\cdot {\widehat{\text{boldmathqunboldmath}}}_e^h:=\{\begin{array}{ll}{\text{boldmathnunboldmath}}_e\cdot {\text{boldmathqunboldmath}}_e^h+{\tau }_e({\mathbb{P}}_0{u}_e^h-u_D)& \text{on ∂Ωe∩ΓD,}\\ {\text{boldmathnunboldmath}}_e\cdot {\text{boldmathqunboldmath}}_e^h+{\tau }_e({\mathbb{P}}_0{u}_e^h-{\hat{u}}^h)& \text{elsewhere,}\end{array}$$

$$-{\int }_{{\Omega }_e}{\text{boldmathqunboldmath}}_e^{h}d\Omega ={\int }_{\partial {\Omega }_e\cap {\Gamma }_D}{u}_D{\text{boldmathnunboldmath}}_{e}d\Gamma +{\int }_{\partial {\Omega }_e\setminus {\Gamma }_D}{\hat{u}}^h{\text{boldmathnunboldmath}}_{e}d\Gamma ,$$

$${\int }_{\partial {\Omega }_e}v{\tau }_e{\mathbb{P}}_0{u}_e^{h}d\Gamma ={\int }_{{\Omega }_e}vsd\Omega +{\int }_{\partial {\Omega }_e\cap {\Gamma }_D}v{\tau }_e{u}_{D}d\Gamma +{\int }_{\partial {\Omega }_e\setminus {\Gamma }_D}v{\tau }_e{\hat{u}}^{h}d\Gamma$$

$$\sum _{e=1}^{{\text{n}}_{\text{el}}}{\int }_{\partial {\Omega }_e\setminus {\Gamma }_D}{\text{boldmathnunboldmath}}_e\cdot {\widehat{\text{boldmathqunboldmath}}}_e^{h}d\Gamma =-\sum _{e=1}^{{\text{n}}_{\text{el}}}{\int }_{\partial {\Omega }_e\cap {\Gamma }_N}td\Gamma ,$$

$$\sum _{e=1}^{{\text{n}}_{\text{el}}}{\int }_{\partial {\Omega }_e\setminus {\Gamma }_D}\left({\text{boldmathnunboldmath}}_e\cdot {\text{boldmathqunboldmath}}_e^h+{\tau }_e({\mathbb{P}}_0{u}_e^h-{\hat{u}}^h)\right)d\Gamma =-\sum _{e=1}^{{\text{n}}_{\text{el}}}{\int }_{\partial {\Omega }_e\cap {\Gamma }_N}td\Gamma .$$

$${\mathcal{A}}_e$$

$${\mathcal{D}}_e$$

$${\mathcal{N}}_e$$

$${\mathcal{B}}_e$$

$$\chi_{\square}(l)=\bigg{\{}\begin{array}[]{ll}1&\text{ if }\ l\in\square\\ 0&\text{ otherwise}.\end{array}$$

$$-|{\Omega }_e|{\mathbf{q}}_e=\sum _{j\in {\mathcal{D}}_e}|{\Gamma }_{e,j}|{\text{boldmathnunboldmath}}_j{u}_{D,j}+\sum _{j\in {\mathcal{B}}_e}|{\Gamma }_{e,j}|{\text{boldmathnunboldmath}}_j{\hat{\text{u}}}_j,$$

$${\mathbf{m}}_e{\mathbf{u}}_e={\mathbf{f}}_e+\sum _{j\in {\mathcal{D}}_e}{\tau }_j{\mathbf{d}}_j+\sum _{j\in {\mathcal{B}}_e}{\tau }_j{\mathbf{r}}_j{\hat{\text{u}}}_j,$$

$$(m_e{)}_{IJ}$$

$$(d_j{)}_I$$

$$(p_{e,i}{)}_l=\frac{1}{{\text{n}}_{\text{fn}}^{e,i}}{\chi }_{{\mathcal{F}}_{e,i}}(l)$$

$${\mathbf{q}}_e=-|{\Omega }_e{|}^{-1}{\mathbf{z}}_e-|{\Omega }_e{|}^{-1}\sum _{j\in {\mathcal{B}}_e}|{\Gamma }_{e,j}|{\text{boldmathnunboldmath}}_j{\hat{\text{u}}}_j,$$

$${\mathbf{u}}_e={\mathbf{m}}_e^{-1}{\mathbf{b}}_e+{\mathbf{m}}_e^{-1}\sum _{j\in {\mathcal{B}}_e}{\tau }_j{\mathbf{r}}_j{\hat{\text{u}}}_j,$$

$${\mathbf{b}}_e:={\mathbf{f}}_e+\sum _{j\in {\mathcal{D}}_e}{\tau }_j{\mathbf{d}}_j,{\mathbf{z}}_e:=\sum _{j\in {\mathcal{D}}_e}|{\Gamma }_{e,j}|{\text{boldmathnunboldmath}}_j{u}_{D,j}.$$

$$\sum_{e=1}^{\texttt{n}_{\texttt{el}}}\Bigl{\{}|\Gamma_{e,i}|\text{\boldmath$n$\unboldmath}_{i}\cdot\mathbf{q}_{e}+|\Gamma_{e,i}|\tau_{i}\mathbf{p}_{e,i}\cdot\mathbf{u}_{e}-|\Gamma_{e,i}|\tau_{i}\hat{\text{u}}_{i}\Bigr{\}}=-\sum_{e=1}^{\texttt{n}_{\texttt{el}}}\Bigl{\{}|\Gamma_{e,i}|t_{i}\,\chi_{\mathcal{N}_{e}}(i)\Bigr{\}}\quad.$$

$$\hat{\mathbf{K}}\hat{\mathbf{u}}=\hat{\mathbf{f}}.$$

$$\displaystyle{\widehat{K}}^{e}_{i,j}:=|\Gamma_{e,i}|\Big{(}\tau_{i}\tau_{j}\mathbf{p}_{e,i}\cdot\left(\mathbf{m}_{e}^{-1}\mathbf{r}_{j}\right)-|\Omega_{e}|^{-1}|\Gamma_{e,j}|\text{\boldmath$n$\unboldmath}_{i}\cdot\text{\boldmath$n$\unboldmath}_{j}-\tau_{i}\delta_{ij}\Big{)},$$

$$\displaystyle{\widehat{f}}^{e}_{i}:=|\Gamma_{e,i}|\Big{(}|\Omega_{e}|^{-1}\text{\boldmath$n$\unboldmath}_{i}\cdot\mathbf{z}_{e}-\tau_{i}\mathbf{p}_{e,i}\cdot\left(\mathbf{m}_{e}^{-1}\mathbf{b}_{e}\right)-t_{i}\,\chi_{\mathcal{N}_{e}}(i)\Big{)},$$

$${\int }_{{\Gamma }_{e,i}}{\mathbb{P}}_0{u}_e^{h}d\Gamma ={\int }_{{\Gamma }_{e,i}}{u}_e^{h}d\Gamma .$$

$${\int }_{{\Gamma }_{e,i}}{\mathbb{P}}_0{u}_e^{h}d\Gamma =|{\Gamma }_{e,i}|{\mathbb{P}}_0{u}_e^h=|{\Gamma }_{e,i}|{\mathbf{p}}_{e,i}\cdot {\mathbf{u}}_e$$

$${\int }_{{\Gamma }_{e,i}}{u}_e^{h}d\Gamma =\sum _{l=1}^{{\text{n}}_{\text{en}}}\left({\int }_{{\Gamma }_{e,i}}{N}_{l}d\Gamma \right)u_l=\sum _{l=1}^{{\text{n}}_{\text{en}}}\left(\frac{1}{{\text{n}}_{\text{fn}}^{e,i}}|{\Gamma }_{e,i}|{\chi }_{{\mathcal{F}}_{e,i}}(l)\right)u_l.$$

$$\left\{\begin{aligned} \text{\boldmath$L$\unboldmath}+\sqrt{\nu}\text{\boldmath$\nabla$\unboldmath}\text{\boldmath$u$\unboldmath}&=\text{\boldmath$0$\unboldmath}&&\text{in $\Omega_{e}$, and for $e=1,\dotsc,\texttt{n}_{\texttt{el}}$,}\\ \text{\boldmath$\nabla$\unboldmath}\cdot\bigl{(}\sqrt{\nu}\text{\boldmath$L$\unboldmath}+p\mathbf{I}_{\texttt{n}_{\texttt{sd}}}\bigr{)}&=\text{\boldmath$s$\unboldmath}&&\text{in $\Omega_{e}$, and for $e=1,\dotsc,\texttt{n}_{\texttt{el}}$,}\\ \text{\boldmath$\nabla$\unboldmath}\cdot\text{\boldmath$u$\unboldmath}&=0&&\text{in $\Omega_{e}$, and for $e=1,\dotsc,\texttt{n}_{\texttt{el}}$,}\\ \text{\boldmath$u$\unboldmath}&=\text{\boldmath$u$\unboldmath}_{D}&&\text{on $\Gamma_{D}$,}\\ \text{\boldmath$n$\unboldmath}\cdot\bigl{(}\sqrt{\nu}\text{\boldmath$L$\unboldmath}+p\mathbf{I}_{\texttt{n}_{\texttt{sd}}}\bigr{)}&=-\text{\boldmath$t$\unboldmath}&&\text{on $\Gamma_{N}$,}\\ \llbracket\text{\boldmath$u$\unboldmath}\otimes\text{\boldmath$n$\unboldmath}\rrbracket&=\text{\boldmath$0$\unboldmath}&&\text{on $\Gamma$,}\\ \llbracket\text{\boldmath$n$\unboldmath}\cdot\bigl{(}\sqrt{\nu}\text{\boldmath$L$\unboldmath}+p\mathbf{I}_{\texttt{n}_{\texttt{sd}}}\bigr{)}\rrbracket&=\text{\boldmath$0$\unboldmath}&&\text{on $\Gamma$.}\\ \end{aligned}\right.$$

Full text

Keywords: finite volume method, face-centred, second-order convergence, hybridisable discontinuous Galerkin

1 Introduction

Cell-centred and vertex-centred second-order finite volume (FV) methods are still the predominant techniques used in commercial and industrial computational fluid dynamics (CFD) solvers due to their robustness, easy implementation and relatively low cost [17, 33, 19, 2, 15, 14]. Both, cell-centred and vertex-centred, FV techniques require a reconstruction of the gradient of the solution to ensure second-order convergence of the unknown and first-order convergence of the fluxes [12, 11, 3, 1]. The accuracy of the scheme is therefore dependent on the accuracy of the reconstruction technique, which in turns depends on the quality of the mesh. In particular, FV methods are known to suffer an important loss of accuracy, and sometimes even a loss of second-order convergence, when unstructured meshes are used or highly stretched and deformed cells are present in the computational mesh [12, 11].

In [28], the authors proposed a novel methodology, called face-centred finite volume (FCFV) method. The technique can be seen as a particular case of the hybridisable discontinuous Galerkin (HDG) method by Cockburn and co-workers [5, 6, 7, 8], by considering an element-by-element constant degree of approximation. Contrary to other FV methods, the resulting FCFV method is insensitive to mesh distortion and stretching. In addition, being a mixed method, it provides first-order convergence of the gradient of the solution without the need of a reconstruction. Contrary to other mixed methods, in the context of incompressible flows, the FCFV method passes the so-called Ladyzhenskaya-Babuška-Brezzi (LBB) condition. However, the FCFV method only provides first-order convergence for the solution.

This paper proposes a second-order FCFV with a computational cost almost identical to the original, first-order, FCFV. The method is also insensitive to mesh distortion and stretching and provides a first-order approximation of the gradient of the solution without the need of a reconstruction. Similar to the first-order FCFV method, the proposed method can be seen as a particular case of the HDG method, in which the space of approximation used for the solution is linear whereas constant degree approximation spaces are used for the gradient of the solution and the solution on the cell faces [24, 25, 26]. Therefore, the FCFV method inherits the convergence properties of HDG and it passes the LBB condition in the context of incompressible flows. It is worth noting that the degrees of freedom of the solution and its gradient inside each cell can be eliminated via a Schur complement procedure. The only global coupled degrees of freedom correspond to the value of the solution on the cell faces. For both the original and the second-order FCFV methods, the solution on the cell faces is approximated with a constant degree approximation. Hence both methods have a global matrix with the same size and same number of non-zero elements. Due to the extra accuracy of the proposed method, compared to the original FCFV [28, 29], this work also proposes a combination of first and second-order methods to produce an error indicator that is used to drive an $h$-adaptivity process.

The remaining of the paper is organised as follows. In Section 2 the proposed second-order FCFV method is presented for the solution of the Poisson equation. Its extension to the Stokes problem is presented in Section 3. The ability to combine first-order and second-order FCFV methods to perform an automatic mesh adaptive process is discussed in Section 4. Section 5 presents a number of numerical experiments to validate the optimal approximation properties of the method and to compare the accuracy of the first-order and second-order FCFV methods in terms of the computational cost. The insensitivity to mesh distortion and stretching is also demonstrated using numerical experiments. Section 6 presents more challenging problems to show the potential of the proposed second-order FCFV method and its application in an automatic mesh adaptive process. Finally, Section 7 summarises the conclusions of the work that has been presented.

2 Second-order FCFV for the Poisson equation

2.1 Problem statement

Let us consider an open bounded domain $\Omega\in\mathbb{R}^{\texttt{n}_{\texttt{sd}}}$ with boundary $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$, $\Gamma_{D}\cap\Gamma_{N}=\emptyset$ and $\texttt{n}_{\texttt{sd}}$ the number of spatial dimensions. The strong form of the Poisson problem is

[TABLE]

where $s$ denotes a source term, $n$ is the outward unit normal vector to the boundary $\partial\Omega$ and $u_{D}$ and $t$ are the Dirichlet and Neumann data.

The domain $\Omega$ is assumed to be partitioned in $\texttt{n}_{\texttt{el}}$ disjoint triangular or tetrahedral cells $\Omega_{e}$ in two and three dimensions respectively

[TABLE]

with boundaries $\partial\Omega_{e}$, defining an internal interface

[TABLE]

The boundary of a cell is also expressed as the union of a set of edges or faces in two and three dimensions respectively, namely

[TABLE]

where $\texttt{n}_{\texttt{fa}}^{e}$ denotes the number of edges/faces of the cell $\Omega_{e}$. For triangular cells $\texttt{n}_{\texttt{fa}}^{e}=3$ and for tetrahedral cells $\texttt{n}_{\texttt{fa}}^{e}=4$.

The proposed FCFV method uses the mixed form of the Poisson problem in the so-called broken computational domain, namely

[TABLE]

where, following [18], the jump operator is defined as the sum of the values from the left and right of an interface, that is $\llbracket\odot\rrbracket=\odot_{e}+\odot_{l}$.

It is worth noting that the last two equations in (5) impose the continuity of the solution and the normal flux across the interface $\Gamma$, respectively.

2.2 Strong form of the local and global problems

As it is usual in HDG methods [8, 21, 20, 31] and FCFV methods, the strong form of the problem is split into the so-called local problem, defined element-by-element,

[TABLE]

for $e=1,\dotsc,\texttt{n}_{\texttt{el}}$, and the global problem, defined over the interface $\Gamma$ and the Neumann boundary

[TABLE]

The local problem of Equation (6) is a pure Dirichlet problem defined on each cell and introduces the value of the solution at the cell faces, $\hat{u}$, as an independent variable. This problem is used to provide an explicit expression of the solution and its gradient, $u$ and $q$, in terms of the solution at the faces, $\hat{u}$. When this explicit expression is introduced in Equation (7) a global problem with the only unknown being $\hat{u}$ is obtained. This means that the dominant cost of the method is associated to solving the global problem. In a second step the solution and its gradient are computed in each cell using the explicit expression derived from the local problem but this step can be easily parallelised and its cost is negligible compared to the cost of the global problem.

Remark 1**.**

The first equation in (7) is automatically satisfied due to the imposition of the Dirichlet boundary condition in the local problem and the unique value of $\hat{u}$ on the interior faces

2.3 Second-order FCFV weak formulation

Let us denote by $\mathcal{V}^{1}(\Omega)$ the space of $\mathcal{L}_{2}(\Omega)$ functions that are, at most, linear in each cell, $\mathcal{V}^{0}(\Omega)$ the space of $\mathcal{L}_{2}(\Omega)$ functions that are constant in each cell and $\mathcal{\hat{V}}^{0}(\Gamma)$ the space of $\mathcal{L}_{2}(\Gamma)$ functions that are constant on each cell face. With these definitions, the discrete weak formulation of the local problem is: find (u_{e}^{h},\text{\boldmathq\unboldmath}_{e}^{h})\in\mathcal{V}^{1}(\Omega_{e})\times[\mathcal{V}^{0}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}} such that

[TABLE]

for all $v\in\mathcal{V}^{1}(\Omega_{e})$ and for $e=1,\dotsc,\texttt{n}_{\texttt{el}}$. It is worth noting that in Equation (8), a constant test function has been arbitrarily chosen in the space $[\mathcal{V}^{0}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}}$ and it has been used that \text{\boldmathq\unboldmath}_{e}^{h}\in[\mathcal{V}^{0}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}}, that is \text{\boldmath\nabla\unboldmath}\cdot\text{\boldmathq\unboldmath}_{e}^{h}=0.

The so-called numerical flux, \widehat{\text{\boldmathq\unboldmath}}_{e}^{h}, is defined as

[TABLE]

where $\mathbb{P}_{0}$ denotes the projection operator over the space of constant functions [24, 25] and $\tau_{e}$ is a stabilisation parameter that is selected to ensure stability, accuracy and convergence of the resulting scheme [8, 4, 21, 22, 20, 23].

Remark 2**.**

The definition of the numerical flux considered in Equation (10) follows the rationale of the hybridised DG method with reduced stabilisation [24]. If the projection operator is not considered in Equation (10), the second-order convergence is lost and the resulting method is only first-order [24], providing no advantages with respect to the original FCFV.

Introducing the expression of the numerical flux in Equation (9) and integrating by parts the first term, leads to the following discrete weak formulation: find (u_{e}^{h},\text{\boldmathq\unboldmath}_{e}^{h})\in\mathcal{V}^{1}(\Omega_{e})\times[\mathcal{V}^{0}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}} such that

[TABLE]

for all $v\in\mathcal{V}^{1}(\Omega_{e})$ and for $e=1,\dotsc,\texttt{n}_{\texttt{el}}$.

The discrete weak form of the global problem is obtained following an analogous process. It becomes: find $\hat{u}^{h}\in\mathcal{\hat{V}}^{0}(\Gamma\cup\Gamma_{N})$ such that

[TABLE]

a constant test function has been arbitrarily chosen in the space $\mathcal{\hat{V}}^{0}(\Gamma\cup\Gamma_{N})$.

Introducing the definition of the numerical flux in Equation (13) leads to the following weak form of the global problem: find $\hat{u}^{h}\in\mathcal{\hat{V}}^{0}(\Gamma\cup\Gamma_{N})$ such that

[TABLE]

2.4 Second-order FCFV discretisation

To simplify the notation, the following sets of faces are introduced

[TABLE]

corresponding to all faces of a cell, the faces on the Dirichlet boundary, the faces on the Neumann boundary and the faces not on the Dirichlet boundary, respectively. It is also convenient to denote the set of nodes of a cell $\Omega_{e}$ belonging to a face $\Gamma_{e,j}$ as $\mathcal{F}_{e,j}$. Finally, the indicator function of a set $\square$ is defined as

[TABLE]

With this notation, the discrete local problem becomes

[TABLE]

where $\tau_{j}$ denotes the value of the stabilisation parameter on the $j$-th face, assumed constant, $\mathbf{q}_{e}$ contains the value of $q$ in the cell and $\mathbf{u}_{e}$ contains the nodal values of the solution in the cell. The matrices and vectors in Equation (17b) are given by

[TABLE]

where $\texttt{n}_{\texttt{en}}$ is the number of cell nodes, the vector $\mathbf{p}_{e,i}$ is introduced to compute the projection of the solution, i.e. the average of the nodal values of $\mathbf{u}_{e}$ on face $\Gamma_{e,i}$. Formally it is defined as

[TABLE]

with $\texttt{n}_{\texttt{fn}}^{e,i}$ being the number of nodes of the face $\Gamma_{e,i}$.

The discrete local problem allows to obtain an explicit expression of both the solution and its gradient in terms of the solution at the cell faces/edges, namely

[TABLE]

where

[TABLE]

It is worth noting that the Equation (21b) involves the solution of a $3\times 3$ system of equations for triangular cells and a $4\times 4$ system for tetrahedral cells. Given the size of the system and the definition of $\mathbf{m}_{e}$, its inverse can be analytically computed to substantially reduce the computational cost of this operation.

The discretisation of the global problem of Equation (14) leads, for $i\in\mathcal{B}_{e}$, to

[TABLE]

By inserting the explicit expressions of Equation (21), in the global problem of Equation (23), a linear system of equations involving only the solution on the faces as an unknown is obtained, namely

[TABLE]

The global matrix $\mathbf{\widehat{K}}$ and right hand side $\mathbf{\hat{f}}$ are obtained by assembling the in the cell contributions given by

[TABLE]

for $i,j\in\mathcal{B}_{e}$ and with $\delta_{ij}$ denoting the Kronecker delta.

As discussed in Remark 2, the use of the projection operator in the numerical flux of Equation (10) is required to obtain a second-order method. If the projection is not introduced, the resulting method is only first-order [24]. This minor difference is exploited in Section 4 to devise an error indicator that can be used to drive an automatic mesh adaptive process. Therefore, it is of interest here to study the difference in the global system of Equation 25 induced by the introduction of the projection operator. The next result shows that only the matrix $\mathbf{m}_{e}$ changes if the projection is not considered.

Lemma 1**.**

Let us consider a face $\Gamma_{e,i}$ of a triangular or tetrahedral cell $\Omega_{e}$ and $u_{e}^{h}\in\mathcal{V}^{1}(\Omega_{e})$. Then, the following equality holds

[TABLE]

Proof.

The first integral of Equation (26) can be written as

[TABLE]

because $\mathbb{P}_{0}u_{e}^{h}$ is constant within each face.

The second integral of Equation (26) can be easily computed using the expression of the linear shape functions used to define the approximation $u_{e}^{h}$, namely

[TABLE]

where $N_{l}$ and $u_{l}$ denote the linear shape function and nodal value associated with the $l$-th node.

Equality (26) follows from the definition of $\mathbf{p}_{e,i}$ introduced in Equation (20). ∎

3 Second-order FCFV for the Stokes equation

3.1 Problem statement

The strong form of the Stokes problem can be written in the partitioned domain and after splitting the second-order momentum conservation equation into a system of two first-order equations, as

[TABLE]

where $\nu>0$ is the viscosity and the last two equations enforce the continuity of the velocity and the normal flux across the interface $\Gamma$ respectively.

Remark 3**.**

To simplify the presentation, this work considers the traditional velocity-pressure HDG formulation of the Stokes equation [32], where the vector $t$ does not correspond to the boundary traction and it is usually called a pseudo-traction [13]. It is worth emphasising that the so-called Cauchy formulation could also be employed here, using the formulation proposed in [16]. This formulation, contrary to other existing HDG methods, guarantees optimal convergence even for low order approximations. This idea was also exploited in the context of linear elasticity to obtain optimal convergence for low order approximations in an HDG context [30, 29, 27].

3.2 Strong form of the local and global problems

As usually done in HDG methods [9, 25, 16] and FCFV methods, the strong form of the problem given by Equation (29) is split into the local and global problems. The local problem

[TABLE]

is a pure Dirichlet problem and, therefore, requires the introduction of a solvability constraint for the pressure [9], namely

[TABLE]

where $\rho_{e}$ is the mean value of the pressure on the boundary of the cell $\Omega_{e}$.

In addition, a compatibility condition is induced by the free-divergence character of the velocity, namely

[TABLE]

As the continuity of the solution is automatically imposed by the introduction of the velocity on the cell faces as an independent variable that is uniquely defined on each face, the global problem imposes the continuity of the normal flux across the interface and the Neumann boundary conditions, that is

[TABLE]

3.3 Second-order FCFV weak formulation

Following the same rationale presented for the Poisson problem, the discrete weak formulation of the local Stokes problem is: find (\text{\boldmathu\unboldmath}_{e}^{h},\text{\boldmathL\unboldmath}_{e}^{h},p_{e}^{h})\in[\mathcal{V}^{1}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}}\times[\mathcal{V}^{0}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}\times\texttt{n}_{\texttt{sd}}}\times[\mathcal{V}^{0}(\Omega_{e})] such that

[TABLE]

for all \text{\boldmathw\unboldmath}\in[\mathcal{V}^{1}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}}. It is worth emphasising that in Equations (34a) and (34c), a constant test function has been arbitrarily chosen in the spaces $[\mathcal{V}^{0}(\Omega_{e})]^{\texttt{n}_{\texttt{sd}}\times\texttt{n}_{\texttt{sd}}}$ and $\mathcal{V}^{0}(\Omega_{e})$ respectively. It is also worth noting that Equation (34c) is exactly the compatibility condition introduced in Equation (32). As done in the standard FCFV method, this equation is then removed from the local problem and imposed only in the global problem.

The weak form of the local problem has been introduced after using the following definition of the numerical flux

[TABLE]

The discrete global problem that accounts for the transmission conditions, the Neumann boundary condition and the incompressibility constraint reads: find (\widehat{\text{\boldmathu\unboldmath}}^{h},\rho_{e})\in[\mathcal{\hat{V}}^{0}(\Gamma\cup\Gamma_{N})]^{\texttt{n}_{\texttt{sd}}}\times\mathbb{R} such that

[TABLE]

It is worth noting that a constant test function in the space $[\mathcal{\hat{V}}^{0}(\Gamma\cup\Gamma_{N})]^{\texttt{n}_{\texttt{sd}}}$ has been used in Equation (36a).

3.4 Second-order FCFV discretisation

Using the notation introduced in Section 2.4, the discrete local Stokes problem provides explicit expressions of the velocity, its gradient and the pressure in terms of the velocity on the faces and the mean pressure. The expressions are

[TABLE]

for $e=1,\dotsc,\texttt{n}_{\texttt{el}}$, where

[TABLE]

and

[TABLE]

Analogously, the discretisation of the global problem of Equation (36) after inserting the expressions of Equation (37), leads to a system of equations

[TABLE]

where the global matrix and right hand side are obtained by assembling the elemental contributions given by

[TABLE]

for $i,j\in\mathcal{B}_{e}$. The matrix $\mathbf{P}_{e,i}$, introduced to account for the projection of the solution on a space of constant functions on the face $\Gamma_{e,i}$, is defined as

[TABLE]

and

[TABLE]

in two and three dimensions respectively. The vector $\mathbf{p}_{e,i}$ was introduced in Equation (20).

4 Mesh adaptivity

As mentioned in Section 2.3, the proposed FCFV can be seen as a particular case of the hybridised DG method with reduced stabilisation [24] and therefore provides second-order convergence for the solution. The key aspect is the projection of the solution over a space of constant functions. Without this projection the method is only first-order accurate. This small difference in the formulation is exploited here to devise an error indicator. Noting $u$ the solution of the proposed FCFV methodology and $\tilde{u}$ the solution of the method where the projection is not performed, the following error indicator is proposed for the Poisson problem

[TABLE]

with a similar definition for the Stokes problem.

To compute the desired cell size, the error indicator of Equation (45) is combined with the a priori local error estimate for elliptic problems [10] given by

[TABLE]

for a constant degree of approximation. By using Richardson extrapolation, the following desired cell size is computed, for a desired error $\varepsilon$,

[TABLE]

where $h_{e}$ is the characteristic cell size used to perform the computation of both $u$ and $\tilde{u}$.

It is worth noting that the difference between $u$ and $\tilde{u}$ is only due to the use of the projection of the solution over a space of constant functions. Therefore, the majority of the calculations required to assemble the global system of equations can be re-used, substantially reducing the computational effort required to compute the error indicator of Equation (45). In fact, as detailed in Remark 1, the only difference between both formulations is in the matrix $\mathbf{m}_{e}$ in Equation (25) and $\mathbf{M}_{e}$ in Equation (42) for the Poisson and Stokes problems respectively.

5 Numerical studies

This Section presents a series of numerical experiments designed to test the optimal convergence properties of the proposed technique and to compare its performance with the recently proposed first-order FCFV [28]. Numerical experiments are also presented to illustrate the accuracy of the method in terms of the stabilisation parameter and the distortion and the stretching of the meshes. A detailed description of the model problems considered for both Poisson and Stokes problems in two and three dimensions can be found in [28] and they are omitted here for brevity.

5.1 Optimal convergence of the second-order FCFV scheme for Poisson equation

A mesh convergence study is performed for the Poisson problem using a series of successively refined triangular meshes. Figure 1 shows the relative $\mathcal{L}_{2}(\Omega)$ norm of the error of the solution and its gradient as a function of the characteristic mesh size $h$.

The results clearly show the optimal second-order convergence of the error of the solution and the first-order convergence of the error of the solution gradient. The results of the original FCFV are also included, clearly showing the gain in accuracy achieved for the solution $u$ for a given spatial discretisation. The gradient of the solution is only marginally more accurate as the approximation space for this variable is not changed with respect to the original FCFV.

The results of the analogous study in three dimensions are shown in Figure 2, demonstrating the optimal convergence of the error for both the solution and its gradient as well as the increased accuracy with respect to the first-order FCFV.

5.2 Optimal convergence of the second-order FCFV scheme for Stokes equation

A mesh convergence is next performed for the Stokes problem in two and three dimensions. Figures 3 and 4 show relative $\mathcal{L}_{2}(\Omega)$ norm of the error of the velocity, its gradient and the pressure as a function of the characteristic mesh size $h$.

The error of the velocity converges with second-order accuracy whereas first-order convergence is observed for its gradient, with an important gain on accuracy in the pressure. For the three dimensional test case, the optimal convergence is again observed for all the variables.

5.3 Computational cost

The convergence studies performed in the previous section show an important gain in accuracy of the proposed second-order FCFV method when compared to the original FCFV in the same mesh. In this section both methods are compared in terms of the computational time.

Figure 5 shows relative $\mathcal{L}_{2}(\Omega)$ norm of the error of the velocity, its gradient and the pressure as a function of the CPU time for the two dimensional Stokes problem.

The results show that the proposed second-order FCFV is able to produce more accurate results for the velocity and the pressure using the same CPU time when compared to the first-order FCFV, whereas similar results are obtained for the veloocity gradient. In three dimensions similar conclusions are obtained, as shown in Figure 6.

The results show that the proposed second-order FCFV provides the same accuracy as the original first-order FCFV with orders of magnitude less CPU time when the velocity is of interest. For instance, an error in the velocity field of the order of 1% is obtained in less than 1 second with the second-order FCFV whereas the first-order FCFV requires 2.7 hours.

The results for the Poisson problem, not displayed here for brevity, show the same advantages for the proposed second-order FCFV.

5.4 Influence of the stabilisation parameter

The next study considers the influence of the stabilisation parameter $\tau$ in the accuracy of the proposed second-order FCFV method. Figure 7 shows the relative error, measured in the $\mathcal{L}_{2}(\Omega)$ norm, as a function of the stabilisation parameter for the two and three dimensional Poisson problem and for two levels of mesh refinement.

The results show that a maximum accuracy in the solution is achieved for a value of the stabilisation parameter of $10^{2}$ or larger, whereas the error of the solution gradient seems insensitive to the choice of this parameter. The behaviour is almost identical in both two and three dimensions.

For the Stokes problem, similar conclusions are obtained, as illustrated in Figure 8.

In this case, the velocity gradient and the pressure are insensitive to the stabilisation parameter $\tau$, whereas the velocity requires a value of $10^{2}$ or larger to provide the maximum accuracy.

It is worth emphasising that the value required to achieve the maximum accuracy of the solution is the same for two and three dimensional problems and for Poisson and Stokes problems. In addition, compared to the results presented in [28] for the first-order FCFV, the proposed second-order FCFV is less sensitive to a particular choice of the stabilisation parameter.

5.5 Influence of the cell distortion and stretching

The last study considers the solution of Poisson and Stokes problems in meshes involving distorted and stretched cells. To illustrate the type of cells tested, Figure 9 shows the mesh corresponding to the third level of refinement where the cells have been randomnly disorted, as explained in [28] and with stretched cells with a stretching factor of 100.

Figure 10 shows a mesh convergence study for the Poisson problem in two and three dimensional meshes that have been distorted by randomly moving the interior nodes (i.e. the nodes not on the boundary of the domain) as described in [28].

The results are almost identical to the ones obtained for regular meshes and displayed in Figures 1 and 2, showing that the proposed method is insensitive to mesh distortion.

The same conclusions are also obtained for the Stokes problem in two and three dimensions as shown in Figure 11.

Next, the influence of the cell stretching is considered. Two cases involving a maximum stretching factor, $s$, of 10 and 1,000 are considered, where this factor is measured in each cell as the ratio between the longest and the shortest edges. Figure 12 shows the relative $\mathcal{L}_{2}(\Omega)$ norm of the error of the solution and its gradient as a function of the characteristic mesh size for the Poisson problem solved in stretched meshes in two and three dimensions.

Almost identical results are observed for both stretching factors. In the two dimensional problem a marginal lower error is observed for the mesh with stretching factor $s=10$ whereas for the three dimensional problem the mesh with stretching factor $s=1,000$ produces a slightly lower error.

The results for the Stokes problem in two and three dimensions are displayed in Figure 13, confirming the conclusions observed for the Poisson problem.

It is worth noting that contrary to other FV methods, the proposed second-order FCFV method not only shows an accuracy that is insensitive to mesh distortion and stretching but also preserves the optimal rate of convergence in all the variables, i.e. the solution and its gradient for the Poisson problem and the velocity, its gradient and the pressure for the Stokes problem.

6 Numerical examples

6.1 Irrotational flow past a full aircraft

To show the ability of the proposed method to efficiently solve large scale problems involving complex geometries, the irrotational flow around a full aircraft is considered. A tetrahedral mesh with 5,125,998 cells is considered, leading to a global system of 11,283,113 equations to find the solution on the cell faces. The magnitude of the velocity, computed from the gradient of the solution, and the pressure, computed from the Bernoulli equation are displayed in Figure 14.

The solution using the proposed second-order FCFV took 5.1 minutes for the computation of all the elemental matrices and the assembly and 6.4 minutes for solving the global system of equations using a direct method. The developed code is written in Matlab and the computation was performed in an Intel*®* Xeon*®* CPU $@$ 3.70GHz and 32GB main memory available.

It is worth noting that the time recorded for assembling the system is slightly higher than the times reported in [28] for the first-order FCFV due to the extra operations required by the second-order method for the computation and assembly of the global matrix. The time required by the proposed second-order FCFV for solving the system is almost identical to the time employed by the first-order method due to the global matrix having the exact same size and structure, with the same number of non-zero elements. Finally, it is worth emphasising that, as shown in the numerical experiments of the previous Section, the extra cost induced by the second-order FCFV for a given spatial discretisation leads to a substantial gain in accuracy and, in general, for a given error, the second-order FCFV is more efficient.

6.2 Stokes flow past a sphere

The next example considers the Stokes flow around a sphere. This classical three dimensional example is used to compare the accuracy and performance of the proposed second-order FCFV against the original first-order FCFV method for a problem involving large three dimensional meshes. The domain of interest is $\Omega=\left([-7,15]\times[-5,5]\times[-5,5]\right)\setminus\mathcal{B}_{1,\mathbf{0}}$, where $\mathcal{B}_{1,\mathbf{0}}$ is the unit ball with centre at the origin. Homogeneous Neumann boundary conditions are imposed on the outlet part of the boundary, corresponding to $x=15$, whereas Dirichlet boundary conditions, corresponding to the exact solution, are imposed on the rest of the boundary.

Six unstructured tetrahedral meshes are considered with 3,107, 10,680, 43,682, 204,099, 686,853 and 2,516,099 cells respectively. The size of the global system of equations to be solved to obtain the velocity on the cell faces is 20,711, 72,249, 299,276, 1,409,916, 4,765,776 and 17,513,075 respectively.

Figure 15(a) shows the convergence of the error for the velocity, pressure and gradient of the velocity as the mesh is refined. The results are compared to the original first-order FCFV and clearly show the advantage of the proposed method by providing second-order convergence on the velocity. In this example the accuracy for the other variables is almost identical, with a marginal gain observed in the computation of the pressure in the finest meshes.

Figure 15(b) shows the convergence of the drag as the mesh is refined. The advantages of the proposed second-order FCFV are observed as the convergence to the exact value is faster than with the original first-order method.

6.3 Mesh adaptivity for the Poisson problem

This Section presents a numerical example to illustrate the strategy described in Section 4 to perform an automatic mesh adaptive process by combining the results of the first-order and second-order FCFV methods. A two dimensional Poisson problem with known analytical solution is considered in $\Omega=[0,1]^{2}$. The source term and Dirichlet boundary conditions are selected so that the analytical solution is given by

[TABLE]

and the desired accuracy is $\varepsilon=10^{-2}$.

The process starts with the coarse mesh represented in Figure 16(a). The approximation with the proposed FCFV on the coarsest mesh is depicted in Figure 17(a).

By comparing the approximation of the second-order FCFV with the approximation of the first-order FCFV method (i.e. without the projection operator), a desired cell size is computed, using Equation (47), and a new mesh is generated. The mesh and the second-order FCFV approximation after two iterations of the mesh adaptive process are displayed in Figures 16(b) and 17(b) respectively. It can be clearly observed how the mesh is coarsened in the regions where the approximation is almost constant, whereas the mesh density is increased in the regions where the approximation changes rapidly. The adaptive process finishes in five iterations, when the desired error is achieved. The final mesh and FCFV approximation are represented in Figures 16(c) and 17(c) respectively.

As described in Section 4, the adaptive process is driven by an error indicator that results from computing the relative difference in each cell, measured in the $\mathcal{L}_{2}(\Omega_{e})$ norm. To illustrate the efficiency of the proposed error indicator, Figure 18 shows both the error indicator and the exact error for the three iterations, $\texttt{n}_{\texttt{i}}$, of the adaptive process corresponding to the meshes and approximations shown in Figures 16 and 17.

To further illustrate the performance of the error indicator and the automatic adaptive process, Figure 19(a) shows the maximum value of the indicator and the exact error over all the cells as a function of the number of iterations of the mesh adaptive process.

The efficiency of the error indicator, defined as the ratio between the exact error and the indicator is also displayed in Figure 19(b), clearly illustrating the suitability of the proposed technique to drive an automatic mesh adaptive process.

6.4 Mesh adaptivity for the Stokes problem

The last example involves the solution of the Stokes equations in three dimensions for the complex geometry, taken from [34], depicted in Figure 20.

The corrugated channel has a height of 0.5$\mu$m and the curved profile is obtained by repeating the expression

[TABLE]

where $f_{n}=0.5\mu$m, $f_{\omega}=4.5\mu$m and $L=12.5\mu$m. Two spheres of radius 0.2$\mu$m are placed inside the channel. The first sphere, with centre $(0,0,0.25)\mu$m, is placed in the middle of the channel, where the cross section is minimum, and it is expected to produce a major flow disturbance. The second sphere, with centre at $(-3.75,1,0.25)\mu$m, is situated in a region with larger cross section and it is expected to disturb much less the flow. This scenario is utilised to show the ability of the proposed method to drive the adaptivity for a problem involving an incompressible flow in a complex geometry.

A Dirichlet boundary condition, corresponding to a velocity inlet given by \text{\boldmathu\unboldmath}_{D}(x,y,z)=64(y^{2}-1/4)(z^{2}-1/16), is introduced at one end of the channel, at $x=-L$, depicted in red in Figure 20. A homogeneous Neumann boundary condition is imposed at the outlet, at $x=L$, depicted in blue in Figure 20. Homogeneous Dirichlet boundary conditions are imposed on the rest of the boundary, corresponding to material walls.

The initial mesh, shown in Figure 21 (a), has 37,415 tetrahedral cells. The initial mesh is generated with a required element size of $0.05\mu$m on the surfaces defining two spheres, to ensure an appropriate geometric representation. An element size of $0.5\mu$m is imposed in the domain, with a smooth transition between these two values.

Next, the automatic adaptive process described in Section 4 is applied, with a desired accuracy of $\varepsilon=5\times 10^{-2}$. Convergence of the estimated error to the desired tolerance is achieved in this example after seven iterations of the adaptive process. The resulting meshes in the fourth and seven iteration are depicted in Figures 21 (b) and (c) respectively. The fourth mesh has 61,871 cells, with a minimum and maximum element size of $0.010\mu$m and $0.594\mu$m respectively. The last mesh has 116,913 cells, with a minimum and maximum element size of $0.003\mu$m and $0.520\mu$m respectively. A detailed view of the three meshes of Figure 21 near the two spheres is shown in Figure 22.

The velocity streamlines and the pressure field obtained using the last mesh are represented in Figure 23.

A detailed view of the velocity streamlines and the pressure field near the sphere surfaces is shown in Figure 24. This Figure offers a qualitative comparison of the major disturbance caused by one sphere compared to the other.

7 Concluding remarks

This paper proposes a second-order FCFV method for the solution of scalar and vector elliptic problems in two and three dimensions.

The proposed method preserves the attractive properties of the original first-order FCFV method, namely the first-order convergence of the gradient of the solution, without the need of a reconstruction, and the insensitivity to mesh distortion and stretching. It also satisfies the LBB condition in the context of incompressible flows. Contrary to the original FCFV, the proposed method guarantees second-order convergence of the solution. Numerical experiments show an increased performance when compared to the first-order method in terms of the CPU time required to achieve a desired accuracy as well as a lower sensitivity to the choice of the stabilisation parameter. A combination of first-order and second-order schemes is used to devise an error indicator that can be used to drive a mesh adaptivity process. An extensive set of numerical experiments has been used to demonstrate the optimal approximation properties of the method and more complex problems demonstrate its potential for large scale three dimensional simulations, including a Stokes problem where an automatic mesh adaptive process is employed.

Acknowledgements

This work was partially supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Actions (Grant numbers: 675919 and 764636). The first author gratefully acknowledges the financial support of the Erasmus+ programme. The second, third and last author were also supported by the Spanish Ministry of Economy and Competitiveness (Grant agreement No. DPI2017-85139-C2-2-R). The second and last authors are grateful for the financial support provided by the Generalitat de Catalunya (Grant agreement No. 2017-SGR-1278).