Numerical simulation of floating bodies in extreme free surface waves

In this paper, we use the in-house Computational Fluid Dynamics (CFD) flow code AMAZON-SC as a numerical wave tank (NWT) to study wave loading on a wave energy converter (WEC) device in heave motion. This is a surface-capturing method for two fluid flows that treats the free surface as contact surface in the density field that is captured automatically without special provision. A time-accurate artificial compressibility method and high resolution Godunov-type scheme are employed in both fluid regions (air/water). The Cartesian cut cell method can provide a boundary-fitted mesh for a complex geometry with no requirement to re-mesh globally or even locally for moving geometry, requiring only changes to cut cell data at the body contour. Extreme wave boundary conditions are prescribed in an empty NWT and compared with physical experiments prior to calculations of extreme waves acting on a floating Bobber-type device. The validation work also includes the wave force on a fixed cylinder compared with theoretical and experimental data under regular waves. Results include free surface elevations, vertical displacement of the float, induced vertical velocity and heave force for a typical Bobber geometry with a hemispherical base under extreme wave conditions.


Introduction
In this paper, we describe developments of the AMAZON-SC 3D numerical wave tank (NWT) to study extreme wave loading of a floating structure (in Heave motion).The extreme wave formulation prescribed as an inlet condition is due to Dalzell (1999) and Ning et al. (2009) is based on a first or second-order Stokes focused wave.
Correspondence to: Z. Z. Hu (z.hu@mmu.ac.uk)The AMAZON-SC 3D code (see e.g.Hu et al., 2009) uses a cell centred finite volume method of the Godunov-type for the space discretization of the Euler and Navier Stokes equations.The computational domain includes both air and water regions with the air/water boundary captured as a discontinuity in the density field thereby admitting the break up and recombination of the free surface.Temporal discretisation uses the artificial compressibility method and a dual time stepping strategy to maintain a divergence free velocity field.Cartesian cut cells are used to provide a fully boundary-fitted gridding capability on a regular background Cartesian grid.Solid objects are cut out of the background mesh leaving a set of irregularly shaped cells fitted to the boundary.The advantages of the cut cell approach have been outlined previously by Causon et al. (2000Causon et al. ( , 2001) ) including its flexibility for dealing with complex geometries whether stationary or in relative motion.The field grid does not need to be recomputed globally or even locally for moving body cases; all that is necessary is to update the local cut cell data at the body contour for as long as the motion continues.The handing of numerical wave paddles and device motion in the AMAZON-SC NWT is therefore straightforward and efficient.
Firstly, extreme design wave conditions are generated in an empty NWT and compared with laboratory measurements as a precursor to calculations to investigate the survivability of the Bobber device operating in a challenging wave climate.Secondly, a fixed submergence horizontal cylinder has been validated in regular waves.Finally, a floating Bobber has been simulated under extreme wave conditions.

The Cartesian cut cell mesh
In this paper, the Euler equations version of AMAZON-SC willbe extended to handle a 3-D floating body in extreme waves.The majority of the flow domain is overlaid with Published by Copernicus Publications on behalf of the European Geosciences Union.520 Z. Z. Hu et al.: Numerical simulation of floating bodies a regular Cartesian mesh.Moving flow boundaries or bodies in relative motion within the flow domain are accommodated by computing local cell-boundary intersections at the boundaries that are in motion.A detailed description of the principles of the cut cell method including the numerical procedures applied at solid boundaries that are either static or in relative motion have been given previously by the authors (Yang et al., 1997;Causon et al., 2000Causon et al., , 2001;;Qian et al., 2003) including extension of the cut cell algorithms to 3-D (Yang et al., 2000;Hu et al., 2009).

The flow solver on a cut cell mesh
The Euler equations for a general moving control volume of fluid expressed in a pseudo-compressible form can be written as where Q is the vector of conserved variables within the volume V , F is the conservative flux through the volume's bounding surface S, whose outward-pointing unit-normal vector is n and B is a source term of body forces.Q, F and B are given by where where u, v and w are the flow velocity components and u b , v b and w b are the velocity components of the (body) boundary S, which are identically zero when the boundary is stationary.ρ is the fluid density, p is the pressure, β is the coefficient of artificial compressibility and g is the gravitational acceleration.The use of a pseudo-compressible form of the describing equations and artificial compressibility parameter β permits the use of efficient modern Riemann-based solution methods developed for compressible flows.
We then discretize Eq. ( 1) over each cell of the flow domain using a finite volume formulation, which gives where Q ij k is the average value of the solution vector Q at cell (i, j , k) stored at the cell centre and V ij k denotes the volume of the cell.F k is the numerical flux across the bounded face k of the finite volume cell, A k is the area of that face and r is an integer identifier for each face of the cell (> 4 in the case of a 2-D cut cell).
The convective flux F k in Eq. ( 5) is then determined by solving a Riemann problem at each cell interface.This involves two stages: first, the left and right state values are reconstructed on the opposite sides of each cell face by projecting the solution data from the stored cell centre values either side of the cell face centre; second, the resulting Riemann problem defined by the left and right state data is solved using an approximate Riemann solver (ARS).
Thus, the required left and right state values corresponding to stored cell centre data Q(x,y,z) can be found anywhere within the cut cell using where r is the normal distance vector from the cell centre to any specific interface or solid boundary.Q i,j,k is the stored computed cell centre data and Secondly, Roe's approximate flux F k is constructed to second order accuracy in each cell.This assumes a 1-D Riemann problem in the direction normal to the cell face and has the form where Q R and Q L are the reconstructed data values on the right and left of face k based on Eq. ( 6) and Ã is the flux Jacobian evaluated using Roe's average state values The quantities R and L are the right and left eigenvectors of Ã.
is the diagonal eigenvalue matrix, where the eigenvalues are given by λ 1,2,3 = un x + vn y + wn z ,λ 4,5 = 0.5 un x + vn y + wn z ± C where C = un x + vn y + wn z

2
+ 4β/ρ.The Jacobian matrix A can be constructed according to its definition Then Q is introduced into Eq.( 9) and the Jacobian matrix Ã constructed.
To achieve a time-accurate solution at each time step of an unsteady flow problem, a first-order Euler implicit difference scheme is used for the time discretisation in Eq. ( 5) as follows, Here, the incompressible flow is numerically modelled by employing the well known artificial compressibility method.Thus, a pseudo time derivative term is added to the equation set in such a manner that the true incompressible flow equations are recovered at each physical time step by a process of marching to convergence in pseudo time until the time derivative term is driven to zero.This gives where τ is the pseudo-time and If the foregoing linearizations are introduced into Eq.( 11), the result can be written as where m is iterated to zero at each physical time step the density and momentum equations are satisfied identically and the divergence of the velocity at time level n+1 is zero as required.The system of equations can be written in matrix form as RHS stands for the right side of Eq. ( 12), D is a block diagonal matrix, L is a block lower triangular matrix and U is a block upper triangular matrix.Each of the elements in D, L and U is a 5 × 5 matrix.An approximate LU factorization (ALU) scheme as proposed by Pan and Lomax (1988) can be adapted in the form wherein within each physical time step of the implicit integration process the pseudo time sub-iterations are terminated when the L 2 norm associated with the iteration process is less than a user specified limit ε, and ε = 10 −4 here.Further details may be found in Qian et al. (2006).

The extreme wave formulation
As is well known, the exact velocity profile for a true physically realisable nonlinear wave under given conditions is not known a priori.Thus, a viable approach is to input reasonable approximate wave conditions along the input boundary to simulate the real phenomenon.This leads to the notion of the extreme wave formulation as a focused wave group in which many wave components in a spectrum are focussed simultaneously at a position in space in order to model the average shape of an extreme wave profile.The derivation here refers to the work of Dalzell (1999) and Ning et al. (2008Ning et al. ( , 2009) ) in which a first or second-order Stokes focused wave can be imposed in such a manner.
A Cartesian coordinate system O-xyz is defined with the origin located at the undisturbed equilibrium free surface, with the z-coordinate vertical and positive upwards.The xcoordinate is zero at the wave-maker located at x = 0.0 m, x 0 is the focus point, t 0 is the focus time and the water depth h.For each wave component i the amplitude A i can be calculated from where N is total number of components, S i (f ) is the spectral density, f is the increment in frequency depending on the number of wave components and band width and a i is the input wave amplitude of the focused wave.
www.nat-hazards-earth-syst-sci.net/11/519/2011/Nat.Hazards Earth Syst.Sci., 11, 519-527, 2011 The corresponding wave elevation η, and horizontal and vertical velocities u and w are expressed as follows: where η (1) , u (1) and w (1) are the linear wave elevation and velocities, η (2) , u (2) and w (2) correspond to the second-order wave elevation and velocities, respectively.Both velocity and wave elevation can be decomposed into N components with different frequencies as follows: where g is the gravitational acceleration, h is the water depth, the wave number k i = ω 2 i /gtanh(k i h), the frequency ω i = 2πf i , the phase slope ε is set to zero for the calculations in this work and At the interface between two immiscible fluids, the present method assumes that the system of equations for nonhomogeneous incompressible flow can treat the free surface numerically as a contact discontinuity in the density field.
The need for special procedures to track the free surface is thus eliminated, since the free surface is captured automatically in the time-marched numerical solution as a discontinuity in the density field.It is asserted that the numerical solution of Eq. (1) for a system containing one or more free surfaces will converge to a physically-correct unique solution.
In this paper, the pressure p can be obtained from p/β in Eq. ( 2).The total force is obtained by integration of the pressure field around the body contour F = − S b pndS, where S b is the body surface as defined approximately by the fitted cut cell surface.

Results
In the following sets of results, the density ratio between water and air is taken as 1000:1.The location of the free surface, which means air and water interface, is determined as the density contour with average value of the two fluids.The value of the gravitational acceleration is taken as g = 9.8 ms −2 .The pseudo time step is set as τ = 5 × 1.0 −3 .The value of the artificial compressibility parameter is used by β = 500.

Extreme wave generation in the tank
In this work, we adopt the wave tank geometry and set up conditions used in the experimental model described by Ning et al. (2009).An experimental tank was used with the dimensions 69m × 3.0m and water depth was set to0.5m.In the study Ning et al. (2009) and Westphalen et al. (2008), four extreme wave cases are investigated with different input amplitudes.Here we reproduce and validate numerically Case 3 only.
In our NWT, the wave characteristics in each case are as shown in Table 1.The length of the numerical tank is the same as the one used by Ning et al. (2009), which is 5 times the characteristic wave length (5λ).In the vertical the tank is set to 1.0 m and the water depth h = 0.5 m is the same as in the physical experiments.The width of the tank is taken with a 3 cell layer thickness as 3D, albeit in a narrow numerical tank.In Ning et al. (2009), the parameters used were the focus point x 0 = 3.27 m and the focus time t 0 = 10.0 s, respectively.The wave maker signals for the simulations were calculated using a first and first plus second order formulation and 16 (=N) wave components as recommended by Ning et al. (2009).The corresponding  wave spectrum can also be obtained as shown in Fig. 1 by a Fast Fourier Transform.

8
Firstly results obtained with the inclusion of first order and first and second order wave components for comparison with the experimental data of Ning et al. (2009) are presented in Fig. 4, which shows the surface elevation at the focus point.The total number of cells in the 3D domain is 471,900 with a uniform mesh spacing of 0.01667m.The results confirm that agreement with the experiments up to the time of wave focussing is satisfactory particularly in the first and second order case, with any discrepancies between experiments and simulations appearing only after the focal time.
A non-uniform mesh was used for the case with 358×70×3 = 75,180 cells and mesh spacing in the refined regions of 0.00875m, which gives 10 cells per wave amplitude in the vertical direction.The results with the inclusion of first order and first and second order simulations for comparison with the experimental data of Ning et al. (2009) are shown in Fig. 5, which illustrates the surface elevation at the focus point.A comparison of the maximum surface elevation between simulation and experiments after the focal point (including the focal point) are shown in Table 2.It can be seen that the first order wave maker signal underestimates the velocities whilst the results for the first plus second order case provide reasonable agreement with a fully nonlinear calculation and by implication with the experimental results.The total CPU time is about 45 hours on one processor of a 600 MHz NEC vector computer.

Wave period T(s)
Frequency Band (Hz) Case 3 1.25 2.18 0.0875 0.6-1.4        2. It can be seen that the first order wave maker signal underestimates the velocities whilst the results for the first plus second order case provide reasonable agreement with a fully nonlinear calculation and by implication with the experimental results.The total CPU time is about 45 h on one processor of a 600 MHz NEC vector computer.

A fixed horizontal cylinder in regular waves
The case next considered is the interaction between regular waves and a half submerged horizontal cylinder in a tank.The purpose of the test case is to provide validation of the wave forces acting on the cylinder compared with the theory based on Morison's equation and experimental results (see Dixon et al., 1979;Easson et al., 1985).To correspond with the physical experiments, first-order regular waves are generated in a tank to interact with the cylinder.The inflow boundary velocity components u, w and the surface elevation  specified using non-reflecting boundary conditions allowing air to leave or enter the domain.The remaining boundaries are set as rigid walls.
The NWT geometry used had outer dimensions 12 m × 1.5 m × 0.21 m and the water depth used was h = 1.0 m.The position of the cylinder was set about one wave period (wave length of 3.90 m) from the wave maker (see Figs. 6 and 7).A non-uniform mesh was used for the case with 475 × 69 × 14 = 458850 cells and mesh spacing in the refined regions = 0.015 m as shown in Fig. 8.The set up parameters are: cylinder diameter, D = 0.25 m, length of cylinder l = 0.12 m, wave amplitude A = 0.125 m, wave frequency ω = 3.817 and k = 1.61.Figure 9 shows time histories of the relative vertical force over one period, which relative force defines to F = F z / gρ πD 2 l/4 and F z is a force on the cylinder resulting from the pressure acting on the surface calculated in the vertical direction.It can be seen that good agreement is achieved with the experimental data and theoretical forces providing satisfactory evidence of the accuracy of the present model.

Wave interaction with the floating Bobber in the tank
The NWT domain for this extreme waves case is 13 m × 1.0 m × 0.48 m with a water depth for the tests of h = 0.5 m.The geometry of the Bobber is as follows: the diameter of the hemispherical base is 0.  The mass of the Bobber geometry is taken as t Bobber geometry is allowed to articulate in he excitation, while all other modes are restrained.A vertical velocity, displacement and heave disp corresponding results for the free surface elevatio is shown in Fig. 12; the heave force on the Bobbe velocity in Fig. 14, and the heave displacement o Fig. 16 illustrates the wave profile around the B incoming wave, diffracted wave, radiated wave a the body. -0. -0.
-0. 0. 0.  .The geometry of the Bobber is as follows: the diameter of the hemispherical base is 0.3m, the vertical sides extend to a flat top 0.15m above the curved section.The initial position of the apex of the Bobber geometry in the tank is at the wave focus point 3.0m×0.24m×0.35mand a non-uniform mesh is used with 425×40×22 = 374,000 cells and refined regions with local mesh spacing of 0.02m around the geometry (see Figs. 10 & 11).The focus time, focus point and the input wave amplitude are set up in the same manner as before in the empty tank.Reflection boundary conditions are used on the Bobber boundary and the other boundaries are specified in the same manner as in the case of the empty NWT.

Relative vertical force
The mass of the Bobber geometry is taken as the volume of the hemispherical base.The Bobber geometry is allowed to articulate in heave motion only, responding to the wave excitation, while all other modes are restrained.As expected, the response is in terms of the vertical velocity, displacement and heave displacement of the Bobber geometry.The corresponding results for the free surface elevation on the front side of the Bobber geometry is shown in Fig. 12; the heave force on the Bobber geometry is shown in Fig. 13; the vertical velocity in Fig. 14, and the heave displacement of the Bobber geometry is shown in Fig. 15.Fig. 16 illustrates the wave profile around the Bobber geometry.These results include the incoming wave, diffracted wave, radiated wave and the wave created by the heave motion of the body.The NWT domain for this extreme waves cas the tests of m h 5 .0 = . The geometry of th hemispherical base is 0.3m, the vertical side section.The initial position of the apex of t focus point 3.0m×0.24m×0.35mand a non-un cells and refined regions with local mesh spa 10 & 11).The focus time, focus point and th manner as before in the empty tank.Reflecti boundary and the other boundaries are spec empty NWT.
The mass of the Bobber geometry is taken Bobber geometry is allowed to articulate in excitation, while all other modes are restrain vertical velocity, displacement and heave corresponding results for the free surface elev is shown in Fig. 12; the heave force on the Bo velocity in Fig. 14, and the heave displaceme Fig. 16 illustrates the wave profile around t incoming wave, diffracted wave, radiated wa the body.The mass of the Bobber geometry is taken as the volume of the hemispherical base.The Bobber geometry is allowed to articulate in heave motion only, responding to the wave excitation, while all other modes are restrained.As expected, the response is in terms of the vertical velocity, displacement and heave displacement of the Bobber geometry.The corresponding results for the free surface elevation on the front side of the Bobber geometry is shown in Fig. 12      heave force on the Bobber geometry is shown in Fig. 13; the vertical velocity in Fig. 14, and the heave displacement of the Bobber geometry is shown in Fig. 15. Figure 16 illustrates the wave profile around the Bobber geometry.These results include the incoming wave, diffracted wave, radiated wave and the wave created by the heave motion of the body.

Conclusions
A numerical tank with free surface capturing and Cartesian cut cell method has been modified and developed for the simulation of wave energy devices under modelled extreme wave conditions.The results show the emerging promise of the NWT for the simulation of nonlinear wave interactions with fixed and floating bodies.Future work will include extensions to other wave energy converter devices admitting a full complement of degrees of freedom as opposed to heaving alone as considered here, and to wave interactions with other floating bodies and fixed structures including aeration of the water component under violent wave impact situations.

Fig. 1
Fig.1 Input wave spectrum with frequency Fig.2 Surface elevation at input position
kh) η = Acos(kx − ωt) with the velocity specification applied in the water component only and the velocity of the air at the inlet boundary set to zero.The top boundary and right far boundary are Nat.Hazards Earth Syst.Sci., 11, 519-527, 2011 www.nat-hazards-earth-syst-sci.net/11/519/2011/ 10 pressure acting on the surface calculated in the vertical direction.It can be se agreement is achieved with the experimental data and theoretical force satisfactory evidence of the accuracy of the present model.

Fig. 6 Fig. 6 .
Fig.6 A horizontal cylinder Fig.7 A horizontal cylinder in the N Fig. 6.A horizontal cylinder.
3 m, the vertical sides extend to a flat top 0.15 m above the curved section.The initial position of the apex of the Bobber geometry in the tank is at the wave focus point 3.0 m × 0.24 m × 0.35 m and a nonuniform mesh is used with 425 × 40 × 22 = 374000 cells and refined regions with local mesh spacing of 0.02 m around the

Fig. 8
Fig.8 Cartesian cut cell mesh around Fig horizontal cylinder

Fig. 8
Fig.8 Cartesian cut cell mesh around Fig.9 Relative vertical forces on horizontal horizontal cylinder cylinder
Fig. 11.The Bobber in the NWT.

Fig. 12 .
Fig. 12.Time history of wave run-up on the front side of the Bobber geometry.

Fig. 13 .Fig
Fig. 13.Time history of the heave force on the Bobber geometry.

Fig. 14 .Fig. 14
Fig. 14.Time history of the vertical velocity of the Bobber geometry.

Fig. 12 Fig. 14 Fig. 16 .
Fig.12 Time history of wave run-up on the Fig.13 Time history of the heave f front side of the Bobber geometry on the Bobber geometry uvn y − uwn z 2un x +vn y +wn z

Table 2 .
It can be seen that underestimates the velocities whilst the results for the reasonable agreement with a fully nonlinear calcula experimental results.The total CPU time is about 45 ho NEC vector computer.

Table 1
Characteristic waves

Table 2 .
Maximum surface elevation at wave gauges.