Encyclopaedia Index

### Contents

1. The problem
2. A multi-physics example
3. The mathematics of the method
4. A simple example of flow-influenced stress
5. References

### 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, it is the stresses which are of major concern, while the fluid and heat flows are of only secondary interest.

(b) Practical occurrence

Engineering examples of fluid/heat/stress interactions include:

• air-cooled gas-turbine blades under transient conditions;
• "residual stresses" resulting from casting or welding;
• thermal stresses induced in liquid-sodium containers and immersed structures in fast-breeder nuclear reactors during emergency shut-down;
• manufacture of bricks, ceramics and other materials;
• flows induced by the contraction of "heat-shrinking" plastics;
• stresses in the pistons and cylinder blocks of diesel engines;
• the failure of steel-frame-supported buildings during fires.

(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 awkwardness.

(d) A better solution

This interfacing is rendered unnecessary by the solid-stress option of PHOENICS, which permits fluid flow, heat flow and solid deformation, and the interactions between them, all to be simulated 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).

### 2. A multi-physics example

(a) Description

In order to demonstrate what this method can achieve, the results from PHOENICS library case s400 will be shown. This concerns the calculation of fluid-flow, conjugate heat transfer and solid stress in the situation illustrated by the sketch below.

```

cooling | air
|
| V |/////// hot radiating wall ///////////|
|   ----------------------------------------
|                       duct              ----->  exit
|-------------                 -------------
|// steel ///|     cavity      |/// steel /|
|-------------------------------------------
|////////////// aluminium /////////////////|
|-------------------------------------------
```

The task is to discover what stresses arise in the metal blocks as a result of the non-uniformities of temperature.

Further details 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 IMMERSOL model, which is economical and fairly accurate for such situations.

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

4. The temperature distributions in the metal blocks are the consequence of radiative heating and convective cooling at their exposed surfaces, and of heat conduction within them.

5. The stresses within them result primarily from the differences in their thermal-expansion coefficients. namely:
• 2.35 e-5 for aluminium, and
• 0.37 e-5 for steel.

(b) Vector and contour plots

• Vectors

The velocity vectors displayed in Fig.1 reveal the nature of the air-flow pattern.

They are calculated at the same time as the displacement vectors shown additionally Fig.2; but, of course, the two sets of vectors have different scales, and indeed dimensions.

The solids are supposed to be confined by a stiff-walled box, but are allowed to slide relative to its walls. This is why the displacement vectors are vertical near the confining-box walls.

They are however not allowed to slide relative to each other; this is what causes the concentrations of stress at their contact surfaces.

• Stresses and strains

The stresses in the x- (horizontal) and y- (vertical) directions are displayed in Fig.3 and Fig.4 respectively.

They have been deduced from the strains shown in
Fig.5
for the x-direction , and in
Fig.6
for the y-direction,

The strains have been deduced from the displacements by differentiation. They are displayed in
Fig. 7
for the x-direction and
Fig. 8
for the y-direction.

However, their implications are less visually dramatic,

• Temperature fields

Fig. 9 displays contours of temperature in the air and the solid, of which the most noticeable feature is that the air is heated by contact with the 80-degree-Celsius top wall and by the metal blocks which have been receiving heat by radiation from above.

Temperature differences within the high-conductivity solids are too small to be discerned visually.

It is interesting to compare Fig. 9 with Fig. 10,
This displays the distribution within the air space of the "radiation temperature", T3, which is computed by IMMERSOL, and which 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.

The solid temperature influences the stresses and strains, of course, primarily through the agency of the temperature-dependent thermal- expansion distribution. However, its variations with position are too slight to be revealed by a contour diagram, as inspection of Fig. 11 will reveal.

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 post-processing.

The results are displayed in
Fig. 12 for the x-direction. and by
Fig. 13 for the y-direction.

The shapes of these contour diagrams, if studied and interpreted in physical terms, will be found to be plausible.

• 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 easily-understood meaning when the walls are near, and nearly parallel, is computed from the solution of the "LTLS" equation.

The distributions of these two quantities are shown by Fig. 14 for the former, and by
Fig. 15 for the latter.

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.

• Contours of auxiliary quantities used in the fluid-flow calculation

The flow field was calculated by means of the LVEL turbulence model, which makes use of the wall-distance field which is also derived from the LTLS distribution.

The contours of WDIS are displayed in Fig. 16.
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.

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

Since the laminar viscosity is of the order of 1.e-5 m**2/s, it is evident that turbulence raises the effective value, far from the walls, by an order of magnitude.

### (a) Similarities between the equations for displacement and velocity

The similarities referred to in section 1 (d) are here decribed for only one cartesian direction; but they prevail for all three directions.

1. The x-direction displacement, U, obeys the equation:
[del**2]* U + [d/dx]* [ D*C1 - Te*C3 ] + Fx*C2 = 0

where:
• Te = local temperature measured above that of the un-stressed solid in the zero-displacement condition, multiplied by the thermal-expansion coefficient;
• D = [d/dx]* U + [d/dy]* V + (d/dz]* W
which is called the "dilatation";
• Fx = external force per unit volume in x-direction;
• V and W = displacements in y and z directions;
• C1, C2 and C3 are functions of Young's modulus and Poisson's ratio.
2. When the viscosity is uniform and the Reynolds number is low, so that convection effects are negligible, the x-direction 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 x-direction,
• c1 = c2 = the reciprocal of the viscosity.

Notes:

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

• [del**2]* U + [d/dx]* [ D*C1 - Te*C3 ] + Fx*C2 = 0

• [del**2]* u - [d/dx]* [ p*c1 ] + fx*c2 = 0 ,
2. 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

3. 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)
4. 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 stress region:
• the linear relation between D ( ie [d/dx]* U + ...) p and Te is introduced by inclusion of a pressure- and temperature- dependent "mass-source" term.

### (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 .

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}
tauxy = {YM / (1 - PR**2)} * {0.5 * (1 - PR)*gamxy}

where:

• gamxy = [d/dy]*u - [d/dx]*v

### (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:

• All the velocity equations are solved first, with the current values of p. (li) The consequent errors in the mass-balance 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.

All that it is necessary to do, in order to modify this algorithm for solving for displacements simultaneously, is to treat the dilatation B 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.

### (a) Description

The following sketch illustrates a combined mechanical-and-thermal-stress situation, in which the fluid flow plays a part. It is the focus of attention in a series of case to be found in the Input-File Library of the solid- stress option,

```
------------------------------------------------------

load V V                          <---------   ^ y
|-----------|                             |
|///////////|                cold-air jet |
|// heated /|                             |
|// block //|                             |
|///////////|                <------z-----O
----------------------------------------------

```

The task is to compute displacements, and thence strains and stresses, in the heated block, subject to a variety of constraints as follows:

• In case s101.htm, there is no load, the block is heated to a uniform temperature, and its

### (b) The computational problem

• The independent variables are y & z. (the flow is steady; and the body is axi-symmetrical, so neither time nor circumferential coordinate x has an influence)
• The dependent variables are:
v & w (velocities or displacements)
P (pressure or dilatation)
T (temperature)
• Derived variables: Te (linear thermal expansion)
ey, ez (strains in y & z directions)
stry, strz (stresses in y & z directions)
• Boundary conditions: as indicated on above sketch.
• The differential equations to be solved are:-
• the Navier-Stokes equations for the fluid region;
• the displacement equations (above) for the solid region;
• the thermal-energy equation for both regions.

### (c) How PHOENICS solves the problem

PHOENICS solves the fluid-flow, displacement and thermal-energy equations simultaneously, being enabled to do so easily because:

1. it does not need to solve for both displacements and velocities at the same place;
2. the displacement equations are so similar to those for velocity that the same solution algorithm will do for both.

The user needs only to set STRA=T in the Q1 file, and to provide appropriate boundary conditions for displacement, etc; then PHOENICS will test, cell-by-cell, whether a solid or a fluid is present, and will solve whichever equation is appropriate there.

This is done whenever the solid-stress option is present.

Only linear stress-strain relations are allowed for; but the Young's Modulus and the thermal- expansion coefficient can be temperature-dependent.

Graphically-displayed results for this case are as follows:

Fig.18 shows simultaneously the displacement vectors and the velocity vectors; the second shows the temperature field which caused the thermal expansion.

The arrows represent velocity vectors in the upper region and displacement vectors in the block, and Fig.19 shows the computed temperature distribution which caused (most of) the displacements.

A library of input-files is supplied with the solid-stress option of the PHOENICS.

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

• CT Yang
Applied Elasticity
McGrawHill, 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 numerical-solution procedures. However, this has been done by:

• JH Hattel and PN Hansen
A Control-Volume-based Finite-Difference 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 Three-Dimensional, Parabolic Flows"
Int J Heat Mass Transfer, vol 15, p 1787, 1972

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:

• KM Bukhari, HQ Qin and DB Spalding
Progress Report (to Rolls-Royce Ltd) on the Calculation of Thermal Stresses in Bodies of Revolution
CHAM Ltd, November, 1990

• KM Bukhari, IS Hamill,HQ Qin and DB Spalding
Stress-Analysis Simulations in PHOENICS. CHAM Ltd, May, 1991

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:

• DB Spalding
Simulation of Fluid Flow, Heat Transfer and Solid Deformation Simultaneously
NAFEMS 4, Brighton 1993
• D Aganofer, Liao Gan-Li 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 Fluid-flow, Heat-transfer and Solid-stress Computation in a Single Computer Code
Helsinki University 4th International Colloquium on Process Simulation, Espoo, 1997
• DB Spalding
Fluid-Structure Interaction in the presence of Heat Transfer and Chemical Reaction
ASME/JSME Joint Pressure Vessels and Piping Conference, San Diego, 1998