On the use of OpenFOAM to model Oscillating wave surge converters

Abstract The computational fluid dynamic (CFD) toolbox OpenFOAM is used to assess the applicability of Reynolds-averaged Navier–Stokes (RANS) solvers to the simulation of oscillating wave surge converters (OWSC) in significant waves. Simulation of these flap type devices requires the solution of the equations of motion and the representation of the OWSC׳s motion in a moving mesh. A new way to simulate the sea floor inside a section of the moving mesh with a moving dissipation zone is presented. To assess the accuracy of the new solver, experiments are conducted in regular and irregular wave traces for a full three dimensional model. Results for acceleration and flow features are presented for numerical and experimental data. It is found that the new numerical model reproduces experimental results within the bounds of experimental accuracy.


Introduction
The oscillating wave surge converter (OWSC) consists of a bottom hinged floating flap as shown in Fig. 1. The waves act on the flap and force it to move back and forth. This motion can be used to generate electricity, for example using an hydraulic power take off system. This type of system is typically installed in shallow water where the horizontal fluid motion is larger than in the deep sea. Further details of this design have been detailed in Folley et al. (2007) and Renzi et al. (2014).
While using numerical simulations of ships in a seaway is by now common engineering practice, the simulation of an OWSC is not straightforward. Qian et al. (2005) presented results for the interaction of a wave driven rotating vane and a shoreline. Simulations were performed using the interface-capturing Cartesian cut cell flow solver AMA-ZON-SC, without considering viscous effects and for a two dimensional case. Schmitt et al. (2012a) compared pressure distributions derived from various numerical tools with experimental data for a fixed OWSC in waves. Results of fully viscous CFD simulations obtained with OpenFOAM showed very good agreement with experimental data. The paper also highlights the issues encountered when applying linearised potential codes like WAMIT to the case of an OWSC. Renzi and Dias (2012) developed a semi-analytical linearised potential solution method and successfully applied it to explain resonance effects encountered during experiments in small amplitude waves. Mahmood and Huynh (2011) presented two dimensional simulations of a bottom hinged vane in oscillating single phase flow. Bhinder et al. (2012) employed the Flow3d CFD code to obtain drag coefficients for an OWSC, oscillating in translational modes only. The body consisted of a cube and was not excited by waves but forced to oscillate. This work highlights the importance of viscosity for these types of devices, they estimated performance reductions of almost 60% when comparing non-viscous and viscous solutions. Rafiee et al. (2013) employed a smoothed particle hydrodynamics (SPH) method to simulate two and three dimensional cases of an OWSC. Viscosity was modelled by a k-ϵ turbulence model and results were compared to experimental data. No quantitative error estimates were given but agreement for flap rotation and pressure at various locations seems to compare well. It should be noted that the cases presented extreme events, that is overtopping waves, are investigated. The wave maker consisted of a moving piston. Results highlight the need of performing three dimensional simulations and thus the importance of the flow around the sides for the motion of the flap. Schmitt et al. (2012b) reviewed the numerical simulation demands of the wave power industry and compared the application of fully viscous CFD solvers to experimental tank tests. Simulation results were shown for cases run in OpenFOAM using a mesh distortion method to accommodate the flap motion and compared well to experimental data. The paper also gives examples of useful applications of CFD tools in the design of an OWSC, while a comparison of run times and cost estimates highlights the necessity of experimental tank testing as a tool in the wave power industry. Recently Palm et al. (2013) presented simulations of a moored wave energy converter. While the fluid forces and motions were solved in OpenFOAM, mooring loads were calculated in a coupled structural code.
Research on OWSCs has thus mainly been based on experimental, model scale tank testing. In tank tests large areas of separation and vortices of the order of magnitude of the flap width can be observed. During a wave cycle these large flow features move around the flap's side and interact with newly created vortices. While RANS CFD methods with wall functions have successfully been applied to turbulent flows, it is not clear whether the aforementioned flow effects and their effect on the flap's motion can be captured with these models. Small design changes to the flap could well affect the separation point, dissipation and other viscous effects. Before numerical tools can be used for shape optimisation or similar research, validation against experimental results is required.
Many floating bodies on a fluid surface can readily be simulated with a mesh distortion method. However, a typical OWSC rotates 7401 during normal operating conditions and up to 7 801 in extreme conditions. Mesh distortion methods usually fail due to highly distorted cells between the bottom and the flap. Remeshing is a possible but very expensive option. In this work we present an algorithm that avoids these issues. The flap moves within a cylindrical mesh zone without distorting any cells. The coupling with the surrounding static mesh is implemented using an efficient arbitrary mesh interface (AMI). The bottom of the tank is simply taken into account by setting a dissipation parameter.
Simulation results are compared to tank tests performed in Queen's University Belfast and show very good agreement.

Numerical model
The fluid solver employed in this numerical study is the interDyMFOAM solver from the OpenFOAM toolbox. The method is based on the volume of fluid algorithm for incompressible flows. A more detailed description can be found in Rusche (2002) and Berberović et al. (2009). The two main extensions to the code are libraries for the equations of motion and mesh motion algorithm. These will be presented in more detail in the following sections. The wavemaker used is based on the method presented in Choi and Yoon (2009). As a turbulence model the standard SST model was used.
This section gives an overview of the interFOAM solver class as provided by the OpenFOAM community and extensions developed for the simulation of WECs. More information on general CFD methods can be found in Versteeg and Malalasekera (2007) and Ferziger and Peric (2002). Detailed explanations of OpenFOAM are given in Weller et al. (1998)  The mass conversation and Navier-Stokes equation are given as where the viscous stress tensor is T ¼ 2μSÀ 2μð∇UÞI=3 with the mean rate of strain tensor S ¼ 0:5½∇U þð∇UÞ T and the body forces f b . In the volume of fluid method only one effective flow velocity exists. The different fluids are identified by a variable γ which is bounded between 0 and 1. A value of 0.5 would thus mean the cell is filled with equal volume parts of both fluids. Intensive properties of the flow like the density ρ are evaluated depending on the species variable γ and the value of each species ρ b and ρ f : The transport equation for γ is: The interface between the two fluids requires special treatment to maintain a sharp interface, numerical diffusion would otherwise 'mix' the two fluids over the whole domain. In OpenFOAM the interface compression treatment is derived from the two-fluid Eulerian model for two fluids denoted with the subscript l and g as given by (Berberović et al., 2009) This equation can be rearranged to an evolution equation for γ, with U r ¼ U l À U g being the relative or compression velocity: The new transport equation for γ now contains a term which is zero inside a single species but sharpens the interface between two fluids. This formulation removes the need of specialised convection schemes as used in other codes.
With n f as the face unit normal flux depending on the gradient of the species ∇γ the relative velocity at cell faces is evaluated with ϕ being the face volume flux: where δ n is a factor to account for non-uniform grids, C γ is a user defined variable to control the magnitude of the surface compression when the velocities of both phases are of the same magnitude. In the present study C γ of one was used, which yields conservative compression. The face volume fluxes are evaluated as a conservative flux from the velocity pressure coupling algorithm and not as usual from cell centre to face interpolation. A wave-maker based on the work presented by Choi and Yoon (2009) was implemented by adding a source term to the momentum equation. In the current implementation the source term is defined as the product of density ρ, the scalar field defining the wave-maker region r and the analytical solution of the wave velocity U ana at each cell centre yielding the adapted impulse equation: The beach is modelled in a similar way by implementing a dissipative source term s Á ρ Á U in the impulse equation (1). The dissipation parameter s can then be set to model the beach and has no effect where set to zero. Tests by Schmitt and Elsaesser (2015) have shown that a beach extending over approximately one wavelength and with a value s increasing from 0 to 5 effectively removes any reflections. Such a beach was used in all simulations. The parameter s is also used to take into account the sea floor in the rotating mesh, as will be explained in detail later. The computational domain consists of two mesh regions, a cylindrical moving mesh surrounding the flap and a static mesh representing the remaining tank geometry, Fig. 2. The boundary conditions used are standard conditions for fixed or moving walls for all outer walls and the flap, that is zeroGradient for pressure and zero flux conditions for velocity. Only the patch describing the top of the domain was set to a fixed pressure and velocity to pressureInletOu-tletVelocity type, which applies a zero-gradient for outflow, while for inflow the velocity is set as the normal component of the internal-cell value. The two domains are coupled via two cylindrical patches, using arbitrary mesh interface (AMI) patches.
The most important term is the convection term in the Navier Stokes equation, the linearLimited discretisation scheme was used for all simulations.

Equations of motion
Under the assumption that the fluid solver gives correct results for the hydrodynamic forces F Hydro on a body, other outer force components like gravity and damping forces can be added to obtain the total outer forces on the body F.
The instantaneous acceleration a can then easily be obtained by dividing the force F by the mass m: F and m stand for generalised forces (including moments) and masses (inertia). Integration of acceleration in time yields velocity, integration of velocity yields the bodies' new position.
In dense fluids the hydrodynamic force changes during one time-step, this effect can be considered as an added mass. Not considering this added mass leads to wrong values for the acceleration, see Bertram (2001). It is possible to use iterative methods to move the body and evaluate the forces within each time step until the value for a converges and the new equilibrium position is found. This implicit method will always yield the correct position for each time step and is unconditionally stable. It could also be expected to be less dependent on the size of the time step.
Interestingly, few people seem to be aware of the physical meaning of this effect, although they do notice that iterative procedures to fulfil Eq. (9) need under-relaxation (Hadẑić et al., 2005).
In this work, the forces on the body are averaged over several time-steps, thus avoiding inner iterations while implicitly taking into account the effect of added mass.
The algorithm used in all simulations of moving flaps in this work is explained in detail in the following section, the coordinate reference system and main variables are shown in Fig. 3.
The total hydrodynamic torque around the hinge M

!
Hydro;n is evaluated as a vector for the current time-step n by integrating pressure and viscous shear forces over the patch describing the flap surface.
The mass moment M ! mass is evaluated as follows: where CoG ! n is the position of the centre of gravity at time-step n and x Hinge is the hinge position.
The total torque for the current time step M Total;n around the hinge is then evaluated as the sum of all components around the hinge axis vector of unit length a ! : M Total;n is then saved for future time-steps and smoothed by averaging over the total moments of up to four preceding time steps: with N w being the number of weights w k larger than zero. In all simulations presented in this work three weights with a value of 1 were used.
The new rotational velocity _ ϕ n þ 1 can now be obtained as with the current time step δt and the flaps inertia around the hinge I Hinge . Similarly the change in rotation angle ϕ can be obtained as The position of the centre of gravity is then updated to the new position, employing Rodrigues' formula (Mebius, 2007): with a ! as the directional unit vector of the flaps hinge axis and x ! Hinge;CoG the vector from the position of the centre of gravity at the start of the simulation CoG ! 0 to the hinge position x ! Hinge . The reason for evaluating the new position of the centre of gravity from the initial position at the start of the simulation, and not from the preceding time-step, is to avoid accumulation of numerical errors.
The algorithm described above was implemented in Open-FOAM as a new body motion solver. The method can be called from any mesh motion solver.

Mesh motion
To adapt the changing computational domain when simulating moving bodies, different algorithms are available. Mesh distortion methods preserve the mesh topology but depending on the motion of the body can result in collapsing or distorted cells. It is also possible to re-mesh all or only parts of the domain to maintain mesh quality but re-meshing is often computationally expensive. In this work, the flap is moved with a cylindrical subset of the mesh. The interface to the static domain is modelled with a sliding interface (Farrell and Maddison, 2011). The representation of the sea floor, which is usually close to the hinge and thus inside the moving domain, is achieved by setting a dissipation parameter in all cells below a defined z-coordinate. The dissipation parameter acts as a negative source term in the impulse equation and reduces the flow velocity. With this new method the mesh quality around the flap is preserved without performing expensive re-meshing even when simulating arbitrary angles of rotation, while enabling the simulation of flaps rotating around a hinge close to the sea floor. The mesh motion method was implemented in the OpenFOAM framework. The actual mesh motion method requires specification of the hinge positions, the moving mesh zone, the height of the sea floor to adapt the dissipation parameter and the body motion solver. In this work the body motion solver described above is used exclusively but other body motion solvers can be used to perform forced oscillation tests for example. Fig. 4 shows two instances during a wave cycle. The flap shape is shown with a longitudinal slice of the tank to illustrate the mesh motion. The sea floor is represented by high dissipation values and can be seen to change inside the moving cylinder while it rotates. This means that the mesh resolution around the bottom must be sufficiently high and the value for the dissipation variable must be set to a high enough value.
Simulations were run for two different mesh refinement levels in the rotating cylinder. Refinement levels close to the flap and in the outer, static mesh are identical, while the rest of the moving cylinder was refined once more, that is all edges were split into half. Fig. 5 shows the rotation angle over time for the coarse and fine meshes. The simulation with the fine mesh shows about 11 larger rotation amplitudes of the flap. The shape and frequency of the rotation traces agree well.
Results of simulations for different values of the dissipation value under the floor level are shown in Fig. 6. The maximum rotation angle for the case with a dissipation coefficient of zero, that is without taking into account the sea floor inside the rotating cylinder, is about 10% or 31 less than for the two cases run with values of 50 or 100. A phase-shift can also be observed. The flap reaches its maximum rotation angle earlier when the floor is not considered, this difference increases over the wave period T displayed. No difference between the two later cases can be observed, all future cases were run with a value of 50.
The accuracy of the solution is affected in two ways by the choice of time step. The solution of the flow field and the solution of the equation of motion of the moving flap are both time step dependent. Only the solution of the flow field is physically related to the Courant number. The accuracy of the solution of the equation of motion can thus not be deemed sufficient for all cases, only because the flow field is solved correctly. For example, a configuration in which the flow velocities are low but the accelerations of the flap high, the time step might be too large for the motion solver. It seems though, that in general the high velocities around the top of the flap and quickly moving fluid interfaces constrain the time-step more than the equations of motion. Fig. 7 shows the rotation angle over time for simulations performed with different Courant numbers. Results show very little variation for Courant numbers smaller than 0.3. In all following simulations a Courant number of 0.2 was used.

Experimental setup
The following section describes the experiments performed in the wave tank at Queen's University Belfast to create data specifically for the comparison with numerical results.
The wave tank at Queen's University's hydraulic laboratory is 4.58 m wide and 20 m long. An Edinburgh Design Ltd. wave-maker with 6 paddles is installed at one end. The bottom is made of two horizontal sections connected by sloped concrete slaps which allow experimental testing 150 mm and 356 mm above the lowest floor level at the wave-maker. A beach consisting of wire meshes is located at the opposite end. An over-view of the bathymetry and the flap location in the experiments described can be seen in Fig. 8.
The flap measures 0.1 m Â 0.65 m Â 0.341 m in x, y and z directions.
Water-levels in the tank are defined with reference to the deepest point in the tank, at the wave-maker.
The flap model consists of three units, the fixed support structure, the hinge and the flap, Fig. 9.
The support structure is made of a 15 mm thick, stainless steel base plate, measuring 1 by 1.4 m, which is fixed to the bottom of the tank by screws. The hinge is held in three bearing blocks. To accommodate an electric drive above water, which was not utilised in the physical experiments shown here, a platform with three cylindrical legs is mounted beside the flap.
The flap itself is made of a foam centre piece, sandwiched by two PVC plates on the front and back face. Three metal fittings connect the flap to the hinge, enabling changes of the flap even without draining the tank.
A 3 axis accelerometer from Kistler, Type 8395A010ATT00 was attached onto the top of the flap. The sensor has a range of 7 10 g. Only the y and z channels were used. With the sensor attached to the top of the flap one channel gives radial a rad , the other tangential accelerations a tan .
It should be noted that accelerations in different directions are measured in slightly different positions inside the sensor. An offset of 4 mm is irrelevant when the complete radius of gyration, that is the distance from hinge to sensor position, is 324.5 mm. To simplify post-processing only this one value was used and it was assumed that the sensor positions are directly above the centre of the hinge.
Mass and inertia data were extracted from 3D CAD models as follows: The wave-probes are standard resistance wave gauges, according to Masterton and Swan (2008) and can also be assumed to be accurate within 70.5 mm.
The accuracy of the accelerometer was not independently assessed. The calibration certificate states a transverse sensitivity of 3% for all three channels. The largest uncertainty is understood to stem from the dynamic bearing friction. Although not directly determined, the value can be assumed to be slightly less than the static bearing friction, which was derived as follows: The flap was left in an upright position (without water in the tank) within the range of about 11.
From the weight and the position of the centre of gravity relative to the hinge it can be calculated that the (static) bearing friction is about 0.01 Nm.
According to numerical results the total hinge moment amplitude in the monochromatic seas is about 1 Nm. The expected error due to bearing friction losses is thus only about 1%.
In the wave series tests simulating the random waves the moment amplitude obtained from numerical simulations is mostly around 0.4 Nm. However at t ¼ 14 s it drops to less than 0.2 Nm. Thus the bearing friction could be a significant part of the total measured value in the physical experiment.

Results
First simulations were run for 20 s in monochromatic seas with a period of 2.0625 s and an amplitude of 0.038 m. This equates approximately a wave of 13 s period and 1.5 m wave height at 40th scale, taking into account the clocking rate of the wave maker. The Ursell number as defined by Fenton (1998)  Numerical results show the start up phase from still water. The second wave crest (5 s) is slightly higher than the preceding, after that the surface elevation settles into a regular pattern with almost constant wave amplitude. While the crest has a smooth sinusoidal shape all troughs indicate some perturbation.
Experimental data shows some slight noise before the first crest. The second crest is the highest in the wave trace, similar to the numerical results. The experimental data shows a distinct drop in the third trough which is not replicated in the numerical data, all following waves have a flat crest. The troughs are always deeper and the crests lower compared to the numerical data. It seems as if a reflected or radiated wave superimposes the original incoming wave. The zero crossing periods match very well. Fig. 10B shows the tangential and radial acceleration component in the accelerometers frame of reference. Numerical rotation data was used to obtain the acceleration components equivalent to the raw experimental results. The skill value as defined by Dias et al. (2009) is a suitable metric to compare the accuracy of numerical models. A value of one would indicate perfect agreement or identical signals. Comparison of numerical and experimental traces yield the following:    The relatively low skill value for the radial acceleration is mainly due to the high frequency noise of the experimental signal, which is not present in the numerical data and obviously not a real feature of the flap motion.
The radial acceleration caused by the flap motion acts against gravity. When the flap moves, radial acceleration drops from the starting level of 9:81 m=s 2 and after settling oscillates with a small amplitude of about 0:5 m=s 2 over each wave cycles around an average of 7:5 m=s 2 .
The tangential acceleration shows much larger amplitudes of up to 7 m=s 2 . The crests show very good agreement in shape and amplitude between numerical and experimental data. Some differences can be observed in the shape of the troughs. While the crests are round, the troughs show a little dip when reaching the highest negative acceleration, the signal than flattens out before rising again. The flat part is much more pronounced in the numerical data, the amplitudes of negative acceleration agree very well between numerical and physical data.
As a second test case a series of waves of similar but varying amplitude and frequency were calibrated in the physical tank, results are shown in Fig. 11. The plot shows results in the same way as previously in Fig. 11.
The skill values are as follows:      amplitude waves around 10 s show some difference, there and at the very beginning and end of the trace high frequency oscillations can be seen on the experimental data. Acceleration data compares very well over all. At around 11 s the numerical data shows higher accelerations. As some friction had been observed in the bearings during the experiments, it seems reasonable to assume that these will be more dominant when exciting forces are smaller, that is the case with smaller wave heights, when these discrepancies occur. As in the monochromatic cases high frequency noise can be observed on the experimental acceleration signals, reducing the skill value especially for the radial acceleration.

Conclusions
A new way of simulating OWSC's in a mesh based RANS CFD code was presented and the solver tested against two experimental benchmark tank tests. The following conclusions can be drawn: The numerical methods presented in Section 2 enable the simulation of an OWSC in normal operating conditions.
The method of using a cylindrical mesh rotating around the hinge point enables efficient simulation of a moving flap.
Modelling the bed with a spatially fixed dissipation zone represents the sea floor well.
Solution of the equations of motion using three weights for smoothing is stable even in significant waves.
Differences between numerical and experimental data are believed to be caused primarily by differences in exciting waves, i.e. reflections and other perturbations.
Further errors are believed to stem from the friction of the bearings used in the experiments.