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 4-8, 2002, Vancouver, British Columbia, Canada

updated March 2004

abstract or contents or update



This is not the only common fallacy that will be exposed. Specifically:

  1. It is often believed that thermal radiation between solids separated by non-participating gases must be computed by view-factor or ray-tracing techniques.



  2. It is often believed that convective heat transfer at low Reynolds number, which prevails for example in electronics equipment, requires an advanced 2-or-more-differential equation model for its prediction.



  3. 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.



To summarise, four misconceptions are to be dispelled, namely that:

  1. solid-stress analysis needs a separate program;

  2. solid-to-solid radiation needs view-factors or ray tracing;

  3. low-Re convection needs (e.g.) Lam-Bremhorst;

  4. wall-distance computation requires trigonometry.

next or back


  1. The problem
    1. Its essential nature
    2. Practical occurrence
    3. The conventional solution
    4. A better solution

  2. A multi-physics example
    1. Stresses resulting from radiation, conduction and convection
    2. Vector and contour plots
    3. How the stress calculations were performed

    next or back or contents

  3. The mathematics of the method
    1. Similarities between the equations for displacement and velocity
    2. Deduction of the associated stresses and strains
    3. The "SIMPLE" algorithm for the computation of displacements
    4. More details of the equations

  4. Details of the auxiliary models
    1. IMMERSOL, for radiation
    2. WGAP, WDIS and LTLS, for radiation and turbulence
    3. LVEL, for turbulence

  5. Discussion of achievements and future prospects

  6. References

next or back or contents

1. The problem

(a) Its essential nature

It is frequently required to simulate fluid-flow and heat-transfer 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:

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 multi-physics example

(a) Description (PHOENICS library case, s400):

The task is to calculate the stresses in radiation-heated 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:

  1. 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.

  2. The radiative heat transfer is represented by the conduction-type IMMERSOL model, which is:
    1. economical and
    2. 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

  3. Both LVEL and IMMERSOL make use of the distributions of:

    1. distance from the wall (WDIS) and

    2. distance between walls (WGAP), both of which are calculated by solving a scalar equation for the

    3. LTLS variable.

  4. The stresses within the metals result primarily from the differences in their thermal-expansion coefficients. namely:

next, back or contents

(b) Vector and contour plots

next, back or contents

(c) How the stresses were calculated

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

  1. The x-direction displacement, U, obeys the equation:


    next, back or contents

next, back or contents


  1. The two equations are now set one below the other, so that they can be easily compared:

  2. The equations can thus be seen to become identical if:

    next, back or contents

  3. The expressions for C1, C2 and C3 are:

    next, back or contents

  4. A solution procedure designed for computing velocities will therefore in fact compute the displacements if:

    1. the convection terms are set to zero within the solid region: and

    2. a suitable linear relation between:
      • D ( ie [d/dx]* U + ...)


      • p
      is introduced by inclusion of a pressure- and temperature-dependent "mass-source" 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.


ex = [d/dx]* U

ey = [d/dx]* V

ez = [d/dx]* W

next, back or contents

Then the corresponding:

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}


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 mass-conservation constraint.

Its essential features are:

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 mass-source 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:

  1. eliminating the convection terms (ie setting Re low); and
  2. making D linearly dependent on p and temperatureT.

The "staggered grid" used as the default in PHOENICS proves to be extremely convenient for solid-displacement 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

next, back or contents


next, back or contents


next, back or contents


next, back or contents

(c) LVEL: summary

next, back or contents


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:

  1. Simultaneous simulation of solid-stress, heat transfer and fluid flow is indeed practicable and economical.

  2. As compared with the alternative, namely the use of distinct methods for each phenomenon with iterative interaction between them, the simultaneous-solution method is very attractive.

  3. It therefore seems possible that, when its attractiveness is fully recognised, SFT (i.e. Solid-Fluid-Thermal) 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:
  1. First, it requires long-held convictions to be questioned and then abandoned. Many people find this uncomfortable.

  2. 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

  3. 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

  4. The method described above has proved to have two deficiencies, namely:
    1. the SIMPLE algorithm converges rather slowly when both the flow and solid-stress situations are complex and intertwined; and

    2. 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

  1. Differentiate:
    the U-displacement equation with respect to y and :
    the V-displacement equation with respect to x.

  2. Subtract one from the other, causing the dilatation and thermal-expansion terms to disappear.

  3. Hence obtain an easy-to-solve transport equation for the z- direction component of vorticity (the CFD word) or rotation (the stress-analyst's word) k of the form:

    [del**2]*k + K*C2 = 0

    k = dU/dy - dV/dx and K = dFx/dy - dFy/dx

    next, back or contents

  4. 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.

  5. 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

  6. Finally, by algebraic manipulation, derive second-order 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

  1. First solve the three (linear, Poisson-type) equations for the vorticity components.

    Often analytical solutions will suffice.

  2. Solve the three (linear, Poisson-type) equations for the displacements (in solids) and velocities (in fluids), which are now linked only at the solid-fluid boundaries.

  3. Iterate as necessary, but many fewer times than for pressure-linked equations (i.e. SIMPLE/SIMPLEST).

  4. 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

next, back or contents

(d) Research opportunities

Single-algorithm SFT presents numerous opportunities, both to numerical analysts and to practically-minded researchers, for example, to:
  1. extend the algorithm to curvilinear grids and to cartesian grids cut by curved and/or moving solid-fluid interfaces;
  2. extend to transient phenomena;
  3. allow deformations which are large enough to affect the flow;
  4. permit plastic as well as elastic deformation of the solid material; and
  5. allow the fluid to be multiphase.

An offer

CHAM will be happy to offer free-of-charge PHOENICS licences to academic researchers willing to collaborate in this important field.


References, back or contents

6. References

  1. The differential equations governing displacements, stresses and strains in elastic solids of non-uniform temperature can be found in numerous textbooks, for example:

    It has not been common to choose the displacements as the dependent variables in numerical-solution procedures. However, this has been done by:

    Their numerical procedure differ from that used here, which was that of

  2. The first use of the present method for solving the solid-displacements and fluid-velocity equations simultaneously appears to have been made by CHAM, late in 1990.

    Reports describing the early work include:

    From that time onwards, the solid-stress option was made available as a (little-advertised) option in successive issues of PHOENICS,

  3. Open-literature and conference publications have been few, but include: