**4. Results and discussion**

The results presented in this section are based on unsteady CFD simulations conducted on a specific tube array geometry over a range of reduced velocities, defined here as *Ur=Ugap/fd* where *f* is the frequency of oscillation of tube 1 and *d* the tube diameter. A completed unsteady simulation provided tube forces as a function of time, which, following a FFT of the data, provided the force coefficients and associated phase information. The choice of time step is very important to the accuracy of the calculations, but also for efficient use of computational resources. The optimal simulation time step was chosen based on previous investigations (Williams, 2004; Omar, 2010; Hassan et al, 2010), which indicated that fifty or more time steps per tube oscillation period is needed. In the present computations on average the time step was based on subdividing the tube period by ~90. The total number of tube oscillations simulated, beginning from a quasi-steady initial solution, was 50 and ensured steady statistics could be obtained when applying the FFT process to the force data. The coefficient data as a function of *Ur* was then processed to obtain an FEI stability map as a function of mass damping parameter, defined here as *mdp=m2 /d2*where *m* is the mass per unit tube length and *2* is the log decrement. In generating the stability maps *mdp* was varied from 0.1 to 200.

Numerical Simulations of Unsteady

compared against the experimental data.

causes a high rate of change in the fluid forces.

coefficient. For example, the phase angle

array in the lift direction.

*C* <sup>1</sup>*XX* ; (b) *C* <sup>1</sup>*YY* (Omar, 2010).

of 

especially for

Fluid Forces in Heat Exchanger Tube Bundles 325

and Takahara (Tanaka & Takahara, 1981). Figure (3) shows the predicted fluid force

As shown in Figure (3), the general behaviour of the fluid force coefficients obtained experimentally and predicted numerically, is similar. All force coefficients decrease sharply from their highest value to a low value in the range of reduced flow velocities between 1.5 and 20, and then become mostly independent of reduced flow velocity. This is because at a low reduced flow velocity, the flow velocity and the tube velocity are comparable. Therefore, the influence of the structural motion on the flow is significant, which in turn

At reduced flow velocities greater than 20, the flow velocity is much higher than the tube velocity; therefore, the influence of the tube motion on the flow is very small. This causes the

The predicted phase angles are shown in Figure (4). The trend of the phase angle as a function of the reduced flow velocity can be divided into two groups. In the first group, as depicted by Figures (4a), the phase angle has a trend opposite to that of the fluid force

approaches a constant value of approximately 150o in the range of reduced flow velocities between 25 and 100. Phase angles with values in the range of 0 to 180o will cause the fluid damping coefficient to have a positive value, which contributes to the instability of the tube

positive value over a range of the reduced flow velocities between 1.5 and 5.

*X X*<sup>1</sup> , *Y Y*<sup>1</sup> )

> *X X*1

*X X*1 increases sharply from a negative value to a

coefficients (*CX X*1 and *CY Y*<sup>1</sup> ) and Figure (4) shows their phase angles (

fluid forces acting on Tube 1 to be almost independent of the flow velocity.

(a) (b)

Fig. 3. Fluid force coefficients for an in-line square tube array with a P/d ratio of 1.33: (a)

*Y Y*<sup>1</sup> , especially in the range of reduced flow velocities between 10 and 80.

An example of the second group is presented by Figure (4b), where the phase angle

*Y Y*<sup>1</sup> decreases to a negative value in the lower ranges of reduced flow velocity, then gradually increases to a value of −20o at a reduced flow velocity of 100. In general, the agreement between the predicted phase angle and the reported experimental results is good

*X X*<sup>1</sup> . The largest deviation from the experimental results is found in the case

An in-line square tube array were considered where the geometry, as given in Figure 1, *P/d*=1.33 *W*=0.288 m. A set of force coefficients and phase shift versus *Ur* curves was created and a stability map generated. For a subset of this data, in particular for the *P/d*=1.33 case, direct comparison to the force coefficients and phase data of Tanaka and Takahara (Tanaka & Takahara, 1981) could be made as will be discussed subsequently.

### **4.1 Mesh sensitivity**

Before performing the complete set of unsteady simulations a mesh sensitivity study was undertaken to determine the best mesh resolution. This assessment was done using a range of meshes ranging in node count from 18,000 to 220,000. In conducting the simulations changes in predicted *CX1X* and *CY1Y* coefficients were monitored as a function of changing mesh resolution. From this study it was determined that 80,000 nodes was a good compromise between accuracy and computational time. In Figure 2a is shown the mesh topology used in the simulations, while in Figure 2b the near wall mesh resolution is highlighted.

Fig. 2. Mesh (80,000 node case) with a) view of tubes 1,2,3,4 and 5 intervening mesh topology (Omar, 2010; Hassan et al, 2010), and b) view of near wall mesh.

#### **4.2 Validation of the predicted fluid force coefficients**

This study presents samples of the predicted fluid force coefficients (*CX X*1 and *CY Y*<sup>1</sup> ) and their phase angles (*X X*<sup>1</sup> and *Y Y*<sup>1</sup> ) and compared against the experimental data of Tanaka

An in-line square tube array were considered where the geometry, as given in Figure 1, *P/d*=1.33 *W*=0.288 m. A set of force coefficients and phase shift versus *Ur* curves was created and a stability map generated. For a subset of this data, in particular for the *P/d*=1.33 case, direct comparison to the force coefficients and phase data of Tanaka and Takahara (Tanaka

Before performing the complete set of unsteady simulations a mesh sensitivity study was undertaken to determine the best mesh resolution. This assessment was done using a range of meshes ranging in node count from 18,000 to 220,000. In conducting the simulations changes in predicted *CX1X* and *CY1Y* coefficients were monitored as a function of changing mesh resolution. From this study it was determined that 80,000 nodes was a good compromise between accuracy and computational time. In Figure 2a is shown the mesh topology used in the simulations, while in Figure 2b the near wall mesh resolution is

(a)

(b)

This study presents samples of the predicted fluid force coefficients (*CX X*1 and *CY Y*<sup>1</sup> ) and

*Y Y*<sup>1</sup> ) and compared against the experimental data of Tanaka

Fig. 2. Mesh (80,000 node case) with a) view of tubes 1,2,3,4 and 5 intervening mesh

topology (Omar, 2010; Hassan et al, 2010), and b) view of near wall mesh.

**4.2 Validation of the predicted fluid force coefficients** 

*X X*<sup>1</sup> and

their phase angles (

& Takahara, 1981) could be made as will be discussed subsequently.

**4.1 Mesh sensitivity** 

highlighted.

and Takahara (Tanaka & Takahara, 1981). Figure (3) shows the predicted fluid force coefficients (*CX X*1 and *CY Y*<sup>1</sup> ) and Figure (4) shows their phase angles (*X X*<sup>1</sup> , *Y Y*<sup>1</sup> ) compared against the experimental data.

As shown in Figure (3), the general behaviour of the fluid force coefficients obtained experimentally and predicted numerically, is similar. All force coefficients decrease sharply from their highest value to a low value in the range of reduced flow velocities between 1.5 and 20, and then become mostly independent of reduced flow velocity. This is because at a low reduced flow velocity, the flow velocity and the tube velocity are comparable. Therefore, the influence of the structural motion on the flow is significant, which in turn causes a high rate of change in the fluid forces.

At reduced flow velocities greater than 20, the flow velocity is much higher than the tube velocity; therefore, the influence of the tube motion on the flow is very small. This causes the fluid forces acting on Tube 1 to be almost independent of the flow velocity.

The predicted phase angles are shown in Figure (4). The trend of the phase angle as a function of the reduced flow velocity can be divided into two groups. In the first group, as depicted by Figures (4a), the phase angle has a trend opposite to that of the fluid force coefficient. For example, the phase angle *X X*1 increases sharply from a negative value to a positive value over a range of the reduced flow velocities between 1.5 and 5. *X X*1 approaches a constant value of approximately 150o in the range of reduced flow velocities between 25 and 100. Phase angles with values in the range of 0 to 180o will cause the fluid damping coefficient to have a positive value, which contributes to the instability of the tube array in the lift direction.

Fig. 3. Fluid force coefficients for an in-line square tube array with a P/d ratio of 1.33: (a) *C* <sup>1</sup>*XX* ; (b) *C* <sup>1</sup>*YY* (Omar, 2010).

An example of the second group is presented by Figure (4b), where the phase angle *Y Y*<sup>1</sup> decreases to a negative value in the lower ranges of reduced flow velocity, then gradually increases to a value of −20o at a reduced flow velocity of 100. In general, the agreement between the predicted phase angle and the reported experimental results is good especially for *X X*<sup>1</sup> . The largest deviation from the experimental results is found in the case of *Y Y*<sup>1</sup> , especially in the range of reduced flow velocities between 10 and 80.

Numerical Simulations of Unsteady

et al, 2010).

**4.4 Separation and attachment** 

Fluid Forces in Heat Exchanger Tube Bundles 327

where *Ms* , *Cs* , and *Ks* are the structural mass, damping, and stiffness, respectively. *Ma* , *Ca* , and *Ka* are the fluid added mass, damping, and stiffness, respectively. Each individual

tube was assumed to have two degrees of freedoms. The size of the system matrices are the same as the total number of degrees of freedoms (18 x 18). The flow added matrices ( *Ma* , *Ca* , and *Ka* ) are functions of added fluid-damping, fluid stiffness coefficients, flow density,

flow velocity and tube diameter. Secondly, the eigenvalue extraction of the above equation was employed to study the stability of the system as a function of the flow velocity. This results in a number of complex eigenvalues and eigen vectors. The reader is forwarded to (Tanaka & Takahara, 1981; Omar, 2010) for more detail. The system is unstable if the real part of the eigenvalue is positive. The critical flow velocity is determined by solving for the flow velocity at which the real part of the eigenvalue becomes zero. Results were located on the stability map and compared to experimental results available in the open literature for the same tube array geometry. The case of Tanaka and Takahara (Tanaka & Takahara, 1981) was selected since it matches the pitch to diameter ratio of the current simulation (*P/d* of 1.33). The predicted reduced critical flow velocity, *Uc*, as a function of mass damping parameter agrees well with the experimental data sets as shown in Figure 5. The agreement

Fig. 5. Comparison of the predicted critical flow velocity with experimental data of Tanaka and Takahara (Tanaka & Takahara, 1981): ◦ experiments, − simulation (Omar, 2010; Hassan

FEI models such as the unsteady analytical model by Lever and Weaver (Lever & Weaver,

respectively in Figure 1. Generally, for either case, an angle of approximately ten degrees is applied. In addition a time delay model is input to account for the delay between stream

*1* and *2*

1981) rely on input of separation and attachment angles for their model, shown as

between the simulations and the experimental results is excellent.

The rest of the fluid force coefficients and the phase angles due to tube motion in the lift and the drag directions are presented in (Omar, 2010).

Fig. 4. Fluid force phase angles for an in-line square tube array with a P/d ratio of 1.33: (a) <sup>1</sup>*XX* ; (b) <sup>1</sup>*YY* (Omar, 2010)

#### **4.3 Prediction of stability map characteristics**

In addition to comparison against fluid force coefficients obtained by experiment, the coefficients were used to simulate a fully flexible 3x3 tube array over a *mdp* range of 0.1 to 200. Firstly, the predicted force coefficients were transformed into added stiffness and added damping parameters and incorporated in the system equations of motion. The equation of motion of the tube system is given by:

$$\left[\left[M\_s + M\_a\right]\right]\left\{\ddot{\mathbf{x}}\right\} + \left[\mathbf{C}\_s + \mathbf{C}\_a\right]\left\{\ddot{\mathbf{x}}\right\} + \left[K\_s + K\_a\right]\left\{\ddot{\mathbf{x}}\right\} = \mathbf{0}$$

The rest of the fluid force coefficients and the phase angles due to tube motion in the lift and

(a)

(b) Fig. 4. Fluid force phase angles for an in-line square tube array with a P/d ratio of 1.33: (a)

In addition to comparison against fluid force coefficients obtained by experiment, the coefficients were used to simulate a fully flexible 3x3 tube array over a *mdp* range of 0.1 to 200. Firstly, the predicted force coefficients were transformed into added stiffness and added damping parameters and incorporated in the system equations of motion. The

[ ][ ][ ] 0 *M Mx CCx KKx s a sa sa*

the drag directions are presented in (Omar, 2010).

<sup>1</sup>*XX* ; (b)

<sup>1</sup>*YY* (Omar, 2010)

**4.3 Prediction of stability map characteristics** 

equation of motion of the tube system is given by:

where *Ms* , *Cs* , and *Ks* are the structural mass, damping, and stiffness, respectively. *Ma* , *Ca* , and *Ka* are the fluid added mass, damping, and stiffness, respectively. Each individual tube was assumed to have two degrees of freedoms. The size of the system matrices are the same as the total number of degrees of freedoms (18 x 18). The flow added matrices ( *Ma* , *Ca* , and *Ka* ) are functions of added fluid-damping, fluid stiffness coefficients, flow density, flow velocity and tube diameter. Secondly, the eigenvalue extraction of the above equation was employed to study the stability of the system as a function of the flow velocity. This results in a number of complex eigenvalues and eigen vectors. The reader is forwarded to (Tanaka & Takahara, 1981; Omar, 2010) for more detail. The system is unstable if the real part of the eigenvalue is positive. The critical flow velocity is determined by solving for the flow velocity at which the real part of the eigenvalue becomes zero. Results were located on the stability map and compared to experimental results available in the open literature for the same tube array geometry. The case of Tanaka and Takahara (Tanaka & Takahara, 1981) was selected since it matches the pitch to diameter ratio of the current simulation (*P/d* of 1.33). The predicted reduced critical flow velocity, *Uc*, as a function of mass damping parameter agrees well with the experimental data sets as shown in Figure 5. The agreement between the simulations and the experimental results is excellent.

Fig. 5. Comparison of the predicted critical flow velocity with experimental data of Tanaka and Takahara (Tanaka & Takahara, 1981): ◦ experiments, − simulation (Omar, 2010; Hassan et al, 2010).

#### **4.4 Separation and attachment**

FEI models such as the unsteady analytical model by Lever and Weaver (Lever & Weaver, 1981) rely on input of separation and attachment angles for their model, shown as *1* and *2* respectively in Figure 1. Generally, for either case, an angle of approximately ten degrees is applied. In addition a time delay model is input to account for the delay between stream

Numerical Simulations of Unsteady

function of pitch-to-diameter ratio and reduced velocity.

Journal of Fluids and Structures 5, 299-322.

of Sound and Vibration 77, 19-37.

Heat Transfer 11, pp. 363-390.

Applications, AIAA Journal 32, No. 8.

**5. Conclusions** 

**6. References** 

Fluid Forces in Heat Exchanger Tube Bundles 329

In this study the use of CFD to generate force coefficients in the unsteady flow models for FEI of Chen (Chen, 1991) and Tanaka and Takahara (Tanaka & Takahara, 1981) has been undertaken. The unsteady CFD simulations, with appropriate care taken for mesh and time step resolution, yield results that enable force coefficients to be well predicted. With comprehensive studies that cover a range of pitch-to-diameter ratios and reduced velocities, a series of FEI stability curves can be generated. The curves generated in this work based on CFD derived data follow the trends of that available in the literature, and reveal expected trends such as increasing critical reduced velocity with increasing *P/d* ratio. The results thus far indicate the CFD based approach for computing stability maps for family of arrays (in this case in-line) is very promising. The use of CFD data also has the added benefit of providing considerable detail on the flow behaviour in the array, and thus enables extracting information that could be used in models such as those of Weaver and Lever (Lever & Weaver, 1981). Along this vein simulations are planned using a larger array in order to isolate the region of interest around the oscillating tube from the downstream vortex shedding. This would allow a more careful investigation of the movement of the stagnation regions on the attachment side, and the degree of mass flow redirected, as a

J. Lever and D. S. Weaver (1982). A theoretical model for the fluid-elastic instability in heat

S. J. Price and M. P.Paidoussis (1984). An improved mathematical model for the stability of cylinder rows subject to cross-flow, Journal of Sound and Vibration 97, 615-640. S.S. Chen (1991). Fluidelastic instabilities in tube bundles exposed to no uniform cross-flow,

H. Tanaka and S. Takahara (1981). Fluid elastic vibration of tube array in cross flow, Journal

D. Weber, R. Brewster, S. Chen, and T. Wei, (2001). Numerical method for fluidelastic

H. Omar, M. Hassan, & A. Gerber (2009). Numerical Estimation of Fluidelastic Instability in Staggered Tube Arrays; ASME, PVP conference 2009 -77472, Prague, Czech. H. Omar (2010). Numerical Simulation of Fluidelastic Instability in Tube Bundles, PhD's

M. Hassan, A. Gerber, H. Omar (2010). Numerical Estimation of Fluidelastic Instability in Tube Arrays, Journal of Pressure Vessel Technology, vol. 132, issue 4, 04130, 2010. G.E. Schneider and M.J. Raw (1987). Control volume finite-element method for heat transfer

F.R. Menter (1994). Two-Equation Eddy-Viscosity Turbulence Models for Engineering

Transactions, SMiRT 16, Washington CD, p. Paper number 1194.

thesis, New Brunswick, Canada, 2010, available from

http://www.synergiescanada.org/theses/unb/1089

exchanger tube bundles, ASME Journal of Pressure Vessel Technology 104, 147-158.

instability evaluation in nuclear reactor system and plant components,

and fluid flow using colocated variables – 1. Computational Procedure, Numerical

lines perturbation and tube motion. The times delay model approximates how quickly flow along the monitoring points shown in Figure 1 is affected by any intervening tube motion.

In models such as that of Lever and Weaver (Lever & Weaver, 1981), it is assumed the flow inside the tube bundle is divided into wake and channel flows. The wake flow appears as large recirculation zones located between tubes in the spaces running transverse to the flow (zones 1 and 2 in Figure 6a). Spaces between tubes aligned with the flow direction allow the flow to proceed relatively freely along lanes (or channels as depicted in Figure 6a). The outer edges of the lanes attach and separate from neighbouring tubes on the basis of *1* and *2*. Streamlines in the channel regions are assumed not to between intervening recirculation zones. Therefore the mass flow at each gap between tubes remains constant. In the present study the flow patterns are examined in light of these approximations.

Figure 6 shows the streamlines prevalent at three instances during one period (of time *T*) of tube motion in the lift direction. At time 0 s (Figure 6a) there exists a significant recirculation zone between tubes 2 and 1 (zone 1), while at time 1/3T (Figure 6b), when the tube has just passed it peak downward position, fluid streamlines are redirected from the channel 2 to the channel 1. Conversely when the tube, at 2/3 T (Figure 6c), is nearing its peak upward position, streamlines are diverted from channel 1 to channel 2. The migration of flow back and forth between channels 1 and 2 through zones 1 and 2 is not entirely periodic as other asymmetries exist in the overall flow structure.

Fig. 6. Flow visualization showing flow conditions proceeding moving tube at period intervals of 0T, 1/3T and 2/3T. (*P/d*=1.33, *Ur*=50, *T*= 3.33s).
