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

The boundaries of numerical domains for free-surface wave simulations with marine structures generate spurious wave reﬂection if no special measures are taken to prevent it. The common way to prevent re-ﬂection is to use dissipation zones at the cost of increased computational effort. On many occasions, the size of the dissipation area is considerably larger than the area of interest where wave interaction with the structure takes place. Our objective is to derive a local absorbing boundary condition that has equal performance to a dissipation zone with lower computational cost. The boundary condition is designed for irregular free-surface wave simulations in numerical methods that resolve the vertical dimension with multiple cells. It is for a range of phase velocities, meaning that the reﬂection coeﬃcient per wave component is lower than a chosen value, say 2%, over a range of values for the dimensionless wave number kh . This is accomplished by extending the Sommerfeld boundary condition with an approximation of the linear dispersion relation in terms of kh , in combination with vertical derivatives of the solution variables. For this article, the boundary condition is extended with a non-zero right-hand side in order to prevent wave reﬂection, while, at the same time, at the same boundary, generating waves that propagate into the domain. Results of irregular wave simulations are shown to correspond to the analytical reﬂection coef-ﬁcient for a range of wave numbers, and to have similar performance to a dissipation zone


Introduction
We employ a numerical method for simulations of extreme, irregular free-surface wave interaction with (floating) structures at sea.A major issue is that waves reflect from computational domain boundaries when no special measures are taken to prevent reflection.Dissipation zones are often used near domain boundaries to resolve the issue.Dissipation zones induce a decrease of wave energy so that when waves reach the boundary and reflection occurs, the wave height has reduced to such an extent that spurious reflection does not interfere with the processes near the structure of interest.
Dissipation zones need to be several typical wave lengths L p long to be effective, where L p is the wave length associated with the peak frequency of the irregular wave spectrum.Consider, for boundary condition that does not require any additional areas attached to the computational domain.Without any additional areas, local absorbing boundary conditions are the more computationally efficient alternative to dissipation zones for preventing wave reflection.Here, a local absorbing boundary condition for long-crested irregular waves is derived with similar spurious wave reflection as a two-wave-lengths-long dissipation zone at a reduced computational effort.
The absorbing boundary condition is based on a local boundary condition operator in combination with a rational approximation of the free surface wave dispersion relation and with second derivates of velocity and pressure in vertical direction, derived in Wellens [26] .This approach extends the range of wave numbers over which spurious wave reflection is low compared to the Sommerfeld boundary condition.The range was further extended in terms of both dispersion (wave number) and direction of the wave components in Düz et al. [7] , by multiplying the boundary condition with an additional boundary condition operator.
What the current article adds to the earlier work is that the absorbing boundary condition is extended with a non-zero righthand side to generate incoming waves, while outgoing waves are absorbed at the same time.Some approximations of the dispersion relation lead to stability issues.A stability analysis is presented to prevent instabilities.The absorbing boundary condition can be combined with any numerical method for free-surface wave simulations that resolves the vertical dimension with multiple grid cells.For the purpose of wave impacts on marine structures, the boundary condition was implemented in ComFLOW [18] , a numerical method for free-surface waves that is based on the Navier-Stokes equations with a Volume of Fluid approach for the free suface.Both the absorbing performance and the combined generating and absorbing performance are verified by comparing the spurious reflection coefficients of irregular wave simulations to the theoretical reflection coefficients.And finally, it is demonstrated that the generating absorbing boundary condition (GABC for short) not only works in theoretical wave propagation simulations, but can be used as a computationally efficient alternative to dissipation zones in a simulation of wave interaction with -and wave impact loads on -a structure in the free surface, which is validated by an experiment.

Governing equations
The numerical method we employ for wave interaction with floating structures solves the Navier-Stokes equations for a fluid in the presence of a free (liquid) surface in a cartesian grid domain, where the grid lines are colinear with the components of the cartesian right-handed axis system with x = (x, y, z) T .The Navier-Stokes equations impose conservation of mass and conservation of momentum, respectively: in which the fluid velocity vector is represented by u = (u, v , w ) T .Scalar p is the pressure, ρ is the fluid density and ν is the kinematic viscosity.f represents gravity (0 , 0 , −g) T acting as an external body force.The indicator function for free surface, S , is displaced in time and space by: in a manner derived from Hirt and Nichols [16] , in combination with a local height function in a 3x3x3 block of cells [18] .The free surface is reconstructed using PLIC [28] .The simulations in Sections 7 and 8 are single phase (water and void) with a boundary condition for the pressure.The following boundary condition for the pressure is imposed by means of interpolation towards the reconstructed position of the free surface: where n s is the direction normal to the free surface and p 0 is the atmospheric pressure.The surface tension is represented by σ and the curvature of the free surface by κ, which is computed from the Near the boundary, sufficiently far away from a structure in the flow, we assume that the flow is irrotational.Near the boundary, we also assume that we can apply a local linearization.This allows for the introduction of a linear potential function , whose spatial derivatives yield the velocity in the direction of the derivative ∇ = u .With the potential, the continuity equation and linearized momentum equation become: The boundary conditions at the free surface simplify to: Both are imposed at the mean free-surface position at z = 0 .ζ is the free-surface position with respect to z = 0 .The propagating periodic solution to the system of equations in (5) to (8) in a finite water depth h , where the bottom is considered impermeable, i.e. ∂ /∂ z = 0 at z = −h, is given by: Here, ζ a,j is the wave amplitude, ω j is the wave frequency, k j is the wave number and j is an indicator for the wave component under consideration.
Eqs. ( 6)-( 8) can be combined with solution (9) to obtain the dispersion relation: which for our purposes is be rewritten into a formulation for the phase velocity c j : By dispersion we mean that each wave frequency has a different phase velocity.This aspect is particularly difficult for boundary conditions for irregular free-surface waves that are considered to be composed of many wave components with each their own amplitude, frequency and, hence, phase velocity.
The system of equations in ( 5) to (8) can also be recombined into a planar wave equation for the potential at the mean free surface: Wave Eq. ( 12) is the starting point for the derivation of wave boundary conditions at domain boundaries.

Background
For a local generating absorbing boundary condition three aspects are relevant: (1) wave directionality, (2) wave dispersion and (3) incoming and outgoing waves, or -in other words -generation and absorption of waves.Literature predominantly discusses absorbing boundary conditions, which are also called non-reflecting boundary conditions or radiation conditions, with an emphasis on wave directionality, see Givoli [9] .It is relevant for the derivation of our boundary condition to discuss all three aspects, starting with wave directionality.
An example of an absorbing boundary condition for wave Eq. ( 12) is the Sommerfeld equation:

∂ ∂t
in which c 0 is a typical value that approximates the apparent phase velocity at the boundary.The apparent phase velocity for waves with phase velocity c that approach the boundary under angle θ with respect to the boundary's normal direction is higher and equal to c /cos ( θ ).If c 0 is the same as the apparent phase velocity, then there is no reflection.The more c 0 is different from the apparent phase velocity of a wave component, the more it is reflected.The amount of reflection is found from the reflection coefficient: which is the ratio of amplitudes of the outgoing and the reflected wave components.One could say that the Sommerfeld boundary condition is only accurate, i.e. has low reflection, for a small range of wave directions where c 0 is close to c /cos ( θ ).The development of boundary conditions is concerned with extending the range of directions over which low reflection occurs.
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 approximated with derivatives of the solution variables in tangential direction τ along the boundary, see Fig. 1 .Engquist and Majda [8] suggested an order expansion of c /cos ( θ ), where including more terms in the expansion leads to higher derivatives along the boundary and a larger range over which the absorbing boundary condition is accurate.The expansion to O ( 1 ) leads to the Sommer- feld equation with c 0 = c.The expansion to O ( 2 ) gives: From Eq. ( 15) we find that higher order approximations of the dispersion relation give higher derivatives along the boundary in the absorbing boundary condition operator.By means of higher derivatives in the ABC, waves with larger incoming angles to the boundary are accurately absorbed.
The other family of boundary conditions approximates c /cos ( θ ) with derivatives of the solution variables in normal direction to the boundary.Introduced by Higdon [15] , it is given by: Wave reflection becomes smaller as the order of the absorbing boundary condition increases.The order of the boundary condition, in this sense, reflects the number of products J used to construct the operator.Higher order absorbing boundary conditions feature higher derivatives.As the order of the absorbing boundary condition increases, it becomes increasingly difficult to implement the numerical equivalent of the higher derivatives at the boundary.Givoli [10] reports that Higdon operators beyond order three are rarely found in the literature.Collino and Joly [5] introduced auxiliary variables to remedy the difficulties with implementing higher derivatives.Instead of N th derivatives, then a system of N + 1 additional equations is solved at the boundary.
Auxiliary variables have been used in N th order absorbing boundary conditions derived by Grote and Keller [12] , Givoli and Neta [11] and Hagstrom and Warburton [14] , among others.The Hagstrom-Warburton formulation of the auxiliary system of recursive equations is: where a j are coefficients and ξ j are auxiliary variables.The set of equations features none higher than first derivatives in normal direction to the boundary.The major difficulty of this way of formulating boundary conditions are the corner conditions required near the ends of domain boundaries.An overview is provided by Hagstrom et al. [13] .We did not choose this approach because it proved difficult to combine auxiliary variable absorbing boundary conditions with incoming waves over the same boundary.Dispersion of free-surface waves, the fact that each wave frequency has a different phase velocity, has not been an explicit subject in most of the literature regarding absorbing boundary conditions.As was discussed for wave directionality, the Sommerfeld equation is accurate in a limited range around phase velocity c 0 .The objective of a better absorbing boundary condition is to extend the range of wave frequencies over which the reflection is low.We found two ways to account for dispersion at the boundary: one is to use a low order boundary condition and combine it with an estimator for the apparent phase velocity at the boundary, the second is to use approximation of a non-local absorbing boundary condition with vertical derivatives of the solution variables.As an example of the former, we mention Orlanski [20] .He suggests the use of the Sommerfeld condition, with a dynamic approximation of the phase velocity, obtained from the solution itself: Fig. 2 shows a representation of the phase velocity, when the Orlanski boundary condition is applied to an irregular wave signal conforming to a JONSWAP spectrum with peak period T p = 15 s.It shows rather large discontinuities in the resolved value for c , coinciding with the times where the denominator in (20) becomes zero, which can and will lead to an unstable simulation.There are stabilization methods, such as weighted moving average filtering of c 0 in [3] , but these methods show large reflection for wave components with an apparent phase velocity that is different from the filtered c 0 .
Dgaygui and Joly [6] discuss absorbing boundary conditions for the simulation of free-surface waves in the nz -plane.Herein, an exact, non-local absorbing boundary condition operator is derived.The exact operator is simplified by means of rational approximations to yield an absorbing boundary condition that is local in space and local in time with vertical derivatives of the solution variables along the boundary.Numerical results are presented for simulations of long waves with a zeroth order and a first order boundary condition.With the latter boundary condition, the reflection was generally small; exact figures were not mentioned.
In our studies in Wellens [26] , we found that using rational approximations of the dispersion relation (10) in combination with vertical derivatives of the solution variables yields a considerable increase of the range of wave frequencies and phase velocities that can be absorbed accurately.In Düz et al. [7] , further extension of the range with low reflection was obtained by combining the vertical derivatives of the solution variables with derivatives in normal direction to the boundary as in Higdon [15] .It proved difficult to turn the latter approach into a boundary condition for an open boundary, where waves from inside the domain are absorbed while waves from outside the domain enter at the same time.
Another way to derive boundary conditions is by means of the method of characteristics discussed by Blayo and Debreu [2] , Verboom and Slob [24] and Van Dongeren and Svendsen [23] .Here, the method of characteristics is applied to the wave equation.Consider again wave Eq. ( 12) , but now in one dimension in the direction normal to the boundary: The term on the left-hand side of Eq. ( 21) can be factorized: The term 22) is constant along lines d n/d t = −c.We call this term the incoming characteristic.
A boundary condition for wave Eq. ( 21) is obtained by assigning a value to the incoming characteristic.When no waves enter the computational domain through the boundary, i.e. by assigning a value of 0 to the incoming characteristic, the Sommerfeld equation is obtained.When the incoming characteristic receives a non-zero time-dependent right-hand side value, we are sending in waves over the boundary, while outgoing waves may leave the computational domain unaffected.In this way, an open, or generating absorbing boundary is obtained.
According to Carpenter [4] and Perkins et al. [21] , the Sommerfeld equation with non-zero right-hand side value should be formulated as follows: ∂ ∂t Below, a generating absorbing boundary condition for open boundaries is derived and evaluated, which combines the method of characteristics with a rational approximation of the dispersion relation and vertical derivatives of the solution variables.The rational approximation and the vertical derivatives of the solution variables give an increased range of wave directions and phase velocities for which reflection is low.The same operator used on the left-hand side of the equation is applied to the right-hand side of the equation that prescribes the incoming characteristic, allowing waves to enter the domain while outgoing waves are absorbed.

Derivation
To arrive at an absorbing boundary condition for free-surface wave simulations, we start from Sommerfeld condition (13) and look for ways to include the effect of wave dispersion in the applied c 0 .As explained in the previous section, application of a solution-dependent value (20) is not general enough.On the other hand, we want to avoid the use of multiple operators as in (16) in order not to complicate the numerical implementation with onesided higher derivatives in normal direction to the boundary.Inspired by Dgaygui and Joly [6] , we consider the Fourier transform of the Sommerfeld equation and replace c 0 by the expression for the linear dispersive phase velocity (11) , approximated by a rational [2,2] Padé approximation.
We start with the Fourier transform of (13) : In our interpretation f 1 : kh → c 0 is a constant function that approximates dispersion relation (10) over a range of phase velocities, or equivalently, a range of values for the dimensionless wave number kh. c 0 intersects the dispersion relation at one value for kh .Here, the difference between the functions is zero and, hence, the reflection coefficient is zero, see Eq. ( 14) .
The main thought behind our boundary condition is that c 0 is replaced by a better approximation of the dispersion relation.A rational polynomial in kh that is a better approximation to the dispersion relation than f 1 is given by The coefficients a 0 , a 1 and b 1 can be chosen such that different kh -ranges of the dispersion relation are approximated well.With a 0 = 1 .0 0 0 , a 1 = 0 .150 and b 1 = 0 .367 , the dispersion relation is approximated near kh = 0 with sixth-order accuracy, as can be shown by means of a series expansion.This is illustrated in Fig. 3 .It shows the dispersion relation together with the approximating functions f 1 and f 2 (with the coefficients stated above), and it also shows the reflection coefficient for both approximations of the dispersion relation.The reflection coefficient for f 2 is lower than for f 1 over a large kh -range.Using f 2 in a boundary condition, therefore, makes it a better boundary condition for irregular waves with wave components in the kh -range where the approximation is sufficiently accurate.It is a 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 is a trade-off between accuracy and range.When the range is increased, also the reflection coefficient increases.Fig. 3 also shows function f 2 with the coefficients a 0 = 1 .040 , a 1 = 0 .106 and b 1 = 0 .289 .With these coefficients, the reflection coefficient is non-zero near kh = 0 , but it is less than 2% in the range between kh = 0 and kh = 6 , which is acceptable in many applications.Finding the coefficients for a specific simulation was automated by means of a brute force algorithm that either gives the lowest reflection coefficient for a given range, or the lowest possible area of the reflection spectrum for a given wave spectrum.
When c 0 = f 2 is substituted into boundary condition (24) , we obtain: Note that Eq. ( 26) is multiplied by 1 Boundary condition ( 26) cannot be transformed back to the time domain, because it is nonlinear in k j .k 2 j will need to be eliminated from the boundary condition.
The wave number k j can be found by taking derivatives of the solution in space.Consider the solution of the linear system at the boundary in Eq. (9) .We find that second derivative of the solution in vertical direction is equal to the square of the wave number times the solution itself Eq. ( 27) is substituted in boundary condition (26) to yield and can be transformed back from Fourier space.In the time domain, the absorbing boundary condition becomes It gives little reflection for wave components in the range of kh -values where rational approximation (25) of the dispersion relation is accurate.To demonstrate how the range of low reflection of our boundary condition has increased, the reflection coefficient is compared to the reflection coefficient for the Sommerfeld boundary condition.The coefficients in rational approximation (25) were given the values a 0 = 1 .040 , a 1 = 0 .106 , b1 = 0 .289 .For the coefficient in the Sommerfeld equations a value of c 0 = 0 .6 gh is chosen.
Fig. 4 shows the reflection coefficient for a range of kh -values and angles of incidence θ for our boundary condition and the Sommerfeld condition.The area where the reflection coefficient is below five percent has been enclosed by a bold line; five percent reflection is considered tolerable in practical wave simulations.From Fig. 4 , we find that the area with less than 5% reflection has increased considerably.
The left-hand side of (29) can be viewed as the incoming characteristic of the wave equation, where the effect of wave dispersion has been approximated using (25) .Following Carpenter [4] and Perkins et al. [21] , the incoming characteristic is combined with a non-zero right-hand side consisting of the same operator as on the left-hand side applied to a known incoming wave potential.With a non-zero right-hand side, absorbing boundary condition (29) becomes: in which: We call this boundary condition a GABC, a Generating Reflecting Boundary Condition.

Stability
The range of kh -values for which boundary condition (30) absorbs outgoing waves with little reflection is the range of kh -values for which the applied underlying dispersion approximation is accurate.That range can be tuned to a particular application by a suitable choice of the parameters.We optimized these coefficients for a reflection of less than 2% from kh = 0 to the largest possible value for kh .However, we found that this set of parameters ( a 0 = 1 .050 , a 1 = 0 .100 , b 1 = 0 .310 ) cannot be applied because simulations became unstable.It appears that a 0 , a 1 , b 1 have to satisfy stability constraints.
So far, we have only considered propagating wave modes, which are oscillatory in time, oscillatory in horizontal space, and exponential in the vertical, see (9) .Other wave modes are oscillatory in time are the evanescent waves, which are exponential in horizontal space and oscillatory in the vertical.It can be shown that they have no effect on the absorption properties of (30) , which is why we have not considered them.However, it appears that Eqs. ( 5)-( 8) support yet another type of solution, which we refer to as spurious modes.These are modes that, like the evanescent modes, are exponential in horizontal space and oscillatory in the vertical.But rather than being oscillatory in time, they are exponential in time: Here, subscript s is used to indicate that it concerns a spurious wave mode that grows exponentially in time.When the exponent is positive, the mode is unstable.With boundary conditions that are normally applied, spurious modes cannot exist, i.e. stability is always satisfied.With our boundary condition, however, it appears that certain combinations of a 0 , a 1 and b 1 leave spurious modes undefined, as a result of which a simulation will become unstable.
When an incoming and outgoing solution mode are substituted into boundary condition (29) , it is possible to formulate a reflection coefficient as the ratio of their amplitudes.For stable solutions, the reflection coefficient necessarily has to be smaller or equal than 1: In (32) , a mechanism is recognized by which the reflection coefficient can become larger than 1.Eq. (32) has singularities where the denominator has roots.The following functions are the expressions from the denominator in (32) .33) are chosen in the range where f 1 has imaginary values, then instability cannot occur.This puts restrictions on the values for the coefficients a 1 and b 1 .In Fig. 5 , the function f 2 is plotted with a 0 = 1 and two different values for a 1 , a 1 = a 0 /π 2 and a 1 = a 0 / 4 π 2 .With these values for a 1 , the roots of f 2 are precisely on the outer limits of the range kh ∈ [ π /2, π ].Therefore, the first requirement for stability is that a 1 has a value between a 0 / π 2 and 4 a 0 / π 2 .
The same line of reasoning applies to the function f 3 in Eq. (33) and the value for b 1 .The second requirement for stability, then, is that the value for b 1 is chosen between 1/ π 2 and 4/ π 2 .
The third and final requirement for stability comes from the limit behaviour of the reflection coefficient in Eq. ( 32) .If we consider the limit kh → ∞ of the reflection coefficient, it is found that b 1 should be larger than a 1 to ensure that R ≤ 1 and, hence, to ensure stability.This is a stronger requirement for stability than that b 1 should be larger than 1/ π 2 : the value for b 1 should in fact be between a 1 and 4/ π 2 .Fig. 6.Reflection coefficient of a spurious wave mode for a stable set and for an unstable set of coefficients.

Summarizing the inequalities for the coefficients
The behaviour of Eq. ( 32) is the same for every interval b 1 are chosen in such a way that the roots of the functions f 2 and f 3 are in these intervals, stability is ensured.Within these constraints, one is free to approximate the dispersion relation for propagating wave modes in any interval.Fig. 6 gives an example of what the reflection coefficient for a spurious wave mode looks like when either a stable set of coefficients -a 0 = 1 .040 , a 1 = 0 .106 and 1 = 0 .289 -or an unstable set of coefficients -a 0 = 1 .050 , a 1 = 0 .100 and b 1 = 0 .310 -is chosen.Here, the stability criterion that a 1 has to be larger than a 0 / π 2 is violated.

Numerical implementation
The absorbing open boundary condition developed in the previous sections is suitable for use in any general-purpose phaseresolving 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 numerical method in ComFLOW approximates the 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 pressure is solved implicitly from a Poisson equation.When the superscripts n that indicate the time level are added to the solution variables, the discrete system of equations becomes: Here, M is a divergence matrix for the discrete continuity equation, C u n d the convective matrix, D a diffusive matrix and G a gradient matrix.V is a diagonal matrix that contains the control volume size, p d is a vector containing the discrete pressures and F d is a vector that accounts for the discrete external force.For the convective term this notation has been chosen to show that it is a nonlinear term and that elements of the vector with the discrete velocities u d have been used to construct the matrix.Now, the predictor velocity ˜ u d is introduced.This auxiliary vector will contain the contributions of convection, diffusion and external forcing at the old time level: With the predictor velocity, the discrete momentum equation becomes: The momentum equation is substituted into continuity Eq. ( 35) .The pressure vector on the new time level t n +1 remains on the left-hand side of the equation.The predictor velocity is shifted to the right-hand side of the equation.With the property that G = −M T we observe that a discrete Poisson equation for the pressure is obtained: When the pressure vector at the new time level has been resolved, the velocity vector at the new time level u n +1 d can be found from Eq. (38) .
ComFLOW solves for velocities and pressures and, therefore, the GABC also needs to be expressed in terms of the same variables.The velocity at the boundary u b is the derivative of the potential in the direction normal to the boundary.The pressure at the boundary p b follows from the linear momentum equation in (6) .
For notation purposes, p b is defined as p / ρ.The GABC in terms of pressure and velocity at the domain boundary then becomes: where z c is the vertical position of the center of the cell.
In ComFLOW , the solution variables are staggered within a cell, see Fig. 7 for the definition of the variables at the boundary.Cell labeling is used to indicate completely filled cells inside the domain (F) and mirror cells beyond the boundary (O).The index i is used to refer to cells in horizontal n -direction.The domain boundary is chosen such that it coincides with the position of the horizontal velocity u b .
It is essential that the velocity and pressure in the GABC are defined at the same position to prevent phase differences that induce reflection.The pressure at the boundary p b is obtained from linear interpolation between the pressures on either side of the boundary: (41) It is equally essential that the pressure and the velocity at the boundary are defined at the same point in time to prevent reflection.For a boundary condition in terms of pressures, the pressure at the boundary is determined at time t n +1 .The horizontal velocity at the new time level u n +1 b can be eliminated by means of the momentum equation at the boundary (38) .
The second derivatives in vertical direction in the are approximated by the operator Z.The operator has been derived for a stretched grid in the vertical direction and is second order accurate.The boundary condition with the approximated vertical derivates is set up an equation for p i,k in a mirror cell outside the domain.It relates three velocities with two times three pressures on either side of the domain boundary.In this way, a discrete equation for the absorbing boundary condition is obtained: in which the matrix coefficients for the pressure are equal to: and the coefficients for the velocities are equal to: In the relations above, these coefficients were used: The boundary condition features only pressures at the new time level t n +1 on the left hand side, and horizontal velocities that include convective and diffusive terms on the old time level t n on the right hand The structure of the discrete GABC is similar to that of Poisson Eq. (39) .
Near the bottom and near free surface, the vertical second derivatives cannot be continued because the required information is not available.Near the bottom, we can use (9) to see that the vertical derivative of velocity and pressure are zero near the bottom at z = −h .Near the free surface such information is not available, see Fig. 8 .In this figure, cell labeling has been used to indicate the cell type.In empty cells (E) no equations are solved, in fluid cells (F) both governing equations are solved and surface cells (S) are partially filled cells in between E-cells and F-cells in which boundary conditions are solved.Hence, no pressures and velocities are determined above the free surface.And a mirror condition, such as that near the bottom is not available.Since the largest velocity variation occurs near the free surface, it seems important for a properly functioning GABC to accurately approximate the second derivative there.
The most obvious solution is to derive an operator for a onesided second derivative.The operator only uses solution variables from below the free surface.Unfortunately, applying a one-sided operator for the second derivative near the free surface resulted in unstable simulations.Instead, a compromise was implemented.In the cell containing the free surface, no second derivative is calculated.In the O-cell nearest to the S-cell inside the domain, an ordinary Sommerfeld equation, such as (13) , is solved with a well chosen value for c 0 , typically the phase velocity that is associated with the peak frequency of the spectrum.This equation is implemented at the cost of accuracy, but at the moment an alternative is not available.

Reflection in irregular wave simulations
The performance of the absorbing boundary condition is tested by means of simulations.The performance is better when there is less reflection.The reflection coefficient is obtained from irregular wave simulations; firstly, with the GABC at the downstream end of the domain for the situation without incoming waves at that end; and secondly, with the GABC at the upstream end of the domain (left in Fig. 9 c) for the situation that incoming waves enter the domain while outgoing waves are absorbed.Reflection coefficients from the simulations are compared to the theoretical reflection coefficient.
There are ways to determine the propagation direction and the frequency content from a set of wave signals, measured at multiple locations [25] .These methods, however, are based on linear theory and perform worse as waves become steeper.Furthermore, they will wrongly attribute numerical effects, such as phase lagging or wave energy dissipation to propagating waves.
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.lation with waves propagating in positive x -direction is performed in an infinitely long domain.'Infinite', in this sense, means that the domain length is chosen such, that during the entire measurement time, reflected waves cannot reach the measurement location.Then another simulation is performed.This simulation is the same as the previous simulation in every respect, except the domain length.The domain is short with a GABC at the outflow end.Registrations of the surface elevation, taken at the same positions, are compared to measurements on the infinite domain.Everything being the same, the difference can only be attributed to the boundary procedure.The infinite domain and the shorter domain with the GABC are shown in Fig. 9 ; the figure also indicates the measurement positions in the middle of the short domain.
One simulation in a short domain is performed with the GABC on the side where only wave absorption takes place (see Fig. 9 b).Another simulation is performed with the GABC on the side where generation and absorption take place (see Fig. 9 c).In the latter simulation, the opposing boundary is fully reflecting to make waves go in two directions.Velocities, pressures and the surface elevation from a column in the infinite domain are used as undisturbed incoming characteristic for use in the boundary condition.
The simulation with the GABC on the side where only absorption takes place is also compared to a simulation with a dissipation zone.A dissipation zone was implemented as in Westhuis [27] , which was further studied in ComFLOW by Meskers [19] .It works by adding a term, a momentum damping function, to the pressure boundary condition at the free surface that counteracts the vertical velocity: Here, p is the pressure at the free surface, β is a coefficient that determines the linear increase of the damping function from a position x where the start of the dissipation zone has been defined and w is the vertical velocity at the free surface.The slope β of the damping function is typically low as to reduce reflection by the damping function itself as much as possible.A mass damping function was not applied.
The setup of the simulations is outlined in Tables 1 and 2 .The domain is 2D, which means that only long-crested waves are considered.The domain sizes and the grid distances are stated in Table 1 .The cells have a uniform size in the horizontal direction and are stretched in the vertical direction.In order to be sure that not only the measurement position in space, but also the measurement positions in time are the same throughout all simulations, the time increment t is kept fixed and small.
The simulations are started from rest: at t 0 there are no waves in the domain and the velocities u are zero.Waves are imposed on the left of the domain by means of linear potential theory.With a linear ramp function the signals of surface elevation and velocities at the inflow boundary are gradually built up over a period of 20 s.The irregular wave signal consists of a superposition of regular wave components that accord with a realistic JONSWAP spectrum with significant wave height H m 0 = 4 m and peak period T p = 15 s, see Holthuijsen [17] .
The coefficients of the outflow boundary conditions are in Table 2 .The coefficients for the GABC are tuned in such a way that the reflection coefficient over the range kh ∈ < 0, 6] is never larger than 2%.The Sommerfeld boundary condition at the free surface has only coefficient η to tune.The best choice for η is the value of the phase velocity associated with the peak period of the spectrum.
Meskers [19] gives the optimal configuration of a dissipation zone when linear theory is assumed.The length of the dissipation zone, expressed in number of wave lengths, will be determined based on the wave length associated with the peak period of the spectrum.The optimal number of wave lengths for this simulation, when a theoretical reflection coefficient of 2% is desired, is two.The slope of the damping function β is based on the peak frequency.
The wave signal at the measurement location in the infinite domain 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 Fourier components.The Fourier components are then converted to both spectra and reflection coefficients.This is shown in Figs.10-12 .Note that the vertical axes of these figures are different.
The results with the dissipation zone in Fig. 10 were unexpected and larger than the theoretical value.This is due to nonlinear effects in the dissipation zone that have transferred energy from higher to lower wave frequencies.Dissipation zones tend not to work well for longer waves with lower frequencies because the vertical velocities are small, which explains the large reflection coefficient at lower frequencies.The dissipation zone does perform well for the short waves with kh -values larger than 3.
The reflection coefficient and the reflection spectrum for the GABC are shown in Fig. 11 .In the range from kh = 0 to kh = 5 , the reflection coefficient of around 5% is higher than the 2% it was designed for.Investigation of the reflection coefficients by means of two sets of regular wave simulations (shown in the inset of Fig. 11 a), one with a maximum steepness H / L of 1 • 10 −2 , the other with a maximum steepness of 5 • 10 −2 , leads to the conclusion that the higher reflection coefficient is mostly due to the Sommerfeld boundary condition imposed in the S-cell at the free surface near the boundary and not due to inaccuracies in the Fourier transform, nor due to nonlinear effects.The Sommerfeld boundary condition affects the shortest components of the irregular wave train, above kh = 5 , most.Another reason why the GABC does not match the theoretical reflection coefficient for the larger kh -values is because the shorter components in the irregular waves are resolved by fewer cells per wave length.
There is a difference between the GABC in Fig. 11 and the dissipation zone in Fig. 10 in the value of the reflection coefficient at kh = 0 .Imposing waves by means of Dirichlet velocity boundary conditions is known to add water to the domain.This is apparent from the non-zero value of the reflection spectrum of the dissipation zone at kh = 0 .The GABC makes the additional water flow out of the domain because it is formulated in a way that, on average, makes the free surface tend to the mean free surface.For that reason, a systematic increase of the water level is not observed with the GABC.Fig. 12 shows the reflection coefficient and the reflection spectrum when the GABC on the upstream side of a short domain is both generating and absorbing.The opposing boundary was fully reflective so that the wave signal at the measurement location in the middle of the domain contains wave components going in downstream and in upstream direction by design.The wave signal from a simulation in a long domain with the incoming boundary 'infinitely' far away from the fully reflective boundary is subtracted from the wave signal in the short domain.Because of this setup the wave signal contains mainly the longer wave components up to kh = 1 .2 , because they propagate faster and arrive earlier at the reflective boundary.The difference in wave signal from the long domain and the short domain is attributed to spurious reflection.From Fig. 12 we find that the reflection spectrum and the associated reflection coefficient of the implemented GABC are somewhat higher than the theoretical reflection.
The GABC, both on the right when only absorbing and on the left when generating and absorbing, behaves as expected, but with more reflection.The amount of wave energy that is reflected is small compared to the input spectrum and at a comparable level as spilling beaches in wave tanks used for experiments.There are two main reasons for the difference between theory and implementation.The first reason is the Sommerfeld condition in the SO-cell at the boundary, see Fig. 8 .This relation was tuned to the peak frequency of the spectrum and will especially cause reflection of the longer and shorter components away from the tuning frequency.The second reason is that nonlinear effects near the boundary are not accounted for in the design of the GABC.

Practical example of a validated simulation
We mean to show a real application of the GABC for a situation with unidirectional incoming waves where a simulation with dissipation zones on the upstream and the downstream side that are The experiment, with identification number 202006, was performed in a long wave tank.Waves were generated with a pivoting wave board at one end of the tank.At the opposing end a spilling beach was present to induce wave breaking and reduce reflection.The model was placed a large distance away from the wave board.
A simulation was performed to compare to the experiment.A better comparison with the experiment is obtained when nonlinear wave kinematics are used for input at the boundary instead of linearized wave kinematics from something like a wave gauge.To obtain the nonlinear kinematics, measurements of the rotation of the wave board were taken and fed into a nonlinear potential flow model [27] with a domain that was so long that wave reflection from the downstream end did not occur.The nonlinear potential flow model propagates waves at a much reduced computational cost compared to ComFLOW .Kinematics at the position of the inflow boundary of ComFLOW were taken from the nonlinear potential flow model.
The inflow boundary in ComFLOW is located at x = 5156 m with respect to the wave board.This is where the GABC is used together with the kinematics from the nonlinear potential flow model as the incoming wave characteristic.The outflow boundary with a GABC is positioned at x = 5456 m. 300 cells were used in xdirection without stretching.In the vertical direction 75 cells were applied with cell size stretching in such a way that every n -th cell away from the cell with the minimum cell size is bigger according to: h n = (1 + γ ) n h min (44) with γ a stretching coefficient and n a counter of how many cells are between it and the cell with grid size h min at the focus point.
In z -direction, γ = 5% and the focus point of the stretching is positioned at the mean free surface level z = 0 .Near the free surface z min is 0.5 m, near the bottom z is nearly 14 m.In y -direction the boundaries of the domain were chosen to coincide with the tank walls, each positioned 100 m away from the center of the structure at y = 0 ; 75 cells are used with a stretching coefficient of 2%, y min being 1.5 m and the largest cells near the y -boundary slightly over 2 m.The simulation time is 120 s The center position of the structure was at x = 5306 m.The GABC was used at the upstream boundary and at the downstream boundary without incoming waves.Fig. 13 b shows a snapshot of the 3D simulation.
Wave probe 2 was positioned at x = 5238 m on the upstream side of the structure in the middle of the tank.Fig. 14 a shows the surface elevation as a function of time for wave probe 2 and compares it to ComFLOW simulations with the GABC, a dissipation zone and the Sommerfeld boundary condition with c 0 = 0 .56 gh .There is good agreement between the GABC simulation, the simulation with a dissipation zone and the experiment with the surface elevation in the GABC simulation marginally higher than in the experiment for the highest waves.The simulation with the Sommerfeld boundary condition gives different results from t = 650 s onwards as a result of spurious reflection.
A similar result is found for the impacts in the registration of pressure sensor P12 on the front column of the structure (positioned at x = 5266 m and at z = −10 m with respect to the mean water surface), which is shown in Fig. 14 b.There is overall good agreement between the pressure in ComFLOW and in the experiment for the simulation with the GABC and with the dissipation zone.The simulation with the Sommerfeld boundary condition demonstrates that spurious reflection can influence the impact pressures.

Conclusion
This article discusses the background, derivation, implementation and results of an open boundary condition for dispersive freesurface waves.We call it a Generating Absorbing Boundary Condition or GABC.The starting point for the derivation of the GABC is the Sommerfeld condition, which is perfectly absorbing for one wave component with one propagation velocity.This velocity is specified in the Sommerfeld condition by means of a tuning parameter.It was found that the range of absorbed wave components can be extended by replacing the tuning parameter with an approximation of the linear dispersion relation in terms of the wave number.When second derivatives of the solution variables in the vertical direction are substituted for the wave number, it yields a GABC with a low 2% reflection for a range of components instead of a single component.Stability issues loom when the GABC is implemented, but the mechanism, by which instabilities occur, is well understood and stability criteria have been formulated.The stability criteria hardly restrain the absorbing performance of the GABC.
In irregular wave simulations, the GABC has similar performance to a dissipation zone that was two typical wave lengths long in the range of wave numbers of interest.It was found in numerical 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 was derived from theory, but comparable to the amount of reflection that is said to be obtained in experimental basins and flumes.The difference between theory and practice here is mainly due to the implementation near the free surface.
It was demonstrated that the GABC can be used in an actual 3D simulation of steep wave interaction with floating structures.The simulation was validated with experimental results.The simulation with the GABC has similar performance to the simulation with the dissipation zone, at a much reduced computational effort.
values of the local height function.Near domain boundaries it is convenient to use an axis system with a component n normal, and a component τ tangential to boundary B, see Fig. 1 .The side walls of the domain are vertical and z is directed upward with respect to the n − τ -plane.The direction of impinging waves on the boundary will be expressed by angle θ with respect the normal direction n to the boundary.

Fig. 2 .
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.

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

Fig. 4 .
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, respectively.

Functions f 1
and f 2 in (33) are plotted in Fig. 5 .The figure shows only the real values of f 1 .It has imaginary values in the range kh ∈ [ π /2, π ].The range of imaginary values gives us the opportunity to ensure stability.If the roots of the functions f 2 and f 3 in Eq. (

Fig. 7 .
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.

Fig. 8 .
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.

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

Fig. 11 .
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 compares 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.

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

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

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

Table 1
Domain and grid.

Table 2
Coefficients of the boundary procedures.