# Beginning at the Beginning … 2

In the first blog in this series I described how the research group at Imperial College, under Professor Spalding, on computational methods for 3D parabolic boundary-layer flows (the “3D Boundary Layer Group”) got started, with me as a relatively junior member. I described how I (re)discovered the need to use a staggered grid when solving for the primitive variables, (pressure and velocities) – and how Suhas Patankar joined the group around mid 1970.

To continue …

To start with I need to get rather technical – anyone who wants to avoid the “technical interlude” can skip the next seven paragraphs.

There is a fundamental problem in any fluid-dynamic computational method that attempts to solve directly for the “primitive variables”, pressure and velocities. That is – how to compute the pressure field?

Let me try to explain. First consider the variables that we need to solve for, and the equations from which they are to be deduced. For 3D flows (with no heat transfer or other processes), to completely define the flow field we need to compute the three velocity components (in x, y and z) and the pressure field – four variables in all. And there are four fundamental equations governing the flow field – mass continuity, and conservation of momentum in each of the x, y, and z directions. So, four variables and four equations – so far so good!

However, then we try to decide which variable to compute from which equation. The x, y and z velocities (u, v, and w, say) are, respectively, the natural dependent variables of the x, y and z momentum equations. This leaves pressure – which we therefore need to compute from the continuity equation. However, for incompressible flows (with which we were concerned) pressure does not appear explicitly in the continuity equation (which is simply a balance of mass inflows and outflows for a computational grid cell). The choice about how to resolve this dilemma is at the heart of any CFD solution method.

In the 3D Boundary Layer Group, we started out using a new computational method known as SIVA (standing for “simultaneous variable adjustment”), devised by Professor Spalding. This addressed the “pressure problem” as follows. Consider (in 2D) the grid cell arrangement shown in the figure below. The continuity equation for the grid cell shown involves the velocities ue, uw, vn, vs. (and also those in the third direction – but for simplicity forget about these for now). In the continuity equation, it is possible to eliminate these velocities, using the finite volume momentum equations, by:- replacing  ue by all the terms on the right-hand side of the finite volume momentum equation “ue = ….”; and so on for uw, vn, vs, and for the third direction velocities.

The SIVA Grid Cluster:- Pressures (and other variables) are computed at grid nodes P, N etc - velocities (u and v) at the staggered locations shown. The grid cell for node P is shown dashed.

This leads to a rather complex equation, in which attention is focussed on the terms involving pressure (arising from the pressure gradient terms in the momentum equations). If the pressure terms are taken to the left hand side, with all the other terms on the right hand side, then what results is the finite volume form of a Poisson Equation for pressure – involving pP, pN, pS, pE and pW.

The SIVA method involved solving this equation, and the momentum equations, at each “cluster of cells”, so as simultaneously to satisfy the (linearised) momentum equations at n, s, e and w (and, in 3D, the two z-direction momentum equations), and the continuity equation at P. Of course, because the momentum equations are linearised, and because changes in values of variables outside the cluster are necessarily neglected, the usual iteration is required to achieve a converged solution.

If this description sounds complex, then it perhaps illustrates the point that SIVA was a complex algorithm to program. It requires the continuity equation to contain a perfect, term-by-term representation of (in 3D) six momentum equations. Any inconsistencies between the continuity equation and the momentum equations as solved for the velocity fields would lead to non-convergence (the two versions would be “pulling against each other”, which would lead to an oscillatory, non-converged solution).

End of “technical interlude” …

So the SIVA solution method that the 3D Boundary Layer Group was developing and using was ingenious – but complex. I was able to get it working for my target application – laminar flow in the inlet region of a rectangular-sectioned passage (as reported in a paper cited later) – and others were applying it to other flow situations. However, from my experiences of working with the method, I was beginning to appreciate the challenges in extending it to more complex situations.

Then Suhas Patankar returned from India to join our research group. For the first few weeks we didn’t see much of him – but he was evidently working hard! He announced the results at one of the weekly progress meetings of the research group. I have a vivid recollection of the occasion – I was tasked with keeping the records of the meetings, so I had to pay particular attention!

Essentially the meeting focussed on the conclusions that Suhas had reached from his few weeks deliberation on a better way to approach the solution of the pressure/velocity equations. And what Suhas presented was the SIMPLE algorithm. I recall Suhas, in his usual quiet way, giving a clear, detailed description of SIMPLE pretty much “fully formed” – exactly as it was later published, and subsequently widely used.

As CFD practitioners will know, the essence of SIMPLE is making pressure adjustments (or corrections) which, at any stage in the solution, are approximate (ie there is no attempt to mimic the full momentum equations, as in SIVA), but which will progressively move the pressure field in the direction of the “correct solution” – that is, a set of pressures which, when the momentum equations are solved, yield velocities that satisfy continuity at every grid cell. When this solution is reached, the pressure corrections become zero, and the solution is converged.

The real beauty of this is that it removes the complexity and the rigidity of methods such as SIVA that attempt to get the pressures “right” at all stages in the solution. It means that approximations and damping can be introduced relatively freely in the pressure correction equation – all that is needed is to ensure that the corrections are always changing the pressures in the right direction by roughly the right amount. This flexibility has, in my view, been crucial in making it possible to apply SIMPLE to the vast range of flow configurations that is has been used for over the last 40 years, and in enabling developers to devise the various extensions, variants, and improvements that have emerged over the years.

I recall that, as I listened to Suhas’s presentation, I could begin to see some of these advantages, particularly from my experiences of working with the SIVA algorithm. My immediate reaction was “what a fantastic idea – why didn’t I think of that?”. But of course I didn’t! Like many great ideas – it is blindingly obvious once it has been explained – but impossible for mere mortals to conceive without prior explanation!

To complete the story very briefly – the research group pretty much adopted SIMPLE straight away, and it was what I used for the later turbulent flow calculations in my PhD. And it is of course what, in one form or another, the majority of the CFD community has used ever since!

I understand that Suhas was at that time a talented amateur magician (he probably still is). It always seems to me that there was a touch of magic in the way that Suhas devised SIMPLE in those few weeks back in 1970! It was privilege to be present at what was perhaps the most crucial breakthrough in the development of CFD (so far at least!).

I was able to make my own modest contribution to improving the SIMPLE algorithm a little later – but this is for the next instalment …

Postscript:- For anyone who is interested, two papers were published by the 3D Boundary Layer Group describing SIVA and the early applications:
Caretto, Curr and Spalding (1972) “Two numerical methods for three-dimensional boundary layers”, Computer Methods in Applied Mathematics and Engineering, Vol 1, pp 39-57.
Curr, Sharma, and Tatchell (1972) “Numerical predictions of some three-dimensional boundary layers in ducts”, Computer Methods in Applied Mathematics and Engineering, Vol 1, pp 143-158.
Plus, there is of course the classic original paper describing SIMPLE:
Patankar and Spalding (1972) “A calculation procedure for heat, mass, and momentum transfer in three-dimensional parabolic flows”, Int J Ht Mass Transfer, Vol 15, pp 1787-1806.
And finally, of more general interest, there is a fascinating write up of Brian Spalding’s career and contributions in a 2008 paper by Akshai Runchal, “Brian Spalding: CFD and Reality”.

## Post Author

Posted August 13th, 2009, by

1 Comment