Experimental and numerical investigation of wave resonance in moonpools at low forward speed

Applied Ocean Research 47 (2014) 28-46

ELSEVIER

Contents lists available at ScienceDirect

Applied Ocean Research

journal homepage: www.elsevier.com/locate/apor

O C B A N

RESEARCH

Experimental and numerical investigation of wave resonance in

moonpools at low forward speed

Arnt G. Fredriksen*, Trygve Kristiansen, Odd M. Faltinsen

Department of Marine Technology, Centre for Ships and Ocean Structures, Norwegian University of Science and Technology, Trondheim, Norway

CrossMark

A R T I C L E I N F O Article history:

Received 10 October 2013

Received in revised form 27 February 2014 Accepted 19 March 2014

Available online 19 April 2014 Keywords: Gap resonance 2D experiments Moonpool Numerical wavetank Domain-decomposition

Harmonic polynomial cell method Potential flow

Viscous flow Flow separation Finite volume method

A B S T R A C T

In o r d e r to s t u d y the b e h a v i o r of r e s o n a n t p i s t o n - m o d e r e s o n a n c e in a m o o n p o o l at l o w f o w a r d / i n c o m i n g c u r r e n t speed, w e p e r f o r m e d a series of e x p e r i m e n t s a n d c o m p a r e d t h e m to a n o n l i n e a r h y b r i d m e t h o d w h i c h couples p o t e n t i a l a n d v i s c o u s flow. T h e s e t t i n g is a 2 D box s e c t i o n w i t h a m o o n p o o l gap i n t h e m i d d l e , forced to oscillate i n h e a v e w i t h a g i v e n a m p l i t u d e and f r e q u e n c y , w h i l e s i m u l t a n e o u s l y t r a v e l l i n g at a g i v e n c o n s t a n t f o r w a r d s p e e d . T h e n u m e r i c a l m e t h o d couples a N a v i e r - S t o k e s ( C F D ) s o l v e r u s i n g the Finite V o l u m e M e t h o d ( F V M ) , w i t h a potential f l o w m e t h o d u s i n g the H a r m o n i c P o l y n o m i a l Cell m e t h o d ( H P C ) . It is f o u n d that the m o o n p o o l b e h a v i o r is slightly r e d u c e d w i t h a l o w f o r w a r d velocity, and the r e d u c t i o n is d e p e n d e n t on the h e a v e f o r c i n g a m p l i t u d e .

© 2 0 1 4 E l s e v i e r Ltd. A l l rights r e s e r v e d .

1. Introduction

Marine operations f r o m ships often involve moonpools to lower or l i f t devices such as subsea modules and ROVs. Resonant piston-mode resonance can be excited by the relative vertical ship motions in the neighborhood of the moonpool and cause strong amplifi-cation of the dynamic wave elevation i n the moonpool. The fact that the resonant piston-mode frequency is typically in the vicin-ity of the heave natural frequency of the ship limits possible sea states for a marine operation. The stronger the shed vorticity due to f l o w separation at the moonpool entrance and inside the moon-pool is, the larger the damping is, and the smaller the maximum resonant piston-mode wave amplitude is for a given ship and main moonpool geometrical parameters i n a given sea condition. It is of practical interest to know the free-surface elevation in the moonpool and the ambient f l o w velocities and accelerations in the vicinity of the moonpool i n order to assess the loads on lifted or lowered devices through the moonpool.

Several authors have studied the importance of including vis-cosity and f l o w separation in the resonant moonpool problem.

* Corresponding author. Tel.: +47 97504522.

£-mcii7 addresses: arnt.g.fredriksen@ntnu.no, arntgunv@gmail.com (A.G. Fredriksen).

0141-1187/$ - see front matter © 2014 Elsevier Ltd. All rights reserved, http://dx.doi.0rg/10.1016/j.apor.2014.03.005

due to the fact that a potential f l o w solution w i l l greatly over-predict the piston-mode amplitude at resonance for sharp-edged lower entrances of the moonpool. Faltinsen et al. [5] investi-gated forced heave of a two-dimensional moonpool section using a domain-decomposition (DD) scheme w i t h i n the framework of linear potential flow theory. Their DD scheme led to a system of integral equations on the transmission interfaces that solved for the piston-mode natural frequency and the steady-state piston-mode amplitude. To improve the potential f l o w models some authors have tried to fit an artificial, empirically based damping to the free-surface condition inside the moonpool. This is known as a numerical damping l i d . Lu et al. [15] investigated the possibil-ity of finding this damping coefficient based on experimental and CFD results. The damping coefficient was i n their work observed not to be sensitive to the variation of moonpool gap width, body draft, breadth-to-draft ratio and body number. Their focus was on wave forces, where Lu et al. [16] used the same setup w i t h focus on the wave elevation i n tlie moonpool. It is not known when using a numerical damping l i d how well the flow i n the vicinity of the moonpool is predicted. Lu and Chen [14] inves-tigated what contributed to the dissipation of the piston-mode amplitude generated f r o m incoming waves. Both the dissipation f r o m the boundary layers inside the moonpool gap, and i n which fluid areas around the moonpool gap the vorticity dissipation was largest.

KG. Fredriksen et al. /Applied Ocean Research 47(2014) 28-46 29

Fig. 1 . Domain overview witti the HPC domain close to the free surface {Qpot) and the FVM domain around the ship edges and down to the bottom of the tanl< ( Q C F D ) . Notice how the grid in the HPC domain follows the free surface.

Kristiansen and Faltinsen [10,11] showed how a hybrid method based on a linear potential solver coupled to a viscous Navier-Stokes (CFD) solver could be used to predict the piston-mode in a 2D moonpool under forced heave oscillations. The viscous domain covered the inlet of the moonpool where the vortic-ity would be shed, and later advected before dissipating, The rest of the domain, including the free surface was discretized using linear potential f l o w theory, where it was solved for the linear acceler-ation potential A direct coupling between ir i n the potential domain and the pressure p i n the viscous domain was used. This guaranteed continuity in pressure and normal velocity across the intersection between the two domains. Both domains where solved using a Finite Volume Method (FVM) w i t h second order accuracy in space. Their thought is that potential flow is best at propagat-ing waves, and that the viscous domain can incorporate vorticity separated f r o m the edges. The method is extremely fast relative to solving the complete viscous f l o w problem w i t h nonlinear free-surface conditions by means of a CFD solver w i t h , for instance, the Volume of Fluid (VOF) method for capturing the free surface. Later, Kristiansen et al. [12] demonstrated that the method gave good results also for the 3D moonpool problem, w i t h modest computa-tional times.

In the previously described works, there was no current or forward velocity. Fredriksen et al. [7] developed the method introduced by Kristiansen and Faltinsen [11] i n order to account for this. A linear perturbation of the free surface is adequate to cap-ture the moonpool f l o w at zero forward velocity and no incoming current. However, a linear perturbation method is insufficient to capture the wave-current interaction. Here in our hybrid coupling method a higher order perturbation method is unsuitable since a single unknown (pressure p i n the viscous domain) can not be cou-pled to several unknown velocity potentials i n the potential flow domain. Fredriksen et al. [7] used nonlinear free-surface conditions in a non-rotating body-fixed coordinate system to overcome this problem. The body-boundary conditions are then exactly satisfied, but a new complication w i t h the free surface is introduced. The matrix system now needs to be updated each time-step to satisfy the nonlinear free-surface conditions on its exact position. The cur-rent work is a further development of this nonlinear hybrid method, where we use the recently developed Harmonic Polynomial Cell (HPC) method in the potential domain instead of the FVM.

The HPC method was first presented for solutions of the 2D Laplace equation in [19], where both computational speed and accuracy were compared against boundary element methods and other field solver methods. It was later verified i n 2D w i t h wave propagation over a submerged trapezoidal bar, see [20]. The first extension to 3D was done by Shao and Faltinsen [21 ], where again, its high efficiency and accuracy was compared w i t h other solvers. Forces on a bottom-mounted free surface piercing vertical circu-lar cylinder i n nonlinear regucircu-lar waves were studied. Nonlinear forces describing up to the fourth harmonic was computed and compared against other numerical results and experimental values with satisfactory results.

There exist several other strategies for coupling viscous f l o w and potential f l o w models. It can simply be done by using a potential flow model to generate initial conditions to a viscous flow model.

An example of this is by using a potential flow model to simulate a wave breaking up to when the free surface intersects itself, then use the potential flow results to generate initial conditions to a viscous flow simulation, see [9]. A stronger coupling strategy w h i c h is similar to ours is summarized in Ref. [8]. Basically they solve a potential flow problem on a large domain using the boundary element method (BEM). On a smaller viscous CFD domain, the NS-equation is split i n an inviscid and a viscous part, u = u'+u^ and p = p ' + p ^ . Since the inviscid part is known from the potential flow BEM calculation, the NS-equation can be solved for and p^. They use this strategy to solve a sediment transport model, where the viscous CFD domain is located close to the sea bottom.

We start by presenting the numerical method i n Section 2, then give an overview of the experimental setup in Section 3. The results are given and discussed i n Section 4. At the end a conclusion about the present w o r k w i l l is given in Section 5.

2. Numerical approach

We present a coupling strategy between a potential flow outer domain w i t h a viscous flow inner domain. Only the water domain is considered (see Fig. 1). The chosen method for the potential flow domain is based on the HPC method and the flow i n the inner viscous domain is solved based on a laminar flow assumption using FVM. The governing equations w i l l be solved i n a body-fixed non-rotating coordinate system, w i t h exact boundary conditions. The numerical method is then capable of simulating free-surface flows w i t h flow separation f r o m free-surface piercing structures. Only flow separadon f r o m sharp corners w i l l be considered w h i c h implies that the details of the boundary-layer flow are not needed. We w i l l here consider the case w i t h forced heave oscillations of a 2D moonpool w i t h l o w forward velocity. The Earth-fixed coordi-nate system is here defined fixed in the initial position of the body, w i t h z = 0 on the mean free surface and the z-axis positive upwards. Correspondingly the body-fixed coordinate system is following the motion of the body in heave and sway.

We w i l l solve the Laplace equation for the absolute velocity potential tp i n the potential flow domain, which is valid for i r r o -tational flow of an incompressible and inviscid fluid, i.e.

V^(p = 0 i n Spot. (1)

The absolute velocity is defined as u = V ^ . We w i l l solve for the absolute velocity potential ^ in a body-fixed and non-rotating coor-dinate system. In the CFD domain we w i l l solve for the relative velocity Ur = ( f r , Wr) in the same body-fixed non-rotating coordi-nate system. The governing equations for mass and momentum conservation i n an incompressible, laminar viscous fluid flow in an accelerated non-rotating coordinate system are given according to Ref. [4] as

V-Ur = 0 i n f2cFD (2)

^ 5 ^ + U r - V U r = - - V p - g k + V . V V - a o in ^2cFD, (3)

dt p

where d'^Urjdt means the time-differendation for a fixed point i n the body-fixed coordinate system, (dvr/drlj + (dwr/dt)k, i.e. we do

30 A.C. Fredriksen et al./Applied Ocean Reseairh 47 (2014)28-46 not time differentiate thie ttnit vectors. However, ttie latter fact

does only matter for a rotating coordinate system. Here ao is the forced body acceleration i n sway and heave. We w i l l introduce P = (P/p) + gz in order to simpHfy the equations. The procedure does not account for the fact that vorticity flow outside boundary layers becomes turbulent for relatively small Reynolds numbers ([18] or [22]). To solve the Navier-Stokes equation numerically the explicit fourth-order Runge-Kutta scheme is used. On each sub-dme-step we apply the well known Chorin's projection method, see [1 ]. One of its key features is that it decouples the computations of the velocity and pressure fields. The method consists of three steps. Step 1 involves advection w i t h o u t the pressure and viscosity terms, resulting i n the temporary artificial velocity field u*. Step 2 involves viscous diffusion terms resulting in the temporary artificial velocity field u" and step 3 solves the pressure and updates to a divergence-free velocity field u'^"^^ The advection step is i n the present code solved using a simple linear upwind method. According to Ref. [11] this is good enough for the zero velocity case, where shed vordcity is important only for about half a period. Srill this does not guaran-tee a good solution when the vorticity separated from the leading edge influences the vorticity separated f r o m the moonpool edges. By leading edge we refer to the bilge in the upstream direction of the body (see point A i n Fig. 1).

2. J. Harmonic polynomial cell (HPC) metYiod

Details i n this section are mostly from [20], but given here to make the details in the coupling at the intersection between poten-tial and viscous flow domains clearer. It is complicated i n a FVIVI to implement a scheme higher than second order spatial accuracy. The main motivation for introducing the HPC is the high spatial accu-racy, which also implies that larger cells can be used compared to a FVIVI. The free-surface waves are propagated w i t h a higher accuracy using the HPC scheme. In particular the wave celerity is captured more accurately than w i t h a conventional FVM. Another feature that makes the HPC method attractive is that the harmonic polyno-mials automatically satisfy the Laplace equation everywhere, then we only have to make sure that the multiplying factors of each poly-nomial are consistent w i t h the boundary conditions both globally and locally.

When describing the HPC method we w i l l operate w i t h a local Cartesian coordinate system. The harmonic polynomials are i n t w o dimensions given by the real and imaginary parts of the complex polynomial

z" = ( x - F i y ) " , (4)

where n is the order of the polynomial and i = In this method we include all harmonic polynomials up to 3rd order and the 4 t h order harmonic polynomial corresponding to real part of Eq. (4). Then we can write 95 as a linear combination of 8 different polyno-mials, where each polynomial satisfy the Laplace equation, i.e.

9'(y.z)

&

9

Fig. 2. An overview of tlie local HPC numbering. Note that this grid will overlap with neighbouring grids, i.e. node 9 in the above grid can be any of the other 8 nodes in the neighbouring cells.

for a definition on the local HPC numbering.) This gives a linear relationship between the coefficients bj and the values tpi,

<Pi (6)

where the element i n d / j = ƒ)• (y,-, z,-) f r o m (5) defines the matrix [D], Assuming 8 points where the value of cp are known, the coefficients

bi are found as,

(7)

Here the matrix [C] is defined by the elements c,j and is the inverse of the matrix [D]. Furthe/, given a ninth point i n the middle of the 8 points (see Fig. 2), we have a basis for constructing a polynomial valid at the middle point. This means

^ö(y,z) = ^ 1=1 _j=i

(8)

Choosing the middle ninth point to be yg = 0 and Zg = 0 the above equation simplifies to only consist of the first constant polynomial. The reason is t h a t / i = 1 a n d ^ = 0 w h e r e j = 2 . . .8, such that

9^9 ( y 9 , Z 9 ) = ^ C , , | ^ ( . ₍₉₎

This implies that ^g can be w r i t t e n as a linear combination of the 8 neighbouring values of (p.

Since w e easily can find analytical expressions for the deriva-tives 9/9y and 9/9z of a polynomial, we can easily find expressions for dcpjdy and 9^/9z,

(5) = bi + bjy + Ö3Z + Ö4 (y2 - z2) + b^yz + be (y^ - 3yz2) + by {3yh - z^) + bg ()/4 _ syZzZ +

Note that z is no longer a complex number, but a coordinate in . the z-direction. The y-axis is the horizontal axis and the z-axis is

the vertical axis, positive upwards. Eq. (5) leads to a linear sys- 9,^9 ^

tem of equations, where we find ^ at 8 different locations <p = (pi, " g ^ = 0,zg = 0) = 2_^C2j(Pi (10)

A.G. Fredriksen et al. /Applied Ocean Research 47 (2014) 28-46 31

^ ( y 9 = 0 , Z 9 = 0) = ^ C 3 , , - ^ , . . ₍₁₁₎

Determining tlie coefficients i n the matrix [C] can be solved as a sub Dirichlet problem for each node. The 8 neighbouring points will span the boundary edge on the local polynomial for the ninth point, and these points w i l l define the elements djj of matrix [D]. The wanted r o w in the matrix [C] can now be found by inverting the matrix [D]. Doing this for all nodes throughout the fluid w i l l lead to a sparse matrix system w i t h at most 9 non-zeros on each row. For global boundary nodes we can choose a different local numbering than i n Fig. 2 and use any of the boundary nodes to be point 9.

2.2. Finite volume method

Details in this section are mostly f r o m [6] and [4], but given here to show the details of the coupling at the intersection. Rectangu-lar cells and a staggered grid arrangement scheme for the velocity and pressure nodes are used. More precisely we put the pressure nodes in the center of the cells, the horizontal y-velocity nodes on the midpoint on the vertical cell edge and the vertical w-velocity nodes are on the midpoints on the horizontal cell edges. This means that the control volumes for the pressure, f-velocity and w-velocity nodes are different. The above described HPC method is coupled to the solution of the pressure when solving the Navier-Stokes equa-tion. Using Chorin's projection method to solve the Navier-Stokes Eq. (3) implies that we solve a Poisson equation for the pressure.

2piM+l _ 1 _

A t V . u * * (12)

where M is the time-step number. To simplify the notation we now use that p = p in the rest of the text. The solution of the pressure is found by using FVM, i.e. we start out by integrating the equation over small cells.

v 2 p d A = / V.u**dA. A JA

(13)

By using the divergence theorem, we can write the discretized Pois-son equation as

i=l ^ I i=1

ds,-. (14)

where the sum is over the four cell edges, and the normal derivative can be further discretized using a finite difference approximation using the two closest nodes. This limits the accuracy on the veloc-ities to first order i n space. To easier explain h o w we transfer information between the FVM domain and the HPC domain, we write out the above equation. Consider a cell w i t h a pressure node P P in the center, it has four neighbour cells i n 2 dimensions, where P N is the pressure node at the midpoint of the north neighbour cell a distance A z f r o m pp, p£ is the pressure node at tlie midpoint of the east neighbour cell a distance A y f r o m pp and similarly for ps (south) and pw (west). On the midpoints of the four edges of the

P E - P ? Ay Az-P Az-P + Ay

w

cell the velocity nodes are located. The horizontal velocities Ve on the east and Vw on the west edge, and the vertical velocities Vn on the north and Vs on the south edge.

PE - PP. Az _ pp

7 Pw ^

A y Ay Az

PP-PS Az 1

= {v*e*Az - v*„'Az + w**Ay -~ vvf A y ) (15)

The equation above is w r i t t e n w i t h constant grid values for clarity, however in the implementation i n the code we can have a varying grid.

2.3. Intersection

On the intersection we use the Bernoulli equation to connect the pressure i n the potential HPC domain to the pressure in the viscous FVM domain. Here Bernoulli's equation is given for a fixed point in the body-fixed non-rotating coordinate system f r o m [4],

nM 1 , . . . I

^M+1 _ ^ a T

^M+1 ₍₁₆₎

where uo = [)72,'?3] is the forced velocity of the body. The corresponding term U Q • V ^ represents the transition from an Earth-fixed coordinate system to a body-fixed coordinate system. The finite-difference approximation i n time used for the time derivative of the velocity potential is more accurate in the middle of two time-steps. Therefore the pressure i n the potential domain and the viscous domain w i l l always be half a time-step "out of phase" w i t h each other.

Let us use the two above Eqs. (15) and (16) to explain w i t h a simple example where the potential domain is on the west (left) side, and the viscous CFD domain is on the east (right) side. On the potential side we solve for q>w, and not pw- The second term we need to consider is how we update the temporary velocity node

u^. We use the following approach: The velocity on the next

time-step (M+1) should be independent on which of the t w o methods are used to obtain the velocity. From potential theory the velocity is updated as u f + ^ = V(/^+^ - U Q , while f r o m the viscous side i t is uj^+^ = u** - A t V p ' ^ + i . Thus, we find the artificial velocity on the intersection as

U** = Ve 1+1 _{- Uo -I- A f V p} M+l ₍₁₇₎ and again we can use the time-discretized Bernoulli Eq. (16) to change variables

U** = V ^ " - A t V Q l V ^ j ' ^ l ^ .- -- U o - V < - U o . (18) Note that all dependency on the next time-step M + l are gone, and we are left w i t h only values f r o m the previous time-step. Instead two higher-order terms V ^ and V (uo • Vcp'^^ appear. How-ever, due to the higher order spacial accuracy of the HPC method we are able to estimate them w i t h high accuracy. On the velocity nodes on the intersection we need to evaluate the following terms.

V V(p = 2 dtp d^cp dtp d^tp dtp d'^tp dtp d^cp'^

dy 9y2 9z dydz' dy dydz dz dz^ (19)

and the other higher order term w i l l look very similar. The resulting equation for our example then becomes.

dz dydz (20) A z A y At 1 \V<Pw\ --ao-V<pv.

Here all values on the right hand side of the equation are based on the previous time-step, while the left hand side is the unknowns to be solved for on the next timestep. Note that we get and A f -terms i n the matrix system, w i t h a'resulting change i n the matrix

32 _{A.G. Fredriksen et al./Applied Ocean Research 47 (2014)28-46} V _{1 1} V> •p V X ¥^ X i - ^< X V V X ^ _{^ 0 ji o} P ' V P o V ^ X , V p O p O N / P 0 p 0 X i V p O f i P 0 O p O p 0

Fig. 3 . Overview of wliere ttie different nodes are located. <p HPC nodes are on the corners of a cell, viscous p nodes are on the center. The velocity nodes are located on the edge (.), the [/-velocity node are on the west/east edge and the w-veloclty node are on the north/south edge, x symbols shows where values need to be interpolated.

system between the sub time-steps in the explicit fourth order Runge-Kutta method.

Since the HPC domain and the FVM domain do not share a com-mon node, we solve this by interpolating to a value of tp on the FVM node i n the potential domain. The problem is that cp is not uniquely defined for any given point outside the HPC nodes, i.e. tp can be found f r o m four different harmonic polynomials. The reason is that four HPC cells (see Fig. 3) overlap each FVM node. Instead of eval-uating the average of four different harmonic polynomials, which would involve in total 16 values of tp, a linear weight function is applied, such that the HPC node closest to the FVM node gets the most influence.

Vw = Wnwtpwnw + Wne(pwne + WseCPWse + Wsw<PWsw (21) where Wnw is the weight applied to the north west HPC neighbour node tpwnw of the FVM node tpw- A similar procedure is followed for the north east tpwnw, south east tpwse and south west tpwsw neighbouring HPC nodes. Here the four weight values should sum up to 1. On the potential f l o w side we need to similarly change vari-ables f r o m ^ to p by using Bernoulli's equation (16). The equation for the HPC cell is given from Eq. (9),

C l , l < O l + C\,2<P2 + C\,3<P3 + C\,4(P4 + Cl,5(p5 + C-i^etpe

+ Cyj(p7 + C]^S(P8-<P9 = 0 (22)

Again we use a simple example to illustrate how the equation sys-tem is changed at the intersection. Now the HPC cells have the node 8 (lower right) inside the viscous domain,

Cl,l(P\ +Ci,2<P2 + C\3(ps + C i , 4 ^ 4 4-Ci,5^£)5 -|-Ci,s(p6

+ C l , 7 « 0 7 - C i , 8 A f p 8 - ( p 9 ₍₂₃₎

= Cl,8 ; - A t Uo • V^ög

Also here we have to interpolate existing values of the pressure to a value at the HPC node pg. Similarly as Eq. (21) we use a weighted average between the four closest FVM nodes.

P8 = WnwPsnw + VV„ePsne + WsePSse + WswPssw (24)

After the coupled matrix system is solved, we know cp on all HPC nodes i n the potential domain andp on all FVM nodes i n the viscous domain. We update on each side of the intersection the nodes to the other value, such that both <p and p are known at the next time-step i n an area close to the intersection. This means that for all practical purposes when running the code the potential f l o w area are extended into the CFD domain, and potential f l o w theory should be valid here. Since we are aiming at solving a resonance problem,

it has been important to get the phase of all terms consistent i n the explicit fourth-order Runge-Kutta time-integration framework.

2.4. Free-surface conditions

We start out w i t h the nonlinear free-surface kinematic and dynamic boundary conditions i n an (non-accelerating) inertial Earth-fixed coordinate system. Our goal is to find the corresponding semi-Lagrangian nonlinear free-surface conditions in a non-inertial non-rotating body-fixed coordinate system. To find the kinematic free-surface condition we use that a fluid particle on the free sur-face remains on the free sursur-face. Then the kinematic free-sursur-face condition is mathematically given in an Earth-fixed coordinate sys-tem as the material derivative of a function Z=z - f ( x , y , t) = 0, see [4] D ( ? - z ) 9? Df ,^dvdl_d(p^^ dt dy dy dz (25) (26) Change of time derivative f r o m an inertial coordinate system to body-fixed (non-inertial) coordinate system is represented by

Here the derivative d''/dt represents change of a value i n time seen f r o m a point fixed in the body-fixed coordinate system, as earlier defined in Eq. (3). The kinematic free-surface condition i n a body-fixed coordinate system when combining Eqs. (25) and (26), and noting that we should exchange cp w i t h the function Z i n Eq. (26) then becomes

d ^ _ 9^0

dt ~ dz

9 ^ 9 f . 9 f

dy dy dy •>13 (27)

On the free surface we require continuity in pressure across the interface, such that the pressure i n the liquid at the free surface is equal to the pressure po i n the air. First Bernoulli's equation is w r i t t e n i n the Earth-fixed coordinate system as

p' dtp 1 „ 2

^ + ^ + 2 l H +SZ^

where p' is the pressure and C(t) is a time-varying constant. Next we evaluate the equation at z = f where p' = Pa = 0, then C(f) = 0 can be determined. The result is combined w i t h Eq. (26) to transform Bernoulli's equation to the dynamic free-surface condition i n a body-fixed coordinate system, , ••

d''<p „ 2 . d t p . d t p

(29) where ^f^/is the wave elevation seen in the Earth-fixed coordinate system. We reach the final expression for the dynamic free-surface boundary condition by rewriting Bernoulli's equation to express the time rate of change of tp on the free surface, as one travels w i t h the wave vertically. Similarly as Eq. (26) we can set up a relation-ship between the time derivative when following a fluid particle in vertical direction and the time derivative of a point fixed i n a body-fixed coordinate system.

d * ^ ^ d V

dt ~ dt ^ dt d z ~ dt ^ dt dz

Inserting for Eqs. (29) and (27) gives

(30)

d*(p

dt

9 ^ 9 ^ 9 f . fdtp d^d^

A.C. Fredriksen et al. / Applied Ocean Research 47(2014) 28-46 33 Here we define tlie derivative d'/dt as tiie cliange of a value in time

wlien moving vertically w i t h the free surface, while 9/9t represents the change of a value in time when fixed to a point in Earth-fixed coordinate system. Here we have used that é'zl<lt=d^l;jdt.

The derivatives of tp are found by evaluating Eqs. (10) and (11). For the y-deri vati ves of the free-surface elevation 9f/9y, we locally construct a polynomial up to fourth order and take derivatives of this.

The free-surface conditions are stabilized by applying a 5 point Chebychev smoothing as originally described by [13]. We did not find it necessary to apply the smoothing each time-step for simu-lations without forward velocity, but found that 10-20 times each period was sufficient. However, smoothing is applied each time-step for simulations w i t h forward velocity, see the discussion about convergence w i t h regards to number of smoothing operations (Appendix A.4). Without a smoothing technique the simulations will break down for cases w i t h forward velocity, due to a growing saw-tooth instability both up- and down-stream. In addition to the normal saw-tooth instability reported in the literature, we here due to the coupling have an inconsistency in the modelling when vortic-ity reaches the intersection. This w i l l occur on the downstream side of the hulls, and normally not inside the moonpool gap. However, vordcity w i l l reach the potential f l o w domain also inside the moon-pool gap for higher heave amplitudes. The smoothing technique will help stabilize this instability, since this for most cases hap-pens downstream i t is thought not to influence the response of the piston-mode amplitude. The height of the potential f l o w domain is set before the computations start, and based on a best guess of the expected wave amplitude. Here we separate between the expected piston-mode amplitude and the expected radiated wave amplitude, such that the potential f l o w domain is smaller outside of the moonpool gap. A commonly used numerical beach known as the Orianski's condition [17], used by [2] among others, is applied here on the kinematic free-surface condition to damp the outgoing waves.

2.5. Body-boundary condition

On the body surface Sg we need to separate between the poten-tial SBP region and the viscous SBV region and therefore also two sets of body-boundary conditions are needed. In the viscous region we also need boundary conditions on the sub-steps in Chorin's projec-tion method. For the tangential velocity at the body we w i l l apply a no-slip condition i n the viscous domain, and a slip condition i n the potential domain. Due to the fact that we do not resolve the details of the boundary layer w i t h our mesh density we have not experienced any problems related to this inconsistency. The rea-sons w h y we do not need the details of the boundary layer are that flow separation occurs f r o m sharp corners and that viscous stresses are secondary. The condition on the normal relative velocity on the body is simply

Ur • n = 0 on S B (32)

We use this condition to generate a condition for the pressure gra-dient and the temporary velocity field u " . From the momentum Eq. (3), we find that the body-boundary conditions are.

dp

= - n ^ u o onbsv u** = - A t i i o onSflv

(33) (34) By remembering that we are solving for the absolute velocity potential in the potential domain we get the following boundary condition

dtp

= n • Uo on S B P

2.6. Bottom and outer wall boundary conditions

In the body-fixed coordinate system we w i l l observe that the bottom and outer walls are moving, which implies that we have a grid that is following the body i n both sway and heave. In order to satisfy the outer wall boundary conditions exactly, as w i t h the body-boundary conditions we would need to re-grid close to the outer walls. This approach is not followed, instead we use that the absolute velocity should be zero at the outer walls

Sw-Ur - n = - n - u o onSw (36)

w h i c h can be transferred into the following condition for the tem-porary velocity u " on the viscous part of the outer walls Swv<

u * * = u ' ^ - A t u o o n S w (37) Further, for the pressure gradient, we impose

= 0 onSwv

dn

Also here we need a separate condition for the absolute velocity potential tp on the potential part of the bottom and outer walls Swp, i.e. (38) 1^ = 0 onSwp dn 2.7. Linear method (39)

The "linear" method has four main differences f r o m the nonlinear method, (1) linear free-surface conditions, (2) linear body-boundary conditions, (3) larger viscous domain inside the moonpool gap and (4) solving for the acceleration potential which is in phase w i t h the pressure inside the viscous domain. Point (3) follows f r o m point (1) as a fixed grid can be applied and the poten-tial domain can then be minimized. Basically the same method as in [11 ] is used here, except that the potential domain is here solved using the HPC method. Because the Navier-Stokes equations are nonlinear, the "linear" computations are not completely linear.

3. Experimental approach

Our motivation was to continue on the experimental campaign started by [11], and perform parameter variations using an auto-matic control system for the job. The experiments were performed in a wave flume at the Marine Technology Centre at NTNU in Trond-heim. The wave flume is 12 m long, 0.6 m w i d e and w i t h a 1.0 m water depth. In both ends there where parabolic beaches w i t h upper position located just below (1 m m ) the free surface. Since the beaches occupied approximately 2.5 m each, the usable rail length for the carriage was 5.6 m. This fact together w i t h that we needed our experiments to reach a steady-state piston-mode oscillation amplitude limited the carriage velocity (U) we could use. We define a Froude number based on the total length of the model including the moonpool gap.

Fn =

Vs{2B + b)

(40)

(35)

where B is the breadth of one hull, and b = 0.18 m is the moonpool gap w i d t h . We have found that using a model of total length 0.9 m, limits the forward carriage velocity to Fn=0.08. The dimensions of each section where B= 0.36 m breadth, 0.585 m wide and w i t h variable draft. See Fig. 4 for a sketch of the experimental setup and Fig. 5 for a picture from the experiments.

The following five parameters were chosen to be varied dur-ing the test campaign: oscillation frequency, oscillation amplitude, draft, carriage velocity and moonpool edge profile. To be able to

AG. Fredril<sen et al. / Applied Ocean Reseai ch 47(2014) 28-46 34

^ r —

A

B T

\

Fig. 4. Principle sketch of model test setup, with a carriage on top of a glass covered wave flume with a parabolic beach at each end. Note that the drawing is not in scale.

perform this extensive test campaign an automated setup was designed and installed, such that experiments could be performed without anyone present. The only thing not automized was the cali-bration process, such that we stopped by at least once a day to check the calibration factor of the wave gauges. The automatic control sys-tem was set up to start the carriage from one side of the tank, then smoothly accelerate the carriage up to the wanted velocity while the heave oscillations also smoothly started. When it reached the other end, it smoothly slowed down and waited until the waves had died out. After waiting for 200 s it returned to the starting position w i t h the same heave amplitude, frequency and carriage velocity. In this manner, the same experiment was repeated twice, before a new heave amplitude, frequency or carriage velocity was tested. The waiting period o f 2 0 0 s was chosen by observation of the measured wave elevation, of the time needed for the waves to die out. However, the circulation i n the tank which was set up by the forward moving model might not have stopped. This means that when the model started moving i n the opposite direction it might be influenced by the residual global flow. We should also consider that the presence of the seiching mode of the tank w i l l result in a horizontal current at the middle of the tank, which also w i l l i n f l u -ence the relative velocity between the the ship and the water. A seiching amplitude of 1 m m w i l l give a horizontal current at the middle of tank o f 2.7% of the forward velocity at Froude number Fn = 0.04. The seiching amplitude for carriage velocities of Froude number Fn = 0.08 was found from the wave gauges to maximum be 0.5 mm, the wave gauges were then around 2 m away f r o m the end of the tank.

The model was equipped w i t h 4 wave gauges, two in the moon-pool gap and one on each side. The wave gauges i n the moonmoon-pool gap were located 6.0 cm f r o m the hull on each side. The wave gauge on the left was 26.5 cm f r o m the model side, and the wave gauge

Fig. 5. Picture of the experimental setup with no appendages, here the model is located close to left beach. To change to a model setup with appendages, only the gray part of the moonpool edges are changed with appendage profiles.

on the right was 21.0 cm f r o m the other model side. All four wave gauges were mounted on the rig and was then forced to move w i t h the forward velocity and heave oscillation. All amplitude results presented in this work are therefore f r o m a body-fixed point of view.

3.1. Experimental error sources

The rails and the glass walls were not perfectly aligned, such that the model had to be smaller than necessary. From trial and error we found that the model had to be 1.5 c m smaller than the w i d t h of the tank to minimize the contact between the model and the glass wall. Due to this gap we have a 3D effect i n the experiments that are difficult to estimate. Visually we could at higher velocities observe vorticity being shed f r o m the small gap between the glass and hull, into the moonpool gap, and behind the second hull. Since we still had some contact between the glass and the hull, the for-ward motion was influenced due to a varying friction i n the length direction of the tank. The control system managed to counteract and minimize this, but i t can not be neglected as an error source.

Another error source is the control system for the carriage and heave actuator. It was not able to reach the desired heave ampli-tude, but in average around 90% of the desired heave amplitude. For a test series w i t h a desired heave amplitude of 10.0 mm, the actual heave amplitude became 9.1 m m . It may not be characterized as an error source since we know the amplitude after the test, but since the control system did not perform as expected we mention it as an potential error source.

A third error source in the experiments is reflections from the beaches. The velocity of the carnage is much lower than the group velocity of the outgoing waves. Such that waves generate upstream may be reflected and influencing the result when the steady-state piston-mode amplitude is evaluated. Note that for both higher car-riage velocity and longer periods the steady-state piston-mode amplitude is found close to the beach at the end of the tank.

When accelerating the structure i n any direction a transient effect w i l l be generated. The transient effect w i l l generate outgoing waves, but w i l l also excite the natural piston-mode and the slosh-ing modes inside the moonpool gap. For the forward velocity case it means that the odd sloshing modes inside the moonpool may be excited, which mainly means that the flrst sloshing mode w i l l be excited at its resonance frequency. For the forced heave oscil-lation case without forward velocity, the transient start-up results in an excitation of the piston-mode at the piston-mode resonance frequency. The piston-mode w i l l decay due to wave radiation and viscous dissipation. The sloshing mode communication w i t h the water outside the moonpool w i l l be low, and not affected by flow separation at the moonpool edges. Since also boundary layer dissi-pation is small, the sloshing mode decay slowly. The first sloshing mode natural period is for the 18 cm long moonpool Ts = 0.48s.

The raw signal f r o m the wave gauges was band-pass filtered to remove frequencies above 1.9 Hz and below 0.5 Hz. This w i l l remove all higher harmonics, including the first natural sloshing mode. An inspection of the frequency spectrums computed by FFT of the raw wave gauge signals, shows negligible traces of nonlin-earities. We choose to tal<e the amplitudes f r o m the experiments as half of the distance f r o m the wave crests to the wave troughs over a steady, or near-steady interval i n time.

4. Results and discussion

For the zero current case w i t h forced heave we can consider our problem as symmetric about the mid-line of the moonpool gap, as long as the excitation amplitude is sufficiently small to avoid asym-metry in the shed vorticity pattern due to instabilities. The problem

A.C. Fredriksen et al. / Applied Ocean Research 47 (2014) 28-46 a ) b ) 35 \'

P

< t ] I I _

I [ ^ ^ I 1 I

a ) b ) c )

Fig. 7. Overview of edge geometries that were tested, (a) sharp edge profile, (b) appendage profile # 1 where each appendage dimensions was I S mm width by 9 m m height and (c) appendage profile # 2 where each appendage dimensions was 27 mm width by I S mm height.

is asymmetric w i t t i a low forward velocity of the hull. A f e w ques-tions and assumpques-tions then arise, w i t h the strength of the shed vorticity being a function of ttie local fluid velocity at any sharp edge. Let us for the moment assume that the shed vorticity is con-centrated i n thin free shear layers w i t h o u t diffusion. Further, we discuss the free shear layer separating from one corner and define the time rate of change of circulation dV/dt = iO.SU^, see [3]. V is the circulation around a closed curve C i n the water domain that encloses the shed vorticity and is equal to the integrated vorticity inside C. Furthermore, Us is the separation velocity just outside of the boundary layer, at the corner. The latter consideration assumes thin boundary layer. The low forward velocity of the body w i l l influence the local velocity at the edges of the moonpool. At the beginning of the oscillation cycle, when the wave amplitude inside the moonpool is increasing, the velocity introduced by the piston-mode oscillation w i l l on the leading edge have the same direction as the undisturbed incoming velocity. Due to the higher local velocity we w i l l expect that more vorticity is being shed f r o m this edge. The opposite w i l l happen on the trailing edge of the moonpool entrance. Here the local induced flow f r o m the increasing piston-mode amplitde w i l l have opposite flow direction as the undisturbed incoming velocity, and weaker vorticity is generated f r o m this edge, see Fig. 6c. Half a period later in the oscillation cycle the situation is turned around. The piston-mode amplitude is now decreasing and the vorticity shed f r o m the leading edge is now lower than that shed f r o m the trailing edge, see Fig. 6d. Due to the mentioned effects, how w i l l the viscous damping of the piston-mode ampli-tude change due to low forward velocity or incoming current? In the above discussion we have neglected the effect of vorticity shed from the leading edge of the hull. This vorticity w i l l influence the local flow around the entrance of the moonpool gap. The assump-tions about the magnitude of vorticity above might not be correct, with a sharp leading corner.

In addition to parameter variation in Figs. 9-14 and Sec-tions 4.1-4.4, a detailed parameter variation of the amplitude are given in Section4.5 (see Fig. 15). At last a velocity variation w i l l be given and results discussed in Section4.6 (see Fig. 16). To check the quality of the numerical simulations a convergence study has been

Table 1

Overview of figure numbers where results from the given combination of edge profile and draft are given. The appendage # 1 has dimensions 18 m m width by 9 m m height and appendage # 2 has dimensions 2 7 m m width by 1 8 m m height. The appendage # 1 cover 20% of the moonpool gap, while appendage # 2 cover 30% of the moonpool gap, see Fig. 7. L is a reference to the left part of the figure, and R is a reference to the right part of the figure.

Draft Edge profile Draft

Sharp corner Appendage # 1 Appendage # 2

15 cm Fig. lOL Fig. 12L Fig. 14L

18 cm Fig. 9 Fig. 11 Fig. 13

21 cm Fig. IOR Fig. 12R Fig. 14R

conducted. Including grid size, size of the potential domain inside the moonpool gap, time-step and smoothing. Details are shown i n Appendix A.

4.1. Results

We w i l l i n this section show results f r o m the parameter study on the 2D moonpool section. We w i l l show comparisons between numerical and experimental results for the parameter variation presented i n Tables 1 and 2. Note that we do not have experimental results for all the variations given in Tables 1 and 2.

Our main objective was to study the influence of a l o w for-ward velocity on the wave amplitude i n the moonpool. Secondly to validate the numerical method presented in this work. The flrst discussion is mostiy based on the experimental results, but w i l l be

Table 2

Overview of how the velocity and heave amplitude variation in Fig. 9 - 1 4 are orga-nized. The letters a - i refer to different parts in the figures.

ri3a Fn ri3a 0.00 ; 0.04 0.08 2.3 mm a b c 4.5 mm d e f 9.1 m m g h i

36 _{A. G. Fredriksen et al. /Applied Ocean Researcli 47 (2014)28-46}

d)

Fig. 8. Examples on how the leading edge wake influence the flow field at the moonpool entrance. The case here is from 9 with Fn = 0.04, >i2a =4.5 m m and r = 1.175 s. (a) is at the start of an oscillation when the boxes are moving upwards, (b) is at the top position, (c) is at middle position moving downwards, (d) is at the bottom position. Velocity arrows are given in the body-fixed coordinate system.

filled i n w i t h results f r o m the numerical work where experimental results are missing.

4.2. Experimental results

We begin by discussing the experimental results. A discussion regarding the comparison w i t h numerical results is given i n the next section. A first check of the quality of the experimental results is to compare the natural period of the piston-mode w i t h results given i n the literature.

Faltinsen et al. [5] provided accurate calculations of the natu-ral periods for the piston-mode, based on non-separated potential flow without current for case a) without appendages. Based on their

results the natural period for the 3 different drafts are T15 = 1.125s, l i s = 1.179s and T21 = 1.233s. The observed natural periods f r o m experiments w i t h o u t forward velcity were Ti5 = 1.13s, Tig = 1.18s and = 1.23s, and thus corresponds well to the results f r o m [5]. The natural periods were found not to change much when introduc-ing a forward velocity to the problem; it might have decreased a f e w percent when comparing the zero Froude number cases w i t h the 0.08 Froude number cases.

By taking Fig. 9 as an example the forward carriage velocity has a minor damping effect on the piston-mode amplitude. A decrease in the piston-mode amplitude between 5—7% f r o m Fn=0.0 to Fn=0.08 is observed, while for other appendage/draft configura-tions no decrease is observed.

A.G. Fredriksen et al / Applied Ocean Research 47 (2014) 28-46 37

a) Fn=0, amp=2.3mm b) Fn=0.04, amp=2.3mm c) Fn=0.08, amp=2.3nnm - Mp NonLin

10

d) Fn=0, amp=4.5mm e) Fn=0.04, amp=4.5mm f) Fn=0.08, amp=4.5mm -Outup NonLin

- Out down NonLin

g) Fn=0, amp=9.1mm ti) Fn=0.04, amp=9.1mm i) Fn=0.08, amp=9.1 mm

1.4

Fig. 9. Resutts for rectangutar side tiults (no appendage) at 18 cm draft, x , >: ; Experimental results for non-dimensional piston-mode amplitudes, o , o : Non-dimensional outgoing waves, dotted ( _ - ) line: Nonlinear numerical results for piston-mode amplitude. Dashed dotted (-..-.) line: Nonlinear upstream waves. Dashed dotted ( ) line: Nonlinear downstream waves. Dashed line with A: Linear numerical results for piston-mode amplitude. Dashed dotted ( ) line: Linear outgoing waves. All results are given in a body-fixed coordinate system.

As expected, the ratio between the piston-mode amplitude and the heave amplitude at resonance decreases as the heave ampli-tude increases. This is due to a quadratic increase i n the strength of the shed vorticity. However, strictly speaking we cannot assume quadratic dependence on heave velocity, but on local fluid velocity on hull edges.

There is a major increase i n damping f r o m sharp corners to appendages, but not much difference between the two appendages tested. The appendages have two effects, it is f r o m the experiments seen that they increase the natural period. From case (a) w i t h o u t appendages to case (c) w i t h large appendages the natural period is increased by 3%. The second effect is the increase i n the strength of shed vorticity, as more water is being pushed through a more narrow entrance.

The effect of changing the draft is i n practice negligible for the maximum value of the piston-mode amplitude, as seen by comparing Figs. 9, 10, 11, 12 and 13, 14. Our simplistic view on this is that the heave morion displaces a certain amount of water which is proportional to the beam of the section. During resonance a large part of that water goes into the moonpool. How large the piston-mode amplitude becomes is then depen-dent on the ratio between the w i d t h of the moonpool and the length of the side hulls, and not much affected by the draft.

The asymptotic value when T^-^Os for both the piston-mode elevation and the outgoing wave amplitude should approach 0, which would be equal to 1 i n the body-fixed results i n Figs. 9-14.

38 A.G. Fredriksen et al./Applied Ocean Research 47 (2014) 28-46

aJ Fn=0, amp=2.3rrm '',s' F"=0-<'4. amprf.Smm c,^) Fn=0.08, amp=2.3mm Fn=0, amp=2.3mm b^,) Fn=0.04, ainp=2.3mm o^^) Fn=0.oe, amp=2.3mm

../'t^v-.i.

d lFn=0,amp=4.5mni e „ ) Fn=0.04, anip=4,5mm f,^) Fn=0.0e. amp=4,5mm dj,) Fn=0, amp=4.Emm e^^) Fn=0.04, amp=4.6mni f^^) Fn=0.08, amp=4.5mm

g,5)Fn=0, amp=9.1mm h,j) Fn=0.04, amp=9.1mni 1^^) Fn=0.08, amp=9.1mm g^^) Fn=0, amp=9.1 mm h^^) Fn=0.04, amp=9.1mm 1^^) Fn=0.08, amp=9.1mm

1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4

T[s]

T[s]

T[s]

T[s]

T[s]

T[s]

Fig. 1 0 . Results for rectangular side hulls (no appendage) at 15 cm (left) and 21 cm (right) draft. See caption and legends in Fig. 9 for description of symbols.

4.3. Numerical setup

The grid has the following properties: In the horizontal direction the grid size around the hull used in all calculations was 0.01 m, w h i c h means that there are 36 cells i n the horizontal direction over each hull and 18 cells across the moonpool. This grid size extends 1.2 hull lengths or 0.3 wave lengths away f r o m the hull depending on what is shortest in upstream direction. In the down-stream direction direction of the hull, the grid size is constant for 1.5 hull lengths or 0.42 wave lengths away f r o m the hull also depend-ing on what is shortest. Outside this area the grid size gradually increases to 30 cells for each wave length. In the damping zone the grid size increases even further. Note that the horizontal grid size are then constant in the vertical column. The wave length used here is computed without current effects, but w i t h water depth effects.

In the vertical direction the grid used depends on the draft and expected piston-mode amplitude. Generally the vertical grid size is 30% smaller than the horizontal grid size around the hull. The grid size is Icept constant until 0.5 hull draft below the bottom of the hull, then gradually increasing until the bottom of the tank, w i t h a total number of 60 cells in vertical direction. Also here the vertical grid size is kept constant for each horizontal row. A problem in using an easy grid generation method as this, is that the aspect ratio of the cells far away f r o m the body (vertically or horizontally) becomes high (or low).

The damping zone starts, at both sides, four wave lengths ( w i t h -out current) away f r o m the body and increases smoothly over one wave length to its maximum value (0.8), and is kept at this value for another three wave lengths before the end of the tank. The total length of the numerical domain is 16 wave lengths i n addition to the length of both hulls including the moonpool.

The time-step A t is chosen to be the lowest of 0.5 times the Courant-Friedrichs-Lewy (CFL) number ( A f = O.SAy/f or A t = 0 . 5 A z / w ) or 120 time-steps per oscillation period A t = 1/120.

The height of the potential f l o w domain inside the moon-pool is for all drafts and appendages set to be 0.03 m for 2.3 m m heave amplitude, 0.06 m for 4.5 m m heave amplitude and 0.075 m for 9.1 m m heave amplitude. These values are based on what is observed f r o m experiments, and included a safety margin to allow over-prediction of the numerical results.

The wave profile in the moonpool gap should be almost flat when the sloshing modes are not present. These sloshing modes w i l l be excited by transient effects in the start-up. These transient sloshing modes are damped out by a numerically added artificial damping around the mean position of the wave profile inside the moonpool gap, this is done i n a similar way as for the damping applied i n the numerical beach and which is done in sloshing anal-ysis by the multimodal potential f l o w method (see [4]). Physically the only damping source of the sloshing modes are through the boundary layers on the glass walls and on the side hulls. The damp-ing ratio ^ on the first sloshdamp-ing mode for the current setup f r o m the

A.C. Fredriksen et al./Applied Ocean Research 47(2014) 28-46

a) Fn=0, amp=2.3mm b) Fn=0.04, amp=2.3mm c) Fn=0.08, amp=2.3mnn

39

d) Fn=0, amp=4.5mm e) Fn=0.04, amp=4.5mm f) Fn=0.08, amp=4.5mm

g) Fn=0, amp=9.1mm h) Fn=0.04, amp=9.1mm i) Fn=0.08, amp=9.1mm

V ^ . .

0 5 ^

T[s]

T[s]

Fig. 1 1 . Results for appendage # 1 (18 mm x9 mm) at 18 cm draft. See caption and legends in Fig. 9 for description of symbols.

boundary layer flow is found f r o m [4] to be ^ = 0.0014. It is believed that this artificial damping inside the moonpool gap is larger than what w i l l be i n reality due to boundary layer damping.

4.4. Numerical results

The validity and limitations of the numerical method presented here are discussed. First we discuss the presented nonlinear body-fixed hybrid method results. W i t h the limitations given i n the numerical method, we believe that the results compare quite well w i t h experiments (see Fig. 9-14). The simulations for higher ampli-tudes somewhat over-predict the damping of the gap amplitude, while it for lower amplitudes are i n good agreement w i t h the exper-imental results. For higher heave amplitudes, the ratio between the piston-mode amplitude and the draft increases. This requires a

larger potential flow domain due to re-gridding. This means that, while the amount of vorticity increases, we have to decrease the viscous domain inside the gap. The consequence is an increase in the probability that vordcity w i l l reach the numerical intersection between the viscous and the potential flow domain. Therefore our numerical method is limited to free shear layers that stays below the level of the wave trough i n the body-fixed coordinate system. (See results from the sensitivity tests of the height of the potential domain inside the moonpool in Fig. 18.)

Also the vorticity separated f r o m the trailing edge on the down-stream side w i l l grow w i t h increasing current velocities, up to a point where the flow separation w i l l cause a "dry transom stern". By "dry transom stern" we make the analogy to what happens for high-speed semi-displacement and planning vessels. However, the dry flow separation case is not relevant here, as

40 h.G. Fredriksen et al. /Applied Ocean Research 47(2014) 28-46

a,j) Fn=0, amp=2.3mm ^ s ' amp=2.3mm oJ Fn=0.08, amp=2.: a^,) Fn=0, amp=2.3rtim b^^) Fn=O.04, amp=2.3mm o^,) Fn=0.08, amp=2.3mm

d,j) Fn=0, amp=4,5mm F"=«-0''. anip=4.5mm f,^) Fn=0.0a, amp=4.5mni dj,) Fn=0, a[iip=4.6mm Fn=0.04, amp=4.5nnm f^,) Fn=0.0e, amp=4.e

9,5) Fn=0, amp=9.1mm h^^) Fn=0.04, amp=9.1mm 1^^) Fn=0.oe, amp=9.1mm g ) Fn=0, amp=9-1 mm h ) Fn=O.04, amp=9.1 mm i ) Fn=0.08, amp=9.1 mm

0.8 1 1.2 1.4 0.8 1 1,2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4

T[s]

T[s]

T[s]

T[s]

T[s]

T[s]

Fig. 1 2 . Results for appendage # 1 ( 1 8 m m x 9 mm) at 15 cm (left) and 21 cm (right) draft. See caption and legends in Fig. 9 for description of symbols.

it w i l l happen for much higher Froude numbers than what we consider (Fn >0.3—0.4 [3]). However, for velocities lower than for dry separation the vorticity w i l l reach the numerical inter-section between the potential flow and viscous f l o w domains, and later the vorticity w i l l i n the experiments reach the free sur-face.

The numerical results for appendage # 1 compare better w i t h experiments than the numerical results for appendage # 2 (see Figs. 11-14). The numerical results and experimental results differ in both amplitude and the predicted period of maximum response. The mesh is almost identical for the t w o cases and the intersec-tion is at the same locaintersec-tion, so there must be some physical effect our numerical method does not capture. For appendage # 2 the edges are closer to each other than they are for appendage # 1 , such that the vorticity created at one edge is more likely to influ-ence to vorticity created at the other edge. Further, our viscous model is not capable of dealing w i t h the turbulent mixing of vor-ticity.

As a first estimate on the influence of turbulent diffusion in our problem, we simulated w i t h our nonlinear method cases where we increased the dynamic viscosity v f r o m 10"^ kg/(ms) to 10-'' lcg/(ms) in the fluid. V\/e increase v to simulate that there exists turbulent diffusion due to eddy viscosity at scales smaller than what we capture w i t h our grid density. An increase i n v w i l l result in an increase in the boundary layer thickness, that w i l l increase the damping of the piston-mode motion. For v = 10"^l<g/(ms) we have seen that the damping contribution f r o m the boundary layer is small, this may not be the case w i t h v = 10-''l<:g/(ms). For the

lowest heave amplitude (2.3 m m ) we see a decrease i n the piston-mode motion, around 10% at resonance. While for the highest heave amplitude (9.1 m m ) the piston-mode response is not much affected by increasing the dynamic viscosity. This could be a result of that we for some areas in the fluid have a turbulent mixing of vorticity, and thus a higher diffusion/cancellation of vorticity w h i c h would lead to a smaller damping of the moonpool piston-mode resonance. Piston-mode amplitude results f r o m the linear hybrid method are higher than the nonlinear results for all zero current cases. Further, it over-predicts the gap response for higher heave ampli-tudes. It is believed that the reason for this is that the potential f l o w domain i n the linear results can be minimized to only contain the top layer close to the free surface. While the potential f l o w domain in the nonlinear method needs to be larger due to re-gridding of the free-surface, since the implementation does not allow the viscous CFD domain to change size i n time. Notice that also the linear results are different f r o m the experiments for appendage # 2 .

By studying velocity plots f r o m our numerical simulations w i t h forward velocity (see examples i n Fig. 8), we see that the lead-ing edge on the upstream side hull creates a wake w h i c h make the relative velocity at the moonpool gap entrance almost zero. In effect, the water i n the moonpool is free to oscillate nearly w i t h o u t the influence of forward velocity. This may explain the negligible influence of low forward velocity on the moonpool behavior. Also the free-vorticity f l o w that develops at the leading edge w i l l easily become turbulent, which again influences the diffusion of vorticity that w i l l not be captured by our laminar model.

A.G. Fredriksen et al. /Applied Ocean Research 47(2014} 28-46 ^'^^'^'^'""^ t>) Fn=0-04, amp=2.3mm o) Fn=0.08, amp=2,3mm 41 -X d) Fn=0, amp=4.5mnn e) Fn=0.04, amp=4.5mm

: # \

T[s]

T[s]

Fig. 1 3 . Results for appendage # 2 (27 m m x l S mm) at 18 cm draft. See caption and legends in Fig. 9 for description of symbols.

4.5. Amplitude variation

For each of the 3 edge profile configurations, we performed an experimental investigation where we varied the heave amplitude, ranging f r o m 2 m m to 24 mm, see Fig. 15. During these test-series' the draft was kept constant to 18 cm and constant forward veloc-ity corresponding to Fn = 0.04. However, the period was varied and each amplitude was tested w i t h 5 different forcing periods around the natural frequency. The results show as expected that the amplitude in the gap is influenced by a nonlinear damping, and further, that the nonlinear damping is higher for the two cases w i t h appendages. It is quite striking that the damping introduced by the two different appendages is so similar. The damping is con-siderably increased relative to square inlet, up to 6 0 — 7 0 % for the largest forcing amplitudes, which is of appreciable practical impor-tance.

4.6. Velocity variation

To isolate the effect of forward velocity on the piston-mode amplitude, we performed an experimental test-series where we I<ept all parameters other than the carriage velocity constant. The tests were done w i t h o u t appendages and using a draft of 0.18 m . The carriage velocity was varied f r o m Fn =0.01 up to Fn =0.11. Note that for velocities above Fn -0,08, i t is somewhat doubtful that the experiments have reached steady-state. The results are presented in Fig. 16 where the forcing heave amplitudes 7730 were fixed at 4.5 m m and 9.1 mm. The period was meant to be kept constant at 1.18 s. However, note that the data analysis after the tests showed that the period had not been constant between tests, w i t h a relative error of up to 2%. A linear curve is fitted to the data between Fn =0.01 and Fn =0.Q8, see Fig. 16. The piston-mode response decays slightly for increasing forward velocity, a 7%

42 A.G. Fredriksen et al. /Applied Ocean Research 47(2014) 28-46

a ) Fn=0, amp=2.3mm b^j) Fn=0.04, amp=2.3mm o^^) Fn=0.06, amp=2.3mm a^,) Fn:=0, anip=2.3mm b^,) Fn=0.04, amp=2.3mni c^^) Fn=O.0e, amp=2.3mm

S

d ) Fn=0, amp=4.5mm e,^) Fn=0.04, amp=4.5mm f^^) Fn=O.OB, amp=4.5nim d^,) Fn=0, amp=4.5mra e^,) Fn=0.04, amp=4.5mnn y Fn=0.08, anip=4.6nim

5 '

^ 4

A '

g,j)Fn=0,amp=9.1mm h^^) Fn=0.04, amp=9.1mm 1^^) Fn=0.08, amp=9.1mm g^^) Fn=0, amp=9.1mm h^,) Fn=0.04, annp=9.1mm 1^,) Fn=0.08, amp=9.1mm

0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4 0.8 1 1.2 1.4

T[s] T[s] T[s] T[s] T[s] T[s]

Fig. 14. Results for appendage # 2 (27 m m x 18 m m ) at 15 cm (left) and 21 cm (right) draft. See caption and legends in Fig. 9 for description of symbols.

o § 0

O O O

Heave variation, Fn=0.04

O Rectangular side tiull Appendage #1 Appendage #2

°.: ^.1 l

Fig. 15. Non-dimensional piston-mode amplitudes in a body-fixed view as a func-tion of forced heave amplitude. For cases w i t h and without appendages. Each amplitude is tested for 5 different forcing penods.

Velocily variation 0 : 0 ' _{o • o} o ; o : • - - Linear fit ^ 11 =9.1 mm '001 0 02 0 03 0 04 0.05 0.06 0.07 0.03 0,09 0.1 0.11 Fn

Fig. 16. Non-dimensional piston-mode amplitudes in a body-fixed view as a func-tion of forward velocity for constant heave amplitude ))3o of 4.5 m m and 9.1 mm, case with no appendage. The data is fitted to a linear curve for Froude numbers between Fn =0.01 and Fn =0.08.

A.G. Fredriksen et al. / Applied Ocean Research 47 (2014) 28-46 43

decay f r o m zero forward velocity to Fn =0.08 is seen for the forcing amplitude ??3n=4.5mm. The decay from forward velocity is somewhat lower for higher amplitudes of oscillation.

5. Conclusions

The main motivation behind this study was to investigate how a low current/low forward velocity influences the resonant piston-mode motion of a moonpool inside a 2D ship section. This was done by both numerical and experimental methods, and for cur-rent velocities up to Fn = 0.08. We found that w i t h i n this range, the low current/low forward velocity had a slightly decreasing effect on the moonpool piston-mode behavior. To what extent this fact depends on the gap w i d t h remains unknown.

In addition to the low forward velocity influence, we varied the moonpool edge profile, the draft and the heave amplitude. We tested 3 different edge profile inside the moonpool gap. First rectan-gular side hulls, then including two appendages covering 20% and 30% of the moonpool gap area. The damping of the piston-mode is increased by using appendages compared to the rectangular side hull. For each of the 3 edge profiles we tested 3 different drafts, where the ratio between the draft and the total hull w i d t h was 1/6, 1/5 and 7/30. However, the change of draft had little influence on the maximum piston-mode response. It did change the period of maximum response of the system, such that for a given period the response was changed.

We demonstrated how to use a direct coupling between a poten-tial f l o w solver and a viscous f l o w Navier-Stokes CFD solver. The coupling was based on the newly developed Harmonic Potential Cell (HPC) method f r o m [20] for the potential f l o w domain, and a viscous solver by the commonly used Finite Volume Method (FVM) by assuming laminar flow. The overall agreement is good between the numerical and experimental results. There are, however, some deviations which are thought to be caused mainly by that vorticity reaches the interface between the potential and viscous domains. An important finding and a fact that has been recognized by many others before, is that f l o w separation at the moonpool gap

entrance provides important damping of the resonant moonpool piston-mode oscillation which necessitates that a viscous solver has to be used. However, a large part of the domain can be described by potential flow.

Acknowledgements

We acknowledge Fredrik Dukan for setting up, implementing and testing the automatic control system for the carriage. This work is supported financially by the Research Council of Norway through Centre for Ships and Ocean Structures (CeSOS).

Appendix A. Convergence

We w i l l here present four parameters we identified as possible candidates that are important for the convergence of the nonlinear hybrid method. The four parameters are grid size, size of the poten-tial domain inside the moonpool gap, time-step size and smoothing size. Convergence is studied for all appendage configurations w i t h draft 18 cm, 4.5 m m heave amplitude and t w o Froude numbers (Fn = 0.00 and Fn = 0.08). We have chosen five different forcing periods around the natural period for each case. The parameters are varied around the inputs used for the calculations i n Figs. 9-14, from here on called the standard case. Figs. 17-20 all have six sub-figures, i n t w o rows and three columns. In the first row are convergence results for Fn =0.00, and second row are results for Fn =0.08. The first column contains results for the rectangular side hull profile (no appendage), second column for appendage # 1 and third column for appendage # 2, see Table 3.

The convergence and sensitivity results w i t h forward velocity should all be viewed l<eeping i n mind that vorticity has reached the intersection between potential and viscous flow. This means that the convergence results could show h o w the local inaccuracy on the intersection is affected by a change i n grid and/or time-step size and how this inaccuracy influence the piston-mode oscillation, and not how the change in grid and time-step size change the global flow field, and then the piston-mode oscillation.

a) No App, Fn=0.00 b) App #1, Fn=0.00 4.6^ ~ ~ + ~, : 4.5 4.4 4.3 4.2 4.1 1 2 1.21 1.22 1.23 1.22 1.23 1,24 1.25 1.26

T[s] T[s] '

Fig. 1 7 . Convergence of the piston-mode amplitude with respect to the number of grid-ceils across one side hull, the number is indicated in the figure. Here 36 is the case corresponding to the "standard case" results in Figs. 9-14. The number of cells in the vertical direction is varied respectly, to keep the aspect ratio A z / A y constant.