/*
|
******************************************************************************
|
Project: OWA EPANET
|
Version: 2.2
|
Module: qualreact.c
|
Description: computes water quality reactions within pipes and tanks
|
Authors: see AUTHORS
|
Copyright: see AUTHORS
|
License: see LICENSE
|
Last Updated: 05/15/2019
|
******************************************************************************
|
*/
|
|
#include <stdlib.h>
|
#include <stdio.h>
|
#include <math.h>
|
#include "types.h"
|
|
// Exported functions
|
char setreactflag(Project *);
|
double getucf(double);
|
void ratecoeffs(Project *);
|
void reactpipes(Project *, long);
|
void reacttanks(Project *, long);
|
double mixtank(Project *, int, double, double ,double);
|
|
// Imported functions
|
extern void addseg(Project *, int, double, double);
|
extern void reversesegs(Project *, int);
|
|
// Local functions
|
static double piperate(Project *, int);
|
static double pipereact(Project *, int, double, double, long);
|
static double tankreact(Project *, double, double, double, long);
|
static double bulkrate(Project *, double, double, double);
|
static double wallrate(Project *, double, double, double, double);
|
|
static void tankmix1(Project *, int, double, double, double);
|
static void tankmix2(Project *, int, double, double, double);
|
static void tankmix3(Project *, int, double, double, double);
|
static void tankmix4(Project *, int, double, double, double);
|
|
|
char setreactflag(Project *pr)
|
/*
|
**-----------------------------------------------------------
|
** Input: none
|
** Output: returns 1 for reactive WQ constituent, 0 otherwise
|
** Purpose: checks if reactive chemical being simulated
|
**-----------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
int i;
|
|
if (pr->quality.Qualflag == TRACE) return 0;
|
else if (pr->quality.Qualflag == AGE) return 1;
|
else
|
{
|
for (i = 1; i <= net->Nlinks; i++)
|
{
|
if (net->Link[i].Type <= PIPE)
|
{
|
if (net->Link[i].Kb != 0.0 || net->Link[i].Kw != 0.0) return 1;
|
}
|
}
|
for (i = 1; i <= net->Ntanks; i++)
|
{
|
if (net->Tank[i].Kb != 0.0) return 1;
|
}
|
}
|
return 0;
|
}
|
|
|
double getucf(double order)
|
/*
|
**--------------------------------------------------------------
|
** Input: order = bulk reaction order
|
** Output: returns a unit conversion factor
|
** Purpose: converts bulk reaction rates from per Liter to
|
** per FT3 basis
|
**--------------------------------------------------------------
|
*/
|
{
|
if (order < 0.0) order = 0.0;
|
if (order == 1.0) return (1.0);
|
else return (1. / pow(LperFT3, (order - 1.0)));
|
}
|
|
|
void ratecoeffs(Project *pr)
|
/*
|
**--------------------------------------------------------------
|
** Input: none
|
** Output: none
|
** Purpose: determines wall reaction coeff. for each pipe
|
**--------------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Quality *qual = &pr->quality;
|
|
int k;
|
double kw;
|
|
for (k = 1; k <= net->Nlinks; k++)
|
{
|
kw = net->Link[k].Kw;
|
if (kw != 0.0) kw = piperate(pr, k);
|
net->Link[k].Rc = kw;
|
qual->PipeRateCoeff[k] = 0.0;
|
}
|
}
|
|
|
void reactpipes(Project *pr, long dt)
|
/*
|
**--------------------------------------------------------------
|
** Input: dt = time step
|
** Output: none
|
** Purpose: reacts water within each pipe over a time step.
|
**--------------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Quality *qual = &pr->quality;
|
|
int k;
|
Pseg seg;
|
double cseg, rsum, vsum;
|
|
// Examine each link in network
|
for (k = 1; k <= net->Nlinks; k++)
|
{
|
// Skip non-pipe links (pumps & valves)
|
if (net->Link[k].Type != PIPE) continue;
|
rsum = 0.0;
|
vsum = 0.0;
|
|
// Examine each segment of the pipe
|
seg = qual->FirstSeg[k];
|
while (seg != NULL)
|
{
|
// React segment over time dt
|
cseg = seg->c;
|
seg->c = pipereact(pr, k, seg->c, seg->v, dt);
|
|
// Update reaction component of mass balance
|
qual->MassBalance.reacted += (cseg - seg->c) * seg->v;
|
|
// Accumulate volume-weighted reaction rate
|
if (qual->Qualflag == CHEM)
|
{
|
rsum += fabs(seg->c - cseg) * seg->v;
|
vsum += seg->v;
|
}
|
seg = seg->prev;
|
}
|
|
// Normalize volume-weighted reaction rate
|
if (vsum > 0.0) qual->PipeRateCoeff[k] = rsum / vsum / dt * SECperDAY;
|
else qual->PipeRateCoeff[k] = 0.0;
|
}
|
}
|
|
|
void reacttanks(Project *pr, long dt)
|
/*
|
**--------------------------------------------------------------
|
** Input: dt = time step
|
** Output: none
|
** Purpose: reacts water within each tank over a time step.
|
**--------------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Quality *qual = &pr->quality;
|
|
int i, k;
|
double c;
|
Pseg seg;
|
Stank *tank;
|
|
// Examine each tank in network
|
for (i = 1; i <= net->Ntanks; i++)
|
{
|
// Skip reservoirs
|
tank = &net->Tank[i];
|
if (tank->A == 0.0) continue;
|
|
// k is segment chain belonging to tank i
|
k = net->Nlinks + i;
|
|
// React each volume segment in the chain
|
seg = qual->FirstSeg[k];
|
while (seg != NULL)
|
{
|
c = seg->c;
|
seg->c = tankreact(pr, seg->c, seg->v, tank->Kb, dt);
|
qual->MassBalance.reacted += (c - seg->c) * seg->v;
|
seg = seg->prev;
|
}
|
}
|
}
|
|
|
double piperate(Project *pr, int k)
|
/*
|
**--------------------------------------------------------------
|
** Input: k = link index
|
** Output: returns reaction rate coeff. for 1st-order wall
|
** reactions or mass transfer rate coeff. for 0-order
|
** reactions
|
** Purpose: finds wall reaction rate coeffs.
|
**--------------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Hydraul *hyd = &pr->hydraul;
|
Quality *qual = &pr->quality;
|
|
double a, d, u, q, kf, kw, y, Re, Sh;
|
|
d = net->Link[k].Diam; // Pipe diameter, ft
|
|
// Ignore mass transfer if Schmidt No. is 0
|
if (qual->Sc == 0.0)
|
{
|
if (qual->WallOrder == 0.0) return BIG;
|
else return (net->Link[k].Kw * (4.0 / d) / pr->Ucf[ELEV]);
|
}
|
|
// Compute Reynolds No.
|
// Flow rate made consistent with how its saved to hydraulics file
|
q = (hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlow[k];
|
a = PI * d * d / 4.0; // pipe area
|
u = fabs(q) / a; // flow velocity
|
Re = u * d / hyd->Viscos; // Reynolds number
|
|
// Compute Sherwood No. for stagnant flow
|
// (mass transfer coeff. = Diffus./radius)
|
if (Re < 1.0) Sh = 2.0;
|
|
// Compute Sherwood No. for turbulent flow using the Notter-Sleicher formula.
|
else if (Re >= 2300.0) Sh = 0.0149 * pow(Re, 0.88) * pow(qual->Sc, 0.333);
|
|
// Compute Sherwood No. for laminar flow using Graetz solution formula.
|
else
|
{
|
y = d / net->Link[k].Len * Re * qual->Sc;
|
Sh = 3.65 + 0.0668 * y / (1.0 + 0.04 * pow(y, 0.667));
|
}
|
|
// Compute mass transfer coeff. (in ft/sec)
|
kf = Sh * qual->Diffus / d;
|
|
// For zero-order reaction, return mass transfer coeff.
|
if (qual->WallOrder == 0.0) return kf;
|
|
// For first-order reaction, return apparent wall coeff.
|
kw = net->Link[k].Kw / pr->Ucf[ELEV]; // Wall coeff, ft/sec
|
kw = (4.0 / d) * kw * kf / (kf + fabs(kw)); // Wall coeff, 1/sec
|
return kw;
|
}
|
|
|
double pipereact(Project *pr, int k, double c, double v, long dt)
|
/*
|
**------------------------------------------------------------
|
** Input: k = link index
|
** c = current quality in segment
|
** v = segment volume
|
** dt = time step
|
** Output: returns new WQ value
|
** Purpose: computes new quality in a pipe segment after
|
** reaction occurs
|
**------------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Quality *qual = &pr->quality;
|
|
double cnew, dc, dcbulk, dcwall, rbulk, rwall;
|
|
// For water age (hrs), update concentration by timestep
|
if (qual->Qualflag == AGE)
|
{
|
dc = (double)dt / 3600.0;
|
cnew = c + dc;
|
cnew = MAX(0.0, cnew);
|
return cnew;
|
}
|
|
// Otherwise find bulk & wall reaction rates
|
rbulk = bulkrate(pr, c, net->Link[k].Kb, qual->BulkOrder) * qual->Bucf;
|
rwall = wallrate(pr, c, net->Link[k].Diam, net->Link[k].Kw, net->Link[k].Rc);
|
|
// Find change in concentration over timestep
|
dcbulk = rbulk * (double)dt;
|
dcwall = rwall * (double)dt;
|
|
// Update cumulative mass reacted
|
if (pr->times.Htime >= pr->times.Rstart)
|
{
|
qual->Wbulk += fabs(dcbulk) * v;
|
qual->Wwall += fabs(dcwall) * v;
|
}
|
|
// Update concentration
|
dc = dcbulk + dcwall;
|
cnew = c + dc;
|
cnew = MAX(0.0, cnew);
|
return cnew;
|
}
|
|
|
double tankreact(Project *pr, double c, double v, double kb, long dt)
|
/*
|
**-------------------------------------------------------
|
** Input: c = current quality in tank
|
** v = tank volume
|
** kb = reaction coeff.
|
** dt = time step
|
** Output: returns new WQ value
|
** Purpose: computes new quality in a tank after
|
** reaction occurs
|
**-------------------------------------------------------
|
*/
|
{
|
Quality *qual = &pr->quality;
|
|
double cnew, dc, rbulk;
|
|
// For water age, update concentration by timestep
|
if (qual->Qualflag == AGE)
|
{
|
dc = (double)dt / 3600.0;
|
}
|
|
// For chemical analysis apply bulk reaction rate
|
else
|
{
|
// Find bulk reaction rate
|
rbulk = bulkrate(pr, c, kb, qual->TankOrder) * qual->Tucf;
|
|
// Find concentration change & update quality
|
dc = rbulk * (double)dt;
|
if (pr->times.Htime >= pr->times.Rstart)
|
{
|
qual->Wtank += fabs(dc) * v;
|
}
|
}
|
cnew = c + dc;
|
cnew = MAX(0.0, cnew);
|
return cnew;
|
}
|
|
|
double bulkrate(Project *pr, double c, double kb, double order)
|
/*
|
**-----------------------------------------------------------
|
** Input: c = current WQ concentration
|
** kb = bulk reaction coeff.
|
** order = bulk reaction order
|
** Output: returns bulk reaction rate
|
** Purpose: computes bulk reaction rate (mass/volume/time)
|
**-----------------------------------------------------------
|
*/
|
{
|
Quality *qual = &pr->quality;
|
|
double c1;
|
|
// Find bulk reaction potential taking into account
|
// limiting potential & reaction order.
|
|
// Zero-order kinetics:
|
if (order == 0.0) c = 1.0;
|
|
// Michaelis-Menton kinetics:
|
else if (order < 0.0)
|
{
|
c1 = qual->Climit + SGN(kb) * c;
|
if (fabs(c1) < TINY) c1 = SGN(c1) * TINY;
|
c = c / c1;
|
}
|
|
// N-th order kinetics:
|
else
|
{
|
// Account for limiting potential
|
if (qual->Climit == 0.0) c1 = c;
|
else c1 = MAX(0.0, SGN(kb) * (qual->Climit - c));
|
|
// Compute concentration potential
|
if (order == 1.0) c = c1;
|
else if (order == 2.0) c = c1 * c;
|
else c = c1 * pow(MAX(0.0, c), order - 1.0);
|
}
|
|
// Reaction rate = bulk coeff. * potential
|
if (c < 0) c = 0;
|
return kb * c;
|
}
|
|
|
double wallrate(Project *pr, double c, double d, double kw, double kf)
|
/*
|
**------------------------------------------------------------
|
** Input: c = current WQ concentration
|
** d = pipe diameter
|
** kw = intrinsic wall reaction coeff.
|
** kf = mass transfer coeff. for 0-order reaction
|
** (ft/sec) or apparent wall reaction coeff.
|
** for 1-st order reaction (1/sec)
|
** Output: returns wall reaction rate in mass/ft3/sec
|
** Purpose: computes wall reaction rate
|
**------------------------------------------------------------
|
*/
|
{
|
Quality *qual = &pr->quality;
|
|
if (kw == 0.0 || d == 0.0) return (0.0);
|
|
if (qual->WallOrder == 0.0) // 0-order reaction */
|
{
|
kf = SGN(kw) * c * kf; //* Mass transfer rate (mass/ft2/sec)
|
kw = kw * SQR(pr->Ucf[ELEV]); // Reaction rate (mass/ft2/sec)
|
if (fabs(kf) < fabs(kw)) kw = kf; // Reaction mass transfer limited
|
return (kw * 4.0 / d); // Reaction rate (mass/ft3/sec)
|
}
|
else return (c * kf); // 1st-order reaction
|
}
|
|
|
double mixtank(Project *pr, int n, double volin, double massin, double volout)
|
/*
|
**------------------------------------------------------------
|
** Input: n = node index
|
** volin = inflow volume to tank over time step
|
** massin = mass inflow to tank over time step
|
** volout = outflow volume from tank over time step
|
** Output: returns new quality for tank
|
** Purpose: mixes inflow with tank's contents to update its quality.
|
**------------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
|
int i;
|
double vnet;
|
i = n - net->Njuncs;
|
vnet = volin - volout;
|
switch (net->Tank[i].MixModel)
|
{
|
case MIX1: tankmix1(pr, i, volin, massin, vnet); break;
|
case MIX2: tankmix2(pr, i, volin, massin, vnet); break;
|
case FIFO: tankmix3(pr, i, volin, massin, vnet); break;
|
case LIFO: tankmix4(pr, i, volin, massin, vnet); break;
|
}
|
return net->Tank[i].C;
|
}
|
|
|
void tankmix1(Project *pr, int i, double vin, double win, double vnet)
|
/*
|
**---------------------------------------------
|
** Input: i = tank index
|
** vin = inflow volume
|
** win = mass inflow
|
** vnet = inflow - outflow
|
** Output: none
|
** Purpose: updates quality in a complete mix tank model
|
**---------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Quality *qual = &pr->quality;
|
|
int k;
|
double vnew;
|
Pseg seg;
|
Stank *tank = &net->Tank[i];
|
|
k = net->Nlinks + i;
|
seg = qual->FirstSeg[k];
|
if (seg)
|
{
|
vnew = seg->v + vin;
|
if (vnew > 0.0) seg->c = (seg->c * seg->v + win) / vnew;
|
seg->v += vnet;
|
seg->v = MAX(0.0, seg->v);
|
tank->C = seg->c;
|
}
|
}
|
|
|
void tankmix2(Project *pr, int i, double vin, double win, double vnet)
|
/*
|
**------------------------------------------------
|
** Input: i = tank index
|
** vin = inflow volume
|
** win = mass inflow
|
** vnet = inflow - outflow
|
** Output: none
|
** Purpose: updates quality in a 2-compartment tank model
|
**------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Quality *qual = &pr->quality;
|
|
int k;
|
double vt, // Transferred volume
|
vmz; // Full mixing zone volume
|
Pseg mixzone, // Mixing zone segment
|
stagzone; // Stagnant zone segment
|
Stank *tank = &pr->network.Tank[i];
|
|
// Identify segments for each compartment
|
k = net->Nlinks + i;
|
mixzone = qual->LastSeg[k];
|
stagzone = qual->FirstSeg[k];
|
if (mixzone == NULL || stagzone == NULL) return;
|
|
// Full mixing zone volume
|
vmz = tank->V1max;
|
|
// Tank is filling
|
vt = 0.0;
|
if (vnet > 0.0)
|
{
|
vt = MAX(0.0, (mixzone->v + vnet - vmz));
|
if (vin > 0.0)
|
{
|
mixzone->c = ((mixzone->c) * (mixzone->v) + win) /
|
(mixzone->v + vin);
|
}
|
if (vt > 0.0)
|
{
|
stagzone->c = ((stagzone->c) * (stagzone->v) +
|
(mixzone->c) * vt) / (stagzone->v + vt);
|
}
|
}
|
|
// Tank is emptying
|
else if (vnet < 0.0)
|
{
|
if (stagzone->v > 0.0) vt = MIN(stagzone->v, (-vnet));
|
if (vin + vt > 0.0)
|
{
|
mixzone->c = ((mixzone->c) * (mixzone->v) + win +
|
(stagzone->c) * vt) / (mixzone->v + vin + vt);
|
}
|
}
|
|
// Update segment volumes
|
if (vt > 0.0)
|
{
|
mixzone->v = vmz;
|
if (vnet > 0.0) stagzone->v += vt;
|
else stagzone->v = MAX(0.0, ((stagzone->v) - vt));
|
}
|
else
|
{
|
mixzone->v += vnet;
|
mixzone->v = MIN(mixzone->v, vmz);
|
mixzone->v = MAX(0.0, mixzone->v);
|
stagzone->v = 0.0;
|
}
|
|
// Use quality of mixing zone to represent quality of
|
// tank since this is where outflow begins to flow from
|
tank->C = mixzone->c;
|
}
|
|
|
void tankmix3(Project *pr, int i, double vin, double win, double vnet)
|
/*
|
**----------------------------------------------------------
|
** Input: i = tank index
|
** vin = inflow volume
|
** win = mass inflow
|
** vnet = inflow - outflow
|
** Output: none
|
** Purpose: Updates quality in a First-In-First-Out (FIFO) tank model.
|
**----------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Quality *qual = &pr->quality;
|
|
int k;
|
double vout, vseg;
|
double cin, vsum, wsum;
|
Pseg seg;
|
Stank *tank = &pr->network.Tank[i];
|
|
k = net->Nlinks + i;
|
if (qual->LastSeg[k] == NULL || qual->FirstSeg[k] == NULL) return;
|
|
// Add new last segment for flow entering the tank
|
if (vin > 0.0)
|
{
|
// ... increase segment volume if inflow has same quality as segment
|
cin = win / vin;
|
seg = qual->LastSeg[k];
|
if (fabs(seg->c - cin) < qual->Ctol) seg->v += vin;
|
|
// ... otherwise add a new last segment to the tank
|
else addseg(pr, k, vin, cin);
|
}
|
|
// Withdraw flow from first segment
|
vsum = 0.0;
|
wsum = 0.0;
|
vout = vin - vnet;
|
while (vout > 0.0)
|
{
|
seg = qual->FirstSeg[k];
|
if (seg == NULL) break;
|
vseg = seg->v; // Flow volume from leading seg
|
vseg = MIN(vseg, vout);
|
if (seg == qual->LastSeg[k]) vseg = vout;
|
vsum += vseg;
|
wsum += (seg->c) * vseg;
|
vout -= vseg; // Remaining flow volume
|
if (vout >= 0.0 && vseg >= seg->v) // Seg used up
|
{
|
if (seg->prev)
|
{
|
qual->FirstSeg[k] = seg->prev;
|
seg->prev = qual->FreeSeg;
|
qual->FreeSeg = seg;
|
}
|
}
|
else seg->v -= vseg; // Remaining volume in segment
|
}
|
|
// Use quality withdrawn from 1st segment
|
// to represent overall quality of tank
|
if (vsum > 0.0) tank->C = wsum / vsum;
|
else if (qual->FirstSeg[k] == NULL) tank->C = 0.0;
|
else tank->C = qual->FirstSeg[k]->c;
|
}
|
|
|
void tankmix4(Project *pr, int i, double vin, double win, double vnet)
|
/*
|
**----------------------------------------------------------
|
** Input: i = tank index
|
** vin = inflow volume
|
** win = mass inflow
|
** vnet = inflow - outflow
|
** Output: none
|
** Purpose: Updates quality in a Last In-First Out (LIFO) tank model.
|
**----------------------------------------------------------
|
*/
|
{
|
Network *net = &pr->network;
|
Quality *qual = &pr->quality;
|
|
int k;
|
double cin, vsum, wsum, vseg;
|
Pseg seg;
|
Stank *tank = &pr->network.Tank[i];
|
|
k = net->Nlinks + i;
|
if (qual->LastSeg[k] == NULL || qual->FirstSeg[k] == NULL) return;
|
|
// Find inflows & outflows
|
if (vin > 0.0) cin = win / vin;
|
else cin = 0.0;
|
|
// If tank filling, then create new last seg
|
tank->C = qual->LastSeg[k]->c;
|
seg = qual->LastSeg[k];
|
if (vnet > 0.0)
|
{
|
// ... inflow quality is same as last segment's quality,
|
// so just add inflow volume to last segment
|
if (fabs(seg->c - cin) < qual->Ctol) seg->v += vnet;
|
|
// ... otherwise add a new last segment with inflow quality
|
else addseg(pr, k, vnet, cin);
|
|
// Update reported tank quality
|
tank->C = qual->LastSeg[k]->c;
|
}
|
|
// If tank emptying then remove last segments until vnet consumed
|
else if (vnet < 0.0)
|
{
|
vsum = 0.0;
|
wsum = 0.0;
|
vnet = -vnet;
|
|
// Reverse segment chain so segments are processed from last to first
|
reversesegs(pr, k);
|
|
// While there is still volume to remove
|
while (vnet > 0.0)
|
{
|
// ... start with reversed first segment
|
seg = qual->FirstSeg[k];
|
if (seg == NULL) break;
|
|
// ... find volume to remove from it
|
vseg = seg->v;
|
vseg = MIN(vseg, vnet);
|
if (seg == qual->LastSeg[k]) vseg = vnet;
|
|
// ... update total volume & mass removed
|
vsum += vseg;
|
wsum += (seg->c) * vseg;
|
|
// ... update remiaing volume to remove
|
vnet -= vseg;
|
|
// ... if no more volume left in current segment
|
if (vnet >= 0.0 && vseg >= seg->v)
|
{
|
// ... replace current segment with previous one
|
if (seg->prev)
|
{
|
qual->FirstSeg[k] = seg->prev;
|
seg->prev = qual->FreeSeg;
|
qual->FreeSeg = seg;
|
}
|
}
|
|
// ... otherwise reduce volume of current segment
|
else seg->v -= vseg;
|
}
|
|
// Restore original orientation of segment chain
|
reversesegs(pr, k);
|
|
// Reported tank quality is mixture of flow released and any inflow
|
tank->C = (wsum + win) / (vsum + vin);
|
}
|
}
|
