A Flux Conserving Meshfree Method for Conservation Laws

1.6 MB·

by  Pratik Suchde

Thisisapreprint
ThefinalversionofthisarticlehasappearedinInternationalJournalforNumericalMethodsinEngineering
Thefinalfulltextisavailableonlineat:http://onlinelibrary.wiley.com/doi/10.1002/nme.5511/full
Pleasecitethisarticleasdoi:10.1002/nme.5511
A Flux Conserving Meshfree Method for Conservation Laws
PratikSuchde1,2∗,Jo¨rgKuhnert1,SimonSchro¨der1andAxelKlar2
1FraunhoferITWM,67663Kaiserslautern,Germany
2DepartmentofMathematics,UniversityofKaiserslautern,67663Kaiserslautern,Germany
7
1
0
2
SUMMARY
n
a Lackofconservationhasbeenthebiggestdrawbackinmeshfreegeneralizedfinitedifferencemethods(GFDMs).
J
In this paper, we present a novel modification of classical meshfree GFDMs to include local balances which

1 produce an approximate conservation of numerical fluxes. This numerical flux conservation is done within the
3
usualmovingleastsquaresframework.UnlikeFiniteVolumeMethods,itisbasedonlocallydefinedcontrolcells,

ratherthanagloballydefinedmesh.Wepresenttheapplicationofthismethodtoanadvectiondiffusionequation
]
A andtheincompressibleNavier–Stokesequations.Oursimulationsshowthattheintroductionoffluxconservation
significantlyreducestheerrorsinconservationinmeshfreeGFDMs.
N
. KEYWORDS: Meshfree methods; Conservation; Finite difference methods; Advection-diffusion equation;
h
Navier–Stokes;FinitePointsetMethod;FPM
t
a
m
[ 1 Introduction

1
v Generationandmanagementofmeshesisoftenthemostdifficultandtimeconsumingpartofnumerical
3
simulation procedures. This is further compounded for complex, time-dependent geometries. The
7
efficiencyofmeshgenerationdeterminestheoverallaccuracyandrobustnessofthesimulationprocess.
9
8 Toavoidthetaskofmeshing,severalclassesofmeshlessormeshfreemethodshavebeendeveloped.
0 Meshfreemethodsusethenumericalbasisofasetofnodes,whichneednotberegularlydistributed,
. tocoverthecomputationaldomain.Thesenodescouldeitherbemass-carryingparticlesornumerical
1
0 points.Foreachnode,theonlyconnectivityinformationrequiredisalocalsetofneighbouringnodes
7 overwhichapproximationsarecarriedout.
1 Smoothed Particle Hydrodynamics (SPH) [21] is one of the first and most widely used meshless
:
v methods.InSPH,thecomputationaldomainisdiscretizedbyparticleswhichcarrymass.Thus,mass
i conservationisdirectlyguaranteed.OneofthebiggestdrawbacksofSPHisthedifficultyinenforcing
X
boundary conditions [18, 24, 29]. While a lot of work has been done to address this issue, the SPH
r
a formulation does not naturally include treatment of most boundary conditions and extra effort is
needed to enforce them. One such method is the addition of ghost or dummy particles outside the
computational domain [18]. Another major drawback of SPH in its classical formulation is the so-
called particle inconsistency issue which results in the absence of a valid approximation order. The
original SPH does not even have C0 particle consistency for irregularly distributed particles and
boundaryparticles.Attemptstosolvethisproblemincludethekernelrenormalization[3]whichcomes
at the price of the loss of certain conservation properties for momentum and energy. In addressing
these issues, the Finite Pointset Method (FPM) evolved out of SPH [16]. FPM [13, 17, 34, 35] is a
meshfreemethodbasedonmovingleastsquares(MLS)approximationssuchthatboundaryconditions
can be directly enforced, without any extra effort. While SPH uses a set of mass-carrying particles
as a numerical basis, FPM and other meshfree GFDMs (for example, [14, 19, 28, 33]) use a set of
∗Correspondenceto:P.Suchde.E-mail:[email protected]
2 P.SUCHDEETAL.
numerical points, referred to as a point cloud. These numerical points do not carry mass, and thus,
pointscaneasilybeaddedorremovedduringasimulation[13],whichisnotpossibleintheparticle-
based SPH. Since point clouds can easily be modified locally, FPM is more adaptive than SPH, in
terms of the discretization of the domain. However, these advantages come with the drawback of
a lack of formal conservation. Despite the lack of conservation, the FPM has been shown to be a
robustmethodwithmanypracticalapplications[7,13,34,37].TheFPMisalsousedasthenumerical
basis of two commercially used meshfree simulation tools: NOGRID [25] and the meshfree module
ofVPS-PAMCRASH[36].Itmustbenotedthatinthemeshfreecontext,theacronymFPMisoftena
confusingone.ItisusedtorepresentnotonlytheaforementionedFinitePointsetMethod,butalsothe
Finite Particle Method [23] and the well-known Finite Point Method [26]. Thus, we henceforth drop
theacronymFPMandrefertotheFinitePointsetMethodundertheumbrellatermofGFDM.
Despite the name, most meshfree methods are not completely devoid of meshes. While meshfree
methods lack the use of a predefined mesh for domain discretization, many meshfree methods use
an easily generable mesh such that the solution is not too heavily dependent on the quality of that
mesh (see [20] for details). Several methods, which solve partial differential equations in their weak
formulation, use a so-called background mesh for numerical integration of the weak-forms. The
element-free Galerkin method (EFG) [2] uses a global background mesh, while the meshless local
Petrov-Galerkinmethod(MLPG)[1]requiresonlylocalbackgroundcells.Meshfreeparticlemethods,
suchasSPH,useapredefinitionofmassparticles,whichoftenrequiressomekindofamesh.Oneof
themethodstosolvetheproblemofparticledistortioninSPHisbyregularization,whichrequiresare-
meshing[39].Severalmethodsinthefamilyofpointinterpolationmethods(PIM)[22]usetessellations
orbackgroundcellsinso-calledT-schemesfortheselectionofsupportdomains.Point-basedmeshfree
GFDMsareoftenreferredtoas“trulymeshfree”.However,theytoousesomenotionofamesh.They
oftenuselocaltriangulationsforaccuratepost-processingoperationsorfortheidentificationofpoints
onthefreesurfaceoffluidflowsimulations.Ineachofthesecases,themeshesusedarealoteasierto
generateanddonotimposetherestrictionsthatmeshesinclassicalmeshedmethodsdo.
A few attempts have been made to overcome the lack of conservation in meshfree GFDMs [4, 5],
but they lose the completely local nature of the differential operators. They introduce a symmetry of
differentialoperatorcoefficientswhichleadstoadiscretedivergencetheorem.However,thatresultsin
aneedtosolvefordifferentialoperatorsinaglobalsystem.Thatincreasessimulationtimesignificantly
and makes those methods extremely inefficient on finely discretized point clouds and moving point
cloudswhichrequireare-computationofnumericaldifferentialoperatorsateverytimestep.
FiniteVolumeMethods(FVMs)havethefeatureoflocalconservationofnumericalfluxesfromone
controlvolumetoitsneighbouringone[8].Inthispaper,weaimtouselocalcontrolcells,inamanner
similar to FVMs, in a meshfree framework to introduce an approximate conservation of numerical
fluxes. A local control cell at a point is formed by the Delaunay tessellation of points in its support
domain. The differential operators are computed using an MLS approach which guarantees usual
monomialconsistencycoupledwithaconservationofnumericalfluxesforspecificfields.Formodeling
fluidflow,weuseafullyLagrangianframework.Sincepointcloudscanbemodifiedlocally,andthe
controlcellsarealsodefinedlocally,thismethoddoesnotfacetheproblemofcelldistortionpresent
in Lagrangian moving mesh methods. It must be noted that the idea of generalizing control volumes
used by FVMs in the context of meshfree methods is not a novel one. The Finite-Volume Particle
Method(FVPM)[10,11]generalizescontrolvolumestoincludethosewhichneednotbedisjoint.In
FVPM, grid generation is replaced by an expensive integration of partition of unity functions. The
Meshless-Finite Volume Method (MFVM) [12] generalizes Voronoi-based moving mesh methods.
MFVMusesavolumepartitioningthatamountstoaVoronoitessellationwithedgessmoothed.Unlike
boththeFVPMandMFVM,themethodpresentedhereisastrong-formmethod.
Thepaperisorganizedasfollows.InSection2,wegiveabriefoverviewofGFDMs.InSection3,we
introduceourproposedmethodbypresentingthedetailsofcontrolcellsinameshfreecontextandthe
construction of flux conserving differential operators. In Section 4 and 5, we present the application
of this method to an advection-diffusion equation and the incompressible Navier–Stokes equations
respectively.Finally,weconcludethepaperinSection6withadiscussionoftheproposedmethod.
2
FLUXCONSERVINGMESHFREEMETHOD 3
2 Meshfree Generalized Finite Difference Methods
FortheFinitePointsetMethodandothersimilarmeshfreeGFDMs,thecomputationaldomainΩ,with
boundary∂Ω,isdiscretizedusingacloudofN numericalpointswithpositions(cid:126)x ,i=1,...,N,which
i
includespointsbothintheinteriorandontheboundaryofthedomain.Thepointsareusuallyirregularly
spaced.Eachnumericalpointcarriesthenecessarynumericaldataoftheproblem.Eachpointihasaset
ofneighbouringpointsS whichcontainsn points,includingitself.TheneighbourhoodorsupportS
i i i
isdeterminedasS ={(cid:126)x :(cid:107)(cid:126)x −(cid:126)x (cid:107)≤βh},whereh=h((cid:126)x,t)istheradiusofthesupport,referred
i j j i
toassmoothinglengthorinteractionradius;andβ ≤1isapositiveconstant.Anumericalpointisnot
amasscarryingparticleandthus,pointscaneasilybeaddedordeletedtoadaptthepointcloudlocally.
The spatial distribution of points is described by three parameters: h, r and r . It is ensured
min max
thatnotwopointsarecloserthanr handthereexistsatleastonepointineverypossiblesphereof
min
radiusr hinthecomputationaldomain.Thus,thesmoothinglengthhalsodeterminesthespatial
max
discretization size. r and r usually have values of approximately 0.2 and 0.45 respectively.
min max
Furtherdetailsaboutpointcloudorganizationandmanagement,includingthesetupoftheinitialpoint
cloud,canbefoundin[7,13].
As the name suggests, numerical derivatives are approximated with a generalized finite difference
approach. For a function f defined at each numerical point i=1,2,...,N, its derivatives are
approximatedas
∂∗f((cid:126)x)≈∂˜∗f = (cid:88) c∗f , (1)
i i ij j
j∈Si
where∗ =x,y,xx,∆,etc.representsthedifferentialoperatorbeingapproximated,∂∗ representsthe
i
continuous ∗-derivative at point i, and ∂˜∗ represents the discrete derivative. For a point i, the stencil
i
coefficients c are found using a weighted least squares approach. The weighted sum of the stencil
ij
coefficientsisminimizedsuchthatthederivativesofmonomialsm∈Muptoacertainorder,usually
2,areexactlyreproduced.
(cid:88)
c∗m =∂∗m ∀m∈M, (2)
ij j i
j∈Si
(cid:88)
minJ = W (c∗)2, (3)
ij ij
j∈Si
whereW isaweightingfunction.Throughoutthispaper,weuseaGaussianweightingfunction
(cid:107)(cid:126)x −(cid:126)x (cid:107)2
W =exp(−4 j i ). (4)
ij h2+h2
i j
Foreachpointi,Eq.(2)andEq.(3)aresolvedtoobtainthenumericaldifferentialoperator,whichcan
thenbesubstitutedinthepartialdifferentialequationsbeingsolvedtoobtainthespatialdiscretization.
Henceforth, the differential operators found by Eq.(2) and Eq.(3) will be referred to as classical
GFDM.
3 Conservation
Consideraconservationlaw
∂φ
+∇·J=0, (5)
∂t
with J=J(φ). For sufficiently smooth φ and J, integrating Eq.(5) over the entire domain, and an
applicationofthedivergencetheoremleadstotheintegralformoftheconservationlaw,
d (cid:90) (cid:90)
φdV =− (cid:126)n·JdA. (6)
dt
Ω ∂Ω
Physically,Eq.(6)canbestatedastherateofchangeofenergywithinthedomainshouldbethesame
astheenergyfluxacrossitsboundary,whennoenergysourceorsinkispresent.Finitevolumemethods
3
4 P.SUCHDEETAL.
Figure1.Controlcellwithinalocalsupport.Centerpointismarkedwithacross.Nearestneighboursl∈I are
i
markedwithcircles.Interiorpoint(left)andboundarypoint(right).
[8] use a local balance on each discretization cell or control volume. FVMs solve Eq.(5) by directly
enforcingalocal,discreteversionofEq.(6).
Meshfree GFDMs, on the other hand, usually directly solve the strong form Eq.(5). The lack of a
discretedivergencetheoremcausestheabsenceofadiscreteformofEq.(6).Thiscombinedwiththe
localnatureofthedifferentialoperators,Eq.(1),resultsinconservationnotbeingensuredatadiscrete
level.WenowpresentanovelmodificationofGFDMdifferentialoperatorstoenforcealocaldiscrete
divergencetheorem.Thisresultsinanapproximatediscreteformofglobalconservation,Eq.(6).
3.1 AMeshfreeControlCell
For each point i, we consider the Delaunay tessellation of the n points in its support domain S .
i i
The Voronoi diagram forms the dual graph to the Delaunay tessellation. Among the tessellations,
the Voronoi cell containing point i is the only one of interest, and is used as the control cell over
whichthefluxbalanceiscarriedout.Figure1showssuchacellwithinthesupportdomain.Thepoint
i is associated with a volume V , taken to be the volume of this cell. V is also used for accurate
i i
post-processing.Thelocaltessellationscanalsobeusedinthecomputationofgeometricparameters
of the point cloud, such as the identification of free surface points in flows with open boundaries.
Alternatively,insteadoftheVoronoidual,acentroidaldualcouldalsobeusedtothesameeffect.
Using the local tessellation for the support S , we obtain a set of nearest neighbours I , in the
i i
Delaunay sense. For l∈Ii, we let i(cid:101)l denote the dual edge (in 2D) or face (in 3D) of the Delaunay
edge il. Further, we let Ail denote the area of i(cid:101)l and(cid:126)nil denote the unit normal of i(cid:101)l, pointing away
from i. For boundary points, the geometric area A and outward pointing unit normal(cid:126)n are used to
i i
truncateandclosethesemi-infiniteVoronoicell,asshowninFigure1.Tosimplifynotation,wedefine
theclosureI¯ ofI toincludethepointiitself,ifiisaboundarypoint.Further,welet(cid:126)n =(cid:126)n and
i i ii i
A =A .
ii i
We note that since the Voronoi cells are defined locally on the support domain of each point, and
not globally, they need not stitch together to form a global tiling of the computational domain. The
uniqueness property of Delaunay tessellations suggests that the Voronoi cells should stitch together
perfectly.However,thatonlyholdsforsufficientlylargesupportdomains.Forsmallsupportdomains,
itispossiblethataninsufficientnumberofpointsineachneighbourhoodleadstounsymmetricVoronoi
cells.Undersuchsituations,symmetryoftheinterfacesvalues,A andA ;(cid:126)n and(cid:126)n ,isviolated.
ij ji ij ji
Figure2showssuchanexampleofadjacentnon-symmetriccells.Wemeasuretheextentofstitching
togethertoformaconsistentglobalmeshbytheerror
(cid:80)N (cid:80) (cid:107)(cid:126)n A +(cid:126)n A (cid:107)
(cid:15) = i=1 l∈Ii il il li li . (7)
mesh 2(cid:80)N V
i=1 i
If i∈/ I , we set (cid:126)n =(cid:126)0 and A =0. If the cells are based on a global tessellation, (cid:126)n A =
l li li il il
−(cid:126)n A ∀i,l, which leads to (cid:15) =0. Table 1 shows the extent of stitching together of locally
li li mesh
4
FLUXCONSERVINGMESHFREEMETHOD 5
defined Voronoi cells for different support sizes β, and the corresponding average number of points
inthesupportdomain.Table1illustratesthatforlargeenoughsupportdomains,whichcorrespondto
high values of β, the intersection of local Voronoi cells is minimum. Throughout this paper, we use
supportsizesintherangeofβ ∈[0.75,1].Thisresultsinapproximately45−55pointsinthesupport
domainofeachpointinthreedimensions,andabout20pointsintwodimensions.
Figure2.A2Dexampleofnon-symmetriclocallydefinedVoronoicells(left),andthesamepointconfiguration
withsymmetriclocallydefinedVoronoicells(right).Thenon-symmetriccasesoccurwhenthesupportdomains
are not sufficiently large. The size of the support domains considered are β =0.6, mean(n )=11 (left) and
i
β =0.7,mean(n )=18(right).
i
Table1.ErrorsinformationofaglobaltilingbystitchingtogetherlocallyobtainedVoronoicells(in3D).The
domainconsideredisaunitspherewithh=0.5,N =1400.
β (cid:15) mean(n )
mesh i
0.5 3.31 16
0.55 1.08 21
0.6 0.33 27
0.65 9×10−3 35
0.7 6×10−15 41
Several meshfree methods, such as EFG [2] and PFEM [27], use a globally defined background
meshforavarietyofpurposes.Suchaglobalbackgroundmeshcouldbeusedforthecontrolcells,but
definingcellsbasedonlocaltessellationsonasmallnumberofpointsissignificantlyfaster,especially
when done in parallel. The concept of shapeless meshfree volumes for each point have also been
proposed by several authors [4, 5, 14, 15], but their use presents several problems in the present
context. Most importantly, the proposed meshless volumes are not closed, which, as we shall show
later,isnecessaryfortheconstructionofourfluxconservingdifferentialoperators.
3.2 FluxConservingDifferentialOperators
Forthedefinitionofournewnumericalfirstderivativeoperators,considerthedivergencetheorem
(cid:90) (cid:90)
∇·(cid:126)gdV = (cid:126)n·(cid:126)gdA, (8)
Ω ∂Ω
forsufficientlysmooth(cid:126)gand∂Ω.Settinggtobe(f,0,0),(0,f,0)and(0,0,f)alternatingly,ascalar
versionofEq.(8)canbeobtainedas
(cid:90) (cid:90)
∂kfdV = fnkdA. (9)
Ω ∂Ω
5
6 P.SUCHDEETAL.
Throughout this paper, the superscript k =x,y,z denotes the spatial dimension. A discrete local
versionofEq.(9)onthecontrolcellscanbewrittenasfollows
V ∂˜kf =(cid:88)F (f), (10)
i i il
l∈I¯i
where F (f) is a numerical flux function which depends on (cid:126)n , A , f and f , and possibly their
il il il i l
derivatives.OnepossibilityforF,toachievealocaldiscretedivergencetheorem,canbe
(cid:40)
f nkA ifi∈∂Ωandl=i,
F (f)= i i i (11)
il f nk A elsewhere,
il il il
with f = fi+fl. (cid:126)n and A are determined locally based on the control cell as mentioned in
il 2 il il
Section3.1. Throughout this paper, we consider only simple symmetric flux functions similar to
Eq.(11). However, more sophisticated flux functions can also be used and can also be coupled with
fluxfunctionsinEq.(1).Theuseoffluxfunctionsinthedefinitionofthedifferentialoperators,Eq.(1),
havebeenconsideredin[4,5];andthosespecifictotheFinitePointsetMethodhavebeenstudiedin
[31].
Fortime-dependentproblems,stencilcoefficientsfornumericaldifferentialoperatorsarefoundsuch
that they satisfy Eq.(10) for conservation of fluxes at the previous time level, in addition to Eq.(2).
Thus, if a field φ needs to be conserved, when proceeding from time-step tn to tn+1, the discrete
first-derivativesarefoundbythefollowingminimization
(cid:88)
ckm =∂km ∀m∈M, (12)
ij j i
j∈Si
(cid:88) ckφ(n) = 1 (cid:88)F (φ(n)), (13)
ij j V il
j∈Si i l∈I¯i
(cid:88)
minJ = W (ck)2, (14)
ij ij
j∈Si
where the bracketed superscript denotes the time-level and k =x,y,z. We emphasize that the
superscript k in Eq.(12) - Eq.(14) is used to denote only the first order derivatives in different
directions, x,y,z, whereas the superscript ∗ in Eq.(1) denotes any arbitrary differential operator.
Eq.(12)-Eq.(14)constituteanenhancementoftheGFDMdifferentialoperatorsbytheconservation
condition Eq.(13). If the monomials m∈M are taken up to second order, Eq.(12) represents 10
conditions in 3D and 6 conditions in 2D. Whereas the number of unknown c values is the same
ij
as the number of points in each support domain, which is typically around 50 and 20, in 3D and
2D respectively, for the support sizes being used. Thus, there is plenty of numerical freedom to
imposetheadditionalfluxconservationcondition,Eq.(13),withoutaffectingtheaccuracyofgradient
reconstructions, Eq.(12). Clearly, the stencil coefficients c can be found locally at each point i,
ij
withouttheneedfortheglobalsystemsusedbypreviousconservativemeshfreeGFDMs[4,5].
For each point i=1,2,...,N, Eq.(12) represents the usual monomial consistency conditions as
described by Eq.(2), which is dependent on all points in the support domain S . On the other hand,
i
Eq.(13) represents the conservation of numerical fluxes at the previous time-level, where the fluxes
aredefinedsolelyonthelocallydefinedcontrolcell,whichdependsonlyonthenearestneighboursI¯.
i
ThisdifferenceisillustratedinFigure1.Inexplicittime-integrationschemeswherespatialderivatives
are computed only at the previous time level, the differential operators defined by Eq.(12) - Eq.(14)
degeneratetoonesbasedsolelyonthelocallydefinedcontrolcell,makingitresembleVoronoi-based
FVMs. This is because the spatial derivatives used in such cases would be completely described by
Eq.(13).Thus,theinclusionofEq.(13)isonlydonewithimplicittime-integrationschemes,whichare
usedthroughoutthispaper.
The additional time taken due to the addition of Eq.(13) is not significant. As mentioned earlier,
tessellations are usually performed in classical GFDMs for prescribing volumes to points, areas to
6
FLUXCONSERVINGMESHFREEMETHOD 7
boundarypointsandthedetectionoffreesurfaces.Thus,noextratessellationsneedtobeperformedto
determinethelocalcontrolcells.Secondly,thesizeofthesystemsbeingsolvedischangedbyavery
smallamount.Forexample,in3D,forsecondorder-accuratedifferentialoperators,classicalGFDMs
solveasystemof10equationsandabout50unknownsateachpoint.TheadditionofEq.(13)changes
thattoasystemof11equationsand50unknowns.Further,whenthedifferentialoperatorsaredefined
locally as done in this paper, the largest portion of time of meshfree GFDM simulations is taken by
theiterativesolversforthelargesparselinearsystemsobtainedbythediscretizationofthePDEs.The
computationofthedifferentialoperatorstakesupamuchsmallerportionofthesimulationtime.Thus,
theadditionofthefluxconservationconstraint,Eq.(13),hasaminimaleffectontheoveralltimetaken
bytheentiresimulation.ThisisshowninthenumericalexamplesinSection5.3.
Henceforth, for the sake of brevity, we denote this variant of GFDMs, in which the differential
operators conserves fluxes of specific numerical fields, as FC-GFDM, where the FC stands for flux
conserving. We note that Eq.(13) does not ensure the conservation of numerical fluxes of the field
determined at level n+1 with respect to the given differential operator stencils. Thus, even with
respect to the field φ, the use of differential operators defined by Eq.(12)-Eq.(14) would only form
anapproximatelyconservativemethod.
ItshouldbeensuredthatthechoiceofthefluxfunctionF isdonesuchthatEq.(12)andEq.(13)are
consistentwitheachotherwhenφintheneighbourhoodofiislinearlydependentonthemonomials
m∈M.ForF definedasinEq.(11),thisisensuredbythefactthattheVoronoicellisclosedandnot
self-intersecting by definition. That guarantees that geometric fluxes are conserved. For 0 order, the
cellisclosed,
(cid:88)
(cid:126)n A =0. (15)
il il
l∈I¯i
Similarly,firstordergeometricconservationensures
(cid:40)
(cid:88)nkxk(cid:48)A = 0 ifk (cid:54)=k(cid:48), (16)
il il il V ifk =k(cid:48),
l∈I¯i i
where k,k(cid:48) =x,y,z denote the spatial dimension, (cid:126)xil is the geometric center of the edge or face i(cid:101)l
forl(cid:54)=i,and(cid:126)xii =(cid:126)xi.BythedefinitionoftheVoronoidiagram,(cid:126)xil = (cid:126)xi+2(cid:126)xl.Eq.(16)isalsousedto
determinethevolumeofthecell.
3.3 HigherOrderDerivatives
Thesameprocedureusedabovecanbeextendedtocomputenumericaldifferentialoperatorsforhigher
orderderivatives.Thisisillustratedwithadiffusionoperator.ConsiderthediffusionoperatorD such
thatDφ=∇·(α∇φ).Thedivergencetheoremappliedtothevectorfieldα∇φleadsto
(cid:90) (cid:90)
∇·(α∇φ)dV = (cid:126)n·(α∇φ)dA. (17)
Ω ∂Ω
AdiscretelocalversionofEq.(17)onthecontrolcellscanbewrittenas
V ∂˜Dφ=(cid:88)G (φ,α), (18)
i i il
l∈I¯i
where G is a numerical flux function. Throughout this paper, G is taken to be symmetric, similar to
Eq.(11).
(cid:40)
α (cid:126)n ·∂˜∇φA ifi∈∂Ωandl=i,
G (φ)= i i i i (19)
il α (cid:126)n ·∂˜∇φA elsewhere,
il il il il
where αil = αi+2αl, ∂˜i∇ =(∂˜ix;∂˜iy;∂˜iz) in 3D, and (cid:126)nil·∂˜i∇l is an approximation of the first order
derivativeinthedirectionof(cid:126)n andisgivenby
il
φ −φ
(cid:126)n ·∂˜∇φ= l i . (20)
il il (cid:107)(cid:126)x −(cid:126)x (cid:107)
l i
7
8 P.SUCHDEETAL.
Once again, this method can easily be extended to use more sophisticated flux functions and less
diffusive approximations than Eq.(20). Using Eq.(18), when proceeding from time-step tn to tn+1,
thediscretediffusionoperatorisfoundbythefollowingminimization
(cid:88)
cDm =∂Dm ∀m∈M, (21)
ij j i
j∈Si
(cid:88) cDφ(n) = 1 (cid:88)G (φ(n),α), (22)
ij j V il
j∈Si i l∈I¯i
(cid:88)
minJ = W (cD)2. (23)
ij ij
j∈Si
The FC-GFDM differential operators can be used on both fixed and moving point clouds. However,
theiruseonfixedpointcloudshasthesignificantdisadvantageofneedingdifferentialoperatorstobe
recomputed at every time-step. On the other hand, since the traditional GFDM differential operators
dependonlyonpointlocations,theycanbecomputedinasinglepre-processingstep.Thisdisadvantage
isnotpresentformovingpointclouds,sincechangingpointlocationsmeansthatdifferentialoperators
alwaysneedtoberecomputed.Wethusrestrictthenumericalstudyinthefollowingsectionstomoving
pointcloudsandLagrangianframeofreferences.
4 Advection Diffusion Equation
Considertheadvection-diffusionequationinLagrangianformulationonafixeddomainΩ⊂R2
D(cid:126)x
=(cid:126)v, (24)
Dt
Dφ
=∇·(α∇φ), (25)
Dt
whereφistheconcentrationortemperature,αisthediffusivity,and(cid:126)vistheadvectionvelocitywhich
isassumedtobedivergence-free.Energyconservation,Eq.(6),leadsto
d (cid:90) (cid:90)
φdV = (cid:126)n·(α∇φ−(cid:126)vφ) dA, (26)
dt
Ω ∂Ω
For the FC-GFDM, the numerical diffusion operator are computed as done in Eq.(21) - Eq.(23).
ThespatialdiscretizationofEq.(25)isobtainedas
Dφi =∂˜Dφ= (cid:88) cDφ . (27)
Dt i ij j
j∈Si
Inthefirststep,movementofthepointcloudisperformedinaccordancewithEq.(24),
(cid:126)v(n)−(cid:126)v(n−1)
(cid:126)x(n+1) =(cid:126)x(n)+(cid:126)v(n)∆t+ i i (∆t)2, (28)
i i i ∆t
for each point i=1,2,...,N. The explicit time-integration for the movement of points results in a
CFL-likeconditiononthetime-stepsize
(cid:18) (cid:19)
h
∆t=C . (29)
∆t (cid:107)(cid:126)v(cid:107)
min
For the advection diffusion test case being considered, we use a coefficient of C =0.01. The
∆t
movementisfollowedbyanimplicitEulertime-integrationofEq.(25)
φ(in+1)−φ(in) = (cid:88) cDφ(n+1), (30)
∆t ij j
j∈Si
8
FLUXCONSERVINGMESHFREEMETHOD 9
10-2
10-3
E
0
10-4
GFDM
FC-GFDM
10-5
10-2 10-1 100
h
Figure3.ConvergenceoferrorfortheAdvectionDiffusionequation.
Table2.Errorsandexperimentalconvergenceratesfortheadvectiondiffusiontestcase.histhesmoothinglength,
N isthenumberofpointsintheentiredomainattheinitialstate,(cid:15) istheerrorinenergyconservation,andris
E
theorderofconvergenceof(cid:15) .
E
GFDM FC-GFDM
h N (cid:15) r (cid:15) r
E E
0.4 752 2.28×10−3 − 1.24×10−3 −
0.2 2778 9.66×10−4 1.24 4.05×10−4 1.61
0.1 10565 4.76×10−4 1.02 1.61×10−4 1.33
0.05 36224 2.89×10−4 0.72 6.73×10−5 1.26
TheimplicitsystemresultingfromEq.(30)canbesolvedusinganiterativesolver.Intheresultsbelow,
we use the BiCGSTAB solver [38]. Homogeneous Neumann boundary conditions are applied for φ.
Giventheseboundaryconditions,errorinenergyconservationcanbemeasuredby
(cid:12) (cid:12)
(cid:12)(cid:82) φ(end)dV −(cid:82) φ(0)dV +(cid:82)tend(cid:2)(cid:82) (cid:126)n·(cid:126)vφdA(cid:3) dt(cid:12)
(cid:12) Ω Ω 0 ∂Ω (cid:12)
(cid:15)E = (cid:82) φ(0)dV , (31)
Ω
where the bracketed superscript end denotes the values at the end of simulation, t=t , and the
end
bracketed superscript 0 denotes the initial condition. The computational domain Ω is taken to be
[−2,2]×[−2,2]. A uniform diffusion coefficient α=0.4 and advection field(cid:126)v =(−y,x) are used.
Initialconditionsaretakentobe
(cid:40)
500 if(cid:107)(cid:126)x−(1,0)(cid:107)2 <0.1,
φ((cid:126)x,0)= (32)
0 elsewhere.
The convergence of the error in energy conservation with a varying smoothing length is shown in
Figure 3. The figure illustrates that the new flux conserving method produces much smaller errors.
Bothmethodstakesimilarsimulationtimes.Foranerror(cid:15)andconsecutivesmoothinglengthsh and
1
h ,therateofconvergenceofthesolutionwithchangingsmoothinglengthismeasuredas
2
(cid:16) (cid:17)
log (cid:15)(h2)
r = (cid:15)(h1) . (33)
(cid:16) (cid:17)
log h2
h1
Theerrorsinenergyconservationandthenumericalordersofconvergenceoftheerrorsaretabulated
inTable2.
9
10 P.SUCHDEETAL.
5 Incompressible Navier–Stokes Equations
ConsiderfluidflowmodeledbytheincompressibleNavier–StokesequationsinLagrangianform
D(cid:126)x
=(cid:126)v, (34)
Dt
∇·(cid:126)v =0, (35)
D(cid:126)v η 1
= ∆(cid:126)v− ∇p+(cid:126)g, (36)
Dt ρ ρ
where(cid:126)visthefluidvelocity,ρisthedensity,ηisthedynamicviscosityand(cid:126)gincludesbothgravitational
accelerationandbodyforces.Weconsiderconservationwithrespecttothevelocityfield,whichleads
to
(cid:90) (cid:90)
∇·(cid:126)vdV = (cid:126)n·(cid:126)vdA, (37)
Ω ∂Ω
(cid:90) (cid:90)
∇·((cid:126)v⊗(cid:126)v)dV = (cid:126)n·((cid:126)v⊗(cid:126)v)dA, (38)
Ω ∂Ω
(cid:90) (cid:90)
∆(cid:126)vdV = (cid:126)n·∇(cid:126)vdA. (39)
Ω ∂Ω
For the FC-GFDM, the numerical Laplace operator is computed in a way similar to the numerical
diffusionoperatorinSection3.3.
(cid:88)
c∆m =∂∆m ∀m∈M, (40)
ij j i
j∈Si
(cid:88) c∆vk,(n) = 1 (cid:88)G (vk,(n),1) k =x,y,z, (41)
ij j V il j
j∈Si i l∈I¯i
(cid:88)
minJ = W (c∆)2, (42)
ij ij
j∈Si
wherethefluxfunctionGisasdefinedearlierandvk,(n) denotesthek-componentofthevelocityat
j
pointj andtime-leveln.Thenumericalfirstderivativeapproximationsarecomputedas
(cid:88)
ckm =∂km ∀m∈M, (43)
ij j i
j∈Si
(cid:88) ckvk,(n) = 1 (cid:88)F (vk,(n)), (44)
ij j V il
j∈Si i l∈I¯i
(cid:88) ckvk,(n)vk(cid:48),(n) = 1 (cid:88)F (vk,(n)vk(cid:48),(n)) k(cid:48) =x,y,z, (45)
ij j j V il
j∈Si i l∈I¯i
(cid:88)
minJ = W (ck)2, (46)
ij ij
j∈Si
wherethefluxfunctionF isasdefinedearlier.WhileEq.(44)addsoneextraconstrainttothediscrete
differentialoperatordefinitions,Eq.(45)addstwoorthree,dependingonthespatialdimension.
Using the above defined numerical differential operators, Eq.(34) - Eq.(36) are solved with a
meshfree projection method, similar to that of Chorin [6]. In the first step, the point cloud is moved
by solving Eq.(34) according to Eq.(28). The projection method begins with the computation of an
intermediatevelocity(cid:126)v∗bysolvingtheconservationofmomentumequation,Eq.(36),
(cid:126)v∗−(cid:126)v(n) η 1
= ∆(cid:126)v∗− ∇p∗+(cid:126)g, (47)
∆t ρ ρ
10