Full Saint-Venant explicit scheme numerical (10 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=10$ m, $S_0=0.0002$, $n=0.04$, $\Delta x=1000$ m, $\Delta t=300$ s, so $\dfrac{\Delta t}{2\Delta x}=0.150$.

Step 1: Areas and friction slope

$A_{i-1}=19.00$, $A_i=18.50$, $A_{i+1}=20.00$ m$^2$.

$$P_i=B+2y_i=13.70\ \text{m},\quad R_i=\frac{A_i}{P_i}=1.3504\ \text{m},\quad S_{f,i}=\frac{n^2Q_i^2}{A_i^2R_i^{4/3}}=0.00452$$

Step 2: Flow depth

$$y_i^{n+1}=\frac{2.0+1.9}{2}-0.01500\,(37.5-40)=1.9500-(-0.0375)=1.9875\ \text{m}$$

Step 3: Discharge

• Average term: $\dfrac{37.5+40}{2}=38.750$
• Convective term: $\dfrac{Q^2}{A}$ at $i+1$ is 70.312, at $i-1$ is 84.211; difference $=-13.8980$; term $=-0.150(-13.8980)=2.085$
• Pressure term: $-gA_i\dfrac{\Delta t}{2\Delta x}(y_{i+1}-y_{i-1})=-9.81(18.50)(0.150)(0.10)=-2.722$
• Gravity and friction: $gA_i\Delta t(S_0-S_f)=9.81(18.50)(300)(0.0002-0.00452)=-235.358$

$$Q_i^{n+1}=38.750+(2.085)+(-2.722)+(-235.358)=-197.25\ \text{m}^3/\text{s}$$

Answer: $y_i^{n+1}=1.9875$ m and $Q_i^{n+1}=-197.25$ m$^3$/s.

Remark: the huge change in $Q$ (even to a negative value) comes from the friction term: $S_f=0.00452$ is about 23 times $S_0$, so the given flow is far from uniform, and this imbalance is multiplied by a long time step ($-235.4$). Also, $\Delta t=300$ s gives a Courant number $(V+\sqrt{gy})\Delta t/\Delta x\approx1.9>1$, so the explicit scheme is outside its stability limit; a time step below about 158 s would be needed to continue the simulation. A negative discharge has no physical meaning here; it is simply the one-step result of the scheme for these data.