Next Article in Journal
A Novel Vision-Based Towing Angle Estimation for Maritime Towing Operations
Next Article in Special Issue
Unsteady RANS CFD Simulations of Sailboat’s Hull and Comparison with Full-Scale Test
Previous Article in Journal
Path Following of a Water-Jetted USV Based on Maneuverability Tests
Previous Article in Special Issue
Numerical Investigation of Unsteady Cavitation Dynamics over a NACA66 Hydrofoil near a Free Surface
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Viscous Damping Identification for a Wave Energy Converter Using CFD-URANS Simulations

by
Marco Fontana
*,†,
Pietro Casalone
*,†,
Sergej Antonello Sirigu
,
Giuseppe Giorgi
,
Giovanni Bracco
and
Giuliana Mattiazzo
Department of Mechanical and Aerospace Engineering, Polytechnic of Turin, C.so Duca degli Abruzzi, 24, 10129 Turin, Italy
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this work.
J. Mar. Sci. Eng. 2020, 8(5), 355; https://doi.org/10.3390/jmse8050355
Submission received: 19 April 2020 / Revised: 12 May 2020 / Accepted: 14 May 2020 / Published: 17 May 2020
(This article belongs to the Special Issue CFD Simulations of Marine Hydrodynamics)

Abstract

:
During the optimization phase of a wave energy converter (WEC), it is essential to be able to rely on a model that is both fast and accurate. In this regard, Computational Fluid Dynamic (CFD) with Reynolds Averaged Navier–Stokes (RANS) approach is not suitable for optimization studies, given its computational cost, while methods based on potential theory are fast but not accurate enough. A good compromise can be found in boundary element methods (BEMs), based on potential theory, with the addition of non-linearities. This paper deals with the identification of viscous parameters to account for such non-linearities, based on CFD-Unsteady RANS (URANS) analysis. The work proposes two different methodologies to identify the viscous damping along the rotational degree of freedom (DOF) of pitch and roll: The first solely involves the outcomes of the CFD simulations, computing the viscous damping coefficients through the logarithmic decrement method, the second approach solves the Cummins’ equation of motion, via a Runge-Kutta scheme, selecting the damping coefficients that minimize the difference with CFD time series. The viscous damping is mostly linear for pitch and quadratic for roll, given the shape of the WEC analysed.

1. Introduction

Developing efficient, cost competitive and survivable wave energy converters (WECs), has proven, in the last years, to be a challenging task [1,2,3]. Fluid dynamic analyses are fundamental at every design stage to evaluate the behaviour of the WEC [4], so different methods for simulating the fluid forces and the interaction between waves and semi-submerged bodies have been developed, with increasing levels of fidelity and associated computational cost. Solving the complete set of Navier–Stokes Equations (NSE) for offshore engineering applications is expensive and not advised because the computational cost is too high compared to the accuracy desired [5]. Therefore, the NSE are simplified to obtain a linear equivalent potential flow theory (PF), where the solution is obtained through linearization, thanks to assumptions of inviscid fluid and irrotational, incompressible flow. In this way, one of the most complex set of equations is transformed into a Laplace problem. Under such hypotheses, the Boundary Element Method (BEM) can be implemented, which is a numerical approach for the solution of linear partial differential equations, expressed as integral equations. BEM represents a valid alternative to Reynolds Average Navier–Stokes (RANS) simulation: Being three (or more) orders of magnitude faster [6], BEM is used in many industrial and research applications, especially when computational resources are limited and accuracy is not of primary importance. Linear approaches produce satisfactory results when modelling linear waves and the WEC is moving around the operative point. In fact, non linearities of fluid forces are less important when oscillations are small and smooth. The success of linear models is justified by their many perks, which include, but not limited to, extreme simplicity, compliance with the superposition principle and great flexibility that makes them compatible with a vast array of mathematical tools. The linearization is always referred to the operating point. In many applications this is a reasonable assumption, since the working condition is not far respect a specific setpoint; in contrast, in the case of wave energy conversion, the objective is to drive the system as far away from equilibrium as possible, in order to maximize the produced energy. This is likely to excite nonlinear dynamics, resulting in non-representative linear models, which usually underestimate energy losses and overestimate the productivity of the WEC [7].
Two main assumptions must be satisfied when using BEMs:
  • The maximum wave steepness should be less than 1/30.
  • The oscillations around the equilibrium point should be small.
The better these hypotheses are met, the higher the accuracy of the results obtained with the linear model [8]. Although linear PF methods have been used successfully in many offshore engineering applications, the linearizing assumptions are often challenged by realistic WEC behaviour, where large amplitude motions result in viscous drag, flow separation, vortex shedding and other nonlinear hydrodynamic effects [5]. In order to improve the linear time domain model of a WEC, described by the Cummins’ equation [9], it is possible to include a viscous damping term as an external force, without significantly increment the computational time required to run the model. In this paper, a notional un-coupled linear and quadratic dependence of the viscous force on the velocity is assumed, hence considering two additive terms. It follows that two viscous drag coefficients should be identified, one for the linear contribution and one for the quadratic contribution. The addition of such non-linear damping force is necessary to avoid divergence in simulations under certain conditions: for example, the model of a forced, undamped system, is likely to explode into very large, non-realistic motion close to resonance conditions. Diverging behaviour can happen without a proper extinction coefficient and is not possible to model some of the conditions the WEC will face at sea [10].
A way to determine damping coefficients is to conduct an experimental campaign, analysing the free decay motion of the WEC [11,12]. For each degree of freedom (DOF), the WEC is released from a set of initial conditions in order to highlight when non-linearities become important; this test is usually performed to identify the roll damping for ships [13]. In order to consider only the significant kinematics of the free decay, the analysed time trace usually covers a few oscillations before the transient elapses; for different starting conditions, the resulting angular acceleration, velocity and rotation angle are post-processed to estimate damping parameters. However, experimental campaigns [14] are expensive and usually not affordable in the early design stage of a WEC. Such a critical obstacle can be overcome by using the use of a CFD Numerical Wave Tank (CNWT). CFD is a Hi-Fi approach to the problem that makes the inclusion of all non-linearities possible. However, CFD is so computationally expensive that it cannot be used for extensive analysis or study of operative conditions: Its most common use is when large non-linearities are expected (survivability, large motion) or for tuning and validating lower-fidelity models, or investigating behaviour of related devices such as tuning tanks [15]. Furthermore, turbulence is a critical feature to handle. There are mainly three techniques to deal with it: Direct Numerical Simulations (DNS), in such simulations, all of the motions contained in the flow are resolved, Large Eddy Simulations (LES) in which NSEs are filtered over the space, finally RANS in which flow quantities are split into the mean value and the perturbation. The resulting equations of LES and RANS are essentially identical, the difference is that in LES most of the turbulence energy is resolved [16]. The URANS approach is chosen because it is the best cost-efficient method considering accuracy and solver time (which could increase easily by several orders of magnitude from RANS to LES and DNS).
The aim of this paper is to define a CFD setup able to emulate a real tank experiment and to describe a methodology to identify the viscous damping coefficients along the two rotational DOFs of a WEC, namely pitch and roll. In particular, the WEC is a prismatic pitching platform that extracts energy from the pitching motion. On the one hand, viscous losses in pitch are detrimental for power extraction, since they increase the dissipation in the DOF exploited to harvest energy; on the other hand, roll response to oblique or short-crested waves is also detrimental, since it absorbs part of the incoming energy, hence decreasing the energy available to the pitch DOF. Therefore, it is important to obtain a reliable prediction of both pitch and roll responses. The goal of this paper is to increase the fidelity of models based on the potential flow theory by adding the effect of viscous non-linearities. It is worthwhile to remark that most of the literature in the wave energy field is concerned with the identification of viscous losses in the DOF used for energy extraction only (usually a smooth surfaces), while ancillary DOFs are often neglected (with sharp edges in the case study herein considered). Hence, this paper also purports to discuss specific challenges and intrinsic differences between different DOFs, despite using the same test-case and identification technique. Results carried out from simulations are consistent with expectations based on the shape of the hull: the linear term of the damping force prevails in the pitch free decay analysis while the quadratic term, that models the viscous non-linearities, is predominant in the roll DOF.
This paper is organized as follows. After a brief presentation of the device considered in Section 2, the mathematical model of the methodology presented in this paper is discussed in Section 3. The CFD simulation and its setup is explained in Section 4 with emphasis on damping techniques to avoid wave reflections and overset mesh. In the final part of CFD section the mesh convergence study is presented and discussed. Results of simulations are shown in Section 5, firstly as regards CFD scenes, then kinematics and finally the application of methodology discussed in the mathematical section. Ultimately the conclusions of Section 6 involve the analysis of the previously shown outcomes and a last post-processing phase by solving the equation of motion in the time domain, including the non-linear damping effect and the optimization of coefficient identification through the goodness of fit.

2. Device

The present case considers the WEC shown in Figure 1. Such device, so-called Inertial Sea Wave Energy Converter (ISWEC), is a pitching gyroscope-based WEC, designed for the Mediterranean Sea. ISWEC absorbs wave energy through a floating hull, that drives a gyroscope within contained by means of reacting inertial effects. In particular, waves force the hull into pitching motion, which is transmitted to the gyroscopic system located on a platform inside the floater; the reaction is due to the gyroscopic effects on the spinning flywheel, which induces a torque that is transformed by the power take off (PTO) into electrical energy.
The hull considered in this work has length L H equals to 10   m and width W H = 5 m .

3. Mathematical Model and Identification Method

The time domain linear model of a floating WEC is based on the Cummins’ equation [9].
( M + A ) X ¨ ( t ) + 0 t h r ( t τ ) X ˙ ( t ) d τ + K X ( t ) = F w ( t )
where the terms of Equation (1) are:
  • A is the matrix of added mass at infinite frequency
  • M is the mass matrix
  • F w is the external force
  • K is the stiffness matrix
  • h r is the radiation impulse response function
The convolution term represents, along with A , the effect of wave radiation in an ideal fluid and it is often referred to as fluid memory effect, because it considers the energy of the radiated waves due to the past motion of the body.
The hydrodynamic parameters in Equation (1) are computed via BEM-based software (e.g., NEMOH [17] and ANSYS Aqwa [18]). It is worth noting that in Equation (1) there are no linearity limitations on the external forces, hence the model is potentially capable to deal with non-linear contributions. Viscous effects are modelled via the addition of a nonlinear term to the right-hand side of Equation (1), divided into linear and quadratic dependence on the pitching velocity δ ˙ .
The identification of such damping parameter is conducted using the logarithmic-decrement decay method [11,12,19], which considers the rate of decay of oscillation starting from a non-equilibrium position. During a free decay test, the system oscillates, for each DOF, at its single natural frequency, thus it is possible to determine both the natural period and the hydrodynamic damping. The non-linear equation for free decay pitch motion is:
( I 55 + A ( ω δ ) ) δ ¨ ( t ) + F vis , 55 NL + K 55 δ ( t ) = 0
where the subscript 55 refers to the pitching DOF, δ is the pitch angle and the frequency-dependent added mass A is calculated at the natural frequency of the system ( ω δ ). F vis , 55 NL is the nonlinear viscous force, defined as
F vis , 55 NL = B 55 , 1 δ ˙ + B 55 , 2 δ ˙ | δ ˙ |
where B 55 , 1 is the coefficient for linear dissipation, and B 55 , 2 is the coefficient for quadratic dissipation. Note that B 55 , 1 accounts for both linear viscous damping and radiation damping. Equation (3) can be directly included in a non-linear time domain model, with the coefficients B 55 , 1 and B 55 , 2 identified as the ones that minimize the error between a fidelity-benchmark (CFD, for instance) and the time-domain model. However, it is often crucial to preserve the linear structure of the lower-fidelity model, for example to retain the ability to simulate the model in the frequency domain. Therefore, it is useful to linearize Equation (3), while preserving the ability to discriminate between linear and quadratic contributions. According to [20], it is possible to assume that for each half-cycle, the oscillation is reasonably sinusoidal. Under such an assumption, the non-linear (quadratic) term of Equation (3) can be linearized using the Fourier series expansion:
δ ˙ | δ ˙ | = 8 3 π ω δ δ i δ ˙
where δ i and ω δ are the amplitude and the frequency of oscillation of the i-th oscillation cycle, respectively. Consequently, the linearized damping force can be defined as:
F vis , 55 L = B 55 , tot δ ˙ = B 55 , 1 + 8 3 π ω δ δ i B 55 , 2 δ ˙
With the inclusion of F vis , 55 L and dividing Equation (2) by the equivalent inertia ( I 55 + A ( ω δ ) ) , it is possible to obtain the non-dimensional linear equation of motion in the canonical form:
δ ¨ + 2 α e q δ ˙ + ω δ 2 δ = 0
where
α e q = α + 4 3 π ω δ δ i β
with α = B 55 , 1 2 ( I 55 + A ( ω δ ) ) and β = B 55 , 2 ( I 55 + A ( ω δ ) ) . Equation (6) describes a linear underdamped system, for which the envelope curve of a decay starting from δ 0 is defined as:
δ = δ 0 e α e q t
Applying Equation (8) for two consecutive peaks i and i + 1 of the decay curve, it is possible to calculate the logarithmic decay:
δ i δ i + 1 = e α e q ( t i + 1 t i )
Thus, identifying the equivalent linear extinction coefficient:
α e q = 1 t i + 1 t i ln | δ i | | δ i + 1 | α + 4 3 π ω δ δ m e a n , i β
where
δ m e a n , i = | δ i | + | δ i + 1 | 2
and the damped natural pitch frequency can be calculated as:
ω δ 0 = ω δ 2 + α e q 2
This method allows to calculate the extinction curve of the α e q as a function of the mean amplitude δ m e a n for each cycle. To exploit all the available experimental data and obtain a more accurate estimation, it is possible to calculate α e q and δ m e a n for each oscillation cycle, considering both maxima and minima peaks and grouping the informations of all the decay tests for the specific DOF. The linear regression fit is then performed on the calculated points, with respect to the curve expression:
α e q = a δ m e a n + b
which allows to identify the linear and quadratic extinction coefficients:
α b
β 3 π 4 ω δ a

4. CFD Model

The analysis has been carried out on the commercial software Star-CCM+ from Siemens [21] and a Unsteady Reynold Average Navier–Stokes (URANS) approach has been used. Particular attention should be paid to the setup of the CFD simulation environment, since different fluid interfaces, fluid-structure interaction and large solid body motion present numerical problems of challenging resolution. Under such conditions, reliable results can only be achieved by an appropriate and comprehensive model definition. The following subsections detail the most relevant aspects to take under consideration when performing pitch and roll free decay tests of a floating pitching platform as the one considered in this paper.

4.1. Domain, Boundary Conditions and Damping Length

Due to its oscillation, the floater generates radiated waves that propagate towards the boundaries of the domain, effectively dissipating energy. In order to correctly represent the amount of energy that is subtracted from the floater due to wave generation, it is crucial to capture the wave pattern generated by the hull during the oscillations. In this regard, the domain needs to be large enough to contain the most energetic part of the wave pattern, freely propagating in each direction. Considering L H and W H , respectively length and width of the hull, domain overall length is 20 L H and its width 10 W H . Since the wave are generated by the hull of length L H , depth of domain is set to 2.5 L H in order to respect the deep water approximation D e p t h / W a v e   L e n g t h > 0.5 .
The choice of the dimensions of the domain also depends on wave reflection phenomena, which are a critical problem in both real and numerical tanks. In fact, reflected waves modify the wave field and ultimately affect the dynamics of the floater; Furthermore, it is difficult, if not impossible, to accurately separate the contribution of radiated and reflected waves, leading to improper measures, in case of real experiments, and to numerical noise and simulation crash, in case of numerical tanks. A possible solution consists in adding a damping zone at the boundaries of the domain in order to decrement the energy and the height of waves before they reach the frontier of the domain. Therefore, damping zones hinder the generation of reflected waves, leaving fluid domain perturbed mainly by the radiated waves for the whole duration of the significant decay.
The damping zone should have an appropriate length (at least twice the wavelength [21]) and intensity of the damping [22]. Damping waves is possible by introducing an additional resistance term S z as a source into the momentum equation relative to the vertical motion w, according to [22]:
S z = ρ f 1 + f 2 | w | e k 1 e 1 w
k = x x s d x e d x s d n
where:
  • x s d is the starting point for wave damping (along x-axis direction)
  • x e d is the end point (i.e., the boundary)
  • ρ is the water density
  • w is the vertical velocity component
  • f 1 , f 2 and n are parameters of the damping model of Choi et al. [22].
In this particular case, considering a wavelength ( λ ) comparable to the hull’s length, considering the results of Peric et al. [23] and after a trial-and-error evaluation stage, based on estimation of the simulation’s residuals, the following values have been imposed:
f 1 = 0 ; f 2 = 10 2 π λ ; n = 2
Bearing this in mind, the waves reflection phenomenon has been addressed imposing a damping factor k in Equation (17) within a distance of 20 m from the boundaries as shown in Figure 2, that represents a top-view of the sea surface of the simulation.
The boundaries conditions are shown in Figure 3. With respect to the boundaries named top and bottom, they are set as velocity inlet, because a unique phase, air and water, respectively, with zero-velocity is imposed. At side, right and left, the boundary conditions are set as pressure outlet, and the pressure is forced to be the atmospheric one for the air section and the hydrostatic one for the water part.
Being the domain a function of wavelength, a parametric approach has been used: For the different free decay tests and different initial conditions the wave produced are quite various and so the domain and mesh refinements, to maximize efficiency of the simulation.

4.2. Mesh Generation

In CFD there are two ways to deal with moving bodies: Mesh morphing/remeshing techniques are best when the movement of the body is very small, while they are too computationally expensive and potentially unstable with large displacements of the floater, since the mesh morphing task becomes a significant portion of the overall computational burden. Alternatively, when large motions are expected, overset methods should be preferred [24,25]. Thus the overset approach has been chosen for this study and the domain is split in two different overlapped parts: background and overset. The overset is the one moving with the body, while the background is fixed. With the overset approach, no remesh operation is needed because the mesh never deforms nor its quality decreases, since it moves with the body and remains unaltered. The overset and background separate regions exchange information due to an overlapping area consisting of some cells, called donor and acceptor, where the solution is interpolated trough Chimera Interpolation techniques according to [26,27]. Equations are solved separately in these different areas and bound together only at the interface between the two zones.
While the dimensions of the background are function of the wavelength, the dimensions of the overset region are constant for all the tests because they are referred to the hull and not to the wave pattern. Another important volumetric refinement is the one around the overset which ensures that, in the zone where information are exchanged, cells belonging to overset and background have the same volume, in order to reduce discretization and interpolation errors. When dealing with trimmed block mesh, it is is important to bear in mind that cells can increase their size only by double.

4.2.1. Tank

A trimmed mesh has been chosen due to the possibility to create anisotropic cells which provide a predominantly hexahedral mesh with minimal cell skewness, it results as the best discretization at the interface between the two fluids. It is worthwhile to remark that wave reflection can occur also when the mesh size suddenly changes. To avoid this problem, usually three different volumetric sea surface refinements are used in order to gradually coarsen the mesh, as shown from a prospective-side view in Figure 4 and side-view in Figure 5.
The wave damping should start before the coarsest grid level is reached since such a change in cell size is a potential source of reflection that should be included in the damping zone [28].

4.2.2. Overset Region

Two main factors influence the overset’s dimensions:
  • Overset does not have to be too large, otherwise the extremities will move too much when it rotates, generating interpolation errors and requiring a lower time step
  • Neither too small, because cells need space to grow from the wall to the interfaces; Moreover this guarantees that the strongest gradients are solved within the overset and not in the zone where information is exchanged.
For the overset sea refinement, we used the shape displayed in Figure 6 and isotropic cells because, otherwise, when the region rotates, the anisotropic sea refinement at the water line will rotate along with the hull.

4.3. Volume of Fluid

The Volume of Fluid (VOF) Multiphase model is an Eulerian numerical model suited to simulate flows of several immiscible fluids on numerical grids, capable of resolving the interface between the different phases [29]. The spatial distribution of each phase at a given time is defined in terms of a variable that is called the volume fraction. A method of calculating such distributions is to solve a transport equation for the phase volume fraction. The method uses the STAR-CCM+ Segregated Flow model which solves each of the momentum equations in turn. The VOF Multiphase model is the standard model to deal with marine environment and wave generation in a RANS simulation [21].

4.4. Turbulence, Wall Treatment, Y+

K-Epsilon is a two-equations model that solves transport equations for the turbulent kinetic energy and the turbulent dissipation rate in order to determine the turbulent eddy viscosity. Such model provides a good compromise between robustness, computational cost and accuracy and are well suited for external flows and when gradients at wall are not too strong [30,31]. In our case a Realizable Two layer K-Epsilon Model is used [32]. This model contains a new transport equation for the turbulent dissipation rate with respect to standard K-Epsilon. In addition, the turbulent viscosity is expressed as a function of mean flow and turbulence properties, rather than assumed to be constant as in the standard model. A two-layer wall treatment is adopted: this method solves for but prescribes and turbulent viscosity algebraically with distance from the wall in the viscosity-dominated near-wall flow regions. In this approach, the boundary layer is divided into two layers. The values specified in the near-wall layer are blended smoothly with the values computed from solving the transport equation far from the wall. The equation for the turbulent kinetic energy is solved across the entire flow domain. The two-layer approach, which resolves correctly for both low and high y + , using respectively standard wall function and blended wall function according to Reichardt’s law [33], fits perfectly with the free decay problem: after few oscillations the velocities are smaller and y + value decreases during the numerical experiment. Based on such considerations, a target value of y + = 1 has been chosen, since the aim of this work is to define damping coefficient accounting for non-linear effects, thus it is important to describe properly the boundary layer zone where viscous stress, separation or vortex shedding that occur. The value of y + = 1 is obtained considering the flat-plate boundary layer theory [34], obtaining a first cell height of 4.222 · 10 4 m , considering a reasonable total thickness of 3 c m and using a number of 10 layers, which guarantee a smooth grow rate of the cells.

4.5. Time Step

In Star-CCM+ the multiphase VOF model requires the use of an implicit unsteady approach. A second order time-marching discretization has been chosen since the first order is numerically dissipative and the property of the wave may not be transported correctly. First order time discretization is the only model unconditionally stable, but it has a drawback: the leading term of the truncation error of convective flux resembles a diffusive flux. This numerical, or artificial, diffusion that is introduced in the solution is magnified in multidimensional problems if the flow is oblique to the grid. A very fine grid would be necessary to obtain accurate solutions. On the other hand, since a second order approach does not require such a refined mesh-discretization, it is faster and more accurate than a first-order time-discretization. Time step is limited by two different conditions:
  • The courant number at the interface between the two fluids has to be at least less than 1 (required by the VOF model), although it is suggested to keep it less than 0.5 in order to have a clean separation of the two fluids.
  • The movement of the overset cells at the interface between the two regions must be less than half the minimum cell size, when using a second order temporal scheme, in order to prevent overset mesh problems due to the exchange of information.
The more severe limitation on the time step comes from the latter condition; therefore, it is convenient to change the time step according to the different initial conditions: since the oscillation period is always the same, the further away from equilibrium is the initial displacement, the smaller the time step needs to be in order to maintain the same courant number.

4.6. Mesh Convergence Study

In order to support the goodness of the numerical scheme created, a grid independence analysis has been carried out. This process proves the stability of the numerical method and ensures that solution does not change when refining the mesh. This convergence study is based on the number of cells used to discretize the wave height that is generated by the oscillations of the hull. In the wave pattern there is the energy lost by the device during the free decay and it is fundamental to capture it accurately. Moreover, prism layer properties are base size independent and computed as absolute values, thus base size has little influence in the overall mesh. The number of layers that are necessary to discretize the wave is very dependent on wave steepness; according to [5], at least 10 layers are suggested and even more for a sharp interface between fluids. Not using enough cells some information may be lost and wave transport and generation may be inaccurate, resulting in overly-damped or magnified waves. Cells count is very dependent on the number of layers of sea refinement because most of the cells come from this mesh control and from the overset zone, as confirmed by results presented in Figure 7.
In order to choose the appropriate mesh sizes CPU time is considered, as shown in Figure 8: passing from 15 cells to 20 means an increase of about 115%. So the configuration of 20 cells per wave height has been excluded because of the prohibitive computational time.
The cells count of the overset is related with the size of the sea refinement, because cells at the interface between background and overset must have same volume, especially at sea surface. Convergence study is done on the kinematics of the free decay with respect to Figure 9, which is also the data used for the post-processing and evaluation of damping coefficients.
The wave height with several probes has also been monitored in order to understand the impact of mesh resolution in areas far from the floater.
The coarsest mesh setups have been compared to the finest ones corresponding to 20 cells in terms of relative errors and the associated root-mean-square as shown in Figure 10 and Figure 11.
Finally, 15 cells per wave height have been chosen, in order to ensure a good compromise between wave capturing and computational cost specially considering the trend of errors in Figure 10, the low related rms (Figure 11) and its computational cost in terms of time (Figure 8).

5. Results

5.1. CFD Scenes

In this section, results are presented in terms of scalar and vector scenes in order to show the validity and the reliability of the conducted CFD simulations for pitch and roll free decay motions. An example of the wave pattern generated by the pitching motion is shown in Figure 12, starting from a non-equilibrium angle equal to 10 . The outcomes are consistent with the physical phenomenon.
The hydrostatic pressure referred to the waterline is well defined as shown in the example of a roll simulation in Figure 13.
The interface is smooth and well-captured, as shown in Figure 14, where the line convolution integral (LIC) of velocity is presented. Such a visualization tends to show a texture where points along the same streamline have the same colour; Displaying LIC of velocity is more computationally efficient the actual streamline scene which the latter is very expensive in terms of computational costs. Generally this technique [35] involves convoluting a white noise image along streamlines computed from the vector field. In the resulting image, the streamlines cover the entire domain of the vector field. This technique has the advantage of being able to visualize large and detailed vector fields in a reasonable display area. Compared with simpler integration-type approaches, which entail following the flow vector at each point to produce a line, LIC produces a whole image at every step.
Lastly, the target value of y + = 1 on the wet surface of the hull, on average, has been correctly obtained, as shown in Figure 15. The scene in Figure 15 refers to an intermediate time step and proves an adequate mesh sizing within the scope of capturing a low value of y + which results into an accurate interpolation by the solver in terms of wall stresses and frictional velocity.

5.2. Post-Processing and Damping Coefficients Identification

The kinematics of pitch and roll motion versus physical time are presented in Figure 16 and Figure 17, respectively. In order to graphically represent and evaluate the importance of non-linearities, each free decay time trace has been normalized with respect to the initial displacement, as reported in Figure 18 and Figure 19, respectively for pitch and roll. While the pitch normalized responses significantly overlap, suggesting linearity, roll normalized responses are diverse, showing a clear nonlinear behaviour. In particular, the larger the initial displacement, the faster the decaying rate, suggesting a higher content of non-linearity. This is hereafter quantified through the identification of the nonlinear damping coefficient β , as discussed in Section 3.
The different behaviour of pitch and roll in terms of damping is clearly presented in kinematics figures, and the roll motion results less damped than pitch. This behaviour is justified considering the radiation curves related to the hull of this work, shown in Figure 20 where the pitch motion has a more significant value with respect to roll damping, on a point of detail greater than roll of about two orders of magnitude. When oscillations become small and non-linearities less important, roll needs more time to reduce the amplitude of motion. Anyway, roll is strongly damped at the beginning, when oscillations are big because the quadratic part, i.e., proportional to the square of velocity, is predominant in roll viscous damping while the linear part is almost negligible.
Each pair of consecutive peaks (considering both maxima and absolute values of minima) defines a α e q and a δ m e a n , i , according to Equations (10) and (11), respectively. Such parameters are shown in Figure 21 and Figure 22 for pitch and roll responses, respectively.
As regards peaks data in Figure 21 and Figure 22, the wide spread in pitch data at low angles is justified by oscillations around the equilibrium point. Such low angles measurements are negligible as far as the identification of viscous damping coefficient is concerned and thus will be neglected.
About the roll motion the previous spread is absent because it is overall less damped; In addition two trends are clear (Figure 22), one until about 2 , the other one afterwards. This behaviour is possibly related to the linearization assumptions done in order to define and obtaine α e q . Therefore, the identification regards only the left side of the plot, where points are well aligned and so the linearization appears to be more valid.
Taking into account the Equation (10) and also referring to the plots in Figure 23 and Figure 24, the more non-linear behaviour of roll is confirmed. Indeed, considering Equation (10), if α represents the intercept and β the slope, it is clear the different pattern both in pitch and roll motion. Pitch kinematics are more damped than the roll one, but the first has a linear behaviour instead of the second one which results less damped, but mainly quadratic. With regards to roll data in Figure 24, points corresponding to δ 0 = 10 differ from the 5 and 7 because the first case involves much larger angles and consequently a more non-linear behaviour.
Numerical results are summarized in Table 1 in terms of coefficient of determination R 2 of the linear fit, the period of oscillations T and finally α and β respectively corresponding to intercept and slope of the linear regression, or in a more physical meaning, the linear and quadratic damping terms.
The computation of values of α and β throughout the logarithmic decrement here proposed can be considered as a starting step for solving Equation (6). Taking into account the previous outcomes, the equation is solved with a Runge-Kutte scheme (RK) of 4th-order accuracy with an adaptative time step [36].
Furthermore, this technique has been enriched considering a parametrization of 100 values α and β centred on the outcomes computed through the linear regression and shown in Table 1. For each value, the goodness of fit has been computed with the cost function Normalized Mean Square Error:
f i t = 1 | | x r e f x | | 2 | | x r e f m e a n ( x r e f ) | | 2
where x r e f and x represents respectively the reference time history (i.e., the one obtained from CFD) and the one computed through RK-solver, the goodness of fit equals to 1 means a perfect match of results. By way of example and for sake of brevity, outcomes are presented in a visual way through plot in Figure 25 and Figure 26, considering pitch and roll free decay with an initial condition of 10 , showing CFD results, Equation (6) solved with values in Table 1 and finally with α and β corresponding to the maximum of goodness of fit (GoF).
For clarity’s sake, zoomed portions of plot are presented in Figure 27 and Figure 28 as far as the first oscillations where large amplitude occur and Figure 29 and Figure 30 as regards the final tail.
A different behaviour of the optimization of α and β is observed as regards pitch and roll. Concerning the first, the maximization of GoF tends to a best fit of larger oscillations as shown in Figure 27 and a poorest fit of the final tail of the decay (Figure 29); The opposite behaviour is evident concerning the roll motion, where the largest amplitudes are not perfectly fitted by the optimal values of α and β (Figure 28), even if the fit is nevertheless better than the guess value of Table 1. The explanation of such trend is analogous to the considerations about the cloud of data in Figure 21 and Figure 22. Results are summarized in Table 2 and Table 3, respectively for pitch and roll motion.
where
  • δ 0 is the initial condition for the decay in degrees
  • α and β are the linear and quadratic terms of the equivalent extinction coefficient α e q in Equation (10)
  • GoF is the the goodness of fit defined in Equation (19)
  • R 2 is the ordinary coefficient of determination of the linear fit model

6. Conclusions

In conclusion, in this work we propose an affordable methodology with the aim of the identification of non-linear viscous damping parameters. Results obtained are in line with the expectations and consistent with the physics and experience about the WEC device considered. The hi-fi CFD techniques could be performed one-time considering few cases (in this work three different initial conditions for the free decay). This first investment in terms of computational costs and effort can lead to the evaluation of those non-linear coefficients, related to significant physical phenomena. This identification provides an important enrichment of low-order models used for next design and analysis steps.

Author Contributions

Conceptualization, M.F. and P.C.; data curation, M.F., P.C., S.A.S. and G.G.; formal analysis, M.F. and P.C.; funding acquisition, G.B. and G.M.; investigation, M.F. and P.C.; methodology, M.F., P.C., S.A.S. and G.G.; project administration, S.A.S., G.G., G.B. and G.M.; resources, M.F., P.C., G.B. and G.M.; software, M.F. and P.C.; supervision, S.A.S., G.G., G.B. and G.M.; validation, M.F., P.C., S.A.S. and G.G.; visualization, M.F. and P.C.; writing—original draft, M.F. and P.C.; writing—review and editing, M.F., P.C., S.A.S. and G.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

Computational resources provided by hpc@polito (http://hpc.polito.it).

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Mattiazzo, G. State of the Art and Perspectives of Wave Energy in the Mediterranean Sea: Backstage of ISWEC. Front. Energy Res. 2019, 7, 114. [Google Scholar] [CrossRef] [Green Version]
  2. Pozzi, N.; Bonetto, A.; Bonfanti, M.; Bracco, G.; Dafnakis, P.; Giorcelli, E.; Passione, B.; Sirigu, S.; Mattiazzo, G. PeWEC: Preliminary design of a full-scale plant for the mediterranean sea. In Proceedings of the NAV International Conference on Ship and Shipping Research, Trieste, Italy, 20–22 June 2018; pp. 504–514. [Google Scholar] [CrossRef]
  3. van Rij, J.; Yu, Y.H.; Edwards, K.; Mekhiche, M. Ocean power technology design optimization. Int. J. Mar. Energy 2017, 20. [Google Scholar] [CrossRef]
  4. Sirigu, A.S.; Gallizio, F.; Giorgi, G.; Bonfanti, M.; Bracco, G.; Mattiazzo, G. Numerical and Experimental Identification of the Aerodynamic Power Losses of the ISWEC. J. Mar. Sci. Eng. 2020, 8, 49. [Google Scholar] [CrossRef] [Green Version]
  5. Windt, C.; Davidson, J.; Ringwood, J.V. High-fidelity numerical modelling of ocean wave energy systems: A review of computational fluid dynamics-based numerical wave tanks. Renew. Sustain. Energy Rev. 2018, 93, 610–630. [Google Scholar] [CrossRef] [Green Version]
  6. Giorgi, G.; Ringwood, J.V. Analytical Formulation of Nonlinear Froude-Krylov Forces for Surging-Heaving-Pitching Point Absorbers. In Proceedings of the ASME 37th International Conference on Ocean, Offshore and Arctic Engineering, Madrid, Spain, 17–22 June 2018. [Google Scholar]
  7. Davidson, J.; Giorgi, S.; Ringwood, J.V. Linear parametric hydrodynamic models for ocean wave energy converters identified from numerical wave tank experiments. Ocean. Eng. 2015, 103, 31–39. [Google Scholar] [CrossRef] [Green Version]
  8. Penalba, M.; Giorgi, G.; Ringwood, J.V. Mathematical modelling of wave energy converters: A review of nonlinear approaches. Renew. Sustain. Energy Rev. 2017, 78, 1188–1207. [Google Scholar] [CrossRef] [Green Version]
  9. Cummins, W. The Impulse Response Function and Ship Motions; Technical Report; David Taylor Model Basin: Washington, DC, USA, 1962. [Google Scholar]
  10. Passione, B. Hydrodynamic Analysis and Mooring Design of a Floating Pitching Wave Energy Converter. Ph.D. Thesis, Politecnico di Torino, Torino, Italy, 2018. [Google Scholar]
  11. Begovic, E.; Bertorello, C.; Prpic Orsic, J. Roll damping coefficients assessment and comparison for round bilge and hard chine hullforms. In Proceedings of the ASME 32nd International Conference on Ocean, Offshore and Arctic Engineering, Nantes, France, 9–14 June 2013. [Google Scholar]
  12. Begovic, E.; Mancini, S.; Day, A.; Incecik, A. Applicability of CFD methods for roll damping determination of intact and damaged ship. In High Performance Scientific Computing Using Distributed Infrastructures: Results and Scientific Applications Derived from the Italian PON ReCaS Project; World Scientific: Singapore, 2017; pp. 343–359. [Google Scholar]
  13. Begovic, E.; Day, A.H.; Incecik, A.; Mancini, S.; Pizzirusso, D. Roll damping assessment of intact and damaged ship by CFD and EFD methods. In Proceedings of the 12th International Conference on the Stability of Ships and Ocean Vehicles, Glasgow, UK, 14–19 June 2015; pp. 13–19. [Google Scholar]
  14. Sirigu, S.A.; Bonfanti, M.; Begovic, E.; Bertorello, C.; Dafnakis, P.; Bracco, G.; Mattiazzo, G. Experimental Investigation of Mooring System on a Wave Energy Converter in Operating and Extreme Wave Conditions. J. Mar. Sci. Eng. 2020, 8, 180. [Google Scholar] [CrossRef] [Green Version]
  15. Sirigu, S.; Bonfanti, M.; Dafnakis, P.; Bracco, G.; Mattiazzo, G.; Brizzolara, S. Pitch Resonance Tuning Tanks: A novel technology for more efficient wave energy harvesting. In Proceedings of the OCEANS 2018 MTS/IEEE Charleston, Charleston, SC, USA, 22–25 October 2018. [Google Scholar] [CrossRef]
  16. Ferziger, J.H.; Perić, M. Computational Methods for Fluid Dynamics; Springer: Berlin/Heidelberg, Germany, 2002. [Google Scholar] [CrossRef]
  17. Penalba, M.; Kelly, T.; Ringwood, J. Using NEMOH for modelling wave energy converters: A comparative study with WAMIT. In Proceedings of the 12th European Wave and Tidal Energy Conference, Cork, Ireland, 27 August–2 September 2017. [Google Scholar]
  18. Ansys, A. Theory Manual; ANSYS® Academic Research; ANSYS: Canonsburg, PA, USA, 2014. [Google Scholar]
  19. Zhao, W.; Efthymiou, M.; McPhail, F.; Wille, S. Nonlinear roll damping of a barge with and without liquid cargo in spherical tanks. J. Ocean. Eng. Sci. 2016, 1, 84–91. [Google Scholar] [CrossRef] [Green Version]
  20. Chakrabarti, S.K. Offshore Structure Modeling; World Scientific: Singapore, 1994; Volume 91994. [Google Scholar]
  21. CD-adapco, S. STAR CCM+ User Guide Version 12.04; CD-Adapco: New York, NY, USA, 2017. [Google Scholar]
  22. Choi, J.; Yoon, S.B. Numerical simulations using momentum source wave-maker applied to RANS equation model. Coast. Eng. 2009, 56, 1043–1060. [Google Scholar] [CrossRef]
  23. Perić, R.; Abdel-Maksoud, M. Reliable damping of free-surface waves in numerical simulations. Ship Technol. Res. 2016, 63, 1–13. [Google Scholar] [CrossRef] [Green Version]
  24. Windt, C.; Davidson, J.; Chandar, D.; Ringwood, J. On the importance of advanced mesh motion methods for WEC experiments in CFD-based numerical wave tanks. In Proceedings of the VIII International Conference on Computational Methods in Marine Engineering, Gothenburg, Sweden, 13–15 May 2019. [Google Scholar]
  25. Davidson, J.; Karimov, M.; Szelechman, A.; Windt, C.; Ringwood, J. Dynamic mesh motion in openfoam for wave energy converter simulation. In Proceedings of the 14th OpenFOAM Workshop, Duisburg, Germany, 23–26 July 2019. [Google Scholar]
  26. Meakin, R.L. Composite overset structured grids. In Handbook of Grid Generation; CRC Press: Boca Raton, FL, USA, 1999; pp. 1–20. [Google Scholar]
  27. Petersson, N.A. Hole-cutting for three-dimensional overlapping grids. SIAM J. Sci. Comput. 1999, 21, 646–665. [Google Scholar] [CrossRef]
  28. Field, A. How the Wave Damping Length Influences the Domain and Wake Refinements in a Hull Performance Analysis. 2017. Available online: https://thesteveportal.plm.automation.siemens.com/articles/en_US/FAQ/How-the-Wave-Damping-Length-can-influence-the-mesh-domain-and-the-wake-refinement-sizes-in-a-Hull-performances-analysis (accessed on 14 May 2020).
  29. Hirt, C.W.; Nichols, B.D. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 1981, 39, 201–225. [Google Scholar] [CrossRef]
  30. Menter, F. Zonal two equation kw turbulence models for aerodynamic flows. In Proceedings of the 23rd Fluid Dynamics, Plasmadynamics, and Lasers Conference, Orlando, FL, USA, 6–9 July 1993; p. 2906. [Google Scholar]
  31. Menter, F.R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 1994, 32, 1598–1605. [Google Scholar] [CrossRef] [Green Version]
  32. Shih, T.H.; Liou, W.W.; Shabbir, A.; Yang, Z.; Zhu, J. A new k-epsilon eddy viscosity model for high Reynolds number turbulent flows: Model development and validation. In NASA Technical Reports Server; Document ID: 19950005029; NASA: Washington, DC, USA, 1994. [Google Scholar]
  33. Reichardt, H. Vollständige Darstellung der turbulenten Geschwindigkeitsverteilung in glatten Leitungen. ZAMM J. Appl. Math. Mech. 1951, 31, 208–219. [Google Scholar] [CrossRef]
  34. White, F. Fluid Mechanics; McGraw-Hill International Editions; McGraw-Hill: New York, NY, USA, 2003. [Google Scholar]
  35. Cabral, B.; Leedom, L.C. Imaging vector fields using line integral convolution. In Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques, Anaheim, CA, USA, 2–6 August 1993; pp. 263–270. [Google Scholar]
  36. Dormand, J.R.; Prince, P.J. A family of embedded Runge-Kutta formulae. J. Comput. Appl. Math. 1980, 6, 19–26. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Schematic representation of the Inertial Sea Wave Energy Converter (ISWEC) device hull, gyroscope, power take off (PTO) and reference system with respect to the incoming wave.
Figure 1. Schematic representation of the Inertial Sea Wave Energy Converter (ISWEC) device hull, gyroscope, power take off (PTO) and reference system with respect to the incoming wave.
Jmse 08 00355 g001
Figure 2. VOF Damping Factor k.
Figure 2. VOF Damping Factor k.
Jmse 08 00355 g002
Figure 3. Domain boundary conditions.
Figure 3. Domain boundary conditions.
Jmse 08 00355 g003
Figure 4. Sea surface refinements in prospective view.
Figure 4. Sea surface refinements in prospective view.
Jmse 08 00355 g004
Figure 5. Sea surface refinements in side-view.
Figure 5. Sea surface refinements in side-view.
Jmse 08 00355 g005
Figure 6. Hull mesh refinement based on maximum rotation angle.
Figure 6. Hull mesh refinement based on maximum rotation angle.
Jmse 08 00355 g006
Figure 7. Number of total cell varying the cells per wave height.
Figure 7. Number of total cell varying the cells per wave height.
Jmse 08 00355 g007
Figure 8. CPU Time (s) for different cells per wave height.
Figure 8. CPU Time (s) for different cells per wave height.
Jmse 08 00355 g008
Figure 9. Different kinematics in terms of pitch angle for various cells per wave height.
Figure 9. Different kinematics in terms of pitch angle for various cells per wave height.
Jmse 08 00355 g009
Figure 10. Relative error for different meshes.
Figure 10. Relative error for different meshes.
Jmse 08 00355 g010
Figure 11. RMS of error compared to 20 cells per wave height.
Figure 11. RMS of error compared to 20 cells per wave height.
Jmse 08 00355 g011
Figure 12. Wave pattern for a pitching free decay simulation, starting from a 10 initial displacement.
Figure 12. Wave pattern for a pitching free decay simulation, starting from a 10 initial displacement.
Jmse 08 00355 g012
Figure 13. Hydrostatic pressure referred to the waterline.
Figure 13. Hydrostatic pressure referred to the waterline.
Jmse 08 00355 g013
Figure 14. Line Integral Convolution of velocity in roll simulation.
Figure 14. Line Integral Convolution of velocity in roll simulation.
Jmse 08 00355 g014
Figure 15. Wall Y+ on wet surface, grey color refers to the surface in air.
Figure 15. Wall Y+ on wet surface, grey color refers to the surface in air.
Jmse 08 00355 g015
Figure 16. Pitch kinematic.
Figure 16. Pitch kinematic.
Jmse 08 00355 g016
Figure 17. Roll kinematic.
Figure 17. Roll kinematic.
Jmse 08 00355 g017
Figure 18. Normalized Pitch kinematic.
Figure 18. Normalized Pitch kinematic.
Jmse 08 00355 g018
Figure 19. Normalized Roll kinematic.
Figure 19. Normalized Roll kinematic.
Jmse 08 00355 g019
Figure 20. Radiation curves.
Figure 20. Radiation curves.
Jmse 08 00355 g020
Figure 21. All pitch peaks data.
Figure 21. All pitch peaks data.
Jmse 08 00355 g021
Figure 22. All roll peaks data.
Figure 22. All roll peaks data.
Jmse 08 00355 g022
Figure 23. Cut pitch data and outcomes.
Figure 23. Cut pitch data and outcomes.
Jmse 08 00355 g023
Figure 24. Cut roll data and outcomes.
Figure 24. Cut roll data and outcomes.
Jmse 08 00355 g024
Figure 25. Pitch results.
Figure 25. Pitch results.
Jmse 08 00355 g025
Figure 26. Roll results.
Figure 26. Roll results.
Jmse 08 00355 g026
Figure 27. Large oscillations in pitch.
Figure 27. Large oscillations in pitch.
Jmse 08 00355 g027
Figure 28. Large oscillations in roll.
Figure 28. Large oscillations in roll.
Jmse 08 00355 g028
Figure 29. Tail pitch results.
Figure 29. Tail pitch results.
Jmse 08 00355 g029
Figure 30. Tail roll results.
Figure 30. Tail roll results.
Jmse 08 00355 g030
Table 1. Numerical results for pitch and roll motion post-processing.
Table 1. Numerical results for pitch and roll motion post-processing.
R 2 T (s) α β
Pitch0.924.050.0890.136
Roll0.824.630.0051.38
Table 2. Outcomes for pitch free decay motion.
Table 2. Outcomes for pitch free decay motion.
Runge-KuttaLinear Regression
δ 0 (deg) α β GoF α β R 2
100.08180.03410.98690.08970.10260.8877
70.08050.03870.96400.08980.12810.8652
50.07910.03840.96150.91530.15380.4568 *
All dataN.A. N.A. N.A. 0.08940.11020.9436
* Low value due to the small amout of points after the angle cut. The solution of differential equation can not be performed considering all the initial conditions.
Table 3. Outcomes for roll free decay motion.
Table 3. Outcomes for roll free decay motion.
Runge-KuttaLinear Regression
δ 0 (deg) α β GoF α β R 2
100.00881.4950.98730.01211.03560.9968
70.00871.53680.99190.01241.39620.9968
50.00511.55770.99610.00141.5680.9810
All dataN.A. 2 N.A. 2 N.A. 2 0.00491.38260.8188

Share and Cite

MDPI and ACS Style

Fontana, M.; Casalone, P.; Sirigu, S.A.; Giorgi, G.; Bracco, G.; Mattiazzo, G. Viscous Damping Identification for a Wave Energy Converter Using CFD-URANS Simulations. J. Mar. Sci. Eng. 2020, 8, 355. https://doi.org/10.3390/jmse8050355

AMA Style

Fontana M, Casalone P, Sirigu SA, Giorgi G, Bracco G, Mattiazzo G. Viscous Damping Identification for a Wave Energy Converter Using CFD-URANS Simulations. Journal of Marine Science and Engineering. 2020; 8(5):355. https://doi.org/10.3390/jmse8050355

Chicago/Turabian Style

Fontana, Marco, Pietro Casalone, Sergej Antonello Sirigu, Giuseppe Giorgi, Giovanni Bracco, and Giuliana Mattiazzo. 2020. "Viscous Damping Identification for a Wave Energy Converter Using CFD-URANS Simulations" Journal of Marine Science and Engineering 8, no. 5: 355. https://doi.org/10.3390/jmse8050355

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop