Nothing
R/manningc.R
In hydraulics: Basic Pipe and Open Channel Hydraulics
Defines functions
return_fcn2
n_adj_func
Qfull
ycfunc
#' Solves the Manning Equation for gravity flow in a circular pipe
#'
#' This function solves the Manning equation for water flow in a circular pipe
#' at less than full. Uniform flow conditions are assumed, so that the pipe slope
#' is equal to the slope of the water surface and the energy grade line. This is a
#' modification of the code prepared by Irucka Embry in his iemisc package.
#' The iemisc::manningcirc function was adapted here for more limited cases
#' commonly used in classroom exercises, additional checks were included to ensure
#' the pipe is flowing less than full, and a cross-section figure is also
#' available. The iemisc::manningcirc and iemisc::manningcircy functions were
#' combined into a single function. Manning n supplied is assumed to be that for
#' full pipe flow; an optional argument may be supplied to account for n varying
#' with depth.
#'
#' The possible applications of this function for solving the Manning equation
#' for circular pipes are:
#' \tabular{ll}{
#' \strong{Given} \tab \strong{Solve for} \cr
#' y_d, Q, Sf, n \tab d \cr
#' d, Sf, Q, n \tab y \cr
#' y, d, Q, n \tab Sf \cr
#' y, d, Sf, n \tab Q \cr
#' d, Q, Sf, y \tab n
#' }
#'
#' @param Q numeric vector that contains the flow rate [\eqn{m^3 s^{-1}}{m^3/s} or \eqn{ft^3 s^{-1}}{ft^3/s}]
#' @param d numeric vector that contains the pipe diameter [\eqn{m}{m} or \eqn{ft}{ft}]
#' @param Sf numeric vector that contains the slope of the pipe [unitless]
#' @param y numeric vector that contains the water depth [\eqn{m}{m} or \eqn{ft}{ft}]
#' @param y_d numeric vector that contains the ratio of depth to diameter [unitless]
#' @param n numeric vector that contains the Manning roughness coefficient (for full flow or fixed).
#' @param n_var If set to TRUE the value of n for full flow is adjusted with depth. [Default is FALSE]
#' @param units character vector that contains the system of units [options are
#' \code{SI} for International System of Units and \code{Eng} for English (US customary)
#' units. This is used for compatibility with iemisc package]
#' @param ret_units If set to TRUE the value(s) returned are of class \code{units} with
#' units attached to the value. [Default is FALSE]
#'
#' @return Returns a list including the missing parameter:
#' \itemize{
#' \item Q - flow rate
#' \item V - flow velocity
#' \item A - cross-sectional area of flow
#' \item P - wetted perimeter
#' \item R - hydraulic radius (A/P)
#' \item y - flow depth
#' \item d - pipe diameter
#' \item Sf - slope
#' \item n - Manning's roughness (for full flow, or as adjusted if n_var is TRUE)
#' \item yc - critical depth
#' \item Fr - Froude number
#' \item Re - Reynolds number
#' \item Qf - Full pipe flow rate
#' }
#'
#' @author Ed Maurer, Irucka Embry
#'
#' @details The Manning equation (also known as the Strickler equation) describes flow conditions in
#' an open channel under uniform flow conditions. It is often expressed as:
#' \deqn{Q = A\frac{C}{n}{R}^{\frac{2}{3}}{S_f}^{\frac{1}{2}}}
#' where \eqn{C} is 1.0 for SI units and 1.49 for Eng (U.S. Customary) units. Critical depth is
#' defined by the relation (at critical conditions):
#' \deqn{\frac{Q^{2}B}{g\,A^{3}}=1}{Q^2B/gA^3=1}
#' where \eqn{B}{B} is the top width of the water surface. Since B equals zero for a full pipe, critical
#' depth is set to the pipe diameter \eqn{d}{d} if the flow \eqn{Q}{Q} exceeds a value that would produce a
#' critical flow at \eqn{\frac{y}{d}=0.99}{y/d=0.99}.
#'
#' @examples
#'
#' #Solving for flow rate, Q: SI Units
#' manningc(d = 0.6, n = 0.013, Sf = 1./400., y = 0.24, units = "SI")
#' #returns 0.1 m3/s
#'
#' #Solving for Sf, if d=600 mm and pipe is to flow half full
#' manningc(d = 0.6, Q = 0.17, n = 0.013, y = 0.3, units = "SI")
#' #returns required slope of 0.003
#'
#' #Solving for diameter, d when given y_d): Eng (US) units
#' manningc(Q = 83.5, n = 0.015, Sf = 0.0002, y_d = 0.9, units = "Eng")
#' #returns 7.0 ft required diameter
#'
#' #Solving for depth, y when given Q: SI units
#' manningc(Q=0.01, n=0.013, Sf=0.001, d = 0.2, units="SI")
#' #returns depth y = 0.158 m, critical depth, yc = 0.085 m
#'
#' #Solving for depth, y when given Q: SI units, and n variable with depth
#' manningc(Q=0.01, n=0.013, Sf=0.001, d = 0.2, n_var = TRUE, units="SI")
#' #returns depth y = 0.174 m, critical depth, yc = 0.085 m
#'
#' @seealso \code{\link{xc_circle}} for a cross-section diagram of the circular channel
#'
#' @importFrom pracma fsolve
#'
#' @name manningc
NULL
# Function for finding critical depth
ycfunc <- function(Q = NULL, d = NULL, g = NULL) {
# Find theta for 99% full (y/D = 0.99)
thetafull <- 2 * acos(1 - (2 * (0.99)))
Qcfull <- sqrt(g * d^5 * ((thetafull - sin(thetafull))^3)/(8^3 * sin(thetafull/2.0)))
if ( Q >= Qcfull ) {
yc <- d
} else {
thetafun <- function(theta) {Q^2 / ( g * d^5 ) - ((theta - sin(theta))^3)/(8^3 * sin(theta/2.0))}
theta <- uniroot(thetafun, interval = c(0.00001, 2*pi))$root
yc <- (d / 2) * (1 - cos(theta / 2))
}
return(yc)
}
Qfull <- function(d = NULL, Sf = NULL, n = NULL, k = NULL) {
Vf <- (k / n) * ((d / 4.0) ^ (2.0/3.0)) * sqrt( Sf )
Qf <- Vf * (0.25 * pi * d^2)
return(Qf)
}
# Calculates n/nf using Eq. 30 from Akgiray, 2005, Can. J. Civ. Eng. 32: 490-499
n_adj_func <- function(yd) {
X <- 1 - yd
n_nf <- 1 - 0.8627*X^5 + 0.4281*X^4 + 0.7626*X^3 - 1.02*X^2 + 0.8057*X
return(max(n_nf, 1.0))
}
# Attaches units to output if specified
units::units_options(allow_mixed = TRUE)
return_fcn2 <- function(x = NULL, units = NULL, ret_units = FALSE) {
if (units == "SI") {
out_units <- c("m^3/s","m/s","m^2","m","m","m","m",1,1,"m",1,1,"m^3/s")
} else {
out_units <- c("ft^3/s","ft/s","ft^2","ft","ft","ft","ft",1,1,"ft",1,1,"ft^3/s")
}
if( ret_units ) {
a <- units::mixed_units(unlist(x), out_units)
#names(a) <- names(x)
x <- a
}
return(x)
}
#' @export
#' @rdname manningc
manningc <- function (Q = NULL, n = NULL, Sf = NULL, y = NULL, d = NULL, y_d = NULL,
n_var = FALSE, units = c("SI", "Eng"), ret_units = FALSE ) {
units <- units
if (length(units) != 1) stop("Incorrect unit system. Specify either SI or Eng.")
#check if any values have class 'units' and change to numeric if necessary
for( i in c("Q", "n", "Sf", "y", "d", "y_d") ) {
v <- get(i)
if(inherits(v, "units")) assign(i, units::drop_units(v))
}
#initial check for missing variables and out of bounds
if(missing(y_d)) {
if (length(c(Q, n, Sf, y, d)) != 4) {
stop("There must be exactly one unknown variable among Q, n, Sf, y, d")
}
if (any(c(Q, n, Sf, y, d) <= 0)) {
stop("Either Q, n, Sf, y, d is <= 0. All of these variables must be positive")
}
if ( ( ! missing(d) ) & ( ! missing(y) ) ) {
if ( y > d ) stop("depth y cannot exceed diameter d.")
}
case <- 1
} else {
if ( (length(c(Q, Sf, n)) != 3) & (! missing(d)) & (! missing(y)) ) {
stop("d, y must be missing when given y_d.")
}
if (any(c(Q, Sf, n, y_d) <= 0)) {
stop("Either Q, n, Sf, y_d is <= 0. All of these variables must be positive")
}
if ( y_d > 1.0 ) {
stop("y_d cannot exceed 1.0")
}
case <- 2
}
if (units == "SI") {
g <- 9.80665 # m / s^2
k <- 1.0
mu <- dvisc(T = 20, units = 'SI')
rho <- dens(T = 20, units = 'SI')
dmax <- 3.5 #m maximum pipe size
} else if (units == "Eng") {
g <- 32.2 #ft / s^2
k <- 1.4859
mu <- dvisc(T = 68, units = 'Eng')
rho <- dens(T = 68, units = 'Eng')
dmax <- 12 #ft
} else if (all(c("SI", "Eng") %in% units == FALSE) == FALSE) {
stop("Incorrect unit system. Must be SI or Eng")
}
#########CASE 1########################
if (case == 1) {
if (missing(Q)) {
nf <- n
if (n_var == TRUE) n <- n * n_adj_func(y/d)
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
theta <- 2 * acos(1 - (2 * (y / d)))
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
Qfun <- function(Q) {Q - ((((theta - sin(theta)) * (d ^ 2 / 8)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((theta * d) / 2) ^ (2 / 3))}
Quse <- uniroot(Qfun, interval = c(0.0000001, 200), extendInt = "yes")
Q <- Quse$root
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
} else if (missing(n)) {
theta <- 2 * acos(1 - (2 * (y / d)))
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
# In this case, if n_var == TRUE, what is returned is assumed to be the adjusted n value (variable with depth)
nfun <- function(n) {Q - ((((theta - sin(theta)) * (d ^ 2 / 8)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((theta * d) / 2) ^ (2 / 3))}
nuse <- uniroot(nfun, interval = c(0.0000001, 200), extendInt = "yes")
n <- nuse$root
# If n_var == TRUE, convert n back to n_full for the next calculation
if (n_var == TRUE) nf <- n / n_adj_func(y/d)
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
} else if (missing(Sf)) {
nf <- n
if (n_var == TRUE) n <- n * n_adj_func(y/d)
theta <- 2 * acos(1 - (2 * (y / d)))
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
Sffun <- function(Sf) {Q - ((((theta - sin(theta)) * (d ^ 2 / 8)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((theta * d) / 2) ^ (2 / 3))}
Sfuse <- uniroot(Sffun, interval = c(0.0000001, 200), extendInt = "yes")
Sf <- Sfuse$root
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
} else if (missing(y)) {
nf <- n
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
if ( Q > Qf ) {
stop("Flow Q exceeds full flow for the pipe.")
}
if (n_var == TRUE) {
# Find the depth with n = nf, used as initial guess in search for y
rh <- (n * Q) / (k * sqrt(Sf))
thetafun <- function (theta) ((theta - sin(theta)) * (d ^ 2 / 8)) * (((theta - sin(theta)) * (d ^ 2 / 8) / ((theta * d) / 2)) ^ (2 / 3)) - rh
thetause <- uniroot(thetafun, c(0.001, 2*pi), extendInt = "yes")
theta <- thetause$root
y_fixedn <- (d / 2) * (1 - cos(theta / 2))
# Find n that would fill pipe at given Q, S, d -- upper bound for n
#n_full <- k / Q * (pi/4)*d^2 * (d/4)^(2/3) * sqrt(Sf)
#set up 2 unknowns (x) and 2 equations (y)
fnc_trial <- function(x) {
n1 <- x[1]
y1 <- x[2]
X <- 1 - y1/d
theta <- 2 * acos(1 - (2 * (y1 / d)))
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
R <- A / P
V1 <- (k / n1) * (R^(2.0/3.0))*sqrt(Sf)
Q1 <- V1 * A
y <- numeric(length(x))
# n/nf for estimated for the depth equals the n/nf
y[1] <- 1 - 0.8627*X^5 + 0.4281*X^4 + 0.7626*X^3 - 1.02*X^2 + 0.8057*X - n1/nf
# flow equals the calculated flow given n and
y[2] <- Q - Q1
return(y)
}
#provide initial guesses for unknowns and run the fsolve command
xstart <- c(nf, y_fixedn)
z <- pracma::fsolve(fnc_trial, xstart)
n <- z$x[1]
y <- z$x[2]
if (y > d) stop("feasible solution not found\n")
theta <- 2 * acos(1 - (2 * (y / d)))
}
else {
rh <- (n * Q) / (k * sqrt(Sf))
thetafun <- function (theta) ((theta - sin(theta)) * (d ^ 2 / 8)) * (((theta - sin(theta)) * (d ^ 2 / 8) / ((theta * d) / 2)) ^ (2 / 3)) - rh
thetause <- uniroot(thetafun, c(0.001, 2*pi), extendInt = "yes")
theta <- thetause$root
y <- (d / 2) * (1 - cos(theta / 2))
}
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
}
}
#########CASE 2########################
if (case == 2) {
nf <- n
if (n_var == TRUE) n <- n * n_adj_func(y_d)
theta <- 2 * acos(1 - (2 * (y_d)))
rh <- (n * Q) / (k * sqrt(Sf))
dfun <- function (d) ((theta - sin(theta)) * (d ^ 2 / 8)) * (((theta - sin(theta)) * (d ^ 2 / 8) / ((theta * d) / 2)) ^ (2 / 3)) - rh
duse <- uniroot(dfun, interval = c(0.001, dmax), extendInt = "yes")
d <- duse$root
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
y <- y_d * d
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
}
}
Try the hydraulics package in your browser
Any scripts or data that you put into this service are public.
hydraulics documentation built on June 22, 2024, 7:38 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions return_fcn2 n_adj_func Qfull ycfunc
#' Solves the Manning Equation for gravity flow in a circular pipe
#'
#' This function solves the Manning equation for water flow in a circular pipe
#' at less than full. Uniform flow conditions are assumed, so that the pipe slope
#' is equal to the slope of the water surface and the energy grade line. This is a
#' modification of the code prepared by Irucka Embry in his iemisc package.
#' The iemisc::manningcirc function was adapted here for more limited cases
#' commonly used in classroom exercises, additional checks were included to ensure
#' the pipe is flowing less than full, and a cross-section figure is also
#' available. The iemisc::manningcirc and iemisc::manningcircy functions were
#' combined into a single function. Manning n supplied is assumed to be that for
#' full pipe flow; an optional argument may be supplied to account for n varying
#' with depth.
#'
#' The possible applications of this function for solving the Manning equation
#' for circular pipes are:
#' \tabular{ll}{
#' \strong{Given} \tab \strong{Solve for} \cr
#' y_d, Q, Sf, n \tab d \cr
#' d, Sf, Q, n \tab y \cr
#' y, d, Q, n \tab Sf \cr
#' y, d, Sf, n \tab Q \cr
#' d, Q, Sf, y \tab n
#' }
#'
#' @param Q numeric vector that contains the flow rate [\eqn{m^3 s^{-1}}{m^3/s} or \eqn{ft^3 s^{-1}}{ft^3/s}]
#' @param d numeric vector that contains the pipe diameter [\eqn{m}{m} or \eqn{ft}{ft}]
#' @param Sf numeric vector that contains the slope of the pipe [unitless]
#' @param y numeric vector that contains the water depth [\eqn{m}{m} or \eqn{ft}{ft}]
#' @param y_d numeric vector that contains the ratio of depth to diameter [unitless]
#' @param n numeric vector that contains the Manning roughness coefficient (for full flow or fixed).
#' @param n_var If set to TRUE the value of n for full flow is adjusted with depth. [Default is FALSE]
#' @param units character vector that contains the system of units [options are
#' \code{SI} for International System of Units and \code{Eng} for English (US customary)
#' units. This is used for compatibility with iemisc package]
#' @param ret_units If set to TRUE the value(s) returned are of class \code{units} with
#' units attached to the value. [Default is FALSE]
#'
#' @return Returns a list including the missing parameter:
#' \itemize{
#' \item Q - flow rate
#' \item V - flow velocity
#' \item A - cross-sectional area of flow
#' \item P - wetted perimeter
#' \item R - hydraulic radius (A/P)
#' \item y - flow depth
#' \item d - pipe diameter
#' \item Sf - slope
#' \item n - Manning's roughness (for full flow, or as adjusted if n_var is TRUE)
#' \item yc - critical depth
#' \item Fr - Froude number
#' \item Re - Reynolds number
#' \item Qf - Full pipe flow rate
#' }
#'
#' @author Ed Maurer, Irucka Embry
#'
#' @details The Manning equation (also known as the Strickler equation) describes flow conditions in
#' an open channel under uniform flow conditions. It is often expressed as:
#' \deqn{Q = A\frac{C}{n}{R}^{\frac{2}{3}}{S_f}^{\frac{1}{2}}}
#' where \eqn{C} is 1.0 for SI units and 1.49 for Eng (U.S. Customary) units. Critical depth is
#' defined by the relation (at critical conditions):
#' \deqn{\frac{Q^{2}B}{g\,A^{3}}=1}{Q^2B/gA^3=1}
#' where \eqn{B}{B} is the top width of the water surface. Since B equals zero for a full pipe, critical
#' depth is set to the pipe diameter \eqn{d}{d} if the flow \eqn{Q}{Q} exceeds a value that would produce a
#' critical flow at \eqn{\frac{y}{d}=0.99}{y/d=0.99}.
#'
#' @examples
#'
#' #Solving for flow rate, Q: SI Units
#' manningc(d = 0.6, n = 0.013, Sf = 1./400., y = 0.24, units = "SI")
#' #returns 0.1 m3/s
#'
#' #Solving for Sf, if d=600 mm and pipe is to flow half full
#' manningc(d = 0.6, Q = 0.17, n = 0.013, y = 0.3, units = "SI")
#' #returns required slope of 0.003
#'
#' #Solving for diameter, d when given y_d): Eng (US) units
#' manningc(Q = 83.5, n = 0.015, Sf = 0.0002, y_d = 0.9, units = "Eng")
#' #returns 7.0 ft required diameter
#'
#' #Solving for depth, y when given Q: SI units
#' manningc(Q=0.01, n=0.013, Sf=0.001, d = 0.2, units="SI")
#' #returns depth y = 0.158 m, critical depth, yc = 0.085 m
#'
#' #Solving for depth, y when given Q: SI units, and n variable with depth
#' manningc(Q=0.01, n=0.013, Sf=0.001, d = 0.2, n_var = TRUE, units="SI")
#' #returns depth y = 0.174 m, critical depth, yc = 0.085 m
#'
#' @seealso \code{\link{xc_circle}} for a cross-section diagram of the circular channel
#'
#' @importFrom pracma fsolve
#'
#' @name manningc
NULL
# Function for finding critical depth
ycfunc <- function(Q = NULL, d = NULL, g = NULL) {
# Find theta for 99% full (y/D = 0.99)
thetafull <- 2 * acos(1 - (2 * (0.99)))
Qcfull <- sqrt(g * d^5 * ((thetafull - sin(thetafull))^3)/(8^3 * sin(thetafull/2.0)))
if ( Q >= Qcfull ) {
yc <- d
} else {
thetafun <- function(theta) {Q^2 / ( g * d^5 ) - ((theta - sin(theta))^3)/(8^3 * sin(theta/2.0))}
theta <- uniroot(thetafun, interval = c(0.00001, 2*pi))$root
yc <- (d / 2) * (1 - cos(theta / 2))
}
return(yc)
}
Qfull <- function(d = NULL, Sf = NULL, n = NULL, k = NULL) {
Vf <- (k / n) * ((d / 4.0) ^ (2.0/3.0)) * sqrt( Sf )
Qf <- Vf * (0.25 * pi * d^2)
return(Qf)
}
# Calculates n/nf using Eq. 30 from Akgiray, 2005, Can. J. Civ. Eng. 32: 490-499
n_adj_func <- function(yd) {
X <- 1 - yd
n_nf <- 1 - 0.8627*X^5 + 0.4281*X^4 + 0.7626*X^3 - 1.02*X^2 + 0.8057*X
return(max(n_nf, 1.0))
}
# Attaches units to output if specified
units::units_options(allow_mixed = TRUE)
return_fcn2 <- function(x = NULL, units = NULL, ret_units = FALSE) {
if (units == "SI") {
out_units <- c("m^3/s","m/s","m^2","m","m","m","m",1,1,"m",1,1,"m^3/s")
} else {
out_units <- c("ft^3/s","ft/s","ft^2","ft","ft","ft","ft",1,1,"ft",1,1,"ft^3/s")
}
if( ret_units ) {
a <- units::mixed_units(unlist(x), out_units)
#names(a) <- names(x)
x <- a
}
return(x)
}
#' @export
#' @rdname manningc
manningc <- function (Q = NULL, n = NULL, Sf = NULL, y = NULL, d = NULL, y_d = NULL,
n_var = FALSE, units = c("SI", "Eng"), ret_units = FALSE ) {
units <- units
if (length(units) != 1) stop("Incorrect unit system. Specify either SI or Eng.")
#check if any values have class 'units' and change to numeric if necessary
for( i in c("Q", "n", "Sf", "y", "d", "y_d") ) {
v <- get(i)
if(inherits(v, "units")) assign(i, units::drop_units(v))
}
#initial check for missing variables and out of bounds
if(missing(y_d)) {
if (length(c(Q, n, Sf, y, d)) != 4) {
stop("There must be exactly one unknown variable among Q, n, Sf, y, d")
}
if (any(c(Q, n, Sf, y, d) <= 0)) {
stop("Either Q, n, Sf, y, d is <= 0. All of these variables must be positive")
}
if ( ( ! missing(d) ) & ( ! missing(y) ) ) {
if ( y > d ) stop("depth y cannot exceed diameter d.")
}
case <- 1
} else {
if ( (length(c(Q, Sf, n)) != 3) & (! missing(d)) & (! missing(y)) ) {
stop("d, y must be missing when given y_d.")
}
if (any(c(Q, Sf, n, y_d) <= 0)) {
stop("Either Q, n, Sf, y_d is <= 0. All of these variables must be positive")
}
if ( y_d > 1.0 ) {
stop("y_d cannot exceed 1.0")
}
case <- 2
}
if (units == "SI") {
g <- 9.80665 # m / s^2
k <- 1.0
mu <- dvisc(T = 20, units = 'SI')
rho <- dens(T = 20, units = 'SI')
dmax <- 3.5 #m maximum pipe size
} else if (units == "Eng") {
g <- 32.2 #ft / s^2
k <- 1.4859
mu <- dvisc(T = 68, units = 'Eng')
rho <- dens(T = 68, units = 'Eng')
dmax <- 12 #ft
} else if (all(c("SI", "Eng") %in% units == FALSE) == FALSE) {
stop("Incorrect unit system. Must be SI or Eng")
}
#########CASE 1########################
if (case == 1) {
if (missing(Q)) {
nf <- n
if (n_var == TRUE) n <- n * n_adj_func(y/d)
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
theta <- 2 * acos(1 - (2 * (y / d)))
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
Qfun <- function(Q) {Q - ((((theta - sin(theta)) * (d ^ 2 / 8)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((theta * d) / 2) ^ (2 / 3))}
Quse <- uniroot(Qfun, interval = c(0.0000001, 200), extendInt = "yes")
Q <- Quse$root
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
} else if (missing(n)) {
theta <- 2 * acos(1 - (2 * (y / d)))
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
# In this case, if n_var == TRUE, what is returned is assumed to be the adjusted n value (variable with depth)
nfun <- function(n) {Q - ((((theta - sin(theta)) * (d ^ 2 / 8)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((theta * d) / 2) ^ (2 / 3))}
nuse <- uniroot(nfun, interval = c(0.0000001, 200), extendInt = "yes")
n <- nuse$root
# If n_var == TRUE, convert n back to n_full for the next calculation
if (n_var == TRUE) nf <- n / n_adj_func(y/d)
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
} else if (missing(Sf)) {
nf <- n
if (n_var == TRUE) n <- n * n_adj_func(y/d)
theta <- 2 * acos(1 - (2 * (y / d)))
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
Sffun <- function(Sf) {Q - ((((theta - sin(theta)) * (d ^ 2 / 8)) ^ (5 / 3) * sqrt(Sf)) * (k / n) / ((theta * d) / 2) ^ (2 / 3))}
Sfuse <- uniroot(Sffun, interval = c(0.0000001, 200), extendInt = "yes")
Sf <- Sfuse$root
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
} else if (missing(y)) {
nf <- n
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
if ( Q > Qf ) {
stop("Flow Q exceeds full flow for the pipe.")
}
if (n_var == TRUE) {
# Find the depth with n = nf, used as initial guess in search for y
rh <- (n * Q) / (k * sqrt(Sf))
thetafun <- function (theta) ((theta - sin(theta)) * (d ^ 2 / 8)) * (((theta - sin(theta)) * (d ^ 2 / 8) / ((theta * d) / 2)) ^ (2 / 3)) - rh
thetause <- uniroot(thetafun, c(0.001, 2*pi), extendInt = "yes")
theta <- thetause$root
y_fixedn <- (d / 2) * (1 - cos(theta / 2))
# Find n that would fill pipe at given Q, S, d -- upper bound for n
#n_full <- k / Q * (pi/4)*d^2 * (d/4)^(2/3) * sqrt(Sf)
#set up 2 unknowns (x) and 2 equations (y)
fnc_trial <- function(x) {
n1 <- x[1]
y1 <- x[2]
X <- 1 - y1/d
theta <- 2 * acos(1 - (2 * (y1 / d)))
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
R <- A / P
V1 <- (k / n1) * (R^(2.0/3.0))*sqrt(Sf)
Q1 <- V1 * A
y <- numeric(length(x))
# n/nf for estimated for the depth equals the n/nf
y[1] <- 1 - 0.8627*X^5 + 0.4281*X^4 + 0.7626*X^3 - 1.02*X^2 + 0.8057*X - n1/nf
# flow equals the calculated flow given n and
y[2] <- Q - Q1
return(y)
}
#provide initial guesses for unknowns and run the fsolve command
xstart <- c(nf, y_fixedn)
z <- pracma::fsolve(fnc_trial, xstart)
n <- z$x[1]
y <- z$x[2]
if (y > d) stop("feasible solution not found\n")
theta <- 2 * acos(1 - (2 * (y / d)))
}
else {
rh <- (n * Q) / (k * sqrt(Sf))
thetafun <- function (theta) ((theta - sin(theta)) * (d ^ 2 / 8)) * (((theta - sin(theta)) * (d ^ 2 / 8) / ((theta * d) / 2)) ^ (2 / 3)) - rh
thetause <- uniroot(thetafun, c(0.001, 2*pi), extendInt = "yes")
theta <- thetause$root
y <- (d / 2) * (1 - cos(theta / 2))
}
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
}
}
#########CASE 2########################
if (case == 2) {
nf <- n
if (n_var == TRUE) n <- n * n_adj_func(y_d)
theta <- 2 * acos(1 - (2 * (y_d)))
rh <- (n * Q) / (k * sqrt(Sf))
dfun <- function (d) ((theta - sin(theta)) * (d ^ 2 / 8)) * (((theta - sin(theta)) * (d ^ 2 / 8) / ((theta * d) / 2)) ^ (2 / 3)) - rh
duse <- uniroot(dfun, interval = c(0.001, dmax), extendInt = "yes")
d <- duse$root
Qf <- Qfull(d = d, Sf = Sf, n = nf, k = k)
y <- y_d * d
A <- (theta - sin(theta)) * (d ^ 2 / 8)
P <- ((theta * d) / 2)
B <- d * sin(theta / 2)
R <- A / P
D <- A / B
V <- Q / A
Re <- (rho * R * V) / mu
if (Re < 2000) {
message(sprintf("Low Reynolds number: %.0f indicates not rough turbulent, Manning eq. not valid\n",Re))
}
Fr <- V / (sqrt(g * D))
yc <- ycfunc(Q = Q, d = d, g = g)
out <- list(Q = Q, V = V, A = A, P = P, R = R, y = y, d = d, Sf = Sf, n = n, yc = yc, Fr = Fr, Re = Re, Qf = Qf)
return(return_fcn2(x = out, units = units, ret_units = ret_units))
}
}
Try the hydraulics package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.