Kurs:Mathematische Modellbildung/Themen/2023-24 Wintersemester/Modellierung eines Vulkanausbruchs

Modellierung eines Vulkanausbruchs

Ätna

Modellierungsthema

Hinführung & Motivation

Jährlich brechen durchschnittlich zwischen 50 und 70 Vulkane weltweit aus. Im Jahr 2023 gibt es Stand Oktober 46 Vulkane im „continuing eruption status“ (Smithsonian Institution: Global Volcanism Program), die also momentan aktiv sind. Auch wenn sich eine Vielzahl dieser, etwa 70%, am pazifischen Feuerring, dem Rand der pazifischen Erdplatte, befinden, so bergen vereinzelt auch aktive Vulkane in unserer Nähe einige Gefahren. Wie oben schon impliziert, liegen Vulkane meistens am Rand von Erdplatten, da hier zwei solche Platten miteinander agieren. Sie falten sich entweder gegenseitig auf, tauchen untereinander ab (Subduktion) oder reiben aneinander vorbei. Diese endogenen Erdkräfte führen nun zu Erbeben, folglich Tsunamis und eben der Bildung von Vulkanen.

Vulkanausbrüche sind, genau wie alle Naturkatastrophen, nicht anthropogen kontrollierbar. Man kann sie zwar voraussagen, aber nicht verhindern. Also kann der Mensch im Hinblick auf diese Naturereignisse nicht agieren, sondern muss reagieren. Durch die Prognostizierbarkeit können jedoch schon im Vorhinein Sicherheitsmaßnamen ergriffen werden.

Vulkanausbrüche sehen oftmals spektakulär aus, denn es explodiert ein Berg, woraufhin Lava austritt und eine riesige Aschewolke in die Stratosphäre befördert wird. Auf diese soll sich in diesem Modellierungsprojekt beschränkt werden. Die Aschewolke besteht aus verschiedensten Partikeln, den Aerosolen, die schädlich für Lebewesen sind, und verbreitet sich je nach Ausbruch mehr oder weniger schnell in Höhe und Breite. Somit entsteht ein gewisser Bereich, in dem im Zeitraum nach dem Ausbruch das Leben unmöglich ist. Dieser Bereich soll durch Berechnungen mittels räumlicher Diffusion erschlossen werden, sodass man im Allgemeinen näherungsweise den Anti-Lebensraum erfassen kann. Das übergeordnete Ziel der Modellierung ist also, Berechnungen anzustellen, mit Hilfe derer man den Raum um einen Vulkan ausfindig machen kann, in dem jegliches Leben unmöglich ist nach einer Eruption, um so Menschen zu schützen.

Wir betrachten in unserem Fall exemplarisch den Ätna.

Didaktischer Hintergrund

Die ersten beiden Zyklen unseres Modellierungsprozesses bewegen sich fachmathematisch auf der Ebene der Sekundarstufe II. Bestenfalls sollte eine solche Modellierung also auch unter Anleitung in der Schule möglich sein. Fachdidaktisch stärkt ein Modellierungsprozess generell die Transferfähigkeit von SuS, deren Grundvorstellungen zu den mathematischen Begriffen und ein generelles Verständnis der Mathematik als Sprache der Wissenschaft. Der Kontext der hier geschaffen wird ist real und die Ergebnisse könnten theoretisch realen Menschen helfen. Das rechtfertigt zum einen alle mathematischen Anwendungen, die genutzt werden, und wirkt zum andern stark motivierend und interessensteigernd. Zudem ist aus der Geographiedidaktik bekannt, dass Naturkatastrophen zu den Themen gehört, die SuS am meisten interessieren. Also bietet dieser Kontext in der Schule nicht nur eine fächerübergreifende Erfahrung, sondern auch ein hochinteressantes Thema für SuS, um die mathematischen Inhalte sowie das zyklische Arbeiten der Modellbildung zu erlernen.

Gruppenteilnehmer

1. Gebhard, Sebastian
2. Leonhardt, Katrin
3. Valentin, Daniel

Zuordnung zu den Nachhaltigkeitszielen

SDG3 UN-Guidelines

• SDG 3: Good Health and Well-being

Vulkane birgen eine große Gefahr für die Menschen, die in ihrer Umgebung leben müssen. Zu den großen Gefahren für die Menschen zählen zum Einen, die Gefahr von Vulkangesteinn getroffen zu werden, zum Anderen birgt die Asche/Aesorol-Wolke eine große Gefahr für die Gesundheit der Menschen. Aus diesem Grund wird modelliert, in welcher Entfernung sich Menschen von Vulkanen aufhalten sollten, um ein gefährdungsfreies Leben führen zu können.

Niveauzuordnung

Zyklus 1/Sekundarstufe II

Im ersten Zyklus der Modellbildung wird der Ausbruch des Vulkans auf die Gefahr durch die Gesteinsbrocken begrenzt. Dabei wird mit Hilfe einer Gleichung eine Parabel beschrieben, die den Kurvenverlauf der Gesteinsbrocken darstellen soll. Das Ziel des ersten Zyklus ist es, einen ungefähren Bereich um den Vulkan zu bestimmen, in dem das Leben von Menschen nicht möglich ist, da die fallenden Gesteinsbrocken während eines Ausbruchs eine zu große Gefahr darstellen.

Zyklus 2/Sekundarstufe II

Im zweiten Zyklus soll Zyklus 1 um das Einwirken des Luftwiderstands und der Masse des Gesteins auf die Flugkurve dessen erweitert werden. Durch die Erweiterung um diese weiteren Parameter wird eine bessere Vorhersage des Bereichs möglich, in dem Menschen auf Grund zu hoher Gefahr durch Vulkangestein nicht leben können.

Zyklus 3/Universitäts-Niveau

Der dritte Zyklus der Modellierung befasst sich dann ausschließlich mit der Aschewolke, die bei einem Vulkanausbruch in die Atmosphäre ausgetoßen wird. Aufgrund dessen, dass diese Aschewolke aus vielen toxischen Aerosolen besteht, birgt sie eine erhebliche Gefahr für die Gesundheit der Menschen. In diesem Zyklus der Modellierung soll dementsprechend bestimmt werden, in welcher Entfernung Menschen zu einem Vulkan leben sollten, damit sie gesundheitlich nicht durch die Luftverschmutzung belastet werden. Zur Berechnug dieser Entfernung wird eine Diffusionsgleichung genutzt.

Modellierungszyklen

Zyklus 1 - Sekundarstufe II

Im ersten Zyklus der Modellierung betrachten wir das beim Vulkanausbruch ausgestoßene Vulkangestein und welche Flugkurve dieses einnimmt, um einen Bereich um den Vulkan zu definieren, indem ein Leben aufgrund der Gefahr von Vulkangestein getroffen zu werden zu hoch ist. Als Beispielvulkan dient der Ätna mit einer Höhe von 3357m. Dieser Teil der Modellierung orientiert sich stark an dem physikalischen Phänomen des schiefen Wurfes.

Grundlegende Voraussetzungen

• Die Vulkanwand als Graph einer linearen Funktion dargestellt wird, um später den Schnittpunkt zwischen dieser und der Flugkurve des Vulkangesteins zu bestimmen.
• Austrittswinkel $\alpha =80^{\circ }$
• Ausbruchsgeschwindigkeit $v_{0}$ ist frei wählbar im Intervall [$400{\frac {km}{h}};1000{\frac {km}{h}}$] bzw. im Intervall [$111,{\overline {1}}{\frac {m}{s}};277,{\overline {7}}{\frac {m}{s}}$]
• Ausbruchshöhe $h=3200m$
• Erdanziehungskraft $g=9,81{\frac {m}{s^{2}}}$

Bemerkung

Es ist wichtig anzumerken, dass die Masse des Vulkangesteins und die Wirkung des Luftwiderstand auf dieses vernachlässigt werden.

Herleitung der Formel

Der Vulkanausbruch kann nicht durch Gleichungen des waagerechten Wurfes beschrieben werden, da das Gestein "schräg nach oben" fliegt und nicht horizontal. Somit erfolgt die Berechnung durch Gleichungen eines schrägen Wurfes nach oben mit Anfangshöhe. Das Gestein wird beim Vulkanausbruch herausgeschleudert, erreicht seinen Scheitelpunkt, in dem Fall auch Hochpunkt, sinkt daraufhin und kommt irgendwann auf der Erdoberfläche auf.

Da in diesem Modellierungszyklus Reibungseffekte vernachlässigt werden, findet das Superpositionsprinzip seine Anwendung. In horizontaler Richtung (x-Richtung) handelt es sich um eine gleichförmige Bewegung des Gesteins mit konstanter Geschwindigkeit, wobei die Geschwinidgkeit gleich der horizontalen Komponente ${\vec {v}}_{x},_{0}$ der Anfangsgeschwindigkeit ${\vec {v}}_{0}$ ist. In vertikaler Richtung (y-Richtung) handelt es sich um eine nach oben gleichmäßig beschleunigte Bewegung des Gesteins, wobei die Startgeschwindigkeit gleich der vertikalen Komponente ${\vec {v}}_{y},_{0}$ der Anfangsgeschwindigkeit ${\vec {v}}_{0}$ ist. Die Beschleunigung ist der Erdanziehungskraft entgegengesetzt und $-g$. Beide Bewegungen behindern sich nicht gegenseitig und die Gesamtbewegung ergibt sich durch Überlagerung der beiden.

Wie bereits erwähnt spielt die Anfangsgeschwindigkeit $v_{0}$ eine Rolle, welche in ihren horizontalen und vertikalen Anteil zerlegt wird. Diese Berechnung wird anschließend durchgeführt, bei der $v_{0}=540{\frac {km}{h}}=150{\frac {m}{s}}$ frei gewählt wurde. Der horizontale Anteil kann durch die Formel "Kosinus gleich Ankathete durch Hypothenuse" bestimmt werden. Es gilt: $cos(\alpha _{0})={\frac {v_{x},_{0}}{v_{0}}}$

Stellt man diese Formel um, sodass sich der horizontale Anteil bestimmen lässt, so ergibt sich: $v_{x},_{0}=cos(\alpha _{0})\cdot v_{0}$

In unserem Fall lässt sich dieser wie folgt berechnen: $v_{x},_{0}=cos(80^{\circ })\cdot 150{\frac {m}{s}}=26,05{\frac {m}{s}}$

Der vertikale Anteil kann durch die Formel "Sinus gleich Gegenkathete durch Hypothenuse" bestimmt werden. Es gilt: $sin(\alpha _{0})={\frac {v_{y},_{0}}{v_{0}}}$

Stellt man auch diese Formel um, sodass sich der vertikale Anteil bestimmen lässt, so ergibt sich: $v_{y},_{0}=sin(\alpha _{0})\cdot v_{0}$

In unserem Fall lässt sich dieser wie folgt berechnen: $v_{y},_{0}=sin(80^{\circ })\cdot 150{\frac {m}{s}}=147,72{\frac {m}{s}}$

Es gibt vier Bewegungsgesetze, mit welchen man zu jedem Zeitpunkt $t$ die Ortskoordinaten $x$ und $y$, sowie die Geschwindigkeitskomponenten $v_{x}$ und $v_{y}$ bestimmen kann. Diese Bewegungsgesetze werden nachfolgend aufgelistet:

1. Zeit-Ort-Gesetz in x-Richtung: Gl.1

$x(t)=v_{x},_{0}\cdot t$

2. Zeit-Geschwindigkeit-Gesetz in x-Richtung: Gl.2

$v_{x}(t)=v_{x},_{0}$

3. Zeit-Ort-Gesetz in y-Richtung: Gl.3

$y(t)=-0,5\cdot g\cdot t^{2}+v_{y},_{0}\cdot t+h$

4.Zeit-Geschwindigkeit-Gesetz in y-Richtung: Gl.4

$v_{y}(t)=-g\cdot t+v_{y},_{0}$

Nachdem die Theorie vorgestellt wurde, beginnen wir mit den ersten Berechnungen. Zunächst wird die Steigzeit $t_{s}$ berechnet. Unter der Steigzeit versteht man die Zeitspanne vom Abwurf des Gesteins bis hin zum Erreichen des Scheitelpunktes. An diesem Punkt ist die y-Koordinate maximal und die y-Komponente der Geschwindigkeit 0, sodass gilt: $v_{y}(t_{s})=0$

Mit Gleichung 4 ergibt sich daraus $-g\cdot t_{s}+v_{y},_{0}=0$. Löst man diese Gleichung nach der Steigzeit auf, so ergibt sich: $t_{s}={\frac {v_{y},_{0}}{g}}$

In unserem Beispiel berechnet sich die Steigzeit zu: $t_{s}={\frac {147,72{\frac {m}{s}}}{9,81{\frac {m}{s^{2}}}}}=15,1s$

Mit Hilfe der Gleichungen 1 und 3 ergibt sich für die Koordinaten des Scheitelpunktes folgendes: $S({\frac {v_{x},_{0}\cdot v_{y},_{0}}{g}}{\bigg \vert }{\frac {(v_{y},_{0})^{2}}{2\cdot g}}+h)$

In unserem Beispiel berechnet sich der Scheitelpunkt zu:

$S({\frac {26,05{\frac {m}{s}}\cdot 147,72{\frac {m}{s}}}{9,81{\frac {m}{s^{2}}}}}{\bigg \vert }{\frac {(147,72{\frac {m}{s}})^{2}}{2\cdot 9,81{\frac {m}{s^{2}}}}}+3200m)=S(392,26m{\bigg \vert }4312,19m)$

Nun wird die Wurfzeit $t_{w}$ berechnet. Dies ist die Zeitspanne vom Abwurf des Gesteins bis zum Auftreffen auf der Erdoberfläche. An diesem Punkt ist die y-Koordinate 0, sodass $y(t_{w})=0$ Mit Gleichung 3 ergibt sich daraus: $-0,5\cdot g\cdot (t_{w})^{2}+v_{y},_{0}\cdot t_{w}+h=0$

Stellt man diese Gleichung nun um, sodass man $t_{w}$ bestimmen kann, so ergibt sich dafür folgendes: $t_{w}={\frac {v_{y},_{0}+{\sqrt {(v_{y},_{0})^{2}+2\cdot g\cdot h}}}{g}}$

In unserem Beispiel ergibt sich als Wurfzeit: $t_{w}={\frac {147,72{\frac {m}{s}}+{\sqrt {(147,72{\frac {m}{s}})^{2}+2\cdot 9,81{\frac {m}{s^{2}}}\cdot 3200m}}}{9,81{\frac {m}{s^{2}}}}}=44,71s$

Anschließend wird noch die Wurfweite $w$ berechnet. Dies ist die Entfernung vom Koordinatenursprung, bei dem $x=0$, bis zum Auftreffen auf der Erdoberfläche.

Da die Wurfweite die Strecke ist, welche das Gestein während der Wurfzeit in horizontaler Lage zurücklegt, gilt: $w=x(t_{w})$.

Daraus ergibt sich mit Gleichung 1: $w=v_{x},_{0}\cdot t_{w}$

In unserem Beispiel ergibt sich als Wurfweite: $w=26,05{\frac {m}{s}}\cdot 44,71s=1164,54$

Die Gleichung der Bahnkurve ergibt sich aus den Gleichungen 1 und 3, bei denen man die Zeit weglässt: $y(x)=-0,5\cdot {\frac {g}{(v_{x},_{0})^{2}}}\cdot x^{2}+{\frac {v_{y},_{0}}{v_{x},_{0}}}\cdot x+h$

Mit Hilfe dieser Gleichung lässt sich zu jeder x-Koordinate die zugehörige y-Koordinate berechnen.

In unserem Beispiel lässt sich diese Gleichung wie folgt bestimmen: $y(x)=-0,5\cdot {\frac {9,81{\frac {m}{s^{2}}}}{(26,05{\frac {m}{s}})^{2}}}\cdot x^{2}+{\frac {147,72{\frac {m}{s}}}{26,05{\frac {m}{s}}}}\cdot x+3200m=-0,007{\frac {1}{m}}\cdot x^{2}+5,67\cdot x+3200m$

Ergebnis

Animation Zyklus 1

In nebenstehendem Video erkennt man die Visualisierung unserer Ergebnisse aus dem Zyklus 1. Die Anfangs allgemein aufgestellten Formeln lassen sich durch Schieberegler verwirklichen, die mithilfe einer Animation variiert werden. Der erste Schieberegler beschreibt die Austrittsgeschwindigkeit $v$ des Steins und der zweite den Austrittswinkel $\alpha$. Somit lassen sich in einem Video nicht nur das Ergebnis unserer Beispielwerte visualisieren, sondern eine Variation von verschiedenen möglichen Vulkanausbrüchen.

Zyklus 2 - Sekundarstufe II

Im Anschluss an die Berechnungen aus Zyklus 1 erweitern wir diese im zweiten Zyklus um die Wirkung des Luftwiderstands und der Masse des Gesteins auf dessen Flugkurve.

Was ist der Luftwiderstand?

Die Richtung der Luftwiederstandskraft ist der Bewegungsrichtung des Vulkangesteins entgegengesetzt. Da sich die Bewegungsrichtung des Gesteins durchgehend ändert, ändert sich auch die Richtung der Luftwiderstandskraft durchgängig. Mithilfe der Newton´schen Luftwiderstandsformel lässt sich der Betrag der Luftwiderstandskraft $F_{L}$ beschreiben.

Die Formel lautet: $F_{L}={\frac {1}{2}}\cdot c_{w}\cdot \rho \cdot A\cdot v^{2}$

wobei $c_{w}$ der Luftwiderstandswert eines Gesteinsbrockens, $\rho$ die Dichte der Luft, $A$ die Querschnittsfläche eines Gesteinsbrockens und $v$ der Betrag der Geschwindigkeit (= Tempo) ist.

Differentialgleichung der Gesteinsflugbahn mit Luftwiderstand

Graphik zur Veranschaulichung der Kräftewirkung

Herleitung der Formeln zur Berechnung der einzelnen Geschwindigkeiten

Rechts sieht man eine Abbildung mit der Kraft des Luftwidrstandes $F_{L}$ und der Erdanziehungskraft $F_{G}$, welche beide während des Fluges des Vulkangesteins auf diesen einwirken. Für die Koordinatenachsen ergeben sich folgende Gleichungen:

$m\cdot a_{x}=-F_{Lx}$

$m\cdot a_{y}=-F_{Ly}-F_{G}$bzw.

$m\cdot {dv_{x} \over dt}=-F_{L}\cdot cos(\alpha )$

$m\cdot {dv_{y} \over dt}=-F_{L}\cdot sin(\alpha )-F_{G}$

wobei gilt

$sin(\alpha )={\frac {v_{y}}{v}}$

$cos(\alpha )={\frac {v_{x}}{v}}$

Die Luftwiderstandskraft lässt sich wie folgt umschreiben:

$F_{L}=K\cdot v^{2}$

K ist Faktor des Luftwiderstandes und berechnet sich mit $K={c_{w}\cdot \rho \cdot A \over 2}$ ,

wobei $c_{w}$ der Luftwiderstandswert eines Gesteinsbrockens, $\rho$ die Dichte der Luft und $A$ die Querschnittsfläche eines Gesteinsbrockens ist.

Da Vulkangestein vielfältige Formen annehmen kann betrachten wir hier einen kugelförmigen Gesteinsbrocken, für den wir den Radius $r$ der Querschnittsfläche frei gewählt haben.

Mit den folgenden Werten lässt sich später der Luftwiderstandsfaktor $K$ berechnen. $r=0,1m$, $c_{w}=0,45$ , $\rho =1,2{\frac {kg}{m^{3}}}$ , $A=0,03m^{2}$

$K={\frac {0,45\cdot 1,2{kg \over m^{3}}\cdot 0,03m^{2}}{2}}=8,1\cdot 10^{-3}{\frac {kg}{m}}$

Damit gilt dann

$m\cdot {dv_{x} \over dt}=-K\cdot v^{2}\cdot {v_{x} \over v}=-K\cdot v\cdot v_{x}$ (1)

$m\cdot {dv_{y} \over dt}=-K\cdot v^{2}\cdot {v_{y} \over v}-m\cdot g=-K\cdot v\cdot v_{y}-m\cdot g$ (2)

Zunächst wird (1) betrachtet, wobei die Geschwindigkeit abhängig von der Zeit t in Sekunden betrachtet wird, es ergibt sich dann:

$m\cdot {d \over dt}\cdot v_{x}(t)=-K\cdot v(t)\cdot v_{x}(t)$ wobei $v(t)={\sqrt {(}}v_{x}(t)^{2}+v_{y}(t)^{2})$

also $m\cdot {d \over dt}\cdot v_{x}(t)=-K\cdot {\sqrt {(}}v_{x}(t)^{2}+v_{y}(t)^{2})\cdot v_{x}(t)$

Nach einem Zeitschritt $t_{n}$zu $t_{n+1}$den wir als $\vartriangle t$ bezeichnen, beträgt die Geschwindigkeit des Gesteinsbrockens in x-Richtung:

$v_{x}(t_{n+1})={(-K\cdot {\sqrt {(}}v_{x}(t_{n})^{2}+v_{y}(t_{n})^{2})\cdot v_{x}(t_{n})\cdot \vartriangle t) \over m}+v_{x}(t_{n})$

Durch einestzen in Gleichung (2) erhalten wir:

$m\cdot {d \over dt}\cdot v_{y}(t)=-K\cdot v(t)\cdot v_{y}(t)-m\cdot g$verwenden wir ebenso, dass $v(t)={\sqrt {(}}v_{x}(t)^{2}+v_{y}(t)^{2})$

Nach einem Zeitschritt $t_{n}$zu $t_{n+1}$den wir als $\vartriangle t$ bezeichnen, beträgt die Geschwindigkeit des Gesteinsbrockens in y-Richtung:

$v_{y}(t_{n+1})={(-K\cdot {\sqrt {(}}v_{x}(t_{n})^{2}+v_{y}(t_{n})^{2})\cdot v_{y}(t_{n})-m\cdot g)\cdot \vartriangle t \over m}+v_{y}(t_{n})$

Zusammengefasst wird hier nun die Geschwindigkeit in x- und y-Richtung nach jedem Zeitschritt $\vartriangle t$ berechnet, indem zu der Geschwindigkeitsänderung, die durch die in y-Richtung wirkende Erdanziehung und den Luftwiderstand zustande kommt, die Geschwindigkeit vor dem jeweiligen Zeitschritt hinzu addiert wird.

Die Gesamtgeschwindigkeit nach $\vartriangle t$ berechnet sich dann mit:

$v(t_{n+1})={\sqrt {(}}v_{x}(t_{n+1})^{2}+v_{y}(t_{n+1})^{2})$

In unserem Modell wählen wir $\vartriangle t$= 0,1s (Sekunden).

Nun können wir nach jedem Intervall $\vartriangle t$ die neue (Gesamt-)Geschwindigkeit bestimmen, nachdem wir die Geschwindigkeiten in x- und y-Richtung ermittelt haben.

Berechnung der zurückgelegten Strecken und die Zeit bis zum Aufprall

Mit Hilfe der oben aufgestellten Formeln lassen sich die vom Gesteinsbrocken zurückgelegten Strecken in x- und y-Richtung für jedes Zeitintervall $\vartriangle t$ berechnen.

Um nun die zurückgelegten Strecken in x- und y-Richtung zu bestimmen, wird die Gleichung $v={s \over t}$ ,wobei $v$ die Geschwindigkeit, $s$ die zurückgelegte Strecke und $t$ die Zeit angibt, so umgestellt, dass man für die zurückgelegten Strecken in x- und y-Richtung folgende Gleichungen erhält.

Die zurückgelegte Strecke des Gesteinsbrockens in x-Richtung lässt sich mit folgender Gleichung berechnen:

$s_{x}(t_{n+1})=s_{x}(t_{n})+v_{x}(t_{n})\cdot \vartriangle t$

und die zurückgelete Strecke in y-Richtung mit folgender Gleichung:

$s_{y}(t_{n+1})=s_{y}(t_{n})+v_{y}(t_{n})\cdot \vartriangle t$

Auszug aus den zur Berechnung benötigten Befehlen 1

Auszug aus den zur Berechnung benötigten Befehlen 2

Im folgenden ist ein Auzug aus unseren Berechnungen, ebenso wie ein Vergleich der beiden Graphen. Diese Graphen zeigen die Flugkurven des Gesteinsbrockens unter Berücksuchtigung des Luftwiderstands und unter Vernachlässigung des Luftwiderstands.

Auszug aus den berechneten Werten

Graphen der Flugbahn

Um zu bestimmen wann der Gesteinsbrocken bei Berücksichtigung des Luftwiderstandes auf dem Boden auftrifft muss das Interval $[t_{i},t_{j}]$ betrachtet werden, indem die berechneten Werte für die zurückgelegte Strecke in y-Richtung ihr Vorzeichen wechseln. Diesen Vorzeichenwechsel sieht man rechts in der Abbildung.

Veranschaulichung Vorzeichenwechsel

Um nun im weiteren Verlauf das genaue Zeitintervall $[t_{i},t_{j}]$ zu bestimmen, indem der Gesteinsbrocken auf der Erdoberfläche auftritt benötigt man die Intervallschachtelung. Dies wurde hier allerdings nicht betrachtet, da ein Zeitschritt von $\vartriangle t=0,1s$ in unserem Fall ausreichend genau ist.

Ergebnis

Die oben angeführten Tabellen spiegeln die von uns berechneten Werte wider. Zudem visualisiert die Abbildung die Flugbahn unter Beachtung des Luftwiderstands im Vergleich zur Flugbahn unter Vernachlässigung des Luftwiderstands. Auf Grundlage unserer Beispielwerte würde ein Gesteinsbrocken mit einer Masse von $12,6kg$ eine Strecke von $730m$ zurücklegen.

Zyklus 3 - Universitätsniveau

Der dritte Zyklus unserer Modellierung behandelt jetzt nicht mehr das beim Vulkanausbruch ausgestoßene Gestein, sondern er fokussiert sich auf die Aerosol/Aschewolke, die bei einem solchen Ausbruch in die Atmosphäre ausgestoßen wird.

Ziel

Das Ziel dieses Zyklusses ist es, eine Verbreitung der Aschewolke zu simulieren um so näherungsweise ein Gebiet um den Vulkan zu bestimmen, in dem die Luftbelastung durch Asche und Aerosole zum Leben zu hoch ist.

Wichtige Zeichen

-Divergenz: $div$

-Gradient: $\nabla$

-Laplace-Operator: $\Delta$

Mathematische Darstellung

Der Begriff Diffusion beschreibt die zufällige Bewegung von Molekülen eines Stoffes in Flüssigkeiten oder Gasen. Diese von Robert Brown entdeckten Bewegungen werden somit auch als Brownsche Bewegung bezeichnet. Hierbei streben die Teilchen, wie so oft in der Natur, einem Gleichgewicht der Konzentrationen im Raum (Flüssigkeit oder Gas) entgegen. Man spricht also von räumlicher Diffusion.

Das Gleichgewicht des Systems wird durch eine partielle Differenzialgleichung beschrieben, die wie folgt aussieht:

$$div(c(x)\nabla u(x))=0$$

Hierbei sind die beiden Funktionen $u$ und $c$ auf ${\mathbb {R} }^{n}\to {\mathbb {R} }$ definiert und beschreiben die unbekannte Konzentration ($u$) und den ortsabhängigen Diffusionskoeffizienten ($c$). Wenn gilt $c(x)=1$, dann ergibt sich hierfür die Laplace Gleichung $\Delta u(x)=0$.

Um die verschiedenen Konzentrationen zu verschiedenen Zeitpunkten darstellen zu können, muss eine zusätzliche zeitliche Komponente in die Gleichung eingebracht werden. Die Funktion $u=u(x,t)$ stellt genau diese Komponente dar. Zusammen mit obiger Gleichung ergibt sich die Diffusionsgleichung

$$\partial _{t}u(x,t)-div[c(x)\nabla u(x,t)]=0.$$

Da in unserem Modell der Vulkan die Stelle mit der hohen Konzentration zum Zeitpunkt 0 darstellt und keine externe Konzentrationsquelle existiert kann in dieser Modellierung mit einer solchen instationären Gleichung gearbeitet werden.

Da es in der Natur nur sehr selten windstill ist, muss mit einer Hintergrundgeschwindigkeit gerechnet werden. Diese wird durch den Vektor ${\vec {v}}$ modelliert und man spricht nun von einer Konvektions-Diffusionsgleichung.

$$\partial _{t}u(x,t)-div[c(x)\nabla u(x,t)]+div[u(x,t){\vec {v}}(x,t)]=0.$$

In unserem Modell wird sowohl beim Diffusionskoeffizienten $c$, als auch bei der Hintergrundgeschwindigkeit ${\vec {v}}$ Konstanz angenommen, wodurch man die Gleichung folgendermaßen vereinfachen kann:

$$\partial _{t}u(x,t)-c\Delta u(x,t)+{\vec {v}}\cdot \nabla u(x,t)=0.$$

Physikalische Herleitung

Um das Phänomen der Diffusion von Orten mit höherer Konzentration hin zu Orten mit niedriger Konzentration mathematisieren zu können, wird ein dünner, horizontal ausgerichteter Stab mit einer vorgegebenen Länge $L$ betrachtet, in dem die Verteilung der Teilchen in vertikaler Richtung homogen, also gleichmäßig ist.

Teilchendiffusion

Man untersucht also in diesem Intervall $[0,L]$ die Konzentration an Molekülen, indem ein bestimmter Punkt $x\in [0,L]$ betrachtet wird, den die Teilchen von links nach rechts überqueren. Dies geschiet in einer Diffusionsflussrate, $q(x,t)$. Außerdem wird die Gesamtkonzentration im blauen Abschnitt (begrenzt durch $x$ und $x+\Delta x$) zum Zeitpunkt $t$ als Integral $\int _{x}^{x+\Delta x}u(s,t)ds$

dargestellt.

Wichtige Physikalische Gesetze

1. Erhaltungsgesetz: Der Stoff in einem Volumen entsteht oder verschwindet nicht, sondern bleibt erhalten (sofern nicht von außen zugefügt oder weggenommen).
2. Fick'sches Gesetz: Die Diffusionsflussrate (Teilchenstromdichte) $q(x,t)$ ist proportional zum Konzentrationsgradient, die Proportionalitätskonstante ist der negative Diffusionskoeffizient -$c(x),$
   $$q(x,t)=-c(x)\partial _{x}u(x,t).$$
   Dieses Gesetz drückt die Tatsache aus, dass sich die Teilchen in der Richtung des Konzentrationsabstiegs bewegen.

(--> Räumliche Diffusion)

Da von außen keine Moleküle hinzugefügt werden, ist nach dem Erhaltungsgesetz davon auszugehen, dass die Konzentrationsveränderungen in $[x,x+\Delta x]$ durch Teilchenbewegungen zustande kommt.

$${\frac {\partial }{\partial t}}\int _{x}^{x+\Delta x}u(s,t)ds=q(x,t)-q(x+\Delta x,t)$$
 

Teilt man nun diese Gleichung durch $\Delta x$ und betrachtet den Grenzübergang $\Delta x\to 0$, erhält man die Kontinuitätsgleichung:

$$\partial _{t}u(x,t)=\partial _{x}q(x,t)$$

Mithilfe des Fick´schen Gesetzes kann man nun die Diffusionsrate $q(x,t)$ durch $-c(x)\partial _{x}u(x,t)$ ersetzen und erhält somit die homogene Diffusionsgleichung

$$\partial _{t}u(x,t)-\partial _{x}[c(x)\partial _{x}u(x,t)]=0.$$

Da bei unserem Vulkanausbruch von einem konstanten Diffusionskoeffizienten ausgegangen wird, kann man die Gleichung auf folgende vereinfachen:

$$\partial _{t}u(x,t)-c\partial _{xx}u(x,t)=0$$

An dieser Stelle beachte man nun den Wind als Hintergrundgeschwindigkeit ${\vec {v}}$, der zusätzlich zum Materialtransport beiträgt, wodurch sich der gesamte Materialfluss $j$ aus Diffusions- und Konvektionsfluss zusammensetzt.

$$j=-c\partial _{x}u+{\vec {v}}u.$$

Mithilfe der oben aus dem Erhaltungsgesetz hergeleiteten Gleichung ergibt sich:

$${\frac {\partial }{\partial t}}\int _{x}^{x+\Delta x}u(s,t)ds=j(x,t)-j(x+\Delta x,t)$$

Auch hier wird durch $\Delta x$ geteilt und auf den Grenzübergang $\Delta x\to 0$ geschaut und somit erhält man die Konvektions-Diffusionsgleichung:

$$\partial _{t}u(x,t)-c\partial _{xx}u(x,t)+\partial _{x}[{\vec {v}}(x,t)u(x,t)]=0.$$

Durchführung mithilfe von Comsol

Annahmen und Grundlagen

Herleitung des Diffusionskoeffizienten

Bei Vulkanausbrüchen werden in kurzer Zeit große Mengen an Gasen und Partikel weit in die Atmosphäre geschleudert. Da beim Vulkanausbruch nicht nur ein Element ausgestoßen wird, haben wir durch Recherche eine repräsentative chemische Zusammensetzung gefunden. Die Chemie des Vulkans bestimmt seine Gefahr. Die Zusammensetzung vulkanischer Gase ist abhängig von der chemischen Zusammensetzung der Lava und von Vulkan zu Vulkan unterschiedlich. In unserem Beispiel mit dem Ätna setzen sich die Vulkangase wie folgt zusammen: Fast die Hälfte der Vulkangase, die dabei frei werden, ist Wasserdampf. Zudem folgt ein Viertel Kohlenstoffdioxid und Schwefeldioxid. Einen relativ kleinen Anteil machen Wasserstoff und Kohlenstoffmonoxid aus. Anschließend haben wir mit Hilfe der einzelnen Diffusionskoeffizienten der Ionen, aus denen die Vulkangase bestehen, den Diffusionskoeffizient der Asche bestimmt. Da Wasserstoff und Kohlenstoffmonoxid zusammen gerade einmal 1% ausmachen, werden diese hier bei der Rechnung vernachlässigt. Schwefeldioxid hat einen einen Diffusionskoeffizienten von ${\frac {1}{10^{9}}}{\frac {m^{2}}{s}}$. Wasserdampf und Kohlenstoffdioxid bilden zusammen Kohlensäure, welche aus zwei Protonen und einem zweifach negativ geladenem Carbonat-Ion besteht. Der Diffusionskoeffizient eines Protons beträgt $9,31\cdot {\frac {1}{10^{9}}}{\frac {m^{2}}{s}}$ und der eines Carbonat-Ions $0,955\cdot {\frac {1}{10^{9}}}{\frac {m^{2}}{s}}$. Die Partikel werden als Asche bezeichnet und besteht hauptsächlich aus Siliciumdioxid. Davon liegt der Diffusionskoeffizient bei ${\frac {1}{10^{12}}}{\frac {m^{2}}{s}}$. Rechnet man alle einzelnen Diffusionskoeffizienten zusammen, so ergibt sich folgendes Ergebnis: $D_{c}=2\cdot 9,31\cdot {\frac {1}{10^{9}}}{\frac {m^{2}}{s}}+0,955\cdot {\frac {1}{10^{9}}}{\frac {m^{2}}{s}}+{\frac {1}{10^{12}}}{\frac {m^{2}}{s}}=2,1\cdot {\frac {1}{10^{8}}}$. Dieses Ergebnis dividiert man anschließend durch 4, da sich der Diffusionskoeffizient aus vier chemischen Verbindungen zusammensetzt, sodass der gesuchte Diffusionskoeffizient $D_{c}=5\cdot {\frac {1}{10^{9}}}{m^{2}}{s}$ beträgt.

Windrichtung

Da es in der Realität selten Windstill ist, haben wir bei unserer Modellierung eine Hintergrundgeschwindigkeit mit einberechnet. Aufgrund der Lage vom Ätna im Gradnetz kommt diese von Südwesten.

Durchführung

Erstellen einer Geometrie

Um in Comsol den physikalischen Vorgang der räumlichen Diffusion simulieren zu können, muss zunächst ein Raum eingegrenzt werden, in dem der Vorgang geschiet. Hierzu wird ein Rechteck erstellt, dass den um den Vulkan liegenden Raum beschreibt und ein darin liegender Kreis stellvertretend für den Vulkan, der als Konzentrationsquelle dient. Das Rechteck hat eine Breite von $15km$ und eine Höhe von $10km$. Also umspannt unser Grundraum $150km^{2}$ um den Vulkan herum. Der für den Krater stehende Kreis hat einen Radius von $0,5km$, also eine Fläche von ca. $0,79km^{2}$.

Einstellen der verschiedenen Parameter

Eigenschaften von Luft

Zunächst muss festgelegt werden, in welchem Stoff die Diffusion erfolgt, hierbei wählen wir natürlich Luft (air), welche schon in Comsol implementiert ist mit ihren Eigenschaften.

Man geht beim Transport verdünnter Spezies von idealer Speziesaktivität, der Konvektion als zusätzlichem Transportmechanismus und einer linearen Konzentration aus.

Als Diffusionskoeffizienten wählen wir $D_{c}=5\cdot {\frac {1}{10^{9}}}{m^{2}}{s}$ (--> siehe Erklärung oben).

Die Reaktionsgeschwindigkeit wird auf $100mol/(m^{3}*s)$ festgelegt.

Beim Netz haben wir uns für das dritt-feinste physikgesteurte Netz entschieden, das "feiner" genannt wird.

Ergebnisse Zyklus 3

Mithilfe obiger Einstellungen kann der PC nun einen zeitabhängigen Verlauf der Ausbreitung unserer Vulkanasche in unserem Gebiet simulieren. Hier sehen wir das Ergebnis:

Verschiebung der Aschewolke unter Windeinfluss aus Südwesten

Verbreitung der Aschewolke in Isolinien

In der folgenden Abbildung sieht man die Ausbreitung der Aschewolke und ihre Konzentration ca. 15 min nach dem Ausbruch.

Ausbreitung der Aschewolke

Ergebnisse

Die Ergebnisse der Zyklen 1 und 2 zeigen stark, dass der Luftwiderstand und die Masse des Gesteins eine große Rolle spielen, denn es wird deutlich, dass die Flugbahn eines Gesteinsbrockens durch die Einwirkung des Luftwiderstandes deutlich verkürzt ausfällt (S.Abb: Zyklus 2). Auch in Zahlen ausgedrückt wird dieser Unterschied gut erkennbar, denn während ein Gesteinsbrocken ohne Beachtung des Luftwiderstands und der Masse über $1100m$ weit fliegt, so fliegt er bei einer doppelt so hohen Austrittsgeschwindigkeit und gleichbleibendem Austrittswinkelwinkel unter Beachtung der Gesteinsmasse ($12,6kg$) und des Luftwiderstandes nur noch etwas über $730m$ weit.

Zyklus 3, der sich mit der Aerosol/Aschewolke beschäftigt hat, hat sehr anschaulich gezeigt, dass sich diese Wolke sehr rasch und sehr weit bewegen kann. Im Laufe der Modellierung hat sich gezeigt, dass es schwierig ist einen genauen Bereich festzulegen, in dem Leben aufgrund einer zu hohen Aerosolbelastung nicht möglich ist. Allerdings ist in obenstehendem Ergebnis die Ausbreitung der Aschewolke nach 15 Minuten veranschaulicht. Hiebei ist zu erkennen, dass die Ausbreitung den Rand unseres festgelegten Gebietes erreicht hat. Also ist davon auszugehen, dass die Wolke, je nach Windrichtung, jeden Punkt in unserem Gebiet von $150km^{2}$ innerhalb von 15 bis 20 Minuten erreichen kann. Der Bereich mit der höchsten Konzentration (rote Farbe im Schaubild) bleibt jedoch relativ nahe am Vulkan.

Generell sind alle Angaben abhängig von verschiedenen Natureinflüssen, wie Windstärke oder Windrichtung, aber auch von der schwere eines Vulkanausbruchs. Je schwerer dieser, desto größer sind die Gefahrenbereiche. Unsere verwendeten Werte basieren auf der Annahme eines mittelschweren Vulkanausbruchs.

Literatur

https://www.aqion.de/site/diffusionskoeffizienten

https://www.eskp.de/grundlagen/naturgefahren/zusammensetzung-vulkanischer-gase-935614/#:~:text=Das%20h%C3%A4ufigste%20vulkanische%20Gas%20ist%20Wasserdampf%20oder%20H,der%20Gase%20variieren%20kann%20%28siehe%20dazu%20obige%20Abbildung%29