**2. Fluid dynamics model**

The Reynolds averaged Navier-Stokes (RANS) equations describing mass and momentum conservation are solved to obtain the time evolution of velocity and pressure in the tube array. The interactions between fluid motion and moving structures are handled by further casting the governing RANS equations in an Arbitrary Lagrangian-Eulerian (ALE). ALE accommodates moving boundaries and any subsequent deformation of the underlying discrete mesh. The mesh velocity, *umj*, is calculated based on the space conservation law:

$$\frac{d}{dt}\int\_{V}dV = \int\_{S} \mu\_{m\bar{\eta}}d\mathfrak{n}\_{\bar{\jmath}}\tag{1}$$

The left hand term represents the rate of change of a volume *V*, while the right hand side represents the integral of the surface (*S*) velocities (where *nj* is the surface normal).

#### **2.1 Conservation equations**

320 Computational Simulations and Applications

theory models such as those of Chen (Chen, 1991) and Tanaka and Takahara (Tanaka & Takahara, 1981). These approaches provide a good estimate of the critical flow velocity, with models based on the unsteady flow theory being the more general. However, the applicability of the unsteady flow theory is restricted to the availability of the force coefficients they depend on. These force coefficients have to be measured experimentally or alternatively, as is investigated in this paper, can be generated using a comprehensive CFD

Several numerical techniques and studies embarked on attempting to simulate unsteady fluid forces (Weber et al, 2001; Omar et al, 2009; Omar, 2010; Hassan et al, 2010). This, in turn, could be used in conjunction with the unsteady flow model. This approach will be presented in this chapter. There has been an attempt to numerically predict the unsteady fluid forces in a tube row by Weber et al (Weber et al, 2001). They utilized a commercial CFD package (STAR-CD) to analyze a tube row with a P/d ratio of 1.35. In this work (Weber et al, 2001), the fluid forces due to the prescribed motion of one tube were predicted. Using these fluid forces, fluid force coefficients equivalent to those of Chen (Chen, 1991) were

The results obtained via this work (Weber et al, 2001)have good agreement with other experiments. Comprehensive numerical studies to simulate the unsteady fluid forces for In line square and normal triangle tube arrays were conducted by Omar et al (Omar et al, 2009); Omar (Omar, 2010); Hassan et al (Hassan et al, 2010). In these studies (Omar et al, 2009; Omar, 2010; Hassan et al, 2010), the ability to extract fluid force coefficients for the unsteady flow theory was undertaken for two tube array configurations, namely i) in-line square tube array and ii) normal triangle tube array. The ability to predict the tube array FEI from CFD derived coefficients for each configuration and for single and fully flexible tube arrays was assessed. The predicted fluid force coefficients and FEI are then compared to available experimental data. The comparison demonstrates the viability of the proposed model to provide fluid force coefficients for the unsteady flow theory. The effect of the P/d

In the present study the accuracy of using CFD generated force coefficients in the unsteady flow theory model is tested by comparing against available experimental data for in-line square tube arrays. The base geometry, depicted in Figure 1a, and test conditions, follows

The Reynolds averaged Navier-Stokes (RANS) equations describing mass and momentum conservation are solved to obtain the time evolution of velocity and pressure in the tube array. The interactions between fluid motion and moving structures are handled by further casting the governing RANS equations in an Arbitrary Lagrangian-Eulerian (ALE). ALE accommodates moving boundaries and any subsequent deformation of the underlying discrete mesh. The mesh velocity, *umj*, is calculated based on the space conservation law:

> *V S <sup>d</sup> dV u dn*

represents the integral of the surface (*S*) velocities (where *nj* is the surface normal).

The left hand term represents the rate of change of a volume *V*, while the right hand side

*mj j*

*dt* (1)

ratio and the Reynolds number on the FEI threshold was also investigated.

the experimental study of Tanaka and Takahara (Tanaka & Takahara, 1981).

model.

obtained.

**2. Fluid dynamics model** 

The RANS equations in ALE form for mass conservation appears as:

$$\frac{\partial \hat{\rho}}{\partial t} + \frac{\partial \hat{\rho}(\mu\_j - \mu\_{mj})}{\partial \mathbf{x}\_j} = \mathbf{0} \tag{2}$$

and momentum equation:

$$\frac{\partial \rho u\_i}{\partial t} + \frac{\partial \rho u\_i (u\_j - u\_{mj})}{\partial \mathbf{x}\_j} = -\frac{\partial P}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left( \mu\_{\epsilon \mathcal{Y}} \frac{\partial u\_i}{\partial \mathbf{x}\_j} \right) + S\_i \tag{3}$$

where *uj* and *ui* are the fluid velocity components,  *is* the fluid density, *xj* and *xi* are cartesian spatial coordinates, *P* is the fluid pressure, *eff* is the effective viscosity (includes laminar and turbulent contributions) , and *Si* is any additional momentum source contributions. This form allows conservative fluid flow calculations with mesh adaptation in time. The discretization process of the governing equations is similar to that applied for finite-volume/finite-element discretization (Schneider & Raw, 1987; Omar, 2010; Hassan et al, 2010).

The *k-* based Shear Stress Transport (SST) model is used (Menter, 1994) to include the influence of turbulent mixing. The turbulence equations are cast in ALE form :

$$\frac{\partial \rho \mathbf{k}}{\partial t} + \frac{\partial \rho \mathbf{k} (\mathbf{u}\_j - \mathbf{u}\_{mj})}{\partial \mathbf{x}\_j} = \mathbf{P}\_k - \mathcal{J} \, ^\* \rho \mathbf{k} \mathbf{o} + \frac{\partial}{\partial \mathbf{x}\_j} \left| \, \Gamma\_{\mathrm{eff}} \, \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_j} \right| \tag{4a}$$

and the -equation:

$$\frac{\partial \rho \alpha o}{\partial t} + \frac{\partial \rho o (u\_j - u\_{mj})}{\partial \mathbf{x}\_j} = \alpha \frac{\alpha o}{k} P\_k - \beta \rho o \rho^2 + \frac{\partial}{\partial \mathbf{x}\_j} \left| \Gamma\_{\epsilon \mathcal{H}} \frac{\partial \alpha o}{\partial \mathbf{x}\_j} \right| \tag{4b}$$

In these equations *Pk* is the production rate of turbulence. The SST model incorporates a transformation between the *k-* model in the free stream regions of the flow, and the *k-* model in the near wall regions of the flow. Further details on the turbulence model implementation can be found in (Menter, 1994).

Based on the oscillation frequency and amplitude, the nodes on the moving tube surface (shaded tube designated as 1 in Figure 1) are moved to the new position and the remaining nodes in the interior are moved according to a set of Laplacian diffusion equations.

The simulation is first run with a static tube, to obtain a quasi steady-state solution to provide initial flow conditions for the transient simulation. Following this quasi steadystate solution a transient simulation is then started but with tube 1 in motion in the lift direction and then in the drag at a specified frequency and amplitude. Reader is forwarded to (Omar et al, 2009; Omar, 2010; Hassan et al, 2010) for more details.

For each time step convergence is achieved to within an RMS residual of 1x10-4 before moving on. The solution of the RANS equations is based on a coupled algebraic multigrid solver (Raw, 1996). Parallel computations are applied in the present case to reduce the solution time required per time step.

Numerical Simulations of Unsteady

Fluid Forces in Heat Exchanger Tube Bundles 323

1,2,3,4,5

1,2,3,4,5

*F U CpCq*

2 *y gap YjX <sup>j</sup> YjY <sup>j</sup>*

where *CX X*1 and *CY Y*1 are fluid force coefficient amplitudes. The fluid forces ( *FX X*<sup>1</sup> , *FY Y*<sup>1</sup> )

three suffixes of the coefficients represent the direction of the force, the tube index, and the direction of vibration, respectively. For example, the lift fluid force component acting on tube 1 due to the motion of tube 4 in the drag direction (Y ) is expressed as *FX Y*<sup>4</sup> . Therefore, the total lift fluid forces consists of 10 different components corresponding to the lift motion effect (*CX X*<sup>1</sup> ...*CX X*<sup>5</sup> ) and to the drag motion (*CX X*<sup>1</sup> ...*CX X*<sup>5</sup> ). Similarly, the overall drag force comprises of 10 components (*CY X*<sup>1</sup> ...*CY X*<sup>5</sup> , *CY Y*<sup>1</sup> ...*CY Y*<sup>5</sup> ). As the centre tube (1) oscillates in Y direction, *FXjY* is equal to zero except for *FX Y*4 and *FX Y*<sup>5</sup> . Similarly, as the centre tube

The reader can refer to (Tanaka & Takahara, 1981; Omar, 2010) for a description. Based on

2 *<sup>X</sup> gap X X X X X Y*

2 *<sup>y</sup> gap Y Y LX Y Y*

*F U C p C p p C q q*

*F U C q C p p C q q*

11 4 4 5 4 4 5

22 33

22 33

*YY YY*

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

*CqCq*

*XX XX*

*C pC p*

11 4 4 5 4 4 5

 */d2*

is the log decrement. In generating the stability maps *mdp*

*j*

*j*

2 *<sup>x</sup> gap XjX <sup>j</sup> XjY <sup>j</sup>*

*F U CpCq*

1

1

acting on tube j lead the displacement of tube j ( *Xj* ,*Yj* ) by phase angles (

oscillates in X direction, the *FYjX* is equal zero except of *FY X*4 and *FX Y*<sup>5</sup> .

2

2

a function of mass damping parameter, defined here as *mdp=m2*

these simplifications, Eq. 5 can be reduced to:

1

1

 

**4. Results and discussion** 

per unit tube length and *2*

was varied from 0.1 to 200.

  <sup>2</sup>

<sup>2</sup>

(5a)

(5b)

*XjX* and (6a)

(6b)

where *m* is the mass

*YjY* ). The
