Unstructured finite volume method for water impact on a rigid body

538Available online at www.sciencedirect.comScienceDirectHT充DJournal of HydrodynamicsEL SEVIER2014,26(4):538-548www.sciencedirect.com/science/journal/l001 6058DOI: 10.1016/S 1001 -6058( 14)60061-5Unstructured finite volume method for water impact on a rigid bodyYU Yan (余艳)College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, ChinaSchool of Science, Harbin University, Harbin 150086, China, E mail: yuyan801 206@ 126.comMING Ping-jian (明平剑)DUAN Wen-yang (段文洋)College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China(Received May 3, 2013, Revised September 11, 2013)Abstract: A new method is presented for the water impact simulation, in which the air- water two phase flow is solved using thepressure-based computational fluid dynamics method. Theoretically, the air effects can be taken into account in the water structureinteraction. The key point of this method is the air-water interface capture, which is treated as a physical discontinuity and can becaptured by a well-designed high order scheme. According to a normalized variable diagram, a high order discrete scheme onunstructured grids is realised, so a numerical method for the free surface flow on a fixed grid can be established. This method isimplemented 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.Key words: air cushion, numerical simulation, unstructured grid, volume of fluid (VOF)Introductionwas on the air cushion and on the period until“touch-Fluid/structure impact problems can be found wi-down". A further example was the case of a dropletdely in marine and ship engineering applications. Iimpact on a thin fluid layer during the initial stages' .many cases, the process in which the body of the fluidHicks and Purvis4I studied the air cushioning effectsbounded by a free surface impacts on approaching st-on a droplet impact on a structure or fluid. Lan et al.!5ructures or other fluids, is affected by air cushioning.studied the wave impact on slab using both experi-Many analytical solutions were. obtained with simpli-mental and numerical methods. Ding et al.'0 discussedfied assumptions. Wood et al.!l focused on the waterthe the characteristics of the wave impact forces onalone with approximations of the pressure-impulsethe undersides of the rectangular structures with va-theory for the impact process. Howison et al.!4I theore-rious length-breadth ratios. The last several decadestically studied the air cushioning effects in the gapsaw a great progress in this field. Initially, the effectsbetween the fluid and the structure. Another exampleof the gravity and the viscosity of the liquid are usua-was the case considered by Smith et al.5 for liquidlly ignored because the impact usually lasts for a verycolumns moving towards each other, but the focusshort period of time. When the compressibility of theliquid is ignored, the velocity potential can be introdu-* Project supported by the National Natural Science Foun-ced. The problem of a two-dimensional wedge ente-dation of China (Grant Nos. 51206031, 51079032), the Chinaring water can thus be studied step by step, from aPostdoctoral Science Foundation funded Project (Grant No.symmetrical wedge to an asymmetrical wedge, fromconstant velocity to free fall, and from a rigid body toBiography: YU Yan (1980-), Female, Ph. D. Candidatean elastic wedge body. So this problem is used to veri-Corresponding author: MING Ping-jian,fy the volume of fluid (VOF) method in this study.E-mail: pingjianming @ hrbeu.edu.cnIn the pre中国煤化工1 method isadopted to stud|YHCNMH Gpact on ast-ructure taking ....。...oning effect.539There are two normal interface reconstruction algori-V.{u[VU +(VU)I]}+ ρb(2)thms for the VOF method. One is based on the simpleline interface construction (SLIC) algorithm, and thewhere U，pb and p are the velocity vector, theother is based on the piecewise linear interface const-body force vector and the pressure, and ρ and μ areruction (PLIC) algorithm. Y ang and James' simula-ted the interaction between an extreme wave and athe density and the dynamic viscosity of the fluid.freely-floating structure using the VOF method basedWhen the fluid domain is discretized into small cells, .on an unstructured grid. Y ang et al.s proposed an un-we define the liquid volume fraction as a in eachsplit Lagrangian advection (ULA) scheme based oncell, which is the ratio of the volume taken by the li-PLIC, which was implemented with two dimension st-quid in the cell and the volume of the cell itself. Thusructured grids. However, these method could not gua-a is 1 foracell full of water and 0 for a cell full ofrantee that the free surfaces calculated on either sideair. For a cell containing both liquid and air 0Zaf[(n)no - ()rum)(13)through the cell face j. The value of the velocity inthe summation of Eq.(9) is taken at the intersection ofeach face with the line linking the centres of the twoThe advantage of this equation is that the term beforeadjacent cells. The accuracy of such an approximationu.。is also positive. This will improve the stability ofwas discussed in Ref.[11]. The required value on thethe solution.cell face is then interpolated from those at the cellcentres, which is another layer of approximation'There are many different interpolation methods. Herewe use .工(u)f = (u.)puD + Aa[(u)Ho - (u})pw]cPwhere 0≤A≤1. The subscripts HO and FUD irthe equation indicate that the value is obtained fromthe high order and the first order upwind methods, re-Fig.1 The grid structure and the normal and tangential direc-spectively, based on the sign of f;. This givestions of the cell faceFor the diffusion term in Eq.(8), we have之(u)pf=之f(n)r + SZzf[(4)mo - ()pum)ou;JupNVu,.dA=二(vVu)oAn=Z(μA。=Ion )。(14)Based on the definition of the first order upwind,the first term on the right hand side of the above equa-To deal with the normal derivative on the face intion can be rewritten asthe above equation, we writeZf(u)pum= 2 f4mp+ 艺fμug=(Vu,)g =Ou( ou,(15)j=l.sgn(f; )=1n; +|0τ )后lxp |艺f,+Ef,|-uZf,+where n; and T; are the normal vector and the unitvector in the tangential direction of the face, respec-tively. Multiplying both sides of the equation by thevector d; linking the centres of the two cells on bothfj+习fsides of the face (see Fig. 1), we obtain=l,sgI(Vu;)ndTd,( ou,(16)(12)n;dn;d;( 0t);where sgn(f;)=1 means that f; > 0, that is, the flowWe can also rewrite Eq.(15) in the following formleaves the cell, marked as cp in Fig.1, from face jOu )(17)and sgn(f;)=-1 means that f; <0, the flow enters中国煤化工the cell from the adjacent cell cj as shown in Fig.1.Because of the mass conservation, the summation inInserting Eq.(11MHCNMHG541 .d」 (0u、The source term and the coefficients at the nth time( On), n,d,\ad), n;d,(Vu)g=step are obtained by using the values taken from theprevious iteration. The initial value of the iteration istaken from the previous time step solution. Eq.(22) isthen solved using a conjugate gradient method.Maj - "icpn;d(Vu)p .(18)The velocity obtained from the solution ofEq.(22), based on the assumed pressure distributionp*, can be written as u*. It does not automaticallySubstituting Eq.(18) into Eq.(14) givessatisfy the mass conservation condition and a correc-tion p' should be applied to the assumed pressure.」a 1Vvn. dA=ZHA-(Ul -up)+Thus we can write:台n, d;u;=u; +ul' .(23);| n;(Vu)A;(19)p=p° +p'(24)The remaining velocity gradient on the right handSubstituting Eqs.(23) and (24) into Eq.(22), we obtainside of Eq.(19) can be obtained from the followingap’' oequationaia ,u=> agug-L(25)(Vu)月=w,(Vu)p +(1 - w,)(Vu)gAs is assumed, the fluid flow is mainly driven by thewhere W, 台xx|xcr-xop| is the weightingpressure gradients, and the neighbour cell's influenceis ignored. This is the baseline of the algorithm of thecoefficient based on the distance. The gradient at thesemi implicit method for the pressure linked equationcell centre is calculated by using the Gaussian method.(SIMPLE). Many modified forms are proposed, e.g.,The first term on the right hand side defines theSIMPLEC.matrix coefficient and the second term can be dealtwith as a source term, known as a deferred term.For the time derivative term in Eq.(8), we useanu'm=-_s(26)[ (u"m)d2= (u)"-(u)"_sTo enforce the mass continuity using the finiteJa Otvolume method, we integrate Eq.(6) over the cell L.The Gauss theorem gives: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 diffu-，f V.UdL=φ UdA=>F;=0(27)sion terms, we obtainF2 aguls +8a(22)F,=Fj+F,(28)=1Equation (27) then giveswhereay =-Min(0,f)+，alaμ;Agρs2 +\t 角ZF'=-Zr;(29)n;'d;where F" is calculated based on the linear averaged(puqp)" sopSai=2+velocity beside the face u;"\tF" =[w,lip中国煤化工"- nd,(Vu), - Ef[(u)n -(u)rup](1-w)usp.FYHCNMH G542F'=U'Ag(30)As an illustration, a one-dimensional problem with auniform grid is used to explain the normalized varia-We can write the velocity correction expressionble diagram (NVD) procedure. Based on the NVD, se-veral total variation diminishing (TVD) schemes canin a similar way as with Eq.(26).be designed. A one-dimensional example is shown inFig.2. The cell face value of a parameter φ，φ, is(u),={9){?'now considered. According to the flow direction onx),the cell face, the values of the same parameter in theupstream and downstream cells are denotedby中andwhereφp, and the value in the further upstream cell (or seco-nd upstream cel) by φsu ，as shown in Fig.3. The nor-(2)= w,+(1-w,)(2malized variable can be defined as .a);a)pa,)(34)Inserting Eqs.(30) and (31) into Eq.(29), we ob-φ -9sutain the pressure correction equation.whichisOatcell SU and 1 at cell D. With the same2) (op'2 (op'Ap=variable value on the downstream and second upwinda),( ax)。a2)。6y ),stream cells, the normalized variable is constant at 1,so we can now just consider the situation near a freesurface where the liquid volume fraction is differenton the downstream and second upwind stream cells.-之r;(32)Therefore, we have at the upstream cell: 中一psuWe rewrite the above equation and insert the ve-(35)locity correction with the pressure correction to themφb-φsuto obtain the pressure correction equationand at the cell interfaceaqpP'=Zap' +s,(33)φ;-su(36)φp-$suFlow directionA，A。_( sa;=|Xg-xqp|(a)。 |Yg-yp|(a2)j中心a=Za，s,=-ZF;j=Fig.2 one-dimensional uniform grid diagram for NVDThe linear equation set, Eq.(33), can be solvediteratively, using the conjugate gradient method. Toaccelerate the solution process, an algebraic multi-gridmethod is adopted. The pressure, the velocity and theflux are updated according to Eqs.(24), (26) and (28),respectively after the pressure correction p' is solved.An algorithm called the compressive interfacecapturing scheme for arbitrary meshes (CICSAM)21Fig.3 Unstructured grid diagram of the upstream cell termsis used for the volume fraction Eq.(3).For the collocated grid finite volume method, allDifferent discrete forms can be expressed by thevariables are stored at the cell centre. The convectionNVD as showI中国煤化工etails can beterm value required on the cell face is obtained throu-found in Ref.[1gh interpolation. There are many schemes available.| YHCNMH G .543in which, β;=(中, -中)/(1-中). The face value can↑UD yUDthen be seen as the usual interpolation using the valuesof the cell centres on both sides of the face, althoughthe weighting cofficient β， depends on the schemeo:used.CDFor a two-dimensional unstructured grid, the linelinking cell centres D and U may not pass throughthe centre of the next cell. The way to choose theFROAMpoint SU is to extend the straight line from DU to05SU, and the distance from U to SU is that fromD to U. The value中sv is then obtained by usingFig.4 NVD diagram for different discrete forms. CD- Centrethe following equationdifference scheme, QUICK- Quadratic upwind interpo-lation for convective kinetics scheme, TUD- Third-orderφsu=φ - 2d.(Vφ)u(41)upwind scheme, FROMM scheme, FUD- First-order up-wind scheme, and SUD- Second-order upwind schemeTo ensure that中su is bounded, we useIn different schemes, different values at the cellcentres are used for SU，U and D to interpolateφsu = Max[Min(Pmx φ3u), φmi]the value at f . When it is normalized, φ is 0 and 1at SU and D. Thus φ, should be a functionof中volume fraction Eq.(3). The gradient of the variationonly, or:is obtained based on an equation similar to Eq.(20).Onno Ubbink (1997) recommended the correla-可=f(q) .(37)tion of the expression for Eq.(11) as follows. This for-mulation is used for the free surface tracking process.再,=克FUD,页=(+h) CD,φ =min |中8c中 +(1-c)(6中u +3)0<元<1,ψ;=号中SUD, ψ,=-+中FROMM,ψ,=克，中<0 or电>1(42),=-(3+64n) QUICK, ,=≤克+→TUDwhere c=unj(t"-t"-)1d, and is the normalcomponent of the velocity on the cell face.Once the value of φ，is obtained from one of theIntegrating Eq.(3) over a cell and using the Gausstheorem for the convection term, we obtainabove schemes, the value of φ，itself can then be ob-tained from Eq.(36) as」, ocds2+f, ,(aU).dA=O(43)Js Atφ=中φ +(1-中,)4s(38)The time derivative term is discretized by using theCrank-Nicolson scheme to obtain a discretized form.From Eq.(35), we also have22,(a"a"-)+ Z(a", +a;")F"=0(44)φsu=_中-中中(39)△=1-中Inserting Eq.(40) into Eq.(44), we obtainInserting this into Eq.(38), we obtaina,a"=a,o"; +s。(45)φj=βrφ+(1- β, )中中国煤化工MHCNMH G544Pressure crrection一Conection[ Calculate fuxu",v",VpCoefficientSIMPLECTransientCNealculationCGSI AMGb=bdr=bUpdate propertiesMoment crrectionρ=ap +(-a)P。PressureL p=p+pμ=a叫+(1-a)luxE-E+Velocity4=码+ujMoment predictionConvectionDiffusionCISccond EulerConverge :implicitGaussVmethod此=bu,v4Fig.5. Flow chart for every time stepwherea°'=βg22.1.3 Solution procedureaop=-Aqpi△t年The above problem is solved by using the follo-wing procedure as shown in Fig.5. First, the volumesa =ax~|2_fraction Eq.(3) is solved. Second, Eqs.(4) and (5) areused to calculate the face density and the viscosity. Fi-nally, the SIMPLE algorithm is used to solve Eqs.(2)and (6).When F" >0, the direction of the fluid flow is .fromcell cp tocell cj,and ap=ag, au=app,so a,=-Bq，aqj=1-βg， a";-=-β， and a")'=1- β;When F;" ≤0, the direction of the fluid flow is7%fromcell cj tocell cp,and ap→ap，a,→aj，so a,=-(1-β,)， a=，al"'=-(1-β;), andFig.6 Sketch of th中国煤化工"TYHCNMH G545the symmetry of the computational domain and theshort period for the impact process, a symmetricalC2C3model is adopted here. The grids are generated byusing the commercial software GAMBIT. Three non-uniform quadrilateral meshes with 7 500, 30 000 andC52 500 cells are used. The time step is adjusted by theB2 |B:/courant number. The shape of the water surface andBAI B.the pressure on the body are given in Fig.8, in whichp' =p/pU2 and s=Ut. As can be seen, differentFig.7 Grid topology for the simulationgrids do not significantly alter the trends of the results.The curve obtained with the grid one deviates fromthat obtained with the other two grids, but the curves2. Simulation modelThe computational domain is shown in Fig.6.obtained with the grids two and three are in good ag-The width and the height are 30 m and 15 m, respec-reement. The computation time seems linearly propor-tively, and the wedge water height is 4 m. The dead-tional to the cell number. Therefore, a mesh similar tothe second grid is now used in this study.，rise anglesare a and β as shown in the figure. Tgenerate a high quality grid, the domain is divided3.2 Comparison with other published resultsinto several sub- domains as in Fig.7. The mesh is refi-The wedge water column impacting on the stru-ned in the sub-domain A, where the wedge water im-cture is extensively found in marine engineering andpacts occur on the structure. In the other sub- domainswas widely studied. Here we use this case to verify thea quadrilateral mesh is generated.present method. The simulation results are shown inFigs.9 and 10. The dimensionless pressure distribu-tions in dimensionless coordinate systems agree wellwith each other for different times and impacting ve-----Grid1locities. This is consistent with the similarity theory of一Grid 3Wwul---- t= 0.5353---t= 0.5555- t= 0.57571.0Duan et a1.[14]0xtZhang et al.15](a) Free surface distribution0.5-1.5[Fig.9 Comparison of the pressure distribution on the body sur-face for different times---- 1.0 m/s-- 1.5 m/sDuan et al[14(b) Dimensionless pressure distributionZhang et al.'sFig.8 Results obtained with different grids0.s-3. Simulation and discussion3.1 Grid-independent studyxs*The case of a symmetric water wedge with thedead-rise angle equal to 60° is used to analyze thFig.10 Compari中国煤化工n the body sur-convergence of results with different grids. Because offace for dif:fYHCNM HG546therefore believed to be due to the approximationbased on an exponential function used for the free sur-一-Presentface profile by Zhang et al.!。Duan et al.142。Zhang et al!15Duan etal.14]Zhang et al!1]5o.5(a) Dead-rise angle of 50° and 70°-1.0.53.1.5[(@) Dead-ise angle of 50f and 70°.0 t2.0空0.5 .a 1.04-20.5(b) Dead-rise angle of40° and 80°-4(b) Dead-rise angle of 40 and 80°3(C) Dead-rise angle of 30° and 90°Fig.11 Wedge water column impact on the structure free sur-face3.3 Wedge water impact on the structuresFig.12 Wedge water column impact on the structure: PressureIn Figs.11 and 12, we show the free surface pro-distributionfile and the pressure distribution on a horizontal wallhit by a liquid wedge defined in the method used in3.4 Simulation with a semi-elliptic water column frontthis study based on the potential theory method fromA water column with a semi -lliptic air cushionDuan et al.[14]and Zhang et a. From Fig.11, it canis now studied by using the numerical methods. Thebe seen that the free surface profiles in this studywidth and the height are 0.1 m and 0.15 m, respec-agree well overally with the data from the literature,tively, while the bottom boundary of the water columnalthough there is some, local discrepancy with the re-is prescribed by the equationsults from Zhang et al."5 especially near or in the jetregion. Good agreements are, achieved between ourstudy and that of Duan et al.', which indicates thaty= 0.03 +0.03//1-0≤x≤0.1( 0.1the nonlinear characteristics affect the free surfaceshape near the jet region. The discrepancy in Fig.12 is中国煤化工MHCNMH G547locity distribution as shown in Fig. 14. The first peak is3「-Elliptic column without cushionat 0.003 s when the bottom right hand corner of the .Elliptic column with cushionwater column approaches the surface of the structure.As shown in Fig.13, the air cushion has an importanteffect on the impact process at this stage. w ithout anyair shear force, the water column remains in a constantshape and the bottom right hand corner hits the struc-ture surface at nearly 0.003 s. A large pressure peakthen occurs while the pressure is much smoother for05the case with an air cushion. The free surface shapesubstantially changes because of the air flow sheart/nsforce in Fig.14. As shown in Fig.13, the case with anelliptic water column front is much more complexFig. 13 Pressure history with an lliptic water column with andthan that with a flat front. The peak pressure is overes-without an air cushiontimated by neglecting the air cushion. The descendingvelocity of the main part of the water column does notchange too much with the air flow and the gravity andthis supports the assumption made by some researcherto ignore the effects of the gravity. However, the va-riation in the water column front shape is crucial tothe impact pressure as shown in Fig.14. The tip hits .the structure surface at 3 ms without an air cushion(ab)and the first pressure peak occurs, while the tip movesa lttle to the right, which causes a delay of the firstpressure peak. The situation is similar for the secondand third pressure peaks.d)4. ConclusionsThis paper presents the development of a finitevolume high order numerical algorithm based on theNVD to simulate the water impact on a rigid structurewith unstructured grids. The basic concept of this al-gorithm is that the free surface can be treated as a phy-sical discontinuity which can be captured by a well-designed high order scheme. Owing to the above me-thods, the water column impact on rigid structures canbe simulated with unstructured fixed grids. The modelin this study is implemented by using an unstructuredgrid finite volume solver.the numerical results with the potential theoretical re-g)(hsults. Both the free surface and the pressure are inFig.14 Snapshot of the velocity and free surface distributiongood agreement.with an elliptical water column at 8 different times. Topwith an air cushion and bottom without an air cushionAcknowledgementThe pressure history at the centre of the structureThis work was supported of the China Scholar-surface is shown in Fig.13. It is shown that the pressu-ship Council (Grant No.2010307245).re on the structure under the action of a semi ellipticbottomed water column is much more complex thanthat for a flat bottomed water column. There are threeReferencesmain peaks during the water column impact processes,1] WOOD D. J.. PEREGRINE D. H. Wave impact on awhile at the same time some spikes also appear beside中国煤化工OorUS hemthe main peaks. These phenomena can be explainedJournal ocean Enginee-by the detail of the free surface movement and the ve-ring, 2000MHCNM HG548[2] HOWISON S. D., OCKENDON J. R. and OLIVER J.[10] REBOUILLAT S., LIKSONOV D. Fluid-structure inte-M. et al. Droplet impact on a thin fluid layer[J]. Jour-raction in partially filled liquid containers: A compara-nal of Fluid Mechanics, 2005, 542: 1-23.tive review of numerical approaches[J]. Computers[3] SMITH F. T. LI L. and WU G. Air cushioning with aand Fluids, 2010, 39(5): 739-746.lubrication/inviscid balance[J]. Journal of Fluid Me-[11] WU G., HU Z. Z. A Taylor series based finite volumechanics, 2003 ,482: 291-318method the Navier- Stokes equations[J]. InternationalHICKS P. D., PURVIS R. Air cushioning and bubbleJournal for NumeriMethods in Fluids, 2008,entrapment:three- dimensional droplet impacts[J].58(12): 1299-1325.Journal of Fluid Mechanics, 2010, 649: 135-163.[12] UBBINK O. Numerical prediction of two fluid systems[5] LAN Ya-mei, GUO Wen-hua and LIU Hua et al. Nume-with sharp interfaces[D]. London, UK: Imperial Collegerical simulation of wave impact on the slab[J]. Journalof Science, Technology and Medicine, 1997.of Hydrodynamics, 2010, 225Suppl.): 986-992.. [13] LEONARD, B. P. The ULTIMATE conservative diffe-[6] DING Zhao-qiang, WANG Guo-yu and REN Bing.rence scheme applied to unsteady one dimensional ad-Three-dimensional numerical simulation of wave slam-vection[J]. Computational method in applied mecha-ming on an open structure[]. Journal of Hydrodyna-nics and engineering, 1991, 88(1): 17-74.mics, 2012, 24(4):526-534.[14] DUAN W., XU G. and WU G. Similarity solution of[7] YANG X., JAMES A. J. Analytic relations for recon-oblique impact of wedge-shaped water column on wed-structing piecewise linear interfaces in triangular andged coastal structures[J]. Coastal I Engineering, 2009,tetrahedral grids[J].Journal of computational physics,56: 400-407.2006, 214(1): 41-54.[15] ZHANGS., YUE D. K. P. and TANIZAWA K. Simula-yulin[8] YANG Wei, LIU Shu-hong andAn unsplittion of plunging wave impact on a vertical wal[J].Lagrangian advection scheme for volume of fluid me-Journal of Fluid Mechanics, 1996, 327: 221-254.thod[J]. Journal of Hydrodynamics, 2010, 22(1): 73-[9] IBRAHIM R., PILIPCHUK V. and IKEDA T. Recentadvances in liquid sloshing dynamics[J]. Applied Me-chanics Review, 2001, 54(2): 133- 199.中国煤化工YHCNMH G