168

2017,29(1):168-171

DOI: 10.1016/S1001-6058(16)60728-X

Numerical simulation of flow through circular array of cylinders using porous media approach with non-constant local inertial resistance coefficient* Jie-min Zhan (詹杰民)1, Wen-qing Hu (胡文清)1, Wen-hao Cai (蔡文豪)1, Ye- jun Gong (龚也君)1, Chi-wai Li (李志伟)2 1. Department of Applied Mechanics and Engineering, School of Engineering, Sun Yat-Sen University, Guangzhou 510275, China, E-mail: [email protected] 2. Department of Civil and Environmental Engineering, Hong Kong Polytechnic University, Hong Kong, China (Received November 29, 2016, Revised December 15, 2016) Abstract: Aquatic vegetation zone is now receiving an increasing attention as an effective way to protect the shorelines and riverbeds. To simulate the flow through the vegetation zone, the vegetation elements are often simplified as equidistant rigid cylinders, and in the whole zone, the porous media approach can be applied. In this study, a non-constant inertial resistance coefficient is introduced to model the unevenly distribution of the drag forces on the cylinders, and an improved porous media approach is applied to one circular array of cylinders positioned in a 2-D flume. The calculated velocity profile is consistent with the experimental data. Key words: Vegetation zone, cylinder array, porous media, inertial resistance coefficient

The aquatic vegetation is widely researched due to the significant influence on ecosystem and hydraulics. On the aspect of laboratory experiments, a series of classic experiments were conducted by Ghisalberti and Nepf[1] and Nezu and Onitsuka[2]. In their studies, the vegetation is simulated by rigid bodies, and the whole vegetation zone is regarded as an array of solid cylinders. The tracking dye is used to track the flow trajectory, showing the complex flow pattern, turbulence and diffusion in the flows through the vegetation. For the numerical simulation, two main methods[3] are widely employed to simulate the flow through the vegetation. The first method, which resolves the flow around each vegetation unit, is called the multi-body model[4]. Each vegetation unit is generally * Project supported by the Special Fund of Marine-Fishery Science-Technology Extension in Guangdong Province (Grant No. A201401B08). Biography: Jie-min Zhan (1963-), Male, Ph. D., Professor Corresponding author: Ye-jun Gong, E-mail: [email protected]

simplified as a cylinder and each one of them is taken into consideration in the model. The second method, the porous model[5], treats the entire vegetation zone as a porous zone and the impact of the vegetation on the flow is governed by the relevant parameters of the porous zone. Compared with the porous model, the multi-body model considers the influence of each single vegetation unit, therefore, is more accurate and time consuming. Relatively speaking, as the porous model regards the vegetation as one porous zone, the computational efficiency is better, therefore, it is more widely used in the real engineering applications. In the porous model, the effect of vegetation is determined by the strength of the resistance source in the momentum equations, hence the determination of a suitable drag coefficient Cd is the key for a successful simulation[3,6]. The drag coefficient Cd

[7]

is

generally defined as Cd = F /(0.5 r U Apro ) , where F 2

is the time-averaged drag force over the vegetable region, U is the characteristic velocity, which is commonly defined as the time-averaged velocity inside the vegetable region, and Apro is the projection area of the vegetable region in the flow direction. Yu et al.[3]

169

proposed an improvement in the definition of the drag force for the multi-cylinder model and obtained a constant value of Cd . However, the drag force acting on the vegetation is not constant in the vegetation zone, and the assumption of a constant drag coefficient may lead to errors. especially, when the vegetation zone is large or the density of plant is high. In this study, we build a new equivalent porous model with a non-constant drag coefficient for the large-scale engineering applications. For the flow of incompressible, homogeneous and viscous fluid, the governing Navier-Stokes equations can be written as ∇⋅v = 0 ,

∂v 1 + (v ⋅ ∇)v = Fb − ∇ p + ν ∆ v ρ ∂t

(1)

where v is the velocity vector, ρ is the fluid density,

C

2ρ k

σΦ

 1 ∂ω ∂ω 1 ∂k ∂k max  2 ,  ω ∂x ∂x k 2 ∂x ∂x j j j j 

   , 0   

(3)

Here L is the length scale of the modeled turbulence L = k /(c1/µ 4ω ) , and the model parameters in Eq.(3) are η2 = 3.51 , σ Φ = 2 / 3 , C = 2 . The von Karman length-scale LvK is a three-dimensional generalization of the classic boundary layer definition LvK = κ

U′ U ′′

(4)

where U ′ = 2Sij Sij , Sij = 0.5(∂U i / ∂x j + ∂U j / ∂xi ) , U ′′ = ∂ 2U i / ∂xk2 ⋅ ∂ 2U i / ∂x 2j and κ = 0.41 is the von

Karman constant. In the software ANSYS Fluent, a limiter is added to compute LvK , which is proportional to the mesh cell size ∆ , i.e. the cubic root of the control volume size ΩCV

p is the pressure, ν is the kinematic viscosity, Fb represents the force acting on the unit mass and it is generally related to the gravitational acceleration g . A hybrid RANS/LES turbulence model, the scale adaptive simulation (SAS) model, is employed. The SAS model is based on the SST model, which behaves very much like the k - ω model away from the wall and serves to control the turbulence level in the near wall region. But the length scale produced by the SST model is too large. To avoid this shortcoming, in the SAS model, the von Karman length scale is introduced into the turbulence scale equation, such that the turbulence length-scale can be resolved correctly. The transport equations for the SAS model are defined as[8]

The porous media are modeled by adding one source term to the momentum equations. The source term is composed of two parts: a viscous loss term and an inertial loss term. For homogeneous porous media, the source term is defined as

∂ρ k ∂ + ( ρρ ui k ) = Gk − cµ k ω + ∂t ∂xi

1 µ  Si = −  vi + C2 ρ v vi  2 α 

∂ ∂x j

α

ω k

 µt  µ + σ k 

 ∂k  ,   ∂x j 

Gk − ρβ ω 2 + QSAS +

(1 − F1 )

2 ρ 1 ∂k ∂ω σ ω , 2 ω ∂x j ∂x j

∂ρω ∂ + ( ρ uiω ) = ∂t ∂xi ∂ ∂xi

 µt  ∂ω   µ + +  σ ω  ∂x j  

(2)

The additional source term QSAS is defined as  L − QSAS = max  ρη2k S 2 Lvk 

LvK

    ′ κη2 U 1/ 3 = max  κ , Cs ⋅ ∆  , ∆ = ΩCV  U ′′  β −a   cm  

(5)

(6)

where α is the permeability and C2 is the inertial resistance factor. This momentum sink contributes to the pressure gradient in the porous cell. In the convectional porous model, the inertial resistance factor is assumed to be constant. However, it is in conflict with the real situation, especially, when the size of the porous region is large. In this study, two assumptions are made in the porous region of the cylinder array: the local inertial resistance is proportional to the local pressure drop, the pressure head loss along the flow direction can be neglected. With these assumptions, the local inertial resistance can be redefined as C2 i =

Ptotal − Pi C2 , Ptotal

170

C2 =

F  Ptotal − Pi 1  Apatch ⋅ ρ vi vxi Si ⋅ ∆ ni   2 ∑ Ai  Ptotal 

(7)

the multi-body and the porous models in this study.

∑  i

where C2i is the inertial resistance coefficient in cell i , and C2 is the averaged coefficient value in the area

of the circular patches. Ptotal is the total pressure on the inlet boundary, Pi is the local cell static pressure, F is the time-averaged drag force of the circular patches, vi is the velocity magnitude, vxi is the velo-

Fig.1 The schematic diagram of computational model

city in the flow direction, Si is the cross sectional area of the local cell normal to the flow direction, Apatch is the area of the circular patches, ∑ Ai is the

The inlet velocity U ∞ is given as 0.098 m/s, and the hydraulic diameter H D = 4S wet / Pwet is used to

summation of the areas of the local cells, and ∆ni is

sectional wet area, and Pwet is the cross sectional wet

the thickness of the local cell. Ptotal is a constant when the inlet velocity remains the same. The main idea of the method is that the total drag force of the circular patches is equal to the summation of the forces on the local cells, equal to the pressure drops multiplied by the cross sectional areas, and the pressure drop can be calculated by the inertial resistance. In order to obtain the local inertial resistance coefficient at each cell, the multi-body model should be adopted first. The flow field variables in Eq.(7) are extracted at the first step calculation of the multi-body model. In this study, an equivalent diameter d =

perimeter. The turbulence intensity I = 0.16( ReDH )−1/ 8 .

determine the inlet turbulence, where S wet is the cross

ReDH is the Reynolds number defined by the hydrau-

lic diameter. The size of the computation mesh used for the multi-body model is 0.12×106, which is about 2 times larger than that for the porous model.

2 Ai / π is employed to calculate the values of Si and ∆ ni , assuming that the local cell as isotropic. In order to validate the feasibility and the rationality of the improved porous method with the adoption of the non-constant local inertial resistance coefficient, the calculated velocity field is compared to the measurements of the experiments conducted by Zong et al.[9]. The computational domain is shown in Fig.1, where the circular vegetation zone with diameter D = 0.22 m is centered at x = 0 . The diameter of each circular dot (vegetation unit) is d = 0.006 m . The flume’s width is 1.2 m, and the water depth is 0.133 m. Additionally, two monitoring lines (marked in red in Fig.1) are set at the appropriate positions to obtain the relative velocity profiles measured in the experiments. The time-averaged horizontal velocity u is measured at line y = 0 , while the time-averaged vertical velocity v is measured at line y = D / 2 , respectively.

The dimensionless constant Φ = N (d / D) 2 (where N represents the number of circular patches) is defined as the ratio of the sum of all circular patch areas to the total vegetation region area, representing the plant density of the vegetation. In the porous model, the value of Φ represents the porosity, and Φ = 0.03 both for

Fig.2 (Color online) Contour maps of the local inertial resistance coefficients

Figure 2 is the contour maps of the local inertial resistance coefficients. It is calculated based on the SAS multi-body model. As is described above, the porous model has a smaller size, that requires an interpolation of the C2 distribution inside the porous zone. Figure 3 shows the time-averaged horizontal velocity profile along the line y = 0 . Similarly, Fig.4 shows the time-averaged vertical velocity profile along the line y = D / 2 . U ∞ is the inlet velocity, D is the diameter of the circular patches. Both results of the multi-body model and the porous model fit well with the experimental data outside the circular region of the cylinder array. The multi-body model performs better in the prediction of the vertical velocity than the porous model away from the cylinders. Moreover, the multi-body model predicts the fluctuation of the flow velocity inside the region of the cylinder array, caused by impacting effect of the flow through the dense array of cylinders. Note that the distribution of the

171

cylinders is unknown in the experiment, so the measured velocity in the region of the cylinder array is not shown in Figs.3, 4.

Fig.6 (Color online) Contours of transient vorticity magnitude near porous zone

Fig.3 (Color online) Time-averaged horizontal velocity profiles along y = 0

Compared to the multi-body method, the improved porous model can give a proper result over the flow field while the simulation cost is much less than the multi-body model, such that it is suitable for the large-scale engineering applications. References

Fig.4 (Color online) Time-averaged vertical velocity profiles along y = D / 2

The vorticity fields simulated by the multi-body and porous models are shown in Figs.5, 6, where Fig.6 shows the details of the vorticity field around the cylinder array. As expected, the porous model can not capture the small scale vortices in the near-field, but still gives some details of the disturbance inside the porous zone, as shown in Fig.6. However, in the farfield, the porous model gives a similar vorticity field as the multi-body model.

Fig.5 (Color online) Vorticity fields of multi-body (a) and porous (b) models

[1] Ghisalberti M., Nepf H. M. Mixing layers and coherent structures in vegetated aquatic flows [J]. Journal of Geophysical Research: Oceans, 2002, 107(2): 1-11. [2] Nezu I., Onitsuka K. Turbulent structures in partly vegetated open-channel flows with LDA and PIV measurements [J]. Journal of Hydraulic Research, 2001, 39(6): 629642. [3] Yu L. H., Zhan J. M., Li Y. S. Numerical investigation of drag force on flow through circular array of cylinders [J]. Journal of Hydrodynamics, 2013, 25(3): 330-338. [4] Zhang H., YANG J. M., XIAO L. F. et al. Large-eddy simulation of the flow past both finite and infinite circular cylinders at Re = 3 900 [J]. Journal of Hydrodynamics, 2015, 27(2): 195-203. [5] Su X., Li C. W. Large eddy simulation of free surface turbulent flow in partly vegetated open channels [J]. International Journal for Numerical Methods in Fluids, 2002, 39(10): 919-937. [6] Yu L. H., Zhan J. M., Li Y. S. Numerical simulation of flow through circular array of cylinders using multi-body and porous models [J]. Coastal Engineering Journal, 2014, 56(3): 1450014. [7] Wang Q., Li M., Xu S. Experimental study on vortex induced vibration (VIV) of a wide-D-section cylinder in a cross flow [J]. Theoretical and Applied Mechanics Letters, 2015, 5(1): 39-44. [8] Menter F. R., Egorov Y. The scale-adaptive simulation method for unsteady turbulent flow predictions. Part 1: Theory and model description [J]. Flow, Turbulence and Combustion, 2010, 85(1): 113-138. [9] Zong L., Nepf H. Vortex development behind a finite porous obstruction in a channel [J]. Journal of Fluid Mechanics, 2012, 691: 368-391.