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

Numerical Simulation of High Velocity Water Impact on the Rigid Plate by Discontinuous Galerkin Method

ZHANG Zhong-yü,YAO Xiong-liang,ZHANG A-man

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

Based on Arbitrary Lagrangian-Eulerian method,an inviscid model is established to investigate high velocity water impact problems.In order to capture the pressure wave,the Arbitrary Lagrangian-Eulerian Discontinuous Galerkin method is proposed to simulate the two-dimensional flow of water jet impact on rigid plates.Discontinuous Galerkin method can be able to achieve high order of spatial accuracy and stability.The physics of wave propagation can be obtained exactly by solving the Riemann problem.For fluid domain deformation,radial basis function method is employed to derive the velocities of the internal fluid nodes given the velocities of the boundary nodes,which needs no grid connectivity information.The water jets with flat or round head impact are investigated.For water jet with flat head,the pressure at the center of rigid plate is in good agreement with those by AutoDyn-2D solver.The results show that compressibility plays an important role during the high speed water impact process,as the pressure wave is generated upon impact and propagated into water domain.A larger pressure can also be observed during the water jet with round head impact process.

Discontinuous Galerkin method;Arbitrary Lagrangian-Eulerian; water jet impact

0 Introduction

Water jet impact problems have been paid great attention since they are common phenomena in many fields.There are two stages in the impact process[1].During the first stage,the compressibility of the liquid is dominant and the shock wave is generated by the impact,which can lead to the peak of impact pressure comparable to the‘water hammer’pressure[2],P0=ρcV, where ρ and c are the density and acoustic velocity of the liquid,and V is the impact velocity. It should be noted that the‘water hammer’pressure is based on one-dimensional approximations.However,based on two-dimensional approximations[3],the expression was introduced for the peak of the impact pressure,P0=ρ c+?（ ）V V,where ? is a constant number.After the release of the impact pressure generated by the contact edge of the jet,the next stage begins and the impact pressure decreases rapidly,until approaching the hydrodynamic pressure,when thesteady state is reached.The effect of fluid compressibility on the evolution of the pressure distribution was firstly analyzed with a two dimensional model[4].Then the characteristic of water jet on the rigid plate was analyzed by Korobkin[5],who adopted the eigenfunctions to solve the governing equation,and then the water impact on the elastic plate was further studied[6].The hydrodynamic pressure and the free-surface deformation caused by the impulsive motion of a rigid plate into a strip of incompressible fluid was also studied[7].Strictly speaking,if the interest in impact is over a longer period,the water jet is well developed and there may be a turbulent wall jet region induced by violent free surface motions,even water jet separation from the wall or the splash occurs[8],which is a complex problem.However,the present study focuses on the initial stage and the compressibility effect.

Numerical simulations of the water jet impact have been the subject of intense research. In the present work,the Arbitrary Lagrangian-Eulerian(ALE)Discontinuous Galerkin method was used to solve this problem.The Discontinuous Galerkin method was firstly proposed by Reed and Hill[9],which has become popular for the solution of various problems.The Discontinuous Galerkin method[10-11]combines two advantages associated to finite element and finite volume methods.One is that the spatial accuracy is improved by means of high-order polynomial approximation within an element as in finite element methods.The other is that the discontinuous solution at element interfaces is assumed to obtain the physics of wave propagation exactly by solving the Riemann problem.Then,the Discontinuous Galerkin method can be able to achieve high order of spatial accuracy,compactness,stability and highly parallel computation.Water jet impact on rigid plates is essentially a moving boundary problem.The ALE method is most commonly used in flow problems involving fluid domain deformation.The Discontinuous Galerkin method in framework of ALE was firstly introduced for compressible Navier-Stokes flows in moving domains[12].The geometric conservation law was always satisfied without introducing a grid velocity source term[13].A high-order Discontinuous Galerkin method on deformable domains was proposed by constructing with mapping between a fixed reference configuration and a time-varying domain[14].An additional equation related to the Jacobian of transformations was introduced to ensure the preservation of constant solutions.

In this work,the goal is to simulate the flow of water jet impact on rigid plates by the ALE Discontinuous Galerkin method.The governing equation is discretized using Discontinuous Galerkin spatial discretization.The deformations of fluid domain are simulated in ALE form,and radial basis functions method is used to interpolate velocities of the boundary nodes to the whole meshes.Simulations of water jets with flat or round head impact on the rigid plates are presented,which aim at providing some reference for the design of marine engineering structures.

1 Mathematical model

1.1 Governing equations

We consider the simulation of water jet impact on the rigid plate as a two-dimensionalunsteady flow for simplicity.The water domain Ω is formed with the width D and the length L,and Γ denotes the boundary of domain. Initially the water moves with the high velocity V to the rigid plate,which is supported rigidly,as shown in Fig.1.

Assume that the water is inviscid and compressible,and the gravity and surface tensor are ignored,we consider the Euler equations in the water domain which are written as a system of conservation laws,

where ρ,p and E denote the density,pressure and specific total energy of the fluid,respectively,and uiis the velocity component of the flow in the coordinate direction i.The Tait EOS is added for the completeness of the governing equations,

For the relation p=p ρ,（ ）e given by the Eq.(3),we can have the acoustic speed c,which is a useful state variable for the unsteady flow of compressible water jet impact on rigid plates,

where peand pρdenote the partial derivatives of pressure p with respect to internal energy per unit mass e and the density ρ.

1.2 ALE formulation

It is obvious that boundaries of water domain deform during the water jet process,and the domain Ω changes with time t,in other words,the motion of mesh must be accounted for in the simulation.Then ALE formulation is adopted,as the method allows us to move and deform the domain.The space discretization of Eq.(1)is carried out with the Discontinuous Galerkin method.

Let us divide the domain Ω（）tinto a collection of nonoverlapping triangular meshes Τ=Thus,the test functions space is defined as

Multiplying Eq.(1)by an arbitrary test function φi,integrating over the each element,integrating by parts the divergence term and substituting Eq.(6)into Eq.(1),we have the weak form of Eq.(1)

where n denotes the unit outward normal vector to the boundary,?K is the boundary of the element.The argumentshave been dropped for the sake of simplicity.

Because the domain varies with time for moving mesh problems,the time derivative in the first term is taken outside the integral to yield the relation in the ALE domain[13]

here,we define w is the grid velocity.

Substituting Eq.(9)into Eq.(8),we have

It can also be obtained for test function φiby Reynolds transport theorem

where?K0represents the internal boundaries and?Kbdenotes the boundaries of water domain.Thenotation indicate the trace of the solution for the interior element and the neighboring element,whiledenotes the trace of those for boundary element.is the Riemann flux function,which is solved by the HLLC approximate Riemann solver in the present work.This HLLC scheme has been successfully used to compute compressible flow on unstructured grids in the ALE domain,which are also found to preserve isolated contact and shear waves exactly.

By assembling together all the elemental contributions,a system of ordinary differential equations governing the evolution in time of the discrete solution can be written as

It should be note that the mass matrix must be recomputed in each step time as the test function changes with the time.In the present work,the test function is chosen as a set of polynomials base on Taylor series expansions at the centroid of element[11].The Eq.(13)is discretized in time by a two-stage Runge-Kutta method.The time step is restricted by the CFL condition of stability,

where Δl is the diameter of each element K,λ1and λ2are the maximum wave velocities in the x-and y-directions,respectively.The GCL condition for Discontinuous Galerkin method was proved mathematically that the ALE formulation satisfies the uniform flow exactly[13].

1.3 Radial basis function method

Radial basis function method[15]can be used to derive the velocity of the internal fluid nodes given the velocity of the boundary nodes,as shown in Fig.2.It has been proven that this method works well even for large mesh deformations,both for 2D and 3D meshes.This method of deforming the mesh is straightforward to implement with a small system of equations for the boundary nodes to be solved.The other superiority is that no grid-connectivity information is needed.

Fig.2 Boundary and internal nodes

The grid velocity function χ（）

X has the form by a sum of basis functions,

where nn is the number of boundary nodes,Xiis the location of boundary nodes.is apolynomial and H is the compactly supported C2function which is respect to the Euclidean distance,

where r is the scaled radius to control the mesh movement.

Then the coefficients αiare determined by the velocities of boundary nodes

and the additionally requiring for the polynomials q（）X of degree less or equal than that of function H

This gives a linear system

With α is the vector of coefficients αiand β contains coefficients of polynomialG is a nn×nn matrix and each unit of which has to compute basis functionsat each boundary node.B is a nn×3 matrix with row j given bySolve the above system and the velocities of the interior nodes can then be derived by evaluating Eq.(16)in the internal nodes.

1.4 Boundary conditions

In the present work,all boundary conditions will be imposed in a weak manner.At the free surface,the pressure is given as atmospheric pressure Pa,and the other variables are equal to the interior element.At the rigid plate,referred to as inviscid walls,flow tangency is enforced to be zero.The density ρband total energy Ebat the integration point are equal to the interior element.Then the boundary state vector at the integration point is computed as

where niis the unit normal to the rigid plate.

2 Results and discussions

2.1 Numerical verification

The case of water jet with flat head impact on the rigid plate is considered here.The width of fluid domain is 2 m and the length is 10 m.The density of water is 1 000 kg/m3and the impact velocity is 100 m/s.The convergence study to this problem is firstly computed on threesets of triangular meshes with 500 elements,3 000 elements and 10 000 elements,respectively.The second order of polynomials is used for Discontinuous Galerkin method in the present work. The pressure histories at the center of rigid plate with different meshes are shown in Fig.3.It can be seen that different meshes do not significantly alter the trends of the results,and the curves obtained with 3 000 elements and 10 000 elements are in good agreement,which means that the results have converged well.In the following sections,the results for water impact with flat head are obtained on the 10 000 triangular elements.

To further validate the present method,the same case is also solved by Autodyn-2D solver[16],and the fluid domain is divided into the 100 000 quadrilateral elements.The SNL model is chosen for water.Fig.4 shows the comparisons of the histories of pressure at the center of the rigid plate obtained by the present method and Autodyn-2D solver and the agreement is quite good.It should note that the maximum pressure is higher than those obtained by the‘water hammer’pressure,or p=ρcV,which is based on one-dimensional approximations. Furthermore,the acoustic velocity c is assumed to be constant as only a small disturbance may exist with the lower impact velocity.However,the acoustic velocity can be up to the 1 800 m/s as the compressibility of water during the high speed impact process,as shown in Fig.5.The same phenomenon has also been founded in the droplet impact on the surface[17].

Fig.4 Histories of pressure at the center by present method and Autodyn-2D solver

Fig.5 Contour of acoustic velocity

2.2 Simulations for water impact with flat head

Fig.6 gives the evolution of shock wave propagating in the water,which is generated by water jet impact with flat head.The case of water jet with flat head impact on the rigid plate is considered here.The width of fluid domain is 2 m and the length is 10 m.The density of wateris 1 000 kg/m3and the impact velocity is 100 m/s.Initially the high speed water moves suddenly onto the rigid plate,and the shock wave is generated only at the area adjacent to the rigid plate.In other words,the fluid away from the rigid plate is not compressed,as shown in Fig.6(a).The compressibility of water is of major significance and the high pressure can be obtained with the high speed water impact.Then the shock wave travels downward in the water, at the same time,the reflection waves are generated by the free surface of water domain,which propagate to the center of the water domain and make the pressure adjacent to the free surface lower.Moreover,the free surface can move along the rigid plate simultaneously,as shown in Fig.6(b).It can be seenfrom Fig.6(c)that the reflection wave has been propagated at center of the rigid plate and makes the pressure at the center decreased,which can also be obtained from pressure curve at the center in Fig.4.At time t=1.1 ms,the shock wave has been departed from the rigid plate completely.

Fig.6 Contours of shock wave generated by water impact with flat head at different time

It is obvious that the free surface moves along the rigid plate,though we only concern the water impact for a short time.Thefree surface moves along the rigid plate consistently.Because of the reflection wave by the free surface,the pulse width of pressure at the center of the rigid plate is comparable to the distance between the center and free surface against acoustic speed c,which can be proved by the pressure curve in Fig.4.It should note that the gravity is ignored during the simulation.However,the gravity may have an influence on the free surface deformation for water jet impact lasting for a long period.

Fig.7 shows pressure distribution along the rigid plate by water impact with flat head. Owing to shock wave generated adjacent to the plate,at time t=0.1 ms,the rigid plate is subjected to large impact pressure,up to 180 MPa,which may cause structure damage.When the reflection wave by free surface propagates to the center of the water domain at time t=0.5 ms,pressure curve becomes the anti-U shape with the interaction between the shock wave and the reflection wave.Then the pressure along the rigid plate gets lower,while the shock wave has fade away from the rigid plate,as shown in Fig.6(c) and 6(d).

2.3 Simulations for water impact with round head

Fig.7 Pressure distribution along the rigid plate by water impact with flat head

The case of water jet impact with round head is also considered.The width of fluid domain is 2 m,and the radius of round head is 1 m.The length is 10 m,the density of water is 1 000 kg/m3and the impact velocity is 100 m/s.The fluid domain is divided into a set of 9 200 triangular elements with 4 721 nodes.The second order of polynomials is used for Discontinuous Galerkin method to solve the problem.Initially,the round head of water domain is tangent to the rigid plate.

Fig.8(a)-(f)shows the evolution of shock wave propagating in the water during water jet impact process.Unlike those of water jet with flat head,the contact area will become larger gradually as the round head arrive at the rigid plate in sequence,which is similar to the phenomenon of water droplet impact on surface.Due to the compressibility of water,a creation of a strong shock wave can be observed(Fig.8(a)).The water adjacent to the contact rigid plate is highly compressed,while the rest of water domain remains unaware of the impact.The two regions are separated by a shock front which travels into water(Fig.9).The intersection between round head and rigid plate is named as‘contact edge’.Then refection waves sweep into the water domain and the shock begins to detach from the contact plate,as shown in Fig.8(b) and 8(c).At the last,the departure of shock wave can be seen in Fig.8(d)-(f),whereas a new shock wave is generated adjacent to the contact rigid plate as the impact lasts in sequence.

Fig.8 Contours of shock wave generated by water impact with round head at different time

The velocity of contact edge Uehas been given in theory for the problem of water droplet impact on surface[17], The contact edge velocity initially has an infinite value by the above equation,remains higher than the shock speed throughout this initial stage.So the shock front should be attached to the contact periphery.After some time,the shock wave velocity is greater than those of contact edge,and the shock wave begins to travel along the free surface.

At time t=0.03 ms,the center of rigid plate is subjected to an acoustic pressure by the generated shock wave,as shown in Fig.10.The contact area will continue to generate secondary shock waves.These shock waves are superimposed on the previously generated shock wave and the rigid plate is subjected to larger impact pressure rather than acoustic pressure.Meanwhile,there is no reflection wave making the pressure lower as the contact edge velocity is faster than shock wave velocity.We can observe that the maximum pressure has been reached up to 234 MPa at the time t=0.1 ms in Fig.10,more beyond to the‘water hammer’pressure.The maximum pressure is also call as contact pressure in the problem of water droplet,which has been studied experimentally.However,when the shock wave begins to travel along the free surface,which can generate the reflection waves,the maximum pressure along the rigid plate decreases rapidly at time t≥0.2 ms.

Fig.9 Sketch of shock front and contact edge

Fig.10 Pressure distribution along the rigid plate by water jet impact with round head

3 Conclusions

The ALE Discontinuous Galerkin method is proposed to investigate the water impact on rigid plates computationally.It has been demonstrated that compressibility of fluid dominates during the early impact process.The simulations show that a shock wave attached to the rigid plate is generated upon impact.Then the shock wave travels downward in the water,while thereflection waves,generated by the free surface,propagate to the center of fluid domain.The maximum pressure has been reached up to 234 MPa during water jet with round head impact process,more beyond to the‘water hammer’pressure,which may cause structure damage.

[1]Hsu C Y,Liang C C,Teng T L,Nguyen A T.A numerical study on high-speed water jet impact[J].Ocean Engineering, 2013,72(1):98-106.

[2]Cook S S.Erosion of water-hammer[J].Pro.R.Soc.Lond.Ser.,1928,119:418-488.

[3]Heymann F J.High-Speed Impact between a liquid drop and a solid surface[J].Journal of Applied Physics,1969,40(13): 5113-5122.

[4]Frankel I.Compressible flow induced by the transient motion of a wavemaker[J].Journal of Applied Mathematics and Physics,1990,4(5):628-655.

[5]Korobkin A A.Global characteristics of jet impact[J].Journal of Fluid Mechanics,1996,307:63-84.

[6]Korobkin A A,Khabakhpasheva T I,Wu G X.Coupled hydrodynamic and structural analysis of compressible jet impact onto elastic panels[J].Journal of Fluids and Structures,2008,24(7):1021-1041.

[7]Needham D J,Billingham J,King A C.The initial development of a jet caused by fluid body and free-surface interaction [J].Journal of Fluid Mechanics,2007,578:67-84.

[8]Chen Z,Xiao X,Wang D Y.Simulation of flat bottom structure slamming[J].China Ocean Engineering,2005,19(3):385-394.

[9]Reed W H,Hill T R.Triangular mesh methods for the neutron transport equation[R].Los Alamos Report,LA-UR-73-479.

[10]Cockburn B,Shu C W.The Runge-Kutta discontinuous Galerkin method for conservation laws V:Multidimensional systems[J].Journal of Computational Physics,1998,141(2):199-224.

[11]Luo H,Baum J D,L?hner R.A discontinuous Galerkin method based on a Taylor basis for the compressible flows on arbitrary grids[J].Journal of Computational Physics,2008,227(20):8875-8893.

[12]Lomtev I,Kirby R M,Karniadakis G E.A discontinuous Galerkin ALE method for compressible viscous flows in moving domains[J].Journal of Computational Physics,1999,155(1):128-159.

[13]Nguyen V T.An arbitrary Lagrangian-Eulerian discontinuous Galerkin method for simulations of flows over variable geometries[J].Journal of Fluids and Structures,2010,26(2):312-329.

[14]Persson P O,Bonet J,Peraire J.Discontinuous Galerkin solution of the Navier-Stokes equations on deformable domains [J].Computer Methods in Applied Mechanics and Engineering,2009,198(7):1585-1595.

[15]Beckert A,Wendland H.Multivariate interpolation for fluid-structure-interaction problems using radial basis functions[J]. Aerospace Science and Technology,2001,5(2):125-134.

[16]ANSYS,ANSYS/Autodyn,User Documentation.Version 6.1[K].2007.

[17]Haller K K,Ventikos Y,Poulikakos D,Monkewitz P.Computational study of high-speed liquid droplet impact[J].Journal of Applied Physics,2002,92(5):2821-2828.

基于間斷有限元方法的可壓縮水高速?zèng)_擊平板特性研究

張忠宇，姚熊亮，張阿漫

（哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001）

基于任意歐拉拉格朗日（ALE）方法，忽略流體的粘性，建立水射流沖擊剛性平板的數(shù)值模型。為了更為精確地捕捉流場(chǎng)中壓力波，該文推導(dǎo)了在ALE框架下的間斷有限元方法。該方法易于提高數(shù)值離散的空間精度，數(shù)值穩(wěn)定性較好，利于精確模擬高速水射流沖擊過程。對(duì)于自由液面的變形，采用徑向基函數(shù)方法確定網(wǎng)格單元的速度。該方法采用邊界網(wǎng)格節(jié)點(diǎn)的信息去推導(dǎo)出內(nèi)部網(wǎng)格節(jié)點(diǎn)的信息，不需要單元信息。平頭水柱射流沖擊的數(shù)值結(jié)果與Autodyn的數(shù)值結(jié)果吻合較好。在驗(yàn)證方法合理性的基礎(chǔ)上，文中對(duì)平頭及圓頭水柱中壓力波的分布特性進(jìn)行了分析。數(shù)值結(jié)果表明：當(dāng)沖擊在平板表面時(shí)，會(huì)在產(chǎn)生一個(gè)壓力波之后向水域內(nèi)傳播，且使得平板受到較大壓力的作用。此外，圓頭水柱的射流沖擊將會(huì)產(chǎn)生一個(gè)更大的壓力峰值。

間斷有限元方法；ALE方法；水射流沖擊

O353.4

A

張忠宇（1988-），男，哈爾濱工程大學(xué)博士研究生；

O353.4

A

10.3969/j.issn.1007-7294.2016.06.004

1007-7294（2016）06-0674-12

姚熊亮（1963-），男，哈爾濱工程大學(xué)教授，博士生導(dǎo)師；

張阿漫（1981-），男，哈爾濱工程大學(xué)教授，博士生導(dǎo)師。

Received date：2016-02-11

Biography:ZHANG Zhong-yu(1988-),male,Ph.D.student in Harbin Engineering University,

E-mail:zhangyu061031@126.com;YAO Xiong-liang(1963-),male,professor/tutor.