Water hammer: rectangular grid and specified intervals

Water hammer is the pressure surge produced by a rapid change of velocity (valve closure, pump trip). It is solved with the characteristic equations of unsteady pipe flow.

For unsteady flow in a closed pipe (area $A_p$, diameter $D$, wave speed $a$, Darcy-Weisbach factor $f$), with head $H$ and discharge $Q$:

Momentum:

$$L_1=\frac{\partial Q}{\partial t}+gA_p\frac{\partial H}{\partial x}+\frac{f\,Q|Q|}{2DA_p}=0$$

Continuity:

$$L_2=\frac{\partial H}{\partial t}+\frac{a^2}{gA_p}\frac{\partial Q}{\partial x}=0$$

Along $C^+$ ($dx/dt=+a$) and $C^-$ ($dx/dt=-a$):

$$C^+:\ \frac{dH}{dt}+\frac{a}{gA_p}\frac{dQ}{dt}+\frac{af}{2gDA_p^2}Q|Q|=0,\qquad C^-:\ \frac{dH}{dt}-\frac{a}{gA_p}\frac{dQ}{dt}-\frac{af}{2gDA_p^2}Q|Q|=0$$

Let $B_p=\dfrac{a}{gA_p}$.

(a) Rectangular grid: characteristics through the diagonal

Divide the pipe into $N$ reaches of length $\Delta x=L/N$ and choose $\Delta t=\Delta x/a$. The $C^+$ and $C^-$ lines through P then pass exactly through the grid points A and B of the previous time level.

 t+dt :         P
               /|\
         C+  /  |  \  C-
 t    :    A----C----B

Integrating $C^+$ from A to P and $C^-$ from B to P (friction at the known end), with $R_p=\dfrac{f\Delta x}{2gDA_p^2}$:

$$H_P=H_A-B_p(Q_P-Q_A)-R_pQ_A|Q_A|,\qquad H_P=H_B+B_p(Q_P-Q_B)+R_pQ_B|Q_B|$$

With $K_1=H_A+B_pQ_A-R_pQ_A|Q_A|$ and $K_2=H_B-B_pQ_B+R_pQ_B|Q_B|$:

$$Q_P=\frac{K_1-K_2}{2B_p},\qquad H_P=\frac{K_1+K_2}{2}$$

No interpolation is needed. Boundaries: upstream reservoir, $H_P=H_{res}$ and $Q_P=(H_{res}-K_2)/B_p$; closed downstream valve, $Q_P=0$ and $H_P=K_1$.

(b) Specified time intervals: characteristics not on the diagonal

When $\Delta t$ and $\Delta x$ are chosen independently (for example to get the same $\Delta t$ in pipes of different length), the Courant condition $a\Delta t\le\Delta x$ must hold. The characteristics through P cut the time line at points R and S between grid points.

 t+dt :              P
                    /|\
              C+  /  |  \  C-
 t    :   A-----R----C----S-----B

With $\theta=\dfrac{a\Delta t}{\Delta x}\le1$, linear interpolation gives

$$Q_R=Q_C-\theta(Q_C-Q_A),\quad H_R=H_C-\theta(H_C-H_A)$$ $$Q_S=Q_C-\theta(Q_C-Q_B),\quad H_S=H_C-\theta(H_C-H_B)$$

The characteristic equations from R and S to P (friction length $a\Delta t$, so $R_p'=\dfrac{f\,a\Delta t}{2gDA_p^2}$) are

$$H_P=H_R-B_p(Q_P-Q_R)-R_p'Q_R|Q_R|,\qquad H_P=H_S+B_p(Q_P-Q_S)+R_p'Q_S|Q_S|$$ $$Q_P=\frac{K_1'-K_2'}{2B_p},\quad H_P=\frac{K_1'+K_2'}{2},\qquad K_1'=H_R+B_pQ_R-R_p'Q_R|Q_R|,\ \ K_2'=H_S-B_pQ_S+R_p'Q_S|Q_S|$$

The interpolation introduces numerical damping (diffusion), which grows as $\theta$ departs from 1, so the diagonal grid of part (a) is preferred when possible.