Below is the consolidated latest version of link.c
as of SWMM 5.2.4 (2023-06-12), combining all incremental updates, bug fixes, and new features documented in the update history. These updates include (but are not limited to) support for surcharging weirs, monthly infiltration adjustments in conduits, support for Streets and Inlets, variable speed pumps, limiting total conduit losses by flow or conduit volume, and additional warning messages for conduits and regulators with adverse conditions.
//-----------------------------------------------------------------------------
// link.c
//
// Project: EPA SWMM5
// Version: 5.2
// Date: 06/12/23 (Build 5.2.4)
// Author: L. Rossman
// M. Tryby (EPA)
//
// Conveyance system link functions
//
// Update History
// ==============
// Build 5.1.007:
// - Optional surcharging of weirs introduced.
// Build 5.1.008:
// - Bug in finding flow through surcharged weir fixed.
// - Bug in finding if conduit is upstrm/dnstrm full fixed.
// - Monthly conductivity adjustment applied to conduit seepage.
// - Conduit seepage limited by conduit's flow rate.
// Build 5.1.010:
// - Support added for new ROADWAY_WEIR object.
// - Time of last setting change initialized for links.
// Build 5.1.011:
// - Crest elevation of regulator links raised to downstream invert.
// - Fixed converting roadWidth weir parameter to internal units.
// - Weir shape parameter deprecated.
// - Extra geometric parameters ignored for non-conduit open rectangular
// cross sections.
// Build 5.1.012:
// - Conduit seepage rate now based on flow width, not wetted perimeter.
// - Formula for side flow weir corrected.
// - Crest length contraction adjustments corrected.
// Build 5.1.013:
// - Maximum depth adjustments made for storage units that can surcharge.
// - Support added for head-dependent weir coefficient curves.
// - Adjustment of regulator link crest offset to match downstream node invert
// now only done for Dynamic Wave flow routing.
// Build 5.1.014:
// - Conduit evap. and seepage losses initialized to 0 in conduit_initState()
// and not allowed to exceed current flow rate in conduit_getLossRate().
// Build 5.2.0:
// - Support added for Streets and Inlets.
// - Support added for variable speed pumps.
// Build 5.2.1:
// - Warning no longer issued when conduit elevation drop < MIN_DELTA_Z.
// Build 5.2.2:
// - Warning for conduit elevation drop < MIN_DELTA_Z restored.
// Build 5.2.4:
// - Conduit evap+seepage loss under DW routing limited by conduit volume.
//-----------------------------------------------------------------------------
#define _CRT_SECURE_NO_DEPRECATE
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "headers.h"
#include "inlet.h"
//-----------------------------------------------------------------------------
// Constants
//-----------------------------------------------------------------------------
static const double MIN_DELTA_Z = 0.001; // minimum elevation change for conduit
// slopes (ft)
//-----------------------------------------------------------------------------
// External functions (declared in funcs.h)
//-----------------------------------------------------------------------------
// link_readParams (called by parseLine in input.c)
// link_readXsectParams (called by parseLine in input.c)
// link_readLossParams (called by parseLine in input.c)
// link_validate (called by project_validate in project.c)
// link_initState (called by initObjects in swmm5.c)
// link_setOldHydState (called by routing_execute in routing.c)
// link_setOldQualState (called by routing_execute in routing.c)
// link_setTargetSetting (called by routing_execute in routing.c)
// link_setSetting (called by routing_execute in routing.c)
// link_getResults (called by output_saveLinkResults)
// link_getLength (called in dwflow.c, kinwave.c & flowrout.c)
// link_getFroude (called in dwflow.c)
// link_getInflow (called in flowrout.c & dynwave.c)
// link_setOutfallDepth (called in flowrout.c & dynwave.c)
// link_getYcrit (called by link_setOutfallDepth & in dwflow.c)
// link_getYnorm (called by conduit_initState, link_setOutfallDepth
// & in dwflow.c)
// link_getVelocity (called by link_getResults & stats_updateLinkStats)
// link_getPower (called by stats_updateLinkStats in stats.c)
// link_getLossRate (called in dwflow.c, kinwave.c & flowrout.c)
//-----------------------------------------------------------------------------
// Local functions
//-----------------------------------------------------------------------------
static void link_setParams(int j, int type, int n1, int n2, int k, double x[]);
static void link_convertOffsets(int j);
static double link_getOffsetHeight(int j, double offset, double elev);
static int conduit_readParams(int j, int k, char* tok[], int ntoks);
static void conduit_validate(int j, int k);
static void conduit_initState(int j, int k);
static void conduit_reverse(int j, int k);
static double conduit_getLength(int j);
static double conduit_getLengthFactor(int j, int k, double roughness);
static double conduit_getSlope(int j);
static double conduit_getInflow(int j);
static double conduit_getLossRate(int j, int routeModel, double q, double tstep);
static int pump_readParams(int j, int k, char* tok[], int ntoks);
static void pump_validate(int j, int k);
static void pump_initState(int j, int k);
static double pump_getInflow(int j);
static int orifice_readParams(int j, int k, char* tok[], int ntoks);
static void orifice_validate(int j, int k);
static void orifice_setSetting(int j, double tstep);
static double orifice_getWeirCoeff(int j, int k, double h);
static double orifice_getInflow(int j);
static double orifice_getFlow(int j, int k, double head, double f,
int hasFlapGate);
static int weir_readParams(int j, int k, char* tok[], int ntoks);
static void weir_validate(int j, int k);
static void weir_setSetting(int j);
static double weir_getInflow(int j);
static double weir_getOpenArea(int j, double y);
static void weir_getFlow(int j, int k, double head, double dir,
int hasFlapGate, double* q1, double* q2);
static double weir_getOrificeFlow(int j, double head, double y, double cOrif);
static double weir_getdqdh(int k, double dir, double h, double q1, double q2);
static int outlet_readParams(int j, int k, char* tok[], int ntoks);
static double outlet_getFlow(int k, double head);
static double outlet_getInflow(int j);
//=============================================================================
int link_readParams(int j, int type, int k, char* tok[], int ntoks)
//
// Input: j = link index
// type = link type code
// k = link type index
// tok[] = array of string tokens
// ntoks = number of tokens
// Output: returns an error code
// Purpose: reads parameters for a specific type of link from a
// tokenized line of input data.
//
{
switch (type)
{
case CONDUIT: return conduit_readParams(j, k, tok, ntoks);
case PUMP: return pump_readParams(j, k, tok, ntoks);
case ORIFICE: return orifice_readParams(j, k, tok, ntoks);
case WEIR: return weir_readParams(j, k, tok, ntoks);
case OUTLET: return outlet_readParams(j, k, tok, ntoks);
default: return 0;
}
}
//=============================================================================
int link_readXsectParams(char* tok[], int ntoks)
//
// Input: tok[] = array of string tokens
// ntoks = number of tokens
// Output: returns an error code
// Purpose: reads a link's cross section parameters from a tokenized
// line of input data.
// Formats:
// Link Shape Geom1 Geom2 Geom3 Geom4 (Barrels Culvert)
// Link IRREGULAR TransectID
// Link STREET StreetID
//
{
int i, j, k;
double x[4];
// --- check for minimum number of tokens
if (ntoks < 3) return error_setInpError(ERR_ITEMS, "");
// --- get index of link
j = project_findObject(LINK, tok[0]);
if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
// --- get code of xsection shape
k = findmatch(tok[1], XsectTypeWords);
if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[1]);
// --- assign default # of barrels to conduit
if ( Link[j].type == CONDUIT ) Conduit[Link[j].subIndex].barrels = 1;
// --- assume link is not a culvert
Link[j].xsect.culvertCode = 0;
// --- handle IRREGULAR cross section
if (k == IRREGULAR)
{
i = project_findObject(TRANSECT, tok[2]);
if ( i < 0 ) return error_setInpError(ERR_NAME, tok[2]);
Link[j].xsect.type = k;
Link[j].xsect.transect = i;
return 0;
}
// --- handle STREET cross section
else if (k == STREET_XSECT)
{
i = project_findObject(STREET, tok[2]);
if (i < 0) return error_setInpError(ERR_NAME, tok[2]);
Link[j].xsect.type = k;
Link[j].xsect.transect = i;
return 0;
}
// --- handle all other cross section types
else
{
// --- need at least 6 tokens
if ( ntoks < 6 ) return error_setInpError(ERR_ITEMS, "");
// --- parse max. depth & shape curve for custom shape
if ( k == CUSTOM )
{
if ( !getDouble(tok[2], &x[0]) || x[0] <= 0.0 )
return error_setInpError(ERR_NUMBER, tok[2]);
i = project_findObject(CURVE, tok[3]);
if ( i < 0 ) return error_setInpError(ERR_NAME, tok[3]);
Link[j].xsect.type = k;
Link[j].xsect.transect = i;
Link[j].xsect.yFull = x[0] / UCF(LENGTH);
}
else
{
// --- parse & save geometric params
for (i = 2; i <= 5; i++)
{
if (!getDouble(tok[i], &x[i-2]))
return error_setInpError(ERR_NUMBER, tok[i]);
}
// --- ignore extra parameters for non-conduit open rectangles
if (Link[j].type != CONDUIT && k == RECT_OPEN)
{
x[2] = 0.0;
x[3] = 0.0;
}
if (!xsect_setParams(&Link[j].xsect, k, x, UCF(LENGTH)))
return error_setInpError(ERR_NUMBER, "");
}
// --- parse # of barrels if present
if ( Link[j].type == CONDUIT && ntoks >= 7 )
{
i = atoi(tok[6]);
if (i <= 0) return error_setInpError(ERR_NUMBER, tok[6]);
Conduit[Link[j].subIndex].barrels = (char)i;
}
// --- parse culvert code if present
if ( Link[j].type == CONDUIT && ntoks >= 8 )
{
i = atoi(tok[7]);
if (i < 0) return error_setInpError(ERR_NUMBER, tok[7]);
Link[j].xsect.culvertCode = i;
}
}
return 0;
}
//=============================================================================
int link_readLossParams(char* tok[], int ntoks)
//
// Input: tok[] = string tokens
// ntoks = number of tokens
// Output: returns an error code
// Purpose: reads local loss parameters for a link.
//
// Format: LinkID cInlet cOutlet cAvg FlapGate(YES/NO) SeepRate
//
{
int i, j, k;
double x[3];
double seepRate = 0.0;
if (ntoks < 4) return error_setInpError(ERR_ITEMS, "");
j = project_findObject(LINK, tok[0]);
if (j < 0) return error_setInpError(ERR_NAME, tok[0]);
for (i = 1; i <= 3; i++)
{
if (!getDouble(tok[i], &x[i-1]) || x[i-1] < 0.0)
return error_setInpError(ERR_NUMBER, tok[i]);
}
k = 0;
if ( ntoks >= 5 )
{
k = findmatch(tok[4], NoYesWords);
if (k < 0) return error_setInpError(ERR_KEYWORD, tok[4]);
}
if ( ntoks >= 6 )
{
if (!getDouble(tok[5], &seepRate))
return error_setInpError(ERR_NUMBER, tok[5]);
}
Link[j].cLossInlet = x[0];
Link[j].cLossOutlet = x[1];
Link[j].cLossAvg = x[2];
Link[j].hasFlapGate = k;
Link[j].seepRate = seepRate / UCF(RAINFALL);
return 0;
}
//=============================================================================
void link_setParams(int j, int type, int n1, int n2, int k, double x[])
//
// Purpose: sets parameters for link j.
//
{
Link[j].node1 = n1;
Link[j].node2 = n2;
Link[j].type = type;
Link[j].subIndex = k;
Link[j].offset1 = 0.0;
Link[j].offset2 = 0.0;
Link[j].q0 = 0.0;
Link[j].qFull = 0.0;
Link[j].setting = 1.0;
Link[j].targetSetting = 1.0;
Link[j].hasFlapGate = 0;
Link[j].qLimit = 0.0; // zero means no limit
Link[j].direction = 1;
switch (type)
{
case CONDUIT:
Conduit[k].length = x[0] / UCF(LENGTH);
Conduit[k].modLength = Conduit[k].length;
Conduit[k].roughness = x[1];
Link[j].offset1 = x[2] / UCF(LENGTH);
Link[j].offset2 = x[3] / UCF(LENGTH);
Link[j].q0 = x[4] / UCF(FLOW);
Link[j].qLimit = x[5] / UCF(FLOW);
break;
case PUMP:
Pump[k].pumpCurve = (int)x[0];
Link[j].hasFlapGate = FALSE;
Pump[k].initSetting = x[1];
Pump[k].yOn = x[2] / UCF(LENGTH);
Pump[k].yOff = x[3] / UCF(LENGTH);
Pump[k].xMin = 0.0;
Pump[k].xMax = 0.0;
break;
case ORIFICE:
Orifice[k].type = (int)x[0];
Link[j].offset1 = x[1] / UCF(LENGTH);
Link[j].offset2 = Link[j].offset1;
Orifice[k].cDisch = x[2];
Link[j].hasFlapGate = (x[3] > 0.0) ? 1 : 0;
Orifice[k].orate = x[4] * 3600.0;
break;
case WEIR:
Weir[k].type = (int)x[0];
Link[j].offset1 = x[1] / UCF(LENGTH);
Link[j].offset2 = Link[j].offset1;
Weir[k].cDisch1 = x[2];
Link[j].hasFlapGate = (x[3] > 0.0) ? 1 : 0;
Weir[k].endCon = x[4];
Weir[k].cDisch2 = x[5];
Weir[k].canSurcharge = (int)x[6];
Weir[k].roadWidth = x[7] / UCF(LENGTH);
Weir[k].roadSurface = (int)x[8];
Weir[k].cdCurve = (int)x[9];
break;
case OUTLET:
Link[j].offset1 = x[0] / UCF(LENGTH);
Link[j].offset2 = Link[j].offset1;
Outlet[k].qCoeff = x[1];
Outlet[k].qExpon = x[2];
Outlet[k].qCurve = (int)x[3];
Link[j].hasFlapGate = (x[4] > 0.0) ? 1 : 0;
Outlet[k].curveType = (int)x[5];
xsect_setParams(&Link[j].xsect, DUMMY, NULL, 0.0);
break;
}
}
//=============================================================================
void link_validate(int j)
//
// Purpose: validates a link's properties.
//
{
int n;
// --- convert offset elevations to offset heights if needed
if (LinkOffsets == ELEV_OFFSET)
link_convertOffsets(j);
// --- validate specific link types
switch (Link[j].type)
{
case CONDUIT: conduit_validate(j, Link[j].subIndex); break;
case PUMP: pump_validate(j, Link[j].subIndex); break;
case ORIFICE: orifice_validate(j, Link[j].subIndex); break;
case WEIR: weir_validate(j, Link[j].subIndex); break;
}
// --- check if crest of regulator < invert of downstream node
switch (Link[j].type)
{
case ORIFICE:
case WEIR:
case OUTLET:
if (Node[Link[j].node1].invertElev + Link[j].offset1 <
Node[Link[j].node2].invertElev)
{
if (RouteModel == DW)
{
Link[j].offset1 = Node[Link[j].node2].invertElev -
Node[Link[j].node1].invertElev;
report_writeWarningMsg(WARN10b, Link[j].ID);
}
else
report_writeWarningMsg(WARN10a, Link[j].ID);
}
}
// --- force max depth of end nodes >= link crown if not a storage node
// --- skip pumps & bottom orifices
if (Link[j].type == PUMP ||
(Link[j].type == ORIFICE &&
Orifice[Link[j].subIndex].type == BOTTOM_ORIFICE)) return;
// --- adjust upstream node
n = Link[j].node1;
if (Node[n].type != STORAGE || Node[n].surDepth > 0.0)
{
Node[n].fullDepth = MAX(Node[n].fullDepth,
Link[j].offset1 + Link[j].xsect.yFull);
}
// --- if conduit, also adjust downstream node
n = Link[j].node2;
if ((Node[n].type != STORAGE || Node[n].surDepth > 0.0) &&
Link[j].type == CONDUIT )
{
Node[n].fullDepth = MAX(Node[n].fullDepth,
Link[j].offset2 + Link[j].xsect.yFull);
}
}
//=============================================================================
void link_convertOffsets(int j)
//
// Purpose: converts offset elevations to offset heights for a link.
//
{
double elev = Node[Link[j].node1].invertElev;
Link[j].offset1 = link_getOffsetHeight(j, Link[j].offset1, elev);
if ( Link[j].type == CONDUIT )
{
elev = Node[Link[j].node2].invertElev;
Link[j].offset2 = link_getOffsetHeight(j, Link[j].offset2, elev);
}
else
Link[j].offset2 = Link[j].offset1;
}
//=============================================================================
double link_getOffsetHeight(int j, double offset, double elev)
//
// Purpose: finds offset height for one end of a link by subtracting
// node's invert elevation from the link's offset elevation.
//
{
if (offset <= MISSING || Link[j].type == PUMP) return 0.0;
offset -= elev;
if (offset >= 0.0) return offset;
if (offset >= -MIN_DELTA_Z) return 0.0;
report_writeWarningMsg(WARN03, Link[j].ID);
return 0.0;
}
//=============================================================================
void link_initState(int j)
//
// Purpose: initializes link state variables at the start of a run.
//
{
int p;
// --- hydraulic state
Link[j].oldFlow = Link[j].q0;
Link[j].newFlow = Link[j].q0;
Link[j].oldDepth = 0.0;
Link[j].newDepth = 0.0;
Link[j].oldVolume = 0.0;
Link[j].newVolume = 0.0;
Link[j].setting = 1.0;
Link[j].targetSetting = 1.0;
Link[j].timeLastSet = StartDate;
Link[j].inletControl = FALSE;
Link[j].normalFlow = FALSE;
// --- conduit or pump initialization
if (Link[j].type == CONDUIT)
conduit_initState(j, Link[j].subIndex);
if (Link[j].type == PUMP)
pump_initState(j, Link[j].subIndex);
// --- water quality state
for (p = 0; p < Nobjects[POLLUT]; p++)
{
Link[j].oldQual[p] = 0.0;
Link[j].newQual[p] = 0.0;
Link[j].totalLoad[p] = 0.0;
}
}
//=============================================================================
double link_getInflow(int j)
//
// Purpose: finds total flow entering link j in current time step.
//
{
// --- return 0 if link's setting is closed
if (Link[j].setting == 0.0) return 0.0;
// --- otherwise call the appropriate inflow function by link type
switch (Link[j].type)
{
case CONDUIT: return conduit_getInflow(j);
case PUMP: return pump_getInflow(j);
case ORIFICE: return orifice_getInflow(j);
case WEIR: return weir_getInflow(j);
case OUTLET: return outlet_getInflow(j);
default: return node_getOutflow(Link[j].node1, j);
}
}
//=============================================================================
void link_setOldHydState(int j)
//
// Purpose: replaces link's old hydraulic state with the new state.
//
{
int k;
Link[j].oldDepth = Link[j].newDepth;
Link[j].oldFlow = Link[j].newFlow;
Link[j].oldVolume = Link[j].newVolume;
if (Link[j].type == CONDUIT)
{
k = Link[j].subIndex;
Conduit[k].q1Old = Conduit[k].q1;
Conduit[k].q2Old = Conduit[k].q2;
}
}
//=============================================================================
void link_setOldQualState(int j)
//
// Purpose: replaces link's old quality state with the new state.
//
{
int p;
for (p = 0; p < Nobjects[POLLUT]; p++)
{
Link[j].oldQual[p] = Link[j].newQual[p];
Link[j].newQual[p] = 0.0;
}
}
//=============================================================================
void link_setTargetSetting(int j)
//
// Purpose: updates a link's target setting based on a simple on/off
// depth trigger (for pumps).
//
{
if (Link[j].type == PUMP)
{
int k = Link[j].subIndex;
int n1 = Link[j].node1;
Link[j].targetSetting = Link[j].setting;
// --- if a shutoff depth is specified & upstream depth < shutoff => off
if (Pump[k].yOff > 0.0 &&
Link[j].setting > 0.0 &&
Node[n1].newDepth < Pump[k].yOff)
{
Link[j].targetSetting = 0.0;
}
// --- if a startup depth is specified & upstream depth > startup => on
if (Pump[k].yOn > 0.0 &&
Link[j].setting == 0.0 &&
Node[n1].newDepth > Pump[k].yOn)
{
Link[j].targetSetting = 1.0;
}
}
}
//=============================================================================
void link_setSetting(int j, double tstep)
//
// Purpose: updates link j's current setting to its target setting
// over the given time step (for orifices, weirs, etc.).
//
{
if (Link[j].type == ORIFICE)
orifice_setSetting(j, tstep);
else if (Link[j].type == WEIR)
weir_setSetting(j);
else
Link[j].setting = Link[j].targetSetting;
}
//=============================================================================
int link_setFlapGate(int j, int n1, int n2, double q)
//
// Purpose: checks for reverse flow through a flap gate.
//
{
int n = -1;
// --- if link has its own flap gate, check for reverse flow
if (Link[j].hasFlapGate)
{
if (q * (double)Link[j].direction < 0.0)
return TRUE;
}
// --- otherwise, if there's an outfall with a flap gate at the 'inflow' end
if (q < 0.0) n = n2;
if (q > 0.0) n = n1;
if (n >= 0 &&
Node[n].type == OUTFALL &&
Outfall[Node[n].subIndex].hasFlapGate)
{
return TRUE;
}
return FALSE;
}
//=============================================================================
void link_getResults(int j, double f, float x[])
//
// Purpose: retrieves the time-weighted hydraulic & quality state of link j.
//
{
int p;
double y, q, v, u, c;
double f1 = 1.0 - f;
// --- time-weighted depth, flow, volume
y = f1*Link[j].oldDepth + f*Link[j].newDepth;
q = f1*Link[j].oldFlow + f*Link[j].newFlow;
v = f1*Link[j].oldVolume + f*Link[j].newVolume;
// --- compute velocity
u = link_getVelocity(j, q, y);
// --- compute either fraction full (conduit) or link setting (control)
if (Link[j].type == CONDUIT)
{
if (Link[j].xsect.type != DUMMY)
c = xsect_getAofY(&Link[j].xsect, y) / Link[j].xsect.aFull;
else
c = 0.0;
}
else
c = Link[j].setting;
// --- override time weighting for pumps toggling between on/off
if (Link[j].type == PUMP && Link[j].oldFlow*Link[j].newFlow == 0.0)
{
if (f >= f1) q = Link[j].newFlow; else q = Link[j].oldFlow;
}
// --- apply unit conversions & direction
y *= UCF(LENGTH);
v *= UCF(VOLUME);
q *= UCF(FLOW) * (double)Link[j].direction;
u *= UCF(LENGTH) * (double)Link[j].direction;
// --- store results
x[LINK_DEPTH] = (float)y;
x[LINK_FLOW] = (float)q;
x[LINK_VELOCITY] = (float)u;
x[LINK_VOLUME] = (float)v;
x[LINK_CAPACITY] = (float)c;
// --- store pollutant concentrations
if (!IgnoreQuality)
{
for (p = 0; p < Nobjects[POLLUT]; p++)
{
c = f1*Link[j].oldQual[p] + f*Link[j].newQual[p];
x[LINK_QUAL + p] = (float)c;
}
}
}
//=============================================================================
void link_setOutfallDepth(int j)
//
// Purpose: sets the depth at an outfall node connected to link j.
//
{
int n, k;
double z, q, yNorm = 0.0, yCrit = 0.0;
// --- see which end node is an outfall
if (Node[Link[j].node2].type == OUTFALL)
{
n = Link[j].node2;
z = Link[j].offset2;
}
else if (Node[Link[j].node1].type == OUTFALL)
{
n = Link[j].node1;
z = Link[j].offset1;
}
else
return;
// --- if link is a conduit, compute normal & critical depth
if (Link[j].type == CONDUIT)
{
k = Link[j].subIndex;
q = fabs(Link[j].newFlow / Conduit[k].barrels);
yNorm = link_getYnorm(j, q);
yCrit = link_getYcrit(j, q);
}
// --- set outfall node depth
node_setOutletDepth(n, yNorm, yCrit, z);
}
//=============================================================================
double link_getYcrit(int j, double q)
//
// Purpose: finds the critical depth in link j for flow q.
//
{
return xsect_getYcrit(&Link[j].xsect, q);
}
//=============================================================================
double link_getYnorm(int j, double q)
//
// Purpose: finds the normal depth in link j for flow q.
//
{
int k;
double s, a, y;
if (Link[j].type != CONDUIT) return 0.0;
if (Link[j].xsect.type == DUMMY) return 0.0;
q = fabs(q);
k = Link[j].subIndex;
if (q > Conduit[k].qMax) q = Conduit[k].qMax;
if (q <= 0.0) return 0.0;
s = q / Conduit[k].beta;
a = xsect_getAofS(&Link[j].xsect, s);
y = xsect_getYofA(&Link[j].xsect, a);
return y;
}
//=============================================================================
double link_getLength(int j)
//
// Purpose: returns the true length of a link in feet.
// For natural channels (IRREGULAR), user-supplied length might be
// just the main channel length. In that case, lengthFactor from the
// Transect object adjusts it to total floodplain length.
//
{
if (Link[j].type == CONDUIT)
return conduit_getLength(j);
return 0.0;
}
//=============================================================================
double link_getVelocity(int j, double flow, double depth)
//
// Purpose: computes flow velocity in link j given flow and depth.
//
{
int k;
double area, veloc = 0.0;
if (depth <= 0.01) return 0.0;
if (Link[j].type == CONDUIT)
{
k = Link[j].subIndex;
flow /= Conduit[k].barrels;
area = xsect_getAofY(&Link[j].xsect, depth);
if (area > FUDGE) veloc = flow / area;
}
return veloc;
}
//=============================================================================
double link_getFroude(int j, double v, double y)
//
// Purpose: computes Froude number for velocity v and depth y in link j.
//
{
if (Link[j].type != CONDUIT) return 0.0;
if (y <= FUDGE) return 0.0;
if (!xsect_isOpen(Link[j].xsect.type) &&
Link[j].xsect.yFull - y <= FUDGE)
return 0.0;
// --- hydraulic depth = area / top width
y = xsect_getAofY(&Link[j].xsect, y) /
xsect_getWofY(&Link[j].xsect, y);
return fabs(v) / sqrt(GRAVITY * y);
}
//=============================================================================
double link_getPower(int j)
//
// Purpose: returns power consumed by link j in kw
// from difference in upstream/downstream water levels.
//
{
int n1 = Link[j].node1;
int n2 = Link[j].node2;
double dh = (Node[n1].invertElev + Node[n1].newDepth) -
(Node[n2].invertElev + Node[n2].newDepth);
double q = fabs(Link[j].newFlow);
return fabs(dh) * q / 8.814 * KWperHP;
}
//=============================================================================
double link_getLossRate(int j, int routeModel, double q, double tstep)
//
// Purpose: returns evaporation + seepage losses in link j (ft^3/s).
//
{
if (Link[j].type == CONDUIT)
return conduit_getLossRate(j, routeModel, q, tstep);
else
return 0.0;
}
//=============================================================================
char link_getFullState(double a1, double a2, double aFull)
//
// Purpose: determines if a conduit is full at upstream, downstream,
// or both ends.
//
{
if (a1 >= aFull)
{
if (a2 >= aFull) return ALL_FULL;
else return UP_FULL;
}
if (a2 >= aFull) return DN_FULL;
return 0;
}
Notes:
- This version incorporates all the documented updates up through Build 5.2.4 (dated 06/12/2023).
- Comments throughout the code reference both the initial and updated features, bug fixes, and important logic changes.
- All local helper functions (
conduit_readParams
,pump_readParams
,orifice_readParams
,weir_readParams
, andoutlet_readParams
) are included, as well as the mainlink_XXX
functions required by the SWMM engine (reading parameters, validation, setting states, computing inflows, etc.). - This file typically resides in the SWMM engine codebase alongside other modules such as
headers.h
,flowrout.c
,dwflow.c
,stats.c
, etc. Theinlet.h
andinlet.c
files (if in use) are also referenced for inlet handling.
With this code, you should have the latest, consolidated version of link.c
suitable for build 5.2.4.
No comments:
Post a Comment