R/twolump.R

In mrke/NicheMapR: R implementation of Niche Mapper software for biophysical modelling

#' Two-lump Transient Heat Budget (for use with deSolve package)
#'
#' Transient, 'two-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 core temperature (°C)
#' @param Tsk_init = 5.1, initial shell (skin) temperature (°C)
#' @param Ts_init = 5.2, initial surface temperature (°C)
#' @param Ww_g = 500, animal weight (g)
#' @param rho_body = 932, animal density (kg/m3)
#' @param x_shell = 0.001, shell thickness, m
#' @param q = 0, metabolic heat production rate W/m3
#' @param c_body_inner = 3073, Specific heat of flesh J/(kg-K)
#' @param c_body_outer = 3073, Specific heat of outer shell J/(kg-K)
#' @param k_inner = 0.5, Thermal conductivity of inner shell (W/mK, range: 0.412-2.8)
#' @param k_outer = 0.5, Thermal conductivity of outer shell (W/mK, range: 0.412-2.8)
#' @param emis = 0.95, Emissivity of animal (0-1)
#' @param alpha = 0.85, solar absorptivity, decimal percent
#' @param geom = 2, Organism shape, 1-2, Determines surface area/volume relationships: 1=cyl, 2=ellips
#' @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 posture = 'n', pointing normal 'n', parallel 'p' to the sun's rays, or 'a' in between?
#' @param orient = 1, does the object orient toward the sun? (0,1)
#' @param fatosk = 0.4, Configuration factor to sky for infrared calculations (-)
#' @param fatosb = 0.4, Configuration factor to substrate for infrared calculations (-)
#' @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 = 101325, air pressure (Pa)
#' @return Tc Core temperature (°C)
#' @return Tsk 'Skin' / shell temperature (°C)
#' @return Ts Outer surface temperature (°C)
#' @return Tcf Final (steady state) temperature (°C), if conditions remained constant indefinitely
#' @usage ode(y = c(Tc_init, Tsk_init, Ts_init), times = time, func = twolump, 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
#' c_body_inner <- c_body
#' c_body_outer <- c_body
#' rho_body <- 1000 # animal density, kg/m3
#' q <- 0 # metabolic rate, W/m3
#' k_inner <- 0.5 # thermal conductivity of core, W/mK
#' k_outer <- 0.5 # thermal conductivity of shell, W/mK
#' geom <- 2 # shape, -
#' x_shell <- 0.005 # thickness of outer shell (m)
#' 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, -
#' 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(Ww_g = Ww_g, x_shell = x_shell, geom = geom, k_inner = k_inner, k_outer = k_outer, q = q, c_body_inner = c_body_inner, c_body_outer = c_body_outer, emis = emis, rho_body = rho_body, alpha = alpha, shape_b = shape_b, shape_c = shape_c, posture = posture, orient = orient, fatosk = fatosk, fatosb = fatosb, alpha_sub = alpha_sub, pdif = pdif, press = press, fluid = fluid, dyn_q = dyn_q)
#'
#'   Tc_init<-Tairf(1) # set inital Tc as air temperature
#'   Tsk_init <- Tc_init
#'   Ts_init <- Tc_init
#'
#'   Tbs_ode <- as.data.frame(ode(y = c(Tc_init, Tsk_init, Ts_init), times = time, func = twolump, parms = indata))
#'   colnames(Tbs_ode) = c("time", "Tc", "Tsk", "Ts", "Tcf")
#'   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(loc, ", ", mons[i], ", ", Ww_g," g")))
#'   with(Tbs_ode, points(Ts ~ time, lty = 2, type = 'l', col = '1'))
#'   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", "Ts", "Tcf", "Tair"), lty = c(1, 2, 1, 2), lwd = c(2.5, 2.5, 2.5, 2.5), col = c("black", "black", "red", "blue"))
#' }
#' @export
twolump<-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 <- as.numeric(y[1]) # core temperature, °C
    Tsk <- as.numeric(y[2]) # surface temperature
    Ts <- as.numeric(y[3]) # skin temperature
    vel[vel < 0.01] <- 0.01 # don't let wind speed go too low
    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
    C <- m * c_body_inner # 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
    V_inner <- (L - x_shell) ^ 3
    V_shell <- V - V_inner
    C_sk <- V_shell * rho_body * c_body_outer
    C_c <- V_inner * rho_body * c_body_inner

    # CYLINDER geometry
    if(geom == 1){
      R1 <- (V / (pi * shape_b * 2)) ^ (1 / 3) # radius, m
      ALENTH <- 2 * R1 * shape_b # length, m
      V_inner <- pi * (R1 - x_shell) ^ 2 * (ALENTH - x_shell)
      V_shell <- V - V_inner
      C_sk <- V_shell * rho_body * c_body_outer
      C_c <- V_inner * rho_body * c_body_inner
      ATOT<- 2 * pi * R1 ^ 2 + 2 * pi * R1 * ALENTH # total surface area, m2
      ATOT_inner <- 2 * pi * (R1 - x_shell) ^ 2 + 2 * pi * (R1 - x_shell) * (ALENTH - x_shell)#(ATOT / V) * V_shell
      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)) ^ (1 / 3) # axis A, m
      B1 <- A1 * shape_b # axis B, m
      C1 <- A1 * shape_c # axis C, m
      A2 <- A1 - x_shell
      B2 <- B1 - x_shell
      C2 <- C1 - x_shell
      V_inner <- (4 / 3) * pi * (A1 - x_shell) * (B1 - x_shell) * (C1 - x_shell)
      V_shell <- V - V_inner
      C_sk <- V_shell * rho_body * c_body_outer
      C_c <- V_inner * rho_body * c_body_inner
      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
      ATOT_inner <- (4 * pi * (((A2 ^ P1 * B2 ^ P1 + A2 ^ P1 * C2 ^ P1 + B2 ^ P1 * C2 ^ 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
    }

    # end geometry section ############################################################

    if (max(Zen) >= 89) {
      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 == 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) {
      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 * (Ts - Tair) / VISDYN ^ 2) # Grashof number
    Raylei <- GR * PR # Rayleigh number

    # get Nusselt for Free Convect
    if (geom == 1) {
      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) {
      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 Transient Equations Derivation vignette
    #R_rad <- (Ts-Trad)/(emis * 0.8 * sigma * ATOT * (Ts^4-Trad^4)) # radiative resistance, eq. 4 of Transient Equations Derivation vignette
    h_rad <- 4 * emis * sigma * ((Ts + Trad) / 2 + 273.15) ^ 3 # radiation heat transfer coefficient, eq. 46 of Transient Equations Derivation vignette
    R_rad <- 1 / (h_rad * ATOT) # radiative resistance
    Q_resp <- 0 # respiratory heat loss
    #m_dot_bl <- 3.5 * (Ww_g / 100) / 1000 / 60 # blood flow g/(min-100g) to kg/s for lizard of weight Ww_g
    #R_bl <- 1 / (m_dot_bl * c_body_inner) # blood resistance
    R_sk <- (x_shell / 2) / (k_outer * ATOT) # resistance of skin
    R_b <- (V_inner ^ (1 / 3)) / (k_inner * ATOT_inner)
    if(geom == 1){
      R_b <- (R - x_shell) / (k_inner * ATOT_inner)
    }
    if(geom == 2){
      R_b <- min(A2, B2, C2) / (k_inner * ATOT_inner)
    }
    F <- 1 / (C_c * R_b) # from eq 70 and 71 of Transient Equations Derivation vignette
    H <- (Q_gen - Q_resp) / C_c # from eq 70 and 71 of Transient Equations Derivation vignette
    A <- (1 / C_sk) * (1 / R_b + 2 / R_sk - (2 / R_sk) / (R_sk / (2 * R_rad) + R_sk / (2 * R_conv) + 1)) # from eq 67 of Transient Equations Derivation vignette
    B <- 1 / (R_b * C_sk) # from eq 68 of Transient Equations Derivation vignette
    D <- (Qabs + Trad / R_rad + Tair / R_conv) / ((R_sk / (2 * R_rad) + R_sk / (2 * R_conv) + 1) * C_sk) # from eq 69 of Transient Equations Derivation vignette
    P1 <- (-1 * (F + A) + ((F + A) ^ 2 - 4 * F * A + 4 * F * B) ^ (1 / 2)) / 2 # from eq 85 of Transient Equations Derivation vignette
    P2 <- (-1 * (F + A) - ((F + A) ^ 2 - 4 * F * A + 4 * F * B) ^ (1 / 2)) / 2 # from eq 85 of Transient Equations Derivation vignette
    alpha <- (Tc + (A + D) / (B - A)) # from eq 91 of Transient Equations Derivation vignette
    beta <- F * Ts + (H * B - F * D) / (B - A) # from eq 92 of Transient Equations Derivation vignette
    C1 <- (alpha * (P2 + F) - beta) / (P2 - P1) # from eq 93 of Transient Equations Derivation vignette
    C2 <- (beta - alpha * (P1 + F)) / (P2 - P1) # from eq 94 of Transient Equations Derivation vignette
    Tcf <- -1 * (F * D + A * H) / (F * (B - A)) # from eq 77 of Transient Equations Derivation vignette
    #Tc <- C1 * exp(P1 * t) + C2 * exp(P2 * t) + Tcf # from eq 86 of Transient Equations Derivation vignette
    dTc <- (Q_gen - Q_resp - (Tc - Tsk) / R_b) / C_c # from eq 61 of Transient Equations Derivation vignette
    dTsk <- ((Tc - Tsk) / R_b - (Tsk - Ts) / (R_sk/2)) / C_sk # from eq 62 of Transient Equations Derivation vignette
    Ts <- (Qabs + Trad / R_rad + Tair / R_conv + 2 * Tsk / R_sk) / (1 / R_rad + 1 / R_conv + 2 / R_sk) # from eq 64 of Transient Equations Derivation vignette
    dTs <- Ts - y[3]
    list(y = c(dTc[1], dTsk[1], dTs[1]), x = Tcf[1])
  })
}