Dynamic and kinematic characterization of the impulsive wavemaker system in a numerical wave tank

We study the dynamical response of a piston-type wavemaker in a numerical wave tank. The two-dimensional, fully viscous unsteady Navier–Stokes equations are solved on a two-phase flow configuration using the volume of fluid method to capture the free surface dynamics. The wavemaker is a moving wall driven by an arbitrary signal waveform. The step response of the wavemaker may generate pulse-like waves similar to an undular bore propagating along the tank. Wave elevation at the piston wall has close similarity to the time response of second order systems found in feedback theory. The scaling found for water elevation at the piston wall for different step velocities and mean still water levels is in agreement with that in the available theory at low Froude numbers. The results along the tank for continuous waves agree with those of potential theory. The power input during the step response was determined during the whole wave generation process showing that net piston forces are predominantly hydrostatic. A power scaling for different mean still water levels and step velocities as a function of the Froude number was obtained. An active absorption strategy based upon a feedback controller driving a secondary piston was implemented. Wave absorption was successfully achieved on regular and irregular waves. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0017376


I. INTRODUCTION
Wave tanks are centerpieces when it comes to studying hydrodynamics and wave structure interaction in offshore and marine engineering. Experimental or numerical, they can reproduce the offshore environmental conditions (waves, currents, winds for some of them) in order to evaluate the dynamic response to these harsh conditions of vessels, offshore platforms, and marine renewable energy structures such as fixed/floating offshore windmills or wave energy converters. For a long time, these phenomena could only be experimentally studied. The wave tank scaling should be done, according to the similarity principle of dimensional analysis, by matching relevant numbers as the Froude and the Reynolds numbers. This is often complicated especially at small scale as explained in Ref. 1. One approach is to test the model through increasing model scales. 2 The most simple configuration for a wave tank is a 2D tank, which is a long tank where waves are generated on one side and absorbed or broken up on the opposite side.
In order to generate waves at a laboratory scale, one should displace a certain amount of water by means of a wavemaker of choice: piston, hinged paddle, double articulated paddle, plunger, duck. First order wavemaker theory has first been studied in Ref. 3. Wave generation alternatives have been theoretically studied in the 1950s 4 including piston, hinged paddle, or double articulated paddle. In this work, we choose to use a piston wavemaker where the piston stroke along the tank is the physical mechanism to create waves. The resulting waves from a piston wavemaker were investigated in Ref. 4 determining a transfer function or relationship between the piston stroke and the wave height. A validation of the piston wavemaker theory was provided by some experiments. 5 Similarly, Madsen 6 developed a theory of wavemaker generation for pistons providing a detailed expression for wave elevation at any position in a wave tank, and Schäffer 7 developed the second-order wavemaker theory for irregular waves. Goring 8 studied solitary wave generation with a piston wavemaker, basing his solution on the theory of Ref. 9. The impulsive wavemaker was experimentally and numerically studied in Ref. 10, notably showing wave profiles from the experiment as well as forces on the wavemaker.
Most of the research on wavemaker theory has been conducted on the basis of potential theory avoiding the determination of viscous effects. At the same time, the necessity for local solutions at the piston wavemaker was considered important in order to understand the physics of the first instants in the wave generation. Local solutions near the piston wall are provided by Joo et al. 11 using the Laplace equation and expansion series for small Froude numbers and basing their work on Ref. 12. They investigated the contact line motion at initial times of the piston wavemaker for ramp, step, exponential, and harmonic velocities. They later on extended their studies to stratified fluid. 13 A study of the steady motion of a ship bow at high Froude number, which is similar to the piston wavemaker problem, was carried out in Ref. 14. Experimental realization of the velocity step case becomes complicated because imposing a prescribed step-like wavemaker motion involves high piston accelerations and therefore peak power inputs. The practical case is to perform a ramp motion as indicated in Ref. 15.
While experimental wave tanks provide a physical platform to implement such systems, testing may be expensive or not suitable to comply with similarity restrictions, whereas numerical studies once performed exclusively with potential codes (irrotational, incompressible, and nonviscous flows) are now commonly carried out taking into account viscous effects that, nevertheless, require more computational resources.
At present, the use of efficient codes solving viscous Navier-Stokes equations has been reported to be useful in order to model extreme wave conditions where the wave breaking process 16 can be captured and identified as the dominant mechanism (instabilities) involved in extreme loading on marine structures and wave energy converter devices. 17 Numerical wave tanks may recreate realistic waves using moving boundaries (which involves a dynamic mesh) or static boundaries with mathematical implementation (in the domain or at the boundaries). In this work, we use the first option that mimics a viscous piston motion in a physical experiment. 18,19 Most of the numerical wave tank models include the volume of fluid (VOF) to simulate two-phase flows 20 as well as the wave generation and propagation problem. [21][22][23][24] Wave tank modeling with active wave absorption was developed in OpenFOAM, 18,25 and recent methods such as Smoothed-Particle Hydrodynamics (SPH) have also been developed to represent two-phase flows. 26,27 In this work, we present a numerical study of a piston wavemaker at a laboratory scale to investigate the wave generation process and characterize the piston wavemaker dynamics, extending Froude number regimes far beyond the ones that were studied before. We perform fully non-linear and viscous numerical simulations of the wave tank where surface waves are produced by the prescribed motion of a piston wavemaker. The moving boundaries allow reproducing the physics of a real facility. A complete mesh independence study is carried out, and validation is performed with theoretical data. To be able to generate waves of high quality, it is necessary to characterize the wavemaker system. It is chosen here to do so by basing the study on the initial-value problem, theoretically studied in Ref. 11. The numerical simulation allows us to study regimes that do not appear in the theoretical study based on potential flow. We measure the water height at the wavemaker and in the tank for a set of prescribed wavemaker motions, allowing one to characterize the near and far field generated waves and determine the power input needs according to the wavemaker dynamics.
The structure of this article is as follows: the initial-value problem and the required equations are presented in Sec. II, the numerical method, mesh independence tests, and validation are presented in Sec. II C, and finally the resulting characterization of the piston wavemaker and the power input analysis are presented in Sec. III.

II. PROBLEM FORMULATION
The problem is described in the schematics in Fig. 1. We consider a 2D wave tank of length L and height d, equipped with a moving piston wavemaker placed at the left wall, and the opposite wall of the tank is situated on the right. The mean still water level is denoted h and can be varied at will. The origin of the coordinate system is located at the left bottom of the wave tank, and all different measurement stations along the wave tank are referred to this coordinate system. In order to characterize the piston wavemaker and the wave tank behavior, a series of tests will be performed to provide the system response to input velocity signals driving the motion of the piston wavemaker. The overall wave tank response will be obtained from Heaviside Θ(t) input signals representing a step response test. 28 We consider the system composed by the wavemaker, whose input signal is its velocity UG and the response is the water level at the piston wavemaker η w , as shown in Fig. 1. We can define a transfer function associated with the wavemaker system, hereafter called G where the fundamental output/input relationship (transfer function) can be expressed as G = η w /UG, which is crucial when implementing feedback controlled wave-absorbers. 29

A. Governing equations
The numerical simulation solves the 2D two-phase laminar Navier-Stokes equations with two incompressible fluids (water and air phases). The mass conservation equation is also solved in primitive variables incorporating the VOF model to deal with each fluid phase. 20 These equations can be written in the vector form as where U is the velocity vector, p * is the pseudo-dynamic pressure (p * = p − ρg ⋅ x), ρ is the density, μ is the dynamic viscosity, g is the gravity acceleration, x is the position vector, and σ is the fluid surface tension coefficient. The volume fraction α is introduced to deal with the two-phase formulation within the volume of fluid (VOF) framework. κ is defined as follows: In the two-phase formulation, density and viscosity on each domain cell are computed as a weighted mean of the form 30,31 The relevant dimensionless number in this work is the Froude number Fr = UG/ √ gh controlling the wave propagation dynamics. Viscous effects become important during the interaction of the starting wave with the piston wavemaker taking place at the beginning of the wave generation as we will discuss in Sec. III C.

B. Boundary conditions
Boundary conditions employed in this work are summarized in Eq. (7). On the piston wavemaker wall, we impose no-slip condition for all velocity components. The initial still water level h is always established before any wavemaker motion by initializing α. The volume fraction α is bounded and may adopt any value within 0 ≤ α ≤ 1 in any place of the physical domain, and Neumann boundary conditions, set to 0, are applied at all boundaries for the α variable as well.
In particular, at the wavemaker, such a condition forces the contact angle of the interface to be perpendicular to the wavemaker wall, which agrees with the experimental results of Refs. 15 and 32 for the case of a continuously accelerated wavemaker. In this work, the piston wavemaker moves according to an input velocity signal displayed in Fig. 3(a). A velocity step UG drives the piston resulting in linear dependency on time displacement XG(t). At the piston wavemaker wall, the velocity matches the velocity of the moving boundary in the x-direction only. A no-slip condition is used at the seabed wall and the right end wall. The pressure is set to a reference pressure (in this case 0) at the atmosphere boundary, and 0-Neumann conditions are used at the other locations. A zero-gradient condition is applied at the atmosphere for outflow, and a velocity u ϕ is assigned for inflow based on the flux in the patch-normal direction. The boundary conditions for the velocity, pressure, and α variables can be summarized as C. Numerical method The governing equations were solved with the open-source software OpenFOAM version 5. 33 OpenFOAM is an object oriented C++ toolbox for solving continuum mechanics problems with the finite volume method. OpenFOAM presents many advantages: as released under the GNU General Public License, it is free and open-source (no licensing fees, unlimited number of jobs, users, and cores). It is also largely used in the scientific community and thus has been validated for many applications. We use interDyM-Foam, a solver of the Navier-Stokes equations for two incompressible isothermal immiscible phases using the volume of fluid (VOF) method. It increases the capabilities of previous solvers allowing one to handle dynamic mesh motion. OpenFOAM solves a single equation of momentum for the two-phase mixture by introducing a volume fraction advection equation of the VOF method used to capture the interface between the phases. Hirt and Nichols 20 presented this method as an efficient and simple way of treating the free surface in numerical simulations, as it stores a minimum amount of information. This method should be carefully used when the surface tension becomes important. Some numerical solvers such as interDyMFoam impose some restrictions in order to keep a sharp interface between both fluid phases. An additional term called artificial compression is introduced here, 31 Uc = min(Cα|U|, max(|U|)).
Cα is a user defined coefficient whose default value is 1. The additional term is only active close to the interface because of the product α(1 − α) and does not impact the solution outside the interface ARTICLE scitation.org/journal/adv region. Its role is to compress the interface and maintain α between 0 and 1 if used with discretization techniques. In the post-processing stage, the value of α = 0.5 is chosen to detect the free surface, which is carried out thanks to linear interpolation. The interDyMFoam solver uses the PIMPLE algorithm, which combines both SIMPLE (Semi-Implicit Method for Pressure Linked Equations) 34 and PISO (Pressure Implicit Split Operator) 35 algorithms and allows for bigger time steps. Simulations were performed on a CPU Xeon E5-2660 v2 cluster running on Simple Linux Utility for Resource Management (Slurm) and based on MPI libraries. The CPU run time for a one second transient simulation and a typical wave was about 2.4 h for a 400 000 element mesh and Δt = 5 × 10 −4 s time step. The geometry of the wave tank required a fine spatial discretization, particularly in zones such as the water-air interface and zones of high velocity gradients such as the wavemaker walls. Explicit schemes are used so that special care is taken when choosing the time step regarding the mesh size and CFL condition below 1. As a result of the rapid input velocity signals during the wave formation and subsequent progression along the wave tank, the temporal discretization requires time steps smaller than 10 −3 s; thus, we choose a time step size Δt = 5 × 10 −4 s for all the simulations.

D. Mesh independence tests
A series of mesh tests are performed to look for mesh independent results. We consider two test cases in order to achieve mesh convergence: (i) the study of the response to an input velocity step characterized by the overshoot and the steady state water height, and the rising, peak, and settling times, as found on the time response of linear dynamical systems, and (ii) the wave propagation of linear waves along the 2D wave tank (wave height, wavelength). In the first case, a set of uniform rectangular cell meshes are generated to look for mesh independence. Most sophisticated meshes, with a uniform zone at the water-air interface, which includes the minimum and the maximum water height during the whole simulation, are also set up. They include refinement at the wavemakers and are denominated "non-uniform." They allow decreasing the computation time as the mesh size is kept in reasonable limits. The mesh, generated thanks to the blockMesh program, is built from three characteristic elements: the cell width at the piston wall Δxw, the cell width in the wave propagation zone (far from the walls) Δx, and the cell height in the wave propagation zone Δz. Δxw will be determined from the uniform mesh study in (i). From the wavemaker, the longitudinal element size Δxi is being computed with the following geometric law: Δxi = Δxwr (i−1) ∀i < nj, where r is the geometric growth rate. After nj cells, the elements reach a constant size Δx in the wave propagation zone. A similar calculation allows us to define the vertical element size from outside the interface zone. Both Δx and Δz are set up considering the number of cells per wavelength and cells per wave height, respectively. These parameters will be set in (ii). Figure 2 shows the mesh in the x-z plane, while one cell size is set up in the y-direction. Mesh type and properties are displayed in Table I. These preliminary tests with different meshes allow an appropriate choice of the mesh size without compromising accuracy and CPU time (cf. Fig. 2). The dynamic mesh is modeled using a mesh expansion/contraction strategy. The mesh uniformly contracts and expands, conserving the global mesh cell volumes as these motions are relatively small compared to the tank length.

Step response
In order to fully test the capacity of the code to represent sudden and rapid water surface elevation, we choose to perform a step response of the water tank. The input signal is a velocity step, as TABLE I. Mesh properties to study the wavemaker response to a velocity step. Δxw is the cell size in the x-direction at the piston wavemaker, Δx is the cell size in the x-direction in the wave propagation zone, Δz is the cell size in the z-direction at the water-air interface, η o is the overshoot water height at the piston wavemaker, η ss is its steady state value, and tr , tp, ts are the rising time, the peak time, and the settling time, respectively. The area a is used as an entry for the GCI study. 36 shown in Fig. 3(a), and given by where Θ(t) is the Heaviside function. The step value is 0 for negative times and UG = U 0 for positive times, resulting in a linear displacement XG(t) of the piston wavemaker. In this part of the study, we set the tank length to L = 4 m, the tank height to 0.25 m, and the mean still water level to h = 0.15 m. Mesh M 1 U is the coarsest and M 6 U the finest. M 7 NU to M 9 NU are non-uniform meshes as previously described. They use a cell size at the wall defined later on in the conclusion of the uniform mesh study and geometric growth rate r = 1.05 in the x-direction and r = 1.2 in the z-direction. For the finest mesh in the z-direction (mesh M 6 U ), it is necessary to reduce the time step to keep the Courant number below 1. This is why the time step is set for all meshes to Δt = 10 −4 s. The measured quantity is the water elevation at the wavemaker η w . The variables of interest are described in Fig. 3(b) and are the following: η o , η ss , the overshoot and steady state values, respectively, and tr, tp, ts, the rising, peak, and settling times, respectively. Results are reported in Table I and shown in Fig. 4. Figure 4(a) shows the time series of the water elevation at the piston wavemaker, η w (t), for each mesh. The results of the four finest meshes are similar. A Grid Convergence Index (GCI) study 36 is carried out with meshes 4, 5, and 6. The local order of accuracy p ranges from 0.18 to 17.57 with a global average of 4.98. This apparent average order is used to assess the GCI error at every point as suggested in Ref. 36. The error made in the last mesh is really low (the maximum GCI error is 0.3%), which allows us to say that the results do not depend on the mesh. Figure 4(b) shows the error for every point, and a zoom around the overshoot is displayed in the inset where errors increase. The mean error is an order of magnitude lower (0.03%). The errors on the overshoot and the steady state values are very low. Figures 4(c) and 4(d) show the convergence of the different parameters composing the typical response to the step (overshoot and steady state water elevation, rising, peak, and settling times). They show that for these parameters, the limit of convergence is mesh M 3 U . It is chosen, for obvious practical reasons, to work with mesh M 3 U at the wavemaker, which allows us to decrease computational times while assuring convergence. Nonuniform mesh shows good agreement with the uniform ones for all variables. In the rest of this work, we make sure that the first cell at the piston wavemaker wall is kept below Δxw = 0.0025 m in the x-direction and Δz = 0.001 m in the z-direction in order to keep the results independent of the mesh. The number of cells per wavelength and height in the wave propagation zone is analyzed in the following paragraph.

Wave propagation
In order to properly study the wave propagation properties as a function of the mesh type and quality, we use an extended wave tank with L = 8 m. The mesh is finer at the water-air interface and is kept uniform in the zone where the wave propagates. At the wavemaker, we set Δxw = 0.001 m, and a transition is made with a 1.05 cell to cell ratio. The piston stroke is set to X 0 = 0.004 m, the wavemaker frequency to f = 1.25 Hz, and the piston wavemaker velocity to UG(t) = X 0 ω/2 sin(ωt + δ) with δ = −π/2. The corresponding wavelength can be estimated from the dispersion relation ω 2 = gk tanh(kh), where k is the wave number, and in this case, λ = 0.82 m. A common discretization is given by 20-25 elements per wave height and 60-70 elements per wavelength in recent works. 37 We conduct our test based on three meshes whose properties are shown in Table II, where M 7 NU is the coarsest mesh and M 9 NU the finest one. The number of cells per wave height ranges from 15 to 60, while the number of cells per wavelength ranges from 60 to 240. The time step is set to Δt = 0.001 s, and the theoretical CFL number is reported. Even if it shows to have a value below 1, the time step for the finest mesh M 9 NU had to be decreased to Δt = 0.0005 s to avoid divergence due to the use of explicit schemes. The simulation end time is 10 s, and two probes are set at x = 2 m and x = 4 m from the origin of the coordinate system (cf. Fig. 1).  Figure 5(a) shows the water level at the wavemaker wall η w (t). No differences between the meshes are observed as expected, since the refinement in the x-direction and in the z-direction is finer than the necessary one studied in (i). Figure 5(b) shows the differences between the meshes that are used at x = 4 m. A zoom over the highest value and for three wave periods is shown in Fig. 5(c). The coarse mesh M 7 NU effectively produces minor differences from the other two, specially at the maxima and minima. The results for the probe at x = 2 m are compared with the theory in Ref. 6 and shown in Fig. 5(d). The errors of wave crests and troughs with the wavemaker theory are shown in Table III. The results are quite accurate for the three meshes (almost all cases with a rms deviation below 0.1 mm) although the coarse mesh lacks accuracy at the maxima and minima. It is chosen to work with the medium mesh M 8 NU in the rest of this work as it allows accurate wave height data and reduces the computational time compared to the fine mesh. The final mesh characteristics are given in Table IV.

III. RESULTS
One of the objectives in this work is to characterize the piston wavemaker system by applying a series of velocity steps. First, a qualitative approach is taken and observations are made about the wave generation at small times and the generated wave pulse is studied at longer times and far away from the wavemaker. Then, we compare the characteristic variables with the theory developed in Ref. 11 and explore higher Froude number regimes. The forces exerted on the wavemaker and the power involved in the wave generation process are studied, and an active wave absorption strategy is finally designed. The step velocity tests are carried out at four mean still water levels: h = 0.050 m, 0.075 m, 0.100 m, 0.150 m, and for velocities ranging from UG = 0.005 m/s to UG = 0.4 m/s. In this problem, the fundamental velocity, length, and timescales are √ gh, h, and √ h/g, respectively. The problem can be written as f (ηw, ρ, g, h, t, UG) = 0, but according to the Buckingham π theorem, 38  t * are the dimensionless water height at the piston wavemaker and time, respectively. In the next sections, as a consequence of this analysis, we will express the results in a dimensionless way using the relevant variables η * w , t * , and Fr.

A. First instants-The overshoot-wave
While applying a velocity step to the wavemaker, one observes the formation of a surface water pulse propagating along the tank as it is shown in Figs. 6 and 7. The wavemaker displaces a given volume of water, and as the fluid is incompressible, this volume is found under the wave profile. The water height at the wavemaker η w (t) first rises along the wall reaching a maximum value, the overshoot. Then, the wave pulse comes off the wall and self-propagates in the positive x-direction, as the phase velocity of the wave becomes higher than the wavemaker velocity Cp > UG. In the time snapshots in Fig. 6, the wave pulse generated by the velocity step is called the overshoot-wave as it is created by the overshoot of the water elevation at the piston wall.

ARTICLE scitation.org/journal/adv
After leaving the wavemaker, the water height at the wavemaker remains constant, as indicated in Fig. 9. The wave pulse profile for times t = 0.05 s, 0.13 s, 0.20 s, 0.31 s, which correspond to the first instants of the pulse formation, is displayed in Fig. 6. In Fig. 6, the wave steepness, defined as |dη/dx|, is superimposed. We observe that the wave steepest shape occurs at the very first instants of its formation, i.e., at t = 0.05 s, where the maximum steepness is higher than 3. When the pulse starts leaving the piston wavemaker, its steepness decreases to values under 1. Such high steepness is associated with the non-linear properties of the overshoot-wave. Another important parameter useful to evaluate when it comes to apply linear theory is the wave height to mean still water level ratio η/h. It is relatively important (values around 0.5 in this example); therefore, Airy theory of linear waves cannot be applied here. The profile of the generated pulse at longer times is shown in Fig. 7, as well as its steepness as a function of the x-coordinate. We can notice the generation of wiggles after the main pulse propagating downstream as already described. The wiggles are trailing waves whose height decreases in the vicinity of the piston wavemaker, while the front wave height increases as it travels. These wave structures are described in the work of Ref. 11 and are a consequence of wave dispersion. The piston motion creates a wave packet in which each wave travels at different velocities due to dispersion effects. We identify the created wave to be an undular bore, which was experimentally and theoretically studied in Refs. 39 and 40, respectively, while Stoker 41 predicted that an impulsive wavemaker would generate an undular bore as in the present case. Recently, similar bore structures were obtained from a different approach where a moving weir at the bottom of a channel may produce the volumetric water displacements necessary to develop such a wave. 42 Undular bores are of particular importance as they appear to more likely represent a real tsunami wave instead of a solitary wave. 43 Important values of the wave steepness are observed during the pulse formation, which decrease rapidly as the pulse propagates without evidence of wave breaking. If the bore strength is important, undular bores may display a wave breaking process. 42 In the time snapshots in Figs. 6 and 7, we observe the vector velocity field of the numerical solution of the Navier-Stokes equations. The vector field becomes intense precisely near the steepness peaks as the overshoot-wave propagates along the tank. The higher the wave height, the higher the intensity of the vector field, which is associated with particle velocity.
The starting overshoot-wave displays a phase celerity Cp as a function of time t, which may be compared with two characteristic properties, the shallow water wave speed Cp = √ g(h + H) and the piston velocity UG, as shown in Fig. 8. The overshoot-wave celerity may be estimated from the mean celerity as Cp = δx/δt between consecutive wave crests or looking for maximum steepness at different time steps during the propagation, as shown in Fig. 6. As the phase celerity of the overshoot-wave is found to be higher than the piston velocity from the first instants of motion, the overshoot-wave travels fast enough to leave the piston wavemaker. As the wave propagates along the tank, its phase celerity increases with time, approaching the shallow water phase speed given by Cp → √ g(h + H). This behavior agrees with that in potential theory and therefore allows us to validate the capabilities of the numerical wave tank for wave propagation.

B. Response to a velocity step
A positive velocity step (UG > 0) creates a water pulse, which leaves the wavemaker wall as it propagates along the tank. The step response is recorded as the water height or water elevation at the wavemaker wall η w (t) and is presented in Fig. 9(a).
The water height first increases, reaching a maximum or overshoot η o at peak time tp, before approaching a lower steady state value η ss . It is of particular interest to note the similarity of this dynamic response to a step response of the second order system [cf. Figs. 9(a) and 3(b)]. The time response of second order linear systems depends on the type of the input signal. For a step input, the system exhibits a characteristic response defined by the rising time tr, the maximum overshoot η o , and the steady state value η ss obtained at a given settling time ts. 28 In Fig. 9(b), we present the normalized time response using the Froude number as in the theory proposed in Ref. 11. We can observe that after reaching its maximum value, the overshoot, the water height slightly oscillates and decreases to its steady state value η ss = hFr. The step response for small Froude numbers found in this work is in agreement with that in the theoretical work of Ref. 11. Nevertheless, some differences with the theory are observed in the overshoot behavior for higher Froude numbers. The overshoot starts to increase beyond the theoretical prediction, and the steady state value approaches a slightly higher value than the expected one from theory, η ss = hFr, as we will discuss in Sec. III B.
The relative deviation from theory 11 is presented in Table V  deviation decreases with the increase in the water level (which actually corresponds, for a given step velocity, to lower Froude numbers). Finally, the deviation for the overshoot value can reach high values (> 40%) for the highest Froude numbers (high velocity, low mean still water level).
The dominant controlling parameter in this work is the Froude number Fr defined as the ratio between inertial (ρh 2 U 2 G ) and gravitational (ρgh 3 ) forces, and thus, Fr = UG/ √ gh. The similarity of the system response to the typical second order response under a step velocity can be analyzed through the characteristic timescales and maximum overshoot of the response. Timescales as well as the amplitude of the water height response must be plotted as a function of the step velocity UG. As in many feedback control systems, the time response of a second order system will be completely determined when one knows the maximum overshoot η o , the steady state amplitude η ss , and the rising, peak, and settling times tr, tp, ts, respectively. If the time response is correctly described by the second order response, we should find a scaling law for the characteristic timescales at different water depths h and step velocities UG otherwise said the Froude number.
The overshoot η o as a function of the step velocity UG is plotted in Fig. 10(a) and the steady state η ss is plotted in Fig. 10(c), both for two mean still water levels h = 0.05 m and h = 0.15 m. The steady state values evolve linearly with the generation step piston velocity as they correspond to the water height of the displaced volume, which increases linearly in time as the step velocity UG is constant for t > 0. However, the overshoot dependence on UG is not linear as the overshoot-wave originates during a rapid transient process. The curve collapse under the proposed scaling is relatively good in the range of mean still water levels shown as 0.05 ≤ h ≤ 0.15 m. The scaling for the overshoot fails when Fr > 0.2 as the overshoot-wave height increases over the linear limit and starts to move faster with higher UG.  The evolution of the characteristic timescales of the response to a velocity step is shown in Fig. 11. As it is indicated in Fig. 3, these timescales are associated with a second order dynamical response, where the water height at the piston wall η w is measured and compared to the steady state height η ss . The rise time tr and the peak time tp are associated with the very first instants of the wave motion, when the overshoot-wave is created. Both time constants appear to be independent of the step velocity and provide an interesting scaling independent of the Froude number where the dimensionless timescale is written as As shown in Fig. 9, after the main wave leaves the piston zone, the water height on the piston wall scales perfectly with the Froude number as η w /(hFr) = 1 because the displaced volume in the steady state regime is entirely determined by the steady state water elevation η ss . The associated settling time ts is computed within a 10% band and slightly grows with the Froude number. In Fig. 11(b), the time scaling indicated in Eq. (11) produces a very tight collapse of each characteristic time as a function of the Froude number.
The scaling seems to confirm the similarity of the water elevation response with a second order response. As the wave pulse is created by the excess or overshoot of the water elevation at the piston wall, we called it the overshoot-wave. The piston motion produces the displacement of a water volume (per unit depth) given by V(t) = UGth, which displays an initial transient peak, the overshoot, superimposed into a water slug rising over the mean still water level h. The overshoot-wave has its own dynamics, moving at shallow water speed, running over the water slug, and therefore leaving the displaced volume faster than the linear waves.

C. Forces involved in the step response
In this section, we present new findings such as the force decomposition and the maximum power as a function of the Froude number. The objective of this section is to determine the forces originated on the piston wavemaker during the step response. As the problem is 2D, there are only two components of the forces projected on the piston wall in the x-direction and z-direction, which are inertial, pressure, and viscous forces. These forces can be calculated from the pressure and velocity fields at the wavemaker using the stress tensor ε , the elementary surface area dS (which is in the 2D case an elementary length), and its normal n, The stress tensor is defined as In Fig. 12, we display the resulting normal and tangential force pro- force Fz(z) in Fig. 14(b). When the overshoot-wave leaves the piston, and the progressing volume pushed by the piston reaches a steady state motion, the shear forces become very small. When the step velocity is increased to UG = 0.3 m/s, the normal force profile f x (z) displays notorious changes in time but finally evolves into an almost hydrostatic vertical profile, as shown in Fig. 13(b). In order to verify the effects of the initial fluid motion on the pressure field during the step response, we recorded the pressure profiles at both step velocities UG = 0.03 m/s and UG = 0.3 m/s. The sudden increase in water elevation at the piston wall is more important at higher UG, which is explained by the higher volume displaced during the initial times. Shear forces change in sign during the formation of the overshoot-wave, as shown in Fig. 12 at both step velocity values. The first water elevation motion produces a positive shear on the piston wall and therefore a positive shear force, which start to decrease, becoming negative as the overshoot-wave leaves the piston wall and the force points downward before vanishing in the steady state regime. The pressure profile is dominated by hydrostatics, as shown in Fig. 12. However, when we subtract hydrostatic pressure due to the initial wave elevation along the z axis, i.e., we plot p(z) − ρg(h + η), we observe traces of the creation of the overshoot-wave on the remanent pressure. As time progresses and the overshoot-wave leaves the piston wavemaker, this remanent dynamic pressure contribution approaches very low values with respect to hydrostatics.
In order to compute the power delivery involved in the process, we must get a good estimate of the resulting forces on the piston. In Fig. 14, we present the normal and tangential forces on the piston wall as a function of time, uncovering the initial transient associated with the overshoot-wave formation and during the steady state regime. The normal and tangential force profiles are obtained integrating the force profiles in Fig. 12  wall during each time step, providing an accurate estimation of the net forces shown in Fig. 14. The resulting normal forces Fx(t) are negative as they oppose the piston motion. If we plot the absolute values, we find a time evolution very similar to that of water elevation η w (t) in Fig. 9. The first rising part is associated with the force excess resulting from the creation of the overshoot-wave followed by a steady state force associated with the progressive motion of the mass of the water slug moving at constant velocity UG. On the other side, the tangential averaged forces Fz(t) shown in Fig. 14 display a change in sense (sign), indicating how the fluid is moving on a boundary layer created on the piston wall. At first, the fluid moves upward, then stops, and then moves downward reaching a local (downward) maximum precisely when the overshoot-wave leaves the piston. The critical time when the resultant shear force is zero t ∼ 0.2 s does not correspond to the peak time as water is still rising at the bottom of the wavemaker, as shown in Fig. 6. X-force characteristic times (rise and peak times) are much larger than those for the kinematic observations. The minimum shear force time corresponds to the maximum force in the x-direction, showing the correlation between both phenomena. We can imagine that a feedback controlled piston wavemaker might also be designed by measuring the vertical force instead of the one in the direction of the piston motion. If the overshoot-wave is going away leading to a fluid flow downward at the interface, the global motion is more complicated, with a flow still going upward at the base of the piston.
During the design of a piston wavemaker, it is important to evaluate the power input during the step response associated with the wave generation process, as it can be particularly useful in determining the scaling. For example, the Flowave facility power demand can creep close to 300 kW when it creates a sea state moving the 168 paddles. 44 Dimensioning the necessary power supply is then of crucial importance. This section is devoted to the evaluation of the energy input required to create a wave pulse resulting from a piston velocity step. The power involved in the step motion of the piston wavemaker is calculated according to the following equation: In this case, the piston wall velocity is Uw = UGx, which is the Heaviside step function. Maximum power delivery as a function of the step velocity is shown in Fig. 15(a). The maximum power corresponds to the maximum water elevation at the piston wall. An expected power increase is found when the step velocity increases, but more impressive is the radical change in the power delivery when the mean still water level is increased. If we look for a scaling of the power delivery, we may use a characteristic force per unit length and velocity to perform the normalization. The normalized maximum instantaneous power can be written as where P is the power per unit length. The involved force is the hydrostatic pressure force, and the velocity is the corresponding shallow wave celerity √ gh. In Fig. 15(b), the normalized maximum power P * vs Froude number collapse for different mean still water levels confirms that the scaling has been properly defined and it follows accurately a quadratic fit (obtained by least square fitting). This law can be used as an entry design tool to define the maximum power that is necessary to generate waves.

D. Active wave absorption
An active wall driven by a feedback controller may be useful not only to cancel wave reflections but also to attenuate wave impacts associated with extreme waves on a vertical wall and reduce their consequences. 45 We have here implemented an active wave absorption strategy using our results from the response of the wavemaker to velocity steps discussed in Sec. III B. First, consider a wave created at the wavemaker, propagating from left to right, whose shape is a leading trough, as shown in Fig. 16. To absorb this wave at the right active wall, a wave crest of nearly opposite phase has to be superimposed, which is generated by the motion of a wave-absorber, a ARTICLE scitation.org/journal/adv FIG. 16. Schematics of the wave absorption problem. The numerical domain is composed of (1) the piston wavemaker, (2) the piston wave-absorber on which the water level η w is measured, (3) the atmosphere, and (4) the seabed. The generated waves (5) at the free-water surface are referenced to the mean still water level (6) h. A second frame of reference, on the wave-absorber, is indicated as (x ′ , z ′ ).
wavemaker situated at the right side of the tank, according to the following strategy. Consider a wave-absorber consisting of an active flat wall with a sensor that measures the water level at the wall η w . On the frame of reference (x ′ , z ′ ) associated with this wave-absorber, the positive motion is from right to left (note that x ′ -direction is the opposite of x). The wall water level η w can be compared to a reference value η ref = 0 in order to attenuate reflections. The error ε = η ref − η w is permanently computed and fed into a proportional controller of gain K, which provides the absorption velocity UA = Kε. The corresponding block diagram of the feedback control strategy using a proportional controller is shown in Fig. 17(a). The efficiency of the control strategy relies on the choice of K where we propose to use the kinematic results of the step response rather than typical methods. 28 In Sec. III B, we have shown that the wall water level at the wavemaker reaches an overshoot value after a short time (as seen in Fig. 9) corresponding to the standard peak time of the system (see Fig. 11). The overshoot η 0 , the maximum water level at the wavemaker during the step, was related to the piston velocity UG by a linear relationship at the lower part (Fr < 0.2) of Fig. 10(b). If we want to absorb a wave of given amplitude at the wave-absorber wall, then we are more efficient if the change in amplitude, η w → η 0 (UG), takes into account the corresponding step velocity UG.
The slope of the linear scaling in Fig. 10(b) (the lower part of the plot) will determine the proportional controller gain K. The absorption velocity is then computed as UA = Kε = −Kη w , where K = 1/1.267 √ g/h represents the inverse of the linear fit slope in Fig. 10(b).
The dimensionless controller gain K * = K √ h/g is reproduced in Fig. 17(b) at different water levels, for a better understanding of the gain selection.
In order to determine the efficiency and limits of the feedback strategy, an irregular wave train, an undular bore, and regular wave cases were tested. The first example of the implementation of such an absorption strategy for irregular waves is shown in Fig. 18. A wave train is generated by the wavemaker on the left with the help of a smooth velocity pulse function defined as with S = 0.077 m, t 0 = 1.30 s, and τ = 0.342 s. This function can be visualized in Fig. 3(a). The wave train is moving in the direction of the wave-absorber at the right end part of the wave tank and is formed of a leading trough followed by wiggles (see time t = 3.5 s). The waves are then absorbed according to the feedback control strategy with an update of the absorption velocity at every time step, where the error at the wave-absorber wavemaker is plotted in  Fig. 19(a). The error alternates between positive and negative values causing the absorption velocity UA, which is also shown in the same graph, to behave similarly. The error decreases to zero while the wave-absorber performs the canceling action. The maximum error corresponds to 21% of the mean water level, that is to say, at the limit of the non-linear portion of the plot in Fig. 10(b). It results in an almost fully absorbed wave state at times t = 7 s and 8 s, where the water surface is calm at every location of the wave tank. Note that a reflected wave will take a time greater than t > 8 s to arrive into the wave-absorber after reflection at the left wavemaker. That means the error ε(t) falls to zero rapidly. In Fig. 18, in the last snapshot at t = 8 s, the black dashed curve is obtained with the controller off, and thus, the wave-absorber is at rest, allowing one to qualitatively compare the efficiency of the absorption strategy. It is possible to compute the energy of the wave tank in order to evaluate (and propose) the absorption efficiency of our method. The kinetic energy and potential energy per unit width (in J/m) for the water phase are defined as where ux and uz are the horizontal and vertical fluid velocities, respectively. The initial energy (at t = 0 s) is equal to the potential energy of the still water level, that is to say, E 0 = ρgLh 2 /2. As the tank length is not constant in time since both generating and absorbing wavemakers move at positions XG(t) and XA(t), a reference energy, which corresponds to the potential energy of the tank for a still water level retrieved by volume conservation considering these new positions, is defined as ).
Results of the energy computations are shown in Fig. 19(b), and they are shown with the initial energy state E 0 as a basis. We observe the decrease in the energy in the system during the generation process (t < 2 s), then a plateau corresponding to the wave propagation stage (2 s ≤ t ≤ 3 s), and finally an increase in the energy (until t = 8 s), which is due to the wave absorption of the incident waves and its convergence until a final energy value. We can compare these results with the case without absorption where we observe a nearly constant value of the total energy once the generating wavemaker stops. It is interesting to note the permanent trade between kinetic and potential energies as the waves reflect on the still wavemakers. The slight decrease in total energy is due to wave attenuation during propagation and reflection against the still walls. An estimation of the reflection coefficient can be made by computing the energy ratio between incident and reflected waves (adapted from Ref. 46), where EI and ER are the incident and reflected wave energies and t f is the time when the error reaches and stays inside 5% of the error band such that |ε(t ≥ t f )| ≤ max(|ε(t)|) × 5/100 is verified. The time t f as well as the 5% error band is indicated in Fig. 19(a). We make sure that at this time, no re-reflection on the absorbing wavemaker has occurred. Perfect wave absorption would lead to Etot = E ref (t f ), but as some reflection happens, this energy level is not reached. It is important to take into account the reference energy as the length of the wave tank is not constant, thus impacting the general level of potential energy. The computation of the reflection coefficient with this method leads to CR = 16%. Another analysis is carried out thanks to the separation of incident and reflected wave fields by means of the Fourier transform of two wave gauge data recordings at different tank locations. 46 It leads to the value of CR = 15%, which is similar to the previous one and shows the attenuation in the wave absorption process. The waterfall plot in Fig. 20(a) helps to visualize the motion of the wavemakers as illustrated by the extension of the wave tank at the left during the wave generation and the oscillations of the waveabsorber wavemaker as the incident waves make contact with the wall sensor. Finally, the efficiency of the process can be observed as

ARTICLE
scitation.org/journal/adv the reflected waves are small compared to the initial waves. A comparison with the case without absorption is also plotted in Fig. 20(b) and shows that the wave height at the (still, controller off) waveabsorber wavemaker fully fluctuates between low (η w < −0.025 m) and high (η w > 0.03 m) values. When the controller is turned on, those fluctuations disappear and η w → 0 as a result of the absorption mechanism.
The second test case consists in an undular bore, which is generated with a velocity square function defined as where U 0 = 0.2 m/s, t 0 = 2.5 s, and Θ is the Heaviside function, for a 4 m long wave tank and a still water level of h = 0.075 m. The undular bore, shown in Fig. 21, is generated on the left, propagates toward the wave-absorber wavemaker on the right, and is absorbed according to the proportional strategy. In this case, the error is almost always negative at the wave-absorber wall, with variations due to the incoming bore wiggles, as shown in Fig. 22(a). The starting error is important due to the high amplitude incoming undular bore, but rapidly the controller action decreases it and makes the stationary error converge to zero. In Fig. 22(b), the energy analysis is shown and a similar behavior is observed as in the previous case, with an increase in the total energy during the wave generation and a constant level during the propagation stage (0 s ≤ t ≤ 2.5 s for the first one and 2.5 s ≤ t ≤ 4 s for the second one). The total energy then decreases as long as the wave-absorber actuates and converges to the final energy, which corresponds to the final still position of both wavemakers. An estimation of the reflection coefficient can be made with the energy analysis and conducts to CR = 10%. The analysis of the reflection coefficient in Ref. 46 leads to CR = 16%. The difference can be explained by the method in Ref. 46, which is not well suited for high amplitude non-linear waves as can be the undular bore. This analysis shows that the strategy is interesting and can effectively absorb non-regular waves, even steep waves.
Finally, tests are carried out on harmonic cases, consisting of the absorption of regular waves generated as for the wave propagation mesh study in Sec. II D 2. The excitation functions and general setups are reported in Table VI. The reflection coefficient is calculated according to Ref. 46 and leads to a value lower than 10% for both cases, showing the efficiency of the proposed strategy in order to absorb waves. A summary of all test cases is reported in Table VI, and the reflection coefficients are given.
When we considered the steady state value of the step response, η ss , it was shown that the water level at the wavemaker reached a constant value after a short time, ts, at which η w /h = Fr, as shown in Fig. 9(b). The absorption velocity and constant coefficient K may be   . . . δ = −π/2, L = 8 m alternatively defined according to UA = Kε(t) with K = √ g/h. This coefficient was already given in Ref. 47 and was deduced from a mathematical study of the problem. The choice of the first strategy based on the overshoot (K = 1/1.267 √ g/h) rather than this last one is justified by the peak time tp associated with the overshoot, which is shorter than the settling time ts associated with the steady state, as shown in Fig. 11, or better said, the wavemaker moving at constant speed generates a transient wave that is used to cancel the incident wave.
The active wall driven by a feedback controller has proven to be useful not only to cancel wave reflections but also to attenuate high amplitude irregular wave impacts as in the undular bore example. In the future work, we will push further the absorption strategy to effectively reduce the consequences of extreme high amplitude waves using controllers for non-linear waves.

IV. CONCLUSION
In this work, we performed numerical simulations of a twodimensional wave tank in order to study the piston-type wavemaker initial-value problem and wave generation using the free and opensource code OpenFOAM. The numerical model reproduced the motion of a solid body piston-type wavemaker by moving a solid boundary driven by an external arbitrary signal waveform. We considered a fully viscous model solving the unsteady Navier-Stokes equations on the basis of a two-phase flow strategy and the volume of fluid method to capture the free surface dynamics. Velocity step signals (Heaviside functions) were applied to the piston-type wavemaker, generating a pulse-like wave that propagated along the tank followed by smaller waves or wiggles, which was identified as an undular bore. Recording of the wave elevation time series at the moving wall and in different tank locations was compared with theoretical data, providing a very good agreement and proving the capabilities of the OpenFOAM solver interDyMFoam to simulate twophase flows with wave propagation involving both free surfaces and moving boundaries. Wave elevation at the piston wall was found to have close similarity to the time response of the second order system found in feedback theory. In particular, the overshoot and rise, peak, and settling timescales were very close to those in the theory. The scaling found for water elevation at the piston wall at different step ARTICLE scitation.org/journal/adv velocities and mean still water levels was in close agreement with that in the theory for low Froude numbers. 11 At higher Froude numbers, the scaling differs considerably from that in the theory, being unable to take into account the main wave dynamics. The resulting main wave pulse is generated and detaches from the piston wall at the same time as the overshoot takes place in the wall elevation signal; thus, we call this wave the overshoot-wave. Results along the tank downstream agree with those of potential theory. The overshootwave propagates faster than piston velocity increasing its velocity and reaching asymptotically the shallow water celerity downstream the tank. As we solved fully viscous equations, we were able to quantitatively determine the power input during the step response associated with the wave generation process using the entire stress tensor at the piston wall. Net piston forces were obtained integrating pressure and shear stresses on the piston wall. A power scaling was found for different mean still water levels and step velocities as a function of the Froude number. Finally, in this work, we proposed a feedback proportional controller driving a secondary piston for wave absorption, where the controller gain was determined from the wavemaker step response. The feedback controlled piston method proved to be very efficient on both regular and irregular wave absorption. This novel approach provided the basis from which more complex active wave generation/absorption strategies can be further implemented on numerical and experimental wave tanks to improve efficiency under the influence of different parameters such as the water depth, the wave steepness, and negative velocity steps.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.