Die Level Set Methode für das

Stefan Problem

CES Master Seminararbeit

Autor:

Marco SchoosRWTH Aachen

marco.schoos@rwth-aachen.deMatrikelnr.: 312553

Betreuerin:

Dr. Julia KowalskiAICES

Aachen, 2016

INHALTSVERZEICHNIS

Inhaltsverzeichnis

1 Einleitung 5

2 Das Stefan Problem 5

3 Die Level Set Methode 6

3.1 Idee der Level Set Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63.2 Update der Level Set Funktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73.3 (Re)Initialisierung der signed distance Funktion . . . . . . . . . . . . . . . . . . . . . . . . . . 83.4 Update des Temperaturfeldes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

4 Diskretisierung 9

4.1 Diskretisierung des Updates der Level Set Funktion . . . . . . . . . . . . . . . . . . . . . . . . 94.1.1 Approximation des Temperatursprungs . . . . . . . . . . . . . . . . . . . . . . . . . . 94.1.2 Stetige Erweiterung des Temperatursprungs . . . . . . . . . . . . . . . . . . . . . . . . 104.1.3 Berechnung der Geschwindigkeitsfunktion F . . . . . . . . . . . . . . . . . . . . . . . . 104.1.4 Diskretisierung der Advektionsgleichung . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4.2 Diskretisierung der (Re)Initialisierung der signed distance Funktion . . . . . . . . . . . . . . . 114.3 Diskretisierung der Temperaturgleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124.4 Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

5 Ergebnisse 14

2

ABBILDUNGSVERZEICHNIS

Abbildungsverzeichnis

1 Darstellung des zweiphasigen Gebietes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Beispielplot einer signed distance Funktion mit Contourplot . . . . . . . . . . . . . . . . . . . 73 Darstellung der vier Koordinatenrichtungen x, y, η und ζ . . . . . . . . . . . . . . . . . . . . 74 Beispiel für ein Feld des Temperatursprungs in einer Dimension . . . . . . . . . . . . . . . . . 95 Beispiel für eine stetige Erweiterung eines Feldes des Temperatursprungs in einer Dimension . 106 Beispiel für die Konstruktion von PL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 Analytische Lösung für das Stefan Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3

Abkürzungsverzeichnis

Lateinische Symbole

Symbol Beschreibung Einheit

c Wärmekapazität J K−1

cfl CFL-Zahl −D Vereinigung der Gebiete der üssigen und festen Phase −d euklidischer Abstand mF Geschwindigkeitsfunktion −G Godunov Funktion −h Gitterweite mk Wärmeleitfähigkeit W m−1 K−1

L Latente Erstarrungswärme Jn Normalvektor der Phasengrenze mS Vorzeichenfunktion −T Temperatur Kt Zeit su Sprung im Wärmeuss −V Geschwindigkeit der Phasengrenze in Normalenrichtung m s−1

x Ortskoordinate my Ortskoordinate m

Griechische Symbole

Symbol Beschreibung Einheit

∆t Zeitdierenz sε Hilfsparameter −η Ortskoordinate im gedrehten Koordinatensystem mΓ Phasengrenze −Ω Gebiet der festen Phase −Ωc Gebiet der üssigen Phase −φ Level Set Funktion −ρ Dichte kg m−3

ζ Ortskoordinate im gedrehten Koordinatensystem m

Indizes

Index Beschreibung

e Eisi Zählvariablej Zählvariablel liquid (üssig)m melt (Schmelz-)s solid (fest)w Wasser

Mathematische Symbole

Symbol Beschreibung

∆ Laplace OperatorR Menge der reellen Zahlen∇ Nabla Operator

2 DAS STEFAN PROBLEM

1 Einleitung

In dieser Arbeit wird eine numerische Methode vorgestellt um das Stefan Problem zu lösen. Betrachtet wirddazu ein zweidimensionales Gebiet und die Level Set Methode ist ein Teil der Diskretisierungsstrategie. Imzweiten Kapitel wird zunächst das Stefan Problem selbst kurz erläutert. Im dritten Kapitel wird die LevelSet Methode beschrieben, während im vierten Kapitel auf die zugehörige Diskretisierung eingegangen wird.Das fünfte und abschlieÿende Kapitel geht kurz auf verschiedene Berechnungsergebnisse ein.Der Groÿteil dieser Arbeit wurde aus [1] entnommen, weitere Verwendung aus der Literatur ist an derentsprechenden Stelle referenziert.

2 Das Stefan Problem

Das Stefan Problem fragt nach der Temperaturverteilung in einem Reinsto, bei welchen ein Phasenüber-gang zu beobachten ist. Insbesondere will man neben der Temperatur auch die Position der Phasengrenzebestimmen. In dieser Arbeit geht es konkret um die Entwicklung der Phasengrenze zwischen Eis und Wasser,sowie das dazugehörige Temperaturfeld. Beide Variablen entwickeln sich durch die Diusion von Wärme. Zu-sätzlich bedingt latente Wärme den Phasenübergang zwischen Eis und Wasser, wobei an der Phasengrenzeein Sprung im Temperaturgradienten, also im Wärmeuss, zu beobachten ist.

Für die Modellierung des Problems werden folgende Annahmen getroen:

• Betrachtet wird ein Stefan Problem, dabei bezeichnet Ω das Gebiet in der sich das Material in der festenPhase bendet. Entsprechend bezeichnet Ωc das Gebiet der üssigen Phase, siehe dazu Abbildung 1.Die entsprechenden Gleichungen werden in D = Ω ∪ Ωc gelöst.

• Die Phasengrenze wird durch Γ angegeben und die Temperatur an der Phasengrenze ist Tm = const,welches der Schmelztemperatur entspricht.

• Zwischen Eis und Wasser ndet keine signikante Dichteänderung statt, das heiÿt ρw ≈ ρe.

• Die Konvektion im Wasser wird vernachlässigt.

Abbildung 1: Darstellung des zweiphasigen Gebietes

Diese Annahmen führen zu den Gleichungen (1) - (3).

5

3 DIE LEVEL SET METHODE

cs∂T

∂t= ∇(ks∇T ) , x ∈ Ω (1)

cl∂T

∂t= ∇(kl∇T ) , x ∈ Ωc (2)

LV = −(kl∂T

∂n− ks

∂T

∂n

), x ∈ Γ (3)

Dabei bezeichnen cs und cl die Wärmekapazität bei konstantem Volumen in der festen (Index s) beziehungs-weise üssigen (Index l) Phase, ks und kl die Wärmeleitfähigkeit, L die latente Erstarrungswärme, n derNormalvektor der Phasengrenze und V die Geschwindigkeit der Phasengrenze in Normalenrichtung.Als weitere Vereinfachung wird angenommen das die Stokonstanten für die Wärmekapazität, die Wärme-leitung und die latente Erstarrungswärme gleich eins sind und somit aus den Gleichungen (1) - (3) eliminiertwerden können. Dies führt zu den Gleichungen (4) und (5), wobei Gleichung (5) als Stefan Bedingung bekanntist.

∂T

∂t= ∆T , x ∈ D (4)

V = −[∂T

∂n

], x ∈ Γ (5)

Gleichung (4) ist eine partielle Dierentialgleichung zweiter Ordnung in D und somit werden zusätzlichzwei Randbedingungen T0 und T∞, sowie eine Anfangsbedingung zur Schlieÿung des Systems benötigt. DieNotwendigkeit der Stefan Bedingung, Gleichung (5), ergibt sich aus der zusätzlichen Variable, die Positionder Phasengrenze Γ.

3 Die Level Set Methode

3.1 Idee der Level Set Methode

Die Phasengrenze Γ wird durch die Nullstellenmenge einer Hilfsfunktion φ : D × R+ → R beschrieben

Γ = x | φ(x, t) = 0 , x ∈ D (6)

Für die Denition der sogenannten Level Set Funktion φ(x, t) gibt es mehrere Möglichkeiten. In dieserArbeit wird eine signed-distance Funktion verwendet. Die signed-distance Funktion ist so deniert, das sieden euklidischen Abstand d eines Punktes x zur Phasengrenze angibt. Das Vorzeichen von φ(x, t) entscheidetdarüber, ob sich der Punkt x in der üssigen Phase Ωc oder in der festen Phase Ω bendet, siehe Gleichung(7).

φ(x, t) =

+d , x ∈ Ωc

0 , x ∈ Γ

−d , x ∈ Ω

(7)

In Abbildung 2 ist für eine kreisförmige Phasengrenze beispielhaft eine signed-distance Funktion dargestellt.Die linke Abbildung zeigt die signed-distance Funktion selbst, während in der rechten Abbildung ein Con-tourplot der linken Abbildung dargestellt ist. In dem Contourplot kann sehr einfach der Abstand zwischenden einzelnen Contourlinien abgelesen werden.

6

3 DIE LEVEL SET METHODE

Abbildung 2: Beispielplot einer signed distance Funktion mit Contourplot

Die Idee hinter der Level Set Methode ist nun φ(x, t) mit der korrekten Geschwindigkeit V an der Phasen-grenze zu bewegen und das Temperaturfeld T (x, t) mit der neuen, implizit gespeicherten, Phasengrenze Γzu aktualisieren.

3.2 Update der Level Set Funktion

Die Level Set Funktion wird mit Hilfe einer Advektionsgleichung verschoben, wobei die Advektionsgeschwin-digkeit der Schmelz- bzw. Erstarrungsrate entspricht, siehe Gleichung (8). Die Geschwindigkeitsfunktion Fist eine stetige Erweiterung von V auf das gesamte Gebiet D.

φt + F |∇φ| = 0 (8)

Aus Gleichung (5) und der Denition von n = ∇φ|∇φ| folgt

V = − [∇T ]n = − [∇T ]∇φ|∇φ|

(9)

wobei der Sprung in ∇T von der üssigen in die feste Phase berechnet wird, vergleiche Gleichung (3).Die Approximation von [∇T ] basiert auf der Approximation von T in vier Koordinatenrichtungen, die kar-tesischen Koordinaten x und y, sowie die um 45 gedrehten Koordinaten η und ζ, siehe Abbildung 3.

Abbildung 3: Darstellung der vier Koordinatenrichtungen x, y, η und ζ

7

3 DIE LEVEL SET METHODE

Jede Approximation des Sprungs in der Ableitung von T in einer der vier Koordinatenrichtungen kann durcheine Advektionsgleichung stetig von der Phasengrenze aus erweitert werden, siehe Gleichungen (10) - (13).

Dabei ist u1 =[∂T∂x

], u2 =

[∂T∂y

], u3 =

[∂T∂η

]und u4 =

[∂T∂ζ

].

u1t + S(φφx)u1x = 0 (10)

u2t + S(φφy)u2y = 0 (11)

u3t + S(φφη)u3η = 0 (12)

u4t + S(φφζ)u4ζ = 0 (13)

S : R→ −1, 0, 1 ist die Vorzeichenfunktion,

S(x) =

1 , x > 0

0 , x = 0

−1 , x < 0

(14)

Um Instabilitäten zu vermeiden denieren wir F als arithmetisches Mittel zwischen den kartesischen Koor-dinaten x und y und den rotierten Koordinaten η und ζ

F =

4∑i=1

uiφxi

|∇φ|(15)

wobei x1 = x, x2 = y, x3 = η und x4 = ζ gilt.

3.3 (Re)Initialisierung der signed distance Funktion

Für den Start der Berechnung stellt sich die Frage, wie eine signed-distance Funktion φ(x, t) zu vorgegebenerRandkurve Γ initialisiert werden kann. Zusätzlich ist es notwendig die signed-distance Funktion zu reinitia-lisieren, falls durch die Verschiebung der Phasengrenze, siehe Kapitel 3.2, φ keine signed-distance Funktionmehr ist.Durch Iteration der Gleichung

φt = S(φ0)(1− |∇φ|) (16)

kann eine Funktion φ0 zu einer signed distance Funktion entwickelt werden. Die Nullstellenmenge, und damitdie Position der Phasengrenze Γ, von φ0 bleibt erhalten.

3.4 Update des Temperaturfeldes

Die Integration von Gleichung (4) ist in einem einphasigen Gebiet durch einfache Methoden, zum Beispiel einfünf Punkte Dierenzen Schema mit Vorwärtseuler, möglich. In einem zweiphasigen Gebiet gibt es jedoch einProblem durch den Sprung des Temperaturgradienten an der Phasengrenze. Damit ist eine Approximationdes Temperaturgradienten in der Nähe der Phasengrenze durch ein Standard Finiten Dierenzen Schemanicht ohne groÿe Fehler möglich. Da φ eine signed distance Funktion ist kann der Wert dieser Funktionverwendet werden um festzustellen ob sich ein Punkt x in der Nähe der Phasengrenze bendet. Zusätzlichkann mit Hilfe von φ der Abstand eines Punktes x zur Phasengrenze in x beziehungsweise y-Richtunginterpoliert werden. Dadurch ist es möglich zwei Interpolationspolynome PL und PR auf der rechten undlinken Seite, beziehungsweise auf der oberen und unteren Seite, der Phasengrenze zu berechnen. ZweifacheDierenzierung dieser Polynome führt zu der gesuchten Approximation von Txx beziehungsweise Tyy in derNähe der Phasengrenze.

8

4 DISKRETISIERUNG

4 Diskretisierung

Für die Diskretisierung der in Kapitel 2 und 3 vorgestellten Gleichungen wird eine quadratische Box mitSeitenlänge L und äquidistanter Gitterweite h = ∆x = ∆y = L

M betrachtet, wobei (M+1)2 die Gesamtanzahlder Gitterpunkte ist. Die Zeitschrittweite wird mit ∆t bezeichnet. Dazu werden die folgenden Denitionenim restlichen Kapitel verwendet

xi,j = ((i− 1)h, (j − 1)h)

φi,j = φ(xi,j)

Ti,j = T (xi,j)

i, j = 1, ...,M + 1

4.1 Diskretisierung des Updates der Level Set Funktion

Aus Kapitel 3.2 ist hier noch einmal kurz die Vorgehensweise beim Update der Level Set Funktion zusam-mengefasst.

1. Approximation des Temperatursprungs in die vier Koordinatenrichtungen, siehe Kapitel 4.1.1.

2. Stetige Erweiterung des Temperatursprungs ausgehend von der Phasengrenze, siehe Kapitel 4.1.2.

3. Berechnung der Geschwindigkeitsfunktion F , siehe Kapitel 4.1.3.

4. Lösen der Advektionsgleichung für φ, siehe Kapitel 4.1.4.

4.1.1 Approximation des Temperatursprungs

Zur Approximation des Temperatursprungs werden vier Felder u1, u2, u3 und u4, Notation siehe Kapitel 3.2,berechnet. An Gitterpunkten auf oder in der Nähe der Phasengrenze approximieren diese vier Felder denSprung in ∂T

∂x ,∂T∂y ,

∂T∂η und ∂T

∂ζ über die Phasengrenze Γ. An Gitterpunkten die weiter von Γ entfernt sind,

sind u1 − u4 normalerweise nahe an null. In Abbildung 4 ist beispielhaft solch ein Feld für eine Dimensionaufgezeigt. Es ist deutlich zu erkennen das der absolute Wert des Temperatursprungs in der Nähe derPhasengrenze am gröÿten ist.

Abbildung 4: Beispiel für ein Feld des Temperatursprungs in einer Dimension

Die Diskretisierung für u1 − u4 sieht wie folgt aus

9

4 DISKRETISIERUNG

u1i,j = −Si,j(φx)((Ti+2,j − Ti+1,j)− (Ti−1,j − Ti−2,j))/hu2i,j = −Si,j(φy)((Ti,j+2 − Ti,j+1)− (Ti,j−1 − Ti,j−2))/h

u3i,j = −Si,j(φη)((Ti+2,j+2 − Ti+1,j+1)− (Ti−1,j−1 − Ti−2,j−2))/(√

2h)

u4i,j = −Si,j(φζ)((Ti+2,j−2 − Ti+1,j−1)− (Ti−1,j+1 − Ti−2,j+2))/(√

2h)

wobei S erneut die Vorzeichenfunktion ist. Die Vorzeichenfunktion ist notwendig um zu garantieren das derSprung konsistent von der festen in die üssige Phase berechnet wird.

4.1.2 Stetige Erweiterung des Temperatursprungs

Für die stetige Erweiterung des Temperatursprungs auf ganz D werden die Advektionsgleichungen (10) -(13) bis zu einem stationären Endwert gelöst (t− > ∞), siehe dazu auch [2]. Jedoch ist es für die korrekteAusbreitung der Level Set Funktion ausreichend die folgenden Gleichungen für einige wenige Zeitschritte miteiner CFL-Zahl von cfl = 0.5 zu berechnen

u1,(new)i,j = u

1,(old)i,j − cfl ∗ (u

1,(old)i,j − u1,(old)i−1,j ) falls Si,j(φφx) > 0

u1,(new)i,j = u

1,(old)i,j + cfl ∗ (u

1,(old)i+1,j − u

1,(old)i,j ) falls Si,j(φφx) < 0

wobei die Gleichungen für u2 − u4 äquivalent gebildet werden. Für den zweiten Term in den beiden obigenGleichungen ist es auch möglich ein ENO Schema gleicher oder höherer Ordnung für gröÿere Stabilität zuverwenden. Genauere Informationen zu den ENO Schemata ndet sich zum Beispiel in [3].In Abbildung 5 ist solch eine stetige Erweiterung von Abbildung 4 dargestellt. Es ist zu erkennen das derWert für x = 0 (also die Position der Phasengrenze) in Abbildung 4 stetig in die linke und rechte Richtungin Abbildung 5 erweitert wurde.

Abbildung 5: Beispiel für eine stetige Erweiterung eines Feldes des Temperatursprungs in einer Dimension

4.1.3 Berechnung der Geschwindigkeitsfunktion F

Die beiden folgenden Gleichungen gelten gleichermaÿen

Fi,j = u1i,j

(φx|∇φ|

)i,j

+ u2i,j

(φy|∇φ|

)i,j

Fi,j = u3i,j

(φη|∇φ|

)i,j

+ u4i,j

(φζ|∇φ|

)i,j

10

4 DISKRETISIERUNG

Bilden des arithmetischen Mittels liefert

Fi,j =1

2(u1i,j

(φx|∇φ|

)i,j

+ u2i,j

(φy|∇φ|

)i,j

+ u3i,j

(φη|∇φ|

)i,j

+ u4i,j

(φζ|∇φ|

)i,j

)

⇒ Fi,j |∇φ|i,j =1

2(u1i,j (φx)i,j + u2i,j (φy)i,j + u3i,j (φη)i,j + u4i,j (φζ)i,j) (17)

4.1.4 Diskretisierung der Advektionsgleichung

Mit den Gleichungen (8) und (17) folgt

φt = −1

2(u1i,j (φx)i,j + u2i,j (φy)i,j + u3i,j (φη)i,j + u4i,j (φζ)i,j)

welches mit einem Runge Kutta Schema dritter Ordnung in der Zeit und einem ENO Schema zweiter Ordnungim Raum gelöst wird. Wichtig ist, dass diese Gleichung nur für einen Hauptzeitschritt ∆t berechnet wird.Das Runge Kutta Scheme hat für die allgemeine Gleichung φt = L(φ) die folgende Form

φ(1) = φ(0) + ∆tL(φ(0))

φ(2) =3

4φ(0) +

1

4(φ(1) + ∆tL(φ(1)))

φ(3) =1

3φ(0) +

2

3(φ(2) + ∆tL(φ(2)))

wobei L die diskrete Form von L ist.

4.2 Diskretisierung der (Re)Initialisierung der signed distance Funktion

In zwei Dimensionen kann Gleichung (16) geschrieben werden als

φt = S(φ0)(1−√φ2x + φ2y)

Eine Möglichkeit zur Diskretisierung dieser Gleichung ist die sogenannte Godunov Methode

φN+1i,j = φNi,j −∆tSε(φ

0i,j)G(φNi,j)

mit

a ≡ D−x φi,j = (φi,j − φi−1,j)/hb ≡ D+

x φi,j = (φi+1,j − φi,j)/hc ≡ D−y φi,j = (φi,j − φi,j−1)/h

d ≡ D+y φi,j = (φi,j+1 − φi,j)/h

Sε(φi,j) =φi,j√φ2i,j + ε2

und

G(φi,j) =

√max((a+)2, ((b−)2)) +max((c+)2, ((d−)2))− 1 , ifφ0i,j > 0√max((a−)2, ((b+)2)) +max((c−)2, ((d+)2))− 1 , ifφ0i,j < 0

0 , otherwise

11

4 DISKRETISIERUNG

Die Verwendung von Sε statt S dient zur Vermeidung von Divisionen durch null, wobei ε = 2h gewählt wird.Um höhere Ordnung zu erzielen können die einseitigen Dierenzen von φx und φy durch ENO Schematahöherer Ordnung ersetzt werden. Zur Zeitdiskretisierung wird das Runge Kutta Schema aus Kapitel 4.1.4verwendet. Da die hier verwendet Zeitschrittweite ∆treinit unabhängig von der Hauptzeitschrittweite ∆t ist,wird ∆treinit = h/5 gewählt. Nähere Informationen zu dem Schema und zur Bedeutung von a+/− etc. ndensich in [4].

4.3 Diskretisierung der Temperaturgleichung

Zur Lösung der Temperaturgleichung wird ein implizites Verfahren verwendet. In Punkten die sich weit vonder Phasengrenze entfernt benden wird ein Standard fünf Punkte Schema verwendet

Tn+1i,j − Tni,j

∆t=Tn+1i+1,j − 2Tn+1

i,j + Tn+1i−1,j

h2+Tn+1i,j+1 − 2Tn+1

i,j + Tn+1i,j−1

h2(18)

Für Punkte die in der Nähe der Phasengrenze liegen wird zunächst der Abstand zur Phasengrenze selbstberechnet.Beispielhaft sei xf ∈ Γ wobei gelten soll xi,j ≤ xf ≤ xi+1,j . Die Abstände

xi+1,j − xf =

(φi+1,j

φi+1,j − φi,j

)h = r1h

xf − xi,j = −(

φi,jφi+1,j − φi,j

)h = r2h

werden verwendet um zwei Interpolationspolynome PL und PR zu konstruieren. Beide Polynome sind lediglichabhängig von x. PL wird konstruiert mit Hilfe der Variablen h, r2, T (xf ), Ti,j , Ti−1,j etc. und PR mit Hilfeder Variablen h, r1, T (xf ), Ti+1,j , Ti+2,j etc. In Abbildung 6 ist dieses Szenario für PL dargestellt.

Abbildung 6: Beispiel für die Konstruktion von PL

Durch PL kann nun die zweite Ableitung an der Stelle xi,j approximiert werden P ′′L(xi,j) ≈ Txx(xi,j).Eine Möglichkeit diese Polynome implizit zu bestimmen ist mit Hilfe der Vandermonde Matrix V.

V =

1 x1 x21 . . . xn−11

1 x2 x22 . . . xn−12

. . . . . . . . .. . . . . .

1 xn x2n . . . xn−1n

12

4 DISKRETISIERUNG

Die Koezienten c eines Polynoms vom Grad n-1 das die Funktion f(x) interpoliert können mit der Gleichung

V c = f (19)

bestimmt werden, wobei f = (f(x1), ..., f(xn)) ist.Angewendet auf das obige Beispiel ergibt sich

(x1, ..., xn) = (xf , xi,j , xi−1,j , ..., xi−n+2,j)

(f(x1), ..., f(xn)) = (T (xf )n+1, Tn+1i,j , ..., Tn+1

i−n+2,j)

Für ein Polynom dritten Grades (n=4) ergibt sich damit folgende Approximation der zweiten Ableitung

Txx(xi,j) ≈ 6c4xi,j + 2c3 (20)

Durch Verwendung der Gleichungen (19) und (20) können die gesuchten Polynome PL und PR nun implizitbestimmt und damit die gesuchten Approximation für Txx und äquivalent Tyy berechnet werden. Dazuwerden diese Gleichungen in das LGS des Standard fünf Punkte Schemas, siehe Gleichung (18), integriert.Sei

Ax = b

das LGS des Standard fünf Punkte Schemas mit

x = (T (x1,1)n+1, T (x2,1)n+1, ..., T (xM+1,M )n+1, T (xM+1,M+1)n+1)

b = (T (x1,1)n, T (x2,1)n, ..., T (xM+1,M )n, T (xM+1,M+1)n)

Dann ergibt sich ein modiziertes LGS

Ax = b

A =

(A1 A2

A3 A4

)x =

(xc

)b =

(bb2

)wobei A1 = A ist. In der Matrix A2 stehen die Koezienten aus der Approximation der Ableitung, sieheGleichung (20), und in der Matrix A3 werden die Koezienten vor den Temperaturwerten Tn+1 geschrie-ben, siehe Gleichung (19). Die Vandermonde Matrix V bendet sich in der Matrix A4. Der Vektor b2 enthältneben vielen Nulleinträgen lediglich die Werte T (xf ).Eine weitere Möglichkeit der Berechnung ndet sich in [2], wobei hier das LGS Ax = b nicht erweitert werdenmuss.Die Zeititeration wird nun so lange durchlaufen bis die Lösung konvergiert, das heiÿt bis die folgende Un-gleichung

M+1∑i,j=1

(Tn+1i,j − Tni,j) < tol mit tol = 1.0 · 10−12

erfüllt ist.

13

5 ERGEBNISSE

4.4 Algorithmus

Der vollständige Algorithmus sieht nun wie folgt aus

1. Initialisierung der Temperatur T (x, 0) und der Level Set Funktion φ(x, 0)

2. Update der Level Set Funktion φ(x, t)

3. Reinitialisierung der Level Set Funktion φ(x, t)

4. Integration des Temperaturfeldes T (x, t)

5. Wiederholung der Schritte 2)-4) bis zum Endzeitpunkt T

5 Ergebnisse

Eine analytische Lösung des Stefan Problems für V = const ist gegeben durch

T (x, t) =

−1 + e−V (x−V t) , x > V t

0 , x ≤ V t

wobei die Phasengrenze parametrisiert ist durch

Γ(t) = x = V t, y = s , s ∈ R

In Abbildung 7 ist diese Lösung für V = 1 und t = 0 dargestellt.

Abbildung 7: Analytische Lösung für das Stefan Problem

In den Tabellen 5 und 6 sind die Ergebnisse eigener Berechnungen aufgezeigt. Die Rechnungen wurden für dieParameter V = 1, T = 5 ·10−4 und ∆t = 1 ·10−5 durchgeführt. In Tabelle 5 wurden Finite Dierenzen ersterOrdnung sowie das Upwind Verfahren erster Ordnung verwendet, während in Tabelle 6 Finite Dierenzenzweiter Ordnung und ENO Schemata zweiter Ordnung verwendet wurden.

14

5 ERGEBNISSE

Gridpoints calc. Interface Position ||T − Tn||L1 Conv. Rate20 4.1054 · 10−4 1.0429 · 10−4 -40 4.5293 · 10−4 5.7189 · 10−5 0.866880 4.7589 · 10−4 2.9889 · 10−5 0.9361160 4.8785 · 10−4 1.5222 · 10−5 0.9734

Tabelle 5: Finite Dierenzen und Upwind Verfahren erster Ordnung

Gridpoints calc. Interface Position ||T − Tn||L1 Conv. Rate20 4.0836 · 10−5 1.0684 · 10−4 -40 4.5383 · 10−5 5.6105 · 10−5 0.929280 4.7925 · 10−4 2.5728 · 10−5 1.1248160 4.9508 · 10−4 6.1668 · 10−6 2.0607

Tabelle 6: Finite Dierenzen und ENO Schemata zweiter Ordnung

In beiden Tabellen lässt sich eine Konvergenz des Fehlers ablesen, wobei in der letzten Berechnung vonTabelle 6 sogar eine Konvergenzordnung von zwei zu erkennen ist. Hingegen scheint in Tabelle 5 die erwarteteKonvergenzordnung von 1 erreicht zu werden. Auällig ist jedoch der relativ groÿe Fehler in der Position derPhasengrenze. Die analytisch berechnete Position der Phasengrenze ist x = V T = 5 · 10−4. Vor allem füreine geringe Anzahl an Gitterpunkten weicht der berechnete Wert jedoch sehr deutlich ab.

15

LITERATUR

Literatur

[1] S. Chen, B. Merriman, S. Osher, and P. Smereka (1997). A Simple Level Set Method for Solving StefanProblems, Journal of Computational Physics 135, 829 (1997)

[2] Stanley Osher and Ronald P. Fedkiwy (2000). Level Set Methods: An Overview and Some RecentResults, Journal of Computational Physics 169, 463502 (2001)

[3] Chi-Wang Shu (1997). Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemesfor Hyperbolic Conservation Laws, NASA/CR-97-206253 ICASE Report No. 97-65 (1997)

[4] M. Sussman, P. Smereka, and S. Osher, (1994). A level set approach for computing solutions to incom-pressible two-phase ow, Journal of Computational Physics 114, 146 (1994)

16