Lehrbrief und Übungen

„Geoströmungstechnik“

Numerische Modellierung

von Strömungs- und Transportprozessen

im geologisch geschichteten Untergrund TU Bergakademie Freiberg Institut für Bohrtechnik und Fluidbergbau Prof. Dr. St. Wagner DM S. Boy Dipl.-Ing. A. Matwijenko Dipl.-Phys. Th. Wagner Dipl.-Geoinformatikerin E. Bennewitz 07 / 2006

1

Lehrbrief "Geoströmungstechnik/Praktikumsanleitung"

für das Lehrgebiet "Geoströmungstechnik" am Institut für Bohrtechnik und Fluidbergbau der Fakultät Geowissenschaften, Geotechnik und Bergbau

an der TU Bergakademie Freiberg

Wagner, St. / Boy, S. / Matwijenko, A. / Wagner, Th. & Studenten

Vorlesung und Übung

Numerische Modellierung von Strömungs- und Transportprozessen im geologisch geschichteten Untergrund

0. Einleitung und Aufgabenstellung Strömung/Diffusion - Transport -Parameteridentifikation (Filterströmung) Symbolverzeichnis / Indizes 1. Lösung der Strömungs- und Transportgleichung 1.1. Grundgleichungen 1.2. Transportmechanismen und ihre mathematische Beschreibung 1.2.1 Allgemeine Grundlagen, Modelle und Modellierungsaufgabe 1.2.2 Transporterscheinungen, Transportprozesse 1.2.3 Mehrphasen-Strömung und Stofftransport 1.3. Analytische Lösungen Laplacetransformation - Rücktransformation mittels Korrespondenzen - Numerische Rücktransformation (Stehfest Algorithmus) 2. Numerische Lösung von partiellen Differentialgleichungen 2.1 Numerische Modellierung von Strömungs- und Transportprozessen - Diskretisierung durch endliche Differenzen oder finite Elemente - Numerische Lösung von Differenzengleichungen - Explizite und implizite Lösungsverfahren 2.2 Numerische Verfahren - Einfache und erweiterte Bilanzmethode (CVM) - Finite - Element- Methode - Methode der Charakteristiken (MOC) - Random-Walk-Methode (RWM) - Randintegralgleichungsmethode (BIEM) 2.3 Deterministische und stochastische Modellierung 2.3.1 Hydrodynamische Dispersion in porösen Medien 2.3.2 Dispersion und Varianz der Geschwindigkeitswerte 2.3.3 Stochastische Prozesse 2.3.4 Die 2 Grundprobleme der Transportmodellierung 2.3.5 Das Konzept der Makrodispersion

2

2.4 Verfahren zur optimalen Prozesssteuerung und Parameteridentifikation - optimale Prozesssteuerung, numerische Suchverfahren (Gradienten, Gauß-Newton...) 2.5 Programmsysteme zur Simulation von Geofiltration und Geomigration für den mehrdimensionalen Fall Literaturverzeichnis 3. Anhang - Lösungsbeispiele - 3.1. Lösung der Grundgleichung (homogene DGL 2.Ordnung, parabolischen Typs, linear mit konstanten Koeffizienten, eindimensional) - explizite Lösungsmethode 3.2. Temperaturausbreitung im Untergrund mit verschiedenen Randbedingungen 3.3. Lösung der Strömungsgleichung für die Grundwasserströmung mit praxisrelevanten Parametern. (explizites Lösungsverfahren) 3.4. Strömungsgleichung - implizites Lösungsverfahren - semi-implizite Lösung (CRANK-NICHOLSON) 3.5. Diffusionsgleichung mit verschiedenen Randbedingungen Laplace-Transformation und numerische Rücktransformation nach STEHFEST 3.6. Lösung der Transportgleichung (eindimensional) - explizites Lösungsverfahren - implizites und semi-implizites Lösungsverfahren - Gauß-Seidel-Iteration - Gaußsches Eliminationsverfahren 3.7. Lösung der Strömungsgleichung für ein nichtlineares Problem 3.8. Parameterbestimmung mit Suchverfahren am Beispiel des Diffusionskoeffizienten - Gauß-Newton-Verfahren - Powell-Verfahren

3

0. Einleitung und Aufgabenstellung Es ist das Anliegen der vorliegenden Praktikumsanleitung, eine Gruppe von physikalischen Bewegungs- und Transportvorgängen, die durch mathematisch ähnliche Modellgleichungen charakterisiert sind, einheitlich darzustellen und für vereinfachende Voraussetzungen numerische Lösungen aufzuzeigen. Dazu gehören: - die Wärmeleitung in Festkörpern, - die Strömung von Flüssigkeiten und Gasen, - die isotherme Diffusion in Festkörpern, Flüssigkeiten und Gasen, - der konvektive Massentransport (Stofftransport) in Flüssigkeiten und Gasen, - Parameteridentifikation (Lösung der inversen Aufgabe) In allen Beispielen wird sich auf laminare Strömungsvorgänge in porösen Medien ohne Trägheitseinfluss im ein- oder mehrdimensionalen Fall beschränkt. Der vorliegende Programm- und Übungsteil soll an die Bearbeitung praxisrelevanter numerischer Simulationsmodelle heranführen. Voraussetzung für eine erfolgreiche Beschäftigung mit den in dieser Anleitung zusammengestellten Übungen und Programmen ist, neben der Kenntnis der mathematischen Grundgleichungen, die Kenntnis der Programmiersprache C bzw. C++. Die Programme und Übungen sind so einfach wie möglich gehalten und sollen zu selbst-ständigen Erweiterungen anregen. Als Grundlage der Praktikumsanleitung dienten die Fachbücher von HÄFNER, SAMES & VOIGT "Wärme- und Stofftransport" /1/, KINZELBACH "Numerische Methoden zur Modellierung des Transportes von Schadstoffen im Grundwasser" /2/ und RICHTER "Modelle für Prozesse im Boden" /3/. In der o.g. Literatur wird jeder Nutzer weiterführende und fundierte Erläuterungen zu den speziellen Aufgabenstellungen finden. Frau Dipl. Math. S. Boy, Herr Dipl. Ing. A. Matwijenko und Diplomphysiker Th. Wagner erstellten und testeten die Übungsprogramme. Die Programme wurden im Rahmen von Lehrveranstaltungen durch Studenten (S. Waage, A. Klauke und E. Bennewitz) ergänzt und sind auf jedem Personalcomputer anwendbar. Frau S. Hengst fertigte die Abbildungen an.

1. Auflage, April 1997 2. korrigierte Auflage, 2000 3. Auflage, 2006

Symbolverzeichnis (Alle Maßeinheiten gelten in SI-Basiseinheiten: m, s, kg,) A Fläche m2

a Temperatur-, Druckleitfähigkeit m2/s allg. Leitzahl C Volumenkonzentration, Teildichte kg/m3

CM Massenkonzentration c spezifische Wärme J/kg⋅K D allgemeine Leitfähigkeit hydrodynamischer Dispersionskoeffizient m2/s Dm molekularer Diffusionskoeffizient m2/s D* mechanischer Dispersionskoeffizient m2/s erf(x) Gaußsche Fehlerfunktion erfc(x) komplementäre Gaußsche Fehlerfunktion erfc(x) = 1 - erf(x) g Erdbeschleunigung : g = 9,80665 m2/s h Wasserspiegelhöhe, Wasserstand k allgemeiner Abbaukoeffizient s-1

K, k (r) Durchlässigkeit, (relative) m2

kf Durchlässigkeitsbeiwert m/s kd Verteilungskoeffizient m3/kg L Länge m M Mächtigkeit m m Masse kg

.m Massenstrom kg/s

.

Am flächenbezogener Massenstrom kg/s⋅m2

Vm.

volumenbezogener Massenstrom (Massenstromdichte) kg/s⋅m3

n Porenanteil, Porosität →

n Normalenvektor p Druck Pa Q allgemeiner Stromterm q allgemeiner Quell/Senkenterm r Radius, allgemeine Koordinatenrichtung m R Retardationsfaktor beim Stofftransport S allgemeiner Speicherkoeffizient T Temperatur K t Zeit s u allgemeine Lösungsfunktion V Volumen m3

.V Volumenstrom m3/s

4

AV.

flächenbezogener Volumenstrom m/s (Geschwindigkeit)

VV.

Volumenstromdichte s-1

w Geschwindigkeit m/s zg Realgasfaktor α Wärmeübergangskoeffizient W/m2⋅K Kolmationsfaktor s-1

Г Rand des Gebietes δ Dispersivität m Ω zwei- oder dreidimensionaler Transport- raum /Gebiet oder Raum) m2 oder m3

ρ Dichte kg/m3

η dynamische Viskosität Pa⋅s λ Wärmeleitfähigkeit W/m⋅K κ isotherme Kompressibilität Pa-1

Indizes Fl Flüssigkeitsphase f Feststoffphase g Gasphase L bei x = L (Rand) R Rand oder bei r = R V auf das Volumen des Transportraumes bezogen 0 bei x = x0 oder r = r0t total, Gesamt- x,y,z Komponente eines Vektors in dieser Koordinatenrichtung - kart. Koordinaten r, ψ, θ - Kugelkoordinaten →

e Einheitsvektor

5

1. Lösung der Strömungs- und Transportgleichung Laminare Strömungs- und Transportvorgänge in homogenen und inhomogenen, porösen Medien. 1.1. Grundgleichungen 1.1.1 Fourier’sches Gesetz der Wärmeleitung

TTgradQ A ∇⋅−=−= λλ.

- Wärmestrom je Flächeneinheit, W/mAQ.

2

λ - Wärmeleitfähigkeit des Mediums, W/m ⋅ k 1.1.2 Die Wärmeleitungsgleichung

( ) VQtTcTdiv

.−

∂∂

=∇⋅ ρλ

- Wärmestrom je Volumeneinheit infolge Quellen / Senken VQ.

c - spezifische Wärme eines inkompressiblen Mediums, J / kg⋅K 1.1.3 Das Darcy-Gesetz der Filterströmung - Geoströmung in porösen Medien

( )zgradgpgradKw ρη

+−=

ρ - Fluiddichte, kg/m3

η - dynamische Viskosität des Fluids, Pa⋅s g - Erdbeschleunigung, m/s2

p - Druck, Pa z - vertikale Koordinate, m K - Durchlässigkeit, m2

- fiktive Geschwindigkeit (Darcy-Geschwindigkeit), m/s →

w

⎮w⎮ = /A ; - Volumenstrom, m.

V.

V 3/s und A - Querschnittsfläche, m2

6

1.1.4 Die Differentialgleichungen der Filterströmung Kontinuitätsgleichung

( )

VA mtnmdiv

..=

∂∂

+⎟⎠⎞

⎜⎝⎛ ρ

mit wm A ρ=.

n - Porosität Strömung einer homogenen, kompressiblen Flüssigkeit unter Vernachlässigung der örtlichen Dichteänderung (grad ρ = 0) : κ - Gesamtkompressibilität, κ = κFl + κf , Pa-1

Isotherme Strömung eines realen Gases

0

00

.

TTp

Vtp

zpnpgrad

zpKdiv

gg

−∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛ κη

(po, To) - Bezugszustand zg - Realgasfaktor 1.1.5 Das Fick’sche Gesetz der Diffusion

CgradDm mA −=.

Dm - molekularer Diffusionskoeffizient, m2/s

Für Volumenkonzentration einer Komponente i: imAi gradDm ρ−=.

1.1.6 Die Diffusionsgleichung

( ) Vii

imi mt

CCgradDdiv

.−

∂∂

=

Vim.

- Massenstrom je Volumeneinheit (Quelle/Senke)

7

1.1.7 Die Stofftransportgleichung Allgemeine Stofftransportgleichung

( )[ ] Viiimi mtCnCwCgradDDndiv

.* −

∂∂

=−+

D* - Dispersionskoeffizient, m2/s Wärmetransportgleichung

( )( ) ( )[ ] ( ) VtFlFl QtTcTwcTgradDcdiv

.* −

∂∂

⋅=⋅−⋅+ ρρρλ

(ρ⋅c)Fl - spezifische Wärmekapazität des Fluids (ρ⋅c)t - totale Wärmekapazität (Fluid/Feststoff) Tabelle1 : Allgemeine Strömungsgleichung – Symbolbedeutung

Dq

tu

augraddiv −

∂∂

=1

Strömungsvorgang u a q

Wärmeleitung Temperatur T Temperaturleit-

fähigkeit, D=λ

Wärmestrom je

Volumeneinheit

Filterströmung

kompressible

Flüssigkeit

Druck p Druckleitfähigkeit

a=k / nηκ

Volumenstrom je

Volumeneinheit

Grundwasser

gespannt

Spiegelhöhe h

(Druckhöhe)

Druckleitfähigkeit

a = kf / S0

Volumenstrom je

Volumeneinheit

Grundwasser

ungespannt

Spiegelhöhe h

(Wasserstand)

a = kf H / n Volumenstrom je

Flächeneinheit

Reales Gas Druck p Druckleitfähigkeit

a=k / nηκ

Quell-Senkenterm

Realgasströmung

Isotherme Diffusion Volumenkonz. C Diffusionskoeffizient

a=D, D=Dm, S=1

Massenstrom je

Volumeneinheit

8

Tabelle 2 : Differentialgleichung der Strömung in drei Koordinatensystemen

Lineare Strömung Dq

tu

axu

−∂∂

=∂∂ 1

2

2

Radialsymmetrische Strömung Dq

tu

aru

rru

−∂∂

=∂∂

+∂∂ 11

2

2

Kugelsymmetrische Strömung Dq

tu

aru

rru

−∂∂

=∂∂

+∂∂ 12

2

2

------------------------------------------------------------------------------------------------------------ Kartesische Koordinaten

qtuS

zuD

zyuD

yxuD

x zyx −∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

∂∂

+⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

∂∂

+⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

∂∂

-------------------------------------------------------------------------------------------------------------- Zylinderkoordinaten

qtuS

zuD

zuD

rruDr

rr zr −∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

∂∂

+⎟⎟⎠

⎞⎜⎜⎝

⎛Φ∂∂

Φ∂∂

+⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

∂∂

Φ2

11

------------------------------------------------------------------------------------------------------------ Kugelkoordinaten

qtuSuD

ruD

rruDr

rr r −∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛Θ∂∂

⋅ΘΘ

+⎟⎟⎠

⎞⎜⎜⎝

⎛Φ∂∂

Θ+⎟⎟

⎠

⎞⎜⎜⎝

⎛∂∂

∂∂

ΘΦ sinsin1

sin11

2222

2

-----------------------------------------------------------------------------------------------------------------

9

Tabelle 3 : Allgemeine Transportgleichung – Symbolbedeutung (HÄFNER u.a. 1992)

( ) qtuSukuwugradDdiv −∂∂

=−−r

Transport-prozeß

u

D

S

k

q

Wärme-transport Bemerkung : w = v∗(ρc)Fl

Temperatur T

Summari-scher Wärme-transport-koeffizient =λ+D*∗(ρc)Fl

Spezifische Wärme-kapazität = (ρc)t

Koeffizient abhängiger Quellen und Senken

Auf Volumen-einheit bezogener Strom unabh. Q/S

Stofftrans-port

Volumen-konzentra-tion C

Hydrodyna-mischer Dispersions-koeffizient =Dm + D*

Porenanteil S = n

Koeffizient abhängiger Quellen und Senken

Auf Volumen-einheit bezogener Strom unabh. Q/S

1. 2. Transportmechanismen und ihre mathematische Beschreibung 1.2.1. Allgemeine Grundlagen, Modelle und Modellierungsmethode Geoströmungstechnik Grundaufgabe : Vorausberechnung (Simulation) gegeben: Geometrie, Struktur, Parameter (z.B. Durchlässigkeit, Porosität) Differentialgleichung ortsabhängige Anfangsbedingungen orts- und zeitabhängige Randbedingungen Lösungsverfahren gesucht : Strömungsgrößen als Funktion des Ortes und der Zeit (Druck- bzw. Standrohrspiegelhöhenverteilung, Geschwindigkeitsverteilung, Dichteverteilung, Berechnung der bestimmten Randbedingungen zugeordneten Randwerte (z.B. Massenstrom oder Spiegelhöhe/ Druck))

10

11

Umkehraufgabe: Parameterbestimmung aus Betriebs- bzw. Feldmessungen oder Labormessungen (Identifikation) gegeben: Geometrie, Strömungsverhalten (orts- und zeitabhängige Anfangs- und Randbedingungen (Potentiale und Ströme) und Potentiale im Strömungsbereich) Lösungsverfahren gesucht: Parameter abh. vom Ort u./o. von der Zeit (Porosität, Durchlässigkeit, Grundwasserneubildung), unbekannte Geometrie, Stoffparameter, Randbedingungen 1.2.2. Transporterscheinungen, Transportprozesse Eine Transportgleichung besteht aus der Erhaltungsgleichung (Erhaltungssatz) mit Speicher-, Strom- und Quellsenkenterm, für die Ansätze mit bezogenen Zustandsvariablen (z.B. Dichte, Druck, Geschwindigkeit, Konzentration und Temperatur) erforderlich sind. Zu ergänzen ist die Zustandsgleichung, eine Beziehung zwischen allen auftretenden Variablen (außer Geschwindigkeit). Wir betrachten ein ruhendes unveränderliches (repräsentatives) Volumenelement (RVE). Die zeitliche Änderung von Masse in einem RVE ist im Zeitbereich von t0 bis t gleich der über die Fläche A des Elementes strömenden Massendifferenz und dem Element durch Quellen oder Senken zugeführten Masse. Speicherterm = Stromterm + Quellsenkenterm Speicherterm - Inhalt der Zustandsvariablen in V Stromterm - Konvektion + Diffusion Quellsenkenterm - Phasenwechsel, Sorption, Reaktionen, Wirkung äußerer Felder (Gravitation) Die Konvektions-Diffusions-Gleichung kann allgemein dargestellt werden in Differenzenform - als örtl.-zeitliche Integraldarstellung - als Differentialgleichung Modellbildung Modell Abbild des zu untersuchenden Originals (zur Modelltheorie siehe z. B.: L. Luckner u.a.1972, 1975...) /4/, /5/ materielle Modelle - Abbild von Strukturen (Sandmodelle)

12

ideelle Modelle - mathematische u. kybernetische Modelle nach der Art des Originals : z.B. geoströmungsmechanisches Modell physikalische Modelle (mit Analogie - funktions- oder verhaltensanalog) Elektro-Analogie : Ohm'sches Gesetz - Darcy'sches Gesetz mathematische Modelle - Widerspiegelung der dynamischen Originalprozesse (ideelle Modelle) Das mathematische Modell dient sowohl der Größenberechnung (z.B. p(x, t)) als auch der Parameterberechnung. Die mathema- tischen Modelle der Grund- und Umkehraufgabe unter- scheiden sich vor allem im Lösungsverfahren. Simulation - wirklichkeitsgetreue Nachahmung eines Vorganges (Strömung in porösen Stoffen) Lagerstättenmodell -------- Simulationsmodell Das Simulationsverfahren muss auf Simulationsmodell anwendbar sein. physikalische Gesetze : Erhaltungssatz (Massenerhaltung) Darcy-Gesetz Zustandsgleichung Differentialgleichung --- Differenzengleichung Lösungsmethoden / Simulationsverfahren - mathematisch - analytische Lösungen - analytisch - numerische Lösungen - diskrete numerische Simulationsverfahren Simulationsmodell - 1, 2, 3 - dimensional Begrenzung - innere u. äußere Ränder RB 1., 2. oder 3.Art Bohrung,Vorfluter - Quelle/Senke Zeitfunktion - instationär/quasistationär/stationär Gitternetzstruktur - orthogonale Vierecke, Dreiecke... Elementstruktur - Quaderelemente......

Simulation eines Strömungsvorganges – Modellanpassung durch Parameteridentifikation Messung, subjektive Parameterwahl, mathematische Verfahren Parameter : Durchlässigkeit/Transmissibilität Porosität/Speicherkoeffizient Speisung/GW-Neubildung/Querströmungsrate - Simulation von Prognosezuständen /(Geschichte) - Ergebnisbewertung und Ergebnisauswertung Stofftransport Strömungsprozess- Transportprozess Wirken von Wärmeleitung bzw. Diffusion + Konvektion physikalische Grundlage : Erhaltungssätze für die Energie bzw. Masse Zwei Gründe für Trennung in Strömung u. Transport : 1. Gleichungen mathematisch unterschiedlich Strömungsgleichung - parabolisch Transportgleichung – parabolisch-hyperbolisch → unterschiedliche Lösungsmethoden 2. deutliche Trennung in den Anwendungen Prozesse der Strömung ( Wärmeleitung, Diffusion und Filterströmung ) gekoppelte Wirkung (Transport) Konvektion und Dispersion konvektiver oder advektiver Transport - Transport von Masse oder Wärme durch eine

Fluidbewegung mit der mittleren Geschwindigkeit . →

wAls Dispersion (Zerstreuung) bezeichnen wir solche Prozesse, die durch die statistische Verteilung der Geschwindigkeit im Inneren des Strömungsraumes hervorgerufen werden. Konvektion 1. erzwungene Konvektion : Druck- oder Schweregradient (Potentialgradient) 2. freie Konvektion : Temperaturgradient erzwungene Konvektion:

iAi Cwm ⋅=.

Mechanische Dispersion : konfigurierter Strömungsraum (poröse Medien) mit statistischer Geschwindigkeitsverteilung.

13

Die mathematische Beschreibung der Stoffdispersion ist analog der Diffusion.

.

m Ai = -D* ⋅grad Ci D* - mechanischer Dispersionskoeffizient, m2/s -- kein Stoffwert !! praktikable Näherung für die Definition der Dispersivität

D* = δ ⋅| |→

w υ allgem. υ =1 Die Stofftransportgleichung Komponente i eines Stoffgemisches aus i =1, 2,.,N Komponenten

Strömungsgeschwindigkeit des Stoffgemisches →

w n < 1 für porösen Strömungsraum n = 1 für den freien Raum (z.B. Transport in einem Rohr)

.

m Ai Summe der Anteile von Konvektion, molekularer Diffusion und der mechanischen Dispersion Daraus ergibt sich die Kontinuitätsgleichung zu:

( )[ ] Viiimi mtCnCwCgradDDndiv

.* −

∂∂

=−+

Im allgemeinen Fall ist zu beachten, dass der Diffusionskoeffizient Dmi und der Dispersionskoeffizient D* Tensorgrößen sind. Für ein zweidimensionales, orthogonales Koordinatensystem, in dem die x-Achse mit der Stromlinie zusammenfällt (Hauptrichtungssystem), gilt:

( ) ( ) Vyyymxxxxm mtCn

yCDDn

yCw

xCDDn

x ii

.** −

∂∂

=⎥⎦

⎤⎢⎣

⎡∂∂

+∂∂

+⎥⎦

⎤⎢⎣

⎡−

∂∂

+∂∂

D*

XX = D*L = δL · wX longitudinaler Dispersionskoeffizient

D*

YY = D*T = δT · wY transversaler Dispersionskoeffizient

Quellen und Senken Zu - bzw. Abführungen von Masse oder Wärme im betrachteten System. Quellen / Senken = Massen- oder Volumenströme Wärmeströme je Volumeneinheit des Transportraumes

14

unabhängige Quellen und Senken z.B. Zustrom Schadstoffkomponente / Konzentration Zuführung von Wärme / Temperatur, Injektion von Fluid

abhängige Quellen / Senken iZiVi Cm ⋅−= λ.

z.B. radioaktiver Stoffzerfall, λZi - Zerfallskonstante der i-ten Komponente lineare und nichtlineare Stoffsenken z.B. Denitrifikation mathematische Darstellung : Stofftransport / Massenstrom

t

VVQiVV

mVmitCVm i ρ

....

=⋅=

Wärmetransport : ( ) QFlVV TcVQ ⋅⋅⋅= ρ..

- Volumenstrom der Quelle je Volumeneinheit VV.

Vim.

- Massenstrom, Ρt - totale Dichte (Gemischdichte) TQ , CQi - Temperatur, Konzentration der Quelle Linearität: Die Gleichung ist linear, wenn die Lösungsfunktion u(x, y, z, t) nur in der 1. Potenz auftritt, d.h. wenn die Koeffizienten D, S, a und q nicht von u abhängen bzw. q nur von u in der 1. Potenz abhängig ist. Homogenität : q = 0 stationär - instationär - quasistationär

Zeitableitung tu∂∂ = 0 stationär

≠ 0 instationär = M quasistationär ( M < ∞ , M = const.) Anfangs - und Randbedingungen Kenntnis von Bedingungen zu Beginn des Prozesses und an den Rändern des betrachteten Raumes. Anfangsbedingungen (AB) homogene AB: u(x, y, z, t) = ua = const. für t=0 durch Transformation der Lösungsfunktion lässt sich die homogene AB stets in der Form schreiben: u(x, y, z, t) = 0 für t=0 inhomogene AB : u(x, y, z, t) = ua (x, y, z) für t=0 (ortsabhängige Verteilung )

15

Randbedingungen (RB) - Randwerte RB stellen die Kenntnis der Lösungsfunktion u bzw. ihrer Ableitungen an den Rändern des betrachteten Raumes dar. Als Randwerte (RW) bezeichnet man die Lösung der jeweiligen Aufgabe auf dem Rand. Wir unterscheiden Randbedingungen 1., 2. und 3.Art. RB 1. Art: Lösungsfunktion u auf dem Rand für alle Zeiten t bekannt: u (x, y, z, t) = uR(x, y, z, t) bei (x,y,z) = Rand ist uR zeitunabhängig, so nennt man die RB stationär RB 2. Art: Im mathematischen Sinne bekannte Normalenableitungen der Lösungs- funktion

( ) [ ] RandzyxbeiFunktionevorgegeben

ntzyxu

==∂

∂ ,,,,,

n - Normalenrichtung Im physikalischen Sinne ist der Strom (Wärmestrom, Massenstrom, Volumenstrom) über den Rand bekannt. Die RB 2. Art für die allgemeine Strömungsgleichung lautet:

( ) ( ) [ ] RandzyxbeitzyxQn

tzyxuDRA

R

==∂

∂− ,,,,,,,,

wobei QA - allgemeiner flächenbezogener Strom Die RB 2. Art für die allgemeine Transportgleichung lautet:

( ) ( ) ( )[ ] [ ] RandzyxbeitzyxuwtzyxQn

tzyxuD RRAR

=⋅−=∂

∂− ,,,,,,,,,,,

RB 3. Art: (gemischte RB) weder Lösungsfunktion noch Normalenableitung sind gegeben, jedoch eine Linearkombination beider. z.B. : Grenzschicht mit Wärmeübergangswiderstand, Wärmeübergangskoeffizient / Strömungswiderstand dünner Schichten – Kolmationsschicht. Die RB 3. Art für die allgem. Strömungsgleichung lautet:

( ) ( ) ( ) [ ] Randzyxbeitzyxutzyxun

tzyxuDRR

R

==+∂

∂− ,,,,,,,,,,,α

16

α- Übergangskoeffizient Die RB 3. Art für die allgem. Transportgleichung lautet:

( ) ( )[ ] ( ) [ ] RandzyxbeitzyxQtzyxuwn

tzyxuDRAR

R

==⋅−∂

∂− ,,,,,,,,,,,

QA als Gesamtstrom über den Rand bekannt (Strom je Fläche) Diskussion der Randbedingungswahl für Transportprozesse: RB 2.Art sinnvoll, wenn Konvektionsterm am Rand vernachlässigbar. RB 3.Art sinnvoll, wenn Konduktion und Konvektion wirksam werden. RB 1.Art sinnvoll, wenn konvektiver Transport überwiegt. Am Ausströmungsrand ist die Vorgabe eines physikalisch sinnvollen Stromes Q im allge-meinen nicht möglich (Ausnahme im Unendlichen). Hier wird die so genannte Transmissions - Randbedingung vorgeschlagen (aus dem Inneren kommender Konduktions- Diffusionsstrom fließt auch über den Rand aus). Randbedingungen im Unendlichen : Keinerlei Einflüsse aus dem Inneren eines Gebietes erreichen den äußeren Rand. Für Strömungs- und Transportprozesse bedeutet es, dass der Strom über den äußeren Rand unveränderlich und gleich dem Wert im Anfangszustand ist: QA(x, y, z, t) = QA(x, y, z, t=0) bei (x, y, z) = Rand im Unendlichen, wobei der Strom endlich groß bleiben muss. Zur Erfüllung dieser Bedingung genügt es, die Endlichkeit der Variablen u (RB1.Art), ihrer Normalenableitung ∂u⁄∂n (RB 2.Art) oder einer Linearkombination (RB 3.Art) zu verlangen. Sowohl für Strömungs- als auch für Transportprobleme kann man für einen Rand im Unendlichen die folgenden RB als äquivalent betrachten: RB 1. Art ( ) ∝<= Rutzyxu ,,,

RB 2. Art ( )∝<=

∂∂

AQn

tzyxuD ,,,

RB 3. Art ( ) ( ) ∝<=+∂

∂⋅− Rutzyxu

ntzyxuD ,,,,,,

α

bzw. ( ) ( ) ( ) ∝<=⋅+∂

∂⋅− tzyxQtzyxuw

ntzyxuD An ,,,,,,,,,

17

wobei uR und |QA| beliebige Konstanten sind, in der Regel = 0 Randwert- und Anfangs-Randwertprobleme Unter einem Randwertproblem versteht man die entsprechende DGL. mit den zugehörigen RB. Randwertaufgaben sind meist stationäre Vorgänge, bei denen eine Anfangsbedingung physikalisch ohne Einfluss ist, z.B. ein Strömungsvorgang nach unendlich langer Zeit. Beispiel: Randwertproblem der eindimensionalen, stationären Strömung mit RB 1.Art

RB: u = u0 bei x = 0 Dq

dxud

−=2

2

u = uL bei x = L Anfangs- und Randwertprobleme ergeben sich für instationäre Strömungs- und Trans-portprozesse. Sie bestehen aus der DGL und den die Aufgabe charakterisierenden AB und RB. Beispiel : eindimensionaler Transport mit RB 1.Art

AB : u = 0 für t = 0 und x > 0 qtuSuk

xuw

xu

x −∂∂⋅=⋅−

∂∂⋅−

∂∂

2

2

RB : u = u0 bei x = 0 und t > 0 u = uL < ∞ bei x = L, t > 0 Insgesamt besteht die Anfangsrandwertaufgabe aus 6 Gleichungen mit 6 Variablen : Massentransport = Kontinuitätsgleichung ρ Impulstransport = Darcysches Gesetz w Zustandsgleichung Fluid p Zustandsgleichung Mischphase fluid-fest n Phasendruckgleichung peff Gesamtdruckgleichung pges den Materialgleichungen für : dynamische Viskosität η Durchlässigkeit k

18und den Anfangs- und Randbedingungen

19

1.2.3 Mehrphasen - Strömung – und Stofftransport

Sedimentgestein / Festgestein : porös, porös-klüftig, klüftig

ungesättigte Zone / gesättigte Zone

Grundfall : Einkomponenten- Einphasen- Systeme, Tracerfall

Erweiterter Grundfall : Kopplung ungesättigter und gesättigter Bereich

Kopplung von Grund- und Oberflächenwasser

Erfassung von Temperatureinflüssen

Dichteabhängige Strömung

Geochemische Prozesse, Zerfallsketten

Mechanische Deformation

Veränderliche Fluideigenschaften

Modellierungsaufgaben

♦ Idealisierung zum hydrogeologischen Modell

♦ Strukturierung zum mathematischen Modell „Konzeptionelles Modell“

♦ Ermittlung der hydraulischen Parameter

♦ Ermittlung der Migrationsparameter

♦ Durchführung der Berechnung „analytisches / numerisches Modell“

♦ Szenariountersuchungen „generisches Simulationsmodell“

♦ Interpretation der Berechnungsergebnisse

♦ Prognoserechnungen

Strömungs- und Transportprozesse in der ungesättigten Zone Basisproblem / Grundfall der Mehrphasenhydraulik :

♦ Versickerndes Niederschlagswasser (Transport der gelösten Schadstoffkomponenten in

der wässrigen Phase Bodenwasser im mit Luft teilgesättigten Untergrund /

Aerationszone)

♦ Transport in einer nichtwässrigen flüssigen Phase (NAPL-nonaqueous phase liquid) –

kohärent verteilte Phasen.

Mehrphasen - Mehrkomponenten-Systeme nichtmischbarer Fluide

( RICHARDSON-Gleichung für Infiltrationsprozesse mit konstantem Druck (Luftdruck) der

Gasphase-Luft : Grundfall der Gas–Wasser-Strömung für ungesättigte Zone)

Richard‘s-Gleichung : p (Gas/Luft) = 0, d.h. pnb (r,t) = pnb,0 (t)

pc = pnb - pb Kapillardruck

( )t

ppCt

ppt

tp cc

c

c

c

∂∂⋅=

∂∂⋅

∂∂

=∂

∂ θθ ),(

C (pc) – Parameterfunktion der (kapillaren) Speicherkapazität der benetzenden Phase (erste

Ableitung der Zustandsfunktion θb = f (pc)).

Mit pnb (t) = const. gilt dpc = -dpb und es ergibt sich:

qt

ppCzgradpgradg

kdiv bbbr −

∂∂⋅=⎟⎟

⎠

⎞⎜⎜⎝

⎛⎟⎟⎠

⎞⎜⎜⎝

⎛+ )(1

ρ

Die Gasphase ist somit nur noch implizit durch die Beziehung θw + θg = 1 gegeben.

Parameterfunktionen Sättigung – Kapillardruck (LUCKNER et al., 1989, VAN GENUCHTEN, 1980,

MUALEM,1976, BROOKS & COREY, 1964)

Hysterese

♦ Unregelmäßige Porengeometrie führt zu Unterschieden in Drainage / Imbibition

♦ Kontaktwinkel zwischen Drainage und Imbibition sind unterschiedlich

20

Durchlässigkeit – Sättigung Mobilität einer Phase kM = spezif. Permeabilität / dynamische Viskosität

sPamin

gkkk f

M ⋅⋅==

2

ρη

Relative Durchlässigkeit kr

= Verhältnis der Durchlässigkeitsbeiwerte bei Teil- und Vollsättigung.

Benetzende Phase : kr (θb) VAN GENUCHTEN (1980)

Die Kurven für die relativen Permeabilitäten und Sättigungen in Abhängigkeit von den

Saugspannungen (Druckhöhen) nach VAN GENUCHTEN, MUALEM und BROOKS &

COREY werden in den Simulationsprogrammen in verschiedener Form genutzt.

Strömung von nicht mischbaren Flüssigkeiten

Mischbar : Salz- und Süßwasser, warmes / kaltes Wasser

Nicht mischbar : Wasser / Öl, Wasser / Chlorkohlenwasserstoffe

NAPL – non aqueous phase liquids

LNAPL – light-NAPL (geringere Dichte als Wasser)

DNAPL – dense-NAPL (höhere Dichte als Wasser)

Zweiphasenströmung Gas-Wasser

Bilanzgleichungen für die Wasser- und Gasphase :

( )q

tnpgrad

Kdiv w

ww

ww −∂∂

=⎭⎬⎫

⎩⎨⎧ θ

ηθ

Wasserphase

21

Gasphase :

( )q

tngrad

ddp

pgradK

div ww

w

capw

g

wg −∂∂

−=⎪⎭

⎪⎬⎫

⎪⎩

⎪⎨⎧

⎟⎟⎠

⎞⎜⎜⎝

⎛+

θθθη

θ

pg = pcap + pw und θg + θw = 1

n – Porosität, Kg,w – relative Permeabilitäten in m2 kf,w = Kw ρw g / ηw

(Korrelationen nach Brooks & Corey, van Genuchten, Mualem, Luckner et al.)

3D- Modellierung der Mehrphasen- und Mehrkomponentenströmung mit

thermodynamischem Gleichgewicht mittels Programm MULTIF – TWOPHASE (Nekrassov,

TU Bergakademie, 2001).

3D- Zweiphasenmodell : Gas/Luft-Wasser mittels TWOPHASE in HOTH (2003).

Die relativen Durchlässigkeiten zeigen eine typische Abhängigkeit vom Verhältnis der

Sättigung der (beiden) Phasen. Eine gleichzeitige Bewegung (beider) Phasen ist demnach nur

möglich, wenn die Wassersättigung größer als die irreduzible Sättigung (Haftwasser ca. 30%)

und die Ölsättigung größer als die so genannte Rest- oder Residualsättigung ist (ca. 10%

Ölsättigung).

Fällt die Sättigung der nicht benetzenden Phase auf den Wert der Restsättigung, geht die

relative Durchlässigkeit gegen Null. Die Phase liegt praktisch fest.

Der Phasenkörper kann jedoch vom Wasser durchströmt werden, wobei eine sehr große

Kontaktfläche entsteht. Dabei bildet die nicht mischbare Phase in Abhängigkeit ihrer

Löslichkeit in Wasser eine permanente Quelle für gelöste Inhaltsstoffe, die vom

Grundwasser aufgenommen und weiter transportiert werden.

Zwei nichtmischbare mobile Phasen sind in einem porösen Medium so verteilt, dass die

besser benetzende Phase die kleineren und die schlechter benetzende Phase die größeren

Poren einnehmen.

22

Parameterbeschreibung und – Ermittlung in der ungesättigten Zone

Bestimmung des relativen K-Wertes als Funktion der Sättigung und die Bestimmung der

Saugspannung als Funktion der Sättigung auf der Basis experimentell gewonnener Mess-

werte: Restwassergehalt / Sättigung / Anpassungsparameter

Kapillarkraft : Tensiometermessungen zur Bestimmung des Gradienten der

Kapillarspannung

Multistep-Flow-Test (Saugverfahren), Infiltrationsexperimente

Optimierungsalgorithmen z.B. LEVENBERG-MARQUARDT-Verfahren

Dichteabhängige Strömung Massenerhaltungsgesetz für das Fluid :

div wt

n m( ) ( ) &ρ ∂∂

ρ= − . rw

www

x

y

z

=⎛

⎝

⎜⎜⎜

⎞

⎠

⎟⎟⎟

Nach Definition der Gesamtkompressibilität χ χ χ= +f F , und unter Voraussetzung der

Boussinesq-Approximation (nichtdeformierbares poröses Medium), ergibt sich:

mtC

ntp

nzgradgpgradkdiv c &−+=⎥⎦

⎤⎢⎣

⎡+

∂∂

βρ∂∂

χρρρη 0)(

Dabei ist die Dichte ρ ( , )p C näherungsweise linear mit:

[ ]ρ ρ χ β( , ) ( )p C p p Co F o c= + − +1

( )

max0

max )0,(,C

pCp ooc ρ

ρρβ

−=

23

Transportgleichung aus der Massenbilanz der jeweiligen Komponente C :

Parametergleichung: Konvektion / Diffusion-Dispersion / Sorption / Stoffabbau

div D grad C w C n R C n RCt

mC[ ] &− − ⋅ ⋅ ⋅ = ⋅ −r λ∂∂

Die Transportgleichung berücksichtigt im Retardationsfaktor R die mögliche Adsorption nach zwei Isothermenkonzepten. Nach der HENRY-Isotherme ist:

R nn

kM d= +−1 1 ρ , nach der LANGMUIR-Isotherme:

( )R C n

nF FF C

Msat d

d

( ) = +− ⋅

+1 1

2ρ .

ρ Dichte, ML-3.

p Druck, MT-2L-1.

k Permeabilität, L2.

χ Kompressibilität des Fluidporenraumes, M-1T2L.

cm& Massenstromdichte der Quellen und Senken, ML-3T-1.

g Schwerebeschleunigung, LT-2.

z Koordinate, L.

Fsat Sättigungskonzentration am Feststoff, M L-3

Fd dyn. Langmuierkoeffizient, M L-3

kd HENRY-Verteilungskoeffizient, L3 / M

R Redardationsfaktor

η dyn. Viskosität, M L-1 T-1

24

1.3. Analytische Lösungen Bei HÄFNER u. a. /1/, van GENUCHTEN & ALVES (1982) /10/ werden eine Vielzahl analytischer Lösungen der Anfangs- und Randwertaufgabe für ein- und zweidimensionale Strömung, Wärmeleitung und Diffusion angeboten. Analytische Lösungen für den Wärme- und Stofftransport bleiben nur für den eindimensionalen Fall überschaubar. In /1/ werden deshalb für mehrdimensionale Transportprobleme nur Spezialfälle betrachtet (Geschwindigkeitsvektor parallel zu einer Koordinatenachse). Für Probleme der Wärmeleitung und des Temperaturausgleiches ist bei TAUTZ /11/ ein umfassendes Angebot analytischer Lösungen gegeben. Auch für komplizierte numerische Simulationsmodelle kann mittels analytischer Lösungen ein physikalisch sinnvoller Vergleich vorgenommen werden. In einer Vielzahl von Fällen wird die analytische Lösung der Problemstellung gerecht und kann der numerischen Lösung vorgezogen werden. In den Übungsaufgaben wird deshalb entweder auf die analytische Lösung verwiesen oder diese wird selbst programmiert. Weil in den letzten Jahren im Ingenieurwesen verstärkt die Methode der Laplace-Transformation angewendet wird, da sie auch für eine Reihe komplizierter Strömungs- und Transportprobleme relativ leicht handhabbar ist, soll der prinzipielle Lösungsweg im folgenden Abschnitt vereinfacht dargestellt werden. Laplace-TransformationDie Methode wird in /1/ S.72-92 ausführlich behandelt. In einem ersten Schritt wird das gesamte Anfangs-Randwertproblem (partielle Differentialgleichung mit Anfangs- und Randbedingungen) in den Laplace-Bereich transformiert.

( )[ ] ( ) ( )∫∝

− =⋅=0

sfdtetftf stL s : Laplace- Variable

f(t) : Zeit- oder Originalfunktion

( )[ ] ( )szyxftzyxf ,,,,,, =L ( )sf : Laplace-Transformierte oder Bildfunktion Haben die mit dieser Methode behandelten Probleme neben der Zeit nur noch eine unab-hängige Variable, so entsteht im Bildbereich aus der partiellen Differentialgleichung eine gewöhnliche Differentialgleichung mit einfacher Lösung. Für die Rücktransformation in den Zeitbereich (inverse Laplace-Transformation) gilt eine komplexe Umkehrformel.

( ) ( )[ ] ∫∝+

∝−

− ⋅==ic

ic

st dsefi

sftfπ211L

Ihre Anwendung führt auf Integrale, die analytisch oder auf numerischem Wege gelöst werden können. Die analytisch auswertbaren Integrale ergeben Korrespondenzen, von denen die wichtigsten in /1/ im Anhang C enthalten sind. In den meisten Fällen ist es jedoch rationeller, die Rücktransformation mittels numerischer Verfahren auszuführen.

25

1.3.1 Mathematische Lösung von speziellen Transportmodellen Randbedingungen: Einseitig begrenztes Gebiet (1. Art bei x = 0 ) RB im Unendlichen bei x = ∞ Lösung der Transportgleichung

26

t

cSckxcw

xcD

∂∂

=−∂∂

−∂∂

2

2

(partielle DGL 2. Ordnung parab. - hyperbol. Typs, inhomogen , instationär, linear) 0,00: >== xtfürcAB

0,0,0:.1 0

>∞=∞<Χ=

>==

txbeictxbeiccArtRB

Laplace –Transformation Lösung :

1. Multiplikation mit und Integration (damit alle Integranden für gegen

Null gehen)

ste− ∫∞

0

dt ∞→t

∫ ∫ ∫∫∞ ∞ ∞

−−−−∞

∂∂

=−∂∂

−∂∂

0 0 002

2

tdetcStdecktde

xcwtde

xcD stststst

Konstante Koeffizienten

( ) ( ) ( ) dtetcStdetxcktdetxc

xwtdetxc

xD stststst −

∞∞ ∞ ∞−−− ∫∫ ∫ ∫ ∂

∂=−

∂∂

−∂∂

00 0 02

2

,,,

∫∞

−

∂∂

=−−

27

02

2

tdetcSck

xdcdw

xdcdD st

Partielle Integration ∫∫ −= xdvuuvdxuv ´´

Partielle Integration des Speicherterms führt auf:

( )[ ] ∫∫∞

−∞−−∞

+=∂∂

00

0

, tdecSsetxcStdetcS ststst

Mit den Integrationsgrenzen und der Anfangsbedingung c = 0 bei t = 0 ergibt sich:

( ) ( )[ ] cSsetxcetxcS +=−∞== −∞ 00,, = 0 = CAnfang = 0

Transformierte DGL

( ) 0

.

2

2

2

2

=+−−

=−−

cSskxdcdw

xdcdD

bzw

csSckxdcdw

xdcdD

Zur Vereinfachung führen wir eine Zeit t´= St ein.

(DGL würde lauten: tcck

xdcdw

xcD

′∂∂

=−−∂∂

2

2

) und erhalten die transformierte Dgl. ohne S:

( ) 02

2

=+−− cskxdcdw

xdcdD

Transformation der Randbedingungen:

∫∞

∞−− =−===0

00

000,0

sc

es

ctdecccx stst

∞<⇒∞==∞=sXXCx

Lösung der transformierten DGL mit eμx – Ansatz :

( ) 02 =⋅+−⋅−⋅ xxx eskeweD μμμ μμ

Division durch führt auf die charakteristische Gleichung: ⇒xeμ

( ) 02 =+−− skwD μμ

Mit den Lösungen:

28

( )[ ]skDwwD

++±= 421 2

2/1μ 00

2

1

<>

μμ

Damit ergibt sich die allgemeine Lösung:

( ) xx eBeBsxc 2121, μμ +=

Einsetzten der RB zur Bestimmung der Koeffizienten:

∞∞ +=∞→ .2

.1

21 μμ eBeBsxx mit der Bedingung: Lösung < ∞ !

( )01 >∞→μ ( )0

0

2 <→μ

ergibt: BB1 = 0 und für x = 0 B2B = c0 / s

Mit der Gesamtlösung : ( ) xes

csxc 20, μ=

( ) xskDwxD

w

eesc ⎟

⎠⎞⎜

⎝⎛ ++−

=420

2

.

Rücktransformation aus dem L-Bereich in den Original-(Zeit)-Bereich mittels Korrespondenzen (z.B. Häfner, Sames, Voigt, 1992, S. 592 ff., Beziehung 9.3)

( ) ( )[ ] 000 ,,,,211 PvtxPvtxPe

sx =+⇒

+−

μ

Dkwv 42 +=

( ) ( )⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

∗

∗++=+ −

−− Dt

vtxerfcvw

DxvtxP

42exp,,0

(Beachte : und erfc – komplementäre Fehlerfunktion) Stt /=∗ Lösung des Transportproblems im Zeitbereich :

( ) ( )vtxPs

CtxC ´,,, 00= mit Dkwv 42 +=

( ) ( ) ( )⎪⎭

⎪⎬

⎫

⎪⎩

⎪⎨

⎧ +⎥⎦⎤

⎢⎣⎡ ++

−⎥⎦⎤

⎢⎣⎡ −=

stDs

tvxerfcvw

Dx

stDs

tvxerfcvw

DxCtxC

42exp

42exp

2, 0

1.3.2 Numerische Verfahren zur Rücktransformation Bei einem im Wesentlichen monoton verlaufenden Funktionsverlauf hat sich das Verfahren von Stehfest als robust und schnell erwiesen. Der Stehfest-Algorithmus ist für die Lösung von Strömungsproblemen sehr, für Transportprobleme nur bedingt geeignet (numerische Dispersion). Alternative Methoden zur numerischen Rücktransformation sind die Verfahren nach Schapery, Zakian und Talbot /1/. Im Lösungsbeispiel 2.2.5 wird die Diffusionsgleichung mittels Laplace-Transformation gelöst und numerisch mit STEHFEST in den Zeitbereich rücktransformiert.

29

2. Numerische Lösung von partiellen Differentialgleichungen Formulierung von diskreten Anfangs-Randwertproblemen und deren numerische Lösung. Diskrete Probleme bestehen aus einer Differenzengleichung mit ihren Anfangs- und/oder Randbedingungen. Die Differenzengleichung ist das diskrete Analogon zur Differential-gleichung und kann durch folgende Betrachtungsweisen gewonnen werden: - durch Taylor-Reihenentwicklung der Differentialquotienten in der jeweiligen Differentialgleichung - durch Diskretisierung der Integraldarstellung der Differentialgleichung - durch Ableitung am endlichen Element (Bilanzmethode) - über das Variationsprinzip Die Lösung des Systems von Differenzengleichungen ist das Kernstück aller numerischen Simulationsverfahren. Für jedes Volumenelement (Gitterpunkt) entsteht eine Differenzengleichung mit 3, 5, oder 7 unbekannten Druckwerten (1-, 2- oder 3-dimensional). Diese Gleichungen sind miteinander verkoppelt; das Gesamtsystem enthält so viele unbekannte Druckwerte, wie es Gitterpunkte im Inneren des Strömungsraumes besitzt. Die einfachste Lösung erfolgt mit Hilfe des Gauß'schen Eliminationsverfahrens, das aufgrund des hohen Berechnungsaufwandes durch eine Vielzahl von Lösungsalgorithmen ersetzt wurde. 2.1 Numerische Modellierung von Strömungs- und Transportprozessen Diskretisierung durch endliche Differenzen oder finite Elemente: Das Geschwindigkeitsfeld als stetige Funktion an jedem Punkt (x, y) erhält man durch Interpolationen zwischen den Zellen oder mittels Finite - Elemente - Methode (- flexible Diskretisierung durch beliebig geformte Elemente). Spiegelhöhe/Druck im Zentrum der Zelle - aus ihnen kann rückwärts mit dem Darcy-Gesetz und der effektiven Porosität die Abstandsgeschwindigkeit auf der Grenze zwischen zwei Zellen bestimmt werden.

( ) ( )ijjiijy

ijyijjiijx

ijx hhyn

kfwundhh

xnkf

w −Δ⋅

−=−Δ⋅

−= ++ 1,,

,,1,

,

Approximation durch : Interpolationsansatz/räumliche Interpolationsfunktion ωi (x, y) Daraus ergibt sich ein Geschwindigkeitsfeld als Funktion von y und x.

( ) ( ) ( ) ( )⎟⎠

⎞⎜⎝

⎛⋅

∂∂

= ∑=

n

iiix yxth

xnyxkftyxw

1,,,, ω und ( ) ( ) ( ) ( )⎟

⎠

⎞⎜⎝

⎛⋅

∂∂

= ∑=

n

iiiy yxth

ynyxkftyxw

1,,,, ω

Bei linearer Interpolation auf Dreieckselementen sind wx und wy elementweise konstant. Die Geschwindigkeitsvektoren sind damit an den Elementgrenzen unstetig. Für die zuverlässige Bestimmung des Geschwindigkeitsfeldes mit dieser Art von Elementen muss eine hinreichend

30

31

feine räumliche Diskretisierung vorausgesetzt werden. Alternativ sind Glättungsverfahren denkbar. Probleme ergeben sich bei der Interpolation sowohl in Differenzenverfahren als auch in Finite-Elemente-Verfahren an Singularitäten des Strömungsfeldes (z.B. Brunnen). Hier wird eine Einbettung analytischer Lösungen in die numerischen Lösungen empfohlen. Die Erstreckung des Strömungsfeldes ist in der Regel größer als die des Transportmodells. Das Transportmodell erfordert i.a. eine feinere Diskretisierung (Zoom-Technik, Interpolationsknoten (Pinch-nodes) oder echte lokale Netzverfeinerung). Ein Strömungsmodell, das Grundlage für ein Transportmodell ist, sollte höheren Genauigkeitsanforderungen genügen als ein Strömungsmodell für regionale Wasserbilanzen. Dispersionsfreie Näherung Der konvektive Transport von Schadstoffen im Grundwasser stellt die mittlere Schadstoffbewegung dar. Diese Betrachtungsweise entspricht einer Vernachlässigung des Dispersionsterms. Sie ist anwendbar, wenn in erster Linie die Richtung des Schadstoff-transports, die Bahnen von Fronten oder Schwerpunkten von Konzentrationsverteilungen und mittlere Ankunftszeiten an bzw. Laufzeiten zu einem Entnahmebrunnen interessieren (z.B. Schutzzonen, hydraulische Abwehr- und Sanierungsmaßnahmen). Bahnlinien und Laufzeiten Die Transportgleichung beschreibt die Entwicklung der Schadstoffkonzentration in einem Kontrollvolumen, das sich längs der Bahnlinie (x(t), y(t)) bewegt. Aus der Integration über die Bahnlinie erhält man die Laufzeit des Schadstoffes zwischen zwei beliebigen Punkten. Durchbruchskurven Bei allen wasserwirtschaftlichen Betrachtungen von Grundwasserverschmutzungen ist eine wesentliche Frage, in welcher Stärke und zeitlichen Abfolge eine Verschmutzung in Wasserfassungen, Quellen oder Vorflutern wieder an die Oberfläche zurück gelangt. Lösung der Transportgleichung Die hyperbolischen Eigenschaften der Transportgleichung bereiten Probleme bei der numerischen Lösung. Sie sind Anlass, über die Standardverfahren der Differenzenmethode und Finite-Elemente-Methode hinaus, weitere Lösungsverfahren anzuwenden, die sich an den Lösungsverfahren der rein hyperbolischen Wellengleichung orientieren. Verfahrensklassen: - Differenzenverfahren (FDM) - Bilanzmethode (Control-Volume-Method) - Finite-Elemente-Verfahren (FEM) - Randintegralgleichungsmethode (BIEM - Boundary Integral Equation Method) - Charakteristikenverfahren (MOC - Method of Characteristics) - Random-Walk-Verfahren (RWM) (-Laplacetransformationsverfahren (Luckner, Schestakow 1985 /6/, Sudicky,1989 /12/, Häfner, Sames, Voigt /1/))

32

Differenzenverfahren - einfachste Lösungsmethode , n Knoten, n Gleichungen Die Bilanzgleichungen können direkt durch Bilanzierung an der Gitterzelle oder durch Ersetzen der Differentialquotienten durch Differenzenquotienten gewonnen werden. Bei überwiegend konvektivem Massenfluss ist es sinnvoll, das Gewicht bei der oberstrom gelegenen Konzentration zu wählen. Man spricht von upwind - Differenzen (Rückwärts-differenzen). Bei überwiegend diffusiv/dispersivem Massenfluss sind zentrale Differenzen im Konvektionsterm geeigneter. Die Maßzahl für das Verhältnis zwischen konvektivem und dispersivem Massenfluss im Element (i,j) ist die Gitter- Pecletzahl Pe* : Pe*

x = wx ∆x / Dxx Pe*y = wy ∆y / Dyy

Dxx und Dyy – hydrodynamische Dispersionskoeffizienten Pe* < 1 überwiegend parabolischer Charakter des Transportproblems Pe* > 1 überwiegend hyperbolischer Charakter Statt einer Festlegung auf Rückwärts- bzw. Zentraldifferenzen ist auch eine Steuerung der Art der Differenzen durch die Pecletzahl oder Flux-Limiter Techniken möglich. Explizites Verfahren Jede Knotengleichung enthält nur einen Term mit der Unbekannten (Vorwärtsdifferenzen in der Zeit). Das explizite Verfahren hat den Nachteil, dass zu seiner Stabilität einschränkende Bedingungen an die Länge des Zeitschrittes erfüllt sein müssen (Bear, 1979 /13/): - Courant - Kriterium (garantiert, dass in einem Zeitschritt die Konzentration in einem Element nicht größer werden kann als die Konzentration in den konvektiven Zuflüssen bzw., dass nicht mehr Schadstoff die Zelle verlassen kann als in ihr zu Beginn des Zeitschrittes enthalten ist.) - Neumann - Kriterium (es lässt sich physikalisch als die Forderung interpretieren, dass während eines Zeitschrittes Konzentrationsgradienten durch dispersive Massenflüsse allein nicht umgekehrt werden können.) - allg. Stabilitätskriterium (am Entnahmeknoten soll nicht mehr Schadstoff das Element verlassen als ursprünglich vorhanden war. Entsprechend darf an Zugabeknoten nicht mehr Schadstoff eingetragen werden als das Element, ohne eine Konzentrationserhöhung über die Zugabekonzentration hinaus, aufnehmen kann.) - Abbaukriterium (Adsorptions- und Reaktionsterme) (im Fall der Abbaureaktion muss der Zeitschritt so klein gewählt werden, dass nicht mehr Schadstoff abgebaut werden kann als in der Zelle vorhanden ist.)

Implizites Verfahren - flexiblere Wahl des Zeitschrittes, Rückwärtsdifferenz in der Zeit oder zentrale Differenzen in der Zeit (Crank, Nicholson, 1947 /63/) - n simultane lineare Gleichungen für n unbekannte Konzentrationen. Zur Lösung des Gleichungssystems stehen direkte und iterative Methoden zur Verfügung. Direkte Gleichungslöser brauchen mehr Speicherplatz als iterative Methoden. Dafür ist bei letzteren, bei schlechter Konditionierung, der Zeitaufwand größer. Das allgemeine Gleichungssystem für den Transport ist im Gegensatz zum Gleichungssystem für die Strömung unsymmetrisch. Damit können spezielle Löser für symmetrische Matrizen zunächst keine Anwendung finden (erst nach geeigneter Transformation - Lösung durch Matrizenzerlegung (Lower und Upper Dreiecksmatrix). Iterative Gleichungslöser : ADI- bzw. IADI-Verfahren (Peaceman & Rachford, 1955 /14/) Blockiterationsverfahren (regelmäßige Rechteckgitter); Punktiterationsverfahren (GAUSS-SEIDEL) mit Überrelaxation (Young, 1954, /15/) Voraussetzung : gute Konditionierung des Gleichungssystems (Systemmatrix diagonal dominant, Pecletzahl < 1); verallgemeinerte konjugierte Gradientenverfahren – PCG-Verfahren (Minimierungsaufgabe) (Schmid, Braess 1987 / 16/, Leismann, Frind 1989 /17/) ORTHOMIN-Verfahren (Northrup, Woo 1988 / 18/, Eisenstatt et al. 1988 /19/) Behandlung nichtlinearer Probleme - Dichteeffekte, nichtlineare Adsorptionsisotherme, Reaktionsterme Iterative Lösungen - wenn die Massenbilanz an jedem Knoten einen vorgegebenen max. Fehlerbetrag nicht mehr überschreitet, kann die Iteration abgebrochen werden. Numerische Lösung von Differenzengleichungen für Strömungsprozesse Zur einfacheren Schreibweise definieren wir die Koeffizienten der Differenzengleichung im Element i nach (eindimensional):

Speicherkoeffizient ( )t

zyxngS i

iifFli ΔΔ⋅Δ⋅Δ

⋅⋅+⋅= 2ρκκ

Leitfähigkeitskoeffizient ( )12

1

2

21

21

++

+Δ+Δ⋅

Δ⋅Δ⋅⎟⎟

⎠

⎞⎜⎜⎝

⎛ ⋅⋅=

iii

ixx

zykgKη

ρ

Quell-/Senkenterm zyxmQ iVi Δ⋅Δ⋅Δ⋅=..

33

Man erhält die Differenzengleichung für eindimensionale Strömung in einfacher Form:

( ) ( ) ( ).

1211

21,, iiiiiiitittii QhhKhhKhhS =−⋅−−⋅−−⋅ +

+−

−Δ+ i=2,...,M-1

Daraus kann leicht die Differenzengleichung für die dreidimensionale Strömung abgeleitet werden. Die Aufgabe besteht nun darin, Methoden zur numerischen Lösung des Gleichungssystems zu entwickeln. Zeitliche Diskretisierung Es ist erforderlich zu entscheiden, welchem der beiden Zeitpunkte tn bzw. tn+1 oder welchem Zwischenwert tn ≤ t ≤ tn+1 die Spiegelhöhen h der örtlichen Strömungsterme zugeordnet werden sollen. Man geht dazu vom physikalischen Sachverhalt aus: Massenänderung im Element während des Zeitintervalls ∆tn

( ) .1 mhhS nn =−⋅ +

Massenstrom über eine Elementoberfläche

Dieser Massenstrom ist zeitlich veränderlich im Intervall ∆tn: .

mhK =Δ⋅ Abbildung 1 : Zeitlicher Verlauf des Massenstromes über eine Elementoberfläche

Die Zuordnung des Massenstromes zu einem Zeitpunkt ist willkürlich: zum Zeitpunkt tn (alle Spiegelwerte h bekannt) = explizites Lösungsverfahren zum Zeitpunkt tn+1 (alle Spiegelwerte h unbekannt) = implizites Lösungsverfahren Mittelwert der Massenströme über die Elementoberfläche = semi-implizites Verfahren (tn+1 + tn) / 2 (jeweils halbe Anteile nach CRANK-NICHOLSON)

allgemein : ( )2

.

1

.

3

.1 ⎟

⎠⎞

⎜⎝⎛⋅+⎟

⎠⎞

⎜⎝⎛⋅−=⎟

⎠⎞

⎜⎝⎛

AAA mmm αα

34

Mittelungsfaktor 0 ≤ α ≤1 = 0 explizites Verfahren = 1 implizites Verfahren = 1⁄2 semi-implizites Verfahren Stabilität, Konsistenz, Konvergenz Nach SMITH (1971) setzt sich der Gesamtfehler einer numerischen Lösung zusammen aus: Gesamtfehler = Diskretisierungsfehler + Stabilitätsfehler ( E - N ) = ( E - U ) + ( U - N ) E - exakte Lösung der Differentialgleichung U - numerisch exakte Lösung des Differenzengleichungssystems (unendlich viele Dezimalstellen) N - numerische Lösung des Differenzengleichungssystems (endlich viele Dezimalstellen) Der Stabilitätsfehler ergibt sich aus der summarischen Wirkung aller Rundungsfehler. Stabilitätsbedingungen : (vgl. Abschnitt 2)

Neumann-Bedingung : 21

2 ≤Δ⋅Δ⋅xStD

Courant-Bedingung : 1≤Δ⋅

Δ⋅→

xS

tw

Abbaukriterium : 2≤Δ⋅ tSk

kombinierte Bedingung nach Bear (1979, /13/) bei "upwind"-Wichtung (sie kombiniert o.g. Bedingungen):

22

1*

*2 ZePexStD

++≤

Δ⋅Δ⋅

Pe* - Gitter-Peclet-Zahl, Ze*- Zerfalls- oder Abbauzahl, auf die Gitterlänge eines ortsdiskreten Gitternetzes bezogen: Ze* = k⋅∆x2/ D D - hydrodynamischer Dispersionskoeffizient m2/s , allg. Leitfähigkeit →

w - Geschwindigkeit m / s, S - Speicherkoeffizient, k - Abbaukoeffizient s-1

Pe - Peclet - Zahl, auf die Gesamtlänge bezogen Pe = |w| ⋅L / D,

35

36

L - Länge m Eine stabile Lösung ist konsistent, wenn sie für kleine Orts- und Zeitschritte gegen die Lösung E der Differentialgleichung konvergiert (durch numerische Testrechnungen aufzeigbar). Konvergenz liegt vor, wenn der Diskretisierungsfehler (E-U) einer numerischen Lösung für kleiner werdende Orts-und Zeitschritte gegen Null strebt. Konsistenz + Stabilität =======⇒ Konvergenz (Äquivalenz - Theorem von LAX (1967) /20/) Bei der numerischen Lösung von praktischen Aufgaben ist es in der Regel relativ einfach, den Stabilitätsfehler relativ klein zu halten. Der Diskretisierungsfehler (Konvergenz) muss durch Variation der Zeit- und Ortsschrittweiten abgeschätzt werden; nur eine konvergente Lösung ist physikalisch sinnvoll und praktisch zu verwerten.

Das explizite Lösungsverfahren Mit der Zuordnung der Spiegelwerte h in den Ausdrücken K⋅∆h zum alten Zeitpunkt tn ergibt sich:

( ) ( ) ( ).

1

211

21

1i

ni

ni

n

ii

ni

n

i

ni

nii QhhKhhKhhS =−⋅−−⋅−−⋅ +

+−

−

+ i = 2,...,M-1

einzige Unbekannte ist hi

n+1 , so dass die Differenzengleichung nach hin+1 umgestellt werden

kann.

LRE ZE ZE ZE RRE • • • • • • • hL i=1 i=2 i i=M-1 i=M hR

Abbildung 2 : Elementdiskretisierung (eindimensional) Für die Randelemente i = 1, M kann obige Gleichung leicht geändert werden. Das explizite Lösungsverfahren beginnt mit der Berechnung im Element i = 1 für den ersten unbekannten Zeitpunkt (n+1), dabei sind im Anfangszeitpunkt (n = 1) alle Spiegelhöhen und Quellen durch die Anfangsbedingung gegeben. Die Berechnung schreitet fort für i = 2,..., M-1 und beginnt wieder in i = 1 für den nächsten Zeitschritt. Das Stabilitätskriterium für das explizite Verfahren lautet

212

1

≤+

i

i

S

K i = 2,....., M-1, d.h. für konstante Größen k, n, η, ρ, κ, ∆x

21

2 ≤Δ⋅⋅⋅

Δ⋅xn

tkκη

Die Beziehung zeigt anschaulich, dass das Verfahren zeitschrittbegrenzt ist, wobei ein engmaschiges Gitternetz zu extrem kleinen Zeitschritten zwingt. Numerische Diskretisierung (ohne Quell-/Senkenterm) DD - dimensionslose diskrete Diffusionskoeffizienten ( Ki+1/2 / Si ) Linkes Randelement i = 1 , Rand i = ½ DD = 0.5 mit ∆s1 = ∆x1 / 2 Rechtes Randelement i = iM , Rand bei iM+1/2 DD = 0.5 mit ∆sM = ∆xM / 2 Zentrale Elemente i =2,..., M-1 DD = 0.25 mit ∆si = ∆xi+1...M-1 (∆s1,M = Abstand vom 1. bzw. letzten Element zum Randwert)

37

Daraus ergeben sich die diskreten Rechengleichungen : Linkes Randelement , i = 1 ( )n

LDnn

i hhhDhh 2111 32 +−⋅+=+

Zentrale Elemente , i = 2,...,M-1 ( )n

iiiDni

ni hhhDhh 11

1 2 +−+ +−⋅+=

Rechtes Randelement , i = iM ( )n

RiMiMDniM

niM hhhDhh 231

1 +−⋅+= −+

(diskreter Diffusionskoeffizient DD = ¼ , zentrale Elemente) Das implizite Lösungsverfahren Das implizite Verfahren ordnet die Terme K ⋅∆h dem jeweils neuen, noch unbekanntem Zeitpunkt tn+1 zu, damit ergibt sich (eindimensional):

( ) ( ) ( ) 111

1

21

11

1

21

1 +++

+

+

+−

+

−

+ =−⋅−−⋅−−⋅ nii

ni

n

ii

ni

n

i

ni

nii QhhKhhKhhS i = 2,...,M - 1

Diese Gleichung enthält drei Unbekannte (hi-1

n+1, hin+1,hi+1

n+1), so dass nur die Aufstellung des gesamten Systems für i = 1,...., M eine Lösung erlaubt ! (Randbedingungen notwendig ! ) Das Gesamtsystem besteht aus (M) Gleichungen und enthält (M) unbekannte h-Werte; es ist eindeutig lösbar. Die Koeffizientenmatrix besitzt eine diagonale Bandstruktur, die Zeilen i=2,..., M-1 besitzen je drei Elemente (eindimensional), man nennt diese Matrix "tridiagonal". Für dieses Gleichungssystem existiert eine sehr effektive Form des Gauß'schen Elimina-tionsverfahrens (PROGONKA). Eliminationsverfahren für tridiagonale Gleichungssysteme Das Gleichungssystem lautet abgekürzt: BB1 ⋅ h1 - C1 h2 = D1 i = 1, LRE -Ai ⋅ hi-1 + Bi ⋅ hi - Ci ⋅ hi+1 = Di i = 2..M-1 , ZE - AM ⋅ hM-1 + BM ⋅ hM = DM i = M , RRE Dabei bedeuten die Koeffizienten für i = 2,..., M-1: Ai = Ki-1/2 BBi = Ki-1/2 + Ki+1/2 + Si Ci = Ki+1/2 Di = Qi

n+1 + Si hin

38

Der Lösungsalgorithmus ist: Berechnung schrittweise für i = 2,...M progressiv w1 = B1 wi = Bi - Ai⋅fi-1f1 = C1 ⁄ B1 fi = Ci ⁄ Bi g1 = D1 ⁄ B1 gi = (Di + Ai⋅qi-1) ⁄ Wi ----------------------------------------------------------- Berechnung schrittweise für i = (M-1),..,1 rekursiv hM

n+1 = gM hi

n+1 = gi + fi ⋅ hi+1 Lösung des Differenzengleichungssystems für mehrdimensionale Probleme Für zwei- und dreidimensionale Probleme (flächenhafte und räumliche Strömungsvorgänge) ist die tridiagonale Form des Gleichungssystems nicht mehr gegeben (5 bzw. 7 Unbekannte in jeder Differenzengleichung). Die Lösung solcher Systeme könnte nach dem Gauß'schen Eliminationsverfahren erfolgen. Bei M N L Gitterpunkten in den drei Koordinatenrichtungen besitzt die Matrix des Systems also (M⋅ N⋅ L)2 - Elemente. Das erfordert einen großen Speicherplatzbedarf und lange Rechenzeiten. Eine effektive Lösung bietet das Sparse- Matrix-Verfahren mit "Geordneter Elimination" an, welches die Tatsache ausnutzt, dass nur wenige Elemente der Matrix verschieden von Null sind (Sames, 1974, /21/). Stabilität und Konvergenz des impliziten Verfahrens Die implizite numerische Lösung ist stabil bei allen Zeitschrittweiten ∆t. Um Oszillationen, d.h. physikalisch unsinnige Schwingungen der Lösung zu vermeiden, ist es notwendig, die Zeitschritte zu beschränken (ca. 10...100 fach größere Zeitschritte als beim expliziten Verfahren). Das semi-implizite Lösungsverfahren

( ) ( ) ( ) ( ) ( ) ( )

( ) 1

111

21

111

211

211

21

1

1

1

+

+++

+

++−

−+

+−

−

+

⋅+⋅−=

⎥⎦

⎤⎢⎣

⎡−+−⋅⋅−⎥

⎦

⎤⎢⎣

⎡−⋅+−⋅⋅−−−⋅

ni

ni

ni

ni

i

ni

ni

i

ni

ni

i

ni

ni

i

ni

nii

QQ

hhKhhKhhKhhKhhS

αα

αα

(α = ½, CRANK-NICHOLSON), i = 2,...., M-1 Die Umstellung ergibt ein analoges Gleichungssystem, es gelten die Ausführungen zur impliziten Lösung. Das semi-implizite Verfahren ist stabil für alle Zeitschritte. Der Diskretisierungsfehler ist gegenüber dem impliziten Verfahren geringer. Demgegenüber steht eine stärkere Neigung zu Oszillationen, d.h. in der Regel erfordert es kleinere Zeitschrittweiten als das implizite Verfahren.

39

Diskretisierung der Randmassenströme mL und mR

( )2/

1.

xhhDA

m LsL

Δ−⋅⋅

= ( )

2/

.

xhhDA

m RiMsR

Δ−⋅⋅

=

hL,R - Spiegelhöhen auf dem Rand DS - spezifischer Diffusions-/Leitfähigkeitskoeffizient A = ∆y ⋅ M Für die Randpunkte i = 1 und i = M kann keine Gleichung aufgestellt werden (Vi = 0 , ∆xi = 0 ). Die RB 1. u. 2.Art werden für die "echten" inneren Randelemente i=2 und i=M-1 geschrieben. RB 1.Art werden in Differenzenform exakt erfasst RB 2.Art können nur näherungsweise erfasst werden (keine Strömung zwischen den Randpunkten i=1 und i=2 bzw. M und M-1). Quellen im Strömungsraum - "innere Ränder" Innere Ränder, deren Abmessung bei der diskreten Betrachtungsweise in endlicher Differenzenform nicht dargestellt werden können, sind z.B. Bohrungen in porösen Schichten, deren Durchmesser sehr viel kleiner als die Elementeabmessung ist. Zur näherungsweisen Erfassung von inneren Rändern wird im Element (i, j) stationäre oder quasistationäre radialsymmetrische Strömung vorausgesetzt, so dass ein einfacher Zusammenhang zwischen Druck (Spiegelhöhe, Temperatur, Konzentration) im Zentrum der Quelle und dem Mittelwert des Druckes (Spiegelhöhe, Temperatur, Konzentration) im Element hergestellt werden kann. Die typische Form der Gleichung ist:

FSrr

hhkconstVk

B

E

BB

−+

−⋅⋅=

ln

.

.

V B - Volumenstrom der Bohrung > 0 Injektion, < 0 Produktion h - Spiegelhöhe im Gitterpunkt (Mittelwert des Volumenelementes) hB - Spiegelhöhe im Bohrloch

rE - Radius der Einzugskontur π

yxrEΔ+Δ

=

rB - Radius der Bohrung SK - Skinfaktor SK > 0 Verringerung der Durchlässigkeit SK < 0 Erhöhung der Durchlässigkeit F - Korrekturfaktor für die Art des Strömungszustandes im Element 0 < F < 1 (GW-Strömung = 1 ) const. = 2 π M GW - Strömung gespannt = 2 π hm ungespannt ( hm Mittelwert aus hB und h)

40 = 2 π M / η Druckdarstellung der Flüssigkeitsströmung

2.2 Numerische Verfahren Einfache und erweiterte Bilanzmethode (Elementeverfahren) (Control Volume Method (CVM)) Bei der Bilanzmethode wird das betrachtete Gebiet V in beliebige Bilanzelemente ∆Vi unterteilt und die partielle Differentialgleichung integriert. Differentialgleichung des Transportes :

qtuSkuuwugradDdiv −∂∂

=−− )(

Bilanzgleichung für das i-te Volumenelement :

[ ] ∫ ∫ ∫∫ ∫ ∫ΔΔ

−∂∂

=−−ii VV

dVqtuSdVkuuwugradDdiv )()(

Die Anwendung des Gaußschen Satzes führt zu einem Oberflächenintegral:

∫ ∫ ∫ ∫ ∫ ∫ ∫ ∫Δ Δ Δ

→

−∂∂

=−−i i iA V V

dVqtuSdVkudAnwuugradD )()(

Für den eindimensionalen Transport hat die Bilanzgleichung im Kontrollvolumen ∆Vi die folgende Form:

∫∫Δ+

Δ−

Δ+

Δ−

Δ+

Δ−−

∂∂

Δ=Δ−−∂∂ 2/

2/

2/

2/

2/

2/)()(

xx

xx

xx

xx

xx

xxx

i

i

i

i

i

i

dxqtuSAdxkuAdAuw

xuD

Die partiellen Ableitungen ∂u / ∂x werden durch die zentralen örtlichen Differenzen ersetzt. Für die Approximation der Unbekannten u an den Grenzen zu den Nachbarelementen zur Berücksichtigung des konvektiven Transportes bestehen mehrere Möglichkeiten:

),(),1(),(),2

( 121

21,

txutxutxxu iixiixi +++−+=

Δ+ δδ

41

bzw.

),(),1(),(),2

(211

21,

txutxutxxu iixiixi−

−−

−+=Δ

− δδ

mit δx = ½ zentrale Wichtung = (1+sign(w))/2 upwind Wichtung („upstream/downstream“) = Kombination bzw. gewichtete Extrapolation (Peclet-Zahl-Wichtung, Flux-Limiter, Front-Limitation) Einfache Bilanzmethode : Alle Parameter und die Funktion u werden durch ihren Mittelwert ui im Volumenelement i repräsentiert (Subgebietskollokation). Mit Hilfe numerischer Standardverfahren ist die Berechnung des DGL-Systems möglich. In den meisten Fällen wird das zeitliche Differential jedoch durch den zeitlichen Differenzenquotienten ersetzt, und man bekommt ein lineares System, das mit dem Gauß-Algorithmus gelöst werden kann. Das Gleichungssystem hat folgende allgemeine Form:

11 11 ++ +

Δ=⎟⎟

⎠

⎞⎜⎜⎝

⎛+

Δnn

n

n

n

buSt

uASt S-Speichermatrix

A-Leitfähigkeitsmatrix b-Vektor der Quellen/Senken, RB (n, n+1 = Zeitschrittindex, ∆tn = tn+1 – tn ) Mit den Globalmatrizen A (Speicher- und Leitfähigkeitsterm) sowie b (Quellen/Senken, RB und Vorrat zur Zeit tn) ergibt sich die Standardschreibweise A x = b . Erweiterte Bilanzmethode : z.B. durch genauere Abbildung der Integrale (z.B. Trapezregel) führt auf äquivalente Formen der FEM-Gleichung und gleiche Lösungsverfahren. Für die Bilanzmethode können folgende Empfehlungen gegeben werden: - Strömungsprobleme sollten aufgrund der langen Beobachtungszeiten implizit gelöst werden. Besondere Beschränkungen für die Orts- und die Zeitschrittweite müssen i.a. nicht beachtet werden. Die Gitterpunktanzahl muss in ausreichendem Maße der physikalischen Aufgabenstellung angepasst sein. - Transportprobleme können explizit gelöst werden (Beobachtungszeit liegt in der Größe des Zeitintervalls ∆ts ). Empfohlen wird das semi-implizite Crank-Nicholson Schema mit zentraler zeitlicher und örtlicher Wichtung bei Einhaltung des Stabilitätskriteriums und des Überwiegens des dispersiven Transportes Pe* < 1.

42

43

Der "Front Limitation" Algorithmus - eine neue Extrapolationstechnik für Advektion-Diffusions Probleme Wie im Kapitel 2 aufgezeigt, bereitet die numerische Lösung der Transportgleichung einige Schwierigkeiten. Sie stellt eine Differentialgleichung gemischten Typs (parabolisch - hyperbolisch) dar und macht Beschränkungen in der Wahl der Ortsschritte (PECLET-Zahl) und Zeitschritte (COURANT-Zahl) erforderlich. Werden diese Beschränkungen ignoriert, zeigen die numerischen Lösungen Oszillationen und numerische Dispersion, die zu unbrauchbaren Ergebnissen führen. Der Front Limitation Algorithmus (Häfner et al. /81/) umgeht diese Schwierigkeiten. Er ist in die Bilanzmethode (CVM) integriert und nutzt die Aufspaltung in eine implizite Lösung für den Diffusion-Dispersions Term (einschließlich Quell-Senken Term) sowie eine explizite Lösung für den Advektionsterm. Eine spezielle Extrapolationstechnik verhindert die numerische Dispersion und erhält scharfe Fronten. Die Methode erfüllt die Massenbilanz in jedem Volumenelement und zeigt keinen Gitterorientierungseffekt. Die Beschränkungen in der Gitterpecletzahl und im NEUMANN- Kriterium entfallen. Die Courant-Zahl sollte < 1 sein. Für die Quell-Senken Terme gelten die allgemeinen Beschränkungen. Durch die Aufhebung der Pecletzahlbeschränkung können die Gitterabstände den geologischen Bedingungen angepasst werden und unterliegen keiner zusätzlichen Verfeinerungen durch das numerische Verfahren. Daraus resultiert eine deutliche Einsparung an Rechenzeit und eröffnet neue, äußerst effektive Möglichkeiten bei der Lösung inverser Probleme. Methode der "Front Limitation" Das Hauptproblem der FD- und FE-Methoden (so genannte Matrixmethoden) liegt in der Diskretisierung des Advektionsterms, der die Ausbreitung einer scharfen Front in eine bevor-zugte Richtung beschreibt. Deshalb wurden von verschiedenen Autoren Korrekturverfahren entwickelt ( Flux-limiter ), die den Massenstrom und seine Richtung an den Gittergrenzen durch Interpolationsfunktionen berücksichtigen (Harten 1983, 1988, Michael & Rubin 1991, Rubin & Blunt 1990, Vasiliev 1981, Zalesak 1979, Zazovskij & Musaev 1987). Der Front Limitation Algorithmus unterscheidet sich in folgenden Punkten: - In Advektion-Diffusions Problemen sind die advektiven Prozesse dominant, d.h. die wesentliche Information kommt mit der Front ("upstream") und muss deshalb auf die Ausströmgrenze des Gitterelementes (i+ 1/2) extrapoliert werden. - Das Ergebnis des Transportprozesses ist eine Konzentrationsverteilung, deshalb wird der Konzentrationswert extrapoliert und nicht der Massenstrom (w⋅ C). - Die gewöhnliche upstream-Wichtung ist eine Extrapolation nullter Ordnung. Die zentrale Wichtung führt, wie andere Wichtungen auch, auf Interpolationsfunktionen. Diese Extra-polation auf die Gittergrenze i+1/2 wird im Front Limitation Algorithmus durch einen Poly-nomialfunktion zweiter Ordnung beschrieben:

2)(p xa = (x)C i

i

p

0=iE 1/2+i

≤∑

44Die Polynomialkoeffizienten ai werden durch die Gitterpunkte "hinter der Front " (upstream)

i, i-1 und i-2 bzw. ihren zugehörigen Zustandsgrößen bestimmt. Der Verlauf der Polynomial-funktion muss monoton sein, das kann für p ≤ 2 leicht gezeigt werden. Differenzenschema eines Advektion-Diffusions Problems. Man erkennt, dass das Differ-enzenschema zu einer symmetrischen, diagonal-dominanten Koeffizientenmatrix führt, die mit bekannten Gleichungslösern gelöst werden kann. Tafel 1 : 2D-Differenzenschema und Polynomialextrapolation

Die Extrapolation wird durch folgende Bedingung beschränkt:

))C ),C,C((),C,C(_ = C +i E1+ii1+ii+i 2/12/1 MINMAXMAXMin

Diese Bedingung realisiert, dass die extrapolierten Konzentrationen Ci+1/2 zwischen den Werten Ci und Ci+1 , als Maximum oder Minimum ("Front Limitation") liegen. Zur Erhaltung der Massenbilanz im Element i liefert die Extrapolation im Element i-1 die Konzentration auf der Gittergrenze CE,i-1/2 . Der Front Limitation Algorithmus wird in jede Koordinatenrichtung angewandt. Ist die Monotonizität der Polynomialfunktion nicht gegeben, muss mit nullter

45

oder erster Ordnung extrapoliert werden. Eine ähnliche Methode wurde 1996 von UNGER et al. / / publiziert. Dort wird ein Van Leer Flux-Limiter als Wichtungsfunktion zwischen upstream- und downstream Gitterpunkten vorgeschlagen und in eine vollständig implizite Lösung integriert. Verifizierung des FL – Algorithmus

1.00E-08

1.00E-07

1.00E-06

1.00E-05

1.00E-04

1.00E-03

1.00E-02

1.00E-01

1.00E+00

0 5 10 15 20 25 30 35 40 45 50

C-ana, Pe=1

C-FL, Pe=1

C-ana, Pe=10

C-FL, Pe=10

C-ana, Pe=100

C-FL, Pe=100

C-ana, Pe=1000

C-FL, Pe=1000

conc

entra

tion

x, m

Abbildung: 1D- Stofftransport für Pecletzahlen von 1 bis 1000, halblogarithmisch

0 500 1000 1500 2000 2500 3000 3500 4000 4500 50000

500

1000

1500

2000

2500

3000

3500

4000

4500

5000

x, m

y, m

analytical solutionnumerical solution with FLsemi-implicit, cno=0.5

Abbildung : Flächenhafte Impulsquelle (DIRAC – Quelle) Gitterpeclet = 10, long./ trans.

Dispersivität = 10, Konzentration nach 10 Jahren

46

FINITE – ELEMENT – Methode - flexibel in der räumlichen Auflösung von Rändern - Einführung unregelmäßiger Gitternetze - Vorteile bei der Einführung tensorieller Konzepte Lit.: PINDER & GRAY 1977 /22/ , ZIENKIEWICZ 1971 /23/ FEM - Konzentration wird an jedem Punkt des modellierten Gebietes durch eine Interpolationsfunktion beschrieben (Stützstellen = Knoten), während bei der Differenzen-methode der Funktionswert an den Netzknoten gesucht wird, der repräsentativ für den zugehörigen Abschnitt ∆ x ist. Die exakte Konzentrationsverteilung wird in einem globalen Sinn durch die Interpolationsfunktion mit Basisfunktionen möglichst gut approximiert. Abstandsgeschwindigkeiten, Dispersionskoeffizienten, Mächtigkeit, Abbaurate und effektive Porosität werden als Elementeigenschaften angesehen, d.h. sie sind über das Element konstant. Sollten diese Eigenschaften an den Knoten gegeben sein, dann kann die Elementeigenschaft durch Mittelung über die Knoteneigenschaften gewonnen werden. Die Interpolationsfunktion hat folgende Form:

47

) ( ) ( ) (∑=

=M

jjj xtutxu

0

,ˆ ω

Galerkin – Verfahren : Beim Einsetzen von in den Differentialoperator L (u) wird sich ein Residuum ergeben. Das Galerkin-Verfahren fordert, dass das Residuum orthogonal zu den Basisfunktionen ist:

u

( ) ),(ˆ txuL ε= → Residuum

( ) ( ) ( ) ( )∫ ∫ ===L L

ii MidxxuLdxxtx0 0

,....2,10ˆ, ωωε

Nach Einsetzen der Funktionen und Anwendung des GREENSCHEN Integralsatzes erhält man ein DGL-System zur Bestimmung von uj (t) folgender Form:

( ) ( )∑ ∑ ∑= = =

⎟⎟⎠

⎞⎜⎜⎝

⎛⎟⎟⎠

⎞⎜⎜⎝

⎛−+=+++

M

j

M

j

M

jj

Lj

jiijjijijij

jij uw

dxd

DxqukwDdt

duS

0 0 00ω

ωω

mit

und ∫=L

jiij dxSS0

ωω ∫=L

jiij dx

dxd

dxd

DD0

ωω

und ∫=L

jiij dxkk0

ωω ∫=L

ji

ij dxdx

dww

0

ωω

Lineare Approximation Basisfunktionen = „Hütchenfunktionen“ (Wert = 1 am Knoten i, Wert = 0 an den benachbarten Knoten) bzw. als Elementefunktion nur im Element verschieden von 0. Mit der Definition der Elementefunktionen sowie der Einführung lokaler Koordinaten lassen sich die Koeffizienten (Integrale) der gewöhnlichen DGL einfach bestimmen. Für den Speicherkoeffizienten S ergibt sich z.B. für die eindimensionale Transportgleichung:

( ) jiiiijiiiijiiiij xSxSxSxSS ,11111,1 61

31

61

+++++− Δ+Δ+Δ+Δ= δδδ

SymbolKronecij −= kerδ

In analoger Weise erhält man die anderen Koeffizienten: Dij , wij , kij und qij . Ergänzt wird das Gleichungssystem durch die Auswertung der Randflüsse (Randbedingungen). Für den einfachsten Fall (eindimensional, konstante Parameter, RB 1.Art bei x=0 und L) der Transportgleichung- erhält man folgende Gleichung:

SquZePeuZeuZePe

xSD

dtdu

dtdu

dtdu

iiii

iii

+⎥⎦

⎤⎢⎣

⎡⎟⎠⎞

⎜⎝⎛ −−+⎟

⎠⎞

⎜⎝⎛ +−⎟

⎠⎞

⎜⎝⎛ −+

Δ

=++

+−

+−

1***

1**

2

11

61

211

322

61

211

61

32

61

i=1,2,....., M-1 AB: ui = ui

0 für t = 0 und i = 1,2,...,M-1 RB: ui = u0 für i = 0 und t > 0 ui = uL für i = M und t > 0 Ze*- Zerfalls- oder Abbauzahl, auf die Gitterlänge eines ortsdiskreten Gitternetzes bezogen: Ze* = k⋅∆x2/D Dieses Gleichungssystem kann mit Standardsoftware (Runge-Kutta-Verfahren) direkt gelöst werden. Meist wird jedoch die Lösung analog zum Vorgehen beim finiten Differenzen-verfahren durch Ersetzen der Ableitungen durch zeitliche Differenzen gefunden.

48

Ersetzen der Zeitableitung : Für die Konzentrationen werden gewichtete Mittel zwischen den entsprechenden Werten zur Zeit t und zur Zeit t + ∆t verwendet. Das Gleichungssystem wird für die Konzentrationen ui(t+∆t) bei Kenntnis der Konzentrationen ui(t) gelöst:

mit ( ) ( )

ttuttu

dtdu iii

Δ−Δ+

= ergibt sich in verallgemeinerter Form:

( ) ( ) ( ) ( ) ( ) ( ) ( )( ) iititijii

ij rtuttuBt

tuttuA =⋅−+Δ+⋅⋅+⎟

⎠⎞

⎜⎝⎛

Δ−Δ+

⋅ δδ 1' ( )

i = 1, 2,...M -1 δt = 1 voll implizites System δt = 0.5 CRANK-NICHOLSON-Schema Letztlich kann das Gleichungssystem allgemein mit den Koeffizientenmatrizen A und B dargestellt werden:

ruBudtdA += bzw. du / dt = A-1 B u + A-1 r und

mit den Globalmatrizen A (Speicher- und Leitfähigkeitsterm) sowie b (Quellen/Senken und Vorrat zur Zeit tn) zur Standardschreibweise A x = b übergeht. Im Gegensatz zur Differenzenmethode ist das Gleichungssystem im Falle δt = 0 nicht explizit lösbar. Durch die Wichtung des Speicherterms muss immer ein gekoppeltes Gleichungssystem gelöst werden. Dazu werden in der Literatur folgende numerische Lösungsmethoden genannt:

• Standard – Galerkin Methode (und Modifikationen) • Petrov – Galerkin Methode

SU – Petrov – Galerkin streamline upwind FU – Full upwind technique SC – Shock-capturing technique PGLS – Petrov – Galerkin least square method .........

49

50

Der Speicherplatzbedarf ist über die Bandbreite der Systemmatrix eine Funktion der Knoten-nummerierung. Die Nummerierung ist längs der kürzeren Rechtecksseite am günstigsten (MARSAL, 1976 /60/). Für einen Vergleich zwischen Punktiterationsverfahren und direkter Gleichungslösung bei der Lösung der Transportgleichung nach dem GALERKIN-Verfahren wird auf YEH (1985) /26/ verwiesen. Das GALERKIN-Verfahren ist auf Modelle beliebiger räumlicher Dimensionalität anwendbar. In drei räumlichen Dimensionen wird das Dreieckselement zum Tetraederelement oder zum Dreiecksprismenelement (GREENKORN 1983 /32/). Numerische Gleichungslösungsverfahren beschreiben VEMURI, KARPLUS 1981 /27, STÖSSINGER (1979) /28/. Die Kopplung über den Dichteterm an die Strömungsgleichung führt zu einem dichteabhängigen Ansatz z.B. für Salzprobleme. NEUMANN et al.(1974) /61/ schlagen bei nichtlinearen Problemen vor, die Koeffizienten-matrix aus numerischen Gründen künstlich auf Diagonalform zu bringen. Man spricht dann von " lumped finite elements ". Bei unzureichender Feinheit der Diskretisierung führt das GALERKIN-Verfahren zu Oszillationen. Bei verschwindenden Dispersivitäten wird das Standard - GALERKIN- Verfahren wegen Oszillationen unbrauchbar. Eine Alternative bietet das TAYLOR-GALERKIN-Verfahren (DONEA 1984 /24/, LÖHNER et.al.1984 /25/). Eine zusätzliche upwind-Gewichtung des Speicherterms führt zum PETROV-GALERKIN-Verfahren (BROOKS & HUGHES 1982 /29/). Das verbesserte Verfahren glättet zwar ebenfalls die anfangs scharfe Front durch numerische Dispersion, liefert aber danach den korrekten Transport der geglätteten Verteilung. In Gebieten großer Konzentrationsgradienten ist eine feine Diskretisierung erforderlich. Da die Konzentrationsgradienten von den Dispersivitäten abhängen, bestimmen diese letztlich die Feinheit der Diskretisierung. Vorteilhaft sind variable Elementnetze. Die genauesten Ergebnisse erzielt man bei der FEM, wenn die Strömungsrichtung mit der Richtung einer Dreiecksseite übereinstimmt. Dies führt zur Methode- principal directions - (FRIND, 1982 /30/) im Falle stationärer Strömung. Prinzipiell können auch Interpolationsfunktionen höherer Ordnung eingesetzt werden (KINZELBACH, FRIND 1986 /31/). Während der Übergang von der eindimensionalen zur mehrdimensionalen Betrachtung bei der Finiten Differenzenmethode sehr einfach ist, sind auf Grund der Konstruktion der Lösung bei der FEM und CVM einige Besonderheiten zu beachten (siehe Literatur). Methode der Charakteristiken EVANS & HARLOW 1957 /33/, ABBOTT 1966 /34/ GARDER et al. (1964) /35/ - für Konvektions-Diffusionsgleichung KONIKOW & BREDEHOEFT 1978 /36/ - Transportmodell für Grundwassersektor Grundidee : Kombination aus - konvektiver Teil der Transportgleichung nach Charakteristikenmethode - Stromterm nach FEM oder FDM Bestimmt wird die zeitliche Änderung der Konzentration in einem Kontrollvolumen, das sich

51

längs einer Bahnlinie x(t), y(t) "Charakteristik" bewegt. Die Änderung der Konzentration wird verursacht durch dispersiven Schadstoffaustausch, Reaktionen und Quellen/Senken, die das Kontrollvolumen auf seiner Bahn antrifft. Populäre Form : Bewegung von Tracerteilchen mit einer Anfangskonzentration entlang einer Charakteristik. Übergang von Zeitpunkt t zu Zeitpunkt t+ ∆t in 3 Schritten: 1. Teilchen werden konvektiv entlang ihrer Bahnlinie bewegt. An jedem Knoten des Rechteckgitters wird eine intermediäre Knotenkonzentration als Mittel der Konzentrationen aller Teilchen bestimmt, die sich zu Ende des konvektiven Schrittes in der 'Zelle' befinden. 2. Betrifft nur die Knotenkonzentrationen: dispersiver Austausch mit Nachbarzellen, Q/S, Reaktion und Verdünnung (FDM oder FEM ohne Konvektionsterm). Daraus folgen die neuen Knotenkonzentrationen zur Zeit t+∆t. 3.Aktualisierung der Teilchenkonzentrationen - jedem Teilchen in einer Zelle wird die Änderung der entsprechenden Knotenkonzentration (nicht die mittlere Knotenkonzentration über die gesamte Zelle, was zur Zerstörung scharfer Fronten führen würde) mitgeteilt. Die Massenbilanzen müssen mit den Knotenkonzentrationen berechnet werden. Die Genauigkeit des Modells wird vom Gitterabstand und von der Dichte der Tracerteilchenbelegung (ca. 4 – 9 im Minimum) bestimmt. Der Programmieraufwand ist hoch. Die MOC neigt sehr stark zu Oszillationen. Random - Walk - Verfahren (Verfahren der Zufallsbahnen) CHANDRASEKHAR 1943 /37/, HAKEN 1983 /38/, GARDINER 1983 /39/ Eine für die praktische Anwendung auf die Modellierung der Schadstoffausbreitung in natürlichen Aquiferen geeignete Form ist in Arbeiten von PRICKETT et al. (1981) /40/, UFFINK (1985) /41/ und AHLSTROM & FOOTE (1976) /42/ gegeben. 'particle - tracking - Verfahren' Im Random-Walk-Verfahren sind die Tracerteilchen mit einer festen Schadstoffmasse verknüpft. Die Dispersion wird dadurch modelliert, dass der konvektiven Bewegung der Einzelteilchen eine Zufallsbewegung überlagert wird, deren statistische Eigenschaften den Eigenschaften des dispersiven Prozesses entsprechen. Erst dadurch, dass viele zufallsbehaftete Einzelbahnen (Random-Walks) simuliert werden, entsteht eine dispergierende Tracerpunkteverteilung. Diese stellt eine Schadstoffmassenverteilung dar, deren Gesamtmasse gleichmäßig auf alle eingesetzten Tracerteilchen verteilt ist. Wie in der MOC ist eine Buchhaltung über Teilchenort und -zahl notwendig. Die Qualität der numerischen Resultate hängt wesentlich von der Anzahl der verwendeten Tracerteilchen ab. In Teilen des Modellgebietes mit geringer Teilchenzahl brauchen die zufallsbehafteten Ergebnisse nicht signifikant zu sein. Auch die RWM neigt zu Oszillationen, die durch die Randomisierung maskiert wird. Die Wahl der Gitterabstände bestimmt wesentlich die Güte der Ergebnisse. Konzentrationen

sollten aus mindestens 25 Tracerpunkten berechnet werden. Die RWM bietet als einzige Lösungsmethode die Möglichkeit der Einführung skalenab-hängiger Dispersivitäten und ist eine robuste Methode, die auf jedes Strömungsmodell aufgepfropft werden kann. Sie liefert keine numerische Dispersion im klassischen Sinne. Der benötigte Speicherplatz ist gering. Da zur Berechnung signifikanter Konzentrations-verteilungen viele Teilchen (>1000) notwendig sind, ist der Rechenzeitbedarf groß. Es wird empfohlen, die Methode zur Berechnung von Momenten einer Schadstoffverteilung, Durchbruchskurven und Aufenthaltszeitverteilungen anstelle von Konzentrationsverteilungen zu verwenden KINZELBACH 1985, 1987 /43/ /44/. Die Randintegralgleichungsmethode (Boundary Integral Equation Method) Die BIEM leitet sich aus der Potentialtheorie ab.

Potenzial U ( )( ) ( ) ( )222

,,ψηξ

μ

−+−+−=

zyxzyxU

P(x, y, z), Q(ξ, η, ψ); ( ) ( ) ( )222 ψηξ −+−+−= zyxr Das Potential U ist im gesamten dreidimensionalen Raum mit Ausnahme des Quellpunktes (P) definiert. Für P = Q verschwindet der Nenner, der Ausdruck μ / r ist deshalb für r = 0 nicht definiert. Das Punktpotential U = μ/r ist eine Lösung der Laplaceschen Differentialgleichung.

Abbildung 5 : Integrationsgebiet Die allgemeine Potentialfunktion: Man nennt nun jede Funktion U(x, y, z), die nach allen drei Veränderlichen zweimal stetig differenzierbar ist und die in einem gewissen Gebiet Ω des Raumes die Gleichung ∆U = 0 erfüllt, eine Potentialfunktion oder auch harmonische Funktion in diesem Gebiet. Wählt man U = 1/r , wobei r der Abstand des Aufpunktes von einem festen Punkt Q bezeichnet, dann ist:

01=⎟

⎠⎞

⎜⎝⎛Δ=Δ

rU

außer in P = Q . Diesen Punkt schließt man aus dem Integrationsgebiet Ω zunächst aus, da U dort eine Singularität hat. Die Lösung von U am Punkt U (P, Q) nennt man singuläre Grundlösung. Man bezeichnet diese singuläre Grundlösung auch als Greensche Funktion des unendlichen Raumes. Diese Lösungen sind bekannt. Will man Q in Ω aufnehmen, so hat man

52

rechnerisch einen Grenzübergang durchzuführen. Bezeichnet man mit u die Lösung der inhomogenen selbstadjungierten Differentialgleichung div ( D ⋅ grad u ) - k ⋅u = q und mit v die Lösung der singulären Grundlösung U(P,Q) von div ( D ⋅ grad v ) - k ⋅ v = 0 dann ergibt sich: - nach Multiplikation der Gleichungen mit v bzw. u - Subtraktion beider Gleichungen - Integration über das Gebiet Ω - Anwendung des Gauß`schen Satzes - Einsetzen der singulären Grundlösung v (P, Q) und Grenzübergang für ε→ 0 die Grundgleichung der Randintegralgleichungsmethode:

( ) ( )∫∫ ∫∫∫Γ Ω

⋅⋅−∇⋅−∇⋅⋅=K

dVvqdAuvvuDPU

Neben der singulären Grundlösung v (Greensche Funktion) sind ebenso die Quellen und Senken bekannt. So stellt obige Gleichung eine Beziehung zwischen den Randwerten u und grad u sowie der Funktion u im Inneren des Gebietes her. Man nutzt diese Eigenschaft zur Bestimmung der Potentialfunktion im Inneren aus den Randwerten, indem man den Auf- punkt P in die Nähe des Randes legt. Die BIEM erniedrigt die Dimension des Problems, denn die Funktion u(P) wird nur von den Randwerten UГ und grad uГ bestimmt. Greensche Funktion für die Strömungsgleichung ∆v = 0

eindimensional ξ−= xr r21

zweidimensional ( ) ( )22 ηξ −+−= yxr rln2

1⋅

⋅π

dreidimensional ( ) ( ) ( )222 ψηξ −+−+−= zyxr r⋅⋅

−π41

Die BIEM wird vorwiegend zur Lösung stationärer Prozesse eingesetzt, aber auch instationäre Prozesse können mit der BIEM gelöst werden, wobei sich dabei der Einsatz der Laplace-Transformation mit numerischer Rücktransformation als günstig erweist.

53

54

2.3 Deterministische Modelle und stochastische Modellierung Der sensitivste Parameter eines Transportmodells ist die Abstandsgeschwindigkeit (bestenfalls auf eine Größenordnung genau). Ein zweiter Unsicherheitsfaktor der Transport-modelle sind die Dispersivitäten (Fick’sches Dispersionsgesetz erst im Fernfeld gültig, Dispersivitäten können sich innerhalb von Prognosezeiträumen ändern). Der dritte Unsicher-heitsfaktor ist der durch Intensität und zeitlichen Verlauf charakterisierte Schadstoffeintrag. Die Vielfalt der Interpretationsmöglichkeiten vergleichbarer Signifikanz kann beträchtlich zunehmen, falls chemische Prozesse wie Reaktion und Adsorption eine Rolle spielen. Die Vergleichbarkeit von Modellrechnungen mit Messwerten aus dem Feld ist beschränkt. Ein regionales zweidimensionales Modell versagt, wenn dreidimensionale Effekte wichtig werden. Eine wesentliche Voraussetzung der Transportmodelle ist die stochastische Homoge-nität des Aquifers. Systematische Berücksichtigung der Unsicherheit von Modellen - Einsatz von Modellen unter Verwendung konservativer Annahmen - Sensitivitätsanalysen zur Abschätzung von Unsicherheiten in den Parametern - Stochastische Modellierung als Kombination der Modellierung von Erwartungswert und Streuungsmaß interessierender Größen. Stochastische Modellierung (Gelhar, et al., 1979 /45/; Gelhar, Axness, 1983 /46/, Dagan, 1984 /47/, Matheron, de Marsily, 1980 /48/) – (siehe Abschnitt 2.3.1 „Hydrodynamische Dispersion“) - stochastische Natur der Durchlässigkeitsverteilung als Ursache der Modellunsicherheit Die stochastische Theorie liefert zwei Dinge, die über die deterministische Modellierung hinausgehen: - den stochastischen Zugang zum Makrodispersionsterm - eine zweite Transportgleichung, die die Unsicherheiten des berechneten Erwartungswertes der Konzentration beschreibt Eine vom Konzept her einfachere, wenn auch rechenintensivere Methode der stochastischen Modellierung des Schadstofftransportes ist die Monte - Carlo - Methode (Smith, Schwartz, 1980,1981a, b /49/, /50/, /51/). 2.3.1 Hydrodynamische Dispersion in porösen Medien Stofftransport – deterministische Beschreibung : Die hydrodynamische Dispersion charakterisiert das „Verwischen“ von Konzentrations-fronten und wird analog zu der aus der Brownschen Molekularbewegung resultierenden Diffusion beschrieben.

( ) qtccwcgradDDdiv m −

∂∂

=⎥⎦

⎤⎢⎣

⎡ −⋅⎟⎠⎞

⎜⎝⎛ +

−

(bei konstanter Porosität)

Daraus wird bei eindimensionaler Betrachtung mit linearer Reaktionskinetik:

qtcckcw

xcD

x−

∂∂

=−⎟⎟⎠

⎞⎜⎜⎝

⎛−

∂∂

∂∂

Für w, D = const. ! erhält man die Differentialgleichung:

qtcck

xcw

xcD −

∂∂

=−∂∂

−∂∂

2

2

(Die Lösungen dieser DGL liegen für bestimmte Modellbedingungen den experimentellen Bestimmungen des Dispersionskoeffizienten zugrunde) Für einen punktförmigen, impulsartigen Stoffeintrag bei x=0 und t=0 ergibt sich z.B.:

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡ −−−=

⎟⎠⎞⎜

⎝⎛

tD

twxtk

tDm

vtxc A

4exp

2),,(

2

1 π (w als Parameter)

Empirisch wurde nachgewiesen (BEAR 1972, HULLA, RAVINGER u.a. 1982, BEIMS 1983), dass der hydrodynamische Dispersionskoeffizient D geschwindigkeitsabhängig ist:

νδ wD = 21 ≤≤ν (in guter Näherung = 1) Die Dispersivität δ ist dann eine charakteristische Länge. Weiterhin kann man den Untersuchungen entnehmen (BEIMS 1983), dass offenbar eine Proportionalität zwischen dem Betrachtungsabstand x und dem Faktor δ besteht:

xαδ = (in Häfner, Sames, Voigt : α = 0,017 bei einer Streuung der Geschwindigkeit im Mittel von 20 % ) Die Proportionalität zwischen δ und x hat Konsequenzen, die nicht mit der Praxis vereinbar sind. Es folgt daraus auch eine lineare Abhängigkeit des Dispersionskoeffizienten D von der Raumkoordinate:

xwD α= oder allgemeiner : 0DxwD +=α (D muss als Materialkonstante invariant gegenüber einer Koordinatentransformation sein,

55

Ergebnisse werden abhängig von der Versuchsanordnung und letztlich besteht ein Widerspruch zur Annahme D = const.). Mit obiger Annahme entsteht die Differenzialgleichung:

( )( ) ( ) qtcck

xcw

xcDxw −

∂∂

=−∂∂

−−∂∂

+ αα 12

2

0

deren Lösung nicht mit der oben angegebenen übereinstimmt. 2.3.2 Dispersion unter Beachtung des Einflusses der Varianz von Geschwindigkeitswerten Setzt man in der DGL des Transportes D = 0, so erhält man eine Konvektionsbeschreibung:

qtcck

xcw −

∂∂

=−∂∂

−

Diese hat unter den gleichen Modellbedingungen (AB und RB) für einen Eintrag eines Schadstoffimpulses (Dirac-Impuls) die Lösung:

( ) ⎟⎠⎞

⎜⎝⎛ −⎟

⎠⎞

⎜⎝⎛−=

wxt

wxk

wmvtxc A δexp,,2

Diese Lösung wird superponiert für verschiedene Geschwindigkeiten w mit einer Wahrscheinlichkeitsdichte p (w) für die Geschwindigkeitsverteilung:

( ) 1=∫+∞

∞−

wdwp und für den Mittelwert der Geschwindigkeit :

( )∧+∞

∞−

=∫ wwdwpw

Setzt man näherungsweise für die Geschwindigkeit eine Normalverteilung, ergibt sich:

( )⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡⎟⎠⎞

⎜⎝⎛ −−

=

∧

2

2

2exp

21

s

ww

swp

π mit Erwartungswert und der Varianz s

∧

w 2.

56

Die Superposition der Konvektionslösung mit der Normalverteilung ergibt:

2

221

3 2,,

⎟⎟⎠

⎞⎜⎜⎝

⎛−−−∧

∧

=⎟⎠⎞

⎜⎝⎛ w

tx

stk

A et

es

mwtxcπ

Aus dem einfachen Vergleich der Lösungen c1 = c3 ergibt sich für den Dispersionskoeffizienten D:

2

2 tsD = und aus 03 =∂∂

xc

findet man die allgemeine Tatsache, dass sich

das Maximum der Konzentration c3 mit der Geschwindigkeit ausbreitet, es gilt: ∧

w

( ) twcxx∧

== maxmax daraus ergeben sich durch Substitution:

max2

2

max

2

2.

2x

w

sbzwxw

sD∧∧ == δ

Diese Beziehungen stellen die von BEIMS (1983) gefundene lineare Abhängigkeit des Dispersionskoeffizienten von der Raumkoordinate x dar. Fazit : Eine Normalverteilung der Geschwindigkeit in einem porösen Medium täuscht einen zeitproportionalen bzw. einen der Weglänge x proportionalen Dispersionskoeffizienten vor. Alle Versuche, die Maßstabsabhängigkeit in die hydrodynamischen Modellrechnungen einzubeziehen, führten auch bei den zugrunde gelegten komplizierten stochastischen Modellen zu keinen prinzipiell anderen Ergebnissen. Identifiziert man bei Säulenversuchen mit unterschiedlichen Längen aus den Durchbrüchen eines Tracers die Dispersivität im herkömmlichen Sinne, so erhält man mit der Säulenlänge wachsende Parameter, auch wenn alle Versuchsbedingungen konstant gehalten werden. Wertet man die Versuche nach der superponierten Lösung aus, bestimmt man Varianzen der Geschwindigkeit, die ihrerseits nur noch zufällig streuen.

57

58

2.3.3 Stochastische Prozesse In stochastischen Prozessen werden Zufallsgrößen untersucht, die außer vom Zufall von einer oder mehreren Veränderlichen abhängen. Der stochastische Prozess X (t, w) ist somit die Gesamtheit von Zufallsgrößen, die von einem reellen Parameter t (z. B. Zeit) abhängen – Realisierung des Prozesses. X (t, w) gibt somit bei festem w z.B. die Anzahl von Teilchen oder Geschwindigkeits-vektoren in Abhängigkeit von der Zeit an. Praktische Bedeutung haben solche stochastischen Prozesse X (t, w) bei denen für alle t der Mittelwert m (t) über alle w und die Varianz / Streuung D² X (t, w) von X (t, w) über alle w existieren, mit: m (t) = E [X (t, w)] und D² X (t, w) = E ( [X (t,w) – m (t)]² ) E (...) = Erwartungs- oder Mittelwert (stochastischer Prozess zweiter Ordnung) Stochastische Modellierung (G. DAGAN, 1995) Parameter :

(1) räumlich verteilt : Permeabilität, Porosität, Retardation, Diffusion / Dispersivität, Sorption ...

(2) diskrete Parameter : Spiegelhöhen, Neubildung, Pumpraten, Ausgangskonzentrationen K (x) – hydraulische Leitfähigkeit als RSF (random space function), d. h. als Ensemble von Realisierungen für den geologischen Untergrund. Vereinfachung: Y = ln K, d.h. const. Mittel- oder Erwartungswert < Y > und der Kovarianz CY (x, y), r = x-y (Abstand zwischen den Punkten) bezüglich der Integrationsskala IY (Anisotropie horizontal /vertikal).

Transportproblem C (x, t) Die analytischen und numerischen Lösungen der Transportgleichung führen mit konstanten Dispersionskoeffizienten (Labormessungen) zu ellipsoiden Schadstoffverteilungen, wobei scheinbar die Dispersionskoeffizienten mit dem Abstand oder der Zeit zunehmen. Praktische Untersuchungen zeigen: (1) Die Konzentrationsfront hat eine irreguläre Form, weit weg von einer

Gauß’schen Glockenkurve.

(2) Die longitudinale Streuung ist weit größer als Mikrodispersion (Porenskala). (3) Die vertikale und transversale Streuung ist viel kleiner als die longitudinale, aber

größer als die aus Laborwerten ermittelte.

2.3.4 Die 2 Grundprobleme der Transportmodellierung -Euler’sche und Lagrang’sche Betrachtungsweise- Grundlagen (Strömungsproblem) Statistische Charakterisierung des Fluid-Geschwindigkeitsfeldes : < w (x, t) > = < w > / n (n-Porosität) und 2-Punkt- Kovarianz uij (x, x´, t, t´) = <ui (x, t) uj (x´, t´)> (i, j = 1,2,3) (für ein räumlich stationäres Feld ist uij eine Funktion von r = x-x´, d. h. u = w - <w> bildet die Fluktuationen der Geschwindigkeit ab.) Euler (Volumenelement bezogen, w und Dd als Zufallsgrößen) → physikalisch sinnvoll bei „feiner“ Diskretisierung ! Langrange (Partikel / Punkt bezogen, wobei eine Vielzahl von Partikel durch < C(x,t)> und σ2

C ersetzt werden kann.)

( ) xdtxCnM ∫= , und ( ) xdtxCxnM

R ∫= ,1

bzw. ( ) ( ) ( ) xdtxCRxRxnM

S jjiiij ,1∫ −−=

(n-Porosität, M-totale Masse, R-zentrales Moment (Ortskoordinate der mittleren Geschwindigkeit <w>), Sij – zweites, räumliches Moment, Varianz sij

2) Das Grundwesen der Lagrang’schen Beschreibung liegt nun in der Abbildung der Partikelbahn x = X (t, t0 , a) durch : dX = w [X (t), t] dt + dXd ; t=t0 , X=a (dXd bildet als „Wiener-Prozess (Mengenvereinigung) “ den Diffusionsprozess des Koeffizienten Dd ab.) Nachteil : Bestimmung der lokalen Konzentrationsverteilung für reaktive Transport- modellierung und Mehrphasenströmung wird komplizierter, da die Massenbilanz der Einzelpartikel eingehalten und bestimmt werden muss.

59

2.3.5 Das Konzept der Makrodispersion Analyse der räumlichen Momente (spatial moments) zur Bestimmung der hydrodynamischen Dispersion (GARABEDIAN, LeBLANC, GELHAR & CELIA, 1991). Grundlagen : SCHWARTZ (1977), GELHAR et al. (1979), SMITH & SCHWARTZ (1980), MATHERON & deMARSILY (1980), DAGAN (1982, 1984), GELHAR & AXNESS (1983), NEUMANN et al. (1987) „ longitudinale Dispersivität erreicht asymptotisch einen konstanten Wert in einem stationären, homogenen Medium begrenzter räumlicher Ausdehnung (correlation scale)“ Eine mögliche Alternative wäre durch einen effektiven- oder Makrodispersions-Tensor als Kombination von Geschwindigkeits-Fluktuationen und porenskaliger Dispersion Kef bzw. Def in Verbindung mit einer mittleren Konzentrations-Verteilung gegeben:

CDCwtC 2∇=∇⋅+∂∂

Nachteil : < C > ist statistisch verteilt und keine messbare Größe, sowohl in Euler’scher als auch in Lagrang’scher Betrachtungsweise. (Der einfache Fall einer deterministischen Konzentration mit konstantem, effektivem Makro-Dispersionstensor wird hier nicht betrachtet.) Konzept der räumlichen Momente (spatial moments) als Zufallsvariable der Zeit : <R>, Rij , Sij , VAR (Sij ), ... ARIS (1956) Zeitableitung für das erste räumliche Moment ( -∞ < x < + ∞) :

321 dxdxdxtc

Mx

ntdxd i

∂∂

= ∫∫∫ und 321 dxdxdxcnM ∫∫∫=

n – Porosität, M – totale, gelöste Masse, i = 1, 2, 3 Nach Einsetzen der Transportgleichung und partieller Integration ergibt sich:

ii w

tdxd

= Die Geschwindigkeit ist gleich der Änderung des ersten Moments

mit der Zeit ( FREYBERG, 1986).

60

Zeitableitung des zweiten räumlichen (zentralen) Moments – Varianz sij2 der Orts– bzw.

Raumpunkte ergibt:

( )( )∫ ∫ ∫

∞+

∞− ∂∂−−

= 321

2

dxdxdxtcn

Mxxxx

tdsd jjiiij

Nach Einsetzen der Transportgleichung und partieller Integration ergibt sich:

wAtdsd

D ijij

ij ==2

21

mit Aij als Tensor für die Makrodispersivität (ARIS, 1956,

FISCHER et al., 1979) Bei konstanter Fluidgeschwindigkeit, kann der Tensor für die Dispersivität auch dargestellt werden als:

twxmitxd

sdA ij

ij ==2

21

i, j = 1, 2, 3

Damit ist der tatsächlich wirksame Dispersionskoeffizient Dij einer Konzentrationsverteilung definiert als halbe Änderungsrate des räumlichen Moments. Für den 3-dimensionalen Transport auf lokaler Skala (transversale Änderung der Konzentration ist größer als die transversale Änderung der Geschwindigkeit (Korre-lationsskala) ergibt sich damit:

( ) ( )td

SdDtDStSw

tdRd

tdRd ij

ijijijij 21;; =≅≅==

Lokale / regionale Skala

• Ausbreitungs- und Prognosezeiten von >10 Jahren – keine Messwerte ! • Die lokal erfüllte Forderung nach: „Ergodizität“ (räumliche Ausdehnung der

Konzentrationsverteilung ist groß genug bzgl. der Integrationsskala) der räumlichen Momente bleibt erhalten, ist auf der regionalen (heterogenen) Skala oft nicht erfüllt.

• Damit ist auch die Forderung nach einer konstanten mittleren Geschwindigkeit <w>

nicht mehr erfüllt : (1) Ist oder wird die initiale transversale Erstreckung L der Konzentrations- verteilung kleiner als die Integrationsskala, verkleinert sich auch die

61

62

longitudinale Streuung (Dispersion begründet sich auf Skalen kleiner L, größere Skalen bevorzugen die Ausbreitung der Schadstoffwolke als Ganzes, als zentrale Zufallsbewegung) (2) Änderungen in Sij (Fluktuationen) können in der regionalen Skala auch zur Abnahme führen, woraus letztlich entsprechend Definition ein negativer Dispersionskoeffizient resultiert. Konsequenzen :

1. Die Annahme der Ergodizität in heterogenen, regionalen Skalen führt zu unrealistisch

großen Streuungen, wobei sich das Zentrum mit konstanter mittlerer Geschwindigkeit bewegt.

2. Ist der Schadstoffpfad kleiner als die log-normale Transmissiblitäts-

Integrationsskala, wächst der Dispersionskoeffizient scheinbar mit der Laufzeit – „Superdiffusion“

Fazit : Der Herleitung der Makrodispersion liegt die korrekte Annahme zugrunde, unabhängig vom Ort zu sein, sie enthält die physikalische Modellgrundlage der Zeitabhängigkeit, wobei das Konzept der räumlichen Momente den Dispersionstensor nicht als Materialkonstante, unab-hängig von Ort und Zeit definiert. Auch die stochastische Beschreibung der Makrodispersion bedarf einer Konditionierung und einer Einbeziehung von Messungen, um letztlich eine praktikable deterministische Charakterisierung von Strömungs- und Transportprozessen vornehmen zu können.

63

2.4 Optimale Prozesssteuerung und Parameteridentifikation Optimale Prozesssteuerung Die optimale Prozesssteuerung ermittelt Parameter, die den Prozess in einer erstrebten, optimalen Weise beeinflussen. Grundlage jeder optimalen Steuerung ist die den Prozess beschreibende Zustands- oder Prozessgleichung (Strömungs- und Transportgleichung). Die abhängige Variable zu der Zustandsgleichung (d.h. ihre Lösung) heißt Zustandsgröße. Bei der Randsteuerung lenken die Randbedingungen den Prozess in optimaler Weise. Anfangswertsteuerung - optimale Anfangsbedingungen Quellsteuerung - optimale örtliche/zeitliche Verteilung von Quellstärken Diese Zielstellung wird mathematisch in einer Zielfunktion formuliert. Aus einer Vielzahl von Lösungen oder unendlich vieler Lösungen nimmt die Zielfunktion einen Extremwert an (Minimum oder Maximum) - Ziel- oder Gütefunktional. Das Zielfunktional kann durch Nebenbedingungen weiter spezifiziert werden. Diese NB stellen jedoch keine Extremalbedingungen dar. Man unterscheidet zwei typische Aufgaben der optimalen Steuerung: - optimale Prozesssteuerung, dabei werden Anfangs- und/oder Randbedingungen und/oder Quellen/Senken gesteuert. - optimale Parametersteuerung, dabei werden solche Parameter ermittelt, die einen gemessenen Prozess optimal nachbilden (Parameteridentifikation, Modelleichung). Eine Zielfunktion kann in folgender Form dargestellt werden: J0 = Min J [ ßi ] , mit Nebenbedingungen ßi - Steuergrößen, i = 1, 2 ,..., N Die Zielfunktion ist mathematisch ein Funktional, da sie nur mittelbar über die Lösungen u von den Steuergrößen abhängt. Zur Minimierung des Funktionals haben sich aus einer Vielzahl von Verfahren für die Zwecke der Prozesssteuerung bewährt: - das Verfahren des steilsten Abstiegs (Gradientenverfahren) - das Gauß - Newton - Verfahren - Verfahren nach Levenberg-Marquard - Verfahren nach POWELL - stochastische Suchverfahren Die Verfahren unterscheiden sich im Programmieraufwand, im notwendigen Rechenaufwand und bezüglich der Genauigkeit bei der Annäherung an das Minimum (siehe Tabelle). Das Problem der Eindeutigkeit einer Optimallösung ist identisch mit der Fragestellung, ob das

Zielfunktional einen globalen Extremwert (Max./Min.) besitzt, dessen zugehörige Steuergrößen die eindeutige Optimallösung darstellen. Ebenso können lokale Extremwerte auftreten, die suboptimale Lösungen darstellen. Abbildung 6 : Verlauf eines Zielfunktionals J [ß] beim eindimensionalen Steuerproblem

Für die praktische Lösung von Steueraufgaben ist die Eindeutigkeit im mathematischen Sinne, d.h. die Existenz eines globalen Minimums (Maximums) in der Regel gegeben. Die Feststellung, ob das ermittelte Extremum globaler oder lokaler Art ist, bereitet Schwierig-keiten. Parameteridentifikation - Ermittlung von Kenngrößen von Strömungs- und Transportvorgängen aus einem gemessenen Prozessverlauf. Kenngrößen : Durchlässigkeit Speichervermögen eines Feststoffes Intensität von Quellen oder Senken Abbauzahl Wechselwirkungsparameter eines Fluid - Feststoffsystems Wärmeleitfähigkeit u.a. Dieses Problem kann als optimale Steueraufgabe betrachtet werden, wobei das Zielfunktional die Übereinstimmung von berechnetem und gemessenem Prozessverlauf bewerten muss. Diese Problemstellung wird deshalb oft auch als optimale Parametersteuerung, Koeffizientenanpassung (Fitting) oder Modelleichung bezeichnet.

64

Mathematisch gesehen ist sie eine inverse Aufgabe oder Umkehraufgabe. Als Lösung für die Parameteridentifikation bei vorgegebenem Typ der Prozessgleichung haben sich zwei grund-sätzliche Möglichkeiten bewährt:

65

- die graphisch-analytische Parameterermittlung (Geradenverfahren, Typkurvenverfahren), basierend auf einer analytischen Lösung des jeweiligen Strömungs- bzw. Transport- problems - (Gegenstand der Vorlesung „Stofftransport in porösen Medien“). - die Parameterbestimmung mit Suchverfahren, basierend auf der Methode der Prozess- steuerung (Gradientenverfahren., Gauß- Newton, Levenberg-Marquard, Powell-Verf.,...) - stochastische Suchverfahren 2.5 Programmsysteme zur Simulation von Geofiltration und Geomigration Zur numerischen Simulation von Strömungs- und Transportprozessen existiert eine Vielzahl von Programmsystemen. Die Einsatzmöglichkeiten gehen in fast allen Fällen über die genannten Aufgaben hinaus. Im Rahmen der Lehrveranstaltungen am Institut und zu lösender Forschungsaufgaben erfolgt eine Beschränkung auf folgende Problemkreise: Tagebauentwässerung / Grundwasserwiederanstieg Grundwassersperren / Infiltrationsanlagen Altlasten / Grundwassersanierung / Trinkwasserschutzgebietsbemessung Nachweis von Grundwasservorräten Sicherheitsbewertung von Endlagern Stofftransport im Grundwasser (Schadstoffausbreitung im Untergrund) Wärmeleitung und Wärmetransport (Geothermie) Mehrphasenströmung und Stofftransport (Bodenzone, Kippen, Öl-Gas-Lagerstätten) Am Institut für Bohrtechnik und Fluidbergbau vorhandene und lizenzierte Programmsysteme: Visual Modflow (Modflow/MT3D), MODCALIF, MODAMO, MULTIF / Twophase, PCGEOFIM, GMS, ECLIPSE, ANSYS / FLOTRAN, FLUENT (TU BAF Lizenz) Institutseigene Entwicklungen : MODCALIF (F. Häfner & S. Boy) MODAMO / MODEW (F. Häfner) MULTIF / „Twophase“ (A. Nekrassow) Lizenzierte Programmsysteme: Generisches Modell „MODFLOW“ -Visual Modflow (Modflow/MT3D) MODFLOW – Programmcode (McDonald & Harbaugh, 1988, 1996) 3-dimensionales, implizite finite Differenzen in FORTRAN, Block- bzw. Elementzentriertes Gitter, stationäre / instationäre Strömung, gespannt / ungespannt oder kombiniert,

66

Solver : streng implizite Prozedur (SIP), gleitende sukzessive Überrelaxation (SSOR), konjugierte Gradientmethode (PCCG), direkter Gleichungslöser (Hill, 1990, Harbaugh, 1995), Pre- und Postprocessing Packages (kommerzielle Software – Visual MODFLOW); Programmtool MODPATH (Pollock, 1989, 1994) – semianalytisches Partikel-Tracking –Verfahren. MODFLOWP (Hill, 1992) – Parameteridentifikationstool mit nichtlinearer Regression, Minimierung (gewichtete Fehlerquadratsumme) nach modifizierter GAUSS-NEWTON-Methode oder einer konjugierten Gradienten Methode. MODFLOW (U.S. Geological Survey, USGS) – Appel & Reilly (1994) MT3D (Zheng, 1990) – Advektions-Dispersions-Gleichung (Transportgleichung) MOC (Konikow & Bredehoeft, 1978) – Grundlage für MOCDENSE (Sanford & Konikow, 1985) – dichteabhängiges Strömungsproblem MOC3D (Konikow et al., 1996) PCGEOFIM ist ein eingetragenes Warenzeichen des Ingenieurbüros für Grundwasser GmbH, Leipzig (IBGW ; D. Sames) Aufgaben: Unterirdische Strömung und Transport in porösen Medien. - Filterströmung / Stofftransport / Kalibrierung Algorithmus: Bilanzmethode (CVM) Random Walk Methode (RWM) Parameteridentifikation Programmiersprache: FORTRAN77 GMS – Groundwater Modeling System (www.scientificsoftwaregroup.com)

• Modellentwicklung / Kalibrierung / Post-processing / Visualisierung • Grundwasserströmung / Stofftransport

2D und 3D-Modellierung für FDM und FEM, enthält MODFLOW 2000, MODPATH, MT3DMS/RT3D, SEAM3D, ART3D, UTCHEM, FEMWATER, PEST, UCODE, MODAEM und SEEP2D ECLIPSE - Black-Oil Simulator ECLIPSE (Schlumberger GeoQuest)

• voll impliziter, 3-Phasen, 3-dimensionaler Black-Oil Simulator • entwickelt Anfang der 80-er Jahre,ursprünglich von ECL, später INTERA • seit Mitte der 90-er Jahre Schlumberger GeoQuest • geschrieben in Fortran 77 nach ANSI Standard, verfügbare Plattformen:

– SUN SPARCstation (OS Solaris), IBM RS/6000 (OS AIX) – SGI Silicon Graphics, Pentium PCs mit Windows NT

• Keywordgesteuert (mit jedem herkömmlichen Editor zu schreiben)

67

ANSYS Swanson Analysis Systems, INC. (SASI), 1970 Aufgaben: Spannungs- und Verformungsverhalten fester Materialien in Abhängigkeit von der Temperatur (linear-elastisches und inelastisches zeitabh. Materialverhalten bei Temperaturbeanspruchung) Strömung, Stoff- Wärme-Transport in Fluiden / Gasen in porösen Gesteinen Algorithmus: Finite-Element-Methode (FEM) Programmiersprache: FORTRAN77 FLOTRAN Compuflo Inc., 1992 (in Verbindung mit Pre-und Postprocessor ANSYS) Aufgabe: Strömung, Stoff- und Wärmetransport in Fluiden Laminare oder turbulente Strömung oder Wärmeübertragung kompressibler und inkompressibler Medien Algorithmus: Finite-Element-Methode (FEM) Programmiersprache: FORTRAN77 FLUENT (FLUENT - Deutschland) Aufgaben : Komplexe Strömungsvorgänge (laminar, turbulent) , Energie - und Stofftransport, chemische Reaktionen, Zweiphasen-Strömung mit Phasenwechsel Algorithmus : Finite - Volumen - Methode Programmiersprache : FORTRAN77

68

Literaturverzeichnis /1/ Häfner,F.; Sames,D.;Voigt,H.-D.: Wärme-und Stofftransport. Springer-Verlag Berlin Heidelberg New York 1992, ISBN 3-540-54665-0 /2/ Kinzelbach,W.: Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser. 2.Auflage, Schriftenreihe gwf Wasser/Abwasser,R. Oldenbourg,Verlag München Wien 1992, Bd.21, ISBN 3-486-26347-1 /3/ Anlauf,R.; Kersebaum,K.Chr.;Liu Ya ping u.a. (eingeleitet und koordiniert von Richter,J.): Modelle für Prozesse im Boden. Programme und Übungen. Ferdinand Enke Verlag Stuttgart 1988, ISBN 3-432-96441-2 /4/ Busch,K.-F.; Luckner,L.: Geohydraulik. Verlag für Grundstoffindustrie Leipzig 1972 /5/ Luckner,L.; Schestakow,W.: Geofiltration. Verlag für Grundstoffindustrie Leipzig 1975 /6/ Luckner,L.; Schestakow,W.:Migrationsprozesse im Boden- und Grundwasserbereich. VEB Deutscher Verlag für Grundstoffindustrie Leipzig 1985 /7/ Busch,K.-F.; Luckner,L.;Tiemer,K.: Geohydraulik. Lehrbuch der Hydrogeologie Band3, 3. neubearbeitete Auflage, Gebrüder Bornträger Berlin Stuttgart 1993, ISBN 3-443-01004-0 /8/ Heeg,W.: Einführung in die Geoströmungstechnik. Lehrbrief 1 : Grundlagen. Bergakademie Freiberg, Sektion Geotechnik Bergbau 1988, BA 0909 01 0 /9/ Heeg,W.: Einführung in die Geoströmungstechnik. Lehrbrief 2: Numerisch diskrete Verfahren für Geoströmungsgleichungen. Bergakademie Freiberg, Sektion Geotechnik und Bergbau, Fachrichtung Bohrtechnik und Fluidbergbau 1987,BA 001702020 /10/ van Genuchten, M.T.; Alves,W.J.: Analytical solutions of the one-dimensional convective-dispersive solute transportequation. U.S.Department of Agriculture Techn.Bull.No.1661,1982 /11/ Tautz, H.: Wärmeleitung und Temperaturausgleich. Akademie-Verlag, Berlin 1971 /12/ Sudicky, E.A.: The Laplace transform Galerkin technique : A time-continuos finite element theory and application to mass transport in groundwater. In: Contaminant Transport in Groundwater, Kobus,H., Kinzelbach, W.(Hrsg.), Balkema, Rotterdam, 1989 /13/ Bear, J.: Hydraulics of groundwater. McGraw Hill Series in Water Resources and Enviromental Engineering, 1979

69

/14/ Peaceman, D.W.; Rachford,H.H.Jr.: The numerical solution of parabolic and elliptical difference equations. Jour. Soc. Industrial and Applied Mathematics, 3 (11), 1955 /15/ Young, D.M.: Iterative methods for solving partialdifferential equations of elliptic type .Trans.Amer.Math.Soc.,76, 1954 /16/ Schmid, G.; Braess,D.: Comparison of fast equation solvers for groundwaterflow problems. In: Grounwater Flow and Quality Modelling. Custodio et al. (Hrsg.), NATO ASI-Series C, Band 224, D.Reidel, Dordrecht, 1987 /17/ Leismann, H.M.; Frind,E.O.: A symmetric-matrix time integration scheme for the efficient solution of advection dispersion problems. Water Resources Research, 25(6), 1989 /18/ Northrup, E.J.; Woo,P.T.: Application of preconditioned conjugate-gradient-like methods in reservoir engineering. Soc. Pet. Eng. Res.J., 1988 /19/ Eisenstatt, S.C.; Elman,H.C.; Schultz,M.H.: Block-preconditioned conjugate-gradient- like-methods for numerical reservoir simulation. Soc. Pet. Eng. Res.J., 1988 /20/ Lax (1967) in: Richtmyer, R.D.; Morton,K.W.: Difference methods for initial-value problems. Wiley-Interscience, New York 1967 /21/ Sames,D.: Programmpaket "Geordnete Elimination". Inst. für Energetik, Leipzig 1974 /22/ Pinder,G.F.; Gray,W.G.: Finite element simulation in surface and subsurface hydrology. Academic Press, New York 1977 /23/ Zienkiewicz, O.C.: The finite element method in engineering science. McGraw-Hill, London 1971 /24/ Donea, J.: A Taylor-Galerkin method for convective transport problems. Int.J.Num.Meth.Eng.,20,1984 /25/ Löhner, R.; Morgan,K.; Zienkiewicz,O.C.: The solution of non-linear hyperbolic equation systems by the finite element method. Institute for Numerical Methods in Engineering, University College of Swansea, Sweansea, U.K. 1984 /26/ Yeh, G.T.: Comparisons of successive iteration and direct methods to solve finite element equations of aquifer contaminant transport. Water Resources Research, 21(3), 1985 27/ Vemuri,V.; Karplus,W.J.: Digital computer treatment of partial differential equations. Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1981

70

/28/ Stössinger,W.: Beschreibung der hydrodynamischen Dispersion mit der Methode der finiten Elemente am Beispiel des instationären Interface zwischen Süß-und Salzwasser in Grundwasserleitern. Mitteilungen Institut für Wasserbau und Wasserwirtschaft, RWTH Aachen, Heft 28, 1979 /29/ Brooks,A.N.; Hughes,T.J.R.: Streamline upwind / Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier- Stokes-equations. Comp.Meth. in Applied Mech. and Eng.,32, 1982 /30/ Frind, E.O.. The principal direction technique: A new approach to groundwater contaminant transport modelling. In: Finite Elements in Water Resources. Proceedings of the 4th International Conference, Hannover, Springer-Verlag, Berlin 1982 /31/ Kinzelbach,W.; Frind,E.O.: Accuracy criteria for advection dispersion models. In: Finite Elements in Water Resources. Proceedings of the 6th International Conference, Lisboa, Springer-Verlag, Berlin 1986 /32/ Greenkorn,R.A.: Flow phenomena in porous media. Marcel Dekker, New York 1983 /33/ Evans, M.E.; Harlow, F.H.: The particle-in-cell method for hydrodynamic calculations. Los Alamos Scientific Lab., Bericht Nr.LA-2139,Los Alamos, New Mexico 1957 /34/ Abbott, M.B.: An introduction to the method of characteristics. American Elsevier, New York 1966 /35/ Garder, A.O.; Peaceman, D.W.;Pozzi,A.L.: Numerical calculation of multidimensional miscible displacement by the method of characteristics. Soc. Petr. Eng. J., 4(1), 1964 /36/ Konikow, L.F.; Bredehoeft,J.D.: Computer model of two-dimensional solute transport and dispersion in groundwater. Techniques of Water Resources Investigations of the United States Geological Survey, Book 7, Chapter C2. United States Government Printing Office, Washington, 1978 /37/ Chandrasekhar, S.: Stochastic problems in physics and astronomy. Rev. of Modern Physics, 15(1), 1943 /38/ Haken, H.: Advanced synergetics. Springer-Verlag, Berlin 1983 /39/ Gardiner, C.W.: Handbook of stochastic methods. Springer Serie Synergetik, 13, Springer-Verlag, Berlin 1983 /40/ Prickett, T.A.; Naymik,T.G.; Lonnquist,C.G.: A "random walk" solute transport model for selected groundwater quality evaluations. Illinois State Water Survey, Bulletin 65, 1981

71

/41/ Uffink,G.J.M.: A random walk method for the simulation of macro-dispersion in a stratified aquifer. In: Relation of Groundwater Quantity and Quality, Proc. IUGG General Assembly, IAHS-Publication Nr.146, Hamburg 1985 /42/ Ahlstrom,S.W.; Foote, H.P.: Transport modeling in the environment using the discrete-parcel-random-walk-approach. Proceedings of the EPA-conference on Environmental Modeling and Simulation, 1976 /43/ Kinzelbach,W.: Modelling of the transport of chlorinated hydrocarbon solvents in groundwater: A case study. Water Science and Technology, 17, 1985 /44/ Kinzelbach,W.: The random walk method in pollutant transport simulation. In: Custodio et al.(Hrsg.) Groundwater flow and Modelling, NATO ASI Series C, Band 224, Reidel, Dordrecht, 1987 /45/ Gelhar,L.W.; Gutjahr,A.L.; Naff,R.L.: Stochastic analysis of macrodispersion in a stratified aquifer. Water Resources Research, 15(6), 1979 /46/ Gelhar,L.W.; Axness, C.L.: Three-dimensional stochastic analysis of macrodispersion in aquifers. Water Resources Research, 19(1), 1983 /47/ Dagan,G.: Solute transport in heterogeneous porous formations. J. Fluid Mech., 145, 1984 /48/ Matheron,G.; de Marsily,G.: Is transport in porous media always diffusive ? A counter-example. Water Resources Research, 16(5), 1980 /49/ Smith,L.; Schwartz, F.W.: Mass transport, 1. A stochastic analysis of macroscopic dispersion. Water Resources Research, 16(2), 1980 /50/ Smith,L.; Schwartz, F.W.: Mass transport, 2. Analysis of uncertainty in prediction. Water Resources Research, 17(2), 1981a /51/ Smith,L.; Schwartz, F.W.: Mass transport, 3. Role of hydraulic conductivity data in prediction. Water Resources Research, 17(5), 1981b /52/ Press,W.H.; Flannery,B.P.; Teukolsky,S.A.; Vetterling;W.T.: Numerical Recipes. Cambridge University Press, Cambridge 1986 /53/ Häfner,F.;Voigt,H.D.; Bamberg, H.F.; Lauterbach;M.: Geohydrodynamische Erkundung von Erdöl-, Erdgas- und Grundwasserlagerstätten. WTI- Wiss.Tech. Information des Zentralen Geologischen Institutes 26(1985) /54/ Autorenkollektiv: Katalog strömungsmechanischer Simulationsmodelle. Bergakademie Freiberg, Sektion Geotechnik und Bergbau, 1980

72

/55/ Hughes, T.R.: The Finite Element Method, Prentice Hall, 1987 /56/ Zienkiewicz,O.C.; Taylor, R.L.: The Finite Element Method . 4th edition, Vol.1, Basic Formulation and Linear Problems McGraw-Hill, London, 1989 /57/ Kiraly, L.: FEM301 - A Three Dimensional Model for Groundwater Flow Simulation, NAGRA, Technischer Bericht, 1985 /58/ Kröhn, K.-P.: Simulation von Transportvorgängen im klüftigen Gestein mit der Methode der Finiten Elemente. Bericht Nr.30/1991, Institut für Strömungsmechanik und Elektron. Rechnen im Bauwesen, Universität Hannover, 1991 /59/ Kinzelbach, W.: Groundwater Modelling - An Introduction with Sample Programs in BASIC. Developments in Water Science, 25, Elsevier Science Publishers, Amsterdam, 1986 /60/ Marsal, D.: Die numerische Lösung partieller Differentialgleichungen Bibliographisches Institut, Mannheim, 1976 /61/ Neuman, S.P., Feddes, R.A., Bresler,E.,: Finite-element simulationof flow in saturated- unsaturated soils considering water uptake by plants. In: Developments of Methods, Tools and Solutions for Unsaturated Flow, Third Annual Report (Part 1), Technion Israel Institute of Technology, Haifa, 1974 /62/ Crank, J.H.: Mathematics of diffusion. Oxford University Press, London, 1956 /63/ Crank, J.H., Nicholson, P.: A practical method for numerical integration of solution of partial differential equations of heat-conduction type. Proc. Cambridge Philos. Soc. 43, 1947