Unstructured finite volume method for water impact on a rigid body*

YU Yan (余艳)

College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China

School of Science, Harbin University, Harbin 150086, China, E-mail: yuyan801206@126.com

MING Ping-jian (明平剑)

College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China

DUAN Wen-yang (段文洋)

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China

Unstructured finite volume method for water impact on a rigid body*

YU Yan (余艳)

College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China

School of Science, Harbin University, Harbin 150086, China, E-mail: yuyan801206@126.com

MING Ping-jian (明平剑)

College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China

DUAN Wen-yang (段文洋)

College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China

(Received May 3, 2013, Revised September 11, 2013)

A new method is presented for the water impact simulation, in which the air-water two phase flow is solved using the pressure-based computational fluid dynamics method. Theoretically, the air effects can be taken into account in the water structure interaction. The key point of this method is the air-water interface capture, which is treated as a physical discontinuity and can be captured by a well-designed high order scheme. According to a normalized variable diagram, a high order discrete scheme on unstructured grids is realised, so a numerical method for the free surface flow on a fixed grid can be established. This method is implemented using an in-house code, the General Transport Equation Analyzer, which is an unstructured grid finite volume solver. The method is verified with the wedge water and structure interaction problem.

air cushion, numerical simulation, unstructured grid, volume of fluid (VOF)

Introduction

Fluid/structure impact problems can be found widely in marine and ship engineering applications. In many cases, the process in which the body of the fluid bounded by a free surface impacts on approaching structures or other fluids, is affected by air cushioning. Many analytical solutions were obtained with simplified assumptions. Wood et al.[1]focused on the water alone with approximations of the pressure-impulse theory for the impact process. Howison et al.[2]theoretically studied the air cushioning effects in the gap between the fluid and the structure. Another example was the case considered by Smith et al.[3]for liquid columns moving towards each other, but the focus was on the air cushion and on the period until “touchdown”. A further example was the case of a droplet impact on a thin fluid layer during the initial stages[2]. Hicks and Purvis[4]studied the air cushioning effects on a droplet impact on a structure or fluid. Lan et al.[5]studied the wave impact on slab using both experimental and numerical methods. Ding et al.[6]discussed the the characteristics of the wave impact forces on the undersides of the rectangular structures with various length-breadth ratios. The last several decades saw a great progress in this field. Initially, the effects of the gravity and the viscosity of the liquid are usually ignored because the impact usually lasts for a very short period of time. When the compressibility of the liquid is ignored, the velocity potential can be introduced. The problem of a two-dimensional wedge entering water can thus be studied step by step, from a symmetrical wedge to an asymmetrical wedge, from constant velocity to free fall, and from a rigid body to an elastic wedge body. So this problem is used to verify the volume of fluid (VOF) method in this study.

In the present research, a numerical method is adopted to study the process of water impact on a structure taking into account the air cushioning effect.There are two normal interface reconstruction algorithms for the VOF method. One is based on the simple line interface construction (SLIC) algorithm, and the other is based on the piecewise linear interface construction (PLIC) algorithm. Yang and James[7]simulated the interaction between an extreme wave and a freely-floating structure using the VOF method based on an unstructured grid. Yang et al.[8]proposed an unsplit Lagrangian advection (ULA) scheme based on PLIC, which was implemented with two dimension structured grids. However, these method could not guarantee that the free surfaces calculated on either side of the cell face would completely match each other. Thus it may be difficult to make the calculated result converge, and the result may even diverge. A number of improved methods were proposed, with good results. The advantages and disadvantages of various numerical methods were reviewed in Ref.[9]. Rebouillat and Liksonov[10]used three methods including the direct numerical method with FEM and FDM, the smooth particle method and the (semi-)analytical method to study the fluid-structure interaction in a partially filled liquid container, but without mentioning the VOF based finite volume method.

In the VOF method, the free surface is reconstructed using a geometric method. However, reconstructing the surface requires ample computation resources, and is based mainly on the advancement of the structured grid. This paper proposes a new concept in which the free surface is treated as a physical discontinuity and can be captured by a well-designed highorder scheme similar to the shock capturing scheme on unstructured grids. Therefore, the sloshing process can be simulated on fixed unstructured grids.

1. Mathematical models and numerical methods

1.1 Mathematical models

We consider the problem of a two-dimensional liquid column in the air, approaching a solid boundary at a high speed, with a large impact. Both the liquid and the air are assumed to be Newtonian fluids and the governing equations are of the Navier-Stokes type. The fluid pressure and velocity are continuous at the interface of the liquid and gas, but the density and the viscosity are discontinuous and in fact change sharply.

In the Cartesian systemoxy , the governing equations based on the mass and momentum conservations can be written as:

where U,ρbandpare the velocity vector, the body force vector and the pressure, andρandµare the density and the dynamic viscosity of the fluid. When the fluid domain is discretized into small cells, we define the liquid volume fraction asαin each cell, which is the ratio of the volume taken by the liquid in the cell and the volume of the cell itself. Thus αis 1 for a cell full of water and 0 for a cell full of air. For a cell containing both liquid and air0＜α＜1. The evolution ofαin each cell is based on the following equation

Within each cell the density and the viscosity in Eq.(2) are obtained by using the following weighting equations

where the subscripts land gindicate the liquid and the gas, respectively. For an incompressible fluid, Eq.(1) can be written as

This is used in this study unless stated otherwise.

1.2 Numerical method

The above problem is solved by using the finite volume method and Eq.(2) can be written as

where (x1,x2)=(x,y), and u1and u2are the components of the velocity vector in thexandy directions.

These equations can be integrated over a cellΩbounded by the boundary ∂Ωwhich consists of faces Afj,j=1,2…m. Thus using the Gauss theorem, we have

Integrating the convection term over the boundary, this can be approximated as

where0≤λ≤1. The subscriptsHO andFUDin the equation indicate that the value is obtained from the high order and the first order upwind methods, respectively, based on the sign offj. This gives

Based on the definition of the first order upwind, the first term on the right hand side of the above equation can be rewritten as

where sgn(fj)=1means that fj＞0, that is, the flow leaves the cell, marked ascpin Fig.1, from face j and sgn(fj)=-1means that fj＜0, the flow enters the cell from the adjacent cellcjas shown in Fig.1. Because of the mass conservation, the summation in the first brackets on the right hand side is zero. Thus the convection term can be written as

The advantage of this equation is that the term before uicpis also positive. This will improve the stability of the solution.

For the diffusion term in Eq.(8), we have

To deal with the normal derivative on the face in the above equation, we write

where njand τjare the normal vector and the unit vector in the tangential direction of the face, respectively. Multiplying both sides of the equation by the vectordjlinking the centres of the two cells on both sides of the face (see Fig.1), we obtain

We can also rewrite Eq.(15) in the following form

Inserting Eq.(17) into Eq.(16), we obtain

Substituting Eq.(18) into Eq.(14) gives

The remaining velocity gradient on the right hand side of Eq.(19) can be obtained from the following equation

The first term on the right hand side defines the matrix coefficient and the second term can be dealt with as a source term, known as a deferred term.

For the time derivative term in Eq.(8), we use

in which the superscript indicates the time step.

Substituting Eqs.(13), (19) and (21) into Eq.(8) with an implicit scheme for the convection and diffusion terms, we obtain

The source term and the coefficients at the nthtime step are obtained by using the values taken from the previous iteration. The initial value of the iteration is taken from the previous time step solution. Eq.(22) is then solved using a conjugate gradient method.

The velocity obtained from the solution of Eq.(22), based on the assumed pressure distribution p∗, can be written as ui∗. It does not automatically satisfy the mass conservation condition and a correctionp′should be applied to the assumed pressure. Thus we can write:

Substituting Eqs.(23) and (24) into Eq.(22), we obtain

As is assumed, the fluid flow is mainly driven by the pressure gradients, and the neighbour cell’s influence is ignored. This is the baseline of the algorithm of the semi implicit method for the pressure linked equation (SIMPLE). Many modified forms are proposed, e.g., SIMPLEC.

To enforce the mass continuity using the finite volume method, we integrate Eq.(6) over the cellΩ. The Gauss theorem gives:

We can write the velocity correction expression in a similar way as with Eq.(26).

Inserting Eqs.(30) and (31) into Eq.(29), we obtain the pressure correction equation.

We rewrite the above equation and insert the velocity correction with the pressure correction to them to obtain the pressure correction equation

The linear equation set, Eq.(33), can be solved iteratively, using the conjugate gradient method. To accelerate the solution process, an algebraic multi-grid method is adopted. The pressure, the velocity and the flux are updated according to Eqs.(24), (26) and (28), respectively after the pressure correctionp′is solved.

For the collocated grid finite volume method, all variables are stored at the cell centre. The convection term value required on the cell face is obtained through interpolation. There are many schemes available. As an illustration, a one-dimensional problem with a uniform grid is used to explain the normalized variable diagram (NVD) procedure. Based on the NVD, several total variation diminishing (TVD) schemes can be designed. A one-dimensional example is shown in Fig.2. The cell face value of a parameterφ,φf, is now considered. According to the flow direction on the cell face, the values of the same parameter in the upstream and downstream cells are denoted byφUand φD, and the value in the further upstream cell (or second upstream cell) by φSU, as shown in Fig.3. The normalized variable can be defined as

which is 0 at cellSU and 1 at cellD. With the same variable value on the downstream and second upwind stream cells, the normalized variable is constant at 1, so we can now just consider the situation near a free surface where the liquid volume fraction is different on the downstream and second upwind stream cells. Therefore, we have at the upstream cell

and at the cell interface

Different discrete forms can be expressed by the NVD as shown below in Fig.4. More details can be found in Ref.[13].

In different schemes, different values at the cell centres are used forSU,U and Dto interpolate the value atf. When it is normalized,φ～is 0 and 1 atSUandD . Thus φ～fshould be a function of φ～Uonly, or:

Once the value of φ～fis obtained from one of the above schemes, the value ofφfitself can then be obtained from Eq.(36) as

Inserting this into Eq.(38), we obtain

For a two-dimensional unstructured grid, the line linking cell centresDandUmay not pass through the centre of the next cell. The way to choose the pointSUis to extend the straight line from DU to SU , and the distance fromU toSU is that from DtoU . The value φSUis then obtained by using the following equation

To ensure that φSUis bounded, we use

where φmaxand φminare 1 and 0, respectively for the volume fraction Eq.(3). The gradient of the variation is obtained based on an equation similar to Eq.(20).

Onno Ubbink (1997) recommended the correlation of the expression for Eq.(11) as follows. This formulation is used for the free surface tracking process.

Integrating Eq.(3) over a cell and using the Gauss theorem for the convection term, we obtain

The time derivative term is discretized by using the Crank-Nicolson scheme to obtain a discretized form.

Inserting Eq.(40) into Eq.(44), we obtain

where

1.3 Solution procedure

The above problem is solved by using the following procedure as shown in Fig.5. First, the volume fraction Eq.(3) is solved. Second, Eqs.(4) and (5) are used to calculate the face density and the viscosity. Finally, the SIMPLE algorithm is used to solve Eqs.(2) and (6).

2. Simulation model

The computational domain is shown in Fig.6. The width and the height are 30 m and 15 m, respectively, and the wedge water height is 4 m. The deadrise angles areαandβas shown in the figure. To generate a high quality grid, the domain is divided into several sub-domains as in Fig.7. The mesh is refined in the sub-domain A, where the wedge water impacts occur on the structure. In the other sub-domains a quadrilateral mesh is generated.

3. Simulation and discussion

3.1 Grid-independent study

The case of a symmetric water wedge with the dead-rise angle equal to 60ois used to analyze the convergence of results with different grids. Because of the symmetry of the computational domain and the short period for the impact process, a symmetrical model is adopted here. The grids are generated by using the commercial software GAMBIT. Three nonuniform quadrilateral meshes with 7 500, 30 000 and 52 500 cells are used. The time step is adjusted by the courant number. The shape of the water surface and the pressure on the body are given in Fig.8, in whichAs can be seen, different grids do not significantly alter the trends of the results. The curve obtained with the grid one deviates from that obtained with the other two grids, but the curves obtained with the grids two and three are in good agreement. The computation time seems linearly proportional to the cell number. Therefore, a mesh similar to the second grid is now used in this study.

3.2 Comparison with other published results

The wedge water column impacting on the structure is extensively found in marine engineering and was widely studied. Here we use this case to verify the present method. The simulation results are shown in Figs.9 and 10. The dimensionless pressure distributions in dimensionless coordinate systems agree well with each other for different times and impacting velocities. This is consistent with the similarity theory of Wu[11].

3.3 Wedge water impact on the structures

In Figs.11 and 12, we show the free surface profile and the pressure distribution on a horizontal wall hit by a liquid wedge defined in the method used in this study based on the potential theory method from Duan et al.[14]and Zhang et al.[15]. From Fig.11, it can be seen that the free surface profiles in this study agree well overally with the data from the literature, although there is some local discrepancy with the results from Zhang et al.[15]especially near or in the jet region. Good agreements are achieved between our study and that of Duan et al.[14], which indicates that the nonlinear characteristics affect the free surface shape near the jet region. The discrepancy in Fig.12 is therefore believed to be due to the approximation based on an exponential function used for the free surface profile by Zhang et al.[15].

3.4 Simulation with a semi-elliptic water column front

A water column with a semi-elliptic air cushion is now studied by using the numerical methods. The width and the height are 0.1 m and 0.15 m, respectively, while the bottom boundary of the water column is prescribed by the equation

The pressure history at the centre of the structure surface is shown in Fig.13. It is shown that the pressure on the structure under the action of a semi elliptic bottomed water column is much more complex than that for a flat bottomed water column. There are three main peaks during the water column impact processes, while at the same time some spikes also appear beside the main peaks. These phenomena can be explained by the detail of the free surface movement and the velocity distribution as shown in Fig.14. The first peak is at 0.003 s when the bottom right hand corner of the water column approaches the surface of the structure. As shown in Fig.13, the air cushion has an important effect on the impact process at this stage. Without any air shear force, the water column remains in a constant shape and the bottom right hand corner hits the structure surface at nearly 0.003 s. A large pressure peak then occurs while the pressure is much smoother for the case with an air cushion. The free surface shape substantially changes because of the air flow shear force in Fig.14. As shown in Fig.13, the case with an elliptic water column front is much more complex than that with a flat front. The peak pressure is overestimated by neglecting the air cushion. The descending velocity of the main part of the water column does not change too much with the air flow and the gravity and this supports the assumption made by some researcher to ignore the effects of the gravity. However, the variation in the water column front shape is crucial to the impact pressure as shown in Fig.14. The tip hits the structure surface at 3 ms without an air cushion and the first pressure peak occurs, while the tip moves a little to the right, which causes a delay of the first pressure peak. The situation is similar for the second and third pressure peaks.

4. Conclusions

This paper presents the development of a finite volume high order numerical algorithm based on the NVD to simulate the water impact on a rigid structure with unstructured grids. The basic concept of this algorithm is that the free surface can be treated as a physical discontinuity which can be captured by a welldesigned high order scheme. Owing to the above methods, the water column impact on rigid structures can be simulated with unstructured fixed grids. The model in this study is implemented by using an unstructured grid finite volume solver.

The proposed method is verified by comparing the numerical results with the potential theoretical results. Both the free surface and the pressure are in good agreement.

Acknowledgement

This work was supported of the China Scholarship Council (Grant No. 2010307245).

[1] WOOD D. J., PEREGRINE D. H. Wave impact on a wall using pressure-impulse theory. II: Porous berm[J]. Journal of Waterway, Port, Coastal, Ocean Engineering, 2000, 126(4): 191-195.

[2] HOWISON S. D., OCKENDON J. R. and OLIVER J. M. et al. Droplet impact on a thin fluid layer[J]. Journal of Fluid Mechanics, 2005, 542: 1-23.

[3] SMITH F. T., LI L. and WU G. Air cushioning with a lubrication/inviscid balance[J]. Journal of Fluid Mechanics, 2003,482: 291-318.

[4] HICKS P. D., PURVIS R. Air cushioning and bubble entrapment in three-dimensional droplet impacts[J]. Journal of Fluid Mechanics, 2010, 649: 135-163.

[5] LAN Ya-mei, GUO Wen-hua and LIU Hua et al. Numerical simulation of wave impact on the slab[J]. Journal of Hydrodynamics, 2010, 22(5Suppl.): 986-992.

[6] DING Zhao-qiang, WANG Guo-yu and REN Bing. Three-dimensional numerical simulation of wave slamming on an open structure[J]. Journal of Hydrodynamics, 2012, 24(4):526-534.

[7] YANG X., JAMES A. J. Analytic relations for reconstructing piecewise linear interfaces in triangular and tetrahedral grids[J]. Journal of computational physics, 2006, 214(1): 41-54.

[8] YANG Wei, LIU Shu-hong and WU Yu-lin. An unsplit Lagrangian advection scheme for volume of fluid method[J]. Journal of Hydrodynamics, 2010, 22(1): 73-80.

[9] IBRAHIM R., PILIPCHUK V. and IKEDA T. Recent advances in liquid sloshing dynamics[J]. Applied Mechanics Review, 2001, 54(2): 133-199.

[10] REBOUILLAT S., LIKSONOV D. Fluid-structure interaction in partially filled liquid containers: A comparative review of numerical approaches[J]. Computers and Fluids, 2010, 39(5): 739-746.

[11] WU G., HU Z. Z. A Taylor series based finite volume method the Navier-Stokes equations[J]. International Journal for Numerical Methods in Fluids, 2008, 58(12):1299-1325.

[12] UBBINK O. Numerical prediction of two fluid systems with sharp interfaces[D]. London, UK: Imperial College of Science, Technology and Medicine, 1997.

[13] LEONARD, B. P. The ULTIMATE conservative difference scheme applied to unsteady one dimensional advection[J]. Computational method in applied mecha- nics and engineering, 1991, 88(1): 17-74.

[14] DUAN W., XU G. and WU G. Similarity solution of oblique impact of wedge-shaped water column on wedged coastal structures[J]. Coastal Engineering, 2009, 56: 400-407.

[15] ZHANG S., YUE D. K. P. and TANIZAWA K. Simulation of plunging wave impact on a vertical wall[J]. Journal of Fluid Mechanics, 1996, 327: 221-254.

10.1016/S1001-6058(14)60061-5

* Project supported by the National Natural Science Foundation of China (Grant Nos. 51206031, 51079032), the China Postdoctoral Science Foundation funded Project (Grant No. 20100471016).

Biography: YU Yan (1980-), Female, Ph. D. Candidate

MING Ping-jian,

E-mail:pingjianming@hrbeu.edu.cn