SWMM 5.2.2 Code for LID Storage Layer

 

Variable Name	Type	Description
thickness	double	the layer thickness (ft)
voidFrac	double	the void volume / total volume
kSat	double	saturated hydraulic conductivity (ft/sec)
clogFactor	double	clogging factor
covered	int	Boolean flag indicating if rain barrel is covered (TRUE if it is covered)

This code defines a struct called TStorageLayer, which represents a storage layer in a Low Impact Development (LID) system. The struct contains several variables:

• thickness: a double that represents the thickness of the layer, in feet.
• voidFrac: a double that represents the void volume of the layer over the total volume of the layer.
• kSat: a double that represents the saturated hydraulic conductivity, the rate at which water can flow through the soil, in feet per second.
• clogFactor: a double that represents the clogging factor, which can be used to account for the potential reduction of the saturated hydraulic conductivity due to clogging or other factors.
• covered: an integer that is a Boolean flag indicating if the rain barrel is covered. This is set to TRUE if the rain barrel is covered and FALSE otherwise.

The struct is used to store information about the LID storage layer, which can then be used in various calculations and simulations related to the LID system. The struct is usually used along with other parameters such as precipitation, infiltration, and evapotranspiration to evaluate the water management of the LID system.

// LID Storage Layer typedef struct { double thickness; // layer thickness (ft) double voidFrac; // void volume / total volume double kSat; // saturated hydraulic conductivity (ft/sec) double clogFactor; // clogging factor int covered; // TRUE if rain barrel is covered } TStorageLayer;

SWMM 5.2.2 Code for Drainage Mat Layer (for green roofs)

 

Variable Name	Type	Description
thickness	double	the layer thickness (ft)
voidFrac	double	the void volume / total volume
roughness	double	Mannings n for green roof drainage mats
alpha	double	slope/roughness term in Manning equation

This code defines a struct called TDrainMatLayer, which represents a Drainage Mat Layer commonly used in green roofs. The struct contains several variables:

• thickness: a double that represents the thickness of the layer, in feet.
• voidFrac: a double that represents the void volume of the layer over the total volume of the layer.
• roughness: a double that represents the Mannings n, a roughness coefficient used in the calculation of flow through open channels, for green roof drainage mats.
• alpha: a double that represents the slope/roughness term in the Manning equation, used to estimate the resistance to flow in open channels.

The struct is used to store information about the Drainage Mat Layer, which can then be used in various calculations and simulations related to green roofs. This struct is usually used along with other layers and other parameters to evaluate the water management of green roofs.

// Drainage Mat Layer (for green roofs)

typedef struct

{

    double    thickness;          // layer thickness (ft)

    double    voidFrac;           // void volume / total volume

    double    roughness;          // Mannings n for green roof drainage mats

    double    alpha;              // slope/roughness term in Manning equation

}  TDrainMatLayer;

SWMM 5.2.2 Underdrain System (part of Storage Layer) Code

 

Variable Name	Type	Description
coeff	double	underdrain flow coefficient (in/hr or mm/hr)
expon	double	underdrain head exponent (for in or mm)
offset	double	offset height of underdrain (ft)
delay	double	rain barrel drain delay time (sec)
hOpen	double	head when drain opens (ft)
hClose	double	head when drain closes (ft)
qCurve	int	curve controlling flow rate (optional)

This code defines a struct called TDrainLayer, which represents an underdrain system as part of a storage layer. The struct contains several variables:

• coeff: a double representing the underdrain flow coefficient, in units of in/hr or mm/hr.
• expon: a double representing the underdrain head exponent, used to calculate the head loss through the underdrain.
• offset: a double representing the offset height of the underdrain, in feet.
• delay: a double that represents the delay time for the rain barrel drain, in seconds.
• hOpen: a double representing the head at which the drain opens, in feet.
• hClose: a double representing the head at which the drain closes, in feet.
• qCurve: an integer that represents the curve controlling the flow rate of the underdrain. This is optional and may not be used in all cases.

The struct is used to store information about the underdrain system, which can then be used in various calculations and simulations related to the storage layer.

// Underdrain System (part of Storage Layer) typedef struct { double coeff; // underdrain flow coeff. (in/hr or mm/hr) double expon; // underdrain head exponent (for in or mm) double offset; // offset height of underdrain (ft) double delay; // rain barrel drain delay time (sec) double hOpen; // head when drain opens (ft) double hClose; // head when drain closes (ft) int qCurve; // curve controlling flow rate (optional) } TDrainLayer;

Runge-Kutta 2nd order method (RK2), also known as the Midpoint Method

The Runge-Kutta 2nd order method (RK2), also known as the midpoint method, is a numerical method for solving ordinary differential equations (ODEs). The basic idea behind this method is to approximate the solution of the ODE at discrete time steps, using a local polynomial approximation of the solution.

The method begins by approximating the solution at the current time step using the initial value and the derivative of the solution at that point, which is given by the ODE. The next step is to calculate the estimate of the solution at the midpoint between the current and next time steps, this is done by using the initial value and the derivative evaluated at the current time step .

The next approximate solution is then calculated by taking a weighted average of the estimates. The weights used in this average are determined by the specific RK method being used.

The RK2 method is a second-order method, meaning that it can achieve an accuracy of O(h^2), where h is the size of the time step. However, the RK2 method is less computationally expensive than higher-order methods like the RK4 method, which is a fourth-order method with an accuracy of O(h^4).

The explicit form of RK2 method is as follows:

y(tn+1) = y(tn) + h * f(tn + h/2, y(tn) + (h/2)*f(tn, y(tn)))

where y is the function to be solved, h is the step size and f(tn, y(tn)) is the derivative of the function.

It's important to note that RK2 method is not always the best choice for solving ODEs, especially if high accuracy is needed or stiff ODEs. The choice of numerical method depends on the specific problem being solved and the desired level of accuracy.

All RK Methods

Runge-Kutta methods are a family of iterative methods for solving ordinary differential equations (ODEs) with a given initial value. The basic idea behind these methods is to approximate the solution of the ODE at discrete time steps, using a local polynomial approximation of the solution.

The most common method in this family is the fourth-order Runge-Kutta (RK4) method, which uses four estimates of the solution at different points within a time step to compute the next approximate solution. The method begins by approximating the solution at the current time step using the initial value and the derivative of the solution at that point, which is given by the ODE.

The first estimate of the solution at the next time step is calculated using the initial value and the derivative evaluated at the current time step. The second estimate is calculated using the initial value and the derivative evaluated at the midpoint between the current and next time steps. The third and fourth estimates are calculated in a similar manner, using the derivative evaluated at other points within the time step.

These four estimates are then used to compute the next approximate solution, by taking a weighted average of the estimates. The weights used in this average are determined by the specific RK method being used.

The RK4 method is a fourth-order method, meaning that it can achieve an accuracy of O(h^4), where h is the size of the time step. However, the RK4 method is more computationally expensive than lower-order methods like the Euler method, which is a first-order method with an accuracy of O(h).

There are other Runge Kutta Methods like RK2, RK3, RK4, RK5 and so on. It's important to note that RK methods are not always the best choice for solving ODEs. The choice of numerical method depends on the specific problem being solved and the desired level of accuracy.

External Source

https://www.researchgate.net/figure/Illustration-of-numerical-discretization-methods-a-FE-b-RK2-c-RK4-and-d-AB2_fig1_332343636

Illustration of numerical discretization methods: (a) FE, (b) RK2, (c) RK4, and (d) AB2 Figure 1(b) shows the application of RK2 method where Predictor and Corrector are the points of the function at each sample based on Equations (4) and (5). From this perspective, the RK2 method can also be considered as a hybrid method with switching factor is equal to one. This mechanism can be seen from Equations (4) and (5) where the the predictor 'measures' value and then the corrector update the values and this switching continues for the complete ODE function. It is obvious that the main challenge with this approach is the high effort for computation is required even the accuracy is improved. The reduced speed calculation due to the finding the solution is the main hindrance of adaption the method in real-time control system.

SWMM 5.2.2 Code for LID Function getEvapRates

This code is for calculating evaporation rates for a low impact development (LID) unit, specifically for the surface layer, pavement, soil, and storage layer. It takes in the volumes/areas of water in each layer (surfaceVol, paveVol, soilVol, and storageVol) as well as the fraction of the surface layer that is pervious (pervFrac) as inputs. It also uses a global variable EvapRate. The function first calculates the available evaporation rate and then calculates the evaporation rate for each layer using the MIN and MAX functions to ensure the evaporation rate does not exceed the available evaporation rate or the volume of water in the layer. The evaporation rate for the surface layer is multiplied by the fraction of the surface layer that is pervious. If the SurfaceInfil variable is greater than zero, the evaporation rates for the pavement, soil, and storage layers are set to zero. The code also uses Tstep variable to convert the rate to total volume over the current time step. The output of this function is updating the member variable SurfaceEvap, PaveEvap, SoilEvap, and StorageEvap and doesn't return any value.

void getEvapRates(double surfaceVol, double paveVol, double soilVol,

    double storageVol, double pervFrac)

//

//  Purpose: computes surface, pavement, soil, and storage evaporation rates.

//  Input:   surfaceVol = volume/area of ponded water on surface layer (ft)

//           paveVol    = volume/area of water in pavement pores (ft)

//           soilVol    = volume/area of water in soil (or pavement) pores (ft)

//           storageVol = volume/area of water in storage layer (ft)

//           pervFrac   = fraction of surface layer that is pervious

//  Output:  none

//

{

    double availEvap;

    //... surface evaporation flux

    availEvap = EvapRate;

    SurfaceEvap = MIN(availEvap, surfaceVol/Tstep);

    SurfaceEvap = MAX(0.0, SurfaceEvap);

    availEvap = MAX(0.0, (availEvap - SurfaceEvap));

    availEvap *= pervFrac;

    //... no subsurface evap if water is infiltrating

    if ( SurfaceInfil > 0.0 )

    {

        PaveEvap = 0.0;

        SoilEvap = 0.0;

        StorageEvap = 0.0;

    }

    else

    {

        //... pavement evaporation flux

        PaveEvap = MIN(availEvap, paveVol / Tstep);

        availEvap = MAX(0.0, (availEvap - PaveEvap));

        //... soil evaporation flux

        SoilEvap = MIN(availEvap, soilVol / Tstep);

        availEvap = MAX(0.0, (availEvap - SoilEvap));

        //... storage evaporation flux

        StorageEvap = MIN(availEvap, storageVol / Tstep);

    }

}

}