Part 1: Solving an Initial Value Problem (Falling Head Test)#
Groundwater Flow and the Hydraulic Head Test#
The falling head test is a test to determine the hydraulic conductivity[1] of a soil sample. The soil sample is put in a box with water. A glass tube in the sample is used to measure the hydraulic head \(h\). Before the test, water is poured in the tube while the valve is closed. At the start of the test, the valve is opened and the decrease of the hydraulic head is measured over time.
The hydraulic head \(h\) in this test can be defined by the following ordinary differential equation:
where
\(h\) = hydraulic head [length] (measured in the tube)
\(K\) = hydraulic conductivity [length/time]
\(A\) = cross section area of the sample [length²]
\(a\) = cross section area of the tube [length²]
\(L\) = sample height [length]
Source: A. Verruijt, Grondmechanica / Soil Mechanics; https://geo.verruijt.net/software/SoilMechBook2012.pdf (to be used, copied and distributed without restriction)
In this notebook we will solve the decreasing hydraulic head (“falling head”) \(h\) in time (\(t\)).
True solution:#
The true solution can be obtained by solving the differential equation analytically: $\(h = h_0 \, \exp \left( \frac{-K A t}{a L} \right) \)$
\(\text{Task 1.1:}\)
Run the cells below to define the parameters and visualise the analytical solution:
import math
import numpy as np
import matplotlib.pyplot as plt
K = 1e-6 # hydraulic conductivity
A = 0.1 # cross section area of the sample
a = 0.001 # cross section area of the tube
L = 0.1 # sample height
h0 = 1.0 # initial groundwater head
t0 = 0.0 # start time
maxtime = 5000
dt = 500
n_steps = round(maxtime / dt)
t = np.linspace(0, maxtime, 1001)
h = [h0 * math.exp(-K * A * t[i] / (a * L)) for i in range(1001)]
plt.plot(t, h, label="Analytical solution")
plt.xlabel("time [s]")
plt.ylabel("h [m]")
plt.legend()
plt.show()
Euler Forward:#
\(\text{Task 1.2:}\)
Formulate the Euler Forward solution and complete the code below:
\(\text{Solution 1.2:}\)
Euler Forward formulation:
t_n = np.linspace(0, maxtime, n_steps + 1)
h_EF = []
h_EF.append(h0)
for i in range(n_steps):
h_new = (1.0 - dt * K * A / (a * L)) * h_EF[-1]
h_EF.append(h_new)
Euler Backward:#
\(\text{Task 1.3:}\)
Formulate the Euler Backward solution and complete the code below:
\(\text{Solution 1.3:}\)
Euler Backward formulation:
h_EB = []
h_EB.append(h0)
for i in range(n_steps):
h_new = h_EB[-1] / (1.0 + dt * K * A / (a * L))
h_EB.append(h_new)
Heun’s Method (= 2-stage Runge-Kutta):#
\(\text{Task 1.4:}\)
Formulate Heun’s solution and complete the code below:
\(\text{Solution 1.4:}\)
Heun’s method:
Hence:
h_RK2 = []
h_RK2.append(h0)
for i in range(n_steps):
fac = -K * A / (a * L)
k1 = fac * h_RK2[-1]
k2 = fac * (h_RK2[-1] + dt * k1)
h_new = h_RK2[-1] + dt / 2 * (k1 + k2)
h_RK2.append(h_new)
4-stage Runge-Kutta:#
\(\text{Task 1.5:}\)
Formulate the 4-stage Runge-Kutta solution and complete the code below:
\(\text{Solution 1.5:}\)
Runge-Kutta’s method:
h_RK4 = []
h_RK4.append(h0)
for i in range(n_steps):
fac = -K * A / (a * L)
k1 = fac * h_RK4[-1]
k2 = fac * (h_RK4[-1] + dt / 2 * k1)
k3 = fac * (h_RK4[-1] + dt / 2 * k2)
k4 = fac * (h_RK4[-1] + dt * k3)
h_new = h_RK4[-1] + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
h_RK4.append(h_new)
Plot results#
\(\text{Task 1.6:}\)
Run the code cell below to show the different solutions:
plt.plot(t, h, label="True solution")
plt.plot(t_n, h_EF, label=f"Explicit Euler (dt={dt})", marker=".")
plt.plot(t_n, h_EB, label=f"Implicit Euler (dt={dt})", marker=".")
plt.plot(t_n, h_RK2, label=f"Heun's method (dt={dt})", marker=".")
plt.plot(t_n, h_RK4, label=f"Runge-Kutta (dt={dt})", marker=".")
plt.xlabel("time [s]")
plt.ylabel("h [m]")
plt.legend()
plt.show()
\(\text{Task 1.7:}\)
Comment on the various results in view of the true solution and their numerical stability.
\(\text{Solution 1.7:}\)
Euler Forward solution is unstable for \(dt > 2000\) and predicts \(h\) to decrease too quickly. Note that for \(dt > 1000\) the solution already starts oscillating and becomes non-physical (showing negative \(h\)-values).
Euler Backward solution is unconditionally stable but predicts \(h\) to decrease too slowly.
Heun’s solution is stable and more accurate than the Euler solutions; still, \(h\) decreases a bit too slowly.
Runge-Kutta’s solution is stable and much more accurate than the other solutions.
By Ronald Brinkgreve, Anna Störiko, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.