Discrete Scheme for SWE

In the ELM-SWE solver, a semi-Lagrangian approach is used to discretize the acceleration terms in the momentum equation. While this approach has the advantage of being stable for large time steps, it can create excessive numerical diffusion of momentum, leading to inaccurate results in lab-scale simulations where strict conservation of momentum is important. For this reason, an alternative Eulerian SWE solver (EM-SWE) is provided. The alternative approach utilizes the momentum-conservative discretization of the acceleration terms suggested by Kramer and Stelling (2008).

Advection

In Kramer and Stelling’s approach, the advective term in the momentum equation is discretized assuming local conservation of momentum about a control volume centered on a cell face, $\begin{array}{l}(\boldsymbol{V} \cdot \nabla ) u_N\end{array}$ , shown in Figure 1 (Perot, 2000).

1)

$\begin{array}{l}\displaystyle \left[ (\boldsymbol{V} \cdot \nabla ) u_N \right]_f \approx \left( \frac{1}{h} \left [ ( \nabla \cdot ( h \boldsymbol{V} u_N ) - u_N \nabla \cdot ( h \boldsymbol{V} ) \right ] \right) _f \approx \\ \displaystyle \frac{\alpha ^L _f}{\overline{h} _f A_L} \displaystyle\sum_{k \in \text{K} ( L ) } s_{L,k} Q_k \left( \boldsymbol{V}^u_k \cdot \textbf{n} _f - u_{N,f} \right) + \frac{\alpha ^R _f}{\overline{h} _f A_R} \displaystyle\sum_{k \in \text{K} ( R ) } s_{R,k} Q_k \left( \boldsymbol{V}^u_k \cdot \textbf{n} _f - u_{N,f} \right)\end{array}$

Figure 1. Definition of Eulerian advection variables.

Variables with $\begin{array}{l}L=\text{L} (f)\end{array}$ indicate quantities for the cell to the left of face $\begin{array}{l}f\end{array}$ and those with $\begin{array}{l}R=\text{R} (f)\end{array}$ indicate the right cell. The weights $\begin{array}{l}\alpha ^L _f\end{array}$ and $\begin{array}{l}\alpha ^R _f\end{array}$ are given to the left and right cells and are based on the area of the control volume.

2)	$\begin{array}{l}\displaystyle \displaystyle \alpha ^L_f = \frac{\Delta x^L_f}{\Delta x^L_f + \Delta x^R_f}, \: \alpha ^R_f = 1- \alpha ^L_f\end{array}$

The water depth at the face () is calculated as a weighted average of the average water depths of the left and right cells:

3)	$\begin{array}{l}\displaystyle \bar{h}_f = \alpha ^R_f h_R + \alpha ^L_f h_L\end{array}$

where $\begin{array}{l}h_{L|R} = \Omega_{L|R} / A_{L|R}\end{array}$ . For both the left and right cells, the flux of momentum into the cell is calculated using upwinded velocities from each edge ( $\begin{array}{l}\boldsymbol{V}_k^u\end{array}$ ) in the direction of the face normal ( $\begin{array}{l}\textbf{n}_k\end{array}$ ) and the face flows $\begin{array}{l}Q_k = A_k u_{N,k}\end{array}$ . Upwinded velocities are taken at the cell centers upwind of the faces and are reconstructed from the face normal velocities using the weighting method of Perot (2000) (Cell Velocity). Because this method is explicit, all acceleration terms can be computed before the start of the outer iterations and do not need to be updated during the iterations.

Coriolis Effect

Discretizations of the Coriolis and eddy viscosity terms in the original SW solver required variables to be evaluated at the backtracked location X. Since this acceleration discretization requires no backtracking of velocities, an alternative discretization of these terms are required. The Coriolis term is discretized with a simple explicit treatment:

4)	$\begin{array}{l}\displaystyle ( f_c \boldsymbol{k} \times \boldsymbol{V}) \cdot \textbf{n}_f \approx - f_c u_{T,f}^n\end{array}$

where $\begin{array}{l}u_{T,f}^n\end{array}$ is the tangential velocity at face k, reconstructed using the normal face velocities from adjacent cells, and $\begin{array}{l}f_c\end{array}$ is the Coriolis parameter.

Barotropic Pressure Gradient

In the EM-SWE model, the barotropic pressure gradient term is treated similar the ELM-SWE model. The term is computed at computational faces utilizing the two-point stencil described above and treats the water levels semi-implicitly. The term is may be written as

5)	$\begin{array}{l}\displaystyle g \nabla z_s \cdot \textbf{n}_f = g \frac{h_f}{\bar{h}_f} \frac{\partial z_s^{n+\theta} }{\partial N}\end{array}$

where , , and is the face hydraulic depth. The pressure gradient term in EM-SWE differs from ELM-SWE in the inclusion of the ratio following Kramer and Stelling (2008). It is noted that this factor is not computed the same here as in Kramer and Stelling (2008) since HEC-RAS uses subgrid bathymetry. However, this adaptation has shown to work well in HEC-RAS.

Solution Procedure

The semi-discrete form of the second fractional step may be written as

6)

$\begin{array}{l}\displaystyle \displaystyle \frac{u_{N,f}^{n+1} - u^n_{N,f}}{\Delta t} + \left [ \left ( \boldsymbol{V} \cdot \nabla \right ) u_N \right ]_f^n = -g \frac{h_f}{\bar{h}_f} \frac{\partial z_{s,f}^{n + \theta}}{\partial N} + \left[ \frac{1}{h} \nabla \cdot \left( \nu_{t} h \nabla u_N \right) \right]_f^n - c_{f,f} u_{N,f}^{n+1} + \frac{\tau _{s,N}}{\rho h_f^n}\end{array}$

The momentum equation above can be rearranged to obtain an expression for the velocity at n+1 as

7)	$\begin{array}{l}\displaystyle \displaystyle u_{N,f}^{n+1} = - \frac{\Delta t \theta g} {1 + \Delta t c_{f,f} } \frac{h_f}{\bar{h}_f} \frac{\partial z_{s,f}^{n+1}} {\partial N} + F_f\end{array}$

where

$\begin{array}{l}F_f = \frac{ B_f^n} {1 + \Delta t c_{f,f} }\end{array}$
$\begin{array}{l}B_f^n = u_{N,f}^n - (1-\theta) \Delta t g \frac{h_f } { \bar{h}_f } \frac{\partial z_s^n } {\partial N} - \Delta t \left [ \left ( \boldsymbol{V} \cdot \nabla \right ) u_N \right ]_f^n + \Delta t \left[ \frac{1}{h} \nabla \cdot \left( \nu_{t} h \nabla u_N \right) \right]_f^n + \Delta t \frac{ \tau_{s,N} } {\rho h_f^n}\end{array}$

To obtain a discrete implicit equation for the water volume, the above equation is inserted into the discrete continuity equation to obtain

8)	$\begin{array}{l}\displaystyle \displaystyle \Omega_i^{n+1} + \sum_{j \in \text C (i)} a_{i,j} z_{s,j}^{n+1} = b_i\end{array}$

where

$\begin{array}{l}a_{i,j} = - \frac{ \Delta t^2 \theta^2 g A_k h_k } { \left ( 1+ \Delta t c_{f,f} \right ) \Delta x_N \bar{h}_k }, j \in \text{N}(i)\end{array}$
$\begin{array}{l}a_{i,i} = - \sum_{j \in \text{N} (i)} a_{i,j}\end{array}$
$\begin{array}{l}b_i = \Omega_i^n + \Delta t Q_i + \frac{1 - \theta}{\theta} \sum_{j \in \text{C}(i)} a_{i,j} z_{s,j}^n + \Delta t \sum_{k \in \text{K}(i)} \left [ \theta F_k + (1-\theta) u_{N,k}^n \right ] A_k\end{array}$

The system of equations is solved using the same Newton-like iterations used for the DWE solver.

Robustness and Stability

The SWE-EM solver has improved momentum conservation compared to the SWE-ELM solver. However, the tradeoff for more accurate momentum conservation is that the method requires the 2D grid be strictly orthogonal, and the time step necessary for stability is limited by the Courant-Friedrichs-Lewy (CFL) condition:

9)	$\begin{array}{l}\displaystyle \displaystyle C = \frac{u \Delta x}{\Delta t} \leq C_{max}\end{array}$

Where $\begin{array}{l}C\end{array}$ is the Courant number, and $\begin{array}{l}C_{max}\end{array}$ is 1. If turbulence (momentum diffusion) is turned on, the explicit treatment of the momentum diffusion results in an additional stability criteria which is approximated as,

10)	$\begin{array}{l}\displaystyle \displaystyle \Delta t < \frac{\Delta x^2}{2v_t}\end{array}$

This CLF condition allows a larger time step than typical stability conditions originating from Eulerian advection schemes.

The linearized scheme is second order accurate in space. The time accuracy depends of the choice of ; for instance, for $\begin{array}{l}\theta\end{array}$ =1 it is first order accurate and for $\begin{array}{l}\theta\end{array}$ =1/2 it is second order accurate.

Discrete Boundary Conditions

The boundary conditions are discretized in manner as in the SWE-ELM solver as:

Water surface elevation: The water surface elevation boundary condition is directly implemented as . The internal cell is then discretized as described above and the terms containing the boundary water surfaces are placed on the right-hand-side of the system of equations.
Normal Depth: The friction slope, , is specified and utilized to compute a flow at each computation face as . Boundary face flows are included in the internal cells as a source term on the right-hand-side of the system of equations.
Flow: The flow boundary condition is specified at each computational face based on the local conveyance. Boundary face flows are included in the internal cells as a source term on the right-hand-side of the system of equations.

Solution Algorithm.

The solution algorithm of the EM-SWE solver proceeds exactly as described in the SWE-ELM solver. However, because of the explicit treatment of the acceleration terms, the time step necessary for stability is limited by the CFL condition in Equation 2-183.

Download PDF

Eulerian Shallow Water Equation Solver