Numerical Simulations of Unsteady Fluid Forces in Heat Exchanger Tube Bundles

Heat transfer is the main concern in designing shell and tube heat exchangers. Flowinduced vibration (FIV) is also a major concern while designing shell and tube heat exchangers. Fluidelastic instability (FEI), turbulence, periodic instability, and acoustic resonance are the FIV mechanisms that could cause to vibrations in heat exchanger tube bundles. FEI is the most dangers mechanism, since it can cause to tube damage in short time when the flow velocity exceeds the critical flow velocity (Ucr). Therefore, intensive researches have been undergoing in the last five decades to predict and under stand this mechanism. Different theories and models are in use to predict the FEI thresholds as function of mass damping parameters (MDP). These theories and models rely on some coefficients and parameters. Experimental approaches are used to predict these parameters for some tube arrays geometries. The experimental approach is expensive and a time consumer. Computational Fluid Dynamics (CFD) is an alternative approach proposed in this study to predict these parameters. This study utilized the CFD model to simulate the unsteady flow and the resulting fluidelastic forces in a tube bundle. Numerical simulations of in-line square tube arrays with a pitch-to-diameter (P/d) ratio of 1.33 utilizing a 2dimensional model are presented. In this model, a single tube was forced to oscillate within an otherwise rigid array. The numerical model solves the Reynolds-Average Navier-Stokes (RANS) equations for unsteady turbulent flow, and is cast in an Arbitrary LagrangianEulerian (ALE) form to handle mesh motion associated with a moving boundary. The fluidelastic instability was predicted for both single and fully flexible tube arrays over a mass damping parameter (MDP) range of 0.1 to 200. Fluid forces acting on the oscillating tube and the surrounding tubes were estimated. The predicted forces were utilized to calculate fluid force coefficients for all tubes. Fluidelastic instability is the most destructive excitation mechanism leading to rapid failure by fatigue or tube-to-tube clashing if the stability threshold is exceeded. Due to this potential for catastrophic failure intensive research has been ongoing for several decades on the topic of predicting and mitigating FEI effects. This has resulted in a vast amount of literature on the topic. Much of the research has been directed towards obtaining a reliable estimate of the critical flow velocity for the purpose of design. There are several models available to analyse FEI problems and the associated critical velocity. These models range from analytical approaches such as the models of Lever and Weaver (Lever & Weaver, 1982) and Paidoussis and Price (Paidoussis & Price, 1984) to the empirically-based unsteady flow


Introduction
Heat transfer is the main concern in designing shell and tube heat exchangers. Flowinduced vibration (FIV) is also a major concern while designing shell and tube heat exchangers. Fluidelastic instability (FEI), turbulence, periodic instability, and acoustic resonance are the FIV mechanisms that could cause to vibrations in heat exchanger tube bundles. FEI is the most dangers mechanism, since it can cause to tube damage in short time when the flow velocity exceeds the critical flow velocity (U cr ). Therefore, intensive researches have been undergoing in the last five decades to predict and under stand this mechanism. Different theories and models are in use to predict the FEI thresholds as function of mass damping parameters (MDP). These theories and models rely on some coefficients and parameters. Experimental approaches are used to predict these parameters for some tube arrays geometries. The experimental approach is expensive and a time consumer. Computational Fluid Dynamics (CFD) is an alternative approach proposed in this study to predict these parameters. This study utilized the CFD model to simulate the unsteady flow and the resulting fluidelastic forces in a tube bundle. Numerical simulations of in-line square tube arrays with a pitch-to-diameter (P/d) ratio of 1.33 utilizing a 2dimensional model are presented. In this model, a single tube was forced to oscillate within an otherwise rigid array. The numerical model solves the Reynolds-Average Navier-Stokes (RANS) equations for unsteady turbulent flow, and is cast in an Arbitrary Lagrangian-Eulerian (ALE) form to handle mesh motion associated with a moving boundary. The fluidelastic instability was predicted for both single and fully flexible tube arrays over a mass damping parameter (MDP) range of 0.1 to 200. Fluid forces acting on the oscillating tube and the surrounding tubes were estimated. The predicted forces were utilized to calculate fluid force coefficients for all tubes. Fluidelastic instability is the most destructive excitation mechanism leading to rapid failure by fatigue or tube-to-tube clashing if the stability threshold is exceeded. Due to this potential for catastrophic failure intensive research has been ongoing for several decades on the topic of predicting and mitigating FEI effects. This has resulted in a vast amount of literature on the topic. Much of the research has been directed towards obtaining a reliable estimate of the critical flow velocity for the purpose of design. There are several models available to analyse FEI problems and the associated critical velocity. These models range from analytical approaches such as the models of Lever and Weaver (Lever & Weaver, 1982) and Païdoussis and Price (Païdoussis & Price, 1984) to the empirically-based unsteady flow www.intechopen.com 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 model. 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 obtained. 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 ratio and the Reynolds number on the FEI threshold was also investigated. 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 experimental study of Tanaka and Takahara (Tanaka & Takahara, 1981).

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, u mj , is calculated based on the space conservation law: 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 n j is the surface normal).

Conservation equations
The RANS equations in ALE form for mass conservation appears as: (2) and momentum equation: where u j and u i are the fluid velocity components,  is the fluid density, x j and x i are cartesian spatial coordinates, P is the fluid pressure,  eff is the effective viscosity (includes laminar and turbulent contributions) , and S i 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 : and the -equation: In these equations P k 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.

Unsteady flow model
The unsteady fluid force includes the inertia, the damping, and the stiffness effects. The inertia effect is assumed to be independent of the flow velocity and can be measured experimentally, or computed using the potential-flow theory. In order to utilize this model, fluid force coefficients are required. The first attempt to obtain these force Coefficients experimentally was done by Tanaka and Takahara (Tanaka & Takahara, 1981). The damping and Stiffness effects were expressed in terms of fluid force amplitudes (C) and phase angles (φ) between fluid forces and tube oscillation. To validate the CFD results with experimental data, a simpler variant of the unsteady flow theory as presented by Tanaka and Takahara (Tanaka & Takahara, 1981) is used. The formulation is presented briefly here, however a complete description and derivation is available in (Tanaka & Takahara, 1981;Omar, 2010).  (Tanaka & Takahara, 1981;Omar, 2010;Hassan et al, 2010), and b) close-up of moving tube (shaded) and related dimensions.
Referring to the square array in Figure 1a, and considering the forces and displacements in x as well as y, and assuming the effects of all surrounding tubes can be summed linearly (and that only the four tubes immediately adjacent the center tube 1 have significant influence) the total force expected on tube 1 can be expressed as

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 U r =U gap /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 U r was then processed to obtain an FEI stability map as a function of mass damping parameter, defined here as mdp=m2 /d 2 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.

www.intechopen.com
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 U r 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.

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 C X1X and C Y1Y 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. and Takahara (Tanaka & Takahara, 1981). Figure (3 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 1 XX  increases sharply from a negative value to a positive value over a range of the reduced flow velocities between 1.5 and 5.  Omar, 2010).
An example of the second group is presented by Figure (4b), where the phase angle 1 YY  decreases to a negative value in the lower ranges of reduced flow velocity, then gradually increases to a value of −20 o 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 1 XX  . The largest deviation from the experimental results is found in the case of 1 YY  , especially in the range of reduced flow velocities between 10 and 80.

www.intechopen.com
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).

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: 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 ( a M , a C , and a K ) 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, U c , 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).

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 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. www.intechopen.com

Conclusions
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 function of pitch-to-diameter ratio and reduced velocity. The purpose of this book is to introduce researchers and graduate students to a broad range of applications of computational simulations, with a particular emphasis on those involving computational fluid dynamics (CFD) simulations. The book is divided into three parts: Part I covers some basic research topics and development in numerical algorithms for CFD simulations, including Reynolds stress transport modeling, central difference schemes for convection-diffusion equations, and flow simulations involving simple geometries such as a flat plate or a vertical channel. Part II covers a variety of important applications in which CFD simulations play a crucial role, including combustion process and automobile engine design, fluid heat exchange, airborne contaminant dispersion over buildings and atmospheric flow around a re-entry capsule, gas-solid two phase flow in long pipes, free surface flow around a ship hull, and hydrodynamic analysis of electrochemical cells. Part III covers applications of non-CFD based computational simulations, including atmospheric optical communications, climate system simulations, porous media flow, combustion, solidification, and sound field simulations for optimal acoustic effects.