Detailed 2D Bridge Modeling

Overview

The detailed modeling approach utilizes multiple rows of 2D cells within the bridge footprint in order to resolve the bridge hydraulics. It transitions between unpressurized flow at low flow conditions, to pressurized flow when the bottom of the bridge deck is wetted, to 2-layer flow when the bridge is overtopped. The approach is more computationally expensive and requires detailed geometry of the bridge, but is more accurate and does not require as much parameterization as the simplified approach.

Note

The detailed bridge modeling approach does NOT automatically modify the terrain to account for the bridge geometry. Therefore, the user must make sure the bridge geometry with the exception of the bridge deck is included in the terrain.

Low Flow Below-Deck

Low flow below the bridge deck is unpressurized. This type of flow is treated the same as any other 2D flow simulation. As noted above, it is important that the bridge geometry including abutments and piers be included in the terrain and NOT in the bridge geometry. This is a major difference between the 1D and simplified 1D/2D bridge approaches, and the detailed 2D bridge approach.

Pressurized Flow Below-Deck

Two-dimensional cells beneath the bridge deck become pressurized as water levels increase past the height of the bottom of the bridge deck. At that point, the computed water surface elevation represents the sum of the height of the bottom of the bridge deck plus the pressure head, $\begin{array}{l}P/\gamma\end{array}$ , where $\begin{array}{l}\gamma\end{array}$ is the specific weight of water.

Once pressurized, wave speeds become infinitely fast (when water compressibility and wall elasticity are ignored) and explicit solvers suffer from prohibitively large time step restrictions. The semi-implicit approach taken in the RAS solver treats the pressure term in the momentum equation implicitly, and is not subject to the same restrictions. Modifications to the semi-implicit computational engine were still necessary to enable modeling of pressurized flow. These include updates to the solution algorithm and the cell property tables and are described in detail in V. Casulli and G. Stelling's 2013 paper "A semi-implicit numerical model for urban drainage systems" (International Journal for Numerical Methods in Fluids, volume 73). The extension of these methods from 1D to 2D, as well as their adaptation to the RAS computational scheme is described in the following sections.

Decomposition of Hydraulic Property Tables

Flow becomes pressurized when rising water encounters a ceiling in the bathymetry. In contrast to open channel flow, this means that the wet surface area within a cell (the derivative of the cell volume) decreases with increasing water surface elevation. This cannot be handled directly in the original implementation of the semi-implicit solver because the Newton's Method-like algorithm requires a non-decreasing surface area-elevation function for stability. Casulli and Stelling showed that this can be handled by decomposing the volume- and area-elevation curves into two non-decreasing functions: one for the positive contribution to volume or area, and one for the negative contribution. Each volume or area contribution can then be linearized independently and solved using a set of nested Newton-like iterations.

The figure below shows an example of a 2D computational cell, with one square unit of surface area and a ceiling one unit high. The total cell surface area-elevation curve (top left plot) is decomposed into positive (top middle plot) and negative (top right plot) contributions. The positive area contribution is constant with elevation through the height of the ceiling. At those elevations, the negative surface area contribution is zero. Above the ceiling, the positive contribution continues at 1.0 (to remain non-decreasing) and the negative surface area contribution jumps to 1.0, so that the total area sums to zero. The volume-elevation curves (bottom row of plots) area the integral of the area curves.

Formally, the decomposition of a function into positive and negative contributions is known as a Jordan decomposition and can be expressed:

$\begin{array}{l}\displaystyle \Omega(H) = \Omega^{+}(H) - \Omega^{-}(H), \hspace{10pt} A(H) = A^{+}(H) - A^{-}(H)\end{array}$

It can be performed by moving through the property table from the lowest to highest elevations and calculating

$\begin{array}{l}\displaystyle A_k^{+} = A_k^{+} + \max(A_k - A_{k-1}, 0), \hspace{10pt} A_k^{-} = A_k^{-} - \min(A_k - A_{k-1}, 0) \newline \Omega_k^{+} = \Omega_{k-1}^{+} + (H_k - H_{k-1}) A_{k-1}^{+}, \hspace{10pt} \Omega_k^{-} = \Omega_{k-1}^{-} + (H_k - H_{k-1}) A_{k-1}^{-}\end{array}$

with $\begin{array}{l}A_0^{+} = A_0\end{array}$ and $\begin{array}{l}A_0^{-} = 0\end{array}$ .

Solver Extension for Pressurized Flow

In 2D open channel flow, the discretizations of the continuity and momentum equations on an unstructured grid result in a matrix equation of the form

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

where

$\begin{array}{l}\boldsymbol{\Omega}\end{array}$ is a vector of volume for each cell,
$\begin{array}{l}\boldsymbol{Z}\end{array}$ is a vector of the unknown water surface elevations 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.

(See Eulerian-Lagrangian Shallow Water Equation Solver, Equation 12).

With consideration of pressurized flow and using the cell volume decomposition, the system of equations can be written

$\begin{array}{l}\displaystyle \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, and
$\begin{array}{l}\boldsymbol{H}\end{array}$ is a vector of the unknown hydraulic heads (water surface elevation + pressure head) for each cell.

The mildly non-linear system is solved by linearizing both positive and negative volume contributions, and solving using a series of nested iterations. That procedure is the same one that is used in the 1D finite-volume pressurized flow solver for pipe networks, and is described in Numerical Methods (Pipe Flow).

Bridge Overtopping

Cloned Cells

When the detailed 2D bridge modeling approach is selected, the set of computational cells and cell faces beneath the footprint of the bridge are flagged. A set of co-located cells and cell faces is then created within the computational engine to calculate and store hydrodynamic state variables for flow above the bridge deck - these cells are referred to as cloned cells, following the terminology of V. Casulli ("Computational grid, subgrid, and pixels", Numerical Methods in Fluids, 2019). The approach of cloning and connecting co-located cells to a 2D grid was shown by Casulli to work successfully with the semi-implicit solution structure.

Cloned cells and faces are given their own unique indices. The property tables - originally computed using only information for the terrain - are modified so that the below-deck cells are given a ceiling at the elevation of the bottom of the bridge deck, and the above-deck cells have a minimum elevation of the bridge deck surface. Modifications to the spatial discretizations for computing average values and gradients at the edge of the bridge were made to account for the connectivity of the added cloned cells.

Overhead view of river flow through a 2D bridge structure within a HEC-RAS 2D computational area.

Oblique view of the connectivity of the cloned cells used to represent flow over the 2D bridge

Plunging Flow Off Bridge Deck

Plunging flow off the downstream edge of a bridge deck can occur when the bridge is overtopped and the downstream water surface is below the bridge deck height. This is handled in HEC-RAS by modifying the discretization of the pressure term in the momentum equation at faces on the edge of the bridge. The approach closely follows that of D-Flow FM, and uses an energy conservation (i.e., Bernoulli equation) approach between the upstream and downstream cells of the bridge edge. In this approach, the water depth at the edge is limited to

$\begin{array}{l}\displaystyle \displaystyle h_{f} = \max ( z_{s,D} - z_{b,f}, \frac{2}{3} E_f )\end{array}$

where

$\begin{array}{l}z_{s,D}, z_{s,U}\end{array}$ are the surface water elevations in the downstream and upstream cells,
$\begin{array}{l}z_{b,f}\end{array}$ is the bed elevation at the bridge edge face,
$\begin{array}{l}E_f = z_{s,U} - z_{b,f} + \dfrac{( \tilde{\textbf{V}}^u_k \cdot \textbf{n} _f )^2}{2g}\end{array}$ is the far-field energy head above the bridge edge, and
$\begin{array}{l}\tilde{\textbf{V}}^u_k\end{array}$ is the upstream cell velocity vector.

The Bernoulli equation is used to derive a correction term to the momentum equation for energy conservation given by

$\begin{array}{l}\displaystyle \displaystyle \mathcal{H}_f = \frac{g}{\Delta x} \max \left [ 0,\frac{2}{3} E_f -(z_{s,D} - z_{b,f} )\right ]\end{array}$

The term $\begin{array}{l}\mathcal{H}_f\end{array}$ is treated explicitly and added to the momentum equation.

Minor Losses

Energy losses occur as flow moves through and past the bridge structure. A proportion of them are associated with phenomena resolvable by a typical HEC-RAS 2D grid resolution, and are implicitly handled by the 2D flow solution. This is provided the RAS 2D grid is constructed using a “good enough" resolution and a proper turbulence model is specified. These resolvable losses include contraction and expansion losses. Other energy losses, however, are dependent on boundary layer dynamics and 3D flow phenomena and are unresolvable with a 2D depth-averaged model using a typical grid resolution. These include bridge pier, deck, and railing drag losses, plunging flow off the bridge deck, and vena contracta past parts of the bridge structure. To account for these subgrid scale losses, minor loss terms should be specified in the momentum equation.

The determination of minor losses in HEC-RAS's detailed bridge approach is based on the TUFLOW 2D bridge treatment, which itself is based on a joint study performed by TUFLOW and the Queensland Department of Transportation. In this study, 3D CFD modeling was performed to simulate flow past submerged bridges, and the results were used to parameterize the TUFLOW 2D model. Total minor losses for the bridge were calculated in a depth-dependent way by splitting the flow vertically into several layers and accounting for the relevant minor loss factors within each layer.

Minor losses are constant in the "pier layer" under low-flow conditions when the bridge deck is not wetted. As water levels rise, the minor losses increase linearly in the "super-structure layer" of the curve. In the TUFLOW study, minor losses were found to peak when the water level was approximately 0.6 times the bridge thickness above the top of the bridge deck. As the bridge becomes more inundated, minor losses decrease inversely to the depth.

The TUFLOW report characterized the peak minor loss coefficient as a function of the ratio of the bridge deck thickness to the height of the bridge deck above the bed.

Bridge deck height to thickness ratio	Minor Loss Value
2	0.42
4	0.28
6	0.2

The minor loss parameterization accounts for the total effect of flow through the bridge structure. Typically, a minor loss coefficient is used to calculate a head loss related to the local velocity head

$\begin{array}{l}\displaystyle H_L = K_L \, \frac{V^2}{2g}\end{array}$

where

$\begin{array}{l}H_L\end{array}$ is the energy head loss, and
$\begin{array}{l}K_L\end{array}$ is the minor loss coefficient.

For 2D bridges in HEC-RAS, these minor losses must be distributed across multiple faces, both above and below the bridge. This is accomplished by dividing the minor loss coefficient by the number of rows of computational faces perpendicular to flow, and within the bridge. The scaled loss coefficient is multiplied by the average velocity below the bridge deck for faces below the deck, or the average velocity above the bridge deck (for faces above). The energy loss is incorporated into the momentum equation as a force term ( $\begin{array}{l}F_{ML}\end{array}$ ) and applied at each face as:

$\begin{array}{l}\displaystyle \dfrac{F_{ML}}{\rho A} = g \dfrac{H_L}{2 \Delta x}\end{array}$

where

$\begin{array}{l}\Delta x\end{array}$ is the cell spacing, and the distance over which the minor losses are applied,
$\begin{array}{l}\rho\end{array}$ is the density of water, and
$\begin{array}{l}A\end{array}$ is the face cross-sectional area.