A non-oscillatory face-centred finite volume method for compressible flows

This work presents the face-centred finite volume (FCFV) paradigm for the simulation of compressible flows. The FCFV method defines the unknowns at the face barycentre and uses a hybridisation procedure to eliminate all the degrees of freedom inside the cells. In addition, Riemann solvers are defined implicitly within the expressions of the numerical fluxes. The resulting methodology provides first-order accurate approximations of the conservative quantities, i.e. density, momentum and energy, as well as of the viscous stress tensor and of the heat flux, without the need of any gradient reconstruction procedure. Hence, the FCFV solver preserves the accuracy of the approximation in presence of distorted and highly stretched cells, providing a solver insensitive to mesh quality. In addition, FCFV is capable of constructing non-oscillatory approximations of sharp discontinuities without resorting to shock capturing or limiting techniques. For flows at low Mach number, the method is robust and is capable of computing accurate solutions in the incompressible limit without the need of introducing specific pressure correction strategies. A set of 2D and 3D benchmarks of external flows is presented to validate the methodology in different flow regimes, from inviscid to viscous laminar flows, from transonic to subsonic incompressible flows, demonstrating its potential to handle compressible flows in realistic scenarios.


Introduction
Finite volume (FV) solvers are the most widespread technology within the aerospace community for the simulation of steady-state compressible flows [1,2]. The success of such methodologies is mainly due to their capability of providing results for large-scale flow problems involving complex geometries by means of overnight simulations. For this reason, FV implementations are accessible in many computational fluid dynamics (CFD) platforms, from open-source to commercial and industrial software [3,4,5,6,7,8,9,10,11]. Nonetheless, these methods require delicate mesh generation procedures, limiting unstructured regions and distorted cells [12,13], in order to construct high-quality meshes suitable for computation.
Traditional FV solvers, see e.g. [2,14,15,16], rely on two paradigms, the so-called cellcentred finite volume (CCFV) method [17] and the vertex-centred finite volume (VCFV) method [18]. The former defines the unknown solution at the centroid of the cells, whereas the latter at the mesh nodes. Recently, the face-centred finite volume (FCFV) paradigm was proposed for a series of linear elliptic partial differential equations (PDEs) [19,20,21,22]. The FCFV method utilises a mixed FV formulation and eliminates the degrees of freedom within each cell via a hybridisation step, leading to a problem defined in terms of the unknowns at the face barycentres only. Finally, the variables inside each cell are retrieved via a computationally inexpensive postprocessing step performed independently cell-by-cell.
The present work discusses the first FCFV formulation for nonlinear hyperbolic PDEs modelling steady-state compressible flows, spanning from viscous compressible Navier-Stokes to inviscid Euler equations. The proposed approach introduces the flow variables, i.e. density, momentum and energy, the deviatoric strain rate tensor and the gradient of temperature as unknowns inside each cell, whereas the conservative quantities on the cell faces define the so-called hybrid unknowns. The resulting method yields first-order accuracy for all above mentioned variables. Contrary to traditional first-order CCFV and VCFV paradigms, the FCFV method provides first-order convergence of the deviatoric strain rate tensor and of the gradient of temperature. In addition, this is achieved without resorting to any flux reconstruction strategy required by traditional second-order CCFV and VCFV approaches. This is especially important as it allows the FCFV method to preserve optimal accuracy in the computation of the stress tensor and the heat flux even in presence of distorted and stretched cells [19,20,21,22], where traditional second-order CCFV and VCFV approaches are known to suffer from a loss of accuracy [12,13,23]. In addition, the FCFV method showed its versatility in devising approximations based on meshes of different cell types and on hybrid meshes [19,21], proving to be a robust methodology insensitive to element type and mesh quality. Although the FCFV method is known to introduce a larger number of unknowns than CCFV and VCFV approaches [19], the simplified procedure required to generate meshes suitable for computation and its capability to avoid gradient reconstruction make the FCFV scheme a competitive alternative to traditional FV solvers from the computational viewpoint.
In the approximation of conservation laws, a critical aspect is represented by the strategy employed to handle the nonlinear convective fluxes. In this context, Riemann solvers have been extensively studied [24,25,26,27]. The present work revisits the traditional Lax-Friedrichs, Roe, HLL and HLLEM Riemann solvers in the context of the FCFV method. An extensive set of numerical examples is employed to validate and compare these approaches. Special attention is devoted to the development of positivity-preserving approximations in order to devise a robust methodology able to compute physically-admissible solutions across a wide range of flow regimes. More precisely, numerical experiments are reported to showcase the capability of the FCFV method to construct non-oscillatory approximations of shock waves and sharp fronts, without the need of shock capturing techniques based on artificial viscosity [28,29,30,31] or flux limiters [32,33,34,35].
It is worth mentioning that FV methods have been reinterpreted in recent years as particular cases of finite element discretisations [14]: the CCFV method can be seen as a discontinuous Galerkin (DG) approach with piecewise constant approximations in each cell [16,36], whereas a VCFV discretisation on a simplicial mesh is equivalent to a continuous Galerkin finite element method with piecewise linear approximations [37,38]. Concerning the FCFV paradigm, it can be interpreted as the lowest-order version of the hybridisable discontinuous Galerkin (HDG) method [39,40,41], in which constant approximations are selected for all the variables. More precisely, the FCFV formulation presented in this work is inspired by the high-order HDG discretisations for compressible flows discussed in [42,43]. The proposed FCFV method thus inherits its properties from HDG, including its stability in the incompressible limit [19,20,21,22], circumventing the Ladyzhenskaya-Babuška-Brezzi (LBB) condition [29,44]. Hence, the FCFV method provides accurate approximations of incompressible flows without the need of introducing specific pressure corrections like the well-known SIMPLE algorithm [45] implemented by commercial and open-source software, see e.g. [3,9].
The remainder of this paper is organised as follows. In section 2, the compressible Navier-Stokes equations are recalled. Section 3 introduces the corresponding FCFV discretisation and the integration of the different Riemann solvers in the definitions of the numerical fluxes. A set of convergence tests to validate the optimal accuracy properties of the method for inviscid and viscous laminar flows is presented in section 4, with special emphasis on the robustness to cell stretching and distortion. Two-and three-dimensional benchmarks of external flows of aerodynamic interest are reported in section 5 to demonstrate the capabilities of the method in different flow regimes, from inviscid to viscous laminar flows, from transonic to subsonic incompressible flows. Finally, section 6 reviews the main results of this work and three appendices provide technical details on the definition of the boundary conditions (A), on the construction of the mixed variable (B) and on the hybridisation procedure (C) in the FCFV solver.

Governing equations for compressible flows
Consider the compressible Navier-Stokes equations, written in non-dimensional conservative form as where Ω ⊂ R nsd denotes an n sd -dimensional open bounded domain with boundary ∂Ω, T end > 0 is the final time of interest, U 0 represents the initial state and B is an operator imposing inlet, outlet or wall boundary conditions as detailed in A.
Remark 1 (Pseudo-time in steady-state flows). The present work focuses on the development of a novel FV spatial discretisation for steady-state compressible flows. In this context, t represents an artificial pseudo-time and time marching algorithms are introduced to speed-up the convergence of the nonlinear solver, as detailed in section 3.5.
The system (1) is written in terms of the vector of conserved quantities U ∈ R nsd+2 and the advection, F , and diffusion, G, flux tensors, namely where ρ is the density, v the velocity vector, E the total specific energy and p denotes the pressure field. The viscous stress tensor σ d and the heat flux vector q are given by where ∇ S := (∇ + ∇ T )/2 is the symmetric part of the gradient operator and T denotes the temperature. The non-dimensional dynamic viscosity is defined according to the Sutherland's law, i.e. µ = (T /T ∞ ) 3/2 (T ∞ + S)/(T + S), being the non-dimensional free-stream temperature and the Sutherland constant, respectively, with S 0 = 110K and T ref = 273K.
In addition, following the assumption of a calorically perfect gas, the pressure is given by p = (γ − 1)ρ (E − v 2 /2), where γ = c p /c v (equal to 1.4 for air) is the ratio of specific heats at constant pressure, c p , and constant volume, c v . Finally, from the ideal gas law it follows that γp = (γ − 1)ρT .
Remark 2 (Non-dimensional quantities). The problem is described in terms of the Mach, Reynolds and Prandtl numbers, defined respectively as c = γp/ρ being the speed of sound, L a characteristic length and κ the thermal conductivity. The subscript ∞ denotes free-stream reference values.

FCFV formulation for compressible flows
The domain Ω is partitioned in a set of n el disjoint cells Ω e satisfying Ω = nel e=1 Ω e . In addition, the boundary of the cell Ω e , denoted by ∂Ω e , is obtained as the union of its faces, namely where n e fa is the total number of faces of the cell Ω e . Finally, the internal interface Γ is defined as

Introducing a set of mixed variables
Following the HDG [39,44,42] and FCFV [19,20] rationales, the second-order problem (1) is written as a system of first-order PDEs via the introduction of a set of socalled mixed variables. In the context of compressible flows, the most common approach, see [43,42,46,47], relies on defining the mixed variable as the gradient of the so-called primal variable U , namely Ψ = ∇U . Nonetheless, this choice leads to the introduction of a mixed variable, associated with the gradient of density, which is redundant since the mass conservation equation is a first-order PDE. In addition, several nonlinearities appear in the resulting expressions to compute the viscous stress and the heat flux starting from U and Ψ. Following [48], in this work only two mixed variables are considered, namely the deviatoric strain rate tensor ε d and the gradient of temperature φ, given by Remark 3 (Deviatoric strain rate). It is worth noticing that the deviatoric strain rate tensor can be expressed as a function of the gradient of velocity as ε d = D∇ S v, where the linear operator D is defined as Interested readers are referred to B for the details concerning the construction of the operator D and its implementation.
Besides reducing the number of mixed variables involved, from (7) it also follows that the viscous stress tensor and the heat flux vector in equation (3) can be obtained using the linear expressions

A mixed hybrid finite volume framework
The FCFV method solves the compressible Navier-Stokes equations in two stages. First, an independent hybrid variable U , representing the vector of conservative variables on the mesh faces Γ ∪ ∂Ω, is introduced. Equation (1) is thus rewritten in each cell Ω e , e = 1, . . . , n el as Equation ( Second, the vector U of conservative variables on the faces is computed by solving the FCFV global problem, which prescribes the continuity of the conservative variables and of the normal fluxes on Γ and the boundary conditions on ∂Ω, namely where n stands for the outward normal vector to the cell face and = + + − denotes the jump operator defined on an internal face as the sum of the values in the neighbouring elements Ω + and Ω − , respectively [49]. The trace boundary operator B(U , U , ε d , φ) imposes the boundary conditions on ∂Ω exploiting the hybrid variable, as detailed in A.
It is worth noticing that the first condition in equation (11) is automatically satisfied due to the Dirichlet boundary conditions imposed in the local problems (10) and because of the unique definition of the hybrid variable U on each face.

Integral form of the FCFV local and global problems
For each cell Ω e , e = 1, . . . , n el , the integral form of the FCFV local problem is obtained by applying the divergence theorem to equation (10). Given U = U 0 at time t = 0, it holds that where v and T denote the velocity and temperature fields on the cell faces ∂Ω e , respectively, and they are defined using the hybrid vector U of conservative variables. This problem corresponds to the hybridisation step of the FCFV method: the goal is to eliminate the unknowns (U , ε d , φ) within each cell by expressing them in terms of the hybrid variable U .
The unknown U is thus computed by means of the global problem (11) whose integral form is The terms F (U )n and G(U , ε d , φ)n appearing in equations (12c) and (13) stand for the convection and diffusion numerical fluxes of the conservation equations, respectively. It is worth recalling that the approximation of the numerical fluxes in the FCFV method plays a crucial role in the accuracy and stability of the computed solution, see [19,20,22,21]. Following the rationale discussed for high-order HDG discretisations in [48,43,50,42,47], the traces of the numerical fluxes on the cell faces are defined as On the one hand, the stabilisation tensor τ a is associated with convection phenomena. Different expressions of τ a are derived from the theory of Riemann solvers for nonlinear hyperbolic PDEs, as described in section 3.4. On the other hand, the term τ d stands for the stabilisation tensor related to viscous effects and is defined by means of the diagonal matrix 1 nsd being an n sd -dimensional vector of ones.
Remark 4 (FCFV method for inviscid Euler equations). The inviscid Euler equations are obtained as the limit of the compressible Navier-Stokes equations when Re → ∞. Setting G = 0, equation (1) thus reduces to the well-known system of first-order hyperbolic PDEs modelling inviscid compressible flows. For each cell Ω e , e = 1, . . . , n el , the corresponding FCFV local problem for the Euler equations is obtained from equation (12) by neglecting the mixed variables and the viscous term, namely Similarly, the global problem follows from the simplification of equation (13) as

Riemann solvers for the FCFV method
In the context of the FCFV method, Riemann solvers are defined implicitly within the convection fluxes, see equation (14a), by means of appropriate expressions of the stabilisation term τ a . In this section, the stabilisation terms leading to the formulation of the Lax-Friedrichs, Roe, HLL and HLLEM numerical fluxes are presented.
· n be the Jacobian of the convective fluxes along the normal direction to a cell face. Moreover, denote by λ max := | v · n| + c the maximum eigenvalue in absolute value of the matrix A n ( U ), v and c being the velocity and the speed of sound evaluated from U , respectively. Following the unified formulation in [48], the Riemann solvers for the FCFV method are detailed below.

Lax-Friedrichs Riemann solver
The FCFV stabilisation tensor inspired by the Lax-Friedriches Riemann solver, see [24], is defined as

Roe Riemann solver
Consider the spectral decomposition A n ( U ) = RΛR −1 , where Λ is a diagonal matrix containing the eigenvalues λ i , i = 1, . . . , n sd +2 of A n ( U ) and R is the corresponding matrix of right eigenvectors. In addition, the diagonal matrix Φ is given by Φ = diag (ϕ 1 , . . . , ϕ nsd+2 ) where ϕ i = max(|λ i |, δ), δ ≥ 0 being a user-defined parameter. The Roe Riemann solver with Harten-Hyman entropy fix [51] is obtained for the FCFV method by setting The parameter δ represents the threshold value of the aforementioned Harten-Hyman entropy fix. Such correction aims to remedy the failure of entropy conditions of the Roe solver, which may produce nonphysical solutions in transonic and supersonic cases. It is worth noticing that for δ = 0, Φ = |Λ| and the traditional Roe Riemann solver is retrieved, namely τ a = |A n ( U )|.

HLL Riemann solver
The HLL Riemann solver is devised to recover the Rankine-Hugoniot condition, for a simplified scenario in which contact discontinuities are neglected, without the need of any user-defined parameter [52]. The resulting positivity-preserving Riemann solver for the FCFV method is given by where s + := max(0, v · n + c) is an estimate of the largest wave speed of the Riemann problem.

HLLEM Riemann solver
In order to exploit both the positivity-preserving properties of the HLL Riemann solver and the capability of Roe's method to capture shear layers, the HLLEM Riemann solver [53,54] is employed to construct the stabilisation tensor where s + is the HLL wave speed estimate and θ( U ) = RΘR −1 replaces the identity I nsd+2 in equation (18c). Note that the definition of θ( U ) exploits the matrix of right eigenvectors arising from the spectral decomposition of A n ( U ), whereas the diagonal matrix Θ is given by

FCFV discrete problem
The discrete form of the FCFV method is obtained by introducing the definition (14) of the numerical fluxes into the local (12) and global (13) problems. In addition, the vector of conservative variables U and the mixed variables ε d and φ are discretised using a constant value at the centroid of each cell, whereas a constant approximation at the barycentre of the faces is employed for the hybrid vector U . Finally, a quadrature rule based on a single integration point is utilised to evaluate the integral quantities on cells and faces.
Moreover, χ Ie and χ Ee are defined to represent the indicator functions associated with the sets I e and E e , respectively.
The semi-discrete form of the FCFV local problem (12) is: for e = 1, . . . , n el , given the initial state U e = U 0 e at t = 0 and the hybrid vector U j on the faces Γ e,j , j = 1, . . . , n e fa , compute (U e , ε d e , φ e ) that satisfy Remark 5 (Symmetry of the mixed variable). The mixed variable ε d is a second-order symmetric tensor commonly represented using a matrix of dimension n sd × n sd . In order to exploit the symmetry property in its discretisation, Voigt notation is employed to store only its m sd = n sd (n sd + 1)/2 non-redundant components. This approach, detailed in B, was first proposed in the context of hybrid discretisation methods for high-order HDG formulations, see [56,57,58,59,60], and later exploited also for FCFV approaches in [20]. For the simulation of compressible and weakly-compressible flows, this approximation of the deviatoric strain rate tensor was employed in [48,61] in the context of high-order HDG methods.
Remark 6 (Time integration scheme). In order to obtain the fully-discrete form of the local problem (20c), an appropriate time integration scheme needs to be introduced. As previously mentioned, the present work proposes a novel FV spatial discretisations and the numerical examples in sections 4 and 5 focus on steady-state flows, whence this term is neglected. Nonetheless, it is worth mentioning that time marching based on an artificial time is a common relaxation approach to speed-up the convergence of a nonlinear solver. In this context, time derivative may be discretised using a backward Euler scheme, that is where ∆t is the artificial time step. Note that similar FCFV discrete problems are obtained using other implicit time integration schemes, e.g. higher-order backward difference formulae (BDF) [62,63,42], providing additional accuracy in the simulation of transient phenomena. Of course, in the case of transient simulations, Newton-Raphson iterations are performed at each time step to solve the nonlinear problem.
In a similar fashion, the discrete form of the FCFV global problem (13) is: compute It is worth noticing that both the local (20) and global (22) problems are nonlinear. More precisely, let Q e = (ε d e , φ e ) be the set of mixed variables introduced by the FCFV formulation in the cell Ω e . The resulting system of algebraic-differential equations arising from the local problem (20) is where U e and Q e are the vectors containing the values of the local and mixed variables, respectively, at the centroid of the cell, whereas the vector U collects the values of the hybrid variable at the barycentres of its faces. On the one hand, equations (20a) and (20b) provide analytical expressions of the mixed variables in terms of the hybrid unknown, see equation (23a). On the other hand, equation (20c) is nonlinear and the residual vector obtained from its spatial discretisation is denoted by R e . Upon linearisation of equation (23b) via a Newton-Raphson procedure, the n el local problems allow to express U e and Q e in each cell Ω e , e = 1, . . . , n el in terms of the unknown U on its faces. The resulting expressions are plugged into equation (22) and all the degrees of freedom inside the cells are eliminated from the global problem, leading to where R e is the nonlinear residual vector related to the unknowns U associated with cell Ω e . The global problem (24), whose structure is detailed in C, is thus solved by means of a Newton-Raphson linearisation.

Inviscid Ringleb flow
The convergence properties of the FCFV method in the inviscid limit are examined through the Ringleb flow problem [63,42]. This 2D example describes a smooth transonic flow, for which an analytical expression of the solution can be computed via the hodograph method [64]. At a given point (x, y), the solution is obtained as result of the nonlinear implicit equation where c is the speed of sound, whereas density ρ, velocity magnitude V , pressure p and J are determined as Finally, the velocity vectorfield is given by sgn (·) being the sign operator and sin(2θ) = 2ρV 2 y.
The computational domain is defined as Ω = [0, 1] 2 and far-field boundary conditions are imposed on ∂Ω. Figure 1 displays two levels of refinement of the domain using uniform meshes of triangular and quadrilateral cells. The corresponding approximation of the Mach number distribution on these meshes is reported in figure 2.
The relative error of the numerical approximation, measured in the L 2 (Ω) norm as function of the characteristic mesh size h, is examined for the four Riemann solvers discussed in section 3.4. The h-convergence study is performed using the sets of meshes

Viscous laminar Couette flow
A Couette flow with source term [65,42] is defined in the domain Ω = [0, 1] 2 to examine the convergence properties of the FCFV method in the viscous laminar regime. The analytical expressions of velocity, pressure and temperature, are where α c = 0.8, β c = 0.85 and the free-stream Mach number is set to M ∞ = 0.15. Assuming constant viscosity, the source term is determined from the exact solution and is given by Finally, the boundary conditions are prescribed on ∂Ω employing the expression of the analytical solution.
The h-convergence study is performed using the meshes of triangular cells in figure 1. Figure 4 reports the results using different Riemann solvers and for Reynolds number Re = 1 and Re = 100. The approximation of the conservative variables displays optimal convergence regardless of the Reynolds number and of the employed Riemann solver. The proposed method is thus able to provide optimal accuracy in the conservative quantities also in case viscous phenomena are considered, with errors between 10 −3 and 10 −5 . Regarding the viscous stress tensor and the heat flux, optimal accuracy is achieved using HLLEM and Roe Riemann solvers independently of the Reynolds number, whereas the Lax-Friedrichs and HLL fluxes appear to be more sensitive to increasing values of the Reynolds number.
Finally, it is worth mentioning that the test case under analysis features an incompressible flow (∇ · v = 0, M ∞ = 0.15). Despite this additional difficulty, the FCFV method is capable of computing an accurate approximation without the need of introducing specific pressure corrections like the well-known SIMPLE algorithm [45]. Thus, an important advantage of the proposed methodology is its robustness in the incompressible limit as futher detailed in section 5.3.

Influence of cell distortion and stretching
In this section, the sensitivity of the FCFV method to cell distortion and stretching is investigated. For the sake of brevity, this study focuses on the viscous case in order to  First, a set of highly distorted meshes is generated by introducing a perturbation on the position of the interior nodes of the regular meshes illustrated in figure 1. In particular, for a given node i, its new position is defined asx i = x i + r i , r i being an n sd -dimensional vector whose components are randomly generated within the interval [− min /3, min /3], where min denotes the minimum edge length of the regular mesh. The third and fifth level of refinement of the meshes featuring distorted quadrilateral and triangular cells are illustrated in figure 5.
The previously introduced Couette flow example for Reynolds number Re = 100 is also employed for the current sensitivity study. More precisely, figure 6 shows the approximation  Second, a set of meshes with stretching near the bottom boundary is produced. For its construction in 2D, the vertical coordinate of the first mesh layer is fixed at the desired stretching factor s. Then, the vertical coordinate of the subsequent layers is defined as where h is the maximum edge length of the corresponding regular mesh, N y is the number of cells in the vertical direction and the growth rate factor β is computed by imposing that the vertical coordinate of the last layer is one, that is by finding the roots of Figure 7 reports the second level of refinement of a set of triangular meshes for different levels of stretching s. A quantitative evaluation of the influence of cell distortion and stretching on the accuracy of the FCFV approximation is performed via an h-convergence study of the error, measured in the L 2 (Ω) norm, using the HLLEM Riemann solver. The results, reported in figures 8a and 8b, respectively, show that optimal convergence of order 1 is achieved for all the variables, independently of the distortion or of the stretching factor of its cells. In addition, the precision of the numerical approximation also results unaffected by the loss of orthogonality and loss of isotropy of the mesh. Indeed, by comparing the results of figure 8 with the ones in figure 4, almost identical levels of accuracy are obtained in the L 2 (Ω) error of the approximate solutions using meshes with uniform, distorted or stretched cells.

Numerical benchmarks
In this section, a set of numerical examples is presented to show the capabilities of the proposed FCFV method to simulate inviscid and viscous compressible flows at different regimes.

Inviscid transonic flow over a NACA 0012 profile
The first test case considers the inviscid transonic flow over a NACA 0012 aerofoil at freestream Mach number M ∞ = 0.8 and angle of attack α = 1.25 • . This classical benchmark for inviscid compressible flows [66,67] is proposed to evaluate the ability of the FCFV solver to capture flow solutions involving shock waves. More precisely, this benchmark is used to demonstrate the importance of the choice of the Riemann solver in the accuracy and stability of the FCFV approximate solution.
Unstructured meshes of triangular cells with non-uniform refinements on the surface of the aerofoil and at the leading and trailing edges are used for the simulation. Figure 9 reports the details of a coarse and a fine mesh featuring 89,250 and 712,164 triangular cells, respectively. The far-field boundary is located at 15 chord units away from the profile and the aerofoil surface is defined as an inviscid wall.   of abrupt variations without the need of any shock capturing or limiting mechanism. This property follows from the result on the monotonicity of first-order schemes in Godunov's theorem [68]. Hence, the FCFV method with the HLL numerical flux ensures the positivity properties of the approximate solution. In addition, the choice of the Riemann solver also controls the amount of numerical diffusion introduced by the FCFV method, influencing the overall accuracy of the computed solution. A qualitative comparison of the different Riemann solvers is performed in figure 11 by illustrating the pressure distribution in the fine mesh at different sections along the vertical body axis. The results display that the Lax-Friedrichs numerical flux provides non-oscillatory solutions. Nonetheless, it introduces excessive numerical dissipation leading to a smeared representation of the shock wave. The Roe Riemann solver is equipped with a Harten-Hyman entropy fix: without this correction, the method fails to converge and nonphysical solutions with localised overshoots appear. On the one hand, using an entropy fix with threshold parameter δ = 0.1, the Roe solver shows insufficient numerical dissipation producing an approximation with oscillations in the vicinity of the shock wave. On the other hand, a threshold value δ = 0.15 for the Roe solver leads to a physically-admissible and accurate solution. It is worth noticing that the parameter δ, which is problem-dependent, needs to be appropriately tuned a priori by the user. Finally, HLL-type Riemann solvers exhibit their ability to produce positivitypreserving and accurate solutions in presence of shocks without the need of any user-defined parameter, thus remedying the aforementioned issue of the Roe solver.
To further analyse the accuracy of the Riemann solvers for the FCFV method, the numerical computation of the pressure coefficient over the aerofoil surface is compared with experimental data [69]. Figure 12 confirms the overdissipative nature of the Lax-Friedrichs solution which shows a smeared representation of the shock wave. The HLL, HLLEM and the Roe solvers (the latter with an entropy fix parameter δ = 0.15) produce nearly identical solutions with a sharp representation of the shock wave showing excellent agreement with the experimental data.
Finally, a quantitative comparison is performed by computing the lift and drag coefficients, reported in table 1. According to experimental data [70], acceptable values lie within the range [0.342, 0.352] for the lift and [0.0217, 0.0227] for the drag coefficient, accounting for a tolerance of 5 lift and drag counts. The reported values for the drag coefficient em-  coefficient, the obtained results show an underestimation of this quantity, regardless of the employed Riemann solver. It is worth noticing that the proximity of the far-field boundary has a strong influence on the precision of the computed quantities, as reported in [71,72]. Indeed, the presented results are in quantitative agreement with references employing a similar domain, see e.g. [66] where the far-field boundary is located at 20 chord units from the aerofoil. In [66], the value of the lift coefficient computed using a first-order stabilised finite element approximation is 0.308, differing between 4 and 6 lift counts from the FCFV solution provided by Roe, HLL and HLLEM Riemann solvers.

Viscous laminar transonic flow over a NACA 0012 aerofoil
The next example consists of the viscous laminar transonic flow over a NACA 0012 profile at free-stream Mach number M ∞ = 0.8 and angle of attack α = 10 • . The Reynolds number, based on the chord length of the aerofoil, is Re = 500. This benchmark is presented to establish the capability of the FCFV method to concurrently capture abrupt variations due to shock waves and viscous effects in boundary layers [73,66,74].
As for the meshes utilised in the previous section, a non-uniform refinement is performed near the aerofoil surface. In addition, exploiting the information on the angle of attack of the free-stream, a priori mesh refinement is introduced in a region surrounding the aerofoil, tilted 10 • from its mean chord line, in order to accurately capture the viscous effects of the flow in the wake of the profile. An unstructured mesh of 1,005,199 triangular cells is displayed in figure 13a and a detail of its refinement on the surface of the aerofoil is reported in figure 13b. The far-field boundary is located at 15 chord units from the profile and the aerofoil surface is considered adiabatic.  The flowfield computed with the HLLEM Riemann solver is depicted in figure 14. The Mach number distribution illustrates the capacity of the method to accurately describe the detached sonic region near the leading edge as well as the appearance of a wake behind the profile. The different Riemann solvers for the FCFV method are compared for this viscous test case in figure 15. The results display the pressure and the skin friction coefficients, computed on the aerofoil surface, as well as the numerical results obtained by Kordulla in [73]. Similarly to the results observed in the inviscid simulation, the Lax-Friedrichs Riemann solver displays discrepancies with respect to the reference curves for both the pressure and the skin friction coefficient. The Lax-Friedrichs results are matched by the ones provided by the HLL Riemann solver which shows an excessive numerical dissipation  in the viscous boundary layer. The overdiffusive nature of the HLL numerical flux, not observed in the inviscid case, is attributed to its misrepresentation of contact and shear waves [53,54]. Concerning Roe numerical flux, this Riemann solver strongly depends upon the choice of the value of the entropy fix also in the viscous case. Without entropy fix (δ = 0), numerical oscillations of the solution near the leading edge appear and larger values of the threshold parameter δ are required to remedy this issue. For sufficiently large values of the entropy fix, the solution computed using the Roe Riemann solver is in good agreement with the reference one. Such an accurate approximation is also achieved by the FCFV method using the HLLEM numerical flux, without the need of tuning any parameter. Table 2 reports the values of the lift and drag coefficients, computed using different Riemann solvers. Reference data from several numerical studies based on various computational methods were collected in [73], reporting values of the lift coefficient in the range [0.415, 0.483] and of the drag coefficient in the interval [0.2430, 0.2868]. The excessive nu- Table 2: Viscous laminar transonic flow over a NACA 0012 profile -Lift, C l , and drag, C d , coefficients computed using different Riemann solvers.

Low Mach number flow over a cylinder
In this section, an incompressible flow over a 2D cylinder at angle of attack α = 0 • is considered, both in the inviscid and viscous laminar case. The objective is to show the robustness of the proposed FCFV solver for compressible flows when low Mach number flows are considered [75,66].
Unstructured meshes of triangular cells are considered for both the inviscid and the viscous simulations. On the one hand, for the inviscid case, the mesh is isotropically refined in the vicinity of the cylinder, for a total of 359, 242 cells, as displayed in figure 16a. The far-field boundary is located at 50 chord units from the cylinder where inviscid wall boundary conditions are imposed. On the other hand, the refinement of the boundary layer and of the wake of the cylinder in the viscous simulation leads to a mesh of 654, 194 triangular cells, reported in figure 16b. In this case, the far-field boundary is placed at 20 chord lengths from the obstacle and the surface of the cylinder is considered adiabatic.  The isolines of the Mach number distribution computed with the HLLEM Riemann solver for the inviscid flow are reported in figure 17 for different values of the far-field condition. The results display the robustness of the FCFV solver for low Mach number simulations, highlighting the capability of the method to devise non-oscillatory solutions even in the incompressible limit. More precisely, the computed solution is not deteriorated by the decrease of the Mach number, even when it approaches zero. Note that the loss of symmetry of the solution is due to the geometric error introduced by the piecewise linear approximation of the surface of the cylinder. This well-known problem, see [76], is related to the production of nonphysical entropy by the low-order discretisation of the curved boundary. It is worth recalling that the objective of this test is to show the robustness of the proposed method in the incompressible limit. In order to remedy the above mentioned issue, several approaches proposed in the literature can be employed within the FCFV paradigm. These include high-order approximation of the geometry [76], appropriate modification of the wall boundary condition [77] or exact treatment of the geometry via the NURBS-enhanced finite element method [78]. The robustness of the FCFV solver in the incompressible limit of the compressible Navier-Stokes equations is studied through a steady-state flow at Re = 30. Figure 18 displays the isolines of the Mach number distribution for far-field conditions at M ∞ = 0.1 and M ∞ = 0.01. In this case, the FCFV solver is able to precisely approximate the flow in the wake of the cylinder, with no loss of accuracy when approaching the incompressible limit. It is worth noticing that in both figure 17 and figure 18, the variation of the farfield boundary condition only affects the scale of the computed Mach number and not its distribution.
Finally, the values of the pressure and skin friction coefficients computed using the HLLEM Riemann solver are compared with the results reported in [66] of a stabilised finite element simulation performed with polynomial approximation of degree k = 3. Good agreement is displayed in figure 19 for both the pressure and the skin friction coefficients, confirming the capability the proposed methodology of accurately simulating viscous laminar flows also in the incompressible limit. In particular, it is worth noticing that the FCFV results computed using the two values of the far-field condition are almost identical and they are not affected by the value of the Mach number going to zero. Hence, the FCFV  This benchmark constitutes a classic CFD validation example for external flows due to its complex flow physics and the availability of experimental results at high Reynolds number [79].
A non-uniform mesh refinement is adopted in the vicinity of the wing surface and towards the leading and trailing edges. Figure 20 details two levels of refinement with meshes consisting of 236,682 and 5,061,252 tetrahedral cells, respectively. The far-field boundary is located at approximately 12 chord lengths from the wing and the aerofoil surface is defined as an inviscid wall.  The FCFV simulation is performed using the HLL Riemann solver, based on its capability of producing positivity-conserving solutions and on its robustness in inviscid supersonic cases. The linear system of equations arising from the FCFV discretisation is solved using a GMRES with restarting parameter 10 and no preconditioner. A detail of the flow computation obtained on the fine mesh is reported in figure 21. The Mach number and pressure distributions clearly show that the FCFV method is able of accurately capturing the characteristic lambda-shock arising in this test benchmark.
The performance of the FCFV method is evaluated by examining the pressure coefficient, computed on the fine mesh, at different sections along the wing span, see figure 22. The obtained results are compared both to experimental data [79] and to computational simulations performed using second-order CCFV and VCFV solvers on the same mesh. More precisely, the commercial CFD software Ansys Fluent [3] is employed as CCFV solver, whereas the VCFV results are obtained using the CFD solver FLITE [7,8]. An upwind scheme is utilised for the treatment of the numerical fluxes in Fluent, whereas a Roe Riemann solver is selected for the FLITE simulation. The numerical results obtained with the FCFV method show excellent agreement both with the experimental data and with the remaining FV schemes. It is worth noticing that both second-order FV schemes produce a small oscillation in the representation of the shock-wave located at the upper mid-chord. On the contrary, the first-order FCFV method based on the positivity-preserving HLL  Riemann solver is capable of computing non-oscillatory solutions, establishing a robust framework for the simulation of 3D problems involving complex flow features.

Concluding remarks
The face-centred finite volume paradigm was proposed for the first time for the approximation of nonlinear hyperbolic PDEs modelling compressible flows. The method is based on a mixed formulation and defines the unknowns, that is, the hybrid vector of conservative variables, at the barycentre of the faces. The unknowns in each cell, i.e. density, momentum, energy, deviatoric strain rate tensor and gradient of temperature, are eliminated via a hybridisation procedure to reduce the global number of degrees of freedom of the problem. In addition, traditional Riemann solvers, i.e. Lax-Friedrichs, Roe, HLL and HLLEM, are devised in the context of FCFV discretisations via appropriate definitions of the numerical fluxes.
The presented methodology provides first-order accuracy of the conservative quantities, as well as of the stress tensor and of the heat flux. More precisely, the first-order convergence of the stress tensor and the heat flux is achieved without the need to perform a reconstruction of the gradients as required by traditional second-order CCFV and VCFV strategies. The FCFV paradigm is thus robust on unstructured meshes and retains optimal accuracy even with highly stretched or distorted cells, avoiding well-known issues of traditional FV schemes.
In addition, the FCFV method is able to construct non-oscillatory solutions of sharp fronts without the need of any shock capturing or limiting technique. The accurate treatment of shocks, expansion fans and shear waves is naturally handled by the Riemann solvers implicitly embedded in the FCFV numerical fluxes.
Finally, the method is robust in the incompressible limit, allowing to seamlessly simulate flows at low Mach number, without the need of introducing specific pressure corrections like the well-known SIMPLE algorithm.
A comprehensive set of two-and three-dimensional numerical examples is employed to demonstrate the optimal convergence properties of the method and its capabilities to solve complex flow problems of aerodynamic interest across various regimes, from inviscid to viscous laminar flows, from transonic to subsonic incompressible flows. Moreover, a detailed comparison of the accuracy and robustness of different Riemann solvers is presented. The FCFV method equipped with HLL-type Riemann solvers thus provides a solution strategy suitable for all flow regimes, which outperforms the traditional Lax-Friedrichs numerical flux in terms of accuracy and the Roe solver in terms of robustness.
An extension of the proposed first-order FCFV approach to achieve second-order accuracy in the primal variable will be investigated in future works starting from the recent results obtained for the second-order FCFV method in the context of second-order elliptic PDEs [22,21]. This approach relies on two ingredients. On the one hand, a piecewise linear approximation is employed for the solution in the cells, while its gradient in the cell and the hybrid solution on the faces are maintained piecewise constant. On the other hand, a projection operator is introduced to propose a new definition of the numerical fluxes. In this context, the number of globally-coupled degrees of freedom of the second-order FCFV method is the same as the original FCFV method. On the contrary, to achieve second-order convergence of the primal variable, a standard HDG method [43,47,48] employs polynomial functions of degree one for all the variables, including the hybrid one, substantially increasing the number of unknowns and the size of the resulting system of equations.

A Boundary conditions for the FCFV formulation
In this appendix, the definitions of the boundary conditions for compressible flow problems are briefly recalled in the context of hybrid discretisation methods [48,42,80,43,47]. Following the HDG framework, the implementation of boundary conditions in the FCFV formulation of compressible flows rely on the exploitation of the hybrid vector U .
Let the boundary ∂Ω be partitioned as ∂Ω = Γ ∞ ∪ Γ out ∪ Γ ad ∪ Γ iso ∪ Γ inv , where the introduced portions are disjoint by pairs. From the modelling viewpoint, Γ ∞ refers to a far-field boundary, Γ out is a subsonic outlet with imposed pressure, Γ ad and Γ iso denote adiabatic and isothermal walls, respectively, whereas Γ inv represents a symmetry boundary or an inviscid wall with slip conditions. The expressions of the trace boundary operator B(U , U , ε d , φ) corresponding to the above mentioned cases are reported in table 3. As classic in the context of compressible flows, inlet and outlet boundaries are identified through a 1D characteristics analysis in the direction of the outward normal to the boundary. More precisely, A ± n := (A n ±|A n |)/2 denote the positive and negative parts of the matrix A n ( U ) and they are defined exploiting the spectral decomposition introduced in section 3.4. Moreover, in table 3, U ∞ denotes the free-stream values of the conservative variables at the far-field, whereas p out and T w are the prescribed values of the outlet pressure and the wall temperature, respectively. Finally, the stabilisation coefficient τ d ρE = 1/ [Re(γ − 1)M 2 ∞ Pr] is obtained extracting the component associated with the energy equation from the diffusive stabilisation tensor (15).

Γ ∞
Far-field, subsonic inlet, supersonic inlet/outlet Subsonic outlet (pressure outlet) Symmetry surface or inviscid wall B = {ρ − ρ, [(I nsd − n ⊗ n)ρv − ρv] T , ρE − ρE} T the matrix N V , accounting for the normal to a surface, being defined as where n k denotes the k-th component of the unit normal vector n.

C Hybridisation of the FCFV solver
This appendix describes the implementation details for the hybridisation procedure and the final system obtained for the hybrid vector U.
It is worth recalling that the FCFV method features two stages, see section 3.3. First, the n el local problems (20) are devised, in order to eliminate the unknowns (U e , ε d e , φ e ) within each cell by expressing them as functions of the hybrid vector U . Stemming from equation (23), the FCFV local problem for Ω e , e = 1, . . . , n el is given by A e U U U e + A e U Q Q e = A e U U U + F e U , where the matrices and vectors above arise from the Newton-Raphson linearisation of nonlinear system of equations (20).
Following from the constant degree of approximation utilised to approximate U e and Q e at the centroid of each cell and U at the barycentre of each face and from the quadrature rule employing a single integration point on cell and faces, the primal, U e , and mixed, Q e , variables are expressed as functions of the hybrid unknown U in a decoupled manner, namely It is worth noticing that the computations in equation (37) are independent cell-by-cell and only involve the inverses of matrices A e U U and A e QQ . The former is a matrix of dimension (n sd + 2) × (n sd + 2), that is, 4 × 4 in 2D and 5 × 5 in 3D. The latter is the identity matrix of dimension (m sd + n sd ) × (m sd + n sd ) (i.e., 5 × 5 in 2D and 9 × 9 in 3D) scaled by the volume of the cell Ω e . Hence, this step requires a reduced computational effort and can be easily performed in parallel.
Similarly, upon linearisation via the Newton-Raphson method, the global problem (22) is expressed as By plugging the expressions obtained from equation (37) into equation (38), the number of unknowns is reduced by eliminating the local unknowns U e and Q e from the global problem. Hence, at each Newton-Raphson iteration, the linear system is solved, where the matrix K and the vector F are obtained from the assembly of the contributions from each cell, namely