A New Method for Constructing the Coefficients of Pressure Correction Equation for Colocated Unstructured Grids

Abstract: One of the important equations in numerical solution of Navier Stokes equations is the pressure correction equation. In this article, a new method for constructing the coefficients of this equation for colocated unstructured meshes is proposed. This method is based on momentum interpolation. The method is compared with the approach adopted in literature. The complete discretization of the Navier Stokes equations using finite volume method is presented. An algorithm similar to SIMPLE is used to evaluate the rate of convergence on two sample problem. The results show that by using rectangular grids, two methods have the same performance, but when triangular meshes are applied, the new method increases the rate of convergence.


INTRODUCTION
In the last two decades, solution of Navier Stokes equations using colocated arrangement has received great interests. On the other hand, unstructured grids are popular for solution of flow field in complex geometries. Colocated arrangement has some obvious advantages over staggered grids especially in case of non-orthogonal meshes. In colocated arrangement, all variables share the same location, hence only one set of volumes is considered. Also, by applying the colocated grids, the convection contribution to the coefficients in discretized equations, are the same for all variables. Finally, Cartesian velocity components can be used in conjunction with non-orthogonal grids for complex geometries which result simpler discretized equations.
After the original work by Rhie and Chow [1] , Peric [2] generalized the same idea for three dimensional flows and calculated several two and three dimensional flow situations. Majumdar [3] studied the role of under relaxation parameter in momentum interpolation for calculation of flow with colocated grids. Lien [4] used the colocated arrangement for unstructured grids successfully. He applied the momentum interpolation as a concept for derivation of pressure correction equation. However, different authors have different ideas to construct the pressure correction equation. For example the approach used by Thomadakis et al. [5] for constructing the coefficients is quiet different from that used by Lien [4] . Here the emphasis is on the fact that these coefficients can affect the numerical solution. Therefore, a new form of the coefficients is proposed and it's rate of convergence and performance is examined.
In this study, two examples are solved by structured and unstructured grids: • Laminar flow in a lid-driven square cavity • Two dimensional parallel flow

FLOW EQUATIONS
For incompressible Newtonian fluid flows, the conservation equation for mass and momentum are as: In above equation, u and v are Cartesian components of velocity vector.
Transport equation for a general variable, such as ϕ, can be written as:

DATA STRUCTURE
As mentioned before, colocated arrangement for velocity components and pressure within an arbitrary finite volume is adopted here. The basic idea of this data structure is shown in Fig. 1. As shown, forming points are the vertices of control volume and the center of control volume is assumed for storage of all variables.
In Fig. 1, each cell has some adjacent cells with two or three shared nodes. For example, in triangular meshes, the adjacent cells have two shared nodes and in tetrahedral meshes, they have three shared nodes.
The forming points and adjacent of a cell are numbered in the counterclockwise direction suitable for divergence theorem.

FINITE VOLUME DISCRETIZATION
By integrating the Eq. 5, over a control volume, using divergence theorem, the general transport equation in integral form is obtained: where, S is the surface vector with the Cartesian components, namely S x , S y . This integral form consists of four parts; transient term, convection term, diffusion and source term.
On triangular meshes, the convective term can be approximated as: The parameter C si is the mass flux over the i'th surface of control volume. The mass fluxes are: Variables u si and v si are the components of velocities at the i th face of the control volume.
The value of ϕ on the faces (i.e. ϕ si ) can be approximated by the following second-order upwind scheme [1] : The vectors of 0 c si r − and ci si r − are shown in Fig. 2.
In Fig. 2, C 0 is the Center of the control volume Ω and C i is the centre of its i'th adjacent cell Diffusion term can be approximated as follows: The normal gradient of any variable can be approximated as below: where, h C0-Ci is the summation of distances h C0-Si and h Ci-Si in Fig. (2). Inserting Eq. 10 in Eq. 9, results: The S i is the magnitude of face vector i S . The pressure source terms can be approximated by using Green's theorem. For example: The transient term is approximated as: In which, superscript (n-1) is used previous time step.
Insertion of the above discretized terms in the integral form of transport equation, gives: The pressure source term can be obtained from: Convection and diffusion source terms are: In above equations si is diffusion coefficient at i'th face and S i is the magnitude of a face vector and S ϕi denotes that component of face vector which is parallel with ϕ. For example in calculation of S pu , the term S ϕi is S xi and for calculation of S pv , the term S ϕi is S yi .

PRESSURE CORRECTION EQUATION
The pressure correction equation is derived from integral form of the continuity equation. This integral form of can be written as: The pressure and velocity link is very important. According to the momentum interpolation method, the face velocities are approximated as: By considering below notations (Eq. 15): Using momentum interpolation concept, the face velocity can be written as: Different approximations can be employed for estimating the face value for si 0 ( ) A Ω in the above equations.
Lien [4] , proposed the face value of The weight factor (β) is defined by: By replacing the above approximations, the integral form of continuity equation can be written as: In the other form: And by introducing perturb mass flux (C' Si ) as: The continuity equation becomes: Therefore, the final form of pressure correction is: Where: And: According to the Lien [4] Proposal, The A ip coefficients are: And according to Eq. 24, The New form for A ip coefficients can be written as: New pressure and velocity field can be calculated from: and where, α is under relaxation factor for p'. In this study, the value of 0.3 is used as relaxation factor for p'. The algorithm of solution is: • Guessing the pressure field p • Calculation of mass flux (Eq. 7) • Calculation of coefficients (Eq. 15, 16, 17) • Calculation of velocities (Eq. 14) • Calculation of pressure correction (Eq. 30, 36) • Calculation of u and v from their starred values using the velocity correction formulas (Eq. 36, 37) • Assuming the corrected pressure as a new guessed pressure and returning to step (2), then repeating the whole procedure until a converged solution is obtained This is quite similar to SIMPLE method used by Patankar [6] for staggered grids.

CONVERGENCE CRITERIA'S AND DEFINITION OF RESIDUALS
After discretization, the residual for a general variable ϕ at a cell can be written as: A different definition is used for residual of continuity equation, which is defined by: To solve the problem, two different meshes are generated (Fig. 3). The problem is solved using two different pressure correction equations mentioned in previous section .
The pressure contours and the streamlines for structured grid are shown in Fig. 4. These results are achieved by using the new method (based on Eq. 34) for pressure correction equation.
In Fig. 5a, dimensionless profiles of u component of velocity vector at the midline of the cavity (at x = 0.5) for both structured and unstructured meshes are compared with Ghia et al. results [7] . The results are achieved using Eq. 34. There is a satisfactory agreement between the results, but in triangular mesh results deviation is more. Some error resources such as non-orthogonality are important for triangular meshes. Using deferred correction term can reduce the effects of non-orthogonality [8] .
In Fig. 5b, profiles of u component of velocity vector are plotted for different methods on structured grids. Figure 5b also shows that the results are equal for two methods.
Total mass imbalance (defined by equation 40), vs. the number of iterations are plotted for different cases in Fig. 6a and b. The important feature of these results is that, the values of total mass imbalance (b) can be interpreted as a good criterion for overall convergence of the problem. The only difficulty in this way is that  [7] (Re = 1000), b: between Lien [4] method and new proposed method the value of total mass imbalance can vary over a wide range for different problems, especially in unstructured grids. As shown in this Fig. 6a, the value of total mass imbalance are equal for two methods in rectangular meshes, but new method needs less number of iterations to converge in triangular meshes (Fig. 6b).
It can be understood that the rate of convergence by using new method (Eq. 34) is more, in the case of triangular meshes. The study of the results for parallel flow between two plates will be presented in the subsequent section as another sample problem. Two-dimensional parallel flow: This case has an analytical solution in fully developed region. The geometry and properties of the fluid are shown in Fig.  7.
For solving the problem, two types of rectangular and triangular meshes are generated. The specifications of the meshes are given in Table 1.  In Fig. 8, contours of u and v components of the velocity are shown. These results are achieved using rectangular mesh and new method pressure correction equation.  [9] . a: at x = 0.5H; b: at x = 2H results [9] in two sections of developing region. This results show the correctness of written code.
Total mass imbalance (defined by equation 40), vs. the number of iterations are plotted for different cases in Fig. 10a, b.
As shown in this Fig. 10a, the value of total mass imbalance are equal for two methods in rectangular meshes, but in triangular meshes new method needs less number of iterations to converge (Fig. 10b).

SUMMARY AND CONCLUSION
By studying two sample problems, it is evident that new method (defined by Eq. 34) can decrease the residuals faster than by Lien [4] proposed method. Two methods are approximately equal when rectangular meshes are used (Fig. 6a, 10a), but by using triangular meshes, new method needs less number of iterations to converge (Fig. 6b, 10b). This phenomenon is due to non-orthogonality of triangular meshes in which, weak pressure correction can't damp the oscillation generated during the solution procedure. Therefore the method of construction of coefficients in Pressure correction equation can affect the rate of convergence. It must be emphasized that, using new method, the calculation time is less, because the residual of continuity equation and total mass imbalance decrease faster for this approach (Fig. 6b, 10b).