Benefits Of Using A Wendland Kernel For Free-surface Flows - CORE

Transcription

Benefits of using a Wendland kernel for free-surfaceflowsE AlaciaA. SoutO'IglesiasNaval Architecture Dept.(ETSIN)Tfeclinical University of Madrid (UPM)28040 Madrid, Spainfabricio.macia,antonio.souto(ifr upm.esM. AntnonoINS BAN00128 Rome, Italy.matteoantouno@gntail.comA ColagrossiINSEAN0012S Rome, ItalyCeSOS, NTNU,Trondheim, Norwaya .col a grossí @ insean. ilAbstract—The aim of this paper Is lo discuss the influenceof the selection of the interpolation kernel in the accuracy ofthe modeling of the internal viscous dissipation in Tree surfaceHows, Simulations corresponding to a standing wave* for whichan analytic solution 1 available, are presented. Wendland andrenormalized Gaussian kernels are considered. The differencesin the flow pattern* and Internal dissipation mechanisms aredocumented for a range of Reynolds numbers. It is shown thatthe simulations with Wendland kernels replicate the dissipationmechanisms more accurately than those with a renormalizedGaussian kernel. Although some explanations are h i n t e d wehave Tailed to clarify which the core structural reasons for Michdifferences are*aim of the present work is to document the benefits of usingsuch a kernel when dealing with viscous free surface flows*The paper is organized as follows; first, the physicalproblem we are interested in is presented* Second, the SPHmodel, emphasizing the properties of the considered kernels isintroduced. Third, Ihe practical problem considered, which hasbeen the dissipation of the kinetic energy in a standing wave.is introduced and the results of the computations carried out,for a range of Reynolds numbers, are presented and discussed*Finally, some conclusions are drawn and future work threadshinted.SYMBOLS11- GOVERNING EQUATIONSFree surface 2D laminar Newtonian incompressible flowsare treated in this papen The viscous effects in these flowsare proportional to the Laplacian of the velocity field, a.The equations that describe these flows are the Navier-Stokes1- INTRODUCTIONIn Colagrossi el al. [1J. a theoretical analysis of the two incompressible ones.In order to close this system of equations it is necessarymost popular SPH formulalions for Newtonian viscous lerms,i.e. Monaghan's [2] and Morris et a/.'s [3], was carried out at to specify the boundary conditions (BC). Let the fluidthe continuous level. It was shown that in tlie presence of a domain be noted as il and its boundary as rííi Suchfree surface and under certain conditions, these viscous terms boundary encompasses afreesurface boundary dilp and solidbecome singular with die inverse of the smoothing length. It boundaries Oils* Along the free surface, both a kinematic andwas also demonstrated that for the Monaghan's [2] viscous a dynamic BC should be in principle fulfilled. They are notterm Ihe singularity does not affect the integralflowquantities* explicitly included in the SPH simulation, since the weaklyAs a consequence, the exact mechanical energy dissipation rate compressible model is considered and such model can beis recovered in the continuous approximation as the smoothing properly shaped in order to being inherently consistent withthe kinematic and pressure conditions [8] and simultaneouslylength goes to zero.providing the correct viscous dissipation [1]. NotwithstandingWliile in practical SPH implementations for pure diffusion that, the dynamic free surface BC is presented (equation 1 forprocesses, the viscous terms behave similarly [4], in order to the normal stress and 2 for the tangent stress) both for the sakeassess the practical implications of the aforementioned results, of completeness and in order to recall that when viscosity isa viscous free surface flow is simulated with SPH, Such relevant, pressure is not in principle zero at the interface andflow is the evolution of a standing wave. This is a classical is in general discontinuos.problem in tl»e scientific lilerature and is of practical intereslsince it is strictly related to the propagation of gravity waves2///i duidn(1)Moreover, an analytic solution for the decay of the kineticD « 0(2)energy is available (see Lighthill [5]). The SPH simulationshave been implemented using free-slip conditions along the In these equations n is the normal vector to the fluid domain,solid boundary of the tank and using a renormalized Gaussian* r is an unitary vector lying on the free surface tangent planeKernel [6]. and afifthorder class 2 Wendland [7] kernel The and D is the rate of strain tensor.WC2RGK30Wendland 5th degree class 2 kernelRenormalized Gaussian kernel

6 international SPHERIC workshopHamburg, Germany. June, 08-10 2011A free slip BC is imposed on the solid boundaries since weare interested in assessing the internal dissipation mechanismsand nol those related to the solid boundaries.11¿rÍ5.1-HI. SPH MODEL[1AA. GeneralA weakly compressible approach is considered for thesimulations [9]* The Monaghan's [2] viscous term is used,which, as mentioned in the introduction, allows in principlethe computation of the exact mechanical energy dissipationrate since it is recovered in the continuous approximation asthe smoothing length goes to zero,11;.,Jt 1.-. - v— R g n o m d l i z e d Gaussian Kernel"- Rr-:;oi»iiAi -ttd &useian K e r n e l 2nd d e r i v a t i v eHendland Kernel-n—Wfcooiendmr— m ? aodderivative1i1 —n'*' JiA0,20,406 0L2141651eFig* I* Wendland and Gaussian Ke/ndFB. KernelsThe kernels considered are the 5* degree class 2 Wendlandkernel [7], [10] and a renorrnalized Gaussian Kernel [4], [6].They are both isotropic, which means that they can be written*te* ¿*Mib**wIS** S . TS¿ -;;;. *# »*., ¿i*¿ 1,.,j: plg(3)- ¿¿!tí"?: /rfwhere x is a vector in R , h is die smoothing length, d is I hedimensionality of the problem, and q \x\fh. They both have'Reaomftlirtd Gau&staii Keen?! d e r i v a t i v ea compact supportiWendland [71 devised a family of radial interpolation1-Z14l*S1*1*t*functions that have compact support and that are definitepositive. We have used the 5th degree class 2 Wendland kernel Fig. 2. Fourier iransfonus, I lie noitnalized wave number ¡ is scaledfioIliac(WC2 from now on) reformulating it so that it has a 2// radius k 1 for the mode with wavelength X hsupport, similarly to Robinson et al- [10], [11].PilpiÉl" Mñ-&(2-g),{l 2í)0for0 9 2for q 2(4)where d is the dimensionality of the kernel and /? is a constantthai depends on d and thai is equal to 7/(64/0 in 2D, whichis how die simulations have been carried out.Together with the Wendland kernel, the renorrnalizedGaussian kernel (RGK from now on) of Molteni andColagrossi [6] has bee also considered According to [6], froma numerical point of view the behavior of the RGK is almostidentical to the classical Gaussian kernel and the followingproperties are well established: (i) amongst IÜ tested kernelshapes, the Gaussian kernel appears to give the best numericalaccuracy in the stable field [12]; (ii) the comparison of theGaussian kernel to classically used spline kernels showedthat the former leads to better stability properties [13]; (¡ii)it presents also a lower computational cost with respect toevolved forms of spline kernels [14];finally*(iv) its gradientcan be straightforwardly obtained from the evaluation of Ihekernel itself.Since the radius of Molteni and Colagrossi RGK [6] is3/J, we have reformulated it (equation 5) so that its radiusis 2/f and thus be possible to compare its performance withthe Wendland kernel one, using Ihe same parameters:exp(-c*2 f 2 ) - m» for 0 q 2WRGKW - Hfor q 2 ! 0(5) :Where a 3/2, and in 2D A/3/2 0.717082 and n% exp(-9). A plot of both kernels and their second derivativesis presented in figure I.Tensile instability may arise, according to [15], when aparticle enters in the negative second derivative region of theneighborhood of another particle and the pressure is positive.As can be appreciated in figure 1, the size of these regionsis very similar for both kernels {q 0.5 for WVc? andq N/2/3 * 0.4714 for fi *). which suggests they shouldbehave similarly in regards to the onset of such instability*We will see this is not the case- Monis [131 showed that theFourier transform of the kernels is relevant for the stability ofSPH simulations. He found that stability is improved if thetransforms fall off rapidly. If we look at the spectrum of bothkernels (figure 2), we find that the RGK spectral componentsare smaller than WC2*s for wave numbers greater than thecorresponding to h. This would suggest that the RGK shouldfilter the high frequency oscillations better than the WC3, Wewill see that this is not the case, Robinson [II] got to asimilarconclusion, when comparing WC2 lo spline kernels,TV, RESUUSA. GeneralSimulations corresponding to the dissipation of kineticenergy in a viscous standing wave have been conducted- There31

6th international SPHERIC workshopHamburg, Germany. June, 08-10 2011is an analytic solulion by Lighthill [5] lo compare to that isalso presented. A range of Reynolds numbers, resolutions andnumber of neighbors has been covered in the simulations. Theresults are doeumenled in this section, showing the diverseperformance of ihe RGK and WC2 kernels.where v is the kinematic viscosity of the fluid The dampingcoefficient is /Í 4 kJP and depends on the wave number andon the kinematic viscosity- Lighthill [5] obtained Ihis result bysubstituting the linear solution (6) into the dissipative term ofthe kinetic energy, that is:g J(D:B)áV J . J ¿ 2 & p%)dV,B. Dissipative effects for Standing WavesIn Ihe following, H denotes the still water depth, L is Ihetank length, A the wave length of the standing wave (which and, then, averaging over a wave period- In this context we arewill be taken as L in our simulations} and k the corresponding neglecting Ihe influence of solid boundaries and of the bottomwave number (Le. k 2jr/-t). The standing wave amplitude friction.In the work of Wu et al. {16J, the linearized Navier-Stokesis indicated by A while denotes the ratio 2A/H- For smallamplitude waves (i.e small e), and small wave steepness equations are solved to inspect the influence of viscosity on{Le 2AfA 0.1) the potential theory gives the following free-surface waves evolving into a two-dimensional lank- Toapproximate analytical solution (see figure 1V-B for a sketch justify the use of the linearized equaticais, Wu et al. [16]assumed e 1 similarly to what was done lo get the solutionof the flow domain and initial velocity vectors)):(6j. Further, they imposed a free-slip condition along thesolid boundaries. Under these hypotheses, the evolution of the fix y.O A (* vj cos(wi)potential energy for the standing wave is:Hg cosh[¿(y tf)] fo(x,y) -ecoa(kx)2 a* 42wcoslniWi*-atCOS5*'{&') t — —sRe* Re2Re,-2««ca .(s[gki)M 1 Rrj110)where Re H \§JJfv is Ihe Reynolds number and a 2 vk .Under Ihe hypotheses:-Re-d.2 -and S*'ReKkny*MS)-di)Ihe potential energy in (10) behaves as:-a.4-M : ; Re«resin( gklllllil!;;;:iEM)!EM * -2ff(SO-(12)This expression shows that, for short times, the dampingrate of gravity waves is /Í 2 x 4vA?. Such a dampingcoefficient coincides with that found by Lighthill [5).-i .i * I0.2 a-*UiULLC. Stnudariom setupThe initial arrangement of the particles is a cartesian latticeFig. 3. Sketch of ite standing wave problem and initial velocity Held(figure 4). Tentative computations with other arrangementssuggest that the conclusions do not depend on the arrangementThe period of oscillation of the standing waves is T 2JT/6 type. Ghost particles (see e.g. [9]) are used to enforce thewhere the circular frequency & is given by Ihe dispersion free slip boundary condition. Simulations have been carriedrelation: u? gk tanh(A-//) and g is the gravity acceleration. out for both kernels covering a range of Reynolds numbers,We set the time origin r 0 when the free surface is horizontal resolutions and in some particular cases different dx/h ratios.and the time derivative of the velocity potential p is zero in Ihe A summary of the simulations parameters is presented inwhole domain. In this way, tlie pressure field can be simply table L The trends observed for Re 500 are maintainedassumed to be hydrostatic with an error of 0(e*}. By definition for Re 1000,2000 and relevant changes appear for Re tlie initial fluid velocity is given by V . For an inviscid 5000- Therefore, due to lack of space, the results for thefluid, the linear theory predicts that the standing wave evolves intermediate Reynolds number cases are skipped in Ihe presentinfinitely without any change of amplitude, i.e* Ihe total energy paper. The simulations have been conducted assuming theis conserved. For viscous fluid an analytical approximation of regime is laminar. This assumption is valid for ihe range ofthe kinetic energy dissipation can be derived following the Reynolds numbers considered, as shown in reference [17],procedure described by Lighthill [5J. This reads:where transitional flows were described for Re 20000Moreover. free slip boundary conditions have been used, which2AHprevent the generation of vortieity close to the walls- In orderft» 0 g—-e-4yk3'[\ cos(2«/)](8)32**iJL(9

6* international SPHERIC workshopHantaig. Germany, lune. 08 10 2011IU-S00; Hto 4O0— ,.«-*/ *» '*;'Fia. 4rInitio) -:ú&Jtti.ia htiice cí tfcr pitiicíei. i U ths epical [utii-fcTABLE IOSentíWCi RGKTfrc:fa" W)ROEi*:.\-WCJ: &JT W'CIKOKft?.**MX Íí*3. A W-vr 7tuÜ'.'. J " ' .*'.ll1W. i*Xi*3iiIOVÍC&4ÍO.ECOJTTiH». Xli. W . t f nÍ5tn* ST*"":;.: ". M/h 05, Rf SO» fffdx Mti 35bocoukf explain the gap in the kinetic energy decay rate curvesobserved In figure 5. TWs dissipation of energy throughnumerical fluctuations was analyzed in other context by Elleroi/.' W*: "4!" ''et al [211, by l o w i n g the analogy with the dissipationT I ft /¿*. i 9nvchanism at mtfecular level.Observing figure 7(b It seems that the instabilities aresonchowrelated to the effect of the bottom boundary* butto check that the decay rates [resented in section IV tí arcthatisnotclear yet at this stage. In addition, t t e stroeture ofvalid for Reynolds nurrters oí the sams order of Ihe onestheRGKkernelvorticity tick! aiggests that tensile instabilitythat will be used in our confutations, simulations with amayhavebeenexcited wlule that is not ths case for thewell established VOP solver. STAR-CCM-f [18). haw beenWC2kernelThisis unexpected since die pressure field is aconducted*compressionowdueto the effect of the hyoVotfatic pressureIn r ards to specific detail* oí rte SPH ím ltínwntationandfordi/h 0-5particlesare not into the negative secondfor this problem (initial condition for Ihe pressure, free slipderivativeregionofthekernel.Robinson 11] showed inboundary conditions* H&) we refer to reference 1191forced turbulence smulatioiis that the WC2 kemd preventedthe clustering of particles, keeping them at relative distances/ . Results for Re 500clos? to the Initial lattice- In our particular case-, ihe initialResults for the decay of ihe Kintflc Energy ai predicted distribution of the particles can be described by figure 8 whereby using tfie Mooaghan's (21 viscosity fotmula and UghtfitlPs *K is a histogram of psvticle pairs ditfarcres for parities Inside[5] reference solution (dashed Mack line) ace presented fur the stfport of the kernel and K . is the radial partirle densityRGK and WC2 kenwls in figure 5 , The curves correspond to function as defined by Robinson- [11]. Tills function is athe highest resolution Hid* 400 and to dxfh 0.5. which good indicator of particte clustering. It is plotted for bothmeans that each particle has around 50 neiglibors* As can be kernels In a later cycle in figura and the differences canappreciated, die agreement with Ihe exact solution is good. be clearly appreciated; the WC2 kernel keeps t t e relativethough SPH tends to underestimate the decay rate, which Is distawe between particles whilst the RGK ones [resent amore evident for the WC2 simulation.significant scatter but not a clustering of neighborly particles*Tte vorticity fieH Is now a n a l y w l W#i*city is conpited by I b i s scatter induces extra inaecurackrs in the evaluation of thedifferentiating the velocity field using MLS correction for the SPH operators» as d t i c i « » l by 1221, 1231 but Its origin andSPH interpolation I2QJ. Remits are presented for one Instant link with the structure of the kernel* is not clear at this stage of the first oscillation cycle in figure 6, This figure shows theonset of the viscous boundary layer at the free-surface dueIn order to check whether the instability Is betng excitedto the unbalance of the zero stress boundary condition. Both for the RGK. the ratio dxfh was mxfcr smaller from 0.5 tokernels f rfonn yet similarly at this stage*#9 0.222. which rwans passing from50 to250 neighbors* soIf the simulation is run for a few more periods, the that at least there are two rings of neighbors inside the negativedirterences between both kernels results can be clearly kernel second derivative region. whose bounds have beenappreciated. as the vorticity field of the RGK one becomes shown in section Hi H. On one side* the noise in the vorticitynoisy (figure 7 . This means thai a certain amount of kinetic rtekl Is significantly reduced» the kinetic erargy increases withenergy is being lost in these numerical fluctuations* which respect to thr resolution d\ih 03 (less spurious dissipation): : wm* r TT55u? "IK I

6th international SPHERIC workshopHamburg, Germany. June, 08-10 2011ioiiu)iWg)rot(u) i t t g f "I I I I III B« I d f f / t f « t t t to* ) « ?Ce *ÍI "L\ .H« 8 « i - 361 1001 'D -0Í -0-1 *** MOS11*)X *H00204DO05(2t WC2 a) WC2HXtUKH Fu 'l 4 4 ' 1 J 4 1Re 503H/ih CO1Í -TVorticiry. & 500. dx/h 0.5. HgfHfi*2xH(b) RGK 0.361, W/rfx 4 0 0and gels ticker IO the one provided by W C 2 kernel for thatresolution (figure 10), which presents less dissipation than iheexact solution. In regards to the radial particle density function,it shows scattering from the initial configuration (figure 11).in a similar manner as in the case with dx/h 0-5- We wouldexpect to llave less scattering due to the vorticily field beingsmoother but it seems this is not the case. Also we wouldexpect the triggering of Ihe tensile instability because there areparticles with positive pressure in an area of negative secondderivative of Ihe kernel, but it seems no clustering of neighborsis taking place. . Results for Re 5000- W -\H(b) RGKFig. 6.I x'HFig. 7. Vorticily. fi* 500, Jx/H 05, Kg/ftf* 2.SS0. ///*/* 400,.*!* -i111rL L lI J* *V* i"i iFig. 8. Particle density graph for t{g!H?5 0, H¡dx 400.The simulations with the WC2 kernel are now moreThe kinetic energy evolution is shown in figure i 2. It is dissipative than the exact solution. This could be due to thenoticeable that due to the reduction in viscosity there is still spurious vorticity that is created mainly in the lateral walls,60% of Ihe kinetic energy alter 20 oscillations, compared to for which we have found yet no explanation (figure 14(a)).less than X% for the Re 500 case (figure 5). As the Reynolds The noise in the RGK simulations is now even more intensenumber gets higher, the gradients of vorticity in the region and the scales of the noise are spread through a w-ide rangeclose to the free surface become larger; tlierefore it becomes (figure 14(b)) including elongated structures; the reason formore difficult to resolve them- This is reflected hi the boundary these patterns is left for future studies- Since we are findinglayers being thinner (figure 13¿ than the ones shown in figure some inaccuracies in the WC2 simulations, it makes sense to6, corresponding to Re 500.observe the radial particle density graphs for a late frame of

6* international SPHERIC workshopFi1 M* *JvvrT*t rTHamburg, Germany. June, 08-10 2011iJii1—i—iItITILV11r11 ii HI i ¡i.11tf \i L T- - - - -i* «ft'-"1-A-*-11JlXA i.i1 iI (ai WC2: ftLiiiiiiiiI1 1llilAlftiMJÜii/xcu:iaii¿14iiIIi(a)« /ff) O5 0CA& G OT\t21aC1Í121-4Ll2LiS1I'eCft-*LC b RGKfig. 9.Particle density. Re 500. dx/hH/d* 400,1 0,5. Hf/H) 0 3 17.641.L21.41-01-S2{b) rí?/fl) M 17-6410.15r.Fig. ILParticle density graph for RGK, dx/h 2/9, U/dx 200.WCJX/IFQ.5RGKdxft 2#Analytic SolutionÍARc 500ü; R'dx 4Ü0L2SPH RGv.i'isio Kr-nrlSPH WenUntd K*m?lI(IS0.05 -0.60.4Analytic Solutionm\M,029Fig. 10.929.49.69.8 10Kinetic eueigy. RGK, dx/h 0 3 , 2 / 9 , Re 500, H / ¿ T 200Fig. 12,these simulations. This is done in ligure 15. where il can beappreciated that Ihe dispersion in the particles position for theWC2 is yet small, while large, as expected, for the RGK one.If we decrease the ratio ttxfh in order to get closer Io Ihecontinuum, the noise in Uie simulations with the WC2 kerneldisappear and the decay rate tends to the exact one (figure16). With ihe RGK the particles cluster in some regions ofthe flow domain (figure 17), in a similar pattern like the oneKinetic energy decay, dx/h 0 5 , Re 5000, ff/tf.v 400corresponding to the tensile instability phenomenon- If we takea look at the radial particle density tunction {figure 17) we seethat the clustering is reflected in the density reaching its peakfor q - 0* This clustering does nol take place under the sameconditions for the WC2 simulations.

6* international SPHERIC workshopHamburg, Germany. June, 08-10 2011 urroHUi H 81fa) WC200Í04ae x.nFig. 13. Vtorticity, Re 5000. WC2. dxlh 0.5. tig/ffp5 0.361,Hjdx 400« G.3QAOfifcS1111.4I .aIBib) RGKFig. 15. Parlkle densilv. Re 500ft///d* 400.M*w 17.641. ix/k 05.R 3 I0Q-1WroJhm4Kc«Ttr[J I Í I » fl« «M MC » ¡ 34" « K OM tt Re-5000Oft- HI,;.* .Fig. 16, Kinetic eneigy decay. WC2.to 5000.fl/rf* 400o-iMevolution of a standing wave lias been simulated using bothkernels. It is a significant test case because a analytical solutionlor the kinetic energy decay rate is available. Simulationshave been conducted varying the Reynolds number, resolutionand number of neighbors- It has been found that WendlandFig. 14. Vbiticily. Be 5000. ¿t/ft O.S.l&/!QM 7.920. H/dx 400. kemel simulations provide better resulls. For Ihe low Reynoldsnumber cases, the renonnalized Gaussian kernel simulationsshow a noisy vorticity field but no clustering takes place. Forthe highest Reynolds number case and for the largest neighborsV, CONCLUSIONS AND FUTURE WORKnumber, the renonnalized Gaussian kemel simulations presentThe aim of this paper has been lo discuss the performance a clustering of neighbors due to the onset of the tensileof two significant kernels in simulating free surface viscous instability* which in the same conditions, does not take placeRows. The kernels considered have been Ihe 5th degree class for the Wendland kemel simulations. We have not been2 Wendland kemel and a renonnalized Gaussian Kernel The able to link the discrepancies found with the properties of'

6th international SPHERIC workshopHamburg, Germany. June, 08-10 2011 * * Jt*ft**0 « i « íff Oft Ü»* 0 » W?and from the Spanish Ministry for Science and Innovalionunder grant T R A 2 0 H M 6 9 8 8 Caracterización Numérica yExperimental de las Cargas Fluido-Dinámicas en el transportede Gas Licuado"*" - SOCOREFERENCES fJft'g?i't'J '¿ii'jj!8Ml'tJg H(b) zoomRg. 17. \ rtícíly. to 5000. á»/ft 2/9. ¡(gfff}"5 W201. ff/Ai' 400.in 1VUiox 1O.t'U1-4lfiLBÍ«7 -i*3A01i«0 AIrLf\v/04*06* .\*1f\*K:*14/**ieisaFig. IS. Particle density. Rt 5000, tígflífi* m 16.201, rfx/A 2/9.tf/dx 400.both kernels- Nonetheless, we think that Wendland kernelprovides belter results and should therefore be used on topof RenormaKzed Gaussian for free surface Bows simulations.ACKNOWLEDGMENTST h e research that backs this paper has received fundingfrom the European Community's FP7 under grant agreementn225967 NextMuSE", from the Centre for Ships and OceanStructures (CeSOSj, NTNU, Trondheim, within the "ViolentWater-Vessel Interactions and Related Structural Load" project[1J A. Colagrossi. M Antuono. A. Son to- Iglesias, D. Le Touzé. andP. Izaguiíre Alza. Iheoretical analysis of SPH in simulating treeSurface viscous flows" in 5th ERCOFTAC SPHERIC workshop, 2010.[2] J. J. Monaghan and R. A. Gingold. "Shock simulation by the particlemethod SPH;* X Com?. Phys*. vol. 52. no, 2. pp. 374-589. 1983*[3] J. R Monis. R J. Fox. and Y. Zhu, "Modeling low Reynolds numberincompressible flows using SPH/* / Comp. Phys . vol. I36 1997.[4J E Macia. J, M. Sánchez, A. SoutoTglesias* and L. M. González,"WCSPH viscosity diffusion processes in vortexflows/1InternationalJournal for Numerical Methods in Fluids, tin press), 5) J. LightliüL Hfrra in fluids. Cambridge University Press. 200!.[6J D. Molteni and A* Colagmssi, 'A ample procedure lo Improve thepressure evaluation in hydrodynamic context using the SPH." ComputerPhysics Communications, vol. ISO, pp. S61-S72, 2009.[7] H. Wendland, "Piecewise polynomial, positive detimte and compacUysupported radial functions of minimal degree/* drfv, CompuL Math.,vol. 4, no. 4. pp. 389-396. 1995. S] A. ColagiossL M. Antuono. and D. L. Touzá, Theoretical considerationson the free-surface role in the Smoothed -particle-hydrodynamics model»"Physical Review Ef vol. 79, no. 5. 2009.[9J A. Colagrossi and M. Landrint. "Numerical simulation of interfacialflows by smoothed particle hydrodynamics;1 J, Compr WIYJ. vol. 19Itpp. 44¿-475. 2003.(10] M. Robinson aud J. Monaghan. "Forced 2D wall-bounded turbulenceusing SPH" in 3nl SPHERIC worishop, 2008. pp. 7S-S4.[Ill M.Robinson, 'Turbulence and Viscous Mixing using Smoothed ParticleHydrodynamics»" Ph.D. dissertation, Department of MathematicalScience. Monash University 2009.112] J. Hongbin and D* Xin. "On cnteaons for smoothed particlehvdrodvnamics kernels in stable field." Journal of ComputationalPhysics, VOL 202. no. 2. pp. 699 - 709. 2005.(131 J. Morris. "Analysis of SPH with applications." PhX). dissertation.Mathematics Dpi. Monash Univk1 Melbourne Australia. 1996.[141 A. ColagrossL "A mesliless lagnutgian method for free-surface andinterface llows with fragmentation." Ph.D. dissertation. Universitá diRoma La Sapienza, 2005.[151 J. W. Swegle. D. L. Hicks, and S. W. Attaway. "Smoothed ParticleHydrodynamics Stability Analysis." Journal of Computational Physics.vol. 116. pp. 123-134. 1995. *(16} G. L Wu. R. Eaiock-T yloi; and D. M. Greaves. "Hie effect or viscosityon the transient free-surface waves in a two dimensional lank" / Engng*Maths., vol. 40, pp. 77-90, 200L(17) A. Cenedese. Transition from Laminar to Turbulent Flow inside anOscillating Tank " Mecamica. vol. IS, pp. 217-224, 1983.[18) I. Derairdzic and S. Muzaferija, "Numerical method for coupled fluidflow, heat transfer and area* analysis using unstructured movingmeshes with ceils of arbitrary topology " Computer Methods in AppliedMechanics and Engineering* vol. 125. no. 1-4. pp. 235 - 255. 1995.(19] M* Antuono. A. ColagrossL S. Marrone. aud C. Lugni, "Propagation ofgravity waves through an SPH «rheme with numerical diffusive terms*"Computer Physics Comm., vol. 182. no. 4. pp. 866 - 877, 2011.(20} A. Colagrossi, G- Colicchio. C Lugnt and M, BroccWni "A study ofviolent sloshing wave impacts using an improved SPH method." Journalof Hydraulic Research, vol- 48. no. Extra Issue, pp. 94*104. 2010,[21] M. Hiero, P. Espafiol. and N. A. Adams. "Implicit atomistic viscositiesin smoothed particle hydrodynamics." Phys, Rev* Et vol. 82, no. 4.2010*[22} N. J. Quinlan. M. Lastiwka* and M. Basa, "Truncation error in meshfree particle methods." International Journal for Numerical Methods inEftgmeering.TO).66. no. 13. pp. 2064-2085. 2006.(23J A. Amicarellt J.-C. Marongiu. F. Leboeut J. Leduc. and J. Caro. "SPHtruncation enor in estimating a 3D function" Computen/ fc Fluids,vol. 44. no. I, pp. 279 - 296. 2011.37

Wendland [71 devised a family of radial interpolation functions that have compact support and that are definite positive. We have used the 5th degree class 2 Wendland kernel (WC2 from now on) reformulating it so that it has a 2// radius support, similarly to Robinson et al- [10], [11].