In [4]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

Assignment 4 Key

Consider flow over an isolated obstacle in a channel, as in class. Assume that the obstacle is a triangle, with height off a flat channel of $h_m = 10\ \mathrm{m}$, and that the triangle's ramp has a slope of 1/100. Assume that the incoming two-d flow transport is a fixed $30\ \mathrm{m^2\,s^{-1}}$, that the flow is in steady state, and that the flow is controlled at the crest.

Q1

What is the thickness of the water column, $d_0$, far upstream of the obstacle? (OK to use a root finder and give a numeric answer)

In [5]:
# Soe variables:

hm = 10 # m
q = 30 # m^2/s
g = 9.8
dhdx = 0.01 # i.e the bump is 1000 m long, 10 m hight

A1

If transcritical then $$ \frac{1}{2} F_0^2 + 1 -\frac{3}{2}F_0^{2/3} = \frac{h_m}{d_0}$$

or in terms of $q$ we get a cubic for $d_o$:

$$d_0^3 - \left(\frac{3}{2}\frac{q^{2/3}}{g^{1/3}} - h_m\right) d_0^2 + \frac{q^2}{2g} = 0 $$
In [7]:
coef = [1., -3/2*q**(2/3)/g**(1/3) - hm, 0, q**2/2/g]
roots = np.roots(coef)
print(roots)
print(q**2/g/roots**3)
d0 = max(roots)
print('d0 = ',d0, 'm')

[ 16.60090986   1.7485285   -1.58191022]
[  2.00733917e-02   1.71790161e+01  -2.31990895e+01]
d0 =  16.6009098571 m

Note that there are two positive roots to this equation: one is for very thin very fast flow that is just transcritical, and one is for very thick slow flow. Lets choose the latter.

So $d_0 = 16.6 \ \mathrm{m}$

Q2:

Knowing the height far upstream you can numerically integrate in x (or calculate the cubic at each point in x) to get the water thickness $d(x)$ as a function of $x$. Plot the water thickness as a function of $x$ from $x=-1000 m$ (upstream) to $x=0\ \mathrm{m}$ (the crest). Check that $F_m = 1$ in your calculation. For precision, make sure that you have a data point every 10 cm or so; include your code, and the expressions you used to get the interface heights. Make the plot as nice as possible (including the obstacle, helps.

A2:

Now lets integrate:

$$ \frac{\partial d}{\partial x} = \left(F^2 - 1 \right)^{-1}\frac{\partial h}{\partial x}$$
In [8]:
dx = 0.01
x = np.linspace(-1000.,0., int(1000/dx))

# init array
d = 0. * x
Fsq = 0. * x
h = 0. * x
d[0] = d0
Fsq[0] = q**2 / g / d0**3
h[0] = 0
for i in range(1,len(x)):
    d[i] = dhdx/(Fsq[i-1] - 1)*dx + d[i-1]
    Fsq[i] = q**2 / g / d[i]**3
    h[i] = dhdx * dx + h[i-1]

u = np.sqrt(Fsq*g*d)

print(h[-1])
um = (g*q)**(1/3)
print('u', u[-1], um)
dm = q / um
print('dm', d[-1], dm)

9.99989999999
u 6.60569378079 6.649399761150975
dm 4.54153658882 4.511685426897414
In [11]:
fig,ax = plt.subplots()
ax.plot(x, Fsq, label='$F^2$')
ax.plot(x, h, label='h')
ax.plot(x, h+d, label='h+d')
ax.set_xlabel('X [m]')
ax.set_ylabel('z [m]')
ax.legend()
print(Fsq[-1])

0.980410562441

Q3

Knowing the water thickness as a function of $x$, use the Momentum Theorem to numerically demonstrate that the momentum balance is satisfied.

A3

Evaluating the momentum theorem for this volume we get that the net flux of moentum is equal to the net forces on the volume:

$$ u_0^2d_0 + \frac{g}{2}d_0^2 - u_m^2 d_m - \frac{g}{2} d_m^2 - g \frac{dh}{dx}\int_0^L d \ \mathrm{dx} = 0$$

where the last term is due to the form drag of the ramp.

In [14]:
print(u[-1]**2 * d[-1])
print(g/2*d[0]**2 )
print(g/2*d[-1]**2 )
print('LHS:')
print(u[0]**2 * d[0] + g/2 * d[0]**2 - u[-1]**2 * d[-1] - g/2 * d[-1]**2) 
print('RHS:')
print(g*dhdx*np.trapz(d, x))

198.170813424
1350.39201961
101.065217479
LHS:
1105.36988468
RHS:
1105.39092161

So the left hand and right hand sides are equal, so the momentum theorem is satisfied.

In [ ]: