Numerical Simulation on Flow Past Two Side-by-Side Inclined Circular Cylinders at Low Reynolds Number

A series of three-dimensional numerical simulations is carried out to investigate the effect of inclined angle on flow behavior behind two side-by-side inclined cylinders at low Reynolds number Re=100 and small spacing ratio T/D=1.5 (T is the center-to-center distance between two side-by-side cylinders, D is the diameter of cylinder). The instantaneous and time-averaged flow fields, force coefficients and Strouhal numbers are analyzed. Special attention is focused on the axial flow characteristics with variation of the inclined angle. The results show that the inclined angle has a significant effect on the gap flow behaviors behind two inclined cylinders. The vortex shedding behind two cylinders is suppressed with the increase of the inclined angle as well as the flip-flop gap flow. Moreover, the mean drag coefficient, root-mean-square lift coefficient and Strouhal numbers decrease monotonously with the increase of the inclined angle, which follows the independent principle at small inclined angles.


Introduction
Multiple cylinders in side-by-side configuration are widely applied in offshore and ocean engineering, such as piers, bridge pilings, risers and pipelines. The complicated wake flow characteristics behind side-by-side cylinders have attracted the attention of many researchers, especially for the gap flow pattern at small spacing ratios.
In the past several decades, comprehensive experimental studies have been conducted to investigate the flow characteristics behind two side-by-side cylinders, especially for the flip-flop flow (Ishigai et al., 1972;Bearman and Wadcock, 1973;Zdravkovich, 1977;Williamson, 1985;Kim and Durbin, 1988). Zdravkovich (1977) observed the flip-flop gap flow at small spacing ratios of 1.2<T/D<2.5 (T is the center-to-center distance between two side-by-side cylinders, D is the diameter of the cylinder). The period of the flip-flop flow is found to be much larger than the vortex shedding period (Williamson, 1985), the non-dimensional vortex shedding frequency and drag coefficient of the cylinder with narrow wakes are larger than those with wide wakes (Sumner et al., 1999;Xu et al., 2003). Wang and Zhou (2005) pointed out that the interactions between dif-ferent-signed vortices lead to the gap flow changeover phenomenon. In recent years, many numerical simulations have been performed on flow past side-by-side cylinders. Meneghini et al. (2001) conducted a two-dimensional numerical simulation on flow past side-by-side cylinders at Re=100-200 using a fractional step method, and elaborated the reason for the repulsive lift force of two cylinders by the pressure contours. Six types of wake patterns behind two side-by-side cylinders are observed by Kang (2003), and the period of flip-flop flow is found to be one order larger than that of the vortex shedding. Afgan et al. (2011) found that the vortex shedding pattern at 1.25≤T/D≤1.75 is dominated by the out-of-phase mode, and would switch to be inphase mode when the changeover phenomenon occurs. Similar conclusion has also be drawn by Alam and Zhou (2013) and Carini et al. (2014), where the phase shift of the vortex shedding is considered to trigger the changeover phenomenon of gap flow.
However, in the practical offshore and ocean engineering applications, the current flow direction is not always perpendicular to the axial of ocean structures, with the existence of an inclined angle between the flow direction and the cylinder axis (Zhao et al., 2009). Many experimental and numerical studies have been performed on flow past a single inclined cylinder, which concentrate on two aspects: the effect of inclined cylinder on the axial flow field and the validation of independence principle (IP). Shirakasi et al. (1986) observed the vortex motion along the axial direction at the back of the inclined cylinder, and found that the larger the inclined angle, the lower the vortex shedding frequency. Later, Li and Gu (2005) analyzed the power spectral density (PSD) of the fluctuating lift coefficients and proposed that the peak values not only existed at the Karman vortex shedding frequency but also the low-frequency components. Hogan and Hall(2010) found that the downstream vortex street will be disturbed if the inclined angle is large enough (α>30°). In recent years, Shao et al. (2013) have found that the inclined orientation has a significant effect on the force coefficients of cylinder. Franzini et al. (2014) reported that there is no regular vortex shedding regime when the inclined angle α=45°. Thapa et al. (2014) stressed that the inclined cylinder would suppress the vortex shedding. Liang et al. (2015) summarized the different types of vortex shedding patterns at different spanwise sections. Hu et al. (2017) further investigated flow past an inclined square cylinder by the PIV technique, and observed the axial flow significantly alters the aerodynamic instabilities of the inclined cylinder. Zhao et al. (2009) stated that the mean drag coefficients and Strouhal number follow the IP, but the root-mean-square lift coefficients (r.m.s lift coefficients) are obviously depended on the inclined angle. Thapa et al. (2015) pointed that the normal Reynolds number should be used in the implementation of the IP.
Although extensive numerical and experimental literatures have investigated flow past two side-by-side vertical cylinders or a single inclined cylinder, there remain some questions to be answered, such as (1) what kind of flow behaviors will be observed when two inclined cylinders are placed in a closely side-by-side configuration? (2) How will the inclined angle influence the force coefficients of two cylinders? The motivation of this three-dimensional numerical study is to answer the above questions for flow past two side-by-side inclined cylinders at a small spacing ratio and low Reynolds number, and special attention is paid on the gap flow and axial flow characteristics.

Governing equations
In this numerical study, the incompressible viscous Navier-Stokes equations are used for the governing equations of flow past two side-by-side inclined cylinders. In the Cartesian coordinate system, the continuity and momentum equations are expressed as: ∇ · U = 0; (1) where U is the velocity of fluid, t is the time, ρ is the density of fluid, p is the pressure, and v is the kinematic viscosity of the fluid. The Navier-Stokes equations are solved by the Open source Field Operation And Manipulation (Open-FOAM) C++ libraries.
The laminar standard solver icoFoam in OpenFOAM is used for the pressure velocity coupling. The time scheme is discretized by the first-order implicit. The Gauss liner scheme is used for the gradient and divergence terms, which has the second-order accuracy. The second-order Gauss liner orthogonal scheme is applied to discretize the Laplacian and pressure terms. Fig. 1 shows the computational domain for the flow past two side-by-side inclined cylinders. The upper and lower cylinders are labeled as C1 and C2, respectively. The length of the computational domain is 30D in the streamwise direction, and the width of the domain in the transverse direction is 21.5D. The distance between the inlet boundary and center of the cylinder is 10D, similar to the cases of Afgan et al. (2007) and Parnaudeau et al. (2008). The height of the cylinder is set as Z=4D (Lei et al., 2001) to ensure the accuracy of the three-dimensional wake flow. The spacing ratio between two side-by-side cylinders is kept to be T=1.5D. The definition of the inclined angle is shown in Fig. 1b, and u n and u t are the orthogonal components of the flow velocity. The range of the inclined angle α is varied from 0° to 60° with an increment of 7.5°. LIU Cai et al. China Ocean Eng., 2019, Vol. 33, No. 3, P. 344-355 345 Similar to the study of Tong et al. (2014), the initial values for the velocity and pressure in the whole domain are set to zero. At the inlet boundary, a uniform streamwise (x-direction) velocity u 0 =0.01 m/s is imposed, the transverse (ydirection) and spanwise (z-direction) velocities are set to v=w=0. The outflow condition is applied on the outlet condition, where both the pressure and velocity gradients are zero. Symmetry boundary condition is applied on two later-al sides, where the normal velocity and pressure gradients are equal to zero. Due to the significant effect of the spanwise velocity for flow past two inclined cylinders, the zero gradient condition is applied at the top and bottom boundaries. The cylinder surfaces are applied as no-slip boundary condition. The details of the boundary conditions are summarized in Table 1. Fig. 2 shows the computational meshes for flow past two side-by-side inclined cylinders at α=30°. The element number on the circumference of each cylinder is 100 and 50 layers of elements are distributed along the spanwise direction, which leads to the element size Δz=0.08D. The minimum mesh size in the O-grid region of the radial direction ∆y ∆y ∆t ∆t close to the cylinder surface is evaluated by y + =u f /v (u f is the friction velocity). In this research, is set to be 0.09D for all the cases and the distance between the nodes next to the wall is in the range of y + <1. The non-dimensional time step is u 0 /D=0.01 ( is the time step) to ensure the Courant-Friedrichs Lewy CFL=u 0 Δt/Δx (Δx is the distance between two adjacent grids).  Tables 2 and 3, respectively. The mean drag coefficient , root-mean-square lift coefficient and Strouhal number St are defined as:

Mesh and time step dependency study ∆y
(3) where is the total mean drag force in the x-direction, is the root-mean-square lift force in the y-direction, H is the cylinder length, f s is the vortex shedding frequency, and is the pressure of the incoming flow. Table 3, the discrepancies in the mean drag coefficient , r.m.s lift coefficient and Strouhal number among the results obtained from three different meshes   Zhou et al. (1999) and Ding et al. (2007). The Strouhal numbers for Mesh 2 and Mesh 3 are 0.163, which is also in agreement with the previous studies. Therefore, the accuracy of the present simulation results is considered to be acceptable using the medium mesh density.  Table 4. It can be seen that the difference of the mean drag coefficient obtained with u 0 Δt/D=0.02 and 0.01 is about 2.88% and it decreases to smaller than 0.5% obtained with u 0 Δt/D=0.005. Moreover, the difference of the r.m.s lift coefficient obtained with three time steps decreases to 0.93%. Therefore, a non-dimensional time step of u 0 Δt/D=0.01 is adopted for all the simulations in the present study.

Results and discussion
3.1 Instantaneous flow field Fig. 3 shows the instantaneous contours of the non-di-mensional streamwise velocity u/u 0 at different inclined angles behind two side-by-side inclined cylinders, where the iso-surfaces of vortex flow wake are identified by the Q criterion (Hunt et al., 1988).
As shown in Fig. 3, for all cases (expect α=60°), the instantaneous contours behave paralleled with the spanwise direction of the cylinders regardless of the increase of the inclined angle, the iso-surfaces of vortex becomes inclined at large inclined angles. Besides, no deformation in the spanwise is observed, and indicates the weak effect of the three-dimensional flow at low Reynolds number. However, for the case of α=60°, there is no vortex shedding behind two cylinders. The flow filed is not fully developed and the vortex shedding is suppressed at large inclined angle. One point should be stressed that the minimum streamwise velocity u/u 0 <0 when α<30°, which indicates the existence of the recirculation zone. At α>30°, the minimum streamwise velocity u/u 0 >0, the recirculation zone would disappear.
In order to further investigate the detailed wake patterns, the non-dimensional instantaneous vorticity contours ω * =ω z D/u 0 (ω z is the vorticity in spanwise) at the middle spanwise section z/D=2.0 for typical cases are presented in Fig. 4. For the cases of α≤37.5°, the direction of gap flow switches from one side to another side randomly, that is to say, the flip-flop flow pattern occurs with irregular time intervals, an asymmetry vortex street with wide wake and narrow wake is observed, consistent with Ishigai et al. (1972), Bearman and Wadcock (1973) and Kim and Durbin (1988). Taking the case of α=22.5° for example, at the instant u 0 t/D=516, the gap flow is deflected toward to the upper cylinder C1, resulting in a narrow wake behind C1 and a wide wake behind C2. Due to the interactions between the different-signed vortices in narrow wakes, the wakes in the inner side of cylinders are approximate parallel with the x- axis at the instant u 0 t/D=542. However, it appears that the current flow pattern is unstable. The gap flow will deflect toward C2 at the instant u 0 t/D=572, resulting in a narrow wake behind C2 and a wide wake behind C1. With further increase of the inclined angle to α≥45°, the changeover phenomenon of gap flow is not significant compared with those at small inclined angles. At α=52.5°, the free shear layers on the lateral side are stretched with large scales and the vortex strength becomes weak. At large inclined angle α=60°, due to the effect of axial flows, no vortex shedding is observed behind the two side-by-side cylinders.
It can be seen that the inclined angle has a significant effect on the vortex shedding pattern behind two side-by-side cylinders. With the increase of the inclined angle, the flipflop flow will be suppressed, and gradually diminish. Fig. 5 presents the three-dimensional streamline topologies behind two side-by-side cylinders at different inclined angles. It can be seen from Fig. 5 that for the vertical case, all the streamline behave paralleled with the xoy plane, and there is no movement in the spanwise direction. When the inclined angle α increases to the range of 7.5° to 22.5°, one portion of the streamlines in the recirculation zone appear to  Fig. 6 shows the time-averaged streamline topologies with the normalized mean streamwise velocity contours on two streamwise sections y/D=±0.75 (through the axis of C1 and C2). The purple dotted line depicts the values of =0, indicating the division between the recirculation zone and outer zone. At α=0°, all the streamlines are approximately paralleled with the x-axis, and the contourline of mean streamwise velocity also behaves paralleled with the cylinders, in good agreement with the aforementioned. The widths of the recirculation zones behind two cylinders are different due to the bistable flow pattern, and no obvious axial flow is observed. When α increases to the range of 7.5° to 22.5°, the streamlines behind the cylinders behave more and more inclined, however, the recirculation zones always exist behind two cylinders. One point should be stated that whether the axial streamlines are in the recirculation zones, it always deflects to the upper side of the cylinder (downstream side). However, increasing the inclined angle to α≥30°, the axial streamlines appear to be inclined significantly and the recirculation zones disappear. In contrast to the previous works, for the inclined cases, the axial streamlines neither gather at downstream side (Hayashi and Kawamura, 1995) nor disturb the vortex street at large inclined angles (Hogan and Hall, 2010), and the reason may be contributed to the low Reynolds number. The free-stream velocity associated with the low Reynolds number (Re=100) is u 0 =0.01 m/s, resulting in very small axial velocity (much smaller than u 0 ) to disturb the vortex street at the downstream side at the large inclined angles (α=30°-52.5°).w /u 0 Fig. 7 shows the time-averaged streamline topologies with the normalized mean spanwise velocity on the middle spanwise section z/D=2. It is clear that, for the vertical case, there are multiple eddies with different scales behind the two cylinders. With the increase of the inclined angle, the scales of eddies decrease gradually and will disappear at α=30°, and the recirculation zones will also disappear simultaneously. However, when α≥30°, the streamline topology behaves like a pattern without separation, and roughly symmetrical about y/D=0. It is observed that the values of the axial velocities increase monotonously with the increase of the inclined angle except for the case of α>45°, where the axial velocity decreases slightly, and the maximum value of the mean spanwise velocity could reach to 40% of the streamwise velocity at α=45°. Fig. 8 shows the time histories of the lift coefficients C l and drag coefficient C d on two side-by-side cylinders at different inclined angles. i=1, 2 represent the force coefficients associated with C1 and C2, respectively.

Force coefficients
For the cases of α≤52.5°, the time histories of the lift coefficients C l behave as a "quasi-periodic" pattern. The magnitudes of the drag coefficients C d for two cylinders change alternatively with the process of flip-flop, while the magnitudes of the lift coefficients C l change frequently. For example, at α=30°, when 420<u 0 t/D<550, the gap flow deflects towards C2, leading to a narrow wake and a larger drag coefficient of C2 than that of C1, and the corresponding amplitudes of the lift coefficients within the aforementioned time duration decrease dramatically. When

350
LIU Cai et al. China Ocean Eng., 2019, Vol. 33, No. 3, P. 344-355 550<u 0 t/D<620, the direction of the gap flow switches from C2 to C1 randomly as shown in Fig. 4e. As a result, the amplitudes of the lift coefficients increase compared with that of 420<u 0 t/D<550. Similar phenomenon is also observed at the other inclined angles. It is found that when the gap flow is biased toward one cylinder with a long-time in-  LIU Cai et al. China Ocean Eng., 2019, Vol. 33, No. 3, P. 344-355 351 terval (deflected pattern), the amplitude of the lift coefficient will decrease dramatically. When the gap flow switches from one side to another side intermittently (flipflop pattern), the amplitude of the lift coefficient will increase gradually. The amplitude of the lift coefficient decreases with the increase of the inclined angles. However, the periodicity of the lift coefficient is more apparent at α>30°, which indicates that the flip-flop is gradually suppressed with the increase of inclined angles. It is interesting to point out that for the case of α=60°, both the time histories of the drag and lift coefficients behave like straight lines, demonstrating the disappearance of the vortex shedding phenomenon.
Moreover, the different-signed lift coefficients reveal that there is a repulsive force acting on the two cylinders.
In order to investigate the effect of the inclined angles on the force coefficients, the variations of the mean drag coefficient and r.m.s. lift coefficient with the inclined angle are shown in Fig. 9. Both and on the two cylinders decrease with the increase of the inclined angle α. Moreover, the amplitudes of and on two cylinders are almost the same, the deviation of between the two cylinders is smaller than 3%, and the largest deviation of is 3.31% at α=22.5°. Compared with the case of α=0°, decreases to 50.55% at the large inclined angle α=60°. Fig. 10 shows the circumferential mean pressure coeffi-

352
LIU Cai et al. China Ocean Eng., 2019, Vol. 33, No. 3, P. 344-355 Cp 0 Cp 2 cient distributions around two cylinders on the middle spanwise Section I-I (z/D=2) at different inclined angles. represents the mean pressure coefficient of a single vertical cylinder, and represent the mean pressure coefficients of the upper and lower cylinders C1 and C2, respectively. As shown in Fig. 10, the high-pressure regions move toward the gap side, as a result, the maximum values of and appear at θ=15° and 345° (-15°) instead of 0° for each case. Consequently, the pressure near the gap is larger than those at other regions, and the repulsive lift coefficient is formed due to the pressure discrepancy. It is observed that and are roughly symmetrical around the center of θ=180° for all cases, implying that the separated points of free shear layers of two cylinders are roughly symmetrical about the x-axis. It can be qualitatively deduced from Fig. 10 that the pressure discrepancy between the base point (θ=180°) and standing point (θ=0°) decreases with the increase of the inclined angle.
The FFT (Fast Fourier Transformation) analysis of the time series of the fluctuating lift coefficients is conducted to obtain the vortex shedding frequencies f. The normalized vortex shedding frequency Strouhal number is calculated by . Fig. 11 shows the PSD (Power Spectral Density) at different inclined angles. It should be noted that the bandwidth of the PSD for α<45° is broader than that for α≥ 45° due to the harmonic vortex shedding caused by the flipflop gap flow. St 1 and St 2 represent the normalized vortex shedding frequencies of wide wake and narrow wake, respectively. The variation of the Strouhal number with the inclined angle is shown in Fig. 12.
At α=0°, the Strouhal numbers of C1 and C2 are close to the results obtained by Kang (2003). With the increase of the inclined angle α to 37.5°, the values of St 1 and St 2 decrease gradually. As a result, compared with the case of α=0°, St 1 and St 2 decrease to 84.21% and 75.98% at the inclined angle α=37.5°, respectively. The decrease of the Strouhal number indicates that the vortex shedding is suppressed gradually with the increase of the inclined angle. The number of harmonic peaks decreases suddenly and two pronounced sharp peaks are observed at α=45° and 52.5°. However, the Strouhal number decreases dramatically at α=52.5°, St 1 and St 2 at α=52.5° are only 44.79% and 39.52% of their counterparts at α=45°, respectively. At α<45°, there is a small peak in the vicinity of St=0, which is considered to be the occurrence frequency of flip-flop gap flow (Kang, 2003). However, the small peak value is not observed due to the disappearance of flip-flop flow when α>45°.
It can be seen from the fact that both the vortex shedding and flip-flop gap flow are suppressed gradually with the increase of the inclined angle (as shown in Fig. 12). At the large inclined angle α=60°, there is no peak observed in the PSD figure, in good agreement with "no vortex shedding pattern" as previously discussed. The possible reasons may be contributed to two aspects: (1) the complete suppression of the vortex shedding pattern at such a large inclined angle; (2) the relatively low normal Reynolds number Re N =50 (Re N =u n D/v). Fig. 11. Power spectrum density of the fluctuating lift coefficient at different inclined angles.
As shown in Fig. 13, when α≤30°, the values of the drag coefficients are close to the IP (cos 2 α) as well as the lift coefficients, the largest deviation between the drag coefficients and IP is 10.9% at α=30°. When α>30°, the deviations of the drag coefficients and lift coefficients will increase dramatically. However, it appears that the deviation between the lift coefficient and IP is small when α≤52.5°. It can be concluded that the mean drag force coefficients follow the IP when α≤30°. However, when α>30°, the mean drag force coefficients will not be sensitive to the IP.
It can be seen from Fig. 14 that for the normalized vortex shedding frequency St 1 in the narrow wakes, the values of two parameters are similar and close to the IP (cosα) when α≤30°, the deviation between two St 1 and IP is about 10.4% at α=30°. When α>30°, the deviation will increase dramatically. For the normalized vortex shedding frequency St 2 in the wide wakes, the values of two parameters are similar and close to the IP (cosα) when α≤45°, the largest deviation is 7.26% at α=15°. As a result, it appears that the variation of St 1 follows the IP when α≤30°, and that of the St 2 follows the IP when α≤45°.

Conclusions
A series of three-dimensional numerical simulations is performed to study the flow past two side-by-side inclined cylinders at low Reynolds number (Re=100) and small spacing ratio (T/D=1.5). The effects of the inclined angle on the flow field and force coefficients are investigated, and the principle findings of this work are summarized as follows.
(1) The inclined angle has a significant effect on the instantaneous and mean flow fields. The vortex shedding patterns behind two side-by-side cylinders are suppressed with the increase of the inclined angle.
(2) The pattern of the axial flow is significantly affected by the inclined angle. For the cases of α<30°, the recirculation zones with reverse flow behind two side-by-side cylinders are observed. With the increase of the inclined angle to α≥30°, the recirculation zones disappear, and the streamline topologies behave like a pattern without separation points. The maximum mean spanwise velocity approaches   to be 40% of that at α=45°.
(3) For all the cases, there is a repulsive lift force acting on the two cylinders. The mean drag coefficient and r.m.s lift coefficient decrease monotonously with the increase of the inclined angle. Both the flip-flop gap flow and vortex shedding are suppressed gradually with the increase of the inclined angle, leading to the disappearance of harmonic vortex shedding at large inclined angles (α≥45°). The force coefficients and Strouhal numbers are found to basically follow the IP when α≤30° except for the Strouhal numbers in the wide wakes, where St 2 follow the IP when α≤45°.