我正在模拟一个带有流动水和温度梯度的环形管，使用deSolve::ode(). 环被建模为一个向量，其中每个元素都有一个温度值和位置。
我正在模拟热扩散公式：
1)
但我也在努力让水沿着环移动。从理论上讲，它只是将管向量中元素i处的温度替换为之前元素s处的温度。由于s可能不是整数，因此可以将其分为整数部分（n）和小数部分（p）：s=n+p。因此，由于水的移动引起的温度变化变为：
2)
问题是s等于在 ode 求解器的每次迭代中评估的dt的水流速度v 。
我的想法是将这些现象视为加法，即首先计算 (1)，然后计算 (2)，最后将它们相加。我害怕时间的影响。带有隐式方法的 ode 求解器自动决定时间步长，并线性按比例缩小酉变化增量。
我的问题是仅在导数函数中返回 (1) + (2) 是否正确，或者我是否应该将这两个过程分开并分别计算导数。在第二种情况下，建议的方法是什么？
编辑：根据@tpetzoldt 的建议，我尝试使用ReacTran::advection.1D(). 我的模型有多种温度变化源：自发对称热扩散；水流；如果传感器（放置在热源之前）附近的温度低于下限阈值时打开，如果高于上限阈值则关闭热源；由周期性外部温度确定的恒定热扩散。
在“动水”部分下面仍然是我以前版本的代码，现在替换为ReacTran::advection.1D(). 该plot_type参数允许可视化水管（“管道”）中温度的时间序列，或传感器（加热器之前和之后）的温度序列。
library(deSolve)
library(dplyr)
library(ggplot2)
library(tidyr)
library(ReacTran)
test <- function(simTime = 5000, vel = 1, L = 500, thresh = c(16, 25), heatT = 25,
heatDisp = .0025, baseTemp = 15, alpha = .025,
adv_method = 'up', plot_type = c('pipe', 'sensors')) {
plot_type <- match.arg(plot_type)
thresh <- c(16, 25)
sensorP <- round(L/2)
vec <- c(rep(baseTemp, L), 0)
eventfun <- function(t, y, pars) {
heat <- y[L + 1] > 0
if (y[sensorP] < thresh[1] & heat == FALSE) { # if heat is FALSE -> T was above the threshold
#browser()
y[L + 1] <- heatT
}
if (y[sensorP] > thresh[2] & heat == TRUE) { # if heat is TRUE -> T was below the threshold
#browser()
y[L + 1] <- 0
}
return(y)
}
rootfun <- function (t, y, pars) {
heat <- y[L + 1] > 0
trigger_root <- 1
if (y[sensorP] < thresh[1] & heat == FALSE & t > 1) { # if heat is FALSE -> T was above the threshold
#browser()
trigger_root <- 0
}
if (y[sensorP] > thresh[2] & heat == TRUE & t > 1) { # if heat is TRUE -> T was below the threshold
#browser()
trigger_root <- 0
}
return(trigger_root)
}
roll <- function(x, n) {
x[((1:length(x)) - (n + 1)) %% length(x) + 1]
}
fun <- function(t, y, pars) {
v <- y[1:L]
# Heat diffusion: dT/dt = alpha * d2T/d2X
d2Td2X <- c(v[2:L], v[1]) + c(v[L], v[1:(L - 1)]) - 2 * v
dT_diff <- pars * d2Td2X
# Moving water
# nS <- floor(vel)
# pS <- vel - nS
#
# v_shifted <- roll(v, nS)
# nS1 <- nS + 1
# v_shifted1 <- roll(v, nS + 1)
#
# dT_flow <- v_shifted + pS * (v_shifted1 - v_shifted) - v
dT_flow <- advection.1D(v, v = vel, dx = 1, C.up = v[L], C.down = v[1],
adv.method = adv_method)$dC
dT <- dT_flow + dT_diff
# heating of the ring after the sensor
dT[sensorP + 1] <- dT[sensorP + 1] + y[L + 1]
# heat dispersion
dT <- dT - heatDisp * (v - baseTemp + 2.5 * sin(t/(60*24) * pi * 2))
return(list(c(dT, 0)))
}
out <- ode.1D(y = vec, times = 1:simTime, func = fun, parms = alpha, nspec = 1,
events = list(func = eventfun, root = T),
rootfunc = rootfun)
if (plot_type == 'sensors') {
## Trend of the temperature at the sensors levels
out %>%
{.[,c(1, sensorP + 1, sensorP + 3, L + 2)]} %>%
as.data.frame() %>%
setNames(c('time', 'pre', 'post', 'heat')) %>%
mutate(Amb = baseTemp + 2.5 * sin(time/(60*24) * pi * 2)) %>%
pivot_longer(-time, values_to = "val", names_to = "trend") %>%
ggplot(aes(time, val)) +
geom_hline(yintercept = thresh) +
geom_line(aes(color = trend)) +
theme_minimal() +
theme(panel.spacing=unit(0, "lines")) +
labs(x = 'time', y = 'T°', color = 'sensor')
} else {
## Trend of the temperature in the whole pipe
out %>%
as.data.frame() %>%
pivot_longer(-time, values_to = "val", names_to = "x") %>%
filter(time %in% round(seq.int(1, simTime, length.out = 40))) %>%
ggplot(aes(as.numeric(x), val)) +
geom_hline(yintercept = thresh) +
geom_line(alpha = .5, show.legend = FALSE) +
geom_point(aes(color = val)) +
scale_color_gradient(low = "#56B1F7", high = "red") +
facet_wrap(~ time) +
theme_minimal() +
theme(panel.spacing=unit(0, "lines")) +
labs(x = 'x', y = 'T°', color = 'T°')
}
}
有趣的是，设置更高的分段数 ( L = 500) 和高速 ( vel = 2) 可以观察到后加热传感器中的尖峰序列。此外，处理时间急剧增加，但更多的是速度增加的影响，而不是管道分辨率的增加。
我现在最大的疑问是ReacTran::advection.1D()在我的上下文中是否有意义，因为我正在模拟水温，而这个函数似乎与流动水中的溶质浓度更相关。