GOLFVOORSPELLING : PARAMETERSCHATTING BIJ HET OPSTELLEN VAN EEN WISKUNDIG MODEL
J. MONBALIU K.U.Leuven Laboratorium voor Hydraulica
Mathematica/ modelling tor wave forecast has become an important tooi, both trom the economical and environmental point of view. To obtain a better understanding of the working of the souree terms in a third generation wind wave model, a framewerk has been build tor tuning unknown parameters in the tormulation of the torces involved in wave generation and wave dissipation through whitecapping. A direct minimization method is applied using a standard optimization routine trom the widely available NAG-Fortran subroutine library to minimize a cost tune-
1. 1NLEIDING
Golfprediktie kan een belangrijk hulpmiddel zijn om vaarroutes voor de scheepvaart aan de weersomstandigheden aan te passen en om zodoende ongelukken zoals die begin dit jaar met de olietanker 'Braer' bij de Shetland eilanden te voorkomen. Het gebruik van een wiskundig model is daarbij een handig en krachtig hulpmiddel (zie Van den Eynde & De Wolf 1990). Maar bij het vertalen van fysisch inzicht naar een mathematisch model moeten meestal een aantal parameters vastgelegd worden. Enerzijds is de preciese waarde van een fysische parameter niet altijd met de gewenste nauwkeurigheid gekend, anderzijds kan men de complexiteit van de fysische werkelijkheid niet altijd precies weergeven, zij het door een gebrek aan kennis, zij het door een gebrek aan middelen (bvb. rekentijd, voldoende metingen voor een stochastisch gebeuren, ...). Zo zal bij het vastleggen van het verband tussen de waterhoogte en het debiet over een overlaat, een debietscoëfficiënt moeten bepaald worden om bepaalde aannames bij het theoretisch model te 'corrigeren'. Men spreekt van ijking of kalibratie. Zo kan ook het verband gezocht worden tussen een variabele en de kombinatie van verschillende andere variabelen, bijvoorbeeld bij evaporatie zullen zowel de luchtdruk als de temperatuur en de windsnelheid een belangrijke rol spelen. In meer complexe gevallen kan het verband tussen input en outWater nr. 69 - maartlapril 1993
tion. The cost tunetion is defined as a weighted sum of squares to take advantage of the special structure of this formulation. Measured physical characteristics such as wave energy and peak frequency of the waves are expressed in tunetion of fetch through the Kahma and Calkoen growth curve laws. In the optimization exercise, we try to reproduce the fetch evafution of the total energy by tuning two elected parameters in the wind input and/or dissipation souree terms. Such a framewerk is useful in increasing our understanding of individua/ souree terms.
put niet door een relatief eenvoudige regressie bepaald worden, maar moet men de evolutie van het systeem volgen. In bepaalde gevallen, zoals bij de simulatie van de stroming in leidingsnetwerken, moet men een stelsel van al dan niet lineaire vergelijkingen oplossen. In andere gevallen moet men de differentiaalvergelijking die het tijdsafhankelijk en/of plaatsafhankelijk systeem beschrijft, oplossen. De waarde die men dan aan een bepaalde parameter in het stelsel van vergelijkingen of in de differentiaalvergelijking geeft, zal de volledige oplossing van het probleem beïnvloeden. Neem de bepaling van de ruwheid van de bodem van een rivier. Men kan dan bijvoorbeeld een waarde toekennen aan de ruwheidscoëfficiënt die de kleinste afwijking geeft tussen de berekende waterhoogte en de gemeten waterhoogte bij de verschillende meetstations langs de rivier over de duur van een meetcampagne. Het systeem van de rivier wordt daarbij beschreven door de formuleri ng van de wetten van behoud van massa en van behoud van momenturn (Houthooft, 1986; Van Erdeghem et al., 1992). Of bij het transport van polluenten zal men een diffusiecoëfficiënt vastleggen, zodanig dat het wiskundig model het verloop van de concentratie in bepaalde meetpunten het best benadert. Merk dus op dat de parameters die men wil kennen niet noodzakelijk meetbaar zijn, maar dat ze wel een invloed hebben op hetgeen gemeten wordt. In het voo rbeeld dat in dit artikel zal be-
schreven worden, gaan we op zoek naar gepaste parameters voor de brontermen van windtoevoer en van energiedissipatie door schuimkopjes, zodanig dat de oplossing van de energietransportvergelijking van door wind opgewekte golven zo goed mogelijk overeenkomt met wat in een fysisch experiment opgemeten werd. De juiste kennis van deze brontermen is belangrijk voor het korrekt voorspellen of bepalen van een windzee al of niet met deining. 2. OPTIMALISERINGSBENADERING
De voorbeelden in de inleiding tonen aan dat men het verschil tussen de waarden die men met een wiskundig model bepaalt en waarden die in werkelijkheid optreden , zo klei n mogelijk wil maken. Daartoe formuleert men meestal een kostfunktie die de som is van kwadraten :
V(A) =
2
n
E
W; (
Y;(X,t,A) - Y;(X,t) )
i•l
Y; en Y, zijn de gemodelleerde en waargenomen toestandsvariabelen, x geeft de ruimtelijke lokatie weer, t is de tijd en A is de vector die de ongekende te bepalen parameters bevat. Daarbij kan aan de toestandsvariabelen, naar gelang het vertrouwen dat men erin heeft, een bepaald gewicht w, toegekend worden . Indien de modelwaarden Y; bekomen worden door 51
spectrum (zie figuur 1) is wellicht de meest gebruikte formulering en wordt gegeven door:
Fig. 1 :Het JONSWAP spectrum (naar Hasselmann et al., 1973)
waarbij
y = Emox ,...-,
PM
aa
EMAX
(/)
r
N
E
if-f.'\2 exp(--r-) 2(flf.p2
L......J
Q) I
0'1
I
'-
I
Q)
c
w
frekwentie
fp
het oplossen van een partiële differentiaalvergelijking, kan het probleem als volgt samengevat worden : parameters A~ (a 1,a2 , .. ,ai, .. ,an)
[H z]
aire golven met verschillende amplitude. Deze superpositie wordt geschreven in de vorm van het verplaatsi ngsspectrum :
V
I
A is invoer partiële differentiaalvergelijking
00
11
V
bereken het minimum van de kostfunktie V
E(w)dw
3.2 De energietransportvergelijking - - - -- -- - -
Naargelang het aantal parameters die als invoer dienen en die moeten bepaald worden, kiest men bepaalde methodes. Is het aantal parameters groot dan gaat de voorkeur naar de zogenoemde 'adjoint models' (Snyder et al., 1992). Deze methode wordt gekenmerkt door een additioneel wiskundig model (i.e. de 'adjoint model') dat gedreven wordt door het verschil tussen de modelresultaten en de waarnemingen en dat omgekeerd in de tijd opgelost wordt. De methode heeft als voordeel dat de gradient van de kostfunctie kan bepaald worden met het equivalent van twee modelruns, maar heeft als nadeel dat een tweede wiskundig model moet opgesteld worden. Indien het aantal te bepalen parameters beperkt is en de dynamische situatie die men beschrijft relatief eenvoudig is, kan beter een rechtstreekse methode gebruikt worden. Men hoeft geen additioneel model te maken, alhoewel het additionele model niet noodzakelijk ingewikkelder is dan het oorspronkelijke. Bovendien kan men gemakkelijk gebruik maken van bibliotheken waarin mathematische subroutines zitten die dergelijke problemen oplossen. Het eigen mathematisch model fungeert dan als subroutine van de kostfunktieberekening. Deze laatste oplossingsmethode werd gekozen voor het hieronder in meer detail weergegeven voorbeeld. 3. TOEPASSING OP DE ENERGIETRANSPORTVERGELIJKING VAN WINDGOLVEN
3.1 Beschrijving van het zeeoppervlak Het wateroppervlak kan gezien worden als een superpositie van een groot aantal line52
_____J
met Ç : verplaatsing van het wateroppervlak (t.o.v. water in rust) E(oo) : frekwentiespectrum ; oo = 2nf en f is de frekwentie F(k) : golfgetalspectrum waarbij k (een vector) het golfgetal is Het spectrum van windgolven behoudt onder invloed van de verschillende effecten die erop inwerken zijn vorm (self-similarity) en kan daardoor beschreven worden met een beperkt aantal parameters. Mees tal werkt men in het frekwentiedomein. Een aantal formuleringen zijn welgekend. Het JONSWAP spectrum (Hasselmann et al. 1973), heeft een t·5 -staart. Het Pierson-Moskowitz (PM) spectrum (Pierson en Moskowitz, 1964) is een bijzonder geval van het JONSWAP spectrum en beschrijft het zeeoppervlak bij volgroeide zeegang. Het Toba spectrum en het equivalente spectrum van Baltjes (Toba, 1973; Baltjes et al. 1987) en het Donelan spectrum (Donelan, 1985) hebben een f·4 staart. Het Toba spectrum heeft daarbij een extra parameter, namelijk de wrijvingsnelheid (u*) in de lucht. Het JONSWAP
I
g = 9.81 = valversnelling [m 2/s] f = frekwentie De schaal parameters zijn : fP= frekwentie van de piek van het spectrum <Xp= Phillips constante De vorm parameters zijn : y = piekfaktor (peak enhancement factor) == 3.3 cra = breedte linkerzijde piek ~ cr = cra = 0.07 indien f < fP crb = breedte rechterzijde piek ~ cr = crb = 0.09 indien f > fP Golven hebben niet enkel een frekwentie maar ook een richting (8). Een aantal formuleringen voor de spreiding van de energie over de verschillende richtingen , kan men vinden in Monbaliu (1992).
voor graviteitsgolven De energie die men terugvindt in de golven van de zee zijn het gevolg van de verschillende processen die erop inwerken. In deze studie zijn we enkel geïnteresseerd in golven opgewekt door wind en er zal dan ook aangenomen worden dat de kracht van de wind de enige bron van energie is. Energie kan verloren gaan door breking van de golven (branding en schuimkopjes), en door middel van wrijving op de bodem. Door de aanname dat het wateroppervlak bestaat uit de superpositie van een groot aantal lineaire (=sinusoïdale) golven, worden hogere orde termen verwaarloosd. Hasselmann (1962, 1963a en 1963b) heeft theoretisch aangetoond dat een aantal combinaties van golven onderling energie kunnen uitwisselen. Deze niet-lineaire wisselwerking, ook wel resonante wisselwerking genoemd, tussen golven heeft geen energieverlies noch energiewinst tot gevolg, enkel een verschuiving van de energie van bepaalde trekwenties en richtingen naar andere trekwenties en richtingen. Bovenvermelde processen, met name de inwerking van de wind, het verlies van golf-
-
Fig. 2 : Het Phillips mechanisme
a 6
uc(kn)
e
-
Water nr. 69 - maart!apri/1993
Fig. 3 : Het Mi/es mechanisme : a. stroomlijnen zonder gradi(jnt in de snelheid. b. met gradiMt in de snelheid en meebewegend met het golfoppervlak
een verhoogde druk (stroomlijnen verder uit elkaar) en op de lijzijde een verminderde druk (stroomlijnen dichter bij elkaar) (zie Figuur 3). De golven groeien en worden steiler. Doordat de golven steiler worden, wordt het mechanisme versterkt, zodat er sprake is van een instabiliteit. Wiskundig formuleert zich dat in een exponentiële groei. Meestal wordt enkel het Miles-mechanisme weerhouden in wiskundige modellen, zodat de windtoevoer kan geschreven worden als :
s. =
öE(j,O) Öt
111
energie door breking en door bodemwrijving en de niet-lineaire wisselwerking tussen verschillende golfcomponenten vormen de basis voor de huidige kennis aangaande··de energiebalans van golven opgewekt door wind. In diep water, waar de waterdiepte voldoende groot is zodat de bodemwrijving en andere (on)diepte-effecten zoals refractie kunnen verwaarloosd worden, wordt de energietransportvergelijking in de x-richting gegeven door :
öE(j,8)
öE(j,8)
öt
öx
SIQI(j,8)
1
van het wateroppervlak. Men krijgt dan een soort resonantie-effekt. Als deze drukschommelingen nu zouden voortbewegen star verbonden aan het wateroppervlak, zou men een voortdurende toevoer van energie krijgen en een exponentiële groei van de golven. Dit is niet het geval en uiteindelijk krijgt men slechts een lineaire groei van de golfenergie in functie van de tijd. Het Phillips-mechanisme is het voornaamste proces voor de initiële groei van windgolven. Een tweede mechanisme is het zogenoemde Miles-mechanisme. Doordat er een gradiënt is van de snelheid ('shear flow') in de atmosferische grenslaag, krijgt men op de loefzijde van het golfoppervlak
=
11
r
w EU,O) V'
Een aantal formuleringen voor de koppelingscoëfficiënt ~~ zijn in de literatuur gegeven. Men kan ze onderverdelen in twee categorieën, met name deze die lineair zijn in een karakteristieke snelheid , aangeduid als het Snyder type naar de metingen van Snyder et al. (1981), en deze die kwadratisch zijn, aangeduid als het Stawart type (Stewart, 1974). Aangezien het tweede type bij de numerieke uitwerking minder goede resultaten gaf (Monbaliu, 1992), weerhouden we hier enkel het eerste type als voorbeeld :
Su.if.9)
=
P. 28u, 0.25 a, Pw (a2 -c-cos9 - I) w E(/,9)
De coëfficiënten a, en a2 hebben in de literatuur een waarde 1, maar zullen hier als parameters beschouwd worden. Merk op
---------- -- ---
= S111(j,8)
+ Sïn(j,8) + Sdisl(j,8)
Fig. 4 : a. schematische weergave van het JONSWAP-experiment. b. evolutie van het verplaatsingspectrum met toenemende strijklengte (fetch-limited)
~meetpunt In woorden betekent dit dat wat in de tijd gewonnen wordt aan energie per eenheid van oppervlakte en per eenheid van gewicht, gelijk is aan wat erin gebracht wordt door de wind (S,n), wat verloren gaat door breking (Sdiss> en wat uitgewisseld wordt door de resonante wisselwerking (Sn1), verminderd met de energiehoeveelheid die weggevoerd wordt aan een snelheid die gelijk is aan de groepsnelheid van de golven (Cgx).
fetc h
12
r-----------------------------~
---- x
10 X
3.3 De brontermen Het effect van de wind op de golven en het daarmee gepaard gaande groeiproces worden voornamelijk door twee mechanismen verklaard (Phillips, 1977). Een eerste mechanisme is het Phillips-mechanisme. De wind wordt voorgesteld als bestaande uit wervels, die worden meegevoerd. Kleine wervels worden dicht bij het wateroppervlak meegevoerd met een kleine snelheid. Hoe groter de wervels, hoe hoger die zich in de atmosferische grenslaag bevinden en hoe sneller die worden getransporteerd (zie Figuur 2). De wervels creëren een drukverdeling op het wateroppervlak. De golflengten van deze drukverdeling kunnen samen vallen met de golflengten Water nr. 69- maartlapri/1993
8
26 km
= 64
--- X
6 4
2
0.10
0.20
0.30
0.40
0 .50
km
128 km
- ..-
X
= 257 km
-
X
1278
km
- - X = 2552
km
-
_ _x_ =
6374
_:j 53
dat hier gekozen werd voor de wrijvingssnelheid als karakteristieke snelheid, in overeenstemming met het idee dat de schuifspanning in de atmosferische grenslaag de drijvende kracht is op het wateroppervlak (Komen et al., 1984 ; Janssen et al., 1987). Voor de dissipatie van golfenergie door schuimkopjes zijn de ideeën van Hasselmann (1974) overgenomen. Dit resulteert in een bronterm die evenredig is met het energiespectrum en evenredig met he t kwadraat van de frekwentie. De evenredigheictsfactor hangt enkel af van het golfgetal en geïntegreerde spectrale grootheden zoals de gemiddelde steilheid van de golven. Komen et al. (1984) werkten deze bronterm in meer detail uit :
Sdi#if,O)
waarbij gemiddelde frekwentie [rad/s] a. : steilheid van de golven a.PM : steilheid van golven bij volgroeide zeegang De coëfficiënt c1 beïnvloedt het algemeen niveau van de demping , de coëfficiënt n verschu ift de demping t.o.v. de piek van het spectrum en de coëfficiënt m geeft de afhankelijkheid van de demping t.o.v. de steilheid van de golven weer. De coëfficiënten c,. m en n werden in deze studie als parameters beschouwd. Dank zij de studie van Komen et al. (1984), kunnen we als richtwaarden 3.33 1o·5 voo r c, , en 2 voor m en n nemen. Tenslotte is er nog de resonante wisselwerking van de golven. Stellen van vier golven kunnen energie uitwisselen, als aan bepaalde resonantievoorwaarden voldaan is. De hoeveelheid energie die uitgewisseld wordt, kan berekend worden met behulp van een vrij ingewikkelde maar vooral rekenintensieve integraal. Daarom ook werd het 'discrete interaction approximation' algoritme van Hasselmann et al.(1985) gebruikt om de rekentijd te beperken.
w:
het JONSWAP spectrum , samen met de uitwerking van de groeicurven voor de parameters die dit spectrum beschreven. Ook anderen hebben experimenten gedaan om de evolutie van het spectrum te bekijken, o.a. Kahma (1981) in de Botnische Golf, en Donelan ( 1985) in het Ontariomeer. Kahma en Calkoen (1992) vonden het daarbij merkwaardig dat de groeicurven voor de energie tussen de verschillende experimenten vrij grote verschillen vertoonden. Daarom brachten ze de data van bovenvernoemde experimen ten bijeen en heranalyseerden deze. Dit resulteerde in nieuwe groeicurven , hier weergegeven in dimensieloze vorm , waarbij het eigenlijk niet duidelijk is wat als karakteristieke snelheid moet genomen worden om gemeten golfgrootheden dimensieloos te maken. Wel is min of meer duidelijk dat het ofwel de windsnelheid buiten de grenslaag moet zijn (uJ, ofwel de wrijvingsnelheid die een maat is voor de schuifspanning op het wateroppervlak. Er is intuïtief een voorkeur voor de wrijvingssnelheid u•. Harde bewijzen daarvoor bestaan echter niet (Janssen et al. 1987):
Dit resulteerde in 1973 tot het JONSWAP rapport (Hasselmann et al., 1973). Belangrij k in dit rapport was de formulering van 54
g
S.6 10-3
1.1 lol «pji
= 0.0081
3.5 De kostfunktie Het criterium om waarden toe te kennen aan bepaalde parameters in de brontermen (vgln. (7) en (8)), bestaat erin dat het verschil tussen gemeten grootheden en waarden die men kan vinden door berekeningen zo klein mogelijk moet zijn. Daarbij werd er vanuit gegaan dat een gewogen kleinste kwadraten benadering een goede keuze is. Volgende kostfunktie werd voorgesteld:
E. = 2.4 I0-3 x~· 78
.r,;
=
0.358
x:0· 244
"• =f"u • --
E. =Eg2
Jp
u.4
g
x =xg •
2
u.
Kahma en Calkoen ( 1992) geven geen groeicurve voor de Phillips constante ap, daarom wordt hier de groeicurve uit het JONSWAP rapport overgenomen:
ap
De subscript i staat voor de strijkl engte (fetch) x,, subscript b voor berekende en subscript g voor gemeten groothede n . Voor de gemeten waarden zullen de gemeten groeicurven (vgln. 9, 10 en 11) genomen wo rden . De berekende waarden zijn de resultaten van het oplossen van de energietransportvergelijking (vgl. (5)). Om een idee te krijgen van het gedrag van
=0.35 x;·22
Voor volgroeide zeegang werden volgende waarden aangenomen (ct. Komen et al.(1984); PM verwijst naar Pierson-Moskowitz (1964)): Fig. 5 : Detail van het verloop van de kostfunktie
3.4 De evolutie van het golfoppervlak onder invloed van wind Om vat te krijgen op de evolutie van het spectrum van windgolven onder invloed van wind , we rd op het eind van de jaren zestig een grote meetcampagne opgezet vóór de kust van het eiland Sylt (grens Duitsland en Denemarken). Deze lokatie was gekozen omwille van de vrij rechte kustlijn. Loodrecht op de kust had men een aantal meetposten uitgezet zodanig dat men bij aflandige wind de evolutie van het energiespectru m met de afstand tot de kust kon volgen. Dit wordt kwalitatief weergegeven in Figuur 4. Merk op dat in dit geval de energietransportvergelijking (vgl. 5) eenvoudiger wordt aangezien de tijdsafhankelijke term ~wegvalt.
~ = /p11-"· = JpJA
0.035 0 .034 0 .033
>
0 .032
CD
0 .031
:0:::: ~
-.. c
0 .030
Cl)
0.029
::I
0
~
0 .0 2 8 0 .027 0.026 0.025 28
29
30
31
32
33
34
35
dissipatie constante c
Wa ter nr. 69- maartlapril 1993
Tabel 1 : Opt1male paramete/Waarden in de brontermen voor de reproduktie van de Kahma en Calkoen groeicurven
parameters
a, 0.61
(a, ,c,)
a2 1
--~-
1.36
~
c 1(*105 ) 2.13
m -- 2 - -
kostfunctie 0.005
-~
~
(a2 ,c 1)
1
...
(c,,m)
10.7 ~ 12.2
1.42
1
Komen et al.
2
0.003
3.20
2.65
0.004
3.33
2
0.060
(1984)
r-
Fig. 6 : Groeicurven bij gebruik van optimale parametercombinaties voor : a. dimensieloze energie. b. dimensieloze piekfrekwentie. c. Phillips constante
G)
"ii ...
1000 -
G)
c
-
Kahma en Calkoen
G) G)
N
0
Q)
ïii c G) E
100
al cl
0.61 2.13
a2 c1
1.42 12.2
c1 m
3 .20
= 2.65
'5
10 10'
10 6
10 7
10 9
10°
dimensieloze fetch
G)
:;;
c
G)
~
-..
:\.,
.:,(. G)
~
. :,(. G)
...~
'
·a. G)
0 .01
'
N
0
Q)
·o; cG)
"- '
" N.
"",, ..,,
...
K a t-ma en CaiKoen
3.6 Beschikbare routines
a1 c1
0 .6 1 2.13
a2 c1 -
1.42 12.2
c1 m
3 .20 2 .65
Alle berekeningen werden uitgevoerd op de IBM -3090 van het rekencentrum van de K.U.Leuven. Het rekencentrum beschikt bovendien over de mathematische subroutine bibliotheek NAG. Daarin hebben wij gekozen voor de routine E04FCF, omdat deze routine een aantal opties had om het optimaliseringsproces van dichtbij te volgen. De routi ne E04FCF is gebaseerd op de Giii-Murray methode (zie Gill en Murray (1976) , Gi ll e t al. (1981 ), en Sca les (1985)). Deze methode maakt gebruik van het feit dat bij een kleinste kwadraten methode de Hessiaan, de matrix met partiële afgeleiden van de tweede orde, kan benaderd worden met een uitdrukking in de Jacobiaan die enkel partiële afgeleiden van de eerste orde nodig heeft. Dit heeft een aanzienlijk rekenvoordeel tot gevolg. De optimaliseringsroutine berekent opeenvolgende benaderingen van de optimale parametervector : A (k+1) = A (lll + fl.lll(l.kl waarbij p(k) : de zoekrichting a (k) : staplengte in die zoekrichting. Het programma bepaalt dus eerst een zoek rich ting en daarna gaat het in de zoekrichting een voorlopig minimum bepalen door de staplengte aan te passen. Vanuit dit voorlopig minimum wordt de procedure dan herhaald. Bij het gebruik van de opgestelde optimaliseringsprogramma's
--, ----------'
E
"'
- ''
'5
0.004 10°
10 6
10 7
10°
109
dimensieloze fetch
-G)
c
.JONSWAP
(IJ
a1 cl
0 .6 1 2.13
a2 c1
1.42 12.2
IQ
c 0 0
(IJ
~
0.01
c1 3.20 m = 2.65
:.2
0.
0.004
L-~~~.......-~~~.......-~~~........-~~~....J
10°
10
6
10
7
dimensieloze fetch
Water nr. 69- maartlapril 1993
101
10 1
__,
de kostfunktie, zijn een aantal gevoeligheidstests uitgevoerd waarbij slechts één parameter varieert. Wij noteren volgende bevindingen : 1 . De kostfunktie waarin enkel de verschillen in de energiehoeveelheden in het spectrum meetellen, is het meest gevoelig, zodat de kostfunktie kan gereduceerd worden tot de eerste term (w2 ,, en w3 ,, gelijk aan 0). 2. De minimale waarden voor een kostfunktie gedefinieerd om de totale energie van het golfspectrum te verkrijgen, liggen bij andere parameterwaarden dan voor een kostfunktie om enkel de piekfrekwentie of de Phillips-constante aP optimaal te benaderen. 3. De kostfunktie is ongevoelig voor veranderingen in de parameter n van de dissipatieterm. Grotere waarden (> 2.5) van de parameter n, leveren echter aanzienlijke problemen op voor het berekenen van de Phillips-constante. 4. De kostfunktie is niet erg gevoelig aan veranderingen in de dissipatie-parameter c 1, zeker niet in de buurt van zijn optimale waarde. 5. De groeicurven blijven bij een verandering in de waarde van om het even welke beschouwde parameter, qua vorm ongeveer dezelfde. 6. Kleine verschillen in een welbepaalde pa rame te r l eiden tot willekeur ige schommelingen (grootte-orde 1 0"3 ) in de kostfunktie (zie Figuur 5). Dit is enerzijds te wijten aan het numeriek schema (aanpasbare integratiestap, staat voor de hoge frekwenties,... ) en anderzijds aan het discrete domein van richtingen en trekwenties waarin gewerkt wordt.
55
zijn volgende zaken belangrijk : 1) het gebruik van schaallaktoren voor de te optimaliseren parameters 2) het gebruik van aanpasbare parameters in de NAG-subroutine bibliotheek 3) het interpreteren van de NAG convergentiecriteria Het gebruik van schaallaktoren voor de parameters is noodzakelijk enerzijds omdat de te kalibreren parameters voor een goede werking van de routines van dezelfde grootte-orde moeten zijn, anderzijds omdat de eerste partiële afgeleiden numeriek moeten benaderd worden. Deze afgeleiden worden op de IBM-3090 van het universitair rekencentrum van de K.U.Leuven via een voorwaartse absolute stap van de grootte-orde 1o-s berekend. Hoger is ook reeds vermeld dat kleine veranderingen in de parameters aanleiding geven tot schommelingen in de kostfunktie. Daarom moet men er zeker van zijn dat de absolute stap voor het berekenen van de afgeleiden een betekenisvol verschil geeft in de waarde van de kostfunktie. In de praktijk komt het voor deze toepassing erop neer dat die absolute stap enkele procenten van de parameterwaarde bedraagt. Het gebruik van de door NAG-bijgeleverde parameters, kunnen wij voor deze toepassing als volgt samenvatten. Bij voorkeur wordt de zoekrichting opnieuw berekend en doet men weinig moeite om het minimum in de huidige zoekrichting te bepalen. Het maximaal aantal keren dat de kostfunktie mag berekend worden, wordt best aan de lage kant gezet (20 tot 50 maal het aantal optimaliseringsparameters). De NAG-routines hebben hun eigen convergentiecriteria. Met een oordeelkundig gebruik van de NAG-tolerantieparameter en een goede keuze voor de schaalfaktor van de kostfunktie kan aan deze criteria voldaan worden. Het is ook mogelijk om zelf een convergentiecriterium in te bouwen en het programma met een extra mededeling stop te zetten. Daarbij moeten wij nog vermelden dat het soms gebeurt dat er geen convergentie optreedt omdat het optimaliseringsprogramma niet in staat is bij een nieuwe zoekricht ing de minimale kostfunktie uit de vorige iteratie te verbeteren. Het optimaliseringsprogramma interpreteert dit als een discontinuïteit en wil daaroverheen springen. Door onze keuze van schaallaktoren wordt die sprong, relatief gezien, veel te groot (= veel groter dan onze parameterwaarden). De gevonden waarden voor de parameters en de gevonden waarde voor de kostfunktie uit de vorige iteratie kunnen meestal als goede resultaten geïnterpreteerd worden.
3.7 Een numeriek voorbeeld Ter illustratie worden hier de parameters bepaald in de wind input (vgl. 7) en in de dissipatie bronterm (vgl. 8) om zodoende de Kahma en Calkoen-groeicurven voor de energie(vgl. 9) en de asymptotisch energieniveaus (vgl. 11) zo goed mogelijk te reproduceren. Dit geeft voor parameters twee aan twee genomen de resultaten van Tabel 1. Ter vergelijking werden de richtwaarden van Komen et al.(1984) mee op-
56
genomen. Deze richtwaarden geven een vrij goede fit voor de JONSWAP-groeicurven (Monbaliu, 1992), doch geven minder goed de groeiende zeegang weer zoals die door Kahma en Calkoen (1992) geïnterpreteerd werd. In Tabel 1 zien wij dat de 'minimale' kostfunkt ie voor de verschillende parameterkombinaties dicht bij elkaar liggen. De verkregen groeicurven zijn voorgesteld in Figuur 6. Alhoewel enkel geprobeerd werd de groeicurve voor de totale energie te benaderen, worden zowel de groeicurve voor de piekfrekwentie als voor de Phillips' constante goed benaderd. De parametercombinatie (a2 ,c,) scoort het best'en in Figuur 6 kunnen wij ook zien dat zijn asymptotisch gedrag heel goed is. Dit geldt eveneens voor de parametercombinatie (c,,m). Voor de parametercombinatie (a 2 ,c 1) zien we bovendien dat afhankelijk van de startpositie andere 'optimale' combinaties eenzelfde waarde voor de kostfunktie geven. Dit wil zeggen dat de kostfunktie voor deze parametercombinatie vrij vlak is in een zone en het minimum bijgevolg niet precies bepaald kan worden . De parametercombinatie (a,,c,) doet het minder goed, wellicht omdat de input door de wind aanzienlijk kleiner is dan Snyder et al. (1981) vooropstelden door metingen. Men krijgt daardoor een verminderde dissipatie, die het evenwicht tussen de verschillende processen naar een andere limiet verschuift.
4. BESLUIT Het is mogelijk om op een snelle en efficiënte manier, doch voor een beperkt aantal parameters die het verloop van de volledige oplossing van een partiële differentiaalvergelijking beïnvloeden, waarden vast te leggen met een dirakte methode. Dit heeft als grote voordeel dat men standaard optimaliseringsroutines kan gebruiken waarbij enkel het bestaand wiskundig model nodig is als subroutine van de kostfunktieberekening. De grootste moeilijkheid is het vastleggen van een kostfunktie die garandeert dat de bekomen parameterwaarden uniek zijn en bovendien of daardoor de fysische werkelijkheid van het probleem weergeven.
J. MONBALIU K. U. Leuven Laboratorium voor Hydraulica de Croylaan 2, bus 4 300 1 Heverlee
REFERENTIES Battjes, J.A., Zitman T.J. en Holthuysen LH., A reanalysis of the spectra observed in JONSWAP, J. Phys. Ocean., 17, 1288-1295, 1987.
Gill, P.E .. Murray, W. en Wright, M.H .. Practical Optimization, Academie Press, 1981. Hasselmann, K., On the non-linear energy transfer in a gravity-wave spectrum. Part 1. General Theory, J. Fluid Mech. 12, 481-500, 1962. Hasselmann. K., On the non-linear energy transfer in a gravity-wave spectrum. Part 2. Conservalion Theorems ; wave- partiele analogy ; irreversibility, J. Fluid Mech., 15,273-281, 1963. Hasselmann, K., On the non-linear energy transfer in a gravity-wave spectrum. Part 3. Evaluation of the energy flux and swell-sea interaction lor a Neumann spetrum, J. Fluid Mech., 15, 385-398, 1963. Hasselmann, K., Barnett, T.P., Bouws, E., Carlson, H., Cartwright, D.E., Ewing, J.A., Gienapp, H., Hasselmann, D.E., Kruseman, P., Meerburg, A. , Müller, P., Olbers, D.J .. Richter, K. Sell, W. en Walden , H., Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project(JONSWAP), Dtsch. Hydrogr. Z., AS, 12, 95 p., 1973. Hasselmann, K., On the speetral dissipation of ocean waves due to whitecapping, Boundary Layer Meteorology, 6, 107-127, 1974. Hasselmann S. en Hasselmann, K., Computations and parametrizations of the nonJinaar energy transfer in a gravity-wave spectrum, Part I : A new method lor efficient computations of the exact nonlinear transfer integral, J. Phys . Ocean., 15, 1369-1377, 1985. Hasselmann S., Hasselmann, K., Allender, J.H. en Barnett, T.P., Computations and parametrizations of !he nonlinear energy transfer in a gravitywave spectrum, Part 11 : Parametrizations of the nonlinear energy transfer lor the application in wave models, J. Phys. Ocean., 15, 1378-1391 , 1985. Houthooft, R., Bepaling van hydraulische parameters in het Scheldebekken, eindwerk Toegepaste Wetenschappen K.U.Leuven, departement Computerwetenschappen, 1986. Janssen, P.A.E.M., Komen. G.J. en De Voogt, W.J.P., Friction velocity sealing in wind wave generation, Boundary Layer Meteorology, 38, 2935, 1987. Kahma, K.K., A study of the growth of the wave spectrum with fetch, J. Phys. Ocean., 11 , 15031515, 1981. Kahma, K.K. en Calkoen. C.J., Reconciling discrepancies in the observed growth of wind-generaled waves, J. Phys. Ocean., 22(12), 13891405, 1992. Komen, G.J., Hasselmann, S. en Hasselmann, K., On the Existence of a Fully Developed WindSea Spectrum, J. Phys. Oceanogr., 14. 8, 1984. Monbaliu, J., Wind and waves : lnvestigation of an optimization approach to parameter estimation, doctoraatsthesis, departement bouwkunde K.U.Leuven. 1992. NAG (Numerical Algorithm Group), lntroductory guide to NAG Fortran Library Mark 14, Oxford, Uniled Kingdom, 1990.
Donelan, M.A.. Hamilton, J.A. en Hui, W.H., Directional spectra of wind generaled waves, Phil. Trans. R. Soc. Lond., A.315, 509-562, 1985.
Phillips, O.M., The dynamics of the upper ocean, second edition, 1977.
Gill, P.E. en Murray, W., Algorithms lor the solution of the non-linear least-squares problem. NPL report NAC 71, National Physical Laboratory, 38 p., 1976.
Pierson, W.J. en Moskowitz , L., A Proposed Speetral Form lor Fully Developed Wind Seas Based on the Similarity Theory of S.A. Kitaigorodskii, J. Geoph. Res.. 69, 24, 5181-5190, 1964.
Water nr. 69 - maartlapril 1993
Scales, L.E. , Introduetion to non-linear optimization, Springer Verlag, New Vork, 1985. Snyder, R.L., Dobson F.W., Elliot J.A. en Long A.B. , Array measurements of atmospheric pressure fluctuations above surface gravity waves, J. Fluid Mech., 102, 1-59, 1981 . Snyder, R.L., Neu W.L., Long A.B. en de Voogt W.J.P., Analysis of field data and Completion of
Inverse-Medelling Experiments, A proposal submilled to the National Science Foundation in renewal of NSF grant OCE-8817273, 1992. Stewart, A.W., The air-sea momentum exchange, Boundary Layer Meteorology, 6, 151-167, 1974. Toba, Y., Local balance in the air-sea boundary process, Part 111, On !he spectrum of wind waves, J. of the Oceanographical Society of Japan, v.
29, 209-220, 1973. Van den Eynde, D. en De Wolf, P., Deiningsprediktie aan de Belgische kust, Water 52 , 163-167, 1990. Van Erdeghem, D., Troch, P.A. en De Troch , F.P., Mathematisch voorspellingsmodel van wasdebieten op de Belgische Maas. Deel 2 : Hydraulische Flood-routing, WATER 63, 1992.
15de Internationaal Congres over Irrigatie en Drainering Water is de bron van alle leven. De aanwezigheld in de wereld is echter niet optimaal verdeeld. In sommige landen is er een tekort, elders heeft men met een teveel te kampen. Soms is het mogelijk, zij het op zeer beperkte schaal, kunstmatig bij een ongunstige toestand in te grijpen. Waterzieke gronden worden gedraineerd en droge gebieden worden met een netwerk van kanaaltjes geïrrigeerd. Sommige landen hebben een rijke ervaring in het droog houden van gronden die soms belangrijk onder het zeeniveau liggen, andere konden door een gepaste irrigatie woestijnstreken tot vruchtbare zones ombouwen. De ICID, "International Commission on lrrigation and Drainage", heeft tot doel kennis en ervaring uit te wisselen. De ICID is een internationale wetenschappelijke en technische organisatie op het gebied van irrigatie, drainering, controle en regeling van stroombekkens. Zij heeft tot doel de ontwikkeling en toepassing van de kennis, de ingenieurswetenschappen, de landbouw, de economie en de sociale gevolgen van een goed waterbeleid te stimuleren; reeds zijn 8 landen aangesloten bij de internationale commissie waarvan de zetel in New Dheli is gevestigd. Elk aangesloten land heeft zijn eigen nationaal comité. Twee derden van de landen kunnen bij de ontwikkelingslanden worden gerekend. Naast de jaarlijkse bijeenkomsten van die internationale commissie wordt om de drie jaar een groot internationaal congres georganiseerd. In 1993 wordt dit congres gehouden in Nederland van 6 tot 11 september. Het gaat door onder voorzitterschap van Prof. Wil Segeren, Rector aan het "lnternationallnstitute for Hydraulic and Environmental Engineering" (IHE) te Delft en heeft het "waterbeheer in de volgende eeuw" als algemeen thema. Bij de deelthema's noteren wij o.a. de systeemanalyses, de micro-irrigatie, het impact van drainering en irrigatie op het milieu en de instrumentatie- en monitoring-systemen. Nederland heeft op het gebied van waterbeheer zeker heel wat te bieden en het congres zal dan ook de gelegenheid zijn om de verwezenlijkingen bij onze noorderburen, nu en in het verleden nader te leren kennen. Naast het eigenlijke congres wordt eveneens een internationale tentoonstelling gehouden in de grote hall van het Nederlands Congrescentrum in Den Haag. Deze tentoonstelling is open voor het grote publiek van 7 tot 11 september 1993 en zal over een oppervlakte van 8.000 m 2 een ruim beeld geven van de ontwikkeling in het huidige waterbeheer. Nadere gegevens in verband met congres en tentoonstelling kunnen worden verkregen bij het Holland Organisation Centre, 16 Lange Voorhout te 2514 EE Den Haag, tel: 0031703657850- fax: 0031703614846. H.R.
Onlangs is de v.z. w. WEL voor een nieuwe periode van drie jaar erkend als culturele vereniging, waardoor giften van 1.000,- BEF al fiscaal aftrekbaar zijn. Het fiscaal attest dat hiertoe door de v.z. w. WEL wordt afgeleverd, is een officieel document dat u bij het aangifteformulier van uw belastingen kunt voegen om te genieten van een belastingsvermindering, die tot meer dan 50 % kan bedragen van het bedrag dat u overmaakt aan WEL. De schenker wordt vermeld in de WEL -tijdschriften. Reeds 12 jaar ijvert de v.z.w. WEL ervoor een zo ruim mogelijke wetenschappelijk verantwoorde informatie te verstrekken aan een zo ruim mogelijk publiek, betreffende alles wat met het water, de energie en het leefmilieu verband houdt in de ruimste zin van het woord. Zodoende heeft WEL bijgedragen tot de bewustwording van de problematiek voor een gezonder leefmilieu.
Water nr. 69 - maartlapril 1993
57
