---
title: "Algorithm for unsteady pipe flow FD equations"
subject: "Computational Techniques in Civil Engineering"
chapter: 6
key: "semester_8/computational_techniques_in_civil_engineering:algorithm-pipe-fd-hp-equations"
source: https://www.bctengineer.com/q/semester_8/computational_techniques_in_civil_engineering/algorithm-pipe-fd-hp-equations
site: BCT Engineering
---

# Algorithm for unsteady pipe flow FD equations

Computational Techniques in Civil Engineering · Chapter 6: Method of Characteristics

- Asked in: 2070 Bhadra, 8 marks
- Answer checked by Claude Opus 5.5 (high effort) on 2026-10-14
- Page: https://www.bctengineer.com/q/semester_8/computational_techniques_in_civil_engineering/algorithm-pipe-fd-hp-equations
- Whole chapter: https://www.bctengineer.com/md/semester_8/computational_techniques_in_civil_engineering/chapter-6

## Question

Prepare an algorithm to compute discharge and head based on the following form of finite difference equations for unsteady pipe flow problem using a rectangular grid.

$H_{P_i} = H_{i-1} - B(Q_{P_i} - Q_{i-1}) - R Q_{i-1}|Q_{i-1}|$

$H_{P_i} = H_{i+1} + B(Q_{P_i} - Q_{i+1}) + R Q_{i+1}|Q_{i+1}|$

where $H$ = head, $Q$ = discharge, $H_{P_i}$ and $Q_{P_i}$ = head and discharge at the point of intersection of two characteristics, $B$ and $R$ = coefficients.

## Short answer

Set up grid with Δt=Δx/a. At interior nodes, equate the C+ and C− equations to solve for Q_{P_i}=[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), then compute H_{P_i} from either. Upstream boundary uses C− with known head; downstream uses C+ with given condition. Update all nodes and march in time.

## Answer

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

```text
 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.

## Check yourself

The answer's key points as multiple-choice questions, with the correct option and why.

1. How is Q_Pi found at an interior node from the C+ and C- equations?
   A. Equate the two H_Pi expressions and solve for Q_Pi, giving a numerator over 2B
   B. Add them and divide by B
   C. Use C+ only and assume H_Pi = H_(i-1)
   D. Take the average of Q_(i-1) and Q_(i+1)
   Answer: A. Equating the expressions gives Q_Pi = [(H_(i-1) - H_(i+1)) + B(Q_(i-1)+Q_(i+1)) - R(...)]/(2B). (Key point: Solve for Q_P at interior node)

2. Which relation must hold between dt and dx for the rectangular grid algorithm?
   A. dt = a dx
   B. dt = dx/a (Courant number 1) so no interpolation is required
   C. dt = 2 dx/a
   D. dt can be any value
   Answer: B. dt = dx/a puts the characteristic feet on the grid nodes of the rectangular grid. (Key point: Time step relation)

3. At the upstream boundary (i = 1), which equation and extra information is used?
   A. C+ equation with the boundary condition
   B. Both C+ and C- without boundary data
   C. Only C- is available, so a boundary condition such as reservoir H_P1 = H_res is needed
   D. The interior formula with i - 1 = 0
   Answer: C. At the end nodes only one characteristic arrives, so a boundary condition supplies the second equation. (Key point: Boundary nodes)

4. What must be done after all nodes are computed in a time step?
   A. Increase B and R
   B. Reset t to zero
   C. Compute Q only
   D. Set Q_i = Q_Pi and H_i = H_Pi for all i, then continue until t >= T
   Answer: D. The new values become the known values for the next step. (Key point: Update step)

## Related in Computational Techniques in Civil Engineering

- [Method of characteristics; FD form for pipe flow](https://www.bctengineer.com/question.md?key=semester_8%2Fcomputational_techniques_in_civil_engineering%3Amoc-definition-fd-pipe-equations) — Computational Techniques in Civil Engineering, chapter 6, asked 4×
- [MOC numerical: discharge and head at P (30 cm pipe)](https://www.bctengineer.com/question.md?key=semester_8%2Fcomputational_techniques_in_civil_engineering%3Amoc-pipe-discharge-head-pc-30cm) — Computational Techniques in Civil Engineering, chapter 6
- [Characteristic form FD solution (velocity, pressure)](https://www.bctengineer.com/question.md?key=semester_8%2Fcomputational_techniques_in_civil_engineering%3Acharacteristic-fd-velocity-pressure) — Computational Techniques in Civil Engineering, chapter 6, asked 2×
- [Water hammer: rectangular grid and specified intervals](https://www.bctengineer.com/question.md?key=semester_8%2Fcomputational_techniques_in_civil_engineering%3Awater-hammer-numerical-solution) — Computational Techniques in Civil Engineering, chapter 6

## Related in other subjects

- [Resection, two point problem, intersection](https://www.bctengineer.com/question.md?key=semester_4%2Fsurveying_ii%3Aresection-two-point-problem-ab-p) — Surveying II, chapter 5
- [Velocity, head loss, discharge in pipes BC, CD, BF, FD (0.1 m)](https://www.bctengineer.com/question.md?key=semester_5%2Fwater_supply_engineering%3Aloop-c72) — Water Supply Engineering, chapter 7
- [Velocity, head loss, discharge in loop pipes BC, CD, BF, FD](https://www.bctengineer.com/question.md?key=semester_5%2Fwater_supply_engineering%3Aloop-a76) — Water Supply Engineering, chapter 7
- [Derive continuity equation for unsteady pipe flow](https://www.bctengineer.com/question.md?key=semester_4%2Fhydraulics%3Acontinuity-equation-unsteady-pipe-flow) — Hydraulics, chapter 4