Hidrológiai Közlöny 1993 (73. évfolyam)
6. szám - Szél Sándor–Gáspár Csaba: A Muskingum-Cunge eljárás felülvizsgálata és továbbfejlesztése
SZÉL S. - GÁSPÁR CS.: A Muskingum-Cunge eljárás 349 - 1 + Cr + C-K = Pe 1 + Cr + Pe CAr 2 C 3A/ 2 (15) 2 ' 2 1 " 6 6 egyenlőség teljesüljön. Ezt Ax-szel osztva nyeijük a [DCAx (Cr -1)] • [(l-e) • Cr + (1-©)] 1 D • Cr - • (Cr 2 - 1) egyenlőséget. Innen, a (13) feltételt figyelembe véve könnyen igazolható, hogy ez az egyenlőség akkor és csakis akkor áll fenn, ha a következő egyszerű feltétel teljesül: 3 Cr 2 = 1 — Ha most Ax, At egyikének értékét valamilyen megfontolás alapján előírjuk, akkor a (16) feltétel egyértelműen meghatározza a másiknak az értékét, éspedig (16)-ból nyilvánvalóan: Megjegyzések: 1. A fenti (ll)-es (Taylor-sorfejtésen alapuló) levezetést szokásos (Cunge, 1969) a (2) egyenlet helyett a jóval egyszerűbb alakú tiszta konvekciós egyenlet megoldásaira levezetni. Ekkor, nagy vonalakban, a meggondolás az, hogy a séma alkalmatlan a tiszta konvekciós problémák kezelésére, mert jelentékeny numerikus diszperzió terheli: ezzel szemben viszont alkalmas olyan transzport esetén, ahol a fizikai diszperzió egyezik az így kiadódó numerikus diszperzióval. Ez a meggondolás elvileg hibás, ui. ebből csak az következik, hogy a differenciaoperátor „jól" (legalábbis másodrendben) approximálja az L konvekciós-diszperziós operátort egy másik differenciálegyenlet (ti. a tiszta konvekciós egyenlet) megoldásain. A (13) feltétel ennek ellenére ugyanilyen formában kiadódik, a magasabb rendű közelítés feltétele (ld. alább) viszont már nem. 2. Miután a (13) feltétel e és 0 értékét külön-külön nem szabja meg, csupán egy lineáris összefüggést ír elő köztük, az e, 0 súlyozóparaméterek szokatlan módon is megválaszthatok, azaz akár negatívak vagy 1-nél nagyobbak is lehetnek, így pl. sokszor választják e értékét 0,5-nek: ekkor a (13) feltétel 0-ra mindig negatív értéket ad, valahányszor Pecl, ami ütközni látszik a fizikai szemlélettel, illetve ellentmondást sejtet (Szilágyi, 1991). Valójában itt szó sincs ellentmondásról, sem pedig hibás alkalmazásról. Ha a séma kellő mértékben pontos és stabil, akkor az konvergens (a matematikailag egzakt tárgyalásra nézve ld. pl. Richtmyer és Morton, 1967), így a gyakorlatban jól használható még akkor is, ha negatív súlyozóparamétereket tartalmaz (ld. Szél, 1988). A numerikus stabilitás kérdésével egy későbbi szakaszban foglalkozunk. A numerikus oszcilláció kiküszöbölése Ha a (13) feltétel teljesül, akkor a fentiek szerint a sémának nincs numerikus diszperziója, és legalább másodrendben közelít. A séma pontossága tovább növelhető, ha sikerül elérni, hogy (11) jobb oldalán a másodrendű hibatag is eltűnjék. Világos ugyanakkor, hogy mivel a numerikus diszperzió eliminálása a-t már meghatározza, a numerikus oszcilláció eltűnését a súlyozóparaméterek megválasztásával már nem lehet elérni, ahhoz a séma további paraméterei (pl. Ax vagy At) speciális megválasztása szükséges. A numerikus oszcilláció eltűnésének feltétele (11) alapján az, hogy a éi.í.^-m* ill. I2zr 1/2 Ar + P? (16) Megjegyzések: 1. A (16) feltétel tiszta konvekciós probléma (pl. kinematikus hullám) esetén a Cr = 1 feltételbe megy át: jól ismert tény, hogy ekkor a séma valóban alakhűen visz át éles frontokat is, numerikus oszcilláció generálása nélkül. 2. A (16) feltétel nem mindig teljesíthető: ehhez nyilván az szükséges, hogy Pe 2 > 3 legyen (az egyenlőség zérus időlépést eredményezne). Ez megszorítást jelent az alkalmazható térlépésre. A (16) feltétel egyébként meglehetősen erős abban az értelemben, hogy ekkor a Courant-szám sohasem nagyobb 1-nél, ami azzal jár, hogy az alkalmazandó időlépés sokszor igen kicsi kell, hogy legyen: ez a szükséges számítási munkát erősen megnövelheti. Az eredményeket összefoglalva: az (5) séma mentes a numerikus diszperziótól és másodrendben közelíti a (3) differenciáloperátort a (2) egyenletnek tetszőleges megoldásán, ha a (13) feltétel teljesül. Ekkor a séma együtthatóit a (14) formulák adják. Ha még ezenfelül a (16) feltétel is teljesül, akkor a séma mentes a numerikus oszcillációtól, és harmadrendben közelíti az L operátort a (2) egyenlet tetszőleges megoldásán. A séma használhatóságához a jó approximáció azonban még kevés: az is szükséges, hogy a séma numerikusan stabil legyen. Ennek vizsgálatát végezzük el a következő szakaszban. A Muskingum-Cunge-séma stabilitása A stabilitás, durván szólva, azt jelenti, hogy bármilyen, térben és időben korlátos (fizikailag nem feltétlen reális) kezdeti, illetve peremfeltétel esetén a séma olyan diszkrét megoldást eredményez, amely szintén korlátos marad, éspedig a rácsháló finomságától függetlenül. Ellenkező esetben az óhatatlanul következő számítási pontatlanságok hatása minden határon túl nő, és tönkreteszi a közelítő megoldást. Ismeretes, hogy a stabilis sémák viszont konvergens módszereket eredményeznek (Marcsuk, 1976; Peyret és Taylor, 1983) mégpedig annál pontosabbakat, minél jobb a séma approximációja a megoldani kívánt differenciálegyenletek megoldásain. A továbbiakban tegyük fel, hogy - mint a numerikus modellezések használatakor mindig - a vizsgált térszakasz, melyen a (2) egyenlet megoldását keressük, nem a félig végtelen (0, + °o) félegyenes, hanem a véges (0, A) szakasz. Tegyük fel továbbá, hogy e szakasz épp N számú részre van osztva az x k rácspontokkal, azaz
