Quantifying the checkerboard problem to reduce numerical dissipation

Journal of Computational Physics 521 (2025) 113537 Available online 29 October 2024 0021-9991/© 2025 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp Quantifying the checkerboard problem to reduce numerical dissipation J.A. Hopmana,∗, D. Santosa, À. Alsalti-Baldelloua,b, J. Rigolaa, F.X. Triasa aHeat and Mass Transfer Technological Center, Technical University of Catalonia, ESEIAAT, c/Colom 11, Terrassa, 08222, Barcelona, Spain bTermo Fluids SL, c/Magí Colet 8, Sabadell, 08204, Barcelona, Spain1 A R T I C L E I N F O A B S T R A C T Keywords: Checkerboarding Collocated grids Conservative discretisation This work provides a comprehensive exploration of various methods in solving incompressible flows using a projection method, and their relation to the occurrence and management of checkerboard oscillations. It employs an algebraic symmetry-preserving framework, clarifying the derivation and implementation of discrete operators while also addressing the associated numerical errors. The lack of a proper definition for the checkerboard problem is addressed by proposing a physics-based coefficient. This coefficient, rooted in the disparity between the compactand wide-stencil Laplacian operators, is able to quantify oscillatory solution fields with a physics-based, global, normalised, non-dimensional value. The influence of mesh and time-step refinement on the occurrence of checkerboarding is highlighted. Therefore, single measurements using this coefficient should be considered with caution, as the value presents little use without any context and can either suggest mesh refinement or use of a different solver. In addition, an example is given on how to employ this coefficient, by establishing a negative feedback between the level of checkerboarding and the inclusion of a pressure predictor, to dynamically balance the checkerboarding and numerical dissipation. This method is tested for laminar and turbulent flows, demonstrating its capabilities in obtaining this dynamical balance, without requiring user input. The method is able to achieve low numerical dissipation in absence of oscillations or diminish oscillation on skew meshes, while it shows minimal loss in accuracy for a turbulent test case. Despite its advantages, the method exhibits a slight decrease in the second-order relation between time-step size and pressure error, suggesting that other feedback mechanisms could be of interest. 1. Introduction The checkerboard problem arises due to the non-trivial pressure-velocity coupling for incompressible Newtonian fluid flows. The governing equations for such flows are given by the momentum and continuity equations: 𝜕𝑡𝐮+(𝐮⋅∇)𝐮=𝜈∇2𝐮−1 𝜌∇𝑝, (1) ∇⋅𝐮=0,(2) * Corresponding author. E-mail address: [email protected] (J.A. Hopman). URL: https://github.com/janneshopman (J.A. Hopman). 1www .termofluids .com. https://doi.org/10.1016/j.jcp.2024.113537 Received 13 June 2024; Received in revised form 24 September 2024; Accepted 24 October 2024 Journal of Computational Physics 521 (2025) 113537 2 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. Fig. 1. Schematic drawing of how the Central Difference Scheme leads to decoupling of nodes. which, for three spatial dimensions, provide only three independent equations and four unknowns. For incompressible flows, the pressure cannot be obtained from the equation of state, and formulating an equation for pressure becomes non-trivial, a topic that has been widely studied in the field of computational fluid dynamics (CFD) [1–4]. Many methods solve this problem iteratively by predicting a velocity field and then finding a pressure field of which the gradient projects it onto a divergence-free space, using the Helmholtz-Hodge theorem [5–9]. If the discrete gradient at node 𝑖is derived through central differencing, its value will only depend on the pressure at neighbouring nodes 𝑗and not on the pressure at node 𝑖itself. Moreover, if done consistently, the discrete Laplacian operator derived from this gradient will be based on a wide stencil. In this stencil, node 𝑖is coupled to nodes 𝑘, which are neighbours to nodes 𝑗, resulting in a decoupling between node 𝑖and directly neighbouring nodes 𝑗, see Fig. 1. The use of this so-called widestencil Laplacian in the Poisson equation and the aforementioned gradient in the velocity correction results in a decoupling of the pressure field between neighbouring cells. This in turn can lead to non-physical oscillations and checkerboard-like patterns in the pressure field, also known as the checkerboard problem [1]. On structured Cartesian meshes, this problem can be circumvented by using a staggered grid arrangement in which the velocities at the cell faces are coupled to a compact-stencil gradient of pressure between directly neighbouring cell-centered nodes [10]. However, for many industrial applications the use of CFD involves complex geometries that require unstructured meshes. Extension of the staggered grid method to such cases is not straight-forward and leads to complex correction schemes and increased computational costs [11,12]. One commonly used strategy to address pressure field oscillations, is the use of the weighted interpolation method (WIM) to evaluate the velocity at the cell faces, in contrast to the direct interpolation method (DIM). This method was originally introduced as the pressure-weighted interpolation method and has allowed the wide-spread usage of the collocated grid arrangement [13,14]. Since its introduction this method has been extended in many ways, e.g. to account for decoupling caused by the use of small timesteps found in transient flows [15–18]and to account for under-relaxation factors [19–21]. Generalisation and unification of the aforementioned problems and accompanying solutions can be found in more recent works [22–24], which offer a comprehensible overview of these methods. One unavoidable consequence that these solutions have in common is the introduction of a numerical error in the form of a non-zero contribution to the kinetic energy balance of either the pressure term or the convective term [25]. Nevertheless, many commercially available and open-source codes favour stability at the cost of accuracy [26–29], by applying a dissipative form of the WIM. This avoids pressure field oscillations while at the same time increasing stability and accepting the consequential numerical error. Usually this is done implicitly through the use of a compact-stencil Laplacian, which directly couples node 𝑖to its neighbouring nodes 𝑗instead of its second-neighbours 𝑘. In this method, divergence-free values for the velocities at the faces follow directly from the Poisson equation, without introducing any numerical error. However, the coupling effect of the WIM is introduced implicitly by the Laplace operator and the numerical error is introduced to the collocated velocities in this case. This method introduces a dissipative pressure error, with the benefit of greatly reducing computational complexity and cost when using unstructured grids. However, with the increase in computational resources available and the ensuing rise of high-fidelity simulations, higher accuracy is desired to describe motion of fluids. Numerical dissipation is limiting this accuracy by disrupting fluid motion, especially at the smaller scales, which are essential to accurately depict turbulent flows, making the use of conservative schemes more important. The symmetry-preserving method eliminates numerical dissipation and conserves the physical properties of the flow, while at the same time warranting unconditional stability, by mimicking the properties of the continuous operators in their discrete counterparts [30]. This method has been extended to collocated grid arrangements [31]and implemented in the open-source code OpenFOAM [32]. The pressure error that is introduced by the application of this method on collocated grids remains as the largest source of numerical dissipation. Other methods to filter pressure oscillations without the introduction of numerical dissipation have been attempted, such as filtering pressure modes that lie on the kernel of the wide-stencil Laplacian operator [33]. However, through examination of the connection between the mesh and this kernel, it was shown that this method only works on Cartesian meshes [34], as the oscillatory part of the kernel vanishes for most unstructured meshes and complex geometries. A method that can eliminate, or at least balance, both the checkerboard problem and numerical dissipation is therefore still sought after. Moreover, the inadequacy of the kernel filtering method illustrated that a broader and clear definition of checkerboarding is lacking, and most published works in literature use a qualitative description of the phenomenon. The present work contributes to this topic by introducing a clear definition of the checkerboard problem. This definition can be used to quantify the level of checkerboarding on any geometry or mesh. Furthermore, a solver algorithm is developed which uses this quantification method to dynamically balance the occurrence of pressure oscillations and numerical dissipation. The structure of this Journal of Computational Physics 521 (2025) 113537 3 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. Fig. 2. Geometric parameters of the operators in Table 1, needed to establish the symmetry-preserving scheme. Table 1 Matrices containing the geometric parameters shown in Fig. 2, which form the building blocks for all symmetry-preserving operators. 𝑚 and 𝑛denote the number of faces and control volumes respectively. Operator Dimensions Description 𝑇𝑓𝑜,𝑇 𝑓𝑛 𝑛×𝑚face-owner and face-neighbour connectivity matrices, containing entry (𝑖, 𝑓) =1if cell 𝑖is connected to face 𝑓as an owner or as a neighbour respectively. 𝑁𝑠𝑚×3𝑚face-normal matrix with diagonal blocks (𝑁𝑠𝑥,𝑁 𝑠𝑦,𝑁 𝑠𝑧)containing the (𝑥, 𝑦, 𝑧)-components of the face-normals, 𝐧𝑓, respectively. 𝐴𝑠𝑚×𝑚diagonal matrix containing face areas, 𝐴𝑓. 𝛿𝑜 𝑛𝑠,𝛿𝑛 𝑛𝑠 𝑚×𝑚face-owner and face-neighbour normal distance diagonal matrices containing 𝛿𝑜 𝑛𝑓 and 𝛿𝑛 𝑛𝑓 respectively, denoting the absolute values of the 𝐧𝑓-projected vectors from face-centroid to owner or neighbour centroid respectively. Ω𝑐𝑛×𝑛cell-volume diagonal matrix. paper is as follows: Section 2gives an overview of the equations and different algorithms that are used and discusses the occurrence of checkerboarding and numerical errors. Section 3discusses methods to quantify checkerboarding and introduces a solver that dynamically balances numerical dissipation and checkerboard oscillations. Section 4shows results for the new solver compared to existing methods. Finally, section 5discusses the conclusions of the work and the outlook for future work. 2. Numerical framework 2.1. Symmetry-preserving method The semi-discretised formulation of equations (1)and (2)for arbitrary collocated grids, using the matrix-vector notation of [31], is given by: Ω𝜕𝑡𝐮𝑐+𝐶(𝐮𝑠)𝐮𝑐=−𝐷𝐮𝑐−Ω𝐺𝑐𝐩𝑐,(3) 𝑀𝐮𝑠=𝟎𝑐.(4) In three dimensions, the discrete collocated kinematic pressure and velocity fields are given by 𝐩𝑐∈ℝ𝑛and 𝐮𝑐=(𝐮𝑇 𝑐,𝑥,𝐮𝑇 𝑐,𝑦,𝐮𝑇 𝑐,𝑧)𝑇 ∈ ℝ3𝑛, whereas 𝐮𝑠∈ℝ𝑚gives the staggered velocities, in which 𝑛and 𝑚give the number of control volumes and faces respectively. A small set of matrices containing geometric information about the mesh, given in Table 1, is enough to derive all the matrix operators needed to from an algebraic symmetry-preserving framework, given in Table 2. Midpoint interpolation in the convective term is necessary to maintain the skew-symmetry of the continuous operator [30], whereas it was shown in [35–38]that employing volumetric interpolation in the other operators leads to unconditional stability, even on highly-distorted meshes. In the following text therefore, if superscript 𝛾is dropped, volumetric interpolation is employed, e.g. 𝐿𝑐= 𝑀𝑐𝐺𝑐=−𝑀Γ𝑐𝑠Ω−1Γ𝑇 𝑐𝑠𝑀𝑇=−𝑀Γ𝑉 𝑐𝑠Ω−1Γ𝑉𝑇 𝑐𝑠 𝑀𝑇. 2.2. Time-stepping algorithm The temporal integration, which was left undiscretised in equation (3), is taken care of by the fractional step method [39], in which a predictor velocity is calculated and corrected through solving a Poisson equation for pressure. Usage of a wide-stencil Laplacian for the pressure Poisson equation in combination with the DIM to calculate the staggered velocities leads to checkerboarding. By applying a compact stencil and/or the WIM, this method can be adjusted to deal with pressure field oscillations. Table 3gives an overview of the four possible fractional step methods found by combining the compactand wide-stencil methods with the DIM and Journal of Computational Physics 521 (2025) 113537 4 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. Table 2 Full set of matrix operators to form an algebraic symmetry-preserving framework. 𝑚 and 𝑛denote the number of faces and control volumes respectively. 𝛾indicates an unspecified interpolation method, for which 𝐿, 𝑀and 𝑉are given as options, indicating linear, midpoint and volumetric interpolation respectively. Operator definition Dimensions Description Ω=𝐼3⊗Ω𝑐3𝑛×3𝑛collocated volumes 𝛿𝑛𝑠 =𝛿𝑜 𝑛𝑠 +𝛿𝑛 𝑛𝑠 𝑚×𝑚face-normal distances Ω𝑠=𝛿𝑛𝑠𝐴𝑠𝑚×𝑚staggered volumes 𝑆𝑠=𝐴𝑠𝑁𝑠𝑚×3𝑚surface vectors 𝑀=(𝑇𝑓𝑜 −𝑇𝑓𝑛)𝐴𝑠𝑛×𝑚divergence 𝐺=𝛿−1 𝑛𝑠 (𝑇𝑓𝑛 −𝑇𝑓𝑜)𝑇=−Ω −1 𝑠𝑀𝑇𝑚×𝑛gradient 𝐿=𝑀𝐺 =−𝑀Ω−1 𝑠𝑀𝑇𝑛×𝑛compact-stencil Laplacian 𝑊𝛾 𝑜=⎧ ⎪ ⎨ ⎪ ⎩ 𝛿−1 𝑛𝑠 𝛿𝑛 𝑛𝑠,𝛾=𝐿 1 2𝐼𝑚,𝛾=𝑀 𝛿−1 𝑛𝑠 𝛿𝑜 𝑛𝑠,𝛾=𝑉 𝑚×𝑚interpolation weights 𝑊𝛾 𝑛=𝐼𝑚−𝑊𝛾 𝑜 Π𝛾 𝑐𝑠 =𝑊𝛾 𝑜𝑇𝑇 𝑓𝑜 +𝑊𝛾 𝑛𝑇𝑇 𝑓𝑛 𝑚×𝑛cell-to-face interpolator Γ𝛾 𝑐𝑠 =𝑁𝑠(𝐼3⊗Π𝛾 𝑐𝑠)𝑚×3𝑛cell-to-face dot-interpolator Γ𝛾 𝑠𝑐 =Ω −1Γ𝛾𝑇 𝑐𝑠 Ω𝑠3𝑛×𝑚face-to-cell interpolator 𝑀𝛾 𝑐=𝑀Γ𝛾 𝑐𝑠 𝑛×3𝑛collocated divergence 𝐺𝛾 𝑐=Γ 𝛾 𝑠𝑐 𝐺=−Ω −1𝑀𝛾𝑇 𝑐3𝑛×𝑛collocated gradient 𝐿𝛾 𝑐=𝑀𝛾 𝑐𝐺𝛾 𝑐=−𝑀Γ𝛾 𝑐𝑠Ω−1Γ𝛾𝑇 𝑐𝑠 𝑀𝑇𝑛×𝑛wide-stencil Laplacian 𝐶𝑐(𝐮𝑠)=𝑀diag(𝐮𝑠)Π𝑀 𝑐𝑠 𝑛×𝑛convective block 𝐶(𝐮𝑠)=𝐼3⊗𝐶 𝑐(𝐮𝑠)3𝑛×3𝑛convective operator 𝐷𝑐=−𝜈𝐿 𝑛 ×𝑛diffusive block 𝐷=𝐼3⊗𝐷 𝑐3𝑛×3𝑛diffusive operator Table 3 Overview of the four possible fractional step methods found by combining the compactand widestencil methods with the DIM and the WIM. Wide stencil Compact stencil DIM WIM DIM WIM 𝐮𝑝 𝑐=(𝐮𝑐,𝐮𝑠)𝐮𝑝∗ 𝑐=𝐮𝑝 𝑐−𝐺𝑐 𝐩𝑝 𝑐 𝐿𝑐 𝐩𝑛+1 𝑐=𝑀𝑐𝐮𝑝 𝑐𝐿 𝐩′ 𝑐=𝑀𝑐𝐮𝑝∗ 𝑐  𝐩𝑛+1 𝑐= 𝐩𝑝 𝑐+ 𝐩′ 𝑐 𝐮𝑛+1 𝑐=𝐮𝑝 𝑐−𝐺𝑐 𝐩𝑛+1 𝑐𝐮𝑛+1 𝑠=Γ 𝑐𝑠𝐮𝑝 𝑐−𝐺 𝐩𝑛+1 𝑐 𝐮𝑛+1 𝑠=Γ 𝑐𝑠𝐮𝑛+1 𝑐𝐮𝑛+1 𝑠=Γ 𝑐𝑠𝐮𝑝 𝑐−𝐺 𝐩𝑛+1 𝑐𝐮𝑛+1 𝑐=Γ 𝑠𝑐 𝐮𝑛+1 𝑠𝐮𝑛+1 𝑐=𝐮𝑝 𝑐−𝐺𝑐 𝐩𝑛+1 𝑐 the WIM. The temporal discretisation is not the main interest in this work and is simply denoted by a function, (𝐮𝑐,𝐮𝑠), which calculates the velocity predictor, 𝐮𝑝 𝑐. As an example, for Forward Euler time-stepping, this term is given by: 𝐮𝑝 𝑐=(𝐮𝑐,𝐮𝑠) =𝐮𝑛 𝑐+Δ𝑡𝑅 (𝐮𝑛 𝑐,𝐮𝑛 𝑠) =𝐮𝑛 𝑐−Δ𝑡Ω−1 (𝐶(𝐮𝑛 𝑠)+𝐷)𝐮𝑛 𝑐. (5) The value for the pressure predictor,  𝐩𝑝 𝑐is usually chosen to be 𝟎𝑐or  𝐩𝑛 𝑐, corresponding to the classical Chorin projection method [40]and the second-order Van Kan projection method [9], respectively. Note that this pressure predictor only has an effect when a compact-stencil Laplacian is used. For the wide-stencil methods, 𝐿𝑐 𝐩𝑝 𝑐could simply be added to both sides of the Poisson equation, resulting in 𝐿𝑐 𝐩𝑛+1 𝑐on the left hand side (LHS) and 𝑀Γ𝑐𝑠𝐮𝑝 𝑐on the right hand side (RHS). Upon closer inspection, the compact-stencil DIM can immediately be discarded for being less accurate than the compact-stencil WIM. This is true because the only difference between these methods is the back-and-forth interpolation of the predictor velocity. This operation can be viewed as the application of a Laplacian filter, as: 𝐮𝑛+1 𝑐=Γ 𝑠𝑐Γ𝑐𝑠𝐮𝑝 𝑐−𝐺𝑐 𝐩′ 𝑐=𝑓(𝐮𝑝 𝑐)−𝐺𝑐 𝐩′ 𝑐. This filter can be understood most easily by considering its effect on a uniform Cartesian mesh, where the classical [1,−2,1] coefficients of the Laplacian are retrieved: 𝑓(𝐮𝑝 𝑐)=𝐿𝑓𝐮𝑝 𝑐(uniform Cartesian),(6) 𝐿𝑓=𝐼+1 4diag(𝐿𝑓,𝑥,𝐿 𝑓,𝑦,𝐿 𝑓,𝑧),(7) [𝐿𝑓,𝑥𝐮𝑝 𝑐,𝑥]𝑖,𝑗,𝑘 =[𝐮𝑝 𝑐,𝑥]𝑖−1,𝑗,𝑘 −2[𝐮𝑝 𝑐,𝑥]𝑖+[𝐮𝑝 𝑐,𝑥]𝑖+1,𝑗,𝑘 ,(8) Journal of Computational Physics 521 (2025) 113537 5 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. where subscripts 𝑖, 𝑗, 𝑘indicate cell numbering in 𝑥, 𝑦, 𝑧-directions respectively. This unnecessary extra filtering is smoothing the velocity field and thereby causing undesired numerical dissipation. This method is therefore not considered any further. The remaining methods each have their own problem, which makes choosing the right method a trade-off between these factors. As discussed in the introduction, the wide-stencil DIM is very prone to the occurrence of checkerboarding. The remainder of this section discusses occurrence of checkerboarding for the WIMs, and the numerical errors that they introduce when trying to diminish the problem. As there is no numerical error for the wide-stencil DIM of this category, this method is discussed no more hereafter. In the wide-stencil WIM, a correction is applied to the staggered velocities, which makes them non-divergence-free as a result. Although carried out in exactly the same way, the calculation of the staggered velocities for the compact-stencil WIM is without correction, and follows directly from the Poisson equation. In this method, however, the correction is applied to the collocated velocities, which in turn makes them non-divergence-free. The divergence of either the staggered or collocated velocities introduces a numerical error to the simulation, which is further discussed in section 2.4. 2.3. Occurrence of checkerboarding Although the WIM, for both the compactand wide-stencil methods, introduces a coupling of the pressure field between neighbouring nodes at the cost of a numerical error, these methods can still show oscillating pressure fields, most commonly caused by the usage of a small time-step in unsteady simulations [25], or by the inclusion of a pressure predictor as  𝐩𝑝 𝑐= 𝐩𝑛 𝑐in case of the compact-stencil WIM [32]. To illustrate this, note that if Δ𝑡 →0+, then 𝐮𝑝 𝑐→𝐮𝑛 𝑐in equation (5), since the effect of the convective and diffusive terms are proportional to Δ𝑡. This removes the coupling established in the wide-stencil WIM, since it is dependent on 𝐮𝑠 and the convective term. To show the decoupling for the compact-stencil WIM, first regard the case of  𝐩𝑝 𝑐=𝟎𝑐, with Forward Euler temporal discretisation as an example, which leads to: 𝐮𝑝∗ 𝑐=𝐮𝑛 𝑐,(9) 𝐮𝑛 𝑐=𝐮𝑛−1 𝑐−𝐺𝑐 𝐩𝑛 𝑐=𝐮0 𝑐−𝐺𝑐 𝑛 ∑ 𝑖  𝐩𝑖 𝑐,(10) 𝐿 𝐩𝑛+1 𝑐=𝑀𝑐𝐮𝑝∗ 𝑐=𝑀𝑐𝐮0 𝑐−𝐿𝑐 𝑛 ∑ 𝑖  𝐩𝑖 𝑐,(11) 𝐿ℙ𝑛+1 𝑐=𝑀𝑐𝐮0 𝑐+(𝐿−𝐿𝑐)ℙ𝑛 𝑐,(12) where ℙ𝑛 𝑐=∑𝑛 𝑖 𝐩𝑖 𝑐and 𝐿ℙ𝑛 𝑐is added on both sides of equation (11)to reach equation (12). Equation (12)gives a stationary iterative method to solve the wide-stencil Poisson equation, to which the solution is a decoupled pressure field. The coupling that the compactstencil Laplacian gave is therefore lost if the time-step becomes too small. Using the Van Kan method, where  𝐩𝑝 𝑐= 𝐩𝑛 𝑐, leads to more severe oscillations in the pressure field. In this case the wide-stencil Laplacian operates on a larger part of the pressure field: 𝐿 𝐩′ 𝑐=𝑀𝑐𝐮𝑝∗ 𝑐=𝑀𝑐𝐮𝑝 𝑐−𝐿𝑐 𝐩𝑝 𝑐,(13) which weakens the coupling that the compact-stencil Laplacian provided. Since the pressure oscillations are “invisible” to the wide-stencil gradient, the collocated velocities remain unaffected and the algorithm advances as if they were not there. The oscillations are retained in this case and could grow until they reach a point in which they cause numerical issues and unstable solutions. Although their retention is evident, their actual origins and method of growth are not entirely understood. In [31]it was argued that the convective term causes spurious modes in the velocity field that lead to checkerboarding, but an analysis of this process was not performed. A relation to the method with which the Poisson equation is solved, including preconditioners, was discussed in [41], but the exact mechanisms remain unclear. 2.4. Numerical errors of the WIM To analyse the numerical errors of the WIM, the global discrete kinetic energy is considered, given by 𝐸𝐾=1 2𝐮𝑇 𝑐Ω𝐮𝑐. The temporal evolution of this term is derived using the product rule and equation (3): 𝑑 𝑑𝑡𝐸𝐾=−1 2⎛⎜⎜⎝ 𝐮𝑇 𝑐(𝐶(𝐮𝑠)+𝐶𝑇(𝐮𝑠))𝐮𝑐 +𝐮𝑇 𝑐(𝐷+𝐷𝑇)𝐮𝑐 +𝐮𝑇 𝑐Ω𝐺𝑐𝐩𝑐+𝐩𝑇 𝑐𝐺𝑇 𝑐Ω𝑇𝐮𝑐⎞⎟⎟⎠ .(14) If 𝐶(𝐮𝑠)is symmetry-preserving, i.e. 𝐶(𝐮𝑠) =−𝐶𝑇(𝐮𝑠), mimicking the continuous operator, then the contribution of the convective term equals zero. Similarly, the contribution of the pressure term equals zero if Ω𝐺𝑐=−𝑀𝑇 𝑐and 𝑀𝑐𝐮𝑐=𝟎𝑐. The global kinetic energy dissipation is therefore reduced to the viscous term, equaling: 𝑑 𝑑𝑡𝐸𝐾=−1 2𝐮𝑇 𝑐(𝐷+𝐷𝑇)𝐮𝑐≤0,(15) which is strictly dissipative if 𝐷is constructed using the symmetry-preserving discretisation, retaining the positive-definiteness of the continuous operator. Journal of Computational Physics 521 (2025) 113537 6 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. 2.4.1. Wide-stencil WIM -convective error For the wide-stencil WIM, 𝑀𝑐𝐮𝑛+1 𝑐=𝟎𝑐and therefore the pressure error equals zero. Conversely, the convective error is non-zero. This is caused by the correction to the staggered velocities and the resulting fact that 𝑀𝐮𝑛+1 𝑠≠0, since it is given by: 𝑀𝐮𝑛+1 𝑠=𝑀(Γ𝑐𝑠𝐮𝑝 𝑐−𝐺 𝐩𝑛+1 𝑐) =  : 𝟎𝑐 𝑀Γ𝑐𝑠 (𝐮𝑝 𝑐−𝐺𝑐 𝐩)+𝑀(Γ𝑐𝑠Γ𝑠𝑐 −𝐼)𝐺 𝐩𝑛+1 𝑐 =(𝐿𝑐−𝐿) 𝐩𝑛+1 𝑐. (16) To see the effect on the evolution of kinetic energy, recall that 𝐶𝑐(𝐮𝑠) =𝑀diag(𝐮𝑠)Π𝑚 𝑐𝑠. The appearance of 𝑀in this definition and the incompressibility constraint of equation (4), is related to the simplified advective form in which equation (1)is written. The off-diagonals of 𝐶𝑐(𝐮𝑠)contain the fluxes between cells, and are skew-symmetric by construction. In contrast the diagonal entries are given by: 𝐶𝑖,𝑖(𝐮𝑠)= 1 2[𝑀𝐮𝑠]𝑖,(17) which are non-zero in this case. Equation (17)is consistent with the situation for compressible flows, see for example [42]. However, in the incompressible framework, this term is not balanced by the mass conservation equation. If 𝐶𝑐(𝐮𝑠)is split into its non-zero diagonal and skew-symmetric off-diagonals as 𝐶𝑐(𝐮𝑠) =1 2diag(𝑀𝐮𝑠)+𝐶𝑂𝐷 𝑐(𝐮𝑠), then the numerical error of the evolution of the global kinetic energy at time-step 𝑛 +1due to the convective error is given by: 𝑑 𝑑𝑡𝐾=− 1 2𝐮𝑛+1𝑇 𝑐⎛⎜⎜⎜⎝ : 𝟎3𝑛×3𝑛 𝐶𝑂𝐷(𝐮𝑛+1 𝑠)+𝐶𝑂𝐷𝑇(𝐮𝑛+1 𝑠)⎞⎟⎟⎟⎠ 𝐮𝑛+1 𝑐 −1 2𝐮𝑛+1𝑇 𝑐[𝐼3⊗diag(𝑀𝐮𝑛+1 𝑠)]𝐮𝑛+1 𝑐. (18) The convective error in the wide-stencil WIM is therefore caused by the divergence of the staggered velocities and proportional to Δ𝑡 (𝐿𝑐−𝐿)𝐩𝑛+1 𝑐. 2.4.2. Compact-stencil WIM - pressure error On the other hand, for the compact-stencil WIM, 𝑀𝐮𝑛+1 𝑠=𝟎𝑐and therefore the convective error equals zero. However, in this case the correction of the collocated velocities leads to 𝑀𝑐𝐮𝑛+1 𝑐≠𝟎𝑐, which is given by: 𝑀𝑐𝐮𝑛+1 𝑐=𝑀Γ𝑐𝑠 (𝐮𝑝 𝑐−𝐺𝑐 𝐩′ 𝑐) =  : 𝟎𝑐 𝑀(Γ𝑐𝑠𝐮𝑝 𝑐)−𝐺 𝐩′ 𝑐+𝑀(𝐼−Γ 𝑐𝑠Γ𝑠𝑐 )𝐺 𝐩′ 𝑐 =(𝐿−𝐿𝑐) 𝐩′ 𝑐. (19) Then the ensuing numerical error in the evolution of the global kinetic energy at time-step 𝑛 +1due to the pressure error is given by: 𝑑 𝑑𝑡𝐾=−𝐮𝑛+1𝑇 𝑐Ω𝐺𝑐𝐩′ 𝑐 =𝐩′𝑇 𝑐𝑀𝑐𝐮𝑛+1 𝑐 =𝐩′𝑇 𝑐(𝐿−𝐿𝑐) 𝐩′ 𝑐. (20) In the compact-stencil WIM, the pressure error is caused by the divergence of the collocated velocities and is, similarly to the convective error, proportional to Δ𝑡 (𝐿−𝐿𝑐)𝐩′ 𝑐. In practice, most numerical codes revert to the compact-stencil WIM. The main reasons for this are threefold: (i) The inclusion of a pressure predictor as  𝐩𝑝 𝑐= 𝐩𝑛 𝑐reduces the order of the pressure error from  (Δ𝑡2)to  (Δ𝑡4)since  𝐩𝑛+1 𝑐∼ (Δ𝑡)and ( 𝐩𝑛+1 𝑐− 𝐩𝑛 𝑐)∼ (Δ𝑡2), making the non-physical contribution quite small. This reduction is not possible with the wide-stencil Laplacian. (ii) The compact-stencil Laplacian reduces the computational complexity and cost when using unstructured grids. And finally, (iii) the convective error term, as shown in equation (18), is not strictly dissipative. Since the global divergence of the staggered velocities will always be equal to zero, there will always be cells with negative divergence to compensate cells that have positive divergence. This means that kinetic energy can be added in some parts of the solution, leading to instabilities. Conversely, the resulting term of equation (20)will be dependent on the term (𝐼−Γ 𝑐𝑠Γ𝑠𝑐 ), which in turn depends on the choice in interpolators and meshing. These factors are much easier to control, and moreover, it was shown in [35–38]that choosing volumetric interpolation for this operator causes this term to be strictly dissipative on circumcenter meshes. In most practical meshes, stability can be warranted in this way. Journal of Computational Physics 521 (2025) 113537 7 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. 3. Defining the checkerboard problem Past works that extensively discuss the topic use different names for the checkerboard problem, such as checkerboarding [43,17, 44], spurious pressure modes [22], odd-even decoupling [45], or zigzagness [46], but none of them give a quantifiable definition for the problem. In those works a constant and uniform solution is applied that does not need quantification of the problem, whereas in this work a solution is sought after that is applied proportionally to the severity of the problem. For this reason a clear quantification method has to be introduced first. 3.1. Using the kernel of the discrete Laplacian to define checkerboarding One way to define checkerboard modes is by performing an eigenvector decomposition of 𝐿𝛾 𝑐and defining the checkerboard modes as 𝐩− 𝑐∈𝐾𝑒𝑟(𝐿𝛾 𝑐). Since 𝐿𝛾 𝑐=𝑀𝛾 𝑐𝐺𝛾 𝑐=𝑀Γ𝛾 𝑐𝑠Γ𝛾 𝑠𝑐𝐺, the 𝐩− 𝑐modes will depend on the choice of interpolator and the mesh. The constant mode vector will always lie on the kernel, in addition to any spurious mode vectors. This method was for instance used by [33], where the method was adequate since midpoint interpolation and Cartesian meshes were used. By looking at the definitions of 𝑀𝛾 𝑐and 𝐺𝛾 𝑐it becomes clear why this null-space exists for these specific conditions, since: [Γ𝛾 𝑠𝑐𝐺𝜙 𝜙 𝜙𝑐]𝑖=1 Ω𝑖∑ 𝑓∈𝐹𝑓(𝑖) 𝑤𝛾 𝑓𝑖Ω𝑓 𝜙 𝜙 𝜙𝑗−𝜙 𝜙 𝜙𝑖 𝛿𝑛𝑓 𝐧𝑓(𝑖) =1 Ω𝑖∑ 𝑓∈𝐹𝑓(𝑖) 𝑤𝛾 𝑓𝑖𝜙 𝜙 𝜙𝑗𝐬𝑓(𝑖)−1 Ω𝑖∑ 𝑓∈𝐹𝑓(𝑖) 𝑤𝛾 𝑓𝑖𝜙 𝜙 𝜙𝑖𝐬𝑓(𝑖), (21) gives the wide-stencil gradient at cell 𝑖. Where 𝜙 𝜙 𝜙𝑐∈ℝ𝑛×1 and 𝜙 𝜙 𝜙𝑖=[𝜙 𝜙 𝜙𝑐]𝑖. Cell 𝑗is neighbouring cell 𝑖through face 𝑓. 𝐹𝑓(𝑖)denotes the set of faces that define cell 𝑖. Ω𝑖=[Ω𝑐]𝑖,𝑖, Ω𝑓=[Ω𝑠]𝑓,𝑓 and 𝛿𝑛𝑓 =[𝛿𝑛𝑠]𝑓. 𝑤𝑖𝑓 equals [𝑊𝛾 𝑜]𝑓or [𝑊𝛾 𝑛]𝑓depending on whether 𝑖is the owner or the neighbour of face 𝑓, respectively. Finally, 𝐧𝑓(𝑖)and 𝐬𝑓(𝑖)respectively give the outward-pointing face-normal vector and outward-pointing surface vector at face 𝑓with respect to cell 𝑖, note that 𝐬𝑓(𝑖)=𝐴𝑓𝐧𝑓(𝑖). The second term of the RHS of equation (21) vanishes in case 𝑤𝛾 𝑓𝑖 is constant over all faces of 𝑖, since ∑𝑓∈𝐹𝑓(𝑖)𝐬𝑓(𝑖)=𝟎. This is the case for midpoint interpolation and for uniform Cartesian meshes. Similarly for the divergence operator: [𝑀Γ𝛾 𝑐𝑠𝜓 𝜓 𝜓𝑐]𝑖=∑ 𝑓∈𝐹𝑓(𝑖)(𝑤𝛾 𝑖𝑓 𝜓 𝜓 𝜓𝑖+𝑤𝛾 𝑗𝑓𝜓 𝜓 𝜓𝑗)⋅𝐬𝑓(𝑖) =∑ 𝑓∈𝐹𝑓(𝑖) 𝑤𝛾 𝑖𝑓 𝜓 𝜓 𝜓𝑖⋅𝐬𝑓(𝑖)+∑ 𝑓∈𝐹𝑓(𝑖) 𝑤𝛾 𝑗𝑓𝜓 𝜓 𝜓𝑗⋅𝐬𝑓(𝑖), (22) in which, again, the first term of the RHS vanishes for midpoint interpolation and for uniform Cartesian meshes. When combining equations (21)and (22)it becomes evident that 𝐿𝛾 𝑐does not connect cell 𝑖to its direct neighbours 𝑗, but only to its second-neighbours 𝑘, when midpoint interpolation is used: [𝐿𝑀 𝑐]𝑖,𝑗 =0,[𝐿𝑀 𝑐]𝑖,𝑘 =𝐴𝑓𝐴𝑔 4Ω𝑗 𝐬𝑓(𝑖)⋅𝐬𝑔(𝑗),(23) in which face 𝑓lies between cells 𝑖and 𝑗, and face 𝑔lies between cells 𝑗and 𝑘. If this odd-even parity is sustained throughout the entire mesh, two disconnected groups of cells will exist. In addition to the constant kernel vector, this gives rise to a spurious kernel vector: [𝐩− 𝑐]𝑖=(−1) 𝑝(𝑖),(24) where 𝑝(𝑖)denotes the parity of cell 𝑖, equal to 0 or 1. In the special case of Cartesian meshes, the dot product in equation (23)will be zero for any diagonal second-neighbour pairing, giving rise to 2𝑁𝑑𝑖𝑚 kernel vectors, with number of dimensions 𝑁𝑑𝑖𝑚. For example, in three-dimensional Cartesian meshes the resulting set of eight vectors is given by [33]: [𝐩− 𝑐(𝐼𝐽𝐾)]𝑖,𝑗,𝑘 =(−1) 𝑖𝐼+𝑗𝐽+𝑘𝐾 ,(25) in which 𝐼, 𝐽, 𝐾∈{0, 1} and indices 𝑖, 𝑗, 𝑘indicate the cell numbering in each of the Cartesian directions. In this notation, 𝐩− 𝑐(000) gives the constant vector. [𝐿𝛾 𝑐]𝑖,𝑗 is generally only zero for midpoint interpolation. Whereas, in general, [𝐿𝛾 𝑐]𝑖,𝑗 ≠0for linear or volumetric interpolation, except on uniform meshes where all interpolators equal midpoint interpolation. However, a set of kernel vectors was derived by [34]for any Cartesian mesh with midpoint, linear or volumetric interpolation. To do so, the fact that 𝐾𝑒𝑟(𝐺𝛾 𝑐) ∈𝐾𝑒𝑟(𝐿𝛾 𝑐) was used, in addition to rewriting 𝐺𝛾 𝑐to 𝐺𝐺Π𝛾 𝑐𝑠. 𝐺𝐺denotes a Gauss-gradient and the overbar denotes swapping the interpolation weights between owners and neighbours of the faces, such that Π𝑀 𝑐𝑠 = Π𝑀 𝑐𝑠 , Π𝑉 𝑐𝑠 = Π𝐿 𝑐𝑠 and Π𝐿 𝑐𝑠 = Π𝑉 𝑐𝑠. For details on this equality, the reader is referred to Appendix A. By taking the interpolation before the Gauss-gradient, cell-centered values can be chosen such that sets of opposing interpolated face values are always equal, for which the cell-centered Gaus-gradient equals zero. The resulting set of eight vectors for three-dimensional Cartesian meshes is then given by: Journal of Computational Physics 521 (2025) 113537 8 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. Fig. 3. Pressure oscillation in a one-dimensional periodic domain, orthogonal to 𝐾𝑒𝑟(𝐿𝑐),butdetectableby𝐶𝑐𝑏. [𝐩− 𝑐(𝐼𝐽𝐾,𝛾)]𝑖,𝑗,𝑘 =(−1) 𝑖𝐼+𝑗𝐽+𝑘𝐾 ([Δ𝑥]𝐼 𝑖[Δ𝑦]𝐽 𝑗[Δ𝑧]𝐾 𝑘)𝛼 ,(26) in which 𝛼={−1, 0, 1} for linear, midpoint and volumetric interpolations respectively. This set is not necessarily mutually orthogonal, especially for linear and volumetric interpolations, however, in all cases the set is linearly independent and therefore spans the nullspace of 𝐿𝛾 𝑐. The derivation of this set of vectors with an example is given in Appendix B. Despite this extension, calculating the kernel of 𝐿𝛾 𝑐for non-Cartesian meshes involves performing a singular value decomposition, for which the computational cost grows exponentially with the number of grid points as (𝑁3 𝑔𝑟𝑖𝑑), which quickly becomes unaffordable [47]. Moreover, for most meshes the rank of the kernel reduces to one if the parity and the orthogonal faces disappear, which is nearly always the case for any unstructured mesh, leaving only the constant mode vector. Therefore using the kernel of 𝐿𝛾 𝑐to define and quantify the checkerboard problem is often insufficient, leading to the search for a more generally applicable definition. 3.2. Applying a general definition of checkerboarding In the previous section a very restrictive definition of checkerboarding was used, for which some easy examples were given in which case the definition would not be fruitful. A more useful definition can be found if the problem is regarded at control volume level. The essential problem of the decoupled control volume is that the pressure gradient over a given face, [𝐺𝐩𝑐]𝑓, might give a significant non-zero value, while the value at the adjacent cell [𝐺𝑐𝐩𝑐]𝑖can be (close to) zero. This problem might occur only in a few cells and therefore lie mostly outside of the kernel of 𝐿𝑐, that is, if the kernel even contains spurious vectors. It turns out that the ratio between the 𝐿2norms of vectors 𝐺𝑐𝐩𝑐and 𝐺𝐩𝑐gives a good global indication of this decoupling. The expression for the 𝐿2norms of cell-centered and face-centered fields is given by ‖𝑎𝑐‖=𝑎𝑇 𝑐Ω𝑐𝑎𝑐and ‖𝑎𝑠‖=𝑎𝑇 𝑠Ω𝑠𝑎𝑠, respectively. Furthermore, in the lower limit the field lies fully inside the kernel, resulting in a ratio of zero, whereas a perfectly smooth field gives the upper limit, resulting in a ratio of one. Since a higher coefficient should indicate more prevalent checkerboarding, the checkerboard coefficient is finally found by subtracting this ratio from one, resulting in: 𝐶𝑐𝑏 (𝐩𝑐)=1−‖𝐺𝑐𝐩𝑐‖ ‖𝐺𝐩𝑐‖=1−𝐩𝑇 𝑐𝐺𝑇 𝑐Ω𝐺𝑐𝐩𝑐 𝐩𝑇 𝑐𝐺𝑇Ω𝑠𝐺𝐩𝑐 =𝐩𝑇 𝑐(𝐿−𝐿𝑐)𝐩𝑐 𝐩𝑇 𝑐𝐿𝐩𝑐 .(27) 𝐶𝑐𝑏 (𝐩𝑐)is defined as zero for the constant pressure field in which case ‖𝐺𝐩𝑐‖ =0. As the 𝐿2norms of both pressure gradient fields are non-negative, the coefficient’s upper bound is 1. The lower bound is 0 as long as 𝐩𝑇 𝑐(𝐿−𝐿𝑐)𝐩𝑐is non-positive, i.e. 𝐿 −𝐿𝑐is negative semidefinite. This holds true for volumetric interpolation and circumcenter meshes [38]. For virtually all simulations in practice, the eigenvalues of 𝐿 −𝐿𝑐were non-positive or at least very close to zero and non-problematic, as long as volumetric interpolation was used, leading to a value of 𝐶𝑐𝑏 in the interval [0,1]. The checkerboard coefficient can be calculated for any collocated scalar field. As pressure field oscillations are the focus of this work, 𝐶𝑐𝑏 (𝐩𝑐)will henceforth simply be denoted by 𝐶𝑐𝑏. Aside from being a non-dimensional and normalised coefficient, rewriting the coefficient as done in equation (27) reveals some other properties which support this definition. Firstly, it reveals that the magnitude of the pressure fields does not matter, since the constant pressure vector lies inside the kernel of 𝐿and 𝐿𝑐and will therefore be filtered in the operation of equation (27). This should always be the case for any definition that is employed, since the value of pressure is not what matters in incompressible flows, rather its gradient. Secondly, the matrix (𝐿−𝐿𝑐)=𝑀(𝐼−Γ 𝑐𝑠Γ𝑠𝑐 )𝐺is innately linked to the fundamental problem of the collocated grid arrangement, which is that the gradient of pressure and the velocities are not defined at the same location. The ensuing interpolators will always apply some smoothing to a field such that only 𝐼≈Γ 𝑐𝑠Γ𝑠𝑐 . Thirdly, this definition has a physical meaning, since the numerical errors in the convective and pressure terms, as discussed in sections 2.4.1 and 2.4.2 respectively, are proportional to the term ± (𝐿−𝐿𝑐)𝐩𝑐. For the wide-stencil WIM, (𝐿𝑐−𝐿)𝐩 =𝑀𝐮𝑠, which is non-zero and leads to the error expressed in equation (18). Similarly, for the compact-stencil WIM, (𝐿𝑐−𝐿)𝐩 =𝑀𝑐𝐮𝑐, which is non-zero and leads to the pressure error as seen in equation (20). One notable difference is that 𝐶𝑐𝑏 is independent of the time-step, which is also a necessary property, since it should be possible to calculate the coefficient without any knowledge of the temporal discretisation. Finally, 𝐶𝑐𝑏 is able to detect local oscillations which lie outside of the kernel of 𝐿𝑐, which it should be able to do, as mentioned before. A simple example to illustrate this point can be provided with an oscillation on a one-dimensional periodic domain, as seen in Fig. 3. The following calculations apply to this configuration: Journal of Computational Physics 521 (2025) 113537 9 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. 𝐩𝑐=[0010 10 ],(28) 𝐩− 𝑐(1,𝛾)=[111111],(29) 𝐺𝐩𝑐=[01 111 0 ],(30) 𝐺𝑐𝐩𝑐=[01 2010 1 2],(31) 𝐩𝑇 𝑐𝐩− 𝑐(1,𝑣)=0,(32) 𝐶𝑐𝑏 =1 − ‖𝐺𝑐𝐩𝑐‖ ‖𝐺𝐩𝑐‖=5 8,(33) where column vectors are represented horizontally for readability. Intuitively, this wiggle should not have a zero value when quantifying the checkerboard problem, therefore 𝐶𝑐𝑏 gives a more desirable outcome than the kernel vector method. 3.3. One possible application of the checkerboard coefficient Since 𝐶𝑐𝑏 gives a global, non-dimensional, normalised coefficient for checkerboarding, it can directly be applied in the solver algorithm to diminish the occurrence of the problem. It could, for example, be used to shift between the possible fractional step method algorithms displayed in Table 3, where choices are made concerning: (i) the width of the Laplacian stencil, (ii) the interpolation method for the velocity correction and (iii) the calculation of the pressure predictor. Since most collocated finite volume fractional step method codes apply the compact-stencil WIM, and since the inclusion of the pressure predictor is a known cause of checkerboarding, one way to use 𝐶𝑐𝑏 is to deliver a negative feedback through this predictor value. To this end, the following expression for 𝐩𝑝 𝑐is used: 𝐩𝑝 𝑐=𝜃𝑝𝐩𝑛 𝑐,(34) 𝜃𝑝=1−𝐶𝑐𝑏 =𝐩𝑛𝑇 𝑐𝐿𝑐𝐩𝑛 𝑐 𝐩𝑛𝑇 𝑐𝐿𝐩𝑛 𝑐 .(35) By doing so, the solver will converge to 𝜃𝑝=1in absence of checkerboarding in the pressure field, benefiting from the lower numerical dissipation this offers. Whereas if the case or mesh is prone to checkerboarding, the algorithm will tend to 𝜃𝑝=0, in which case the increase in numerical dissipation can damp the oscillations. This feedback onto the pressure predictor offers a dynamical balance between two problems, so that the user of the code does not have to decide which problem will be more prevalent, creating a unified solver that should perform well in both situations. For this initial attempt at establishing such a negative feedback, a simple linear relation is used to explore its effects, other relations might also work and could give better results. 4. Results In this section, the new solver with the predictor pressure given by equations (34)and (35)is tested and compared to the Chorin and Van Kan methods, which use 𝜃𝑝=0and 𝜃𝑝=1, respectively. The values of 𝜃𝑝for the resulting three methods, labelled 𝜃0, 𝜃1 and 𝜃𝑝, are given in Table 4. These three methods are tested and compared in this section. To do so, firstly, numerical dissipation is measured using a two-dimensional Taylor-Green vortex. Secondly, a temporal and spatial convergence study is performed using a two-dimensional lid-driven cavity. And finally, a turbulent channel flow case is used to measure overall accuracy of the solver and to study the behaviour of the checkerboard coefficient in unsteady flows. All cases were run with OpenFOAM, using the symmetrypreserving spatial discretisation and Runge-Kutta 3 temporal integration, which are both implemented in the RKSymFoam solver developed for [32,48]. Table 4 Settings of the tested solvers. 𝜃0𝜃1𝜃𝑑𝑦 𝜃𝑝011−𝐶𝑐𝑏 In each case the newly implemented solvers are compared to a readily implemented solver available in OpenFOAM, which is either icoFoam for laminar flows or pisoFoam for turbulent flows [28]. Apart from some choices in spatial discretisation, the biggest difference of the OpenFOAM solvers is how they remove time-step dependency of the compact-stencil WIM, which is similar to and derived from the methods used by [16,21]. For explicit time-stepping it is given by: 𝐮𝑝 𝑠=Γ 𝛾 𝑐𝑠𝐮𝑝 𝑐−𝐶𝑢𝑠𝐮𝑠,𝑐𝑜𝑟𝑟,(36) 𝐮𝑠,𝑐𝑜𝑟𝑟 =𝐮𝑛 𝑠−Γ 𝛾 𝑐𝑠𝐮𝑛 𝑐,(37) [𝐶𝑢𝑠]𝑓,𝑓 =1−𝑚𝑖𝑛 (||[𝐮𝑠,𝑐𝑜𝑟𝑟]𝑓|| ||[𝐮𝑛 𝑠]𝑓||+𝜖,1),(38) where 𝜖is a very small number to avoid division by zero. This term stabilises the results at the cost of introducing a sizable amount of numerical dissipation, which was shown by [49]. Journal of Computational Physics 521 (2025) 113537 16 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. 𝐷𝜖 𝑘=−𝜈⟨𝜕𝑗𝑢𝑖𝜕𝑗𝑢𝑖⟩+𝜈⟨𝜕𝑗𝑢𝑖⟩⟨𝜕𝑗𝑢𝑖⟩, 𝐶𝑇 𝑘=𝐶𝑘−𝐶𝑃 𝑘+⟨𝑢𝑗⟩𝜕𝑗𝑘, 𝐷𝑣 𝑘=𝐷𝑘−𝐷𝜖 𝑘,(48) such that the only choice in discretisation that remains to be made, is how to take a cell-centered gradient of the velocity, 𝜕𝑗𝑢𝑖, as ⟨𝑢𝑗⟩𝜕𝑗𝑘 =0for the channel configuration. In this work, the cell-centered gradient is taken for each component of 𝐮𝑐separately using 𝐺𝑐, creating a tensor field with 9 components at the cell-centers. The derivation of the equations and the implementation of the method in OpenFOAM can be found in the GitHub repository of [55]. The implemented methods rely on a Runge-Kutta based solver, whereas pisoFoam relies on the temporal integration methods implemented in OpenFOAM. To make a fair comparison between these methods the Backward Euler temporal integration was chosen, since it is one of few schemes available to both solvers. Therefore, the results on the coarse mesh, shown in Figs. 12a through 13e, contain a fairly large error, which is mainly caused by the Backward Euler time scheme. The emphasis of these figures, however, is to show the difference between the implemented methods and pisoFoam, not the overall quality of the solution. To show the influence of the temporal error, the results for the 𝜃0method with a Runge-Kutta 3 scheme are shown in Figs. 12a, 12band 12c. From these figures it can also be seen that even on a coarse mesh, the implemented solvers are able to attain fairly decent results. From Fig. 12a it can be seen that there the mean stream-wise velocity for all solvers is generally much too high. pisoFoam shows the highest mean stream-wise velocity, indicating that it is less able to develop turbulence on this mesh, which may be due to its dissipative nature. This finding is confirmed by Fig. 12b, in which the viscous shear stress, 𝜌𝜈 𝑑⟨𝑢𝑥⟩ 𝑑𝑦 , and the Reynolds stress, 𝜌⟨𝑢′ 𝑥𝑢′ 𝑦⟩are plotted. The highest value of the viscous stress and the lowest value of the Reynolds stress for the pisoFoam solver confirms that this solver is least able to develop turbulence on this specific mesh. The root-mean-square (RMS) velocities are depicted in Fig. 12c, which shows that the values of 𝑅𝑀𝑆(𝑢𝑥)are too high and the RMS velocities in other directions are too low, the least accurate solution again given by pisoFoam. These inaccuracies can be caused by coarse mesh resolution, especially in the stream-wise direction. The turbulent kinetic energy budgets are given in Figs. 13a through 13e. In general, the implemented solvers show behaviour resembling the reference solution, however, because of the coarse mesh they have not attained an accurate solution. The lower turbulence of pisoFoam can be seen in these figures by the lower absolute values of the peaks in production, transport, viscous diffusion and pressure diffusion. The dissipation budget also shows lower absolute values, and a lesser expression of the characteristic bump in the line around 𝑦+=15. For all solvers, the location of the peak in production seems to occur more or less at the location where the Reynolds shear stress and the viscous stress are equal, as is expected [56]. The value of the peak is slightly too high or too low for the implemented methods, with the 𝜃0method having the lowest production of the three. pisoFoam has the lowest production, in line with what was seen before. The inaccuracies in production seem to be mainly compensated by the viscous diffusion around the peak and the dissipation closer to the wall. The under-estimation of the transport and viscous diffusion terms further away from the wall also seem to be compensated by the dissipation, of which the values are too high in this region. The characteristic bump in the dissipation term is not well resolved on this mesh, as the profiles are too smooth. Finally, the pressure diffusion term shows relatively large deviations from the reference value, since pressure diffusion is small in general, the deviations appear more visibly. In general, the implemented solvers show similar results, and the deviations from the reference values are expected to be caused by the coarse mesh resolution. Finally, the results for the 𝑋180 mesh using the Runge-Kutta 3 temporal integration method are shown in Figs. 14a through 14c. Fig. 14a shows that the mean stream-wise velocity is very accurately predicted, indicating that the turbulence is properly developed. The RMS velocities also depict this, although a slight deviation can still be seen in the values of 𝑅𝑀𝑆(𝑢𝑥). This can be explained by some under-resolution in the stream-wise direction, since the mesh is not of full direct numerical simulation (DNS) quality. Finally, Fig. 14c shows that the implemented solvers converge to one another and lie very closely to the reference values. Convergence to the reference is not yet attained for all values, but the right trend is visible and suggests that additional mesh refinement would further increase convergence of the results. Especially mesh refinement in the stream-wise direction is needed to improve the 𝑅𝑀𝑆(𝑢𝑥) values. Overall, with mesh refinement, the implemented solvers seem to converge to the same solution, and checkerboarding seems to diminish. The range in which the 𝜃𝑑𝑦 method operates is therefore also smaller. On coarse meshes and in laminar flows however, the method still has significant benefits and its ability to balance dynamically between adding numerical dissipation or allowing checkerboarding. In flows that show transition from laminar flow to turbulence over time, the method should perform well with respect to the other methods. Finally, if local quantification and local balancing is developed, flows with turbulence transition areas could benefit from this method as well. 5. Concluding remarks Numerical analysis This work has presented an elaborate overview of the different methods in obtaining and avoiding checkerboard oscillations when using the projection method to solve incompressible flows on collocated grids, by using the symmetry-preserving framework of which the derivations and implementations of the discrete operators were also explained in detail. Distinction between usage of a wide-stencil or a compact-stencil Laplacian combined with a direct (DIM) or weighted (WIM) interpolation method lead to a total of four methods. The wide-stencil DIM is the method to which one would arrive when using a central differencing scheme, which is very prone to checkerboarding and therefore usually avoided. The compact-stencil DIM introduces an unnecessary dissipative smoothing Journal of Computational Physics 521 (2025) 113537 17 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. Fig. 12. (a) Mean stream-wise velocity on the 𝑋40 mesh. (b) Profiles of the viscous and Reynolds shear stress on the 𝑋40 mesh, showing that pisoFoam is not able to fully develop turbulence. (c) RMS velocities on the 𝑋40 mesh. Influence of the temporal integration scheme is shown by plotting the 𝜃0method using the RK3 scheme. of the velocity fields, and is therefore also undesirable. Between the WIMs, the wide-stencil WIM leads to an error in the convective term, which is not strictly-dissipative and can lead to unstable solutions. The pressure error of the compact-stencil WIM is favorable because its order can be reduced by employing the Van Kan method. Moreover, the compact-stencil Laplacian is computationally favorable. Quantifying checkerboarding The problem of the lack of definition of checkerboarding in literature was addressed and a solution to this problem was suggested by quantifying the phenomenon as 𝐶𝑐𝑏 =1 −𝐩𝑇 𝑐(𝐿𝑐)𝐩𝑐 𝐩𝑇 𝑐𝐿𝐩𝑐 , which is a physics-based, global, non-dimensional, normalised coefficient. The basis of this coefficient is tied to the difference between the wide-stencil and compact-stencil Laplacian operators. This difference is linked to the fundamental problem of the collocated grid arrangement, which is unavoidably accompanied by interpolations, of which there cannot be an inverse such that the original field is reestablished, i.e. Γ𝑐𝑠Γ𝑠𝑐 ≈𝐼𝑠. This difference also forms the basis of the numerical errors that are introduced by the different fractional step methods to counter-act checkerboard oscillations. Application of the checkerboard coefficient and performance This coefficient was tested using laminar and turbulent flows, and compared to qualitative and other quantitative assessments, in which it showed to give an intuitive and seemingly correct estimation of the levels of checkerboarding. An example of a possible usage of the coefficient was also introduced, by employing it as a parameter that gradually includes a pressure predictor into the momentum predictor equation. Since the inclusion of a pressure predictor itself can cause checkerboarding, this effectively establishes a dynamic balancing between checkerboarding and numerical dissipation through negative feedback. In the performed test cases, this solver showed its possibility in achieving this balance without requiring any user input. By doing so it was able to achieve low levels of numerical dissipation in cases without any oscillations, whereas it was also able to reduce the amount of checkerboarding in more Journal of Computational Physics 521 (2025) 113537 18 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. Fig. 13. (a) Turbulent kinetic energy production on the 𝑋40 mesh. (b) Turbulent kinetic energy transport on the 𝑋40 mesh. (c) Viscous diffusion on the 𝑋40 mesh. (d) Dissipation on the 𝑋40 mesh. (e) Pressure diffusion on the 𝑋40 mesh. Journal of Computational Physics 521 (2025) 113537 19 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. Fig. 14. (a) Mean stream-wise velocity on the 𝑋180 mesh. (b) RMS velocities on the 𝑋180 mesh. (c) Turbulent kinetic energy budgets on the 𝑋180 mesh. challenging cases. In the transient test case, the solver showed no significant difference in accuracy, while also showing convergence to the same solution as the Chorin and Van Kan methods. Scope of applicability Currently, the method has its greatest use for laminar flows, flows on coarse meshes, cases with transition between laminar and turbulent flows or more generally, simulations that are run without knowing the resulting flow a priori. For refined meshes and turbulent flows, the checkerboarding seems to diminish in general, and the operating window of the dynamic solver becomes smaller, while the accuracy of the results are similar to the Van Kan and Chorin methods, making the method less useful. Critical considerations One drawback of this method to consider is the loss of its second-order relation between time-step size and pressure error, for which it could be interesting to consider a different feedback mechanism instead of the proposed linear feedback, such that the method stops utilising a pressure predictor only at certain levels of checkerboarding. Another consideration of the quantification method is that the level of checkerboarding decreases greatly when the mesh is refined. Therefore, a single measurement of the level of checkerboarding without any context can be misleading or present little use. However, it still remains as a valid quantification method where no other methods are widely-known or -used. Future work An interesting suggestion for future works would be the measurement of checkerboarding on a local or even cell level. This could provide more insight into problematic areas, suggesting mesh refinement or even feedback mechanisms that act on a local level, for example by local inclusion of a pressure predictor or hybrid usage of a compactand wide-stencil scheme depending on the level of checkerboarding. This could then also make this method useful in cases that have laminar to turbulent flow transitional regions. Another potential area of interest is using 𝐶𝑐𝑏 in other feedback loops, such as for the determination of the time-step size, since it Journal of Computational Physics 521 (2025) 113537 20 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. affects the level of checkerboarding. However, since a decreasing time-step size increases the levels of checkerboarding, the negative feedback would tend to increase the time-step. Generally speaking, a larger time-step is already sought after, but is limited by other factors, such as the CFL condition. CRediT authorship contribution statement J.A. Hopman: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. D. Santos: Writing – review & editing, Methodology, Formal analysis. À. Alsalti-Baldellou: Writing – review & editing, Methodology, Formal analysis. J. Rigola: Supervision, Project administration, Funding acquisition. F.X. Trias: Writing – review & editing, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Appendix A. A note on collocated gradients Using the definition 𝐺𝛾 𝑐=Γ 𝛾 𝑠𝑐𝐺=Ω −1Γ𝛾𝑇 𝑐𝑠 Ω𝑠𝐺can be quite cumbersome in numerical analyses and in terms of implementation into codes. In this section this term is rewritten to a more intuitive notation which can be easily implemented into most codes. The starting point for this derivation is a general collocated gradient, 𝐺𝛾 𝑐, for which the interpolation method is not yet specified. Using the definitions from 2, the following derivation is made: 𝐺𝛾 𝑐=Γ 𝛾 𝑠𝑐𝐺 =Ω −1Γ𝛾𝑇 𝑐𝑠 Ω𝑠𝐺 =−Ω −1 [𝐼3⊗Π𝛾𝑇 𝑐𝑠 ]𝑁𝑇 𝑓𝑀𝑇 =Ω −1 [𝐼3⊗(𝑇𝑓𝑜𝑊𝛾 𝑜+𝑇𝑓𝑛𝑊𝛾 𝑛)]𝑆𝑇 𝑓(𝑇𝑇 𝑓𝑛 −𝑇𝑇 𝑓𝑜) =Ω −1 ⎛⎜⎜⎜⎜⎝ −[𝐼3⊗(𝑇𝑓𝑜 (𝐼𝑚−𝑊𝛾 𝑛))]𝑆𝑇 𝑓𝑇𝑇 𝑓𝑜 [𝐼3⊗(𝑇𝑓𝑜𝑊𝛾 𝑜)] 𝑆𝑇 𝑓𝑇𝑇 𝑓𝑛 −[𝐼3⊗(𝑇𝑓𝑛𝑊𝛾 𝑛)] 𝑆𝑇 𝑓𝑇𝑇 𝑓𝑜 [𝐼3⊗(𝑇𝑓𝑛 (𝐼𝑚−𝑊𝛾 𝑜))]𝑆𝑇 𝑓𝑇𝑇 𝑓𝑛 ⎞⎟⎟⎟⎟⎠ =Ω −1 [𝐼3⊗(𝑇𝑓𝑜 −𝑇𝑓𝑛)]𝑆𝑇 𝑓(𝑊𝛾 𝑛𝑇𝑇 𝑓𝑜 +𝑊𝛾 𝑜𝑇𝑛𝑓 ) =Ω −1 [𝐼3⊗𝑀]𝑁𝑇 𝑓(𝑊𝛾 𝑛𝑇𝑇 𝑓𝑜 +𝑊𝛾 𝑜𝑇𝑛𝑓 ) =𝐺𝐺Π𝛾 𝑐𝑠, (A.1) where at each cell 𝑖the following simplification is made: [−[𝐼3⊗𝑇 𝑓𝑜𝐼𝑚]𝑆𝑇 𝑓𝑇𝑇 𝑓𝑜 +[𝐼3⊗𝑇 𝑓𝑛𝐼𝑚]𝑆𝑇 𝑓𝑇𝑇 𝑓𝑛]𝑖=− ∑ 𝑓∈𝐹(𝑖) 𝐬𝑓(𝑖)=𝟎,(A.2) which follows from the fact that the sum of all outward-pointing surface vectors of a closed geometry equals zero, which holds for any number of dimensions. In the final line, two new operators are introduced, the Gauss gradient operator, 𝐺𝐺, and the inverse-weighted cell-to-face interpolator, Π𝛾 𝑐𝑠. The Gauss gradient can be explained as assigning the normal direction of each face to the scalar value at the face multiplied by the face area, summing these face-vectors to the cell-center, then dividing by the cell volume, which is a common method to take gradients in the finite volume method. The inverse-weighted cell-to-face interpolator operates exactly as the normal cell-to-face interpolator, but the weights of neighbour and owner are swapped. This has the noteworthy property that: Π𝑀 𝑐𝑠 = Π𝑀 𝑐𝑠 ,Π𝑉 𝑐𝑠 = Π𝐿 𝑐𝑠,Π𝐿 𝑐𝑠 = Π𝑉 𝑐𝑠,(A.3) such that finally: 𝐺𝑐=Γ 𝑉 𝑠𝑐𝐺=𝐺𝐺Π𝐿 𝑐𝑠.(A.4) In conclusion, using the derivation provided by equation (A.1), in combination with the definitions of equation (A.3), it becomes easy to rewrite any wide-stencil gradient to a combination of a scalar interpolation and a Gauss gradient, which can be helpful in numerical analysis. Moreover, as discussed below, it greatly simplifies implementation into numerical algorithms. Journal of Computational Physics 521 (2025) 113537 21 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. A.1. A note on implementation Aside from simplifying numerical analyses, equation (A.4)greatly simplifies the implementation of the collocated gradient into code. For example, in the popular open-source finite volume code OpenFOAM, implementing 𝐺𝑐as Γ𝑠𝑐𝐺is not straight-forward, especially in symmetry-preserving codes. For these codes, the choice of interpolation has to be hard-coded into the solver, to remove any degree of freedom in choice of interpolator. Instead of writing a function for the face-to-cell interpolator, Γ𝑠𝑐 , equation (A.4) can be used, resulting in more conventional functions that are usually readily available. For example, the symmetry-preserving OpenFOAM solver RKSymFoam uses the following syntax to hard-code 𝐺𝑉 𝑐𝐩𝑛 𝑐as 𝐺𝐺Π𝐿 𝑐𝑠𝐩𝑛 𝑐[48]: const volVectorField gradpn ( fv :: gaussGrad<scalar >::gradf ( linear<scalar>(mesh). interpolate(pn), "gradpn" ) ); which uses the gradf() and interpolate() functions that are provided in the source code of OpenFOAM. Appendix B. Kernel vectors for arbitrary Cartesian meshes In this section, a set of kernel vectors for arbitrary Cartesian meshes is derived for the linear, midpoint and volumetric wide-stencil Laplacian operators, 𝐿𝐿 𝑐, 𝐿𝑀 𝑐, 𝐿𝑉 𝑐. The derivation is shown for the two-dimensional case. In this case, four linearly independent kernel vectors are required to span the null-space; the constant kernel vector, 𝐩− 𝑐(00,𝛾), the vertical kernel vector, 𝐩− 𝑐(10,𝛾), the horizontal kernel vector, 𝐩− 𝑐(01,𝛾)and finally, the checkered kernel vector, 𝐩− 𝑐(11,𝛾). Since the case for 𝛾=𝑀was given in [33], only the cases for 𝛾=𝐿 and 𝛾=𝑉will be derived here. For the derivations, the fact that 𝐾𝑒𝑟(𝐺𝛾 𝑐) ∈𝐾𝑒𝑟(𝐿𝛾 𝑐)was used in addition to the equality 𝐺𝛾 𝑐=𝐺𝐺Π𝛾 𝑐𝑠, which was derived in equation (A.1). After applying these equalities, a set of vectors, 𝜙 𝜙 𝜙𝑐needs to be found such that: 𝐺𝛾 𝑐𝜙 𝜙 𝜙𝑐=𝐺𝐺Π𝛾 𝑐𝑠𝜙 𝜙 𝜙𝑐=𝟎𝑐.(B.1) Since the constant kernel vector is a trivial solution to this equation, the kernel vector that is considered first is the vertical linear kernel vector, 𝐩− 𝑐(10,𝐿). If the values of this vector are chosen as: [𝐩− 𝑐(10,𝐿)]𝑖,𝑗 =(−1) 𝑖(Δ𝑥𝑖)−1 ,(B.2) a vertically striped pattern will be obtained, as seen in Fig. B.15. Since the vertically neighbouring cells have equal values, the values at faces 𝑠and 𝑛are trivially given by: [Π𝐿 𝑐𝑠𝐩− 𝑐(10,𝐿)]𝑠=[Π𝐿 𝑐𝑠𝐩− 𝑐(10,𝐿)]𝑛=−1 Δ𝑥1 .(B.3) Given the fact that Π𝐿 𝑐𝑠 =Π 𝑉 𝑐𝑠, the value at face 𝑤is calculated as: [Π𝐿 𝑐𝑠𝐩− 𝑐(10,𝐿)]𝑤=[Π𝑉 𝑐𝑠𝐩− 𝑐(10,𝐿)]𝑤=Δ𝑥0𝜙𝑊+Δ𝑥1𝜙𝑃 Δ𝑥0+Δ𝑥1 =Δ𝑥0∕Δ𝑥0−Δ𝑥1∕Δ𝑥1 Δ𝑥0+Δ𝑥1 =0. (B.4) Similarly, [Π𝐿 𝑐𝑠𝐩− 𝑐(10,𝐿)]𝑒=0. Since this leads to opposing face pairs with equal values, it becomes immediately evident that [𝐺𝐺Π𝐿 𝑐𝑠𝐩− 𝑐(10,𝐿)]𝑃=0. This equality holds throughout the whole mesh and therefore the vector given by equation (B.2)lies on the kernel of 𝐿𝐿 𝑐. The horizontal linear kernel vector is derived similarly by swapping the axes. This will be shown by deriving the values for the horizontal volumetric kernel vector, 𝐩− 𝑐(01,𝑉 ), in which, additionally, the values at each cell are replaced by their inverse. This kernel vector has its values given by: [𝐩− 𝑐(01,𝑉 )]𝑖,𝑗 =(−1) 𝑗Δ𝑦𝑗,(B.5) so that: Journal of Computational Physics 521 (2025) 113537 22 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. Fig. B.15. Vertical linear kernel vector, 𝐩− 𝑐(10,𝑉 ), with 𝐼=0,𝐽=0and 𝛾=𝑉. [Π𝑉 𝑐𝑠𝐩− 𝑐(01,𝑉 )]𝑤=[Π𝑉 𝑐𝑠𝐩− 𝑐(01,𝑉 )]𝑒=−Δ𝑦1,(B.6) [Π𝑉 𝑐𝑠𝐩− 𝑐(01,𝑉 )𝑐]𝑠=[Π𝐿 𝑐𝑠𝐩− 𝑐(01,𝑉 )]𝑠=Δ𝑦1𝜙𝑆+Δ𝑦0𝜙𝑃 Δ𝑦0+Δ𝑦1 (B.7) =Δ𝑦1Δ𝑦0−Δ𝑦0Δ𝑦1 Δ𝑦0+Δ𝑦1 =0, [Π𝑉 𝑐𝑠𝐩− 𝑐(01,𝑉 )]𝑛=0,(B.8) which also leads to equal opposite face pair values. The checkered kernel vectors are slightly more difficult to see directly. For 𝐿𝑉 𝑐, the values of kernel vector 𝐩− 𝑐(11,𝑉 )are given by: [𝐩− 𝑐(11,𝑉 )]𝑖,𝑗 =(−1) 𝑖+𝑗Δ𝑥𝑖Δ𝑦𝑗,(B.9) leading to: [Π𝑉 𝑐𝑠𝜙 𝜙 𝜙𝑐]𝑤=[Π𝐿 𝑐𝑠𝜙 𝜙 𝜙𝑐]𝑤=Δ𝑥1𝜙𝑊+Δ𝑥0𝜙𝑃 Δ𝑥0+Δ𝑥1 (B.10) =Δ𝑥1(Δ𝑥0Δ𝑦1)−Δ𝑥0(Δ𝑥1Δ𝑦1) Δ𝑥0+Δ𝑥1 =0, [Π𝑉 𝑐𝑠𝜙 𝜙 𝜙𝑐]𝑒=[Π𝑉 𝑐𝑠𝜙 𝜙 𝜙𝑐]𝑠=[Π𝑉 𝑐𝑠𝜙 𝜙 𝜙𝑐]𝑛=0.(B.11) Which gives similar results for the checkered volumetric kernel vector. Extending these arguments for all combinations of patterns and interpolators and to the three-dimensional case, a full set of linearly independent kernel vectors for arbitrary Cartesian meshes can be derived. For the wide-stencil Laplacian operators 𝐿𝛾 𝑐this set is given by: [𝐩− 𝑐(𝐼𝐽𝐾,𝛾)]𝑖,𝑗,𝑘 =(−1) 𝑖𝐼+𝑗𝐽+𝑘𝐾 ([Δ𝑥]𝐼 𝑖[Δ𝑦]𝐽 𝑗[Δ𝑧]𝐾 𝑘)𝛼 ,∀𝐼,𝐽,𝐾 ∈{0,1}, 𝛼=⎧ ⎪ ⎨ ⎪ ⎩ −1 if 𝛾=𝐿 0if 𝛾=𝑀 1if 𝛾=𝑉 . (B.12) In conclusion, a set of vectors that spans the kernel of 𝐿𝛾 𝑐can always be derived for any Cartesian mesh and for any of the discussed interpolators, using equation (B.12). Since the kernel is sometimes used in describing the checkerboard phenomenon, using equation (B.12)saves a lot of computational effort compared to performing a singular value decomposition. Data availability Data will be made available on request. Journal of Computational Physics 521 (2025) 113537 23 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. References [1] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Company, 1980. [2] P. Wesseling, Principles of Computational Fluid Dynamics, vol. 29, Springer Science & Business Media, 2009. [3] J.H. Ferziger, M. Perić, Computational Methods for Fluid Dynamics, Springer-Verlag, Berlin, Heidelberg, New York, 2002. [4] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Pearson Education, 2007. [5] A.J. Chorin, Numerical solution of the Navier–Stokes equations, Math. Comput. 22 (104) (1968) 745–762. [6] R. Temam, Sur l’approximation de la solution des equations de Navier—Stokes par la methode des pas fractionnaires ii, Arch. Ration. Mech. Anal. 33 (1969) 377–385. [7] S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows, in: Numerical Prediction of Flow, Heat Transfer, Turbulence and Combustion, Elsevier, 1983, pp. 54–73. [8] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier–Stokes equations, J. Comput. Phys. 59 (2) (1985) 308–323. [9] J. Van Kan, A second-order accurate pressure-correction scheme for viscous incompressible flow, SIAM J. Sci. Stat. Comput. 7(3) (1986) 870–891. [10] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (12) (1965) 2182–2189. [11] B. Perot, Conservation properties of unstructured staggered mesh schemes, J. Comput. Phys. 159 (1) (2000) 58–89. [12] A. Pascau, Cell face velocity alternatives in a structured colocated grid for the unsteady Navier–Stokes equations, Int. J. Numer. Methods Fluids 65 (7) (2011) 812–833. [13] C.M. Rhie, W.L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21 (11) (1983) 1525–1532. [14] M. Perić, R. Kessler, G. Scheuerer, Comparison of finite-volume numerical methods with staggered and colocated grids, Comput. Fluids 16 (4) (1988) 389–403. [15] W.Z. Shen, J.A. Michelsen, J.N. Sørensen, Improved Rhie-Chow interpolation for unsteady flow computations, AIAA J. 39 (12) (2001) 2406–2409. [16] S.K. Choi, Note on the use of momentum interpolation method for unsteady flows, Numer. Heat Transf., Part A, Appl. 36 (5) (1999) 545–550. [17] B. Yu, Y. Kawaguchi, W.Q. Tao, H. Ozoe, Checkerboard pressure predictions due to the underrelaxation factor and time step size for a nonstaggered grid with momentum interpolation method, Numer. Heat Transf., Part B, Fundam. 41 (1) (2002) 85–94. [18] B. Yu, W.-Q. Tao, J.-J. Wei, Y. Kawaguchi, T. Tagawa, H. Ozoe, Discussion on momentum interpolation method for collocated grids of incompressible flow, Numer. Heat Transf., Part B, Fundam. 42 (2) (2002) 141–166. [19] S. Majumdar, Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids, Numer. Heat Transf. 13 (1) (1988) 125–132. [20] T.F. Miller, F.W. Schmidt, Use of a pressure-weighted interpolation method for the solution of the incompressible Navier–Stokes equations on a nonstaggered grid system, Numer. Heat Transf. 14 (2) (1988) 213–233. [21] B. Yu, W.Q. Tao, J.J. Wei, Y. Kawaguchi, T. Tagawa, H. Ozoe, Discussion on momentum interpolation method for collocated grids of incompressible flow, Numer. Heat Transf., Part B, Fundam. 42 (2) (2002) 141–166. [22] J.L. Guermond, P. Minev, J. Shen, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Eng. 195 (44–47) (2006) 6011–6045. [23] S. Zhang, X. Zhao, S. Bayyuk, Generalized formulations for the Rhie–Chow interpolation, J. Comput. Phys. 258 (2014) 880–914. [24] P. Bartholomew, F. Denner, M.H. Abdol-azis, A. Marquis, B.G.M. Van Wachem, Unified formulation of the momentum-weighted interpolation for collocated variable arrangements, J. Comput. Phys. 375 (2018) 177–208. [25] F.N. Felten, T.S. Lund, Kinetic energy conservation issues associated with the collocated mesh scheme for incompressible flow, J. Comput. Phys. 215 (2) (2006) 465–484. [26] J.E. Mattson, An Introduction to Ansys Fluent 2023, 1st edition, SDC Publications, 2023. [27] S.D.I. Software, Simcenter STAR-CCM+ User Guide v. 2306, 2023. [28] C. Greenshields, OpenFOAM v11 User Guide, London, UK, https://doc .cfd .direct /openfoam /user -guide -v11, 2023. [29] F. Archambeau, N. Méchitoua, M. Sakiz, Code saturne: a finite volume code for the computation of turbulent incompressible flows-industrial applications, Int. J. Finite Vol. 1(1) (2004). [30] R.W. Verstappen, A.E. Veldman, Symmetry-preserving discretization of turbulent flow, J. Comput. Phys. 187 (1) (2003) 343–368. [31] F.X. Trias, O. Lehmkuhl, A. Oliva, C.D. Pérez-Segarra, R.W. Verstappen, Symmetry-preserving discretization of Navier–Stokes equations on collocated unstructured grids, J. Comput. Phys. 258 (2014) 246–267. [32] E.M. Komen, J.A. Hopman, E.M. Frederix, F.X. Trias, R.W. Verstappen, A symmetry-preserving second-order time-accurate piso-based method, Comput. Fluids 225 (2021) 104979. [33] J. Larsson, G. Iaccarino, et al., A co-located incompressible Navier–Stokes solver with exact mass, momentum and kinetic energy conservation in the inviscid limit, J. Comput. Phys. 229 (12) (2010) 4425–4430. [34] J.A. Hopman, F. Trias, J. Rigola, On a conservative solution to checkerboarding: examining the discrete laplacian kernel using mesh connectivity, in: ERCOFTAC Workshop Direct and Large Eddy Simulation, Springer, 2023, pp. 306–311. [35] J.A. Hopman, F.X. Trias Miquel, J. Rigola Serrano, Symmetry-preserving discretisation methods for magnetohydrodynamics, in: World Congress in Computational Mechanics and ECCOMAS Congress, Oslo, Norway, 2022. [36] D. Santos, F.X. Trias, G. Colomer, C.D. Pérez-Segarra, An energy-preserving unconditionally stable fractional step method on collocated grids, in: World Congress in Computational Mechanics and ECCOMAS Congress, Oslo, Norway, 2022. [37] D. Santos, F.X. Trias, J.A. Hopman, C.D. Pérez-Segarra, Pressure-velocity coupling on unstructured collocated grids: reconciling stability and energy-conservation, in: Proceedings of the Tenth International Symposium on Turbulence, Heat and Mass Transfer, Rome, Italy, 2023, pp. 259–262. [38] D. Santos, J.A. Hopman, C.D. Pérez-Segarra, F.X. Trias, On an energy-preserving unconditionally stable projection method on collocated unstructured grids, Available at SSRN: https://doi .org /10 .2139 /ssrn .4788673, 2024. [39] N.N. Yanenko, The Method of Fractional Steps, Springer-Verlag, New York, 1971. [40] A.J. Chorin, A numerical method for solving incompressible viscous flow problems, J. Comput. Phys. 2(1) (1967) 12–26. [41] J.A. Hopman, A. Alsalti-Baldellou, F.X. Trias, J. Rigola, On a conservative solution to checkerboarding: examining the causes of non-physical pressure modes, in: Proceedings of the 14th International ERCOFTAC Symposium on Engineering Turbulence Modelling and Measurements, Barcelona, Spain, 2023, pp. 599–603. [42] A.E.P. Veldman, A general condition for kinetic-energy preserving discretization of flow transport equations, J. Comput. Phys. 398 (2019) 108894. [43] M.M. Larmaei, J. Behzadi, T.F. Mahdi, Treatment of checkerboard pressure in the collocated unstructured finite-volume scheme, Numer. Heat Transf., Part B, Fundam. 58 (2) (2010) 121–144. [44] C.M. Klaij, On the stabilization of finite volume methods with co-located variables for incompressible flow, J. Comput. Phys. 297 (2015) 84–89. [45] P. Rauwoens, J. Vierendeels, B. Merci, A solution for the odd-even decoupling problem in pressure-correction algorithms for variable density flows, J. Comput. Phys. 227 (1) (2007) 79–99. [46] A.W. Date, Fluid dynamical view of pressure checkerboarding problem and smoothing pressure correction on meshes with colocated variables, Int. J. Heat Mass Transf. 46 (25) (2003) 4885–4898. [47] G.H. Golub, C.F. Van Loan, Matrix Computations, Johns Hopkins University Press, 1996. [48] J.A. Hopman, E.M.A. Frederix, RKSymFoam GitHub page, https://github .com /janneshopman /RKSymFoam, 2023. [49] V. Vuorinen, J.-P. Keskinen, C. Duwig, B.J. Boersma, On the implementation of low-dissipative Runge–Kutta projection methods for time dependent flows using OpenFOAM®, Comput. Fluids 93 (2014) 153–163. Journal of Computational Physics 521 (2025) 113537 24 J.A. Hopman, D. Santos, À. Alsalti-Baldellou et al. [50] G.I. Taylor, A.E. Green, Mechanism of the production of small eddies from large ones, Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci. 158 (895) (1937) 499–521. [51] A.W. Vreman, J.G. Kuerten, Comparison of direct numerical simulation databases of turbulent channel flow at 𝑅𝑒𝜏= 180, Phys. Fluids 26 (1) (2014). [52] E.M.J. Komen, E.M.A. Frederix, T.H.J. Coppen, V. D’alessandro, J.G.M. Kuerten, Analysis of the numerical dissipation rate of different Runge-Kutta and velocity interpolation methods in an unstructured collocated finite volume method in openfoam, Comput. Phys. Comun. (2020). [53] H. Zhang, F.X. Trias, A. Gorobets, Y. Tan, A. Oliva, Direct numerical simulation of a fully developed turbulent square duct flow up to 𝑅𝑒𝜏= 1200, Int. J. Heat Fluid Flow 54 (2015) 258–267. [54] P. Durbin, B. Reif, Statistical Theory and Modeling for Turbulent Flows, Wiley & Sons, Ltd, 2011. [55] J.A. Hopman, runTimeChannelBudgets GitHub page, 2023, https://github .com /janneshopman /runTimeChannelBudgets. [56] S.B. Pope, Turbulent Flows, Cambridge University Press, 2000.