Methode der Finite-Differenzen

Mit dem analytischen Strömungsmodell können bereits aussagefähige Modelle gerechnet werden. Es setzt jedoch konstante Parameter voraus. Um beispielsweise auch unterschiedliche Durchlässigkeiten zu berücksichtigen, sind numerische Verfahren notwendig. Ein leicht verständliches Verfahren ist die Methode der Finite-Differenzen, die hier am Beispiel der eindimensionalen Strömung erläutert wird.

Dreidimensionale Kontinuitätsbedingung

Die Kontinitätsbedingung wird anhand der Durchströmung eines Würfels erläutert (Abb. 1).

Abb. 1: Einheitswürfel.

Bei einer stationären Strömung ist die Summe der Ein- und ausfließenden Grundwassermengen gleich Null:

$$Q_{x1}+Q_{x2}+Q_{y1}+Q_{y2}+Q_{z1}+Q_{z2}=0$$

Die Strömung in x-Richtung ist

$$Q_{x1}=q_{x1}\cdot \mathrm{\Delta y}\cdot \mathrm{\Delta z}\text{\hspace{1em}\hspace{1em}\hspace{1em}und\hspace{1em}\hspace{1em}\hspace{1em}}{Q}_{x2}=q_{x2}\cdot \mathrm{\Delta y}\cdot \mathrm{\Delta z}$$ $$\Rightarrow Q_x=(q_{x2}-q_{x1})\cdot \mathrm{\Delta y}\cdot \mathrm{\Delta z}={\mathrm{\Delta q}}_x\cdot \mathrm{\Delta y}\cdot \mathrm{\Delta z}$$

wobei q der spezifische Durchfluss für eine Fläche von 1 m₂ ist. Die Bilanzgleichung mit dem spezifischen Durchfluss

$${\mathrm{\Delta q}}_x\cdot \mathrm{\Delta y}\cdot \mathrm{\Delta z}+{\mathrm{\Delta q}}_y\cdot \mathrm{\Delta x}\cdot \mathrm{\Delta z}+{\mathrm{\Delta q}}_z\cdot \mathrm{\Delta x}\cdot \mathrm{\Delta y}=0$$

wird durch das Volumen V = Δx · Δy · Δz · geteilt und man erhält die dreidimensionale Kontinuitätsbedingung:

$$\frac{{\mathrm{\Delta q}}_x}{\mathrm{\Delta x}}+\frac{{\mathrm{\Delta q}}_y}{\mathrm{\Delta y}}+\frac{{\mathrm{\Delta q}}_z}{\mathrm{\Delta z}}=0$$

Für eine eindimensionale Betrachtung ist die Strömung in y- und z-Richtung gleich Null und die Gleichung reduziert sich:

$$\frac{{\mathrm{\Delta q}}_x}{\mathrm{\Delta x}}=0$$

Laplace-Gleichung

Die Grundwasserströmung basiert auf dem Gesetz von Darcy und der Kontinuitätsbedingung, die zusammengeführt werden.

Gesetzt von Darcy		Kontinuitätsbedingung
$$q=k_f\cdot \frac{\mathrm{\Delta h}}{\mathrm{\Delta l}}$$		$$\frac{\mathrm{\Delta q}}{\mathrm{\Delta x}}=0$$
↘		↙
	$$\frac{\Delta (k_f\cdot \frac{\mathrm{\Delta h}}{\mathrm{\Delta l}})}{\mathrm{\Delta x}}=0$$

Die Zusammenführung der Darcy-Gleichung und der Kontinuitätsbedingung ist die Laplace-Gleichung. Es handelt sich um eine Differentialgleichung, die durch Diskretisierung gelöst wird, wobei der spezifische Durchfluss auf der diskreten Länge Δx gemittelt wird.

Diskretisierung

Ein zweidimensionales Modell wird in Form eines Rechteck-Gitters diskretisiert, wie es bereits im Kapitel „Analytisches Strömungsmodell“ verwendet wurde.

In Abb. 2 ist die Diskretisierung eines eindimensionalen Aquifers dargestellt, der in drei Abschnitte unterteilt wird. Für jeden Abschnitt wird ein spezieller k_f-Wert definiert. Die Grundwasserhöhen h_i werden für vier Stützstellen x_i berechnet.

Abb. 2: Eindimensionaler gespannter Grundwasserleiter mit drei k_f-Wert-Zonen.

Gleichungssystem

Die Grundwasserhöhen an den Stellen x₁ und x₄ seien bekannt. Für die Stellen x₂ und x₃ werden die k_f-Werte und Grundwasserhöhen in die Laplace-Gleichung eingesetzt (Tab. 1).

Tab. 1: Laplace-Gleichung an den Stutzstellen x₂ und x₃

x2	x3
$$\frac{k_{f2}\cdot \frac{h_3-h_2}{x_3-x_2}-k_{f1}\cdot \frac{h_2-h_1}{x_2-x_1}}{x_3-x_1}=0$$	$$\frac{k_{f3}\cdot \frac{h_4-h_3}{x_4-x_3}-k_{f2}\cdot \frac{h_3-h_2}{x_3-x_2}}{x_4-x_2}=0$$

Damit wird ein Gleichungssystem aufgestellt, das aus zwei Gleichungen mit zwei Unbekannten besteht. Durch Auflösen nach h₂ und h₃ können die Grundwasserhöhen für die Stützstellen x₂ und x₃ berechnet werden. Jedoch ist die Lösung des Gleichungssystems rechenintensiv. Daher wird das Gleichungs­system vereinfachend für konstante k_f-Werte und Gitterabstände gelöst.

Prinzipiell kann der Gitterabstand beliebig sein. Hier wird ein konstanter Gitterabstand Δx gewählt,

$$k_{f2}{h}_3-h_2(k_{f2}-k_{f1})+k_{f1}{h}_1=0k_{f3}{h}_4-h_3(k_{f3}-k_{f2})+k_{f2}{h}_2=0$$

so dass das Gleichungssystem leicht nach h₂ und h₃ aufgelöst werden kann.

$$h_2=\frac{k_{f2}{h}_3+k_{f1}{h}_1}{k_{f1}+k_{f2}}{h}_3=\frac{k_{f3}{h}_4+k_{f2}{h}_2}{k_{f2}+k_{f3}}$$

Das Gleichungssystem wird weiter vereinfacht, indem eine konstante Durchlässigkeit vorgegeben wird:

$$h_2=\frac{h_3+h_1}{2}{h}_3=\frac{h_4+h_3}{2}$$

Beispielsweise für bekannte Grundwasserhöhen von h₁ = 13 m und h₄ = 10 m an den Modellrändern sind die übrigen Grundwasserhöhen h₂ und h₃

$$h_2=\frac{h_4+2h_1}{3}=\frac{10+2\cdot 13}{3}=12\text{ m}$$ $$h_3=\frac{h_1+2h_4}{3}=\frac{13+2\cdot 10}{3}=11\text{ m}$$

An den Modellrändern x₁ und x₄ wurden bekannte Grundwasserhöhen vorgegeben. Stattdessen kann auch ein Grundwasserfluss definiert werden: Beispielsweise ist der spezifische Grundwasserfluss für den Knoten x₁

$$q=k_{f1}\cdot \frac{\mathrm{\Delta h}}{\mathrm{\Delta x}}=\frac{h_2-h_1}{\mathrm{\Delta x}}\Rightarrow h_1=h_2-\frac{q}{k_{f1}}\mathrm{\Delta x}$$

In diesem Fall wird für die Stützstelle x₁ eine weitere Gleichung aufgestellt, so dass sich das Gleichungssystem auf drei Gleichungen mit drei Unbekannten erweitert. Für k_f1 = 10^-3 m/s, Δx = 100 m und q = 10^-5 m/s ist die Grundwasserhöhe h₂:

$$h_2=h_4-2\cdot \frac{q}{k_{f1}}\cdot \mathrm{\Delta x}=10-2\cdot \frac{{10}^{-5}}{{10}^{-3}}\cdot 100=12\text{ m}$$

In beiden Fällen handelt es sich um sogenannte Randbedingungen, die bereits im Kapitel „Grundwasserströmung (Darcy 1856)“ anhand der Durchströmung einer Säule erläutert wurden.

Gleichungslöser

Wenige Gleichungen können noch durch Umstellen aufgelöst werden. In der Regel werden umfangreiche Gleichungssysteme aufgestellt, die mit numerischen Verfahren gelöst werden. Grundsätzlich wird zwischen einem direkten und einem iterativen Gleichungslöser unterschieden:

Der direkte Gleichungslöser löst das Gleichungssystem exakt. Er wird bei großen Matrizen ungenau, da die Genauigkeit des Datentyps überschritten werden kann, so dass Rundungsfehler die Folge sind.

Der iterative Gleichungslöser verbessert sukzessiv vorgegebene Startparameter, bis je nach Abbruchkriterium die Iteration beendet wird. Der einfachste iterative Gleichungslöser wurde bereits im Kapitel „Grundwasserabsenkung in einer Baugrube“ angewendet. Es handelt sich dabei um die Jacobi-Iteration, die die Werte aus der vorherigen Iteration übernimmt. Verbesserte Iterationsverfahren nutzten die Differenzen zwischen dem letzten und dem aktuellen Schritt, die zusätzlich mit einem Wichtungsfaktor multipliziert werden.

 

Beispiel: Lateraler Fazieswechsel

Der Aquifer weist entlang der Fließrichtung einen lateralen Fazieswechsel Sand-Kies-Sand auf. Der Parameter k_f und die Randbedingungen, die als Festpotentiale vorgegeben werden, sind Abb. 3 zu entnehmen, wobei Δx = 100 m ist.

Abb. 3: Diskretisierung und Modellvorgaben für einen lateralen Fazieswechsel.

Die Grundwasserhöhen werden mit der Jacobi-Iteration ermittelt. Für die Grundwasserhöhen werden Startwerte von 10 m vorgegeben. Nach 40 Iterationen betragen die Differenzen zur vorherigen Berechnung weniger als 1 cm, so dass die Iteration abgebrochen wird (Abb. 4):

Abb. 4: Berechnete Grundwasserhöhe im Profil.

Bei der Erhöhung der Durchlässigkeit um den Faktor 10 nimmt der Gradient ebenfalls um den Faktor 10 ab:

$$I_{\mathrm{1,2}}=\frac{h_1-h_2}{\mathrm{\Delta x}}=\frac{\mathrm{13,00}-\mathrm{11,57}}{100}=\mathrm{1,4}\cdot {10}^{-2}$$ $$I_{\mathrm{2,3}}=\frac{h_2-h_3}{\mathrm{\Delta x}}=\frac{\mathrm{11,57}-\mathrm{11,43}}{100}=\mathrm{1,4}\cdot {10}^{-3}$$ $$I_{\mathrm{3,4}}=\frac{h_3-h_4}{\mathrm{\Delta x}}=\frac{\mathrm{11,43}-\mathrm{10,00}}{100}=\mathrm{1,4}\cdot {10}^{-2}$$