{'fraction_non_alphanumeric': 0.04365966323945761, 'fraction_numerical': 0.010480305965330552, 'mean_word_length': 4.023453093812376, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
ON THE INTERPOLATION PROBLEM FOR THE POISSON EQUATION ON COLLOCATED MESHES

The appearence of unphysical velocities in highly distorted meshes is a common problem in many simulations. In collocated meshes, this problem arises from the interpolation of the pressure gradient from faces to cells. Using an algebraic form for the classical incompressible Navier-Stokes equations, this problem is adressed. Starting from the work of F. X. Trias et. al. [FX.Trias et al. JCP 258: 246-267, 2014], a new approach for studying the Poisson equation obtained using the Fractional Step Method is found, such as a new interpolator is proposed in order to found a stable solution, which avoid the appearence of these unpleasant velocities. The stability provided by the interpolator is formally proved for cartesian meshes and its rotations, using fully-explicit time discretizations. The construction of the Poisson equation is supported on mimicking the symmetry properties of the differential operators and the Fractional Step Method. Then it is reinterpreted using a recursive application of the Fractional Step Method in order to study the system as an stationary iterative solver. Furthermore, a numerical analysis for unstructured mesh is also provided.


INTRODUCTION
The continuity and Navier-Stokes equations are the equations which describe the mass and momentum conservation. Even though they were formulated long time ago, they still do not have an analytical solution for general cases. This is why numerical methods are so important, because they give us numerical approximations of these solutions. One family of numerical methods that can be applied to the resolution of this equations are the Finite Volume Methods. These methods are characterised by their apparent mathematical simplicity and intuitive description of the fundamental physics of the equations. Despite the fact that a staggered configuration would be the natural way in which the FVM should be applied to discretise the equations [1,2,3], a collocated grid arrangement is often preferred due to its simplicity in order to deal with complex geometries and its computational cost [4,5]. Indeed, this approach is found in the most popular codes, such as ANSYS-FLUENT or OpenFOAM.
One method used to solve the velocity-pressure coupling is the classical fractional step method [6,7,8].
However, working with collocated meshes and a wide stencil Poisson equation in such a method lead us to the well-known checkerboard problem. This problem can be solved by means of working with a compact stencil Poisson equation [9,10] or by regularizing of the convective term such as [4].
Furthermore, other problems or unstabilities can be found in the numerical solution of these equations. One of these problems is the appearence of spurious velocities in highly distorted meshes, coming from the interpolation from faces to cells of the discrete pressure gradient, which is strongly related to the resolution of the Poisson equation itself. Many reconstructions of the discrete gradient have been proposed [11,12,13,14]. Some of them were proposed in terms of improving accuracy, but they can also fix the interpolation problem. However, they do not usually respect the underlying symmetries of the continuos gradient and divergence operator. This work is focused on this problem from the point of view of respecting these symmetries.
2 Symmetry-preserving finite volume discretization on collocated grids The well-known dimensionless Navier-Stokes equations for the simulation of laminar or turbulent, Newtonian and incompressible flows in primitive variables are where Re is the dimensionless Reynolds number. A fully-conservative finite-volume discretization respecting the symmetries of the differential operators for structured and unstructured meshes was presented in [4]. Using this discretization approach, the properties of the continuous equations are preserved, as far as the symmetries of the convective and diffusive operators.
As a collocated formulation is preferred over a staggered one due to its simplicity, pressure and velocity are computed at the control volume center. Following the notation used in [4], let us define p c = (p 1 , p 2 , . . . , p n ) * ∈ R n and u c ∈ R 3n as the cell-centered pressure and velocity fields, where n is the number of control volumes. The subindices c and s refer to whether the variables are cell-centered or staggered at the faces. For simplicity u c is defined as a column vector and arranged as u c = (u 1 , u 2 , u 3 ) * , are the vectors containing the velocity components corresponding to the x i -spatial direction. In order to assure mass conservation in the control volumes, a secondary velocity field is defined at the faces u s = (u s ) 1 , (u s ) 2 , (u s ) 3 , . . . , (u s ) m * ∈ R m where m is the number of faces on the computational domain. Defining a linear shift transformation from cell to faces Γ c→s ∈ R m×3n , both velocities can be related as follows The matrices Ω ∈ R 3n×3n , C (u s ) ∈ R 3n×3n and D ∈ R 3n×3n are block diagonal matrices given by where I 3 ∈ R 3×3 is the identity matrix and Ω c ∈ R n×n is a diagonal matrix with the cell-centered control volumes. C c (u s ) ∈ R n×n and D c ∈ R n×n are the cell-centered convective and diffusive operators for a discrete scalar field, respectively. Finally, G c ∈ R 3n×n represents the discrete gradient operator and the matrix M ∈ R n×m is the face-to-center discrete divergence operator.
In this way, the equations on an arbitrary collocated mesh can be written as follows The construction of the discrete operators is adressed in section 3 of [4] and it is not the object of this work.

Solving the pressure-velocity coupling. Checkerboard problem
In order to simplify the notation, semidiscretized momentum equation (4a) can be rewritten as follows where R (u c ) ≡ −Ω −1 (C (u s ) u c + Du c ) represents the convective and diffusive terms. For the temporal discretization, a fully explicit scheme is assumed whereas the incompressibility constraint is treated implicitly. Introducing the functions f and R, which represent a general fully explicit discretization, the fully-discretized NS equations read In the case presented in [4], Observe that, in general, we can split f into two parts. The first one will contain the new velocity u n+1 c and the second one will contain all the previous ones. This split can be represented as where α is a general parameter and g is a function containing the contribution of all previous velocites.
To solve the velocity-pressure coupling, a classical fractional step method [6,7,8] is used. For the staggered velocity field, u s , this projection is 'naturally' derived from the Helmholtz-Hodge vector decomposition theorem [16], whereby a velocity u p s can be uniquely decomposed into a solenoidal vector, u n+1 s , and a curl-free vector, expressed as the gradient of a scalar field, Gp c . This decomposition is written as where G ∈ R m×n is the center-to-face staggered gradient operator which is related with the divergence operator via [4] G where Ω s ∈ R m×m is a diagonal matrix with the staggered control volumes. Then, taking the divergence of Eq.(7) yields a discrete Poisson equation forp c Finally, using the definition of G given in Eq.(8) the previous equation becomes where the discrete Laplacian operator L ∈ R n×n is, by construction, a symmetric negative-definite matrix.
Once the solution is obtained, u n+1 s results from the correction (7) u n+1

Constraints on the shift operators
Linear shift operators are needed to relate the cell-centered velocity fields to the staggered ones and viceversa. Namely, the linear shift operator, Γ c→s ∈ R m×3n , transforms a cell-centered velocity field into a staggered one whereas the cell-centered fields are related to the staggered ones via the linear shift transformation Γ s→c ∈ R 3n×m , u c = Γ s→c u s .
Notice, that in general Γ s→c Γ c→s = I holds only approximately i.e. u c ≈ Γ s→c Γ c→s u c . More importantly, recalling the definition for the staggered gradient operator given in Eq. (8), the cell-centered discrete gradient operator results and, therefore, the face-to-cell shift operator Γ s→c is restricted as follows to force that the pressure gradient contribution to the global kinetic energy exactly vanishes (see section 2.2 of [4]).
3 Analysis on the Poisson equation obtained using the FSM 3.1 Computing the cell-center predictor velocity u p c Let us introduce the following notation: In the projection method, assuming we are at the n + 1 time step, the predicted velocity, u p c , is computed as follows wherep c p is a 'prediction' of the pressure at the next time step,p c n+1 =p c p +p c . The definition of p c introduce an additional error to the momentum equation in a collocated formulation. Furthermore, as we are working with collocated meshes, we will face the checkerboard problem. This problem is widely adressed in sections 2.3.3 and 2.4 of [4]. So equation (17) can be rewritten as: Introducing the last result in equation (10), we obtain: where L c = M c G c ∈ R n×n is the collocated Laplacian operator. Finally, adding Lp c p at both sides, we arrive at: 3.2 Reinterpretation of the recursive application of the FSM as a stationary iterative solver of the wide-stencil Poisson equation

Particular case: a first-order forward Euler
For the first-order forward Euler scheme, we can specify the previous functions and constants as α = 1, h n = −u n c and R n = R (u n c ). Assuming we are in n + 1 time step, the velocity centered field at time n can be computed as follows: Introducing (21) in (20) leads to: Now, we can definep c n = ∑ n i=0pc i . Observe thatp c n+1 =p c n+1 −p c n . This allows us to rearrange the previous equation as: This can be viewed as a stationary iterative solver. The stability of this process will depend on the eigenvalues of L −1 L c and (L − L c ). This is going to be analysed deeply for the general case in the next section.
3.2.2 The general case: assuming a general fully-explicit scheme Recalling we are working with collocated meshes and assuming we are in n + 1 time step, the velocity centered field at time n can be computed as follows: Recalling that h is a function which contains a linear combination of the previous velocities, we can apply this result recursively up to u 0 c and rearrange the results in new functions: where K is a constant, and F and H are functions that depend onp c i and R i respectively, and also in the election of the time step ∆t at each iteration.
Observe that, in the same way we procedeed for (25), we can rewrite equation (20) in terms of all the previous time steps using the last result as follows: Lp c n+1 =K + L cF (p c n ,p c n−1 , ...,p c 0 ) +H(∆t n+1 R n , ∆t n R n−1 , ..., Here, we have done an arrangement in terms of functions as we did in (25). This can be viewed as a stationary iterative solver. Note that inside R i we also have the velocities at previous time steps and, consequently, the pressures at previous time steps, but, if the time step ∆t is chosen wisely, the stability of this process will depend on the eigenvalues of L −1 L c and (L − L c ), which subsequently depend on the interpolator operators. This will be carefully analysed in the next sections.

Conditions to solve the system
Inverting L in (26), and taking into account thatK and the function H can be bounded by means of a wisely chosen time step, then, let us define T 1 = L −1 (L − L c ) and T 2 = L −1 L c . If µ T 1 i and µ T 2 i are their eigenvalues, the process will converge if |µ T 1 i | < 1 and |µ T 2 i | < 1. Notice that One the one hand, L −1 L c is the product of two negative definite matrices (see Corollary 1), then, by Corollary 2, L −1 L c is a positive definite matrix. So µ T 1 i = 1 − µ T 2 i < 1 (it is a well-known result in solving iterative systems). On the other hand, if L − L c is negative definite, then L −1 (L − L c ) is positive definite, so µ T 1 i > 0. We can conclude that µ T 1 i ∈ (0, 1) and it implies that µ T 2 i ∈ (0, 1) and the process converge. In conclusion, we have to study the definiteness of L − L c .

Relation of L − L c and the interpolator operators
We have to proof when L − L c is negative definite: According to Lemma 1, the defineteness of −M(Ω −1 s − Γ c→s Ω −1 Γ * c→s )M * is given by the defineteness of Ω −1 s − Γ c→s Ω −1 Γ * c→s , so, as we want L − L c to be negative definite, we should impose Ω −1 s − Γ c→s Ω −1 Γ * c→s to be positive definite. Using the shift operators relation (15): Multiplying each side of the expression by Ω s and applying again Lemma 1, we conclude that should be positive definite. As we can see here, the defineteness of L − L c are directly related to the interpolation from faces to cells matrix. Ideally, as it is easier to define the interpolation from cells to faces via Γ c→s , we would like some expression involving this matrix. In order to do that, using (15) and multiplying each side of the expression by Ω s , we will need ρ(B) ≤ 1, being ρ(·) the spectral radius (observe that B is positive definite).
s , we can construct the previous matrix B = D T D and another matrix C = DD T . A general theorem states that eigen(D T D) = eigen(DD T ) except for λ = 0 (and its multiplicity) eigenvalue. So, So, we need I − Ω − 1 2 Γ T c→s Ω s Γ c→s Ω − 1 2 to be positive definite. Finally, we can multiply the previous expression by Ω 1 2 and we will obtain that Finally, we have proven that L − L c being negative definite is equivalent to Ω − Γ T c→s Ω s Γ c→s being positive definite. An interpretation of this in terms of a Cartesian mesh can be found in the next section.

Relation between the interpolator matrix and the defineteness
Taking into account that Γ c→s = N s Π = N s,1 N s,2 N s,3 (I 3 ⊗ Π c→s ), where the matrices N s,i ∈ R m×m are diagonal matries containing the x i -spatial components of the face normal vectors, and Π c→s ∈ R m×n is the operator that interpolates a cell-centered scalar field to the faces. The calculation of the matrix and the coefficients are straightforward: c→s Ω s Γ c→s =   Π T c→s N s,1 Ω s N s,1 Π c→s Π T c→s N s,1 Ω s N s,2 Π c→s Π T c→s N s,1 Ω s N s,3 Π c→s Π T c→s N s,2 Ω s N s,1 Π c→s Π T c→s N s,2 Ω s N s,2 Π c→s Π T c→s N s,2 Ω s N s,3 Π c→s Π T c→s N s,3 Ω s N s,1 Π c→s Π T c→s N s,3 Ω s N s,2 Π c→s Π T c→s N s,3 Ω s N s,3 Π c→s and where a i j are the coefficients of Π c→s . For more details in the notation, the reader is referred again to [4]. Now, in order to relate all this to the grid, let us recall the Greshgorin circle theorem. It tells us that, if A ∈ R n×n , R i = ∑ j =i |a i j |, and D(a ii , R) is the closed disc of radius R and center a ii , then every eigenvalue of A lies within at least one D(a ii , R). Given a general scalar interpolation matrix from cells to faces, and applying the Greshgorin circle theorem to our previous matrix Ω − Γ T c→s Ω s Γ c→s , we will find that, if V r − ∑ n s=1 ∑ j=x,y,z | ∑ m l=1 a lr V stagg l a ls n i l n j l | ≥ 0 ∀r ∈ {1, ..., n}, i ∈ {x, y, z}, then Ω − Γ T c→s Ω s Γ c→s is positive definite and consequently, L − L c is negative definite. We will understand a Cartesian mesh as a mesh whose faces are parallel to Cartesian axis. For our purpose, we can have elements with different volumes. Let us assume we apply the previous result to this kind of mesh, in particular to the highlighted control volume. If we take i = x:

A stable interpolation for Cartesian grids
Here, we are changing the numbers for the cardinal points N, S,W, E for simplicity, s is the subindex referred to control volumes and l is the subindex referred to faces. For example, a eP is the interpolator coefficient from cell P to face e. The faces n, s, e, w are referred to the control volume P.
Now, let us define the following interpolator, being i the index used for faces and j the index used for control volumes: whereṼ i j is the quantity of volume of V stagg i which is inside the control volume j. This quantity can be computed asṼ i j = d i j * S i , where d i j is the distance from control volume center j to the face i, and S i is the surface of i. The interpolator is called volume weighted. Observe that, if j and k are two adjacent control volumes which share face i, then a i j + a ik = 1. Applying this interpolator to equation (37) we obtain: Last equality holds becauseṼ eP −Ṽ wP = d eP * S e + d wP * S w = V P , independently of the cell center location (at least if we define d i j as the distance from the cell node to the plane which contains the face). However, this suggests putting the cell center right in the center of the rectangle, which is the common way. Furthermore, if we have a non-periodic case, and we are in the adjacent cell of the (for example) left boundary, observe that we would get V P −Ṽ eP ≥ 0. So, the interpolator gives a stable solution for the Poisson equation in Cartesian grids even in non-periodic cases.
Observe that the results for the y-coordinate are analogous to the previous ones. Furthermore, due to the same symmetries, this result can be extended to a 3-dimensional case (or to a general n-dimensional case). The extension is straightforward. In addition, the stability can be also proved for rotations of the Cartesian meshes, but the proof is not presented here due to its extension.

Comments on unstructured grids
As no theoretical result could be found for triangular unstructured meshes, a numerical test has been carried out: • Checking the sign of the eigenvalues of L − L c .
• Computing the global error and the rate of convergence of the new interpolator.
When placing the cell center in the orthocenter, some eigenvalues of L − L c lie on the positive axis, even using the weighted volume scheme. This problem seems to be avoided if the cell center is placed in the circumcenter and it is placed inside the triangle. However, when the circumcenter lies outside the triangle the previous problem appears again. The possibility of working with negative staggered volumes has been tested, but it introduces other problems such as L not being negative definite. A possible solution would be to work with meshes constructed with acute-angled triangles, because they have the circumcenter inside the triangle.
For computing the global error and the rate of convergence, • 60 random meshes have been used.
• Average step-size is defined as trace(Ω c )/n, where n is the number of control volumes.
• The next mapping function has been used to stretch the meshes:  The global error is understood as ||p s c − p c c || 2 , where: • p s c is obtained solving Lp s c = Mu s . • p c c is obtained solving Lp c c = MΓ c→s u c . So, we understand the global error as the mean difference between the result of an iteration using the staggered velocity and the collocated velocity.
The rate of convergence, p conv , is found fitting the points assuming a relation y = kx p conv , where k is a constant. As we can see in the next table, the rate of convergence value is around 1 for both schemes. Furthermore, the order of the error found with both schemes is similar, and it becomes practically the same for stretched meshes.

Conclusions
Apart from the well-known checkerboard problem, the appearance of unphysical velocities is a common problem found in highly distorted meshes. The volume weighted scheme solves this problem, at least for Cartesian meshes. A theoretical result could not be found for unstructured triangular meshes, but the accuracy for the resolution of one iteration of the Poisson equation has been tested. The results show that we have a first-order error solving the Poisson equation when approaching the staggered velocity using the volume weighted scheme. Furthermore, problems with extending the methodology are found when placing the cell-center outside the control volume, such as breaking the negative definiteness of L. Finally, it is interesting to remark that an interpolator that makes (33) positive definite (thus, guaranteeing stability) could not be unique, so other possibilites (also for unstructured meshes) may be found.