Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo DOKTORSKI ŠTUDIJSKI PROGRAM III. STOPNJE GRAJENO OKOLJE Kandidatka: ELVIRA DŽEBO, univ. dipl. inž. vod. in kom. inž. SIMULACIJE VALOV ZARADI PORUŠITEV PREGRAD Z BREZMREŽNO NUMERIČNO METODO HIDRODINAMIKE ZGLAJENIH DELCEV Doktorska disertacija štev: 3/GO SIMULATIONS OF DAM-BREAK WAVES WITH MESHLESS NUMERICAL METHOD OF SMOOTHED PARTICLE HYDRODYNAMICS Doctoral thesis No.: 3/GO Soglasje k temi doktorske disertacije je dala Komisija za doktorski študij na 29. redni seji, 13. junija 2012. Za mentorja je bil imenovan prof. dr. Matjaž Četina. Ljubljana, 3. oktober 2013 Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo Komisijo za oceno ustreznosti teme doktorske disertacije v sestavi: prof. dr. Matjaž Četina, - doc. dr. Dušan Žagar, - prof. dr. Matjaž Mikoš, - izr. prof. dr. Anton Bergant, Litostroj Power in UL FS. - je imenoval Senat Fakultete za gradbeništvo in geodezijo na 29. redni seji, 28. marca 2012. Poročevalce za oceno doktorske disertacije v sestavi: doc. dr. Dušan Žagar, - prof. dr. Matjaž Mikoš, - izr. prof. dr. Anton Bergant, Litostroj Power in UL FS. - je imenoval Senat Fakultete za gradbeništvo in geodezijo na 2. redni seji, 26. junija 2013. Komisijo za zagovor doktorske disertacije v sestavi: - prof. dr. Matjaž Mikoš, dekan UL FGG, predsednik in član, - prof. dr. Matjaž Četina, mentor, - doc. dr. Dušan Žagar, - izr. prof. dr. Anton Bergant, Litostroj Power in UL FS, je imenoval Senat Fakultete za gradbeništvo in geodezijo na 3. redni seji, 25. septembra 2013. Univerza v Ljubljani Fakulteta za gradbeništvo in geodezijo IZJAVA O AVTORSTVU Podpisana ELVIRA DŽEBO, univ. dipl. inž. vod. in kom. inženirstva, izjavljam, da sem avtorica doktorske disertacije z naslovom: »SIMULACIJE VALOV ZARADI PORUŠITEV PREGRAD Z BREZMREŽNO NUMERIČNO METODO HIDRODINAMIKE ZGLAJENIH DELCEV«. Izjavljam, da je elektronska različica v vsem enaka tiskani različici. Izjavljam, da dovoljujem objavo elektronske različice v repozitoriju UL FGG. Ljubljana, 3. oktober 2013 ……………………………….. (podpis) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev VII Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. POPRAVKI Stran z napako Vrstica z napako Namesto Naj bo VIII Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev IX Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. BIBLIOGRAFSKO – DOKUMENTACIJSKA STRAN IN IZVLEČEK UDK: 519.61:532.511.011.5(043.3) Avtor: Elvira Džebo Mentor: prof. dr. Matjaž Četina Naslov: Simulacije valov zaradi porušitev pregrad z brezmrežno numerično metodo hidrodinamike zglajenih delcev Obseg in oprema: 121 str., 12 pregl., 107 sl., 68 en. Ključne besede: porušitev pregrade, tok s prosto gladino, hidrodinamika zglajenih delcev, metoda SPH Izvleček Poplavni valovi lahko povzročijo velike naravne katastrofe ter ogrozijo lastnino in življenja ljudi. Med najnevarnejšimi poplavnimi valovi so porušitveni valovi, ki nastanejo zaradi morebitnih porušitev pregrad, zato so vnaprejšnje simulacije tovrstnih valov zelo pomembne. Rezultati simulacij vplivajo na odločitev glede morebitne potrebe po povečanju poplavne varnosti. Simulacije takšnih primerov običajno izvajamo z uporabo matematičnih modelov. Način reševanja osnovnih enačb matematičnega modela pa je odvisen od izbrane numerične metode. Običajno so v uporabi mrežne Eulerjeve metode, ki pa imajo nekatere slabosti (npr. numerična difuzija). Zato se razvijajo nove brezmrežne Lagrangeove metode. Ena izmed najbolj priljubljenih brezmrežnih metod je metoda hidrodinamike zglajenih delcev (SPH). V okviru disertacije smo uporabili matematični model Tis Isat, ki osnovne enačbe rešuje po tej metodi. Model smo ustrezno dopolnili z novimi vrstami robnih pogojev in vmesnikom, ki nam omogoča povezovanje med modeli različne resolucije ter njegovo natančnost preverili na različnih testnih primerih, kjer so bili za primerjavo na voljo eksperimentalni podatki in/ali referenčni numerični rezultati iz literature. Uporabnost razvitega modela smo nato preverili še na praktičnem primeru porušitve dela nasipa predvidene akumulacije črpalne elektrarne Kolarjev vrh. Na tem primeru smo na osnovi primerjave z rezultati fizičnega modela preverili različne možnosti definiranja hrapavosti terena, primerjali rezultate različnih velikosti delcev ter izvedli parametrično analizo modela. Razvit in v okviru disertacije dopolnjen model Tis Isat se je pokazal kot uporabno orodje za simulacijo primerov toka s prosto gladino, kjer prihaja do izrazitih sprememb hitrosti in gladin v toku. Uspešno je bil uporabljen tudi za praktični primer simulacije porušitvenega vala na realni topografiji, kar je eden redkih tovrstnih primerov v svetu. Pri nadaljnjem delu bo za doseganje še večje natančnosti pri računu gladin potrebno uporabiti manjše delce, kar bo omogočil hitri razvoj vedno zmogljivejših računalnikov. X Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. BIBLIOGRAPHIC-DOCUMENTALISTIC INFORMATION UDC: 519.61:532.511.011.5(043.3) Author: Elvira Džebo Supervisor: prof. Matjaž Četina, PhD Title: Simulations of dam – break waves with meshless numerical method of smoothed particle hydrodynamics Notes: 121 p., 12 tab., 107 Figure, 68 eq. Key words: dam – break, free surface flow, smoothed particle hydrodynamics, SPH method Abstract Flood waves can cause large property damage and loss of human lives. Some of the most dangerous flood waves are destructive waves caused by dam failure. Accurate numerical simulations of such events are therefore of the most importance. Simulation results affect our decision about increasing flood protection. Such simulations are usually performed using mathematical models. The way of solving the basic equations of a mathematical model depends on the numerical method used. Mesh-based Eulerian methods are usually applied, but they have some disadvantages (e.g. numerical diffusion), so new meshfree Lagrangian methods have been developed. The most popular meshfree method for such simulations is Smoothed Particle Hydrodynamics (SPH). In this doctoral thesis, the Tis Isat model is used, which solves the basic equations by using the SPH method. This model was upgraded with a new boundary condition and a new interface for coupling between models of different resolutions. The new model was validated against results of laboratory experiments and/or available numerical model results from literature. Our model gave good results for all dam- break simulations. We subsequently used it on a real case study simulation of a dam break at the upper storage reservoir of the pumped storage HPP Kolarjev vrh in Slovenia. Different ways of defining terrain roughness were also verified. We show that the new model is applicable to simulate dam-break flows with rapidly changing velocities and free surface. Tis Isat has also successfully been used to simulate a real-world problem. Such simulations are relatively rare. In order to achieve greater accuracy, simulations must be performed with smaller particle size. This will eventually be practicable when more powerful CPU processors are available. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev XI Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. ZAHVALA Za pomoč pri nastajanju doktorske disertacije se zahvaljujem mentorju prof. dr. Matjažu Četini. Iskrena hvala domačim za vso pomoč in skrb za malo Adejo. Posebna zahvala pa gre dr. Gregorju Petkovšku in podjetju CGS plus d.o.o. za vodstvo in pomoč pri dopolnjevanju računalniškega programa Tis Isat. XII Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. KAZALO VSEBINE 1 UVOD 1 1.1 Splošen opis problema 1 1.2 Zasnova doktorske disertacije 4 2 OSNOVNE ENAČBE IN NUMERIČNE METODE REŠEVANJA 7 2.1 Osnovne enačbe 7 2.1.1 Predstavitev osnovnih enačb 8 2.1.2 Modeli turbulence 9 2.2 Konvencionalne numerične metode reševanja 10 2.3 SPH metoda reševanja 11 2.3.1 Predstavitev osnovnih enačb metode SPH 11 2.3.2 Modeli turbulence po metodi SPH 17 2.3.3 Upoštevanje modelov turbulence in trenja ob ostenju v modelu Tis Isat 21 2.3.4 Diskretizacija po času 22 2.3.5 Robni pogoji 23 2.3.6 Časovni koraki 25 2.3.7 Izračun globine po metodi SPH 25 2.3.8 Natančnost rezultatov 26 3. RAZVOJ IN DOPOLNITVE MODELA TIS ISAT 27 3.1 Dopolnitev robnih pogojev 28 3.2 Povezave med modeli različne resolucije 29 4. PRIMERI TESTIRANJA MODELA TIS ISAT NA OSNOVI REZULTATOV S FIZIČNIH MODELOV 30 4.1 Porušitev dvodimenzionalnega vodnega stolpca 30 4.2 Porušitev tridimenzionalnega vodnega stolpca 42 4.3 Porušitev vodnega stolpca v kanalu z razširitvijo ter v kanalu z zožitvijo in razširitvijo 45 4.4 Porušitev vodnega stolpca v kanalu s hipno razširitvijo 62 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev XIII Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 5. UPORABA MODELA TIS ISAT ZA PRAKTIČNI PRIMER PORUŠITVE NASIPA ZGORNJE AKUMULACIJE ČRPALNE ELEKTRARNE KOLARJEV VRH 72 5.1 Akumulacija 73 5.2 Dolina dolvodno od porušenega dela nasipa 78 5.2.1 Gladek teren 79 5.2.2 Hrapav teren 90 6. ZAKLJUČKI IN NAPOTKI ZA NADALJNJE DELO 105 6.1 Zaključki in ugotovitve 105 6.2 Napotki za nadaljnje delo 108 7 POVZETEK 110 8 SUMMARY 112 VIRI 114 XIV Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. KAZALO PREGLEDNIC Preglednica 1: Začetni parametri 31 Preglednica 2: Začetni parametri simulacij v laboratorijskem kanalu (Džebo in sod., 2011) 39 Preglednica 3: Začetni parametri 3D simulacij v kanalu 44 Preglednica 4: Začetne globine vode v kanalih a in b s suhim in mokrim dnom (Džebo in sod., 2013 a) 47 Preglednica 5: Parametri simulacij v kanalih a in b (Džebo in sod., 2013 a) 47 Preglednica 6: Začetni parametri pri primeru hipne razširitve kanala (Džebo in sod., 2012 b) 63 Preglednica 7: Časi napredovanja čela vala; PCFLOW2D – ORTHOCURVE. 78 Preglednica 8: Časi napredovanja čela vala; Tis Isat (gladek teren, konstanten parameter viskoznosti ob stenah) 79 Preglednica 9: Čas napredovanja čela vala; Tis Isat (gladko dno – spremenljiv parameter viskoznosti ob stenah) 86 Preglednica 10: Čas napredovanja čela vala; Tis Isat (hrapav teren, dvignjena mrežna vozlišča) 92 Preglednica 11: Čas napredovanja čela vala z modelom Tis Isat (hrapav teren - točkovni elementi) 98 Preglednica 12: Računski časi simulacij 104 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev XV Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. LIST OF TABLES Table 1: Initial parametres 31 Table 2: Initial parametres of simulations in the laboratory channel (Džebo et al., 2011) 39 Table 3: Initial parametres of 3D simulations in the channel 44 Table 4: Water depths in channels a and b with dry and wet bottom (Džebo et al., 2013 a) 47 Table 5: Parametres of simulations in channels a and b (Džebo et al., 2013 a) 47 Table 6: Initial parametres in the case of sudden expansion of the channel (Džebo et al., 2012 b) 63 Table 7: Propagation time of the surge front; PCFLOW2D – ORTHOCURVE 78 Table 8: Propagation time of the surge front; Tis Isat model (smooth terrain – constant coefficient for wall – particle interaction) 79 Table 9: Propagation time of the surge front; Tis Isat (smooth terrain – variable coefficient for wall – particle interaction) 86 Table 10: Propagation time of the surge front; Tis Isat (rough terrain – elevated mesh nodes) 92 Table 11: Propagation time of the surge front; Tis Isat (rough terrain - point elements on the bottom of the slopes) 98 Table 12: Computational times 104 XVI Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. KAZALO SLIK Slika 1: Eulerjev (levo) in Lagrangeov (desno) pristop 7 Slika 2: Jedrna funkcija 15 Slika 3: Izračun globine po metodi modela Tis Isat: a) tlorisni prikaz 2D SPH modela; b) podolžni profil 2D SPH modela; c) izračun vrednosti jedrne funkcije (Džebo in sod., 2013 a) 26 Slika 4: Povezovanje 2D/3D in 3D/2D modela (Džebo in sod., 2013 a) 28 Slika 5: Točkovni elementi na ravnem terenu. 28 Slika 6: Povezovanje med modeli različne resolucije (Džebo in sod., 2012 b). 29 Slika 7: Fizični model (Martin in Moyce, 1952) 31 Slika 8: Propagacija čela vala; Ho/ Wo = 1 (Petkovšek in sod., 2010 a) 33 Slika 9: Relativna višina vodnega stolpca; Ho/ Wo = 1 (Petkovšek in sod., 2010 a) 33 Slika 10: Propagacija čela vala; Ho/ Wo = 2 (Petkovšek in sod., 2010 a) 34 Slika 11: Relativna višina vodnega stolpca; Ho/ Wo = 2 (Petkovšek in sod., 2010 a) 34 Slika 12: Propagacija čela vala; Ho/ Wo = 4 (Petkovšek in sod., 2010 a) 35 Slika 13: Relativna višina vodnega stolpca; Ho/ Wo = 4 (Petkovšek in sod., 2010 a) 35 Slika 14: Primerjava napredovanja čela vala 0,08 s po porušitvi pregrade; Ho/ Wo = 2, število delcev je 50 x 100 (Petkovšek in sod., 2010 a) 36 Slika 15: Primerjava napredovanja čela vala 0,3 s po porušitvi pregrade; Ho/ Wo = 2, število delcev je 50 x 100 (Petkovšek in sod., 2010 a) 37 Slika 16: Primerjava napredovanja čela vala 0,3 s po porušitvi pregrade; Ho/ Wo = 2, število delcev je 50 x 100 (Petkovšek in sod., 2010 a) 37 Slika 17: Kanal brez zapore (zgoraj) in kanal z zaporo (spodaj) (Džebo in sod., 2011) 38 Slika 18: Napredovanje čela vala v kanalu brez pregrade; t = 1 s (Džebo in sod., 2011) 39 Slika 19: Napredovanje čela vala v kanalu brez zapore (Džebo in sod,. 2011) 40 Slika 20: Višina vodnega stolpca v kanalu brez zapore (Džebo in sod,. 2011) 40 Slika 21: Napredovanje čela vala v kanalu z zaporo (Džebo in sod,. 2011) 41 Slika 22: Višina vodnega stolpca v kanalu z zaporo (Džebo in sod,. 2011) 42 Slika 23: Laboratorijski eksperiment (levo) ter laboratorijski eksperiment z lociranimi vertikalnimi sondami (desno) (Kleefsman in sod., 2005) 43 Slika 24: Primerjava med merjenimi in simuliranimi globinami vode v točki H4 44 Slika 25: Primerjava med merjenimi in simuliranimi globinami vode v točki H2 45 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev XVII Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 26: Kanal z razširitvijo (kanal a) ter kanal z zožitvijo in razširitvijo (kanal b) (Rajar, 1972) 46 Slika 27: Pojava, ki nastaneta v tovrstnih kanalih: a) kanal z razširitvijo b) kanal z zožitvijo (Rajar, 1972) 46 Slika 28: Napredovanje čela vala – primer 1 (Džebo in sod., 2013 a) 49 Slika 29: Primerjava med meritvami in rezultati simulacij 15,5 s po porušitvi pregrade – primer 1 (Džebo in sod., 2013 a) 49 Slika 30: Primerjava med meritvami in rezultati simulacij 27 s po porušitvi pregrade – primer 1 (Džebo in sod., 2013 a) 50 Slika 31: Primerjava med meritvami in rezultati simulacij 45 s po porušitvi pregrade – primer 1 (Džebo in sod., 2013 a) 50 Slika 32: Aksonometrični pogled dela kanala z razširitvijo 45 s po porušitvi pregrade, simuliran s povezanim 2D/3D modelom Tis Isat z večjimi delci. Delci v roza barvi prikazujejo območja z višjo, modri delci pa območja z nižjo globino (Džebo in sod., 2013 a) 51 Slika 33: Aksonometrični pogled dela kanala z razširitvijo 45 s po porušitvi pregrade, simuliran s povezanim 2D/3D modelom Tis Isat z manjšimi delci. Delci v roza barvi prikazujejo območja z višjo, modri delci pa območja z nižjo globino (Džebo in sod., 2013 a). 52 Slika 34: Aksonometrični pogled dela kanala z razširitvijo 45 s po porušitvi pregrade, simuliran z 3D modelom Tis Isat z večjimi delci. Delci v roza barvi prikazujejo območja z višjo, modri delci pa območja z nižjo globino (Džebo in sod., 2013 a) 53 Slika 35: Aksonometrični pogled dela kanala z razširitvijo 45 s po porušitvi pregrade, simuliran z 3D modelom Tis Isat z manjšimi delci. Delci v roza barvi prikazujejo območja z višjo, modri delci pa območja z nižjo globino (Džebo in sod., 2013 a) 54 Slika 36: Prečni profil; Y = 32,22 m (Džebo in sod., 2013 a) 55 Slika 37: Prečni profil; Y = 33,11 m (Džebo in sod., 2013 a) 55 Slika 38: Prečni profil; Y = 34,00 m (Džebo in sod., 2013 a) 56 Slika 39: Napredovanje čela vala – primer 2 (Džebo in sod., 2013 a) 56 Slika 40: Primerjava med meritvami in rezultati simulacij 16,5 s po porušitvi pregrade – primer 2 (Džebo in sod., 2013 a) 57 XVIII Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 41: Primerjava med meritvami in rezultati simulacij 35 s po porušitvi pregrade – primer 2 (Džebo in sod., 2013 a) 57 Slika 42: Napredovanje čela vala – primer 3 (Džebo in sod., 2013 a) 58 Slika 43: Primerjava med meritvami in rezultati simulacij 20 s po porušitvi pregrade – primer 3 (Džebo in sod., 2013 a) 58 Slika 44: Primerjava med meritvami in rezultati simulacij 30 s po porušitvi pregrade – primer 3 (Džebo in sod., 2013 a) 59 Slika 45: Primerjava med meritvami in rezultati simulacij 58 s po porušitvi pregrade – primer 3 (Džebo in sod., 2013 a) 59 Slika 46: Napredovanje čela vala – primer 4 (Džebo in sod., 2013 a) 60 Slika 47: Primerjava med meritvami in rezultati simulacij 19,5 s po porušitvi pregrade – primer 4 (Džebo in sod., 2013 a) 60 Slika 48: Primerjava med meritvami in rezultati simulacij 26,5 s po porušitvi pregrade – primer 4 (Džebo in sod., 2013 a) 61 Slika 49: Primerjava med meritvami in rezultati simulacij 40 s po porušitvi pregrade – primer 4 (Džebo in sod., 2013 a) 61 Slika 50: Kanal z razširitvijo (Popovska, 1988) 62 Slika 51: Simulirana in merjena globina v točki zožitve in razširitve kanala 63 Slika 52: Prečni profil (od leve stene do sredine kanala) na razdalji 0,6 m od razširitve; t = 5 s. 64 Slika 53: Prečni profil (od leve stene do sredine kanala) na razdalji 1 m od razširitve; t = 5 s. 64 Slika 54: Prečni profil (od leve stene do sredine kanala) na razdalji 2 m od razširitve; t = 5 s. 65 Slika 55: Prečni profil (od leve stene do sredine kanala) na razdalji 0,6 m od razširitve; t = 7 s. 65 Slika 56: Prečni profil (od leve stene do sredine kanala) na razdalji 1 m od hipne razširitve; t = 7 s 66 Slika 57: Prečni profil (od leve stene do sredine kanala) na razdalji 2 m od hipne razširitve; t = 7 s. 66 Slika 58: Prečni profil (od leve stene do sredine kanala) na razdalji 3 m od hipne razširitve; t = 7 s 67 Slika 59: Prečni profil (od leve stene do sredine kanala) na razdalji 4 m od hipne razširitve; t = 7 s 67 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev XIX Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 60: Prečni profil (od leve stene do sredine kanala) na razdalji 0,6 m od hipne razširitve; t = 10 s 68 Slika 61: Prečni profil (od leve stene do sredine kanala) na razdalji 1 m od hipne razširitve; t = 10 s 68 Slika 62: Prečni profil (od leve stene do sredine kanala) na razdalji 2 m od hipne razširitve; t = 10 s. 69 Slika 63: Prečni profil (od leve stene do sredine kanala) na razdalji 3 m od hipne razširitve; t = 10 s 69 Slika 64: Prečni profil (od leve stene do sredine kanala) na razdalji 4 m od hipne razširitve; t = 10 s 70 Slika 65: Prečni profil (od leve stene do sredine kanala) na razdalji 5 m od hipne razširitve; t = 10 s 70 Slika 66: Kolarjev vrh (Google Earth) 72 Slika 67: Fizični model akumulacije Kolarjev vrh in južne doline pod pregrado (Legiša in Rajar, 1980) 73 Slika 68: Akumulacija Kolarjev vrh s sondami (Džebo in sod., 2013 b) 74 Slika 69: Nihanje gladin v akumulaciji na sondi 1 (Džebo in sod., 2013 b) 75 Slika 70: Nihanje gladin v akumulaciji na sondi 2 (Džebo in sod., 2013 b) 75 Slika 71: Nihanje gladin v akumulaciji na sondi 3 (Džebo in sod., 2013 b) 76 Slika 72: Nihanje gladin v akumulaciji na sondi 4 (Džebo in sod., 2013 b) 76 Slika 73: Nihanje gladin v akumulaciji na sondi 5 (Džebo in sod., 2013 b) 77 Slika 74: Nihanje gladin v akumulaciji na sondi 6 (Džebo in sod., 2013 b) 77 Slika 75: Podolžni profil doline s petimi sondami (Krzyk, 2004) 78 Slika 76: Gladek teren 79 Slika 77: Iztok iz akumulacije (gladko dno – konstanten parameter viskoznosti ob stenah) 80 Slika 78: Nihanje gladin v dolini na sondi 1 (gladko dno – konstanten parameter viskoznosti ob stenah) 81 Slika 79: Nihanje gladin v dolini na sondi 2 (gladko dno – konstanten parameter viskoznosti ob stenah) 81 Slika 80: Nihanje gladin v dolini na sondi 3 (gladko dno – konstanten parameter viskoznosti ob stenah) 82 Slika 81: Nihanje gladin v dolini na sondi 4 (gladko dno – konstanten parameter viskoznosti ob stenah) 82 XX Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 82: Nihanje gladin v dolini na sondi 5 (gladko dno – konstanten parameter viskoznosti ob stenah) 83 Slika 83: Tlorisni prikaz toka vode 27 s po porušitvi pregrade ( d = 2,5 m; bvis = 0,014) 84 Slika 84: Prečni profil; t = 27 s ( d = 2,5 m; bvis = 0,014) 84 Slika 85: Tlorisni prikaz toka vode 62 s po porušitvi pregrade ( d = 2,5 m; bvis = 0,014) 85 Slika 86: Prečni profil; t = 62 s ( d = 2,5 m; bvis = 0,014) 85 Slika 87: Iztok iz akumulacije (gladko dno – spremenljiv parameter viskoznosti ob stenah) 87 Slika 88: Nihanje gladin v dolini na sondi 1 (gladko dno – spremenljiv parameter viskoznosti ob stenah) 88 Slika 89: Nihanje gladin v dolini na sondi 2 (gladko dno – spremenljiv parameter viskoznosti ob stenah) 88 Slika 90: Nihanje gladin v dolini na sondi 3 (gladko dno – spremenljiv parameter viskoznosti ob stenah) 89 Slika 91: Nihanje gladin v dolini na sondi 4 (gladko dno – spremenljiv parameter viskoznosti ob stenah) 89 Slika 92: Nihanje gladin v dolini na sondi 5 (gladko dno – spremenljiv parameter viskoznosti ob stenah) 90 Slika 93: Hrapav teren (dvignjena mrežna vozlišča) 91 Slika 94: Različne gostote hrapavosti 91 Slika 95: Iztok iz akumulacije (hrapavo dno – dvignjena mrežna vozlišča) 94 Slika 96: Nihanje gladin v dolini na sondi 1 (hrapavo dno – dvignjena mrežna vozlišča) 95 Slika 97: Nihanje gladin v dolini na sondi 2 (hrapavo dno – dvignjena mrežna vozlišča) 95 Slika 98: Nihanje gladin v dolini na sondi 3 (hrapavo dno – dvignjena mrežna vozlišča) 96 Slika 99: Nihanje gladin v dolini na sondi 4 (hrapavo dno – dvignjena mrežna vozlišča) 96 Slika 100: Nihanje gladin v dolini na sondi 5 (hrapavo dno – dvignjena mrežna vozlišča) 97 Slika 101: Hrapav teren (točkovni elementi na dnu pobočij) 98 Slika 102: Iztok iz akumulacije (hrapavo dno – podajanje točkovnih hrap) 100 Slika 103: Nihanje gladin v dolini na sondi 1(hrapavo dno – podajanje točkovnih hrap) 101 Slika 104: Nihanje gladin v dolini na sondi 2 (hrapavo dno – podajanje točkovnih hrap) 101 Slika 105: Nihanje gladin v dolini na sondi 3 (hrapavo dno – podajanje točkovnih hrap) 102 Slika 106: Nihanje gladin v dolini na sondi 4 (hrapavo dno – podajanje točkovnih hrap) 102 Slika 107: Nihanje gladin v dolini na sondi 5 (hrapavo dno – podajanje točkovnih hrap) 103 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev XXI Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. LIST OF FIGURES Figure 1: Eulerian (left) and Lagrangian (right) approach 7 Figure 2: Smoothing kernel 15 Figure 3: Water depth calculation with the 2-D Tis Isat model: a) bottom view of the channel b) longitudinal profile of the channel c) smoothing kernel influence (Džebo et al., 2013 a) 26 Figure 4: 2D/3D and 3D/2D coupling (Džebo et al., 2013 a) 28 Figure 5: Point elements on flat terrain 28 Figure 6: Coupling between models of different resolution (Džebo et al., 2012 b) 29 Figure 7: Physical model (Martin and Moyce, 1952) 31 Figure 8: Surge front propagation; Ho/ Wo = 1 (Petkovšek et al., 2010 a) 33 Figure 9: Relative height of the water column; Ho/ Wo = 1 (Petkovšek et al., 2010 a) 33 Figure 10: Surge front propagation; Ho/ Wo = 2 (Petkovšek et al., 2010 a) 34 Figure 11: Relative height of the water column; Ho/ Wo = 2 (Petkovšek et al., 2010 a) 34 Figure 12: Surge front propagation; Ho/ Wo = 4 (Petkovšek et al., 2010 a) 35 Figure 13: Relative height of the water column; Ho/ Wo = 4 (Petkovšek et al., 2010 a) 35 Figure 14: Comparison of fluid propagation 0,08 s after a dam break; Ho/ Wo = 2, number of particles is 50 x 100 (Petkovšek et al., 2010 a) 36 Figure 15: Comparison of fluid propagation 0,3 s after a dam break; Ho/ Wo = 2, number of particles is 50 x 100 (Petkovšek et al., 2010 a) 37 Figure 16: Comparison of fluid propagation 0,44 s after a dam break; Ho/ Wo = 1, number of particles is 50 x 50 (Petkovšek et al., 2010 a) 37 Figure 17: Channel without barrier (top) and channel with barrier (bottom) (Džebo et al., 2011) 38 Figure 18: Surge front propagation in the channel without barrier; t = 1 s (Džebo et al., 2011) 39 Figure 19: Surge front propagatin in channel without barrier (Džebo et al., 2011) 40 Figure 20: Height of the water column in channel without barrier (Džebo et al., 2011) 40 Figure 21: Surge front propagatin in channel with barrier (Džebo et al., 2011) 41 Figure 22: Height of the water column in channel with barrier (Džebo et al., 2011) 42 Figure 23: Laboratory experiment (left) and laboratory experiment with positions of vertical height probes (right) (Kleefsman et al., 2005) 43 XXII Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Figure 24: Comparison between measured and simulated water depths in H4 44 Figure 25: Comparison between measured and simulated water depths in H2 45 Figure 26: Channel with expansion (channel a) and channel with contraction and expansion (channel b) (Rajar, 1972) 46 Figure 27: Phenomena, which occur in such channels: a) channel with expansion b) channel with contraction (Rajar, 1972) 46 Figure 28: Surge front propagation – case 1 (Džebo et al., 2013 a). 49 Figure 29: Comparison between measurements and simulation results 15,5 s after the dam break – case 1 (Džebo et al., 2013 a) 49 Figure 30: Comparison between measurements and simulation results 27 s after the dam break – case 1 (Džebo et al., 2013 a) 50 Figure 31: Comparison between measurements and simulation results 45 s after the dam break – case 1 (Džebo et al., 2013 a) 50 Figure 32: Axonometric view of the part of the channel with an expansion 45 s after the dam break, simulated with the coupled 2D/3D Tis Isat model with larger particles. Pink particles show areas with higher and blue particles show areas with lower water depths (Džebo et al., 2013 a) 51 Figure 33: Axonometric view of the part of the channel with an expansion 45 s after the dam break, simulated with the coupled 2D/3D Tis Isat model with smaller particles. Pink particles show areas with higher and blue particles show areas with lower water depths (Džebo et al., 2013 a) 52 Figure 34: Axonometric view of the part of the channel with an expansion 45 s after the dam break, simulated with the fully 3D Tis Isat model with larger particles. Pink particles show areas with higher and blue particles show areas with lower water depths (Džebo et al., 2013 a) 53 Figure 35: Axonometric view of the part of the channel with an expansion 45 s after the dam break, simulated with the fully 3D Tis Isat model with smaller particles. Pink particles show areas with higher and blue particles show areas with lower water depths (Džebo et al., 2013 a) 54 Figure 36: Cross – section; Y = 32,22 m (Džebo et al., 2013 a) 55 Figure 37: Cross – section; Y = 33,11 m (Džebo et al., 2013 a) 55 Figure 38: Cross – section; Y = 34,00 m (Džebo et al., 2013 a) 56 Figure 39: Surge front propagation – case 2 (Džebo et al., 2013 a) 56 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev XXIII Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Figure 40: Comparison between measurements and simulation results 16,5 s after the dam break – case 2 (Džebo et al., 2013 a) 57 Figure 41: Comparison between measurements and simulation results 35 s after the dam break – case 2 (Džebo et al., 2013 a) 57 Figure 42: Surge front propagation – case 3 (Džebo et al., 2013 a) 58 Figure 43: Comparison between measurements and simulation results 20 s after the dam break – case 3 (Džebo et al., 2013 a) 58 Figure 44: Comparison between measurements and simulation results 30 s after the dam break – case 3 (Džebo et al., 2013 a) 59 Figure 45: Comparison between measurements and simulation results 58 s after the dam break – case 3 (Džebo et al., 2013 a) 59 Figure 46: Surge front propagation – case 4 (Džebo et al., 2013 a) 60 Figure 47: Comparison between measurements and simulation results 19,5 s after the dam break – case 4 (Džebo et al., 2013 a) 60 Figure 48: Comparison between measurements and simulation results 26,5 s after the dam break – case 4 (Džebo et al., 2013 a) 61 Figure 49: Comparison between measurements and simulation results 40 s after the dam break – case 4 (Džebo et al., 2013 a) 61 Figure 50: Channel with expansion (Popovska, 1988) 62 Figure 51: Simulated and measured depth in the point of contracted and expanded channel 63 Figure 52: Cross – section (from the left edge to the middle of the channel) at distance 0,6 m from the expansion ; t = 5 s 64 Figure 53: Cross – section (from the left edge to the middle of the channel) at distance 1 m from the expansion ; t = 5 s 64 Figure 54: Cross – section (from the left edge to the middle of the channel) at distance 2 m from the expansion ; t = 5 s 65 Figure 55: Cross – section (from the left edge to the middle of the channel) at distance 0,6 m from the expansion ; t = 7 s 65 Figure 56: Cross – section (from the left edge to the middle of the channel) at distance 1 m from the expansion ; t = 7 s 66 Figure 57: Cross – section (from the left edge to the middle of the channel) at distance 2 m from the expansion ; t = 7 s 66 XXIV Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Figure 58: Cross – section (from the left edge to the middle of the channel) at distance 3 m from the expansion ; t = 7 s 67 Figure 59: Cross – section (from the left edge to the middle of the channel) at distance 4 m from the expansion ; t = 7 s 67 Figure 60: Cross – section (from the left edge to the middle of the channel) at distance 0,6 m from the expansion ; t = 10 s 68 Figure 61: Cross – section (from the left edge to the middle of the channel) at distance 1 m from the expansion ; t = 10 s 68 Figure 62: Cross – section (from the left edge to the middle of the channel) at distance 2 m from the expansion ; t = 10 s 69 Figure 63: Cross – section (from the left edge to the middle of the channel) at distance 3 m from the expansion ; t = 10 s 69 Figure 64: Cross – section (from the left edge to the middle of the channel) at distance 4 m from the expansion ; t = 10 s 70 Figure 65: Cross – section (from the left edge to the middle of the channel) at distance 5 m from the expansion ; t = 10 s 70 Figure 66: Kolarjev vrh (Google Earth) 72 Figure 67: Physical model of Kolarjev vrh reservoir and south valley below the dam (Legiša and Rajar, 1980) 73 Figure 68: Accumulation of Kolarjev vrh with gauges (Džebo et al., 2013 b) 74 Figure 69: Depth oscillations in accumulation at Gauge 1 (Džebo et al., 2013 b) 75 Figure 70: Depth oscillations in accumulation at Gauge 2 (Džebo et al., 2013 b) 75 Figure 71: Depth oscillations in accumulation at Gauge 3 (Džebo et al., 2013 b) 76 Figure 72: Depth oscillations in accumulation at Gauge 4 (Džebo et al., 2013 b) 76 Figure 73: Depth oscillations in accumulation at Gauge 5 (Džebo et al., 2013 b) 77 Figure 74: Depth oscillations in accumulation at Gauge 6 (Džebo et al., 2013 b) 77 Figure 75: Longitudinal profile with positions of five gauges (Krzyk, 2004) 78 Figure 76: Smooth terrain 79 Figure 77: Outlet from the accumulation (smooth terrain – constant coefficient for wall – particle interaction) 80 Figure 78: Depth oscillations in valley at Gauge 1 (smooth terrain – constant coefficient for wall – particle interaction) 81 Figure 79: Depth oscillations in valley at Gauge 2 (smooth terrain – constant coefficient for wall – particle interaction) 81 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev XXV Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Figure 80: Depth oscillations in valley at Gauge 3 (smooth terrain – constant coefficient for wall – particle interaction) 82 Figure 81: Depth oscillations in valley at Gauge 4 (smooth terrain – constant coefficient for wall – particle interaction) 82 Figure 82: Depth oscillations in valley at Gauge 5 (smooth terrain – constant coefficient for wall – particle interaction) 83 Figure 83: Ground view of water flow 27 s after the dam break ( d = 2,5 m; bvis = 0,014) 84 Figure 84: Cross – section; t = 27 s ( d = 2,5 m; bvis = 0,014) 84 Figure 85: Ground plan of water flow 62 s after the dam break ( d = 2,5 m; bvis = 0,014) 85 Figure 86: Cross – section; t = 62 s ( d = 2,5 m; bvis = 0,014) 85 Figure 87: Outlet from the accumulation (smooth terrain – variable coefficient for wall – particle interaction) 87 Figure 88: Depth oscillations in valley at Gauge 1 (smooth terrain – variable coefficient for wall – particle interaction) 88 Figure 89: Depth oscillations in valley at Gauge 2 (smooth terrain – variable coefficient for wall – particle interaction) 88 Figure 90: Depth oscillations in valley at Gauge 3 (smooth terrain – variable coefficient for wall – particle interaction) 89 Figure 91: Depth oscillations in valley at Gauge 4 (smooth terrain – variable coefficient for wall – particle interaction) 89 Figure 92: Depth oscillations in valley at Gauge 5 (smooth terrain – variable coefficient for wall – particle interaction) 90 Figure 93: Rough terrain (elevated mesh nodes) 91 Figure 94: Different rough density 91 Figure 95: Outlet from the accumulation (rough terrain – elevated mesh nodes) 94 Figure 96: Depth oscillations in valley at Gauge 1 (rough terrain – elevated mesh nodes) 95 Figure 97: Depth oscillations in valley at Gauge 2 (rough terrain – elevated mesh nodes) 95 Figure 98: Depth oscillations in valley at Gauge 3 (rough terrain – elevated mesh nodes) 96 Figure 99: Depth oscillations in valley at Gauge 4 (rough terrain – elevated mesh nodes) 96 Figure 100: Depth oscillations in valley at Gauge 5 (rough terrain – elevated mesh nodes) 97 Figure 101: Rough terrain (point elements on the bottom of the slopes) 98 Figure 102: Outlet from the accumulation (rough terrain – point elements on the bottom of the slopes) 100 XXVI Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Figure 103: Depth oscillations in valley at Gauge 1 (rough terrain – point elements on the bottom of the slopes) 101 Figure 104: Depth oscillations in valley at Gauge 2 (rough terrain – point elements on the bottom of the slopes) 101 Figure 105: Depth oscillations in valley at Gauge 3 (rough terrain – point elements on the bottom of the slopes) 102 Figure 106: Depth oscillations in valley at Gauge 4 (rough terrain – point elements on the bottom of the slopes) 102 Figure 107: Depth oscillations in valley at Gauge 5 (rough terrain – point elements on the bottom of the slopes) 103 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 1 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 1 UVOD 1.1 Splošen opis problema V okviru disertacije smo obravnavali porušitvene valove. Tovrstni valovi nastanejo zaradi porušitve pregrade ter lahko povzročijo ogromno škodo. Gibanje porušitvenih valov opišemo z nelinearnimi parcialnimi diferencialnimi enačbami (St. Venantove enačbe), ki pri poljubnih robnih pogojih niso analitično rešljive in jih rešujemo s pomočjo numeričnih metod. Izbira numerične metode, s pomočjo katere bomo reševali osnovne enačbe, je zelo pomembna. Naša odločitev bo vplivala na točnost in zanesljivost rezultatov, saj vse numerične metode niso enako primerne za vse vrste primerov. Poznamo dva osnovna načina izpeljevanja osnovnih enačb. To sta Eulerjev in Lagrangeov način. V praksi se najpogosteje uporablja Eulerjev način izpeljevanja osnovnih enačb. Ta način je že preverjen in daje dobre rezultate. Pri Eulerjevem pristopu je predmet preučevanja natančno določen del prostora, v katerem se nahaja tekočina. V tem določenem delu prostora se izračunavajo različne fizikalne količine, kot sta hitrost in tlak, ter so funkcija koordinat in časa. Najpopularnejše numerične metode, ki so osnovane na Eulerjevem načinu, so metoda končnih razlik, metoda končnih elementov in metoda končnih volumnov. Največji slabosti teh metod sta numerična difuzija in vpliv mreže na rezultate. Numerična difuzija lahko znatno vpliva na kakovost rezultatov. Njen vpliv lahko zmanjšamo s povečanjem gostote numerične mreže, ne moremo se ji pa popolnoma izogniti. Vpliv mreže na rezultate pa lahko omilimo z uporabo kompleksnejših mrež, katerih priprava je zahtevno in zamudno delo. Tem slabostim konvencionalnih metod pa se lahko izognemo z uporabo brezmrežnih Lagrangeovih metod. Pri Lagrangeovem pristopu obravnavamo masne delce tekočine, ki se premikajo po prostoru. Fizikalne količine (hitrost, tlak, itd.) določenega delca z znanimi začetnimi koordinatami so samo funkcija časa, zato konvekcijskega pospeška (enačba 2, tretji člen na desni strani enačaja), ki določa spremembo hitrosti zaradi spremembe položaja delcev, ni potrebno izračunavati. Najpopularnejša Lagrangeova metoda, ki se najpogosteje uporablja za simulacije tokov s prosto gladino, je metoda hidrodinamike zglajenih delcev ali smoothed particle hydrodynamics (SPH). To je razmeroma nova metoda, ki še ni dobro raziskana in je še vedno v fazi preverjanja. Metoda SPH je tako imenovana metoda delcev, ki je zaradi svoje narave zelo primerna za simulacije hitrih sprememb vodne gladine (npr. simulacije porušitvenih valov). Za razliko od Eulerjevega pristopa, kjer računska mreža prekriva celotno območje, 2 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. kjer bi se voda lahko med simulacijo nahajala, pa v SPH metodi podamo samo lokacijo delcev vode na začetku simulacije. To nam omogoča uporabo večje resolucije delcev, po drugi strani pa premikajoči se delci niso omejeni na velikost računskega območja, ampak se lahko premikajo po celotni domeni. Zgoraj naštete prednosti metode SPH so nam dale spodbudo za raziskovanje in uporabo le-te. Hkrati pa smo tudi pričakovali, da nam bo metoda SPH dala točnejše rezultate za simulacije porušitvenih valov. V okviru disertacije smo to poskusili tudi dokazati na takšen način, da smo rezultate simulacij po metodi SPH primerjali z meritvami na fizičnih modelih ter z rezultati simulacij konvencionalnih modelov. Poleg tega pa do sedaj še ni bila izvedena raziskava simulacij porušitvenega vala na realni topografiji po metodi SPH, kjer bi pridobljene rezultate po metodi SPH primerjali z meritvami ter rezultati drugega konvencionalnega modela. Na omenjenem primeru smo preverili tudi vpliv različne hrapavosti terena na rezultate. Tudi te raziskave še niso bile opravljene. Metodo so razvili Gingold in Monaghan (1977) ter Lucy (1977), prvotno z namenom simulacij astrofizikalnih pojavov v tridimenzijskem (3D) odprtem prostoru. Metoda se je razvijala, izpopolnjevala ter preverjala in sčasoma postala zanimiva tudi za simulacije problemov na drugih področjih. Uporabo metode SPH za simulacije dinamičnega odziva trdnih snovi so predlagali Libersky in Petsheck (1993). Metoda se je izkazala kot zelo primerna za simulacije pojavov, kjer pride do hipnih sprememb, zato je svoje mesto našla tudi pri simulacijah eksplozij (Liu in sod., 2003 a). Metoda je bila prvič predlagana za simulacije toka s prosto gladino leta 1994 (Monaghan, 1994). Monaghan (1994) je prilagodil osnovne enačbe metode SPH, ki so se uporabljale za simulacije stisljivih plinov v astrofiziki, za simulacije toka s prosto gladino. Voda se običajno v hidromehaniki obravnava kot nestisljiva tekočina. Metoda SPH pa je bila razvita za simulacije stisljivih plinov, zato je Monaghan (1994) vodo predpostavil kot tekočino z majhno stisljivostjo, pri čemer razlika v gostoti vode ne sme biti večja kot 1 %. Takšen način obravnavanja vode je še danes v uporabi. Obstajajo tudi razne dopolnitve metode SPH, ki vodo obravnavajo kot nestisljivo, vendar se je s tem računski čas že tako ali tako zamudnih SPH simulacij samo še povečal. Metoda SPH je postala v zelo kratkem času razmeroma priljubljena za simulacije toka s prosto gladino. Opravljene so bile razne študije valovanja ter vpliva le-tega na obalne zgradbe in plovila na morju (Dalrymple in Knio, 2000; Gómez-Gesteira in Dalrymple, 2004; Rogers Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 3 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. in sod., 2010). Z metodo SPH so bile opravljene simulacije delovanja hidrotehničnih objektov (Gatti in sod., 2007, Lopez in sod., 2009; Lopez in sod., 2010), izračunani večfazni tokovi (Colagrossi in Landrini, 2003; Falappi in sod., 2007; Sibilla, 2007), narejene študije porušitve pregrad (Prakash in sod., 2001; Roubtsova in Kahawita, 2006) ter toka ne-Newtonske tekočine (Shao in Lo, 2003). V metodi SPH je predvsem potrebno večjo pozornost posvetiti definiranju robnih pogojev in razvijanju novih pristopov za zniževanje računskega časa izredno zamudnih SPH simulacij. Obstajata dva načina definiranja robnih pogojev po metodi SPH. Prvi način je ta, da se rob predpostavi z navideznimi delci ob robu. Robni pogoj je lahko dinamični z uporabo mirujočih ali gibljivih delcev (Crespo in sod., 2007; Dalrymple in Knio, 2000) ali odbojni (Monaghan, 1994; Monaghan in Kos, 1999). Pri drugem načinu pa se rob izračuna s pomočjo integracije ob robu med stikom tekočine s trdnim telesom (Kulasegaram in sod., 2003; Petkovšek in sod., 2010 a). Najpogosteje se uporablja prva možnost, ki rob nadomesti z navideznimi delci ob robu, pa čeprav ta robni pogoj povzroča nekatere fizikalno neobstoječe pojave ob steni: delce vode potiska stran od roba oz. se le-ti zalepijo na njega (Petkovšek in sod., 2010 a). Integracija ob robu izboljša obnašanje delcev ob steni. Iz tega razloga je razvoj novih robnih pogojev po metodi SPH zelo pomemben, saj ustrezno definirani robovi lahko precej pripomorejo k izboljšavi rezultatov. Pri metodi SPH definicija terena ter določanje hrapavosti stene ni popolnoma samoumevno. Hrapavost sten se definira s pomočjo robnih pogojev. Predlaganih je bilo nekaj načinov definiranja trdnih robov. Prakash in sod. (2001) so simulirali porušitvene valove na realni topografiji, ki bi nastali zaradi porušitve pregrade. Teren so definirali s pomočjo robnih delcev ter interpolacijsko dolžino 10 m (simulacije grobe resolucije) ali 5 m (simulacije fine resolucije). Roubtsova in Kahawita (2006) sta simulirala dobro znano katastrofo ob prelitju pregrade Vaiont, ki se je zgodila leta 1963. Tudi v tem primeru je bil rob definiran s pomočjo delcev ob robu in drsnega robnega pogoja (Roubtsova in Kahawita, 2006). Lee in sod. (2009) so simulirali širjenje valov, ki bi nastali zaradi porušitev rečne pregrade v Gouloirsu v Franciji. Tudi v tem primeru so bili uporabljeni navidezni delci (Violeau in Issa, 2007). Vpliv robnega pogoja na delce tekočine se običajno opiše z modeli turbulence. Najpogosteje se uporablja model umetne viskoznosti, ki nam opiše obnašanje delcev in pri katerem je 4 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. vrednost turbulentne viskoznosti navadno enaka tako za delce tekočine kot za robne delce. Hitrost robnih delcev je običajno enaka vrednosti 0. Poleg razvoja robnih pogojev pa nam skrajševanje časa SPH simulacij predstavlja ogromen izziv. Metoda SPH je izrazito zamudna in zato se razvijajo novi pristopi, ki nam skrajšujejo računski čas simulacij. SPH modele običajno povezujemo z drugimi računsko manj zahtevnimi modeli. Predlagani so bili različni načini povezovanja modelov: povezovanje SPH modelov s konvencionalnimi modeli (Narayanaswamy in sod., 2010), povezovanje med modeli različne resolucije (Børve in sod., 2001; Børve in sod., 2004; Lastiwka in sod., 2005; Feldman in Bonet, 2007; Vacondio in sod., 2013) in povezovanje med modeli različne dimenzionalnosti (Džebo in sod., 2012 a; Džebo in sod., 2013 a). 1.2 Zasnova doktorske disertacije Kot osnovo našega dela smo uporabili model Tis Isat, ki je tipični SPH model, ki vodo obravnava kot malo stisljivo in je bil razvit na Katedri za mehaniko tekočin z laboratorijem. Model uporablja dva različna koeficienta turbulentne viskoznosti, to sta koeficient turbulentne viskoznosti med delci tekočine in koeficient turbulentne viskoznosti med steno in tekočino. Osnovni model je bil preverjen na številnih testnih primerih (Petkovšek in sod., 2010 a in b). Osnovni cilji disertacije so: - razvoj, verifikacija in uporaba metode SPH za račun porušitvenih valov. - nadgradnja modela Tis Isat z novo vrsto robnega pogoja ter preveriti vpliv hrapavosti tal na širjenje poplavnega vala. Robni pogoji predstavljajo izziv vsem raziskovalcem metode SPH, saj robni pogoji niso samoumevni. Poleg tega pa so študije širjenja poplavnega vala v odvisnosti od hrapavosti ostenja redke in bo ta disertacija pripomogla k natančnejšemu razumevanju vpliva hrapavosti tal na širjenje poplavnega vala. - preveriti, kako vrednost parametra viskoznosti ob stenah in velikost delcev vplivata na rezultate simulacij po metodi SPH. - izdelati vmesnik za povečanje resolucije z manjšimi delci na določenih območjih, kjer se iskane količine hitro spreminjajo. Metoda SPH je časovno zelo zamudna. Omenjeni vmesnik nam bo omogočil izvajanje simulacij z manjšimi delci samo na območju, kjer je tok izrazito spremenljiv in so manjši delci za zagotovitev zadostne natančnosti nujno potrebni. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 5 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. - narediti simulacijo primera porušitvenega vala na realni topografiji, kar bo dalo disertaciji poleg raziskovalnega tudi praktični pomen. Zastavljene cilje smo dosegli na naslednji način. Dopolnjen model Tis Isat smo v okviru disertacije najprej preverili na dveh 2D in enem 3D testnem primeru. Rezultate 2D simulacij modela Tis Isat, ki smo jih imeli na voljo, smo najprej primerjali z rezultati drugega SPH modela SPHysics, katerega koda je brezplačno dostopna na internetni strani (https://wiki.manchester.ac.uk/sphysics/index.php/Downloads) ter z meritvami. Na drugem 2D primeru smo preverili, kako se rezultati simulacij izboljšujejo z manjšanjem delcev. Rezultate smo primerjali z meritvami. Nato smo preverili še delovanje 3D matematičnega modela Tis Isat. Rezultate smo primerjali z rezultati drugih dveh SPH modelov (Lee in sod., 2010), rezultati konvencionalnega modela (Kleefsman in sod., 2005) ter z meritvami (Kleefsman in sod., 2005). Ko smo se prepričali, da naš model daje dobre rezultate, smo preverili še vmesnik, ki nam omogoča povezovanje med 2D in 3D različico modelov Tis Isat. Ta vmesnik lahko uporabimo za simulacije toka v kanalih, kjer je tok nekaj časa širinsko povprečen, nato pa se tloris kanala spremeni in povzroči spremembe toka po širini kanala. V tem primeru je možno v ravnem delu kanala uporabiti 2D model, pred območjem tlorisnih razlik pa simulacije nadaljujemo s 3D modelom. Vmesnik smo preverili na primeru toka vode v kanalu z razširitvijo ter v kanalu z zožitvijo in razširitvijo. Omenjeni vmesnik je pozitivno vplival na zmanjšanje računskega časa simulacij. V nadaljevanju smo razvili še en vmesnik, ki nam je omogočil povezovanje med modeli različne resolucije. Tudi v tem primeru vplivamo na računski čas simulacij, saj simulacije izvajamo z večjimi delci na celotnem območju, kjer ni izrazitih nestacionarnih pojavov. Na območju, kjer bi radi dobili natančnejše rezultate, pa delce pomanjšamo in posledično povečamo število le-teh. Vmesnik smo preverili na primeru toka vode v kanalu z zožitvijo in razširitvijo. V nadaljevanju smo simulirali valove, ki bi nastali zaradi morebitne porušitve pregrade zgornje akumulacije črpalne elektrarne Kolarjev vrh in ta primer nam je dokončno potrdil uporabnost modela Tis Isat tudi za praktične primere. V tem primeru smo hrapavost terena definirali na dva načina. Najprej smo predpostavili gladek teren, kjer smo hrapavost terena definirali s pomočjo parametra hrapavosti ob stenah, ki smo ga podali na dva načina, kot konstanten ali spremenljiv v odvisnosti od naklona terena. Nato smo teren opisali kot realno hrapav, kjer smo hrapavost podali s pomočjo dvigovanja mrežnih vozlišč. Model smo dopolnili še z možnostjo podajanja točkovnih elementov na pobočjih ter izvedli simulacije tudi s tako dopolnjenim modelom. 6 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Izkazalo se je, da je naš model primerljiv z ostalimi SPH in konvencionalnimi modeli ter nam v primerjavi z meritvami daje zelo dobre rezultate. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 7 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 2 OSNOVNE ENAČBE IN NUMERIČNE METODE REŠEVANJA Kot smo že omenili, osnovne enačbe rešujemo s pomočjo numeričnih metod. Razdelimo jih v dve skupini, ki se med seboj razlikujeta po tem, na kakšen način diskretiziramo računsko območje. Eulerjev pristop računsko območje določi s pomočjo uporabe numerične mreže, v kateri izračunava različne fizikalne količine (slika 1, levo). Lagrangeova metoda pa računsko območje določi s pomočjo masnih delcev, ki se premikajo po prostoru. V tem primeru se fizikalne količine izračunavajo v samem delcu, ki se premika po prostoru (slika 1, desno). Slika 1: Eulerjev (levo) in Lagrangeov (desno) pristop Figure 1: Eulerian (left) and Lagrangian (right) approach 2.1 Osnovne enačbe Zaradi v uvodu že omenjenih prednosti smo se v okviru doktorske disertacije odločili za uporabo Lagrange-ovega pristopa oz. metode zglajenih delcev SPH. V nadaljevanju najprej podajamo osnovne enačbe hidrodinamike tekočin, nato prirejene enačbe za metodo SPH in postopek njihovega reševanja. 8 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 2.1.1 Predstavitev osnovnih enačb Osnovne enačbe, ki se nahajajo v vseh modelih, ki simulirajo hidrodinamiko tekočine, so kontinuitetna in dinamična enačba ter enačba stanja. Splošna oblika kontinuitetne enačbe je: ௗఘ + ݀݅ݒሺߩݒറሻ = 0 (1) ௗ௧ kjer je ρ gostota vode, ݒറ pa vektor hitrosti vode. V splošnem velja, da je voda praktično nestisljiva (ρ = konst), zato je prvi člen na levi strani enačaja enak 0. Metoda SPH pa obravnava tekočino kot malo stisljivo, zato člena na levi strani enačbe ostane. Splošna oblika dinamične enačbe je: ௗ௩ሬറ = − ଵ ݃ݎܽ݀ ݌ + ܨറ + ߥ∆ݒറ + జ ݃ݎܽ݀ ݀݅ݒ ݒറ (2) ௗ௧ ఘ ଷ kjer je p tlak vode, ܨറ so sile, ki delujejo na obravnavano enoto mase oziroma delec vode, ν pa kinematični koeficient turbulentne viskoznosti. Tlake izračunamo s pomočjo enačbe stanja. Običajno se le-ta v Eulerjevem modelu zapiše kot: ߩሺܶ, ݏ, ݌ሻ = 0 (3) Temperaturo T in slanost s izračunavamo s pomočjo advekcijsko-difuzijske enačbe: డ஼ + ݑ డ஼ = ܦ డమ஼ డ௧ ௜ డ௫ మ (4) ೔ డ௫೔ kjer je C koncentracija, D pa koeficient difuzije. Advekcijsko-difuzijski enačbi, s katerima izračunavamo T in s, se zapišeta kot: డ் + డሺ௨்ሻ + డሺ௩்ሻ + డሺ௪்ሻ = డ ቀܦ డ்ቁ + డ ቀܦ డ்ቁ + డ ቀܦ డ்ቁ (5) డ௧ డ௫ డ௬ డ௭ డ௫ ்ೣ డ௫ డ௬ ்೤ డ௬ డ௭ ்೥ డ௭ డ௦ + డሺ௨௦ሻ + డሺ௩௦ሻ + డሺ௪௦ሻ = డ ቀܦ డ௦ቁ + డ ቀܦ డ௦ቁ + డ ቀܦ డ௦ቁ (6) డ௧ డ௫ డ௬ డ௭ డ௫ ்ೣ డ௫ డ௬ ்೤ డ௬ డ௭ ்೥ డ௭ Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 9 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. kjer so ܦ் , ܦ in ܦ koeficienti turbulentne difuzije v ೣ ்೤ ்೥ x, y in z smeri, u, v in w pa komponente hitrosti v x, y in z smeri. Pri metodi SPH pa se običajno uporablja naslednjo obliko enačbe stanja: ݌ = ܿଶሺߩ − ߩ଴ሻ (7) Kjer je c hitrost zvoka, ρ je gostota, ρ pa začetna gostota tekočine. 0 2.1.2 Modeli turbulence Intenzivnost turbulence vpliva na hitrost, tlak in višino vode. Povezana je z dodatnimi strižnimi silami, ki se pojavijo med posameznimi sloji tekočine ter med robom in tekočino. Vpliv turbulence določimo z modeli turbulence. Poznamo jih več vrst, delimo jih glede na število vsebovanih transportnih enačb za količine, ki karakterizirajo stanje turbulentnosti toka v določeni točki. Različni modeli turbulence so na kratko opisani v nadaljevanju. - model brez transportnih enačb: V tem primeru podamo konstantne vrednosti turbulentnih koeficientov. To je najpreprostejši turbulentni model, ki se zaradi razvoja zmogljivejših računalnikov opušča. - model mešalne dolžine: Vrednosti turbulentne viskoznosti se določijo s pomočjo mešalne dolžine (dolžina mešanja v turbulentnem toku) in gradienta hitrosti. - model Smagorinsky: Tudi ta model je osnovan na eni enačbi. V modelu moramo določiti brezdimenzijski koeficient Smagorinsky, ki pa ni nujno, da je konstanten in se lahko v posameznih točkah računskega področja spreminja. Poleg koeficienta Smagorinsky pa na izračun vrednosti koeficienta turbulentne viskoznosti vplivajo tudi gradienti srednjih hitrosti toka. 10 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. - model z eno transportno enačbo ( k model): Intenzivnost turbulence se določa s turbulentno kinetično energijo ( k). - model z dvema transportnima enačbama (npr. k-ε model): Model z dvema transportnima enačbama je natančnejši od modela z eno transportno enačbo, saj določa intenzivnost turbulence s pomočjo turbulentne kinetične energije ( k) in stopnje disipacije (ε). 2.2 Konvencionalne numerične metode reševanja S pomočjo numeričnih metod rešujemo opisane osnovne nelinearne parcialne diferencialne enačbe toka, za katere za poljubne robne pogoje analitične rešitve ne obstajajo. Poznamo več vrst numeričnih metod za reševanje osnovnih enačb, zapisanih v Eulerjevem načinu, najpopularnejše pa so naslednje tri: metoda končnih elementov, metoda robnih elementov in metoda končnih razlik. - metoda končnih elementov: Računsko področje diskretiziramo na končne prostorske elemente, ki so poljubno oblikovani. Ta metoda se odlično prilagaja zapletenim robovom računskega področja, vendar je časovno zelo zamudna in zahteva relativno veliko spominske kapacitete računalnika. - metoda robnih elementov: Računsko področje pri prostorskem 3D primeru diskretiziramo na površinske ploskovne like. Po načinu reševanja je ta metoda podobna metodi končnih elementov, le da se prostorski elementi nadomestijo s ploskovnimi liki. S tem se potreba po računskem času in predvsem potreba po spominskih kapacitetah računalnikov bistveno zmanjša. - metoda končnih razlik: Ta metoda je relativno enostavna in jasna ter časovno razmeroma hitra. Slabost metode je ta, da se Kartezijev koordinatni sistem slabše prilagaja robovom. Tej pomanjkljivosti se lahko izognemo z vpeljavo novega, bolj kompleksnega krivočrtnega koordinatnega sistema, ki se bolje prilagaja lokaciji robov, vendar je priprava tovrstnih mrež običajno nekoliko zamudnejša. Ker po drugi strani obstaja že veliko izkušenj z reševanjem nelinearnih Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 11 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. konvekcijskih členov v dinamični enačbi, je metoda končnih razlik še vedno ena izmed najpopularnejših konvencionalnih numeričnih metod. Metoda, ki združuje dobre lastnosti metode končnih razlik in metode končnih elementov, je metoda končnih volumnov (Patankar, 1980; Četina, 1988; Četina, 1992). Princip diskretizacije računskega področja je naslednji: Računsko področje se razdeli na kontrolne volumne. Tlak, gostota, temperatura, slanost in koncentracija se običajno izračunavajo v središču volumna, hitrost pa na premaknjenih pozicijah. Osnovne enačbe nato pretvorimo v diferenčne na takšen način, da jih integriramo znotraj kontrolnih volumnov po prostoru in času. Za diskretizacijo pa uporabimo izbrano numerično shemo (na primer: hibridna numerična shema, quick, quickest, ultimate itd.). Pri reševanju osnovnih enačb s pomočjo konvencionalnih numeričnih metod reševanja se pojavi problem numerične difuzije. Določena napaka nastane v vsaki iteraciji aproksimacije konvekcijskih členov in se skozi celoten izračun povečuje. Numerično napako je možno zmanjšati z zgoščevanjem numerične mreže ali uporabo numeričnih shem višjega reda točnosti, vendar je ne moremo popolnoma odpraviti. Lahko pa se ji izognemo z uporabo Lagrangeovih brezmrežnih numeričnih metod, kot je metoda SPH, ki je podrobneje opisana v nadaljevanju. 2.3 SPH metoda reševanja 2.3.1 Predstavitev osnovnih enačb metode SPH Metoda SPH je Lagrangeova brezmrežna numerična metoda. Vrednost poljubne fizikalne količine A poljubnega delca na lokaciji z radijvektorjem ݎറ se izračuna po spodnji enačbi: ܣሺݎറሻ ≈ ׬ ܣ൫ݎఫሬറ൯ܹ൫ݎറ − ݎఫሬറ൯݀ݎఫሬറ (8) ఆೝሬറ kjer ߗ௥റ označuje obravnavano območje okoli delca na lokaciji z radijvektorjem ݎറ, W pa je tako imenovana jedrna funkcija, katere lastnosti so podrobneje opisane v nadaljevanju. V SPH metodi je celoten sistem predstavljen s končnim številom delcev, ki imajo določeno maso in zasedejo določeno mesto v prostoru. Masa delca j je: 12 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. ݉௝ = ∆ܸ௝ߩ௝ (9) kjer je ∆ V končni volumen delca j j ter ߩ௝ gostota delca j. Integracija v enačbi (8) se nato zamenja z vsoto vseh pripadajočih vrednosti sosednjih delcev v dosegu nekega obravnavanega območja, ki ga določa jedrna funkcija W. S tem dobimo diskretizirano obliko enačbe, ta proces diskretizacije pa se imenuje aproksimacija delcev. SPH jedrna aproksimacija se lahko zapiše v naslednji diskretizirani obliki: ܣሺݎറሻ = ׬ ܣ൫ݎఫሬറ൯ܹ൫ݎറ − ݎఫሬറ൯݀ݎఫሬറ (10) ఆ ≈ ∑௝ ܣ൫ݎఫሬറ൯ܹ൫ݎറ − ݎఫሬറ൯∆ܸ௝ = ∑ ௠ೕ ௝ ܣ൫ݎఫሬറ൯ܹ൫ݎറ − ݎఫሬറ൯ ೝሬറ ఘೕ oz. jedrna aproksimacija za delec i na lokaciji z radijvektorjem ݎറ௜ se zapiše kot: ܣሺݎሬపሬറሻ = ∑ ௠ೕ ௝ ܣ൫ݎ ሬറ − ݎ (11) ఘ ఫ ሬറ൯ ∙ ܹ൫ݎప ఫሬറ൯ ೕ Gradienti se izračunajo po spodnji enačbi: grad ܣሺݎറሻ ≈ ׬ ܣ൫ݎఫሬറ൯ ∙ grad ܹ൫ݎറ − ݎఫሬറ൯݀ݎఫሬറ = ׬ ܣ൫ݎఫሬറ൯ ∙ ܹ൫ݎറ − ݎఫሬറ൯݀ܵ − ׬ ܣ൫ݎఫሬറ൯ ∙ grad ܹ൫ݎറ − ݎఫሬറ൯ ݀ݎറ (12) ఆ ௌ ఆ ೝሬറ ೝሬറ Prvi integral na desni strani predstavlja integracijo po robu med stikom tekočine z robom, kjer je S obravnavana površina roba, ki jo določa jedrna funkcija W. Vrednost tega integrala je znotraj tekočine enaka 0. Če se delec nahaja v neposredni bližini roba, pa ta člen predstavlja vpliv roba na obravnavani delec tekočine. Vpeljava diskretizirane oblike enačbe v metodo SPH je pravzaprav ključnega pomena, saj s tem postane metoda SPH uporabna tudi brez vpeljave numerične mreže. Diskretizirana oblika vpeljuje maso in gostoto delcev v samo enačbo. Najverjetneje je ta lastnost metode SPH vplivala na to, da je le-ta postala v zelo kratkem času izredno priljubljena za reševanje hidrodinamičnih problemov, pri katerih je gostota ključnega pomena. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 13 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Lastnosti jedrne funkcije W Znanih je več različnih vrst jedrnih funkcij, vendar imajo vse enake glavne lastnosti, ki sta jih opisala Liu in Liu (2003 a): - Jedrna funkcija mora biti normalizirana skozi celotno obravnavano območje: ׬ ܹ൫ݎറ − ݎఫሬറ൯݀ݎఫሬറ = 1 (13) ఆೝሬറ - Jedrna funkcija naj bo omejena: ܹ൫ݎറ − ݎఫሬറ൯ = 0, za หݎറ − ݎఫሬറห > ݇ ∙ ℎ (14) Omejitev je določena s pomočjo dolžine jedrne funkcije h in s pomočjo normiranega faktorja k ( k določa razporeditev določene jedrne funkcije). Z neenačbo หݎറ − ݎఫሬറห ≤ ݇ ∙ ℎ definiramo obravnavano območje delca, ki se nahaja na lokaciji ݎറ. - Vrednost jedrne funkcije naj bo pozitivna na celotnem obravnavanem območju delca, ki se nahaja na lokaciji ݎറ: ܹ൫ݎറ − ݎఫሬറ൯ ≥ 0 (15) Ta pogoj nam zagotovi fizikalno smiselnost oz. stabilnost nekega fizičnega pojava. Obstajajo tudi takšne jedrne funkcije, ki imajo negativne vrednosti, vendar se pri simulacijah hidrodinamike takšne jedrne funkcije ne uporabljajo, saj bi povzročile nekatere fizično nelogične parametre, kot so negativne vrednosti gostote. - Z naraščanjem razdalje med delci, vrednost jedrne funkcije upada. Bližje so si delci, večji je vpliv med njimi. 14 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. - Jedrna funkcija naj bo soda funkcija (simetrične lastnosti). To pomeni, da imajo vsi delci, ki so enako oddaljeni od obravnavanega delca, enak vpliv na njega, ne glede na njihovo pripadajočo lokacijo. - Jedrna funkcija naj bo gladka, to pomeni, da je zvezno odvedljiva. S tem zagotovimo boljšo aproksimacijo in rezultate. Znano je, da jedrna funkcija z zglajenimi vrednostmi funkcij in odvodov običajno prinese boljše rezultate (Monaghan, 1992; Fulk, 1994). Vse funkcije, ki imajo zgoraj naštete lastnosti, so lahko SPH jedrne funkcije. V nadaljevanju so predstavljene nekatere izmed njih. Lucy (1977) je predlagal uporabo naslednje zvonaste jedrne funkcije: ܹ൫ݎറ − ݎఫሬറ൯ = ܹሺݍሻ = ଵ ൜ሺ1 + 3ݍሻሺ1 − ݍሻଷ ݍ ≤ 1 ௪ 0 ݍ > 1 (16) kjer je ݍ = ௟, ௛ l je razdalja med dvema delcema ter w je normalizacijski faktor. Vrednost normalizacijskega faktorja je ߨℎଶ 5 ⁄ za 2D in 16ߨℎଷ 1 ⁄ 05 za 3D simulacije. Gingold in Monaghan (1977) sta predlagala naslednjo obliko jedrne funkcije: ܹ൫ݎറ − ݎఫሬറ൯ = ܹሺݍሻ = ଵ ݁ି௤మ (17) ௪ kjer je ݓ = ߨℎଶ za 2D in ݓ = ߨଷ ଶ ⁄ ℎଷ za 3D simulacije. Monaghan in Lattanzio (1985) sta predlagala naslednjo jedrno funkcijo: 1 − ଷ ݍଶ + ଷ ݍଷ 0 ≤ ݍ < 1 ଶ ସ ܹ൫ݎറ − ݎఫሬറ൯ = ܹሺݍሻ = ଵ ൞ ଵ ௪ ሺ2 − ݍሻଷ 1 ≤ ݍ < 2 (18) ସ 0 ݍ ≥ 2 kjer je ݓ = 7ߨℎଶ 1 ⁄ 0 za 2D in ݓ = ߨℎଷ za 3D simulacije. Ta oblika jedrne funkcije se najpogosteje uporablja (slika 2). Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 15 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 2: Jedrna funkcija Figure 2: Smoothing kernel Pogosto se uporablja tudi jedrna funkcija, ki jo je predlagal Wendland (1995): ସ ܹ൫ݎറ − ݎఫሬറ൯ = ܹሺݍሻ = ଵ ቀ1 − ௤ቁ ሺ2ݍ + 1ሻ 0 ≤ ݍ ≤ 2 (19) ௪ ଶ kjer je ݓ = 4ߨℎଶ/7 za 2D in ݓ = 8ߨℎଷ/7 za 3D simulacije. Morris (1996) je predlagal funkcijo višjega reda: ۓሺ3 − ݍሻହ − 6ሺ2 − ݍሻହ + 15ሺ1 − ݍሻହ 0 ≤ ݍ < 1 ሺ3 − ݍሻହ − 6ሺ2 − ݍሻହ 1 ≤ ݍ < 2 ܹ൫ݎറ − ݎఫሬറ൯ = ܹሺݍሻ = ଵ (20) ௪ ۔ ሺ3 − ݍሻହ 2 ≤ ݍ < 3 ە 0 ݍ ≥ 3 kjer je ݓ = 478ߨℎଶ 7 ⁄ za 2D in ݓ = 359ߨℎଷ 3 ⁄ za 3D simulacije. Kontinuitetna enačba (1) se v SPH metodi integrira na obeh straneh enačbe znotraj območja dosega delca ݎఫሬറ in pomnoži z jedrno funkcijo: − ׬ ଵ ௗఘ ఆ ܹ൫ݎറ − ݎఫሬറ൯݀ݎఫሬറ = ׬ ݀݅ݒሺݒറሻ ܹ൫ݎറ − ݎఫሬറ൯݀ݎఫሬറ (21) ೝሬറ ఘ ௗ௧ ఆೝሬറ ali 16 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. − ׬ ଵ ௗఘ (22) ఆ ܹ൫ݎറ − ݎఫሬറ൯݀ݎఫሬറ = ׬ ൣ݀݅ݒ൫ܹ൫ݎറ − ݎఫሬറ൯ݒറ൯ − ݒറ ∙ grad ܹ൫ݎറ − ݎఫሬറ൯൧݀ݎఫሬറ ೝሬറ ఘ ௗ௧ ఆೝሬറ Z upoštevanjem Gaussovega teorema dobimo: ଵ ௗఘ = ׬ ݒറ ∙ grad ܹ൫ݎറ − ݎ (23) ఘ ௗ௧ ఫ ሬറ൯݀ݎఫሬറ − ׬ ܹ൫ݎറ − ݎఫሬറ൯ݒറ ∙ ݁റ ݀ܵ ఆೝሬറ ௌ Zadnji integral je površinski integral, kjer je ݁റ smerni vektor, ki je pravokoten na gladino in kaže navzven. Če je delec, ki se nahaja na lokaciji ݎറ daleč stran od gladine, se lahko ta integral zanemari. Dobimo naslednjo obliko kontinuitetne enačbe: ଵ ௗఘ = ׬ ݒറ ∙ grad ܹ൫ݎറ − ݎ (24) ఘ ௗ௧ ఫ ሬറ൯݀ݎఫሬറ ఆೝሬറ ki se v diskretizirani obliki zapiše kot: ଵ ௗఘ೔ = ∑ ௠ೕ ݒሬሬറܹ ఘ పఫ ௜௝′ ௝ (25) ೕ ௗ௧ ఘೕ ali ௗఘ೔ = ∑ ݉ ሬሬറܹ ௗ௧ ௝ݒపఫ ௜௝′ ௝ (26) kjer je v razlika med hitrostima med delcema ij i in j. V primeru, da zanemarimo viskoznost (viskoznost je obravnavana v poglavju 2.3.2), se dinamična enačba v Lagrangeovem pristopu zapiše kot: ௗ௩ሬറ = − ଵ ∙ grad ݌ + ܨറ ௗ௧ ఘ (27) Člen ଵ ∙ grad ݌ se izračuna po spodnji enačbi: ఘ ଵ ∙ grad ݌ = ௣ ∙ grad ߩ + grad ቀ௣ቁ (28) ఘ ఘమ ఘ Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 17 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. oziroma če obe strani enačbe integriramo po celotnem vplivnem območju jedrne funkcije W v okolici delca ݎఫሬറ, pomnožimo z jedrno funkcijo ter nato integral zamenjamo z vsoto vseh vrednosti v dosegu obravnavanega območja, dobimo naslednjo obliko enačbe: ቀ௚௥௔ௗ ௣ቁ = ∑௝ ݉ ఘ ௝ ൬௣೔మ + ௣ೕమ൰ grad௜ ܹ௜௝ (29) ௜ ఘ೔ ఘೕ Z upoštevanjem enačbe 27 dobimo naslednjo diskretizirano obliko dinamične enačbe po metodi SPH: ௗ௩ሬഢሬറ = − ∑ ݉ ௗ௧ ௝ ൬௣೔ ఘమ + ௣ೕమ൰ grad௜ ܹ௜௝ + ܨറ ௝ (30) ೔ ఘೕ Tlake izračunamo s pomočjo spodnje enačbe stanja: ݌ = ܿଶሺߩ − ߩ଴ሻ (31) kjer je c hitrost zvoka in je običajno enaka desetkratniku največje predvidene hitrosti delcev tekočine (Monaghan, 1994). Ta predpostavka izhaja iz naslednje povezave: ܯଶ = ௩మ = ∆ఘ (32) ௖మ ఘబ Če bi podali realno vrednost hitrosti zvoka, bi posledično dobili majhno vrednost Machovega števila M in sprememba gostote bi bila skoraj nepomembna. Ta predpostavka pa velja za nestisljive tekočine. Kot smo že povedali, pri metodi SPH vodo obravnavamo kot malo stisljivo tekočino (~1 %) in zato moramo podati manjšo vrednost hitrosti zvoka od realne vrednosti. Podrobnejši opis in izpeljave osnovnih enačb najdemo v literaturah (Gómez- Gesteira in sod., 2010; Violeau, 2012; Liu in Liu, 2003 b). 2.3.2 Modeli turbulence po metodi SPH Ko se je metoda SPH razvijala, se je turbulenca modelirala s pomočjo t.i. modela umetne viskoznosti (Monaghan, 1992). To je zelo enostaven in preprost model, ki je kljub temu dal zelo dobre rezultate. Dinamična enačba se lahko v SPH obliki zapiše kot: 18 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. ௗ௩ሬഢሬറ = − ∑ ൬௣೔ ௗ௧ ఘమ + ௣ೕమ + ߎ௜௝൰ ݉௝ grad ∙ ܹ′௜௝ + ܨറ ௝ (33) ೔ ఘೕ kjer je ߎ௜௝ viskozni člen, ki je enak: ିఈ௖തഢതണതఓ೔ೕ ݒሬሬറ ݎሬሬറ < 0 ߎ ఘതതതത పఫ పఫ ௜௝ = ൝ ഢണ (34) 0 ݒሬపሬఫሬറ ݎሬపሬఫሬറ > 0 kjer je ߤ௜௝ = ௛௩೔ೕ௥೔ೕ, ݒሬሬറ = ݒሬറ − ݒሬറ, ݎሬሬറ = ݎሬറ − ݎ ௥మ పఫ ప ఫ పఫ ప ఫ ሬറ, ߟଶ = 0,01ℎଶ, α je parameter, katerega ೔ೕାఎమ vrednost se umerja, ܿതపఫ തത je povprečna vrednost hitrosti zvoka ter h dolžina jedrne funkcije. Morris in sod. (1997) so predlagali uporabo modela z laminarno viskoznostjo. Dinamična enačba se zapiše: ௗ௩ሬറ = − ଵ ∙ grad ݌ + ܨറ + ߥ∆ݒറ (35) ௗ௧ ఘ kjer je: ሬሬሬറ∙ ୥୰ୟୢ ௐ ሺߥ∆ݒറሻ ೔ೕ′ ௜ = ∑௝ ݉௝ ቆସజ௥ഢണ ሬሬറ (36) ൫ఘ೔ାఘೕ൯ห௥ሬഢሬണሬറหమ ቇ ݒపఫ ߥ je kinematična viskoznost laminarnega toka (ߥ = 10ି଺ ݉ଶ ݏ ൗ ), člen ଵ ∙ grad ݌ pa se ఘ izračuna po enačbi 28. Gotoh in sod. so leta 2001 predlagali uporabo laminarne viskoznosti v kombinaciji z SPS (sub-particle scale) modelom turbulence. Dinamična enačba je v tem primeru: ௗ௩ሬറ = − ଵ ∙ grad ݌ + ܨറ + ߥ∆ݒറ + ଵ ∙ grad ߬̅ (37) ௗ௧ ఘ ఘ kjer se člen ߥ∆ݒറ izračuna s pomočjo enačbe (36), člen ଵ ∙ grad݌ s pomočjo enačbe 28, ߬̅ pa je ఘ SPS napetostni tenzor. V SPH metodi ga izračunamo kot: Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 19 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. ఛ ೔ೕ = 2ߥ ݇ߜ ܥ ఘ ௧ܵ௜௝ − ଶଷ ௜௝ − ଶଷ ூ߂ଶߜ௜௝หܵ௜௝หଶ (38) kjer je ߬௜௝ SPS napetostni tenzor, ߥ௧ je kinematični koeficient turbulentne viskoznosti in se izračuna po enačbi 39, k je SPS turbulentna kinetična energija, vrednost konstante C = I 0,0066. Uporabila sta standardni Smagorinsky model (Smagorinsky, 1963): ்߭ = ሺܥ௦ߜሻଶ|ܵ̅| (39) kjer je ܥௌ Smagorinsky-jeva konstanta in njena vrednost je 0,12, ߜ je numerični parameter, ki je običajno enak dvakratniku dolžine jedrne funkcije, |ܵ̅| je stopnja deformacije hitrosti, ki se poda kot: |ܵ̅| = ට2ݏ௜௝ ݏ (40) ௜ ௜௝௜ kjer je ݏ௜௝ = − ଵ ൤డ௩೔ + డ௩ೕ൨ (41) ଶ డ௥ೕ డ௥೔ Pozneje so se razvili pravi modeli turbulence SPH metode. k model sta v SPH modelih predlagala Violeau in Issa (2007). Turbulentna viskoznost delca i ࣏ࢀ,࢏ je : మ ߭ ௞೔ ்,௜ = ܥఓ (42) ఌ೔ ݇௜ … turbulentna kinetična energija ߝ௜ … stopnja disipacije k transportna enačba je zapisana kot (Violeau in Issa, 2007): జ ௞ ௗ௞೔ = ∑ ௠ೕ ೅,೔ାజ೅,ೕ ೔ೕ௥೔ೕ ௝ ∙ ݓ + ܲ ௗ௧ ఘ మ ௛′൫ݎ௜௝൯ ௜ − ߝ௜ (43) ೕ ఙೖ ௥೔ೕାఎమ kjer je ݇௜௝ = ݇௜ − ݇௝. ߪ௞ … Schmidtovo število (ߪ௞ = 1) 20 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. ܲ௜ zapišemo kot: మ ܲ ௞೔ ௜ = 2ܥఓ ݏതതത ݏതതത (44) ೔ ఌ పఫ పఫ ೔ ௜ ௜ ali ܲ ௞೔ ௜ = ݉݅݊ ቀඥܥఓ, 2ܥఓ ݏതതത ݏതതത ቁ 2݇ തതത ݏതതത (45) ఌ పఫ పఫ ௜ݏపఫ పఫ ೔ ௜ ௜ ௜ ௜ kjer je: ଷ య మ ൗ ߝ ସ ൗ ௞೔ ௜ = ܥఓ (46) ௅೘,೔ L je mešalna dolžina. Tudi m,i k – ε model sta predlagala (Violeau in Issa, 2007). Enačbe 42 – 45 ostanejo nespremenjene, medtem ko se stopnja disipacije izračuna s pomočjo enačbe 47: జ ఌ ௗఌ೔ = ∑ ௠ೕ ೅,೔ାజ೅,ೕ ೔ೕ௥೔ೕ ௝ ݓ ሺܥ (47) ௗ௧ ఘ మ ௛′൫ݎ௜௝൯ + ఌ೔ ఌଵܲ௜ − ܥఌଶߝ௜ሻ ೕ ఙഄ ௥೔ೕାఎమ ௞೔ kjer je ߝ௜௝ = ߝ௜ − ߝ௝. Vrednosti konstant so: ܥఓ=0,09 ܥ௄=0,09 ߪఌ =1,30 ܥఌ=0,069 ܥఌଵ=1,44 ܥఌଶ=1,92 Vrednost k pa se izračuna po spodnji enačbi: మ ݇௜ = ௨∗,೔ (48) ඥ஼ഋ ݑ∗,௜ … trenjska hitrost delca i Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 21 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 2.3.3 Upoštevanje modelov turbulence in trenja ob ostenju v modelu Tis Isat Model Tis Isat uporablja naslednjo obliko dinamične enačbe: ሬሬሬറ ሬሬሬറ ሬሬሬറ ௗ௩ሬഢሬറ = ∑ ቆ൬− ൬௣೔ ௩ഢണ ௥ഢണ൰݁ሬሬറ + జ′௩ഢണቇ݉ ௗ௧ ఘమ + ௣ೕమ൰ + జ′ పఫ ௝ܹ′௜௝ + ܨറ ௝ (49) ೔ ఘೕ ଷ ௥ሬഢሬണሬറ ௥ሬഢሬണሬറ ห௥ሬഢሬണሬറห kjer je ν ' enako kinematični viskoznosti ν deljeno s povprečno gostoto delcev i in j, p oz je i pj tlak v delcih i oz j, ݎሬଓሬଔሬറ je razdalja med delcema i in j, W' je odvod jedrne funkcije ter ݁റ௜௝ je smerni vektor. Za kinematično viskoznost v primeru turbulentnega toka upoštevamo vrednosti po spodaj opisanem principu. Model Tis Isat izračunava dve vrsti turbulentne viskoznosti in sicer viskoznost med delci (50): ሬሬሬറห ߥ௔ = ܽݒ݅ݏ ∙ ݀ଶ ห௩ഢണ (50) ௟ೌ in viskoznost ob steni (51): ሬሬሬറห ߥ௕ = ܾݒ݅ݏ ∙ ݀ଶ ห௩ഢണ (51) ௟್ kjer je avis brezdimenzijski parameter viskoznosti v tekočini, bvis je brezdimenzijski parameter viskoznosti ob stenah, d je velikost delcev, ݒሬపሬఫሬറ je razlika v hitrosti med obravnavanima delcema, l je razdalja med delcema in je razdalja med delcem in steno. a lb Vrednost parametra avis se običajno giblje okoli 0,01. To vrednost so predlagali Gómez- Gesteira in sod. (2009). Vrednost parametra bvis pa običajno določimo z umerjanjem modela. Ta vrednost je lahko konstantna ali pa variabilna. V okviru disertacije smo preverili konstantne in variabilne vrednosti parametra bvis v odvisnosti od naklona terena. V tem primeru smo variabilne vrednosti parametra bvis določili s pomočjo enačbe, ki jo je izpeljal Krzyk (2004) z uporabo poenostavljene Rickenmannove enačbe (Rickenmann, 1996; Rickenmann, 2000): ݊ = ܫఈ ܥ ൗ (52) 22 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Manningov koeficient hrapavosti n se je izračunaval za vsako računsko celico posebej. I je padec terena na območju obravnavane računske celice, konstanti C in α pa sta bili določeni z umerjanjem matematičnega modela. Konstanta C je bila posredno izračunana po enačbi: ܥ = ܫఈ ݊ ൘ (53) kjer je ܫ povprečni padec terena, ݊ je izbrani povprečni Manningov koeficient hrapavosti, α pa je izbrani eksponent. V okviru študije porušitve pregrade na realni topografiji, smo parameter bvis izpeljali na podoben način: ܾݒ݅ݏ = ܫఈ ܥ ⁄ ௕ (54) kjer se vrednost parametra C izračuna na naslednji način: b ܥ௕ = ܫ̅ఈ ܾ ⁄ ݒ݅ݏ_ܽ (55) kjer je bvis_a izbrana povprečna vrednost parametra viskoznosti ob stenah. Ta parameter se umerja. 2.3.4 Diskretizacija po času Diskretizacija po času v SPH metodi (Gómez-Gesteira, 2009) običajno sledi eni izmed naslednjih dveh shem za izračun vrednosti spremenljivk v naslednjem časovnem koraku: - Shema Verlet: Vrednost spremenljivke A v naslednjem časovnem koraku k+1 se določi s pomočjo spodnje enačbe: ܣ௞ାଵ = ܣ௞ିଵ + ቀௗ஺ቁ 2 ∙ ݀ݐ (56) ௗ௧ ௞ Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 23 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. kjer je A vrednost spremenljivke v prejšnjem časovnem koraku, je odvod k-1 ቀௗ஺ቁ ௗ௧ ௞ spremenljivke v trenutnem časovnem koraku. - Prediktor – Korektor shema: Pri tej shemi pa se vrednost spremenljivke A v naslednjem časovnem koraku k+1 določi na podlagi enačbe: ܣ௞ାଵ = 2 ∙ ܣ௞ାଵ ଶൗ − ܣ௞ (57) kjer je A vrednost spremenljivke v trenutnem časovnem koraku. k 2.3.5 Robni pogoji Verjetno največji izziv v metodi SPH predstavljajo robni pogoji. Obravnavanje robnih pogojev v SPH metodi ni samoumevno. Poznamo dve vrsti definiranja robnih pogojev: a) Rob se nadomesti s pomočjo navideznih delcev ob robu. V tem primeru rob vpliva na tekočino s pomočjo predpisanih robnih količin delcev in se gradienti jedrne funkcije izračunavajo brez upoštevanja integracije po robu (prvi člen na desni strani enačbe 12 odpade). Vpliv robnih delcev na tekočino lahko upoštevamo na dva načina. - Če predpostavimo dinamični robni pogoj (Crespo in sod., 2007; Dalrymple in Knio, 2000), se lastnosti robnih delcev izračunavajo s pomočjo dinamične in kontinuitetne enačbe, vendar se le-ti ne premikajo v odvisnosti od izračunanih hitrosti. Hitrosti robnih delcev so vnaprej določene in so lahko enake 0 (fiksirani robovi) ali pa so znane (gibljivi robovi). - Če predpostavimo odbojni robni pogoj (Monaghan, 1994; Monaghan in Kos, 1999), potem na delce tekočine vplivajo odbojne sile, ki začnejo vplivati na delce, ko se le-ti približajo robu in preprečijo, da bi delci prečkali rob. b) Rob se izračuna s pomočjo integracije med stikom tekočine s trdnim telesom (Kulasegaram in sod., 2003; Petkovšek in sod., 2010 a), kjer predpostavimo nekatere vrednosti na robu. V tem primeru prvi člen na desni strani enačbe 12 ostane in se izračuna s pomočjo integracije. 24 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Hitrost na robu je predpostavljena (ݒറ = 0), tlak se pa izračuna: ݌௦ = ݉ܽݔ ቀ0; ݌௝; ௖௩೛ቁ (58) ௥೗ kjer je ݌௝ tlak delca j, ݒ௣ je pravokotna hitrost, ki kaže od roba obravnavanega delca proti steni. Člen ܿݒ௣ ݎ ⁄ ௟ preprečuje, da bi delci prečkali rob in pride v poštev samo, ko se delci približajo robu (Monaghan in Kos, 1999). Kontinuitetna in dinamična enačba ob steni dobita obliko: ௗఘ ೕ,ೞ = ߩ ሬሬሬറ ݁ሬሬറܹܫ (59) ௗ௧ ௝ݒఫ௦ ఫ௦ ሬሬሬሬറ ሬሬሬറ ሬሬሬറ ሬሬሬറ ௗ௩ണ,ೞ = ቆ൬௣ೞା௣ೕ − జೞ జണೞ ௘ണೞ൰ ݁ሬሬറ − జೞ௩ണೞቇ ܹܫ (60) ௗ௧ ఘ ௡ ೕ ଷ ௟ ௥೗ Indeks s označuje steno, WI je jedrna funkcija ob steni, ݁ሬ௡ሬറ je pravokotnica na rob, ݒ௦ je viskoznost ob steni roba. ܹܫ = ׬ ܹ൫ݎറ − ݎఫሬറ൯݀ܵ (61) ali ۓ 0,7 − ݍଶ ቀ1 − ଷ ݍଶ + ଷ ݍଷቁ 0 ≤ ݍ < 1 ସ ଵ଴ ܹܫሺݍሻ = ଵ (62) ௛ ۔0,8 − ݍଶ ቀ2 − 2ݍ + ଷ ݍଶ − ଵ ݍଷቁ 1 ≤ ݍ < 2 ସ ଵ଴ ە 0 ݍ ≥ 2 Obstaja več predlaganih metod, kako preprečiti nihanja tlakov ob steni. Ena izmed možnosti je opisana v literaturi (Molteni in Colagrossi, 2008). Ta varianta je vgrajena tudi v model Tis Isat: ߩி = 0,4 ∙ ܿ ∙ ℎ ∙ ఘ೔ೕ ∙ ݀ଷ ∙ ܹ ห௥ሬ ௜௝′ (63) ഢሬണ ሬറห kjer je ߩி gostota toka, d velikost delcev, c zvočna hitrost, h pa dolžina jedrne funkcije. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 25 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 2.3.6 Časovni koraki Za doseganje numerične stabilnosti, mora biti velikost časovnega koraka po metodi SPH manjša od: ∆ݐ ≤ ߣ ௗ (64) ௖ kjer je d velikost delcev, c je hitrost zvoka, λ pa konstanta, ki je običajno enaka vrednosti 0,25. 2.3.7 Izračun globine po metodi SPH Metoda SPH je osnovana na delcih, zaradi tega je točnost izračunane globine odvisna od velikosti delcev. Džebo in sod. (2013 a) so predlagali naslednji enačbi za izračun globine vode: ℎௗ௘௣ = ݀ଶ ∑௕ ܹ൫ݎሬ௕ሬറ − ݎఫሬറ൯ (65) za 2D in ℎௗ௘௣ = ݀ଷ ∑௕ ܹ൫ݎሬ௕ሬറ − ݎఫሬറ൯ (66) za 3D simulacije, kjer je ݎሬ௕ሬറ − ݎఫሬറ enako horizontalni razdalji med obravnavano točko na lokaciji z radijvektorjem ݎሬ௕ሬറ in delcem na lokaciji z radijvektorjem ݎఫሬറ. Grafični prikaz 2D izračuna globine po metodi SPH je prikazan na sliki 3. 26 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 3: Izračun globine po metodi modela Tis Isat: a) tlorisni prikaz 2D SPH modela; b) podolžni profil 2D SPH modela; c) izračun vrednosti jedrne funkcije (Džebo in sod., 2013 a) Figure 3: Water depth calculation with the 2-D Tis Isat model: a) bottom view of the channel b) longitudinal profile of the channel c) smoothing kernel influence (Džebo et al., 2013 a) 2.3.8 Natančnost rezultatov Natančnost ujemanja simuliranih rezultatov se lahko določi s pomočjo različnih metod. V okviru doktorske disertacije smo natančnost simuliranih rezultatov določili vizualno ali s pomočjo spodnje enačbe (mean absolute error): ೙ ܯܣܧ = ∑ ห ೔సభ ீ௜ೞ೔೘ೠ೗೔ೝೌ೙ೌିீ௜೘೐ೝೕ೐೙ೌห (67) ௡ Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 27 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 3. RAZVOJ IN DOPOLNITVE MODELA TIS ISAT Numerični model Tis Isat je napisan v programskem jeziku C++ in je tipični SPH model, ki obravnava vodo kot malo stisljivo. Delovno okolje modela Tis Isat je program Scilab 4.1.2. Program Tis Isat je sestavljen iz numeričnega modela tisisat.dll in vmesnika tisisat.sci, preko katerega se program poganja in je napisan za okolje Scilab. Model Tis Isat uporablja jedrno funkcijo (enačba (18)), ki sta jo predlagala Monaghan in Lattanzio (1985) in se pri tovrstnih simulacijah najpogosteje uporablja. Robne pogoje izračunava s pomočjo integracije med stikom tekočine z robom. Diskretizacija po času pa sledi Verletovemu algoritmu (enačba (56)). V model je vgrajen vmesnik, ki omogoča povezovanje med 2D in 3D modeloma Tis Isat. V kanalih, kjer je tok širinsko povprečen (to so na primer relativno široki prizmatični pravokotni kanali), se običajno lahko uporabljajo računsko manj zahtevni 2D modeli. V primeru, da obravnavano območje spreminja tlorisno obliko, pa pridejo v poštev samo 3D modeli. Velikokrat pa se v praksi srečamo s primerom, kjer je tok v kanalu dolgo časa širinsko povprečen in šele nato pride do spremembe tlorisne oblike. Ker je SPH metoda časovno zelo zamudna, bi bila uporaba polnega 3D modela na daljšem oz. večjem računskem območju težko uresničljiva, saj bi bile simulacije časovno predolge. Iz tega razloga bi morali izvajati simulacije z večjimi delci, kar bi negativno vplivalo na točnost rezultatov. V tem primeru je najbolj smiselna uporaba vmesnika za povezovanje modelov Tis Isat različne dimenzionalnosti. Možno je povezati modela v obe smeri, 2D/3D in 3D/2D. Slika 4 prikazuje primer uporabe vmesnika za povezovanje med modeloma različne dimenzionalnosti na primeru toka v kanalu z razširitvijo oz. zožitvijo. Leva stran slike 4 prikazuje primer povezovanja 2D modela s 3D modelom. V ravnem ožjem delu kanala uporabimo 2D model, ki je časovno razmeroma hiter in daje dobre rezultate. Malo pred razširitvijo se začne povezovanje 2D/3D modela. Cel postopek povezovanja se odvija na razdalji cs· dt ( cs je zvočna hitrost, dt pa časovni korak) gorvodno od mesta razširitve kanala. To je najmanjša razdalja, ki je predlagana v modelu Tis Isat in predstavlja razdaljo, ki jo motnja s hitrostjo cs v časovnem koraku dt prepotuje na tem območju. Na tem območju se delci pomnožijo v vrste po celotni širini kanala. Hkrati pa se jim prenesejo lastnosti o hitrostih in tlakih posameznih 28 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. delcev iz 2D modela. Simulacijo nadaljujemo s 3D modelom. Povezovanje lahko poteka tudi v drugi smeri. Desna stran slike 4 prikazuje primer povezovanja 3D modela z 2D modelom. Na območju pred zožitvijo uporabljamo 3D model. Na razdalji cs· dt pred zožitvijo pa do zožitve se začne postopek povezovanja. V tem postopku uporabimo delce, ki se nahajajo na razdalji |0,5 ∙ ݀| od sredine kanala. Model Tis Isat te delce prestavi v os kanala ter jim pripiše podatke o hitrostih in tlakih iz 3D modela. Simulacije se nadaljujejo z 2D modelom. Slika 4: Povezovanje 2D/3D in 3D/2D modela (Džebo in sod., 2013 a) Figure 4: 2D/3D and 3D/2D coupling (Džebo et al., 2013 a) 3.1 Dopolnitev robnih pogojev Upor toku vode je lahko hrapavostni in oblikovni. Hrapavostni upor se v modelu Tis Isat definira s pomočjo parametra viskoznosti med tekočino in steno. Model smo nadgradili z dodatno možnostjo, ki nam omogoča upoštevanje oblikovnega upora toku vode. Na poljubnih lokacijah na terenu lahko postavimo točkovne elemente, ki so med seboj poljubno oddaljeni (slika 5). Ti elementi lahko ponazarjajo npr. razne zarasti ali večje geometrijske ovire na brežinah, katerih vpliv na tok samo s pomočjo hrapavostnega upora ne moremo ustrezno in dovolj natančno definirati. Slika 5: Točkovni elementi na ravnem terenu. Figure 5: Point elements on flat terrain Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 29 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 3.2 Povezave med modeli različne resolucije V okviru doktorske disertacije smo osnovni 3D model Tis Isat dopolnili tudi z vmesnikom, ki nam omogoča povezovanje med modeli različne resolucije. Leva stran slike 6 prikazuje primer povezovanja modela z večjimi delci z modelom z manjšimi delci v kanalu z razširitvijo. Na razdalji cs· dt gorvodno od razširitve pa do same razširitve kanala smo obravnavane delce pomanjšali za 2x ( d/2). Da smo zagotovili obravnavani volumen vode, smo morali povečati število delcev za 8x. Delcem smo prepisali povprečne hitrosti in tlake v odvisnosti od globine. Simulacije smo nadaljevali z manjšimi delci. Desna stran slike 6 prikazuje primer povezovanja modela z manjšimi delci z modelom z večjimi delci. Na razdalji cs·dt gorvodno od zožitve pa do zožitve kanala smo delce povečali za 2x ter posledično število delcev zmanjšali za 8x. Tudi tem delcem smo prepisali povprečne hitrosti in tlake v odvisnosti od globine. Simulacije smo nato nadaljevali z večjimi delci. Slika 6: Povezovanje med modeli različne resolucije (Džebo in sod., 2012 b). Figure 6: Coupling between models of different resolution (Džebo et al., 2012 b) 30 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 4. PRIMERI TESTIRANJA MODELA TIS ISAT NA OSNOVI REZULTATOV S FIZIČNIH MODELOV Rezultati s fizičnih modelov, ki smo jih na podlagi objav iz literature in lastnih meritev uporabili za verifikacijo matematičnega modela, so bili pridobljeni z umerjenimi merilnimi inštrumenti. Pri merjenju pa vedno naredimo tudi neko napako, ki je posledica človeške napake in nenatančnosti instrumenta, zato imajo meritve neko merilno negotovost. Negotovost pomeni razliko med dejansko in izmerjeno vrednostjo meritve oziroma je realna vrednost meritve vedno nekoliko nižja ali višja od izmerjene vrednosti. Poleg tega se moramo zavedati, da tudi fizični modeli ter na njih izvedene meritve, zaradi vpliva modelnega merila ne prikazujejo popolnoma realnega stanja v naravi. Izvajanje meritev na že zgrajenih objektih v naravi pa je razmeroma drago, zahtevno in zato tudi redko, poleg tega pa nas običajno posledice našega posega v naravo zanimajo še predno je objekt zgrajen. Zato je za testiranje matematičnih modelov uporaba meritev s fizičnih modelov kljub nekaterim slabostim, še vedno najprimernejša. 4.1 Porušitev dvodimenzionalnega vodnega stolpca Najprej smo primerjali rezultate modela Tis Isat z rezultati drugega SPH modela SPHysics, katerega koda je brezplačno dostopna na strani (http://wiki.manchester.ac.uk/sphysics/index.php/Main_Page). Študijo smo opravili na primeru porušitve dvodimenzionalnega vodnega stolpca. Rezultati tega laboratorijskega eksperimenta se pogosto uporabljajo za preverjanje SPH modelov. Številni avtorji (Jones in Belton, 2006; Koshizuka in Oka, 1996; Monaghan, 1994; Roubtsova in Kahawita, 2006) so dosegli v začetnih sekundah po porušitvi pregrade dobro ujemanje z meritvami, nato se je točnost slabšala. Meritve na fizičnem modelu (slika 7) sta opravila Martin in Moyce (1952) na primerih porušitev dvodimenzionalnega vodnega stolpca s tremi različnimi razmerji med višino in širino vodnega stolpca Ho/ Wo: 1 ( Ho = 5,715 cm in Wo = 5,715 cm), 2 ( Ho = 11,43 cm in Wo = 5,715 cm) in 4 ( Ho = 11,43 cm in Wo = 2,8575 cm). Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 31 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 7: Fizični model (Martin in Moyce, 1952) Figure 7: Physical model (Martin and Moyce, 1952) Meritve na fizičnem modelu smo primerjali z rezultati simulacij 2D modela Tis Isat (Žagar in sod., 2008) in z rezultati modela SPHysics. Vhodni parametri so prikazani v preglednici 1. V modelu SPHysics smo izbrali enako velikost in število delcev, kot so bili uporabljeni v modelu Tis Isat. V raziskavi smo preverili obe možnosti definiranja robnih pogojev, ki jih ima na razpolago model SPHysics (dinamični in odbojni robni pogoj navideznih delcev ob robu) ter oba algoritma za izračun integracij po času (metodi Verlet in Prediktor – korektor). Preglednica 1: Začetni parametri Table 1: Initial parametres Parameter Kratica Vrednost Model Tis Isat 2500 (50 x 50) ali 10000 (100 x 100) za H / = 1; 0 W0 5000 (50 x 100) ali 20000 (100 x 200) Število delcev np za H / = 2; 0 W0 10000 (50 x 200) ali 40000 (100 x 400) za H / = 4; 0 W0 Velikost delcev d 0,001143 m Parameter viskoznosti avis 0,01 v tekočini 32 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Parameter viskoznosti bvis 0,01 ob stenah Model SPHysics 2500 (50 x 50) ali 10000 (100 x 100) za H / = 1; 0 W0 Število delcev np 5000 (50 x 100) ali 20000 (100 x 200) za H / = 2; 0 W0 10000 (50 x 200) za H / = 4; 0 W0 Velikost delcev d 0,001143 m Viskoznost i_viscos Laminarna viskoznost + SPS Primerjava meritev in računov za različna razmerja H / je podana v slikah 8-13. Prikazani 0 W0 so samo najboljši rezultati, ki smo jih dobili z umerjanjem parametra bvis oz. s primerjanjem simuliranih rezultatov in meritev propagacije čela vala in relativne višine vodnega stolpca pri treh različnih razmerjih med višino in širino vodnega stolpca ( Ho/ Wo) 1, 2 in 4. Vrednosti so normalizirane na začetno višino Ho ali širino Wo vodnega stolpca. Tudi čas je normaliziran na časovno merilo, ki sta ga predlagala Martin in Moyce (1952): ܶ = ݊ݐට݃ ܹ ൗ ݋ (68) kjer je ݊ = ඥܪ݋ ܹ ⁄ ݋, t je dejanski čas in g je težnostni pospešek. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 33 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 8: Propagacija čela vala; Ho/ Wo = 1 (Petkovšek in sod., 2010 a) Figure 8: Surge front propagation; Ho/ Wo = 1 (Petkovšek et al., 2010 a) Slika 9: Relativna višina vodnega stolpca; Ho/ Wo = 1 (Petkovšek in sod., 2010 a) Figure 9: Relative height of the water column; Ho/ Wo = 1 (Petkovšek in sod., 2010 a) 34 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 10: Propagacija čela vala; Ho/ Wo = 2 (Petkovšek in sod., 2010 a) Figure 10: Surge front propagation; Ho/ Wo = 2 (Petkovšek et al., 2010 a) Slika 11: Relativna višina vodnega stolpca; Ho/ Wo = 2 (Petkovšek in sod., 2010 a) Figure 11: Relative height of the water column; Ho/ Wo = 2 (Petkovšek in sod., 2010 a) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 35 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 12: Propagacija čela vala; Ho/ Wo = 4 (Petkovšek in sod., 2010 a) Figure 12: Surge front propagation; Ho/ Wo = 4 (Petkovšek et al., 2010 a) Slika 13: Relativna višina vodnega stolpca; Ho/ Wo = 4 (Petkovšek in sod., 2010 a) Figure 13: Relative height of the water column; Ho/ Wo = 4 (Petkovšek in sod., 2010 a) 36 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Rezultati simulacij modela SPHysics kažejo hitrejše napredovanje čela vala v primerjavi z meritvami in rezultati simulacij modela Tis Isat. Relativna višina vodnega stolpca pa se je z obema SPH modeloma zmanjševala v skladu z meritvami. Slike 14 - 16 prikazujejo primerjavo napredovanja čela vala, simuliranega z modeloma Tis Isat in SPHysics. Model Tis Isat Model SPHysics Slika 14: Primerjava napredovanja čela vala 0,08 s po porušitvi pregrade; Ho/ Wo = 2, število delcev je 50 x 100 (Petkovšek in sod., 2010 a) Figure 14: Comparison of fluid propagation 0,08 s after a dam break; Ho/ Wo = 2, number of particles is 50 x 100 (Petkovšek et al., 2010 a) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 37 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Model Tis Isat Model SPHysics Slika 15: Primerjava napredovanja čela vala 0,3 s po porušitvi pregrade; Ho/ Wo = 2, število delcev je 50 x 100 (Petkovšek in sod., 2010 a) Figure 15: Comparison of fluid propagation 0,3 s after a dam break; Ho/ Wo = 2, number of particles is 50 x 100 (Petkovšek et al., 2010 a) Model Tis Isat Model SPHysics Slika 16: Primerjava napredovanja čela vala 0,3 s po porušitvi pregrade; Ho/ Wo = 2, število delcev je 50 x 100 (Petkovšek in sod., 2010 a) Figure 16: Comparison of fluid propagation 0,44 s after a dam break; Ho/ Wo = 1, number of particles is 50 x 50 (Petkovšek et al., 2010 a) 38 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Model SPHysics pokaže določene netočnosti v simulacijah. Slike napredovanja čela vala z modelom SPHysics so pokazale, da se delci ponekod odlepijo od stene. Netočnosti se pokažejo tudi na površini vode (slika 15). Medtem ko se nekateri delci odlepijo od glavnine in odtečejo naprej (slika 16), drugi ostanejo zalepljeni na robu stene (slika 16). Ti pojavi nimajo fizikalnega ozadja. Kot smo že omenili, model SPHysics rob nadomesti z navideznimi delci na robu, kar povzroči omenjene netočnosti v simulacijah. Z uvedbo robnih pogojev s pomočjo integracije in trenjem ob robu v model Tis Isat pa se tem netočnostim izognemo. Model Tis Isat da vsaj tako dobre rezultate kot model SPHysics, s tem da je obnašanje ob robovih bistveno boljše. Rezultate te študije smo objavili v (Žagar in sod, 2009; Petkovšek in sod., 2010 a in b). V nadaljevanju smo vpliv različne velikosti delcev na rezultate preverili še s simulacijami z 2D matematičnim modelom Tis Isat v kanalu brez in z zaporo (slika 17), kjer smo rezultate primerjali z lastnimi meritvami. V prvem primeru je voda prosto odtekala iz kanala, v drugem pa se je poplavni val odbil od zapore, zato voda ni mogla odteči iz kanala. Slika 17: Kanal brez zapore (zgoraj) in kanal z zaporo (spodaj) (Džebo in sod., 2011). Figure 17: Channel without barrier (top) and channel with barrier (bottom) (Džebo et al., 2011) Rezultate 2D modela Tis Isat smo primerjali z meritvami, izvedenimi v laboratoriju Katedre za mehaniko tekočin UL FGG. Stranske stene kanala so bile steklene. Akumulacija je bila 2,20 m dolga in 0,20 m široka. V različnih vertikalnih vzdolžnih presekih kanala, ki so bili Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 39 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. oddaljeni od prednje stene od 3,1 cm do 18,3 cm, so bile izvedene meritve. Preseki so bili osvetljeni s halogenim svetlobnim virom, ki je omogočil vizualizacijo širitve poplavnega vala v teh presekih (Tropea in sod., 2007). Uporabljena je bila hitroslikovna digitalna kamera Casio EX-F1 (300 slik na sekundo). Meritve so bile posnete v vzdolžnih presekih od pregrade do zapore v kanalu oz. od pregrade 80 cm dolvodno v kanalu. Natančnost posnetkov je bila 1/300 s. Slika 18 prikazuje posnetek laboratorijske porušitve vodnega stolpca v kanalu brez zapore v času t = 1 s po trenutnem dvigu pregrade. Slika 18: Napredovanje čela vala v kanalu brez pregrade; t = 1 s (Džebo in sod., 2011) Figure 18: Surge front propagation in the channel without barrier; t = 1 s (Džebo et al., 2011) V preglednici 2 so prikazani začetni parametri simulacij. Preglednica 2: Začetni parametri simulacij v laboratorijskem kanalu (Džebo in sod,. 2011) Table 2: Initial parametres of simulations in the laboratory channel (Džebo et al., 2011) Parameter Kratica Vrednost 3520 - 352000 (brez zapore); Število delcev np 1540 - 154000 (z zaporo) Velikost delcev d 0,001 – 0,01 m Parameter viskoznosti avis 0,01 v tekočini Parameter viskoznosti bvis 0,005 – 0,03 ob stenah 40 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Kot smo že omenili, so bile meritve narejene v različnih podolžnih profilih kanala, ki so bili različno oddaljeni od stene. Na meritve, ki so bile izvedene bližje steni, ima trenje ob steni večji vpliv, zato je iz slik 19 – 22 razvidno, da meritve, ki so bile posnete bližje steni, kažejo počasnejše napredovanje čela vala. Slika 19: Napredovanje čela vala v kanalu brez zapore (Džebo in sod,. 2011) Figure 19: Surge front propagatin in channel without barrier (Džebo et al., 2011) Slika 20: Višina vodnega stolpca v kanalu brez zapore (Džebo in sod,. 2011) Figure 20: Height of the water column in channel without barrier (Džebo et al., 2011) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 41 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Simulacije so bile izvedene z večjimi ( d = 0,01 m) in manjšimi delci ( d = 0,001 m). Dokazali smo, da rezultati simulacij z manjšimi delci bolje sledijo meritvam v bližini sredine kanala. Slika 20 prikazuje spreminjanje višine vodnega stolpca v kanalu brez zapore. Pokazala je, da simulacije, ki so bile izvedene z uporabo manjših delcev, bolje sledijo meritvam ob zadnji steni kanala. Sliki 21 in 22 prikazujeta napredovanje čela vala in spreminjanje višine vodnega stolpca v odvisnosti od časa v kanalu z zaporo. Preverili smo rezultate simulacij, ki smo jih izvedli z večjimi ( d = 0,01 m) in manjšimi ( d = 0,001 m) delci. Iz slik je razvidno, da dobimo boljše rezultate z uporabo manjših delcev. Dobili smo dobro ujemanje rezultatov modela Tis Isat z meritvami. Slika 21: Napredovanje čela vala v kanalu z zaporo (Džebo in sod,. 2011) Figure 21: Surge front propagatin in channel with barrier (Džebo et al., 2011) 42 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 22: Višina vodnega stolpca v kanalu z zaporo (Džebo in sod,. 2011) Figure 22: Height of the water column in channel with barrier (Džebo et al., 2011) Kot smo pričakovali, smo dobili boljše rezultate z uporabo manjših delcev. S tem pa smo znatno podaljšali računski čas simulacij. Simulacije so bile izvedene na sistemu s štirjedrnim procesorjem Intel® CoreTM i7 s frekvenco 3,33 GHz. Za 7 s dolgo simulacijo v kanalu brez zapore smo potrebovali 8 min ( d = 0,01 m) oziroma 125,5 ur ( d = 0,001 m) računskega časa ter za 5 s dolgo simulacijo v kanalu z zaporo smo potrebovali 4 min ( d = 0,01 m) oziroma 77 ur ( d = 0,001 m) računskega časa. Ugotovili smo, da z 10 x zmanjšanjem velikosti delcev računski čas podaljšamo za približno 1000 x. 4.2 Porušitev tridimenzionalnega vodnega stolpca V nadaljevanju smo izvedli simulacije porušitev tridimenzionalnega vodnega stolpca v kanalu z oviro, za katerega obstajajo rezultati meritev iz literature (slika 23). Tudi to je osnovni primerjalni test, ki nam je še dodatno potrdil delovanje modela Tis Isat. Eksperimenti so bili izvedeni na pomorskem raziskovalnem inštitutu na Nizozemskem (Maritime Research Institute Netherlands - MARIN). Uporabljen je bil kanal z dimenzijami 3,22 x 1,00 x 1,00 m. Volumen vode v akumulaciji je bil 1,22 m3, globina vode pa 0,55 m. V kanalu je bila locirana škatla z dimenzijami 0,161 x 0,403 x 0,161 m. Od štirih vertikalnih sond so bile tri locirane v kanalu, ena pa v akumulaciji. Več informacij o eksperimentu je dostopnih v člankih (Kleefsman in sod., 2005; Issa in Violeau., 2006; Lee in Han., 2010). Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 43 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 23: Laboratorijski eksperiment (levo) ter laboratorijski eksperiment z lociranimi vertikalnimi sondami (desno) (Kleefsman in sod., 2005) Figure 23: Laboratory experiment (left) and laboratory experiment with positions of vertical height probes (right) (Kleefsman et al., 2005) Primerjava rezultatov računa in meritev je bila opravljena v točkah H4 (slika 24) in H2 (slika 25). Rezultate simulacij modela Tis Isat smo primerjali z meritvami, izvedenimi na fizičnem modelu, rezultati numerične metode »Volume of Fluid« (VOF) (Kleefsman in sod., 2005), rezultati drugega SPH modela, ki obravnava tekočino kot malo stisljivo (WCSPH) (Lee in sod., 2010) in še enega SPH modela, ki obravnava tekočino kot nestisljivo (ISPH) (Lee in sod., 2010). VOF metoda je tipična Eulerjeva metoda. Za to študijo so Kleefsman in sod. (2005) uporabili 236 x 76 x 68 mrežnih celic. Rezultati drugega SPH modela, ki tekočino obravnava kot malo stisljivo (WCSPH) oz. nestisljivo (ISPH) so na voljo v članku (Lee in sod., 2010). V obeh primerih je bil robni pogoj definiran s pomočjo navideznih delcev na robu. Vhodni parametri so podani v preglednici 3. 44 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Preglednica 3: Začetni parametri 3D simulacij v kanalu Table 3: Initial parametres of 3D simulations in the channel Parameter Kratica Vrednosti WCSPH ISPH Matematični model (Lee in sod., (Lee in sod., Tis Isat 2010) 2010) Število delcev np 108540 108540 108540 Velikost delca d 0,55/30 m 0,55/30 m 0,55/30 m Dolžina jedrne funkcije h 1,5 d 1,0 d 1,2 d Konstantna viskoznost ν 10-6 m2/s 10-6 m2/s / Parameter viskoznosti v avis / / 0,01 tekočini Parameter viskoznosti ob bvis / / 0,004 stenah Slika 24: Primerjava med merjenimi in simuliranimi globinami vode v točki H4 Figure 24: Comparison between measured and simulated water depths in H4 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 45 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 25: Primerjava med merjenimi in simuliranimi globinami vode v točki H2 Figure 25: Comparison between measured and simulated water depths in H2 Na slikah 24 in 25 je prikazana primerjava med merjenimi in simuliranimi globinami vode v točkah H4 in H2. V točki H4 dobimo v času od 0 do 2,75 s z modelom Tis Isat najboljše ujemanje z meritvami. Po tem času je opazno rahlo odstopanje od meritev, ki ga pokažejo vsi uporabljeni modeli. V točki H2 pa smo z modelom Tis Isat dobili enako dobre rezultate kot ostali avtorji s svojimi modeli. 6 s dolga simulacija, izvedena na sistemu s štirjedrnim procesorjem Intel® CoreTM i7 s frekvenco delovanja 3,33 GHz je trajala 14,5 h. 4.3 Porušitev vodnega stolpca v kanalu z razširitvijo ter v kanalu z zožitvijo in razširitvijo Delovanje povezanega 2-D/3-D modela Tis Isat smo preverili na primeru porušitve vodnega stolpca v kanalu z razširitvijo (slika 26 a, kanal a) ter v kanalu z zožitvijo in razširitvijo (slika 26 b, kanal b). Širina akumulacije in kanala do zožitve oz. razširitve je bila 0,4 m, naklon dna kanala pa je bil 0,087 %. Meritve, ki so bile izvedene na fizičnem modelu in rezultati 1D modela, ki je osnovan po metodi končnih razlik, so dostopni v literaturi (Rajar, 1972; Rajar, 1978). 46 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 26: Kanal z razširitvijo (kanal a) ter kanal z zožitvijo in razširitvijo (kanal b) (Rajar, 1972) Figure 26: Channel with expansion (channel a) and channel with contraction and expansion (channel b) (Rajar, 1972) V takšnih kanalih nastaneta dva tipična pojava (slika 27). V kanalu z razširitvijo (kanal a) se vodna gladina po razširitvi hipno zniža. V območju, kjer je globina nizka, se pojavi deroči tok z zelo neenakomerno hitrostjo vode. Na koncu tega območja pa nastane poševni vodni skok. V kanalu z zožitvijo se del vala odbije od zožene stene in napreduje v nasprotno smer kot tok vode. Slika 27: Pojava, ki nastaneta v tovrstnih kanalih: a) kanal z razširitvijo b) kanal z zožitvijo (Rajar, 1972) Figure 27: Phenomena, which occur in such channels: a) channel with expansion b) channel with contraction (Rajar, 1972) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 47 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Povezani 2-D/3-D model Tis Isat smo preverili v kanalih a in b s suhim in mokrim dnom. Začetne globine vode so prikazane v preglednici 4. Preglednica 4: Začetne globine vode v kanalih a in b s suhim in mokrim dnom (Džebo in sod., 2013 a) Table 4: Water depths in channels a and b with dry and wet bottom (Džebo et al., 2013 a) Začetna globina vode v Začetna globina vode v Primeri Kanal akumulaciji [cm] kanalu [cm] Kanal a - Primer 1 38,88 0,00 suho dno Kanal a - Primer 2 43,72 6,99 mokro dno Kanal b - Primer 3 29,22 0,00 suho dno Kanal b - Primer 4 30,67 5,83 mokro dno Uporabljeni vhodni parametri so prikazani v preglednici 5. Preglednica 5: Parametri simulacij v kanalih a in b (Džebo in sod., 2013 a) Table 5: Parametres of simulations in channels a and b (Džebo et al., 2013 a) Parameter Kratica Vrednost Primer 1 Primer 2 Primer 3 Primer 4 Velikost 0,039 m in d 0,0437 m 0,029 m 0,031 m delcev 0,02 cm Število 8.350 9.252 7.911 Povezan delcev na 2D delcev 2D delcev 2D delcev 2D/3D 10.968 začetku np in in + model 2D delcev simulacije 33.403 46.326 37.773 Tis Isat ( t = 0 s) 2D delcev 3D delcev 3D delcev 48 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 83.663 3D 3D delcev 129.710 3D 187.806 116.969 model in delcev 3D delcev 3D delcev Tis Isat 668,487 3D delcev Povezan 0,001 2D/3D in 0,001 0,001 0,0013 model Parameter 0,0011 Tis Isat viskoznosti bvis ob stenah 3D 0,004 model in 0,003 0,001 0,0015 Tis Isat 0,0041 V vseh štirih primerih so bili delci 10 x manjši kot je bila začetna globina vode v akumulaciji. S tem smo dobili v vseh primerih 10 vrstic delcev vode po višini v akumulaciji. Za primer 1 smo izvedli še simulacije z manjšimi delci. S tem smo želeli pokazati, kako se z manjšanjem velikosti delcev izboljšujejo rezultati in hkrati zvišuje računski čas. S tem smo tudi dokazali, da je razvoj novih pristopov, ki zmanjšujejo računski čas simulacij, pri metodi SPH zelo pomemben. V preglednici je prikazano začetno število delcev. V 3D matematičnem modelu Tis Isat se število delcev ne spreminja skozi celotno simulacijo, medtem ko pa je število delcev v povezanem 2D/3D modelu spremenljivo. Dolžina jedrne funkcije h je bila enaka velikosti delcev. Vrednost parametra viskoznosti v tekočini avis je bila pri vseh primerih 0,01. Vrednosti parametra viskoznosti ob stenah bvis smo določili z umerjanjem (simulirano napredovanje čela vala smo primerjali z meritvami). Primer 1 Rezultate povezanega 2D/3D modela Tis Isat smo primerjali z meritvami, rezultati 3D modela Tis Isat in rezultati 1D mrežnega modela. Slika 28 prikazuje napredovanje čela vala. Vsi modeli se dobro ujemajo z meritvami. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 49 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 28: Napredovanje čela vala – primer 1 (Džebo in sod., 2013 a). Figure 28: Surge front propagation – case 1 (Džebo et al., 2013 a). 15,5 s po porušitvi pregrade dosežemo boljše ujemanje globin z meritvami z modeloma Tis Isat kot pa z 1D mrežnim modelom (slika 29). Ta model prikazuje nekoliko počasnejše napredovanje čela vala, kot to kažejo meritve. 27 s (slika 30) in 45 s (slika 31) po porušitvi pregrade pa dobimo z modeloma Tis Isat enako dobre rezultate, kot z 1D mrežnim modelom. Slika 29: Primerjava med meritvami in rezultati simulacij 15,5 s po porušitvi pregrade – primer 1 (Džebo in sod., 2013 a) Figure 29: Comparison between measurements and simulation results 15,5 s after the dam break – case 1 (Džebo et al., 2013 a) 50 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 30: Primerjava med meritvami in rezultati simulacij 27 s po porušitvi pregrade – primer 1 (Džebo in sod., 2013 a) Figure 30: Comparison between measurements and simulation results 27 s after the dam break – case 1 (Džebo et al., 2013 a) Slika 31: Primerjava med meritvami in rezultati simulacij 45 s po porušitvi pregrade – primer 1 (Džebo in sod., 2013 a) Figure 31: Comparison between measurements and simulation results 45 s after the dam break – case 1 (Džebo et al., 2013 a) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 51 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. SPH modela lahko simulirata hipni padec vodne gladine tik po razširitvi (slike 32 - 35), kar 1D mrežni model ni sposoben simulirati. Delci roza barve prikazujejo območja z višjo, modri delci pa območja z nižjo globino. Slika 32: Aksonometrični pogled dela kanala z razširitvijo 45 s po porušitvi pregrade, simuliran s povezanim 2D/3D modelom Tis Isat z večjimi delci. Delci v roza barvi prikazujejo območja z višjo, modri delci pa območja z nižjo globino (Džebo in sod., 2013 a) Figure 32: Axonometric view of the part of the channel with an expansion 45 s after the dam break, simulated with the coupled 2D/3D Tis Isat model with larger particles. Pink particles show areas with higher and blue particles show areas with lower water depths (Džebo et al., 2013 a) 52 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 33: Aksonometrični pogled dela kanala z razširitvijo 45 s po porušitvi pregrade, simuliran s povezanim 2D/3D modelom Tis Isat z manjšimi delci. Delci v roza barvi prikazujejo območja z višjo, modri delci pa območja z nižjo globino (Džebo in sod., 2013 a). Figure 33: Axonometric view of the part of the channel with an expansion 45 s after the dam break, simulated with the coupled 2D/3D Tis Isat model with smaller particles. Pink particles show areas with higher and blue particles show areas with lower water depths (Džebo et al., 2013 a) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 53 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 34: Aksonometrični pogled dela kanala z razširitvijo 45 s po porušitvi pregrade, simuliran z 3D modelom Tis Isat z večjimi delci. Delci v roza barvi prikazujejo območja z višjo, modri delci pa območja z nižjo globino (Džebo in sod., 2013 a) Figure 34: Axonometric view of the part of the channel with an expansion 45 s after the dam break, simulated with the fully 3D Tis Isat model with larger particles. Pink particles show areas with higher and blue particles show areas with lower water depths (Džebo et al., 2013 a) 54 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 35: Aksonometrični pogled dela kanala z razširitvijo 45 s po porušitvi pregrade, simuliran z 3D modelom Tis Isat z manjšimi delci. Delci v roza barvi prikazujejo območja z višjo, modri delci pa območja z nižjo globino (Džebo in sod., 2013 a) Figure 35: Axonometric view of the part of the channel with an expansion 45 s after the dam break, simulated with the fully 3D Tis Isat model with smaller particles. Pink particles show areas with higher and blue particles show areas with lower water depths (Džebo et al., 2013 a) Slike 32 - 35 lepo prikazujejo pojav, ki se pojavi v kanalih z razširitvijo. Nižja gladina vode se pojavi v sredini kanala, višja pa ob stenah. Iz slik 36 – 38 je razvidno, da simulacije z manjšimi delci dajo boljše rezultate in zato je razvoj tovrstnih pristopov, ki skrajšujejo računski čas simulacij, zelo pomemben. Prečni profili so še dodatno potrdili sposobnost SPH modelov za simulacije hipnih sprememb gladine, ki se pojavijo v kanalu z razširitvijo. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 55 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 36: Prečni profil; Y = 32,22 m (Džebo in sod., 2013 a) Figure 36: Cross – section; Y = 32,22 m (Džebo et al., 2013 a) Slika 37: Prečni profil; Y = 33,11 m (Džebo in sod., 2013 a) Figure 37: Cross – section; Y = 33,11 m (Džebo et al., 2013 a) 56 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 38: Prečni profil; Y = 34,00 m (Džebo in sod., 2013 a) Figure 38: Cross – section; Y = 34,00 m (Džebo et al., 2013 a) Primer 2 Slika 39 prikazuje napredovanje čela vala za drugi primer. Vsi modeli dajo dobre rezultate. Slika 39: Napredovanje čela vala – primer 2 (Džebo in sod., 2013 a) Figure 39: Surge front propagation – case 2 (Džebo et al., 2013 a) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 57 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Sliki 40 in 41 prikazujeta, da vsi uporabljeni modeli dajejo dobre rezultate. Prednost SPH modelov pred 1D mrežnim modelom je v tem, da je sposoben simulirati območje z nižjo globino vode, ki nastane po razširitvi kanala. Slika 40: Primerjava med meritvami in rezultati simulacij 16,5 s po porušitvi pregrade – primer 2 (Džebo in sod., 2013 a) Figure 40: Comparison between measurements and simulation results 16,5 s after the dam break – case 2 (Džebo et al., 2013 a) Slika 41: Primerjava med meritvami in rezultati simulacij 35 s po porušitvi pregrade – primer 2 (Džebo in sod., 2013 a) Figure 41: Comparison between measurements and simulation results 35 s after the dam break – case 2 (Džebo et al., 2013 a) 58 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Primer 3 Slika 42 prikazuje napredovanje čela vala za tretji primer. Vsi modeli se dobro ujemajo z meritvami. Slika 42: Napredovanje čela vala – primer 3 (Džebo in sod., 2013 a) Figure 42: Surge front propagation – case 3 (Džebo et al., 2013 a) Tudi v tem primeru dobimo enako dobre rezultate z vsemi modeli (slike 43 - 45). SPH modela lahko simulirata območje z višjo globino vode, ki nastane zaradi odboja vala od zoženo steno in območje z nižjo gladino vode, ki nastane takoj po razširitvi. Slika 43: Primerjava med meritvami in rezultati simulacij 20 s po porušitvi pregrade – primer 3 (Džebo in sod., 2013 a) Figure 43: Comparison between measurements and simulation results 20 s after the dam break – case 3 (Džebo et al., 2013 a) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 59 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 44: Primerjava med meritvami in rezultati simulacij 30 s po porušitvi pregrade – primer 3 (Džebo in sod., 2013 a) Figure 44: Comparison between measurements and simulation results 30 s after the dam break – case 3 (Džebo et al., 2013 a) Slika 45: Primerjava med meritvami in rezultati simulacij 58 s po porušitvi pregrade – primer 3 (Džebo in sod., 2013 a) Figure 45: Comparison between measurements and simulation results 58 s after the dam break – case 3 (Džebo et al., 2013 a) 60 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Primer 4 Slika 46 prikazuje napredovanje čela vala za četrti primer. Vsi modeli dajo dobre rezultate. Slika 46: Napredovanje čela vala – primer 4 (Džebo in sod., 2013 a) Figure 46: Surge front propagation – case 4 (Džebo et al., 2013 a) V tem primeru pa dobimo zelo dobro ujemanje SPH modelov z meritvami (slike 47 - 49). Prav tako pa SPH modela pokažeta oba pojava, ki se pojavita v tovrstnih kanalih. Slika 47: Primerjava med meritvami in rezultati simulacij 19,5 s po porušitvi pregrade – primer 4 (Džebo in sod., 2013 a) Figure 47: Comparison between measurements and simulation results 19,5 s after the dam break – case 4 (Džebo et al., 2013 a) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 61 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 48: Primerjava med meritvami in rezultati simulacij 26,5 s po porušitvi pregrade – primer 4 (Džebo in sod., 2013 a) Figure 48: Comparison between measurements and simulation results 26,5 s after the dam break – case 4 (Džebo et al., 2013 a) Slika 49: Primerjava med meritvami in rezultati simulacij 40 s po porušitvi pregrade – primer 4 (Džebo in sod., 2013 a) Figure 49: Comparison between measurements and simulation results 40 s after the dam break – case 4 (Džebo et al., 2013 a) Vse simulacije smo izvajali na sistemu s štirjedrnim procesorjem (Intel® CoreTM i7 3,33 GHz). Z uporabo povezanega 2D/3D SPH modela smo računski čas skrajšali za 74 % (primer 62 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 1), 56 % (primer 2), 73 % (primer 3) in 33 % (primer 4) glede na računski čas polnega 3D modela. Izsledke te študije smo objavili v člankih (Džebo in sod., 2012 a; Džebo in sod., 2013 a). 4.4 Porušitev vodnega stolpca v kanalu s hipno razširitvijo Simulirali smo porušitev vodnega stolpca v kanalu z razširitvijo (slika 50). Voda je bila 0,35 m globoka in se je nahajala v 8,0 m dolgi in 1,2 m široki akumulaciji. Po porušitvi pregrade je voda odtekala skozi ožji kanal, ki je bil 4,0 m dolg in 0,4 m širok. Kanal se je hipno razširil na 2,8 m. Slika 50: Kanal z razširitvijo (Popovska, 1988) Figure 50: Channel with expansion (Popovska, 1988) Rezultate simulacij dopolnjenega modela Tis Isat z možnostjo uporabe različne resolucije smo primerjali z meritvami (Popovska, 1988), rezultati 2D mrežnega modela (Popovska, 1988) in rezultati 3D modela Tis Isat. Dvakrat večje delce smo uporabili do razširitve, nato smo delce zmanjšali na polovično velikost in simulacije nadaljevali z manjšimi delci. V preglednici 6 so podani začetni parametri simulacij. V 3D matematičnem modelu Tis Isat so bili delci 20 x manjši, kot je bila globina vode v akumulaciji. V dopolnjenem 3D modelu pa so bili delci v akumulaciji in zoženem kanalu 10 x manjši, kot je bila začetna globina vode v akumulaciji. Tik pred razširitvijo se delci zmanjšajo na polovično velikosti začetnih delcev. Število delcev se v tem delu posledično poveča za 8x. V 3D matematičnem modelu Tis Isat je ostalo število delcev nespremenjeno skozi celotno simulacijo, medtem ko se število delcev v dopolnjenem 3D modelu spreminja. Dolžina jedrne funkcije h je bila enaka velikosti delcev. Vrednost parametra viskoznosti v tekočini avis je bila 0,01. Vrednosti parametra viskoznosti ob stenah smo določili z umerjanjem (rezultate simulacij smo primerjali z meritvami). Dopolnjeni 3D model z možnostjo uporabe različne resolucije uporablja poleg manjših tudi večje delce. Po Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 63 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. enačbi (51) je trenje poleg vrednosti parametra bvis odvisno tudi od kvadrata velikosti delcev. Zato je umerjena vrednost parametra dodatne viskoznosti ob stenah v dopolnjenem 3D modelu nižja, kot je v osnovnem 3D modelu. Preglednica 6: Začetni parametri pri primeru hipne razširitve kanala (Džebo in sod., 2012 b) Table 6: Initial parametres in the case of sudden expansion of the channel (Džebo et al., 2012 b) Parameter Kratica Vrednost 3D model Tis Dopolnjen 3D model Tis Isat Isat Število delcev np 630660 76062 Velikost delcev d 0,0175 m 0,035 m; 0,0175 m Parameter viskoznosti ob bvis 0,001 0,0005 stenah V nadaljevanju smo vizualno primerjali simulirane in merjene globine v osi kanala v točki zožitve in razširitve v odvisnosti od časa (slika 51) ter globine v različnih prečnih prerezih kanala v različnih časih po porušitvi pregrade (slike 52 - 65). Vsi modeli pokažejo razmeroma dobro ujemanje z meritvami, le prvi dve sekundi po porušitvi pregrade pokaže 2D mrežni model v točki zožitve preveliko globino vode (slika 51). Zožitev Razširitev Slika 51: Simulirana in merjena globina v točki zožitve in razširitve kanala Figure 51: Simulated and measured depth in the point of contracted and expanded channel 64 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slike 52 - 54 prikazujejo prečne profile v času t = 5 s od porušitve pregrade. 0,6 m od hipne razširitve dobimo z SPH modeloma slabše rezultate kot z mrežnim 2D modelom. Na razdalji x = 1,0 m od razširitve pokaže 2D mrežni model preveliko globino med levo steno in sredino razširjenega kanala. SPH modela pa pokažeta preveliko globino v sredini kanala in premalo ob levi steni kanala. Na razdalji 2 metrov od razširitve so rezultati SPH modelov boljši kot rezultati mrežnega modela. Slika 52: Prečni profil (od leve stene do sredine kanala) na razdalji 0,6 m od razširitve; t = 5 s. Figure 52: Cross – section (from the left edge to the middle of the channel) at distance 0,6 m from the expansion ; t = 5 s Slika 53: Prečni profil (od leve stene do sredine kanala) na razdalji 1 m od razširitve; t = 5 s. Figure 53: Cross – section (from the left edge to the middle of the channel) at distance 1 m from the expansion ; t = 5 s Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 65 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 54: Prečni profil (od leve stene do sredine kanala) na razdalji 2 m od razširitve; t = 5 s. Figure 54: Cross – section (from the left edge to the middle of the channel) at distance 2 m from the expansion ; t = 5 s Slike 55 - 59 prikazujejo prečne profile v času t = 7 s od porušitve pregrade. Na razdalji 0,6 m od razširitve kaže mrežni model boljše ujemanje z meritvami na sredini kanala. Modela Tis Isat dosežeta veliko boljše ujemanje z meritvami na razdaljah 1, 2, 3 in 4 m od razširitve. Ob levi steni dobimo z dopolnjenim modelom Tis Isat boljše rezultate kot z osnovnim modelom Tis Isat. Slika 55: Prečni profil (od leve stene do sredine kanala) na razdalji 0,6 m od razširitve; t = 7 s. Figure 55: Cross – section (from the left edge to the middle of the channel) at distance 0,6 m from the expansion ; t = 7 s 66 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 56: Prečni profil (od leve stene do sredine kanala) na razdalji 1 m od hipne razširitve; t = 7 s Figure 56: Cross – section (from the left edge to the middle of the channel) at distance 1 m from the expansion ; t = 7 s Slika 57: Prečni profil (od leve stene do sredine kanala) na razdalji 2 m od hipne razširitve; t = 7 s. Figure 57: Cross – section (from the left edge to the middle of the channel) at distance 2 m from the expansion ; t = 7 s Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 67 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 58: Prečni profil (od leve stene do sredine kanala) na razdalji 3 m od hipne razširitve; t = 7 s Figure 58: Cross – section (from the left edge to the middle of the channel) at distance 3 m from the expansion ; t = 7 s Slika 59: Prečni profil (od leve stene do sredine kanala) na razdalji 4 m od hipne razširitve; t = 7 s Figure 59: Cross – section (from the left edge to the middle of the channel) at distance 4 m from the expansion ; t = 7 s Tudi 10 s po porušitvi pregrade na razdalji 0,6 m od razširitve modela SPH ne kažeta zadovoljivega ujemanja z meritvami (slika 60). 2D mrežni model kaže, v primerjavi z 68 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. meritvami, preveliko globino vode v vseh prečnih profilih. SPH modela kažeta boljše rezultate na razdalji 1, 2, 3, 4 in 5 m od razširitve (slike 61 - 65). Slika 60: Prečni profil (od leve stene do sredine kanala) na razdalji 0,6 m od hipne razširitve; t = 10 s Figure 60: Cross – section (from the left edge to the middle of the channel) at distance 0,6 m from the expansion ; t = 10 s Slika 61: Prečni profil (od leve stene do sredine kanala) na razdalji 1 m od hipne razširitve; t = 10 s Figure 61: Cross – section (from the left edge to the middle of the channel) at distance 1 m from the expansion ; t = 10 s Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 69 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 62: Prečni profil (od leve stene do sredine kanala) na razdalji 2 m od hipne razširitve; t = 10 s. Figure 62: Cross – section (from the left edge to the middle of the channel) at distance 2 m from the expansion ; t = 10 s Slika 63: Prečni profil (od leve stene do sredine kanala) na razdalji 3 m od hipne razširitve; t = 10 s Figure 63: Cross – section (from the left edge to the middle of the channel) at distance 3 m from the expansion ; t = 10 s 70 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 64: Prečni profil (od leve stene do sredine kanala) na razdalji 4 m od hipne razširitve; t = 10 s Figure 64: Cross – section (from the left edge to the middle of the channel) at distance 4 m from the expansion ; t = 10 s Slika 65: Prečni profil (od leve stene do sredine kanala) na razdalji 5 m od hipne razširitve; t = 10 s Figure 65: Cross – section (from the left edge to the middle of the channel) at distance 5 m from the expansion ; t = 10 s Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 71 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Simulacije smo izvedli na sistemu s štirjedrnim procesorjem Intel® CoreTM i7 3,33 GHz. Z osnovnim modelom Tis Isat smo potrebovali 43 h in 40 minut za 11 s dolgo simulacijo ter 24 h in 30 min z dopolnjenim modelom Tis Isat. Ugotovimo lahko, da z uporabo dopolnjenega modela Tis Isat praktično nismo vplivali na točnost rezultatov, smo pa računski čas skrajšali za 44 %. 72 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 5. UPORABA MODELA TIS ISAT ZA PRAKTIČNI PRIMER PORUŠITVE NASIPA ZGORNJE AKUMULACIJE ČRPALNE ELEKTRARNE KOLARJEV VRH Nazadnje smo uporabnost modela Tis Isat preverili še na primeru porušitvenega vala na realni topografiji. Simulirali smo porušitev nasipa zgornje akumulacije črpalne elektrarne Kolarjev vrh. Njena izgradnja je bila predvidena že leta 1979 in čeprav do danes še ni bila izvedena, postaja v zadnjem času ideja o njeni izgradnji spet aktualna. Slika 66: Kolarjev vrh (Google Earth) Figure 66: Kolarjev vrh (Google Earth) Predvideno je, da bo akumulacija postavljena na vrh izravnanega Kolarjevega vrha (slika 66). Površina akumulacije bo 160 000 m2 in v njej naj bi bilo 3,1 x 106 m3 vode. Kota dna akumulacije bo locirana na koti 975 m nadmorske višine, najvišja globina vode v njej bo 20 m. Krona nasipa bo široka 20 m in locirana 1,5 m nad najvišjo globino vode, torej na koti 996,50 m n. m. Na notranji strani akumulacije bo znašal naklon brežin 1:2, na zunanji pa 1:1,5. Na voljo smo imeli meritve, izvedene na fizičnem modelu (slika 67), ki je bil zgrajen v Vodnogospodarskem inštitutu iz Ljubljane (Legiša in Rajar, 1980) in rezultate simulacij 2D matematičnega modela PCFLOW2D-ORTHOCURVE (Krzyk, 2004). Fizični model je bil zgrajen v merilu 1:200 in je obravnaval samo južno dolino pod akumulacijo, ki je bila Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 73 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. najkrajša. Na fizičnem modelu so bile locirane merske sonde, ki so beležile višino vode v odvisnosti od časa in s pomočjo katerih smo lahko dobili višino vala in čas njegovega napredovanja. Slika 67: Fizični model akumulacije Kolarjev vrh in južne doline pod pregrado (Legiša in Rajar, 1980) Figure 67: Physical model of Kolarjev vrh reservoir and south valley below the dam (Legiša and Rajar, 1980) 5.1 Akumulacija Najprej smo v uporabljenem 3D modelu SPH umerili parameter viskoznosti ob stenah v akumulaciji. Kot že omenjeno smo imeli na voljo meritve, pridobljene na fizičnem modelu (Legiša in Rajar., 1980). Slika akumulacije s položajem merskih sond je prikazana na sliki 68. 74 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 68: Akumulacija Kolarjev vrh s sondami (Džebo in sod., 2013 b) Figure 68: Accumulation of Kolarjev vrh with gauges (Džebo et al., 2013 b) Slike 69 – 74 prikazujejo nihanje gladin na šestih sondah v akumulaciji. Izvajali smo simulacije z večjimi ( d = 5 m) in manjšimi ( d = 2,5 m) delci. Število delcev je bilo 21890 ( d = 5 m) oz. 174884 ( d = 2,5 m). Vrednost parametra viskoznosti v tekočini je bila avis_aku = 0,01. Parameter viskoznosti ob stenah v akumulaciji bvis_aku smo umerili na takšen način, da smo simulirano globino primerjali z merjeno globino vode. Po umerjanju se je izkazalo, da rezultati z vrednostjo parametra bvis_aku = 0,001 večinoma najbolje sledijo meritvam. Pri simulacijah z večjimi delci na sondi 3 pa se izkaže, da najbolje sledijo meritvam izračuni z vrednostjo parametra bvis_aku = 0,03. Simulacije z manjšimi delci pa, prav tako na sondi 3, bolje sledijo meritvam pri vrednosti parametra bvis_aku = 0,001. Poleg tega ta simulacija pokaže naglo zvišanje gladine v času 100 s po porušitvi dela nasipa, prav tako kot to pokažejo meritve. Po pričakovanjih pokažejo slike, ki prikazujejo primerjavo časovnega poteka izračunanih in merjenih gladin, na vseh sondah boljše ujemanje pri simulacijah, izvedenih z manjšimi delci. Iz zgoraj naštetih razlogov smo pri vseh nadaljnjih simulacijah predpostavili, da je vrednost parametra viskoznosti ob stenah v akumulaciji bvis_aku = 0,001. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 75 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 69: Nihanje gladin v akumulaciji na sondi 1 (Džebo in sod., 2013 b) Figure 69: Depth oscillations in accumulation at Gauge 1 (Džebo et al., 2013 b) Slika 70: Nihanje gladin v akumulaciji na sondi 2 (Džebo in sod., 2013 b) Figure 70: Depth oscillations in accumulation at Gauge 2 (Džebo et al., 2013 b) 76 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 71: Nihanje gladin v akumulaciji na sondi 3 (Džebo in sod., 2013 b) Figure 71: Depth oscillations in accumulation at Gauge 3 (Džebo et al., 2013 b) Slika 72: Nihanje gladin v akumulaciji na sondi 4 (Džebo in sod., 2013 b) Figure 72: Depth oscillations in accumulation at Gauge 4 (Džebo et al., 2013 b) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 77 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 73: Nihanje gladin v akumulaciji na sondi 5 (Džebo in sod., 2013 b) Figure 73: Depth oscillations in accumulation at Gauge 5 (Džebo et al., 2013 b) Slika 74: Nihanje gladin v akumulaciji na sondi 6 (Džebo in sod., 2013 b) Figure 74: Depth oscillations in accumulation at Gauge 6 (Džebo et al., 2013 b) 78 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 5.2 Dolina dolvodno od porušenega dela nasipa Nato smo preverili še dogajanje v dolini dolvodno od porušenega dela nasipa. Na fizičnem modelu je bilo postavljenih pet merskih sond (slika 75). Slika 75: Podolžni profil doline s petimi sondami (Krzyk, 2004) Figure 75: Longitudinal profile with positions of five gauges (Krzyk, 2004) V preglednici 7 so prikazani najboljši rezultati, ki so bili pridobljeni z modelom PCFLOW2D – ORTHOCURVE (Krzyk, 2004) z uporabo konstantne ter spremenljive vrednosti Manningovega koeficienta hrapavosti n (Krzyk, 2004). Vrednost spremenljivega koeficienta n je bila odvisna od naklona terena, kar prikazuje enačba 52. Natančnost ujemanja simuliranih časov prihoda vala do posameznih sond z merjenimi rezultati pa smo določili s pomočjo enačbe 67 ( MAE). Preglednica 7: Časi napredovanja čela vala; PCFLOW2D – ORTHOCURVE Table 7: Propagation time of the surge front; PCFLOW2D – ORTHOCURVE Eksperiment Sonda 1 Sonda 2 Sonda 3 Sonda 4 Sonda 5 MAE meritve 30 s 62 s 99 s 134 s 175 s / n = konst. = 0,05 24 s 54 s 96 s 132 s 178 s 4,4 n = 0,045; α = 0,2 26 s 56 s 96 s 132 s 176 s 3,2 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 79 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. V SPH modelu Tis Isat smo nato uporabili dva načina definiranja hrapavosti terena. Najprej smo teren opisali kot gladek in nato kot hrapav. 5.2.1 Gladek teren Najprej smo teren obravnavali kot gladek (slika 76). Pri tem načinu smo hrapavost dna definirali s pomočjo parametra viskoznosti ob stenah bvis. Slika 76: Gladek teren Figure 76: Smooth terrain a) Konstantna vrednost parametra viskoznosti ob stenah Tako kot pri umerjanju vrednosti parametra viskoznosti ob stenah v akumulaciji bvis_aku smo tudi pri umerjanju parametra viskoznosti ob stenah bvis za tok dolvodno pojav simulirali z večjimi ( d = 5 m) in manjšimi ( d = 2,5 m) delci. Rezultati so prikazani v preglednici 8. Preglednica 8: Časi napredovanja čela vala; model Tis Isat (gladek teren, konstanten parameter viskoznosti ob stenah) Table 8: Propagation time of the surge front; Tis Isat model (smooth terrain – constant coefficient for wall – particle interaction) Eksperiment Sonda 1 Sonda 2 Sonda 3 Sonda 4 Sonda 5 MAE meritve 30 s 62 s 99 s 134 s 175 s / d = 5,0 m bvis = 0,0070 26 s 59 s 95 s 131 s 176 s 3,0 bvis = 0,0073 27 s 59 s 96 s 133 s 178 s 2,6 80 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. bvis = 0,0074 27 s 59 s 96 s 134 s 179 s 2,6 bvis = 0,0075 27 s 59 s 97 s 134 s 179 s 2,4 bvis = 0,0076 27 s 59 s 97 s 135 s 179 s 2,6 bvis = 0,0080 28 s 60 s 99 s 139 s 181 s 3,0 d = 2,5 m bvis = 0,010 24 s 53 s 90 s 128 s 175 s 6,0 bvis = 0,014 25 s 57 s 99 s 137 s 179 s 3,4 bvis = 0,015 27 s 58 s 101 s 138 s 181 s 3,8 bvis = 0,016 27 s 59 s 103 s 139 s 183 s 4,6 bvis = 0,020 28 s 61 s 106 s 145 s 190 s 7,2 Najboljše rezultate dobimo pri velikosti delcev d = 5,0 m in vrednosti parametra bvis = 0,0075, kjer je vrednost MAE = 2,4 ter pri velikosti delcev d = 2,5 m in vrednosti parametra bvis = 0,014, kjer je vrednost MAE = 3,4. Zato smo pri teh dveh vrednostih parametra bvis preverili ujemanje simuliranih globin v posameznih sondah v primerjavi z merjenimi globinami ter globinami, pridobljenimi z modelom PCFLOW2D – ORTHOCURVE (slike 78 - 82). Krivulje smo še dodatno zgladili s pomočjo uporabe MO Excell-ove funkcije TREND(). Primerjali smo tudi simuliran pretok z modelom Tis Isat z izračunanim pretokom z modelom PCFLOW2D – ORTHOCURVE na iztoku iz akumulacije (slika 77). Slika 77: Iztok iz akumulacije (gladko dno – konstanten parameter viskoznosti ob stenah) Figure 77: Outlet from the accumulation (smooth terrain – constant coefficient for wall – particle interaction) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 81 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 78: Nihanje gladin v dolini na sondi 1 (gladko dno – konstanten parameter viskoznosti ob stenah) Figure 78: Depth oscillations in valley at Gauge 1 (smooth terrain – constant coefficient for wall – particle interaction) Slika 79: Nihanje gladin v dolini na sondi 2 (gladko dno – konstanten parameter viskoznosti ob stenah) Figure 79: Depth oscillations in valley at Gauge 2 (smooth terrain – constant coefficient for wall – particle interaction) 82 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 80: Nihanje gladin v dolini na sondi 3 (gladko dno – konstanten parameter viskoznosti ob stenah) Figure 80: Depth oscillations in valley at Gauge 3 (smooth terrain – constant coefficient for wall – particle interaction) Slika 81: Nihanje gladin v dolini na sondi 4 (gladko dno – konstanten parameter viskoznosti ob stenah) Figure 81: Depth oscillations in valley at Gauge 4 (smooth terrain – constant coefficient for wall – particle interaction) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 83 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 82: Nihanje gladin v dolini na sondi 5 (gladko dno – konstanten parameter viskoznosti ob stenah) Figure 82: Depth oscillations in valley at Gauge 5 (smooth terrain – constant coefficient for wall – particle interaction) Potrebno je poudariti, da delamo z razmeroma velikimi delci. Simulacije z manjšimi delci bi bile še dolgotrajnejše. Na sondi 1 dobimo dobre rezultate v prvih 80 s po porušitvi pregrade. Po 80 s nam simulacije z modelom Tis Isat v primerjavi z meritvami pokažejo prenizko globino vode. Med 120 in 170 s pokažejo simulacije z manjšimi delci nekoliko višjo gladino kot simulacije z večjimi delci. Na sondi 2 se pojavi voda nekoliko kasneje, kot to kažejo meritve. Razlog za takšno dogajanje je v tem, da delci vode ne potujejo po sredini doline, ampak voda nekaj časa teče po bregu, nato se približa sondi (sliki 85 in 86). Simulacije z manjšimi delci se bolje ujemajo z meritvami. Podobno dogajanje se pojavi tudi na sondi 1 (sliki 83 in 84). Na sondi 3 pokažeta obe simulaciji z modelom Tis Isat nekoliko nižjo gladino kot meritve. Na sondah 4 in 5 pa dobimo dobro ujemanje z meritvami. Model PCFLOW2D- ORTHOCURVE pokaže na prvih štirih sondah previsoko gladino vode, ki se nato naglo zniža. Iz teh primerov je razvidno, da simulacije z manjšimi delci dajo bolj točne rezultate. 84 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 83: Tlorisni prikaz toka vode 27 s po porušitvi pregrade ( d = 2,5 m; bvis = 0,014) Figure 83: Ground view of water flow 27 s after the dam break ( d = 2,5 m; bvis = 0,014) Slika 84: Prečni profil; t = 27 s ( d = 2,5 m; bvis = 0,014) Figure 84: Cross – section; t = 27 s ( d = 2,5 m; bvis = 0,014) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 85 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 85: Tlorisni prikaz toka vode 62 s po porušitvi pregrade ( d = 2,5 m; bvis = 0,014) Figure 85: Ground plan of water flow 62 s after the dam break ( d = 2,5 m; bvis = 0,014) Slika 86: Prečni profil; t = 62 s ( d = 2,5 m; bvis = 0,014) Figure 86: Cross – section; t = 62 s ( d = 2,5 m; bvis = 0,014) 86 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. b) Parameter viskoznosti ob stenah v odvisnosti od naklona terena Parameter viskoznosti ob stenah je odvisen od naklona terena (poglavje 2.3.3). Preglednica 9: Čas napredovanja čela vala; Tis Isat (gladko dno – spremenljiv parameter viskoznosti ob stenah) Table 9: Propagation time of the surge front; Tis Isat (smooth terrain – variable coefficient for wall – particle interaction) Eksperiment Sonda 1 Sonda 2 Sonda 3 Sonda 4 Sonda 5 MAE meritve 30 s 62 s 99 s 134 s 175 s / d = 5,0 m bvis = 0,0073; α = 0,0001 27 s 58 s 96 s 134 s 179 s 2,8 bvis = 0,0073; α = 0,001 27 s 58 s 96 s 134 s 178 s 2,6 bvis = 0,0073; α = 0,01 27 s 58 s 96 s 133 s 178 s 2,8 bvis = 0,0073; α = 0,1 27 s 58 s 96 s 133 s 178 s 2,8 bvis = 0,0073; α = 0,2 27 s 58 s 96 s 132 s 176 s 2,6 bvis = 0,0073; α = 0,3 28 s 59 s 97 s 132 s 178 s 2,4 bvis = 0,0073; α = 0,4 28 s 59 s 96 s 133 s 179 s 2,6 bvis = 0,0073; α = 0,5 26 s 57 s 95 s 131 s 174 s 3,4 bvis = 0,0073; α = 0,8 26 s 57 s 95 s 130 s 172 s 4,0 bvis = 0,0073; α = 0,9 26 s 57 s 95 s 129 s 172 s 4,2 bvis = 0,0073; α = 1 26 s 57 s 95 s 129 s 173 s 4,0 bvis = 0,0073; α = 2 24 s 56 s 99 s 135 s 175 s 2,6 bvis = 0,0073; α = 3 24 s 56 s 102 s 137 s 179 s 4,4 bvis = 0,008; α = 0,01 28 s 59 s 100 s 141 s 182 s 4,0 bvis = 0,008; α = 0,1 27 s 59 s 99 s 137 s 180 s 2,8 bvis = 0,008; α = 0,2 27 s 59 s 100 s 139 s 181 s 3,6 bvis = 0,008; α = 0,3 27 s 59 s 98 s 137 s 180 s 3,0 bvis = 0,008; α = 0,4 27 s 59 s 97 s 135 s 180 s 2,8 bvis = 0,008; α = 0,5 29 s 61 s 102 s 142 s 183 s 4,2 bvis = 0,01; α = 0,3 29 s 61 s 102 s 142 s 183 s 4,2 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 87 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. bvis = 0,02; α = 0,3 29 s 61 s 102 s 142 s 183 s 4,2 bvis = 0,03; α = 0,3 29 s 61 s 102 s 142 s 183 s 4,2 d = 2,5 m bvis = 0,014; α = 0,3 25 s 56 s 97 s 136 s 178 s 3,6 bvis = 0,015; α = 0,2 25 s 58 s 100 s 138 s 180 s 3,8 bvis = 0,015; α = 0,3 26 s 58 s 100 s 137 s 179 s 3,2 bvis = 0,015; α = 0,4 26 s 57 s 100 s 137 s 179 s 3,4 bvis = 0,015; α = 0,6 25 s 57 s 99 s 137 s 178 s 3,2 bvis = 0,015; α = 0,8 25 s 57 s 99 s 137 s 178 s 3,2 bvis = 0,016; α = 0,3 26 s 58 s 102 s 138 s 182 s 4,4 bvis = 0,020; α = 0,3 27 s 60 s 105 s 144 s 189 s 7,0 Najboljše rezultate pri simulacijah z večjimi delci ( d = 5 m) dobimo pri vrednosti parametra bvis = 0,0073 ter α = 0,3, kjer je vrednost MAE = 2,4 ter pri simulacijah z manjšimi delci ( d = 2,5 m) in vrednosti parametra bvis = 0,015 ter α = 0,3, kjer je vrednost MAE = 3,2 in največje odstopanje rezultatov od meritev manjše ali enako 4 s (preglednica 9). Simulirane globine omenjenih dveh primerov smo primerjali z meritvami in simuliranimi globinami, pridobljenimi z modelom PCFLOW2D – ORTHOCURVE (slike 88 - 92). Primerjali smo tudi simulirane pretoke na iztoku iz akumulacije (slika 87). Slika 87: Iztok iz akumulacije (gladko dno – spremenljiv parameter viskoznosti ob stenah) Figure 87: Outlet from the accumulation (smooth terrain – variable coefficient for wall – particle interaction) 88 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 88: Nihanje gladin v dolini na sondi 1 (gladko dno – spremenljiv parameter viskoznosti ob stenah) Figure 88: Depth oscillations in valley at Gauge 1 (smooth terrain – variable coefficient for wall – particle interaction) Slika 89: Nihanje gladin v dolini na sondi 2 (gladko dno – spremenljiv parameter viskoznosti ob stenah) Figure 89: Depth oscillations in valley at Gauge 2 (smooth terrain – variable coefficient for wall – particle interaction) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 89 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 90: Nihanje gladin v dolini na sondi 3 (gladko dno – spremenljiv parameter viskoznosti ob stenah) Figure 90: Depth oscillations in valley at Gauge 3 (smooth terrain – variable coefficient for wall – particle interaction) Slika 91: Nihanje gladin v dolini na sondi 4 (gladko dno – spremenljiv parameter viskoznosti ob stenah) Figure 91: Depth oscillations in valley at Gauge 4 (smooth terrain – variable coefficient for wall – particle interaction) 90 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 92: Nihanje gladin v dolini na sondi 5 (gladko dno – spremenljiv parameter viskoznosti ob stenah) Figure 92: Depth oscillations in valley at Gauge 5 (smooth terrain – variable coefficient for wall – particle interaction) Tudi pri tem primeru, kjer je teren gladek in parameter viskoznosti ob stenah odvisen od naklona terena, dobimo precej podobne rezultate kot pri prejšnjem primeru, kjer je bil teren gladek in parameter bvis konstanten. Pri obeh primerih z gladkim terenom je pri večjih delcih najboljša vrednost MAE enaka 2,4, medtem ko pa je pri manjših delcih najboljša vrednost MAE za prvi primer (konstanten parameter bvis) enaka 3,4, za drugi (parameter bvis v odvisnosti od naklona) pa 3,2. Iz tega sklepamo, da se z uporabo parametra viskoznosti ob stenah v odvisnosti od naklona terena rezultati le malo izboljšajo. Tudi nihanje gladin v dolini na sondah 1, 2, 3, 4 in 5 je precej podobno prejšnjemu primeru. 5.2.2 Hrapav teren V nadaljevanju smo predpostavili, da je teren realno hrapav. Hrapavost smo opisali na dva načina. V prvem primeru smo dvigovali mrežna vozlišča dna pobočij (Žagar in sod., 2009) v drugem pa podajali točkovne elemente na dnu pobočij. V obeh primerih je dolina ostala gladka. Takšna hrapavost terena je bila definirana tudi na fizičnem modelu, kjer je bil gozd nad dnom doline ponazorjen s kamenčki različnih velikost, dno doline pa je bilo narejeno z gladkim betonom. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 91 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. a) Dvigovanje mrežnih vozlišč dna pobočij Hrapav teren smo najprej opisali s pomočjo dvigovanja mrežnih vozlišč (slika 93). Tudi na tem primeru smo uporabili dve velikosti delcev, in sicer d = 5,0 m in d = 2,5 m. Slika 93: Hrapav teren (dvignjena mrežna vozlišča) Figure 93: Rough terrain (elevated mesh nodes) Časi napredovanja čela vala do posameznih sond so podani v preglednici 10. Uporabili smo različne vrste hrapavosti (slika 94) in višine hrap ehe. Slika 94: Različne gostote hrapavosti Figure 94: Different rough density 92 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Preglednica 10: Čas napredovanja čela vala; Tis Isat (hrapav teren, dvignjena mrežna vozlišča) Table 10: Propagation time of the surge front; Tis Isat (rough terrain – elevated mesh nodes) Eksperiment Sonda 1 Sonda 2 Sonda 3 Sonda 4 Sonda 5 MAE meritve 30 s 62 s 99 s 134 s 175 s / d = 5,0 m Vrsta hrapavosti: 1 bvis = 0,006; ehe =0,50 m 24 s 54 s 89 s 123 s 163 s 9,4 bvis = 0,007; ehe =0,50 m 25 s 57 s 92 s 126 s 168 s 6,4 bvis = 0,001; ehe =2,00 m 24 s 51 s 96 s 125 s 157 s 9,4 bvis = 0,002; ehe =2,00 m 24 s 54 s 95 s 127 s 164 s 7,2 bvis = 0,001; ehe =2,20 m 25 s 55 s 100 s 135 s 172 s 3,4 bvis = 0,002; ehe =2,20 m 25 s 57 s 103 s 142 s 181 s 5,6 bvis = 0,003; ehe =2,00 m 25 s 57 s 98 s 137 s 178 s 3,4 bvis = 0,003; ehe =2,01 m 25 s 57 s 99 s 138 s 179 s 3,6 bvis = 0,003; ehe =2,02 m 25 s 57 s 99 s 137 s 179 s 3,4 bvis = 0,003; ehe =2,03 m 25 s 57 s 100 s 138 s 180 s 4,0 Vrsta hrapavosti: 2 bvis = 0,003; ehe =2,00 m 25 s 55 s 97 s 129 s 173 s 4,2 bvis = 0,003; ehe =2,05 m 24 s 55 s 98 s 135 s 176 s 3,2 bvis = 0,003; ehe =2,10 m 25 s 57 s 98 s 138 s 180 s 4,0 Vrsta hrapavosti: 3 bvis = 0,003; ehe =2,00 m 25 s 57 s 99 s 137 s 178 s 3,2 bvis = 0,003; ehe =2,05 m 25 s 57 s 100 s 141 s 181 s 4,8 Vrsta hrapavosti: 4 bvis = 0,003; ehe =2,00 m 25 s 57 s 97 s 135 s 177 s 3,0 bvis = 0,003; ehe =2,03 m 25 s 57 s 98 s 137 s 180 s 3,8 Vrsta hrapavosti: 5 bvis = 0,003; ehe =2,00 m 25 s 57 s 97 s 134 s 176 s 2,6 bvis = 0,003; ehe =2,03 m 24 s 58 s 99 s 141 s 181 s 4,6 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 93 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Vrsta hrapavosti: 6 bvis = 0,003; ehe =2,00 m 24 s 55 s 98 s 133 s 176 s 3,2 bvis = 0,003; ehe =2,03 m 24 s 55 s 99 s 137 s 178 s 3,8 Vrsta hrapavosti: 7 bvis = 0,003; ehe =2,00 m 24 s 56 s 98 s 130 s 175 s 3,4 bvis = 0,003; ehe =2,10 m 24 s 58 s 99 s 141 s 181 s 4,6 Vrsta hrapavosti: 8 bvis = 0,006; ehe =0,50 m 24 s 54 s 89 s 121 s 163 s 9,8 bvis = 0,007; ehe =1,00 m 24 s 54 s 95 s 132 s 173 s 4,4 bvis = 0,001; ehe =2,00 m 26 s 56 s 101 s 137 s 169 s 4,2 d = 2,5 m Vrsta hrapavosti: 5 bvis =0,0001; ehe=2,00 m 20 s 46 s 87 s 116 s 146 s 17,0 bvis =0,0005; ehe=2,00 m 20 s 47 s 88 s 121 s 149 s 15,0 bvis = 0,001; ehe =2,00 m 21 s 49 s 89 s 120 s 150 s 14,2 bvis = 0,001; ehe =2,30 m 22 s 52 s 95 s 124 s 167 s 8,0 bvis = 0,001; ehe =2,50 m 24 s 56 s 99 s 130 s 173 s 3,6 bvis = 0,001; ehe =2,55 m 23 s 57 s 99 s 140 s 178 s 4,2 bvis = 0,001; ehe =2,57 m 23 s 57 s 100 s 139 s 175 s 3,6 bvis = 0,001; ehe =2,58 m 23 s 58 s 100 s 141 s 180 s 4,8 bvis = 0,001; ehe =2,59 m 23 s 57 s 99 s 134 s 175 s 2,4 bvis = 0,001; ehe =2,60 m 23 s 57 s 101 s 139 s 175 s 3,8 bvis = 0,001; ehe =2,61 m 23 s 57 s 103 s 143 s 180 s 6,0 bvis = 0,001; ehe =2,62 m 23 s 58 s 103 s 142 s 178 s 5,2 bvis = 0,001; ehe =3,00 m 24 s 64 s 112 s 146 s 183 s 8,2 bvis = 0,001; ehe =4,00 m 31 s 68 s 120 s 153 s 192 s 12,8 bvis = 0,002; ehe =2,00 m 21 s 51 s 92 s 130 s 166 s 8,0 bvis = 0,003; ehe =2,30 m 22 s 55 s 98 s 138 s 181 s 5,2 bvis = 0,004; ehe =2,00 m 22 s 54 s 96 s 134 s 173 s 4,2 bvis = 0,005; ehe =2,00 m 23 s 55 s 96 s 135 s 177 s 4,0 bvis = 0,005; ehe =2,05 m 23 s 55 s 97 s 135 s 179 s 4,2 bvis = 0,005; ehe =2,10 m 23 s 56 s 98 s 137 s 183 s 5,0 bvis = 0,006; ehe =2,00 m 23 s 56 s 97 s 137 s 182 s 5,0 94 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. bvis = 0,007; ehe =2,00 m 24 s 57 s 99 s 140 s 188 s 6,0 bvis = 0,008; ehe =1,00 m 21 s 50 s 89 s 123 s 166 s 10,2 bvis = 0,008; ehe =2,00 m 24 s 58 s 102 s 144 s 193 s 8,2 Najboljše rezultate pri simulacijah z večjimi delci smo dobili pri vrsti hrapavosti 5, vrednosti parametra bvis = 0,003 in višini hrap ehe = 2 m. Vrednost MAE je pri tej simulaciji enaka 2,6. Ker nam je ta simulacija dala najboljše rezultate, smo izvajali simulacije z manjšimi delci samo pri tej vrsti hrapavosti. Najboljše rezultate smo dobili pri vrednosti bvis = 0,001 in višini hrap ehe = 2,59. Pri teh pogojih je vrednost MAE enaka 2,4. Globine vode v posameznih sondah so prikazane na slikah 96 – 100, iztok iz akumulacije za ta dva primera pa na sliki 95. Slika 95: Iztok iz akumulacije (hrapavo dno – dvignjena mrežna vozlišča) Figure 95: Outlet from the accumulation (rough terrain – elevated mesh nodes) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 95 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 96: Nihanje gladin v dolini na sondi 1 (hrapavo dno – dvignjena mrežna vozlišča) Figure 96: Depth oscillations in valley at Gauge 1 (rough terrain – elevated mesh nodes) Slika 97: Nihanje gladin v dolini na sondi 2 (hrapavo dno – dvignjena mrežna vozlišča) Figure 97: Depth oscillations in valley at Gauge 2 (rough terrain – elevated mesh nodes) 96 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 98: Nihanje gladin v dolini na sondi 3 (hrapavo dno – dvignjena mrežna vozlišča) Figure 98: Depth oscillations in valley at Gauge 3 (rough terrain – elevated mesh nodes) Slika 99: Nihanje gladin v dolini na sondi 4 (hrapavo dno – dvignjena mrežna vozlišča) Figure 99: Depth oscillations in valley at Gauge 4 (rough terrain – elevated mesh nodes) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 97 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 100: Nihanje gladin v dolini na sondi 5 (hrapavo dno – dvignjena mrežna vozlišča) Figure 100: Depth oscillations in valley at Gauge 5 (rough terrain – elevated mesh nodes) Na sondi 1 pokaže simulacija z manjšimi delci zelo dobro ujemanje z meritvami, saj nam model Tis Isat pokaže podobno nihanje gladine kot so to pokazale meritve. Tudi na sondi 2 se rezultati, v primerjavi z prejšnjimi primeri, bistveno izboljšajo. To smo tudi pričakovali, saj hrapavo dno zadržuje vodo v dolini in ji ne dovoli, da bi odtekala po bregovih. Na ostalih sondah dobimo dobro ujemanje z meritvami. b) Podajanje točkovnih elementov V nadaljevanju smo podajali točkovne elemente na dnu pobočij (slika 101). Simulacije smo izvajali z uporabo dveh različnih velikosti delcev ( d = 5,0 m in d = 2,5 m). 98 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 101: Hrapav teren (točkovni elementi na dnu pobočij) Figure 101: Rough terrain (point elements on the bottom of the slopes) Pridobljeni rezultati so podani v preglednici 11. Preglednica 11: Čas napredovanja čela vala z modelom Tis Isat (hrapav teren - točkovni elementi) Table 11: Propagation time of the surge front; Tis Isat (rough terrain - point elements on the bottom of the slopes) Eksperiment Sonda 1 Sonda 2 Sonda 3 Sonda 4 Sonda 5 MAE meritve 30 s 62 s 99 s 134 s 175 s / d = 5,0 m Točkovne hrape so postavljene na razdaljah 25 m bvis = 0,001; ht = 0,50 m 22 s 47 s 85 s 114 s 145 s 17,4 bvis = 0,003; ht = 0,50 m 24 s 53 s 97 s 132 s 177 s 4,2 bvis = 0,003; ht = 1,00 m 24 s 55 s 103 s 139 s 185 s 6,4 bvis = 0,003; ht = 1,50 m 26 s 57 s 103 s 145 s 196 s 9,0 bvis = 0,004; ht = 0,50 m 24 s 55 s 103 s 139 s 184 s 6,2 bvis = 0,005; ht = 0,50 m 25 s 57 s 103 s 142 s 187 s 6,8 Točkovne hrape so postavljene na razdaljah 50 m bvis = 0,001; ht = 2,50 m 22 s 42 s 71 s 100 s 135 s 26,0 bvis = 0,003; ht = 2,50 m 24 s 51 s 89 s 125 s 162 s 9,8 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 99 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. bvis = 0,005; ht = 2,50 m 25 s 56 s 95 s 131 s 171 s 4,4 bvis = 0,005; ht = 3,00 m 25 s 57 s 98 s 134 s 176 s 2,4 bvis = 0,005; ht = 3,50 m 26 s 58 s 100 s 137 s 178 s 3,0 bvis = 0,005; ht = 4,00 m 26 s 58 s 100 s 138 s 178 s 3,2 bvis = 0,0055; ht = 2,50 m 26 s 58 s 99 s 136 s 180 s 3,0 bvis = 0,0055; ht = 5,00 m 26 s 58 s 101 s 140 s 181 s 4,4 bvis = 0,006; ht = 2,50 m 26 s 58 s 100 s 138 s 180 s 3,6 bvis = 0,006; ht = 5,00 m 26 s 58 s 103 s 141 s 183 s 5,4 bvis = 0,0065; ht = 2,50 m 27 s 59 s 103 s 140 s 182 s 4,6 bvis = 0,0073; ht = 0,00 m 26 s 59 s 96 s 137 s 180 s 3,6 bvis = 0,0073; ht = 2,50 m 27 s 60 s 106 s 142 s 187 s 6,4 bvis = 0,0073; ht = 5,00 m 28 s 61 s 106 s 144 s 190 s 7,0 Točkovne hrape so postavljene na razdaljah 100 m bvis = 0,005; ht = 2,50 m 26 s 55 s 91 s 122 s 162 s 8,8 bvis = 0,006; ht = 2,50 m 26 s 57 s 93 s 127 s 168 s 5,8 bvis = 0,007; ht = 2,00 m 27 s 59 s 96 s 131 s 176 s 2,6 bvis = 0,007; ht = 2,50 m 28 s 60 s 96 s 132 s 178 s 2,4 bvis = 0,007; ht = 3,00 m 28 s 60 s 99 s 133 s 178 s 1,6 bvis = 0,0071; ht = 2,50 m 28 s 60 s 96 s 133 s 179 s 2,4 bvis = 0,0072; ht = 2,50 m 27 s 59 s 96 s 133 s 178 s 2,6 bvis = 0,0074; ht = 2,50 m 27 s 59 s 97 s 135 s 180 s 2,8 bvis = 0,0076; ht = 2,50 m 27 s 59 s 99 s 137 s 180 s 2,8 bvis = 0,008; ht = 2,50 m 27 s 60 s 100 s 140 s 181 s 3,6 bvis = 0,01; ht = 2,50 m 29 s 62 s 108 s 146 s 190 s 7,4 d = 2,5 m Točkovne hrape so postavljene na 100 m bvis = 0,001; ht = 3,00 m 19 s 37 s 57 s 88 s 123 s 35,2 bvis = 0,0015; ht = 3,00 m 20 s 38 s 59 s 91 s 123 s 33,8 bvis = 0,002; ht = 3,00 m 20 s 40 s 62 s 93 s 128 s 31,4 bvis = 0,01; ht = 3,00 m 24 s 54 s 90 s 129 s 173 s 6,0 bvis = 0,011; ht = 3,00 m 24 s 54 s 92 s 130 s 174 s 5,2 bvis = 0,012; ht = 1,00 m 24 s 55 s 94 s 132 s 175 s 4,0 bvis = 0,012; ht = 3,00 m 25 s 55 s 93 s 132 s 175 s 4,0 100 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. bvis = 0,013; ht = 0,90 m 25 s 57 s 97 s 136 s 177 s 3,2 bvis = 0,013; ht = 1,00 m 25 s 57 s 97 s 135 s 177 s 3,0 bvis = 0,013; ht = 2,90 m 25 s 57 s 97 s 136 s 177 s 3,2 bvis = 0,013; ht = 3,00 m 26 s 57 s 97 s 135 s 177 s 2,8 bvis = 0,013; ht = 3,10 m 26 s 57 s 97 s 136 s 177 s 3,0 bvis = 0,014; ht = 1,00 m 25 s 57 s 99 s 137 s 179 s 3,4 bvis = 0,014; ht = 3,00 m 26 s 58 s 100 s 137 s 179 s 3,2 bvis = 0,015; ht = 1,00 m 26 s 58 s 101 s 138 s 181 s 4,0 bvis = 0,02; ht = 3,00 m 28 s 61 s 106 s 146 s 191 s 7,6 Najboljše rezultate smo dobili pri simulaciji z velikimi delci ter vrednosti parametra bvis = 0,007 in višini lociranih točk nad terenom ht = 3 m. Pri teh pogojih je vrednost MAE = 1,6. Pri simulacijah z manjšimi delci smo najmanjšo vrednost MAE = 2,8 dobili pri vrednosti parametra bvis = 0,013 in višini lociranih točk nad terenom ht = 3 m. Pri teh dveh simulacijah smo preverili nihanja na sondah v dolini (slike 103 - 107) in iztok iz akumulacije (slika 102). Slika 102: Iztok iz akumulacije (hrapavo dno – podajanje točkovnih hrap). Figure 102: Outlet from the accumulation (rough terrain – point elements on the bottom of the slopes) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 101 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 103: Nihanje gladin v dolini na sondi 1(hrapavo dno – podajanje točkovnih hrap) Figure 103: Depth oscillations in valley at Gauge 1 (rough terrain – point elements on the bottom of the slopes) Slika 104: Nihanje gladin v dolini na sondi 2 (hrapavo dno – podajanje točkovnih hrap) Figure 104: Depth oscillations in valley at Gauge 2 (rough terrain – point elements on the bottom of the slopes) 102 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 105: Nihanje gladin v dolini na sondi 3 (hrapavo dno – podajanje točkovnih hrap) Figure 105: Depth oscillations in valley at Gauge 3 (rough terrain – point elements on the bottom of the slopes) Slika 106: Nihanje gladin v dolini na sondi 4 (hrapavo dno – podajanje točkovnih hrap) Figure 106: Depth oscillations in valley at Gauge 4 (rough terrain – point elements on the bottom of the slopes) Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 103 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Slika 107: Nihanje gladin v dolini na sondi 5 (hrapavo dno – podajanje točkovnih hrap) Figure 107: Depth oscillations in valley at Gauge 5 (rough terrain – point elements on the bottom of the slopes Tudi v tem primeru opazimo, da čelo vala v bližini sonde 2 potuje po bregovih in pozneje doseže dno doline, kjer se nahaja sonda. Globine na sondi 1 se dobro ujemajo z meritvami prvih 80 s, potem pa simulirane globine, v primerjavi z meritvami, padajo. Na sondi 3 dobimo nekoliko nižje globine, kot to kažejo meritve. Na ostalih sondah se simulirane globine dobro ujemajo z meritvami. Z uporabo metode SPH smo dobili vsaj tako dobre rezultate, kot je to pokazala uporaba konvencionalne metode. Odločitev o tem, katero numerično metodo bomo uporabili za izračun dogajanja v nekem primeru, je zelo pomembna, saj bo naša odločitev vplivala na točnost rezultatov. Zato je priporočljivo, da se za podobne izračune najprej izvrši natančno umerjanje in šele nato odloči, katera numerična metoda je primernejša za dani primer. Simulacije so bile izvedene na sistemu s štirjedrnim procesorjem Intel® CoreTM i7 3.33 GHz. V preglednici 12 so podani računski časi, ki smo jih potrebovali za 200 s dolge simulacije, izvedene z osnovnim in dopolnjenim modelom Tis Isat, ki omogoča upoštevanje oblikovnega upora toku vode s pomočjo dodajanja točkovnih elementov. 104 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Preglednica 12: Računski časi simulacij Table 12: Computational times Model Tis Isat Računski čas d = 5,0 m d = 2,5 m Osnovni model 1,6 h 54,0 h Dopolnjeni model 4,0 h 138,0 h Z 2x zmanjšanjem velikosti delcev, smo računski čas povečali za približno 34 x. Dopolnjeni model potrebuje približno 2,5 x več računskega časa kot osnovni model. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 105 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 6. ZAKLJUČKI IN NAPOTKI ZA NADALJNJE DELO 6.1 Zaključki in ugotovitve Poplavni valovi, ki nastanejo zaradi porušitev pregrad, običajno povzročijo veliko škodo. Gibanje vala pa nam opišejo osnovne nelinearne parcialne diferencialne osnovne enačbe, ki so rešljive samo s pomočjo uporabe numeričnih metod. Izbira numerične metode za reševanje osnovnih enačb je zelo pomembna. Nobena od že znanih metod pa ne zadovoljuje popolnoma in zato se razvijajo vedno novejše metode, ki omilijo ali odstranijo slabosti prejšnjih metod. Metoda SPH se je izkazala kot zelo uspešna metoda za simulacije hidrodinamičnih pojavov. Zelo dobro lahko simulira izrazito dinamične pojave s hipnimi spremembami vodne gladine s spreminjajočimi smermi in hitrostmi toka. Metoda je bila zato dobro sprejeta in v kratkem času preizkušena tudi na številnih primerih toka s prosto gladino. Za simulacije toka s prosto gladino v okviru doktorske disertacije smo uporabili osnovni matematični model Tis Isat, ki je namenjen simulacijam toka vode po metodi SPH. Najprej smo preverili rezultate modela Tis Isat z rezultati drugega SPH modela SPHysics ter jih primerjali z meritvami, pridobljenimi na fizičnem modelu. Študija je bila opravljena na primeru porušitve dvodimenzionalnega vodnega stolpca z različnimi razmerji med širino in višino vodnega stolpca. Model Tis Isat je dal vsaj tako dobre rezultate, kot smo jih dobili z modelom SPHysics, v daljšem časovnem obdobju po porušitvi pregrade pa znatno boljše rezultate. Pokazali smo, da se pri simulacijah z modelom SPHysics ponekod pojavijo motnje v gladini in fizikalno nelogični pojavi ob robovih. Omenjeni model definira steno s pomočjo navideznih delcev ob robu, model Tis Isat pa s pomočjo integracije med območjem stika vode z robom. Ta študija nam je še dodatno pokazala, kako pomembno je razvijati nove robne pogoje v metodi SPH. V nadaljevanju smo preverili delovanje modela Tis Isat na osnovnih testnih primerih, ki se običajno uporabljajo za preverjanje SPH modelov. To so simulacije toka v kanalu s hipno porušitvijo vodnega stolpca. Najprej smo izvajali simulacije v kanalu z 2D modelom Tis Isat. Meritve na fizičnem modelu so bile izvedene v laboratoriju Katedre za mehaniko tekočin z laboratorijem. V prvem primeru je voda prosto odtekala iz kanala, v drugem pa se je val odbil 106 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. od zapore. Na tem primeru smo preverili, v kolikšni meri se rezultati izboljšujejo z zmanjševanjem velikosti delcev. Pričakovano so se rezultati pri uporabi manjših delcev bistveno izboljšali. Nato smo preverili še delovanje 3D matematičnega modela Tis Isat v kanalu z oviro. Rezultate našega SPH modela smo primerjali z meritvami, rezultati drugih dveh SPH modelov ter rezultati konvencionalnega modela. Rezultati našega SPH modela so se izkazali kot primerljivi z rezultati drugih modelov. Slabost metode SPH je ta, da so simulacije časovno zelo zamudne. Natančnost simulacij se dokazano povečuje z zmanjševanjem velikosti delcev, s tem pa se računski čas simulacij bistveno poveča. Iz tega razloga se razvijajo različni pristopi, ki skrajšujejo računski čas simulacij metode SPH. V praksi se velikokrat srečamo s primeri prizmatičnih kanalov, kjer je tloris obravnavanega območja nespremenljive oblike. V tem primeru je smiselno uporabiti 2D SPH model, ki je časovno manj potraten. V kanalih z razširitvijo oz. zožitvijo, kjer se tloris kanala spreminja, pa pride v poštev samo 3D SPH model. Velikokrat pa se zgodi, da je kanal pred razširitvijo enostavne tlorisne geometrije. V tem primeru je na ravnem območju kanala smiselno uporabiti časovno hitrejši 2D, za razširitvijo pa zamuden 3D SPH model. V model Tis Isat smo za takšne primere vgradili vmesnik, ki omogoča povezovanje med modeli različne dimenzionalnosti. Delovanje tega vmesnika smo preverili na primeru toka v kanalu z razširitvijo ter v kanalu z zožitvijo in razširitvijo. Rezultate povezanega 2D/3D Tis Isat modela smo primerjali z rezultati 3D Tis Isat modela, rezultati konvencionalnega modela ter z meritvami. Z uporabo omenjenega vmesnika za povezovanje med modeli različne dimenzionalnosti praktično nismo vplivali na zmanjšanje točnosti rezultatov, smo pa bistveno skrajšali čas simulacij. V praksi se velikokrat srečamo s primeri, kjer nas bolj natančni rezultati toka zanimajo na samo enem delu obravnavanega območja. Ker je metoda SPH časovno izredno zamudna, bi bilo izvajanje simulacij z manjšimi delci na celotnem območju težko izvedljivo. Iz tega razloga smo v okviru doktorske disertacije razvili vmesnik, ki nam omogoča povezovanje med modeli različne resolucije. Preizkusili smo ga na primeru toka v kanalu s hipno razširitvijo. Rezultate povezanega 3D modela Tis Isat, ki omogoča povezovanje med modeli različne resolucije, smo primerjali z rezultati 3D modela Tis Isat, rezultati konvencionalnega modela in meritvami. Pokazali smo, da nam omenjeni vmesnik znatno skrajša računski čas simulacij, s tem da praktično ne zmanjša točnosti rezultatov. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 107 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Nazadnje smo uporabnost modela Tis Isat preverili še na praktičnem primeru toka iz narave. Simulirali smo poplavni val, ki bi nastal pri morebitni porušitvi zgornje akumulacije črpalne elektrarne Kolarjev vrh. Takšne elektrarne so ene izmed najbolj gospodarsko upravičenih in hkrati naravi najbolj prijaznih elektrarn. Na tem primeru smo preverili različne možnosti definiranja robnih pogojev. Teren smo definirali na dva načina. Najprej smo predpostavili, da je teren gladek, kjer smo hrapavost podali s pomočjo parametra viskoznosti ob stenah. Nato smo podali še realno hrapav teren, kjer smo hrapavost terena definirali bodisi s pomočjo dvigovanja mrežnih vozlišč ali pa s podajanjem točkovnih elementov na pobočjih. Rezultati simulacij so nam pokazali, kako zelo pomembno je pravilno definiranje hrapavosti terena, saj so se pri uporabi različnih pristopov pojavile razmeroma velike razlike v rezultatih. Najbolj točne rezultate so dale simulacije, kjer smo teren predpostavili kot realno hrapav z dvigovanjem mrežnih vozlišč. Na tem primeru smo analizirali vpliv parametra viskoznosti ob stenah in velikosti delcev na rezultate. S pomočjo parametra viskoznosti ob steni ( bvis) pa vplivamo na hitrost premikanja vala po terenu in sicer večja kot je vrednost parametra bvis, počasneje se val premika po terenu. Pokazali smo, da rezultati simulacij z manjšimi delci bolje sledijo meritvam kot rezultati simulacij z večjimi delci. Poudariti je potrebno, da so uporabljeni delci še vedno preveliki in da bi se z zmanjševanjem delcev šele lahko pokazala resnična prednost metode SPH pri opisu gladine. Z obravnavanim primerom, ki predstavlja doslej v svetu enega redkih primerov simuliranja porušitvenega vala na realni topografiji z metodo SPH, smo dokazali uporabnost te napredne metode tudi za praktične aplikacije. Tovrstne simulacije omogočajo točnejšo določitev poplavnih območij in s tem pripomorejo k zagotovitvi večje varnosti ljudi in njihove lastnine. Pokazalo se je, da je naš razviti in dopolnjeni model Tis Isat, ki osnovne enačbe toka s prosto gladino rešuje z numerično metodo SPH, uporaben za tovrstne simulacije. Dobili smo rezultate, ki so primerljivi z rezultati drugih SPH modelov in konvencionalnih modelov, ki uporabljajo numerične metode končnih volumnov, končnih elementov ali končnih razlik. Prednost metode SPH pa se je pokazala predvsem pri simulacijah pojavov, kjer pri toku s prosto gladino prihaja do izrazitih dinamičnih sprememb hitrosti in gladin, do česar pride npr. v kanalih z razširitvijo in zožitvijo. Razviti SPH model Tis Isat je zelo dobro simuliral tudi tok vode, ki nastane zaradi porušitve dela nasipa na realni topografiji dolvodne doline, kar je bil tudi eden osnovnih ciljev doktorske disertacije. 108 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Bistvene prednosti metode SPH za obravnavane primere toka lahko strnemo v naslednjih točkah: - metoda SPH je brezmrežna metoda in nima numerične mreže, ki bi lahko vplivala na točnost rezultatov. - metoda SPH je delčna metoda in zato tudi zelo primerna za simulacije izrazito dinamičnih pojavov, kot je gibanje porušitvenih valov. V okviru disertacije smo dokazali, da je metoda primerna za tovrstne simulacije. - metoda SPH računsko območje razdeli s pomočjo masnih delcev, ki se premikajo po prostoru. Takšna delčna metoda pa je zelo primerna za prikaz gibanja vode pri izrazitih dinamičnih spremembah gladine in omogoča izdelavo animacij toka vode. Manjši so uporabljeni delci, bolj nazorna je animacija. Metoda ima tudi nekatere slabosti: - podajanje robnih pogojev ni vedno urejeno na fizikalno logičen način. Običajno se rob opiše s pomočjo navideznih delcev ob robu, ki pa povzročajo nekatere fizikalno nelogične pojave (poglavje 4.1). - z manjšanjem delcev dokazano izboljšujemo rezultate, hkrati pa negativno vplivamo na računski čas simulacij. Metoda SPH je časovno zamudna, vendar pričakujemo, da se bo z nadaljnjim razvojem računalništva kmalu dalo simulirati pojave s še manjšimi delci v računsko sprejemljivem času. Nekateri raziskovalci metode SPH so že razvili paralelno kodo, s pomočjo katere lahko izvajajo bistveno zahtevnejše simulacije glede števila in velikosti delcev. - v metodo SPH so vgrajene nekatere empirične relacije (npr. hitrost zvoka je omejena na desetkratnik maksimalne predvidene hitrosti delcev tekočine), zato jo nekateri teoretiki popolnoma zavračajo 6.2 Napotki za nadaljnje delo Pri nadaljnjem razvoju numeričnega SPH modela Tis Isat je potrebno dati poseben poudarek razvoju novih pristopov za skrajšanje računskega časa simulacij in nadaljnjim izboljšavam robnih pogojev. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 109 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Model Tis Isat je že imel razvit vmesnik, ki je omogočal povezovanje med modeloma Tis Isat različne dimenzionalnosti. V okviru doktorske disertacije smo razvili vmesnik, ki omogočata povezovanje med Tis Isat modeloma različne resolucije. V nadaljevanju bi bilo smiselno razviti tudi vmesnik, ki bi nam omogočal povezovanje med SPH in konvencionalnimi modeli. V tem primeru bi lahko v primeru prizmatičnih kanalov na tlorisno nespremenljivem obravnavanem območju izvajali simulacije s hitrejšimi konvencionalnimi modeli ter na območjih, kjer bi prihajalo do hitrih sprememb globin in hitrosti vode, uporabili SPH model. Podatke o toku vode in premeščanju plavin moramo poznati npr. pri simulacijah delovanja hidroenergetskih objektov. Voda erodira zemljo in sedimenti zasipavajo akumulacijo, kar nam zmanjša razpoložljivi volumen vode v akumulaciji. Če hočemo spet povečati volumen, namenjen za vodo, moramo sedimente odstraniti z izpiranjem. Za simulacije takšnih primerov bi bilo smiselno model Tis Isat nadgradili z za vodo in sedimente prepustnimi robnimi pogoji, kar bi omogočilo simulacijo dotoka oziroma iztoka vode in plavin iz akumulacije. 110 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 7 POVZETEK Doktorska disertacija obravnava simulacije valov, ki bi nastali zaradi morebitnih porušitev pregrad po metodi hidrodinamike zglajenih delcev oziroma metodi Smoothed Particle Hydrodynamics (SPH). Takšni problemi se običajno simulirajo s pomočjo matematičnih modelov, ki osnovne enačbe rešujejo s konvencionalnimi numeričnimi metodami. Te že preverjene metode so zasnovane na Eulerjevem načinu diskretiziranja računskega območja, kjer se računsko območje definira s pomočjo numerične mreže, v kateri se izračunavajo fizikalne količine. Največja slabost konvencionalnih numeričnih metod je numerična difuzija, ki nam lahko bistveno poslabša točnost rezultatov. Numerični difuziji se lahko izognemo z uporabo Lagrangeove metode hidrodinamike zglajenih delcev. Lagrangeove metode diskretizirajo računsko območje s pomočjo samih delcev, ki se premikajo po prostoru. To je novejša metoda, ki je v zelo kratkem času postala priljubljena za simulacije toka s prosto gladino. Metoda je še zlasti primerna za simulacije izrazitih dinamičnih pojavov, kjer prihaja do hitrih sprememb gladine in hitrosti vode, kot so npr. lomi valov. V okviru disertacije smo uporabili osnovni matematični model Tis Isat, ki je bil razvit na Katedri za mehaniko tekočin z laboratorijem UL FGG in se uporablja za simulacije toka vode po metodi SPH. Model smo dopolnili in ga preverili na številnih primerih toka s prosto gladino. Najprej smo izvajali simulacije z 2D modelom Tis Isat. Rezultate našega modela smo primerjali z rezultati drugega SPH modela SPHysics ter z meritvami. Nato smo izvedli simulacije še z 3D modelom Tis Isat, ki je delovanje modela še dodatno potrdil. Znano je, da so simulacije po metodi SPH časovno zelo zamudne. V praksi pa se velikokrat srečamo s primeri, kjer nas boljša natančnost zanima samo na določenem delu računske domene. Za takšne primere ima model že izdelan vmesnik, ki nam omogoča povezovanje med modeli različne dimenzionalnosti. V okviru doktorata pa smo izdelali vmesnik, ki nam omogoča povezovanje med modeli različne resolucije. Oba vmesnika smo preverili na primerih toka vode v kanalih z razširitvijo in/ali zožitvijo. S takšnim pristopom smo bistveno znižali računski čas simulacij, nismo pa vplivali na točnost rezultatov. V osnovni model Tis Isat smo vgradili tudi nove robne pogoje. Upor toku vode je običajno hrapavostni in oblikovni. Hrapavostni upor se pri metodi SPH običajno določi s pomočjo Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 111 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. parametrov viskoznosti. Z dopolnitvijo modela, ki nam omogoča dodajanje točkovnih elementov na dnu pobočij, pa lahko opišemo oblikovni upor, ki ga predstavljajo razne zarasti na terenu. Dopolnjeni model smo preverili na realni topografiji. Simulirali smo poplavne valove, ki bi nastali zaradi morebitne porušitve dela nasipa akumulacije črpalne elektrarne Kolarjev vrh. Preverili smo obnašanje vode na gladkem in hrapavem terenu. Gladek teren smo definirali s pomočjo parametra viskoznosti ob stenah, ki je bil konstanten ali spremenljiv, odvisno od naklona terena. Hrapav teren pa smo podali na dva načina in sicer najprej s pomočjo dvigovanja mrežnih vozlišč, nato pa s podajanjem točkovnih elementov. Simulacije so nam pokazale razmeroma veliko razliko v točnosti rezultatov v odvisnosti od uporabljenega načina podajanja hrapavosti. S tem smo pokazali, da je pravilno definiranje hrapavosti terena ključnega pomena pri simulacijah toka porušitvenih valov na realni topografiji. Model Tis Isat se je izkazal kot zelo uporabno orodje za simulacijo porušitvenih valov po metodi SPH in bo pripomogel k povečanju točnosti tovrstnih računov v hidrotehnični praksi. 112 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. 8 SUMMARY This thesis presents numerical simulations of flood waves caused by dam failure using the meshless numerical method of smoothed particle hydrodynamics (SPH). Flood waves are usually simulated by mathematical models which mostly use conventional numerical techniques to solve the set of governing equations. These methods are based on the Eulerian approach that defines the computational domain by using a numerical mesh. Eulerian methods are most commonly used for fluid flow simulations and they give good results, but they also have disadvantages. One of the greatest weaknesses of mesh-based methods is numerical diffusion that can significantly decrease the quality of the results. Numerical diffusion can be completely eliminated by using the meshless Lagrangian SPH method. Lagrangian approach discretizes the computational domain into moving mass particles. To date, the new SPH approach has become increasingly popular for simulations in fluid mechanics. The SPH method is appropriate for simulating phenomena with rapidly changing velocities and free surface (for example: breaking waves). In this thesis, the basic mathematical model Tis Isat was used. This model was developed at the Chair of fluid mechanics with laboratory at the Faculty of civil and geodetic engineering of the University of Ljubljana and designed for simulating free surface flows based on the SPH method. The simulations of the upgraded Tis Isat model were compared to measurements of laboratory experimental testing. First, we verified the 2D Tis Isat model and simulated a collapse of a two-dimensional rectangular fluid column. Results of our simulations were compared to measurements and results of another SPH model SPHysics. Thereafter, the 3D Tis Isat model was tested. Both models gave satisfactory results. SPH method is very time-consuming. Often we come across cases where simulations with high resolution are required on certain parts of the computational domain. The 2D and 3D Tis Isat models can already be coupled in one- and two-way coupling. Within this thesis, we developed a new interface for coupling between models of different resolutions. Both interfaces were used to simulate fluid flows in channels with a contraction and/or an expansion. By using these two coupling approaches, the computational time was significantly reduced without decreasing the accuracy of results. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 113 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Drag is usually accounted by roughness and form. In the SPH method, roughness can be described by using a coefficient to calculate eddy viscosity on the walls. The Tis Isat model was upgraded with a new boundary condition that enables us to put point elements on selected locations. By using these elements, flow between obstacles can be simulated. The upgraded Tis Isat model was verified on a real case scenario. Flood waves due to a dam break at the upper storage reservoir of the pumped storage HPP Kolarjev vrh in Slovenia were simulated. Different ways of defining terrain roughness were used and flows over smooth and rough terrain were studied. The roughness of the smooth terrain was controlled using constant or space-variable eddy viscosity coefficients to calculate the dam viscosity between wall-particle interactions. Rough terrain was controlled by elevating mesh-nodes or loading the point elements. The simulations show the remarkably improved results by using rough terrain. The Tis Isat model has proved to be a useful tool for flood wave simulations using the SPH method and will help to improve the accuracy of flood wave simulations. 114 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. VIRI Børve, S., Omang, M., Trulsen, J. 2001. Regularized smoothed particle hydrodynamics: a new approach to simulating magnetohydrodynamic shocks. The Astrophysical Journal 561: 82–93. Børve, S., Omang, M., Trulsen, J. 2004. Regularized smoothed particle hydrodynamics with improved multi-resolution handling. Journal of Computational Physics 208: 345–367. Colagrossi, A., Landrini, M. 2003. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. Journal of Computational Physics 191: 448-475. Crespo, A.J.C., Gómez-Gesteira, M., Dalrymple, R.A. 2007. Boundary conditions generated by dynamic particles in SPH methods. Computers, Materials & Continua 5, 3: 173-184. Četina, M. 1988. Matematično modeliranje dvodimenzionalnih turbulentnih tokov. Magistrsko delo. Ljubljana, Fakulteta za gradbeništvo in geodezijo: 126 str. Četina, M. 1992. Tridimenzionalni matematični baroklini model za izračun totov v jezerih in morju. Doktorska disertacija. Ljubljana, Fakulteta za gradbeništvo in geodezijo: 127 str. Dalrymple, R.A., Knio, O. 2000. SPH modelling of water waves. Coastal Dynamics: 141-147. Džebo, E., Žagar, D., Četina, M., Jeromel, M., Bajcar, T., Širok, B., Petkovšek, G. 2011. Analiza porušitve pregrade po metodi SPH. V: Eberlinc, M. (ur.), Širok, B. (ur.). Kuhljevi dnevi, 22. september 2011, Mengeš. Zbornik del. Mengeš: SDM – Slovensko društvo za mehaniko: str. 25-32. Džebo, E., Žagar, D., Četina, M., Petkovšek, G. 2012 a. Simulation of dam-break flow in channel expansion with coupled 2-D/3-D SPH model. 7th International SPHERIC Workshop, 29. – 31. maj 2012, Prato. Zbornik del. Prato: SPHERIC – the SPH European Research Interest Community: str. 403-408. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 115 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Džebo, E., Žagar, D., Četina, M., Petkovšek, G. 2012 b. Simulacije valov zaradi porušitve pregrade v kanalu z razširitvijo s povezanimi modeli različne resolucije po metodi SPH. V: Hriberšek, M. (ur.), Ravnik, J. (ur.). Kuhljevi dnevi, 26. in 27. september 2012, Rogaška Slatina. Zbornik del. Rogaška Slatina: SDM – Slovensko društvo za mehaniko: str. 25-32. Džebo, E., Žagar, D., Četina, M., Petkovšek, G. 2013 a. Reducing the computational time of the SPH method with a coupled 2-D/3-D SPH model. Journal of Mechanical Engineering: sprejeto v objavo. Džebo, E., Žagar, D., Četina, M., Petkovšek, G. 2013 b. Simulacije vodne gladine v akumulaciji po metodi SPH. Kuhljevi dnevi, 25. in 26. september 2013, Rogaška Slatina. Zbornik del. Rogaška Slatina: SDM – Slovensko društvo za mehaniko: sprejeto v objavo. Falappi, S., Galatti, M., Maffio, A. 2007. SPH simulation sediment scour in reservoir sedimentation problems. 2nd International SPHERIC Workshop, 23. – 25. maj 2007, Madrid. Zbornik del. Madrid: SPHERIC – the SPH European Research Interest Community: str. 9-12. Feldman, J., Bonet, J. 2007. Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems. International Journal for Numerical Methods in Engineering 72: 295–324. Fulk, D.A. 1994. A numerical analysis of smoothed particle hydrodynamics. Doktorska disertacija. Hobson Way, Air Force Institute of Technology: 352 str. Gatti, D., Maffio, A., Zucala, D., Di Monaco, A. 2007. SPH simulation of hydrodynamics problems related to dam safety. 32nd Congress of IAHR, 1. – 6. julij 2007, Venice. Zbornik del. Venice: IAHR – International Association for Hydro – Environment Engineering and Research. Gingold, R.A., Monaghan, J.J. 1977. Smoothed particle hydrodynamics. Theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society 181: 375-389. 116 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Gómez-Gesteira, M., Dalrymple, R.A. 2004. Using a Three-Dimensional Smoothed Particle Hydrodynamics Method for Wave Impact on a Tall Structure. Journal of Waterway, Port, Coastal and Ocean Engineering 130: 63-69. Gómez-Gesteira, M., Rogers, B.D., Dalrymple, R.A., Crespo, A.J.C., Narayanaswamy, M. 2009. User Guide for the SPHysics code: 73 str. Gómez-Gesteira, M., Rogers, B.D., Dalrymple, R.A., Crespo, A.J.C. 2010. State-of-the-art of classical SPH for free-surface flows. Journal of Hydraulic Research 48: 6-27. Gotoh, H., Shao S., Memita, T. 2004. SPH-LES model for numerical investigation of wave interaction with partially immersed breakwater. Coastal Engineering Journal 46: 39-63. Issa, R., Violeau, D. 2006. 3D dambreaking (Test-case 2). SPH European Research Interest Community: 8 str. Jones, D.A., Belton, D. 2006. Smoothed particle hydrodynamics: applications within DSTO. Victoria, Department of Defence: 47 str. Kleefsman, K.M.T., Fekken, G., Veldman, A.E.P., Iwanowski, B., Buchner, B. 2005. A Volume-of-Fluid based simulation method for wave impact problems. Journal of Computational Physics 206: 363–393. Koshizuka, S., Oka, Y. 1996. Moving Particle Semi-Implicit Method for Fragmentation of Incompressible Fluid. Nuclear Science and Engineering 123: 421-434. Krzyk, M. 2004. Dvodimenzijsko matematično modeliranje toka v strmih strugah. Doktorska disertacija. Ljubljana, Fakulteta za gradbeništvo in geodezijo: 126 str. Kulasegaram, S., Bonet, J., Lewis, R.W., Profit, M. 2003. A variational formulation based contact algorithm for rigid boundaries in two-dimensional SPH applications. Computer Mechanics 33, 4: 316-325. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 117 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Lastiwka, M., Quinland, N., Basa, M. 2005. Adaptive particle distribution for smoothed particle hydrodynamics. International Journal for Numerical Methods in Fluids 47: 1403– 1409. Lee, E-S., Violeau, D., Issa, R., Ploix, S., Marc, R. 2009. Simulating a Real Dam Spillway Flow with 3-D SPH. 4th International SPHERIC Workshop, 26. – 29. maj 2009, Nantes. Zbornik del. Nantes: SPHERIC – the SPH European Research Interest Community: str. 339- 345. Lee, H., Han, S. 2010. Solwing the shallow water equations using 2D SPH partlices for interactive applications. The Visual Computer 26: 865–872. Legiša, D., Rajar, R. 1980. Modelna raziskava vala ob hipni porušitvi nasipa akumulacije Kolarjev vrh (ČE Kozjak). Strokovna naloga. Ljubljana, Vodnogospodarski inštitut in Univerza v Ljubljani, Fakulteta za arhitekturo, gradbeništvo in geodezijo, Laboratorij za mehaniko tekočin: 27 str. Libersky, L.D., Petsheck, A.G., Carney, T.C., Hipp, J.R., Allahdadi, F.A. 1993. High strain Lagrangian hydrodynamics – a three-dimensional SPH code for dynamic material response. Journal of Computational Physics 109: 67-75. Liu, M.B., Liu, G.R., Zong, Z., Lam, K.Y. 2003 a. Computer simulation of the high explosive explosion using smoothed particle hydrodynamics methodology. Computers & Fluids 32, 3: 305-322. Liu, G.R., Liu, M.B. 2003 b. Smoothed Particle Hydrodynamics - a mesh free particle method, World Scientific: 449 str. Lopez, D., Marivela, R., Arranda, F. 2009. SPH model calibration using data from pressure of the prototype stilling basin of Villar del Rey dam, Spain. 33rd Congress of IAHR, 9. – 14. avgust 2009, Vancouver. Zbornik del. Vancouver: IAHR – International Association for Hydro – Environment Engineering and Research: str. 187-198. 118 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Lopez, D., Marivela, R., Garrote, L. 2010. Smoothed particle hydrodynamics model applied to hydraulic structures: a hydraulic jump test case. Journal of Hydraulic Research 48: 142- 158. Lucy, L.B. 1977. A numerical approach to the testing of the fission hypothesis. The Astronomical Journal 82, 12: 1013-1024. Martin, J.C., Moyce, W.J. 1952. An experimental study of the collapse of liquid columns on a rigid horizontal plane. Philosophical Transactions of the Royal Society 244: 312 – 324. Molteni, D., Colagrossi, A. 2009. A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH. Computer Physics Communications 180, 6: 861-872. Morris, J.P. 1996. Analysis of Smoothed Particle Hydrodynamics with Applications. Doktorska disertacija. Melbourne, Monash University: 338 str. Morris, J.P., Fox, P.J., Zhu, Y. 1997. Modeling Low Reynolds Number Incompressible Flows Using SPH. Journal of Computational Physics 136: 214-226. Monaghan, J.J., Lattanzio, J.C. 1985. A refined particle method for astrophysical problems. Astronomy and Astrophysics 149: 135-143. Monaghan, J.J. 1992. Smooth Particle Hydrodynamics. Annual Review of Astronomy and Astrophysics 30: 573-574. Monaghan, J.J. 1994. Simulating Free Surface Flows with SPH. Journal of Computational Physics 110: 399-406. Monaghan, J.J., Kos, A. 1999. Solitary Waves on a Cretan Beach. Journal of Waterway, Port, Coastal and Ocean Engineering 125, 3: 145-154. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 119 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Narayanaswamy, M., Crespo, A.J.C., Gómez-Gesteira, M., Dalrymple, R.A. 2010. SPHysics - FUNWAVE hybrid model for coastal wave propagation. Journal of Hydraulic Research 48: 85–93. Patankar, S.V. 1980. Numerical heat transfer and fluid flow. McGraw-Hill book company: 197. Petkovšek, G., Džebo, E., Četina, M., Žagar, D. 2010 a. Application of Non-Discrete Boundaries with Friction to Smoothed Particle Hydrodynamics. Journal of Mechanical Engineering 56, 5: 307 – 315. Petkovšek, G., Džebo, E., Četina, M., Žagar, D. 2010 b. Simulacija laboratorijske porušitve pregrade po metodi SPH z analizo modela trenja. V: Širok, B. (ur.), Eberlinc, M. (ur.). Kuhljevi dnevi, 23. september 2010, Ljubljana. Zbornik del. Ljubljana: SDM – Slovensko društvo za mehaniko: str. 159–166. Popovska, C. 1988. Numericko i eksperimentalno istražuvanje na dvodimenzionalno nestacionarno tecenje vo otvoreni tokovi. Doktorska disertacija. Skopje, Fakulteta za gradbeništvo. Prakash, M., Debroux, F., Gleary, P. 2001. Three-dimensional modelling of dam-break induced flows using Smoothed Particle Hydrodynamics. 14th Australasian fluid mechanics conference, 10. – 14. december 2001, Adelaide. Zbornik del. Adelaide: AFMS – Australasian Fluid Mechanics Society: str. 379-382. Rajar, R. 1972. Recherche théorique et expérimentale sur la propagation des ondes de rupture de barrage dans une vallée naturelle. Doktorska disertacija. Toulouse, Univerza Paul Sabatier: 342 str. Rajar, R. 1978. Mathematical Simulation of Dam-Break Flow. Journal of the Hydraulics Division 104: 1011-1026. Rickenmann, D. 1996. Fliessgeschwindigkeit in Wildbächen und Gebrigsflüssen. Wasser, Energie, Luft – eau, énergie, air 88, 11-12: 298-304. 120 Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….…hidrodinamike zglajenih delcev Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Rickenmann, D. 2000. Dynamics of sediments and water in alpine catchments – processes and prediction. Scientific report WSL, Birmensdorf. Rogers, B.D., Darlymple, R.A., Stansby, P.K. 2010. Simulation of caisson breakwater movement with 2-D SPH. Journal of Hydraulic Research 48: 6-27. Roubtsova, V., Kahawita, R. 2006. The SPH technique applied to free surface flows. Computers and Fluids 35: 1359-1371. Shao, S., Lo, E.Y.M. 2003. Incompressible SPH method for simulating Newtonian and non- Newtonian flows with a free surface. Advenaces in Water Resources 26: 787-800. Sibilla, S. 2007. SPH simulation of local scour processes. 2nd International SPHERIC Workshop, 23. – 25. maj 2007, Madrid. Zbornik del. Madrid: SPHERIC – the SPH European Research Interest Community: str. 169-172. Smagorinsky, J. 1963. General circulation experiments with the primitive equations, I. the basic experiment. Monthly Weather Review 91: 99-164. Tropea, C., Yarin, A.L., Foss, J.F. 2007. Springer Handbook of Experimental Fluid Mechanics. Berlin, Springer-Verlag: 1557 str. Vacondio, R., Rogers, B.D., Stansby, P.K., Mignosa, P., Feldman, J. 2013. Variable resolution for SPH: A dynamic particle coalescing and splitting scheme. Computer Methods in Applied Mechanics and Engineering 256: 132-148. Violeau, D., Issa, R. 2007. Influence of turbulence closure on free surface flow modelling with SPH. 32nd Congress of IAHR, 1. – 6. julij 2007, Venice. Zbornik del. Venice: IAHR – International Association for Hydro – Environment Engineering and Research. Violeau, D. 2012. Fluid Mechanics and the SPH Method – Theory and Applications, Oxford University Press: 594 str. Džebo, E. 2013. Simulacije valov zaradi porušitev pregrad z ……….… hidrodinamike zglajenih delcev 121 Doktorska disertacija. Ljubljana, UL FGG, Doktorski študijski program Grajeno okolje, smer Gradbeništvo. Wendland, H. 1995. Piecewiese polynomial, positive definite and compactly supported radial functions of minimal degree. Advances in computational Mathematics 4, 1: 389—396. Žagar, D., Džebo, E., Četina, M., Petkovšek, G. 2008. Effects of boundary friction in SPH flow simulations. SPHERIC newsletter 7: str. 3. Žagar, D., Džebo, E., Četina, M., Petkovšek, G. 2009. Simulations of dam break and flow through a steep valley using SPH. 33rd Congress of IAHR, 9. – 14. avgust 2009, Vancouver. Zbornik del. Vancouver: IAHR – International Association for Hydro – Environment Engineering and Research: str. 171-178.