M. I. Alamayreh, A. Fenocchi, G. Petaccia, S. Sibilla, E. Persi

1. Alternative Energy Department, Al-Zaytoonah University of Jordan, Jordan, Amman

2. Department of Civil Engineering and Architecture, University of Pavia, Pavia, Italy

Abstract: The drag and side force coefficients of a half-submerged cylinder in a free-surface flow were calculated through numerical simulations, with the aim of supporting the numerical modelling of log transport in rivers. The variability of these coefficients with the yaw angle with respect to the flow direction and with the ratio between the flow depth and the diameter of the cylinder were investigated. Simulations were performed with the three-dimensional code ANSYS/CFX, employing the volume of fluid multiphase technique to reproduce the critical interaction between the free surface and the cylinder. The numerical tests, showing the rise of the drag force coefficient for increasing yaw angles passing from flow-parallel to flow-perpendicular cylinder and the peak of the side force coefficient for flow-oblique cylinder, were validated by comparison with the results of laboratory experiments. The simulations were then extended to conditions with significant blockage in the vertical direction which had not been previously experimented, revealing a strong increase in the force coefficients for decreasing ratios between the flow depth and the cylinder diameter. A detailed description of the reproduced flow features in the proximity of the cylinder for the different cases was furthermore obtained. Such report, in addition to the analysis of the force coefficients, could serve a much wider range than that of log transport, i.e., any case in which a floating cylinder interacts with free-surface flow.

Key words: Large woody debris (LWD) transport, floating objects, drag and side force coefficients, vertical blockage, free-surface interaction

Introduction

Large woody debris (LWD) transport is a serious risk during floods, as it can cause the clogging of river sections, especially where inline structures such as bridges are located[1-3]. Numerical models of the hydrodynamic transport of floating logs[4-5]can help developing strategies to reduce such risk, by identifying deposition-prone areas[6-7]. To properly model LWD transport, the hydrodynamic forces acting on logs must be correctly evaluated. An appropriate description of the flow around each log under all conditions is hence needed.

This work is part of an experimental and numerical effort[5,8]to characterise the open-channel flow dynamics around a half-submerged yawed cylinder,considered representative of a real LWD[6-7], aimed at supporting two-dimensional numerical modelling of LWD transport[5]. Though some works[9-10]attempted to simulate the fluid flow dynamics around trees considering their real complex geometry, a complete description of the variation of the hydrodynamic force coefficients for simple floating cylinders with the yaw angle is not available yet. Such knowledge is essential for a numerical model of LWD transport, as it allows proper computation of forces and torques acting on the logs, which determine their path and orientation.Furthermore, during the initial and final phases of flood events, the ratio between the flow depth and the diameter of logs is small. High blockage in the vertical direction may therefore be involved, causing relatively stronger hydrodynamic forces to be exerted on the logs than those acting when blockage is negligible.

The force coefficients of cylinders depend on many parameters such as their relative density with respect to water, yaw angle, aspect ratio, submergence and blockage ratio. Past literature on force coefficients focused on fully submerged cylinders with negligible blockage, overlooking the partially submerged case,which is typical of log transport[8], and therefore the relative subcase with shallow flow relative to the log diameter which is typical of the incipient motion and arrest of logs. In the case of partially submerged cylinders, the interaction between the free surface and the body plays a key role, vortex shedding and flow separation occurring behind the cylinder[10]. Chu et al.[11]determined through numerical simulations the values of the drag and lift coefficients of a fully submerged horizontal cylinder close to the water surface and perpendicular to the flow for different cylinder Froude numbers, with U being the mean undisturbed flow velocity, g the gravity acceleration and D the cylinder diameter)and low submergence. Chu et al.[11]showed that the amplitude of the water surface alterations increases with the cylinder Froude number, a hydraulic jump occurring downstream of the cylinder forregardless of the submergence within their tested range, due to the obstruction caused by the body.When the cylinder is large compared with the flow depth, the hydrodynamics are further complicated by flow constriction below the body, which also interacts with the phenomena occurring at the surface.

An accurate evaluation of the force coefficients around a yawed cylinder can be obtained either experimentally or through computational fluid dynamics (CFD). The latter provides an insight into flow dynamics which complements experiments and helps validating their results. Various numerical methods have been adopted in literature to simulate flow around cylinders to determine the force coefficients,ranging from Eulerian, mesh-based methods as finite volume (FV)[11-12]or finite difference (FD)[13], to Lagrangian, meshless methods as smoothed particle hydrodynamics (SPH)[14-15]or discrete vortex(DV)[16-17]. Some of the mentioned studies have been devoted to the analysis of laminar flow[14-15], while others have addressed the turbulent regime, either adopting the Reynolds-averaged Navier-Stokes(RANS) equations with a turbulence model[12]or performing large eddy simulation (LES)[13].

The aims of this work are:

(1) To improve the knowledge on the flow dynamics around a half-submerged horizontal cylinder in an open channel, including the case in which the ratio between the flow depth and the cylinder diameter is small, through numerical FV solution of the incompressible RANS equations under steady conditions.

(2) To validate through numerical simulations the laws for the variation of the drag and side force coefficients with the yaw angle obtained from the laboratory experiments in Persi et al.[8], which were performed under negligible blockage, further extending them to the shallow flow case.

In this paper, the numerical background of the performed simulations, the tested conditions and the methodology employed to evaluate the force coefficients are first presented, the simulated flow behaviours around the cylinders and the computed force coefficients for the different cases being then reported and analysed.

1. Numerical methods

1.1 Simulation of the flow field

The numerical simulations of flow around a half-submerged cylinder in an open channel were performed with the three-dimensional FV RANS code ANSYS/CFX[18]. Both water and air phases were considered to solve the free surface through the volume of fluid (VOF) multiphase Eulerian-Eulerian approach[19-20]. Using the l and g subscripts to denote the liquid and gaseous phases respectively, the continuity and momentum equations for the phases are:

where u=(ux, uy, uz) is the velocity vector, ρ is the density, p is the pressure, μ is the kinematic viscosity and e is the rate of deformation tensor. If n is the unit normal vector pointing towards the positive curvature region, the local curvature of the water-air interface is κ= ∇ ·n and the dynamic condition expressing the balance of surface forces at the interface can be written as

where σ is the surface tension between water and air, set to 0.0728 N/m[21]. The kinematic interface condition implies that at the interface. Water and air fractions coexist at mesh elements in the proximity of the free surface due to the interference between phases. The evolution of the water fraction α is described by

The computational domain is 2.70 m long, 0.49 m wide and 0.40 m high, matching the dimensions of the laboratory flume in the experiments by Persi et al.[8],whose results are taken as comparison. The same flow conditions with undisturbed flow depthh =0.15 m and discharge Q =29.5l/s were reproduced.

The geometry (Fig. 1) and the numerical meshes(Fig. 2) were created within the ANSYS Workbench and Mesh packages, respectively. Each mesh is made up of approximately 6×105tetrahedral elements,refined in the proximity of the cylinder to more accurately reproduce the flow close to it and across the average free-surface elevation throughout the domain to determine a sharper interface (Fig. 2).Preliminary tests were performed to check that the selected mesh resolution ensured mesh-independence of the solution, especially as regards the force coefficients.

Fig. 1 (Color online) Schematics of the computational domain with the adopted boundary conditions

Fig. 2 (Color online) Slice view of a tetrahedral grid used in the numerical simulations

Within the ANSYS/CFX flow solver, the convective terms are discretised with a 2nd-order-accurate bounded scheme and the equations are integrated in time with an implicit 2nd-order-accurate scheme. The FV method is coupled with a finite element (FE)strategy to describe the variation of flow variables within each mesh element[18]. The solver adopts a multigrid (MG) approach to reduce long-wavelength errors in the solution of the coupled linear set of algebraic equations resulting from the discretisation by performing the initial iterations over a coarser mesh, accelerating and improving the convergence of the system[18]. Turbulence was modelled using a shear stress transport (SST) k -ω two-equation model[22].

Simulations were initially performed to reproduce the experiments in Persi et al.[8], in which a half-submerged cylinder was considered. Exact half submergence, i.e., with the cylinder axis at the undisturbed water surface elevation, was selected to simplify the positioning operations in the flume. Log submergence under real conditions slightly varies according to wood density and water content, even though we do not expect such small variation to significantly affect the force coefficients. The length of the cylinder was L=0.15 m and its diameter was D=0.025 m, hence the aspect ratio and the flow depth to cylinder diameter ratio are L/ D =6 and h/ D = 6, respectively. The blockage ratio, i.e.,LD over the undisturbed wetted area A=0.0735 m2, is about 5% for the experimental campaign, which should not influence the hydrodynamic coefficients[23-24].Yaw angles within the range 0°＜ θ ＜90° with 10°increments were tested, in addition to θ = 45°.

The influence of the ratio between the flow depth and the cylinder diameter, not assessed in the laboratory, was then investigated numerically testing half-submerged cylinders with diameters D=0.15 m and D=0.05 m. As the length of the cylinder is still L=0.15 m , the aspect ratios are L/ D =1 and L/ D =3, respectively. The ratios of flow depth to cylinder diameter are h/ D = 1 and h/ D=3.Under such conditions, the vertical blockage, defined as the reciprocal of the h/ D ratio, would become the most significant parameter for the hydrodynamic coefficients, as the flow deviation should affect both the flow above and below the cylinder, increasing the influence of the free-surface and bottom interactions.

Proper boundary conditions were set (Fig. 1) to reproduce the laboratory experiments[8]. A uniform normal water velocity U = Q / A =0.401m/s was set at the water inlet boundary, while atmospheric pressure was assumed at the water outlet boundary.Unit water volume fraction was set at both boundaries.Atmospheric pressure and unit air volume fraction were set at the three air boundaries. No-slip wall conditions were set at the cylinder boundary and at flume bottom and side walls. The complete geometric and hydraulic parameters of the study are collected in Table 1.

Table 1 Geometric and hydraulic parameters of the study

Simulations were performed starting from waterat-rest initial conditions with h =0.15 m up to the achievement of the steady state. Convergence of solutions was ensured by requiring the residuals of the discretised equations to be lower than 10-4times the residuals at the first iteration across the whole computational domain. Each run was performed in parallelusing an Intel® Core i7-7700 CPU with 8 threads,resulting in a computational time of ～45 min per simulation.

1.2 Evaluation of the force coefficients

To calculate the drag and side forces on the cylinder in the simulations, the momentum conservation equation was applied to a suitable control volume enclosing the body. We describe the approach, which was adopted e.g., by Allen and Smith[9], for the drag force coefficientCD(Fig. 3), related to the horizontal component of the hydrodynamic force in the main flow direction, the extension to the side force coefficientCS, relative to the force normal to the main flow, being straightforward. The streamwise momentum conservation equation in the control volume of a channel enclosing the partially submerged cylinder can be written, considering its upstream and downstream cross sections 1 and 2 respectively, as

Fig. 3 Scheme of the momentum conservation approach for the calculation of the drag force on the cylinder considering the control volume enclosing a partially submerged cylinder for θ = 90° (side view)

wherepgis the pressure at the centre of gravity of the cross section,Acis the actual wetted area of the cross section,FDis the drag force on the cylinder andβis the Boussinesq coefficient, which corrects the mean undisturbed velocityUto contemplate the departure from uniform flow conditions owing to boundary and wake effects. For the calculations, the upstream cross section was placed 1 0Dupstream of the centre of mass of the cylinder, while for the downstream one, influenced by wake effects, a position 20Ddownstream of the centre of mass of the cylinder was determined as optimal through preliminary simulations, the resulting force coefficients being the closest to the experimental ones by Persi et al.[8]. The general expression of theβcoefficient for theicross section is

whereuis the local velocity component in the main flow direction. The Boussinesq coefficients were herein calculated from the numerical solution. If theicross section is discretised in the numerical model byNicell surfaces,βcan be computed as

whereujandAjare the local velocity component and area in the main flow direction for cell surfacej.Reformulating Eq. (8), the drag force can be expressed as

The drag force coefficient is eventually obtained as

The side force coefficient for the side forceFSis similarly

2. Results and discussion

2.1 Flow field analysis for different yaw angles and ratios between the flow depth and the cylinder diameter

The observed behaviour of the drag and side force coefficients depending on the geometrical parameters can be better understood analysing first the simulated flow field around the cylinder. Figs. 4, 5 display the velocity, static pressure and vorticity fields in the proximity of the half-submerged cylinder, both in the vertical plane at the centreline of the channel and in the horizontal plane at the undisturbed water surface elevation for h/ D = 6 and θ= 0°, 45°,90°. It can be seen from the planar views of the velocity fields (Fig. 4) that recirculation occurs in the wake behind the cylinder for θ= 90°, while for θ = 45°water flows along the yawed cylinder. This is confirmed by the vorticity fields (Fig. 5), which display quite high values in the wake regions, albeit much lower than those in the boundary layer of the cylinder, null velocity occurring at the cylinder boundary owing to the no-slip condition. The static pressure (Fig. 5) is the highest in front of the cylinder for θ = 90° due to stagnation, which determines a rise in the local water level. Complete stagnation does not occur for θ= 45°, the upstream pressure and water level rise being consistently lower than those in the previous case. F or θ = 0° , water flows along the cylinder and only a thin wake is present downstream of it, entraining high vorticity arising from the boundary layer along the sides. Stagnation forms upstream of the cylinder, while a local pressure minimum appears in the wake immediately downstream of the cylinder, leading to a local drop of the free-surface elevation.

Fig. 4 (Color online) Simulated velocity magnitude and streamlines in the vertical plane at the centreline of the channel and in the horizontal plane at the undisturbed water surface elevation for h / D = 6 and θ = 0°, 45°, 90°

Fig. 5 (Color online) Simulated vorticity in the horizontal plane at the undisturbed water surface elevation and static pressure in the vertical plane at the centreline of the channel for h / D = 6 and θ = 0°, 45°, 90°

Figures 6-8 portray the water volume fraction,velocity and static pressure fields in the vertical plane at the centreline of the channel for h / D =1, 3, 6 and θ=45°, 90°. The behaviour of the free surface (Fig.6) in the proximity of the floating cylinder strongly depends on the ratio between the flow depth and the cylinder diameter and on the yaw angle. Maximum superelevation results for θ = 90°. The ratio between the superelevation and the cylinder diameter increases with h/ D . The highest velocities (Fig. 7) occur just below the cylinder due to flow contraction, the maximum values being attained for h / D =1, owing to the vertical blockage of the cylinder caused by its smaller distance from the bottom. Complete stagnation (Fig.8) occurs only for θ = 90°, at a lower elevation than the undisturbed free surface, due to the deformation of the latter with water flowing also above the cylinder.The static pressure distribution along the cylinder is hence asymmetric, with low-pressure regions below and above it. This pressure field structure in the vertical plane is substantially independent of the relative shallowness of the flow, though for lower h/ D values the pressure minimum below the cylinder is enhanced by bottom proximity. The local pressure minimum just downstream of the cylinder for h/ D = 3and h / D =6 is caused by vortices in the wake, being absent for h/ D =1 as in that case the boundary layer does not separate. Separation is likely prevented for h/ D =1 by the downward deflection of the flow caused by the low-pressure region below the cylinder. The θ= 45° cases display lower maximum static pressure values (Fig. 8) upstream of the cylinder, with only partial stagnation occurring.The missing separated wake is related to the oblique flow along the cylinder and strongly reduces the downstream pressure minimum compared to θ = 90°.

Fig. 6 (Color online) Simulated water volume fraction in the vertical plane at the centreline of the channel for h / D =1, 3, 6 and θ=45°, 90°

Fig. 7 (Color online) Simulated velocity magnitude and streamlines in the vertical plane at the centreline of the channel forh/ D =1, 3, 6 and θ =45°, 90°

Fig. 8 (Color online) Simulated static pressure in the vertical plane at the centreline of the channel for h / D = 1, 3, 6 and θ=45°, 90°

2.2 Computed drag and side force coefficients

The computed values of CDand CSas function of the yaw angle for h / D =6 are shown in Fig.9, compared with the correspondent experimental values by Persi et al.[8]. Successful agreement is found between simulations and experiments. The relative mean absolute errors of the numerically found values against the laboratory ones over the tested yaw angles are 4.03% for CDand 2.89% for CS, while the maximum discrepancies are 8.46% and 6.88%,respectively. However, the force coefficients obtained from the simulations are generally smaller than the experimentally measured ones. This is because the model cannot fully reproduce the complex local turbulence dissipation phenomena, especially those ascribed to the free-surface interaction near the cylinder, which cause additional force to be exerted on the body. The CDcurve as function of the yaw angle(Fig. 9(a)) displays: (1) an almost constant CD≈0.2,equal to the known experimental value of the drag force coefficient for a cylinder aligned with the stream direction in the tested flow regime[25], for moderate yaw angles θ ≤ 1 0°, (2) a monotonic, almost linear increase for higher yaw angles, (3) a plateau at maximum drag for 80°≤ θ ≤ 9 0°. T he values of CSas function of the yaw angle (Fig. 9(b)) show instead an almost symmetric behaviour, peaking around θ= 45°.

Fig. 9 (Color online) Comparison between simulated and experimental[8] force coefficients for h/ D =6 as function of the yaw angle

The impact of the ratio between the flow depth and the cylinder diameter h / D on the values of the drag and side force coefficients is presented in Fig. 10.Both force coefficients increase more than linearly for decreasing h/ D ratios, the difference between h/ D = 3and h/ D=1 conditions being strikingly relevant. Such behaviour can be certainly ascribed to the relevant vertical blockage for h/ D =1, the computed values being even larger than those predicted by the empirical relation by Gippel et al.[26].The latter expression was however based on experiments on fully submerged cylinders and hence does not consider the effect of the free surface.

Fig. 10 (Color online) Simulated force coefficients for h/ D =1, 3, 6 as function of the yaw angle

3. Conclusion

Force coefficients for a yawed half-submerged cylinder inside an open channel with different ratios between the flow depth and the cylinder diameter were computed through FV numerical simulations and compared against laboratory experiments. The model allowed describing the flow dynamics under the tested conditions, improving the knowledge on the interaction between woody debris and water flows needed for the development of numerical models of log transport in rivers. The drag and side force coefficients calculated from the simulations well matched the experiments, yet being mostly smaller than the laboratory ones due to the inability of the model to fully capture turbulent dissipation processes at the free surface near the cylinder. Drag force coefficients peak at θ = 90° due to maximum flow obstruction and diminish with decreasing yaw angles following an S-like curve. Side force coefficients culminate for θ = 45°, the largest flow deviation being attained, and reduce around the maximum following an almost symmetric bell-shaped curve. Force coefficients increase with a strongly more-than- linear behaviour for decreasing ratios between the flow depth and the cylinder diameter, so that the relative shallowness of the flow should be definitely taken into account in numerical models of log transport to properly describe motion during the initial and final phases of flood events.

Acknowledgements

The authors greatly appreciate the financial support received by Dr. M. I. Alamayreh from the Centre for International Cooperation and Development of the University of Pavia (CICOPS) to carry out his stay as visiting researcher.