Vývoj numerického 2D-ADT transportního modelu plavenin Jiří Lindner České vysoké učení technické v Praze, Fakulta stavební, Katedra hydrotechniky e-mail:
[email protected]
Abstrakt V tomto příspěvku je popsán úspěšný vývoj nového simulačního prostředku pro modelování transportu drobných plaveninových částic (popř. znečistění) ve vodních tocích. Modelový prostředek byl pojmenován 2D-ADT (dvojdimenzionální - advekčně difúzní transportní) model a jedná se o transportní model s integrovanou sedimentační a erozní složkou. Tento numerický prostředek vznikl na základě aktuální vodohospodářsko ekologické problematiky. Vývoj modelu byl zaměřen na mapování šíření plaveninových částic ve vodních tocích, na kterých může být vázáno chemické znečištění (např. AOX, ΣPCB, NEL, ΣDDT). 2D-ADT model je založen na advekčně difúzní rovnici (1), která byla řešena metodou konečných diferencí, explicitním QUICKEST schématem. Programově byl 2D-ADT model vytvořen v programovacím jazyce C# a pro vynesení výstupů z tohoto modelu bylo naprogramováno příslušné rozhraní k postprocesoru, užívaném na pracovišti Katedry hydrotechniky, ČVUT v Praze.
1
Úvod
Přiblížit se k vyřešení problému transportu plavenin (popř. znečistění) ve vodních tocích můžeme jednak za pomocí fyzikálního modelování nebo s použitím numerického modelu. Zaměříme-li se na druhou možnost, tj. numerické modelování, zjistíme, že po správném určení řídící diferenční rovnice tohoto advekčně difúzního jevu nám vyvstane otázka správného aplikování optimální numerické metody a jejího diferenčního schématu. Cílem této práce proto bylo nalezení optimálního řešení 2D advekčně difúzní rovnice, dle kritérií přesnosti, stability a náročnosti řešení, které by bylo možno posléze aplikovat v novém simulačním numerickém prostředku pro modelování transportu drobných plaveninových částic ve vodních tocích. Pro řešení 2D úloh šíření látky v toku byla vybrána jako nejvhodnější metoda konečných diferencí se dvěma diferenčními schématy. Jako první byl problém řešen explicitním Upwind schématem, následovalo řešení za pomocí rovněž explicitního QUICKEST schématu. 2D advekčně difúzní (2D A-D) rovnice v otevřených korytech pro transport lze psát v následujícím tvaru [2]: ∂c ∂ 2c ∂ 2c ∂c ∂c = Dx 2 + D y 2 − v x − vy ∂t ∂x ∂y ∂x ∂y c (x,y,t) koncentrace Dx, Dy, difúzní koeficient vx, vy, (u, v) rychlosti ve směru x, y
2
(1)
Provedené simulace
Pro posouzení a porovnání těchto explicitních schémat byl sestaven jednoduchý numerický model, na kterém byly posléze provedeny základní simulace transportu látek při rovnoměrném ustáleném proudění. Za základ modelu byla navržena čtvercová síť, pole s 625
uzly, 25 x-směru a 25 y-směru, s jednorázovým nadávkováním znečištění c(x,y,t) = 1 kg m-2 v centru sítě. Při simulaci byly použity koeficienty podélné a příčné difúze Dx=Dy= 0.2 m2s-1. Vzdálenost mezi jednotlivými uzly ∆x = 2 m a ∆y = 2 m. Časový krok byl zvolen velice krátký (z důvodu použitých explicitních schémat) ∆t = 2,5 s. Celková doba simulace byla při 10 časových krocích rovna 25 s. Pro každé diferenční schéma byly provedeny dvě sady simulací. Jednak bylo modelováno šíření při nulových rychlostech u a v pro nasimulování transportu látek pouze difúzní složkou šíření. Posléze byly provedeny simulace transportu za nenulových rychlostí, v tomto případě za rychlostí kladných ve směru x a y. Zde potom transport probíhal za přispění jak advekční, tak difúzní složky transportu. Na závěr byl proveden analytický výpočet pro ověření obou explicitních schémat. Výstupy z výše popsaných simulací byly graficky zpracovány (Obr.1, 2, 4). Na těchto grafických výstupech dominuje znázornění koncentrace pomocí grafu proloženého v podélných a příčných osách sítě (osy x a y). V levém horním rohu je pozičně znázorněno šíření koncentrace po uplynutí celkového času simulace. Model po programové stránce spolu s dvojrozměrným grafickým výstupem byl vytvořen v programovacím jazyce C# (.net).
3
Explicitní Upwind (UP) schéma
Upwind schéma užívá pro výpočet standardu “upstream” (protiproudu) a 2D advekčně difúzní rovnici (1) je možné podle tohoto schématu rozepsat následovně [2]:
cin,+j1 = cin, j − Cxa (cin, j − cín−1, j ) − Cya (cin, j − cin, j−1) + Dxd (cin−1, j − 2cin, j + cin+1, j ) + Dyd (cin, j−1 − 2cin, j + cin, j+1 ) Advekční Courantovo číslo-
Crxa = Cxa =
Difúzní Courantovo číslo -
D xd =
vx ∆t ∆x
,
D xx ∆t , ∆x 2
Crya = Cya =
D yd =
vy ∆t ∆y
D yy ∆t ∆y 2
(2)
,
(3)
,
(4)
Obr.1: Simulace advekčně difúzního šíření látky explicitní Upwind metodou po 25 sekundách a při u = v = 0 ms-1 a u = v = 0,1 ms-1. Upwind schéma je stabilní, když jsou v časovém kroku splněny následující podmínky [2]: ∆t ≤
(
∆x 2
)
4 D + v x + v y ∆x
(5)
4
Explicitní QUICKEST schéma
Algoritmus QUICKEST (Quadratic Upstream Interpolation for Convective Kinematics with Estimated Streaming Terms) schématu je založen na formulaci “conservative control volume”, to znamená, že výpočet změny hodnoty jednoho uzlu v čase dosáhneme kvadratickou (tj. parabolickou) interpolací hodnot dvou sousedních uzlů, spolu s hodnotou nejblíže následujícího uzlu, nacházejícího se “proti proudu”. Rozepsání 2D A-D rovnice (1) pro QUICKEST schéma bude vypadat následovně [1]: 2 1 1 1 1 cin, +j 1 = cin, j − Crxa (cin, j + cín+1, j ) − Crxa (cín+1, j − cin, j ) − − D xd − (Crxa ) ⋅ (cín+1, j − 2cin, j + cín−1, j ) 2 6 6 2 2 1 1 1 1 + Cr xa (c in, j + c ín−1, j ) − Cr xa (c ín, j − c ín−1, j ) − − D xd − (Cr xa ) ⋅ (c ín− 2 , j − 2 c in−1, j + c ín, j ) 2 2 6 6
1 − Cr ya (c in, j + c ín, j +1 ) − 2 1 + Cr ya (c in, j + c ín, j −1 ) − 2
2 1 1 1 Cr ya (c ín, j +1 − c ín, j ) − − D yd − (Cr ya ) ⋅ (c ín, j +1 − 2 c in, j + c ín, j −1 ) 2 6 6 2 1 1 1 Cr ya (c ín, j − c ín, j −1 ) − − D yd − (Cr ya ) ⋅ (c ín, j − 2 c in, j −1 + c ín, j − 2 ) 2 6 6
+ D xd (c ín+1, j − 2c in, j + c ín−1, j ) + D yd (c ín, j +1 − 2c in, j + c ín, j −1 ) ,
(6)
Obr.2:Simulace advekčně difúzního šíření látky explicitní metodou QUICKEST po 25 sekundách a při u = v = 0 m s-1 a u = v = 0,1 m s-1. Numerický výpočet explicitního QUICKEST schématu je stabilní, pokud se nachází v tzv. “oblasti stability”, to znamená, pokud se Courantovo číslo Cr a difúzní parametr γ nachází v oblasti viz Obr.3. [1] u ⋅ ∆t v ⋅ ∆t = ∆x ∆y D ⋅ ∆t D ⋅ ∆t γ = x 2 = y 2 ∆x ∆y Cr =
Obr.3: Oblasti stability pro QUCKEST schéma [1] .
(7) (8)
5
Porovnání s analytickým řešením
Pro verifikaci obou numerických, explicitních schémat je zapotřebí ověřit získané výsledky i za pomocí analytického řešení provedeného na stejném příkladu. Poskytne nám to jednak možnost ověření správnosti numerického řešení, ale i možnost verifikování celého numerického modelu. Jestliže bude okamžité dávkování M = 1 kg provedeno ve středu sítě a při t = 0, pak šíření koncentrace v otevřeném recipientu bude dle Rutheforforda [3] dáno vztahem: ( x − v x ⋅ t )2 ( y − v y ⋅ t )2 M c (x, y, t ) = ⋅ exp− + 4 ⋅ D y ⋅ t 4 ⋅ π ⋅ t ⋅ Dx ⋅ D y 4 ⋅ D x ⋅ t (9)
Obr.4: Analytické řešení advekčně difúzního šíření látky po 25 sekundách a při u = v = 0 m s-1 a u = v = 0,1 m s-1.
6
Zhodnocení
Z dosažených výsledků a grafického znázornění advekčně difúzního transportu na Grafu 1. je dostatečně průkazné, že pro řešení tohoto a podobných případů je vhodnější explicitní QUICKEST schéma, které se více přibližuje analytickému řešení. Vhodné je jak svou dostatečnou přesností výpočtu (neprojevuje se u něj v takové míře numerická difúze jako u Upwind schématu), tak především svou nízkou náročností (v porovnání s implicitními schématy). Z těchto důvodů má tato metoda nejlepší předpoklady na využití v nově vznikajícím simulačním prostředku pro modelování transportu drobných plaveninových částic ve vodních tocích nazvaného 2D-ADT model.
Graf 1: Znázornění postupu koncentrace v ose x pro schémata Upwind, QUICKEST a pro analytické řešení.
7
Modelování transportu se sedimentačně erozním členem
Při modelování reálného šíření látek (např. v tocích se vyskytujících plaveninových částic) je pro výpočet transportu zapotřebí použít advekčně difúzní rovnici (1) se sedimentačně erozním členem Si(D,E). Hodnota sedimentačně erozního členu Si(D,E) bude rovna rozdílu do rovnice přidané koncentrace (vlivem eroze), a z rovnice odebrané koncentrace (vlivem sedimentace):
Si(D,E) = SE - SD
(10)
SD …..sedimentační člen vyjadřuje úbytek koncentrace způsobený sedimentací. Pro výpočet sedimentace byla dle Krone - Methy [5], použita rovnice: τ S D = w S ⋅ c ⋅ 1 − b g m-2s-1, (11) τ cd kde ws je sedimentační rychlost. K sedimentaci bude docházet v případě, kdy bude splněn níže popsaný poměr mezi tečným napětím a kritickým tečným napětím pro sedimentaci: τ b < τ cd . (12) SE…..erozní člen vyjadřuje nárůst koncentrace způsobený účinky eroze. Pro výpočet eroze byla použita dle Lumborga [6] rovnice: (α ⋅ τ b −τ ce ) SE = E ⋅ e g m-2s-1, (13) kde E je erozní konstanta. K erozi bude docházet v případě, kdy bude platit: τ ce < τ b , (14) -2 τ b je tečné napětí u dna N m a dle autorů Methy a Partheniadese [7] bude počítáno následovně: g τ b = ρ ⋅ 2 w2 , (15) C kde C je Chezyho rychlostní součinitel vztažený k hloubce a součiniteli drsnosti dle Manninga, a w je velikost vektoru průměrné svislicové rychlosti proudění. τ cd je kritické tečné napětí pro sedimentaci N m-2 a pro jeho výpočet může být použita rovnice dle Krone Methy [5] : (16) τ = θ ⋅ (ρ − ρ ) ⋅ g ⋅ d cd
s
w
,
τ ce je kritické tečné napětí pro erozi N m-2, ve většině případů jsou ovšem τ ce , τ cd , určeny pomocí laboratorních experimentů. Z výše uvedených vzorců (1), (10), (11), (12), (13), (14), (15), byl sestaven algoritmus pro řešení 2D advekčně difúzního transportu látek v tocích, viz. Obr. 5, který byl také použit v nově vyvinutém 2D-ADT transportním modelu.
Obr. 5: Řešící algoritmus pro 2D-ADT transportní model.
8
Vstupy a výstupy z 2D-ADT transportního modelu
Pro sestavení 2D-ADT transportního modelu pro šíření plaveninových látek v otevřeném korytě jsou zapotřebí tyto vstupní údaje: - kompletní návrh výpočetní sítě, se svými rozměry, s počtem výpočetních uzlů, se vzdálenostmi mezi jednotlivými uzly dx, dy, a jasně vymezenými okrajovými podmínkami - střední profilové rychlosti ve směru x, y, a hloubky v jednotlivých uzlech sítě - mapové podklady modelované oblasti ve formě DTM (digitální model terénu) - součinitel podélné a příčné difúze Dx, Dy v základních jednotkách m2s-1 - hodnoty počátečních koncentrací c pro každý bod sítě, v základních jednotkách kg m-3 - vstupní hodnoty pro sedimentační a erozní složku modelu (τ ce, ,τ cd , sedimentační rychlosti atd.) Výstupem z 2D-ADT transportního modelu po ukončení doby simulace jsou: - hodnoty výsledných koncentrací (po simulovaném čase) pro každý uzel v síti, které zobrazují šíření látky v toku - množství sedimentovaných popřípadě erodovaných částic v jednotlivých uzlech
9
Určení součinitelů difúze Dx a Dy
Přesné určení součinitelů difúze lze provést pouze na základě laboratorních pokusů. Proto se ve většině případů přistupuje k přibližnému určení těchto součinitelů, nebo k jejich odhadu. Pro model 2D-ADT byly součinitelé difúze získány z několika zdrojů a následně navzájem porovnány. Modelem 2D-ADT bylo ve většině případů simulováno šíření za izotropních podmínek, proto oba součinitelé, jak podélné tak příčné difúze, Dx, Dy byly navzájem rovny. Pro určení součinitelů difúze byly použity následující metody: 1. odhad součinitelů difúze dle Fischera [8], který sestavil tabulku Dx, Dy pro pokusné kanály. 2. určení součinitelů difúze pomocí empirických vztahů, dle Fischera [8]: u∗ = g ⋅ d ⋅ s , (17) Dx = Dy = 0,15 ⋅ d ⋅ u ∗ . (18) 3. určení difúzních součinitelů v podélném a v příčném směru z empirických rovnic popsaných Marvanem [3]: Dx = d L ⋅ u ∗ ⋅ h , DY = d T ⋅ u ∗ ⋅ h , kde je třecí rychlost u∗ vyjádřena pomocí následujícího vzorce:
u∗ =
τ0 . ρ
(19)
(20)
10 Praktické použití 2D-ADT modelu a zhodnocení Pro úspěšné otestování transportního 2D-ADT modelu byl řešen praktický úkol, šíření plaveninových částic, které pocházejí z resuspendovaných kohezivních sedimentů, na toku Třebovka nad soutokem s Tichou Orlicí. Ověření sestaveného modelu bylo provedeno prostřednictvím simulací pohybu těchto plavenin na fyzikálním modelu. Mezi výstupy z obou simulací byla konstatována podobnost a pro numerický model byly po vizuálním srovnání obou modelů vybrány součinitelé difúze Dx, Dy = 0,01 m2 s-1 jako optimálně se hodící pro dané fyzikální podmínky.
Druhá praktická úloha byla řešena na Labi v oblasti pod přehradou Les Království, kde bylo pomocí 2D-ADT modelu simulováno šíření sedimentů, při extrémní povodni, z usazených starých ekologických zátěží. Při těchto simulacích bylo prokázáno, že při modelovaném průtoku Q500 se znečištění dostane prakticky do celého říčního profilu, včetně inundací, jako je např. Jaroměřský rybník, kde se posléze plaveniny nejvíce ukládají. Rovněž bylo ověřeno, že část znečištění, vázané na plaveninové částečky, postupuje dále po toku. Obě tyto praktické aplikace ukázaly, že 2D-ADT model je schopen řešit šíření jemných plaveninových částic v tocích a v inundačních územích s reálnou topografií, které se s vysokou pravděpodobností podobá reálnému transportu plaveninových látek v českých tocích. Simulační prostředek 2D-ADT byl pomocí těchto praktických úloh úspěšně odzkoušen a je připraven na řešení široké škály inženýrských úloh, ať již v oblasti ekologie, transportu plavenin nebo ve stokování při modelování usazovacích nádrží. Pro komerční využití 2DADT modelu by bylo zapotřebí do programového vybavení tohoto simulačního prostředku doplnit uživatelské rozhraní.
Reference [1] Leonard, B. P., 1979. A Stable Accurate Convective Modelling Procedure Based On Quadratic Upstream Interpolation. Computer Methods in Applied Mechanics and Engineering, Vol. 19, pp 59-98. [2] Lindner, J., 2006. Two Dimension Equations (parabolic) for diffusion transport. HeriotWatt University of Edinburgh, UK [3] Marvan, F. G., 2001. A two-dimensional numerical transport model for organic-rich cohesive sediments in estuarine waters. Heriot-Watt University of Edinburgh, UK. [4] Rutherford, J. C., 1994. River Mixing. John Wiley & Sons, Chichester, England [5] Krone, R. B, Mehta, A. J.: Flume studies on the transport of sediment in estuarine shoaling processes, Hydraulik Engineering Laboratory, University of Berkeley, California, USA, 1962 [6] Lumborg, U.: Cohesive sediment transport modelling – application to the Lister Dyb tidal area in the Danish Wadden Sea, Journal of Coastal Research, Special Issue 41, pp 114123., 2004. [7] Mehta, A. J., Patheniades, E., An investigation of the depositional properties of flocculated fine sediment. Journal od hydraulic research, international Association of Hydraulic Research. Vol 13, No. 14, pp 361-381, 1975. [8] Fischer, B. H., Imberger, J., List, E. J., Koh, R. C. Y., Brooks, h. N.: Mixing in Inland and Coastal Waters. Academic Press, A Subsidiary of Harcourt Brace Jovanovich, Publishers, 1979
