International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer
journal homepage: www.elsevier.com/locate/ijhmt

A tribute to D.B. Spalding and his contributions in science and engineering
V. Artemov a, S.B. Beale b,*, G. de Vahl Davis c, M.P. Escudier d, N. Fueyo e, B.E. Launder f, E. Leonardi c,1,
M.R. Malin g, W.J. Minkowycz h, S.V. Patankar i, A. Pollard j, W. Rodi k, A. Runchal l, S.P. Vanka maMoscow Power Engineering Institute, Moscow, Russia
National Research Council of Canada, Montreal Road, Ottawa, ON, Canada K1A 0R6c
The University of New South Wales, School of Mechanical and Manufacturing Engineering, Sydney, NSW 2052, AustraliadUniversity of Liverpool, Department of Engineering, Liverpool L69 7ZE, UKe
University of Zaragoza, Numerical Fluid Dynamics Group, 50018 Zaragoza, Spainf
The University of Manchester, School of Mechanical, Aerospace and Civil Engineering, P.O. Box 88, Manchester M60 1QD, UK
g
Concentration Heat and Momentum Limited, 40 High St., Wimbledon Village, London SW19 5AU, UKhUniversity of Illinois at Chicago, Department of Mechanical and Industrial Engineering, Chicago, IL 60607, USAi
University of Minnesota, Department of Mechanical Engineering, Minneapolis, MN 55455-0111, USAj
Queen’s University, Department of Mechanical and Materials Engineering, Kingston, ON, Canada K7L 3N6
k
University of Karlsruhe, Kaiserstr.12, D-76128 Karlsruhe, GermanylAnalytic & Computational Research, Inc., 1931 Stradella Road, Bel Air, CA 90077, USAm
University of Illinois at Urbana-Champaign, Department of Mechanical Science and Engineering, 1206 W Green Street, Urbana, IL 61801, USAba r t i c l e

i n f oArticle history:
Received 11 March 2009
Accepted 23 March 2009
Available online 20 May 2009
Keywords:
D.B. Spalding
Fluid dynamics
Heat transfer
Mass transfer
Combustiona b s t r a c t
This paper presents a summary of some of the scientiﬁc and engineering contributions of Prof. D.B. Spalding up to the present time. Starting from early work on combustion, and his unique work in mass transfer theory, Spalding’s unpublished ‘‘uniﬁed theory” is described brieﬂy. Subsequent to this, developments
in algorithms by the Imperial College group led to the birth of modern computational ﬂuid dynamics,
including the well-known SIMPLE algorithm. Developments in combustion, multi-phase ﬂow and turbulence modelling are also described. Finally, a number of academic and industrial applications of computational ﬂuid dynamics and heat transfer applications considered in subsequent years are mentioned.
Crown Copyright Ó 2009 Published by Elsevier Ltd. All rights reserved.

1. Introduction
In May 2008, an International Symposium (CHT’08: Advances in
Computational Heat Transfer) was held in Marrakech, Morocco under the auspices of the International Centre for Heat and Mass
Transfer in honour of Prof. Brian Spalding’s 85th birthday. The conference included a keynote lecture by Prof. Spalding himself andsome lectures describing his research contributions. At the conference, it was decided that a special issue of this journal would be
published in Prof. Spalding’s honour. Whereas this issue will haveinvited papers on a variety of topics, it seemed appropriate to include also an extensive review of Spalding’s research contributions.

* Corresponding author.
E-mail addresses: artemov@itf.mpei.ac.ru (V. Artemov), steven.beale@nrccnrc.gc.ca (S.B. Beale), g.devahldavis@unsw.edu.au (G. de Vahl Davis), escudier@liv.ac.uk (M.P. Escudier), norberto.fueyo@unizar.es (N. Fueyo), brian.launder@manchester.ac.uk (B.E. Launder), mrm@cham.co.uk (M.R. Malin), wjm@uic.edu (W.J.
Minkowycz), patankar@inres.com (S.V. Patankar), pollard@me.queensu.ca (A.
Pollard), rodi@uka.de (W. Rodi), runchal@acricfd.com (A. Runchal), spvanka@uiuc.edu (S.P. Vanka).
1
Eddie Leonardi passed away on 14 December 2008.

With this thought, the authors of this paper (most of whom haveclosely worked with Prof. Spalding) had a meeting at the conference to plan the nature and the structure of this paper. The discussion at that meeting and several months of work thereafter haveculminated in this tribute.
Here, the authors have not just given a factual summary of Spalding’s publications. They have tried to take the reader on a historical journey to show the development of his key ideas and theirsubsequent impact. The authors have shown how the theme ofgeneralization or uniﬁcation runs through all his work. It is thisgrand vision that enabled Spalding to make so many break-throughcontributions. In this paper, the reader will get an insight into thisgreat man’s way of thinking and his extraordinary impact on science and engineering. The authors have all greatly enjoyed theactivity of writing this paper and, in the process, re-lived the experience of working with Prof. Spalding.
2. Nomenclature
In a paper of this nature, the usual Nomenclature section is not
required. Spalding’s passions for science and poetry are well0017-9310/$ - see front matter Crown Copyright Ó 2009 Published by Elsevier Ltd. All rights reserved.
doi:10.1016/j.ijheatmasstransfer.2009.03.038

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

known. It is perhaps the interplay between these seemingly disparate pursuits that has given him the talent to invent wonderfully
descriptive acronyms that describe many of his scientiﬁc contributions. The list is long; only a few are provided here.
PHOENICSCHAMCORA
ESCIMO
FLASH
FLIRT
IPSA
MOSIE
PICALOSIMPLETHIRBLEaParabolic, Hyperbolic Or Elliptic Numerical
Integration Code Series
Originally Combustion, Heat and Momentum Ltd.;
later (since 1974) Concentration, Heat andMomentum Ltd.
COmbustion and Radiation Analyser
Engulfment, Stretching, Coherence, Inter-diffusion,and Moving Observer
FLow Around Ships’ Hulls
FLow In Right angled Tee junctions
InterPhase Slip Algorithma
Movement Of Smoke In Enclosures
PIston and Cylinder CALculatOr
Semi-Implicit Method for Pressure-Linked Equations
Transfer of Heat In Rivers, Bays, Lakes and Estuaries

Sometimes also referred to as InterPhase Slip Analyser.

ously, and within the separate disciplines of either heat, mass or
momentum transfer, varying approaches had been used to formulate the governing equations. As a result he showed that generalized transport equations could be applied to essentially all thestate-variables in ﬂuid mechanics and heat/mass transfer, and indeed to ‘synthesized’ variables such as mixture fraction. This generalization of the equations of motion was pre-requisite to thesubsequent development of generalized solution methodologies
to engineering mechanics problems. Spalding’s mass transfer formulation was ﬁrmly based on mass fractions, as opposed to molar-based balances, it being reasoned that mass (not moles) was
conserved in chemical reactions, and moreover mass-weighted
mixture velocities are typically obtained in solutions of the Navier–Stokes equations. Indeed, the general form of the transport
equation he derived became the natural equation into which allothers were forced and extra terms that did not ﬁt into the substantive derivative, diffusion, production or dissipation terms werelumped as sources (or sinks).
In 1960, Spalding published a pioneering article [10] in vol. 1 ofthis journal, in which he postulated the ‘Ohm’s law’ of mass transfer, namely_ 00 ¼ gB
mOthers are sprinkled throughout this paper.
3. Early work
The origins of Spalding’s life-long contributions to heat andmass transfer and to combustion go back to his Ph.D. Thesis atCambridge University [1]. This thesis introduced the elements thatwere later to govern his professional career. His goal was to understand the combustion of liquid fuels. Spalding came to the conclusion that aerodynamic ﬂow, heat transfer and mass transfer all play
key roles in a synergistic manner. To understand these phenomena,
Spalding then used key hydrodynamic concepts from von Kármán
[2], heat transfer concepts from Kruzhilin [3], and mass transfer
concepts from Eckert and Lieblein [4]. Of course the interactionof chemical kinetics and ﬂow processes plays a key role in combustion. Spalding also beneﬁted from the work of Zeldovich and FrankKamenetskii [5,6], and of Semenov [7].
Subsequently, he went on to make the then unforeseen prediction that the chemical-reaction-rate constants had no inﬂuence on
combustion until a critical rate of mass transfer was reached. He
deduced these critical rates from his theory that were borne out
by experiments. He also developed a methodology to predict ﬂame
extinction from this theory. This work was later expanded andpublished in book form, Spalding [8].
Until the early 1960s Spalding continued to work on combustion and burning of droplets and made a number of innovative contributions. These included the concept of the mass transfer driving
force B, deﬁned below, also known as the Spalding number. These
efforts culminated in his remarkable book on convective heat andmass transfer [9]. Moreover, during this time, he also expanded his
efforts to include general engineering ﬂows in which wall shear
plays an important role. This would lead to his outstanding insights
and the proposal of what came to be known as ‘‘The Uniﬁed Theory”. Prior to the development of ﬁnite-volume schemes, solutions
to problems in convective heat and mass transfer involved thereplacement of diffusion terms with rate equations, based on thenotion of a heat/mass transfer coefﬁcient, which Spalding denoted
by g. The ability to obtain mathematical solutions to the resulting
ordinary differential equations was thus central to the science ofheat and mass transfer. Spalding’s major contribution to masstransfer was to unify (generalize) the problem, so that essentiallythe same systematic approach could be applied across a range of
disciplines from cooling towers to ablation to combustion. Previ-

3885

_ 00

ð1Þ

where m is the mass ﬂow (‘current’), g is a ‘conductance’ or masstransfer coefﬁcient which depends only on aerodynamic and notthermodynamic factors, and B is a mass-transfer driving force (potential) given by B = (PG  PS)/(PS  PT). PG and PS are the values ofsome conserved quantity, P, in the bulk and at the wall surface ofthe ﬂuid, respectively, and PT represents the ‘transferred substancestate’ or ‘T-state’. In later years the symbol ‘P’ was replaced by ‘u’.
For simple mass transfer 0 6 PT 6 1, whereas if chemical reactions
occur 1 6 PT 6 1. Eq. (1) may be applied to enthalpy andmomentum as well as chemical species. A variety of laminar andturbulent mass transfer problems, with or without heat transfer(including radiation), phase change, and chemical reactions wereconsidered in Spalding’s book on convective mass transfer [9],which was published 3 years later. The methodology was subsequently adopted by many others, and remains a standard approachin the engineering curriculum to this date; see, for example, [11,12].
An important point is that the methodology is appropriate for all
convection–diffusion systems, i.e., a general system of conserved
properties, and combinations of such variables.
Around 1960, Spalding had a simple but far-reaching idea [13].
Superﬁcially this may have seemed a relatively unimportant contribution, yet it actually underlined his ability to apply a fresh
direction of thinking to everything he touched. From the 1930s, research workers had been attempting to provide a continuous formula for the variation of mean velocity in the vicinity of asmooth surface where, for simple shear ﬂows, u+ = f(y+). Immediately next to the wall, in the viscous sublayer, the variation was assuredly linear and in the fully turbulent region for y+ > 30 alogarithmic variation was accepted. But what of the region in between, sometimes called the buffer layer? The usual route at thetime was to employ piecewise ﬁts to different segments of the region. It was not easy to come up with a single analytical expression
for the whole universal velocity proﬁle. Spalding’s contribution
was to express the dimensionless distance y+ as a function of thedimensional velocity u+ (instead of the usual practice of writingu+ in terms of y+). The logarithmic part could be easily inverted
as y+ = (1/E)exp(ju+), where E and j are constants. To complete
the expression, he added to the right-hand side the laminar termu+ and then subtracted a few terms from the power-series expansion of the exponential function. The resulting expression did havethe correct asymptotic behaviour (linear and logarithmic) at verysmall and very large values of u+ and did beautifully represent
the transition layer as well.

3886

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

A useful outcome of the above-mentioned analytical expression
was that we could obtain an analytical expression for dy+/du+,which could then be easily related to the dimensionless turbulent(or effective) viscosity. Thus,"

leff dyþ
jðjuþ Þ2 ðjuþ Þ3¼ þ ¼1þexpðjuþ Þ  1  juþ 
2
6
l du
E

#
ð2Þ

where the ﬁrst term on the right-hand side represents the laminar
contribution, while the rest is due to turbulence. Once the effective
viscosity was expressed this way, it was possible to extend the idea
to get an expression for the effective thermal conductivity

Ceff

l

"

¼

leff =l 1 1 jðjuþ Þ2 ðjuþ Þ3¼ þexpðjuþ Þ  1  juþ reff2
6
r rt E#

ð3Þ
where r and rt are the laminar Prandtl number and the turbulent
Prandtl number, respectively, and reff is the effective Prandtl
number.
Using this framework, Spalding proposed that a thermal boundary layer growing through an already established velocity boundary layer (in a turbulent ﬂow) can be calculated downstream of adiscontinuity in wall temperature. This was published in Ref.
[14], where he initially set the turbulent Prandtl number to unity.
The key idea in this paper was used by Kestin and Persen [15],Gardner and Kestin [16], and Smith and Shah [17] to produce
extensive numerical solutions for a range of Prandtl numbers andfor both wall-temperature and heat-ﬂux boundary conditions.
In a later paper [18], Spalding examined these solutions, discovered a regularity in them, provided approximate analytical formulas for them, and introduced the idea of a non-unity turbulent
Prandtl number. He suggested that rt = 0.9 was more appropriate.
This paper made another signiﬁcant contribution. Since the inﬂuence of the laminar Prandtl number is limited to the region near
the wall, the solutions for different values of the laminar Prandtl
number behaved identically once the thermal boundary layer became thicker than the laminar sublayer. Thus, Spalding proposedthat the only effect of the laminar Prandtl number was in providingan ‘‘extra resistance of the laminar sublayer”. This resistance resulted from the fact that the value of the laminar Prandtl number
was different from the value of the turbulent Prandtl number.
In terms of the relevant dimensionless variables, it was possible
to write dT+/du+ = reff. Therefore,T þ ¼ rt

Z 1
0



reff
þ
 1 durt

ð4ÞThis expression was considered to be valid for u+ values beyond thelaminar sublayer, where the effect of the laminar Prandtl number
has all but vanished. This consideration allows the upper limit of
the integral to be set to 1, since the integrand falls to zero for large
values of u+.
Spalding identiﬁed this integral as the resistance of the laminar
sublayer and called it the P function, which depends only on the ratio r/rt. Jayatilleke [19] working with Spalding, developed a comprehensive set of formulae for P functions for smooth and roughwalls in turbulent ﬂows. The most commonly used of these expressions is

P ¼ 9:24

 3=4rrt!
 1 ð1 þ 0:28 expð0:007r=rt ÞÞð5Þ

For moderate Prandtl numbers, this expression is sometimes simpliﬁed to: P ¼ 9:0ðr=rt  1Þðr=rt Þ0:25 . This function is widely used
to this day to provide the wall boundary condition in turbulent ﬂow
calculations.

Fig. 1. Illustration of the type of process which was considered. From Ref. [22] withpermission.

4. The Uniﬁed Theory
In his book Convective Mass Transfer [9], Spalding sought tounify heat transfer and mass transfer. He also considered simultaneous heat and mass transfer (including evaporation, ablation, andcombustion) in a single framework. Around the same time, he became fascinated by the Russian book by Kutateladze and Leontev,
which he translated [20].
In the Preface to the translation, he wrote:‘‘No previous textbook authors have been bold enough to recommend methods for calculating drag and heat transfer under conditions of high Mach number, pressure gradient, and simultaneousmass transfer. . . The authors have set a standard of comprehensiveness which will be both an inspiration and a challenge to later
workers.”
It seems that Spalding himself accepted the challenge and set
out to create a single calculation method for turbulent shear layers(including boundary layers, jets, wakes, wall jets, and duct ﬂows)
in the presence of pressure gradient and mass transfer. This is whatcame to be known as the Uniﬁed Theory. Uniﬁcation of Flow, Heatand Mass Transport was to stay as the dominant theme of his research throughout his professional career and has had a signiﬁcant
impact on the CFD community.
In 1964, Spalding published an Aeronautical Research Council
Report [21] with the title ‘‘A Uniﬁed Theory of Friction, Heat Transfer and Mass Transfer in the Turbulent Boundary Layer and WallJet”. In December 1964 the report was retyped and annotated[22]: there were 87 typed foolscap pages and 29 ﬁgures. It was
the true start of the extraordinary contributions made by Spalding
to the prediction of turbulent shear ﬂows. As if the title were notambitious enough in scope, Fig. 1 of the paper (reproduced here)
revealed the extent of his vision although, as he makes clear on
p. 1 of the paper, not all the features represented are usually present simultaneously. Remarkably, the paper was never published ina refereed journal even in an abbreviated form. A possible explanation is that developments were moving so rapidly, both within his
own group and elsewhere, that Spalding had decided the Uniﬁed
Theory would be superseded very quickly and was not worth pursuing further.2

2
Spalding himself commented recently that he felt the paper was ‘‘too short for abook and too long for a journal paper”.

3887

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

On p. 4 we read‘‘The new theory rests on two main postulates. The ﬁrst is that thevelocity, temperature and concentration proﬁles can be describedby formulae having two main components, one accounting foreffects of momentum, heat and mass transfer to the wall, and theother accounting for interactions with the main stream; thus thegeneral proﬁles have both ‘‘boundary-layer” and ‘‘jet” components.
The second is that ﬂuid is entrained into the wall layer in the samemanner that it is into a turbulent jet and in accordance with similar
quantitative laws.”
It is interesting to note that prior to the Uniﬁed Theory (UT),
Spalding’s publications indicated little interest in turbulent ﬂows,
although the 1961 paper by Ricou and Spalding [23], which wasconcerned with entrainment by turbulent jets, is a notable exception. Spalding went on to say that once the ‘‘new” theory had been
conceived it soon became apparent that its elements were not newat all. He suggested that it might be helpful to regard the Uniﬁed
Theory as a combination of the velocity-proﬁle formulae of Coles[24] and others with the entrainment hypothesis of Head [25](although that was not the genesis of the theory).
The theory was developed in detail in Section 2 of the report.
The integral forms of the boundary-layer equations for momentum
and a conserved property were written and the velocity proﬁle represented ‘‘with sufﬁcient accuracy” by the relation1

z ¼ s2 uþ ðyþ Þ þ ð1  zE Þð1  cos pnÞ=2

ð6Þ

where z = u/uG is non-dimensional velocity, uG being the freestream velocity, with an appropriate form for the wall function
u+(y+) to include the effects of mass transfer. The key feature of thisequation, which is very close to the well-known relation proposedby Coles, is the appearance of the parameter zE which, it was said,‘‘will assume great importance later”. An analogous equation was
introduced for the conserved property u which could stand forthe composition of a chemically inert component of a mixture, stagnation enthalpy, etc. A preliminary entrainment law was then introduced guided by information about the plane free turbulent mixing
layer. Head had attempted to correlate the entrainment rate with ashape factor Hl, deﬁned as (yG  d1)/d2, yG being the velocity boundary-layer thickness whereas Spalding, displaying his usual physicalinsight, realized that the entrainment rate was more likely to be
proportional to (1  zE), the strength of the jet or wake componentof the velocity proﬁle. In fact in the original Uniﬁed Theory [21] hedeveloped a more complicated pair of equations
1
zE < 1 : mG ¼ 0:1023ð1  zE Þ 1 þ zE3
1
zE > 1 : mG ¼ 0:09ðzE  1Þ 1 þ zE =ð1 þ zE Þ3ð7Þ

wherein mG is the non-dimensional entrainment rate. Both equations were later simpliﬁed by Spalding [22] as they were thought
to be over-elaborate and found to over-estimate entrainment rates
and were modiﬁed further by him shortly thereafter [26]. The Uniﬁed Theory was shown to be consistent with many other aspectsof turbulent boundary-layer theory and, indeed, showed how a
uniﬁed approach could be developed to encompass such inﬂuences as heat and mass transfer, roughness, and even property
variations.
The responsibility for developing the purely hydrodynamic aspects of the Uniﬁed Theory for constant-property boundary layers
was given to two research students, Escudier and Nicoll. Analysis ofexperimental data available in the literature for turbulent boundary layers (zE < 1) was carried out by the former while the latter
was concerned with wall jets (zE > 1) where few data were avail-able and so extensive experimental work was required. Entrainment rates were derived using the formula

mG ¼

dRm
dRx

ð8Þ
R yG

Rx

where Rm ¼ 0 ðudy=mÞdy and Rx ¼ 0 ðuG dx=mÞdx. Although yG was
estimated systematically by ﬁtting the velocity-proﬁle equation tomeasured data, the uncertainty inherent in this procedure was
obvious. A computer program was developed to predict boundary-layer and wall-jet development based on the Uniﬁed Theory.
This code was programmed, originally on punched ﬁve-hole tapeand later IBM cards, and ran on the University of London Atlas3
computer. The outcome of the Escudier–Nicoll collaboration waspresented in [27], in which the ﬁnal version for the entrainment rateequations was

zE < 1 : mG ¼ 0:075ð1  zE Þ
zE > 1 : mG ¼ 0:03zE  1

ð9Þ

Extensive calculations were presented and discussed with the ﬁnal
conclusion that ‘‘The boundary-layer predictions are in reasonableagreement with experiment, except where separation conditions
are approached” and ‘‘The wall-jet predictions are good, both forzero pressure gradients and for adverse pressure gradients”. On
the other hand Escudier and Nicoll gave the usual explanation forthe poor performance for adverse-pressure-gradient boundary layers, namely, to blame the experimental data for not satisfying thetwo-dimensional condition. There is some truth here; see Fellouah
and Pollard [28].
It has already been mentioned that in the year following publication of the Uniﬁed Theory, Spalding modiﬁed the equation relating the non-dimensional entrainment rate, mG, to the wakeparameter zE [26]. In fact, it was in the latter paper that Spaldingﬁrst considered an approach to boundary-layer prediction based
upon solution of the momentum-integral equation and the kinetic-energy-deﬁcit equation of the turbulent boundary layer. He
showed that there was a formal relation between the entrainment
function, mG, and the dissipation integral, s, which appeared in thewell-established theories of Rotta [29] and Truckenbrodt [30]. Hewent on to recommend improved, but provisional, expressionsfor s valid for zE < 1, including near-separating boundary layers,
and for zE > 1. Spalding stated ‘‘The author’s present opinion is that
it will ultimately be preferable to use the ‘‘kinetic-energy” method
rather than the ‘‘mass-conservation” method”. His reason was ‘thegreater directness with which s can be connected with density-variation, radius-of-curvature and other effects’. It is clear that theentrainment approach was to all intents and purposes history in
the eyes of Spalding by 1965: he says ‘‘it should not be forgotten thatthe role of the entrainment-based theory, although perhaps brief, wasa valuable one”, but there was no hint that the dissipation-integral
approach was also doomed to be short lived.
5. From uniﬁed theory to computational ﬂuid dynamics
By 1965, UT was well established and it was Spalding’s prime
focus and his expectation that UT would be generally valid for avery broad range of ﬂows beyond the so called ‘‘boundary” layer
ﬂows. UT was a ‘‘grand design” built upon the insights of Taylor[31] to have a single theory that covered the canonical ﬂows ofboundary layers, wakes and jets. The essence was the use of classical proﬁle methods which had been perfected by Polhausen andothers. Spalding was aware that the proﬁle methods, then in vogue,
had severe limits. An analytical proﬁle had to be devised that represented the ﬂow well, and had the right number of parameters to3

http://www.chilton-computing.org.uk/acl/technology/atlas/overview.htm.

3888

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

satisfy the initial and boundary conditions. This was a demanding
and complex task. Spalding therefore proposed to generate universal proﬁles by representing them with piece-wise segments andthen derive the matching constants from the initial and boundary
conditions and other mathematical constraints. Although in principle any proﬁle could be used, Spalding decided to rely on simple
linear or polynomial segments. Such a piece-wise proﬁle could
approximate to a given accuracy any ideal proﬁle that might describe the ﬂow. He worked towards it with a series of students
and had a fair amount of success. His previous students had already determined optimal entrainment functions, log-law constants, and heat and mass transfer resistance required to describe
a wide range of ﬂows. A research student, S.V. Patankar, had a positive outcome on the theoretical side in building a general purpose‘‘integral-proﬁle” computer code based on piece-wise linear
segments.
Prediction of the hydrodynamic aspects of boundary-layer
ﬂows, based on solution of the momentum and kinetic energy
equations, was later extended to the temperature equation byusing two parameters for the temperature proﬁle: the thermalboundary-layer thickness and a wake-like parameter. These wereobtained from the energy equation and a ‘‘T-squared” equation,
which was the counterpart of the kinetic energy equation. Thiswas published at the 1966 International Heat Transfer Conferencein Chicago [32]. At this time, Spalding effectively started using theidea of weighting functions. The kinetic energy equation and the‘‘T-square” equation did not require a physically understandablebasis; they could be considered as integrals of the momentum
and energy equations after multiplication by a weighting function(which could be the dependent variable itself). In the search for
greater ﬂexibility and wider applicability, many families of weighting functions were considered. For example, the use of powers of u
(or T) led to the kinetic energy and related equations. Powers of y(as weighting functions) led to various ‘‘moment of momentum”
equations. Another option was to use subdomain weighting functions. With so many integral equations available, the proﬁles could
now have a large number of parameters; they could be polynomials, or better still, piece-wise linear. During this time of exploration, Spalding also suggested a variety of independent variables
for the cross-stream coordinate in the boundary layer. The distancey or the dimensionless distance y/d were the obvious candidates;
but he suggested the stream function (and its dimensionless counterpart), which simpliﬁed the representation of the cross-stream
velocity or ﬂow rate.
Spalding was at that time conﬁdent that many ﬂows of engineering interest could be represented by his uniﬁed theory andthe piece-wise proﬁles. Although the theory had produced a number of remarkable advances, ﬂows with strong favourable or adverse pressure gradients remained outside the scope of UT. Suchﬂows included important classes of industrial ﬂows such as wall
jets and separated ﬂows, including ﬂow behind a pipe step or overa cavity. The major problem with the uniﬁed theory was the needto calculate proﬁle parameters, using a set of algebraic equations.
These equations could develop singularities, e.g., in the transitionof a wall jet to a boundary layer.
Shear plays a key role in separated ﬂows, including those where
the boundary layer is destroyed by, say, an adverse pressure gradient, or a geometry that induces separation. Since the core of UT relied on the parabolic version of the Navier–Stokes equations, suchﬂows were outside the scope of UT. Spalding’s idea was that proﬁles would be found that could evolve and represent these ﬂows
which otherwise had no ‘‘self-similar” behaviour. He therefore assigned two new students to tackle the problem of bringing suchﬂows within the scope of the uniﬁed theory: Wolfshtein had already joined the group in August 1964 and Spalding asked him
to tackle the problem of the impinging jet on a ﬂat plate; ﬂow

behind a backward-facing step in a channel and that in a drivenlid square cavity were to be the focus of attention for Runchal.
These extensions would have ﬁrmly established the Uniﬁed Theory
not only for parabolic ﬂows such as the boundary layers but also
for elliptic ﬂows with strong pressure gradients, recirculation andimpingement. By the end of 1965, from the preliminary work ofRunchal and Wolfshtein, it became apparent that UT had two serious shortcomings as far as separated or high pressure gradientﬂows were concerned. The underlying structure of UT was the parabolic Navier–Stokes equations and included an inherent assumption that ‘‘self-similarity” existed. Spalding realized that the UT andthe proﬁle method had a fatal ﬂaw. There was no easy way to represent the role of the axial diffusion that destroys the self-similarity and plays a key role in separated and reversed ﬂows. The UT
could not be extended to such ﬂows. During this attempt to extend
the UT to strong pressure gradient ﬂows, it also became apparent
that the only general method then available to tackle the full (elliptic) form of the Navier–Stokes equations was the method of ﬁnite
differences. In one way the ﬁnite difference method (FDM) was
similar to the central idea of the Uniﬁed Theory: both used
piece-wise polynomial proﬁles. But there was an essential difference also. While UT relied on solving for the matching constants
that completed a proﬁle, the FDM used a simple universal piecewise proﬁle to solve directly for the value of the governing variable.
In essence there is no explicit proﬁle beyond the immediate neighbourhood of the point of interest (or node) in the computational
domain. This may be thought of as a locally evolving proﬁle rather
than a generic proﬁle.
Finite difference methods for Navier–Stokes equations havebeen around for a long time: Thom [33] had used them well before
the advent of electronic computers and Burggraf [34] had more recently applied them to a separated ﬂow in a square cavity. It became apparent that it would be simpler and more general to use
the FDM method rather than to modify the UT. In January 1966,
Runchal and Wolfshtein combined their individual research goals
and started writing a general purpose FDM computer program.
Since their immediate interest was in solving two-dimensional
problems, the stream function–vorticity (w–x) form of the Navier–Stokes equations was the focus of their attention. The w–x
approach eliminates pressure from the equations and also results
in a set of two rather than three coupled partial differential equations. Early on, it became apparent that there were limitations of
the standard, and then preferred, option of using central differences for FDM. The solution failed to converge for high Reynolds
number ﬂows. Stability analysis showed that this was because
the matrix coefﬁcients were no longer positive-deﬁnite at highReynolds numbers. Spalding then made the analogy of how thewind from the pigsty always stinks – or that northerly winds (inthe Northern hemisphere) generally bring the cold. These discussions led to the development and use of the ‘‘upwind” concept. It
later turned out that one-sided differences had been previously
used by Courant et al. [35] and others. The important contribution
made by Spalding, was the adoption of the concept of a ﬁnite-volume method (FVM) as the core of mathematical development
rather than the use of the traditional FDM.4 Spalding likened thecomputational grid to a series of ‘tanks’ (that hold the nodes) andthe ‘tubes’ (that convey the ﬂuid along the grid). These improvements rapidly led to success. By June of 1966, Runchal and Wolfshtein reported results both for the wall jet and for the separated ﬂows
behind a backward-facing step and in a driven-lid square cavity for
low and high Reynolds number ﬂows. This was reported in a seriesof papers from Imperial College [36–39]. By late 1966, the method4
It was not until some time later that the terminology ‘‘ﬁnite volume” wasadopted.

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

was well established and began to be used widely by other researchers all over the world. It was subsequently published in book form by
Gosman et al. [40].
One shortcoming of the w–x method was that it was essentiallylimited to two-dimensional ﬂows and could not easily be extended
to three dimensions in any obvious manner. The possible extension
using the vector vorticity and the vector velocity potential (whichis the counterpart of the stream function) was thought at the time
to be far too complicated; it was feared that the delightful simplicity in the two-dimensional arena would simply go away and as aresult further advances (as documented below) were based on
the so-called primitive variables, namely the velocity components
and pressure. However, the vorticity–vector potential formulation
was in fact later successfully used by Mallinson and de Vahl Davis[41] and many others.
With the development and the success of the FVM and theinability to extend the UT to ﬂows with strong pressure gradients,
the progression of UT (based on integral equations) came to an
abrupt halt in the summer of 1966 when Spalding travelled with
Patankar to the 1966 International Heat Transfer Conference in Chicago to present a paper on their work on the UT [32] and to work on
a consultancy project for the summer break. During this period,
with the already demonstrated success of the FDM for elliptic ﬂows
at the back of their minds, they decided to essentially abandon the‘‘piece-wise-proﬁle” method for UT and switched over to FDM forthe parabolic form of Navier–Stokes equations for two-dimensionalboundary layers. They introduced the innovation of replacing the
distance normal to the wall by a non-dimensional stream function
based variable – essentially the von-Mises transformation – as the
cross-stream coordinate, a piecewise-linear proﬁle for velocity (andtemperature), and the subdomain weighting function. Suitableentrainment formulas were used at free boundaries to determine
the growth of the calculation domain into the surrounding freestream. This method had the speed, ﬂexibility, accuracy, and simplicity needed for a general procedure for two-dimensional boundary layers, free-shear ﬂows, wall jets, and duct ﬂows.
Thus by the end of 1966, both the elliptic and parabolic FDMs
were ﬁrmly established and the entire focus of the group shiftedto application and improvement of this technology. With the publication of the papers [42,43] and later the books by Patankar and
Spalding [44,45] and Gosman et al. [40] both the elliptical and theparabolic methods gained considerable following and were widelyused by researchers at Imperial College, Stanford University, andelsewhere. A number of Ph.D. theses applied the methods to a variety of parabolic and elliptic ﬂows. This was also the time when research on mathematical models of turbulence was in full swing.
The CFD codes provided a convenient vehicle for testing turbulence
and combustion models. The two-ﬂux and six-ﬂux models for radiation [46] and formulations for particle concentrations in different
size ranges were also implemented in the framework of thesemethods.
By 1968, evidence began building up that for many applications
the FDM and FMV with upwind differences led to undue false(numerical) diffusion. The effort therefore shifted to ﬁnding superior methods for representing the convective terms. A number ofhigher order schemes, including the concept of vector-differencing
following the stream-lines, were explored, although none proved
entirely satisfactory. Spalding [47] then proposed an exponential
method to replace upwinding, but eventually the Imperial College
group settled on a ‘‘hybrid” difference scheme [48] that blended
the central and upwind difference methods based on the local
Péclet number. For quite some time this [40,43] became the standard approach for solving high Reynolds number ﬂows, but was later augmented by higher order upwinding, QUICK [49] andvectored upwinding schemes, and by the use of ﬂuxes, not proﬁlesof variables.

3889

It was realized that the ‘‘tank-and-tube” approach could be replaced by the Gauss Theorem. This avoids the use of the Taylor series and leads to an elegant formalism to derive the discrete set of
equations in any arbitrary shaped control volume. This generalized
FVM revolutionized the engineering practice by giving rise to the
new ﬁeld of computational ﬂuid dynamics (CFD). General and ﬂexible computer codes could be written to solve ﬂuid ﬂow problems
in any complex geometry at any Reynolds number. In an important
way the FVM also changed the basic thinking behind computational ﬂuid dynamics. Unlike the mathematical approach used by
FDM and the ﬁnite element method (FEM), an understanding of
the FVM was enhanced by a physical approach to the problem.
The focus of interest was not the variable and its interrelationshipparameters with other nodes; the focus now was on the control
volumes and the ﬂuxes crossing their boundaries and the sourcesor sinks contained within them.
Further research, and the continued creativity of Spalding, led toa new and improved version of the parabolic boundary layer method, which was published in the 1970 edition of the book by Patankar and Spalding [45]. This contained a new computer program
called GENMIX. In 1977 Spalding published a book [50] which isin the nature of a user manual for GENMIX. The method now included an improved way of calculating the entrainment and a procedure for determining the unknown pressure gradient in conﬁned
ﬂows. The new method now contained the hybrid scheme developed earlier for the elliptic model for high Reynolds number ﬂows.
The parabolic boundary layer procedure of Patankar–Spalding
was based on a non-iterative scheme that employed a purely‘‘marching” procedure. Spalding explained the non-iterative determination of the pressure gradient in a conﬁned ﬂow via a ‘carsteering analogy’. The calculation of the entrainment at the freeboundaries also remained a source of oscillations and instability
in many applications of the method. Formally, the equations are
nonlinear and an iterative solution would have addressed this in
a forthright manner. (Later in the general-purpose CFD code PHOENICS, described below, Spalding [51] addressed these issues anddid employ iteration at the forward step.) The entrainment rates(and the unknown pressure gradient in a conﬁned ﬂow) shouldhave been recognized as eigenvalues in the equations to be determined from certain overall properties of the resulting solution.
With what we know now, a much more robust and efﬁcient procedure for two-dimensional boundary layers could now be constructed; but there is hardly an incentive to develop such amethod in this day and age.
In 1967 Spalding made a contribution that on the face of it
seems trivial but had a profound inﬂuence on later thinking anddevelopment of CFD tools. He introduced the concept of a ‘‘General
Transport Equation” to express all second order convective-diffusive transport equations. He proposed to write them as a single
equation for a generalized variable u that represented velocity
components, enthalpy, or chemical species. This generalized transport equation has four components: (1) transient or accumulation,2) convection, (3) diffusion, and (4) source term. Thus the Navier–Stokes, energy and mass (species) transport equations could all be
expressed by this single equation:@
ðq/Þ þ divðqu/Þ ¼ divðC grad /Þ þ S
@t

ð10Þ

Early examples of this formulation appear in the 1969 book by Gosman et al. [40]. This was his effort at ‘‘uniﬁcation” of mathematicalexpression. This approach led to the fast development of CFD methodology since one can concentrate efforts at solving this single
equation rather than look at the speciﬁcs of each variable.
Another innovation that Spalding introduced around 1966 wasthat of the ‘‘compass” notation to express the FD or FV stencil. This
simpliﬁed the visualization and expression of the ﬂuxes and3890

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

meshed in with his ‘‘tank-and-tube” approach to FV. The six points(or control volumes) surrounding a central point (or CV) were designated ‘‘north”, ‘‘south”, ‘‘east”, ‘‘west”, ‘‘high”, and ‘‘low”. Thus,the ‘‘standard” form of the linear algebraic equations, written in
the so-called zero-residual form, isaW ð/W  /P Þ þ aE ð/E  /P Þ þ aS ð/S  /P Þ þ aN ð/N  /P Þþ aL ð/L  /P Þ þ aH ð/H  /P Þ þ S ¼ 0

ð11Þ

By 1969, the interest of the Imperial College group began toshift to three-dimensional ﬂows. The stream function–vorticity
variables had distinct advantages for two-dimensional ﬂows but,as previously mentioned, their use for three-dimensional ﬂows
was problematic. After some preliminary (w–x) extension to suchﬂows, Spalding started to look for methods based directly on theprimitive form of the Navier–Stokes equations. A number of alternatives were tried but none proved entirely satisfactory. The ﬁrst
hint of success came with the coupled SIVA scheme of Caretto
et al. [52]. However the coupled scheme was deemed to be computationally expensive and the focus shifted to the search for an efﬁcient sequential or segregated method of solving the governing
equations. By this time Harlow and co-workers [53] had demonstrated the advantages of a staggered grid algorithm and Chorin[54] had published work on pressure splitting (or pressure projection) algorithms. Patankar and Spalding [55] combined these ideas
with the FVM technology and hybrid-differencing developed earlier to arrive at a robust and general 3D Navier–Stokes solver called
SIMPLE (Semi-Implicit Method for Pressure-Linked Equations). The
SIMPLE algorithm went on to achieve a classic status, and led to thecommercialization of the new CFD technology.

tions for the two-dimensional boundary layer contain the same
inconsistency about pressure? Spalding was able to point out thatthe issue was present all along in the two-dimensional boundary
layer, but remained invisible; it is the three-dimensional boundary
layer that brought the issue to the forefront.
In a two-dimensional boundary layer, we obtain the streamwise
velocity by solving the streamwise momentum equation, in which
we use a single value of dp/dx. Then we obtain the cross-stream
velocity from the continuity equation. The cross-stream momentum equation does not get used at all! If we were to use it, we
would be able to calculate the variation of pressure in the crossstream direction, thus exposing the inconsistency with the valueof dp/dx used.
In a three-dimensional boundary layer, there are two crossstream velocities. The continuity equation alone cannot determine
these two velocity components. How the ﬂow distributes in thetwo cross-stream directions is decided by the variation of pressure
over the cross-stream plane. We are forced to solve the crossstream momentum equations and calculate the variation of pressure over the cross-stream plane. The pressure inconsistency cannot be hidden under the carpet!
The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) was described in the paper by Patankar and Spalding [55].
The main ingredients of the method are:
 calculation of the streamwise velocity component by solving thestreamwise momentum equation with a value of dp/dx thatremains the same across the cross-stream plane calculation of the cross-stream velocity components and crossstream pressure variation by solving the cross-stream momentum equations and continuity.

6. The SIMPLE algorithm
Around 1970, there was an effort to develop a procedure forthree-dimensional boundary layers (or three-dimensional parabolic ﬂows). It was quickly realized that this was not a simple or
straightforward extension of the two-dimensional boundary layer
procedure. The three-dimensional boundary layer contained within it a parabolic behaviour in the streamwise direction and an elliptic behaviour in the cross-stream plane. For this ‘‘elliptic” problem,the stream function–vorticity procedure could not be used because
the three-dimensionality of the ﬂow precluded the use of thestream function as a variable. The obvious choice was to use thevelocity components directly.
The breakthrough contribution by Spalding was the identiﬁcation of the role of the pressure gradient. To preserve the marching
nature of the solution procedure for this parabolic problem, Spalding required that the streamwise momentum equation was driven by a single pressure gradient dp/dx across the wholeboundary layer at a given x. This pressure gradient could be obtained from the pressure gradient in the external free stream (forunconﬁned ﬂows) or from the overall mass conservation (for conﬁned ﬂows). Then the cross-stream velocity components were governed by their momentum equations and a pressure variation in
the cross-stream plane that was needed to satisfy the local continuity equation. Although this pressure variation in the crossstream plane was calculated, it was not used in modifying thestreamwise pressure gradient dp/dx used for the streamwise
momentum equation. Thus, an inconsistency was deliberately accepted in the treatment of pressure. Spalding maintained that thisinconsistency was a small price to pay for the beneﬁt of using themarching procedure for the parabolic ﬂow. In practice, the crossstream variation in pressure would be rather small and no signiﬁcant error was involved.
If this was an important issue, why did we not know about it
while studying the two-dimensional boundary layer? Do the solu-

The latter step required the development of a new algorithm,
which is the heart of SIMPLE. It was recognized that we could always solve the momentum equations for any speciﬁed pressureﬁeld; but then the resulting velocity ﬁeld would not normally satisfy the continuity equation. We could then propose ‘‘corrections”
to the speciﬁed pressure ﬁeld so that the resulting (corrected)
velocities would satisfy continuity. This led to the birth of what
was called the ‘‘pressure correction equation”. It was actually arewritten form of the continuity equation, which gave a directmeans of calculating the pressure corrections (and hence thepressure).
Although in the 1972 paper, Patankar and Spalding presented
the SIMPLE algorithm in the context of three-dimensional parabolic ﬂows, it was obvious that the method contained within it aprocedure for two-dimensional elliptic ﬂows, since that is how
the cross-stream velocities were treated. Further, the method
could easily be extended to three-dimensional elliptic ﬂows, as it
did not rest upon any concepts (such as the stream function) thatwere speciﬁc to a two-dimensional situation. Thus, SIMPLE provided a general calculation procedure for solving the coupledmomentum and continuity equations in any dimensionality. Onceagain, this method became very widely used all over the world. Improved variants of the method, such as SIMPLEC [56], SIMPLER[57], and SIMPLEST [58] were later developed by different
researchers, in part to strengthen the coupling between velocity
and pressure and provide more robust algorithms with fasterconvergence.

7. Multi-phase ﬂow
During the mid and late 1970s, Spalding began extending SIMPLE-like methods to multi-phase ﬂows. The InterPhase Slip Algorithm, IPSA, is the two-phase Eulerian–Eulerian algorithm which

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

Spalding pioneered [59,60]. It is a natural and elegant extension of
the single-phase SIMPLE algorithm. Although in its original formulation and the description below, only two phases were identiﬁed,
all the main ideas are easily extendable to multi-phase (i.e., morethan two) ﬂows. It ﬁrst appeared in 1977 at an ICHMT Seminarin Dubrovnik, see Spalding [60].
Spalding considered a two-phase ﬂow as composed of two
interpenetrating and intermingling continua. Each location (andas a result, each ﬁnite-volume cell) potentially had both phases(e.g., liquid and gas) present. The amount of each phase was characterized by its volume fraction. Such a model of two-phase ﬂow
does not describe how the two phases are locally distributed. The
volume-fraction characterization cannot distinguish between large
chunks of liquid and ﬁne droplets. It is this characterization that
leads to the simplicity of the formulation and also to its main
weakness. Spalding later presented extensions for removing thisweakness [61].
For the two interpenetrating continua, he proposed two sets ofconservation equations for the phases. These were similar to thesingle-phase equations, except for an additional pair of inter-phasetransport terms. Thus, the velocity difference between the twophases at a point would lead to a force term due to the shear force
exerted by one on the other, and temperature difference between
the two phases would lead to an inter-phase heat transfer. The inter-phase transport strongly depends on the surface area available
for shear or heat transfer. This in turn requires the size of the droplets of the particulate phase. The volume-fraction characterization
does not provide this information. Indeed the multi-phase ﬂow
problem is a classic example of a closure problem, of which turbulence is, perhaps, the example best known to the engineering
student.
It seemed practical to solve the conservation equations for thetwo phases in a segregated manner; that is, to solve for the velocities and temperatures of one phase followed by those for the otherphase. Thus, while solving for one phase, the values of velocity andtemperature for the other phase were assumed to be known fromthe previous iteration. When inter-phase transport was dominant,this procedure would converge very slowly since the assumed values of the other phase would dominate the behaviour of the current phase. Spalding had already seen this problem in anothercontext. In a heat exchanger analysis, discussed elsewhere in thisarticle, the temperatures of the shell-side ﬂuid and the tube-side
ﬂuid inﬂuence each other in the same manner. Spalding thus invented the Partial Elimination Algorithm (PEA) precisely to handle
such a situation. The trick was to solve the equations for temperature in an algebraically transformed form that eliminated the direct contribution of the ‘other’ temperature. (There are still someterms that indirectly contain the neighbour values of the ‘other’temperature, which is why the transformation is partial and not total.) The reduction of the inﬂuence of the ‘other’ temperature (orvelocity) gave an impressive speedup of the iteration process. Spalding immediately identiﬁed the two-phase ﬂow problem as a candidate for PEA. It turns out that the use of the PEA is essential forthe success of most two-phase ﬂow calculations. A ‘sum-to-unity’requirement was also imposed and ensured that the coupling between the phase volume fractions is preserved at all times during
the solution procedure.
If the overall continuity equation is written as a mass-conservation equation, the terms for the denser ﬂuid (liquid) excessively
dominate the terms for the lighter ﬂuid (gas). To overcome this difﬁculty, Spalding wrote a ‘‘volume continuity equation” [62] whichgave comparable weighting to the two phases. The Gas And Liquid
Analyser (GALA) employs volumetric rather than mass continuity,and is advantageous for computing the movement of ﬂuids containing discontinuities of density, especially those of large magnitude such as between a gas and a liquid. This overall continuity

3891

equation forms the basis of creating the pressure or pressure correction equation for two-phase ﬂows. Although the notion of a single ‘shared’ pressure is fundamental to the procedure, IPSA permits
there to be two different ﬂuid pressures, provided an appropriate
algebraic relationship between the two individual pressures isknown. Thus, the phase continuity equations are used to obtainthe volume fractions of the individual phases and, by analogy to
SIMPLE, overall continuity is converted to a pressure correction
equation.
Like SIMPLE, IPSA was rapidly adopted by researchers aroundthe world, and variants of IPSA form the basis for two-phase Eulerian–Eulerian formulations in almost all multi-phase ﬁnite-volume
CFD codes today. The basic framework for two-phase ﬂows was later extended, conceptually, to multi-phase ﬂows of various kinds,
using the ‘algebraic-slip method’ (sometimes called the mixture
model or drift-ﬂux model). In essence, the algebraic slip model iscomprised of continuity and momentum equations for the mixture,continuity equations for the dispersed phases, and algebraic relationships for the slip velocities. One could consider liquid dropletsin different size ranges, gas packets of different concentrations, andso on. In the subsequent implementation of IPSA in the PHOENICS
code, a phase-diffusion term was also included in the phase continuity equations. Spalding modelled this term by means of a gradient diffusion assumption so as to represent the turbulent ﬂux
associated with correlations between ﬂuctuating velocity and volume fraction. Some later researchers argued for excluding thephase-diffusion term in favour of a turbulent dispersion force inthe momentum equations, e.g., Lahey et al. [63]. The phase-diffusion terms proved to be an important element in some variantsof Spalding’s two-ﬂuid model of turbulence, Malin and Spalding[64]. Later, (i.e., post 1980), Spalding used the two-phase concepts
to devise advanced turbulence models [64–67] (considering the
turbulent ﬂow as a two-phase mixture of laminar and turbulent
chunks) and advanced combustion models [68–70] (considering a
mixture of reacted and unreacted gases).

8. Turbulence modellingWith the development of the Patankar and Spalding [42,44,45]
ﬁnite-volume solver for two-dimensional thin shear ﬂows, theneed (so far as turbulence modelling was concerned) shifted toproviding a route for determining the effective turbulent stresses
and heat ﬂuxes throughout the computational domain, a problem
that Spalding simpliﬁed to that of obtaining the local turbulent viscosity. The initial model in the Patankar–Spalding method was Prandtl’s mixing-length hypothesis, a scheme employed in several
early papers including his group’s submission for the Stanford Conference, Kline et al. [71].
Spalding’s vision, however, was already on providing a corresponding computational procedure for recirculating ﬂows and here
he recognized that the mixing-length scheme was quite inadequate. In Ref. [72], he sets out to explain why, in a separated ﬂow,
the wall heat ﬂux often reaches a maximum at the reattachment
point. As he remarked [72]‘‘This fact assumes particular interest when it is recalled that, at thereattachment point, the (time-mean) shear stress is zero; and mostformulae for calculating heat transfer rates which are based on aphysical theory of ﬂow in a boundary layer would predict that,
where the shear stress is zero, the heat transfer must also vanish.”
The apparently paradoxical behaviour of separated ﬂows he explained as due to the diffusional transport of turbulence energy, k,
towards the wall. He provided evidence for this assertion bycarrying out an approximate analytical solution of the turbulenceenergy equation, showing that, with the empirical coefﬁcients

3892

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

appearing in the diffusion and dissipation terms chosen from quite
different ﬂow types, the observed power-law dependency between
Stanton and Reynolds numbers was correctly reproduced. The turbulence energy equation had been employed over the preceding
twenty years by several turbulence modellers, e.g., Prandtl [73],but Spalding’s paper was the ﬁrst to show an application of theequation to an idealization of a separated ﬂow and how it resolved
the question of why heat transfer in separated ﬂows followed different rules of behaviour from those in boundary layers. It was presumably from this research that the decision was taken to embed
the Prandtl–Glushko turbulence energy model into the slightly later numerical scheme for recirculating ﬂows, Gosman et al. [40]).
A second major contribution, linked to the foregoing, concerned
improving the near-wall logarithmic velocity proﬁle (onto which
the numerical solution for the outer region was ﬁxed to avoid
what, at the time, was the impossible numerical requirement –not to mention the modelling problems – of resolving the viscous
sublayer in an elliptic solver). As noted above, Spalding recognizedthat the usual formula

u

þ

u 

1

ðsw =qÞ2

¼

1

j

"
ln

1

Eyðsw =qÞ2

#
ð12Þ

m

and the corresponding temperature distribution
12

þT  qðsw =qÞ cp DT=qw ¼ ðrt =jÞ ln"#1Eyðsw =qÞ2mþ Pðrt =rÞð13Þ

were deﬁcient in several respects in the vicinity of stagnation points
1
simply because there, the friction velocity ðsw =qÞ2 was an inappropriate velocity scale. Having shown the necessity of including the
turbulence energy in computations of separated ﬂows, he proposedthat the dimensionless parameters appearing in the above law-ofthe-wall formulae should be changed to:u " 1 1#
1 1
uc4l k2
Ec4l k2
1¼ ln
ðsw =qÞ j
m
1
4

1
2

T  qcl k cp DT=qw ¼ ðrt =jÞ ln

ð14Þ
" 1 1#
Ec4l k2mþ Pðrt =rÞð15Þ

The starred quantities introduce the turbulence energy, k. (Inequilibrium, u*, T* are the same as u+, T+, but not otherwise.)
The change from Eqs. (12)1 and (13) to Eqs. (14) and (15) simply
involved replacing sw/q by c2l k. Indeed, the coefﬁcient cl was chosen as 0.09 so that, in an equilibrium ﬂow near a wall (where thewall shear stress is proportional to the turbulence energy in the‘‘log-law” region), the formula gives the same results as the more
traditional form, Eqs. (12) and (13). A further important featureof Eq. (14) is that the wall shear stress becomes directly proportional to the wall-parallel velocity which is important when theﬂow is close to separation. Spalding’s version of the law of the wall
was employed in software at Imperial College from the beginningof the 1970s and was ﬁnally published in the review paper by
Launder and Spalding [74] to which further reference is made below. The form became, and remains, the standard wall treatment
for most of the commercial CFD codes that have emerged overthe past 30 years. Unlike the original law-of-the-wall, it always
gives at least qualitatively correct heat transfer patterns in recirculating ﬂows, although more reliable wall-function schemes arenowadays employed in research-level CFD.
Despite the very real beneﬁts of introducing the turbulence energy to compute recirculating or other far-from-equilibrium ﬂows,
Spalding realized that, because the Prandtl–Glushko model (likethe mixing-length hypothesis) needed an empirical speciﬁcationof the turbulent length scale, it was extremely limited in its widthof applicability. He saw that more general and more complete

models required, in addition to an equation for k, a further transport equation for calculating the distribution of the turbulent
length scale over the ﬂow ﬁeld. Inspired by the work of Kolmogorov [75] and Rotta [76] who proposed length-scale equations forthe ﬂuctuation frequency, x, and the product of k and the lengthscale, L (i.e., a kL-equation) respectively, Spalding pushed forward
the development of models using such an equation, i.e., the development of two-equation models. In an early report [77] he used
Rotta’s length-scale equation in integral form to show that various
free turbulent shear ﬂows (notably the mixing layer, the plane jet,and the fan jet) could be calculated with the same set of constants,
while, in contrast, the mixing-length model required different ratios of mixing length to ﬂow width to give the correct rate of
spread. Then, to support these ﬁndings by numerical calculations
and to develop a general model for the computer codes his CFD
team had created, he gave two Ph.D. students the task of developing and testing, with the aid of the Patankar–Spalding boundarylayer code, a model using Rotta’s length-scale equation together
with a transport equation for k and the Prandtl–Kolmogorov
eddy-viscosity assumption. Rodi performed the development andtesting of the model for free ﬂows while Ng took on the development and testing for near-wall ﬂows. This research led to theestablishment of the k–kL model. In Rodi and Spalding [78], thesuperiority of this two-equation model over the mixing-length
model was conﬁrmed for free ﬂows. This work also brought out
the so-called plane/round jet anomaly, i.e., the fact that the plane
and the round jets in stagnant surroundings cannot be calculatedwith a model with the same constants. The near-wall calculations,
Ng and Spalding [79], showed that the model also performed fairly
well with these ﬂows, provided that an extra term was added to
the kL length-scale equation to produce the correct decay of lengthscale near the wall.
In parallel with this work, Spalding developed his k–W model inwhich he employed an equation for W  k/L2 (which is the squareof Kolmogorov’s frequency x) as the length-scale-determining
equation where the quantity W was interpreted as the time-meansquare of the vorticity ﬂuctuations. This work was ﬁrst reported inRef. [80] and then in a paper in which he also calculated the concentration ﬂuctuations in a round jet, Spalding [81], a paper which
was to serve as a basis for many other papers, especially in the
combustion community. An extended account of work on the k–
W model was provided in an internal report, Spalding [82], which
also appeared in a commemorative volume, Spalding [83]. While it
may seem strange that Spalding should adopt a different dependent variable from his students, he recognized that this was anarea where a range of choices was available and that the advantages or weaknesses of different routes needed to be explored. Indeed, he was supportive of the fact that a group led by one of
the lecturers, B.E. Launder, explored the use of an equation forthe dissipation rate e  k3/2/L as the length-scale equation, Hanjalic
and Launder [84], Jones and Launder [85].
Spalding made some comparison of various length-scale equations and addressed the question of choice in his 1971 report [82]on the k–W model. He stated in this report that ‘‘it cannot yet be said
that any model is obviously superior” – and that the choice could bebased upon aesthetic considerations. However, already in 1970 in
an internal Imperial College note addressed to his coworkers
Spalding analyzed the general length-scale equation for the variablem

Z ¼ k Lnð16Þand considered which choice of length-scale variable (i.e., which
exponents m and n) would best satisfy the requirements for bothfree and near-wall ﬂows. He came to the conclusion that the variable e (m = 3/2, n = 1), which was then being used by Hanjalic
and Jones, satisﬁed the requirement better than the other lengthscale variables investigated at Imperial College at the time. A similar

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

analysis was published later in Launder and Spalding [74] which
established the basis for giving the preference to the e-equation.
The 1970 note states, however, that Kolmogorov’s [13] variable x,
which is the frequency of the turbulent ﬂuctuations, could also bea suitable variable, an early ﬁnding that is supported by the more
recent success of the k–x model, Wilcox [86], and its ‘Shear–Stress
Transport’ (SST) variant, Menter [87].
In free turbulent shear ﬂows, the models using different lengthscale equations produced nearly identical results, as observed invarious papers by Spalding; but, of the models investigated at
the time at Imperial College, only the e-equation could dispensewith wall-effect corrections (additional terms or variation of theconstants). This fact was responsible for Spalding being converted
to giving his preference to the e-equation in the coming years.
Among his earlier involvements with the k–e model were the calculations for the 1972 NASA-Langley conference on free-shear
ﬂows [88], to which the Imperial College group contributed calculations of 23 different free-shear ﬂows with a number of turbulence models ranging from the mixing-length hypothesis to asecond-moment-closure, those involving a scale-determining
equation all using the e-equation, Launder et al. [88]. These calculations ﬁrmly brought out to a largely American audience the superiority and greater generality of the two-equation model overmodels using an empirical length-scale prescription for free-shear
ﬂows. In the years following, Spalding was also involved in calculations of elliptic ﬂows with the k–e model and this early work
was summarized in the paper of Launder and Spalding [74], which
for a long time has served as one of the main references for the k–emodel, with numerous citations in the literature.
In the 1980s Spalding returned to his k–W model and made thesame version applicable to both near-wall and free ﬂows by adding
a source term involving the length-scale gradient in the W-equation. The model was tested successfully for a variety of ﬂows, Malinand Spalding [64,89], Ilegbusi and Spalding [90]; but by then the k–e model was so widely accepted that it remained the main modelin use. Around the same time Spalding also introduced two-ﬂuid(turbulent/non-turbulent) models involving the calculation of theintermittency. One version was a modiﬁed form of k–e model, Malin and Spalding [64], while the other, Ilegbusi and Spalding [65–67], did not use an eddy viscosity but rather a diffusion ﬂux basedon the direct calculation of the cross-stream velocity of the two ﬂuids, while reverting, for the length-scale determination, to an
empirical prescription. While the ideas were novel, these models
did not have any great impact, compared to the advance of the original k–e model.
Indeed, the k–e model rapidly became the most widely used
two-equation model and, while now superseded at the research level in favour of more general routes to determining the stress–
strain inter-linkage, it is still the main workhorse in industrial
CFD calculations. While Spalding was neither the developer nor
the main advocate of the model, and in fact later expressed doubts
about the wisdom of the choice of the e-equation [91], he played amajor and decisive role in developing, promoting and disseminating models of the two-equation eddy-viscosity type which he recognized very early on as the simplest useful level for use within
general CFD computer codes since such models required no turbulence quantities to be prescribed within the computation domain.
This enabled non-specialists in turbulence modelling to embark on
the calculation of turbulent ﬂows for practical conﬁgurations overa very wide range of ﬁelds.

9. Turbulent combustion and chemical reaction
Brian Spalding has made many very signiﬁcant contributions tothe ﬁeld of turbulent combustion modelling. Indeed, many of the3893

currently prevailing key ideas in combustion modelling were already present in his pioneering developments, including the often-controlling role of mixing at the molecular level, and themulti-scale nature of the problem with the crucial part played bysurface phenomena at the ﬂame front. In Spalding’s own words,[92]:‘‘Few painters invent new colours: it is their arrangement on thecanvas which constitutes novelty.”
Spalding’s most cited contribution to turbulent combustion isperhaps his 1971 paper in the 13th Combustion Symposium [93].
From the experimental observations of turbulent, premixed, conﬁned, bluff-body-stabilized ﬂames available at the time, he reﬂected on the relative independence of the angle of spread of theﬂame from the experimental conditions (such as fuel type, equivalence ratio, mean speed of the approaching ﬂow or temperatureand pressure levels). From this evidence, he concluded that, for awide range of ﬂame conditions, the chemical kinetics were indeed
subsidiary to aerodynamic processes. In a much earlier, pre-CFD
paper published in 1967 [94], Spalding had analyzed the implications of assuming that the fresh mixture was immediately burned
as it was entrained into the shear layer (then known as the ‘‘entrained-is-burned” model). The assumption provided at the time
results in the correct direction; but, from the vantage point afforded by the ﬁrst CFD calculations years later, it was deemed to be
crude. The entrainment-is-burned model was thus replaced in
the 1971 paper with the ‘‘mixed-is-burned” model, which would
be later known as the ‘‘eddy-break-up”, or EBU, model. The essential feature of EBU is that turbulent burning proceeds at the rate at
which large parcels of fresh gas break into smaller ones, thus creating interface area for the molecular processes to act; since it ismainly on this interface area that combustion takes place, at sufﬁciently-large Damköhler numbers, it was postulated that the rateof combustion can be modelled by recourse to the rate of eddy
break-up.
The seminal 1971 paper [92] addressed ﬂames in turbulent
shear ﬂows, where, using an equilibrium assumption, the rate of
kinetic energy decay can be estimated to be equal to the production. Thus, the reaction rate per unit volume in this mixing-controlled regime was postulated to be:

 
@u
 
_r000
mix ¼ Cð1  sÞq
@y

ð17Þ

where s is the mixture local reactedness (or reaction progress). An
Arrhenius rate was combined with this mixing rate in a harmonic
average to cater for situations where chemical kinetics is the ratecontrolling process.
The model was further reﬁned in a paper published in 1975 [95]that is one of the ﬁrst endeavours at the comprehensive, multidimensional modelling of a combustion chamber. The calculation
was fully elliptic (as opposed to previous mainly parabolic simulations of laboratory ﬂames); the two-equation k–e model was used
for turbulence; and a four-ﬂux model of radiation was included.
Advantage was taken of the k and e predictions to re-cast theeddy-break-up model into one of its now widespread forms, viz.
1=2
ðqe=kÞ
r_ 000mix ¼ Cg

ð18Þ

where g is the scalar (mixture-fraction) variance [81].
The eddy-break-up model, as described above and used for several decades (with minor, non-fundamental variants) by many
researchers was largely regarded by Spalding as a ﬁrst, crude attempt [93,95] to account for the often-controlling inﬂuence of
ﬂame aerodynamics on the combustion rate, and of the relevanceof the micro-scale phenomena of interface-surface creation andof molecular mixing across this interface. The original EBU model

3894

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

was therefore much reﬁned in subsequent research by Spalding
and co-workers. Thus, in [96,97] he outlined the ideas for a‘‘slide-and-stretch” model of turbulent combustion, in which‘‘coherent bodies of gas are squeezed and stretched as they travel
through the ﬂame”. The idea was further elaborated in [98], andwas subsequently given its own character as the ESCIMO model.
ESCIMO stands for Engulfment, Stretching, Coherence, Inter-diffusion and Moving Observer, i.e., the main constituents of the model. Spalding, drawing from the ideas of many previous authors,
formulated a model [100] that can be regarded as the ﬁrst multiscale interpretation of turbulent combustion suitable for embodiment in a CFD code. Chemical reaction was, as for EBU, assumed
to take place in the interface area separating the reactants (in diffusion ﬂames) or reactants and products (in premixed ﬂames). Theinterface area was arranged in what were equivalently termed‘‘folds”, ‘‘parcels”, ‘‘sandwiches”, or ‘‘layers” (Fig. 2); these areformed in the main shear regions of the ﬂow; they stretch as they
are convected; and then ‘die’ when they are re-engulfed by othereddies. The task of ESCIMO was therefore twofold: ﬁrst, to calculate the statistics of the folds (i.e., the local distribution of theirsizes, or widths), which was called the ‘‘demographical part”[100]; and, second, to compute the (one-dimensional) property
distribution within each fold (which was termed the ‘‘biographicalpart”). The ﬁrst is an Eulerian point of view; the second is Lagrangian, because the transient, one-dimensional proﬁles depend (barring simpliﬁcations) on the fold history as it moves through theﬂow (this is the Moving Observer part of the model name). Also
the ﬁrst describes the ﬂow in the macro-scales, where convection
and stretching occur, while the second focuses on the microscales,
where molecular diffusion and chemical reaction take place. The
ESCIMO model was applied by Spalding and co-workers to several
canonical conﬁgurations, such as stirred reactors, bafﬂe ﬂames andfree jets [100,101], for which analytical solutions for the fold populations could be calculated with suitable simplifying assumptions.
It was later successfully applied to turbulent, non-premixed hydrogen–air ﬂames [100].
The uptake of ESCIMO by the scientiﬁc community at large wasscanty, despite the fact that it embodied, in a framework which
rendered the problem solvable with the late 1970s computers,

the key concepts later found in widely-accepted models. Thus,
ﬂamelet models of non-premixed combustion make use of thesame notion of a distribution of thin reaction layers, which areparameterized using not a length scale (the fold width) but the inverse of a time-scale (the dissipation rate). The time-dependence ofthe biographical part of the fold (e.g., the species one-dimensional
proﬁle within the layer, Fig. 2), is dispensed with in the so-called
steady-state ﬂamelets. This is calculated locally, irrespective of
the parcel (or ﬂamelet) history. The more recent, and independently derived, linear-eddy model (LEM) of Kerstein [102] remarkably shows some common modelling principles with ESCIMO.
Notably, LEM resorts also to a one-dimensional calculation,embedded in the three-dimensional one, to resolve the coupling
between mixing and chemical reaction in the small scales that cannot be resolved by the 3D mesh. Unlike Spalding’s ESCIMO, thesmall scales are calculated in an Eulerian framework; and the effectof turbulent mixing is simulated through stochastic stirring of the
one-dimensional ﬁeld, rather than by using a local ‘‘eddy
distribution”.
The development of Eulerian multiphase methods such as IPSA
[60,103], discussed above, opened up new opportunities to add
physical realism to turbulent combustion modelling. Spalding soonrecognized that the ﬂame front, where burned (burning) and unburned matter intermingle as in a ‘‘land-and-lagoons landscape”[68], and could be modelled as a two-phase ﬂow, with one phaserepresenting the burning ﬂuid and the other the unburned one.
Combustion is therefore regarded in this framework primarily as
an interface phenomenon taking place on the parcel boundary,and triggered by the entrainment of fresh ﬂuid into the burningone. The combustion rate was primarily governed by this entrainment rate, which required modelling. The models for entrainment
often involve the fragment surface area per unit volume, which is
the inverse of a length scale and can be interpreted as a characteristic fragment size; several algebraic and transport models havebeen suggested for this fragment size. Kinetic control of combustion is also provided within the burning ﬂuid. It is interesting tonote that a similar concept underlies the single-ﬂuid eddy-dissipation concept of Magnussen [104], where the ‘‘ﬂuids” (using Spalding’s term) are the ‘‘surroundings” or mean mixture and theFig. 2. The ESCIMO model: (a) the ‘‘folds”, from [99]; (b) ‘‘sliding and stretching”. From Ref. [92], with permission.

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905‘‘ﬁne structure”; in the latter, the reactants are mixed at the molecular level. Spalding’s interphase-transfer rate is the rate of creationof ﬁne structure in the eddy-dissipation model.
Because of its true two-phase character, Spalding’s two-ﬂuid
model provided the beneﬁt of allowing for two different local densities (for the fresh and burning ﬂuids, respectively). Fluid fragments with different densities react differently to a commonpressure gradient, thus originating a ‘‘non-gradient” (and, sometimes, ‘‘counter-gradient”) diffusion, which conventional models
cannot reproduce, but two-ﬂuid ones can [69,70].
The ‘‘unburned-burning” (or colder–hotter) model of turbulent
combustion is just an instance of the new class of two-ﬂuid turbulence models pursued by Spalding. While ‘‘heavier–lighter” models, such as the combustion one outlined above, allowed for arepresentation of the ‘‘sifting” phenomenon (i.e., fragments of ﬂuid
with a given density moving through fragments with another), inother models a property other than density was selected as aﬂuid-distinguishing property, as indicated in the turbulence modelling section above.
The latest member in this lineage of increasingly sophisticatedmodels which build upon similar principles are the multi-ﬂuid
models (MFMs) of combustion and chemical reaction [105,106].
MFMs are often introduced by Spalding as an evolution of theeddy-break model [107]. While in the EBU model only two thermochemical states are permitted (burned and unburned), in MFMs anarbitrary number of intermediate states are allowed by deﬁning
population deﬁning attributes (PDAs), such as mixture fraction
and/or reactedness. These PDAs are discretized, and several ﬂuids
are allowed to co-exist in the domain, each characterized by aset of values for the PDAs. The multiﬂuid model solves conservation equations for the local mass fraction of each ﬂuid. In the turbulent ﬂow, fragments of these ﬂuids with different compositions
collide (the collision among fragments being an idea originally behind the theories on turbulence of Reynolds and Prandtl); and,after some contact time, separate, thus leading to a certain fractionof intermediate ﬂuids at the interface, which Spalding terms ‘‘theoffspring”. The population distribution is in this way changed.
The theory allows also for ‘‘continuously-varying attributes” inthe ﬂuids, i.e., properties which are not discretised.
The intuitively developed MFMs share many analogies with themathematically derived probability density function (PDF) equation solved with a Monte Carlo method [108,109]. Fluid fragments
are akin to ﬂuid particles; however, while the deﬁning propertiesof the former are discretised as population deﬁning attributes,
the properties of the Monte Carlo particles are allowed to vary continuously. In both approaches, micro-mixing is modelled when two
ﬂuid fragments collide: in the Monte Carlo method, the so-calledmolecular mixing models, which frequently address mixing as aparticle–particle interaction are used; in MFMs, mixing is represented through the rate of mass transfer between the colliding
ﬂuid fragments as they enter in contact with each other, resulting
in the ‘‘offspring”. In MFMs, the PDF can be reconstructed from themass fraction of the ﬂuid population.
CFD has been applied to the prediction of the performance of
combustion processes since the 1960s, and Spalding and the Imperial College team pioneered the use of the technology for combustion calculations. First, Patankar and Spalding [42] devised themilestone GENMIX CFD code for the calculation of parabolic ﬂows,
and during 1969 Spalding incorporated combustion into GENMIX,[45,50], which could simulate laminar and turbulent ﬂame jets,
and was later used for the prediction of laminar ﬂame propagation
by Spalding and Stephenson [110]. Spalding proposed a CFD modelfor Diesel Engines in 1969 [111,112], but even as early as 1967,
Spalding had employed the simple chemical reacting system(SCRS) concept into a CFD code for the simulation of turbulent diffusion ﬂames in furnaces, Pun and Spalding [113]. An underlying

3895

feature of the model is that the combustion process is represented
by a global, one-step, inﬁnitely fast chemical reaction in which fueland oxidant react to form products at stoichometric proportions.
The intermediate species are ignored, and the mass fractions of
the reactants, product, and inert species are expressed as ﬁxed
algebraic relationships in terms of a conserved scalar variable
called the mixture fraction. This approach was to become the basisof all CFD modelling of furnaces and combustors for some time tocome, and it was soon extended by Spalding and co-workers to allow fuel and oxidant to co-exist at the same location, although at
different times, by the use of a presumed probability density function (pdf) and the solution of a transport equation for the varianceof the mixture fraction. Thus, the presumed-pdf model probably
started in CFD with Spalding’s ‘‘two-spike” double-delta-function
presumption [114], and subsequently, Lockwood and Naguib[115] made the ‘‘clipped-Gaussian” presumption; and they were
followed by many others, including Kent and Bilger [116]. The
third type of combustion model pioneered by Spalding’s Imperial
College team was to employ the SCRS with a ﬁnite reaction rate
represented by the minimum of the Arrhenius and Spalding’s
eddy-break-up model in the transport equation for the fuel massfraction. These three combustion models are documented, and applied to a number of furnace arrangements in a paper by Khalil
et al. [95], which reports on the modelling work performed undercontract to the Atomic Energy Research Establishment (AERE), Harwell (UK). This paper dealt with two-dimensional ﬂows, however
three-dimensional furnace ﬂows had, in fact, already been simulated by Patankar and Spalding [117,118], and applications tothree-dimensional gas–turbine combustors were soon to be published by Serag-Eldin [119].
10. Partially parabolic and complex ﬂows
One of the early applications of the three-dimensional parabolic
procedure was to ducts with curvature. In 1972, immediately after
the publication of the parabolic calculation procedure, Pratap5[120], working under the supervision of Spalding, attempted thecomputation of ﬂow and heat transfer in helically coiled pipes. This
work consisted of the study of both the developing and fully developed regions, the fully developed being the spatially asymptotic
state. Flow in coiled pipes is of importance in a variety of applications such as heat exchangers, river bends, curved arteries, pipe elbows, etc. When ﬂuid ﬂows through a curved duct the curvatureof the duct induces secondary ﬂows which are directed to the outside in the central region and towards the inside in the near-wall region. These secondary ﬂows are referred to as Dean vortices[121,122]. In the developing region, the streamwise ﬂow evolves
from a parabolic proﬁle to a distorted proﬁle with a peak on the outside region. This is a result of the momentum transfer from the inner
region to the outer region of the bend.
Prior to 1972, there had been no computations of the developing
region of the helical coil. Pratap et al. [123] published the ﬁrst paperon laminar developing ﬂow in a curved pipe. In this paper, the developing region was solved by the parabolic calculation procedure andthe streamwise proﬁles were successfully compared with the dataof Austin [124]. Subsequently, the heat transfer enhancement wasalso compared with experimental data of Mori and Nakayama[125] and Dravid et al. [126] for Pr = 4.0 and 0.71, respectively. This
paper revealed the ability of the parabolic duct ﬂow procedure to
predict industrially relevant complex ﬂows. This procedure was
next extended from laminar ﬂow to turbulent ﬂow. Here, themomentum and continuity equations for a curvilinear coordinate
system were augmented with the k–e turbulence model [127]. The

5

V.S. Pratap is currently known as S.P. Vanka.

3896

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

above two works considered large radius ratios. However, when theduct curvature is large the parabolic assumption can break down.
Spalding proposed a new procedure, called the ‘semi-elliptic’, or‘partially parabolic’ procedure, which accounted for such crossstream pressure variations in accelerating the streamwise ﬂow. Thisidea was followed up by Pratap and Spalding [128] who demonstrated the possibility of capturing such effects. This procedure
was applied to the ﬂow in a duct with strong curvature, and experimental data taken by Pratap [129] were used for comparison.
Flows in ducts rotating about an axis perpendicular to the duct
axis present another situation with body forces on the ﬂuid.
Majumdar et al. [130] applied the three-dimensional parabolic calculation procedure to solve the governing three-dimensional ﬂow
equations in their parabolic form. The momentum equations contain extra terms due to the centrifugal and Coriolis forces. In a subsequent study, Majumdar and Spalding [131] applied the partiallyparabolic procedure to ﬂow in rotating ducts. Agreement with data
was generally good; however certain discrepancies were observed.
This algorithm recognizes that the three-dimensional scalar pressure ﬁeld is indeed elliptic and that there is ﬂuid recirculation in
the cross-plane of the ﬂow, which demands the allocation of scarcecomputer memory to pressure and velocity components, but since
there is no downstream to upstream ﬂow in the main ﬂow direction, then only the pressure need be stored in memory as the ﬂow
then essentially acts parabolically. Pollard and Spalding [132] took
this idea and developed a partially elliptic method.
Majumdar et al. [133] studied the local and mean heat-transfer
characteristics for air ﬂowing turbulently in the fully developed region of a circular and of a square duct rotating around a horizontal
axis parallel to the duct axis. Results showed that the heat transfer
was increased by rotation. The results were compared with experimental data with good agreement. A prior study by Majumdar[134] also considered radial diffusers that were rotating about aperpendicular axis.
The development of the partially parabolic algorithm was a testament to Spalding’s grasp of the physics of ﬂuid mechanics, and oftransport phenomena in general. The above works on ﬂows withbody forces – duct ﬂows inﬂuenced by rotation and curvature –were path-breaking. The newly developed partially parabolic (aswell as the parabolic) procedures were successfully applied. These
led to a new set of computational procedures to predict complex
three-dimensional ﬂows which, by 1975, were being consideredin the Heat Transfer Section at Imperial College (a ‘complex ﬂow’
then being deﬁned as one with externally applied rates of strain
or in which two or more basic turbulent ﬂows interact [135], as opposed to ‘basic ﬂows’: jets, wakes, channels, boundary layers). One
element of the research focus during this period was to establish
the competency to predict three-dimensional turbulent ﬂows
using the turbulence models, discussed elsewhere in this paper.
In many cases, experimental data were also required, and if theequipment existed in other parts of the world, students were sent
there. A review of the 75 or so papers that emerged from the groupduring the period of 1970–1980 reﬂect the vast number of topics
that engaged Spalding, for example: free and forced turbulent convection, Abdelmeguid and Spalding [136]; ﬂow over and aroundvarious shapes and obstacles, Caretto et al. [137] combustion, Serageldin and Spalding [138]; low-Prandtl number heat transfer, ElHadidy et al. [139], ﬂows in T-junctions, Pollard [140]; two-phaseﬂow, Majumdar et al. [141]. A signiﬁcant feature of many of theseinvestigations was the combination of the development of computational methods and models that better captured the main features of the physics of the phenomena under study.
High performance computers are widely prevalent nowadays;
however, in those days, computing devices were quite modest interms of both memory and speed. Thus, what would today ﬁt intoa small handheld device, then required a huge infrastructure andlarge numbers of support personnel. Restrictions on memory sizeand data-storage did not always permit the desire to solve the Navier–Stokes equations for large-scale problems to be satisﬁed.
Thus, innovative algorithm development was always an integral
part of the CFD community activities. In the early to mid 1970s,
two-dimensional ﬂows were being computed, but on what today
seem like very small numbers of control volumes. With the addition of new hardware with more memory and storage, and thedevelopment of ﬂoating point operations, the idea of solving equations in three-dimensional space and even time could be contemplated; however, with that, there was the need to invent even moreelegant algorithms.
11. Heat exchangers
Patankar and Spalding [142] ﬁrst proposed to perform 3D calculations for ﬂuid ﬂow and heat transfer in shell-and-tube heatexchangers by means of a porous medium or ‘distributed resistance’ analogy. This was essentially a very early example of a multi-scale analysis, where detailed computations around the bafﬂesof the heat exchanger were performed but drag and heat transfer
were ‘modelled’ based on prescribed volumetric friction and heattransfer factors. These factors were to be obtained either analytically, or from detailed ﬁne-scale numerical calculations, LeFeuvre[143], Beale and Spalding [144], or from correlations of experimental data. This approach was quite different from conventional heatexchanger design approaches based on ‘presumed ﬂow’ methods,
such as the Delaware method [145]. As previously mentioned,
the partial elimination algorithm (PEA) proved most beneﬁcial inthis context, since in later applications (see below) the volumeaveraged equations for shell-side, and tube-side ﬂuids closely
mimicked those for two-phase ﬂow and heat transfer. The Patankar–Spalding approach was subsequently adopted and reﬁned bya number of other researchers. A complete discussion of Spalding’s
philosophical approach to heat exchanger analysis may be found in
his articles in the Heat Exchanger Design Handbook [146]. Thissubject has remained one of enduring interest and Spalding continues to lecture on the subject to this day [147].
12. Stress analysis
What numerical method is typically used by engineers to calculate stresses and strains in solid objects? The usual answer is the‘‘Finite Element Method” (FEM), which is routinely applied to thesimulation of structural analysis problems, and forms the basis
for many computer codes. There are, however, numerous multiphysics problems, where solid stresses must be computed at thesame time as, say, turbulent ﬂuid ﬂow with heat and mass transfer,chemical reactions and so forth. Such ﬂuid–solid interaction (FSI)
problems are of much scientiﬁc interest and practical necessity.
Some examples include thermal and centrifugal stresses in a rotating turbine, pressure wave propagation in ﬂuid ﬂow in an elastic
tube such as blood-ﬂow in human arteries, or ﬂow in a printerhead chamber containing a piezoelectric membrane, to name but
a few. In these and other cases, the choice of methodology is lessclear since FEM algorithms are arguably not as sophisticated for
simulation of coupled systems of hydrodynamic equations, compared to the ﬁnite-volume method (FVM), which is well developed
and currently employed in almost all major CFD codes. Moreover,such CFD codes already contain numerous sub-component models
for physical phenomena (turbulence, radiation, chemical reactions,multi-phase ﬂow, etc.) which are necessary constituents for modelling complex FSI problems.
Spalding was among the ﬁrst to propose that the FVM may
readily be adapted to obtain solutions for elasticity problems

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

and, to suggest that it is probably simpler to develop the equations
in a FVM-based code, than to port the large set of existing physicalﬂuid mechanics and heat transfer models from established FVM
based CFD codes over to FEM-based codes. Beale and Elias[148,149] showed how the FVM code PHOENICS could be applied
to stress analysis problems using the analogy between the stream
and stress functions in creeping-ﬂuid-ﬂow and solid mechanics.
That approach was notionally similar to the w–x methodology ofRunchal et al. [43], and was thus limited to two-dimensional planestrain/stress problems. Spalding subsequently considered a displacement, rather than a force analogy which has the major advantage of being quite general in three dimensions.
The governing equations for linear thermo-elasticity are basedon the equation of equilibrium, which combined with the generalized Hooke’s law and the strain–displacement relations, result inNavier’s displacement equations. Spalding originally considered
the displacement equations in the following form,

G div grad u þ f  grad p ¼ 0

ð19Þ

where

p ¼ ðk þ GÞD þ HðT  T 0 Þð20Þ

where u is displacement, and f is a body-force. The ‘pressure’ p isthus proportional to the cubical dilatation, D ¼ divu, where k and
G and are Lamé constants, and H is proportional to the coefﬁcientof thermal expansion. Eq. (19) strongly resembles the (Stokes)
momentum equation for creeping ﬂow. Thus, Spalding [150] ﬁrst
proposed in 1993 the FVM solution of FSI problems, based on thisanalogy between the equations for displacement (solids) and velocity (ﬂuids), with the following modiﬁcation to the usual FVM: inﬂuid regions u = (u, v, w) represents velocity whereas in solid regions, it represents displacement. In the latter, the convection terms
are set to zero. The linear relation between pressure p, dilatation D
and T, Eq. (20), is introduced by the inclusion of a pressure and temperature dependent ‘mass-source’ term. Strains and stresses may
readily be computed from the displacements. Later, Spalding[151] replaced the pressure/dilatation term in Eq. (19) with anisotropic exchange coefﬁcients, C, for the solid-phase-diffusion terms.
He also noted that other forms of the governing equations, forexample based on rotation (the stress-equation counterpart to vorticity) are possible, and may ultimately prove superior to displacements. For the solution algorithm it is possible to use the wellknown SIMPLE algorithm, but in [151] Spalding speculated that
coupled schemes such as the SIVA algorithm [52] may be more
effective in this context. Combined solid-thermal ﬂuids problems,for example, in porous media, should prove readily amenable to thisclass of analysis in the future. The uniﬁcation of stresses and strains
into the general framework of ﬂuid ﬂow and heat transfer withchemical reactions and turbulence, is another logical step forwardin the development of a general-purpose engineering methodology
founded upon physical reasoning.
13. Concentration Heat and Momentum Ltd.
There can be little doubt that Spalding and his associates pioneered the commercial CFD industry by founding the company,
originally known as Combustion Heat and Momentum Ltd. (CHAM)
in 1969. The company, which became Concentration Heat andMomentum Ltd. in 1974, provided commercial CFD services basedon the pioneering technology that had emerged from Spalding’s research group at Imperial College. While the activities of CHAM andthe genesis of the CFD code PHOENICS (see below) are irrevocablyintertwined, well before the advent of PHOENICS, CHAM wasessentially the world’s only commercial CFD company, and remained so until the early 1980s when Creare Inc. began marketing
the FLUENT code, the early origins of which can be traced to

3897

Swithenbank et al. [152]; see also Runchal [153]. Between 1974
and 1980, CHAM developed many application-speciﬁc CFD computer codes for industrial clients across a wide range of industrial
sectors, including the aerospace, automotive, defence, chemical,
environmental, ﬁre and safety, marine, manufacturing and process,nuclear-power, and fossil-fuel power industries. These CFD codes(some of which found their way into the public domain) also provided the means for CHAM to undertake CFD consultancy contracts
aimed at solving practical problems for industry. CHAM’s clientswere for the most part, large engineering companies in the UK,
Europe, and the USA. During this period CHAM scored numerous‘ﬁrsts’, under Spalding’s direction and supervision. For example,
it created the ﬁrst CFD code for the steam generators of the nuclear
industry: for this, CHAM had to devise CFD programs which could
simulate the separate ﬂows of steam and water within the same
space. As another example, CHAM produced the ﬁrst Electrolytic
SmelTER (ESTER) CFD code, for the aluminium-production industry, to simulate multi-anode electrolytic smelters of the Hall-cell
type, in which electromagnetic effects and gravity waves on themolten-metal/electrolyte interface can interact to limit the performance of the equipment. Later on, ESTER was superseded by a
PHOENICS-based product, Rosten [154], that was also used extensively by the industry [155–157]. Coal- and gas-ﬁred furnaces,
heat-exchangers, diesel and petrol (gasoline) engines, power condensers, steam generators, and natural-draught cooling towers
were all among the equipment items for which CHAM produced
the ﬁrst CFD packages.
During the 1970s the partially parabolic method, described earlier, provided CHAM with a more practical route than fully elliptic
codes for some 3D industrial applications. For example, the partially-parabolic CFD code FLASH was employed to compute 3D ﬂow
around ship and submarine hulls using a curvilinear coordinate
system. This work was largely funded by Ishikawa-Harima Heavy
Industries (IHI) (Japan), National Maritime Institute (NMI) (Teddington), and Admiralty Experimental Works (AEW) (Haslar), andresulted in a substantial number of publications by Spalding andco-workers, e.g., Abdelmeguid et al. [158] and Markatos et al.
[159]. As computer technology became more powerful, FLASH
was superseded by an in-house elliptic CFD code, and then ultimately by a PHOENICS-based CFD model, with a general nonorthogonal body-ﬁtted coordinate system, Malin et al. [160]. The
partially parabolic concept also proved useful in turbomachinery
applications, and in 1976 the SEMFT partially-parabolic CFD code
was developed by Singhal and Spalding [161] for predicting 2D
steady, turbulent ﬂow in turbomachinery cascades. This code, supplied to Sulzer Brothers in 1977, was capable of handling incompressible, subsonic, supersonic, or transonic ﬂows. The CATHY3
code was another partially parabolic turbomachinery code thatwas supplied to a number of companies, including the Westinghouse Electric Corporation in 1978. A curvilinear coordinate system was employed for the prediction of steady, 3D turbulent
ﬂow in rotating turbomachine passages. The code was derived
from earlier work by Pratap and Spalding [128] and Majumdar
and Spalding [162]. It was replaced in 1979 by a fully elliptic codeof the same name, Malin et al. [163]. This elliptic CFD code included a cavitation model, and the software was supplied to NASA
Marshall for the prediction of steady, 3D turbulent ﬂow in centrifugal impellers.
As was mentioned earlier, Patankar and Spalding [142] pioneered the use of the distributed-resistance concept in CFD modelling of 3D ﬂow and heat transfer in a shell-and-tube heatexchanger. This approach has been used extensively by CHAM
since the early 1970s for conducting CFD simulations of heat-exchanger, steam-generator, and condenser applications in the defence, marine, nuclear, and power industries. These applications
were funded by numerous companies including Combustion

3898

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

Engineering (USA), Rolls Royce (UK), National Nuclear Corporation(UK), Electrical Power Research Institute (USA), European Atomic
Energy Community (Italy), Northern Engineering Industries Parsons (UK), General Electric Company (UK), and Vattenfall (Sweden). The publications resulting from these and other studies by
Spalding and co-workers were numerous, see, for example, Singhal
et al. [164,165], Al-Sanea et al. [166], Jureindini et al. [167], andMalin [168]. Some applications included the simulation of twophase ﬂow by means of homogeneous, algebraic-slip or two-ﬂuid
modelling via IPSA. Eventually PHOENICS became the vehicle forall such computations at CHAM, but notable pre-PHOENICS CFD
codes developed by Spalding and other CHAM personnel were
ATHOS, CALIPSOS, FOCS, LOBSTER, STELLA, and URSULA.
The origin of many industrial multi-phase computations can betraced to GALA [61,62] and IPSA [60,61], both devised by Spaldingin the mid seventies. These were described earlier in the manuscript, but it should be mentioned here that Spalding’s remarkable
paper [61] also describes different interface-tracking techniques asalternatives to using GALA for the simulation free-surface ﬂows.
One method involves tracking particles, while the other involves
the solution of a conservation equation for the liquid height. Interface tracking techniques were employed successfully both by
CHAM as well as the group at Imperial College in the late 1970s
and early 1980s, see, for example, Maxwell [169], Awn [170], andCastrejon and Andrews [171]. IPSA has made a tremendous impact
on the numerical computation of industrial multi-phase ﬂows. It
formed the basis of the bulk of the two-phase ﬂow simulations
made by CHAM and its clients with PHOENICS and earlier CHAM
codes, Moult et al. [172], Baghdadi et al. [173], Markatos et al.
[174,175], Marchand et al. [176], Al-Sanea et al. [177], Boisson
and Malin [178]. Inevitably, the method later found its way into
competing commercial CFD products.
Spalding was content to employ the algebraic slip model formany suitable multi-phase-ﬂow applications. This model was inuse at CHAM during the late 1970s for steam-generator applications, Singhal et al. [179,180], but its use became much more widespread during the 1980s for use in the chemical and materials
processing industries. Several variants of the model were in use
during this time, mostly restricted to considerations of gravitational and/or centrifugal force ﬁelds, for example: Pericleous andDrake [181] (particulate separation in cyclones); Tacke and Ludwig[182] (inclusion separation in tundishes); Pericleous and Patel[183] (solid suspensions in stirred tank reactors). Eventually, Spalding proposed a more general formulation that employed GALA forthe satisfaction of volumetric continuity. The method was implemented in PHOENICS in 1992. In more recent times, multi-phasecapabilities using algebraic slip modelling began to appear in othergeneral-purpose CFD software.
During the mid 1970s to early 1980s another novel applicationarea was the prediction of ﬁre and smoke movement in buildings.
At that time, the Fire Research Station, Borehamwood UK (FRS)
funded CHAM’s pioneering work in the development and application of the 2D and 3D elliptic MOSIE CFD codes, Markatos et al.
[184]. These special-purpose ﬁre-and-smoke simulation codes
were superseded by a PHOENICS-based model, which was subsequently used by the FRS and later by the Building Research Establishment (BRE), to create the well-known CFD ﬁre-simulation code
JASMINE, Markatos & Cox [185], Kumar et al. [186], which makes
use of an early version of the PHOENICS code as its CFD engine.
In this way, Spalding’s contribution can be seen as providing FRS
and BRE with a basic CFD capability from which they were able
to develop and validate the CFD technology embodied in JASMINE
to such an extent that CFD now forms an accepted component ofﬁre safety engineering.
In the early years at CHAM, combustion modelling formed a signiﬁcant part of the business, and by 1973 the company had created

for the Rocket Propulsion Establishment (Westcott, UK) a very advanced elliptic CFD combustion code called BAFL, see Tatchell and
Spalding [187] and Jensen et al. [188], for predicting ﬂow, combustion and heat transfer in the region immediately downstream of arocket moving at supersonic or subsonic speeds. BAFL was ahead ofits time, and it incorporated a number of sophisticated mathematical models which enabled it to handle compressible, turbulent,recirculating ﬂow with two-way chemical reactions between gaseous species, ﬂuctuations of temperature and concentration, radiative heat transfer in particular wave-length bands (Lockwood andSpalding [46]), and reaction with particulate species in distinct size
ranges (Spalding [189]). Indeed, some of the ﬁrst CFD codes sold by
CHAM to gas–turbine manufacturers were speciﬁcally for combustor simulation. Several of those combustor codes remained in use
for well over a decade, having of course also been signiﬁcantly further developed by their users; and at least one of them entered the
public domain by way of the US Army, enabling competing CFDcode vendors to start business; which they did with alacrity. In
1976, CHAM supplied the NASA Lewis Research Centre with a combustion CFD code using a chemical equilibrium model for characterising hydrocarbon reactions, and a chemical kinetics model for
predicting the oxides of nitrogen, Elghobashi et al. [190]. Thiswas perhaps the ﬁrst CFD model for predicting pollutant formation. Another notable CHAM combustion code of the 1970s was
CORA3, which was delivered to a number of industrial clients. This
CFD combustion code, among others, made use of Spalding’s ‘eddybreak-up’ [93,191] and ‘presumed-pdf’ methods [81,191], which
took account of turbulent ﬂuctuations. The impact of Spalding’s
work on CFD in the power generation industry was mentioned
by Fiveland, see Hanna [192], who started CFD research work inthe boiler manufacturing industry during the 1970s at Babcock
and Wilcox in the USA:‘‘At the time I was aware of Professor Spalding’s pioneering CFD
work in the UK, and we saw value in using this ‘‘new” CFD technology to improve our boiler and burner designs.”
Computers during the 1970s had small memories and slow execution times, and so ingenious programming was needed to enableindustrial problems to be solved. Moreover, there was much skepticism in some circles as to whether computer simulation of complex turbulent and multi-phase ﬂows would ever be of any
practical use. Some time in 1978, Spalding conceived the idea ofa single CFD code capable of handling all ﬂuid-ﬂow processes. Consequently, CHAM abandoned the policy of developing individual
application-speciﬁc CFD codes, and during late 1978 the company
began creating the world’s ﬁrst general-purpose CFD code:
PHOENICS.

14. The PHOENICS code
Many scientists and engineers distinguish themselves by writing a monograph on their ﬁeld of expertise, at some point in theircareers. Spalding however devoted his time and energy instead to apractical piece of CFD software, which in turn helped spawn awhole industry-ﬁeld, and created jobs and wealth, where none previously existed. CFD as a sub-set of computer aided engineering
can thus be directly traced to the developments of the Parabolic,
Hyperbolic Or Elliptic Numerical Integration Code Series, PHOENICS[51,193–195].
PHOENICS was launched commercially in 1981 as a FORTRANinput program and later, in 1984, as an interactive-input CFD code.
Here for the ﬁrst time a single code was to be used for all thermoﬂuids problems. This represented a major philosophical change. Indeed, at the time, there was substantial resistance to the idea that ageneral-purpose CFD code was even possible. Credit for PHOENICS

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

is due also to a number of other people, especially H.I. Rosten, andit would be remiss not to mention the signiﬁcant and ongoing contributions made since the code’s inception by J.C. Ludwig. However, PHOENICS was essentially Spalding’s brainchild, andremains the subject of his attention to the present day. PHOENICS
represented a paradigm shift in the advance of CFD; engineers
were now able to apply the algorithms and closure assumptions
previously developed at Imperial College and elsewhere to a whole
host of real-world industrial problems. Once this step was reached
CFD advanced irreversibly beyond the ivory tower of academia andinto the real world of engineering products and processes. The
requirement to provide a working software product to a client
led to a new level of professionalism, in terms of software design,
memory management/optimization, portability, documentation,and training, not typically found in a university environment. Because of the general-purpose nature of the code, CFD advanced beyond thermo-ﬂuids to a whole range of multi-physics problems, forexample, chemical-vapour deposition to give one example.
In the early 1980s, PHOENICS was the only general-purpose CFD
code available, and represented an embodiment of many of the scientiﬁc principles previously published by Spalding’s group, manyof which were meticulously reproduced in a case-library for execution by the user. The concept of writing a single- and multi-phasecode, for laminar or turbulent ﬂow, with complex chemical kinetics, heat and mass transfer (including thermal radiation) and variable properties was quite unique. The code was designed with acentral core known as ‘EARTH’. This contained all established utilities which were hard-coded and to which the user did not have access (in order to avoid damage by the inexperienced and thus tofacilitate maintenance and support). Open source components
were also provided wherein physical models such as turbulence
closures, wall functions, and solution algorithms were coded andcould be inspected, altered or replaced by the user. ‘Hooks’(threads) were provided at different points in the software forusers to attach their own coding, thus allowing the extension of
the code capabilities and enabling it to be used for many problems
for which it was not speciﬁcally designed. This ‘open-source’ concept, which was pioneered in 1981, was later adopted by someother CFD-code vendors. PHOENICS was ﬁrmly based on the ﬁnite-volume method. At the heart of code was the single-phase
SIMPLEST algorithm [58], a variant of SIMPLE [55,137]. Multiphase ﬂows were accommodated using IPSA [60] (Eulerian) and
GENTRA (Lagrangian). Users could select from a number of linear
equation ‘solvers’ such as point-by-point, slab and whole-ﬁeld (orindeed provide their own), and control relaxation using inertial
or linear methodologies on a per-variable basis. Originally a staggered velocity scheme [53] was employed, although later collocated schemes [196] were added as an option both in rectilinear
and body-ﬁtted coordinates. Boundary conditions and the source
terms themselves were coded as linear source terms, S = GC(V  uP)
in Eq. (11), where G is a geometric factor, C is a source-term ‘coefﬁcient’ and V is a source term ‘value’, over which the user could exert substantial control; thus they were programmable. The code
was based on a segregated scheme with a number of ‘iterations’being performed for each ﬁeld-variable, together with multiple‘sweeps’ over the entire computational ﬁeld, until convergence
was reached and the solver stopped. Spot-values and residuals
were originally printed and, in later versions of the code, plotted
dynamically by the graphical user interface (GUI), to provide theuser with a measure of conﬁdence in the results.
In early versions of PHOENICS, users were provided the opportunity to tailor the code to their needs using FORTRAN either in the
pre-processor, ‘SATELLITE’, or at run-time, in ‘GROUND’. Through
GROUND, the user could in fact access all EARTH variables and arrays which were stored in a single ‘F-array’ by means of functions
and sub-routines which performed a variety of operations. The

3899

nature of the code allowed for extensive programmability, so thatusers could if they wished, completely by-pass all the default solvers and so forth, in favour of their own utilities. It is important toemphasise the enormous number of academic and industrial CFD
computations that have been based on the PHOENICS code. The
proceedings of the First International PHOENICS Users Conferencein 1985 [197] reveals that from the beginning there were a large
number of applications of PHOENICS across a very wide range ofindustries. Subsequent User Conferences, held biennially to date,
together with publications in the PHOENICS Journal of Computational Fluid Dynamics, and in the open literature, served to enlarge
the number of industrial applications to truly immense proportions. Numerous PhD theses were conducted by Spalding’s students at the Computational Fluid Dynamics Unit at Imperial
College [198–215] and elsewhere, based on PHOENICS. While nowadays the notion of post-graduate students utilizing a commercial
code in their research work is scarcely controversial, at the time it
was quite novel.
Some of the code developments included: body-ﬁtted grids(1983); biological ﬂows – bladder, lung, etc. (1987); programmable
input language and conjugate heat transfer (1989); expert-based
system to improve convergence (1991); solid stresses/strains(1992); parallel implementation (1993); and ﬁne-grid embedding(1997). Later, new programming features were introduced whichallowed users to create their own physical models without theneed to resort to FORTRAN. As the user-base grew, so did the library of cases. Over the years, innumerable developments by users,
both industrial and academic were adopted into the code. In addition, special versions of PHOENICS were developed for a wide
range of applications: electrolytic smelters, chemical-vapour deposition, electronic cooling, blast furnaces, cooling towers, external
aerodynamics. Special-purpose PHOENICS User Groups were
established by CHAM.
PHOENICS was the world’s ﬁrst general-purpose CFD code; andby the time of its release in 1981, the underlying scientiﬁc principles had been extensively published. Therefore, once CHAM
showed that a general-purpose CFD code could be created andmarketed proﬁtably, competitors soon followed the lead. The practices and principles embodied in PHOENICS have been much imitated; and there are now numerous such CFD packages in themarket offered by competing vendors. The lineage of most of thesecodes can be traced back either to CHAM itself or to the Imperial
College team from which CHAM emerged, and the inﬂuence ofCHAM and PHOENICS was sometimes affected by staff movement.
Certainly, it can be argued that CHAM’s innovation in creating andmarketing PHOENICS founded a worldwide industry which isbased on the selling of CFD software and services, which although
no ﬁrm ﬁgures are available, now has an annual turnover in excess
of six ﬁgures.
Within a few years of the ﬁrst release of PHOENICS, Spalding
initiated work at CHAM on the creation and marketing of PHOENICS-based, special-purpose products (SPPs) that employed application-speciﬁc GUIs. The ﬁrst of these was ESTER, which wasdeveloped speciﬁcally for the aluminium industry. Other SPPs
quickly followed in the late 1980s and early1990s, for example,
HOTBOX (for electronics cooling), FLAIR (for heating, ventilationand ﬁre and smoke movement in buildings), TACT (for cooling towers), PHOENICS-VWT (for motorsport vehicle aerodynamics) and
PHOENICS-CVD (for chemical-vapour-deposition). Spalding argued
for the application of CFD to CVD reactors as early as 1979 [216],
including a proposal to use embedded ﬁne-grids to resolve moreaccurately the ﬁne structure of the ﬂuid ﬂow and deposition process in the vicinity of the wafer.
PHOENICS-CVD, Heritage [217], remains the most technicallyadvanced of the SPPs, being able to simulate the behaviour of awide range of CVD reactors by modelling ﬂuid ﬂow and heat

3900

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

transfer in a multi-component gas, including both gas-phase(homogeneous) and surface (heterogeneous) chemical reactions,
thermal radiation and plasma effects. During the 1980s and1990s, PHOENICS was heavily used by CHAM and other organizations for metallurgical- and materials-processing applications.
Such applications are often challenging and can involve combinations of ﬂuid-ﬂow, heat and mass transfer, phase change, electromagnetics, free-surface ﬂow, multiple phases, and chemistry. The
successful use of PHOENICS-ESTER by the aluminium industry
has already been mentioned. In addition Alcan International collaborated with CHAM to use IPSA in PHOENICS for simulating themixing and solidiﬁcation of submerged liquid metal jets, see Enright et al. [218]. CHAM also used IPSA to create PHOENICS-based
CFD models for the NASA Lewis Research Centre to simulate solidliquid phase change in binary materials, see Prakash [219,220]. In
the mid-1980s CHAM also began work on developing a PHOENICS-based CFD model of the HIsmelt process, which involves
the direct smelting of iron ore by injecting the iron and coal ﬁnes
directly into a molten bath, plus the subsequent release and burning of the smelting gas in the top space. Staff movement from
CHAM in the early 1990s facilitated the continued reﬁnement of
the model by Rio Tinto (Australia), Davis et al. [221]. PHOENICSbased CFD models for materials processing were developed formany other companies, including Arbed Recherche (Luxembourg),
Centro Sviluppo Materiali (Italy) and Hoogevens (The Netherlands). During the 1980s and early 1990s, PHOENICS was also used
extensively at Massachusetts Institute of Technology by Szekely
and co-workers for the development of CFD models of materialsprocessing operations, A large number of papers were publishedduring this period, for example: Ilegbusi and Szekely [222,223](tundishes and continuous casting), Cartwright et al. [223] (Czochralski systems) Choo and Szkeley [224] (arc welding). and Saluja
et al. [225] (stirred melts). PHOENICS has also been in use at Metallurgical Research Institute AB (Sweden) and the Center for Advanced Study and Research of the National Polytechnic Institute(Mexico) for many years to this day, to conduct research into metals processing, see, for example, Solhed and Jonsson [226] andRodriguez et al. [227].
Cartesian cut-cell methods for representing complex geometries were in use at CHAM for speciﬁc applications during the early
1980s, but it was not until a decade later that Spalding began toinitiate a program of work to equip PHOENICS with a generally
applicable facility called PARSOL (partial solids). The motivation
was to treat imported CAD geometries implicitly and simplify considerably the grid generation associated with using body-ﬁtted
non-orthogonal meshes. PARSOL has never been formally published, although some discussion was given in Spalding [228].
The method employs a background Cartesian grid for the majorityof the ﬂow domain with special treatments being applied at cells
which are cut by solid bodies, thus retaining a boundary-conforming grid. PARSOL is also partially unstructured in that it allows local
grid reﬁnement by means of embedded ﬁne-grid regions. The
method computes the fractional areas and volumes, and employs
a collection of special algorithms for computing interfacial areas,
evaluating wall stresses, and for computing advection and diffusion near solid boundaries, etc. Cartesian cut-cell methods havenow become a recognized alternative to both structured
and unstructured body-conforming meshes, and they are now used in at least three other commercial CFD codes.
During 1993 CHAM implemented into PHOENICS a number of
low-Reynolds-number turbulence models that employed normalwall-distance-based damping functions. It was quickly recognizedthat the use of conventional search-based geometrical proceduresfor determining the wall-distance distribution would be far too
computationally expensive for handling the arbitrary complex
geometries used in industrial applications of PHOENICS. It wasthen that Spalding proposed the ingenious and extremely economical approach of computing approximate values of the wall distance by solving a Poisson-type transport equation for a scalarvariable. Spalding’s proposal was mainly intuitive, being basedon an analogy with thermal radiation between two plates. This differential-equation based wall-distance algorithm was immediately
put to use in PHOENICS for all relevant turbulence closures. The
method was only ever formally presented in a poster session at
the 10th International Heat Transference Conference, Brighton,
by Spalding himself, but thereafter it was exploited, reﬁned andpublished by Tucker [229,230] and others, and eventually themethod or some variant of it found its way into most commercial
CFD codes.
The invention of the wall-distance method emerged simultaneously with Spalding’s proposal for a zero-equation low-Reynolds-number turbulence model, LVEL, to determine theturbulent kinematic viscosity for situations in which ﬂuid ﬂows
through spaces cluttered with many solid objects [231]. In suchcases, the grid density between nearby solids is often too coarse
for any more advanced turbulence model, such as the two-equation k–e model, to be meaningfully employed. The turbulent kinematic viscosity was calculated via Spalding’s law of the wall [232],
which covers the entire laminar and turbulent regimes of the ﬂow.
A key element of this model is the Poisson-based differential equation for wall distance mentioned earlier. The LVEL model was
implemented in PHOENICS, and used to great practical advantage
in electronics-cooling applications, so much so that almost inevitably the model began to appear in other commercial CFD codes, andespecially ones dedicated to the thermal design of electronic
components and systems. Spalding wrote in his one of his papers[231]‘‘The LVEL model is to be regarded as providing a practical solution
to heat-transfer-engineers’ problems, and not as one providing newscientiﬁc insight.”
Certainly, the practical usefulness of the model has been demonstrated by many workers, including Rodgers et al. [233], Dhinsa
et al. [234], and Choi et al. [235]. The latter report that the LVEL
model was as effective as the more advanced two-equation models
for simulating rack-mounted servers, with the big advantage beingthat ‘‘signiﬁcant computation time can be saved (factor of three orhigher), especially when conducting transient CFD simulations orwhen making steady-state runs with different rack settings”.
Spalding’s ability to recognize the need for engineers in industry to get an adequate solution within reasonable computing time
led to the development of the IMMERSOL radiation model. Thismodel was proposed by Spalding in 1994, mainly for economically
representing radiation in congested spaces, as found in electronicscooling applications, and in any other problems where many conducting solids are distributed within the ﬂuid. IMMERSOL, which
stands for IMMERsed SOLids, was implemented into PHOENICS
by CHAM as an economic option for representing radiation in complex geometries for both optically thick and thin media. The model
remains unpublished by Spalding in the open literature, but adescription has been given by Rasmussen [236] and Yang et al.
[237]. The model involves the solution of a diffusion equation for
the radiosity, in which the diffusivity is proportional to the absorption and scattering coefﬁcients of the medium, and the distance between solid walls. The inter-wall distances are computed
automatically from a Poisson equation, as explained earlier in connection with computation of wall distances for turbulence modelling. The radiosity equation was devised by simplifying the six-ﬂux
model of Hamakar and Schuster, see Spalding [238], by ignoring
the directional effects of radiation, but the extension to transparent
media was based on physical intuition. IMMERSOL provides an
economically realisable approximation to the radiative transfer

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905

equation that gives the exact solution for thermal radiation between two inﬁnitely long parallel plates; and in more complex situations, its predictions are always plausible, and of the right orderof magnitude. Successful applications of the model for combustion
applications have been reported by both Yang et al. [237,239] andRasmussen [236], and also by Hien and Istiadji [240] for building
ventilation studies with solar radiation. Thus far, IMMERSOL appears to have made little impact on industrial CFD, other than on
those organizations who use PHOENICS. One exception is foundin the work of Rasmussen [236], who implemented IMMERSOL intoanother commercial CFD code, STAR-CD, via user programming,and then extended the model to handle wavelength dependence.
PHOENICS stood in the vanguard of practical engineering CFD;
for many years it set the standard which the competition attempted to replicate. Inevitably it was imitated, and eventually another CFD code overtook it in terms of sales-volume. Today the CFD
practitioner can choose from a number of code-suites, many ofwhich were written by former students of Spalding at Imperial College, or by entrepreneurs who began their initial careers at CHAM
and left to form new companies. However PHOENICS also still enjoys widespread use world-wide, and especially so in research circles. It is still the only general-purpose code for which the ‘P’(Parabolic) and the ‘H’ (Hyperbolic) options are appropriate. Indeed
the parabolic procedure in PHOENICS was extended by Malin and
Spalding [241] to handle under-expanded free jets discharging into
subsonic and supersonic free streams, under contract to S&C Thermoﬂuids UK, for use in rocket-exhaust-plume software. This produced solutions much more economically than the elliptic
solution procedure. Probably no other commercial code to date
has the extensive programmability available to the user. PHOENICS
thus is Spalding’s own personal, continuously-evolving, living
monograph to CFD.
At the 1995 PHOENICS CFD User Conference in Tokyo, Spalding
made the following light-hearted observation in terms of CFD Tshirt slogans, on the transition of CFD technology from its researchorigins in Imperial College to a tool for industrial analysis anddesign:3901

contribution Brian Spalding has made to progress in the ﬁelds ofmechanical and chemical engineering. Spalding’s physical insights,
mathematical ability and extraordinary work ethic enabled him to
make so many advances in the creation and development of mathematical models to describe turbulence, heat transfer, mass transfer and combustion. His work is characterized by a superb
creativity and an extraordinary vision. His ideas often produced aquantum leap in our understanding and predictive capability,
and opened the door to new scientiﬁc opportunities. He had thedesire and ability to generalize a given concept, formula, or procedure. This is what led to far-reaching consequences from his research contributions. The extraordinary depth of knowledge ismatched by the enormous breadth of the applications of his work.
In addition to his technical work, he was able to communicate his
ideas to others effectively through his clear and lucid lectures andwritings. His communications skills, combined with being ﬂuent inGerman and Russian among other languages, enabled him to havea very large impact around the world. His leadership, vision, passion and single-minded obsession also provided industry withmulti-purpose tools for the practice of CFD and the developmentof a suite of tools and services for commercial application. These
contributions have had a profound impact in the world of engineering software. He was always a strong advocate for international cooperation and exchange, and had a major inﬂuence on
the lives and careers of so many of his students and colleagues,
among whom we are fortunate to count ourselves. We have tried
here to acknowledge some of the achievements of this remarkable
man and his remarkable life. Doubtless there will be more in the
future.
Acknowledgements
The authors gratefully acknowledge the assistance of Colleen
King in the preparation and collection of material used to complete
this paper and of Micha Wolfshtein for his comments on an earlyversion.
References

DecadeSloganWhat was happening

1960s
1970s
1980s
1990s
2000s

We think we can
We can!
You can too
Everyone canEveryone will

Imperial College researchFirst applications
Start of code-vending industry
Commercial over-sell
Web-based services

Spalding’s pioneering role in industrial CFD was summarized by F.C.
Lockwood, an original member of the Imperial College team, as follows [242]:‘‘Spalding was able to combine mathematical ﬂair with a copious
amount of physical intuition which led to a scientiﬁc methodology
which was universally appreciated by the engineer at the ‘‘sharpend” of real industry problems,”
and‘‘Commercial CFD software developers are not always noted fortheir generosity of acknowledgement so there is merit in remarking
that virtually all of the commercial codes being marketed today
owe much to the pioneering work of Spalding and co-workers.”

15. Conclusion
The preparation of this article has been a labour of love on thepart of the authors. It is intended as a small tribute to the very large[1] D.B. Spalding, The Combustion of Liquid Fuels, Ph.D. thesis, Cambridge
University, 1951.
[2] T. von Kármán, Über Laminare Und Turbulente Reibung, Z. Agnew. Math.
Mech. 1921 (1) (1921) 233–252.
[3] G.N. Kruzhilin, Issledovanie Teplovogo Pogranichnogo Sloia, Zh. Tekh. Fiz. 2(1936) 183.
[4] E.R.G. Eckert, V. Lieblein, Berechnung Des Stoffüberganges an Einer Ebenen,
Längs Angeströmten Oberﬂäche Bei Großem Teildruckgefälle, Forschung auf
dem Gebiete des Ingenieurwesens 16 (2) (1949) 33–42.
[5] Y.B. Zeldovich, V.A. Frank-Kamenetskii, Teoriia Teplovogo Rasprostraneniia
Plameni, Zh. Khim. Fiz. (J. Phys. Chem. (USSR)) 12 (1) (1938) 100–105.
[6] Y.B. Zeldovich, Selected Works of Yakov Borisovich Zeldovich (A. Granik, E.
Jackson, Trans.), Princeton University Press, 1992–1993.
[7] N.N. Semenov, Teplovaia Teoriia Goreniia I Vzryvov: Iii. Teoriia Normalnogo
Rasprostraneniia Plameni, Usp. Fiz. Nauk (Progress Phys. Sci. (USSR)) 24 (4)(1940) 433–486.
[8] D.B. Spalding, Some Fundamentals of Combustion, Butterworth, London,
1955.
[9] D.B. Spalding, Convective Mass Transfer; an Introduction, Edward, Arnold,
London, 1963.
[10] D.B. Spalding, A standard formulation of the steady convective mass transfer
problem, Int. J. Heat Mass Transfer 1 (1960) 192–207.
[11] W.M. Kays, M.E. Crawford, B. Weigand, Convective Heat and Mass Transfer,
fourth ed., McGraw-Hill, New York, 2005.
[12] A.F. Mills, Mass Transfer, Prentice Hall, Upper Saddle River, NJ, 2001.
[13] D.B. Spalding, A single formula for the law of the wall, ASME J. Appl. Math. 28(3) (1961) 455–458.
[14] D.B. Spalding, Heat transfer to a turbulent stream from a surface with a stepwise discontinuity in wall temperature, in: Proceedings of the International
Developments in Heat Transfer, Part II, 1961, pp. 439–446.
[15] J. Kestin, L.N. Persen, Application of Schmidt’s method to the calculation of
Spalding’s function and the skin-friction coefﬁcient in turbulent ﬂow, Int. J.
Heat Mass Transfer 5 (1962) 143–152.
[16] G.O. Gardner, J. Kestin, Calculation of the Spalding function over a range of
Prandtl numbers, Int. J. Heat Mass Transfer 6 (1963) 289–299.

3902

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905[17] A.G. Smith, V.L. Shah, The calculation of wall and ﬂuid temperatures for theincompressible turbulent boundary layer, with arbitrary distribution of wallheat ﬂux, Int. J. Heat Mass Transfer 5 (1962) 1170–1189.
[18] D.B. Spalding, Contribution to the theory of heat transfer across a turbulent
boundary layer, Int. J. Heat Mass Transfer 7 (1964) 743–761.
[19] C.L.V. Jayatillake, The inﬂuence of Prandtl number and surface roughness on
the resistance of the laminar sub-layer to momentum and heat transfer, Prog.
Heat Mass Transfer 1 (1969) 193–329.
[20] S.S. Kutateladze, A.I. Leontev, Turbulent Boundary Layers in Compressible
Gases (D.B. Spalding, Trans.), Arnold Press, London, 1964.
[21] D.B. Spalding, A Uniﬁed Theory of Friction, Heat Transfer, and Mass Transfer
in the Turbulent Boundary Layer and Wall Jet, ARC No. 25 925, Her Majesty’s
Stationery Ofﬁce, London, 1964.
[22] D.B. Spalding, A Uniﬁed Theory of Friction, Heat Transfer and Mass Transfer in
the Turbulent Boundary Layer and Wall Jet ARC CP No. 829, Her Majesty’s
Stationery Ofﬁce, London, 1965.
[23] F.P. Ricou, D.B. Spalding, Measurements of entrainment by axisymmetricalturbulent jets, J. Fluid Mech. 11 (1) (1961) 21–32.
[24] D. Coles, The law of the wake in the turbulent boundary layer, J. Fluid Mech. 1(1956) 191–226.
[25] M.R. Head, Entrainment in the Turbulent Boundary Layer, ARC R&M No. 3152,
Her Majesty’s Stationery Ofﬁce, London, 1960.
[26] D.B. Spalding, The Kinetic-Energy Deﬁcit Equation of the Turbulent Boundary
Layer, AGARDograph No. 97, Part 1, 1965, pp. 191–244.
[27] M.P. Escudier, W.B. Nicoll, The entrainment function in turbulent-boundarylayer and wall-jet calculations, J. Fluid Mech. 25 (2) (1966) 337–366.
[28] H. Fellouah, A. Pollard, Reynolds number effects within the development
region of a turbulent round free jet, Int. J. Heat Mass Transfer, in press,
doi:10.1016/j.ijheatmasstransfer.2009.03.029.
[29] J. Rotta, Uber die Theorie der Turbulenten Grenzschichten Ström. Forsch. Nr.
1. Transl As: On the Theory of the Turbulent Boundary Layer, NACA TM 1344,
1953, 1950.
[30] E. Truckenbrodt, Ein Qudraturverfahren Zur Berechnung der Laminaren Und
Turbulenten Reibungschicht Bei Ebener Und Rotationssymmetrischer
Strömung, Ing. Arch. 20 (1952) 211.
[31] B.R. Morton, G.I. Taylor, J.S. Turner, Turbulent gravitational convection frommaintained and instantaneous sources, Proc. R. Soc. Lond. A 234 (1956) 1–23.
[32] S.V. Patankar, D.B. Spalding, A calculation procedure for heat transfer byforced convection through two-dimensional uniform-property turbulent
boundary layers on smooth impermeable walls, in: Proceedings of the 3rdInternational Heat Transfer Conference, Chicago, vol. II, 1966, pp. 50–63.
[33] A. Thom, An Investigation of Fluid Flow in Two Dimensions, ARC R&M No.
1194, Her Majesty’s Stationery Ofﬁce, 1928.
[34] O.R. Burggraf, Analytical and numerical studies of steady separated ﬂows, J.
Fluid Mech. 24 (1966) 113–151.
[35] R. Courant, E. Isaacson, M. Rees, On the solution of non-linear hyperbolic
differential equations by ﬁnite-differences, Commun. Pure Appl. Math. 5(1952) 243–255.
[36] A.K. Runchal, Three ﬁnite difference methods for Navier–Stokes equations, in:Proceedings of the 2nd Applied Mechanics Conference, University ofStrathclyde, Glasgow, 1967.
[37] A.K. Runchal, D.B. Spalding, M. Wolfshtein, The Numerical Solution of the
Elliptic Equations for the Transport of Vorticity, Heat and Matter in Two
Dimensional Flows, SF/TN/2, Imperial College, Mech. Eng. Dept., 1967.
[38] A.K. Runchal, M. Wolfshtein, A Finite-Difference Procedure for the Integrationof the Navier–Stokes Equations, SF/TN/1, Imperial College, Mech. Eng. Dept.,
1966.
[39] A.K. Runchal, M. Wolfshtein, A Fortran IV Computer Program for the Solutionof Steady-State, Two-Dimensional Equations of Motion, Energy andConcentrations, SF/TN/10, Imperial College, Mech. Eng. Dept., 1967.
[40] A.D. Gosman, W.M. Pun, A.K. Runchal, D.B. Spalding, M. Wolfshtein, Heat andMass Transfer in Recirculating Flows, Academic Press, 1969.
[41] G.D. Mallinson, G. de Vahl Davis, Three-Dimensional Natural Convection in a
Box – A Numerical Study, J. Fluid Mech. 83 (1) (1977) 1–31.
[42] S.V. Patankar, D.B. Spalding, A ﬁnite-difference procedure for solving theequations of the two-dimensional boundary layer, Int. J. Heat Mass Transfer
10 (1967) 1389–1412.
[43] A.K. Runchal, D.B. Spalding, M. Wolfshtein, Numerical solution of the elliptic
equations for transport of vorticity, heat, and matter in two-dimensionalﬂow, Phys. Fluids 12 (1969) 21–28.
[44] S.V. Patankar, D.B. Spalding, Heat and Mass Transfer in Boundary Layers,
Morgan-Grampian, London, 1967.
[45] S.V. Patankar, D.B. Spalding, Heat and Mass Transfer in Boundary Layers: A
General Calculation Procedure, second ed., Intertext Books, London, 1970.
[46] F.C. Lockwood, D.B. Spalding, Prediction of Turbulent Reacting Duct Flow with
Signiﬁcant Radiation, Report No. CCK/TN/A/6, Mech. Eng. Dept., Imperial
College, London, 1971.
[47] D.B. Spalding, Convergence and accuracy of three ﬁnite difference schemes
for a two dimensional conduction and convection problem, Int. J. Numer.
Methods Eng. 4 (1972) 540–550.
[48] A.K. Runchal, Convergence and accuracy of three ﬁnite difference schemes fora two dimensional conduction and convection problem, Int. J. Numer.
Methods Eng. 4 (4) (1972) 541–550.
[49] B.P. Leonard, A stable and accurate convective modeling procedure based on
quadratic upstream interpolation, Comput. Methods Appl. Mech. Eng. 19 (1)(1979) 59–98.

[50] D.B. Spalding, GENMIX: A General Computer Program for Two-Dimensional
Parabolic Phenomena, 1st ed., Pergamon Press, Oxford, New York, 1977.
[51] D.B. Spalding, PHOENICS a general-purpose computer program for multidimensional one- and two-phase ﬂow, Math. Comput. Simul. XXIII (1981)
267–276.
[52] L.S. Caretto, R.M. Curr, D.B. Spalding, Two numerical methods for threedimensional boundary layers, Comput. Methods Appl. Mech. Eng. 1 (1) (1972)
39–57.
[53] F.H. Harlow, J.E. Welch, Numerical calculation of time-dependent viscousincompressible ﬂow of ﬂuid with free surface, Phys. Fluids 8 (12) (1965)
2182–2189.
[54] A.J. Chorin, The numerical solution of the Navier–Stokes equations for an
incompressible ﬂuid, Bull. Am. Math. Soc. 73 (1967) 928–931.
[55] S.V. Patankar, D.B. Spalding, A calculation procedure for heat, mass andmomentum transfer in three-dimensional parabolic ﬂows, Int. J. Heat Mass
Transfer 15 (10) (1972) 1787–1806.
[56] J.P. Van Doormaal, J.P. Raithby, Enhancements of the SIMPLE method forpredicting incompressible ﬂuid ﬂows, Numer. Heat Transfer B Fundam. 7 (2)(1984) 147–163.
[57] S.V. Patankar, A calculation for two-dimensional elliptic problems, Numer.
Heat Transfer 4 (1981) 409–426.
[58] D.B. Spalding, Mathematical Modelling of Fluid-Mechanics, Heat-Transfer andChemical-Reaction Processes: A Lecture Course, HTS/80/1, Computational
Fluid Dynamics Unit, Imperial College, University of London, London, 1980.
[59] D.B. Spalding, IPSA 1981, New Developments and Computed Results, HTS/81/
2, Imperial College, London, 1981.
[60] D.B. Spalding, Numerical computation of multi-phase ﬂow and heat transfer,
in: C. Taylor (Ed.), Recent Advances in Numerical Methods in Fluids, Pineridge
Press, 1980, pp. 139–168.
[61] D.B. Spalding, The calculation of free-convection phenomena in gas–liquid
mixtures, in: D.B. Spalding, N. Afgan (Eds.), Heat Transfer & Turbulent
Buoyant Convection, vol. 2, Hemisphere, Washington, DC, 1977.
[62] D.B. Spalding, A Method for Computing Steady and Unsteady Flows
Possessing Discontinuities of Density, Report 910/2, CHAM Ltd.,
Wimbledon, London, 1974.
[63] R.T. Lahey, M. Lopez de Bertodano, O.C. Jones, Phase distribution in complex
geometry conduits, Nucl. Eng. Des. 141 (1993) 177–201.
[64] M.R. Malin, D.B. Spalding, A two-ﬂuid model of turbulence and its application
to heated plane jets and wakes, Physicochem. Hydrodyn. 5 (5-6) (1984) 339–
361.
[65] O.J. Ilegbusi, D.B. Spalding, A two-ﬂuid model of turbulence and its
application to near-wall ﬂows, Physicochem. Hydrodyn. 9 (1-2) (1987)
127–160.
[66] O.J. Ilegbusi, D.B. Spalding, Application of a two-ﬂuid model of turbulence toturbulent ﬂow in conduits and free shear layers, Physicochem. Hydrodyn. 9(1–2) (1987) 161–181.
[67] O.J. Ilegbusi, D.B. Spalding, Prediction of ﬂuid ﬂow and heat transfer
characteristics of turbulent shear ﬂows with a two-ﬂuid model of
turbulence, Int. J. Heat Mass Transfer 32 (4) (1989) 767–774.
[68] D.B. Spalding, Towards a Two-Fluid Model of Turbulent Combustion in Gases
with Special Reference to the Spark Ignition Engine, CFD/82/18,
Computational Fluid Dynamics Unit, Imperial College, London, 1982.
[69] D.B. Spalding, The two-ﬂuid model of turbulence applied to combustion
phenomena, AIAA J. 24 (6) (1986) 876–884.
[70] D.B. Spalding, A turbulence model for buoyant and combusting ﬂows, Int. J.
Numer. Methods Eng. 24 (1) (1987) 1–23.
[71] S.J. Kline, M.V. Morkovin, G. Sovran, D.J. Cockrell, Computation of turbulent
boundary layers, in: AFOSR-IFP-Stanford Conference, 1968.
[72] D.B. Spalding, Heat transfer from turbulent separated ﬂows, J. Fluid Mech. 27(1967) 97–109.
[73] L. Prandtl, Uber Ein Neues Formelsystem Für die Ausgebildete Turbulenz,
Nachr. Akad. Wiss. Göttingen (1945) 6–19.
[74] B.E. Launder, D.B. Spalding, The numerical computation of turbulent ﬂows,
Comput. Methods Appl. Mech. Eng. 3 (2) (1974) 269–289.
[75] A.N. Kolmogorov, Equations of turbulent motion of an incompressible ﬂuid,
Izv. Akad. Nauk. SSSR Ser. Phys. V (1–2) (1942) 56–58.
[76] J. Rotta, Statistische Theorie Nichthomogener Turbulenz I and II, Z. Phys.
129, 131 (1951) 547–572, 51–77 (Translated into English by W. Rodi as
Imperial College, Mech. Eng. Dept. Technical Notes Nos. TWF/TN/38 andTWF/TN/39).
[77] D.B. Spalding, The calculation of the length scale of turbulence in someturbulent boundary layers remote from walls, in: T. Irvine (Ed.), Progress in
Heat & Mass Transfer, Pergamon Press, Oxford, 1967, pp. 255–266.
[78] W. Rodi, D.B. Spalding, A two-parameter model of turbulence, and its
application to free jets, Wärme Stoffübertrag. 3 (2) (1970) 85–95.
[79] K.H. Ng, D.B. Spalding, Some applications of turbulence models for boundary
layers near walls, Phys. Fluids 15 (1) (1972) 20–30.
[80] D.B. Spalding, The Prediction of Two-Dimensional, Steady, Turbulent Elliptic
Flows, Heat and Mass Transfer in Turbulent Flows with Separated Regions
and Measurement Techniques, Paper Delivered at International Center for
Heat and Mass Transfer Seminar, Herceg-Novi, Yugoslavia, Imperial College,
Mech. Eng. Dept. Report No. EF/TN/A/16, London, 1969.
[81] D.B. Spalding, Concentration ﬂuctuations in a round turbulent free jet, Chem.
Eng. Sci. 26 (1) (1971) 95–107.
[82] D.B. Spalding, The kW Model of Turbulence, Technical Note No. TM/TN/A16,
Mech. Eng. Dept., Imperial College, London, 1971.

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905[83] D.B. Spalding, A two-equation model of turbulence, commemorative lecture
for Prof. F. Bosnajakovic, VDI Forschungsheft 549 (1972) 5–16.
[84] K. Hanjalic, B.E. Launder, A Reynolds-stress model of turbulence and its
application to asymmetric shear ﬂows, J. Fluid Mech. 52 (1972) 609–638.
[85] W.P. Jones, B.E. Launder, The prediction of laminarization with a twoequation model of turbulence, Int. J. Heat Mass Transfer 15 (1972) 301–314.
[86] D.C. Wilcox, Reassessment of the scale-determining equation for advanced
turbulence models, AIAA J. 26 (1988) 1299–1310.
[87] F.R. Menter, Two-equation eddy-viscosity turbulence models for engineering
applications, AIAA J. 32 (1994) 1598–1605.
[88] B.E. Launder, A. Morse, W. Rodi, D.B. Spalding, The prediction of free shear
ﬂows – a comparison of the performance of six turbulence models, in:Proceedings of the NASA Conference on Free Shear Flows, Langley VA, 1972.
[89] M.R. Malin, D.B. Spalding, The prediction of turbulent jets and plumes by theuse of the k-w model of turbulence, Physicochem. Hydrodyn. 5 (2) (1984)
153–198.
[90] O.J. Ilegbusi, D.B. Spalding, An improved version of the k-w model of
turbulence, ASME J. Heat Transfer 107 (1) (1985) 63–69.
[91] D.B. Spalding, The Vorticity-Fluctuations (kW) Model of Turbulence; Early
Papers, Report CFD/82/17, Computational Fluid Dynamics Unit, Imperial
College, London, 1982.
[92] D.B. Spalding, The ESCIMO Theory of Turbulent Combustion, HTS/76/13,
Department of Mechanical Engineering, Imperial College, London, 1976.
[93] D.B. Spalding, Mixing and chemical reaction in steady conﬁned turbulent
ﬂames, in: Proceedings of the 13th Symposium (International) on
Combustion, Pittsburgh, PA, 1971, pp. 649–657.
[94] D.B. Spalding, The spread of turbulent ﬂames in conﬁned ducts, in:Proceedings of the 11th Symposium (International) on Combustion, 1967,
pp. 807–815.
[95] E.E. Khalil, D.B. Spalding, J.H. Whitelaw, The calculation of local ﬂow
properties in two-dimensional furnaces, Int. J. Heat Mass Transfer 18 (6)(1975) 775–791.
[96] D.B. Spalding, Mathematical-models of turbulent ﬂames, Combust. Sci.
Technol. 13 (1–6) (1976) 3–25.
[97] D.B. Spalding, Some general-comments on papers concerned with theory ofreactive turbulence, Combust. Sci. Technol. 13 (1–6) (1976) 245–255.
[98] D.B. Spalding, Development of the eddy-breakup model of turbulent
combustion, in: Proceedings of the 16th Symposium (International) on
Combustion, Pittsburgh, PA, 1977, pp. 1657–1663.
[99] D.B. Spalding, Computer Modelling Techniques for Laminar and Turbulent
Combustion, HTS/77/12, Department of Mechanical Engineering, Imperial
College, London, 1977.
[100] D.B. Spalding, A general theory of turbulent combustion, J. Energy 2 (1)(1978) 16–23.
[101] D.B. Spalding, The inﬂuences of laminar transport and chemical kinetics onthe time-mean reaction rate in a turbulent ﬂame, in: Proceedings of the 17th
Symposium (International) on Combustion, 1979, pp. 431–440.
[102] A.R. Kerstein, Linear-eddy modeling of turbulent transport. II: Application toshear layer mixing, Combust. Flame 75 (1989) 397–413.
[103] D.B. Spalding, Developments in the IPSA procedure for numericalcomputation of multiphase-ﬂow phenomena with interphase slip, unequal
temperatures, etc, in: T.M. Shih (Ed.), Proceedings of the 2nd International
Conference on Numerical Properties and Methodologies in Heat Transfer,
Maryland University, 1981, pp. 421–436.
[104] B.F. Magnussen, On the structure of turbulence and a generalized eddy
dissipation concept for chemical reaction in turbulent ﬂow, in: 19th AIAA
Aerospace Science Meeting, St. Louis, Missouri, 1981.
[105] D.B. Spalding, Multi-ﬂuid models of turbulence progress and prospects, in:Proceedings of the CFD96 4th Annual Conference of the CFD Society ofCanada, Ottawa, 1996, pp. 27–56.
[106] D.B. Spalding, Turbulent mixing and chemical reaction; the multi-ﬂuid
approach, in: International Symposium on the Physics of Heat Transfer in
Boiling and Condensation, Moscow, 1997.
[107] D.B. Spalding. Connexions between the multi-ﬂuid and ﬂamelet models of
turbulent combustion, A lecture delivered at the University of Leeds, England.
<http://www.cham.co.uk/phoenics/d_polis/d_lecs/mfm/ﬂamelet.htm>, 1998
(Accessed 3 March 2008).
[108] C. Dopazo, Probability density function approach for a turbulent axisymmetricheated jet centerline evolution, Phys. Fluids 22 (1) (1979) 20–30.
[109] S.B. Pope, Monte-Carlo method for the PDF equations of turbulent reactive
ﬂow, Combust. Sci. Technol. 25 (1981) 159–174.
[110] D.B. Spalding, P.L. Stephenson, Laminar ﬂame propagation in hydrogen andbromine mixtures, Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 324 (1558)(1971) 315–337.
[111] D.B. Spalding, Predicting the performance of diesel engine combustion
chambers, Proc. IMechE 184 (1969) 241–251.
[112] S.V. Patankar, D.B. Spalding, Mathematical models of ﬂuid ﬂow and heattransfer in furnaces, in: Proceedings of the 4th Symposium of Flames &Industry; Predictive Methods for Industrial Flames, 1972, pp. 1–5.
[113] W.M. Pun, D.B. Spalding, A procedure for predicting the velocity andtemperature distribution in a conﬁned, steady, turbulent, gaseous diffusionﬂame, in: Proceedings of the International Astronautics Federation Meeting,
Belgrade, 1967, pp. 3–21.
[114] D.B. Spalding, Turbulent, Physically-Controlled Combustion Processes, Report
RF/TN/A/4, Mech. Eng. Dept. Imperial College, University of London, London,
1971.

3903

[115] F.C. Lockwood, A.S. Naguib, The prediction of the ﬂuctuations in theproperties of free, round-jet, turbulent, diffusion ﬂames, Combust. Flame 24(1975) 109–124.
[116] J.H. Kent, R.W. Bilger, The prediction of turbulent diffusion ﬂame ﬁelds andnitric oxide formation, in: Proceedings of the 16th Symposium (International)
on Combustion Pittsburgh, PA, 1976, pp. 1643–1656.
[117] S.V. Patankar, D.B. Spalding, A computer model for three-dimensional ﬂow in
furnaces, in: Proceedings of the 14th International Symposium on
Combustion, 1972, pp. 605–614.
[118] S.V. Patankar, D.B. Spalding, Simultaneous predictions of ﬂow pattern andradiation for three-dimensional ﬂames, in: N. Afgan, J. Beer (Eds.), Heat
Transfer in Flames, Hemisphere, Washington, 1974, pp. 73–94.
[119] M.A. Serag-Eldin, The Numerical Prediction of the Flow and Combustion
Processes in a Three-Dimensional Combustion Chamber, Ph.D. thesis,
Imperial College, London, 1977.
[120] V.S. Pratap, Flow and Heat Transfer in Curved Ducts, Ph.D. thesis, Imperial
College, University of London, London, 1975.
[121] W.R. Dean, Note on the motion of ﬂuid in a curved pipe, Phil. Mag. 14 (1927)
208–223.
[122] W.R. Dean, Stream-line motion of a ﬂuid in a curved pipe, Phil. Mag. 5 (1928)
673–695.
[123] S.V. Patankar, V.S. Pratap, D.B. Spalding, Prediction of laminar ﬂow and heattransfer in helically coiled pipes, J. Fluid Mech. 62 (1974) 539–551.
[124] L.R. Austin, The Development of Viscous Flow within Helical Coils, Ph.D.
thesis, University of Utah, 1971.
[125] Y. Mori, W. Nakayama, Study on forced convective heat transfer in curved
pipes (1st report, laminar region), Int. J. Heat Mass Transfer 8 (1965) 67–82.
[126] A.N. Dravid, K.A. Smith, E.W. Merrill, P.L.T. Brian, Effect of secondary ﬂuid
motion on laminar ﬂows and heat transfer in helically coiled tubes, AIChE J.
17 (1971) 1017–1122.
[127] S.V. Patankar, V.S. Pratap, D.B. Spalding, Prediction of turbulent ﬂow incurved pipes, J. Fluid Mech. 67 (1975) 583–595.
[128] V.S. Pratap, D.B. Spalding, Fluid ﬂow and heat transfer in three-dimensional
duct ﬂows, Int. J. Heat Mass Transfer 19 (10) (1976) 1183–1188.
[129] V.S. Pratap, D.B. Spalding, Numerical computations of ﬂow in curved ducts,
Aeronaut. Q. XXVI (1975) 219–228.
[130] A.K. Majumdar, V.S. Pratap, D.B. Spalding, Numerical computation of ﬂow inrotating ducts, ASME J. Fluids Eng. 99 (1977) 148153.
[131] A.K. Majumdar, D.B. Spalding, A numerical investigation of three-dimensional
ﬂows in a rotating duct by a partially-parabolic procedure, in: Proceedings of
the ASME Paper No. 77-WA/FE-7, 1977.
[132] A. Pollard, D.B. Spalding, The Prediction of the three-dimensional turbulent
ﬂow ﬁeld in a ﬂow splitting tee-junction, Comput. Math. Appl. Mech. Eng. 13(1978) 293–306.
[133] A.K. Majumdar, W.D. Morris, D. Skiadaressis, D.B. Spalding, Heat transfer inrotating ducts, Mech. Eng. Bull. (Durgapur) 8 (4) (1978) 87–95.
[134] A.K. Majumdar, D.B. Spalding, Numerical investigation of ﬂow in rotatingradial diffusers, Mech. Eng. Bull. (Durgapur) 8 (1–2) (1977) 19–28.
[135] P. Bradshaw. Complex turbulent ﬂows, in: W.T. Koiter (Eds.), Proceedings ofthe 14th Iutam Congress, Theoretical and Applied Mechanics, North Holland,
Delft, 1976, pp. 103–113.
[136] A.M. Abdelmeguid, D.B. Spalding, Turbulent ﬂow and heat transfer in pipes
with buoyancy effects, J. Fluid Mech. 94 (2) (1979) 383–400.
[137] L.S. Caretto, A.D. Gosman, S.V. Patankar, D.B. Spalding, Two Calculation
procedures for steady, three-dimensional ﬂows with recirculation, in:Proceedings of the 3rd International Conference on Numerical Methods in
Fluid Mechanics, vol. 2, 1973, pp. 60–68.
[138] M.A. Serageldin, D.B. Spalding, Computations of three-dimensional gasturbine combustion chamber ﬂows, Trans. ASME J. Eng. Power 101 (3) (1979)
326–336.
[139] F. El-Hadidy, F. Gori, D.B. Spalding, Further results on the heat transfer to lowPrandtl number ﬂuids in pipes, Numer. Heat Transfer A Appl. 5 (1) (1982)
107–117.
[140] A. Pollard, The numerical calculation of partially elliptic ﬂows, Numer. Heat
Transfer A Appl. 2 (3) (1979) 267–278.
[141] A.K. Majumdar, A.K. Singhal, D.B. Spalding, Numerical modeling of wet
cooling towers. I. Mathematical and physical models, ASME J. Heat Transfer
105 (4) (1983) 728–735.
[142] S.V. Patankar, D.B. Spalding, A calculation procedure for the transient andsteady-state behavior of shell-and-tube heat exchangers, in: N. Afgan, E.U.
Schlünder (Eds.), Heat Exchangers: Design and Theory Sourcebook, Scripta
Book Company, Washington, DC, 1974, pp. 155–176.
[143] R.F. Le Feuvre, Laminar and Turbulent Forced Convection Processes
through in-Line Tube Banks, Ph.D. thesis, Imperial College, University of
London, 1973.
[144] S.B. Beale, D.B. Spalding, Numerical study of ﬂuid ﬂow and heat transfer intube banks with stream-wise periodic boundary conditions, Trans. CSME 22(4A) (1998) 394–416.
[145] R.W. Serth, Process Heat Transfer: Principles and Applications, Academic
Press, 2007.
[146] G.F. Hewitt (Ed.), Heat Exchanger Design Handbook, vol. 1, Begell House, New
York, 2008.
[147] D.B. Spalding, Numbers count: SFT analysis of heat-transfer equipment. The
James P. Hartnett Lecture (Paper KN-15), in: G. de Vahl Davis, E. Leonardi(Eds.), Proceedings of the IHTC-13: International Heat Transfer Conference,
Sydney, 2006.

3904

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905[148] S.B. Beale, S.R. Elias, Stress distribution in a plate subject to uniaxial loading,
PHOENICS J. Comput. Fluid Dyn. Appl. 3 (3) (1990) 255–287.
[149] S.B. Beale, S.R. Elias, Numerical Solution of Two-Dimensional Elasticity
Problems by Means of a SIMPLE-Based Finite-Difference Scheme, TR-LT-020,
National Research Council, Ottawa, 1991.
[150] D.B. Spalding, Simulation of ﬂuid ﬂow, heat transfer and solid deformation
simultaneously, in: NAFEMS Conference No. 4, Brighton, 1993.
[151] D.B. Spalding, Extending the frontiers of computational heat transfer, Paper
Cht-08-001, in: Proceedings of the CHT-08 ICHMT International Symposium
on Advances in Computational Heat Transfer, Marrakech, 2008.
[152] J. Swithenbank, A. Turan, P.G. Felton, D.B. Spalding, Fundamental Modellingof Mixing, Evaporation and Kinetics in Gas Turbine Combustors, Report No.
HIC 320, Department of Chemical Engineering & Fuel Technology, Universityof Shefﬁeld, 1979.
[153] A.K. Runchal, Brian spalding: CFD & reality, Int. J. Heat Mass Transfer, in press,
doi:10.1016/j.ijheatmasstransfer.2009.03.058.
[154] H.I. Rosten, The Mathematical Foundation of the ESTER Computer Code,CHAM TR/84, Wimbledon, 1982.
[155] W.E. Wahnsiedler, Hydrodynamic modeling of commercial Hall-Heroult cells,Light Metals (1987) 269–287.
[156] V. Potocnik, Modeling of metal–bath interface waves in Hall-Heroult cells
using ESTER/PHOENICS, in: Proceedings of the Light Metals, 1989, pp. 227–235.
[157] M. Segatz, D. Vogelsang, C. Droste, P. Baekler, Modeling of transient magnetohydrodynamic phenomena in Hall-Heroult cells, Light Metals (1993) 361–
368.
[158] A.M. Abdelmeguid, N.C. Markatos, K.A. Muraoka, D.B. Spalding, A method ofpredicting three-dimensional turbulent ﬂows around ships’ Hulls, in:Proceedings of the International Symposium on Ship Viscous Resistance,
SSPA, Gothenborg, Sweden, vol. 3, 1978, pp. 1–24.
[159] N.C. Markatos, M.R. Malin, D.G. Tatchell, Computer analysis of threedimensional turbulent ﬂows around ships’ Hulls, Proc. IMechE 194 (24)(1980) 239–248.
[160] M.R. Malin, H.I. Rosten, D.B. Spalding, D.G. Tatchell, Application of PHOENICS
to ﬂow around ship’s Hulls, in: L. Larsson, C. S. Ohlsson (Eds.), Proceedings ofthe 2nd International Symposium on Ship Viscous Resistance, Gothenburg,
Sweden, 1985, pp. 16:1–16:23.
[161] A.K. Singhal, D.B. Spalding, A Two-Dimensional Partially-Parabolic Procedure
for Axial-Flow Turbomachinery Cascades, ARC R&M No. 3807, Her Majesty’s
Stationery Ofﬁce, London, 1978.
[162] A.K. Majumdar, D.B. Spalding, A Parametric Computation of Flows andPerformance of Rotating Radial Impellers, ARC R&M 3810, Her Majesty’s
Stationery Ofﬁce, London, 1977.
[163] M.R. Malin, H.I. Rosten, D.G. Tatchell, Three-dimensional computations ofﬂows in centrifugal pumps and compressors; in performance prediction ofcentrifugal pumps and compressors, in: Proceedings of the 25th Annual
International Gas Turbine Conference and Exhibit and 22nd Annual Fluids
Engineering Conference, New Orleans, LA, 1979, pp. 33–45.
[164] A.K. Singhal, D.B. Spalding, Mathematical modelling of multi-phase ﬂow andheat transfer in steam generators, in: Proceedings of the 2nd Multi-phase
Flow Transport: Fundamentals, Reactor-Safety Applications, 1980, pp. 373–
406.
[165] A.K. Singhal, L.W. Keeton, A.J. Prekwas, J.S. Weems, D.B. Spalding, ATHOS – A
Computer Program for Thermal Hydraulic Analysis of Steam Generators,
1982.
[166] S.A. Al-Sanea, N. Rhodes, D.G. Tatchell, T.S. Wilkinson, A computer model for
detailed calculation of the ﬂow in power-station condensers, in: Proceedings
I. Chem. E. Symposium Series No. 75, Condensers: Theory & Practice, 1983,
pp. 70–88.
[167] R.H. Jureidini, M.R. Malin, M.J. Lord, K.K. Yau, A three-dimensional model for
power condenser design, in: Proceedings of the 4th International PHOENICS
User Conference, Miami, FL, 1991.
[168] M.R. Malin, Modelling ﬂow in an experimental marine condenser, Int.
Commun. Heat Mass Transfer 24 (5) (1997) 597–608.
[169] T.T. Maxwell, Numerical Modelling of Free-Surface Flows, Ph.D. thesis,
Imperial College, University of London, 1977.
[170] A.G. Awn, Numerical Modelling of Flows with Moving Interfaces, Ph.D. thesis,
Imperial College, University of London, London, 1979.
[171] A. Castrejon, M.A. Andrews. A procedure for calculating moving-interface
ﬂows with PHOENICS, in: C.A. Brebbia, S.A. Orszag (Eds.), Proceedings of the
1st International PHOENICS Users Conference, Dartford, Published in Lecture
Notes in Engineering, 18 (Springer-Verlag, London, 1985), pp. 433–443.
[172] A. Moult, V.S. Pratap, D.B. Spalding, Finite-difference predictions of waves instratiﬁed gas–liquid ﬂows, in: Proceedings of the 1st International
Conference on Physicochemical Hydrodynamics, Washington, DC, 1978.
[173] A.H.A. Baghdadi, H.I. Rosten, A.K. Singhal, D.B. Spalding, D.G. Tatchell,
Numerical modelling of two-phase ﬂow with slip, in: Proceedings of the
ICHMT seminar on two-phase ﬂows, Dubrovnik, Yugoslavia, 1978.
[174] N.C. Markatos, A. Moult, P.J. Phelps, D.B. Spalding, The calculation of steady,three-dimensional, two-phase ﬂow and heat transfer in steam generators, in:
N. Afgan, F. Durst, G. Tsiklauri (Eds.), Proceedings of the ICHMT Seminar on
Two-Phase Flows, Dubrovnik, Yugoslavia, vol. 1, 1978, pp. 485–502.
[175] N.C. Markatos, S.M. Rawnsley, D.B. Spalding, Heat transfer during a smallbreak loss-of-coolant accident in a pressurized water reactor-a parametricstudy for a 4-in lower-plenum break, Int. J. Heat Mass Transfer 27 (8) (1984)
1379–1394.

[176] E.O. Marchand, A.K. Singhal, D.B. Spalding, Numerical modelling of two-phaseﬂow in steam generators, in: 2nd Multi-phase Flow & Heat Transfer
Symposium, Miami, FL, 1978.
[177] S.A. Al-Sanea, N. Rhodes, T.S. Wilkinson, Mathematical modelling of twophase condenser ﬂows, in: Proceedings of the BHRA 2nd International
Conference on Multi-Phase Flow, London, 1985, pp. 169–182.
[178] N. Boisson, M.R. Malin, Numerical prediction of two-phase ﬂow in bubble
columns, Int. J. Numer. Methods Fluids 23 (1996) 1289–1310.
[179] A.K. Singhal, D.B. Spalding, Mathematical Formulation and Solution
Procedure of URSULA2 Computer Program, H3201/1, CHAM, 1979, 1978.
[180] A.K. Singhal, L.W. Keeton, D.B. Spalding, Predictions of thermal hydraulics of a
PWR steam generator by using the homogeneous, the two-ﬂuid, and thealgebraic-slip model, in: Proceedings of the AIChE Symposium Series No.199,
vol. 76, 1980, pp. 45–55.
[181] K.A. Pericleous, S.N. Drake, An algebraic-slip model for multi-phaseapplications, numerical simulation of ﬂuid ﬂow and heat/mass transfer
processes, in: C.A. Brebbia, S.A. Orszag (Eds.), Lecture Notes in Engineering,
Springer- Verlag, London, 1986.
[182] K.H. Tacke, J.C. Ludwig, Steel ﬂow and inclusion separation in continuous
casting tundishes, Steel Res. 58 (1987) 262–270.
[183] K.A. Pericleous, M.K. Patel, The modelling of tangential and axial agitators inchemical reactors, Physicochem. Hydodyn. 8 (2) (1987) 105–123.
[184] N.C. Markatos, M.R. Malin, G. Cox, Mathematical modelling of buoyancyinduced smoke ﬂow in enclosures, Int. J. Heat Mass Transfer 25 (1982) 63–75.
[185] N.C. Markatos, G. Cox, Hydrodynamics & heat transfer in enclosurescontaining a ﬁre source, Physicochem. Hydrodyn. 5 (1) (1984) 53–66.
[186] S. Kumar, N. Hoffmann, G. Cox, Some validation of JASMINE for ﬁres in
hospital wards, in: C.A. Brebbia, S.A. Orszag (Eds.), Lecture Notes inEngineering, Springer-Verlag, 1986, pp. 159–169.
[187] D.G. Tatchell, D.B. Spalding, The Rocket Based-Flow Program – BAFL, CHAM
Technical Report 640/1, CHAM Ltd., Wimbledon, 1973.
[188] D.E. Jensen, D.B. Spalding, D.G. Tatchell, A.S. Wilson, Computation of
structures of ﬂames with recirculating ﬂow and radial pressure gradients,
Combust. Flame 34 (1979) 309–326.
[189] D.B. Spalding, Notes on the PARMIX Program, Report BL/TN/A/35, Imperial
College, Mech. Eng. Dept., London, 1970.
[190] S. Elghobashi, D.B. Spalding, S.K. Srivatsa, Prediction of Hydrodynamics andChemistry of Conﬁned Turbulent Methane–Air Mixtures with Attention to
the Formation of Oxides of Nitrogen, NASA CR-135170, CHAM TR 951/2,
1977.
[191] D.B. Spalding, Turbulence Models & Their Experimental Veriﬁcation: – 12.
Turbulence Effects on Combustion, Report HTS/73/27, Imperial College,
London, 1973.
[192] K. Hanna, Visions of the Future: CFD in the Power Generation Industry,
FLUENT News (Fall Edition), 2004.
[193] D.B. Spalding, PHOENICS a General-Purpose Computer Program for MultiDimensional One- and Two-Phase Flow, HTS/81/11, Computational Fluid
Dynamics Unit, Imperial College, University of London, 1981.
[194] D.B. Spalding, Four Lectures on the PHOENICS Code, CFD/82/5, Computational
Fluid Dynamics Unit, Imperial College, University of London, 1982.
[195] D.B. Spalding, PHOENICS 1984. A Multi-Dimensional Multi-Phase GeneralPurpose Computer Simulator for Fluid Flow, Heat Transfer and Combustion,
CFD/84/18, Computational Fluid Dynamics Unit, Imperial College, Universityof London, 1984.
[196] C.M. Rhie, W.L. Chow, Numerical study of the turbulent ﬂow past an airfoil
with trailing edge separation, AIAA J. 21 (11) (1983) 1525–1532.
[197] N.C. Markatos, D.G. Tatchell, M. Cross, N. Rhodes, Numerical simulation of
ﬂuid ﬂow and heat/mass transfer processes, in: C.A. Brebbia, S.A. Orszag(Eds.), Proceedings of the 1st International PHOENICS Users Conference,
Dartford, London (1985), Published in Lecture Notes in Engineering, vol. 18,
Springer-Verlag, 1986.
[198] G.W. Carroll, Experimental and Theoretical Investigation of ThreeDimensional Turbulent Mixing in Jets & Ducts, Ph.D. thesis, Imperial
College, University of London, London, 1980.
[199] P.A. Tennant, A Numerical Analysis of the Inﬂuence of Internal Nozzle
Geometry on Diesel Fuel Injection, Ph.D. thesis, Imperial College, Universityof London, 1981.
[200] J.O. Ilegebusi, A Revised Two-Equation Model of Turbulence, Ph.D. thesis,
Imperial College, University of London, 1983.
[201] A. Castrejon, Unsteady Free-Convection Heat Transfer inside and Annular
Cavity, Ph.D. thesis, Imperial College, University of London, 1984.
[202] H.Q. Qin, Computer and Water-Model Simulation of Flow through Poppet
Valves, Ph.D. thesis, Imperial College, University of London, 1984.
[203] S.M. Lo, Development of a New Interphase Friction Law for Air–Water
Flows in Inclined Pipes, Ph.D. thesis, Imperial College, University of London,
1985.
[204] G.T. Robert, Shock Tube of Convective Heat Transfer, Ph.D. thesis, Imperial
College, University of London, 1985.
[205] M.J. Andrews, Turbulent Mixing by Rayleigh–Taylor Instability, Ph.D. thesis,
Imperial College, University of London, 1986.
[206] L. Jun, Computer Modelling of Flow with Free Surface, Ph.D. thesis, Imperial
College, University of London, 1986.
[207] M.R. Malin, Turbulence Modelling for Flow and Heat Transfer in Jets, Wakes
and Plumes, Ph.D. thesis, Imperial College, University of London, 1986.
[208] J.Z.Y. Wu, Studies in Unsteady Flame Propagation, Ph.D. thesis, Imperial
College, University of London, 1986.

V. Artemov et al. / International Journal of Heat and Mass Transfer 52 (2009) 3884–3905[209] S.T. Xi, Transient Turbulent Jets of Miscible and Immiscible Fluids, Ph.D.
thesis, Imperial College, University of London, 1986.
[210] K.A. Fenech, Analysing Aspects of the Performance of an Iron Blast Furnace,
Ph.D. thesis, Imperial College, University of London, 1987.
[211] C.A. Miller, An Evaluation of the Performance of the PHOENICS Code, Ph.D.
thesis, Imperial College, University of London, 1988.
[212] D. Radosavljevik, The Numerical Simulation of Direct-Contact NaturalDraught Cooling Tower Performance under the Inﬂuence of Cross-Wind,
Ph.D. thesis, Imperial College, University of London, 1988.
[213] F. Villasenor, The Kelvin–Helholtz Instability in a Stratiﬁed Shear Flow, Ph.D.
thesis, Imperial College, University of London, 1988.
[214] N. Fueyo, Two-Fluid Models of Turbulence for Axisymmetrical Jets andSprays, Ph.D. thesis, Imperial College, University of London, 1989.
[215] S.B. Beale, Fluid Flow and Heat Transfer in Tube Banks, Ph.D. thesis, Imperial
College, University of London, 1993.
[216] D.B. Spalding, Numerical Modelling of CVD Reactors, HTS/79/6, Mech. Eng.
Dept., Imperial College, University of London, London, 1979.
[217] J.R. Heritage, PHOENICS-CVD: a code for the design and development ofchemical vapour deposition equipment and processes, PHOENICS J. Comput.
Fluid Dyn. Appl. 7 (4) (1994) 165–180.
[218] P.G. Enright, L. Katgerman, J.C. Ludwig, S.C. Rogers, Mixing and solidiﬁcationof a turbulent liquid jet in a co-ﬂowing stream, Int. J. Numer. Methods Eng. 24(1987) 231–249.
[219] C. Prakash, Two-phase model for binary solid-liquid phase change, part I,
Numer. Heat Transfer B Fundam. 18 (1990) 131–145.
[220] C. Prakash, Two-phase model for binary solid-liquid phase change, part II,
Numer. Heat Transfer B Fundam. 18 (1990) 147–167.
[221] M.P. Davis, R.J. Dry, M.P. Schwarz, Flow simulation of the HIsmelt process, in:Proceedings of the 3rd International Conference on CFD in the Minerals &Process Industries, Melbourne, 2003, pp. 306–311.
[222] O.J. Ilegbusi, J. Szekely, Fluid Flow and tracer dispersion in shallow tundishes,
Steel Res. 59 (9) (1988) 399–405.
[223] R. Cartwright, J.O. Ilegbusi, J. Szekely, A comparison of order-of-magnitudeand numerical analyses of ﬂow phenomena in Czochralski and magnetic
Czochralski systems, J. Crystal Growth 94 (1989) 321–333.
[224] R.T.C. Choo, J. Szekely, Vaporization kinetics and surface temperature in amutually coupled spot gas tungsten arc weld and weld pool, Weld. J. 71 (3)(1992) 77–93.
[225] N. Saluja, J.O. Ilegbusi, J. Szekely, On the calculation of the electromagnetic
force ﬁeld in the circular stirring of metallic melts, J. Appl. Phys. 68 (11)(1990) 5845–5850.
[226] H. Solhed, L. Jonsson, An investigation of slag ﬂotation and entrapment in acontinuous-casting tundish using ﬂuid-ﬂow simulations, sampling andphysical metallurgy, Scand. J. Metall. 32 (2003) 15–32.
[227] E. Rodriguez, A.H. Castillejos, F.A. Acosta, Experimental and numericalinvestigation of ﬂuid ﬂow and mixing in pachuca tanks, Metall. Mater.
Trans. B 38B (2007) 641–656.

3905[228] D.B. Spalding, CAD to SFT, with Aeronautical applications, Plenary Lecture, in:Proceedings of the 38th Israel Annual Conference on Aerospace Sciences,1998.
[229] P.G. Tucker, Assessment of geometric multilevel convergence robustness anda wall distance method for ﬂows with multiple internal boundaries, Appl.
Math. Modell. 22 (1998) 293–311.
[230] P.G. Tucker, Differential equation-based wall distance computation for DES
and RANS, J. Comput. Phys. 190 (1) (2003) 229–248.
[231] D. Aganofer, L. Gan-Li, D.B. Spalding, The LVEL turbulence model for
conjugate heat transfer at low Reynolds numbers, in: Proceedings of the
EEP Application of CAE/CAD to Electronic Systems, ASME International
Mechanical Engineering Congress and Exposition, Atlanta, GA, 1996.
[232] D.B. Spalding, A new analytical expression for the drag of a ﬂat plate valid forboth the turbulent and laminar regimes, Int. J. Heat Mass Transfer 5 (12)(1962) 1133–1138.
[233] P. Rodgers, J. Lohan, V. Eveloy, C.M. Fager, J. Rantala, Validating numericalpredictions of component thermal interaction on electronic printed circuitboards in forced convection air ﬂows by experimental analysis, Adv. Electron.
Packag. 1 (1999) 999–1008.
[234] K.K. Dhinsa, C.J. Bailey, K.A. Pericleous, Turbulence Modelling and Its Impact
on CFD Predictions for Cooling of Electronic Components, in: Proceedings of
the Ninth Intersociety Conference Thermal and ThermomechanicalPhenomena in Electronic Systems, 2004.
[235] J. Choi et al., A CFD-based tool for studying temperatures in rack-mounted
servers, IEEE Trans. Comput. 57 (8) (2008) 1129–1142.
[236] N.B.K. Rasmussen, The composite radiosity and gap (CRG) model of thermalradiation, in: Proceedings of the 6th European Conference on Industrial
Furnaces and Boilers (INFUB-6) 2002 Conference, Estoril, Lisbon, Portugal,
2002. (Also published as Danish Gas Technology Centre Report No. CO201,
Hørsholm, Demark, 2003.)[237] Y. Yang, R.A. de Jong, M. Reuter, Use of CFD to Predict the Performance of aheat treatment furnace, in: Proceedings of the 4th International Conference
on CFD in the Oil and Gas, metallurgical and process Industries, Trondheim,
Norway, 2005, pp. 1–9.
[238] D.B. Spalding, Idealisations of Radiation, HTS/80/1, Mech. Eng., Dept.,
Imperial College, University of London, 1980.
[239] Y. Yang, R.A. de Jong, M. Reuter, CFD prediction for the performance of a heat
treatment furnace, Prog. Comput. Fluid Dyn. Int. J. 7 (2–4) (2007) 209–218.
[240] W.N. Hien, A.D. Istiadji, Effects of external shading devices on daylighting andnatural ventilation, in: Proceedings of the 8th International IBPSA
Conference, Eindhoven, The Netherlands, 2003, pp. 475–482.
[241] M.R. Malin, D.B. Spalding, Extensions to the PHOENICS Parabolic Solver forunder-Expanded Jets, Report C/4366/1 & 2, CHAM Ltd., London, 1998.
[242] F.C. Lockwood, T. Abbas, N.H. Kandamby, V. Sakthitharan, CFD experience on
industrial combustors, Prog. Comput. Fluid Dyn. Int. J. 1 (1/2/3) (2001) 1–13.