Evaluation of the Viscous Drag for a Domed Cylindrical Moored Wave Energy Converter

Viscous drag, nonlinear in nature, is an important aspect of the fluid–structure interaction modelling and is usually not taken into account when the fluid is assumed to be inviscid. Potential flow solvers can competently compute radiation damping, which is related to the radiated wave field. However, the drag damping primarily related to the viscous effects is usually neglected in the radiation/diffraction problems solved by the boundary element method (BEM), also known as the boundary integral element method (BIEM). This drag force can have a significant impact in the case of structures extending much deeper below the free surface, or for those that are completely submerged. In this paper, the drag coefficient Cd was quantified for the heave and surge response of a structure which consists of a moored horizontally oriented domed cylinder with two surface piercing square columns located at the top surface. The domed cylinder is the primary part and is submerged. The drag coefficient is estimated using the experimental measurements related to harmonic monochromatic wave–structure interaction. Finally, this estimated drag coefficient was used in the modified time domain model, which includes the nonlinear viscous correction term, and the resulting device response in heave and surge directions is presented for an irregular incoming wave field. The comparison of the numerical model and the experiments validates the estimated Cd values obtained earlier. Prior to the time domain model, frequency-dependent parameters such as added mass, radiation damping, and excitation force were computed using three mainstream potential flow packages (that is, ANSYS AQWA, WAMIT, and NEMOH), and a comparison is presented. The effect of free surface on the drag coefficient is investigated through differences in Cd values between heave and surge modes.


Introduction and Background
A part of the resistance to structures moving in fluid, for example, marine energy devices, is caused by the fluid viscosity. Displaced structures set up waves and generate eddies. Resistance due to waves is sometimes termed as gravity resistance as the waves build up potential energy. However, eddy formation and their resistance (or dissipation) is more related to the viscosity of the fluid. Therefore, fully submerged structures (such as submarines) experience relatively significant viscous resistance, while the generated waves and eddies ultimately dissipate due to the viscosity of the fluid.
The seakeeping and/or manoeuvring time domain mathematical model [1] relies on the frequency domain hydrodynamic parameters, which are computed through linear potential theory, where fluid is assumed nonviscous (inviscid).
Inviscid potential flow solvers provide robust solution of the wave-structure interaction response for a wave energy device and in order to account for the viscous effects, the use of the Morison equation has been reported by a number of studies [2][3][4]. However, identification of the drag coefficient used in this approach is quite challenging. A value for this coefficient can be found empirically or using numerical modelling, such as the RANS (Reynolds-averaged Navier-Stokes) solvers [5][6][7][8][9].
The authors of [10] examined two-dimensional nonlinear time-domain solutions of the viscous effects for a twin rectangular hull under heave oscillation and compared results with existing linear, inviscid frequency domain data. In the case for low frequencies, the viscous effects are reported to be significant even for small amplitudes of the body oscillation. Estimate of viscous effects is crucial for an accurate numerical modelling of the wave energy structures. In order to include the nonlinear drag force in such numerical modelling, the CFD (computational fluid dynamics) solvers can be employed to deduce the additional damping coefficients for a floating wave energy structure. In Reference [8], the viscous damping of a three-dimensional model of the float-over barge was estimated using a RANS (Reynolds-averaged Navier-Stokes) solver and a comparison with the experimental data showed applicability of the CFD for evaluation of the viscous damping associated with the low frequency waves.
The viscous drag's impact on the power production of a generic surging type wave energy converter was studied using a nonlinear time domain model with viscous correction [2,7]. In relation to the viscous drag, the power output for a heaving wave energy converter was analysed in Reference [6], where it was reported that the viscous drag force caused a 4% reduction in the annual power production of such a point absorber device. A similar study for a generic surging type wave energy converter reported that the numerical power prediction can result in significant overestimation (more than 100%) if viscous effects are not taken into account [2]. CFD simulations can provide drag values (coefficient) that can be used in the time domain dynamic numerical model, as reported by References [6,7,11].
In addition to the viscous drag, the incoming wave nonlinearity, moorings scopes, and the power take off offer added nonlinear contributions, as has been reviewed in Reference [3]. In Reference [4], CFD modelling of the viscous effects was reported for a generic model scale buoy moored in waves, and issues and limitations associated with the CFD analysis were also discussed, such as being computationally demanding for long term power prediction. Recently, the authors of [12] reported a study of the nonlinear forces on a cylindrical generic moored point-absorbing wave energy converter using a RANS solver, the Euler equations solution, and the linear radiation-diffraction method. Viscous effects in relation to device scale are discussed. The scale effects and viscous drag relationship is usually seen as drag verses Reynolds number problem; however, for oscillatory flow, this can be a challenging task, in particular when RANS solvers are used, because results can be very sensitive to the mesh structure and the turbulence model used. The viscous force can have a significant impact on the control strategy considered. For a heaving wave follower type wave energy device, the identification of the drag using an extended Kalman filter was reported in Reference [13]; the study showed the importance of the viscous correction in relation to the control system performance.
In relation to multidegree of freedom motion, in Reference [14], the heading stability of a turret-moored ship was assessed by varying damping in different modes of motion. The study showed that viscous effects in pitch mode are stronger than heave or roll. However, for wave energy converters, there is a need to investigate the variation of the drag coefficient in relation to the multidegree of freedom response. In this paper, the drag coefficient was estimated for a moored wave energy converter, and a comparison of the nonlinear time domain numerical model against experimental data is presented. This study shows that small-scale monochromatic tests data can be used to identify the drag coefficient needed for the Morison equation's nonlinear viscous term. Further, this paper shows that for a device that has an irregular shape, there can be considerable difference in the identified drag coefficient in relation to the motion mode (such as heave and surge response).

Introduction
The significance of the viscous drag part of the hydrodynamic interaction that a marine structure experiences is case-dependent, due to a varying encounter (wave) frequency and the response amplitudes, which are usually classified in the form of the Reynolds number (Re) and the Keulegan-Carpenter (KC) number (Equation (2)). By definition, Re is given by Equation (1), where the numerator is the inertial part and the denominator is the viscous contribution: where, ρ = fluid density, -U = maximum velocity, -D = relevant characteristic dimension of the rigid structure, µ = dynamic viscosity of the fluid.
For the oscillatory flow phenomenon, the problem becomes more complex as the drag force also depends on another nondimensional measure of the KC number (for example, see [15]): with: -U m = amplitude of the sinusoidal velocity, -T = time period of the sinusoidal velocity.
For sinusoidal motion ( or for deep water waves), the KC number can be rewritten as ( [15]): where a is the amplitude of the wave or of structure displacement. The viscous drag force evaluation for simple geometrical shapes like cylinders and/or squares (with smooth and sharp corners) has been widely studied.
The viscous parameter β is defined as the ratio of Re to the KC number. The amplitude of oscillation and the structures geometry are the relevant properties that quantify the viscous force in a particular situation.
Marine renewable devices such as wave energy converters are designed to exhibit relatively much larger motion amplitudes compared to conventional tension leg platforms (TLPs). Therefore, viscous effects are of primary importance, as their influence will have a direct effect on the energy production. Preliminary design assessment of such devices has implemented the seakeeping mathematical models where the viscous force (due to skin friction, flow separation, and eddies contribution) can be included using an additional viscous drag term [1,16]. As viscous force incorporation in this paper is achieved through the Morison equation, therefore, a brief introduction of the Morison equation is presented next.
The Morison equation ( [17]) provides a semi-empirical formulation to model the unsteady force on rigid structures in oscillatory flow. According to this expression, the total in-line force on an immersed object within an oscillatory viscous fluid with velocityẊ and accelerationẌ is expressed as a summation of two components: force due to the inertia; an effect of the irrotational (potential ) assumption, i.e., ρVC IẌ ; -force due to the viscous drag; effect of the skin friction and flow separation, i.e., 1 2 ρA r C dẊ |Ẋ | For a three-dimensional structure moving in-line with the oscillatory flow, the Morison equation becomes [15]: where t is the time index and: The device is composed of two primary parts: (i). a submerged horizontal cylindrical part which has domed ends; (ii). two rectangular columns attached at the top surface of the domed cylinder.
The rectangular columns are in surface piercing mode, and the draft depth is such that only half the length of these columns is submerged. The device is moored using two vertical mooring lines attached to the submerged cylinder and a free floating clump mass. The clump mass is then connected to the bottom (floor) of the tank using two additional horizontal mooring lines. For visual description of this setup, a picture of the meshed geometry and the schematic showing device dimensions are presented in Figure 1.
Frequency-dependent hydrodynamic properties are computed using a radiation-diffraction boundary element method solver. Three state of the art solvers-WAMIT (Software package of WAMIT Inc. (http://www.wamit.com/)), NEMOH (Open source code from LHEEA lab, Ecole Centrale de Nantes, France, version 2.0.) and AQWA (Software product of ANSYS Inc. (http://www.ansys.com/) version14.5.))-were employed for this purpose, and a comparison of these packages is then reported.
NEMOH is an open source code that is based on Aquaplus, which is a BEM solver for calculation of first-order hydrodynamic coefficients.
Utilising the frequency domain data of NEMOH, a time domain model was used to compute the time series response of the floating structure, and the results are then compared against the experimental measurements. The body geometry and input parameters, such as mass properties, stiffness, and damping, are sourced from the panchromatic wave tests for configuration B of Reference [18] (see Figure 1), and these parameters are shown in Table 2. The time domain solution was derived using the convolution integral for the radiation force according to Cummins relation [19], and the viscous drag of the fluid structure interaction is included in the time domain model using nonlinear quadratic drag damping of Morison's equation [17]. Following NEMOH's results, a time domain model with viscous correction ( [2]) was then implemented. In Reference [2] such a model was referred to as a potential time domain viscous (PTDV) model, and a comparison of the viscous correction was presented against CFD computations. It was shown that for a floating wave energy converter, viscous contributions show a significant impact on the power output of the surging wave energy device. Only surge mode was reported in that study with C d values derived from curve fitting of the Morison equation and the radiation force using CFD [7].
However, in the study presented here, the structure is a small-scale model which is composed of a submerged and a semisubmerged part, and the experimental measurements data were used for the estimate of the C d values.

Frequency Domain Model And Results
Hydrodynamic parameters from frequency domain analysis using WAMIT, AQWA, and NEMOH were computed using a slightly varying mesh structure and a number of tested frequency points. These frequency domain results of the hydrodynamic parameters are shown in Figure 2a-c. Compared to RANSE solvers, mesh independence for BEM solvers of the potential theory is relatively usually not an issue unless the considered panels are extremely coarse [2]. In the modelling reported here, it was ensured that a fairly dense mesh had been achieved, for example, in NEMOH and AQWA, the body was discretised into 1629 panels. It can be seen that results from all three packages are in reasonably good agreement; however, a small (negligible) percentage of the discrepancy in these results can be attributed to the differences of mesh structures, input frequency points, and the computational resource. Furthermore, the range of the frequency steps used in this computation is considerably large for NEMOH (300 frequency steps) and WAMIT (300 frequency steps) compared to the AQWA, where only 50 frequencies were used.
In Figure 2a, it can be seen that the added mass due to surge motion is approximately twice the added mass associated with the heave motion. This shows that the mass of the water that is accelerated due to surge motion of the structure is much higher (nearly twice) compared to the heave motion. Similar, from Figure 2b, radiation damping experienced by the structure due to radiated wave field is also considerably higher compared to the one observed in heave motion. This shows that the rigid floating structure is radiating significantly more waves due to surge motion compared to the heave motion. However, for waves with a wave period higher than 1.8 s (i.e., frequency lower than 3.5 rad/s), the radiation damping drops to zero for both of the cases, i.e., the surge and the heave. It can be seen in Figure 2c that the magnitude of the complex amplitude of the wave excitation force component in surge direction is significantly higher than the corresponding component in the heave direction.
Following the computation of the frequency-dependent hydrodynamic parameters (such as added mass matrix A, radiation damping matrix B, and the excitation force vector Fe), the frequency domain response amplitude of the device is obtained from the complex magnitude of (Z(ω)) using Equation (5): where: The magnitude of thez is the response amplitude operator (RAO) as a function of the wave-frequency of unit wave-amplitude. C is the additional linear damping. For frequency domain analysis, a value for this linear damping was adopted from free decay tests of surge and heave motions shown in Reference [18] (Figure 7). A value of 3 N/m/s for heave and 7 N/m/s for the surge motion were used for the linear frequency domain results.

Time Domain Model
Time domain analysis (based on frequency domain parameters of NEMOH) relies on the equation of motion Equation (6): where: µ ∞ = added mass at infinite frequency, -Z(t) = acceleration of structure, -F ex (t) = wave excitation force , -F d (t) = non-linear viscous drag force.
Representation of the term ∑ N prony i=1 I i in a state-space model is given below, in Equation (7): where: ,N prony = complex velocity-dependent variables, α = real part of the Prony coefficient, β = imaginary part of the Prony coefficient.
In Equation (6), F d (t) is the additional quadratic drag damping term in accordance with the Morison equation. i.e.,: where V r = (Ż(t) b − V w (t)) is the relative velocity, with V w (t) being the instantaneous velocity of the incident wave.Ż(t) b is the body velocity. ρ is the fluid density. A r is the relevant cross-sectional area of the structure. The excitation force due to irregular wave signal is given by: where denotes the real part of the complex variable, and: -A wave = wave amplitude of j th frequency wave, -Fe(ω j ) = frequency-dependent wave excitation force (complex variable), ϕ j = phase variable used to generate random wave field' -N wave = total number of wave frequencies .
Here (in Equation (6)), the convolution product of radiation impulse response function times the velocity has been replaced by additional state variable I i using the approach described in Reference [20]. The use of the prony method is implied for the computation of coefficients α and β using the method shown in Reference [21]. In Equation (7), N prony is the total number of prony coefficients. Finally, the equation of motion (Equation (6)) is solved using the ode45 (based on an explicit Runge-Kutta formula) solver function in MATLAB.

Evaluation of the Drag Coefficient (C d )
The drag coefficient of the Morison equation (Equation (4)) can be determined by means of experiments or numerical computations. The process involves measurements/computations of the force time history from the experiments or numerical computations. This force series is then utilised in the curve fitting against the Morison equation yielding appropriate coefficient values based on the criteria of the curve fitting approach. For such a purpose, a number of methods can be considered, for example [22]: Morison's method, the Fourier series approach, least squares method, and weighted least squares method.
For the present case study, the experimental normalised motion response, as reported in Reference [18], was used in order to estimate the C d values, and this was achieved by curve fitting of these experimental readings against the time domain model of Equation (6) through the least squares approach.
The Morison viscous force term is proportional to the cross-sectional area of the moving structure. From the geometry of the structure (Figure 1), it can be seen that this area in heave mode is the cross-sectional area of the domed cylinder orthogonal to the heave direction; however, for the surge mode, the cylindrical cross-sectional area orthogonal to the surge direction has an additional contribution from the two individual surface-piercing rectangular columns.

The Drag Coefficient (C d )
In Figure 3, normalised heave-response is shown against experiments for the estimated C d value of 0.5, whereas for surge-response, a similar comparison is shown in Figure 4, where the estimated C d value is 2.4.
Corresponding KC and Re values can be computed using Equations (3) and (1). For the KC number, the amplitude of motion a is obtained through multiplication of the experimental values of the normalised response (say RAO n ) per unit wave amplitude (shown in Figures 3 and 4) by the corresponding wave amplitude. Considering the sinusoidal response, for the Re number, the velocity amplitude U can be obtained by the following relation (Equation (10)): where ω is the wave frequency in rad/s. The KC and Re values are shown in Table 1.

Time Domain Model Results And Analysis
The instantaneous time history of the heave and surge displacement for the given irregular random wave time series are computed using the time domain model. The time series and the corresponding spectral plot of the incoming wave is shown in Figure 5; the peak period of the random wave is approximately 2.3 s, and the significant wave height is ≈72.6 mm. It can be seen that the the resonance period of the device was within the range of the considered wave spectrum.  Table 2. Key properties in this Table 2 are the total stiffness contributions experienced by the float due to attached mooring lines [18]. For the heave mode response, an extract of the time series results from Equation (6) are shown against the measured heave response in Figure 6. A similar comparison for the surge case is shown in Figure 7. The root mean squared error (RMSE) of the computed and experimental results can be computed using Equation (11). RMSE values of time series results are shown in Table 3.
Here, x j refers to the computational results and y j to the experimental measurements, j and N are the index and the total number of data points, respectively.
A comparison between the experiment data and the computed response using the time domain model (Equation (6)) was also done in terms of the spectral densities. For this, the frequency components of the fast Fourier transform were used to obtain the spectral density using Equation (12). In order to eliminate the initial conditions influence, from the 512s of time series data, results from 256 s onward were used for this analysis.
S(ω i ) = (y.yNdt)/π, (12) where y is the fast Fourier transform of the input signal, y is the complex conjugate of y, N is the number of data points, and dt the sampling time. The corresponding frequency in this case is given by Equation (13): The spectral density of heave is shown in Figure 8, where the results show a good match with the experiment; however, it can be seen that the time domain model carries some overprediction, which is shown at two frequency points. For the surge case, the spectral density is presented in Figure 9. Again, it is seen that time domain model response, although in good level of accuracy, still carries some minor overprediction, as is shown by some spikes at a few frequencies.    From the irregular wave response, the RAO or normalised response RN shown in Figures 10  and 11 was obtained via Equation (14): -S r = Spectral density of the response, -S w = Spectral density of the wave signal. Figure 10 shows the heave response of the device in relation to the unit wave amplitude at various wave periods. Here, experiments from the regular wave of varying amplitudes are plotted alongside computational results of the frequency and time domain analysis using NEMOH's hydrodynamic parameters. For this graph, the resulting response is divided by the wave amplitude in order to produce the normalised position. Figure 11 shows the frequency domain response amplitude of the surge motion of the device alongside the normalised response from the panchromatic waves and measured data from experiments for monochromatic waves of fixed amplitudes. It is seen that the device response to the panchromatic wave force signal is in accordance with the trend predicted by the frequency domain approach and the experiments.
The RMSE of the time history responses (surge and heave) is shown in Table 3. Table 3. Root mean squared error (RMSE) of the two time series comparison; the computational and the measured experimental data (presented in Figures 6-9).

Discussion
Following the RMSE shown in Table 3 , it is noted that incorporation of the viscous force produced very good agreement with the experiments; however, the heave response of the device is relatively in better agreement with the experiments compared to the surge response's RMSE. From analysis of the spectral curves shown in Figures 8 and 9, it can be noted that for the surge case, experimental results carry significant contributions for frequencies below 2.5 rad/s as well as for frequencies above 3.5 rad/s compared to the heave case (Figure 8).
For a measure of correlation between two signals (or series) (namely x j and y j , for i = 1,2,3,...,n), the correlation coefficient r can be computed using Equation (15): wherex andȳ are the mean values of x and y, respectively, and n is the total number of data points. The value of the correlation coefficient between the prediction and the experimental data is shown in Table 4. As the value of this correlation coefficient is above 0.9, this shows that the results are in excellent agreement. From Figures 10 and 11, it is worth noticing that the frequency domain approach is the most robust and offered a fairly good indication of the float's response to waves of differing frequencies per unit wave amplitude. Overall, a normalised position from the panchromatic wave field showed a quite similar trend to that of the linear frequency domain results; however, there are instances where nonlinear drag force became significant.
It can be seen that the viscous correction in the time domain model-which is based on linear frequency domain data-offered a robust method compared to more rigorous computational methods, such as Reynolds-averaged Navier-Stokes equations (RANSE) models, where computation of a monochromatic wave requires a relatively much longer computation time; however, for the quantification of the viscous drag coefficients and for the validation of the modified time domain models, CFD can be employed, as has been demonstrated [23].
One important feature observed in this modelling is the magnitude of the drag coefficient for the surge and heave motions. As shown in Table 2, the drag coefficient obtained for the heave motion is nearly 50% less than surge motion. This implies that the significance of the drag force for heave is considerably less compared to that of the surge case. The cross-sectional area or the water-plane area for the surge mode is larger than that for the heave mode due to, firstly, the additional rectangular columns and, secondly, the fact that the device is surface piercing and a part of the extra damping is, possibly, coming from the nonlinearity of the free surface interaction.
For heave motion, because the surface-piercing profile, as part of the device, is above the water plane and there is no motion against the direction of the incoming waves, there is relatively less drag offered by the free surface interaction. Furthermore, a much higher total stiffness (hydrostatic plus moorings) force can counteract part of the drag force in heave mode; hence, the resultant drag damping value is less than surge.
Finally, an analysis of the time domain model with and without an additional nonlinear viscous force (F d ) is presented. An extract of the time series response is shown in Figure 12. The case with C d = 0 refers to the scenario where F d is not included in the model. It is observed that in the absence of the additional drag force, the RMSE value (corresponding to the experiment measurements) for the heave case (where instead of nonlinear damping term, the linear damping, similar to the frequency domain model, of 3 N/m/s was employed in order to avoid unrealistic motion amplification) increases by 38% . Similarly, for C d = 0, the surge response RMSE value increases by 22.6%. This shows the relative importance of the nonlinear features that can be competently captured using drag force F d .

Conclusions
It was demonstrated in this paper that results for the frequency domain parameters computed using three radiation-diffraction panel codes-WAMIT, AQWA and NEMOH-show very good agreement.
Simulations for an irregular (panchromatic) wave were performed using the nonlinear time domain model, and a comparison against the experimental results was presented. The statistical measures, such as RMSE and correlation coefficient of the computed and the experimental data, were reported, which demonstrated a very good agreement. Thus, nonlinear features of wave structure interactions can be taken into account through the additional nonlinear viscous drag damping in accordance with the Morison equation.
Preliminary hydrodynamic evaluation using the frequency domain linear approach was shown to be reasonably accurate, and this produced a robust estimate of the structure's response for a range of wave frequencies. However, from the normalised response of time series results shown against the RAO of the linear frequency domain model, it was observed that there are instances where the nonlinear viscous effects were significant. Consequently, it is concluded that viscous effects can cause a considerable influence on the resultant response.
For the particular device, it was shown that the magnitude of the estimated drag coefficient C d for the heave motion is ≈5 times less than the corresponding C d value of the surge motion.