-
Notifications
You must be signed in to change notification settings - Fork 10
/
plumerise.R
175 lines (170 loc) · 6.58 KB
/
plumerise.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
#' Calculate plume rise height.
#'
#' @description Calculate the maximum height of rise based on Brigs (1975), the height is calculated using different formulations depending on stability and wind conditions.
#'
#' @format data.frame with the input, rise (m) and effective higt (m)
#'
#' @param df data.frame with micrometeorological and emission data
#' @param imax maximum number of iteractions
#' @param ermax maximum error
#' @param Hmax use weil limit for plume rise, see details
#' @param verbose display additional information
#'
#' @return a data.frame with effective height of emissions for pointSource function
#'
#' @export
#'
#' @references The plume rise formulas are from Brigs (1975):"Brigs, G. A. Plume rise predictions, Lectures on Air Pollution and Environmental Impact Analyses. Amer. Meteor. Soc. p. 59-111, 1975." and Arya 1999: "Arya, S.P., 1999, Air Pollution Meteorology and Dispersion, Oxford University Press, New York, 310 p."
#'
#' The limits are from Weil (1979): "WEIL, J.C. Assessmet of plume rise and dispersion models using LIDAR data, PPSP-MP-24. Prepared by Environmental Center, Martin Marietta Corporation, for Maryland Department of Natural Resources. 1979."
#'
#' The example is data from a chimney of the Candiota thermoelectric powerplant from Arabage et al (2006) "Arabage, M. C.; Degrazia, G. A.; Moraes O. L. Simulação euleriana da dispersão local da pluma de poluente atmosférico de Candiota-RS. Revista Brasileira de Meteorologia, v.21, n.2, p. 153-160, 2006."
#'
#' @details The input data.frame must contains the folloging colluns:
#'
#' - z: height of the emission (m)
#'
#' - r: source raius (m)
#'
#' - Ve: emission velocity (m/s)
#'
#' - Te: emission temperature (K)
#'
#' - ws: wind speed (m/s)
#'
#' - Temp: ambient temperature (K)
#'
#' - h: height of the Atmospheric Boundary Layer-ABL (m)
#'
#' - L: Monin-Obuhkov Lench (m)
#'
#' - dtdz: lapse ration of potential temperature, used only for stable ABL (K/m)
#'
#' - Ustar: atriction velocity, used only for neutral ABL (m/s)
#'
#' - Wstar: scale of convectie velocity, used only for convective ABL (m/s)
#'
#' Addcitionaly some combination of ws, Wstar and Ustar can produce inacurate results, Weil (1979) propose a geometric limit of 0.62 * (h - Hs) for the rise value.
#'
#' @examples
#' candiota <- matrix(c(150,1,20,420,3.11,273.15 + 3.16,200,-34.86,3.11,0.33,
#' 150,1,20,420,3.81,273.15 + 4.69,300,-34.83,3.81,0.40,
#' 150,1,20,420,3.23,273.15 + 5.53,400,-24.43,3.23,0.48,
#' 150,1,20,420,3.47,273.15 + 6.41,500,-15.15,3.48,0.52,
#' 150,1,20,420,3.37,273.15 + 6.35,600, -8.85,3.37,2.30,
#' 150,1,20,420,3.69,273.15 + 5.93,800,-10.08,3.69,2.80,
#' 150,1,20,420,3.59,273.15 + 6.08,800, -7.23,3.49,1.57,
#' 150,1,20,420,4.14,273.15 + 6.53,900,-28.12,4.14,0.97),
#' ncol = 10, byrow = TRUE)
#' candiota <- data.frame(candiota)
#' names(candiota) <- c("z","r","Ve","Te","ws","Temp","h","L","Ustar","Wstar")
#' row.names(candiota) <- c("08:00","09:00",paste(10:15,":00",sep=""))
#' candiota <- plumeRise(candiota,Hmax = TRUE)
#' print(candiota)
#'
plumeRise <- function(df, imax = 10, ermax = 1/100, Hmax = TRUE, verbose = TRUE)
{
if(imax < 2) imax = 2
g <- 9.81 # m / s2
rise <- rep(0,nrow(df))
He <- rep(0,nrow(df))
for(j in 1:nrow(df)){
# source parameters
Hs <- df$z[j] # source height
r <- df$r[j] # source radius
Vi <- df$Ve[j] # exhaust velocity
Ti <- df$Te[j] # exhaust temperature
# micrometeorological data
U <- df$ws[j] # mean wind speed
Ta <- df$Temp[j] # temperature
h <- df$h[j] # atmospheric boundary layer height
L <- df$L[j] # monin-obukov length
# boyance term
Flu <- g * Vi * r^2 * abs(Ti - Ta)/Ti
i <- 1
# strong convection
if(h/abs(L) > 10 & L < 0){
if(verbose)
cat(paste("strong convective, h/L =",h/L,"\n"))
Wstar<- df$Wstar[j] # w*
deltaH <- 4.3 * (Flu / (U * Wstar^2))^(3/5) * h^(2/5)
}else
# slytly convective
if(h/abs(L) <= 10 & h/abs(L) > 1 & L < 0){
if(verbose)
cat(paste("convective, h/L =",h/L,"\n"))
Wstar<- df$Wstar[j] # w*
Wd <- 0.4 * Wstar # downdrafts mean speed
a <- ( Flu / U*Wd^2 )^(3/5)
b <- 2 * Hs
x <- Hs / 2
err <- Inf
while(err > ermax & i < imax){
i <- i+1
f <- x^3 - a*(x^2) -2*a*b*x -a*(b^2)
fl <- 3*(x^2) -2*a*x -2*a*b
old <- x
deltaH <- x - (f/fl)
err <- abs( (deltaH - old )/old )
x <- deltaH
# if(verbose){
# print(paste("Dh old=",old))
# print(paste("err=",err))
# print("=====================================================")
# }
}
}else
# neutral
if(abs(h/L) <= 1.0 & abs(h/L) >= 0.0){ # & L < 0.0
if(verbose)
cat(paste("neutral, h/L =",h/L,"\n"))
Ustar<- df$Ustar[j] # u*
a <- ((1.3 * Flu) /(U * Ustar))^(3)
x <- Hs / 4
err <- Inf
while(err > ermax & i< imax){
i <- i+1
f <- x^5 -a*x^2 -2*a*x*Hs -2*a*Hs^2
fl <- 5*x^4 -2*a*x -2*a*Hs
old <- x
deltaH <- x - (f/fl)
err <- abs( (deltaH - old )/old )
x <- deltaH
# if(verbose){
# print(paste("Dh old=",old))
# print(paste("err=",err))
# print("=====================================================")
# }
}
}else
# stable
if(h/L > 1.0){
dtdz <- df$dtdz[j] # temperature vertical gradient
s <- (g/Ta) * dtdz
if(U <= 1){
if(verbose)
cat(paste("stable, h/L =",h/L,"- calm,","U=",U,"m/s","\n"))
deltaH <- 5 * Flu^(1/4) * s^(3/5)
}else{
if(verbose)
cat(paste("stable, h/L =",h/L,"- windy,","U=",U,"m/s","\n"))
deltaH <- 2.6 * (Flu/(U*s))^(1/3)
}
}
if(i == imax)
cat("* max iterations reached!","\n")
# Weil (1979) limit for plume rises
if(Hmax){
weil <- 0.62 * (h - Hs)
rise[j] = min(deltaH,weil)
if(verbose & weil < deltaH)
cat(paste("using weil max=",weil,"\n"))
}else{
rise[j] = deltaH
}
He[j] <- rise[j] + Hs
}
df <- cbind(df,rise)
df <- cbind(df,He)
return(df)
}