Simultaneous Prediction of Solid stress, Heat Transfer and
Fluid flow by a Single Algorithm
By Brian Spalding
Invited Lecture presented at 2002 ASME Pressure Vessels and Piping
Conference,
August 48, 2002, Vancouver, British Columbia, Canada
updated March 2004
abstract or contents or
update
Abstract
 It is often believed that fluidflow and solidstress
problems must
be solved by different methods and different computer
programs.
 This is not true, if the solidstress problems are formulated
in terms of displacements.
 The lecture exemplifies and explains how both displacements and
velocities can be calculated at the same time.
next
This is not the only common fallacy that will be exposed.
Specifically:
 It is often believed that thermal radiation between solids
separated by nonparticipating gases must be computed by
viewfactor or raytracing techniques.
But:
 These become hopelessly expensive when many solids are
present, as in a furnished room or an engine compartment; and
 it is not true because the IMMERSOL conductivitytype
method handles the problem inexpensively and well.
next
 It is often believed that convective heat transfer at low
Reynolds number, which prevails for example in electronics
equipment, requires an advanced 2ormoredifferential equation
model for its prediction.
However:
 these are also computationally expensive; and
 the nodifferentialequation LVEL model does just as well at
much lower cost.
next
 Both kinds of turbulence model require knowledge of the
"distance from the wall", and sometimes "distance between walls";
and it is often believed that complicated trigonometrical computations
must be conducted so as to acquire it.
However:
 these are expensive as well as complicated; and
 they are also not needed, because the easytosolve
LTLS equation provides all that is needed at
negligible cost.
next
To summarise, four misconceptions are to be dispelled, namely
that:
 solidstress analysis needs a separate program;
 solidtosolid radiation needs viewfactors or ray tracing;
 lowRe convection needs (e.g.) LamBremhorst;
 walldistance computation requires trigonometry.
next or back
Contents
update
 The problem
 Its essential nature
 Practical occurrence
 The conventional solution
 A better solution
 A multiphysics example
 Stresses resulting from radiation, conduction and
convection
 Vector and contour plots
 How the stress calculations were performed
next or back or
contents
 The mathematics of the method
 Similarities between the equations for displacement and velocity
 Deduction of the associated stresses and strains
 The "SIMPLE" algorithm for the computation
of displacements

More details of the equations
 Details of the auxiliary models
 IMMERSOL, for radiation
 WGAP, WDIS and LTLS, for radiation and turbulence
 LVEL, for turbulence
 Discussion of achievements and future
prospects
 References
next or back or
contents
1. The problem
(a) Its essential nature
It is frequently required to simulate fluidflow and
heattransfer
processes in and around solids which are, partly as a consequence of the
flow, subject to thermal and mechanical stresses.
Often, indeed, the stresses are the major concern, while
the fluid and heat flows are of only secondary interest.
next, back or contents
(b) Practical occurrence
Engineering examples of fluid/heat/stress interactions include:
 gasturbine blades under transient conditions;
 "residual stresses" resulting from casting or welding;
 thermal stresses in
nuclear reactors during emergency shutdown;
 manufacture of bricks and ceramics;
 stresses in the cylinder blocks of diesel engines;
 the failure of steelframe buildings during fires.
next, back or contents
(c) The conventional solution
It has been customary for two computer codes to be used for
the solution of such problems, one for
the fluid flow and the other for the stresses
Iterative interaction between the two codes is then employed, often
with considerable inconvenience.
next, back or contents
(d) A better solution
It is however possible for fluid flow, heat flow and
solid deformation, and the interactions between them, all to be
calculated at the same time.
The method of doing so exploits the similarity between the
equations governing velocity (in fluids) and those governing
displacement (in solids).
In the present lecture, the results of such a calculation
will be shown first.
The explanation of how it was conducted will then follow.
next, back or contents
2. A multiphysics example
(a) Description (PHOENICS library case, s400):
The task is to calculate the stresses in radiationheated
solids cooled by air.
20 deg C air
 80 deg C
 V /////// hot radiating wall ///////////
 
 duct > exit
 
// copper // cavity //copper //
 ? temperature ?
////////////// aluminium /////////////////

next, back or contents
Details of the calculation are:
 The Reynolds number (based on the inflow velocity and horizontal duct
width) is 1000.
Therefore the
LVEL model
is used for simulation of the
turbulence.
 The radiative heat transfer is represented by the
conductiontype
IMMERSOL model,
which is:
 economical and
 fairly accurate
for such situations.
The absorptivity of the air is taken as 0.01 per meter;
the scattering coefficient as 0.0;
and the solid surface emissivity as 0.9 .
next, back or contents
 Both LVEL and IMMERSOL make use of the distributions of:
 distance from
the wall (WDIS) and
 distance between walls (WGAP), both of which are
calculated by solving a scalar equation for the
 LTLS variable.
 The stresses within the metals result primarily from the differences in their
thermalexpansion coefficients. namely:
 2.35 e5 for aluminium, and
 1.67 e5 for copper.
next, back or contents
(b) Vector and contour plots
Vectors
The vectors displayed in
Fig.1 show simultaneously:
 the nature of the airflow pattern (upper part of the domain), and
 the displacements of the solid material (lower part of the domain).
Both are calculated at the same time. Of course, the two sets of vectors
have different scales, and indeed dimensions. namely m/s for velocity and m
for displacements.
The displacement vectors are shown on their own in
Fig. 1a
next, back or contents
The solids are supposed to be confined by a stiffwalled box, but are
allowed to slide relative to its walls. This is why the displacement
vectors are vertical near the confiningbox walls.
They are however not allowed to slide relative to each other; this is
what causes the concentrations of stress at their contact surfaces.
next, back or contents
 Stresses and strains
The stresses in the x (horizontal) and y (vertical) directions are
displayed in
Fig.2 and
Fig.3
respectively.
They have been deduced from the strains shown in
Fig.4 for the xdirection , and in
Fig.5 for the ydirection,
next, back or contents
The strains have been deduced from the displacements by
differentiation.
The displacements are however not here displayed as contour plots because these do not
easily display small variations.
That is the end of the stressstrain results.
There will now be
shown some of the other variables which had to be computed.
next, back or contents
 Temperature fields
Fig. 6 displays contours
of temperature in the air and the solid, and reveals that:
 the air is heated by contact with:
 the 80degreeCelsius top wall, and
 the metal blocks, which have been receiving heat by radiation
from the top wall;
 temperature differences within the highconductivity solids are too small
to be discerned visually.
next, back or contents
It is interesting to compare Fig. 6 with
Fig. 7.
This displays
the distribution within the air space of:
 the "radiation temperature", T3, which:
 is computed by IMMERSOL, and
 is defined as the temperature which would be taken up
by a probe which was affected only
by radiation.
Obviously, and understandably, T3 and TEM1 have very different
values, unless the absorptivity is very great (as in solids).
next, back or contents
The solid temperature influences the stresses and strains, of course,
primarily through the agency of the temperaturedependent
thermalexpansion distribution.
However, its variations with position, within a single material, are too
slight to be revealed by a contour diagram, as inspection of
Fig. 8 will reveal.
next, back or
contents
 Radiationflux contours
The IMMERSOL model, of which the solution of the T3 equation is the
major feature, enables the radiant heat fluxes in the coordinate
directions to be established by postprocessing.
The results are displayed in
Fig. 9 for the
xdirection. and by
Fig. 10 for the
ydirection.
The values and patterns displayed, if studied and interpreted in physical
terms, will be found to be plausible.
Where calculation by hand is easy, namely for the parallel surfaces,
they will be found to be correct.
next, back or
contents
 Contours of auxiliary quantities used by IMMERSOL
A crucial feature of the IMMERSOL model is its use of the distribution of
the "distance between the walls", WGAP.
This quantity, which has an
easilyunderstood meaning when the walls are near, and nearly parallel, is
computed from the solution of the "LTLS" equation;
this will be explained later in the lecture.
The distributions of these two quantities are shown by
Fig. 11
for the former, and by
Fig. 14
for the latter.
next, back or
contents
It will be seen that WGAP has a uniform value in the region of
between the top of the duct and the tops of the upper metal slabs,
between which the actual distance is 0.008 meters.
Further, it has approximately twice this value near the convex
corners; and it becomes zero in the concave corners.
next, back or
contents
 Contours of auxiliary quantities used in the fluidflow
calculation
The flow field was calculated by means of the LVEL turbulence
model,
which makes use of the walldistance (WDIS) field.
This, like WGAP, is also
derived from the LTLS distribution.
The contours of WDIS are displayed in
Fig. 13.
which exhibits:
 the expected maximum of 0.004 between the parallel
horizontal walls, and
 a somewhat greater value near the cavity,
where the true distance from the wall depends on the direction in
which it is measured.
next, back or
contents
LVEL, like IMMERSOL, is a "heuristic" model, by which is meant that
it is incapable of rigorous justification, but is nonetheless
useful.
WDIS is calculated once for all, at the start of the computation.
From it, and from the developing velocity distribution, the
evolving distribution of ENUT, the effective (turbulent) viscosity
is derived.
The resulting contours of ENUT are shown in
Fig. 14.
Since the laminar viscosity is of the order of 1.e5 m**2/s, it is
evident that turbulence raises the effective value, far from the
walls, by an order of magnitude.
next, back or
contents
(c) How the stresses were calculated
 As will be shown below, the equations governing the displacements are
very similar to those governing the velocities.
 The CFD code PHOENICS, like many others, can calculate velocities in
fluids; but this
ability is not
needed in the solid region; so such codes are usually idle there.
 However, PHOENICS can be "tricked" into calculating what it "thinks"
are velocities everywhere; whereas what it actually calculates in the
solid regions are displacements.
 The details of the "trickery" now follow.
next, back or
contents
3. The mathematics of the method
(a) Similarities between the equations for displacement and velocity
The similarities already referred to are here described for only one
cartesian direction, x; but they prevail for all three directions.
next, back or
contents
 The xdirection displacement, U, obeys the equation:
where:
 Te = local temperature measured above that of the unstressed
solid in the zerodisplacement condition, multiplied by
the thermalexpansion coefficient;
 D = [d/dx]* U + [d/dy]* V + (d/dz]* W
which is
called the "dilatation";
 Fx = external force per unit volume in xdirection;
 V and W = displacements in y and z directions;
 C1, C2 and C3 are functions of Young's modulus and
Poisson's ratio.
next, back or
contents
 When the viscosity is uniform and the Reynolds number is low, so that
convection effects are negligible,
the xdirection velocity, u, obeys the equation:
[del**2]* u  [d/dx]* [ p*c1 ] + fx*c2 = 0 ,
where
 p = pressure,
 fx = external force per unit volume in xdirection,
 c1 = c2 = the reciprocal of the viscosity.
next, back or
contents
Notes:
 The two equations are now set one below the other, so that they
can be easily compared:
 The equations can thus be seen to become identical if:
 p*c1 = D*C1  Te*C3
which implies:
D = [p*c1 + Te*C3 ] / C1
 and
fx * c2 = Fx * C2
next, back or
contents
 The expressions for C1, C2 and C3 are:
 C1 = 1/(1  2*PR)
 C2 = 2*(1 + PR) / YM
where
 PR = Poisson's Ratio. and
 YM = Young's Modulus
and
 C3 = 2 *(1 + PR)/(1  2*PR)
next, back or
contents
 A solution procedure designed for computing velocities will
therefore in fact compute the displacements if:
 the convection terms are set to zero within the solid
region: and

a suitable linear relation between:
 D ( ie [d/dx]* U + ...)
and
 p
is introduced by inclusion of a
pressure and temperaturedependent "masssource" term.
next, back or
contents
(b) Deduction of the associated stresses and strains
The strains (ie extensions ex, ey and ez) are
obtained from differentiation of the
computed displacements.
Thus:
ex = [d/dx]* U
ey = [d/dx]* V
ez = [d/dx]* W
next, back or
contents
Then the corresponding:
 normal stresses, sx, sy, sz, and
 shear stresses tauxy, tauyz, tauzx,
are obtained from the strains by way of equations such as:
sx = {YM / (1  PR**2)} * {ex + PR*ey} and
tauxy = {YM / (1  PR**2)} * {0.5 * (1  PR)*gamxy}
where:
 gamxy = [d/dy]*U  [d/dx]*V
next, back or
contents
(c) The "SIMPLE" algorithm for the computation of displacements
PHOENICS employs (a variant of) the "SIMPLE" algorithm of Patankar &
Spalding (1972) for computing velocities from pressures, under a massconservation
constraint.
Its essential features are:
 All the velocity equations are solved first, with the
current values of p.
 The consequent errors in the massbalance equations are computed.
 These errors are used as sources in equations for
corrections to p.
 The corresponding corrections are applied, and the process is
repeated.
next, back or
contents
All that it is necessary to do, in order to solve
for
displacements simultaneously, is, in solid regions, to treat the
dilatation D as the masssource error and to ensure
that p obeys the above linear relation to it.
Therefore a CFD code based on SIMPLE can be made to solve the displacement equations
by:
 eliminating the convection terms (ie setting Re low); and
 making D linearly dependent on p and
temperatureT.
The "staggered grid" used as the default in PHOENICS proves to be extremely
convenient for soliddisplacement analysis; for the velocities and
displacements are stored at exactly the right places in relation to
p.
next, back or
contents
4. Details of the auxiliary models
(a) IMMERSOL: summary
 The solved differential equation is:
div( effective_conductivity * T3 ) + source = 0
 effective conductivity =
0.75 * sigma * T3**3 / (abso + scat + 1/WGAP)
 source = abso * sigma * ( T1**4  T3**4 )
 in solids, abso = large, so T3 > T1
 surface resistances account for nonunity emissivities
next, back or
contents
Notes:
 The main novelty is the inclusion of WGAP, ie the
distance
between the walls, in the formulae.
 This enables a conductiontype model to be used even with
nonparticipating media.
 Of course, an economical means of calculating WGAP is needed.
 This is provided by the LTLS equation (see below).
 IMMERSOL gives quantitatively correct predictions in geometrically
simple circumstances and plausible ones in complex ones.
 It is economical enough to be generalised for wavelengthdependent
radiation.
next, back or
contents
(b) WGAP, WDIS and LTLS
 The solved differential equation is:
div (grad LTLS) = 1
 The boundary conditions are:
LTLS = 0 , in all solids.
 The solution for plane channel flow is:
LTLS = WDIS * ( WGAP  WDIS ) / 2 , and
grad LTLS = WGAP / 2  WDIS
 These relations are supposed to prevail also in two and
threedimensional circumstances.
 That is all!
next, back or
contents
Notes:
 The LTLS equation is very simple, and therefore easy to solve.
 Its solution yields values of LTLS and grad LTLS at
all points in the field.
 WDIS and WGAP are then deduced from them.
 Their values are quantitatively correct predictions in geometrically
simple circumstances and plausible in complex ones.
 The method is especially useful, indeed the only practicable one,
when the space in question contains many solids of arbitrary shapes.
next, back or
contents
(c) LVEL: summary
next, back or
contents
Notes:
 The LVEL model is very simple, and therefore easy to implement.
 The predicted effective viscosities are quantitatively correct
in geometrically simple circumstances and plausible in complex
ones.
 The method is especially useful, indeed often the only practicable one,
when the space in question contains many solids of arbitrary shapes.
 LVEL handles the complete Reynoldsnumber range: laminar, transitional
and fully turbulent.
 LVEL can be easily extended so as to improve its accuracy in locations
far from walls.
Click
here for an SFT example involving LVEL in natural convection.
next, back or
contents
5. Discussion of current achievements and future prospects
(a) Preliminary conclusions
The following conclusions appear to be justified:
 Simultaneous simulation of solidstress, heat transfer and fluid flow
is indeed practicable and economical.
 As compared with the alternative, namely the use of distinct methods
for each phenomenon with iterative interaction between them, the
simultaneoussolution method is very attractive.
 It therefore seems possible that, when its attractiveness is fully
recognised, SFT (i.e. SolidFluidThermal) analysis may become as
popular as CFD.
next, back or
contents
(b) Why has this recognition not already become widespread ?
The author now proposes several answers to this question, as follows:
 First, it requires longheld convictions to be questioned and then
abandoned. Many people find this uncomfortable.
 Secondly, those who are in principle willing to try new methods prefer
to do so only after very many others have led the way.
"Never do anything first" commends itself to those who have
observed how seldom pioneers actually prosper.
next, back or
contents
 Thirdly, in the present case, The pioneers ( the author and
his colleagues) have been reprehensibly backward in documenting,
illustrating and generally promoting the SFT method.
True: it has been a part of the PHOENICS package for several years, but
with too many limitations to be easily usable; and finally
 The method described above has proved to have two
deficiencies, namely:
 the SIMPLE algorithm converges rather slowly when both the flow and
solidstress situations are complex and intertwined; and
 additional 'source terms' must be provided in order that
'bending' phenomena can be properly represented.
How to overcome these deficiencies will now be described.
next, back or
contents
The better algorithm; revival of an old trick
Mathematical aspects
Click here for full details
 Differentiate:
the Udisplacement equation with respect to y
and :
the Vdisplacement equation with respect to x.
 Subtract one from the other, causing the dilatation and
thermalexpansion terms to disappear.
 Hence obtain an easytosolve transport equation for the z direction component of
vorticity (the CFD word) or rotation (the stressanalyst's
word) k of the form:
[del**2]*k + K*C2 = 0
where:
k = dU/dy  dV/dx and
K = dFx/dy  dFy/dx
next, back or
contents
 Do the same for the V and W equations and for the
W
 and U equations,
thus obtaining equations for the other two vorticity
(i.e. rotation) components, i and j.
 Differentiate the definition of the displacement, D with respect
to x, y and z, and thence obtain 3 expressions (only 1 is shown)
for its gradients, of the following form:
dD/dx = [del**2]*U + dj/dz  dk/dy
 Finally, by algebraic manipulation, derive secondorder differential
equations for U, V and W of the form:
[del**2]* U + [dj/dz  dk/dy]*C4 +[d/dx]* [ Te*C3 ] + Fx*C2 = 0
which is simpler than the earlier equation for
[del**2]* U because [dj/dz  dk/dy] is known from the
solution for vorticity.
next, back or
contents
New solution procedure
 First solve the three (linear, Poissontype) equations for the vorticity
components.
Often analytical solutions will suffice.
 Solve the three (linear, Poissontype) equations for the displacements
(in solids) and velocities (in fluids), which are now linked only at the
solidfluid boundaries.
 Iterate as necessary, but many fewer times than for pressurelinked
equations (i.e. SIMPLE/SIMPLEST).
 The results of such a calculation can be seen
here, where the
bending of the beam is clearly visible.
Click here for
more details.
next, back or
contents
Present stage of development
 The new algorithm has been found to 'converge' much more rapidly than
the old.
 If appropriate expressions for 'vorticity sources'
(I, J and K) are introduced so as to express the
'bending moments' exerted by the applied forces, bending phenomena
are well simulated.
 The forces exerted by the fluid on the solid are easily taken into
account they are calculated simultaneously.
 There appear to be no difficulties of principle in the way of
extension to:
 timedependent phenomena;
 deformations which are large enough to affect the fluid flow;
 and perhaps even plastic deformations.
 SFT therefore appears to be ready to take off like
this or
this.
next, back or
contents
(d) Research opportunities
Singlealgorithm SFT presents numerous opportunities, both to numerical analysts and to
practicallyminded researchers, for example, to:
 extend the algorithm to curvilinear grids and to cartesian grids
cut by curved and/or moving solidfluid interfaces;
 extend to transient phenomena;
 allow deformations which are large enough to affect the flow;
 permit plastic as well as elastic deformation of the solid material;
and
 allow the fluid to be multiphase.
An offer
CHAM will be happy to offer freeofcharge PHOENICS licences to academic
researchers willing to collaborate in this important field.
END of LECTURE
References, back or
contents
6. References
 The differential equations governing displacements, stresses and
strains in elastic solids of nonuniform temperature can be found
in numerous textbooks, for example:
 CT Yang
Applied Elasticity
McGraw Hill, 1953
 BA Boley and JH Weiner
Theory of Thermal Stresses
John Wiley, 1960
 PP Benham, RJ Crawford and CG Armstrong:
Mechanics of Engineering Materials
Longmans, 2nd edition, 1996
It has not been common to choose the displacements as the
dependent variables in numericalsolution procedures. However, this
has been done by:
 JH Hattel and PN Hansen
A ControlVolumebased FiniteDifference Procedure for
solving the Equilibrium Equations in terms of
Displacements
Applied Mathematical Modelling, 1990
Their numerical procedure differ from that used here, which was that of
 SV Patankar and DB Spalding
"A Calculation Procedure for Heat, Mass and Momentum
Transfer in ThreeDimensional, Parabolic Flows"
Int J Heat Mass Transfer, vol 15, p 1787, 1972
 The first use of the present method for solving the soliddisplacements
and fluidvelocity equations simultaneously appears to have been
made by CHAM, late in 1990.
Reports describing the early work include:
 KM Bukhari, HQ Qin and DB Spalding
Progress Report (to RollsRoyce Ltd) on the Calculation of
Thermal Stresses in Bodies of Evolution
CHAM Ltd, November, 1990
 KM Bukhari, IS Hamill,HQ Qin and DB Spalding
StressAnalysis Simulations in PHOENICS.
CHAM Ltd, May, 1991
From that time onwards, the solidstress option was made available
as a (littleadvertised) option in successive issues of
PHOENICS,
 Openliterature and conference publications have been few, but
include:
 DB Spalding
Simulation of Fluid Flow, Heat Transfer and Solid Deformation
Simultaneously
NAFEMS 4, Brighton 1993
 D Aganofer, Liao GanLi and DB Spalding
The LVEL Turbulence Model for Conjugate Heat Transfer at
Low Reynolds Numbers
EEP6, ASME International Mechanical Engineering Congress and
Exposition, Atlanta, 1996
 DB Spalding
Simultaneous Fluidflow, Heattransfer and Solidstress
Computation in a Single Computer Code
Helsinki University 4th International Colloquium on Process
Simulation, Espoo, 1997
 DB Spalding
FluidStructure Interaction in the presence of Heat
Transfer and Chemical Reaction
ASME/JSME Joint Pressure Vessels and Piping Conference, San
Diego, 1998