IHE, Delft, 2017-12-20
@T.N.Olsthoorn
The analytical solution is given by (see syllabus)
$$s = A e^{-a x} \sin(\omega t - a x + \theta)$$Give the envelopes between which the head will fluctuate due to only this wavy head fluctuations at $x=0$
Show the wave in the aquifer at a number of times using the following variable values:
$kD = 500\, m2/d$, the aquifer's transmissivity
$S = 0.001$, the aquifer's storage coefficient
$A = 2 \, m$, the wave's amplitude
$\omega = 1/\pi \, d^{-1}$, the wave's frequency
$\theta = 0.2/\pi$, the wave's initial time delay
In [1]:
import numpy as np
import matplotlib.pyplot as plt
In [2]:
def newfig(title='?', xlabel='?', ylabel='?', xlim=None, ylim=None,
xscale='linear', yscale='linear', size_inches=(14, 8)):
'''Setup a new axis for plotting'''
fig, ax = plt.subplots()
fig.set_size_inches(size_inches)
ax.set_title(title)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.set_xscale(xscale)
ax.set_yscale(yscale)
if xlim is not None: ax.set_xlim(xlim)
if ylim is not None: ax.set_ylim(ylim)
ax.grid(True)
return ax
In [4]:
# define the variables and their values
Tc = 0.5 # [d] cycle time
omega = 2 * np.pi / Tc # [1/d] angle velocity randians/day
theta = 0.2 * np.pi # [-] initial time delay in radians
S = 0.01 # [-] storage coefficient of aquifer
kD = 1000 # m2/d, transmissivity
A = 2 # m, sine ampletude
a = np.sqrt(omega / 2 * S / kD) # 1/m, damping and delay factor
x = np.linspace(0, 500, 501) # m, choose values to compute and show
# Top and bottom envelopes
env1 = +A * np.exp(-a * x) # top envelope
env2 = -A * np.exp(-a * x) # bottom envelop
# Setting up a nice plot
ax = newfig('Sinusoidal fluctation at x=0', '$x$ [m]', '$s = [h - h_0]$ [m]')
# plot the envelopes
ax.plot(x, env1, 'b', label='top envelope')
ax.plot(x, env2, 'b', label='botom envelope')
# choose values for time
time = np.arange(0, 0.5, 0.05)
# do for each time, plot a wave
for ti in time:
s = A * np.exp(-a * x) * np.sin(omega * ti - a * x + theta)
ax.plot(x, s, label=f't = {ti:.2f}')
ax.legend(loc='right')
Out[4]: