Mzir Gholm i Korzni*,Sergio A.Glindo-Torresb,Alexnder Scheuermnn Dvid J.W illims
aSchool of Civil Engineering,The University of Queensland,St Lucia,Brisbane,QLD 4072,Australia
bSchool ofMathematics and Physics,The University of Queensland,St Lucia,Brisbane,QLD 4072,Australia
Parametric study on smoothed particle hydrodynam ics for accurate determ ination of drag coefficient for a circular cylinder
Maziar Gholam i Korzania,*,Sergio A.Galindo-Torresa,b,Alexander Scheuermanna, David J.W illiamsa
aSchool of Civil Engineering,The University of Queensland,St Lucia,Brisbane,QLD 4072,Australia
bSchool ofMathematics and Physics,The University of Queensland,St Lucia,Brisbane,QLD 4072,Australia
Abstract
Simulationsof two-dimensional(2D)flow pasta circular cylinderw ith the smoothed particle hydrodynamics(SPH)methodwere conducted in order to accurately determ ine the drag coefficient.The fluid wasmodeled asa viscous liquid w ith weak compressibility.Boundary conditions, such as a no-slip solid wall,inflow and outflow,and periodic boundaries,were employed to resemble the physical problem.A sensitivity analysis,which has been rarely addressed in previous studies,was conducted on several SPH parameters.Hence,the effects of distinct parameters,such as the kernel choicesand the domain dimensions,were investigated w ith the goalof obtaining highly accurate results.A range of Reynolds numbers(1-500)was simulated,and the results were compared w ith existing experimental data.Itwas observed that the domain dimensions and the resolution of SPH particles,in comparison to the obstacle size,affected the obtained drag coefficient significantly.Other parameters,such as the background pressure,influenced the transient condition,but did not influence the steady state atwhich the drag coefficientwas determ ined.
?2017 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http:// creativecommons.org/licenses/by-nc-nd/4.0/).
Keywords:Smoothed particle hydrodynamics;Drag coefficient;Reynolds number;Sensitivity analysis;Viscous flow
Recent advancements in high-performance computing technologies have had significant impacts on mesh-free methods such as smoothed particle hydrodynam ics(SPH)in computational fluid dynamics(CFD).Although SPH is computationally expensive due to intensive interactions between integration points(in am isleadingmanner also denoted as particles),it may be advantageous in several situations, including moving boundaries(e.g.,free-surface flow),complex boundary geometries,and shock wave simulations. Consequently,SPH hasgrown in popularity since its inception in 1977 for astrophysics purposes(Gingold and Monaghan, 1977;Lucy,1977).Furthermore,recent developments in this method have extended its applicability to several engineering and scientific areas(Liu and Liu,2005;Monaghan and Gingold,1983)including fluid and solid dynam ics(Gray et al.,2001),coastal management(M irmohammadi and Ketabdari,2011),and geotechnical engineering(Bui et al., 2007).In order to ensure the best performance,the accuracy, stability,and validity of SPH still need to be investigated thoroughly,since it is rather a new numerical approach.
Despite the grow ing popularity of SPH,there have been few studies on the sensitivity of results to SPH parametersand the influence of SPH variables on the accuracy of results.On the other hand,this method is better known for simulatingfluidmotion in games(Gourlay,2014;Nie etal.,2015).As a result,the qualitative visualization of fluid flow has almost always been the main concern,rather than the accuracy in fluid simulations in engineering applications.Takeda et al. (1994)simulated flow past a circular cylinder for Reynolds numbers less than 55,and their results perfectly matched the results from the finite differencemethod(FDM)and experimental data.Although they developed a no-slip boundary and a new viscosity equation,they carried out sensitivity analyses only on the smoothing length and domain shape.Morris etal. (1997)used SPH to simulate the same problem for very low Reynolds numbers(Re≤1).They proposed an artificial velocity for boundary particles to fulfi ll the no-slip boundary condition.They also used a simplified viscosity equation based on a finite difference approximation for very low Reynoldsnumbers.Marrone etal.(2013)obtained good results for the same problem using theδ-SPH method(Antoci et al., 2007).They also compared results for circular and square shapesw ith those obtained from FDM,andmostly studied the effectsofobstacleshapeand vorticesgeometries.Nonetheless, a series of questions,including selection of a smoothing function,smoothing length,and different viscosity equations, as well as impacts of background pressure and the speed of sound,have remained unsolved.
The capability of thismethod should be examined quantitatively in relation to awell-documented problem such as flow over a bluff body(Anderson,2007),to provide an appropriate benchmark to verify the accuracy and validity of a numerical method as well as a new ly developed code.In this study,a weakly compressible SPH method was adopted(Monaghan, 2005).The drag coefficient,as themain output,was studied for different SPH variables,including viscosity equations, kernels,smoothing lengths,speed of sound,and background pressures,as well as the domain parameters,including dimensions and resolution numbers.This study provides a practical approach to increasing the accuracy of future SPH simulations.
This paper is structured as follows:the SPH method is described in the subsequent section,the third section briefl y explains the developed code,flow past a circular cylinder is discussed intensively in the fourth section,and the paper closes w ith a discussion of this work aswell as planned extensions of it.
The SPH method,originally developed for astrophysical purposes(Gingold and Monaghan,1977;Lucy,1977),is basically an interpolation technique.A comprehensive review of this method is presented in Liu and Liu(2005)and Monaghan(1994).In SPH,the computational domain is discretized into a finite number of particles(or integration points).These particles carry material properties such as velocity,density,and stress,and move w ith the material velocity according to the governing equations.The material propertiesof each particle are then calculated through the use of an interpolation process over its neighboring particles (integration domain)(Bui and Fukagawa,2013).The interpolation process is based on the integral representation of a field function.Numerous scientists have recently shown interest in thismethod and introduced details on the derivation and formulation of SPH(Li and Liu,2004;Liu and Liu, 2005,2010).Below,the SPH method is introduced in detail against the background of this study.As is generally known,fluid motion is governed by the continuity and momentum(Navier-Stokes)equations,which are formulated in SPH as follows:
The continuity equation:
Themomentum equation:
wheretis time;Wabis the kernel(smoothing function),which should have particular characteristics in order to approximate Dirac's delta function(Li and Liu,2004);?adenotes the gradient with respect to the coordinates of particle a;the subscriptsa and b denote the integration pointand particles in the neighborhood of particle a,respectively;and the particle field variables are v,P,m,μ,andρ,which represent the velocity,pressure,mass,dynam ic viscosity,and density, respectively.The last term in themomentum equation denotes the viscosity,which w ill be discussed in Section 2.2.
2.1.Equation of state
In a weakly compressible SPH(WCSPH)method,an equation of state(EOS)mustbe used to correlate density and pressure.As shown in Eq.(2),the pressure plays an essential role,so the estimation of pressure from the density field is of paramount importance in theWCSPH method.Three options are available:
whereCsis the speed of sound and should be greater than 10U,w ithUbeing the upstream velocity of the flow (Monaghan,1988);andρ0andP0are the density and background pressure at rest(initial condition),respectively.
EOSalso hasasignificanteffecton tensile instability,which isawell-documented issue in SPH,bymaintaining a persistent positive pressure.SPH particles repel each other when the pressure ispositive,sim ilarly to atoms,and attracteach other in the case of negative pressure,unlike atoms.However,the attraction causes SPH particles to form clumps,resulting in tensile instability(Monaghan,2000;Swegle etal.,1995).
The ideal gas EOS(Eq.(3))(Morris et al.,1997)avoids clustering of particles since it always maintains a positive pressure,but it requires very small time steps due to the possibility of localized high pressure.An alternative approach is to choose larger time steps in Eq.(4),which is computationally more efficient(Fang et al.,2006).However,Eq.(4) does suffer from particle clustering due to the development ofnegativepressure.For the flow pastan obstacle,thepressure becomes negative in thewake of the bluff body.This problem can be overcome through introduction of the background pressureP0in Eq.(4)(Marrone et al.,2011,2013;Morris, 1996).Consequently,Eq.(5)was chosen for this study in order to keep the pressure field positiveand to avoid the tensile instability.It should be noted that the background pressure is only applicable to problems w ith a confined domain.In the case of a free surface problem,there is no need to use this term,since thegravity acceleration can almostalwayskeep the pressure field positive.Alternatively,an artificial pressure (stress)should be used to deal w ith the tensile instability (Monaghan,2000).
2.2.Viscosity equation
In SPH,the chosen viscosity not only determines the behavior of the liquid,but also influences the performance of the numerical calculation.Therefore,an artificial viscosity is often employed in SPH instead of the real viscosity component of the Navier-Stokes equation(the last term of Eq.(2)) (Gholam i Korzani et al.,2016).The artificial viscosity was proposed by Monaghan and Gingold(1983)to dampen nonphysical numerical oscillations and to simulate shock waves (Monaghan,2012).A lthough this viscosity conserves total linear and angular momentum exactly,it is not a suitable approach tomodeling viscous fluids.Therefore,an appropriate viscosity term,which should be properly adapted to the SPH scheme,isneeded to simulate the shear flow more realistically. Several viscosity equations from the literature were used and are compared below:
(1)Morris et al.(1997):
whereraband vabare the distance and velocity vector(vab=va-vb)between particles a and b,respectively.?W/?radenotes the value of the derivative of the kernel function w ith respect to the location of particle a(ra)forrab. This formulation conserves linear momentum exactly,while angularmomentum is only approximately conserved.
(2)Shao and Lo(2003):
Eq.(7)is am ixture of a standard SPH fi rst derivative,w ith a finite difference approximation for the fi rst derivative.
(3)Takeda etal.(1994):
wherenis the dimension of the problem,and xabis the position vector,which is defined as xab=xa-xb.A lthough the solution of this equation employs the second derivative of the kernel,it is the exact derivation of viscosity from the Navier-Stokesequations.Themost importantassumption in Eq.(8)is that the fluid is compressible,and the bulk viscosity,ζ,can be considered negligible.
(4)Incompressible viscosity:
Thisequationwasused in this study as itconsiders the fluid incompressible by assum ing the divergence term(?$v)to be equal to zero.
2.3.Kernel
The smoothing function,also called the kernel,is of utmost importance since the performance of SPH critically depends on it.Kernels must satisfy numerous conditions, such as unity(the volume integral is equal to one),compact support(the kernel must not have an infinite reach),and positivity(Li and Liu,2004).Kernels depend on the smoothing length,h,and a non-dimensional distance between particles,q=r/h.In this study,some well-known kernels were tested:
(1)Quadratic(Johnson et al.,1996):
whereαdis2/(πh2)in two dimensionsand 5/(4πh3)in three dimensions.This kernel was used to simulate high-velocity impactproblems.However,the second derivative of the kernel is a constant value,which is not suitable for Eqs.(8)and(9). (2)Cubic spline(Monaghan,1988):
whereαdis 10/(7πh2)in two dimensions and 1/(πh3)in three dimensions.Thiskernelhasbeen,so far,themostw idely used kernel in SPH literature.
(3)Quintic(Wendland C2)(Wendland,1995)
whereαdis 7/(4πh2)in two dimensions and 7/(8πh3)in three dimensions.In general,the higher the order of the kernel is,the greater the accuracy of the SPH scheme w ill be. Therefore,it provides the best comprom ise between accuracy and computational time.
(4)Quintic spline(Morris et al.,1997;Schoenberg,1946):
whereαdis 7/(47πh2)in two dimensionsand 1/(120πh3)in three dimensions.This approach resembles a Gaussian kernel w ith a compactsupport.It is important to note that this kernel has a radius of influence of 3hand,therefore,more particles interact w ith one another,potentially reducing the computational efficiency.
Fig.1 shows the calculated functions for the presented kernels and their fi rstand second derivatives,aswell as the Laplacian,asa function ofqwhenh=1 andn=2.In addition,the impacts of the introduced kernels on the problem are further discussed in Section4.2.1.Severaltheoreticalandmathematical studieshave been conducted on kernels to investigate their influenceson SPH performance(Dehnen and Aly,2012;Violeau and Leroy,2014).As a result,the kernel standard deviation is the most important parameter in comparing distinct kernels. However,the aim of the present study was to practically investigate the effectsof the chosen kernelon the resultof the calculations.
2.4.Boundary conditions
2.4.1.Solid boundary
Immovable solid bodies are usually built w ith parallel layers of fluid particles fixed in place but not allowed to move,and thus acting as a solid boundary(Dalrymple and Knio,2001).In other words,Eq.(2)is calculated for each boundary particle,but the location and velocity of each particle are not updated at the end of each time step. Nevertheless,forces acting on the solid body can be obtained by summ ing Eq.(2)over all fixed boundary particlesmultiplied by theirmasses.Themost significant advantage of this approach is that the computational treatmentof the system is considerably simplified,since no special considerations are necessary for boundary particles,except for the constraintof nomovement.Several layers of particles are needed to avoid the kernel support domain truncation,since the truncation produces some errors in estimating the key variables,such as the density.
The essential characteristic of a solid boundary is the enforcementof the no-slip condition.Morrisetal.(1997)used an antisymmetric approximationmethod,which was sim ilar to thatof Takeda etal.(1994),to extrapolate the velocity of free particles through fixed boundary particles.In simple terms,the fixed boundary particleshave an artificial velocity thatmirrors their free particle counterparts to ensure that the overall fluid velocity on the surface is zero.
2.4.2.Inflow/outflow boundary
The inflow/outflow algorithm proposed by Federico et al. (2012)was used to generate steady continuous flow.As shown in Fig.2,two new sets of boundary particles(inflow and outflow particles)are defined at the domain boundaries in order to impose distinct upstream and downstream flow conditions.The inflow/outflow particles affect the free fluid particlesw ithin the domain.The inverse isnot true,as the inflow/ outflow particles are unaffected by the free ones.The inflow/ outflow region covered by these particlesmust be at least as w ide as the kernel domain.
At the inflow boundary,the desired velocity and density are assigned to the inflow particles.They are distributed on a regular grid and move according to their assigned velocity. When an inflow particle crosses the inflow threshold,it becomes a free fluid particle and w ill move according to the governing equations.At the outflow boundary,it is possible to impose specific outflow conditions similar to the inflow condition.Nonetheless,a fluid particle that crosses the outflow threshold becomes an outflow particle and all its field variables are maintained constant w ith time,except for the location,which evolves according to its velocity.Further details of the computational procedure are reported in Federico et al.(2012).
2.4.3.Periodic boundary
For this type of boundary condition,the cubical simulation box is replicated throughout space to form an infinite lattice, which is often used to elim inate surface effects from the computation.This boundary condition was applied perpendicular to the flow direction(Fig.2)and implemented by considering that particles located at a boundary are linked to particles at the opposite boundary.Therefore,a particle that leaves a boundary immediately re-enters at the opposite boundary with the same velocity.
Fig.1.Comparison of kernels and derivative values for h=1 and n=2.
Fig.2.Initial sketch of boundary conditions.
An open source code,called PersianSPH(Gholam iKorzani etal.,2015),hasbeendeveloped inC++tosolve theintroduced equations.All calculations and neighbor searching were conducted in parallel using the OpenMP library.Therefore,the whole domain was divided into several segments,and then the cell-linked list approach(Hockney and Eastwood,1988)was used to find neighboring particlesinorder to calculate the forces between them in parallel.There isnoparticular lim itation to the numberof particles,unless the numberof threads is restricted. In addition,themodified Verletscheme(Monaghan,1994)was used for time integration.The time step size,Δt,of the integrationmethod used has to be lim ited forstability reasonsbased on several criteria,which can be summarized as follows:
Courant-Friedrichs-Lewy(CFL)condition:
Condition on the viscous diffusion:
Condition on the force per unitmass:
whereaiis theacceleration of particlei,andνdenotes the fluid kinematic viscosity.The chosen time step sizemustbe smaller than them inimum of these conditions.
4.1.Modeling case
Themain problem studied was the calculation of the drag coefficient for a two-dimensional(2D)circular cylinder under the influence of various SPH parameters.As shown in Fig.3, the periodic boundary was used in the vertical direction,and the inflow/outflow boundarieswere used to generate a steady flow in the horizontal direction.A solid boundary wasused to simulate the circular cylinder.
Parameters studied can be divided into two main categories:SPH factors such as kernel choices,and domain factors such as the domain shape.A llparametersand coefficients w ill be introduced and discussed thoroughly in Section 4.2. However,all the parameters considered are defined concisely in Section 2 and Fig.3,except for the resolution number, which is the ratio of the cylinder diameter to the initial distance between SPH particles at time zero(refer to Section 4.2.7).
Fig.3.Schematic view of domain and boundaries.
4.2.Numerical results
As described earlier,the primary goal of this study was to achieve accurate results for the drag coefficient,which can be compared w ith experimental data.To achieve the best result, severalmodifications are required from a base setup of input parameters defining the numerical model.Inevitably,some initialassumptionsmade to simulate the problem mustbeused for the investigated parameters at early stages.A fterwards,an appropriate value or equation is chosen to replace the initial assumptions based on the results of each modification stage. By modifying these parameters step by step,the drag coefficient,as themain outputof this study,gradually converges to the experimental value.Finally,the drag coefficient is calculated for numerous Reynolds numbers based on themodified and calibrated models.
For the initial calculations in this study,the Reynolds numberwas fixed at60 for all stages,unless otherw ise noted. The target drag coefficientCd,extracted from the available experimental data(Anderson,2007),was 1.41 for this Reynolds number.A square w ith side lengths eight times the cylinder diameterwasmodeled as the whole domain,and the cylinderwas located in themiddle.In addition,the smoothing length,the speed of sound,and background pressurewere 1.1 times the initial distance of particles,ten times the imposed flow speed,and 0˙07C2sρ0,respectively.The resolution number was 30,and the Morris viscosity equation(Eq.(6))was employed.The impacts of these factors and selection criteria w ill be discussed in Sections 4.2.1 through 4.2.7.
4.2.1.Kernel selection
Various kernels have been introduced in Section 2.3.As shown in Fig.1,there are significant differences between the kernels and their derivatives,so it seems obvious that the results for different kernels must be very different as well. However,it should be noted that the optimum smoothing lengths for each kernel are not identical.Therefore,a fair comparison of kernels requires the selection of the individual optimum smoothing length.
Simulations for each kernel were conducted using the properties introduced in Section 4.2.Thesmoothing lengthwas not changed,so the number of neighboring particles was maintained constant,except for the case of the quintic spline kernel.Asshown in Fig.4,thequintic spline isthemostsuitable kernel,w ithCd=1˙79,because this kernel shows the best match w ith the experimental results and itsCdtime history curve is smoother(w ithout any fluctuations)than that of the other kernels in the steady state condition.A lthough thiskernel iscomputationally lessefficient,itisthemostreliablekernel for theproblem considered.Therefore,thequinticsplinekernelwas used for the restof the study.The accuracy of resultsusing this kernelhasbeen investigated by Gholam iKorzanietal.(2014).
4.2.2.Smoothing length
After choosing a proper kernel,finding the optimum smoothing length,h,is the next critical step.The number of neighboring particles grows significantly w ith increasingh, which decreases the computational efficiency.In SPH,the smoothing length is dependent on particle packing and the initial distance between particles.Hence,the problem was studied w ith varioussmoothing lengths.The smoothing length equalsαΔx,whereΔxis the initial distance between particles at the initial time step,andαvaried from 0.9 to 1.3.
As shown in Fig.5,the simulations forα=0.9 and 1.0 ended earlier since they were highly unstable.The values of drag coefficientCdforαvalues of 0.9,1.0,1.1,1.2,and 1.3 were 1.46,1.59,1.79,1.79,and 1.79,respectively.Cdwas exactly the same forαvalues of 1.1,1.2,and 1.3,demonstrating that a state had been reached in which the smoothing length did not change the results.Since the computation time increased with increasingh,the optimal value forαwas 1.1, which wasmaintained for the rest of the study.It is worth noting that although this optimum value results in aCdvalue that is quite different from the experimental value of 1.41,the difference can be reduced further w ith the introduction of othermodifications.
4.2.3.Viscosity equation
As already explained in Section 2.2,all equations for modeling the viscosity effect have advantages and disadvantages.Four simulationswere conducted to determine the impactsof theseequationson the resulting valuesofCd.Thedrag force is usually divided into two components:frictional drag and pressure drag.Frictional drag is derived from the friction between the fluid and the surface and involves viscous flow.Pressure drag is derived from eddyingmotions and is associated w ith the formation of awake.Theviscousdrag forcealso individually shows the effectof the viscosity equation.Hence, two coefficients were calculated:the total and viscous drag coeffi cients.They are compared for the different equations used in Table 1.
Fig.4.Drag coefficient curve vs.time for different kernels.
Fig.5.Drag coefficient curve vs.time for various smoothing lengths.
As illustrated in Fig.6 and Table 1,there is no significant difference between the resulting drag coefficients.However, some discrepancieswere observed in the transient state condition,show ing that the results of Eqs.(8)and(9)were slightly smoother.Eqs.(6)and(7)were proposed for low Reynolds numbers.Since the results for the equation of Takeda et al. (Eq.(8))and theequation of incompressibleviscosity(Eq.(9)) were closer to theexperimentalvalueofCd,Eq.(8)waschosen for later simulations.There was no difficulty in obtaining the second derivative of the kernel,which is needed for theapplication of the viscosity equation of Takeda et al.(1994), because the quintic spline kernel as a high-order kernel was used.
Table 1Drag coeffi cients for different viscosity equations.
4.2.4.Background pressure
As already discussed in Section 2.1,the value for the background pressure,P0,should be defined in such away that no negative pressure is generated during simulations.In order to explore this,P0can be defined asβCsρ0,whereβis a coefficientused to adjustP0.Consequently,P0w illbecome zero ifβapproaches zero(i.e.,no background pressure). Conversely,the EOSwillbe transformed into thatof the ideal gas ifβapproaches unity.Fig.7 illustrates the fluctuations ofCdw ith increasingβ,but in all cases it reached the same steady statevalue forCd.Thisstudy attempted tomaintainβas low as possible to m inim ize fluctuations and maintain a positive pressure.A value ofβ=0˙07 was chosen for the remainder of the study.
4.2.5.Speed of sound
The speed of sound in the fluid,Cs,as another important component of the EOS,influences the selection of the time step and resulting density fluctuations.With the increaseof the speed of sound,the time stepmustbe reduced,and the density variations become negligible.In contrast,density variations w illbe increased by reducing the speed of sound,whichmeans that the fluid displays a higher compressibility.Therefore,a balance should be established between the time step and weak compressibility by optim izing the speed of sound in the fluid.
The speed of sound in the fluid is proportional to the fluid velocityU,leading toCs=λU.This is an appropriate approach since the idea is not to use realistic values for the speed of sound,but to keep theMach numberU/Cssmalland to reduce compressibility effects.In thisstudy,λvaried from 5 to 20.Fig.8 shows that all curves reached a particular convergence value,but thosew ith higher values for the speed of sound had greater fluctuations in the transient state.As a result,the speed of sound was chosen to be ten times the fluid velocity as suggested by Monaghan(1988).
Fig.6.Drag coefficient curve vs.time for various viscosity equations.
Fig.7.Drag coefficient curve vs.time for various background pressures.
4.2.6.Problem domain geometrical shape
Two geometrical shapes were considered in this study to model the whole domain of the problem:a square(L=W, whereLis the length andWis the w idth),and a rectangle (L=2W)(Fig.3).Six different side lengths,proportional to the diameter of the cylinder,were investigated for the square domain shape.As demonstrated in Fig.9,the drag coefficient decreased gradually w ith increasing side length,and finally reached a plateau,signaling the point when the results were independent of the domain size.Table 2 shows the drag coefficient and relative difference between the calculated and experimental(1.41)drag coefficients,and also shows that the square domain shape resulted in an accurate drag coefficient. Since the number of particles grows significantly w ith increasing side length of the problem domain,it may be preferable to reach a comprom ise in accuracy by reducing the domain size.Although the obtainedCdvalue for the square w ith a side length greater than 20Dis closer to the experimental value,it is not computationally beneficial.Hence,to establish a balance between accuracy and computational time, theCdvalue can be in an acceptable range for a squarew ith a side length of 16D.
On the other hand,a rectangular domain shape is advantageous,since the numberof particles is reduced to half thatof the square domain shape.To study the impact of the rectangular domain on the accuracy of the results,two different flow velocities were considered.The results show that the drag coefficient values for Reynolds numbers of 60 and 200 were 1.48 and 1.26,respectively,for the square domain,and 1.58 and 1.32,respectively,for the rectangular domain.Therefore, the discrepancies between the drag coefficients calculated for the square and rectangular domain shapes were less than 7% and 4%forReynoldsnumbersof60 and 200,respectively.The difference between the rectangular domain shape results and experimental data was around 13%for both flow conditions.
Fig.8.Drag coefficient curve vs.time for various speeds of sound.
Fig.9.Drag coefficient curve vs.time for square domain shapew ith various side lengths.
Table 2Drag coeffi cient for square domain shapew ith various side lengths.
In addition,under the rectangular domain shape,a shifted cylinder condition was studied.In this condition,the cylinder was shifted about 5.5Dto the left.Consequently,the calculated drag coefficient(for a Reynolds number of 200)was 1.49,which was 13%greater than that calculated for the centered cylinder.Thus,the best condition for both domain shapes is the centered cylinder.
The results for the square domain shapeweremoreaccurate than those for the rectangular domain shape.However,the rectangular domain shape is preferred for reducing the computational time.Accordingly,the square domain shape was used for Reynolds numbers less than and equal to 200, since these simulationswere not computationally demanding. In contrast,the rectangular domain shape was selected forRe=500.As w ill be seen later,even w ith the rectangular domain shape,the results are accurate enough compared w ith the experimental data.
4.2.7.Resolution number(spatial resolution)
As stated previously,the resolution number,D/Δx,is defined as the ratio of the cylinder diameter to the initial distance between particles.To demonstrate how this factor affected the results,some simulationswere carried out through variation of this number.Two different flow velocities(corresponding to Reynolds numbers of 60 and 200)were investigated and the resolution number varied from 10 to 80.It should be noted that the case ofRe=60 is a flow condition preceding the generation of vortices,and the case ofRe=200 is a flow condition follow ing the generation of vortices in the wake of the cylinder.
Fig.10.Drag coefficient curve vs.time for various resolution numbers.
Fig.10(a)illustrates that lower resolution numbersproduced some noise in the results prior to the generation of vortices.In contrast,the drag coefficientrosegently to reach a plateau after thegenerationofvortices,asillustrated in Fig.10(b).Therefore, thedrag coefficientwasstrongly dependenton this factor,while this dependency was more significant for higher Reynolds numbers.A generalresultof thisstudy istheobservation thatthe resolutionnumberisdependenton theReynoldsnumber.Fig.10 shows that the drag coefficient is free and independentof the resolution numberwhen the resolutionnumber isgreater than a certain value(e.g.,20 forRe=60 and 60 forRe=200).That means thevalue of the resolution numbershow ing no influence on the drag coefficientmustbe determ inedmanually for every Reynolds number,and is henceforth denoted as the optimum resolution number.
4.2.8.Reynolds number
Finally,the model calibrated using the above-mentioned seven steps(Sections 4.2.1 to 4.2.7)was used to study the dependency of the drag coefficienton the Reynoldsnumber.In this section,results obtained from SPH were compared to available experimental data(Anderson,2007).
Fig.11 shows a comparison between the results of this study,the results of Marrone et al.(2013),and the available experimental data.A lthough there were some discrepancies between the acquired numerical results and the experimental data for Reynolds numbers exceeding 200,a satisfying agreementw ith experimental resultswas achieved.The deviation of the calculations from experimental results can be explained by the chosen rectangular domain shape and the smaller domain size(refer to Section 4.2.6).
Fig.12 shows the velocity field on the left-hand side and velocity vectorson the right-hand side,around the cylinder for different Reynolds numbers,demonstrating agreement with experimental observations and previous studies(Marrone et al.,2013).In addition,Table 3 shows the m inimum required resolution numbers for eachRevalue exam ined in this study,demonstrating that the drag coefficient is independent of the particle size.It is important to note that the maximum resolution number(D/Δx=80)w illwork for everyRevalue considered in this study,and the choice of the individual values was made based on computational efficiency issues only.
The paper presents a comprehensive parametric study of the flow past a circular cylinder in two dimensions.The impacts of several SPH variables on the calculation of drag coefficients in the 2D problem were studied for Reynolds numbers ranging from 1 to 500.Through comparison w ith existing experimental data and results from other numerical tools,the results of this study showed strong agreementw ith the existing experimental data.
The significance of the studied parameters can be summarized as follows:
(1)The quintic spline kernelwas themost reliable kernel for this study,even though it was computationally less efficient.
Fig.11.Drag coefficient curve vs.Reynolds numbers for a circular cylinder obstacle.
Fig.12.Snapshots from velocity field(left)and corresponding velocity vectors(right)for different Reynolds numbers.
(2)Different viscosity equations did not have any noticeable impacts on the results.Therefore,the simplest equation (Eq.(6))can be used to improve the computationalefficiency. It should be noted that Eqs.(8)and(9)are beneficial for reproducing vortices as in the present study no turbulence modelwas adopted.
(3)Thebackground pressure played a vital role in avoiding the clustering of particles due to tensile instability,but its drawback was some numericalnoise in the results.Therefore, itmust bemaintained as low as possible to keep the pressure positive.
(4)Variation of the speed of sound did notaffect the steady state results.It only produced fluctuations in the transient condition.
(5)The square shape of the problem domain was more accurate than the rectangular shape.However,the rectangular domain shape was preferred for higher Reynolds numbers to reduce the computational time w ithout significantly compromising the accuracy.
(6)The most suitable side length for the square domain shape was 26D,but a square w ith a side length of 16Dwas used to establish a balance between the accuracy and the computational time.
(7)As SPH is computationally expensive,it is important to keep themodel size as smallas possible.Therefore,obtaining efficient problem domain sizes or shapes,and using proper boundary conditions(e.g.,the periodic boundary condition) are absolutely important.
(8)The best location of the obstacle for achievementof the most accurate results was in the center of both considered domain shapes.
(9)The resultswere significantly sensitive to the resolution number for different flow velocities.Therefore,the resolution numberwas individually optim ized for each Reynoldsnumber in order to reach a threshold independentof particle size.This was done to reduce simulation time,but it is clear that the highestused resolution numberw ill cover the complete range of Reynolds numbers considered in this study.
In conclusion,the presented results and discussion provide guidance and indications for the selection of parameters used to optimize numerical studies using SPH.
Table 3Resolution number for each Reynolds number.
Anderson,J.D.,2007.Fundamentals of Aerodynamics,fourth ed.M cGraw-Hill Higher Education,Boston.
Antoci,C.,Gallati,M.,Sibilla,S.,2007.Numerical simulation of fluidstructure interaction by SPH.Comput.Struct.85(11),879-890.http:// dx.doi.org/10.1016/j.compstruc.2007.01.002.
Bui,H.H.,Sako,K.,Fukagawa,R.,2007.Numerical simulation of soil-water interaction using Smoothed Particle Hydrodynamics(SPH)method.J. Terramech.44(5),339-346.http://dx.doi.org/10.1016/j.jterra.2007.10.003.
Bui,H.H.,Fukagawa,R.,2013.An improved SPHmethod for saturated soils and its application to investigate themechanisms of embankment failure: Case of hydrostatic pore-water pressure.Int.J.Numer.Anal.Methods Geomech.37(1),31-50.http://dx.doi.org/10.1002/nag.1084.
Dalrymple,R.A.,Knio,O.,2001.SPH modelling of water waves.In:Proceedings of the Fourth Conference on Coastal Dynamics.ASCE,Lund Sweden.
Dehnen,W.,A ly,H.,2012.Improving convergence in Smoothed Particle Hydrodynam ics simulations w ithout pairing instability.Mon.Notices R. Astronom.Soc.425(2),1068-1082.http://dx.doi.org/10.1111/j.1365-2966.2012.21439.x.
Fang,J.,Owens,R.G.,Tacher,L.,Parriaux,A.,2006.A numericalstudy of the SPH method for simulating transient viscoelastic free surface flows.J. New t.Fluid Mech.139(1),68-84.http://dx.doi.org/10.1016/ j.jnnfm.2006.07.004.
Federico,I.,Marrone,S.,Colagrossi,A.,Aristodemo,F.,Antuono,M.,2012. Simulating 2D open-channel flows through an SPH model.Eur.J.Mech. B/Fluids 34,35-46.http://dx.doi.org/10.1016/j.euromechflu.2012.02.002.
Gholam i Korzani,M.,Galindo-Torres,S.A.,W illiams,D.,Scheuermann,A., 2014.Numerical simulation of tank discharge using Smoothed Particle Hydrodynam ics.Appl.Mech.Mater.553,168-173.http://dx.doi.org/ 10.4028/www.scientific.net/AMM.553.168.
Gholam i Korzani,M.,Galindo-Torres,S.A.,Scheuermann,A.,2015.PersianSPH:AMulti-physical Smoothed Particle Hydrodynamics Code.http:// korzani.wixsite.com/persiansph/persiansph[Retrieved 29 October 2016].
Gholami Korzani,M.,Galindo-Torres,S.A.,Scheuermann,A.,Williams,D., 2016.Smoothed Particle Hydrodynam ics into the fluid dynamics of classicalproblems.Appl.Mech.Mater.846,73-78.http://dx.doi.org/10.4028/ www.scientific.net/AMM.846.73.
Gingold,R.A.,Monaghan,J.J.,1977.Smoothed particle hydrodynamics: Theory and application to non-sphericalstars.Mon.Notices R.Astronom. Soc.181(3),375-389.http://dx.doi.org/10.1093/mnras/181.3.375.
Gourlay,M.J.,2014.Fluid Simulation for Video Games.https://software.intel. com/en-us/articles/fluid-simulation-for-video-games-part-1[Retrieved 29 October 2016].
Gray,J.P.,Monaghan,J.J.,Sw ift,R.P.,2001.SPH elastic dynam ics.Comput. Methods Appl.Mech.Eng.190(49),6641-6662.http://dx.doi.org/ 10.1016/S0045-7825(01)00254-7.
Hockney,R.W.,Eastwood,J.W.,1988.Computer Simulation Using Particles. CRC Press,Bristol.
Johnson,G.R.,Stryk,R.A.,Beissel,S.R.,1996.SPH for high velocity impact computations.Comput.Methods Appl.Mech.Eng.139(1),347-373. http://dx.doi.org/10.1016/S0045-7825(96)01089-4.
Li,S.,Liu,W.K.,2004.Meshfree Particle Methods.Springer-Verlag,Berlin.
Liu,G.R.,Liu,M.B.,2005.Smoothed Particle Hydrodynam ics:A Meshfree Particle Method.World Scientifi c,New Jersey.
Liu,M.B.,Liu,G.R.,2010.Smoothed Particle Hydrodynamics(SPH):An overview and recent developments.Archives Comput.Methods Eng. 17(1),25-76.http://dx.doi.org/10.1007/s11831-010-9040-7.
Lucy,L.B.,1977.A numerical approach to the testing of the fission hypothesis.Astronom.J.82,1013-1024.
Marrone,S.,Antuono,M.,Colagrossi,A.,Colicchio,G.,Le Touz■e,D., Graziani,G.,2011.δ-SPH model for simulating violent impact flows. Comput.Methods Appl.Mech.Eng.200(13),1526-1542.http:// dx.doi.org/10.1016/j.cma.2010.12.016.
Marrone,S.,Colagrossi,A.,Antuono,M.,Colicchio,G.,Graziani,G.,2013. An accurate SPH modeling of viscous flows around bodies at low and moderate Reynolds numbers.J.Comput.Phys.245,456-475.http:// dx.doi.org/10.1016/j.jcp.2013.03.011.
M irmohammadi,A.,Ketabdari,M.J.,2011.Numerical simulation of wave scouring beneath marine pipeline using Smoothed Particle Hydrodynam ics.Int.J.Sediment Res.26(3),331-342.http://dx.doi.org/10.1016/ S1001-6279(11)60097-8.
Monaghan,J.J.,Gingold,R.A.,1983.Shock simulation by the particlemethod SPH.J.Comput.Phys.52(2),374-389.http://dx.doi.org/10.1016/0021-9991(83)90036-0.
Monaghan,J.J.,1988.An introduction to SPH.Comput.Phys.Commun. 48(1),89-96.http://dx.doi.org/10.1016/0010-4655(88)90026-4.
Monaghan,J.J.,1994.Simulating free surface flows with SPH.J.Comput. Phys.110(2),399-406.http://dx.doi.org/10.1006/jcph.1994.1034.
Monaghan,J.J.,2000.SPH w ithout a tensile instability.J.Comput.Phys. 159(2),290-311.http://dx.doi.org/10.1006/jcph.2000.6439.
Monaghan,J.J.,2005.Smoothed particle hydrodynamics.Rep.Prog.Phys. 68(8),1703-1759.http://dx.doi.org/10.1088/0034-4885/68/8/R01.
Monaghan,J.J.,2012.Smoothed particle hydrodynam ics and its diverse applications.Annu.Rev.Fluid Mech.44,323-346.http://dx.doi.org/ 10.1146/annurev-fl uid-120710-101220.
Morris,J.,1996.Analysis of Smoothed Particle Hydrodynamics w ith Applications.Ph.D.Dissertation.Monash University,Melbourne.
Morris,J.,Fox,P.,Zhu,Y.,1997.Modeling low Reynolds number incompressible flows using SPH.J.Comput.Phys.136(1),214-226.http:// dx.doi.org/10.1006/jcph.1997.5776.
Nie,X.,Chen,L.,Xiang,T.,2015.Real-time incompressible fluid simulation on the GPU.Int.J.Comput.Games Technol.http://dx.doi.org/10.1155/ 2015/417417.
Schoenberg,I.J.,1946.Contributions to the problem of approximation of equidistant data by analytic functions,Part B:On the problem of osculatory interpolation,a second class of analytic approximation formulae.Q. Appl.Math.4(2),112-141.
Shao,S.,Lo,E.Y.M.,2003.Incompressible SPH method for simulating Newtonian and non-New tonian flows with a free surface.Adv.Water Resour. 26(7),787-800.http://dx.doi.org/10.1016/S0309-1708(03)00030-7.
Swegle,J.W.,Hicks,D.L.,Attaway,S.W.,1995.Smoothed particle hydrodynam ics stability analysis.J.Comput.Phys.116(1),123-134.http:// dx.doi.org/10.1006/jcph.1995.1010.
Takeda,H.,M iyama,S.M.,Sekiya,M.,1994.Numericalsimulation of viscous fl ow by smoothed particle hydrodynam ics.Prog.Theor.Phys.92(5), 939-960.
Violeau,D.,Leroy,A.,2014.On the maximum time step in weakly compressible SPH.J.Comput.Phys.256,388-415.http://dx.doi.org/ 10.1016/j.jcp.2013.09.001.
Wendland,H.,1995.Piecew ise polynom ial,positive definite and compactly supported radial functions ofminimal degree.Adv.Comput.Math.4(1), 389-396.http://dx.doi.org/10.1007/BF02123482.
Received 16 November 2016;accepted 21 March 2017 Available online 13 June 2017
This work was supported by the Australian Research Council Discovery Project(Grant No.DP120102188).
*Corresponding author.
E-mailaddress:m.gholam ikorzani@uq.edu.au(MaziarGholam iKorzani).
Peer review under responsibility of Hohai University.
http://dx.doi.org/10.1016/j.wse.2017.06.001
1674-2370/?2017 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http:// creativecommons.org/licenses/by-nc-nd/4.0/).
Water Science and Engineering2017年2期