Title: | Steady and Unsteady Open-Channel Flow Computation |
---|---|
Description: | A tool for undergraduate and graduate courses in open-channel hydraulics. Provides functions for computing normal and critical depths, steady-state water surface profiles (e.g. backwater curves) and unsteady flow computations (e.g. flood wave routing) as described in Koohafkan MC, Younis BA (2015). "Open-channel computation with R." The R Journal, 7(2), 249–262. <doi: 10.32614/RJ-2015-034>. |
Authors: | Michael C Koohafkan [aut, cre] |
Maintainer: | Michael C Koohafkan <[email protected]> |
License: | GPL (>= 3) |
Version: | 1.2-3 |
Built: | 2025-03-02 04:34:15 UTC |
Source: | https://github.com/mkoohafkan/rivr |
This package is designed as an educational tool for students and instructors of undergraduate courses in open channel hydraulics. Functions are provided for computing normal and critical depths, steady (e.g. backwater curves) and unsteady (flood wave routing) flow computations for prismatic trapezoidal channels. See the vignettes to get started.
Compute geometry relations for trapezoidal channels.
channel_geom(y, B, SS)
channel_geom(y, B, SS)
y |
Flow depth [ |
B |
Channel bottom width [ |
SS |
Channel sideslope [ |
Channel geometry relations are routinely calculated in numerical solutions of steady, gradually-varied and unsteady flows. This function is used extensively by internal procedures and is made accessible to the user for convenience.
Named vector:
A |
Flow area [ |
P |
Wetted perimeter [ |
R |
Hydraulic radius [ |
dAdy |
Water surface width [ |
dPdy |
First derivative of wetted perimeter w.r.t. flow depth. |
dRdy |
First derivative of hydraulic radius w.r.t. flow depth. |
DH |
Hydraulic depth [ |
ybar |
Vertical distance from water surface to centroid of flow area [ |
channel_geom(1.71, 100, 0) # rectangular channel channel_geom(5.79, 6.1, 1.5) # trapezoidal channel with sideslope 3H:2V
channel_geom(1.71, 100, 0) # rectangular channel channel_geom(5.79, 6.1, 1.5) # trapezoidal channel with sideslope 3H:2V
Compute the gradually-varied flow profile of a prismatic channel.
compute_profile( So, n, Q, y0, Cm, g, B, SS, z0 = 0, x0 = 0, stepdist, totaldist )
compute_profile( So, n, Q, y0, Cm, g, B, SS, z0 = 0, x0 = 0, stepdist, totaldist )
So |
Channel slope [ |
n |
Manning's roughness coefficient. |
Q |
Flow rate [ |
y0 |
The water depth at the control section [ |
Cm |
Unit conversion coefficient for Manning's equation. For SI units, Cm = 1. |
g |
Gravitational acceleration [ |
B |
Channel bottom width [ |
SS |
Channel sideslope [ |
z0 |
Elevation reference datum at control section [ |
x0 |
Distance reference at control section [ |
stepdist |
The spatial interval used in the Standard step method [ |
totaldist |
The total distance upstream (or downstream) to compute the profile [ |
Computes the longitudinal water surface profile of a prismatic channel using the standard step method by solving the non-linear ODE
The standard-step method operates by stepping along the channel by a constant distance interval, starting from a cross-section where the flow depth is known (the control section). The flow depth is computed at the adjacent cross-section (target section). The computed value at the target is then used as the basis for computing flow depth at the next cross-section, i.e. the previous target section becomes the new control section for each step. A Newton-Raphson scheme is used each step to compute the flow depth and friction slope. Technically, the average friction slope of the control and target section is used to compute the flow depth at the target section.
data.frame with columns:
x |
Along-channel distance. |
z |
Elevation. |
y |
Flow depth. |
v |
Flow velocity. |
A |
Flow area. |
Sf |
Friction slope. |
E |
Total energy. |
Fr |
Froude Number. |
# example M1 profile compute_profile(0.001, 0.045, 250, 2.7, 1.486, 32.2, 100, 0, stepdist = 10, totaldist = 3000) # example M2 profile compute_profile(0.001, 0.045, 250, 0.64, 1.486, 32.2, 100, 0, stepdist = 10, totaldist = 3000) # example S2 profile compute_profile(0.005, 0.01, 250, 2.65, 1.486, 32.2, 10, 0, stepdist = 10, totaldist = 2000) # example S3 profile compute_profile(0.005, 0.01, 250, 0.5, 1.486, 32.2, 10, 0, stepdist = 10, totaldist = 2000)
# example M1 profile compute_profile(0.001, 0.045, 250, 2.7, 1.486, 32.2, 100, 0, stepdist = 10, totaldist = 3000) # example M2 profile compute_profile(0.001, 0.045, 250, 0.64, 1.486, 32.2, 100, 0, stepdist = 10, totaldist = 3000) # example S2 profile compute_profile(0.005, 0.01, 250, 2.65, 1.486, 32.2, 10, 0, stepdist = 10, totaldist = 2000) # example S3 profile compute_profile(0.005, 0.01, 250, 0.5, 1.486, 32.2, 10, 0, stepdist = 10, totaldist = 2000)
Calculate the channel conveyance.
conveyance(n, A, R, Cm)
conveyance(n, A, R, Cm)
n |
Manning's roughness coefficient (dimensionless). |
A |
Flow area [ |
R |
Hydraulic radius [ |
Cm |
Unit conversion coefficient for Manning's equation. For SI units, Cm = 1. |
Channel conveyance is routinely calculated in numerical solutions of steady, gradually-varied and unsteady flows. This function is used extensively by internal procedures and is made accessible to the user for convenience.
The channel conveyance.
Calculate the critical depth.
critical_depth(Q, yopt, g, B, SS)
critical_depth(Q, yopt, g, B, SS)
Q |
Flow rate [ |
yopt |
Initial guess for normal depth [ |
g |
Gravitational acceleration [ |
B |
Channel bottom width [ |
SS |
Channel sideslope [ |
The critical depth is the water depth at which a channel flow regime will transition from supercritical to subcritical (or vice versa). Calculation of the critical depth is based on a specific energy formulation, i.e.
where is the flow depth,
is
the elevation relative to some datum (assumed to be 0), and the last term
represents kinetic energy. More specifically, the function operates
by finding the point where the derivative of specific energy w.r.t.
is zero, i.e.
when
.
The critical depth [
].
critical_depth(250, 2, 32.2, 100, 0) # rectangular channel critical_depth(126, 1, 9.81, 6.1, 1.5) # trapezoidal channel with sideslope 3H:2V
critical_depth(250, 2, 32.2, 100, 0) # rectangular channel critical_depth(126, 1, 9.81, 6.1, 1.5) # trapezoidal channel with sideslope 3H:2V
Demonstrate package functionality via Shiny apps
demo_shiny(ex)
demo_shiny(ex)
ex |
Example to run. |
Demonstrations available:
"gvf"
Gradually-varied flow.
## Not run: # get list of available demos demo_shiny() # run the gradually-varied flow demo demo_shiny("gvf") ## End(Not run)
## Not run: # get list of available demos demo_shiny() # run the gradually-varied flow demo demo_shiny("gvf") ## End(Not run)
Calculate the Froude Number.
froude(Q, g, A, DH)
froude(Q, g, A, DH)
Q |
Flow rate [ |
g |
Gravitational acceleration [ |
A |
Flow area [ |
DH |
Hydraulic depth [ |
The Froude number is a dimensionless measure of bulk flow characteristics that represents the relative importance of inertial forces and gravitational forces. For open channel flow, the Froude number of open channel flow is defined as
where is the flow velocity,
is the gravitational
acceleration and
is the hydraulic depth. The Froude number is related
to the energy state of the flow and can be used to identify flows as
either supercritical (
) or subcritical (
).
The Froude Number (dimensionless).
froude(250, 32.2, 171, 1.71) # subcritical flow froude(250, 32.2, 57.9, 0.579) # critical flow froude(250, 32.2, 45, 0.45) # supercritical flow
froude(250, 32.2, 171, 1.71) # subcritical flow froude(250, 32.2, 57.9, 0.579) # critical flow froude(250, 32.2, 45, 0.45) # supercritical flow
Calculate the normal (equilibrium) depth using Manning's equation.
normal_depth(So, n, Q, yopt, Cm, B, SS)
normal_depth(So, n, Q, yopt, Cm, B, SS)
So |
Channel slope [ |
n |
Manning's roughness coefficient. |
Q |
Flow rate [ |
yopt |
Initial guess for normal depth [ |
Cm |
Unit conversion coefficient for Manning's equation. For SI units, Cm = 1. |
B |
Channel bottom width [ |
SS |
Channel sideslope [ |
The normal depth is the equilibrium depth of a channel for a given flow rate, channel slope, geometry and roughness. Manning's equation is used to calculate the equilibrium depth. Manning's equation for normal flow is defined as
where is the channel flow,
is the channel slope,
is the
cross-sectional flow area,
is the hydraulic depth and
is a conversion factor
based on the unit system used. This function uses a Newton-Raphson root-finding approach
to calculate the normal depth, i.e.
when
.
The normal depth [
].
normal_depth(0.001, 0.045, 250, 3, 1.486, 100, 0) # rectangular channel normal_depth(0.0008, 0.013, 126, 5, 1, 6.1, 1.5) # trapezoidal channel with sideslope 3H:2V
normal_depth(0.001, 0.045, 250, 3, 1.486, 100, 0) # rectangular channel normal_depth(0.0008, 0.013, 126, 5, 1, 6.1, 1.5) # trapezoidal channel with sideslope 3H:2V
Route a flood wave down a prismatic channel.
route_wave( So, n, Cm, g, B, SS, initial.condition, boundary.condition, downstream.condition, timestep, spacestep, numnodes, monitor.nodes, monitor.times, engine = c("Dynamic", "Kinematic"), scheme = c("MacCormack", "Lax"), boundary.type = c("QQ", "Qy", "yQ", "yy") )
route_wave( So, n, Cm, g, B, SS, initial.condition, boundary.condition, downstream.condition, timestep, spacestep, numnodes, monitor.nodes, monitor.times, engine = c("Dynamic", "Kinematic"), scheme = c("MacCormack", "Lax"), boundary.type = c("QQ", "Qy", "yQ", "yy") )
So |
Channel slope [ |
n |
Manning's roughness coefficient. |
Cm |
Unit conversion coefficient for Manning's equation. For SI units, Cm = 1. |
g |
Gravitational acceleration [ |
B |
Channel bottom width [ |
SS |
Channel sideslope [ |
initial.condition |
The initial flow rate [ |
boundary.condition |
Vector specifying the upstream boundary condition
for the full duration of the model. If |
downstream.condition |
Only used if |
timestep |
Temporal resolution of the model. Also the assumed time interval [ |
spacestep |
the spatial resolution of the model, interpreted as the distance [ |
numnodes |
The number of nodes used to discretize the channel. The total channel extent is
computed as |
monitor.nodes |
the nodes to be monitored every time step. Specified as a vector of node
indices, with 1 being the upstream boundary and |
monitor.times |
the time steps at which to monitor every node. Specified as a vector of
indices of |
engine |
The engine to be used for routing the flood wave. May be either "Kinematic" or "Dynamic". |
scheme |
Only used if |
boundary.type |
Only used if |
Provides implementations of a Kinematic Wave Model (KWM) and a Dynamic Wave Model (DWM) with the choice of two numerical schemes. The MacCormack scheme is a second-order accurate predictor-corrector scheme that provides efficient flood wave routing. The Lax diffusive scheme can be used to obtain smooth solutions for problems with discontinuities in the boundary conditions, e.g. sudden gate closures. The DWM implementation uses the Method of Characteristics (MOC) to compute the flow regime at the model boundaries, and allows the user to specify boundaries in terms of depths and/or flows. the KWM implementation assumes the normal depth at the upstream boundary and is only first-order accurate.
data.frame with columns:
step |
Time step. |
node |
Node index. |
time |
Time since start. |
distance |
Downstream distance. |
flow |
Flow rate. |
depth |
Flow depth. |
velocity |
Flow velocity. |
area |
Flow area. |
monitor.type |
Row refers to a monitored node ("node") or timestep ("timestep"). |
## Not run: # kinematic wave routing times = seq(0, 30000, by = 25) floodwave = ifelse(times >= 9000, 250, 250 + (750/pi)*(1 - cos(pi*times/(60*75)))) route_wave(0.001, 0.045, 1.486, 32.2, 100, 0, initial.condition = 250, boundary.condition = floodwave, timestep = 25, spacestep = 50, numnodes=301, monitor.nodes = c(1, 101, 201, 301), monitor.times = seq(1, length(times), by = 10), engine = "Kinematic") # dynamic wave routing with zero-gradient downstream condition using MacCormack scheme route_wave(0.001, 0.045, 1.486, 32.2, 100, 0, initial.condition = 250, boundary.condition = floodwave, downstream.condition = rep(-1, length(times)), timestep = 25, spacestep = 500, numnodes = 31, engine = "Dynamic", scheme = "MacCormack", monitor.nodes = c(1, 11, 21, 31), monitor.times = seq(1, length(times), by = 10)) # mixed boundary conditions (sudden gate closure) using Lax scheme lax = route_wave(0.00008, 0.013, 1, 9.81, 6.1, 1.5, initial.condition = 126, boundary.condition = rep(5.79, 2001), downstream.condition = rep(0, 2001), timestep = 1, spacestep = 10, numnodes = 501, monitor.nodes = c(1, 151, 251, 301, 501), monitor.times = c(1, 501, 1001, 1501, 2001), engine="Dynamic", scheme="Lax", boundary.type="yQ") # extract data for a monitored point require(dplyr) filter(lax, monitor.type == "node", node == 151) ## End(Not run)
## Not run: # kinematic wave routing times = seq(0, 30000, by = 25) floodwave = ifelse(times >= 9000, 250, 250 + (750/pi)*(1 - cos(pi*times/(60*75)))) route_wave(0.001, 0.045, 1.486, 32.2, 100, 0, initial.condition = 250, boundary.condition = floodwave, timestep = 25, spacestep = 50, numnodes=301, monitor.nodes = c(1, 101, 201, 301), monitor.times = seq(1, length(times), by = 10), engine = "Kinematic") # dynamic wave routing with zero-gradient downstream condition using MacCormack scheme route_wave(0.001, 0.045, 1.486, 32.2, 100, 0, initial.condition = 250, boundary.condition = floodwave, downstream.condition = rep(-1, length(times)), timestep = 25, spacestep = 500, numnodes = 31, engine = "Dynamic", scheme = "MacCormack", monitor.nodes = c(1, 11, 21, 31), monitor.times = seq(1, length(times), by = 10)) # mixed boundary conditions (sudden gate closure) using Lax scheme lax = route_wave(0.00008, 0.013, 1, 9.81, 6.1, 1.5, initial.condition = 126, boundary.condition = rep(5.79, 2001), downstream.condition = rep(0, 2001), timestep = 1, spacestep = 10, numnodes = 501, monitor.nodes = c(1, 151, 251, 301, 501), monitor.times = c(1, 501, 1001, 1501, 2001), engine="Dynamic", scheme="Lax", boundary.type="yQ") # extract data for a monitored point require(dplyr) filter(lax, monitor.type == "node", node == 151) ## End(Not run)
Digitized results from the California Water Olympics. The variables are as follows:
t The time (in seconds) since the start of the model run.
Q The flow rate [].
x The distance downstream [] at which the hydrograph was recorded.
The data can be used to validate numerical solutions to flood wave routing for a channel under the following conditions:
Channel width is 100 feet.
Channel slope is 0.001.
Channel extent is 150,000 feet.
Channel roughness (Manning's n) is 0.045.
Channel sideslope is 0 (rectangular channel).
Initial flow rate is 250 cfs.
Upstream boundary condition is defined as
data(waterolympics)
data(waterolympics)
A data frame with 40 rows and 3 variables
Sobey, Rodney. "H11: Hydrograph Routing." Review of One-Dimensional Hydrodynamic and Transport Models. Bay-Delta Modeling Forum, 15 June 2001. Web. 13 Mar. 2015. <http://www.cwemf.org/1-DReview/>.