Strojniški vestnik - Journal of Mechanical Engineering 53(2007)2, 105-113 UDK - UDC 519.65:004.021 Kratki znanstveni prispevek - Short scientific paper (1.03) Vpliv pravokotnosti mreže na konvergenco programa SIMPLE za reševanje Navier-Stokes-ovih enačb The Influence of Grid Orthogonality on the Convergence of the SIMPLE Algorithm for Solving Navier-Stokes Equations Ivo Džijan - Zdravko Virag - Hrvoje Kozmar (University of Zagreb, Croatia) Razvili smo metodo končnih volumnov za reševanje Navier-Stokesovih enačb na lokalno pravokotni nestrukturirani mreži z uporabo algoritma SIMPLE. Razvito metodo smo primerjali s podobno metodo na strukturirani, ne nujno pravokotni mreži, v členih pretekle konvergence in obsegu pod-relaksacijskih faktorjev, pri katerih se metodi približujeta. Kadar je strukturirana mreža pravokotna, sta stopnji približevanja obeh metod podobni. V primerih, kadar strukturirana mreža ni pravokotna, se pokaže prednost predlagane metode pri lokalno pravokotni mreži v razmerah pretekle konvergence. V teh primerih je obseg pod-relaksacijskih faktorjev, pri katerih je predlagana metoda zadovoljivo konvergentna, mnogo večji kot pri metodi na nepravokotni mreži. © 2007 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: Navier-Stokesove enačbe, metode končnih volumnov, algoritmi SIMPLE, nestrukturirane mreže) A finite-volume method for solving the Navier-Stokes equations on a locally orthogonal unstructured grid using the SIMPLE algorithm has been developed. The developed method was compared with a similar method on a structured, not necessarily orthogonal grid, in terms of convergence history and the range of under-relaxation factors in which the methods converge. When the structured grid is orthogonal, the convergence rates of the two methods are similar. For the cases when the structured grid is non-orthogonal, the superiority of the proposed method on the locally orthogonal grid is demonstrated in terms of convergence history. In these cases, the range of under-relaxation factors in which the proposed method shows satisfactory convergence is much wider than for the method on the non-orthogonal grid. © 2007 Journal of Mechanical Engineering. All rights reserved. (Keywords: Navier-Stokes equations, finite volume methods, SIMPLE algorithm, unstructured grid) 0 INTRODUCTION The rapid development of computers has brought about rapid developments in the field of computational fluid dynamics. Calculation domains are now more complex, which increases the need to use an unstructured grid for their discretization. Finite volume methods are widely applied for solving fluid flow problems. Initially, these methods were used on structured staggered grids. Nowadays they are used on unstructured collocated grids, on which segregate algorithms with the pressure-based approach are applied for incompressible flows. The most popular algorithm based on pressure correction is the SIMPLE (Semi-Implicit Method for the Pressure- Linked Equation) algorithm, Caretto et al. [1] and Patankar and Spalding [2]. In the pressure-velocity correction relation the effects coming from velocity corrections in neighboring nodes on the pressure correction in the central node are neglected. The consequence of this neglecting is the overestima-tion of the pressure correction, which can cause the divergence of the numerical process. To ensure the stability of the numerical process, the under-relax-ation factor for the pressure is introduced. The optimal value of this factor cannot be estimated in advance since it depends on the grid’s characteristics and the nature of the problem. The SIMPLE algorithm is originally defined on a staggered grid where the pressure is calculated 105 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)2, 105-113 in the cell centre and the velocity components are calculated on the cell faces. On a collocated grid, the pressure field and velocity components are calculated in the cell centre. The application of the SIMPLE algorithm on a collocated grid started with Rhie and Chow [3]. The grid non-orthogonality is one of the factors that increases the number of iterations of the SIMPLE algorithm. If the connecting line of two neighboring nodes is not perpendicular to the cell face, some terms that appear due to non-orthogonality are usually neglected. This is the case with the CAFFA public-domain computer code [4]. It is believed that this neglecting slows down the rate of convergence of the numerical method. Therefore, the modification of the finite volume method on an unstructured locally orthogonal grid is proposed. The rate of convergence of the SIMPLE algorithm on that grid will be compared with the rate of convergence on a structured, not necessarily orthogonal grid. 1 MATHEMATICAL MODEL AND NUMERICAL PROCEDURE The mathematical model of steady, laminar, incompressible fluid flow with constant viscosity and without mass forces is adopted. The model is described with the following Navier-Stokes equations: d OX (pVj) = 0 d t dp d v. ------I/7V.V,. ) =-----— + JU-------- ôx dx ôxôx (1) (2) Fig. 1. A part of calculation domain and a typical cell of locally-orthogonal unstructured grid where p, v., p, ju and x. are the fluid density, velocity, pressure, 'viscosity and coordinates, respectively. These equations will be numerically solved on an unstructured locally orthogonal grid. A part of such a computational grid is shown in Fig. 1. The main nodes, C and N, at which the velocity and pressure fields are calculated, are placed within their respective cells. The connecting line CN is perpendicular to the cell face and the nodes C and N are at an equal distance from an auxiliary node n, which enables a simple formulation of the high-order interpolation. It is clear that such a grid is possible in every 2D case, because the cell vertex a in Fig. 1 is the circumcenter of the triangle formed by the nodes C, M and K. Such a grid generator is described by Džijan [5], which includes a grid-smoothing procedure that forces the main nodes to be close to the cell centroids and the auxiliary nodes to be close to the cell-face centroids. Discretization of the equations starts with integrating over the cell volume V, according to Fig. 1. By using the Gauss theorem and the mean-value theorem, the integrated governing equations take the form: m L [F] =0 (3) Fv\ mA cV. dx -V dp dx mA dv dx (4) where F = p A v. «. = p A vn is the mass flow through the cell face and J, = F vt n - juA (dv, / dXj ) n «y. is the momentum flux through the cell face. A and « are the cell-face area and its outward normal vector, V is the cell volume, while k denotes the cell-face index and m is the number of cell faces on the considered volume. The scope of the differencing schemes is to define the velocity v;|n and its normal derivatives at the auxiliary node n in terms of velocity values at the main nodes. Since the adopted grid is locally orthogonal, these values are defined by using only the values at nodes C and N. A blending scheme of the central differencing scheme (CDS) and of the first-order upwind differencing scheme (UDS) is used in the CAFFA computer code. Therefore, the same scheme will also be used in the proposed method. In the case of the locally orthogonal grid, the diffusion part of the flux vector is modeled with the following equation: Jd -mA dvi dx -mA -mA 2s (5). 106 Džijan I. - Virag Z. - Kozmar H. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)2, 105-113 In the case of a non-orthogonal grid, an additional term appears. In the CAFFA computer code this term is implemented by using the deferred correction approach, i.e., by using the velocity values from the previous iteration. In the first-order UDS, the convective flux is modelled by: F viC for F>0 F viN for F<0 (6), and in the CDS (for the case when the node n lies in the middle of the CN connection line) by: •^1MC+4N ) (7). By introducing the mixing factor g, the final expression for the momentum flux is: T t1 \ rUDS , xCDS , xd (8) where, for g = 0 the result is the UDS, and for g = 1 the CDS. Introducing the expressions for the fluxes into (4) results in: where: E[aN^|N_ max) -F,0] -V dp dx mA 2s m aC = Xk] + bi (9) (10), (11) and b^r\l\12Fvi] -E[max(-*,0)]\lC + E[max(-*,0)v.Nî} L" " " (12). It is obvious that all the terms coming from the CDS are treated as a deferred correction, the same as in the CAFFA code. To reduce the possibility that the numerical process diverges, this equation is under-relaxed in the following form: dx u aCvi\C L>,+--------- (13) which was proposed by Patankar [6]. The last term on the right-hand side is calculated from the previous iteration, and auv is the under-relaxation factor for the velocity. According to Rhie and Chow [3], the mass flow through the cell face is defined as follows: F = pAvn=pA{vn)-pA {n} dp dn (dp n] (14) where the line above a symbol indicates the linear interpolation between the values at nodes C and N, as follows: ( v ) = — (v.I +V.I h and Udp 2{dn dp dn (15), (16) (17). C CNy In the case of a locally orthogonal grid, the normal derivative of the pressure is defined by using the CDS, as follows: K dp dn pC 2s (18). In the CAFFA code, where the grid is non-orthogonal, additional terms emerge and are treated explicitly by using the values from the previous iteration. Solving the momentum equation with a given pressure field/?* results in the velocity field v *, and the mass flow F, which does not necessarily satisfy the continuity equation. For that reason, the velocity corrections v ’ and the corresponding F’ and pressure corrections ' are searched, so that the corrected velocity field v. = v.* + v ’ and corrected mass flow F = F + F’ satisfy the continuity equation. According to Equation (14), the corrected mass flow is approximated as follows: F= *EäK 2s ^an PC F -< Ä-fC (19). Introducing the corrected mass flows in the continuity equation (3) results in the following equation for the pressure correction: where and m i "c/'c-IKK] =b"' aC =Z[flN m bp'=-np k (20) (21) (22). The solving of this equation results in the pressure-correction field. Therefore, the pressure Vpliv pravokotnosti mreže - The Influence of Grid Orthogonality 107 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)2, 105-113 field is corrected using the following equation: pC = p*C+app'C (23) where a is the under-relaxation factor for the pressure. The velocities in the main nodes are corrected as follows: viC =vi C aC dXj (24). The gradients of the physical values in the main nodes are calculated using the Gauss formula as follows: dx 1v±WA 2> k=1 (25) where f can stand for vi, p or p’. The steps in the SIMPLE algorithm for solving the Navier-Stokes equations can be summarized as follows: 1. Guess the pressure field p* and the velocity field vi. 2. Solve the momentum equation (13) to obtain vi*. 3. Solve the p’ equation (20). Correct the pressure according to (23), correct the velocity according to (24) and the mass flow according to (19). 4. Treat the corrected pressure as p* and return to Step 2. Repeat the whole procedure until a converged solution is obtained. The converged solution is obtained when the normalized residuals for the continuity and momentum equations become smaller than some small number, e. In this paper e = 10-6 was used. The residual for the continuity equation is: R M ,_ and the residual for the momentum is: *,=z LMn]V (26) (27). In the above expressions, l denotes the cell index and M the total number of cells. The values of the variables in the above formula are from the current iteration, and the coefficients are prepared for the next iteration. The following residuals are usually normalized: the mass residuals with the inlet mass-flow rate and the residuals for the momentum equation with the inlet momentum flow rate. 2 RESULTS The described numerical method is implemented in the FVM computer code. In this code the residuals are defined and normalized in the same way as in the CAFFA code. The rate of convergence of the described method and of the method used in the CAFFA code will be compared by varying the grid’s non-orthogonality and the differencing scheme. Also, the range of under-relaxation factors in which the numerical procedure converges will be analyzed. 2.1 Laminar flow in a lid-driven cavity with inclined side walls In this test the 2D laminar fluid flow is calculated in a closed cavity whose lid is moving in a tangential direction with velocity v Peric [7]. The Reynolds number based on the side length a is Re = p.v .a/'u= 1000. The calculation is performed for different inclination angles of the side walls, /?= 90°, 67.5° and 45°. In this problem the residuals defined by (26) and (27) are not normalized. Fig. 2 shows the qualitative picture of the streamlines for ß= 90° and 45°. It is obvious that the initially assumed constant-velocity field will be very different from the final solution. The problem is solved using the CAFFA numerical code on structured grids of size 40x40 cells, and with the FVM code on unstructured grids with approximately 1600 cells. Fig. 3 a shows a part of the unstructured grid for ß= 45° that is used in the FVM code. The borders of the finite volumes are presented, and the main nodes are marked. Fig. 3 b shows a part of the geometric grid for the same case, which is used in the CAFFA code. The displayed lines connect the main nodes at which the pressure and velocity fields are calculated. In the SIMPLE algorithm, two under-relaxation factors should be given. The rate of convergence depends on the values of these two factors. Their optimal values are not known in advance, so that the described problem will be solved for a range of under-relaxation factors by varying a from 0.5 to 0.95 with a step of 0.025, and a from 0.1 to 0.6 with a step of 0.1. The comparison criteria will be the number of iterations needed for the residuals to fall below s= 106. In the CAFFA and FVM codes different solvers for linear algebraic equations are used. For this reason, a sufficient number of inner iterations is given at every iterative step to be sure that the systems are solved equally well in both codes. Fig. 4 shows the numbers of outer iterations N required to reduce the residual levels to s as a function of the under-relaxation factors a and a, 108 Džijan I. - Virag Z. - Kozmar H. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)2, 105-113 a) b) Fig. 2. Streamlines in laminar flow in a lid-driven cavity a) b = 90°, b) b = 45° f/ // / / // / /, J / / / / /////// /////// //////// ///////// //////// j// / //////////// // / / / //////// ////// //////// //////, '///////. v / / / / / //////// JY ////// //////// /z//////////////, ///////////////// // / / ////// //////// // / / / ////// '//////// /////, '///// ////////. ////// ////// //////// ///////////// //////// a) b) Fig. 3. A part of the grid for the lid-driven cavity problem for b = 45° a) FVM, b) CAFFA 1000 800 600 400 200 1000 800 600 400 200 0.5 0.6 0.7 0.8 0.9 0.5 0.6 0.7 0.8 0.9