A CFD simulation of wave loads on a pile-type oscillating-water-column device ＊

Zhenhua Huang, Conghao Xu, Shijie Huang

Department of Ocean and Resources Engineering, School of Ocean and Earth Science Technology, University of Hawaii at Manoa, Honolulu, USA

Abstract: Oscillating-water-column (OWC) devices provide a practical and convenient way to convert wave energy to electricity.Wave loads on bottom-sitting WECs are important design parameters for their foundation safety. This study considers an OWC integrated into a bottom-sitting pile structure (an OWC pile). Results from three-dimensional numerical simulations are used to examine the wave loads on an OWC pile. It appears that viscous loads and lateral loads are relatively small. Both the inline and vertical wave forces and the resulting overturning moment are not sensitive to the variability in wave period, but increase with wave height.Scales suitable for the wave forces and overturning moment are proposed for estimating wave loads on the OWC pile.

Key words: Numerical simulation, wave energy conversion, OpenFOAM, foundation safety

Introduction

Foundation design of a bottom-sitting wave energy converter (WEC) begins with site-specific wave loads on the structure, including wave force and the resulting overturning moment. Most existing studies of wave energy converters have focused on wave power extraction efficiency, few studies have been dedicated to wave loads on WECs.

The ocean wave is one of the most promising renewable energy sources. Many types of WECs have been proposed in the past to extract energy from ocean waves for electricity generation. Oscillating water columns (OWC) are one type of the WECs that have been widely investigated and tested so far. A typical OWC-type WEC consists of an air chamber above the water surface with one turbine (power take-off device)mounted to the air chamber for extracting wave energy. The bottom of the chamber is open and submerged so that the incident waves can cause the water column inside the chamber to oscillate, which generates a fluctuation in the air pressure inside the chamber, causing the air flow to drive the power take off (PTO) device for wave energy extraction. One advantage of the OWC-type WECs is its simplicity of wave energy extraction mechanism.

Many studies on the hydrodynamics and energy extraction of various OWC devices can be found in the literature. The methods used in these studies include wave-flume tests, frequency-domain analysis based on potential flow theory, numerical simulations based on potential flow theory, and computational fluid dynamics (CFD) simulations. In experimental studies and CFD simulations, the power-take-off(PTO) is usually modeled by an orifice. Existing experimental studies have focused mainly on the conversion efficiency of various OWC devices[1-6],hydrodynamic characteristics such as vortex shedding,wave scattering and motion responses[2-5,7-8].

A limited number of studies on wave loads on wave energy converters can be found in the literature.These include a theoretical analysis of wave loading on a 2-D heaving WEC (The Berkeley Wedge)[9], an experimental study of wave forces on a fixed 2-D OWC device[10], an experimental study of wave loads on a hemispherical point absorber (Wavestar)[11], a CFD simulation of the wave force on a fixed 2-D rectangular OWC device[12], a CFD simulation of the wave forces on a 3-D buoy-type WEC[13]. For 2-D wave energy converters, which can also serve as breakwaters, there exists a correlation between the reflection coefficient and the dynamic horizontal wave force and a correlation between the air pressure in the OWC chamber and the vertical dynamic loading[12].Dynamic horizontal force is usually much larger than the vertical dynamic force[10,12].

CFD simulations of OWC devices can provide important information about the complex flow field around an OWC device and the spatial distribution of the water surface inside the OWC chamber, which is otherwise very difficult to obtain in laboratory experiments. Most of existing CFD simulations of OWC devices focused on 2-D configurations such as a rectangular shaped OWC[14]. We present a numerical study on the wave loads on a circular OWC-type WEC integrated into a large pile[4]. The numerical model is an OpenFOAM-based two-phase flow model and the PTO is modeled using an orifice mounted to the top cover of the OWC chamber. OpenFOAM is an open-source CFD library for solving CFD problems using a volume of fluid method to track the free surface[15].

1. Description of numerical model

1.1 The pile-type OWC model

The oscillating-water-column device studied here is a circular OWC supported by a C-shaped structure[4]. See Fig. 1 for a 3-D view of the OWC device.The model has an overall dimension of 40 cm in total height. The distance from the lower tip of the OWC chamber to the flume bed is DS= 244 mm. The inner diameter of the cylinder is D =1 25 mm and the thickness of all walls is 3 mm. On the top of the tube sector an orifice of diameter DO= 14 mm is used to simulate a nonlinear power takeoff (PTO) device. The opening ratio of the orifice area to the cross-sectional area of the air chamber is ε = 1.25%.

Fig. 1 A sketch of the circular OWC model. Not drawn to scale

To mathematically model the OWC pile, a Cartesian coordinate system is defined. Referring to Fig. 1, the x- axis points in the direction of wave propagation, the y- axis points vertically upward,and the direction of the z- axis is obtained by the right-hand rule. y- axis coincides with the line of symmetry of the air chamber. The origins of all three coordinates are at point O, which is on the seabed.

1.2 Numerical wave tank

A relaxation-zone-based wave-generation method[15]is used for wave generation and absorption. In this method, the computational domain is divided into three sections: a wave generating relaxation section, a test section and the wave-absorbing relaxation section(numerical beach), see Fig. 2 for the numerical setup and the three sections. In this study, the length of the wave generating relaxation section is 4 m long, the test section is 6 m long, and the wave absorbing section is 4 m long. Numerical tests have shown that the reflection coefficient from the wave absorbing relaxation zone is less than 0.05 for the waves examined in this study. The turbulence model used is a k- ω model. The air-water surface is tracked using a VOF method. The thin wall of the OWC chamber and the high velocity through the orifice requires very fine mesh in the vicinity of the orifice and the thin walls. As shown in Fig. 3, the simulated surface displacement of the stabilized waves agrees well with that of linear wave theory, confirming the quality of the waves generated in the numerical wave tank.

In the wave generating relaxation zone, the value of of a variable φ (e.g., velocity u or fluid saturation α) is a blend of the actual computed value by the RANS formulationCφ and the value solved theoretically using a wave theorytφ , with a relaxation factorRα . Therefore, the equation for φ in the wave generation zone can be written as

The value of the relaxation factorRα is computed by

The numerical wave tank was set up with the boundary conditions resembling the actual test conditions[4]. These conditions are: wave inlet boundary conditions at x =0 m, wall boundary conditions on all rigid walls, and atmospheric boundary condition at the top boundary of the computational domain. All wall boundaries are assumed to be hydraulically smooth because the walls of the wave flume in the experiment are glass and the walls of the OWC model are stainless steel.

Fig. 2 Numerical setup of the numerical wave tank, not drawn to scale. 1G, 2G and 3G are the three locations where waves are measured by wae gauges

Fig. 3(a) (Color online) Times series of the waves measured at three locations in the numerical wave tank with x being the distance from the left end of the numerical wave tank

Fig. 3(b) (Color online) A snapshot of the waves in the numerical wave tank with red line being the simulation result and the black line being the theoretical airy wave solution

An unstructured nested mesh is used to cover the entire computational domain. In order to reduce the computational load, a coarse mesh of resolution 16.4 mm×10 mm×50 mm ( Δx ×Δ y ×Δ z) is used in the wave generation and absorption zones. In the test section(from x= 4 m to 10 m), the mesh is refined to a resolution of 8.2 mm×10 mm×25 mm. The mesh within ±50 mm around the still water surface is further refined to are solution of 8.2 mm×2.5 mm×12.5 mm.In the vicinity of the OWC model, a mesh of 20.5 mm×12.5 mm×15.9 mm is used. near the thin walls of the OWC pile mode, the mesh resolution is further refined to 1 mm×1 mm×1 mm, which closely snaps the surface of the OWC model. The total cell count of the final mesh is about 1.4×106.

Fig.4 (Color online) Comparisons between the simulated and measured surface displacements at three locations and the air pressure. From top to bottom, the plots show the measured surface displacements at the locations 1G,and 3G and the relative air pressure. The circles are measurement and the solid lines are numerical simulations

2. Model verification and validation

Representative results are presented in this section to show the agreement between the simulated and measured results. These results include surface displacements at three locations ( 1G,2G and3G as shown in Fig. 2) and the relative pressure of the air in the pneumatic chamber.

Comparisons between the time series of the measured and simulated relative pressure and the time series of the measured and simulated surface displacements at three locations are shown in Fig. 4(a)for wave period T = 0.8 s , wave height H=0.04 m and water depth h = 0.31 m, and in Fig. 4(b) for wave period T =1.3 s, wave height H = 0.04 m ,and water depth h = 0.29 m. The simulated and measured surface displacements agree very well at all three locations.

The relative pressure of the air in the pneumatic chamber, which is basically the pressure drop across the PTO, is not sinusoidal, even though the incident waves are. It can be modeled by

where u (t) is the air flow velocity,aρ the air density, cfa quadratic loss coefficient, andgL an empirical length scale related to inertia effect[4]. For sinusoidal waves, the relative pressure p (t) in the air chamber is not sinusoidal, which means the PTO used in the study is nonlinear. As shown in Fig. 4, the numerical model can simulate the relative pressure of the air inside the pneumatic chamber very well,implying that the numerical simulation can simulate the behavior of this nonlinear PTO very well.

The bottleneck in simulating OWC-type WECs is associated with the high air velocity through the PTO.The air flow velocity through the orifice is mainly controlled by the opening ratio ε, which is 1.25% for the OWC-pile WEC studied here. Even though the wave height is normally small, the maximum velocity of the air flow through the PTO device can be amplified by a factor of 1/ε due to conservation of mass for incompressible fluids. For waves of a height H = 0.04 m and a wave period T =1s, the velocity of the water surface inside the air chamber can be estimated by 2 πH / T =0.125 m/s , for which the air flow velocity through the orifice can reach about 10 m/s.The above estimation of the air flow velocity has been confirmed by our numerical simulations (see the top panel of Fig. 6 for an example). The fine mesh and high velocity speed in the vicinity of the PTO requires very small-time step, which creates a bottleneck in 3-D simulations of OWC devices.

3. Results on wave loads acting on the OWC pile

3.1 Calculation of pressure force and overturning moment

The air pressure inside the OWC air chamber does not contribute to the horizontal wave force and the resulting overturning moment because the air pressure can be regarded as spatially uniform inside the air chamber. Like in wave loads on slender cylinders, the viscous components of the wave loads are orders of magnitude smaller than the pressure components, and thus can be safely ignored in the wave load analysis. Using the pressure p (t) for the fluid phase (air and water) provided by the OpenFOAM simulation, the wave force F and the overturning moment M on the structure are calculated by:

and

where S consists of the outer and inner surfaces of the structure, ndA is the differential surface vector and r is the length of the arm of rotation at the location (x , y , z), which is the location of ndA . The arm of rotation is determined using a center of rotation set at the center of the circular tube at y=0 . The wave force F (t) has three components,with l, J and k being the unit vectors in the x, y and z directions, respectively.is the inline force,zF is the lateral force andyF is the vertical force. Similarly, the overturning moment M(t) has three components and can be expressed as

Both F(t) , M(t) are calculated using an OpenFOAM built-in function object “forcesIncompressible” available in version 4.x. If the simulation is performed using FOAM-Extend 3.2 environment, the data format needs to be converted to OpenFOAM 4.x format in order to use the function object “forcesIncompressible” available in OpenFOAM 4.x.

3.2 Hydro-static loads

In the absence of waves, the loads on the OWC pile is purely the hydro-static loads due to buoyancy.The static force has only the vertical component,which comes from the pressure acting on the surface of the lower tip of the OWC chamber. For h = 0.31 m, the theoretical value ofis 0.37 N. In numerical simulations, the mesh size will affect the accuracy of the computed values of the hydro-static loads. Let the numerical error of the wall thickness be the mesh size, the numerical error in the surface area of the lower tip is estimated as δ A =πr δ r , where r is the inner diameter of the OWC chamber and rδ is the numerical error in the wall thickness. For the present OWC model and a mesh size of 0.00125 m,which gives δ A =2.45 × 10-4m2. The draft is Dr= 0.06 m for h = 0.31 m , which results in a possible numerical error of 0.144 N. Therefore, the mean vertical force should be in the range of 0.37 ± 0.1 44 N . For the model in still water, the simulatedis 0.44 N for the mesh used in this study, which is in the range of 0.37 ± 0.1 44 N . The hydro-static overturning moment is due toThis overturning moment is in the z-direction, and denotedasThe simulated value ofis-1.19×10-2N⋅m. Based on the estimation of the numerical error in, the numerical error inshould be about 2.78×10-3N⋅m.

The focus of this study is on the wave loads,which exclude the hydro-static loads. In the following, the hydro-static forces and moments are not included in the results for wave-induced forces and overturning moments.

Fig. 5(a) (Color online) Time series of the calculated total pressure force on the OWC-pile ( T = 1.4 s , H =0.04 m and h = 0.31 m )

3.3 Characteristics of wave loads

The three components of the force and overturning moment acting on the OWC-pile are shown in Fig. 5 for wave period T = 1.4 s , wave height H=0.04 m and water depth h = 0.31 m . The inline wave forcexF dominates the total wave force.

Fig. 5(b) (Color online) Time series of the corresponding overturning moment about the sea bed ( T = 1.4 s , H=0.04 m and h = 0.31 m )

The inline wave forcexF is much larger than the other two components. This force has a period same as the incident wave period. The air pressure inside the air chamber does not contribute toxF because the air pressure is spatially uniform. The maximum is 2.75 N, while the minimum -2.55 N.However, the mean ofxF is almost zero because the duration of positivexF is slightly different from that of negativexF . The asymmetry ofxF can be attributed to the asymmetry of the supporting structure of the OWC pile.

The vertical wave forceyF has a zero mean because the non-zero hydro-static force is not included.In general, the vertical wave force on the OWC pile can be decomposed into two components

The lateral hydrodynamic forcezF is at least two orders of magnitude smaller than the inline force.Physically,zF is directly related to the vortex shedding from the OWC pile, which may have a frequency higher than the wave frequency. Vortex shedding from the OWC pile is related to a Keulegan-Carpenter(KC) number defined by KC = UT/ D with D being the diameter of the OWC chamber, T the wave period and U a velocity scale. If we take U=ωH/ 2, we have KC = πH / D . For H = 0.04 m,T = 1s, we obtain KC =1, which is smaller than the critical value of 7, meaning the vortex does not shed from the outer surface of the OWC chamber and its main support structure. The small lateral force is related to the vortex shedding from sharp edges of the supporting structure, which is related to a wall KC number defined by KCwall= UT/ d with d being the thickness of the wall of the supporting structure.For the present problem, this wall KC number is much larger than 7, which guarantees vortex shedding from the sharp edges of the supporting structure. The above analysis of vortex shedding is confirmed by the computed vorticity distribution shown in the bottom panel of Fig. 6, which shows that the vortex shedding mainly occurs at the lower tip of the OWC skirt and there is weak vortex shedding at the sharp edges of the support structure (see Fig. 6). This explains why the lateral force due to vortex shedding is weak compared to the inline wave force.

Fig. 6(a) (Color online) An example of the simulated air velocity and the air-water interface ( T = 1.0 s , H=0.04 m and h = 0.31 m )

Fig. 6(b) (Color online) An example of the simulated vorticity( T = 1.0 s , H = 0.04 m and h = 0.31 m)

The z-component of the overturning moment Mzis much larger than the other two components.BothxF,yF contribute to Mz. The vertical wave force contributes to the overturning moment Mzbecause of the geometry of the OWC pile. The lateral wave forcezF is at least two orders of magnitude smaller than the inline wave forcexF and the vertical wave forceyF, and thus can be ignored when analyzing the overturning moment. In the following,we will only discuss the z-component of the overturning moment.

3.4 Changes of the pressure force and overturning moment with wave period

The changes of the pressure force and overturning moment with wave period are demonstrated by fixing the wave height at 0.04 m and varying the wave period from 0.7 s to 1.6 s in water of depth h = 0.31m. Hereinafter, all results will be presented using dimensionless variables defined below.

Because wave loads are due to the wave-induced pressure, which is scaled by ρwgA with ρwbeing the density of water and A the wave amplitude, we normalize the wave force by ρwg HDh/2 with D being the diameter of the OWC chamber and h the water depth. Similarly, we normalize the overturning moment by ρwgH Dh /2 × H /2. Wave period will be normalized byand wave height H will be normalized by the diameter of the OWC pile D.Therefore, the normalized wave force and overturning moment are denoted by:

Changes of the maximum and minimum of the inline wave forceand the vertical wave forcewith dimensionless wave periodare shown in Fig. 7(a). Changes of the maximum and minimum of the overturning moment mz(t) with dimensionless wave periodare shown in Fig. 7(b).

As shown in Fig. 7(a), the inline wave force is about several times as large as the vertical wave force.Both the maximum and minimum values of the inline and vertical wave forces are not very sensitive to the variation in wave period within the tested range, even though the wave forces are slightly smaller for long waves. This means thatis a reasonable scale for estimating both the inline and vertical wave forces.

Fig. 7(a) (Color online) Change of maximum and minimum pressure forces with ( H = 0.04 m, h=0.31m, lines show averages)

Fig. 7(b) (Color online) Changes of maximum and minimum overturning moments with ( H = 0.04 m,m, lines show averages)

The overturning moment mz(t) shown in Fig.7(b) contains contributions from both fx(t ) and fy(t ). For practical purpose, the variations of max( mz), min( mz) withare mild andis a reasonable scale for estimating the overturning moment. An average value of the overturning moment can be used for a first-cut estimation of overturning moment.

The difference between the maximum and minimum values of both wave forces and overturning moment are within the range of the variation of the wave forces and overturning moment withIf we neglect the variation of the wave loads withwe should also neglect the difference between the maximum and minimum values of the wave loads. We define the magnitude of each normalized wave load as:

It can be hypothesized that the magnitudes of the wave forces and overturning moment can be estimated by:

where the coefficients α, β and γ are basically average values ofandrespectively. Knowing the coefficients α, β and γ helps estimate the wave loads on an OWC pile for a first-cut design of the OWC pile foundation. The coefficients α, β and γ are expected to be dependent of H/ D due to the weakly nonlinear nature of the wave loads. For the data in Fig. 7, we have α= 0.378 ± 0.0 14, β = 0.106 ± 0.0 09 and γ =0.448 ± 0.0 23.

For some WECs, there is a correlation between wave loads and wave energy extraction efficiency[9].For the WEC studied here, the wave loads given in Fig. 7 do not show peak values at certain wave lengths as the wave power extraction efficiency does[4],indicating that energy removed by the nonlinear PTO-simulating orifice is not a main factor affecting the wave loads on the OWC pile. Because the wave loads are not directly related to wave power extraction efficiency which is affected by the nonlinear nature of the PTO-simulating orifice, they are not likely affected by the air flow through the orifice in any significant way. i.e., the wave loads seem to be affected more by wave scattering than by wave radiation.

3.5 Changes of the pressure force and overturning moment with wave height

Fig. 8(a) (Color online) Changes of maximum and minimum pressure forces with H /D ( T=1.0 s, h = 0.31m )

Fig. 8(b) (Color online) Changes of maximum and minimum overturning moments with H/ D ( T=1.0 s, h=0.31m)

Changes of the pressure force and overturning moment with wave height are demonstrated by fixing the wave period at 1.4 s and varying wave height from 0.02 m to 0.05 m in water of depth 0.31 m.Changes of the maximum and minimum of the inline wave force fx(t) and the vertical wave force fy(t ) with dimensionless wave height H /D are shown in Fig. 8(a), changes of the maximum and minimum of the overturning moment mz(t ) with dimensionless wave height H /D are shown in Fig.8(b).

As shown in Fig. 8(a), the inline wave force is several times as large as the vertical wave force. For both the inline and vertical wave forces, they increase with increasing wave height. For weak waves, the vertical wave force is almost zero. The overturning moment also increases with increase wave height, as shown in Fig. 8(b). The coefficients α, β and γ increase with increasing wave height because they are the average values ofanddefined in Eqs. (10)-(12).

where k is the wave number, ηmaxis the maximum surface displacement and umaxis the maximum wave orbital velocity[16]. The terms in the square brackets on the right-hand side of Eq. (16) represent the effect of wave direction. For the present problem of a slender OWC pile, kD /2 ≤ 1 and thus the diffraction effect is not important. As a result, the wave runup can be estimated by

If we adopt the following estimates for ηmax,andwe have the following estimate for the wave runup on the OWC pile

For a fixed wave height H = 0.04 m, T=0.7-1.6 s and R≈ 0.0 2 m, which means the wave runup is not sensitive to wave period for a fixed wave height. Therefore, wave period does not have a significant effect on the wave loads through wave runup, as shown in Fig. 7. When the wave period is fixed at T = 1.4 s , however, the wave runup on the OWC pile can be estimated by

which means the wave runup increases with wave height. Because larger wave runups generate larger wave loads, larger wave forces and overturning moments are expected for larger wave heights, as shown in Fig. 8. The wave runup effects on wave loads cannot be revealed by frequency domain solvers based on linear potential theory.

Recently Xu and Huang[17]experimentally examined a row of OWC piles, with each OWC pile be identical to that studied here. CFD simulations based on OpenFOAM are underway and the results will have reported elsewhere.

4. Conclusions

Wave loads on a bottom-sitting circular pile integrated with an OWC have been studied using OpenFOAM simulations. The following main conclusions can be drawn: (1) The inline wave force is several times as large as the vertical wave force, (2)The dominant overturning moment is due to the inline wave force and the vertical wave force, (3) The lateral wave loads are relatively weak and can be ignored, (4)The magnitudes of both the wave forces and the overturning moment are not very sensitive to variations in wave period, but increases with increasing wave height. It seems that the wave forces can be scaled byand and the overturning moment can be scaled by

Acknowledgements

This work was supported by the US National Science Foundation (Grant No. CBET-1706938), the Extreme Science and Engineering Discovery Environment (XSEDE) (Grant Nos. OCE170015,ENG180008). XSEDE is supported by National Science Foundation (Grant No. ACI-1548562). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. This is SOEST contribution number 10621.