Numerische Modellierung von Flachwasserwellen — Besondere Lernleistung im Fach Physik — Jan Repenning Max-Planck-Schule Kiel 18. Ma¨rz 2011

2

Danksagung Zun¨achst m¨ochte ich mich bei all denjenigen bedanken, die diese Arbeit unterst¨ utzt und erst erm¨oglicht haben. Insbesondere bedanke ich mich bei Dr. Joke L¨ ubbecke und Jonathan Durgadoo vom IfM-Geomar in Kiel die mich bei meiner Aufgabe betreut haben und mir immer zur Seite standen, wenn es irgendwo nicht voran ging oder etwas erkl¨art werden musste. Weiterhin bedanke ich mich bei Herrn Kuckhoff von der Max-Planck-Schule Kiel, der die Betreuung dieser Arbeit von schulischer Seite u ¨bernommen hat. Dank geb¨ uhrt außerdem Dr. Joachim Dengg, der diese Kooperation zwischen Wissenschaft und Schule im Rahmen des NaT-Working-Projektes erst erm¨oglicht hat.

Jan Repenning

Inhaltsverzeichnis

3

Inhaltsverzeichnis 1. Einf¨ uhrung

5

1.1. Ozeanzirkulationsmodelle . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.2. Approximationen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2. Experimentelle Untersuchung

8

2.1. Versuchsaufbau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.2. Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

3. Mathematische Beschreibung

12

4. Numerische Simulation

13

4.1. Diskretisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.1.1. Funktionen mit einer Variable . . . . . . . . . . . . . . . . . . . 13 4.1.2. Mehrere Variablen . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.2. Gitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.3. Anfangs- und Randbedingungen . . . . . . . . . . . . . . . . . . . . . . 16 4.4. Numerischer Programmcode . . . . . . . . . . . . . . . . . . . . . . . . 18 5. Ergebnisse der Simulation

18

6. Herausforderungen & Probleme

22

6.1. Wahl der Parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 6.2. Diskretisierungsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . 23 6.3. Anfangsbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 7. Zusammenfassung & Fazit

24

A. Programmcode

26

B. Quellen

31

B.1. Literatur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

Max-Planck-Schule Kiel

4

Abbildungsverzeichnis B.2. Internetseiten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 B.3. Bildquellen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Abbildungsverzeichnis 1.

Versuchsaufbau . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.

Vertiefung der Wasseroberfl¨ache beim Einschlag des Tropfens . . . . . . 10

3.

Konzentrische Wellen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4.

Verdeutlichung der Wellen . . . . . . . . . . . . . . . . . . . . . . . . . 11

5.

Zwei Schemata zur Diskretisierung von Differenzialgleichungen . . . . . 15

6.

B-Gitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

7.

Ansicht der Welle nach 70 Sekunden . . . . . . . . . . . . . . . . . . . 19

8.

Ausbreitung der Welle . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

9.

Diagonaler Schnitt durch die Welle nach 80 Sekunden . . . . . . . . . . 21

10.

Unbrauchbare Simulation aufgrund schlecht gew¨ahlter Parameter . . . 22

11.

Einseitige Ausbreitung der Welle . . . . . . . . . . . . . . . . . . . . . 23

12.

Ungegl¨attete Welle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

Jan Repenning

5

1. Einfu ¨hrung Diese besondere Lernleistung behandelt die numerische Simulation anhand des Beispiels des Wassertropfens, der ins Wasser f¨allt. Dieses Beispiel klingt zun¨achst banal, denn welchen Nutzen k¨onnte man aus der genauen Kenntnis der Art und Weise der Ausbreitung dieser Wellen ziehen? Doch n¨aher betrachtet gibt es daf¨ ur sogar zwei Gr¨ unde. Erstens l¨asst sich mit diesem einfachen Beispiel der grunds¨atzliche Ablauf bei der Erstellung eines Modells nachvollziehen und zweitens gibt es neben den kleinen Wellen auf dem Teich auch noch andere Wellenarten, die den gleichen Gesetzen gehorchen. Besonders wichtig sind in diesem Zusammenhang Tsunamis, die sich somit auch modellieren lassen. Doch zun¨achst soll das allgemeine Vorgehen behandelt werden. Am Anfang jedes Modells steht die Beobachtung. Zuerst muss verstanden werden, wie sich Objekte verhalten, wie sie interagieren und warum dies so ist. Wenn das verstanden ist, lassen sich allgemeine Gesetzm¨aßigkeiten aufstellen, die physikalischen Gesetze. Diese sind die mathematische Repr¨asentation der Prozesse, die sich beobachten lassen. Viele dieser Gesetze sind recht einfach und man ist durchaus in der Lage, sie St¨ uck f¨ ur St¨ uck selbst zu berechnen. Doch einige Gleichungen lassen sich kaum noch leicht berechnen, stattdessen ben¨otigt man Computer, die dies tun. Dazu muss allerdings das Gesetz in eine Form gebracht werden, die vom Computer verarbeitet werden kann. An dieser Stelle kommen die numerischen Simulationen ins Spiel. Die verschiedenen numerischen Methoden erm¨oglichen es, schwierige Berechnungen mit dem Computer auszuf¨ uhren und so Aussagen u ¨ber deren Ergebnisse zu treffen. Doch warum sollte man dies tun? Diese Frage stellt sich erst einmal. Der Rechenaufwand ist groß und wirkliche Ver¨anderungen bringt es auch nicht. Doch dabei u ¨bersieht man, dass die Ergebnisse solcher Simulationen l¨angst unseren Alltag bestimmen. Man findet sie im Wetterbericht, der f¨ ur drei Tage eine sehr genaue Prognose liefert, aber auch im Auto oder in der Wirtschaft. Der Wetterbericht verarbeitet die atmosph¨arischen Daten, die gemessen werden k¨onnen und kann dann die zuk¨ unftige Entwicklung vorhersagen. In der Autoindustrie werden noch vor den Prototypen Si-

Max-Planck-Schule Kiel

6

Einf¨ uhrung

mulationen angefertigt, um beispielsweise die Aerodynamik zu optimieren und in der Wirtschaft versucht man zuk¨ unftige Trends zu erkennen und die Entwicklung vorherzusagen. Außerdem bilden sie die Grundlage f¨ ur wichtige politische Entscheidungen, so bieten sie eine Prognose f¨ ur den Klimawandel und man kann versuchen, vorherzusagen, welche Maßnahmen in Zukunft die gr¨oßte Wirkung haben werden.

Abgesehen vom Alltag spielen numerische Modelle vor allem in der Forschung, insbesondere in der Ozeanographie, eine entscheidende Rolle. Denn neben den Vorhersagen lassen sich vor allem Aussagen dar¨ uber treffen, wie sich teilweise gigantische Systeme verhalten. So ist es deutlich g¨ unstiger, Werte nur an einigen Stellen im Ozean zu messen und den Rest nachzumodellieren ( Hindcasting“), als jeden Quadratzenti” meter einzeln auszumessen. Daneben sinken die Preise f¨ ur die Rechenleistung st¨andig, was einen weiteren Preisvorteil bietet. Außerdem haben Modelle den Vorteil der Anpassbarkeit. W¨ahrend man den Ozean meist nur in Echtzeit beobachten kann, ist es bei Modellen m¨oglich, Ver¨anderungen einzubringen und ihre Auswirkungen zu beobachten, die Detailgenauigkeit nach Belieben einzustellen und zeitlich sowohl vor als auch zur¨ uck gehen kann und das in unterschiedlichen Geschwindigkeiten.

Numerische Modelle haben also vor allem in letzter Zeit erheblich an Bedeutung gewonnen und sind heute ein essentieller Bestandteil des Lebens und der Forschung. Doch neben einem allgemeinem Verst¨andnis bietet die Untersuchung von Flachwasserwellen auch ein Verst¨andnis von Tsunamis. Diese Riesenwellen sind auf dem offenen Meer noch sehr klein, besitzen aber eine große Wellenl¨ange. Sobald sie auf flacheres Wasser treffen, t¨ urmen sie sich auf und werden so zu einer großen Gefahr. Dass Tsunamis eine große Opferzahl fordern wurde bei den großen Tsunamis 2004 in Indonesien und zuletzt 2011 in Japan deutlich. Fr¨ uhwarnsysteme und Modelle, die zeigen, wie sich die Welle ausbreiten wird, k¨onnen helfen Menschen zu retten, auch wenn bei einer zu kurzen Vorlaufzeit auch das nicht hilft.

Jan Repenning

Ozeanzirkulationsmodelle

7

1.1. Ozeanzirkulationsmodelle Die Grundlage der meisten Ozeanzirkulationsmodelle bilden die Navier-StokesGleichungen. Sie beschreiben, basierend auf Massen-, Impuls- und Energieerhaltung, die hydrodynamischen Vorg¨ange. Da es sich dabei jedoch um ein recht komplexes Gleichungssystem handelt, ist es n¨otig, sofern man nur begrenzte Ressourcen zur Verf¨ ugung hat, diese allgemein g¨ ultigen Gleichungen zu vereinfachen, auch wenn dies auf Kosten der Allgemeing¨ ultigkeit geschieht. Auch die Flachwassergleichungen, die die Grundlage dieser Arbeit bilden, lassen sich aus den Navier-Stokes-Gleichungen ableiten, sind aber auch eigenst¨andig anerkannt. Grundlage f¨ ur jegliche Form der Vereinfachung bilden die Approximationen.

1.2. Approximationen ¨ So erh¨alt man unter Ausnutzung der Ublichen N¨aherungen“ die primitiven Gleichun” ¨ gen. Die Ublichen N¨aherungen“ bestehen dabei aus folgenden Approximationen: ” • Boussinesq Approximation Im Ozean sind die Dichteschwankungen meist sehr klein. Daher nimmt man die Dichte in allen Gleichungen als konstant an. Dies gilt jedoch nicht im Zusammenhang mit der Auswirkung der Schwerkraft, da sie f¨ ur gr¨oßere Ph¨anomene sorgt. So basiert beispielsweise der Golfstrom auf dem Absinken und Aufsteigen der Wassermassen an seinen Wendepunkten. • Flachwasserapproximationen 1. Sph¨arische Approximation Die Erde wird als Kugel angenommen und die Abweichungen von dieser Idealisierung ignoriert. 2. D¨ unne-Schale-Approximation Diese Approximation nutzt aus, dass die Wasserschicht (maximal 11 Kilo-

Max-Planck-Schule Kiel

8

Experimentelle Untersuchung meter dick) im Vergleich zur gesamten Erde (6 300 Kilometer Radius) sehr d¨ unn ist. 3. Hydrostatische Approximation Hierbei wird ausgenutzt, dass die horizontalen Skalen, die die Ausdehnung der Weltmeere umfassen, deutlich gr¨oßer sind als die vertikalen Skalen, also die Tiefe. ~ · ~u sondern • Traditionelle Approximation Die Corioliskraft wird nicht als 2Ω als f~ × ~u geschrieben.

2. Experimentelle Untersuchung Um die Ausbreitung der Flachwasserwellen unter kontrollierten Bedingungen zu beobachten, nutzt man ein Experiment im verkleinerten Maßstab. Die daraus folgenden Erkenntnisse lassen sich zur Formulierung allgemeiner Gesetzm¨aßigkeiten nutzen oder in diesem Fall zur sp¨ateren Verifizierung der Simulationsergebnisse.

2.1. Versuchsaufbau F¨ ur das Experiment wird ein rundes Gef¨aß mit einem Durchmesser von 23 Zentimetern genutzt. Dieses ist bis zu einer H¨ohe von 5 Zentimetern mit Wasser gef¨ ullt. Aus einer Tropfb¨ urette fallen dann gleichm¨aßig Tropfen (Volumen: ∼0,04 ml) in das Wasser. Zu Beginn ist die Wasseroberfl¨ache glatt. Das Gef¨aß wird mit einer zentral montierten Lampe von oben beleuchtet, wodurch ein R¨ uckschluss auf die optisch schwer beobachtbaren Wellen, mithilfe der von ihnen geworfenen Schatten, m¨oglich wird. Dieser Versuchsaufbau ist in Abbildung 1 dargestellt. W¨ahrend des Versuches wird die Fallh¨ohe der Tropfen variiert und der Versuch mit einer Kamera aufgezeichnet.

Jan Repenning

Ergebnisse

9

Abbildung 1: Versuchsaufbau

2.2. Ergebnisse Bei allen Fallh¨ohen zeigen sich einige charakteristische Eigenschaften. In Abbildung 2 ist eine Seitenansicht des Gef¨aßes dargestellt. Man erkennt, wie der Tropfen, hier aus einer H¨ohe von 25 Zentimetern, auf die Wasseroberfl¨ache trifft. Es kommt zu einer Verdr¨angung des Wassers und eine Eindellung bildet sich f¨ ur kurze Zeit aus. Dabei steigt die Tiefe der Delle mit zunehmender Fallh¨ohe an. Ausgehend von dieser urspr¨ unglichen Auslenkung breitet sich das Signal daraufhin in konzentrischen Kreisen aus. Abbildung 3 zeigt eine Aufsicht des Gef¨aßes nach kurzer Zeit. Darin wird das

Max-Planck-Schule Kiel

10

Experimentelle Untersuchung

Abbildung 2: Vertiefung der Wasseroberfl¨ache beim Einschlag des Tropfens konzentrische Wellenausbreitungsmuster deutlich. Außerdem zeigt sich, dass die Wellen, aus dem Zentrum propagierend, an den R¨andern des Gef¨aßes reflektiert werden. Nach einiger Zeit stellt sich wieder eine glatte Wasseroberfl¨ache ein.

Die Wellen, die von Wassertropfen aus verschiedenen Fallh¨ohen hervorgerufen werden, unterscheiden sich haupts¨achlich in ihrer Amplitude und Wellenl¨ange. An dieser Stelle soll nun keine genaue Untersuchung dieser Unterschiede folgen sondern nur eine exemplarische Berechnung anhand eines Beispiels. Da die Maximalamplituden jedoch nur schwer zu ermitteln sind, wird eine ganze Wellenfront betrachtet. Dabei wird die Gesamtl¨ange eines bestimmten Wellenzuges gemessen und durch die Anzahl der Maxima geteilt. Eine solche Untersuchung ist in Abbildung 4 dargestellt. Der schwarz markierte Bereich umfasst dabei etwa 4 Zentimeter beziehungsweise 6 Maxima und somit 5 Wellenperioden. Die Wellenl¨ange betr¨agt daher etwa

4 cm 5

= 0, 8 cm.

Jan Repenning

Ergebnisse

11

Abbildung 3: Konzentrische Wellen

Abbildung 4: Verdeutlichung der Wellen

Max-Planck-Schule Kiel

12

Mathematische Beschreibung

3. Mathematische Beschreibung Physikalische Prozesse k¨onnen durch mathematische Gleichungen beschrieben werden. Diese k¨onnen kurz und einfach sein wie das newtonsche Grundgesetz F~ = m · ~a oder sich als umfangreiche Gleichungssysteme darstellen wie es bei den Navier-StokesGleichungen der Hydrodynamik der Fall ist, die bereits in Kapitel 1.1 behandelt wurden. Hier werden die Flachwassergleichungen verwendet, die unter der Vorraussetzung gelten, dass die Wellenl¨ange der untersuchten Welle deutlich gr¨oßer als die Wasserh¨ohe ist. Dies ist nicht nur bei dem Versuch mit dem fallenden Wassertropfen der Fall, sondern gilt auch f¨ ur Tsunamis. Bei der Flachwassergleichung handelt es sich um eine partielle Differenzialgleichung zweiter Ordnung: ∂ 2f ∂ 2f ∂ 2f = gH( 2 + 2 ) ∂t2 ∂x ∂y

(1)

Diese Gleichung l¨asst sich allerdings zur einfacheren Verwendung in drei partielle Differenzialgleichungen erster Ordnung aufteilen: ∂f ∂u = −g ∂t ∂x ∂v ∂f = −g ∂t ∂y ∂f ∂u ∂v = −H( + ) ∂t ∂x ∂y

(2) (3) (4)

In diese Gleichungen fließen mehrere Variablen und Konstanten ein. g ist die Normfallbeschleunigung und betr¨agt etwa 9, 81m · s2 , H ist die Wassertiefe. f ist dann die Auslenkung der Wasseroberfl¨ache am Ort (x|y) zu der Zeit t. u bezeichnet die Geschwindigkeit in x-Richtung, v die Geschwindigkeit in y-Richtung. Die Gleichung erm¨oglicht Aussagen zur Ver¨anderung der Wasserh¨ohe im n¨achsten Zeitschritt. Dazu werden die aktuellen Auslenkungen der Umgebung des betrachteten Ortes herangezogen und so Vorraussagen gemacht. Dass die Aufteilung der Differenzialgleichung korrekt ist l¨asst sich einfach zeigen. Differenziert man partiell Gleichung

Jan Repenning

13 2 nach x und Gleichung 3 nach y und Gleichung 4 nach t, so ergibt sich: ∂ 2u ∂ 2f = −g 2 ∂t∂x ∂x

(5)

∂ 2v ∂ 2f = −g 2 ∂t∂y ∂y

(6)

∂ 2v ∂ 2u ∂ 2f + ) = −H( ∂t2 ∂t∂x ∂t∂y

(7)

Nun setzt man in die Gleichung 7 die Gleichungen 5 und 6 ein und erh¨alt: ∂ 2f ∂ 2f ∂ 2f = −H(−g − g ) ∂t2 ∂x2 ∂y 2

(8)

Klammert man nun −g aus, ergibt sich Gleichung 1. Stattdessen teilt man den untersuchten Raum und die Zeit durch ein Rechengitter in kleine Teile, die dann durch Austausch des Differentialquotienten mit dem Differenzenquotienten mit der Methode der finiten Differenzen berechnet werden k¨onnen. Diese Diskretisierung erm¨oglicht erst das numerische L¨osen des Problems.

4. Numerische Simulation 4.1. Diskretisierung 4.1.1. Funktionen mit einer Variable Bei Funktionen die nur eine Variable besitzen gibt es grunds¨atzlich drei verschiedene M¨oglichkeiten den Differentialquotienten durch einen Differenzenquotienten zu ersetzen: Vorw¨ arts Man w¨ahlt sowohl den aktuellen Punkt als auch den n¨achsten Punkt als Referenz und erh¨alt so die Ableitung.   fi+1 − fi df = dx i ∆x R¨ uckw¨ arts Statt des n¨achsten Punkts w¨ahlt man hierbei den vorherigen Punkt aus.   fi − fi−1 df = dx i ∆x

Max-Planck-Schule Kiel

14

Numerische Simulation

Zentral Bei der zentralen Differenz werden sowohl der vorherige als auch der nachfolgende Punkt betrachtet. 

df dx

 = i

fi+1 − fi−1 2∆x

4.1.2. Mehrere Variablen Im Beispiel gibt es jedoch mehrere Variablen: zwei r¨aumliche Dimensionen und eine zeitliche. Dabei ergeben sich die verschiedenen Diskretisierungsans¨atze durch Kombination der oben erl¨auterten, einfachen Prototypen. Zwei bekannte Ans¨atze sind: Leapfrog Hierbei handelt es sich um eine Zentral/Zentral-Methode. Um den n¨achsten Zeitschritt zu berechnen werden die aktuelle Umgebung und der letzte Wert genutzt. Upstream Beim Upstream-Ansatz wird der aktuelle Wert und ein Umgebungswert in einer bestimmten Richtung verwendet. Dabei haben die unterschiedlichen Schemata verschiedene Vor- und Nachteile. Das Leapfrog-Schema wirkt numerisch dispersiv, dass heißt neben den korrekten L¨osungen, gibt es noch einen oszillierenden L¨osungsanteil, der zu einer Verf¨alschung der Ergebnisse f¨ uhrt. Zudem entkoppeln sich die verschiedenen Zeitschritte leicht, was jedoch mit der Anwendung eines Zeitfilters, beispielsweise des Robert-Asselin Zeitfilters, im Rahmen gehalten werden kann. Hingegen neigt das Upstream-Schema zu einer D¨ampfung der berechneten Werte, das Schema wirkt numerisch diffusiv. Maxima verlieren an Sch¨arfe und Amplitude, da es zu einer Auff¨acherung kommt. In Abbildung 5 sind diese beiden Schemata noch einmal symbolisch dargestellt. Sie zeigt ein zweidimensionales Gitter, in welchem die i-Achse eine r¨aumliche Dimension und die n-Achse die zeitliche Dimension darstellt. Dementsprechend bezeichnet i den betrachteten Ort und i + 1 sowie i − 1 die Umgebung. Die aktuellen Daten liegen f¨ ur den Zeitpunkt n vor, n − 1 stellt daher die direkte Vergangenheit dar, w¨ahrend n + 1 f¨ ur die weitere Entwicklung steht, die berechnet werden soll. In der Zeichnung steht

Jan Repenning

Diskretisierung

15

an dieser Stelle ein weißer Knotenpunkt, wohingegen die schwarzen Knotenpunkte jeweils Daten kennzeichnen die f¨ ur die Kalkulation dieses Punktes herangezogen werden. Im Modell kommt das Leapfrog-Schema zum Einsatz, da nur dieses im Vergleich

Abbildung 5: Zwei Schemata zur Diskretisierung von Differenzialgleichungen

zum Upstream-Schema eine annehmbare Stabilit¨at der L¨osung liefert. Ausgangspunkt sind dabei die drei Differentialgleichungen erster Ordnung:

∂u ∂f = −g ∂t ∂x ∂v ∂f = −g ∂t ∂y ∂u ∂v ∂f = −H( + ) ∂t ∂x ∂y

Max-Planck-Schule Kiel

(9) (10) (11)

16

Numerische Simulation

Diskretisiert man diese nun mit dem Leapfrog-Schema erh¨alt man: ui,j,n+1 − ui,j,n−1 fi+1,j,n − fi−1,j,n = −g 2∆t 2∆x vi,j,n+1 − vi,j,n−1 fi,j+1,n − fi,j−1,n = −g 2∆t 2∆y

(12) (13)

fi,j,n+1 − fi,j,n−1 ui+1,j,n − ui−1,j,n vi,j+1,n − vi,j−1,n = −H( + ) 2∆t 2∆x 2∆y

(14)

Wobei i die Koordinate auf der x-Achse und j die Koordinate auf der y-Achse ist, w¨ahrend n den aktuellen Zeitschritt bezeichnet. Im Modell wird zus¨atzlich vorrausgesetzt, dass gilt: ∆x = ∆y.

4.2. Gitter Zur numerischen Berechnung muss die kontinuierliche L¨osung aufgeteilt werden. Da die verschiedenen Variablen in zwei Dimensionen definiert sind nutzt man f¨ ur diese Aufteilung Gitter. Bei der einfachsten Variante, dem A-Gitter, liegen alle Variablen auf demselben Gitter bei jeweils gleichen r¨aumlichen Koordinaten. Dies ist jedoch problematisch. Im Beispiel wird daher ein B-Gitter verwendet. Hierbei ist das Gitter, auf dem die Variablen f¨ ur die Geschwindigkeit in die beiden Richtungen liegen, um jeweils 21 ∆x in x- und y-Richtung relativ zum Gitter, auf dem die Wasserh¨ohe berechnet wird. Dies hat den Vorteil, dass die Aufl¨osung gr¨oßer ist, da doppelt soviele Koordinaten betrachtet werden, und auch die Gitterdispersion klein gehalten wird. Außerdem gibt es auch noch die M¨oglichkeit f¨ ur jede Variable ein eigenes versetztes Gitter zu nutzen, dabei handelt es sich dann um ein C-Gitter. Ein B-Gitter ist in Abbildung 6 gezeigt. Dabei ist das Gitter f¨ ur die Auslenkung blau gehalten, w¨ahrend das Gitter f¨ ur die Geschwindigkeiten rot ist.

4.3. Anfangs- und Randbedingungen Der simulierte Bereich hat eine r¨aumliche Begrenzung, daher muss definiert werden, was an den R¨andern der Simulation geschiet. In diesem Modell kommt es zu keinerlei

Jan Repenning

Anfangs- und Randbedingungen

17

Abbildung 6: B-Gitter Transport u ¨ber die R¨ander hinweg, ein Bereich um den Simulationsbereich herum besitzt immer weder Auslenkung noch Geschwindigkeit. Außerdem m¨ ussen die Anfangsbedingungen festgelegt werden. Da es sehr schwer ist, u ¨ber die Fallh¨ohe und die Beschaffenheit des Tropfens auf die Initialauslenkung zu schließen, m¨ ussen die empirisch gewonnenen Daten aus dem zuvor beschriebenen Versuch herangezogen werden. Die Auslenkung wird auf Grundlage dieser Beobachtungen mit einem Gauss-H¨ ugel angen¨ahert. Die Funktion dazu lautet: f (x) = −e−0,5·x

2

(15)

Beziehungsweise in zwei Dimensionen: 2 +y 2 )

f (x, y) = −e−0,5·(x

(16)

Allerdings muss diese Funktion noch normiert werden, um zu gew¨ahrleisten, dass die Anfangsauslenkung nicht den Betrag 1 hat, sondern einstellbar ist. Dabei sei K der Betrag der Auslenkung. Es gilt dann: f (x, y) = −e−0,5·(x

2 +y 2 )

·K

(17)

Abgesehen von dieser anf¨anglichen Auslenkung, ist der Rest der Wasseroberfl¨ache glatt.

Max-Planck-Schule Kiel

18

Ergebnisse der Simulation

4.4. Numerischer Programmcode Nachdem nun alle Vorraussetzungen f¨ ur eine computergest¨ utzte Berechnung gekl¨art sind, muss man die bereits vorbereiteten Gleichungen und Bedingungen in einer Programmiersprache formulieren. Speziell f¨ ur derartige Probleme, wurde das Programm MATLAB mit einer eigenen high-level Skriptsprache entwickelt. Auf diese Weise lassen sich die Gleichungen sehr direkt umsetzen. Als freie Alternative, die auch in diesem Projekt genutzt wurde, gibt es das Open- Source-Programm GNU Octave (http://www.gnu.org/software/octave/), das weitgehend kompatibel zu MATLAB-Code ist. F¨ ur die erste Differenzialgleichung (siehe Gleichung 2) ergibt sich dann auf einem B-Gitter folgender Code: dfdx(i,j)=0.5*((f(i,j,2)-f(i+1,j,2)/dx)+(f(i,j+1,2)-f(i+1,j+1,2))/dx); dudt(i,j)=-g*dfdx(i,j) - r*u(i,j,2); u(i,j,3)=dudt(i,j)*2*dt+u(i,j,1); Zun¨achst berechnet man anhand der u- und v-Werte

∂f . ∂x

Dadurch verbindet man den

¨ Schritt der Ubertragung der Werte des u/v-Gitters auf das f -Gitter und die Berechnung des Gradienten. Anschließend wird nach Gleichung 2 gerechnet um schließlich den n¨achsten Wert von u zu erhalten. Der restliche Code ist im Anhang zu finden.

5. Ergebnisse der Simulation Schon auf den ersten Blick zeigt sich, dass das Modell eine konzentrische Welle, wie sie im Experiment beobachtet werden konnte, berechnet. Im Verlauf der Zeit wird deutlich, dass auch das Modell charakteristische Eigenschaften echter Wellen widerspiegelt. So kommt es an den R¨andern des Simulationsbereiches zu einer Reflexion und sp¨ater zu Interferenz der Wellen. Außerdem verliert die Wasseroberfl¨ache in sp¨aterer Zeit an Struktur, was der glatten Oberfl¨ache in der Wirklichkeit entspricht. In der

Jan Repenning

19

(a) Dreidimensionale Seitenansicht

(b) Aufsicht

Abbildung 7: Ansicht der Welle nach 70 Sekunden

Max-Planck-Schule Kiel

20

Ergebnisse der Simulation

(a) Ausgangslage

(b) nach 60 Sekunden

(c) nach 110 Sekunden

(d) nach 170 Sekunden

(e) nach 210 Sekunden

(f) nach 350 Sekunden

Abbildung 8: Ausbreitung der Welle

Jan Repenning

21 Simulation kommt es jedoch nicht zu einer Gl¨attung der Wasseroberfl¨ache sondern zu einer gleichm¨aßigen Oszillation, da die zugrundeliegenden Gleichungen unter anderem keine Reibung ber¨ ucksichtigen. Betrachtet man einen Schnitt senkrecht zur Wassero-

Abbildung 9: Diagonaler Schnitt durch die Welle nach 80 Sekunden

berfl¨ache, l¨asst sich zus¨atzlich noch eine Aussage zu der Wellenl¨ange machen. Diese betr¨agt im Beispiel in den Außenbereichen etwa 4 Meter. Daran wird deutlich, dass es sich hierbei um eine Flachwasserwelle handelt, die eine sehr große Wellenl¨ange besitzt. Alerdings ergeben sich auch ¨ahnliche Ergebnisse bei Parametern die dem Experiment entsprechen. Es gibt also eine sehr große Diskrepanz zwischen Simulation und Experiment, wobei man anmerken muss, dass das Experiment mit den beobachteten Messwerten nicht in die Definition der Flachwasserwellen passt.

Max-Planck-Schule Kiel

22

Herausforderungen & Probleme

6. Herausforderungen & Probleme 6.1. Wahl der Parameter Das gesamte Modell reagiert ¨außerst sensibel auf die Auswahl der Parameter. Es ist schwierig eine Kombination zu finden, die eine sinnvolle Simulation erm¨oglicht und dabei stabil bleibt. Die meisten Parameter sorgen daf¨ ur, dass sich numerische Fehler nicht im weiteren Verlauf ausgleichen sondern immer weiter potenzieren. Dies sorgt dann f¨ ur gigantische Werte. Ein Beispiel f¨ ur solches Verhalten findet sich in Abbildung 10. Man erkennt deutlich, dass sich nur noch ein symmetrisches Muster aus Schwarz und Weiß ergibt, das heißt, dass die berechneten Werte so groß, beziehungsweise so klein sind, dass sie die Anzeigem¨oglichkeiten u ¨bersteigen.

Abbildung 10: Unbrauchbare Simulation aufgrund schlecht gew¨ahlter Parameter

Jan Repenning

Diskretisierungsverfahren

23

6.2. Diskretisierungsverfahren Alternativ zu zentralen Differenzen kann auch das Upstream-Schema gew¨ahlt werden. Dies hat jedoch den Nachteil, dass es oft f¨ ur einseitig propagierende Wellen sorgt. Dieser Vorgang ist in Abbildung 11 gezeigt. Die Welle propagiert haupts¨achlich nach oben rechts. Zus¨atzlich zeigt sich, dass die Werte im inneren Bereich unverh¨altnism¨aßig wachsen, w¨ahrend sie an den R¨andern wenigstens noch im Skalierungsbereich liegen.

Abbildung 11: Einseitige Ausbreitung der Welle

6.3. Anfangsbedingungen Die anf¨angliche Auslenkung der glatten Wasseroberfl¨ache muss in Form eines gaußschen H¨ ugels aufbereitet werden. Nutzt man stattdessen die M¨oglichkeit, nur ein einziges Feld des Gitters zu ver¨andern, kommt es zu einer stufenartigen Ausbreitung der Welle. Dies ist in Abbildung 12 gezeigt. Es kommt zwar auch zu einer Ausbreitung, allerdings unterscheiden sich benachbarte Felder teilweise sehr stark. Die Welle hat somit keinerlei Stetigkeit.

Max-Planck-Schule Kiel

24

Zusammenfassung & Fazit

Abbildung 12: Ungegl¨attete Welle

7. Zusammenfassung & Fazit Insgesamt zeigt dieses Projekt, dass die allgemeine Vorgehensweise zur Erstellung von Modellen gut funktioniert. Es ist m¨oglich, anhand von mathematischen Gleichungen, Prozesse, die im Experiment beobachtet werden k¨onnen, nachzuvollziehen. Allerdings gibt es daf¨ ur eine Vielzahl von M¨oglichkeiten, insbesondere in Bezug auf Diskretisierung und genutztem Gitter. Diese unterscheiden sich in ihrer Qualit¨at und m¨ ussen abh¨angig vom betrachteten Fall eingesetzt werden, um ein bestm¨ogliches Ergebnis zu gew¨ahrleisten. Leider war es in diesem Fall nicht m¨oglich eine ausreichend gute Deckung mit der Realit¨at zu erreichen, so dass eine Vergleichbarkeit der erhaltenen Werte nicht m¨oglich und das Experiment daher leider nicht ad¨aquat repr¨asentiert werden kann. Dies liegt jedoch auch in der Natur des Modells, abh¨angig von den gemachten Vereinfachungen. Ein Modell kann nat¨ urlich niemals eine vollst¨andige Deckung erm¨oglichen und das ist zum Teil auch gar nicht erw¨ unscht. Um jedoch verl¨assliche Ergebnisse zu erhalten

Jan Repenning

25 werden in der Forschung deutlich umfangreichere Modelle genutzt, um deren innere Zusammenh¨ange sich der Forscher, der nur seine Daten nutzen m¨ochte, nicht k¨ ummern muss. Das Modell wird somit immer mehr zu einer Blackbox von denen Einzelne nur noch Teile kennen. Um sie auszuf¨ uhren werden allerdings auch deutlich gr¨oßere Rechenzentren ben¨otigt. Darin zeigt sich, dass man in diesem Bereich noch lange nicht am Ende angekommen ist und Zukunftsvisionen, wie hypothetische Quantencomputer, lassen auf eine ebenso rasante Weiterentwicklung in der Zukunft hoffen. Mich pers¨onlich hat dieses Projekt auch weitergebracht. Zum einen konnte ich einen tiefen fachlichen Einblick in diese Materie nehmen, wobei mir vor allem meine beiden Betreuer geholfen haben. Zum anderen konnte ich dieses Wissen auch gleich praktisch verarbeiten und erg¨anzen. Dabei hat es mir vor allem Spaß gemacht die Formeln zu verarbeiten und die Ergebnisse der Simulation zu analysieren und zu verbessern. Dieses Projekt hat mir aber auch gezeigt, wie anf¨allig solche Modelle sein k¨onnen. Die richtige Gleichung zu nutzen wird dabei manchmal fast nebens¨achlich. Dieses Verhalten hat mich sehr u ¨berrascht, da ich davon ausgegangen war, dass sich zwar Ergebnisse von unterschiedlichen Modellen unterscheiden, aber zumindest der Großteil ein realit¨atsnahes Ergebnis liefert. Insgesamt stellen die in dieser Arbeit erworbenen F¨ahigkeiten f¨ ur mich außerdem eine gute Grundlage f¨ ur die Zukunft dar, da diese Technik, teilweise Schl¨ usselelement der heutigen Forschung ist.

Max-Planck-Schule Kiel

26

Programmcode

A. Programmcode Skript zur Berechnung der Daten und Speicherung dieser in Textdateien. Die Dateien werden in des Arbeitsverzeichnis geschrieben. Das Skript ist aus Octave mit dem Befehl wellensimulation aufzurufen. 1 % 2 % Simulation von Flachwasserwellen 3 % Jan Repenning 4 % 5 % Kommentare werden mit % eingeleitet 6 % Dieses nicht-interaktive Skript wird innerhalb der Octave7 % Kommandozeile mit "wellensimulation" aufgerufen 8 9 % Zun¨ achst werden eventuell offene Graphiken geschlossen, 10 % Variablen freigegeben und das Programm zur¨ uckgesetzt 11 12 clear all; close all 13 14 % Setzen der Parameter 15 % g: Erdbeschleunigung (auf 10 gerundet) 16 % H: Wassertiefe 17 % dx: Schrittweite im Raum (gilt auch f¨ ur y-Achse) 18 % dt: Schrittweite in der Zeit 19 % K: Anfangsauslenkung 20 % r: D¨ ampfungsfaktor 21 22 % Semikola werden ben¨ otigt, um Zuweisungen oder Befehle 23 % in einer Zeile zu trennen 24 % Ein Befehl pro Zeile ist auch ohne Semikolon legitim 25 % Das Semikolon unterdrueckt allerdings die Ergebnisausgabe 26

Jan Repenning

27

27 g=10; H=0.05; dx=1; dt=0.1; K=0.005 r=0; 28 29 % Setzen des Rechengitters 30 % Von -49 bis 49, da der Nullpunkt zentral sein soll und 31 % sonst ein 150x150 Gitter geplottet (nicht berechnet) wird 32 33 [X,Y] = meshgrid(-49:1:49, -49:1:49); 34 35 % Gauss-Huegel mit der H¨ ohe 1 36 % Setzen der Gr¨ oße der zweidimensionalen Matrix 37 38 Z = exp(-(0.5*X.^2+0.5*Y.^2)); [x,y]=size(Z); 39 40 % Erzeugung von mit Nullen gef¨ ullten Matrizen 41 u0 = zeros(x,y); v0 = zeros(x,y); 42 43 % Setzen der Anfangswerte f¨ ur die Geschwindigkeiten 44 u(:,:,1)=u0; u(:,:,2)=u0; v(:,:,1)=v0; v(:,:,2)=v0; 45 46 % Anfangswerte f¨ ur die Auslenkung 47 % Normierung der Auslenkung 48 f(:,:,1)=u0; f(:,:,2)=Z*(-K); 49 50 % Beginn der for-Schleife 51 % Festlegung der Anzahl der Zeitschritte 52 for n=1:5000 53 54

% Festlegen der Randbedingungen

55

u(1,:,:)=0; u(x,:,:)=0; u(:,1,:)=0; u(:,y,:)=0;

56

v(1,:,:)=0; v(x,:,:)=0; v(:,1,:)=0; v(:,y,:)=0;

57

u(2,:,:)=0; u((x-1),:,:)=0; v(:,2,:)=0; v(:,(y-1),:)=0;

Max-Planck-Schule Kiel

28

58

Programmcode

dfdx(1,:)=0; dfdx(x,:)=0; dfdy(:,1)=0; dfdy(:,y)=0;

59 60 61

% Beginn der beiden Schleifen, die ¨ uber jedes Feld

62

% des Gitters iterieren

63

% Eigentliche Berechnung findet hier statt

64

for i=2:(x-2)

65

for j=2:(y-2)

66 67

% Berechnung der ersten Gleichung

68

% Zun¨ achst wird dfdx anhand der Werte aus dem u/v-Gitter berechnet

69

dfdx(i,j)=0.5*((f(i,j,2)-f(i+1,j,2)/dx)+(f(i,j+1,2)-f(i+1,j+1,2))/dx);

70

% Berechnung von dudt ¨ uber dfdx

71

dudt(i,j)=-g*dfdx(i,j) - r*u(i,j,2);

72

% Berechnung des der tats¨ achlichen Geschwindigkeit

73

u(i,j,3)=dudt(i,j)*2*dt+u(i,j,1);

74 75

% Berechnung der zweiten Gleichung (Analog zu der ersten)

76

dfdy(i,j)=0.5*((f(i,j,2)-f(i,j+1,2)/dx)+(f(i+1,j,2)-f(i+1,j+1,2))/dx);

77

dvdt(i,j)=-g*dfdy(i,j) - r*v(i,j,2);

78

v(i,j,3)=dvdt(i,j)*2*dt+v(i,j,1);

79 80

% Berechnung der dritten gleichung

81

% dudx und dvdy werden ¨ uber die Werte des f-Gitters berechnet

82

dudx(i,j)=0.5*((u(i-1,j-1,2)-u(i,j-1,2)/dx)+(u(i-1,j,2)-u(i,j,2))/dx);

83

dvdy(i,j)=0.5*((v(i-1,j-1,2)-v(i-1,j,2)/dx)+(v(i,j-1,2)-v(i,j,2))/dx);

84

% Berechnung von dfdt

85

dfdt(i,j)=-H*(dudx(i,j)+dvdy(i,j));

86

% Berechnung der tats¨ achlichen Auslenkung

87

f(i,j,3)=dfdt(i,j)*2*dt+f(i,j,1);

88

Jan Repenning

29

89

end

90

end

91

% Ende der for-Schleife ¨ uber dem Gitter

92 93

% Zur¨ ucksetzen der Zeitscheibe

94

% n+1 wird zu n, n wird zu n-1

95

% Auf diese Weise kann der Code oben mit den gleichen Variablen

96

% weiterrechnen

97

u(:,:,1)=u(:,:,2); u(:,:,2)=u(:,:,3);

98

v(:,:,1)=v(:,:,2); v(:,:,2)=v(:,:,3);

99

f(:,:,1)=f(:,:,2); f(:,:,2)=f(:,:,3);

100 101

% Anwendung des Robert-Asselin Zeitfilters

102

% Dadurch wird die Entkoppelung der Zeitschritte ged¨ ampft

103

f(:,:,2)=f(:,:,2)+0.1*(f(:,:,3)-2*f(:,:,2)+f(:,:,1));

104

u(:,:,2)=u(:,:,2)+0.1*(u(:,:,3)-2*u(:,:,2)+u(:,:,1));

105

v(:,:,2)=v(:,:,2)+0.1*(v(:,:,3)-2*v(:,:,2)+v(:,:,1));

106 107

% Setzen der Variable f¨ ur den Export

108

fplot=f(:,:,2);

109 110

% Dateiname wird definiert mit f¨ unf Ziffern

111

savename=sprintf(’%05d.txt’,n);

112 113

% Speichern der Daten in einer Textdatei

114

save(savename, "fplot", "-ascii", "-tabs");

115 end 116 % Ende der for-Schleife zur zeitlichen Iteration 117 118 % Ende des Skripts

Max-Planck-Schule Kiel

30

Programmcode

Skript zur grafischen Darstellung der berechneten Daten. Der Aufruf erfolgt u ¨ber plot , die Ausgabe erfolgt im Arbeitsverzeichnis. 1 % 2 % plot.m 3 % Jan Repenning 4 % 5 % Nutzung: plot 6 % Beispiel: plot 0.005 120 7 % 8 9 % Auslesen der Kommandozeilenparameter 10 arg_list=argv (); 11 n=str2num(arg_list{2}); 12 r=str2num(arg_list{1}); 13 14 % Die Graphik wird nicht angezeigt 15 set(0, ’defaultfigurevisible’, ’off’); 16 17 % Definition des Dateinamens zum Laden der Daten 18 lade=sprintf(’%05d.txt’,n); 19 plot=load(lade, "-ascii"); 20 21 % Erstellen der dreidimensionalen Ausgabe 22 mesh(plot); 23 24 % Alternativ kann untenstehender, auskommentierter Code 25 % zur zweidimensionalen Ausgabe genutzt werden 26 % 27 %pcolor(plot); 28 %shading flat; 29 %

Jan Repenning

31

30 31 % Setzen der Skalierung ¨ uber den Parameter r 32 set(gca,’zlim’,[-r r]); 33 caxis([-r r]); 34 35 % Definition des Titels der Ausgabe 36 titel=sprintf(’%d Zeitschritte’,n); 37 title(titel) 38 39 % Gestaltung des seitlichen Farbschl¨ ussels 40 colorbar; 41 colormap ocean; 42 43 % Dateiname f¨ ur die Bildausgabe 44 % Hier sind verschiedene Formate m¨ oglich 45 filename=sprintf(’%d.png’,n); 46 47 % Ausgabe der Graphik 48 print(filename);

B. Quellen B.1. Literatur • Haidvogel, D.B. & A. Beckmann, Numerical ocean circulation modeling“, Im” perial College Press, 1999 • K¨ampf, J., Ocean Modelling for Beginners“, Springer Verlag, 2009 ” • Zus¨atzlich stammt der Großteil meiner Informationen aus Vorlesungsskripten von meinen Betreuern.

Max-Planck-Schule Kiel

32

Quellen

B.2. Internetseiten • Alle Internetseiten wurden zuletzt abgerufen am 9.3.2011. • http://www.fluvial.ch/m/2D_Flachwasser_Skript.pdf • http://www.climate.unibe.ch/~stocker/papers/stocker07skriptEKM. pdf • http://edoc.ub.uni-muenchen.de/7052/1/Erhardt_Gabriele.pdf • http://en.wikipedia.org/wiki/Shallow_water_equations • http://kiwi.atmos.colostate.edu/group/dave/pdf/ShallowWater.pdf • http://users.ices.utexas.edu/~arbogast/cam397/dawson_v2.pdf • http://kluedo.ub.uni-kl.de/volltexte/2003/1581/pdf/diss.pdf • http://www.mathworks.com/moler/exm/chapters/water.pdf • http://physics.nmt.edu/~raymond/classes/ph332/notes/shallowgov/ shallowgov.pdf • http://www.mtnmath.com/whatrh/node66.html • http://srnwp.met.hu/workshops/BadOrb_2009/Presentations/01_ Montag/01_Numerical_Developments/03_Williams/Williams.pdf

B.3. Bildquellen • Alle Graphiken in dieser Arbeit sind selbst erstellt. Nur das Titelbild enth¨alt ein fremdes Werk: http://annisleung.com/IntroEMAC/wp-content/ uploads/2009/12/water_drop.jpg (9.3.2011)

Jan Repenning