Acta hydrotechnica 32/57 (2019), Ljubljana
ISSN 1581-0267
Open Access Journal
Odprtodostopna revija
UDKIUDC: 626.882:532.51(078.9)
Izvirni znanstveni članek - Original scientific paper
Prejeto/Received: 10.03.2020 Sprej etolAccepted: 16.04.2020
Določitev koeficienta upora potopljenega telesa z uporabo metode SPH Evaluation of the drag coefficient of a fully submerged body using SPH
Gorazd Novak1'*, José M. Domínguez2, Angelo Tafuni3, Matjaž Četina1, Dušan Žagar1
1 Fakulteta za gradbeništvo in geodezijo, Univerza v Ljubljani, Jamova 2, Ljubljana, Slovenija 2 EPHYSLAB Environmental Physics Laboratory, University of Vigo, Ourense, Spain 3 School of Applied Engineering and Technology, New Jersey Institute of Technology, Newark, USA
Odprtokodno programsko opremo DualSPHysics smo uporabili za trodimenzijsko modeliranje vodnega toka s prosto gladino v ozkem vodoravnem odprtem kanalu pravokotnega prečnega prereza z namenom določitve sil, ki delujejo na različne mirujoče potopljene objekte. Simulacije so obravnavale obtekanje mirujoče krogle, kocke in kapsuli podobnega telesa, sestavljenega iz treh osnovnih geometrijskih teles. Obravnavali smo različne globine potopitve in različne hitrosti vodnega toka. Izračunane vzdolžne sile naraščajo s kvadratom hitrosti, prečne so enake nič, navpične pa konstantne. Z umerjenim in verificiranim modelom smo izračunali koeficient upora za potopljeno sestavljeno telo. Raziskava je pokazala na nekatere omejitve pri uporabi orodja DualSPHysics in bo pripomogla k nadaljnjemu razvoju kode in nadgradnji raziskav toka v ribji stezi z vertikalnimi prekati.
Ključne besede: CFD, FSI, 3-D model, metoda SPH, DualSPHysics, hidrodinamika, sile, potopljeno telo.
Abstract
The open-source software DualSPHysics was used so set up a three-dimensional model of a free-surface water flow in a narrow horizontal open channel with a rectangular cross section, and to determine the forces acting on various submerged fixed objects. Simulations included a sphere, a cube, and a composite capsule-like object made out of three basic bodies. Objects were located in the channel's central axis at various distances from the channel bed, and were exposed to various flow velocities. Calculating the hydrodynamic forces acting upon the objects resulted in longitudinal forces increasing with the square of the velocity, transverse forces equal to zero, and constant vertical forces. The calibrated and validated model was used to evaluate the drag coefficient of the submerged composite capsule-like object. This research reveals certain limitations of DualSPHysics and facilitates future improvements of the code and the upgrade of studies related to vertical slot fishways.
Keywords: CFD, FSI, 3-D model, SPH method, DualSPHysics, hydrodynamics, forces, submerged body
* Stik / Correspondence: gorazd.novak@fgg.uni-lj .si
© Novak G. et al.; Vsebina tega članka se sme uporabljati v skladu s pogoji licence Creative Commons Priznanje avtorstva -Nekomercialno - Deljenje pod enakimi pogoji 4.0.
© Novak G. et al.; This is an open-access article distributed under the terms of the Creative Commons Attribution - NonCommercial - ShareAlike 4.0 Licence.
https://doi.org/10.15292/acta.hydro.2019.08
Izvleček
107
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
1. Uvod
V prispevku so prikazani rezultati numeričnih simulacij za določitev sil na različna telesa, potopljena v toku. Ti rezultati so uporabni na različnih področjih, zlasti pri raziskavah interakcij med objekti in tekočino (angl. fluid-structure-interaction, FSI); v našem primeru omogočajo nadgradnjo raziskave toka v ribji stezi, ki je opisana v Novak et al. (2019).
Ribje steze imajo velik ekološki pomen, saj premoščajo prekinitve migracijskih poti, ki jih povzročajo različne pregrade na vodotokih. Ribje steze so bile in še vedno ostajajo deležne velike pozornosti, tako z vidika terenskih raziskav (Bombač et al., 2015), fizičnih hidravličnih modelov (Tan et al., 2018), globinsko povprečenih dvodimenzionalnih (2-D) numeričnih modelov (Bermudez et al. 2010, Bombač et al. 2017) in različnih trodimenzionalnih (3-D) numeričnih modelov (Duguay et al. 2017, Sanagiotto et al. 2019). Med slednje sodi tudi raziskava Novak et al. (2019), kjer je bila ribja steza z vertikalnimi prekati uspešno modelirana z uporabo metode hidrodinamike zglajenih delcev (metoda SPH, angl. smoothed-particle hydrodynamics). Funkcionalnost ribjih stez je največkrat obravnavana in modelno optimizirana z vidika hidravličnih parametrov, kot so hitrostna polja, tokovnice, turbulentna kinetična energija in vrtinčnost, le izjemoma (oz. samo posredno) pa z vidika hidrodinamičnih sil. Ena izmed takih raziskav je bila npr. Lauder in Drucker (2002), ki pa je bila izrazito eksperimentalna in ni vključevala niti numeričnega modela niti metode SPH.
Da bi zapolnili to vrzel s 3-D modelom, temelječim na metodi SPH, smo izvedli simulacije, ki so prikazane v nadaljevanju. Na osnovi teh rezultatov bo možno izvesti podrobnejšo analizo hidrodinamičnih sil, ki delujejo na ribo med njenim gibanjem v ribji stezi. Poleg tega bo model uporaben tudi za druge aplikacije na področju FSI, npr. za sile na objekte v rečnem ali morskem toku in interakcijo med drobirskimi tokovi ter okolico.
Novost te raziskave je uporaba metode SPH za račun hidrodinamičnih sil na različno oblikovana
telesa, ki so popolnoma potopljena v vodni tok. Metoda se je doslej sicer uporabljala tudi za določitev sil, vendar predvsem v povezavi s porušitvami pregrad in morskimi valovi. Tako npr. raziskava Aureli et al. (2015) obravnava silo porušitvenega vala, ki pljuskne ob objekt; Westphalen et al. (2009) so z različnimi modeli, med katerimi je bil tudi SPH, analizirali delovanje pretvornikov energije valov; Aghili et al. (2014) so obravnavali interakcijo med valovi in potopljeno vodoravno ploščo; Tafuni et al. (2016) pa so s to metodo računali tlake, ki jih na morskem dnu povzroča plovba hitrega čolna.
Namen prispevka je dvojen: 1) preveriti, ali je 3-D model, s katerim so bila uspešno simulirana hitrostna polja v omenjeni ribji stezi (Novak et al., 2019), primerno orodje tudi za določitev sil, ki na ribo delujejo v času njenega gibanja vzdolž ribje steze, in 2) numerično določiti še neznani koeficient upora pri obtekanju popolnoma potopljenega telesa, sestavljenega iz osnovnih geometrijskih teles.
2. Sile, ki delujejo na telesa v toku
»Določanje sile, s katero deluje tekočina na telo, ki ga ta tekočina obliva, ali pa, s katero mirujoča tekočina deluj e na telo, ki se giblje v njej, je osnoven inženirski problem, ki v praksi zelo pogosto nastopa (npr. delovanje sile vetra na stavbe in dimnike, upor avtomobilov v zraku, upor podvodnega dela ladij, upor snežnih plazov na zaščitne zgradbe itd.)« (cit. po Rajar, 1997, str. 216). V dani raziskavi gre za primer povsem potopljenega mirujočega telesa, ki ga obteka voda, ki teče v vzdolžni smeri x. Pri tem na telo v vzdolžni smeri deluje sila upora Fx, ki je v splošnem kombinacija upora oblike (vpliv normalnih napetosti) in običajno veliko manjšega upora trenja (vpliv strižnih napetosti). Poleg skupne sile upora na telo delujeta še njegova teža G in vzgon Fz, obe v navpični smeri z, kot za tok s hitrostnim poljem v* = (u, 0,0) prikazuje slika 1. Sila v prečni smeri Fy je enaka nič, saj gre za vpliv hidrostatičnega tlaka, ki je ob bokih telesa nasprotno enak. Teža telesa je dana z enačbo (1):
G = Pteio • Vteio • 9	(1)
108
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
kjer je ptelo gostota telesa [kg/m3], Vteo je celotni volumen telesa [m3], g pa je gravitacijski pospešek [kg m/s2]. Vzgon, ki deluje navpično navzgor, je posledica razlike hidrostatičnih tlakov in je enak teži izpodrinjene tekočine (enačba 2):
Fz Pvoda • Vpotopljeno • d
(2)
kjer je pvoda gostota vode [kg/m3] (tu privzamemo konstantno vrednost 1000 kg/m3), Vpotopljeno pa je volumen potopljenega dela telesa (kar je v primeru povsem potopljenega telesa kar volumen telesa) [m3].
u
r
I_^x
G
Slika 1: Sile pri obtekanju telesa.
Figure 1: Forces on a body submerged in a flow.
Privzeli smo, da je gostota telesa enaka gostoti vode (saj gre za aproksimacijo ribe), zato je teža telesa G po enačbi (1) enaka sili vzgona Fz po enačbi (2). Upor oblike in upor trenja nastopata hkrati in ju je navadno težko izvrednotiti ločeno, zato se celotni upor Fx najpogosteje izrazi z enačbo (3):
Fx = -• Cd
S • p • u2
(3)
kjer je S prečni presek telesa, p je gostota tekočine, u pa hitrost tekočine. Koeficient upora Cd je dobljen eksperimentalno in je funkcija Reynoldsovega števila Re, danega z enačbo (4):
u-d
Re = —
v
(4)
kjer je d karakteristična dolžina, v pa viskoznost tekočine. Vrednosti Cd so tabelirane za nekaj osnovnih teles (npr. Rajar, 1997), v primeru drugače oblikovanih potopljenih teles, kakršno nastopa tudi v tej raziskavi, pa moramo Cd določiti eksperimentalno ali z numeričnimi simulacijami. Primer take določitve je raziskava Marinho et al. (2009), kjer je bil hidrodinamični upor plavalca pod vodo določen s komercialnim 3-D numeričnim modelom.
Račun upora lahko upošteva tudi upor trenja in v ta namen uporabljamo bolj zapletene empirične enačbe, v katerih kot prispevna površina S nastopa omočeni del telesa. Primer takega pristopa je raziskava Granville (1976), ki obravnava obtekanje torpeda, a podobno, kot že omenjeno delo Marinho et al. (2009), kaže na to, da je upor oblike prevladujoč, enačba (3) pa torej dovolj natančna (tj. kot S se upošteva samo prečni prerez).
3. Metodologija
Za račun sil in določitev koeficienta Cd smo uporabili 3-D numerični model DualSPHysics, ki temelji na metodi SPH.
3.1 Glavne značilnosti metode SPH
Metoda SPH je ena izmed metod na področju računalniške dinamike tekočin (angl. computational fluid dynamics, CFD). Gre za brezmrežno Lagrangeevo metodo, ki tekočino aproksimira z masnimi delci, njihov medsebojni vpliv pa modelira s t. i. jedrno funkcijo W. Za vsak delec se vpliv sosednjih delcev upošteva z določeno utežjo znotraj domene, ki je definirana z značilno jedrno dolžino h, imenovano dolžina glajenja. Jedrna funkcija predpisuje, kako vpliv z oddaljenostjo od posameznega delca pada. Jedrna funkcija Wje lahko različnih oblik, najpogosteje se uporabljajo kubična (Monaghan in Lattanzio, 1985), Wendlandova (Wendland, 1995) in Gaussova.
Podobno kot drugi modeli dinamike tekočin tudi metoda SPH rešuje Navier-Stokesovi enačbi, tj. kontinuitetno enačbo (ohranitev mase) in dinamično enačbo (ohranitev gibalne količine), dopolnjuje pa ju enačba stanja, ki povezuje tlake in gostoto. Pri tem SPH obravnava vodo kot tekočino z majhno stisljivostjo. Podrobnosti so podane v Crespo et al. (2015).
Podobno kot v ostalih modelih je tudi v SPH ključnega pomena ločljivost modela, ki pa je tu izražena z začetno razdaljo med delci, tj. parametrom dp. SPH temelji na preračunavanju medsebojnih interakcij masnih delcev (v 3-D primerih velike ločljivosti število delcev pogosto presega milijon) v zelo kratkih računskih korakih (ki so potrebni, ker je metoda eksplicitna), zato so
109
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
simulacije zelo zahtevne z vidika računalniških virov. Ta pomanjkljivost se v zadnjem obdobju vse uspešneje rešuje z uporabo zmogljivih grafičnih kartic (GPU) oz. s paralelnim programiranjem (Domínguez et al., 2013).
3.2 DualSPHysics
Glavni poudarki teoretičnega ozadja programske opreme DualSPHysics so podrobneje predstavljeni v Crespo et al. (2015). DualSPHysics je odprtokodna programska oprema, ki je namenjena modeliranju inženirskih problemov, v katerih nastopajo nelinearni in večfazni tokovi ter tokovi z zelo razgibano prosto gladino. DualSPHysics je bil doslej uporabljen predvsem na področju porušitvenih valov, morskega valovanja in na področju FSI (npr. pljuskanje v cisternah). Z nedavno dodano možnostjo modeliranja odprtih robnih pogojev se vse bolj uveljavlja tudi pri reševanju problemov s turbulentnim tokom, kakršen nastopa denimo v ribji stezi (Novak et al., 2019). Podobno kot v omenjeni raziskavi so tudi v tem delu robni pogoji simulirani z uporabo koncepta t. i. navideznih vozlišč (angl. ghost nodes), ki je podrobneje opisan v Tafuni et al. (2018).
DualSPHysics omogoča modeliranje viskoznosti tekočine na enega izmed dveh načinov: 1) s konceptom umetne viskoznosti a, kjer je tekočina modelirana kot neviskozna, numerična viskoznost pa je namenjena stabilizaciji računske sheme, ali 2) s kombinacijo fizikalne viskoznosti ter modela turbulence SPS (angl. sub-particle scale). Slednji se je v tej raziskavi, podobno kot tudi že v Novak et al. (2019), izkazal za manj ustreznega, zato je bila uporabljena umetna viskoznost a.
Pri obeh načinih obravnavanja viskoznosti ima uporabnik možnost izbrati vrednost koeficienta, imenovanega visco-bound, s katerim se pomnoži viskoznost na stiku med tekočinskimi in trdnimi delci. V vseh simulacijah pričujoče raziskave je bil ta faktor nastavljen na vrednost 0, s čimer so bila potopljena telesa in stene kanala efektivno modelirani kot gladki objekti.
Paket DualSPHysics poleg istoimenskega solverja vključuje različna orodja za pred- in postprocesiranje. V raziskavi smo uporabili
naslednja: GenCase, PartVTK in ComputeForces. GenCase iz uporabniško določene geometrije modela generira domeno s tekočinskimi, trdnimi in/ali praznimi delci; orodje PartVTK je namenjeno vizualizaciji glavnih rezultatov; orodje ComputeForces pa na podlagi pospeškov masnih delcev izračuna vse tri komponente sil, s katerimi tekočina deluje na izbrano trdno telo.
GenCase za generiranje delcev uporablja 3-D Kartezijev sistem, tj. delci so generirani v vozliščih 3-D mreže. S pomočjo že vgrajenih ukazov je možno generirati različne oblike (npr. točke, črte, trikotnike, večkotnike, piramide, prizme, krogle, valje in elipsoide) ter telesom predpisati, ali so tekočinska, trdna ali prazna, kakšna je njihova masa, kako se gibljejo ipd. Za generiranje kompleksnih oblik se v model lahko uvozi tudi t. i. stereo-litografske (STL) objekte, ki so narejeni z namenskimi orodji CAD računalniško podprtega oblikovanja.
Program ComputeForces izračuna sile tekočine na telo z uporabo dinamične enačbe, v kateri nastopa Wendlandova jedrna funkcija W. Za izbrane tekočinske in trdne delce na stiku med tekočino in trdnim telesom se numerični pospešek izračuna z upoštevanjem interakcij med delci (enačba 5):
% = - lb mb ( J + | + na6) Vawab + g (5)
kjer indeksa a in b označujeta sosednje delce, P je tlak tekočinskih delcev, nab pa je viskoznostni člen (DualSPHysics, 2018). Sila se nato izračuna kot vsota produktov mase in pospeškov (enačba 6):
F = mJ1d^
^ dt
(6)
Za časovno integracijo lahko v DualSPHysics uporabimo računsko enostavnej šo Verletovo shemo in numerično bolj stabilno, a računsko zahtevnejšo simplektično (angl. symplectic) metodo. Pri eksplicitnih shemah integracije je časovni korak odvisen od Courant-Friedrichs-Lewyjevega (CFL) pogoja. V tem modelu je uporabljena simplektična metoda z variabilnim časovnim korakom, pri čemer je začetni časovni korak At,zač = 0,0001 s, minimalni pa At,min = 0,000001 s.
110
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
4. Simulacije
Simulacije so potekale v več fazah. V preliminarni smo najprej za primer mirujoče tekočine in navpične trdne stene, omočene z ene strani, potrdili, da je izračunana sila na ploskev enaka rezultanti hidrostatičnega tlaka. Poleg tega smo za obtekanje kocke potrdili, da pri dani konstantni hitrosti toka globina potopitve nima vpliva na izračunane sile, če je telo dovolj stran od ostenja. V tem modelu je razmerje med prečno dimenzijo telesa (premer krogle ali stranica kocke) in dimenzijo odprtega kanala (širina je enaka globini vode) enako 1 : 10.
Sledila je faza umerjanja, v kateri smo s simulacijami obtekanja krogle najprej določili vpliv izbire vrednosti dp in a. Potrdili smo, da so rezultati simulacij zelo odvisni tako od izbire dp (tj. od ločljivosti modela) kakor tudi od izbire a.
Sledila je faza verifikacije, v kateri smo umerjene vrednosti parametrov uporabili pri modeliranju obtekanja kocke.
Nazadnje smo določili še neznani koeficient upora Cd za sestavljeno telo. Pregled faz je prikazan v preglednici 1.
Preglednica 1: Pregled faz simulacij.
Table 1: Overview of simulations.
faza	opis	telo	z [m]	u [m/s]
1	preliminarna	stena, kocka	0,2-0,8	0-2, korak 0,5
2	umerjanje	krogla	0,5	0,25-2, korak 0,25
3	verifikacija	kocka	0,5	0,25-2, korak 0,25
4	določitev Cd	sestavljeno	0,5	0,25-2, korak 0,25
Geometrija modeliranega kanala je predstavljala ravni vodoravni kanal s pravokotnim prečnim prerezom, dolžine L = 2 m, širine B = 1 m in višine H = 1,1 m. V kanalu je tekel vodni tok z globino h = 1 m. Tok smo modelirali z odprtimi robnimi pogoji, tako da smo na dotočnem in iztočnem delu predpisali konstantno globino h in enakomerne hitrosti u.
Telo je v vsakem posameznem primeru mirovalo na lokaciji x = L/2, y = B/2, (tj. v osi kanala na sredini dolžine), izpostavljeno pa je bilo različnim hitrostim toka. Raziskani so bili naslednji parametri:
1)	Ločljivost modela (začetna razdalja med delci; manjša razdalja pomeni večjo ločljivost): dp = 0,6; 0,8; 1; 1,5; 2 cm.
2)	Oblika potopljenega telesa:
-	krogla s polmerom rkrogla = 5 cm,
-	kocka s stranico a = 10 cm,
-	kapsuli podobno telo, sestavljeno iz dveh krogel in valja, prikazano na sliki 2: r = 5 cm, l = 20 cm. Telo je bilo generirano in obravnavano kot enotno.
Koeficient upora Cd za kroglo oz. kocko znaša 0,47 oz. 1,05 (Rajar, 1997), za kapsuli podobno telo pa ni tabelirano. Obravnavani krogla in kocka sta imeli primerljivo velik prečni presek, kar je olajšalo primerjavo vpliva posameznih relevantnih parametrov.
r
Slika 2: Kapsuli podobno sestavljeno telo (vzdolžni prerez). Ni v merilu.
Figure 2: Composite object (longitudinal section). Not to scale.
3)	Pozicija telesa (v osi, na različnih oddaljenostih od dna) z = 0,2; 0,4; 0,5; 0,6; 0,8 m. Globine, različne od 0,5 m, so bile obravnavane le v preliminarni fazi.
4)	Hitrost toka u = 0,001; 0,25; 0,5; 0,75; 1; 1,25; 1,5; 1,75; 2 m/s. Hitrost u = 0,001 m/s je bila modelirana le v preliminarni fazi in je predstavljala hidrostatični primer, saj je bil pretok tedaj zanemarljiv. V vseh primerih je bilo predpisano, da se voda že od samega začetka in ves čas trajanja simulacije giblje s hitrostjo u, tako na vtoku, iztoku, kakor tudi vzdolž kanala (tj. stalni tok). Časovno spreminjanje sil za primer obtekanja krogle pri u = 1 m/s je prikazano na sliki 3.
111
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
0	1	2	3	4	5
t[s]
Slika 3: Časovno spreminjanje sil za primer krogle pri u = 2 m/s. Figure 3: Forces acting on the sphere at u = 2 m/s.
Iz slike 3 je razvidno, da sile zanihajo le na samem začetku simulacije (približno do časa 0,2 s), nato pa ostajajo konstantne, zaradi česar je čas simulacij lahko kratek. Omejili smo se na 2 s, kot rezultat pa smo upoštevali povprečje zadnje sekunde simulacije. Ker smo izbrali, da model izračuna izhodne podatke vsakih 5 stotink sekunde, je bil posamezni rezultat sile povprečje 20 vrednosti.
Za simulacijo z ločljivostjo dp = 1 cm, ki zajema približno N = 2,5 milijona začetnih delcev, sta bili na delovni postaji z eno grafično kartico NVidia GTX 1080 potrebni približno 2 uri računskega časa.
5. Rezultati in razprava
5.1. Preliminarne simulacije
S preliminarnimi simulacijami obtekanja kocke smo potrdili naslednje: 1) v dovolj velikem kanalu ostajajo sile Fx, Fy in Fz konstantne neodvisno od globine potopitve. Ker oddaljenost telesa od dna ne vpliva na velikost sil, so bile vse nadaljnje simulacije izvedene za z = 0,5 m. 2) Sila Fx narašča s kvadratom hitrosti u, sila Fy ostaja enaka nič, sila Fz pa ostaja enaka teži izpodrinjene tekočine. To je skladno s teorijo, vendar pa so izračunane sile po pričakovanjih zelo odvisne od ločljivosti modela, ki je dana z začetno razdaljo med delci dp.
5.2 Vpliv izbire vrednosti dp in a
V DualSPHysics sta med glavnimi parametri, s katerimi je treba umeriti model, vrednosti dp in a.
Za določitev vpliva izbire dp in a smo izvedli simulacije obtekanja krogle. Sila upora, ki deluje na kroglo pri u =1 m/s in u = 2 m/s - oboje pri priporočeni vrednosti a = 0,01 -, je za različne dp prikazana na sliki 4. Modra grafa prikazujeta računsko (tj. teoretično) vrednost sile in modelni rezultat pri u = 2 m/s, rdeča grafa pa pri u = 1 m/s.
Iz slike 4 je razvidno, da je pri veliki hitrosti (u = 2 m/s) vpliv dp zelo izrazit, medtem ko je pri srednji hitrosti (u = 1 m/s) vpliv dp manjši. Privzeli smo, da so z vidika umerjanja za večino primerov merodajnejši rezultati za srednjo hitrost, sile pri velikih hitrostih pa smo korigirali, kot je opisano v nadaljevanju.
Iz grafov za u = 1 m/s je razvidno, da je modelna sila za vse dp < 0,01 m praktično enaka in da se dobro ujema z računsko, kar pomeni, da večanje ločljivosti modela (tj. zmanjševanje dp), ki bistveno poveča računski čas, od določene vrednosti dp dalje ni več smiselno.
Kot optimalni dp smo izbrali dp = 0,01 m in to vrednost uporabili v vseh nadaljnjih simulacijah. Za modeliranje sile upora v večini tu obravnavanih primerov zadostuje taka ločljivost modela, da je prečna dimenzija modeliranih teles enaka 10 dp.
S tako izbranim dp smo za obe omenjeni hitrosti analizirali še vpliv izbire vrednosti a, kot je prikazano na sliki 5.
112
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
	Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana	
Slika 4: Sila upora za kroglo pri različnih dp. Figure 4: Drag force for the sphere at various dp values.
Slika 5: Sila upora za kroglo pri različnih a. Figure 5: Drag force for the sphere at various a values.
Iz slike 5 je razvidno, da je pri veliki hitrosti najboljše ujemanje med modelom in teorijo doseženo pri a = 0,03. Pri srednji hitrosti je zelo dobro ujemanje sil doseženo v območju a = 0,005 do a = 0,01. V tem razponu daje boljše ujemanje pri u = 2 m/s vrednost a = 0,01 in to vrednost smo uporabili v vseh nadaljnjih simulacijah. To je skladno s priporočilom, da se za večino aplikacij v odprtih kanalih lahko uporabi a = 0,01 (Altomare et
al., 2015). S tako izbranimi vrednostmi dp in a smo za obtekanje krogle dobili rezultate, ki so prikazani na sliki 6.
Iz slike 6 je razvidno zelo dobro ujemanje sil Fx za u <1 m/s, zelo dobro ujemanje sil Fy za vse obravnavane u in sprejemljivo ujemanje sil Fz, medtem ko so izračunane sile Fx za u > 1,25 m/s z naraščanjem u čedalje bolj prevelike. To verjetno
113
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
lahko pripišemo bolj turbulentnim razmeram (tj. bolj kompleksni tokovni sliki v neposredni okolici telesa), ki jih brez ustreznega modela turbulence ne moremo modelirati dovolj natančno.
Povečanje turbulentnosti je razvidno iz slike 7, ki prikazuje hitrosti u v vzdolžni osi kanala. Ker je bil
cilj raziskave določitev koeficienta upora, smo se omejili na sile Fx. Če določimo enačbo trendne črte, ki jo določajo rezultati Fx,m za u < 1 m/s (enačba je polinom drugega reda) in jo uporabimo za ekstrapolacijo sil pri u > 1 m/s, dobimo zadovoljivo ujemanje med Fx,m in Fx,r tudi v območju večjih hitrosti (slika 8).
Slika 6: Sile na kroglo pri različnih u. Figure 6: Forces on the sphere at various u values.
Slika 7: Hitrosti pri obtekanju krogle (levo) in kocke (desno), za u = 1 m/s (zgoraj) in u = 2 m/s (spodaj). Figure 7: Velocities for the sphere (left) and the cube (right) at u = 1 m/s (above) and u = 2 m/s (below).
114
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
10 9 8 7 6 5 4 3 2 1 0
						
						m
						
					_i	• /
						
Y =	2,36x2 -1	Q,946x +	0,5825			
	R2 =	0,9918				
						
						
						
0,25
0,5
0,75
-Fx,r
- Fx,m
1	1,25
u [m/s]
Fx,m po enačbi
1,5
1,75
Polinomska (Fx,m)
Slika 8: Sila upora za kroglo pri različnih u. Figure 8: Drag force for the sphere at various u values.
Slika 9: Sile na kocko pri različnih u. Figure 9: Forces on the cube for various u values.
Pristop z ekstrapolacijo smo preverili z modeliranjem obtekanja kocke, kot je opisano v nadaljevanju.
5.3 Verifikacija
Ustreznost opisanega postopka umerjanja smo preverili z modeliranjem obtekanja kocke. Rezultati so prikazani na sliki 9, kjer je razvidno, da so trendi
sil zelo podobni kot pri obtekanju krogle: sile Fx se lepo ujemajo do u = 1,5 m/s, sile Fy se skladajo z računskimi, sile Fz pa so za vse u približno enako prevelike.
Tako kot pri umerjanju so nas tudi tu zanimale predvsem sile Fx. Če jih ekstrapoliramo tako kot pri umerjanju, dobimo opazno boljše ujemanje tudi pri večjih hitrostih, kot kaže slika 10.
115
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
30
25
20
15
10
						
						
						
y = 6	,12x2-1,	334x + 0,	6925	1		
	R2 = 0	,9998	1	l/		
						
0,25 0,5 0,75	1	1,25 1,5 1,75	2
u [m/s]
>—Fx,r --■--Fx,m ■ Fx,m po enačbi	Polinomska (Fx,m)
Slika 10: Sila upora za kocko pri različnih u. Figure 10: Drag force for the cube at various u values.
10 9 8 7 6 . 5 4 3 2 1 0
						
						
						
						
						
						
y = l,88x - 0,33x + 0,842 R2 = 0,9888			5			
		S				
	•j." ."»'T					
___						
0,25 0,50 0,75 1,00 1,25 1,50 1,75 2,00 U [m/s]
--•--Fx,m —•—Fx,m po enačbi .........Polinomska (Fx,m)
Slika 11: Sile za določitev koeficient upora za sestavljeno telo. Figure 11: Forces for the evaluation of the drag coefficient for the composite body.
5.4 Koeficient upora za sestavljeno telo
Za obravnavano sestavljeno telo ne poznamo vrednosti Cd, saj telo ni med tipičnimi tabeliranimi primeri. Koeficient Cd smo določiti po naslednjem postopku: 1) z umerjenim in verificiranim modelom simuliramo obtekanje kapsuli podobnega telesa, 2) dobljene sile Fx za u < 1 m/s uporabimo za določitev
enačbe, s katero izračunamo sile Fx za u > 1 m/s, 3) iz tako dobljenih modelnih sil (slika 11) izračunamo Cd z uporabo enačbe (3).
Koeficient Cd se običajno izrazi v odvisnosti od Reynoldsovega števila Re. Ob privzeti kinematični viskoznosti vode u = 10-6 m2/s smo dobili rezultat, prikazan na sliki 12.
116
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
1,4						
						
1,2	_•_					
						
1,0						
						
0,8						
			M			
^ 0,6			_1_			
				• m		
0,4					» w	
						
0,2						
						
0,0						
						
0,00		0,50 1,00 1,50 2,00				
			Re [106]			
Slika 12: Koeficient upora za sestavljeno telo. Figure 12: Drag coefficient for the composite body.
Iz slike 11 je razvidno, da je koeficient pri nižjih hitrostih nerealno visok, z večanjem Re pa se ustali približno na Cd = 0,50. To je skladno s pričakovanji, saj je tako nekoliko večji kot pri krogli enakega premera (Cd = 0,47) in manjši kot pri horizontalnem valju z enakim razmerjem premera in dolžine (Cd = 0,85). Odstopanje od tabeliranih teoretičnih vrednosti predstavlja previsoka vrednost Re, pri kateri se Cd ustali, saj je v literaturi ta vrednost že pri Re = 10000.
6. Zaključki
3-D model DualSPHysics, ki temelji na metodi SPH, smo uspešno uporabili za simulacije obtekanja različnih potopljenih teles z namenom računa hidrodinamičnih sil in koeficienta upora za sestavljeno telo z netipično geometrijo.
Med umerjanjem modela smo pokazali, da za račun sil zadošča taka ločljivost modela (izražena z začetno razdaljo med delci dp), da je prečna dimenzija objekta, ki je odločilna za račun sile upora, enaka 10 dp.
Nadalje smo pokazali, da je za večino primerov najbolje uporabiti vrednost umetne viskoznosti a = 0,01, kar je skladno s priporočili v literaturi.
Za obravnavano obtekanje krogle in kocke je bilo ujemanje med modelnimi in teoretičnimi silami upora zelo dobro za hitrosti u < 1 m/s, pri večjih hitrostih pa so bile modelne sile izrazito prevelike. To odstopanje verjetno lahko pripišemo kompleksnejši tokovni sliki na območju stika med telesom in tekočino. Gre za problematiko, ki terja ustrezen model turbulence, kakršnega v DualSPHysics zaenkrat še ni, zato bomo to vprašanje obravnavali v nadaljnjem delu.
Kot rezultate za u > 1 m/s smo upoštevali vrednosti, ki smo jih ekstrapolirali na podlagi rezultatov za u < 1 m/s.
Iz tako dobljenih sil upora smo izračunali koeficient upora Cd za sestavljeno telo netipične oblike. Dobljene vrednosti Cd so skladne s pričakovanji, saj kažejo, da je upor takega telesa večji kot pri krogli enakega premera in manjši kot pri primerljivo velikem horizontalnem valju.
Poudariti velja, da gre v našem primeru za raziskavo, ki se loteva zelo kompleksnega področja (določitev Cd za potopljena telesa praviloma terja kombinacijo zahtevnih eksperimentov in numeričnih modelov) in da metoda SPH doslej ni bila uporabljena za račun hidrodinamičnih sil na povsem potopljena telesa v toku.
117
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
Na podlagi pričujoče raziskave bomo nadaljevali simulacije toka v ribji stezi (Novak et al., 2019) v smeri določitve hidrodinamičnih sil, ki delujejo na ribe v času njihovega potovanja po stezi. Z uporabo programske opreme DualSPHysics bo možno izvesti podrobne hidrodinamične analize brez uporabe neracionalno zahtevnega fizičnega modeliranja.
Ta model se lahko uporabi za račun koeficienta upora za različna sestavljena telesa, pozicionirana nesimetrično na tok (zasuka in/ali nagib telesa glede na smer toka), in za telesa, ki se gibljejo proti toku.
Viri
Aghili, M., Ghadimi, P., Maghrebi, Y. F., Nowruzi, H. (2014). Simulating the interaction of solitary wave and submerged horizontal plate using SPH method. Int. J. Phys.	Res.	2(2),	16-26.
https://doi.org/10.14419/ijpr.v2i2.2451.
Altomare, C., Domínguez, J.M., Crespo, A.J.C., Suzuki, T., Caceres, I., Gómez-Gesteira, M. (2015). Hybridisation of the wave propagation model SWASH and the meshfree particle method SPH for real coastal applications. Coastal Engineering Journal, 57(4) https://doi.org/10.1142/S0578563415500242.
Aureli, F., Dazzi, S., Maranzoni, A., Mignosa, P., Vacondio, R. (2015). Experimental and numerical evaluation of the force due to the impact of a dam-break wave on a structure. Advances in Water Resources 76, 29-42. https://doi.org/10.1016/j.advwatres.2014.11.009.
Bermúdez, M., Puertas, J., Cea, L., Pena, L., Balairón, L. (2010). Influence of pool geometry on the biological efficiency of vertical slot fishways, Ecol. Eng. 36, 13551364. https://doi.org/10.1016/j.ecoleng.2010.06.013.
Bombač, M., Novak, G., Mlačnik, J., Četina, M. (2015). Extensive field measurements of flow in vertical slot fishway as data for validation of numerical simulations, Ecol.	Eng.	84,	476-484.
https://doi.org/10.1016/j.ecoleng.2015.09.030.
Bombač, M., Četina, M., Novak, G. (2017). Study on flow characteristics in vertical slot fishways regarding slot optimization, Ecol. Eng. 107, 126-136. https://doi.org/10.1016/j.ecoleng.2017.07.008.
Crespo, A.J.C., Domínguez, J. M., Rogers, B. D., Gómez-Gesteira, M., Longshaw, S., Canelas, R., Vacondio, R., Barreiro, A., García-Feal, O. (2015). DualSPHysics: Open-source parallel CFD solver based
on Smoothed Particle Hydrodynamics SPH, Comput. Phys.	Commun.	187,	204-216.
https://doi.org/10.1016/j.cpc.2014.10.004.
Domínguez, J. M., Crespo, A. J. C., Gómez-Gesteira, M. (2013). Optimization strategies for CPU and GPU implementations of a smoothed particle hydrodynamics method, Comput. Phys. Commun. 184, 617-627. https://doi.org/10.1016/i.cpc.2012.10.015.
DualSPHysics (2018). Users Guide for DualSPHysics code. DualSPHysics v4.2. May 2018. 146 p.
Duguay, L., Lacey, R., Gaucher, J. (2017). A case study of a pool and weir fishway modeled with OpenFOAM and FLOW-3D, Ecol. Eng. 103, 31-42. https://doi.org/10.1016/j.ecoleng.2017.01.042.
Granville, P.S. (1976). Elements of the drag of underwater bodies. David W. Taylor Naval Ship Research and Development Center, Bethesda, Md., USA
Lauder, G. V., Drucker, E. G. (2002). Forces, Fishes, and Fluids: Hydrodynamic Mechanisms of Aquatic Locomotion, News Physiol. Sci. 17, 235-240. https://doi.org/10.1152/nips.01398.2002.
Marinho, D. A., Reis, V. M., Alves, F. B., Villas-Boas, J. P., Machado, L., Silva, A. J., Rouboa, A. I. (2009). Hydrodynamic Drag During Gliding in Swimming, J. App.	Biomech.	25,	253-257.
https://doi.org/10.1123/íab.25.3.253.
Monaghan, J. J., Lattanzio, J. C. (1985). A refined method for astrophysical problems. Astron. Astrophys. 149, 135-143.
Novak, G., Tafuni, A., Domínguez, J. M., Četina, M., Žagar, D. (2019). A Numerical Study of Fluid Flow in a Vertical Slot Fishway with the Smoothed Particle Method, Water, 2019, 11, 1928. https://doi.org/10.3390/w11091928.
Rajar, R. (1997). Hidromehanika. UL FGG, Ljubljana, 235 str.
Sanagiotto, D. G., Rossi, J. B., Bravo, J. M. (2019). Applications of Computational Fluid Dynamics in The Design and Rehabilitation of Nonstandard Vertical Slot Fishways, Water, 2019, 11, 199; https://doi.org/10.3390/w11020199.
Tafuni, A., Sahin, I., Hyman, M. (2016). Numerical investigation of wave elevation and bottom pressure generated by a planing hull in finite-depth water. Applied Ocean	Research	58,	281-291.
https://doi.org/10.1016/j.apor.2016.04.002.
118
Novak G. et al.: Določitev koeficienta upora potopljenega telesa z uporabo metode SPH - Evaluation of the drag
coefficient of a fully submerged body using SPH
_Acta hydrotechnica 32/57 (2019), 107-119, Ljubljana_
Tafuni, A., Domínguez, J. M., Vacondio, R., Crespo, A. J. C. (2018). A versatile algorithm for the treatment of open boundary conditions in Smoothed particle hydrodynamics GPU models, Comput. Methods Appl. Mech.	Engrg.	342,	604-624.
https://doi.org/10.1016/j.cma.2018.08.004.
Tan, J., Tao, L., Gao, Z., Dai, H., Shi, X. (2018). Modeling Fish Movement Trajectories in Relation to Hydraulic Response Relationships in an Experimental Fishway, Water, 2018, 10, 1511. https://doi.org/10.3390/w10111511.
Wendland, H. (1995). Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Advances in Computational Mathematics	4,	389-396.
http://dx.doi.org/10.1007/BF02123482.
Westphalen, J., Greaves, D. M., Williams, C. K., Taylor, P. H., Causon, D. M., Mingham, C. G., Hu, Z. Z., Stansby, P. K., Rogers, B. D., Omidvar, P. (2009). Extreme Wave Loading on Offshore Wave Energy Devices using CFD: a Hierarchical Team Approach. Proceedings of the 8th European Wave and Tidal Energy Conference, Uppsala, Sweden, 2009.
119