3D effects on vortex-shedding flow and hydrodynamic coefficients of a vertically oscillating cylinder with a bottom-attached disk

Vortex-shedding flow induced by the vertical oscillation of a cylinder with bottom-attached disks of different diameter ratio Dd/Dc and thickness ratio td/Dc is studied by a 3D (three-dimensional) numerical model developed in this paper, and compared with the results obtained through 2D (two-dimensional) numerical model. The high-order upwind scheme is applied to stabilize the computation, and convergence is accelerated by the multi-grid method. Qualitative and quantitative analyses of the differences between the 2D and 3D simulation results reveal the 3D effect on the flow field characteristics and hydrodynamic coefficients of the vertically oscillating cylinder with a bottom-attached disk. The 3D effect on the fluid field is mainly reflected in the significance of three vortex-shedding patterns: ωx has a greater effect on the flow fields around the sharp edges relative to the vortices generated in the 2D simulation. In the slice along the axial orientation, the vortex effect of ωy along the radial axis is smaller than that of ωx along the circumferential direction, indicating the radial effect on the velocity more pronounced than the circumferential effect around the sharp edges of the disk. The rotational interaction ωz of the fluid in the horizontal plane during the heave motion is insignificant. Based on the 2D and 3D simulation results, the turning points that separate the increasing regimes of the added mass coefficient and damping ratio are identified. The dependence of the turning point on the diameter ratio Dd/Dc and thickness ratio td/Dc are discussed in detail.


Introduction
To enhance the damping mechanism during the resonant oscillation of a TLP (Tension Leg Platform) or Spar in the vertical direction, additional mechanical damping devices or other active damping systems are introduced externally. A typical example of these damping devices is a heave plate, which is attached to the bottom of a TLP cylinder or Spar (Fig. 1). The additional disk at the bottom of the cylinder enhances the vortex-shedding process, which consequently changes the hydrodynamic properties by introducing extra damping and increasing the added mass of the system.
The problem of interest can be modelled as a cylinder with a disk attached at the bottom. The system vertically oscillates in an infinite fluid as demonstrated in Fig. 1. The flow around this system is governed by two fundamental β non-dimensional parameters, the Keulegan-Carpenter number (KC) and the frequency parameter ( ), which are defined as follows: (1) β where D c is the diameter of the cylinder. a and f are the amplitude and frequency of the oscillation, respectively.
And v is the kinematic viscosity of the fluid. The parameters KC and represent the dimensionless amplitude and frequency of the oscillation, respectively, as are related to Reynolds number as follows: China Ocean Eng., 2017, Vol. 31, No. 4, P. 428-437 DOI: 10.1007/s13344-017-0049-7, ISSN 0890-5487 http://www.chinaoceanengin.cn/ E-mail: coe@nhri.cn is the velocity amplitude of the oscillation, and is the circular frequency of the oscillation. β The hydrodynamic performance of the vertical oscillation of a cylinder with a disk attached at the bottom has been experimentally investigated by Thiagarajan and Troesch (1998) with the particle image velocimetry (PIV) technique. They find that the vortex-shedding pattern is depended on both flow parameters (KC and ) and the geometry of the structure. Thiagarajan (2003a, 2003b) found that three distinct vortex-shedding modes/regimes (independent, interactive and unidirectional vortex shedding) occur due to different KC values and the thickness of the disk (t d ). Diffusion or convection is found to be predominant in each vortex-shedding regime. Consequently, distinct hydrodynamic heave damping behaviours are found to be associated with different vortex-shedding regimes, i.e., a lower, weakly KCdependent damping was found in the independent vortexshedding regime, whereas the damping was found to be non-linearly dependent on KC in the latter two vortex-shedding regimes as a result of the strong convection effect. Li et al. (2013) conduct systematic experimental studies on the hydrodynamic characteristics of heave plates under forced oscillation. The effects of variables such as the Keulegan-Carpenter (KC) number, frequency of oscillation, plate depth, thickness-to-width ratio, shape of the edge, perforation ratio and hole size on the hydrodynamic coefficients are analysed. Using 3D direct numerical simulations, Yang et al. (2014) investigate the hydrodynamic characteristics of a heave plate under steady current conditions. The changes of flow field and the parameter dependence of the hydrodynamic coefficients on the flow field are revealed in the study.
This study aims to investigate the 3D effect on the flow characteristics and associated hydrodynamic coefficients for an axial oscillating cylinder with a disk attached at the bottom. Numerical simulations are performed for systematic variation of the disk geometry. The diameter ratio and thickness ratio t d /D c are listed in Table 1. The 3D effect are revealed based on quantitative analysis of the differences between the 2D and 3D simulation results. A finite difference solution of the viscous fluid structure interaction problem is obtained by solving Navier-Stokes equations for an incompressible viscous flow. The present numerical method is validated against several benchmark problems and available experimental results.

Governing equations and boundary conditions
The fluid in Fig. 1 can be modelled by the unsteady incompressible Navier-Stokes equations and continuity equation: ∇ · u = 0; (4) where u is the velocity vector of the flow field. p is the dynamic pressure of the flow field, and is the numerical Reynolds number deduced by the nondimensionalization in the numerical simulation. The cylinder with an attached disk is forced to oscillate sinusoidally along its vertical axis as: ω where a and are the amplitude and circular frequency of the oscillation, respectively.
No-flux or no-slip velocity boundary condition is imposed on the oscillating cylinder and disk surfaces: where S 0 is the wet surface of the oscillating cylinder and the disk. As Thiagarajan and Troesch (1998) observed that the disturbance of the free surface is very small due to the small-amplitude and high-frequency heave oscillation of the floating structures. Therefore, nonlinear free surface effects are not included in the present numerical modelling. A freestream flow boundary condition is assumed, and the pressure on the free surface is constant. Free-stream velocity conditions are also applied to the entire open boundary, and static pressure is assumed at the far open boundaries. The Neumann boundary condition for the pressure on the cylinder and disk surfaces is obtained by applying the momentum equations in the direction normal to the cylinder and disk surfaces.  WANG Quan et al. China Ocean Eng., 2017, Vol. 31, No. 4, P. 428-437 429 2.2 Hydrodynamic force calculation The hydrodynamic force that acts on the structure that oscillates along the vertical direction is calculated by integrating the vertical component of the stress vector, which is the vertical projection of the stress tensor along the wet surface of the cylinder and disk, where S 0 is the wet surface of the cylinder and the disk. p is the dynamic pressure, and are the components of the unit-normal vectors along the structure wet surface. The right-hand-side terms of the above expression represent three different components of the viscous force, the dynamic pressure component the viscous shear stress component and the viscous normal stress component The hydrodynamic damping and added mass characteristics of the structure that oscillates vertically are calculated from the time histories of F(t), which were obtained through numerical simulation using a Fourier average method. The viscous drag force is represented by a similar expression to the Morison equation, in which the drag force is proportional to the square of the velocity, ∀ where S is the cross section area of the disk. is the immersed volume of the cylinder and disk; C d and C a are the drag and added mass coefficient, respectively. The added mass and drag coefficients are calculated using the Fourier average method (Sarpkaya, 2010), (13) is the amplitude of the oscillation velocity, and T is the period of oscillation. Alternatively, the drag force can be expressed with the equivalent linear damping coefficient B, (15) The equivalence between the Morison-like formula of the drag force and its equivalent linear counterpart results in the following relation between the drag coefficient C d and equivalent linear damping coefficient B, μ where is the dynamic viscosity. D d is the diameter of the disk, and D c is the diameter of the cylinder. The linear damping coefficient B is further normalized to give the nondimensional damping ratio, where is the critical damping of the system. m and m a represent the mass and added mass, respectively. By substituting the mass m+m a with and B with Eq. (16), the damping ratio is, where T d and t d represent the draft of the cylinder and the thickness of the disk, respectively.

Numerical scheme
A finite difference method based on the projection method (Brown, 2001) is used to solve the governing equations. The pressure field is obtained by solving the following Poisson equation, δt where n is the time step. is the time increment. D n+1 is assumed to be zero. But D n is retained as a correction term.
For the temporal integration of the Navier-Stokes equations, the Crank-Nicolson implicit scheme is used. To linearize the convection term, the former step value is used for the convection velocity, (20) Eqs. (19) and (20) are iteratively solved at each time step by successive over relaxation (SOR) coupled with the multi-grid method, which is effective to accelerate the convergence speed.
(21) To handle the complex geometry, the boundary-fitted curvilinear grids are adopted in the numerical simulation. The relations in the former sections, which are described in the Cartesian coordinates, are transformed to the curvilinear coordinate form based on the chain rules of derivations (Thompson et al., 1985).

Model validation β
To obtain the numerical results for the present investigations, the mesh configurations for different computational domains were tested with KC=1.5 and =89236, which produces the most severe flow field change in the studied cases. Because the vortex formation around the structure depends on the local flow field around the sharp edges of the structure (Tao et al., 2000), the vortex-shedding and damping coefficients are affected by the restricted domain around the structure. A computational domain of (11R×15.4R) was ad-opted based on the domain size with the far-field effect minimized. The convergence trends of different spatial and time discretization combinations shown in Table 2 were tested.
The numerical results of the tested cases are sensitive to the grid numbers/spatial discretization scheme. The time step is not the main factor as long as the spatial grid number is within the convergence range. More details of the testing results are listed in Table 2. Based on these convergence tests, the minimum mesh size Δx min =0.05 with the grid number n x ×n y ×n z =128×128×32 and time step Δt=T/1000, where T is the dimensionless oscillation period, was adopted for further investigation.
Before investigation cases, the code was validated against the experimental data reported by Tao and Dray (2008). The validation results are shown in Fig. 2. The variation trend of the numerical results is consistent with experimental test results with small constant variation. To maintain consistency with the result of Tao and Dray (2008), the added mass coefficients in the model are further rewritten in a nondimensional form A as follows: is the theoretical ideal fluid added mass, and m is the displacement mass of the solid part. It m = (1/4) ρπD 2 t should be noted that in the experiment because the disk is the focus of the experiment. In the numerical simulation, the displacement mass of the vertically oscillating cylinder with a disk is the primary concern. The constant variation between the numerical and experimental results may come from the simplification of the numerical simulation. The diameter of the rod, which drives the oscillation of the structure, is unclear from the available data. Meanwhile, the effect of the rod/cylinder is not included in the deduction of the experimental results. Those factors result in inequality between the numerical results and the experimental results. If those secondary effects are ignored, the validation results are acceptable for this study.  Fig. 3-Fig. 6. The 2D and 3D results show the same vortex generation trend, where the forced oscillation of the structure generates, sheds, develops and dissip-   ates vortices around the two sharp edges of the disk. In the 2D simulation, only one vortex-shedding pattern is observed, whereas in the 3D simulation, and are also obtained in addition to , which reveal more information in the flow field of the forced oscillation of the structure. Here, , and are three vortex components of the vortex vector in the Cartesian coordinates, which is the coordinate system adopted in the numerical simulation to unify and simplify the programming procedure. In fact, the vortex vector can be expressed in both Cartesian and cylindrical coordinates: θ ω r ω θ ω ′ z where i, j and k are the base frame of the Cartesian coordinate system. r, and z are the base frame of the cylindrical coordinate system, and , and are the vortex components in the cylindrical coordinates. It should be noted that in the local system of the slice in Fig. 4-Fig. 6, the following equivalent relations exist, then, To unify the coordinate system, all results presented in this paper are based on the Cartesian coordinates, although the discussions about the vortex-shedding pattern in this section are based on the cylindrical coordinates. The 3D effects on the fluid field characteristics are revealed from the vortex-shedding patterns as follows. As demonstrated in Fig. 4, among the vortex-shedding patterns of the 3D results, show nearly the same vortexshedding pattern as the 2D results with the opposite sign, which is the result of the coordinate transformation. The same trend between the 2D results and may originate from the identical direction perpendicular to the presented plane. However, the vortex generated in the 2D simulation tends to concentrate in the local area of the disk edges, whereas affects a larger area around the sharp edges. This may indicate a larger pressure change along the sur-face of the longitudinally vertical oscillating structure. indicates that the effect of the vortex along the radial axis is smaller than that of its counterpart . In other words, the radial effect on the fluid velocity is stronger than the circumferential one around the sharp edges of the disk. The affected area of the vertical vortex-shedding pattern is smaller than that by . The value of is closer to zero with the effect negligible, as shown in Fig. 6. The negligible value of indicates that the rotational interaction of the fluid in the horizontal plane during the heave motion is insignificant. This confirms the rationality of the 2D asymmetric assumption from another perspective. The formation of the vertical vortex-shedding pattern , i.e., the rotational interaction of fluid in the horizontal plane, is mainly caused by the induced reverse flow that fills the vacancy left by the disk during the oscillation.
The rotational interaction of the fluid in the horizontal plane is relatively restricted in the limited area near the sharp edges of the disk. However, the advection of fluid in the horizontal plane affects the vortex formation and the determination of the hydrodynamic coefficients, which will be revealed in the following sections.

Added mass coefficient C a
The added mass coefficient C a of D150 with a thickness of 0.254 m obtained in the 3D simulation is presented in Fig. 7 with the counterpart 2D simulation presented for comparison. Generally, the 3D results C a3D are larger than the 2D result C a2D , with two clear increasing regimes: the linear increasing regime and the constant increasing regime. In the linear increasing regime, the 2D and 3D results increase linearly with the increase of the oscillation amplitude KC. In the constant increasing regime, the variation between the 2D and 3D results remains constant regardless of the variation of the oscillation amplitude KC.
The presence of the two increasing regimes with different slopes of the added mass coefficient variation may be interpreted as being influenced by the intensity of the 3D effect. In the low KC range (KC<0.6), when the oscillation amplitude is small, the disturbance of the flow field is restrained to a small area. The fluid can be considered a 2D flow, as results in minor difference between the 2D and 3D simulation results. With the increase of the oscillation amplitude, the influence of the flow in the heave direction enhances, and the lateral flow begins affecting and being affected by the local flow. This mutual influence is strengthened with the increase in the oscillation amplitude KC, which results in the linear increasing regime of the added mass coefficient C a . It should be noted that the increase of the mutual influence is accompanied by the increase in flow turbulence in the oscillation direction, as is a dominant factor that determines the added mass coefficient C a . In the high KC range, the flow turbulence in the oscillation is more intense than the stabilizing mutual influence in the lateral flow. As a result, the increase of added mass coefficient from 2D to 3D, which reflects the 3D effect that is mainly caused by the mutual influence, stabilizes as a constant value.
The existence of two different increasing regimes (linear and constant increasing regimes) indicates the existence of a turning point, which is a specific oscillation amplitude KC that divides the two increasing regimes. This indication is confirmed by the quantitative analysis of the difference between 2D and 3D results. Fig. 9 shows the increment rate from the 2D results to the 3D ones. Two obvious regimes can be indentified, a linear increasing regime in the low KC range and a constant increasing ratio with the oscillation amplitude larger than the specific oscillation amplitude KC, i.e., the turning point.
In the 2D simulations, the added mass coefficient depends linearly on the oscillation amplitude KC in the entire investigated KC range, which is different from the 3D simulations. With the two increasing regimes, two linear depend-ent ranges with different slope rates are observed in the 3D results in Fig. 7. Fig. 8 demonstrates the added mass coefficients of the investigated configurations with a thickness of 0.254 m and different diameter ratios obtained through both 2D and 3D simulations. The magenta dotted line is the connecting line of the turning points, which divide two regimes. The emergence of the turning point varies with the increase in the diameter ratio D d /D c and exhibits a nonlinear trend. In addition to the increase of the added mass coefficient C a , two KC-dependent increasing regimes are also observed in all demonstrated results.
The turning points of the configurations with a thickness of 0.254 m and different diameter ratios of D d /D c are shown in Fig. 9, which is the result of quantitative analyses on the disparities between 2D and 3D simulations. The percentages shown in Fig. 9 are the augments from the 2D results to the 3D results as a percentage of the 2D results at the corresponding KC. Quantitative analyses based on different diameter ratios D d /D c reveal two distinct linear and constant increasing regimes. The turning point, which is the specific oscillation amplitude KC that separates the two increasing regimes, shows a nonlinear correlation with the diameter ratio D d /D c , whereas the constant increasing ratios appear to converge to a fixed value of 10%, as shown in Fig. 9 for an increase of the diameter ratio D d /D c .   Fig. 9. Quantitative analysis of the variation of the added mass coefficient C a between the 2D and 3D simulated results.
To systematically investigate the dependence of the turning point on the thickness ratio t d /D c and diameter ratio D d /D c , numerical simulations were carried out for the scenarios in Table 1 with all presented thickness values. The results are shown in Fig. 10. It can be observed from Fig. 10 that with the same diameter ratio D d /D c , the turning point tends to stabilize to a fixed value as the thickness ratio t d /D c increases. The thickness ratio t d /D c is the primary factor that influences the turning point for each specified diameter ratio D d /D c . Furthermore, the turning point tends to converge with the decrease of the diameter ratio D d /D c . Although the occurrence of the turning point varies with the diameter ratio D d /D c , the augment at the turning point, tends to stabilize to a fixed value as demonstrated in Fig. 11. Based on the investigated models, the augments from the 2D results to the 3D results tended to be stable when D d /D c >1.5, and this stabilization tended to converge as the thickness ratio t d /D c increased. For the investigated configurations of D d /D c >1.5, an increase of 10% from the 2D results is a conservative approximation of the 3D results for the constant augment.

Damping ratio Z
Generally, the damping ratios obtained through the 3D model are higher than those obtained through the 2D model, as demonstrated in Fig. 12. The results shown in Fig. 12 are simulated by the model with thickness of 0.254 m and have different diameter ratios D d /D c . The distinct increasing regimes observed for the added mass coefficient are not obvious in Fig. 12, although it appears in the 3D results.
The quantitative analysis of the variation between 2D and 3D damping ratios reveals that the regime division is similar to those in the added mass coefficients, as shown in Fig. 13. The turning points that separate the regimes are consistent with those of the added mass coefficients. However, the factors or modes are entirely different from those of the added mass coefficients. Two distinct regimes are observed in Fig. 13, the linear increasing regime in the low KC range and the "constant" regime in the high KC range. It is similar to the added mass coefficient. However, the "constant" increasing regimes in Fig. 13 show some downtrend, which might be described as the pseudo-constant regime. The configuration of the damping ratio with   the lower KC value is in contrast with the added mass coefficient. Unlike the various increasing slopes in the linear increasing regime of added mass coefficients, the increasing slope of the damping ratio in the linear increasing regime appears to be constant. This implies the key role of the thickness ratio t d /D c in the determination of the increasing slope.
The quantitative analysis of the damping ratio variation on all models results in Table 1 is summarized in Fig. 14. The occurrence of the turning point of the damping ratio is identical to that of the added mass coefficients. The dependences of the damping ratio on the diameter ratio D d /D c and thickness ratio t d /D c show different patterns from those of the added mass coefficient (see Fig. 11).The damping ratio tends to converge with the decrease in either the thickness ratio t d /D c , or the diameter ratio D d /D c , whereas the added mass coefficient augment ratio tends to converge with the increase in the thickness ratio t d /D c or diameter ratio D d /D c .
Although there is an augment ratio of the damping ratio from the 2D results to the 3D results, as shown in Fig. 13 and Fig. 14, the absolute increase of the results is small. This small variation of 2D and 3D damping ratios is mainly due to the scaling factor between the drag coefficient and damping ratio, which cancel out the 3D effect. However, the discussion of the 3D effect on the drag coefficient or the linear damping coefficient is insignificant because the damping ratio z is required in practical application.
The 2D approximation of the damping ratio is adequate for the simplification if vertical oscillation is the only factor that causes the flow field changes. However, to reveal the 3D effect on the flow regime and added mass coefficient, the 2D approximation is imprecise. The proper approximation of the augment ratio would be a practical simplification.

Conclusions
Numerical studies on the vertical oscillation of a cylinder with a bottom-attached disk were conducted with different diameter ratios D d /D c and thickness ratios t d /D c . A numerical model has been developed for these studies with high-order upwind scheme applied to stabilize the computation. The convergence is accelerated by the multi-grid method. Qualitative and quantitative analyses on the variation between 2D and 3D results reveal the 3D effect of the flow influenced by the vertical oscillation of the structure.
The 3D effect on the fluid field is mainly reflected by the significance of three vortex-shedding patterns, shows the large affection flow field around the sharp edges of the disk compared with the vortex generated in the 2D simulation. In the slice along the axial orientation, the effect of along the radial direction is smaller than along the circumferential direction, i.e., the radial effect on the fluid velocity is larger than the circumferential one around the sharp ω z edges of the disk. The rotational interaction of the fluid in the horizontal plane during the heave motion is insignificant. However, the advection interaction of fluid in the horizontal plane affects the formation of the vortex and determines the hydrodynamic coefficients.
Quantitatively, the 3D effect on hydrodynamic coefficients is mainly reflected in the increase of the added mass coefficient obtained through 3D numerical simulation compared with that obtained through 2D simulation. Two types of increasing regimes are revealed, the linear increasing and constant increasing regime. The regime is separated by the turning point, which is related to the oscillation amplitude KC. It is affected by the diameter ratio D d /D c and thickness ratio t d /D c . The linear increasing trend is observed in the linear increasing regime of the added mass coefficient. When KC is larger than the turning point, the augment ratio of the added mass coefficient is constant and only depends on the configuration. The turning point depends on the configura-  tion and tends to converge with the increase in the thickness ratio t d /D c and diameter ratio D d /D c . The augment ratio of the added mass coefficient tends to converge with the increase in either the thickness ratio t d /D c or the diameter ratio D d /D c . In the constant regime, an approximation of 10% increase of the 2D results is a practical simplification.
The damping ratio of the 3D effect is relative small. The increase of 1% from the 2D results is a conservative approximation. Considering the minor impact of the 3D effect on the damping ratio, the 2D results are an acceptable approximation for the simplification. Turning points exist for trends of the added mass coefficient and damping ratios in all investigated models, the two hydrodynamic coefficients have different augment modes. The augment ratio of the damping ratio shows a converging tendency with the decrease in either the thickness ratio t d /D c or the diameter ratio D d /D c .
The entirely different dependence of the augment ratio indicates the key role of the cylinder wall in the disturbance of the flow pattern, when the diameter ratio D d /D c approaches one, and this requires further investigation.