Algorithm for water hammer using MOC

Problem: a reservoir at the upstream end, a pipe of length $L$ with a valve at the downstream end that closes in a given time.

Algorithm

1. Input data: $L$, $D$, wave speed $a$, friction factor $f$, reservoir head $H_{res}$, initial steady discharge $Q_0$, valve closure law $\tau(t)$ (relative opening), total simulation time $T$.
2. Grid: choose the number of reaches $N$; set $\Delta x=L/N$ and $\Delta t=\Delta x/a$ (characteristics along the diagonal).
3. Constants: $A_p=\pi D^2/4$, $B_p=a/(gA_p)$, $R_p=f\Delta x/(2gDA_p^2)$.
4. Initial conditions ($t=0$): $Q_i=Q_0$ at all nodes $i=1,\dots,N+1$; $H_i=H_{res}-R_pQ_0^2(i-1)$ (steady friction loss along the pipe).
5. Time loop: set $t=t+\Delta t$.
6. Interior nodes $i=2,\dots,N$:
   • $K_1=H_{i-1}+B_pQ_{i-1}-R_pQ_{i-1}|Q_{i-1}|$
   • $K_2=H_{i+1}-B_pQ_{i+1}+R_pQ_{i+1}|Q_{i+1}|$
   • $Q_{P_i}=(K_1-K_2)/(2B_p)$, $H_{P_i}=(K_1+K_2)/2$
7. Upstream node ($i=1$, reservoir): $H_{P_1}=H_{res}$; $K_2=H_2-B_pQ_2+R_pQ_2|Q_2|$; $Q_{P_1}=(H_{P_1}-K_2)/B_p$.
8. Downstream node ($i=N+1$, valve): $K_1=H_N+B_pQ_N-R_pQ_N|Q_N|$. With the valve law $Q_P=Q_0\tau\sqrt{H_P/H_0}$ the solution is $$Q_P=-B_pC_v+\sqrt{(B_pC_v)^2+2C_vK_1},\qquad C_v=\frac{(Q_0\tau)^2}{2H_0},\qquad H_P=K_1-B_pQ_P$$ For complete closure ($\tau=0$): $Q_P=0$, $H_P=K_1$.
9. Update: replace old values by the new ones, $Q_i\leftarrow Q_{P_i}$, $H_i\leftarrow H_{P_i}$.
10. Store/print $H_i$ and $Q_i$ (especially the maximum $H$ at the valve and mid-pipe).
11. If $t<T$ go to step 5; otherwise stop. Plot $H$ against $t$ to read the maximum pressure rise.

 Start -> read data -> grid, B, R -> initial Q, H
   -> t = t + dt
   -> interior nodes (C+ and C-)
   -> reservoir node (C-) , valve node (C+ and valve law)
   -> update Q, H -> output
   -> t < T ?  yes: back to "t = t + dt"   no: stop