A generating and absorbing boundary condition for dispersive waves in detailed simulations of free-surface flow interaction with marine structures

Delft University of Technology

A generating and absorbing boundary condition for dispersive waves in detailed

simulations of free-surface flow interaction with marine structures

Wellens, Peter; Borsboom, Mart

DOI

10.1016/j.compfluid.2019.104387

Publication date

2020

Document Version

Final published version

Published in

Computers and Fluids

Citation (APA)

Wellens, P., & Borsboom, M. (2020). A generating and absorbing boundary condition for dispersive waves

in detailed simulations of free-surface flow interaction with marine structures. Computers and Fluids, 200,

[104387]. https://doi.org/10.1016/j.compfluid.2019.104387

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

Computers

and

Fluids

journalhomepage:www.elsevier.com/locate/compfluid

A

generating

and

absorbing

boundary

condition

for

dispersive

waves

in

detailed

simulations

of

free-surface

flow

interaction

with

marine

structures

Peter

Wellens

a,∗

_,

_Mart

_Borsboom

b

a Delft University of Technology, Mekelweg 2, Delft 2628CD, the Netherlands b Deltares, PO Box 177, Delft 2600MH, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 10 January 2019 Revised 10 November 2019 Accepted 13 November 2019 Available online 14 November 2019

Keywords:

Free surface waves Wave dispersion Boundary condition Volume-of-Fluid (VoF)

a

b

s

t

r

a

c

t

Theboundariesofnumericaldomainsforfree-surfacewavesimulationswithmarinestructuresgenerate spuriouswavereflectionifnospecialmeasuresaretakentopreventit.Thecommonwaytoprevent re-flectionistousedissipationzonesatthecostofincreasedcomputationaleffort.Onmanyoccasions,the sizeofthedissipationareaisconsiderablylargerthantheareaofinterestwherewaveinteractionwith thestructuretakesplace.Ourobjectiveistoderivealocalabsorbingboundaryconditionthathasequal performancetoadissipationzonewithlowercomputationalcost.Theboundaryconditionisdesignedfor irregularfree-surfacewavesimulationsinnumericalmethods thatresolvetheverticaldimension with multiplecells.Itisforarangeofphasevelocities,meaningthatthereflectioncoefficientperwave com-ponentislowerthanachosenvalue,say2%,overarangeofvaluesforthedimensionlesswavenumber

kh.ThisisaccomplishedbyextendingtheSommerfeldboundaryconditionwithanapproximationofthe lineardispersionrelationintermsofkh,incombinationwithverticalderivativesofthesolutionvariables. Forthisarticle,theboundaryconditionisextendedwithanon-zeroright-handsideinordertoprevent wavereflection,while,atthesametime,atthesameboundary,generatingwavesthatpropagateintothe domain.Resultsofirregularwavesimulationsareshowntocorrespondtotheanalyticalreflection coef-ficientforarangeofwavenumbers,and tohavesimilarperformancetoadissipation zoneatalower cost.

© 2019TheAuthor(s).PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Weemployanumericalmethodforsimulationsofextreme, ir-regular free-surface wave interaction with(floating) structures at sea. A major issue is that waves reflect from computational do-main boundarieswhen nospecial measures are takento prevent reflection. Dissipation zonesare often used near domain bound-aries to resolve theissue.Dissipation zonesinduce a decreaseof wave energyso that whenwavesreach theboundary and reflec-tion occurs, the wave height hasreduced to such an extent that spurious reflectiondoesnot interferewiththeprocessesnearthe structureofinterest.

Dissipation zones need to be several typical wave lengths Lp

long tobe effective,whereLp isthe wave lengthassociated with

the peak frequencyof the irregularwave spectrum.Consider, for

∗ _{Corresponding author.}

E-mail address: p.r.wellens@tudelft.nl (P. Wellens).

instance,thesetupofthedomaininRivera-Arreba etal.[22]:for thissituationwithunidirectional incomingwaves,thedissipation zoneinfrontofthestructure(withrespecttothewavedirection) isoneLplongandthedissipationzonebehindthestructuretakes

uptwoLp.Withanareaofinterestthat isoneLp,the dissipation

zonescomprisethelargerpartofthedomainsizeand,hence,the computationaleffort.Inasettingwithmultidirectionalwaves,with dissipation zoneson all sidesofthe area ofinterest, the ratioof areaofinteresttodissipationzoneareabecomesevenworse. Tech-niques such as stretching the grid with increasing distance from theareaofinterestcanbeemployedtomitigatepartofthe com-putationaleffortandmemoryrequirements,butdissipationzones take up a considerable amount of the total computational effort requiredforthesimulationnonetheless.

There are two common alternatives for dissipation zones: 1. Perfectly Matched Layers (PMLs) [1], and 2. local boundary con-ditions.PMLs areimproved dissipation zonesandrequiresmaller zones to be added to the computational domain. We will not consider PMLs, because we wish to investigatea local absorbing

https://doi.org/10.1016/j.compfluid.2019.104387

boundaryconditionthat doesnot requireanyadditionalareas at-tachedtothecomputationaldomain.Withoutanyadditionalareas, localabsorbingboundaryconditionsarethemorecomputationally efficientalternativetodissipationzonesforpreventingwave reflec-tion.Here, a local absorbing boundary conditionfor long-crested irregularwavesisderivedwithsimilarspuriouswavereflectionas atwo-wave-lengths-long dissipation zone at a reduced computa-tionaleffort.

The absorbing boundarycondition isbased on a local bound-aryconditionoperatorincombinationwitha rational approxima-tionofthefreesurfacewave dispersionrelationandwithsecond derivatesofvelocity andpressureinvertical direction, derived in Wellens [26]. This approach extends the rangeof wave numbers overwhichspuriouswave reflectionislowcomparedtothe Som-merfeld boundary condition. The range was further extended in termsofbothdispersion(wavenumber)anddirectionofthewave componentsinDüzetal.[7],by multiplyingtheboundary condi-tionwithanadditionalboundaryconditionoperator.

What the current article adds to the earlier work is that the absorbingboundary conditionisextended withanon-zero right-handside togenerateincomingwaves,while outgoingwavesare absorbedat the same time. Some approximations of the disper-sion relation lead to stability issues. A stability analysis is pre-sentedto prevent instabilities.The absorbing boundary condition canbecombinedwithanynumericalmethodforfree-surfacewave simulationsthatresolvestheverticaldimensionwithmultiplegrid cells. Forthe purposeofwave impacts onmarine structures, the boundaryconditionwasimplementedinComFLOW[18],a numer-ical method for free-surface waves that is based on the Navier-StokesequationswithaVolumeofFluidapproachforthefree su-face.Boththeabsorbingperformanceandthecombinedgenerating andabsorbing performance are verified by comparing the spuri-ousreflectioncoefficientsofirregularwavesimulationstothe the-oreticalreflectioncoefficients. Andfinally,it isdemonstratedthat thegeneratingabsorbingboundarycondition(GABCforshort)not onlyworksintheoreticalwavepropagationsimulations,butcanbe usedasacomputationallyefficientalternativetodissipationzones ina simulationofwave interactionwith-andwaveimpact loads on-astructureinthefreesurface,whichisvalidatedbyan exper-iment.

2. Governingequations

The numerical method we employ for wave interaction with floating structures solves the Navier–Stokes equations for a fluid inthe presence ofa free (liquid) surfacein a cartesian grid do-main,wherethegridlinesarecolinearwiththecomponentsofthe cartesianright-handedaxissystemwithx=

(

x,y,z

)

T_{.TheNavier–}

Stokesequationsimposeconservationofmassandconservationof momentum,respectively:

∇

· u=0 (1)

∂

u

∂

t +

(

u·

∇

)

u+ 1

ρ ∇

p−

∇

·

(

ν∇

u

)

− f=0, (2)

inwhichthefluidvelocity vectorisrepresentedbyu₌

(

u_,

v

,w

)

T_.

Scalarpisthe pressure,

ρ

isthe fluiddensityand

ν

isthe kine-maticviscosity.frepresentsgravity

(

0,0,−g

)

T _actingasan

exter-nalbodyforce.

The indicator function forfree surface, S, isdisplaced in time andspaceby:

DS Dt =

∂

S

∂

t +u·

∇

S=0, (3)

ina manner derived from Hirt andNichols [16], in combination withalocalheightfunctionina3x3x3blockofcells[18].Thefree surfaceisreconstructedusingPLIC[28].

Fig. 1. Axis system at the boundary.

ThesimulationsinSections7and8aresinglephase(waterand void) with a boundary conditionfor the pressure. The following boundaryconditionforthepressureisimposedbymeansof inter-polationtowardsthereconstructedpositionofthefreesurface:

−p+2

μ∂

uns

∂

ns =−p0+

σ κ

,

(4)

wherensisthedirectionnormaltothefree surfaceandp0 isthe atmosphericpressure.Thesurfacetensionisrepresentedby

σ

and thecurvatureofthefreesurfaceby

κ

,whichiscomputedfromthe valuesofthelocalheightfunction.

Near domain boundaries it is convenient to use an axis sys-tem witha component nnormal, anda component

τ

tangential toboundaryB,seeFig.1.Thesidewallsofthedomainarevertical andzisdirectedupwardwithrespecttothen−

τ

-plane.The di-rectionofimpinging wavesontheboundarywillbe expressedby angle

θ

withrespectthenormaldirectionntotheboundary.

Neartheboundary,sufficientlyfarawayfromastructureinthe flow, we assume that theflow isirrotational. Nearthe boundary, wealsoassumethatwecanapplyalocallinearization.Thisallows forthe introduction ofa linearpotential function

, whose spa-tialderivativesyield thevelocity inthedirectionofthederivative

∇

=u.Withthepotential,thecontinuityequationandlinearized momentumequationbecome:

∇

2

_{=0 (5)}

∂

∂

t +

p

ρ

− gz=0. (6)

Theboundaryconditionsatthefreesurfacesimplifyto:

∂ζ

∂

t −

∂

∂

z =0 (7)

p=g

ζ

. (8)

Bothareimposed atthemeanfree-surfacepositionatz=0.

ζ

is thefree-surfacepositionwithrespecttoz₌0.

Thepropagatingperiodicsolutiontothesystemofequationsin

(5)to(8)inafinitewaterdepthh,wherethebottomisconsidered impermeable,i.e.

∂

_/

∂

z₌0atz_=−h,isgivenby:

j=

ζ

a, jg

ω

j cosh

(

kj

(

h+z

))

cosh

(

kjh

)

sin

(ω

jt− kjn

)

. (9)

Here,

ζ

a,jisthewaveamplitude,

ω

jisthewavefrequency,kjisthe

wavenumberandjisanindicatorforthewavecomponentunder consideration.

Eqs. (6)–(8) can be combined with solution (9) to obtain the dispersionrelation:

ω

2

j=kjgtanh

(

kjh

)

, (10)

whichforourpurposesis berewrittenintoa formulationforthe phasevelocitycj: cj=

ω

j kj =



gh



tanh

(

kjh

)

kjh . (11)

By dispersion we mean that each wave frequency has a differ-ent phase velocity.This aspect isparticularly difficult for bound-aryconditionsforirregularfree-surfacewavesthatareconsidered to be composed ofmany wave components witheach their own amplitude,frequencyand,hence,phasevelocity.

The systemof equationsin (5)to (8)can alsobe recombined intoaplanarwaveequationforthepotentialatthemeanfree sur-face:



∂

2

∂

t2− c 2 j

∇

2



j=0 atz=0. (12)

Wave Eq. (12) is the starting point for the derivation of wave boundaryconditionsatdomainboundaries.

3. Background

Fora local generatingabsorbing boundary conditionthree as-pectsarerelevant:(1)wavedirectionality,(2)wavedispersionand (3)incomingandoutgoingwaves,or-inotherwords-generation and absorption of waves. Literature predominantly discusses ab-sorbing boundaryconditions,which arealso callednon-reflecting boundaryconditionsorradiationconditions,withanemphasison wave directionality, see Givoli [9]. It is relevant for the deriva-tionofourboundaryconditiontodiscussallthreeaspects,starting withwavedirectionality.

An example of an absorbing boundary condition for wave

Eq.(12)istheSommerfeldequation:



∂

∂

t +c0

∂

∂

n



=0, (13)

inwhichc0isatypicalvaluethatapproximatestheapparentphase velocity atthe boundary. The apparent phase velocity forwaves with phase velocity c that approach the boundary under angle

θ

withrespect to the boundary’s normal directionis higher and equaltoc/cos(

θ

).Ifc0isthesameastheapparentphasevelocity, thenthereisnoreflection.Themorec0isdifferentfromthe appar-entphase velocityofa wave component,the moreitisreflected. Theamountofreflectionisfoundfromthereflectioncoefficient: R=c0− c/cos

(θ

)

c0+c/cos

(θ

)

, (14)

whichistheratioofamplitudesoftheoutgoingandthereflected wave components.One could say that the Sommerfeld boundary conditionisonlyaccurate,i.e.haslowreflection,forasmallrange ofwavedirectionswherec0isclosetoc/cos(

θ

).Thedevelopment ofboundary conditionsis concernedwithextending therangeof directionsoverwhichlowreflectionoccurs.

There are two families of absorbing boundary conditions for wave directionality. One is the family of boundary conditions where the c/cos(

θ

) behaviour in the apparent phase velocity is approximatedwithderivativesofthesolutionvariablesin tangen-tial direction

τ

along the boundary,see Fig.1. Engquistand Ma-jda[8]suggestedanorderexpansionofc/cos(

θ

),whereincluding moretermsintheexpansionleadstohigherderivativesalongthe boundary anda larger rangeover whichthe absorbing boundary

conditionisaccurate.TheexpansiontoO

(

1

)

leadstothe Sommer-feldequationwithc0=c.TheexpansiontoO

(

2

)

gives:



c

∂

2

∂

t

∂

n−

∂

2

∂

t2+ c2 2

∂

2

∂τ

2



=0. (15)

FromEq.(15)wefindthat higherorderapproximationsofthe dispersion relation give higherderivatives along the boundary in the absorbing boundary condition operator. By means of higher derivativesin theABC, waveswithlarger incomingangles to the boundaryareaccuratelyabsorbed.

Theotherfamilyofboundaryconditionsapproximatesc/cos(

θ

) withderivativesofthesolutionvariablesinnormaldirectiontothe boundary.IntroducedbyHigdon[15],itisgivenby:

J  j=1



∂

∂

t +cj

∂

∂

x



=0. (16)

Wavereflectionbecomessmallerastheorderoftheabsorbing boundaryconditionincreases.Theorderoftheboundarycondition, inthissense,reflects thenumberof productsJused toconstruct theoperator. Higher orderabsorbing boundaryconditions feature higherderivatives.As theorderoftheabsorbingboundary condi-tionincreases, itbecomes increasingly difficultto implement the numerical equivalent of the higher derivatives at the boundary. Givoli[10] reportsthat Higdon operators beyondorder threeare rarelyfoundintheliterature.CollinoandJoly[5]introduced auxil-iaryvariablestoremedythedifficultieswithimplementinghigher derivatives.Instead ofNthderivatives,thenasystemofN+1 ad-ditionalequationsissolvedattheboundary.

Auxiliary variables have been used in Nth order absorbing boundaryconditions derived by Grote andKeller [12],Givoli and Neta [11] andHagstrom and Warburton [14], among others.The Hagstrom-Warburtonformulationoftheauxiliarysystemof recur-siveequationsis:



a0

∂

∂

t +c

∂

∂

n



=a0

∂

∂

t

ξ

1 (17)



aj

_∂

∂

_t +c

_∂

∂

_n



ξ

j=



aj

_∂

∂

_t +c

_∂

∂

_n



ξ

j+1 for j=1,2,...,N, (18)

ξ

N+1=0, (19)

whereaj arecoefficientsand

ξ

jare auxiliaryvariables.The setof

equationsfeaturesnonehigherthanfirstderivativesinnormal di-rection to the boundary. The major difficulty of thisway of for-mulatingboundary conditions are the cornerconditions required neartheends ofdomainboundaries. An overviewis providedby Hagstrometal.[13].We didnot choosethisapproachbecause it proveddifficult tocombineauxiliary variableabsorbing boundary conditionswithincomingwavesoverthesameboundary.

Dispersionof free-surface waves,the fact that each wave fre-quencyhasadifferentphasevelocity,hasnotbeenanexplicit sub-ject in mostof the literature regarding absorbing boundary con-ditions.As wasdiscussed forwave directionality,theSommerfeld equation is accurate ina limitedrange around phase velocity c0. The objective ofa better absorbing boundary conditionis to ex-tend the range of wave frequencies over which the reflection is low.Wefoundtwowaystoaccountfordispersionattheboundary: oneistousealoworderboundaryconditionandcombineitwith anestimatorfortheapparentphasevelocity attheboundary, the secondistouseapproximationofanon-localabsorbingboundary conditionwithverticalderivativesofthesolutionvariables.Asan exampleoftheformer,wementionOrlanski[20].Hesuggeststhe

Fig. 2. The resolved phase velocity c as a function of time, when the Orlanski boundary condition is used. The discontinuities result from the zero crossings of the denominator in the equation for c . The irregular wave signal was generated with a JONSWAP spectrum with T p = 15 s.

useof theSommerfeld condition, witha dynamicapproximation ofthephasevelocity,obtainedfromthesolutionitself:

c0=−

∂

/

∂

t

∂

/

∂

n. (20)

Fig. 2shows arepresentation ofthe phase velocity, whenthe Orlanskiboundaryconditionisappliedtoanirregularwavesignal conformingtoaJONSWAPspectrumwithpeakperiodTp=15s.It

showsrather largediscontinuitiesinthe resolved valueforc, co-incidingwiththe timeswhere the denominatorin (20)becomes zero,whichcanandwillleadtoanunstablesimulation.Thereare stabilizationmethods,suchasweightedmovingaveragefilteringof

c0 in[3],butthesemethodsshow large reflectionforwave com-ponentswithanapparentphasevelocitythatisdifferentfromthe filteredc0.

DgayguiandJoly[6]discussabsorbingboundaryconditionsfor the simulation of free-surface waves in the nz-plane. Herein, an exact,non-localabsorbingboundaryconditionoperatorisderived. Theexact operatoris simplifiedby means ofrational approxima-tions to yield an absorbing boundary condition that is local in space andlocal in time with vertical derivatives of the solution variablesalong theboundary.Numerical resultsare presentedfor simulations of long waves with a zeroth order and a first order boundarycondition.Withthelatterboundarycondition,the reflec-tionwasgenerallysmall;exactfigureswerenotmentioned.

InourstudiesinWellens[26],wefoundthatusingrational ap-proximations ofthe dispersion relation (10) incombination with verticalderivatives ofthesolution variables yields a considerable increaseoftherangeofwavefrequenciesandphasevelocitiesthat canbeabsorbedaccurately.InDüz etal.[7],furtherextension of therangewithlowreflectionwasobtainedby combiningthe ver-ticalderivativesof thesolution variables withderivatives in nor-maldirectiontotheboundaryasinHigdon[15].Itproveddifficult toturnthelatterapproachintoaboundaryconditionforanopen boundary,wherewavesfrominsidethedomainareabsorbedwhile wavesfromoutsidethedomainenteratthesametime.

Anotherwaytoderiveboundaryconditionsisbymeansofthe methodofcharacteristicsdiscussed byBlayoandDebreu[2], Ver-boomandSlob[24] andVan Dongeren andSvendsen [23]. Here, themethodofcharacteristicsisappliedtothewaveequation. Con-sideragainwaveEq.(12),butnowinonedimensioninthe direc-tionnormaltotheboundary:



∂

2

∂

t2− c 2

∂

2

∂

n2



=0. (21)

Thetermontheleft-handsideofEq.(21)canbefactorized:



∂

∂

t − c

∂

∂

n



∂

∂

t +c

∂

∂

n



. (22)

The term



_∂∂_t+c_∂∂_n



in Eq. (22) is constant along lines

dn/dt=−c.Wecallthistermtheincomingcharacteristic. AboundaryconditionforwaveEq.(21)isobtainedbyassigning a value tothe incomingcharacteristic. Whenno wavesenter the computational domain through the boundary, i.e. by assigning a valueof0totheincomingcharacteristic,theSommerfeldequation isobtained.Whentheincomingcharacteristicreceivesanon-zero time-dependent right-hand side value, we are sending in waves over the boundary, while outgoing wavesmay leave the compu-tationaldomainunaffected.Inthisway,anopen,orgenerating ab-sorbingboundaryisobtained.

AccordingtoCarpenter[4]andPerkinsetal.[21],the Sommer-feldequation withnon-zeroright-hand side value should be for-mulatedasfollows:



∂

∂

t +c

∂

∂

n



=



∂

∂

t +c

∂

∂

n



in_{. (23)}

Below, a generating absorbing boundary condition for open boundariesisderivedandevaluated,whichcombinesthemethod ofcharacteristics witha rational approximation ofthe dispersion relation andverticalderivatives ofthe solution variables.The ra-tional approximation and the vertical derivatives of the solution variablesgiveanincreasedrangeofwavedirectionsandphase ve-locitiesforwhichreflectionislow.Thesameoperatorusedonthe left-handsideoftheequation isapplied totheright-handsideof theequation that prescribes theincomingcharacteristic,allowing wavestoenterthedomainwhileoutgoingwavesareabsorbed. 4. Derivation

To arrive at an absorbing boundary conditionfor free-surface wave simulations, we start from Sommerfeld condition (13) and lookforways toincludethe effectofwave dispersion inthe ap-plied c0. As explained in the previous section, application of a solution-dependentvalue(20)isnotgeneralenough.Ontheother hand,wewanttoavoidtheuseofmultipleoperatorsasin(16)in ordernot to complicatethe numericalimplementation with one-sidedhigher derivativesin normaldirection tothe boundary. In-spiredby DgayguiandJoly[6],weconsidertheFouriertransform ofthe Sommerfeldequation andreplace c0 bythe expressionfor thelineardispersivephase velocity(11),approximatedbya ratio-nal[2,2]Padé approximation.

WestartwiththeFouriertransformof(13):



−

ω

j+c0kj



φ

j=0for j=1,2,· · · ,N. (24)

Inourinterpretationf1:kh→c0isaconstantfunctionthat ap-proximates dispersion relation(10) over a rangeof phase veloci-ties,orequivalently,arangeofvaluesforthedimensionlesswave number kh.c0 intersects thedispersion relationat one value for

kh.Here,thedifferencebetweenthefunctionsiszeroand,hence, thereflectioncoefficientiszero,seeEq.(14).

Themain thoughtbehindourboundaryconditionisthatc0 is replaced by a better approximation of the dispersion relation. A rationalpolynomialinkhthatisabetterapproximationtothe dis-persionrelationthanf1 isgivenby

f2:kh→



gha0+a1

(

kh

)

2 1+b1

(

kh

)

2 . (25)

The coefficientsa0,a1 andb1 canbe chosen such that differ-ent kh-ranges of the dispersion relation are approximated well. With a0=1.000, a1=0.150 andb1=0.367,the dispersion rela-tionisapproximatednearkh=0withsixth-orderaccuracy,ascan be shown by means of a series expansion. This is illustrated in

Fig.3.It showsthedispersion relation togetherwiththe approx-imating functions f1 and f2 (with the coefficients stated above),

Fig. 3. Approximation of the dispersion relation in (a). In (b), the reflection coeffi- cient associated with these approximations.

anditalsoshowsthereflectioncoefficientforbothapproximations ofthedispersionrelation.Thereflectioncoefficientforf2 islower thanforf1overalargekh-range.Usingf2inaboundarycondition, therefore,makesitabetterboundaryconditionforirregularwaves withwave componentsin thekh-rangewhere theapproximation is sufficiently accurate. It isa strength of the approach that one can choose the range over which the approximation is accurate depending on what range of wave numbers is simulated. There isa trade-off betweenaccuracyandrange.Whentherangeis in-creased,alsothereflectioncoefficient increases.Fig.3 alsoshows function f2 withthe coefficients a0=1.040, a1=0.106 andb1= 0.289.Withthesecoefficients,thereflectioncoefficientisnon-zero nearkh=0,butitislessthan2%intherangebetweenkh=0and

kh₌6_,which isacceptable inmanyapplications.Finding the co-efficients fora specific simulation wasautomated by meansof a bruteforcealgorithmthateithergivesthelowestreflection coeffi-cientforagivenrange,orthelowestpossibleareaofthereflection spectrumforagivenwavespectrum.

When c0=f2 is substituted into boundarycondition (24), we obtain:

−

1+b1



kjh



2

ω

j+



gh

a0+a1



kjh



2

kj

φ

j=0 for j=1,2,· · · ,N. (26)

NotethatEq.(26)ismultipliedby

1+b1



kjh



2

.

Boundary condition (26) cannot be transformed back to the timedomain,becauseitisnonlinearinkj.k2j willneedtobe

elim-inatedfromtheboundarycondition.

Thewavenumberkjcanbe foundbytakingderivativesofthe

solutioninspace.Considerthesolutionofthelinearsystematthe boundary in Eq.(9). We find that second derivative of the solu-tioninverticaldirectionisequaltothesquareofthewavenumber timesthesolutionitself

∂

2

φ

j

∂

z2 =k

2

j

φ

j. (27)

Eq.(27)issubstitutedinboundarycondition(26)toyield



−



1+b1h2∂ 2 ∂z2



ω

j+



gh



a0+a1h2∂ 2 ∂z2



kj

φ

j=0for j=1,2,· · · ,N, (28)

andcanbetransformed backfromFourierspace.Inthetime do-main,theabsorbingboundaryconditionbecomes



1+b1h2

∂

2

∂

z2



∂

∂

t +



gh



a0+a1h2

∂

2

∂

z2



∂

∂

n



=0. (29)

It gives little reflection for wave components in the range of

kh-valueswhererationalapproximation(25)ofthedispersion rela-tionisaccurate.Todemonstratehowtherangeoflowreflectionof ourboundaryconditionhasincreased,thereflectioncoefficient is comparedtothe reflectioncoefficient fortheSommerfeld bound-arycondition.Thecoefficientsinrationalapproximation(25)were giventhevaluesa0=1.040,a1=0.106,b1=0.289.Forthe coeffi-cientintheSommerfeldequationsa valueofc0=0.6



ghis cho-sen.

Fig.4showsthe reflectioncoefficient fora rangeofkh-values andanglesofincidence

θ

forourboundaryconditionandthe Som-merfeldcondition.Thearea wherethereflectioncoefficientis be-lowfivepercenthasbeenenclosedbyaboldline;fivepercent re-flectionisconsideredtolerableinpracticalwavesimulations.From

Fig.4, we findthat the area withless than5% reflection has in-creasedconsiderably.

Theleft-handsideof(29)canbeviewedastheincoming char-acteristicofthewaveequation,wheretheeffectofwavedispersion has been approximated using (25). Following Carpenter [4] and Perkinsetal.[21],theincomingcharacteristiciscombinedwitha non-zeroright-handsideconsistingofthesameoperatorasonthe left-handside appliedtoa knownincomingwavepotential. With anon-zeroright-handside,absorbingboundarycondition(29) be-comes:



1+b1h2

∂

2

∂

z2



∂

∂

t +



gh



a0+a1h2

∂

2

∂

z2



∂

∂

x



=Rin_{, (30)} inwhich: Rin₌



1+b1h2

∂

2

∂

z2



∂

∂

t +



gh



a0+a1h2

∂

2

∂

z2



∂

∂

x



in_.

We call thisboundary conditiona GABC, a Generating Reflecting BoundaryCondition.

5. Stability

Therangeof kh-valuesforwhichboundary condition(30) ab-sorbsoutgoingwaveswithlittlereflectionistherangeofkh-values forwhichthe appliedunderlying dispersion approximationis ac-curate. That range can be tuned to a particular application by a suitablechoiceoftheparameters.Weoptimizedthesecoefficients for a reflection of less than 2% from kh=0 to the largest pos-sible value for kh. However, we found that this set of parame-ters(a0=1.050,a1=0.100,b1=0.310)cannotbeappliedbecause

Fig. 4. Reflection coefficient in percentages as a function of incidence angle θ and dimensionless wave number kh . The Sommerfeld condition in (a) is compared to the boundary condition with a rational approximation in (b). In the boundary conditions, c 0 = 2 / 3 gh and a 0 = 1 . 040 , a 1 = 0 . 106 , b1 = 0 . 289 were used, respec- tively.

simulationsbecameunstable.Itappearsthata0,a1,b1 haveto sat-isfystabilityconstraints.

Sofar,wehaveonlyconsideredpropagatingwavemodes,which are oscillatory in time, oscillatory in horizontalspace, and expo-nential in the vertical, see (9). Other wave modes that are os-cillatory in time are the evanescent waves, which are exponen-tial in horizontal space and oscillatory in the vertical. It can be shown that they have no effect on the absorption properties of

(30),whichiswhywe havenot consideredthem.However,it ap-pearsthatEqs.(5)–(8)supportyetanothertypeofsolution,which we refer to as spurious modes. These are modes that, like the evanescentmodes,areexponentialinhorizontalspaceand oscilla-toryinthevertical.Butratherthanbeingoscillatory intime,they areexponentialintime:

s=

ζ

a g

ω

s cos

(

ks

(

h+z

))

cos

(

ksh

)

exp

(ω

st− ksn

)

, (31)

Here, subscript s is used to indicate that it concerns a spurious wavemodethatgrows exponentiallyintime. Whentheexponent

Fig. 5. Roots of the reflection coefficient’s denominator.

is positive, the mode is unstable.With boundary conditions that arenormallyapplied, spuriousmodescannot exist,i.e.stabilityis alwayssatisfied.Withourboundarycondition,however,itappears that certain combinations of a0,a1 andb1 leave spurious modes undefined,asaresultofwhichasimulationwillbecomeunstable. Whenanincomingandoutgoingsolutionmodearesubstituted intoboundarycondition(29),itispossibletoformulateareflection coefficientastheratiooftheiramplitudes.Forstablesolutions,the reflectioncoefficientnecessarilyhastobesmallerorequalthan1:

Rsp= a0−a1(kh)2 1−b1(kh)2 −

tan(ksh) ksh a0−a1(kh)2 1−b1(kh)2 +

tan(ksh) ksh ≤ 1 (32)

In(32),amechanismisrecognizedby whichthereflection co-efficientcanbecomelargerthan1.Eq.(32)hassingularitieswhere thedenominatorhasroots.Thefollowingfunctionsarethe expres-sionsfromthedenominatorin(32).

f1:ksh→

tan(ksh) ksh f2:ksh→a0− a1

(

ksh

)

2 f3:ksh→1− b1

(

ksh

)

2 (33)

Functionsf1andf2in(33)areplottedinFig.5.Thefigureshows only the real values of f1. It has imaginary values in the range

kh∈[

π

/2,

π

].The rangeofimaginaryvaluesgivesusthe oppor-tunitytoensure stability.Iftherootsofthefunctionsf2 andf3 in Eq. (33)are chosen in the rangewhere f1 has imaginaryvalues, then instability cannot occur. Thisputs restrictions on thevalues for the coefficients a1 and b1. In Fig. 5, the function f2 is plot-ted with a0=1 andtwo different valuesfor a1, a1=a0/

π

2 and

a1=a0/4

π

2.Withthesevaluesfora1,therootsoff2areprecisely ontheouterlimitsoftherangekh∈[

π

/2,

π

].Therefore,thefirst requirementforstabilityisthata1hasavaluebetweena0/

π

2 and 4a0/

π

2.

The same line of reasoning applies to the function f3 in Eq.(33)andthevalueforb1.Thesecondrequirementforstability, then,isthatthevalueforb1ischosenbetween1/

π

2 and4/

π

2.

The third and final requirement for stability comes from the limitbehaviour ofthereflectioncoefficientinEq.(32).Ifwe con-siderthelimitkh→∞ofthereflectioncoefficient,itisfoundthat

b1shouldbelargerthana1toensurethatR≤ 1and,hence,to en-surestability.Thisisastrongerrequirementforstabilitythanthat

b1 shouldbe largerthan 1/

π

2:thevalue forb1 shouldinfact be betweena1 and4/

π

2.

Fig. 6. Reflection coefficient of a spurious wave mode for a stable set and for an unstable set of coefficients.

Summarizingtheinequalitiesforthecoefficients

a0

π2 <a1<4_πa20 and

a1<b1<_π42.

(34)

The behaviour of Eq. (32) is the same for every interval

kh ∈ < n

π

/2, n

π

] for n=1,3,5,· · · . If the coefficients a1 and

b1 are chosen in such a way that the roots of the functions f2 and f3 are in these intervals, stability is ensured. Within these constraints,one isfree toapproximatethe dispersionrelation for propagatingwavemodesinanyinterval.Fig.6givesanexampleof whatthereflectioncoefficientforaspuriouswavemodelookslike when either a stable set of coefficients – a0=1.040, a1=0.106 and b1=0.289 – or an unstableset of coefficients– a0=1.050,

a1=0.100andb1=0.310– ischosen.Here,thestabilitycriterion thata1hastobelargerthana0/

π

2isviolated.

6. Numericalimplementation

The absorbingopen boundaryconditiondeveloped inthe pre-vious sections is suitable for use in any general-purpose phase-resolving free-surface wave model. Here, we have implemented the boundary condition in ComFLOW_. ComFLOW _{was developed} for offshore applications with extreme waves interacting with structures. The numericalmethodin ComFLOW approximatesthe Navier-Stokes equation on a fixed 3D cartesian grid. The free surface is advected by means of an improved Volume-of-Fluid method.The velocity is solved explicitly,whereas the pressureis solvedimplicitlyfromaPoissonequation.

Fig. 7. Definition of the solution variables at the boundary. The boundary condition is applied to solve for p i . It is positioned in the center of the mirror cell outside the domain.

Whenthesuperscriptsnthatindicatethetimelevelareadded to the solution variables, the discrete system of equations be-comes: Mun+1 d =0, (35) Vu n+1 d − und



t =−C



un d



un d− 1

ρ

Gpnd+1+

ν

Du n d+F n d. (36)

Here,Misadivergencematrixforthediscretecontinuity equa-tion, C



un

d



the convective matrix, D a diffusive matrix and G a

gradient matrix.V is a diagonal matrixthat contains thecontrol volume size,p_d isa vector containingthe discrete pressures and

Fd isavectorthataccountsforthediscreteexternalforce.Forthe

convectivetermthisnotationhasbeenchosentoshowthatitisa nonlineartermandthat elementsofthevector withthediscrete velocitiesudhavebeenusedtoconstructthematrix.

Now,thepredictorvelocityu˜d isintroduced.Thisauxiliary

vec-torwillcontainthecontributions ofconvection,diffusionand ex-ternalforcingattheoldtimelevel:

˜ un_d=un d−



tV−1



C



un d



un d−

ν

Du n d− F n d



. (37)

Withthe predictor velocity, the discrete momentum equation becomes: un+1 d =u˜ n d−



tV−1 1

ρ

Gpnd+1. (38)

ThemomentumequationissubstitutedintocontinuityEq.(35). The pressure vector on the new time level tn+1 _{remains on the} left-handsideoftheequation. Thepredictorvelocityisshifted to the right-hand side ofthe equation. With the property that G= −MT _{weobservethatadiscretePoissonequationforthepressure}

isobtained: MV−1_MT_pn+1 d =

ρ



tMu˜ n d. (39)

Whenthe pressure vector atthenewtime level hasbeen re-solved,thevelocityvectoratthenewtimelevelun+1

d canbefound

fromEq.(38).

ComFLOW_{solvesforvelocitiesandpressuresand,therefore,the} GABC alsoneeds tobe expressed interms ofthesame variables. The velocity at the boundary ub is the derivative of the

poten-tialinthe directionnormalto theboundary. The pressure atthe boundary p_b followsfrom thelinear momentum equation in (6). Fornotationpurposes,pb isdefinedasp/

ρ

.TheGABCintermsof

pressureandvelocityatthedomainboundarythenbecomes:

−



gh



a0+a1h2

∂

2

∂

z2



ub+



1+b1h2

∂

2

∂

z2



pb=gzc+Rin, (40)

wherezc istheverticalpositionofthecenterofthecell.

InComFLOW,thesolutionvariablesarestaggeredwithinacell, seeFig.7forthedefinition ofthevariablesatthe boundary.Cell

labelingis used to indicate completely filled cells inside the do-main(F)andmirrorcellsbeyondtheboundary(O).Theindexi is usedtorefertocellsinhorizontaln-direction.Thedomain bound-aryischosen suchthatitcoincides withthepositionofthe hori-zontalvelocityub.

ItisessentialthatthevelocityandpressureintheGABCare de-finedatthesamepositiontopreventphasedifferencesthatinduce reflection.Thepressureattheboundarypbisobtainedfromlinear

interpolationbetweenthepressures on eitherside ofthe bound-ary:

pb= 1

2

(

pi−1+pi

)

. (41)

It isequallyessentialthat thepressureandthevelocity atthe boundaryaredefinedatthesamepointintimetoprevent reflec-tion.Foraboundaryconditionintermsofpressures, thepressure attheboundaryisdeterminedattimetn+1_{.Thehorizontal} veloc-ityatthenewtimelevelu_bn+1 canbeeliminated bymeansofthe momentumequationattheboundary(38).

The second derivatives in vertical direction in the GABC are approximated by the operator Z. The operator has been derived fora stretched grid in the vertical directionand issecond order accurate.The boundary conditionwith the approximatedvertical derivates is set up an equation for pi,k in a mirror cell outside

thedomain.Itrelates threevelocitieswithtwo timesthree pres-suresoneithersideofthedomainboundary.Inthisway,adiscrete equationfortheabsorbingboundaryconditionisobtained:



Li−1 ,k−1 Li−1 ,kLi−1 ,k+1 pi−1 +



Li,k−1 Li,kLi,k+1 pi =gzc+



Ni−1 ,k−1 Ni−1 ,kNi−1 ,k+1 u˜i−1 +Rin, (42) inwhichthematrixcoefficientsforthepressureareequalto: Li−1,k−1=

−

χτ

+1 2

ψ

Zk−1 Li−1,k=−

ϕτ

+1₂−

χτ

−1 2

ψ

Zk Li−1,k+1=

−

χτ

+1 2

ψ

Zk+1 Li,k−1 =

χτ

+1 2

ψ

Zk−1 Li,k=

ϕτ

+1₂+

χτ

+1 2

ψ

Zk Li,k+1 =

χτ

+1 2

ψ

Zk+1,

andthecoefficientsforthevelocitiesareequalto: Ni−1,k−1=

χ

Z1

Ni−1,k=

ϕ

−

χ

Z2 Ni−1,k+1=

χ

Z3.

Intherelationsabove,thesecoefficientswereused:

ϕ

=



gha0,

χ

=



gha1h2,

ψ

=b1h2and

τ

=



t



x.

Theboundaryconditionfeaturesonlypressuresatthenewtime leveltn+1 _{ontheleft handside,andhorizontalvelocitiesthat} in-cludeconvective and diffusive terms onthe old time level tn _on

therighthandside. Thestructure ofthediscrete GABC issimilar tothatofPoissonEq.(39).

Near thebottomandnearthefreesurface, theverticalsecond derivativescannot be continuedbecause therequiredinformation isnot available. Nearthe bottom,we can use(9) tosee that the verticalderivative ofvelocityandpressure arezeronearthe bot-tomatz=−h.Nearthefreesurfacesuchinformationisnot avail-able,see Fig. 8.In this figure, cell labeling hasbeen used to in-dicatethecelltype.Inemptycells(E)noequationsaresolved,in

Fig. 8. At the free surface near the boundary no second derivative is calculated. Instead, a Sommerfeld equation is solved at the SO-cell boundary.

fluidcells(F)bothgoverningequationsaresolvedandsurfacecells (S)arepartiallyfilledcellsinbetweenE-cellsandF-cellsinwhich boundary conditions are solved. Hence, no pressures and veloci-tiesaredeterminedabovethefreesurface.Andamirrorcondition, suchasthatnearthebottomisnotavailable.Sincethelargest ve-locityvariationoccursnearthefreesurface,itseemsimportantfor aproperlyfunctioningGABCtoaccuratelyapproximatethesecond derivativethere.

The mostobvioussolution isto derive an operator fora one-sidedsecond derivative.The operatoronlyuses solutionvariables frombelow the free surface. Unfortunately, applying a one-sided operator for the second derivative near the free surface resulted inunstablesimulations.Instead,a compromisewasimplemented. Inthecellcontainingthefree surface,nosecond derivativeis cal-culated. In theO-cell nearest to theS-cell inside the domain, an ordinarySommerfeld equation,suchas(13),issolved witha well chosen valueforc0,typicallythephasevelocity thatisassociated withthepeakfrequencyofthespectrum. Thisequationis imple-mented atthecost ofaccuracy,butatthemomentan alternative isnotavailable.

7. Reflectioninirregularwavesimulations

Theperformanceoftheabsorbingboundaryconditionistested bymeansofsimulations.Theperformance isbetterwhenthereis lessreflection. Thereflectioncoefficientisobtainedfromirregular wave simulations; firstly, withthe GABC atthe downstream end of the domain for the situation without incoming waves at that end;andsecondly,withtheGABCattheupstreamendofthe do-main (leftinFig.9c) forthe situationthat incomingwavesenter thedomain whileoutgoingwavesare absorbed.Reflection coeffi-cientsfromthesimulationsarecomparedtothetheoretical reflec-tioncoefficient.

Therearewaystodeterminethepropagationdirectionandthe frequencycontentfromasetofwavesignals,measuredatmultiple locations[25].Thesemethods,however,arebasedonlineartheory and perform worse aswaves become steeper. Furthermore,they willwronglyattributenumericaleffects,such asphaselaggingor waveenergydissipationtopropagatingwaves.

Therefore, another method to assess the performance of the GABC in practical circumstances is proposed. First a wave

simu-Fig. 9. Simulation in a long domain, simulation in a short domain with a GABC on the right where only wave absorption takes place, and simulation in a short domain with a GABC on the side where generation and absorption take place and full reflection on the other side.

lationwithwavespropagating inpositivex-directionisperformed in an infinitely long domain. ’Infinite’, in this sense, means that thedomainlengthischosensuch,thatduringtheentire measure-ment time, reflected waves cannot reach the measurement loca-tion.Thenanothersimulationisperformed.Thissimulationisthe same astheprevious simulation inevery respect, exceptthe do-mainlength.ThedomainisshortwithaGABCattheoutflowend. Registrationsofthesurfaceelevation,takenatthesamepositions, arecomparedtomeasurementsontheinfinitedomain.Everything beingthesame,thedifferencecanonlybeattributedtothe bound-ary procedure. The infinitedomain andthe shorter domain with the GABC areshown inFig.9; thefigure alsoindicates the mea-surementpositionsinthemiddleoftheshortdomain.

OnesimulationinashortdomainisperformedwiththeGABC ontheside whereonlywave absorptiontakesplace(seeFig.9b). AnothersimulationisperformedwiththeGABConthesidewhere generation and absorption take place (see Fig. 9c). In the lat-ter simulation, the opposingboundaryis fullyreflecting to make waves go in two directions. Velocities, pressures and the surface elevationfromacolumnintheinfinitedomainareusedas undis-turbedincomingcharacteristicforuseintheboundarycondition.

ThesimulationwiththeGABC onthesidewhereonly absorp-tiontakesplaceisalsocomparedtoasimulationwithadissipation zone. A dissipation zone was implemented as in Westhuis [27], whichwasfurtherstudiedinComFLOW _{byMeskers[19].Itworks} byaddingaterm,amomentumdampingfunction,tothepressure boundaryconditionatthefreesurfacethatcounteractsthevertical velocity:

p=min

(

0,

β(

x− xstart

)

)

w. (43)

Here, pisthe pressureatthe free surface,

β

isa coefficient that determinesthelinearincreaseofthedampingfunctionfroma po-sition xwherethe startofthe dissipationzone hasbeen defined andw isthe vertical velocity atthe free surface. The slope

β

of

Table 1 Domain and grid.

Parameter Symbol Value

Water depth h 100 m

Domain length  ∞ 10000 m

Domain length  400 m

Measurement position xm 200 m Horizontal grid distance x 1 m Vertical grid distance zmin 0.2 m Vertical grid distance zmax 4.95 m Vertical stretch factor ξ 1.05 [-]

Time step t 0.01 s

Simulated time T 600 s

Table 2

Coefficients of the boundary procedures. Parameter Symbol Value GABC coefficient a0 1.040 [-] GABC coefficient a1 0.106 [-] GABC coefficient b1 0.289 [-] Outgoing phase velocity c/ gh 0.58 [-] Dissipation zone length  zone 600 m Slope damping function β 0.05 [-]

the damping function is typically low as to reduce reflectionby thedampingfunctionitselfasmuchaspossible.Amassdamping functionwasnotapplied.

The setup of the simulations is outlined in Tables 1 and 2. Thedomainis2D, whichmeansthatonly long-crestedwavesare considered.Thedomain sizesandthegriddistances arestatedin

Table1. Thecellshave auniformsize in thehorizontaldirection andarestretchedintheverticaldirection.Inordertobesurethat notonlythemeasurementpositioninspace,butalsothe measure-ment positions in time are the same throughoutall simulations, thetimeincrement



tiskeptfixedandsmall.

Thesimulationsarestartedfromrest:att0therearenowaves inthedomainandthevelocitiesuarezero.Wavesareimposedon theleft of thedomain by means oflinear potential theory.With alinearramp functionthesignals ofsurfaceelevationand veloc-ities atthe inflowboundary are gradually builtup over a period of 20 s. The irregular wave signal consists of a superposition of regular wave components that accord with a realistic JONSWAP spectrumwithsignificantwaveheightHm0=4mandpeakperiod

Tp=15s,seeHolthuijsen[17].

The coefficients of the outflow boundary conditions are in

Table2.ThecoefficientsfortheGABCaretunedinsuchawaythat thereflectioncoefficientovertherangekh_{∈ <}0,6]isneverlarger than 2%. The Sommerfeld boundary conditionatthe free surface hasonly coefficient

η

totune. The bestchoice for

η

is thevalue ofthephasevelocityassociatedwiththepeakperiodofthe spec-trum.

Meskers [19] gives the optimal configuration of a dissipation zone when linear theory is assumed. The length of the dissipa-tion zone, expressed in number of wave lengths, will be deter-minedbasedonthe wavelength associatedwiththepeak period ofthespectrum.Theoptimalnumberofwavelengthsforthis sim-ulation,whenatheoreticalreflectioncoefficientof2%isdesired,is two. The slope of the dampingfunction

β

is basedon the peak frequency.

Thewavesignalatthemeasurementlocationintheinfinite do-main is subtracted from the wave signals in the domains where boundary procedures were applied, so that a reflection signal is obtained. The wave and reflection signals are decomposed into their Fouriercomponents.The Fourier componentsare then con-vertedtobothspectraandreflectioncoefficients.Thisisshownin

Fig. 10. Reflection for the dissipation zone. In (a) the reflection coefficient and in (b) the reflection spectrum.

Figs.10–12.Notethatthe verticalaxesofthesefiguresare differ-ent.

The results with the dissipation zone in Fig. 10 were unex-pectedand largerthan the theoretical value. This isdue to non-lineareffects inthedissipationzone thathavetransferred energy fromhighertolowerwavefrequencies.Dissipationzonestendnot toworkwell forlongerwaveswithlowerfrequenciesbecausethe verticalvelocitiesaresmall,whichexplainsthelargereflection co-efficientat lower frequencies. The dissipation zone doesperform wellfortheshortwaveswithkh-valueslargerthan3.

The reflection coefficient and the reflection spectrum for the GABC are shownin Fig. 11. In the range fromkh₌0 to kh₌5_, thereflectioncoefficientofaround5%ishigherthanthe2%itwas designedfor.Investigation ofthe reflectioncoefficientsby means of two sets of regular wave simulations (shown in the inset of

Fig.11a),onewithamaximumsteepnessH/Lof1· 10−2_,theother withamaximumsteepnessof5· 10−2_{,leadstotheconclusionthat} thehigherreflectioncoefficient is mostlydueto theSommerfeld boundaryconditionimposed intheS-cellatthefreesurfacenear the boundary and not due to inaccuracies in the Fourier trans-form,norduetononlineareffects.TheSommerfeldboundary con-ditionaffectstheshortestcomponentsoftheirregularwavetrain, abovekh=5,most.AnotherreasonwhytheGABCdoesnotmatch

Fig. 11. Reflection for the GABC. In (a) the reflection coefficient of the irregular wave simulation. The inset magnifies the range between kh = 0 and 5 and also com- pares to two sets of reflection coefficients from regular wave simulations, one with a maximum steepness of 1 · 10 −2 , the other with a maximum steepness of 5 · 10 −2 . In (b) the reflection spectrum.

thetheoreticalreflectioncoefficientforthelargerkh-valuesis be-causetheshortercomponentsintheirregularwavesare resolved byfewercellsperwavelength.

ThereisadifferencebetweentheGABCinFig.11andthe dis-sipation zone in Fig. 10 in the value of the reflection coefficient atkh₌0.ImposingwavesbymeansofDirichletvelocityboundary conditionsisknowntoadd watertothedomain.Thisisapparent fromthenon-zerovalueofthereflectionspectrumofthe dissipa-tionzoneatkh=0.TheGABCmakestheadditionalwaterflowout ofthedomainbecauseitisformulatedinawaythat,onaverage, makesthefreesurfacetendtothemeanfreesurface.Forthat rea-son,asystematicincreaseofthewaterlevelisnotobservedwith theGABC.

Fig.12showsthereflectioncoefficientandthereflection spec-trumwhen theGABC on theupstream side of ashort domainis both generatingand absorbing.The opposingboundary wasfully reflective so that the wave signal at the measurement location inthe middleofthe domaincontainswave componentsgoing in downstreamandinupstreamdirectionbydesign.Thewave signal froma simulationina longdomain withthe incomingboundary ’infinitely’farawayfromthefullyreflectiveboundaryissubtracted fromthe wave signal in the shortdomain. Because ofthis setup thewave signal containsmainly the longerwave components up tokh=1.2,becausetheypropagatefasterandarriveearlieratthe

Fig. 12. Reflection for the GABC in the presence of incoming wave. In (a) the re- flection coefficient and in (b) the reflection spectrum.

reflective boundary. The difference in wave signal from the long domainandthe shortdomainisattributedtospurious reflection. From Fig.12we findthatthe reflectionspectrumandthe associ-atedreflectioncoefficientoftheimplementedGABCaresomewhat higherthanthetheoreticalreflection.

The GABC,both on therightwhen onlyabsorbing andonthe leftwhengeneratingandabsorbing,behavesasexpected,butwith more reflection. The amount of wave energy that is reflected is smallcomparedtotheinputspectrumandatacomparablelevelas spillingbeachesinwavetanksusedforexperiments.Therearetwo main reasonsforthedifferencebetweentheory and implementa-tion.ThefirstreasonistheSommerfeldconditionintheSO-cellat theboundary, seeFig.8.Thisrelationwastunedto thepeak fre-quencyofthespectrumandwillespeciallycausereflectionofthe longer and shortercomponents away fromthe tuning frequency. Thesecond reasonisthatnonlineareffectsneartheboundaryare notaccountedforinthedesignoftheGABC.

8. Practicalexampleofavalidatedsimulation

WemeantoshowarealapplicationoftheGABCforasituation withunidirectional incomingwaveswhereasimulationwith dis-sipationzonesontheupstreamandthedownstreamsidethatare

Fig. 13. Photo of the experiment and snapshot of the simulation.

two typical wave lengths long,would be prohibitively expensive. Thesimulationisvalidatedbymeansofanexperiment.

During the development of ComFLOW a specifically designed experiment was performed at the Maritime Research Institute Netherlands (MARIN) to validate the simulations. It was for a schematized model of a restrained semi-submersible with two columnsandapontoonunderwater,seeFig.13a.Theexperiments wereperformedatscale1:50.

The experiment, withidentification number202006, was per-formedinalongwavetank.Wavesweregeneratedwithapivoting waveboardatoneendofthetank.Attheopposingendaspilling beachwaspresenttoinducewave breakingandreducereflection. Themodelwasplacedalargedistanceawayfromthewaveboard. Asimulation wasperformedtocompare tothe experiment. A bettercomparisonwiththeexperimentis obtainedwhen nonlin-ear wave kinematics are used forinput atthe boundary instead oflinearizedwave kinematicsfromsomething likea wave gauge. Toobtainthenonlinearkinematics,measurements oftherotation ofthe wave boardwere takenandfed into a nonlinearpotential flowmodel[27]withadomainthatwassolongthatwave reflec-tion fromthe downstream end did not occur. The nonlinear po-tentialflowmodelpropagateswavesatamuchreduced computa-tionalcost compared to ComFLOW_{.Kinematics at theposition of} theinflowboundaryofComFLOW were takenfromthe nonlinear potentialflowmodel.

The inflow boundary in ComFLOW _{is located at x}=_{5156 m} withrespect to thewave board.This iswhere the GABC is used together with the kinematics from the nonlinear potential flow modelastheincomingwavecharacteristic.Theoutflowboundary withaGABCispositionedatx=5456m.300cellswereusedinx -directionwithoutstretching.Intheverticaldirection75cellswere appliedwithcellsizestretchinginsuchawaythateveryn-thcell

Fig. 14. Numerical results compared with the experimental results of test 202006.

awayfromthecellwiththeminimumcellsizeisbiggeraccording to:

hn=

(

1+

γ

)

nhmin (44)

with

γ

astretchingcoefficientandnacounterofhowmanycells arebetweenitandthecellwithgridsizehmin atthefocuspoint. Inz-direction,

γ

=5%andthefocuspointofthestretchingis po-sitionedatthemeanfreesurfacelevelz=0.Nearthefreesurface



zmin is0.5m,nearthebottom



zisnearly 14m.Iny-direction the boundariesof the domain were chosen to coincide withthe tank walls, each positioned 100 m away from the center of the structureaty₌0;75cellsareusedwithastretchingcoefficientof 2%,



ymin being 1.5mandthe largest cellsnearthey-boundary slightlyover2m.Thesimulationtimeis120s

The center position of the structure was at x₌5306 m. The GABCwasusedattheupstreamboundaryandatthedownstream boundarywithout incomingwaves. Fig.13bshows asnapshot of the3Dsimulation.

Wave probe 2waspositioned atx=5238mon theupstream side of the structure in the middle of the tank. Fig. 14a shows thesurfaceelevationasa functionoftime forwave probe 2and comparesittoComFLOWsimulationswiththeGABC,adissipation zoneandtheSommerfeld boundaryconditionwithc0=0.56



gh. ThereisgoodagreementbetweentheGABCsimulation,the simu-lationwithadissipationzoneandtheexperimentwiththesurface

elevationintheGABCsimulationmarginallyhigherthaninthe ex-perimentforthehighestwaves.Thesimulationwiththe Sommer-feldboundaryconditiongivesdifferentresultsfromt=650s on-wardsasaresultofspuriousreflection.

Asimilar resultisfound fortheimpacts in theregistrationof pressure sensor P12 on the front column of the structure (posi-tionedatx₌5266mandatz₌₋₁₀mwithrespecttothemean watersurface), whichis shownin Fig.14b. There isoverall good agreement betweenthe pressure in ComFLOW _{andin the} exper-iment for the simulation with the GABC and with the dissipa-tion zone. The simulation with theSommerfeld boundary condi-tiondemonstratesthatspuriousreflectioncaninfluencetheimpact pressures.

9. Conclusion

Thisarticle discusses thebackground, derivation, implementa-tionandresultsofanopenboundaryconditionfordispersive free-surfacewaves. We call it a Generating Absorbing Boundary Con-ditionorGABC. Thestartingpoint forthe derivationofthe GABC istheSommerfeldcondition, whichisperfectlyabsorbingforone wave component with one propagation velocity. This velocity is specified in theSommerfeld conditionby means of a tuning pa-rameter.Itwasfoundthattherangeofabsorbedwavecomponents can be extended by replacing the tuning parameter withan ap-proximationofthelineardispersionrelationintermsofthewave number.When second derivativesofthesolutionvariables inthe verticaldirectionare substitutedforthewave number,ityields a GABC witha low2% reflectionforarangeofcomponentsinstead ofasinglecomponent.StabilityissuesloomwhentheGABCis im-plemented,butthemechanism,bywhichinstabilitiesoccur,iswell understoodandstabilitycriteriahavebeenformulated.The stabil-itycriteriahardlyrestraintheabsorbingperformanceoftheGABC. In irregular wave simulations, the GABC has similar perfor-mance to a dissipation zone that was two typical wave lengths longintherangeofwavenumbersofinterest.Itwasfoundin nu-merical simulations with the GABC that the reflection coefficient can be as low as 5% for mildly steep waves. This is somewhat higher than the 2% that wasderived from theory, but compara-bleto theamountofreflectionthat issaid to beobtainedin ex-perimentalbasinsandflumes.The differencebetweentheoryand practice hereis mainly due to theimplementation near the free surface.

ItwasdemonstratedthattheGABCcanbeusedinanactual3D simulationof steepwave interaction withfloating structures. The simulationwasvalidatedwithexperimentalresults.Thesimulation withtheGABChassimilarperformancetothesimulationwiththe dissipationzone,atamuchreducedcomputationaleffort.

DeclarationofCompetingInterests

Theauthorsdeclarethattheyhavenoknowncompeting finan-cialinterestsorpersonalrelationshipsthatcouldhaveappearedto influencetheworkreportedinthispaper.

The authors declare the following financial interests/personal relationshipswhichmaybeconsideredaspotential competing in-terests:

CRediTauthorshipcontributionstatement

PeterWellens:Conceptualization,Methodology,Software, Vali-dation,Formalanalysis,Writing-originaldraft,Visualization.Mart Borsboom:Conceptualization,Formalanalysis,Writing -review & editing.

Acknowledgement

This research has been supported by the Dutch Technology FoundationSTW, applied science division of NWOand the technol-ogyprogrammeoftheMinistryofEconomicAffairsinThe Nether-lands (GWI.6433). MARIN andespecially Tim Bunnik is acknowl-edgedfororganisingtheComFLOWJIPs.Andwewouldliketo ex-tendourgratitudetoArhurVeldmanfromRijksuniversiteit Gronin-genforhisfeedbackonthedevelopmentinthisarticle.

References

[1] Berenger JP . A perfectly matched layer for the absorption of electromagnet- ic-waves. J Comput Phys 1994;114(2):185–200 .

[2] Blayo E , Debreu L . Revisiting open boundary conditions from the point of view of characteristic variables. Ocean Modell 2005;9(3):231–52 .

[3] Borsboom M , Groeneweg J , Doorn N , van Gent M . Flexible boundary conditions for a boussinesq-type wave model. In: Edge BL, Hemsley JM, editors. Ocean wave measurement and analysis : proc. of the 4th int. symposium, WAVES 2001, San Francisco, CA. American Society of Civil Engineers; 2001. p. 884–93 . [4] Carpenter KM . Radiation conditions for the lateral boundaries of limited-area

numerical-models. Q J R Meteorolog Soc 1982;108(457):717–19 .

[5] Collino F , Joly P . New absorbing boundary-conditions for the finite-element solution of 3d maxwells equations. IEEE Trans Magn 1995;31(3):1696–701 . [6] Dgaygui K , Joly P . Absorbing boundary-conditions for linear gravity-waves.

SIAM J Appl Math 1994;54(1):93–131 .

[7] Düz B , Borsboom MJA , Veldman AEP , Wellens PR , Huijsmans RHM . An ab- sorbing boundary condition for free surface water waves. Comput Fluids 2017;156:562–78 .

[8] Engquist B , Majda A . Absorbing boundary-conditions for numerical-simulation of waves. Math Comput 1977;31(139):629–51 .

[9] Givoli D . Non-reflecting boundary conditions. J Comput Phys 1991;94(1):1–29 . [10] Givoli D . High-order local non-reflecting boundary conditions: a review. Wave

Motion 2004;39(4):319–26 .

[11] Givoli D , Neta B . High-order non-reflecting boundary conditions for dispersive waves. Wave Motion 2003;37(3):257–71 .

[12] Grote MJ , Keller JB . Nonreflecting boundary conditions for time-dependent scattering. J Comput Phys 1996;127(1):52–65 .

[13] Hagstrom T , Mar-Or A , Givoli D . High-order local absorbing conditions for the wave equation: extensions and improvements. J Comput Phys 2008;227(6):3322–57 .

[14] Hagstrom T , Warburton T . A new auxiliary variable formulation of high-order local radiation boundary conditions: corner compatibility conditions and ex- tensions to first-order systems. Wave Motion 2004;39(4):327–38 .

[15] Higdon RL . Absorbing boundary-conditions for difference approximations to the multidimensional wave-equation. Math Comput 1986;47(176):437–59 . [16] Hirt CW , Nichols BD . Volume of fluid (VOF) method for the dynamics of free

boundaries. J Comput Phys 1981;39(1):201–25 .

[17] Holthuijsen LH . Waves in oceanic and coastal waters. Cambridge University Press; 2007 .

[18] Kleefsman K . Water impact loading on offshore structures - a numerical study. University of Groningen, Netherlands; 2005. Ph.D. thesis .

[19] Meskers G . Realistic Inflow Conditions for Numerical Simulation of Green Wa- ter Loading. Delft University of Technology; 2002. Master’s thesis .

[20] Orlanski I . Simple boundary-condition for unbounded hyperbolic flows. J Com- put Phys 1976;21(3):251–69 .

[21] Perkins AL , Smedstad LF , Blake DW , Heburn GW , Wallcraft AJ . A new nested boundary condition for a primitive equation ocean model. J Geophys Res-O- ceans 1997;102(C2):3483–500 .

[22] Rivera-Arreba I , Bruinsma N , Bachynski E , Viré A , Paulsen B , Jacobsen N . Mod- eling of a semisubmersible floating wind platform in severe waves. In: Pro- ceedings of the 37th international conference on offshore mechanics and arctic engineering - 2018 - Vol 9; 2018 .

[23] Van Dongeren AR , Svendsen IA . Absorbing-generating boundary condi- tion for shallow water models. J Waterway Port Coastal Ocean Eng-ASCE 1997;123(6):303–13 .

[24] Verboom GK , Slob A . Weakly-reflective boundary-conditions for 2-dimensional shallow-water flow problems. Adv Water Resour 1984;7(4):192–7 .

[25] Voogt A , Bunnik T , Huijsmans R . Validation of an analytical method to calcu- late wave setdown on current. In: Proceedings of the 24th international con- ference on offshore mechanics and arctic engineering - 2005 - Vol 1, Pts A and B; 2005. p. 875–9 .

[26] Wellens P . Wave simulation in truncated domains for offshore applications. Delft University of Technology; 2012. Ph.D. thesis .

[27] Westhuis JH . The Numerical Simulation of Nonlinear Waves in a Hydrodynamic Model Test Basin. University of Twente; 2001. Ph.D. thesis .

[28] Youngs DL . An Interface Tracking Method for a 3D Eulerian Hydrodynamics Code. Technical Report. Atomic Weapons Research Establishment; 1987 .