Method of Characteristics

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.

Governing equations (velocity $V$, pressure $p$)

For a pipe of diameter $D$ inclined at angle $\theta$ to the horizontal (water density $\rho$, wave speed $a$, friction factor $f$), neglecting the small convective terms:

Momentum:

$$L_1=\frac{\partial V}{\partial t}+\frac1\rho\frac{\partial p}{\partial x}+g\sin\theta+\frac{f}{2D}V|V|=0$$

Continuity:

$$L_2=\frac{\partial p}{\partial t}+\rho a^2\frac{\partial V}{\partial x}=0$$

Characteristic equations

Form $L_1+\lambda L_2=0$:

$$\left[\frac{\partial V}{\partial t}+\lambda\rho a^2\frac{\partial V}{\partial x}\right]+\lambda\left[\frac{\partial p}{\partial t}+\frac{1}{\lambda\rho}\frac{\partial p}{\partial x}\right]+g\sin\theta+\frac{f}{2D}V|V|=0$$

For both brackets to be total derivatives, $\dfrac{dx}{dt}=\lambda\rho a^2=\dfrac{1}{\lambda\rho}$, so $\lambda=\pm\dfrac{1}{\rho a}$ and $\dfrac{dx}{dt}=\pm a$:

$$C^+\ \left(\frac{dx}{dt}=+a\right):\quad \frac{dV}{dt}+\frac{1}{\rho a}\frac{dp}{dt}+g\sin\theta+\frac{f}{2D}V|V|=0$$ $$C^-\ \left(\frac{dx}{dt}=-a\right):\quad \frac{dV}{dt}-\frac{1}{\rho a}\frac{dp}{dt}+g\sin\theta+\frac{f}{2D}V|V|=0$$

Finite difference form

Use the rectangular grid with A, C, B at time $t$ and P at $t+\Delta t$, with $\Delta x=a\Delta t$. Multiply by $\rho a$, integrate along each characteristic (friction at the known point), and put $\Delta x\sin\theta=\Delta z$ (rise of the pipe over $\Delta x$):

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

$$C^+:\quad p_P=p_A-\rho a\left(V_P-V_A\right)-\rho g\left(z_P-z_A\right)-\frac{\rho f\Delta x}{2D}V_A|V_A|$$ $$C^-:\quad p_P=p_B+\rho a\left(V_P-V_B\right)+\rho g\left(z_B-z_P\right)+\frac{\rho f\Delta x}{2D}V_B|V_B|$$

Solution

Equating the two expressions for $p_P$:

$$V_P=\frac{(p_A-p_B)+\rho a(V_A+V_B)-\rho g(z_B-z_A)-\dfrac{\rho f\Delta x}{2D}\left[V_A|V_A|+V_B|V_B|\right]}{2\rho a}$$

and $p_P$ is then obtained from either equation. At each time step, the interior points are found this way; boundary points use one characteristic plus the boundary condition (e.g. $p=$ reservoir pressure, or $V=0$ at a closed valve).

Idea

When the characteristics do not lie on the diagonal of the grid, $a\Delta t<\Delta x$: the $C^+$ and $C^-$ lines through P (at $t_2$) meet the line $t_1$ at points R and S that lie between the grid points, not at A and B. $Q$ and $H$ at R and S are found by linear interpolation along the time line $t_1$.

 t2 :              P
                  /|\
            C+  /  |  \  C-
 t1 :   A-----R----C----S-----B
        (x-dx)         (x+dx)
        |<-a dt->|   |<-a dt->|

Here R lies on the $C^+$ line between A and C, and S on the $C^-$ line between C and B, each at a distance $a\Delta t$ from C.

Interpolation factor

$$\theta=\frac{a\,\Delta t}{\Delta x}=\frac{850\times1}{1000}=0.85$$

Point R (between A and C, $C^+$ foot):

$$Q_R=Q_C-\theta\,(Q_C-Q_A),\qquad H_R=H_C-\theta\,(H_C-H_A)$$ $$Q_R=0.64-0.85(0.64-0.6)=0.606\ \text{m}^3/\text{s},\qquad H_R=20.4-0.85(20.4-20)=20.06\ \text{m}$$

Point S (between C and B, $C^-$ foot):

$$Q_S=Q_C-\theta\,(Q_C-Q_B),\qquad H_S=H_C-\theta\,(H_C-H_B)$$ $$Q_S=0.64-0.85(0.64-0.65)=0.6485\ \text{m}^3/\text{s},\qquad H_S=20.4-0.85(20.4-20.6)=20.57\ \text{m}$$

Answer: at $t_1=1$ s, $Q_R=0.606$ m$^3$/s, $H_R=20.06$ m on the $C^+$ line, and $Q_S=0.6485$ m$^3$/s, $H_S=20.57$ m on the $C^-$ line. These replace $Q_A,H_A$ and $Q_B,H_B$ in the $C^+$ and $C^-$ equations

$$H_P=H_R-B_p(Q_P-Q_R)-R_pQ_R|Q_R|,\qquad H_P=H_S+B_p(Q_P-Q_S)+R_pQ_S|Q_S|$$

with $B_p=a/(gA_p)$ and $R_p=f\,a\Delta t/(2gDA_p^2)$ (the friction length is $a\Delta t=850$ m).

Idea

When the characteristics do not lie on the diagonal of the grid, $a\Delta t<\Delta x$: the $C^+$ and $C^-$ lines through P (at $t_2$) meet the line $t_1$ at points R and S that lie between the grid points, not at A and B. $Q$ and $H$ at R and S are found by linear interpolation along the time line $t_1$.

 t2 :              P
                  /|\
            C+  /  |  \  C-
 t1 :   A-----R----C----S-----B
        (x-dx)         (x+dx)
        |<-a dt->|   |<-a dt->|

Here R lies on the $C^+$ line between A and C, and S on the $C^-$ line between C and B, each at a distance $a\Delta t$ from C.

Interpolation factor

$$\theta=\frac{a\,\Delta t}{\Delta x}=\frac{800\times1}{1000}=0.80$$

Point R (between A and C, $C^+$ foot):

$$Q_R=Q_C-\theta\,(Q_C-Q_A),\qquad H_R=H_C-\theta\,(H_C-H_A)$$ $$Q_R=0.74-0.80(0.74-0.7)=0.708\ \text{m}^3/\text{s},\qquad H_R=20.4-0.80(20.4-20)=20.08\ \text{m}$$

Point S (between C and B, $C^-$ foot):

$$Q_S=Q_C-\theta\,(Q_C-Q_B),\qquad H_S=H_C-\theta\,(H_C-H_B)$$ $$Q_S=0.74-0.80(0.74-0.76)=0.756\ \text{m}^3/\text{s},\qquad H_S=20.4-0.80(20.4-20.6)=20.56\ \text{m}$$

Answer: at $t_1=1$ s, $Q_R=0.708$ m$^3$/s, $H_R=20.08$ m on the $C^+$ line, and $Q_S=0.756$ m$^3$/s, $H_S=20.56$ m on the $C^-$ line. These replace $Q_A,H_A$ and $Q_B,H_B$ in the $C^+$ and $C^-$ equations

$$H_P=H_R-B_p(Q_P-Q_R)-R_pQ_R|Q_R|,\qquad H_P=H_S+B_p(Q_P-Q_S)+R_pQ_S|Q_S|$$

with $B_p=a/(gA_p)$ and $R_p=f\,a\Delta t/(2gDA_p^2)$ (the friction length is $a\Delta t=800$ m).

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.

Use the finite difference form of the characteristic equations on the rectangular grid (A and B at time $t$, P at $t+\Delta t$):

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

$H$ is the piezometric head (it already contains the elevation), so the 1 m elevation difference between A and P is not needed to find $H_P$ and $Q_P$.

Step 1: Constants

$$A_p=\frac{\pi}{4}(0.3)^2=0.07069\ \text{m}^2$$ $$B_p=\frac{1200}{9.81\times0.07069}=1730.53\ \text{s/m}^2,\qquad R_p=\frac{0.02\times500}{2\times9.81\times0.3\times0.07069^2}=340.03\ \text{s}^2/\text{m}^5$$

(Here $a\Delta t=1200\times0.4=480$ m is close to $\Delta x=500$ m, so the characteristics are taken through A and B as in the figure.)

Step 2: Known quantities

$$K_1=H_A+B_pQ_A-R_pQ_A|Q_A|=26.5+1730.53(0.4)-340.03(0.16)=664.31$$ $$K_2=H_B-B_pQ_B+R_pQ_B|Q_B|=27.5-1730.53(0.45)+340.03(0.2025)=-682.38$$

Step 3: Solve for P

$$Q_P=\frac{K_1-K_2}{2B_p}=\frac{664.31-(-682.38)}{2\times1730.53}=0.3891\ \text{m}^3/\text{s}$$ $$H_P=\frac{K_1+K_2}{2}=-9.038\ \text{m}$$

Answer: $Q_P\approx0.3891$ m$^3$/s and $H_P\approx-9.038$ m. The negative head means the data give a pressure below the vapour pressure at P (column separation would occur), because the large impedance $B_p$ makes the $Q$ difference between A and B produce a big head change.

Method of characteristics

The method of characteristics replaces the two partial differential equations of unsteady flow by ordinary differential equations along the characteristic lines $dx/dt=\pm a$ in the $x$-$t$ plane (the $C^+$ and $C^-$ lines). These are integrated by finite differences from the known time level to the next. It is widely used for water hammer.

Numerical diffusion

Smoothing of the solution (loss of amplitude) caused by the truncation error of the scheme. In the MOC it arises when the characteristic foot falls between grid points and $Q,H$ are found by linear interpolation: sharp pressure peaks are averaged and flattened.

Numerical dispersion

Spreading of waves caused by phase errors: different frequency components travel at different speeds, so a steep wave front becomes distorted and may show spurious oscillations. It is a property of the discretisation, not of the physical pipe.

Stability

A scheme is stable if errors made at one step do not grow in later steps. For the MOC with interpolation, the Courant condition must hold:

$$C=\frac{a\,\Delta t}{\Delta x}\le1$$

so that the characteristic foot lies within the grid interval used for interpolation. $C=1$ (characteristics along the grid diagonal) gives no interpolation error, no numerical diffusion and the most accurate result.

Set-up

One reach means two nodes: A at the reservoir ($x=0$) and B at the dead end ($x=L$). The time step follows from $\Delta x=L$ and $\Delta t=L/a$:

$$\Delta t=\frac{500}{1250}=0.40\ \text{s}$$ $$A_p=\frac{\pi}{4}(0.45)^2=0.15904\ \text{m}^2,\quad B_p=\frac{a}{gA_p}=801.17\ \text{s/m}^2,\quad R_p=\frac{f\,L}{2gDA_p^2}=44.78\ \text{s}^2/\text{m}^5$$

Equations

$C^+$ (from A at the old time to B) and $C^-$ (from B at the old time to A):

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

• Node A (reservoir): $H_{PA}=100+2\sin(\pi t)$ is known; use $C^-$ from B:

$$Q_{PA}=\frac{H_{PA}-K_2}{B_p},\qquad K_2=H_B-B_pQ_B+R_pQ_B|Q_B|$$

• Node B (closed end): $Q_{PB}=0$; use $C^+$ from A:

$$H_{PB}=K_1=H_A+B_pQ_A-R_pQ_A|Q_A|$$

Initial condition ($t=0$)

No flow and static head: $Q_A=Q_B=0$, $H_A=H_B=100$ m.

First step ($t=0.40$ s)

$H_{PA}=100+2\sin(0.40\pi)=101.9021$ m; $K_2=100-0+0=100$; $Q_{PA}=(101.9021-100)/801.17=0.00237$ m$^3$/s; $K_1=100$, so $H_{PB}=100$ m.

Marching in time

Repeating the two boundary calculations at every step (each uses the values of the previous step):

($Q_B=0$ at all times.)

Answer: at $t=2$ s, the discharge from A is $Q_A=0.00768$ m$^3$/s and the HGL elevation at B is $H_B=93.846$ m.

Characteristic curve

In the $x$-$t$ plane, a characteristic curve is a line $dx/dt=\pm a$ along which the partial differential equations of unsteady pipe flow reduce to ordinary differential equations. It is the path of a small pressure disturbance: the $C^+$ curve ($dx/dt=+a$) carries information downstream and the $C^-$ curve ($dx/dt=-a$) upstream. The solution at a point $P$ is found from the values at the two points A and B where the $C^+$ and $C^-$ curves through $P$ meet the previous time line.

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

Numerical part

Equations used (velocity and pressure form, $\rho=1000$ kg/m$^3$, $\Delta x=500$ m so $a\Delta t=1000\times0.5=500$ m):

$$C^+:\ p_P=p_A-\rho a(V_P-V_A)-\rho g(z_P-z_A)-\frac{\rho f\Delta x}{2D}V_A|V_A|$$ $$C^-:\ p_P=p_B+\rho a(V_P-V_B)+\rho g(z_B-z_P)+\frac{\rho f\Delta x}{2D}V_B|V_B|$$

Assumption: the pipe has a uniform slope, so P is 2.50 m above A and 2.50 m below B ($z_B-z_A=5.0$ m, $z_P-z_A=2.5$ m, $z_B-z_P=2.5$ m).

Constants

• $\rho a=1000\times1000=1.0\times10^{6}$ kg/(m$^2$s)
• $\dfrac{\rho f\Delta x}{2D}=\dfrac{1000\times0.02\times500}{2\times0.35}=14285.71$ kg/m$^3$ (friction coefficient)
• $\rho g(z_B-z_A)=9810\times5.0=49050$ N/m$^2$

Velocity at P

$$V_P=\frac{(p_A-p_B)+\rho a(V_A+V_B)-\rho g(z_B-z_A)-\frac{\rho f\Delta x}{2D}\left[V_A|V_A|+V_B|V_B|\right]}{2\rho a}$$ $$V_P=\frac{(102000-124000)+10^6(12.25)-49050-14285.71(36+39.0625)}{2\times10^6} =\frac{11106629}{2\times10^6}=5.5533\ \text{m/s}$$

Pressure at P (from $C^+$)

$$p_P=102000-10^6(5.5533-6)-9810(2.5)-14285.71(36)=9875\ \text{N/m}^2=9.88\ \text{kN/m}^2$$

Check with $C^-$: $p_P=124000+10^6(5.5533-6.25)+9810(2.5)+14285.71(39.0625)=9875$ N/m$^2$ (same).

(The same result follows from the head form: $B_p=a/g=101.94$ s, $H_A=10.398$ m, $H_B=17.640$ m, giving $H_P=3.507$ m and $p_P=\gamma(H_P-z_P)$.)

Answer: $V_P\approx5.5533$ m/s and $p_P\approx9.88$ kN/m$^2$.

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

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. Equivalent: $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$, number of reaches $N$, total time $T$, 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 the boundary condition (e.g. 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|$ together with the boundary condition (closed end $Q_P=0$, so $H_P=H_N+BQ_N-RQ_N|Q_N|$; or 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 time.
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 required.