Numerical Prediction of Hydrodynamic Loading on Circular Cylinder Array in Oscillatory Flow Using Direct-Forcing Immersed Boundary Method

Abstract

Cylindrical structures are commonly used in offshore engineering, for example, a tension-leg platform (TLP). Prediction of hydrodynamic loadings on those cylindrical structures is one of important issues in design of those marine structures. This study aims to provide a numerical model to simulate fluid-structure interaction around the cylindrical structures and to estimate those loadings using the direct-forcing immersed boundary method. Oscillatory flows are considered to simulate the flow caused by progressive waves in shallow water. Virtual forces due to the existence of those cylindrical structures are distributed in the fluid domain in the established immersed boundary model. As a results, influence of the marine structure on the fluid flow is included in the model. Furthermore, hydrodynamic loadings exerted on the marine structure are determined by the integral of virtual forces according to Newton’s third law. A square array of four cylinders is considered as the marine structure in this study. Time histories of inline and lift coefficients are provided in the numerical study. The proposed approach can be useful for scientists and engineers who would like to understand the interaction of the oscillatory flow with the cylinder array or to estimate hydrodynamic loading on the array of cylinders.

1. Introduction

Marine structures are commonly utilized to undertake activities in the offshore region. For example, to obtain oil in the offshore region, a tension-leg platform (TLP) is used. Another example is the wave power station proposed by a Norwegian company [1]. The TLP frame is used for the wave power station. A platform with four circular cylinders is used as the main frame of the proposed wave power station. Since those marine structures are in the offshore region, hydrodynamic loading from waves and currents should be considered while one designs this sort of marine structures. Circular cylinders are often utilized as the frame of the marine structure, so prediction of the loadings on the circular cylinders is important. The aim of this study is to establish a numerical model to predict hydrodynamic loading exerted on those circular cylinders in marine structures.

Oscillatory flows are regarded as the flows due to a progressive wave train in shallow water, so interaction of oscillatory flows with cylinders have been studied by several researchers. Sumer and Fredsoe [2] reviewed published papers regarding studies of flow around a single cylinder in oscillatory flows. Keulegan-Carpenter (KC) number was used to categorize the flow variations into a variety of flow regimes. The study of hydrodynamic loadings exerted on a single circular cylinder in oscillatory flows can be found in numerous publications. Williamson [3] conducted a series of experiments to measure inline force and lift of an oscillating cylinder. Flows were visualized at various KC numbers. Sarpkaya [4] provided theoretical and experimental results of inertia coefficient and flow visualization. Obasaju et al. [5] measured forces on a circular cylinder in oscillating flow. They categorized the flow variations into asymmetric, the transverse, the diagonal, the third vortex, and the quasi-steady regimes as KC varied from 4 to 55.

Numerical studies on the force prediction for a single cylinder in oscillatory flows have been reported in the past two decades. Lin et al. [6] employed a hybrid Lagrangian/Eulerian discrete vortex method to simulate flow around a circular cylinder in oscillatory flow up to KC numbers of 30. Iliadis and Anagnostopoulos [7] considered viscous oscillatory flows interacting with a circular cylinder by solving streamfunction and vorticity formulation of the Navier-Stokes equation. The finite element method was adopted in their study. They visualized flow variations using vorticity contours at various KC numbers. The inline and the transverse forces on the circular cylinder were determined and reported. Zheng and Dalton [8] utilized the finite difference method to study the interaction between a rectangular cylinder and oscillatory flows. Variations of flow due to the change of KC number were visualized in their study. Prediction of forces on the square cylinder was also predicted in their model. An et al. [9] investigated the oscillatory flow around a circular cylinder at low KC numbers. Honji vortices were simulated and explained in their study.

Interaction of oscillatory flows with a number of cylinders has been reported in several papers. Williamson [3] studied flow around a pair of circular cylinders in his flow visualization experiments. Those two cylinder were either side-by-side or at 45°. The effect of gap between two cylinders on flow variation and forces were reported. Chern et al. [10] simulated oscillatory flows around a pair of side-by-side cylinders up to KC of 15. Phase diagrams of force components were revealed to observe the route to chaos in the flow. The relationship between two vortical systems behind cylinders were analyzed. Anagnostopoulos and Dikarou [11] established a numerical model to simulate viscous oscillatory flow past four cylinders. The cylinders were in a square arrangement. Inline and transverse forces of the cylinders were predicted and the relationship among them and KC was given in the study. Also, the effect of pitch ratio on the viscous oscillatory flows was discussed.

An immersed boundary method which added a virtual force in the momentum equations to simulate the existence of solid has been receiving more and more attentions. One of immersed boundary methods is the so-called direct forcing method proposed by Yusof [12]. The direct forcing method determine a forcing term by calculating the difference between the interpolated velocities on the boundary points and the desired solid boundary velocities. As a result, it does not require to fit grids for a complex configuration. A Cartesian grid can be used in the direct-forcing immersed boundary method. The principle of the direct forcing method has been adopted and applied to several fluid-structure interactions in Fadlun et al. [13], Tseng and Ferziger [14], and Verzicco et al. [15]. The present model based on the direct forcing method was also validated by the uniform flow past cylinders in the work by Noor et al. [16]. In terms of those studies, it turns out that the direct-forcing immersed boundary method is able to simulate fluid-solid interactions and to predict hydrodynamic loadings on cylinders properly.

The aim of this study is to utilize the direct-forcing immersed boundary method to study the oscillatory flow around a cylinder array in a square arrangement. Since the geometric configuration of the cylinder array is complex in the flow domain, one can take the advantage of the propose direct-forcing immersed boundary method to simulate the oscillating flow with the complex cylinder configuration in Cartesian grids. The other advantage of the direct-forcing immersed boundary method is the prediction of forces on a solid immersed in fluids. The virtual force in the model can be used to determine forces. It is common that forces on a single cylinder in an oscillatory flow are divided to two parts, the drag and inertia forces (see Williamson [3] and Sarpkaya [4]). Nevertheless, the resultant forces of the cylinder in this study can be calculated by the integral of the virtual force in the solid region. Consequently, the oscillatory flow interacting with the cylinder array and the hydrodynamic loadings on the cylinder array can be obtained using the proposed model.

2. Mathematical Formulae and Numerical Model

In the present work, the direct forcing method is used to establish the proposed immersed boundary method. A virtual force is added to the incompressible Navier-Stokes equations in order to accommodate the interaction between solids and fluids. A solid body is identified by a volume-of-body function 𝜂 which denotes a fraction of solid within a cell. For a cell full of solids, 𝜂 is equal to 1 while it becomes 0 for a cell full of fluids. It is fractional in a boundary cell which consists of solids and fluids both.

2.1. Equations of Fluid Motion

An incompressible fluid is considered in the study. The motion of fluids should satisfy the conservation laws of mass and momentum. Mathematically, those physical laws can be denoted in dimensionless form and shown as𝜕𝐮∇⋅𝐮=0,(2.1)𝜕𝑡+∇⋅(𝐮𝐮)=−∇𝑝+R1e∇²𝐮+𝜂𝐟,(2.2) where 𝐮 and 𝑝 are dimensionless velocity and pressure, respectively. Equations (2.1) and (2.2) are nondimensionalized by the characteristic velocity 𝑈_𝑚 which is the amplitude of the oscillating velocity and the characteristic length 𝐷 which is the diameter of a circular cylinder. Meanwhile, time is nondimensionalized by the characteristic time 𝐷/𝑈_𝑚. Re is Reynolds number and defined by 𝑈_𝑚𝐷/𝜈, where 𝜈 is kinematic viscosity. The proposed immersed boundary model adds the virtual force 𝐟 to include the effect of solid in the viscous flow.

The momentum equation (2.2) is solved by three steps. First, the velocity is stepped from the 𝑚th time level to the first intermediate level 𝐮^∗ by solving (2.2) without the pressure gradient and virtual force at the beginning of each time step. This step is implemented by the following formula: 𝐮^∗−𝐮^𝑚Δ𝑡=𝑆^𝑚,(2.3) where 𝑆 includes the convective and diffusive terms of (2.2). Apparently, the first intermediate velocity 𝐮^∗ in (2.3) does not satisfy (2.1), so 𝐮^∗ is marched to the second intermediate velocity 𝐮^∗∗ by including the pressure gradient, 𝐮^∗∗−𝐮^∗Δ𝑡=−∇𝑝^𝑚+1/2,(2.4) at the second step. Taking the divergence of (2.4) gives ∇⋅𝐮^∗∗−∇⋅𝐮^∗Δ𝑡=−∇²𝑝^𝑚+1/2.(2.5) In principle, the second intermediate velocity 𝐮^∗∗ should satisfy the continuity equation, that is, ∇⋅𝐮^∗∗=0.(2.6) Substitution of (2.6) to (2.5) results in the Poisson equation of pressure: ∇²𝑝^𝑚+1/2=1Δ𝑡∇⋅𝐮^∗.(2.7) Once (2.7) is solved, 𝐮^∗∗ will be obtained by (2.4). So far, the effect of solid on the viscous fluid flow has not been included. The virtual force 𝐟 caused by a solid should be involved at the third step, so 𝐮^∗∗ is updated to the (𝑚+1)th time level by imposing the virtual force term, that is, 𝐮^𝑚+1−𝐮^∗∗Δ𝑡=𝜂𝐟^𝑚+1.(2.8) The dimensionless virtual force 𝐟 reveals the existence of a force to hold or to drive a solid body which is either stationary or moving. Given that the velocity of a solid is 𝐮_𝑠 and may be not the same as the calculated velocity 𝐮^∗∗,  𝐟 acting on the solid will have to ensure that the fluid velocity 𝐮^𝑚+1 is equal to the solid velocity 𝐮_𝑠 at the (𝑚+1)th time step, that is, 𝐮^𝑚+1=𝐮^𝑠𝑚+1. Hence, the virtual force is defined as the rate of momentum change of a solid body and proportional to the difference between the solid velocity at the (𝑚+1)th time step and the fluid velocity at the 𝑚th time step. The force exists at the fluid domain, where the solid body is immersed and zero elsewhere. Furthermore, it can be simply written as 𝐟^𝑚+1=𝐮^𝑚+1−𝐮^∗∗Δ𝑡=𝐮^𝑠𝑚+1−𝐮^∗∗Δ𝑡.(2.9) Presumably, cylinders are fixed so 𝐮_𝑠 is always zero for all cylinders.

2.2. Oscillatory Flow Boundary Condition

Oscillatory flows are considered in this study. Transient velocity boundary conditions are imposed at four boundaries of the computational domain to simulate oscillatory flows. Consider an oscillatory flow of dimensionless period 𝑇. The dimensionless horizontal velocity component of the oscillatory flow varies according to the condition 𝑢=sin2𝜋𝑡𝑇.(2.10) This boundary condition has been used by Zheng and Dalton [8] and Chern et al. [10] for simulations of oscillatory flows with square cylinders. The dimensionless vertical velocity component of the oscillatory flow at four boundaries of the computational domain is set as zero during simulations.

2.3. Calculation of Hydrodynamic Force on Cylinder

The integration of the virtual force will be a good approximation of the resultant dimensionless force exerted on a single circular cylinder, 𝐅=−_Ω𝜂𝐟𝑑𝐴,(2.11) where 𝐅 is the resultant hydrodynamic force vector and 𝐴 is the area occupied by the circular cylinder. The inline force coefficient 𝐶_𝑓 in the oscillatory flow direction can be denoted as 𝐶_𝑓=−2𝐹_𝑥(2.12) and the life or transverse force coefficient 𝐶_𝑙 can be determined by 𝐶_𝑙=−2𝐹_𝑦.(2.13)

2.4. Numerical Procedure to Predict Fluid-Solid Interaction

The numerical procedure to predict viscous oscillatory flows and to trace the falling object at each time step can be summarized in the following steps.(1)Determine 𝜂 through the positions and diameters of cylinders at the first time step.(2)Calculate 𝐮^∗ via (2.3).(3)Solve the pressure Poisson equation (2.7) and advance to 𝐮^∗∗ via (2.4).(4)Compute the virtual force required via (2.9).(5)Determine the velocity via (2.8) and back to Step 2 and for the next time step.

Uniform Cartesian grids are adopted in this study. 250 × 250 and 430 × 430 grids are used for cases for a single cylinder and four cylinders, respectively. The finite volume method is used to solve the momentum equations. The advective scheme is discretized by the so-called third QUICK scheme. The Adams-Bashforth scheme is used for the temporal derivative. The dimensionless time increment Δ𝑡 is set as 10⁻⁴ which satisfies the CFL condition so the stability can be retained in the model. The total dimensionless time for simulations is 230. It takes about 2 and half days for each simulation of the 2-D oscillatory flow around a single cylinder at a PC cluster consisting of Intel Xeon CPUs.

3. Results and Discussion

3.1. Validation of Numerical Model

Although the proposed numerical immersed boundary model was validated in our previous study (Noor et al. [16]) which considered a uniform past circular cylinders, it will be validated by the study of the interaction of oscillatory flow with a single circular cylinder again. Figure 1 shows the schematic of the benchmark test problem. Consider a circular cylinder of diameter 𝐷 in a fluid domain. The oscillatory flow is featured by Keulegan-Carpenter (KC) number which is determined by the formula 𝑈_𝑚𝑇/𝐷. Cases at KC = 2 and 10 and Re = 200 are simulated in the validation test. Figure 2 presents the evolution of vorticity contours during a cycle at KC = 2. It is found that a pair of symmetric vortices appearing in two sides of the cylinder alternatively due to the change of flow direction. The results agree with Iliadis and Anagnostopoulos’ study [7]. Furthermore, the development of the symmetric wake in half a cycle is calculated and shown in Figure 3. The wake is elongated as time increases. The predicted wake length in the model is very close to Iliadis and Anagnostopoulos’ result. As KC increases, the symmetry in the pair of vortices breaks, Figure 4 shows the evolution of shedding vortices in a cycle. It is found that vortices are not attached to the cylinder anymore. They shed from the cylinder and travel downstream. Finally, they are damped in the far field region. In addition, the inline force coefficient 𝐶_𝑓 is estimated using the proposed approach in the previous subsection. The time histories of 𝐶_𝑓 at KC = 2 and 10 are shown in Figures 5(a) and 5(b), respectively. It is found that 𝐶_𝑓 reaches a maximum while the oscillatory flow changes its direction. The comparison of the predicted 𝐶_𝑓 reveals that the present immersed boundary model can simulate the interaction of oscillatory flow with a single cylinder and the hydrodynamic loading reasonably.

Figure 1 

Schematic of interaction of an oscillatory flow with a single circular cylinder.

Figure 2 

Evolution of vorticity distribution at KC = 2 and Re = 200. Dash lines refer to the negative vorticity contours.

Figure 3 

Evolution of wake during half a cycle at KC = 2 and Re = 200.

Figure 4 

Evolution of vorticity distribution at KC = 10 and Re = 200. Dash lines refer to the negative vorticity contours.

(a) KC = 2

(b) KC = 10

Figure 5 

Time histories of 𝐶𝑓 at (a) KC = 2 and (b) KC = 10.

3.2. Interaction of Oscillatory Flow with Cylinder Array

A cylinder array is commonly used in marine structures. For example, a tension-leg platform which is normally used for the offshore production of oil or gas consists of a group of vertical cylinders as mentioned in Introduction. Those cylinders have to endure hydrodynamic loadings from waves and currents. To design those cylindrical structures, stress distribution caused by those hydrodynamic loadings must be provided in advance. Hence, the hydrodynamic loadings due to oscillatory flows are calculated using the current immersed boundary model. An array of four circular cylinders are considered in this study as shown in Figure 6. Four cylinders are allocated in a square arrangement. The distance between two centers of two adjacent cylinders is denoted as 𝑃. The pitch ratio P/D is one of parameters to affect the hydrodynamic behavior around those cylinders. To simplify the problem, P/D is set as 2 in this study. Also, Re is fixed at 200 in this subsection. To investigate the effect of KC number, three various KC numbers 2, 5, and 10 are considered to observe its effect on the oscillatory flow. Results are illustrated in the following subsections.

Figure 6 

Schematic of interaction of an oscillatory flow with a cylinder array. 𝑃/𝐷 = 2.

3.2.1. Evolution of Four Vortical Systems

Figure 7 shows the evolution of vorticity contours around four cylinders in a cycle. Each cylinder has its symmetric vortical system as mentioned in the case with a single cylinder. Those vortical systems do not interact with others since 𝑃/𝐷 of 2 is wide. Nevertheless, the symmetry in those vortical systems vanishes as KC is increased to 5. Figure 8 presents the snapshots of vorticity contours for the 11th cycle at KC = 5. It is found that the gap flow influence the vortical systems, so vortices are pushed downstream. The vortices disappear after they leave from the cylinders. Figure 9 presents vorticity contours at KC = 10. It is found that vortices become irregular and shed from cylinders. They are not damped until they travel far away from cylinders.

(a)

(b)

(c)

(d)

Figure 7 

Snapshots of vorticity contours in the oscillatory flow interacting with four cylinders during a cycle at KC = 2 and Re = 200.

(a)

(b)

(c)

(d)

Figure 8 

Snapshots of vorticity contours in the oscillatory flow interacting with four cylinders during a cycle at KC = 5 and Re = 200.

(a)

(b)

(c)

(d)

Figure 9 

Snapshots of vorticity contours in the oscillatory flow interacting with four cylinders during a cycle at KC = 10 and Re = 200.

3.2.2. Variation of 𝐶_𝑓 with KC

Prediction of hydrodynamic forces are important for the design of the cylinder array. To understand the variation of those forces with respect to KC, the time histories of 𝐶_𝑓 of the first cylinder as shown in Figure 6 are provided. Figure 10 demonstrates the time histories of 𝐶_𝑓. 𝐶_𝑓 behaves sinusoidaly at KC = 2 as shown in Figure 10(a). It looks just like the case with a single cylinder. As KC increases to 5, the sinusoidal form of 𝐶_𝑓 is not regular any more as shown in Figure 10(b). Also, the amplitude of 𝐶_𝑓 decreases from 10 to 4 at KC = 5. As KC increases to 10, 𝐶_𝑓 becomes more irregular and also the amplitude decreases again.

(a)

(b)

(c)

Figure 10 

Time histories of 𝐶𝑓 of the first cylinder at KC = 2, 5, and 10.

3.2.3. Variation of 𝐶_𝑙 with KC

The other force component is along the transverse direction. It is called the life or transverse force. When KC = 2, as Figure 7 shows four vortical systems are symmetric during a cycle. As a result, lifts or transverse forces on four cylinders are very little in comparison with 𝐶_𝑓. Nevertheless, as KC increases to 5, the vortical systems are not symmetric. The transverse forces on cylinders are obvious. Figure 11 shows time histories of 𝐶_𝑙 at KC = 5 and 10. 𝐶_𝑙 at KC = 5 is not regular is smaller than 𝐶_𝑓. While KC increases to 10, the frequency of variation of 𝐶_𝑙 becomes faster. Another feature at KC = 10 is those spikes. They do not appear regularly but randomly. This phenomenon does not occur in 𝐶_𝑓.

(a)

(b)

Figure 11 

Time histories of 𝐶𝑙 of the first cylinder at KC = 5 and 10.

3.2.4. Phase Diagram of 𝐶_𝑓 versus 𝐶_𝑙

The resultant hydrodynamic forces can be regarded as the responses of the vortical systems around those four cylinders. Trajectories of 𝐶_𝑓 versus 𝐶_𝑙 at successive instants reveal the states of the vortical systems. Figure 12 shows phase diagrams of 𝐶_𝑓 and 𝐶_𝑙 of the first cylinder at various KC numbers. While KC = 2, Figure 12(a) shows that the vortical system around the first cylinder is almost periodic. Therefore, the trajectory at KC = 2 is close to a straight line and the vortical system follows the line up and down. As KC increases to 5 and the transverse force becomes significant, the trajectory as shown in Figure 12(b) looks like a twisted “8”. Nevertheless, the trajectory shown in Figure 12(c) becomes completely chaotic as KC increases to 10. In terms of those sub figures, it illustrates the way of the vortical system changes from a periodic state to a chaotic state due to the increase of KC.

(a)

(b)

(c)

Figure 12 

Phase diagrams of 𝐶𝑓 versus 𝐶𝑙 of the first cylinder at KC = 2, 5, and 10.

4. Conclusions

A direct forcing immersed boundary model has been established to predict interactions between oscillatory flows with a cylinder array. The proposed immersed boundary model was validated by an oscillating flow interacting with a single cylinder at Keulegan-Carpenter (KC) number = 2 and 10 and Reynolds number = 200. The development of the symmetric wake predicted by the proposed model agreed with other researchers’ study at KC = 2. Meanwhile, time histories of the inline force coefficient 𝐶_𝑓 at KC = 2 and 10 were compared with the same published study. Good agreements were found between the present results and the other study. The established numerical model was further applied to simulate oscillatory flows around four cylinders at moderate KC numbers. As KC was 2, the symmetry was found in all vortical systems around four cylinders. The resultant 𝐶_𝑓 of one of cylinders was sinusoidal and the lift coefficient 𝐶_𝑙 was very small at KC = 2. As KC increased to 5, the symmetry in vortical systems broke. 𝐶_𝑓 became irregular and its amplitude decreased and 𝐶_𝑙 was significant. While KC was 10, the vortical systems were more complex. The change of 𝐶_𝑓 became more frequent and a number of spikes were found in 𝐶_𝑙. Phase diagrams of 𝐶_𝑓 versus 𝐶_𝑙 were used to demonstrate the state vortical system around the cylinder. It showed that the trajectory changed from a straight line to a twisted “8” as KC changed from 2 to 5. Finally, when KC was 10, the trajectory in the phase diagram was completely chaotic. It meant that the vortical system around the cylinder was completely disordered.

Nomenclatures

Latin Symbols

𝐶𝑓:	Inline force coefficient
𝐶𝑙:	Lift coefficient
𝐷:	Diameter of cylinder
𝐟:	Dimensionless virtual force
𝐅:	Total dimensionless force acting on solids
KC:	Keulegan-Carpenter number, 𝑈𝑚𝑇/𝐷
𝑃:	Dimensionless distance between centers of two adjacent cylinders
𝐿:	Dimensionless length of wake
𝑝:	Dimensionless pressure
𝑅:	Dimensionless radius of cylinder
Re:	Reynolds number, 𝑈𝑚𝐷/𝜈
𝑇:	Dimensionless period of oscillating velocity
𝑡:	Dimensionless time
Δ𝑡:	Dimensionless time increment
𝑈𝑚:	Amplitude of oscillating velocity, m⋅s−1
𝐮(𝑢,𝑣):	Dimensionless velocity
𝑥,𝑦:	Dimensionless horizontal and vertical cartesian coordinates.

Greek Symbols

𝜂:	The volume of solid function
𝜈:	Kinematic viscosity of fluid, m2⋅s−1.

Subscripts

𝑠:

Solid.

Superscripts

𝑚:	Time step level
∗:	Intermediate time step level.

Acknowledgment

The authors would like to express their gratitude for the financial support from National Science Council, Taiwan (Grant no.: NSC 100-2212-E-011-163).