Department of Mechanical
Engineering, Brunel University, Uxbridge, Middlesex, UB8 3PH, UK.
Email: Vincent.Anjorin@brunel.ac.uk
Keywords: SIMPLE algorithm, Flow solver, Pressurecorrection equation, NavierStokes equations
The SIMPLEV (SIMPLEVincent) algorithm is an improved version of the
standard SIMPLE algorithm where the underrelaxation and temporal terms
are removed from the pressure correction equation of the SIMPLE algrithm.
This paper focuses on the methodologies of the SIMPLEV algorithm and how
it is used to solve the discretized momentum equations that yield the pressure
correction terms of the pressure correction equation. The SIMPLE and SIMPLEV
algorithms rely on underrelaxation to avoid divergence and there is the
problem that by using a very small velocity underrelaxation factor or
a very small time step, the convergence rate of the SIMPLE and SIMPLEV
algorithms becomes extremely slow since a small pressure underrelaxation
factor must be employed. The SIMPLE algorithm is assessed with a version
that has the terms for the pressure correction equation that increases
the pressure correction as the velocity underrelaxation factors or the
time step terms decrease. It is shown that the SIMPLEV algorithm gives
the reverse effect of this.
A method for solving the momentum equations for laminar flow problems has been developed by Patankar and Spalding (1972). They proposed the SIMPLE (SemiImplicit Method for PressureLinked Equations) algorithm to iteratively solve the momentum equations in their discretized form. The methodology of the SIMPLE algorithm has been recently discussed by Versteeg and Malalasekera (1995), explaining how to calculate the velocity field for a twodimensional control volume. The real difficulty in calculating the velocity field lies in determining the unknown pressure field. There is no obvious equation to solve thus pressure is used to satisfy the condition for continuity. To solve this the pressure field is indirectly specified via the continuity equation. This indirect specification is achieved by obtaining a whole set of discretized equations from the momentum and continuity equations and solving the discretized equations by a direct solution.
2.1. Governing Equations
Twodimensional incompressible laminar constantdensity flow (Melaaen, 1993) is governed by a set of partial differential equations. The momentum and continuity equations in their primitive form are shown in equations (13) where the equation for conservation of mass is given by
(1)
The conservation of momentum in the x and y directions are governed by the umomentum equation expressed as
2.2. Discretisation
In order to numerically solve
the velocity and pressure fields that obey the discretized momentum and
continuity equations, the finitevolume method was applied. This method
involves integrating the continuity and momentum equations over a twodimensional
control volume on a staggered differential grid (Patankar
and Spalding, 1972; Harlow and Welch, 1965) shown
in Figure 1.
This yields the governing equations in their discretized form as shown in equations (46). The staggered grid evaluates the scalar variables, in this case only the pressure, which are stored at the scalar nodes marked (·), and located at the intersection of two unbroken grid lines. Such points are indicated by the capital letters P, W, E, N and S. The uvelocity components are stored at the east and west cell faces of the scalar control volume and are indicated by the lower case letters e and w. The vvelocity components are stored at the north and south cell faces of the scalar control volume which are indicated by the lower case letters n and s. These velocity components are located at the intersection of a dashed and unbroken line that construct the scalar cell faces and are indicated by arrows. The horizontal arrows shown in Figure 1 indicate the locations of the u_{e} and u_{w} velocity components and the vertical arrows indicate the locations of the v_{n} and v_{s} velocity components. Forward staggered velocity grids are used. The uniform grids are forward staggered since the uvelocity, u_{e}, is at a distance of 1/2 _{u }from the scalar node P_{p}. Similarly, the location for the vvelocity, v_{n}, is at a distance of 1/2 _{v} from the scalar node. After the process of discretization, the discretized continuity equation becomes:
, (4)
. (8)
The second terms on the righthand side of equations (7) and (8) are the temporal terms where (k1) implies values of the velocity components from the previous time step. The first terms on the righthand side of equations (7) and (8) stand for the constant part of the average value of the source terms. The method of discretizing the momentum and continuity equations is fully explained by Versteeg and Malalasekera (1995) where they show how to employ differencing schemes to successfully interpolate between the node points.
The velocity and pressure fields
of the SIMPLE and SIMPLEV algorithms are all solved on a staggered differential
grid arrangement. This arrangement is used to prevent the pressurevelocity
(PV) coupling that links the mass conservation equation and the momentum
equations. A further discussion of this is presented by Simoneau
and Pollard (1994). In the past, a majority of applications use the
staggered grid arrangement to solve all the flow variables for numerical
modelling purposes. As an example of a typical problem, Stathopoulos
and Baskaran (1996) carried out simulations of the mean wind environmental
conditions around buildings for the assessment of the dispersion of pollutants.
The differential equations for the computation of turbulent wind flow conditions
around buildings were discretized and represented on a nonuniform staggered
grid arrangement containing 235,000 nodes. However more recently, finitevolume
flow calculations have often used the Rhie and Chow interpolation method
(Oliveira et al, 1998; Rhie and
Chow, 1983) which is based on a nonstaggered grid. This approach cures
the problem of oddeven coupling between the pressure and velocity fields
by employing the technique of interpolating the cellface velocities via
momentum interpolation.
3. Solution Procedure of the SIMPLE Algorithm
The standard SIMPLE algorithm (Patankar, 1980) can be broken down into the following steps:
First solve the discretized u and v momentum equations (9) and (10) by using the current guessed pressure field p^{*}to yield the intermediate velocity fields , .
(14)
. (17)
4. Solution Procedure of the SIMPLEV Algorithm
The pressure correction equation of the SIMPLEV algorithm is derived by first determining if there is a connection between the underrelaxation factor (a_{u} and a_{v}) and the timestep . Using the underrelaxed discretized umomentum equation (22), the temporal term is disregarded so that a comparison can be made of the effect of underrelaxation and the original temporal term. Taking the lefthand side of equation (22) and removing the temporal term gives
. (35)
. (37)
. (41)
Equations (40) and (41) are substituted into the continuity equation (4) to give
(42)
. (46)
. (48)
The problem of the flow around a square cylinder was chosen because the geometry is simple. It is an easy problem in which the SIMPLEV algorithm can be applied to investigate how the pressure field on a control volume can be improved in order for the velocity fields to satisfy continuity and hence how the converged solution of the pressurecorrection equation can be determined. An inlet Reynolds number of 800 was used. The Reynolds number for this case is defined using the diameter of the square cylinder, the freestream velocity and the viscosity as mentioned by Barton (1998). A schematic showing the geometry and boundary conditions is shown in Figure 3. The grid used was compressed towards the upper and lower boundaries, as well as the solid boundaries of the square cylinder to better resolve the higher gradients in those regions. The upper and lower boundaries are situated 4D away from the cylinder, where D is known as the diameter of the cylinder. The length scales are nondimensionalised by D. The inlet and outlet boundaries are placed 6D from the cylinder.
The predicted velocity profiles using the SIMPLEV scheme are shown in Figure 4. The uvelocity profiles were predicted at the positions x = 0, 8, 12, 14, 20 & 25 downstream from the inlet boundary. It can be shown that the uvelocity profiles are symmetric about the centreline of the flow field just as the flow approaches the cylinder. Behind the cylinder the Uvelocity profile changes due to the flow separation in front of the cylinder and the formation of recirculation regions behind the cylinder. As the flow moves towards the outlet boundary, the Uvelocity increases further away from the recirculation regions.
Figure 5
shows the flow around a square cylinder where the flow approaches, separates
and is forced away from the square cylinder forming two recirculation regions
behind the cylinder. The upstream flow slowly recovers downstream returning
to its original profile. Using nondimensionalized timesteps at particular
time intervals the recirculation regions increase with respect to time
(Barton, 1998).
5.2. Flow over a BackwardFacing Step
The problem of the flow over a backwardfacing step was chosen because it is fundamental in design and geometry, and is used in a variety of engineering applications. The sudden changes in pressure in which our SIMPLEV algorithm numerically evaluates can be studied using backwardfacing step configurations. The configuration of the flow over a backward facing step is shown in Figure 6 which has an inlet channel of height (h) and a parabolic inlet flow profile.
The channel expansion number was set to 2 and this is defined as the ratio of the total height of the main channel to the height of the inlet channel (h). The inlet Reynolds number was set at Re = 800, where the Reynolds number is based on the total channel height (2h), the average inlet velocity and the dynamic viscosity .
5.2.1. Flow Predictions for the BackwardFacing Step Problem
Shown in Figure
7 is the flow configuration of the backwardfacing step problem where
the flow moves downstream from the inlet channel forming a recirculation
region close to the step (Barton, 1998). It has
a main recirculation point where the flow reattaches. As the flow recovers
downstream, a recirculation region is formed that separates and later recovers
as the pressure recovers. Initially the flow has four recirculation regions.
The two recirculation regions furthest downstream eventually disappear
with increasing time because the vorticity of the flow decreases finally
enabling the flow to have two recirculation regions.
Shown in Table 1 are the predictions
of the growth of the main reattachment and separation positions with time.
The table shows that the SIMPLEV algorithm gives the longest reattachment
and separation positions for iterations expressed in nondimensional time
units.









T

X1

X2

X3

X1

X2

X3


























































5.3. Performance of the SIMPLEV Algorithm
The performance of the SIMPLEV algorithm
was compared with the SIMPLE algorithm to search for the minimum number
of iterations and the optimum underrelaxation factors that enable the
solution to converge. During the computational process, the pressurecorrection
equations of the SIMPLE and SIMPLEV algorithms were iteratively solved
until the tridiagonal matrix algorithm (TDMA) solver (Anderson
et al, 1984) terminates operation. This algorithm solves the matrix
of the pressure correction equations of the SIMPLE and SIMPLEV algorithms.
The convergence criterion that was used for the TDMA solver was 5x10^{5}
for the backwardfacing step problem and 7x10^{4
}for the problem
of laminar flow around a square cylinder. For both algorithms, the convergence
criteria remained the same so that there can be a comparison between the
two schemes. This happens when the velocity and pressure residuals reach
their desired convergence criteria.
In order to obtain convergence, smaller values of underrelaxation factors should be used. However this poses a problem that an increasing number of iterations are required to yield convergence and the solutions are therefore obtained with more CPU time than in the case of the SIMPLE method. An alternative way to consider the algorithms is by introducing a timestep instead of underrelaxation. The results of Figure 9 presents the convergence rates of the SIMPLE and SIMPLEV schemes expressed in nondimensional time intervals. For both schemes it can be shown that the converged solution is easily obtained as the time step decreases. However using very small number of time steps leads to an increasing number of iterations that are required to enable the solution of the pressure correction equation to converge. For the case of laminar flow around a square cylinder of Figure 9 an initial time step of 0.01 units was applied for both schemes in which convergence was obtained after 782 iterations for the SIMPLEV scheme and 812 iterations for the SIMPLE scheme. This was marched in time of 0.01 units until a maximum time step of 0.085 units was employed. For the 60x60 grid system, the SIMPLE and SIMPLEV algorithms require more iterations than the 30x30 grid system to obtain convergence. This is because by having more grid points in the computational domain increases the number of algebraic equations that need to be solved for the pressure and velocity fields. Because both methods are iterative, more equations almost always means slower convergence measured in terms of the total number of iterations required. In addition, each individual TDMA step requires more CPU time, therefore causing slow convergence of both schemes. With more number of node points in the computational domain, the SIMPLEV algorithm is therefore not able to improve the efficiency of computations. The computational tests for the two grid systems of Figure 9 were performed in order to investigate the influence of the number of nodes on the SIMPLE and SIMPLEV algorithms. Above all the SIMPLEV scheme yields better convergence as fewer iterations are required for the converge solution of the pressure correction equation.
Figure 10
shows the results of the number of iterations of the iterative versions
of the SIMPLE and SIMPLEV algorithms as well as the changes of the residual
values of the velocities and pressures during the iterative computation
of the AFLOW code. It can be shown that the SIMPLEV algorithm attains the
required value of the convergence criterion in less time than the SIMPLE
algorithm. It should be noted that the results shown in Figure
10 were obtained for the 30x30 and 60x60 grid systems. When the SIMPLEV
algorithm is employed, there is a reduction in the residual values for
every iteration and the convergence criterion in all cases is attained
fastest than the SIMPLE method. However it can be observed from Figure
10 that the SIMPLEV algorithm is not able to improve the efficiency
of computations if a grid system with large number of nodes is employed.
This makes it quite difficult for the converging solutions to be obtained.
The set of optimal underrelaxation coefficients for which there is higher
computational efficiency is shown in Table 2. This set was determined for
the 30x30 grid system. A comparison of the minimal number of iterations
and the set of optimal time steps for the SIMPLE and SIMPLEV methods for
which there is higher computational efficiency is shown in Table 3. They
indicate the necessary number of time steps for the solution of the computational
process to reach convergence.
A special thanks to Mr Callum Downie for his help in setting up the AFLOW code on the UNIX workstations.
1. Alvaro Valencia., (1995) , Heat Transfer Enhancement in a Channel with a Builtin Square Cylinder?, International Communications in Heat and Mass Transfer, 22, Issue 1, 4758.
2. Anderson, D. A. Tannehill, J. C. & Pletcher, R. H., (1984), Computational Fluid Mechanics and Heat Transfer, Hemisphere Publishing Corporation, Taylor & Francis Group, New York.
3. Barton, I. E., (1995), A Numerical Investigation of Incompressible Dilute Particulate Flows, Ph.D. Thesis, Aerospace Division, School of Engineering, University of Manchester.
4. Barton, I. E., (1998), Comparison of SIMPLE and PISOtype Algorithms for Transient Flows, International Journal for Numerical Methods in Fluids., 26, 459483.
5. Barton, I. E., (1998), Improved Laminar Predictions using a Stabilised TimeDependent SIMPLE Scheme, Int. J. Numer. Meth. Fluids., 28, 841857.
6. Barton, I. E. & Kirby, R., (1998), FiniteDifference Scheme for the Solution of Fluid Flow Problems on NonStaggered Grids, Int. J. Numer. Meth. Fluids., 33, 939959.
7. Breuer, M. Bernsdorf, J. Zeiser, T. and Durst F., (2000), Accurate Computations of the Laminar Flow Past a Square Cylinder Based on Two Different Methods: LatticeBoltzmann and FiniteVolume, International Journal of Heat and Fluid Flow, 21, 186196.
8. Fletcher, C. A. J., (1988), Computational Techniques for Fluid Dynamics. Vol. II, Springer, New York.
9. Guo T. Chew Y.T. Luo S.C. and Su M.D., (1998), A New Numerical Simulation Method of High Reynolds Number Flow around a Cylinder, Computer Methods in Applied Mechanics and Engineering, 158, Issue 34, 357366.
10. Harlow, F. H. and Welch, J. E., (1965), Numerical Calculations of TimeDependent Viscous Incompressible Flow of a Fluid with a Free Surface, Physics of Fluids., 8, 21822189.
11. Huang, P. G. and Launder, B. E. and Leschziner, M .A., (1985), Discretization of NonLinear Convection Processes: A Broad Range Comparison of Four Schemes, Comput. Methods Appl. Mech. Eng., 48, 124.
12. Iwai, H. Nakabe, K. Suzuki, K., (2000), Flow and Heat Transfer Characteristics of BackwardFacing Step Laminar Flow in a Rectangular Duct, International Journal of Heat and Mass Transfer, 43, Issue 3, 457471.
13. Jang, D. S. Jetli, R. and Acharya, S., (1986), Comparison of the PISO, SIMPLER, and SIMPLEC Algorithms for the Treatment of the PressureVelocity Coupling in Steady Flow Problems, Numer. Heat Transfer, 19, 209228.
14. Johansson, P. and Davidson, L., (1997), Full Multigrid Method Applied to Turbulent Flow using the SIMPLEC Algorithm together with a Collocated Arrangement?, Doktorsavhandlingar vid Chalmers Tekniska Hogskola.
15. Jun Zhang., (1996), Acceleration of FivePoint RedBlack GaussSeidel in Multigrid for Poisson Equation, Applied Mathematics and Computation, 80, Issue 1, 7393.
16. Keskar, J. Lyn, D.A., (1999), Computations of a Laminar BackwardFacing Step Flow at Re = 800 with a Spectral Domain Decomposition Method, International Journal for Numerical Methods in Fluids, 29, Issue 4, 411427.
17. Leonard, B. P., (1979), A Stable and Accurate Convection Modelling Procedure Based on Quadratic Upstream Interpolation, Comput.Methods Appl. Mech. Eng., 19, 5998.
18. Melaaen, M. C., (1993), Nonstaggered Calculation of Laminar and Turbulent Flows using Nonorthogonal Coordinates, Numer. Heat Transfer, Part A., 24, 375392.
19. Oliveira P.J. Pinho F.T. Pinto G.A., (1998), Numerical Simulation of NonLinear Elastic Flows with a General Collocated FiniteVolume Method, Journal of NonNewtonian Fluid Mechanics, 79, Issue 1, 143.
20. Patankar, S. V. and Spalding, D. B., (1972), A Calculation Procedure for Heat, Mass and Momentum Transfer in Threedimensional Parabolic Flows. International Journal for Heat and Mass Transfer., 15, 17871806.
21. Patankar, S. V., (1980), Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, Taylor & Francis Group, New York.
22. Ranade, V. V., (1997), An Efficient Computational Model for Simulating Flow in Stirred Vessels: A Case of the Rushton Turbine, Chemical Engineering Science, 52, Issue 24, 44734484.
23. Rhie, C. M. And Chow, W. L., (1983), Numerical Study of the Turbulent Flow past an Airfoil with Trailing Edge Separation?, AIAA Journal, 21, 15251532.
24. Robert R. and Hwang., (2000), Ground Effect on Vortical Flows behind a Square Cylinder, Proceedings of the International Offshore and Polar Engineering Conference, 3, 467470.
25. Simoneau, J. and Pollard, A., (1994), FiniteVolume Methods for Laminar and Turbulent Flow using a Penalty Function Approach, International Journal for Nunerical Methods in Fluids., 18,
26. Spalding, D. B., (1972), A Novel Finitedifference Formulation for Differential Expressions Involving Both First and Second Derivatives, Int. J. Numer. Methods. Eng., 4, 551.
27. Stathopoulos,T. Baskaran, B. A.,(1996), Computer Simulation of Wind Environmental Conditions around Buildings, Engineering Structures, 18, Issue 11, 876885.
28. Versteeg, H. K and Malalasekera, W., (1995), An Introduction to Computational Fluid Dynamics: The Finite Volume Method, McGrawHill, Loughborough University.
29. Wanik, A. & Schnell, U., (1989), Some Remarks on the PISO and SIMPLE Algorithms for Steady Turbulent Flow Problems, Computers & Fluids., 17, No.4, 555570.
30. Wengle, H. Huppertz, A. Barwolff, G. Janke, G., (2001), ?Manipulated Transitional BackwardFacing Step Flow: An Experimental and Direct Numerical Simulation Investigation, European Journal of Mechanics, B/Fluids, 20, Issue 1, 2546.
31. Zaiyong, S. Hongqing, H. and Timin, C., (1998), Solving Viscous Flow in Nozzles using Improved SIMPLE Algorithm, Tuijin Jishu/Journal of Propulsion Technology, 19, Issue 3, 4851.