Variable resolution smoothed particle hydrodynamics schemes for 2-D and 3-D viscous flows
Abstract
Smoothed Particle Hydrodynamics (SPH) is a Lagrangian particle-based method for the numerical solution of the partial differential equations that govern the motion of fluids. The main aim of this thesis work is to better enable the applicability of SPH to problems involving multi-scale fluid dynamics. In the first part of the thesis, the capability of the SPH method to simulate three-dimensional isotropic turbulence is investigated with a detailed comparison of Lagrangian and Eulerian SPH formulations. The main reason for this first investigation is to provide an assessment of the error introduced by the particle disorder on the SPH discrete operators when being purely Lagrangian. When the free decay of isotropic turbulence in a triple periodic box is studied, the Eulerian SPH formulation achieves a very good agreement with other well validated reference solutions, whereas Lagrangian SPH yields an inaccurate prediction of turbulent energy spectra. When considering linearly forced isotropic turbulence, the use of a Godunov-type SPH scheme becomes essential for the achievement of a stable solution. The efficacy of the particle shifting technique applied to turbulent SPH flows is also studied in this part of the thesis and numerical findings indicate that corrective terms derived from the arbitrary Lagrangian–Eulerian theory are essential for a proper estimation of turbulence characteristics. Subsequently, numerical analyses of a decaying isotropic turbulent flow are carried out for the first time using SPH schemes based on high-order kernels. A dramatic increase in the accuracy of the results is observed when high-order SPH is employed, especially for the description of the vorticity dynamics. Motivated by the findings of the computational investigations above, the second part of this thesis focuses on the implementation and testing of a novel SPH variable-resolution algorithm. A domain-decomposition approach is adopted to partition the computational domain into regions having different particle resolutions. Each numerical sub-problem is then closed by appending buffer regions to every sub-domain, and populating these regions with particles whose physical quantities are obtained by means of interpolations over adjacent sub-domains. These interpolations are carried out using a second-order kernel correction procedure to ensure the proper consistency and accuracy of the interpolation process. The mass transfer among sub-domains is modeled by evaluating the Eulerian mass flux at the domain boundaries. Particles that belong to a specific zone are created/destroyed in the buffer regions and do not interact with fluid particles that belong to a different resolution zone. The algorithm is implemented in the DualSPhysics open-source code [62] and optimized thanks to DualSPHysics’ parallel framework. The algorithm is tested on a series of different fluid dynamics problems: a 2-D hydrostatic tank case, a flow past cylinder for different values of the Reynolds number, a flow past an oscillating cylinder in the cross-flow direction, and the propagation of regular waves across a rectangular tank. The present algorithm is able to simulate efficiently fluid dynamics problems characterized by a wide range of spatial scales, achieving a ratio between the coarsest and the finest resolution up to a factor equal to 256.The investigation is then extended to 3-D fluid dynamics problems, such as the flow past a sphere with a Reynolds Number equal to 300 and 500, for which a SPH solution using a uniform resolution is unfeasible due to the high computational cost, showing a good agreement between the results obtained with the variable-resolution algorithm herein presented and relevant numerical investigations in the literature. The work is then concluded with the simulation and validation of a 3-D dam-breaking flow impacting a cubic obstacle.
Full text
ABSTRACT
VARIABLE RESOLUTION SMOOTHED PARTICLE
HYDRODYNAMICS SCHEMES FOR 2-D AND 3-D VISCOUS FLOWS
by
F ancesco Ricci
Smoo hed Pa icle Hyd odynamics (SPH) is a Lag angian pa icle-based me hod o
he nume ical solu ion o he pa ial di e en ial equa ions ha go e n he mo ion
o luids. The main aim o his hesis wo k is o be e enable he applicabili y
o SPH o p oblems in ol ing mul i-scale luid dynamics. In he i s pa o he
hesis, he capabili y o he SPH me hod o simula e h ee-dimensional iso opic
u bulence is in es iga ed wi h a de ailed compa ison o Lag angian and Eule ian SPH
o mula ions. The main eason o his i s in es iga ion is o p o ide an assessmen
o he e o in oduced by he pa icle diso de on he SPH disc e e ope a o s when
being pu ely Lag angian. When he ee decay o iso opic u bulence in a iple
pe iodic box is s udied, he Eule ian SPH o mula ion achie es a e y good ag eemen
wi h o he well alida ed e e ence solu ions, whe eas Lag angian SPH yields an
inaccu a e p edic ion o u bulen ene gy spec a. When conside ing linea ly o ced
iso opic u bulence, he use o a Goduno - ype SPH scheme becomes essen ial o he
achie emen o a s able solu ion. The e icacy o he pa icle shi ing echnique applied
o u bulen SPH lows is also s udied in his pa o he hesis and nume ical indings
indica e ha co ec i e e ms de i ed om he a bi a y Lag angian–Eule ian heo y
a e essen ial o a p ope es ima ion o u bulence cha ac e is ics. Subsequen ly,
nume ical analyses o a decaying iso opic u bulen low a e ca ied ou o he i s
ime using SPH schemes based on high-o de ke nels. A d ama ic inc ease in he
accu acy o he esul s is obse ed when high-o de SPH is employed, especially o
he desc ip ion o he o ici y dynamics.
Mo i a ed by he indings o he compu a ional in es iga ions abo e, he
second pa o his hesis ocuses on he implemen a ion and es ing o a no el
SPH a iable- esolu ion algo i hm. A domain-decomposi ion app oach is adop ed o
pa i ion he compu a ional domain in o egions ha ing di e en pa icle esolu ions.
Each nume ical sub-p oblem is hen closed by appending bu e egions o e e y
sub-domain, and popula ing hese egions wi h pa icles whose physical quan i ies a e
ob ained by means o in e pola ions o e adjacen sub-domains. These in e pola ions
a e ca ied ou using a second-o de ke nel co ec ion p ocedu e o ensu e he
p ope consis ency and accu acy o he in e pola ion p ocess. The mass ans e
among sub-domains is modeled by e alua ing he Eule ian mass lux a he domain
bounda ies. Pa icles ha belong o a speci ic zone a e c ea ed/des oyed in he bu e
egions and do no in e ac wi h luid pa icles ha belong o a di e en esolu ion
zone.
The algo i hm is implemen ed in he DualSPhysics open-sou ce code [62] and
op imized hanks o DualSPHysics’ pa allel amewo k. The algo i hm is es ed on
a se ies o di e en luid dynamics p oblems: a 2-D hyd os a ic ank case, a low
pas cylinde o di e en alues o he Reynolds numbe , a low pas an oscilla ing
cylinde in he c oss- low di ec ion, and he p opaga ion o egula wa es ac oss a
ec angula ank. The p esen algo i hm is able o simula e e icien ly luid dynamics
p oblems cha ac e ized by a wide ange o spa ial scales, achie ing a a io be ween
he coa ses and he ines esolu ion up o a ac o equal o 256.
The in es iga ion is hen ex ended o 3-D luid dynamics p oblems, such as he
low pas a sphe e wi h a Reynolds Numbe equal o 300 and 500, o which a SPH
solu ion using a uni o m esolu ion is un easible due o he high compu a ional cos ,
showing a good ag eemen be ween he esul s ob ained wi h he a iable- esolu ion
algo i hm he ein p esen ed and ele an nume ical in es iga ions in he li e a u e.
The wo k is hen concluded wi h he simula ion and alida ion o a 3-D dam-b eaking
low impac ing a cubic obs acle.
VARIABLE RESOLUTION SMOOTHED PARTICLE
HYDRODYNAMICS SCHEMES FOR 2-D AND 3-D VISCOUS FLOWS
by
F ancesco Ricci
A Disse a ion
Submi ed o he Facul y o
New Je sey Ins i u e o Technology
in Pa ial Ful illmen o he Requi emen s o he Deg ee o
Doc o o Philosophy in Mechanical Enginee ing
Depa men o Mechanical and Indus ial Enginee ing
Augus 2023
Copy igh ©2023 by F ancesco Ricci
ALL RIGHTS RESERVED
APPROVAL PAGE
VARIABLE RESOLUTION SMOOTHED PARTICLE
HYDRODYNAMICS SCHEMES FOR 2-D AND 3-D VISCOUS FLOWS
F ancesco Ricci
D . Angelan onio Ta uni, Disse a ion Ad iso Da e
Assis an P o esso o Mechanical Enginee ing, NJIT
D . Samaneh Fa okhi ad, Commi ee Membe Da e
Assis an P o esso o Mechanical Enginee ing, NJIT
D . Samuel Liebe , Commi ee Membe Da e
Assis an P o esso o Mechanical Enginee ing Technology, NJIT
D . Simone Ma as, Commi ee Membe Da e
Assis an P o esso o Mechanical Enginee ing, NJIT
D . Jos´e Manuel Dom´ınguez Alonso, Commi ee Membe Da e
Assis an P o esso o Applied Physics, Uni e sidade de Vigo, Ou ense, Spain
BIOGRAPHICAL SKETCH
Au ho : F ancesco Ricci
Deg ee: Doc o o Philosophy
Da e: Augus 2023
Unde g adua e and G adua e Educa ion:
•Doc o o Philosophy in Mechanical Enginee ing,
New Je sey Ins i u e o Technology, Newa k, NJ, US, 2023
•Mas e o Science in Compu a ional Fluid Dynamics,
C an ield Uni e si y, C an ield,UK, 2018
•Mas e o Science in Mechanical Enginee ing,
Poli ecnico di Ba i, Ba i, I aly 2016
•Bachelo o Science in Mechanical Enginee ing,
Poli ecnico di Ba i, Ba i, I aly, 2013
Majo : Mechanical Enginee ing
P esen a ions and Publica ions:
F., P. A.S.F. Sil a, P. Tsou sanis, . F. An oniadis, “Ho e ing o o solu ions by high-
o de me hods on uns uc u ed g ids,,” Ae ospace Science and Technology,
Volume 97, 2020.
F. Ricci, R. Vacondio, A. Ta uni; Di ec nume ical simula ion o h ee-dimensional
iso opic u bulence wi h smoo hed pa icle hyd odynamics.” Physics o
Fluids, 2023; 35 (6): 065148.
F. Ricci, R. Vacondio and A. Ta uni, “A a iable esolu ion SPH scheme based
on independen domains coupling,” P oceedings o he 17 h SPHERIC
In e na ional Wo kshop, Rhodes, 2023.
F. Ricci, R. Vacondio and . Ta uni, “High-o de SPH schemes o DNS o u bulen
lows,” P oceedings o he 2022 SPHERIC In e na ional Wo kshop, Ca ania,
2022.
i
To my g and a he and my nephew
ACKNOWLEDGMENT
Fi s and o emos , I exp ess my g a i ude o my ad iso s, P o . Angelo Ta uni and
P o . Rena o Vacondio, o allowing me o pu sue his Ph.D. p og am and o hei
cons an suppo and guidance h oughou all hese yea s.
Special hanks also o my de ense commi ee, in he pe son o P o . Fa okhi ad,
P o . Liebe , P o . Ma as and P o . Dominguez o hei eedback on my esea ch.
I would also like o hanks he inancial suppo ecei ed by Gene al Mo o s
unde G an No. GAC3794, he Na ional Science Founda ion unde G an No.
2209793 and he depa men o Mechanical Enginee ing Technology.
I’m also g a e ul o he help I ecei ed om all he esea che s o he
DualSPHysics g oup, in pa icula , he people o he SPH esea ch g oup o he
Uni e si y o Manches e o hei cons uc i e c i icism o his esea ch du ing ou
mee ings and o he ePhysLab o hei echnical suppo and o he oppo uni y o
spending a pe iod o s ay a he Uni e si y o Vigo.
Mos impo an ly, I wan o men ion my pa en s o all he sac i ices hey ha e
made o me since I was bo n and my belo ed sis e o being my con idan . Also,
special hanks o my iends in I aly o showing hei lo e and suppo e en wi h an
ocean be ween us and o he people in he US wi h which I sha ed hese yea s a
om home.
i
LIST OF FIGURES
(Con inued)
Figu e Page
4.19 Vo ici y con ou s o ω= 1,5,10,20,30 a x=−0.5 o he (a) 2nd, (b)
4 h and (c) 6 h o de schemes o he Taylo -G een Vo ex a Re=1,600.
(d) Re e ence solu ion in [253] . ...................... 79
4.20 Tu bulen ene gy spec um 2nd, 4 h and ) 6 h o de schemes o he
Taylo -G een Vo ex a Re=1,600 wi h N= 2563pa icles. Nume ical
esul s a e compa ed o he e e ence solu ion in [253]. ......... 80
5.1 (a) Example o wo sub-domains Γ1and Γ2. (b) Bu e egions ∂Γ2
1and
∂Γ1
2wi h wid hs l∂Γ2
1= 2h1and l∂Γ1
2= 2h2a e appended o hei
espec i e sub-domains. .......................... 82
5.2 Coupling p ocedu e be ween sub-domains Γ1and Γ2. Bu e pa icles
(o ange) in e pola e hei p ope ies o e he luid pa icles (ligh blue)
in he coupled subdomain. ........................ 83
5.3 A bu e pa icle ha mo es in o he luid domain is ans o med in o a
luid pa icle (p ocess 1). A luid pa icle ha en e s he bu e egion
is ans o med in o a bu e pa icle (p ocess 2). A bu e pa icle ha
mo es ou side he ex ended subdomain ∂Γj
i∪Γiis dele ed (p ocess 3). 84
5.4 Pa icle inse ion p ocedu e: a each ime s ep, he no mal mass lux
a he ou e bounda y o he sub-domain is calcula ed and added o
he mass accumula ion poin s ( ed squa es). When he mass a he
accumula ion poin s eaches he e e ence pa icle mass, new pa icles
(g een) a e c ea ed. ............................ 86
5.5 Ske ch o he egula iza ion p ocedu e o bu e pa icles. The shi ing
co ec ion is applied only in he di ec ion angen ial o he in e ace,
while neglec ed in he no mal di ec ion n. In he co ne egion, his
p ocedu e is deac i a ed. ......................... 87
5.6 Call unc ion o he DualSPHysics GPU code using a symplec ic in eg a o 92
5.7 Call unc ion o he main loop in he new mul i- esolu ion algo i hm. .96
5.8 Compu a ional domain o he hyd os a ic ank case ........... 98
5.9 Hyd os a ic ank case: (a) densi y con ou s and (b) p essu e dis ibu ion
agains he hyd os a ic solu ion a = 20s ............... 99
5.10 Compu a ional domain o 2-D low pas a ci cula cylinde ....... 100
5.11 Dimensionless p essu e o low pas a cylinde wi h (a) Re = 100 and (b)
Re = 200 .................................. 101
xiii
LIST OF FIGURES
(Con inued)
Figu e Page
5.12 Dimensionless o ici y o low pas a cylinde wi h (a) Re = 100 and (b)
Re = 200 .................................. 103
5.13 S eamlines o low pas a cylinde wi h (a) Re = 100 and (b) Re=200 .103
5.14 Time his o y o he d ag and li coe icien o a low pas cylinde wi h
Re = 100 and Re = 200 .......................... 105
5.15 Time his o y o he li coe icien o low pas a cylinde wi h (a) Re =
100 and (b) Re = 200. CLsindica es he li coe icien ob ained om a
single esolu ion simula ion while CLmbelongs o he mul i- esolu ion
simula ion. ................................. 105
5.16 Ske ch o he di e en SPH sub-domains c ea ed o esol e he low pas
a cylinde a a ious Reynolds numbe s ................. 107
5.17 D ag coe icien o low pas a cylinde wi h (a) Re = 1000, (b) Re =
3000 and (c) Re = 9500. These SPH solu ions a e compa ed agains
nume ical esul s in [118]. ......................... 108
5.18 Vo ici y con ou s o low pas a cylinde a Re = 1000 ......... 109
5.19 S eamlines o low pas a cylinde a Re = 1000 ............. 110
5.20 Vo ici y con ou s o low pas a cylinde a Re = 3000 ......... 111
5.21 S eamlines o low pas a cylinde a Re = 3000 ............. 112
5.22 Vo ici y con ou s o low pas a cylinde a Re = 9500 ......... 113
5.23 S eamlines o low pas a cylinde a Re = 9500 ............. 114
5.24 Time his o y o he li coe icien o low pas an oscilla ing cylinde
a Re = 100 o di e en ampli ude Aand equency F a ios: (a)
(A, F) = (0.25,0.9), (b) (A, F) = (0.25,0.5), (a) (A, F) = (0.25,1.5),
(a) (A, F) = (1.25,1.5) .......................... 117
5.25 Powe Spec al Densi y (PSD) o he low pas an oscilla ing cylinde a
Re=100 o di e en ampli ude Aand equency F a io ........ 118
5.26 Con ou s o dimensionless o ici y o low pas an oscilla ing cylinde a
Re = 100 wi h di e en ampli ude Aand equency F a ios ..... 119
5.27 Compa ison o SPH ae odynamic o ces agains esul s in [194] o low
pas an oscilla ing cylinde wi h A= 0.25 ................ 120
5.28 Compu a ional domain o he p opaga ion o egula wa es case .... 120
xi
LIST OF FIGURES
(Con inued)
Figu e Page
5.29 Compa ison o ee-su ace ele a ion (a, b) and o bi al eloci ies (c– )
be ween he mul i- esolu ion and uni o m esolu ion SPH simula ions
egula wa es p opaga ion ......................... 121
5.30 (a) Densi y and (b) eloci y con ou s o he mul i- esolu ion simula ion
o he wa e p opaga ion es case .................... 122
5.31 Time his o y o he mass a ia ion in he mul i- esolu ion simula ion . . 123
6.1 Longi udinal- and c oss- sec ions o he compu a ional domain o he low
pas a sphe e. ............................... 125
6.2 (a) Time his o ies and (b) no malized powe spec al densi ies o he
ae odynamic coe icien s o he low pas a sphe e a Re = 300. .... 127
6.3 (a) Time his o y and (b) no malized powe spec al densi y o he alue
o he s eamwise coe icien a a x/D = 5.75 on he wake cen e line o
he low pas a sphe e a Re = 300. ................... 128
6.4 Compa ison o he a e age and oo -mean-squa e alues o he s eamwise
eloci y along he wake cen e line wi h he nume ical esul s in [242] .128
6.5 Flow isualiza ion by he Q c i e ion, colo ed by he eloci y magni ude
U, a e e y qua e o pe iod om a iew no mal o he (x, z) plane o
he low pas a sphe e a Re = 300. ................... 130
6.6 Flow isualiza ion by he Q c i e ion, colo ed by he eloci y magni ude
U, a e e y qua e o pe iod om a iew no mal o he (x, y) plane o
he low pas a sphe e a Re = 300. ................... 131
6.7 (a) Time his o ies and (b) no malized powe spec al densi ies o he
ae odynamic coe icien s o he low pas a sphe e a Re = 500. .... 132
6.8 No malized powe spec al densi y o he s eamwise eloci y and
p essu e a (a)-(b) (x, y, z) = (2.5D, 0,0), (c)-(d) (x, y, z) = (2D, 0.3,0),
(c)-(d) (x, y, z) = (2D, 0.0,0.3), o he low pas a sphe e a Re = 500. 133
6.9 Flow isualiza ion by he Q c i e ion, colo ed by he o ici y ωx, a e e y
pe iod T1= 1/S 1 om a iew no mal o he (x, y) plane o he low
pas a sphe e a Re = 500. ........................ 135
6.10 Con igu a ion o he expe imen al se up o he h ee-dimensional dam-
b eaking es case [113]. .......................... 136
6.11 Veloci y magni ude con ou s o he h ee-dimensional dam b eak using
he mul i- esolu ion algo i hm. The e inemen egion is highligh ed in
ed. ..................................... 137
x
LIST OF FIGURES
(Con inued)
Figu e Page
6.12 Compa ison o he ime his o y o he p essu e a he on side measu emen s
gauges be ween he p esen simula ion and he expe imen al da a [113]. 138
6.13 Compa ison o he ime his o y o he p essu e a he op side measu emen s
gauges be ween he p esen simula ion and he expe imen al da a [113]. 139
6.14 Compa ison o he ime his o y o he wa e ele a ion a he measu emen s
gauges be ween he p esen simula ion and he expe imen al da a [113]. 140
6.15 Snapsho s o he densi y [kg/m3] con ou s ac oss he cen e -line sec ion in
he longi udinal di ec ion du ing he i s impac o dam-b eaking lows
agains he obs acle. ........................... 142
x i
CHAPTER 1
INTRODUCTION
1.1 Backg ound and Mo i a ion
The equa ions ha model he mo ion o an incomp essible iscous low, i.e., he
Na ie -S okes equa ions, a e cha ac e ized by non-linea e ms ha make hei
ma hema ical ea men s ill an open challenge 100 yea s a e hei i s de i a ion.
Analy ical solu ions ha e been de eloped o speci ic cases by linea izing e ms
o educing dimensionali y, p o iding quali a i e desc ip ions o low in simple
geome ies. Howe e , in o de o ob ain quan i a i e esul s, sol ing he Na ie -S okes
equa ions nume ically is c ucial.
Compu a ional Fluid Dynamics (CFD) is he compu a ional science encom-
passing all nume ical me hods used o analyze luid mo ion and se es wo main
pu poses. F om a scien i ic pe spec i e, i helps be e unde s and complex
phenomena like u bulence, mul iphase lows, and low-induced noise. In indus y, i
educes cos s associa ed wi h expe imen al aspec s o he design p ocess by na owing
he ange o a iables in expe imen al uns, op imizing designs, and sho ening he
ime om concep ual design o p oduc ion in a ious sec o s such as au omo i e,
ae ospace, and ene gy.
The ea lies bu s ill mos popula nume ical me hods in CFD a e mesh-based
me hods such as he Fini e Di e ence Me hod (FDM) [92], he Fini e Volume Me hod
(FVM) [129], and he Fini e Elemen Me hod (FEM) [284]. In mesh-based me hods,
he physical domain is disc e ized by compu a ional nodes, which a e opologically
connec ed.
Gene a ing he compu a ional g id is a complex and c ucial pa o he nume ical
compu a ion wo k low. I is o en he mos ime-consuming ask and demands
1
expe ise o ensu e he p ope placemen o compu a ional nodes, conside ing he
physics o speci ic phenomena while main aining g id quali y. High-aspec a io
o en angled g id elemen s can nega i ely impac he accu acy and s abili y o he
nume ical solu ion, pa icula ly in he p esence o in ica e geome ical ea u es [73].
Be o e applying he disc e iza ion speci ic o a pa icula nume ical echnique,
he i s choice conce n he kinema ic desc ip ion o he con inuum. Usually,
wo di e en o mula ions a e employed [64]: in he Eule ian desc ip ion, he
compu a ional nodes emain ixed in space and ime, and con ec i e e ms a e added
in o he go e ning equa ion, while he con inuum mo es and de o ms wi h espec
o he mesh, as opposed o he Lag angian desc ip ion, whe e each compu a ional
elemen is associa ed o a ma e ial pa icle hus ollowing i du ing i s mo ion. Bo h
o mula ions ha e hei ad an ages and disad an ages, depending on he pa icula
physics o he p oblem: o luid dynamics lows, he Eule ian desc ip ion is usually
p e e ed because i is concep ually simple . Howe e , i is gene ally unable o
desc ibe accu a ely mo ing in e aces. The Lag angian o mula ion is p e e ed when
dealing wi h s uc u al mechanics because i is implici ly able o ack in e aces
and deal wi h a ime-dependen s ess-s ain ela ionship. Howe e , as opposed o
he Eule ian desc ip ion, when he e is la ge ma e ial de o ma ion, he mo ion o he
compu a ional nodes can esul in dis o ed o en angled elemen s ha can de e io a e
he accu acy and he con e gence o he nume ical me hods.
Due o he successes in p edic ing ae odynamical lows, he CFD echnique
has also been ex ended o luid-s uc u e in e ac ion (FSI) p oblems, whe e mo able
o de o ming objec s in e ac wi h su ounding luid lows. This class o p oblems
is cha ac e ized by s ong non-linea i ies and a mul i-physics na u e, making he
nume ical solu ions in a single ma hema ical amewo k a he challenging. The ange
o applica ions o FSI p oblems spans di e en enginee ing a eas, including ae oelas-
2
ici y [104], bio-medical lows [96], biological lows [241], s uc u al enginee ing [124],
coas al and ma ine applica ions [278] among he o he s.
Many di e en CFD me hods ha e been p oposed o ackle FSI p oblems.
Among hose, he wo mos popula echniques a e A bi a y Lag angian-Eule ian
(ALE) [93] and he Imme se Bounda y Me hod (IBM)[193].
The idea behind he A bi a y-Lag angian-Eule ian (ALE) desc ip ion is o
e ain he ad an ages o bo h he Eule ian and he Lag angian o mula ion while
mi iga ing hei d awbacks. Usually, in ALE me hods, a body- i ed g id disc e izes
he solid domain, while an Eule ian desc ip ion is employed o he luid a
egion. Ins ead, nea he in e ace, an a bi a y eloci y is imposed on he mesh
o a oid excessi e mesh dis o ion, and con ec i e luxes accoun o his a bi a y
mo ion in he go e ning equa ions. Howe e , in he p esence o la ge displacemen ,
he e-meshing and emapping p ocess is una oidable, along wi h he associa ed
compu a ional cos and he challenges in emapping he solu ion o e he new mesh,
which can de e io a e he accu acy o he nume ical solu ion.
In he Imme sed Bounda y Me hod, non-con o mal meshes a e used o essella e
he luid and he solid egion, which a e usually desc ibed using a di e en kinema ic
desc ip ion [159]. The abili y o employ o e lapping g ids, a oid ALE me hods’
compu a ionally cumbe some e-meshing p ocess and g ea ly simpli y he g id
gene a ion p ocess. The coupling be ween he non-con o mal meshing is usually
achie ed by adding a local olume ic o ce in o he go e ning equa ions o he luid
pa , using a smoo hing unc ion o edis ibu e he e ec o he in e ace o e a
ange o compu a ional nodes. Howe e , di e en s a egies ha e been p oposed
[87]. Ano he coupling app oach is he Cu -Cell Fini e-Volume app oach [275], which
doesn’ use a olume ic o cing and has be e mass and momen um conse a ion
p ope ies. Ne e heless, wo main d awbacks cha ac e ize he IBM: he i s one,
as discussed in [25], conce ns he ”added-mass” e ec when dealing wi h a high
3
luid- o-solid densi y a io; he second one ega ds he simula ion o a high Reynolds
numbe lows, o which a mesh e inemen p ocedu e is equi ed in o de o ensu e
an adequa e esolu ion nea he solid bounda ies [254].
Ano he class o p oblems ha in ol e mo ing in e aces is ee-su ace and
mul iphase lows p oblems. In his case, wo main s a egies a e used in mesh-based
me hods, on -cap u ing and on - acking app oaches [243].
The Volume o Fluid (VOF) echnique [94] is he mos popula on -cap u ing
me hod. In he VOF, he in e ace be ween wo luids is ep esen ed by a colo
unc ion, which alue is based on he olume ac ion o each phase in a pa icula
compu a ional elemen . The me hod is composed o wo s eps: in he i s one, he
in e ace is econs uc ed om he alue o he colo unc ion in each cell. Ea ly
s udies used a Simple Line In e ace Cons uc ion (SLIC) [181, 94], ep esen ing
he in e ace by segmen s aligned wi h he mesh. Al hough i is e y simple,
his me hod leads o la ge in e ace smea ing. Mo e accu a e app oaches, such as
he Piecewise Linea In e ace Cons uc ion (PLIC) [207], signi ican ly imp o e he
me hod’s accu acy. The colo unc ion is ad ec ed using he eloci y ield in he
second pa o he app oach. Despi e he ma hema ical o mula ion ensu ing good
conse a ion p ope ies, he discon inuous econs uc ion can lead o ins abili ies in
he p esence o high-cu a u e in e aces.
Ano he app oach is he Le el-Se me hod [190], in which he in e ace is
ep esen ed by a smoo h unc ion ha mo es acco ding o an ad ec ion equa ion.
The ad an age o he Le el-se me hod is ha wi h espec o he VOF, he unc ion
is smoo h; howe e , i has a wo se conse a ion mass p ope y wi h espec o he
o me me hod.
In on - acking echniques [244], ins ead, he bounda ies be ween di e en
phases a e ep esen ed by a se o ma ke poin s ha a e connec ed and mo e along
wi h he luid. The d awbacks o hese me hods a e he addi ional da a s uc u e o
4
desc ibe he on and he explici ea men ha equi e he opology change o he
in e ace.
Opposed o g id-based me hods, in meshless-based me hods he app oxima ion
o he go e ning equa ions is buil upon a se o compu a ional nodes o which
g id connec i i y is no speci ied. Among hem, he e a e he Smoo hed Pa icle
Hyd odynamics (SPH) me hod [80], he Meshless local Pe o -Gale kin me hod [12],
he Di use Elemen Me hod [180], he Elemen -F ee Gale kin me hods [19], and he
Mo ing-pa icle semi-implici me hod [116].
The SPH me hod is a guably he mos popula o simula ing incomp essible
iscous lows. In he SPH me hod, he Na ie -S okes equa ions a e app oxima ed
upon a se o disc e e pa icles ha ca y he physical p ope ies o he luid. The
in e pola ion is based on he con olu ion wi h a ke nel smoo hing unc ion. One o
he ad an ages o he SPH me hod is i s obus ness; in ac , as demons a ed in [24],
he SPH o mula ion is consis en wi h a a ia ional app oach, ensu ing conse a ion
o all ele an physical quan i ies. Mo eo e , due o he Lag angian o mula ion, he
luid p ope ies a e ad ec ed exac ly and able o implici ly desc ibe in e aces, suach
as ee-su ace.
The SPH Esea ch and Enginee ing In e na ional Communi y (SPHERIC) is
an in e na ional o ganiza ion ha g oups he communi y o SPH esea che and
indus ial p ac i ione s. The main objec i e o SPHERIC is o s ee he esea ch
ocus on he SPH me hod. In ac , he s ee ing commi ee has iden i ied i e di e en
aspec s o SPH ha need o be add essed in o de o encou age he widesp ead
adop ion o he me hod o CFD s udy. The SPH G and Challenges a e [245]:
Con e gence, consis ency and s abili y, Bounda y condi ions, Adap i i y, Coupling o
o he models, Applicabili y o indus y. One opic ha has been poo ly add essed by
he SPH esea ch and limi s he ange o applicabili y and he ideli y o he me hod
conce n he simula ion o u bulen lows. In his wo k is s udied he issue o he
5
inclusion o u bulence e ec s in he SPH me hod. In pa icula , a majo ocus is
gi en o he ela ionship be ween u bulence modeling and he issue o adap i i y
wi hin he SPH me hod.
1.2 Aim and Objec i es
The main goal o he p esen disse a ion is o ex end he ange o applicabili y o
he Smoo hed Pa icle Hyd odynamics me hod, wi h a ocus on u bulen lows. To
his end, he p esen objec i es a e se :
1. Gain insigh in o he pe o mance o he SPH me hod in he compu a ion o
u bulen lows.
2. Analyze he e ec o he disc e iza ion e o due o he disc e e ope a o , he
densi y di usion e m and he pa icle shi ing on he nume ical compu a ion
o iso opic u bulence.
3. Assess he pe o mance o high-o de ke nel scheme in he nume ical solu ion
o iso opic u bulence
4. De elop a no el and highly-e icien a iable- esolu ion app oach o u bulence
p oblems whe e high local esolu ion is equi ed.
5. Implemen he new algo i hm in he DualSPHysics open-sou ce code and
alida ing ac oss di e en es cases.
6. Ex ension and alida ion o he app oach o he nume ical compu a ion o h ee-
dimensional lows.
1.3 S uc u e o he Disse a ion
This disse a ion is s uc u ed as ollows:
1. In Chap e 1 is p esen ed an o e iew o he nume ical app oaches o he
compu a ion o lows cha ac e ized by mo ing in e aces
2. Chap e 2 p esen a li e a u e e iew o e he s a e o he a o he SPH me hod
wi h a pa icula ocus o e he u bulence and i s modeling wi hin he SPH
me hod and he issue o adap i i y.
3. Chap e 3 p esen he ma hema ical basis o he SPH along wi h he nume ical
echnique adop ed in his wo k.
4. In Chap e 4 he nume ical compu a ion o homogenous iso opic u bulence
wi hin he SPH me hod is ad essed. The e ec o he disc e e ope a o o e
he accu acy is s udied, and he e ec o a Densi y Di usion Te m and he
Pa icle Shi ing Technique is assessed. The nume ical simula ion esul s wi h
high-o de SPH schemes o decaying iso opic u bulence a e p esen ed.
6
app oach whe e he low is cha ac e ized by complex mo ing in e ace,e.g., ee-su ace
o mo ing objec s.
As p e iously discussed, one o he p ima y limi a ions o he Smoo hed Pa icle
Hyd odynamics (SPH) me hod, which ini ially hinde ed i s b oad applica ion in
enginee ing, is i s high compu a ional cos . This is due o a la ge compu a ional
s encil ( ypically on he o de o 30+ and 300+ pa icles o 2-D and 3-D simula ions,
espec i ely) compa ed o Fini e Volume Me hod (FVM) o Fini e Elemen Me hod
(FEM). Howe e , wi h he ise o massi ely pa allel a chi ec u es, such as hose based
on G aphics P ocessing Uni s (GPUs), nume ous in-house o open-sou ce ([62, 21, 34])
codes ha e apidly de eloped.
O iginally, GPUs we e u ilized in g aphics applica ions, such as image/ ideo
p ocessing and ideo gaming. Thei a chi ec u e is buil a ound he S eaming
Mul ip ocesso uni , composed o se e al A i hme ic Logic Uni s (also known as
CUDA co es). The la es gene a ions o GPUs ha e hund eds o s eaming mul ip o-
cesso s, enabling hem o pe o m housands o a i hme ic ope a ions simul aneously.
The SPH me hod, pa icula ly in i s Weakly-Comp essible o mula ion u ilizing
an explici ime-scheme, is especially sui ed o such a chi ec u es. This is due o he
high a i hme ic in ensi y o pa icle-pa icle in e ac ion calcula ions which, as no ed
in Dominguez e al. [59], gene ally ep esen s he mos compu a ionally expensi e
ope a ion in an SPH simula ion.
In he SPH me hod, he compu a ional s encil is compu ed h ough he c ea ion
o a neighbo lis . Dominguez e al. [59] ha e discussed he p ima y echniques o
c ea ing his neighbo lis , emphasizing he impo ance o pa icle eo de ing o ensu e
op imal coalesced memo y access. Addi ionally, in Dominguez e al.[61], a ious
op imiza ions a e de ailed o u he enhance he e iciency o a GPU implemen a ion
o an SPH model.
13
2.2.2 Con e gence and accu acy o SPH me hod
The disc e iza ion e o in he SPH me hod is composed o wo con ibu ions: he
i s one is due o he smoo hing p ocedu e ha , o a con en ional smoo hing ke nel
unc ion, has an o de o O(h2), whe e his he smoo hing leng h. The second sou ce
o e o is due o he disc e e app oxima ion and is a unc ion o bo h he smoo hing
leng h and he numbe o neighbo s Nbincluded in he ke nel suppo .
Monaghan [164] conjec u ed ha he SPH me hod because he disc e e p essu e
g adien ends o a ange he pa icles in a glass-like con igu a ion, p esen s mo e
a o able con e gence p ope ies han Mon e-Ca lo me hods. In Zhu e al. [283]
is no ed ha o main ain cons an he a e o con e gence, one mus ha e h→0,
Nb→ ∞ and N→ ∞ simul aneously, and p opose a powe -law o ela e hand Nb o
he o al numbe o pa icles Nin he compu a ional domain o a ypical 2nd o de
smoo hing ke nel unc ion as:
h∝N−1/6, Nb∝N0.5(2.1)
The same conclusions we e eached in [199], whe e he disc e iza ion e o o a 1D SPH
app oxima ion was analyzed h ough a second Eule -McLau in summa ion o mula.
The immedia e consequence o hese indings is ha o p ese e he second-o de
accu acy, he compu a ional s encil, which is al eady la ge in compa ison o mesh-
based me hods, mus g ow as hsh inks, inc easing he compu a ional cos .
Besides, ea ly SPH p ac i ione s un in o he so-called ”pai ing ins abili y”
o e a ce ain h eshold o neighbo s. As demons a ed in Dehnen and Aly[55], his
phenomenon is caused by nega i e alues in he Fou ie ans o m o he smoo hing
ke nel. Fo his eason, Wendland ke nels [265] ha e become he s anda d in he SPH
me hod.
Ano he widesp ead echnique o dec ease he inaccu acies due o a diso de
dis ibu ion is he Pa icle Shi ing Technique (PST) [136], which aims o es o e a
14
mo e egula dis ibu ion by mo ing he posi ion o he pa icles, ypically modeling
he displacemen wi h a Fick’s law based on he concen a ion g adien . A simila
app oach has also been p oposed wi hin he SPH-ALE o mula ion [187].
The PST has soon become a co ne s one when dealing wi h incomp essible
iscous lows, al hough i equi es ca e ul ea men in he p esence o a ee su ace
due o he ze o- h e o in oduced in he SPH g adien ope a o by he unca ed
suppo . In his case, he o mula ion o he PST is usually modi ied by elimina ing
he no mal componen o he ee su ace o he shi ing ec o [136]. Di e en
p ocedu es ha e been p oposed in o de o iden i y he luid pa icles ha belong o
he ee-su ace: in he me hod p oposed in Lee e al.[125], he ee-su ace pa icles
a e iden i ied h ough he alue o he di e gence o he posi ion ec o . A mo e
compu a ionally cos ly bu accu a e app oach has been p oposed in Ma one e al.
[152], whe e he minimum alue o he eigen alues o he eno maliza ion ma ix
[202] and an addi ional p ocedu e, based on scanning he ”umb ella-shaped” egions
a e used o de e mine he pa icles belonging o he ee-su ace. Recen ly, di e en
wo ks ha e add essed he inconsis ency in oduced a he ee-su ace om he shi ing
algo i hm [109, 262, 145, 120].
Addi ionally, i e a i e explici and implici shi ing algo i hms [203] ha e also
been p oposed o ensu e a be e egula iza ion o he pa icle dis ibu ion.
Ano he aspec closely ela ed o he con e gence p ope y is he consis ency
o he SPH me hod, in pa icula in he p esence o bounda ies ha unca e he
suppo domain o he smoo hing ke nel. In ha case, nei he he ze o no he
i s -o de consis ency is ensu ed. Di e en nume ical echniques ha e been p oposed
o co ec he inconsis ency: app oach based on he eno maliza ion ma ix [202, 101],
Mo ing Leas -Squa e schemes [58], he Co ec i e Smoo hed Pa icle Hyd odynamics
Me hod (CSPH) [36], he Fini e Pa icle Me hod [140] (FPM), he modi ied Smoo hed
Pa icle Hyd odynamics me hod (MSPH) [16].
15
The downside o hese app oaches is he compu a ional cos associa ed wi h he
solu ion o he linea sys em, which inc eases s eeply wi h he dimensionali y and he
o de o consis ency equi ed. Besides, mos o hese me hods b eak he symme ici y
o he pa icle in e ac ion, wi h he loss o he conse a ion p ope ies o he me hod,
al hough some au ho s [185] sugges ha his p ope y can be elaxed.
Di e en s a egies ha e been p oposed in he las yea s o achie e a highe
con e gence a e. Using a Riemann-SPH scheme, in A esani e al. [13] has been
p oposed a WENO- econs uc ion whe e he polynomials a e ob ained h ough an
MLS in e pola ion. The econs uc ion s encils a e de ined by pa i ioning he
suppo domain o he ke nel in o di e en sec o s. This app oach has been u he ly
imp o ed in An ona e al. [7], whe e he FPM me hod is used in place o he MLS
scheme. This scheme has also been ex ended in A esani e al. [14] o include
a high-o de space- ime econs uc ion wi h an ADER-WENO-SPH scheme. The
d awbacks o his app oach a e ha he con e gence a e is s ill limi ed by he SPH
app oxima ion in he pa icle in e ac ion and he compu a ional cos associa ed wi h
he MLS econs uc ion ha equi es a ma ix in e sion o each s encil
Lind and S ansby [135] ha e shown ha i is possible o achie e high-o de
con e gence a es by employing high-o de ke nels. These ke nels a e ob ained by
elaxing he non-nega i e p ope y o he smoo hing ke nel. The sho coming o his
app oach is he ke nel mus be e y well sampled, es ic ing he ange o applicabili y
o a uni o m dis ibu ion o pa icles. Following his app oach in Nasa e al. [178],
high-o de Di ichle and Von Neumann bounda y condi ions a e p oposed, while in
Nasa e al. [177], a new o mula ion based on a ke nel consis ency co ec ion has
been p oposed o limi he smoo hing e o due o he disc e e SPH ope a o when
high-o de ke nels a e employed.
16
2.2.3 Weakly-Comp essible s. incomp essible SPH
The e a e wo main app oaches o ea ing incomp essible lows in SPH: he Weakly-
Comp essible SPH (WCSPH) and he Incomp essible SPH (ISPH).
The la e app oach has been p oposed i s ly in Cummins and Rudman [52],
whe e a Cho in’s p ojec ion me hod is used o en o ce an incomp essibili y o he low
h ough he solu ion o a p essu e Poisson equa ion ha a ises om he di e gence-
ee eloci y ield assump ion. Al e na i e algo i hms ha e been also p oposed in Shao
and Lo [221] and Hu and Adams[97]. No ably, i has been shown in [230, 231] ha he
ISPH and he mo ing pa icle semi-implici (MPS) me hod a e equi alen . The ISPH
has been applied o ee-su ace lows [108, 107, 228, 128], o sho e applica ion [133],
and mul i-phase low [132]. The ad an ages o he ISPH agains he WCSPH a e he
abili y o ob ain a smoo he p essu e ield, sol e he well-known p oblems ha a lic
he WCSPH, and ha is able o use a la ge ime s ep. Howe e , while he la e
app oach is compu a ionally e icien due o he explici scheme ha is mo e easily
pa allelizable, especially by exploi ing a pa allel amewo k wi h a SIMD pa adigm,
such as OpenMP and CUDA, he ISPH p esen s bo lenecks ha limi he e iciency
o he me hod. The lack o opological connec i i y be ween compu a ional nodes,
one o he dis inc i e ai s o he SPH me hod, o ces he cons uc ion o he PPE
ma ix a e e y ime s ep. Fu he mo e, because o he la ge compu a ional s encil,
he memo y equi emen s o s o ing he PPE ma ix a e conside ably la ge han he
WCSPH, and also, wi h espec o he la e o mula ion ha en o ces he ee-su ace
bounda y condi ion implici ly, he ISPH mus ely on ee-su ace de ec ion me hod
o co ec ly apply he Di ichle bounda y condi ions o close he PPE.
In he WCSPH he densi y and he p essu e a e coupled h ough a s i equa ion
o s a e, usually Tai ’s Equa ion. To a oid an excessi e es ic ion o he ime-s ep
due o he CFL condi ion, he Mach numbe is aken a ound 0.1, which bound he
17
a ia ion o he densi y wi hin he 1%. The ad an age o his app oach is he abili y
o use an explici ime-s epping scheme.
One o he d awbacks o he WCSPH o mula ion is he p esence o high-
equency oscilla ion ha a ec s he smoo hness o he densi y ield. This aspec
has been s udied by se e al au ho s [112, 68], and i s s ems om he combina ion o
wo di e en ac o s: on he one hand, he e is he employmen o a s i equa ion;
on he o he hand, he e is he colloca ion na u e o he SPH scheme.
Va ious echnique ha e been p oposed du ing he pas wo decades o add ess
his issue. Se e al au ho s [255, 192, 98] ha e p oposed an app oach based on he
de ini ion o a local Riemann p oblem o sol e he pa icle in e ac ion. The Riemann-
SPH me hod has been used wi h di e en limi e s([200, 174, 99, 210, 117], howe e , o
iolen ee-su ace low, his app oach has e ealed o in oduce oo much dissipa ion.
In Fe a i e al. [72], has been p oposed a di usi e e m ob ained applying a
Rusano Flux o he Riemann-SPH app oach, bu in oducing less dissipa ion wi h
espec o he la e app oach. One o he d awbacks o his e m is ha i is nei he
consis en , so ha doesn’ anish o h→0, no p ese e he hyd os a ic solu ion. To
es o e he consis ency in Mol eni and Colag ossi[162] is p oposed a densi y di usion
e m based on he disc e iza ion o he laplacian o he densi y ield wi h he Mo is
o mula, which has been imp o ed in An uono e al. [10], he so-called ”δ-SPH
scheme, in o de o p ese e he hyd os a ic solu ion wi h a consis ency co ec ion a
he ee-su ace based on he calcula ion o he eno maliza ion ma ix. In G een e
al. [85] is p esen ed a di usi e e m based on he applica ion o a Roe’s app oxima ing
Riemann sol e , and also i is shown ha he δ-SPH can be iewed as a pa icula
case o his model. Mo e ecen ly, in Fou akas e al. [75], a new di usi e model is
p oposed, based on he neglec ion o he hyd os a ic p essu e in he calcula ion o he
Laplacian, which is able o p ese e he hyd os a ic solu ion wi hou elying on he
calcula ion o he eno maliza ion ma ix, educing he compu a ional cos .
18
2.2.4 Bounda y condi ion in SPH
The imposi ion o bounda y condi ions o close he nume ical p oblem is a challenging
opic wi hin he SPH me hod, and i has been lis ed as one o he SPH ”G and
Challenges” by he SPHERIC commi ee. Among he easons a e he lack o he Del a
K onecke p ope y and he inconsis ency o he SPH in e pola ion in he p esence
o bounda ies ha unca e he suppo domain o he smoo hing leng h.
Fo he de ini ion o inle /ou le bounda y condi ions, he mos popula
app oaches a e based on he ex ension o he compu a ional domain by ”bu e
egions” in which he luid p ope ies a e imposed o en o ce Di ichle BC o
ex apola ed by he luid domain o he de ini ion o Neumann BC. They also
model he in low and ou low by gene a ing o dele ing luid pa icles ha en e
he compu a ional domain. Di e en s a egies ha e been p oposed based on his
app oach, mos no ably [122, 247, 69, 238].
Rega ding solid bounda y modeling, a ious app oaches ha e been p oposed
in he li e a u e, each wi h ad an ages and sho comings. In Monaghan [165], he
solid BC was imposed by means o solid pa icles ha exe ed a epulsi e o ce o e
he luid pa icles o en o ce he no-pene a ion condi ion. The epulsi e o ce was
modeled based on he Lenna d-Jones po en ial; howe e , his o mula ion was unable
o simula e a smoo h su ace, imposing an implici oughness wi h a spa ial scale equal
o he pa icle dis ance ha esul s in a diso de ed con igu a ion o he luid pa icles
close o he in e ace. A u he modi ica ion o add ess his issue was p oposed in
[171, 170] based on he de ini ion o he no mals o he solid in e ace o ensu e ha
a pa icle mo es in he pa allel di ec ion o he solid bounda y expe ience a cons an
o ce. The en o cemen o he no-slip condi ion is hen implici ly accoun ed o by
including he iscous e m in he calcula ion o he in e ac ion o ces be ween he
luid and solid pa icles.
19
A di e en s a egy, he ”Dynamic Bounda y Condi ions” me hod, has been
p oposed in Dal ymple and Knio [53], whe e ypically one laye o pa icles is placed
o desc ibe he solid bounda ies. The ad an age o his me hod is i s compu a ional
e iciency and capabili y o disc e ize complex domain because hese solid pa icles
beha e as luid pa icles when calcula ing hei densi y alue. As s udied in C espo e
al. [50], he solid pa icles exe a o ce ha depends on he dis ance and he p essu e
o he inciden luid pa icles. This app oach has been used o simula e he in e ac ion
be ween inciden wa es and coas al s udy. The sho comings o his app oach a e,
howe e , he unphysical gaps be ween luid and solid pa icles and he gene a ion o
la ge oscilla ions in he densi y ield. Besides, he no-pene a ion condi ion is no
explici ly en o ced.
In he ”Ghos Pa icle App oach” [149], as in he DBC, se e al laye s o pa icles
a e c ea ed a he beginning o he simula ion o ep esen he solid in e ace. To
gene a e hese bounda y pa icles, he solid in e ace is ep esen ed by a piecewise
linea unc ion, usually a spline, and disc e ized by a se o pa icles wi h a spacing
equal o he cha ac e is ic pa icle size o he p oblem. A e he no mal and he
angen o his in e ace a e calcula ed, he i s laye o solid pa icles and he
associa ed in e pola ion poin s in he luid domain a e c ea ed by a simul aneous
con ac ion and expansion o he solid su ace. This p ocess is epea ed ecu si ely
o ensu e a su icien numbe o laye s based on he wid h o he suppo domain o
he ke nel smoo hing unc ion. The luid p ope ies o he solid pa icles a e hen
e ie ed, in e pola ing o e he luid domain a he in e pola ion poin s de ined in
he p ocedu e ou lined abo e. In Ma one e al. [149], he in e pola ion is ca ied ou
using an MLS echnique, while in English e al. [66], whe e he mDBC is p oposed
o add ess he issue wi h he DBC, he co ec ion echnique o Liu and Liu [140] o
es o e pa icle consis ency is employed.
20
One sho coming o he Ghos Pa icle me hod is ha complex geome ical
ea u es, i.e., sha p angles, mus be ca e ully ea ed o de ine he in e pola ion poin s
in he luid domain co ec ly. Mo eo e , in he case o subme ged hin elemen s, his
app oach equi es placing su icien laye s on bo h sides o he elemen , which can
lead o an unaccep able numbe o pa icles wi hou a a iable- esolu ion app oach.
In Adami e al. [1], he p essu e is assigned based on he o ce balance a he solid
in e ace.
Ano he popula app oach is he ”mi o ing ghos pa icles”, whe e bounda y
pa icles a e gene a ed by mi o ing wi h espec o he solid in e ace, he posi ion
o luid pa icles close o he con ou s, ha ei he ca y he ield p ope ies o ob ain
hei alues h ough in e pola ion. Howe e , his app oach is mo e compu a ionally
cumbe some wi h espec o he ”Ghos Pa icles app oach” because he mi o ing
p ocedu e mus be execu ed a each ime s ep. Mo eo e , i is di icul o handle
complex 3-D geome ies.
The Vi ual Bounda y Pa icle [72] uses a di e en s a egy o ensu e he
consis ency o he SPH in e pola ion a he bounda ies. The solid in e ace is
disc e ized by a se o bounda y pa icles whose pu pose is pu ely geome ic. An
in e io luid pa icle close o he solid con ou gene a ed a se o ic i ious pa icles
using a local poin -symme y ins ead o a plane-symme y employed in he mi o ed
pa icle app oach. This p ocedu e has been u he imp o ed in [246, 76] o ensu e
ha he local s encil esembles ha o an in e io pa icle wi h ull suppo and o
ha e a be e de ini ion o he local s encil nea co ne egions.
Based on his app oach, in Fou akas e al. [75], he Local-Uni o m-S encil
bounda y condi ion has been p esen ed: o each in e io pa icle, a he beginning
o he simula ion, a uni o m s encil o i ual pa icles is c ea ed. This local s encil
mo es alongside he luid pa icles associa ed. T iangula elemen s hen disc e ize
he solid in e ace, and a each ime s ep, using a aycas ing algo i hm, he pa icles
21
o he uni o m local s encil a e iden i ied in he bounda y egion. These pa icles
a e hen ac i a ed and a e used o en o ce he solid bounda y condi ion. A uni o m
s encil ensu es a ze o- h and i s -o de consis ency.
Ano he al e na i e is based on accoun ing o he unca ed ke nel a he
bounda y using su ace in eg als. These in eg als calcula e a co ec i e ac o ha
en e s he go e ning equa ions. This app oach, called ”Semi-analy ical wall bounda y
condi ion,” was i s p oposed in Kulaseg am e al. [121] and hen de eloped in
[148, 71, 155]. Howe e , his app oach is no sui ed o an e icien implemen a ion
on GPUs.
2.3 Tu bulence and i s Modeling in SPH
2.3.1 Tu bulence and i s s a is ical modeling
The s udy o u bulence is one o he mos complex aspec s o luid dynamics. I s
s a is ical a he han ma hema ical modeling is also c ucial due o he ubiqui ousness
o u bulen lows in any eal-li e lows o p ac ical in e es .
Despi e ha he complex na u e o u bulence makes a o mal de ini ion
di icul , an a emp can be made o cha ac e ize a u bulen low by he ollowing
p ope ies:
1. chao ic
2. has a la ge and con inuous spec um o spa ial and ime scales
3. h ee dimensional
4. e y sensible o change in ini ial and bounda y condi ions
5. in e mi en
6. mixing and dissipa ion a e enhanced wi h espec o lamina lows.
I mus be s essed ha he wo ds ”chao ic” and ” andom” a e no synonymous
o highligh ha u bulence a ises om a non-linea dynamical sys em, he Na ie -
S okes equa ions, which is de e minis ic.
22
nume ical esul s wi h he log-wall o he u bulen bounda y laye . In Violeau and
Issa (2007) [257], app oaches such as he k-ϵand he Explici Algeb aic Reynolds
S ess Models (EARSM) ha e been in es iga ed o add ess he simula ion o dam-
b eaking lows. The ein, he au ho s epo a good quali a i e ag eemen agains
expe imen al esul s, especially o he EARSM model, which is mo e sui ed o s udy
egions wi h high dis o ion nea he ee su ace. Fu u e con ibu ions led o a
k-ϵmodel combined wi h he semi-analy ical wall bounda y condi ions o Weakly
Comp essible SPH (WCSPH) in Fe and e al. (2013) [71] and Incomp essible SPH
(ISPH) in Le oy e al. (2014) [128], showing good ag eemen wi h he Fini e-Volume
Me hod (FVM) o a ish-pass low. The k-ϵmodel has also been used wi h he
ISPH o in es iga e wa e o e opping [219] and b eaking [220] and mo e ecen ly
soli a y[260] and pe iodic wa es[261].
In la e yea s, a swi ch o La ge Eddy Simula ion (LES) has been obse ed,
ying o exploi he analogy be ween he LES il e ing p ocedu e and he SPH
in e pola ion. One o he pionee ing con ibu ions o LES applied o SPH can be
ound in Lo and Shao (2002) [142], whe e u bulen soli a y beach wa es ha e been
s udied wi h a ocus on u bulence du ing he b eaking phase. This model has hen
been ex ended o WCSPH by Dal ymple and Roge s (2006)[54], p o iding alida ion
agains nume ical expe imen s o cases o wa e o e opping and beach wa es in 2-D
and dam b eaking low in 3-D.
In May ho e (2015) [156], a LES app oach in SPH has been coupled wi h
he semi-analy ical wall bounda y condi ions o s udying u bulen channel low,
howe e , esul s he ein ha e shown an o e p edic ion o he s eamwise eloci y.
The au ho s poin ed ou ha his is p obably due o insu icien esolu ion when
cap u ing he o ex. Mo e ecen ly, he LES app oach has been employed in he
simula ion o u bulen open channel lows o e and wi hin na u al po ous g a el
beds [105] and in he modeling o oil spill[225] using he ISPH o mula ion.
29
In Di Mascio e al. (2017) [57] and An uono e al. (2021) [9], a new app oach
o LES modeling in SPH is p esen ed. The key idea o hese con ibu ions is o use a
ime-space il e ing p ocedu e, whe e he SPH ke nel ope a o ac s as a spa ial il e
while ime il e ing is implici ly accoun ed o wi h addi ional e ms in he go e ning
equa ions. This s a egy o e s a consis en app oach o LES while in e p e ing he
δ-SPH by Mol eni and Colag ossi (2009) [162] om a LES pe spec i e, i.e., he del a
coe icien as he de ia o ic s ain o he mean low. The δ-LES model has also been
employed o s udying g a i y wa es [158] and dam-b eak low [158].
F om a Di ec Nume ical Simula ion (DNS) pe spec i e, Robinson and Monaghan
(2012) [208] ha e a emp ed he DNS o wo-dimensional decaying u bulence,
al hough 2-D u bulence is undamen ally di e en om 3-D because o he in e se
ene gy cascade phenomena [22, 119]. No able con ibu ions ha e been made in
May ho e e al. (2015) [156], whe e a wall-bounded low is simula ed, and a lowe
bound in e ms o he numbe o pa icles pe o ex is sugges ed.
F om he p esen ed li e a u e e iew, i is e iden ha se e al c ucial aspec s
emain o be add essed conce ning u bulence modeling in SPH. The ini ial consid-
e a ion in ol es de e mining he mos sui able app oach o u bulence modeling.
While he RANS modeling appea s o be a mo e a o able choice, gi en he high
compu a ional cos associa ed wi h he SPH me hod, i seems o lack he necessa y
heo e ical ounda ions o i s applica ion in he ypical SPH domain, which includes
ee-su ace lows wi h conside able non-s a iona i y and dis o ions.
LES modeling appea s o be he mos sui able app oach, bu he e a e s ill
p oblema ic aspec s ela ed o he SPH me hod ha equi e u he in es iga ion.
The i s issue conce ns he me hod’s o de o accu acy: LES modeling demands
a high deg ee o accu acy o ensu e ha he nume ical dissipa ion in oduced by
he nume ical scheme does no comp omise he ideli y o he esul s. Howe e ,
his con lic s wi h he cu en con e gence o de o he SPH me hod, which is
30
a mos second-o de in op imal si ua ions, while unde no mal condi ions, due
o he unca ion e o caused by he i egula dis ibu ion o pa icles, i anges
be ween 1 and 2. Al hough app oaches o achie ing a con e gence o de highe han
second-o de ha e been p oposed, as p e iously discussed, hey do no ye seem o be
ma u e enough o b oade applica ion.
The second aspec pe ains o he compu a ional cos o he me hod. As
p e iously discussed, a u bulen low is cha ac e ized by a wide spec um o spa ial
and empo al scales. In g id-based me hods, his aspec is add essed by inc easing
he esolu ion in a eas whe e o ex de elopmen is expec ed, pa icula ly nea
solid su aces. The in oduc ion o a iable esolu ion in SPH encoun e s challenges
a ising om he me hod’s pa icle-based Lag angian o mula ion. In he ollowing
sec ion, hese aspec s will be discussed, and an o e iew o he s a e-o - he-a o
mul i- esolu ion in he SPH me hod will be p o ided.
2.4 Adap i i y wi hin he SPH Me hod
One weakness o he SPH app oach is ha adop ing a mul i- esolu ion o mula ion
is mo e challenging han in mesh-based me hod. The e a e se e al easons:
•As p e iously illus a ed, in o de o p ese e he con e gence p ope y, he
a io be ween he smoo hing leng h and he mean pa icle size mus be kep
a leas cons an . This means ha in o de o inc ease he nume ical accu acy,
he alue o he ke nel smoo hing leng h canno be educed wi hou inc easing
he o al numbe o pa icles in he simula ion.
•The in e ac ion be ween pa icles wi h di e en smoo hing leng hs mus be
ea ed ca e ully. Di e en app oaches ha e been p oposed in he li e a u e,
and some o hem a e de i ed ying o p ese e he a ia ional consis ency
[24, 26]. Ne e heless, he smoo hed o mula ion o he me hod p e en s a
sha p a ia ion o he pa icle size.
•The symme ici y o he smoo hing ke nel unc ion, a co ne s one o he SPH
o mula ion, implies he iso opic dis ibu ion o he compu a ional nodes, as
opposed o g id-based me hods ha can exploi he aniso opy o he low by
using di e en spacing depending on he spa ial coo dina e. A ypical example
is he in la ion laye s used o disc e ize he nea -wall egion, whe e he spa ial
esolu ion in he wall-no mal di ec ion is usually one o de o magni ude lowe
han in he s eamwise o spanwise di ec ion.
31
His o ically, ea ly a emp s o he in oduc ion o adap i i y in he SPH
me hod we e ocused on he in oduc ion o a a iable smoo hing leng h o mula ion
coupled wi h he de ini ion o egions wi h di e en pa icle sizes a he beginning
o he simula ions, o example in Bone and Paz [24] o he collapse o a ci cula
dam o e a su ace, cylind ical wa e blas [26], wedge wa e en ies [184], hea ing
cylinde and cone in a eling wa es [188][189]. Howe e , hese app oaches we e
limi ed o p oblems wi h a sho ime scale and we e un easible o cases in which
he mo ion o he compu a ional nodes highly dis o ed he ini ial con igu a ion o
pa icles.
Ins ead, la e e o s aimed o dynamically inc ease he pa icle esolu ion by
in oducing a dynamic e inemen p ocess unde he cons ain s o conse ing mass,
momen um, and angula eloci y and minimizing he e o in he es ima ion o he
densi y wi h di e en spli ing pa e ns [70, 204]. The basic concep is o minimize, in
a leas squa es sense, he e o in he SPH es ima ion o he densi y ield be ween he
o iginal pa icle dis ibu ion and he e ined one. A pa ame ic s udy on he op imal
spli ing pa e n, along wi h he op imal alue o he weigh o each pa icle size λi,
hei smoo hing leng h alue hi, and he dis ance wi h espec o he o iginal pa icle
ϵi, has been p esen ed in Vacondio e al. [248]. The dynamic spli ing p ocedu e has
been coupled in Vacondio e al. [249] o a de- e inemen p ocess based on coalescing
pai s o pa icles wi h simila sizes and ex ended o 3-D in Vacondio e al. [250]. A
simila app oach is used in Yang e al. [272] o he compu a ion o mul i-phase [273]
and ee-su ace [274] lows, whe e he spli ing c i e ion is no based on he geome ic
de ini ion o a ixed e inemen egion, bu ins ead on he posi ion wi h espec o he
ee-su ace o in e ace be ween di e en phase.
Howe e , one weakness o his app oach is he loss o compu a ional e iciency
wi h he inc ease o he e inemen a io. In ac , as discussed in Vacondio e al. [248],
in o de o minimize he e o in oduced by he spli ing p ocedu e, he alue o
32
he smoo hing leng h be ween he coa se and he ine pa icles emains almos equal,
while he pa icle size dec eases due o he spli ing p ocedu e ( ypically hexagonal o
squa e). This means ha he compu a ional s encil inc ease along wi h he esolu ion,
in oducing an impo an sou ce o ine iciency o he me hod, which is e en mo e
se ious in h ee-dimensional applica ions. Mo eo e , he coalescing p ocedu e is only
pe o med pai wise, so less equen wi h espec o he spli ing p ocedu e, which
causes an unnecessa y o e head in e ms o he numbe o pa icles. In [175] [89], his
issue is ackled by aking as ke nel smoo hing leng h, ins ead o he op imal alue
p esc ibed by he leas squa e minimiza ion, he a e age alue in he suppo domain
o he ke nel, and esul s a e p esen ed o he compu a ion o low pas blu bodies.
A di e en app oach is p oposed in Ba ca olo e al. [15]: as in he spli ing-
coalescing p ocedu e, a pa icle is spli in o ine pa icles when en e ing he e inemen
egion. Howe e , in his app oach, he o iginal pa icle, ins ead o being dele ed, is
e ained and i is ipo e ically ad ec ed. A weigh unc ion go e ns he ansi ion
be ween he wo zones o a oid p essu e discon inui ies a he in e ace. The same
app oach has been imp o ed in Chi on e al. [37] by esol ing he in e ac ion
be ween coa se and e ined pa icles wi h he de ini ion o bu e zones ha a oid
he in e ac ion be ween pa icles o di e en sizes. Using his app oach, in Sun e
al. [235], esul s a e p esen ed o lows pas bodies wi h a ious shapes, coupling
he mul i esolu ion algo i hm wi h a Tensile Ins abili y Con ol (TIC) e m, while
in Sun e al. [236] he me hod is employed o s udy wa e en y o ci cula
cylinde s. In Gao e al. [78], a block-based adap i e pa icle e inemen algo i hm
is p oposed o dynamically changing he pa icle e inemen domain coupled wi h a
egula iza ion p ocess simila Pa icle Shi ing Technique (PTS) [136] o ob ain an
iso opic dis ibu ion o he e ined pa icle in he ansi ion zone.
Howe e , one o he weaknesses o he Adap i e Pa icle Re inemen app oach
is he decoupling o he posi ion be ween coa se and ine (daugh e ) pa icles, which
33
in highly dis o ed lows, as de ailed in Chaneac e al. [35], esul s in an o e -c ea ion
o pa icles which e en ually leads o uns able solu ions.
In he domain-decomposi ion me hod, he compu a ional domain is subdi ided
in o many compu a ional sub-p oblem, and a e ad anced in ime wi h an app op ia e
coupling s a egy. An app oach based on his o mula ion has been p oposed in Bian
e al. [20]. Howe e , he esul s p esen ed only wo le els o e inemen , and he
densi y o mula ion wasn’ able o ea ee-su ace lows. In Shiba a e al. [224]
a simila s a egy based on inle /ou le bounda y condi ions is in oduced o he
coupling o di e en esolu ion zones in he con ex o he MPS me hod. Mul i-
esolu ion app oaches, mainly de o ed o Fluid-s uc u e p oblems in which di e en
esolu ions a e de ined be ween he luid and he solid phase, ha e been p oposed in
Zhang e al. [281] and Khayye e al. [110]. Howe e , hese app oaces a e limi ed
only o luid-s uc u e p oblem and allow a e y small a ia ion in he pa icle size
be ween he di e en phases.
34
CHAPTER 3
NUMERICAL FRAMEWORK
This chap e p esen s he basis o he ma hema ical ea men in he SPH me hod
alongside he compu a ional echniques ha ep esen he s a e-o - he-a o he SPH
me hodology and ha ha e been used in his wo k.
3.1 Basis o he SPH Me hodology
3.1.1 SPH con inuous in e pola ion
The ma hema ical ea men o he SPH in e pola ion me hod s a s om he
con olu ion o a ield unc ion and he Di ac’s Del a unc ion δ(x−x′):
(x) = ZΩ
(x′)δ(x−x′)dx′.(3.1)
An app oxima ion o he p e ious iden i y is ob ained by subs i u ing he δ unc ion
wi h a smoo hing ke nel unc ion W(x−x′, h):
⟨ (x)⟩=ZΩ
(x′)W(x−x′, h)dx′,(3.2)
whe e his he smoo hing leng h, a pa ame e ha de ines he size o he ke nel suppo
Ω.
The p ope ies o he ke nel smoo hing unc ion Win luence he con e gence,
accu acy, and s abili y o he SPH in e pola ion and will be discussed in he nex
sec ion.
Exp essing he de i a i e o (3.2) wi h espec o x′and using in eg a ion by
pa s, i can be de i ed he exp ession o he SPH g adien ope a o :
∇⟨ (x)⟩=ZΩ∇ (x′)W(x−x′, h)dx′−ZΩ
(x′)∇W(x−x′, h)dx′=
=Z∂Ω
(x′)W(x−x′, h)·¯
ndS −ZΩ
(x′)∇W(x−x′, h)dx′.
(3.3)
35
The i s in eg al is ob ained by applying he Gauss heo em o pass om a olume
o a su ace in eg al and, assuming ha he ke nel unc ion has compac suppo and
his is no unca ed, his e m is equal o ze o, and he ollowing iden i y is ob ained:
∇⟨ (x)⟩=−ZΩ
(x′)∇W(x−x′, h)dx′.(3.4)
Mo eo e , i he smoo hing ke nel W(x−x′, h) is an e en unc ion, he las exp ession
can be ew i en as:
⟨∇ (x)⟩=ZΩ
(x′)∇W(x′−x, h)dx′.(3.5)
whe e he di e en ial ope a o ∇is e e ed o x.
F om he e on, he b acke no a ion o iden i y he SPH app oxima ion will be
omi ed.
3.2 P ope ies o he Smoo hing Ke nel Func ion
The e is a se o p ope ies he e a e desi able o he ke nel smoo hing unc ion:
1. Uni y:
ZΩ
W(x−x′, h)dx′= 1.(3.6)
2. E en unc ion:
W(x−x′, h) = W(x′−x, h).(3.7)
3. Compac ely suppo ed:
W(x−x′, h) = 0 i |x−x′|> ah. (3.8)
4. Posi i i y:
W(x−x′, h)≥0,∀x′.(3.9)
36
5. Del a unc ion:
lim
h→0W(x−x′, h) = δ(x−x′).(3.10)
6. Mono onically dec easing
7. Smoo hness
The i s condi ion ensu es ha he SPH in e pola ion has a C0consis ency when
he suppo domain o he ke nel unc ion is a om he bounda ies and is able o
ep oduce a cons an unc ion exac ly .
Toge he wi h he equi emen ha he smoo hing ke nel is an e en unc ion,
i can be demons a ed ha he SPH in e pola ion in Equa ion (3.2) is con e ging
wi h a 2nd o de con e gence a e. In ac , expanding a ield unc ion (x′) in Taylo
se ies:
(x′) = ((x) + ′(x)(x′−x) + 1
2 ′′(x)(x′−x)2+O(x′−x)3,(3.11)
and mul iplying he p e ious equa ion by Equa ion (3.2), i is ob ained:
(x) = (x)ZΩ
W(x−x′, h)dx′− ′(x)ZΩ
(x−x′)W(x−x′, h)dx′
+1
2 ′′(x)ZΩ
(x′−x)2W(x−x′, h)dx′+ZΩO(x′−x)3dx′.
(3.12)
F om he p e ious iden i y i can be seen ha in o de o ensu e ha he
app oxima ion has a 2nd o de o accu acy, he i s wo momen Mkmus be equal
o 0:
M0=ZΩ
W(x−x′, h)dx′= 1,(3.13)
M1=ZΩ
(x−x′)W(x−x′, h)dx′= 0.(3.14)
I can be obse ed ha hese wo condi ions a e ul illed i he smoo hing ke nel
W(x−x′, h) e i ies he Uni y and he E en p ope ies in Equa ions (3.6) and (3.7).
37
Mo eo e , i can also be demons a ed ha all odd momen Mka e iden ically equal
o 0.
Rega ding he o he p ope ies: he smoo hing unc ion is chosen o ha e
compac suppo o educe he nume ical s encil and he compu a ional cos ; he
posi i i y condi ion, a he han a ma hema ical equi emen , is based on he physical
admissibili y o hyd odynamical s a es, i.e., densi y and ene gy. This condi ion also
p e en s he anishing o he e en-o de momen s o he smoo hing unc ion so ha
he SPH app oxima ion has, a mos , a 2nd-o de accu acy. The Del a unc ion in
Equa ion (3.10) p ope ies ensu e ha he smoo hing unc ion eco e s he Di ac
dis ibu ion as he smoo hing leng h h ends o ze o. The six h p ope y is based on
he assump ion ha close alue mus ha e a bigge in luence o e he SPH es ima ion,
while he las p ope y imp o es he accu acy o he in e pola ion [199].
Ano he in e es ing p ope y ha de i es om he symme ic condi ion, is ha
he ke nel unc ion can be exp essed as a unc ion o he dis ance be ween spa ial
poin s:
W(x−x′, h) = W(|x−x′|, h).(3.15)
Fu he mo e, he ke nel unc ion can be w i en in dimensionless o m as:
W(|x−x′|, h) = W(q),(3.16)
whe e qis ob ained ope a ing he ollowing change o a iable:
q=|x−x′|
h.(3.17)
The e o e, a gene al exp ession o ke nel unc ion is:
W(q) = αd
hd (q),(3.18)
38
discussed in [198], his disc e iza ion ensu e a mo e iso opic pa icle a angemen
and a oid he clumping o pa icles. Besides, as shown in [259], he se o disc e e
ope a o s chosen is consis en wi h a a ia ional o mula ion, ensu ing he p ope
conse a ion o he linea and angula momen um.
Ne e heless, some au ho s ha e p oposed al e na i e o mula ion [235, 185],
based on using bo h he an isymme ic and symme ic SPH ope a o in he
disc e iza ion o he momen um equa ion, depending on he sign o he p essu e.
Howe e , hese app oaches don’ explici ly p ese e momen um conse a ion.
3.5 Viscosi y Models
3.5.1 A i icial iscosi y
The a i icial iscosi y model has been in oduced in he SPH me hod o ensu e he
s abili y o he nume ical solu ion in he p esence o s ong shocks [164]. In his
model, a dissipa i e e m Φa, simila o he Von Neumann-Rich mye iscosi y, is
added in o he momen um equa ion:
Γa=
−αh¯cabµab
¯ρab ∇aWab,i ab · ab ≤0
0,o he wise
(3.51)
wi h µab:
µab = ab · ab
ab
(3.52)
The pa ame e αis a unable coe icien usually aken in he ange be ween 0.1−0.01,
and ¯cab and ¯ρab de ine he a e age alues o he speed o sound and he densi y.
The ad an ages o his model is ha i p ese e he conse a ion o he angula
momen um. As demons a ed in [67], he alue o αcan be ela ed o he physical
kinema ic iscosi y nu as:
ν=αhc0
2(D+ 2) (3.53)
45
whe e Dis he dimensionali y o he p oblem.
3.5.2 Lamina iscosi y
The lamina iscosi y model s a om he SPH disc e iza ion o he iscous e m
o an incomp essible low
Γ = µ∇2 .(3.54)
Ins ead o disc e izing using an SPH spa ial ope a o which would be e y sensi i e o
he pa icle diso de [256], an SPH and a ini e-di e ence i s de i a i e ope a o a e
mixed [173]in o de o ob ain he ollowing disc e iza ion o he iscous s ess enso
[142]:
(µ∇2 )a= 4mb
ν ab ·∇aWab
(ρa+ρb) 2
ab
ab.(3.55)
Wi h espec o he a i icial iscosi y model, his model de i es di ec ly om he
disc e iza ion o he iscous e m in he go e ning equa ions, al hough i doesn’
conse e he angula momen um.
3.5.3 Densi y di usion e ms
One o he d awbacks o he Weakly-Comp essible Smoo hed Pa icle Hyd odynamics
is he p esence o spu ious oscilla ion in he densi y ield a he spa ial scale o
pa icles.This issue has been ad essed by adding a di usion e m in he con inui y
equa ion. In his sec ion, he Densi y Di usion Te ms used in he p esen wo ks a e
p esen ed.
δ-SPH model [8] Mol eni and Colag ossi [162] ha e p oposed a di usi e e m
based on he disc e iza ion o he Laplacian o he densi y ield in he o m :
Da= 2hc0δX
b
ϕab
ab ·∇Wab
2
ab
Vb,(3.56)
46
whe e ϕab = (ρa−ρb) and δis a unable pa ame e ha usually is aken in he
ange 0.3−0.1. This e m is consis en , because i goes o 0 as he smoo hing leng h
h→0 and p ese es he mass conse a ion. Howe e , because o he singula i y o
he Mo is o mula in he case o unca ed suppo , i di e ges in he p esence o a
ee su ace. As emedy, in [8] he δ-SPH model has been p oposed, sugges ing he
ollowing exp ession o ϕab:
ϕab = (ρa−ρb)−1
2∇ρL
a+∇ρL
b· ba, (3.57)
wi h he eno malized densi y ∇ρL
agi en by:
∇ρL
a=X
b
(ρa−ρb)La∇aWab.(3.58)
The eno maliza ion ma ix is de ined La:
La="X
b
(xb−xa)⊗Wab#−1
,(3.59)
and es o e he i s -o de consis ency also in he p esence o bounda ies o ee-su ace
ha unca e he suppo . Ne e hless, i equi e he in e sion o a ma ix 2x2 in
2-D and 3x3 in 3-D o e e y pa icles in he compu a ional domain, so he addi ional
compu a ional cos is no negligible.
G een e al. DDT [86] A di e en app oach o s abilize he densi y ield has
been p oposed in [84] and is used in his wo k. The basic idea behind his model is
o apply a Goduno SPH scheme o he con inui y equa ion and employing a Roe’s
app oxima e sol e o he Riemann p oblem a he in e ace o each neighbo ing
pa icle. A limi e unc ion is hen applied o alues o he densi y ob ained in his
ashion, so o p ese e he mono onici y o he scheme and a oid oscilla ions in he
densi y ield. The di usi e e m can ben w i en as:
Da= 2c0X
b
ϕab
||xa−xb||
∂W
∂q Vb.(3.60)
47
ϕab is now de ined as:
ϕab =Bab (ρa−ρb)−1
2(ξab∇ρC
a+ξba∇ρC
b)(xa−xb),(3.61)
whe e:
Bab =||xa−xb||
h
2ρa
ρ
a+ρ
b
.(3.62)
Wi h espec o Equa ion (3.57), he uning pa ame e δis now absen and he
magni ude o he di usion in oduced in he go e ning equa ions is now au oma ically
adjus ed by he de ini ion o econs uc ed densi y alues:
ρ
a=ρa+1
2ξ(η)∇ρC
ab ·(xb−xa),(3.63)
whe e ∇ρC
ab is a co ec ed SPH app oxima ion o he densi y g adien and ξ(η) is he
limi e unc ion, wi h ηde ined as:
η=∇ρab ·(xa−xb)
ρb−ρa
.(3.64)
Following he indica ions in [84], he limi e unc ion chosen in his wo k is he Van
Albada [252]:
ξ(η) = η2+η
η2+ 1.(3.65)
I can be seen ha , assuming ξab =ξba = 1.0 and no icing ha Bab ≈0.5, his
o mula ion is e y simila o Equa ion 3.57.
Fou akas e al. DDT [75] Fou akas e al. [75] ha e p oposed a new densi y
di usion model ha is able o main ain he hyd os a ic solu ion while a oiding he
compu a ion o he eno malized densi y g adien . The main concep is o compu e
he Laplacian o he densi y ield by neglec ing he hyd os a ic con ibu ion. In ac :
Da=δhc0X
b
ψab ·∇aWabVb(3.66)
48
whe e ψab is equal o:
ψab = 2(ρT
ba −ρH
ab)xab
||xab||,(3.67)
whe e he supe sc ip Tand H e e s o he o al and o he hyd os a ic pa o he
densi y. The hyd os a ic densi y di e ence can be ob ained by:
ρH
ab =ρ0
γ
sPH
ab + 1
CB−1
(3.68)
whe e PH
ab is he hyd os a ic p essu e di e ence:
PH
ab =ρ0gzab (3.69)
and CB=c2
0ρ0
γ.
3.6 Pa icle Shi ing
An impo an issue o he SPH me hod is he o ma ion o oid egions, especially in
he p esence o s ong o ex s uc u es, ha can a ec he s abili y and accu acy o
he nume ical solu ion.
This p oblem has been add essed in [134], whe e he Pa icle Shi ing Technique
(PST) has been p oposed: he basic concep is o shi he posi ion o he pa icles
o ensu e a mo e uni o m dis ibu ion by modelling he shi ing δxapplied o he
pa icle posi ion wi h a Fickian law based on he g adien o he concen a ion C, in
he o m:
δx=−D∇C(3.70)
whe e Dis he di usion coe icien de ined as:
D=Ah2.(3.71)
49
He e, Ais a dimensionless cons an ha is uned based on he pa icula p oblem,
and his he smoo hing leng h.
The g adien o concen a ion ∇Cis calcula ed using an SPH disc e iza ion o
he g adien gi en by:
∇C=X
b
mb
ρb∇aWab.(3.72)
In he p esence o a ee su ace, he SPH disc e iza ion in Equa ion (3.72) is
inaccu a e due o he unca ed suppo , esul ing in an upwa d mo ion o he
pa icles om he bulk o he low. To add ess his issue, in [134], only he angen ial
componen is e ained a he ee su ace while he pa icle shi ing is neglec ed in
he no mal di ec ion:
δx=
0∇· ≤AF SM ,
(¯
¯
I−¯
n⊗¯
n)δx AF ST ≤ ∇· ≤AF SM ,
δx ∇· ≥AF ST ,
(3.73)
In he p e ious exp ession, ¯
¯
Iis he second o de iden i y enso and ¯
nis he no mal
a he ee su aces es ima ed as:
¯
n=−∇C
||∇C||.(3.74)
Fo de ec ing he ee-su ace in e ace, he me hod p oposed by [125], based on he
pa icle posi ion di e gence, is used:
∇· a=X
b
Vb ab ·∇aWab (3.75)
wi h AFSM and AFST equal espec i ely o 1.1 and 1.7 in 2-D and 2.1 and 2.8 in 3-D.
50
ALE co ec i e e ms The shi ing o mula ion wi h he ALE co ec i e e ms is
gi en by [237]:
dρa
d =X
bρa
ρb
mb( ab +δ ab)·∇Wab +mb
ρb
(ρaδ a+ρbδ b)·∇Wab−Da,(3.76)
d a
d =X
bhmbPb+Pa
ρaρb∇Wab + ( a⊗δ a+ b⊗δ b)·∇Wab+
+ a(δ a−δ b)·∇Wabi+ Γa,
(3.77)
dxa
d = a+δ a.(3.78)
3.7 Solid Bounda y T ea men
In his wo k a e used he modi ied Dynamic Bounda y Condi ion (mDBC), p oposed
in [66], o he de ini ion o solid bounda y condi ions. This app oach is a modi ica ion
o he Dynamic Bounda y Condi ion (DBC) o add ess he issue o he unphysical gap
be ween he dummy pa icles and he luid pa icles. In he DBC o mula ion, a se
o dummy pa icles ha ep esen he solid in e ace a e placed in he compu a ional
domain. Thei eloci y is se o ze o, while hei densi y is e ol ed by using he SPH
con inui y equa ion. In his way, when luid pa icles come close o he solid pa icles,
he densi y and he p essu e o he o me inc ease, gene a ing a epulsi e o ce o e
he luid pa icles. Howe e , as discussed in [60], his epulsi e mechanism c ea es a
gap in he o de o he smoo hing leng h so ha he e is an inco ec de ini ion o
he solid in e ace ha ei he mus be conside ed p io o placing he solid bounda y
pa icles, o a e when he in he pos -p ocessing o he esul s.
In he mDBC, as in he Ghos Pa icle app oach [149], o each bounda y
pa icle, a ghos node is c ea ed by mi o ing he posi ion o he solid bounda y
along he solid- luid in e ace (Figu e 3.2). Once he posi ion o he ghos node is
de ined, he densi y and i s g adien a e compu ed a he ghos node adop ion a
co ec ed SPH ope a o [140]. This ope a o consis in he solu ion op he ollowing
51
Figu e 3.2 Mi o ing o ghos nodes (c osses) and he ke nel adius a ound he ghos
nodes o bounda y pa icles in a la su ace (a) and a co ne (c). Fluid pa icles
(pink) included in he ke nel sum a ound ghos nodes o bounda y pa icles in a la
su ace (b) and a co ne (d) [66].
52
linea sys em
A =b
bm=X
b
bwmVb
Amn =X
b
wm nVb
wm=Wgb Wx
gb Wy
gb Wz
gb
n=1xgb ygb zgb
=ρgρx
gρy
gρz
g
(3.79)
whe e ρgand ρi
ga e he densi y and i s g adien a he ghos node g.
Then, he densi y a he bounda y pa icle is ound h ough:
ρb=ρg+ ( b− g)·ρgρx
gρy
gρz
g(3.80)
In case he ma ix Ais ill-condi ioned, due o an insu icien numbe o pa icles,
ins ead o Equa ion (3.79), he densi y is calcula ed wi h a Shepa d co ec ion:
ρg=PbρgWgbVb
PbWgbVb
(3.81)
Fo he eloci y ub, a e en o ce no-slip condi ion by an symme ic e lec ion along he
solid in e ace:
ub= 2us−ug(3.82)
whe e usis he eloci y o he solid in e ace, and ugis he eloci y a he ghos node
calcula ed as:
ug=PbugWgbVb
PbWgbVb
.(3.83)
53
3.8 Time S epping Scheme
In his wo k p esen ed la e , a second-o de symplec ic p edic o -co ec o [127] is
employed as ime scheme:
ρn+1
2
a=ρn
a+∆
2Mn
a,
n+1
2
a= n+∆
2Fn
a,
xn+1
2
a=xn+∆
2 n
a,
ρn+1
a=ρn
a2−ϵ
2 + ϵ;ϵ=−∆ Mn
a
ρa1
2
,
n+1
a= n+∆
2Fn+1
2
a,
xn+1
a=xn+∆
2( n
a+ n+1
a) + δx,
(3.84)
whe e δx is he pa icle shi ing ob ained wi h Equa ion (3.70).
The ime-s ep is chosen acco ding o he ollowing CFL condi ion:
∆ =CCFL min(∆ ,∆ c ),(3.85)
whe e:
∆ = min
a h
Fa
,(3.86)
∆ c = min
a
h
c0+ maxb|h( b− a)·(xb−xa)
(xb−xa)2|.(3.87)
54
ini ial pa icle dis ibu ion gene a ed by a p elimina y simula ion o he same es
case. This app oach is simila o he one p oposed by Colag ossi e al. [44]. This
ini ializa ion, shown in Figu e 4.3(b), allows o bypass he a o emen ioned p oblems
and o ob ain a much smoo he eloci y ield [Figu e 4.3(d),4.3( )]. Speci ically, when
using a andom-like ini ial dis ibu ion o e a Ca esian one, a smalle dissipa ion
is obse ed a he beginning o he simula ion as shown in Figu e 4.4(a). This
di e ence be ween he wo ini ializa ions seems o lose signi icance wi h inc easing
alues o he ke nel suppo . Howe e , using a la ge ke nel suppo leads o a la ge
smoo hing e o o he SPH spa ial ope a o s, he e o e his could be de imen al
in he la e s ages when u bulen s uc u es become much smalle and a e a
isk o being smoo hed ou a i icially. To e i y his and o ul ima ely choose an
op imal ini ial se up, a sensi i i y analysis is pe o med o assess he e ec o he
smoo hing-leng h- o-pa icle-spacing a io, h/dp, on he nume ical solu ion. The
esul s a e depic ed in Figu e 4.4(b) and indica e ha alues o h/dp > 1.5 p oduce
wo se esul s o a gi en dp, especially o > 5 when he ini ial o ices a e b oken
in o smalle s uc u es. Fo his eason, a alue o h/dp = 1.5 is chosen and kep
ixed in all simula ions he ea e . Mo eo e , o sa is y he weakly comp essibili y
assump ion, he speed o sound is chosen o be c0= 10 · max, limi ing he densi y
a ia ions o 1% o he e e ence densi y.
4.3.2 Eule ian SPH
Figu e 4.5(a) shows he his o y o he kine ic ene gy o Re = 1600 wi h he Eule ian
SPH model dep i ed o any di usi e e ms in he con inui y equa ion, i.e., Equa ion
(3.49), and o 643, 1283, 2563and 5123compu a ional nodes. I is ema ked ha
o his Reynolds numbe , a ough es ima ion o he a io be ween he Kolmogo o
leng h [Equa ion (2.3)] scale and he e e ence leng h gi es L/η ≈250, indica ing
ha he wo ines esolu ions adop ed in his s udy a e able o esol e he smalles
61
(a) (b)
Figu e 4.4 (a) E ec o he ini ial pa icle dis ibu ion on he decay o he kine ic
ene gy o he 3-D Taylo -G een Vo ex a Re = 1600 wi h wo alues o he
smoo hing-leng h- o-pa icle-spacing a io, i.e., h/dp = 1.5 and 2.0. (b) E ec o he
smoo hing-leng h- o-pa icle-spacing a io on he decay o he kine ic ene gy o he
3-D Taylo -G een Vo ex a Re = 1600 wi h a andom-like pa icle ini ial dis ibu ion.
scales o mo ion. Resul s a e in good ag eemen wi h he high-o de pseudo-spec al
solu ion by Van Rees e al. (2011) [253], especially o he ines esolu ion, o which
he g id independence is almos achie ed. In Figu e 4.5(b), he ime his o y o he
ens ophy is epo ed o he same se o simula ions. I is possible o obse e ha
he nume ical solu ion is con e ging owa ds he e e ence solu ion wi h he inc ease
o he esolu ion.
4.3.3 Lag angian SPH
The same es case as in he p e ious sec ion is simula ed he e using Lag angian
SPH wi hou any di usi e e ms added o he con inui y equa ion. Looking a he
kine ic ene gy decay in Figu e 4.6(a), he solu ion is also con e ging, howe e , an
excessi e dissipa ion is obse ed when he p ocess o o ex oll-up s a s, i.e., o
3≤ ≤6. A compa ison wi h he p e ious Eule ian SPH esul s sugges s ha his
o e -di usion could be due o he diso de in he SPH pa icle dis ibu ion, which
leads o low accu acy o he SPH spa ial ope a o s. This o e -di usion causes he
62
(a) (b)
Figu e 4.5 (a) Time his o y o he kine ic ene gy and (b) ime his o y o he
ens ophy o he 3-D Taylo -G een Vo ex a Re = 1600 simula ed wi h Eule ian
SPH using 643, 1283, 2563and 5123pa icles. SPH esul s a e compa ed o he
e e ence solu ion in Van Rees e al. (2011) [253].
ea ly b eakdown o he ini ial o ices, and a disc epancy wi h he e e ence solu ion
is hus obse ed [Figu e 4.6(a)]. As can be seen in Figu e 4.6(b), his also causes a
poo ag eemen be ween SPH and he DNS solu ion in Van Rees e al. (2011) [253]
o wha conce ns he ime e olu ion o he ens ophy. When compa ing he Eule ian
and Lag angian esul s om Figu es 4.5(b) and 4.6(b), espec i ely, he la e se e ely
unde es ima es he peak alue, ema king he lowe o de o accu acy o Lag angian
SPH due o he poo e pa icle dis ibu ion.
4.3.4 In luence o he dissipa ion model
Addi ional simula ions o he same 3-D Taylo -G een Vo ex case a e p esen ed he e
wi h he adop ion o he densi y di usion e m (DDT) by G een e al. (2019) [84]
in he Lag angian SPH simula ions. In Figu e 4.7(a), he solu ions wi h and wi hou
he DDT a e compa ed o 5123pa icles in he compu a ional domain, e ealing a
negligible e ec o he di usion e m om he s andpoin o he ime his o y o he
low kine ic ene gy. Howe e , when looking a he ime e olu ion o he ens ophy
[Figu e 4.7(b)], he addi ion o he G een e al. (2019) [84] DDT in he con inui y
63
(a) (b)
Figu e 4.6 (a) His o y o he kine ic ene gy and (b) o he ens ophy o he Taylo -
G een Vo ex a Re = 1600 o esolu ions o 643, 1283, 2563and 5123pa icles wi h
he Lag angian SPH o mula ion. Nume ical esul s a e compa ed o he e e ence
solu ion in Van Rees e al. (2011) [253].
equa ion is able o yield a highe peak alue wi h espec o he s anda d Lag angian
SPH, likely due o he less noisy densi y ield which leads o an imp o emen o he
accu acy o his es case.
This is con i med when looking a he con ou s o he densi y ield in Figu e
4.8(a), whe e i can be clea ly seen ha when using a Lag angian SPH model
wi hou any dissipa ion, he densi y ield appea s domina ed by s ong nume ical
noise. Ins ead, when he G een e al. (2019) [84] DDT is enabled [Figu e 4.8(b)], a
much smoo he densi y ield is ob ained. This e ec can also explain he di e ences
be ween he u bulen ene gy spec a displayed in Figu e 4.9, whe e o wa e numbe s
in he ange o he ke nel adius size, he Lag angian SPH wi hou any DDT shows a
spu ious la ening, which is ins ead signi ican ly educed when he DDT is enabled.
Mo eo e , he Eule ian SPH is capable o compu ing he spec um co ec ly ac oss
he whole ange o equencies.
64
(a) (b)
Figu e 4.7 (a) Time his o y o he kine ic ene gy and (b) ime his o y o he
ens ophy o he 3-D Taylo -G een Vo ex a Re = 1600 simula ed wi h 5123pa icles
and wi h a Lag angian SPH o mula ion wi h and wi hou he G een e al. (2019)
[84] di usi e e m. Nume ical esul s a e compa ed agains he e e ence solu ion in
Van Rees e al. (2011) [253].
(a) (b)
Figu e 4.8 Densi y ield a = 7 o 3-D Taylo -G een Vo ex a Re = 1600: (a)
Lag angian SPH wi h no dissipa ion; (b) Lag angian SPH wi h he G een e al. (2019)
[84] di usi e e m enabled.
65
Figu e 4.9 Tu bulen ene gy spec a o he 3-D Taylo -G een Vo ex a Re = 1600
wi h 5123pa icles in he domain and wi h Eule ian SPH, and Lag angian SPH wi h
and wi hou he G een e al. (2019) [84] dissipa ion e m. Nume ical esul s a e
compa ed agains he e e ence solu ion in Van Rees e al. (2011) [253].
66
4.4 Fo ced Iso opic Tu bulence
The pe o mance o he abo e models is e alua ed he e o an iso opic u bulen low
in a iple pe iodic box wi h a linea o cing e m added o he momen um equa ion:
d
d =−∇P
ρ+ν∇2 +A ,(4.7)
whe e Ais he o cing cons an . This pa icula es case is chosen o assess
he s abili y and obus ness o he abo e schemes when simula ing long pe iods
o physical ime, which has been possible o sol e mainly hanks o he highly
pa allel pe o mance o he DualSPHysics code [62]. When compa ed o a mo e
adi ional band-limi ed o cing, a linea o cing injec s ene gy a all scales o mo ion
bu is capable o achie ing a s a iona y s a e wi h an ene gy spec um as good as
band-limi ed o cing [212]. A solenoidal eloci y ield peaking a a wa enumbe o
k0= 2 has been chosen as he ini ial condi ion [212]. As discussed by Rosales and
Mene eau (2005) [212], he ini ial eloci y ield has no in luence o e he s a iona y
solu ion which depends only on he Acons an o Equa ion (4.7) and on he size do he
domain. As o he p e ious es case, simula ions ha e been pe o med wi h bo h he
Lag angian and Eule ian SPH schemes wi h and wi hou he densi y di usion e m
o Equa ion (3.60). A ough es ima ion o he kolmogo o scales gi es L/η ≈200, so
a esolu ion wi h 2563is chosen. A lis o he di e en simula ions ha ha e been un
oge he wi h de ails abou he SPH o mula ion, densi y di usion scheme and alue
o he coe icien A used a e summa ized in Table 4.1. A esolu ion o 2563pa icles
is used o all cases and he luid kinema ic iscosi y is se o ν= 9.9471 ×10−5
4.4.1 Role o he densi y di usion e m
Upon ackling he o ced iso opic u bulence p oblem wi h an Eule ian SPH scheme
wi hou any densi y di usion e m, simula ion esul s ha e shown o be uns able due
o he insu gence o s ong oscilla ions in he densi y ield [Figu e 4.10(a)]. This
67
Table 4.1 Simula ion pa ame e s o he Fo ced Iso opic Tu bulence P oblem
Case SPH Fo mula ion DDT A
1a Eule ian G een e al. (2019)[84] 0.1
2a Lag angian G een e al. (2019)[84] 0.1
3a Eule ian G een e al. (2019)[84] 0.3
4a Lag angian G een e al. (2019)[84] 0.3
4b Lag angian G een e al. (2019)[84] + Shi .
ALE
0.3
4c Lag angian G een e al. (2019)[84] + Shi . 0.3
nume ical noise was no de ec ed in he p e ious TGV case, and i is p obably
due o he longe simula ion imes which b ing o ligh he s abili y issues o
a cen e ed-colloca ed nume ical scheme. To o e come his issue, wo di e en
dissipa ion models ha e been included in he con inui y equa ion and in es iga ed:
he δ-SPH by An uono e al. (2010)[8] and he G een e al. (2019) [84] models.
The con ou s o he densi y ield ob ained wi h he δ-SPH model[8] wi h δ= 0.1 a e
depic ed in Figu e 4.10(b): a highe le el o noise can be seen when compa ed o he
solu ion ob ained wi h he G een e al. (2019)[84] model [Figu e 4.10(c)]. This is
also highligh ed by he p obabili y dis ibu ions o densi y alues in Figu e 4.11(a),
while he densi y ene gy spec a o hese h ee solu ions in Figu e 4.11(b) indica e
ha δ-SPH is unable o ep oduce he low-mid ange scales. Al hough his issue
could be pa ially mi iga ed by inc easing he alue o he δpa ame e acco ding o
he magni ude o he o cing e m, he G een e al. (2019) [84] model appea s o be
he be e choice due o i s abili y o au oma ically adjus he amoun o in oduced
dissipa ion based on local low condi ions. Fo his eason, his model is included in
he con inui y equa ion o all Eule ian and Lag angian SPH simula ions hence o h.
68
(a) (b) (c)
Figu e 4.10 Con ou s o he densi y ield a = 9 o he o ced iso opic u bulen
cases wi h A= 0.1: (a) Eule ian SPH; (b) Eule ian SPH wi h he An uono e
al. (2010)[8] dissipa ion model; (c) Eule ian SPH wi h he G een e al. (2019)[84]
dissipa ion model.
(a) (b)
Figu e 4.11 Fo ced iso opic u bulen p oblem: (a) P obabili y dis ibu ion o he
densi y ield, and (b) ene gy spec a a = 9 o Eule ian SPH, Lag angian SPH wi h
he An uono e al. (2010)[8] dissipa ion model and Lag angian SPH wi h he G een
e al. (2019)[84] dissipa ion model.
69
(a) (b)
Figu e 4.12 (a) Time his o y o 2
ms, and (b) ime his o y o ϵ/3 2
ms o Eule ian
SPH esul s o he o ced iso opic u bulen p oblem (case 1a).
4.4.2 Di e ences be ween Eule ian and Lag angian app oaches
Two global quan i ies a e de ined o in es iga e he s a iona y beha io o he
p oblem: he mean squa e alue o he eloci y luc ua ions, 2
ms =⟨ · ⟩/3, and
he eddy u no e ime, τ= 2
ms/ϵ, whe e ϵ=−ν⟨ ·∇2 ⟩is he mean dissipa ion
a e. As can be seen in Figu e 4.12(a), 2
ms eaches a s a iona y alue o /τ > 5 up
o 20 eddy u no e imes in all he di e en cases. F om he kine ic ene gy balance
o a s eady s a e p oblem, i is known ha he a io ϵ/3 2
ms mus lie a ound he
alue o he pa ame e Aspeci ied in he go e ning equa ions. This is co obo a ed
by he ime his o y o ϵ/3 2
ms in Figu e 4.12(b), whe e he co ec alue o A= 1 is
e ie ed.
In Figu es 4.13(a) and 4.13(b), he esul s o he same se o simula ions a e
shown when Lag angian SPH is employed. I can be seen ha 2
ms is sligh ly less
oscilla o y o he Lag angian scheme. This no wi hs anding, he di e ences be ween
he wo app oaches a e no signi ican and he inal alues o he mean squa e alue o
he eloci y luc ua ions a e almos iden ical when he s a iona y solu ion is eached.
As wi h he Eule ian SPH, also Lag angian SPH p edic s he co ec alue o ϵ/3 2
ms
well, wi h a s able esul o /τ > 8 Figu e 4.13(b)).
70
(a)
(b)
(c)
Figu e 4.17 His o y o he kine ic ene gy o he (a) 2nd, (b) 4 h and (c) 6 h o de
schemes o he Taylo -G een Vo ex a Re=1,600. Nume ical esul s a e compa ed
o he e e ence solu ion in [253].
77
(a) (b)
(c) (d)
(e) ( )
Figu e 4.18 Time e olu ion o he ens ophy and kine ic ene gy dissipa ion a e o
he (a)-(b) 2nd, (c)-(d) 4 h and (e)-( ) 6 h o de schemes o he Taylo -G een Vo ex
a Re=1,600. Nume ical esul s a e compa ed o he e e ence solu ion in [253].
78
(a) (b)
(c) (d)
Figu e 4.19 Vo ici y con ou s o ω= 1,5,10,20,30 a x=−0.5 o he (a) 2nd,
(b) 4 h and (c) 6 h o de schemes o he Taylo -G een Vo ex a Re=1,600. (d)
Re e ence solu ion in [253] .
79
Figu e 4.20 Tu bulen ene gy spec um 2nd, 4 h and ) 6 h o de schemes o he
Taylo -G een Vo ex a Re=1,600 wi h N= 2563pa icles. Nume ical esul s a e
compa ed o he e e ence solu ion in [253].
80
CHAPTER 5
MULTI-RESOLUTION ALGORITHM
In he p esen chap e , an adap i e esolu ion algo i hm o SPH is p esen ed.
The scheme p esen ed he ein is based on he decomposi ion o he compu a ional
domain in o di e en sub-domains, each wi h i s own cha ac e is ic pa icle size
and smoo hing leng h. A each ime s ep, he compu a ional p oblem is sol ed in
e e y sub-domain independen ly, and he sub-domain closu e is p o ided by a bu e
egion ha ac s as a Di ichle bounda y condi ion. Coupling be ween he di e en
sub-domains is ob ained by in e pola ing he physical quan i ies in he bu e egion
using in o ma ion a ailable a he luid pa icles lying in adjacen domains.
5.1 Mul i-Resolu ion Algo i hm
The main idea behind he a iable esolu ion algo i hm o his s udy is based on a
decomposi ion o he compu a ional domain, Γ, in o a se o Nsub-domains, Γiwi h
i= 1, . . . , N , such ha Γ = SN
i=1 Γi[Figu e 5.1(a)]. Each sub-domain is cha ac e ized
by i s own cha ac e is ic pa icle size, dpi, and smoo hing leng h, hi. To sol e he
compu a ional p oblem om n o n+1, a closu e Di ichle bounda y condi ion mus
be p o ided a he bounda ies o each sub-domain. To his end, each sub-domain
Γiis ex ended by appending a bu e egion, ∂Γi, wi h wid h l∂Γi= 2hiand such
ha ∂Γi∈Γjand ∂Γj∈Γi. The la e condi ion allows es ablishing a bijec i e
opological connec ion be ween wo sub-domains, o malized wi h he no a ion ∂Γj
i,
i.e., he bu e egion o sub-domain iis coupled wi h he sub-domain j.
E e y bu e egion is popula ed by special SPH pa icles ollowing a s a egy
simila o he one adop ed in [238] o open bounda y condi ions. The “bu e ”
pa icles a e dis inc om egula luid pa icles as hey a e no upda ed using he
go e ning equa ions o mo ion. Ins ead, a he beginning o each ime s ep, a bu e
81
(a) (b)
Figu e 5.1 (a) Example o wo sub-domains Γ1and Γ2. (b) Bu e egions ∂Γ2
1
and ∂Γ1
2wi h wid hs l∂Γ2
1= 2h1and l∂Γ1
2= 2h2a e appended o hei espec i e
sub-domains.
pa icle a∈∂Γj
iob ains i s physical p ope ies om a spa ial in e pola ion o e luid
pa icles in Γj, as shown in Figu e 5.2. Bu e pa icles loca ed close o he in e ace
ha e a unca ed suppo , he e o e he co ec ed SPH in e pola ion p oposed in [139]
is employed o ensu e consis ency.
The p ocedu e o es o ing pa icle consis ency p oposed in [139] s a s om
a mul i-dimensional Taylo se ies expansion o a ield unc ion (x) mul iplied by he
ke nel unc ion and i s de i a i e up o he desi ed o de o consis ency. In his wo k,
a second-o de consis ency condi ion is en o ced, which in wo dimensions esul s in
82
Figu e 5.2 Coupling p ocedu e be ween sub-domains Γ1and Γ2. Bu e pa icles
(o ange) in e pola e hei p ope ies o e he luid pa icles (ligh blue) in he coupled
subdomain.
he ollowing linea sys em o he gene ic bu e pa icle a∈∂Γj
i:
A =b(5.1)
bm=X
b∈Γj
bwmVb(5.2)
Amn =X
b∈Γj
wm nVb(5.3)
w=Wab Wx
ab Wy
ab Wxx
ab Wxy
ab Wyy
ab (5.4)
=1xab zab
x2
ab
2xabyab
y2
ab
2(5.5)
= a x
a y
a xx
a xy
a yy
a(5.6)
whe e Wis he ke nel smoo hing unc ion, xab = [(xa−xb),(ya−yb)] is he in e -
pa icle dis ance, and Vis he olume. Once he p ope ies a e ob ained h ough he
83
Figu e 5.3 A bu e pa icle ha mo es in o he luid domain is ans o med in o a
luid pa icle (p ocess 1). A luid pa icle ha en e s he bu e egion is ans o med
in o a bu e pa icle (p ocess 2). A bu e pa icle ha mo es ou side he ex ended
subdomain ∂Γj
i∪Γiis dele ed (p ocess 3).
in e pola ion p ocedu e, he posi ion o he bu e pa icle ais upda ed in ime using
a simple Eule ime in eg a ion scheme:
xn+1
a= ∆ n
a.(5.7)
whe e xn+1
ais he posi ion o he pa icles a ime-s ep n+1, and n
ais he eloci y
a n.
To handle he exchange o mass among he di e en sub-domains, he ollowing
condi ions a e checked a he beginning o each ime s ep:
1. I a bu e pa icle mo es in o he luid domain, i is con e ed o a luid pa icle
(Figu e 5.3, p ocess 1).
2. I a luid pa icle en e s he bu e egion, i becomes a bu e pa icle (Figu e
5.3, p ocess 2).
3. I a bu e pa icle mo es ou side he ex ended subdomain ∂Γj
i∪Γi, his pa icle
is dele ed (Figu e 5.3, p ocess 3).
84
Figu e 5.4 depic s he p ocedu e o pa icle inse ion in o a gi en sub-domain. The
ex e nal bounda y o he gene ic sub-domain Γiis di ided in o mass segmen s o
leng h dpiand he Eule ian mass lux ac oss he segmen s is calcula ed a each ime
s ep using he mid-poin ule:
˙ma= max(0,−ρma( ma− b )·ndpid ) (5.8)
whe e ˙mais he mass lux du ing he ime s ep d ,ρmaand maa e he densi y
and he eloci y calcula ed a he mass accumula ion poin using he same co ec ed
in e pola ion employed o bu e pa icles, and b is he eloci y o he in e ace in
case o a mo ing subdomain. A e ˙mais added o ma, i ma≥(dp)D/ρ0, whe e
Dis he numbe o p oblem dimensions, a pa icle is inse ed in he bu e egion
a a dis ance dp/2 in he no mal di ec ion o he in e ace, as shown in Figu e 5.4.
In p esence o ee-su ace, his p ocedu e has o be adjus ed o p e en a non-ze o
mass lux a he accumula ion poin abo e he wa e le el. The ollowing o mula ion
based on he ee-su ace de ec ion me hod p oposed in [125] is employed o he
compu a ion o he mass- lux a he in e ace:
˙ma=
max(0,−ρma( ma− b )·ndpid ) i ∇· ≥ ∇· h
0 i ∇· <∇· h;
(5.9)
whe e ∇· =Pb
mb
ρb·∇Wab and ∇· h is a h eshold alue equal o 1.5 in 2-D and
2.75 in 3-D.
The e o e, because bu e pa icles mo e wi h a Lag angian eloci y which is no
ob ained om he go e ning equa ions di ec ly bu om an in e pola ion o e luid
pa icles in adjacen sub-domains, hey do no bene i om he egula iza ion e ec o
he p essu e g adien disc e iza ion [198]. These issues can yield an i egula pa icle
dis ibu ion in he bu e egions, which e en ually a ec s he co ec en o cemen
o bounda y condi ions o he sub-domains. A solu ion o his p oblem is ound
85
Figu e 5.4 Pa icle inse ion p ocedu e: a each ime s ep, he no mal mass lux
a he ou e bounda y o he sub-domain is calcula ed and added o he mass
accumula ion poin s ( ed squa es). When he mass a he accumula ion poin s eaches
he e e ence pa icle mass, new pa icles (g een) a e c ea ed.
by using he Pa icle Shi ing Technique (PST) in he bu e egions, wi h special
a en ion o a oid shi ing he bu e pa icles owa d he edges o he sub-domain.
The la e isk is mi iga ed by su ounding he edges o he bu e egions wi h laye s
o “ ixed pa icles” which ha e a ixed posi ion in space and cons an densi y ρ0.
O he han hese wo physical a iables, ixed pa icles do no ha e any o he physical
quan i y associa ed and hey in e ac wi h bu e pa icles only wi h he sole pu pose
o compu ing he shi ing co ec ion. In con as o he shi ing o mula ion adop ed
o luid pa icles, bu e pa icles a e shi ed only in he di ec ion angen ial o he
in e ace be ween he bu e and he luid egion o a oid inaccu acies when compu ing
he mass lux a he in e ace. Mo eo e , because no mal ec o s a e ill-de ined in he
co ne egions, he shi ing is disabled in hese a eas (Figu e 5.5).
Algo i hm 1 shows a pseudo-code o he p oposed mul i- esolu ion s a egy
including all he salien s eps o he simula ion.
86
Table 5.3 Main Loop Func ions and hei Desc ip ion
Func ion Desc ip ion
P edic o s ep.
In e ac ion Fo ces Call o pa icle in e ac ion (PI).
P eIn e ac ion Fo ces P epa es a iables and assigns memo y o PI.
In e ac ion Fo ces Compu es pa icle in e ac ion.
D Va iable Compu es he alue o he new a iable ime s ep.
Compu eSymplec icP e Compu es Sys em Upda e using
PosIn e ac ion Fo ces Memo y elease o a ays in GPU.
RunCellDi ide Gene a es neighbou lis .
CORRECTOR Co ec o s ep.
In e ac ion Fo ces Call o pa icle in e ac ion.
P eIn e ac ion Fo ces P epa es a iables and assigns memo y o PI.
In e ac ion Fo ces Compu es pa icle in e ac ion.
D Va iable Compu es he alue o he new a iable ime s ep.
Compu eSymplec icCo Compu es sys em upda e using
symplec ic=co ec o .
PosIn e ac ion Fo ces Memo y elease o a ays in GPU.
RunCellDi ide Gene a es Neighbou Lis .
FinishRun Shows and s o es inal o e iew o execu ion.
o di e en esea ch g oups, he DualSPHysics package is cons an ly upda ed wi h
new unc ionali ies.
To keep he ange o applicabili y o DualSPHysics in ac and possibly o ex end
i , he a iable esolu ion algo i hm was implemen ed, a o ing he p ese a ion o
93
he base s uc u e o he code, a oiding majo changes ha could lead o long- e m
main enance issues.
The e a e wo obse a ions ha can be made by looking a he algo i hm
p oposed in Sec ion 5.1 o a iable esolu ion and o he DualSPHysics code s uc u e
ou lined p e iously:
•Each compu a ional sub-p oblem is dependen on he o he ones only when he
physical p ope ies o he bu e pa icles a e in e pola ed
•To each compu a ional sub-p oblem co esponds a JSphSingle class ins ance.
Conside ing hese wo aspec s, he implemen a ion o he new algo i hm has
been s uc u ed by ope a ing a he highes le el o abs ac ion. A new class,
named JSphGpuW appe , has been implemen ed: in his new objec , an a ay o
JSphGpuSingle objec s is alloca ed, and a each objec co esponds a nume ical sub-
domain. In his way, each compu a ional sub-p oblem can be managed independen ly
and synch onized o he o he when needed.
•The numbe o sub-domains is ead om he con igu a ion ile
•The JSphGpuSingle a ay is alloca ed
•Each JSphGpuSingle is ini ialized and con igu a ed
•The main loop o ad ancing he simula ion is execu ed.
•Sub ou ine o comple ing he simula ion.
The new main loop is s uc u ed in he same way as in he s anda d DualSPHysics
implemen a ion. S ill, a he end o he ime-s ep, he ou ines o upda ing he bu e
egions a e execu ed, ollowing he pseudocode ou lined in Algo i hm 1:
1. Fo each subdomain, he solu ion is ad anced in ime by calling he same
sub ou ines as in he o iginal DualSPHysics code.
94
Table 5.4 Main loop Func ions and hei Desc ip ion
Func ion Desc ip ion
In e ac ion Bu e Ex apFlux Call o compu ing he lux a he accumu-
la ion poin .
Compu eS ep P ocedu e o ans o ma ion be ween luid
and bu e pa icles.
Bu e Lis C ea e De ini ion o he poin in which bu e pa icles
mus be c ea ed.
Bu e C ea eNewPa C ea ion o bu e pa icles.
RunCellDi ide Gene a es neighbou lis .
In e ac ion Bu e Ex apFlux Call o in e pola ing he bu e pa icles.
2. Each sub-domain compu es he coupling pa o he new mul i- esolu ion
algo i hm.
A isualiza ion o he s uc u e o he main loop o he mul i- esolu ion implemen-
a ion is shown in Figu e 5.7, while in Table 5.4 a e desc ibed he unc ions o he
coupling sec ion o he mul i- esolu ion algo i hm is di ided in o hese main s eps,
which can be di ided in h ee main s eps:
1. Compu a ion o he lux a he accumula ion poin .
2. E alua ion o he mass segmen and pa icle c ea ion and dele ion p ocedu e.
3. Pa icle eo de ing and in e pola ion o ob aining he physical p ope ies o he
bu e pa icles.
The main ad an age o he p esen implemen a ion s a egy is ha he main s uc u e
o he DualSPHysics code is le unchanged as he changes in oduced o he a iable
95
Figu e 5.7 Call unc ion o he main loop in he new mul i- esolu ion algo i hm.
96
Table 5.5 New Files Added o he Sou ce Code
Files Desc ip ion
JSphGpuW appe .cpp/.h Implemen s he class JSphGpuW appe .
JSphGpuSingleBu e .cpp De ine he hos unc ion o he coupling
algo i hm.
JSph Gpu Bu e .cu De ine he Cuda ke nel unc ion o he
coupling algo i hm.
JSphBu e .cpp/.h Implemen s he class JSphBu e .
JSphBu e Zone.cpp/.h Implemen s he class JSphBu e Zone.
In e ac ion Bu e Ex apFlux Call o in e pola ing he bu e pa icles.
esolu ion algo i hm a e addi i e. In Table 5.5 a e lis ed he new iles in oduced in
he sou ce code. Two new classes a e in oduced in o he code:
JSphBu e Zone In his class, he geome ical de ini ion o he bu e egion is
de ined.
JSphBu e De ine an a ay o JSphBu e Zone objec s, one o each bu e egion
o he sub-domain. I also implemen s he ou ine needed o e ie e he lis o
pa icles in he bu e egion and label hem as bu e pa icles.
5.3 Resul s and Discussion
5.3.1 Hyd os a ic ank
Despi e being a simple es case, he hyd os a ic ank is o en e y di icul o he
SPH me hod, because o he inconsis ency due o he p esence o a ee-su ace and
he ank wall ha in oduce nume ical e o in he solu ion. Mo eo e , is an ideal es
case o e i y he accu acy o he coupling be ween sub-domain and any inaccu acies
in oduced by he in e pola ion o bu e egions.
97
Figu e 5.8 Compu a ional domain o he hyd os a ic ank case
Fo hese easons, an hyd os a ic ank is chosen as a s udy case o alida e he
p oposed a iable- esolu ion algo i hm. The compu a ional domain consis s o a 2-D
ec angula ank, wi h wid h L= 1m, illed wi h wa e up o a heigh H= 1m.
The wa e densi y is se equal o ρ0=ρ∞= 1000 kg/m3and he g a i y g=
9.81m/s2. Rega ding he nume ical pa ame e s, an a i icial iscosi y model is used
wi h α= 0.01, while he speed o sound is equal o c0= 10 max, whe e max =√gH.
The Fou akas DDT in Equa ion 3.66 is used o o de o p ese e he hyd os a ic
solu ion.
Two e inemen zone a e employed, he coa se one wi h a esolu ion equal o
H/∆x= 50, and he ines one wi h H/∆x= 100. The e inemen egion is a squa e
box wi h leng h h1= 0.5m, cen e ed a he midline o he ec angula ank. In Figu e
5.8 he geome ical de ini ion o he compu a ional domain and o he e inemen
egions is shown.
98
(a) (b)
Figu e 5.9 Hyd os a ic ank case: (a) densi y con ou s and (b) p essu e dis ibu ion
agains he hyd os a ic solu ion a = 20s
Figu e 5.9(b) epo s he p essu e dis ibu ion agains he e ical coo dina e
a = 20s: as can be seen, he p esen mul i- esolu ion app oach is able o p ese e
he hyd os a ic solu ion, and no discon inui ies a e isible a y/H = 0.5, whe e he
ho izon al in e ace be ween he coa se and he ine esolu ion is placed. Mo eo e ,
despi e a low alue o he a i icial iscosi y coe icien α, in Figu e 5.9(a) he densi y
con ou s a = 20s e eal a egula dis ibu ion o he pa icles, ha doesn’ show
any spu ious mo ion, especially a he in e ace be ween he sub-domains.
5.3.2 Flow pas a ixed, ci cula cylinde
Flow pas a ci cula cylinde is simula ed o se e al Reynolds numbe s (Re) as a
i s s udy case o benchma k he new mul i- esolu ion algo i hm agains consolida ed
li e a u e esul s. This p oblem has been in es iga ed ex ensi ely bo h nume ically
and expe imen ally, hus i is sui able o demons a ing he ad an ages o he p esen
mul i- esolu ion app oach wi h espec o using a uni o m SPH pa icle esolu ion.
Since he Reynolds numbe based on he smoo hing leng h, Reh=U∞hν−1, mus
be in he o de o 1 o esol e he low nea he cylinde p ope ly, his es case is
99
Figu e 5.10 Compu a ional domain o 2-D low pas a ci cula cylinde
qui e challenging o esol e wi h a uni o m pa icle esolu ion, e en a low o mode a e
Reynolds numbe s, due o he p ohibi i e numbe o pa icles equi ed o he domain
disc e iza ion. Mo eo e , he le el o esolu ion needed in he a - ield is subs an ially
lowe han a ound he cylinde and hus he p oposed mul i- esolu ion algo i hm is
adop ed o educe he a e age pa icle esolu ion in he a - ield wi hou a ec ing he
global accu acy o he simula ion.
The compu a ional domain shown in Figu e 5.10 is chosen ollowing he wo k
in [238], whe e a cylinde wi h diame e D= 0.1 m is cen e ed a (x, y) = (0,0)
and he o e all domain dimensions a e se o 25D×20D o minimize he change o
any blockage e ec . A no-slip solid bounda y condi ion is applied o he cylinde ,
while ee-slip condi ions a e de ined a he bo om and uppe walls. In low-ou low
condi ions a e imposed using he o mula ion in [238]: a he inle , he eloci y and
100
(a) Re=100 (b) Re=200
Figu e 5.11 Dimensionless p essu e o low pas a cylinde wi h (a) Re = 100 and
(b) Re = 200
he densi y a e p esc ibed, whe eas a he ou le , he eloci y is ex apola ed om
he luid o he bu e egion, while he densi y is p esc ibed o he e e ence alue.
The luid is ini ialized wi h U∞(x, y) = (1,0) m/s while he densi y has an ini ial
alue equal o ρ0=ρ∞= 1000 kg/m3. The Reynolds numbe Re=U∞Dν−1is a ied
by changing he alue o he kinema ic iscosi y ν. The smoo hing leng h is se o
h= 2∆x, cons an ac oss he di e en sub-domains, and he DDT in Equa ion (3.66)
is ac i a ed in o de o di use he oscilla ions a ec ing he densi y ield due o he
cen e ed colloca ed SPH scheme. Finally, PST is applied o a oid he c ea ion o
emp y egions due o he o ical s uc u es ha de elop in he wake.
Fo he se up o he di e en esolu ion sub-domains, a minimum esolu ion
co esponding o D/∆x= 25 is chosen o he a ield. Then, new sub-domains a e
nes ed inside he a - ield, each wi h a esolu ion ha doubles as hey ge close o he
cylinde . The numbe o subdomains c ea ed depends on he low Reynolds numbe ,
wi h highe Re equi ing highe o e all esolu ion, he e o e mo e sub-domains. Table
5.6 summa izes he numbe o sub-domains used and hei espec i e dimensions and
pa icle esolu ions as a unc ion o he Reynolds numbe .
101
Table 5.6 Numbe o Sub-Domains Used and hei Respec i e Dimensions and
Pa icle Resolu ion as a Func ion o he Reynolds Numbe
Case Re Numbe o zones D/∆xmax
1 100 3 100
2 200 3 100
3 1000
4 200
5 400
6 800
4 3000
5 400
6 800
7 1600
5 9500
6 800
7 1600
8 3200
Re = 100 and 200 Fo he i s wo Reynolds numbe s Re = 100,200, only
h ee sub-domains a e u ilized as shown p e iously, and hese e inemen egions
a e cen e ed a he cylinde and ex ended downs eam o be e esol e he wake
egion. The con ou s o dimensionless p essu e P⋆=P(x/D, y/D)ρ−1U2
∞shown
in Figu e 5.11 depic he o ices’ co es, clea ly isible as low-p essu e egions
de eloping in he wake. These o ical s uc u es ha e highe in ensi y wi h inc easing
Reynolds numbe . No discon inui ies a e isible h ough he in e ace be ween
he ine and he medium e inemen egions, demons a ing he obus ness o he
coupling p ocedu e be ween sub-domains. Fu he mo e, he dimensionless o ici y
ζ⋆=ζ(x/D, y/D)DU−1can be obse ed in Figu e 5.12, colo ed such ha clockwise
and coun e -clockwise luid o a ions end owa ds blue and ed, espec i ely. The
ypical on Ka man s ee associa ed wi h he pe iodic o ex shedding is cap u ed
co ec ly and he o ical s uc u es c oss he in e ace o he di e en esolu ion
102
(a) ∗= 1 (b) ∗= 2
(c) ∗= 3 (d) ∗= 4
(e) ∗= 5 ( ) ∗= 6
Figu e 5.18 Vo ici y con ou s o low pas a cylinde a Re = 1000
109
(a) ∗= 1 (b) ∗= 2
(c) ∗= 3 (d) ∗= 4
(e) ∗= 5 ( ) ∗= 6
Figu e 5.19 S eamlines o low pas a cylinde a Re = 1000
110
(a) ∗= 1 (b) ∗= 2
(c) ∗= 3 (d) ∗= 4
(e) ∗= 5 ( ) ∗= 6
Figu e 5.20 Vo ici y con ou s o low pas a cylinde a Re = 3000
111
(a) ∗= 1 (b) ∗= 2
(c) ∗= 3 (d) ∗= 4
(e) ∗= 5 ( ) ∗= 6
Figu e 5.21 S eamlines o low pas a cylinde a Re = 3000
112
(a) ∗= 1 (b) ∗= 2
(c) ∗= 3 (d) ∗= 4
(e) ∗= 5 ( ) ∗= 6
Figu e 5.22 Vo ici y con ou s o low pas a cylinde a Re = 9500
A ound ∗≈2 [Figu es 5.22(b) and 5.23(b)], he p ima y o ex pai de aches om
he cylinde , as can be seen also by looking a he ime his o y o he d ag coe icien
113
(a) ∗= 1 (b) ∗= 2
(c) ∗= 3 (d) ∗= 4
(e) ∗= 5 ( ) ∗= 6
Figu e 5.23 S eamlines o low pas a cylinde a Re = 9500
114
in Figu e 5.17(c). Meanwhile, a new o ex pai is o med a a highe angle, and
he same in e play be ween p ima y and seconda y o ici y is obse ed, inc easing
he p essu e d ag up o ∗= 3. In e es ingly, a ∗= 4, he second o ex pai
me ges wi h he p e ious, p ima y o ices a e being sepa a ed om he body.
The e olu ion a la e imes is domina ed again by he complex in e ac ion among
o ices in he bounda y laye , esul ing in a highly uns eady low. The con e gence
s udy pe o med o his las Reynolds numbe has esul ed in using 6, 7, and 8
di e en esolu ion zones, eaching a maximum esolu ion nea he cylinde equal o
D/∆x= 800,1600,3200, espec i ely. The quan i a i e e ec o hese h ee di e en
esolu ions can be assessed in Figu e 5.17(c), whe e he ime his o y o he d ag
coe icien is shown. While all h ee esolu ions cap u e he d ag coe icien well un il
∗= 4, only he simula ion wi h D/∆x= 3200 is a good ma ch wi h he e e ence
solu ion o ∗>4, when he low is cha ac e ized by a high deg ee o uns eadiness.
Mo eo e , wi h he highes esolu ion adop ed, he low emains almos pe ec ly
symme ical as is also he case o he e e ence solu ion in [118].
Con a y o wha has been shown o Re = 100 and 200 in Sec ion 5.3.2, no
compa isons be ween mul i- esolu ion and uni o m esolu ion SPH simula ions a e
shown he e o Re = 1000, 3000 and 9500. This is mainly due o he p ohibi i e cos o
unning uni o m esolu ion simula ions o he h ee Reynolds numbe s selec ed when
he highes esolu ion has o be used. Fo example, he mul i- esolu ion simula ion o
Re = 9,500 wi h D/∆x= 3200 equi es 8 sub-domains, esul ing in a o al numbe
o SPH pa icles deployed o N≈10 ×106. Simila ly, a uni o m esolu ion SPH
simula ion wi h his le el o pa icle spacing would equi e N≈9×109pa icles,
which can be achie ed only wi h sophis ica ed memo y-dis ibu ed pa alleliza ion
[63].
115
5.3.3 Flow pas an oscilla ing, ci cula cylinde
The low pas a ci cula cylinde oscilla ing in a ans e se di ec ion o he s eamwise
di ec ion is in es iga ed o Re = 100 o di e en alues o he oscilla ion ampli ude
and equency. This es case gene a es complica ed lows and i has been in es iga ed
wi h di e en echniques, including expe imen s compu a ional s udies [194] and
expe imen al in es iga ions [32, 33, 268, 5]. He ein, low pas an oscilla ing cylinde
has been chosen o assess he capabili y o he p oposed mul i- esolu ion scheme o
simula e lows wi h mo ing bounda ies. The compu a ional se up is he same as he
one employed o s udying he low pas a ixed cylinde case (see Sec ion 5.3.2) and
he domain has been disc e ized wi h h ee di e en esolu ion zones, wi h pa icle
size equal o D/∆x= 25,50,100 (see Figu e 5.10).
A sinusoidal mo ion y( ) is applied o he cylinde in he c oss- low di ec ion:
y( ) = ADsin(2πF s ) (5.11)
whe e A=ymax/D, wi h ymax equal o he maximum displacemen and F= 0/ s,
whe e sis he equency o he o ex shedding when he cylinde is ixed. The
e inemen egions also mo e acco ding o Equa ion (5.11), so ha he ela i e
posi ion o he cylinde wi h espec o he e inemen a eas does no change in ime.
Fou di e en con igu a ions o ampli ude and equency ha e been conside ed which
can be ound in Table 5.8.
In Figu e 5.24(a), he ime his o y o he li coe icien CLis shown o [A, F] =
[0.25,0.9]. This is de ined as a “locked con igu a ion” because he o ex shedding
p esen s a dominan equency equal o 0, also co obo a ed by he Powe Spec al
Densi y (PSD) o he li coe icien ime signal shown in Figu e 5.25(a), whe e only
one peak is isible.
Fo an unlocked case, ins ead, he li signal con ains wo o mo e equencies,
and his is achie ed by changing he equency a io o F= 0.5 and F= 1.5. In he
116
Table 5.8 Ampli ude and F equency Values Chosen o he Simula ion o Flow Pas
an Oscilla ing Cylinde
Case A F
1 0.25 0.9
2 0.25 0.5
3 0.25 1.5
4 1.25 1.5
(a) (A, F ) = (0.25,0.9) (b) (A, F ) = (0.25,0.5)
(c) (A, F ) = (0.25,1.5) (d) (A, F ) = (1.25,1.5)
Figu e 5.24 Time his o y o he li coe icien o low pas an oscilla ing cylinde a
Re = 100 o di e en ampli ude Aand equency F a ios: (a) (A, F) = (0.25,0.9),
(b) (A, F) = (0.25,0.5), (a) (A, F) = (0.25,1.5), (a) (A, F) = (1.25,1.5)
i s case, he pe iodic signal is cha ac e ized by wo modes wi h di e en magni udes
and pe iod b= 2 0, whe e 0is he pe iod o he imposed sinusoidal mo ion. This is
117
(a) (A, F ) = (0.25,0.9) (b) (A, F ) = (0.25,0.5)
(c) (A, F ) = (0.25,1.5) (d) (A, F ) = (1.25,1.5)
Figu e 5.25 Powe Spec al Densi y (PSD) o he low pas an oscilla ing cylinde
a Re=100 o di e en ampli ude Aand equency F a io
con i med by looking a he PSD in Figu e 5.25(b), whe e wo peaks associa ed wi h
= 0and = 2 0can be obse ed in ag eemen wi h indings in [194]. The same
bea ing phenomena [194] is also p esen o he F= 1.5 case, al hough in his case
he li coe icien has a pe iod b= 8 0.
The las case simula es a mo e challenging con igu a ion ob ained by using
(A, F) = (1.25,1.5). Figu es 5.24(d) and 5.25(d) highligh a dominan equency
= 0coupled wi h a second = 2 0and hi d = 3 0 equency, in ag eemen wi h
he esul s epo ed by [194]. These small seconda y equencies a e no associa ed
wi h bea ing phenomena, a he hey in luence he o ex shedding pa e ns, esul ing
118
Figu e 6.1 Longi udinal- and c oss- sec ions o he compu a ional domain o he
low pas a sphe e.
he sphe e. A he same ime, in he downs eam di ec ion, he compu a ional
domain is ex ended by 20Din o de o esol e a leas 3 e ical s uc u es. A squa e
c oss-sec ion wi h size 10Dis used, which esul s in a blockage a io app oxima ely
equal o 0.7%. The de ini ion o he bounda y condi ions ollows he se up ollowed in
subec ion 5.3.2: a no-slip bounda y condi ion is applied o he sphe e, while ee-slip
condi ions a e de ined a he c oss-sec ion walls. Inle -ou le bounda y condi ions
and he ini ial condi ions a e imposed again in he same manne as in sec ion 5.3.2.
Howe e , di e en om he wo-dimensional s udy, he alue o he smoo hing leng h
is h= 1.5∆xin o de o educe he compu a ional cos . The DDT in Equa ion 3.66
and he PST a e applied o he en i e domain.
Fo he se up o he e inemen egions, a minimum esolu ion equal o D/∆x=
6.25 is chosen. As o he wo-dimensional s udy on he low pas a cylinde , new
sub-domains a e nes ed inside he a ield, each wi h a esolu ion ha doubles as
hey ge close o he cylinde , wi h a maximum esolu ion equal o D/∆x= 100,
which co esponds o a a io be ween he bounda y laye δbl and he pa icles size
equal δbl/∆x≈6−5 o a ange Re = 300 −500. I is no ed ha he o al numbe
o pa icles in hese simula ions is Np ≈11 ×106, while wi h a uni o m esolu ion,
he numbe o pa icles would ha e been equal o Np= 3 ×109.
125
Table 6.1 Nume ical Resul s o he D ag Coe icien CD, Li Coe icien CLand
S ouhal Numbe S Values o he Flow Pas a Sphe e a Re=300
Case CD∆CDCL∆CLS
P esen 0.667 0.0032 0.066 0.017 0.134
Cos an inescu and Squi es [46] 0.655 0.0032 0.065 – 0.136
Johnson and Pa el [102] 0.656 0.0035 0.069 0.016 0.138
Tomboulides and O szag [242] 0.671 0.0028 – – 0.136
Ploumhans e al. [195] 0.683 0.0025 0.061 0.014 0.13
6.1.2 Re=300
He ein, he i s case conside ed is o a Reynold Numbe Re = 300. The low is
simula ed o an o e all du a ion o ∗= 250 ime uni , whe e ∗= U/D, and only
he las 120 ime uni s a e conside ed o collec ing he low s a is ics.
In Figu e 6.2(a) he ime his o y o he ae odynamical coe icien s is shown: he
a e age alue o he d ag CD=Fx/(0.5ρU2
∞πD2/4) and li CL=Fy/(0.5ρU2
∞πD2/4)
a e espec i ely 0.667 and 0.066 wi h an oscilla ion ampli ude ∆CDand ∆CL
equal o 0.0032 and 0.17. These alues ag ee wi h o he nume ical in es iga ions,
as seen in Table 6.1, whe e he p esen esul s a e summa ized and compa ed
agains he li e a u e. F om Figu e 6.2(a) is also e iden ha he side coe icien
CS=Fy/(0.5ρU2
∞πD2/4) is equal o 0, indica ing ha he plane (x, y) is a plane o
symme y o he p esen low, as expec ed om p e ious expe imen al and nume ical
in es iga ions wi h Re = 300 [242, 102].
F om he equency analysis in Figu e 6.2(b), whe e he Powe Spec al Densi ies
o CDand CLa e epo ed, i can be ound he alue o he non-dimensional o ex
shedding equency S = 0.134, which is again in ag eemen wi h o he nume ical
in es iga ion, o example wi h he esul s in [242] whe e is p edic ed a a alue
S = 0.137. In e es ingly, as in [102], he d ag coe icien p esen s a peak in he
spec um a a alue equal o 2S , which is no p esen in he li coe icien spec um.
126
(a) (b)
Figu e 6.2 (a) Time his o ies and (b) no malized powe spec al densi ies o he
ae odynamic coe icien s o he low pas a sphe e a Re = 300.
To u he in es iga e his aspec , Figu e 6.3 displays he ime his o y o he
s eam-wise eloci y and i s powe spec al densi y a a poin dis an 5.75D om he
sphe e in he downs eam di ec ion. The only peaks p esen a e a a alue equal o
S = 0.134, along wi h i s h ee supe -ha monics a 2S , 3S , and 4S , as in [242].
Figu e 6.4 p esen s he a e age s eamwise eloci y Ua g along he wake
cen e line, wi h a compa a i e analysis agains he nume ical ou comes om [242].
A good co ela ion is obse ed up o a poin x/D = 7, beyond which he displayed
esul s indica e lowe alues o Ua g. The oo -mean-squa e alue o he s eamwise
eloci y on he same cen e line is also highligh ed in he igu e, wi h he cu en
simula ion p edic ing a diminished peak alue a x/D = 1.7. Fu he downs eam
he e is a easonable ag eemen be ween he wo se s o esul s.
The low isualiza ion displayed in Figu es 6.5-6.6, whe e he Q c i e ion, a
o ex iden i ica ion me hod, is shown a e e y qua e o he o ex shedding pe iod,
helps unde s and he de elopmen o he o ical s uc u es in he wake. In Figu es
6.6(a)-6.5(a), a o ex is abou o sepa a e om he o ical s uc u es su ounding
he sphe e. This cohe en s uc u e is con ec ed downs eam, as can be seen in
Figu es 6.6(b)-6.5(b)), whe e a e now clea ly isible he legs o he hai pin o ex.
127
(a) (b)
Figu e 6.3 (a) Time his o y and (b) no malized powe spec al densi y o he alue
o he s eamwise coe icien a a x/D = 5.75 on he wake cen e line o he low pas
a sphe e a Re = 300.
Figu e 6.4 Compa ison o he a e age and oo -mean-squa e alues o he s eamwise
eloci y along he wake cen e line wi h he nume ical esul s in [242]
128
While he shed-wake s uc u e mo es u he away om he sphe e [Figu es 6.6(c)-
6.5(c)], a new s uc u e, esul ing om he in e ac ion be ween he wake and he ou e
low, is o med a ound he legs o he hai pin o ex (6.6(d)-6.5(d)). The las igu e
also shows he de elopmen o a new cohe en s uc u e in he nea -wake egion. To
summa ize, he o ex-shedding mechanism is cha ac e ized by wo hai pin s uc u es
o ien ed in opposi e di ec ions; one esul s om he wake-shed, he o he om he
in e ac ion be ween he wake and he ou e low.
6.1.3 Re=500
To u he alida e he p esen mul i- esolu ion algo i hm, he Reynolds numbe is
aised up o 500, whe e he low pas a sphe e is cha ac e ized by he loss o he
plana symme y and a mo e complex o ex shedding. The nume ical simula ion is
compu ed o a o al o ∗= 300 ime uni , and he las ∗= 180 a e used o collec
he s a is ics.
In Figu e 6.7(a), he ime his o y o he o ce coe icien is epo ed. The a e age
alue o CDo he p esen simula ion is equal o 0.566, which is in pe ec ag eemen
wi h he esul s in [160]. Wi h espec o he case wi h Re = 300, he loss o plana
symme y can be disce ned om he e olu ion o he side coe icien .
The spec al analysis o he ae odynamical coe icien s is epo ed in Figu e
6.7(b), e eals he mo e complex na u e wi h espec o he case wi h Re = 300,
whe e only one dominan peak, co esponding wi h he o ex shedding equency,
was p esen in he spec um o CD. Ins ead, in his case, he powe spec al densi y
o he d ag coe icien is cha ac e ized by di e en peaks, wi h he dominan one a
S 2= 0.44, which is close o he alue ound in [126, 51]. As o he p e ious case, o
u he analyze he low’s spec al cha ac e is ic, in Figu e 6.8, he no malized powe
spec al densi ies o he s eam-wise eloci y a e shown, and he p essu e a h ee
di e en loca ions in he nea wake egion, chosen o ma ch he same measu emen
129
(a)
(b)
(c)
(d)
Figu e 6.5 Flow isualiza ion by he Q c i e ion, colo ed by he eloci y magni ude
U, a e e y qua e o pe iod om a iew no mal o he (x, z) plane o he low pas
a sphe e a Re = 300.
130
(a)
(b)
(c)
(d)
Figu e 6.6 Flow isualiza ion by he Q c i e ion, colo ed by he eloci y magni ude
U, a e e y qua e o pe iod om a iew no mal o he (x, y) plane o he low pas
a sphe e a Re = 300.
131
(a) (b)
Figu e 6.7 (a) Time his o ies and (b) no malized powe spec al densi ies o he
ae odynamic coe icien s o he low pas a sphe e a Re = 500.
Table 6.2 Nume ical Resul s o he D ag Coe icien CD, S ouhal Numbe S 1and
S 2Values o he Flow Pas a Sphe e a Re=500
Case CDS 1S 2
P esen 0.566 0.0167 0.044
Mi al and Najja [160] 0.56 0.150 0.05
Lee [126] 0.54 0.164 0.045
Tomboulides and O szag [242] - 0.167 0.045
C i ellini e al. [51] 0.558 0.157 0.044
poin s in [242, 51]. The i s poin on he wake cen e line a (x, y, z) = (2.5D, 0,0)
shows he p esence o wo di e en equencies, one a S 2= 0.44, he second a S 1=
0.166, while a loca ions 2 and 3, he dominan equency is he S 1= 0.166. The
alues o he low and high equencies ound a e in good ag eemen wi h he esul s
in [242], as can be seen om Table 6.2, whe e he alue o he d ag coe icien and
he S ouhal numbe s ound in his s udy a e compa ed agains o he expe imen al
and nume ical in es iga ions.
132
(a) S eamwise eloci y (b) P essu e
(c) S eamwise eloci y (d) P essu e
(e) S eamwise eloci y ( ) P essu e
Figu e 6.8 No malized powe spec al densi y o he s eamwise eloci y and
p essu e a (a)-(b) (x, y, z) = (2.5D, 0,0), (c)-(d) (x, y, z) = (2D, 0.3,0), (c)-(d)
(x, y, z) = (2D, 0.0,0.3), o he low pas a sphe e a Re = 500.
133
In Figu e 6.9 he iso-su aces o he Q c i e ion a e e y pe iod T1= 1/S 1a e
displayed: i is clea ha he S 1is he equency o he o ex shedding mechanism,
which esembles he same shedding p ocess seen a Re = 300, al hough in his case,
he o ex s uc u e sequence is mo e i egula . Mo eo e , i can be seen ha he
o ex o ien a ion changes om cycle o cycle.
6.2 3-D Dam-B eaking Flow
He ein, he abili y o he p oposed mul i- esolu ion algo i hm o deal wi h h ee-
dimensional ee-su ace lows is alida ed by simula ing a dam-b eaking low ha
impac s an obs acle. This es case is e y popula in he SPH communi y and has
been chosen as SPHERIC benchma k case #2. The con igu a ion o he nume ical
expe imen [113] is epo ed in Figu e 6.10.
A baseline spa ial esolu ion ∆xmax = 0.005 is chosen, and a e inemen egion
has been de ined a ound he obs acle, whe ein a ine esolu ion o ∆xmin = 0.0025
is applied, as shown in Figu e 6.11. A uni o m alue ac oss he di e en sub-domain
o he smoo hing leng h o pa icle dis ance is se equal o h/∆x= 1.5.
The e e ence densi y has a alue equal o ρ0= 1000 kg/m3, while he
g a i a ional accelle a ion g=−9.81 m/s2. The sound speed is chosen equal o
c0= 10√gH, whe e H= 0.55mis he heigh o he wa e column. To ensu e
he nume ical s abili y, an a i icial iscosi y model wi h α= 0.01 is used, which
co esponds o a physical kinema ic iscosi y ν= 8.6×1−−5, along wi h he DDT
e m in Equa ion 3.66. No-slip wall bounda y condi ions a e en o ced on he walls.
Figu e 6.12 shows he compa ison be ween he nume ical simula ion and he
expe imen al esul s o [113] a he p essu e p obes placed on he on side o he
obs acles. A good ag eemen is ound a gauges P1 and P2, whe e he nume ical
solu ion is able o p edic qui e well he e olu ion o he p essu e a he i s impac .
Some disc epancies a e ins ead ound a P3 and P4, whe e he SPH simula ion unde -
134