Algorithm for unsteady pipe flow FD equations

Given: the two finite difference (characteristic) equations on a rectangular grid, with $\Delta t=\Delta x/a$:

$$H_{P_i}=H_{i-1}-B\left(Q_{P_i}-Q_{i-1}\right)-R\,Q_{i-1}|Q_{i-1}|\qquad(C^+)$$ $$H_{P_i}=H_{i+1}+B\left(Q_{P_i}-Q_{i+1}\right)+R\,Q_{i+1}|Q_{i+1}|\qquad(C^-)$$

with $B=\dfrac{a}{gA_p}$ and $R=\dfrac{f\Delta x}{2gDA_p^2}$.

Solution of the two equations at an interior node

Equate the two expressions for $H_{P_i}$ and solve for $Q_{P_i}$:

$$Q_{P_i}=\frac{H_{i-1}-H_{i+1}+B\left(Q_{i-1}+Q_{i+1}\right)-R\left(Q_{i-1}|Q_{i-1}|+Q_{i+1}|Q_{i+1}|\right)}{2B}$$

Then $H_{P_i}$ follows by substituting $Q_{P_i}$ in either equation. Equivalently, $H_{P_i}=\tfrac12\left(K_1+K_2\right)$, with $K_1=H_{i-1}+BQ_{i-1}-RQ_{i-1}|Q_{i-1}|$ and $K_2=H_{i+1}-BQ_{i+1}+RQ_{i+1}|Q_{i+1}|$.

Algorithm

1. Read $L$, $D$, $a$, $f$, the number of reaches $N$, the total time $T$, the initial $Q_i^0$ and $H_i^0$ for $i=1,\dots,N+1$, and the boundary conditions (upstream and downstream).
2. Compute $\Delta x=L/N$, $\Delta t=\Delta x/a$, $A_p=\pi D^2/4$, $B$ and $R$.
3. Set $t=0$ and store the initial values.
4. Repeat while $t<T$:
   1. $t=t+\Delta t$.
   2. For $i=2,\dots,N$ (interior nodes):
      • $Q_{P_i}=\dfrac{(H_{i-1}-H_{i+1})+B(Q_{i-1}+Q_{i+1})-R(Q_{i-1}|Q_{i-1}|+Q_{i+1}|Q_{i+1}|)}{2B}$
      • $H_{P_i}=H_{i-1}-B(Q_{P_i}-Q_{i-1})-RQ_{i-1}|Q_{i-1}|$
   3. Upstream boundary ($i=1$): only the $C^-$ equation is available. Use it with the boundary condition (e.g. a reservoir, $H_{P_1}=H_{res}$):
      • $Q_{P_1}=\dfrac{H_{P_1}-H_2}{B}+Q_2-\dfrac{R}{B}Q_2|Q_2|$
   4. Downstream boundary ($i=N+1$): only the $C^+$ equation is available: $H_{P}=H_N-B(Q_P-Q_N)-RQ_N|Q_N|$, used together with the boundary condition (for a closed end $Q_P=0$, so $H_P=H_N+BQ_N-RQ_N|Q_N|$; otherwise a valve, pump or reservoir relation).
   5. Update: $Q_i\leftarrow Q_{P_i}$, $H_i\leftarrow H_{P_i}$ for all $i$.
   6. Output $H_i$ and $Q_i$ at the required nodes and times.
5. Stop when $t\ge T$.

 read data; compute dx, dt, B, R; set initial Q, H; t = 0
 while t < T:
     t = t + dt
     for i = 2..N:   Q_P(i) from the C+ / C- pair, then H_P(i)
     upstream node:    C-  + boundary condition
     downstream node:  C+  + boundary condition
     Q(i) = Q_P(i);  H(i) = H_P(i)   for all i
     print / store results
 end while

The time step must satisfy $\Delta t=\Delta x/a$ (Courant number 1) so that the characteristics pass through the grid points and no interpolation is needed.