VD-PQ; A Velocity-Dependent Viscous Damping Model for Wave-Structure Interaction Analysis

: For the analysis and design of coastal and offshore structures, viscous loads represent one of the most influential parameters that dominate their response. Very commonly, the potential flow theory is used for identifying the excitation wave loads, while the viscous damping loads are taken into consideration as distributed drag type loads and/or as linear and quadratic damping loads ap-proximated with the use of motion decay curves of the structure in specific degrees of freedom. In the present paper, is developed and proposed a numerical analysis method for addressing wave-structure interaction effects through a velocity-dependent viscous damping model. Results derived by a computational fluid dynamics model are coupled with a model that uses the boundary element method for the estimation of the viscous damping loads iteratively in every time-step of the analysis. The computational fluid dynamics model solves the Navier–Stokes equations considering incompressible flow, while the second model solves the modified Cummins Equation of motion of the structure in the time domain. Details about the development of the coupling method and the veloc-ity-dependent viscous damping (VD-PQ) are presented. The coupling between the different models is realized through a dynamic-link library. The proposed coupling method is applied for the case of a wave energy converter. Results derived with the use of the developed numerical analysis method are compared against experimental data and relevant numerical analysis predictions. The importance of considering the instantaneous velocity of the structure in estimating the viscous damping loads is demonstrated. The proposed numerical analysis method for estimating the viscous damping loads provides good accuracy compared to experimental data and, at the same time, low computational cost.


Introduction
Renewables represent the largest source of energy and are set to penetrate the global energy system more quickly than any fuel in human history; offshore wind turbines are contributing significantly to this growth and is expecting to dominate into renewables up to 2040 [1]. However, since the higher energy wind resource exists in deep-sea areas (approximately deeper than 70 m), the development of floating wind turbines in deep seas has been targeted. Different floating support platforms are possible for use with offshore wind turbines [2,3]. On the other hand, the wave energy converters (WECs) technology is in the pre-commercialization phase. Most of the types of the WECs are designed to operate and produce power close to their resonance in different degrees of freedom since large amplitudes in motions result in larger produced power [4]. In addition to marine renewable energy structures, the offshore oil and gas structures are still being developing at a growing rate. A huge development of coastal and offshore energy structures and systems is expected for the years to come in different energy technologies. Analysis and design methods covering the whole life-cycle range of coastal and offshore energy structures are continuously redeveloped and reassessed. Analysis and design of coastal and offshore energy structures is an equally demanding and challenging task. Accuracy of the calculated response is required at a maximum level [5], while the computational cost is very important, especially when we approach the design phase of a coastal and offshore energy structure and analysis for many operational and extreme environmental conditions is required to be performed. Unfortunately, the response calculation accuracy contradicts the computational cost (Figure 1). Wave-structure interaction effects and resulting excitation wave loads should be appropriately addressed with rationality but within the limitations of computational capacity and cost. It will be ideal if numerical methods with low computational cost and larger fidelity can be used in the design phase of marine energy structures, but with higher accuracy results. For the analysis and design of coastal and offshore energy structures, wave loads and rest hydrodynamic coefficients are commonly calculated by using potential flow theory [6]. Assumptions of inviscid fluid, incompressible and irrotational flow and small amplitude oscillations are adopted with the potential theory that makes the physical problem under study linear. The aforementioned assumptions are challenged by real sea nonlinear hydrodynamic conditions (e.g., vortex shedding, resonance, viscous drag, flow separation) that are not addressed by the potential flow theory. The calculated wave loads and rest hydrodynamic coefficients are used as input from finite element model tools for the solution of the equation of motion of a floating structure, and consequently, for its analysis and design. This method is very frequently applied in marine structure designs since low computational cost (capacity and time) is required. For the design of WECs, potential flow theory has been used and reported in [7][8][9][10]. Usually and when the potential flow theory is adopted for the analysis and design of WECs, in order to consider viscous effects during analysis, distributed drag type loads are added [11][12][13].
Alternatively, viscous damping effects for a specific motion of a floating structure can be simulated with linear and quadratic damping coefficients that can be calculated with the use of free decay curves. The free decay curves can be generated by experimental tests or by the use of a computational fluid dynamics (CFD) model. In [6] and [14], different mathematical methods for calculating the linear and quadratic coefficients are proposed with the use of decay curves. The PQ method [14] has been applied for the case of a semi-submersible wind turbine [15], while the method proposed in [6] has been used for a heaving type WEC [16], both with efficient results. It is noted that for both methods, a constant value of the linear and quadratic coefficients is used during the analysis of the floating structure in the time domain irrespectively of the velocity of the structure.
On the other hand, the modeling of the dynamics of the fluid flow with the use of CFD methods provides high fidelity accuracy and calculation of the excitation wave loads but is still considered to be relatively computationally expensive, especially when a big number of analyses for different type of environmental loadings are required. CFD analysis is the most powerful tool and a mainstream method that can be used for the analysis of different types of WECs. Up to today, CFD methods and tools have been used successfully for the analysis of different types of WECs, providing robust and reliable results. The required extensions that should be used with the CFD toolbox OpenFOAM in order to simulate efficiently oscillating wave surge converters are discussed in [17]. In [18], a CFD based numerical wave tank model has been developed for the validation study of the 1:5 scale Wavestar point-absorber WEC by including the power-take-off the system in the CFD model. The required process of coupling a high-fidelity power-take-off model with a CFD-based numerical wave tank for performing a robust wave-to-wire simulation is discussed in [19]. The validation of the numerical model of the 1:20 scale Wavestar WEC considering different test cases of increasing complexity is presented in [20]. In addition, in [21], a CFD-based model is coupled with different control strategy systems (active and passive) in order to examine the performance of a three-dimensional oscillating wave surge WEC. Many critical steps have been accomplished, but further research work should be made towards the modeling of different components of WECs with CFD-based models [22,23]. Viscous damping loads are very important for the analysis of WECs and should be accounted in the relevant analysis.
For the analysis and design of WECs, viscous loads may be accounted in wave-structure interaction problems as constant loads when the potential flow theory is adopted. In the present paper, a process for accounting viscous damping loads during the analysis of WECs in the time domain is developed and presented. The velocity-dependent viscous damping model (VD-PQ) for wave-structure interaction analysis is proposed. The model is realized through a Dynamic-link library that is coupling the different numerical models. The importance of considering the instantaneous velocity of the structure in estimating the viscous damping loads for every time-step of the analysis is demonstrated through comparisons of RAOs of motions of a WEC in two degrees of freedom against experimental data. VD-PQ is generic, provides good accuracy compared to experimental data and, at the same time, low computational cost, especially when dealing with specific types of analysis (e.g., full long-term analysis).

Description of the Numerical Analysis Method
The numerical analysis method used in the present paper is generic and can be used for any possible design of coastal and offshore floating structures. The required information that is needed as input in order for the method to be used is a decay curve of the structure for a specific degree of freedom calculated either with the use of a CFD model or with the use of relevant experimental data. CFD or experiments are considered equally reliable and accurate, and their possible uncertainties are out of the scope of the present paper. In the present paper, the decay curves of two motions of the WEC are calculated with the use of a CFD model.
The numerical analysis method of the present paper consists of: (a) a 3D numerical model based on potential theory (PTM) capable for the calculation of linear hydrodynamic coefficients (added mass, radiation damping, excitation loads) in the frequency domain, (b) a 3D computational fluid dynamics model (CFDM) capable for the calculation of decay curves of the oscillation of the floating structure in specific degrees of freedom, (c) a numerical analysis process (NAP) capable for the calculation of velocity-dependent linear and quadratic viscous damping coefficients in every time-step of the solution of the equation of motion and (d) a 3D model (FEM) capable for the solution of the equation of motion of the structure in the time domain and the estimation of the response of all components of the structure and wave field, as well. PTM and CFDM are participating statically during the numerical analysis since specific outputs derived with their use are utilized as input by the FEM and NAP. FEM and NAP are directly coupled, and in every time-step of the analysis, NAP calculates and provides the velocity-dependent viscous damping values that are used for the solution of the equation of motion of the structure. In Figure 2, an outline and coupling interconnections between the different models utilized are presented. A description of the different models used in the present paper is presented below. Initially, a PTM is used in order the linear hydrodynamic coefficients of the floating structure under study to be calculated. Those hydrodynamic coefficients are the added mass, radiation damping and excitation wave loads. All the hydrodynamic coefficients are calculated for a specific grid of the wet surface panels of the floating structure and after an appropriate grid size study. Furthermore, the hydrostatic stiffness is calculated. The flow is assumed irrotational and incompressible, and the fluid inviscid while its motion can be described with the use of the velocity potential. The velocity potential, φ, is described as below: where φ0 is the potential of the incident waves: where φD is the diffraction potential, φs is the scattered potential, φj, j=1,…,6, is the radiation potential of each rigid body motion associated with the waves that are radiated due to the forced motions of the floating structure, ω is the wave frequency, g is the gravitational acceleration and k is the wavenumber. The boundary value problem is solved based on the three-dimensional panel method utilizing Green's theorem, and appropriate boundary conditions are used on the free surface, on the sea bottom and on the floating body. Moreover, the radiation condition for the outgoing waves is adopted [24,25]. With regard to the hydrodynamic coefficients, namely, added mass, Aij, i, j=1,…,6, radiation damping, Bij, i, j=1,…,6, and wave excitation loads, Xi, i=1,…,6, those are calculated after the solution of the first-order boundary value problem as below: where ni, i=1,…,6, is the normal component of the i-th rigid body mode on the mean body wetted surface, SB, and ρ is the water density. As far as the coefficients of the hydrostatic stiffness matrix, Cij, i, j=1,…,6, those are calculated with the use of the following equation: where d is the draft of the floating structure and Di is the divergence of the motion displacement.
As far as the computational fluids dynamic model, the CFDM is developed in order to solve the Navier-Stokes equations and estimate realistic viscous damping loads in specific degrees of freedom of the structure. By applying an initial displacement on the floating structure in a specific degree of freedom and letting the structure oscillate, the decay curve of this motion is calculated with the use of the CFDM. Incompressible viscous flow is considered by using the following set of equations: where V is the velocity of the fluid, p is the pressure, μ is the viscosity, ρ is the water density, g is the gravitational acceleration, and F represents the force applied due to the presence of the structure in the computational domain. The Navier-Stokes equations are solved in the whole numerical domain on a structured grid by second-order finite volume discretization techniques in space and time. Capturing the free surface of the waves in the two-phase flow problem under study is modeled through the level set method, while water-structure interaction is modeled with the use of the immersed boundary method. The immersed boundary method has been recently adopted in offshore and ocean engineering problems [26][27][28][29][30][31] for the robust and efficient numerical modeling of physical problems related to the interaction between two fluids and moving rigid structures, unstructured grids, high-density flows, WECs operating in multiple degrees of freedom and water entry/exit problems. In [32][33][34], the slip boundary condition is imposed on all sides of the numerical domain, and no penetration boundary condition is imposed on the solid-fluid interface. It should be noted that the effect of the modeling details of the CFD model is out of the scope of the present paper. As explained previously, the CFD model should be able to provide the decay curves of the floating structure in specific degrees of freedom with a good level of accuracy.
As far as the third model, the FEM model is developed in order to solve the equation of motion of the floating structure and, consequently, calculate all relevant responses (e.g., motions, the tension of mooring lines). With the use of Newton's second law, the equation of motion of the floating body in the time-domain is as below [35]: where M is the structural mass, α is the added mass that corresponds to infinite frequency, B1 and B2 are the linear and quadratic damping coefficients of the hydrodynamic damping, C is the summation of hydrostatic stiffness and mooring lines stiffness, and X is the excitation wave load. It is noted that the hydrodynamic damping can be reproduced accurately with linear and quadratic damping terms. Moreover, ( ) ( ) ( ) are the displacement, velocity and acceleration of any translational or rotational degree of freedom of the structure and h(t-τ) is the retardation function. The calculation of the coefficients of the matrix of the retardation function and the added mass in infinite frequency is not straightforward due to the semi-infinite integral and also due to the possibility of instability for the numerical calculation of the coefficients.
With the use of the CFDM, decay curves are generated for specific degrees of freedom of the floating structure ( Figure 3). It is noted that the decay curve in Figure 3 is an example of a possible curve; the cycles of oscillations and rest characteristics depend upon the wave-structure interaction of a specific structure. With the use of the data of the decay curve and by following the methods proposed in [14] (PQ method) or in [6] linear, B1, and quadratic, B2, damping coefficients can be calculated and inserted in the equation of motion of the floating body in time-domain (Equation (8)). Based on the PQ method [14], successive positive or negative amplitudes of the decay curve in the whole decay time are considered in order to determine the relative decrement of the decay curve. The B1 and B2 coefficients are calculated with the following equations: where p and q are coefficients calculated with the use of data of the decay curve, M is the total mass of the floater (e.g., structural mass and added mass in infinite frequency), and Tn is the natural frequency of the floater for the examined degree of freedom that we want to calculate the B1 and B2 coefficients. For the estimation of the p and q coefficients, points with coordinates ( ) ( ) ( )  are calculated and plotted where Φ is the peak amplitude (positive or negative) of the decay curve. A line can be fitted through those points to calculate the intercept (p value) with the vertical axis and the slope (q value) of the line. Again, it is stated that decay curves can be generated either with the use of a CFD model or with relevant experimental tests; both types of data can be used for the estimation of the B1 and B2 coefficients. If we look carefully at Figure 3 (or any possible decay curve), we can see that at the beginning of the decay curve, the velocity of the structure is large and afterward becomes very small as the number of oscillations increases. For the case where the damping force has a large dependence on the Reynolds number or on the Keulegan-Carpenter number and for the case that the floating structure is not lightly damped, the velocity changes significantly for each subsequent oscillation. For those conditions, the use of constant B1 and B2 coefficients may be conservative and not representative for the total decay time since B1 and B2 are calculated considering the whole decay time irrespectively of the velocity of the structure.
In order to take into account the velocity of the floating structure and its effect on the viscous damping, a velocity-dependent viscous damping model (VD-PQ) is developed where Φ is the amplitude of the i th peak of the motion (positive or negative) of the decay curve, T is the relevant time that the i th peak occurs, and k is the total number of positive or negative peaks of the decay curve. Based on all the calculated values of ). It should be stressed that the t x  value is calculated with the use of Equation (8). For the two different regions, we can calculate with Equations (9) and (10)  is assigned as the linear and quadratic damping coefficients. Afterward (this is made for every time-step, t, of the analysis), an iteration is made as below in order the correct linear and quadratic damping coefficients pair to be used in every time-step. For a specific time-step, t, the linear and quadratic damping coefficients from the previous time-step, t-1, are randomly assigned. With the use of Equation (8) and after appropriate numerical integration techniques, the velocity of the floating structure x  is calculated, and its value is compared with the mean velocity value, M k,decay x  , as calculated by the decay curve (Equation (12)), as below: If Equation (13) If Equation (14) x  is compared with the M k,decay x  as below: If Equation (15) holds true, then ( ) ( ) ( ) are calculated with Equation (8) and the analysis moves to the next time-step, t+1. The iterative process that is performed at every time-step of the analysis (Equations (8)- (15)) is realized through a Dynamic-link library (DLL) and is coupled with the FEM model for solving Equation (8). In this way, a velocity-dependent viscous damping model (VD-PQ) for wave-structure interaction analysis was developed and implemented.

Results and Discussion
The numerical analysis method that was developed in the present paper was applied for a heaving type WEC. The WEC consists of two rigid bodies connected with two vertical tendons. The first body is a fully submerged horizontal cylinder with domed ends and rectangular-shaped surface piercing columns, while the second one is a fully submerged mass [36,37]. A system of mooring lines is used to keep the WEC in place by providing horizontal and heave stiffness. The predictions of the numerical analysis method were compared with relevant experimental data performed at the Kelvin Hydrodynamics Laboratory in the University of Strathclyde, Glasgow, UK [36]. It is noted that the numerical analysis method was developed for the dimensions and properties of the test model without any upscale. In Figure 4, the examined WEC in the present paper during experiments is presented. Images of the different numerical models, namely, the PTM, CFDM and FEM, that were developed and used for the purposes of the present paper are presented in Figure 5. It is noted that the details and rest characteristics of the components of WEC as presented in Table 1 and [36]. For both the PTM and CFDM models grid convergence study was performed for selecting the final grid mesh; more details can be found in [16]. The FEM model is a multi-body model consisting of two rigid bodies, cylinder and mass, which are connected with vertical tendons. The second rigid body is a mass and is modeled as a point mass. Hydrodynamic effects of the second rigid body are ignored, as well as the hydrodynamic interaction between the two rigid bodies. The two bodies are free to undergo all possible motions (three translational and three rotational).  With the use of the CFDM, the decay curves of the heave and surge motions of the first rigid body are calculated and presented in Figure 6. The in-house CFD tool developed by [32,33] was used for the estimation of the decay curves in surge and heave. The heave natural frequency of the WEC predicted by the CFDM is Theave = 2.28 sec, while the surge natural period is Tsurge = 3.87 sec. Based on the decay curves, the mean velocity value M k x  is 12.62 mm/sec and 3.41 mm/sec for heave and surge, respectively. For every motion and based on the peak (positive and negative) values, Φ, of the decay curves, points with co-   With the use of the VD-PQ damping coefficients for both velocity ranges, analysis is performed for regular waves with different wave periods and heights in order the RAOs of the heave and surge motions to be calculated. Three different wave heights were examined equally to 12 mm, 25 mm and 75 mm. Similarly, analysis is performed for the case where only the radiation damping calculated by potential theory is used. The numerical analysis predictions from the two different models are compared against experimental data as reported in [36]. In Figures 7 and 8, comparisons of heave and surge of the first rigid body of the examined WEC between experimental data and numerical predictions are presented. The predictions are calculated with the use of the FEM model with the developed iterative viscous damping VD-PQ and with the use of radiation damping. For both damping models, the results follow the trends of the relevant experimental data with differences. For heave motion and irrespectively of the examined wave height, the use of the VD-PQ model results in RAOs values closer to the experimental data compared to the case that the radiation damping is only used during the analysis. This is more intense for examined wave periods in the wave period range close to the natural period in heave, Theave = 2.28 sec, where the model that uses only the radiation damping completely fails to predict the response. VD-PQ results to very efficient results for predicting the response of the WEC especially close to the resonance of the body. On average, the relative difference between numerical predictions and experimental data is 8.41% and 43.25% for VD-PQ and radiation potential models, respectively. Especially for examined wave periods close to the heave resonance of WEC, the model which uses only radiation damping results to unrealistic predictions that are up to four times larger than the experimental data values. Similar observations can be found as far as the RAOs of the surge motion. The use of the VD-PQ model results in a better estimation of the responses in all the examined wave periods and especially close to the resonance of the WEC in the surge (Tsurge = 3.87 sec). The relative difference for surge motion is 6.85% and 20.72% for VD-PQ and radiation potential models, respectively. Based on the results, it is clear that a correction of the viscous damping model should take place when studying WECs since WECs operate close to theirs.

Conclusions
In the present paper, the boundary element method is coupled with results derived by a computational fluid dynamics model for the estimation and update of the viscous damping loads in every time-step of the analysis. A coupling numerical analysis method for addressing wave-structure interaction effects through a velocity-dependent viscous damping model (VD-PQ) in the time domain is developed and proposed. The coupling between the different models is realized through a Dynamic-link library. The numerical analysis method is applied for the case of a heaving type WEC, and comparisons are made against experimental data and relevant numerical analysis predictions. The basic findings of the research are highlighted below: A method for accounting velocity-dependent viscous damping loads for every timestep of the analysis of a WEC in a boundary element model is developed and applied. VD-PQ model provides very good accuracy compared to experimental data and, at the same time, low computational cost, especially for specific types of analysis (e.g., full long-term analysis).
The consideration of instantaneous velocity of the structure in estimating the viscous damping loads for every time-step of the solution of the equation of motion is important and should be accounted especially when the wave period is close to the natural period of motions of the structure.
Viscous damping loads should be appropriately addressed and accounted when dealing with the analysis and design of WECs.
The use of radiation damping, without any possible correction, for the analysis of WECs, may result to over-predictions in motions and consequently in the produced power of the WEC.
The combination and/or coupling of different numerical analysis models may result to increase of response calculation accuracy and decrease in computational cost.
The VD-PQ viscous damping model that is proposed in the present paper can be appropriately further developed (e.g., in the elastic degrees of freedom) for the case of WECs that hydroelasticity is important and dominates their response since the coupling of the structural deformations with CFD-based models is computationally expensive. Moreover, the proposed numerical analysis model will be further developed and assessed for the case of floating structures that are free to move in all six degrees of freedom and for the case of floating wind turbines. The aforementioned issues are prioritized as future work.
Funding: This research received no external funding. Data Availability Statement: Data sharing not applicable.