De turbulente stroming over de zeereep

HILL MDL 2.0

DE T U R B U L E N T E S T R O M I N G

OVER DE Z E E R E E P

F

LANDSCAPE A N D ENVIRONMENTAL RESEARCH GROUP UNIVERSITEIT V A N A M S T E R D A M

Rapport nr:

Titel en subtitel:

HILL MDL 2.0

De turbulente stroming over de zeereep

Schrijvers:

J.H. van Boxel, S.M. Arens

Naam en adres opdl'achtnemer:

Vakgroep Fysische Geografie en Bodemkunde Universiteit van Amsterdam

Serie nr.:

Nieuwe Prinsengracht 130, 1018 VZ Amsterdam

Naam en adl'es opdl'achtgever:

Rijkswaterstaat, Dienst Weg- en Waterbouwkunde tbv TAW-C

Postbus 5044, 2600 GA Delft

Opmel'ldngen:

Referaat:

Ontvanger catalogus nummer:

Datum Rapport: Maart 1997

Kode uitvoerende Organisatie: FGBL-UvA

Nr. Rapport uitvoerende organisatie:

Projektnaam:

Fysische Modellering zeereep-ontwikkeling

Kontraktnummer:

Type rapport: Project rapportage

Kode andere opdrachtgever:

Doel van het project "Fysische modellering zeereepontwikkeling" is de ontwikkeling van een model waarmee effecten van beheersmaatregelen op de morfologische ontwikkeling van de zeereep kunnen worden voorspeld, teneinde daarmee het zeereepbeheer te optimaliseren. Dit model bestaat uit een tweetal deelmodellen. Een stromingsmodel berekend de luchtstroming over de zeereep en de effecten van topografie en vegetatie op deze turbulente stroming. Een zandtransportmodel berekend op grond van de resultaten van het stromingsmodel de mate van zandtransport, de plaatsen waar erosie en waar depositie optreedt en de nieuwe vorm van de zeereep. Ais de vorm van de zeereep significant veranderd wordt het stromingsmodel opnieuw aangeroepen.

Op grond van een eerdere haalbaarheidsstudie was aan het licht gekomen dat het stromil)gsmodel op een aantal punten verbeterd moest worden: Het model diende niet meer vast te lopen op recirkulatiewervels, de snelle redistributie van turbulente kinetische energie moest verbeterd worden, een fout in de berekening van de laterale wind moest verwijderd worden, de effecten van de kromming van de stroomlijnen op de turbulentie moesten opgenomen worden, de uitvoer moest geschikt gemaakt worden voor het SCOPE-model en het interpolatiealgoritme dat gebruikt werd voor de ruwheidslengte moest gewijzigd worden. Bovendien diende een validatie plaats te vinden op grond van meetgegevens.

Aile genoemde wijzigingen zijn uitgevoerd en hebben de betrouwbaarheid van het model aanzienlijk verbeterd. Een gevoeligheidsstudie laat zien dat het model voor een aantal parameters gevoelig is. Het zou aan te raden zijn een meer gedetailleerde gevoeligheidsstudie uit te voeren. Bij deze studie bleek de impulsflux op enige hoogte (tussen 10 en 1000 m) onderschat werd door het model. De oorzaak hiervan is zeer waarschijnlijk een procedure voor het vereffenen van het impulsfluxprofiel. Dit zal moeten worden gewijzigd.

Bij de validatie op grond van over Nederlandse zeerepen gemeten windprofielen bleek dat de windprofielen aan de zeewaartse zijde en op de top van het duin redelijk worden benaderd. Het is niet duidelijk in hoeverre afwijkingen te wijten zijn aan de onjuist gemodelleerde impulsflux. De vertraging van de stroming aan de lijzijde van de zeereep wordt door het model onderschat. Hieraan is gezien de geringe relevante voor het berekende zandtransport weinig aandacht besteed.

De gevoeligheidsstudie zien dat voor een goede validatie van het model ook gegevens nodig zijn omtrent de impulsflux of de wrijvingssnelheid, daar deze laatste de belangrijkste invoerparameter is voor de berekening van het zandtransport.

Trefwool'den:

SAFE, zeereep, model, wind, turbulentie, zandtransport

Classificatie: Classificatie van dcze pagina:

Distributiesysteem:

INHOUDSOPGAVE

V O O R W O O R D 3

1 I N L E I D I N G 5

2 T H E O R I E 7 2.1 Effecten van overgangen en topografie op het windprofiel 7

2.2 Reynoldse splitsing 9 2.3 De hydrodynamische drukverstoring 10

2.4 Vergelijking voor de horizontale wind I I 2.5 Tweede orde vergelijkingen voor de varianties 12

2.6 De vergelijking voor de impulsflux 16

2.7 De dissipatievergelijking 17 2.8 Aannamen voor sluiting van de 2''''-orde vergelijkingen 17

2.9 Randvoorwaarden 18

3 M O D E L B E S C H R I J V I N G 21

4 A A N P A S S I N G E N I N H E T M O D E L 25

4.1 Recirculatiewervels 25 4.2 Redistributie binnen één simulatiestap 25

4.3 Fout laterale wind 26 4.4 Effecten van de kromtestraal op de stroming 28

4.5 Invoer en uitvoer geschikt maken voor S C O P E model 30

4.6 Slingering ruwheidslengte, Zg 30

4.7 Geluidssignaal 31 4.8 Schaling initiële profielen 31

4.9 De parameterfile H I L L P A R . P A R 31

5 G E V O E L I G H E I D S A N A L Y S E E N P A R A M E T E R I S A T I E 33 5.1 Het Zeman en Jensen model en Askervan Hill metingen 33 5.2 Vergelijking met Z&J87 bij verwaarlozing van krommingseffecten 35

5.3 Het windprofiel met inbegrip van krommingseffecten 37

5.4 Gevoeligheid voor de parameter Delta-Alfa 37 5.5 Gevoeligheid voor windsnelheid; schaling 39

6 V A L I D A T I E 41 6.1 Vergelijking met de metingen in Nieuw-Haamstede 41

6.2 Vergelijking met de metingen in Groote Keeten 42 6.3 Vergelijking met de metingen op Schiermonnikoog 44

6.4 Conclusie validaties 46

7 C O N C L U S I E S E N A A N B E V E L I N G E N 47

8 V E R W I J Z I N G E N 49

B I J L A G E N : 1 H O O F D P R O G R A M M A V A N H E T H I L L M O D E L 51 2 D E P A R A M E T E R F I L E H I L L P A R . P A R 53 3 D E C O N F I G U R A T I E F I L E H I L L _ C F G . C F G 55 4 D E P R O C E D U R E H I L L F L O W 57 5 P R O C E D U R E S V O O R T E R M E N 7^' O R D E V E R G E L I J K I N G E N 59 6 V O O R B E E L D V A N D E I N V O E R F I L E S 61 7 V O O R B E E L D V A N D E U I T V O E R F I L E S C O P E I N . P R N 63 8 V O O R B E E L D V A N D E F I L E I N I T I A L . I N 65 9 B E R E K E N I N G V A N D E K R O M T E S T R A A L 67

VOORWOORD

Het voor u liggende rapport vormt onderdeel van een serie rapporten die zijn verschenen in het kader van het project "Fysische Modellering Zeereepontwikkeling. De project streeft naar de ontwikkeling van een simulatiemodel SAFE (Simulation of Aeolian Foredune Evolution", een model dat ontwilckeld wordt ter ondersteuning van zeereepbeheer ders.

In dit rapport wordt het deelmodel "H1LL_MDL 2.0" besproken dat de luchtstroming over het duin simuleert. Dit stromingsmodel, gebaseerd op fysische principes, is een 2-dimensionaal 2''^ orde sluitingsmodel (StuU 1988). De benodigde theorie is voor een groot deel ontleend aan Zeman en Jensen, 1987. In een eerder project was reeds een voorlopige versie van dit model geproduceerd (Van Dijk et al. 1995). Dit rapport behandelt de theorie en de implementatie daarvan in het model, de aanpassingen t.o.v. de eerdere versie, een gevoeligheidsstudie en een validatiestudie.

Dit rapport vloeit voort uit een opdracht van de Dienst Weg- en Waterbouwkunde van Rijkswaterstaat voor de Technische Adviescommissie voor de Waterkeringen, werkgroep C "Zandige kusten" (overeenkomst TAWC/EOMODl).

De begeleidingsgroep voor het project "Fysische Modellering Zeereepontwikkeling" is als volgt samengesteld:

- Ir. A.P. de Looff (Dienst Weg- en Waterbouwkunde) - Dr.ir. H.J. Steetzel (Alkyon)

- Ir. H.J. Verhagen (IHE)

- Ir. J. van der K o l f f (Hoogheemraadschap van Delfland)

Voorpagina: De kleuren geven de v/indsneiheid aan in m/s als functie van de afstand tot de top van het duin en de hoogte, berekend voor een symmetrisch kaal duin van 10 m hoog. De isolijnen geven de hydrodynamisctie drukverstoring in Pa.

1 INLEIDING

Voor de berekening van de hoeveelheid zand die onder bepaalde omstandigheden getrans-porteerd kan worden is de windsnelheid dicht bij de grond of de wrijvingssnelheid een omnisbaar gegeven. In eerste benadering is de transportcapaciteit van de lucht evenredig met de derde macht van de wrijvingssnelheid (Bagnold, 1941). Ook recentere publikaties (Hsu

1971, Lettau & Lettau 1977, White 1979 en Williams 1964) voorspellen een transportcapaciteit die voor windsnelheden ruim boven de kritische windsnelheid het derde-machtsverband goed benadert.

Boven vlak en homogeen terrein kan, onder neutrale omstandigheden, de verandering van de windsnelheid met de hoogte worden beschreven met het logaritmisch windprofiel. In dat geval kan de wrijvingssnelheid worden afgeleid uit de windsnelheid op een hoogte van bijvoorbeeld twee meter en de ruwheid van het oppervlak. Het systeem dat we hier beschouwen, strand, zeereep en duin, is echter niet vlak en niet homogeen. Voor een eenvoudige ruwheidsovergang is nog wel enigszins bekend hoe de wind zich aanpast aan het nieuwe oppervlak. Zijn er meerdere ruwheidsovergangen na elkaar dan wordt het profiel te ingewikkeld voor een analytische beschrijving. Voor niet vlakke terreinen is het nog veel moeilijker een analytische oplossing te geven.

Een analytische beschrijving van het windprofiel of het verloop van de wrijvingssnelheid over strand en zeereep is onmogelijk te geven. Er tevens zijn geen eenvoudige empirische verbanden beschikbaar, die gebruikt zouden kunnen worden om een schatting te maken voor het verloop van de wrijvingssnelheid. Daarom is gekozen voor een fysisch deterministisch numeriek model, waarmee de luchtstroming over het duin, en dus ook de wrijvingssnelheid gesimuleerd wordt.

In een haalbaarheidsstudie (Van Dijk et al. 1995, Steetzel 1995) is onderzocht of een combinatie van een stromingsmodel en een model voor het aeolisch transport van zand een zinvolle bijdrage levert aan ons inzicht in de ontwikkeling van de zeereep en kan worden gebruikt als beheersinstrument. Gezien de positieve uitkomst van die studie is besloten de modellen verder te ontwikkelen. Dit rapport beschrijft de uitwerking van het stromingsmodel. Uit de haalbaarheidsstudie (Van Dijk et al. 1995) bleek bovendien dat het wenselijk was het model op een aantal punten aan te passen.

In de opdracht (TWAC/EOMODl) worden de belangrijkste punten waarop het model aangepast moet worden expliciet genoemd:

- Recirculatiewervels dienen tijdig te worden gesignaleerd, zodat het model daar niet meer op vast loopt,

- Er dienen voorzieningen te worden getroffen voor de redistributie van turbulente kinetische energie binnen één tijdstap.

- De termen die de effecten van de kromming van de stroomlijnen op de produktie en/of destructie van turbulentie dienen te worden opgenomen.

- In de configuratiefile dient een "set-variabele" te komen, waarmee de in/uitvoer wordt toegespitst op het gebruik van het SCOPE-model.

- Het verwijderen van de automatisclie slingering (a.g.v. het toegepaste spline-algoritme) in het verloop van de ruwheidslengte.

- Verder dient een validatie van het model plaats te vinden op grond van meetgegevens. De theorie in hoofdstuk 2 komt voor een deel overeen met de theorie beschreven in een eerder rapport (Van Dijk et al. 1995), maar is uitgebreid. Deze uitbreidingen zijn verwerkt in het model waarvan de implementatie in hoofdstuk 3 wordt besproken. Hoofdstuk 4 behandeh de belangrijkste wijzigingen die zijn aangebracht naar aanleiding van eerdere ervaringen (Van Dijk et al. 1995, Steetzel 1995, diverse besprekingen en de opdracht). Ten slotte komt in hoofdstuk 5 een gevoeligheidstudie en de parameterisatie aan de orde en in hoofdstuk 6 een validatie aan de hand van windmetingen over Nederlandse zeerepen.

2 T H E O R I E

Het gebruikte model is gebaseerd op het Zeman en Jensen (1987) model. Dit is een tweedimensionaal model. De X-coördinaat is in de windrichting en de Z-coördinaat verticaal. In de Y-richting (lateraal) wordt verondersteld dat er geen gradiënten zijn. Het coördinaat stelsel volgt de stroomlijnen. Er wordt voorwaarts geïntegreerd, de wind volgend (Lagrangiaanse benadering). Het model maakt gebruik van tweede orde sluiting (zie paragraaf 2.8 of StuU 1988). Het bevat dus niet alleen een vergelijking voor de gemiddelde wind maar ook voor alle relevante turbulente variabelen.

Hieronder wordt na een korte inleiding over het logaritmisch windprofiel en de aanpassing van het windprofiel bij ruwheidsovergangen een beschrijving gegeven van de belangrijkste vergelijkingen die in het model gebruikt worden. Voor de details wordt verwezen naar Zeman en Jensen (1987).

Omdat de fortran code van het oorspronkelijke Zeman en Jensen model geen commentaar bevatte, nauwelijks gedocumenteerd was en op een aantal plaatsen afweek van de in de literatuur beschreven theorie is het model volledig opnieuw geïmplementeerd. Daarbij is gekozen voor een implementatie in TurboPascal.

De hier besproken theorie is grotendeels ook al besproken in Van Dijk et al. 1995 (Modellering zandtransport zeereep. Tekenprocedures SAFE 0.0). Op een aantal punten is het model en de gebruikte theorie echter gewijzigd. Voor de duidelijkheid is besloten de theorie hier weer integraal te behandelen.

2.1 Effecten van overgangen en topografie op het windprofiel

Boven een vlak en homogeen terrein kan men onder neutrale condities de wrijvingssnelheid afleiden uit het logaritmisch windprofiel:

U{z) =

K

ln (z-d

\ ^ 0 )

[m s-'] (1)

Waarin U(z) = de windsnelheid op hoogte z, U» = de wrijvingssnelheid,

d = de nulvlaksverplaatsings t.g.v. vegetatie,

ZQ = de ruwheidslengte voor het oppervlak,

K = de Von Karman constante (0.40).

Als de atmosferische stabiliteit niet neutraal is wijkt het windprofiel af van dit logaritmisch verband. De stabiliteit wijkt meer af van neutraal als de verticale gradiënt van de temperatuur groot is en de windsnelheid laag. Aangezien verstuiving slechts optreedt bij hoge windsnelheden zullen we met deze afwijking geen rekening houden.

2.1.1 Niet homogene oppervlakken, m.n. ruwheidsovergangen

Als het terrein niet homogeen is, bijvoorbeeld door de aanwezigheid van een ruwheidsovergang, zal het windprofiel zich na de overgang gaan aanpassen aan het nieuwe oppervlak. De laag waar de invloed van de nieuwe ruwheid voelbaar is heet de interne grenslaag, de hoogte van deze laag, H, is een functie van de afstand (fetch) na de overgang, X, (Elliot 1958, Kroon 1985):

H = z,,a [xlz4-' [m] (2)

met

a = 0.75 + O.03hi(zo2/Zoi) H (3)

Hierin is ZQ, de ruwheidslengte voor de overgang en de ruwheidslengte van het nieuwe terrein. Alleen het onderste deel van deze interne grenslaag heeft weer een logaritmisch profiel. De dikte van deze aangepaste grenslaag bedraagt ongeveer 15% van de interne grenslaag (Monteith 1973). In Figuur 1 is te zien dat de dikte van de interne grenslaag bij deze overgang ongeveer 10% is van de fetch. Alleen het onderste deel van de

interne grenslaag heeft weer een logaritmisch F i g u u r l.- Hoogte van de interne grenslaag en aangepaste

windprofiel. Deze aangepaste grenslaag grenslaag (z^, = o.i mm en z „ = lo mm).

beslaat slechts ongeveer 15% van de interne

grenslaag, zodat 100 m na de overgang de dikte van de aangepaste laag slechts 1.5 m is. In het duingebied, waar vele ruwheidsovergangen voorkomen, zal dus slechts een dunne laag nabij het oppervlak een logaritmisch windprofiel vertonen. Willen we de wrijvingssnelheid bepalen uit het windprofiel, dan zijn dus windgegevens nodig van dicht bij de grond. Zeer dicht bij de grond, in de zogenaamde ruwheidslaag, (z < d + 10 ZQ) zal het gemeten windprofiel ook afwijken van het logaritmische windprofiel (Jacobs en Van Boxel 1988, Tennekes 1982). Daar zijn de effecten van de afzonderlijke ruwheidselementen voelbaar.

2.1.2 De stroming over een heuvel

Wanneer zich in de stroming zich over een duin beweegt zal zich voor het duin een overdruk opbouwen. De wind perst zich als het ware tegen het duin aan. Wordt de helling convex, dan is er een kracht nodig die moet zorgen dat de wind het duin volgt, er ontstaat een lichte onderdruk. Aan de lijzijde van het duin waar de helling weer concaaf wordt ontstaat weer een overdruk (Figuur 2).

Waar de druk in de windrichting toeneemt 2: Hydrodynamische drukverstoring en mndsnelheid

wordt de wind vertraagd en waar de druk in over een klokvormig duin zonder riiwheidsovergang.

de windrichting afneemt versnelt de wind

FLOW OVER A BELL-SHAPED DUNE

ROUGHNESS LENGTH CONSTANT

I D S P Ë E D F H E E I D S P Ë E D — \ F H E E ^ „,L,.._„^ ! / Ï / :^Nj„^L^„v,-;-..f — \ F H E E ^ „,L,.._„^ ! / Ï / " -f -•• • / y.. j

' ' i l Y D R 2DyH2DyH..P.BES.S. ..P.BES.S.

(Figuur 2). De berekening van de liydrodynamisclie drukverstoringen wordt beschreven in paragraaf 2.3.

2.1.3 Windmetingen over de zeereep en schaling

Schaling is een principe dat in de meteorologie veel wordt toegepast (Tennekes & Lumley 1972, Tennekes 1982, Stull 1988, Jacobs & van Boxel 1991, Arens et al. 1995, Jacobs et el. 1995a, 1995b). Als geschikte schalingsparameters kunnen worden gevonden krijgen de afgeleide verbanden een meer algemene geldigheid. Wanneer bijvoorbeeld in vergelijking (1) de hoogte (z-d) wordt geschaald met de ruwheidslengte (ZQ) en de windsnelheid met de wrijvingssnelheid (U»), dan krijgen alle logaritmische windprofielen dezelfde vorm.

Door Arens (1992a, 1992b, 1992c, 1994) zijn in de loop van 1990 en 1991 windmetingen gedaan op het strand en de zeereep. Deze metingen zijn verricht op drie lokaties: Schiermonnikoog, Groote Keeten en Nieuw-Haamstede. Steeds zijn v i j f tot zes masten opgesteld en is op iedere mast op drie tot vier hoogten de windsnelheid gemeten (voor meer details zie Arens 1992a, 1992b, 1992c, 1994). Na het verwijderen van evidente fouten en de metingen bij windsnelheden van minder dan 2 m/s op 5 m hoogte op het strand resteerden data over nog ruim 4000 uren. Deze zijn geanalyseerd.

Een van de belangrijkste conclusies die uit deze gegevens volgde was dat voor iedere windrichting de windsnelheid op een bepaalde hoogte geschaald kon worden met de windsnelheid op dezelfde hoogte op het strand (Arens et al. 1995). D.w.z. dat als de windsnelheid op het strand 2 maal zo groot is dat dan ook op een willekeurig punt boven de zeereep de windsnelheid 2 maal zo groot is. Voor de verschillende windrichtingen waren de resultaten wel verschillend. Van Kaam-Peters (1992) analyseerde met name de metingen gedaan op Schiermonnikoog en ontdekte dat de schaling minder goed opging voor de windmetingen aan de lijzijde van de zeereep. Dit word mogelijk veroorzaakt door het soms optreden van een recirkulatiewervel (Arens et al. 1995). Er is niet onderzocht onder welke omstandigheden deze lij-vortex optreedt. Van Kaam-Peters ontdekte wel dat de relatieve windsnelheid (windsnelheid gedeeld door de windsnelheid op dezelfde hoogte op het strand) hier een zeer zwakke negatieve correlatie vertoonde met de windsnelheid, ü p andere punten is deze correlatie afwezig.

Dat windschaling succesvol kan zijn is op een geheel verschillende lokatie bevestigd door metingen van Baas (1994) in een verstuiving in een fossiele zeereep op Ameland.

2.2 Reynoldse splitsing

Bij de beschrijving van een turbulente stroming kan men niet volstaan met een beschrijving van de gemiddelden, de vergelijkingen moeten ook uitdrukkingen bevatten voor het turbulente deel. Dit turbulente deel wordt vooral gekenmerkt door grillige fluctuaties. Op een bepaald moment en bepaalde plaats is de windsnelheid meestal niet gelijk aan de gemiddelde windsnelheid, maar hij fluctueert om dat gemiddelde. In de vergelijkingen wordt dat als volgt tot uiting gebracht:

U = U + u

W = W + w

Hierin staan de variabelen met een ~ er boven voor de momentane waarden, U , V en W voor het gemiddelde van de longitudinale (in de windrichting), laterale, respectievelijk verticale wind en u, v, w voor de turbulente afwijkingen van het gemiddelde.

De gemiddelde waarde van u-u (aangegeven door üïï) is de variantie van de longitudinale windsnelheid (Stull 1988). Evenzo stellen vv en ww de variantie van de laterale en de verticale windsnelheid voor.

Wordt de covariantie van de horizontale wind en de verticale wind (üw) vermenigvuldigd met de dichtheid van de lucht p, dan geeft dit de verticale flux van horizontale impuls. Aan de grond komt dit overeen met de wrijving. De wrijvingssnelheid is aan het oppervlak dan ook gelijk

aan V (-uw)

De uitdrukkingen uu, vv, ww en uw worden tweede orde termen genoemd. De variaties zijn uiteraard altijd positief uw is in de regel negatief; normaliter is de impulsflux naar beneden gericht omdat de windsnelheid met de hoogte toeneemt.

2.3 De hydrodynamische drukverstoring

Zeman en Jensen (1988) beschrijven de hydrodynamische drukverstoring als:

Pix,z) = -Ul^ 7 V ( ^ ) {X-S)^ ^^2,-2^ (5)

Hierin is Ug de windsnelheid op een referentiehoogte, z de hoogte, x de horizontale afstand en s een integratievariabele. r|'(s) is een vormfunctie die afhangt van de vorm van het duin. Zeman en Jensen gebruiken als referentiesnelheid de windsnelheid op een hoogte die gerelateerd is aan de horizontale lengteschaal van het duin. In hun toepassingen was dat meestal tussen de 100 en 200 m. De horizontale lengteschaal van het duin is echter niet zonder meer vast te stellen. Er zijn hiervoor verschillende definities (Arens et al. 1995), welke vooral voor onregelmatig gevormde duinen verschillende resultaten geven. Omdat bovendien de windsnelheid op hoogte groter dan 100 m niet meer sterk verandert is ervoor gekozen de windsnelheid aan de top van het model (boven de grenslaag) als referentiesnelheid te gebruiken bij de berekening van de hydrodynamische drukverstoring.

Hoe de vormfunctie ri'(s) wordt berekend, wordt in Zeman en Jensen (1987) niet vermeld. Uit de fortrancode van het oorspronkelijke model blijkt dat r|'(s) wordt berekend als de tangens van de hellingshoek ter plaatse van s gedeeld door de maximale hoogte van het duin. Dit delen door de maximale hoogte maakt dat het Zeman en Jensen model niet getest kan worden boven vlak terrein, omdat dat de drukverstoring dan oneindig groot wordt. Ook uit dimensie-overwegingen is dit delen dor de hoogte van het duin onwenselijk, daar dan de

eenheden van de berekende drukverstoring niet meer kloppen. Daarnaast lijkt de keuze van de tangens van de hoek geen gelukkige, omdat de tangens bij steile hellingen zeer groot wordt. In onze implementatie van het model is daarom gekozen voor een vormfactor die evenredig is met de sinus van de hellingshoek. De sinus is immers voor hellingshoeken kleiner dan 30° vrijwel gelijk aan de tangens, maar is voor een hoek van 90° gelijk aan 1.0 i.p.v. oneindig groot. Als evenredigheidsfactor is niet de hoogte van het duin gekozen, maar een dimensieloze parameter, die onafhankelijk is van het duin. Voor de waarde van deze parameter is voorlopig 0.20 gekozen. De juiste waarde kan pas vastgesteld worden als het model uitgebreid getest is. De vormfunctie luidt nu dus:

^'(s) = Par sm(a)

Waarin: a = hellingshoek ter plaatse van s. Par = Evenredigheidsparameter (0.20).

Op de voorpagina is de hydrodynamische drukverstoring weergegeven in Pa als functie van de afstand tot de top van het duin en de hoogte, die wordt berekend bij een windsnelheid van 9 m/s (op 10 m hoogte) voor een symmetrisch kaal duin van 5 m hoog.

2.4 Vergelijking voor de horizontale wind

l n prognostische meteorologische modellen wordt voor de berekening van de windsnelheid en de verandering daarvan vaak gebruik gemaakt van de Navier-Stokes vergelijking:

^ = - 1 ^ + fV - ^ [m s-^] (7) dt p dx dz

Hierin is p de dichtheid, dp/dx de grootschalige drukgradiënt in de x-richting, f de Coriolisparameter en düw/dz de verticale gradiënt van de impulsflux (wrijvingskracht). De eerste term geeft het effect van de horizontale drukgradiënt op de wind, de tweede het effect van de Corioliskracht en de derde term het effect van de wrijvingskracht.

Het gebruikte model is een Lagragiaans model, waarin de lucht wordt gevolgd over een afstand dx. Daarom is dU/dt vervangen door U-dU/dx. Door het gebruik van stroomlijncoördinaten en de kromming van de stroomlijnen worden extra termen geïntroduceerd (Zeman en Jensen 1987):

IJ dU ^ _ dP ^ . d^ _ d^ _ W d^ [in (8) dx dx ^ ^' dz dx u dz

I I I I I I I V V

De eerste term beschrijft de reactie van de wind op de hydrodynamische drukverstoring, de horizontale drukgradiënt die ontstaat door de aanwezigheid van het duin.

De tweede term bevat de grootschalige drukgradiënt (9p/9x = f V g , met Vg = de component van de geostrofische wind in de laterale richting). Deze zal grotendeels gecompenseerd worden door de Corioliskracht. Echter vooral dicht bij de grond zal de wind niet geostrofisch zijn en is een correctie noodzakelijk.

De derde term bestaat uit de verticale divergentie van uw. puw staat voor de wrijving. Als de wrijvingskracht met de hoogte verandert zal dat gevolgen hebben voor de windsnelheid. De voorlaatste term komt niet voor in modellen die ontwikkeld zijn voor homogene oppervlakken, daar in dat geval in de horizontale richting geen gradiënten voorkomen. De laatste term is het gevolg van het feit dat de stroming zich niet horizontaal over de zeereep kan bewegen en dat in stroomlijncoördinaten wordt gerekend. Voor de transformatie van een orthogonaal stelsel naar een stroomlij nstelsel wordt de Von Mises transformatie toegepast (Milne-Thompson 1949). Voor details omtrent de transformatieregels wordt verwezen naar Zeman en Jensen (1987).

2.5 Tweede orde vergelijkingen voor de varianties

l n het model komen een vijftal tweede orde vergelijkingen voor, namelijk drie voor de respectievelijke varianties, één voor de impulsflux (uw) en één voor de dissipatie.

De varianties worden berekend met de volgende vergelijkingen:

U U U duu dw dww = - luw du 2uu dU dz dx prod w 2uw U R I-T C ' uu- w-.— du w-.— u c + 2ww— + 4MW— - — dx R T I- ww-- RD. RD. RD. duuw dz 2 —e 3 I I I I I I IV dww dwww dz V I 3 V I I (9) (10) (11)

Hierin is R de kromtestraal van de stroming (negatief voor convexe stroomlijnen), T de dissipatie tijdschaal (T = qq/s), C een evenredigheidsconstante (C = 3.25; Zeman & Jensen 1987), qq de turbulente kinetische energie (deze is gedefinieerd als qq s üu+vv+ww) en s de snelheid waarmee deze turbulente kinetische energie wordt gedissipeerd door wrijving aan de rand van de kleinste turbulente wervels. Voor de dissipatie bevat het model een aparte vergelijking (zie paragraaf 2.7).

In de fortrancode van het Zeman en Jensen model wordt de kromtestraal R benaderd als de reciproke van de tweede afgeleide van de hoogte van de stroomlijn naar x. Hoewel dit formeel niet juist is, vormt het een goede benadering, zolang de helling van de stroomlijn niet steiler is dan 20°. In ons model is gebruikt gemaakt van een nauwkeurigere benadering (Zeman & Jensen 1987):

fdz^^"^'"

R 1 +

dx) _dx'

[m] (12)

Waarin z de hoogte van de stroomlijn is ten opzichte van een vast niveau (bijv. t.o.v. NAP). De eerste drie termen in (9), (10) en (11) zijn produktietermen. Zij geven aan welke factoren van belang zijn voor de produktie van turbulente fluctuaties in de drie richtingen. De eerste

term is in ieder tweede orde stromingsmodel aanwezig, de tweede term komt alleen voor bij inhomogene oppervlakken en de derde alleen bij gekromde oppervlakken. Bij de beoordeling van de verschillende produktietermen termen moet men vooral letten op het teken van de tweede orde termen. De varianties zijn uiteraard ahijd positief, üw is in de regel negatief De produktie van turbulentie vindt voornamelijk plaats in de longitudinale richting. Waar de wind vertraagt (9U/9x negatief) vindt extra produktie plaats in de longitudinale richting, terwijl de bijdrage aan de verticale variantie juist negatief is. Versnellingen hebben een omgekeerd effect. Ook concave hellingen (R is positief) geven een produktie van longitudinale variantie en een negatieve produktie van verticale variantie, terwijl dit op convexe hellingen weer juist omgekeerd is.

De vierde term is de redistributieterm. Omdat de produktie van turbulente kinetische energie voornamelijk plaatsvindt in de longitudinale richting zal de longitudinale variantie groter zijn dan de laterale en de verticale. De snelheidsfluctuaties geven echter aanleiding tot drukfluctuaties. Deze drukfluctuaties produceren weer turbulente snelheidsfluctuaties. Daar druk echter alzijdig werkt zullen deze fluctuaties gelijkelijk over alle richtingen verdeeld zijn. Dit wordt de redistributie van turbulente kinetische energie genoemd. T/C is een tijdschaal die aangeeft hoe snel het redistributieproces verloopt. In navolging van Zeman en Jensen (1987) is hier gekozen voor een lokale schaling. In een eerdere versie van het model werd gebruik gemaakt van een mengwegschaling (Wilson & Shaw 1977), waarbij de mengweglengte recht evenredig met de hoogte boven het oppervlak werd genomen.

De vijfde term is de rapid distortion, een herverdeling van turbulente kinetische energie over de drie richtingen doordat de wervels vervormd worden in een stroming waar een snelheidsgradiënt bestaat. De parameterisering van de rapid distortion wordt in de paragraaf 2.5.1, 2.5.2 uitgewerkt.

De zesde term is een transport term. Net zoals uw het verticale transport van horizontale impuls door turbulente wervels voorstelt, kan iluw opgevat worden als het verticale turbulente transport van longitudinale fluctuaties. Komt er aan de ene kant meer binnen dan er aan de andere kant uitgaat, dan levert dat een bijdrage aan de toename van de longitudinale variantie. Hetzelfde geldt voor 9vvw/9z en dwww/dz.

De laatste term in de vergelijkingen geeft de dissipatie, de snelheid waarmee deze turbulente kinetische energie wordt gedissipeerd door wrijving aan de rand van de kleinste turbulente wervels. Omdat verondersteld wordt dat de turbulentie op de schaal van de kleinste wervels isotroop is wordt de dissipatie gelijkelijk over de drie varianties verdeeld.

De turbulente kinetische energie qq is hier gedefinieerd als qq = uu+vv+ww (soms wordt ook gedefinieerd TKE = 0.5qq). Als de vergelijkingen voor de varianties worden opgeteld ontstaat de vergelijking voor de turbulente kinetisch energie qq:

U = - luw — - l(uu-ww) — + luw— - - le K^J)

dx dz dx R dz

I I I I I I V I V I I

De redistributieterm (IV) en de rapid distortion (V) komen in deze vergelijking niet voor omdat deze alleen zorgen voor een herverdeling van de TKE over de verschillende componenten, maar geen verandering brengen in de totale hoeveelheid TKE.

Het effect van de kromtestraal op de totale turbulente kinetische energie is niet evident. Term 111 suggereert dat op convexe hellingen (duinvoet) de TKE afneemt, omdat uw steeds negatief is en R > 0. De kromtestraal heeft echter ook een effect op de waarde van uw (zie ook par. 2.6). Door de krommingseffecten neemt op een convexe helling de absolute waarde van uw toe. Term 1 zorgt dus tegelijkertijd voor een extra produktie van TKE. Tot dusverre is niet onderzocht welk effect overheerst. Op concave hellingen (bijv. duintop) zijn de effecten uiteraard omgekeerd.

2.5.1 Rapid distortion

Snelheidsfluctuaties geven aanleiding tot drukfluctuaties, welke als gevolg hebben dat de turbulente kinetische energie, die voornamelijk wordt geproduceerd in de longitudinale en in iets mindere mate in de verticale richting, wordt herverdeeld over alle richtingen. In een stroming waar snelheidsgradiënten voorkomen zullen de turbulente wervels vervormd worden. Ook dat gaat gepaard met drukfluctuaties die voor een herverdeling van de turbulente kinetische energie zorgen. Dit laatste proces wordt rapid distortion genoemd. Rapid omdat het een proces is dat erg snel werkt. Binnen zeer korte tijd wordt een deel van de produktie herverdeeld.

Hieronder wordt de rapid distortion steeds gesplitst in de rapid distortion van de normale produktie en de rapid distortion van de krommingsproduktie:

RD^^ = RD!^ + RD' (14)

De parameterisatie van de rapid distortion die beschreven is in Zeman en Jensen (1987) is gebaseerd op de rapid distortion theorie van Zeman en Tennekes (1975).

2.5.2 Rapid distortion in variantievergelijkingen

In de parameterisatie van de rapid distortion van de normale produktie termen zijn vaak de produktie termen herkenbaar:

RD. RD. (2 2 a, — dU uw — dz 4 — ac/ + — a. uw — 3 ^ dz RD p _ a. - 2a^ — dU uw — dz 4 / .dU — a. [uu-w) — 3 ' dx 4 / ^ dU — aAww-uu) — 3 dx 4 , .dU — a, ( w - w w — 3 ' dx 2 dU 5 dx 2 dU 5 dx (15) (16) (17)

De rapid distortion parameters a, en zijn beide positief en van dezelfde orde van grootte (circa 0.3). Boven een homogeen terrein (dU/dx = 0) leidt een grote verticale windsnelheidsgradiënt (dU/dx) en/of een grote wrijving (-üw) tot een rapid distortion die groot en positief is voor de longitudinale variantie ü ü en negatief voor de laterale variantie vv en de verticale variantie ww. Aangezien de rapid distortion termen in de variantievergelijkingen voorkomt met een minteken betekent dit dat TKE, die voornamelijk wordt geproduceerd in

de longitudinale richting via rapid distortion wordt omgezet in snelheidsfluctuaties in de laterale en verticale richting.

Voor de parameters a, en a2 geven Zeman en Jensen (1987) aan dat zij op theoretische gronden beiden gelijk zouden moeten zijn aan 0.30. Op een andere plaats in hun publikatie geven ze echter aan dat ze gewerkt hebben met a,=0.375 en a2=0.225. Na een aantal substituties en het optellen van termen in hun model bleken deze waarden ook verwerkt te zijn in de fortrancode van hun model. Ook in H I L L _ M D L wordt standaard gerekend met a,=0.375 en a2=0.225. Als men w i l rekenen met andere waarden dan kan men in H1LL_PAR.PAR waarden opnemen voor de parameters ALFA_0 en DEL_ALFA, waarbij a,=ALFA_0+DEL_ALFA en a2=ALFA_0+DEL_ALFA.

Omdat rapid distortion zorgt voor een herverdeling van TKE over de drie richtingen geldt uiteraard dat RD^^ + RD^^ + RD^„^ = O en komt de rapid distortion term niet voor in de TKE-vergelijking (13).

2.5.3 Rapid distortion in impulsfluxvergelijking

Ook voor de impulsfluxvergelijking geven Zeman en Jensen (1987) een rapid distortion term:

Waarin b,, = uu-qV3 en b33 = ww-qV3.

De bijdrage van de eerste term (-1/5 qq 9U/9z) aan de rapid distortion (gevolg van het vervormen van wervels) is negatief, aangezien qq positief is en 9U/5z over het algemeen ook. De rapid distortion staat met een minteken in de impulsfluxvergelijking, zodat dit leidt tot een "grotere" uw. Aangezien üw echter negatief is zal de absolute waarde afnemen. Deze term in de rapid distortion leidt dus tot een afname van de wrijving. Dat wordt gedeeltelijk teniet gedaan door de laatste term omdat b33 over het algemeen negatief is. Als a, en gelijk zijn aan hun theoretische waarde van 0.3 zal de tweede term, die b,, bevat, wegvallen.

In het HILL-model is voor de rapid distortion in de impulsfluxvergelijking gebruik gemaakt van een uitdrukking die is ontleend aan Wilson en Shaw (1977):

= - C , q q ^ (19)

Wilson & Shaw (1977) gebruikten voor CR een waarde van 0.07. Op grond van de rapid distortion theorie van Zeman en Jensen verwachten we een waarde die iets onder de 0.20 ligt. Deze waarden lopen nogal uiteen. Een onjuiste keuze van deze parameter maakt dat het model boven een homogeen terrein geen stationaire toestand meer zal berekenen.

Van de eigenschap dat een onjuiste keuze van de parameter C^ boven homogeen terrein geen stationair windprofiel geeft is gebruikt om een goede schatting te verkrijgen van de waarde van de parameter C. In de initialisatiefase (een stroming over een homogeen terrein) past het model deze parameter zodanig aan dat de stroming stationair wordt. Opmerkelijk is dat steeds een waarde wordt gevonden van ongeveer 0.14, onafhankelijk van de gekozen ruwheid of

windsnelheid. De waarde van 0.14 ligt dichtbij de waarde die volgt uit de rapid distortion theorie van Zeman en Jensen (1987). Overigens zal de exacte waarde nog beïnvloed worden door de keuze van de andere parameters in het model, zoals de initiële anisotropie in de TKE-verdeling.

Misschien zou het aan te bevelen zijn om ook voor de impulsfluxvergelijking de rapid distortion formulering van Zeman en Jensen (1987) te gebruiken en dan de parameter DEL_ALFA aan te passen om het windprofiel boven homogeen terrein stationair te krijgen.

2.5.4 Rapid distortion kromniingsproduktie

Ook de produktie (of destructie) van turbulente kinetische energie die plaatsvindt als gevolg van de kromming van de stroomlijnen wordt via rapid distortion snel herverdeeld over de drie richtingen. Zeman en Jensen (1987) geven daarvoor de volgende uitdrukkingen:

RD',^, = -(16/3) ^ U/R (20)

RD^ = - (4/3) « 2 U/R (21)

RD^l = +(20/3) ^ U/R (22)

i?D;,, = 2 a^im - ^ / 2 ) U/R (23)

Ook deze uitdrukkingen zijn verwerkt in het model. Bij de implementatie van deze theorie is om redenen van duidelijkheid gepoogd de source code van het model zo goed mogelijk te laten aansluiten bij de theorie. Als illustratie hiervan worden in bijlage 5 de source code getoond van de procedures die gebruikt worden voor het berekenen van de termen van de difi'erentiaalvergelijkingen voor ïïu, vv en ww.

2.6 De vergelijldng voor de impulsflux

De impulsflux (wrijving) wordt beschreven door de vergelijking:

JJ duw dU — \ U C — duww

U = - WW— + (2MM - ww) — — uw RD„^,

-dx dz ^ ' R T ""^ dz

I I I I I V V V I

De betekenis van de termen is grotendeels vergelijkbaar met de vergelijkingen voor de varianties. De eerste produktieterm is weer in ieder tweede orde sluitingsmodel aanwezig. De tweede en derde term hangen samen met de convexiteit van de helling. Op concave hellingen wordt extra wrijving geproduceerd, terwijl op convexe hellingen de wrijving juist vermindert. Dit draagt bij aan de hoge snelheidsgradiënten die kunnen ontstaan op de top van de zeereep. Dit laatste effect lijkt op wat er gebeurd in een stabiel gestratificeerde grenslaag, ook daar is de wrijving gering maar zijn de windsnelheidsgradiënten groot.

De vijfde term is de rapid distortion. Voor de rapid distortion van de normale produktietermen is gebruik gemaakt van een formulering die ontleend is aan Wilson & Shaw (1977). Deze formulering wijkt enigszins af van de formulering van Zeman en Jensen. De waarde van de parameter CR wordt door het model zelf bepaald tijdens de initialisatie fase (zie paragraaf 2.5.3) en bedraagt steeds circa 0.14.

De vijfde term is weer een transportterm.

2.7 De dissipatievergelijking

ln sommige tweede orde modellen komt geen expliciete vergelijking voor de dissipatie voor, maar worden bepaalde aannamen gemaakt voor het verband tussen de dissipatie en sommige andere variabelen in het model (bijv. de turbulente kinetische energie) (Stull 1988). In het Zeman en Jensen (1987) model wordt deze methode van sluiting niet toegepast, dit model bevat wel een vergelijking voor de dissipatie:

Hierin is IJ de totale produktie aan turbulente kinetische energie (de som van de eerste drie termen in de variantievergelijkingen), p een dimensieloze model parameter die volgens Zeman en Jensen (1987) de waarde 0.75 heeft, en de turbulente uitwisselingscoëfficiënt voor dissipatie. Ook het getal 3.8 is in feite een modelparameter.

De eerste term is een soort produktieterm. Is de dissipatie veel groter dan de produktie van turbulente kinetische energie, dan zal deze term de dissipatie doen afnemen. Is de dissipatie daarentegen relatief klein dan beschrijft deze term een toename van de dissipatie. De andere term is een transportterm, die vergelijkbaar is met de triple produkten in de andere tweede orde vergelijkingen. De turbulente uitwisselingscoëfficiënt, K^,, die hier gebruikt wordt is echter volgens Zeman en Jensen (1987) 0.68 maal zo groot als de turbulente uitwisselingscoëfficiënt voor de andere termen.

Het fysisch fundament onder de dissipatievergelijking is niet erg sterk, zeker niet in vergelijking tot de andere tweede orde vergelijkingen. In feite komt het er op neer dat de hier gebruikte vergelijking niet veel beter onderbouwd is dan de sluiting die in vele andere tweede orde modellen wordt toegepast voor de dissipatie.

2.8 Aannamen voor sluiting van de 2''^-orde vergelijkingen

In het voorgaande bleek dat in de vergelijkingen voor eerste orde variabelen (i.e. U) tweede orde termen verschijnen. Het is mogelijk voor de tweede orde variabelen weer vergelijkingen af te leiden. Er zijn echter meer tweede orde variabelen dan eerste orde variabelen. Bovendien zijn de tweede orde vergelijkingen ingewikkelder en bevatten ze termen die moeilijk te interpreteren zijn. Een nog groter probleem is dat in de tweede orde vergelijkingen weer derde orde termen verschijnen. Als voor deze termen weer vergelijkingen worden afgeleid zijn dat er nog meer, ze zijn nog ingewikkelder, nog moeilijker te interpreteren en ze bevatten vierde orde termen. Dit probleem is theoretisch onoplosbaar (Stull 1988). Men houdt ahijd meer onbekenden dan vergelijkingen. Op enig niveau zullen dus aannamen gemaakt moeten worden

voor de hogere orde termen. Dit heet het sluiten van het stelsel van vergelijkingen. In het Zeman en Jensen model gebeurt dat op het tweede orde niveau. Vandaar dat gesproken wordt van een tweede orde sluitingsmodel.

In de vergelijkingen zoals genoemd in de vorige paragrafen zitten al enkele aannamen voor sluiting verwerkt. Deze hebben voornamelijk betrekking op drukfluctuaties: de redistributieterm in de vergelijkingen voor de varianties en term IV in de impulsflux-vergelijking (vergelijkbaar met redistributieterm) en term V (rapid distortion).

De triple produkten worden geïnterpreteerd als transporttermen. Net zoals üw het verticale transport van horizontale impuls door turbulente wervels voorstelt, kan üüw, vvw, respectievelijk www, opgevat worden als het verticale turbulente transport van longitudinale, laterale respectievelijk verticale fluctuaties. Het is in de natuurkunde heel gebruikelijk de flux te berekenen als het produkt van een gradiënt en een uitwisselingscoëfficiënt.

duu -—• , , 3 w uuw = - K , WW = - K dz dz (26) dww r^duw www = - K , uww = - K dz dz

ln het geval van laminaire stromingen is deze uitwisselingscoëfficiënt, K, de diffusieconstante. Bij een turbulente stroming is de uitwisselingscoëfficiënt afhankelijk van de hoogte. Vaak wordt een benadering toegepast waarbij de turbulente uitwisselingscoëfficiënt wordt berekend uit de mengweglengte. Hier is er echter gekozen voor de bepaling van de turbulente uitwisselingscoëfficiënt op grond van lokale parameters (Zeman & Jensen 1987):

K = 0.075 '^qq l e [m's"'] (27)

Wanneer deze turbulente uitwisselingscoëfficiënt wordt vergeleken met de turbulente uitwisselingscoëfficiënt die wordt berekend met de mengwegtheorie (Wilson & Shaw 1977), dan blijken de verschillen niet erg groot te zijn. Bovendien is meestal de transportterm ook niet de belangrijkste term in de variantie- en impulsfluxvergelijkingen.

2.9 Randvoorwaarden

Het model heeft voor alle variabelen randvoorwaarden nodig voor de onderzijde (bodem of vegetatie), en bovenzijde.

De bovenzijde van het model is gekozen ruim boven de atmosferische grenslaag. Hier zal de windsnelheid geostrofisch zijn. Deze blijft gedurende de simulatie constant. Omdat boven de grenslaag de turbulente fluctuaties verwaarloosbaar zijn ten opzichte van de fluctuaties in de grenslaag wordt voor de varianties en de impulsflux aangenomen dat deze aan de bovenkant van het model gelijk zijn aan 0.0. Als er geen turbulentie is aan de bovenkant zal uiteraard ook de dissipatie nul zijn.

De onderrand van het stromingsmodel ligt op de topografische hoogte plus de verplaatsingshoogte (d) plus de ruwheidslengte (zg). Als onder-randvoorwaarde voor de varianties en de impulsflux is aangenomen dat de gradiënt tussen het onderste gridpunt en de

bodem nul is. Voor de windsnelheid wordt aangenomen dat deze nul is aan de grond. De windsnelheid op het eerste gridpunt wordt bepaald uit het logaritmisch windprofiel, de lokale ruwheidslengte en de wrijvingssnelheid aan de grond (U» = ^ ( - u w ) ) . Voor uu, vv en ww is aangenomen dat ze recht evenredig zijn met -uw. De evenredigheidsconstanten (4.4, 2.5 resp. 1.6) werden ook gebruikt bij het berekenen van de initiële profielen van uu, vv en ww. Het startprofiel voor de gemiddelde wind is nabij het oppervlak een logaritmisch windprofiel, hoger in de grenslaag is de verandering van de windsnelheid met de hoogte minder dan het logaritmisch profiel voorspelt. Boven de grenslaag is de windsnelheid constant.

Voor de beginprofielen van de overige variabelen worden betrekkelijk willekeurige aannamen gemaakt. In de initialisatiefase van het model wordt de stroming echter over een uniform vlak terrein geleid, dat dezelfde eigenschappen heeft als het terrein aan de bovenwindse kant van de te modelleren zeereep. Hierdoor krijgen de profielen de gelegenheid zich aan te passen aan het terrein.

3 MODELBESCHRIJVING

3.1 structuur van het model

De source code van het model is verdeeld over vier Pascal files: - PIILLTYPE.PAS Type definities en globale variabelen, - HILLSOUN.PAS Geluidseffecten,

- HILLCALC.PAS Enkele algemene rekenprocedures, - H1LL_MDL.PAS De rest (het grootste deel dus).

3.1.1 Het hoofdprogramma

Bijlage 1 bevat een volledige listing van het hoofdprogramma. De belangrijkste punten zijn samengevat in Tabel 1.

Eerst wordt met de procedures Read_Param en Read_Config een aantal instellingen ingelezen uit de files H1LL_PAR.PAR en HILL_CFG.CFG. De parameterfile H1LL_PAR.PAR bevat parameters welke in principe niet door de gebruikers gewijzigd moeten worden. Is de parameterfile niet beschikbaar dan wordt gerekend met de standaardwaarden van de parameters. In bijlage 2

w o r d t beschreven w e l k e parameters het windmodel.

H1LL_PAR.PAR kan bevatten en wat het format van de parameterfile moet zijn.

De configuratiefile bevat parameters welke wel standaard door de gebruiker gewijzigd worden, zoals bijvoorbeeld de wrijvingssnelheid of de windsnelheid op 10 m hoogte aan de bovenwindse zijde van het model. Bijlage 3 wordt het format van de configuratiefile en de betekenis van de op te nemen parameters beschreven.

Na het inlezen van de parameterfile en de configuratiefile worden de gegevens omtrent de topografie (Read_Dune) en de ruwheidslengte (Read_Roughness) ingelezen en worden deze op vaste gridafstanden weggeschreven (Write_Dune).

Vervolgens worden eerst vrij willekeurige startprofielen berekend voor alle variabelen (lnitial_Profiles).

In de file H I L L MDL.CFG (zie bijlage 3) staat ook of de initialisatiefase doorlopen moet worden. Deze initialisatiefase bestaat uit twee procedures. In de procedure Time_Steps loopt het model over een vlak en uniform terrein met vaste tijdstappen. Tijdens deze eerste initialisatiefase ligt het windprofiel nog vast en worden de profielen van de turbulente

PROCEDURE DO_MAIN; BEGIN Read Param; Read Config; Read Dune; Read Roughness; Write Dune; I n i t i a l P r o f i l e s ; I f ( I n i t i a l i s e = TRUE) Then Begin Time Steps; Flat Surface; Generate Output ( ' I n i t i a l i n ' ) ; End Else R e a d _ V a r i a b l e s ( ' I n i t i a l i n ' ) ; H i l l j l o w ; END;

grootheden aangepast. Daarna gebeurt vrijwel hetzelfde maar dan met vaste afstandstappen (Nu verandert ook het windprofiel). De tijdstap neemt in dit laatste geval met de hoogte af omdat de windsnelheid met de hoogte toeneemt. Tijdens de initialisatiefase wordt ook de rapid distortion parameter, C^, uit de impulsfluxvergelijking voortdurend aangepast (zie paragraaf 2.5.3). Als de initialisatiefase doorlopen wordt dan worden de resultaten hiervan weggeschreven in de file I N I T I A L . I N (Generate_Output).

Wordt de initialisatiefase niet doorlopen dan worden dat de resultaten hiervan gelezen uit de file INITIAL.IN (Read_Variables).

De echte simulatie over de zeereep, waarvan de topografie en de ruwheid met de afstand verandert, gebeurt vervolgens met de procedure Hill_Flow.

3.1.2 De procedures TimeJSteps, Flat_Surface en Hill_Flow

De procedures Time_Steps, Flat_Surface en Hill_Flow hebben vrijwel dezelfde opbouw. In bijlage 4 is een listing van de procedure Hill_FIow opgenomen. Tabel 2 geeft een samenvatting van de procedures Time_Steps, Flat_Surface en Hill_Flow.

Nadat in het begin van de procedure een aantal variabelen een waarde hebben gekregen begint de while loop (in Time_Steps "While (t < tend) do"). Zolang de eindtijd of de eindafstand nog niet bereikt is wordt deze lus doorlopen.

Eerst wordt eventueel uitvoer weggeschreven naar file. PROCEDURE H111 Flow; BEGIN While (X < Xend) do Begin Maybe_Output Sound Signal Calc_Pressure Calcuiate_ü_over_R Prepare Calculations Hi gher_Örder_Equati ons Calculate New_Values

## Correct_da't_uw; ## niet in Hi 11 flow Max_Ti nie_or_Di stance

Change_dt_or_dx Calculate_New Values Smoothing;

Boundary_Conditions Calculate_WindSpeed; Copy New t o Old End; " " . END

I n H i l l _ F l o w worden vervolgens de

hydrodynamische drukverstoring en de gradiënt Xabel 2 ; SamenvaUingvan de procedures Tunejteps.

hiervan berekend (Calculate_Pressure) en de Flat_Surface en HUI Flow. §#..## wil zeggen niet

k r o m m i n g v a n de s t r o o m l i j n e n HUIFIOW.

(Calculate_U_over_R). I n initialisatiefase

(procedures Time_Steps en FIat_Surface) wordt dit deel overgeslagen, omdat daar wordt gerekend met vlak en homogeen terrein.

Dan worden alle termen van de tweede orde vergelijkingen berekend (Higher_Order_Equations).

Tijdens de initialisatie procedures wordt zonodig de rapid distortion parameter CR in vergelijking (19) en (24) aangepast. Als dat nodig was wordt ook de verandering van üw gecorrigeerd voor de in de berekeningen gebruikt onjuiste waarde van de parameter CR. Omdat vooral Time_Steps (maar in mindere mate ook Flat_Surface) wel eens geconfronteerd kan worden met onrealistische beginprofielen wordt gekeken over welke tijd of afstand maximaal vooruit geïntegreerd kan worden (Max_Time_or_Distance) en wordt zonodig de stapgrootte aangepast (Change_dt_or_dx). Ook in Hill_Flow wordt deze procedure

aangeroepen, omdat ook ruwheidsovergangen en/of de topografie kunnen zorgen voor snelle veranderingen in de stroming, die het noodzakelijk maken de stapgrootte te reduceren. Aan de hand van de berekende termen van de tweede orde vergelijkingen worden met de al dan niet aangepaste stapgrootte de nieuwe waarden van de tweede orde variabelen berekend (Calculate_New_Values). Voor de randvoorwaarden wordt de procedure Boundary_Conditions aangeroepen.

Pas daarna wordt de windsnelheid berekend (Calculate_WindSpeed). In Time_Steps wordt nog geen nieuwe windsnelheid berekend, omdat eerst het turbulente deel (de hogere orde termen) zich goed moeten instellen. De procedure Copy_New_To_01d bewaart de berekende resultaten in aparte arrays, welke in de volgende loop weer gebruikt kunnen worden.

Tinie_Steps/F1atSurface/Hin_Flow -> Maybe Output -> Generate Output -> Sound Signal -> Calc Pressure -> Cubic Eval -> Calculate U over R -> Prepare Calculations -> Calculate ddz

-> Calculate Eddy -> _{Set_Fi rst_and_Last}

-> Calculate Time Scales

-> Higher Order Equations

-> Calculate Terms uu -> Calculate Transport -> Calculate Terms" "vv -> Calculate Transport -> Calculate Terms" "ww -> Calculate Transport -> Calculate Terms' _"qq -> Calculate Transport -> Calculate Terms" uw -> Calculate Transport -> Calculate Terms" eps -> Calculate Transport -> Calculate ddx

-> _{Correct ddt uw} -> Max Time or Distance

-> Change dt or dx

-> Calculate New Turbulence -> Calculate New

-> _{V e r t i c a l Smoothing} -> Boundary_Condi t i ons

-> Calculate Terms Wind

-> _{Calculate~ddx (Rind)} -> Calculate New (Wind)

-> _{Copy New t o Old}

T a b e l 3 ; Schema van de belangrijkste procedure aanroepen vanaf het niveau van Time_Steps, Flat Surface of Hill_Flow.

3.1.3 De overige procedures

Het zou te ver voeren om in het kader van dit rapport voor iedere procedure een gedetailleerde beschrijving te geven. Een schema van de belangrijkste procedure-aanroepen vanaf de procedures Time_Steps, Flat_Surface en Hill_Flow is gegeven in Tabel 3.

In de procedures Calculate_Terms_uu/vv/ww/qq/uw/eps is gepoogd om de Pascal code zo goed mogelijk aan te laten sluiten bij de theorie die in het voorgaande besproken is (dit in tegenstelling tot de fortran code van Zeman en Jensen). Als voorbeeld zijn enkele van deze procedures opgenomen in bijlage 5.

3.2 Invoer en uitvoer files

Het model gebruikte de volgende invoerfiles: - H1LL_PAR.PAR Parameterfile (optioneel), - H1LL_CFG.CFG Configuratiefile (verplicht),

- PROFDUNE.IN Topografische hoogte (verplicht), - PROF_Z0.1N Ruwheidslengte (verplicht).

De opbouw van de parameterfile HILL PAR.PAR en de configuratiefile H I L L CFG.CFG is toegelicht in de bijlagen 2 en 3.

De files PROFDUNE.IN en PROF ZO.IN bevatten op iedere regel een X-coördinaat en een hoogte respectievelijk ruwheidslengte. De afstand dient op te lopen, er wordt niet gesorteerd. De files worden afgesloten met een regel "9999 9999".

De file INITIAL.IN wordt geschreven als het model wordt gerund met de parameter INITIALIZE=TRUE. Is deze parameter FALSE, dan wordt I N I T I A L . I N gelezen. In deze file staan de initiële profielen van U , uu, vv, ww, qq, uw en de dissipatie. Deze profielen zijn geschaald met de juiste macht van U* (U* voor U , U»^ voor uu, vv, ww, qq en uw, U*^ voor de dissipatie). Door de schaling zijn alle grootheden (behalve de hoogte en de dissipatie) dimensieloos. Door de waarden in INITIAL.IN te schalen hoeft de initiële fase niet opnieuw te worden doorlopen als het model wordt gerund met een andere waarde van U , of 10 m windsnelheid. Het initialiseren moet wel overgedaan worden als de bovenwindse ruwheidslengte veranderd of de verticale verdeling van de gridpunten.

De uitvoerfile SCOPEJN.PRN bevat de informatie die nodig is voor het SCOPE model. Iedere regel in deze file bevat de afstand, gevolgd door het quotiënt van de wrijvingssnelheid en de bovenwindse wrijvingssnelheid (Z, U./U*o). Deze uitvoer wordt iedere 1.0 m geproduceerd. Een voorbeeld van de uitvoerfile is opgenomen in bijlage 7.

Ais de parameter SCOPE op FALSE wordt gezet (SCOPE=FALSE opnemen in HILL_CFG.CFG) dan worden meer uitvoerfiles geproduceerd (DUNESLOP.PRN, PRESSURE.PRN, PROFILE.PRN en HILL_OUT.PRN). De inhoud en het format van deze uitvoerfiles wordt in dit rapport niet besproken. Normaal zal het model worden gerund met de parameter SCOPE=TRUE. SCOPEJN.PRN is de enige informatie die het SCOPE model nodig heeft. De extra uitvoer in deze files kan gebruikt worden om de resultaten meer gedetailleerd te evalueren. Op deze manier zijn de gegevens verkregen voor de produktie van de grafieken in de hoofdstukken 5 en 6.

4 AANPASSINGEN IN H E T MODEL

4.1 Recirculatiewervels

Omdat het HILL-model een voorwaarts integrerend model is vormen eventuele recirculatiewervels een probleem. In een recirculatiewervel waait de wind namelijk de andere kant op.

Het HILL-model begint aan één kant te rekenen (normaal de zeewaartse zijde) en van daar uit wordt stapsgewijs steeds uit de condities op één punt de stroming op het volgende punt benedenstrooms berekend. In het geval de windrichting omdraait zou het punt van waaruit gerekend wordt aan de landwaartse zijde moeten liggen. Aangezien de berekening is gestart vanaf de andere kant zijn de condities in het uitgangspunt nog niet bekend en kan dus ook de stroming in het nieuwe punt niet berekend worden.

Voor de toepassing van het model op de zeereepontwikkeling is een exacte modellering van de stroming in de recirkulatiewervel niet erg belangrijk. In die wervel zijn de windsnelheden sowieso zo laag dat het transport stopt. De oude versie van het HILL-model crashte echter als er recirculatiewervels optraden.

In het nieuwe model worden recirculatiewervels vermeden door voor alle grootheden slechts een beperkte variatie toe te staan rond de waarde die berekend is in het initiële profiel (boven het strand). Voor de windsnelheid is de toegestane variatie een factor 4, voor de turbulente grootheden wordt een variatie van een factor 16 toegestaan.

Aangezien de windsnelheid slechts een factor 4 kan variëren rondom het initiële profiel kunnen negatieve waarden niet meer voorkomen en is in het model recirculatie uitgesloten. Aangezien de impulsflux (uw) van de belangrijke variabelen de enige is die negatief is, wordt in het Var Record voor uw in de arrays uw.New en uw. Old niet de waarde van uw opgeslagen, maar -1.0*uw. Dit maakte een aantal procedures die moeten controleren of variabelen niet te groot of te klein zijn eenvoudiger. Waar de impulsflux uw gebruikt wordt gaat dat meestal via een hulpvariabele, uwo, die dan gelijk gemaakt wordt aan -1.0*uw.old. De +/- tekens voor de betreffende termen in de source code moeten dus nog steeds overeenstemmen met de theorie.

Nadat deze wijzigingen waren ingevoerd bleek het model zeer stabiel. Het crashte niet meer, ook niet achter steile heuvels.

4.2 Redistributie binnen één simulatiestap

Eigenlijk geldt voor alle processen dat ze veranderen als binnen een tijdstap de variabelen te veel veranderen. De redistributie van turbulente kinetische energie (TKE) over de verschillende richtingen is het proces met de kleinste tijdschaal, daarom mogen daar eerst problemen verwacht worden.

Om onstabiliteiten die hiervan het gevolg zouden kunnen zijn is het model op twee manieren aangepast. De stapgrootte die vast was is variabel gemaakt en de drukterm is opgesplitst in twee termen, namelijk de rapid distortion en de redistributie ten gevolge van turbulente snelheidsfluctuaties.

4.2.1 Variabele stapgrootte

Het probleem dat modellen zich onstabiel gaan gedragen als de simulatiestap te groot wordt doet zich bij alle processen voor. Daarom is gekozen voor een algemene oplossing. Zodra een van de parameters (U, uu, vv, ww, qq, uw, eps) significant verandert (hier is significant gedefinieerd als meer dan 2%), dan wordt de stapgrootte aangepast.

Het model rekent dus niet meer met een vaste stap, maar met een variabele stapgrootte. In H1LL_CFG.CFG kan hiervoor een maximale waarde worden opgegeven, bijv. "STEPSIZE = 0.20" (eenheden in meter). De opgegeven stapgrootte is een maximum. Zonodig wordt de stapgrootte verkleind.

De variabele NSTEP is verwijderd uit de configuratiefile, daar het aantal stappen nu per definitie variabel is.

4.2.2 Rapid Distortion

De rapid distortion is de produktie of destructie van turbulentie (uu, vv, en/of ww) en van impulsfiux (üw) doordat dicht bij de grond waar een flinke windsnelheidsgradiënt aanwezig is de turbulente wervels vervormd worden. De hiervoor benodigde theorie is ontleend aan Zeman & Jensen (1987) en Zeman & Tennekes (1975) en is besproken in het hoofdstuk theorie (paragraaf 2.5.1, 2.5.2).

Voor de turbulente kinetische energie (uu, vv, ww) heeft rapid distortion tot gevolg dat een deel van de geproduceerde turbulente zeer snel herverdeeld wordt over de verschillende richtingen. In het model gebeurt dat binnen één simulatiestap. Voor de impulsfiux (üw) wordt een deel van de produktie teniet gedaan.

NB; Voor de impulsflux was al een rapid distortion term opgenomen. Deze is niet aangepast aan de door Zeman & Jensen voorgestelde parameterisatie van de rapid distortion. De belangrijkste reden is dat deze term een parameter ( C Rapid) bevat die door het model automatisch aangepast wordt om in de ongestoorde situatie een goed windprofiel te genereren, dat bovendien in overeenstemming is met het gemodelleerde üw-profiel. Aanpassing van deze term vergt daarom aan aanzienlijk ingreep in het model.

4.3 Fout laterale wind

In het model speelt de laterale component van de windsnelheid een betrekkelijk geringe rol, omdat in de richting dwars op de lengterichting waarin gerekend wordt homogeniteit wordt verondersteld. Toch zijn er enkele termen waarin de laterale component terugkomt. Daarbij gaat het steeds om de afwijking van de laterale windsnelheid van de component van de geostrofische windsnelheid dwars op de modelrichting (V-Vg).

ZJ87 beschreven de geostrofische afwijking als functie van de Coriolisparameter ( f ) , de wrijvingssnelheid (U»), de hoogte (z) en de grenslaaghoogte (H):

/ ( F - F ^ ) = 3 uil

H

_z \ 2 [m (28)

Vaak wordt aangenomen dat de grenslaaghoogte evenredig is met de windsnelheid en omgekeerd evenredig met de Coriolisparameter (H ~ U/f, Tennekes & Lumley 1972). ZJ87 veronderstelde een evenredigheid met U , en de reciproke van de Coriolisparameter, waarbij de Von Karman constante de evenredigheidsconstante is.

/f = 0.4 [m] (29)

ƒ

F i g u u r 3 : Verband tussen grenslaaghoogte en

wrijvingssneiheidvoor de neutrale grenslaag Wanneer de uitdrukking voor dc grcuslaaghoogte

wordt ingevuld en vergelijking (28) wordt enigszins herschreven volgt:

F - F -7.5 1

-0.4

[m s-'] (30)

Dit schaalt met U*, in ieder geval in de onderste helft van de grenslaag. Op wat grotere hoogte is de schaling alleen correct als de hoogte ook wordt geschaald, en wel met de grenslaaghoogte (H = 0.4 U»/f).

In het Zeman en Jensen model is bij de implementatie van de bovenstaande theorie de van afliankelijke grenslaaghoogte vervangen door een constante grenslaaghoogte z, (inversiehoogte). Zodoende verandert (28) in:

F - F ( 1

\ 2

m s -21 (31)

Het is eenvoudig te zien dat dit niet schaalt met U«, maar met U*^. Dit heeft als gevolg dat bij hoge een te grote afwijking van de geostrofische wind wordt berekend. De vergelijking voor vv bevat in de produktieterm de gradiënt van de geostrofische afwijking (d(V-Vg)/dz). Deze wordt nu erg groot. In de eerste plaats wordt de geostrofische afwijking te groot berekend, maar bovendien wordt gerekend met een te kleine grenslaaghoogte.

In het HILL-model is gekozen voor een mengvorm. Omdat ook alle andere termen werken met een vaste grenslaaghoogte is wordt de hoogte nog steeds geschaald met de constante inversiehoogte, maar voor de grenslaaghoogte die buiten de haken voorkomt is wel de theoretische grenslaaghoogte gebruikt. De uitdrukking voor de geostrofische afwijking wordt dan als volgt:

\ 2

10000:

U* = 0.60 m/G

Zi = 1500 m 1000:

Figuur 4 toont liet verschil tussen de drie benaderingen. Te zien is dat dicht bij de grond vergelijking (32) hetzelfde resuhaat geeft als de theoretische benadering. Op grotere hoogte komen de waarden van V-Vg door de schaling van de hoogte met een te geringe grenslaaghoogte op een iets lager niveau terecht. De gradiënten zullen hierdoor iets te groot zijn, maar veel minder dan bij de Zeman en Jensen variant.

Figuur 5 laat zien wat een onjuiste berekening van de laterale component van de geostrofische afwijking voor gevolgen heeft

voor de laterale component van de turbulente F i g u u r 4: Geostrofische afivijking als functie van de hoogte

kinetische energie, vv. De verschillen zijn het Theoretisch (vette lijn) zoals gebruikt in Zeman&Jensen

grootste in de bovenste helft van de (stippellijn) en zoals gebruikt in het HILL-model (dunne

grenslaag.

-6.0 -5.0 -4.0 -3.0 -2.0 -1.0 Geoolrofioche afwijking V-Vg (m/o)

getrokken lijn).

Dicht bij de grond wordt de produktie van turbulente kinetische energie voornamelijk bepaald door de windsnelheidsgradiënt (dU/dz) en de impulsflux (üw) en ontvangt uu de meeste energie door redistributie vanuit uu (deze produktie van TKE vindt namelijk plaats in de uu component). Op iets grotere hoogte begint de gradiënt van de geostro-fische afwijking (d(V-Vg)/dz) mee te tellen in de produktie van vv. Het is dan ook niet verwonderlijk dat de Z&J variant veel te grote waarde voor vv voorspelt in de bovenste helft van de grenslaag, met name omdat daar ook de dissipatie (een verliesterm) klein is.

u * = 0.80 m/o Zi = 1500 m

0.5 1.0 1.5 LatQratQ variantiQ v v (m2/Q2)

F i g u u r 5 : Verticaal profiel van de berekende laterale variantie, vv, voor correctie in de laterale wind (stippellijn) en na correctie (getrokken lijn).

Een te hoge waarde van vv in de bovenste helft van de grenslaag werkt via de redistributie ook door in de andere termen van de turbulente kinetische energie (uu en ww). Omdat alle vergelijkingen gekoppeld zijn heeft dit artefact ook gevolgen voor de impulsflux (uw) en voor de windsnelheid (U). Via de transportterm in de diverse vergelijkingen en termen die afhangen van verticale gradiënten diffunderen de fouten vervolgens naar beneden. Vandaar dat de correctie nodig was.

4.4 Effecten van de kromtestraal op de stroming

De meeste grenslaagmodellen zijn ontworpen voor vlak terrein. Is het terrein niet vlak, dan zullen de stroomlijnen een kromming gaan vertonen. Finnigan (1988) laat zien dat het effect van een concave kromming (bijv. duinvoet) lijkt op het effect van een onstabiele temperatuuropbouw in de grenslaag, nl de produktie van extra turbulentie. Een convexe

kromming (dus bij de top van het duin) heeft een effect heeft dat vergelijkbaar is met een stabiele temperatuuropbouw, namelijk een onderdrukking van de turbulentie.

In de in Van Dijk et al. (1995) besproken versie van het HILL-model waren de effecten van de kromtestraal op de stroming en de turbulentie wel besproken in de theorie, maar niet meegenomen bij de implementatie van het model. Het opnemen van de krommingseffecten gaf aanvankelijk problemen.

Allereerst is de theorie van Zeman en Jensen (1987) in detail vergeleken met hetgeen Finnigan (1988) over dit onderwerp publiceerde. Doordat er aanzienlijke verschillen in notatie waren was deze vergelijking niet echt eenvoudig. Nadat echter de theorie van Finnigan was omgewerkt naar de notatie van Zeman en Jensen bleken beide identiek te zijn.

Alle termen waarin de kromtestraal zijn nu opgenomen in het model (zie ook paragraaf 2.5 en 2.6). Aangezien door de kromtestraal effecten extra turbulentie wordt geproduceerd, komt een deel hiervan weer in aanmerking voor een snelle herverdeling via de rapid distortion (zie paragraaf 2.5.4). Ook dit effect is opgenomen.

Aanvankelijk gaf de invoering van de krommingseffecten nog problemen. Daarom kunnen in HILL_PAR.PAR twee termen opgenomen worden (PROD_CURVE en RAPID_CURVE), waarmee de kromtestraaleffecten op de produktie van turbulente kinetische energie en op de rapid distortion afzonderlijk verminderd kumien worden. Normaal zijn beide parameters gelijk zijn aan 1.0. Door deze parameters gelijk te maken aan 0.0 wordt het effect van de kromtestraal uitgeschakeld (waarden tussen 0.0 en 1.0 zijn ook mogelijk). Is het bestand HILL_PAR.PAR niet aanwezig, of zijn daarin geen waarden opgegeven voor de parameters PROD CURVE en RAPID_CURVE, dan wordt gerekend met de default waarde van 1.0 voor beide parameters.

Voor de berekening van de kromtestraal wordt steeds gebruik gemaakt van drie verticale profielen PRl, PR2 en PR3, die ieder de hoogte van alle stroomlijnen bevatten op de punten x l , x2 en x3. De hoogte van een willekeurige stroomlijn op deze punten is nu z l , z2 respectievelijk z3. We kunnen nu over beide trajecten de eerste afgeleide berekenen als dz/dx,2 = (z2-zl)/x2-xl) en dz/dxj^ = (z3-z2)/x3-x2). Vervolgens is de eerste afgeleide dz/dx gedefinieerd als het gemiddelde van deze 2. De tweede afgeleide is gedefinieerd als d^z/dx^ = (dz/dx32-dz/d,2)/((x3-xl)/2). Daarna is de kromtestraal berekend met (Zeman en Jensen

1987):

[ [dx)

Als de profielen P R l , PR2 en PR3 erg dicht bij elkaar liggen wordt de berekening van de kromtestraal gevoelig voor afrondfouten. Daarom worden de profielen genomen op punten die tenminste 10 cm uiteen liggen.

De procedure waarin de kromtestraal wordt berekend is bijgevoegd in bijlage 9.

3/2 - 1