19 JUNI
1985
BEREKENINGEN
VAN
BANEN
EN DISTRIBUTIES
VAN
DEELTJES IN
TWEE-~~-DRIEDIMENSIONALE
.
:nnm!!!;T
SNELHEI
o
QUIUMDAG
.
.
'; I \ J rFA 8506
BEREKENINGEN VAN BANEN EN DISTRIBUTIES VAN DEELT,JES IN TWEE- EN DRIE DIMENSIO-NALE SNELHEIDSVELDEN J.J. Bosman
*
G.C. van Dam**
H. Gerritsen***
W.A.M. de Jong**
R.A. Pasmanter**
L. Postma***
D. Riepma****
*
Waterloopkundig Laboratorium, de voorst**
Fysische Afdeling, Rijkswaterstaat***
Waterloopkundig Laboratorium, Delft3
-*
***
*
*
***
***
Deze bundel is een neerslag van een colloquiumdag gehouden op 19 ju-ni 1985, georgaju-niseerd door de Fysische Afdeling van de Directie Wa-terhuishouding en Waterbeweging van de Rijkswaterstaat, in samenwer-king met de OVerleggroep Transportverschijnselen van de Raad van Overleg van het fysisch-oceanografisch onderzoek van de NOordzee.
Enkele auteurs hebben volstaan met het inzenden van een in ander ka-der voorbereide tekst, naar inhoud overeenstemmend met de gehouden voordracht. Hierdoor komen enkele afwijkende titels voor; op een voorafgaand blad is steeds de oorspronkelij ke titel van de voor-dracht vermeld: ook in de inhoudsopgave (blz. 5) worden eventuele
afwijkende titels voorafgegaan door de oorspronkelijke.
De teksten zijn afgedrukt in de volgorde van de voordrachten op de colloquiumdag. Aan het eind volgt een kort hoofdstuk onder de titel "Elementen uit de diskussies". Daarbij is alleen verwerkt wat de ge-schreven teksten verduidelijkt of aanvult en wat zonder tijdrovende navraag en kontrole met redelijke zekerheid uit aantekeningen kon worden gerekonstrueerd, waarbij meer gestreefd is naar leesbaarheid en duidelijkheid dan naar een zo letterlijk mogelijke weergave. De inhoud van de toegeleverde teksten is uiteraard voor rekening van de sprekers/auteurs; zij dragen echter geen verantwoordelijkheid voor de weergave van de diskussie.
***
***
*
*
***
*
Inhoud
Inleiding
blz. G.C. v. Dam 7
Afleiding van snelheidsvelden uit op roosters gege-ven componenten
W.A.M. de Jong 13
Aspekten van tweedimensionale transportberekeningen voor oppervlaktewater door middel van (passieve) deeltjessimulatie
G.C. van Dam 29
Deterministic chaos: a new mechanism for dispersion in shallow seas
(Dynamical systems, Deterministic-Chaos and Dispersion in Shallow Tidal Seas)
R.A. Pasmanter 61
Seizoengemiddelde verplaatsingen in de centrale Noordzee
(Mean residual displacements during the winter season in the Central North Sea and possible year to year fluc-tuations. COmparison of two models (first results»)
H.W. Riepma 97
Eerste resultaten met een deeltjesbanenmodel t.b.v. het onderzoek naar subgrid-scale-effekten
(Some first particle track computations with fine grid and coarse grid velocity fields of the Southern North Sea)
H. Gerritsen 113
Enige bijzondere aspekten van het stoftransport voor de Zuid-Nederlandse Noordzeekust
L. Postma 137
Berekeningen van deeltjesdistributies voor H.W. Riepma 161 de Noordzee in drie dimensies
(Current meter records and the problem of the simulation of particle motions in the North Sea near the Dutch coast)
Transporttheorie via statistische aanpak van deterministische bewegingsvergelijking
J.J. Bosman 173
-
7-Algemene inleiding
Het lijkt passend, als inleiding enkele opmerkingen te maken bij de formulering van het thema.
In de formulering van het thema van deze dag is kortheidshalve weg-gelaten dat al het te rapporteren werk iets te maken heeft met op-pervlaktewateren.
Bepaalde uitspraken zullen ook wel een algemenere (vloeistofmechani-sche) strekking hebben. Sommige ideeën, technieken of konklusies zullen wellicht ook ruimer toepasbaar zijn, b.v. op verschijnselen in de atmosfeer.
Het kader oppervlaktewateren brengt mee dat in het vlak van de toe-passingen het veelal gaat om "vreemde" deeltjes, meer dan over de waterdeeitjes zelf.
Onderstaand enkele opmerkingen bij de "trefwoordenn van het thema BEREKENINGENI VAN BANENIV EN DISTRIBUTIESV VAN DEELTJESII IN
TWEE- EN DRIEDIMENSIONALE SNELHEIDSVELDENIII in relatie tot het programma en in de volgorde van de Romeinse cijfers.
Daarbij is de toedeling van de opmerkingen aan de trefwoorden hier en daar enigszins willekeurig.
I. BEREKENINGEN
Het zal inderdaad gaan over berekeningen en de daarbij gebruikte hulpmiddelen, i.h.b. de (wiskundige) modellen.
Maar daarbij speelt de empirie voortdurend een rol. Het gaat er ten-slotte om, de fysische werkelijkheid te begrijpen en er voorspellin-gen over te kunnen doen. De (transport-)modellen worden daarom geka-libreerd met behulp van meetresultaten: de rekenresultaten worden getoetst aan de werkelijkheid: snelheidsvelden worden afgeleid uit waterbewegingsmodellen die op hun beurt met waarnemingen gekali-breerd zijn, of de snelheidsvelden worden rechtstreeks uit waarne-mingen afgeleid.
8
-II. DEELTJES
Dit woord slaat vooral op de technieken waartoe we ons vandaag zul-len beperken. We kunnen er ons verschilzul-lende dingen bij voorstelzul-len, eventueel betrekkelijk grote waterpakketten, maar soms, in het
bij-zonder in de bijdrage van BOSMAN, gaat het om deeltjes die zich bij waarneming ook als zodanig aan ons voordoen en waarbij we ook in hun gedrag als zodanig geïnteresseerd zijn.
v66r- en nadelen van deeltjes-technieken t.o.v. andere methoden
ko-men in het algemeen niet aan de orde, misschien in de diskussies.
lIl. SNELHEIDSVELDEN
De snelheidsvelden worden afgeleid uit uitkomsten van
waterbewe-gLnqsmodelLerrê (DE JONG, VAN DAM, RIEPMA-l, GERRITSEN, POSTMA) of
uit waarnemingenb (RIEPMA-2, BOSMAN) , of ze worden in analytische
vorm gepostuleerdc (PASMANTER, VAN DAM, BOSMAN) •
In de eerste twee gevallen is sprake van "afgeleid uit" omdat het
bij deeltjes-technieken essentieel is dat de snelheid in elk gewenst
punt gegeven is. Daarom moet er in beide gevallen geïnterpoleerd
worden. Dit is geen triviale zaak, het geïnterpoleerde veld moet aan
bepaalde eisen voldoen, vooral bij stationaire situaties luistert
dit erg nauw. De voordracht van DE JONG is geheel aan dit onderwerp
gewijd.
De analytische velden zijn naar hun aard reeds als kontinue funkties
gegeven. Aan bepaalde eisen is voldaan, i.h.a. is gezorgd voor over-eenstemming met de kontinuïteitsvoorwaarde voor de vloeistof maar is
niet expliciet rekening gehouden met de dynamika. Maar struktuur,
frekwenties, golflengten etc. zijn wel gekozen met een blik op de
empirie.
Bij de snelheden berekend met waterbewegingsmodellen is zowel met de
dynamika (bewegingsvergelijking) als met de kontinuïteitswet
reke-ning gehouden.
De in analytische vorm gepostuleerde veldenc kunnen op twee manieren
9
-(c1) als een schematische weergave van de werkelijkheid (PASMANTER, VAN DAM)
(c2) als een aanvulling op andere velden (VAN DAM), met name als een representatie van "subgrid"-details van de op eindige roosters bere-kende snelheden van de waterbewegingsmodellen.
("Subgrid" moet hierbij niet te letterlijk worden opgevat. Het is voldoende bekend dat ook op wat grotere schalen nog wel detailstruk-turen in de op de roosters berekende snelheden kunnen ontbreken) •
(d ) Er zal ook een tweede manier van representatie van "subgrid"-snelheidsstrukturen aan de orde komen, namelijk de techniek van de randomverplaatsingen (dronkemans-stap, Monte Carlo-techniek). Men kan aan deze verplaatsingen ook een snelheid toekennen (random-snel-heid), maar deze snelheden vormen geen konsistent veld in de zin van een kontinu medium: twee deeltjes in identieke uitgangsposities kun-nen gelijktijdig totaal verschillende, bijvoorbeeld tegengestelde verplaatsingen ondergaan: algemener: er is geen korrelatie tussen de verschillende verplaatsingen/snelheden, ook als de onderlinge af-stand tot nul nadert. Hierin zit een onrealistisch element van deze techniek. TOch kan het effekt van subgridsnelheidsstrukturen er tot op zekere hoogte mee gesimuleerd worden.
IV. BANEN
Men kan deeltjesbanen onderscheiden die volledig bepaald zijn door het snelheidsveld van de vloeistof en dus samenvallen met banen van vloeistofdeeltjes, - en deeltjesbanen die daarvan afwijken. Alleen de eerste worden uitsluitend uit het snelheidsveld berekend. In het
laatste geval is sprake van niet-passief gedrag t.o.v. de vloeistof: uitwendige krachten spelen een rol.
Men kan echter nog in een ander geval van niet-passief gedrag spre-ken, namelijk wanneer de beschouwde ("vreemde") deeltjes door hun (massale) aanwezigheid de vloeistofbeweging beïnvloeden. De banen behoeven zich niet te onderscheiden van de banen van de vloeistofbe-weging, maar zouden er wel anders uitzien zonder de aanwezigheid van de vreemde deeltjes.
- 10
-In het geval dat er géén afwijking is t.o.v. de banen van de vloei-stofbeweging en er bovendien geen merkbare beïnvloeding is van de vloeistofbeweging zou men van zuiver passiefa gedrag kunnen spre-ken. De meeste sprekers van vandaag beperken zich tot dit geval. We zijn dan, mits het snelheidsveld volledig gegeven is, eigenlijk al-leen bezig met kinematika. In de beide andere gevallen komt de dyna-mika expliciet aan de orde.
Is er alleen sprake van terugkoppeling naar de waterbewegingb1, dan beperkt het dynamische probleem zich tot een bijstelling van de be-wegingsvergelijking voor de vloeistof (dichtheidsstromen, onderdruk-king van vertikale uitwisseling van impuls). Een belangrijke konse-kwentie is dat de ontkoppeling tussen waterbewegingsmodel en trans-portmodel niet meer mogelijk is. Deze ontkoppeling biedt veel voor-deel en wordt tegenwoordig dan ook toegepast waar dit mogelijk is. Het andere geval waarbij de dynamika om de hoek komt kijken is dat, waarbij, bij gegeven vloeistofbeweging, de dynamika een rol speelt bij het berekenen van de deeltjesbaan: uitwendige krachten werken op de beschouwde deeltjes op andere wijze dan op de omringende vloei-stofb2• Dit geval wordt door BOSMAN aan de orde gesteld. De kracht waar het om gaat is de zwaartekracht. Een duidelijke differentiële werking daarvan treedt pas op als het deeltje niet alleen een afwij-kende dichtheid heeft t.o.v. de vloeistof, maar ook een bepaalde mi-nimum-afmeting. Is het deeltje zeer klein, dan overheersen de inwen-dige (viskeuze) krachten.
Tenslotte kan nog het geval genoemd worden waarin differentiële wer-king van uitwendige krachten en beïnvloeding van de vloeistofbewe-ging beide een rol spelenb1+b2• Men denke aan bewegende sliblagen. Geen akademisch geval dus.
v.
DISTRIBUTIESBij afwezigheid van bronnen en putten wordt de overgang van de ene distributie in de andere, bepaald door de deeltjesverplaatsingen in het tussenliggende tijdinterval.
- 11
-In de rekenmodellen, zeker in de toegepaste sfeer, is er bij de be-rekening van de verplaatsingen die de ene distributie in de andere overvoeren in het algemeen sprake van twee afzonderlij ke verplaat-singen, korresponderend met een meer grootschalig en een meer klein-schalig deel van het veld. Het grootschalige deel wordt veelal als systematisch en reproduceerbaar beschouwd en het kleinschalige als ongeordend en niet-reproduceerbaar en daarom alleen op een statisti-sche wijze weer te geven. Men dient zich hierbij te realiseren dat de scheidslijn tussen de twee regimes in het algemeen niet met een duidelijk fysisch onderscheid in het prototype korrespondeert (zoals b.v. het geval zou zijn tussen verplaatsingen door molekulaire warm-tebeweging enerzijds en alle overige bewegingen anderzijds). De scheiding heeft daarmee uit fysisch oogpunt een willekeur ig karak-ter. Komen de snelheden uit roostermodellen, dan verschuift de grens tussen de twee regimes in de richting van de kleinere schalen als de maaswijdte van het rooster fijner wordt gekozen.
Soms is een afzonderlijke representatie van een kleinschalig deel van het veld in de modelsimulatie voor het beoogde doel overbodig. De betrekkelij ke willekeur van de scheiding tussen de twee reg imes betekent dat een stochastische modelbenadering van de kleinere scha-len meer te maken heeft met een gebrek aan deterministische kennis beneden een bepaalde schaal dan met een essentieel stochastisch ka-rakter.
Het modelmatige onderscheid tussen de beide regimes is meestal nauw verwant met het onderscheid tussen advektie- en diffusietermen bij berekeningen in termen van koncentraties.
Is een koncentratiedistributie in een gebied stationair, maar niet homogeen, dan betekent dit onder meer dat het grootschalige c.q. de-terministische snelheidsveld stationair is en bovendien dat er (in het inwendige en/of aan de randen van het gebied) zowel (één of meer) bronnen als (één of meer) "putten" aanwezig zijn.
Bij een deeltjes-representatie blijven er in een evenwichtssituatie als gevolg van het eindige aantal deeltjes fluktuaties rond de even-wichtskoncentraties bestaan.
- 12
-Tracht men een stationaire distibutie te berekenen uitsluitend met behulp van een stationair snelheidsveld dan krijgt men, hoe fijn de struktuur van het veld ook is, altijd een onrealistisch resultaat. neeltjesbanen en stroomlijnen vallen in dit geval samen. Ze verlaten of het gebied bij een put (en zijn dan bij een bron het gebied bin-nengekomen), of ze vormen een gesloten kontour • Ligt een kontinue bron van deeltjes (met een zekere afmeting) in een stroombuis m~t: een 'doorgaand karakter', dan verlaten de deeltjes bij de desbetref-fende put de uitgang het gebied. Hun koncentratie is (bij konstant deel tjesdebiet) langs de stroombuis konstant. IS de buis gesloten, dan worden de deeltjes permanent in deze buis gevangen en hun kon-centratie blijft onbeperkt 6plopen. Het onrealistische karakter van de beschreven distibuties illustreert de essentiële rol van de ver-anderlijke snelheden, ook in de natuur. Het veranderlijke element is blijkbaar altijd aanwezig (eventueel alleen in de kleine schalen) en zal in een model met meer dan één dimensie ook altijd vertegenwoor-digd moeten zijn als het om transportverschijnselen gaat. In de praktijk van het modelleren kan het inzicht wel eens worden vertroe-beld doordat bij aktuele berekening zo'n noodzakelijk mechanisme soms onbedoeld door numerieke artefakten wordt verzorgd.
Hoewel met het voorgaande een duidelijke beperking van de betekenis van stationaire snelheidsvelden gegeven is, is het toch van groot belang zich van de korrekte struktuur van deze velden (o.m. de ge-noemde eigenschappen van de stroomlijnen) te overtuigen. Anders is men al spoedig bezig bepaalde artefakten met andere modelmechanismen
te kompenseren of, erger nog, met andere artefakten. Beide leidt na-tuurlijk in het algemeen niet tot bevredigende resultaten.
Hiermee is de aandacht gerechtvaardigd die in de voordracht van
oe
Jong wordt besteed aan stationaire velden in een praktij ksituatie, i.h.b. aan de interpolatietechnieken waarmee deze velden uit gegeven snelheidmatrices worden afgeleid.- 13
-Afleiding van snelheidsvelden uit op roosters gegeven komponenten
- 15
-Afleiding van snelheidsvelden uit op roosters gegeven komponenten(*)
W.A.M. de Jong
Het model waarmee bovenstaand probleem is onder zocht is een imple-mentatie van het waterbewegingsmodel WAQUA voor het IJsselmeer. WAQUA geeft de waarden van de variabelen in de punten van een
roos-ter. Dit rooster is weergegeven in figuur 1.
n+2 + m+~ + + + -' u-snelheid.
0 0 0
,
: v-snelheid.n+1 + + + + +: waterstand.
0 0 0 0: diepte.
n + + + +
figuur 1. Het WAQUA-rooster.
Wat we zoeken, is een procedure, waarmee uit deze gegeven waarden voor elk willekeurig punt de snelheid kan worden berekend. Met name voor de berekening van stationaire concentratieverdelingen is het belangrijk een verantwoorde keuze te doen; de gekozen procedure kan grote invloed hebben op de uiteindelijke evenwichtstoestand, die met het deeltjesmodel berekend wordt.
In het deeltjesmodel wordt het gebied onderverdeeld in vakken:
+ + + + +
.
I+ + + + +
+ + + + +
+ + + + +
figuur 2. De vakken in het deeltjesmodel.
- 16
-De WAQUA-snelheden zijn dan de snelheden door de grensvlakken heen. De meest eenvoudige interpolatie is bilineaire interpolatie tussen de hoekpunten van de vakken. Van tevoren moeten dan de
snelheidskom-ponenten in deze hoekpunten worden benaderd, bijvoorbeeld met line-aire interpolatie.
û = ~(u +u )
m+~,n+~ m+~,n m+~,n+l
v
= ~ (v +v )m+~,n+~ m,n+~ m+l,n+~
met: u,v= de WAQUA-snelheden1
û,v=
de geïnterpoleerde snelheidskomponenten.Met deze methode vinden we de volgende interpolatie van een snel-heidsprofiel:
figuur 3. Bilineaire interpolatie met van tevoren middeling.
De blokken in figuur 3 geven de WAQUA-snelheden weer. ZOals uit deze figuur blijkt, worden door deze methode de snelheidsfluctuaties af-gevlakt. Nu zijn het juist de snelheidsfluctuaties, die voor de dis-persie zorgen.
We kunnen de afvlakking verminderen door direkt tussen de WAQUA-snelheden te interpoleren:
- 17
-figuur 4. Bilineaire interpolatie tussen de WAQUA-snelheden.
Dit maakt wel de interpolatie iets ingewikkelder, omdat de u-kompo-nent en de v-komponent in verschillende vakken wordt geïnterpo-leerd. Dit wordt geïllustreerd door figuur 5:
+ + + + + + +
I
,
f T •••f ...T • • • C . I...:...I + + + + +figuur 5. Interpolatievakken voor een punt C bij bilineaire interpolatie tussen de WAQUA-snelheden.
We zijn met deze laatste methode echter nog niet tevreden. We zouden namelijk van de interpolatiemethode willen eisen, dat deze de conti-nuïteitsvergelijking niet verstoort.
In het geval van een stationaire situatie is de continuïteitsverge-lijking van het WAQUA-differentieschema:
H u -H u +H v -H v =0
m+~,n m+~,n m-~,n m-~,n m,n+~ m,n+~ m,n-~ m,n-~
met: H
=
~(h +h +' +~ )- 18
-Deze vergelijking stelt de som van de debieten door de grensvlakken van een vak gelijk aan nul.
+ + + + +
n + +
Cf
+ +A
+ + + + +
rn
figuur 6. Een vak van het deeltjesmodel.
We stellen de roosterafstand gelijk aan één, zodat het debiet over AB gelijk is aan Hm+~,nUm+~,n • We eisen nu dat het debiet van de interpolatie over AB gelijk is aan het WAQUA-debiet. Verder eisen we dat de overgang in A en in B op de interpolatie over de aangrenzende vakken continu is.
Dit zijn drie eisen en dus moet het interpolatiepolynoom tenminste van tweede orde zijn. We zoeken een interpolatiepolynoom voor AB van de vorm:
(1) û(y)
=
ay2+by+cAan de continuïteit in A en B wordt voldaan door de waarde in deze punten voor te schrijven:
(2a) (2b) Û (0) Û (1) = ~(u +u ) m+~,n rn+~,n-l
=
~(u +u ) rn+~,n m+~,n+lWe zullen in eerste instantie uitgaan van een konstante diepte over de vakgrens. Dan wordt de eis m.b.t. het debiet:
H m+~,n
ofwel:
- 19
-De diepte valt dus weg uit de vergelijking.
Met de methode (1), (2) en (3) vinden we de volgende interpolatie van een snelheidsprofiel:
figuur 7. Interpolatie volgens (1), (2) en (3).
Zoals in deze figuur te zien is, kiezen we aan de randen van het ge-bied de snelheid gelijk aan de WAQUA-snelheid.
Op de methode (1), (2) en (3) hebben we een verfijning toegepast voor de punten, waarin de samenkomende parabolen afgeleiden met te-gengesteld teken hebben. In die punten vervangen we de waarden (3) door het gemiddelde van de maxima van de twee parabolen. vervolgens berekenen we de parabolen voor een tweede maal. Dit geeft als
resul-taat:
figuur 8. Interpolatie volgens (1), (2) en (3) met korrektie.
20
-Op deze manier vinden we interpolatiepolynomen voor de u-snelheden op de vertikale lijnen in figuur 9 en voor de v-snelheden op de
ho-rizontale lijnen. u u u u + + + + + __._ + + + + +
.
I + + + + +•
• • + + + + + v v v figuur 9. De polynoomlijnen.voor een punt tussen de lijnen berekenen we de snelheid met lineaire interpolatie.
We vergelijken de eerste methode (zie figuur 3) met de laatste me-thode (figuur 8) aan de hand van deeltjesbanen in een stationair snelheidsveld in het IJsselmeer. In het bijzonder bekijken
we
deel-tjesbanen, die via het Ketelmeer het gebied binnen komen. In een stationair snelheidsveld zouden deze banen via een van de uitgangen weer naar buiten moeten gaan.Figuur 10 laat een aantal deeltjesbanen zien, die berekend zijn met de eerste methode.
De deeltjes vertrekken uit het Ketelmeer en komen in de wervel A te-recht of bij B op het land. (De deeltjes kunnen over het land bewe-gen, doordat de snelheidskomponenten in een hoekpunt als B, na de middeling vooraf, ongelijk zijn aan nul. In het deeltjesmodel kunnen de deeltjes niet over land bewegen, omdat daar een reflektieprocedu-re wordt toegepast, waardoor een deeltje, dat op het land terecht komt, direkt terug in het water wordt geplaatst. Bij de berekening van de banen in figuur 10 is echter geen gebruik gemaakt van deze reflektieprocedure.)
21
-I -I
1 I
B-
.-::-
i--
,_
J
11
~c·]"
iJ
...
~-
I
I -I -I " J r--r-)
r-~r-Ir.
A,_
r--r-t; ~ir."
./ ~-
--
~ ~'(
r--(!_
\\~~
~f
~ Jr:-~
....,_,
.__~ L-I-. figuur 10. figuur 11.Figuur 11 laat een aantal deeltjesbanen zien in hetzelfde snelheids-veld als figuur 10, maar nu berekend met de laatste methode. Een ge-deelte van de banen loopt vast in kleine werveltjes1 een gedeelte loopt vast in grotere wervels als A of Br een ander gedeelte verlaat echter wel het gebied via de uitgangen. Figuur 11 laat een duidelij-ke verbetering zien ten opzichte van figuur 10.
22
-Voor een ander snelheidsveld blijkt de interpolatiemethode uit fi-guur 8 echter geen verbeter ing te geven. Figuur 12 toont de deel-tjesbanen in het betreffende snelheidsveld , berekend met de eerste methode, en figuur 13 laat de banen zien, die berekend zijn met de
laatste methode. In dit geval zijn de resultaten even slecht.
I I I (. II I figuur 12. I I I I I
I
I
r
figuur 13.23
-De oorzaak van deze slechte resultaten wordt duidelijk als we de
geometrie bij de ingang van het gebied bekijken. Figuur 14 geeft de
vier vakken bij de ingang weer:
figuur 14. Lineaire interpolatie
WAQUA-snelheden.
v
po
sitief
/ /:
U
n
egati
ef
van de
De pijlen geven de grootte van de WAQUA-snelheden aan en de
gear-ceerde gebieden geven het teken van de snelheidskomponenten aan, als
er lineair geïnterpoleerd wordt tussen de WAQUA-snelheden. Als een
deeltje het Ketelmeer verlaat, dan kan het niet verder stromen dan
de lijn AB. vervolgens beweegt het naar links en naar beneden en
komt tenslotte weer in het vak rechtsonder terecht.
De getallen bij de snelheden geven de diepten ter plaatse aan. Deze
diepten vertonen grote verschillen. De debieten staan daardoor in
een heel andere verhouding tot elkaar dan de snelheden. Figuur 15
toont dezelfde situatie, maar nu geven de pijlen debieten aan
- 24
-"\\
d,
posdiei
/ /: d
une~~
:
ltj
~ / / /// //1,.,1' I /
.
~>
.
'
"
;'tIrr«1".1/figuur 15. Lineaire interpolatie van de debieten.
Als we het deeltje de debieten laat volgen, dan stroomt het wel
ver-der het gebied in.
uit figuur 14 en figuur 15 blijkt, dat we alleen betere resultaten
kunnen verwachten als we bij de interpolatie expliciet rekening
hou-den met de diepte. Daarvoor moet dan eerst worhou-den afgesproken hoe we
de diepte zullen interpoleren. Een eerste mogelijkheid is:
bilineai-re interpolatie tussen de hoekpunten van het vak. Deze hoekpunten
vallen samen met de dieptepunten van WAQUA.
Ter bepaling van de interpolatiepolynomen stellen we weer de eis,
dat het debiet van de interpolatie over een grensvlak gelijk moet
zijn aan het WAQUA-debiet, ofwel:
1
25
-met:
n-
ht gezochte interpolatiepolynoom voor de u:H-
de interpolatie van de waterdiepte.Bovendien leggen we weer de voorwaarden (2) op. Met deze voorwaarden ligt het polynoom vast.
Verder houden we rekening met het diepteverloop bij de interpolatie tussen de polynoomlijnen. + + + + _!. i-X +
yl
è
+ + figuur 16.voor de interpolatie van de snelheid naar het punt C interpoleren we eerst het debiet naar C, waarna we dit debiet delen door de diepte:
û
=
c(1-x)B1 (y)Ol (y) + xB2(y)02(y) (1-x)B1(y) + xB2(y)
De zojuist beschreven methode voor het bepalen van het interpolatie-polynoom en de interpolatie tussen de polynoomlijnen geeft echter nogal wat problemen in de randvakken. Bovendien is de methode nogal omslachtig: eerst wordt een interpolatiepolynoom voor de snelheid omgerekend naar debieten en deze debieten worden na de interpolatie weer omgezet in snelheden.
De volgende methode is minder omslachtig en vermijdt de problemen in de randvakken. In plaats van een interpolatiepolynoom voor de snel-heid, zoeken we een interpolatiepolynoom voor het debiet:
Hu(y)= ay2+by+c
- 26 -(4) 1
rI
Hu(y)dy = H u m+~,n m+~,n Hu(0) Hu(1) = ~(H u +H u ) m+~,n-l m+~,n-l m+~,n m+~,n = ~(H u +H u ) m+~,n m+~,n m+~,n+l m+~,n+lNadat we onder deze voorwaarden het pol1noom hebben berekend, passen we weer de korrektie toe, die in figuur 8 geïllustreerd werd.
voor het punt C uit het figuur 15 bepalen
we
de snelheid, door eerst het debiet lineair te interpoleren en vervolg~ns het debiet door de diepte te delen.oe
diepte bepalen we door te interpoleren tussen de dieptewaarden uit de vergelijking (4). Voor de u-komponenten wordt dus geïnterpoleerd tussen de waarden van de waterdiepte in de omrin-gende WAQUA-punten voor de u-snelheid.Een randvak heeft een aparte behandeling nodig~
Hu
+ +
+
+ +
figuur 17.
Het debiet van de v-komponent kan gewoon lineair worden geïnterpo-leerd. Als x
<
i,
dan wordt de diepte voor de v-komponent bilineair geïnterpoleerd. Als x ~± '
dan wordt de diepte lineairgeïnterpo-leerd tussen de twee snelheidspunten A en B.
Voor de u-komponent wordt geïnterpoleerd tussen de punten
Co
en Cl· In Cl is het debiet van de u-komponent gelijk aan nul. In dat geval geeft de interpolatiemethode onrealistische resultaten. Stel bij-voorbeeld, dat de diepte in Cl' aan de rand, gelijk is aan nul. Dan vinden we voor de snelheid in C:= (l-X)HU(CO) (1-x) H (CO)
27
-In dat geval is over de gehele lijn COCl de snelheid gelijk aan de snelheid in CO. Dit is in strijd met de voorwaarde voor de gesloten rand in Cl. Daarom zullen we in een randvlak, naar de rand toe niet het debiet, maar de snelheid lineair interpoleren:
Figuur 18 toont deeltjesbanen in hetzelfde snelheidsveld als de fi-guren 12 en 13, maar nu berekend met de zojuist beschreven methode. De deeltjes stromen nu wel verder het gebied in. Een gedeelte van de banen komt nog wel in een wervel terecht, maar de rest verdwijnt, in overeenstemming met de theorie, via de uitgangen.
I
I
I
1
"""'",
r
~I I- t-I-I--n
~l...
~.J
"
figuur 18.- 28
-Conclusies.
De eerste methode, bilineaire interpolatie tussen de hoekpunten, is het eenvoudigst uitte voeren en kost het minste rekentijd. Deze methode heeft echter het nadeel dat de snelheidsfluctuaties worden
afgevlakt.
Deze afvlakking kan vermeden worden door direkt tussen de WAQtTA -snelheden te interpoleren. Deze interpolatie is iets ingewikkelder
en kost meer rekenwerk.
Voor een stationaire berekening is het belangrijk dat de
continuï-teitsvergelijking niet wordt verstoord. Bovendien dient rekening te
worden gehouden met het diepteverloop, aangezien gebleken is, dat
anders de resultaten zeer onrealistisch kunnen worden.
De methode, waarbij eerst de debiet wordt geïnterpoleerd, houdt re-kening met het diepteverloop. Verder wordt er globaal aan de conti-nuïteitsvergelijking voldaan: d.w.z.: de som van de debieten over de vakgrenzen is gelijk aan nul.
De implementatie van deze laatste methode is aanzienlijk ingewikkel-der dan die van de aningewikkel-dere methoden en het rekenen kost meer
Aspekten van tweedimensionale transportberekeningen voor oppervlaktewater door middel van passieve deeltjessimulatie
- 31
-Aspekten van tweedimensionale transportberekeningen voor
oppervlaktewater door middel van passieve deeltjessimulatie (*) G.C. van Dam
1. Inleiding
De beperking tot twee dimensies betreft in dit geval de horizontale richtingen (2 DH)1 2 DV-modellen blijven buiten beschouwing. Voor de praktische ("geografische") toepassingen wordt uitgegaan van 2DH-waterbewegingsmodellen. De daarvan afkomstige snelheidsvelden zijn wel in tweedimensionale vorm gegeven, maar ze dragen toch het stempel van een derde dimensie: van waterdiepteverschillen1 dit zijn in het niet-stationaire geval ook de veranderingen van de waterdiepte met de tijd. DOOr deze snelheidsvelden en de bijbehorende dieptevelden als uitgangspunt te nemen zijn de 2DH-deeltjesmodellen in dezelfde zin feitelijk driedimensionaal als de 2DH-waterbewegingsmodellen waarvan wordt uitgegaan.
In een aantal nevenstudies wordt zuiver tweedimensionaal ge~rkt, zo-als bij de bestudering van analytisch geformuleerde velden en tweedi-mensionale random walk-studies, c.q. kombinaties van beide. Brengt men dergelij ke zuiver tweedimensionale velden of mechanismen in ongewij-zigde tweedimensionale vorm over naar het transportmodel, dan ontstaan fouten, die veelal niet te verwaarlozen zijn (paragraaf 5). Eerst iets over de aanvullende velden en andere suppletiemethoden.
2. Random-verplaatsingen
Het meest bekende middel om berekende snelheidsvelden te suppleren is wel het gebruik van random-verplaatsingen. De gelij kwaardigheid met gradiênt-diffusie van de gewone, per stap (gemiddeld) konstante ver-plaatsing (random in richting), bij voldoend aantal deeltjes, stappen enz., is bekend. Einstein7 gaf in 1905 een bewijs voor het lD-geval1 de gevallen voor meer dimensies volgen hier in feite uit als men zich realiseert dat men de random-verplaatsingen opgebouwd mag denken uit onafhankelijke komponenten in de koBrdinaatrichtingen.
Maier-Reimer9 konstateerde dat dit mechanisme in zijn simulaties met snelheidsvelden uit 2DH-getijmodellen voor de NOOrdzee, niet toerei-kend was. In de meeste gevallen groeiden zijn deeltjeswolken (uit mo-mentlozingen) in "diameter" ongeveer met t~, dus in feite net als
wan-neer het snelheidsveld afwezig of homogeen was geweest. Slechts in
en-kele gevallen was de exponent van t iets hoger. Hij loste dit op door
de lengte van zijn random-stappen bij een momentlozill9 met diameter
nul met ta te laten groeien (t gerekend vanaf het moment van lozing).
Op deze wijze kon Maier-Reimer goede overeenstemming verkrijgen met
experimenten, zoals die van Joseph, Sendner en Weidemann8,11. Dat hij
geen behoefte had aan verdere verfijningen, komt doordat zijn wolken
niet of nauwelijks groter werden dan de maas van het rooster waarop
het waterbewegingsJROdel de snelheden berekende en toeleverde. In
fi-guur 1 ziet men een voorbeeld waarin de wolkjes waaruit de
desbetref-fende kontinue lozing is opgebouwd zelfs klein blijven ten opzichte
van de roostermaas. Dezelfde figuur laat meteen een aantrekkelijkheid
van het deeltjes-koncept zien: het snelheidsveld in het beschouwde
ge-bied is weliswaar praktisch homogeen, maar tegelijkertijd sterk
tijds-afhankelijk en deze tijdstijds-afhankelijkheid roept een gedetailleerde
sub-gridstruktuur op, z6nder dat men hiervoor eerst een lokaal deelmodel
32
-hoeft in te nesten met enkele honderden malen zoveel roosterpunten binnen eenzelfde oppervlakte.
DOordat Ma,ier-Reimer zich tot kleinschalige numer ieke exper imenten
be-perkte, bleek niet dat de inbomogenitei ten van het gegeven
snelheids-veld steeds meer tot de totale verspreiding bijdragen naarmate het
diffunderende systeem groter is. Daar het in de toepassingen bij de
Rij kswaterstaat zeer duidelij k de bedoeling is de schaal van één of
twee maas-lengten royaal te passeren, moet met het genoemde effekt
re-kening gehouden worden. Voor de wijze waarop dit gerealiseerd is in de
modellen waarin met random-verplaatsingen wordt gewerkt, wordt
verwe-zen naar lit. 3 (N.B. de beschouwing aldaar is strikt tweedimel"sio
-naal) •
Hoewel de lijn van Maier-Reimer dus wel is opgepakt en ·verlengd",werd
tegelijkertijd het onbevredigende ervan beseft. DOOr de
random-ver-plaatsingen van de tijd te laten afhangen, kennen we de deeltjes een
geheugen toe. Fysisch is het konoept dus ongerijmd. Men kan zich
mis-schien voorstellen of wensen dat dit niet tot (merkbaar) onjuiste
re-sultaten behoeft te leiden, maar werkelijk korrekt kan het model toch
nooit zijn. vooral bij kontinue lozing, waarbij elk deeltje volgens
dit koncept een andere leeftijd heeft en de verschillende
leeftijds-groepen zich ook met elkaar kunnen vermengen, vraagt men zich af
hoe-ver men kan gaan. In de praktij k zal het vaak meevallen, als gevolg
van het afremmen en vervolgens tot stilstand komen van het aangroeien
van de random-stap, beschreven in lito 3.
De fysische ongerijmdheid van het koncept heeft het zoeken naar
alter-natieven gestimuleerd. Op deze wijze ontstond de interesse voor
analy-tisch geformuleerde kunstmatige snelheidsvelden.
3. Kunstmatige snelheidsvelden
In figuur 2 ziet men een aanschouwelijke voorstelling van het
eenvou-digst denkbare tweedimensionale stationaire "wervelveld" (zie de
for-mules). Aan de kontinuïteitsvergelijking voor het zuiver
tweedimensio-nale geval (waarin verschillen in waterdiepte of waterstand niet
be-staan) is voldaan. Wil men de zijden van de vierkanten géén hoek laten
maken met de koördinaatassen, dan heeft men per komponent twee
harmo-nische funkties nodig (figuur 3). Het wervelveld heeft overigens
pre-cies dezelfde eigenschappen. De stroomlijnen c.q. deeltjesbanen binnen
de vierkanten hebben dezelfde vorm (figuur 4).
De getoonde velden zijn niet alleen sterk gestyleerd; ze bevatten
bo-vendien slechts één ·golflengte". De eenvoudigste manier om tot een
"spektrum" te komen is een aantal van deze velden, met verschillende
golflengten en willekeurige fasen (ruimtelijk; het veld is nog steeds
stationair gedacht) gewoon lineair bij elkaar op te tellen. Hiermee
wordt niet bereikt dat kleine wervels in grotere worden meegevoerd,
zoals in de natuur wel gebeurt. Een dergelijke situatie is trouwens
onverenigbaar met stationariteit. De sommatie levert meer
gekompli-ceerde velden op, zoals in figuur 5 weergegeven in vektorvorm en in
samen-- 33
-vallend met deeltjesbanen) . Er zijn voor figuur 6 ongeveer 60 deeltjes losgelaten, zodanig dat het veld min of meer gelijkmatig werd gevuld, doch overigens willekeurig. Hierbij valt op dat slechts twee van de 60 banen een groter gebied omspannen dan overeenkomt met de grootste golflengte van de deelvelden: mogelijk gaan deze banen zelfs naar on-eindig. Er bestaat grote twijfel over de vraag of dergelijke banen in-derdaad v66rkomen of dat, ondanks de zeer kleine tijdstap waarmee ge-rekend is (waardoor de overige banen netjes sluiten) hier toch nog van een incidenteel numer iek artefakt sprake is. Hiervoor bestaat naast het incidentele voorkomen van deze banen in het rekenresultaat nog een tweede aanvi jZlng. Zodra men namelij k de rekenstap iets te groot neemt, zodat de deeltjes in de berekening uit de in werkelijkheid ge-sloten banen langzaam naar buiten "spiralen", verzamelen ze zich in banen die zich op dezelfde wijze tussen de grootste "wervels" d66r be-wegen als de twee bijzondere banen van figuur 6. Deeltjes van uiteen-lopende herkomst komen in deze "straten" samen en blijven er. De schijnbare aantrekking is in feite een afstoting door de wervels ter weerszijden: de numerieke afwijking heeft een "middelpuntvliedend" ka-rakter.
uitsluitsel over het mogelijke bestaan van grotere banen in velden van dit type is nog niet verkregen, hoewel externe deskundigheid is inge-roepen. uit de illustratie is echter wel duidelijk dat een deeltjes-verzameling van eindige afmetingen, zeker wanneer deze afmetingen kleiner zijn dan die van de grootste wervel, in het algemeen opgeslo-ten zal blijven binnen een eindig gebied (bij voorbeeld in de "ring" tussen twee "koncentrische" stroomlijnen), zolang het veld "bevroren" blijft. Dit is bij de numerieke experimenten ook bevestigd. Zodra er
echter een veranderlijk element wordt ingebracht dat zich overal in
het veld manifesteert, breidt een deeltjesdistributie zich op den duur onbegrensd in alle richtingen uit. Eén van de manieren om dit te rea-liseren is een kleinschalige random walk, met verplaatsingen, kleiner
dan de kleinste wervel in het spektrum (figuur 7, overgenomen uit
lito 1). Bij gebruikmaking van een kleine random walk kan men in het numerieke model alle deeltjes bij een momentlozing zelfs in één punt loslaten, iets wat fyisch uiteraard orunogelijk is. Vervangt men de
random walk door een zekere verander lijkheid in één of meer van de
deelvelden (golflengten), dan ontstaat eveneens een onbeperkt
voort-gaande verspreiding. uiteraard moet dan de deeltjesverzameling een
eindige begin-afmeting hebben, wat ook fysisch realistisch is. Men be-schrijft op deze wijze het dispersieproces geheel door advektie, zon-der dat hierdoor vanuit "makroskopisch" gezichtspunt essentiële
ele-menten verloren gaan. Bovendien kan men door juiste keuze van het
spektrum in principe elk schaaleffekt in het dispersieproces
verkrij-gen dat men in het prototype aantreft als gevolg van de spektrale
struktuur die dáár heerst. Dit laatste zou voor de groei van
individu-ele wolken 66k nog wel moge lijk zijn met een gemodificeerde random
walk-methode (namelijk door in de faktor ta de exponent a op meer in-gewikkelde wijze van de tijd te laten afhangen dan tot op heden), maar
de onrealistische elementen van die methode vervallen daarmee niet.
Dit komt niet alleen tot uiting bij kontinue lozing. Bij de advektieve methode verkr ijgt men gr illiger vlekvormen die voortdurend verande-ren. Een zichtbaar effekt bij kontinue lozing is het meanderen van de
pluim als men de advektieve methode (kunstmatig wervelveld) toepast
(figuur 16).
Het spektrum van de kunstmatige wervels is uit een eindig aantal dis-krete golflengten opgebouwd ("lijnenspektrum")• De dispersie blijkt betrekkelijk ongevoelig voor de onderlinge spatiëring naar golflengte
- 34
-"energie") over de golflengten. Een nauwkeurige analyse van de invloed van de spatiëring kon nog niet worden voltooid, omdat dit onderzoek nogal reken-intensief is. Wel bleek dat bij verschillende systemen van spatiëring toekenning van gelijke snelheden aan alle golflengten steeds bij goede benader ing het bekende resultaat oplevert van een maximumkoncentratie die met t-2 afneemt, c.q. een "diameter" die groeit met t1 (als t
=
0 het tijdstip is waarop de diameter ongeveer nul is). DOor een zekere progressie in de snelheid aan te brengen werd een meer progressieve dispersie verkregen. Interessant is ook dat de door verschillende theorieën voor speLde afwijking van de ruimtelijke verdeling t.o.v. de normale distributie, uit de simulatie duidelijk te voorschijn komt en de voorspelde richting heeft (figuur 7, inzet). Een aanschouwelijke voorstelling van het advektieve dispersieproces wordt verkregen door bij de simulaties zonder random walk niet met een min of meer willekeurige deeltjesverzameling te beginnen, maar met een ver zameling die geheel op een cirkelomtrek is gelegen, een methode waarmee bij de Fysische Afdeling beg in 1983 is begonnen. En passant ziet men daarbij meteen wat er met de "watermassa" binnen de kontour gebeurt. Zolang het systeem nog klein is t.o.v. de grootste wervel van het gekozen spektrum, maakt het niet veel uit of men werkt met een stationair veld of met veranderlijke komponenten. Immers, aanvankelijk verplaatsen de grotere wervels het systeem voortdurend naar andere de-len van het veld, wat hetzelfde effekt sorteert als een veld dat (te-vens) in de tijd verandert. Zo is figuur 8 met een stationair veldbe-rekend, waarbij de grootste golflengte gelijk is aan ca. 20 keer de diameter van de begincirkel, dus nog groot is ten opzichte van de li-neaire afmetingen die op het laatste tijdstip ("t
=
20") bereikt zijn. Men kan de in beeld gebrachte ontwikkeling dus in feite beschou-wen als de beginfase van een "algemeen" geval (met tijdafhankelij k veld) •De lineaire afmetingen van zoln systeem nemen op den duur onbegrensd toe: de kontour , begonnen als cirkelomtrek, wordt oneindig lang en neemt in lengte veel sneller toe dan de "overall·-afmeting van het systeem. Het omsloten oppervlak verandert echter niet: desalniettemin komt de kontour uiteindelij k "overal dicht" te liggen, evenals de deeltjes die zich binnen de kontour bevinden, hoewel een geheel binnen de kontour gelegen oppervlakje altijd hetzelfde aantal deeltjes blijft bevatten: de ·koncentratie" binnen de kontour blijft altijd dezelfde, al wordt dit, wanneer men het inwendige van de cirkel bij de aanvang gelijkmatig met een eindig aantal deeltjes zou beleggen, op een gege-ven moment wel een onmeetbare grootheid: de diameter van alle cirkel-tjes die men nog binnen de kontour kan tekenen nadert tot nul en daar-mee ook het omsloten aantal deeltjes per cirkeltje bij eindig totaal aantal deeltjes. Ondanks het konstant blijven van de deeltjesaantallen per oppervlakte-eenheid, nemen de afstanden tussen deeltjes langs kon-touren enorm toe: er vinden geweldige uitrekkingen plaats waardoor al gauw op sommige plaatsen wegens het eindige aantal deeltjes de kontour moeilijk te rekonstrueren is, terwijl op andere plaatsen de uitrekking nog niet of nauwelijks zichtbaar is. SOms lopen twee kontourfragmenten met een vrij grote en een relatief kleine onderlinge deeltjesafstand op korte afstand evenwijdig aan elkaar (figuur 9). Naast verspreiding is ook sprake van een zekere ordening. Iedere deeltjesverzameling bin-nen de begincirkel is immers veroordeeld om binnen de kontour te blij-ven zoals deze zich verder ontwikkelt. omgekeerd kan men rond iedere
35
-ver zameling van eindige afmet ingen een gesloten kontour trekken. Men ziet dan ook, dat wanneer het snelheidsveld wordt losgelaten op een deeltjeswolk met bij voorbeeld een gauss-achtige distr ibutie, deze wolk na enige tijd de gedaante krijgt van een veelpotige spinachtige draadfiguur (figuur 10).
De berekeningen voor figuur 8 werden tot en met het voorlaatste tijd-stip (nt
=
10") uitgevoerd met de UNIVAC 1100 met 2000 deeltjes. Met dit deeltjes-aantal zijn sommige kontourfragmenten voor tijdstip 10 al moeilijk te konstrueren. Met een CRAY (te Reading) werd de berekening herhaald met 10000 deeltjes, waarmee een redelijke rekonstruktie voor t=
10 geen problemen meer oplevert. De berekening werd met dit deel-tjesaantal doorgezet tot t=
20. Hoewel een zeer illustratieve figuur werd verkregen, levert de rekonstruktie van bepaalde kontourgedeelten grote problemen op, groter dan bij 2000 deeltjes voor t=
10. Het aan-tal benodigde tijdstappen hangt af van het gebruikte rekenschema. Met een tweede-orde-schema kan voor het tijdstip t=
20 met 100 tijdstap-pen worden volstaan. Voor deze 100 stappen van 10000 deeltjes is op de CRAY ca. 13 s rekentijd nodig (gevektoriseerde versie van het program-ma). Daarbij bedroeg het aantal harmonische komponenten voor beidekoördinaatrichtingen tien.
Berekeningen om het gedrag van dichtheid of omvang van deeltjeswolken in de tijd te onderzoeken (met een tijdafhankelijk element in het snelheidsveld of een ander "ovarstap=-mechan Ieme) strekken zich uit over veel meer tijdstappen (men wil in beginsel doorrekenen tot de deeltjeswolk groot ten opzichte van de grootste wervel is).
om
tot aanvaardbare rekentijden te komen wordt met kleinere aantallen deel-tjes gewerkt dan wanneer men kontouren zichtbaar wil maken.Het is interessant om na te gaan in hoeverre eigenschappen van onze in eerste instantie als kleinschalig bedoelde kunstmatige velden, zoals deze tot uitdrukking komen in hun werking op cirkelvormige kontouren, op grotere schaal zijn terug te vinden in de met 2DH-waterbewegingsmo-dellen verkregen velden. We kunnen eerst kijken naar de geschiedenis van een aantal kontouren in een snelheidsveld verkregen uit een model
van de zuidelijke Noordzee met een maaswijdte van 10 km (figuur
17a). De diameter van de begincirkels is 7600 m, na verloop van tijd hebben sommige kontouren een grootste lineaire afmeting van ruim 4 maaslengten bereikt, maar er gebeurt verder heel weinig mee in verge-lijking met wat we zagen in de kunstmatige velden, hoewel is doorgere-kend over een tijdsduur van 240 getijperioden (M2) en er niet met tijdsgemiddelde velden maar met de volledige getijbeweging is ge-werkt. Interessanter is nog de vergelijking met het gedrag van kontou-ren met diameters van de aanvankelijke cirkels van 800 en 4000 m (2 en 10 keer de maaswijdte van het waterbewegingsmodel) in snelheidsvelden uit het WAQUA-model van de Westerschelde (figuren 11 en 12). De ver-vorming van de kleinere cirkel komt wat langzamer op gang dan die van de grotere, maar in beide gevallen verloopt het proces razendsnel in vergelijking met het Noordzeemodel (zie de aangegeven aantallen getij-perioden) • Het verschil in homogeniteit tussen de beide (tijdafhanke-lijke) velden is opvallend en is niet louter toe te schrijven aan het verschil in maaswijdte tussen de beide modellen: ook wanneer men reke-ning houdt met de maaswijdte, is de stroomsnelheidsverdeling in grote
36
-delen van de Noord zee aanzienlij k vlakker dan in het estuar ium met zijn geaccidenteerde bodem en grillige begrenzing.
Beperken we ons verder tot het estuarium, dan blijkt uit de figuren duidelijk hoe groot de rol is van het berekende snelheidsveld in het totale dispersieproces. Reeds na 8 getijperioden (M2) zijn de grootste lineaire afmetingen van de (als kleine cirkel begonnen) deeltjesverza-meling wanneer alleen het berekende snelheidsveld gewerkt heeft, niet veel kleiner meer dan wanneer ook de suppletie voor de kleine schalen is toegepast (figuur 13). Aan de andere kant is de struktuur van de verzameling in het ene geval toch wel zo sterk verschillend van die in het andere geval, dat dezelfde illustratie tevens kan dienen om de be-langr ijke rol van "de kleine schalen" voor de "doormenging" te den.on -streren.
Dat sterke vervorming in tamelijk korte tijden ook door een
tijdonaf-hankelijk veld kan worden veroorzaakt moge blij~en uit enkele
resulta-ten verkregen met een berekend snelheidsveld voor het IJsselmeer
(fi-guur 14).
Dat een waterbewegingsmodel met een kleinere maaswijdte een grotere
neigen bijdrage" aan de dispersie geeft, wordt al enigszins
geïllu-streerd door vergelijking van de figuren 15a en 15b verkregen met
mo-dellen met maaswijdten van respektievelij k 2 en 1 kilometer. In het
fijnere model is in dit geval met dezelfde toevoeging aan het
bereken-de snelheidsveld gewerkt als in het grovere, met het gevolg dat in de
afgebeelde periode een grotere totale dispersie optreedt in het
fijne-re model. De gebruikte parameter-waarden zijn die, welke bij de
afre-geling van het grove model bepaald zijn door vergelijking met
waarne-mingen. Door de parameters niet bij te stellen verloopt de dispersie
in het fijnere model dus sneller dan met deze waarneningen
overeen-komt.
4. Kombinatie van kunstmatige velden met berekende
Men zou zich kunnen voorstellen dat de kunstmatige snelheidsvelden
weliswaar uit fysisch oogpunt de voorkeur verdienen als "subgrid"-aan-vulling op dynamisch berekende snelheidsvelden, maar dat de toepassing
uit een oogpunt van rekentijd bezwaarlij k zou zijn. Erhter blij kt in
de praktij k dat men met een klein aantal golflengten kan volstaan
(aangevuld met een kleine random-stap als representant van de
aller-kleinste schalen). De praktische toepassing van de kunstmatige velden
is tot dusverre om geheel andere redenen beperkt gebleven. In het
ge-geven (berekende) snelheidsveld komen meestal snelheden voor die een
orde groter zijn dan de benodigde snelheden in de toe te voegen
wer-vels. Geeft men deze laatste gefixeerde posities, dan wordt hun
dis-persieve werking onderdrukt doordat de snelheid van het gegeven veld
de deeltjes zo snel door de wervels heen "trekt", dat hun invloed op
de deeltjesbaan sterk vermindert. Is het gegeven veld in het
interes-segebied, bij voorbeeld een zekere omgeving van een bron, redelijk ho-mogeen, dan is het aangeduide probleem op te lossen door het
wervel-veld met de (globale) snelheid van het grootschalige veld in het
inte-ressegebied te laten meebewegen • Deze oplossing is eenvoudig en is dan
ook in de programmatuur voorzien en als zodanig toegepast. Een
resul-taat voor een beperkt gebied (in het IJsselmeer) met een betrekkelijk
redelij 37 redelij
-ke homogeniteit niet meer voldaan, dan blijkt inderdaad dat de disper-sie merkbaar van de plaats gaat afhangen door het van gebied tot ge-bied variërende snelheidsverschil tussen het gegeven veld en het aan-vullende veld. Men zou het aanvullende veld dus overal willen laten bewegen met de lokale snelheid. Dit is een eis waaraan met de ruimte-lijk "starre" velden die in het voorafgaande ten tonele gevoerd werden (ook waar over de tijdafhankelij ke var iant werd gesproken) niet kan worden voldaan. Een oplossing 1 ijkt mogel ijk door de reeds genoemde tijdafhankelij kheid door middel van wervels van eindige levensd uur
(wat ook fysisch realistisch is), te kombineren met onderlinge onafhankelij kheid van wervels, d .w.z. de wervels ontstaan lokaal en worden door het gegeven veld (dus onafhankelijk van elkaar) geadvekteerd. Ook dit laatste is fysisch realistisch. Deze en andere
oplossingen worden nader onderzocht.
5. Driedimensionale aspekten
We komen nu terug op de kwestie dat er weliswaar tweedimensionaal wordt gerekend, maar een driedimensionaal systeem wordt beschreven. De verschillen in waterdiepte drukken hun stempel op de aangeboden snel-heidsvelden.
om
uit de tweedimensionale deeltjesdistributies (3D-)kon-centraties te berekenen moeten de deeltjesaantallen per oppervlak-te-element door de lokale en momentane waterdiepte worden gedeeld. Een zuiver tweedimensionaal ontworpen mechanisme, zoals 2D-random-walk past niet in dit wezenlijk driedimensionale geheel. Als de dieptever-schillen in het interessegebied gering zijn kunnen een 2D-random-walk of een snelheidsveld dat gebaseerd is op een zuiver tweedimensionale kontinuïteitswet zonder bezwaar worden toegepast. Een relatief vlakke bodem zoals die van het IJsselmeer vertoont echter reeds te grote diepteverschillen om zonder korrekties een 2D "diffusie"-mechanisme te kunnen toepassen. Voor een stationaire koncentratieverdeling leidt dit tot koncentratieverschillen van dezelfde relatieve grootte als de diepteverschillen • Immers, een 2D diffusiekoncept leidt in situaties waar homogenisering optreedt tot een gelijkmatige 2D-verdeling: gelij-ke aantallen deeltjes binnen gelijgelij-ke oppervlakgelij-ken. De (3D)koncentratie is in zo'n geval dus omgekeerd evenredig met de waterdiepte, in strijd met de realiteit.Het is uiteraard niet toelaatbaar om voor deze afwijking een korrektie achteraf toe te passen, bij voorbeeld door de volumekoncentratie maar evenredig aan de oppervlaktekoncentratie te nemen. Men kan geen ver-trouwen hebben in een dergelij ke eindsituatie als men "tussentijds" steeds met onjuiste verdelingen heeft gewerkt. De korrektie zal moeten worden aangebracht in de simulatie van het proces van het uitwisse-lingstransport. Wanneer de "diffusie" wordt bewerkstelligd door middel van stappen ter grootte ds (2 D hor., random in richting) dan moet elk deeltje bij elke tijdstap een extra verplaatsing a ds2
4H
in hellingafwaartse richting gegeven worden om het gewenste 3D-gedrag te krijgen (a= absolute waarde van de lokale gradient van H; H= lokale hoogte van de waterkolom). Voor de afleiding van deze korrektie wordt verwezen naar lit. 4, evenals voor de overeenkomstige procedure als i.p.v. random-verplaatsingen een supplementair snelheidsveld wordt toegepast •
Behalve bij geringe diepteverschillen is een 2D-aanpak van het uitwis-selingstransport wellicht ook toelaatbaar in gevallen waarin de
deel-- 38
-tjesdistr ibutie zeer snel var ieert met de tijd, althans ten opzichte van de geometrie van het gemodelleerde gebied. Men kan hierbij denken aan de momentlozingen in het "kalamiteitenmodel" voor de Westerschel-de5,6. De overweging is, dat voor een bepaald waterpakketje (en de
zich daarin bevindende deeltjes) grootte en richting van de bodemhel-ling zodanig varieert in de tijd, dat de zojuist beschreven fouten el-kaar voor een groot deel opheffen.
6. Artefakten door middeling van snelheden in ruimte en tijd
Voor problemen in samenhang met ruimtelijke middeling c.q. interpola -tie wordt verwezen naar de bijdrage van W.A.M. de Jong. De problemen vloeien geheel voort uit de diskretisering, de fouten naderen tot nul als de maaswijdte van de modellen tot nul nadert.
Voor problemen samenhangend met middeling of interpolatie in de tijd geldt in wezen hetzelfde. Ze worden bestreden door de tijdstap vol-doende klein te nemen en/of hogere-orde-rekenschema's toe te passen. Echter, wanneer het snelheidsveld periodiek is en een groot aantal pe-rioden moet worden doorgerekend, bestaat vaak de wens om de netto-ver-plaatsing over een periode te bepalen met behulp van één snelheids-waarde voor de gehele periode. Dit is in feite alleen mogelijk als het snelheidsveld homogeen is. IS dit niet het geval, dan is er sprake van een benadering. uit recent onderzoek2 is gebleken dat benadering met behulp van zogenaamde Eulerse grootheden reeds voor relatief homogene getijstroomvelden tot ontoelaatbaar grote fouten kan leiden (fig.
17). Tracht men dit te ondervangen door voor een gegeven periodiek veld (zoals een getij stroming, eventueel met een stationaire windin-vloed) de periodegemiddelde verplaatsingen a.h.w. eenmalig uit te re-kenen, dan blijft men toch nog met een probleem zitten: de verplaat-singen hangen af van de fase waarop de deeltjes uit de gekozen punten vertrekken. In lito 2 worden voorbeelden gepresenteerd (fig. 18) waar-in deze fase-afhankelijkheid klein is t.o.v. de afwijkingen die ont-staan bij gebruik van Eulerse grootheden. Dit behoeft echter niet te betekenen dat de fase-afhankelijkheid in het algemeen een v erwaarloos-baar effekt is.
7. Technische bijzonderheden
voor een bespreking van een aantal min of meer technische zaken zoals nauwkeurigheid, rekentijden, speciale technieken (m.n. superpositie), randbehandeling, afregeling van rekenmodellen, afbraak, reakties, na-bewerking en presentatie van rekenresultaten, wordt verwezen naar de bijdrage over 2 DH-deeltjesmodellen aan de in september te houden "WAQUA-vervolgkursus"4.
- 39
-8. Literatuur
1. Dam, G.C. van, Distinct-particle simulations.
Paragraaf VIII van Chapter 2 in .R:>llutant transfer and transport in the sea, G. Kullenberg ed., vol. I.
CRC Press Inc., Soca Raton, Florida, 1982.
2. Dam, G.C. van, Reststromen en resttransport in modelberekeningen. Rijkswaterstaat, Directie Waterhuishouding en Waterbeweging, Fys
i-sche Afdeling, nota FA 8402, 1984.
3. Dam, G.C. van, Konstante en aangroeiende random-stap in deeltj es-modellen.
Rijkswaterstaat, Directie Waterhuishouding en Waterbeweging, Fy si-sche Afdeling, notitie 85-FA-219, 1985.
4. Dam, G.C. van, Deeltjesmodellen(tweedimensionaal-horizontaal). Bijdrage aan (syllabus bij) de kursus Waterkwaliteitsmodelering in relatie tot WAQUA-gebruik.
Rijkswaterstaat en Waterloopkundig Laboratorium, september 1985 (Tevens afzonderlij k als nota FA 8507, Rij kswaterstaat, Directie Waterhuishouding en Waterbeweging, Fysische Afdeling).
5. Dam, G.C. van, A particle model applied to the Western Scheldt. Liverpool Conference on Mixing and Dispersion in Estuar ies, 24/25 September 1985 (Abstract: notitie 85-FA-226, Rijkswaterstaat, Di-rectie Waterhuishouding en waterbeweging , Fysische Afdeling, 1985) •
6. Dekker, L., en F.O.B. Lefèvre, Presentatie calamiteitenmodel
Westerschelde.
Rijkswaterstaat, Directie Waterhuishouding en waterbeweging , Di s-trict Kust en Zee, Adviesdienst Vlissingen, notitie WWKZ-85.V300,
1985.
7. Einstein, A., Ober die von der molekular kinetische Theor ie der Wärme geforderte Be wegung von in ruhenden Flüssig keiten suspen-dierten Teilchen.
Ann. Phys. (Leipzig), 4, 17, 1905.
8. Joseph, J., H. Sendner, en H. Weidemann, Untersuchungen liberdie horizontale Diffusion in der Nordsee.
Deutsche Hydrographische Zeitschift, 17, 2, 1964.
9. Maier-Reimer, E., Hydrodynamisch-numerische Untersuchungen zu ho-rizontalen Ausbreitungs- und Transportvorgangen in der Nordsee. Dissertatie, Mitt. Inst. für Heereskunde der Universität Hamburg, XXI, 1973.
10. Pasmanter, R.A., Dynanical systens, determinist ic chaos and dis-persion in shallow tidal seas.
In: Berekeningen van banen en distributies van deeltjes in
twee-en
driedimensionale snelheidsvelden. Colloquiumdag (19 juni 1985)Rijkswaterstaat, Directie Waterhuishouding en Waterbeweging ,
Fy-sische Afdeling (nota FA 8506) in samenwer king met OVerleggroep Transportverschijnselen van de Raad van OVerleg voor het fysisch-oceanografisch onderzoek van de Noordzee (Tevens afzonderlijk als
nota FA 8505, RWS,
ww,
Fys. Afd.), 1985.11. Weidemann, H. (editor), The ICES diffusion experiment RHENO 1965. Conseil International pour l'EXploration de la Her, Charlottenlund Slot, Danemark; Rapports et Procès Verbaux des Réunions, vol. 163,
- 41
-Abb.ZO.Z Transport und Diffusion naeb Monte-Carlo-Verfahren aua konstanter Quelle bei MZ-Gezeit und 14 m/.ee Nord-West-Wind über der Nords.e in ZO Perioden P-O .Q-2
a sin
(
2~
X
+~X)
cos(
2~
t
+W
y)vy =
-S
a cos(2ïTX
+'l(i
)
,
(2ïT
)
Isotrope geval (t=l)À x sin ~ +
