Numerical simulation of flow past circular duct

Ze-gao YIN*, Xian-wei Cao, Hong-da SHI, Jian MA

1. Ocean Engineering Key Laboratory of Shandong Province, Ocean University of China,Qingdao 266100, P. R. China

2. Civil Engineering Department, Zhejiang University, Hangzhou 310027, P. R. China

1 Introduction

Pipelines are the main way that oil and gas are transported through the sea. The exterior flow past horizontal circular ducts is an important factor that many scholars and engineers are concerned about. In order to investigate flow near the circular cylinder, physical models have been intensively employed. Bearman and Zdravkovich (1978)used a wind tunnel test to investigate flow around a circular cylinder when the Reynolds numbers (Re)were 2.5 × 104and 4.8 × 104, and the gap ratio between the circular cylinder and the plate wall varied from 0 to 3.5. Then, the smoke-wire method was used to express the flow visualization and pressure distribution, and it was found that the Strouhal number remained constant when the gap ratio was greater than 0.3. On the basis of Bearman’s research, Freds?e and Hansen (1987)conducted flume experiments to study flow separation, drag force, and lift force of the cylinder.According to the variable diameter cylinder, flow past the cylinder near the bottom boundary was investigated with the experiments. Qi et al. (2006)analyzed the characteristics of the vortex of flow past the cylinder on the basis of physical model data when the gap ratio was 0.5.

With the increase of computational power, mathematical models have been utilized to simulate flow past the circular duct. Williamson (1996), Travin et al. (2000), Kravchenko and Moin (2000), and Sarpkaya (2004)used the Reynolds-averaged Navier-Stokes equations to study the flow near a circular cylinder. Kalro and Tezduyar (1997), Breuer (2000), Jordan(2002), and Catalano et al. (2003)used the large eddy simulation method to study the flow around a cylinder. This research focused mainly on the hydraulic characteristics and the fluid structure near the cylinder.

As far as direct numerical simulation for Reynolds numbers from 10 to 1 000 is concerned,Henderson (1997)improved the capabilities of research on flow transition from two dimensions to three dimensions. Mittal and Balachandar (1995)investigated the three-dimensional effect on simulation accuracy at Re = 525.

Although these scholars have concentrated on flow past the cylinder, there have been few studies on the flow force on circular ducts with various upper and bottom gap ratios. In this study, the Renormalization Group (RNG)k-ε turbulence model and Volume of Fluid (VOF)method were employed to calculate the hydraulic parameters around a circular duct with various upper and bottom gap ratios, and the forces on the circular duct were also computed and analyzed.

2 Mathematical model

The circular duct boundary causes the streamline to curve intensively. Compared with the standard k-ε model, the RNG k-ε model more accurately simulates complex boundary flow. Therefore, the RNG k-ε model was employed in this study. The tensor form of the Navier-Stokes equations and the RNG k-ε equations can be written as follows (Liu et al.2004; Baumal et al. 1998; Bazdidi et al. 2004):

Continuity equation:

Momentum equation:

k equation:

ε equation:

where ρ is fluid density; μ is molecular viscosity; uiis the velocity component in the xi(i =1, 2)coordinate direction, denoting x and y coordinates, respectively, when i is equal to 1 and 2;is pressure in the xidirection; k is turbulent kinetic energy; ε is the turbulent dissipation rate;Fiis the body force in the xidirection: F1=0 and F2=g; δijis the Kronecker delta; when i = j,and when

3 Numerical methods and boundary conditions

3.1 Numerical methods

The VOF method was adopted to handle the free surface. Its concept can be expressed as follows: a water volume fraction function (WVFF), Fw, is defined to present the relative water volume fraction in every computational grid cell. In each cell, the WVFF plus the volume fraction of air, Fa, is 1:

For a given cell, three situations exist: Fw=1, when the cell is completely filled with water; Fw=0, when the cell is completely filled with air; and 0＜Fw＜1, which indicates that the cell is filled with both water and air, and an interface between the two exists in the cell.Free surface flow conforms to the third situation. The control equation of Fwcan be expressed as follows:

where ρwand ρa(bǔ)are the densities of water and air, respectively; μwand μaare the molecular viscosity coefficients of water and air, respectively; and ρ and μ can be obtained by Eqs. (7)and (8)after Fwhas been computed.

After the governing equations of the turbulence model consisting of the coupled RNG k-ε model and the VOF method, including Eqs. (1)through (4)and Eqs. (6)through (8), are solved, the hydraulic factors can be obtained.

In this study, the mesh generation was performed with GAMBIT 2.2.30. In the Fluent software, a two-dimensional, single-precision, implicit, segregated solver was selected. The PLIC-VOF algorithm was used for the reconstruction of the interface, and the semi-implicit method for a pressure-linked equation algorithm was adopted for the solution of discrete equations. The solutions were considered to have converged when residuals were less than 10-3. The simulations were performed on a PC with a 2.8 GHz processor and one GB of RAM.

3.2 Boundary conditions

The inlet boundary was divided into a water velocity inlet and an air pressure boundary.The water velocity at the inlet boundary was defined as U0. The air boundary was regarded as atmospheric pressure. The k and ε at the inlet boundary were evaluated according to the following formulas:

where H0is the water depth at the inlet boundary. The outlet boundary was regarded as the atmospheric pressure boundary. At the wall, a no-slip boundary was specified, and fluid velocities were set to zero. The wall function was adopted to treat the viscous sub-layer.

4 Numerical results and analysis

4.1 Model validation

In Fig. 1, we can see that D is the diameter of the circular duct; e is the distance between the lower end of the circular duct and the bottom plane; H is the water depth; U is the mean flow velocity; Re=UDν, with ν being the kinematic viscosity; E is defined as the bottom gap ratio: E=e D; and G is defined as the upper gap ratio: G = (H - D- e)D.

In order to validate the mathematical model, a physical model and a mathematical model were established where Re = 5 × 105, E = 1, G = 3, and D = 0.2 m. Velocity measurement was conducted with an acoustic Doppler velocimeter (ADV), as in Fig. 2. Pressure was measured by pressure sensors. The computational area was partitioned by unstructured grids, whose cell size was 0.1 m, and the total number of grid cells was 3 477. Due to the hydraulic complexity near the circular duct boundary, the area near the duct was partitioned with refined grids (Fig. 3).Then, the data of the physical model were compared with computed results.

Fig. 1 Sketch map of computational case

Fig. 2 Velocity measurement sketch of physical model with ADV

Fig. 3 Computational grids of case where E = 1 and G = 3

From Fig. 4, we can see that the measured free surface agrees with the computational results for the most part. From Table 1, we can see that all the relative errors between physical model velocities (Up)and computational velocities (Uc)are less than 7%, and measured data approach computational results. From Table 2, we can see that all the relative errors between physical model pressures (pp)and computational pressures (pc)are less than 6%, and measured data also approach computational results. Thus, computations of flow past the circular duct with the RNG k-ε turbulent model and the VOF method are valid and reliable.

Fig. 4 Comparison between computational free surface and measured free surface (E = 1, G = 3, x is the horizontal distance from the duct center, and H is the water depth)

Table 1 Comparison between Up and Uc

Table 2 Comparison between pp and computational pc

4.2 Computational velocities and pressure

When Re = 5 × 105, E varies among the values 0, 1, and 2; G varies among the values 1, 2,3, and 4; and we obtain 12 conditions with different values of E and G (Table 3). The mathmatical model was used to compute the hydraulic factors associated with the 12 conditions.The results of the 7th condition, where E = 1 and G = 3, were chosen to be analyzed.

Table 3 Twelve computational conditions

From Fig. 5 we can see that the velocity direction is significantly deflected near the front of the circular duct boundary owing to the squeezing action that the circular duct boundary exerts on the flow, and the magnitude of the velocity becomes larger. The velocity distribution is symmetrical in the wake area, and the magnitudes of velocities are less than those of other areas.

From Fig. 6 we can see that the pressure along the duct boundary is unevenly distributed.The pressure on the duct close to the inflow direction and the bottom plane is the largest, and the pressure on the duct close to the free surface is the smallest. The pressure recovers faster than the velocity does. Pressure contours approximately parallel each other in the section that is about five times the duct diameter away from circular duct, and the contour direction is almost horizontal, showing that the influence of the duct on pressure can be ignored here.

Fig. 5 Computational flow velocities (E = 1, G = 3)

Fig. 6 Computational pressure (Pa)(E = 1, G = 3)

4.3 Computational turbulent kinetic energy and turbulence dissipation rate

It can be seen in Fig. 7 and Fig. 8 that the distribution of turbulent kinetic energy approximately agrees with that of the turbulence dissipation rate. The largest area lies on a duct close to the inflow and the free surface, and the turbulent stress is the largest here. In the wake area, the turbulent intensity is the lowest. The turbulent kinetic energy and the turbulence dissipation rate are distributed symmetrically, and the horizontal line across the duct center is the symmetry axis.

Fig. 7 Turbulent kinetic energy (m2/s2)(E = 1, G = 3)

Fig. 8 Turbulence dissipation rate (m2/s3)(E = 1, G = 3)

4.4 Analysis of force on duct

The force on the duct can be decomposed into a horizontal drag force and a vertical lift force. The drag force coefficient (CD), lift force coefficient (CL), composite force (F), and azimuth (α)can be defined as follows:

where Fxand Fyare the horizontal drag force and vertical lift force, respectively.

In Fig. 9, Fig. 10, Fig. 11, and Fig. 12, we can see that, when E = 0, CD, CL, and F are at their largest values, but α is at its smallest. If we do not consider the duct dead weight and the internal flow force, the circular duct is most likely to break when E = 0. With an increase of E from 0 to 1, CDand F decrease sharply, and CLdoes not decreases so much, but α increases dramatically. With an increase of E from 1 to 2, CD, CL, F, and α vary slightly.With an increase of G from 1 to 4, CD, CL, and F decrease, but α increases when E = 0.More specifically, with the increase of G from 1 to 2 when E = 0, CD, CL, and F decrease sharply, but α increases quickly, showing that lift force contributes more and more to composite force. With the continuing increase of G from 2 to 4 when E = 0, CD, CL, and F decrease slowly. The value of α, however, also increases slowly. With an increase of G from 1 to 4 when E = 1 and E = 2, CD, CL, F, and α vary slightly. More specifically, CDvaries from 1.5 to 1.7, CLvaries from 60 to 61, F varies from 7 550 N to 7 620 N, and α varies from 88.3° to 88.6°, showing that the variation of G has little influence on CD, CL, F, and α when E = 1 and E = 2.

Fig. 9 Relation curve of CD and G

Fig. 10 Relation curve of CL and G

Fig. 11 Relation curve of F and G

Fig. 12 Relation curve of α and G

5 Conclusions

(1)Using Fluent software, the RNG k-ε turbulence model and the VOF method were employed to simulate the flow past a circular duct. The computational results regarding related hydraulic parameters, including surface position, velocities, pressure, turbulent kinetic energy,and the turbulence dissipation rate, were obtained and analyzed.

(2)The forces on the duct were computed and analyzed with different bottom gap ratios and upper gap ratios. When E = 0, the drag force coefficient, lift force coefficient, and composite force are at their largest, but the azimuth is at its smallest, compared to other conditions. When E = 0, the drag force coefficient, lift force coefficient, and composite force decrease, but the azimuth increases along with the upper gap ratio. With an increase of the bottom gap ratio from 0 to 1, the drag force coefficient and the composite force decrease sharply, and the lift force coefficient does not decrease so much, but the azimuth increases dramatically.

(3)The bottom gap ratio is the most important factor influencing the force on the circular duct. When the bottom gap ratio is less than 1, the upper gap ratio has a significant influence on the force of the circular duct. When the bottom gap ratio is greater than 1, the variation of the upper gap ratio has little influence on the force of the circular duct.

Baumal, A. E., McPhee, J. J., Calamai, P. H., and Mompean, G. 1998. Numerical simulation of a turbulent flow near a right-angled corner using the Speziale non-linear model with RNG k-ε equations.Computers and Fluids, 27(7), 847-859. [doi:10.1016/S0045-7930(98)00004-8]

Bazdidi, T. F, Shahmir., A, and Haghparast, K. A. 2004. Numerical analysis of a single row of coolant jets injected into a heated crossflow. Journal of Computational and Applied Mathematics, 168(1-2), 53-63.[doi:10.1016/j.cam.2003.05.012]

Bearman, P. W, and Zdravkovich, M. M. 1978. Flow around a circular cylinder near a plane boundary. Journal of Fluid Mechanics, 89(1), 33-47. [doi:10.1017/S002211207800244X]

Breuer, M. 2000. A challenging test case for large eddy simulation: High Reynolds number circular cylinder flow. International Journal of Heat and Fluid Flow, 21(5), 648-654.

Catalano, P., Wang, M., Iaccarino, G., and Moin, P. 2003. Numerical simulation of the flow around a circular cylinder at high Reynolds numbers. International Journal of Heat and Fluid Flow, 24(4), 463-469. [doi:10.1016/S0142-727X(03)00061-4]

Freds?e, J., and Hansen, E. A. 1987. Lift forces on pipelines in steady flow. Journal of Waterway, Port,Coastal and Ocean Engineering, 113(2), 139-155. [doi:10.1061/(ASCE)0733-950X(1987)113:2(139)]

Henderson, R. D. 1997. Nonlinear dynamics and pattern formation in turbulent wake transition. Journal of Fluid Mechanics, 352(1), 65-112. [doi:10.1017/S0022112097007465]

Jordan, S. A. 2002. Investigation of the cylinder separated shear-layer physics by large-eddy simulation.International Journal of Heat and Fluid Flow, 23(1), 1-12.

Kalro, V., and Tezduyar, T. 1997. Parallel 3D computation of unsteady flows around circular cylinders.Parallel Computing, 23(9), 1235-1248. [doi:10.1016/S0167-8191(97)00050-1]

Kravchenko, A. G., and Moin, P. 2000. Numerical studies of flow over a circular cylinder at ReD= 3900.Physics of Fluids, 12(2), 403-417. [doi:10.1063/1.870318]

Liu, P. Q., Qu, Q. L., Wang, Z. G., and Zhang, H. M. 2004. Numerical simulation on hydrodynamic characteristics of bifurcation pipe with internal crescent rib. Journal of Hydraulic Engineering, 35(3),42-46. (in Chinese)

Mittal, R., and Balachandar, S. 1995. Effect of three-dimensionality on the lift and drag of nominally two-dimensional cylinders. Physics of Fluids, 7(8), 1841-1865.

Qi, E. R., Li G. Y., Li, W., Wu, J., and Zhang, X. 2006. Study of vortex characteristics of the flow around a horizontal circular cylinder at various gap-ratios in the cross-flow. Journal of Hydrodynamics, Series B,18(3), 334-340.

Sarpkaya, T. 2004. A critical review of the intrinsic nature of vortex-induced vibrations. Journal of Fluids and Structures, 19(4), 389-447. [doi:10.1016/S0889-9746(04)00035-0]

Travin, A., Shur, M., Strelets, M., and Spalart, P. 2000. Detached-eddy simulations past a circular cylinder.Flow, Turbulence and Combustion, 63(1-4), 293-313. [doi:10.1023/A:1009901401183]

Williamson, C. H. K. 1996. Vortex dynamics in the cylinder wake. Annual Review of Fluid Mechanics, 28(1),477-539. [doi:10.1146/annurev.fl.28.010196.002401]