Građevinski fakultet u Subotici Univerzitet u Novom Sadu
Miodrag Spasojević
NUMERIČKA HIDRAULIKA - OTVORENI TOKOVI -
1996.
NUMERIČKA HIDRAULIKA -OTVORENI TOKOVIMiodrag Spasojević
RECENZENTI: Prof. Dr. Marko Ivetić Prof. Dr. Miodrag Jovanović
Subotica, 1996
SADRŽAJ UVOD 1.
PREGLED NUMERIČKIH METODA 1.1. Metod konačnih elemenata 1.2. Metod karakteristika 1.3. Metod konačnih razlika (priraštaja) 1.4. Konsistencija, konvergencija, stabilnost
1 6 10 19 24
2.
JEDNODIMENZIONALNI NEUSTALJENI TRANSPORT ZAGAĐIVAČA 2.1. Jednačina održanja mase zagađivača (jednačina transporta zagađivača) 2.2. Jednodimenzionalna advekcija 2.2.1. Metod konačnih razlika 2.2.2. Metod karakteristika 2.3. Jednodimenzionalna difuzija 2.3.1. Metod konačnih razlika 2.4. Jednodimenzionalna advekcija i difuzija 2.4.1. Metod razlomljenih koraka (metod razdvojenih operatora) 2.4.2. Hibridni metod
37
3.
JEDNODIMENZIONALNO NEUSTALJENO TEČENJE U OTVORENIM TOKOVIMA 3.1. Integralni oblik jednačina 3.2. Diferencijalni oblik jednačina 3.3. Granice važenja jednodimenzionalnih jednačina. Složeni preseci 3.4. Kratak pregled šema metoda konačnih razlika 3.5. Preissmann-ova šema 3.5.1. Primena Preissmann-ove šeme 3.5.2. Granični uslovi
ZADACI
37 47 47 54 69 69 77 78 81
88 88 99 104 109 116 125 129
UVOD
Ovaj materijal nastao je kao posledica predavanja predmeta Numerička hidraulika tokom poslediplomske nastave i dela (otprilike polovine) predmeta Hidraulika 2 tokom redovne nastave. Cilj pomenutih predmeta je da se stekne teorijsko znanje i praktično iskustvo u pravljenju i korišćenju matematičkih modela u inženjerskoj hidraulici (uglavnom) otvorenih tokova. Odnosno, cilj je da se što bolje izuči veza između razmatranih fizičkih pojava, jednačina kojima se te pojave opisuju i numeričkih metoda kojima se jednačine rešavaju, kako bi se rezultati proračuna što bolje tumačili i analizirali. Ovaj pisani materijal podeljen je na tri poglavlja. Poglavlje 1. predstavlja uvod u numeričke metode za rešavanje parcijalnih diferencijalnih jednačina kojima se opisuju fizički procesi koji se izučavaju u predmetima Hidraulika i Mehanika fluida. U poglavlju 2. izvedena je jednačina održanja mase zagađivača za ciljem da se razjasne mehanizmi transporta zagađivača u fluidnoj struji. Nadalje su u poglavlju 2. razmatrani numerički metodi za rešavanje jednačine kojom se opisuje jednodimenzionalni neustaljeni tranport zagađivača. U poglavlju 3. izvedene su jednačine kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima (tzv. de St. Venant-ove jednačine) i proučavani numerički metodi za njihovo rešavanje. Izučavanje transporta zagađivača ovde ima specifičnu ulogu. Mada je to važna (i trenutno veoma popularna) oblast sama za sebe, ovde je modelisanje transporta zagađivača korišćeno i kao značajno pedagoško sredstvo. Naime, transport zagađivača se opisuje jednom parcijalnom diferencijalnom jednačinom, koja je linearna, a uz to se može podeliti na dva dela od kojih oba imaju jasan matematički karakter, pa stoga i jasno definisane potrebe za početnim i graničnim uslovima. Smatralo se da je pogodnije da se principi numeričkog rešavanja prvo razjasne na primeru jednačine transporta zagađivača, pre nego što se pređe na razmatranje strujanja vode, koje se opisuje sistemom nelinearnih parcijalnih diferencijalnih jednačina. Zato se ovde transport zagađivača izučava u poglavlju 2, uz pretpostavku da je strujno polje poznato, dok se samo strujno polje odnosno rešavanje jednačina kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima izučava naknadno u poglavlju 3. Posebna pažnja posvećena je sticanju praktičnog iskustva u pravljenju i korišćenju matematičkih modela transporta zagađivača i jednodimenzionalnog neustaljenog tečenja u otvorenim tokovima. Zadaci, pridodati na kraju ovog materijala, omogućavaju da se delovi izložene materije savladaju i kroz praktičan rad. Zadaci su sastavljeni tako da prate redosled izložene materije, kao i da se rezultati prethodnih zadataka mogu koristiti u narednim zadacima.
Praktični problemi na čije rešavanje se ovde izučavana materija može primeniti obuhvataju: 1) Neustaljeno tečenje u otvorenim tokovima, kao što su: - reke (prirodni vodotoci), kod kojih neustaljenost potiče od poplavnih talasa (proleće/jesen), rada hidroelektrana itd., - akumulacije, kod kojih neustaljenost takođe potiče od poplavnih talasa ali i od upravljanja akumulacijom, - kišna kanalizacija, gde neustaljenost potiče od padavina odnosno oticanja, - kanalski sistemi za navodnjavanje, gde neustaljenost izaziva rad crpki, upravljanje sistemom (pomoću ustava i ostale kontrolne opreme), - kanali u sastavu hidroenergetskih postrojenja, gde je neustaljenost obično izazvana radom postrojenja, - ušća reka u more, gde neustaljenost potiče od od plime i oseke, - itd. 2) Transport zagađivača u otvorenim tokovima ili sistemima pod pritiskom, ali i transport suspendovanog nanosa, mešanje slatke i slane vode na ušću reka u more, itd., jer se sve ove pojave opisuju istom transportnom jednačinom. Za bolje razumevanje i dalje izučavanje izložene materije preporučuje se izbor literature, koja je delimično korišćena i prilikom sastavljanja ovog materijala: Cunge, Holly, Verwey (1980), Practical Aspcets of Computational River Hydraulics, Abbott (1979), Computational Hydraulics, Lai (1986), Numerical Modeling of Unsteady Open Channel Flow, Richtmyer, Morton (1967), Difference Methods for Initial Value Problems, Mahmood, Yevjevich (1975), Unstedy Flow in Open Channels, Carnahan, Luther, Wilkes (1969), Appliad Numerical Methods.
1. PREGLED NUMERIČKIH METODA Neustaljeno tečenje u otvorenim tokovima i transport zagađivača opisuju se parcijalnim diferencijalnim jednačinama, pa je korisno podsetiti se klasifikacije parcijalnih diferencijalnih jednačina, jer različiti tipovi jednačina zahtevaju različite početne i granične uslove. Takođe, pokazalo se da postoji jaka veza između matematičkog karaktera (tipa) jednačine i njenog numeričkog rešenja. Kod nekih tipova jednačina se relativno lako dolazi do pouzdanog numeričkog rešenja, dok je kod drugih to skopčano sa mnogim problemima. Stoga nije svaka numerička metoda podjednako pogodna za svaki tip jednačina. Kod onih tipova jednačina koji prouzrokuju probleme pri numeričkom rešavanju mora se mnogo više voditi računa o izboru metode, a najbolji rezultati se dobijaju ako je numerička metoda usklađena sa matematičkim karakterom (tipom) jednačine.
Klasifikacija parcijalnih diferencijalnih jednačina Opšta parcijalna diferencijalna jednačina drugog reda može se napisati u sledećem obliku: ∂ 2u ∂ 2u ∂ 2u ∂ u ∂ u , a +c 2 + 2b 2 = f x , y , u, ∂ x∂ y ∂ x ∂ y ∂x ∂y gde su: x, y - nezavisno promenljive, a u(x, y) - zavisno promenljiva. Napisana jednačina je: − linearna, ako su koeficijenti a, b, c = fx, y ∂ u ∂ u − kvazilinearna, ako su koeficijenti = f x , y , u, , , ∂ x ∂ y ∂ u ∂ u ∂ 2u ∂ 2u − nelinearna, ako su koeficijenti = f x , y , u, , , , . ∂ x ∂ y ∂ x2 ∂ y2 Napisana parcijalna diferencijalna jednačina je: − hiperboličnog tipa, ako je diskriminanta b2 - ac > 0, − paraboličnog tipa, ako je diskriminanta b2 - ac = 0, − eliptičnog tipa, ako je diskriminanta b2 - ac < 0. U nastavku će se, koristeći kao primere jednačine kojima se opisuju poznati fizički procesi, ukratko navesti osnovne odlike sva tri tipa jednačina.
(
Hiperbolične jednačine b2 − ac > 0
)
Kao primer navešće se jednačina kojom se opisuje jednodimenzionalna advekcija zagadjivača (transport zagađivača brzinama fluida): 1
∂ C ∂ C +u =0 ∂ t ∂ x gde su: x , t - nezavisno promenljive (koordinate po prostoru i vremenu), u(x, t) - parametar (brzina fluida poznata iz proračuna strujanja), C(x,t) - zavisno promenljiva (nepoznata koncentracija zagadjivača). Uočava se da je navedena jednačina linearna. Ako se zbog jednostavnosti pretpostavi da je brzina u konstantna, pa se prethodna jednačina diferencira po vremenu, dobija se:
∂ 2C ∂ 2C + =0 u ∂ t2 ∂x ∂ t odakle se, poređenjem sa prethodno napisanim opštim oblikom parcijalne diferencijalne u u2 jednačine, zaključuje da je: a = 1, b = , c = 0, pa je diskriminanta b 2 − ac = − 0 uvek 2 2 pozitivna, odnosno jednačina je hiperbolična. Kao drugi primer navešće se sistem pojednostavljenih (linearizovanih) jednačina kojima se opisuje jednodimenzionalno tečenje u kanalima (propagacija talasa): jednačina kontinuiteta:
∂ h ∂u + h0 =0 ∂t ∂x
dinamička jednačina:
∂ u ∂h +g =0 ∂t ∂x
gde su: x , t - nezavisno promenljive (koordinate po prostoru i vremenu), u(x,t), h(x,t) zavisno promenljive (brzina i dubina vode), h0, g - parametri (konstantna dubina i gravitaciono ubrzanje). Ako se prva od prethodne dve jednačine diferencira po vremenu, a druga po prostoru, pa se iz tako dobijene dve jednačine eliminiše mešoviti izvod, sledi:
gh0
∂ 2h ∂ 2h − =0 ∂ x2 ∂ t 2
odakle se, poređenjem sa prethodno napisanim opštim oblikom parcijalne diferencijalne jednačine, zaključuje da je: a = gh0, b = 0, c = -1, pa je diskriminanta b 2 − ac = 0 + gh0 uvek pozitivna, odnosno jednačina je hiperbolična. Kao što se iz navedenih primera vidi, hiperbolične jednačine su vezane za neustaljene procese kao što su propagacija talasa, advekcija zagađivača itd. Pošto se jednačinama obično opisuju neustaljene pojave, rešenje napreduje kroz vreme. Hiperbolične jednačine zahtevaju početni uslov i granične uslove na “uzvodnim” granicama (na onim granicama odakle nailazi talas).
2
(
Parabolične jednačine b2 − ac = 0
)
Kao primer navešće se jednodimenzionalni transport zagađivača difuzijom:
∂ C ∂ 2C =D ∂ t ∂ x2 gde su: x , t - nezavisno promenljive (koordinate po prostoru i vremenu), D(x, t) - parametar (koeficijent difuzije, D > 0), C(x,t) - zavisno promenljiva (nepoznata koncentracija zagadjivača). Uočava se da je navedena jednačina linearna. Poređenjem sa prethodno napisanim opštim oblikom parcijalne diferencijalne jednačine, zaključuje da je: a = D, b = c = 0, pa je diskriminanta b 2 − ac identički jednaka nuli, odnosno jednačina je parabolična. Parabolične jednačine su vezane za neustaljene pojave kao što su neustaljeni transport zagađivača difuzijom, neustaljeno provođenje toplote itd. Pošto se jednačinama obično opisuju neustaljene pojave, onda rešenje napreduje kroz vreme. Pararbolične jednačine zahtevaju početni uslov i granične uslove na svim granicama.
(
Eliptične jednačine b2 − ac < 0
)
Kao primer će se navesti ustaljeno dvodimenzionalno provođenje toplote (koje se opisuje poznatom Laplace-ovom jednačinom):
∂ 2T ∂ 2T + =0 ∂ x2 ∂ y2 gde su: x , y - nezavisno promenljive (koordinate po prostoru), T(x,y) - zavisno promenljiva (nepoznata temperatura). Poređenjem sa prethodno napisanim opštim oblikom parcijalne diferencijalne jednačine, zaključuje da je: a = c = 1, b = 0, pa je diskriminanta b 2 − ac = 0 − 1 uvek negativna, odnosno jednačina je eliptična. Eliptične jednačine su vezane za pojave kao što su ustaljeno dvodimenzionalno provođenje toplote, strujanje podzemnih voda itd. Eliptične jednačine zahtevaju samo granične uslove, i to na svim granicama.
* * *
3
Može se lako videti da su u prethodnim slučajevima, gde je matematički katakter (tip) jednačina bilo moguće lako i jednoznačno odrediti, korišćene pojednostavljene jednačine odnosno delovi kompletnih jednačina. U većini slučajeva kompletne parcijalne diferencijalne jednačine, kojima se opisuju pojave u hidraulici, se ne mogu jednostavno klasifikovati. Na primer, transport zagađivača uključuje i advekciju i difuziju zagađivača. Kompletna jednodimenzionalna jednačina transporta zagađivača glasi:
∂ C ∂C ∂ 2C +u =D ∂ t ∂x ∂ x2 Teoretski, navedena jednačina je parabolična, jer je diskriminanta b 2 − ac = 0 . Međutim, u slučaju tzv. “slabe” difuzije (ako D →0 ) tj. dominantne advekcije, jednačina postaje hiperbolična ( b 2 − ac 〉 0) . Praktično, navedena jednačina ima i parabolične i hiperbolične karakteristike, što može da stvara probleme pri numeričkom rešavanju: advektivni deo zahteva granične uslove samo na “uzvodnoj” granici, dok difuzioni deo zahteva granične uslove na svim granicama; neke numeričke metode su pogodnije za hiperbolične jednačine, dok su druge pogodnije za parabolične jednačine itd.
Postojeći numerički metodi Najpoznatiji (najčešće korišćeni) numerički metodi u hidraulici odnosno mehanici fluida su: − metod konačnih razlika (ili konačnih priraštaja), − metod karakteristika, − metod konačnih elemenata, i − metod graničnih elemenata. Pobrojaće se nekolko osnovnih karakteristika svakog od navedenih metoda, ali će se ovde nadalje uglavnom izučavati metod konačnih razlika i metod karakteristika. Metod konačnih razlika je najčešće korišćen u numeričkoj hidraulici i mehanici fluida, zbog svoje jednostavnosti i jasne veze između numerike i fizičkog značenja članova u jednačinama. Metod karakteristika se uglavnom koristi za izrazito hiperboličke jednačine - kojima se opisuju advekcija, propagacija strmih talasa itd.- jer najbolje odgovara karakteru jednačina. Ovakvi fizički procesi, odnosno ovaj tip jednačina je zapravo najkomplikovaniji za numeričko rešavanje. Metod konačnih elemenata je “inženjerski” metodi (prvo uveden od strane inženjera, a zatim proučavan od strane matematičara), koji se prvo pojavio, a i danas se najčešće koristi, u teoriji konstrukcija. Osnovna prednost ovog metoda je što, u odnosu na ostale metode, omogućava bolji opis oblasti strujanja. Stoga se ovaj metod u hidraulici i mehanici fluida koristi za modelisanje dvodimenzionalnog i trodimenzionalnog tečenja (lako je zamisliti kako se i najkomplikovanija npr. dvodimenzionalna oblast strujanja moze lako izdeliti na niz trougaonih elemenata). U slučaju jednodimenzionalnog tečenja metod konačnih elemenata 4
nema nikakvih posebnih prednosti u odnosu na, recimo, metod konačnih razlika. Iako će se ovde isključivo izučavati jednodimenzionalno tečenje, ipak će se napraviti kratak osvrt i na metod konačnih elemenata, tek toliko da se vidi koliko se ovaj metod razlikuje od npr. metoda konačnih razlika. Metod graničnih elemenata se uglavnom koristi za dvodimenzionalno i trodimenzionalno modelisanje strujanja podzemne vode. Ovo je integralni metod, koji koristi Green-ove funkcije da svede integraciju po trodimenzionalnoj oblasti strujanja na integraciju po (dvodimenzionalnoj) površini kojom je oblast strujanja ograničena, odnosno da svede integraciju po dvodimenzionalnoj oblasti strujanja na integraciju po (jednodimenzionalnoj) konturi oblasti. Pošto nije poznata primena ovog metoda u drugim oblastima osim strujanja podzemnih voda, to se ovaj metod ovde neće izučavati.
5
1.1. Metod konačnih elemenata Za razliku od diferencijalnih metoda kao što je metod konačnih razlika, koji se zasniva na približnom predstavljanju (aproksimiranju) izvoda u parcijalnim diferencijalnim jednačinama, metod konačnih elemenata je integralni metod. Postoji više načina da se objasni metod konačnih elemenata. Čini se da je najjednostavnije da se metod konačnih elemenata opiše kao metod otežanih ostataka (weighted-residual method) odnosno metod zaostale greške (residual error). Za to će se koristiti jednostavni primer. Razmatraće se jednodimenzionalna konzola (kao na skici), čiji se slobodni krajzagreva, tako da je uspostavljen konstantan (ustaljen) rasporedtemperature duž konzole. Toplota se prenosi provođenjem. Traži se raspored temperatura duž konzole.
Raspored temperature opisan je običnom diferencijalnom jednačinom: d 2T =3 dx 2 Zadati su sledeći granični uslovi:na slobodnom kraju (u x = 0), koji se konstantno zagreva, održava se konstantna temperatura: T = 1 ; na uklještenom kraju (u x = 3), koji se nalazi na zidu sa termičkom izolacijom, poznat je izvod temperature:
∂T ∂x
= 0.
Tačno (analitičko) rešenje:
T=
3 2 x − 9x + 1 2
se dobija jednostavno, tako što se napisana diferencijalnajednačina dvaput integrali, pa se integracione konstante odrede na osnovu graničnih uslova. Nadalje će se odrediti numeričko rešenje pomoću metoda konačnih elemenata. Konzola će se podeliti na elemente. U ovom primeru koristiće se svega tri elementa, zbog jednostavnosti. Pošto je problem jednodimenzionalan, to su i elementi jednodimenzionalni (u slučaju npr.dvodimenzionalnog problema elementi mogu biti trougaoni, pravougaoni itd.). Elementi se dodiruju u tzv. čvorovima ili računskim čvorovima (u ovom slučaju, sa tri elementa, postoje četiri čvora). Računski čvorovi ili računske tačke definišu tzv. računsku mrežu.
6
Pretpostaviće se da se približno rešenje T * može predstaviti kao linearna kombinacija vrednosti temperature T u čvorovima pomnoženih tzv. baznim funkcijama ϕ:
T ∗=
m
∑T ϕ j =1
j
j
gde je: m - broj čvorova; Tj - vrednost nepoznate temperature u čvoru j; ϕj - bazna funkcija za čvor j. Bazne funkcije su lokalne funkcije, različite od nule u okolini određenog čvora, a jednake nuli u svim ostalim čvorovima. Pod uslovom da zadrže lokalni karakter, bazne funkcije se mogu izabrati tako da budu linearne, kvadratne itd. Bazne funkcije takođe moraju zadovoljiti “osnovni” granični uslov T1=1.Na skici su prikazane izabrane funkcije za razmatrani primer. Uočava se da oblik približnog rešenja T * između čvorova zavisi od izbora ϕ. Zadatkom se traže nepoznate vrednosti Tj u čvorovima.
Definisaće se tzv. zaostala greška (residual error), tako što će se uvrstiti približno rešenje T * u polaznu jednačinu:
d 2T ∗ R( x) = −3≠ 0 dx 2 Ne može se zahtevati da zaostala greška bude jednaka nuli u svakom x (to bi bilo tačno rašenje), ali se može zahtevati da suma (t.j. integral) svih lokalnih zaostalih grešaka bude 7
jednak nuli (t.j. zahtevaće se da zaostala greška bude jednaka nuli, u osrednjenom ili “otežanom” smislu): 3
∫ R( x) W ( x) dx = 0 i
i = 1, 2, ...., m (*)
0
gde su Wi(x) tzv. težinske ili “test” funkcije. Nadalje će se koristiti najjednostavniji pristup, tzv. Galerkin-ov metod, koji za težinske funkcije W koristi bazne funkcije ϕ : ϕj =Wj Jednačina (*) postaje: 3
∫ 0
d2 m 2 ∑ T jϕ j ( x ) − 3ϕ i ( x )dx = 0 i = 1, 2, ..., m dx j =1
Na ovaj način je dobijen sistem od m jednačina - po jedna jednačina za svako i (t.j. ϕi) - sa m nepoznatih T (po jedno T u svakom čvoru). Preostaje da se izaberu bazne funkcije ϕ (koje u ovom slučaju služe i kao težinske funkcije). Zbog integraljenja bi bilo pogodno da funkcije ϕ budu što je moguće jednostavnije, d 2ϕ , što znači da bi što znači linearne funkcije. Međutim, u poslednjoj jednačini se javlja dx 2 funkcija ϕ morala biti barem kvadratna funkcija (zbog diferenciranja). Ovo poslednje se može izbeći ako se koristi parcijalna integracija (što je još jedan ključni element metode konačnih elemenata). Ako se se poslednja jednačina preuredi uvođenjem
d 2ϕ j
3 m
∫ ∑T 0 j =1
j
dx 2
d2 pod Σ: dx 2
3
ϕ i dx − ∫ 3ϕ i dx = 0
i = 1, ... m
0
pa zatim primeni parcijalna integracija: b
b
a
a
∫ udv =[ uv] ba − ∫ vdu
( d ( uv) = udv + vdu)
na prvi integral u jednačini, tako što se izabere da je: u = ϕi du =
dv =
dϕ i dx dx
8
d 2ϕ j
dx dx 2 dϕ j v= dx
dobija se: 3 m dϕ j dϕ i dϕ j 3 m dx − ∫ 3ϕ i dx = 0 i=1,...m (**) ϕ − ∫ ∑T ∑ Tj dx i 0 0 j =1 j dx dx j =1 0 3
dϕ , pa se za ϕ može izabrati linearna funkcija dx (čime se obezbeđuje da je integraljenje jednostavno): Sistem jednačina (**) sada sadrži samo
0 x > x j+ 1, x < x j− 1 x − x j− 1 ϕj = x j− 1 ≤ x ≤ x j x j − x j− 1 x j+ 1 − x x j ≤ x ≤ x j+ 1 x j+ 1 − x j
Postupak rešavanja se može sažeti u nekoliko koraka: 1. Definisati ϕ funkcije u svakom čvoru. 2. Napisati jednačine (**) za svaki čvor (nepoznate su temperature Tj u čvorovima). 3. Obaviti neophodna diferenciranja i integraljenja tako da se dobije sistem od m linearnih algebarskih jednačina sa m nepoznatih Tj. 4. Primeniti granične uslove. 5. Rešiti dobijeni sistem algebarskih jednačina. Napominje se da je nepotrebno pisati jednačinu za čvor j=1. Umesto toga sedirektno koristi granični uslov T1=1, kao jednačina za taj čvor.
9
1.2. Metod karakteristika
Metod karakteristika se može nazvati “pre-kompjuterskim” metodom. Belgijski hidrauličar Massau je još 1889.g primenio ovaj metod pri rešavanju problema rečne hidraulike (propagacije talasa u otvorenim tokovima). Naravno, pri tome je korišćena grafička integracija. Izložiće se, uprošćeno, osnovna ideja metoda karakteristika, koristeći kao primere jednačine kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima i transport zagađivača. Same jednačine će se izvesti u poglavljima 2. i 3. Metod karakteristika se najčešće vezuje za rešavanje parcijalnih diferencijalnih jednačina hiperboličnog tipa, pa će se kao prvi primer koristiti advektivni deo jednodimenzionalne jednačine transporta zagađivača (za ovu jednačinu je već ranije u poglavlju 1. utvrđeno da ima nedvosmislen hiperbolični karakter):
∂C ∂C +u =0 ∂t ∂x gde su: x , t - nezavisno promenljive (koordinate po prostoru i vremenu), u(x, t) - parametar (brzina fluida poznata iz proračuna strujanja), C(x,t) - zavisno promenljiva (nepoznata koncentracija zagadjivača). Kao početni uslov zadaje se poznato C(x,t0), što će se obeležiti sa C0(x). Kao granični uslov zadaje se poznato C(x1,t), što će se obeležiti sa C1(t). Ovde je sa x1 obeležena tzv. “uzvodna” granica (kao na skici), pa je C1(t) tzv “uzvodni” granični uslov (hiperbolične jednačine zahtevaju granični uslov samo na “uzvodnoj” granici).
Osnovna ideja metoda karakteristika sastoji se u sledećem: pretpostavlja se da u ravni (x,t) postoje takve linije duž kojih navedena parcijalna diferencijalna jednačina postaje obična diferencijalna jednačina. Te linije (koje se nazivaju i karakteristične linije ili prosto karakteristike) se takođe opisuju običnim diferencijalnim jednačinama, tako da se početna parcijalna diferencijalna jednačina svodi na sistem od dve obične diferencijalne jednačine, čime je rešavanje umnogome olakšano. Za navedeni primer, ako se trajektorija (fluidnog delića) definiše kao:
dx =u dt pa se definicija trajektorije uvrsti u razmatranu parcijalnu diferencijalnu jednačinu, dobija se:
∂ C dx ∂ C + =0 ∂ t dt ∂ x
10
Lako se uočava da je leva strana poslednje jednačine zapravo totalni (materijalni, Lagranžijanski) izvod DC/Dt, tako da se polazna parcijalna diferencijalna jednačina svodi na sistem od dve obične diferencijalne jednačine:
DC dx = 0 što važi duž trajektorije = u. dt Dt
dx = u, a dt to je duž trajektorije delića, polazna parcijalna diferencijalna jednačina postaje obična DC diferencijalna jednačina = 0. Vidi se da su u ovom slučaju karakteristične linije zapravo Dt trajektorije fluidnih delića. Izraz “što važi duž trajektorije” ovde označava da svuda gde važi jednačina
Rešenje, dobijeno jednostavnom integracijom, glasi:
C = const. duž trajektorije
dx =u dt
što je prikazano i na skici:
dx = const.), trajektorije prave i dt paralelne linije. U opštem slučaju, ako je u ≠ const., trajektorije nisu ni prave ni paralelne. Vidi se takođe da su, ukoliko je brzina u = const. (
dx = u , čime se zapravo određuje dt trajektorija koja prolazi kroz tačku u kojoj se traži rešenje. Ako se, kao na sledećoj skici, dx pretpostavi da je u > 0 ( > 0), pa traži rešenje u tački B, onda treba odrediti trajektoriju koja dt prolazi kroz B, što se može uraditi unapred, jer je brzina u poznata. U tom slučaju rešenje je trivijalno: CB = CA, gde je CA poznato iz početnog ili graničnog uslova: Problem je sveden na integraciju jednačine
11
Jedna od prednosti metode karakteristika je i u tome što ova metoda omogućava i dobro fizičko tumačenje problema. Kao što se sa poslednje skice vidi, tačka A je uvek “uzvodno” u odnosu na tačku B, tako da “informacija uvek putuje” od “uzvodnog” ka “nizvodnom” kraju oblasti strujanja. Rešenje na nizvodnoj granici se takođe dobija na osnovu “uzvodne” informacije. Odavde je jasno zašto navedena hiperbolična jednačina ne zahteva nizvodni granični uslov. Metoda karakteristika se može primeniti i na jednačine koje nisu isključivo hiperbolične. Kao sledeći primer razmatraće se sistem od dve parcijalne diferencijalne jednačine kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima: jednačina kontinuiteta:
∂h A∂u ∂h + +u =0 ∂t B∂x ∂x
dinamička jednačina:
∂u ∂u ∂h +u +g + g (S f − S0 ) = 0 ∂t ∂x ∂x gde je: B(x,t) - širina vodnog ogledala (skica), A(x,t) - površina poprečnog preseka, h(x,t) dubina, u(x,t) - brzina, Sf - nagib trenja (linije energije), S0 - nagib dna, g - gravitaciono ubrzanje.
Navedene jednačine se prvo moraju svesti na tzv. karakteristične jednačine ili “karakteristični” oblik jednačina, tj. oblik na koji se može primeniti metoda karakteristika.
12
Pretpostaviće se, zbog jednostavnosti, da se radi o prizmatičnom kanalu pravougaonog poprečnog preseka (mada nema nikakvih principijelnih teškoća i ako se ova pretpostavka ne A uvede, samo je izvođenje nešto duže). Iz navedene pretpostavke sledi da je = h . B Uvešće se nova promenljiva c =
gh odnosno: c 2 = gh
gde je nova promenljiva c zapravo brzina propagacije talasa. Odrediće se izvodi prethodnog izraza po x i t:
∂c ∂ ∂h ∂ ( c 2 ) = 2c = ( gh) = g ∂x ∂x ∂x ∂x i:
∂c ∂ ∂h ∂ (c 2 ) = 2c = ( gh) = g ∂t ∂t ∂t ∂t
odakle sledi:
∂h c ∂c =2 ∂x g ∂x i:
∂h c ∂c =2 ∂t g ∂t Uvrstiće se dobijeni izrazi za
∂h ∂ h i u polazne jednačine, čime se zapravo ∂t ∂ x
nepoznata dubina h zamenjuje novom promenljivom, nepoznatom brzinom propagacije talasa c, tako da se dobije: jednačina kontinuiteta: 2c
∂c ∂u ∂c +c + 2u =0 ∂t ∂x ∂x
dinamička jednačina:
∂u ∂u ∂c +u + 2c + g(S f − S0 ) = 0 ∂t ∂x ∂x Sabiranjem poslednje dve jednačine, odnosno oduzimanjem jednačine kontinuiteta od dinamičke jednačine, dobija se:
∂u ∂c ∂ ( u + 2 c ) + ( u + c) + 2 ( u + c) + g( S f − S0 ) = 0 ∂t ∂x ∂x
13
odnosno:
∂u ∂c ∂ ( u − 2 c ) + ( u − c) − 2 ( u − c) + g( S f − S0 ) = 0 ∂t ∂x ∂x ili, posle sređivanja:
∂ ∂ ( u + 2c ) + ( u + c) ( u + 2 c) + g ( S f − S 0 ) = 0 ∂t ∂x odnosno:
∂ ∂ ( u − 2c ) + ( u − c) ( u − 2 c) + g ( S f − S 0 ) = 0 ∂t ∂x
Dobijenesu tzv. karakteristične jednačine ili karakteristični oblik jednačina, koji je u potpunosti ekvivalentan sa polaznim parom parcijalnih diferencijalnih jednačina, jednačinom kontinuiteta i dinamičkom jednačinom. Definisaće se takozvana: pozitivna (C+) karakteristika (trajektorija):
dx =u+c dt
negativna (C-) karakteristika (trajektorija):
dx =u−c . dt
i:
Kada se uvedene definicije uvrste u karakteristične jednačine, jednačine se svode na:
D+ dx (u + 2c) = − g ( S f − S 0 ) što važi duž C+ karakteristike = u+ c Dt dt i:
D− dx (u − 2c) = − g ( S f − S 0 ) što važi duž C- karakteristike = u− c Dt dt gde “+” i “-” u totalnim izvodima simbolički označavaju da se jednačine odnose na pozitivnu i negativnu karakteristiku. Kao što se vidi, početni sistem od dve parcijalne diferencijalne jednačine se sveo na sistem od četiri obične diferencijalne jednačine sa četiri nepoznate: u, c, x, t (gde su x, t nepoznate koordinate tačke u kojoj se računaju nepoznate u i c). U ovom primeru karakteristike zapravo predstavljaju trajektorije (putanje) propagacije poremećaja (talasa) u (x,t) ravni. Razmotriće se kako izgledaju karakteristike u slučaju mirnog odnosno burnog tečenja. Ako se pretpostavi da je u > 0, vidi se da pozitivna karakteristika uvek ima pozitini nagib, bez obzira da li je tečenje mirno ili burno (jer je brzina propagacije talasa c po definiciji pozitivan broj):
14
+
C :
c=
dx = u + c > 0 (uvek) dt
u2 U slučaju mirnog tečenja je Fr = < 1, pa je u < gh , odnosno gh gh ), tako da negativna karakteristika ima negativan nagib:
C− :
u < c (jer je
dx =u−c<0 dt
kao što se vidi na skici:
u2 > 1, pa je u > gh tako da negativna karakteristika takođe ima pozitivan nagib: U slučaju burnog tečenja je Fr =
C−:
gh , odnosno u > c ( c =
gh ),
dx = u−c > 0 dt
kao što se vidi na skici:
Razmotriće se postupak rešavanja dobijenog sistema običnih diferencijalnih jednačina na primeru mirnog tečenja (mada isti postupak važi i za burno tečenje). Potrebno je rešiti sistem od četiri obične diferencijalne jednačine:
D+ dx (u + 2c) = − g ( S f − S 0 ) što važi duž C+: =u+c Dt dt
15
D− dx (u − 2c) = − g ( S f − S 0 ) što važi duž C-: =u−c Dt dt Sve četiri jednačine su istovremeno zadovoljene u tački M (skica) u kojoj se seku pozitivna i negativna karakteristika.
Jasno je da se u ovom slučaju trajektorije (karakteristike) ne mogu odrediti unapred, jer su i u i c deo rešenja. Trajektorije koje polaze iz tačaka L i D (kao na skici) seći će se u nekoj tački M, ali su koordinate tačke M (xM, tM) nepoznate unapred. Vrednosti u tačkama L i D , tj. xL, tL, uL, cL i xD, tD, tD, uD, cD, su poznate ili iz početnog ili iz graničnih uslova. Integraljenjem
D+ duž C+ dobija se: Dt tM
( u M + 2c M ) − ( u L + 2c L ) = −
∫ g( S
f
− S 0 ) dt
tL
Integraljenjem jednačine pozitivne karakteristike C+,
dx = u + c , dobija se: dt
tM
x M − x L = ∫ ( u + c) dt tL
Integraljenjem
D− Dt
duž C- dobija se: tM
(u M − 2c M ) − ( u D − 2c D ) = − ∫ g ( S f − S 0 ) dt tD
Integraljenjem jednačine negativne karakteristike C-, tM
x M − x D = ∫ (u − c)dt tD
16
dx = u − c , dobija se: dt
Dobijen je sistem od četiri jednačine sa četiri nepoznate: uM, cM, xM, tM. Pri tome je nagib dna S0 poznat, dok je nagib linije energije (nagib trenja) Sf kvadratna funkcija brzine Sf = Sf (u2). Ako se za preostale integrale koristi npr. trapezno pravilo za integraciju dobiće se sistem od četiri nelinearne algebarske jednačine (nelinearne zbog Sf = Sf (u2)) sa četiri nepoznate. Rešavanjem sistema dobija se rešenje u proizvoljnoj tački M. Rešavanje se nastavlja tako što nove karakteristike C+ i C- započinju u bilo koje dve poznate tačke, i traži se rešenje u njihovom preseku.
Opisani postupak rešavanja je veoma tačan. Naime, jedina numerička greška uvodi sa numeričkom integracijom. Nažalost, veliki nedostatak opisanog postupka je u tome da se položaj računskih tačaka u kojima se dobija rešenje ne zna unapred, tj. ne postoji kontrola nad računskom mrežom. Za praktične zadatke je veoma važno da se rešenje dobije u izabranoj tački u prostoru i vremenu. Zbog toga su razvijene varijante metoda karakteristika sa tzv. “fiksnom” (unapred određenom, izabranom) računskom mrežom. U ovom slučaju se računske tačke unapred definišu, tako što se (x,t) prostor unapred izdeli na niz tačaka (kao na skici), pa se rešenje traži samo u poznatim računskim tačkama.
U ovom slučaju poznati su xM, tM, tL = tD, a nepoznati xL, xD. Međutim, kod “fiksne” računske mreže vrednosti u i c nisu poznati u tačkama L i D, nego samo u okolnim fiksnim računskim tačkama, pa se mora koristiti interpolacija, što unosi dodatnu numeričku grešku. Uopšte, primena metoda karakteristika za tečenje u otvorenim tokovima ima niz prednosti. Metod karakteristika daje tačnije rešenje u poređenju sa drugim metodama. Primena metoda ne zavisi od režima tečenja (mirno - burno). Konačno, metod karakteristika pomaže pri identifikaciji graničnih uslova i matematičkog karaktera jednačina za neustaljeno 17
tečenje u otvorenim tokovima. Poznato je da jednačine kojima se opisuje jednodimenzionalno neustaljeno tečenje zahtevaju dva granična uslova (o tome se detaljno raspravlja u poglavlju 3.). Kao što se vidi sa sledeće skice, da bi se našlo rešenje u tački M u slučaju mirnog tečenja, potrebno je zadati jedan granični uslov na uzvodnoj, a drugi na nizvodnoj granici, što potvrđuje da u slučaju mirnog tečenja jednačine imaju parabolični karakter. Da bi se našlo rešenje u tački M u slučaju burnog tečenja, oba granična uslova se moraju zadati na uzvodnoj granici, što potvrđuje da jednačine u slučaju burnog tečenja imaju hiperbolični karakter.
Nažalost metod karakteristika se retko primenjuje za proračun neustaljenog tečenja u otvorenim tokovima, pre svega zbog praktičnih poteškoća. Metod je komplikovan za programiranje (pogotovo za prirodne vodotoke), nepogodan za komplikovanije slučajeve kao što su mreže vodotoka, uključivanje hidrauličkih objekata - ustava, preliva i sl. - u model itd. Stoga je uobičajeno da se metod karakteristika koristi kao “etalon” pri testiranju nekog drugog metoda koji je pogodniji za praktičnu upotrebu (pri čemu se koriste pojednostavljeni test slučajevi), pa se onda taj drugi metod upotrebljava za rešavanje komplikovanijih praktičnih zadataka.
18
1.3. Metod konačnih razlika (priraštaja) Osnovna ideja metoda konačnih razlika je da se parcijalni izvodi u parcijalnim diferencijalnim jednačinama aproksimiraju (približno predstavljaju) kao konačni priraštaji koristeći vrednosti funkcija u diskretnim tačkama. Ravan (oblast) x - t (prostor - vreme) se izdeli u računsku mrežu diskretnih računskih tačaka, kao na skici:
gde indeks i označava prostor, a indeks n vreme. Pogodno je razmatrati opštu tačku xi, tn (bilo koje n i i), tako da postupci objašnjeni za opštu tačku važe i za sve ostale računske tačke:
gde je: ∆xi = xi +1 − xi ∆xi −1 = xi − xi −1 Funkcije, koje su inače neprekidno definisane u čitavoj oblasti, zamenjuju se tzv. funkcijama mreže, definisanim jedino u diskretnim računskim tačkama. Ako je f neprekidna funkcija, onda će se sa f$ obeležiti tzv. funkcija mreže. Uvešće se sledeće obeležavanje:
f$i n = f$ ( xi , t n ) ≅ f ( xi , t n ) Koristiće se razvijanje funkcija u Taylor-ov red da bi se izabrale aproksimacije izvoda i procenila greška koja se pri tome pravi.
19
Razmatraće se, za početak, prvi izvod (po prostoru). Funkcija f će se razviti u Taylorov red oko tačke (xi, tn), od (xi, tn) do (xi+1, tn):
∂ f f ( xi +1 ) = f ( xi + ∆xi ) = f ( xi ) + ∂x
xi
∂ 2f ∆xi + ∂ x2
xi
∆xi2 ∂ 3 f + 2! ∂ x3
xi
∆xi3 +⋅⋅⋅ 3!
gde je: ∆x i = x i +1 − x i
, xi +1 = xi + ∆xi
Odavde je:
∂f ∂x
xi
f ( x i + ∆x i ) − f ( x i ) ∂ 2 f = − 3 ∆xi ∂x
xi
∆x i ∂ 3 f − 3 2! ∂x
xi
∆xi2 −⋅⋅ 3!
Ako se na desnoj strani prethodnog izraza zanemare (odbace) svi članovi osim prvog, može se pisati:
∂f ∂x
= xi
f ( x i + ∆x i ) − f ( xi ) + O(∆xi ) ∆xi
gde je sa O(∆xi) obeležena greška zanemarenja (odbacivanja). Red aproksimacije se definiše kao eksponent priraštaja ∆xi kojim se množi prvi zanemareni (odbačeni) član. U ovom slučaju:
∂ 2f O( ∆xi ) = − ∂ x2
xi
∆xi ∂ 3 f − 2! ∂ x 3
xi
∆xi2 −⋅⋅⋅ 3!
pa se zanemarivanjem O(∆xi) u ovom slučaju dobija aproksimacija prvog reda. Ako ∆xi→0, aproksimacija se svodi na tačnu definiciju izvoda:
∂ f ∂x
= lim xi
∆xi → 0
f ( xi + ∆xi ) − f ( xi ) ∆xi
Međutim, u praktičnim inženjerskim zadacima ∆x nije nikada nula (niti čak blisko nuli), pa je zato ispravno pisati:
∂ f ∂x
≅ xi
f$i +1 − f$i (aproksimacija prvog reda) ∆xi
Prethodno napisano je jedna od (a ima ih više) mogućih aproksimacija prvog izvoda izražena pomoću diskretnih vrednosti u računskim tačkama. Prikazaće se još jedna moguća aproksimacija prvog izvoda. Razviće se f oko xi, prvo od xi do xi+1 (∆xi =xi+1 - xi):
20
f ( xi +1 ) = f ( xi ) +
∂ f ∂x
∆xi + xi
∂ 2f ∂ x2
∆xi2 ∂ 3 f + 2! ∂ x3
xi
xi
∆xi3 +L 3!
(1)
a zatim od xi do xi-1 (∆xi-1=xi - xi-1):
∂ f f ( xi −1 ) = f ( xi ) − ∂x
xi
∂ 2f ∆xi −1 + ∂ x2
xi
∆xi2−1 ∂ 3 f − 2! ∂ x3
xi
∆xi3−1 +L 3!
(2)
Oduzimanjem (2) od (1) dobija se:
∂ f f ( xi +1 ) − f ( xi −1 ) + ∂x
∂ 2f ( ∆xi + ∆xi−1 ) + ∂ x 2 x i
( ∆x xi
2 i
− ∆xi2−1 ) 2!
∂ 3f + ∂ x3
( ∆x xi
3 i
+ ∆xi3−1 ) 3!
+L
odnosno:
∂ f ∂x
∆xi2 − ∆xi2−1 ∆xi3 + ∆xi3−1 f ( xi +1 ) − f ( xi −1 ) ∂ 2 f ∂ 3f − − −L 2 3 ∆ ∆ 2 ∆ ∆ 3 ∆ ∆ x + x !( x x ) !( x x ) + + x x ∂ ∂ i i i i i i − 1 − 1 − 1 xi xi xi 1444444444 424444444444 3 O∆x) gde je O (∆x) greška zanemarivanja. =
Prvi član greške pokazuje da je i ovo aproksimacija prvog reda, osim ako računska mreža nije uniformna (ekvidistantna) po x (∆x =const. t.j. ∆xi =∆xi-1). U tom slučaju član sa ∂ 2f 2 2 nestaje, a prvi zanemareni član množi se sa ∆x , pa je ovo aproksimacija drugog reda: ∂x O(∆x2). Može se dakle pisati:
∂ f ∂x
≅ xi
f$i +1 − f$i −1 ∆xi + ∆xi −1
što je aproksimacija prvog reda za ∆x≠const., odnosno drugog reda za ∆x=const. Na osnovu reda aproksimacije može se utvrditi kako se smanjuje greška ako se smanjuje ∆x. Smatra se da je najveći deo greške zanemarenja sadržan u prvom zanemarenom članu. Zato, ako se ∆x prepolovi, kod aproksimacije prvog reda, gde je greška odbacivanja O(∆x), prvobitna greška se množi sa jednom polovinom, a kod aproksimacije drugog reda, gde je greška odbacivanja O(∆x2), sa jednom četvrtinom. Iz ovoga se može zaključiti da se što je aproksimacija višeg reda, sa smanjenjem ∆x prvobitna greška brže smanjuje. Međutim, red aproksimacije ništa ne govori o veličini prvobitne greške (sa prvobitnim ∆x). Uopšte, red aproksimacije nije jedini kriterijum za ocenu valjanosti određene aproksimacije. Tako aproksimacija prvog reda ne mora uvek biti lošija od aproksimacije drugog reda, jer o tome odlučuju i drugi kriterijumi (stabilnost, konvergencija), o kojima se govori u poglavlju 1.4.
21
Prvi izvod po vremenu (slično izvodu po prostoru) se može aproksimirati sa:
∂ f ∂t
≅ tn
f i n+1 − f i n ∆t
što je aproksimacija prvog reda O(∆t). Za prve izvode (po prostoru i vremenu) mogu se koristiti i opštije aproksimacije:
∂ f f$i +n1+1 − f$i n+1 f$i +n1 − f$i n ≅θ + (1 − θ ) ∆xi ∆xi ∂x gde je 0 ≤ θ ≤ 1 težinski faktor po vremenu, odnosno:
∂ f f i n +1 − f i n f i +n1+1 − f i +n1 ≅ψ + (1 − ψ ) ∆t ∆t ∂t gde je 0 ≤ Ψ ≤ 1 težinski faktor po prostoru. Napisana uopštena aproksimacija izvoda po prostoru je zapravo osrednjena vrednost aproksimacija izvoda po prostoru u dva susedna vremenska nivoa, dok je napisana uopštena aproksimacija izvoda po vremenu zapravo osrednjena vrednost aproksimacija izvoda po vremenu u dve susedne računske tačke (prostorne lokacije), što se vidi i na skici:
Postoji takođe i više mogućih aproksimacija za drugi izvod. Na primer, ako se saberu jednačine (1) i (2), uz pretpostavku da je ∆x=const., dobija se sledeća aproksimacija drugog izvoda:
∂ 2f f$i +1 − 2 f$i + f$i −1 ≅ ∆x 2 ∂ x2 što je, uz pretpostavku da je ∆x=const., aproksimacija drugog reda. 22
Može se napisati i uopštena verzija prethodne aproksimacije:
∂ 2f f$i +n1+1 − 2 f$i n +1 + f$i −n1+1 f$i +n1 − 2 f$i n + f$i −n1 ≅θ + (1 − θ ) ∂ x2 ∆x 2 ∆x 2 Pokazana aproksimacija drugog izvoda je poznata kao tzv. CrankNicholson-ova aproksimacija. Crank-Nicholson-ova aproksimacija se može izvesti i na osnovu proste logike, koristeći kombinaciju aproksimacija za prve izvode po prostoru, kao što se vidi sa skice:
∂ f ∂ x2 2
= i
∂ ∂ f ≅ ∂ x ∂ x i
∂f ∂x
− i+
1 2
∂f ∂x
i−
1 2
∆x
f i +1 − f i f − f i −1 − i f − 2 f i + f i −1 ∆x ∆x = i +1 ≅ ∆x ∆x 2
Primena metode konačnih razlika (uz korišćenje navedenih ili sličnih aproksimacija za parcijalne izvode) se može sumirati u nekoliko tačaka: 1. Parcijalni izvodi u parcijalnoj diferencijalnoj jednačini se zamene odgovarajućim aproksimacijama, koristeći diskretne vrednosti u računskim tačkama; 2. Napišu se rezultujuće “diskretne” (algebarske) jednačine za svaku računsku tačku (dobija se sistem algebarskih jednačina sa nepoznatim vrednostima u računskim tačkama); 3. Reši se sistem algebarskih jednačina uz korišćenje graničnih uslova.
23
1.4 Konsistencija, konvergencija, stabilnost Uvešće se pojmovi konsistencije, konvergencije i stabilnosti kao kriterijumi za ocenu valjanosti približnog rešenja dobijenog metodom konačnih razlika. Pri tome će se koristiti relativno grube “inženjerske” definicije. Na početku je neophodno definisati šta se obično podrazumeva pod nazivom “šema metode konačnih razlika”. Šema metode konačnih razlika podrazumeva usvojenu kombinaciju aproksimacija izvoda primenjenu na određenu parcijalnu diferencijalnu jednačinu (ili sistem jednačina). Kada se izabrana kombinacija aproksimacija izvoda primeni na određenu jednačinu, dobija se tzv. približna (algebarska) jednačina. Kao primer će se koristiti jednačina jednodimenzionalne advekcije zagađivača:
∂C ∂C +u =0 ∂t ∂x gde su: x , t - nezavisno promenljive (koordinate po prostoru i vremenu), u(x, t) - parametar (brzina fluida poznata iz proračuna strujanja), C(x,t) - zavisno promenljiva (nepoznata koncentracija zagadjivača). Pretpostaviće se zbog jednostavnosti da je u=const. Napisana jednačina zahteva:
− početni uslov tj. poznato C(x,t0=0)=C0(x), i − granični uslov na “uzvodnoj” granici tj. poznato C(x1=0,t)=C1(t). Usvojiće se, takođe kao primer, šema metode konačnih razlika kod koje se za izvod po vremenu koristi ranije izvedena aproksimacija prvog reda (zanemareni (odbačeni) članovi su obeleženi sa O(∆t)):
∂C Cin+1 − Cin ≅ ∂t i ∆t dok se za izvod po prostoru koristi ranije izvedena aproksimacija drugog reda (zanemareni (odbačeni) članovi su O(∆x2), pod uslovom da je ∆x=const.):
∂C Cin+1 − Cin−1 = 2∆x ∂xi Skica prikazuje uobičajenu grafičku predstavu usvojene šeme u (x,t) ravni:
Konsistencija znači sledeće: smatra se da je određena šema metode konačnih razlika konsistentna ako se izabrane aproksimacije parcijalnih izvoda konačnim razlikama svode na 24
parcijalne izvode (tj. ako se približna jednačina svodi na tačnu jednačinu) kada ∆t, ∆x → 0. Konsistencija se istovremeno i najlakše demonstrira. Kada se izabrana šema primeni na razmatranu parcijalnu diferencijalnu jednačinu dobija se sledeća približna (algebarska) jednačina: Cin+1 − Cin Cin+1 − Cin−1 +u =0 ∆t 2 ∆x Lako je pokazati da:
Cin+1 − Cin ∂ C = ∆t → 0 ∂t ∆t n n Ci +1 − Ci −1 ∂ C i ako ∆x → 0 onda O(∆x2) → 0 odnosno lim = ∆x → 0 2 ∆x ∂x tj. približna jednačina se svodi na tačnu jednačinu, čime je dokazano da je izabrana šema konsistentna. Zapravo, lako je pokazati da je svaka šema kod koje su aproksimacije izvoda dobijene iz Taylor-ovog reda konsistentna. ako ∆t → 0 onda O(∆t) → 0 odnosno lim
Stabilnost (ili numerička stabilnost) se može definisati na sledeći način: približno rešenje dobijeno primenom metode konačnih razlika je stabilno ako greška nastala usled aproksimacija izvoda konačnim razlikama ostaje ograničena (ne raste neprekidno) tokom vremena. Konvergencija se (opet grubo) može definisati na sledeći način: približno rešenje dobijeno metodom konačnih razlika mora da se postepeno pribličava (konvergira) tačnom rešenju kada ∆t, ∆x → 0. U protivnom numerička šema nije konvergentna. Za razliku od konsistencije, konvergencija i stabilnost se teže dokazuju (demonstriraju). U slučaju linearnih jednačina konvergencija je osigurana ako su zadovoljeni uslovi koje zahteva poznata Lax-ova teorema: “Za dobro postavljen početno-granični problem i numeričku šemu koja zadovoljava uslov konsistencije, stabilnost predstavlja potreban i dovoljan uslov za konvergenciju”. Prema tome, da bi se utvrdilo da li je zadovoljena konvergencija, neophodno je sprovesti analizu stabilnosti. Ovde će se izučavati tzv. Von Neumann-ova analiza stabilnosti, koja daje tačan uslov stabilnosti za linearne parcijalne diferencijalne jednačine sa periodičnim početnim uslovima. Za ostale slučajeve Von Neumann-ova analiza stabilnosti daje neophodan ali ne i dovoljan uslov stabilnosti, mada je u praktičnom smislu ova analiza je korisna za sve slučajeve.
Von-Neumann-ova analiza stabilnosti Analiza će se prvo primeniti na tačno rešenje C(x,t) razmatrane linearne parcijalne jednačine:
∂C ∂C +u =0 ∂t ∂x 25
(u=const.)
a zatim na približno (numeričko) rešenje, za izabranu šemu metode konačnih razlika. Osnovna ideja analize (za oba slučaja, i tačno i približno rešenje) se svodi na sledeće: Prvo, početna vrednost C(x,t0=0) (početno rešenje odnosno početni uslov) će se prikazati preko Fourier-ovog reda, u t0=0, koristeći diskretne računske tačke duž x. Zatim će se posmatrati se priraštaj (propagacija) tokom vremena (kako vreme t raste) funkcije C(x,t) koja se svodi na pomenuti Fourier-ov niz za t =0. Fourier-ov niz se može izraziti preko trigonometrijskih funkcija (sinusa i kosinusa), na primer: M
C ( x , o) =
∑
m =− M
a m sin(m
2π x ) ili 2L
M
∑
m =− M
bm cos(m
2π x) 2L
gde je: m - jedna komponenta niza, am - amplituda m-te komponente, 2L - interval (duž x) duž koga je funkcija definisana. Ako se sa λm =
2L m
definiše talasna dužina m-te komponente, onda je za komponentu
m=1: λ1=2L (što je zapravo λmax), a za komponentu m=M: λM=2∆x (što je zapravo λmin). Vidi se da je λmin zapravo definisano računskom mrežom. Naime, ukupni interval 2L je izdeljen na 2M računskih tačaka tačaka (kao na skici) sa korakom ∆x, pa je 2L 2L ∆x = odnosno λ min = = 2 ∆x . 2M M
Pošto je razmatrana parcijalna diferencijalna jednačina linearna, C(x,0) se može izraziti kao linearna kombinacija trigonomtrijskih funksija (sinusa i cosinusa). Dalje izvođenje je jednostavnije ako se u Fourier-ovom nizu trigonometrijske funkcije (sunus odnosno cosinus) zamene ekvivalentnim eksponencijalnim kompleksnim funkcijama. Radi podsećanja, kompleksna funkcija (kao i kompleksan broj) ima svoj realni i imaginarni deo, pa se kompleksna funkcija f može pisati kao f =Re(f) + jIm(f) , gde je j = − 1 26
Prema tome:
∑ am sin(m
2π 2L
∑ bm cos(m
odnosno
x)
2π 2L
x)
će se zameniti sa: M
∑A
m =− M
m
e
2π x 2L
jm
gde je Am - amplituda m -te komponente, a j = − 1 . Ako m
se
poslednji
izraz
preuredi,
koristeći
λ m=
2L , m
odnosno
uvodeći
2π 2π = = σ m , gde je σm tzv. talasni broj (zamena za talasnu dužinu), dobija se: 2L λ m M
C ( x ,0) =
∑A
m =− M
m
e
j σm x
što je zapravo početni uslov (početna vrednost) izražen preko Fourier-ovog reda, koristeći eksponencijalne kompleksne funkcije. Pošto je razmatrana parcijalna diferencijalna jednačina linearna, tj. pošto se C(x,t) može predstaviti kao linearna kombinacija komponenti Fourier-ovog niza, to je dovoljno posmatrati ponašanje (propagaciju tokom vremena) jedne komponente (proizvoljne m-te komponente), pri čemu analiza važi za bilo koju komponentu, a ukupno rešenje se dobija sabiranjem svih komponenti. Da bi se proučila propagacija m-te komponente Am e j σm x tokom vremena, definisaće se Cm(x,t) -m-ta komponenta rešenja u prostoru i vremenu - kao: Cm ( x , t ) = Am e j σm x e − j β m t gde prva eksponencijalna funkcija definiše raspored po prostoru, a druga predstavlja tzv. faktor amplifikacije (porasta) po vremenu, što se može shvatiti i kao Fourier-ov niz po vremenu, analogan onome po prostoru. Na osnovu analogije sa σ, β se može definisati kao tzv. ugaona frekvencija β m = gde je Tm perioda. Međutim, za razliku od talasnog broja σ m =
2π
λm
2π , Tm
, koji je definitivno realan
broj (jer je λm - talasna dužina - zasigurno realan broj), za βm se ništa ne zna unapred - pa βm u principu može biti i kompleksan broj. Prvo će se odrediti veza između βm i σm, tako što će zameniti: Cm ( x , t ) = Am e u jednačinu: 27
j σm x
e− j β m t
∂C ∂C +u =0 ∂t ∂x Dobija se:
A e j σm x e − j β m t ⋅ (− jβ m ) + u Am e j σm x ⋅ ( jσ m ) ⋅ e − j β m t = 0 1m4444244443 14444244443 ∂C ∂t
∂C ∂x
odnosno, nakon skraćivanja:
βm = uσm što zapravo predstavlja uslov da Cm(x,t) bude rešenje jednačine
∂C ∂C +u = 0 (ako je ovaj ∂t ∂x
uslov zadovoljen, Cm(x,t) zadovoljava jednačinu, tj. leva strana jednačine je jednaka nuli). U slučaju tačnog rešenja, pošto su i σm i brzina u realni brojevi, sledi da je i njihov proizvod βm=uσm takođe realan broj. Onda je odnosom:
βm =u σm zapravo definisana brzina propagacije komponente tačnog rešenja Cm(x,t) kroz vreme. Takođe, u slučaju tačnog rešenja
βm = u ne zavisi od komponente m, t.j. brzina σm
propagacije je ista (=u) za sve komponente Cm rešenja tačne jednačine. Prema tome, sve komponente Cm tačnog rešenja se propagiraju kroz vreme (odnosno putuju nizvodno) istom brzinom u. To znači da tačno rešenje zadržava isti oblik kao na početku, tj. nema distorzije rešenja tokom vremena, jer je u svakom trenutku raspored komponenti koje formiraju rešenje isti kao i na početku. Nadalje će se posmatrati faktor amlifikacije (porasta) po vremenu e − j β m t koji ukazuje da li će rešenje (t.j. njegova amplituda, pikovi) rasti (pojačavati se) kroz vreme ili će se rešenje rasplinjavati, odnosno opadati, kroz vreme. Odrediće se apsolutna vrednost faktora amplifikacije: e− j β m t
gde se e − j β m t može napisati kao e-jθ . Radi podsećanja, apsolutna vrednost kompleksne funkcije f je:
f =
Re2 ( f ) + I m2 ( f )
pa sledi da je:
28
e− j
θ
=
Re2 ( e − j θ ) + Ι m2 ( e − j θ )
Ako se koristi Euler-ova formula: e-jθ = cosθ - j sinθ
gde je u ovom slučaju: Re (e-jθ) = cosθ
Ιm(e-jθ) = - sinθ
dobija se:
e − j β m t = e − j θ = cos2 θ + sin 2 θ = 1 Dobijeno je da je apsolutna vrednost faktora amplifikacije identički jednaka jedinici za bilo koju komponentu Cm. Ovo znači da se bilo koja komponenta, pa prema tome i njihov zbir (ukupno rešenje) neće tokom vremena ni amplificirati (pojačavati) ni rasplinjavati. Iz svega prethodnog, za tačno rešenje jednačine
∂C ∂C +u = 0 može se zaključiti ∂t ∂x
sledeće: početna vrednost C(x,0) putuje nizvodno (propagira se kroz vreme) brzinom u, bez promene oblika (bez distorzije) i bez amplifikacije odnosno rasplinjavanja (kao na skici).
U nastavku će se razmatrati približno rešenje dobijeno primenom izabranih aproksimacija izvoda:
∂ C C$ in+1 + C$ in ≅ ∆t ∂t
∂ C C$ in+1 − C$ in−1 ≅ 2∆x ∂x
na razmatranu jednačinu, gde C$ in označava približnu vrednost funkcije u račinskoj tački (xi,tn). Rezultujuća približna (aproksimativna) jednačina glasi: C$ in+1 − C$ in−1 C$ in+1 − C$ in +u =0 ∆t 2 ∆x
Približno rešenje C$ ( xi , t n ) = C$ in će se predstaviti pomoću Fourier-ovog niza, pri čemu će se koristiti xi = i∆x i tn = n∆ t , pa se dobija da je:
29
M
C$ in =
∑ A$
m =− M
m
e
j σ$m i ∆x
$
e − j β m n ∆t
Proizvoljna m - ta komponenta se onda može izraziti kao:
( C$ )
n i m
= A$ m e
j σ$m i ∆x
$
e − j β m n ∆t
(1) $
gde je: A$ m - amplituda, e j σm i ∆x - početni raspored po prostoru (u t = 0), e − j β m n ∆t - faktor aplifikacije po vremenu. $
Zanemariće se, zbog jednostavnosti, indeks m, pa se za svako C$ u približnoj jednačini može pisati: $ $ C$ in = A$ e j σ i ∆x e − j β n ∆t $ $ C$ n+1 = A$ e j σ i ∆x e − j β ( n +1) ∆t i
C$ in+1 = A$ e C$ in−1 = A$ e
j σ$ ( i +1) ∆x j σ$ ( i −1) ∆x
(1a)
$ n ∆t
e− j β
$ n ∆t
e− j β
Kao i kod tačnog rešenja, prvo će se odrediti veza između βm i σm. Približna jednačina se može preurediti u sledeći oblik:
u∆t $ n C$ in+1 − C$ in = − ( C − C$in−1 ) 2 ∆x i +1 gde je
(2)
u∆t = Cr tzv. Courant-ov broj. ∆x Ako se uvede (1) t.j. (1a) u (2) (koristeći ea+b = eaeb) dobija se: A$ e =−
) j σ i ∆x
)
)
)
)
e − j β n ∆t e − j β ∆t − A$ e j σ i ∆x e − j β n ∆t =
(
Cr $ j σ) i ∆x j σ) ∆x − j β) n ∆t $ j σ) i ∆x − j σ) ∆x − j β) n ∆t − Ae Ae e e e e 2
)
odnosno, nakon skraćivanja: )
e − j β ∆t = 1 −
Pošto je: sin θ =
e j θ − e− j θ 2j )
e− j β
) Cr j σ) ∆x ( e − e − j σ ∆x 2
)
, poslednja jednačina se može napisati kao: ∆t
= 1 − j Cr sin( σ$ ∆x )
(3)
što je zapravo uslov koji mora biti zadovoljen da bi C$ m bilo rešenje približne jednačine (2). Odmah se zapaža da veza β i σ kod tačnog rešenja (β = uσ) sadrži jedino u = const., 30
dok veza β i σ kod približnog rešenja (jednačina (3)) pored u = const . sadrži ∆x i ∆t odnosno u∆t Cr = . ∆x Razmatraće se nadalje brzina propagacije m-te komponente približnog rešenja. Već je 2π realan broj (jer je λ - talasna dužina - realan broj), dok β, međutim, u utvrđeno da je σ =
λ
principu može biti kompleksan broj. Kod tačnog rešenja jednostavna veza β i σ (β = uσ) ukazuje da je u ovom slučaju i β- realan broj (jer su i u i σ - realni brojevi). Kod približnog rešenja se, međutim, iz jednačine (3) (veza β i σ) ne vidi da li je β realan ili kompleksan broj, pa se u ovom slučaju mora ostaviti mogućnost da β može biti kompleksan broj. Zato će se, pošto je sama brzina propagacije rešenja realan broj, ta brzina definisati kao:
Re ( β$ ) σ$
β$ ) σ$
(umesto
Iz jednačine (3) sledi:
e − j ( Re ( β ) + j I m ( β ) ) ∆t = 1 − jCr sin( σ$ ∆x ) )
)
odnosno, pošto je − j 2 = − ( − 1 ) = 1 : 2
e − j Re ( β ) ∆t e I m ( β ) ∆t = 1 − jCr sin( σ$ ∆x ) )
)
Ako se primeni Euler-ova formula, e-jθ = cosθ - j sinθ , dobija se:
[
]
cos ( Re ( β$ ) ∆t ) − j sin ( Re ( β$ ) ∆t ) e I m ( β ) ∆t = 1 − jCr sin(σ$ ∆x) )
Kompleksni brojevi su jednaki ako su im jednaki realni i imaginarni delovi, tj. ako je:
cos ( Re ( β$ ) ∆t ) e I m ( β ) ∆t = 1 )
− sin ( Re ( β$ ) ∆t ) e I m ( β ) ∆t = − Cr sin (σ$ ∆t ) )
Ako se poslednje dve jednačine podele (druga sa prvom) dobija se:
tg ( Re ( β$ ) ∆t ) = Cr sin(σ$ ∆x ) odakle se jednostavno može dobiti Re ( β$ ) :
1 Re ( β$ ) = arc tg[ Cr sin(σ$ ∆t )] ∆t Brzina propagacije m-te komponente približnog rešenja je onda jednaka:
Re ( β$m ) 1 = arc tg[ Cr sin(σ$m ∆x )] σ$m σ$m ∆t 31
što ako se pomnoži i podeli sa ∆x i koristi da jeσ$m = Re ( β$m ) = σ$
∆x 2π
λm
2π
λm
daje:
λm 2π = f Cr , arc tg Cr sin , ∆x , ∆t λ m ∆x ∆t ∆x
∆x
Iz poslednjeg izraza se vidi da će se svaka komponenta C$ m približnog rešenja propagirati različitom brzinom (≠ u), jer je brzina propagacije između ostalog funkcija i talasne dužine, koja je različita za svaku komponentu. Prema tome, u t > 0 različite komponente neće zauzimati istu poziciju kao u početnom rešenju, odnosno postojaće tzv. “fazno pomeranje” ili “fazna greška”. Tako izazvana distorzija rešenja manifestuje se kroz tzv. “numeričke oscilacije”, kao na skici.
Konačno, razmatraće se faktor amplifikacije m-te komponente približnog ) − j β m n ∆t e . Za ovo je dovoljno razmatrati samo prvi vremenski korak n = 1, t.j. rešenja ) − j β m ∆t e , jer se zaključci mogu primeniti na bilo koji vremenski korak. Za prvi vremenski korak jednačina (3) postaje: )
e − j β m ∆t = 1 − j Cr sin(σ$m ∆x ) Apsolutna vrednost faktora amplifikacije, uz podsećanje da je e − j θ =
Re2 + I m2 ,
iznosi: e
) − j β m ∆t
= 1 + Cr2 sin 2 (σ$m ∆x ) > 1 uvek 1442443 >0
Apsolutna vrednost faktora amplifikacije uvek veća od jedinice znači da će m-ta (bilo koja, odnosno svaka) komponenta približnog rešenja biti amplificirana (pojačana) u svakom vremenskom koraku, bez ograničenja. Iz ovoga sledi da je razmatrana šema uvek nestabilna. Kao zaključak, razmatrana šema ne samo da ima faznu grešku, već je i bezuslovno (uvek) nestabilna, pa stoga i neupotrebljiva.
Grafički prikaz (portret) “fazne” i “amplitudne” greške
32
Definisaće se “amplitudna” greška R1 kao odnos apsolutne vrednosti faktora amplifikacije približnog rešenja i apsolutne vrednosti faktora amplifikacije tačnog rešenja. Za razmatranu jednačinu i izabranu šemu dobija se da je: )
R1 = Ako se uvrsti σ m =
2π
λm
e − j β m ∆t 1
= 1 + Cr 2 sin 2 (σ$m ∆x )
, sledi da je:
2π λm R1 = 1 + Cr 2 sin 2 = f Cr , ∆x λ m ∆x λ m za Grafički prikaz “amplitudne” greške se dobija tako što se crtaju se krive R1 ∆x različite Courant-ove brojeve Cr. Ilustracija grafičkog prikaza “amplitudne” greške, prikazana na skici, ne odnosi se na neku konkretnu šemu.
Ako je R1 > 1 to ukazuje da se amplitude tokom vremenskog koraka amplificiraju (pojačavaju, rastu) pa je prema tome numerička šema nestabilna. Ako je R1 < 1, to znači da se tokom vremenskog koraka amplitude smanjuju (opadaju), pri čemu je numerička šema stabilna, ali se približno rešenje “rasplinjava” u poređenju sa stvarnim rešenjem. Ovo “rasplinjavanje” numeričkog rešenja u odnosu na pravo rešenje se obično naziva “numerička difuzija”, po analogiji sa stvarnom fizičkom difuzijom, koja ima slične efekte na rešenje. Konačno, ako je R1 = 1, to znači da nema ni amplifikacije ni rasplinjavanja rešenja, tj. amplitude ostaju nepromenjene. Šema može biti bezuslovno (uvek) nestabilna, bezuslovno (uvek) stabilna, i uslovno stabilna tj. stabilna samo za određene vrednosti Courant-ovih Cr brojeva. Definisaće se “fazna” greška R2 kao odnos brzine propagacije m-te komponente približnog rešenja i brzine propagacije m-te komponente tačnog rešenja. Za razmatranu jednačinu i izbranu šemu je:
33
R2 = Ako se uvrsti σ m =
Re ( β$m ) σ$m u 2π
λm
=
1 1 arc tg[ Cr sin (σ$m ∆x )] u σ$m ∆t
, sledi da je:
λm
R2 =
2π ∆x arc tg Cr sin 2π Cr λ m
λm = f Cr , x ∆x
λ m za različite Grafički prikaz “fazne” greške se dobija tako što se crtaju krive R2 ∆x Courant-ove brojeve Cr. Ilustracija grafičkog prikaza “fazne” greške, prikazana na skici, ne odnosi se na neku konkretnu numeričku šemu.
Prethodno je objašnjeno da se fazna greška manifestuje kroz tzv. numeričke oscilacije. Očigledno je da kada R2 teži jedinici, to znači da brzina približnog rešenja teži brzini tačnog rešenja, pa će i intenzitet numeričkih oscilacija biti manji. Razmatrana šema metode konačnih razlika je tzv. “eksplicitna” šema. Približna jednačina za računsku tačku i: C$ in+1 − C$ in C$ in+1 − C$ in−1 +u =0 ∆t 2 ∆x sadrži samo jednu nepoznatu: C$ in+1 (vrednosti iz vremenskog nivoa n su poznate ili iz početnog uslova ili iz prethodno sračinatog vremenskog koraka). Nepoznata C$ in+1 se može izraziti (i sračunati) direktno t.j. eksplicitno iz približne jednačine, pa se takve šeme nazivaju eksplicitne šeme. U slučaju eksplicitnih šema, jednačina za bilo koju računsku tačku može se rešiti zasebno (nezavisno) od ostalih računskih tačaka, što je jednostavno za programiranje. Takozvana “implicitna” verzija razmatrane šeme koristi sledeće aproksimacije za izvode po vremenu i prostoru:
34
∂ C C$ in+1 − C$ in ≅ ∆t ∂t ∂ C C$ in++11 − C$ in−+11 = ∂x 2 ∆x što se simbolično može prikazati skicom:
Kada se navedene aproksimacije primene na parcijalne izvode u razmatranoj jednačini, dobija se približna (algebarska) jednačina: C$ in+1 − C$ in C$ in++11 − C$ in−+11 +u =0 ∆t 2 ∆x Uočava se da približna jednačina napisana za tačku i sadrži nepoznate n +1 $ $ Ci −1 , Cin +1 , C$ in++11 , tj. sadrži nepoznate u tački i i susednim tačkama (i + 1), (i - 1), pa se ne može se direktno rešiti. Rešenje je moguće samo ako se istovremeno napišu i reše jednačine za sve računske tačke, odnosno ako se formira i reši sistem algebarskih jednačina. Dakle, kod tzv. “implicitnih” šema mora se napisati i rešiti sistem jednačina, što je komplikovanije i za programiranje. Međutim, slična anliza kakva je sprovedena za “eksplicitnu” verziju razmatrane šeme pokazuje da je “implicitna” verzija razmatrane šeme uvek (tj. bezuslovno) stabilna. Naime, apsolutna vrednost faktora amplifikacije je u ovom slučaju: 1
)
e − j β m ∆t =
1 + Cr2 sin 2 (σ$m ∆x ) 1442443 >0 14442444 3
<1
uvek
>1
pa je i amplitudna greška, za slučaj razmatrane jednačine i “implicitne” verzije razmatrane šeme, uvek manja od jedinice: R1 < 1
uvek
Ako je “amplitudna” greška uvek manja od jedinice, to znači da će svaka komponenta rešenja biti rasplinuta u svakom vremenskom koraku, u odnosu na tačno rešenje.
35
Analize pokazuju da su tzv. “eksplicitne” šeme obično “uslovno” stabilne (samo u granicama gde je zadovoljen uslov stabilnosti, odnosno samo za određene vrednosti Courantovog broja), dok su “implicitne” šeme obično “bezuslovno” stabilne. To je istovremeno i osnovni razlog za korišćenje “implicitnih” šema, mada su one komlikovanije za programiranje od “eksplicitnih” šema. Na kraju će se dodati i nekoliko praktičnih napomena: Prilikom analize stabilnosti, odnosno interpretacije rezultata te analize, mora se voditi računa o tome da svaka komponenta približnog rešenja mora biti stabilna, odnosno da jedna nestabilna komponenta može uništiti celokupno rešenje. U prethodnom tekstu je proučavana jednodimenzionalna advekcija zagađivača. Za drukčije probleme (drukčije jednačine) tačno rešenje ne mora da ima faktor amplifikacije identički jednak jedinici i brzinu propagacije identički jednaku brzini strujanja vode u, ali se može primeniti ista vrsta analize.
λm
se može analizirati ponašanje različitih (svih) komponenti Fourier∆x ovog niza. Ako se zamisli da je ∆x = const., onda se različito λm (za isto ∆x) može shvatiti kao shvatiti kao razmatranje različitih komponenti (sa različitim talasnim dužinama). Promenom
Takođe, promenom
λm ∆x
se može analizirati i ponašanje samo jedne komponente za
različito ∆x. Na primer, za isto λm (istu komponentu) povećanje smanjenje ∆x tj. kao korišćenje finije računske mreže.
36
λm ∆x
se može shvatiti kao
2. JEDNODIMENZIONALNI NEUSTALJENI TRANSPORT ZAGAĐIVAČA Proučavanje transporta zagađivača zahteva da se prethodno poznaje strujno polje (brzine), što se može postići numeričkim rešavanjem jednačina kojima se opisuje strujanje fluida. U poglavlju 2 se pretpostavlja da je strujno polje već sračunato i poznato. Jednačine kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima i njihovo numeričko rešavanje izučavaju se u poglavlju 3. Cilj ovog poglavlja je da se prouče: − jednačina transporta zagađivača, tj. poreklo i značenje oba transportna mehanizma zagađivača, a to su advekcija i difuzija, − numeričko rešavanje jednačine advekcije zagađivača primenom različitih šema metoda konačnih razlika, − numeričko rešavanje jednačine advekcije primenom metoda karakteristika, − numeričko rešavanje jednačine difuzije primenom različitih šema metoda konačnih razlika, − numeričko rešavanje kompletne jednačine transporta zagađivača (sa oba transportna mehanizma, advekcijom i difuzijom), pri čemu će se koristiti postupci (npr. hibridni metod ili metod razdvojenih operatora) koji omogućuju da se kombinuju različiti numerički metodi za različite transportne mehanizme, kao na primer metod karakteristika za advekciju i metod konačnih razlika za difuziju.
2.1. Jednačina održanja mase zagađivača (jednačina transporta zagađivača) Razmatraće se zapremina ∀, ograničena površinom S, koja sadrži vodu i zagađivač, kao na skici.
Uvode se sledeće oznake: d∀ - elementarna zapremina (voda i zagađivač) dM - ukupna masa (voda i zagađivač) u d∀ dMZ - masa zagađivača u d∀ dM ρ= gustina mešavine (vode i zagađivača) d∀ Koncentracija zagađivača može da se definiše kao:
C=
dM Z dM
bezdimenzionalna koncentracija [-],
ili: 37
C∀ =
dM Z d∀
dimenzionalna (zapreminska) koncentracija [M/L3].
Veza između dve definicije koncentracije je sledeća:
ρC=
dM dM Z dM Z = = C∀ d∀ dM d∀
Masa zagađivača se može izraziti (izračunati) pomoću dimenzionalne (zapreminske) koncentracije C∀:
dM Z =
dM Z d∀ = C∀ d∀ d∀
ili pomoću bezdimenzionalne koncentracije C:
dM Z =
dM Z dM = Cρ d∀ dM
Zbog daljih izvođenja će se ukratko ponoviti:
Reynolds - ova transportna teorema Pratiće se fluid (voda i zagađivač) mase M, koji se kreće. Zapremina ∀f, koju zauzima fluid mase M, može se tokom vremena deformisati (smanjivati, povećavati, menjati oblik). Ukupna masa zagađivača unutar posmatrane mase fluida M u ∀f (u bilo kom vremenu t) iznosi: MZ =
∫ dM
Z
=
∀f
∫ ρCd∀
∀f
Promena ukupne mase zagađivača unutar posmatrane mase fluida se može izraziti kao:
DM Z D = Dt Dt
∫ ρ Cd∀ ∀f
D označava materijalni izvod, t.j. ukazuje da se prati izabrana masa fluida. Ovo je tzv. Dt Lagranžijanski (Lagrange) pristup problemu, veoma pogodan za izvođenje jednačina kojima se opisuju zakoni održanja. gde
Uvešće se nadalje sledeće obeležavanje: ∀f (t) - zapremina koju je posmatrana masa fluida zauzimala u t ∀f (t+δt) - zapremina koju posmatrana masa fluida zauzima u t+δt.
38
Sada se ukupna promena mase zagađivača unutar posmatrane mase fluida može pisati kao:
DM Z D = Dt Dt
1 ∫ ρ Cd∀ = δlim ∫ ρ Cd∀ − ∫ ρ Cd∀ t →0 δ t ∀f ∀ f (t ) ∀ f ( t + δt )
pri čemu napisano predstavlja klasičnu definiciju izvoda. Za rešavanje praktičnih zadataka je, međutim, nepraktično da se prati izabrana masa fluida. Mnogo praktičnije da se izabere (fiksna) nepokretna zapremina i da se izučava promena mase zagađivača unutar nepokretne kontrolne zapremine ∀k ograničene površinom Sk , što je tzv. Ojlerijanski(Euler) pristup problemu. Za kontrolnu zapreminu će se izabrati zapremina koju je fluid zauzimao na počerku izučavanja - u trenutku t pa je:
∀ k (t) = ∀ f ( t) Napisanom materijalnom izvodu će se dodati i oduzeti:
1 δt → 0 δ t lim
∫ ρ Cd∀
∀ k ( t + δt
)
tako da se dobije:
1 1 DM Z = lim ∫ ρ Cd∀ − ∫ ρ Cd∀ + lim ∫ ρ Cd∀ − ∫ ρ Cd∀ (*) δt → 0 δ t Dt ∀ f ( t + δt ) ∀ k ( t + δt ) ∀ f ( t ) ≡∀ k ( t ) δt → 0 δ t ∀ k ( t + δt ) Prvo će se razmatrati prvi član sa desne strane jednačine (*). Razlika integrala u razmatranom članu se svodi na razliku (mase) fluidnih delića posmatrane mase M koji napuštaju Vk kroz (deo površine) Sk, da zauzmu novu zapreminu (na skici obeleženu sa “plus”) i (mase) “novih” fluidnih delića koji ulaze kroz (deo površine) Sk da zauzmu zapreminu (na skici obeleženu sa “minus”) koju su delići mase M napustili. Stoga napisana razlika integrala nije ništa drugo nego izlaz-ulaz tj. protok delića (koji nose neku masu zagađivača) kroz Sk:
39
1 r r 1 ρ CV ⋅ n dS δ t lim ∫ ρ Cd∀ − ∫ ρ Cd∀ = lim ∫ δt → 0 δ t ∀ f ( t + δt ) ∀ k ( t + δt ) δt → 0 δ t Sk ( t + δt ) gde r je: V - vektor brzine fluidnog delića, r nr - jedinični vektor (ort) spoljašnje normale na dS, r V ⋅ n - vektorskir proizvod koji označava komponentu upravnu na dS r r r (tamo gde je V ⋅ n > 0 delići napuštaju zapreminu, a tamo gde je V ⋅ n < 0 delići osvajaju zapreminu). Drugi član sa desne strane jednačine (*) zapravo predstavlja klasičnu definiciju izvoda:
1 ∂ lim ∫ ρ C d∀ − ∫ ρ C d∀ = ρ C d∀ δt →0 δ t ∂ t ∀ ∫(t ) ∀k (t +δt ) ∀k (t ) k što je lokalna promena (unutar nepokretne kontrolne zapremine) mase zagađivača. Ako je ∀k ne samo nepokretna već i nedeformabilna zapremina, tj. ako je ∀k(t)=∀k(t+δt)=const., onda je:
∂ ( ρ C) ∂ d∀ ρ C d∀ = ∫ ∫ ∂ t ∀ (t ) ∂t ∀ (t ) k
k
Konačno, ako je ∀k nepokretna nedeformabilna zapremina, jednačina (*) se može pisati kao: DM Z D = Dt Dt
∫ ρ Cd∀ = ∫ ∀f
∀k
r r ∂ ( ρ C) d∀ + ∫ ρ CV ⋅ n dS ∂t S
Lagranžijanski
k
Ojlerijanski
što zapravo predstavlja prelaz sa Lagranžijanskog na Ojlerijanski pristup izučavanja problema. Ovaj prelaz je poznat kao Reynolds-ova transportna teorema.
Jednačina održanja mase zagađivača Ako se sada za trenutak pretpostavi da se masa zagađivača, unutar posmatrane mase fluida M koju pratimo kroz vreme, ne menja i ako se koristi Lagranžijanski pristup, jednačina održanja mase zagađivača je veoma jednostavna:
DM z =0 Dt Ako se pomoću Reynolds-ove transportne teoreme pređe na Ojlerijanski pristup, prethodna jednačina postaje: 40
r r ∂ ( ρ C) d∀ = − ∫ ρ CV ⋅ n dS ∂t ∀ S
∫
gde leva strana jednačine predstavlja lokalnu promenu mase zagađivača u ∀, dok desna strana jednačine predstavlja ulaz-izlaz zagađivača kroz S odnosno transport zagađivača kroz S brzinama fluida, što se još naziva i transport advekcijom (konvekcijom) ili prosto advekcija (konvekcija). Kao što se vidi, koncept advekcije se javlja kada se sa Lagranžijanskog pređe na Ojlerijanski pristup izučavanju problema, tj. kada se pređe sa izučavanja promena unutar zapremine koju zauzima razmatrana masa fluida koji se kreće na izučavanje promena unutar fiksne zapremine. Nadalje će se razmatrati još jedan mehanizam (u ovom slučaju fizička pojava) koji može da promeni masu zagađivača, a to je difuzija. Započeće se sa tzv. molekularnom difuzijom. Ovde se razmatra zagađivač koji je sa fluidom (vodom) pomešan na molekularnom nivou (razne hemikalije, boje i sl.). Čak i kada fluid miruje (kada su brzine jednake nuli), slučajno (haotično) molekularno kretanje, zvano Brown-ovo kretanje, će izazvati difuziju zagađivača, odnosno transport zagađivača iz zone visoke u zonu niske koncentracije. Ako se razmatra početni raspored koncentracije zagađivača (crvene boje u vodi) kao na skici (sa početno oštrom granicom između “čiste” i “ofarbane” boje):
uočava se da se tokom vremena oštra granica ubrzo gubi (usled razmene “čiste” i “ofarbane” vode kroz granicu) tako da raspored koncentracije posle nekog vremena izgleda kao na sledećoj skici:
r Sa q m će se obeležiti protok (pronos, transport) mase zagađivača usled molekularne difuzije. Protok zagađivača usled difuzije ili, kako se još naziva, difuzioni fluks, definisan je poznatim Fick-ovim zakonom.
Fick-ov zakon glasi: protok mase zagađivača usled molekularne difuzije proporcionalan je negativnom gradijentu koncentracije:
r q m = − Dm ∇ C∀ = − Dm ∇ ( ρ C ) → →
41
gde je sa Dm obeležen faktor proporcionalnosti, tzv. koeficijent molekularne difuzije. Ako se koristi Lagranžijanski način izučavanja problema (ako se prati razmatrana masa fluida) i pošto, osim difuzije, drugih uzroka promene mase zagađivača nema, kompletan zakon održanja mase zagađivača glasi:
DM Z r r = − ∫ q m ⋅ n dS Dt Sf gde Sf - označava površinu koja okružuje zapreminu ∀f koju zauzima razmatrana masa fluida, r r r r q m ⋅ n > 0 - označava masu zagađivača koja napušta ∀f kroz Sf, q m ⋅ n < 0 - označava masu zagađivača koja ulazi u ∀f kroz Sf. Kada se sa Lagranžijanskog pređe na Ojlerijanski pristup, koristeći Reynold-ovu transportnu teoremu, zakon održanja mase zagađivača postaje: r r ∂ ( ρ C) r r d∀ = − ∫ ρ CV ⋅ n dS − ∫ q m ⋅ n dS ∂t ∀ S S
∫
gde član sa leve strane jednačine predstavlja lokalnu promenu mase zagađivača u ∀, prvi član sa desne strane predstavlja ulaz-izlaz zagađivača kroz S usled kretanja fluidnih delića r brzinama V , tj. transport zagađivača advekcijom, a drugi član sa desne strane jednačine predstavlja ulaz-izlaz zagađivača kroz S usled molekularne difuzije. Prelaskom na tenzorsko obeležavanje (kartezijanske koordinate), jednačina održanja mase zagađivača postaje:
∂ ( ρ C) r d∀ = − ∫ ρ C u j n j dS − ∫ q m j n j dS ∂t ∀ S S
∫
Površinski integrali se mogu zameniti zapreminskim, koristeći Gauss-ovu teoremu:
∫a
j
∂ aj
n j dS =
∫∂ x
∀
S
d∀
j
tako da jednačina održanja mase zagađivača postaje:
(
)
∂ ρ Cuj ∂ qm j ∂ ( ρ C) d∀ = − ∫ d∀ − ∫ d∀ ∂t ∂ xj ∀ ∀ ∀ ∂ xj
∫
Napisano predstavlja integralni oblik jednačine održanja mase zagađivača. Poslednja jednačina važi za bilo koju zapreminu ∀ pa i za elementarnu zapreminu d∀, pa se može napisati diferencijalni oblik jednačine održanja mase zagađivača:
∂ qm j ∂ ( ρ C) ∂ ρ Cuj = − + ∂t ∂ xj ∂ xj
(
42
)
Pošto je iz Fick-ovog zakona poznato da je:
r q m = − Dm ∇ ( ρ C ) → onda se komponenta molekularnog difuzionog fluksa u tenzorskom obliku može pisati kao:
q m j = − Dm
∂ ( ρ C) ∂ xj
Ako se još uvede da je ρC = C∀ - zapreminska koncentracija, jednačina održanja mase zagađivača u diferencijalnom obliku može se pisati kao:
∂ C∀ ∂ ∂ ∂ C∀ + C∀ u j = D ∂ t ∂ xj ∂ xj m ∂ xj
(
)
gde prvi član sa leve strane predstavlja lokalnu promenu, drugi član sa leve strane predstavlja advekciju, a član sa desne strane predstavlja molekularnu difuziju. U poslednjoj jednačini u i C∀ su trenutne vrednosti brzine fluida i koncentracije zagađivača. Poslednja jednačina će se prilagoditi turbulentnom tečenju, tako što će se trenutne vrednosti razdvojiti na osrednjene veličine i na fluktuacije: C∀ = C∀ + C∀′ u = u + u′
gde su: C∀ , u - veličine osrednjene po vremenu C∀′ , u ′ - fluktuacije (odstupanja trenutnih vrednosti od C∀ , u ). Osrednjena vrednost npr. brzine je definisana kao: 1 u= T
t 0+ T
∫ udt
t
0
gde je T - kratak vremenski interval, ali dovoljno dug u poređenju sa periodom fluktuacija. Na skici je prikazana promena trenutnih i osrednjenih brzina kroz vreme.
43
Ako se u jednačinama C∀ zameni sa C∀ + C∀′ i u sa u + u ′ , pa onda primene pravila osrednjavanja:
f ′ = 0,
f g = f g,
f ′ g = 0,
∂ f ∂ f = , ∂x ∂x
f g = f g + f ′g′
(gde su f, g bilo koje dve turbulentne veličine) dobija se jednačina održanja mase zagađivača za osrednjeno strujanje:
∂ C∀ ∂ ∂ ∂ ∂ C∀ + C∀ u j + C∀′ u ′j = D ∂ t ∂ xj ∂ xj ∂ xj m ∂ xj
(
)
(
)
Član: C∀′ u ′j se obično naziva član turbulentne difuzije ili protok mase zagađivača r usled turbulentne difuzije, jer ima isti efekat kao i molekularna difuzija, samo sada za V ≠ 0 . Da bi se prikazao efekat turbulentne difuzije, zanemariće se za trenutak molekularna difuzija, pa će se razmatrati kako se zamišljeni početni raspored koncentracije, prikazan na skici, transportuje nizvodno.
Ako se, zbog jednostavnosti, pretpostavi da je tečenje ustaljeno i jednoliko (uniformno), očigledno je da bi se bez člana turbulente difuzije u′ C ′ (i bez molekularne difuzije), očuvao početni oblik rasporeda zagađivača, sa oštrim granicama kao na skici, bez obzira na njegovo putovanje niz struju (transport brzinama fluida). Razmotriće se sada uticaj člana turbulentne difuzije. Na skici je prikazano stvarno ponašanje trenutnih veličina, u ovom slučaju brzine, u turbulentnom tečenju.
Očigledno je da će u turbulentnom tečenju različiti fluidni delići biti transportovani do različitih pozicija - neki preko, a neki ispod originalno oštre granice između C = 0 i C > 0. Pri tome je efekat isti kao kod molekularne difuzije - originalno oštra granica između C = 0 i C > 0 se gubi. Ovo je potpuno razumljivo ako se ima u vidu da su u ′, C ′ takođe slučajnog karaktera kao i Brown-ovo kretanje. Pri tome se, kao i kod molekularne difuzije, zagađivač transportuje od zone veće koncetracije ka zoni manje koncentracije zagađivača. 44
Zbog svega navedenog, član turbulentne difuzije se definiše po analogiji sa molekularnom difuzijom kao:
u ′ C ′ = q t = − Dt
∂C ∂x
odnosno usvaja se da je protok mase usled turbulentne difuzije proporcionalan gradijentu (osrednjene) koncentracije, gde je faktor proporcionalnosti Dt tzv. koeficijent turbulentne difuzije. Jednačina održanja mase zagađivača za osrednjeno strujanje:
∂ C∀ ∂ ∂ ∂ ∂ C∀ + C∀ u j = − C∀′ u ′j + D ∂ t ∂ xj ∂ xj ∂ xj m ∂ xj
(
)
(
)
uz korišćenje:
C∀′ u′j = q t j = − Dt
∂ C∀ ∂ xj
i spajanje članova molekularne i turbulentne difuzije u jedinstveni difuzioni član:
∂ C∀ ∂ ∂ C∀ ∂ ∂ C∀ D Dm = Dt + ∂ xj ∂ xj ∂ xj ∂ xj ∂ xj ∂ xj ∂
(gde je D = Dm + Dt - koficijent difuzije) i uz izostavljanje oznaka za osrednjene u, C∀, postaje:
∂ C∀ ∂ ∂ ∂ C∀ + C∀ u) = ( D ∂ t ∂ xj ∂ xj ∂ xj Uvodeći zamenu C∀ = ρC, jednačina održanja mase zagađivača za osrednjeno strujanje postaje:
∂ ( ρ C) ∂ ∂ ∂ ( ρ C) D ρ Cuj = + ∂t ∂ xj ∂ xj ∂ xj
(
)
gde je prvi član sa leve strane član lokalne promene, drugi član sa leve strane član advekcije, a član sa desne strane član difuzije (zbirno molekularne i turbulentne). Poslednja jednačina se može dalje preurediti kao:
ρ
(
)
∂ ρ uj ∂C ∂ ρ ∂C ∂ ∂ ( ρ C) D +C + ρ uj +C = ∂t ∂t ∂ xj ∂ xj ∂ xj ∂ xj 45
odnosno:
(
∂ ρ ∂ ρ uj C + ∂ xj ∂ t
)
∂ C ∂C ∂ ∂ ( ρ C) D + ρ + uj = ∂ xj ∂ xj ∂ xj ∂t
Na osnovu jednačine kontinuiteta (jednačine održanja mase fluida):
∂ ρ ∂ ρ uj = 0 + ∂ t ∂ xj
(
)
vidi se da je prvi član prethodno napisane jednačine jednak nuli, pa se jednačina održanja mase zagađivača može pisati i u sledećem obliku:
∂ C ∂C ∂ ∂ ( ρ C) D + uj = ∂ xj ∂ xj ∂ xj ∂t
ρ
U većini praktičnih zadataka može se pretpostaviti da koncentracija zagađivača ne menja značajno gustinu fluida, t.j. da je ρ = const, pa su najčešće praktično korišćeni oblici jednačine održanja mase zagađivača:
∂C ∂C ∂ +uj = ∂t ∂ xj ∂ xj
∂C D ∂x j
odnosno:
∂C ∂ ∂ ∂C + Cuj = D ∂ t ∂ xj ∂ xj ∂ xj
(
)
Dobijena jednačina je trodimenzionalna linearna parcijalna diferencijalna jednačina (u, D - poznate funkcije), sa jednom nepoznatom C(x, y, z, t). Jednačina ima kombinovani hiperboličko-parabolički karakter, pri čemu advekcija prouzrokuje hiperbolički, a difuzija parabolički karakter. U oba slučaja potreban je početni uslov. Što se tiče graničnih uslova, advekcija zahteva granični uslov samo na uzvodnoj granici, dok difuzija zahteva granične uslove na svim granicama. Uočava se da u jednačini održanja mase zagađivača postoje dva mehanizma transporta zagađivača, različita ne samo sa fizičkog, već i sa matematičkog i numeričkog aspekta. Advekcija, koja jednačini daje hiperbolički karakter, prouzrokuje niz problema pri numeričkom rešavanju. Difuzija, koja jednačini daje parabolički karakter, ne stvara skoro nikakve probleme pri numeričkom rešavanju (mada je izvođenje tzv. “radnih” jednačina, onih koje se konačno programiraju, komplikovanije). Pošto nijedan numerički metod nije dobar za oba mehanizma, to će se advekcija i difuzija u početku razmatrati odvojeno.
46
2. 2. Jednodimenzionalna advekcija Razmatra se jednačina jednodimenzionalne advekcije:
∂C ∂C +u =0 ∂t ∂x odnosno:
∂C ∂ (uC ) + =0 ∂t ∂x U poglavlju 1.4 analizirano je kako se ponaša tačno rešenje jednodimenzionalne jednačine advekcije. Pokazano je da se, za u = const., početni raspored koncentracije zagađivača (početni uslov) transportuje nizvodno bez promene (neizmenjen) brzinom u (kao na skici).
2.2.1. Metod konačnih razlika Do sada su već pominjane dve šeme metoda konačnih razlika primenjene za rešavanje jednačine jednodimenzionalne advekcije. Ove dve razmatrane šeme prikazane su grafički na sledećoj skici:
47
Prva od dve pomenute šeme razmatrana je u poglavlju 1.4. Pokazano je da je to eksplicitna šema drugog reda po prostoru. Na osnovu analize stabilnosti zaključeno je da je šema bezuslovno nestabilna. Druga šema, koja je razmatrana u zadatku na kraju poglavlja 1.4, je tzv. "upwind" šema. Ovo je takođe eksplicitna šema, ali prvog reda po prostoru. Analiza stabilnosti pokazuje da je “upwind” šema uslovno stabilna (za Courant-ov broj Cr ≤ 0 ), kao i da prouzrokuje veliku amplitudnu grešku, tj. veliku numeričku difuziju (odnosno rasplinjavanje približnog rešenja u odnosu na tačno rešenje). Ova dva primera ilustruju tipične probleme koji se javljaju prilikom primene različitih šema metoda konačnih razlika za rešavanje jednačine advekcije:
− šeme prvog reda su obično stabilne, ali unose veliku amplitudnu grešku, tj. pruzrokuju numeričku difuziju odnosno rasplinjavanje rešenja; − šeme drugog reda su često nestabilne, a ako su stabilne onda obično imaju malu amplitudnu grešku, ali unose veliku faznu grešku (različite komponente rešenja putuju različitim brzinama), što se manifestuje kao tzv. numeričke oscilacije (skica).
Postoji veliki broj šema čije je ponašanje negde između dva pomenuta ekstrema. Pokazaće se, radi ilustracije, još dva primera.
Eksplicitna "leap-frog" šema Ova šema koristi sledeće aproksimacije za izvode po vremenu i prostoru:
∂ C Cin+1 − Cin−1 ≅ 2∆t ∂t ∂ C Cin+1 − Cin−1 ≅ 2∆x ∂x što se može šematski prikazati kao na skici:
48
“Leap-frog” šema je šema drugog reda i po vremenu i po prostoru (pod uslovom da su ∆x, ∆t = const.). Analiza stabilnosti pokazuje da je šema uslovno stabilna, odnosno stabilna za Cr ≤ 1, a nestabilna za Cr > 1 (što je slučaj kod većine eksplicitnih šema). Takođe, šema nema amplitudnu grešku unutar stabilne zone, tj. za Cr ≤ 1, amplitudna greška je identički jednaka jedinici (R1 ≡ 1). Da i ova šema “boluje” od problema uobičajenih za metod konačnih razlika primenjen na advekciju, pokazaće se poređenjem fazne greške za “leap-frog” šemu (dobijene iz Von Neumann-ove analize stabilnosti): Cr sin( σ ∆t ) 1 R2 = arctg u σ ∆t ± 1 − C 2 sin 2 ( σ ∆x ) r sa faznom greškom za "upwind" šemu: C sin( σ ∆t ) 1 r R2 = arctg u σ ∆t 1 − 2C sin 2 σ ∆x r 2 gde je: uσ ∆t = Cr 2u
∆x
λm
.
Ako se za obe šeme sračuna R2 , za npr.
λm ∆x
= 0, dobijaju se sledeće vrednosti za
različite Courant-ove brojeve:
Cr =
1 0.5 0.75
upwind 1.0 1.0 1.01
leap-frog 1.0 0.95 0.96
Iz prethodne tabele se vidi da “upwind” šema (za koju se zna da ima veliku amplitudnu grešku, tj. veliku numeričku difuziju) ima malu faznu grešku, a da “leap-frog” šema (koja nema amplitudnu grešku unutar stabilne zone) ima značajnu faznu grešku, što se manifestuje kroz tzv. numeričke oscilacije. Kao rezultat numeričkih oscilacija, numeričko rešenje često “proizvodi” negativne koncentracije zagađivača, koje nemaju fizičko značenje.
49
Nadalje će razmotriti neki praktični problemi pri primeni “leap-frog” šeme. Približna jednačina, dobijena primenom “leap-frog” šeme na jednačinu advekcije, ima direktno rešenje za tačku (i, n+1) pod uslovom da su poznate sve vrednosti u dva prethodna vremenska nivoa: (n) i (n-1). Ovo zapravo predstavlja problem pri primeni šeme. Naime, kao što se vidi sa sledeće skice, potrebna su dva početna uslova (u n i n-1) da se započne proračun. Kao rešenje uobičajeno je da se koristi neka druga šema metoda konačnih razlika za prvi računski korak po vremenu (n = 1), a onda “leap-frog” šema za ostale računske korake (n = 2 i n > 2). Nedostatak ovakvog dovijanja je u tome što je potrebno koristiti dve različite šeme tokom istog proračuna. Ako je “leap-frog” šema već izabrana zbog male amplitudne greške, onda druga šema može da unese previše grešaka u rešenje za početni vremenski korak n = 1.
Postoji još jedan problem koji se javlja pri praktičnoj primeni “leap-frog” šeme. Primena “leap-frog” šeme zahteva da se poznaje nizvodni granični uslov (što se vidi i u slučaju nizvodne granice na skici), dok sama jednačina advekcije to ne traži. Zbog toga ne postoji fizičko tumačenje šta bi taj uslov mogao biti. Za ovaj praktični problem obično se koristi jedno od dva sledeća rešenja: ili se produži oblast proračuna dovoljno daleko nizvodno i koristi uslov CI+1 = 0 (što zahteva dodatno računarsko vreme i prostor), ili se na nizvodnoj granici uslovi da je gradijent koncentracije jednak nuli CI = CI-1 (što je očigledno pogrešan granični uslov koji može da "pokvari" rešenje unutar oblasti). U zaključku, može se reći da je “leap-frog” šema, za koju je utvrđeno da nema amplitudnu grešku, ali zato ima nezanemarljivu faznu grešku (numeričke oscilacije), još i donekle nezgodna za primenu.
Kombinovana implicitno-eksplicitna šema Implicitne šeme, koje se odlikuju time da približna jednačina nema direktnog rešenja za jednu računsku tačku (već se mora napisati i rešiti sistem jednačina u svakom vremenskom koraku), su obično bezuslovno stabilne. Nažalost, implicitne šeme često imaju značajnu 50
amplitudnu grešku, odnosno prouzrokuju tzv. numeričku difuziju. Nasuprot tome, eksplicitne šeme su obično uslovno stabilne. Uobičajeni uslov u∆t ≤ 1 , često nameće veoma mali računski korak po vremenu (∆t). stabilnosti, Cr = ∆x Računski korak po prostoru (∆x) je obično diktiran potrebama problema. Nagla promena morfologije korita, potreba za izučavanjem detalja u toku i sl., mogu nametnuti malo ∆x, pa se u∆t ≤ 1 , dobija malo ∆t. Kao rezultat, potrebno je puno iz uslova stabilnosti, Cr = ∆x vremenskih koraka da se stigne od t = 0 do zahtevanog krajnjeg vremena t = T. Sa druge strane (kao što je već rečeno), eksplicitne šeme drugog reda obično imaju malu amplitudnu grešku (malu numeričku difuziju), ali unose veliku faznu grešku, što prouzrokuje tzv. numeričke oscilacije. Razmatrana implicitno-eksplicitna šema nastala je kao pokušaj da se iskoriste dobre, a eliminišu loše strane implicitnih i eksplicitnih šema. Kombinovana imlicitno-eksplicitna šema sastoji se iz dva koraka (odnosno polukoraka). Implicitni (polu)korak koristi sledeće aproksimacije za izvode po vremenu i prostoru: n+
1 2
∂ C Ci − Cin ≅ ∆t / 2 ∂t n+
1 2
n+
∂ C Ci +1 − Ci −1 ≅ 2∆x ∂x
1 2
kao što prikazano na skici:
Kada se navedene aproksimacije primene na jednačinu jednodimenzionalne advekcije, n+
1 2
n+
dobija se približna (algebarska) jednačina za tačku i koja sadrži tri nepoznate Ci −1 , Ci (u međukoraku).
1 2
n+
, Ci +1
1 2
Amplitudna greška za implicitni (polu)korak (ovde obeležena indeksom I), dobijena iz Von Neumann-ove analize stabilnosti:
51
R1 I =
1 C2 1 + r sin 2 ( σ ∆x ) 4
≤ 1 uvek
Cr2 ∆t (gde dolazi od ) ukazuje da je implicitni korak bezuslovno stabilan, ali da unosi 4 2 numeričku difuziju (rasplinjavanje rešenja). Eksplicitni (polu)korak koristi sledeće aproksimacije za izvode po vremenu i prostoru: n +1 i
n+
∂ C C − Ci ≅ ∆t / 2 ∂t n+
1
1 2
n+
1
∂ C Ci +12 − Ci −12 ≅ ∂x 2∆x kao što je prikazano na skici:
Kada se aproksimacije primene na izvode u jednačini jednodimenzionalne advekcije, dobija se približna (algebarska) jednačina za tačku i koja sadrži samo jednu nepoznatu Cin+1 . Očigledno je da je eksplicitni (polu)korak zapravo ona ista nestabilna šema koja je analizirana u poglavlju 1.4. Kao što je pokazano u poglavlju 1.4, iz Von Neumann-ove analize stabilnosti dobija se da je amplitudna greška za eksplicitni (polu)korak (ovde obeležena indeksom E):
R1
E
= 1+
Cr2 sin 2 ( σ ∆x ) 4
≥ 1 uvek
što znači da je eksplicitni korak bezuslovno (uvek) nestabilan. Međutim, za dva uzastopna vremenska koraka (polukoraka) ukupna amplitudna greška R1 je proizvod zasebnih R1 za svaki od (polu)koraka. Stoga se, za jedan vremenski korak sračunat primenom kombinovane implicitno-eksplicitne šeme, koja unutar svakog vremenskog koraka ima implicitni i eksplicitni polukorak, dobija da je ukupna amplitudna greška:
52
R1 = R1 I ⋅ R1
E
≡1
tj. ukupna amplitudna greška je identički jednaka jedinici za sve Courant-ove (Cr) brojeve i sve σ∆x. Iz ovoga se zaključuje da je kombinovana implicitno-eksplicitna šema bezuslovno stabilna i da nema amplitudne greške odnosno “rasplinjavanja” rešenja (R1 ≡ 1), što su veoma dobre karakteristike. Nažalost, ukupna fazna greška za kombinovanu implicitno-eksplicitnu šemu (koja se dobija kao zbir faznih grešaka za svaki od zasebnih (polu)koraka) iznosi: R2 = R2 I + R2
E
=
Cr 2 arctg sin ( σ ∆x ) 2 u σ ∆t
λm
= 10 daje vrednost R2 = 0,89. Ako se uporede odgovarajuće ∆x vrednosti, može se zaključiti da je, u pogledu fazne greške, kombinovana implicitnoeksplicitna šema mnogo gora nego nego i “leap-frog” i “upwind” šema. Zapravo, kombinovana implicitno-eksplicitna šema “proizvodi” značajne numeričke oscilacije za male
što za, na primer, Cr = 1 i
vrednosti
λm ∆x
. ***
Uopšte, primenom Von Neumann-ove analize stabilnosti se može pokazati zašto skoro sve šeme metoda konačnih razlika (bilo da su implicitne ili eksplicitne, prvog ili drugog reda) daju slabe rezultate za čistu advekciju i male vrednosti vrednosti” znači
λm ∆x
λm ∆x
(iskustvo pokazuje da “male
≤ 15). Zbog toga se u sledećem poglavlju razmatra metod karakteristika
za advekciju.
53
2.2.2. Metod karakteristika Osnovna ideja metoda karakteristika je već ukratko objašnjena u poglavlju 1.2. Teži se da se iskoristi hiperbolički karakter razmatrane parcijalne diferencijalne jednačine, tako da se parcijalna diferencijalna jednačina svede na sistem običnih diferencijalnih jednačina koje se lakše rešavaju. U poglavlju 1.2 je čak, kao primer, korišćena jednačina jednodimenzionalne advekcije:
∂C ∂C +u =0 ∂t ∂x Rečeno je da, ako se trajektorija (fluidnog delića) definiše kao:
dx =u dt i ova definicija uvrsti u razmatranu parcijalnu diferencijalnu jednačinu, dobija se:
∂ C dx ∂ C + =0 ∂ t dt ∂ x Lako se prepoznaje da je leva strana poslednje jednačine zapravo totalni (materijalni) izvod koncentracije C(x,t):
dC =
∂C ∂C dt + dx ∂t ∂x
odnosno:
dC ∂ C ∂ C dx DC = + = dt ∂ t ∂ x dt Dt Prema tome, razmatrana parcijalna diferencijalna jednačina:
∂C ∂C +u =0 ∂t ∂x svodi se na dve obične diferencijalne jednačine:
DC =0 Dt
što važi duž trajektorije
dx =u dt
gde je:
DC = 0 - linearna obična diferencijalna jednačina prvog reda, a Dt dx = u ( x , t ) - kvazi-linearna obična diferencijalna jednačina, takođe prvog reda. dt
54
Ako se razmatra računska mreža, kao na skici:
onda se nepoznata C u tački (i,n+1) može odrediti na sledeći način: 1) Prvo se odredi trajektorija koja stiže u tačku i sa koordinatom xi u vreme tn+1 (tzv. "dolazna" tačka D), a polazi iz tačke P sa (nepoznatom) koordinatom xP u vreme tn (tzv. "polazna" tačka P). Ovo se postiže tako što se integrali jednačina trajektorije:
dx = u ( x, t ) dt od polazne (P) do dolazne (D) tačke trajektorije: xi
t n +1
xP
tn
∫ dx = ∫ u( x, t ) dt
odakle se dobija: t n +1
x P = xi −
∫ u ( x, t ) dt
tn
gde je xP - nepoznata x-koordinata polazne tačke P trajektorije koja stiže u xi u vreme tn+1. 2) Zatim se integrali jednačina
DC = 0 duž trajektorije. Dt
Pošto je
DC = 0 duž Dt
D
trajektorije, tj. ∫ DC = 0 , tj. C = const. duž trajektorije, sledi da je: P
CD = CP odnosno: Cin+ 1 = CPn
Kao što se vidi, problem rešavanja razmatrane parcijalne diferencijalne jednačine jednodimenzionalne advekcije zamenjen je novim problemom: naći C xnP ako su poznati samo Cin− 2 , Cin−1 , Cin , Cin+1 itd. (vrednosti Cn su poznate samo u diskretnim računskim tačkama, bilo iz početnog uslova ili iz prethodnog računskog koraka). Problem znači postaje interpolacioni problem.
55
Razmatraće se (za početak) linearna interpolacija, jer je najjednostavnija. Pretpostaviće se, zbog jednostavnosti, da je: u = const. pa je integral u izrazu za xP trivijalan, odnosno dobija se da je: xP = xi - u∆t Pretpostaviće se, takođe, da je: u∆t ≤ ∆x odakle se dobija da je: u∆t ≤ 1 tj. Cr ≤ 1 ∆x tako da trajektorija izgleda kao na skici:
Odrediće se sada CxnP linearnom interpolacijom između (i) i (i-1): C Pn − Cin−1 Cin − Cin−1 = x P − xi −1 xi − xi −1 odnosno: C Pn = Cin−1 +
x P − x i −1 n ( C − Cin−1 ) x i − x i −1 i
Pošto je xP = xi - u∆t i CPn = Cin+ 1 , iz poslednje jednačine se dobija: Cin+1 = Cin−1 +
a pošto je xi − xi −1 = ∆x , odnosno
xi − u∆t − xi −1 n ( Ci − Cin−1 ) xi − xi −1
∆x − u∆t u∆t = 1− , dobija se: ∆x ∆x
Cin+1 = Cin−1 + ( 1 − Cr )( Cin − Cin−1 ) 56
gde je: Cr =
u∆t , ili: ∆x C in +1 = C in−1 + C in − C in−1 −
u∆t n C i − C in−1 ∆x
(
)
tj.: C in +1 − C in C in − C in−1 = −u ∆x ∆t Poslednja jednačina je približna (algebarska) jednačina dobijena primenom metoda karakteristika sa linearnom interpolacijom na jednačinu jednodimenzionalne advekcije. Uočava se da se istovetna približna (algebarska) jednačina dobija i primenom eksplicitne “upwind” šeme metoda konačnih razlika na jednodimenzionalnu jednačinu advekcije (što je analizirano u zadatku na kraju poglavlja 1.4.). Može se zaključiti da je primena metoda karakteristika sa linearnom interpolacijom na jednačinu jednodimenzionalne advekcije potpuno ekvivalentna primeni “upwind” šeme metoda konačnih razlika na istu jednačinu. Razmotriće se sada poznate karakteristike "upwind" šeme metoda konačnih razlika (utvrđene na osnovu Von Neumann-ove analize stabilnosti u zadatku na kraju poglavlja 1.4.), ali u novom svetlu (sa stanovišta metoda karakteristika sa linearnom interpolacijom). Iz Von Neumann-ove analize stabilnosti poznato je da je za Cr > 1 "upwind" šema nestabilna. Sa stanovišta metoda karakteristika, Cr > 1 znači da se C Pn određuje ekstrapolacijom (skica), što opet znači da vrednost C Pn nije ograničena vrednostima Cin , Cin− 1 , pa i greška može biti neograničena, što je zapravo i definicija numeričke nestabilnosti.
Analiza stabilnosti (analiza amplitudne greške) pokazuje da za Cr < 1 "upwind" šema uvodi numeričku difuziju koja "raspli njava" rešenje. Ako se problem numeričke difuzije razmatra u svetlu primene metoda karakteristika sa linearnom interpolacijom, može se zaključiti da linearna interpolacija "oduzima" masu kod lokalnih maksimuma, a "dodaje" masu kod lokalnih minimuma - t.j. teži da preraspodeli masu od maksimuma ka minimumu, što se vidi sa sledeće skice:
57
a to proizvodi numerički efekat analogan fizičkom dejstvu difuzije:
Takođe, analiza stabilnosti pokazuje da za Cr = 1 "upwind" šema daje tačno rešenje, (jer su i amplitudna i fazna greška identički jednake jedinici, R1 = R2 ≡ 1 , što znači da nema ni amplitudne ni fazne greške). Tačno rešenje za Cr = 1 se sa stanovišta metode karakteristika sa interpolacijom objašnjava time što za Cr = 1 nema interpolacione greške. Za Cr = 1 je u∆t = ∆x, odnosno: x P = x i − u ∆t = x i −1
pa nema interpolacije, tj. dobija se tačno rešenje: C Pn = Cin−1
Očigledno je da linearna interpolacija kod metoda karakteristika izaziva niz numeričkih problema. Zato je logičan pokušaj da se metod karakteristika poboljša tako što će se C Pn odrediti interpolacijom višeg reda. Razmatraće se, na primer, kvadratna interpolacija. Pretpostaviće se da je Cn kvadratna funkcija između i, i-1, i-2 (kao na skici):
tj.: C n = A1 + A2 x + A3 x 2
Odavde se može pisati da je:
58
Cin = A1 + A2 xi + A3 xi2 Cin−1 = A1 + A2 xi −1 + A3 xi2−1 Cin− 2 = A1 + A2 xi − 2 + A3 xi2− 2
gde su Cin , Cin−1 , Cin− 2 poznati, pa se mogu sračunati nepoznati koeficijenti A1, A2, A3 i odrediti: C Pn = A1 + A2 x P + A3 x P2
Na sličan način se može koristiti kubni interpolacioni polinom, sa četiri konstante određene tako što se zahteva da polinom prolazi kroz četiri tačke: i-2, i-1, i, i+1. Ista ideja (polinomi višeg reda u vremenskom nivou n) se koristi da se izvedu eksplicitne šeme višeg reda (po prostoru) metoda konačnih razlika. Prvo se definiše, recimo, kvadratni polinom u okolini tačke i tako što se koriste diskretne vrednosti Cn u i i okolnim tačkama, a zatim se definišu aproksimacije izvoda po prostoru u tački i kao izvodi polinoma. Na primer, ako je: C n = A1 + A2 x + A3 x 2
gde su:
A1 , A2 , A3 = f ( Cin , Cin−1 , Cin−2 )
onda se aproksimacija prvog izvoda po prostoru dobija kao:
∂C ≅ A2 + 2 A3 x ∂x Slično kao što je pokazano da je primena metoda karakteristika sa linearnom interpolacijom ekvivalentna primeni “upwind” šeme metoda konačnih razlika, tako se može pokazati da je primena metoda karakteristika sa interpolacijom višeg reda ekvivalentna primeni eksplicitnih šema višeg reda metoda konačnih razlika. Nažalost, pokazalo se da eksplicitne šeme višeg reda metoda konačnih razlika uglavnom smanjuju amplitudnu grešku (numeričku difuziju), ali zato povećavaju faznu grešku (numeričke oscilacije, negativne koncentracije). Ovakvo ponašanje šema višeg reda metoda konačnih razlika se, u svetlu metoda karakteristika sa interpolacijom višeg reda, može objasniti time da se korišćenjem više tačaka za konstruisanje polinoma višeg reda, koriste i informacije koje su sve dalje (fizički) od C Pn . Postoji, međutim, način da se konstruiše lokalni interpolacioni polinom višeg reda između dve tačke, ako se pri tome koriste ne samo poznate koncentracije, već i izvodi koncentracija u te dve tačke. Na ovoj ideji je zasnovan takozvani:
Holly-Preissmann-ov metod Holly-Preissmann-ov metod (ili, kako se još naziva metod “dve tačke”) je ekvivalentan sa šemom četvrtog reda metoda konačnih razlika, što će se kasnije i pokazati.
59
Razmatraće se dve susedne računske tačke, i i i-1, kao na skici:
Predstaviće se Cn kao kubni polinom između i, i-1: C n = A1 + A2 x + A3 x 2 + A4 x 3
odakle se može dobiti izvod:
∂C n = Cx n = A2 + 2 A3 x + 3 A4 x 2 ∂x ∂ Cn = CX n . gde je zbog jednostavnosti uvedeno novo obeležavanje: ∂x Nepoznati koeficijenti polinoma A1, ... , A4 se mogu odrediti koristeći poznate koncentracije i njihove izvode u dve susedne tačke: Cin , Cin−1 , CX in , CX in−1 . Može se pisati da je: Cin = A1 + A2 xi + A3 xi2 + A4 xi3 Cin−1 = A1 + A2 xi −1 + A3 xi2−1 + A4 xi3−1 CX in = A2 + 2 A3 xi + 3 A4 xi2 CX in−1 = A2 + 2 A3 xi −1 + 3 A4 xi2−1
Napisano predstavlja sistem od četiri jednačine sa četiri nepoznate A1, ... , A4 (uz pretpostavku da su C, CX poznati u vremenskom nivou tn). Za dalja izvođenja je jednostavnije da se koristi lokalni koordinatni sistem, kao na skici:
60
xi − x . Vidi se da je za x = xi lokalna koordinata r = 0, a za x = xi −1 lokalna x i − x i −1 koordinata r = 1.
gde je: r =
Pokazaće se da lokalna koordinata ima i određeno fizičko značenje. koordinata tačke P (skica) je: rP =
Lokalna
xi − x P x i − x i −1
gde je xi − xi −1 = ∆x , pa se dobija da je: xi − x P = rP ∆x
(a)
Prethodno je već pokazano da se, za u = const., integraljenjem
dx = u od P do D dobija: dt
x P = x i − u ∆t
odnosno: x i − x P = u ∆t
(b)
Iz jednačina (a) i (b) se dobija da je: r P ∆x = u ∆t
t.j. da je:
rP =
u∆t = Cr . ∆x
Znači, za u = const., lokalna koordinata tačke P je isto što i lokalni Courant-ov broj ( rP ≡ Cr ). Nadalje će se raspored koncentracije Cn predstaviti kao kubni polinom u lokalnim koordinatama (odakle se može dobiti i izvod koncentracije u lokalnim koordinatama) :
61
C n ( r ) = A1′ + A2 ′ r + A3′ r 2 + A4 ′ r 3
∂ Cn = A2 ′ + 2 A3′ r + 3 A4 ′ r 2 ∂r
(*)
gde su sa A1′ ,..., A4 ′ obeleženi koeficijenti polinoma u lokalnom koordinatnom sistemu, da bi se ukazalo sa su ovi koeficijenti različiti od koeficijenata A1, ... , A4. Prethodno je pretpostavljeno je da je poznat izvod
∂C ∂C , a ne , pa se zato mora se ∂x ∂r
naći veza između ova dva izvoda. Pošto je r = r(x) t.j. x = x(r), onda je:
∂ C ∂ C ∂ r ∂ C 1 = = − ∂ x ∂ r ∂ x ∂ r xi − xi −1 tj.:
∂C ∂C = − ( x i − x i −1 ) ∂r ∂x
odakle se za poznato
∂C ∂C može odrediti . ∂x ∂r
Prema tome, poznato je: za r = 0 ( x = xi ) :
C ( r ) = Cin
za r = 1 ( x = xi −1 ) :
C ( r ) = Cin−1
n
n
i
∂ Cn = −( xi − xi −1 ) CX in ∂r
i
∂ Cn = −( xi − xi −1 ) CX in−1 . ∂r
Kada se poznate vrednosti uvrste u jednačine (*) dobija se sistem od četiri jednačine sa A1′ ,..., A4 ′ . Rešenjem sistema dobijaju se koeficijenti četiri nepoznate n n n n A′ = f ( Ci , Ci −1 , CX i , CX i −1 ) . Polinom se može preurediti u sledeći oblik: C n ( r ) = a1Cin−1 + a 2 Cin + a 3 CX in−1 + a 4 CX in
gde su koeficijenti: a1 = r 2 ( 3 − 2 r ) a 2 = 1 − a1 a 3 = r 2 ( 1 − r ) ( x i − x i −1 ) 2 a 4 = − r ( 1 − r ) ( x i − x i −1 )
62
Iz prethodnih izraza se vidi da su za r = 1 koeficijenti a1 = 1, a2 = a3 = a4 = 0, pa je Cn(1) = Cin−1 . Znači da, kada se tačka P poklapa sa tačkom i-1 (odnosno kada je rP = 1, što znači da je Cr =1 ), iz metode karakteristika sledi da je Cin+1 = C Pn = Cin−1 , što se i očekivalo. Dobijeni polinom je polinom četvrtog reda koji koristi lokalne informacije u svega dve susedne tačke, pod uslovom da su u vremenskom nivou tn poznate ne samo koncentracije Cn već i izvodi koncentracija CXn. Polinom omogućava da se sračuna Cn+1, ali da bi se sračunalo Cn+2 potrebno je znati i CX . Prema tome, očigledno je potrebno znati početno CXn da bi se započeo proračun, ali je takođe potrebno sračunati CXn+1 da bi se proračun nastavio. Odnosno, potrebno je računati ne samo C već i CX. n+1
Da bi se sračunalo CX, koristiće se ista osnovna ideja metode karakteristika. Jednačina jednodimenzionalne advekcije:
∂C ∂C +u =0 ∂t ∂x će se diferencirti po x, pri čemu se mogu permutovati
∂ ∂ i u prvom članu, tako da se ∂x ∂t
dobija:
∂ ∂ C ∂ u ∂ C ∂ ∂ C +u + =0 ∂ t ∂ x ∂ x ∂ x ∂ x∂ x Koristeći ranije uvedeno obeležavanje CX =
∂C , poslednja jednačina postaje: ∂x
∂ ( CX ) ∂ ( CX ) ∂u +u =− CX ∂t ∂x ∂x
(*)
Ako se trajektorija (fluidnog delića) definiše kao:
dx =u dt onda se leva strana jednačine (*) može prepoznati kao totalni (materijalni) izvod izvoda koncentracije CX:
∂ ( CX ) dx ∂ ( CX ) D ( CX ) + = dt ∂ x Dt ∂t Prema tome, parcijalna diferencijalna jednačina (*) se svodi na dve obične diferencijalne jednačine:
63
D ( CX ) ∂u =− CX Dt ∂x
što važi duž trajektorije
dx =u dt
Jednačina trajektorije:
dx =u dt je prethodno već integraljena od polazne tačke P do dolazne tačke D trajektorije, čime je određena koordinata tačke P, tj. xP (ili rP u lokalnom koordinatnom sistemu). Kada se duž trajektorije integrali: d ( CX ) ∂u =− CX dt ∂x dobija se: xi
tn+ 1
xP
tn
∂u
∫ d ( CX ) = − ∫ ∂ x CXdt
odnosno: t n +1
CX in+1 − CX Pn = −
∂u
∫ ∂ x CXdt .
tn
Integral na desnoj strani poslednje jednačine zahteva numeričku integraciju. Ako se za približnu integraciju koristi, na primer, trapezno pravilo, dobija se:
CX
n +1 i
n +1 n ∂ u 1 ∂ u n +1 n − CX = − CX P + CX i ( t n+1 − t n 2 ∂ x P ∂ x i n P
)
odakle se lako može sračunati CX in+1: n
CX in+1 = CX Pn
U poslednjem izrazu izvod brzine
t n+1 − t n ∂ u 1− 2 ∂ x P n +1
t n+1 − t n ∂ u 1+ 2 ∂ x i
∂u je poznat, tj. može se jednostavno odrediti ako je ∂x
poznata brzina u u računskim tačkama, koristeći jednostavnu tehniku konačnih razlika. Ostaje da se odredi CX Pn . Za to će se koristiti isti interpolacioni polinom kao za C Pn . Izvod koncentracije u lokalim koordinatama je:
∂ Cn = A2′ ′ + 2 A3′ r + 3 A4 ′ r 2 ∂r 64
Pošto je prethodno već pokazano da je:
∂ Cn ∂ Cn = −( xi − xi −1 ) ∂r ∂x sledi da je:
∂ Cn −( xi − xi −1 ) = A2 ′ + 2 A3′ r + 3 A4 ′ r 2 ∂x gde su A2 ′ , A3′ , A4 ′ već sračunati koeficijenti. Napisani polinom se može preurediti u sledeći oblik: CX n = b1Cin−1 + b2 Cin + b3 CX in−1 + b4 CX in gde su koeficijenti: b1 =
6r ( r − 1) xi − xi −1
b2 = −b1 b3 = r ( 3r − 2) b4 = ( r − 1) ( 3r − 1) Praktična primena Holly-Preissmann-ovog metoda se može ukratko prikazati pomoću sledećeg algoritma: 1) Primeniti poznati početni uslov da se odredi Cin , CX in za sve i (izvode ili sračunati ili proceniti konačnim razlikama). 2) Petlja po vremenu. 3) Primeniti uzvodni granični uslov da se odredi C n+1 , CX n+1 u i = 1. 4) Petlja po tačkama i = 2, ... , I (I - poslednja tačka). 5) Sračunati koordinatu xP polazne tačke P trajektorije koja stiže u i u vreme tn+1 , t n +1
tj. x P = xi −
∫ udt , što za u = const. daje x
P
= xi − udt .
tn
6) Konstruisati interpolacioni polinom četvrtog reda koristeći C n , CX n , u dve tačke koje uokviruju tačku P (tj. sračunati koeficijente a1, ... , a4, b1, ... , b4). Ako postoji mogućnost da Courant-ov broj bude veći od jedinice (Cr > 1)
65
može se koristiti uopšteni postupak. Tačke koje uokviruju P obeležiće se indeksima d (desno od P) i l (levo od P), kao na skici.
xd − x . Pri xd − xl tome treba voditi računa da se prilikom sračunavanja koeficijenata polinoma koriste Cdn , Cln , , CX dn , CX ln .
Za uopšteni slučaj je ∆x = x d − x l , a lokalna koordinata r =
7) odrediti C Pn , CX Pn iz interpolacionog polinoma. 8) Sračunati Cin+1 , CX in+1 integraljenjem obične diferencijalne jednačine duž trajektorije (koristeći pri tome, na primer, trapezno pravilo)
C
n +1
Uzvodna granica zahteva poseban tretman. Na uzvodnoj granici treba odrediti , CX n+1 u i = 1. Koncentracija Cin=+11 se određuje direktno iz poznate funkcije uzvodnog
graničnog uslova: Ci =1 ( t ) . Izvod koncentracije CX in=+11 se može odrediti iz jednačine jednodimenzionalne advekcije:
∂C ∂C +u =0 ∂t ∂x odakle sledi: CX =
∂C 1∂ C =− ∂x u ∂t
što važi svuda pa i na uzvodnoj granici. Prema tome: n +1
CX
n +1 i =1
1 ∂ C = − n+1 ui =1 ∂ t i =1
n +1
∂ C gde je izvod poznate funkcije Ci =1 ( t ) . ∂ t i =1 Ako se desi da je Curant-ov broj veći od jedinice (Cr > 1), jedna ili više trajektorija mogu seći uzvodnu granicu, kao na skici:
66
U ovakvom slučaju, jedina razlika u praktičnoj primeni metode je ta da će se umesto C Pn i CX Pn koristiti C PnG , CX PnG , gde je sa PG obeležena polazna tačka trajektorije, koja se nalazi na uzvodnoj granici (kao na skici). Prema tome, u ovom slučaju je potrebno: − odrediti PG iz proračuna trajektorije, − odrediti C PG pomoću interpolacije (linearna interpolacija je obično dovoljo dobra, mada se može koristiti i interpolacija višeg reda), − odrediti CX PG iz zahteva da jednačina advekcije:
∂C ∂C +u =0 ∂t ∂x bude zadovoljena na granici, odakle sledi: CX PG = −
1 ∂ C u PG ∂ t P
G
∂ C gde je u PG određeno interpolacijom iz poznatog polja brzina, a određeno iz poznate ∂ t P funkcije uzvodnog graničnog uslova Ci =1 ( t ) .
G
Na nizvodnoj granici granični uslov nije potreban, kao što se može zaključiti i sa skice:
Prethodno je već rečeno da je Holly-Preismann-ova metoda ekvivalentna korišćenju šeme četvrtog reda metoda konačnih razlika, što se može lako pokazati. Po Holly-Preissmann-ovom metodu, raspored koncentracije između (i) i ( i − 1) , tj. u lokalnim koordinatama između r = 0 i r = 1, predstavljen je kubnim polinomom: 67
C n = A1′ + A2 ′ r + A3′ r 2 + A4 ′ r 3 čiji izvodi su:
∂ Cn = A2 ′ + 2 A3′ r + 3 A4 ′ r 2 ∂r ∂ 2C n = 2 A3′ + 6 A4 ′ r ∂ r2 ∂ 3C n = 6 A4 ′ ∂ r3 Iz polinoma i njegovih izvoda se za r = 0 jednostavno dobijaju vrednosti koeficijenata: A1′ = Cin
∂ Cn ′ A2 = ∂r
i
1 ∂ 2Cn ′ A3 = 2 ∂ r2
i
1 ∂ 3C n ′ A4 = 6 ∂ r3
i
Koncentracija Crn , za bilo koje 0 ≤ r ≤ 1, se sada može izraziti tako što će se u polinomu r zameniti sa ∆r = r − 0 (što je moguće jer je r lokalna koordinata koja ima vrednosti između nule i jedinice): Crn = A1′ + A2 ′ ∆r + A3′ ∆r 2 + A4 ′ ∆r 3 Kada se uvrste izrazi za koeficijente A1′ ,..., A4 ′ , dobija se:
∂ Cn ∂ 2 C n ∆r 2 ∂ 3 C n ∆r 3 C =C + ∆r + + ∂r i ∂ r 2 i 2! ∂ r 3 i 3! n r
n i
Istovetan izraz za Crn može se dobiti razvijanjem funkcije u Taylor-ov red i zanemarenjem svih članova višeg reda počev od člana pomnoženog sa ∆r 4 . Iz ovoga se može zaključiti da je korišćenje Holly-Preissmann-ovog metoda ekvivalentno korišćenju šeme četvrtog reda metoda konačnih razlika. Praktična primena Holly-Preissmann-ovog metoda pokazuje da je ovaj metod veoma ∂C . Naime, uticaj pogrešnog početnog CX se tolerantan na greške u proceni početnog CX = ∂x izgubi posle par računskih koraka po vremenu. Holly-Preissmann-ov metod se može relativno lako uopštiti, tako da bude ekvivalentan korišćenju šeme šestog reda metoda konačnih razlika, ako se koristi lokalni interpolacioni polinom (između dve tačke) petog stepena, što zahteva poznate koncentracije, kao i prve i druge izvode koncentracija (C n , CX n , CXX n ) u i i i-1.
68
2.3. Jednodimenzionalna difuzija Već je rečeno da jednačina difuzije ne stvara nikakve posebne probleme pri numeričkom rešavanju, kao i da se primenom metoda konačnih razlika za rešavanje jednačine difuzije dobijaju dovoljno tačni rezultati. Stoga će se ovde razmatrati jedino primena metoda konačnih razlika za rešavanje jednačine difuzije.
2.3.1. Metod konačnih razlika Pretpostaviće se, zbog jednostavnosti, da je Dx = const. (gde je Dx koeficijent difuzije koji uključuje i molekularnu i turbulentnu difuziju, a u slučaju osrednjavanja jednačine po dubini toka ili po čitavom poprečnom preseku i uticaje tog osrednjavanja - tzv. disperziju), pa se jednačina jednodimenzionalne difuzije može pisati u sledećem obliku:
∂C ∂ 2C = Dx ∂t ∂ x2 Ako se za
∂C primeni aproksimacija prvog reda (tzv. Euler-ova aproksimacija): ∂t ∂ C Cin+1 − Cin ≅ ∆t ∂t
∂ 2C a za primeni tzv. Crank-Nicholson-ova aproksimacija (aproksimacija drugog reda pod ∂ x2 uslovom da je ∆x = const.):
∂ 2 C C in+1 − 2C in + C in−1 ≅ ∂x 2 ∆x 2 ili, još bolje, uopštena verzija Crank-Nicholsonove aproksimacije:
∂ 2C Cin++11 − 2Cin+1 + Cin−+11 Cin+1 − 2Cin + Cin−1 ( ) θ θ ≅ + 1 − ∂ x2 ∆x 2 ∆x 2 (što je u principu implicitna aproksimacija, odnosno puna implicitna aproksimacija za θ = 1, ali može postati i eksplicitna aproksimacija ako je θ = 0 ) dobija se sledeća približna jednačina: Cin+1 − Cin Cin++11 − 2Cin+1 + Cin−+11 Cin+1 − 2Cin + Cin−1 ( ) = Dx θ + D 1 − θ (*) x ∆t ∆x 2 ∆x 2 Ako je θ = 0 , jednačina (*) ima jednu nepoznatu (Cin+1 ), pa je predložena šema eksplicitna (postoji rešenje za svaku pojedinačnu računsku tačku, nezavisno od ostalih računskih tačaka). Ako je θ = 1 , jednačina (*) ima tri nepoznate (Cin−+11 , Cin+1 , Cin++11 ), pa je predložena šema 69
implicitna (mora se napisati i rešiti sistem jednačina). Detaljna analiza predložene šeme se može naći npr. u knjizi Richtmyer-a i Morton-a: "Difference Methods for Initial Valne Problem" (strane 9-22, strana 189, itd.). Greška zanemarivanja (odbacivanja) za predloženu šemu je O ( ∆t ) + O ( ∆x 2 ) , pod 1 uslovom da je ∆x = const. Ako je θ = (tj. ako je šema centrirana i u vremenu i u prostoru), 2 greška zanemarivanja jeO ( ∆t 2 ) + O ( ∆x 2 ) . Za analizu stabilnosti predložene šeme za jednačinu difuzije, može se koristiti Von Neumann-ov postupak, ali i drugi načini analize stabilnosti, s tim što je problem difuzije nešto drukčiji od problema advekcije, jer je rasplinjavanje rešenja kod difuzije normalna fizička pojava. Pošto jednačina difuzije ne stvara probleme pri nnumeričkom rešavanju, to će se ovde samo prikazati rezultati analize stabilnosti, bez ulaženja u detalje. Richtmyer i Morton su analizu stabilnosti radili na više načina, ali je zaključak bio isti. 1 Za 0 ≤ θ ≤ predložena šema je uslovno stabilna, a uslov stabilnosti je: 2
Dx
dok je za
∆t 1 2 ≤ 2 − 4θ ∆x
1 1 ≤ θ ≤ 1 šema bezuslovno stabilna. Vidi se da je zaθ = : 2 2 Dx
∆t ≤∞ ∆x
Richtmyer i Morton analiziraju i veliki broj drugih šema za jednačinu difuzije, ali je za potrebe proračuna transporta zagađivača uopštena Crank-Nicholson-ova šema dovoljno tačna. Nadalje će se razmotriti primena predložene šeme za rešavanje jednačine difuzije. Za svaku računsku tačku (i) se može napisati po jedna jednačina (*) sa tri nepoznate: Cin−+11 , Cin+1 , Cin++11 , osim za prvu (i = 1) i poslednju (i = I) računsku tačku. Znači, ako postoji I računskih tačaka (i = 1, ..., I), može se napisati I-2 jednačina (*) (za tzv. “unutrašnje” tačke i = 2, ... , I-1). Prema tome, ukupan broj jednačina je I-2, dok je broj nepoznatih I (po jedna nepoznata koncentracija u svakoj od računskih tačaka). Da bi se zatvorio sistem, dodaju se dva granična uslova, po jedan na svakoj granici. Kao što je pokazano u poglavlju 1, jednačina difuzije je parabolična jednačina i zahteva početni uslov i granične uslove na svim granicama, što u slučaju jednodimenzionalnog tečenja znači na uzvodnoj i nizvodnoj granici. Prvo će se razmatrati unutrašnje tačke i = 2, ... , I-1. Jednačina (*):
Cin+1 − Cin Dxθ Dx ( 1 − θ ) n n +1 n +1 = + Cin−+11 ) + ( Ci +1 − 2Cin + Cin−1 ) 2 ( Ci +1 − 2Ci 2 ∆t ∆x ∆x 70
se, ako se pomnoži sa ∆t i pregrupišu članovi, može pisati kao:
D ∆t n +1 D ∆t n +1 D x ∆t n +1 C i −1 + − 1 − 2θ x 2 C i + θ x 2 C i +1 = θ 2 ∆x ∆x ∆x D ∆t = −(1 − θ ) x 2 C in+1 − 2C in + C in−1 − C in ∆x
(
)
odnosno:
Pi Cin−+11 + Qi Cin+1 + Ri Cin++11 = Si (**) gde su Pi , Qi , Ri , Si poznati koeficijenti u jednačini za tačku i. Mada su u ovom slučaju koeficijenti Pi i Ri isti, u opštem slučaju kada ∆x nije konstantno, onda je i Crank-Nicholsonova aproksimacija nešto drukčija, pa je i Pi ≠ Ri . Jednačine graničnih uslova se uopšteno mogu pisati (što će se naknadno detaljnije objasniti) kao: uzvodni granični uslov (u i = 1): Qi Cin+1 + Ri Cin++11 = S i nizvodni granični uslov (u i = I):
Pi Cin−+11 + Qi Cin+1 = S i
Kada se napišu jednačine (**) za sve tačke i = 2, .. , I-1 i dodaju granični uslovi, dobija se sistem jednačina: i=1 i=2 i-1 i i+1 I-1 I
Q1C1n+1 + R1C2n+1 P2 C1n+1 + Q2 C2n+1 + R2 C3n+1
= S1 = S2
Pi −1Cin−+21 + Qi −1Cin−+11 + Ri −1Cin+1 Pi Cin−+11 + Qi Cin+1 + Ri Cin++11 Pi +1Cin+1 + Qi +1Cin++11 + Ri +1Cin++21
= S i −1 = Si = S i +1
PI −1C In−+21 + QI −1C In−+11 + R I −1C In+1 = S I −1 PI C In−+11 + QI C In+1 = SI
Da bi se rešio sistem u principu se može direktno invertovati matrica koeficijenata, ali je to jedna od najzahtevnijih operacija u pogledu vremena rada računara, koja bi se uz to morala obavljati u svakom vremenskom koraku. Vidi se, međutim, da je matrica koeficijenata tzv. tri-dijagonalna matrica. Tridijagonalna matrica koeficijenata omogućava da se izbegne invertovanje matrice tako što će se koristiti tzv.:
Thomas-ov (ili “double-sweep”) algoritam Krenuće se sa uzvodnog kraja. Ako se jednačina za i = 1 napiše u sledećem obliku:
71
C1n+1 = −
R1 n+1 S1 C + Q1 2 Q1
i uvede obeležavanje: E1 = −
R1 Q1
F1 =
S1 Q1
dobija se: C1n+1 = E1C2n+1 + F1
Dobijeno će se uvrstiti u jednačinu za i = 2, tako da se dobije:
P2 ( E1C2n+1 + F1 ) + Q2 C2n +1 + R2 C3n +1 = S 2 što se može preurediti kao:
C2n+1 =
− R2 S 2 − P2 F1 C3n +1 + P2 E1 + Q2 P2 E1 + Q2
odnosno, ako se uvede obeležavanje:
E2 =
− R2 P2 E1 + Q2
F2 =
S 2 − P2 F1 P2 E1 + Q2
kao:
C 2n+1 = E 2 C 3n +1 + F2 Uopšteno, može se pretpostaviti da se jednačina za i-1 može pisati u sledećem obliku:
Cin−+11 = E i −1Cin+1 + Fi −1 Prethodni izraz će se uvrstiti u sledeću jednačinu, za i, tako da se dobije
Pi ( E i −1 Cin+1 + Fi −1 ) + Qi Cin+1 + Ri Cin++11 = Si što se može preurediti kao:
Cin+1 = −
Ri Pi E i −1 + Qi
Cin++11 +
Si − Pi Fi −1 Pi E i −1 + Qi
Fi =
Si − Pi Fi −1 Pi E i −1 + Qi
odnosno, ako se uvede obeležavanje:
Ei =
− Ri Pi E i −1 + Qi
kao: 72
Cin+1 = E i Cin++11 + Fi Kao što se vidi, jednačine za sve tačke i = 1, 2, ... , I-1 se mogu napisati u obliku:
Cin+1 = E i Cin++11 + Fi gde je:
E1 = −
R1 S F1 = 1 Q1 Q1
za i = 1
i:
Ei = −
Ri
F1 =
Pi E i −1 + Qi
S i − Pi Fi −1 Pi E i −1 + Qi
za i = 2, ... , I-1
Uočava se rekurzivni karakter izraza za E, F: da bi se sračunale vrednosti koeficijenata za razmatranu tačku, potrebno je znati vrednosti koeficijenata za prethodnu tačku u nizu. Na ovome se i zasniva ideja Thomas-ovog (ili “double-sweep”) algoritma. Thomas-ov (“double-sweep”) algoritam se sastoji iz sledećih koraka: 1) Primeniti uzvodni granični uslov da se sračuna E 1 , F1 (ovo će se detaljnije prikazati naknadno), 2) Petlja od i = 2, ... , I-1 da se sračuna E i , Fi . Uočava se da se za I-1 dobija:
C In−+11 = E I −1C In+1 + FI −1 čime je okončan tzv. “forward sweep” ili proračun napred. 3) Koristiti nizvodni granični uslov:
PI C In−+11 + QI C In+1 = S I zajedno sa jednačinom za I-1:
C In−+11 = E I −1C In+1 + FI −1 da se sračuna C In+1 . 4) Petlja od i = I-1, ... , 1 da se sračuna Cin+1 iz:
Cin+1 = E i Cin++11 + Fi čime se završava tzv. “backward sweep” ili proračun unazad. Nije naodmet napomenuti da Thomas-ov ili “double-sweep” algoritam nije posebno vezan za Crank-Nicholson-ovu šemu ili za difuziju, već je to opšti postupak za rešavanje sistema jednačina sa tri-dijagonalnom matricom koeficijenata. 73
Za θ = 0 , koeficijenti postaju: P=R=0
Q = −1
odakle sledi da su: Ei = 0
Fi = − Si = Cin +
Dx ∆t n ( Ci −1 − 2Cin + Cin+1 ) ∆x 2
Cin+1 = Fi = Cin +
Dx ∆t n ( Ci −1 − 2Cin + Cin+1 ) ∆x 2
pa se dobija da je:
što je zapravo eksplicitna verzija razmatrane šeme primenjena na jednačinu difuzije.
Granični uslovi za jednačinu difuzije Postoje dva opšta tipa graničnih uslova za jednačinu difuzije (poznati granično-početni problem): 1) Poznata (zadata) koncentracija, što je tzv. Dirichlet-ov uslov, i 2) Poznati (zadati) izvod (gradijent) koncentracije, što je tzv. Neumann-ov uslov (uobičajeni ∂C primer je = 0 na zatvorenoj granici). ∂x Razmotriće se svi mogući slučajevi graničnih uslova: 1) Zadato C(t) u i = 1. Jednačina: Q1C1n+1 + R1C2n+1 = S1 se svodi na: C1n+1 = C( t n+1 ) odnosno vednost poznate funkcije određena za t n +1 . Jednačina graničnog uslova se svodi na oblik korišćen u algoritmu, ako su koeficijenti: Q1 = 1
R1 = 0
odakle sledi da su:
74
S1 = C( t n+1 )
E1 = −
R1 =0 Q1
F1 =
S1 = C( t n+1 ) Q1
2) Zadato C(t) u i = I. U ovom slučaju se jednostavno odredi: C In+1 = C( t n+1 )
kao vrednost poznate funkcije određene za t n +1 i započne proračun unazad (tzv. “backward sweep”): C In−+11 = E I −1C In+1 + FI −1
itd. 3) Zadato
∂C = CX ( t ) u i = 1. Da bi se zadati granični uslov mogao iskoristiti, ∂x
razviće se funkcija C( x ) u Taylor-ov niz u okolini tačke i = 1 (skica), što važi za bilo koje vreme t:
C2 = C1 +
Odavde se dobija aproksimacija za
∂ 2C ∂ x2
gde je
∂ C ∂x
∂C ∂x
x1
∆x ∂ 2 C + 1! ∂ x 2
x1
∆x 2 +... 2!
∂ 2C ui = 1: ∂ x2
x1
∂C 2 ≅ 2 C2 − C1 − ∂x ∆x
∆x x1
= CX ( t ) poznato iz graničnog uslova. x1
Novodobijena aproksimacija za
∂ 2C se primenjuje umesto Crank-Nicholson-ove ∂ x2
aproksimacije u jednačini difuzije:
∂C ∂ 2C = Dx ∂t ∂ x2 za tačku i = 1, pri čemu se dobija sledeća približna jednačina:
75
Cin+1 − Cin 2 = θ Dx C n+1 − C1n+1 − CX ( t n+1 ) ∆x + ∆t ∆x 2 2 2 C n − C1n − CX ( t n ) ∆x +( 1 − θ ) D x ∆x 2 2
(
)
(
)
Članovi u poslednjoj jednačini se mogu pregrupisati tako da se dobije: 2 n+1 2 n+1 1 + θ Dx C = + −θ Dx 2 C1 ∆t ∆x ∆x 2 2 C1n 2 n n = Dx −θ CX ( t n+1 ) ∆x + ( 1 − θ ) C2 − C1 − CX ( t n ) ∆x + ∆t ∆x 2
[
)]
(
što se lako svodi na poznati oblik jednačine uzvodnog graničnog uslova: Q1C1n+1 + R1C2n+1 = S1 odakle se dalje računaju:
E1 = −
R1 Q1
F1 =
S1 Q1
itd. Za najčešći slučaj, CX ( t ) = 0, dobija se da je koeficijent: 2 S1 = D x (1 − θ ∆x 2
C1n )( C − C ) + ∆t n 2
n 1
dok su koeficijenti Q1 i R1 isti kao i za bilo koju drugu tačku. 4) Zadato
∂C = CX ( t ) u i = I . Ista vrsta analize kao u slučaju 3) daje koeficijente ∂x
PI , R I , S I u jednačini: PI C In−+11 + R I C In+1 = S I koja se onda kombinuje sa poslednjom jednačinom “forward sweep-a”: C In−+11 = E I −1C In+1 + FI −1 da se sračuna C In+1 i započne proračun unazad (tzv. “return sweep”).
76
2.4. Jednodimenzionalna advekcija i difuzija Kompletna jednačina transporta (održanja mase) zagađivača:
∂C ∂C ∂ 2C +u = Dx ∂t ∂x ∂ x2 se može jednostavno rešiti primenom metoda konačnih razlika. Na primer, ako se kombinuje "upwind" šema za advekciju sa (uopštenom) CrankNichoslon-ovom aproksimacijom za difuzioni član, dobija se približna jednačina: Cin+1 − Cin Cin + Cin−1 Cin++11 − 2Cin+1 + Cin−+11 Cin+1 − 2Cin + Cin−1 +u = Dxθ + Dx ( 1 − θ ) ∆t ∆x ∆x 2 ∆x 2 iz koje se, ako je još i θ = 0, čak dobija i eksplicitno rešenje za Cin+ 1 . Međutim, ako se pogleda kakav efekat na krajnje rešenje ima difuzija kao fizički proces, a kakav "upwind" šema primenjena na samu advekciju (što je pokazano na skici):
postavlja se pitanje kako u ukupnom rešenju razmatrati uticaj stvarne fizičke difuzije, ako je on pomešan sa veštačkom numeričkom difuzijom, izazvanom primenom “upwind” šeme na samu advekciju. Zato će se ovde razmatrati mogućnosti da se neki poboljšan metod za advekciju (npr. Holly-Preissmann-ov metod) kombinuje sa Crank-Nicholson-ovom šemom metoda konačnih razlika za difuziju. Proučiće se dva pristupa ili metoda: 1) Metod razlomljenih koraka (fractional-step method) ili metod razdvojenih operatora (splitoperator method), i 2) Hibridni metod.
77
2.4.1. Metod razlomljenih koraka (metod razdvojenih operatora) Ako se na
∂C u jednačini advekcije-difuzije: ∂t ∂C ∂C ∂ 2C = −u + Dx ∂t ∂x ∂ x2
(*)
primeni aproksimacija prvog reda (tzv. Euler-ova aproksimacija), jednačina se može pisati kao: Cin+1 − Ci∗ Ci∗ − Cin ∂C ∂ 2C + = −u + Dx ∂x ∆t ∆t ∂ x2
(**)
gde je Ci∗ “međuvrednost” nepoznate koncentracije C. Poslednja jednačina se može razdvojiti na dve jednačine koje se onda rešavaju u dva uzastopna koraka: Ci∗ − Cin ∂C ∂C ∂C advektivni korak: = −u = −u (što kad ∆t→0 postaje ) i ∂x ∆t ∂t ∂x Cin+1 − Ci∗ ∂ 2C ∂C ∂ 2C difuzioni korak: = Dx (što kad ∆t→0 postaje = Dx ) ∆t ∂t ∂ x2 ∂ x2
pri čemu zbir tačno daje jednaačinu (**). Ovaj metod se može objasniti i na drugi način. Prvo se razmatra
∂C kao lokalna ∂t
promena C usled dejstva (učinka, delovanja) čiste advekcije: a
∂ C ∂C = −u ∂x ∂t
(a)
gde a označava advekciju. Ako se na izvod koncentracije po vremenu u jednačini (a) primeni Euler-ova aproksimacija, dobija se: a
( )
n +1
∂C Ca −Cn ≅ ∆t ∂t Zatim se razmatra
∂C kao lokalna promena usled delovanja difuzije dodate na dejstvo ∂t
(učinak) advekcije:
∂ C ∂t
a +d
∂C ∂ 2C = −u + Dx ∂x ∂ x2
78
(b)
gde a+d označava kombinovanu advekciju i difuziju. Ako se na izvod koncentracije po vremenu u jednačini (b) primeni Euler-ova aproksimacija, dobija se:
∂C ∂t
a+d
(C ) ≅
a + d n +1
−Cn
∆t
Oduzimanjem jednačine (a) od (b) dobija se:
∂ C ∂t
a +d
a
∂ C ∂ 2C − = D x ∂ x2 ∂t
gde je:
∂ C ∂t
a +d
∂ C ( C a + d ) − C n ( C a ) − C n ( C a +d ) − ( C a ) − − = ≅ ∆t ∆t ∆t ∂t n +1
a
n +1
odakle se vidi da je “međuvrednost” C ∗ zapravo vrednost ( C a ) advekcije.
n +1
n+1
n +1
usled dejstva čiste
Opet se može reći da se metod svodi na razdvajanje kompletne jednačine na dva dela, koji se onda rešavaju u dva uzastopna koraka: advektivni korak:
∂C ∂C = −u ,i ∂t ∂x
difuzioni korak:
∂C ∂ 2C = Dx ∂t ∂ x2
Kombinacija Holly-Preissmann-ovog metoda za advekciju i Crank-Nicholson-ove šeme za difuziju po metodu razlomljenih koraka Tokom svakog vremenskog koraka postoje dve faze (međukoraka). U prvoj fazi se rešava advekcija koristeći Holly-Preissmann-ov metod, što je detaljno razmatrano u poglavlju 2.2.2. Tamo je rečeno da se po ovom metodu računa se ne samo ∂C = CX , koji su potrebni zbog lokalne interolacije advekcija koncentracije C, već i izvoda ∂x višeg reda. Jednačina advekcije za CX:
∂ ( CX ) ∂ ( CX ) ∂u +u =− CX ∂t ∂x ∂x je dobijena tako što je jednačina advekcije za C diferencirana po x. U drugoj fazi se rešava difuzija koristeći Crank-Nicholson-ovu šemu. Pri tome se ∂C takođe mora računati izvod = CX , da bi proračun izvoda bio kompletan, tj. da se odredi i ∂x 79
promena CX usled difuzije. Jednačina difuzije izvoda koncentracije CX se dobija tako što se jednačina difuzije koncentracije:
∂C ∂ 2C = Dx ∂t ∂ x2 diferencira po x. Ako se zbog jednostavnosti predpostavi da je Dx = const., onda se diferenciranjem dobija:
∂ ( CX ) ∂ 2 ( CX ) = Dx ∂t ∂ x2 Jednačina difuzije za izvod koncentracije CX rešava se na isti način kao i jednačina difuzije za samu koncentraciju C (što je detaljno razmatrano u poglavlju 2.3. odnosno 2.3.1.). Ukoliko je D x ≠ const., diferenciranjem jednačine difuzije koncentracije C dobijaju se dodatni članovi u jednačini difuzije izvoda koncentracije CX, što ovu jednačinu čini nešto komplikovanijom, ali ne prouzrokuje nikakav fundamentalni problem pri njenom rešavanju. Problem može da predstavlja zadavanje graničnih uslova. Teorijski, kompletna jednačina advekcije-difuzije koncentracije C je paraboličnog karaktera, što znači da zahteva kao granični uslov ili koncentracije C ili izvode koncentracije CX, i to po svim granicama. Jednačina same advekcije (koja se rešava u advektivnom koraku) je hiperboličnog karaktera, što znači da zahteva granične uslove samo na uzvodnoj granici. Prirodni (najčešće korišćeni) granični uslov za jednačinu advekcije je zadata koncentracija C na uzvodnoj granici. Iz ovoga uslova se može dobiti i granični uslov koji na uzvodnoj granici zahteva jednačina advekcije izvoda koncentracije CX, (kao što je pokazano u poglavlju 2.2.2). Jednačina difuzije koncentracije C je paraboličnog karaktera, što znači da kao granični uslov zahteva ili koncentracije C ili izvode koncentracije CX, i to po svim granicama. Ukoliko je kao granični uslov za jednačinu same advekcije korišćena poznata koncentracija C na uzvodnoj granici, taj isti uslov se može zadati i kao granični uslov za jednačinu difuzije koncentracije C na uzvodnoj granici. Međutim, jednačina difuzije izvoda koncentracije CX, koja je takođe paraboličnog karaktera, zahteva kao granični uslov ili prve izvode koncentracije CX ili druge izvode koncentracije CXX, i to po svim granicama. Problem može nastati ako je koncentracija C već zadata na uzvodnoj granici kao prirodni granični uslov za advekciju koncentracije C. Za ovaj problem ne postoji čisto (rigorozno, strogo) teorijsko rešenje, ali praksa pokazuje da se i C i CX dobijeni na kraju advektivnog koraka uvek mogu koristiti kao potrebni granični uslov za proračun difuzije.
80
2.4.2. Hibridni metod Kombinacija Holly-Preissmann-ovog metoda za advekciju i Crank-Nicholson-ove šeme za difuziju po hibridnom metodu Razmatra se kompletna jednačina advekcije-difuzije:
∂C ∂C ∂ 2C +u = Dx ∂t ∂x ∂ x2
dx = u , leva strana jednačine (*) postaje dt
Ako se trajektorija (delića) definiše kao totalni (materijalni) izvod
(*)
DC i jednačina (*) se svodi na dve obične diferencijalne jednačine: Dt
∂ 2C DC = Dx Dt ∂ x2
što važi duž trajektorije
dx =u dt
dx = u , tako da se odredi položaj dt polazne tačke P trajektorije koja dolazi u tačku i u trenutku tn+1 (dolazna tačka D, kao na skici) : Trajektorija se definiše integraljenjem jednačine
Integraljenjem se dobija: t n +1
x P = xi −
∫ udt
tn
gde je xP prostorna koordinata polazne tačke trajektorije, u poznatom vremenu tn, kao na skici. Integraljenjem jednačine:
∂ 2C DC = Dx Dt ∂ x2 duž trajektorije dobija se: D
Cin+1 = C Pn + ∫ Dx P
81
∂ 2C dt ∂ x2
Ako se za integral na desnoj strani prethodne jednačine koristi trapezno pravilo, dobija se: n +1
n
Cin+1 = C Pn + Dx
∂ 2C ∂ 2C 2 + 2 ∂ x P ∂ x i 2
(t
n +1
− tn )
n +1
∂ 2C Za izvod u prethodnoj jednačini se prirodno može koristiti Crank2 ∂ x i Nicholson-ova aproksimacija, jer je to aproksimacija drugog izvoda u tački i, kao što je i n ∂ 2C potrebno. Ostaje da se vidi kako aproksimirati izvod 2 , što je aproksimacija u tački P . ∂ x P Jedna mogućnost je da se koristi drugi izvod kubnog interpolacionog polinoma između i i i-1, ali praksa pokazuje da ovaj pristup vodi u uslovnu stabilnost i velike numeričke greške. Druga mogućnost (za šta se pokazalo da daje mnogo tačnije rezultate) je da se i za n
∂ C jednostavno koristi Crank-Nicholson-ova aproksimacija u tački i, samo ovog puta u 2 ∂ x P vremenskom nivou tn. 2
∂ 2C u poslednjoj jednačini koristi CrankPrema tome, ako se za oba izvoda ∂ x2 Nicholson-ova aproksimacija (skica), i ako se pri integraljenju koristi uopšteno trapezno 1 pravilo (sa težinskim faktorom θ umesto ), dobija se: 2
n +1 i
C
Cin++11 − 2Cin+1 + Cin−+11 Cin+1 − 2Cin + Cin−1 = C + Dx θ + (1 − θ ) ( t n+1 − t n ∆x 2 ∆x 2 n P
)
gde je C Pn dobijeno iz kubnog interpolaciong polinoma (po Holly-Preissmann-ovom metodu):
C Pn = a1Cin−1 + a 2 Cin + a 3 CX in−1 + a 4 CX in Pretposlednja jednačina se može napisati u obliku: Pi ′ Cin−+11 + Qi ′ Cin +1 + Ri ′ Cin++11 = Si ′ 82
gde su Pi ′ , Qi ′ , Ri ′ , Si ′ poznati koeficijenti u jednačini napisanoj za tačku i ("prim" ovde ukazuje da su koeficijenti izvedeni za hibridni metod, ne za čistu difuziju). Koeficijent S ′ i
sadrži C
n P
koje se određuje interpolacijom.
Kubna interpolacija zahteva i izvode koncentracije CX, koji treba da odslikavaju kompletan proces, odnosno promenu izvoda CX ne samo usled advekcije, već i usled difuzije. Jednačina advekcije-difuzije za izvode koncentracije CX se dobija tako što se jednačina advekcije-difuzije koncentracije C:
∂C ∂C ∂ 2C +u = Dx ∂t ∂x ∂ x2 diferencira po x. Pretpostaviće se, zbog jednostavnosti, da su u, Dx = const., pa se integraljenjem dobija:
∂ ( CX ) ∂ ( CX ) ∂ 2 ( CX ) +u = Dx ∂t ∂x ∂ x2 Ako se trajektorija (delića) definiše kao postaje totalni (materijalni) izvod
dx = u , leva strana prethodne jednačine dt
D ( CX ) i prethodna jednačina se svodi na dve jednačine: Dt
D ( CX ) ∂ 2 ( CX ) = Dx Dt ∂ x2
što važi duž trajektorije
dx =u dt
D ( CX ) ∂ 2 ( CX ) = Dx duž trajektorije (koristeći isti Dt ∂ x2 postupak kao i za odgovarajuću jednačinu koncentracije C), dobija se: Integraljenjem
CX
n +1 i
jednačine
CX in++11 − 2CX in+1 + CX in−+11 CX in+1 − 2CX in + CX in−1 = CX + Dx θ + (1 − θ ) ( t n+1 − t n ∆x 2 ∆x 2 n P
)
što se može preurediti u: PX i ′ CX in−+11 + QX i ′ CX in +1 + RX i ′ CX in++11 = SX i ′ gde su PX i ′ , QX i ′ , RX i ′ , SX i ′ poznati koeficijenti u jednačini za CX. Ovde "prim" ponovo označava hibridni metod. Treba voditi računa o tome da koeficijent SX i ′ sadrži CX Pn koje se određuje kubnom interpolacijom (po Holly-Preissmann-ovom metodu). Kao konačni rezultat, dobijena su dva sistema linearnih algebarskih jednačina: jedan za koncentracije C, a drugi za izvode koncentracija CX. Thomas-ov (ili “double-sweep) algoritam je primenljiv na oba, odnosno jednačine koncentracije C se mogu svesti na oblik: Cin+1 = E i ′ Cin++11 + Fi ′ 83
gde su E ′, F ′ poznati koeficijenti, a jednačine izvoda koncentracije CX se mogu svesti na oblik: CX in+1 = EX i ′ CX in++11 + FX i ′ gde su EX ′, FX ′ poznati koeficijenti.
Granični uslovi Pobrojaće se svi zahtevani granični uslovi: − kompletna jednačina advekcije-difuzije za koncentraciju C je paraboličnog karaktera, pa zahteva granične uslove na svim granicama, što mogu biti ili koncentracije C ili izvode koncentracija CX, − kompletna jednačina advekcije-difuzije za izvod koncentracije CX je takođe paraboličnog karaktera, pa takođe zahteva granične uslove na svim granicama, što mogu biti ili prvi izvod koncentracije CX ili drugi izvod koncentracije CXX. − vrednost koncentracije C Pn u polaznoj tački trajektorije određuje se kubnom interpolacija po Holly-Preissmann-ovom metodu, što zahteva poznato C ( t ) na uzvodnoj granici, dok nizvodni granični uslov nije potreban, − vrednost izvoda koncentracije CX Pn u polaznoj tački trajektorije se takođe određuje kubnom interpolacija po Holly-Preissmann-ovom metodu, što nameće da se odredi i CX n+1 na uzvodnoj granici, dok nizvodni uslov nije potreban. Izvod koncentracije CX n+1 na uzvodnoj granici se određuje iz jednačine čiste advekcije:
∂C ∂C +u =0 ∂t ∂x odakle sledi da je na uzvodnoj granici:
CX
n +1
1 ∂ C = − n+1 u ∂t
n +1
Postavlja se pitanje kako ispuniti sve navedene zahteve. Razmotriće se zasebno uzvodna i nizvodna granica.
Uzvodna granica Iskustvo pokazuje da je, kada advekcija dominira (što obično znači u rekama i kanalima gde na uzvodnoj granici advekcija dominira u odnosu na difuziju, dok difuzija dominira u jezerima, kod strujanja izazvanog vetrom, na nepropusnim granicama, itd.) na uzvodnoj granici najbolje zadati tzv. “advektivni” granični uslov, tj. zadati koncentraciju C pa odrediti izvod koncentracije CX iz jednačine čiste advekcije: 84
CX = −
1∂C u ∂t
Zadata koncentracija C se onda može koristiti i za kubnu interpolaciju po HollyPreissmann-ovom metodu (da se odredi C Pn ) i kao uzvodni granični uslov za jednačinu advekcije-difuzije koncentracije C. Izvod koncentracije CX dobijen iz jednačine čiste advekcije na uzvodnoj granici se može koristiti i za kubnu interpolaciju po Holly-Preissmannovom metodu (da se odredi CX Pn ) i kao uzvodni granični uslov za jednačinu advekcijedifuzije izvoda koncentracije CX.
Nizvodna granica Kao što je poznato, čista advekcija ne zahteva granični uslov na nizvodnoj granici, dok kompletna jednačina advekcije-difuzije koncentracije C zahteva kao nizvodni granični uslov ili koncentraciju C ili izvod koncentracije CX, a kompletna jednačina advekcije-difuzije izvoda koncentracije CX zahteva kao nizvodni granični uslov ili sam izvod CX ili drugi izvod koncentracije CXX. Logično da se na nizvodnoj granici koristi tzv. "difuzioni" granični uslov CX, jer je on zajednički za advekciju-difuziju i koncentracije C i izvoda koncentracije CX. Pri tome se, kao ∂ 2C , iz kod čiste difuzije, koristi Taylor-ov niz da se nađe odgovarajuća aproksimacija za ∂ x2 koje sledi jednačina koja vezuje C I i C I −1 (gde je indeksom I obeležena poslednja nizvodna računska tačka), pa se ta jednačina kombinuje sa poslednjom jednačinom proračuna unapred (tzv. “forward sweep-a”) i sračuna C In+1 . U praktičnim problemima je, međutim, izvod koncentracije CX obično nepoznat, pa se kao praktično rešenje ukazuju dve mogućnosti. Jedna mogućnost je da se prosto zanemari difuzija na nizvodnoj granici, odnosno da se sračunaju C In+1 , CX In+1 samo usled advekcije, po HollyPreissmann-ovom metodu, odnosno da se ne piše kompletna jednačina advekcije-difuzije za tačku I (skica).
Druga mogućnost je da se nizvodna granica pomeri dovoljno daleko tako da se može koristiti uslov da je izvod koncentracije CX = 0 . U slučaju zatvorene (nepropusne) granice (skica), jasno je da je difuzija, odnosno transport zagađivača difuzijom, kroz takvu granicu nula, pa se kao granični uslov zadaje CX = 0 . 85
Postupak rešavanja se može prikazati pomoću sledećeg algoritma: 1) Primeniti početni uslov da se odredi Cin , CX in , za sve tačke i=1,...,I; 2) Odrediti položaj polazne tačke trajektorije: D
x P = xi − ∫ udt P
pa sračunati C Pn , CX Pn koristeći kubnu interpolaciju (prethodno sračunati HollyPreissman-ove koeficijente a1 ,..., a 4 , b1 ,..., b4 ), za tačke i=2,...,I; 3) Koristeći uzvodni granični uslov sračunati koeficijente: Q1′ , R1′ , S1′ odnosno E1′ , F1′ i QX ′ , RX ′ , SX ′ odnosno EX ′ , FX ′ 1
1
1
1
1
za tačku i=1; 4) Sprovesti proračun unapred (tzv. “forward sweep”), tj. sračunati koeficijente: Qi ′ , Ri ′ , S i ′ odnosno E i ′ , Fi ′ i QX ′ , RX ′ , SX ′ odnosno EX ′ , FX ′ i
i
i
i
i
za tačke i=2,.....,I; 5) Primeniti nizvodni granični uslov da se sračuna C In+1 , CX In+1 ; 6) Sprovesti proračun unazad (tzv. “backward sweep”), tj. sračunati Cin+1 , Cxin+1 : C n+1 = E ′ C n+1 + F ′ i
i
i +1
i
CX = EX i ′ CX in++11 + FX i ′ za sve tačke i=I-1,.....1. n +1 i
Napomene: Kada je u ≠ const., Dx ≠ const, neophodno integraljenje u prethodno ispisanim jednačinama hibridnog metoda je unekoliko složenije, ali to ne menja osnovne procedure metoda. Bilo koja razložna integracija daje dovoljno dobre rezultate, jer ovde drugi problemi imaju više uticaja na tačnost. Na primer, ako je u ≠ const., pa se za integraciju u jednačini:
86
D
x p = xi − ∫ udt P
koristi trapezno pravilo, dobija se:
x p = xi −
uin+1 + u Pn ∆t 2
Ovde je u Pn nepoznato, jer je u poznato samo u čvorovima, pa ga treba odrediti interpolacijom. Ako se koristi linearna interpolacija, što za udređivanje u daje zadovoljavajuće rezultate, sledi da je:
u Pn = uin−1 + ( uin − uin−1 )
x P − xi −1 xi − xi −1
Kada se poslednji izraz uvrsti u pretposlednju jednačinu, dobija se:
x p = xi −
x P − xi −1 ∆t n+1 n n n u + u + u − u ( ) i −1 i i −1 2 i xi − xi −1
iz čega se može dobiti eksplicitni izraz za x P .
87
3. JEDNODIMENZIONALNO NEUSTALJENO TEČENJE U OTVORENIM TOKOVIMA 3.1. Integralni oblik jednačina Integralni oblik jednačina kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima će se izvesti polazeći od integralnog oblika osnovnih jednačina mehanike fluida: jednačine održanja mase i jednačine održanja količine kretanja. Ako se koristi Lagranžijanski (Lagrange) pristup, tj. ako se prati izabrana masa fluida koja zauzima zapreminu ∀f, onda se jednačina održanja mase (jednačina kontinuiteta) jednostavno formuliše kao:
D Dt
∫ ρ d∀ = 0 ∀f
gde je: ρ - gustina fluida. Prelaskom na Ojlerijanski (Euler) pristup (Reynolds-ova transportna teorema), t.j. na razmatranje promena unutar nepokretne zapremine ∀ ograničene površinom S, jednačina održanja mase postaje:
∂ ∫ ρ d∀ + ∫ ρ u j n j dS = 0 ∂t∀ S r gde je: uj - komponenta vektora brzine v u pravcu j, a nj komponenta jediničnog vektora r r r spoljne normale n u pravcu j, tako da uj nj zapravo predstavla vektorski proizvod v ⋅ n . U prethodnoj jednačini prvi član sa leve strane označava lokalnu promenu mase u ∀, a drugi r r član označava izlaz-ulaz mase kroz S (vektorski proizvod v ⋅ n > 0 predstavlja (pozitivan) r r izlaz mase, a proizvod v ⋅ n < 0 predstavlja (negativan) ulaz mase, što se može zaključiti i sa skice).
Takođe, ako se koristi Lagranžijanski (Lagrange) pristup, tj. prati izabrana masa fluida(koja zauzima zapreminu ∀f ograničenu površinom Sf), pa na tu masu primeni drugi Njutnov zakon (koji tvrdi da je promena količine kretanja razmatrane mase jednaka sumi svih sila koje na tu masu deluju), onda se jednačina održanja količine kretanja (napisana za pravac i) može napisati u sledećem obliku:
88
D Dt
∫ ρu
i
∀f
d∀ = ∫ ρ f i d∀ − ∀f
∫ pn
i
Sf
dS + ∫ σ jid n j dS Sf
gde je: ui - komponenta vektora brzine u pravcu i, fj - komponenta vektora zapreminske sile (po jedinici mase) u pravcu i, p - pritisak (sferni deo napona), σ jid - komponenta tenzora devijatorskog dela napona, koja deluje u pravcu i na ravan sa normalom u pravcu j. U prethodnoj jednačini član sa leve strane predstavlja promenu količine kretanja, a članovi sa desne strane predstavljaju sile, i to prvi član zapreminsku silu a drugi i treći površinske sile (razdvojene na površinsku silu usled pritiska i površinsku silu usled devijatorskog dela napona). Prelaskom na Ojlerijanski pristup (Reynolds-ova transportna teorema), tj. na razmatranje promena unutar nepokretne zapremine ∀ ograničene površinom S, jednačina održanja količine kretanja postaje:
∂ ∫ ρ u d∀ = − ∫ ρ u j n j ui dS + ∫ ρ f i d∀ − ∫ p ni dS + ∫ σ jid n j dS ∂t∀ i S ∀ S S gde prvi član sa leve strane označava lokalnu promenu količine kretanja u ∀, a drugi član sa leve strane označava izlaz-ulaz količine kretanja kroz S. Članovi sa desne strane imaju isto značenje kao u prethodnoj jednačini. Razmatraće se kako se masa i količina kretanja menjaju tokom konačnog vremenskog intervala od t1 do t2, pa će se zato polazne jednačine integraliti i po vremenu: jednačina održanja mase:
∂ ∫ ∂ t ∫ ρ d∀dt = − ∫ ∫ ρ u j n j dSdt S ∀
t2
t2
t1
t1
jednačina održanja količine kretanja (za pravac i):
∂ ∫ ∂ t ∫ ρ ui d∀dt = − ∫ ∫ ρ u j n j ui dSdt + ∫ ∫ ρ f i d∀dt − ∫ ∫ p ni dSdt + ∫ ∫ σ jid n j dSdt t ∀ t S t ∀ t S t S
t2
t2
t2
t2
t2
1
1
1
1
1
Prethodno napisane jednačine su trodimenzionalne, kao što je to uostalom i tečenje u otvorenim tokovima. Jednodimenzionalno tečenje je aproksimacija koja se može primeniti samo u slučajevima kada važe određene pretpostavke. Francuski hidrauličar de St. Venant (Jean-Claude Barre de Saint-Venant) je još 187071. formulisao uslove pod kojima se neustaljeno tečenje u otvorenim tokovima može opisati jednodimenzionalnim jednačinama i izveo odgovarajuće jednačine (koje se obično nazivaju i de St. Venant-ove jednačine). Navešće se poznate de St. Venant-ove pretpostavke približno onako kako ih je sam autor formulisao (a u zagradi će se dodati i dopunska objašnjenja): 1) Nivo slobodne površine se ne menja previše naglo (što znači da strujnice nisu previše zakrivljene, ondosno da je vertikalno ubrzanje zanemarljivo, tj. da važi pretpostavka o hidrostatičkom rasporedu pritisaka). 89
2) Promena nivoa slobodne površine poprečno na struju nema značajnog uticaja na propagaciju talasa (što znači da se može pretpostaviti da je nivo slobodne površine upravno na struju horizontalan). 3) Neuniformnost brzina po poprečnom preseku struje ne utiče značajno na propagaciju talasa (što znači da se može pretpostaviti konstantna brzina po poprečnom preseku). 4) Gubici na trenje u neustaljenom tečenju ne razlikuju se značajno od onih u ustaljenom tečenju (što znači da se mogu koristiti formule za trenje koje važe za ustaljeno tečenje). 5) Prosečni pad dna u podužnom pravcu je mali (tj. cos α ≅ 1 , odnosno sin α ≅ t g α , gde je α ugao koji dno zaklapa sa horizontalom, kao na skici).
U nastavku će se izvesti jednačine za jednodimenzionalno neustaljeno tečenje u otvorenim tokovima (de St. Venant-ove jednačine), pri čemu se neće slediti originalno izvođenje, ali će se koristiti originalne de St. Venant-ove pretpostavke, kao i dodatna pretpostavka: 6) Kod većine praktičnih problema može se smatrati da je gustina ρ = const. po poprečnom preseku. Razmatraće se kontrolna zapremina ∀ između x = x1 i x = x2 i t = t1 i t = t2, kao na skici:
gde je: u - komponenta brzine upravna na poprečni presek struje (što je pod generalnom pretpostavkom jednodimenzionalnog strujanja jedina komponenta brzine različita od nule), S površina koja ograničava zapreminu ∀, A - površina poprečnog preseka struje upravno na pravac strujanja. Izvođenje će krenuti od polaznih trodimenzionalnih jednačina, integraljenih po vremenu. Prvo će se razmatrati jednačina održanja mase, član po član. Lokalna promena mase u ∀ , u intervalu t1 do t2:
90
t2
∂
∫ ∂ t ∫ ρ d∀dt ∀
t1
je, ako se pretpostavi da je ρ = const po A, ekvivalentna sa:
∂ ∫ ∂ t ∫ ρ A dxdt t x
t2
x2
1
1
odakle se (trivijalnim) integraljenjem po vremenu dobija: x2
x2
∫ ρ A dx
x1
− ∫ ρ A dx x1
t2
t1
što se sažetije može napisati kao: x2
∫ ( ρ A)
x1
t2
− ( ρ A) t dx 1
(Izlaz-ulaz) mase kroz S, u intervalu t1 do t2: t2
∫∫ ρ u
j
n j dSdt
t1 S
je, pošto su uz generalnu pretpostavku jednodimenzionalnog tečenja sve komponente brzine osim one upravne na poprečni presek struje jednake nuli, odnosno pošto od čitave površine S kod jednodimenzionalnog strujanja protok postoji samo kroz poprečne preseke A1 i A2, jednak:
∫ ∫ ρ u dA − ∫ ρ u dA dt t1 A1 A2
t2
što se, ako se koriste pretpostavke da su ρ i u konstantni po poprečnom preseku pa je ∫ ρ udA = ρ u ∫ dA = ρ u A , može pisati kao: A
A
t2
∫ ( ρ u A)
t1
x2
− ( ρ u A) x dt 1
Konačno, integralni oblik jednačine održanja mase za jednodimenzionalno neustaljeno tečenje u otvorenim tokovima glasi: x2
t2
x1
t1
∫ ( ρ A) t2 − ( ρ A) t1 dx + ∫ ( ρ u A) x2 − ( ρ u A) x1 dt = 0
Razmatraće se nadalje jednačina održanja količine kretanja, takođe član po član. 91
Lokalna promena količine kretanja u ∀, u intervalu t1 do t2: t2
∂
∫ ∂ t ∫ ρ u d∀dt i
∀
t1
je, pod pretpostavkom da su ρ i u = const po A (gde je u = ui, i - pravac upravan na presek struje, za koji se jadnačina i piše) , jednaka:
∂ ∫ ∂ t ∫ ρ u A dxdt t x
t2
x2
1
1
odakle se (trivijalnim) integraljenjem po vremenu dobija: x2
x2
∫ ρ u A dx − ∫ ρ u A dx
x1
x1
t2
t1
što se sažetije može pisati kao: x2
∫ ( ρ u A)
x1
t2
− ( ρ u A) t dx 1
(Izlaz-ulaz) količine kretanja kroz S, u intervalu t1 do t2: t2
∫∫ ρ u
j
n j ui dSdt
t1 S
je, pošto su uz generalnu pretpostavku jednodimenzionalnog tečenja sve komponente brzine osim one upravne na poprečni presek struje jednake nuli, odnosno pošto od čitave površine S kod jednodimenzionalnog strujanja protok postoji samo kroz poprečne preseke A1 i A2, jednak:
2 2 dt ρ u dA − ρ u dA ∫ ∫ ∫ t1 A A 2 1
t2
što se, ako se koriste pretpostavke da su ρ i u konstantni po poprečnom preseku pa je ∫ ρ u 2 dA = ρ u 2 ∫ dA = ρ u 2 A , može pisati kao: A
A
t2
∫ ( ρ u A) 2
t1
x2
− ( ρ u 2 A) x dt 1
Dejstvo zapreminskih sila na fluid u ∀ , tokom intervala t1 do t2: t2
∫∫ ρ f
i
t1 ∀
92
d∀dt
će se izraziti pod pretpostavkom da se od zapreminskih sila ovde razmatra jedino sila gravitacije (težina). Onda je vektor zapreminske sile:
r r f =g r gde je g - gravitaciono ubrzanje koje deluje vertikalno naniže. Pošto se kod jednodimenzionalnog tečenja jednačine pišu za pravac tečenja, to je ovde potrebna komponenta gravitacionog ubrzanja u pravcu x: f i = gi = g x
što se, uz pretpostavku da je podužni nagib dna mali ( cos α ≅ 1 ), može pisati kao: g x = g sin α ≅ g tg α = − g
gde su uvedene oznake: Zdna - kota dna, − g
∂ Zd n a = gS 0 ∂x
∂ Zd n a = S 0 - lokalni nagib dna. ∂x
Dejstvo zapreminskih sila na fluid u ∀ tokom intervala t1 do t2 (pisano za pravac tečenja) se onda može izraziti kao: t2
∫∫ ρ gS
0
d∀dt
t1 ∀
što je, pod pretpostavkom da je ρ = const. po A, ekvivalentno sa: t 2 x2
∫∫ρ
g S 0 A dxdt
t1 x1
Dejstvo sile pritiska na fluid u ∀ , tokom intervala t1 do t2: t2
− ∫ ∫ p ni dS dt t1 S
se, u jednodimenzionalni jednačinama napisanim za pravac strujanja, svodi na dejstvo sile pritiska u pravcu strujanja x. 93
Pogodno je za dalje izvođenje da se dejstvo sile pritiska u x pravcu razdvoji na dva dela: a) dejstvo sile pritiska na poprečne preseke u x1 i x2, i b) reakcija na dejstvo sile pritiska na omotač (čvrstu granicu) usled neprizmatičnosti korita. a) Sila pritiska u (bilo kom) poprečnom preseku se računa na osnovu pretpostavke o hidrostatičkom rasporedu pritisaka. Za poprečni presek kao na skici će se uvesti lokalni koordinatni sistem sa koordinatom ζ umesto z.
Elementarna sila pritiska, tj. sila pritiska na elementarnu površinu (b dζ ) se dobija množenjem elementarne površine sa pritiskom u težištu površine:
ρ g( h − ζ ) b dζ Sila pritiska na ceo poprečni presek dobija se integraljenjem elementarnih sila po čitavom preseku (uz korišćenje pretpostavke da je ρ = const. po preseku): h
ρ g ∫ ( h − ζ ) b dζ = ρ g I 1 0 h
gde je: I 1 =
∫ ( h − ζ ) b dζ
- statički moment površine preseka u odnosu na nivo slobodne
0
površine. Dejstvo sile pritiska u presecima x1 i x2, tokom intervala t1 do t2 se sada može izraziti kao: t2
g ∫ ( ρ I 1 ) x − ( ρ I 1 ) x dt 1 2 t1
b) Reakcija na dejstvo sile pritiska na omotač (čvrstu granicu) usled neprizmatičnosti će se izvesti tako što će se prvo razmotriti elementarna sila (sila na elementarnu površinu) usled neprizmatičnosti pri istoj dubini: h( x1 ) = h( x 2 ) = h0 .
94
∂b dx predstavlja povećanje širine elementarne površine (b dζ ) ∂x ∂b dx dζ predstavlja priraštaj elementarne površine na istom Onda ∂x
Sa skice se vidi da na rastojanju dx. rastojanju.
Pošto je jasno da se sile pritiska u dva susedna preseka na istu elementarnu površinu pri istoj dubini poništavaju, onda će, ako su elementarne površine različite, pri istoj dubini ukupna elementarna sila pritiska usled neprizmatičnosti biti jednaka sili na priračtaj elementarne površine. Prema tome, reakcija na elementarnu silu pritiska usled neprizmatičnosti pri istoj dubini (čiji je smer suprotan od smera same sile pritiska) se može izraziti kao: ∂ b ρ g( h − ζ ) dx dζ ∂ x h =h 0
gde
∂b ∂b > 0 znači da je sila pozitivna (kao na prethodnoj skici), a < 0 znači da sila deluje ∂x ∂x
suprotno od pozitivnog smera ose x. Reakcija na silu pritiska usled neprizmatičnosti na čvrsti omotač između dva preseka na rastojanju dx dobija se kada se integrale po dubini elementarne sile reakcije (za istu dubinu): ∂ b ρ g ∫ ( h − ζ ) dx dζ = ρ g I 2 dx ∂ x h= h o h
0
∂ b gde je sa I 2 = ∫ (h − ζ ) dζ ∂ x h= h0 o na nivo slobodne površine. h
obeležen statički momenat promene površine u odnosu
Ukupna reakcija na silu pritiska usled neprizmatičnosti na čvrst omotač između preseka x = x1 i x = x2, tokom perioda t1 do t2 iznosi: t 2 x2
g ∫ ∫ ρ I 2 dx dt t1 x1
Očigledno je da izvedeni izraz važi samo za postepenu promenu širine. Dejstvo površinske sile od devijatorskog dela napona na fluid u ∀ , od t1 do t2: 95
t2
∫ ∫σ
d ji
n j dSdt
t1 S
će se izvesti pod pretpostavkom da je u dejstvu devijatorskog dela napona dominantna dejstvo napona trenja na dodiru sa čvrstom granicom. Ova pretpostavka opravdava se time što je devijatorski deo napona propocionalan izvodima brzina, a ti izvodi su najveći uz čvrstu granicu. Sila trenja na dodiru sa čvrstom granicom, za poprečni presek kao na skici:
se može izraziti kao: − ∫ τ d na dO O
gde τdna označava napon trenja na dodiru vode sa čvrstom granicom (dnom), O označava okvašeni obim, a znak minus označava smer suprotan od smera tečenja. Ako se razmatra element okvašenog obima preseka (prikazan na prethodnoj skici) vidi se da je:
cos ψ =
dy dy odnosno dO = dO cos ψ
Stoga se sila trenja na dodiru sa čvrstom granicom za poprečni presek može izraziti kao: B
−∫ 0
τ d na cos ψ
dy
a ukupna sila trenja na dodiru sa čvrstom granicom, na rastojanju od x = x1 do x = x2, tokom perioda t1 do t2 može izraziti kao: t 2 x2 B
−∫ ∫ ∫ t1 x1 0
τ d na cos ψ
dydxdt
Pošto je raspored napona trenja po okvašenom obimu većinom nepoznat, to se u 96
praktičnim zadacima obično koristi srednja (prosečna) vrednost napona trenja za presek: B 1 τ d na ~ τ = ∫ dy O 0 cos ψ
odakle se dobija da je: B
∫ 0
τ d na cos ψ
dy = τ~ O
Ako se iskoristi poslednji izraz, ukupna sila trenja na dodiru sa čvrstom granicom, na rastojanju od x = x1 do x = x2, tokom perioda t1 do t2 se može izraziti kao: t 2 x2
− ∫ ∫ τ~ Odxdt t1 x1
Prosečni (srednji) napon trenja za presek može da se izrazi i preko prosečne brzine za presek:
τ~ = Cτ
ρ u2 2
gde je: u - prosečna (srednja) brzina za presek, a Cτ - koeficijent trenja. Uz pomoć prethodnog izraza za prosečni (srednji) napon trenja za presek, ukupna sila usled napona trenja postaje: t 2 x2
− ∫ ∫ Cτ t1 x1
ρ u2 2
Odxdt
Prema de St. Venant-ovim pretpostavkama, za τ~ se mogu koristiti izrazi koji važe za ustaljeno tečenje. Na primer, iz poznatog Darcy-Weisbach-ovog izraza: E1 − E 2 = Ei z g
L u2 = Cτ R 2g
odnosno: E1 − E 2 E i z g u2 = = S f = Cτ L 2 gR L gde je: R =
A - hidraulički radijus, Sf - nagib linije energije ili tzv. nagib trenja, dobija se da O
je: Cτ
u2 = gRS f 2
odnosno: 97
Cτ
ρ u2
= ρ g RS f
2
Ako se iskoristi prethodni izraz, ukupna sila usled napona trenja postaje: t 2 x2
− ∫ ∫ ρ gRS f Odxdt t1 x1
ili, pošto je R =
A : O t 2 x2
− ∫ ∫ ρ gAS f dxdt t1 x1
Za nagib trenja Sf, u prethodnom izrazu za ukupnu silu usled napona trenja, mogu se koristiti uobičajene empirijske veze srednje brzine i nagiba trenja za ustaljeno tečenje. Na primer, ako se koristi Manning-ova formula: 1 23 u = R Sf n gde je n Manning-ov koeficijent hrapavosti, dobija se: Sf =
u2n2 4
R3
=
Q2n 4
A2 R 3
gde je sa Q obeležen proticaj. Ako se koristi Chezy-jeva formula:
u = C RS f gde je sa C obeležen Chezy-jev koeficijent hrapavosti, dobija se: Sf =
Q2 u2 = 2 2 CR C A R
Kada se uvrste svi prethodno izvedeni članovi, dobija se integralni oblik jednačine održanja količine kretanja za jednodimenzionalno neustaljeno tečenje u otvorenim tokovima: x2
t2
x1
t1
∫ ( ρ u A) t2 − ( ρ u A) t1 dx + ∫ ( ρ u 2 A) x2 − ( ρ u 2 A) x1 dt = t 2 x2
t2
t 2 x2
t 2 x2
t1 x1
t1
t1 x1
t1 x1
= g ∫ ∫ ρ A S 0 dxdt + g ∫ ( ρ I 1 ) x − ( ρ I 1 ) x dt + g ∫ ∫ ρ I 2 dxdt − g ∫ ∫ ρ A S f dxdt 1 2 Uočava se da izvedene integralne jednačine dozvoljavaju postojanje diskontinuiteta (strmih talasa, hidrauličkog skoka itd.) u toku - jer nikakve pretpostavke o kontinualnosti funkcija nisu uvedene tokom izvođenja. Takođe, tokom izvođenja nije zahtevano da ( x 2 − x1 ) ili ( t 2 − t 1 ) budu ograničeni (mali).
98
3.2. Diferencijalni oblik jednačina Diferencijalne jednačine će se izvesti od integralnih, uz pretpostavku da su zavisno promenljive neprekidne i diferencijabilne funkcije. Koristeći razvijanje funkcije u Taylor-ov red može se, na primer, pisati:
( f )x
∂ f + ∂x
= ( f )x
2
1
x1
∂ 2f ( x2 − x1 ) + ∂ x 2
(x
( f )t
= ( f )t
2
1
∂ f + ∂x
t1
∂ 2f ( t 2 − t1 ) + ∂ t 2
(t t1
− x1 )
2
+...
2!
x1
ili:
2
2
− t1 )
2
+...
2!
gde f može biti ( ρ A) ili ( ρ u A) , ( ρ u 2 A) , itd. Ako se u prethodnim izrazima zanemare svi članovi sa izvodima višim od prvog izvoda, onda se sledeći integrali mogu pisati kao:
∫ [( f )
t2
t1
x2
−( f )x
1
]
∂ f dt = ∫ ( f ) x + 1 ∂x t1 t2
t2 ∂ f ( x2 − x1 ) − ( f ) x1 dt = ∫ ∂ x t1 x1
x1
x2 ∂ f ( t 2 − t1 ) + ( f ) t1 dx = ∫ ∂ t x1 t1
t1
(x
2
− x1 ) dt
odnosno:
∫ [( f )
x2
x1
t2
−( f )t
1
]
∂ f dx = ∫ ( f ) t + 1 ∂x x1 x2
(t
2
− t1 ) dt
Kada se pusti da x 2 → x1 i t 2 → t 1 i nađu limesi prethodnih integrala, dobija se:
∫ [( f ) x − ( f ) x ] dt = xlim ∫ ∂ x ( x 2 − x1 ) dt = ∫ ∫ ∂ x dxdt x →x →x t2
t2
lim 2
1
2
t1
1
2
1
∂ f
t1
t 2 x2
x1
∂ f
t1 x1
gde je (videti skicu):
∂f lim x →x ∂ x 2
1
(x 2 − x1 ) = x1
( f )x
2
− ( f ) x1
x 2 − x1
odnosno, dobija se: 99
(x 2 − x1 ) = ∫ ∂ x2
x1
f dx ∂x
∫ [( f )
x2
lim
t 2 → t1
x1
t2
− ( f )t
1
]
∂ f dx = lim ∫ t →t x ∂ t x2
2
1
1
(t
− t1 ) dx =
2
x2 t 2
∂ f
∫∫ ∂ t
dtdx .
x1 t1
t1
Ako se prethodno dobijeni izrazi primene na odgovarajuće članove u integralnim jednačinama, odgovarajući članovi u jednačini održanja mase postaju: x2
x2 t 2
x1
x1 t1
lim ∫ ( ρ A) t − ( ρ A) t dx = t 2 → t1 2 1
∫∫
∂ ( ρ A) dtdx ∂t
t2
t 2 x2
t1
t1 x1
lim ∫ ( ρ u A) x − ( ρ u A) x dt = x2 → x1 2 1
∂ ( ρ u A) dxdt ∂x
∫∫
a odgovarajući članovi u jednačini održanja količine kretanja postaju: x2
x2 t 2
x1
x1 t1
lim ∫ ( ρ u A) t − ( ρ u A) t dx = t 2 → t1 2 1 t2
∫∫
lim ∫ ( ρ u 2 A) x − ( ρ u 2 A) x dt = x2 → x1 2 1 t1
∂ ( ρ u 2 A) ∫ ∫ ∂ x dxdt t x
t 2 x2
1
t2
t 2 x2
t1
t1 x1
lim ∫ ( ρ I 1 ) x − ( ρ I 1 ) x dt = x2 → x1 2 1
∂ ( ρ u A) dtdx ∂t
∫∫
1
∂ ( ρ I1 ) dxdt ∂x
Dobijeni članovi će se uvrstiti u integralne jednačine, pa jednačina održanja mase postaje: t x ∂ ( ρ A) ∂ ( ρ u A) ∫ ∫ ∂ t dtdx + ∫ ∫ ∂ x dxdt = 0 x t t x
x2 t 2
2
2
1 1
1
1
a jednačina održanja količine kretanja postaje: t x ∂ ( ρ u A) ∂ ( ρ u 2 A) ∫ ∫ ∂ t dtdx + ∫ ∫ ∂ x dxdt = x t t x
x2 t 2
2
1 1
1
2
1
t x t x ∂ ( ρ I1 ) = g ∫ ∫ ρ A S 0 dxdt − g ∫ ∫ dxdt + g ∫ ∫ ρ I 2 dxdt − g ∫ ∫ ρ A S f dxdt ∂x t x t x t x t x t 2 x2
1
1
t 2 x2
1
1
2
2
2
2
1
1
1
1
Napisane jednačine važe za proizvoljnu oblast integracije ( x1 , x 2 ) i ( t1 , t 2 ) , pa prema tome važe i u diferencijalnom obliku: jednačina održanja mase:
100
∂ ( ρ A) ∂ ( ρ u A) + =0 ∂t ∂x jednačina održanja količine kretanja: ∂ ( ρ u A ) ∂ ( ρ u 2 A) ∂ ( ρ I1 ) + = −g + g ρ I 2 + g ρ A S0 − S f ∂t ∂x ∂x
(
)
Ako se pretpostavi da je gustina ρ = const. i koristi Q = uA, dobija se sledeći oblik jednačina: jednačina održanja mase:
∂ A ∂Q + =0 ∂t ∂x jednačina količine kretanja: ∂ Q ∂ Q2 + + g I1 = g I 2 + g A S0 − S f ∂ t ∂ x A
(
)
Izvedeni oblik jednačina je tzv. konzervativni ili divergentni oblik jednačina (jednačine su izvedene direktno iz osnovnih zakona održanja u integralnom obliku, a ako je desna strana jednaka nuli sledi da je divergencija mase i količine kretanja po bilo kojoj zatvorenoj konturi u ( x , t ) ravni jednaka nuli - tj. masa i količina kretanja su očuvani). Iako je tokom izvođenja divergentnog oblika diferencijalnih jednačina pretpostavljeno da su zavisno promenljive neprekidne (što isključuje diskontinuitete, strme talase, hidraulički skok i sl.), može se pokazati da divergentni oblik još uvek omogućava tzv. "slabo" (weak) rešenje u okolini diskontinuiteta (tj. obezbeđuje očuvanje mase i količine kretanja u okolini diskontinuiteta). U slučajevima kada se ne očekuje pojava strmih talasa, hidrauličkog skoka i sličnih diskontinuiteta, obično se eliminišu šlanovi sa integralima I1 i I2 , jer je ove integrale teško sračunati. Ako je f = f ( x , y , z, t ) bilo kakva funkcija, onda se na osnovu Leibnitz-ove teoreme može pisati da je:
∂ f ∂ ∂ ∫ ∂ z dz = ∂ x ∫ f dz − f ( x, y, z, t ) ∂ z z
z2
z2
1
1
Ako se prethodno primeni na član:
∂ I1 ∂ h = ∫ ( h − ζ ) b dζ ∂x ∂x0 kao krajnji rezultat se dobija: 101
z2
z x z
1
∂ I1 ∂h =A +I ∂x ∂x 2 Kada se prethodni izraz uvrsti u poslednji napisani oblik jednačine održanja količine kretanja, diferencijalne jednačine postaju: jednačina održanja mase:
∂ A ∂Q + =0 ∂t ∂x dinamička jednačina:
∂ h ∂ Q ∂ Q2 + − S0 + g A S f = 0 + g A ∂ t ∂ x A ∂ x Poslednja jednačina nije više u tzv. divergentnom obliku i ne omogućava očuvanje količine kretanja u slučaju pojave diskontinuiteta, pa se zato za nju koristi izraz dinamička jednačina a ne više jednačina održanja količine kretanja. Dobijene su dve jednačine sa dve zavisno promenljive Q ( x , t ) i h ( x , t ) (dok je npr. površina A = A ( h, x , t ) , itd.). Jednačine se inače mogu formulisti tako da zavisno promeljive budu: a) Q ( x , t ) i h ( x , t ) , ili b) Q ( x , t ) i Z ( x , t ) , ili c) u ( x , t ) i h ( x , t ) , ili d) u ( x , t ) i Z ( x , t ) gde je sa Z obeležen nivo slobodne površine vode. Na primer, ako se sa Z d na obeleži kota dna, a sa S 0 obeleži lokalni nagib dna, onde se može pisati da je:
h = Z − Zd n a S0 = −
∂ Zd n a ∂x
Takođe, sa skice se vidi da je priraštaj površine poprečnog preseka usled priraštaja kote slobodne površine:
102
dA = BdZ
pa je:
dA =B dZ Konačno, pošto je A = A ( Z ) , uz pomoć prethodnog izraza dobija se da je:
∂ A ∂ A∂ Z ∂Z = =B ∂t ∂Z ∂t ∂t Ako se prethodni izrazi uvrste u diferencijalne jednačine, dobija se sledeći oblik jednačina: jednačina kontinuiteta:
∂Z 1∂Q + =0 ∂t B∂x dinamička jednačina:
∂ Q ∂ Q2 ∂Z + + g ASf = 0 + g A ∂ t ∂ x A ∂x Ovo je i najčešće korišćeni oblik diferencijalnih jednačina, i to iz praktičnih razloga. Naime, lokalni pad dna S0 i dubinu h je lako odrediti za kanal, ali ne i za reku, što se vidi i sa skice. Zato je pogodnije koristiti oblik jednačina sa kotom slobodne površine vode Z.
103
3.3. Granice važenja jednodimenzionalnih jednačina. Složeni preseci Pod složenim presecima ovde se, pre svega, misli na složene preseke u prirodnim vodotocima, ali se slični primeri mogu naći i kod kanala, npr. u mešovitoj kanalizaciji. Kod prirodnih vodotoka se razlikuje tzv. glavno korito i inundacija:
Voda većim delom godine teče kroz glavno korito, osim u kišnim periodima kada ispuni i inundacije. Otpor trenja u inundacijama je obično veći od onog u glavnom koritu (zbog vegetacije) Postavlja se pitanje kako formulisati nagib trenja (nagib linije energije) za složeno korito? U poglavlju 3.1. je rečeno da se nagib trenja (nagib linije energije) izražava pomoću istih impirijskih izraza koji se koriste i za ustaljeno tečenje, kao što su npr.: - Manning-ov (odnosno Strickler-ov) izraz: 2 1 Q = AR 3 S f n
( Q = K S t r AR
2 3
Sf )
gde je n (KStr) - Manning-ov (Strickler-ov) koeficijent hrapavosti, - Chezy-jev izraz:
Q = CA RS f gde je C - Chezy-jev koeficijent hrapavosti, - Darcy-Weisbach-ov izraz:
Ei z g
L u2 f L u2 = Cτ = R 2g 4 R 2g
gde je f - Darcy-Weisbach-ov koeficijent hrapavosti. Darcy-Weisbach-ov izraz se može preformulisati kao:
104
Ei z g
u2 u2 = S f = Cτ = f L 2 gR 8 gR
odakle sledi da je:
8 gRS f
u=
f
odnosno:
8 gR Sf f
Q= A
Kao što se vidi, svi prethodni izrazi mogu se napisati u opštem obliku:
Q = K Sf tj., nagib trenja se uopšteno može izraziti kao: Q2 Sf = 2 K gde je, ako se koristi Manning-ov (Strickler-ov) izraz: K=
2 1 AR 3 n
2
( K = K S t r AR 3 )
ako se koristi Chezy-jev izraz: K = CA R a ako se koristi Darcy-Weisbach-ov izraz:
K=A
8 gR f
U prethodnim izrazima K (conveyance) definiše propusnu moć korita, jer u sebi sadrži i koeficijent trenja i geometrijske karakteristike preseka: površinu porečnog preseka A i hidraulički radijus R. Kod složenih korita se propusna moć glavnog korita i inundacija značajno razlikuju. Tečenje u inundaciji je, zbog velikog trenja, znatno sporije od tečenja u glavnom toku. Jedan način da se razlike u propusnoj moći glavnog korita i inundacija uključe u jednačine je da se napravi razlika između širine toka koja se koristi u dinamičkoj jednačini (tzv. Baktivno) i širine toka koja se koristi u jednačini kontinuiteta (tzv Bakumulisanja), kao što je pokazano na skici.
105
Baktivno ili tzv. "aktivna" širina se onda koristi u dinamičkoj jednačini uz pretpostavku da su brzine u inundacijama toliko male da ne daju značajan doprinos dinamičkoj jednačini. Pri tome u dinamičku jednačinu ulazi propusna moć K za "aktivnu" površinu (skica). Bakumulisanja ili "akumulaciona" širina se koristi u jednačini kontinuiteta uz pretpostavku da je dominantna ulogs inundacije da akumuliše višak vode. Drugi (bolji) način da se u jednačinama uzmu u obzir razlike u propusnoj moći glavnog korita i inundacija je da se definiše propusna moć K za složeno korito i cela unese u dinamičku jednačinu. Propusna moć za celo korito definiše se kao: K=
∑K
i
i
gde je, ako se koristi Manning-ova formula: 5 2
1 1 Ai3 K i = Ai Ri3 = n ni 23 O ili, ako se koristi Chezy-jeva formula:
K i = Ai C i Ri = C i
3 2 i 1 2
A O
i gde "i" označava element preseka (kao na skici), pri čemu je zanemareno trenje na dodiru površina A1 i A2, A2 i A3 itd.
Onda se (za ceo presek) može pisati: Q = ∑ Qi =
∑K
i
i
106
i
Sfi
Uobičajeno je, međutim, da se pretpostavi da je pad linije energije konstantan po poprečnom preseku, pa se može pisati da je: Q=
∑K
i
Sf
i
odnosno da se u jednačinama koristi: Q2 Sf = ∑ Ki2 i
Razmotriće se koliko je ova poslednja pretpostavka dobra. Obeležiće se, za trenutak, brzina u jednoj tački preseka sa u, brzina osrednjena po dubini (srednja brzina u jednoj vertikali poprečnog preseka) sa u~ , a srednja brzina za ceo ~ presek sa u~ . Raspored brzina po dubini toka, posebno u razvijenom turbulentnom tečenju, je takav da nema velikog odstupanja brzina u tačkama duž vertkale od u~ :
h
1 gde je: u~ = ∫ udz . h0 Međutim, raspored brzina u~ (tj. brzina osrednjenih po dubini) u pravcu upravno na tok može biti veoma nejednolik u slučaju složenih preseka, kao što se vidi sa skice:
U ovakvom slučaju je pretpostavka da je nagib trenja (nagib linije energije) konstantan po poprečnom preseku u suprotnosti sa jednom od osnovnih de St. Venant-ovih pretpostavki 107
∂E konstantno ∂x ∂ Π ∂ u2 po preseku, onda pošto brzina odnosno varirati po varira po preseku, mora i ∂ x 2g ∂x da je nivo slobodne površine konstantan po preseku (ako se pretpostavi da je
preseku). Ako se pak želi zadržati de St. Venant-ova pretpostavka da je nivo slobodne površine konstantan po preseku, onda se u slučaju neravnomernog rasporeda brzina po širini toka mora menjati i nagib trenja po širini toka, tj. trebalo bi u jednačinama koristiti:
Q = ∑ Qi =
∑K
i
Ovo poslednje je, jednodimenzionalni model.
i
i
međutim,
veoma
Sfi teško
(ili
nemoguće)
uključiti
u
Zato se zadržava pretpostavka da je nagib linije energije ( Sf ) uniforman po preseku, što je tačno jedino ako je raspored brzina uniforman po preseku. U slučaju složenih preseka sa značajnom neravnomernošću brzina po preseku, pretpostavkom da je nagib linije energije ( Sf ) uniforman po preseku se u proračun unosi greška. Ovo su istovremeno i situacije na granici važenja jednodimenzionalnih modela. Da bi se ipak donekle uzeo u obzir uticaj neravnomernosti rasporeda brzina po širini toka, u drugi član u dinamičkoj jednačini se uvodi tzv. koeficijent neravnomernosti rasporeda brzine α. Umesto:
∂ ∂ Q2 ( uQ) = ∂x ∂ x A piše se:
∂ Q2 α ∂ x A gde je:
1
α =~ ~2
~ 1 B u~ = ∫ u ( y ) h ( y ) dy A0
B
∫ u~ 2 ( y ) h ( y ) dy
u A0
U praktičnim zadacima koeficijent α postaje parametar kalibracije: podešavanjem vrednosti koeficijenta α postiže se slaganje merenih i sračunatih vrednosti nivoa slobodne površine vode. U većini slučajeva kalibracijom se dobija da je α ≅ 1 . Ako se dobije da je α = 2 ili 3, to znači da je raspored brzina po širini toka toliko neravnomeran da St. Venant-ove pretpostavke više ne važe - tečenje više nije jednodimenzionalno. Na ovaj način koeficijent α se može koristiti i za procenu granice važenja jednodimenzionalnih jednačina.
108
3.4. Kratak pregled šema metoda konačnih razlika Uobičajeno je da se različite šeme metoda konačnih razlika, koje se žele primeniti na jednačine kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima, prvo primene na pojednostavljene (linearizovane) jednačine tečenja o otvorenim tokovima:
∂h ∂u + h0 = 0 h0 = const . ∂t ∂x ∂u ∂h +g =0 ∂t ∂x Pojednostavljene (linearizovane) jednačine tečenja imaju jasan matematički karakter (jednačine su hiperboličnog tipa), kao što je pokazano na početku poglavlja 1, pa se jasno mogu definisati potrebe za početnim i graničnim uslovima. Takođe, pojednostavljene jednačine se mogu svesti na tzv. “karakteristični” oblik (slično kao što je to urađeno za kompletne de St. Venant-ove jednačine u poglavlju 1.2), pa onda rešiti metodom karakteristika, koja, ako se ne koristi interpolacija, daje tačno rešenje. Konačno, pošto su test jednačine linearne, to se na njih i predloženu šemu metoda konačnih razlika može primeniti Von-Neumann-ova analiza stabilnosti, prikazana u poglavlju 1.4. Tek kada se predložena šema testira na pojednostavljenim jednačinama, prelazi se na njenu primenu na pune jednačine za jednodimenzionalno neustaljeno tečenje u otvorenim tokovima. Pune jednačine su komplikovanije, a uz to i nelinearne, pa se na njih ne može primeniti rigorozna matematička analiza. Stoga se testiranje predložene šeme primenjene na pune jednačine većinom svodi na numeričke eksperimente i poređenje sa rezultatima merenja. Predložen je i koristi se veliki broj šema metoda konačnih razlika za rešavanje jednačina kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima. Ne postoji "najbolja" šema. Među najčešće korišćenim (najpoznatijim) šemama nema značajne razlike u tačnosti. Zato je praktična primenljivost važan parametar u izboru šeme. Šeme metode konačnih razlika za rešavanje jednačina kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima se mogu razvrstati u dve grupe: 1) šeme zasnovane na diferencijalnim jednačinama, koje mogu biti: a) eksplicitne šeme, ili b) implicitne šeme, i 2) šeme zasnovane na integralnim jednačinama.
Šeme zasnovane na diferencijalnim jednačinama Ove šeme su zasnovane na direktnim aproksimacijama izvoda u parcijalnim diferencijalnim jednačinama pomoću konačnih razlika, kao što je pokazano u poglavlju 1.3. Eksplicitne šeme su one kod kojih postoji eksplicitno rešenje za nepoznate u bilo kojoj tački (ne rešava se sistem jednačina za sve tačke).
109
Ključne prednosti eksplicitnih šema sastoje se u tome da su eksplicitne šeme jednostavne za programiranje i da, potencijalno, daju dobru tačnost u slučaju proračuna propagacije strmih talasa. Jedan od osnovnih nedostataka eksplicitnih šema je u tome da su eksplicitne šeme većinom uslovno stabilne. Uobičajeni uslov stabilnosti je:
c∆t gde je: Cr = - Courant-ov broj, a c = ∆x
Cr ≤ 1
gh - brzina propagacije talasa.
Uslov stabilnosti često zahteva mali računski korak po vremenu ∆t, posebno ako problem zahteva malo ∆x, pa je potrebno veliko računarsko vreme da bi se sračunao ukupan period propagacije talasa. Takođe, za šeme kod kojih aproksimacija izvoda po prostoru uključuje tri tačke po prostoru (čime se postiže da šema bude drugog reda), nije sasvim očigledno kako treba računati zavisno promenljive na granicama, jer se jedna od zahtevanih tačaka u aproksimaciji izvoda po prostoru nalazi van razmatrane oblasti strujanja, što predstavlja praktičnu poteškoću u primeni šeme. Ukratko će se prikazati par primera eksplicitnih šema.
Lax-ova šema kao aproksimacije izvoda po prostoru i vremenu koristi sledeće aproksimacije za izvode po prostoru i vremenu:
∂ f f n − f i −n1 ≅ i +1 2 ∆x ∂x ∂ f f n − f 0n ≅ i ∆t ∂t što je prikazano i na skici:
gde je:
f i n1 + f i −n1 f 0n = α f i n + ( 1 − α ) + 2
110
Šema daje tačno rešenje linearizovanih (pojednostavljenih) jednačina tečenja za slučaj kada je α = 0 i
∆x = ∆t
gh0 .
Ukoliko je α = 0 i ∆x = const., Lax-ova šema je šema drugog reda, a ukoliko je α > 0 šema je prvog reda. Ukoliko je α = 1 šema je nestabilna. Ako je α = 1, aproksimacije izvoda po prostoru i vremenu postaju:
∂ f f i +n1 − f i −n1 ≅ ∂x 2 ∆x ∂ f f i n+1 − f i n ≅ ∂t ∆t što je zapravo ona ista šema čija je primena na čistu advekciju analizirana u poglavlju 1.4 (gde je takođe zaključeno da je šema nestabilna).
Leap-frog šema koristi sledeće aproksimacije za izvode po prostoru i vremenu:
∂ f f i +n1 − f i −n1 ≅ 2 ∆x ∂x ∂ f f i n+1 − f i n−1 ≅ 2 ∆t ∂t kao što je prikazano i na skici:
Leap-frog šema je šema drugog reda za ∆x = const., a inače je prvog reda. I ova šema daje tačno rešenje linearizovanih (pojednostavljenih) jednačina tečenja za slučaj kada je ∆x = ∆t
gh0 .
Osim problema sa aproksimacijama izvoda po prostoru na granicama, primena ove šeme stvara i dodatnu praktičnu poteškoću: potrebna su dva početna uslova (dva vremenska nivoa) da se započne proračun.
Implicitne šeme su one kod kojih se ne može dobiti eksplicitno rešenje u jednoj tački, već se mora mora napisati i rešavati sistem jednačina, što je i osnovna poteškoća kod primene implicitnih šema. Osnovna prednost implicitnih šema je u tome što su implicitne šeme 111
većinom bezuslovno stabilne, što dozvoljava da se u proračunima koristi relativno veliki računski korak po vremenu ∆t. Ukratko će se prikazati par primera implicitnih šema. Abbott-Ionescu šema koristi tzv. smaknutu ili pomerenu računsku mrežu (staggared grid), tako da se zavisno promenljive Z i Q ne računaju u istim, već u različitim ili "smaknutim" tačkama, kao na skici:
U jednačini kontinuiteta koriste se aproksimacije izvoda po prostoru koje su centrirane oko tačaka u kojima se računa kota Z (skica), dok su u dinamičkoj jednačini aproksimacije izvoda po prostoru centrirane oko tačaka u kojima se računa protok Q. Abbott-Ionascu šema je šema drugog reda. Šema je bezuslovno stabilna, nema amplitudnu, ali ima faznu grešku, što se manifestuje kroz numeričke oscilacije, kao u slučaju hidrograma na skici:
Stoga je kod primene ove šeme uobičajeno da se u jednačine doda tzv. disipativni član ili član veštaške difuzije, čime se postiže da se "izglade" oscilacije. Šema koristi tri tačke u aproksimacijama izvoda po prostoru, što stvara probleme pri primeni na granicama. Takođe, pošto se zavisno promenljive Z i Q računaju u različitim tačkama, to stvara dodatne praktične poteškoće pri primeni šeme na granicama. Šema ne daje tačno rešenje pojednostavljenih (linearizovanih) jednačina tečenja ni za jednu kombinaciju parametara ∆x , ∆t , gh0 .
Vasiljev šema koristi sledeće aproksimacije izvoda po vremenu i prostoru:
∂ f f i n+1 − f i n ≅ ∂t ∆t ∂ f f i +n1+1 − f i −n1+1 ≅ 2 ∆x ∂x što je prikazano i na skici:
112
Šema koristi tri tačke u aproksimacijama izvoda po prostoru, što stvara probleme pri primeni na granicama. Takođe, šema ima značajnu amplitudnu grešku (tj. prouzrokuje značajnu numeričku difuziju). Ova ista šema, ali primenjena na jednačinu čiste advekcije, analizirana je i u poglavlju 2.2.1 - tamo je to tzv. implicitni polukorak kombinovane implicitno-eksplicitne šeme - gde je takođe zaključeno da šema prouzrokuje značajnu numeričku difuziju.
Šeme zasnovane na integralnim jednačinama Prikazaće se osnovna ideja od koje polaze šeme zasnovane na integralnim jednačinama. Razmatraće se integralni oblik de St. Venant-ovih jednačina. Koristiće se, zbog jednostavnosti, jednačina kontinuiteta. Pretpostaviće se, zbog jednostavnosti, da je ρ = const i koristiće se da je uA = Q, tako da jednačina kontinuiteta u integralnom obliku glasi: x2
∫( A
t2
x1
)
t2
(
)
− At1 dx + ∫ Qx2 − Qx1 dt = 0 t1
Ako se posmatra element računske mreže, kao na skici:
vidi se da je: x2
∫A
t2
dx - integral duž BC, od x1 do x2 u t2 (što predstavlja ukupnu masu, po jedinici gustine,
x1
unutar kontrolne zapremine u tranutku t2), x2
∫ A dx - integral duž AD, od x1 do x2 u t1 (što predstavlja ukupnu masu, po jedinici gustine, t1
x1
unutar kontrolne zapremine u trenutku t1),
113
t2
∫Q
x2
dt - integral duž DC, od t1 do t2 u x2 (što predstavlja ukupan izlaz mase, po jedinici
t1
gustine, tokom razmatranog vremenskog intervala kroz presek u x2), t2
∫Q
x1
dt - integral duž AB, od t1 do t2 u x1 (što predstavlja ukupan ulaz mase, po jedinici
t1
gustine, tokom razmatranog vremenskog intervala kroz presek u x1). Raspored funkcija duž pojedinih segmenata konture ABCD će se aproksimirati koristeći "podatke" u računskim tačkama: duž BC: A( x , t n+1 ) ≅ ψ Ain+1 + ( 1 − ψ ) Ain++11 duž AD: A( x , t n ) ≅ ψ Ain + ( 1 − ψ ) Ain+1 duž DC: Q( xi +1 , t ) ≅ θ Qin++11 + ( 1 − θ ) Qin+1 duž AB: Q( xi , t ) ≅ θ Qin+1 + ( 1 − θ ) Qin gde su 0 ≤ ψ ≤ 1 i 0 ≤ θ ≤ 1 težinski faktori po prostoru i vremenu. U prethodnim izrazima je vrednost određene funkcije duž odgovarajućeg segmenta predstavljena kao prosek vrednosti iz računskih tačaka na krajevima segmenta, pomoću tzv. “težinskih” faktora ψ i θ . Kada se prethodne aproksimacije uvrste u jednačinu kontinuiteta u integralnom obliku, dobija se: x2
∫ [ψ A
n +1 i
+ ( 1 − ψ ) Ain++11 − ψ Ain + ( 1 − ψ ) Ain+1 ] dx
x1 t2
+ ∫ [θ Qin++11 + ( 1 − θ ) Qin+1 − θ Qin+1 + ( 1 − θ ) Qin ] dt = 0 t1
Vrednosti integrala u prethodnoj jednačini se sada mogu jednostavno odrediti. Integracija je trivijalna, pošto se vrednosti funkcija u tačkama mogu izvući ispred integrala. Da bi se krajnji rezultat napisao u opštijem obliku, zameniće se x1 i x2 sa xi i xi+1, odnosno t1 i t2 sa tn i tn+1 (kao što je to i obeleženo na poslednjoj skici). Kada se uvede obeležavanje xi +1 − xi = ∆xi i t n +1 − t n = ∆t n , pa cela jednačina podeli sa ∆xi ∆t n i preurede članovi, dobija se:
Ain+1 − Ain Ain++11 − Ain+1 Qin++11 − Qin+1 Qin+1 − Qin ψ + (1 − ψ ) +θ + (1 − θ ) =0 ∆t ∆t ∆x ∆x Kao što se vidi, pri izvođenju se pošlo od integralne jednačine, a došlo se do algebarske (približne) jednačine koja izgleda kao da je dobijena primenom najopštije šeme metode konačnih razlika na jednačinu kontinuiteta u diferencijalnom obliku:
114
∂ A ∂Q + =0 ∂t ∂x gde se kao “najopštija šema metode konačnih razlika” smatra šema koja koristi uopštene aproksimacije za izvode po vremenu i prostoru kakve su prikazane u poglavlju 1.3. Dva parametra, ψ - težinski faktor po prostoru i θ - tažinski faktor po vremenu, služe da se odrede osrednjene vrednosti izvoda za element računske mreže (kontrolnu zapreminu), kao na skici:
Uočava se da se pri izvođenju pošlo od integralnih jednačina, odnosno da nisu uvođene pretpostavke o postojanju izvoda. Patametri ψ , θ definišu čitavu familiju šema. U principu, šeme zasnovane na integralnim jednačinama su implicitne. U ovom slučaju nepoznate su: Ain+1 , Ain++11 , Qin+1 , Qin++11 odnosno: Zin+1 , Zin++11 , Qin+1 , Qin++11 jer je A = A(Z). Šeme zasnovane na integralnim jednačinama imaju niz povoljnih osobina sa stanovišta praktične primene. Kao što se vidi, približne (aproksimirane, algebarske) jednačine uključuju samo dve tačke u prostoru, pa nema prektičnih poteškoća pri primeni na granicama, kao kod pojedinih šema zasnovanih na diferencijalnim jednačinama koje uključuju tri tačke u prostoru. Pošto se jednačine pišu za kontrolnu zapreminu (računsku deonicu od x1 do x2) a ne za računsku tačku, onda je postupak rešavanja za svaku računsku deonicu je (pa i za graničnu deonicu) isti. Takođe, susedne kontrolne zapremine mogu imati različito ∆x, a da to ne utiče na aproksimacije. Obe zavisno promenljive se računaju u istim tačkama, što je takođe pogodno sa stanovišta primene graničnih uslova. Konačno, približne jednačine zadržavaju isto fizičko značenje kao i polazne integralne jednačine, tj. njima se izražavaju zakoni održanja (mase, količine kretanja) unutar kontrolne zapremine (elementa računske mreže), što omogućava i jasnije fizičko tumačenje rezultata. Najpoznatija šema metode konačnih razlika zasnovana na integralnim jednačinama (istovremeno i jedna od uopšte najčešće korišćenih šema metode konačnih razlika) je tzv. Preissmann-ova šema ili, kako se još naziva šema 4 tačke, box šema itd.
115
3.5. Preissmann-ova šema Preissmann-ova šema je zasnovana na integralnim jednačinama, napisanim za kontrolnu zapreminu (element računske mreže). Šema se naziva još i šema četiri tačke jer četiri računske tačke, kao na skici, definišu kontrolnu zapreminu (element računske mreže, “računsku ćeliju”).
Preissmann-ova šema predlaže sledeće aproksimacije (približne vrednosti) izvoda funkcije f ( x , t ) po prostoru i vremenu, kao i aproksimaciju same funkcije f ( x , t ) :
∂ f f i +n1+1 − f i n+1 ≅θ + (1 − θ ∆x ∂x
)
f i +n1 − f i n ∆x
∂ f f i n +1 − f i n f i +n1+1 − f i +n1 ≅ψ + (1 − ψ ) ∆t ∆t ∂t f i +n1+1 + f i n +1 f i +n1 + f i n f ≅θ + (1 − θ ) 2 2 pri čemu se, kako je to pokazano na kraju prethodnog poglavlja, predložene aproksimacije funkcije i njihenih izvoda po prostoru i vremenu odnose na računsku ćeliju, a ne na računsku tačku. Funkcija f ( x , t ) može biti bilo koja funkcija koja se javlja u jednačinama (Q, Z, A, 1 K,....). Težinski faktor po vremenu ima vrednosti θ ≥ , dok težinski faktor po prostoru po 2 1 originalnoj šemi ima vrednost ψ = . Pogodnije je, međutim, da se pri izvođenju približnih 2 jednačina težinski faktori ψ i θ ne zamenjuju konkretnim vrednostima, već da se ostane pri opštim oznakama, jer to omogućava veću fleksibilnost pri kasnijem korišćenju šeme.
1 Preissmann-ova šema je šema prvog reda, osim ako je θ ≡ ψ ≡ , kada je šema 2 drugog reda. Šema daje tačno rešenje za pojednostavljene linearizovane jednačine tečenja (za posebnu kombinaciju ∆x i ∆t). De St. Venant-ove jednačine: jednačina kontinuiteta:
∂ A ∂Q + =0 ∂t ∂x 116
i dinamička jednačina:
∂ Q ∂ Q2 ∂Z + + g ASf = 0 α + gA ∂ t ∂ x A ∂x Q2 ) čine sistem od dve nelinearne parcijalne K2 diferencijalne jednačine (nelinearnost se očigledno vidi kod dinamičke jednačine). Kada se navedene aproksimacije primene na funkcije i njihove izvode u de St. Venant-ovim jednačinama dobija se sistem od dve nelinearne algebarske jednačine koje se mogu linearizovati i rešiti primenom Newton-Raphson-ovog iterativnog postupka. Primena Preissmann-ovih aproksimacija, kao i linearizacija tako dobijene nelinearne algebarske jednačine će se detaljno prikazati na primeru jednačine kontinuiteta, jer je jednačina kontinuiteta relativno jednostavna. Primena Preissmann-ovih aproksimacija na dinamičku jednačinu (i linearizacija tako dobijene nelinearne algebarske jednačine) se neće detaljno prikazivati, jer bi to odnelo previše prostora, već će se samo izneti sugestije za pojedine članove jednačine. (gde se član trenja može pisati kao S f =
Jednačina kontinuiteta Kada se na jednačinu kontinuiteta:
∂ A ∂Q + =0 ∂t ∂x primene Preissmann-ove aproksimacije, dobija se sledeća približna (algebarska) jednačina:
Ain+1 − Ain Ain++11 − Ain+1 Qin++11 − Qin+1 Qin+1 − Qin ψ + (1 − ψ ) +θ + (1 − θ ) = 0 (*) ∆t ∆t ∆x ∆x Za linearizaciju i rešavanje dobijene nelinearne algebarske jednačine će se koristiti Newton-Raphson-ov iterativni postupak. Leva strana približne jednačine kontinuiteta (*) je funkcija četiri nepoznate (Qin+1 , Qin++11 , Ain+1 , Ain++11 ), pa se približna jednačina kontinuiteta može uopšteno pisati kao: F ( Qin+1 , Qin++11 , Ain+1 , Ain++11 ) = 0 Iterativni postupak podrazumeva postepeno približavanja tačnom rešenju dok se konačno u poslednjoj iteraciji ne dobije rešenje koje zadovoljava jednačinu. Razmatraće se dve uzastopne iteracije ( m) i ( m + 1) (indeks “m” će se koristiti da obeleži iteraciju). Sa m+1Qin+1 , m+1Qin++11 , m+1Ain+1 , m+1Ain++11 će se obeležiti vrednosti nepoznatih u “sledećoj” ili ( m + 1) iteraciji za koje se traži da zadovolje približnu jednačinu kontinuiteta, tako da je: F ( m+1Qin+1 , m+1Qin++11 , m+1Ain+1 , m+1Ain++11 ) = 0
117
ili, skraćeno: m+ 1
F=0
Sa m Qin+1 ,m Qin++11 , mAin+1 ,m Ain++11 će se obeležiti vrednosti nepoznatih u prethodnoj ili ( m) iteraciji koje još uvek ne zadovoljavaju približnu jednačinu kontinuiteta, tako da je: F ( m Qin+1 , m Qin++11 , m Ain+1 , m Ain++11 ) ≠ 0 ili skraćeno: m
F≠0
Vrednost funkcije F u sledećoj (m+1) iteraciji može se izraziti razvijanjem funkcije u Taylor-ov red u okolini vrednosti iz prethodne (m) iteracije, uz zanemarenje izvoda višeg reda od prvog: m+1
F = mF +
∂F n +1 i m
∂Q
∆Qi +
∂F n +1 i +1 m
∂Q
∆Qi +1 +
∂F n +1 i m
∂ A
∆Ai +
∂F ∂ Ain++11
∆Ai +1 m
odnosno, pošto je A = A( Z ) : m+1
F = mF +
∂F ∂ Qin+1
∆Qi + m
∂F ∂ Qin++11
∆Qi +1 + m
∂F ∂ Ain+1
∂ Ain+1 ∂ F ∂ Ain++11 ∆ Z + ∆Z i +1 i n +1 n +1 n +1 Z A Z ∂ ∂ ∂ i + 1 i i + 1 m m m m
gde su ∆Q i ∆Z priraštaji (promene) Q i Z između dve uzastopne iteracije, npr:
∆Qi =
m+1
Qin+1 − mQin+1
∆Z i =
m +1
Z in+1 − mZ in+1
Pošto se zahteva da vrednost funkcije F određena sa vrednostima nepoznatih iz sledeće (m+1) iteracije bude jednaka nuli, dobija se:
∂F ∂ Qin+1
∆Qi + m
∂F ∂ Qin++11
∆Qi +1 + m
∂F ∂ Ain+1
∂ Ain+1 ∂ F ∂ Ain++11 ∆ Z + ∆Z i +1 =− mF i n +1 n +1 n +1 Z A Z ∂ ∂ ∂ i +1 m i i +1 m m m
što je zapravo linearizovana približna (algebarska) jednačina kontinuiteta, sa nepoznatima ∆Q i ∆Z. Izvodi funkcije na levoj strani prethodne jednačine , kao i funkcija na desnoj strani jednačine određeni su sa vrednostima promenljivih poznatim ili iz poznate prethodne iteracije (m) tekućeg vremenskog nivoa (n+1) ili iz poznatog prethodnog vremenskog nivoa (n), tj. to su poznati koeficijenti u jednačini. Prema tome, ostaje da se odrede izvodi:
∂F n +1 i m
∂Q
=−
118
θ ∆x
∂F
=
n +1 i +1 m
∂Q
∂F n +1 i m
∂ A ∂F
n +1 i +1 m
∂ A
=
=
θ ∆x
ψ ∆t
1− ψ ∆t
∂ Ain+1 m n+1 = bi ∂ Z in+1 m ∂ Ain++11 = m bin++11 n +1 ∂ Z i +1 m Poslednja dva izvoda su određeni na osnovu uobičajene aproksimacije. Pošto je (videti skicu) dA priraštaj površine poprečnog preseka (približno) dA = bdZ , to se može pisati da je = b. dZ
Kada se uvrste dobijeni izvodi, linearizovana jednačina postaje: 1 − ψ m n +1 bi +1 ∆Z i +1 = ∆x ∆x ∆t ∆t m n +1 m n +1 m n +1 m n +1 Ai +1 − Ain+1 Qi +1 − Qi Qin+1 − Qin Ai − Ain = −ψ − (1 − ψ ) −θ − (1 − θ ) ∆t ∆t ∆x ∆x
−
θ
∆Qi +
θ
∆Qi +1 +
ψ
m n +1 i
b
∆Z i +
što se može pisati kao: A∆Zi +1 + B∆Qi +1 = C∆Z i + D∆Qi + G gde su ∆Z i +1 , ∆Qi +1 , ∆Z i , ∆Qi nepoznate, a:
A=
1 − ψ m n+1 bi +1 ∆t B=
θ ∆x
119
C=−
ψ ∆t
D= m
G = −ψ
Ain+1 − Ain − (1 − ψ ∆t
)
m
m n +1 i
b
θ ∆x
Ain++11 − Ain+1 −θ ∆t
m
Qin++11 − mQin+1 Qin+1 − Qin − (1 − θ ) ∆x ∆x
poznati koeficijenti čije se vrednosti mogu sračunati iz vrednosti zavisno promenljivih poznatih ili iz poznate prethodne iteracije (m) tekućeg vremenskog nivoa (n+1) ili iz poznatog prethodnog vremenskog nivoa (n).
Dinamička jednačina Drugi član u dinamičkoj jednačini će se napisati kao: Q2 ∂ A ∂ Q2 Q∂ Q −α 2 α = 2α ∂ x A A∂x A ∂x tako da dinamička jednačina postaje:
∂Z ∂Q Q∂ Q Q2 ∂ A + 2α −α 2 +gA + g ASf = 0 ∂t A∂x ∂x A ∂x U nastavku će se izneti sugestije kako da se primene Preissmann-ove aproksimacije na svaki član dinamičke jednačine.
∂Q , je jednostavan, pa se može direktno primeniti Preissmann∂t ∂ f ova aproksimacija za izvod po vremenu . ∂t Prvi član jednačine,
∂ f Q∂Q , može da se piše kao g , gde g i f označavaju odgovarajuće A∂x ∂x funkcije, pa se mogu koristiti Preissmann-ove aproksimacije: Drugi član, 2α
gin++11 + g in+1 g in+1 + gin g ≅θ + (1 − θ ) 2 2
∂ f f i +n1+1 − f i n+1 ≅θ + (1 − θ ∆x ∂x
)
f i +n1 − f i n ∆x
∂ f Q2 ∂ A Treći član, α 2 , se takođe može pisati kao g , gde g i f označavaju ∂x A ∂x odgovarajuće funkcije, pa se mogu koristiti iste Preissmann-ove aproksimacije kao u drugom članu. Ovde, međutim, treba uočiti da se: 120
Q2 2 2 = u A može aproksimirati na dva načina, ili kao kvadrat srednje brzine: 2 ui + ui +1 2 u = 2 ili kao srednja vrednost kvadrata brzina: u2 =
ui2 + ui2+1 2
Tek nakon detaljnije analize (jer to nije unapred očigledno) može se pokazati da je prva od navedene dve aproksimacije bolja, jer daje tačnu energetsku jednačinu za ustaljeno tečenje, dok je druga lošija jer za ustaljeno tečenje i naglu promenu preseka veštački "proizvodi" energiju. Zato je bolje koristiti: 2
Q 2 θ Qin++11 Qin+1 1 − θ Qin+1 Qin ≅ + + + 4 Ain+1 Ain A 2 4 Ain++11 Ain+1 Četrvrti član, g A piše kao g
2
∂Z , je takođe jednostavan, jer je g = const., pa takođe može da se ∂x
∂ f , gde g i f označavaju odgovarajuće funkcije, tako da se mogu direktno ∂x
primeniti Preissmann-ove aproksimacije, kao u drugom članu.
Peti član, ili član trenja gASf, se može aproksimirati na više načina, pa će mu se ovde posvetiti i više pažnje. Pre svega, nagib trenja: Q2 Sf = 2 K se zapravo mora pisati kao: Sf =
Q Q K2
jer jedino tako član tranja može imati odgovarajući znak u zavisnosti od smera tečenja. Nadalje, pošto je dinamička jednačina (kao i jednačina kontinuiteta) napisana za deonicu između računskih tačaka ( i ) i ( i + 1) , to Sf predstavlja nagib energetske linije za čitavu deonicu. Istovremeno, Q i K su definisani samo u računskim tačkama (presecima) ( i ) i
( i + 1) . Znači, mora se koristiti neka vrsta osrednjenih Q i K u ( i ) i ( i + 1) , da se izrazi Sf za deonicu
121
Još uvek ne postoji odgovor na pitanje koji način osrednjavanja je najbolji (najbolje odslikava fizičko značenje člana trenja i sl.), pa se u modelima mogu sresti različiti načini osrednjavanja. Međutim, treba imati u vidu da različiti načini osrednjavanja mogu dati značajno različite rezultate. Da bi se ilustrovalo prethodno rečeno, razmatraće se, kao primer, ustaljeno tečenje, gde je: Qi +1 = Qi = Q Posmatraće se četiri uobičajena načina osrednjavanja: Sf1 =
2
1) osrednjeno K :
2) osrednjeno Sf :
Sf2
Q2
β Ki2 + ( 1 − β ) Ki2+1
Q2 =β + (1 − β Ki2 2
3) "otežana" geometrijska sredina K :
4) osrednjeno K:
Sf 4 =
S f3 =
)
Q2 Ki2+1 Q2
(K ) ⋅ (K )(
1− β ) 2 i +1
2 β i
Q2
[β
Ki + ( 1 − β ) Ki +1 ]
2
gde je β - težinski faktor. Ako je Ki ≅ Ki +1 , onda sva osrednjavanja daju približno isto Sf . Međutim, ako je Ki ≠ Ki +1 , onda za usvojeno β, vrednost Sf puno zavisi od načina osrednjavanja. Na primer, 1 za Ki = 2 Ki +1 i β = , dobija se da je S f 2 ≅ 1,5S f 1 , dok S f 3 i S f 4 daju vrednosti između prve 2 dve. Zbog toga u praktičnim primenama β postaje parametar koji se kalibriše (podešavanjem vrednosti koeficijenta β postiže se slaganje merenih i sračunatih rezultata), tako da korisnici programa naviknu da kalibrišu β zavisno od vrste osrednjavanja u programu. Ovde se predlaže da se koristi osrednjeno Sf , tj.: S f = β S f i + ( 1 − β ) S f i +1 = β
Qi Qi K
2 i
+ (1 − β )
Qi +1 Qi +1 Ki2+1
odnosno:
Q n+1 Q n+1 Qin++11 Qin++11 i i Sf =θ β + (1 − β ) ( Kin+1 ) 2 ( Kin++11 ) 2
+ (1 − θ
Qn Qn Qin+1 Qin+1 i i ) β n 2 + (1 − β ) n 2 ( Ki ) ( Ki +1 )
Takođe, član gASf sadrži i A, pa treba i na A primeniti Preissmann-ovu aproksimaciju i pomnožiti sa aproksimiranim Sf.
122
Nakon primene Ppreissmann-ovih aproksimacija, dinamička jednačina postaje algebarska nelinearna jednačina, koja se uopšteno može pisati kao: F ( Qin+1 , Qin++11 , Z in+1 , Zin++11 , Ain+1 , Ain++11 , Kin+1 , Kin++11 ) = 0
Kada se primeni Newton-Raphson-ov iterativni postupak, dobija se:
∂F ∂ Qin+1 +
∆Qi + m
∂F ∂ Ain++11
∂F ∂ Qin++11
n +1 i +1 n +1 i +1 m
∂ A ∂Z m
∆Qi +1 + m
∆Z i +1 +
∂F ∂F ∂ F ∂ Ain+1 ∆ + ∆ + ∆Z Z Z ∂ Z in+1 m i ∂ Z in++11 m i +1 ∂ Ain+1 m ∂ Z in+1 m i
∂F ∂ Kin+1
∂ Kin+1 ∂ F ∂ Kin++11 ∆ + ∆Z i +1 = Z n +1 i n +1 n +1 ∂ ∂ ∂ Z K Z i i i + + 1 1 m m m m
= − F ( m Qin+1 , m Qin++11 , m Z in+1 , m Z in++11 , m Ain+1 , m Ain++11 , m Kin+1 , m Kin++11 ) odnosno: ∂F ∂ Z in++11
+ m
∂F n +1 i +1
∂A
∂ Ain++11 ∂ Z in++11 m
+ m
∂F ∂K
n +1 i +1
∂ K in++11 ∂ Z in++11 m
∆Z i +1 + ∂ F ∂ Qin++11 m
∆Qi +1 = m
∂F ∂F ∂ K in +1 ∂ F ∂ Ain +1 ∂F =− + + ∆ − Z i n +1 n +1 n +1 n +1 ∂ Z in +1 ∂ Qin +1 A Z K Z ∂ ∂ ∂ ∂ i i i i m m m m m m m n +1 m n +1 m n +1 m n +1 m n +1 m n +1 m n +1 m n +1 − F Qi , Qi +1 , Z i , Z i +1 , Ai , Ai +1 , K i , K i +1
(
)
∆Qi − m
Poslednja jednačina je zapravo linearizovana približna (algebarska) dinamička jednačina, napisana u obliku: A′∆Zi +1 + B ′∆Qi +1 = C ′∆Z i + D′∆Qi + G ′ gde su A′, B ′, C ′, D′, G ′ poznati koeficijenti.
Praktične napomene Može se usvojiti da je vrednost funkcije f u prvoj (početnoj, nultoj) iteraciji tekućeg vremenskog nivoa (n+1) jednaka vrednosti funkcije f iz prethodnog vremenskog nivoa (n): 0
f
n +1
= f
n
gde je f bilo koja zavisno promenljiva. Uz ovu pretpostavku, koeficijent G u linearizovanoj približnoj jednačini kontinuiteta postaje: G = −θ
Qin+1 − Qin Qin+1 − Qin Qin+1 − Qin − (1 − θ ) =− ∆x ∆x ∆x
123
Očigledno je da se za ustaljeno tečenje koeficijent G takođe svodi na prethodni oblik, pa pošto je za ustaljeno tečenje ∆Z = ∆Q = 0 (jer nema promene po vremenu, pa ni između dve iteracije), onda se linearizovana približna jednačina kontinuiteta svodi na: G=−
Qin+1 − Qin =0 ∆x
t.j. Qin+1 = Qin čime je zapravo provereno da li se linearizovana približna (algebarska) jednačina kontinuiteta svodi na odgovarajući oblik u slučaju ustaljenog tečenja. Slična provera se može napraviti i za linearizovanu približnu (algebarsku) dinamičku jednačinu, koja u slučaju ustaljenog tečenja treba da se svede na Bernoulli-jevu jednačinu. Još jedna uobičajena provera da li su linearizovane približne (algebarske) jednačine dobro izvedene je provera dimenzija članova jednačina. Poslednja napomena odnosi se na određivanje geometrijskih karakteristika poprečnog preseka struje. Za svaku poznatu kotu Z (bilo da je to Z n poznato iz prethodnog vremenskog nivoa, ili m Z n+1 poznato iz prethodne iteracije tekućeg vremenskog nivoa) mogu se uvek ∂ A odrediti odgovarajuće geometrijske karakteristike preseka (širina b = , površina A, ∂ Z ∂K ). Kod jednostavnih preseka okvašeni obim O, hidraulički radijus R, propusnost K, izvod ∂Z (trougaoni, trapezni, itd., lao na skici) to je moguće uraditi i analitički, dok se kod proizvoljnih preseka prirodnih vodotoka to uglavnom radi grafički, tako što se presek izdeli na niz elemenata (kao na skici), pa se geometrijske karakteristike kao što su površina A ili okvašeni obim O, za određenu kotu Z sračunavaju kumulativno.
124
3.5.1. Primena Preissmann-ove šeme Pošto su parcijalne diferencijalne jednačine kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovoma (jednačina kontinuiteta i dinamička jednačina) diskretizovane primenom Preissmann-ovih aproksimacija, i pošto su tako dobijene nelinearne algebarske jednačine linearizovane primenom Newton-Raphson-ovog iterativnog postupka, problem je sveden na rešavanje sistema linearnih algebarskih jednačina u svakoj iteraciji svakog vremenskog koraka. Ako se razmatrana jednodimenzionalna oblast tečenja podeli na I računskih tačaka, kao na skici:
pri čemu dve susedne računske tačke definišu tzv. računsku deonicu, onda se za svaku računsku deonicu mogu napisati dve jednačine: linearizovana približna (algebarska) jednačina kontinuiteta: A∆z i +1 + B∆Qi +1 = C∆z i + D∆Qi + G i linearizovana približna (algebarska) dinamička jednačina: A′∆z i +1 + B' ∆Qi +1 = C ′∆z i + D ′∆Qi + G ′ Pošto je ukupan broj tačunskih deonica (I-1), to se ukupno može napisati 2(I-1) linearizovanih algebarskih jednačina za sve računske deonice. Pošto su nepoznate ∆Q i ∆Z u svakoj računskoj tački, to je ukupan broj nepoznatih u napisanim jednačinama 2I. Znači potrebna su dva granična uslova da bi se "zatvorio" sistem jednačina. Jednačine graničnih uslova, takođe u linearizovanom algebarskom obliku, se uopšteno mogu napisati kao:
α ∆ Z + β ∆Q = γ što je zapravo uopšten izraz za vezu Q i Z, tj. ∆Q i ∆Z, ali se pod ovaj opšti oblik mogu podvesti i jednostavniji granični uslovi, kao što su poznati hidrogram Q(t) ili poznati nivogram Z(t). Kako se pojedini tipovi granični uslova svode na napisani uopšteni oblik se detaljno razmatra u poglavlju 3.5.2. U ovom poglavlju se isključivo razmatra postupak za rešavanje sistema algebarskih jednačina. Nadalje će se razmatrati mirno tečenje, što znači da je potreban po jedan granični uslov na uzvodnom i nizvodnom kraju (postupak za rešavanje burnog tečenja je principijelno isti, mada bi algoritam bio nešto drukčiji, jer se kod burnog tečenja oba granična uslova zadaju na uzvodnoj granici). 125
Napisaće se kompletan sistem algebarskih jednačina, zajedno sa jednačinama graničnih uslova. U napisanom sistemu indeks uz koeficijente A, B , C , D, G i A′, B ′, C ′, D ′, G ′ odnosi se na računske deonice (jer su i jednačine pisane za računske deonice), dok se indeks uz nepoznate ∆Q i ∆Z odnosi (kao i same nepoznate) na računske tačke. U jednačinama uzvodnog i nizvodnog graničnog uslova, uz odgovarajuće koeficijente koriste se indeksi uz i niz. Jednačine su obeležene brojevima, tako da je jednačina (1) zapravo linearizovana približna jednačina kontinuiteta za prvu računsku deonicu, a jednačina (1’) linearizovana približna dinamička jednačina takođe za prvu računsku deonicu, itd.
α u z ∆Z1 + β u z ∆Q1
= γ u z (0) = G1 (1)
− C1 ∆Z1 − D1 ∆Q1 + A1 ∆Z 2 + B1 ∆Q2 − C ′ ∆Z − D ′ ∆Q + A ′ ∆Z + B ′ ∆Q 1
1
1
1
1
2
1
2
− C2 ∆Z 2 − D2 ∆Q2 + A2 ∆Z 3 + B2 ∆Q3 − C2 ′ ∆Z 2 − D2 ′ ∆Q2 + A2 ′ ∆Z 3 + B2 ′ ∆Q3 . . . − Ci ∆Zi − Di ∆Qi + Ai ∆Zi +1 + Bi ∆Qi +1 − C ′ ∆Z − D ′ ∆Q + A ′ ∆Z + B ′ ∆Q i
i
i
i
i
i +1
i
i +1
= G1′ (1’) = G2 (2) = G2′ (2’)
= Gi
(i)
= Gi′
(i’)
. . . − C I −1 ∆Z I −1 − DI −1 ∆QI −1 + AI −1 ∆Z I + B I −1 ∆QI = G I −1 (I-1) − C I −1′ ∆Z I −1 − DI −1′ ∆QI −1 + AI −1′ ∆Z I + B I −1′ ∆QI = G I′−1 (I-1)’ α n i z ∆Z I + β n i z ∆QI = γ n i z (I) Napisani sistem algebarskih jednačina se, naravno, može rešiti tako što bi se direktno invertovala matrica koeficijenata, ali je sa stanovišta rada računara, direktno invertovanje matrice jedna od “najskupljih” operacija. Računarsko vreme potrebno za direktno invertovanje matrice koeficijenata proporcionalno je sa I3. U ovom slučaju se, međutim, može uočiti da je matrica koeficijenata tzv. trakasta dijagonalna matrica, pa se i ovde može koristiti Thomas-ov (ili double-sweep) algoritam, čije je potrebno računarsko vreme proporcionalno samo sa I.
Thomas-ov (“double-sweep”) algoritam Jednačina uzvodnog graničnog uslova:
α u z ∆Z1 + β u z ∆Q1 = γ u z se može napisati kao:
∆Q1 = E1 ∆z1 + F1 gde su: 126
(a)
E1 = −
α uz γ , F1 = uz β uz β uz
Takođe, iz jednačina ( 1) i ( 1′ ) : − C1 ∆Z1 − D1 ∆Q1 + A1 ∆Z 2 + B1 ∆Q2 = G1
(1)
− C1′ ∆Z1 − D1′ ∆Q1 + A1′ ∆Z 2 + B1′ ∆Q2 = G1′
(1’)
se može eliminisati ∆Qi, tako da se dobije: ∆Z1 = L1 ∆Z 2 + M 1 ∆Q2 + N 1
(b)
gde su: A1 D1′ − A1′ D1 L1 = C D ′ − C ′D 1
M1 =
1
1
1
B1 D1′ − B1′ D1 C D ′ − C ′D 1
1
1
1
D1G1′ − D1′ G1 N1 = C D ′ − C ′D 1
1
1
1
Ako se izrazi (a) i (b) uvreste u jednačinu (1), dobija se: ∆Q2 = E 2 ∆Z 2 + F2
gde su:
E2 =
F2 =
L1 ( C1 + D1 E1 ) − A1 B1 − M 1 ( C1 + D1 E1 )
N 1 ( C1 + D1 E1 ) + D1 F1 + G1 B1 − M 1 ( C1 + D1 E1 )
Uopšteno, može se pretpostaviti da postoji veza: ∆Qi = E i ∆Zi + Fi (*)
Ako se zatim eliminiše ∆Qi iz ( i ) i ( i ′ ) dobija se: ∆zi = Li ∆Zi +1 + M i ∆Qi +1 + N i (**)
gde su: 127
Ai Di ′ − Ai ′ Di Li = C D ′ − C ′D i
Mi =
i
i
i
Bi Di ′ − Bi ′ Di C D ′ − C ′D i
i
i
i
Di Gi ′ − Di ′ Gi Ni = C D ′ − C ′D i
i
i
i
Zamenom izraza (*) i (**) u jednačinu (i), dobija se: ∆Qi +1 = E i +1 ∆Zi +1 + Fi +1
gde su:
E i +1 =
Fi +1 =
Li ( Ci + Di E i ) − Ai Bi − M i ( Ci + Di E i )
N i ( Ci + Di E i ) + Di Fi + Gi Bi − M i ( Ci + Di E i )
Kao i u poglavlju 2.3, uočava se rekurzivni karakter izraza za E, F: da bi se sračunale vrednosti koeficijenata za narednu tačku, potrebno je znati vrednosti koeficijenata za prethodnu tačku u nizu. Thomas-ov (“double-sweep”) algoritam, koji je zasnovan na rekurzivnom karakteru izraza za E i F, sastoji se iz dve faze: 1) proračun unapred (tzv. “forward sweep”), pri čemu se računaju koeficijenti L, M , N , E , F i 2) proračun unazad (tzv. “backward sweep”) pri čemu se: a) koristi jednačina (**) da se sračuna ∆Zi za poznate ∆Zi+1, ∆Qi+1, i b) koristi jednačina (*) da se sračuna ∆Qi za poznato ∆Zi. Pošto su Li, Mi, Ni koeficijenti koji omogućavaju da se sračuna ∆Z na jednom kraju deonice ako su poznati ∆Z, ∆Q na drugom kraju deonice, može se reći da ovi koeficijenti određuju uticaj računske deonice na rešenje. Koeficijenti Ei, Fi omogućavaju da se sračuna ∆Q za poznato ∆Z u istoj računskoj tački, pa se može reći da ovi koeficijenti određuju uticaj računske tačke na rešenje. Postupak rešavanja za jednu iteraciju sastoji se iz sledećih koraka: 1) Primeniti uzvodni granišni uslov da se sračuna E1, F1. 2) Sračunati koefisijente A,...G, i A′, ... , G′ za računsku deonicu 1, između računskih tačaka 1 i 2. Nema potrebe da se ovi koeficijenti obeležavaju indeksima, jer će se odmah upotrebiti za računanje koeficijenata L,M,N,E,F. 128
3) Sračunati L1, M1, N1 za deonicu 1. 4) Sračunati E2, F2 . 5) Ponoviti korake 2 - 4 za svaku deonicu i, sve do I-1. 6) Kraj proračuna unapred. Kao poslednja jednačina proračuna unapred obijeno je: ∆QI = E I ∆Z I + FI
7) Kombinovati jednačinu nizvodnog graničnog uslova:
α n i z ∆Z I + β n i z ∆QI = γ n i z sa poslednjom jednačinom dobijenom iz proračuna unapred, tako da se sračuna ∆ZI, ∆QI. 8) Sračunati: ∆Z I −1 = LI −1 ∆Z I + M I −1 ∆QI + N I −1
9) Sračunati: ∆QI −1 = E I −1 ∆Z I −1 + FI −1
10) Ponoviti korake 8 - 9 za sve tačke i = I-2, ... , 1. 11) Sračunati: m+1
m+1
Zin+1 = mZin+1 + ∆Zi
Qin+1 = mQin+1 + ∆Qi
koji će se u sledećoj iteraciji koristiti kao poznate vrednosti iz prethodne iteracije. Ovim je završen proračun unazad. Tokom svakog vremenskog koraka prethodno opisani postupak se ponavlja za sledeću iteraciju itd., sve dok promene proticaja i kote između dve iteracije ∆Q i ∆Z u svim računskim tačkama ne budu manji od unapred zadatih malih vrednosti ε ∆Q i ε ∆Z . Jasno je da prilikom sračunavanja koeficijenata L,M,N odnosno E,F imenilac u odgovarajućim izrazima mora biti različit od nule. Pokazuje se da je, ako tečenje nije blizu kritičnog, uvek:
( CD′ ≠ C′D )
CD ′ − C ′D ≠ 0 B − M i ( C + DE i ) ≠ 0
( B ≠ M ( C + DE ) )
129
i
i
tako da pri sračunavanju koeficijenata Li, Mi, Ni, Ei, Fi ne dolazi do deljenja sa nulom. Ali ako tečenje teži kritičnom, može se javiti problem deljenja sa nulom, pa za takav slučaj treba razviti poseban postupak rešavanja.
3.5.2. Granični uslovi U poglavlju 1.2., gde je izučavana metoda karakteristika, kao primer su korišćene jednačine za jednodimenzionalno neustaljeno tečenje u otvorenim tokovima. Tamo je zaključeno da: a) Mirno tečenje (Fr < 1) zahteva po jedan granični uslov na uzvodnom i nizvodnom kraju, kao što se vidi sa skice.
b) Burno tečenje (Fr > 1) zahteva oba granična uslova na uzvodnom kraju (skica):
Kao što je pokazano u prethodnom poglavlju, u slučaju mirnog tečenja postupak rešavanja zahteva i proračun unapred (tzv. “forward sweep”) i proračun unazad (tzv. “backward sweep”). U slučaju burnog tečenja, postupak rešavanja za mirno tečenje, prikazan u prethodnom poglavlju, se mora modifikovati, pri čemu je dovoljan samo proračun unapred (“forward sweep”). Nadalje se razmatra samo slučaj mirnog tečenja. Postoje tri tipa graničnih uslova: 1) Poznati (zadati) protok kao funkcija vremena Q( t ) (što je čest uzvodni granični uslov - uzvodni hidrogram), 2) Poznata (zadata) kota slobodne površine kao funkcija vremena Z ( t ) (što je čest nizvodni granični uslov - nizvodni nivogram, ušće kanala u jezero i sl.). 3) Poznata (zadata) veza protoka i kote Q( Z ) (što je takođe čest nizvodni granični uslov - kriva protoka, tečenje preko preliva na nizvodnom kraju itd.). 130
Pošto su osnovne jednačine linearizovane, moraju se linearizovati i jednačine graničnih uslova. Uopšteno će se granični uslov u linearizovanom obliku formulisati kao:
α ∆ Z + β ∆Q = γ pa će se za svaki tip graničnog uslova izvesti koeficijenti α, β, γ. Na taj način postupak rešavanja, opisan u prethodnom poglavlju, ostaje nepromenjen, a za različite granične uslove menjaju se samo vrednosti koeficijenata α, β, γ. U nastavku će se razmotriti sva tri tipa graničnih uslova. 1) Koristeći poznati hidrogram Q( t ) , na granici treba zadati vrednost protoka u nepoznatoj ( m + 1) iteraciji u razmatranom vremenskom nivou t n+ 1 . Granični uslov se onda može pisati kao: m+1 n +1 Q = Q(t) t n +1
gde je Q ( t )
t n +1
vrednost protoka određena iz poznate funkcije Q( t ) za t = t n+ 1 .
Pošto je: m+1
Q n+1 = mQ n+1 + ∆Q
dobija se: m
Q n+1 + ∆Q = Q ( t )
tn +1
pa se granični uslov može pisati u opštem obliku:
α ∆Z + β ∆Q = γ ako je:
α =0
β =1
γ = Q(t)
tn +1
− m Q n +1
2) Koristeći poznati nivogram Z ( t ) , na granici treba zadati vrednost kote u nepoznatoj (m+1) iteraciji u razmatranom vremenskom nivou t n+ 1 . Granični uslov se može pisati kao: m+1
Z n+1 = Z ( t )
tn +1
gde je Z ( t ) tn +1 vrednost kote određena iz poznate funkcije Z ( t ) za t n+ 1 . Pošto je: m+ 1
Z n+1 = mZ n +1 + ∆Z
dobija se: m
Z n+1 + ∆Z = Z ( t ) 131
t n +1
pa se granični uslov može pisati u opštem obliku:
α ∆ Z + β ∆Q = γ ako su:
α =1
γ = Z (t)
β =0
t n +1
− m Z n +1
3) Koristeći poznatu zavisnost Q = f ( Z ) , na granici treba zadati vezu protoka i kote u nepoznatoj ( m + 1) iteraciji u razmatranom vremenu t n+ 1 . Granični uslov se može pisati kao: m+ 1
gde je f
( m+1Z n+1 ) vrednost funkcije
Q n +1 = f
( m+1Z n+1 )
f ( Z ) određena za Z = m+1Z n+1 .
Nelinearna funkcija f ( Z ) se može razviti u Taylor-ov red u okolini m Z n+1 :
f
( m+1Z n+1 ) = f ( m Z n+1 ) +
∂ f ∂Z
∆Z + .... m
Z
n +1
gde je: ∆Z = m+1Z n+1 − mZ n+1 Ako se zanemare članovi sa izvodima drugog i višeg reda, i pošto je:
m+1
Q n+1 = mQ n+1 + ∆Q
dobija se: m
Q n +1 + ∆Q = f
(
m
)
Z n +1 +
∂f ∂Z
⋅ ∆Z m
Z n +1
što se može pisati u opštem obliku:
α ∆Z + β ∆Q = γ ako je:
α =
∂ f ∂Z
β = −1 m
γ = mQ n+1 − f
( m Z n+1 )
Z n +1
Uočava se da će u sva tri slučaja koeficijent γ imati vrednost približnu nuli, ali ga ne treba nasilno izjednačavati sa nulom, jer bi se time akumulirala numerička greška. U nastavku će se ukazati na par praktičnih problema (i sugerisati rešenje) koji se javljaju pri korišćenju graničnog uslova tipa poznato Z(t).
132
Započinjanje proračuna unapred (određivanje koeficijenata E1, F1) Uzvodni granični uslov u opštem obliku:
α u z ∆Z1 + β u z ∆Q1 = γ u z napisan je kao: ∆Q1 = E1 ∆Z1 + F1 gde je:
E1 = −
αuz β uz
F1 =
γ uz β uz
Pošto je za poznato Z ( t ) koeficijent β = 0, javiće se problem deljenja sa nulom, pa se postupak rešavanja za ovaj slučaj mora modifikovati. Postupak rešavanja će se modifikovati tako što će se lokalno izmeniti uloge ∆Q1 i ∆Z1, tj. granični uslov će se napisati kao:
(e
∆Z1 = e1 ∆Q1 + f 1
1
≠ E1
f 1 ≠ F1 )
gde je:
e1 = −
β uz α uz
f1 =
γ uz αuz
(zapravo e1 = 0 jer je βuz = 0, tj. u ovom slučaju ∆Z1 ne zavisi od ∆Q1). Iz jednačina (1) i (1′): − C1 ∆Z1 − D1 ∆Q1 + A1 ∆Z 2 + B1 ∆Q2 = G1 − C1′ ∆Z1 − D1′ ∆Q1 + A1′ ∆Z 2 + B1′ ∆Q2 = G1′
(1) (1’)
će se eliminisati ∆Z1, tako da se dobije: ∆Q1 = l1 ∆Z 2 + m1 ∆Q2 + n1
(l
1
≠ L1
gde su:
l1 =
C ′A − CA′ C ′D − CD′
m1 =
C ′B − CB ′ C ′D − CD ′ 133
m1 ≠ M 1
n1 ≠ N 1 )
n1 =
CG ′ − C ′G C ′D − CD′
Takođe, ako se ∆Z1 = e1 ∆Q1 + f 1 uvrsti u jednačine (1) i (1′), pa onda iz te dve jednačine eliminiše ∆Q1, dobija se: ∆Q2 = E 2 ∆Z 2 + F2 gde su:
E2 =
F2 =
A′D − AD′ BD′ − B ′D
(Cf 1 + G )D ′ − (C ′f1 + G ′)D BD ′ − B ′D
Nadalje je postupak rešavanja isti kao i kod originalnog postupka, izloženog u prethodnom poglavlju. Kao i kod originalnog postupka, osim za slučaj kada tečenje teži kritičnom, uvek je BD ′ − B ′D ≠ 0 t.j. BD′ ≠ B ′D , pa nema deljenja sa nulom.
Započinjanje proračuna unazad (sračunavanje ∆ZI, ∆QI) Proračun unapred završava se određivanjem koeficijenata EI, FI u jednačini: ∆QI = E I ∆Z I + FI
(a)
Jednačina nizvodnog graničnog uslova u opštem obliku glasi:
α n i z ∆Z I + β n i z ∆QI = γ n i z (b) Eliminisanjem ∆QI iz jednačina (a) i (b) dobija se:
γ ni z − FI β ni z ∆Z I = α ni z EI + β ni z
(*)
Uočava se da je i ovde, za granični uslov tipa poznato Z ( t ) , koeficijent β = 0, pa bi se dobilo deljenje sa nulom. Međutim u ovom slučaju originalna jednačina graničnog uslova glasi: m
Z In+1 + ∆Z I = Z ( t )
134
t n +1
pa se umesto jednačine (*), u slučaju graničnog uslova tipa zadato Z(t) na nizvodnoj granici jednostavno može pisati: ∆Z I = Z ( t )
t n +1
− mZ In+1
Kada se iz prethodnog izraza sračuna ∆ZI, jednostavno se odredi: ∆QI = E I ∆Z I + FI i započne proračun unazad.
135
ZADACI
136
Zadatak 1 Rešiti problem provođenja toplote duž konzole, opisan u poglavlju 1.1., korišćenjem Galerkin-ove varijante metoda konačnih elemenata sa linearnim baznim funkcijama ϕ . Koraci: 1. Obaviti neophodna integraljenja tako da se dobije sistem linearnih algebarskih jednačina sa nepoznatim temperaturama u računskim čvorovima. 2. Rešiti sistem algebarskih jednačina, tako da se dobiju rezultujuće vrednosti temperature u računskim čvorovima. 3. Nacrtati približno rešenje dobijeno metodom konačnih elemenata, kao i tačno analitičko rešenje u istoj razmeri (duž iste koordinatne ose x ). 4. Uporediti približno rešenje dobijeno metodom konačnih elemenata sa tačnim (analitičkim) rešenjem i analizirati (objasniti) slaganje ili neslaganje. Posebno objasniti kako i gde je pri rešavanju upotrebljen granični uslov na uklještenom kraju konzole.
137
Zadatak 2 1. Napraviti algoritam za približno rešavanje problema provođenja toplote duž konzole iz Zadatka 1, koristeći metod konačnih razlika. Za drugi izvod koristiti aproksimaciju drugog reda, odnosno tzv. Crank-Nicholson-ovu aproksimaciju. Posebnu pažnju obratiti na primenu „prirodnog“ graničnog uslova na uklještenom kraju. 2. Sračunati približno rešenje prvo za četiri računske tačke (kao u Zadatku 1), a onda sa sedam računskih tačaka (tj. podeliti dužine računskih deonica na pola). Zadatak se može, ali i ne mora rešiti pomoću računara. 3. Nacrtati oba rešenja dobijena primenom metoda konačnih razlika, rešenje iz Zadatka 1 dobijeno primenom metode konačnih elemenata, kao i tačno (analitičko) rešenje, sve u istoj razmeri (duž iste koordinatne ose x ). 4. Uporediti približna rešenja (dobijena metodom konačnih razlika i metodom konačnih elemenata) sa tačnim (analitičkim) rešenjem. Analizirati slaganje ili neslaganje približnih rešenja sa tačnim rešenjem. Oceniti odnos tačnosti i uloženog truda za dva približna rešenja.
138
Zadatak 3 Razmatra se takozvana „upwind“ šema metoda konačnih razlika, primenjena na jednodimenzionalnu jednačinu advekcije: ∂C ∂C +u⋅ =0 ∂t ∂x Navedena šema koristi sledeće aproksimacije za izvode po vremenu i prostoru:
∂C C in +1 − C in ≅ ∂t ∆t n ∂C C i − C in−1 ≅ ∂x ∆x 1. Utvrditi da li je razmatrana šema konsistentna. 2. Sprovesti Von-Neumann-ovu analizu stabilnosti razmatrane šeme. 3. Sračunati vrednosti i napraviti grafički prikaz „amplitudne“ i „fazne“ greške za Courant-ove brojeve Cr = 0.25, 0.5, 0.75, 1.0 i 2.0 (i
λ
= 2, 4, 6,..., 30 ). ∆x 4. Na osnovu Von-Neumann-ove analize stabilnosti odnosno prikazane „amplitudne“ i „fazne“ greške izvesti zaključke o uslovima stabilnosti razmatrane šeme.
139
Zadatak 4 Primeniti „upwind“ šemu metoda konačnih razlika na linearnu jednodimenzionalnu jednačinu advekcije. Rešiti jednačinu advekcije:
∂C ∂C +u⋅ =0 ∂t ∂x za 0 < x < X i 0 < t < T , poznato C (0, t ) kao granični uslov i poznato C ( x, 0) kao početni uslov. „Upwind“ šema metoda konačnih razlika koristi sledeće aproksimacije za izvode po vremenu i prostoru: ∂C C in +1 − C in ≅ ∂t ∆t
∂C C in − C in−1 ≅ ∂x ∆x
1. Napisati program koji rešava jednačinu advekcije koristeći „upwind“ šemu metoda konačnih razlika, za konstantne ∆x i ∆t . 2. Koristiti program da se proračuna advekcija početnog rasporeda koncentracije C , zadatog Gauss-ovom krivom:
− ( x − x 0 )2 C = 10 ⋅ exp . 2 2 ⋅ σ x (gde se sa x0 obeležena koordinata maksimalne koncentracije) u kanalu dužine X = 10 km , za T = 160 min . U trenutku t = 0 maksimalna vrednost koncentracije zadate Gauss-ovom krivom se nalazi u x0 = 0 , tako da se tačno polovina (simetrične)
Gauss-ove krive nalazi unutar kanala, iz čega se može naći početni uslov C ( x, 0) . Za t > 0 preostala polovina Gauss-ove krive ulazi u kanal, iz čega se može naći granični uslov C (0, t ) . Usvojiti σ x = 264 m , u = 0.5 m / s , ∆x = 200 m . Sprovesti proračune za ∆t = 100 s , 200 s , 300 s , 400 s i 800 s (tj. za odgovarajuće Courant-ove brojeve). Nacrtati sve sračunate krive u istoj razmeri (duž iste ose x ), zajedno sa tačnim rešenjem, zbog poređenja. 3. Analizirati rezultate koristeći grafički prikaz „fazne“ i „amplitudne“ greške „upwind“ šeme iz Zadatka 3. Da li su se rezultati mogli predvideti na osnovu „fazne“ i „amplitudne“ greške „upwind“ šeme metode konačnih razlika. Zašto rezultati za ∆t = 100 s i ∆t = 300 s nemaju istu grešku u T = 160 min ? Ako bi se koristila „upwind“ šema za rešavanje praktičnog problema, kako bi se mogla smanjiti greška?
140
Zadatak 5 Testirati Holly-Preissmann-ov metod za čistu advekciju. 1. Izvesti radne jednačine (tj. one jednačine koje će se programirati) za HollyPreissmann-ov metod za čistu advekciju. 2. Napisati program i ponoviti proračune iz Zadatka 4 koristeći Holly-Preissmann-ov metod. 3. Porediti grafički rezultate dobijene primenom „upwind“ šeme metode konačnih razlika sa rezultatima dobijenim primenom Holly-Preissmann-ovog metoda; diskutovati upoređene rezultate. 4. Ponoviti jedan od proračuna primenom Holly-Preissmann-ovog metoda, ali ovog puta sa početnim (u t = 0 ) izvodima koncentracije po prostoru ( CX ) jednakim nuli. Diskutovati uticaj ove nedoslednosti na rezultat.
141
Zadatak 6 Primeniti metod razlomljenih koraka, tj. razdvojenih operatora na jednodimenzionalnu jednačinu transporta zagađivača (jednačinu advekcije-difuzije) tako što će se kombinovati Holly-Preissmann-ov metod za advekciju i Crank-Nicholson-ova šema za difuziju. 1. Razmatra se kombinovana advekcija i difuzija početne raspodele zagađivača koji je u početnom trenutku t 0 ispušten u jednodimenzionalnu struju fluida. Analitičko rešenje kombinovane advekcije i difuzije dati su sledećom funkcijom: C ( x, t ) =
− ( x − u ⋅ t − x 0 )2 ⋅ exp D x ⋅ π ⋅ (t 0 + t ) 4 ⋅ D x ⋅ (t 0 + t ) C0
gde je t 0 = 3484.8 s ; D x = 10 m 2 / s ; C 0 = 3308.75 jedinica (proveriti dimenziju);
2.
3.
4. 5.
u = 0.5 m / s ; x0 = 1400 m . Odrediti početni raspored zagađivača iz analitičkog rešenja za t = 0 i pretpostaviti da zagađivač ne ulazi u razmatranu oblast kroz uzvodnu granicu. Izvesti radne jednačine za metod razlomljenih koraka (razdvojenih operatora) zasnovan na korišćenju Holly-Preissmann-ovog metoda za advekciju i uopštene Crank-Nicholson-ove šeme za difuziju. Napisati program za metod razlomljenih koraka i koristiti ga za proračun advekcije i difuzije opisane raspodele zagađivača tokom perioda od 9600 s . Sprovesti proračune za θ = 0, 0.5, 1 . Usvojiti računsku mrežu sa konstantnom dužinom računskog koraka po prostoru od 200 m i koristiti konstantan računski korak po vremenu od 200 s . Porediti grafički približno rešenje sa analitičkim u t = 9600 s , za sve tri vrednosti θ . Ponovit proračune, ali ovoga puta za slabu difuziju, sa koeficijentom difuzije D x = 1 m 2 / s . Ponovo grafički porediti približno i tačno rešenje. Za ovaj koeficijent difuzije parametri analitičkog rešenja su: t 0 = 34848 s ; C 0 = 3308.75 jedinica
(proveriti dimenziju); u = 0.5 m / s ; x0 = 1400 m . 6. Diskutovati tačnost metoda, osetljivost metoda na vrednost θ , kao i osetljivost na jačinu difuzije. Da li u ovom slučaju potencijalna nestabilnost Crank-Nicholson-ove šeme za male vrednosti θ predstavlja problem?
142
Zadatak 7 Izvesti linearizovane jednačine tečenja za Preissmann-ovu šemu. 1. Proveriti (ponovo izvesti) koeficijente A, B, C , D, G u linearizovanoj jednačini kontinuiteta, kao što je to urađeno u poglavlju 3.5. 2. Izvesti koeficijente A' , B ' , C ' , D ' , G ' u linearizovanoj dinamičkoj jednačini. Za nagib trenja koristiti osrednjavanje predloženo u poglavlju 3.5. Voditi računa da se u članu:
α⋅
Q 2 ∂A ⋅ A 2 ∂x
(poglavlje 3.5.) koristi diskretizacija koja daje dobru energetsku jednačinu pri ustaljenom tečenju. 3. Sprovesti sledeće provere koeficijenata: a. Pri ustaljenom tečenju koeficijenti moraju da se svedu na odgovarajuću jednačinu kontinuiteta i energetsku jednačinu. b. Koeficijenti moraju biti dimenzionalno konsistentni (članovi u jednačinama moraju imati iste dimenzije).
143
Zadatak 8 Primeniti Preissmann-ovu šemu, Newton-Raphson-ov iterativni postupak i Thomas-ov algoritam (double-sweep) algoritam za računanje jednodimenzionalnog neustaljenog tečenja u neprizmatičnom kanalu pravougaonog poprečnog preseka. 1. Napraviti program za rešavanje jednačina kojima se opisuje jednodimenzionalno neustaljeno tečenje u otvorenim tokovima, koristeći Preissmann-ovu šemu, NewtonRaphson-ov iterativni postupak i Thomas-ov algoritam (double-sweep) algoritam. Program treba da omogući da se kao uzvodni granični uslov može zadati hidrogram Q(t ) , a kao nizvodni granični uslov bilo Q(t ) , Z (t ) ili Q(Z ) . 2. Koristiti napravljen program za računanje propagacije talasa u dugačkom hrapavom kanalu obloženim kamenom. Kanal je dugačak 24 km , Manning-ov koeficijent hrapavosti je 1 15 m −1 / 3 s , nagib dna kanala je konstantan i iznosi 0.0005 . Širina kanala se menja linearno od 8 m na uzvodnom kraju do 20 m na nizvodnom kraju. Kao uzvodni granični uslov koristiti sledeći hidrogram Q(t ) : t<0 t (min) 0 15 30 45 60 75 90 105 3 Q (m /s) 100 100 100 100 250 350 300 250 200
120 150
135 100
t > 135 100
Koristiti korak po prostoru 1 km , i parametre β = 0.5 , θ = 0.55 . Početni uslov (u t = 0 ) za sve proračune zadat je sledećim analitičkim izrazima: Q = 100 m 3 / s
Z ( x ) = 23.741 − 0.9653 ⋅ x + 0.0097 ⋅ x 2
gde je x u km, Z u m, x = 0 na uzvodnom kraju kanala. Usvojiti da je kota dna na nizvodnom kraju jednaka nuli. Napraviti proračune za period od t = 0 do t = 420 min , za sledeće uslove: a. Računski korak po vremenu je ∆t = 15 min , kao nizvodni granični uslov koristiti vezu protoka i kote koja važi za jednoliko tečenje: Q = K ⋅ S0 b. Računski korak po vremenu je ∆t = 15 min , kao nizvodni granični uslov koristiti nivogram Z (t ) napravljen korišćenjem rezultata proračuna pod a. c. Računski korak po vremenu je ∆t = 15 min , kao nizvodni granični uslov koristiti hidrogram Q(t ) napravljen korišćenjem rezultata proračuna pod a. d. Ponoviti proračun pod a., ali sa ∆t = 5 min . e. Ponoviti proračun pod a., ali sa ∆t = 30 min . Napomene: • Ako program radi kako treba, onda proračuni pod a., b. i c. treba da daju istovetne rezultate. • Tokom prvih 30 minuta proračuni treba da pokažu približno ustaljeno tečenje, jer zadati početni uslov (približno) predstavlja liniju nivoa za zadati kanal pri ustaljenom tečenju od Q = 100 m 3 / s . 3. Pripremiti dva crteža koji porede sračunate nivograme i hidrograme na nizvodnom kraju kanala za proračune pod a., d., i e. Koji računski korak po vremenu je najpogodniji za zadati problem? 144
Zadatak 9
Glavni tok 8 km
16 km
QHtL
Retenzija
QHZL
2 km 3 km
Prag Pritoka
10 km
QHtL Primeniti Preissmann-ovu šemu na granatu mrežu otvorenih tokova sa „unutrašnjim“ graničnim uslovima. 1. Izvesti „radne“ jednačine kojima se opisuje ušće kritoke u glavni tok, preliv odnosno prag sa pravougaonim poprečnim presekom i bočnu retenziju. 2. Modifikovati program iz Zadatka 8 tako da može da simulira konfiguraciju vodotoka i „unutrašnjih“ graničnih uslova kao na slici. 3. Polazeći od približnog početnog uslova, zadatog u nastavku, sprovesti proračune neustaljenog tečenja sa konstantnim graničnim uslovima, tako da se nakon tzv. perioda stabilizacije, dobije stvarno ustaljeno stanje kao stvarni početni uslov u trenutku t = 0 . 4. Sprovesti proračun neustaljenog tečenja od vremena 0 h do 10 h , koristeći θ = 1.0 , β = 0.5 , α = 1.0 i ∆t = 15 min . 5. Proveriti da li su svi „unutrašnji“ i spoljašnji granični uslovi zadovoljeni, tako što će se pokazati da sračunati protoci i nivoi slobodne površine vode zadovoljavaju odgovarajuće hidrauličke zavisnosti. 6. Pripremiti i prikazati sledeće grafike: a. Hidrograme na oba uzvodna kraja modela, na nizvodnom kraju modela, kao i u računskim tačkama neposredno uzvodno i neposredno nizvodno od bočne retenzije. b. Nivograme na oba uzvodna kraja modela i na nizvodnom kraju modela. c. Podužne profile slobodne površine vode u glavnom toku za vremenske trenutke 0 h , 1 h i 5 h . d. Podužne profile slobodne površine vode u pritoci za vremenske trenutke 0 h , 1h i 5 h. 7. Diskutovati pojavu i posledice neustaljenog tečenja u modelisanom sistemu (punjenje i pražnjenje bočne retenzije, režim prelivanja preko praga, ublaženje talasa itd.). Konfiguracija vodotoka i graničnih uslova.
A. Glavni tok je isti kao u Zadatku 8. 145
B. Pritoka: izlaz iz energetskog postrojenja, kanal pravougaonog poprečnog preseka, 10 km dug i 20 m širok, sa betonskom oblogom (Manning-ov koeficijent hrapavosti n = 0.02 m −1 / 3 s ), koji se uliva u glavni tok na 16 km od uzvodnog kraja glavnog toka (kao na slici). Pritoka ima računske tačke zadate u sledećoj tabeli.
Stacionaža (km) 0.0 (uzvodni kraj) 2.0 3.5 5.5 8.0 10.0 (ušće)
Koda dna (m) 9.0 7.0 6.5 6.0 5.0 4.0
C. Bočna retenzija je locirana 3 km od uzvodnog kraja glavnog toka. Ima konstantnu površinu od 5 ha , i kotu dna 10.75 m . D. Prag se nalazi u pritoci, na 2 km uzvodno od ušća u glavni tok. Kota krune praga je Z p = 16 m , širina B = 20 m , a koeficijent prelivanja je µ = 1.0 . E. Uzvodni hidrogram (granični uslov) na pritoci je: Q (m3/s) 150.0 550.0 150.0
Vreme (h) −∞
Napomena: Linearno povećanje odnosno smanjenje protoka tokom uključenja i sključenja, koji pojedinačno traju trideset minuta. F. Granični uslovi na glavnom toku su isti kao u Zadatku 8 pod tačkom a. G. Približni početni uslov: o protoci zadati na osnovu graničnih uslova za trenutak t = 0 , o glavni tok: konstantna dubina vode od 12 m , o pritoka: konstantni nivo slobodne površine vode od 18.8 m uzvodno od praga, odnosno 17.0 m nizvodno od praga.
146
