In previous blogs in this series I have recounted my recollections of the first development of 3D CFD methods at Imperial College in the very early 1970s - and how, after a “false start” with the SIVA algorithm, I adopted Suhas Patankar’s SIMPLE algorithm (with a refinement of my own) for my PhD work.
As will have been apparent from the earlier blogs, the main focus of the 3D Boundary Layer Group, of which I was part, had been on developing the numerical methods and algorithms themselves - leading to the successful application of SIMPLE to a range of laminar flow situations, including my own work on developing flows in square and rectangular sectioned ducts. But of course, as we all know, this only takes us part way towards a practical CFD solution for real engineering flows - there is still the troublesome matter of turbulence modelling!
During the years I have referred to in these blogs (essentially 1965 to 1972), as well as the groundbreaking work on solution methods that I have described, equally important work was going on by others in Spalding’s group at Imperial College on turbulence modelling, and on new measurement techniques for the required supporting experimental work.
These complementary strands of work were led, respectively, by Brian Launder and Jim Whitelaw - both relatively junior lecturers when I first knew them. Launder, and Spalding himself, with a series of talented research students (Wolfgang Rodi, Bill Jones, Kemo Hanjelic, and others), laid the foundations of practical turbulence modelling during this period, culminating in the emergence of the two-equation k-epsilon model in the form still used today - and Whitelaw, initially with Franz Durst and then other research students, developed the Laser Doppler Anemometer for measuring turbulent flows into the practical methodology that is widely used today.
I was myself only an interested observer of these activities, taking place literally “just down the corridor” - so I won’t attempt any kind of detailed account. Suffice it to say that, looking back on all this, it does seem remarkable that so much progress was made in these three fields (numerical methods, turbulence modelling, and experimental methods) by this small group of people in such a short time. I think that it just shows what a hotbed of talent Spalding had assembled at Imperial College - and how, within a small team, success can be infectious. Achievement in one area can inspire success and further achievement elsewhere within the team - and so on - a sort of positive feedback effect. That was certainly the way it felt at the time!
Turning to my own work - when I had completed my laminar flow computations, the next step was to move on turbulent flows - initially developing flow in square ducts. By this time (I guess during 1971) the two-equation k-epsilon model was already established as the “standard” modelling approach for turbulent flows. The model was originally devised by Harlow’s group at Los Alamos (first published in 1965) - but the main development and refinement into a practical engineering tool had taken place at Imperial College. As well as being used (for 2D flows) within Spalding’s group, it was beginning to be publicised more widely, via a lecture course which led to the (at that time) standard text on the subject - Launder and Spalding, “Mathematical Models of Turbulence”, Academic Press, 1972). So this was naturally the model that I adopted.
However, this didn’t quite do the job. For turbulent flow in a square duct (actually for any non-circular sectioned duct), some strange things happen far downstream…
Consider the limiting case of fully-developed flow. For fully-developed laminar flow the axial flow establishes a fixed velocity profile (parabolic in a circular duct - other profiles for other shapes), with no lateral motion at all. It is fairly easy to see that this satisfies the mass continuity equation. In the axial direction, by definition (of fully-developed flow) the gradient in axial velocity is zero. This is one of the three terms in the continuity equation, The others involve lateral gradients of the lateral velocities - and if these velocities are zero, then the gradients will be zero - so the continuity equation is satisfied. And, for laminar flows, in the absence of any forces acting on the fluid in the lateral directions (buoyancy, rotation, electromagnetic, etc), then this is the observed physical behaviour - zero lateral velocities - whatever the duct cross sectional shape.
However, for turbulent flows the behaviour is different, Even in the absence of any forces acting on the fluid, for any non-circular duct the fully-developed flow exhibits a pattern of lateral motions. For a square dust you get essentially eight triangular vortices, with flow from the duct centre to the corner along the diagonals, along the walls from corner to centre planes, and then back out to the duct centre along the centre planes (as illustrated in the figure below). The velocities involved are small (a few percent of the axial flow) but they can have an appreciable effect (of order 10%) on the frictional pressure drop and heat transfer - so they are of some practical engineering significance.
Fully Developed Flow in a Square Duct - Predictions at Re 83,000. Left hand quadrant shows axial velocity contours and lateral velocities - right hand quadrant shows normalised turbulent kinetic energy and length scale.
The reason that I’ve referred to this somewhat esoteric phenomenon, is to make a point about the status of turbulence modelling knowledge in 1971. The point is that it was clear there is no way that the k-epsilon model on its own was going to predict the existence of these lateral motions. Indeed, no model based on the “turbulent viscosity” approach would predict them. It can be shown that, whatever turbulent viscosity distribution one might presume, there is no way that a variation of turbulent viscosity alone can lead to the initiation of lateral motions of the kind. In the absence of source terms of some kind in the lateral momentum equation, the only possible (wrong) fully-developed solution is zero lateral velocities.
However, in parallel with the creation of the k-epsilon model, the work by Spalding and Launder (and particularly Launder with his student Kemo Hanjelic) had led to the development of a full Reynolds stress model - in which, rather than representing the turbulent Reynolds stresses via the turbulent viscosity, the stresses (six in all for a 3D flow) are solved for directly, via their own differential transport equations. At this stage it was impracticable to contemplate solving a full 3D flow using such models (the available computer resources were insufficient), but the models had been used in reduced form to address some particular turbulence phenomena - one of which was these lateral motions for fully developed flow in a square duct.
Launder and Ying had solved the 2D situation (ie only perfectly fully developed flow) for flow in a square duct using a hybrid model, in which the k-epsilon model was used in conjunction with an algebraic form of the Reynolds stress model for the stresses affecting the lateral velocities. Here “algebraic” means that, by making presumptions about the convection and diffusion terms, the Reynolds stress equations were algebraic rather than differential - which made solution more practicable. The result of this is that the normal and shear Reynolds stress terms in the lateral momentum equations appear as source terms, which (Ying showed) led to the prediction reasonably well of these odd lateral motions.
So this was the model that I used - what I was effectively doing was (among other things) investigating the development of these lateral motions up to the fully developed regime that Ying had considered. And what I was able to show is that the model did indeed predict the gradual development of these lateral motions as the flow developed - but that it was only pretty far downstream (more than 30 duct widths) that they became really established and had a real influence on axial velocity, and hence on friction or heat transfer. Which is, I believe, why these lateral motions are really only of academic interest. It is only for situations involving perfectly straight ducts of 30-plus duct widths in length, with no other disturbances likely to create lateral motions themselves, for which these turbulence-induced motions are significant - and such cases are relatively rare in practical engineering situations. But they are an interesting - and challenging - test of turbulence modelling.
Looking back on this work, two things stand out for me. Firstly, it is remarkable that by 1971 the basic k-epsilon model had already reached the form in which it is known and widely used today - and was already a “standard model“ ready for adoption by “users“ such as myself. And, secondly, it is perhaps even more remarkable that, for the “special cases” for which the k-epsilon model fails, there was already a fully-formulated Reynolds-stress modelling approach available to be applied.
My square duct calculations may have been the very first 3D application of the k-epsilon model anywhere - certainly all the prior work I referred to in my thesis is 2D - but if anyone has a counter-claim that’s fine by me. What is more important to me is to recall being a very early “early adopter” of the outcome of the pioneering work on understanding and modelling turbulence by Spalding, Launder and Whitelaw during that remarkable period at the end of the 1960s and early 1970s!