{'fraction_non_alphanumeric': 0.04210790149758757, 'fraction_numerical': 0.020490005639451095, 'mean_word_length': 4.6, 'pattern_counts': {'":': 0, '<': 0, '<?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': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fourth Order Compact Formulation of Navier-Stokes Equations and Driven Cavity Flow at High Reynolds Numbers

A new fourth order compact formulation for the steady 2-D incompressible Navier-Stokes equations is presented. The formulation is in the same form of the Navier-Stokes equations such that any numerical method that solve the Navier-Stokes equations can also be applied to this fourth order compact formulation. In particular in this work the formulation is solved with an efficient numerical method that requires the solution of tridiagonal systems using a fine grid mesh of 601x601. Using this formulation, the steady 2-D incompressible flow in a driven cavity is solved up to Reynolds number of 20,000 with fourth order spatial accuracy. Detailed solutions are presented.


INTRODUCTION
High Order Compact (HOC) formulations are becoming more popular in Computational Fluid Dynamics (CFD) field of study. Compact formulations provide more accurate solutions in a compact stencil.
In finite differences, a standard three point discretization provides second order spatial accuracy and this type of discretization is very widely used. When a high order spatial discretization is desired, ie. fourth order accuracy, then a five point discretization have to be used. However in a five point discretization there is a complexity in handling the points near the boundaries.
High order compact schemes provide fourth order spatial accuracy in a 3×3 stencil and this type of compact formulations do not have the complexity near the boundaries that a standard wide (five point) fourth order formulation would have. Dennis and Hudson [1], MacKinnon and Johnson [2], Gupta et al. [3], Spotz and Carey [4] and Li et al. [5] have demonstrated the efficiency of the high order compact schemes on the streamfunction and vorticity formulation of 2-D steady incompressible Navier-Stokes equations.
In the literature, it is possible to find numerous different type of iterative numerical methods for the Navier-Stokes equations. These numerical methods, however, could not be easily used in HOC schemes because of the final form of the HOC formulations used in [1], [2], [3], [4] and [5]. This fact might be counted as a disadvantage of HOC formulations that the coding stage is rather complex due to the resulting stencil used in these studies. It would be very useful if any numerical method for the solution of Navier-Stokes equations described in books and papers could be easily applied to high order compact (HOC) formulations.
In this study, we will present a new fourth order compact formulation. The difference of this formulation with References [1][2][3][4][5] is not in the way that the fourth order compact scheme is obtained. The main difference, however, is in the way that the final form of the equations are written. The main advantage of this formulation is that, any iterative numerical method used for Navier-Stokes equations, can be easily applied to this new HOC formulation, since the final form of the presented HOC formulation is in the same form with the Navier-Stokes equations. Moreover if someone already have a second order accurate (O∆x 2 ) code for the solution of steady 2-D incompressible Navier-Stokes equations, using the presented formulation, they can easily convert their existing code to fourth order accuracy (O∆x 4 ) by just adding some coefficients into their existing code. With this new compact formulation, we have solved the steady 2-D incompressible driven cavity flow at very high Reynolds numbers using a very fine grid mesh to demonstrate the efficiency of this new formulation.

FOURTH ORDER COMPACT FORMULATION
In non-dimensional form, steady 2-D incompressible Navier-Stokes equations in streamfunction (ψ) and vorticity (ω) formulation are given as 1 Re where x and y are the Cartesian coordinates and Re is the Reynolds number. For first order and second order derivatives the following discretizations are fourth order accurate where φ x and φ xx are standard second order central discretizations such that If we apply the discretizations in Equations (3) and (4) to Equations (1) and (2), we obtain the following equations 1 Re In these equations we have third and fourth derivatives (∂ 3 /∂x 3 , ∂ 3 /∂y 3 , ∂ 4 /∂x 4 and ∂ 4 /∂y 4 ) of streamfunction and vorticity (ψ and ω) variables. In order to find an expression for these derivatives we use Equations (1) and (2). For example, when we take the first and second x-derivative (∂/∂x and ∂ 2 /∂x 2 ) of the streamfunction equation (1) we obtain the following equations And also, by taking the first and second y-derivative (∂/∂y and ∂ 2 /∂y 2 )o ft h e streamfunction equation (1) we obtain the following equations Using standard second order central discretizations given in Table I, these equations ((9), (10), (11) and (12)) can be written as the followings When we substitute Equations (14) and (16) into Equation (7) we obtain the following finite difference equation.
We note that, the solution of Equation (17) is also a solution to streamfunction equation (1) with fourth order spatial accuracy. Therefore if we numerically solve Equation (17), the solution we obtain will satisfy the streamfunction equation up to fourth order accuracy.
In order to obtain a fourth order approximation for the vorticity equation (2), we follow the same procedure. When we take the first and second derivatives of the vorticity equation (2) with respect to x-and y-coordinates we obtain the followings If we substitute Equations (18) and (20) for the third derivatives of vorticity (∂ 3 ω/∂x 3 and ∂ 3 ω/∂y 3 ) into Equations (8), (19) and (21) and also if we substitute Equations (13) and (15) for the third derivatives of streamfunction (∂ 3 ψ/∂x 3 and ∂ 3 ψ/∂y 3 ) into Equations (8), (19) and (21) and finally if we substitute Equations (19) and (21) for the fourth derivatives of vorticity (∂ 4 ω/∂x 4 and ∂ 4 ω/∂y 4 ) into Equation (8), then we obtain the following equation Again we note that, the solution of Equation (22) satisfy the vorticity equation (2) with fourth order accuracy.
As the final form of our HOC scheme, we prefer to write Equations (17) and (22) as the following We note that, the finite difference Equations (23) and (24) are fourth order accurate O(∆x 4 , ∆x 2 ∆y 2 , ∆y 4 ) approximation of the streamfunction and vorticity equations (1) and (2). In Equations (23) and (24), however, if A, B, C, D, E and F are chosen to be equal to zero then the finite difference Equations (23) and (24) simply become 1 Re Equations (26) and (27) are the standard second order accurate O(∆x 2 , ∆y 2 ) approximation of the streamfunction and vorticity equations (1) and (2). When we use Equations (23) and (24) for the numerical solution of 2-D steady incompressible Navier-Stokes equations, we can easily switch between second and fourth order accuracy just by using homogeneous values for the coefficients A, B, C, D, E and F or by using the expressions defined in Equation (25) in the code.
In Equations (23), (24) and (25) instead of finite difference discretizations, if we substitute for partial derivatives we obtain the following differential equations We note that the numerical solutions of Equations (28) and (29), strictly provided that second order discretizations in Table I are used and also strictly provided that a uniform grid mesh with ∆x and ∆y is used, are fourth order accurate to streamfunction and vorticity equations (1) and (2). We prefer to call Equations (28) and (29) Fourth Order Navier-Stokes (FONS) equations. The only difference between FONS equations (28) and (29) and Navier-Stokes (NS) equations (1) and (2) are the coefficients A, B, C, D, E and F . In fact the NS equations are a subset of the FONS equations. We note that FONS equations (28) and (29) are in the same form with Navier-Stokes (NS) equations (1) and (2), therefore any iterative numerical method (such as SOR, ADI, factorization schemes, pseudo time iterations and etc.) used to solve streamfunction and vorticity equations (1) and (2) can also be easily applied to fourth order equations (28) and (29). Moreover, any existing code that solve the streamfunction and vorticity equations with second order accuracy can easily be modified to provide fourth order accuracy just by adding the coefficients A, B, C, D, E and F into the existing code to obtain the solution of FONS equations. Of course, when the coefficients A, B, C, D, E and F are added into a second order accurate code to obtain fourth order accuracy, evaluating these coefficients would require extra CPU work. This might be considered as the cost of increasing accuracy from second order to fourth order.

NUMERICAL METHOD
Recently Erturk et al. [6] have presented a new, stable and efficient numerical method that solve the streamfunction and vorticity equations. The numerical method solve the governing steady equations through iterations in the pseudo time. In this study, we will apply the numerical method Erturk et al. [6] have proposed, to FONS equations (28) and (29) and solve the steady driven cavity flow with fourth order accuracy. For details about the numerical method, the reader is referred to Erturk et al. [6]. When we apply the numerical method to Equations (28) and (29), we obtain the following equations 1 − ∆t ∂ 2 ∂x 2 1 − ∆t ∂ 2 ∂y 2 ψ n+1 = ψ n +∆tω n − ∆tA n + ∆t The solution methodology of these two equations are quite simple. First the streamfunction equation (31) is solved in two steps. For streamfunction equation, a new variable f is defined as the following Using this variable in Equation (31) we obtain the following equation In this equation, the only unknown is the variable f . We first solve this equation for f by solving a tridiagonal system. After this, when we obtain the value of f at every grid point we solve Equation (33) for streamfunction (ψ n+1 ) by solving another tridiagonal system.
After solving the streamfunction equation (31), we solve the vorticity equation (32). For this, similarly, we introduce a new variable g which is defined as the following Using this variable in Equation (32), we obtain the following equation In this equation the only unknown is the variable g. By solving a tridiagonal system, we obtain the value of g at every grid point. Then we solve Equation (35) for vorticity (ω n+1 )b y solving another tridiagonal system.
In a compact formulation, the stencil have 3×3 points. The solution at the first diagonal grid points near the corners of the cavity would require the vorticity values at the corner points. However, the corner points are singular points for vorticity. Gupta et al. [7] have introduced an explicit asymptotic solution in the neighborhood of sharp corners. Similarly, Störtkuhl et al. [8] have presented an analytical asymptotic solutions near the corners of cavity and using finite element bilinear shape functions they also have presented a singularity removed boundary condition for vorticity at the corner points as well as at the wall points. We follow Störtkuhl et al. [8] and use the following expression for calculating vorticity values at the wall where V is the speed of the wall which is equal to 1 for the moving top wall and equal to 0 for the three stationary walls. For corner points, we use the following expression for calculating the vorticity values where again V is equal to 1 for the upper two corners and it is equal to 0 for the bottom two corners. The reader is referred to Störtkuhl et al. [8] for details.

RESULTS AND DISCUSSIONS
The schematics of the driven cavity flow is given in Figure 1. In this figure the abbreviations BR, BL and TL refer to bottom right, bottom left and top left corners of the cavity, respectively. The number following these abbreviations refer to the vortices that appear in the flow, which are numbered according to size. For every Reynolds number considered, we have continued our iterations until, in the computational domain both the maximum residual of Equations (23) and (24), which are given as are less than 10 −10 . Such a low value is chosen to ensure the accuracy of the solution. At these residual levels, the maximum absolute change in streamfunction value between two time steps, (max(|ψ n+1 − ψ n |)), was in the order of 10 −16 and for vorticity, (max(|ω n+1 − ω n |)), it was in the order of 10 −14 . Obviously these convergence levels are far more less than satisfactory, however such low values demonstrate the efficiency of the numerical method used in this study which was presented by Erturk et al. [6].
Using an efficient numerical method, Erturk et al. [6] have clearly shown that numerical solutions of driven cavity flow is computable for Re >10,000 when a grid mesh larger than 256×256 is used. With a grid mesh of 601×601 Erturk et al. [6] have solved the cavity flow up to Re=21,000 using the numerical method also used in this study. In order to be able to obtain solutions at high Reynolds numbers, following Erturk et al. [6], in this study we have used a large grid mesh with 601×601 grids. With this many number of grid points we obtained steady solutions of the cavity flow up to Re=20,000 with fourth order accuracy. Figures 2 to 6 show the streamfunction and vorticity contours of the driven cavity flow between Re=1,000 and Re=20,000. These figures show the vortices that are formed in the flow field as the Reynolds number increases. From these contour figures, we conclude that the fourth order compact formulation provides very smooth solutions.
In Figure 7 we plot a very enlarged view of the top right corner (where the moving wall moves towards the stationary wall) of the streamfunction contour plot for the highest Reynolds number considered, Re=20,000. In this figure the dotted lines show the grid lines. As it is seen in this enlarged figure, fourth order streamfunction contours are very smooth even at the first set of grid points near the corners. Table II tabulates the streamfunction and vorticity values at the center of the primary and secondary vortices and also the location of the center of these vortices for future references. This table is in good agreement with that of Erturk et al. [6].
Using Richardson extrapolation on the solutions obtained with different grid meshes, Erturk et al. [6] have presented theoretically fourth and sixth order accurate (O∆x 4 and O∆x 6 ) streamfunction and vorticity values at the center of the primary vortex. Table III and IV compares the forth order compact scheme solutions of the streamfunction and the vorticity values at the center of the primary vortex with the fourth order (O∆x 4 ) Richardson extrapolated solutions tabulated in Erturk et al. [6]. The present solutions and the solutions of Erturk et al. [6] agree with each other.

CONCLUSIONS
In this study a new fourth order compact formulation is presented. The uniqueness of this formulation is that the final form of the HOC formulation is in the same form of the Navier-Stokes equations such that any numerical method that solve the Navier-Stokes equations can be easily applied to the FONS equations in order to obtain fourth order accurate solutions (O∆x 4 ). Moreover with this formulation, any existing code that solve the Navier-Stokes equations with second order accuracy (O∆x 2 ) can be altered to provide fourth order accurate (O∆x 4 ) solutions just by adding some coefficients into the code at the expense of extra CPU work of evaluating these coefficients.
In this study, the presented fourth order compact formulation is solved with a very efficient numerical method introduced by Erturk et al. [6]. Using a fine grid mesh of 601×601, as it was suggested by Erturk et al. [6] in order to be able to compute for high Reynolds numbers, the driven cavity flow is solved up to Reynolds numbers of Re=20,000. The solutions obtained agree well with previous studies. The presented fourth order accurate compact formulation is proved to be very efficient.