Experimental and Numerical Investigations of Turbulent Flow in a Staggered Tube Bundle

Abstract

Experimental and numerical investigations were carried out on turbulent flow in a staggered tube bundle. In the experiment, a hot-wire anemometer was employed to investigate the developments of the flow at the Reynolds number of 9160. In microscopic numerical study, the standard k-ε model was used with a boundary fitted coordinate system to calculate velocity, pressure drop and kinetic energy of turbulence. The new macroscopic turbulence model with volume averaging and additional source terms was proposed for turbulent flows in tube bundles or packed beds. The turbulence constants were determined by the microscopic numerical experiments. The microscopic and macroscopic numerical results were in good agreement with the experimental results, or the empirical equation in the axial velocity and kinetic energy of turbulence.

1. Introduction

Tube bundles have been widely used in many industry processes such as heat exchangers, cooling systems for nuclear power plants and so on because of its good heat transfer and mixing characteristics. It is very important to study turbulent flows in tube bundles for the effective design. Turbulent flows in the tube bundles have been investigated both experimentally and numerically by many researchers. It is, however, difficult to obtain the experimental data on the turbulent flows in tube bundles and to calculate the turbulence flows due to the complex geometries of tube bundles.

Most of the earlier experimental studies on flows in the tube bundles focused on the measurement of pressure drop. Jakob¹⁾ studied the pressure drop in tube bundles by changing the transverse and longitudinal pitches in the tube bundles and the Reynolds number, and proposed the correlations for pressure drop through the tube bundles. The correlations are simple to apply and used well.

Recently, turbulence characteristics are focused on in the tube bundles. Although it is difficult to measure fluid velocities and turbulence characteristics in such a complicated and high turbulent flow field, some experimental studies have been reported. LDA (Laser Doppler Anemometry) measurements were performed by Simonin and Barcouda²⁾ and Balabani and Yianneskis.³⁾ Also, PIV (Particle Image Velocimetry) measurements were carried out by Iwaki et al.⁴⁾ and Paul et al.⁵⁾ There are a few studies of hot-wire anemometry measurements in tube bundles. Hot-wire is the most common method to measure fluid velocities, because it has small probes which offer high spatial resolution and excellent frequency response characteristics. Aiba et al.⁶⁾ employed hot-wire anemometry to measure mean fluid velocities and turbulent intensities in the tube bundles. They investigated the detailed turbulence intensities along the transverse direction, but they reported the turbulence characteristics at only two locations along the axial direction.

Numerical modeling of turbulent flows provides a useful means to predict the fluid flows and their characteristics such as the kinetic energy of turbulence. The Reynolds-Averaged Navier-Stokes (RANS) with different microscopic turbulence models has been used to predict the turbulent flows in staggered tube bundles. Zijlema et al.⁷⁾ employed the standard k-ε turbulence model^8,9) with a boundary fitted coordinate system in the periodic flows in the tube bundle. Watterson et al.¹⁰⁾ used the low Reynolds number turbulence model¹¹⁾ on unstructured grids. Wang et al.¹²⁾ used the Reynolds stress model to investigate the turbulent flows in the tube bundles with different wall treatments. Those numerical studies were carried out for periodic turbulent flows in the tube bundle with 7 rows at the maximum. The effects of the inlet kinetic energy of turbulence on the development of the turbulent flows in the tube bundles were not investigated in those numerical studies.

On the other hand, the macroscopic model is to regard the tube bundle as the packed bed and to calculate the fluid flows in the backed bed macroscopically without any complex geometry along the surface of the tube. Some macroscopic turbulence models with volume averaging and additional source terms for fluid-particle interactions were reported to calculate turbulent flows in packed beds. As compared with the microscopic turbulence models, the macroscopic turbulence models have advantages in the generation of grids and the cost of computations. Takeda and Lockwood¹³⁾ proposed the macroscopic k-l one-equation turbulence model in packed beds. They determined the turbulence model constant by setting the mixing length, l, and using experimental data in packed beds. Nakayama and Kuwahara¹⁴⁾ and Pedras and Lemos¹⁵⁾ proposed the macroscopic k-ε two-equation turbulence models. The turbulence constants in both models are determined by the numerical experiments of fully developed flows in two-dimensional square or circular tube bundles. However, the turbulence characteristics predicted by those models were not compared with experimental data for the model validation.

In the present work, the tube bundle is considered as the simplified model of the packed bed since fluid flows and their characteristics in the packed bed are very complex ones. Experimental and numerical studies are carried out to investigate the turbulent flows in a staggered tube bundle. 15 rows in a staggered tube bundle are prepared, that seems to be much enough to investigate the characteristics of the fluid flow especially in the downstream and to confirm that the kinetic energy of turbulence converge to the constant value or not, in other words, to understand the effect of the inlet kinetic energy of turbulence. In the experimental study, the development of the flow is examined at the Reynolds number, Re_∞, of 9160 with a hot-wire anemometer measurement. In the numerical study, firstly the microscopic turbulence model is used to calculate fluid velocities, pressure drops and kinetic energies of turbulence in the tube bundle, which are compared with the experimental results or the empirical equation. Secondly the macroscopic turbulence model is proposed by using the results of the microscopic turbulence model and is employed to calculate fluid velocities, pressure drops and kinetic energies of turbulence in the tube bundle. The effects of the inlet kinetic energy of turbulence on the development of the turbulent flows in the tube bundle are investigated in the numerical studies.

2. Experiment

Experiments were carried out in a wind tunnel with air at approximately 20°C as the working fluid. Figure 1 shows a schematic diagram of the wind tunnel test section with a staggered tube bundle. The dimensions of the test section are 200 mm in width, 100 mm in height and 600 mm in length. The test section consists of 15 staggered rows of acrylic tubes with an outside diameter of D = ϕ20 mm and a length of 100 mm with each row consisting of 4 or 5 full tubes. Half tubes are fixed on the side walls to simulate infinite array of tube bundle and minimize boundary layer effects. The transverse and longitudinal pitches, S_T and S_L, are 40 mm and 20 mm, respectively. The void fraction in the tube bundle is 0.607.

Fig. 1.

Schematic diagram of the test section with a staggered tube bundle.

The mean and fluctuating velocity components in x and y directions were measured by a hot-wire anemometer with a crossed-wire probe (Dantec. 55P62) at the inlet of the test section and 8 points along a centerline of the tube bundle as showed in Fig. 1. The probe with two ϕ5 mm diameter, 1.25 mm long platinum-plated tungsten wire sensors was placed vertical to the flow in the test section. The sampling frequency and the sampling time were fixed at 1 kHz and 8.192 s, respectively. The probe was calibrated for the velocity and the angular sensitivity at 20°C before the measurements. The hot-wire probe has the advantages of small size and high response to measure fluctuating velocity components in the tube bundle. The experimental conditions are summarized in Table 1. The steady and uniform flow was approaching the test section. The superficial velocity, u_∞, defined as the mean velocity at the cross section of the inlet is 6.88 m/s. Reynolds numbers, Re_∞, based on the tube diameter and the superficial velocity is 9160. The Reynolds number is defined as:   

$$R{e}_{\infty }=\frac{\rho u_{\infty }D}{\mu }.$$(1)

Kinetic energy of turbulence is obtained as:   

$$k=\frac{1}{2}\left(u^{\prime 2}+v^{\prime 2}\right).$$(2)

Table 1. Experimental conditions.

Fluid		Air
Temperature, T	[K]	293
Density, ρ	[kg/m3]	1.204
Viscosity, μ	[Pa·s]	18.09×10−6
Superficial velocity, u∞	[m/s]	6.88
Reynolds number, Re∞	[–]	9160
Void fraction, α	[–]	0.607

3. Microscopic Turbulence Model

3.1. Governing Equation

The standard k-ε turbulence model^8,9) is used to calculate the turbulent flows in the staggered tube bundle with a boundary fitted coordinate system. Figure 2 shows typical control volumes of a physical plane and a transformed plane. A collocated grid system in which Cartesian velocities and pressure (as well as other scalars) are stored at the cell center locations is employed in this calculation. The physical plane presented by the boundary fitted coordinate system must be transformed into the rectangular plane. Transport equations for steady two-dimensional flows can be expressed in the physical plane as:   

$$\frac{\partial }{\partial x}\left(\rho u\varphi \right)+\frac{\partial }{\partial y}\left(\rho v\varphi \right)=\frac{\partial }{\partial x}\left({\mathrm{\Gamma }}_{\varphi }\frac{\partial \varphi }{\partial x}\right)+\frac{\partial }{\partial y}\left({\mathrm{\Gamma }}_{\varphi }\frac{\partial \varphi }{\partial y}\right)+S_{\varphi },$$(3)

where ϕ represents the dependent variables that denote mass (1), momentum (u, v), kinetic energy of turbulence (k) and eddy dissipation rate of kinetic energy of turbulence (ε). Γ_ϕ and S_ϕ are the exchange coefficient and the source term associated with the variable, ϕ. After applying the transformation, the transport equations result in the following equations:   

$$\begin{array}{l}{\displaystyle \frac{\partial }{\partial \xi }}\left(\rho U\varphi \right)+{\displaystyle \frac{\partial }{\partial \eta }}\left(\rho V\varphi \right)={\displaystyle \frac{\partial }{\partial \xi }}\left[{\displaystyle \frac{{\mathrm{\Gamma }}_{\varphi }}{J}}\left(g^{11}{\displaystyle \frac{\partial \varphi }{\partial \xi }}+g^{12}{\displaystyle \frac{\partial \varphi }{\partial \eta }}\right)\right]\\ +{\displaystyle \frac{\partial }{\partial \eta }}\left[{\displaystyle \frac{{\mathrm{\Gamma }}_{\varphi }}{J}}\left(g^{21}{\displaystyle \frac{\partial \varphi }{\partial \xi }}+g^{22}{\displaystyle \frac{\partial \varphi }{\partial \eta }}\right)\right]+S{\prime }_{\varphi },\end{array}$$(4)

where U and V are the contravariant velocity components in ξ and η directions, respectively and obtained as:   

$$U={\beta }_{\xi }^{x}u+{\beta }_{\xi }^{y}v,$$(5)

  

$$V={\beta }_{\eta }^{x}u+{\beta }_{\eta }^{y}v,$$(6)

where β_i^j represent metrics which are given in the following equations:   

$${\beta }_{\xi }^x=\frac{\partial y}{\partial \eta },$$(7)

  

$${\beta }_{\xi }^y=-\frac{\partial x}{\partial \eta },$$(8)

  

$${\beta }_{\eta }^x=-\frac{\partial y}{\partial \xi },$$(9)

  

$${\beta }_{\eta }^y=\frac{\partial x}{\partial \xi }.$$(10)

Contravariant velocity components are defined at the cell face locations by interpolating Cartesian velocity components, as showed in Fig. 2. J is Jacobian, and g^ij are the contravariant metrics components which are shown in Eqs. (11), (12) and (13).   

$$g^{11}={\beta }_{\xi }^x{\beta }_{\xi }^x+{\beta }_{\xi }^y{\beta }_{\xi }^y$$(11)

  

$$g^{12}=g^{21}={\beta }_{\xi }^x{\beta }_{\eta }^x+{\beta }_{\xi }^y{\beta }_{\eta }^y$$(12)

  

$$g^{22}={\beta }_{\eta }^x{\beta }_{\eta }^x+{\beta }_{\eta }^y{\beta }_{\eta }^y$$(13)

Γ_ϕ and S’_ϕ are the exchange coefficient and the source term associated with the variable, ϕ, which are shown in Table 2. G_k in the source term of the transport equation for kinetic energy of turbulence is a product of kinetic energy of turbulence.

Fig. 2.

Typical control volume.

Table 2. Exchange coefficient and source terms for dependent variable, ϕ.

variable mass		ϕ	Γϕ	Sϕ
		1	0	0
				$-\left({\beta }_{\xi }^x\frac{\partial p}{\partial \xi }+{\beta }_{\eta }^x\frac{\partial p}{\partial \eta }\right)$
u-momentum		u	μeff	$+\frac{\partial }{\partial \xi }\left[\frac{{\mu }_{\text{eff}}}{J}{\beta }_{\xi }^x\left({\beta }_{\xi }^x\frac{\partial u}{\partial \xi }+{\beta }_{\eta }^x\frac{\partial u}{\partial \eta }\right)\right]+\frac{\partial }{\partial \eta }\left[\frac{{\mu }_{\text{eff}}}{J}{\beta }_{\eta }^x\left({\beta }_{\xi }^x\frac{\partial u}{\partial \xi }+{\beta }_{\eta }^x\frac{\partial u}{\partial \eta }\right)\right]$
				$+\frac{\partial }{\partial \xi }\left[\frac{{\mu }_{\text{eff}}}{J}{\beta }_{\xi }^y\left({\beta }_{\xi }^x\frac{\partial v}{\partial \xi }+{\beta }_{\eta }^x\frac{\partial v}{\partial \eta }\right)\right]+\frac{\partial }{\partial \eta }\left[\frac{{\mu }_{\text{eff}}}{J}{\beta }_{\eta }^y\left({\beta }_{\xi }^x\frac{\partial v}{\partial \xi }+{\beta }_{\eta }^x\frac{\partial v}{\partial \eta }\right)\right]$
				$-\left({\beta }_{\xi }^y\frac{\partial p}{\partial \xi }+{\beta }_{\eta }^y\frac{\partial p}{\partial \eta }\right)$
v-momentum		v	μeff	$+\frac{\partial }{\partial \xi }\left[\frac{{\mu }_{\text{eff}}}{J}{\beta }_{\xi }^x\left({\beta }_{\xi }^y\frac{\partial u}{\partial \xi }+{\beta }_{\eta }^y\frac{\partial u}{\partial \eta }\right)\right]+\frac{\partial }{\partial \eta }\left[\frac{{\mu }_{\text{eff}}}{J}{\beta }_{\eta }^x\left({\beta }_{\xi }^y\frac{\partial u}{\partial \xi }+{\beta }_{\eta }^y\frac{\partial u}{\partial \eta }\right)\right]$
				$+\frac{\partial }{\partial \xi }\left[\frac{{\mu }_{\text{eff}}}{J}{\beta }_{\xi }^y\left({\beta }_{\xi }^y\frac{\partial v}{\partial \xi }+{\beta }_{\eta }^y\frac{\partial v}{\partial \eta }\right)\right]+\frac{\partial }{\partial \eta }\left[\frac{{\mu }_{\text{eff}}}{J}{\beta }_{\eta }^y\left({\beta }_{\xi }^y\frac{\partial v}{\partial \xi }+{\beta }_{\eta }^y\frac{\partial v}{\partial \eta }\right)\right]$
kinetic energy of turbulence		k	$\frac{{\mu }_{\text{eff}}}{{\sigma }_k}$	J(GK−ρε)
Eddy dissipation rate		ε	$\frac{{\mu }_{\text{eff}}}{{\sigma }_{\epsilon }}$	$J\frac{\epsilon }{k}\left(C_{\epsilon 1}{G}_K-C_{\epsilon 2}\rho \epsilon \right)$
${\mu }_{\text{eff}}=\mu +\rho {\text{C}}_{\mu }\frac{k^2}{\epsilon }$
$G_k={\mu }_{\text{eff}}\left\{2{\left[\left(\frac{{\beta }_{\xi }^x}{J}\frac{\partial u}{\partial \xi }+\frac{{\beta }_{\eta }^x}{J}\frac{\partial u}{\partial \eta }\right)+\left(\frac{{\beta }_{\xi }^y}{J}\frac{\partial v}{\partial \xi }+\frac{{\beta }_{\eta }^y}{J}\frac{\partial v}{\partial \eta }\right)\right]}^2+{\left[\left(\frac{{\beta }_{\xi }^y}{J}\frac{\partial u}{\partial \xi }+\frac{{\beta }_{\eta }^y}{J}\frac{\partial u}{\partial \eta }\right)+\left(\frac{{\beta }_{\xi }^x}{J}\frac{\partial v}{\partial \xi }+\frac{{\beta }_{\eta }^x}{J}\frac{\partial v}{\partial \eta }\right)\right]}^2\right\}$
Cμ	Cε1				Cε2	σk		σε
0.09	1.44				1.92	1.0		1.3

3.2. Numerical Method

The finite volume method applied on the collocated grid system is employed to discretize the governing equations. Convective and diffusive terms are discretized by the power law scheme¹⁶⁾ corresponding with the first-order upwind or the second-order central difference scheme for the former and the second-order central difference scheme for the latter. SIMPLE¹⁷⁾ algorithm is used to deal with the pressure-velocity coupling between the continuity and the momentum equations to avoid checkerboard pressure field.

The numbers of grids for ξ and η directions are 1202 × 100 for the test section, and 144 × 20 for a unit of the tube bundle as showed in Fig. 3. The calculation is sufficiently accurate for this number of grids. The wall function model⁹⁾ is applied to the bottom and top of the test section and tube walls for transportation equations of momentum and kinetic energy of turbulence. The inlet condition for the velocity is prescribed by setting a uniform superficial velocity profile. The effects of the inlet kinetic energy of turbulence on the flow development in the tube bundle are examined in this study. Normalized inlet kinetic energy of turbulence, k^*_inlet, is set to 0.1, 1.0 and 10 in the calculation and defined as:   

$$k_{\text{inlet}}^{\text{*}}=\frac{k_{\text{inlet}}}{k_{\text{inlet},\text{exp}}},$$(14)

where k_inlet,exp is the inlet kinetic energy of turbulence obtained in the present experiments. The inlet eddy dissipation rate of kinetic energy of turbulence is obtained from the following equation:   

$${\epsilon }_{\text{inlet}}=\frac{k_{\text{inlet}}^{3/2}}{1.52D_{\text{eq}}},$$(15)

where D_eq is the equivalent diameter of the inlet. In this calculation D_eq is set equal to 0.133 m due to the dimensions of cross section of the inlet, 200 mm × 100 mm.

Fig. 3.

Grid arrangement of a unit of the tube bundle (144×20).

4. Microscopic Model Results and Discussion

4.1. Velocity and Pressure Drop

The microscopic model results will be shown with the experimental results. Figure 4 shows the axial distributions of normalized axial velocity, u^*, along centerline. u^* is defined as:   

$$u^{*}=\frac{u}{u_{\infty }},$$(16)

where u_∞ is the superficial velocity. The axial velocities in both experiment and calculation increase as the fluid enters into the tube bundle, since the cross-sectional area decreases in the tube bundle. The calculated maximum axial velocity is found at the center between each tube along the axial distance. The flow is fully developed in the tube bundle. The calculation results are in good agreement with the experimental results. The small differences at the entrance and downstream in the tube bundle would come from the well-known weakness in the original turbulence model for wake flow and momentum dissipation, causing small discrepancies in following pressure drop and kinetic energy of turbulence.

Fig. 4.

Axial distributions of axial velocity component along the centerline by microscopic turbulence model.

Table 3 shows the calculated pressure drop, Δp, across the tube bundle at k^*_inlet = 1.0 and the pressure drop given by the empirical equation of Jakob,¹⁾ Δp_Jakob, obtained from the relationship between the pressure drop and the friction coefficient respect to the number of the rows for the flow direction, the distance between the rows, and Reynolds number. (Error) is defined as:   

$$\left(\text{Error}\right)=\frac{\mathrm{\Delta }p-\mathrm{\Delta }{p}_{\text{Jakob}}}{\mathrm{\Delta }{p}_{\text{Jakob}}}.$$(17)

The calculated pressure drop is in good agreement with that by the empirical equation. It can be seen that the microscopic turbulence model predicts the pressure drop across the tube bundle well.

Table 3. Pressure drops across the tube bundle by microscopic turbulence model.

	Pressure drop, Δp [Pa]	(Error) [–]
Calc. (k*inlet = 1.0)	808	0.092
Eq. of Jakob1)	740	–

4.2. Kinetic Energy of Turbulence

Figure 5 shows the developments of the calculated normalized kinetic energy of turbulence, k^*, at k^*_inlet = 0.1–10 with the experimental results along the centerline. k^* is defined as:   

$$k^{*}=\frac{k}{u_{\infty }^2}.$$(18)

When the fluid enters into the tube bundle, the kinetic energies of turbulence increase rapidly in the calculations and experiment. The kinetic energy of turbulence obtained from the experiment developed more slowly than the one obtained from the calculation within a short axial distance from the entrance of the tube bundle (x/D < 10) whereas the kinetic energy of turbulence obtained from the calculation increased rapidly and then decreased. The further discussion of the difference around the entrance seems to be difficult since the number of the available experimental data is small due to the difficulties of the measurements in such a narrow space. In the calculation, the increase in the inlet kinetic energies of turbulence, k_inlet, results in an increase in the kinetic energies of turbulence, k, at the entrance of the tube bundle. At the downstream rows in the tube bundle (10 < x/D < 18), the kinetic energy of turbulence shows the similar distribution pattern between each tube and reaches a constant value. It can be found that the value of the kinetic energy of turbulence with fully developed flow in the tube bundle is barely affected by the inlet kinetic energy of turbulence. This is because the kinetic energy of turbulence generated by the presence of the tube bundle is nearly equal to its dissipation rate. The calculated results are in good agreement with the experimental results. After the fluid comes out from the tube bundle (x/D > 18), the kinetic energies of turbulence drop sharply.

Fig. 5.

Axial distributions of kinetic energy of turbulence along the centerline by microscopic turbulence model.

5. Macroscopic Numerical Model

5.1. Macroscopic Turbulence Model

Takeda and Lockwood¹³⁾ proposed a macroscopic turbulence model in packed beds with volume averaging and additional source terms for fluid-particle interactions. The volume-averaged transport equations for the kinetic energy of turbulence, k, and its eddy dissipation rate, ε, based on the standard k-ε model^8,9) are given in the following equations:   

$$\nabla \cdot \left(\alpha \rho uk\right)=\nabla \cdot \left(\alpha \frac{{\mu }_{\text{eff}}}{{\sigma }_k}\nabla k\right)+\alpha G_k+C_{k1}\beta {\left|\mathit{u}\right|}^2-\alpha \rho \epsilon ,$$(19)

  

$$\begin{array}{l}\nabla \cdot \left(\alpha \rho u\epsilon \right)=\nabla \cdot \left(\alpha {\displaystyle \frac{{\mu }_{\text{eff}}}{{\sigma }_{\epsilon }}}\nabla \epsilon \right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \frac{\epsilon }{k}}\left(\alpha C_{\epsilon 1}{G}_k+C_{\epsilon 2}{C}_{k1}\beta {\left|\mathit{u}\right|}^2-\alpha C_{\epsilon 2}\rho \epsilon \right),\end{array}$$(20)

where α is the void fraction in packed beds. These macroscopic turbulence equations, in contrast with the microscopic k-ε model l,^8,9) require the extra source terms. The terms with G_k in Eqs. (19) and (20) represent the production of kinetic energy of turbulence and its eddy dissipation rate due to the interactions between fluid and particles, respectively. C_k1 is a turbulence model constant. The inter-phase momentum transfer coefficient, β, is given by:   

$$\beta =\frac{\alpha }{\left|\mathbf{u}\right|}\frac{\mathrm{\Delta }p}{L},$$(21)

where L is the length of packed beds. These source terms correspond to the net production rate converted from pressure drop, Δp. If there is no packing material present, the values of these sources vanish, and then the standard k-ε model equation for free fluid flow is recovered.

Takeda and Lockwood¹³⁾ used the k-l one-equation model with Eq. (22) involving a turbulence mixing length, l, related to the packing structure instead of applying the transport equation for the eddy dissipation rate, Eq. (20).   

$$\epsilon =\frac{k^{2/3}}{l}$$(22)

They determined the constant, C_k1, in the comparison with the measurements for heat and mass transfer to obtain Peclet number of 11 for high Reynolds number flow. However, it is difficult to determine the mixing length in packed beds. Moreover, the constant, C_k1, is dependent of the mixing length, that is to say, the k-l model is not useful for any packed beds with various geometric structures.

Hence, we propose a new macroscopic k-ε two-equation model based on the k-l one-equation model.¹³⁾ The volume-averaged transport equation for the kinetic energy of turbulence is given by Eq. (20) using the constant, C_k1. On the other hand, the volume-averaged transport equation for the eddy dissipation rate is proposed as:   

$$\begin{array}{l}\nabla \cdot \left(\alpha \rho u\epsilon \right)=\nabla \cdot \left(\alpha {\displaystyle \frac{{\mu }_{\text{eff}}}{{\sigma }_{\epsilon }}}\nabla \epsilon \right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \frac{\epsilon }{k}}\left(\alpha C_{\epsilon 1}{G}_k+C_{\epsilon 2}{C}_{k1}\beta \left|\mathbf{u}\right|\sqrt{k}-\alpha C_{\epsilon 2}\rho \epsilon \right),\end{array}$$(23)

where C_k2 is a turbulence model constant. The third term on the right hand side of Eq. (23) is the production of eddy dissipation rate of kinetic energy of turbulence due to the presence of packed beds. These model constants, C_k1 and C_k2, can be determined from numerical experiments.

5.2. Constants in the Present Turbulence Model

Once the turbulence model constants, C_k1 and C_k2, are determined, the macroscopic turbulence model can be solved as the similar manner of the standard k-ε two-equation model. When the flow is fully developed in the tube bundle, the production rate of kinetic energy of turbulence is equal to its eddy dissipation rate, as showed in the microscopic calculations and the experiments in the tube bundle (Fig. 5). As this phenomenon can be seen in packed beds, Eqs. (19) and (23) of the macroscopic model can be reduced as:   

$$C_{k1}\beta {\left|\mathbf{u}\right|}^2=\alpha \rho \epsilon ,$$(24)

  

$$C_{\epsilon 2}{C}_{k1}\beta \left|\mathbf{u}\right|\sqrt{k}=\alpha C_{\epsilon 2}\rho \epsilon .$$(25)

By using Eq. (21), Eqs. (24) and (25) can be written in Eqs. (26) and (27), respectively.   

$$\rho \epsilon =C_{k1}\frac{\mathrm{\Delta }p}{L}\left|\mathbf{u}\right|$$(26)

  

$$\rho \epsilon =C_{k2}\frac{\mathrm{\Delta }p}{L}\sqrt{k}$$(27)

To determine the turbulence model constants, C_k1 and C_k2, numerical experiments are carried out in the two-dimensional tube bundles with the boundary fitted coordinate system by using the microscopic k-ε turbulence model.^8,9) Figure 6 shows a schematic diagram of the staggered tube bundle. The tube bundle consists of 15 staggered rows of tubes with an outside diameter of D = ϕ20 mm. The working fluid is air at 20°C. The flow is fully developed with the upstream length of 2000 mm before the fluid enters the tube bundle. Arrangements of the tube bundles in these numerical experiments, including the longitudinal and transverse pitches, S_L and S_T, and the void fraction, α, are summarized in Table 4. The void fraction, α, ranges from 0.564 to 0.836 by varying the longitudinal and transverse pitches, S_L and S_T.

Fig. 6.

Schematic diagram of the tube bundle for the determination of C_k1 and C_k2.

Table 4. Arrangements of the tube bundles for the determination of C_k1 and C_k2.

		Case 1	Case 2	Case 3	Case 4	Case 5
Tube diameter, D	[m]	0.020	0.020	0.020	0.020	0.020
Longitudinal pitch, SL	[m]	0.020	0.020	0.020	0.040	0.040
Transverse pitch, ST	[m]	0.036	0.040	0.048	0.036	0.048
Length of the tube bundle, L	[m]	0.300	0.300	0.300	.0580	.0580
Void fraction, α	[–]	0.564	0.607	0.673	0.782	0.836
Number of tube rows, N	[–]	15	15	15	15	15

The superficial velocity, u_∞, is varied from 5.04 m/s to 11.0 m/s for the determination of C_k1 and C_k2, and then Reynolds number, Re_∞, is varied from 6710 to 14600. The wall function model⁹⁾ is applied to tube walls for the momentum and energy transport equations. Symmetry conditions are applied to the bottom (y = 0) and the top (y = S_T) of the tube bundle because the microscopic calculations uses the isotropic turbulence model and provides symmetry results.

After the microscopic numerical experiments were conducted, the volume-integrated quantities in the tube bundles for α = 0.546–0.836 and Re_∞ = 6710–14600 are substituted into Eqs. (26) and (27) to determine the constants, C_k1 and C_k2, respectively. The relationships between (ρε) and [Δp/L·|u|] for C_k1, and (ρε) and [Δp/L·k^1/2] for C_k2 are plotted in Fig. 7, respectively. The values of C_k1 = 0.650 and C_k2 = 1.41 can be determined by fitting a straight line to all data. The effects of the eddy dissipation rate, ε, are considered in the constant, C_k1, used in the transport equation for the kinetic energy of turbulence. On the other hand, the effects of the kinetic energy of turbulence, k, and its eddy dissipation rate, ε, are considered in the constant, C_k1, used in the transport equation for the eddy dissipation rate.

Fig. 7.

Determination of value for C_k1 and C_k2.

5.3. Governing Equation

The proposed macroscopic turbulence model is used to calculate the turbulent flow in the tube bundle. The void fractions, α, in a domain of the tube bundle and the free space in Fig. 1 are set equal to 0.607 and 1.0, respectively. The macroscopic governing equations for the flow can be obtained by volume averaging the corresponding microscopic equations. The transport equations used in the turbulent flow can be expressed in a two-dimensional rectangular coordinate system as:   

$$\frac{\partial }{\partial x}\left(\alpha \rho u\varphi \right)+\frac{\partial }{\partial y}\left(\alpha \rho v\varphi \right)=\frac{\partial }{\partial x}\left(\alpha {\mathrm{\Gamma }}_{\varphi }\frac{\partial \varphi }{\partial x}\right)+\frac{\partial }{\partial y}\left(\alpha {\mathrm{\Gamma }}_{\varphi }\frac{\partial \varphi }{\partial y}\right)+S_{\varphi },$$(28)

where ϕ represents the dependent variables that denote mass (1), momentum (u, v), kinetic energy of turbulence (k) and eddy dissipation rate of kinetic energy of turbulence (ε). Γ_ϕ and S_ϕ are the exchange coefficient and the source term associated with the variable, ϕ, as showed in Table 5. The last terms in the source terms of the momentum equations are derived from the fluid-particle interactions. The pressure drop, Δp, across the tube bundle, which is used in the inter-phase momentum transfer coefficient, β, as showed in Eq. (21) is given by:   

$$\mathrm{\Delta }p=4f_{\text{TB}}\left(\frac{\rho u_{\text{max}}^2}{2}\right)N,$$(29)

where f_TB is the friction coefficient in the tube bundle, u_max is the mean velocity in the minimum cross-sectional area and N is the number of tube rows along the flow direction. Empirical equation of Jakob¹⁾ can be used for the friction coefficient, f_TB, for the turbulent flow though the tube bundle.

Table 5. Exchange coefficient and source terms for dependent variable, ϕ.

variable mass		ϕ		Γϕ	Sϕ
		1		0	0
u-momentum		u		μeff	$-\alpha \frac{\partial p}{\partial x}+\frac{\partial }{\partial x}\left(\alpha {\mu }_{\text{eff}}\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left(\alpha {\mu }_{\text{eff}}\frac{\partial v}{\partial x}\right)-\beta u$
v-momentum		v		μeff	$-\alpha \frac{\partial p}{\partial y}+\frac{\partial }{\partial x}\left(\alpha {\mu }_{\text{eff}}\frac{\partial u}{\partial y}\right)+\frac{\partial }{\partial y}\left(\alpha {\mu }_{\text{eff}}\frac{\partial v}{\partial y}\right)-\beta u$
kinetic energy of turbulence		k		$\frac{{\mu }_{\text{eff}}}{{\sigma }_k}$	$\alpha \left(G_k+C_{k1}\frac{\beta }{\alpha }{\left|\mathbf{u}\right|}^2-\rho \epsilon \right)$
Eddy dissipation rate		ε		$\frac{{\mu }_{\text{eff}}}{{\sigma }_{\epsilon }}$	$\alpha \frac{\epsilon }{k}\left(C_{\epsilon 1}{G}_k+C_{\epsilon 2}{C}_{k2}\frac{\beta }{\alpha }\left|\mathbf{u}\right|\sqrt{k}-C_{\epsilon 2}\rho \epsilon \right)$
${\mu }_{\text{eff}}=\mu +\rho {\text{C}}_{\mu }\frac{k^2}{\epsilon }$
$G_k={\mu }_{\text{eff}}\left\{2\left[{\left(\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]+{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2\right\}$
Cμ	Ck1		Ck2			Cε1	Cε2	σk	σε
0.09	0.650		1.41			1.44	1.92	1.0	1.3

5.4. Numerical Method

The finite volume method applied on a staggered grid system in which velocities are stored at the cell face locations and pressure and other scalars are stored at the cell center locations is employed to discretize the governing equations. Convective and diffusive terms are discretized by the power law scheme¹⁶⁾ corresponding with the first-order upwind or the second-order central difference scheme for the former and the second-order central difference scheme for the latter. SIMPLE¹⁸⁾ algorithm is used to deal with the pressure-velocity coupling between the continuity and the momentum equations.

The numbers of grids for x and y directions are 360 × 40 for the test section. The calculation is sufficiently accurate for this number of grids. The symmetry condition is applied to the bottom and top of the test section for the momentum and energy transport equations, since the macroscopic turbulence model can be used for one-dimensional and fully developed flow. The inlet conditions for the velocity and the eddy dissipation rate are prescribed by setting a uniform superficial velocity profile and by using Eq. (15), respectively, as well as the microscopic numerical method. The effects of the inlet kinetic energy of turbulence on the flow development in the tube bundle are examined in this study as well as the microscopic numerical study. Normalized inlet kinetic energy of turbulence, k^*_inlet, is set to 0.1, 1.0 and 10 in the calculation.

6. Macroscopic Model Results and Discussion

6.1. Velocity and Pressure Drop

The macroscopic model results will be shown for the investigation of turbulent flows in the tube bundle and the model validation by comparing with the microscopic model results and the experimental results. Figure 8 shows the axial distributions of normalized axial velocities, u^*, along the centerline by the macroscopic turbulence model, the microscopic turbulence model and the experiment. As the void fraction is set to 0.607 uniformly in the tube bundle (expect for at the entrance and exit of the tube bundle), a uniform distribution is obtained for the axial velocity in the macroscopic model result. Note that the microscopic result is shown only along the centerline of the tube bundle. Actually, a non-uniform distribution is achieved for velocity components in x and y directions. Also, the experimental results are obtained at the measuring points between each tube along the centerline, as shown in Fig. 1.

Fig. 8.

Axial distributions of axial velocity components along the centerline.

Table 6 shows the calculated pressure drop, Δp, across the tube bundle by the macroscopic turbulence model at k^*_inlet = 1.0 and the pressure drop given by the empirical equation of Jakob, Δp_Jakob.¹⁾ The calculated pressure drop is in good agreement with that by the empirical equation. It can be seen that the macroscopic turbulence model predicts the pressure drop across the tube bundle well.

Table 6. Pressure drops across the tube bundle by macroscopic turbulence model.

	Pressure drop, Δp [Pa]	(Error) [–]
Calc. (k*inlet = 1.0)	757	0.023
Eq. of Jakob1)	740	–

6.2. Kinetic Energy of Turbulence

Figure 9 shows the axial distributions of the calculated normalized kinetic energy of turbulence, k^*, at k^*_inlet = 0.1–10 along the centerline by the macroscopic turbulence model. Microscopic model result (at k^*_inlet = 1.0) and the experimental ones are also shown in Fig. 9 for the comparison. After the fluid enters into the tube bundle, the kinetic energies of turbulence increase sharply in the macroscopic calculations. As the fluid goes downstream, the kinetic energies of turbulence reach a constant value at any case of the inlet kinetic energies of turbulence. It can be found that the inlet kinetic energy of turbulence has no effect on the kinetic energy of turbulence in the fully developed flow region (x/D > 10) by using the macroscopic turbulence model as well as the microscopic turbulence model as showed in Subsection 5.2. The macroscopic model results appear to be in reasonably good agreement with the microscopic model result and the experimental ones. This agreement indicates the macroscopic turbulence model is valid for the prediction of kinetic energy of turbulence in fully developed flow in the tube bundle. The turbulence model constants, C_k1 and C_k2, are contributed to the good prediction. This is because C_k1 includes the effect of the eddy dissipation rate in the tube bundle, which is used in the transport equation for the kinetic energy of turbulence. C_k2, on the other hand, includes the effect of kinetic energy of turbulence in the tube bundle, which is used in the transport equation for the eddy dissipation rate. The interaction between the kinetic energy of turbulence and its dissipation rate should be considered for the prediction of turbulence flows in tube bundles. The reason why the macroscopic model overpredicts the kinetic energy of turbulence at a short axial distance from the entrance of the tube bundle where the flow is not fully developed is that the constants are determined when the fully developed flow is assumed. Considering the fluid flow in the packed bed, it may be applied to packed beds to predict the turbulent flows if an equation proposed by Ergun¹⁹⁾ is used for the inter-phase momentum transfer coefficient, β, instead of equation of Jakob¹⁾ for tube bundles.

Fig. 9.

Axial distributions of kinetic energy of turbulence along the centerline in the tube bundle by macroscopic turbulence model.

7. Conclusions

Experimental and numerical studies of turbulent flows were performed in a staggered tube bundle. Experimental data of axial velocity and kinetic energy of turbulence were obtained in the tube bundle by a hot-wire anemometer. The standard k-ε turbulence model was used for a microscopic numerical investigation in the tube bundle by using a boundary fitted coordinate system. The microscopic numerical results reveal that the predictions of velocity, pressure drop and kinetic energy of turbulence distributions were in good agreement with the measured values or the empirical equation. As the flow is fully developed in the tube bundle, kinetic energy of turbulence remains a constant value. The effects of the inlet kinetic energy of turbulence on the development of kinetic energy of turbulence were examined in the tube bundle. The kinetic energy of turbulence in the tube bundle is not affected by the inlet kinetic energy of turbulence.

The macroscopic numerical model was proposed with volume averaging and additional source terms for fluid-particle interactions. The constants in the macroscopic model were determined by conducting the microscopic calculations. As well as the microscopic results, the macroscopic numerical results show that the predictions of velocity, pressure drop and kinetic energy of turbulence are in good agreement with the measured values or the empirical equation. Also, the kinetic energy of turbulence in the tube bundle is not affected by the inlet kinetic energy of turbulence in the fully developed flow region. Therefore, the macroscopic turbulence model proposed in this study may be useful for predicting the turbulent flows in the tube bundles or packed beds.

Acknowledgement

This work was supported by JSPS KAKENHI Grant Number JP15560648.

Nomenclature

C_k1, C_k2: constants in k-ε two-equation model in a packed bed [–]

C_μ, Cε₁, Cε₂: constants in standard k-ε two-equation model [–]

D: tube diameter [m]

D_eq: equivalent diameter of inlet [m]

f_TB: friction coefficient in a tube bundle [–]

g^ij: contravariant metrics component [–]

G_k: product of kinetic energy of turbulence k [kg/(m·s³)]

J: Jacobian [–]

k: kinetic energy of turbulence [m²/s²]

L: length of a packed bed or a tube bundle [m]

l: mixing length [m]

m: mass fraction [–]

N: the number of tubes along flow direction [–]

p: pressure [Pa]

Δp: pressure drop [Pa]

Re_∞: Reynolds number based on superficial velocity [–]

S_L, S_T: longitudinal and transverse pitches in a tube bundle [m]

S_ϕ, S’_ϕ: volumetric source terms for dependent variable ϕ

T: temperature [K]

U, V: contravariant velocity component [m/s]

u, v: Cartesian velocity component [m/s]

u_max: mean velocity in minimum cross-sectional area [m/s]

u_∞: superficial velocity [m/s]

u’, v’: fluctional velocity component [m/s]

x, y: Cartesian coordinate [m]

Greek symbols

α: void fraction [–]

β: inter-phase momentum transfer coefficient [kg/(m³·s)]

β_ξ^x, β_ξ^y, β_η^x, β_η^x: metrics [–]

Γ_ϕ: exchange coefficient for dependent variable ϕ [kg/(m·s)]

ε: eddy dissipation rate [m²/s³]

ξ, η: boundary fitted coordinate [m]

μ: fluid viscosity [Pa·s]

ρ: fluid density [kg/m³]

σ_ϕ: Prandtl number [–]

ϕ: dependent variable

Subscripts

eff: effective

exp: experiment

inlet: inlet

l: laminar

t: turbulent

Superscript

*: normalized

References

• 1)   M.  Jacob: Mech. Eng. (Am. Soc. Mech. Eng.), 60 (1938), 381.
• 2)   O.  Simonin and  M.  Barcouda: Proc. in 4th Int. Symp. on Applications of Laser Anemometry to Fluid Mechanics, Lisbon, (1988), 5.23
• 3)   S.  Balabani and  M.  Yianneskis: Proc. Inst. Mech. Eng. Part C, J. Mech. Eng. Sci., 210 (1996), 317.
• 4)   C.  Iwaki,  K. H.  Cheong,  H.  Monji and  G.  Matsui: Exp. Fluids, 37 (2004), 350.
• 5)   S. S.  Paul,  M. F.  Tachie and  S. J.  Ormiston: Int. J. Heat Fluid Flow, 28 (2007), 441.
• 6)   S.  Aiba,  H.  Tsuchida and  T.  Ota: Bull. JSME, 25 (1982), 927.
• 7)   M.  Zijlema,  A.  Segal and  P.  Wesseling: Comput. Fluids, 24 (1995), 209.
• 8)   W. P.  Jones and  B. E.  Launder: Int. J. Heat Mass Transf., 15 (1972), 301.
• 9)   B. E.  Launder and  D. B.  Spalding: Comput. Methods Appl. Mech. Eng., 3 (1974), 269.
• 10)   J. K.  Watterson,  W. N.  Dawes,  A. M.  Savill and  A. J.  White: Int. J. Heat Fluid Flow, 20 (1999), 581.
• 11)   Z.  Yang and  T.  Shih: Near-Wall Turbulent Flows, ed. by R. M. C. So et al., Elsevier, Amsterdam, (1993), 165.
• 12)   Y. Q.  Wang,  P.  Jackson and  T.  Phaneuf: J. Thermophys. Heat Transf., 20 (2006), 738.
• 13)   K.  Takeda and  F. C.  Lockwood: Tetsu-to-Hagané, 82 (1996), 486.
• 14)   A.  Nakayama and  F.  Kuwahara: J. Fluids Eng., 121 (1999), 427.
• 15)   M. H. J.  Pedras and  M. J. S.  de Lemos: Int. J. Heat Mass Transf., 44 (2001), 1081.
• 16)   S.  Patankar: Numerical Heat Transfer and Fluid Flow, Series in Computational Methods in Mechanics and Thermal Sciences, CRC Press, New York, (1980), 79.
• 17)   C. M.  Rhie and  W. L.  Chow: AIAA J., 21 (1983), 1525.
• 18)   S. V.  Patankar and  D. B.  Spalding: Int. J. Heat Mass Transf., 15 (1972), 1787.
• 19)   S.  Ergun: Chem. Eng. Prog., 48 (1952), 89.