Hamiltonian Boussinesq formulation of wave–ship interactions

Contents lists available at ScienceDirect

Applied

Mathematical

Modelling

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

Hamiltonian

Boussinesq

formulation

of

wave–ship

interactions

E.

van

Groesen

a, b, ∗

_,

_Andonowati

c, b

a Department of Applied Mathematics, University of Twente, Netherlands b LabMath-Indonesia, Bandung, Indonesia

c Department of Mathematics, Institut Teknologi Bandung, Indonesia

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 21 December 2015 Revised 10 July 2016 Accepted 6 October 2016 Available online 18 October 2016

Keywords: Wave–ship interaction Hamiltonian formulation Boussinesq description Momentum exchange

a

b

s

t

r

a

c

t

Inthispaperanewapproachisdescribedforthefullynonlineartreatmentofthedynamic wave–shipinteractionforpotentialflows.Amajorreductionofcomputationalcomplexity isobtainedbydescribingthefluidmotioninhorizontalvariablesonly,thesurface eleva-tionandthe potentialatthesurface.InsuchBoussinesqtypeofequations,theinternal fluidmotionisnotcalculated,butmodeledinaconsistentapproximativeway.The equa-tionsforthewave–shipinteractionarebasedonaLagrangianvariationalprinciple,leading totheformulationofthecoupledsystemasaHamiltoniansystem.Withtheshipposition andorientationas canonicalcoordinates,thecanonicallyconjugatemomentumvariables arethe sumoftheship momemtaand thefluidmomenta.A beneficialconsequenceof thisisthatthe momentumexchangebetweenfluidand shipwillbedescribed without theneedtocalculatethepressure,whichsimplifiesthenumericalimplementationofthe equationsconsiderably.ProvidedthatthepotentialswithmixedDirichlet–Neumanndata canbe calculated,the presentedship dynamicscanbe insertedinexistingfree surface flowsolvers.

© 2016ElsevierInc.Allrightsreserved.

1. Introduction

In this paper the Hamiltonian description of inviscid free surface waves will be extended in a straightforward consistent way to include ship–wave interaction. For inviscid, irrotational flows with a free surface, the Hamiltonian formulation of the continuity and Bernoulli’s equation at the surface uses the potential

φ

and the surface elevation

η

as canonical conjugate variables; this result based on contributions by Luke [1], Zakharov [2], Broer [3]and Miles [4]is summarized in Section2.2. The dependence of

φ

,

η

on horizontal variables alone makes that the Hamiltonian formulation is of Boussinesq-type, i.e. dimension reduced.

For the coupling of a rigid ship with waves, Lagrangian and Hamiltonian methods have been used before, and the general equations are being used in numerical implementations. The use of Lagrangian principles for solid–inviscid fluid interaction started with Thomson and Tait, and Kirchhoff, see Chapter 6 of Lamb [5]. Miloh [6]seems to be the first to include the radiation potential to describe linearized surface elevation caused by forced objects. In Van Daalen [7]a concise derivation of the Lagrangian and Hamiltonian formulation is given for the water–ship interaction in incoming waves. Nowadays, the

∗ _{Corresponding author at: Department of Applied Mathematics, University of Twente, Netherlands}

E-mail addresses: groesen@math.utwente.nl , e.w.c.vangroesen@utwente.nl (E. van Groesen). http://dx.doi.org/10.1016/j.apm.2016.10.018

Lagrangian–Hamiltonian formulation seems standard, and has been extended to the coupling of water motion and flexible bodies by Xing and Price [8]; see this paper also for references to the recent literature.

Yet, listing the total set of the dynamic equations combined with the various Laplace problems that have to be solved, may somewhat hide the variational structure behind the equations, which may, at best, cause that not full advantage is taken of the peculiarities of the special structure. This is particularly evident in the treatment of the momentum exchange between ship and fluid, and the role the pressure is given. In the Hamiltonian formulation of this paper, it is found that

thedynamicchangeofthecanonicalmomentum,whichisthesumoftherigidbodyandthefluidmomentum,equalsthestatic hydrodynamic force. Written in a Newtonian way for one motion component, the canonical momentum is ps+ pf where

psis the ship momentum; pfis the fluid momentum, i.e. the integral over the wetted ship hull of the fluid potential. The Hamilton equation is then of the form, writing

τ

H _{for the hydrostatic forces:}

d

dt

(

ps+pf

)

=

τ

H_{. (1)}

Note that dpf/ dt is (minus) the dynamic part of the pressure; in this Hamilton equation it is forced by the hydro-static part of the pressure. In the explicit time stepping procedure that we will derive these two momentum parts are treated differently.

However, in contrast, characteristic of seemingly all treatments of wave–ship interaction is the fact that this balance law is described by a formulation that thechangeintheshipmomentumequalsthefluidforcesactingontheship; see for instance the formulation for the ‘Dissel’ software by Lin and Kuang [9–11]. In the simplified notation this reads:

d dt

(

ps

)

=

d

dt

(

−pf

)

+

τ

H_{, (2)}

where the right hand side is now the total pressure. With this formulation, an explicit time stepping gives a different result; an implicit scheme will lead to the known problems to update the pressure to the same time level as the potential itself. It will be shown that the Hamiltonian formulation with explicit time stepping avoids the calculation of the pressure, and will consistently respects the order in which forces work.

The description of the Hamiltonian formulation of the coupled system in Section3finds its basis in work of Van Daalen [7], see also Van Daalen et al. [12]. The main modifications presented here concern the explicit expression of the equations in the Hamiltonian variables instead of using the pressure, which leads to the time evolution that avoids the calculation of the pressure.

Concerning the interior flow that has to respect Dirichlet and Neumann conditions, the total potential will be split in a way as described e.g. by Cummins [13]and Ogilvie [14]. One component of the total potential is the fluid potential with normal derivative vanishing at the ship and Dirichlet value at the free surface; it describes the tangential flow along the ship. The other part is the instantaneous radiation potential that is the sum of normalized instantaneous radiation potentials multiplied by the corresponding velocities for each of the six degrees of freedom. The evolution of the impulsive radiation potential is part of the Hamiltonian evolution, without the necessity to consider the convolution as described in Cummins, see Section4.3.

Solving the interior potentials to obtain their effect on the free surface and at the ship hull can be done by using any reliable Laplace solver. The capabilities of the Laplace solver will determine directly the applicability of the wave–ship sys- tem, such as simulations above bathymetry, in harbors, and for ships with small under keel clearance. Results of the Laplace solver can be used directly in the governing Hamilton equations. The wave potential has to be updated each time step; for simplified ‘linearized’ ship interaction, it is sufficient that the normalized radiation potentials (and related added-mass coefficients) are calculated only once with the ship in its hydrodynamic equilibrium position. When aiming to keep the advantage of efficient calculations based on the dimension reduction in the Boussinesq description, the calculation of the effect of the potentials at the free surface and at the ship will have to be based on an approximation of the interior flow, which requires an extension of the methods used in Boussinesq models for free surface flows. In forth coming publications we will show the performance of the coupled system using Boussinesq reduction in two models of HAWASSI software [15]. The organization of this paper is as follows. In the next section we set notation and briefly describe the derivation of the Hamiltonian formulation for the free surface flows. The generalization to include water-body interaction will be described in Section3. We restrict to the case of one free floating ship; the case of more or moving ships is a rather direct extension. In Section 4the potentials and momenta that have to be calculated are summarized, and the time stepping for the Hamiltonian evolution is described. Conclusions are formulated in the final section.

2. Hamiltonianwavedynamics

The dynamical system for the combined wave–ship motions will be an extension of the Hamiltonian description of free- surface waves which will be shortly described in this section. Just as in Classical Mechanics for interacting rigid bodies, the total energy is the main quantity, expressed in canonical variables that form the state variable that has to be evolved in time. In the first section we set notation, followed by the Hamiltonian formulation for waves only. The Lagrangian description for the wave–ship interaction, from which the Hamiltonian formulation for the water-body interaction will be derived in Section 3, is a generalization of the free surface Lagrangian. In the last section we deal with the practical aspect how to obtain the required kinetic energy as explicit expression in the canonical variables for waves above flat and varying bottom.

2.1. Notation

The fixed earth Cartesian coordinate system is denoted by x=

(

x,y,z

)

with x,y the horizontal coordinates perpendicular to the z axis in the opposite direction of gravity; e_jfor j₌1 _,2 _,3 denote the unit vectors along these axes. The bathymetry is described by the depth z= D

(

x,y

)

. The surface elevation is described by

η

( x,y,t), the fluid potential by



( x,y,z,t) and the potential at the free surface by:

φ

(

x,y,t

)

≡

(

x,y,z=

η

(

x,y,t

)

,t

)

.

In the following we assume for ease of presentation that the mass density

ρ

of water is constant, taken to be

ρ

=1 . Furthermore, the word ‘potential’ will be used for quantities that have to satisfy the Laplace equation in the fluid domain and that have vanishing normal derivative at the bottom, or vanishing gradient in the limit of infinite depth. Depending on the case, vanishing gradients are required at lateral boundaries far away, or specific inflow conditions are prescribed there.

The notation to be used is that

∇

2is the gradient in horizontal variables, and that n=N/



1 ₊

|

∇

2

η|

2and N=

(

−

∇

h

η

,1

)

are the normalized and non-normalized normal, with

∂

n



and

∂

N



the corresponding notation for normal derivatives.

2.2.Freesurfaceformulation

Starting point is the description of free surface flows. Bateman [16]remarked already that variations of the potential in the integrated pressure:

− 



∂

t



+ 1 2

|

∇

3

|

2 + gz



dxdydz, (3)

produces the condition that



should be a potential. Luke [1]considered the expression in its dependence on both



and the free surface

η

variable, and showed that the correct two nonlinear free surface conditions are obtained in addition. In fact, Miles [4] observed that by taking the partial time derivative of the dynamic pressure part outside the depth integration:

_η

−D

∂

t



dz=

∂

t _η

−D



dz−

φ∂

t

η

,

the irrelevant total time derivative can be discarded because it will not produce a contribution to the governing equations since variations at the initial and final time are supposed to vanish in the formulation as a canonical action principle. This then leads to the Lagrangian:

L= 

φ∂

t

η

dxdy− H

(

φ

,

η

)

, (4)

where H is the Hamiltonian, the sum of kinetic and potential energy:

H

(

φ

,

η

)

≡ K

(

φ

,

η

)

+P

(

η

)

. (5)

The kinetic energy expressed in the surface variables is determined by the potential satisfying the Dirichlet condition at the surface and impermeability of the bottom (and lateral conditions):

K

(

φ

,

η

)

=



_ 1 2

|

∇

3

|

2 dxdydz









=0 for z<

η



=

φ

at z=

η

∂

N



=0 at z=−D

(

x,y

)



. (6)

The potential energy depends on the surface elevation as:

P

(

η

)

= 

1 2g

η

2_{dxdy. (7)}

The free surface equations in Hamiltonian form are then given by the Hamilton equations as were derived earlier in different ways by Zakharov for infinite depth [2]and Broer [3]for finite depth:

∂

t

η

=

δ

_φH

∂

t

φ

=−

δ

ηH. (8)

Here and in the following

δ

_φH denotes the variational derivative of H with respect to

φ

, and similarly for

δ

_ηH.

To show that this is equivalent with the usual nonlinear boundary conditions at the free surface, the variational derivatives of the kinetic and potential energy have to be calculated. For variations

δη

the variation of the potential energy is found from:

P

(

η

+

εδη

)

=P

(

η

)

+

ε



g

η

.

δη

dx+O

(

ε

2

)

,

from which

δ

_ηP₌g

η

_. The kinetic energy is a bit more complicated. At fixed

η

, consider a variation

δφ

in the surface potential. Let

δ



be the corresponding variation of the potential in the interior. The we get for small

ε

:

Using Green’s identity: 

∇

3



·

∇

3

δ



dxdydz=− 

3



·

δ



dxdydz+

∂

n



·

δ



,

this reduces for potential functions



3



= 0 in the interior, to the boundary integral which vanishes on bottom and lateral boundaries. Hence there remains the surface integral:

K

(

φ

+

εδφ

,

η

)

=K

(

φ

,

η

)

+

ε



∂

N

|

z=η·

δφ

dxdy+O

(

ε

2

)

,

so that:

δ

φK

(

φ

,

η

)

=

∂

N

|

z=η.

For the variational derivative with respect to

η

, the result is:

δ

ηK

(

φ

,

η

)

=



1 2

|

∇

3

|

2 −

∂

z



∂

N





z=η,

the first term being obtained by varying the surface for fixed



, and the second term as correction because the surface potential has to retain the same value at the deformed surface.

Written in full the Hamilton equations become:

∂

t

η

=

δ

_φH=

∂

N

|

z=η

∂

t

φ

=−

δ

ηH=−

g

η

+1 2

|

∇

3

|

2₋

_∂

z



∂

N





z=η. (9)

The first equation is the usual continuity equation. The second equation is the Bernoulli equation at a pressure free surface. This can be better recognized in its common form which is obtained by using the fact that at z=

η

it holds that

∂

t



=

∂

t

φ

−

∂

z



∂

t

η

, so that grouping the terms together and using the continuity equation there results:

∂

t



+g

η

+

1 2

|∇

3

|

2_=0atz=

_η

_.

Having recovered the correct continuity and momentum equation at the free surface, this justifies the Hamiltonian formulation as given above.

2.3. Kineticenergyapproximations

The variational and Hamiltonian description in the previous section requires the expression of the kinetic energy to be explicit in terms of the canonical variables

φ

,

η

. Dirichlet’s principle (see e.g. Treves [17,Chapter 3]) expresses that the kinetic energy can be obtained from solving the minimization problem:

K

(

φ

,

η

)

=min 

{

D

()

|



=

φ

atz=

η}

, (10) with: D

()

=  1 2

|

∇

3

|

2 dxdydz, (11)

the Dirichlet functional. Under mild conditions for

η

,

φ

and the depth, mainly smoothness of the surface z₌

η

and the bottom D( x) and decay of

∇

3



at lateral boundaries, a unique solution of the minimization problem exists which then satisfies the Laplace equation



= 0 in the interior, the bottom and lateral boundary conditions and the Dirichlet condition at the surface. Hence K₌K

(

φ

_,

η

)

is well defined for finite varying depth and also for infinite depth. Even for the special case of a flat bottom, an explicit solution of the Laplace boundary value problem can only be obtained in closed form in the case of a flat surface (for the limit of vanishingly small waves) in the form of a Fourier superposition of Airy depth profiles:

(

x,z

)

=  ˆ

φ

(

k

)

coshk

(

z+D coshkD e ikx_{dk, (12)}

with

φ

ˆ the spatial Fourier transform of

φ

. For infinite depth, the quotient in the integrand is replaced by its limiting value into e|k|z_.

For practical purposes an explicit expression for K(

φ

,

η

) is therefore only possible in an approximate way. There seems to be essentially one way to achieve that in a consistent way, and that is to look for a good approximation of the fluid po- tential



explicitly expressed in

φ

,

η

(and depth D) and then substitute that approximation in the Dirichlet functional. This guarantees that the approximate expression remains positive definite and can serve as an explicitly defined approximate kinetic energy.

In the past years this methodology has been used to obtain so-called VariationalBoussinesqequations that turned out to be well capable to deal with nonlinear dispersive waves over varying bottom. As an illustration, starting with the simplest approximation, taking

(

x,z

)

=

φ

(

x

)

independent of the vertical coordinate, there results K

(

φ

,

η

)

= 1

2  

(

D+

η

)

|∇φ|

2

dx

one or more depth dependent profiles, the dispersive properties can be improved and optimized for the case under investigation, see Klopman et al. [18,19], Adytia and Van Groesen 2012 [20]and Lakhturov et al. [21].

Somewhat different, preserving the Hamiltonian structure, but possibly losing the positive definiteness, is to use Taylor expansions in wave height and inverse wave length to find approximations, see for the lowest order of approximation Broer [22] and Broer et al. [23]. Craig and Sulem [24] approximated Zakharov’s formulation up to fifth order accuracy by a Taylor expansion of the Dirichlet-to- Neumann operator that maps the fluid potential at the fluid surface to the normal derivative of fluid potential at the surface; approximations up to fifth order in the equations are also given in Kurnia and van Groesen [25,26]that lead to pseudo-spectral implementations with Fourier integral operators that can also deal with varying bottom. Simulations compared to data of measurements show remarkably good performance of these modern codes in a large variety of nonlinear waves above flat and varying bottom as shown in the recent references mentioned above.

3. Wave–shipinteraction

The dynamical system for the combined wave–ship motions is an extension of the Hamiltonian description of free- surface waves in the previous section. In the first section we set notation, followed by the description of the wave–ship interaction by Lagrangian principles. Then the Hamiltonian formulation for the water-body interaction will be derived.

The generalization to include a ship in the Hamiltonian wave formulation proceeds along a number of steps.

Rewriting the pressure principle in the previous section led to the canonical action principle and to the Hamiltonian formulation of surface waves. A ship as a rigid body in the gravitational force filed, is a standard example of a system from Classical Mechanics for which the Lagrangian formulation is well known. To couple the wave dynamics with the ship dynamics the two action principles will be combined, where the fluid dynamic part of the free surface formulation is adjusted for the presence of the ship. The free surface F given by z=

η

(

x,y,t

)

being pierced by a rigid body, the wetted part of the ship hull S has to be taken into account. The term in the canonical action principle  

φ∂

t

η

d xd y is modified at the position of the ship by replacing the normal free surface velocity

∂

t

η

by the ship velocity in the normal direction which has to be compatible with the normal derivative of the fluid motion at the ship hull.

This coupled action principle then leads to the correct equations. A Hamiltonian form is obtained by defining the canonical momentum which will be the sum of the ship momentum and the total fluid momentum. The time derivative of the total fluid momentum, which is the dynamic pressure, is in this formulation not a driving force, but is updated in a time step as part of the canonical momentum under influence of the hydrostatic force.

After setting notation, the Lagrangian formulation and the introduction of the momentum will be detailed, after which the fully coupled Hamiltonian system is derived. In the next section, details of the Hamiltonian evolution and time stepping are described.

3.1. Notation

The ship has six degrees of freedom; the position of the ship and its orientation are determined by the center of gravity X=

(

X,Y,Z

)

and three angles

θ

, forming together the 6-dimensional position vector

ξ

=

(

X,

θ

)

. Points at the hull are specified by a body vector r with respect to the center of gravity as xS =X+r; variations in the ship position and orientation affect points at the ship hull as:

δ

xS=

δ

X+

δ

r=

δ

X+r×

δθ

, in particular x˙S = X˙ + r× ˙

θ

.

Explicitly, we assume that the wetted part of the ship hull S can be described as a function z=

ζ

(

x,y,t

)

; the character- istic function

χ

S (equal to 1 in the ship area, zero outside) and its complement for the free surface

χ

F = 1 −

χ

S are then well defined. We will write S0, F0 for the corresponding areas in ( x,y)-space. We use for the normal at the ship hull the same notation n=N/

|

N

|

with N=

(

−

∇

2

ζ

,1

)

, and introduce the six-dimensional generalized normal

ν

:

ν

=

(

n,n× r

)

. (13)

The potential restricted to the wetted part of the ship



| S will be denoted by

φ

, occasionally also denoted by

φ

S to distinguish from the potential at the free surface

φ

F =



|

F. It should be noted, however, that these potentials are restrictions at the ship hull and at the free surface of the same interior potential



, and are therefore related and cannot be determined by a division of the interior domain because of the global dependence through Laplace’s equation.

3.2.Lagrangianformulation

The integral over the dynamic part of the pressure can be rewritten as:

 dxdy _η −D

∂

t



dz=

∂

t  dxdy _η −D



dz−  F0

φ∂

t

η

dxdy−  S

φ

˙ xS.ndS.

Since x˙S =X˙ +r× ˙

θ

, the normal velocity of the ship can be expressed with a 6-dimensional velocity vector

ξ

˙ as:

˙

Therefore the force integral _S

φ

x˙S.nd S is linear in the velocity vector

ξ

˙:



S

φ

˙

xS.ndS=

ξ

˙·

γ

, (14)

with

γ

the six-dimensional momentum vector defined as:

γ

= 

S

φν

dS. (15)

The potential energy of the water is the sum of the free surface contribution:

PF

(

η

)

= 1 2g  F0

η

2_{dxdy, (16)}

and the part at the wetted ship hull:

PS

(

ξ

)

= 1 2g  S0

ζ

2_{dxdy. (17)}

The dynamics of the rigid body without water interaction is determined by the body Lagrangian in the gravitational field as:

Lb=

1

2M

ξ

˙· ˙

ξ

− gM

ξ

· e3. (18)

Here M is the mass–inertia matrix; the diagonal elements are given by diag

(

M

)

=

(

M,M,M,I1,I2,I3

)

with M the total mass of the body, and Ikthe moments of inertia around the three coordinate axes.

Adding the canonical action of the free water part and the ship Lagrangian produces the total Lagrangian density:

LSF = 1 2M

ξ

˙· ˙

ξ

+

ξ

˙·

γ

− PSF

(

ξ

,

φ

S

)

+  F0

φ

F

∂

t

η

dxdy− HF

(

φ

F,

η

)

. (19) Here HF is the free surface Hamiltonian:

HF

(

φ

F,

η

)

=KF

(

φ

F,

η

)

+PF

(

η

)

, (20)

and PSF(

ξ

,

φ

S) is the total ship potential energy including the kinetic energy of the fluid under the ship:

PSF

(

ξ

,

φ

S

)

=gM

ξ

· e3+KS

(

φ

,

ζ

)

+PS

(

ζ

)

. (21) The wave–ship Lagrangian (19)will be used to derive the wave–ship Hamiltonian formulation.

3.3. Hamiltonianformulation

To arrive at the first order in time Hamiltonian dynamics, the momentum variable



canonically conjugate to

ξ

is introduced as



₌

∂

_ξ˙L (see Van Daalen [7]); explicitly:



=

∂

_ξ˙L=M

ξ

˙+

γ

, (22)

is the sum of the rigid body momentum and the fluid momentum

γ

. The total Lagrangian can then be rewritten by eliminating

ξ

˙ in favor of



:

˙

ξ

=M−1

(



−

γ

)

, (23)

to obtain the canonical action density:

L=



· ˙

ξ

− HS

(



,

ξ,

φ

S

)

+



F0

φ∂

t

η

dxdy− HF

(

φ

,

η

)

, (24)

with a resulting ship Hamiltonian HS:

HS

(



,

ξ

,

φ

S

)

=

1 2M

−1

₍

_

₋

_γ

₎

_·

₍

_

₋

_γ

₎

_+P

SF

(

ξ

,

φ

S

)

. (25)

In the canonical variables (

η

F,

φ

,

ξ

,



) the total Hamiltonian becomes:

H

(

η

,

φ

,

ξ

,

)

=HS

(



,

ξ

,

φ

S

)

+HF

(

φ

,

η

)

. (26) From (24) all Hamilton equations follow. As in the previous section, the free surface equations are obtained from variations with respect to the free surface elevation

η

F and the potential at the free surface

φ

F:

∂

t



η

F

φ

F



=



δ

φFHF

(

φ

,

η

)

−

δ

ηHF

(

φ

,

η

)



=



∂

N



inF0 −

g

η

+

δ

ηKF

(

φ

,

η

)



in F0



. (27)

The coupled ship fluid interaction is given by the Hamilton equations obtained from variations with respect to the canonical variables

ξ

,



leading to:

∂

t



ξ





=



∂

HS

(



,

ξ

,

φ

S

)

−

∂

ξHS

(



,

ξ

,

φ

S

)



, (28)

and written in full as:

∂

t

ξ

=M−1

(



−

γ

)

, (29)

∂

t



=−

(

∂

_ξPSF

ξ

,

φ

S

)



. (30)

Eq.(29)is a reformulation of the definition of the canonical momentum



. Eq.(30)needs some more work to recover the usual description with the appearance of the pressure as driving force for the ship acceleration.

First observe the variation of the potential energy for a variation

δξ

:

δ

PS

(

ξ

)

=gM

δξ

· e3+g 

S0

ζ

.

δξ

·

ν

dxdy,

so that the derivative of the potential energy term becomes:

∂

ξPS

(

ξ

)

=gMe3+g 

S0

ζν

dxdy.

For the variation of the kinetic energy with

δξ

a similar result as for the free surface is obtained (keeping the surface value of the potential the same):

δ

KS

(

φ

,

ξ

)

=  S



₁ 2

|

∇|

2 −

∂

n



.

∂

z





δξ

·

ν

dS, so that:

∂

ξKS

φ

,

ξ



=  S



₁ 2

|

∇|

2₋

_∂

n



.

∂

z





ν

dS.

Expanding the change of canonical momentum leads to:

∂

t



=M

ξ

¨+

∂

t  S

φν

dS=M

ξ

¨+  S

(

∂

t

φ

)

ν

dS.

Taking all contributions to the momentum equation together, the result can be written as:

M

ξ

¨=−gMe3+  S



−

∂

t

φ

− 1 2

|

∇|

2 − g

ζ

+

∂

n



.

∂

z





ν

dS, or as: M

ξ

¨=−gMe3+  S p

ν

dS, (31)

with the pressure:

p=



−

∂

t



− 1 2

|

∇|

2 − gz



z=ζ.

This is the usual expression with the explicit appearance of the integrated pressure and justifies the ship momentum equation in the Hamiltonian formulation.

Finally, variations of the Lagrangian with respect to the potential at the ship hull

φ

Sthat appears in the fluid momentum

γ

and in PSF(

ξ

,

φ

S) leads to:

∂

n



=M−1

(



−

γ

)

·

ν

atS,

which, using the dynamic Eq.(23)for

ξ

˙ produces the compatibility between the fluid velocity under the ship and the ship velocity:

∂

n



=

ξ

˙·

ν

atS. (32)

4. Hamiltonianevolution

The Hamiltonian system of the fully nonlinear dynamic wave–ship interaction in fluid variables that are defined on the free surface only, constitute a very compact description for which it is, however, not obvious how to perform a time stepping procedure. The main aim of this section is to define a time stepping procedure for which it will be required for two different mixed boundary-value problems to transfer the Dirichlet data to Neumann data and vice versa. Briefly stated, the procedure will be as follows.

At the start of a new update in an explicit time stepping procedure, the new ship position will be known, determined by

ξ

; hence the projections of the new free surface and ship area can be updated. Since the ship velocity

ξ

˙ is not known, nor the potential or its normal derivative at the ship hull, the compatibility condition at the ship hull

∂

n



=

ξ

˙ ·

ν

at S cannot be used immediately. In order to proceed, we use the splitting of the interior potential as the sum of two parts as described by Cummins [13]and Ogilvie [14]. The so-called fluid potential



_fl will be required to have the Dirichlet value at the free surface, and be tangential along the wetted part of the ship (usually associated with solving the diffraction problem). Then the fluid momentum can be calculated. After having calculated the added mass matrix for the ship in the new position as the values of the normalized radiation potentials, the dynamic equation for

ξ

˙ can be reformulated to provide, using the known updated value of the canonical ship momentum



, the desired ship velocity

ξ

˙ at the new time, after which the instantaneous radiation potential



rad can be calculated, and the updated total potential is known. Then enough information is available to determine the right-hand-sides of all the dynamic equations.

In the first section we describe the splitting of the total potential in an instantaneous radiation potential and a fluid potential by formulating the mixed boundary value problems that have to be satisfied at each instant. It is important to note that in a time dynamic code, the evolution (‘memory’) of each instantaneous ship impulse is not calculated separately, but becomes part of the fluid potential that will furthermore include incoming and diffracted waves. Yet at each time step the instantaneous radiation problem and the diffraction problem will have to be distinguished, which is consistent with the order in which forces act. The time stepping procedure is described in detail in the second section. In the last section the Hamiltonian ship dynamics is linked to the often used forced damped harmonic oscillator model. It will also be shown that the kernel of the convolution to calculate the history of successive instantaneous impulses can be recovered and linked to the frequency dependent added mass and damping coefficients.

4.1. Boundaryvalueproblemsforpotentials

In order to be able to describe the time stepping algorithm in the next section, the underlying fluid motion, with the mixed boundary value problems for the various distinguishable fluid potentials, will now be detailed.

From the Hamilton equations it follows that at each time the total potential



tot, satisfying the Laplace equation and the bottom and lateral conditions, has to satisfy mixed boundary conditions at the surface: a Dirichlet condition at the free surface and a Neumann condition at the ship:



tot=

φ

atF

∂

n



tot=

ξ

˙·

ν

atS.

(33)

The total potential



will be split in two contributions (see Cummins[13]and Ogilvie[14]), the fluid part



flof the potential and the radiation potential



_rad:



tot=



f l+



rad . (34)

Note that in the time stepping of the fully nonlinear wave–ship problem with updated values of the canonical variables, and the ship in its updated position, the calculation of the total potential is a genuine linear superposition of the constituent fluid and the radiation potential: the linearity of the boundary value problem makes the splitting possible and exact.

The fluid potential



flis defined to take care of the potential at the free surface, and to have a tangentialflow under the ship. Hence it satisfies the mixed Dirichlet–Neumann boundary value problem:

∂

n



f l=0atS



f l=

φ

atF. (35)

The vanishing Neumann condition at the ship expresses that



fl does not directly affects the ship’s motion. The fluid po- tential determines part of the wave momentum vector

γ

₌_S

φν

d S that will be denoted by:

β

= 

S



f l

ν

dS. (36)

Note that if, as in boundary integral methods, the potential of incoming waves is given in the absence of the ship, this ‘incoming’ potential



0has to be compensated by a diffraction potential



diff (usually denoted by



7) that satisfies:

∂

n



di f f =−

∂

n

0

atS

so that the sum



0 +



di f f satisfies (35); the evolution of previous instantaneous radiation potentials has then to be calculated separately in BI panel methods (see Section 4.3). The fluid potential



fl is different since now the sum of



0 (in the free surface area) and the history of previous instantaneous radiation potentials is taken as ‘incoming’ wave; the diffraction potential is added to get tangential flow under the ship.

The instantaneous radiation potential (instantaneous impulse potential) has to satisfy:

∂

n



rad=

ξ

˙·

ν

atS



rad=0atF,

(37)

so that



rad is only the direct effect of the instantaneous impulse

ξ

˙ ·

ν

and does not include the effect of impulses at previous times. The instantaneous radiation potential is linear in the velocity vector

ξ

˙ . Since in the time stepping procedure that will be described in the next section,

ξ

˙ is not known, the instantaneous radiation potential has to be split in normalized instantaneous potentials for each component of

ν

. Then the linear dependence on

ξ

˙ makes it possible to write the potential



rad with the vector of normalized instantaneous potentials

(1) as:



rad=

ξ

˙·

(1), with

∂

n

(1)=

ν

. (38)

Hence, each component



(_j1) satisfies a mixed Dirichlet–Neumann problem with

ν

jas the specific Neumann value at the ship boundary:

∂

n



(_j1)=

ν

jatS



(1) j =0atF. (39)

Having obtained the fluid potential and the normalized instantaneous radiation potential vector, the total momentum vector is given by:

γ

=  S



tot

ν

dS=  S



w

ν

dS+  S



rad

ν

dS=

β

+Ma

_ξ

_{˙. (40)}

Here the so-called added mass–inertia matrix Ma _{is introduced that is determined by the normalized instantaneous} radiation potentials as:

Ma=  S

ν

(

(1)

₎

T_{dS, (41)} i.e., Majk=  S

ν

j



(1) k dS=  S



(1) k

∂

n



(1) j dS.

Note that Ma_jk are genuine constants that may depend on time if in a fully nonlinear treatment the normalized instan- taneous radiation potentials are calculated at each updated ship position. But in a linear theory, where all potentials are calculated at the equilibrium position of the ship, Ma_jk are genuine constants, not depending on frequency, different from added mass quantities that take into account the history of the radiation potentials in time domain–frequency coupled methods, see Section4.3.

A final observation is that the splitting of the total potential leads to a splitting of the kinetic energy in the following way. From Green ’s identity:



∇



f l·

∇



raddxdydz=



F0∪S0



rad

∂

n



f ldxdy,

and the facts that

∂

n



f l = 0 at S and



rad =0 at F, the fluid potential and the instantaneous radiation potential are ‘orthogonal’:



∇



f l·

∇



raddxdydz=0, which holds, in fact, for each component:



∇



f l·

∇



(j1)dxdydz=0. (42)

Hence the kinetic fluid energy is the sum of the separate energies:



|

∇

tot

|

2dxdydz= 

∇

f l



2 dxdydz+ 

|

∇

rad

|

2dxdydz, and the kinetic energy over the ship area can be written as:

KS

(

φ

,

ζ

)

=KS

(

φ

f l,

ζ

)

+

1 2M

4.2. Timestepping

The structure of the equations in Hamiltonian form determine the fully nonlinear time marching. The task that has to be executed at each time step is to calculate the right-hand side (the input for the next update) from the data at that time step. The data are the values of the state variable, the canonical variables in the Hamiltonian formulation:

η

F,

φ

F,

ξ

and



. The calculation of the right-hand sides should be done in a specific order, which is related to the order in which forces work.

First the ship is repositioned to the new position and orientation according to the vector

ξ

. Together with

η

F this defines the boundaries of the fluid in which the surface potentials have to be calculated.

Given the interior domain and ship positioning, the fluid potential and the normalized instantaneous radiation potentials can be calculated by solving the mixed Dirichlet–Neumann problems. With this information and the available value of the total momentum



, the ship velocity

ξ

˙ can be determined. The total potential is then obtained as the sum of the fluid potential and the actual instantaneous radiation potential. Then sufficient information is available to calculate all terms in the right-hand sides of the Hamilton equations.

In more detail, in each time step, the following program has to be executed.

1. For given

ξ

determine the new position and orientation of the ship, the free surface F and the wetted part of the ship hull S. For the updated ship position, the hydrostatic contribution

∂

_ξPSFin the right-hand side of

∂

t



can be calculated. 2. Determine, with the ship in the new position, the fluid part of the potential,



fl, calculate the fluid momentum vector

β

and determine the normalized instantaneous radiation potentials

(1)_{, and the added mass matrix M}a_.

3. The total potential is given by the sum of the fluid potential and the instantaneous radiation potential



tot =



f l+



rad with



rad =

ξ

˙ ·



(1). Then with ( 40) the actual value of

ξ

˙ can be determined by using the known value of



from



=M

ξ

˙₊

γ

₌

(

M₊Ma

)

ξ

˙₊

β

as:

˙

ξ

=

(

M+Ma

)

−1

(



−

β

)

. (44)

This calculated value needs to be used for the right hand side in the equation of the canonical ship position; although equivalent to the Hamiltonian formulation

ξ

˙ ₌M−1

(



₋

γ

)

it had to be rephrased to make it possible to calculate

ξ

˙ without the knowledge of

γ

in the time stepping.

4. Knowing

ξ

˙, the instantaneous radiation potential can be calculated, and the updated total potential



totis determined. Then all remaining parts of the right-hand sides in the Hamilton equations can be calculated, which completes the Hamiltonian evolution procedure.

4.3. Shipdynamics

In this section we will give an interpretation of the Hamiltonian ship dynamics and relate it to other well known formulations in the literature. The change of the canonical momentum will be written for brevity as:

∂

t



=

τ

H,

where

τ

H _{are the hydrostatic force terms. With}

_

_=M

_ξ

˙₊

_γ

_{this is:}

∂

t

(

M

ξ

˙+

γ

)

=

τ

H,

which is mathematically equivalent to:

∂

t

(

M

ξ

˙

)

=−

∂

t

γ

+

τ

H.

In this last non-Hamiltonian formulation, the dynamic pressure term

∂

t

γ

at the right hand side is treated as a force in a numerical time stepping procedure. This is cumbersome, since then the dynamic part of the pressure has to be evaluated at the same instant that the potential is updated, requiring a backward differentiation in time or solving a separate dynamic equation for the pressure. In the Hamiltonian formulation the time stepping algorithm as described above can be used without the need to evaluate the dynamic part of the pressure

The fully coupled Hamiltonian dynamics has clear advantages compared to frequency Boundary Integral panel methods. In the dynamic formulation the splitting of the total potential is given by



tot =



f l+



rad, with



rad the instantaneous radiation potential. With the corresponding splitting of the total wave momentum as

γ

=

β

+Ma

ξ

˙ we can write the dynamics as:

∂

t

((

M+Ma

)

ξ

˙

)

=− ˙

β

+

τ

H. (45)

The momentum part of the instantaneous radiation is given by

∂

t

(

Ma

ξ

˙

)

, where the added mass matrix Ma is ‘constant’, possibly depending on time if the changing position and orientation of the ship are taken into account.

In the dynamic code the term

β

˙ is the dynamic fluid force caused by the fluid potential



fl which is determined by the updated surface potential and the scattered waves. Hence the history of the previous instantaneous radiation potentials has been assembled in the updated potential



fl. In particular, when the ship is forced to move (oscillate) in otherwise

still water, the fluid potential



_flat a certain instant consists only of the history of previous instantaneous radiation pulses. Denoted by



R, this potential caused by the ship motion can be written at time t with a convolutions as [14]:



R

(

x,t

)

=

ξ

˙

(

t

)

·

(1)

(

x,t

)

+

t

−∞

˙

ξ

(

σ

)

·

χ

(

x,t−

σ

)

d

σ

.

Here it is assumed that the radiation potential is calculated with the ship in a fixed position, so that the impulse response function

χ

( x,

τ

) only depends on the time difference

τ

= t −

σ

.

The corresponding momentum

γ

is obtained by integrating over the ship area for each component of

ν

:

γ

= 



R

(

x,t

)

ν

dS=

ξ

˙

(

t

)

· 

(1)

₍

_x

₎

_ν

_dS+ t 0 ˙

ξ

(

σ

)

· 

χ

(

x,t−

σ

)

ν

dSd

σ

, (46) i.e.,

γ

=Ma

ξ

˙+

β

, with,

β

= t 0 ˙

ξ

(

σ

)

· P

(

t−

σ

)

d

σ

andP

(

t

)

= 

χ

(

x,t

)

ν

dS.

Note that Ma

ξ

˙ is the instantaneous momentum, while the convolution

β

describes the history of all the previous instanta- neous impulses. When so desired, for instance for use in frequency domain methods, the function P can be obtained with a dynamic code by taking a unit pulse in a certain direction μ so

ξ

˙₌

δ

_Dir

(

t

)

μ

which leads to:

γ

=

δ

Dir

(

t

)

Ma

μ

+P

(

t

)

μ

.

In frequency domain Boundary Integral panel methods the convolution has to be calculated at each time step, which may be cumbersome [27]. Then the convolution is avoided by adjusting the dynamic equation with frequency depending coefficients, leading to damped mass-spring models for linearized ship dynamics. To illustrate this, taking the temporal Fourier transform of (46), the convolution becomes a product and there results for each frequency:

ˆ

γ

=Ma

ξ

˙

(

ω

)

+

_β

ˆ _with

_β

ˆ

₍

_ω

₎

₌

_ξ

˙

₍

_ω

₎

_{· ˆP}

₍

_ω

₎

_.

Writing the real and imaginary part of Pˆ

(

ω

)

explicitly:

ˆ

P

(

ω

)

=a

(

ω

)

+ib

(

ω

)

,

insertion in the Fourier transformed dynamic equation there results:

−ω2

₍

_M+Ma_+Pˆ

₍

_ω

₎₎



_ξ

₍

_ω

₎

_=−τH k.

The equation for a single harmonic mode can be rewritten in the usual form of an oscillator equation like:

M

ξ

¨+A

(

ω

)



ξ

¨+B

(

ω

)



ξ

˙

(

ω

)

=−τH k,

with A and B the frequency dependent response functions, the added mass and damping coefficient respectively:

A

(

ω

)

=Ma+a

(

ω

)

,B

(

ω

)

=

ω

b

(

ω

)

.

Because the function P( t) above is a caustic function, i.e. vanishing for t< 0, A(

ω

) and B(

ω

) are related by the Kramers– Kronigrelation, given explicitly as:

A

(

ω

)

=Ma+a

(

ω

)

=Ma+

_π

2 _∞ 0 B

(

ν

)

ω

2₋

ν

2d

ν

B

(

ω

)

= 2

ω

2

π

 _∞ 0 a

(

ν

)

ω

2₋

ν

2d

ν

. 5. Conclusions

Starting with an action principle with the integrated Bernoulli pressure as Lagrangian, it was shown how to obtain in a consistent way the extension of the known Hamiltonian formulation for free surface waves to a Hamiltonian formulation of the fully nonlinear wave–ship interaction, while retaining the dimension reduction of dealing only with the fluid potential at the free surface and at the wetted ship hull. A nontrivial explicit time stepping procedure was described that is consistent with the order in which forces work. In particular, a main observation made in this paper is that the Hamiltonian formulation avoids the calculation of the dynamic part of the pressure which becomes part of a canonical momentum variable that is updated by the hydrostatic forces.

In the time stepping procedure, mixed Dirichlet–Neumann boundary value problems have to be dealt with. The data are prescribed at the free surface and at the ship hull which have to be inverted from Dirichlet to Neumann and reversely. This makes it possible to use the coupled wave–ship interaction formulation in (Boussinesq or non-Boussinesq) wave models that can perform this task by extending the dynamics with the 6-degrees of freedom for the ship motion.

In successive publications we will describe and show the performance of fully Boussinesq implementations of the wave–ship interaction by extending free surface codes developed earlier, [20,25,28,29], as part of HAWASSI software [15].

Acknowledgment

Comments of Dr. Ed van Daalen (MARIN) that improved the presentation in a preliminary version of the manuscript are gladly acknowledged. Thanks to anonymous referees whose comments further improved the presentation.

References

[1] J.C. Luke , A variational principle for a fluid with a free surface, J. Fluid Mech. 27 (1967) 395–397 .

[2] V.E. Zakharov , Stability of periodic waves of finite amplitude on the surface of a deep fluid, J. Appl. Mech. Tech. Phys. 9 (2) (1968) 190–194 . [3] L.J.F. Broer , On the Hamiltonian theory of surface waves, Appl. Sci. Res. 29 (1974) 430–446 .

[4] J.W. Miles , On Hamilton’s principle for surface waves, J. Fluid Mech. 83 (1977) 153–158 . [5] H. Lamb , Hydrodynamics, Cambridge University Press, 1932 .

[6] T. Miloh , Hamilton’s principle, Lagrange’s method and ship motion theory, J. Ship Res. 28 (1984) 229–237 .

[7] E.F.G. Van Daalen , Numerical and Theoretical Studies of Water Waves and Floating Bodies, University of Twente, 1993 Ph.D. thesis) . ISBN 90–9005656-4 [8] J.T. Xing , W.G. Price , The theory of non-linear elastic ship-water interaction dynamics, J. Sound Vib. 230 (20 0 0) 877–914 .

[9] R.-Q. Lin , W. Kuang , Nonlinear ship-wave interaction model, part 2; ship boundary condition, J. Ship Res. 2 (2006) 181–186 .

[10] R.-Q. Lin , W. Kuang , A fully nonlinear, dynamically consistent numerical model for solid-body ship motion. i. ship motion with fixed heading, Proc. R. Soc. A 2010 (0310) (2010) 1–17 .

[11] R.-Q. Lin , W. Kuang , A fully nonlinear, dynamically consistent numerical model for ship maneuvering in a seaway, Model. Simul.Eng. 2011 (2011) 1–10 . ID 356741

[12] E.F.G. Van Daalen , E. van Groesen , P.J. Zandbergen , A Hamiltonian formulation for nonlinear wave-body interactions, in: Proceedings of the Eight International Workshop on Water Waves and Floating Bodies, 1993 . Canada

[13] W.E. Cummins , The impulse response function and ship motions, Schiffstechnik 9 (1962) 101–109 .

[14] T.F. Ogilvie , Recent progress toward the understanding and prediction of ship motions, in: Proceedings of the Fifth Symposium on Naval Hydrodynam- ics, 1964, pp. 3–79 . Reprinted as David Taylor Model Basin Report No. 4.

[15] HAWASSI software, Variational Boussinesq Model VBM and Analytic Boussinesq AB, see http://hawassi.labmath-indonesia.org .

[16] H. Bateman , Notes on a differential equation which occurs in the two-dimensional motion of a compressible fluid and the associated variational problems, Proc. R. Soc. A 125 (1929) 598–618 .

[17] F. Treves , Basic Linear Partial Differential Equations, Academic Press, 1975 .

[18] G. Klopman , M. Dingemans , E.V. Groesen , Propagation of wave groups over Bathymetry using a variational Boussinesq model, in: Proceedings of the International Workshop on Water Waves and Floating Bodies, Plitvice, Croatia, 2007, pp. 125–128 .

[19] E. Klopmnan G , van Groesen , M.W. Dingemans , A variational approach to Boussinesq modelling of fully non-linear water waves, J. Fluid Mech. 657 (2010) 36–63 .

[20] D. Adytia , E. van Groesen , Optimized variational 1d Boussinesq modelling of coastal waves propagating over a slope, Coast. Eng. 64 (2012) 139–150 . [21] I. Lakhturov , D. Adytia , E. van Groesen , Optimized variational 1d Boussineq modelling for broad-band waves over flat bottom, Wave Motion 49 (2012)

309–322 .

[22] L.J.F. Broer , Approximate equations for long water waves, Appl. Sci. Res. 31 (1975) 377–395 .

[23] L.J.F. Broer , E.W.C. van Groesen , J.M.W. Timmers , Stale model equations for long water waves, Appl. Sci. Res. 32 (1976) 619–636 . [24] W. Craig , C. Sulem , Numerical simulation of gravity waves, J. Comput. Phys. 108 (1993) 73–83 .

[25] R. Kurnia , E. van Groesen , High order Hamiltonian water wave models with wave-breaking mechanism, Coast. Eng. 93 (2014) 55–70 .

[26] R. Kurnia , E. van Groesen , Spatial-spectral Hamiltonian Boussinesq wave simulations, in: Proceedings of the International Conference on Advances in Computational and Experimental Marine Hydrodynamics, vol.2, 2015, pp. 19–24 . ISBN: 978–93-80689–22-7

[27] E. Kristiansen , A. Hjulstad , O. Egeland , State-space representation of radiation forces in time-domain vessel motions, Model. Identif. Control 27 (2006) 23–41 .

[28] Lawrence, D. Adytia, Variational Boussinesq model for dispersive strongly nonlinear waves, in: Proceedings of the International Conference on Offshore Mechanics and Arctic Engineering, OMAE2016, 2016, Busan, Korea.

[29] R. Kurnia , E. van Groesen , et al. , Localization in spatial-spectral method for water wave applications, in: R. Kirby (Ed.), Spectral and High Order Methods for Partial Differential Equations, ICOSAHOM 2014, Lecture Notes in Computational Science and Engineering, vol. 106, Springer, 2015, pp. 305–313 .