Drag, added mass and radiation damping of oscillating vertical cylindrical bodies in heave and surge in still water

Forced heave and surge motion of axisymmetric vertical cylindrical bodies with flat and rounded bases are simulated using the advanced CFD software STAR-CCM+ where the overset method is used so that the mesh local to the body moves within a stationary outer mesh. Viscous effects generating drag, and also influencing added mass and radiation damping, are determined. The Reynolds-Averaged Navier–Stokes (RANS) equations are adopted with different turbulence closure models and the water surface is captured by the volume of fluid (VOF) method. These results are of basic interest but the main motive is to assess appropriate drag coefficients for use with linear diffraction models of wave energy converters (WEC) and we have particular interest in the multi-body WEC M4. A basic dynamical model is set up so that drag, added mass and radiation damping coefficients may be obtained from the CFD results. Added mass and radiation damping coefficientswerealsoobtainedfromthepotentialflowsolverWAMITforcomparison.Mesh convergence studies were undertaken and while mesh independence was achieved for total force it was not possible for the very small shear force. In the laboratory a free heave decay test was undertaken without mechanical contact for bodies with rounded bases and the inferred drag and added mass coefficients were very close to those from CFD. Some general observations are possible for motion in heave. For the hemispherical base the drag coefficient C d is very low for small amplitudes but this increases as amplitude is increased. For the rounded base with a flat central area the C d is larger and for the wholly flat base it is larger again but with values less than 0.35. For a larger geometric scale (times 32) for the hemisphere and round base cases the C d are generally somewhat reduced. For surge motion the C d show less variation and are always greater than heave values by at least a factor of 2 which is indicative of effects due to separation and wake generation. (http://creativecommons.org/licenses/by/4.0/).


Introduction
Floating cylindrical bodies are widely used in offshore engineering, e.g. spar, tension leg and other semi-submersible platforms for oil production and point absorbers for wave energy conversion (Drew et al., 0000;Falcão, 2010;Babarit, 2015). We have particular interest in the multi-body WEC M4 composed of several connected cylindrical bodies (Stansby et al., 2015b, a;EatockTaylor et al., 2016;Stansby et al., 2016;Sun et al., 2016Sun et al., , 2017Stansby et al., 2017). Linear diffraction analysis in the frequency and time domains is quite standard practice for operational conditions and extreme conditions require model testing (Drew et al., 0000;Falcão, 2010;Santo et al., 2017) or full CFD simulation for nonlinear conditions. Linear (and second-order) diffraction analysis can be undertaken with standard modern computing while CFD requires massive parallel processing with numerical convergence, and hence accuracy, uncertain. It is thus attractive to add viscous effects to the linear diffraction analysis based on potential flow theory (Stansby et al., 2015b;EatockTaylor et al., 2016;Stansby et al., 2016;Sun et al., 2016Sun et al., , 2017Stansby et al., 2017). This may be achieved through drag coefficients with body motion relative to the fluid defining an additional drag force. Viscous effects also have a small influence on added mass and possibly radiation damping and this is a secondary consideration.
In spite of the widespread use of cylindrical bodies offshore there is little information on drag coefficients, at least in the open literature. While CFD is computationally demanding or even possibly unattainable for general configurations undergoing multi-mode motion in waves, which are generally irregular and multi-directional, high resolution CFD is possible with simple geometries in single modes of motion. The aim is to determine drag coefficients to enable drag forces to be added to linear diffraction analyses. This strictly requires knowledge of the body motion relative to the fluid motion, due to incident waves with additional diffraction and radiation effects in these cases. This is difficult (or impossible) to represent as a single velocity and we consider only still water cases. This might appear consistent with the linear diffraction analysis for moving bodies where the forces due to fixed and moving body are superimposed. There is no evidence to suggest that viscous effects may be simply superimposed in this way although for oscillating bodies in line with steady flow (defined by a single velocity) this can give accurate results, with different drag coefficients for steady and oscillatory components (Verley and Moe, 1979). This is unlikely to be possible in waves with a spatially varying velocity field without an obvious representative single velocity. The aim of this study is thus to determine drag coefficients for cylindrical bodies with flat and rounded bases for heave and surge motion in still water. This is of fundamental interest and informs suitable values to be included in linear diffraction analysis; these values are basically tuning parameters accounting for complex boundary layer and wake effects. We determine the prominent component for viscous interaction, that due to body motion.
In Section 2 the dynamical equation is described which enables viscous effects to be determined from the CFD results. This also enables some experimental free decay tests to be analysed. In Section 3 the CFD methodology including mesh configuration is described. Results for forced heave and surge motion follow in Section 4 including comparison with coefficients from the free decay tests. Results for near full-scale conditions are included. There is a discussion in Section 5 and conclusions are drawn in Section 6.

Dynamical equation
The forces on an oscillating body may be categorised as added mass in phase with acceleration, radiation damping and drag both in phase with velocity and buoyancy in phase with vertical displacement if the body motion is small, i.e. the problem is assumed linear. Nonlinear effects associated with free surface movement may become significant as amplitudes increase. If viscous effects are negligible, which is often assumed for large bodies, potential flow may be assumed and linear radiation/diffraction panel methods may be applied for small motions. Here we use WAMIT version 6.3 (WAMIT, 2004) to evaluate the associated added mass and radiation damping forces. In general for sinusoidal body motion forces may be defined by constant added mass, radiation damping and drag coefficients to a close approximation. For vertical sinusoidal oscillation of a body with immersed volume V in still water level, with frequency ω and amplitude z 0 , displacement z (positive upwards) is defined by: giving body velocity and acceleration as: If F is the total body force with static buoyancy ρgV subtracted, equal to the mass of the body where ρ is water density and g is gravitational acceleration, according to Newton's law F = mz, where m = ρV . This fluid force on the body may be represented using added mass, radiation damping and drag coefficients, C a , C rd , and C d , with a buoyancy or hydrostatic stiffness term. With a linearized drag term, e.g. Sarpkaya and Isaacson (1981), where A is cross sectional area at the water plane (also frontal area normal to relative velocity in drag term), andż rms = z 0 ω √ 2/2. To determine C a the hydrostatic stiffness term is subtracted from F, giving F ′ , and the component in phase with acceleration (and displacement) is determined from To determine C d the radiation damping (assumed to be due to potential flow only) is subtracted from F, giving new F ′ where C rd = ρωC WAMIT (C WAMIT is the WAMIT output) and components in phase with velocity are determined from Hence C a and C d may be obtained from Eqs. (6) and (8) if the moving body force is known from either CFD or experiment. Since it is assumed that radiation damping is unaffected by viscosity any viscous influence is absorbed in the drag coefficient. Whilst C rd is calculated from the radiation damping coefficient of WAMIT, the equivalent drag coefficient C drd (for comparison with C d ) is given by: On the other hand if the coefficients are known substituting the moving body force = ρVz, into Eq. (4), gives the following which can be solved numerically. While the analysis above is for heave the analysis for surge, with x replacing z, is identical.
In this study drag and added mass coefficients are determined directly from CFD for motion in heave and surge. Analysis of free decay tests is also undertaken for motion in heave with drag and added mass coefficients determined by curve fitting the response to experimental measurements. Direct response prediction by CFD is also undertaken.

The CFD model
The commercial software STAR-CCM+ is applied. For the free surface flow of interest here the two-phase (water-air) solver is required. The governing equations are the Reynolds-Averaged Navier-Stokes (RANS) equations with a turbulence closure model; the k-ω SST model, the standard k-ε low-Re model and the V2F k-ε turbulence model are assessed. Low Re wall boundary treatment is applied with y + (u * y ν ) ∼ 1 for adjacent mesh size, where u * is friction velocity, y is normal distance from wall and ν is kinematic viscosity. The VOF method is used to capture the free surface. Water density is set to 1000 kg/m 3 , and air density is set to 1.205 kg/m 3 . Dynamic viscosity is 0.001 Pa s and 1.82 × 10 −5 Pa s for water and air respectively. To assess the influence of air its density was varied between 1.205 kg/m 3 and 0.05 kg/m 3 showing negligible influence on force which is thus hydrodynamic. The overset method is used with a local mesh fixed relative to the body and a stationary outer mesh connected by the overset interface where velocity and pressure are interpolated. No effect of the interface was discernible. For heave motion an axisymmetric mesh is used reducing the problem to two dimensions to save computational time. The vertical centreline plane of the computational domain is shown in Fig. 1; this is a wedge shape with periodic boundary conditions in the azimuthal direction. The pressure on the upper (air) surface is set to zero and the base and outer boundary are free slip with zero pressure gradient. To absorb radiated waves at the outer boundary a damping zone of one wavelength width is employed where a resistance term is added to the momentum equation for vertical velocity: in which w is the vertical velocity, x sd and x ed are the start and end points of the wave damping zone respectively; other coefficients can be referred to the manual of STAR-CCM+. Fig. 2 shows the mesh used for the hemisphere base case and similar meshes were used for other cases. An outer mesh, uniform in the radial direction, is applied and a non-uniform inner mesh with mesh size increasing from the body surface with y + ∼1 to the overset region where the mesh size is equal to the outer mesh. Between 40 000 and 72 000 cells were required to give convergence of the total force. For surge motion a 3D domain is set up with the same inner mesh dimensions but now the outer domain is a cuboid with equal dimensions in the horizontal directions.

Free decay tests in heave
First sensitivity to turbulence model is tested based on the hemispherical base configuration used in the WEC Wakes project (Stratigaki et al., 2014), as shown in Fig. 3. Mechanical friction was found to contaminate free-decay displacement results but sensitivity to turbulence model was undertaken for this configuration (considered again later) without comparison to experiment. The motion is obtained from F = ρVz where F is obtained from CFD. Results are shown in Fig. 4 with the laminar flow as a reference, the standard low Re k-ε turbulence model, the k-ω SST model and the k-ε V2F model. Numerical convergence was obtained with a solid boundary mesh size of 0.0025 m. The results with the latter two more sophisticated models are almost coincident while laminar flow and the standard low Re k-ε model give smaller damping and larger displacements.
The widely used k-ω SST model is chosen to compare with further free decay tests of a free-floating float with position measurement without mechanical contact. These tests were undertaken in the Plymouth COAST basin following tests on the M4 wave energy machine (Sun et al., 2016Stansby et al., 2017). Floats with a hemispherical (float 1) and rounded base (float 2) were used as shown in Fig. 5 with a water depth of 0.8 m. The mass of float 1 was 9.75 kg and that of float 2 is 23.42 kg. Two initial offsets of −0.75 m and +0.076 m are set for float 1, and one offset of −0.083 m for float 2. The natural period in heave is 0.87 s for float 1 and 0.81 s for float 2. Displacement in three dimensions was measured using a Qualisys camera system and each float was released manually from its initial offset. The motion was shown to remain effectively in heave only for at least four oscillations. The motion is calculated directly by CFD with the k-ω model and radial domain size of 8 m with an outer damping region of 1 m; the mesh size adjacent to the solid boundary was again 0.0025 m. This value will be seen to be below that to provide mesh independent forces. The C d and C a giving the best fit to experimental displacement with C rd from WAMIT were also obtained from Eq. (9). Note that C a determines the natural frequency and C d the decay rate.
Figs. 6a-6c shows that agreement from this solution and CFD with experimental measurement is very close. C d and C a were also determined on a cycle-by-cycle basis for the first four cycles using Eqs. (6) and (8) and the mean amplitude for each    Table 1 to be very close with little cycle-to-cycle variation and values are very close to those obtained by curve fitting using Eq. (9) with single values of C d and C a ; differences in peak displacement are less than 5%. The C a from WAMIT was 0.38 for float 1 and 0.36 for float 2.   Table 1) and solved by Eq. (9) with best fit C a = 0.38  Table 1) and solved by Eq. (9) with best fit C a = 0.38 C d = 0.35.  Table 1) and solved by Eq. (9) with best fit C a = 0.36 C d = 0.45.

Table 2
Values of C a and C d for each cycle of forced heave oscillation with T = 1.0 s and z 0 = 0.06 m for the hemispherical base with different outer radius for wave damping: KC = 1.26; β = 9 × 10 4 ; kD =1.21; kd = 3.1.

Forced heave motion
The comparisons with free decay tests indicate that the k-ω turbulence model performs well. Forced sinusoidal oscillations are computed with amplitudes of 0.03 m and 0.06 m and a period of 1 s. These values are representative of those occurring in the M4 tests. The case of the flat base cylinder (float 3 in Fig. 5) is also considered since this is common for some offshore applications. The aim is to assess magnitudes of pressure force and shear force and sensitivity to outer boundary radius and mesh size. Results for pressure and shear force for the hemispherical base (Float 1) with different mesh resolutions are shown in Fig. 7 for motion in heave with amplitude of 0.03 m. The outer radius R is 4.8 m (approximately 3 wavelengths). It can be seen that total force is mesh independent but the very small shear force is only approaching independence. Similar behaviour is observed with amplitude of 0.06 m in Fig. 8 with shear force appearing closer to convergence. Although the shear force is small it is predominantly in phase with velocity and needs to be compared with the total drag force in phase with velocity. First we consider the influence of the outer radius, tested for the hemispherical base undergoing motion in heave with amplitude of 0.06 m and period 1.0 s.
A constant absorption zone is adopted with 1.56 m. The effect of outer radius on the coefficients for a period 1.0 s for several cases is listed in Table 2. Here C a is about 0.45, close to the value of 0.43 from WAMIT. C d decreases with increasing outer radius, When R/L (outer radius/wave length) is greater than 3 (4.8 m with a period of 1 s) C d changes are negligible. Fig. 9 shows the variation with C d calculated from total force to be about 30% greater than that calculated from pressure force. Table 3 shows C d and C a for shear, pressure and total force for the three shapes; values of added mass and radiation damping from WAMIT may be compared. There are some trends. For the hemisphere base C d is small with amplitude of 0.03 m, about 1% compared to that due to radiation damping (C drd ) and the shear component is larger than the pressure component, but this increases to about 6% with an amplitude of 0.06 m, and the pressure component is now larger than the shear component. There is little information in the literature on the shear and pressure forces but for the case of a circular cylinder in oscillatory cross flow (pre-separation) these two components may be shown analytically to be equal with laminar or turbulent flow (Cobbin et al., 1995). For the round base the pressure component is larger than the shear component and the drag coefficient is smaller for the larger amplitude (in contrast to the hemisphere base). However the C d is about 10% of the equivalent for radiation damping C drd for the smaller amplitude and about 14% for the larger amplitude. For the flat base the C d is similar for both amplitudes and the pressure component is much greater than the shear component, C d is now about 40% of the equivalent for radiation damping. These changes will be related to the size of separation zones, negligible for the hemisphere and significant for the flat base. Extending the amplitude range for the flat base case it is shown in Table 3 that C d changes little with amplitude up to 0.1 m.
The flow is defined generally by the non-dimensional parameters Keulegan Carpenter number KC = 2π z 0 /D where z 0 is amplitude and D is diameter of a float, the period parameter β = D 2 /νT where ν is kinematic viscosity and T is period (convenient alternative to Reynolds number), the depth parameter kd where k is wave number and d is water depth, and diameter to wavelength parameter kD. For kd >3 waves are deep water and this parameter has no influence.

Forced surge motion
Forced surge motion is investigated as for heave motion with the three floats of Fig. 5 with amplitudes of 0.03 m and 0.06 m and a period of 1 s. Mesh settings are the same as for the heave simulations and non-uniform for overset mesh with a fine grid with y + ∼1, and with uniform radial grid for outer region. The outer domain is 9.6 m × 9.6 m in horizontal plane and 4 m vertically with a water depth of 2 m.
Results for C d and C a are shown in Table 4 as for heave motion. A clear difference is that C d is always significant with a minimum of 10% of equivalent radiation damping C drd for the hemisphere base with an amplitude of 0.03 m. For this case the shear force is 13% of the pressure force and in other cases is a smaller proportion. A surprising result is that C a is always slightly reduced by viscosity in contrast to generally being increased in heave motion. The larger C d values are indicative of separated flow.

Forced heave and surge motion at full scale
The model scale is considered to be 1:32 so a 1 s period represents 5.66 s at full scale, applying Froude scaling. This is just one case and a range of values of period would occur at full-scale. However the intention is to assess the effect of a large realistic scale principally on drag coefficient. Results equivalent to those at laboratory scale are shown in Table 5 for heave and Table 6 for surge. Mesh settings are adjusted from those for the small-scale test, with minimum cell size of about 0.00001 m giving y + <5. The parameters, kD, KC and kd are the same as laboratory while β is increased by a factor of 181.
For the hemisphere base the C d values are similar to laboratory scale although for the small amplitude the value is now very small. For the round base the values are again very similar at both scales. For surge motion results are shown in Table 7 and C d values are slightly lower than those at laboratory scale. A general observation is that shear force is less than at laboratory scale although values are always small in both cases.

Discussion
This study has shown that it is difficult to obtain fully converged results for drag effects associated with shear force even for simple, single modes of oscillation of heave and surge; however the shear force is a very small component of total force. The STAR CCM+ software was run on 4 to 8 processors and each run typically required 36 to 90 h. Longer runs are impractical and full parametric studies for a wide range of KC and β are also impractical. However for the prediction of free heave decay for two shapes best fit coefficients were quite accurately aligned with coefficients obtained by CFD. This suggests the drag coefficients obtained provide a useful indication for practical application. The added mass coefficient on the other hand is generally close to the WAMIT value with only a small viscous effect that may be positive or negative, depending on mode. It should also be borne in mind that any effect of viscosity on radiation damping is absorbed into the drag coefficient although this is expected to be relatively small as it was for added mass.
Some trends are significant. For heave motion C d is very low with the hemispherical base for small amplitudes but this increases as amplitude (KC ) increases. For the rounded base the C d is larger and for the flat base it is larger again but with values less than 0.35. This may be thought surprising for an almost square edge. In steady flow for example a square cylinder has a C d of about 2 and in oscillatory flow with slightly round corners was also close to 2 (Smith and Stansby, 1991). However the latter was for laminar flow with β = 172. This case was run using STAR CCM+ which gave a similar result. However with β = 90 000 as in the present study, C d of about 0.2 was obtained, less than that for the flat base float. For the larger geometric scale (times 32) the C d are generally slightly reduced for both the hemisphere and round base cases. For surge motion the C d shows less variation and are always greater than heave values by at least a factor of two which is indicative of some effect of separation and wake generation.
It is also of interest to consider the effect of body shape on radiation damping. For wave energy converters it is well known that float shapes that are good radiators are good converters so high radiation damping is advantageous. Values of radiation damping coefficient from WAMIT for heave are compared in Table 7 for the five shapes considered here (including the WEC Wakes case). It is clear that as the ratio of vertical side length to diameter decreases the radiation damping coefficient increases.
In relation to the performance of the multi-float WEC M4 it was shown experimentally that rounded bases gave improved energy capture markedly over flat bases (Stansby et al., 2015a). This was considered to be due to reduced drag in heave.  While significant changes to drag are observed for the different float shapes changes in radiation damping may also have an effect. Clearly the interactions are complex. Linear diffraction modelling in the time and frequency domains for floats with rounded bases gave good power and response predictions with zero drag coefficient. Interestingly the predictions improved as wave fields became more complex, as they changed from regular to irregular and irregular, multi-directional waves and this is thought to be due to wave basin reflections averaging out for the more complex cases. For the case with flat-base  floats linear diffraction modelling with drag effects defined by C d ≈ 0.4 gave the best agreement with experiment, slightly smaller than but consistent with these CFD results. This study is intended to inform values for drag coefficients in linear and second-order diffraction modelling. The actual flows will be complex and superposition of drag for different modes of oscillation is not strictly valid as it is in linear potential flow theory. Ultimately values should be calibrated against experiment, as has been done for M4, or against CFD for the complete system. However given the difficulty in obtaining mesh convergence for these simple modes and single bodies such full CFD is presently not practical.

Conclusions
Forces on cylindrical bodies with rounded and flat bases undergoing heave and surge motion have been investigated using the CFD software STAR CCM+ with the overset mesh facility and free-surface modelled with VOF method. The k-ω SST turbulence closure model was mainly applied. The aim is to determine how viscous effects determine drag coefficients and modify added mass to inform their magnitude for superimposing in linear diffraction modelling. The investigation was prompted by experimental analysis of the multi-float WEC M4 showing 60% power capture improvement by changing from flat to rounded bases and subsequent linear diffraction giving excellent predictions assuming zero drag. With forced heave oscillations numerical convergence was not achieved for the very small shear force although it was for total force. Even for these idealised motions and simple shapes computations were demanding, using parallel processing. Despite these limitations very good predictions of free decay in heave experiments was obtained for a hemispherical ended float, and another rounded base cylinder. C d and C a producing the best fit to response were in close agreement with those from CFD on a cycle-by-cycle basis.
There are some significant trends. For heave motion C d is very low with the hemispherical base for small amplitudes but this increases as amplitude (KC ) increases. For the rounded base the C d is larger and for the flat base it is larger again but with values less than 0.35. For surge motion the C d show less variation and are always greater than heave values by at least a factor of 2 which is indicative of some effect of separation and wake generation. The added mass coefficient on the other hand is generally close to the WAMIT value with a small viscous effect that may be positive in heave or negative in surge. This is consistent with the excellent predictions of energy capture and response of the multi-float WEC M4 with hemi-spherical and rounded bases by linear diffraction modelling assuming zero drag Stansby et al., 2017). Importantly at a larger geometric scale (times 32) the trends are similar although viscous effects are somewhat reduced.