Full Saint-Venant explicit scheme numerical (50 m)

A simple explicit scheme is the Lax diffusive scheme: forward in time, central in space, with the old value at $i$ replaced by the average of its two neighbours. The equations for a rectangular channel are

$$B\frac{\partial y}{\partial t}+\frac{\partial Q}{\partial x}=0,\qquad \frac{\partial Q}{\partial t}+\frac{\partial}{\partial x}\left(\frac{Q^2}{A}\right)+gA\frac{\partial y}{\partial x}-gA(S_0-S_f)=0$$

Finite difference form (grid $(i,n+1)$ from $(i-1,n),(i,n),(i+1,n)$):

$$y_i^{n+1}=\frac{y_{i+1}^n+y_{i-1}^n}{2}-\frac{\Delta t}{2B\Delta x}\left(Q_{i+1}^n-Q_{i-1}^n\right)$$ $$Q_i^{n+1}=\frac{Q_{i+1}^n+Q_{i-1}^n}{2}-\frac{\Delta t}{2\Delta x}\left[\frac{(Q^2)_{i+1}}{A_{i+1}}-\frac{(Q^2)_{i-1}}{A_{i-1}}\right]-gA_i\frac{\Delta t}{2\Delta x}\left(y_{i+1}-y_{i-1}\right)+gA_i\Delta t\,(S_0-S_{f,i})$$

Data

$B=50$ m, $S_0=0.0002$, $n=0.035$, $\Delta x=1000$ m, $\Delta t=600$ s, so $\dfrac{\Delta t}{2\Delta x}=0.300$.

Step 1: Areas and friction slope

$A_{i-1}=93.50$, $A_i=91.00$, $A_{i+1}=98.50$ m$^2$.

$$P_i=B+2y_i=53.64\ \text{m},\quad R_i=\frac{A_i}{P_i}=1.6965\ \text{m},\quad S_{f,i}=\frac{n^2Q_i^2}{A_i^2R_i^{4/3}}=9.475\times10^{-5}$$

Step 2: Flow depth

$$y_i^{n+1}=\frac{1.97+1.87}{2}-0.00600\,(35-38)=1.9200-(-0.0180)=1.9380\ \text{m}$$

Step 3: Discharge

• Average term: $\dfrac{35+38}{2}=36.500$
• Convective term: $\dfrac{Q^2}{A}$ at $i+1$ is 12.437, at $i-1$ is 15.444; difference $=-3.0073$; term $=-0.300(-3.0073)=0.902$
• Pressure term: $-gA_i\dfrac{\Delta t}{2\Delta x}(y_{i+1}-y_{i-1})=-9.81(91.00)(0.300)(0.10)=-26.781$
• Gravity and friction: $gA_i\Delta t(S_0-S_f)=9.81(91.00)(600)(0.0002-9.475\times10^{-5})=56.374$

$$Q_i^{n+1}=36.500+(0.902)+(-26.781)+(56.374)=66.99\ \text{m}^3/\text{s}$$

Answer: $y_i^{n+1}=1.9380$ m and $Q_i^{n+1}=66.99$ m$^3$/s.

Remark: the large jump in $Q$ comes mainly from the gravity-friction term ($56.4$): the given flow is far from uniform ($S_f$ is less than half of $S_0$), and this imbalance is multiplied by a long time step. Also, $\Delta t=600$ s gives a Courant number $(V+\sqrt{gy})\Delta t/\Delta x\approx2.8>1$, so the explicit scheme is outside its stability limit; a time step below about 216 s would be needed to continue the simulation. One step is computed here, as asked.