Summary of dwflow.c in SWMM5

Below is a detailed overview of dwflow.c, which is the portion of SWMM5 that handles dynamic wave hydraulic computations inside a single conduit. It specifically focuses on how SWMM uses the momentum equation (St. Venant) with various boundary conditions (subcritical, supercritical, surcharged) to update the conduit’s flow each time step. This file is tightly integrated with the rest of SWMM’s hydraulic engine (found in dynwave.c, flowrout.c, etc.), but here we focus on the finite‐difference solution for one conduit under dynamic wave routing.

---

1. Purpose of dwflow.c

• Dynamic wave flow is SWMM’s most detailed hydraulic mode.
• dwflow.c implements the specialized code to:
  1. Solve the momentum equation in a conduit (pipe or channel) under dynamic wave conditions.
  2. Perform one update of the conduit’s flow, depth, and flow classification at each iteration/time step.
  3. Account for:
     • Friction slope (Manning or Darcy-Weisbach for force mains).
     • Local (entrance/exit) losses.
     • Inertial (acceleration) terms.
     • Preissmann slot approach for surcharged pipes.
     • Evaporation/infiltration losses in the conduit barrel.

The dynamic wave engine calls these routines for each conduit, iterating until a system-wide flow solution converges.

---

2. Key Functions

2.1 dwflow_findConduitFlow(...)

void dwflow_findConduitFlow(int j, int steps, double omega, double dt)

• Inputs:
  • j: the link (conduit) index.
  • steps: how many internal iteration steps have already been taken for the current time step.
  • omega: an under‐relaxation factor (0–1).
  • dt: the current routing time step in seconds.
• Outputs:
  • It updates Link[j].newFlow, Link[j].newDepth, Link[j].newVolume, etc.
• Algorithm:
  1. Gathers upstream/downstream node heads $`h1`, `h2`$ and conduit geometry.
  2. Calculates “effective” flow depths (y1, y2) at each end, ensuring they are $\ge$ FUDGE to avoid zero area.
  3. Finds cross-sectional areas, hydraulic radii, midpoint area, checks if the conduit is full.
  4. Applies the finite difference form of the momentum equation: $$q = \frac{q_{\mathrm{old}} + [\text{inertial, slope, friction, local loss, seepage}]}{\text{denominator}}$$
     • Terms (dq1 through dq6) represent friction slope, gravity slope, inertial difference, local losses, etc.
     • Some inertial damping might apply based on the Froude number or user‐chosen inertial damping option.
  5. Possibly limit the flow to normal flow or inlet‐controlled culvert flow if that is smaller (depending on user’s normal flow or Froude criteria).
  6. Applies under‐relaxation: q_new = (1 - omega)*q_last + omega*q.
  7. Updates various variables: conduit’s mid-area, the new flow depth, volume, flow, etc.

In short, dwflow_findConduitFlow is the core routine that:

• Estimates the new flow in a conduit,
• Balances the momentum equation’s terms,
• Enforces flow or velocity constraints (e.g., normal flow limitation, infiltration, culvert, flap gates).

---

2.2 findSurfArea(...)

static void findSurfArea(int j, double q, double length,
    double* h1, double* h2, double* y1, double* y2)

• Purpose: Divides the conduit’s surface area between its upstream and downstream nodes (for stability in dynamic wave).
• How:
  1. Based on the conduit’s flow classification (subcritical, supercritical, critical, or partially dry), it calculates an “effective top width” and how much portion belongs to each node’s “surface area.”
  2. Adjusts h1, h2, y1, y2 if a node is at critical or normal depth.
  3. Helps the dynamic wave solver by distributing the conduit’s free surface to the node surcharges.

---

2.3 getFlowClass(...)

static int getFlowClass(
    int link, double q, double h1, double h2,
    double y1, double y2, double* yC, double* yN, double* fasnh)

• Classifies the conduit’s flow regime as:
  • SUBCRITICAL, SUPCRITICAL,
  • UP_CRITICAL, DN_CRITICAL,
  • UP_DRY, DN_DRY,
  • or DRY.
• It checks depths and heads at each end, compares them to normal and critical depths (found via link_getYnorm, link_getYcrit).
• fasnh is the fraction between normal & critical depth used to distribute free surface area.

---

2.4 findLocalLosses(...)

static double findLocalLosses(int link, double a1, double a2, double aMid, double q)

• If the user sets local minor loss coefficients (cLossInlet, cLossOutlet, cLossAvg), then for each end or midpoint, we add a term $\sim \text{coefficient} \times q/a$.
• Summed up as part of the momentum equation.

---

2.5 getSlotWidth(...), getWidth(...), getArea(...), getHydRad(...)

• Various helpers to derive:
  • The Preissmann slot top width if surcharged (SLOT method),
  • The top width of the cross‐section at a given depth,
  • The area of flow at a given depth,
  • The hydraulic radius.
• These rely heavily on cross‐section geometry (TXsect) and the presence (or not) of the Preissmann slot.

---

2.6 checkNormalFlow(...)

static double checkNormalFlow(int j, double q, double y1, double y2, double a1, double r1)

• Goal: Possibly override the dynamic wave flow with normal flow if:
  1. The user has set “normal flow limit” (SLOPE, FROUDE, or BOTH).
  2. A slope or Froude condition is met, meaning the flow is “supercritical or steep enough” that normal flow may apply.
• If normal flow is less than the dynamic wave flow, we reduce flow to that normal capacity.
• We do: $$q_{\mathrm{norm}} = \beta \times a_1 \times (r_1)^{2/3}$$ (where $\beta$ is basically $1/n$ or a force main friction factor).
• Then, if qNorm < q, we set q = qNorm.

---

3. Key Data Fields & Variables

• Conduit[k].q1 and Conduit[k].q2: store new flows at the start/end or iteration.
• Conduit[k].barrels: number of parallel barrels. Flow in SWMM is computed per barrel, then multiplied.
• Link[j].oldFlow, Link[j].newFlow: flows in cfs at the end of the last time step vs. current iteration result.
• Link[j].flowClass: the flow classification code (dry, subcritical, supercritical, etc.).
• Link[j].dqdh: partial derivative $\partial Q / \partial H$ used by the dynamic wave solver’s global matrix.
• Link[j].froude: Froude number to detect supercritical conditions.

---

4. Typical Flow of Computation

1. At each iteration in the dynamic wave solver (in dynwave_execute(...)), SWMM calls: \text{dwflow_findConduitFlow}(j, \dots) \quad \text{for each conduit } j.
2. The function:
   1. Reads previous iteration flow (qLast) and node depths (h1, h2).
   2. Determines partial or full flow, calculates cross‐section areas.
   3. Summons the momentum equation in finite difference form: $$q = \frac{q_\mathrm{old} - \text{(gravity slope term)} + \text{(inertial term)} + \dots}{1 + \text{(friction)} + \text{(local losses)}}$$
   4. Possibly modifies q for:
      • Culvert inlet control,
      • Normal flow limits,
      • Flap gates disallowing reverse flow,
      • User flow limit.
   5. Under‐relaxes with omega.
3. The updated flow is stored in Conduit[k].q1, and Link[j].newFlow is q * barrels.
4. Once each conduit is updated, the dynamic wave solver checks for convergence and possibly repeats more iteration passes.

---

5. Key Takeaways

• dwflow.c is specialized code for one conduit’s momentum equation solution under dynamic wave.
• The code ensures we capture:
  • Full or partial pipe, plus the possibility of Preissmann slot for surcharging.
  • Inertial terms that can be damped (none, partial, or full).
  • Local and friction losses.
  • Flow constraints (inlet control, normal flow, flap gates).
• Ultimately, the rest of SWMM orchestrates calls to dwflow_findConduitFlow() within a global loop, each iteration refining flows in all conduits until the entire network converges at a time step—thus giving a robust unsteady flow simulation.