System analysis: Irreversible change

David Kneis (firstname.lastname @ tu-dresden.de)

Intro

  • Many systems have a tendency to return to a well-known steady state after perturbation.

  • But this is not necessarily so. Non-linear systems sometimes change dramatically and irreversibly in response to even minor perturbations.

  • Dramatic changes like system collapses are typically the result of positive feedback.

About the example

  • Reservoirs are equipped with special outlet structures to safely handle extreme inflow rates.

  • These structures are typically designed for floods with a return period of 1000 to 10000 years.

  • What if …

About the example

Figure 1: Hirakud dam, Mahanadi River, India.

Load the model

  • The model is part of the rodeoEasy package. You may open the respective workbook on your local computer.
rm(list=ls())                                     # initial clean-up

#remotes::install_github("dkneis/rodeoEasy")      # if not installed yet
library("rodeoEasy")

workbook <- system.file("models/damBreak.xlsx",   # contains the model
  package="rodeoEasy")

m <- build(workbook)                              # build model object

State variables

name unit description default
V m3 reservoir volume 9500000
H_crest m above lake bottom height of dam 10

Model equations

name unit description rate V H_crest
inflow m³ / s inflow to reservoir q_in 1 0
flow_weir m³ / s controlled outflow over the weir structure overflow(mu, w_weir, V/area - h_weir) -1 0
flow_crest m³ / s uncontrolled flow over the dam’s crest overflow(mu, w_crest, V/area - H_crest) -1 0
erosion m / s erosion of the dam in response to uncontrolled overflow ke * max(0, V/area - H_crest) * H_crest 0 -1

Functions

name code
max # intrinsic function
overflow # classical formula to compute the flow rate over a weir
overflow # mu: coefficient for energy dissipation, depends on weir shape
overflow # width: width of the weir crest (m)
overflow # flow depth over the weir crest (m)
overflow overflow <- function (mu, width, depth) {
overflow 2/3 * mu * sqrt(2 * 9.81) * width * max(0, depth)^(3/2)
overflow }

Parameters

name unit description default
area m2 lake surface area 1e+06
mu - weir overflow coefficient 7e-01
h_weir m above lake bottom height of weir base 9e+00
w_weir m width of weir 5e+00
w_crest m width of dam crest 1e+02
q_in m3 / s inflow rate 5e+00
ke 1 / m / s erosion rate constant 1e-02

Predictions for different inflows

scen <- list(                             # inflow scenarios
  low=      c(q_in=  1),
  moderate= c(q_in= 10),
  high=     c(q_in= 20)
)

x <- run.scenarios(m,                     # simulations
  times=seq(0, 24, .1) * 3600,
  scenarios=scen, plot.vars=F)

x[,"time"] <- x[,"time"] / 3600           # let time be in hours
x <- cbind(x, total_outflow=with(x,       # add total outflow rate
  flow_weir + flow_crest))

Simulation outcomes

oma <- par("mar")                         
par(mar=c(4,4,1,1), mfrow=c(2,2), cex=1.2)

key <- show.scenarios(x, "V", xlab="time (hours)")
key <- show.scenarios(x, "H_crest")
key <- show.scenarios(x, "total_outflow", ylog=TRUE)

plot(0, 0, type="n", ann=F, axes=F)
legend("center", bty="n", lty=1, col=key,
  legend=names(key), title="Inflow rate")

par(mfrow=c(1,1), cex=1, mar=oma)

Simulation outcomes