Function to calculate temperature change around borehole heat exchanger as a function of depth and time, based on the Finite Line Source (Eskilson, 1987).

BHETemp.R

library("sf")
library("terra")
library("viridis")

# Finite Line Source (Eskilson, 1987)  
# FLS
# s = term to be integrated, ranges from D to D+L (m)
# aB = temperaturleitfähigkeit (m2/s)
# D = Überdeckung der Sonde (m)
# L = Sondenlänge (m)
# rB = Radius (m)
# z = Tiefe (m)
# not defined by r = 0 and s = z
dFLS <- function(s, aB, tt, rad, z){
  # rad[rad < 0.1] <- 0.1
  rad2 <- rad * rad
  sqrt_4abtt <- sqrt(4 * aB * tt)
  r_p <- sqrt(rad2 + (z - s)^2)
  r_m <- sqrt(rad2 + (z + s)^2)
  y <- erfc(r_p / sqrt_4abtt )/r_p - erfc(r_m / sqrt_4abtt )/r_m
  return(y)
}

# FLS <- function(tt, r, z, L, D = 0, aB){
#   integrate(dFLS, lower = D, upper = D + L, aB = aB, 
#             tt = tt, r = r, z = z)$value
# }

# for sapply in TBHE
TBHEi <- function(rad, tt, z, L, D = 0, aB, wlf, qs){
  tst <- integrate(dFLS, lower = D, upper = D + L, aB = aB, 
            tt = tt, rad = rad, z = z)
  # return(-0.5*tst$value * qs/(2*pi*wlf))
  return( tst$value )
}

# for TBHEi_rz
TBHEi_rz0 <- function(rad, z, tt, L, D = 0, aB, wlf, qs){
  # tst <- pracma::quadcc(dFLS, a = D, b = D + L, aB = aB, 
  #           tt = tt, rad = rad, z = z)
  # return(-0.5*tst * qs/(2*pi*wlf))
  tst <- integrate(dFLS, lower = D, upper = D + L, aB = aB, 
            tt = tt, rad = rad, z = z)
  #return(-0.5*tst$value * qs/(2*pi*wlf))
  return( tst$value )
}

# for outer(rad, z) in TBHE
TBHEi_rz <- function(rad, z, tt, L, D = 0, aB, wlf, qs){
 mapply(TBHEi_rz0, rad, z, MoreArgs = list(tt = tt, L = L, D = D, aB = aB, wlf=wlf, qs=qs))
}

TBHE <- function(rad, tt, z, L, D = 0, aB, wlf, qs){
  rad[rad < 0.1] <- 0.1
  if(length(rad) > 1 && length(tt) == 1 && length(z) == 1){
    TT <- sapply(rad, TBHEi, tt, z, L, D = D, aB, wlf, qs)
  }else if(length(rad) > 1 && length(tt) == 1 && length(z) > 1){
    TT <- outer(rad, z, TBHEi_rz, tt = tt, L = L, D = D, aB = aB, wlf = wlf, qs = qs)
  }
  return(-0.5*TT * qs/(2*pi*wlf) )
}

## and the so-called 'complementary error function'
erfc <- function(x) 2 * pnorm(x * sqrt(2), lower = FALSE)

# x -> sf object
localDT <- function(x, r, aB, wlf, D = 0, tt = 50, resolution = NULL, radmin = 1, 
                    method = c("linear", "nearest", "pchip", "cubic", "spline"),
                    n = NULL){
  method <- match.arg(method)
  #print(method)
  Xvect <- terra::vect(x)
  Dr <- terra::distance(r, Xvect)
  if(is.null(resolution)) resolution <- min(terra::res(r))
  #v <- seq(from = 0, to = max(Dr[]), by = resolution)
  #radius <- seq(from = 0, to = max(Dr[]), by = resolution)
  if(is.null(n)){
    n <- round(max(Dr[])/resolution)
  }
  radius <- expspace(radmin, max(Dr[]), n )
  
  DT <- TBHE(rad = radius, tt = tt, z = D + x$L/2, L = x$L, D = D, aB = aB, wlf = wlf, qs = x$q)
  
  
  # DT <- numeric(length(v))
  # for(u in seq_along(v)){
  #   # tst <- integrate(FLS, lower = D, upper = D + x$L, aB = aB, 
  #   #                  tt = tt, r = v[u], z = D + x$L/2)
  #   # DT[u] <- -0.5*tst$value * x$q/(2*pi*wlf)
  #   DT[u] <- TBHE(tt, r = v[u], z = D + x$L/2, L = x$L, D = D, aB = aB, wlf = wlf, qs = x$q)
  #   
  # }
  # plot(v, DT)
  rk <- r
  rk[] <- signal::interp1(x = radius,
                          y = DT,
                          xi = Dr[],
                          interp = method,
                          extrap = TRUE)
  return(r + rk)
}

# generate n samples exponentially distributed between u0 and un
expspace <- function(u0, un, n){
  v <- (un/u0)^(1/(n - 1))
  u <- u0 * v^(0:(n- 1))
  return(u)
}