Numerical Performance of Higher-Order Semicompact Scheme for Arbitrary Triangular Cavity

An efficient fourth-order semicompact finite difference scheme has been developed to solve steady incompressible Navier-Stokes (N-S) equations in stream function and vorticity formulation in a triangular cavity of arbitrary geometry. The governing equations are transformed into curvilinear coordinates by a simple linear transformation to handle the nonregular geometry of the problem. The main feature of the new higher-order semicompact scheme is that it can calculate a triangle flow with arbitrary shape for high Reynolds numbers. It is found that the solutions obtained with the present scheme are in good agreement with the analytical results or with the existing results depending on the availability.


Introduction
Numerical simulations for the solution of steady incompressible viscous flow within a driven cavity are a complex and significant topic, and a lot of researchers contributed to this subject [1][2][3][4].Although there are still some minor discrepancies in the results, the square cavity flow has been essentially characterized [5].As pointed out in McQuain et al. [6], the results for the square cavity may not be applied to other important geometries such as a triangular cavity.Also the latter shapes are more common in practice.
The problem under consideration is the steady motion of an incompressible viscous flow in a triangular cavity of arbitrary geometry.This flow was studied numerically by McQuain et al. [6], Jyotsna and Vanka [7], Li and Tang [8], and Gaskell et al. [9].McQuain et al. [6] applied the Batchelor's mean square law to triangular cavity flow and analytically obtained the inviscid core vorticity for infinite Reynolds number.Recent calculations of the steady problem in an equilateral triangular cavity have been given in Ribbens et al. [10] for Re ≤ 500.Erturk and Gokcol [11] have presented high accurate, fine grid solutions of 2D steady incompressible flow in triangle cavities.The governing equations are solved up to very low residuals at various Reynolds numbers.The results of these studies, however, show some discrepancies.Furthermore, the scheme in [11] is only second-order spatial accuracy.This constitutes the main motivation of this study.
The main object of this study is the development of an accurate and efficient scheme for solving the N-S problems in triangular geometries.In the present paper, we have attempted an efficient fourth-order semicompact scheme (compactness of the scheme is relaxed for few terms of the governing equations) for stream function vorticity form of incompressible Navier-Stokes equations in curvilinear coordinates inside a triangular cavity.A geometric transformation handling triangles of arbitrary shape is also presented.The developed equations are solved numerically subject to appropriate boundary conditions by a fourth-order accurate finite difference method.Overall, besides opening up new possibilities, the method may be considered an efficient one for computation of flow for this physical configuration and the results thus obtained are not only useful for engineers in the process design; they also supplement the existing literature.

Governing Equations
A triangular driven cavity of general shape is given by locating its three vertices at Â(  , ℎ), B(  , ℎ), and Ĉ(  ,   ), with the International Journal of Engineering Mathematics upper side moving to the right via a constant velocity  0 ; see Figure 1(a).The 2D steady incompressible flow inside a triangular cavity is governed by the N-S equations.We consider the N-S equations in streamfunction and vorticity formulation in Cartesian coordinates, such that [5] where where  is the streamfunction,  is the vorticity,  and V are the components of the velocity in and -directions, respectively, and Re is Reynolds number.We note that these equations are non dimensional and a length scale of (ℎ−  )/3 and a velocity scale of  0 , that is, the velocity of the lid, are used to nondimensionalize the parameters and Reynolds number is defined accordingly.By a simple linear transformation the triangle  Â BĈ is transformed to a rightangled triangle with vertices (0, 0), (1, 0), and From these relations, we can calculate the transformation metrics as follows: and also Using the chain rule, the governing N-S equations (1) in general curvilinear coordinates are as follows: where Substituting for the transformation metrics, we obtain the equations that govern the flow in a triangular cavity shown in Figure 1(b) as follows: where and we obtain The boundary conditions in -plane are follows [11]: (, V) ⋅  = { 1 on the top side, 0 on the other two sides, (, V) ⋅  = 0 on all three sides, ( where  is the outside normal unit vector and  is the tangential unit vector in clockwise direction.In the -plane condition (11) remains the same form, while (13)-( 14) need to be converted as follows.On the top side Â B, substituting (11) into ( 13)- (14), with  = (1, 0) and  = (0, 1) yields On side B Ĉ we have Combining this with (11) we obtain Similarly on side Ĉ Â by using we obtain   = 0,   = 0 on side .
3. Higher-Order Semicompact Scheme where  is any known function, Δℎ is the constant step length in  and  directions, and If we define and then approximate the stream function (8) in the form of (25) we get where A 2D finite difference scheme is compact if it only involves at most the 8 nearest neighboring grid points (of the center point) in the approximation formula [14].Then in (23), first and second terms on the left side of the equation can be discretized on the compact nine point stencil as in any other compact scheme.However, discretization of the third term γ needs some special attention.We refer this term as 1 and the discretization of this term is considered after the discussion of the vorticity equation.

Discretization of the Vorticity Equation. Equation (9)
with fourth-order approximation can be written as The first and second terms on the left hand side of (25) can be discretized within the compact stencil.We call the third term γ as 2 which needs a larger stencil to approximate to fourth-order.The term on the right hand side of (25) contains The first two terms on the right hand side of (26) again can be discretized within the compact stencil to the fourth-order; International Journal of Engineering Mathematics however, the last four terms, call them as 3, 4, 5, and 6, need larger stencil.

Numerical Treatment of the Terms TERM1-TERM6
. We obtain from (8) that Therefore, Similarly, from (9) we obtain In (29), we have third derivatives such as   ( 2  +  2  ) and   ( 2  +  2  ) which present some problems to the fourthorder compact formulation.In order to find an expression for these derivatives, we use ( 8)-( 9) to obtain the following equations: where  = 2, 3, . . .,  and Employing Taylor series expansion, we get where  +  is the first-order forward difference operator in the -direction and Using (32) in (34), we get the fourth-order accuracy (Δℎ 4 ) expression as follows: where  +  is the first-order forward difference operator in the -direction,  2   is the first-order cross difference operator, that is Substituting (36) into (35), we get where Θ = ( ,2 − Δℎ) + (4 + ) ,2 ,  = 24 2 /Δℎ 2 .Similarly on side  the vorticity is equal to International Journal of Engineering Mathematics where  = 2, 3, . . ., .On side , it is equal to where  =  + 2 − ,  = 2, 3, . . ., .It is found that the point (, +1) is outside the 9-point domain being required and the value of  ,+1 is not available; we replace  ,+1 by using the following equation: At the three singular corners, the values of vorticity need some special attention in the process of discretization.We use the average of vorticity on the two adjacent points with the singular corners.

Results and Discussion
Following the numerical procedure described in the previous section, the new schemes (23) and (25) are fourth-order accurate; therefore, they need to be solved in an iterative manner such as SOR iteration with the relaxation parameter  = 2/3.In all of the cases considered, we start the iterations from a homogeneous initial guess and continue until a certain condition of convergence is satisfied.As convergence criteria we decided to use the difference of the - variables between two steps normalized by the previous value of the corresponding variable, such that max (       +1 −        ) < 10 −8 , max (       +1 −        ) < 10 −8 . (41) At this convergence level, this would indicate that the variables  and  are changing less than 0.000001% of their value between two iterations at every grid point in the mesh.These residuals indicate the degree to which the numerical solution has converged to steady state.We considered different triangle geometries with different Reynolds numbers.Numerical tests for a variety of triangular geometries have been investigated, with Re up to 2000 for an equilateral cavity and 1500 for scalene cavities.It should be pointed out that Reynolds number is based on ℎ, which is consistent with the definition in the case of the equilateral cavity.We solved the flow in this triangle cavity at various Reynolds numbers ranging between 1 and 2000.We note that if the length of one side of the triangle 2 √ 3 was used in nondimensionalization, as it was used by Erturk and Gokcol [11], then our Reynolds number of 2000 would 2 √ 3-fold such that it would correspond to a Reynolds number of 6928.

Equilateral Triangular Cavity. Let us first consider a nondimensional equilateral triangle with coordinates of corner points
which was also considered by McQuain et al. [6], Jyotsna and Vanka [7], Li and Tang [8], Gaskell et al. [9], and Ribbens et al. [10].We solved the flow in this triangle cavity at various Reynolds numbers ranging between 1 and 2000.
In order to verify the accuracy of the present numerical study, we consider arbitrary triangular cavity flow with the following analytical solutions [15]: the grid points on the boundaries.This way we would be able to avoid any effect of numerical boundary condition approximations on the numerical solution and concentrate on the accuracy of the solution of both formulations in the interior domain for given analytical values at the boundaries.We note that in (43), the vorticity  is independent of Re and the solution of the streamfunction  at different Reynolds numbers looks almost the same in a contour plot.We solve this model problem using different grid meshes, 11 × 11/2, 21 × 21/2, 41 × 41/2, 81 × 81/2, respectively.In these higherorder solutions, the average absolute (  ) error and the the maximum norm ( MAX ) error between the exact solution Exa  given in (43) and the numerical solution Sol  obtained from ( 23) and (25) are defined as where  being the total number of grid points.Since the grid size is decreased by a factor of 2, we can calculate the convergence rate  using the following formula [16]: As observed by Wan and Zhou [17], the accuracy of the present results is getting better as the grid number increases, and we can see that when the grid spacing is decreased progressively by half, the scheme maintains fourth-order of spatial accurate; the convergence rate is very close to  ≈ 4. In Figure 2, the average absolute error   and the the maximum norm error  MAX are plotted with respect to the grid spacing in a log-log scale [16].From Figure 2, we can again clearly see that the fully HOC formulation indeed provide fourth-order accurate solutions; the slope between log  and log ℎ is close to  ≈ 4.
Figures 3 and 4 show the streamline and vorticity contours of the triangular cavity flow for a variety of Reynolds numbers obtained with using (128 × 128)/2 grids points.These figures show the streams and vortices that as the Reynolds number increases.From these contour figures, we conclude that the fourth-order compact formulation provides very smooth solutions and it is seen that as Re becomes large, the location of the center of the primary eddy and its streamfunction value seem to have converged.In terms of quantitative analysis, Table 1 tabulates the center of the primary eddy and the streamfunction and vorticity values at the core, together with results found in the literature.Our results are in good agreement with the results in [6] and [9] up to the maximum Reynolds number (Re = 500).
In Figure 5(a), we plot the vorticity values at the center of the primary vortex tabulated in Table 1, with respect to Reynolds number.The results of Li and Tang [8] start to behave differently starting from Re = 100 from the rest of the results.Also in Figure 5(a), the results of McQuain et al. [6] start to deviate from present results and Gaskell et al. 's [9] results after Re = 200.Having a different behavior, the vorticity value of McQuain et al. [6] shows an increase such that their vorticity value at Re = 500 is greater than the vorticity value at Re = 350, whereas the present results and Erturk et al. 's [11] results show a continuous decrease until Re ≤ 1500.We believe that these behaviors are due to the coarse grids used in those studies and in order to resolve these behaviors we decided to solve the same flow problem using several coarse grid meshes.For a given grid mesh we have solved the flow for increasing Reynolds number until a particular Re, where the solution was not convergent but oscillating.For this particular Re when the number of grid points was increased, the convergence was recovered and we were able to obtain a solution.
According to the mean square law [18], the value of vorticity is approximately constant at the primary vortex.For an equilateral cavity with length of side 2 √ 3 this constant is International Journal of Engineering Mathematics found to be 1.054 by McQuain et al. [6].This infinite Reynolds number value of the vorticity of 1.054 is shown with the dotted line in Figure 5(a).We note that especially the eddies at the bottom corner occupy some portion of the corner as Reynolds number increases.However, the increase in the size of the portion of the bottom corner eddies almost stops after Re = 500 and the size of the primary eddy remains almost constant beyond and so is the value of the vorticity at the core of the primary eddy.It looks like, for an equilateral triangular cavity flow at high Reynolds numbers, the mean square law predicts the strength of the primary eddy within an error due to the effect of the secondary eddies.For circular or elliptic boundaries Ribbens et al. [10], for square cavity flow Erturk and Gokcol [11], and for rectangle cavities McQuain et al. [6] have shown that the mean square law is approximately valid and successful in predicting the strength of the primary eddy at high Reynolds numbers.Using this triangle geometry, we could carry out the calculations for Reynolds numbers as high as 800.In Figure 5(b), we plot the vorticity values at the center of the primary vortex obtained with present study and Li and Tang.[8].We can see from Figure 5(b) that the vorticity values are getting more and more close to the theoretical value as the Reynolds number increases.In Figure 6, we plot the streamlines and vorticity contours for Re = 350 and 500.Again it is seen that as Re becomes large, the location of the center of the primary eddy and its streamfunction value seem to have converged.

Isosceles Triangle Cavity.
We first considered an isosceles triangle which was also considered by Jyotsna and Vanka [7], Gaskell et al. [9], and Erturk and Gokcol [11], where We note that our definition of Reynolds number is equivalent to one fourth of the Reynolds number definition used by Jyotsna and Vanka [7]. Figure 7 shows the streamline contours of the flow in this triangle at various Reynolds numbers.Comparing the location of the primary eddy center with that of Jyotsna and Vanka [7], Gaskell et al. [9], and Erturk and Gokcol [11], we believe that our results are more accurate.

International Journal of Engineering Mathematics
We also considered an isosceles right triangle with the 90 ∘ corner being at the top right corner, such that with corner points [11]   = 0,   = 1, ℎ = 1,   = 1,   = 0. (49) Figure 9 shows the flow topology as a function of Reynolds numbers and at last, we also considered an isosceles right triangle with the 90 ∘ corner being at the top left corner, such that with corner points [11]   = 0,   = 1, ℎ = 1,   = 0,   = 0.
Using this triangle geometry, we were able to obtain the solutions and Figure 10 shows the flow topology as a function of Reynolds numbers.For the considered triangle geometry, we can see that the eddy is closest to the moving lid.As it is obvious from these figures, in both cases, the flow behaves very differently as Reynolds number increases, which shows that the flow structures in a triangle cavity are greatly affected by the triangle geometry.

Conclusion
In this work we have developed a new semicompact fourthorder scheme for the time-independent - form of the 2D, incompressible N-S equations governing the fluid flow in a triangular driven cavity.Our numerical scheme has been proved robust for a wide range of Reynolds numbers and applicable to triangular cavities with arbitrary shape.The key point with the present scheme is that it allows direct iteration for low-to-high Reynolds numbers.It is the success of the current method with the wide range of Re and mesh sizes that indicates the potential of this method as an accurate and stable numerical method applicable to a wide range of problems.We have tested the present method for both the equilateral triangular cavity problem and the scalene triangular cavity problem, and excellent agreement is found in all the cases, both qualitatively and quantitatively.

Figure 1 :
Figure 1: Geometric transformation of the triangular cavity: Physical domain (a) and Computational domain (b).

Figure 2 :
Figure 2: The average absolute error   and the maximum norm error  MAX in streamfunction and vorticity variables with respect to grid spacing.(a) streamfunction and (b) vorticity.

Figure 3 :
Figure 3: Streamline contours of triangular cavity flow at various Reynolds numbers.

Figure 4 :
Figure 4: Vorticity contours of triangular cavity flow at various Reynolds numbers.

Figure 5 :
Figure 5: Comparison of the vorticity at the center of the primary eddy in the equilateral triangular cavity (a) and right-oriented triangular cavity (b).

Figure 6 :
Figure 6: Streamlines distributions for a right-oriented triangular cavity at different Reynolds numbers.

Figure 8 :
Figure 8: Contour figures of isosceles triangle with various corner angle.

Figure 9 :Figure 10 :
Figure 9: Contour figures of right hand side aligned right triangle at different Reynolds numbers.

Table 1 :
Properties of the center of the primary eddy, located at (, ) with streamfunction value and vorticity for equilateral triangle.