NASA CONTRACTOR 

REPORT 

I " 

NASA CR-178996 

A CRITICAL EVALUATION OF VARIOUS METHODS FOR THE ANALYSIS 

OF FLOW-SOLID INTERACTION IN A NEST OF THIN CYLINDERS 

SUBJECTED TO CROSS FLOWS 

By Sang-Wook Kim 

Universities Spac Research Association 

Systems Dynamics 

Science and Engin 

Interim Report 

January 1987 

Prepared for 

NASA-Marshall Space Flight Center 

Marshall Space Flight Center, Alabama 35812 

(EASE-CH-1789S6) B CRIIICAZ EQALUATIOI CI E187-1L7El 

VABIOCS PETEOCS FCB IEE ALALFSlS CP 

€LCW-SCL;ID INSEBACTICL 'IN A EESl Cf THLN 

Etc~ort (Universities Siace EiezEarch 63/34) 43649 

CYLlNfrEBS SUBJEC'JED TC CBCSS €IC%$ Interim UnCldS 

Ill 
ACKNOWLEDGNENTS 

The author expresses his thanks to Dr. N. C. Costes for reviewing the manuscript and for many helpful discussions. 

This work was supported by NASA contract NAS8-35918, SSPJE. 

i 
TABLE OF CONTENTS 

Page 

I . INTRODUCTION ....................................................... 

I1 . CROSS FLOW INDUCED VIBRATIONS .................................. 

2.1 Vortex Induced Vibration ......................................... 2.2 Fluctuating Turbulent Wall Pressure 

2.3 Fluidelastic Instability ............................................ 2.4 Other Vibration Inducing Mechanisms ............................. 

I- .............................. 

2.4.1 Aeroacoustic Noise ......................................... 2.4.2 Jet Switching .............................................. 2.4.3 Oscillating Mean Flows and Swirling Flows ................. Semi- Analytical Solutions Including Combined Excitation 

Mechanisms ....................................................... 2.5 

111 . NUMERICAL ANALYSIS METHODS FOR FLOW-SOLID INTERACTIONS ... 

3.1 Laminar Flow- Solid Interactions ................................... 3.2 Turbulent Flow-Solid Interactions ................................. 3.3 Turbulence Models ................................................ Extended k- E: Turbulence Models and Algebraic 

Reynolds Stress Models ................................... 3.3.1 

3.3.2 

Structural Analysis Methods Including Large and/or 

Permanent Deformations ........................................... 

3.5.1 Control Volume Finite Difference Analyses of 

3.5.2 Finite Element Analysis of Steady Two Dimensional 

Laminar Flows through Tube Banks ........................ 3.5.3 Body Fitted Grid Finite Difference Analysis of 

Steady Flows through Tube Banks ......................... 

Large Eddy Simulation Method ............................. 

Numerical Solutions of Flows through a Nest of Cylinders ......... 

Steady Laminar Flows through Tube Banks ................ 

3.4 

3.5 

IV . CONCLUSIONS AND RECOMMENDATIONS .............................. 

REFERENCES ................................................................ 

1 

1 

2 

4 

5 

10 

10 

12 

13 

13 

17 

18 

21 

25 

26 

31 

36 

44 

44 

46 

52 

52 

57 

Y 

iii 
LIST OF ILLUSTRATIONS 

Figure Title Page 

1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

10. 

11. 

12. 

13. 

14. 

15. 

16. 

17. 

18. 

19. 

Parameters used for a nest of cylinders ............................. 3 

Idealized response of a cylinder in a nest of cylinders subjected 

to cross flow 3 

Mean and lower boundary curves for fluidelastic instability .......... 9 

........................................................ 

Stability diagram for triangular arrays of cylinders .................. 9 

Possible tube arrangements to reduce aeroacoustic noise ............. 11 

Coupling of jets in the wake region of a row of cylinders 

subjected to cross-flow .............................................. 12 

Space Shuttle Main Engine (SSME) power head assembly ............. 14 

Flow field inside of the SSME-Main Injector Assembly ................ 15 

Finite element analysis of a nuclear reactor safety problem.. ......... 22 

Geometry of compliant surfaces, flow and coordinate system for 

finit e elem en t transient response simulation .......................... 23 

Turbulent recirculating flow in a nuclear reactor outlet plenum ...... 32 

Comparison of measured vertical velocity with computed velocity ..... 33 

The energy cascade process ......................................... 33 

An elasto-viscoplasticity analysis of a thick cylinder subjected 

to internal pressure ................................................. 41 

Control volume finite difference analysis of flow through a 

porous media.. ...................................................... 45 

A finite element solution of one dimensional convectiondiffusion problem .................................................... 47 

A Petrov-Galerkin finite element solution of a two dimensional .i 

steady laminar flow through a tube bank ............................ 50 

J 

Flow over a backward-facing step ................................... 51 

INS3D computation of flow through cylinders.. ....................... 53 

iv 
LIST OF TABLES 

Table Title Page 

1. Values of Kfe and a for equation (2.4) ............................... 7 

2. Values of coefficients for equation (2.5) .............................. 8 

3. Lower bound on critical flow velocity ................................. 8 

. 

4. Comparison of different numerical methods for the numerical 

analysis of flow-solid interactions. .................................... 55 

lV 

, 
NOMENCLATURE 

time rate of change of 3 

a matrix relating strain vector to nodal displacements, B = *= LN 

body force vector (includes acceleration force) 

inertia coefficient for displaced fluid due to vibration of structure 

a coefficient used in LES subgrid-scale model 

elasto-viscoplastic material property matrix 

diameter of a cylinder 

Young’s modulus 

yield function 

load vector in a system of equations 

dominant frequency of turbulent buffeting spectrum 

natural frequency of a cylinder 

vortex shedding frequency 

filtering function for LES method 

area moment of inertia 

identity matrix 

finite element stiffness matrix 

an experimentally determined constant for fluidelastic instability equation 

tangential stiffness matrix for elasto-viscoplastic finite element equation 

artificial viscosity in Petrov - G alerkin finite element method 

differential operator relating strain to displacement 

length scale of energy containing eddies 

mass matrix 

added mass, m = m +C vpd /4 

mass per unit length of a cylinder 

matrix of interpolating polynomials 

2 

Oa 

vi 
Ne 

'e 

Pm 

PS 

Rt 

Red 

'b 

s1 

St 

svs 

V 

g 

V 

ct 

ct 

P 

a 

Y 

yP 

n 

A tn 

6 

" 6 .. 

13 

6U z 

E 

number of nodes in an element 

Peclet number 

mode shape function for sound pressure 

sound pressure 

Reynolds number, p Ud / p 

turbulent Reynolds number 

Strouhal number based on dominant frequency of turbulent buffeting 

spectrum 

cylinder spacing in flow direction in a tube bank 

cylinder spacing in normal to the flow direction in a tube bank 

Strouhal number based on vortex shedding frequency 

displacement vector in spatial coordinate directions 

free stream flow velocity 

flow velocity between cylinders in a nest of cylinders 

velocity vector in spatial coordinate direction 

a constant coefficient used in fluidelastic instability equation 

penalty parameter for penalty finite element method 

a constant coefficient used in fluidelastic instability equation 

a constant coefficient used in fluidelastic instability equation 

fluidity of plastic flows 

twice of the computational mesh size in LES method 

time increment, ~t" = t n+l - tn 

dimensionless logarithmic decrement of amplitude due to fluid and structural 

damping 

Kronecker delta (=1 for i = j, = 0 for i # j) 

variation of displacements 

turbulent dissipation function 

F - a column vector of strain (or total strain) 

c 

-e 

a column vector of elastic strain 

vii 
E 

- VP 

0 

1-I 

vt 

Tr 

P 

w 

- w 

aR 

aR1 

2 

a column vector of viscoplastic strain 

interpolation parameter 

molecular viscosity of a fluid 

eddy viscosity 

a constant, T = 3.141592654 

fluid density 

radian frequency 

vorticity vector 

boundary of domain 

part of boundary on which Dirichlet boundary condition is specified 

part of boundary on which Neumann boundary condition is specified 

finite element basis polynomial, k = 1, Ne 

viii 
CONTRACT REPORT 

A CRITICAL EVALUATION OF VARIOUS METHODS FOR THE ANALYSIS OF 

FLOW-SOLID INTERACTION IN A NEST OF THIN CYLINDERS 

SUBJECTED TO CROSS FLOWS 

I. INTRODUCTION 

Cross flows over a cylinder or a nest of cylinders can be found in chimney 

stacks, aircrafts and missiles at high angles of attack, marine tow cables, offshore 

oil rigs and platforms, heat exchanger tube banks, and in the Space Shuttle Main 

Engine. These flows induce the structures to vibrate and the vibrating structures 

alter the flow field so that the fluid and the structural motions are strongly coupled 

together. It is imperative to understand the flow induced vibrations and to have 

predictive ability to avoid failure of structures. 

uation of various methods of analysis for flow induced vibrations of a nest of cylinders 

in cross flows as well as a critical evaluation of existing numerical methods for flowsolid interactions. 

This report contains a critical evalTopics included in Section 2 are : description of discrete vibration inducing 

mechanisms for a nest of cylinders ; experimental correlation equations derived through 

physical arguments or analytical modelings ; recent developments in semi-numerical 

modeling of the flow-induced vibrations due to multiple number of vibration inducing 

mechanisms; and accumulated experimental data which are used as design guides for 

tube banks. 

Flow-solid interactions have been studied for the last two decades. Nevertheless, application of major numerical analysis methods such as the finite difference 

method and the finite element method to these problems is still in its early development 

stage, and there exists only a limited number of publications related to these topics. 

In-depth prediction methods to be developed in the future for flow-solid interactions 

of a nest of cylinders in a cross flow, a problem which is geometrically and materially 

nonlinear in nature, would require major improvements in numerical analysis methods, 

unsteady turbulence models, and structural analysis methods including large and /or 

permanent deformations. Recent developments on these topics are discussed in 

Section 3. 

In addition to reviewing the flow-solid interaction analysis methods for tube 

banks, most recently available experimental data as well as experimental correlation 

equations are included in this report in order to provide a stand-alone guideline for 

design of tube banks. Also some of the finite element methods for flows as well as 

for structures, which may be used advantageously in the future development of a 

numerical analysis method for flow-solid interaction of tube banks, are discussed in 

some detail. 

* 

11. CROSS FLOW INDUCED VIBRATIONS 

Vibration of a cylinder or a nest of cylinders in cross flows is induced by 

vortex shedding, fluctuating turbulent wall pressure (buffeting), fluidelastic coupling, 

jet switching, oscillating mean free stream flows, and aerodynamic noise generated in 
the flow field. 

act synergistically rather than separately. Consequently , detailed analytical modeling 

of the fully coupled flow-solid interaction is difficult. But depending on the arrangement of cylinders, flow velocities, and Reynolds numbers, a few of the vibration 

generating mechanisms usually dominate the rest so that rather simple analytical 

models can be developed. Most of the available analytical models are for these simple 

flow-solid interaction cases; yet these are also important tools in preliminary design 

analysis of tube banks. 

lation equations, obtained by curve-fitting of experimental data, are used as design 

criteria. 

experimental correlation equations , experimental data , and semi-numerical solutions 

are discussed below. 

In many practical cases, all of the vibration generating mechanisms 

For practical design of tube banks, the experimental correEach of the above vibration generating mechanisms , analytical models , 

2.1 Vortex Induced Vibration 

In the frontal area of a cylinder, boundary layers are formed on the top and 

bottom surfaces of the cylinder. 

toward the rear end, the linear momentum of fluid particles is dissipated in overcoming the adverse pressure gradient in the rear part of the cylinder. When the 

momentum is completely lost, the boundary layers separate and form two free shear 

layers which bound the wake. Since the fluid particles near the center line travel 

much slower than those near the free stream, these shear layers roll up into discrete 

vortices. 

number , Re is changed [ 11 . Discrete laminar vortices begin to appear at approximately Red = 40. 

boundary layer flows become turbulent for Re 

narrower and disorganized up to Red = 3.5 x 10 ; and organized turbulent vortex 

street is formed again for Red higher than 3.5 x 10 . 

only at low Reynolds numbers, vortex shedding is the major source of mechanical 

excitation. 

stream is turbulent (i.e. , beginning from turbulence intensities 2 to 4 percent) , the 

fluctuating turbulent wall pressure (buffeting) becomes a source of mechanical excitation which is of equal importance as the vortex shedding. 

velocity is increased or decreased so that the vortex shedding frequency approaches 

the natural frequency of the cylinder, resonant vibration with large amplitude is 

produced. 

As the fluid particles in the boundary layers travel 

The pattern of vortex shedding for a rigid cylinder changes as the Reynolds 

The vortex flow becomes turbulent at Red as low as 150; the 

d’ 

> 3 x lo5; the turbulent wake becomes ds 

6 

For flow- solid interaction occurring when a cylinder is subjected to cross flow , 

As soon as turbulence begins to appear in the flow field or if the free 

As the free stream flow 

For a single row of cylinders placed perpendicular to the flow direction, vortex 

shedding is the major source of mechanical excitation only for small values of V /f d, 

gn where V is the flow velocity between the cylinders and fn is the natural frequency 

of a cylinder. As soon as V /f d is increased, fluctuating turbulent wall pressure 

becomes an equally important source of excitation as the vortex shedding; and for 

V /f d > 75, jet switching dominates the rest of vibration inducing mechanisms [3,41. 

g 

gn I 

gn 

Cylinders in a closely spaced tube bank do not respond as single cylinders; 

rather, interaction with the flow field causes coupled motion in the group of cylinders. 

The vortex shedding is confined to a range of moderate to large cylinder spacing 

ratios in tube banks [21. For St/d < 1.5, vortex shedding would be suppressed due 
to proximity to the cylinders [ 11, where S is the cylinder spacing in normal to the 

flow direction as shown in Figure 1. 

cylinders (tube bank) to increasing flow velocity is shown in Figure 2 [ 51. It can be 

seen in Figure 2 that the resonance to fluid periodicity, i.e., resonance to fluctuating 

turbulent wall pressure or vortex shedding, and fluidelastic instability are the most 

important mechanisms that can cause failure of tube banks. 

t 

Idealized response of a cylinder in a nest of 

FLU I DE LAST IC 

I NSTAB I L ITY 

RESPONSE TO FLOW 

RONRESONANT 

"PER IODIC I TY" 1 1 

Figure 1. Parameters used for a nest of cylinders. 

- 

v 

Figure 2. Idealized response of a cylinder in a nest of cylinders 

subjected to cross flow. 

3 
Vibrations due to vortex shedding and buffeting are classified as forced vibrations since these two vibration inducing mechanisms do not require flow-solid interactions and can exist for flow around rigid bodies. 

equations, used by tube bank designers to avoid resonance to vortex shedding, are 

usually developed on the basis of the contribution of Chen [ 61. On the other hand, 

correlation equations developed on the basis of the theory formulated by Owen [7] 

provide design means to avoid resonance to fluctuating turbulent wall pressure 

(buffeting). 

for resonance represent the same physical phenomenon and the only difference depends 

on the physical interpretation of the mechanism that induces resonance. In fact, as 

shown in Figure 2, the vibration due to vortex shedding, if it exists, and buffeting 

coexist in the low free stream velocity region. 

The experimental correlation 

Later, Paidoussis [ 81 showed that the two seemingly different conditions 

The correlation equation derived by Chen [6], for resonance in vortex induced 

vibration, is given as: 

f d/V = Svs , ng 

where Svs = fvsd/V 

and the rest of the notations are the same as before. 

geometrical arrangement of cylinders can be found in Reference 6. 

vortex shedding can be avoided by rendering the left hand side of equation (2.1) to 

be different, by 40 to 50 percent, from the Strouhal number which will cause resonance 

for a particular arrangement of cylinders. 

is the Strouhal number, fvs is the vortex shedding frequency, 

g A table of Svs for various 

Resonance to 

2.2 Fluctuating Turbulent Wall Pressure 

Turbulent flows induce vibration of structures through the fluctuating turbulent 

I 

wall pressure acting on the surface of structures submerged in the flow field. 

A nest of cylinders subjected to cross flow exhibit a peak amplitude at a flow 

speed whose dominant frequency of fluctuating turbulent wall pressure spectrum coincides with the natural frequency of the submerged cylinders. 

mental correlation equations used for predicting resonance which is caused by buffeting are based on the expression developed by Owen [7]. 

Most of the experiThe correlation equation due to Owen [7] is based on the assumptions that: the 

conditions for fully turbulent flow have been developed through regions sufficiently 

deep inside the tube bank; the rates of production and dissipation of the turbulent 

kinetic energy are in equilibrium; the length scale of the energy containing eddies is 

of the same order of magnitude as the cylinder spacing in the flow direction; the 

source of vibration is associated with the randomly fluctuating pressure forces imposed 

on the cylinders by the turbulent eddies; the flow velocity is low enough so that the 

fluidelastic vibration would not occur, yet the Reynolds number is high enough so that 

the effect of molecular viscosity is negligible; the tube spacing normal to the flow 

direction, St/d, is appreciably greater than unity so that a cylinder preserves some 

of the features of a circular cylinder in isolation but not large enough for these 

features to dominate; and the rate of dissipation of turbulent energy and the work 

done by the mean pressure gradient across any row of cylinders are in balance as has 

I 

i 
been observed in experiments. The correlation equation due to Owen [ 71 states that 

buffeting resonance would occur if: 

where V 

is a Strouhal number, fb is the dominant frequency of the turbulent buffeting spectrum, and the rest of the notations are the same as before. 

f is given as: 

is defined to be the flow velocity at the minimum gap in a row, Sb = fbd/V 

g g 

According to Owen 171, 

b 

2 d [3.05 (1 - -) + 0.281 1 V 

fb = 

R (1 - d/St) St Y 

where R is the length scale of the energy containing eddies, and the rest of the notations are the same as before. It was mentioned that equations (2.1) and (2.2) represent the same phyiscal resonance conditions at low free stream velocities [ 51. Nevertheless, there exists significant quantitative difference between these two equations if 

they are plotted on the same graph. Therefore, equations (2.1) and (2.2) can provide only rough criteria for avoiding resonance to forced vibration in the design of 

a tube bank. 

A prediction method, which utilizes the random vibration theory [ 9,101 for the 

maximum amplitude in resonant forced vibration was developed by Chen and Wambsganss 

[ 11,121. 

of fluctuating turbulent wall pressure. 

distribution of turbulent wall pressure as function of time nor the spectra and frequency composition are known a priori. 

plate measurements have been used in References 11 and 12 to predict the maximum 

amplitude. The specific example problem in References 11 and 12 will not be discussed 

herein. 

in References 9 and 10. 

In this method, it is required to know the spectra and frequency composition 

For preliminary design analysis, neither the 

Therefore, a set of data obtained from flat 

General procedures for flow induced random vibration analysis can be found 

Recently, methods for predicting fluctuating turbulent wall pressures [ 13,141 

and more detailed experimental data [ 151 have begun to appear in the open literature. 

Use of such a prediction method in a finite element analysis of compliant plates has 

been reported in References 16 and 17. These will be discussed in Section 3.2. 

2.3 Fluidelastic Instability 

Tube banks subjected to increasing flow velocity begin to vibrate with large 

amplitude at a critical flow velocity with a definite frequency or a group of frequencies which are different from that of the natural vibration mode of the cylinders. 

Above the critical flow velocity, the amplitude of vibration increases as the square of 

the flow velocity, apparently without limit. In this type of vibration, called fluidelastic or aeroelastic vibrations, the flow field excites the structure to vibrate and 

in turn the vibrating structure disturbs the fluctuation inducing flow field. Accordingly, this is a fully coupled flow-solid interaction phenomenon and is classified as a 

self-excited vibration. 

5 
Fluidelastic instability is of utmost concern in tube bank designs because it can 

cause total failure in the structural system. On the other hand, resonance to forced 

vibrations, i.e., vibrations due to vortex shedding and/or buffeting is of less concern, because it can be avoided by keeping the flow velocities away from the expected 

resonance velocity by a factor of 40 to 50 percent. 

The analytical expression for the onset of fluidelastic instability for a row of 

cylinders placed perpendicular to the flow direction was first proposed by Connors 

[ 41 , and is given as: 

n 

where the gap velocity V is defined as V 

g g 

or a function of cylinder spacing to be determined from experiments; a is a constant 

to be derived from theoretical analysis, yet it can be modified to fit experimental data 

more precisely; p is the density of the fluid; m is the added mass defined as rn = mo + 

Canpd /4; mo is the mass per unit length of the cylinder; Ca is an experimentally 

determined inertia coefficient for displaced fluid due to structural vibration [ 11 ; IT = 

3.14159; 6 is the logarithmic decrement of amplitude due to fluid and structural 

damping which is dimensionless; and the rest of the notations are the same as before. 

Equation (2.4) was derived by assuming that the energy supplied to the structural 

vibration would be equal to the energy dissipated by the damping force to maintain 

equilibrium in fluidelastic vibration. 

= V[ ( St/d) /{( St/d)- 111 ; Kfe is a constant, 

2 

Following the line of Connors [ 41, Blevins [ 181 formulated an analytical model 

The underlying assumptions in his modeling 

the vortex shedding frequency of the cylinders is well above the natural frewhich is more rigorous mathematically. 

are: 

quency of the cylinders so that the cylinders are not excited by vortex shedding or 

buffeting; flow velocity is low enough so that jet switching is inactive; and the forces 

due to asymmetry of the flow pattern caused by the displacement of cylinders from 

equilibrium is the only source of vibration. Blevins [ 181 sub-divided the fluidelastic 

instability criterion into three cases depending on vibration frequencies and damping 

forces along and transverse to the flow direction. These are: a symmetric frequency 

and symmetric damping force case, an asymmetric frequency and asymmetric damping 

force case. The first case of Blevins [ 181 reproduced exactly the same result as that 

obtained by Connors [ 41. Later, Heinecke [ 191 showed that the same equation, equation (2.4), can be derived by simple dimensional anaIysis and energy consideration of 

a dampled, hormonically excited oscillator. 

Blevin's [ 181 and Heinecke's [ 191 models, it might be suspected that the Connors' 

criterion for fluidelastic instability is not conservative. 

of tube banks, which satisfied the Connors' fluidelastic instability criterion, can be 

found in Reference 8. 

From the underlying assumptions of 

In fact, a list of failure cases 

Numerous experimental correlation equations have been derived following the line . of Connors' work by incorporating more experimental data into equation (2.4). 

order to provide a stand-alone guideline for tube bank designs, different values of 

the coefficient Kfe and exponent ~1 in equation (2.4), proposed by many different 

groups of workers in the area, are summarized in Table 1 [ 201. 

In 

6 
TABLE 1. VALUES OF Kfe AND a [20] 

Investigators 

Zonnors [4] 

Blevins [181 

Y. N. Chen 1211 

Gross [221 

Gorman [23 

Savkar [24] 

Connors [251 

Pettigrew et al. [ 261 

Weaver and Grover [27] 

Chen and Jendrzejczyk 

Tanaka et al. [29] 

fe 

Re-0.25 

4/k 

3.3 

4.95 (St/d)2 

0.37 + 1.76 St/d 

3.3 

7.1 

2.49 to 6.03 

3.0 

0.5 

0.5 

0.5 

1.0 

0.5 

0.5 

0.5 

0.5 

0.21 

0.2 to 1.08 

0.75 

Remarks 

Tube row with St/d = 1.42 

Cx and C are fluidY 

elastic- stiffness force 

coefficients 

Re = Reynolds number, 

f3 = constant 

For square array, and k 

determined from fluid 

force 

Suggested design guideline 

For triangular arrays 

For square array 

1.41 5 St/d I 2.12 

Suggested design guideline 

Rotated triangular array 

St/d = 1.375 

arrays and mixed array 

in water flow 

For various rectangular 

For square array, 

St/d = 2.0 

Another type of empirical correlation equation has been developed by means of 

similitude analysis [ 271 . This equation is given as: 

where a, 8, and y are exponents to be determined from experiments; B is a constant, 

or might be a function of other nondimensional parameters; the gap velocity, V 

defined in the same way as in equation (2.4); and the rest of the notations are 

is 

g’ 
Investigators B c1 B Y 

Paidoussis [81 2.3 0.5 0.25 1 

Paidoussis 1301 5.8 0.4 0.4 1 

TABLE 3. LOWER BOUND ON CRITICAL FLOW VELOCITY 1201 

Remark (Based on 

published data) 

Includes all data 

Excludes some data 

Array 

Tube row 

Parameter Range for 6m 

0.05 < 6m < 0.3 

0.3 < 6m .C 4.0 

4.0 < 6m < 300 

Square 

Rotated square 

Triangular 

Rotated triangular 

Vg/fnd 

0.06 1.35 (St/d - 0.375)6, 

2.30 (St/d - 0.375)6O,*5 

6.00 (St/d - 0.375)6:*5 

e 

9 00 

450 

30° 

6 Oo 

0.03 < dm < 0.7 

0.7 < 6m < 300 

0.15 

2.35 6 $50 

2.106 m 

I 3.54 (St/d - 0.5)6k5 

0.1 

6.53 (St/d - 0.9)6k5 

3.58 (St/d - 0.9)Sm 

0.17 2.86 m 

2.86,!p 

* e is defined in Figure 1. 

2 ** = m6/pd . 

8 

1- 
0 

100 I I- 

- 

- 

- 

- 

9 

IS 

9 10- 

r 

- c - 

UNSTABLE 0 

W - 

- > - 

n - 

W 

0 

3 

W 

K 

- 

- STABLE 

- 

n 

0.001 0.01 0.1 1 10 100 1000 1 

MASS DAMPING, rn(Zn{)/pd* 

Figure 3. Mean and lower boundary curves for fluidelastic instability 1311. 

100 

L - 

0 

0 

J 

W 

> 

J 

3 

k 10 

K 

0 

1 

Figure 4. Stability diagram for triangular arrays of cylinders [ 331 . 

9 
2.4 Other Vibration Inducing Mechanisms 

Other vibration-inducing mechanisms are : ( 1) aeroacoustic noise ; (2) swirling 

flows; and (3) oscillating mean free-stream velocity. These vibration sources have 

received less attention than the others considered earlier in this section, and are 

more difficult to treat, both analytically as well as experimentally, especially for tube 

banks. Another vibration inducing source is jet switching which is a type of fluidelastic vibration and can be observed most prominently for a single row of cylinders 

subjected to cross flows. Jet switching is included in this section. 

2.4.1 Aeroacoustic Noise 

Aeroacoustic sound is generated by vortex shedding and turbulence. If the 

Furthermore , if these two frequencies approach 

vortex shedding frequency coincides with the acoustic frequency , then the two 

mechanisms reinforce each other. 

the natural frequency of cylinders, then large amplitude vibration which can Bead to 

failure of tube banks can be caused. 

A review of numerous analytical and finite difference numerical solution of 

acoustic wave propagation in both homogeneous and inhomogeneous medium can be 

found in Candel [ 341. 

ence 34 are a direct wave propagation equation and the Helmholtz equation. The 

direct wave propagation equation is given as: 

The acoustic wave propagation equations considered in Referwhere ps is the sound pressure, t is time, c is the speed of sound, c = c(z), x_ is 

the spatial coordinates , and v is the Laplacian operator. The Helmholtz equation can 

be obtained from equation (2.6) by setting 

iwt Ps = Pm e Y 

where p,(x) is the mode shape function, e is the natural base of logarithm, i = fi, 

and u is the radian frequency. 

optima1 design of acoustic chamber can be found in Bernhard [351. 

A finite element analysis of Helmholtz equation for 

A review of sound induced by vortex shedding from cylinders is given by 

Blevins 1361. 

equation given as: 

The governing equation considered in the review is the Lighthill’s 

T.. , 22 a2 

ap c v p= -- axi a~ 11 

a t2 j 

10 
where T.. = pv-v - T~~ - c pAij, T~~ is the fluid stress tensor, 6ij is the Kronecker 

delta, and the rest of the notations are the same as before. 

2 

11 1j 

Recently, Quinn and Howe [37] presented an analytical solution for acoustic 

wave propagation through a tube bank. 

flat plates at zero angle of attack. 

in their equation by inserting aeroacoustic dipole sources in the forcing function term 

of the classical velocity potential equation for acoustic wave equation. 

decreasing frequency in such a way that onset of resonance is delayed. 

They modeled the cylinders as rigid finite 

The effect of vortex shedding was incorporated 

The result of 

. this analysis showed qualitatively that both attenuation and sound speed increase with 

Blevins [ 361 concluded that the degree to which self-sound synchronizes vortex 

shedding from various tubes is not known and that the role of acoustic damping had 

not been well understood until now. 

!Lore recently, an experimental study to eliminate aerodynamic noise in a tube 

bank was reported by Zdravkovich and Nuttall [ 381. 

with a tube bank of three rows of cylinders, with seven cylinders per each row. 

The array of cylinders is shown in Figure 5. 

resonance could be eliminated either by spacing unequally two successive rows of 

cylinders or by removing certain cylinders from the tube bank. Removal of a few 

cylinders from the first row was found to be the most effective way to eliminate 

acoustic resonance. For the second row, removal of a few cylinders was effective 

only at the pressure nodes; while removal in the third row was not effective at all. 

Experiments were performed 

It was shown that the acoustic 

Other approaches to reduce aeroacoustic noise was also found by Baird 1391 and 

Cohan and Dean [ 401. 

and higher modes of standing waves, but this technique may not be used advantageously for densely populated cylinders because of high pressure drop in the tube 

bank [40]. 

They recommended to use baffles in order to avoid the first 

+ + + 

++++ ++++ ++++ ++++ ++++ ++++ ++++ 

++++ ++++ ++++ ++++ ++++ ++++ ++++ I ++++ 2++++3++++4++++ 5++++ 6++++ ?++++ 

'+ 2 3+ 4+ 5+ 6+ 7+ 

I 

. \ 

PRESSURE NODE'- -- 

Figure 5. Possible tube arrangements to reduce 

aeroacoustic noise [ 371 . 

11 
2.4.2 Jet Switching 

In this phenomenon, fluid forces excite the cylinders to vibrate and the displacement of cylinders cause the jet pairs emerging from closely spaced cylinders to 

switch. Jet switching is a self-excited vibration which requires fully coupled flowsolid interaction. 

in Figure 6. 

Coupling of jets in the wake region of a row of cylinders is shown 

I f- / 2z- - 

Figure 6. Coupling of jets in the wake region of a row of 

cylinders subjected to cross-flow [ 11 . 

It was claimed in Reference 3 that jet switching may not occur for V /f d < 75, rSn 

and it was reported in Reference 4 that jet switching was not observed for V /f d < 

30. gn 

It was observed experimentally that the flow slightly inside the tube bank is 

suppressed so that jet switching is confined to a few cylinders at the flow inlet 

region. 

force as well as to jet switching. 

to damage due to fatigue and inelastic deformation due to large displacements. 

form analytical solutions cannot take this fact into account. Also numerical analysis 

for this type of problem has not been attempted as yet. 

correlation equations include this phenomenon implicitly , for experiments are concerned 

with failure of tube banks in an integrated sense. 

The cylinders at the flow inlet region are subjected to the greatest fluid 

Consequently , these cylinders are most susceptible 

Closed 

However, experimental 
2.4.3 Oscillating Mean Flows and Swirling Flows 

Oscillating mean flows and swirling flows excite a nest of cylinders submerged 

in the flow field to vibrate. 

oscillating mean flows is mostly based on the approach proposed by Morison [41], in 

which a cylinder is represented by an Euler beam and the fluid force acting on the 

cylinder is modeled as forcing function in the Euler beam equation. 

recent analyses have concentrated on modeling the fluid force due to oscillating mean 

flow, but flow-solid interaction has seldom been considered. 

Analytical analysis of flow induced vibrations due to 

Most of the 

It is natural to expect the swirling flows would induce vibration of cylinders submerged in the flow field. More importantly, the swirling component of the 

flow can alter the entire flow field so that the vibration mode of submerged cylinders 

may be different from those without swirl velocity. For example, it is well known 

that the flow field of jets with swirling component is quite different from those without 

swirling component [ 421. Swirling flow through a nest of cylinders can be found in a 

swirl crossflow stirling engine heater [43] or in the Space Shuttle Main Engine (SSME) 

Main Injector Assembly (MIA) (Fig. 7). A finite difference numerical analysis of the 

flow field inside MIA showed that there exists the swirling flow [44] (Fig. 8). As a 

remark, the vibration of cylinders has not been considered in Reference 43 or 44. 

Only the heat transfer rate was considered experimentally in Reference 43. Neither 

the details of flow field inside tube banks nor the effect of swirl flow on vibration 

of cylinders are known yet. 

2.5 Semi- Analytical Solutions Including Combined Excitation Mechanisms 

The vibration of cylinders in a tube bank is different from a collection of 

vibrating cylinders in isolation. 

coupled vibration of cylinders in tube banks were proposed by Blevins [ 181, Chen 

[ 451 , and Price and Paidoussis [ 461 . In all of these models, cylinders were represented by Euler Beams and the interaction between vibrating cylinders was included 

in the governing equations by specifying the forcing functions to linearly depend on 

the displacements and velocities of cylinders. 

that the amplitude of vibration is small enough so that the flow field is not disturbed 

due to vibration of cylinders. These models successfully describe the coupled vibration mode between cylinders, but can not be used with confidence for problems with 

large amplitude vibrations which lead to mechanical wear, fretting, corrosion, and 

fatigue failure of tube banks. 

More complicated mathematical models to account for 

It was also assumed in these models 

Blevins [18] considered "whirling" of a row of cylinders subjected to flows with 

low reduced velocities so that whirling vibration of cylinders is mainly due to lift 

force (caused by the asymmetry of flow pattern due to displacement of cylinders) and 

drag force. 

The flow velocity for the onset of whirling instability could be obtained from the 

analysis. 

Other vibration inducing mechanisms were excluded from his analysis. 

Chen [45] considered vibration of a row of cylinders subjected to a multiple 

number of excitation forces such as vortex shedding, fluidelastic coupling, drag and 

lift forces, and fluid inertia coupling. Apparently his model can be extended to 

vibration of tube banks with multiple rows once the fluid force coefficients for all the 

cylinders in the tube bank are provided by experiments or by other means, for 

example, computational fluid dynamics. 

discussed below, for its generality and completeness. 

The mathematical model due to Chen 1451 is 

13 
14 

k 

0 

a 

; 

h 

W 

a, 

E: 

a, 

0 

cd 

a 

cn 
F: FUEL SIDE, OX: OXIDIZER SIDE 

F 

0 = 160a 0 = 200 

Figure 8. Flow field inside the SSME Nain Injector Assembly [ 441 . (a) Horizontal velocity profile on the middle plane. (b) Vertical 

velocity profile at 30 deg plane from center line. 

15 
The equation of motion for a single cylinder, which has the same flexural rigidities in both flow direction and transverse to the flow direction, can be written as [45] : 

2 a U. a U-i -1 

at2 -' 

4 a U. 

a& 

+ mi - = f. Y Ei Ii - 4 + 'i at 

where the subscript i denotes the ith cylinder in the tube bank, R is measured along 

the cylinder axis, E is the modulus of elasticity, I is the area moment of inertia, g 

(ul,u2), u1 is the displacement in the flow direction (hereafter denoted as x-coorT- - 

dinate direction), u2 is the displacement in the transverse direction (hereafter denoted 

as y-coordinate direction), fT = (f ,f ), f 

difficulty in flow- solid interaction analysis lies in determining the forcing functions 

which depend on the displacement of cylinders. Chen 1451 decomposed the forcing 

function as: 

and f are forcing functions. The major XY x Y 

f. = f. C +_fig -I- f.e + f.h + f.O 

-1 -1 -1 -l Y 

where 

r) 

a &U -4 f.C = - N;' a.. - 

'1 at2 -1 

R sin (ust + $i 

[l + ci sin (2 wst + qi )] d R1 

-1 

C. 1 

1 2 f.e = 2 p V d (pij u.. + v.. u ) 

-1 11 11 2j 

f.h=-p Vd(s.. - G..+t..G ) 

-1 13 13 11 2j 

(2. lo) 

(2.11) 

(2.12) 

(2.13) 

(2.14) 

In equations (2.10) through (2.14), $: is the inertia force of fluid displaced by 

vibration of cylinders, M' is mass of displaced fluid, a.. is added mass coefficients for 

ith cylinder due to displacement of jth cylinder; fi is the column vector of drag force 

in x-direction and lift force in y-direction, c ' is the lift coefficient, c: is the drag 

coefficient, is a small number to be determined experimentally, us, us = 2~rS/d, is 

the vortex shedding frequency, $.' and +.d are phase delay angles for lift and drag 

11 R 

i 

1 1 

16 
is the fluid- ’ ij forces, respectively; L: is a column vector of fluidelastic forces, 

elastic force coefficients; ji 

hydrodynamic damping coefficients ; and i: is a column vector of forcing functions 

due to the rest of vibration inducing mechanisms. Equations (2.9) and (2.10) are 

quite general so that response of cylinders due to a specific vibration mechanism or 

any set of combined mechanisms can be analyzed by simplifying the load term, equation (2.9). But the quality of analysis will depend on the physical coefficients used 

in equations (2.9) and (2.10). 

analysis method, for natural frequencies and mode shapes of coupled vibrations, 

critical flow velocities for different excitation forces , and response to vortex shedding, 

lift, and drag forces for a row of cylinders have been presented in Reference 45. 

h is the hydrodynamic damping forces, s.. and t.. are 

11 11 

A set of analytical solutions, using the Eigen value 

Prices and Paidoussis [46] considered vibration of two rows of cylinders in cross 

flow. 

tests [ 47,481. 

the same as that of Chen [ 451. But they elaborated on modeling the effect of fluid 

retardness and time delay. 

cylinder in a tube bank of rigid cylinders yielded reasonably good approximation of 

critical flow velocity for tube banks with all flexible cylinders for m6 /pd2 less than 300, 

but not for m6 /pd greater than 300. Also, based on their analytical analysis, they 

recommended to separate 6 from m/pd in the experimental correlation equations, as 

in equation (2.5). 

Fluid force ciefficients used in the analysis were obtained from wind tunnel 

The form of governing equation of motion used in their analysis is 

The results showed that: analysis of a single flexible 

2 

2 

111. NUMERICAL ANALYSIS METHODS FOR 

FLOW- SOLID INTERACTIONS 

Geometrical and material nonlinearity inherent to flow- solid interaction phenomena 

precludes the possibility of obtaining closed form solutions to these problems except 

for simplified cases. 

numerically by solving unsteady turbulent flow equations , i.e. , unsteady NavierStokes equations with appropriate turbulence models. 

tribution and skin friction on the surface of the structure, one could then solve 

structural dynamics equations for vibration (and stability analysis) iteratively. 

vibration inducing mechanisms such as vortex shedding, fluctuating turbulent wall 

pressure, and jet switching may be obtained as part of numerical solution. 

It would be desirable to predict the flow-solid interactions 

By obtaining the pressure disThe 

Significant progress in numerical modeling of flow-solid interactions, turbulence 

models, and nonlinear structural analysis methods has been made recently. 

literature covers : 

tainers ; fully coupled laminar flow- solid interactions; and nonlinear stability, vibration, 

and elasto-viscoplastic analysis of structures subjected to known applied loads. For 

turbulent flow-solid interaction cases, the applied loads due to fluid flow have been 

obtained by solving the turbulent flow equations separately or by incorporating 

analytical or experimental correlation equations for fluctuating wall pressures induced 

by turbulence. Fully coupled numerical analysis of flow-solid interaction in cross 

flows, either laminar or turbulent, for a nest of cylinders has not been made yet. 

Existing 

unsteady turbulent flows around rigid bodies and /or rigid conSome of the numerical analysis methods for laminar flow-solid interactions, turbulent flow- solid interactions, mathematical models of turbulence, and structural 

17 
analysis methods which may be relevant to developing numerical methods for a class 

of flow-solid interaction of a nest of cylinders subjected to cross flows are discussed 

below. 

3.1 Laminar Flow-Solid Interactions 

Most of the recent investigations using numerical analysis for flow- solid interactions have been intended for safety analysis of nuclear reactors exposed to transient 

fluid loading. 

surrounding or submerged structures are of the same order of magnitude, and the 

time period of interest is very short (a few milliseconds for some cases). The effect 

of turbulent viscosity has not been included in these models. Researchers in this 

area developed a special mathematical model and numerical analysis methods which are 

suitable for this special class of flow-solid interaction problems. 

tions are described in Lagrangian-Eulerian coordinate system which was originally 

proposed by Noh [49] and further developed by Hirt [50] to be used for ail tlow 

velocities. The advantages of the Lagrangian-Eulerian description are : excessive 

distortion of computational grids can be avoided; and grid points can be displaced 

independently of the fluid motion so that implementation of fluid-structure interface 

condition can be simplified. 

In these problems , displacements excursed by fluid particles and 

The governing equaLet v.(j=1,2,3) be the velocity vector of a fluid particle, v.~ be the grid 3 I 

velocity vector, and v be the difference velocity vector defined as: 

j 

D G v. =v -v 

I jj (3.1) 

where j denotes the spatial coordinate directions. 

the fluid flow in the Lagrangian-Eulerian coordinate system are obtained to be: 

Then the governing equations for 

6.. - pbi = 0 a a V. 

+ vjD &) - - ax 11 

j j 

a vi 

j 

'ij =o j ax 

(3.3) 

(3.4) 

where equations ( 3.2) through ( 3.5) represent conservation of mass, linear momentum, 

energy, and constitutive equations respectively, p is the density, 0.. is the fluid 

13 

ia 
stress tensor, b. is the component of the body force along the ith coordinate direction, e is the internal energy, p is the molecular viscosity, and p is the pressure. 

Finite element representation of equations (3.2) through (3.5) are all based on 

the same method of weighted residuals, but detailed solution techniques used by 

different groups of workers, i.e., Donea et al. [51], Belytschko et al. [ 52,531, and 

Liu [ 541, are somewhat different. 

In the finite element method, each variable is interpolated using shape functions 

and nodal values of these variables. These are: 

where the shape functions $k(x), om(x), and x (x) are functions of spatial coordenotes number of these variables for each element. 

polating polynomials and inter-element continuity requirements are different depending 

on the differential equations to be solved and the numerical method employed. The 

finite element system of equations are obtained by multiplying the differential equations with appropriate test functions, substituting the continuous variables with 

interpolated expressions, integrating over an element, and assembling the resulting 

element system of equations into a global system of equations. 

S 

dinates; nodal values of these variables are functions of time; and N k9 Nm9 and Ns 

In general, the order of interDonea et al. [ 511 used linear shape functions for the velocities, and constant 

Consequently, p 

elements for density and internal energy. In this method, the pressure, p, was 

determined using an explicit relationship of the form p = p(p ,e). 

is constant in each element and the shape function x was not used in his case. 

a special form of interpolation as given in equation (3.9) was used for internal energy 

in [51]. The finite element system of equations due to Donea et al. [51] is as 

follows : 

Also 

Dp = fP - (3.11) 

(3.12) V Mv. = f. 

Ll -1 

19 
in which each of the entries in the above matrices are given by: 

(3.13) 

(3.14) 

(3.15) 

(3.16) 

(3.17) 

(3.19) 

In equations (3.10) through (3.13), D and E are diagonal matrices due to the shape 

functions used and the mass matrix M was made to be a diagonal matrix using a 

lumping method [ 551 . 

a R 

tion, and the rest of the notations are the same as before. The dissipation term due 

to molecular viscosity , was neglected from the conservation of 

linear momentum equation equation ( 3.3). The resulting discrete system of equations 

for this inviscid flow field was integrated in time direction using a predictor-corrector 

method. 

In equations (3.14) through (3.19) fie denotes an element , 

denotes the boundary of the same element , Ti denotes the applied surface trac- e 

( a v. / a x. + a v. / a x.) 13 11 

The structural domain was discretized into a number of conical shell elements 

[ 571 and a time dependent elasto-plasticity finite element method developed by 

Belytschko and Hsieh [57] was used for structural analysis in Donea et al. [511. The 

structural analysis part will be discussed separately in Section 3.4. 

20 
One of the example problems considered in Reference 51 was a thin cylinder 

vessel with hemispherical bottom almost completely filled with water, subjected to 

blast loading from inside. 

few milliseconds. 

Figure 9 together with moving finite element meshes. 

in Reference 51, a slipping flow boundary condition was used for tangential velocity; 

and the normal velocity was prescribed to be the same as the normal component of 

the moving boundary. This again is obtained as a part of numerical solution at each 

time-level. 

The system of equations were integrated in time up to a 

The physical problem and the computational results are shown in 

Since no viscosity was included 

The finite element method for flow analysis used by Belytschko et al. [51,53] 

is basically the same as the one used in Reference 51, but he used an averaging 

method for each of the variables in an element so that the element mass matrix would 

be a diagonal matrix. He also used an hourglass control algorithm to avoid zeroenergy mode resulting from the averaging method. In this case, the load vector fV 

contains an extra term contributed by the hourglass control algorithm in addition to 

the usual load vectors due to convection, molecular diffusion, surface tensions, and 

body forces. 

explicit time integration method, up to a few milliseconds in a finite element modeling 

of hypothetical nuclear reactor core disruptive accident. Their computational results 

compared favorably with the experimental displacement data. In the experiment 

permanent plastic deformation and severe buckling of steel columns were reported. 

In order to model the buckling of steel columns in the numerical analysis, a slight 

imperfection had been included in the geometry data of the steel column. The 

structural analysis method used can be found in Reference 57. 

The system of equations was integrated in time-direction, using an 

Liu solved the same set of flow equations using the penalty method to exclude 

the conservation of mass equation and pressure variable from the system of equations. 

The penalty method has been used very frequently for incompressible laminar flow 

problems [ 581 and the Stoke's problem [59]. In his method, the stress term, equation (3.5), was considered as an extra variable and was interpolated separately, 

which is called the mixed finite element method. 

model in the structural analysis part to cover strain rate sensitive solid materials. 

The elasto-viscoplasticity model is covered separately in Section 3.4. 

He included an elasto-viscoplasticity 

As can be noticed, there exists significant difference between the nuclear reactor 

In the latter problem, it is very important to include the turbulent viscosity 

Even at the onset of failure due to 

safety analysis problems and the flow-solid interaction of a nest of cylinders in cross 

flows. 

into the flow equations. 

vibrate with respect to its equilibrium position. 

large deformation, the distance excursed by any cylinder is small compared to the 

distance traveled by a fluid particle. 

Eulerian description may not be necessary for flow-solid interaction analysis of a nest 

of cylinders in cross flows. 

be used advantageously to obtain convergent solutions for convection dominated flows 

as has been pointed out by Zienkiewicz [ 581 . 

Any cylinder of a tube bank subjected to cross flows will 

For these reasons, direct use of the LangrangianBut the concept of Lagrangian-Eulerian description can 

3.2 Turbulent Flow - Solid Interactions 

Coupled and uncoupled turbulent flow-compliant material interactions were studied 

numerically by Buckingharn et al. [ 16,171 in an effort to reduce drag and aeroacoustic 

noise of submersible hulls. Uncoupled analysis was used to eliminate a number of 

21 
WATER LEVEL 

E TAIR GAP 

(h = 2.5 cm) 

.THIN VESSEL 

VI \ 1 

t=lms 

. . . - . . . . . 

t=2ms t=3ms 

(C) 

t=4ms 

8 

EMENTS 

1 

% 

L2 

9) 

0 

-2 

0 120 240 360 480 530 

HEIGHT FROM TOP (mm) 

(d) 

Figure 9. Finite element analysis of a nuclear reactor safety problem C511. 

domain. (c) Deformed configurations at t = 1, 2, 3, 

(a) Description of the physical problem. (b) Discretized computational 

and 4 ms. (d) Computational result. 

22 
candidate compliant materials and internal structures (Fig. 10) which do not exhibit 

the desired vibration response from the more elaborate and time consuming fully 

coupled flow-solid interaction analysis. Different materials considered include polyvinylchloride, rubbers, and many others. For all the cases of coupled and uncoupled 

turbulent flow-compliant material interactions considered in References 16 and 17, the 

structural dynamics equation was solved by the finite element codes due to Hallquist 

[60,61]. 

model the compliant materials. 

method for structural dynamics are discussed in Section 3.4. 

An elastic or a linear viscoelastic constitutive equation [62] were used to 

These constitutive equations and the finite element 

THIN, RELATIVELY STIFF 

NEOPRENEOVERLAYER 

CTANGU LAR 

RIGID SUBSTRAT 

Figure 10. Geometry of compliant surfaces, flow and coordinate system for 

finite element transient response simulation [ 16 , 171 . 

In the uncoupled analysis , the fluctuating turbulent wall pressure was generated 

using the Ash code [ 13,141. 

sure beneath fully turbulent boundary layer flows. 

pressure at a spatial and temporal location is obtained by summing up all the pressure 

events affecting the location. Information required to generate each of the pressure 

events include: 

rms 

(or turbulent wall shear stress and its relationship to prmS) ; probability distribution 

of time between pressure events : probability distribution of frequencies for the 

pressure events; convection speed of the pressure events; decay ratio of amplitude 

of pressure events; wave form of pressure carriers; and etc. In the work of Ash 

[ 13,141, all this information was obtained from several different sets of experimental 

data for turbulent boundary layer flows; and the statistical results of the generated 

fluctuating turbulent wall pressure compared favorably with experimental data. 

Universality of these experimental correlation equations for different flow situations 

need to be examined further. 

pressure, much statistical information must be known priori for the flow field, therefore his method is of limited use for general flow situations [ 131. 

Ash code can generate fluctuating turbulent wall presIn the method, the instantaneous 

root mean square value of fluctuating turbulent wall pressure, p 

In order to simulate the fluctuating turbulent wall 

In the fully coupled case [ 16,171, both the two- and three-dimensional NavierStokes equations were solved by a pseudo-spectral method [ 63,641 developed by 

Orszag and Kells [ 63,641 . The pseudo-spectral method [ 63,641 was originally 

developed to solve the Navier-Stokes equations for finite amplitude disturbances in 

order to predict transition criterion from laminar to turbulent flow of plane Poiseuille 

and plane Couette flows. The flow field, considered in References 16 and 17, is 

governed by the incompressible Navier-Stokes equations : 

- v.y=o (3.20) 

23 
av 

at 

- 

- + (v, . v) v = -1P + vv2 ,v z- (3.21) 

OiXlX , osysy 9 and 0szr;Z (3.22) 

I where equation (3.20) is the conservation of mass equation; equation (3.21) is the 

conservation of linear momentum equation; equation (3.22) is the domain of the prob1em;y - (vl,v2, v ) is a column vector of velocities in the three spatial coordinate Tdirections; 1 is the gradient operator; and X, Y, and Z denote extension of the 

domain. The domain in the z-direction was set at a distance of several boundary 

layer thicknesses away from the wall so that the free stream boundary condition can be used at the outer edge of the domain. 

compliant material interface were set to move together in the normal direction, i.e. ) 

the flow velocity in the normal to the surface of compliant material direction is 

dictated by the normal velocity of the compliant material which is obtained from the 

structural dynamics solution but the flow velocities in both of the tangential to the 

surface directions were set equal to zero. 

cyclic boundary condition such that v(O,y,z) = v(X,y,z) was used. In Orszag et al. 

[63,64], the Navier-Stokes equation was first rewritten in a rotational form given as: 

The fluid particles at the flowIn the free stream flow direction, a 

( 3.23) 

where g(x,t), dx,t) = yXy(x,t), is the curl of the velocity vector, and the rest of 

the notations aFe-the same as before. Equation (3.23) was solved by a pseudospectral method. 

in the large eddy simulation method which will be discussed in the next section. 

Therefore, a brief explanation and a discussion on the advantages and disadvantages 

of the method is given below. 

The pseudo-spectral method has also been used quite frequently 

The pseudo-spectra method they used belongs to the classical collocation 

method. 

are interpolated usir,g basis polynomials and constant coefficients assigned to each 

of the polynomials. 

In the collocation method, solution to the governing differential equations 

In their case, the velocity field was represented by: 

where m, n, and p are iztegers so that a set of (m,n,p) designate each collocation 

point; (ZM), (2N), and (P+1) are the number of collocation points in each of the 

spatial coordinate directions respectively ; a( m ,n ,p , t) is function of time only, and 

T (z) , T (z) = ~os(pcos-~z), is the Chebyshev polynomial of degree p. 

equation (3.24) into equation (3.2.4) and collocating at each of the collocation points 

Substituting P P 

I 
yields a system of nonlinear ordinary differential equations, which can be solved by 

any nonlinear ordinary differential equation solution method. 

computation of the collocation method requires enormous computational time as the 

number of collocation points is increased. 

method due to Orszag et al [ 63,641 are: 

utilized by performing most of the computations in the Fourier transformed space and 

then mapped back onto the physical space using the fast Fourier Transformation 

method [65] so that the required number of computations are reduced by orders of 

magnitude ; and the Adams-Bashforth-Crank-Nicholson time itnegration method devised 

to be used in the pseudo-spectral method is one of the best time integrators for 

ordinary differential equations with regard to truncation errors and numerical stability 

requirement. 

spectral method can be found in a number of literatures [ 66,67,68] ; therefore, it is 

not discussed herein. 

But a straight forward 

The special features of the pseudo-spectral 

orthogonality of the basis polynomials is best 

Detailed explanation on the computational procedure of the pseudoIn order to use the pseudo-spectral method to solve any ordinary or partial 

differential equations, basis polynomials which span entire domain and satisfy the 

prescribed boundary conditions must be available a priori. Due to this difficulty, the 

pseudo-spectral method has been used mainly for channel flows or axisymmetric 

annular flows where cyclic boundary conditions can be specified for the upstream and 

the downstream locations. The validity of cyclic boundary conditions has been 

questioned in References 67, 68, and 80. The boundary conditions are believed to 

be approximately applicable for fully developed channel or pipe flows [ 63,64,92] . Accordingly, the method may be suitable for problems with extremely limiting cases of 

boundary conditions defined on simple regular domains. 

to flow-solid interaction analysis of a nest of cylinders subjected to cross flows would 

be extremely difficult, if not impossible, due to the complicated domains and boundary 

conditions. 

Application of the method 

3.3 Turbulence Models 

In numerical analysis of a nest of cylinders subjected to cross flows, 

pressure straining and redistribution of Reynolds stresses near the wall region may 

become important due to a number of cylinders in the cross flow. 

interactions, fluid force is transmitted to the structure through wall pressure and 

skin friction ; therefore, evaluation of fluctuating turbulent wall pressure is very 

important. 

would require an appropriate unsteady turbulence model to be incorporated into the 

flow analysis, a few turbulence models are discussed below. 

Also in flow-solid 

Since successful numerical analysis of flow-solid interaction of tube banks 

Significant progress has been made in modeling steady turbulent flows. But due 

to pressure straining, redistribution of Reynolds stresses, and the complicated geometry of tube banks, some of the steady state turbulence models with proven performance 

may not even be applicable for steady turbulent flows through a nest of rigid cylinders. For example, low Reynolds number turbulence models which include the distance from the walls explicitly in the turbulence equations may hardly be applicable 

for flows through tube banks, since there exist so many cylinder walls in any tube 

bank, e.g., there are 600 tubes in the SSME Main Injector Assembly. On the other 

hand, numerical modeling of unsteady turbulent flows has shorter history than the 

steady case; consequently, there exist only a limited number of publications in this 

area and physical understanding of unsteady turbulent flows is not complete at the 

present time. Therefore, this discussion is limited to the most successful turbulence 

models which may be further extended to unsteady turbulent flows and which can be 

25 
I 

used for turbulent flow-solid interaction of tube banks. 

number k- E turbulence models, a few different versions of the algebraic Reynolds 

stress models, and the large eddy simulation method (LES) . be used for turbulent flows through tube banks at its present state of development as 

can be seen in the following discussions, but the method has been originally developed 

for and primarily used for unsteady turbulent flows. 

various turbulence models can be found in Nallasamy [ 691 , Martin [ 701 , Lakshminarayana 

[ 711 , and the references appearing in these literatures. 

These are a few low Reynolds 

The LES method cannot 

Comprehensive reviews of 

3.3.1. Extended k- E Turbulence Models and Algebraic Reynolds Stress Models 

For a class of complex flow problems, anisotropy of Reynolds stresses induce 

significant differences in the predicted flow field, e. g. , the existence of secondary 

flow in curved square ducts [ 721. For these flows, the conventional (or standard) 

k- E turbulence model yields unsatisfactory computational results. The Reynolds stress 

model [ 731 can describe the convection, diffusion, and anisotropy of each component 

of the Reynolds stresses. But use of the Reynolds stress model is prohibited in many 

cases of practical importance due to the number of variables which result from six 

differential equations for Reynolds stresses and one equation for turbulence dissipation 

function for three-dimensional flows. 

In order to account for the anisotrophy of Reynolds stresses, the algebraic 

stress model, which retains most of the basic features of the original Reynolds stress 

model, was developed [ 741 . The algebraic relationship between the six Reynolds stress 

components are obtained from the Reynolds stress equation by neglecting the convection and diffusion terms [76] or by assuming the convective transport of Reynolds 

stresses to be proportional to the transport rate of turbulent kinetic energy [ 741. 

Hence in algebraic Reynolds stress models, the k- E turbulence models are enhanced 

with this algebraic relationship instead of using the classical Boussinesq eddy viscosity 

assumption [ 761 . 

usually applied in two different ways. 

boundary of flow domain near the wall is slightly separated from the wall and the 

wall functions are applied at the first grid point. 

functions are not applicable near the separation points on the walls; nevertheless, 

they have been used successfully in many cases. 

grid point is located on the wall, no slip boundary condition is applied at the grid 

point, and the necessary modifications are included in the turbulence models in order 

to account for the low Reynolds number effect near the wall. 

In the numberical analysis of turbulent flows, the wall boundary conditions are 

In the first approach, the computational 

Strictly speaking, these wall 

In the second approach, the first 

In flow-solid interaction of a nest of cylinders subjected to cross flows, the 

vibration of cylinders is a three-dimensional motion. The capability of wall function 

methods for the analysis of flows with moving boundaries has not been established as 

yet, and need to be established so that the method can be used for the analysis of 

flow-solid interactions. Alternatively, the low Reynolds number turbulence model may 

be used so that the no-slip flow boundary conditions can be specified at the flowsolid interfaces. 

For unsteady turbulent flows, the mean and fluctuating quantities are defined 

The most familiar averaging technique for unsteady turbulent 

differently than those for steady turbulent flows. Therefore , it would be appropriate 

to include this topic. 

I flows would be the ensemble averaging method in the sense that it yields exactly the 
same flow equations as the time averaging method for steady turbulent flows. 

average of any fluctuating quantity, Y (x ,t) , is defined as: 

Ensemble 

N 

(3.25) 

n=O 

where T is the period of fluctuation, n is an integer, and F denotes ensemble average 

of Y . On the other hand, the time average for steady turbulent flows is defined as 

t+T 

(3.26) 

where denotes the averaged quantity. For both of the averaging techniques, the 

velocity component is decomposed as 

( 3.27) 

where yT = (vl,v2,v3) is a column vector of the three velocity components in the 

three spatial coordinate directions, and 

fluctuating velocities. 

flows are given as: 

is the corresponding column vector of the 

The time dependent turbulent flow equations for incompressible 

a G. 

J =o 

j 

ax (3.28) 

(3.29) 

7 

The Reynolds stress, v.'v ?, in equation (3.29) is modeled using the Boussinesq 1j viscosity assumption givcn as : 

(3.30) 

where v is the turbulent viscosity; 6.. is the Kronecker delta such that 6.. = 1 for 

1 = j, and 6ij = 0 for i # j; and k is the turbulent kinetic energy, k = 1/2 vilvit . t 11 lL 

27 
The standard two equation (k-e) turbulence model due to Launder and Spalding 

[ 781 is given as: 

where E is the dissipation rate of the turbulent kinetic energy, 

(3.31) 

(3.32) 

P is the production rate of the turbulent kinetic energy, 

Y 

2 = c k /E, is the same turbulent viscosity as used in equation (3.20); and vt' 9 1-I 

c 

1-I 

78. 

near the wall, where the molecular viscosity is dominant over the eddy viscosity, the 

wall function methods are used to supply wall boundary conditions for the differential 

equations [ 781. 

Reynolds number effects. 

Poreh 1791 is given as 

= 0.09, cEl = 1.44, cE2 = 1.92, ak = 1.0, and uE = 1.3 have been used in Reference 

Since equations (3.31) and (3.32) are not valid in the viscous sub-layer region 

There exist a few modified k- E turbulence models to include low 

The extended k- E turbulence model due to Hassid and 

k 

2 

Y 

at j ax ax 

2 

t k 

2 - [l - 0.3 exp(-R 11 5- - 2v 

28 

(3.33) 

(3.34) 

I 
2 2 where vt = c k /E, c = c f = 1 - exp (-0.0015Rt), Rt = k /(vE) is the turbulent Reynolds number, and c = 0.09. The model by Jones and Launder [ 801 is 

given as 

f 

u u lJw u' P 

lJm 

(3.35) 

(3.36) 

2 where c = c exp[-2.5/(1+Rt/50)], cE2 = c [1-0.3exp(-Rt )], c = 0.09, cEl = 

1.45, and cEw = 2.0, and rest of the notations are the same as in equations (3.33) 

and (3.34). The wall boundary conditions for the low Reynolds number turbulence 

equations are given as: 

lJ uw Em PO3 

k=o E=O 1 (3.37) 

Various algebraic relationships between Reynolds stresses can be obtained by 

introducing simplifying assumptions into the Reynolds stress equation [ 73,771 . simplifying assumptions are : production and dissipation of turbulent kinetic energy 

are in equilibrium, transport of Reynolds stresses are proportional to the transport 

of turbulent kinetic energy, transport and diffusion of Reynolds stress are negligible, 

and any appropriate combination of these assumptions. 

[73,77] and the algebraic Reynolds stress model due to Rodi [74] are given below. 

The Reynolds stress equation [73,77] is given as: 

These 

The Reynolds stress equation 

a v.1~ I 

ax 

avi' a v.1 

j - v.'v.'v 11 k ' - p (6 jki v' + Aik v+ 2v - axk - axk 1 , (3.38) 

k 

where (I) is the time rate of change, (11) is the convective transport, (111) is the 

stress production, (IV) is the pressure straining term, (V) is the diffusive transport, 

and (VI) is the dissipation of the Reynolds stresses. The closure model for higher 

29 
order terms for low Reynolds number flows proposed by Hanjalic and Launder [77] 

are given as 

a v.' a v.' c, + 8 

(P.. - - 6.. k) 11 3 9 6.. k) - 111 P(1 + L) = - c1 E E- (v.'v' - - 

ij 3 11 

j p axi ax 

- (30c2 55 - 2, k (F avi av. 8c2 - 2 2 

11 (D.. 11 - 3 PS..) 13 + o..,~ 11 

(3.39) 

E- 2 0.125 E (v.'v.' - - k6 ..) + 0.015 (P.. - Dij) 9ij , w 11 3 11 11 

1 a T. 

P.. = - v.'v.' - 1 + v.'v.' - 

11 I- 1 axk 1 1 axk 

+ v.'v 1 k ' - aTkxi 1 

I 11 

P 

P - - (6jkVi' -t bikVj') = 0 

(3.40) 

( 3.41) 

(3.42) 

(3.43) 

(3.44) 

(3.45) 

(3.46 

where P = -v 'v 

mean strains; fs = (1+Rt/l0)-'; Rt = k /(vE) is the turbulent Reynolds number; and 

c1 = 1.5, c2 = 0.4, and cs = 0.11 have been used in Reference 77. 

dissipation function, E, in equations (3.39) through (3.45) is determined by equation 

(aF,/ axk) is the rate of generation of turbulence energy by the 

The turbulence 

2 Rk 

30 
(3.46) which is the convection-diffusion equation of the turbulence dissipation function. The pressure straining due to wall proximity, equation (3.40), was proposed in 

Reference 73 but was disregarded in Reference 77. The form of wall proximity equation, equation (3.40), is not appropriate for numerical analysis of turbulent flows 

through the tube banks due to the explicit appearance of the distance from the wall. 

The algebraic Reynolds stress model derived by Rodi [74] is given below, where 

the convective transport term of the Reynolds stresses is assumed to be proportional 

to the transport of turbulent kinetic energy, and the diffusion of Reynolds stresses 

is neglected. Accordingly, 

(3.47) 

where y % 0.6, and the rest of the notations are the same as in the Reynolds stress 

equation. Successful application of the algebraic Reynolds stress models can be found 

in a number of publications [74,81,82]; hence, it is not elaborated herein. But its 

extension to unsteady turbulent flows remains to be tested. Computational time for 

using algebraic Reynolds stress models is comparable to the k- E turbulence models 

for most of the cases. 

Extension of the k-E turbulence models to unsteady flows has been limited to 

flows with simple geometries [ 833 and boundary layer flows. Even for these cases, 

straight forward extension of the k- E turbulence models has been more unsuccessful 

than successful. These model inadequacies motivated a group of researchers to devise 

a new triple decomposition model of turbulence [84], in which any fluctuating quantity 

is decomposed into three parts, i.e., time mean quantity, organized oscillations, and 

random fluctuations. However, this method also remains to be further studied to be 

useful for engineering calculations. 

A successful application of the standard k-E type turbulence model to unsteady 

recirculating turbulent flow can be found in Boyle and Golay [ 851. Measured velocity 

field and computed velocity field at 80 sec into transient of the scaled nuclear reactor 

outlet plenum are shown in Figure 11 [ 871. 

velocity components, at various locations inside of the plenum chamber, which vary in 

time are shown in Figure 12 [87]. The computational results compare favorably with 

experimental data during the flow development period, but the same results begin to 

diverge from experimental data as steady state is approached. 

Experimental and computational 

3.3.2 Large Eddy Simulation (LES) Method 

In turbulent flows, instability of mean flow generates large eddies and these 

large eddies break up into smaller eddies. After many cascade processes, see Figure 

13 [86], the fine scale eddies are dissipated by viscous forces. 

scale turbulence, or motion of large eddies, strongly depends on the type of flow and 

geometry of the flow field; whereas, fine scale turbulence, or motion of fine scale 

eddies, is almost universal for any turbulent flows. Therefore, fine scale turbulence 

can be modeled more easily than the large scale turbulence due to its universality. 

The large eddy simulation (LES) method best utilizes these physical observations. 

Until now, the LES method has been used for channel flows and axisymmetric annular 

flows. 

The scale of large 

31 
$*e. 

.... 

.... 

* 

. 

. Dbbbbb 

VELOCITY SCALE 

0- 1 .O rnlsec 

0.5 

(b) 

Figure 11. Turbulent recirculating flow in a nuclear reactor outlet plenum [851. 

(a) Measure mean velocity profile after 80 sec into transient. (b) Computed 

mean velocity profile after 80 sec into transient. 
20 60 100 140 

TIME (sec) 

20 60 100 140 

TIME (sec) 

Figure 12. Comparison of measured vertical velocity with computer 

velocity 1851, a through h in the Figure denotes spatial 

locations specified in Figure 10(a>. 

ENERGY PRODUCTION 

QOUOW-Q~OQ 

oo~o~oooooeoooooo~ 0 00 

DlSl PAT ION 

Figure 13. The energy cascade process [86]. 

33 
In the large eddy simulation method, any flow variable is averaged over a finite 

volume which is determined by the grid size used in computations. 

called filtering. Let Y be some flow variable, then Y can be decomposed as: 

This process is 

y=y!+y' (3.48) 

where Y ' is called the residual field, and is the large scale field defined as: 

(3.49) 

where AV = AxlAx2Ax 

direction, and the integration is performed over the volume AV. 

tion, G(_X,~'), in equation (3.49), was first introduced by Leonard [87] in order to 

account for nonnegligible variation of large scale motion inside of the averaging volume, 

AV. 

Axk(k=1,2,3) is the grid size in the k-th spatial coordinate 

The filtering func3' 

The filtering function is given as: 

(3.50) 

where Gi(i=l, 2,3) is the distribution function. 

function is the Gaussian distribution function given as : 

The most frequently used distribution 

22 {-6(xi-xi'> /A I 

G.(x.,x.') = 4; e 111 (3.51) 

where A is twice the computational mesh size in xi-coordinate direction. 

ent forms of filtering functions can be used in practice. 

was used in the first application of the LES method in Deardroff [88]. 

the Navier-Stokes equation used in LES is given as [88] : 

Many differA flux form of 

For example, G(x,,x') = 1 

(3.52) 

Filtering equation (3.52), using any filtering function G(x ,x') , yields : 

34 
I a T~ a t. 

- 

at ax 

+- a (E)=-- ax 11 a X. j 1 

a T. J=o 

j 

ax 

I (3.53) 

-- --- - - where v.v -- - = v.v.+nij, qij - v.'v. + v.v.' + v.'v.', pj = -(nij-nkk6ij/3), and 6.. is the 

Kronecker delta. 

residual stress, v.'v.', due to Deardroff [88] is given by: 

1j 11 13 13 11 11 

A subgrid-scale-model (closure model for LES equation) of the 

- 13 

- 

L- A.. = v.v - 7.T = o 

9 1j 1j 

a X. 

(3.54) 

(3.55) 

- 

where 7.v ' = 

[ 891 proposed cA = 0.17, Deardorff [ 881 used cA = 0.1, and Moin and Kim [ 901 used 

= 0, and cA is a constant to be determined experimentally. Lilly 1j Ji 

c 

term to compensate for the same term appearing in the T.. expression so that equation 

(3.55) holds when the subscripts i and j are contracted. in equation (3.53) need to be modeled separately, and Lilly [89] proposed: 

= 0.065. In equation (3.53), nkk6ij/3 has been added to the pressure gradient A 

11 

Consequently n kk 

(3.57) 

where c1 = 0.094. 

The governing equation for pressure distribution can be obtained as follows. 4 

By letting P = p/p+nkk/3 and taking the divergence of the momentum equation, equation ( 3.53), yields: 

(3.58) 

where dissipation due to molecular viscosity has been neglected and the time derivative 

term can be dropped out using the conservation of mass condition. 

35 
Several different methods have been used by different workers [ 88,91,92] to 

solve the LES equation, and the form of the LES equations were slightly modified to 

be suitable for each of the numerical methods used. In Deardorff [ 881 , the classical 

finite difference method was used to solve the LES equation as given in equation 

(3.53); in Schumann [ 911, the control volume finite difference method [go] was used 

and any time derivative term appearing in the Navier-Stokes equations was filtered 

over the control volume as before whereas the rest of the terms appearing as fluxes 

in the control volume finite difference method was filtered over each of the surfaces 

of the control volume where fluxes are evaluated. 

spectral method [63,64] and the LES equation of the rotational form of the NavierStokes equations, i.e., equation (3.23), were used. Particularly in Moin and Kim [92] , the LES equation was refined to include the exact treatment of the Leonard's stress, 

i. e. , equation (3.54), and molecular viscosity. 

In Moin and Kim [ 921 , the pseudoThe computational results given in these references [ 86,88,92] compare 

excellently with experimental data. But it should be noted that the constant coefficient cA , used in equation (3.56), as proposed by Lilly [ 891 , is larger than twice the 

value used in Moin and Kim [ 921 . Hence, the eddy diffusity, which is proportional 

to the square of c 

different. Conceptually, the subgrid-scale-model should not dominate the computational results when fine grids are used for numerical analysis, and the residual stress 

itself should vanish as computational grid becomes infinitely small. But the aforementioned computational results, using practical grid size, showed that the subgridscale model is important to render the computational results comparable to experimental data. 

made for the LES method to be applicable to practical engineering calculations. 

used in these calculations for the similar channel flows is quite A' 

Therefore, further refinement of the subgrid-scale-model need to be 

3.4 Structural Analysis Methods Including Large and /or Permant Deformations 

Due to fluidelastic instability, thin cylinders in tube banks may go through 

large deformations which may cause viscoplastic deformations of the cylinders. 

fore, the structural analysis method to be used for flow-solid interaction of tube banks 

need to be capable of modeling elastic vibrations as well as large and/or permanent 

viscoplastic deformations, in order to predict Bauschinger effect, hysteretic effect, 

and cyclic straining. But to date, no specific numerical structural analysis method 

designed for this purpose has been developed. Structural analysis methods which 

can be further developed to meet these requirements are discussed in this section. 

ThereMost of the structural analysis methods are based on the finite element displacement method. 

mization of energy functional or derived from the virtual work principle. 

problems are usually described on the frame work of Eulerian coordinate system; and 

hyperelasticity, large deformations, and viscoplastic deformation problems are usually 

described on the convected coordinate system in order to account for large displacements. 

principle is described on the Lagrangian coordinate system which is updated at each 

of the solution time levels. The virtual work principle and the minimization of energy 

functionals are identical for elasticity problems ; but the virtual work principle is more 

convenient for geometrically and /or materially nonlinear problems, since derivation of 

energy functionals for nonlinear problems can be avoided. 

the virtual work principle and the convected coordinate system are used. 

The finite element system of equations are derived through the miniElasticity 

In the methods utilizing the convected coordinate system, the virtual work 

In the following discussions, 
For time dependent deformation of structures, consider a point initially at 

xi(i=l, 2,3) moved to ui( 1=1,2,3) in the same Cartesian coordinate system. In Lagrangian formulation, the deformations are described in terms of the convected coordinates 

u.(i=l, 2,3). The force equilibrium equation, or the momentum equation, for a structurd system is given as 

1 

0.. + b. = 0 

11 ,j 1 

c u(5,t) = u_,(t) on aR1 

(3.59) 

(3.60) 

0.. n = Ti(t) on aR2 (3.61) 11 j 

where bi is the body force; 52 is the structural domain; a Rl is the boundary on which 

the Dirichlet boundary condition is specified; a R2 is the boundary on which surface 

traction is specified; aRln aR2 is a null set, aRIU an2 = an, and aR is the entire 

boundary of the domain; both the surface traction Ti(t) and the prescribed displacements uo(t) are functions of time in general; and 2 is an outward normal vector on 

the boundary. 

momentum equation as a negative body force. 

The inertial force due to acceleration is usually included into the 

The virtual work statement is given as: 

(3.61) 

where 6 ui(i=l, 3) is the virtual displacement which satisfies the Dirichlet boundary 

condition, equation (3.60), on the boundary 8 Rl. Integrating by parts the stress 

term and using the fact that 

since 6 ui = 0 on a R in equation (3.61), yields : 

(3.62) 

which is the virtual work expression to be used in the finite element method. 

finite element analysis, it is more convenient to write the stress tensor, 

In 

“ij, and the 

37 
T- I, and 5 - {“xx ? “yy “zz Oxy 9“yz OZX 

strain tensor, Eij, into column vectors zT = 

I. The shearing strain tensor, E..(i#j), is not equal to the iEXX “yy , Ezz’YXy’YyZ ’YZX 

engineering shearing strain, y..(i#j), hence adequate modification need to be included 

13 

in the material property coefficients in the stress-strain relationship. 

matrix notations, a..6u. can be expressed as 65 2. Let 

13 

Using the 

T 

11 1,j 

(3.63) 

where N is a matrix of interpolating polynomials, 5 is the column vector of nodal displacements, & is a differential operator, and ,B = SY. 

problems, L and B depend on displacements. Introducing equations (3.63) into equation (3.62) and taking the variation of equation (3.62) with respect to the nodal displacements yields, in matrix notation, 

gs 

For geometrically nonlinear 

s 5=z 

/BTodx_-/3T>dx- c f NTJdS=O c= . 0 aR2 aC 

(3.64) 

For geometrically and/or materially nonlinear problems, equation ( 3.64) is not 

generally satisfied at any computational stage. 

incremental form yields 

Hence, writing equation ( 3.64) in 

(3.65) 

where dz = 0 if the surface traction is not incremented at any computational stage. 

Equation (3.62), or equation (3.65), represents the most general form of the 

finite element equation for any spatial dimensions; however, its specific forms for 

different problems and different solution methods are so diverse that they cannot be 

covered in this review. These different problems include : different strain/stress 

assumptions ; geometrical differences such as one dimensional bars, two dimensional 

plates, and quasi three dimensional shells, etc. ; and the geometrical and/or material 

nonlinearities . interpolation methods such as 

Lagrangian interpolating polynomials and the Hermite interpolating polynomials ; 

different formulations such as mixed formulation and displacement formulation ; and 

even time integration methods. 

easily for the flow-solid interaction analysis are covered herein. 

elasto-viscoplasticity analysis method of Hughes and Taylor [ 931 and the nonlinear 

Different solution methods include : 

Only a couple of algorithms that can be extended 

These are the 

38 
vibration analysis method, which includes large deformation and elasto-viscoplastic 

deformations, due to Rodal et al. [99,100]. 

The most rigorous mathematical treatment of viscoplastic deformation can be 

found in Hughes and Taylor 1931 and Owen and Hinton [94]. In References 93 and 

94, a linear strain-displacement relationship has been used and the acceleration force 

has been neglected, on the assumption that it is negligible compared with forces which 

cause plastic deformation. Hence, the method is suitable for problems in which only 

the permanent plastic deformation is important, but not suitable for problems in which 

the vibrational aspect is as important as the permanent plastic deformations. 

general form numerical analysis method for viscoplasticity due to Hughes and Taylor 

[ 931 and Owen and Hinton [ 941 is given below. 

The 

In elasto-viscoplastic analysis, it is assumed that the total strain rate, g, can 

be decomposed into an elastic part, %, and a viscoplastic part, zvp, i.e., 

(3.66) 

where the over-dot, ('1 notation denotes time rate of change. 

rate due to Perzyna [95] is given as: 

The viscoplastic strain 

(3.67) 

where y is the fluidity of the plastic flow; <$(F)> implies <$(F)> = $(F) if F > 0, and 

<$(F)> = 0 if F < 0; F is the yield function such as the Huber-von Mises yield function [ 551 , and Q is the plastic potential function. 

identical to Q , but different for nonassociated plasticity. 

to the elastic strain rate as: 

For associated plasticity, F is 

The stress rate is related 

or, in incremental form, 

( 3.68) 

(3.69) 

where the superscript n denotes the time level t = tn, A denotes increment, and the 

viscoplastic strain increment, A E can be interpolated as 

The increment is defined for a time interval from tn to tn+l, i.e., Atn - tn+l-tn. 

n 

VP ' 

(3.70) 

39 
where e is an interpolation parameter such that 9 = 0 corresponds to the fully explicit 

time integration method, 8 = 1 corresponds to the fully implicit method, and e = 1/2 

corresponds to the Crank-Nicholson type method. The time integration method for the 

present elasto-viscoplasticity analysis is unconditionally stable for o 2 1/2. 

remark, numerical stability does not imply that convergent solutions can always be 

obtained for many nonlinear problems. 

of stress state, equation (3.70) can be rewritten as 

As a 

Since the viscoplastic strain rate is a function 

(3.71) 

where (E .n+l - E -n ) = (aEvp/aa)nA%n has been used, and cn = BAt (a&vp/3u)n. The 

incremental strain-incremental displacement relationship can be written as : 

-vp "vp n- 

(3.72) n n As =B Aa 

Fwhere ,B is a constant matrix for infinitesimal strain-displacement relationship. 

stitutincg equations (3.71) and (3.72) into equation (3.69) yields: 

Sub- 

(3.73) 

n -1 where I! = (D-l + C ) . Equation (3.73) is the incremental constitutive equation 

to be <sed inqhe inGrementa1 virtual work principle, i.e., equation (3.65). 

material nonlinearity only, the incremental virtual work equation is given as : 

For 

Inserting equation (3.73) into equation (3.74) yields 

(3.74) 

(3.75) 

The incremental displacement, A an is obtained by matrix inversion of equation (3.75) 

as 

BD kn At dx - AT n -1 Aan = [E,] + e 

n where i,KT = BT sn B is called the tangential stiffness matrix. 

40 

F: 9s 

(3.76) 
In Hughes and Taylor [931, the above methodology was used to solve an elastoviscoplastic deformation of a thick wall cylinder subjected to constant internal pressure. 

The physical problem was reduced to a one-dimensional problem by assuming an 

axisymmetric plane strain state. The physical problem and the computational results 

are shown in Figure 14. 

2 1,000 

INTERNAL PRESSURE VS TIME 

THICK-WALLED CYLINDER 

flur - -0.3 

2Y -0.4 0.4 

-0.5 

-0.6 

-0.05 

RESULTS AT GAUSS -0.10 

POINTS 1.0 1.4 1.8 

a = 1.0 t -= t = 106 

RESULTS FOR TIME STEPS 

r 

FINITE ELEMENT 

r 

2-12 ARE IDENTICAL 

(b) 

Figure 14. An elasto-viscoplasticity analysis of a thick cylinder subjected 

to internal pressure [93]. (a) Description of physical problem and 

a finite element discretization. (b) Comparison with 

exact plasticity solution. 

41 
The numerical analysis method given in equations (3.66) through (3.76) are so 

Yet , specific forms 

For 

general that it can be used for any elasto-viscoplastic problems. 

of strain-displacement relationships for individual problems need to be derived. 

flow-solid interaction of tube banks, a general shell element would be preferred. 

Refinement of the lasto-viscoplasticity analysis methods in References [ 931 and 

[ 941 can be made in several different ways. These refinements would be: incorporation of temperature dependency of the viscoplastic strain rate into equation ( 3.67) ; 

inclusion of geometrical nonlinearity ; inclusion of work hardening into the yield function; and inclusion of the acceleration force term into the governing equation so that 

Bauschinger effect , hysteretic effect , and cyclic strain hardening can be predicted. 

Examples of temperature dependent viscoplastic strain rate can be found in References 

96 and 97, in which flow formulation was used and the effect of elasticity was completely neglected. Examples of geometrical nonlinearity and work hardening can be 

found in References 98, 99, and elsewhere. 

Rodal et al. [99,100] developed a finite element structural analysis method which 

includes finite strain, large deformation , and elasto-viscoplasticity so that Bauschinger 

effect, cyclic straining, and hysteretic effect could be predicted. In this method, all 

the nonlinearities such as geometrical and material nonlinearities were incorporated into 

the load vector term so that the contributions of the nonlinearities appear only as 

modifying forces to the applied load. 

et al. [ 99,1001 is given as 

The structural dynamics equation due to Rodal 

where M is the usual mass matrix, - K is the elastic stiffness matrix, 

1 vector of externally applied loads, and is a column vector of modifying loads due 

to geometrical and material nonlinearities. The computational procedures used in 

References 99 and 100 are simpler than the one used in References 93 or 94, in which 

the viscoplasticity has been treated more rigorously. Undoubtedly , the numerical 

method due to Rodal et al. [ 99,1001 would be adequate for modeling short period 

vibrations such as those associated with impact problems; but the accuracy of the 

method as applied to structures vibrating for a longer period of time than the impact 

duration need to be tested further. 

is a column 

cc - 

The time dependent elasto-plastic finite element method [57] used in the laminar 

flow-solid interaction analysis of Donea et al. [51] and Belytschko et al [52,53] are 

described below together with the Newmark time integration method [loll. 

implicit Newmark time integration method is suitable to achieve a rigorous force equilibrium condition in integrating a system of second order ordinary differential equations in the time coordinate direction. 

to be tested yet. 

publications including References 55 and 102, and are not included herein. 

The 

But its adequacy for flow-solid interaction has 

Elasto-plasticity constitutive equations can be found in a number of 

Belytschko and Hsieh [ 571 developed a time dependent elasto-plasticity finite 

element method. 

displacement relationship) was used and the large dispalcement was taken care of by 

using the convected coordinate system. 

57 is given as: 

In this method , infinitesimal strain assumption (a linear strainThe finite element equation used in Reference 

42 
(3.78) 

which is equivalent to equation ( 3.64) , but the acceleration term appears explicitly. 

The semi-discretized (fully discretized in spatial coordinates, but left continuous with 

respect to the time variable) finite element equation becomes 

In Belytschko and Hsieh [57], the mass matrix was lumped into a diagonal matrix and 

the elasto-plastic stiffness term was incorporated into the load vector term as a modifying load vector. 

tions was solved explicitly in time. 

methods are suitable to solve short-period structural responses such as the one for 

impact problems, or to solve transient problems which asymptotically approach steady 

state. The force equilibrium condition may not be satisfied rigorously at each time 

step. The implicit Newmark time integration method which can satisfy force equilibrium 

condition more rigorously is briefly described below. 

The resulting system of second order ordinary differential equaAs discussed before, explicit time integration 

In the Newmark time integration method [loll, the displacement vector in equation (3.79), which is a continuous function of time, is expanded into a parametric 

series form such that 

-n+l - n+l n .- n M a - co(g - c. a - c1 kn - c2a, 

.n+l -n .- n+ 1 

a = c a + [(l - Cl) in + 152 I Atn 

where 

c2 = 1/27 - 1 

( 3.80) 

and 01 and y are constant such that (a,y) = (1/2,1/4) corresponds to the unconditionally stable, constant-average-acceleration time integration method. 

equations ( 3.80) into equation ( 3.79) yields 

Substituting 

43 
(3.81) 

which can be solved for displacements for each new time-level. 

3.5 Numerical Solutions of Flows Through a Nest of Cylinders 

Only very recently, a group of researchers began to numerically model the fluid 

flow through tube banks. The recent developments of flow analysis method for tube 

banks are discussed in this section. These are: control volume finite difference 

analyses due to Singhal et al. [44], Wambsganss et al. [103], and Sha and Yang [lo41 ; 

two-dimensional finite element analysis due to Zienkiewicz et al. [105]; and body fitted - 

grid finite differeme analyses due to Kwak and Chang [ 1061 , Chang and Kwak 

11071, and Rogers et al. [ 1081. In all of the above studies, the flow fields were 

assumed to be steady and the viscosity is constant. 

3.5.1 Control Volume Finite Difference Analyses of Steady Laminar Flows Through 

Tube Banks 

For tube banks with a large number of cylinders, detailed numerical modeling of 

the flow field is prohibited due to the present computer capability including available 

memory and computational speed. 

the concept of flow through porous media were proposed by Sha et al. 11041, 

Wambsganss et al. [ 1031 , and Singhal et al. [ 441. In these methods, the flow domain 

is discretized into a number of cells where each cell contains several cylinders (Fig. 

15). Then the resistance of each cell to fluid flow is represented by volume porosity 

(fraction of volume occupied by the fluid in the cell) and surface permeability (fraction of open projected flow area in the direction of fluid flow in the cell). 

governing equations are the same Navier- Stokes equations , but any volumetric quanquantity needs to be multiplied by the volume porosity; and each velocity component, 

by the surface permeability in the direction of each velocity component. 

field shown in Figure 8 [44] was obtained by solving the flow equations for porous 

media using the control volume finite difference method [ 901. 

Approximate numerical analysis methods based on 

The 

The flow 

In the fluidelastic instability analysis of Wambsganss et al. [1031, the computed 

flow velocity normal to the cylinder walls was taken as reference velocity for predicting the fluidelastic instability using the instability diagram shown in Figure 4. The 

validity of these methods remains in question; since fluidelastic instability is caused 

primarily by flow-solid interaction, and in most of the cases, the most serious fluidelastic instability is confined to the flow inlet region of tube banks according to 

experimental observations. 

More recently, a new computational method, called the control-volume based 

finite elements, was developed by Baliza and Patankar [ 1161 , Prakash and Patankar 

[ 1171 , and Prakash [ 1181 among many others. 

finite elements to discretize the domain and to interpolate the solution vectors. The 

discrete system of equations is derived by considering the flux through the element 

sides as in the control-volume based finite difference method. 

the standard finite element methods, the discrete system of equations are derived 

using the method of weighted residuals (MWR) to be introduced in the following section. The mathematical foundations of the two different methods are quite different ; 

and consequently, the final discrete systems of equations are also quite different. 

The method utilizes isoparametric 

On the other hand, in 

44 
:UT 

"2 

Figure 15. Control volume finite difference analysis of flow through 

a porous media [ 1031. (a) Discretization of a tube bank. 

(b) Configuration of a cell. 

Until now, the control-volume based finite element method utilized only the linear isoparametric elements. 

lines connecting the nodal points. 

methods for structural analysis utilizes quadratic or higher order elements to best 

approximate the geometry as well as solution vectors. 

finite element method for flow-solid interaction analysis in the future remains as an 

open question. 

In this case, the curved geometry is approximated by straight 

On the other hand, most of the finite element 

Use of the control-volume based 

45 
3.5.2 Finite Element Analysis of Steady Two-Dimensional Laminar Flows Through 

Tube Banks 

The finite element method has certain advantages over the finite difference 

method. capability to include various types of boundary conditions more naturally and rigorously in the formulation ; capability to show convergent 

nature of the method mathematically ; and capability to model complicated domains more 

precisely than the finite difference method. 

These advantages are: 

Application of the finite element method for flows has shorter history than the 

finite difference method, and the finite element method is still in its early development 

stage as far as application to fluid flows is concerned. Much of the development time 

for the finite element method for flows has been spent in overcoming numerical instability arising from convection in high Reynolds number flows. 

development of the Petrov-Galerkin method [ 1091 for convection dominated flows, more 

rapid development of the method for flows is expected in the near future. 

Petrov-Galerkin method and its application to two-dimensional steady laminar flows 

through tube banks are discussed in this section. 

But with the recent 

The 

The Bubnov-Galerkin method [lo91 (in which test functions are selected from 

the same space as the trial functions) yields exactly the same discrete equation as 

the central differencing finite difference method for one- dimensional problems discretized using linear elements. As a remark, for most of the cases, the discrete equations obtained through the two different methods are quite different. 

it was believed that the Bubnov-Galerkin method (or equivalently, the central differencing scheme) is the most accurate discrete representation of the continuous differential equations. But, Brooks and Hughes [ 1091 showed that : the Bubnov-Galerkin 

method underestimates the numerical diffusion ; the conventional Petrov-Galerkin method 

(or equivalently, the backward differencing scheme) overestimates the numerical 

diffusion; and the optimal Petrov-Galerkin Method (this should be considered as the 

consistent Petrov-Galerkin method when source terms exist in the differential equation) 

yields the most accurate numerical results. A finite element solution of one-dimensional 

convection-diffusion equation of the form : 

For a long time, 

in xe(0,l) a 2T 

ax 

koT =O aT 

v-- 

ax (3.82) 

with boundary conditions T( 0) = 0 and T( 1) = 1, is shown in Figure 16, where T is 

the temperature. 

is given as [lo51 

The basis of test functions for the consistent Petrov-Galerkin method 

a 9i 2 Oi*(X) = Cpi(X> + - ax v koT'/lvI , i = 1, N e (3.83) 

(3.84) 

46 
EXACT 

-----_-- CENTRAL DIFFERENCE 

--- UPWIND DI FFER ENCE 

Figure 16. A finite element solution of one-dimensional 

convection-diffusion problem. 

where a = coth(Pe/2)-2/Pe, Pe = Ivlh/k is the grid Peclet number, v is the convection velocity, @i(x) is the shape function, m = 2 for linear elements, and Ne is the 

number of nodes in an element. 

flows as shown in the two-dimensional flows through tube banks. 

method to three-dimensional cases was suggested to be the same [ 1051. 

This method naturally extends to two-dimensional 

Extension of the 

For a two-dimensional case, the basis of the test fucntions of the consistent 

Petrov-Galerkin method is given as [ 1051 : 

(3.85) 

where zT = (x,y), kO" is the same as in equation (3.84) but IvI = (v12+~22)1/2 at 

each of the Gauss points used in numerical integration, and h is the length scale of 

the element. 

The finite element equation for Navier- Stokes equations is obtained by the 

method of weighted residuals. In the method: the flow domain is discretized into a 

number of elements; the Navier-Stokes equations are multiplied by the test functiqn; 

the second order diffusion term is integrated by parts; continuous variables are interpolated using the nodal values of the variables and interpolating polynomials; the 

resulting equations are integrated over each element using the Gauss numerical quadrature to obtain element system of equations; and the element system of equations are 

assembled together to obtain a global system of equations. A standard finite element 

system of equations for an element is given in the following. 

47 
K 

2K +K 

2E 11+5 22 s 12 

s 22 P11 K p21 

where 

r 

Yl 

V -2 

(3.87) 

(3.88) 

In equation (3.88), i = 1, Ne; j = 1, Ne; k = 1, Ne; N 

nodes in an element; n = 1, N 

integration is performed over each element, and $. 

is the number of velocity e 

N is the number of pressure nodes in an element, P' P 

denotes 3 +i/3x, etc. 

1 ,x 

In practice, the penalty finite element method [119] has been used quite frequently for flow problems (or more frequently for some problems such as polymer 

extrusion problems [llO]). In the penalty method, the unknown pressure is preeliminated from the Navier-Stokes equations. 

obtained in the post process if necessary. 

shown in this section was also solved using the penalty method. 

pressure variable, we set : 

The pressure field solution can be 

The laminar flow through tube banks 

To eliminate the 

v v = Pla 

P e- 

(3.89) 

48 
where p is the pressure, and CY is the penalty parameter. As CY approaches infinity, 

2 y % 0 and, hence, the continuity condition is satisfied approximately. In these 

7 computations, ~1 is taken to be 10 through 10" times the molecular viscosity to avoid 

ill-conditioned system of equations. Finite element discretization of equation (3.89) 

yields 

P P 

P 

(3.90) 

where Q 

(l,N ), and N is the number of pressure nodes in an element as before. Solving 

for pressure, p, in equation (3.90) and substituting the result into equation (3.86) 

yields penalize3 discrete Navier- Stokes equation to be solved for nodal velocities 

only. 

and s2 are the same as in equation (3.88), = j$R$ndx, R = (1,N ) , n = 4 P 

P P 

The penalized discrete Navier-Stokes equations are given by 

V 

FF -2 -2 

(3.91) 

Equation (3.91) is (2Ne by 2Ne) system of equations, and the necessity for solving 

the conservation of mass equation and the pressure variable has been eliminated. 

desired, the pressure field can be recovered by solving equation (3.90) in the post 

process. One advantage of the penalty method is that the number of unknowns are 

reduced so that a slight computational efficiency is achieved. But Zienkiewicz [ 1101 

has pointed out that solving the standard discrete Navier-Stokes equations can achieve 

comparable computational efficiency if the frontal solution technique [ 1111 is used and 

the pressure variables are eliminated at element levels. But in this case, the pressure 

needs to be discontinuous across element interfaces. The standard solution method 

would be preferred in the flow-solid interaction analysis since the pressure solution is 

required for most of the cases. 

dimensional steady laminar flow through a tube bank is shown in Figure 17, where 

the Reynolds number is 500. 

If 

A Petrov-Galerkin finite element solution of twoMuch progress for the Petrov-Galerkin type analysis methods has been made 

recently for unsteady two-dimensional flows. the Taylor-Galerkin 

method due to Zienkiewicz et al. [ 1051, Donea 11121, Lohner et al. [ 1131, and Donea 

et al. [ 1141, in which Taylor series expansion of the time derivative term is made 

prior to the spatial discretization, and the Taylor series expanded Navier-Stokes equations are solved by the Bubnov-Galerkin method; and direct extension of the PetrovGalerkin method for the two-dimensional steady flows into unsteady two-dimensional 

flows by considering the time-coordinate as an extra spatial coordinate due to Yu and 

Heinrich [ 1051. 

dimensional unsteady flows can be extended to three-dimensional unsteady flows. 

Among these are : 

It is believed that most of the upwinding techniques for two49 
Figure 17. A Petrov-Galerkin finite element solution of a two-dimensional 

steady laminar flow through a tube bank. (a) Pressure contour 

lines. (b) Velocity vectors. (c) Stream lines. 

Finite element computations of turbulent , separated flows also begin to appear 

in the literature. Among these are Taylor et al. [ 1221, Smith [ 1231 , and Benim and 

Zinser [ 1241. In these papers, the standard k- E turbulence model for high Reynolds 

number flows has been used. In Taylor et al. [ 1221, the flow domain was discretized 

using quadratic quadrilateral isoparametric elements as shown in Figure 18a. 

numerically obtained mean velocity and the turbulence intensity profiles are compared 

with experimental data [125] as well as the finite difference solutions [125] in Figures 

18b and 18c. It can be seen that both the finite element method and the finite 

difference method predict the mean velocity profile accurately ; however, it can also 

be seen that the finite element method predicts the turbulence intensity profile better 

than the finite difference method for this particular case. 

The 

The boundary element method, which has been tested for incompressible laminar 

flows [ 1271, also exists. In this method, the residual function, defined in the same 

way as in the finite element methods, is weighted against the test function and integrated over the domain. 

tion of the governing differential equation, and hence the domain integration is reduced 

into a boundary integration for Laplace's equation. 

But in this method, the test fucntion is the particular soluBut still one needs to integrate 

50 
20.0 

17.0 

15.5 

10.0 

I 1 II I I I J+ 

-30.0 -25.0 0.5 10.0 30.0 00.0 120.0 180.0 

(a 1 

0 2 4 6 8 12 

h 

0 2 4 6 8 12 

0: EXPERIMENTAL DATA [125], X: FINITE ELEMENT COMPUTATIONAL RESULTS, - -- 

FINITE DIFFERENCE SOLUTION 

Figure 18. Flow over a backward-facing step. (a) Finite element 

discretization. (b) Velocity profile. (c) Turbulent 

kinetic energy profile. 

51 
over the domain for Poisson type equations. 

only when the particular solution for the governing differential equation is available. 

For nonlinear problems, the nonlinear terms can be considered as forcing functions in 

the iterative solution procedure. 

exhibiting singularities as those associated with fracture mechanics. 

the method for turbulent flows is not certain, not to mention the flow-solid 

interactions. 

Originally, the method was applicable 

The method is particularly suitable for problems 

But the use of 

3.5.3 Body-Fitted Grid Finite-Difference Analysis of Steady Flows Through Tube Banks 

In the body-fitted grid finite-difference method, the flow domain is mapped onto 

an orthogonal computational grid; the flow equations on the computational grid are 

obtained by transforming the conservation of mass equation, the Navier-Stokes equations, and the energy equation using the coordinate transformation relationship. The 

resulting system of partial differential equations are solved on the computational grid 

system. 

If the physical flow domain is complicated or highly curved, then the numerical 

results may not be accurate due to large truncation error involved in the grid transformation relationship. This difficulty was overcome in Kwak and Chang [ 1061 and 

Chang and Kwak [lo71 by dividing the flow domain into several sub-domains which 

overlap at the interface regions. 

slight computational efficiency can be achieved by reducing the size of the system of 

equations to be solved. 

quite successful for laminar flows, but its extension to turbulent flows is not so 

obvious as yet. 

may accumulate in an iterative solution procedure of coupled flow and turbulence equations so that the computational results may not be so accurate. 

fitted grid finite-difference method for complicated turbulent flows may depend on 

overcoming the truncation error, one method for doing so is the sub-domain method 

[ 106,1071 . 

Thus, grid transformation error can be reduced and 

The body-fitted grid finite-difference method proved to be 

It may be expected that the truncation error in grid transformation 

Success of the bodyIn order to achieve further computational efficiency, a pseudocompressibility 

method was used in Kwak and Chang [ 1061, Chang and Kwak [ 1071, and Rogers et al. 

[log], so that a time dependent solution method can be used to obtain a steady state 

solution. It was shown in References 106, 107, and 108 that the pseudocompressibility 

method takes less computational time than solving the steady state Navier-Stokes 

equations directly. 

two rows of cylinders, using the above described method, are shown in Figure 19 11081. 

Numerically solved grids and flow fields through a cylinder and 

The body-fitted grid finite-difference method is also intended to be further 

developed for flow-solid interaction analysis of tube banks [ 106,107,1081 . 

IV. CONCLUSIONS AND RECOMMENDATIONS 

This report provides a review of various methods used for design of tube banks 

as well as computational methods to analyze flow-solid interactions of tube banks subjected to cross flows. 

Some conclusions can be made regarding the present state of understanding of 

the flow field in tube banks and the flow-solid interactions in such tube banks. These 

are listed below. 

52 
N 

0 

. ..+-- 4- 

. ... *A '. . 

Figure 19. INS3D computation of flow through cylinders. (a) Grid for a single 

post on a flat plate. (b) Grid spacing near post-plate junction for a single 

post. (c) Velocity vectors for multiple post in z/d = 0.01 plane. 

(d) Grid for calculating flow around two rows of cylinders. 

(e) Velocity vectors for single post in z/D = 0.01 plane. 

53 
1) No major experiment has ever been made for the flow field in the region 

between tube banks and all the understanding of the flow field depends on conjectures 

by researchers in this area. 

2) Analytical models for flow-solid interaction in tube banks are based on too 

Consequently, practical design of tube banks mostly 

many restricting assumptions so that these analytical models cannot fully describe the 

practical problems of interest. 

depend on experimental data. 

3) Most of the analytical models use the Euler beam equation to represent thin 

cylinders. The validity of modeling thin cylinders by the Euler beam equation need 

to be studied further if these models are going to be extended to include large 

deformations and material nonlinearity. 

4) The mathematical model developed by Chen [45] is one of the most complete 

models relating to this subject, in the sense that it includes multiple number of 

vibration inducing mechanisms. Many of the coefficients, such as the fluidelastic 

stiffness coefficients used in the model, need to be supplied by experiments or by 

computational fluid dynamics as suggested by Chen [ 451. But the model can predict 

only weakly coupled flow- solid interactions. Its usefulness for predicting failure 

mechanisms such as cyclic straining or hysteretic effects is limited in its present 

form because of the reasons stated in (3). 

5) Detailed numerical modeling of the flow field in the region between tube 

banks has just begun and the numerical solutions obtained until now are for steady, 

constant-viscosity flows only. 

Solving flow-solid interaction problems requires, in general, capabilities to 

solve fluid flow problems, structural dynamics problems, and coupling the flowstructure interface conditions. Therefore, the state of development for turbulence 

models, numerical methods such as the finite element method and the finite difference 

method, nonlinear structural analysis methods, and existing analytical methods applicable to laminar as well as to turbulent flow-solid interactions have been included in 

the discussions in order to evaluate the possibility of initiating a development of 

major numerical analysis method for cross flows. 

It was found that the existing numerical analysis methods for flow-solid interactions are designed for specific problems of each researcher's interest. At this 

time, it is not known whether these methods can be naturally extended to flow-solid 

interactions for a nest of cylinders subjected to turbulent flow at high Reynolds 

numbers. 

The development stage, advantages, and disadvantages of the various numerical 

methods discussed in this report are summarized in Table 4 as they are related to the 

flow - solid interaction analysis. 

Based on these observations, some recommendations can be made. These are 

listed below. 

1) 

element method. 

enables the method to be applciable to high Reynolds number flows which in the past 

had been solved mostly by the finite difference method. It is, therefore, recommended 

that the use of the finite element method for flow-solid interaction analysis in tube 

banks be pursued further. 

Presently, most of the structural analyses have been based on the finite 

Recent developments in the Petrov-Galerkin finite element method 

54 
TABLE 4. COMPARISON OF DIFFERENT NUMERICAL METHODS FOR THE 

NUMERICAL ANALYSIS OF FLOW-SOLID INTERACTIONS 

(Smaller number denotes higher rating) 

CV-FDM' BFG-FDM2 CV-FEM3 FEM4 BEM5 FDM6 

Capability to handle 2 3 2 1 4 1 

nonlinear equations 

Capability to handle 3 

complicated geomet ry 

Number of grids 3 

required 

3 

3 

3 2 1 3 

3 2 1 3 

Experience gained 1 2 3 2 NI~ 1 

for turbulent flow 

analysis 

Experience gained NI NI NI 1 2 2 

for structural 

analysis 

Experience gained NI NI NI 1 NI 2 

for flow - solid 

interaction 

Computer time to 

solve the same 

turbulent flows 

Computer time to 

solve the same 

nonlinear structural 

problems 

1 

NI 

1 

NI 

2 3 NI 1 

NI 1 2 2 

1. CV-FDM - Control-Volume Based Finite-Difference Method 

2. BFG-FDM - Body-Fitted Grid Finite Difference Method 

3. CV-FEM - Control-Volume Based Finite Element Method 

4. FEM - Finite Element Method 

5. BEM - Boundary Element Method 

6. FDM - Finite Difference Model 

7. NI - No Informsiion Available. 

55 
2) One of the best performing elements for shell structures are the Ahmad 

element [ 551 and the semiloop element 11261 . Use of these elements for geometrically 

and /or materially nonlinear structures has been reported elsewhere , including References 55 and 126. But use of these elements for vibration problems which require 

rigorous mathematical modeling of elasto-viscoplasticity [ 931 , such as flow-solid interaction analysis of tube banks, has not appeared as yet. 

two types of shell elements need to be further tested and developed for the present 

problem. 

3) 

It is recommended that these 

Many of the structural elements including both of the Ahmad element and 

the semiloop element are based on the eight-node serendipity element. 

of the serendipity element for flow calculations is recommended. The serendipity 

element or the complete quadratic element can approximate the curved and smooth 

boundary far better than linear elements. 

Therefore , use 

4) Flow-solid interaction in a nest of thin cylinders subjected to cross-flows 

can be characterized by the large deformation of the structures which entails continuously changing flow domain. It would be advantageous, therefore, to use the 

Lagrangian formulation for the structural analysis and to use the Euler-Lagrangian 

formulation , introduced in Section 3.4 , for the turbulent flow. 

5) The body-fitted grid finite-difference method has already made significant 

progress toward flow-solid interaction analysis in tube banks. 

this method be further developed because it can be used together with the finite 

element structural analysis codes. 

would be the complexities associated with a computer program structure originating 

from combining two numerical methods which are quite different in their nature. 

It is recommended that 

The only drawback envisioned with this approach 

6) Regardless of whether the finite element or the finite difference method is 

used for flow-solid interaction analysis in tube banks, difficulties may be encountered 

as a result of the use of unsteady turbulence models which are not well established 

as yet. There exist a few turbulence models well established for steady state turbulence models as has been discussed in Section 3.3. There also appeared a multiple 

scale k- E turbulence model proposed by Hanjalic et al. [ 1201 and an extended multiple 

scale turbulence model by Kim [ 121. At this time, it is not clear which 

turbulence model would be most suitable to be extended to unsteady turbulent flows. 

It should be noted that the straight forward use of the standard k-E turbulence model 

to unsteady flows has almost always been unsuccessful. Accordingly , any turbulence 

model to be used for flow-solid interaction needs to be tested against available experimental data. It is also recommended that tests for unsteady, turbulent flows be made 

of as many turbulence models as possible. 

A successful numerical modeling of the flow-solid interaction would enhance the 

efficiency and confidence in the design of tube banks. For example, fretting of 

cylinders due to large deformation can be predicted directly , viscoplastic deformations 

which will shorten the lifetime of tube banks can be predicted directly, and any 

instability of tube banks subjected to cross-flows could be predicted using the nonlinear stability criterion analysis method [ 551 . 

56 
REFERENCES 

1. Blevins, R. D.: Flow Induced Vibrations. Van Nostrand Co., New York, 1977. 

2. Chen, Y . N . : 

for Fluid-Elastic Coupling in Cross Flows. 

Structure Interaction Phenomena in Pressure Vessel and Piping System, ASME 

Winter Annual Meeting, Atlanta, Georgia, 1977. 

The Sensitive Tube Spacing Region of Tube Bank Heat Exchangers 

Presented at the Symposium on Fluid 

3. Roberts, B . W. : Low Frequency Aeroelastic Vibrations in Cascade of Circular 

Cylinders. Mechanical Engineering Science Monograph No. 4, September 1966. 

4. Connors, H. J. : Fluidelastic Vibration of Tube Arrays Excited by Cross Flow. 

Presented at the Symposium on Flow Induced Vibration in Heat Exchangers, 

ASME Winter Annual Meeting, December 1970. 

5. Paidoussis, M. P. : Fluidelastic Vibration of Cylinder Arrays in Axial and Cross 

Flow. Proceedings of the Joint Pressure Vessels and Piping, Materials, Nuclear 

Engineering and Solar Conference, pp . 11- 49, Denver, Colorado, June 1981. 

6. Chen, Y. N.: 

Due to von Karman Streets. ASME J. Engineering for Industry, Vol. 90, pp. 

Flow Induced Vibration and Noise in Tube Bank Heat Exchanges 

134- 146, 1968. 

7. Owen, P. R. : Buffeting Excitation of Boiler Tube Vibration. J. Mechanical 

Engineering Science, Vol. 7, pp. 431-439, 1965. 

8. Paidoussis, M. P. : Flow Induced Vibrations in Nuclear Reactors and Heat 

Exchangers : Practical Experiences and State of Knowledge. Practical Experiences with Flow Induced Vibrations, pp. 1- 81, Springer-Verlag , Berlin, 1980. 

9. Crandal, S. H., and Mark, W. D. : Random Vibration in Mechanical Systems. 

Academic Press, New York, 1963. 

10. Richard, E. J., and Mead, D. J. : Noise and Acoustic Fatigue in Aeronautics. 

J. Wiley and Sons, New Yorl, 1968. 

11. Wambsganss, M. W., Jr., and Chen, S. S. : Tentative Design Guide for Calculating the Vibration Response of Flexible Cylindrical Elements in Axial Flow. 

Argonne National Laboratory Report ANL-ETD- 71-07, 1971. 

12. Chen, S. S., and Wambsganss, M. W., Jr. : Parallel Flow Induced Vibration of 

Fuel Rods. Nucl. Engr. Design, Vol. 18, pp. 253-278, 1972. 

13. Ash, R. L. : Simulation of Turbulent Wall Pressure. NASA CR-2958, 1978. 

14. Ash, R. L. : Simulation of Turbulent Wall Pressure. NASA CR- 146534, 1976. 

15. Wlezien, R. W., and Way, J. L. : Techniques for the Experimental Investigation 

of the Near Wake of a Circular Cylinder. J. AIAA, Vol. 17, No. 6, pp. 563- 

570, 1979. 

16. Buckingham, A. C., Hall, M. S., and Chun, R. C.: Numerical Simulation of 

Compliant Material Response to Turbulent Flow. 

January 1984. 

AIAA Paper AIAA-84-0537, 

57 
17. 

18. 

19. 

20. 

21. 

22. 

23. 

24. 

25. 

26. 

27. 

28. 

29. 

30. 

58 

Buckingham, A. C., Hall, M. S., and Chun, R. C.: Numerical Simulation of 

Compliant Material Response to Turbulent Flow. J. AIAA Vol. 23, No. 7, pp. 

1046- 1052, 1985. 

Blevins, R. D. : Fluidelastic Whirling of a Tube Row. J. Pressure Vessel 

Technology, pp. 263-267, 1974. 

Heinecke, E. P. : 

Cross Flow. 

Nuclear Plant, Keswick, May 1978. 

Fluidelastic Vibrations in Heat Exchangers with Tubes in 

Proceedings of BNES International Conference on Vibration in 

Chen, S. S. : Guidelines for the Instability Flow Velocity of Tube Arrays in 

Crossflow. J. Sound and Vibration, Vol. 93, No. 3, pp. 439-445, 1984. 

Chen, Y. N.: The Orbital Movement of and the Damping of the Fluidelastic 

Vibration of Tube Banks due to Vortex Formation, Part 2, Criterion for the 

Fluidelastic Orbital Vibration of Tube Arrays. 

96, Series B, pp. 1065-1071, 1974. 

Transactions of the ASME, Vol. 

Gross, H . : 

Application in Design of Tubular Heat Exchangers. Ph .D . Thesis , Technical 

University of Hanover , 1975. 

Investigations in Aeroelastic Vibration Mechanisms and Their 

Gorman, D. J. : 

Cross Flow Induced Vibration in Nuclear Reactor Heat Exchange Equipment. 

Nuclear Science and Engineering, Vol. 61, pp. 324-336, 1976. 

Experimental Development of Design Criteria to Limit Liquid 

Savkar, S. D. : A Brief Review of Flow Induced Vibration of Tube Arrays in 

Cross Flow. J. Fluid Engineering, Vol. 99, pp. 517-519, 1977. 

Connors, H . J. : Fluidelastic Vibration of Heat Exchanger Tube Arrays. J. 

Mechanical Design, Vol. 100, pp. 347-353, 1977. 

Pettigrew, M. J. , Sylvestre, Y. , and Campagna, A. 0. : 

of Heat Exchanger and Steam Generator Designs. 

Design, Vol. 48, pp. 97-115, 1978. 

Vibration Analysis 

Nuclear Engineering and 

Weaver, D. S. and Grover, L. K.: Cross Flow Induced Vibrations in a Tube 

Bank - Turbulent Buffeting and Fluidelastic Instability. J. Sound and Vibration, Vol. 59, pp. 277-294, 1978. 

Chen, S. S., and Jendrzejczyk , J. A. : Experiments on Fluidelastic Instability 

in Tube Banks Subjected to Liquid Cross Flow. J. Sound and Vibration, Vol. 

78, pp. 355-381, 1981. 

Tanaka, H., Takahara, S. , and Ohta, K. : Flow Induced Vibration of Tube 

Arrays with Various Pitch-to-Daimeter Ratios. ASME, Flow Induced Vibration 

of Circular Cylindrical Structures, Eds. Chen, S. S., et al. , pp. 45-56, New 

York, 1982. 

Paidoussis, M. P. : Fluidelastic Vibration of Cylinder Arrays in Axial and Cross 

Flow, State of the Art, ASME, Flow Induced Vibration Design Guidelines, Ed. 

Chen, P. Y., pp. 11-46, New York, 1980. 
31. 

32. 

33. 

34. 

35. 

36. 

37. 

38. 

39. 

40. 

41. 

42. 

43. 

44. 

45. 

Blevins, R. D. : Discussions of Guidelines for the Instability Flow Velocity of 

Tube Arrays in Crossflow. J. Sound and Vibration, Vol. 97, No. 4, pp. 641- 

643, 1984. 

Chen, S. S.: 

643- 644, 1984. 

Chen, S. S.: 

Proceedings of 

Author's Reply. J. Sound and 

The Instability Flow Velocity of 

the International Conference on 

Vibration, Vol. 97, No. 4, pp. 

Tube Arrays in Cross Flow. 

Flow Induced Vibration in Fluid 

Engineering, Reading, England, September 1982. 

Candel, S. M. : A Review of Numerical Methods in Acoustic Wave Analysis. 

ONERA TP NO. 1983-83, 1983. 

Bernhard, R. J. : A Finite Element Method for Synthesis of Acoustical Shapes. 

J. Sound and Vibration, Vol. 98, pp. 55-65, 1985. 

Blevins, R. D. : 

J. Sound and Vibration, Vol. 92, No. 4, pp. 415-520, 1984. 

Quinn, M. C., and Howe, M. S. : The Influence of Mean Flow on the Acoustic 

Properties of a Tube Bank. 

NASA Aeroacoustics Conference, Williamsburg , Virginia, October 1984. 

Zdravkovich, M. M., and Nuttall, J. A. : On the Elimination of Aerodynamic 

Noise in a Staggered Tube Banks. J. Sound and Vibration, Vol. 34, No. 2, 

Review of Sound Induced by Vortex Shedding from Cylinders. 

Paper No. AIAA-84-2312, presented at the 9th AIAA/ 

pp. 173-177, 1974. 

Baird, R. C. : Pulsation Induced Vibration in Utility Steam Generation Units. 

Transactions of the ASME, Vol. 25, pp. 38-44, 1954. 

Cohan, L. J., and Dean, W. J. : Elimination of Destructive, Self-excited Vibration in Large, Gas and Oil-Fired Utility Units. 

Vol. 87, pp. 223-228, 1965. 

Morison, J. R. : The Force Exerted by Surface Waves on Piles. Petroleum 

Transactions, Vol. 189, pp. 149-157, 1950. 

Srinivasan, R., and Mongia, H. C. : Numerical Computation of Swirling Recirculating Flow. NASA CR- 165196, 1981. 

Bankston, C. P., and Back, L. H. : Heat Transfer from Combustion Gases to 

a Single Row of Closely Spaced Tubes in a Swirl Crossflow Stirling Engine 

Heater. J. Heat Transfer, Vol. 104, pp. 55-61, 1982. 

Singhal , A. K . , et al. : Analysis of Physical-Chemical Processes Governing 

SSME Internal Fluid Flow. Cham of North America, Report No. 4045/22, 

Alabama, June 1985. 

Chen, S. S. : A Mathematical Model for Cross Flow Induced Vibrations of Tube 

Rows. Technical Memorandum, ANL-CT- 77-4, Argonne National Laboratory, 

Illinois, 1976. 

Chemical Process and Engineering, 

59 
46. 

47. 

48. 

49. 

50. 

51. 

52. 

53. 

54. 

55. 

56. 

57. 

58. 

59. 

60 

Price, S. J., and Paidoussis, M. P. : An Improved Mathematical Model for the 

Stability of Cylinder Rows Subjected to Cross Flow. J. Sound and Vibration, 

Vol. 97, NO. 4, pp. 615-640, 1984. 

Price, S. J., and Paidoussis, M. P. : Fluidelastic Instability of a Double Row 

of Circular Cylinders Subjected to a Cross Flow. 

Acoustics, Stress and Reliability in Design, Vol. 105, pp. 59-66, 1983. 

ASME J. of Design, Vibration, 

Simpson, A,, and Flower, J. W. : An Improved Mathematical Model for the 

Aerodynamic Forces on Tandom Cylinders in Motion with Aeroelastic Applications. 

J. Sound and Vibration, Vol. 51, pp. 183-217, 1977. 

Noh, W. F. : A Time-Dependent , Two- Space Dimensional, Coupled EulerianLagrangian Code. Methods in Computational Physics, B. Alder, et al., Eds. , Academic Press, New York, 1964. 

Hirt, C. W., Amsden, A. A., and Cook, J. L.: An Arbitrary LagrangianEulerian Computing Method for All Flow Speeds. J. Computational Physics, 

Vol. 14, pp. 227-253, 1974. 

Donea, J., Giuliani, S., and Halleux, J. P. : An Arbitrary Lagrangian-Eulerian 

Finite Element Method for Transient Dynamic Fluid- Structure Interactions. 

Comput. Meths. Appl. Mech. Engrg., Vol. 33, pp. 689-723, 1982. 

Belytschko, T., Flangan, D. P., and Kennedy, J. M. : Finite Element Methods 

with User-Controlled Meshes for Fluid-Structure Interactions. Comput. Meths. 

Appl. Mech. Engrg., Vol. 33, pp. 669-688, 1982. 

Belytschko, T . , Kennedy, J. M. , and Schoeberle , D . F. : Quasi-Eulerian 

Finite Element Formulation for Fluid- Structure Interactions. Transactions of the 

ASME, J. Pressure Vessel Technology, Vol. 102, pp. 62-69, 1980. 

Liu, W. K. : Development of Finite Element Procedures and Computer Implementation Aspects in Fluid-Structure Interactions. New Concepts in Finite Element 

Analysis. Eds. Hughes, et al., pp. 247-268, ASME, New York, 1981. 

Zienkiewicz, 0. C . : The Finite Element Method. 3rd Edition, McGraw-Hill, 

New York, 1977. 

Halleux, J. P. : Modeling of Elasto-Plastic Material Behaviour. Report EUR7329, Commissions of the European Communities , 1981. 

Belytschko, T., and Hsieh, B . J. : Nonlinear Transient Finite Element Analysis 

with Convected Coordinates. 

1973. 

Int. J. Num. Meth. Engr., Vol. 7, pp. 255-271, 

Zienkiewicz, 0. C., Loehner, R., Morgan, K., and Nakazawa, S.: Finite 

Elements in Fluid Mechanics - A Decade of Progress. Finite Elements in Fluids, 

Vol. 5, Eds. R. H. Gallagher, et al., pp. 1-26, J. Wiley and Sons, New York, 

1984. 

Oden, J. T.: RIP Methods for Stokesian Flow. Finite Elements in Fluids, Vol. 

4, Eds. R. H. Gallagher, et al., pp. 305-318, J. Wiley and Sons, New York, 

1982. 
60. 

61. 

62. 

63. 

1 

64. 

65. 

66. 

67. 

68. 

69. 

70. 

71. 

72. 

74. 

75. 

Hallquist , J. 0. : User's Manual for DYNA-2D, An Explicit Two-Dimensional 

Hydrodynamics Finite Element Code with Interactive Rezoning. Report No. UCID18756, Lawrence Livermore National ILaboratory , 1982. 

Hallquist , J. 0. : User's Manual for DYNA3D and DYNAP (Nonlinear Dynamic 

Analysis of Solids in Three Dimensions). Report No. UCID-19156, Lawrence 

Livermore National Laboratory, 1981. 

Hallquist, J. 0. : Theoretical Manyal for DYNA3D. Report No. UCID-19401, 

Lawrence Livermore National Laboratory, 1983. 

Orszag, S. A., and Kells, L. C. : Transition to Turbulence in Plane Poiseuille 

and Couette Flow. J. Fluid Mechanics, Vol. 96, pp. 159-205, 1980. 

Orszag, S. A., and Patera, A. T. : Subcritical Transition to Turbulence in 

Plane Channel Flows. Physical Review Letters, Vol. 45, No. 12, pp. 989-993, 

1980. 

Cooley, J. W., and Tukey, J. W.: An Algorithm for the Machine Computation 

of Complex Fourier Series. Mathematics of Computation, Vol. 19, pp . 297- 301, 

1965. 

Gottlieb, D., and Orszag, S. A. : Numerical Analysis of Spectral Method: 

Theory and Application. NSF-CMBS Monograph 26, SIAM, 1977. 

Orszag, S. A. : Numerical Simulation of Incompressible Flow within Simple 

Boundaries, I. Galerkin (Spectral) Representations. Studies in Applied Mathematics, Vol. L., No. 4, 1971. 

Orszag, S . A. : 

Studies in Applied Mathematics, Vol. LI, No. 3, 1972. 

Comparison of Pseudospectral and Spectral Approximation. 

Nallasamy , M. : 

to Internal Fluid Flow. NASA TP-2474, 1985. 

A Critical Evaluation of Various Turbulence Models as Applied 

Marvin, J. G . : Turbulence Modeling for Computational Aerodynamics. J . AIAA , Vol. 21, NO. 7, pp. 941-955, 1983. 

Lakshminarayana, B . : Paper No. 

AIAA- 85- 1652, AIAA 18th Fluid Dynamics, Plasmadynamics, and Laser Conference, 

Cincinatti, Ohio, July 1985. 

Turbulence Modeling for Complex Flows. 

Launder, B. E., and Yang, W. M. : Secondary Flows in Ducts of Square Cross 

Section. J. Fluid Mechanics, Vol. 54, pp. 289-295, 1972. 

Hanjalic, K., and Launder, B. E.: A Reynolds Stress Model of Turbulence and 

Its Application to Thin Shear Flows. J. Fluid Mechanics, Vol. 52, pp. 609-638, 

1972. 

Rodi, W. : A New Algebraic Relation for Calculating the Reynolds Stresses. 

EAMM, Vol. 56, pp. 219-221, 1976. 

Launder, B. E., Reece, G. J., and Rodi, W. : Progress in the Development of 

a Reynolds Stress Turbulence Closure. J. Fluid Mechanics, Vol. 68, pp. 537- 

566, 1975. 

61 
76. 

77. 

78. 

79. 

80. 

81. 

82. 

83. 

84. 

85. 

86. 

87. 

88. 

89. 

90. 

62 

Rodi, W. : Examples of Turbulence Models for Incompressible Flows. J . AIAA , Vol. 20, NO. 7, pp. 872-879, 1982. 

Hanjalic, K., and Launder, B. E. : 

Closure for Low Reynolds Number Turbulence. J. Fluid Engineering, Vol. 74, 

pp. 593-610, 1976. 

Contribution Towards a Reynolds Stress 

Launder, B . E., and Spalding, D. B . : The Numerical Computation of Turbulent Flow. Comp. Meth. Appl. Mech. Engrg., Vol. 3, pp. 269-289, 1974. 

Patel, V. C., Rodi, W., and Scheuerer, G. : Turbulence Models for Near-Wall 

and Low Reynolds Number Flows: A Review. J. AIAA, Vol. 23, No. 9, pp. 

1308- 1319, 1985. 1 

Jones, W. P., and Launder, B. E. : The Prediction of Laminarization with a 

Two-Equation Model of Turbulence. Int. J. Heat Mass Transfer, Vol. 5, pp. 

301-314, 1973. 

Naot, D., and Rodi, W. : Numerical Simulation of Secondary Currents in Open 

Channel Flow with an Algebraic Stress Turbulence Model. ASCE J. Hydraulic 

Div., Vol. 108, pp. 948-969, 1982. 

Gibson, M. M. : Prediction of Curved Free Shear Layers with a Reynolds Stress 

Model of Turbulence. Proc. the 2nd Symposium on Turbulent Shear Flows, 

London, 1979. 

Celik, I., and Rodi, W. : Calculation of Oscillatory Turbulent Flows in Open 

Channels. 

Louisiana, February, 1984. 

Proc . Energy Sources Technology Conference, New Orleans, 

Telionis, D. P. : Unsteady Viscous Flows. Springer-Verlag , New York, 1981. 

Boyle, D. K., and Golay, M. W. : Measurement of a Recirculating TwoDimensional Turbulent Flow and Comparison to Turbulence Model Predictions. 

11. Transient Case. ASME Winter Annual Meeting, Boston, MA., November 1983. 

Schumann, U . , Grotzback, G . , and Kleiser , L. : Direct Numerical Simulation of 

Turbulence. Prediction Methods for Turbulent Flows, Ed. W. Kollmann , Hemisphere Publishing Co., New York, 1980. 

Leonard, A. : 

Flows. Advances in Geophysics, Vol. 18A, pp. 237-248, 1974. 

Energy Cascade in Large Eddy Simulations of Turbulent Fluid 

Deardroff, J. W. : A Numerical Study of Three-Dimensional Turbulent Channel 

Flow at Large Reynolds Number. 

1970. 

J. -Fluid Mechanics, Vol. 41, pp . 453- 480, 

Lilly , D. K . : 

Sciences, IBM Form No. 320-1951, 1967. 

Proceedings IBM Scientific Computing Symposium on Environmental 

Patankar, S. V. : Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Co., 1980. 
91. Schumann, U. : Subgrid Scale Model for Finite Difference Simulations of Turbulent Flows in Plane Channels and Annuli. J. Computational Physics, Vol. 18, 

pp. 376-404, 1975. 

92. Moin, P., and Kim, J. : Numerical Investigation of Turbulent Channel Flow. 

J. Fluid Mechanics, Vol. 118, pp. 341-377, 1982. 

93. Hughes, T. J. R., and Taylor, R. L. : Unconditionally Stable Algorithms for 

Quasi Static Elasto-Viscoplastic Finite Element Analysis. Computers and Struc- 

- tures, Vol. 8, pp. 169-173, 1978. 

94. Owen, D. R. J., and Hinton, E. : Finite Elements in Plasticity. Pineridge 

t Press Ltd., Swansea, U. K., 1980. 

95. Perzyna , E. : Fundamental Problems in Viscoplasticity . Recent Advances in 

Applied Mechanics, Academic Press, New York, 1966. 

96. Dawson, P. R. : Development of Models for Simulating Warm and Hot Forming 

Processes. Proc. Manufacturing Systems Research Conference, SME , 1985. 

97. Anand, L. : Constitutive Equations for the Rate Dependent Deformaiton of 

Metals at Elevated Temperatures. ASME, J. Engineering Materials and Technology, Vol. 104, pp. 12-17, 1982. 

98. Bodner , S. R. , and Partom, Y . : Constitutive Equations for Elastic-Viscoplastic 

Strain Hardening Materials. ASME, J. Applied Mechanics, Vol. 42, pp. 385- 

389, 1975. 

99. Rodal, J. J. A., Steigmann, D. J., and Witmer, E. A.: Numerical Simulation 

of Transient Finite Deformations of Thin Beams. J. Applied Mechanics, Vol. 

50, pp. 765-769, 1983. 

100. Rodal, J. J. A., and Witmer, E. A. : Finite Strain Large Deflection ElasticViscoplastic Finite Element Transient Response Analysis of Structures. 

CR-159874, 1979. 

NASA 

101. Bathe, K. J., and Wilson, E. J. : Numerical Methods in Finite Element Analysis. 

Prentice-Hall, New Jersey, 1976. 

102. Schreyer, H. L., Kulak, R. F., and Kramer, J. M.: Accurate Numerical Solutions for Elastic-Plastic Models. 

101, NO. 3, pp. 226-234, 1979. 

ASME, J . Pressure Vessel Technology, Vol. 

103. Wambsganss, M. W., Yang, C. I., and Halle, H.: Fluidelastic Instability in 

Shell and Tube Heat Exchangers - A Framework for a Prediction Method. 

ASME Winter Annual Meeting, New Orleans, December 1984. 

Proc. 

104. Sha, W. T., Yang, C. I., Kao, T. T., and Cho, S. M. : Multidimensional 

Numerical Modeling of Heat Exchangers. ASME, J. Heat Transfer, Vol. 104, 

pp. 417-425, 1982. 

105. Zienkiewicz, 0. C., Loehner, R., Morgan, K., and Nakazawa, S.: Finite 

Elements in Fluid Mechanics - A Decade of Progress. 

Vol. 5, J. Wiley and Sons, New York, 1984. 

Finite Elements in Fluids, 

63 
106. 

107. 

108. 

109. 

110. 

111. 

112. 

113. 

114. 

115. 

116. 

117. 

118. 

119. 

64 

Kwak, D., and Chang, J. L. C. : 

Stokes Flow Solver, Part I. - INS3D Code. CFD Workshop, University of 

Tennessee Space Institute, Tullahoma, Tennesses, June 1985. 

A Three Dimensional Incompressible NavierChang, J. L. C., and Kwak, D. : 

Stokes Solver, Part 11. 

tion. CFD Workshop, University of Tennessee Space Institute, Tullahoma, 

Tennessee, June 1985. 

A Three Dimensional Incompressible Navier- - Space Shuttle Main Engine (SSME) Power Head SimulaRogers, S. E., Kaul, U. K., and Kwak, D.: 

Multiple LOX Posts and Its Application to the Space Shuttle Main Engine (SSME). 

Proceedings AIAA Aerospace Science Meeting, Reno, Nevada, January 1986. 

A Numerical Study of Single and 

Brooks, A. N., and Hughes, T. J. R. : Streamline Upwind/Petrov-Galerkin 

Formulations for Convection Dominated Flows with Particular Emphasis on the 

Incompressible Navier- Stokes Equation. 

1 

Comput. Meth . Applied Mech. Engrg . , Vol. 32, pp. 199-259. 

Zienkiewicz, 0. C. : 

cesses. 

J. Wiley and Sons, New York, 1984. 

Flow Formulation for Numerical Solution of Forming ProNumerical Analysis of Forming Processes, Eds. J. F. T. Pittman, et al., 

Taylor, C., and Hughes, T. G. : 

Stokes Equation. Pineridge Press, Swansea, U . K . , 1980. 

Finite Element Programming of the NavierDonea, J. : 

Int. J. Nume. Meth. Engrg., Vol. 20, pp. 101-119, 1984. 

Lohner, R., Morgan, K., and Zienkiewicz, 0. C. : 

Hyperbolic Equation Systems by the Finite Element Method. 

Meth. Engrg., Vol. 4, pp. 1043-1063, 1984. 

A Taylor-Galerkin Method for Convective Transport Problems. 

The Solution of Nonlinear 

Int. J. Nume. 

Donea, J., Giuliani, S., Laval, H., and Quartapelle, L.: Time-Accurate Solution of Advection-Diffusion Problems by Finite Elements. Comp . Meth. Applied 

Mech. Engrg., Vol. 45, pp. 123-145, 1984. 

Yu, C. C., and Heinrich, J. C. : Petrov-Galerkin Methods for the Time 

Dependent Convective Transport Equation. 

Engrg., 1986. 

Balija, B. R., and Patankar, S. V. : 

Convection-Diffusion Problems. Numerical Heat Transfer, Vol. 3, pp . 393-409, 

1980. 

will appear in Int . J . Nume. Meth. 

A New Finite Element Formulation on 

Prakash, C., and Patankar, S. V. : A Control Volume Based Finite Element 

Method for Solving the Navier-Stokes Equations Using Equal-Order VelocityPressure Interpolation. Numerical Heat Transfer, Vol. 8, pp . 259- 280 , 1985. 

Prakash, C. : 

Mass Transfer and for Fluid Flow Using Equal-Order Velocity-Pressure Interpolation. Numerical Heat Transfer, in press, 1986. 

Kikuchi, N., Oden, J. T., and Song, Y. J.: 

Methods and Smoothing Schemes of Pressure for Stokes' Flow Problems. 

Elements in Fluids, Vol. 5, Eds. R. H. Gallagher, et al., pp. 107-126, J. 

Wiley and Sons, New York, 1984. 

An Improved Control Volume Finite Element Method for Heat and 

Convergence of Modified Penalty 

Finite 
120. Hanjalic, K., Launder, B . E., and Schiestel, R. : Multiple-Time-Scale Concepts 

in Turbulent Transport Modeling. In Turbulent Shear Flows, Vol. 2, Eds. L. 

J. S. Bradbury, et al., Springer-Verlag, pp. 36-49, 1980. 

121. Kim, S.-W. : Finite Element Computation of Turbulent Flows and an 

Extended Multiple-Scale Turbulence Closure Model. In preparation. 

122. Taylor, C., Thomas, C. E., and Morgan, K.: F. E. M. and the Two Equation 

Model of Turbulence. Proc. the 2nd Int. Conf. on Numerical Methods in Laminar 

and Turbulent Flow , Eds. C . Taylor and B . A. Schrefler, Venice, Pineridge 

Press, 1981. 

123. Smith, R. PI.?.: On the Finite Element Calculation of Turbulent Flow Using the s 

k-E Model. Int. J. Nume. Fluids, Vol. 4, pp. 303-319, 1984. 

124. Benim, A. C., and Zinser, W. : Investigation into the Finite Element Analysis 

of Confined Turbulent Flows Using a k- E Model of Turbulence. 

Appl. Mech. Engrg., Vol. 51, pp. 507-523, 1985. 

Comput. Meth. 

125. Atkins, D. J., Maskell, S. J., and Patrick, M. A.: Numerical Prediction of 

Separated Flows. Int. J. Nume. Meth. Engrg., Vol. 15, No. 1, pp. 129-144, 

1980. 

126. Irons, B., and Ahmad, S. : Techniques of Finite Elements. J. Wiley and Sons, 

New Yorl, 1980. 

127. Brebbia, C. A., and Walker, S. : Boundary Element Techniques in Engineering, 

Newnes-Butterworths , London, 1980. 

65 
APPROVAL 

A CRITICAL EVALUATION OF VARIOUS METHODS FOR THE 

ANALYSIS OF FLOW-SOLID INTERACTION IN A NEST OF 

THIN CYLINDERS SUBJECTED TO CROSS FLOWS 

By Sang-Wook Kim 

The information in this report has been reviewed for 

technical content. Review of any information concerning Department 

of Defense or nuclear energy activities or programs has been made 

by the MSFC Security Classification Officer. This report, in its 

entirety, has been determined to be unclassified. 

G. F. McDonough 

Director, Systems Dynamics Laboratory 

. 

U.S. GOVERNMENT PRINTING OFFICE 1987-730-067/40078 
1. REPORT NO. 2. GOVERNMHT ACCESSION NO. 

NASA CR - 178996 

3. RECIPIENT'S CATALOG NO. 

4. TITLE AN0 SUBTITLE 

A Critical Evaluation of Various Methods for the Analysis 

Subjected to Cross Flows 

Sang-Wook Kim * 

of Flow-Solid Interaction in a Nest of Thin Cylinders 

7. AUTHOR(S) . 

9. PERFORMING ORGANIZATION NAME AND ADDRESS 

The report presents a review of various experimental, analytical, and numerical 

analysis msthods for flow-solid interaction of a nest of cylinders subjected to crossflows. 

5. REPORT DATE 

January 1987 

6. PERFORMING ORGANIZATION CODE 

B.PERFORMING ORGANIZATION REPORr A! 

(0. WORK UNIT NO. 

A nest of cylinders subjected to cross-flows can be found in numerous 

engineering applications including the Space Shuttle Main Engine - Main Injector 

Assembly (SSME-MIA) and nuclear reactor heat exchangers. Despite its extreme 

importance in engineering applications, understanding of the flow-solid interaction 

process is quite limited and design of the tube banks are mostly dependent on 

experiments and /or experimental correlation equations. 

George C. Marshall Space Flight Center 

Marshall Space Flight Center , Alabama 35812 

12. SPONSORING AGENCY NAME AND ADDRESS 

National Aeronautics and Space Administration 

Washington, D .C . 20546 

For future development of major numerical analysis methods for the flow-solid 

interaction of a nest of cylinders subjected to cross flow, various turbulence models, 

numerical analysis methods for flows, the finite element computational method of 

nonlinear structural dynamics, and existing laminar flow - solid interaction analysis 

methods are included in this report. 

11. CONTRACT OR GRANT NO. 

NASS- 35918 

13. TYPE OF REPOR;' & PERIOD COVERE[ 

Contractor Rep or t 

1.; SPONSORING AGENCY CODE 

17. KEY WORDS 

Flow- Solid Interaction, Finite Element 

Method, Finite Difference Method, 

Turbulent Flow, Elasto-Viscoplasticity , Cross-Flow , Tube Bank 

DWPC - Form 3292 (*Y 1969) 

le. D I STRl B UTlON STATEMENT 

Unclassified - Unlimited 

For sale by National Tecwcal Information Service, Springfield. Virginia 22161 

19. SECURITY CLASSIF. (d thl.'nesr) 20. SECURITY CLASSIF. (of thh page) 21. NO. OF PAGES 

Unclassified Unclassified 74 

22. PRICE 

NTIS