Method of characteristics; FD form for pipe flow

Method of characteristics (MOC)

The method of characteristics converts the two partial differential equations of unsteady flow (a hyperbolic system) into ordinary differential equations that hold only along special lines in the $x$-$t$ plane, called characteristic lines. These ordinary equations are then integrated by finite differences.

Why it is necessary:

• The PDEs (continuity and momentum) are non-linear because of the friction term $Q|Q|$ and cannot be solved analytically for real boundary conditions.
• Disturbances (pressure waves, water hammer) travel at finite speed $a$; the characteristics follow exactly these waves, so the numerical solution keeps the physical wave behaviour and sharp fronts are captured.
• It handles boundary conditions (reservoir, valve, pump) conveniently, and gives accurate results for water hammer.

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$$

Conversion to characteristic (ordinary differential) equations

Combine the two equations linearly with an unknown multiplier $\lambda$: $L_1+\lambda L_2=0$.

$$\left[\frac{\partial Q}{\partial t}+\lambda\frac{a^2}{gA_p}\frac{\partial Q}{\partial x}\right]+\lambda\left[\frac{\partial H}{\partial t}+\frac{gA_p}{\lambda}\frac{\partial H}{\partial x}\right]+\frac{fQ|Q|}{2DA_p}=0$$

If the two brackets are total derivatives $\dfrac{dQ}{dt}$ and $\dfrac{dH}{dt}$ along a line $dx/dt$, then

$$\frac{dx}{dt}=\lambda\frac{a^2}{gA_p}=\frac{gA_p}{\lambda}\ \Rightarrow\ \lambda=\pm\frac{gA_p}{a},\qquad \frac{dx}{dt}=\pm a$$

This gives two pairs of equations:

Multiplying by $a/(gA_p)$:

$$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$$

Finite difference form (rectangular grid, $\Delta x=a\Delta t$)

 t
 t+dt :          P            unknown H_P, Q_P
                / \
          C+  /   \  C-
 t    :     A --C-- B         known
         (x-dx)  (x)  (x+dx)        --> x

Integrate $C^+$ from A to P and $C^-$ from B to P, taking the friction term at the known point (first-order) and using $a\Delta t=\Delta x$:

$$H_P=H_A-B_p\left(Q_P-Q_A\right)-R_pQ_A|Q_A|\qquad(C^+)$$ $$H_P=H_B+B_p\left(Q_P-Q_B\right)+R_pQ_B|Q_B|\qquad(C^-)$$ $$B_p=\frac{a}{gA_p},\qquad R_p=\frac{f\,\Delta x}{2gDA_p^2}$$

Define the known quantities

$$K_1=H_A+B_pQ_A-R_pQ_A|Q_A|,\qquad K_2=H_B-B_pQ_B+R_pQ_B|Q_B|$$

Then $H_P=K_1-B_pQ_P$ and $H_P=K_2+B_pQ_P$, so

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

With $H_P$ and $Q_P$ known for all interior points, the boundary points are found from one characteristic plus the boundary condition (e.g. a reservoir: $H_P=H_{res}$; closed valve: $Q_P=0$), and the procedure is repeated for the next time step.