Search or ask a question
Journal ArticleDOI

# A Volume-of-Fluid based simulation method for wave impact problems

10 Jun 2005-Journal of Computational Physics (ACADEMIC PRESS INC ELSEVIER SCIENCE)-Vol. 206, Iss: 1, pp 363-393
TL;DR: In this paper, the authors considered some aspects of water impact and green water loading by numerically investigating a dambreak problem and water entry problems, based on the Navier-Stokes equations that describe the flow of a viscous fluid.
About: This article is published in Journal of Computational Physics.The article was published on 2005-06-10 and is currently open access. It has received 618 citations till now. The article focuses on the topics: Finite volume method & Volume of fluid method.

## Summary (5 min read)

### 1. Introduction

• During the last decades increasing attention has been paid to the development of calculation methods for wave impact problems.
• In high or steep waves offshore structures face the problem of wave impact on their bow or bottom.
• A good prediction method is based on the Navier–Stokes equations, which describe the motion of an incompressible, viscous fluid, with a free liquid surface.
• The physical properties such as surface tension can be applied easily.

### 2. Governing equations

• Flow of a homogeneous, incompressible, viscous fluid is described by the continuity equation and the Navier–Stokes equations, describing conservation of mass and momentum, respectively.
• Further, l denotes the dynamic viscosity and F = (Fx,Fy,Fz) is an external body force, for example gravity.

### 2.1. Boundary conditions

• At the solid walls of the computational domain and at the objects inside the domain a no-slip boundary condition is used.
• Some of the domain boundaries may let fluid flow in or out of the domain.
• In their method the wave on the inflow boundary can be prescribed as a regular linear wave or a regular 5th order Stokes wave [39].
• At the outflow boundary a non-reflecting boundary condition is needed to prevent the waves from reflecting from the boundary into the domain.
• A Sommerfeld condition (see, e.g. [40]) is very appropriate in cases where a regular wave is used.

### 2.2. Free surface

• At the free surface boundary conditions are necessary for the pressure and the velocities.
• Here, un and ut are the normal and tangential component of the velocity, respectively, p0 is the atmospheric pressure, r is the surface tension and j denotes the total curvature of the free surface.

### 3. Numerical model

• To solve the Navier–Stokes equations numerically the computational domain is covered with a fixed Cartesian grid.
• The variables are staggered as in the original MAC method [20], which means that the velocities are defined on cell faces, whereas the pressure is defined in cell centres.
• The body geometry is piecewise linear and cuts through the fixed rectangular grid.
• Note that these cells do not have to be completely filled with fluid.
• The velocities that are positioned at the cell faces are labelled using the cell labels.

### 3.1. Discretisation of the continuity equation

• The continuity and momentum equations are discretised using the finite volume method.
• The natural form of the equations when using the finite volume method is the conservative formulation as given in Eqs. (1) and (2).
• Axe AxwÞ; dxðAyn AysÞÞ=kdyðAxe AxwÞ; dxðAyn AysÞk: ð8Þ Recognising l in the denominator of nb, the discrete continuity equation (Eq. (7)) can be written as ueA x edy þ vnAyndx uwAxwdy vsAysdxþ ubðAxe AxwÞdy þ vbðAyn AysÞdx ¼ 0: ð9Þ.

### 3.2. Discretisation of the Navier–Stokes equations

• In the case of uncut cells this simply means that the control volume consists of half of the neighbouring cells to the left and to the right of the velocity (see the left of Fig. 3).
• The discretisation is explained for the momentum equation in x-direction.
• Here, uc is the central velocity around which the control volume is placed and F b cdxcdy is the volume of the control volume (see Fig. 3 for the notation).
• The part of the fluxes due to the moving body through the upper and lower cell boundaries contains a max-function that distinguishes between the situation that the body is moving into or out of the cell.
• This method has been chosen to prevent instabilities, which can occur in the discretisation when a division is needed by the volume of a cut cell that can be arbitrary small.

### 3.3. Temporal discretisation

• The equations of motion are discretised in time using the forward Euler method.
• This first order method is accurate enough, because the order of the overall accuracy is already determined by the first order accuracy of the free-surface displacement algorithm (see Section 3.5).
• Using superscript n for the time level, the temporal discretisation results in: M0unþ1h ¼ Mbunþ1b ; ð24Þ X unþ1h unh dt ¼ Cðunh; ubÞunh 1 q ðM0ÞTpnþ1h lDunh þ Fnh: ð25Þ.
• The continuity equation is discretised at a new time level to ensure a divergence free velocity field.
• The spatial discretisation is written in matrix notation, where M is the divergence operator with M0 working on the interior velocities and Mb on the boundary velocities, X contains cell volumes, C contains the convective coefficients (which depend on the velocity vector) and D contains diffusive coefficients.

### 3.4. Solution method

• First an auxiliary vector field ~unh is calculated from Eq. (27).
• From this equation, which does not require boundary conditions, the pressure is solved using the SOR (Successive Over Relaxation) method where the optimal relaxation parameter is determined during the iterations [2].
• Once the pressure field is known, the new velocity field is calculated from Eq. (26) using the pressure gradient.

### 3.5. Free-surface displacement

• After the new velocity field has been calculated, the free surface can be displaced.
• This is done using an adapted version of the VOF method first introduced in [23].
• In the neighbourhood of the free surface, the calculation of the fluxes is somewhat more complicated.
• The initial condition has also been plot.
• When combining VOF with a local height function, the amount of flotsam and jetsam has decreased considerably as can be seen in the right of Fig.

### 3.6.1. Pressure at the free surface

• A boundary condition for the pressure is needed in surface cells.
• The term containing the viscosity is neglected, which leaves p ¼ p0 rj: ð31Þ.
• To calculate the contribution of the surface tension in the pressure at the free surface, the total curvature of the free surface has to be determined in every S-cell.
• The local height function is defined based on the orientation of the free surface.
• If, for example, the orientation of the free surface is mainly horizontal (as in Fig. 11) with the fluid below the free surface, the fluid cell below the surface cell is used for the interpolation.

### 3.6.2. Velocities at the free surface

• Velocities in the neighbourhood of the free surface can be grouped in different classes (see Fig. 12).
• The first class contains the velocities between two F-cells, between two S-cells and between an S- and F-cell.
• These velocities are determined by solving the momentum equation, so they are called momentum velocities.
• These velocities are determined using boundary conditions, which will be described below.
• These are determined using the tangential free-surface condition, Eq. (6), as described in [12].

### 3.6.3. SE-velocities

• The choice for the method to determine the velocities at the cell faces between surface and empty cells (SE-velocities) turns out very important for the robustness and the accuracy of the flow model.
• In the case that more E-cells surround the S-cell, the net mass flux through FS-boundaries is divided over the SE-boundaries such that $Æ u = 0 is satisfied. • The SE-velocities are determined by extrapolating interior velocities, also known as Method 2. • The first column of Table 1 indicates that Method 1 gives less accurate results in wave simulations than Method 2. • When using constant extrapolation, the method is much more robust. ### 4. Stability of the method • In the case of uncut cells with fixed objects the stability of the equation containing the time integration term and the convective term is given by the CFL-restriction, which in one dimension reads dtjuj/h 6 1 (h is the size of the uncut cell). • When cut cells are present, for the chosen convective discretisation this criterion is not changed. • This result is not directly straightforward when looking at the equation containing the time derivative and the convective term ou ot ¼ X 1Cðu; ubÞu: ð37Þ. • When the object is moving tangential to its boundary, the eigenvalues of the matrix C(u,ub) can again be estimated by O(Xu/h), which means that stability is guaranteed when the CFL-restriction is used. • They can become arbitrary large due to the factor X 1. ### 5. Convergence of the method • To check the numerical convergence of the method, which is first order in space and time, simulations have been performed of a regular wave. • The waves are generated using 5th order Stokes theory. • The number of time steps has been varied for a convergence study in time, more specifically, 63, 125 and 250 time steps per period have been used. • The upper part of Fig. 16 shows the wave elevation at the same moment in time as when the error is calculated. ### 6. Interactive body–fluid motion • In the current method the objects are either moving with a prescribed motion, or the motion is calculated from the interaction between the object and the liquid dynamics. • The acceleration of the body can be computed by Newton s second law Abody ¼ I 1F body; ð40Þ where I is the inertia matrix and Fbody contains the forces and moments exerted by the fluid. • Inserting this into the explicit integration gives Anþ1 ¼ ma m An: ð41Þ. • For problems with moving bodies that are heavy, the coupling will be stable with just a bit under-relaxation resulting in faster convergence. • This slows down the method, because in every sub-iteration step the Poisson equation for the pressure has to be solved. ### 7. Dambreak simulation • At the Maritime Research Institute Netherlands experiments have been performed for breaking dam flows. • Four vertical height probes have been used; one in the reservoir and the other three in the tank. • In Fig. 19 time histories of the water height at two locations are shown: in the reservoir, and in the tank just in front of the box. • In the bottom graphs of Fig. 20, where the time history of pressure transducers at the top of the box are shown, a clear difference occurs between simulation and experiment. • The magnitude of the impact is better predicted on the finer grids. ### 8. Water entry of wedge, cone and circular cylinder • Circular cylinders and of a three-dimensional cone. • The tests have been performed with prescribed constant entry velocities. ### 8.2. Cone entry • As the three-dimensional equivalent of the wedge, the entry of a cone has also been studied. • The slamming coefficient, as a non-dimensional measure of the total hydrodynamic force, has been compared with the theoretical prediction by Schiffman and Spencer made in 1951 [38]. • The non-dimensional parameter k(b) is considered to give most accurate results at value 1.6 (see [1]). • Three different entry velocities have been chosen, leading to the same results in a suitable set of scaled variables. • From this it can be seen that the jets at the side of the cone, which are clearly present in the entry of a wedge (see Fig. 22), are not well resolved. ### 8.3. Entry of a circular cylinder • The entry of a circular cylinder has also been studied. • The free-surface shape observed in the experiment is very well resolved by COMFLOW. • Not all the details of the droplets of the splash are captured by the simulation, but the jets that appear at the sides of the cylinder are well predicted. • In Fig. 26 the slamming coefficients of the cylinder entry with different entry velocities have been plotted versus the non-dimensional penetration depth. • Besides the experimental result of Campbell and Weynberg, also the theory of Von Karman (1929), reported by Faltinsen in [9], has been included. ### 9. Free falling wedge • Finally, results are shown of a wedge falling freely into initially calm water. • The measurements have also been used by Muzaferija et al. [32] to compare their method with. • They use a finite volume method with polyhedral control volumes incorporating a freesurface capturing model. • Muzaferija et al. presented some three-dimensional simulations, which also showed a decrease in the vertical force and better agreement with the experiment. • An initial velocity of 6.15 m/s is prescribed after which the wedge is falling freely. ### 10. Conclusions • In the present paper results are shown of the simulation of engineering problems using a Navier–Stokes solver with a VOF-based free-surface displacement. • The combination of the original VOF method with a local height function has shown to improve the free-surface treatment by examining the rotation of a slotted disk. • Finally, a free falling wedge has been simulated, where a full coupling between the liquid motion and the body motion has been established. • The results of the simulations give much confidence in the performance of the method. • For validation, simulations of a moving vessel in high waves resulting in green water on the deck will be performed and compared with experiments. Did you find this useful? Give us your feedback ##### Citations More filters Journal ArticleDOI TL;DR: In this article, a smoothed particle hydrodynamics model with numerical diffusive terms is used to analyze violent water flows and boundary conditions on solid surfaces of arbitrary shape are enforced with a new technique based on fixed ghost particles. 535 citations Journal ArticleDOI TL;DR: In this paper, the authors proposed modified MPS methods for the prediction of wave impact pressure on a coastal structure by introducing new formulations for the pressure gradient and a new formulation of the source term of the Poisson Pressure Equation (PPE). 288 citations Journal ArticleDOI TL;DR: In this paper, the authors conduct experimental measurements on a dam break flow over a horizontal dry bed in order to provide a detailed insight, with emphasis on the pressure loads, into the dynamics of the dam break wave impacting a vertical wall downstream the dam. 236 citations ### Cites background or methods from "A Volume-of-Fluid based simulation ..." • ...(2010) (sensor 2), 30 mm Kleefsman et al. (2005) (sensor 3) and 80 mm Zhou et al. (1999) (sensor 4) above the tank bed.... [...] • ..., 2004)(Kleefsman et al., 2005) were approved as two benchmark tests for validation of numerical codes by Smoothed Particle Hydrodynamics European Research Interest Community (interest group, ????).... [...] • ...However, description of applied pressure transducers was missing in both (Kleefsman et al., 2005)(Wemmenhove et al., 2010).... [...] • ...The MARIN experimental setup was also used in the work of Wemmenhove et al. (Wemmenhove et al., 2010) and Kleefsmann et al. (Kleefsman et al., 2005).... [...] • ...However, description of applied pressure transducers was missing in both (Kleefsman et al., 2005)(Wemmenhove et al.... [...] Journal ArticleDOI 13 Jun 2011-PLOS ONE TL;DR: Both the achieved speed-ups and the quantitative agreement with experiments suggest that CUDA-based GPU programming can be used in SPH methods with efficiency and reliability. Abstract: Smoothed Particle Hydrodynamics (SPH) is a numerical method commonly used in Computational Fluid Dynamics (CFD) to simulate complex free-surface flows. Simulations with this mesh-free particle method far exceed the capacity of a single processor. In this paper, as part of a dual-functioning code for either central processing units (CPUs) or Graphics Processor Units (GPUs), a parallelisation using GPUs is presented. The GPU parallelisation technique uses the Compute Unified Device Architecture (CUDA) of nVidia devices. Simulations with more than one million particles on a single GPU card exhibit speedups of up to two orders of magnitude over using a single-core CPU. It is demonstrated that the code achieves different speedups with different CUDA-enabled GPUs. The numerical behaviour of the SPH code is validated with a standard benchmark test case of dam break flow impacting on an obstacle where good agreement with the experimental results is observed. Both the achieved speed-ups and the quantitative agreement with experiments suggest that CUDA-based GPU programming can be used in SPH methods with efficiency and reliability. 219 citations ### Cites background or methods from "A Volume-of-Fluid based simulation ..." • ...Right snapshots correspond to figures from [44].... [...] • ...The experiment [44] provides water heights and pressure measurements at different locations.... [...] • ...The experiment for validation The experiment described in [44] consists of a dam break flow impacting with an obstacle.... [...] • ...Experimental configuration of the [44] experiment and measeurement positions for the experimental data: Side view, top view and location of pressure sensors.... [...] Journal ArticleDOI Martin Ferrand TL;DR: Ferrand et al. as discussed by the authors proposed a new approach based on a renormalising factor for writing all boundary terms, which significantly improved traditional wall treatment in smoothed particle hydrodynamics, for pressure forces, wall friction and turbulent conditions. Abstract: Computational / Numerical Methods > International Journal for Numerical Methods in Fluids > Early View > AbstractJOURNAL TOOLS Get New Content Alerts Get RSS feed Save to My Profile Get Sample Copy Recommend to Your Librarian JOURNAL MENU Journal HomeFIND ISSUES Current Issue All IssuesFIND ARTICLES Early View Most Accessed Most CitedGET ACCESS Subscribe / Renew FOR CONTRIBUTORS OnlineOpen Author Guidelines Submit an ArticleABOUT THIS JOURNAL News Overview Editorial Board Permissions Advertise ContactSPECIAL FEATURES La Tex Class File Numerical Engineering Backfile Collection FREE Engineering sample issues Other Numerical Engineering Journals Biomedical Engineering Zienkiewicz Award 2010 Keyword Tag Cloud Wiley Job Network To our Authors: newsletterNumerical Engineering JournalsWiley Job NetworkWiley Science Jobs Hydraulics Product Support Engineer Salary: �35000 - �40000 per annum + package Location: Berkshire Product Support Engineer - Hydraulics Permanent Salaried with company benefits package Working as part of a dedicated Sales? Mechanical Engineer - Sensor Components - Oil & Gas - Cambridge Salary: Negotiable Location: Cambridge Mechanical Engineer - Sensor Components - Oil & Gas - Cambridgeshire One of the country's leading Scientific Consultancy's f? Usability Design Engineer - Medical Devices - Creative Design Salary: Negotiable Location: Hertfordshire Usability Design Engineer - Medical Devices - Creative Design - Hertfordshire Are you an experienced Design Engineer with a ? Product Designer - Draughtsperson - SolidWorks - Hertfordshire Salary: Negotiable Location: Hertfordshire Product Designer - Draughtsperson - SolidWorks - Hertfordshire A well established Medical Devices Company, based in Hertford? Research Engineer - Medical Imaging - Edinburgh - Photonics Salary: Negotiable Location: Edinburgh Research Engineer - Medical Imaging - Edinburgh - Photonics A growing Medical Imaging company, based in Edinburgh, are curre?Find more Wiley Science Jobs �Research ArticleUnified semi-analytical wall boundary conditions for inviscid, laminar or turbulent flows in the meshless SPH method M. Ferrand1,*, D.R. Laurence2, B.D. Rogers2, D. Violeau3, C. Kassiotis3Article first published online: 20 MAR 2012DOI: 10.1002/fld.3666Copyright � 2012 John Wiley & Sons, Ltd.IssueCover image for Vol. 70 Issue 6International Journal for Numerical Methods in FluidsEarly View (Online Version of Record published before inclusion in an issue)Additional Information(Show All)How to CiteAuthor InformationPublication HistorySEARCHSearch ScopeSearch String Advanced > Saved Searches >ARTICLE TOOLS Get PDF (12213K) Save to My Profile E-mail Link to this Article Export Citation for this Article Get Citation Alerts Request PermissionsMore Sharing ServicesShare|Share on citeulikeShare on connoteaShare on deliciousShare on www.mendeley.comShare on twitter Abstract Article References Cited ByView Full Article (HTML) Get PDF (12213K)Keywords: wall boundary conditions; Navier?Stokes; smoothed particles hydrodynamics; Lagrangian; laminar flow; viscous flows; turbulent flow; meshfree; free surfaceSUMMARYWall boundary conditions in smoothed particle hydrodynamics (SPH) is a key issue to perform accurate simulations. We propose here a new approach based on a renormalising factor for writing all boundary terms. This factor depends on the local shape of a wall and on the position of a particle relative to the wall, which is described by segments (in two-dimensions), instead of the cumbersome fictitious or ghost particles used in most existing SPH models. By solving a dynamic equation for the renormalising factor, we significantly improve traditional wall treatment in SPH, for pressure forces, wall friction and turbulent conditions. The new model is demonstrated for cases including hydrostatic conditions for still water in a tank of complex geometry and a dam break over triangular bed profile with sharp angle where significant improved behaviour is obtained in comparison with the conventional boundary techniques. The latter case is also compared with a finite volume and volume-of-fluid scheme. The performance of the model for a two-dimensional laminar flow in a channel is demonstrated where the profiles of velocity are in agreement with the theoretical ones, demonstrating that the derived wall shear stress balances the pressure gradient. Finally, the performance of the model is demonstrated for flow in a schematic fish pass where both the velocity field and turbulent viscosity fields are satisfactorily reproduced compared with mesh-based codes. 208 citations ##### References More filters Journal ArticleDOI TL;DR: The PSC algorithm as mentioned in this paper approximates the Hamilton-Jacobi equations with parabolic right-hand-sides by using techniques from the hyperbolic conservation laws, which can be used also for more general surface motion problems. 13,020 citations Journal ArticleDOI TL;DR: In this paper, the concept of a fractional volume of fluid (VOF) has been used to approximate free boundaries in finite-difference numerical simulations, which is shown to be more flexible and efficient than other methods for treating complicated free boundary configurations. 11,567 citations ### "A Volume-of-Fluid based simulation ..." refers methods in this paper • ...This is done using an adapted version of the VOF method first introduced in [23].... [...] • ...The procedure explained by Hirt and Nichols in [23] is used to calculate these fluxes.... [...] • ...6 the resulting free surface after one rotation is shown with on the left standard VOF of Hirt–Nichols and on the right the current method.... [...] • ...The pressure in surface cells is now calculated from the pressure at the free surface and the pressure in an adjacent fluid cell as in the method of Hirt and Nichols [23].... [...] • ...In the VOF method, first introduced by Hirt and Nichols [23], a VOF function F is introduced with values between zero and one, indicating the fractional volume of a cell that is filled with a certain fluid.... [...] Journal ArticleDOI TL;DR: In this paper, a new technique is described for the numerical investigation of the time-dependent flow of an incompressible fluid, the boundary of which is partially confined and partially free The full Navier-Stokes equations are written in finite-difference form, and the solution is accomplished by finite-time step advancement. Abstract: A new technique is described for the numerical investigation of the time‐dependent flow of an incompressible fluid, the boundary of which is partially confined and partially free The full Navier‐Stokes equations are written in finite‐difference form, and the solution is accomplished by finite‐time‐step advancement The primary dependent variables are the pressure and the velocity components Also used is a set of marker particles which move with the fluid The technique is called the marker and cell method Some examples of the application of this method are presented All non‐linear effects are completely included, and the transient aspects can be computed for as much elapsed time as desired 5,841 citations ### "A Volume-of-Fluid based simulation ..." refers background or methods in this paper • ...The variables are staggered as in the original MAC method [20], which means that the velocities are defined on cell faces, whereas the pressure is defined in cell centres.... [...] • ...Method 1: the divergence of every S-cell is set to zero [20].... [...] Book 01 Jan 2007 3,984 citations Journal ArticleDOI TL;DR: In this article, the authors extended previous work on the solution of the Navier-Stokes equations in the presence of moving immersed boundaries which interact with the fluid and introduced an improved numerical representation of the δ-function. 2,968 citations ### "A Volume-of-Fluid based simulation ..." refers methods in this paper • ...One is the immersed boundary method developed by Peskin [34] since the 1970s and used, for e.g. [26,27].... [...] • ...One is the immersed boundary method developed by Peskin [34] since the 1970s and used, for e.... [...] ##### Frequently Asked Questions (17) ###### Q1. What contributions have the authors mentioned in the paper "University of groningen a volume-of-fluid based simulation method for wave impact problems" ? In this paper, some aspects of water impact and green water loading are considered by numerically investigating a dambreak problem and water entry problems. To track the free surface the VOF function Fs is used, which is 0 if no fluid is present in the cell, 1 if the open part of the cell is completely filled with fluid and between 0 and 1 if the cell is partly filled with fluid. To investigate convergence of the method under grid refinement, the circular cylinder entry simulations have also been run with different grids. In the case of uncut cells with fixed objects the stability of the equation containing the time integration term and the convective term is given by the CFL-restriction, which in one dimension reads dtjuj/h 6 1 (h is the size of the uncut cell). The second series of simulations consists of drop tests, where wedges with different dead-rise angles, a circular cylinder and a cone are entering the water. when performing wave simulations an inflow boundary is needed where the incoming wave is prescribed, and at the opposite boundary a non-reflecting outflow condition should be used. Using this discretisation, the discrete gradient operator in the pressure term is the negative transpose of the discrete divergence operator Eq. (9), which is also an analytic property ($ = (\$Æ)T) [43].

At the outflow boundary a non-reflecting boundary condition is needed to prevent the waves from reflecting from the boundary into the domain.

The weight factor K has been chosen such that the stabilising term is only used when the body is moving; note that it equals unity for fixed objects.

The boundary condition that defines the pressure at the free surface is given by Eq. (5), which describes the continuity of normal stresses.

For the initial state of the entry of a circular cylinder, the hydrodynamic slamming force can be estimated byF ¼ V q p 2 ð2VR 2V 2tÞ; ð46Þwhere t denotes time with t = 0 the moment of first impact.

To assess the performance of the displacement method, the standard advection test of rotation of a slotted disk, defined by Zalesak [46], has been performed.

The X in this estimation cancels the contribution of X 1 in X 1C, leaving the stability criterion for cut cells the same as for uncut cells [10].

In two dimensions the integrand can be written as the sum of the second order horizontal and verticalderivatives of the horizontal velocity.

In the current method the objects are either moving with a prescribed motion, or the motion is calculated from the interaction between the object and the liquid dynamics.

The slamming coefficient is given by Cs = F/qRV2, with F the total vertical hydrodynamic force, R the radius of the circular cylinder and V the entry velocity.

The direction in which the function is defined is the direction of the coordinate axis that is most normal to the free surface (which is the positive z-direction in Fig. 5).