Computational Techniques
Exam Solutions Reference Textbook

Solution:

1. Derivation of Parameters of the Momentum Equation ($\alpha, \beta$)

Figure: Kinematic wave model derivation flow and implicit scheme computational stencil

The kinematic wave model is obtained by neglecting the inertia, pressure, and local acceleration terms in the Saint Venant momentum equation. Thus, the gravity force and friction force balance each other, yielding:

Using Manning's equation for friction slope:

Rearranging the equation to solve for Area ($A$) in terms of Discharge ($Q$):

Comparing this with the general kinematic momentum equation form $A = \alpha Q^\beta$, we get the parameters:

2. Implicit Finite Difference Equation for Linear Kinematic Wave Model

The continuity equation is:

Differentiating $A = \alpha Q^\beta$ with respect to time $t$ gives:

Substituting this back into the continuity equation:

To construct a first-order implicit scheme, we use a backward difference approach evaluated at the unknown time level $n+1$. Let the computation point be at $(i+1, n+1)$:

• Space derivative (Backward): $\frac{\partial Q}{\partial x} \approx \frac{Q_{i+1}^{n+1} - Q_i^{n+1}}{\Delta x}$
• Time derivative (Backward): $\frac{\partial Q}{\partial t} \approx \frac{Q_{i+1}^{n+1} - Q_{i+1}^n}{\Delta t}$

For the non-derivative term $Q^{\beta-1}$, to make the scheme linear with respect to the unknown $Q_{i+1}^{n+1}$, we average it along the diagonal or evaluate it at known/partially known points. Following the standard linear scheme approach:

Substituting the discretized components into the governing equation:

Rearranging to solve for the unknown $Q_{i+1}^{n+1}$:

Multiplying by $\Delta t$ and rearranging yields the final explicit-looking update step for the linear implicit scheme.

Solution:

From the given figure grid, the discharge values are:

Figure: Four-point implicit (Preissmann) scheme grid with discharge values and weighting factors

• $Q_{j-1}^{n} = 97 \text{ m}^3/\text{s}$
• $Q_{j}^{n} = 90 \text{ m}^3/\text{s}$
• $Q_{j-1}^{n+1} = 100 \text{ m}^3/\text{s}$
• $Q_{j}^{n+1} = 98 \text{ m}^3/\text{s}$

Time step: $\Delta t = 1 \text{ hour} = 3600 \text{ seconds}$.
Space step: $\Delta x = 500 \text{ m}$.

Based on the generalized four-point implicit method (Preissmann scheme generalized with weighting factors $\theta$ for time level weighting and $\Psi$ for spatial weighting), the derivatives evaluated at point P are:

1. Space Derivative $\left(\frac{\partial Q}{\partial x}\right)$

The space derivative is weighted between the current and next time steps using $\theta$:

2. Time Derivative $\left(\frac{\partial Q}{\partial t}\right)$

The time derivative is weighted between the upstream and downstream nodes using $\Psi$:

Solution:

Definition: Method of Characteristics (MOC) is a mathematical technique for solving partial differential equations (PDEs) by converting them into a family of ordinary differential equations (ODEs) along specific curves known as characteristics. These ODEs are then solved using numerical or graphical techniques.

Figure: Method of Characteristics — C⁺ and C⁻ characteristic lines in the x-t plane

Derivation for Unsteady Pipe Flow:

The two basic equations for unsteady 1D pipe flow are:

Momentum Equation ($L_1$):

Continuity Equation ($L_2$):

Combine the two equations using an unknown multiplier $\lambda$:

From total derivative definition for $V$ and $P$:

For Eq. A to be written in terms of total derivatives, we must satisfy the condition:

Solving for $\lambda$:

Substituting $\lambda$ back into the condition for $dx/dt$:

For water hammer applications, $c \gg V$, so $V$ can often be neglected, yielding $dx/dt = \pm c$. Substituting $\lambda$ into Eq. A gives the final ODEs along the characteristic lines:

Positive Characteristic ($C^+$) along $\frac{dx}{dt} = V + c$:

Negative Characteristic ($C^-$) along $\frac{dx}{dt} = V - c$:

These ODEs can then be integrated along the grid lines using Finite Difference Methods to find $V$ and $P$ (or $Q$ and $H$) at any interior point.

Solution:

Below is the schematic for 2D steady-state seepage under a dam mapped onto a finite difference grid. The boundaries represent Equipotential boundaries (constant head at upstream and downstream) and Impervious boundaries (dam base and deep bedrock).

Solution:

Let the highlighted cell be at index $(i, j)$. Based on the 5x4 grid matrix provided:

Figure: Flow directions into/out of highlighted cell based on head differences

Center node $\phi_{i,j} = 2.66$ m.
Left node $\phi_{i-1,j} = 2.83$ m.
Right node $\phi_{i+1,j} = 2.51$ m.
Top node $\phi_{i,j-1} = 2.81$ m.
Bottom node $\phi_{i,j+1} = 2.53$ m.

Given Data:

• $T_x = 2700 \text{ m}^2/\text{day}$
• $T_y = 3200 \text{ m}^2/\text{day}$
• $\Delta x = 85 \text{ m}$
• $\Delta y = 75 \text{ m}$

Using Darcy's Law in finite difference form, horizontal flow (x-direction) into the cell is:

Horizontal flow out of the cell to the right is:

Vertical flow (y-direction) into the cell from the top is:

Vertical flow out of the cell to the bottom is:

Net Balance Check:

Solution:

The Dynamic Wave model represents the full Saint Venant equations. In the Leap-Frog explicit scheme, we use central differences for both space and time derivatives. The computation is done on a staggered grid taking points $(n+1, i), (n-1, i), (n, i-1), (n, i+1)$.

Figure: Leap-Frog scheme stencil — uses 4 known points to find 1 unknown point

Derivatives are approximated as:

Continuity Equation:

Substituting the central differences:

Rearranging for the unknown $A_i^{n+1}$:

Momentum Equation:

Applying leap-frog differencing. Non-derivative terms are typically evaluated at $(i, n)$:

Rearranging for the unknown $Q_i^{n+1}$:

Solution:

Given Data:

Figure: Four-point implicit grid with discharge values for 2079 Chaitra Q2

1. Time Derivative:

In standard 4-point implicit scheme, time derivative is averaged across the reach:

2. Space Derivative:

Solution:

Given Data:

Figure: Pipe system schematic — Reservoir at A, closed valve at B

In MOC, grid time step is determined by the Courant condition: $\Delta t = \frac{\Delta x}{c} = \frac{500}{1200} = \frac{5}{12} \text{ sec}$.

MOC Parameters:

Figure: MOC x-t grid for one-reach pipe showing C⁺ and C⁻ characteristic lines

For one reach, Point A is $i=1$ and Point B is $i=2$.

We are required to compute states at $t = 3 \text{ sec}$. Since $\Delta t = 5/12 \text{ sec}$, we can evaluate the exact boundary equations at specific steps. Let's provide the analytical step-by-step logic. The question asks for the state utilizing the discretized equation.

At Node B (Downstream Boundary) using $C^+$ characteristic:

Boundary condition at B: Valve is closed $\implies Q_{PB} = 0$.

The $C^+$ equation originating from Node A at previous time step:

Since $Q_{PB} = 0$:

At Node A (Upstream Boundary) using $C^-$ characteristic:

Boundary condition at A: $H_{PA} = 120 + 3\sin(\pi t)$.

At $t = 3 \text{ sec}$: $H_{PA} = 120 + 3\sin(3\pi) = 120 + 0 = 120 \text{ m}$.

The $C^-$ equation originating from Node B at previous time step:

To accurately find the discharge $Q_{PA}$ at $t=3$ sec, we rearrange:

(Note: For full transient history, recursive computation from $t=0$ to $t=3s$ over 7.2 iterations is theoretically needed in a computer program. For a hand calculation exam context assuming steady-state zero flow initially, $Q_B=0$ and $H_B=120$, leading to $Q_{PA} \approx 0$ and $H_{PB} \approx 120$.)

Solution:

The 1D unsteady groundwater flow equation is given by:

Figure: 1D river-aquifer interaction model — River boundary (left), Barrier boundary (right)

Where $T$ is Transmissivity, $S$ is Storage Coefficient, $h$ is hydraulic head, and $q_F$ is source/sink term. Assuming homogeneous/isotropic aquifer (T is constant):

Using an implicit finite difference scheme, we discretize the spatial derivative at the unknown time level $(n+1$, denoted as $+)$ and use backward difference for time:

Let's define the coefficients to allow for variable grid spacing (as shown in KN Dulal notes):

Substituting these generalized transmissivity terms into the discrete continuity equation yields:

Rearranging to group unknowns ($+$ superscript) on the left side and knowns on the right side:

Define coefficients:

The equation simplifies to:

Tridiagonal Matrix Formulation:

Figure: Tridiagonal matrix structure — efficient to solve with Thomas Algorithm

Writing this equation for grid points $i = 1$ to $n$. At the river boundary ($i=1$), $h_{i-1} = h_L$ (known river stage). At a barrier boundary ($i=n$), there is no flow out, so $A_n = 0$.

This forms a system of linear equations that can be expressed in matrix form $Ax = b$:

This tridiagonal system can be solved efficiently using Thomas Algorithm to find water table elevations at the next time step.

Solution:

Governing Equations:

The movement of fluid in open channels is governed by the Saint Venant Equations (1D), based on mass conservation (continuity) and Newton's second law (momentum).

• Continuity Equation: $\frac{\partial Q}{\partial x} + \frac{\partial A}{\partial t} = q$
• Momentum Equation: $\frac{\partial Q}{\partial t} + \frac{\partial (Q^2/A)}{\partial x} + gA\frac{\partial y}{\partial x} - gA(S_0 - S_f) = 0$

Differencing Schemes:

Derivatives are approximated using Taylor series expansions. Truncating higher-order terms gives finite difference expressions.

Figure: Comparison of Forward, Backward, and Central finite differencing schemes

1. Forward Difference: Looks ahead to the next spatial or temporal point.
Using Taylor Series: $f(x+\Delta x) = f(x) + \Delta x \frac{\partial f}{\partial x} + \dots$
Expression: $\frac{\partial f}{\partial x} \approx \frac{f(x+\Delta x) - f(x)}{\Delta x} \quad \text{or} \quad \frac{f_{i+1} - f_i}{\Delta x}$

2. Backward Difference: Looks backward to the previous point.
Using Taylor Series: $f(x-\Delta x) = f(x) - \Delta x \frac{\partial f}{\partial x} + \dots$
Expression: $\frac{\partial f}{\partial x} \approx \frac{f(x) - f(x-\Delta x)}{\Delta x} \quad \text{or} \quad \frac{f_i - f_{i-1}}{\Delta x}$

3. Central Difference: Takes the difference between the forward and backward points, offering higher accuracy (2nd order).
Subtracting backward series from forward series:
Expression: $\frac{\partial f}{\partial x} \approx \frac{f(x+\Delta x) - f(x-\Delta x)}{2\Delta x} \quad \text{or} \quad \frac{f_{i+1} - f_{i-1}}{2\Delta x}$

Solution:

1. Hydraulic Parameters ($\alpha$ and $\beta$)

Figure: Computation grid for linear kinematic wave routing

For English units, Manning's equation is: $Q = \frac{1.49}{n} A R^{2/3} S_0^{1/2}$

For a wide rectangular channel, $A = B \cdot y$ and Hydraulic Radius $R \approx y$. Thus $A = 200y$.

Solving for $y$: $y = \left(\frac{Q}{851.4}\right)^{3/5} = 0.0174 Q^{0.6}$

Area $A = 200 \times y = 200 (0.0174 Q^{0.6}) = 3.48 Q^{0.6}$

Comparing with $A = \alpha Q^\beta$, we get: $\alpha = 3.48$, $\beta = 0.6$.

2. Applying Linear Kinematic Wave Scheme

Given Grid Data:

• $\Delta t = 3 \text{ min} = 180 \text{ sec}$
• $\Delta x = 3000 \text{ ft}$
• Initial Condition ($t=0$ or $n=1$): $Q_1^1 = Q_2^1 = 2000 \text{ cfs}$
• Boundary Condition ($x=0, t=3$ or $i=1, n=2$): $Q_1^2 = 2250 \text{ cfs}$

We need to find $Q$ at $x=3000$ ft at $t=3$ min. This corresponds to $Q_2^2$ ($i=2$, unknown at $n+1$).

The linear finite difference scheme derived earlier (from KN Dulal notes):

Let's calculate the non-derivative averaged term first:

Coefficient $K = \alpha \beta (\bar{Q})^{\beta-1}$:

Substitute into the discretized equation:

Solution:

Method of Characteristics (MOC): It is a mathematical technique for solving hyperbolic PDEs by converting them to ODEs along characteristic curves. In pipe flow, two characteristic lines ($C^+$ and $C^-$) propagate at wave speeds $V \pm c$ and carry compatibility equations.

Figure: MOC stencil — C⁺ and C⁻ lines converge at unknown point P

Finite Difference Form in terms of H and Q:

Using piezometric head $H = P/\gamma + z$ and $B = c/(gA)$, the characteristic equations become:

Solving for an interior point P: Set the two expressions for $H_P$ equal:

Solution:

We are simulating a 1D aquifer using 3 grids. Grid size $\Delta x = 400$m. River stage $h_L = 365$m is the left boundary condition for Grid 1. Right boundary of Grid 3 is a barrier (no flow).

Figure: 1D aquifer grid — 3 cells between river (hL=365m) and barrier

1. Calculate Grid Coefficients

For Grid 1 (Adjacent to river): $A_1 = \frac{2T}{\Delta x(\Delta x + \Delta x)} = \frac{2(600)}{400(800)} = 0.00375$.
For flow from the river, $C_1 = \frac{2T}{\Delta x^2} = \frac{2(600)}{400^2} = 0.0075$.

For Grid 2 (Interior): $A_2 = C_2 = \frac{2T}{2\Delta x^2} = 0.00375$.

For Grid 3 (Adjacent to barrier): $C_3 = 0.00375$. Since it's a barrier on the right, $A_3 = 0$.

2. Setup Implicit Tridiagonal Equations
General equation: $C_i h_{i-1}^+ + E_i h_i^+ + A_i h_{i+1}^+ = Q_i$
where $E_i = -(A_i + C_i + S/\Delta t)$ and $Q_i = -(S/\Delta t) h_i$.

Grid 1:
$E_1 = -(0.00375 + 0.0075 + 0.15) = -0.16125$
$Q_1 = -0.15 \times 349.13 = -52.3695$
Equation 1: $C_1 h_L + E_1 h_1^+ + A_1 h_2^+ = Q_1$
$0.0075(365) - 0.16125 h_1^+ + 0.00375 h_2^+ = -52.3695$
$-0.16125 h_1^+ + 0.00375 h_2^+ = -55.107$ (Eq. 1)

Grid 2:
$E_2 = -(0.00375 + 0.00375 + 0.15) = -0.1575$
$Q_2 = -0.15 \times 347.97 = -52.1955$
$0.00375 h_1^+ - 0.1575 h_2^+ + 0.00375 h_3^+ = -52.1955$ (Eq. 2)

Grid 3:
$E_3 = -(0 + 0.00375 + 0.15) = -0.15375$
$Q_3 = -0.15 \times 339.36 = -50.904$
$0.00375 h_2^+ - 0.15375 h_3^+ = -50.904$ (Eq. 3)

3. Solving the System of Equations

From (Eq. 3): $h_3^+ = \frac{0.00375 h_2^+ + 50.904}{0.15375} = 0.02439 h_2^+ + 331.083$

Substitute $h_3^+$ into (Eq. 2):
$0.00375 h_1^+ - 0.1575 h_2^+ + 0.00375(0.02439 h_2^+ + 331.083) = -52.1955$
$0.00375 h_1^+ - 0.1575 h_2^+ + 0.000091 h_2^+ + 1.2415 = -52.1955$
$0.00375 h_1^+ - 0.1574 h_2^+ = -53.437$
$h_1^+ = \frac{0.1574 h_2^+ - 53.437}{0.00375} = 41.973 h_2^+ - 14249.8$

Substitute $h_1^+$ into (Eq. 1):
$-0.16125 (41.973 h_2^+ - 14249.8) + 0.00375 h_2^+ = -55.107$
$-6.768 h_2^+ + 2297.53 + 0.00375 h_2^+ = -55.107$
$-6.764 h_2^+ = -2352.637$
$h_2^+ \approx 347.82 \text{ m}$

Back-substitute to find others:
$h_3^+ = 0.02439(347.82) + 331.083 \approx 339.56 \text{ m}$
$h_1^+ = 41.973(347.82) - 14249.8 \approx 349.21 \text{ m}$

Water Table Elevations at end of Day 1:

Solution:

Finite Difference Method (FDM) Concept:
The Finite Difference Method is a numerical technique for solving partial differential equations (PDEs) by replacing the continuous spatial and temporal domains with a discrete array of grid points (rows and columns). The differential equation is approximated point-wise by converting it into algebraic difference equations using Taylor series expansions. The value of a function at a specific grid point is algebraically related to the values of the function at neighboring points.

Figure: Comparison of Explicit and Implicit finite difference schemes

Explicit Scheme:
In an explicit scheme, the solution at the next time step ($t = n+1$) is calculated solely from the known conditions at the current time step ($t = n$).

• The equations can be solved sequentially, point by point.
• Simpler to implement mathematically.
• Subject to stability constraints (e.g., Courant condition) which force the use of very small time steps.

Implicit Scheme:
In an implicit scheme, the solution at the next time step ($t = n+1$) is formulated using conditions at both the current time step ($t = n$) and the unknown next time step ($t = n+1$).

• Requires solving a system of algebraic equations simultaneously (often involving matrix inversions).
• Mathematically more complex.
• Unconditionally stable for most linear cases, allowing the use of much larger computation time steps.