を使用して、流水と温度勾配のあるリング チューブをシミュレートしていdeSolve::ode()ます。リングは、各要素が温度値と位置を持つベクトルとしてモデル化されます。
私は熱拡散式をモデル化しています：
1)
しかし、リングに沿って水を動かすのにも苦労しています。理論的には、チューブ ベクトルの要素iの温度を要素s の場所の温度に置き換えるだけです。sは整数ではない可能性があるため、整数部分 ( n ) と小数部分 ( p )に分けることができます: s=n+p . したがって、水の移動による温度変化は次のようになります。
2)
問題は、ode ソルバーの各反復で評価されるdtによって、sが水の速度vに等しいことです。
私の考えは、現象を加法的に扱うことです。つまり、最初に (1) を計算し、次に (2) を計算し、最後にそれらを足し合わせます。時間の影響が怖いです。陰解法を使用する ode ソルバーは、時間ステップを自動的に決定し、ユニタリ変化デルタを線形に縮小します。
私の質問は、微分関数で (1) + (2) を返すだけが正しいのか、それとも 2 つのプロセスを分割して微分を別々に計算する必要があるのかということです。2 番目のケースでは、どのようなアプローチが推奨されますか?
編集: @tpetzoldt の提案に従って、 を使用して水の流れを実装しようとしましたReacTran::advection.1D()。私のモデルには、温度の変動の複数の原因があります。自発対称熱拡散です。水の流れ; （熱源の前に配置された）センサー付近の温度が下限しきい値を下回った場合にオンになり、上限しきい値を超えた場合にオフになる熱源。周期的な外部温度によって決定される一定の熱分散。
「Moving water」セクションの下には、以前のバージョンのコードがまだあり、現在は に置き換えられてい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()、私が水温をモデリングしているので、私のコンテキストで意味があるかどうかですが、この関数は流れる水の溶質の濃度に関連しているようです.