Modeling Low Reynolds Number Incompressible Flows

ARTICLE NO.CP975776

Modeling Low Reynolds Number Incompressible Flows Using SPH

Joseph P.Morris,Patrick J.Fox,and Yi Zhu

School of Ci v il Engineering,Purdue Uni v ersity,West Lafayette,Indiana47907

E-mail:jpmorris@ecn.purdue.edu

Received December16,1996;revised June12,1997

those obtained using grid-based methods tailored for spe-

cific problems.Moreover,SPH can be more computation-The method of smoothed particle hydrodynamics(SPH)is ex-

tended to model incompressibleflows of low Reynolds number.ally expensive than alternative techniques for a given appli-For suchflows,modification of the standard SPH formalism is re-cation.

quired to minimize errors associated with the use of a quasi-incom-For typical astrophysical applications,SPH is used to pressible equation of state.Treatment of viscosity,state equation,

model compressiblefluids at high Reynolds number(usu-kernel interpolation,and boundary conditions are described.Simu-

ally ReՆ103[4]).Although SPH has been used to model lations using the method show close agreement with series solu-

tions for Couette and Poiseuilleflows.Furthermore,comparison free surfaceflows of inviscid incompressiblefluids[3],and withfinite element solutions forflow past a regular lattice of cylin-some research has been performed using SPH for com-ders shows close agreement for the velocity and pressurefields.pressible gases with Reynolds numbers down tofive[5], The SPH results exhibit small pressurefluctuations near curved

low Reynolds number(ReՅ1)incompressibleflows have boundaries.Further improvements to the boundary conditions may

not been modeled using SPH.The objective of this paper be possible which will reduce these errors.A similar method to that

used here may permit the simulation of otherflows at low Reynolds is to present modifications of the standard SPH formalism numbers using SPH.Further development will be needed for cases needed to simulate suchflows.Treatment of viscosity, involving free surfaces or substantially different equations of equation of state,kernel interpolation,and boundary con-state.ᮊ1997Academic Press

ditions are described.Simulations using the method show

close agreement with series solutions for Couette and

Poiseuilleflows and withfinite element solutions forflow

1.INTRODUCTION

past a regular lattice of cylinders.

The study of incompressiblefluidflows in which viscous

forces are either comparable with or dominate inertial

2.THE METHOD

forces has applications to many physical problems.Indus-

trial,biological,and environmental processes often involve 2.1.Background

flows of low Reynolds number(Re).Many applications

Using SPH thefluid is represented by particles,typically to thefields of environmental,mechanical,chemical,and

offixed mass,which follow thefluid motion,advect contact petroleum engineering likewise involve slow viscous in-

discontinuities,preserve Galilean invariance,and reduce compressibleflows throughfilters,substrates,porous mate-

computational diffusion of variousfluid properties includ-rials,and other potentially deformable structures.

ing momentum.The equations governing the evolution of Smoothed particle hydrodynamics(SPH),while origi-

thefluid become expressions for interparticle forces and nally developed for astrophysical applications[1,2],has

fluxes when written in SPH form.Using the standard ap-been applied successfully to a wide range of problems.

proach to SPH[6,7],the particles(which may also be SPH is a fully Lagrangian technique in which the numerical

regarded as interpolation points)move with the localfluid solution is achieved without a grid.An advantage of SPH

velocity.Each particle carries mass m,velocity v,and other is the relative ease with which new physics may be incorpo-

fluid quantities specific to a given problem.The equations rated into the formulation.It is also straightforward to

governing the evolution offluid quantities are expressed allow boundaries to move or deform and to model the

as summation interpolants using a kernel function W with interaction of severalfluid phases bounded by a free sur-

smoothing length h.For example,the density at particle face.SPH has certain advantages over otherfluid dynami-

a,␳a,may be evaluated using

cal methods,which may,for example,encounter difficulty

with deformable boundaries,multiphasefluids,free sur-

faces,and the extension to three dimensions.While SPH␳

aϭ͸b m b W ab,(1) is a versatile method,errors can sometimes be larger than

214

0021-9991/97$25.00

Copyright©1997by Academic Press

All rights of reproduction in any form reserved.

MODELING LOW RE INCOMPRESSIBLE FLOWS USING SPH

215

where W ab denotes

usually only important for faster flows involving shocks.The simulations presented here,therefore,use (6)to evolve density.

W ab ϭW (r ab ,h )

(2)

2.3.Dynamic Pressure

and

For low Reynolds number flows,local variations in the r ab ϭr a Ϫr b ,

(3)pressure gradient which force the fluid motion can be very small,compared with the hydrostatic pressure gradient.where r a denotes the position of particle a .The kernel This is of special significance to SPH,since pressure is typically takes the form

obtained using an explicit function of density and is only accurate to about 1%(see the following section).For many problems it is simpler to use SPH to model the dynamic W (r ab ,h )ϭ

1h ␴f ͉ͩr ab ͉h

ͪ

,(4)

pressure,p d ,defined as

where ␴is the number of dimensions and the function f

p d ϭp t Ϫp h ,

(7)

is discussed in Section 2.7.Other expressions for quantities at the particles are obtained by summation involving the where p t and p h are the total and hydrostatic pressures,

kernel or its derivatives.For example,the most common respectively.Since pressure appears in the Navier–Stokes SPH expression for the pressure gradient term is

equations only as a gradient,the effect of p h is that of a body force,

Ϫ

ͩ1␳

ٌp ͪa

ϭϪ͸b

m b

ͩ

p a ␳2a

ϩp b

␳2b ͪ

ٌa W ab ,(5)

Ϫ1␳ٌp t ϭϪ1␳ٌp d Ϫ1

␳

ٌp h ,(8)

where p a is the pressure at particle a and ٌa denotes the gradient with respect to the coordinates of particle a .For a ϭϪ1

␳

ٌp d ϩF ,

kernel of the form in (4),this pressure gradient formulation conserves momentum exactly,since forces acting between individual particles are antisymmetric.

where F is a force per unit mass.Using this approach,The following sections consider specific aspects of the pressure gradient driven flow through a periodic lattice can SPH formalism and the approaches necessary to simulate be easily simulated (see Section 3.3).Periodic boundary incompressible flows at low Reynolds number.

conditions are applied to all fluid quantities and the flow is driven by the effective body force.The dynamic pressure 2.2.Evolution of Density

is modeled using the equation of state (see (10)in Section 2.4).For simplicity,p is used in the following sections to If (1)is used when modeling incompressible free-surface denote the dynamic pressure p d .

flows,particle density is smoothed out at the edge of the fluid and spurious pressure gradients are induced at the 2.4.Equation of State

surface.To avoid this problem,Monaghan [3]initially set the density to a reference value and evolved particle densi-Most grid-based techniques model the flow of water ties according to the following SPH equation for continuity:

as incompressible since the speed of sound in water is usually large compared with bulk fluid motions (i.e.,a very low Mach number).Using SPH,fluid pressure is an explicit d ␳a

dt

ϭ͸b

m b v

ab иٌa W ab .

(6)

function of local fluid density and particle motions are driven by local density gradients.Therefore,it is necessary to use a quasi-incompressible equation of state to model Our test cases do not involve free surfaces,and thus (1)may be used to evolve particle densities.One disadvantage such flows with SPH.Since the actual equation of state for water would require a prohibitively small time step for of this approach is that ␳must be evaluated by summing over the particles before other quantities may be interpo-stability (by the CFL condition [8]),an artificial state equa-tion is used.The chosen sound speed is low enough to be lated.However,(6)allows ␳to be evolved concurrently with particle velocities and other field quantities,thus sig-practical,yet high enough to limit density fluctuations to within 3%.A similar approach has been used in grid-based nificantly reducing the computational effort.While (6)does not conserve mass exactly ((1)does,provided the techniques to model steady incompressible flow [9–11].Previous applications of SPH to incompressible flows

total number and mass of particles are constant),this is

216

MORRIS,FOX,AND ZHU

[3–12]have employed a state equation suggested by Bat-d v dt ϭϪ1

␳

ٌp ϩ␯ٌ2v ϩF ,(11)

chelor [13],

where ␯is the kinematic viscosity and F is a body force p ϭp 0

ͭͩ␳

␳0

ͪͲ

Ϫ1ͮ

,

(9)

per unit mass.Substituting a velocity scale V 0,a length scale L 0,and a body force per unit mass of magnitude F into (11),we find that the square of the sound speed should where Ͳϭ7and a zero subscript denotes reference quanti-be comparable with the largest of

ties.The choice of Ͳϭ7in (9)causes pressure to respond strongly to variations in density.Thus,perturbations to the density field remain small,even at high Reynolds numbers.c 2

ȁ

V 20ͳ,␯V 0L 0ͳ,FL 0

ͳ

,(12)

However,as the density fluctuations increase,small errors in density correspond to increasingly larger errors in pres-where

sure.For lower Reynolds numbers,more accurate pressure estimates are obtained using SPH if Ͳis taken to be unity (as it is in [9–11]),since errors in density and pressure ͳϭ

⌬␳␳0

.(13)

remain proportional.

In previous work involving incompressible fluids,the The first term in (12)corresponds to that derived by subtraction of 1in (9)was introduced to remove spurious Monaghan [3].The second and third terms ensure that boundary effects at a free surface.It is well established pressure forces are comparable with viscous and body that SPH is unstable when attractive forces act between forces,respectively.These conditions should only be re-particles [14–17,29].Consequently,for the simulations garded as a guide to an estimate of an appropriate sound presented in Section 3,this subtraction was found to lead to speed.After a simulation has been run initially at low numerical instabilities in regions of sustained low pressure.resolution and the actual variation in p is known,the value Since the test simulations (and many applications)have of c can be changed to produce the desired density vari-particles filling all space,this work uses

ation.

p ϭc 2␳,

(10)

2.5.Boundary Conditions

Initial applications of quasi-incompressible SPH in-where c is the speed of sound.

The sound speed must be chosen carefully to ensure volved high Reynolds number simulations of free surface flows interacting with free-slip boundaries.Such work em-both an efficient and accurate solution of a given problem.The value of c must be large enough that the behavior of ployed boundary particles which exerted strong repulsive forces to prevent SPH particles from penetrating solid sur-the corresponding quasi-incompressible fluid is sufficiently close to that of the real fluid,yet it should not be so large faces [3,19,20].To realistically model flows at lower Reyn-olds numbers,no-slip boundary conditions are needed.In as to make the time step prohibitively small.Using a scale analysis,Monaghan [3]argued that,for density to vary by addition,for the free surface flows considered by Mo-naghan [3],boundary particles do not contribute to the at most 1%,the Mach number of the flow should be 0.1or less.In fact,for typical smoothing lengths used with density of the free SPH particles,thus permitting the fluid to freely leave a solid boundary with no pressure-driven SPH,kernel interpolation is only accurate to within ap-proximately 1%.The principal cause of this variation is restoring force.For our work,boundary particles contrib-ute to the density of fluid particles such that pressure de-small fluctuations in density which inevitably occur as parti-cles move past one another.Thus,pressure gradients ob-creases when fluid and boundary particles diverge.It is possible to implement such a boundary condition using tained using a high sound speed are potentially noisy.Nev-ertheless,the velocities obtained are accurate if smoothed image particles [21].These images are created by reflecting fluid particles across the boundary with opposite velocities.either by XSPH [18]or viscosity.For many applications,a close estimate of bulk fluid motion is sufficient.However,This procedure works well for straight channels,but intro-duces density errors for curved surfaces.Takeda et al .[5]for other problems,accurate estimates of the pressure field are needed.We found that computed pressures are in close achieved a no-slip boundary condition using special bound-ary terms which mimic a half-space filled with SPH image agreement with other techniques when c is chosen such that the density varies by at most 3%.

particles.While these approaches have proved useful for compressible and moderate to high Reynolds number It is possible to estimate an appropriate value of c by considering the balance of forces in the Navier–Stokes flows,we found they did not give sufficiently stable results for our simulations.

momentum equation,

In this work,actual SPH particles are used to represent a no-slip boundary surface.These particles contribute to the usual SPH expressions for density and pressure gradients, and their own densities are also evolved.Evolving the densi-ties of boundary particles was found to better capture peak pressures than if boundary densities were kept constant. Obstacles within theflow are created by placing all the parti-cles on a regular lattice throughout the computational do-main and designating those particles falling within a solid object to be boundary particles.The advantage of this ap-proach is that a‘‘quiet’’start is guaranteed since the particle number density is initially constant throughout the compu-tational domain.However,it can lead to an imperfect

repre-

sentation of a curved boundary at low resolution.An alter-

native would be to place boundary particles on the surface FIG.1.Construction of artificial velocity for boundary particles to

simulate a no-slip boundary condition.

of each obstacle.One difficulty with this method is‘‘filling’’

thefluid space with a distribution of free particles corre-

sponding to constant density.In this case,the initial config-

uration can be‘‘relaxed’’to a quiet state before driving

ͱϭminͩͱmax,1ϩd B d aͪ.(15) forces are applied.The former approach was used in this

work for simplicity.

Afirst-order no-slip condition can be created if boundary

Numerical simulations have shown that good results are particles are moved with the velocity of the object to which

obtained ifͱmax is approximately1.5.If the boundary is in they are attached.However,this approach can produce sig-

motion,v a in(14)should be replaced by thefluid velocity nificant errors in velocity near boundaries.Antisymmetry

relative to the boundary.The artificial velocity v B is used to can be approximated by extrapolating the velocity of free

calculate viscous forces(see Section2.6)),but it is not used particles through the boundary to the position of the bound-

to evolve boundary particle positions.Consequently,the ac-ary particles.Ideally,a local estimate of the velocity gradi-

tual boundary velocity is used in(6)such that densities re-ents at the surface of the boundary would be used to assign

main consistent with(1).For concave surfaces,a similar these artificial velocities to interior points;however,such

method could be used in which the tangent plane is con-estimates would require a second summation over the parti-

structed by considering the nearest point on the curve to cles and,hence,a substantial increase in the computa-

each boundary particle(rather thanfluid particle).The two tional effort.

methods give identical results for a plane surface.

This work uses a simpler approach which is stable,accu-

There are alternative techniques which take boundary rate,and requires little extra computation.The method is

curvature into account.For example,Monaghan[19,20]in-similar to that used by Takeda et al.[5]to obtain a functional

troduced special boundary particle formulations in the con-form for the viscous force due to a solid boundary.Figure1

text of high Reynolds number free-surfaceflows which ex-illustrates the concept for a portion of a curved boundary.

plicitly account for local curvature effects when calculating Foreachfluid particle a,thenormal distance d a to thebound-

the contribution from boundary particles.This approach ary is calculated.This normal is used to define a tangent

can reduce the noise otherwise introduced by a discrete plane(a line in two dimensions)from which the normal dis-

boundary representation.Similar improvements to the tance d B to each boundary particle B is calculated.The ve-

boundary conditions employed here,incorporating ap-locity of particle a is extrapolated across the tangent plane,

proaches taken by other authors for high Reynolds number assuming zero velocity on the plane itself,thus giving each

problems,are the focus of further investigation. boundary particle velocity v BϭϪ(d B/d a)v a.In practice,the

discrete arrangement of boundary particles may permit a

2.6.Viscosity

fluid particle to closely approach the nominal curve describ-

ing the boundary.In such circumstances,the magnitude of Most implementations of SPH employ an artificial vis-v B must be restricted.Accordingly,the following formula is cosityfirst introduced to permit the modeling of strong used to calculate relative velocities betweenfluid and shocks[6,7].This viscosity is incorporated into the momen-boundary particles,tum equation

v abϭͱv a,(14)d v

a

dt ϭϪ͸b m bͩp a␳2aϩp b␳2bϩ⌸abٌͪa W ab,(16)

where

218MORRIS,FOX,AND ZHU

where

d v a

dt

ϭϪ͸b

m b ͩp a ␳

2a

ϩp

b ␳

2b

ٌͪa

W ab

(22)

ϩ͸b

m b

(Ȑa

ϩȐb

)v ab

␳a ␳b

ͩ1r ab

ѨW ab

a

ͪϩF a

,

⌸ab ϭ

Ά

ϪͰc Ȑ˜ab ϩͱȐ˜2

ab ␳ab ,if v ab иr ab Ͻ0;

0,

otherwise,

(17)

where F a is the body force evaluated at particle a .

Ȑ˜ab ϭh v ab иr ab

r 2ab

ϩ0.01h 2,

(18)

2.7.The Kernel

There are many possible choices of the function f in (4).and ␳ab is the average density of particles a and b .The Most SPH simulations employ a cubic spline kernel [27,28],

0.01h 2term is included to keep the denominator nonzero.Although this formulation has been used to model real viscosity [4,22],it produced inaccurate velocity profiles for our simulations.The advantage of (16)is that it guarantees f (s )ϭ

107ȏ

Ά

1Ϫ3s 2/2ϩ3s 3/4,if 0Յs Ͻ1;

(2Ϫs )3/4,

if 1Յs Ͻ2;0,

if s Ն2

(23)

conservation of angular momentum,which is important for applications involving relatively large fluid velocities or an unbounded fluid edge.Since this work involves low velocities and SPH particles fill all space,a more realistic (here normalized for two dimensions),since it resembles form of viscosity has been adopted.Previous expressions a Gaussian while having compact support.However,it has used to calculate realistic viscous forces have typically in-been shown that SPH can be unstable to transverse modes volved nested summations over the particles [23,24](and,when kernels with compact support are used [14,17,29].hence,twice the computational effort),or have directly As higher-order splines more closely approximating a employed second derivatives of the kernel [5].The disad-Gaussian are employed,these instabilities are reduced.vantage of using second derivatives is that interpolation is One reason for the poor performance of lower-order much more susceptible to error at low resolution (espe-splines is that the stability properties of SPH depend cially for low-order spline kernels [25]).Our method em-strongly upon the second derivative of the kernel.The ploys an SPH estimation of viscous diffusion which is simi-second derivative of the cubic spline is a piecewise-linear lar to an expression used in [26]to model heat conduction,

function,and,accordingly,the stability properties are infe-rior to those of smoother kernels.For example,when the quintic spline [27]is used,

ͭͩ

1

␳

ٌиȐٌͪv ͮa

ϭ

͸b m b (Ȑa ϩȐb )r ab

иٌa W ab

␳a ␳b

(r 2

ab

ϩ0.01h 2

)

v ab ,(19)

where Ȑϭ␳␯is the dynamic viscosity.This hybrid expres-sion combines a standard SPH first derivative with a finite f (s )ϭ

7

478ȏ

Ά

(3Ϫs )5Ϫ6(2Ϫs )5ϩ15(1Ϫs )5,

0Յs Ͻ1;(3Ϫs )5Ϫ6(2Ϫs )5,1Յs Ͻ2;(3Ϫs )5,

2Յs Ͻ3;0,

s Ն3,

difference approximation of a first derivative.By taking a Taylor expansion about particle a ,it can be shown that this expression is appropriate [26].This formulation conserves linear momentum exactly,while angular momentum is only (24)

approximately conserved.If the kernel takes the form of (4)then

transverse mode instabilities are negligible [14,29].For most applications using a cubic spline,these instabilities are not important since the resulting perturbations in density ٌa W ab ϭ

r ab ͉r ab ͉ѨW ab Ѩr a

(20)

typically peak at about 1%.However,for a quasi-incom-pressible equation of state,such variations are significant (see Section 2.3).Thus,for simulations involving very low and we can simplify (19)to

Reynolds numbers,it was found that a cubic spline rapidly produced significant noise in the pressure and velocity fields,whereas simulations employing the quintic spline ͭͩ

1

␳

ٌиȐٌͪv ͮa

ϭ

͸b m b (Ȑa ϩȐb )v ab ␳a ␳

b ͩ1r ab ѨW ab

Ѩr a

ͪ.

(21)

remained stable.Although the quintic spline is computa-tionally more intensive (by approximately a factor of 2),it was found to be the most reliable for our simulations.

Thus,the momentum equation is written as

MODELING LOW RE INCOMPRESSIBLE FLOWS USING SPH

219

Morris [29]presents a discussion of alternative kernels velocities are evolved according to (22).To ensure stability of the integration scheme,the time step is limited according which have superior stability properties to standard cubic splines with no extra computational expense.It is not clear,to the conditions set forth in Section 2.8.No-slip boundary conditions are simulated by assigning artificial velocities to however,if these kernels represent a significant improve-ment over the quintic spline for a given region of support the boundary particles using (14).

for this application.

3.MODEL TESTING AND VERIFICATION

2.8.Time Integration

3.1.Couette Flow

A modified Euler technique was used to perform the time integration.For stability,several time step criteria The first test case is Couette flow between infinite plates located at y ϭ0and y ϭL.The system is initially at rest.must be satisfied,including a CFL condition [8],

At time t ϭ0,the upper plate moves at constant velocity V 0parallel to the x -axis.The series solution for the time-⌬t Յ0.25

h c

,(25)

dependent behavior of this flow is

and additional constraints due to the magnitude of particle v x (y ,t )ϭV 0L y ϩ͸

ȍn ϭ1

2V 0n ȏ(Ϫ1)n

sin ͩn ȏL y ͪexp ͩϪ␯n 2ȏ2L 2t ͪ

,

accelerations f a [7],

(28)

⌬t Յ0.25min

a

ͩh

f a ͪ

1/2

,(26)

where v x is the fluid velocity in the x -direction.The flow was simulated using SPH for ␯ϭ10Ϫ6m 2s Ϫ1,L ϭ10Ϫ3m,and viscous diffusion,

␳ϭ103kgm Ϫ3,V 0ϭ1.25ϫ10Ϫ5ms Ϫ1,and with 50particles spanning the channel.This corresponds to a Reynolds number of 1.25ϫ10Ϫ2,using

⌬t Յ0.125h 2

␯

.

(27)

Equation (27)is based upon the usual condition for an ex-Re ϭ

V 0L

␯

.(29)

plicit finite difference method simulating diffusion.At suf-ficiently high resolution (small h )or large viscosity,(27)is Figure 2shows a comparison between velocity profiles typically the dominant time constraint.The choice of kernel obtained using (28)and SPH at several times including the and the arrangement of particles influences the coefficients steady state solution (t ϭȍ).The results are in close in (25)–(27).In particular,different splines can have differ-agreement (within 0.5%),confirming the accuracy of the ent ‘‘effective’’resolution lengths for the same value of h .approach used to evaluate viscous and boundary forces For example,use of a cubic spline (which is ‘‘narrower’’than with SPH.Lower resolution simulations completed with a quintic for the same smoothing length)may require 20particles spanning the channel were found to agree to slightly smaller coefficients in the above expressions.

within approximately 2%of the series solution values.

2.9.Summary of the Method

3.2.Poiseuille Flow

A summary of the method to simulate quasi-incompress-The second test case is Poiseuille flow between stationary ible flow at low Reynolds number using SPH can now be infinite plates at y ϭ0and y ϭL .The fluid is initially at presented.The particles are placed on a hexagonal lattice rest and is driven by an applied body force F parallel to filling all space and assigned a constant initial density (␳0).the x -axis for t Ն0.The series solution for the transient Boundary particles are defined as those initially lying inside behavior is

obstacles within the flow field and beyond solid walls.The equation of state for the fluid is given by (10)and the sound speed c is chosen such that density fluctuations are at most v x (y ,t )ϭF

2␯y (y ϪL )ϩ͸

ȍ

n ϭ04FL 2␯ȏ3(2n ϩ1)

3

(30)

3%.Inthe caseof periodicflow drivenby an appliedpressure gradient,the dynamic pressure (7)is used and the applied pressure differential is simulated by an applied body force sin

ͩ

ȏy L (2n ϩ1)ͪexp ͩϪ(2n ϩ1)2ȏ2␯

L 2

t ͪ

.(8).For stability at low Reynolds numbers,a quintic spline kernel is used with a smoothing length of 1.5times the initial nearest neighbor distance.Particle densities are evolved us-SPH was used to simulate Poiseuille flow for ␯ϭ10Ϫ6m 2s Ϫ1,

L ϭ10Ϫ3m,␳ϭ103kgm Ϫ3,F ϭ10Ϫ4ms Ϫ2,and 50particles

ing (6)to reduce the computational effort,and the particle

FIG.2.Comparison of SPH and series solutions for Couetteflow(Reϭ1.25ϫ10Ϫ2).

spanning the channel.This corresponds to a peakfluid spanning the channel agreed to within approximately2% velocity V0ϭ1.25ϫ10Ϫ5msϪ1and a Reynolds number of

of the series solution values.

1.25ϫ10Ϫ2using(29).A comparison between velocity

profiles obtained using(30)and SPH appears in Fig.3.

3.3.Flow through a Periodic Lattice of Cylinders The results are again in close agreement,with the largest

discrepancy being about0.7%for the steady state solution.The Couette and Poiseuille simulations tested the inter-Lower resolution simulations completed with20particles

action between viscous and body forces and the effective-

FIG.3.Comparison of SPH and series solutions for Poiseuilleflow(Reϭ1.25ϫ10Ϫ2).

FIG.5.Paths for comparison of SPH and FEM solutions.

FIG.4.Single cylinder within a periodic lattice.

were obtained by interpolating the particle quantities to a

50by50array of grid points using the quintic kernel.

Smoothing lengths of1and3grid spacings were used for ness of the no-slip boundary condition in the SPH model.the velocity and pressure,respectively.A greater amount However,theseflows are one dimensional and do not of smoothing was needed to remove smallfluctuations produce variations in dynamic pressure.A more challeng-from the pressurefield.Contours generated using this ing test of the method involvesflow through a square lattice method are inaccurate in the immediate vicinity of the cyl-of cylinders.This particular configuration has been studied inder.

extensively[30–32]as a simple model offlow through Figure6shows a comparison of velocity profiles obtained fibrous porous media.To solve this problem with SPH,a using SPH and FEM for paths1and2defined in Fig.5. single cylinder and its associated volume within the lattice The results obtained using SPH are in close agreement is considered(Fig.4).Flow is driven by a pressure gradient with those from the FEM throughout theflow domain. (modeled using an effective body force F),and periodic Corresponding contour plots of velocity magnitude are boundary conditions are applied to model an infinite peri-shown in Fig.7.Good agreement is obtained for the bulk odic array.of theflow,although the contour smoothing method inac-

curately represents SPH velocities near the cylinder.

3.3.1.Low Reynolds Number

Periodicflow past a cylinder was simulated using SPH

for Lϭ0.1m,␯ϭ10Ϫ6m2sϪ1,aϭ2ϫ10Ϫ2m,Fϭ1.5ϫ

10Ϫ7msϪ2,and cϭ5.77ϫ10Ϫ4msϪ1.Replacing L with a in

(29)and taking the velocity scale to be V0ϭ5ϫ10Ϫ5msϪ1

gives Reϭ1.The SPH simulation was run using approxi-

mately3000particles placed on a hexagonal lattice with a

nearest neighbor distance of0.002m.The cylinder was

modeled by considering all particles within its perimeter

to be of the type described in Section2.5.The particles

started from rest and steady state was reached after ap-

proximately1500steps.To investigate long-term behavior,

the simulation was continued for another6000steps such

that particle arrangements became disordered.The prob-

lem was also modeled using afinite element method(FEM)

program for steady incompressible viscousflow.Velocity

and pressure distributions from the two solutions were

compared by plotting values within one nearest neighbor

distance of the four paths described in Fig.5.The results

were also compared using contour plots.As the FEM em-

ploys a mesh to obtain a solution,it is relatively easy to FIG.6.Comparison of SPH and FEM velocity profiles along paths obtain contour plots.The corresponding plots for SPH

1and2for Reϭ1.

FIG.7.Contour plots of velocity magnitude using (a)FEM and (b)SPH for Re ϭ1(contour lines are labeled in units of 10Ϫ4m/s).

Figure 8shows the dynamic pressure along paths 3FEM solutions are in close agreement.The peaks in the pressure obtained using SPH on the boundary itself fell and 4(Fig.5).The arc of path 3was taken 0.002m (one SPH nearest neighbor distance)beyond the cylinder short of the FEM results by approximately 8%.The FEM better captures pressure extrema since grid-stretching boundary since the SPH boundary particles may not necessarily lie on the boundary surface itself.The SPH increases resolution in the vicinity of the cylinder.Corre-sponding contour plots of pressure given in Fig.9show dynamic pressure profile shows small local fluctuations in the vicinity of the cylinder.Elsewhere,the SPH and

that good agreement is again obtained for the bulk of the flow.This simulation was repeated with twice the particle resolution (approximately 11,000particles)and peak pressures were reproduced to within 5%.Pressure

values in the immediate vicinity of the boundary,how-ever,still exhibited small fluctuations.3.3.2.Very Low Reynolds Number

Flow through a periodic lattice of cylinders was also solved for L ϭ0.1m,␯ϭ10Ϫ4m 2s Ϫ1,a ϭ2ϫ10Ϫ2m,F ϭ5ϫ10Ϫ5ms Ϫ2,and c ϭ1ϫ10Ϫ2ms Ϫ1.Taking V 0ϭ1.5ϫ10Ϫ4ms Ϫ1,this gives Re ϭ0.03.Once again,the simulation was initialized with zero velocity and the steady state was reached after approximately 1500steps.However,the ini-tial lattice was relatively unchanged at this point.To dem-onstrate the long-term behavior of the method,the simula-tion was run for another 300,000steps.The final particle configuration,shown in Fig.10,was typical of those ob-served in these simulations once the initial lattice had bro-ken up.A comparison of velocities along paths 1and 2appears in Fig.11and velocity contour plots are shown in Fig.12.The results obtained using SPH are in close FIG.8.Comparison of SPH and FEM pressure profiles along paths 3and 4for Re ϭ1.

agreement with those obtained by the FEM.The flux

FIG.9.Contour plots of pressure using (a)FEM and (b)SPH for Re ϭ1(contour lines are labeled in units of 10Ϫ6Pa).

through the lattice is 1.70ϫ10Ϫ4m/s,which agrees to within were not fully captured by SPH,the discrepancies were somewhat smaller (about 5%).The simulation was re-3%with the analytical solution of Drummond and Tahir [31]for Stokes flow.A comparison of pressure fields pre-peated for fewer time steps with twice the resolution (ap-proximately 11,000particles)and peak pressures were ob-sented in Figs.13and 14also shows close agreement for the bulk of the flow,with similar fluctuations as observed tained with less than 4%error.Once again,however,small pressure fluctuations were observed near the boundary.

for Re ϭ1.Although pressure extrema on the boundary

FIG.10.Final particle positions corresponding to the results pre-FIG.11.Comparison of SPH and FEM velocity profiles along paths sented for Re ϭ0.03.Fluid and boundary particles are shown in black 1and 2for Re ϭ0.03.

and gray,respectively.

FIG.12.Contour plots of velocity magnitude using (a)FEM and (b)SPH for Re ϭ0.03(contour lines are labeled in units of 10Ϫ4m/s).

4.CONCLUSIONS AND FUTURE WORK

kernel result in a method which is stable and gives accurate estimates of velocity for problems of the type considered.Extensions of the method of smoothed particle hydrody-Peak pressures on solid boundaries are smoothed out at namics (SPH)have been presented which allow the simula-lower resolutions,but are achieved to within 5%for simula-tion of incompressible flows for low Reynolds numbers.tions involving over 10,000particles.However,pressures Results suggest that the proposed equation of state,viscos-in the vicinity of a circular boundary exhibited local fluctu-ity formulation,boundary conditions,and interpolation

ations of several percent.Further improvements to the boundary conditions may be possible which will reduce these errors (see Section 2.5).Additional tests involving convex and concave boundary surfaces are planned.

Straightforward extensions of the method are possible which will increase its utility.For example,it is possible to simulate the dispersion of a solute through the lattice of cylinders described in Section 3.3.A quantity which does not influence the flow (e.g.,the concentration of a dilute solute)may be evolved independently from one cell to the next while still assuming the velocity and density fields are periodic.As most of the computational expense is associated with locating and summing over nearest neighbors,the increase in work is moderate.A method similar to that used here may permit the simulation of other flows at low Reynolds numbers using SPH.However,further development will be needed for cases involving free surfaces or substantially different equations of state.The extension to three-dimensional problems is,in theory,straightforward.Obtaining steady-state results is relatively inexpensive,however,extended simulations in three-dimensions involving many crossing times at very low FIG.13.Comparison of SPH and FEM pressure profiles along paths 3and 4for Re ϭ0.03.

Reynolds numbers will probably necessitate the use of

FIG.14.Contour plots of pressure using (a)FEM and (b)SPH for Re ϭ0.03(contour lines are labeled in units of 10Ϫ3Pa).

zengleichungen der mathematischen physik,Math.Ann.100,32supercomputers.A modification of the method to include (1928).

surface tension is currently being undertaken to simulate 9.A.J.Chorin,A numerical method for solving incompressible flow the flow of multiphase fluids in porous media.

problems,J.Comput.Phys.2,12(1967).

10.E.Turkel,Preconditioned methods for solving the incompressible and

ACKNOWLEDGMENT/DISCLAIMER

low-speed compressible equations,J.Comput.Phys.72,277(1987).11.P.Tamamidis,G.Zhang,and D.N.Assanis,Comparison of pressure-This work was sponsored by the Air Force Office of Scientific Research,based and artificial compressibility methods for solving three-dimen-USAF,under Grant F49620–96-1–0020.The views and conclusions con-sional steady incompressible viscous flows,J.Comput.Phys.124,tained herein are those of the authors and should not be interpreted 1(1996).

as necessarily representing the official policies or endorsements,either 12.J.J.,Monaghan,Simulating gravity currents with SPH lock gates,

expressed or implied,of the Air Force Office of Scientific Research or Applied Mathematics Reports and Preprints 95/5,Monash Univer-the U.S.Government.

sity,1995.(unpublished)13.G.K.Batchelor,An Introduction to Fluid Dynamics.(Cambridge

REFERENCES

Univ.Press,Cambridge,UK,1967).

14.J.P.Morris,A study of the stability properties of SPH,Applied

1.L.B.Lucy,A numerical approach to the testing of the fission hypothe-Mathematics Reports and Preprints 94/22,Monash University,Mel-sis,Astron.J.83,1013(1977).bourne,Australia,1994.(unpublished)

2.R.A.Gingold and J.J.Monaghan,Smoothed particle hydrodynamics:15.D.S.Balsara,Von-Neumann stability analysis of Smoothed Particle

Theory and application to non-spherical stars,Mon.Not.R.Astron.Hydrodynamics:Suggestions for optimal algorithms,J.Comput.Phys.Soc.181,375(1977).121,357(1995).

3.J.J Monaghan,Simulating free surface flows with SPH,J.Comput.16.J.W.Swegle,D.L.Hicks and S.W.Attaway,Smoothed Particle

Phys.110,399,(1994).Hydrodynamics stability analysis,J.Comput.Phys.116,123(1995).4.P.Artymowicz and S.H.Lubow,Dynamics of binary-disk interaction.17.J.P.Morris,Stability Properties of SPH,Publ.Astron.Soc.Aust.13,

1.Resonances and disk gap sizes,Astrophys.J.421,651(1994).97(1996).

5.H.Takeda,S.M.Miyama and M.Sekiya,Numerical simulation of 18.J.J.Monaghan,On the problem of penetration in particle methods,

viscous flow by Smoothed Particle Hydrodynamics,Prog.Theor.J.Comput.Phys.82,1(19).

Phys.92(5),939(1994).19.J.J.Monaghan,Improved modelling of boundaries,Applied Mathe-6.W.Benz,Smooth particle hydrodynamics:A review,in NATO Work-matics Reports and Preprints 95/30,Monash University,Melbourne,shop,Les Arcs,France,19.Australia,1995.(unpublished)7.J.J.Monaghan,Smoothed Particle Hydrodynamics,Annu.Re v .20.J.J.Monaghan,Simulating gravity currents with SPH:III boundary

Astron.Astrophys.30,543(1992).

forces,Applied Mathematics Reports and Preprints 95/11,Monash University,Melbourne,Australia,1995.(unpublished)

8.R.Courant,K.Friedrichs,and H.Lewy,U

¨ber die partiellen differen-

226MORRIS,FOX,AND ZHU

21.L.D.Libersky,A.G.Petschek,T.C.Carney,J.R.Hipp,and F.A.Applied Mathematics Reports and Preprints95/18,Monash Univer-

Allahdadi,High strain Lagrangian hydrodynamics,J.Comput.Phys.

sity,Melbourne,Australia,1995.(unpublished) 109,67(1993).27.Schoenberg,I.J.,Contributions to the problem of approximation of

equidistant data by analytic functions,Q.Appl.Math.4,45(1946).

22.S.T.Maddison,J.R.Murray,and J.J.Monaghan,SPH simulations

of accretion disks and narrow rings,Publ.Astron.Soc.Aust.13,28.Monaghan,J.J.and Lattanzio,J.C.,A refined particle method for 66(1996).astrophysical problems,Astron.Astrophys.149,135(1985).

23.O.Flebbe,S.Mu¨nzel,H.Herold,H.Riffert,and H.Ruder,Smoothed29.J.P.Morris,Analysis of Smoothed Particle Hydrodynamics with Appli-

Particle Hydrodynamics:Physical viscosity and the simulation of ac-cations,Ph.D.thesis,Monash University,Melbourne,Australia,1996.

cretion disks,Astrophys.J.431,754(1994).30.N.Epstein and J.H.Masliyah,Creepingflow through clusters of 24.S.J.Watkins,A.S.Bhattal,N.Francis,J.A.Turner,and A.P.spheroids and elliptical cylinders,Chem.Enj.J.3,169(1972).

Whitworth,A new prescription for viscosity in Smoothed Particle31.J.E.Drummond and M.I.Tahir,Laminar viscousflow through Hydrodynamics,Astron.Astrophys.Suppl.Ser.119,177(1996).regular arrays of parallel solid cylinders,Int.J.Multiphase Flow10, 25.L.Brookshaw,A method of calculating radiative heat diffusion in515(1984).

particle simulations,Proc.Astron.Soc.Aust.6,207,(1985).32.G.W.Jackson and D.F.James,The permeability offibrous porous 26.J.J.Monaghan,Heat conduction with discontinuous conductivity,

media,Can.J.Chem.Eng.,3(1986).