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

Cylindrical structures are commonly used in offshore engineering, for example, a tensionleg 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.


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 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.

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.

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 In principle, the second intermediate velocity u * * should satisfy the continuity equation, that is, Substitution of 2.6 to 2.5 results in the Poisson equation of pressure: Once 2.7 is solved, u * * will be obtained by 2.4 .So far, the effect of solid on the viscous fluid flow has not been included.The virtual force f caused by a solid should be involved at the third step, so u * * is updated to the m 1 th time level by imposing the virtual force term, that is, The dimensionless virtual force f 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 u s and may be not the same as the calculated velocity u * * , f acting on the solid will have to ensure that the fluid velocity u m 1 is equal to the solid velocity u s at the m 1 th time step, that is, . 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 m 1 th time step and the fluid velocity at the mth 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 Presumably, cylinders are fixed so u s is always zero for all cylinders.

Oscillatory Flow Boundary Condition
Oscillatory

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, where F is the resultant hydrodynamic force vector and A is the area occupied by the circular cylinder.The inline force coefficient C f in the oscillatory flow direction can be denoted as 12 and the life or transverse force coefficient C l can be determined by

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.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 Δt is set as 10 −4 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.

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 D in a fluid domain.The oscillatory flow is featured by Keulegan-Carpenter KC number which is determined by the formula U m T/D.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 C f is estimated using the proposed approach in the previous subsection.The time histories of C f at KC 2 and 10 are shown in Figures 5 a and 5 b , respectively.It is found that C f reaches a maximum while the oscillatory flow changes its direction.The comparison of the predicted C f reveals that the present immersed boundary model can simulate the interaction of oscillatory flow with a single cylinder and the hydrodynamic loading reasonably.

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 P .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.

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 P/D 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.

Variation of C f 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 C f of the first cylinder as shown in Figure 6 are provided.Figure 10 demonstrates the time histories of C f .C f 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 C f is not regular any more as shown in Figure 10 b .Also, the amplitude of C f decreases from 10 to 4 at KC 5.As KC increases to 10, C f becomes more irregular and also the amplitude decreases again.

Variation of C l 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 C f .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 C l at KC 5 and 10.C l at KC 5 is not regular is smaller than C f .While KC increases to 10, the frequency of variation of C l becomes faster.Another feature at KC 10 is those spikes.They do not appear regularly but randomly.This phenomenon does not occur in C f .

Phase Diagram of C f versus C l
The resultant hydrodynamic forces can be regarded as the responses of the vortical systems around those four cylinders.Trajectories of C f versus C l at successive instants reveal the states of the vortical systems.Figure 12 shows phase diagrams of C f and C l 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.

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 C f 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 C f of one of cylinders was sinusoidal and the lift coefficient C l was very small at KC 2. As KC increased to 5, the symmetry in vortical systems broke.C f became irregular and its amplitude decreased and C l was significant.While KC was 10, the vortical systems were more complex.The change of C f became more frequent and a number of spikes were found in C l .Phase diagrams of C f versus C l 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 Diameter of cylinder f: Dimensionless virtual force F: Total dimensionless force acting on solids KC: Keulegan-Carpenter number, U m T/D P : Dimensionless distance between centers of two adjacent cylinders L: Dimensionless length of wake p: Dimensionless pressure R: Dimensionless radius of cylinder Re: Reynolds number, U m D/ν T : Dimensionless period of oscillating velocity t: Dimensionless time Δt: Dimensionless time increment U m : Amplitude of oscillating velocity, m • s −1 u u, v : Dimensionless velocity x, y: Dimensionless horizontal and vertical cartesian coordinates.
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 T .The dimensionless horizontal velocity component of the oscillatory flow varies according to the condition u sin 2πt T .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.

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

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

Figure 3 :
Figure 3: Evolution of wake during half a cycle at KC 2 and Re 200.

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

Figure 5 :Figure 6 :
Figure 5: Time histories of C f at a KC 2 and b KC 10.

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

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

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

Figure 10 :
Figure 10: Time histories of C f of the first cylinder at KC 2, 5, and 10.

Figure 11 :
Figure 11: Time histories of C l of the first cylinder at KC 5 and 10.

Figure 12 :
Figure 12: Phase diagrams of C f versus C l of the first cylinder at KC 2, 5, and 10.
ηf, 2.2 where u and p are dimensionless velocity and pressure, respectively.Equations 2.1 and 2.2 are nondimensionalized by the characteristic velocity U m which is the amplitude of the oscillating velocity and the characteristic length D which is the diameter of a circular cylinder.Meanwhile, time is nondimensionalized by the characteristic time D/U m .Re is Reynolds number and defined by U m D/ν, where ν is kinematic viscosity.The proposed immersed boundary model adds the virtual force f 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 mth time level to the first intermediate level u * 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: