Prolog

rm(list=ls()) # alle evtl. verbliebenen Variablen aus dem Arbeitsspeicher loeschen

Aufgabe 1a

Wir suchen die Loesung fuer die selbe Situation wie in Aufgabenblatt 1. Es handelt sich um einen ungespannten Grundwasserleiter. In der Vorlesung hatten wir allerdings die numerische Loesung fuer den gespannten Grundwasserleiter hergeleitet. Wie schon vorher bei der analytischen Loesung diskutiert, unterscheidet sich die ungespannte Loesung von der gespannten nur marginal. Die ungespannte eindimesionalte Gleichung lautet:

\[ K^* \cdot \frac{\partial^2(h^*)}{\partial{x}^2} = j_{w,top} \]

mit \(K^* = \frac{K}{2}\) und \(h^* = h^2\).

Es muss also das folgende Gleichungssystem geloest werden:

\[\frac{K^*}{\Delta x^2} \begin{bmatrix} -2 & 1 & 0 & 0 &\dots & 0 \\ 1 & -2 & 1 & 0 &\dots & 0 \\ 0 & 1 & -2 & 1 &\dots & 0 \\ \vdots & & & & \ddots & \vdots \\ 0 & 0 & 0 &\dots & 1 & -2 \end{bmatrix} \cdot \begin{bmatrix} h_2^* \\ h_3^* \\ h_4^* \\ \vdots \\ h_{n-1}^* \end{bmatrix} = \begin{bmatrix} j_{w,top} - H_o \cdot \frac{K^*}{\Delta x^2} \\ j_{w,top} \\ j_{w,top} \\ \vdots \\ j_{w,top} - H_L \cdot \frac{K^*}{\Delta x^2} \end{bmatrix} \]

Aufgabe 1b: Numerische Loesung implementieren

Parameter, Geometrie, Randbedingungen festlegen

# Aquifereigenschaften und diffuse Zufl??sse
K = 5e-5 # 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

# Geometrie initialisieren
L = 420 # m -> aus der Karte, dies entspricht x1 aus der Vorlesung
d = 1 # m Tiefe des Flie??querschnitts
deltax=20
x=seq(0, 420, by=deltax) # m Array, der die x-Koordinate darstellt
Num=L/deltax+1 # Anzahl der Gitterzellen

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

# Transformierte Variablen erstellen
Ho = ho^2 # transformiert RB Ostrand
HL = hw^2 # transformierte RB Westrand
Ktran=K/2 # transformierte Leitfaehigkeit, siehe oben

Jetzt wird die Systemmatrix und der Vektor b auf der rechten Seite erstellt. Fuer die Erstellung der Systemmatrix muessen die Diagonalen und Nebendiagonalen belegt werden.

# Systemmatrix belegen
S=diag(-2,(Num-2)) # Initialisieren einer quadratischen Matrix, die auf der Diagonalen mit -2 belegt ist
diag(S[,-1]) <- matrix(1,Num-3) # obere Nebendiagonale hinzufuegen
diag(S[-1,]) <- matrix(1,Num-3) # untere Nebendiagonale hinzufuegen
S=S*(Ktran/deltax^2) # Jetzt noch diese Matrix mit dem Faktor versehen.

# Vektor b auf der rechten Seite (Spaltenvektor!) belegen
b=matrix(j_wn,Num-2) # Spaltenvektor initialisieren, er enthaelt jetzt ueberall j_wn
b[1]=b[1]-Ho*(Ktran)/deltax^2 # Randbedingung links
b[Num-2] = b[Num-2]-HL*(Ktran)/deltax^2 # Randbedingung rechts

# Vektor h mit Grundwasserhoehen initialisieren und Raender einsetzen
h_num=matrix(,Num) # initialisiert einen Vektor mit leeren Elementen, der Vektor schliesst die Randzellen mit ein. Aktive Zellen sind nur j=2..Num-1.
h_num[1]=ho # linker Rand (wird nicht berechnet, siehe unten)
h_num[Num]=hw # rechter Rand (wird nicht berechnet, siehe unten)

Jetzt das Gleichungssystem loesen lassen:

h_num[2:(Num-1),1]=sqrt(solve(S,b)) # Aktive Zellen sind nur j=2..Num-1

Das Ergebnis gemeinsam mit der analytischen Loesung (aus dem letzten Aufgabenblatt) plotten.

plot(x,h_num+h_Aqbot, ylab="h in m NN")

# Vergeleich mit analytischer L??sung aus Aufgabenblatt 1
xa=1:420
h_star=j_wn/K*xa^2+((hw^2-ho^2)/L-j_wn*L/K)*xa + ho^2 # analytische Loesung
h_ana=sqrt(h_star) # das ist die Ruecksubstitution

lines(xa,h_ana+h_Aqbot) # f??gt eine Linie zu einem bereits existierenden plot hinzu

Aufgabe 1c: Einen Brunnen bauen

Hierfuer muss der Vektor b nur einen zusaetzlichen Senkenterm an irgendeiner Stelle erhalten. Ich waehle 300 m, also \(\Delta x \cdot 15\). Ich waehle weiterhin, dass die Brunnenentnahme 20 mal so hoch sei wie die Grundwasserneubildung sein soll. Dieser Faktor kann spaeter veraendert werden.

Jo=-20*j_wn # Ich ziehe Wasser aus 15 benachbarten Zellen

Der Vektor b muss neu belegt werden, die Systemmatrix bleibt unveraendert.

b2 = b
b2[15] = b[15]+Jo

Ein neues h muss initialisiert werden:

h_num_br=matrix(0,Num)
h_num_br[1]=ho # linker Rand 
h_num_br[Num]=hw # rechter Rand

Gleichungssystem loesen

h_num_br[2:(Num-1),1]=sqrt(solve(S,b2)) # Aktive Zellen sind nur j=2..Num-1

Plotten der numerischen und analytischen Loesung zum Vergleich

plot(x,h_num_br+h_Aqbot, col="red", ylab="h in m NN")
lines(x,h_num+h_Aqbot) # analytische Loesung

Stellen Sie fest: Das Wasser fliesst nur aus dem linken Rand in den Brunnen!