Numerical Methods (Pipe Flow)

The pipe flow solver utilizes the same semi-implicit, finite-volume algorithms used in the 2D solver and the 1D finite volume solver. These methods are detailed for the 1D finite volume solver in Numerical Methods (1D FV). In short, velocities in the continuity equation and the pressure gradient term in the momentum equation are discretized using the $\begin{array}{l}\theta\end{array}$ -method. Frictional terms are treated implicitly, and the advection term is treated explicitly using the Eulerian-Lagrangian method. The equations are discretized on the 1D unstructured pipe network grid. Cell property tables relate water levels to cell volume and surface area, and face property tables relate water levels to cross-sectional and wetted areas. Equations for each computational cell are aggregated into matrix form and solved using an iterative, Newton's method-like solution algorithm.

Two major differences are present in the pipe network solution, the first in the discretization of the advection term and the second in the matrix solution method.

Advection Discretization

Acceleration terms at a computational face, $\begin{array}{l}j\end{array}$ , are discretized using a Lagrangian approach.

$\begin{array}{l}\displaystyle \dfrac{\partial V_j}{\partial t} + V_j \dfrac{\partial V_j}{\partial x} \approx \dfrac{V_j^{n+1} - V_{X(j)}^n }{\Delta t}\end{array}$

where $\begin{array}{l}V_{X(j)}^n\end{array}$ is the velocity at the start of the time step evaluated at a location, $\begin{array}{l}X(j)\end{array}$ . This location is found by integrating the velocity field back in time, starting from the face location. Since $\begin{array}{l}X(j)\end{array}$ does not necessarily correspond to a cross-section, an interpolation technique is applied. For a positive velocity at the face

$\begin{array}{l}\displaystyle V_{X(j)}^n = V_j^n - V_j^* \Delta t \dfrac{V_j^n - V_{j-1}^n}{x_j - x_{j-1}}\end{array}$

where $\begin{array}{l}V_j^*\end{array}$ is the advective speed, and $\begin{array}{l}V_j^* \Delta t\end{array}$ is the distance traveled back into the upwind cell.

In areas of rapidly varying flow conditions, the choice of the advective velocity is important (Stelling and Duinmeijer, 2003). In strong flow contractions, a discretization consistent with constancy of energy head is applied.

$\begin{array}{l}\displaystyle V_j^* = \dfrac{V_j^n + V_{j-1}^n}{2}\end{array}$

In areas without strong contractions (e.g., flow expansions), a discretization consistent with momentum conservation is applied instead.

$\begin{array}{l}\displaystyle V_j^* = \dfrac{A_j^n V_j^n + A_{j-1}^n V_{j-1}^n}{2A_j^n}\end{array}$

A dynamic choice between the two discretizations is made based on local flow conditions near the face. This prevents unphysical increases in energy through flow contractions, and unphysical conservation of energy through flow expansions.

Matrix Solution Technique

The second difference involves the Newton's method-like matrix solution technique, and it arise from the decomposition of the cell property tables into positive and negative contributions. Using the cell volume decomposition, equation 10) in Numerical Methods (1D FV) is rewritten:

$\begin{array}{l}\displaystyle \boldsymbol{\Omega}^{+}(\boldsymbol{H}) + \boldsymbol{\Omega}^{-}(\boldsymbol{H}) + \boldsymbol{\Psi}\boldsymbol{H} = \boldsymbol{b}\end{array}$

where

$\begin{array}{l}\boldsymbol{\Omega}^+\end{array}$ and $\begin{array}{l}\boldsymbol{\Omega}^-\end{array}$ are vectors of the the decomposed positive and negative volume contributions of each cell,
$\begin{array}{l}\boldsymbol{H}\end{array}$ is a vector of the unknown hydraulic heads (water surface elevation + pressure head) for each cell,
$\begin{array}{l}\boldsymbol{\Psi}\end{array}$ is the coefficient matrix from the semi-implicit discretization of the continuity and momentum equations, and
$\begin{array}{l}\boldsymbol{b}\end{array}$ is the right hand side vector.

In the methods outlined by Casulli and Stelling (2013), the negative volume contribution is linearized with a series of iterates, $\begin{array}{l}m\end{array}$ ,

$\begin{array}{l}\displaystyle \boldsymbol{\Omega}^{-}(\boldsymbol{H}^{m+1}) = \boldsymbol{\Omega}^{-}(\boldsymbol{H}^m) + \boldsymbol{A}^{-}(\boldsymbol{H}^m})(\boldsymbol{H}^{m+1} - \boldsymbol{H}^m)\end{array}$

and the positive volume contribution can be linearized with a series of iterates, $\begin{array}{l}k\end{array}$ ,

$\begin{array}{l}\displaystyle \boldsymbol{\Omega}^{+}(\boldsymbol{H}^{m+1,k+1}) = \boldsymbol{\Omega}^{+}(\boldsymbol{H}^{m+1,k}) + \boldsymbol{A}^{+}(\boldsymbol{H}^{m+1,k})(\boldsymbol{H}^{m+1,k+1} - \boldsymbol{H}^{m+1,k})\end{array}$

where

$\begin{array}{l}\boldsymbol{A}^+\end{array}$ and $\begin{array}{l}\boldsymbol{A}^-\end{array}$ are vectors of the the decomposed positive and negative wetted area contributions of each cell.

This leads to a set of nested iteration loops, where an initial guess at $\begin{array}{l}H\end{array}$ is iterated on until convergence in the $\begin{array}{l}k\end{array}$ iterates, and then the $\begin{array}{l}m\end{array}$ iterates are updated, and the $\begin{array}{l}k\end{array}$ iterates are run again. This process continues until both the $\begin{array}{l}m\end{array}$ and $\begin{array}{l}k\end{array}$ iterates have converged to within a given tolerance. Both iteration loops are formulated to solve for the differences in successive iterations (e.g., $\begin{array}{l}\boldsymbol{H}^{m+1} - \boldsymbol{H}^m\end{array}$ ) as shown below.

$\begin{array}{rl} \left[ \boldsymbol{\Psi} + \boldsymbol{A}^+(\boldsymbol{H}^m) - \boldsymbol{A}^-(\boldsymbol{H}^m) \right] ( \boldsymbol{H}^{m+1}-\boldsymbol{H}^m ) &= -( \boldsymbol{\Omega}^+(\boldsymbol{H}^m) + \boldsymbol{\Psi} \boldsymbol{H}^m - \boldsymbol{b} ) + \boldsymbol{V}^-(\boldsymbol{H}^m) \\ \left[ \boldsymbol{\Psi} + \boldsymbol{A}^+(\boldsymbol{H}^{m+1,k}) - \boldsymbol{A}^-(\boldsymbol{H}^m) \right] ( \boldsymbol{H}^{m+1,k+1}-\boldsymbol{H}^{m+1,k} ) &= -( \boldsymbol{\Omega}^+(\boldsymbol{H}^{m+1,k}) + \boldsymbol{\Psi} \boldsymbol{H}^{m+1,k} - \boldsymbol{b} ) + \boldsymbol{V}^-(\boldsymbol{H}^m) + \boldsymbol{A}^-(\boldsymbol{H}^m)(\boldsymbol{H}^{m+1,k}-\boldsymbol{H}^m) \end{array}$

References

Casulli, V. and Stelling, G.S. 2013. A semi-implicit numerical model for urban drainage systems. International Journal for Numerical Methods in Fluids. 73:600-614.
Stelling, G.S. and Duinmeijer, S.P.A. 2003. A staggered conservative scheme for every Froude number in rapidly varied shallow water flows. International Journal for Numerical Methods in Fluids. 43:1329-1354.