Royal Society Publishing

A free-surface capturing method for two fluid flows with moving bodies

L Qian, D.M Causon, C.G Mingham, D.M Ingram


A two-fluid solver has been developed for flow problems with both moving solid bodies and free surfaces. The underlying scheme is based on the solution of the incompressible Navier–Stokes equations for a variable density fluid system with a free surface. The cut cell method is used for tracking moving solid boundaries across a stationary background Cartesian grid. The computational domain encompasses fully both fluid regions and the fluid interface is treated as a contact discontinuity in the density field, which is captured automatically without special provision as part of the numerical solution using a time-accurate artificial compressibility method and high resolution Godunov-type scheme. A pressure-splitting algorithm is proposed for the accurate treatment of the normal pressure gradient at the interface in the presence of a gravity term. The Cartesian cut cell technique provides a highly efficient and fully automated process for generating body fitted meshes, which is particularly useful for moving boundary problems. Several test cases have been calculated using the present approach including a moving paddle as a wave generator and the initial stages of entry into still water of rigid wedges. The results compare well with other theoretical results and experimental data. Finally, test cases involving the entry into water and subsequent total immersion of a two-dimensional rigid wedge-shaped body as well as the inverse problem of wedge egress have been calculated to demonstrate the ability of the current method to tackle more general two fluid flows with interface break-up, reconnection, entrapment of one fluid into the other, as well as handling moving bodies of complex geometry.

1. Introduction

Numerical simulations of flow problems with a moving free surface separating two immiscible fluids, which may involve interface break-up and/or reconnection, or entrapment of one fluid into the other, are complicated and difficult, but of huge significance to many engineering sectors such as aerospace, civil and process engineering. The situation is further complicated if the boundaries are in motion or moving bodies are present in the fluid system.

Traditional approaches to flow problems with free surfaces are surface-fitting methods and surface-tracking methods which treat the free surface as a sharp interface whose motion is followed explicitly. Surface-fitting methods (Glimm et al. 1988; Lemos 1992; Farm et al. 1994) solve the flow equations in the liquid region only with the free surface treated as a moving upper boundary of the computational domain. This method is very efficient for simple free surface problems, but its validity is tested by the skewness of the resulting grid. If the free surface becomes highly distorted a new mesh would have to be generated to maintain the accuracy of the solution. The other limitation is that cases of wave breaking and fluid–fluid interactions such as the trapping of one fluid into the other cannot be predicted, as in this case no explicit boundary conditions can be specified at the interface. Surface-tracking methods, however, investigate both fluid regions on a fixed grid system with the free surface identified by a marker function such as the volume fraction in the widely used volume-of-fluid (VoF) method (Hirt & Nichols 1981; Youngs 1982). In this method, a transport equation is solved for the volume fraction at each time step whereupon the shape of the free surface is reconstructed explicitly using the distribution of the volume fraction function. This method can define sharp interfaces and is robust. However, tracking and reconstruction of free surfaces remains complicated and difficult, especially in three dimensions (Ubbink 1997) (unpublished work).

More recently, another approach, referred to here as the free surface capturing method, has been proposed initially by Kelecy & Pletcher (1997) and Pan & Chang (2000), based on the artificial compressibility method with high resolution Riemann solvers. In this method, the free surface is treated as a contact discontinuity in the density field in a manner similar to shock-capturing in compressible flow. The moving interface is captured automatically as part of the numerical solution along with other flow variables such as pressure and velocity. The enforcement of a conservation law eliminates the need for complex surface tracking and reconstruction procedures. Difficulties encountered by previous methods such as interface break-up, reconnection and entrapment of one fluid into the other can be dealt with relatively easily. The other advantage of this method is that its extension to three-dimensional flow problems is straightforward. The method is robust and simple and it has since been developed by the current authors for hydraulic flow problems (Qian et al. 2001). Recent versions of VoF methods such as that of Lafaurie et al. (1994) or Ubbink (1997), (unpublished work) have abandoned the explicit interface reconstruction phase instead treating the free surface as a transitional layer between the two fluids. Special solvers independent of the main flow solver are used to maintain the sharpness of the interface. These solvers take advantage of the stability of upwind schemes and the anti-diffusive property of downwind schemes blending them to yield a discretization which is both bounded and has lower levels of numerical diffusion. The VoF methods of Lafaurie et al. (1994) and Ubbink (1997), (unpublished work) along with the level set method in which the zero level set of a smooth function, e.g. the signed distance from the interface, is used to locate the free surface (Sussman et al. 1994), are sometimes, therefore, also categorized as interface capturing methods.

The coupling of moving free surface and moving boundaries or bodies in motion arises in many engineering problems. However, owing to the complexity of handling a moving and possibly breaking free surface and/or the lack of a suitable grid system for accommodating moving boundaries/bodies, many applications of free surface methods have only been to simpler fluid systems with stationary bodies. In this paper, the free surface capturing method (Kelecy & Pletcher 1997; Qian et al. 2003) is extended to free surface flows with moving boundaries or moving solid bodies using the Cartesian cut cell approach (Causon et al. 2000; Qian et al. 2003). As an alternative to structured and unstructured grid generation methods, the Cartesian cut cell method provides a fully automated process for generating a boundary-fitted mesh. Grid generation in the classical established sense is eliminated. Solid regions are simply cut out of a stationary background Cartesian mesh, their boundaries represented by different types of cut cells so as to define a boundary-fitted mesh. Moving boundaries or internal moving bodies are accommodated easily by up-dating a small number of cut cells locally at the moving boundary as long as it remains in motion. There is no requirement to re-mesh globally or even locally as with other methods. The developed cut cell solver thus has the capacity to tackle both free surfaces and moving bodies and can, therefore, be applied to a variety of engineering problems. As the flow equations are solved in both fluids, complex topographical distortions of the free surface such as breaking and recombinations are admitted. An alternative cut cell VoF projection model applied to flows with static embedded solid boundaries is described by Almgren et al. (1995, 1997). Trebotich & Colella (2001) describe an extension of a projection method for moving deformable boundaries. In the following sections, the mathematical formulation of the authors' method is reviewed, followed by a detailed description of the numerical implementation of the interface capturing scheme on a Cartesian cut cell grid with moving bodies. Improvements to the original method include a novel pressure-splitting scheme to better calculate the normal pressure gradient in the presence of the gravity term. Several test cases including wave generation in a wave tank by a moving paddle, as well as initial water entry (slamming) of two-dimensional wedges are calculated to validate the new cut cell method. Test cases involving initial water entry and total immersion of a two-dimensional rigid wedge-shaped body and the inverse problem of its egress from water to air are also used to demonstrate the potential of the current method to tackle general two fluid flows with interface break-up, reconnection, entrapment of one fluid into the other, as well as moving bodies of complex geometry. Finally, some conclusions are drawn from the study, and a brief discussion of future work is given.

2. Numerical method

(a) Governing equations and boundary conditions

The integral form of the two-dimensional incompressible Euler equations for a fluid system with variable density field can be written asEmbedded Image(2.1)where Ω is the domain of interest or control volume, S is the boundary surrounding Ω, n is the unit normal to S in the outward direction, Q is the vector of conserved variables, F is the vector of flux function through S and B is the source term for body forces. For two-dimensional problems, this equation set consists of four scalar equations which, in turn, are the mass conservation (density) equation, x-direction momentum equation, y-direction momentum equation and an incompressibility constraint. It can be shown that the density equation is in fact the basis for deriving the volume fraction equation in the well known VoF method (Hirt & Nichols 1981; Ubbink 1997) (unpublished work). By using the artificial compressibility method (Chorin 1967; Soh & Doodrich 1988; Beddhu et al. 1994) (i.e. by introducing a fictitious time derivative of pressure into the constraint equation for incompressibility) and assuming the only body force is the gravity, Q, F and B are given as Q=[ρ, ρu, ρv, p/β]T, F=fInx+gIny and B=[0,0,−ρg,0], where fI=[ρ(uub), ρu(uub)+p, ρv(uub), u]T, gI=[ρ(vvb), ρu(vvb), ρv(vvb)+p,v]T, nx and ny are the unit vectors along x- and y-directions, respectively, u and v are the flow velocity components, ρ is the density, p is the pressure, β is the coefficient of artificial compressibility and g is the gravitational acceleration. In fI and gI, ub and vb are the velocities of the boundary S, which are zero when the boundary is stationary.

Introducing a fictitious time derivative of pressure into the continuity equation produces a system of hyperbolic equations which can then be solved by any of the recently developed upwind finite volume techniques such as the characteristics-based Godunov-type schemes (Rogers & Kwak 1990). Clearly, from the above formulation, any meaningful solution can only be achieved when a divergence-free velocity field is recovered, i.e. ∂p/∂t→0. For a steady-state calculation this is not a problem. For unsteady flow problems, to achieve a fully time accurate solution, a divergence-free velocity must be attained at every physical time step. This can be achieved by using a dual-time stepping technique, subiterating the equations in a pseudotime domain to achieve a steady-state solution at each physical time step.

The boundary conditions implemented in the present study can be classified as

  1. Solid stationary or moving boundary: At this boundary the no-penetration condition can be applied to the velocity and the density is assumed to have a zero normal gradient. For the pressure, if the boundary is stationary or the acceleration of the boundary is small, we have ∇p=[0,−ρg].

  2. Outlet or open boundary: The pressure at this boundary is fixed and a zero gradient condition is applied to the velocity and density. This definition allows fluid to enter or leave the computational domain freely according to the local flow velocity and direction.

(b) The Cartesian cut cell mesh

The Cartesian cut cell mesh is an alternative to the traditional structured or unstructured body-fitted mesh. The approach is conceptually quite simple and in fact there is no grid generation in the conventional sense. This is replaced by straightforward calculations for the boundary segment intersections with a stationary background Cartesian grid. The intersection points between the boundary contour and the background grid define cut cells which have one of their sides coincident with a boundary segment so that the resulting cut cell grid is boundary-fitted. If the flow domain is multiply-connected, i.e. there is more than one solid region, these regions are specified in a similar manner with additional sets of data points. The majority of the flow domain is overlaid with a regular Cartesian mesh. Moving flow boundaries or bodies in relative motion within the flow domain are accommodated by recomputing local cell-boundary intersections at the boundaries that are in motion. Unlike conventional grid generation methods where local or even global remeshing would be required only those cells cut by the moving boundary are affected. Furthermore, the amplitude of the boundary motion is unrestricted and solid objects may traverse the flow domain. In practice, instances may arise where a cut cell becomes arbitrarily small which would obviously adversely impact on the accuracy of the solver. However, this can be dealt with quite simply using a consistent cell-merging procedure. In areas where the mesh resolution is too coarse accurately to resolve local geometric or flow features, the cut cell method lends itself well to established methods of mesh adaptation using cell subdivision techniques such as quadtree/octree in two/three dimensions or adaptive mesh refinement (AMR). A detailed description of the principles of the cut cell method including the numerical procedures applied at solid boundaries that are either static or in relative motion have been given previously by the authors (Yang et al. 1997; Causon et al. 2000; Causon et al. 2001; Qian et al. 2001) including extensions of the cut cell algorithms to three dimensions (Yang et al. 2000) and the use of mesh adaptation. Further useful discussion of techniques to handle the ‘small cut cell problem’ can be found in Almgren et al. (1995, 1997).

(c) Numerical solution

(i) Finite volume discretization of the flow equations

In the present study, a cell centred finite volume approach has been adopted to discretize the governing equations on Cartesian cut cell grids, see figure 1. For each control volume (i,j), equation (2.1) can be written asEmbedded Image(2.2)where Qi,j are the average quantities at cell (i,j) stored at the cell centre, and ∂Ci,j and Ai,j denote the boundary of the cell and area of cell i,j, respectively. The surface integration on the right side of equation (2.2) is evaluated by summing the flux vectors over each edge of a cell, and the discrete form of the integral isEmbedded Image(2.3)where m is the number of the edges of cell (i,j), Fk is the numerical flux through edge k of cell (i,j), and Δlk is the length of the edge. In the present study each fluid cell has four edges (m=4), and each cut cell has three to five edges as shown in figure 1.

Figure 1

Calculation of gradients at a cut cell.

To evaluate the inviscid numerical fluxes Embedded Image, Roe's flux function is adopted locally at each cell edge, assuming a one-dimensional Riemann problem in the direction normal to the cell edge, as follows:Embedded Image(2.4)where Embedded Image are the reconstructed right and left states at edge k of cell (i, j) and A is the flux Jacobian evaluated by Roe's average state. The quantities R, L and Λ are the right and left eigenvectors of A and the eigenvalues of A. The flux Jacobian A isEmbedded Image(2.5)and its eigenvalues are given byEmbedded Image(2.6)where Embedded Image.

(ii) Calculation of gradients and reconstruction technique on moving boundaries

To achieve second-order accuracy, a piecewise linear representation of the cell variables must first be reconstructed from the solution data before the two Riemann states at each cell edge are computed. This can be applied directly to the conserved variables Q or other characteristic variables. In the present case, it is advantageous to employ the primitive variables, U=[ρ,u,v,p]T. Choosing the primitive variables allows a highly compressive limiter to be applied to the density to help minimize numerical diffusion at the interface, while less compressive limiters or no limiting can be applied separately to the velocity components and pressure. For a given cell with centre (i,j), for example, this requires the construction of the cell variables in the formEmbedded Image(2.7)where r is the position vector from the cell centre to any point (x,y) within cell (i,j), Ui,j is the stored cell centre data, and ΔUi,j is the gradient of solution data at cell (i,j). The gradient calculation is straightforward on flow cells away from a solid boundary. First, a central difference gradient of primitive variables is estimated using the known values at neighbouring cell centres. Then, a slope limiter function is applied to prevent spurious over- or under-shoots (unbounded solutions) and/or to maintain a sharp interface. For cells near a solid boundary, however, a different gradient calculation is needed (Causon et al. 2001) to take the boundary condition for a moving body into account. If reflection boundary conditions are used on a solid boundary, the variables in the fictional cell R can be obtained by (see figure 1)Embedded Image(2.8)Embedded Image(2.9)Embedded Image(2.10)where v=unx+vny and vb=ubnx+vbny. The gradient on cut cell (i,j) may be of two types: fluid and solid. We calculate the fluid gradients and solid gradients separately, i.e.Embedded Image(2.11)Embedded Image(2.12)where Δxi+1/2,j=xi+1,jxi,j, Δyi,j+1/2=yi,j+1yi,j andEmbedded Image(2.13)Embedded Image(2.14)where Δxi,R=xRxi,j, Δyj,R=yi,jyR and G is a slope limiter function that is used to prevent over- or under-shoots and to maintain a sharp interface. The limiter function, among others, may take one of the following forms.

The van Leer limiterEmbedded Image(2.15)The k limiterEmbedded Image(2.16)where s=sign(b) and 1≤k≤2. In the present study, unless otherwise stated, the k limiter has been used for the calculations and it has been found that using k=2 (which corresponds to the superbee limiter) for density and k=1 (which corresponds to the minmod limiter) for velocity and pressure gives good performance in terms of the stability of the solution and sharpness of the interface.

Another limiter which is derived from and analogous to the high resolution interface capturing (HRIC) scheme used by Muzaferija & Peric (1998) is also included as a possible alternative to the superbee limiter for the density field. This limiter (which is referred to herein as the hyperbee limiter) is similar to the superbee limiter but is more compressive under certain flow conditions.Embedded Image(2.17)where s=sign(b) and k=2.

Once the two types of gradients have been calculated, a length average technique is used to obtain unique gradients in the cut cell:Embedded Image(2.18)Embedded Image(2.19)where Δxf=|AB|, Δxs=|BC|, Δys=|CD| and Δyf=|DE|. Δx and Δy are the uncut cell side lengths in the x and y directions, respectively. Finally we have ∇Ui,j=Uxnx+Uyny.

(iii) Pressure splitting scheme

In the current formulation, the method for calculating the normal pressure gradient term ∂p/∂y using the pressure values at the neighbouring cells is accurate for flows with constant density, but is not accurate for flows with a variable density field such as free-surface problems under the influence of gravity. Taking free surface flow problems with water and air (where the density ratio is taken as 1000:1) as a example, if both fluids are stationary, such as in a closed container as shown in figure 2, the normal pressure gradient should balance the gravity term, i.e. ∂p/∂y=−ρg. As there is a jump (discontinuity) in the normal pressure gradient across the interface due to the density gradient, calculating this gradient term using the pressure values at the neighbouring cells will introduce errors in the pressure values on either side of the interface as shown in figure 3. The interpolated pressure values using various limiters will lie at some point between A and B, where point B represents the correct pressure value at the interface. As a result, the gravity term cannot be exactly balanced by the normal pressure gradient within each cell near the free-surface. Furthermore, as the gravity term appears as a source term in the governing equations, the error may accumulate with time and have a serious impact on the accuracy and stability of the underlying two-fluid interface capturing method. In this paper, a new scheme is proposed to tackle this problem. The normal pressure gradient term is split into hydrostatic and kinematic terms asEmbedded Image(2.20)in which the hydrostatic term at each cell is calculated directly using local density values as (∇py)hydro=−ρg and a second order scheme for the kinematic pressure gradient term can be devised, for example, for cell j on a uniform grid,Embedded Image(2.21)Embedded Image(2.22)andEmbedded Image(2.23)Embedded Image(2.24)and the limited kinematic pressure gradient term can be obtained by applying equations (2.11)–(2.14) to equations (2.21) and (2.22) or equations (2.23) and (2.24).

Figure 2

Stationary water and air in a closed container.

Figure 3

Normal hydrostatic pressure distribution across the interface of water and air.

Thus, the gravity source term is exactly balanced by the normal hydrostatic pressure gradient term within each cell and, in contrast to existing implementations (Kelecy & Pletcher 1997; Qian et al. 2003), no errors will be introduced at the interface with this method of calculating the pressure gradient term. The proposed scheme for calculating the normal pressure gradient term is also general, in the sense that it can be used in any numerical method for flows with large normal density gradients and gravity effects.

To show the effectiveness of the approach, a simple test case was set up. The initial conditions for this problem are shown in figure 2, where a closed square container is half occupied by water and half by air. No disturbance is applied to the container and the water and air are only under the influence of gravity. Under these assumptions, the water and air should remain stationary. This problem is used here as a test to determine whether the numerical scheme can maintain this simple situation correctly. The computational domain was discretized with 50×50 equally spaced cells with Δxy=0.002 m. The time step used for the solution was Δt=0.000 05 s and the β value was 200. The problem was firstly calculated using the conventional method of calculating the normal pressure gradient term. The velocity vectors after 2000 time steps clearly showed that a spurious velocity field had been generated near the free surface with a large downward velocity just beneath the water surface. The same problem was recalculated using the new method for calculating the normal pressure gradient and the velocities remained identically zero after 2000 time steps. This test illustrated well the origin of errors arising from an unbalanced treatment of the normal pressure gradient and the effectiveness of the new method.

(iv) Implicit time stepping and linearization of the discretized equations

By discretizing equation (2.2) in time and omitting the subscripts for simplicity, the first-order Euler implicit difference scheme, for example, can be employed:Embedded Image(2.25)where A is the cell area. To achieve a time-accurate solution at each physical time step in unsteady flow problems, equation (2.25) must be further modified to obtain a divergence free velocity field. This is accomplished by introducing a pseudotime derivative into the system of equations, asEmbedded Image(2.26)where τ is the pseudotime and Ita=diag[1,1,1,0]. The right side of equation (2.26) can be linearized using Newton's method at the m+1 pseudotime level to yieldEmbedded Imagewhere Im=diag[1/Δτ+1/Δt, 1/Δτ+1/Δt, 1/Δτ+1/Δt, 1/Δτ]. When Δ(Qn+1)m=Qn+1,m+1Qn+1,m is iterated to zero, the density and momentum equations are satisfied, and the divergence of the velocity at time level n+1 is zero. The system of equations can be written in matrix form asEmbedded Image(2.27)where D is a block diagonal matrix, L is a block lower triangular matrix, and U is a block upper triangular matrix. Each of the elements in D, L and U is a 4×4 matrix. An approximate LU factorization (ALU) scheme as proposed by Pan & Lomax (1988) can be adopted to obtain the inverse of equation (2.27) in the formEmbedded Image(2.28)Within each time step of the implicit integration process, the subiteration is terminated when the L2 norm of the change in successive sub-iterates,Embedded Image(2.29)is less than a specified limit ϵ. In the present study ϵ=10−4.

3. Results

Two test cases have been calculated initially to show the feasibility and accuracy of the method for free surface flow problems with moving bodies. They are wave generation in a wave tank by a moving paddle and the entry into water of a two-dimensional rigid wedge. To demonstrate the applicability of the current method to complex free surface flow problems with interface break-up, reconnection, entrapment of one fluid into the other, as well as to moving bodies of relatively complex geometry, test cases involving water entry and total immersion of a two-dimensional rigid body and the inverse problem of the body's egress from water into air have also been calculated. The time step Δt used for advancing the solution was within the range of 5×10−5–5×10−4 s and the pseudotime step was normally of the order of 100 times the real physical time step in order to accelerate convergence. It was found that a β value between 200 and 2000 produced good results for the current test cases.

The density ratio of water to air is taken throughout as 1000:1 and viscous effects have been ignored in all calculations. The value of the gravitational acceleration was taken as 9.8 m s−2. The location of the free surface (air/water interface) was determined as the density contour with the average value of the two fluid densities. All the calculations presented here were performed on a 600 MHz DEC-α computer or a NEC SX-6i vector machine with 4 nodes. For a typical case using 100×100 grid points and 20 000 time steps, the computational time is around 4 h on the DEC-α machine.

(a) Wave generation in a wave tank

In wave tank tests, piston paddles oscillating horizontally are predominantly used for generating waves (wavemaker) in shallow water flow regimes. As this application involves a moving body and motion of the free surface, it is a good case to use for testing the current flow solver.

In our calculations, waves are numerically generated using a moving paddle situated at the seaward boundary, much like one used in a wave tank. The velocity of the moving paddle can be simply sinusoidal or specified as the actual velocity from the wave tank physical tests. In the first example, regular waves are generated in a 6 m long numerical wave tank with a still water depth of 0.3 m by prescribing the velocity of the paddle as Uin=−0.2 sin(2πt) (period of the wave T=1 s). A mesh with 120×40 cells was used for this simulation, the computational domain was extended 0.3 m above the still water surface into the air region and time step used was Δt=0.0001 s. The pseudo time step Δτ and β values were 0.005 s and 500, respectively. Figure 4 shows wave patterns at times t=5.0, 5.2, 5.4, 5.6 and 5.8 s during one period of development. In this figure, the moving paddle is clearly seen at the left edge of each frame and moves transversely across three grid cells.

Figure 4

Plots of the wave pattern during one period development at t=5.0, 5.2, 5.4, 5.6, 5.8 s.

In the second example, the experiments conducted by Gao (2003) were reproduced numerically. In this test, regular waves are generated in a 8.85 m long wave tank with a still water depth of 0.28 m by prescribing the velocity of the paddle to be the same as that used in the experiments as shown in figure 5 (period of the wave T=1 s). The values of the time step and β were the same as the previous test case. A total of 10 s was simulated and figure 6 shows several snapshots of the location of the free surface and velocity vectors (for both water and air) for one period of development (from t=8.5 to t=9.5 s) using Δxy=0.02 m. The free surface motion is one of the sources of vorticity creation and this can be seen in the results around the free surface contour. It can be seen from the first and last pictures that the flow pattern almost repeats itself after one period. In figure 7, the free surface elevation at three locations (x=0.55, 3.75 and 5.45 m from the initial position of the wave paddle) is compared with the corresponding experimental results (Gao 2003). To show grid convergence in the calculations, three sets of grids were used. The size of each grid was Δxy=0.04 m, Δxy=0.02 m and Δxy=0.014 m, respectively. Unsurprisingly, while the coarsest grid either under-predicted or over-predicted the peak values of the surface elevation, the values predicted on the two finer grids lie close to each other and are in satisfactory agreement with the experimental data which indicates that an accurate grid independent solution can be achieved by the present numerical method.

Figure 5

Paddle displacement versus time used in wave tank test (Gao 2003).

Figure 6

Plots of the wave pattern and velocity vectors during one period development at (a) t=8.5 s (Vmax=0.39 m s−1), (b) t=8.75 s (Vmax=0.42 m s−1), (c) t=9.00 s (Vmax=0.41 m s−1), (d) t=9.25 s (Vmax=0.41 m s−1), (e) t=9.5 s (Vmax=0.42 m s−1).

Figure 7

Comparison with experimental results for free surface elevation. (a) Position 1, x=0.55 m, (b) position 2, x=3.55 m, (c) position 3, x=5.45 m.

As expected, the particular type of slope limiter employed influences the sharpness of the interface. In figure 8, the density distribution across the interface at a distance x=2.5 m from the initial position of the paddle and time T=5.0 s is plotted for three different limiters on a grid of size Δxy=0.02 m. Of these limiters, the van Leer limiter is the most diffusive and the hyperbee limiter is the most compressive wherein the interface is only spread across 2–3 grid cells after a relatively long period of calculation. The density distribution across the interface at two different times (T=5.0 and 10.0 s) at the same position (x=2.5 m) was compared for the case of the hyperbee limiter. This showed that compressive limiters such as hyperbee and superbee have the ability to maintain a sharp interface (anti diffusive property) during long time calculations over many time steps.

Figure 8

Effect of different limiters on the thickness of the interface.

(b) Water entry of a two-dimensional rigid wedge

The initial stages of water impact (slamming) of a two-dimensional rigid wedge has been a subject of extensive studies both theoretically and experimentally (Greenhow & Lin 1983, Greenhow 1987, Zhao & Faltinsen 1993; Mei et al. 1999), owing to its practical importance in ship design and ocean engineering. As a large body of data is available from the theoretical studies, especially for the water entry of a rigid two-dimensional wedge with constant velocity, this flow problem was used as a test case in the present study.

Three deadrise angles (the angle between the wedge slope and the horizontal or undisturbed water surface), namely α=30, 45 and 81° were selected for the calculations in a square container of size 2×2 m half occupied by water and half by air.

The body was then caused to move vertically downwards to the initially calm water surface at a prescribed unit velocity V=1. In theoretical works, based on potential flow assumptions (Mei et al. 1999), it was concluded that if V.T is small (T is the duration of the impact, and thus V.T represents the penetration depth (Z) of the apex of the wedge into the water), a similarity solution can be found for this flow problem. Thus, in the early stages of water entry, the solutions are expected to be self-similar for different T if the variables are properly non-dimensionalized. The computational domain was discretized using a cut cell mesh with 160×160 cells (Δxy=0.0125 m). This grid size was chosen after conducting grid convergence tests for the α=30° case using three different grids, as shown in figure 9, where the data plotted corresponds to the time T=125 m s. The time step Δt and the pseudo time step Δτ were 0.00005 and 0.005 s, respectively. The β value was 500.

Figure 9

Grid convergence test for α=30°: pressure coefficient distribution.

On the left in figure 10, the distribution of pressure coefficient Cp=p/(1/2)ρV2 along the base of the wedge is given at three different times after impact for each of the three deadrise angles. For each deadrise angle, the numerical solutions show good correlation (self similarity). A similarity solution and analytic solution (Mei et al. 1999) are also given in each plot on this figure, which are in close agreement with the present numerical results. The corresponding free surface profiles calculated with the present numerical method at the same times during the impact are shown on the right of figure 10 and compared with the analytic solutions (Mei et al. 1999). In these figures, the apex of the wedge is located at the point (0,−1). Except for the case with α=81°, wave overturning (jet flow) is predicted where the free surface contacts the body surface. This phenomenon is confirmed by data from other studies using fully nonlinear simulations (Zhao & Faltinsen 1993) and flow visualization (Greenhow & Lin 1983). In the analytic solutions a single-valued height function is used to locate the free surface so wave overturning cannot be modelled, which accounts for the small discrepancies in the free surface profiles between the analytical and numerical results for the α=30 and 45° cases.

Figure 10

Pressure coefficient distribution along the bottom of the wedge (left) and the free surface profile (right) during the initial stage of water impact of rigid wedges. (a,b) α=30°; (c,d) α=45°; (e,f) α=81°.

In figure 11, the free surface profiles for the water entry of a rigid wedge with a deadrise angle of 30 degrees are compared with the photograph taken by Greenhow & Lin (1983), for a similar penetration depth of the wedge. Although the details of the water spray are not reproduced in the numerical simulation due to neglecting the surface tension and possible turbulence effects, the shape of the free surface and the two main jets are in close agreement with the flow visualization.

Figure 11

Comparison of free surface profiles with experiments: wedge deadrise angle=30°; number of cells: 320×320; grid size: Δxy=0.006 25 m.

(c) Water entry, total immersion and egress of a two-dimensional rigid wedge-shaped body

In this test case, firstly, the complex process of a two-dimensional wedge-shaped rigid body (α=45°) entering initially calm water with a constant prescribed velocity until totally immersed has been simulated. The geometry, initial position of the body and the calm water elevation are shown in the first frame of figure 12. The Froude number based on the downward velocity and the length of the base of the body is 0.71. A computational domain of 6 m×6 m was discretized using a uniform cut cell mesh with 120×120 cells and the time step used was Δt=0.0001 s. The pseudo time step Δτ and the β values were 0.05 s and 500, respectively.

Figure 12

Plots of the interface position and velocity vectors for a two-dimensional rigid body entering the water surface.

Figure 12 shows several snapshots of the flow development during the initial water entry and its subsequent total immersion in water. As the body pierces the surface (T=0.25 and 0.35 s), displacement waves are formed on each side of the body originating at the vertex of the wedge as the displaced water rises up the sides of the body. At the same time, the pocket of air previously trapped between the body and the water surface is forced around the sides of the wedge and separates at the corner to form a pair of symmetric vortices (T=0.5 s). A second pair of vortices is created at the rear points of detachment forming a distinct recirculating base flow region behind the body in air. As the body travels further into the water, the displacement waves created on the free surface as the body penetrates the water surface rise (T=0.75, 1.0 and 1.5 s) and subsequently overturn towards each other (T=1.75 and 2.0 s), closing the gap between them. Finally, they are reconnected (T=2.23 s), leaving air bubbles trapped behind the body and creating a jet of air directed upwards above the reconnection point on the free surface.

Secondly, the inverse problem of the above test case is also simulated, i.e. emergence of a wedge-shaped rigid two-dimensional body from beneath the water surface into air. The geometry and the initial position of the body, as well as the calm water height are shown in the top left-hand frame of figure 13. The Froude number based on the upward velocity and the length of the base of the body is 0.32. The other calculation specifics are the same as the previous test case. Figure 13 shows several snapshots of the flow development as the body emerges from beneath the water surface and subsequently progresses completely out of the water into air. Once again, these pictures show the complex flow patterns during interactions between the moving body, free surface and induced vortices and after T=2.8, the free surface has nearly returned to its equilibrium state.

Figure 13

Plots of the interface position and velocity vectors for a two-dimensional rigid body emerging from beneath the water surface.

Although there is no direct experimental data available for comparison, these computations have shown that the present method is capable of handing complex two fluid flows with moving bodies, interface break-up and reconnection, as well as entrapment of one fluid into the other. It may also serve as a test case for other numerical or experimental works.

4. Conclusions

Some two-dimensional two-fluid free surface flow problems have been studied numerically using a newly-developed interface-capturing Cartesian cut cell flow solver. These include a moving paddle used as a wavemaker in a wave tank, the initial stages of entry into water of a rigid wedge moving at a constant prescribed velocity and water entry and subsequent total immersion of a wedge-shaped rigid body and the inverse problem of its egress from water into air. For the cases of the paddle as a wavemaker and the initial stages of water entry of a rigid wedge, the comparisons of mesh independent numerical solutions with data and analytical solutions are satisfactory. The case of water penetration and immersion of a wedge-shaped rigid body and the inverse egress problem show the future promise of the method as a tool for studying free surface flow problems with interface break-up and recombination and entrapment of one fluid by another in the presence of moving bodies. The method uses an efficient dual time-stepping artificial compressibility algorithm and modern Riemann-based upwind solvers of the Godunov-type which capture the moving free surface accurately as part of the numerical solution without interface-tracking or other special provisions commonly used in VoF methods. Mesh generation in the conventional sense is eliminated in favour of defining local cut cell data on a stationary background Cartesian grid. In cases where the boundaries are in motion or moving bodies are to be included, the only required changes to the meshing procedures consist of local updates to the cut cell data at the boundaries that are in motion for as long as the motion continues. These updates involve only a small percentage of the total number of cells in the grid, the majority being unaffected. This is in contrast to conventional grid generation methods which must be re-meshed in the field at least locally. A new pressure-splitting algorithm is proposed for the accurate treatment of the normal pressure gradient at the free surface in the presence of gravity terms. This prescription can also be used with other solvers.

Future work will include a facility for local dynamic mesh adaptation to enable the method to handle geometries of arbitrary complexity. The extension of the interface capturing two-fluid solver to three-dimensional flow problems is straightforward and existing three-dimensional Cartesian cut cell algorithms can be incorporated from previous work. The method is currently in use in a study of a novel wave power device and for numerically modelling wave attack at porous seawall structures. Although the reported cases are inviscid, these other applications have necessitated the incorporation into the new method of viscous effects including turbulent transport. Results will be published shortly. A three-dimensional version of the present method will be applicable to a number of engineering problems involving free surfaces, interface break-up and recombination with static or moving boundaries such as green water overtopping of vessels and offshore structures with appropriate wave climate boundary conditions. Finally, the solver is not restricted to two fluids and the potential exists to extend it in a straightforward manner to multi-fluid applications. Also, there is no restriction that the body should be rigid.


The authors would like to thank colleague Dr F. Gao for providing the experimental data of wave generation. The data from analytic and similarity solutions for the problem of water entry of wedges were kindly provided by Dr Y. Liu and Professor Dick K. P. Yue of MIT, USA. The photograph shown in figure 11 was provided by Dr M. Greenhow of Brunel University. This work was financially supported by the EPSRC (UK) under the grants GR/N24162/01 and GR/S12333/01.

  • Received July 22, 2004.
  • Accepted June 10, 2005.


View Abstract