Download PDF
Download page Eulerian Shallow Water Equation Solver.
Eulerian Shallow Water Equation Solver
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 (2008) approach, the advective term in the momentum equation is discretized assuming local conservation of momentum about a control volume centered on a cell face, (\boldsymbol{V} \cdot \nabla ) u_N, shown in the Figure below.
Figure 1. Definition of Eulerian advection variables.
The advection term is given by
1) | \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 ) } Q^-_{L,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 ) } Q^-_{R,k} \left( \boldsymbol{V}^u_k \cdot \textbf{n} _f - u_{N,f} \right) |
where
h_i = \Omega_i / A_i^W
Q^-_{i,k} = \min \left ( s_{i,k} Q_k,0 \right ) : inflow at face k to cell i
A_k : face vertical area
A_{\text{L}} : left cell horizontal area
A_{\text{R}} : right cell horizontal area
\textbf{V}_j : cell-average current velocity vector
\textbf{V}_k^u = \textbf{V}_{L} \,\, \text{for} \,\, u_{N,k} > 0
\textbf{V}_k^u = \textbf{V}_{R(k)} \,\, \text{for} \,\, u_{N,k} < 0
\textbf{n}_f : face-normal unit vector
Variables with L=\text{L} (f) indicate quantities for the cell to the left of face f and those with R=\text{R} (f) indicate the right cell. The weights \alpha ^L _f and \alpha ^R _f are given to the left and right cells and are based on the area of the control volume.
2) | \displaystyle \alpha ^L_f = \frac{\Delta x^L_f}{\Delta x^L_f + \Delta x^R_f}, \: \alpha ^R_f = 1- \alpha ^L_f |
The water depth at the face () is calculated as a weighted average of the average water depths of the left and right cells:
3) | \bar{h}_f = \alpha ^R_f h_R + \alpha ^L_f h_L |
where h_{L|R} = \Omega_{L|R} / A_{L|R} . For both the left and right cells, the flux of momentum into the cell is calculated using upwinded velocities from each edge (\boldsymbol{V}_k^u ) in the direction of the face normal (\textbf{n}_k ) and the face flows Q_k = A_k u_{N,k}. 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) | ( f_c \boldsymbol{k} \times \boldsymbol{V}) \cdot \textbf{n}_f \approx - f_c u_{T,f}^n |
where u_{T,f}^n is the tangential velocity at face k, reconstructed using the normal face velocities from adjacent cells, and f_c 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) | g \nabla z_s \cdot \textbf{n}_f = g \frac{\partial z_s^{n+\theta} }{\partial N} |
where z_s is the water surface elevation and N is the face-normal direction.
Solution Procedure
The semi-discrete form of the second fractional step may be written as
6) | \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{\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} |
The momentum equation above can be rearranged to obtain an expression for the velocity at n+1 as
7) | \displaystyle u_{N,f}^{n+1} = - \frac{\Delta t \theta g} {1 + \Delta t c_{f,f} } \frac{\partial z_{s,f}^{n+1}} {\partial N} + F_f |
where
F_f = \frac{ B_f^n} {1 + \Delta t c_{f,f} }
B_f^n = u_{N,f}^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 \left ( \frac{ \tau_{s,N} } {\rho h_f^n} - \frac{1}{\rho g} \frac{\partial p_a}{\partial N} \right )
To obtain a discrete implicit equation for the water volume, the above equation is inserted into the discrete continuity equation to obtain
8) | \displaystyle \Omega_i^{n+1} + \sum_{j \in \text C (i)} a_{i,j} z_{s,j}^{n+1} = b_i |
where
a_{i,j} = - \frac{ \Delta t^2 \theta^2 g A_k } { \left ( 1+ \Delta t c_{f,f} \right ) \Delta x_N }, j \in \text{N}(i)
a_{i,i} = - \sum_{j \in \text{N} (i)} a_{i,j}
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
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) | \displaystyle C = \frac{u \Delta x}{\Delta t} \leq C_{max} |
Where C is the Courant number, and C_{max} 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) | \displaystyle \Delta t < \frac{\Delta x^2}{2v_t} |
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 \theta=1 it is first order accurate and for \theta =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.