#' One-lump Transient Heat Budget (for use with deSolve package)
#'
#' Transient, 'one-lump', heat budget for computing rate of change of temperature
#' under environmental conditions that vary with time, using interpolation functions to
#' estimate environmental conditions at particular time intervals.
#' Michael Kearney, Raymond Huey and Warren Porter developed this R function and example in September 2017.
#' @param t time intervals (s) at which output is required
#' @param Tc_init = 5, initial temperature (°C) Organism shape, 0-5, Determines whether standard or custom shapes/surface area/volume relationships are used: 0=plate, 1=cyl, 2=ellips, 3=lizard (desert iguana), 4=frog (leopard frog), 5=custom (see details)
#' @param Ww_g = 500 weight (g)
#' @param rho_body = 1000, animal density (kg/m3)
#' @param q = 0, metabolic rate (W/m3)
#' @param c_body = 3073, specific heat of flesh (J/kg-C)
#' @param k_flesh = 0.5, conductivity of flesh (W/mK)
#' @param emis = 0.95, emissivity of skin (-)
#' @param alpha = 0.85, animal solar absorptivity (-)
#' @param geom = 2, organism shape, 0-5, Determines whether standard or custom shapes/surface area/volume relationships are used: 0=plate, 1=cyl, 2=ellips, 3=lizard (desert iguana), 4=frog (leopard frog), 5=custom (see parameter 'shape_coeffs')
#' @param shape_b = 1/5, proportionality factor (-) for going from volume to area, represents ratio of width:height for a plate, length:diameter for cylinder, b axis:a axis for ellipsoid
#' @param shape_c = 1/5, proportionality factor (-) for going from volume to area, represents ratio of length:height for a plate, c axis:a axis for ellipsoid
#' @param shape_coefs = c(10.4713,.688,0.425,0.85,3.798,.683,0.694,.743), Custom surface area coefficients. Operates if posture = 5, and consists of 4 pairs of values representing the parameters a and b of a relationship AREA=a*Ww_g^b, where AREA is in cm2 and Ww_g is in g. The first pair are a and b for total surface area, then a and b for ventral area, then for silhouette area normal to the sun, then silhouette area perpendicular to the sun
#' @param posture = 'n', pointing normal 'n' or parallel 'p' to the sun's rays, or average 'a'?
#' @param orient = 1, does the object orient toward the sun? (0,1)
#' @param fatosk = 0.4, solar configuration factor to sky (-)
#' @param fatosb = 0.4, solar configuration factor to substrate (-)
#' @param dyn_Q = 0, dynamic metabolic heat generation as a function of temperature, based on approxfun called qf (1) or constant based on parameter q (0)
#' @param alpha_sub = 0.2, substrate solar reflectivity, decimal percent
#' @param pdif = 0.1, proportion of solar energy that is diffuse (rather than direct beam)
#' @param fluid = 0, fluid type, air (0) or water (1)
#' @param Tairf air temperature function with time, generated by 'approxfun' (°C)
#' @param Tradf radiant temperature function with time, generated by 'approxfun'(°C), averaging ground and sky
#' @param velf wind speed function with time, generated by 'approxfun' (m/s)
#' @param Qsolf radiation function with time, generated by 'approxfun' (W/m2)
#' @param Zenf zenith angle of sun function with time, generated by 'approxfun' (90 is below horizon), degrees
#' @param press air pressure (Pa)
#' @return Tc Core temperature (°C)
#' @return Tcf Final (steady state) temperature (°C), if conditions remained constant indefinitely
#' @return tau Time constant (s)
#' @return dTc Rate of change of core temperature (°C/s)
#' @usage ode(y = Tb_init, t = times, func = onelump_var, parms = indata)
#' @examples
#' library(deSolve) # note due to some kind of bug in deSolve, it must be loaded before NicheMapR!
#' library(NicheMapR)
#'
#' # get microclimate data
#' loc <- c(133.8779, -23.6987) # Alice Springs, Australia
#' Usrhyt <- 0.05 # height of midpoint of animal, m
#' micro <- micro_global(loc = loc, Usrhyt = Usrhyt) # run the model with specified location and animal height
#' metout <- as.data.frame(micro$metout) # above ground microclimatic conditions, min shade
#' soil <- as.data.frame(micro$soil) # soil temperatures, minimum shade
#'
#' # append dummy dates
#' days <- rep(seq(1, 12), 24)
#' days <- days[order(days)]
#' dates <- days + metout$TIME / 60 / 24 - 1 # dates for hourly output
#' dates2 <- seq(1, 12, 1) # dates for daily output
#' metout <- cbind(dates, metout)
#' soil <- cbind(dates, soil)
#'
#' # combine relevant input fields
#' microclimate <- cbind(metout[, 1:5], metout[, 8], soil[, 4], metout[, 13:15], metout[, 6])
#' colnames(microclimate) <- c('dates', 'DOY', 'TIME', 'TALOC', 'TA1.2m', 'VLOC', 'TS', 'ZEN', 'SOLR', 'TSKYC', 'RHLOC')
#'
#' # define animal parameters - here simulating a 1000 g ellipsoid
#' c_body <- 3342 # specific heat of flesh, J/kg-C
#' rho_body <- 1000 # animal density, kg/m3
#' q <- 0 # metabolic rate, W/m3
#' k_flesh <- 0.5 # thermal conductivity of flesh, W/mK
#' geom <- 2 # shape, -
#' posture <- 'n' # pointing normal 'n' or parallel 'p' to the sun's rays, or average 'a'?
#' orient <- 1 # does the object orient toward the sun? (0,1)
#' shape_b <- 1/5 # shape coefficient a, -
#' shape_c <- 1/5 # shape coefficient b, -
#' shape_coefs <- c(10.4713, 0.688, 0.425, 0.85, 3.798, 0.683, 0.694, 0.743)
#' fatosk <- 0.4 # solar configuration factor to sky, -
#' fatosb <- 0.4 # solar configuration factor to substrate, -
#' alpha <- 0.9 # animal solar absorptivity, -
#' emis <- 0.95 # emissivity of skin, -
#' Ww_g <- 1000 # weight, g
#' alpha_sub <- 0.8 # substrate solar absorptivity, -
#' press <- 101325 # air pressure, Pa
#' pdif <- 0.1 # proportion of solar energy that is diffuse, -
#' dyn_q <- 0 # not dynamically varying q, -
#' fluid <- 0 # air, -
#'
#' # loop through middle day of each month
#' DOYs = c(15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349)
#' mons = c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December")
#'
#' for (i in 1:length(DOYs)) {
#' simday = DOYs[i]
#' microclim <- subset(microclimate, microclimate$DOY == simday)
#'
#' # use approxfun to create interpolations for the required environmental variables
#' time <- seq(0, 60 * 24, 60) #60 minute intervals from microclimate output
#' time <- time * 60 # minutes to seconds
#' hours <- time/3600 # seconds to hours
#' Qsolf <- approxfun(time, c(microclim[, 9], (microclim[1, 9] + microclim[24, 9]) /
#' 2), rule = 2)
#' # approximate radiant temperature as the average of sky and substrate temperature
#' Tradf <- approxfun(time, rowMeans(cbind(c(microclim[, 7], (microclim[1, 7] + microclim[24, 7]) / 24), c(microclim[, 10], (microclim[1, 10] + microclim[24, 10]) / 24)), na.rm = TRUE), rule = 2)
#' velf <- approxfun(time, c(microclim[, 6], (microclim[1, 6] + microclim[24, 6]) / 2), rule = 2)
#' Tairf <- approxfun(time, c(microclim[, 4], (microclim[1, 4] + microclim[24, 4]) / 2), rule = 2)
#' Zenf <- approxfun(time, c(microclim[, 8], 90), rule = 2)
#'
#' t = seq(1, 3600 * 24, 60) # sequence of times for predictions (1 min intervals)
#' indata<-list(alpha = alpha, emis = emis, alpha_sub = alpha_sub, press = press, Ww_g = Ww_g, c_body = c_body, rho_body = rho_body, q = q, k_flesh = k_flesh, geom = geom, posture = posture, shape_b = shape_b, shape_c = shape_c, shape_coefs = shape_coefs, pdif = pdif, fatosk = fatosk, fatosb = fatosb, fluid = fluid, dyn_q = dyn_q)
#'
#' Tc_init<-Tairf(1) # set inital Tc as air temperature
#'
#' Tbs_ode <- as.data.frame(ode(y = Tc_init, times = time, func = onelump_var, parms = indata))
#' colnames(Tbs_ode) <- c('time', 'Tc', 'Tcf', 'tau', 'dTdt')
#' Tbs_ode$time <- Tbs_ode$time / 3600 # convert to hours
#'
#' with(Tbs_ode, plot(Tc ~ time, type = 'l', col = '1', ylim = c(-10, 80), xlim = c(0, 23), ylab='Temperature, °C',xlab = 'hour of day', main = paste0(round(loc[1], 1), round(loc[2], 1), ", ", mons[i], ", ", Ww_g," g")))
#' with(Tbs_ode, points(Tcf ~ time, type = 'l', col = '2'))
#' points(Tairf(time) ~ hours, type = 'l', col = 'blue', lty = 2)
#' legend(0,70, c("Tc", "Tcf", "Tair"), lty = c(1, 1, 2), lwd = c(2.5, 2.5, 2.5), col = c("black", "red", "blue"))
#' }
#' @export
onelump_var <- function(t, y, indata) {
with(as.list(c(indata, y)), {
sigma <- 5.67e-8 # Stefan-Boltzman, W/(m.K)
Tair <- Tairf(t) # interpolate air temp to current time step
vel <- max(0,velf(t)) # interpolate wind speed to current time step
Qsol <- max(0,Qsolf(t)) # interpolate solar to current time step
Trad <- Tradf(t) # interpolate radiant temp to current time step
Zen <- min(max(0,Zenf(t)),90) # interpolate zenith angle to current time step
Zenith <- Zen * pi / 180 # zenith angle in radians
Tc <- y
Tskin <- Tc + 0.01
vel[vel < 0.01] <- 0.01 # don't let wind speed go too low - always some free convection
S2 <- 0.0001 # shape factor, arbitrary initialization, because not defined except for ellipsoid
if(fluid == 0){
DENSTY <- press / (287.04 * (Tair + 273.15)) # air density, kg/m3
THCOND <- 0.02425 + (7.038 * 10 ^ -5 * Tair) # air thermal conductivity, W/(m.K)
VISDYN <- (1.8325 * 10 ^ -5 * ((296.16 + 120) / ((Tair + 273.15) + 120))) * (((Tair + 273.15) / 296.16) ^ 1.5) # dynamic viscosity of air, kg/(m.s)
DIFVPR <- 2.26e-05 * (((Tair + 273.15) / 273.15) ^ 1.81) * (1e+05 / press)
}else{
WATERPROP.out <- WATERPROP(Tair)
DENSTY <- WATERPROP.out$DENSTY
THCOND <- WATERPROP.out$THCOND
VISDYN <- WATERPROP.out$VISDYN
DIFVPR <- 0
}
VISKIN <- VISDYN / DENSTY
# geometry section ############################################################
m <- Ww_g / 1000 # convert weight to kg
C <- m * c_body # thermal capacitance, J/K
V <- m / rho_body # volume, m3
if(dyn_q == 1){
Q_gen <- qf(Tc) * V # total metabolic heat, J
}else{
Q_gen <- q * V
}
L <- V ^ (1 / 3) # characteristic dimension, m
# FLAT PLATE geometry
if (geom == 0) {
AHEIT <- (V / (shape_b * shape_c)) ^ (1 / 3) # length, m
AWIDTH <- AHEIT * shape_b # width, m
ALENTH <- AHEIT * shape_c # height, m
ATOT <- ALENTH * AWIDTH * 2 + ALENTH * AHEIT * 2 + AWIDTH * AHEIT * 2 # total area, m2
ASILN <- ALENTH * AWIDTH # max silhouette area, m2
ASILP <- AWIDTH * AHEIT # min silhouette area, m2
L <- AHEIT # characteristic dimension, m
if (AWIDTH <= ALENTH) {
L <- AWIDTH
} else{
L <- ALENTH
}
R <- ALENTH / 2 # 'radius', m
}
# CYLINDER geometry
if (geom == 1) {
R1 <- (V / (pi * shape_b * 2)) ^ (1 / 3) # radius, m
ALENTH <- 2 * R1 * shape_b # length, m
ATOT <- 2 * pi * R1 ^ 2 + 2 * pi * R1 * ALENTH # total surface area, m2
AWIDTH <- 2 * R1 # width, m
ASILN <- AWIDTH * ALENTH # max silhouette area, m2
ASILP <- pi * R1 ^ 2 # min silhouette area, m2
L <- ALENTH # characteristic dimension, m
R2 <- L / 2
if (R1 > R2) {
# choose shortest dimension as R
R <- R2
} else{
R <- R1
}
}
# Ellipsoid geometry
if (geom == 2) {
A1 <- ((3 / 4) * V / (pi * shape_b * shape_c)) ^ 0.333 # axis A, m
B1 <- A1 * shape_b # axis B, m
C1 <- A1 * shape_c # axis C, m
P1 <- 1.6075 # a constant
ATOT <- (4 * pi * (((A1 ^ P1 * B1 ^ P1 + A1 ^ P1 * C1 ^ P1 + B1 ^ P1 * C1 ^ P1)) / 3) ^ (1 / P1)) # total surface area, m2
ASILN <- max(pi * A1 * C1, pi * B1 * C1) # max silhouette area, m2
ASILP <- min(pi * A1 * C1, pi * B1 * C1) # min silhouette area, m2
S2 <- (A1 ^ 2 * B1 ^ 2 * C1 ^ 2) / (A1 ^ 2 * B1 ^ 2 + A1 ^ 2 * C1 ^ 2 + B1 ^ 2 * C1 ^ 2) # fraction of semi-major and minor axes, see Porter and Kearney 2009 supp1
#k_flesh <- 0.5 + 6.14 * B1 + 0.439 # thermal conductivity of flesh as a function of radius, see Porter and Kearney 2009
}
# Lizard geometry - DESERT IGUANA (PORTER ET AL. 1973 OECOLOGIA)
if (geom == 3) {
ATOT <- (10.4713 * Ww_g ^ .688) / 10000. # total surface area, m2
AV <- (0.425 * Ww_g ^ .85) / 10000. # ventral surface area, m2
# NORMAL AND POINTING @ SUN SILHOUETTE AREA: PORTER & TRACY 1984
ASILN <- (3.798 * Ww_g ^ .683) / 10000. # Max. silhouette area (normal to the sun), m2
ASILP <- (0.694 * Ww_g ^ .743) / 10000. # Min. silhouette area (pointing toward the sun), m2
R <- L
}
# Frog geometry - LEOPARD FROG (C.R. TRACY 1976 ECOL. MONOG.)
if (geom == 4) {
ATOT <- (12.79 * Ww_g ^ 0.606) / 10000. # total surface area, m2
AV <- (0.425 * Ww_g ^ 0.85) / 10000. # ventral surface area, m2
# NORMAL AND POINTING @ SUN SILHOUETTE AREA: EQ'N 11 TRACY 1976
ZEN <- 0
PCTN <- 1.38171E-06 * ZEN ^ 4 - 1.93335E-04 * ZEN ^ 3 + 4.75761E-03 * ZEN ^ 2 - 0.167912 * ZEN + 45.8228
ASILN <- PCTN * ATOT / 100 # Max. silhouette area (normal to the sun), m2
ZEN <- 90
PCTP <- 1.38171E-06 * ZEN ^ 4 - 1.93335E-04 * ZEN ^ 3 + 4.75761E-03 * ZEN ^ 2 - 0.167912 * ZEN + 45.8228
ASILP <- PCTP * ATOT / 100 # Min. silhouette area (pointing toward the sun), m2
R <- L
}
# user defined geometry
if (geom == 5) {
ATOT <- (shape_coefs[1] * Ww_g ^ shape_coefs[2]) / 10000. # total surface area, m2
AV <- (shape_coefs[3] * Ww_g ^ shape_coefs[4]) / 10000 # ventral surface area, m2
# NORMAL AND POINTING @ SUN SILHOUETTE AREA: PORTER & TRACY 1984
# User must define Max. silhouette area (normal to the sun)
ASILN <- (shape_coefs[5] * Ww_g ^ shape_coefs[6]) / 10000 # Max. silhouette area (normal to the sun), m2
# User must define Min. silhouette area (pointing toward the sun)
ASILP <- (shape_coefs[7] * Ww_g ^ shape_coefs[8]) / 10000 # Min. silhouette area (pointing toward the sun), m2
R <- L
}
# end geometry section ############################################################
if (max(Zen) >= 90) {
Q_norm <- 0
}else{
if(orient == 1){
Q_norm <- (Qsol / cos(Zenith))
}else{
Q_norm <- Qsol
}
}
if (Q_norm > 1367) {
Q_norm <- 1367 #making sure that low sun angles don't lead to solar values greater than the solar constant
}
if (posture == 'p') {
Qabs <- (Q_norm * (1 - pdif) * ASILP + Qsol * pdif * fatosk * ATOT + Qsol * (1 - alpha_sub) * fatosb * ATOT) * alpha
}
if (posture == 'n') {
Qabs <- (Q_norm * (1 - pdif) * ASILN + Qsol * pdif * fatosk * ATOT + Qsol * (1 - alpha_sub) * fatosb * ATOT) * alpha
}
if (posture == 'a') {
Qabs <- (Q_norm * (1 - pdif) * (ASILN + ASILP) / 2 + Qsol * pdif * fatosk * ATOT + Qsol * (1 - alpha_sub) * fatosb * ATOT) * alpha
}
Re <- DENSTY * vel * L / VISDYN # Reynolds number
PR <- 1005.7 * VISDYN / THCOND # Prandlt number
if(fluid == 0){
SC <- VISDYN / (DENSTY * DIFVPR) # Schmidt number
}else{
SC <- 1
}
if (geom == 0) {
NUfor <- 0.102 * Re ^ 0.675 * PR ^ (1. / 3.)
}
if (geom == 3 | geom == 5) {
NUfor <- 0.35 * Re ^ 0.6
}
if (geom == 1) {
# FORCED CONVECTION OF A CYLINDER
# ADJUSTING NU - RE CORRELATION FOR RE NUMBER (P. 260 MCADAMS,1954)
if (Re < 4) {
NUfor = 0.891 * Re ^ 0.33
} else{
if (Re < 40) {
NUfor = 0.821 * Re ^ 0.385
} else{
if (Re < 4000) {
NUfor = 0.615 * Re ^ 0.466
} else{
if (Re < 40000) {
NUfor = 0.174 * Re ^ 0.618
} else{
if (Re < 400000) {
NUfor = 0.0239 * Re ^ 0.805
} else{
NUfor = 0.0239 * Re ^ 0.805
}
}
}
}
}
}
if (geom == 2 | geom == 4) {
NUfor <- 0.35 * Re ^ (0.6) # Nusselt number, forced convection
}
h_conv_forced <- NUfor * THCOND / L # convection coefficent, forced
GR <- abs(DENSTY ^ 2 * (1 / (Tair + 273.15)) * 9.80665 * L ^ 3 * (Tskin - Tair) / VISDYN ^ 2) # Grashof number
Raylei <- GR * PR # Rayleigh number
# get Nusselt for Free Convect
if (geom == 0) {
NUfre = 0.55 * Raylei ^ 0.25
}
if (geom == 1 | geom == 3 | geom == 5) {
if (Raylei < 1.0e-05) {
NUfre = 0.4
} else{
if (Raylei < 0.1) {
NUfre = 0.976 * Raylei ^ 0.0784
} else{
if (Raylei < 100) {
NUfre = 1.1173 * Raylei ^ 0.1344
} else{
if (Raylei < 10000.) {
NUfre = 0.7455 * Raylei ^ 0.2167
} else{
if (Raylei < 1.0e+09) {
NUfre = 0.5168 * Raylei ^ 0.2501
} else{
if (Raylei < 1.0e+12) {
NUfre = 0.5168 * Raylei ^ 0.2501
} else{
NUfre = 0.5168 * Raylei ^ 0.2501
}
}
}
}
}
}
}
if (geom == 2 | geom == 4) {
Raylei = (GR ^ 0.25) * (PR ^ 0.333)
NUfre = 2 + 0.60 * Raylei
}
Nutotal <- (NUfre^3 + NUfor^3)^(1./3.)
sh <- Nutotal * (SC / PR) ^ (1 / 3)
h_conv <- Nutotal*(THCOND/L) # combined convection coefficient
#R_conv <- 1 / (h_conv * ATOT) # convective resistance, eq. 3 of Kearney, Huey and Porter in prep. Appendix 1
h_rad <- 4 * emis * sigma * ((Tc + Trad) / 2 + 273.15) ^ 3 # radiation heat transfer coefficient, eq. 46 of Transient Equations Derivation vignette
if(geom == 2){ # ellipsoid
j <- (Qabs + Q_gen + h_conv * ATOT * ((q * S2) / (2 * k_flesh) + Tair) + h_rad * ATOT * ((q * S2) / (2 * k_flesh) + Trad)) / C #based on eq. 52 of Kearney, Huey and Porter in prep. Appendix 1
}else{ # assume cylinder
j <- (Qabs + Q_gen + h_conv * ATOT * ((q * R ^ 2) / (4 * k_flesh) + Tair) + h_rad * ATOT * ((q * R ^ 2) / (2 * k_flesh) + Trad)) / C #based on eq. 52 of Kearney, Huey and Porter in prep. Appendix 1
}
theta_Tc <- ATOT * (Tc * h_conv + Tc * h_rad) / C #based on eq. 48 of Transient Equations Derivation vignette
theta <- ATOT * (h_conv + h_rad) / C #based on eq. 48 of Transient Equations Derivation vignette
Tcf <- j / theta # final Tc = j/theta, based on eq. 23 of Transient Equations Derivation vignette
Tci <- Tc # initial temperature
Tc <- (Tci - Tcf) * exp(-1 * theta * t) + Tcf # Tc at time t, based on eq. 33 of Transient Equations Derivation vignette
tau <- 1 / theta # time constant
dTc <- j - theta_Tc # rate of temperature change (°C/sec)
list(y = dTc, x = Tcf, z = tau, zz = dTc)
})
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.