A numerical piston-type wave-maker toolbox for the open-source library OpenFOAM ＊

Dong-xu Wang, Jia-wen Sun, Jin-song Gui, Zhe Ma, De-zhi Ning, Ke-zhao Fang

1. Department of Ocean and Civil Engineering, Dalian Ocean University, Dalian 116023, China 2. National Marine Environment Monitoring Center, State Oceanic Administration, Dalian 116023, China 3. Deepwater Engineering Research Center, Dalian University of Technology, Dalian 116024, China 4. State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024,China

Abstract: A novel numerical piston-type wave-maker toolbox for the OpenFOAM is developed and demonstrated in this paper. This toolbox is implemented in C++ for improving the solutions of nonlinear wave problems in the field of hydrodynamics. As a toolbox that generates waves using the piston-type method only, it contains several frequently used functions, including the generation and the absorption of waves of various forms, an active absorption system and the porous media flow. Furthermore, to illustrate the operability of the toolbox, some validations and applications are presented, including the regular waves, the irregular waves, and the solitary waves.In each case, a satisfactory agreement is obtained in comparison with the published experimental or theoretical results,so this toolbox may be applied with a considerable confidence.

Key words: Piston-type wave-maker, OpenFOAM, nonlinear wave, hydrodynamics

Introduction

In the field of the ocean researches, both the experiments[1]and the numerical simulations[2]are widely used. A verified numerical model used to simulate waves of various forms may save a great amount of resources. The studies of the wave propagation[3]and the wave breaking[4]as well as the wave-structure interactions[5]are often carried out by numerical simulations. The OpenFOAM is an opensource library with many libraries and codes for solving complex CFD problems. Because of its high flexibility and zero cost, users can develop their own solvers based on idiographic problems. For these reasons, the OpenFOAM is one of the most widely used CFD toolboxes in the world and is used by practising engineers as well as researchers.

For ocean engineering, a two-phase incompressible fluid model named the InterFoam/InterDyMFoam is already available in the OpenFOAM, and most numerical wave flumes (NWFs) were developed based on it. With this model, coupled with the standard volume of fluid (VOF) scheme, the RANS equations are solved with the PIMPLE loop, a combination of the PISO and SIMPLE algorithms[6].Its main structure is inherited from the original PISO,and the outer loop of the PIMPLE algorithm is to handle the linearization of the convection term in the momentum equations. The detailed explanation of these algorithms in the application using the VOF method can be found in Jasak[7]. However, as far as we know, in this two-phase model, the generation and the absorption of different waves are not considered in the official version until now.

Therefore, the ability as well as the accuracy of generating and absorbing waves is the main concern for the further development. In most NWFs, the static boundary is selected for the wave generation, such as in Refs. [8-9]. This static boundary is a Dirichlet-type boundary in which the velocity and free surface elevations are updated at every time step according to the wave theory. The advantage of the static boundary wave generation is that it does not increase the computational cost of simulation. However, this method cannot handle problems with a strong nonlinearity, which means there is always an imbalance between the velocities below the wave crests and the wave troughs[10]. If no corrections are introduced, the mean water level in the wave flumes will increase as the simulation progresses. Consequently, to obtain a better result, a correction must be linked with the fixed boundary wave generation. To solve this problem, Jacobsen et al.[11]established a well-known NWF named the Waves2Foam using a method named the “wave relaxation zones”. The waves are generated by imposing the wave profile inside the relaxation zone, which could also absorb the reflection wave simultaneously. However, its computational domain has always to be increased at least by two wave lengths because of the relaxation zones at the inlet and outlet boundaries.

Unlike the static boundary wave generation, the moving boundary wave generation shares the same mechanism as the model test in the laboratory, which means that there is no flux across the boundary, and there will never be the problem of non-conservation of mass. Furthermore, it allows the complete simulations of laboratory experiments whether under shallow water or deep water. Because of these considerations,Higuera et al.[12]developed a three-dimensional pistontype numerical model with the volume-averaged Reynolds-averaged Navier-Stokes equations (VARANS)as the governing equations[13]. This model can not only simulate the interaction between the water and the porous coastal structure with excellent results but also absorb the secondary reflection at the wave paddle. However, when generating waves with its piston-type wave-maker, the movements of each paddle at every time step are needed[10], which may be slightly complicated if the simulation time is long.Additionally, the only active absorption method of their model is developed after Schäffer and Klopman[14], in which the shallow water regime is applied. Thus, the evanescent modes are insignificant while absorbing the secondary reflection in the vicinity of the paddle, but they may be significant under certain conditions. Other numerical models were developed[15-16], but their codes are not available.

Therefore, an open-source toolbox using the platform OpenFOAM-4.1 is desirable and is developed in this paper, specialized in the numerical piston-type wave-maker with various functions, which is easily applicable with a high accuracy. The piston-type numerical wave-maker is linked with an active absorption system (AAS), in which three methods to calculate the correct velocity are used under different conditions. Additionally, a limiter is added into the AAS to avoid the “slow drift”[14]of the wave paddle,as will be demonstrated in the next section. By accounting for the wet, dry and partial cells, the water level of each time step required for the AAS is calculated automatically rather than being input,which is one advantage of this toolbox. When generating and absorbing waves, only a few wave properties and alternatives are needed for input according to the various requirements, which means that this model is easier to use and extend. Another advantage is that this model could handle problems involving a strong nonlinearity, such as for the solitary wave whose nonlinearity (ε=H/d) is 0.5, which cannot be simulated using the static wave generation method precisely. Detailed discussions and comparisons of the abovementioned advantages are in Section 2.2.

1. Numerical model

1.1 Governing equations

1.1.1 Governing equations of the wave field

The governing equations of the wave field are the combination of the continuity equation and Reynoldsaveraged Navier-Stokes equations (RANS) as well as the VOF equation:

whereuis the velocity vector,αstands for the volume of the fluid in every cell, i.e., 0 for the air, 1 for the water, with values between 0, 1 indicating cells containing a mixture of the two fluids,*Sis the source term for the wave absorption in the damping zone and will be discussed in the next section, andρ,μare the weighted density and the dynamic viscosity of two phases, respectively, calculated by

Cκ∇αin Eq. (2) is the surface tension term, in whichCis the surface tension coefficient,κis the curvature of the interface.gis the acceleration of gravity,Xis the position vector andPρghis the pseudo-dynamic pressure

which does not have any real physical meaning but is used as a convenient numerical technique. It represents the dynamic pressure only if the free surface is located atZ=0. The third term on the left-hand side of Eq. (3) is the artificial compression term, which is only effective at the free surface, andUris the artificial compression velocity

wherecais a factor with a value of 1 by default and its value can be specified by users to enhance or eliminate (ca=0) the compression of the interface.

1.1.2Governing equations of the porous mediaflow field

In the porous media flow field, the governing equations are the VARANS equations[17], which also include the conservation of mass, the conservation of momentum and the VOF function advection equation:

whereφis the porosity,νis the kinematic viscosity andIis the hydraulic gradient

whereA,Bandcare coefficients depending on the physical properties of the porous material and on the flow regime. In the present toolbox,A,Baccording to Van Gent[18]are expressed as:

whereDstands for the average diameter of the material. The two factors1α,1βare equal to 500,2[19], respectively, in the present model.KCis the Keulegan-Carpenter number, which is defined as follows

whereTis the wave period,Uois the maximum oscillatory velocity. The additional mass coefficientcis shown to be less significant thanAorBto the results of the model and a value ofc=0 is applied by default. Users can also set other values.

After some transformation with both sides of Eq.(9) being multiplied byρand then divided byφ,the following equation can be obtained

A new variable is introduced for the convenience of presentation

and Eqs. (8)-(10) can be rewritten as follows:These three equations are similar to Eqs. (1)-(3), with only minor changes required to adapt the OpenFOAM to solve the VARANS equations[17].

1.1.3 Turbulence

Turbulence is not easy to define. Thus, the present model does not develop any new turbulence models but, rather, supports a large number of the turbulence models (e.g.,k-ε, SSTk-ωand LES)that are already included in the OpenFOAM. Thek-ε, SSTk-ωmodels have been applied inside the porous media[17].

1.2 Wave generation

The mesh motion in the OpenFOAM consists of a type of dynamic mesh, the solver and the diffusivity model[20]. Because there are many modules in the official version, selecting the appropriate modules for the development is important.

Clearly, the first step of the development is to determine the type of dynamic mesh. There are two types of dynamic mesh in the OpenFOAM: the dynamicFvMesh, the topoChangerFvMesh. The difference between them is whether the topology changes during the simulation or not. For simulating a wave paddle, the motion of the boundary is small, which means that the topology does not need to be changed.Therefore, the dynamicFvMesh is used in the present model. Among the classes of the dynamicFvMesh,there is an original class called the dynamicMotion-SolverFvMesh in which the mesh motion is obtained by solving a mesh motion equation with the boundary motion treated as a boundary condition to determine the position of other mesh points. The moving points of the other mesh points require models and the corresponding mesh motion equations to be solved

wherekdis the diffusivity, and it is used to determine how the points should be moved after solving the cell motion equation in each time step. In the present model, a diffusivity model named the inverseDistance is selected. When using that model,one or more boundaries are specified and the diffusivity of the field is based on the inverse of the distance from that boundary. For example, in terms of a two-dimensional wave flume, if the inlet boundary is specified as the moving boundary, the inverse distance of the outlet boundary will be the smallest, which means that the movement of the mesh points at the outlet boundary will also be the smallest.Umis the deformation velocity of the mesh points. The displacement of the points during each time step Δt is determined as follows

In the above scenario, a new class inherited from the fixedValuePointPatchField＜vector＞ is developed and adopted to simulate the displacement of the moving boundary. This new class of moving boundary is named the pistonWaveGeneration (Fig. 1), and the displacement of the wave paddle at every time step can be treated as a vector (x,y,z), that is to say, this vector should be (x,0,0) if the wave propagates along thex-direction. Clearly, the main job of the class is calculating that vector under different target wave conditions, which will be discussed in subsequent sections.

Fig. 1 Inheritance diagram for pistonWaveGeneration

1.2.1 Regular wave

The generation of regular waves is always a basic function of a wave-maker and plays a pivotal role in generating other waves. In the linear theory for a piston wave-maker, for regular waves, there is a linear relationship between the wave height and the stroke of the paddle

whereHsis the transfer function,His the wave height,Sdenotes the stroke of the paddle,kis the wave number anddis the initial water depth. If the free surface of the target wave is

the movement and the velocity of the wave paddle(x=0) can be expressed as

wheretis the simulated time,ωstands for the angular frequency.

1.2.2 Irregular wave

The actual waves in the ocean are irregular waves. Irregular waves are described by superimposing a series of regular wave components

whereiA,ki,iωandiθdonate the amplitude, the wave number, the angular frequency and the random initial phase of the componenti, respectively. The initial phase is distributed randomly in the interval[0,2π], therefore, this method is also called the Random Phase Spectrum Method (RPHM). The Jonswap, one of the most widely used wave spectra, is implemented in the present model for the wave generation:

wherefis the frequency,fpis the peak frequency, which is the inverse of the peak periodTp.γis the spectrum peak lifting factor, whose value is equal to 3.3 by default in the present model and could be modified by the users. The number of regular componentsNis 100 in the present model; thus, a discretizing function is needed to divide the total energy into 100 components. As shown in Fig. 2, the representative frequency of theith componentiωis the average of,

The shaded area stands for the energy of theith component. Since the wave spectrum is divided by the rule of the equivalent energy, the amplitude of every component is the same and calculated by

Thus, the velocity of the piston wave paddle is the linear superposition of every component

Fig. 2 The sketch of variables in the method RPHM

1.2.3 Solitary wave

In the wave theory, a solitary wave is the unidirectional permanent wave solution of the Kdv equation. The period of this wave is considered to be infinite. All of the free surface of the solitary wave is above the still water surface and can be expressed as

The integral of Eq. (36) is the displacement function of the wave paddle

Because of the infinite period, the stroke of the wave paddle is

However, the motion time for the paddle cannot actually be infinite, thus, a truncation of time is applied

Through the implementation of Eq. (39), the motion time of the wave paddle for the solitary wave generation is

So far, the displacement of the paddle is clear enough, but in the present model, a minor adjustment is made. The initial position of the paddle is set to(-T/2,-S/2), which means that the wave paddle will perform a unidirectional motion within the motion timeTto generate a solitary wave and stop if the simulation time is larger thanT. The final function is

1.3 Wave absorption

1.3.1 Absorption for reflection

The wave absorption, like the wave generation,plays a crucial role in the development of NWFs. In the present model, the reflection is absorbed by a damping zone, as shown by a source termS*in Eq.(2). The general definition ofS*is as follows whereχis a function for the damping layer. Linear types of damping layers are used in the present model:1

χis a parameter to adjust the wave damping effect,as was studied in detail by Perić[21],x0is the starting point of the damping zone,Lsdenotes the length of the damping zone. In the subsequent cases in this paper,1χis equal to 5.

1.3.2 Active absorption for secondary reflection

Due to the presence of the wave-maker, the secondary reflection, which always occurs at the wave paddle, becomes a major problem in both the physical model test and the NWF. Using a large wave flume is a widely adopted traditional way to reduce the influence of the secondary reflection. However, in a numerical model, this traditional approach may require a large amount of calculations. Therefore, the capability of the active absorption is often required to provide a better result with a lower cost. In the present model, the moving boundary for the wave generation is controlled to absorb waves impinging on it. Ifucorrdenotes the correct velocity of the active absorbing wave paddle, the corresponding displacement in every time step Δtis

The sum of these displacements constitutes the final correction of the active absorbing wave-maker.Thus, the final displacement of the active absorbing wave paddle is the superposition of the normal displacement and the correction part

For regular waves, three types ofucorrare applied in the present model and users can select any one of them to obtain a good result. The first, by Schäffer and Klopman[14], is the simplest method

heregis the gravitational acceleration, with the value of 9.81 m/s2.dis the water depth.mη,0ηare the measured, theoretical surface elevations at the wave paddle, respectively. In the second method, as the most widely used one,ucorrcan be expressed as

hereω, as noted above, is the angular frequency of the target wave,Hsis the linear relationship between the wave height and the paddle stroke. Other variables are the same as in Eq. (46). The last method,ucorrtakes the following form

In Eq. (48),xdenotes the displacement of the wave paddle; thus, the iteration in every time step is needed to solve this equation. In contrast to the linear transfer function for the progressive wave given in Eq. (22),Hsnis the transfer function of the evanescent wave and is calculated as

whereknis the wave number of the evanescent modes, due to the mismatch between the shape of the progressive wave velocity profile and the type of the wave paddle motion,knis real,kn＞0. This wave number is calculated by the dispersion relation for the evanescent modes:

Clearly, there are an infinite number of solutions ofknto this dispersion relation. Returning to Eq. (48), a definite value ofnis needed to make this numerical model effective. Zhang[22]derived a relationship between the transfer function andn(Fig. 3). In the figure, the abscissa axis is the non-dimensional frequencyand the ordinate axis is the value of the transfer function.is represented byR, and the number at the right side of the black curves is the value ofn. In the figure, the evanescent transfer functionRincreases continuously with the increasing frequency and the value ofRatn=20 is almost the same as that atn=200. Thus, in the present model,nis set to 20 by default. It is worth pointing out that if the non-dimensional frequency is larger than 3.4, the sum of the evanescent wave transfer function is larger than that of the linear one,which means that the local disturbance at the wave paddle exceeds the progressive wave.During the calculation, the third term on the right-hand side of Eq. (48) always makes the wave paddle have a “slow drift”[14], due to a linear growth of deviations in every period of the wave paddle, and it will cause a saturation, i.e., the paddle will end up with an extreme position outside of the stroke.Wang[23]analysed this problem, and proposed two modifications, i.e., the linear compensation and the PID control, to minimize the unnecessary motion of the paddle. The results seem very good, but the universal values of the parameters in the two modifications are not obtained, which means that if the wave condition changes, the parameters may need to be redefined. Therefore, in the present model, for convenience as well as accuracy, a limiter is added to this model to make sure that the displacement of the paddle is within a reasonable range. The corrected velocity is changed to

Fig. 3 The transfer function Hs and the total evanescentmode transfer function R with different n for a piston wave paddle[22]

Obviously, after adding this limiter, the method becomes a combination of the second and the third methods.Sis the stroke of the paddle, andμis the limiting factor, which is equal to 0.7 by default.The validation and the discussion of this method will be presented in the next section.

Because the irregular wave is considered as the sum of several regular waves, the active absorption is similar to regular situations. Notably, with the evane-scent model under the condition of irregular waves the calculation amount will increase because of 100 components. Thus, according to Eq. (47), the final velocityuof the wave paddle is directly given as follows

Fig. 4 Schematic diagram of the NWF of regular waves (wave gauge (WG))

hereiω,Hsiandη0iare the angular frequency,the linear transfer function and the theoretical surface elevation at the paddle of the componenti,respectively,mηis the measured surface elevation at the wave paddle,iλis a weight coefficient[24]calculated by

In the present model, as noted above,Nis equal to 100 and every component has the same energy,therefore,iλis equal to 0.01.

2. Validation

2.1 Regular/irregular waves

2.1.1 Wave generation

The regular wave generation is discussed in this subsection. A two-dimensional NWF is set up to simulate the generation of regular waves. The laminar flow model is used throughout Section 2. The initial water depthdis 0.50 m, and the wave heightHis 0.05 m. The period of the target wave is 1.6 s. The wavelengthL, according to the linear wave theory, is approximately equal to 3 m. The length of the NWF is 18 m, which is six times of the wavelength. The start point of the damping zone is set tox=12 m, so that the length of the damping layer is twice of the wavelength. The factor of the damping zone (χ1) is set to 5. Five wave gauges, whose positions in thexdirection are 1L, 2L, 3L, 4Land 6L, are labelled from #1 to #5, respectively. A schematic diagram of this NWF is shown in Fig. 4.

Fig. 5 Non-dimensional time histories of surface elevations at wave gauge #1-#4 of mesh I, mesh II and mesh III

First, the sensitivity of the model to the refinement and the density of the mesh is assessed using three different meshes, as labelled I, II and III, with the same background mesh and refined in the vicinity of the free surface fromx=0 m tox=12 m. Cells in the damping zone are not refined to enhance the effect of the numerical viscosity. Basic parameters of the three meshes are listed in Table 1. The total computation time is 20T, which is equal to 32 s. The initial time step is set to 0.001 s and is modified automatically according to the maximum Courant number, which is set to 0.5 in the present model. The non-dimensional time histories of the surface elevations at the wave gauges #1-#4 are shown in Fig.5. In general, the three curves are almost coincident,which means that the convergent results are obtained by these meshes. However, the results of mesh I are different from the others in detail at some moments(e.g.,t=16s, WG #3). Compared to mesh III, mesh II involves less cells to obtain the same result.Because of the abovementioned reasons, mesh II is selected to simulate the following examples of regular waves in the latter part of this subsection.

After validating the generation and the absorption of regular waves, the validation of irregular waves is discussed in this section. The wave properties are the same as the regular waves. The initial water depth is set to 0.50 m. The significant wave height is 0.05 m, and the peak period is 1.6 s,which means that the peak frequency is 0.625 Hz. The wavelength of the significant wave is still approximately equal to 3 m. Therefore, the NWF setup here is the same as the one used for regular waves (Fig. 4).However, because there must be many components for a smaller wave height, mesh III in Table 1 is selected for the simulation of the irregular waves,which means that Δznear the free surface isH/20.The linear form with a damping factor of 5 is used to absorb the target wave. The total computation time is 240 s, Δt=0.001s.

Table 1 Mesh parameters of the three meshes

Table 2 The statistical information of the irregular wave at WGs #1-#4

The time histories of the surface elevations measured by the WGs #1-#4 and their corresponding spectrums are shown in Figs. 6, 7, respectively. The spectrums at different positions in Fig. 12 almost coincide with each other as well as the target spectrum,which is a desirable result. The statistical information of the irregular wave is listed in Table 2, where it is shown that the significant wave heights of the WGs#1-#4 are 0.052 m, 0.053 m, 0.054 m and 0.052 m,respectively, which are slightly larger than the target height (H=0.005 m). The corresponding peak periods are 1.549 s, 1.518 s, 1.555 s and 1.533 s,which are slightly smaller than the target period(T=1.60 s). This may be due to the errors during the discrete process (Eq. (31)).

2.1.2 Active absorption system

The AAS of regular waves is tested in this section. Most arrangements are the same as those for the regular waves(2.1.1), but the damping zone is removed so there is only a vertical wall at the right side of the NWF, which means that there will be a total reflection if the AAS is off. The length of the NWF becomes 12 m accordingly. Three distinct types of corrected velocities, by the methods #1-#3 and defined by Eqs. (46)-(48), are applied and calculated.The displacements of the wave paddles are shown in Fig. 8. From the figure, we find that the displacement of the total reflection (with AAS off) is a sinusoidal motion, which can be confirmed by Eq. (24). With the methods #1, #2, very similar displacements are obtained because the relative water depthkdis small.Notably, with the method #3, a different displacement is obtained, and if it is not modified by the limiter, it will continue to decrease fromt=10s, which means that the system will be out of operation.

Fig. 8 Displacements of wave paddles in method #1, method #2 and method #3 as well as total reflection

The non-dimensional surface elevations obtained with these three methods, as well as the result in the case of the total reflection (with AAS off) at the wave paddle, are compared in Fig. 9. Thex-axis and they-axis stand for the time and the non-dimensional surface elevation, respectively. Results from the WG#1, which are similar to those from the WG #2, are not shown in Figs. 9-11 because of the limitation of space.The elevations calculated by different AASs see almost no difference. The final stable amplitudes at the WGs #2-#4 are twice of the target wave height.Thus, the stable standing waves can be observed in comparison with an increasing wave obtained by the total reflection. This means that the secondary reflection at the wave paddle is absorbed by the AAS.

Fig. 9 Non-dimensional time histories of surface elevations at WGs #2-#4 of method #1, method #2, method #3 and total reflection

For a quantitative analysis of this stability, a variable called the “stability factor” is considered. In the time range [t-Δt,t+Δt], a stability factor is defined as

hereηmax,ηminare the maximum, minimum surface elevations in the time range, respectively,Hmis the average wave height of this wave profile,r(t)should be approximately equal to 1 and does not change with time sharply if the wave field in the NWF is stable. The time range is [10T,18T] with Δt=2T,which means that the data in the time range [8T,20T]are used. The stability factors in those four cases at the WGs #2-#4 are shown in Fig. 10. The stability factors obtained by the three methods take approximately the value of 1.0 or approaching this value. The stability factor of the total reflection is always increasing and will continue to increase as the time progresses. This condition indicates that the stable wave fields are in those NWFs when the AAS is on. In other words, the AAS in the present model is effective in eliminating the secondary reflection.

Fig. 10 Stability factors at WGs #2-#4 of method #1, method#2, method #3 and total reflection

The AAS validation for irregular waves is performed as follows. With the same arrangements as in the case of the irregular waves(1.1.1), but the damping zone is also eliminated so there is a vertical wall at the right side of the NWF and the length of the NWF is 12 m.Two cases, with and without the AAS, are calculated.In contrast to the regular waves, there is no stable standing wave in the irregular NWF. Therefore, it is essential to separate the incident and reflected waves.A time-domain two-point method is employed to separate the incident and reflected waves. Because of the method, there is one more wave gauge located atx=2 m for the separation. Figure 11 shows the time histories of the surface elevations when the AAS is on/off. Figure 12 is the separated spectrum calculated by the present model. From those two figures, we can see that lower surface elevations are present with the AAS. Moreover, the separated spectrum agrees with the target spectrum, which means that this AAS has a satisfactory effect.

Fig. 11 The time histories of surface elevations with the AAS on /off

2.2 Solitary wave

The validation of the solitary wave is performed in this subsection. As mentioned before, the pistontype wave-maker is more suitable to handle a wave with a strong nonlinearity (ε＞0.3) than the other numerical methods. To show this, the production and the propagation of a solitary wave with a strong nonlinearity (ε=0.5) are not only simulated but also compared with other popular solvers (e.g., the Waves2Foam and the OLAFoam). With the OLAFoam, the solitary wave is generated using the static boundary method. It is a Dirichlet-type boundary in which the velocity and the free surface elevations are updated at each time step following the wave theory. The wave theory for the solitary waves involved in the OLAFoam is the Boussinesq theory and the Grimshaw theory. Unlike the OLAFoam, with the Waves2Foam, the solitary waves are not generated at the boundary. A solitary wave, based on the Boussinesq theory and calculated by the wave properties we set, is generated in the NWF in the first time step. In addition, the velocities at different positions are calculated and passed to the field to make the solitary wave propagate. An NWF is set up for the solitary wave simulations using the present model as well as other solvers. The length, the height of the NWF is 30 m, 2 m, respectively. Δx,Δzare equal to 0.04 m, 0.02 m, respectively, which indicates a total of 75 000 cells. The water depth is 1.0 m and the amplitude of the solitary wave is 0.5 m, thus, the n onlinearity (ε=H/d) is 0.5. The total computation time is 10.0 s and the time step is 0.001 s. Figure 13 shows the stable non-dimensional surface elevations obtained by different solvers as well as the theoretical value att=3.45s . Two theories used in the OLAFoam show diverse results although with the same wave properties. With the Boussinesq theory, a large wave is generated with a non-dimensional wave crest of approximately 1.35. A smaller wave, with a non-dimensional wave crest of 0.85, is obtained by the Grimshaw theory. In comparison with the OLAFoam,the non-dimensional wave crest obtained by the Waves2Foam is approximately 0.9, despite the good agreement in the first time step. The results of the present model agree well with the theoretical value not only in this time step but also in other time steps,which means that this piston-type wave-maker is more suitable to handle a wave with a strong nonlinearity.

Fig. 12 The separated spectrum calculated by the present model:Target spectrum, spectrums before and after separation

2.3 Porous dam break

A succinct validation of the porous media flow is presented in this subsection. The porous dam break experiments by Liu et al.[25]are reproduced. The twodimensional NWF is 0.89 m long, 0.40 m high. A porous dam of 0.29 m long, 0.40 m high is placed at the centre of the NWF,x=0.30 m tox=0.59 m.The mean size and the porosity of the porous material is 0.0159 m, 0.49, respectively. A reservoir is set at the left side of the NWF with a gap (0.02 m) between the reservoir and the porous dam. Two initial

Fig. 13 The non-d imen sional surface elevatio ns from the present modelandotherfamoussolversaswellasthe theory value at t=3.45s . The present model, theoretical solution, Boussinesq model in OLAFoam, Grimshaw model in OLAFoam, Waves2Foam

Fig. 14 (Color online) The numerical result comparisons of the high water level and the corresponding experimental data at different time steps

water levels labelled as high (d=0.355 m)and low(d=0.136 m)are considered. There is water at a height of 0.02 m at the bottom of the NWF in the first time step. Both Δx, Δzare equal to 0.005 m, which indicates a total of 14 240 cells. The total computation time in the two cases is 2.5 s, and the time step is 0.001 s. The comparisons of the numerical results for the High water level and the corresponding laboratory data in different time steps are shown in Fig. 14. The water rushes to the dam at the beginning of simulation.Then, a strong impact occurs and the water is accumulated to form a small upward jet on the surface of the porous dam (T=0.35s). As the time progresses, the water flows out from the dam and, with the combined action of the reflection from the right wall of the NWF, to form a vortex (T=1.55s).These phenomena are not only observed in the laboratory but also in the numerical simulation. The comparisons of the numerical results of the Low water level and the corresponding laboratory data in different time steps are shown in Fig. 15. Most of the phenomena are the same as the High water level, but there are no vortexes observed. The overall agreement observed in Figs. 14, 15 is very good, which means that the porous media flow in the present model is effective.

3. Conclusions

An open-source toolbo x of a numerical pistontype wave-maker is utilized in this paper. The boundary for the wave generation, as well as the AAS,is developed based on the class inheritance. The wave absorption is achieved by adding a damping zone at the end of the wave flume as well as an AAS linked with the moving-wall wave-generation boundary. This toolbox is tested and discussed in Sections 2, with typical examples. Conclusions of the paper are as follows:

This toolbox has a satisfactory capability to simulate several common types of waves in the ocean researches, such as the regular waves and the irregular waves. Additionally, compared to other available solvers, the present model is more suitable for the simulation of strongly nonlinear waves. Therefore, the present model can be used to simulate completely the actual piston wave-maker in the laboratory so that the nonlinear hydrodynamic problems in the ocean engineering can be handled expediently.

The AAS is used to eliminate the secondary reflection, which always occurs at the wave paddle.The “stability factor” method is used in Section 2.1.2.The results show that the AAS has a good effect under the conditions of regular waves. For irregular waves, a time-domain method is used to separate the incident and reflected waves. The spectrum in the present model agrees with the target spectrum.

The above results show the model's good capability for handling strongly nonlinear problems, which is important during the development of numerical wave flumes. In general, the above results indicate that the present model is well suited for handling many ocean engineering problems under a wide range of nonlinear wave conditions. Future work may focus on the wave interaction with floating structures under various wave conditions.

Fig. 15 (Color online) The numerical result comparisons of the low water level and the corresponding experimental data at different time steps

Acknowledgement

The work was supported by the Dalian Ocean University and the Deepwater Engineering Research Center of Dalian University of Technology as well as the National Marine Environment Monitoring Center.