body interactions

European Journal of Mechanics B/Fluids 38 (2013) 1–17

Contents lists available at SciVerse ScienceDirect

European Journal of Mechanics B/Fluids journal homepage: www.elsevier.com/locate/ejmflu

Development of a free surface flow solver for the simulation of wave/body interactions H.B. Gu ∗ , D.M. Causon, C.G. Mingham, L. Qian Centre for Mathematical Modelling and Flow Analysis, School of Computing, Mathematics and Digital Technology, The Manchester Metropolitan University, Manchester M1 5GD, UK

article

info

Article history: Received 19 March 2012 Received in revised form 7 September 2012 Accepted 29 September 2012 Available online 12 October 2012 Keywords: Interface tracking Level set method Free surface flow Two-phase flow Partial cell Fixed and moving bodies

abstract A two-fluid free surface flow model, that can simulate wave interactions with stationary or moving bodies, has been developed. The incompressible Navier–Stokes (NS) equations are solved using a finite difference solver with an efficient particle level set method to track the interface between water and air. A semiLagrangian approach and the fast marching method have been applied to solve the level set equation and, with the aid of the locally deployed particles, to reinitialize the level set function. To represent solid stationary and moving boundaries in the flow field effectively, a partial cell treatment combined with a local relative stationary method has been implemented. A new simple level set boundary condition for both fixed and moving bodies is proposed. The free surface flow model has been validated for a number of 2D and 3D benchmark test cases for applications to complex flow situations involving the violent motions of water, air and their interfaces, as well as their interaction with fixed or moving bodies. Detailed comparisons between the numerical results and available analytical solutions or experimental data show that the model is generic, accurate and holds promise for the numerical simulation of such flow problems. © 2012 Elsevier Masson SAS. All rights reserved.

1. Introduction Numerical simulation of flows with a moving free surface separating two immiscible fluids is complicated and difficult especially those involving interface break-up, mixing or entrapment of one fluid within another and the situation is further complicated if fixed or moving bodies are present. However, this is hugely important in many applications in aerospace, civil, coastal and offshore engineering. To tackle this challenging computational problem adequately, an efficient flow solver along with robust and accurate methods for tracking the location of the free surface and for representing solid boundaries must be developed. Fortunately, advances have been made in each aspect over the past few years that can be used as a firm foundation for the present work. Essentially, three types of approach have been used most frequently for flow problems with free surface. One is the surfacefitting method [1,2], which solves the flow equations only in the liquid region with the free surface treated as a moving upper boundary. This method is very efficient for simple free surface problems but its validity is limited by the requirement that the free surface ordinate is single-valued which limits wave steepness to waves that do not overturn. An obvious implication is that water–air interactions such as the trapping of air bubbles as

∗

Corresponding author. E-mail address: [email protected] (H.B. Gu).

0997-7546/$ – see front matter © 2012 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.euromechflu.2012.09.010

waves overturn cannot be predicted. Another approach is the socalled free surface capturing method [3,4]. This method views the free surface as a contact discontinuity in the density field, and its position is captured automatically as part of the numerical solution, along with other flow variables such as pressure and velocity, in a manner similar to shock capturing in gas dynamics, by enforcing a conservation law thereby eliminating the need for complex surface tracking and reconstruction of the free surface. However, the spatial resolution of the free surface can be compromised by numerical dissipation particularly in a long calculation if the grid is not sufficiently fine around the free surface. The third approach involves surface tracking such as in early volume of fluid (VOF) methods [5] where a transport equation is solved for the volume fraction function associated with each finite volume cell at each time step. The shape of the free surface is then reconstructed from updated values of the volume fraction function. These methods have improved steadily to the point where a sharply defined free surface is generally possible and the methods are robust. However, the tracking and reconstruction algorithms are complicated and difficult, especially in three dimensions. The level set method [6,7] solves a similar transport equation to the VOF method for a level set function, which is, unlike the volume fraction, a continuous function with the zero level set defined as the free surface. This method can in principle provide an exact representation of the free surface using the zero level set at each time step whilst limiting numerical dissipation in long calculations, and has become increasingly popular owing to its ability to deal with complex distortions of the free surface including stretching,

2

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

crossing, breaking and recombination of the interface. The level set method along with some variants of the VOF method [8] can also be categorized as a free surface capturing method in the sense that no explicit reconstruction of the interface is needed. Compared with other interface tracking/capturing methods in computational fluid dynamics (CFD), such as the classical marker and cell (MAC) [9] and VOF [5] methods, one of the most significant advantages of the level set method is its ability to capture easily the full topological evolution of the interface, especially in three dimensions. However, discretization of the level set equation can lead to significant numerical dissipation, as was found with volume fraction equation in VOF methods. This usually manifests itself as a loss of mass or volume in areas of high curvature or other underresolved regions. To overcome this problem, spatially high-order difference schemes such as essentially non-oscillatory (ENO) [10] or weighted ENO (WENO) [11] schemes may be used for the advective terms including the level set equation itself with a TVD-Runge–Kutta scheme for the time discretization. In recent developments, Enright et al. [12] was one of the first to present a particle level set method in order to further enhance the resolution of the free surface calculation. In this method, by introducing massless particles around the free surface and by comparing their relative positions to the interface, the location of the interface can be corrected at each time step, especially for regions where the resolution of the mesh is not fine enough to capture all the details of the interface. Enright’s particle level set method [13] has been shown to be a method of high accuracy for the re-construction of the interface. A 5th-order HJ-WENO and 3rd-order TVD Runge–Kutta scheme were adopted in combination with massless particles to correct computational errors in the level set function. This particle level set method, however, is subject to a restrictive limitation on the allowable time step according to the CFL criterion. In order to improve the efficiency of the level set method, Strain [14] put forward a semi-Lagrangian advection algorithm. However, whilst free of the limits of the CFL condition, the method suffers from serious numerical dissipation. Thus, in the present study, and following this work, a semi-Lagrangian advection scheme coupled with a characteristic-based particle method [14,15] is adopted to accurately and efficiently track a passively advected interface. The characteristic-based particle method using the local particle level set method improves mass conservation and limits numerical dissipation. An important step in the level set method is the reinitialization of the level set function φ during a calculation, which is usually completed by an iterative solution of the Eikonal equation, discretized, for example, by the HJ-WENO scheme [16] or a Godunovtype scheme. To reduce the required CPU time for the iteration process the fast marching method [17–19] has been adopted in the current work. Although the accuracy of this scheme is only formally 1st-order, the use of particles allows more accurate numerical results to be obtained [15] and the computational time is much reduced. For the related flow solver any suitable finite difference, finite element or finite volume method can be used to discretize the Navier–Stokes equations. To handle an obstacle on a fixed-grid, there are a number of methods for the treatment of boundary conditions. When a virtual force is introduced to represent the reaction of the obstacle to the ambient flow, the method is usually termed the immersed boundary method (IBM) or fictitious domain method (FDM). Such methods eliminate the need to impose boundary conditions explicitly at a solid surface and instead virtual forces are computed. Boundary conditions may be invoked at a solid surface, using the cut-cell method where background Cartesian grid cells are simply cut by solid bodies. The cut-cell (also called ‘‘partial cell’’) method has been widely employed in finite

volume methods [20]. Within a finite difference method, which is used for the present work, the partial cell technique also renders an obstacle very easy to simulate. To simulate a moving obstacle, a moving boundary condition may be directly imposed on the obstacle surface. However, this is difficult to implement directly on a fixed-grid, particularly when the obstacle passes across the cells of the grid. Fortunately, using the local relative stationary method combined with the partial cell technique it is easy to simulate flow problems with moving bodies [21] and therefore we adopt this method in the current work where wave interactions with fixed or moving bodies are a major focus. When a solid surface comes into contact with an interface, the boundary value at the interface also presents a challenge even in a fixed grid system. Guendelman et al. [22] and Yang et al. [23] put forward methods for the treatment of the boundary value for the level set function in such cases. The method of Guendelman et al. appears complicated whilst on the other hand the method of Yang et al. may induce non-physical behavior of the interface near the contact point with solid surface. In this paper, a simple method based on one sided interpolation is described for treating the problem, which is suitable for both fixed and moving bodies. In the following sections, the essential details of the mathematical formulation of the level set free surface flow solver are outlined. More details relating to the implementation of the particle level set method can be found as indicated in the original referenced sources. Here, the main focus will be on the implementation of the authors’ finite difference flow solver that combines these elements specifically for applications to free surface flow problems relevant to study of wave interactions with fixed and floating bodies. The developed solver is firstly validated by computing a 2D problem involving the collapse of a water column. To validate the partial cell treatment of solid obstacles, the collapse of a water column and its impingement on an obstacle is then simulated in 2D followed by 3D simulations of an impulse wave interaction with a bottommounted circular cylinder and a collapsing water column and its action on a down-stream cube. Grid refinement tests are also conducted for these test cases showing that the present method is capable of producing results of comparable accuracy to alternative methods using much coarser meshes. Subsequently, a numerical wave tank (NWT) with a piston-type wave maker oscillating foreand-aft is simulated to test the ability of the present model for free surface flows with a moving body. Finally, a complex flow problem involving two moving bodies is simulated to demonstrate the robustness of the numerical model for simulating wave/moving body interactions. In the last section conclusions are drawn with reference to future extensions. 2. Numerical method 2.1. Basic governing equations Here, we are solving the flow equations in both air and water separated by a free surface. The relevant equations of motion assuming incompressible air and water can be written as

∂ ui =0 ∂ xi ∂ ui ∂ ui 1 ∂p 1 ∂τij + uj =− + gi + ∂t ∂ xj ρ (φ) ∂ xi ρ (φ) ∂ xj

(1) (2)

where i = 1, 2 for 2D and i = 1, 2, 3 for 3D. ui is the velocity component in i direction, p is pressure, ρ is density of air or water, and τij = 2µ(φ)σij , σij is defined as the rate of deformation tensor

σij =

1 2



∂ uj ∂ ui + ∂ xj ∂ xi



.

(3)

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

3

In Eq. (2) φ is the level set function defined as a signed distance function that is positive in water and negative in air, (Fig. 1). Accordingly, the level set transport equation can be derived as

In general, the tentative velocities do not satisfy the continuity equation (8). Thus, these must be updated by an iterative pressure correction method using,

∂φ ∂φ = 0. + ui · ∂t ∂ xi

θ uin+1 − θ u˜ i 1 ∂ pn+1 = −θ + θ gi 1t ρ (φ) ∂ xi

(4)

The zero level set is initialized as the position of the free surface between water and air. From the level set function, the density ρ(φ) and the viscosity µ(φ) are written as

ρ (φ) = ρa (1 − H (φ)) + ρw H (φ) µ (φ) = µa (1 − H (φ)) + µw H (φ)

 0,    φ 1 2π φ H (φ) = 0.5 + + sin ,  3h 2π 3h   1,

∂

(6)

xi

φ < −1.5h |φ| ≤ 1.5h

(7)

φ > 1.5h

in which h is the local minimum size of the spatial step. 2.2. Partial cell technique Naturally, at solid bodies inside the computational domain boundary conditions are normally required. However, to simplify the treatment of internal obstacles and boundaries the partial cell technique (PCT) is used that does not require explicit boundary conditions (even if the obstacles are a complex shape). The main idea of the PCT is to define an openness function θ as the ratio of fluid (water) volume to computational cell volume [13,24]. That is, the whole computational domain is divided into fluid cells, solid cells and partial cells. When θ = 1 the cell is a fluid cell; θ = 0 a solid cell; and for 0 < θ < 1 a partial cell, see Fig. 1. With the openness function so defined, the equations of motion for incompressible two-phase flow become

∂ (θ ui ) =0 ∂ xi ∂ (θ ui ) ∂ (θ ui ) θ ∂p 1 ∂τij + θ uj =− + θ gi + θ . ∂t ∂ xj ρ (φ) ∂ xi ρ (φ) ∂ xj



θ

∂ p n +1 ρ (φ) ∂ xi 1

 =

1 ∂θ u˜ i

1 t ∂ xi

+

∂ (θ gi ) . ∂ xi

(12)

Eq. (12) is obtained by enforcing Eq. (11) at the time step n + 1, ∂(θ uin+1 )/∂ xi = 0, ensures the velocities calculated from Eq. (11) are divergence free. 2.4. Local relative stationary method for a moving boundary There are two major issues that require special attention in the numerical implementation of the partial cell method for a moving body. One is tracking the moving body and another is ensuring conservation of mass and momentum [21]. In this study only rigid bodies are considered but the procedures are generalizable. The body surface is expressed by the zero set of a combination of quadratic functions e.g. Fm (x, y, t n ) = 0. The updated position of the body surface on the grid can be determined by simple geometric principles at time step n+1, e.g. Fm (x, y, t n+1 ) = 0. The openness function θ n+1 can also be updated by a combination of quadratic functions. To ensure that the conservation laws are satisfied during body movement from one time step to another the following treatment is adopted. We assume that during the small time interval 1t the body moves with constant velocities uiB , and recast the moving body problem into an equivalent problem where the body remains stationary and the fluid motion is superimposed through an additional velocity with the same magnitude as that of the moving body but in the opposite direction. This is referred to as the local relative stationary method. Details can be found in [21]. Thus, the continuity equation (8) in the partial cell is revised to be,

∂ (θ ui − θ uiB ) ∂ (θ ui ) ∂ (θ uiB ) =0⇒ = . ∂ xi ∂ xi ∂ xi

(13)

(8)

Substituting the above modified continuity equation into the twostep method, a modified Poisson pressure equation is obtained,

(9)

∂ xi

When θ = 1 the above equations return to Eqs. (1) and (2). 2.3. Numerical solution in fluid domain In the spatial domain, a staggered-grid system is used within which scalars are defined at cell centers and vectors are defined at the centers of cell faces, see Fig. 1. The finite difference method is adopted to discretize spatial derivatives on a non-uniform mesh. When a value is required at a location where it is not defined linear interpolation is used to obtain the value. A two-step projection method is used to solve Eqs. (8) and (9) with a combined upwind and central scheme for the convection terms and a central difference scheme to discretize the diffusion terms [24]. After the convection and diffusion terms have been calculated the tentative velocities u˜ i are obtained from

θ u˜ i − θ uni ∂ (θ ui ) 1 ∂τij = −θ uj +θ . 1t ∂ xj ρ (φ) ∂ xj

where the pn+1 is obtained by solving the Poisson pressure equation,

(5)

where H (φ) is a Heaviside-type function, ρw and µw are density and viscosity of water respectively and ρa and µa are density and viscosity of air. The Heaviside function is used to define a smooth function which ensures that the density or viscosity vary smoothly from ρa to ρw or from µa to µw within a narrow band around the interface. There are several choices for the Heaviside function. We adopt the following formula,

(11)

(10)



θ

n+1

∂ pn+1 ρ (φ) ∂ xi 1

1 ∂θ n+1 u˜ i − uiB



 =

1t

∂ xi

 +

∂ (θ gi ) . (14) ∂ xi

It is seen that a moving body changes the pressure field around the body, which consequently affects the flow around it. Apart from the pressure field, the moving body affects the viscous stresses. This is considered in the first step of the projection method where the convection and diffusion terms are evaluated. The local relative stationary method is still applied to calculate the relative velocity gradients near the body. Additional terms result that modify the normal and shear stress around the moving body. The advantages of the local relative stationary method are that it is very easy to implement numerically and easy to deal with any number of multiple bodies since a moving coordinate system is established at the level of an individual cell. No re-meshing is required and the flow computation remains based on a fixed-grid system. The movement of the body is updated in a Lagrangian manner and thus is free of numerical diffusion. Extensions to three-dimensions and more complicated motion than shown here e.g. with rotation are also straightforward (see [21]).

4

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

Fig. 1. Illustration of partial cell technique, level set definition and the storage location of the dependent variables.

2.5. Level set method There are many methods [12,14–16,25–28] to solve the level set Eq. (4), amongst these the semi-Lagrangian method has high computational efficiency. We modify the original formula to be found in the source references for use on non-uniform grids. The distance function in 2D at the time step n + 1 can be expressed as, 1 φin,+ = αβφrn+1,s+1 + (1 − α)βφrn,s+1 + α(1 − β)φrn+1,s j

+ (1 − α)(1 − β)φrn,s

(15)

where r and s relate to the left bottom corner index of the grid cell and the new grid point value comes from that at the time step n according to the velocities ui,j and vi,j . Parameters α and β are the coefficients for the bilinear interpolation in the x and y directions respectively. The formula can easily be extended to 3D. This scheme is unconditionally stable, so the time step is no longer limited by the CFL condition. In our implementation, we only calculate the distance function within a narrow band around the interface, which is normally taken as 3–6 times the local spatial step on each side of the interface. Thus, at each time step the zero level set must not be allowed to pass out of the narrow band. This places a modest but acceptable restriction on the maximum allowable time step. Unfortunately, as noted by Sussman [29], the level set quickly ceases to be a signed distance function especially for flows undergoing extreme topological changes. Re-initialization algorithms can maintain the signed distance property by directly solving the equation

|∇φ| = 1,

(16)

or by solving the steady state equation

φτ + sgn(φ0 )(|∇φ| − 1) = 0,

(17)

where τ is a fictional time, sgn(φ0 ) is a one-dimensional signum function approximated numerically as sgn(φ0 ) = 

φ0 φ + ε2 2 0

,

(18)

and ε is a small value that can be equal to 1x. Using the highly efficient fast marching method [30] we implement a re-initialization process by solving Eq. (16) with the computation also performed

in a narrow band around the interface. The narrow band approach makes a worthwhile saving in computational time. To improve the conservation of the level set method, Enright et al. [12] proposed a hybrid particle level set method in which particles are used to assist the level set method to accurately track the flow characteristics in under-resolved regions. We implement this method to correct the level set values in the level set calculation. During the calculation, a particle reseeding operation [12,15] is used. This operation does not need to be implemented every time step, especially where the interface curvature is not very high. A heap structure is used for sorting the data within the fast marching method and during particle reseeding, which ensure that the interface tracking procedure remains very efficient. 2.6. Level set boundary condition for obstacle Within the numerical model, the level set function at the boundary of a fixed or moving body must be carefully treated. Several researchers have described treatments for this problem by considering the angle between the obstacle surface normal direction and water surface normal direction [23] or by interpolation from a nearby cell (Guendelman et al. [22]). In our model the angle between the water surface and the obstacle surface can remain arbitrary consistent with the manner in which the water surface comes into contact with the obstacle surface. By assuming a constant slope for the free surface near the solid surface, the level set function values at the solid boundary can be approximated. Using the formula for calculating the distance between a point and line, the level set function of a point m with coordinate (xm , ym ) which is the closest point to a boundary point n can be expressed as

ϕm = axm + bym + c

(19)

see Fig. 2 where point n is a cell center point where the level set function is defined, and the cell is crossed by a body surface. Around the cell only the cell with cell center point m is fully occupied by fluid and is closest to the free surface. If ϕm is its level set value and c is a constant, (a, b) is the normal unit vector which can be calculated from ∇ϕ/ |∇ϕ|, the boundary level set value of point n with coordinate (xn , yn ) can then be obtained by

ϕn = axn + byn + c .

(20)

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

5

Fig. 2. Illustration of body boundary level set condition.

Fig. 3. Illustration of the narrow band level set approach involving moving solid boundaries. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Eqs. (19) and (20) can be easily extended into 3D case. The boundary condition for level set is also confined within the narrow band. If the body is moving, the boundary level set value outside the narrow band must also be considered. Outside the narrow band we set the level set value to +106 in water area and −106 in air area. If the moving body is initially located only in air, when it moves into the water region, the moving boundary condition for the level set in the narrow band is calculated by (17) and (18); outside the narrow band we set the moving boundary level set value to that of the cell adjacent to it in the fluid area. This method keeps

the level set value within the body at +106 when it moves into the water area and allows the body to pass through the interface smoothly. If the moving body is initially only in the water region, when it moves into the air region, the same method is employed. To illustrate the moving body boundary treatment we show a case of a circular cylinder descending into initially calm water at constant velocity. Fig. 3 shows the layout of the local mesh, the locations of the cylinder, the free surface, the positive and negative particles and the narrow band during the simulation. From the figure we can see how the circular cylinder passes cleanly through the water free surface.

6

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

Fig. 4. Initial status of water column breaking up.

Fig. 5. Particles distribution with or without reseeding at t = 0.2 s.

Fig. 6. Comparison of water column height at left wall with experimental data. Fig. 7. Comparison of water front location with experimental data.

3. Validation 3.1. Collapse of 2D water column In free surface hydrodynamics, the sudden collapse of a rectangular column of fluid onto a horizontal surface is a classical problem. The dam-break problem is a popular validation case for various surface tracking and capturing methods, because of the relatively simple geometry and initial conditions. Martin and

Moyce [31] performed a benchmark laboratory experiment using a high-speed camera to capture the fluid motion at selected time intervals and obtain data for the free surface position at the walls as functions of time. Many investigators have repeated the case experimentally or numerically [8,14,32,33]. In this test the water column height and width and the tank length is set according to the experiment as shown in Fig. 4. The domain is 0.6 m × 0.6 m, the height of the water column is 0.3 m and its width is 0.15 m. The

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

7

Fig. 8. Water column break up (mesh 80 × 80).

density of water is 1000 kg/m3 and for air 1.2 kg/m3 , and the acceleration due to gravity is taken to be 9.81 m/s. The viscous terms in the model are set to zero. Initially the velocity and pressure is zero everywhere. All boundaries of the tank are set as free slip boundaries. Firstly the effect of particle reseeding based on a mesh with 40 × 40 cells is tested. Fig. 5 shows the distribution of the particles around the free surface at t = 0.2 s for cases with or without reseeding. If particles are transported without reseeding negative particles originally in the top-right corner of the water column will remain far away from the free water surface, and cannot correct the local level set within the narrow band. If every 1, 5 or 10 time steps the particles are reseeded, they all remain closely congregated around the free surface. This shows that particle reseeding is important if particles are to be used to correct the local level set function.

A quantitative comparison with the experiment is made by using the reduction in water column height at the left wall and position of the leading-edge of the water front at the tank bottom. The non-dimensional height of the collapsing water column at the left wall versus the non-dimensional time is shown in Fig. 6. The numerical results for three successively refined meshes of sizes 20 × 20, 40 × 40 and 80 × 80 are presented and in all calculations the predicted water column heights at the left wall are in close agreement with the experimental data. The nondimensional positions of the water front for the same test case are shown in Fig. 7. The middle and the finest meshes give good results compared with the experimental data. This shows that the present model can use much coarser meshes and provide results of comparable accuracy to those methods in other published work [14,32,34].

8

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

Fig. 9. Comparison of interface pattern with particles and without particles.

The free surface patterns during the flow development are also compared with other numerical results [20,32] for the mesh with 80 × 80 cells. In this case the initial time step is set to 0.001 s and is adjusted automatically with the Courant number set to 0.3 in the calculations. For a simulation of 1.2 s of real physical time the CPU time is about 54 min (PC Intel(R) Core(TM) 2 Duo CPU E6750 @2.66 GHz 2.66 GHz, 3.24 GB of RAM). In the calculations, particles are reseeded at every time step. Fig. 8 shows snapshots of the free surface at times t = 0.2, 0.6, 0.7, 0.8, 0.9, 1.2 s. The water moves to the right from the initial state under gravity and at t = 0.2 s approximately 75% of the tank base is covered with water. The water front reaches the right wall at about t = 0.3 s and then the water moves up the right wall under its own inertia. When its potential energy reaches a maximum, the water column at the right wall starts falling downwards. A snapshot at t = 0.6 s confirms this instant. For t = 0.7 s to t = 1.05 s the returning free

surface has turned over and some air has been trapped in the body of water at around t = 1.0 s. Afterwards the backward moving water hits the left wall and moves upwards. Around t = 1.2 s another maximum in the potential energy is reached and the water column along the left wall attains its highest position. For this test case, if no particles are used in the level set solution procedure, some of the detail of air entrapment in the water at around t = 0.8 s and t = 1.0 s would be lost. Fig. 9 shows a comparison of the results from simulations with particles and without particles. The mesh of 80 × 80 is used. Without particles a clear pattern can be obtained which is similar to results with particles, but air entrapment cannot be observed. With particles, which improve the conservation of the level set function, the interface can be better captured so that more intricate details of the interface connection, separation and air entrapment are observed. This shows that particles can enhance the numerical

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

9

Fig. 10. Water free surface pattern: numerical (left) and experimental (right) results at various times.

accuracy but at the cost of longer computing times. For example, without the use of particles about 20 min is required to run the test case for a real time of 1.2 s (PC Intel(R) Core(TM) 2 Duo CPU E6750 @2.66 GHz 2.66 GHz, 3.24 GB of RAM), but

54 min are required if particles are used. So in practice, a balance should be struck between computational costs and numerical accuracy when deciding if particles should be deployed in the calculation.

10

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

Fig. 11. Free surface pattern (from top-left to bottom-right t = 0, 10, 15, 30 s).

Fig. 12a. Water elevation on surface of inner cylinder.

3.2. Water column collapsing and impinging upon a rectangular obstacle In this test case, a water column collapses and impinges upon a rectangular obstacle. A water column 0.3 m high and 0.15 m wide is initially confined to the left of a tank. The tank is 0.6 m long in the horizontal direction and a rectangular obstacle 0.024 m long and 0.048 m high is set in the middle of the tank. At time t = 0 the confinement is suddenly withdrawn and the water moves to the right and strikes the object. Subsequently, a tongue of water moves up and over the top surface of the object, and impinges on the back wall. In this process the

complex phenomena of air entrapment, wave breaking and mixing are observed. Koshizuka [33] reported a laboratory experiment with photographs of the free surface pattern. In the present 2D numerical tests, the domain is 0.6 m × 0.45 m in the horizontal and vertical directions and the corresponding mesh has 100 × 75 cells. The real physical elapsed time for the events described is 1.0 staking about 3.5 h running on a PC Pentium(R) 4 CPU 2.80 GHz 2.79 GHz, 2.00 GB of RAM. Fig. 10 shows a comparison between the water free surface at different times and the corresponding free surface patterns from the experiments. Fig. 10(a) shows the initial conditions at t = 0 s. Fig. 10(b) shows a tongue of water moving over the rectangular box at t = 0.2 s. Fig. 10(c) captures

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

11

Fig. 12b. Horizontal force acting on inner cylinder.

the instant where the tongue has impinged upon the back wall, with a large block of air entrapped behind the box and under the tongue of water. Fig. 10(d) shows that the advancing water front has separated into two parts at t = 0.5 s with one part moving up the back wall and another part moving down and forwards along the bottom right corner. The body of water connected to the top right part of the box then moves downwards. Fig. 10(e) shows the flow pattern at t = 1.0 s at which time several blocks of air are enclosed behind the object and a body of water is again moving backwards towards it. The predicted positions and pattern of the free surface compare very well with Koshizuka’s results. 3.3. An impulse wave interaction with a bottom-mounted circular cylinder An impulse wave in a circular cylindrical tank interaction with an inner circular cylinder placed at its center is next simulated. The combination of incident, reflected and diffracted waves leads to rather complex free surface behavior. Bai et al. [35] have studied this case using a higher-order boundary element method and Chern et al. [36] have studied it using a pseudo-spectral matrix element method. We compare our results with those of Bai et al. simulating the flow interaction with the obstacle by calculating the force and free surface run-up on the cylinder surface. This case is a fully nonlinear free surface problem. The initial wave profile is defined as

   2 r + rd 2 −y η(t = 0) = a exp − x − 2

where r and rd are the radii of the inner cylinder and the outer cylindrical tank, which are set to 1 m and 10 m respectively, and a is the amplitude of the initial wave profile, set equal to 0.1 m. The impulse wave peak amplitude is initially located half way along the radial line between the walls of the inner cylinder and outer cylinder. The water depth is 1.0 m in this simulation. Both the boundary of the cylindrical tank and the inner circular cylinder are described by the partial cell method in the numerical simulations. Three meshes are used to discretize the problem with 56 × 56 × 37, 106 × 106 × 37 and 206 × 206 × 37 cells uniform in the horizontal direction whilst non-uniform in the vertical direction with 15 cells of size 0.02 m covering the expected range of free surface variation, even in the case of violent deformation. Fig. 11 are plots of the free surface elevation from t = 0 s to t = 30 s, showing the process of impulsive wave propagation. In the initial

impulse the wave starts half way between two cylinders. The body of water falls under gravity and propagates in every direction reaching the inner cylinder from the right and the outer cylinder from the left (at around t = 10 s). Reflected waves are generated at both cylinders and diffracted waves at the inner cylinder (t = 15 s). As time progresses, the incident wave, reflected waves and diffracted waves are superposed and complex wave profiles are formed in the tank (t = 30 s). The free surface patterns are similar to those predicted by Bai et al. [35] and Chern et al. [36]. A numerical wave gauge is positioned at the front of the inner cylinder to record the wave elevation and the horizontal force on the inner cylinder is calculated by interpolation. Figs. 12a and 12b show the numerical results on different meshes compared with the results of Bai et al. [35]. It can be seen that the results from the middle mesh are close to those on the fine mesh, confirming mesh convergence and are also in satisfactory agreement with the results of Bai et al.

3.4. 3D dambreak simulation

A classical dambreak and subsequent action on a box is used to test the capability of our model to simulate a complex 3D free surface flow interaction with an obstacle. The experiment is conducted in a large tank of dimensions 3.22 m × 1 m × 1 m with an open roof [37]. A water column 1.22 m wide and 0.55 m high is enclosed in the right part of the tank by a door. Positioned in front of the door is a box of width 0.4 m and height 0.16 m centered at (0.75, 0) on the floor of the tank. When the door is opened impulsively water escapes confinement into the left part of the tank and a complex interaction between the water, box and left vertical wall is observed. This case is a 3D extension of the dam break simulation discussed in 3.2. In the reported laboratory experiments [37], water height and pressure were measured at the positions shown in Fig. 13. Four water height sensors were located at H1 : (x, y, z ) = (0.5, 0.0, 0.0 − 1.0); H2 : (x, y, z ) = (1.0, 0.0, 0.0 − 1.0); H3 : (x, y, z ) = (1.5, 0.0, 0.0 − 1.0); H4 : (x, y, z ) = (2.66, 0.0, 0.0 − 1.0). Eight pressure sensors were positioned at the front and top of the box,

12

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

Fig. 13. Dambreak setup and measurement positions. Source: From Kleefsman et al. [37].

e.g. P1 : (x, y, z ) = (0.831, −0.026, 0.025); P3 : (x, y, z ) = (0.831, −0.026, 0.099); P5 : (x, y, z ) = (0.806, 0.026, 0.165); P7 : (x, y, z ) = (0.733, 0.026, 0.165). In the numerical simulations, two meshes were used with 80 × 25 × 25 and 160 × 50 × 50 cells in the x, y and z directions, and uniform spacing. The time to 6 s was simulated with the Courant number set to 0.3. In the calculations, a coarse mesh calculation took about 30 h, and fine mesh calculation about 9 days on a Pentium(R) 4 CPU 2.80 GHz PC with 2.00 GB of RAM. Fig. 14 shows the water free surface elevation at t = 0.0 s, 0.5 s, 1.0 s and 2.0 s respectively. The complicated topology of the water free surface impacting on the cube, separating, breaking and mixing can be clearly observed. The time history of water height at the designated four locations is compared with experimental data in Fig. 15, and the pressure time history in Fig. 16. It is shown that the predictions of water height and pressure on the two meshes are in close agreement with the experimental observations. This test confirms that with an appropriate mesh the present numerical model can be used to accurately simulate complicated three dimensional free surface flows and interactions with an obstacle. 3.5. Numerical wave generation in a wave tank A piston paddle oscillating horizontally is often used for generating waves in a physical wave tank. This is an appropriate case to use for testing the ability of the present model to generate waves in a tank numerically using a numerical wave paddle. In this test, the laboratory experiment conducted by Gao [38] is reproduced numerically in 2D. A piston paddle is located at the left hand end of a tank 8.85 m long with a water depth of 0.28 m. Fig. 17 shows the piston paddle displacement, from which the velocity of the paddle can be obtained for use as a boundary condition in the numerical simulation. Three wave gauges are set at distances of 0.55, 3.55 and 5.45 m from the paddle to record the wave history, the latter gauge placed roughly three quarters of the distance down the tank. There is 10 s of recorded wave history data available from starting the paddle moving. The setup of the numerical wave tank is the same as in the experiment. To test grid convergence three uniform meshes were used with 1x = 1z = 0.01 m, 1x = 1z = 0.014 m and 1x = 1z = 0.02 m, respectively. Fig. 18 shows the numerical results from the three meshes compared with the experimental data. Unsurprisingly, the coarsest mesh under-predicts the wave peaks and wave troughs while the values predicted by the two finer meshes are in closer agreement. The wave peaks generated in the numerical wave tank are very close to the experimental data but the wave troughs are slightly different. The results do not show any loss of wave peak amplitude along the tank as a result of numerical dissipation, as observed with other numerical methods. In general terms, the fine mesh results are in satisfactory agreement with the laboratory data. The present results indicates that the numerical wave tank is suitable for the study of wave interactions with fixed or floating bodies including cases with breaking waves, which will be reported separately.

Fig. 14. Water free surface at different times.

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

Fig. 15. Water height time history at four gauge sensors.

13

Fig. 16. Pressure time history at four pressure sensors.

3.6. Motion of two elliptical cylinders near a free surface This case demonstrates the capabilities of the numerical model to simulate an intricate free surface flow interaction involving two moving elliptical-shaped bodies. The present test case with two ellipses is based on a case with one ellipse that moves through still water in a tank reported by Lin [21]. In the present case, as one ellipse moves horizontally right-to-left below the surface through still water, a second ellipse is set to move down and up in the middle of the tank. The initial position of the two bodies is shown in Fig. 19. The tank is 10 m long and 3 m high with water depth 2.5 m. The center of the right hand ellipse is initially set at (8 m, −0.2 m) and its surface is defined by (x − 8)2 /0.22 + (z + 0.2)2 /0.12 = 1. The body is set moving to the left from stationary going through a periodic acceleration and deceleration with ax = − sin(ω1 t ) where ω1 is the angular frequency and period 6.0 s. The body also moves in the vertical direction with velocity v1 = 0.16 cos(ω1 t ). The second ellipse moves only in the vertical direction with velocity v2 = 0.48 cos(ω2 t ), and period ω2

Fig. 17. Paddle displacement used in the wave tank (Gao [38]).

of 12.0 s. Its center is initially set at (4.95 m, 0.2 m) and its surface defined by (x − 4.95)2 /0.22 + (z − 0.2)2 /0.12 = 1. In the simulation a mesh with 400 × 67 cells is used to discretize the domain. A uniform spacing of 0.025 m is used in the horizontal

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

eta (m)

0.04 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 0.04 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04

eta (m)

eta (m)

14

0.04 0.03 0.02 0.01 0 -0.01 -0.02 -0.03 -0.04

experiments dx =0.01m dx =0.014m dx =0.02m

0

2

4

6

8

10

times (s) Fig. 18. Free surface elevation comparison between numerical and experimental results (from top to bottom wave gauges at 0.55, 3.55 and 5.45 m).

direction. The same uniform spacing of 0.025 m is used from −0.5 to 0.5 m in the vertical direction around the free surface, and then gradually increased with distance from −0.5 to −2.5. Before simulating two bodies moving at the same time, each body moving separately is first simulated as shown in Figs. 20 and 21. Then the two elliptical cylinders moving in accordance with the described motions is simulated. Vorticity induced by the body motion is shown in Figs. 20–22 with free surface elevations at different times. Fig. 20 with single body motion corresponds to the test case from Lin [21] based on solving single phase flow with a VOF method to track the motion of the free surface. The present results are in close agreement with those of Lin. The body induced vortices interact with the free surface in a very complicated manner due to the relatively long period of body motion. The body motion also induces waves in front and behind it that lie close to the body when the body speed is relatively large (t < 3 s). Wave steepness is very large in front and behind the body and the strongest vorticity occurs at t = 3.0 s when the body reaches its highest speed (Fig. 20(d)). For t > 3 s the body is gradually slowing down and the wave steepness gradually decreases (Fig. 20(e)). Before the body stops the vertical speed is increased inducing a small wave with higher wave steepness at the free surface (Fig. 20(f)). As the body slows to a stop the waves on the free surface gradually decrease in amplitude and the vorticity induced by the body vanishes (Fig. 20(g)). Fig. 21 shows the vertical motion of the second elliptical cylinder in the tank. At the beginning the elliptical cylinder pushes the water free surface down and induces the adjacent water free surface to move upwards (Fig. 21(b)). Waves are generated above the cylinder and transported to the left and right of the tank. After the body moves into the water vortices are generated along with the body and behind it as it continues its descent (Fig. 21(c) and (d)). At t = 3.0 s the body reaches its lowest position and then starts moving upwards. For

z 0.5 0

x -2.5 0

10 Fig. 19. Setup of two elliptical cylinders in a tank.

t > 3.0 s vortices are again induced by the body behind it as it moves upwards (Fig. 21(e) and (f)), whilst waves move outwards horizontally. Finally, the body stops at its initial position and the water moves out from it bilaterally downwards (Fig. 21(g)). In the extended case where two elliptical cylinders move simultaneously in the tank, the interaction between the respective bodies and the water free surface is even more complicated (Fig. 22). This test is used to show that the present numerical model has the ability to treat complicated free surface flow phenomena. Fig. 22 shows that the resulting flow phenomena, as expected, is a superposition of the two separate cases with the single elliptical cylinders. At the beginning the separation distance between the two cylinders is relatively far, so their behavior is largely uncoupled (Fig. 22(b)). When the right cylinder moves towards the left the waves generated by the two bodies are superposed. The free surface and the vorticity induced by two bodies is obviously more complicated than the individual cases, and Fig. 23 shows the local details of the flow at t = 3.0 s. From these tests we conclude that the capability of the present method to treat complex free surface flow interactions with moving bodies appears promising.

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

15

Fig. 20. An elliptical cylinder moving from right to left, colors represent vorticity (s−1 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 21. An elliptical cylinder moving down and up, colors represent vorticity (s−1 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

4. Conclusions

very effective for conserving volume and mass when the interface is in motion. The method can track a moving interface undergoing complex topological changes including shifting, stretching, tearing, breaking and merging. Several benchmark test cases involving complicated free surface flow interactions with fixed or moving bodies have been used to validate the new numerical model and generally good agreement has been achieved between the numerical results and analytical solutions or experimental data. The combined flow solver is being extended to cases involving wave interactions with floating bodies in two and three dimensions and may be linked to models that include compressibility within the flow calculations to accommodate aeration within a composite flow model.

A free surface solver has been developed for the simulation of wave interactions with fixed or moving bodies in which a narrow band Semi-Lagrangian fast marching particle level set method has been incorporated into an efficient finite difference incompressible Navier–Stokes flow solver in conjunction with the partial cell and local relatively stationary methods to simulate wave interaction with fixed and moving bodies. A level set advective solver is implemented within a narrow band around the free surface and in principle the Semi-Lagrangian method has no formal time step restriction due to any stability constraint. The fast marching technique is very efficient and when used with particles it is found to be

16

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17

References

Fig. 22. Two elliptical cylinders moving in the tank, colors represent vorticity (s−1 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 23. Local vorticity at t = 3.0 s.

[1] C.M. Lemos, A simple numerical technique for turbulent flows with free surfaces, Internat. J. Numer. Methods Fluids 15 (1992) 127–146. [2] J. Glimm, J. Grove, B. Lingquist, O.A. McBryan, G. Tryggvason, The bifurcation of tracked scalar waves, SIAM J. Sci. Stat. Comput. 9 (1) (1988) 61–79. [3] F.J. Kelecy, R.H. Pletcher, The development of a free surface capturing approach for multidimensional free surface flows in closed containers, J. Comput. Phys. 138 (1997) 939–980. [4] L. Qian, D.M. Causon, C.G. Mingham, D.M. Ingram, A free-surface capturing method for two fluid flows with moving bodies, Proc. R. Soc. A 462 (2006) 21–42. [5] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [6] S. Osher, R. Fedkiw, Level set methods: an overview and some recent results, J. Comput. Phys. 169 (2) (2001) 463–502. [7] S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, in: Applied Mathematical Sciences, vol. 153, Springe-Verlag, New York, 2002. [8] O. Ubbink, Numerical prediction of two fluid systems with sharp interfaces, Thesis, Department of Mechanical Engineering, Imperial College of Science, Technology & Medicine, January 1997. [9] F. Harlow, J. Welch, Numerical calculation of time-dependent viscous incompressible flow, Phys. Fluids 8 (1965) 2182–2189. [10] A. Harten, S. Osher, B. Engquist, S.R. Chakravarthy, Some results on uniformly high order accurate essentially non-oscillatory schemes, Appl. Numer. Math. 2 (1986) 347–377. [11] X.D. Liu, S. Osher, T.F. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200–212. [12] D. Enright, R. Fedkiw, J. Ferziger, I. Mitchell, A hybrid particle level set method for improved interface capturing, J. Comput. Phys. 183 (2002) 83–116. [13] D. Enright, D. Nguyen, F. Gibou, R. Fedkiw, Using the particle level set method and second order accurate pressure boundary condition for free surface flows, in: Proceedings of FEDSM’03 2004, 4th ASME JSME Joint fluids Engineering Conference, Honolulu, Hawaii, USA, July 6–11, 2003. [14] J. Strain, Semi-Lagrangian methods for level set equations, J. Comput. Phys. 151 (2) (1999) 498–533. [15] D. Enright, F. Losasso, R. Fedkiw, A fast and accurate semi-Lagrangian particle level set method, Comput. Struct. 83 (2005) 479–490. [16] G.S. Jiang, D.P. Peng, Weighted ENO schemes for Hamilton–Jacobi equations, SIAM J. Sci. Comput. 21 (6) (2000) 2126–2143. [17] J. Sethian, Level Set Methods, Cambridge University Press, Cambridge, UK, 1996. [18] J. Sethian, Level Set Methods and Fast Marching Methods, Cambridge University Press, Cambridge, United Kingdom, 1999. [19] J. Sethian, Evolution, implementation, and application of level set and fast marching methods for advancing fronts, J. Comput. Phys. 169 (2) (2001) 503–555. [20] L. Qian, D.M. Causon, D.M. Ingram, C.G. Mingham, A Cartesian cut cell twofluid solver for hydraulic flow problems, ASCE J. Hydraul. Eng. 129 (9) (2003) 688–696. [21] P. Lin, A fixed-grid model for simulation of a moving body in free surface flows, Comput. & Fluids 36 (2007) 549–561. [22] E. Guendelman, A. Selle, F. Losasso, R. Fedkiw, Coupling water and smoke to thin deformable and rigid shells, in: Proceedings of ACM SIGGRAPH 2005 TOG Homepage, ACM Trans. Graph. (TOG) 24 (3) (2005) 973–981. [23] J. Yang, F. Stern, Sharp interface immersed-boundary/level set method for wave-body interactions, J. Comput. Phys. 228 (2009) 6590–6616. [24] P. Lin, Numerical modeling of breaking waves, A Dissertation Presented to the Faculty of the Graduate School of Cornell University, 1998. [25] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 83 (1) (1989) 32–78. [26] G.S. Jiang, C.C. Wu, A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics, J. Comput. Phys. 150 (2) (1999) 561–594. [27] J. Strain, A fast modular semi-Lagrangian method for moving interfaces, J. Comput. Phys. 161 (2000) 512–536. [28] W. Mulder, S. Osher, J. Sethian, Computing interface motion in compressible gas dynamics, J. Comput. Phys. 100 (2) (1992) 209–228. [29] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible two-phase flow, J. Comput. Phys. 114 (1) (1994) 146–159. [30] J. Sethian, Fast marching methods, Dept. of Mathematics, University of California, Berkeley, 1998. [31] J.C. Martin, W.J. Moyce, An experimental study of the collapse of liquid column on a rigid horizontal plane, Phil. Trans. Roy. Soc. Lond. Ser. A 244 (1952) 312–324. [32] D. Greaves, Simulation of interface and free surface flows in a viscous fluid using adapting quadtree grids, Internat. J. Numer. Methods Fluids 44 (2004) 1093–1117. [33] S. Koshizuka, H. Tamako, Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation, Int. J. Comput. Fluid Dyn. 4 (1) (1995) 29–46. [34] X. Lv, Q.P. Zou, D. Reeve, Z.Y. Wang, An unstructured 3D LES solver for free surface flow and breaking waves, in: Proceedings of the Eighteenth (2008) International Offshore and Polar Engineering Conference, Vancouver, BC, Canada, July 6–11, 2008, pp. 62–68. [35] W. Bai, R. Eatock Taylor, Higher-order boundary element simulation of fully nonlinear wave radiation by oscillating vertical cylinders, Appl. Ocean Res. 28 (2006) 247–265.

H.B. Gu et al. / European Journal of Mechanics B/Fluids 38 (2013) 1–17 [36] M.J. Chern, A.G.L. Borthwick, R. Eatock Taylor, Simulation of non-linear free surface motions in a cylindrical domain using Chebyshev–Fourier spectral collocation method, Internat. J. Numer. Methods Fluids 36 (2001) 465–496.

17

[37] K.M.T. Kleefsman, G. Fekken, A.E.P. Veldman, A volume-of-fluid based simulation method for wave impact problems, J. Comput. Phys. 206 (2005) 363–393. [38] F. Gao, An efficient finite element technique for free surface flow, Ph.D. Thesis, Brighton University, UK, 2003.