Arbitrary Lagrangian–Eulerian finite element analysis
of free surface flow using a velocity–vorticity formulation
D.C. Lo, D.L. Young
*Department of Civil Engineering and Hydrotech Research Institute, National Taiwan University, Taipei 10617, Taiwan Received 26 July 2002; received in revised form 13 August 2003; accepted 30 September 2003
Abstract
This paper describes the application of velocity–vorticity formulation of the Navier–Stokes equations for two-di-mensional free surface flow using an arbitrary Lagrangian–Eulerian method. The velocity Poisson equations and the vorticity transport equations are solved using a finite element method to obtain the velocity and the vorticity fields in the interior region of the computational domain. The boundary-fitted coordinates system is adopted to solve the boundary equations for kinematic and dynamic conditions at the free surface using a finite difference method. The numerical model for the velocity–vorticity formulation is validated for a square cavity flow at Re¼ 400 and 1000. The solitary wave reflected from a vertical wall is chosen as a test case for comparison and validation of the free surface flow model. Then the proposed numerical model is used to obtain flow results for the following free surface flow cases: (i) interaction between two opposite solitary waves, (ii) seiche phenomenon in a rectangular reservoir, and (iii) solitary wave through a submerged rectangular structure in a viscous fluid. The efficiency of the present numerical model for numerical treatment of free surface flows is discussed. Furthermore the advantage of this formulation with respect to primitive variables formulation is addressed from the computational point of view.
2003 Elsevier Inc. All rights reserved.
Keywords: Arbitrary Lagrangian–Eulerian; Free surface flow; Velocity–vorticity formulation
1. Introduction
Flows with free surfaces and moving interfaces find wide range of engineering applications in many interdisciplinary fields. The numerical treatment of free surface problems is highly interesting but complex because the computational domain is continuously moving. Although the underlying physics of free surface flows is understood, the numerical treatment differs with respect to handling the moving boundaries at the free surface. Irrespective of the method used to handle the free surface movement, the kinematic and dy-namic boundary conditions have to be enforced at the free surface. The important aspect of the free surface www.elsevier.com/locate/jcp
*Corresponding author. Tel.: +886-2-2362-6114; fax: +886-2-2363-9258. E-mail address:[email protected](D.L. Young).
0021-9991/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2003.09.019
flow problem is the description of the mesh movement, in other words, recalculating the coordinates of the grid points in the computational domain. This can be achieved using one of the following methods: Eulerian, Lagrangian and arbitrary Lagrangian–Eulerian.
In the Eulerian method the coordinates of the grid points are fixed throughout the computation during the transients. The main advantage of this method is that the fluid can undergo large deformations without loss of accuracy. But the great disadvantage of the Eulerian method is the difficulty in tracking the fluid interface as the fluid moves through the mesh. The most popular method in this category is the volume of fluid method (VOF) [1], which was originally developed for the finite volume method (FVM) and the marker and cell (MAC) method [2,3]. In the Lagrangian method, on the other hand the elevation of the free surface is tracked by a coordinates system that moves with the fluid particles. It qualitatively gives accurate free surface shapes and it is easy to implement the free surface boundary conditions. However, the mesh may be overly twisted due to the change of fluid domain, which may introduce numerical errors. As far as the Lagrangian approach is concerned, the implementation of this scheme in three dimensions is very difficult indeed. The Lagrangian formulation is well studied by many investigators for solving fluid flow problems [4,5]. The advantages of both the Eulerian and the Lagrangian methods are combined to yield the arbitrary Lagrangian–Eulerian (ALE) technique, in which the grid points can be moved independently of the fluid motion. In the ALE method, the nodes in the interior domain are displaced in an arbitrarily prescribed way to obtain a mesh of proper shape and thus to avoid mesh crossing. The ALE method can be considered as a moving deformable control volume problem, in which, the node velocity is different from the fluid velocity. The ALE formulation is well developed by various investigators for solving moving boundary problems [6–11].
Numerical simulations of viscous free surface flows are well established for the primitive variables pressure–velocity formulation and results have been obtained using various numerical schemes such as the finite difference method (FDM) [12,13], the finite element method (FEM) [14–19] and the FVM [20,21]. The velocity–vorticity form of the Navier–Stokes equations, pioneered by Fasel [22], has been established as an alternate to the pressure–velocity formulation for the solution of incompressible viscous flow problems. The pressure field can be eliminated by taking the curl of the Navier–Stokes equations. The velocity–vorticity formulation is well developed and established for viscous incompressible flows by many researchers for the case of fixed boundary problems [23–27]. Traditionally, most of the investigators have solved the free surface flow problems using the governing equations either in primitive variable form or vorticity– streamfunction form [28,29]. To the best of our knowledge, the present study is the first attempt to solve free surface flows using velocity–vorticity formulation. Moreover, with the use of velocity–vorticity for-mulation, the computations of the pressure fields both in the interior regions and at the free surface of the computational domain are completely eliminated. Thus the specification of boundary conditions for pressure, which is not an explicit variable in the governing equations, is not required.
Flows with free surface find an important place in many engineering applications [30–33]. Moreover, many numerical schemes for potential flows with fully non-linear free surface boundary conditions are available in the literature [34–36]. In order to comprehend the complexity of free surface flow including the viscous effect, a more general approach for handling the viscous free surface flow is needed. However the main difficulties associated with the free surface flows are: (a) the governing equations and the boundary conditions are highly non-linear, (b) difficulty in maintaining a good control over the moving mesh in the presence of the ALE method arises and (c) the vorticity boundary conditions and the free surface boundary conditions are not known a priori. These issues pose great challenge for numerical treatment of the free surface flows. In our computational methodology for unsteady incompressible viscous fluid flows, the fractional step methods have been found to be effective. In the present study, a new computational procedure is developed to solve the incompressible viscous free surface flow problems using the velocity–vorticity form of the Navier–Stokes equations. At each time step the free surface position and the velocity are satisfied with the free surface kinematic boundary condition (FSKBC) and
the free surface dynamic boundary condition (FSDBC). During the free surface flow simulation, the mesh update for the interior domain is controlled by the boundary positions by adopting an iterative procedure. Using a boundary-fitted coordinates system, and a moving mesh technique, the mesh update for the interior domain can be easily implemented. As far as the accuracy is concerned, the second order Crank– Nicolson FDM is used to solve the free surface boundary conditions. Thus the use of FDM has enabled us to make use of the boundary fitted coordinates system to track the free surface positions easily. As far as the interior domain is concerned, the FEM is found to be an effective tool to compute the flow variables at the displaced positions of the grid points in the interior domain using the ALE scheme. Thus the present numerical procedure exploits the advantages of the FDM at the free surface and the FEM in the interior domain.
The contents of this paper are organized as follows: Section 2 presents the governing equations for the free surface flow; the main variables of the equations are velocity, vorticity and the free surface elevation. In Section 3, the numerical scheme to solve the governing equations using a semi-implicit FEM is discussed. The implementation of boundary-fitted coordinates system to track the moving boundary is also men-tioned. Section 4 outlines the iterative numerical solution procedure for the entire computational domain. The validation of the present numerical model and the results obtained for different cases of free surface flow problems are delineated in detail in Section 5. The main conclusions of the present study are given in Section 6.
2. Governing equations
2.1. Incompressible viscous flow
The governing equations for a viscous incompressible fluid are given by the Navier–Stokes equations. Let X be the fluid domain, which is surrounded by a piecewise smooth boundary C. By using the vector notation, the corresponding non-dimensional basic equations of conservation of mass and momentum are [37,38] Continuity r ~u¼ 0: ð1Þ Momentum o~u otþ ~u r~u¼ rp þ 1 Rer 2~u 1 Fr2~j; ð2Þ
where ~u– the velocity vector, p – the pressure, Re – the Reynolds number, Fr – the Froude number and t – the time. Eqs. (1) and (2) represent the Navier–Stokes equations in the primitive variable (pressure–ve-locity) form.
By definition, the vorticity vector ~xcan be expressed as ~
x¼ r ~u: ð3Þ
By taking the curl on both sides of Eq. (2) and using Eqs. (1) and (3), we obtain the vorticity transport equation as follows: o~x ot þ ~u r~x¼ ~x r~uþ 1 Rer 2~x: ð4Þ
By taking the curl of Eq. (3) and using Eq. (1) we get,
r2~u¼ r ~x ð5Þ
which is the vector form of the Poisson equation for ~u.
Eqs. (4) and (5) with ~u¼ ðu; v; wÞ as velocity vector and ~x¼ ðxx;xy;xzÞ as vorticity vector are known as
velocity–vorticity form of the Navier–Stokes equations and can replace Eqs. (1) and (2) in which ~uand p are the primitive variables. Since the velocity–vorticity form of the governing equations is used for the first time to solve free surface flow problems, only a dimensional case is considered as a first attempt. The two-dimensional velocity–vorticity equations for the incompressible flows can be derived from Eqs. (4) and (5) as follows: ox ot þ u ox ox þ v ox oy ¼ 1 Rer 2x; ð6Þ r2u¼ ox oy; ð7Þ r2v¼ox ox: ð8Þ
The solution of the vorticity transport equation given by Eq. (6) along with the velocity Poisson equations given by Eqs. (7) and (8) gives the velocity and the vorticity fields in the computational domain for a given time step.
We seek a solution in the domain X, which satisfies the initial conditions, u¼ u0; v¼ v0at t¼ 0, with
no-slip boundary conditions of velocity on the solid boundary C of X. 2.2. Free surface boundary conditions
In the present study, the free surface flow is considered as a two-phase flow. It is assumed that the adjacent fluids are impermeable at the interface with constant density and viscosity. An interesting aspect of the free surface flow is that the position of the free surface boundary is not known a priori. The FSKBC can be derived from the following equation, which describes the surface that constitutes the boundary. The mathematical expression for a moving surface in two-dimension can be written as follows:
Fðx; y; tÞ ¼ y gðx; tÞ ¼ 0 ð9Þ
in which the coordinates x and y represent the horizontal and vertical directions respectively and g is the local coordinate in the vertical direction.
According to the Lagrangian statement of a material surface, a particle on the surface must remain on the surface itself. By satisfying the above kinematic condition, Eq. (9) can be rewritten as
og
ot ¼ v u og
ox; ð10Þ
where g is the height of the free surface.
The FSDBC represents the continuity of stress on the free surface and is obtained by means of a force balance. By assuming a constant pressure (pa) to act on the free surface, the dynamic boundary condition
along the free surface (a n-line) leads to dpa dn ¼ opa ox ox onþ opa oy oy on¼ 0: ð11Þ
Neglecting the surface tension effect on the free surface, the FSDBC can be obtained by combining Eqs. (2) and (11) to yield utþ uuxþ vuy 1 Reðuxxþ uyyÞ þ gxðvtþ uvxþ vvyþ 1 Fr2 1 Reðvxxþ vyyÞÞ ¼ 0: ð12Þ
This equation is valid for rotational and viscous flows.
The above Eqs. (10) and (12) have three unknown variables, two velocities (u and v) and the height of free surface (g). In order to solve the velocities (u and v) and the height of the free surface (g), the continuity equation has to be included to close the system of u, v and g on the free surface. The continuity equation can be written as follows:
ou oxþ
ov
oy¼ 0: ð13Þ
These above three Eqs. (10), (12) and (13) become one-dimensional since they are to be satisfied only at the free surface. As far as the vorticity boundary condition is concerned, it is not known a priori. In general, when the velocity field is determined, the velocity distribution at the boundary can be approximated using a polynomial function in the spatial coordinates and can be written as follows:
uðyÞ ¼ a0þ a1yþ a2y2þ a3y3þ OðDyÞ 4 ; vðxÞ ¼ b0þ b1xþ b2x2þ b3x3þ OðDxÞ 4 : ð14Þ
From the vorticity definition, x¼ov
ox ou
oy: ð15Þ
The vorticity values at the boundary can be computed using Eq. (14). Since the pressure term does not appear explicitly in the governing equations, the question about the influence of pressure on the free surface flow arises. However, the effect of pressure missing in the velocity–vorticity method may appear in the guessed boundary conditions or initial conditions. In the meantime, the effect of pressure term appears implicitly in the FSDBC in this study as shown in Eq. (12).
2.3. Arbitrary Lagrangian–Eulerian description
The ALE method is discussed in detail by Hughes et al. [9]. In this method the mesh is neither connected to the material as in the case of the Lagrangian method nor fixed to the spatial coordinates system as in the case of the Eulerian method. But the mesh can be prescribed in an arbitrary manner by defining a mesh velocity. The remapping of the state variables is required during computation since the mesh is not connected to the material. Hence the computation of the mesh velocity and remapping of the state variables are the two important steps for the ALE method. The work of Ramaswamy and Kawahara [10] can be referred for a detailed description of the application of ALE for free surface flow using the FEM.
After introducing the mesh velocities from the ALE scheme, the vorticity transport equation given by Eq. (6) becomes as ox ot þ ðu ^uÞ ox ox þ ðv ^vÞ ox oy ¼ 1 Rer 2x; ð16Þ
where u and ^u are the particle and mesh velocities in the x-direction, v and ^vare the particle and mesh velocities in the y-direction.ðu ^uÞ and ðv ^vÞ represent the convective velocities in the x and y directions respectively.
The application of ALE scheme modifies the FSKBC as og
ot ¼ v ðu ^uÞ og
ox; ð17Þ
where g is the height of the free surface. When the mesh in the x-direction is fixed, the mesh velocity ^u becomes equal to zero. Otherwise, when the boundary is moving in the x-direction, mesh velocity in Eq. (17) will have a non-zero value.
Similarly the FSDBC given by Eq. (12) is modified to utþ ðu ^uÞuxþ ðv ^vÞuy
1
Reðuxxþ uyyÞ þ gxðvtþ ðu ^uÞvxþ ðv ^vÞvyþ 1 Fr2
1
Reðvxxþ vyyÞÞ ¼ 0: ð18Þ The above non-linear equation represents the dynamic boundary condition at the free surface. The FSDBC combines the components of unsteady term, convective term, gravitational term and the viscous term.
3. Numerical formulation
As mentioned earlier, in the solution of two-dimensional free surface flow using the velocity–vorticity formulation, the vorticity transport and the Poisson-type velocity equations are solved using the FEM. The free surface boundary conditions and vorticity transport boundary condition for the free surface are solved by the FDM. In this section, both the numerical schemes are briefly described.
3.1. Finite element formulation
In the present study the computational domain is discretized using the linear triangular elements. In order to improve the accuracy of calculation and to maintain the element symmetry, the arbitrary quadrilateral element is divided into four triangular elements. The central node of a quadrilateral el-ement that is common to all the four triangular elel-ements, is eliminated by substituting the central node value in terms of the four corner node values. The stiffness equation (19) for the four triangular ele-ments can be written using the stiffness equation (20) of a quadrilateral element as shown in the fol-lowing:
Kij/j¼ fj i¼ 1; 5 j ¼ 1; 5; ð19Þ
Kij0/j¼ fj0 i¼ 1; 4 j ¼ 1; 4: ð20Þ
The weak finite element formulation of Eq. (16) is obtained by using the standard Galerkin method as
0¼ Z Z Ni ox ot þ ðu ^uÞox ox þ ðv ^vÞ ox oy 1 Re o2x ox2 þo 2x oy2 dX ð21Þ
and the above equation can be further simplified as, 0¼X Z Z NiNj dxj dt dXþ X Z Z Nkðuk h ^ukÞNi oNj ox þ Nkðvk ^vkÞNi oNj oy xjdX 1 Re X Z qnNidC 1 Re X Z Z oNi ox oNj ox þoNi oy oNj oy xjdX; ð22Þ where qn¼oxoxnxþoxoyny.
For two-dimensional flows, if qn¼ 0 or x is specified on the boundary, then the surface integral term in
the above equation vanishes. Finally Eq. (22) can be written in a matrix form as, X Mije h i xj n o þ Le ij h i xj 1 Re K e ij h i xj ¼ 0f g; ð23Þ where Mije¼ Z Z NiNjdX; Leij ¼ Z Z Nkðuk ^ukÞNi oNj ox þ Nkðvk ^vkÞNi oNj oy dX and Kije ¼ Z Z oN i ox oNj ox þoNi oy oNj oy dX; xj ¼x nþ1 j x n j Dt :
Fractional step methods have been proven to be effective schemes. In the meantime semi-implicit FEM method has rendered some advantages with respect to stability and simplicity. Eq. (23) can be written in a more general form with respect to the treatment of the time derivative term, as
1 Dt M ðeÞ ij h i þ h LðeÞij 1 ReK ðeÞ ij fxjgnþ1 ¼ 1 Dt M ðeÞ ij h i ð1 hÞ LðeÞij 1 ReK ðeÞ ij fxjgn: ð24Þ
In the present work the computations for the above equation are carried out using h¼ 0:5, which yields the second order accurate Crank–Nicolson scheme.
Similarly, the application of the GalerkinÕs weighted residual method to the velocity Poisson equations given by Eqs. (7) and (8) results in the finite element formulation for the velocity fields.
3.2. The treatment of free surface boundary equations
The equations for the free surface boundary conditions are solved in velocity–vorticity form using a FDM based on the boundary-fitted coordinates system. Any arbitrary physical space (x; y; t) can be mapped to a uniform computational space (n; g; s) with the help of the boundary-fitted coordinates system. In this study, the boundary conditions are imposed on the boundary-fitted coordinates system. However, for complicated geometries, the grid sizes Dx and Dy in the x and y directions are both uniform, and hence the grid is non-orthogonal. If the algorithm is designed for curvilinear orthogonal grids, then the non-orthogonal grids cannot be used. In the finite difference scheme, the n-direction is discretized by a central difference scheme, however the time s and space g are discretized using a backward difference scheme. The time derivative term in the FSKBC equation given by Eq. (17), is discretized using the Crank–Nicolson scheme which provides second order accuracy. As far as the FSDBC is concerned, the components of Eq. (18) can be separated into four terms as (i) time derivative term, (ii) convective term, (iii) gravity term, and (iv) the viscous term. These terms are rewritten in terms of the boundary-fitted coordinates system as follows:
(i) the time derivative (unsteady) term: ut xnþ vt yn¼ us xnþ vs yn
(ii) the convective term:
¼ ððu ^uÞuxþ ðv ^vÞuyÞxnþ ððu ^uÞvxþ ðv ^vÞvyÞyn
¼ ððu ^uÞuxþ ðv ^vÞðvx xÞÞxnþ ððu ^uÞðx þ uyÞ þ ðv ^vÞvyÞyn
¼ ðu ^uÞunþ ðv ^vÞvnþ xððu ^uÞyn ðv ^vÞxnÞ
(iii) the gravity term: 1
Fr2yn
(iv) the viscous term:
1 Reðuxxþ uyyÞ xn 1 Reðvxxþ vyyÞ yn¼ 1 Reðxy xnþ xx ynÞ ¼ 1 Re x2 nþ y 2 n J xg xnxgþ ynyg J xn
Now the final form of the FSDBC equation in the boundary-fitted coordinates system can be written as usxnþ vsynþ u h ^uunþ v ^vvn i þ xh u ^uyn v ^vxn i þ yn 1 Fr2 þ 1 Re x2 nþ yn2 J xg xnxgþ ynyg J xn ¼ 0; ð25Þ
where J¼ ðxnyg ynxgÞ is the Jacobian for the coordinate transformation.
The continuity equation given by Eq. (13) is also transformed to the boundary-fitted coordinates system as
ðygun ynugÞ þ ðvgxn vnxgÞ ¼ 0: ð26Þ
The vorticity boundary condition given by Eq. (15) can be written in terms of the transformed coordinates for the free surface boundary as
x¼ov ox ou oy¼ ðygvn ynvgÞ ðugxn unxgÞ J : ð27Þ
The free surface boundary equations are treated by using the transient boundary-fitted coordinates system. It is very important to note that the FSDBC is a fully non-linear equation. Hence, the non-linear coupling of the flow variables in the free surface boundary equations has to be treated carefully to get accurate results. This can be achieved by using an iterative solution procedure as described in Section 4.
4. Solution procedure
The governing equations for free surface flow in velocity–vorticity form are established in the previous section. The final equations to be solved for the velocity field are given by Eqs. (7) and (8). The vorticity field is obtained by solving Eq. (24). The boundary conditions for velocity are known and are specified at all the boundaries. The vorticity boundary condition at all the boundaries except at the free surface are
computed using the expression given by Eq. (15), whereas the vorticities at the free surface are computed using Eq. (27). The free surface height is computed using Eq. (17) and the velocities at the free surface are computed with Eqs. (25) and (26). The velocity Poisson equations and the vorticity transport equation are coupled and the equation governing the FSDBC is non-linear. Hence an iterative numerical solution procedure is employed to obtain the field variables and the important steps are given in the following:
For a given time step
1. Get the initial and the boundary conditions.
2. Solve the free surface height, g using FSKBC (Eq. (17)) through FDM and allocate new mesh distribu-tion.
3. Calculate the mesh velocity.
4. Combine the continuity equation (Eq. (26)) and FSDBC (Eq. (25)) to solve the unknown velocity at the free surface using FDM until convergence for the velocity is obtained.
5. Solve the velocity Poisson equations (Eqs. (7) and (8)) using FEM. 5.1. Calculate the velocity distribution at all the nodal points.
5.2. Determine the new vorticity values at the boundaries (Eqs. (15) and (27)). 6. Solve the vorticity transport equation (Eq. (24)) using FEM.
7. Check for convergence of the velocity and the vorticity components for the present iteration. In the present study the following convergence criteria are used for the flow variables.
unþ1kþ1 unþ1 k 6 106 vnþ1kþ1 vnþ1 k 6 106 xnþ1kþ1 xnþ1k 6 106
When the convergence criteria are satisfied, the iteration loop is stopped for the given time step. In the successive time step, we use the velocity and the vorticity components from the previous time step as the initial guess and use the iterative steps 2–6 to get the values for the new time step. The above solution procedure is repeated until the prescribed number of time steps is reached.
5. Numerical examples
5.1. Validation of the numerical model for a lid-driven square cavity flow
In order to verify the velocity–vorticity formulation, the lid-driven cavity flow is considered as a test case. Flow results are obtained for Re¼ 400 and 1000. The test problem consists of a square cavity of size 1:1 (dimensionless) with the top lid moving with a constant velocity of unity (dimensionless) in the x-di-rection. The cavity is assumed to be filled with an incompressible viscous fluid. In the viscous fluid flow problem, the intensity of the non-linear effect and the convective effects are related to the magnitude of the Reynolds number. For flows with high Reynolds numbers, the use of a coarse mesh often yields either non-convergent solution or a non-convergent solution without physical meaning. Furthermore, the coarse mesh is not sufficient to capture all the flow regimes accurately in a cavity flow. By refining the computational mesh, it is possible to improve the accuracy of the solution for the flow variables. A non-uniform mesh is adopted to calculate the viscous flow regimes in the present case. The boundary values of vorticity are calculated by assuming a quadratic variation of velocity normal to the wall. The computational results are obtained initially using a mesh size of 41 41 for Re ¼ 400. The velocity profiles at the vertical and the horizontal centerlines for Re¼ 400 are compared with the results of Ghia et al. [39] in Fig. 1. The present results show some deviations from the results of Ghia et al. [39] and hence the mesh is refined to 51 51. The results
obtained using the finer mesh are very close to the benchmark solution by Ghia et al. [39]. It is to be noted that Ghia et al. [39] used the vorticity–streamfunction formulation and solved their governing equations by the FDM using a multigrid numerical scheme with 129 129 mesh. But the present study used a different formulation, a different numerical scheme and a different mesh size to obtain the solution.
Fig. 2 shows the comparison of the present results obtained using a mesh size of 51 51, for the dis-tribution of velocities along the vertical and the horizontal centerlines, with the results of Ghia et al. [39]
-0.40 -0.20 0.00 0.20 0.40 0.60 0.80 1.00 0.00 0.20 0.40 0.60 0.80 1.00 Present (51* 51) Ghia et al. Botella and Peyret
u y 0.00 0.20 0.40 0.60 0.80 1.00 -0.60 -0.40 -0.20 0.00 0.20 0.40 x v Present (51* 51) Ghia et al. Botella and Peyret
(b) (a)
Fig. 2. Velocity profiles for Re¼ 1000 on vertical and horizontal centerlines (a) u ðx ¼ 0:5; yÞ, (b) v ðy ¼ 0:5; xÞ. -0.40 -0.20 0.00 0.20 0.40 0.60 0.80 1.00 0.00 0.20 0.40 0.60 0.80 1.00 Ghiaet al. .(FDM, 129*129) (51*51) u y Present (41*41) 0.00 0.20 0.40 0.60 0.80 1.00 -0.60 -0.40 -0.20 0.00 0.20 0.40 x v (b) Present Ghiaet al. .[F DM, 129*129] (51*51) Present (41*41) Present (a)
and Botella and Peyret [40] for Re¼ 1000. It can be observed that the present predictions agree well with those of Ghia et al. [39] and Botella and Peyret [40] as well. Using an entirely different formulation for the governing equations and using a different numerical scheme, our predictions agree well with the above authors for Re¼ 400 and 1000. Hence the present numerical model is validated for two-dimensional cavity flows for Re¼ 400 and 1000.
5.2. Validation of free surface flow model for solitary wave propagation
In this section, the present numerical model developed for the simulation of free surface flow is tested with the classical example of a solitary wave reflected from a vertical wall assuming an inviscid flow. The test problem considered is a solitary wave reflected from a vertical wall, for which analytical and numerical results are available in the literature [41,42]. The propagation of solitary wave is the problem of a wave returning to its initial position after colliding with the wall of the right-hand side. In the present compu-tations, a domain with finite depth, d equal to 1 and length equal to 20 is assumed. The computational domain is discretized with 16 000 linear triangular elements resulting in 4221 discrete nodes. A time step of 0.02 is used in the simulation. Calculations are carried out for different initial heights, H of wave given by the expression,H
d¼ 0:1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5. Fig. 3 shows the comparison of the present
numerical results for the maximum run-up height with the analytical solution of Byatt-Smith [41], and the numerical results of Sung et al. [42] and the experimental results of Chang [28] and Maxworthy [43]. Byatt-Smith [41] obtained the maximum run-up height analytically by the method of superposition given by the following expression: 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.00 0.40 0.80 1.20 1.60 present study
Byatt-Smith (1971):(1st order approx) Byatt-Smith (1971):(2st order approx) Sung (1999):numerical
Maxworthy (1976):exp. Chang (1997):exp.
R/d
(H: wave height , d: water depth)
(R
: Max run up-height)
H/d
Rmax d ¼ 2 H d þ 1 2 H d 2 þ o H d 3 : ð28Þ
The present numerical results agree exactly with the results of the previous investigators especially with the analytical results of Byatt-Smith [41].
5.3. Numerical results for solitary wave reflected from a vertical wall for viscous flow
After validating the present numerical models for velocity–vorticity formulation as well as for the free-surface flow problem, numerical simulation results are obtained for a solitary wave reflected from a vertical wall usingH
d¼ 0:2 . We used the central difference scheme for obtaining the results for Re ¼ 50 000 without
using upwinding. The same number of elements and nodes that was used for the case in Section 5.2 is used in this case also. Fig. 4 depicts the numerical simulation of a solitary wave propagation for time, t¼ 0–25 for Re¼ 1000, 50 000 and Fr ¼ 1 and as well as for an inviscid flow. For Re ¼ 50 000 the decay rate of the amplitude was found to be the smallest, and the wave height approaches the case of the inviscid flow. This is as expected because at such high Reynolds number values, the viscous effects become insignificant near the boundary. The wave damping is caused by the presence of no-slip solid boundary condition and the viscous effect in the fluid. The vorticity contours for Re¼ 1000 are shown in Fig. 5 during the wave propagation between the end walls. The presence of the vortex near the free surface as observed in the above figure illustrates the fact that the viscous force is confined only to the wall boundary layer and the free surface. The concentration of vorticity on the free surface is an interesting finding in our study. The use of potential flow theory, generally used in wave mechanics studies, is not capable of predicting the details of the free surface viscous flow. In the meantime, the wave height decays and the observed phase lag are attributed to the viscous damping. It is to be noted that the Reynolds number for this example plays a major role, as far as the damping characteristics of the solitary wave are concerned.
5.4. The Interaction of two opposite solitary waves
In this section, the results for collision between two opposite solitary waves using the present numerical scheme are discussed. The results are obtained assuming the fluid as inviscid. Su and Mirie [44] used a perturbation method to study the head-on collision of two solitary waves generated by an inviscid, in-compressible fluid. In a similar work, Mirie and Su [45] calculated the collision of two solitary waves to
0.00 4.00 8.00 12.00 16.00 20.00 1.00 1.20 1.40 y t=25 t=20 t=15 t=12 t=20 t=25 0.00 4.00 8.00 12.00 16.00 20.00 1.00 1.20 1.40 Re=1000 Re =50000 invisc id t=0 t=1 t=5 t=10 x y x Re=1000 Re =50000 invisc id
confirm the generation of secondary waves as predicted in their previous work [44]. They found that the solitary waves experience a phase lag change and were trailed by a dispersive wave train when the two solitary waves collide head-on. The maximum run-up height during the collision is obtained by their calibrated formula given as
0 5 10 15 x 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 y -0.6469 94 -0.2386 77 -0.1 70 625 -0.170625 -0.10 2572 -0.03 4518 9 -0.0345189 -0.0345189 Re =1000 t= 5 0 5 10 15 x 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 y -0.442836 -0.3067 3 -0.238677 -0.170-0.106252572 -0.0345189 -0.0 3451 89 -0.0345189 Re =1000 t= 10 0 5 10 15 x 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 y -0.30673 -0.238677 -0.170625 -0.102572 -0.102572 -0 .03 45 18 9 -0.0 3451 89 -0.0345189 Re =1000 t=15
Maximum run-up¼ ALþ BRþ 1 2ALBRþ 1 4 ALB 2 R þ BRA2L ; ð29Þ
where AL, BR are the initial amplitudes of the solitary waves on the left and the right-hand sides,
respectively.
We consider the collision of two solitary waves with different amplitudes. The initial amplitudes of the two solitary waves are assumed to be equal to 0.4 on the left-hand side and 0.2 on the right-hand side. The grid size used was Dx¼ 0:2, with Dt ¼ 0:02 and the domain consisted of 4221 grid points. In this moving grid, the interval of Dx is fixed but Dy is changed by the FSKBC. The interaction of the two solitary waves is shown in Fig. 6 for t¼ 0–20. The interaction is observed to be highly non-linear during the collision be-tween the two solitary waves. Therefore, the coupling of the boundary equations should be taken into consideration for the iteration loop to get a convergent solution. The maximum run-up calculated using the present numerical model is 0.641 while the value obtained using Eq. (29) is 0.652. Hence the present nu-merical model could predict the maximum run-up with less than 2% error with second order accuracy. Figs. 7 and 8 illustrate the velocity distributions during the interactive process of the two solitary waves for t¼ 5– 14. The results obtained using the present numerical model could predict well the underlying physics during the collision of the solitary waves as discussed by Su and Mirie [44]. We also found that after the collision, the solitary waves recuperate almost all of their original amplitudes in the time range of the transient considered in our computation.
5.5. Water oscillation and seiche phenomenon in a rectangular reservoir
The computation for an oscillating free surface flow is carried out for a uni-node standing wave in a rectangular basin. We assume the height of the mean water level as 5 m and length as 98 m. The initial conditions of the wave velocity are equal to zero while the wave position is represented by a curvilinear oscillating elevation given by
ys¼ A cos
px L
þ d; ð30Þ
where A is the initial amplitude, L is the length of half wave, d is depth of the mean water level, and ysis the
free surface position. The flow is considered to be inviscid at this moment so that the exact linear solution can be easily verified. The flow velocities u0, and v0 obtained from the analytical solution of the shallow water theory [46] can be written as
0.00 10.00 20.00 30.00 40.00 0.80 1.00 1.20 1.40 1.60 t=0 t=2 t=4 t=6 t=8 t=9 t=9.5 y x t=10 t=11 t=12 t=13 t=14 t=16 t=18 t=20 x y 0.80 1.00 1.20 1.40 1.60 1.80 1.80 0.00 10.00 20.00 30.00 40.00
u0¼A ffiffiffiffiffiffi gd p d sin px L sin p ffiffiffiffiffiffi gd p t L ; v0¼ A ffiffiffiffiffiffi gd p py dL cos px L sin p ffiffiffiffiffiffi gd p t L : ð31Þ
In this study, the governing equations are all dimensionless. However, in the present case, water oscil-lation in a reservoir with a dimensional scale is used. In order to convert the dimensionless quantities to the
10 15 20 25 30 x 0 0.5 1 1.5 2 y t= 5 10 15 20 25 30 x 0 0.5 1 1.5 2 y t= 8 10 15 20 25 30 x 0 0.5 1 1.5 2 y t= 10
dimensional quantities, the programming output data must satisfy the following conversion system (Table 1 shows the transformation of dimensionless to dimensional physical quantities).
Characteristic time t¼pdffiffiffiffiffiffigd¼t 0 t ) t ¼ t0 ðd=pffiffiffiffiffiffigdÞ: 10 15 20 25 30 x 0 0.5 1 1.5 2 y t=11 10 15 20 25 30 x 0 0.5 1 1.5 2 y t=12 10 15 20 25 30 x 0 0.5 1 1.5 2 y t=14
If d¼ 5 m, g ¼ 9:8 m/s2, then t¼ 1:4t0. Characteristic length H¼d h. Characteristic velocity u¼pffiffiffiffiffiffigd¼u 0 u; v ¼pffiffiffiffiffiffigd¼v0 v;
where d is the actual water level (m), and h is the non-dimensional static water level, t0is the actual time (s), t
is the non-dimensional time, u0is the actual velocity (m/s), u is the non-dimensional velocity. Similar
for-mulation is performed for v0and v.
In the present calculations, a mesh of size of 21 197 and Dt ¼ 0:02 are used. In order to verify the accuracy of the model, the fluid is assumed to have small viscosity (Re¼ 50 000, Fr ¼ 1), so that the os-cillation in the reservoir is in agreement with the linear analysis and has a little decay of amplitude as shown in Fig. 9. While the initial amplitude increases from 0.1 to 0.5 m and then from 0.5 to 2 m, the non-linear effect in the free surface flow is more significant for the latter case as shown in Figs. 10 and 11, as compared to the case of 0.1 m amplitude. Fig. 11 displays the time evolutions of free surface positions of water wave oscillation at the middle, left and right walls, respectively. Since the initial wave shape is very steep, the non-linear effect in the flow is very conspicuous.
Fig. 12 shows the moving boundary and the moving mesh used in the present numerical model for the seiche phenomenon in a rectangular reservoir for t¼ 42, 68 and 84 s. A close look at the moving com-putational mesh clearly indicates that the present numerical model has accurately coupled the flow variables obtained using the FDM at the free surface and the flow variables obtained using the FEM at the interior of the domain. This is achieved by making a one to one correspondence between the moving free surface nodes handled by the FDM and the moving mesh adopted in the FEM. The smooth displacement of the grid points along the length and height of the reservoir highlights the efficiency of the present numerical model. Fig. 13 depicts the free surface profile and the velocity vector during a surge motion in the rectangular reservoir with an initial amplitude equal to 2 m. In the present numerical simulation, it is very interesting to notice that there is no gravitational force term present in the governing equations for the interior region of
0.0 0 14.00 28.00 42.00 56.00 70.00 84.00 4.8 0 4.9 0 5.0 0 5.1 0 t (sec) y (m) Right (x=98 m) Left (x=0 m) middle
Present numerical: (non-linear effect) Linear analysis solution
98.00
Fig. 9. The time history of surface elevation for 0.1 m amplitude. Table 1
The transformation of dimensionless to dimensional physical quantities
Length (m) Time (s) Period, T (s) Velocity, u (m/s) Velocity, v (m/s)
Dimensionless 1 1 39.2 1 1
the computational domain. However the gravity term comes through the FSDBC. Due to the effect of the gravity from the dynamic condition of the free surface, the gravity wave is generated vividly by the gravitational and the inertial forces.
5.6. Solitary wave passing over a submerged structure
The interaction of a solitary wave with a submerged structure is discussed in this section. Fig. 14 il-lustrates the geometry of the computational domain and the boundary conditions of the problem. For an accurate computation of the flow variables at the rear of the solid structure, a non-uniform mesh shown in Fig. 15 is used. The initial condition of the computational domain is assumed to have a solitary wave starting from x0¼ 10 to pass over the submerged rectangular block. The total computational domain is
restricted within20 6 x 6 20. The dimensionless velocity distribution of the free surface flow, for an initial wave height of A0¼ 0:4 and Re ¼ 1000 and Fr ¼ 1 is calculated using a total number of 12 131 grid points.
In order to reduce the computational time and increase the accuracy of the solution, a non-uniform grid is adopted near the solid and at the free surface boundaries. In the upper region, the mesh distribution is changed in accordance with the free surface in the vertical direction, while in the lower region the mesh is fixed in space. When the free surface position is determined by the FSKBC, the grid in the upper region has
0.00 14.00 28.00 42.00 56.00 70.00 84.00 98.00 4.00 4.50 5.00 5.50 right(x=98 m) left(x=0 m) middle (x=49 m) t (sec) y (m)
Present numerical: (non-linear effect) Linear analysis solution
Fig. 10. The time history of surface elevation for 0.5 m amplitude.
0.00 14.00 28.00 42.00 56.00 70.00 84.00 98.00 2.00 4.00 6.00 8.00 10.00 right(x=98 m) left(x=0 m) middle x=49 m y (m) t (sec)
to be allowed to move conformably with the free surface in the vertical direction. Fig. 15 depicts the non-uniform adaptive grid distribution for the calculation of a solitary wave that passes over the submerged structure. In order to understand the vortex motion generated at the rear of the structure, the local streamline contours near the submerged structure are plotted and are shown in Fig. 16. The velocity dis-tribution of a solitary wave passing over the structure at t¼ 16 is given in Fig. 17. Fig. 18 also depicts the
0 25 x 0 1 2 3 4 5 6 7 8 9 y 42 sec x 0 1 2 3 4 5 6 7 8 9 y 68 x 0 1 2 3 4 5 6 7 8 9 y 3T = 8 4 50 75 sec 0 25 50 75 100 sec 0 25 50 75 100
gradual development of the recirculating flow regimes at the downstream of the structure and is clearly predicted by the present numerical study.
For the case of Re¼ 10 000, the vorticity is generated and confined around the rear of the submerged structure as seen in Fig. 19. The recirculating flow generated at the trailing edge of the structure is clearly observed in Fig. 20. With passage of time it is noticed that the vorticity grows in magnitude, keeping in contact with the structure. We can conclude to a first approximation that the viscosity plays an important role in the interaction between the solitary wave and the solid structure. When the solitary wave passes over
0 x 0 1 2 3 4 5 6 7 8 9 y 5T /8 = 17.5 sec x 0 1 2 3 4 5 6 7 8 9 y 3T/ 4 = 21 sec x 0 1 2 3 4 5 6 7 8 9 y 1T = 28 sec 50 25 75 0 25 50 75 0 25 50 75
the submerged structure at high Reynolds number, the vortical flow always increases on the rear surface of the structure. From the simulated results, it is observed that the present numerical model captures well the more significant eddy structure in this case. Fig. 20 displays the separated-flow pattern at different time levels during the wave–structure interaction. While the transmitted wave moves downstream for a distance away from the structure, the primary vortex is always attached to the structure and grows continuously. It is also found that the small solitary wave generated at the upstream advances from the submerged structure at t¼ 15. The highly non-linear growth of the present free surface flow has resulted in a more complex flow field with steep velocity gradients within the boundary layer due to the presence of the submerged structure. All the above results have been obtained without involving pressure as one of the flow variables because of the velocity–vorticity form of the Navier–Stokes equations used in the present study. It is well known from the available literature that the presence of the pressure term in the governing equations requires a special treatment since the pressure is not a variable in the continuity equation. Though the predictor-corrector or
x y
H=1
B=0.5
L=2
Fig. 14. A schematic diagram of the geometry and boundary conditions for solitary wave passing over a submerged structure.
-2 0 2 4 x -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 y t=9
the velocity correction method can be used to handle the pressure–velocity form of the Navier–Stokes equations, the numerical procedure adds further complexity for free surface flow problems.
5.7. Extension to three-dimensional free surface flows
The proposed numerical model can be extended to study three-dimensional free surface flow problems. Using a similar procedure followed for the two-dimensional Navier–Stokes equations, the vorticity transport and the velocity Poisson equations can be obtained from the three-dimensional Navier–Stokes equations and the definition of the vorticity. The application of the ALE scheme is also straightforward; in
-3 -2 -1 0 1 2 x -1 -0.5 0 0.5 y t=15 -3 -2 -1 0 1 2 x -1 -0.5 0 0.5 y t=16
Fig. 16. Streamline contours near the submerged structure.
0 5 x -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 y t=16 10
which case, there will be three mesh velocities to be computed. The complexity of the problem arises from the management of the moving mesh for the free surface. The boundary-fitted coordinates system can be used to deal with the movement of the two-dimensional free surface. But this will make the finite difference solution procedure at the free surface more cumbersome. An alternative will be to use an efficient mesh generator such as transfinite interpolation method to compute the mesh displacement at the free surface and the interior of the computational domain. The present research group is carrying out further research in this direction.
6. Conclusions
A numerical model for free surface flow problems is proposed using the Navier–Stokes equations in velocity–vorticity form for the first time. An effective numerical algorithm has been developed which combines the finite difference and the FEMs to solve the governing equations at the boundary and the interior of the computational domain using an ALE scheme. The numerical test results show that the velocity–vorticity formulation gives accurate solutions with 2% error for flow problems with free surface involving the interaction of two solitary waves.
• The test results obtained for a square cavity lid-driven flow are accurate and show good agreement with the bench mark solutions for Re¼ 400 and 1000.
• The test results for a solitary wave reflected from a vertical wall also show good agreement with the re-sults obtained by experimental, numerical as well as analytical solution. The maximum run-up height could be well predicted as compared to the analytical, experimental and other numerical solutions. • The use of boundary-fitted coordinates system in conjunction with the FDM is effective in predicting the
viscous and non-viscous free surface flow regimes characterized by the value of the Reynolds number. • The velocity–vorticity formulation predicts accurately the non-linear free surface wave behavior for sol-itary waves reflected from a vertical wall, for two opposite solsol-itary waves, for the seiche phenomenon in a rectangular basin and for the solitary wave passing over a submerged structure. The present numeri-cal model treats the free surface flows as fully non-linear and hence produces more accurate results compared to other investigators who assume weak non-linear wave theory.
-3 -2 -1 0 1 2 x -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 y t=16
x y -2 -1 0 1 -1 -0.5 0 0.5 1 2.71412 2.13703 1.55995 0.982865 0.405782 -0.171301 -0.748385 -1.32547 -1.90255 -2.47964 -3.05672 -3.6338 -4.21089 -4.78797 -5.36505 -5.94214 t=12 x y -2 -1 0 1 -1 -0.5 0 0.5 1 3.83821 3.1932 2.54819 1.90318 1.25817 0.613155 -0.0318555 -0.676866 -1.32188 -1.96689 -2.6119 -3.25691 -3.90192 -4.54693 -5.19194 -5.83695 t=15 x y -2 -1 0 1 -1 -0.5 0 0.5 1 3 2.12425 1.23107 0.5 0.337889 0.2 -0.1 -0.2 -0.5 -1 -1.44847 -2.34166 -3.23484 -4.12802 -5.0212 -5.91439 t=20
• The numerical predictions of surface elevations for an oscillating wave in a rectangular basin exactly co-incide with the analytical solutions for the linear wave with small amplitude 0.1m. However the non-lin-ear wave prevails when the amplitude is increased from 0.5 to 2 m.
• The boundary-fitted coordinates system used for treating the free surface movement produces a smooth variation of the moving mesh and the moving boundary. The use of non-uniform mesh is highly essential to capture the vortex generation at the rear of the solid structure in the case of wave propagation over a submerged structure.
Acknowledgements
The work reported in this article was supported by the National Science Council, Taiwan. It is greatly appreciated. We would also like to thank the reviewers for their very constructive review comments and suggestions. A special thank is extended to Dr. K. Murugesan for his instructive review of the manuscript.
x y -2 -1 0 1 2 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 t=10 x y -2 -1 0 1 2 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -1 x y -2 -1 0 1 2 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 x y -2 -1 0 1 2 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 t=12 t=15 t=20
References
[1] C.W. Hirt, B.D. Nicholas, Volume of fluid method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201. [2] F.H. Halow, J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface, Phys.
Fluids 12 (1965) 2182.
[3] J.A. Viecelli, A method for including arbitrary external boundaries in the MAC incompressible fluid computing technique, J. Comput. Phys. 4 (1969) 543.
[4] P. Hansbo, Lagrangian incompressible flow computations in three dimensions by use of space–time finite elements, Int. J. Numer. Meth. Fluids 20 (1995) 989.
[5] B. Ramaswamy, M. Kawahara, Lagrangian finite element analysis applied to viscous free surface flow, Int. J. Numer. Meth. Fluids 7 (1987) 953.
[6] C.W. Hirt, A.A. Amsden, J.L. Cook, An arbitrary Lagrangian Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 227.
[7] R.K.C. Chan, A generalized arbitrary Lagrangian–Eulerian method for incompressible flows with sharp interfaces, J. Comput. Phys. 17 (1975) 311.
[8] W.E. Pracht, Calculating three-dimensional fluid flows at all speeds with an Eulerian–Lagrangian computing mesh, J. Comput. Phys. 17 (1975) 132.
[9] T.J.R. Hughes, W.K. Liu, T.K. Zimmerman, Lagrangian Eulerian finite element formulation for incompressible viscous flows, Comput. Meth. Appl. Mech. Eng. 29 (1981) 329.
[10] B. Ramaswamy, M. Kawahara, Arbitrary Lagrangian–Eulerian finite element method for unsteady convective incompressible viscous free surface flow, Int. J. Numer. Meth. Fluids 7 (1987) 1053.
[11] M. Souli, J.P. Zolesio, Arbitrary Lagrangian–Eulerian and free surface methods in fluid mechanics, Comput. Meth. Appl. Mech. Eng. 191 (2001) 451.
[12] D.H. Zhang, A.T. Chwang, Numerical study of nonlinear shallow water waves produced by a submerged moving disturbance in viscous flow, Phys. Fluids 8 (1996) 147.
[13] D.H. Zhang, A.T. Chwang, On solitary waves forced by underwater moving objects, J. Fluid Mech. 389 (1999) 119.
[14] H.G. Choi, A study on segregated finite element algorithms for the Navier–Stokes equations, Ph.D. Thesis, Seoul National University, 1996.
[15] D.L. Young, Y.S. Liu, Finite element analysis of free surface flows, in: W.L. Hogarth, B.J. Noye (Eds.), Computational Techniques and Applications, Hemisphere Pub, Washington, DC, USA, 1990, p. 461.
[16] A. Masud, T.J.R. Hughes, A space–time Galerkin/least-squares finite element formulation of the Navier–Stokes equations for moving domain problems, Comput. Meth. Appl. Mech. Eng. 146 (1997) 91.
[17] S.E. Navti, K. Ravindran, R.W. Lewis, D.C. Taylor, Finite element modeling of surface tension effects using a Lagrangian– Eulerian kinematic description, Comput. Meth. Appl. Mech. Eng. 147 (1997) 41.
[18] X. Cai, H.P. Langtangen, B.F. Nielsen, A. Tveito, A finite element method for fully nonlinear water waves, J. Comput. Phys. 143 (1998) 544.
[19] S.R. Idelsohn, E. Onate, C. Sacco, Finite element solution of free-surface ship-wave problems, Int. J. Numer. Meth. Eng. 45 (1999) 503.
[20] P.J. Zwart, G.D. Raithby, M.J. Raw, The integrated space–time finite volume method and its application to moving boundary problems, J. Comput. Phys. 154 (1999) 497.
[21] B.T. Helenbrook, L. Martinelli, C.K. Law, A numerical method for solving incompressible flow problems with a surface of discontinuity, J. Comput. Phys. 148 (1999) 366.
[22] H. Fasel, Investigation of the stability of boundary layers by a finite-difference model of the Navier–Stokes equations, J. Fluid Mech. 78 (1976) 355.
[23] D.L. Young, S.K. Yang, T.I. Eldho, Solution of the Navier–Stokes equations in velocity–vorticity form using a Eulerian– Lagrangian boundary element method, Int. J. Numer. Meth. Fluids 34 (2000) 627.
[24] D.L. Young, Y.H. Liu, T.I. Eldho, A combined BEM–FEM model for the velocity–vorticity formulation of the Navier–Stokes equations in three dimensions, Eng. Anal. Bound. Elem. 24 (2000) 307.
[25] G. Guj, F. Stella, A vorticity–velocity method for the numerical solution of 3D incompressible flows, J. Comput. Phys. 106 (1993) 286.
[26] C.J. Ho, F.H. Lin, Numerical simulation of three-dimensional incompressible flow by a new formulation, Int. J. Numer. Meth. Fluids 23 (1996) 1073.
[27] H.J.H. Clercx, A spectral solver for the Navier–Stokes equations in the velocity–vorticity formulation for flows with two nonperiodic directions, J. Comput. Phys. 137 (1997) 186.
[28] J.H. Chang, Interaction of solitary waves with structures in viscous fluid, Ph.D Thesis, National Cheng Kung University, Tainan, Taiwan,1997.
[30] T.Y. Wu, Generation of upstream advancing solitons by moving disturbances, J. Fluid Mech. 184 (1987) 75.
[31] C.C. Mei, H.S. Choi, Forces on a slender ship advancing near the critical speed in a wide canal, J. Fluid Mech. 179 (1987) 59. [32] X.N. Chen, S.D. Sharma, A slender ship moving at a near-critical speed in a shallow channel, J. Fluid Mech. 291 (1995) 263. [33] C. Katsis, T.R. Akylas, On the excitation of long nonlinear water waves by a moving pressure distribution, Part 2:
Three-dimensional effects, J. Fluid Mech. 177 (1987) 49.
[34] A. Clement, Coupling of two absorbing boundary conditions for 2D time-domain simulations of free surface gravity waves, J. Comput. Phys. 126 (1996) 139.
[35] D. Ambrosi, L. Quartapelle, A Taylor–Galerkin method for simulating nonlinear dispersive water waves, J. Comput. Phys. 146 (1998) 546.
[36] I. Robertson, S. Sherwin, Free-surface flow simulation using hp/spectral elements, J. Comput. Phys. 155 (1999) 26. [37] G.K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 1967.
[38] J.A. Liggett, Fluid Mechanics, McGraw-Hill, New York, 1994.
[39] U. Ghia, K.N. Ghia, C.T. Shin, High-Re solutions for incompressible flow using the Navier–Stokes equations and a multigrid method, J. Comput. Phys. 48 (1982) 387.
[40] O. Botella, R. Peyret, Benchmark spectral results on the lid-driven cavity flow, Comput. Fluids 27 (1998) 421.
[41] J.G.B. Byatt-Smith, An integral equation for unsteady surface waves and a comment on the Boussinesq equation, J. Fluid Mech. 49 (1971) 625.
[42] J. Sung, H.G. Choi, J.Y. Yoo, Numerical study on fluid flow and heat liquid film involving a hydraulic jump, Heat Transfer 28 (1999) 18.
[43] T. Maxworthy, Experiments on collisions between solitary wave, J. Fluid Mech. 76 (1976) 177. [44] C.H. Su, R.M. Mirie, On head-on collision between two solitary waves, J. Fluid Mech. 98 (1980) 509.
[45] R.M. Mirie, C.H. Su, Collisions between two solitary waves. Part 2. A numerical study, J. Fluid Mech. 115 (1982) 475. [46] G. Neumann, W.J. Pierson, Principle of physical oceanography, Prentice Hall, New York, 1966.
