Sebastian Gatzka
Numerische Simulation der periodischen Blasenbildung an einer Einströmöffnung mit
dem Programm OpenFOAM
Numerical simulation of the periodic bubble formation at an orifice applying the program
OpenFOAM
ISBN: 978-3-8366-2709-2
Herstellung: Diplomica® Verlag GmbH, Hamburg, 2009
Zugl. Technische Universität Darmstadt, Darmstadt, Deutschland, Studienarbeit, 2009
Dieses Werk ist urheberrechtlich geschützt. Die dadurch begründeten Rechte,
insbesondere die der Übersetzung, des Nachdrucks, des Vortrags, der Entnahme von
Abbildungen und Tabellen, der Funksendung, der Mikroverfilmung oder der
Vervielfältigung auf anderen Wegen und der Speicherung in Datenverarbeitungsanlagen,
bleiben, auch bei nur auszugsweiser Verwertung, vorbehalten. Eine Vervielfältigung
dieses Werkes oder von Teilen dieses Werkes ist auch im Einzelfall nur in den Grenzen
der gesetzlichen Bestimmungen des Urheberrechtsgesetzes der Bundesrepublik
Deutschland in der jeweils geltenden Fassung zulässig. Sie ist grundsätzlich
vergütungspflichtig. Zuwiderhandlungen unterliegen den Strafbestimmungen des
Urheberrechtes.
Die Wiedergabe von Gebrauchsnamen, Handelsnamen, Warenbezeichnungen usw. in
diesem Werk berechtigt auch ohne besondere Kennzeichnung nicht zu der Annahme,
dass solche Namen im Sinne der Warenzeichen- und Markenschutz-Gesetzgebung als frei
zu betrachten wären und daher von jedermann benutzt werden dürften.
Die Informationen in diesem Werk wurden mit Sorgfalt erarbeitet. Dennoch können
Fehler nicht vollständig ausgeschlossen werden und der Verlag, die Autoren oder
Übersetzer übernehmen keine juristische Verantwortung oder irgendeine Haftung für evtl.
verbliebene fehlerhafte Angaben und deren Folgen.
© Diplomica Verlag GmbH
http://www.diplomica.de, Hamburg 2009
Inhaltsverzeichnis
Nomenklatur
II
1 Einleitung
1
2 Physikalische Grundgleichungen
5
2.1 Massenerhaltung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.2 Impulserhaltung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.3 Drei-Phasen-Kontaktwinkel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
3 Numerische Verfahren
8
3.1 Numerisches Gitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
3.2 Finite-Volumen-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3.3 Volume-of-Fluid-Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3.1
Grundlagen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.3.2
Vorgehensweise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.3.3
Diskretisierung der Transportgleichung für
. . . . . . . . . . . . . . . . . . . . . . . . . 16
3.4 Berücksichtigung der Grenzflächenspannung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.5 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.6 Zeitschrittweitensteuerung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.7 Lösungsalgorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4 Numerische Testprobleme
20
4.1 Der gebrochene Damm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2 Einfluss der Grenzflächenspannung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2.1
Parasitäre Strömungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
4.2.2
Durchgeführte Simulationen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.3 Aufsteigende Gasblase in ruhender Flüssigkeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.4 Zusammenfassung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5 Blasenbildung an einer Einströmöffnung
43
5.1 Problembeschreibung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.2 Ablauf der Simulation und Simulationsparameter . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.3 Problemgebiet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.4 Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.5 Erweiterung des Modells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6 Zusammenfassung
76
Literaturverzeichnis
78
Anhang
80
I
Nomenklatur
Generell sind im folgenden Text vektorielle Größen durch Fettdruck charakterisiert. Somit ist z.B. zwi-
schen dem vektoriell geschriebenem Spannungstensor
und der skalaren Grenzflächenspannung zu
unterscheiden.
Lateinische Buchstaben
a
Breite des Wasservolumens
[m]
B
Breite des Rechengebietes
[m]
Bo
Bond-Zahl, Eötvös-Zahl
[-]
C o
Courant-Zahl
[-]
d
Durchmesser
[m]
d
e
äquivalenter Durchmesser
[m]
f
Kraft infolge der Grenzflächenspannung
N
m
3
f
Frequenz
1
s
g
Gravitationskonstante
m
s
2
H
Höhe des Rechengebietes
[m]
H
dimensionslose Position der Wasserhöhe
[-]
h
Höhe des Wasservolumens
[m]
h
Zellengröße, äquidistante Zellengröße
[m]
I
Einheitsmatrix
[-]
K
dimensionslose Konstante,
Qualitätsmerkmal zur Beurteilung des CSF-Modells
[-]
M o
Morton-Zahl
[-]
n
Normalenvektor
[-]
n
Stabilisierter Normalenvektor
[-]
n
Zeitschritt
[-]
n
Anzahl an Unterzyklen des MULES-Schemas
[-]
N
Anzahl
[-]
p
mittlerer Druck
N
m
2
II
p
Druck
N
m
2
p
d
modifizierter Druck
N
m
2
Q
Einströmrate
ml
min
r
Radius
[m]
R
0
Blendenradius
[m]
r
e
äquivalenter Radius
[m]
Re
Reynolds-Zahl
[-]
S
f
Normalenvektor der Zellenseite
m
2
S
Oberfläche
m
2
s
Standardabweichung, empirische Varianz
hier:
N
m
2
t
Zeit
[s]
t
B
Blasenbildungszeit
[s]
T
Temperatur
[K]
u
Geschwindigkeit
m
s
U
dimensionslose Geschwindigkeit
[-]
u
Betrag des Geschwindigkeitsfeldes
m
s
2
u
r
künstliche Kompressiongeschwindigkeit
m
s
u
ma x
maximale Geschwindigkeit der parasitären Strömungen
m
s
V
Volumen
m
3
t
Tangentialvektor
[-]
X
dimensionslose Position der Wasserfront
[-]
x
k
Position der Drei-Phasen-Konstaktlinie beim Simulationsstart
[m]
III
Griechische Buchstaben
Drei-Phasen-Kontaktwinkel
[
o
]
Stabilisierungsparameter
1
m
Übergangsbereich
[m]
t
Veränderung der Blasenbildungszeiten
[%]
V
Veränderung der Blasenvolumen
[%]
VOF-Variable
[-]
Krümmung
1
m
µ
dynamische Viskosität
kg
ms
kinematische Viskosität
m
2
s
r
volumetrischer Fluss der künstlichen
Kompressionsgeschwindigkeit ?
m
3
s
Drei-Phasen-Kontaktwinkel
[
o
]
Dichte
kg
m
3
Spannungstensor
N
m
2
Grenzflächenspannung
N
m
Reibungsspannungstensor
N
m
2
dimensionslose Geschwindigkeit bezogen auf U
[-]
dimensionslose Zeit bezogen auf H
[-]
T
dimensionslose Zeit bezogen X
[-]
IV
Indizes
B
Blase
c
numerische Zelle
E
Ethanol
f
Mittelpunkt der Zellenseite
G
Ergebnisse aus dem Modell von Gerlach et. al. [7]
l
Flüssigkeit
M
Mittelpunkt
N
Mittelpunkt der Nachbarzelle
P
Mittelpunkt der numerischen Zelle
s
Festkörper
sim
Simulation
v
Dampf
W
Wasser
b
Rechengebietsrand
i
Grenzfläche
R
Simulation mit Einlasskanal
sat
Sättigungsbedingungen
w
Wand
1
x
1
-Richtung
1
Fluid 1
2
x
2
-Richtung
2
Fluid 2
V
Abkürzungen
CDS
Central Differencing Scheme
CFD
Computational Fluid Dynamics
CLSVOF
Coupled Leve-Set- and Volume-of-Fluid-Method
CSF
Continuum Surface Force
FVM
Finite-Volumen-Methode
MULES
Multidimensional universal limiter with explicit solution
OpenFOAM
Open Field Operation and Manipulation
PISO
Pressure-Implicit Split-Operator
TVD
Total variation diminishing
VOF
Volume-of-Fluid
Superskripte
0
Einheitsvektor
VI
1 Einleitung
Die hohe Nachfrage nach leistungsfähigen Wärmeübertragern für die Kühlung von elektronischen Bautei-
len, zur Effizienzsteigerung von Kreisprozessen für Wärmekraft- und Kälteanlagen, sowie für die chemi-
sche oder verfahrenstechnische Industrie oder der Umwelttechnik haben die Zweige der Thermodynamik
sowie der Wärme- und Stoffübertragung, die sich mit Phasenwechselphänomenen beschäftigen, in den
letzten Jahren deutlich angekurbelt.
Bei einem Phasenwechsel von Flüssigkeit zu Dampf (Verdampfen) oder von Dampf zu Flüssigkeit (Kon-
densation) findet man sehr hohe Wärmeübergangskoeffizienten vor, da das Fluid sehr viel Wärme durch
seine Verdampfungsenthalpie aufnehmen kann. Verdampfungs- und Kondensationsprozesse sind deshalb
viel effizienter als einphasige Prozesse.
Am effizientesten sind Verdampfungsprozesse, wenn Blasensieden vorliegt. Dabei werden die höchsten
Wärmeübertragungskoeffizienten beobachtet.
Das Blasensieden ist ein Phänomen beim Behältersieden. Dabei wird ein offener und zum Teil mit Flüs-
sigkeit gefüllte Behälter von unten beheizt. Zur Verdampfung der Flüssigkeit kommt es, wenn die Wand-
temperatur T
w
größer als die Sättigungstemperatur T
sat
der Flüssigkeit zuzüglich einer möglicherweise
notwendigen Überhitzung ist [24]. Der sich einstellende Wärmetransportmechanismus hängt maßgeb-
lich davon ab, wie hoch diese Wandüberhitzung
T = T
w
- T
sat
(1.1)
ist. Die verschiedenen Siedebereiche sind in der Nukijama-Kurve (Abb. 1.1) aufgetragen.
Abbildung 1.1: Nukijama-Kurve für Wasser [24]
1
Wie zu erkennen ist, findet man Blasensieden bei Wandüberhitzungen im Bereich zwischen fünf und 30
Kelvin vor. In diesem Bereich steigt der übertragene Wärmestrom mit größer werdender Wandüberhit-
zung steil an. Wo im Gegensatz zum stillen Sieden die Verdampfung noch an der Flüssigkeitsoberfläche
stattgefunden hat, findet sie beim Blasensieden direkt an der Heizfläche statt. Die Dampfblasen entste-
hen dabei fast ausschließlich in kleinen Vertiefungen in der Heizoberfläche. Da alle technischen Ober-
flächen eine gewissen Rauhigkeit aufweisen ist es nicht auszuschließen, dass kleine Gaseinschlüsse in
den oft mikroskopisch kleinen Rauhigkeitsvertiefungen nicht vollständig ausgespült werden [24]. Abb.
1.2 zeigt eine schematische Darstellung solcher Gaseinschlüsse. Da die Dampfblasen später aus diesen
Gaseinschlüssen ,,keimen", spricht man oft auch von Keimstellen.
Abbildung 1.2: Stark vergrößerte Darstellung der Oberflächenstruktur einer typischen metallischen Heiz-
fläche [16]
Betrachtet man das thermodynamische Gleichgewicht einer solchen Keimstelle, kann gezeigt werden,
dass bei einer konstanten Wandüberhitzung nur Keimstellen wachsen können, die eine bestimmte Grö-
ße überschreiten. Man spricht in diesem Zusammenhang von Keimstellenaktivierung. Auf die Herleitung
dieser Tatsache wird an dieser Stelle nicht weiter eingegangen. Es wird vielmehr auf die einschlägige
Literatur ([1], [24]) verwiesen.
Nach der Aktivierung wächst das Volumen der Dampfblasen mit der stetigen Verdampfung solange an,
bis die Auftriebskräfte nicht mehr im Gleichgewicht mit den Trägheits-, Grenzflächen- und Haftkräften
stehen. Die Blase beginnt aufgrund des Dichteunterschieds zur umgebenden Flüssigkeit aufzusteigen.
Vorher schnürt sich die Blase am unteren Teil (dem Blasenfuß) ein, bis sie sich dort durchtrennt. Die
Blase kann nun infolge der Auftriebskräfte ungehindert, entgegen der Gravitation, nach oben steigen
(vgl. Abb. 1.3).
Abbildung 1.3: Abreißvorgang nach Mitrovi´
c. a Form der haftenden Blase;
b Einschüren und Abreißen der Blase; c Kondensation des Dampfrests [1]
2
Manglik [19] betont in einer aktuellen Veröffentlichung, dass aufgrund deren Komplexität noch kei-
ne ganzheitliche Theorie für Siedevorgänge, Mehrphasenströmungen und Grenzflächenphänomene in
Aussicht ist. Es existieren zwar einige Korrelationen, die als Berechnungsgrundlage dienen können, die
allerdings auf spezielle Anwendungsfälle beschränkt sind. Einen kompakten Überblick über die Menge
an Parametern, die einen Einfluss auf die Modellierung des Blasensiedens haben zeigt Abb. 1.4.
Abbildung 1.4: Parameter beim Blasensieden [19]
Die vorliegende Arbeit kann sich nicht mit der gesamten Komplexität des Blasensiedens beschäftigen.
Vielmehr soll hier nur der Vorgang des periodischen Blasenwachstums bis zum Abriss und das anschlie-
ßende Aufsteigen in der umgebenden Flüssigkeit betrachtet werden. Dazu wird eine numerische Si-
mulation mit dem CFD-Berechnungsprogramm OpenFOAM
1
Version 1.4.1 erstellt. Die zentrale Verein-
fachung, die das Gesamtproblem in ein geeignetes Modell überführt, ist das Tatsache, dass der Ver-
dampfungsprozess selbst nicht mitsimuliert wird. Vielmehr wird die Verdampfung dadurch modelliert,
dass ein Gas von unten durch eine Blendenöffnung in das mit Flüssigkeit gefüllte Rechengebiet ein-
strömt. In diesem Fall muss auch die Energiegleichung nicht gelöst werden, da die Temperatur keine
Rolle spielt. Für einen solchen Anwendungsfall mit zwei inkompressiblen Phasen steht in OpenFOAM
der Löser interFoam zur Verfügung.
In einem sehr viel komplexeren Modell [6] wurde bereits das periodische Anwachsen und Abreißen
einer Einzelblase an einer Heizwand erfolgreich simuliert. Dabei wurde jedoch ein randangepasstes nu-
merische Gitter verwendet. Dabei fällt die Grenzfläche immer mit Elementseiten zusammen. Bei einer
solchen Vorgehensweise können allerdings Topologieänderungen (z.B. durch Koaleszenz zweier Blasen)
nicht realisiert werden. Um solche Vorgänge berücksichtigen zu können, kann ein raumfestes Gitter
eingesetzt werden. Zur Verfolgung der Grenzfläche wird dann ein zusätzliches Verfahren benötigt. In
OpenFOAM ist zu diesem Zweck die Volume-of-Fluid-Methode (VOF) verfügbar.
In diesem Sinne ist die vorliegende Arbeit also Vormodell für ein numerisches Siedemodell zu verstehen,
in dem untersucht wird, ob die dynamischen Blasenwachstums- und Abrissvorgänge ohne Einbeziehung
der eigentlichen Verdampfung akkurat simuliert werden können.
1
OpenFOAM: Open Field Operation And Manipulation.
3
Im Folgenden werden die mathematischen Grundlagen vorgestellt und deren numerische Umsetzung
und Implementierung in OpenFOAM genauer beleuchtet. Für eine erste Beurteilung der Leistungsfähig-
keit des Berechnungsprogrammes werden drei Testfälle simuliert, für die Referenzergebnisse existieren.
Im Einzelnen sind das Testrechnungen mit dominierenden Grenzflächen-, Trägheits- und Auftriebskräf-
ten, den Kräftegruppen wie man sie typischerweise auch beim Blasensieden vorfindet.
Als zentraler Punkt dieser Arbeit werden verschiedene Simulationen zum Blasenwachstum mit unter-
schiedlichen Parametern berechnet. Von besonderem Interesse sind dabei die Zusammenhänge zwischen
dem einströmenden Gasmassenstrom und dem Kontaktwinkel der Drei-Phasen-Kontaktlinie auf die Bla-
senwachstumszeit und das entstehende Blasenvolumen. Die Ergebnisse werden anschließend sowohl
qualitativ als auch quantitativ mit vorliegenden numerischen Ergebnissen aus dem vergleichbaren Mo-
dell von Gerlach et al. [7] verglichen.
4
2 Physikalische Grundgleichungen
In der Kontinuumsmechanik bilden die Erhaltungssätze für Masse, Impuls und Energie die mathemati-
sche Grundlage zur Beschreibung von Strömungsvorgängen. Dies gilt auch für die Klasse an Mehrpha-
senströmungen mit zwei inkompressiblen Phasen, wie sie in der vorliegenden Arbeit betrachtet werden.
Wie bereits in der Einleitung erwähnt, wird im betrachteten Modell des Blasensiedens die eigentliche
Verdampfung nicht mitsimuliert. Deshalb muss bei dieser isothermen Betrachtungsweise die Energieglei-
chung nicht gelöst werden. Im Folgenden beschränken wir uns also auf die Darstellung der Massen- und
Impulserhaltung. Die Gleichungen werden dabei in ihrer konservativen Form in Operatorschreibweise
angegeben. Es sei dabei angemerkt, dass bei der Verwendung der konservativen Formen in einer Finite-
Volumen-Methode die globalen Erhaltungsprinzipien gesichert bleiben [5]. Bei der Diskretisierung von
nicht-konservativen Gleichungen kann es z.B. vorkommen, dass sich die Massenflüsse über die Kontroll-
volumenseiten nicht gegenseitig aufheben. In Folge dessen wäre die lokale und globale Massenerhaltung
nicht erfüllt [11].
2.1 Massenerhaltung
Der Erhaltungssatz der Masse lässt die Interpretation zu, dass die zeitliche Änderung der Masse in einem
Volumen gleich der über die Oberfläche ein- und austretenden Masse ist [22]. Die differentielle Form der
Massenerhaltungsgleichung wird auch Kontinuitätsgleichung genannt:
t
+ · u = 0
(2.1)
Für ein inkompressibles (volumenbeständiges) Fluid wird die Dichte als konstant angenommen. Die
Kontinuitätsgleichung vereinfacht sich in diesem Fall zu
· u = 0
(2.2)
5
2.2 Impulserhaltung
Nach dem Impulssatz ist die zeitliche Änderung des Impulses eines Körpers gleich der Summe aller auf
den Körper wirkenden Oberflächen- und Volumenkräfte [22].
(u)
t
+ · uu = · + g + f
(2.3)
Darin ist der Spannungstensor
enthalten. Er setzt sich zusammen aus Komponenten des Druckes p
und der Reibungsspannung
= -pI +
(2.4)
Der Term f
in (2.3) beschreibt die Kraft infolge der Grenzflächenspannung an der Grenzfläche zwischen
den zwei Phasen. Auf Ihn wird in Kapitel 3.4 noch einmal eingegangen.
Setzt man den Spannungstensor (2.4) in Gleichung (2.3) ein folgt daraus
(u)
t
+ · uu = -p + · + g + f
(2.5)
Weiterhin gilt für inkompressible Newtonsche Fluide der Reibungsspannungstensor [23]
= µ (u) + (u)
T
(2.6)
Bindet man das Materialgesetz (2.6) in Gleichung (2.5) ein, erhält man die die Navier-Stokes-Gleichungen
(u)
t
+ · uu = -p + · µ u + (u)
T
+ g + f
(2.7)
Implementierung in interFoam
Die Navier-Stokes-Gleichungen (2.7) wird in interFoam in einer veränderten Form implementiert. Für
eine effizientere numerische Berechnung, wird die Divergenz des Reibungsspannungstensors · durch
geeignete Operationen umgeformt [21]
· µ u + (u)
T
= · µu + (u) · µ
(2.8)
Aus der Differenz des Druckes p und des hydrostatischen Druckes
g · x wird ein modifizierter Druck
p
d
= p - g · x
(2.9)
bzw. dessen Gradient
p
d
= p - g - g · x
(2.10)
in (2.7) eingeführt. Die implementierten Navier-Stokes-Gleichungen lauten damit
(u)
t
+ · uu = -p
d
+ · µu + (u) · µ - g · x + f
(2.11)
6
2.3 Drei-Phasen-Kontaktwinkel
Trifft die Phasengrenze auf eine Wand können folgende Phänomene beobachtet werden (vgl. [25]). Die
Kräfte zwischen den Molekülen der Flüssigkeit und des Feststoffes werden als Adhäsionskräfte (Haft-
kräfte) bezeichnet. Wenn diese Adhäsionskräfte größer als die Kohäsionskräfte
1
sind, drängen sich die
Moleküle in Richtung des Feststoffes. Die Kontaktfläche zwischen Flüssigkeit und Feststoff vergrößert
sich. Man sagt die Flüssigkeit benetzt den Feststoff. Das Gegenteil trifft zu, wenn die Kohäsionskräfte
größer als die Ahäsionskräfte sind. In diesem Fall verkleinert sich das Kontaktgebiet zwischen Flüssigkeit
und Feststoff.
Abbildung 2.1: Grenzflächenspannungen an der Drei-Phasen-Kontaktlinie.
Das Kräftegleichgewicht (vgl. Abb. 2.1) an der Drei-Phasen-Kontaktlinie führt auf die Young-Laplace-
Gleichung
l v
cos
() =
sv
-
sl
(2.12)
Darin tauchen neben dem Drei-Phasen-Kontaktwinkel
und der Grenzflächenspannung zwischen Flüs-
sigkeit und Gas
l v
auch die Grenzflächenspannungen zwischen Flüssigkeit und Festkörper
sl
und die
Grenzflächenspannung zwischen Festkörper und Gas
sv
auf. Die Kontaktwinkel zwischen drei Stoffpaar-
kombinationen werden experimentell bestimmt.
Bei
= 90
o
steht die Phasengrenze senkrecht auf der Wand. Für
90
o
benetzt die Flüssigkeit die
Wand.
Abbildung 2.2: Beispiele für verschiedene Kontaktwinkel (
= ) [13]
1
Die intermolekularen Bindekräfte innerhalb der Flüssigkeit.
7
3 Numerische Verfahren
Zur numerischen Simulation der Testfälle und der Blasenbildung an der Einströmöffnung wird die Soft-
ware OpenFOAM verwendet. Sie wird von OpenCFD Limited entwickelt, und ist im Gegensatz zu kom-
merzieller Software wie FLUENT oder StarCD im Rahmen freier Softwarelizensierung kostenfrei verfüg-
bar. Die Diskretisierung der beschreibenden Gleichungen erfolgt in OpenFOAM mit der Finite-Volumen-
Methode (FVM). Um die Position der Phasengrenze zu berechnen, wird die Volume-of-Fluid-Methode
(VOF) verwendet. Die beiden Verfahren finden ihre Anwendung im interFoam-Löser, der Teil des Open-
FOAM Programmpakets ist und speziell für die Berechnung von inkompressiblen Zweiphasenströmungen
entwickelt wurde.
3.1 Numerisches Gitter
Grundlage der genannten Verfahren (FVM VOF) bildet ein numerisches Gitter. Im Rahmen der Gitter-
generierung wird das betrachtete kontinuierliche Problemgebiet in eine endliche Anzahl einfacher Teil-
gebiete - die sogenannten Kontrollvolumen oder Zellen - zerlegt. Bei Problemen der Strömungsmechanik
werden dabei bevorzugt strukturierte Gitter aus Hexaedern verwendet [22]. Auch in der vorliegenden
Arbeit wird bis auf eine Ausnahme dieser Gittertyp verwendet.
In OpenFOAM werden jedoch alle Gitter aus Polyedern, also Zellen mit einer unbegrenzten Anzahl an
Seitenflächen, erzeugt. Durch diese Allgemeingültigkeit der verwendeten Netze lassen sich natürlich
auch Spezialfälle wie Hexaeder oder Prismen erzeugen. Diese Freiheit im Umgang mit numerischen Net-
zen macht es auch möglich, gängige Netzformate der kommerziellen Konkurrenz nach OpenFOAM zu
konvertieren.
Alle Gitter werden darüber hinaus ausschließlich in einem dreidimensionalen kartesischen Koordinaten-
system erzeugt. Um zweidimensionale oder achsensymmetrische Berechnungen durchführen zu können,
müssen spezielle Randbedingungen für die Dimension angegeben werden, für die keine Lösungen benö-
tigt werden [18]. Darauf wird in konkreten Fällen noch einmal eingegangen.
Weiterhin werden in OpenFOAM alle Variablen in den Mittelpunkten der Gitterzellen abgespeichert [17].
Man spricht in diesem Fall auch von zellenorientierer Anordnung der Variablen. Diese Wahl der Varia-
blenspeicherung hat sich als vorteilhaft hinsichtlich des rechnerischen Aufwandes und bei komplexen
Rechengebieten erwiesen ([21], [22]).
Abbildung 3.1: Zwei numerische Zellen in OpenFOAM [17]
8
3.2 Finite-Volumen-Verfahren
Liegt das numerische Gitter vor, werden die Erhaltungsgleichungen über alle Kontrollvolumen sowie die
Zeit integriert. Die meisten auftretenden Volumenintegrale werden unter Anwendung des Gaußschen
Integralsatzes in Oberflächenintegrale überführt. Die Volumen- und Oberflächenintegrale werden dar-
aufhin numerisch approximiert und liefern für jedes Kontrollvolumen eine algebraische Gleichung. Die
einzelnen Gleichungen aller Kontrollvolumen bilden jeweils einen Beitrag zu einem globalen algebrai-
schen Gleichungssystem, welches wiederum durch ein geeignetes numerisches Verfahren gelößt werden
kann. Nach der Lösung liegt für jede Variable in den Kontrollvolumen ein diskreter Wert vor. Man spricht
deshalb bei der oben genannten Vorgehensweise auch von Diskretisierung. Bei transienten Problemen
wird dazu noch der betrachtete Zeitraum in eine endliche Anzahl kleiner Zeitintervalle
t unterteilt, für
die jeweile eine Lösung berechnet wird.
Abbildung 3.2: Finite-Volumen-Diskretisierung (örtlich und zeitlich) [17]
Für die Behandlung von inkompressiblen Zweiphasenströmungen (so wie für alle anderen Anwendun-
gen auch) stellt OpenFOAM eine beispielhafte Berechnung zur Verfügung. Darin sind die für diesen
Fall empfohlenen Einstellungen an Diskretisierungsschemata, Stabilitätsbedingungen, Zeitschrittweiten-
steuerung etc. voreingestellt. Diese Einstellungen wurden für alle folgenden Berechnungen übernommen
und sollen hier kurz dargestellt werden.
9
Diskretisierungsverfahren
Im Folgenden werden die verwendeten Diskretisierungsverfahren kurz vorgestellt.
Interpolationsverfahren
Für eine Finite-Volumen-Methode ist es typisch, dass nach der Überführung von Volumenintegralen in
Oberflächenintegrale die in den Gleichungen auftauchenden Variablen an den Zellenseitenmittelpunk-
ten gesucht sind. Da die Variablen bzw. Integranden prinzipiell jedoch in den Mittelpunkten der Zellen
definiert sind (vgl. Kap. 3.1), müssen die Werte an den Zellenseitenmitten aus den Mittelpunkten der
Zellen interpoliert werden, die sich die entsprechende Zellenseite teilen.
1
Für diese sogenannte Cell-to-
face-Interpolation steht eine Vielzahl an Möglichkeiten zur Verfügung. In der vorliegenden Arbeit wird
jedoch, wenn nicht anders erwähnt, ein lineares Zentraldifferenzenverfahren (CDS) verwendet [18].
Gradient
Die Oberflächenintegrale können in die Summen über alle Kontrollvolumenseiten aufgeteilt werden. Am
Beispiel des Druckgradienten in (2.11) werden die einzelnen Schritte noch einmal deutlich.
V
p
d
d V
=
S
p
d
dS
f
S
f
· p
d
, f
(3.1)
Wie zu erkennen ist, werden im approximierten Integral nicht mehr die Variablenwerte in den Zellenmit-
ten, sondern auf den Zellenseiten (Index f ) evaluiert. Multipliziert mit dem (nach außen gerichteten)
Vektor der Kontrollvolumenseite S
f
(vgl. Abb. 3.1) ergibt sich damit ein Fluss der Variable (hier p
d
)
durch die Kontrollvolumenseite. Zur Interpolation dieser gesuchten Größe auf der Kontrollvolumenseite
kommt dann das oben genannte CDS-Verfahren zum Einsatz.
Laplace-Term
Bei der Diskretisierung des Laplace-Terms
V
· µu dV =
S
dS
· µu =
f
µ
f
S
f
· (u)
f
(3.2)
muss sowohl für die dynamische Viskosität
µ
f
an der Kontrollvolumenseite also auch für die Interpo-
lation des Gradienten normal zur Kontollvolumenseite S
f
· (u)
f
ein Interpolationsverfahren gewählt
werden.
Für
µ
f
wird dabei das oben erwähnte Zentraldifferenzenverfahren verwendet. Für S
f
· (u)
f
wird auch
ein Zentraldifferenzenverfahren, aber mit einer Korrekturmöglichkeit im Fall nicht-orthogonaler Gitter
gewählt.
1
Es sei an dieser Stelle darauf hingewiesen, dass zur Unterscheidung der Variablen in der Zellenmitte
P
und an den Zellen-
seitenmitten
f
unterschiedliche Indizes verwendet werden.
10
Konvektionsterm
V
· uu dV =
S
dS
· uu =
f
S
f
· uu
f
(3.3)
Für die Interpolation von
uu
f
wird ein lineares TVD-Schema
2
verwendet. Es dient zur Sicherung der
Beschränktheit und der Vermeidung von ungewollten Oszillationen in der Lösung [11].
Zeitableitung
t
V
udV =
uV
n
- uV
n
-1
t
(3.4)
Die Zeitableitung (3.4) wird mit dem impliziten Euler-Verfahren approximiert. n bezeichnet dabei einen
Zeitschritt,
t die Zeitschrittweite. Es ist von 1. Ordnung und stabil für alle Zeitschrittweiten.
Quellterme
V
gdV =
P
g
P
· V
c
(3.5)
Quellterme der Form (3.5) werden durch die Mittelpunktsregel auf den Wert im Mittelpunkt des Kon-
trollvolumens approximiert.
2
TVD = T otal variation d iminishing.
11
3.3 Volume-of-Fluid-Methode
Zweiphasenströmungen zeichnen sich durch das Vorhandensein einer Grenzfläche oder Phasengrenze zwi-
schen den beiden Phasen aus. Bei der numerischen Behandlung solcher Probleme ist es Teil der Lösung
die Lage der Grenzfläche zu berechnen. Die numerischen Verfahren zur Lokalisierung der Grenzfläche
lassen sich grob in zwei Hauptkategorien einteilen:
· Interfaceverfolgungsmethoden. Dabei liegt die freie Oberfläche als scharfe Trennung zwischen den
zwei Phasen vor. Dafür wird ein randangepasstes Gitter verwendet, welches sich in jedem Zeit-
schritt mit der Grenzfläche bewegt (Abb. 3.3(a)).
· Interfaceerfassungsmethoden. Auf einem festen Rechengitter wird die Position der Grenzfläche aus
dem Volumenverhältnis einer der beiden Phase in den Kontrollvolumen bestimmt. Die Grenzfläche
liegt in Folge dessen nicht scharf vor, sondern ist verschmiert. Um trotzdem ausreichend dünne
Grenzflächen berechnen zu können, muss die Gitterweite entsprechend kleiner gewählt werden,
als bei Interfaceverfolgungsmethoden notwendig wäre. Zu den häufig verwendeten Interfaceerfas-
sungsmethoden gehören die Marker-and-Cell-Methode (MAC) und die Volume-of-Fluid-Methode
(VOF). Bei der MAC-Methode werden masselose Partikel zu Beginn der Berechnung auf der Grenz-
fläche verteilt. Deren Bewegung während der Berechnung dient dann zur Erfassung der Grenz-
fläche (Abb. 3.3(b)). Bei der VOF-Methode wird die Position der Grenzfläche durch Lösen einer
Transportgleichung für den Volumenanteil einer Phase gefunden (Abb. 3.3(c)).
Abbildung 3.3: Methoden zur Berechnung der freien Oberfläche: a) Randangepasstes Netz; b) MAC-
Methode; c) VOF-Methode [21]
Die MAC-Methode bringt sehr hohe Anforderungen hinsichtlich des Speicherplatzes mit sich, da eine
große Anzahl an Koordinaten der Marker-Partikel abgespeichert werden müssen [12]. Wie bereits in der
Einleitung erwähnt wurde, können mit einem konturangepasstem Netz keine Topologieänderungen be-
rücksichtigt werden. Da jedoch bereits der Abreißvorgang einer Blase mit einer Topologieänderung ver-
bunden ist, wird für die folgenden Berechnungen die VOF-Methode verwendet, die im Löser interFoam
implementiert ist.
12
3.3.1 Grundlagen
In einer numerischen Zelle wird der Volumenanteil, den die Phase 1 bzw. 2 einnimmt, durch die Variable
beschrieben:
=
V
1
V
c
= 1 -
V
2
V
c
(3.6)
Für
sind drei Zustände möglich:
·
= 1. Die Zelle ist vollständig mit Fluid der Phase 1 gefüllt.
·
= 0. Die Zelle ist vollständig mit Fluid der Phase 2 gefüllt.
· 0
1. Die Zelle enthält die Grenzfläche zwischen Phase 1 und 2 sowie anteilsmäßig Fluid aus
beiden Phasen.
In der Regel tendieren die numerischen Verfahren zur Lokalisierung der Phasengrenze, die darüber hin-
aus die Beschränktheit zwischen 0 und 1 garantieren sollen, dazu, die Phasengrenze über mehrere Zellen
zu verschmieren [25]. Die Phasengrenze liegt somit nicht als unendlich dünne Grenzfläche vor, sonderen
zeichnet sich durch einen Übergangsbereich
aus, in dem sich kontinuierlich verändert.
Fluid 1
Fluid 2
= 0
= 1
= 0.25
= 0.5
= 0.75
b
a
(a) Grenzfläche mit
-Konturlinien
0.5
1
0
0.25
0.75
b
a
(b)
entlang der Strecke ab
Abbildung 3.4: Schematische Darstellung des Übergangsbereiches der Phasengrenze.
Nach einer Abbildung in [25].
13
Abbildung 3.5: Verschmierte Grenzfläche mit drei Konturlinien für
Die Impulsbilanz (2.11) wird bei Verwendung der Volume-of-Fluid-Methode nicht für jede Phase einzeln
gelöst. Dichte
und Viskosität µ = werden vielmehr als gewichtete Mittelwerte aufgefasst
=
1
+ (1 - )
2
(3.7)
µ = µ
1
+ (1 - )µ
2
(3.8)
und die Gleichung mit diesen Zustandswerten gelöst.
14
3.3.2 Vorgehensweise
Um die Position der Phasengrenze zu finden, haben Hirt Nichols [12] die Volume-of-Fluid-Methode
entwickelt. Die zugrundeliegende Überlegung dieser Methode ist die Tatsache, dass die Normale zur ge-
suchten Oberfläche in der Richtung liegt, in der sich
am stärksten ändert. Für diese Information kann
der Gradient von
herangezogen werden (Abb. 3.6(a)). Unter Hinzunahme der entsprechenden Werte
für
in den numerischen Zellen kann somit der Verlauf der Phasengrenze in der Zelle gefunden werden
(Abb. 3.6(b)).
(a) Tatsächliche
Grenzflächenposition,
Fluidverteilung und Grenzflächennorma-
len
(b) Diskret vorliegende Grenzfläche
Abbildung 3.6: Schematische Darstellung zur Positionsbestimmung der Phasengrenze bei Verwendung
einer VOF-Methode. (Nach einer Abbildung in [25])
Das
-Feld und somit auch die Bewegung der Phasengrenze im Rechengebiet muss nach Hirt Nichols
durch die Lösung einer zusätzlichen Transportgleichungen, die ausschließlich konvektive Terme enthält,
berechnet werden:
t
+ · (u) = 0
(3.9)
In OpenFOAM wird (3.9) jedoch durch einen weiteren Term zur künstlichen Komprimierung der Grenz-
fläche erweitert.
t
+ · (u) + · u
r
(1 - ) = 0
(3.10)
Darin stellt u
r
ein geeignetes Geschwindigkeitsfeld zur Komprimierung der Phasengrenze dar, um deren
Verschmierung zu verringern. Durch
(1 - ) wird sichergestellt, dass u
r
nur im Bereich der Phasengren-
ze aktiv ist und das Strömungsfeld ausserhalb dieser Region nicht beeinflusst wird ([10], [21]).
15
Es sei bemerkt, dass (3.10) eine nicht-lineare Gleichung für
darstellt. Das macht ein iteratives Lösungs-
verfahren für diese Gleichung notwendig. Für größere Courant-Zahlen
C o
=
|u|t
h
(3.11)
wurden dabei allerdings Konvergenzprobleme beobachtet [21]. Bei den transienten Berechnungen, wie
sie in der vorliegenden Arbeit zur Anwendung kommen, ist die Courant-Zahl allerdings klein genug um
akkurate Lösungen zu erhalten. Konvergenzeprobleme bei der Berechnungn von (3.10) sind somit nicht
zu erwarten.
Hirt Nichols haben bereits angemerkt, dass ihre Betrachtungen zur Volume-of-Fluid-Methode ohne
eine geeignete Lösungsmethode für die Transportgleichung (3.10) nutzlos sind. Gleichung (3.10) wird
in interFoam mit dem MULES-Schema
3
gelöst. MULES ist ein Verfahren zur Lösung von Transportglei-
chungen mit rein konvektiven Termen. Es wurde speziell für den Einsatz mit interFoam entwickelt, um
die Beschränktheit von
unabhängig von allen anderen verwendeten numerischen Verfahren, der Netz-
struktur o.ä. zu garantieren. Damit ist die Wahl des numerischen Verfahren für die Diskretisierung der
Konvektionsterme nicht zwangsläufig auf beschränkte oder strikt stabile Verfahren limitiert [18].
3.3.3 Diskretisierung der Transportgleichung für
Die Diskretisierung der konvektiven Terme
V
· u dV =
S
u
dS =
f
S
f
· u
f
(3.12)
V
· u
r
(1 - ) dV =
S
u
r
(1 - )dS =
f
S
f
· u
r
(1 - )
f
(3.13)
bedarf besonderer Beachtung. Für die Interpolation von S
f
· u
f
in (3.12) kommt ein Van-Leer-Schema
zum Einsatz. Es handelt sich dabei um ein spezielles Upwind-Verfahren 2. Ordnung, in dem numerische
Oszillationen minimiert werden [11].
Für die Berechnung des Fluss-Terms S
f
·u
r
, f
=
r
in (3.13) gibt es viele Möglichkeiten [21]. In interFoam
wird
r
explizit aus dem Geschwindigkeitsfeld berechnet, dabei aber durch die maximale Geschwindig-
keit (bzw. den maximalen Fluss) im gesamten Rechengebiet limitiert.
Der Einfluss der künstlichen Kompression kann mit einem Parameter c
gesteuert werden. Mit c
= 0
wird dabei vollständig auf die künstliche Kompression verzichtet. c
= 1 entpricht konservativer Kom-
pression und bei c
1 wird durch höhere Kompressionsgeschwindigkeiten verstärkt komprimiert. Dabei
verbleibt die maximale Geschwindigkeit im gesamten Rechengebiet weiterhin als obere Grenze. Bei den
vorliegenden Berechnungen wurde jeweils die konservative Komprimierung (c
= 1) genutzt.
Die diskretisierten Flüsse werden daraufhin dem MULES-Schema übergeben, das die Zeitintegration
übernimmt.
3
MULES: multidimensional universal l imiter with explicit solution. Entwickelt von OpenCFD.
16
3.4 Berücksichtigung der Grenzflächenspannung
Ist die Verteilung von
im gesamten Problemgebiet bekannt, können auch Effekte durch Grenzflächen-
spannungen berücksichtig werden. In den Navier-Stokes-Gleichungen (2.11) muss die Kraft infolge der
Grenzflächenspannung
als geeignete Volumenkraft f
umgesetzt werden. Auf physikalischer Ebene
erzeugt die Grenzflächenspannung über einer Phasengrenze mit der Krümmung
einen Drucksprung
[24]. Die Grenzflächenkräfte haben dabei eine Komponente in Normalenrichtung der Grenzfläche.
Wenn die Fluide in Ruhe sind, stehen die Grenzflächenkräfte und der Drucksprung
p =
(3.14)
mechanisch im Gleichgewicht. Anderenfalls würde sich die Grenzfläche bewegen. Die Schwierigkeit be-
steht darin, diesen Drucksprung als geeigneten Druckgradienten auszudrücken, um die Einbindung in
die Impulsbilanz (2.11) möglich zu machen. Brackbill et al. [2] haben dazu das CSF-Modell
4
entwickelt.
Es lautet
f
()
(3.15)
Die Krümmung wird dabei aus der Divergenz des Normaleneinheitsvektors der Phasengrenze bestimmt
= - · n
0
i
= - ·
||
(3.16)
Die Definition der Krümmung in (3.16) bedeutet, dass für
0 das Fluid 1 auf der konkaven Seite der
Phasengrenze liegt, für
0 Fluid 2. Der Vektor · () zeigt jedoch immer in Richtung der konkaven
Seite der Phasengrenze [25] und somit auch in Richtung des Drucksprungs.
f
= () = - ·
||
(3.17)
4
CSF: continuum surface f orce
17
Diskretisierung von
f
Numerisch wird der Ausdruck für f
(3.17) wie ein Quellterm behandelt. Die Krümmung
wird dabei al-
lerdings nicht mit dem Normaleneinheitsvektor in (3.16) berechnet, sondern mithilfe eines stabilisierten
Normaleneinheitsvektor
n
0
i
=
()
|| +
.
(3.18)
Er unterscheidet sich von einem herkömmlichen Normaleneinheitsvektor durch den Skalar
. Er dient
zur Stabilisierung der Berechnung ausserhalb der Phasengrenze, wo || 0 gilt.
=
10
-8
1
N
c
N
c
c
=1
3
V
c
(3.19)
3.5 Randbedingungen
Die Berücksichtigung der unterschiedlichen Randbedingungstypen für Druck und Geschwindigkeit er-
folgt im üblichen Rahmen bei der Verwendung von Finite-Volumen-Methoden. Auf diese Randbedin-
gungstypen wird hier nicht eingegangen, es wird lediglich auf die ausführliche Literatur über Finite-
Volumen-Methoden verwiesen ([5], [11], [22]).
Interessanter ist der Randbedingungstyp für das
-Feld bei einen vorgeschriebenen Drei-Phasen- Kontakt-
winkel
, der nach der Einführung des Normaleneinheitsvektor der Phasengrenze (3.18) erklärt werden
kann. Wird eine fester Drei-Phasen-Kontaktwinkel als Randbedingung verlangt, muss der Normalenein-
heitsvektor der Phasengrenze gemäß
n
0
i
= n
w
cos
() + t
w
sin
()
(3.20)
vorgeschrieben werden. In (3.20) ist n
w
der Normaleneinheitsvektor der Wand, der in Richtung der
Wand zeigt, und t
w
der Einheitsvektor tangential zur Wand, der in Richtung der Flüssigkeit zeigt (vgl.
hierzu [25]).
3.6 Zeitschrittweitensteuerung
Die Zeitschrittweitensteuerung erfolgt in OpenFOAM adaptiv. Als Grundlage für stabile und akkurate
Berechnungen wird die Courant-Zahl für eine Zelle auf einen maximalen Wert von 0.5 limitiert.
C o
=
|u|t
h
0.5
(3.21)
Darin ist h die Zellengröße in Geschwindigkeitsrichtung [18]. Die Zeitschrittweite wird entsprechend
(3.21) zu
t
0.5h
|u|
(3.22)
gewählt.
18
3.7 Lösungsalgorithmus
Es gilt nun das inkompressible Zweiphasensystem, welches durch die oben dargestellten Gleichungen
beschrieben wird, zu lösen.
t
+ · (u) + · u
r
(1 - ) = 0
(u)
t
+ · uu = -p
d
+ · µu + (u) · µ - g · x + f
· u = 0
mit den zwischen den zwei Phasen gemittelten Größen
=
1
+ (1 - )
2
µ = µ
1
+ (1 - )µ
2
Nach der Auswahl einer geeigneten Zeitschrittweite wird zuerst die Transportgleichung für
mit dem
bereits erwähnten MULES-Schema gelöst. Die Transportgleichung wird dabei in jedem Zeitschritt mehr-
mals gelöst. Das hat den großen Vorteil, dass stabile Lösungen für
erhalten werden, ohne eine Zeit-
schrittweitenreduzierung auf Kosten der Simulationszeit inkauf nehmen zu müssen [18]. Ein Parameter
(n
), der auch vom Benutzer gesteuert werden kann, gibt an wie oft die Transportgleichung in jedem
Zeitschritt gelöst werden soll. In den folgenden Berechnungen wurde n
= 2 gewählt. Das bedeutet, dass
die Transportgleichung in jedem Zeitschritt zwei Mal jeweils mit der Hälfte der tatsächlichen Zeitschritt-
weite gelöst wird.
Aus der Lösung für
werden mit Hilfe des beschriebenen CSF-Modells die Grenzflächenkräfte berechnet
und stehen somit in den Impulserhaltungsgleichungen zur Verfügung.
Die Lösung von Impulserhaltungs- und Kontinuitätsgleichung erfolgt daraufhin mit dem PISO
5
-Algorihtmus,
der die Druck-Geschwindigkeits-Kopplung verarbeitet.
5
PISO: Pressure-Implicit Split-Operator
19
4 Numerische Testprobleme
Um die Zuverlässigkeit des zur Verfügung stehenden Lösers interFoam beurteilen zu können, werden im
Folgenden drei Testprobleme mit einer Grenzfläche zwischen zwei inkomperessiblen Fluiden gelöst und
mit experimentellen oder numerischen Referenzdaten verglichen.
4.1 Der gebrochene Damm
Der gebrochene Damm ist ein Problem zur Validierung von numerischen Berechnungsprogrammen zur
Simulation von Mehrphasenströmungen. Es wird in einer Vielzahl von numerischen Arbeiten als Testbe-
rechnung herangezogen. Da das Problem relativ einfache Randbedingungen, aber sehr große Verände-
rungen innerhalb des Rechengebietes beinhaltet, eignet es sich hervorragend, um die Leistungsfähigkeit
eines Rechenprogrammes zu untersuchen. Anhand der Ergebnisse kann eine Aussage über die korrekte
Berechnung der Grenzflächenposition und der Implementierung von Trägheitskräften getroffen werden.
Ausgangspunkt der Berechnung ist eine Wassersäule der Breite a und Höhe h, die zwischen zwei Wänden
eingeschlossen ist. Zum Zeitpunkt t
= 0 wird eine der Wände plötzlich entfernt, und die Wassersäule er-
gießt sich unter dem Einfluss der Gravitation auf einer ebenen Fläche jenseits der bis dahin vorhandenen
Wand (vgl. Abb. 4.1).
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
0000000000000000000000000000
0000000000000000000000000000
1111111111111111111111111111
1111111111111111111111111111
x
g
II
III
I
IV
a
h
x
2
x
1
B
H
Abbildung 4.1: Problemgeometrie des gebrochenen Dammes
Die Diskretisierung soll auf einem zweidimensionalen Rechengebiet durchgeführt werden. Da Open-
FOAM jedoch aussschließlich dreidimensionale Netze erzeugt, wird in die dritte Dimension auch eine
Lage Zellen erzeugt. Auf den Flächen, deren Normalen in x
3
-Richtung zeigen, müssen dann spezielle
Randbedingungen gewählt werden (Randbedingungstyp empty), damit in dieser Richtung keine Lösung
berechnet wird.
20
Die Größe des Rechengebietes (B
= 1 m, H = 0.2 m) wurde so gewählt, dass die Wasserfront am Simu-
lationsende das rechte Ende des Rechengebietes noch nicht erreicht hat. Bei der Diskretisierung kamen
dabei 50000 Kontrollvolumen zum Einsatz. Da zu erwarten ist, dass sich das Wasservolumen und somit
die Grenzfläche für vorgeschrittene Simulationszeiten vorwiegend im unteren Bereich der Problemgeo-
metrie befinden wird, wurden die Zellen in negativer x
2
-Richtung verdichtet. Dadurch ergibt sich im
relevanten Berechnungsbereich eine höhere Zellendichte.
Entlang der x
2
-Achse wurde dabei ein Expansionverhältnis von 10 gewählt. Das bedeutet, dass die letzte
Zelle entlang des Rechengebietsrandes (Bereich I) 10mal höher ist als die erste Zelle. In x
1
-Richtung
wurde auf eine Verdichtung verzichtet.
Abbildung 4.2: Numerisches Gitter für den gebrochenen Damm
Die Simulation wurde mit einem Wasservolumen der Breite a
= 0.0572 m und der Höhe h = 0.1144 m
gestartet. Das entspricht einem Höhen-Breiten-Verhältnis von h
/a = 2. Durch die verschmierte Phasen-
grenze entspricht das initialisierte Wasservolumen allerding nicht exakt dem erwarteten Volumen bei
diesem Höhen-Breiten-Verhältnis. Bei dieser Diskretisierungsstufe liegt der Initialisierungsfehler aller-
dings nur bei 1.3 %. Als Stoffwerte werden die Werte für Wasser und Luft bei 20
o
C
gewählt.
[kg/m
3
] [10
-6
m
2
/s] [N/m]
Luft
1.188
15.35
-
Wasser
998.21
1.004
0.07
Tabelle 4.1: Stoffwerte für den gebrochenen Damm [26]
Als Randbedingungen wurden die Geschwindigkeitskomponenten an den beiden Wänden (Bereiche I
und IV) zu Null gesetzt. Auf den beiden anderen Seiten des Problemgebiets (Bereiche II und III) wurden
Ausstromrandbedingungen gewählt. Das Druckniveau wird im Bereich II mit p
= 0 festgelegt. Der Drei-
Phasenkontaktwinkel zwischen Wand, Wasser und Luft wurde im Bereich I und IV mit 90
o
gewählt.
21
Bereich
u
p
I
u
= 0
p/ x
1
= 0 / x
1
= 0
II
u
1
= 0, u
2
/ x
2
= 0
p
= 0
/ x
2
= 0
III
u
2
= 0, u
1
/ x
1
= 0 p/ x
1
= 0 / x
1
= 0
IV
u
= 0
p/ x
2
= 0 / x
2
= 0
Tabelle 4.2: Randbedingungen für den gebrochenen Damm
Martin Moyce [20] führten mit den oben genannten Werten und Stoffen Experimente durch. Die
Beurteilung der numerischen Ergebnisse soll deshalb anhand der charakteristischen experimentellen
Werte erfolgen. Diese sind die Position der Wasserfront x und die Höhe der restlichen Flüssigkeit
(vgl.
Abb. 4.1). Martin Moyce überführten diese Messwerte in dimensionslose Größen und fassten sie in
Tabellen zusammen. Als dimensionslose Größen wurden dabei gewählt:
X
=
x
a
H
=
2a
T
=
2g
a
t
=
g
a
t
X
und H sind dabei jeweils dimensionslose Längen, die der Position der Wasserfront bzw. der Wasserhöhe
entsprechen. T und
sind dagegen dimensionslose Zeiten. T gehört zu den Messwerten der Wasserfront,
zu den Messwerten der Wasserhöhe. Die Simulationsergebnisse für die Position der Wasserfront und
die Höhe der Wassersäule wurden auf diese dimensionslosen Kenngrößen reduziert, um den Vergleich
mit dem Experiment zu ermöglichen.
Als Wasseroberfläche wurde dabei die Isolinie mit
= 0.5 interpretiert. In Abbildung 4.3 ist diese Kontur
durch die schwarze Linie visualisiert.
Abbildung 4.3:
-Feld für t = 0.15s
22
Den Vergleich der experimentellen Daten mit den Simulationsergebnissen zeigen Abb. 4.4 und Abb. 4.5.
0
1
2
3
4
5
6
7
8
9
10
0
5
10
15
T
X
Martin Moyce
Simulation
Abbildung 4.4: Zeitlicher Verlauf der dimensionslosen Position der Wasserfront
0
1
2
3
4
5
6
7
8
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
H
Martin Moyce
Simulation
Abbildung 4.5: Zeitlicher Verlauf der dimensionslosen Position der Wasserhöhe
Etwaige Unterschiede der Absolutwerte lassen sich durchaus darauf zurückführen, dass in [20] der An-
fangszeitpunkt des Experimentes nicht genau bestimmt werden konnte. Die Zeitdifferenz zwischen den
Messpunkten ist also genauer als die absoluten Zeitangaben. Die Übereinstimmung der Daten ist - auch
im Hinblick auf diese Tatsache - ausgesprochen gut.
23
4.2 Einfluss der Grenzflächenspannung
Mit dem folgenden Testfall soll die Implementierung der Obeflächenspannung in OpenFOAM überprüft
werden. Dazu wird ein Fluidvolumen betrachtet, welches sich zu Beginn der Simulation nicht in sei-
nem durch die Grenzflächenkräfte bestimmten Gleichgewichtslage befindet. Für fortgeschrittene Berech-
nungszeiten sollte das Fluidvolumen dann eine Gleichgewichtslage einnehmen.
In der Simulation wird ein anfänglich würfelförmiges Flüssigkeitsvolumen mit konstanter Grenzflächen-
spannung in einem gasförmigen Umgebungsmedium unter Vernachlässigung der Gravitation als Anfangs-
konfiguration gewählt (Abb. 4.6(a)).
Die stabile Form eines solchen Flüssigkeitstropfens ist unter Annahme verschwindender Gravitation die
Kugelgestalt (Abb. 4.6(b)). Das anfänglich ,,gestörte" Flüssigkeitsvolumen sollte sich nach einer endli-
chen Anzahl an Zeitschritten unter dem Einfluss der Grenzflächenspannung in dieser stabilen Gleichge-
wichtslage befinden. Dabei kann es durchaus zu Schwingungen um die Gleichgewichtslage kommen, die
aber aufgrund der inneren Reibung der Flüssigkeit sukzessive abklingen sollten.
Betrachet man das Problem in einer zweidimensionalen Ebene, ist die stabile Gleichgewichtslage keine
Kugelform, sondern eine Kreisform mit konstantem Radius R. Aufgrund der gekrümmten Phasengrenze
kommt es zu einer Erhöhung des Druckes innerhalb des Tropfens gemäß
p =
R
(4.1)
Als relevante Bewertungskriterien der Berechnung sollen sowohl die Form der stabilen Gleichgewichts-
lage als auch der konstant erhöhte Druck innerhalb der Flüssigkeit dienen.
B
L
2
1
(a) ,,Gestörtes" Fluidvolumen
L
2
1
R
(b) Gleichgewicht
Abbildung 4.6: Problemgeometrie
Abbildung 4.6 zeigt eine schematische Darstellung des zweidimensionalen Rechengebietes sowohl für
das anfänglich gestörte Fluidvolumen und das Fluidvolumen unter Gleichgewichtsbedingungen. Das
quadratische bzw. kreisförmige Flüssigkeitsgebiet (Bereich 1) mit der Kantenlänge B bzw. dem Radi-
us R ist dabei in einem quadratischen Rechengebiet der Kantenlänge L zentriert und von einem anderen
Fluid (Bereich 2) umgeben.
24
Die Simulation eines solchen Testfalles hat sich in OpenFOAM als unerwartet schwierig dargestellt. Zwar
existieren viele numerische Ergebnisse unterschiedlicher Autoren aus vergleichbaren Berechnungen im
Sinne von Abb. 4.6, die jedoch alle nicht in ihrer von den Autoren vorgeschlagenen Weise in OpenFOAM
nachgerechnet werden konnten. Hier soll eine Zusammenfassung der durchgeführten Berechnungen fol-
gen. Die dabei aufgetretenen Probleme werden dabei dargestellt und wenn möglich erklärt. Letztendlich
wird gezeigt unter welchen Vorraussetzungen eine halbwegs befriedigende Lösung für das Problem ge-
funden werden konnte.
Bevor jedoch auf die einzelnen Berechnungen eingegangen wird, soll ein Phänomen näher erläutert wer-
den, welches im Zusammenhang mit der Berechnung von Grenzflächenspannungen in VOF-Methoden
beobachtet wird: die sogenannten parasitären Strömungen.
4.2.1 Parasitäre Strömungen
Als parasitäre Strömungen werden unphysikalische Geschwindigkeitsfelder im Bereich der Phasengrenze
bezeichnet. Lafaurie et. al. [14] beobachteten bei einem unter den oben genannten Bedingungen ma-
kroskopisch in Ruhe befindlichen Tropfen ein Geschwindigkeitsfeld im Bereich der Phasengrenze und
führten als Erste den Begriff von ,,parasitären" Strömungen ein.
Abbildung 4.7: Parasitäres Geschwindigkeitsfeld trotz makroskopischem Gleichgewicht [14]
Weiterhin beobachteten sie, dass die maximale Geschwindigkeit der parasitären Strömung bei konstan-
tem Radius von den verwendeten Stoffwerten abhängt
u
ma x
= K
(4.2)
Die maximale Geschwindigkeit wächst nach (4.2) mit der Grenzflächenspannung und kann bei Simula-
tionen mit dominanten Grenzflächenkräften die Berechnung des Geschwindigkeitsfeldes negativ beein-
flussen, oder sogar die Auflösung der Phasengrenze zur Folge haben [9]. In (4.2) ist K eine Konstante,
die von Vincent Caltagirone [27] als Qualitätsmerkmal des numerischen Modells zur Modellierung der
Grenzflächenkräfte interpretiert wird. Je kleiner die Dimension dieser Konstante, desto kleiner skaliert
sich auch das parasitäre Strömungsfeld, und desto ,,besser" ist das zugrundeliegende Modell. Während
bei Lafaurie et. al. [14] vor 15 Jahren Konstanten im Bereich von K 10
-2
beobachtet wurden, wird
in der aktuelleren Veröffentlichung von Vincent Caltagirone [27] aus 2004 ein kleinerer Bereich von
10
-3
K 10
-10
angegeben.
25
Die Ursachen des parasitären Strömungsfeldes sind auf Fehler in der Berechnung des Normalenvektors
der Grenzfläche (3.18) bzw. der Krümmung (3.16) zurückzuführen. Betrachtet man einen Tropfen in
seiner kreisrunden Gleichgewichtsform unter Verwendung einer ,,idealen" VOF-Methode, wäre der Nor-
malenvektor der Grenzfläche ein Vektor, der in radialer Richter in das Innere des Tropfens zeigen würde.
Die Krümmung wäre mit
= 1/R konstant, und die Grenzflächenkraft (3.17) wären somit ein Vektor in
radialer Richtung.
Diese Idealvorstellung kann in realen VOF-Methoden jedoch nicht umgesetzt werden. Der Normalenvek-
tor der Grenzfläche besitzt immer auch ein Komponente in Umfangsrichtung und die Krümmung variiert
sowohl in radialer als auch in Umfangsrichtung geringfügig. Die Grenzflächenkraft hat demnach auch
eine fehlerhafte ,,rotatorische" Umfangsrichtung [9].
In der Impulsbilanz (2.11) wären unter idealen Bedingungen alle geschwindigkeitsabhängigen Kompo-
nenten gleich Null, und die Grenzflächenkräfte wären im Gleichgewicht mit dem Druckgradienten. Unter
realen Bedingungen kann die rotatorische Komponente der Grenzflächenkraft jedoch nicht vom Druck-
gradient kompensiert werden. Vielmehr werden zur Erfüllung der Gleichung rotatorische Geschwindig-
keitskomponenten - die parasitären Strömungen - erzeugt. Neben den rotatorischen Geschwindigkeits-
komponenten erzeugen reale VOF-Methoden auch nicht-rotatorische fehlerhafte Geschwindigkeitskom-
ponenten. Sie tragen in der Regel nicht zu den parasitären Strömungen bei, sondern verfälschen das
Ergebnis des Druckfeldes [9].
Die im folgenden vorgestellten numerischen Testfälle dienen dazu, die Fehler durch parasitäre Strömun-
gen des in interFoam verwendeten CSF-Modells näher zu beurteilen.
26
4.2.2 Durchgeführte Simulationen
Reibungsfreie Fluide
Wie im Kapitel 3.4 ausgeführt wurde, wird in interFoam die Kraft infolge der Grenzflächenspannung
mit dem CSF-Modell von Brackbill et al. [2] berechnet. Zusammen mit dem Modell in [2] wird auch eine
vollständige Testfallbeschreibung angegeben, wie sie auch hier zur Anwendung kommen soll.
Dabei wird der Tropfen bereits in seiner Gleichgewichtslage in Kreisform (R
= 2 cm, Abb. 4.6(b)) auf
einem Rechengebiet (L
= 6 cm) mit unterschiedlichen Diskretisierungsstufen (30x30, 60x60 Kontroll-
volumen) zentriert initialsiert. Zur Anwendung kommen dabei zwei reibungsfreie Fluide (
1
=
2
= 0)
mit unterschiedlichen Dichten, wobei der Tropfen die doppelte Dichte wie die umgebende Flüssigkeit
aufweist (
1
= 2
2
= 1000 kg/m
3
). Als konstante Grenzflächenspannung wird
= 23.61 · 10
-3
kg
/s
2
angenommen. Das entspricht der Grenzflächenspannung, wie man sie zwischen Ethanol und Wasser
vorfindet. Die Dichten und Viskositäten entsprechen allerdings nicht diesem Stoffpaar. Über Randbedin-
gungen wurden in [2] keine Angaben gemacht. Es wurden daher Neumann-Randbedingungen für alle
Variablen auf allen Rechengebietsrändern gewählt.
Fluid
[kg/m
3
] [m
2
/s]
[kg/s
2
]
1
1000
0
23.61
· 10
-3
2
500
0
-
Tabelle 4.3: Verwendete Stoffwerte bei der Simulation reibungsfreier Fluide [2]
Den konstant erhöhten Druck innerhalb des Tropfens würde man unter den oben genannten Vorrausset-
zungen nach (4.1) mit
p =
R
= 11.805
N
m
2
(4.3)
erwarten.
Obwohl der Tropfen zu Beginn der Simulation bereits in seiner Gleichgewichtslage initialisiert wurde,
baut sich bereits mit den ersten Zeitschritten eine Schwingung der Phasengrenze um diese Gleichge-
wichtslage auf. Man spricht dabei auch von Kapillarwellen. Die Schwingungsamplituden klingen im Ver-
lauf der Simulation jedoch nicht ab. Es kann davon ausgegangen werden, dass aufgrund der reibungs-
freien Flüssigkeiten auch keine dämpfenden Kräfte auf diese Schwingung wirken. Ab einem gewissen
Zeitpunkt wurde beobachtet, dass das Geschwindigkeitsfeld des schwingenden Tropfens asymmetrisch
wird und er von seiner zentrierten Position im Rechengebiet abdriftet und sogar das Rechengebiet ver-
lässt. In diesem Fall ist ein Vergleich von Drucksprung (4.3) oder statischer Form unmöglich.
Gerlach et al. [8] haben den Testfall von Brackbill et al. um den realitätsnäheren Fall eines höheren Dich-
teverhältnisses (
1
/
2
= 1000), wie man es bei einem Alkohol-Luft-Stoffpaar erwarten würde, erweitert.
Alle anderen Parameter blieben unverändert. In diesem Fall verhielt sich die Simulation allerdings noch
instabiler. Das Aufbauen der Schwingung und das Abdriften des Tropfens erfolgte schneller.
27
Um dem Abdriften des Tropfens entgegenzuwirken, wurde versucht die Symmetrie des Problems auszu-
nutzen. Dafür wurde nur das obere rechte Viertel der Problemgeometrie in Abb. 4.6 diskretisiert. Durch
diese Maßnahme konnte das Abdriften des Tropfens zwar verhindert werden, aber es wurde beobach-
tet, dass die Kapillarwellen auch für sehr hohe Simulationszeiten (t
1000 s) nicht abklingen. Somit
war auch unter diesen Vorraussetzungen kein Vergleich des Drucksprunges bzw. der statischen Form des
Tropfens möglich.
Die Instabilität dieser Simulationen ist offensichlich auf die in Kap. 4.2.1 beschriebenen parasitären Strö-
mungen zurückzuführen. Nach Gleichung (4.2) würden für reibungsfreie Fluide (
0) unendlich hohe
parasitäre Strömungsfelder erwartet werden. Es ist also nicht verwunderlich, dass die Phasengrenze im
oben beschriebenem Maße destabilisert wird.
Wie oben beschrieben, können reibungsfreie Fluide wie bei Brackbill et al. oder Gerlach et al. nicht er-
folgreich simuliert werden. Sowohl physikalisch als auch numerisch ist die Annahme eines reibungsfreien
Fluides generell fragwürdig. Zwar werden für den Einsatz als Wärmeträger (z.B. in Wärmerohren) Fluide
mit sehr niedrigen Viskositäten bevorzugt [24], numerisch besteht jedoch zwischen einem reibungsfrei-
em und einem Fluid mit geringer Viskosität ein bedeutender Unterschied hinsichtlich der Stabilit der
Simulationen, wie sich im folgenden zeigen wird.
Reales Alkohol-Luft-Stoffpaar
Vincent Caltagirone [27] verwenden die praxisrelevantere Stoffwertkombination eines Ethanoltrop-
fens in Luft. Stoffwerte in dieser Größenordnung werden eher in Anwendungsbereichen von Phasen-
wechselprozessen angetroffen.
Fluid
[kg/m
3
] [10
-6
m
2
/s]
[kg/s
2
]
1
797.88
1.5040
23.61
· 10
-3
2
1.1768
8.4967
-
Tabelle 4.4: Verwendete Stoffwerte bei der Simulation mit realen Stoffwerten [27]
Auf einem Rechengebiet mit L
= 75 mm und 120x120 Kontrollvolumen wird ein rechteckig zentrierten
Ethanolvolumen mit B
= 40 mm initialisert. Die Diskretisierung wurde dabei so durchgeführt, dass die
Grenzfläche des Ethanolvolumens zu Beginn der Berechnung mit Elementseiten zusammenfällt. Damit
liegt zu Beginn der Simulation nur eine minimale Verschmierung der Phasengrenze vor. Wie zu erwarten
war, bilden zu sich Beginn der Simulation Kapillarwellen um die Gleichgewichtslage aus, die jedoch
im Gegensatz zur Simulation nach Brackbill et al. [2] langsam abklingen. Dies ist offensichtlich auf
das Vorhandensein der Viskosität zurückzuführen. Trotz merklich abnehmender Schwingungsamplituden
wird das Geschwindigkeitsfeld mit fortlaufender Simuationszeit auch in dieser Simulation asymmetrisch.
Der Tropfen beginnt auch diesmal seine zentrierte Position im Rechengebiet zu verlassen und driftet bis
zum Kontakt mit dem Rechengebietsrand ab. Ein Vergleich von Drucksprung oder statischer Form ist
somit auch diesmal nicht möglich.
28
Alkohol-Luft-Stoffpaar mit künstlich erhöhten Viskositäten
Da möglicherweise ein starkes parasitäres Strömungsfeld durch die niedrige Viskosität induziert wird
1
,
wurde versucht, mit künstlich erhöhten Viskositäten ein stabiles Ergebnis zu erhalten. Ab einer Erhöhung
der Viskositäten um den Faktor 10000 auf
1
= 1.5040 · 10
-2
m
2
/s und
2
= 8.4967 · 10
-2
m
2
/s wurden
stabile Lösungen erzielt. Die weiteren Stoffwerte blieben unverändert. Mit diesen künstlich erhöhten
Zähigkeitseigenschaften konnte das Abdriften der Blase aus der zentrierten Position im Rechengebiet
verhindert werden, sodass Vergleiche für Kontur und Drucksprung möglich waren.
Fluid
[kg/m
3
] [10
-2
m
2
/s]
[kg/s
2
]
1
797.88
1.5040
23.61
· 10
-3
2
1.1768
8.4967
-
Tabelle 4.5: Verwendete Stoffwerte bei der Simulation mit künstlich erhöhten Stoffwerten
Initialisiert man das
-Feld mit der oben genannten Diskretisierung und den Abmessungen für L und
B
befinden sich 28.4444 % Ethanol im Rechengebiet. Das entspricht in 2D einer Fläche von A
E
= 1.6 ·
10
-3
m
2
. Aufgrund der feinen Diskretisierung entspricht das auch der Fläche eines Rechteckts mit der
Kantenlänge B. Der Initialisierungsfehler wie er im Fall des gebrochenen Dammes erwähnt wurde, ist im
vorliegenden Fall also vernachlässigbar klein. In der Gleichgewichtslage würde man dann entsprechend
einen Kreis mit dem Radius
R
E
=
A
E
= 0.02256 m
(4.4)
erwarten, sodass der theoretische Drucksprung im Tropfen
p
E
=
R
E
= 1.0462
N
m
2
(4.5)
beträgt. Wegen der erhöhten Viskosität muss man jedoch in Kauf nehmen, dass physikalische Diffusions-
vorgänge etwas langsamer ablaufen als in der Realität. Das zieht im vorliegenden Fall nach sich, dass
das Erreichen der stabilen Gleichgewichtsform der Blase und die Druckerhöhung etwas längere Simula-
tionszeiten benötigt, um konstante Werte zu erreichen.
Bei der Auswertung der numerischen Ergebnisse, werden die Drücke in Zellen mit
0.99·
1
summiert
und durch die Anzahl an betrachteten Zellen N
c
geteilt. Somit ergibt sich ein mittlerer Druck
p
=
1
N
c
N
c
i
p
i
(4.6)
Zusätzlich wird der Druck in der Mitte des Rechengebiets - also auch in der Mitte des Tropfens - durch
Auswertung des Druckes in einer der vier Zellen, die sich um den Rechengebietsmittelpunkt verteilen
p
M
p(x = y = 37.5 mm)
(4.7)
protokolliert.
1
vgl. Gleichung (4.2)
29
Ende der Leseprobe aus 176 Seiten
- nach oben
Details
- Arbeit zitieren
- , 2009, Numerische Simulation der periodischen Blasenbildung an einer Einströmöffnung mit dem Programm OpenFOAM, Hamburg, Bedey Media GmbH, https://www.diplom.de/document/226629
Ähnliche Arbeiten
Deutschlands größter Sharing Community
Jetzt Zusammenfassungen, Skripte und Klausuren kostenlos downloaden!
Uniturm.de
Ihre Vorteile als Autor
Geld verdienen Für jedes verkaufte Exemplar erhalten Sie Autorenhonorar - bis zu 40%.Kostenlose Buchveröffentlichung
Als eBook-Studie im Original und / oder auf Wunsch zusätzlich als hochwertiges Fachbuch in unserem renommierten Buchverlag.
Renommee als Fachbuchautor
Präsentieren Sie sich als Fachfrau oder Fachmann in Ihrem Fachgebiet und machen Sie sich bekannt.
Persönliche Betreuung
Unser Lektorat wird Sie persönlich betreuen und ist auch telefonisch unter 040-655992-14 erreichbar.
Weltweit im Buchhandel
Ihr Buch ist weltweit im Buchhandel und Online-Buchhandel wie z.B. amazon erhältlich.
Jetzt Autor werden!
Kommentare