Numerische Simulation der periodischen Blasenbildung an einer Einströmöffnung mit dem Programm OpenFOAM

Gatzka, Sebastian

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

von Sebastian Gatzka (Autor:in)

Ingenieurwissenschaften - Maschinenbau

Zusammenfassung

Inhaltsangabe:Einleitung:
Die hohe Nachfrage nach leistungsfähigen Wärmeübertragern für die Kühlung von elektronischen Bauteilen, zur Effizienzsteigerung von Kreisprozessen für Wärmekraft- und Kälteanlagen, sowie für die chemische 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 (Kondensation) 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.
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 Simulation mit dem CFD-Berechnungsprogramm OpenFOAM Version 1.4.1 erstellt. Die zentrale Vereinfachung, die das Gesamtproblem in ein geeignetes Modell überführt, ist das Tatsache, dass der Verdampfungsprozess 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 einströ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 numerische
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 […]

Leseprobe

Inhaltsverzeichnis

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.

http://www.diplomica.de, Hamburg 2009

Inhaltsverzeichnis

Nomenklatur

1 Einleitung

2 Physikalische Grundgleichungen

2.1 Massenerhaltung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2 Impulserhaltung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3 Drei-Phasen-Kontaktwinkel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 Numerische Verfahren

3.1 Numerisches Gitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2 Finite-Volumen-Verfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

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

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

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

Literaturverzeichnis

Anhang

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

Breite des Wasservolumens

[m]

Breite des Rechengebietes

[m]

Bond-Zahl, Eötvös-Zahl

[-]

C o

Courant-Zahl

[-]

Durchmesser

[m]

äquivalenter Durchmesser

[m]

Kraft infolge der Grenzflächenspannung

Frequenz

Gravitationskonstante

Höhe des Rechengebietes

[m]

dimensionslose Position der Wasserhöhe

[-]

Höhe des Wasservolumens

[m]

Zellengröße, äquidistante Zellengröße

[m]

Einheitsmatrix

[-]

dimensionslose Konstante,

Qualitätsmerkmal zur Beurteilung des CSF-Modells

[-]

M o

Morton-Zahl

[-]

Normalenvektor

[-]

Stabilisierter Normalenvektor

[-]

Zeitschritt

[-]

Anzahl an Unterzyklen des MULES-Schemas

[-]

Anzahl

[-]

mittlerer Druck

Druck

modifizierter Druck

Einströmrate

min

Radius

[m]

Blendenradius

[m]

äquivalenter Radius

[m]

Reynolds-Zahl

[-]

Normalenvektor der Zellenseite

Oberfläche

Standardabweichung, empirische Varianz

hier:

Zeit

[s]

Blasenbildungszeit

[s]

Temperatur

[K]

Geschwindigkeit

dimensionslose Geschwindigkeit

[-]

Betrag des Geschwindigkeitsfeldes

künstliche Kompressiongeschwindigkeit

ma x

maximale Geschwindigkeit der parasitären Strömungen

Volumen

Tangentialvektor

[-]

dimensionslose Position der Wasserfront

[-]

Position der Drei-Phasen-Konstaktlinie beim Simulationsstart

[m]

III

Griechische Buchstaben

Drei-Phasen-Kontaktwinkel

[

]

Stabilisierungsparameter

Übergangsbereich

[m]

Veränderung der Blasenbildungszeiten

[%]

Veränderung der Blasenvolumen

[%]

VOF-Variable

[-]

Krümmung

dynamische Viskosität

kinematische Viskosität

volumetrischer Fluss der künstlichen

Kompressionsgeschwindigkeit ?

Drei-Phasen-Kontaktwinkel

[

]

Dichte

Spannungstensor

Grenzflächenspannung

Reibungsspannungstensor

dimensionslose Geschwindigkeit bezogen auf U

[-]

dimensionslose Zeit bezogen auf H

[-]

dimensionslose Zeit bezogen X

[-]

Indizes

Blase

numerische Zelle

Ethanol

Mittelpunkt der Zellenseite

Ergebnisse aus dem Modell von Gerlach et. al. [7]

Flüssigkeit

Mittelpunkt

Mittelpunkt der Nachbarzelle

Mittelpunkt der numerischen Zelle

Festkörper

sim

Simulation

Dampf

Wasser

Rechengebietsrand

Grenzfläche

Simulation mit Einlasskanal

sat

Sättigungsbedingungen

Wand

-Richtung

Fluid 1

-Richtung

Fluid 2

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

Einheitsvektor

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

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

- 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]

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]

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

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.

OpenFOAM: Open Field Operation And Manipulation.

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.

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:

+ · 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)

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)

+ · 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)

+ · uu = -p + · + g + f

(2.5)

Weiterhin gilt für inkompressible Newtonsche Fluide der Reibungsspannungstensor [23]

= µ (u) + (u)

(2.6)

Bindet man das Materialgesetz (2.6) in Gleichung (2.5) ein, erhält man die die Navier-Stokes-Gleichungen

(u)

+ · uu = -p + · µ u + (u)

+ 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)

= · µu + (u) · µ

(2.8)

Aus der Differenz des Druckes p und des hydrostatischen Druckes

g · x wird ein modifizierter Druck

= p - g · x

(2.9)

bzw. dessen Gradient

= p - g - g · x

(2.10)

in (2.7) eingeführt. Die implementierten Navier-Stokes-Gleichungen lauten damit

(u)

+ · uu = -p

+ · µu + (u) · µ - g · x + f

(2.11)

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

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

() =

(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

und die

Grenzflächenspannung zwischen Festkörper und Gas

auf. Die Kontaktwinkel zwischen drei Stoffpaar-

kombinationen werden experimentell bestimmt.

Bei

= 90

steht die Phasengrenze senkrecht auf der Wand. Für

benetzt die Flüssigkeit die

Wand.

Abbildung 2.2: Beispiele für verschiedene Kontaktwinkel (

= ) [13]

Die intermolekularen Bindekräfte innerhalb der Flüssigkeit.

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]

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.

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.

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.

d V

· p

, 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

(vgl. Abb. 3.1) ergibt sich damit ein Fluss der Variable (hier p

)

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

· µu dV =

· µu =

· (u)

(3.2)

muss sowohl für die dynamische Viskosität

an der Kontrollvolumenseite also auch für die Interpo-

lation des Gradienten normal zur Kontollvolumenseite S

· (u)

ein Interpolationsverfahren gewählt

werden.

Für

wird dabei das oben erwähnte Zentraldifferenzenverfahren verwendet. Für S

· (u)

wird auch

ein Zentraldifferenzenverfahren, aber mit einer Korrekturmöglichkeit im Fall nicht-orthogonaler Gitter

gewählt.

Es sei an dieser Stelle darauf hingewiesen, dass zur Unterscheidung der Variablen in der Zellenmitte

und an den Zellen-

seitenmitten

unterschiedliche Indizes verwendet werden.

Konvektionsterm

· uu dV =

· uu =

· uu

(3.3)

Für die Interpolation von

wird ein lineares TVD-Schema

verwendet. Es dient zur Sicherung der

Beschränktheit und der Vermeidung von ungewollten Oszillationen in der Lösung [11].

Zeitableitung

udV =

- uV

-1

(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

gdV =

· V

(3.5)

Quellterme der Form (3.5) werden durch die Mittelpunktsregel auf den Wert im Mittelpunkt des Kon-

trollvolumens approximiert.

TVD = T otal variation d iminishing.

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.

3.3.1 Grundlagen

In einer numerischen Zelle wird der Volumenanteil, den die Phase 1 bzw. 2 einnimmt, durch die Variable

beschrieben:

= 1 -

(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

(a) Grenzfläche mit

-Konturlinien

0.5

0.25

0.75

(b)

entlang der Strecke ab

Abbildung 3.4: Schematische Darstellung des Übergangsbereiches der Phasengrenze.

Nach einer Abbildung in [25].

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 - )

(3.7)

µ = µ

+ (1 - )µ

(3.8)

und die Gleichung mit diesen Zustandswerten gelöst.

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:

+ · (u) = 0

(3.9)

In OpenFOAM wird (3.9) jedoch durch einen weiteren Term zur künstlichen Komprimierung der Grenz-

fläche erweitert.

+ · (u) + · u

(1 - ) = 0

(3.10)

Darin stellt u

ein geeignetes Geschwindigkeitsfeld zur Komprimierung der Phasengrenze dar, um deren

Verschmierung zu verringern. Durch

(1 - ) wird sichergestellt, dass u

nur im Bereich der Phasengren-

ze aktiv ist und das Strömungsfeld ausserhalb dieser Region nicht beeinflusst wird ([10], [21]).

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

(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

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

· u dV =

dS =

· u

(3.12)

· u

(1 - ) dV =

(1 - )dS =

· u

(1 - )

(3.13)

bedarf besonderer Beachtung. Für die Interpolation von S

· u

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

·u

, f

in (3.13) gibt es viele Möglichkeiten [21]. In interFoam

wird

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.

MULES: multidimensional universal l imiter with explicit solution. Entwickelt von OpenCFD.

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

entwickelt.

Es lautet

()

(3.15)

Die Krümmung wird dabei aus der Divergenz des Normaleneinheitsvektors der Phasengrenze bestimmt

= - · n

= - ·

(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.

= () = - ·

(3.17)

CSF: continuum surface f orce

Diskretisierung von

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

()

|| +

(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.

-8

(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

cos

() + t

sin

()

(3.20)

vorgeschrieben werden. In (3.20) ist n

der Normaleneinheitsvektor der Wand, der in Richtung der

Wand zeigt, und t

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

0.5

(3.21)

Darin ist h die Zellengröße in Geschwindigkeitsrichtung [18]. Die Zeitschrittweite wird entsprechend

(3.21) zu

0.5h

|u|

(3.22)

gewählt.

3.7 Lösungsalgorithmus

Es gilt nun das inkompressible Zweiphasensystem, welches durch die oben dargestellten Gleichungen

beschrieben wird, zu lösen.

+ · (u) + · u

(1 - ) = 0

(u)

+ · uu = -p

+ · µu + (u) · µ - g · x + f

· u = 0

mit den zwischen den zwei Phasen gemittelten Größen

+ (1 - )

µ = µ

+ (1 - )µ

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

), 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

-Algorihtmus,

der die Druck-Geschwindigkeits-Kopplung verarbeitet.

PISO: Pressure-Implicit Split-Operator

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).

0000000000000000000000000000

1111111111111111111111111111

III

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

-Richtung zeigen, müssen dann spezielle

Randbedingungen gewählt werden (Randbedingungstyp empty), damit in dieser Richtung keine Lösung

berechnet wird.

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

-Richtung verdichtet. Dadurch ergibt sich im

relevanten Berechnungsbereich eine höhere Zellendichte.

Entlang der x

-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

-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

gewählt.

[kg/m

] [10

-6

/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

gewählt.

Bereich

= 0

p/ x

= 0 / x

= 0

= 0, u

/ x

= 0

/ x

= 0

III

= 0, u

/ x

= 0 p/ x

= 0 / x

= 0

p/ x

= 0 / x

= 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:

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

Den Vergleich der experimentellen Daten mit den Simulationsergebnissen zeigen Abb. 4.4 und Abb. 4.5.

Martin Moyce

Simulation

Abbildung 4.4: Zeitlicher Verlauf der dimensionslosen Position der Wasserfront

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

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.

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 =

(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.

(a) ,,Gestörtes" Fluidvolumen

(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.

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

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

-3

K 10

-10

angegeben.

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.

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 (

= 0)

mit unterschiedlichen Dichten, wobei der Tropfen die doppelte Dichte wie die umgebende Flüssigkeit

aufweist (

= 2

= 1000 kg/m

). Als konstante Grenzflächenspannung wird

= 23.61 · 10

-3

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

] [m

/s]

[kg/s

]

1000

23.61

· 10

-3

500

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 =

= 11.805

(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 (

= 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.

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

] [10

-6

/s]

[kg/s

]

797.88

1.5040

23.61

· 10

-3

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.

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

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.5040 · 10

-2

/s und

= 8.4967 · 10

-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

] [10

-2

/s]

[kg/s

]

797.88

1.5040

23.61

· 10

-3

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

befinden sich 28.4444 % Ethanol im Rechengebiet. Das entspricht in 2D einer Fläche von A

= 1.6 ·

-3

. 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

= 0.02256 m

(4.4)

erwarten, sodass der theoretische Drucksprung im Tropfen

= 1.0462

(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·

summiert

und durch die Anzahl an betrachteten Zellen N

geteilt. Somit ergibt sich ein mittlerer Druck

(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(x = y = 37.5 mm)

(4.7)

protokolliert.

vgl. Gleichung (4.2)

Details

Seiten

Erscheinungsform

Originalausgabe

Jahr

2009

ISBN (eBook)

9783836627092

Dateigröße

4.8 MB

Sprache

Deutsch

Institution / Hochschule

Technische Universität Darmstadt – Maschinenbau

Erscheinungsdatum

2014 (April)

Note

1,0

Schlagworte

openfoam numerik blasenbildung thermodynamik

Autor

Sebastian Gatzka (Autor:in)