ARTIFICIAL COMPRESSIBILITY WITH RIEMANN SOLVERS: CONVERGENCE OF LIMITERS ON UNSTRUCTURED MESHES

Free-surface flows and other variable density incompressible flows have numerous important applications in engineering. One way such flows can be modelled is to extend established numerical methods for compressible flows to incompressible flows using the method of artificial compressibility. Artificial compressibility introduces a pseudo-time derivative for pressure and, in each real-time step, the solution advances in pseudo-time until convergence to an incompressible limit—a fundamentally different approach than SIMPLE, PISO, and PIMPLE, the standard methods used in OpenFOAM® . Although the artificial compressibility method is widespread in the literature, its application to freesurface flows is not. In this paper, we apply the method to variable density flows on 3D unstructured meshes, implementing a Godunov-type scheme with MUSCL reconstruction and Riemann solvers, where the free surface gets captured automatically by the contact wave in the Riemann solver. The critical problem in this implementation lies in the slope limiters used in the MUSCL reconstruction step. It is well-known that slope limiters can inhibit convergence to steady state on unstructured meshes; the problem is exacerbated here as convergence in pseudo-time is required not just once, but at every realtime step. We compare the limited gradient schemes included in OpenFOAM® with an improved limiter from the literature, testing the solver against dam-break and hydrostatic pressure benchmarks. This work opens OpenFOAM® up to the method of artificial compressibility, breaking the mould of PIMPLE and harnessing high-resolution shock-capturing schemes that can scale better in parallel.


Introduction
Variable density incompressible flows are of much practical interest in engineering. One important class of variable density incompressible flows is free-surface flows of water and air, for example, dam breaks [1], hydraulic jumps [2], hydraulic structures [3], and nature-based flood defences, both coastal [4] and fluvial [5]. These flows in these papers [1,2,3,4,5] were all modelled using the volume-of-fluid (VOF) solver interFoam. As with most OpenFOAM ® solvers, interFoam uses the PIMPLE algorithm to update the pressure field. In this paper, we harness the flexibility of OpenFOAM ® to implement an alternative free-surface solver based on the method of artificial compressibility [6], which has a key computational advantage compared to PIMPLE.
Both artificial compressibility and PIMPLE address the problem of the incompressible Navier-Stokes equations having no equation for pressure. In artificial compressibility, the idea is to add a pseudo-time derivative for pressure to the incompressibility constraint, and the solution advances in this pseudo-time until convergence to an incompressible limit [7,8,9]. This makes the governing equations hyperbolic, and thus solvable with the wealth of methods developed for compressible flows, for example, Godunov-type schemes involving MUSCL (monotonic upstream-centered scheme for conservation laws) reconstruction and Riemann solvers. While the method was originally developed for steady-state cases, it can be generalised to transient cases by employing dual time stepping, where the solution converges to the incompressible limit each real-time step. This is a fundamentally different approach to the PIMPLE algorithm. PIMPLE, a combination of SIM-PLE (semi-implicit method for pressure-linked equations) [10] and PISO (pressure-implicit with splitting of operators) [11], is a predictor-corrector method that relies on matrix inversion. While such methods require fewer total iterations than artificial compressibility and therefore would be more efficient in serial [12], artificial compressibility is easier and more efficient to parallelise [12,13,14] and can lead to almost linear speed-up [15]. This is because solving local wave propagation problems requires less communication between processors than matrix inversion [12]. Therefore, it is worthwhile to depart from the convention of using PIMPLE in OpenFOAM ® and investigate implementing artificial compressibility instead because of its scaling potential on the massively parallel architectures of today.
Capitalising on the powerful tools provided by OpenFOAM ® , we generalise the specific variabledensity artificial-compressibility scheme in [22] to 3D unstructured meshes, testing the solver against dam-break and hydrostatic pressure benchmarks. Even though the new solver does not use much of the numerics in OpenFOAM ® , it is useful to develop the solver in OpenFOAM ® because the lowlevel structure does not have to be built from scratch. A downside is that the form of parallelisation in OpenFOAM ® , domain decomposition, is not optimal for artificial compressibility as the waves propagate locally [12]. However, we focus on the numerics in this paper, leaving the implementation of efficient parallelisation to future work.

Implementation
The variable-density artificial-compressibility equations are where ρ is density, u velocity, p pressure, µ dynamic viscosity, g gravity, β the artificial compressibility coefficient, t real-time, and τ pseudo-time. Each real-time step requires convergence in pseudo-time, and so two nested time loops are needed for the dual time stepping, as well as a Runge-Kutta loop. If the solution converges in pseudo-time, then ∂p/∂τ ≈ 0 and the artificial compressibility coefficient β should not impact the solution. However, when the maximum density is around ρ = 10 3 as for water, setting β < 10 3 can result in instability and setting β > 10 4 can inhibit convergence, meaning ∂p/∂τ ≈ 0 is not achieved [17]. Therefore some trial and error is required to pick a suitable value of β. In practice, we have found that setting β = 1100 is just large enough to avoid instability, and so will have minimal impact on convergence. We solve the governing equations (1)-(3) by implementing in OpenFOAM ® the numerical scheme described in detail and implemented on 2D Cartesian meshes by [22]. In that study, a Godunov-type scheme was developed that was explicit in pseudo-time and point-implicit in real-time, especially easy to parallelise, and had a low memory footprint. Specifically, the explicit Runge-Kutta pseudo-time update equation is given by where Q is the state vector, Ω the computational cell, F the flux vector, I C the cancellation matrix, B the body force term, α S the Runge-Kutta coefficient, s max the total number of Runge-Kutta stages, α P I the point-implicit scaling coefficient, n the current real-time iteration, m the current pseudo-time iteration, and s the current Runge-Kutta stage. The real-time derivative is treated as a source term, calculated point-implicitly with a backward-differencing scheme, and the fluxes are approximated using a Riemann solver. For full details, see [22]. Note that we consider the inviscid flux only (i.e. we set µ = 0), but the viscous flux could be added easily using standard finite difference methods [17,43]. The present paper goes further than [22], implementing the numerical scheme (4)- (6) in OpenFOAM ® so that it can be used on 3D unstructured meshes and in parallel. The new implementation uses the Roe Riemann solver and the new pressure gradient calculation developed in [22]. It has a structure inspired by foam-extend's dbnsFoam [61]. However, it was written from scratch to take into account the different governing equations and the fact that several cell-limited gradient schemes were added to OpenFOAM ® since dbnsFoam was released. Therefore, it has the potential to be forward-compatible with any limiters added to OpenFOAM ® in the future.
Switching to an unstructured mesh does require slightly changing the scheme from [22]. The difference is in the MUSCL reconstruction step, which is where the primitive variables are extrapolated to cell faces to provide the left and right states for the Riemann solver. This extrapolation relies upon a limited gradient to avoid the introduction of new extrema, Gibbs-type oscillations. Unstructured meshes require a different approach than structured meshes to do MUSCL reconstruction, and there are two components to this difference: calculating the gradient and limiting the gradient.
First, the gradient calculation for Cartesian meshes is very simple, but unstructured meshes require a method such as Green-Gauss or least-squares. These well-known standard methods are already implemented in OpenFOAM ® , so using them is simply a case of writing fvc::grad() in the code and specifying Gauss or leastSquares in fvSchemes. The gradients can be limited by specifying a limiter alongside Gauss or leastSquares in fvSchemes. Note that this use of fvc::grad() is the only time we use OpenFOAM ® numerics here, but it is not the gradient calculation itself which is problematic.
The problem lies in the limiters, a problem that did not exist in the implementation for 2D Cartesian meshes in [22]. The difference is that, while the limiter calculation for structured meshes is straightforward, involving constructing left-and right-sided gradients at each cell and applying a standard limiter to these, it is not obvious what the left-and right-sided gradients would be on an unstructured mesh [62]. The standard way to overcome this problem is the framework of Barth and Jespersen [63], which is the method implemented in the cell-limited gradient schemes in OpenFOAM ® .

Barth and Jespersen.
Let W be a scalar variable or one of the components of a vector variable. In the context of a Godunov-type scheme, the idea is to construct a scalar limiter function Φ ∈ [0, 1] to be used as where ∇W i is gradient of the variable before limiting and W f,i is the Riemann state at face f on the side belonging to cell i. Then the Riemann states either side of the face can be put into the Riemann solver to get the numerical flux. Barth and Jespersen [63] introduced a method that ensures the interpolated variable W f,i does not exceed that of the neighbouring cells in magnitude. Note that this means all the neighbouring cells of cell i, not just the cell on the other side of face f . Practically, the calculation loops over all the cell interfaces twice [64]. There is a loop to calculate the maximum and minimum values in the neighbouring cells, and then a loop to determine how much limiting is required, where ϕ(y) = min(1, y) (11) and Thus, the quantity ∆ i,f is a measure of how much W changes between the centre of cell i and the centre of the cell interface f . The function ϕ clips ∆ i,f so that it does not exceed ∆ i,min or ∆ i,max in magnitude.
In the uniform 1D case, the Barth-Jespersen limiter can be recast into the Spekreijse [65] form where the limiter is and it is applied to the ratio as shown in [66]. Both ϕ from (11) and ψ from (17) are plotted in Figure 1.
The practical problem with this method is that the non-differentiability of (11) can inhibit convergence in steady-state cases, instead leading to bounded odd-even modes. The simulation is not unstable; it does not crash or blow up, but neither does it make any progress towards convergence. It would continue indefinitely if it was not stopped. This is a well-known problem [48,64,66,67,68,69,70,71,72,73]. Fortunately, there are ways to overcome it, including replacing the non-differentiable function (11) with a differentiable alternative and switching off the limiter in areas of relatively uniform flow. In the present paper, we investigate these techniques in the context of artificial compressibility, where pseudo-time convergence is required not just once, but at every real-time step, and thus the limiters must converge more reliably than for a simple steady-state simulation.
2.2. Venkatakrishnan. Venkatakrishnan [66] addressed the convergence problem by replacing the nondifferentiable function (11) with the differentiable function ϕ(y) = y 2 + 2y As shown in Figure 1, this function is much smoother than (11), and the fact that ϕ > 1 for y > 2 does not affect ψ, which is bounded by the Barth-Jespersen limiter for all r.
A slight modification, not plotted in Figure 1, ensures that the limiter is also not activated in smooth regions: where and V i is the cell volume [64]. The parameter K is a threshold that marks the largest size of oscillations untouched by the limiter. If the gradient is small, then which means there is no limiting. This stops residuals stalling due to numerical noise. By a similar argument, increasing K increases Φ i and so decreases the amount of limiting, and while this is more conducive to convergence to a steady-state, it also increases the potential for Gibbs-type oscillations and thus instability. The more convergent options are more unstable, that is, there is a trade-off between convergence and stability.

2.3.
Michalak and Ollivier-Gooch. While Barth-Jespersen and Venkatakrishnan are the two standard limiters [64], there are more options in the literature. For example, Michalak and Ollivier-Gooch [67] replaced the function (11) with where P (y) is the cubic polynomial with P 0 = 0 (23) and 1 ≤ y t ≤ 2 is a threshold. The polynomial itself is not explicitly stated in [67], but a simple derivation gives As shown in Figure 1, the function ϕ is differentiable everywhere and does not exceed 1. Like Venkatakrishnan, Michalak and Ollivier-Gooch also proposed switching off the limiter in uniform regions of flow. However, they used a different method for this, smoothly transitioning to switching off the limiter when This is done by definingφ (22) andφ i now plugged into (10). The indicator-type function is given by with the smooth transition function s given by s(y) = 2y 3 − 3y 2 + 1.
2.4. Cell-limited schemes in OpenFOAM ® . The above three limiters are already implemented in OpenFOAM ® , and can be easily accessed by specifying cellLimited, cellLimited<Venkatakrishnan>, or cellLimited<cubic> respectively alongside the gradient calculation in fvSchemes. However, there are some major problems with the limiters as implemented in the OpenFOAM ® source code. First, cellLimited<Venkatakrishnan> is not the same as the original limiter put forward by Venkatakrishnan, and so it cannot be expected to have the same convergence properties shown in the original study [66]. One difference is that it uses the unmodified version (19) instead of (20), and so may still be active in uniform regions of flow. Another difference is that, regardless of the limiter, OpenFOAM ® always clips Φ so that it never exceeds 1. As indeed noted in the source code documentation, this clipping makes this particular limiter non-differentiable, and so it "no longer conforms to the basic principles of this kind of limiter function" (VenkatakrishnanGradientLimiter.H, lines 53-54). The whole reason Venkatakrishnan developed the limiter was so it could be differentiable, therefore it seems that cellLimited<Venkatakrishnan> is of little practical use.
Second, cellLimited<cubic> is also not the same as the original limiter put forward by Michalak and Ollivier-Gooch, and so cannot be expected to have the same convergence properties shown in the original study [67] either. One difference is that it does not use (31), and therefore does not stop activation in uniform regions of flow. Another difference is that the limiter uses an incorrect cubic polynomial, one with a very large slope discontinuity at ϕ(y t ), not the one stated in (27)- (29). * In this study, we do not attempt to fix cellLimited<Venkatakrishnan>; the automatic clipping of Φ makes this too problematic. However, we do implement an improved version of cellLimited<cubic> that corresponds to the limiter in the original paper [67], both by using the correct cubic polynomial and by stopping activation in uniform regions of flow. This improved version is called cellLimited<Michalak>.
Note that there are three other limited gradient schemes in OpenFOAM ® : cellMDLimited, faceLimited, and faceMDLimited. In the cellLimited limiters outlined above, one scalar limiter is applied to the x, y, and z components of the gradient equally, irrespective of which neighbouring cells have the minimum and maximum values. This can lead to excessive limiting. The cellMDLimited limiter takes an alternative approach. Consider a cell with faces f = 1, 2, ..., F . Set (∇W ) i,0 = (∇W ) i to be the gradient before limiting. Loop through the faces f , calculating the extrapolate and then clipping the gradient in the direction from the cell centre to that face centre: for f = 1, 2, ..., F . The final iteration gives us the cellMDLimited gradient. Therefore, the limiter "is applied to the gradient in each face direction separately" (cellMDLimitedGrad.H, lines 36-38) not, as is sometimes suggested, to each coordinate direction separately. Meanwhile, the faceLimited and faceMDLimited limiters are like the cellLimited and cellMDLimited limiters but clip the gradient between the face-neighbour values rather than the cell-neighbour values. None of these allow differentiable functions as a run-time selectable option like cellLimited does. However, simulations using them are included for comparison purposes in the following section, alongside cellLimited, cellLimited<Venkatakrishnan>, cellLimited<cubic>, and the newly implemented cellLimited<Michalak>.

Application
The new solver was put through three benchmark tests: a simple dam break to illustrate the importance of limiter choice, a more complicated dam break to compare the new solver against interFoam, and a hydrostatic case to demonstrate that the solver can be run on arbitrary 3D unstructured meshes.
3.1. Convergence of limiters. First, we replicated from [22] the first real-time step of a dam break on a uniform 2D Cartesian mesh. The mesh was a metre length in each direction, divided into 20 × 20 cells. Despite its simplicity, the mesh had to be stored as a 3D unstructured mesh in OpenFOAM ® . However, due to its simplicity, the mesh was useful for isolating the effect of different limiters on convergence. Initial conditions were set to be ρ = 1000 in the water and ρ = 1 in the air, with u = 0 and p = 0 everywhere. The solver was run for one real-time step of size ∆t = 0.01 seconds with the following settings: three Runge-Kutta stages, an artificial compressibility coefficient of β = 1100, and a Courant number for pseudo-time of 0.48 to satisfy the explicit stability constraint. The Green-Gauss gradient calculation was chosen for ρ and u, and the corresponding limiter changed in each simulation to show its effect on convergence.
Recall that the solution at each real-time step is reached when the solution converges in pseudotime. This is measured with residuals, here calculated as the maximum absolute difference between a conserved variable in the current and previous pseudo-time steps. Figure 2 shows that, of all the limiters available in OpenFOAM ® as standard, only faceLimited converged in this case, and not very smoothly. It might seem like the convergence problem was because all the OpenFOAM ® limiters are non-differentiable, as discussed in Section 2.4. However, Figure 2 shows that, despite being differentiable, cellLimited<Michalak> (K = 0) did not converge either. It was only when the limiter was switched off in areas of uniform flow that convergence was achieved (K = 1). Clearly, when K = 1, the newly implemented limiter cellLimited<Michalak> has improved convergence properties compared to the options already present in OpenFOAM ® .
It is important to note that, although the residuals stalled for nearly all of the other limiters, there was no instability and therefore the simulations did not blow up. Only a few cells failed to converge, and they switched between very similar states, which meant the final results did not differ substantially between the limiters. Consider Figure 3. The residuals only stalled in the top-left corner above the water column. Once K was changed to 1 for cellLimited<Michalak>, the limiter switched off in this area of uniform flow.

3.2.
Comparison with interFoam. The new solver was tested with the damBreak tutorial to compare its performance with interFoam. Initial conditions were set to be ρ = 1000 in the water and ρ = 1 in the air, with u = 0 and p = 0 everywhere. The simulations were run until t = 1.0 or until the residuals stalled, whichever was sooner, with the following settings: one Runge-Kutta stage, an artificial compressibility coefficient of β = 1100, and a Courant number for pseudo-time of 0.48 to satisfy the explicit stability constraint. At each real-time step, the absolute tolerances for pseudo-time convergence were 0.0001 for ρ, 0.01 for ρu, and 0.01 for p. The Green-Gauss gradient calculation was chosen for ρ and u, and the corresponding limiter changed along with the Courant number for real-time. Recall that real-time is treated point-implicitly, and so there is no explicit stability constraint for the real-time Courant number. Figure 4 shows how far each simulation could go. The residuals were prone to stalling when the tail of the column of water hit either the obstacle or the far wall and, as expected, cellLimited<Michalak> converged better than cellLimited. However, sometimes cellLimited<Michalak> stalled, and not always with the same parameters. In this particular benchmark, when the real-time Courant number was 1, the choice y t = 1.5 converged but y t = 2.0 stalled, and it was the other way around if the Courant number was 2 or 5. Both choices of y t stalled if the Courant number was 0.5, perhaps due to the increased number of real-time steps meaning there were more chances to stall. Moreover, the effect of switching off the limiter in uniform regions of flow (that is, setting K > 0) was not noticeable. This suggests that the standard slope limiters for unstructured meshes are not sufficiently robust for this solver, whether the limiters already available in OpenFOAM ® or the limiter implemented in this paper. Again, this is nothing to do with stability. At no point did the simulations blow up or crash; they simply stopped making progress towards convergence.
Despite this lack of convergence, we could use parameters that worked for this particular benchmark to compare the new solver, as it stands, with interFoam. Figures 5-6 show the results for cellLimited<Michalak> with y t = 2.0, K = 1, and real-time Courant numbers of 2 and 5. Figure  5 includes results for the interFoam simulations, where the settings were as in the standard laminar damBreak tutorial, but with zero viscosity and surface tension, and the atmosphere patch a wall. Now, interFoam keeps the free surface sharp using a compression term with coefficient cAlpha. Although the compression term has a physical basis [74, p.117], in practice cAlpha is usually arbitrarily set to 1. Setting cAlpha to 0 is equivalent to removing the interface compression term-not a practice acceptable for practical simulations but useful nonetheless for this comparison. Figure 5 shows that the new solver with a real-time Courant number of 2 had very similar behaviour to interFoam with a cAlpha of 0. Surprisingly, when the real-time Courant number was increased to 5, the new solver still managed to capture the highly transient behaviour of the dam break, with only slightly less detail. This is encouraging Often called artificial compression, completely unrelated to artificial compressibility. because fewer total iterations are required for this larger real-time Courant number, as shown by Figure  6.
While the solver has potential, Figure 5 shows that the use of a high-resolution Godunov-type scheme is not sufficient by itself to keep the interface sharp. Activating the interface compression term does indeed keep the interface sharp when using interFoam, but this comes at the cost of a deformed free surface (see the t = 0.6 frame). It was hoped that a Godunov-type scheme might circumvent this problem, but clearly that is not the case, and something like the interface compression term in interFoam is required to keep the interface sharp. This has been done before in the context of Godunov-type schemes for artificial compressibility coupled with VOF [57,55], and would be a good next step once the convergence problem has been solved.
3.3. Mesh from snappyHexMesh. The meshes in the previous benchmarks appeared structured, even if they were not stored that way in OpenFOAM ® . Therefore, to demonstrate that the solver can be run on arbitrary meshes, it was also tested on a much less uniform mesh generated by snappyHexMesh. A simple hydrostatic case was chosen to avoid the convergence problem and because the solution is known, allowing us to test the impact of the pressure gradient calculation on the simulation. The mesh from the iglooWithFridges tutorial was used, with all the boundary conditions changed to walls, and the domain filled with water up to z = 2, as shown in Figure 7. As before, the initial conditions for velocity and pressure were set to zero. The solver was run for one real-time step of size ∆t = 0.01 seconds with the following settings: one Runge-Kutta stage, an artificial compressibility coefficient of β = 1100, and a Courant number for pseudo-time of 0.48 to satisfy the explicit stability constraint. The absolute tolerances for pseudo-time convergence were 0.01 for ρ, 0.01 for ρu, and 0.01 for p. The Green-Gauss gradient calculation was chosen for ρ and u, and the corresponding limiter set to cellLimited<Michalak> with y t = 1.5 and K = 1.
Recall that the pressure gradient calculation used throughout this paper is from [22]. It is based on a rearrangement of the momentum equation, and so it is automatically well-balanced, that is, is balances exactly with the gravity source term when the velocities are zero. While there are other well-balanced methods [19,44,75,76,77], this one is very easy to implement on unstructured meshes. Therefore, we investigate its effectiveness on unstructured meshes by comparing it to when the pressure gradient is calculated instead by cellLimited<Michalak> with y t = 1.5 and K = 1. Figure 8 shows that the method from [22] indeed performs better along the direction of gravity than the limiter cellLimited<Michalak>, which is too limiting. The resulting difference in the pressure and velocity in the direction of gravity is notable. The small but non-zero velocities in the other directions can be attributed to the free surface not starting out entirely flat due to the irregular shape of the mesh. Figure 9 shows that the improved performance of the well-balanced limiter comes at the cost of slower convergence, but this is primarily due to the fact that it takes longer for the pressure field to converge to the correct magnitude from p = 0. In any case, this test demonstrated that the new solver can be run on 3D unstructured meshes and that the pressure gradient calculation of [22] retains its well-balancedness, although not as cleanly as on a 2D Cartesian mesh.

Conclusion
In OpenFOAM ® , the traditional way to model free-surface flows is with the VOF solver interFoam, which uses the PIMPLE algorithm based on matrix inversion. However, there are other ways to deal with the lack of a pressure equation in the incompressible Navier-Stokes equations. In this paper, we focus on the method of artificial compressibility, which scales better than matrix inversion as it requires less communication between processors. Indeed, although artificial compressibility is not new, it is currently very relevant due to its suitability for implementation on massively parallel architectures [12,13,14].
Despite its potential for efficient parallelisation, artificial compressibility has not been investigated extensively for variable density flows (including free-surface flows), having only been implemented on 2D and structured meshes so far [17,19,20,22]. The present study harnesses OpenFOAM ® to implement artificial compressibility for variable density incompressible flows on 3D unstructured meshes. We found that, since pseudo-time convergence is required at every real-time step, the slope limiter used in the MUSCL reconstruction step is critical, and the limiters currently in OpenFOAM ® are not fit for this purpose. Therefore, we implemented Michalak and Ollivier-Gooch's [67] limiter fully to avail of its  improved convergence properties. When the results converge, they both compare well with interFoam and are well-balanced, but even the improved limiter is not sufficiently robust to ensure convergence. While the solver requires a mechanism to keep the free surface sharp as in interFoam, the convergence problem is more urgent. Without convergence, there is no solution.
Consequently, further research should first focus on alternative routes to second-order accuracy than Barth and Jespersen's [63] framework for MUSCL. This could involve smoothing out the multidimensional limiter cellMDLimited rather than cellLimited, or implementing the WENO (weighted essentially nonoscillatory) [78] method. Of course, convergence is necessary but not sufficient for a good solution, and so more work would need to be done after achieving robust convergence: a mechanism could be employed to keep the interface sharp, the solver could be compared thoroughly with interFoam, and it could be applied to practical cases. Ultimately, this will allow us to use Riemann solvers to capture the free surface automatically and in a way that is particularly suited to the massively parallel architectures of today.