Dieses Skript gibt eine graphische Darstellung der analytischen Loesung eines einfachen Grundwasserstroemungsproblems mit folgenden Eigenschaften:

Prolog

# Prolog
rm(list=ls()) # alle evtl. verbliebenen Variablen aus dem Arbeitsspeicher loeschen
dev.off() # entfernt alle offenen Plots 
## null device 
##           1
# Variablen definieren:
K = 1e-4 # m/s -> aus der Tabelle im Vorlesungsskript
j_wn=-0.5/(3600*24*1000) # m/s = 0.5 mm/d, negativ, weil Fluss nach unten

Bemerkung zu den gegebenen Potentialen an den Raendern:

Die Potentiale an den Raendern sind gegeben in “m ueber Normalnull”. Wir koennen die Referenzhoehe jedoch nicht auf NN belassen, da die Potentialhoehe h im ungespannten GWL auch eine Information ??ber die GW Maechtigkeit (also den durchstroemten Querschnitt) enthaelt. Letzterer betraegt nur zwischen 6.5 m - 8 m und nicht 118.5 m - 110 m.

Deswegen gibt es keine andere Moeglichkeit, als die Referenzhoehe auf die Grundwassersohle zu verlegen. Dies muessen wir breits bei der Festlegung der Randbedingungen beachten.

# Randbedingungen:
h_Aqbot=102 # m ueber NN -> Lage des Aquifersohle
ho=110-h_Aqbot # m ueber GW-Sohle -> h am Ostrand
hw=108.5-h_Aqbot # m ueber GW-Sohle -> h am Westrand = h(x=420 m)


# Geometrie initialisieren
L = 420 # m -> aus der Karte, dies entspricht x1 aus der Vorlesung
d = 1 # m Tiefe des Flie??querschnitts
deltax=1
x=seq(0, 420, by=deltax) # m Array, der die x-Koordinate darstellt (deltax ist automatisch ein m)

Aufgabe 2a - Druckgleichung

Die Druckgleichung fuer dieses Stroemungsproblem hatten wir in der Vorlesung hergeleitet. Die Gleichung muss hier nur eingegeben werden. Beachten Sie die Substitution: \(h^{*} = h^2\) fuer den ungespannten Grundwasserleiter.

h_star=j_wn/K*x^2+((hw^2-ho^2)/L-j_wn*L/K)*x + ho^2
h=sqrt(h_star) # das ist die Ruecksubstitution

# Ergebnisse graphisch darstellen
plot(x,h+h_Aqbot)

Aufgabe 2b Berechnung des Flusses

Die Berechnung des Flusses erfolgt nun unter Beruecksichtigung der Potentialgradienten. Der Fluss (\(J_x\) in m\(^3\) s\(^{-1}\)) berechnet sich aus der Fliessgeschwindigkeit (\(j_x\) in m s\(^{-1}\)) und dem Fliessquerschnitt (\(A= m \cdot d = h(x) \cdot d\) in m\(^2\)): \[ J_x = A \cdot j_x \] Dabei ist d (in m) die Breite eines Streifens Aquifer in y-Richtung. Setzt man f??r die Fliessgeschwindigkeit das Darcy gesetz ein, so ergibt sich fuer den Fluss: \[ J_x = -K \cdot h(x) \frac{dh}{dx} \approx -K \cdot h(x) \frac{\Delta h}{\Delta x} \]

Da fuer die Berechnung des Gradienten jeweils zwei benachbarte Potentiale benoetigt werden, erfolgt die Berechnung zwischen den Stuetzstellen (gegeben im Array x), also an der Stelle i+0.5. Praktischerweise erzeugt man dafuer eine neue x-Koordinate (siehe x1 unten). Weiterhin wird auch die durchflossene Maechtigkeit (m=h(x)) an den Stellen i+0.5 benoetigt, ist aber nur an den Stellen i bekannt. Hierfuer gibt es mehere Moegleichkeiten:

Ich habe die zweite Option implementiert:

# Festlegen einer weiteren x-Achse, welche die Punkte zwischen den St??tzstellen 
x1=seq(deltax/2, L-deltax/2, by=deltax)

# Die Maechtigkeit ist der durchflossene Querschnitt. Diesen erhalten wir als lineare 
#Interpoltation zwischen den beiden Stuetzstellen
h1=rowMeans(cbind(h[1:length(h)-1], h[2:length(h)]))

# Jetzt folgt die Berechnung des Flusses entlang des Potentialgradienten
Jx=-K*h1*diff(h)/diff(x) # Darcy-Fluss in m^3/s ->
# j_x = K*m*d*dh/dx, m entspricht der Grundwasserh??he h, also m=h
# Fluss wird an den Stellen x(i+0.5) berechnet

plot(x1, Jx)

Im Ergebnis erhalten wir eine lineare Funktion. Das war zu erwarten, denn schliesslich handelt es sich um die erste Ableitung einer quadratischen Gleichung, nur dass diese hier durch Differenzenbildung und nicht analytisch herbeigefuehrt wurde.

Aufgabe 2c - Berechnung der Flussdivergenz

Die Flussdivergenz berechnet sich nun aus den Differenz von Zu- und Abfluss an den Stuetzstellen. Dies ergibt je einen Wert in der Mitte zwischen zwei Stuetzstellen auf x1. Daher wird wiederum eine neue x-Koordinate benoetigt.

x2 = seq(deltax, L-deltax, by=deltax)

Div_Jx=diff(Jx)/diff(x1)

plot(x2, Div_Jx, xlim = c(0,420), xlab="x [m]", ylab = "div(J_x) [m^3 s^-1]")

Der Plot sieht aus wie ein Sternenhimmel .. Allerdings zeigt ein Blick auf die y-Achse, dass die Werte ueber einen sehr kleinen (\(<10^{-15}\)) Bereich um den Wert \(j_{wn}\) streuen. Das Ergebnis ist also eine Konstante und der Wert ist aus der ursprueenglichen Stroemungsgleichung ablesbar - er war also zu erwarten.

Die Streung kommt durch Rundungsfehler bei der Berechnung zustande. Diese sind vernachlaessigbar klein und werden nur duch sehr geringe Aufloesung der y-Achse, welche automatisch gewaehlt wurde, sichtbar.

Zusatz: Den Plot verschoenern:

Hier verwende ich das Kommando “layout”, welches es gestattet Graphen auf dem selben Blatt sehr flexibel anzuordnen. Um in den Achsen- und Titelbezeichnungen auch mathematische Zeichen verwenden zu koennen (z.B. Indices und Potenzen), verwende ich die Funktion “expression”. Beide sind zu Beginn gewoehnungsbeduertfig, aber sehr hilfreich. Links mit Starthilfen gibt es auf der Seite zur R-Starthilfe auf der Vorlesungsseite.

layout(matrix(c(1,2,3), 3, 1, byrow = TRUE))
plot(x,h+h_Aqbot, xlim = c(0,420), main="Verlauf der Potentialhoehen" , xlab="x [m]" , ylab=expression(h==h[p]+h[GWS] ~ "[m NN]")) # Die Tilde erlaubt Kombination mit normalem Text
plot(x1, Jx, xlim = c(0,420), main="Verlauf des Grundwasserflusses" , xlab="x [m]", ylab = expression(J[x] ~ "[" ~ m^3 ~ s^{-1} ~ "]"))
plot(x2, Div_Jx, xlim = c(0,420), main="Verlauf der Divergenz des Grundwasserflusses" ,ylim = c(-j_wn*0.99, -j_wn*1.01) , xlab="x [m]", ylab = expression("div(" ~J[x] ~ ")  [" ~ m^3 ~ m^{-1} ~ s^{-1} ~"]"))