- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
Advances in Mechanical Engineering
Volume 2013 (2013), Article ID 308436, 17 pages
Finite Element Based Viscous Numerical Wave Flume
1Research and Development Center of the Civil Engineering Technology, Dalian University, Dalian 116622, China
2School of Hydraulic Engineering, Dalian University of Technology, Dalian 116024, China
3Centre for Deepwater Engineering, Dalian University of Technology, Dalian 116024, China
Received 2 September 2013; Accepted 12 October 2013
Academic Editor: Nao-Aki Noda
Copyright © 2013 Jianmin Qin et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
A two-dimensional numerical wave flume (NWF) for viscous fluid flows with free surface is developed in this work. It is based on the upwind finite element solutions of Navier-Stokes equations, CLEAR-volume of fluid method for free surface capture, internal wave maker for wave generation, and sponge layer for wave absorbing. The wave generation and absorption by prescribing velocity boundary conditions along inlet and radiation boundary condition along outlet are also incorporated. The numerical model is validated against several benchmarks, including dam-breaking flow, liquid sloshing in baffled tank, linear water wave propagation and reflection from vertical wall, nonlinear solitary wave fission over sharp step, and wave-induced fluid resonance in narrow gap confined by floating structures. The comparisons with available experimental data, numerical results, and theoretical solutions confirm that the present numerical wave flume has good performance in dealing with complex interface flows and water wave interaction with structures.
Water wave motion, as one kind of interfacial flows, is widely encountered in ocean engineering. The great challenges remained nowadays for the water wave problems include the numerical predictions of the complex wave propagation, transformation, and wave-structure interactions, which usually involve difficulties such as the violent free surface, turbulence flow, viscous effects, and complexity of boundary configurations. Under these situations, the viscosity and vortices play rather important role, and hence the viscous fluid model based on the numerical solutions of the Navier-Stokes equations becomes necessary since the convectional potential flow model assumes that the fluid is inviscid and the flow is irrotational. Promoted by the great demands of engineering practices and attributed to the fast development of computational facilities and numerical techniques, computational fluid dynamics (CFD) becomes a powerful tool for the ocean hydrodynamics, which consequently leads to the growth of two-dimensional numerical wave flume (NWF) and three-dimensional numerical wave tank (NWT). Although the theoretical and experimental investigations play still the most important role in ocean engineering, the NWF/NWT technique has become more and more popular due to its inspiring advantages over the convectional laboratory tests and theoretical analysis. The limitations in the convectional approaches include that, for example, the theoretical predictions of hydrodynamic characteristics are usually limited to simple geometries and physics and the physical modeling is extremely difficult in matching the model scale, and it is time-consuming and costly.
Generally speaking, the existed numerical wave flume/tank can be divided into two main groups depending on the basic mathematical model, that is, the potential flow theory model and the viscous fluid model. The former is preferred for the situations that the fluid viscosity and vortical flow are negligible and works with high numerical efficiency. A great number of numerical models based on the potential flow theory have been proposed, including the linear and nonlinear models in both frequency and time domains. However, due to the primary assumptions made for the traditional potential flow theory, that is, the irrotational flows of inviscid fluids, together with the limitations of single-value height function used for the free surface predictions, the potential theory based NWF/NWT is not available for the complex water wave problems with strong nonlinearity of interface and significant mechanical energy dissipations. On the other hand, the NWF/NWT based on viscous fluid model performs much better with great generality and remarkable ability although the computational efforts are expected to be increased.
The viscous numerical wave flume is highly dependent on the numerical solution of the Navier-Stokes equations. As an important numerical method, the finite element method is widely used in solid mechanics and structural analysis in contrast to the popularities of the finite difference and finite volume methods in the community of computational fluid dynamics . However, the merits of finite element method in CFD should be addressed, for example, the flexible adaptation to complex boundaries in conjunction with irregular mesh partition, the high numerical accuracy, and the convenience in modularization programming . The most important potential benefit of using the finite element method for the NWF/NWT consists in the unified simulations for the fluid and structure interactions (FSI), in which the dynamics responses of the structures should be considered as well. Under this situation, the finite element method can be used as the uniform numerical method for both fluid flow and structure analysis in the frame of arbitrary Lagrangian-Eulerian reference system .
In addition to the necessary flow solver, the free surface prediction is another critical issue for a numerical wave flume/tank. In general, the methods developed so far for predicting the flow interface fall into two categories, that is, the Lagrangian method and the Eulerian method. In Lagrangian approach, the discrete grid nodes can be regarded as material points to move with the fluid flow, giving rise to the moving grid method . In this sense, it is also called the interface tracking type method. For the simple cases, the interface tracking is a desirable choice attributed to its less computational demands and high numerical precision. However, if the extremely large deformation of free surface is involved, the tracking method becomes invalid, leading to severe distortion of the computational meshes and finally the failure of computations. In addition to the moving grid method, usually used with viscous fluid models, the height function method, for example, the wave profile function defined by the Taylor expansion with respect to the still water level, can be also grouped into this class, which generally works in conjunction with the potential flow models. By using the height function method, it can bring great convenience since the computational meshes are allowed to remain unchanged and the free surface boundary conditions are directly incorporated in the definition of wave profile. However, the height function is a singlevalue function, bearing limitations for complex interfacial problems. The Eulerian method, on the other hand, is free from either the mesh distortion as the moving grid method or the limitation of singlevalue as the height function. Hence, it plays more important role in the numerical wave flume/tank in the context of viscous fluid flows. In Eulerian frame, the computational grids are fixed throughout the numerical calculation while the free surface is captured by means of a marker function. The Eulearian interface capturing methods include the Marker and Cell (MAC) , volume of fluid method (VOF) , and level set method (LS)  among others. The marker and cell method is the ancestor of volume of fluid method, which captures the free surface with the marker particles assigned in the computational cells. The main drawback of MAC is the huge requirement of computer capacity. Hirt and Nichols  raised the idea of VOF instead of MAC, in which the free surface is captured by a fractional volume function. The volume of fluid method is theoretically available to any complex free surface flows even involving the folding, splashing, and breaking processes. Several numerical wave flumes/tanks have been established based on the volume of fluid method, for example, the COMBRAS model by  and the numerical models developed, respectively, in references [9–13].
The original volume of fluid method  was actually proposed in the context of finite difference method. Therefore the orthogonal rectangular computational meshes are required, not only for the solution of governing equations but also for the implementation of the volume fluid advection, which brings difficulties in dealing with the problems of wave action on structures with rather complex boundary configurations. The combination of VOF with irregular mesh partition remains still great challenges, although some achievements have been made, for example, the coupling of VOF with finite volume method by using the triangular mesh cells , where the concept of flux across cell face is adopted, in addition to the other work by [15, 16]. As for the finite element method, Ashgriz et al.  proposed a novel and efficient volume of fluid method, referred to as the computational Lagrangian-Eulerian advection remap volume of fluid method (CLEAR-VOF), in which the volume of fluid can be advected in consistence with any irregular mesh partition and the interface reconstruction employs a piecewise linear interface construction (PLIC) strategy. In CLEAR-VOF, the fluid transport is conducted by means of the Lagrangian movement of a “fluid polygon”, and then the fractional volume of fluid is remapped on the Eulerian meshes, which is essentially distinct from the idea behind the concept of donoracceptor  or flux crossing cell side . Since the CLEAR-VOF method is applicable to the irregular computational meshes, the collaboration of VOF and FEM consequently becomes feasible. It should be addressed here that the level set method, as a new interface capture technique, was also introduced into the field of wave hydrodynamics. However, the level set method requires the numerical computations to be carried out for the two phases on the both sides of the interface. Generally speaking, the level set method is more time-consuming than the volume of fluid method, since the single-phase simulations are allowed for the VOF.
Analogue to the facilities of wave flume/tank in laboratory, the NWF/NWT should have the function to produce desirable wave trains. A straightforward way for the wave generation in a numerical model is to impose the velocity boundary conditions together with the wave profile along the inlet of the computational domain. The required velocity components and wave profile can be evaluated according to an appropriate wave theory. This method is the common approach for wave generation and proves to be efficient. However, under some complex conditions, for example, the water waves that interact strongly with structures and involve significant reflection waves, the method of wave generation by describing inlet boundary conditions may be problematic because the reflected wave will travel upstream and eventually collides with the inlet boundary. These interactions further result in the so-called severe second reflection, the same as the situation encountered in laboratory tests. The increase in the size of wave flume/tank seems to be a solution. However, it is practice impossible for the laboratory facility building and will give rise to the excessive increase in computational costs for numerical modeling. In order to overcome this difficulty, the technique of internal wake maker (IWM) was proposed. As for the numerical wave flume based on Navier-Stokes equations, Lin and Liu  devised a promising approach by means of adding a source term to the continuum equation. With internal wave maker, the target waves can be produced numerically. And the undesirable refection wave is possible to be reduced greatly and even eliminated completely, although the second reflection may still occur if structures are involved.
Similar to the wave generation, both numerical and physical wave flumes/tanks require the function of wave absorption. Inspired by the artificial beach used in the laboratory tests, the damping zone is most often incorporated in the NWF/NWT, which can be realized numerically through introducing a damping force into the Navier-Stokes equations or a damping term into the free surface boundary conditions for the potential flow models. The wave damper devised for numerical model is usually dependent on the velocity field and its robust performance in eliminating the undesirable reflection wave can be expected.
As far as the simple harmonic water waves are concerned, the application of Sommerfeld-Orlanski radiation boundary condition on the outlet may produce satisfying results without giving rise to wave reflection. Note that the usage of radiation boundary condition in the numerical model serves as an advantage over the physical model test since the radiation boundary does not exist in any laboratory facility. However, this method loses its attractiveness if the wave phase celerity becomes unclear, for example, the irregular waves or the regular waves disturbed significantly by structure or uneven bottom.
Although the wave damper is preferred by taking consideration its generality, it does not mean that the approach with damping zone can be applied to any cases. If the wave propagation is accompanied by the mass transportation, such as the typical high-order Stokes waves, the introducing of artificial damping zone may lead to sustaining accumulation of fluid bulk. Under this case, the coupled strategy of using together the damping zone and Sommerfeld-Orlanski boundary condition at outlet is suggested, where the local averaged horizontal velocity component can be adopted as the phase speed.
This work focuses on developing a comprehensive finite element based two-dimensional viscous numerical wave flume (2D NWF). The numerical implementations associated with the previously mentioned necessary four components, namely, the flow solver, interface predictor, wave maker, and wave absorber, will be described in detail. Several numerical validations against benchmarks will be presented in Section 3, and finally the conclusions are drawn in Section 4.
2.1. Governing Equations
The governing equations for the motion of homogeneous incompressible viscous Newtonian fluids can be described by the following continuity equation and Navier-Stokes equations, being of the form where is the velocity component in the th direction, is the pressure, is the fluid density, is the fluid viscosity, is the external force and , a source term used as one of the options to generate the desirable target waves in this work. The source term is activated only in the source region under necessary situation. Otherwise, it remains zero, either outside of the source region or with the other wave generation methods. Various formulations have been designed for the source term in order to generate desirable target waves . As for the linear Airy wave, the source function reads where is the wave height, is the area of source region, is and the wave frequency with being the wave period.
In addition to the wave generation, the efficient method to reduce the reflection waves plays the same important role for a numerical wave flume. This is implemented in this work by introducing an artificial damping force into the momentum equation, serving as the artificial beach in the laboratory tests to dissipate the wave energy. The body-force term in (2) is accordingly assumed to be where is the gravitational acceleration and denotes the artificial resistance used to damp out the reflection waves. Similar to Kim et al. , the damping function used in this work reads where is the horizontal coordinate at the start point of sponge layer, is the total length of sponge layer, at least longer than two times of the incident wave length, is the elevation of flume bottom, is the prescribed maximum elevation of the free surface, and is the damping coefficient. For the sake of numerical stability, the artificial damping force works only in the vertical direction. The damping zones are generally arranged at both ends of the computational domain. A sketch definition of the numerical wave flume with internal wave maker for wave generation and sponge layer for wave absorption is shown in Figure 1.
2.2. Numerical Discretization and Implementation
The Navier-Stokes equations are typical advection-diffusion equation. It is well known that the standard Galerkin finite element method is analogue to the central finite difference scheme, which may lead to the numerical instability when it is used to solve the advection-dominated problems. In order to overcome this difficulty, the upwind FEM schemes have been developed, such as the Taylor-Gelarkin (TG) , Streamline Upwind Petrov Galerkin (SUPG) , Galerkin Least-Squares (GLS) , and Characteristic-Based Split (CBS) . The Taylor-Gelarkin scheme considers the upwind characteristics of fluid flows by means of representing the time derivatives of with spatial derivatives. A three-step Taylor-Gerlarkin finite element method was proposed [24, 25], which proves to be with third-order accuracy in time but without introducing the higher order derivatives as the normal Tayler-Galerkin scheme. Another advantage of the three-step finite element method is that it allows using relatively large time step and saves the computational efforts since it satisfies the uniform Courant-Friedrichs-Lewy (CFL) condition. Employing the three-step finite element method, the momentum equation can be discretized in time as follows: where the superscripts , , and represent the time instants at , , and , respectively, with being the time interval and .
It can be seen that the unknowns of velocity can be obtained by explicitly solving the above equations in sequence. However, according to (8), the pressure at the time level has to be solved beforehand in order to obtain the velocity . By applying the divergence operation on the both sides of (8) and incorporating the continuity equation (incompressibility) attime level, that is, or , we may have the pressure Poisson equation, that is, (9).
Since the fractional strategy is used in the three-step finite element method, the standard Galerkin weighted residual method can be adopted for the spatial discretization. Thus the corresponding finite element equations in elemental level can be obtained, referring to (10)~(11), where is the weight function and shares the same form with the trial function in Galerkin method, and denotes the component of outward normal unit vector in theth direction. The corresponding weak form of (9) in a finite element is formulated in (13):
Note that since the projection procedure is involved in three-step finite element method, that is, the pressure is decoupled from the momentum equation and is solved separately, the same order interpolation function (the bilinear trial function used in this work) can be adopted for the pressure as the velocity, which satisfies with the Ladyzenskaja-Babuska-Brezzi (LBB) condition.
As the elemental coefficient matrix is obtained according to (10)~(13), the global matrix should be assembled according to the accumulation rule, which finally leads to the algebraic system in terms of the unknowns of velocity and pressure. As for the momentum equations, due to the explicit scheme used in the time discretization, the coefficient matrix of velocity presents a diagonal dominated linear equation set. Therefore, the lumped mass matrix is used to solve the velocity, which proves to be in high computational efficiency but with little numerical errors. On the other hand, the solution of pressure Poisson equation is very timeconsuming. In this study, the preconditioned BiConjugate Gradient Stabilized (BICGSTAB) method is adopted . Although BICGSTAB solver is originally devised for asymmetric linear system, it is helpful to the convergent acceleration and numerical stability when it is applied to the sparse symmetric problems.
2.3. Free Surface Capture
The computational Lagrangian-Eulerian advection remap volume of fluid method (CLEAR-VOF) is adopted in this work for the free surface prediction, which is a novel interface capture approach with inherent consistence with the unstructured irregular mesh partition in the context of finite element method. The main ideas behind the CLEAR-VOF method are based on the concepts of Lagrangian transportation of “fluid polygon” for the fluid advection and piecewise linear interface construction (PLIC) scheme for the free-surface reconstruction.
As for the fluid advection step, the notion of “fluid polygon” is introduced in CLEAR-VOF. As is well known, the fractional function of fluid volume, denoted by herein, varies from 0.0 to 1.0, defined by the occupied fluid volume over the volume of computational cell. In CLEAR-VOF, corresponds to a “full cell,” for which the “fluid polygon” coincides with the mesh cell, while presents a void cell, indicating that there is no fluid in the computational cell, and hence the “fluid polygon” does not exist. As for the case of , a “partial cell” is defined in a mother cell (the computational mesh cell contains fluids), for which the “fluid polygon” should be reconstructed according to the local distribution of VOF function.
Assuming that the “fluid polygon” has been determined, the motion of the “fluid polygon” can then be described approximately in the Lagrangian sense by using the velocities at its vertices in a tiny time step. As for the case that the “fluid polygon” is just the full cell with , the velocity solutions obtained on the elemental nodes can be used directly for the fluid advection. As for the partial cell, the velocity components on the vertices of the “fluid polygon” have to be interpolated based on the nodal velocities of the mother cell. A detailed illustration of the advection of fluid volume in CLEAR-VOF is shown in Figure 2.
Figure 2 shows that, for example, a fluid polygon ABJK is defined in the mother cell ABCD. With the available velocity vector on its vertex , the displacement after one time step can be evaluated by which finally leads to the advection of “fluid polygon” from ABJK to abjk in a Lagrangian sense, where the superscripts and stand for the two sequent time levels. In order to determine the new fractional volume at thetime level, the fluid volume in abjk should be distributed to the background Eulerian meshes. This is achieved by a computational geometry procedure, that is, to calculate the intersections of the fluid polygon abjk with the possible computational cells. For example, the present fluid polygon abjk will insect with the adjacent cells ABCD, BEFC, and DCGH and gives the intersection areas of apqrs, pbjq, and srk, respectively.
On the other hand, as for a particular mesh cell, for example, considering again the typical mesh cell ABCD, part of the fluid volume will be remained in it, that is, the volume of apqrs. Meanwhile the adjacent fluid polygons may also bring some fluids into it. The summations of the fluid volume remained in and flowing into the mesh ABCD can then be used to update the new fractional function of the volume of fluid: whereis the area of the current computational mesh cell and means the intersection area with theth fluid polygon. After the fluid volume advection and accumulation are conducted throughout the total computational cells, the new distribution of can be obtained for the next time step and the interface reconstruction can be followed.
The significance of interface reconstruction is not only sharpening the configuration of free surface but also saving computer costs since the very fine meshes are not necessary in order to predict the free surface with the same precision. That is the meanings of using the VOF method to capture the interface. Moreover, it is also a requirement from the construction of fluid polygon as mentioned above. For the present numerical model, the popular PLIC method is adopted and a brief description is included as follows. Firstly, the local normal outward unit vector for an interface cell should be calculated by using the local gradient of , which is defined at the centre of the corresponding computational cell: where is the fractional function of volume of fluid and is the differential operator. The details on calculating the local gradient , requiring the appropriate arrangements of the imaginary elements, are referred to . With the obtained direction of free surface, a line segment passing the interfacial element can be defined as where is a constant unknown requiring to be determined by an iterative procedure. The convergent criterion for the iteration is devised as that the fluid volume in the current element beneath the interface line approaches to the determined fluid volume in it. In this work, the uniform criterion is used. The accomplishment of fluid advection and interface reconstruction means that the computational domain and the necessary boundary information for the next step numerical calculation are available.
It should be mentioned that the fractional volume of fluid obtained numerically as above might be with unphysical values attributed to the numerical errors. During the computations, the real-time numerical monitor is necessary. If an exception of occurs in an element, it will be reset to unit. For convenience, if the fraction volume is approaching to be zero or unit, the numerical truncation is applied as follows: where the fuzzy number is adopted in the present study. Moreover a partial element will be forced to be full element if it completely surrounded by full cells. And it will be reset to be void if it is not adjacent to at least one full element. By applying these implementations, the unphysical error can be eliminated completely. In addition to these local adjustments, a global correction concerning the mass conservation is also required. This is realized through increasing or decreasing the fractional volume for the free surface cells: where is the fractional volume in a free surface element, is the area of the free surface element, and is the total volume of imbalance in the computational domain. Note that a full element may be a free surface element and should be considered in the above formula, which is different from that used in , where the global adjustment is only conducted for partial elements. The present revision proves to be more robust according to the numerical calculations included in this paper. Note that the previous local adjustment must be carried out again after the global mass correction since the fractional volume may be changed during the global correction, especially for the possibility involving unphysical value of fractional volume. This procedure is actually implemented by a numerical iteration.
2.4. Initial and Boundary Conditions
In the previous sections, the numerical discretization of governing equation and the interface capture with CLEAR-VOF method are described. In order to obtain proper solutions for a particular physical problem, the governing equations must be subjected to appropriate initial and boundary conditions.
As for the water wave problems concerned in this work, it is convenient to start the simulations from the still water state. Therefore, the initial conditions can be applied as follows, that is, the zero velocity and the static water pressure, where represents the position vector and represents the vertical coordinate.
For a two-dimensional numerical wave flume, the boundaries are generally composed of the inlet, outlet, solid wall, and free surface. The inlet boundary often serves as the wave generator in addition to the internal wave maker. Generally speaking, if the problem is not involved in significant wave reflections, the method for wave generation with the inlet boundary conditions is preferred due to its simplicity in numerical implementation. In practice, the water wave solutions based on the potential flow theory in terms of the velocity components and wave profiles along the inlet are often adopted as the boundary conditions: where the time-dependent wave profile should be converted into the fractional volume fluid.
On the other hand, if the internal wave maker is utilized. The inlet boundary disappears, leaving one/two outlets for the two-dimensional numerical wave flume. The corresponding boundary conditions along the outlets depend on the types of wave absorption. For the simple case, it is convenient to adopt the Sommerfeld-Orlanski boundary condition , which reads where presents the physical variables of velocity and wave profile, presents the local outward normal vector at outlet, and the phase speed of water wave.
For the situations that the phase speed is not clear, as an approximation, the averaged velocity along the vertical outlet can be specified. The above radiation boundary condition is generally used together with the prescribed pressure, where and is the local water depth and time-dependent wave profile, respectively.
As a more general method, the artificial sponge layers can be introduced, leading to that the velocity components and wave profile along the outlets approach to be zero. Under this circumstance, the exterior ends of the sponge layers can be simply regarded as solid wall. For some cases that the travelling waves in the sponge layer cannot be damped out, the Sommerfeld-Orlanski boundary conditions are also incorporated along the exits and work together with the sponge layer. Such a coupling strategy enables the user to shorten the length of sponge layer and hence receive considerable improvement of the computational efficiency.
As far as the solid wall is concerned, including the flume bottom, structure surface, and the outer ends of the sponge layers, the on-slip boundary conditions are applied: where is the time-dependent velocity of the solid wall surface and the zero gradient of pressure is used for the pressure Poisson equation since it should be solved separately in the present numerical model.
The most important for a water wave problem is the treatment of free surface boundary conditions. According to the requirement of stress balance, both pressure and velocity boundary conditions should be imposed along the free surface. As for the present numerical wave flume, the surface tension is neglected, thus the interface dynamic boundary condition reads where and are the pressures in the water phase and air phase, respectively,is the dynamic viscosity of the air, and the subscripts and denote the normal and tangential directions, respectively. By neglecting further the viscous effect, a simplified normal dynamic free surface boundary condition can be used directly for the pressure Poisson equation. However, the velocity boundary conditions on the interface are not able to be incorporated in a straightforward manner. Alternatively, the velocity extrapolations [28, 29] are adopted in this work.
2.5. Numerical Procedure
The overall numerical procedure for the present numerical wave flume is summarized as follows:(1)generating initial computational mesh covering a fairly large computational domain: either regular mesh or irregular mesh is compatible with the present model;(2)specifying appropriate initial and boundary conditions;(3)flow simulations with three-step finite element method including the following:(i)employing (6) to obtain the velocity based on the known variables at the th time level;(ii)calculating with (7);(iii)solving pressure Poisson equation, that is, (9), to obtain at the next time step;(iv)solving (8) to obtain the velocity at time level, that is, ;(4)predicting interface with CLEAR-VOF:(i)fluid bulk advection in Lagrangian sense with available velocities atand time levels;(ii)interface reconstruction with PLIC algorithm;(iii)conducting necessary numerical corrections for the fractional volume;(5)updating the computational domain with obtained distribution of fractional volume;(6)going back to Step (2) if continue, otherwise, program stops.
In addition to the above main procedures, several important issues associated with the numerical implementation should be clarified.
(1) For the sake of numerical stabilities excessive large time step should be avoided. In this study, the time step is not only restricted to the CFL condition but also the advection of fluid bulk in the CLEAR-VOF frame. For the latter, it means that the displacement of fluid polygon in one time step is not permitted to go beyond the region covered by the mother cell and its immediate elements. To guarantee this rule, the following formula is used to automatically determine the time step increment: where is the characteristic length of the computational element, which is evaluated by with being the individual mesh area, is the absolute velocity at the mesh centre, and is a safe coefficient, normally used as 0.2 in this work.
(2) For water wave problems, it is convenient to consider only the water phase with much higher density, neglecting the inertia influence from the air, which is helpful to reduce the computer demands. Therefore the void cells are excluded in the present model. As for the full elements they are treated as usual, whereas the partial elements have to be treated with particular technique since they are generally involved in strong discontinuity in terms of the fluid properties of density and viscosity. In this work, the density and viscosity in partial/interfacial elements are averaged:
In the present model, the turbulence is not considered by means of introducing the turbulent closure. The numerical computations are carried out with very fine mesh resolution, which approaches to the two-dimensional direct numerical simulation. The above viscosity in (2) is therefore restricted to the molecular viscosity. Otherwise, if the turbulence model is incorporated, the viscosity should be interpreted as the effective viscosity.
(3) Considering the small quantity of the fluid viscosity , the integral terms related to in (2) can be neglected without introducing noticeable numerical errors. And the numerical tests of this work show that the following three terms in the pressure Poisson equation can also be neglected with little influence on the numerical accuracy; that is, , and .
(4) Generally speaking, the advection of fluid volume should be carried out throughout the computational domain for the interfacial flows. However, for some special cases of water wave problems, the evolution of free surface is definitely restricted to a certain region. For example, it is enough for us to consider only a fluid layer with the thickness of twice wave height in order to capture the free surface associated with the propagation of linear water wave. Within this subdomain, the CLEAR-VOF procedure is activated. Otherwise, out of the sub-domain, the fractional volume of fluid is directly set to unit since the computational cells are always occupied by water. By introducing the appropriate sub-domain for free surface capture, the computational efforts can be reduced considerably.
3. Numerical Validations
3.1. Dam-Breaking Flow
The classic dam breaking problem is firstly considered to validate the preset numerical model. This case mainly concerns the numerical implements of free surface capture with CLEAR-VOF. Martin and Moyce  conducted a benchmark experiment study, as shown in Figure 3, where the initial depth and breadth of the water column are with the same dimension . In the experiment the gate is removed suddenly and the water column collapses under the action of gravity . The fluid density and dynamic viscosity are and , respectively.
The computation starts from the initial still state with zero velocities and hydrostatic pressure. No-slip boundary conditions are imposed along the solid bottom and dynamic boundary condition is specified at the free surface. The numerical results of velocity field and free surface (denoted by red solid line) at several typical time instants are shown in Figure 4. The numerical results confirm that the present numerical model is capable of producing sharp predictions for the evolution of violent free surface. Especially, Figure 4 suggests that the fluid advection by CLEAR-VOF along the bottom solid wall is satisfying, not leading to the appearance of the numerical flotsam.
For the purpose of further comparisons, the developments of surface front on the bottom and water elevation on the left vertical wall with time are compared with experimental data  and available numerical results , as shown in Figure 5, in which the dimensionless parameters are used; that is, , , and . It is observed that the present numerical results agree well with both the experimental and numerical data, indicating that the present numerical model based on the three-step finite element method and CLEAR-VOF technique has good accuracy in capturing the flow interface.
3.2. Liquid Sloshing in Baffled Tank
The liquid sloshing in a rectangular container with horizontal baffles is adopted as the second case to validate the numerical model. The tank is partially filled with water with a depth of m and the tank width is m. A pair of rigid baffles with thickness 0.005 m is installed on the two opposite side-walls. The twin baffles are submerged in water with the same depth and the same breadth . Various parameters of and are considered in the computations. A sketch definition for the numerical setup is referred to in Figure 6.
The container is initially placed on a horizontal plane and then is forced to oscillate as a sinusoidal motion, which results in the following acceleration: where is the oscillatory amplitude of the tank and is the oscillatory frequency with being the natural frequency of the counterpart tank without baffles. For the present sloshing problem, it is convenient to utilize the noninertial reference coordinate, which is fixed on the moving container with the same acceleration of the container with respect to the ground. Accordingly, the governing equation of (2) should be reformulated by introducing the acceleration of the reference system into the body force term, leading to the horizontal component being
The numerical results along the right wall of the forced oscillating baffledtank by using the present non-inertial based FEM viscous fluid model are compared with those obtained from the linear and non-linear potential flow models . The comparisons are restricted to the numerical results because of the absence of experimental data. As shown in Figure 7, the present numerical results are in general agreements with the previous numerical solutions. However, the discrepancies between the results of potential flow model and the present viscous flow model are noticeable. The potential flow models produce larger amplitudes compared with the present viscous fluid model. This is because that the potential flow theory is not able to consider properly the physical dissipation, especially the strong flow shear and significant vortex shedding in the flow field due to the existence of baffles. Therefore, it can be speculated that the present viscous fluid model is able to produce more reasonable predictions for the reality than the potential flow models. The comparisons in Figure 7 suggest clearly that it is necessary to use the viscous fluid model for the water wave problems where the physical dissipations play important role.
3.3. Linear Wave Propagation and Reflection from of Vertical Wall
As mentioned previously, a numerical wave flume has to be provided with two fundamental functions, that is, the wave generation and the wave reduction. Firstly, the strategy with the collaboration of internal wave maker and sponge layer for the present numerical wave flume is validated. A typical example with wave height m and wave length m is considered. The two-dimensional numerical wave flume is set to be with a water depth of= 0.5 m and a total length of m. The corresponding wave steepness is and the ratio of wave height over water depth is 0.048, which are consistent with the Airy wave theory, giving the wave period s according to the linear dispersion relation. On the both ends of the computational domain, sponge layers are located there. The left sponge layer covers the region from −6.0 m to 0.0 m, while the right one is from 18.0 m to 24.0 m. The centroid of the source region is located at (9.0 m, 0.379 m) with height of 0.104 m and length of 0.184 m, which are in accordance with the principals suggested by Lin and Liu .
The computations are subjected to the no-slip boundary conditions along the flume bottom and the both ends of the sponge layers together with the free surface boundary condition on the wave profile. The numerical result of the wave surface history above the source region is compared with the theoretical solution, as shown in Figure 8, where the wave surface and time are normalized by the wave amplitude and wave period , respectively. It can be seen that the present numerical results agree well with the theoretical solution of linear wave theory. This comparison shows that the present wave generator works well and is able to accurately produce the target wave trains.
The performance of wave reduction with sponge layer is examined in Figure 9 by considering the wave profiles in the left sponge layer during with a time interval of. It can be seen from this figure that the wave amplitude in the dissipative zone is reduced gradually to be nearly zero. It shows that the sponge layer works fairly well. The good performance of wave damping in the right sponge layer is also observed in this study. A number of numerical trials in this work show that the desirable wave reduction can be achieved as the length of the artificial sponge layer is greater than twice the incident wave length.
The propagation of linear water wave over even bed and its reflection from vertical wall are also simulated by using the present numerical wave flume. The computational parameters include the water depth m, wave length m, wave height m, and the corresponding wave period s. The internal wave maker is adopted again in this case, located at (−8.0 m, 0.38 m) with 8 cm in length and 6 cm in height. At the left boundary of the computational domain m, a vertical wall is built. Different from the previous case, the radiation boundary is assumed at the other end of the computational domain at . Therefore this case serves as another validation focusing on the combination of internal wave maker and free radiation of water wave. Accordingly, the Sommerfeld-Orlanski boundary condition is imposed along the right boundary since the linear harmonic wave with determined phase speed is simulated in this case. On the left vertical wall and the flume bottom, the no-slip boundary condition is specified.
The typical wave profiles at, , , and are depicted in Figure 10. It can be seen that the standing wave trains are fully developed after 20 wave periods in the numerical wave flume and the numerical results remain stable solutions after the long-time calculation. It is observed that the distance between the wave node and antinode is 0.75 m, equal to 0.25, while the standing wave height is 6 cm, twice the incident wave height. These quantities are in exact agreement with the theoretical analysis, indicating that the present numerical wave flume has correctly simulated the action of linear water wave on vertical wall. Especially, the stable wave profiles are observed at the right radiation boundary, confirming the good performance of radiation boundary condition implemented in the present numerical model. Considering that the distance between the source region and vertical wall is 6.0 m, just three times of the incident wave length, the reflection waves have the same phase as the waves produced by the internal wave maker, which finally leads to the wave height in the right half domain is twice the incident wave height. However, Figure 10 shows that the wave profiles in the neighborhood of source region involve noticeable distortion. This is induced by the interactions between the reflection wave and the internal wave maker, which may lead to the occurrence of unphysical incident waves if the source region is located very close to the structure. Generally speaking, the distance more than two times of the incident wave length is acceptable.
3.4. Solitary Wave Passing Sharp Step
As a typical nonlinear dynamic process, the propagation of solitary wave and its interaction with structures have been paid much attention during the past decades. When a solitary wave travels from deep water to shallow water, it disintegrates into two or more solitons. The fission of solitary wave is often used as a standard test case for the validation of numerical wave flume.
According to the theoretical solution , the solitary wave profile can be represented by whereis the wave profile, the wave height, with being the water depth and withbeing the horizontal coordinate,is the time, and , is the wave celerity defined as
From the Eulerian point of view, the velocity components in the horizontaland verticaldirections, denoted by and, respectively, read where the dimensionlessandare introduced. The derivatives associated with the dimensionless wave profileare
The sketch definition of the numerical setup is shown in Figure 11, which is in accordance with the experiment setup in . The deep water depth is m and the shallow water depth is m. The edge of step is located at . The computational domain in the wave travelling direction covers the region from to . The incident wave amplitude is cm. Four wave gauges are arranged in the computational domain, which are located at , and, respectively.
In the computations, the no-slip boundary condition is applied on the bottom, while the free surface boundary condition is treated as (25). At the outlet , the radiation boundary condition is used. The incident wave is generated by describing the wave profile and velocity components along the inlet at , referring to (30) and (33), respectively. In order to be compatible with the desirable initial condition at with the still water state in the computational domain , the free parameter included in (30) is set to be −15.0 m together with a relative timeto be determined, which gives rise toand the nearly zero vertical velocity component at the inlet with .
The numerical results of dimensionless free surface evolution with respect to the initial still water level, normalized by the deep water depth , at four locations ( m, 7.0 m, 10.0 m and 13.0 m) are compared with the experimental observations  in Figure 12. For the purpose of further comparisons, the numerical data based on the other two different Navier-Stoke solvers by Liu and Cheng  and Shen and Chan  are also included. The former is the COBRAS model, developed by Lin and Liu  with the finite difference method, SOLA-VOF, and RANS turbulent model, while the latter is a combined immersed boundary (IB) and volume of fluid (VOF) model in the context of finite difference discretization. Different mesh resolutions for the present numerical computations have been tested and the numerical convergence was confirmed. The comparisons in Figure 12(a), for the position that the numerical wave gauge is located at m, just over the step, suggest that the present model works well. The predictions are in good agreement with the available experimental and numerical results, indicating that the target solitary wave is correctly generated and the wave propagation is accurately predicted. It can be seen from Figure 12(b) that the fission of solitary wave occurs after it passes the shelf, leading to two distinguishable solitons. The primary soliton seems to has a larger amplitude compared with that observed at Figure 12(a) while the secondary one follows the leading soliton with a lower amplitude. The comparisons with experiments and the other numerical simulations confirm the good performance of the present numerical model. Figure 12(c) presents the comparisons at the third wave gauge located at= 10 m where the separated soliton develops more significantly. Meanwhile the wave height of the primary soliton is observed to increase further. Comparing with the other two numerical results, the wave phases for the two solitons predicted by the present model are closer to the experimental observations. However, the wave heights of this work are slightly smaller than the experimental data, while the wave heights for the leading soliton predicted by [35, 36] are larger than the experiment. As for the comparisons at the last wave gauge at= 13 m, the discrepancies among the numerical results and experimental data are seen. The discrepancies between the numerical results by the three different models and the experimental data consist in both the wave phase and wave height. For the leading soliton, the results of  are in agreement with the experiment, while the model in  predicts an earlier appearance of the leading soliton and the present model produces a phase lag with respect to the experiment. As far as the phase of the second soliton is concerned, the present model gives the best result compared with the laboratory test. The discrepancies between the numerical and experimental results were speculated to be attributed to the errors in the experimental observations , which should be further clarified. Generally speaking, the present numerical model is able to produce acceptable predictions on the problem of solitary wave passing shelf with non-linear fission.
3.5. Water Wave-Induced Fluid Resonance in Narrow Gap
The advantage of volume of fluid method is the capability to deal with extremely violent free surface flows, which is further demonstrate in this section by taking consideration of the wave-induced resonance in narrow gap. It was known that the very large wave response occurs in the gaps confined by adjacent floating bodies as the incident wave frequency approaches to the resonant frequency.
Saitoh et al.  investigated the gap resonance problem by experiments. In the laboratory tests, two identical floating boxes were fixed in a wave flume. The linear water waves with constant wave height cm and various wave periods were generated. The experimental setup is shown in Figure 13, in which the water depth is, the same as the box breadth. The typical case with the same draft m for the twin bodies and a narrow gap in width cm is used as the benchmark for the numerical validation.
In the numerical computations, the internal wave maker is used for the wave generation, and two sponge layers are located at the both ends of the numerical wave flume to eliminate the reflection waves. The no-slip boundary conditions are applied along the solid walls and the free-surface boundary conditions are treated as (25). The comparisons between the numerical results and the experimental data are presented in Figure 14, where the vertical axisis the relative wave height in the narrow gap and horizontal axis represents the dimensionless wave frequency with being the wave number. It can be seen that the present numerical results are in fairly good agreement with the experimental data, both in the predictions of the resonant frequency and the extremely large resonant wave height. The resonant wave frequency and the resonant wave height in the gap are and , respectively.
The previous examinations based on the convectional potential flow model produced much overprediction on the resonant wave height, although the potential flow model is able to predict the resonant frequency with good accuracy [38, 39]. This is because that the potential flow model fails to consider the physical dissipation related to the fluid viscosity, vortex motion, and shear flow since the assumptions of perfect fluid and irrotational flow are involved. It is believed that the mechanical energy dissipation plays a very important role as the resonance occurs in the narrow gap. In order to demonstrate this aspect, the velocity field in the narrow gap during one period of fluid oscillation is provided in Figure 15, for the typical case with, , namely, the resonant condition. Figure 15 shows clearly that the complicated vortex structures are fully developed in the gap together with the significant vortex shedding from the sharp corners, which is expected to be accounted for the majority of energy dissipation . The numerical results presented above suggest that the numerical wave flume developed in this work has good performance in predicting the complex problems associated with wave-structure interaction. And the development and application of the numerical wave flume based on viscous fluid theory become absolutely necessary for the cases with significant physical dissipations.
A comprehensive two-dimensional viscous numerical wave flume is developed in this work. The present numerical model is featured by the flow solver with upwind finite element solution of the Navier-Stokes equations, interface predictor with CLEAR-VOF method, wave generator with internal wave maker and describing boundary conditions at inlet, and wave absorber with sponge layer and applying Sommerfeld-Orlanski boundary condition at outlet. The collaboration of FEM and VOF in consistence with the irregular mesh partition for the complex water wave problems and the various combinations of wave generations and wave absorptions are addressed. The 2D NWF developed in this work is validated against available theoretical, experimental, and numerical results, including the cases of dam-breaking flow, liquid sloshing in baffled tank, linear water wave propagation and reflection from vertical wall, non-linear solitary wave fission passing shelf, and water wave-induced fluid resonance in narrow gap confined by floating structures. Good agreements were obtained, confirming that the present numerical wave flume performs well in predicting complex interface flows and water wave interaction with structures. Comparing with the other existing viscous numerical wave flumes, the present model has advantage in dealing with the problems with complex boundary configuration since the finite element method with unstructured mesh partition is adopted, which further allows the more complicated fluid-structures interaction to be simulated in the unified frame of the finite element method. And the free surface capture by using the CLEAR-VOF method, inherently consistent with the irregular computational meshes, is rather simple in the numerical implement, which can be fulfilled by means of the straightforward computational geometry techniques. Several typical numerical examples show that the necessity of using the numerical wave flume based on viscous fluid model for the situations involving significant dissipative effects should be paid special attention.
This work is supported by the programs from NSFC of China with Grant nos. A020317 and 51279029. The financial support from the Fundamental Research Funds for the Central Universities (DUT13JS05) is also acknowledged.
- O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, Butterworth-Heinemann, Oxford, UK, 5th edition, 2000.
- T. J. Chung, Finite Element Analysis in Fluid Dynamics, McGraw-Hill, New York, NY, USA, 1978.
- J. Donea, S. Giuliani, and J. P. Halleux, “An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions,” Computer Methods in Applied Mechanics and Engineering, vol. 33, no. 1-3, pp. 689–723, 1982.
- M. Zhao, B. Teng, and L. Tan, “Finite element solution of wave forces on submerged horizontal circular cylinders,” China Ocean Engineering, vol. 18, no. 3, pp. 335–346, 2004.
- J. E. Welch, F. H. Harlow, I. P. Shannon, and B. J. Daly, “The MAC method: a computing technique for solving viscous incompressible, transient fluid flow problems involving free surface,” Report LA-3425, Los Alamos Scientific Report, Los Alamos, Calif, USA, 1966.
- C. W. Hirt and B. D. Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries,” Journal of Computational Physics, vol. 39, no. 1, pp. 201–225, 1981.
- S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations,” Journal of Computational Physics, vol. 79, no. 1, pp. 12–49, 1988.
- P. Lin and P. L.-F. Liu, “A numerical study of breaking waves in the surf zone,” Journal of Fluid Mechanics, vol. 359, pp. 239–264, 1998.
- H. Miyata, “Finite-difference simulation of breaking waves,” Journal of Computational Physics, vol. 65, no. 1, pp. 179–214, 1986.
- E. D. Christensen and R. Deigaard, “Large eddy simulation of breaking waves,” Coastal Engineering, vol. 42, no. 1, pp. 53–86, 2001.
- J. C. Park, Y. Uno, T. Sato, H. Miyata, and H. H. Chun, “Numerical reproduction of fully nonlinear multi-directional waves by a viscous 3D numerical wave tank,” Ocean Engineering, vol. 31, no. 11-12, pp. 1549–1565, 2004.
- P. Lubin, S. Vincent, S. Abadie, and J.-P. Caltagirone, “Three-dimensional large eddy simulation of air entrainment under plunging breaking waves,” Coastal Engineering, vol. 53, no. 8, pp. 631–655, 2006.
- P. Adams, K. George, M. Stephens, K. A. Brucker, T. T. O'Shea, and D. G. Dommermuth, “A numerical simulation of a plunging breaking wave,” Physics of Fluids, vol. 22, no. 9, Article ID 091111, 2010.
- R. Löhner, C. Yang, and E. Oñate, “On the simulation of flows with violent free surface motion,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 41–43, pp. 5597–5620, 2006.
- J. García and E. Oñate, “An unstructured finite element solver for ship hydrodynamics problems,” Journal of Applied Mechanics, vol. 70, no. 1, pp. 18–26, 2003.
- J. Garcia-Espinosa, A. Valls, and E. Oñate, “ODDLS: a new unstructured mesh finite element method for the analysis of free surface flow problems,” International Journal for Numerical Methods in Engineering, vol. 76, no. 9, pp. 1297–1327, 2008.
- N. Ashgriz, T. Barbat, and G. Wang, “A computational Lagrangian-Eulerian advection remap for free surface flows,” International Journal for Numerical Methods in Fluids, vol. 44, no. 1, pp. 1–32, 2004.
- P. Lin and P. L.-F. Liu, “Internal wave-maker for Navier-Stokes equations models,” Journal of Waterway, Port, Coastal and Ocean Engineering, vol. 125, no. 4, pp. 207–215, 1999.
- M. H. Kim, J. M. Niedzwecki, J. M. Roesset, J. C. Park, S. Y. Hong, and A. Tavassoli, “Fully nonlinear multidirectional waves by a 3-D viscous numerical wave tank,” Journal of Offshore Mechanics and Arctic Engineering, vol. 123, no. 3, pp. 124–133, 2001.
- J. Donea, “A Taylor-Galerkin method for convective transport problems,” International Journal for Numerical Methods in Engineering, vol. 20, no. 1, pp. 101–119, 1984.
- A. N. Brooks and T. J. R. Hughes, “Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations,” Computer Methods in Applied Mechanics and Engineering, vol. 32, no. 1-3, pp. 199–259, 1982.
- A. Masud and T. J. R. Hughes, “A space-time Galerkin/least-squares finite element formulation of the Navier-Stokes equations for moving domain problems,” Computer Methods in Applied Mechanics and Engineering, vol. 146, no. 1-2, pp. 91–126, 1997.
- O. C. Zienkiewicz, P. Nithiarasu, R. Codina, M. Vázquez, and P. Ortiz, “The characteristic-based-split procedure: an efficient and accurate algorithm for fluid problems,” International Journal for Numerical Methods in Fluids, vol. 31, no. 1, pp. 359–392, 1999.
- C. B. Jiang and M. Kawahara, “Analysis of unsteady incompressible flows by a three-step finite element method,” International Journal for Numerical Methods in Fluids, vol. 16, no. 9, pp. 793–811, 1993.
- C. B. Jiang and M. Kawahara, “A three-step finite element method for unsteady incompressible flows,” Computational Mechanics, vol. 11, no. 5-6, pp. 355–370, 1993.
- H. A. van der Vorst, “CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear system,” SIAM Journal on Scientific and Statistical Computing, vol. 13, no. 2, pp. 631–644, 1992.
- I. Orlanski, “A simple boundary condition for unbounded hyperbolic flows,” Journal of Computational Physics, vol. 21, no. 3, pp. 251–269, 1976.
- C. Yang, R. Löhner, and H. Lu, “An unstructured-grid based volume-of-fluid method for extreme wave and freely-floating structure interactions,” Journal of Hydrodynamics, Ser. B, vol. 18, no. 3, pp. 415–422, 2006.
- R. Löhner, C. Yang, and E. Oñate, “Simulation of flows with violent free surface motion and moving objects using unstructured grids,” International Journal for Numerical Methods in Fluids, vol. 53, no. 8, pp. 1315–1338, 2007.
- J. C. Martin and W. J. Moyce, “An experimental study of the collapse of liquid columns on a rigid horizontal plane,” Philosophical Transactions of the Royal Society of London Series A, vol. 244, pp. 312–324, 1952.
- M. S. Kim and W. I. Lee, “A new VOF-based numerical scheme for the simulation of fluid flow with free surface. Part I: new free surface-tracking algorithm and its verification,” International Journal for Numerical Methods in Fluids, vol. 42, no. 7, pp. 765–790, 2003.
- K. C. Biswal, S. K. Bhattacharyya, and P. K. Sinha, “Non-linear sloshing in partially liquid filled containers with baffles,” International Journal for Numerical Methods in Engineering, vol. 68, no. 3, pp. 317–337, 2006.
- J. J. Lee, J. E. Skjelbreia, and F. Raichlen, “Measurement of velocities in solitary waves,” Journal of the Waterway, Port, Coastal and Ocean Division, vol. 108, no. 2, pp. 200–218, 1982.
- F. J. Seabra-Santos, D. P. Renouard, and A. M. Temperville, “Numerical and experimental study of the transformation of a solitary wave over a shelf or isolated obstacle,” Journal of Fluid Mechanics, vol. 176, pp. 117–134, 1987.
- P. L.-F. Liu and Y. Cheng, “A numerical study of the evolution of a solitary wave over a shelf,” Physics of Fluids, vol. 13, no. 6, pp. 1660–1667, 2001.
- L. Shen and E.-S. Chan, “Application of a combined IB-VOF model to wave-structure interactions,” Applied Ocean Research, vol. 32, no. 1, pp. 40–48, 2010.
- T. Saitoh, G. P. Miao, and H. Ishida, “Theoretical analysis on appearance condition of fluid resonance in a narrow gap between two modules of very large floating structure,” in Proceedings of the 3rd Asia-Pacific Workshop on Marine Hydrodynamics, pp. 170–175, Shanghai, China.
- B. Li, L. Cheng, A. J. Deeks, and B. Teng, “A modified scaled boundary finite-element method for problems with parallel side-faces. Part II. Application and evaluation,” Applied Ocean Research, vol. 27, no. 4-5, pp. 224–234, 2005.
- L. Sun, R. Eatock-Taylor, and P. H. Taylor, “First- and second-order analysis of resonant waves between adjacent barges,” Journal of Fluids and Structures, vol. 26, no. 6, pp. 954–978, 2010.
- O. M. Faltinsen, O. F. Rognebakke, and A. N. Timokha, “Two-dimensional resonant piston-like sloshing in a moonpool,” Journal of Fluid Mechanics, vol. 575, pp. 359–397, 2007.