deSolve环形水管温度模拟：扩散平流叠加逻辑与ReacTran适配疑问

问题背景 #

我正在使用deSolve::ode()模拟带温度梯度的流动水环形管道，管道被建模为向量，每个元素对应一个位置的温度值。

初期建模遇到的问题 #

我首先对热扩散公式进行建模：

1. 热扩散公式
   
   但我在叠加环内水流运动的过程中遇到了问题。理论上只需将管道向量中第i个单元的温度替换为前s个位置单元的温度即可，s可能不是整数，可拆分为整数部分（n）和小数部分（p）：s=n+p，因此水流运动带来的温度变化公式如下：
2. 水流运动温度变化公式
   
   问题在于s等于ode求解器每次迭代时的水流速度v乘以dt。我原本打算将两种过程的影响叠加，先计算（1）的结果再计算（2）的结果，最后将二者相加，但我担心时间步长的影响：采用隐式方法的ode求解器会自动决定时间步长，且会对单位变化量delta进行线性缩放。

初期疑问 #

• 在导数函数中直接返回（1）+（2）的结果是否正确？
• 还是应该将两个过程拆分后分别计算导数？如果需要拆分，建议的实现方式是什么？

后续更新 #

根据@tpetzoldt的建议，我尝试使用ReacTran::advection.1D()实现水流运动的模拟。我的模型包含多个温度变化源：

• 自发对称热扩散
• 水流运动
• 阈值控制的热源（热源前的传感器温度低于下限阈值时开启，高于上限阈值时关闭）
• 随外部温度周期变化的恒定热耗散
  下方代码段里保留了我最初的水流运动实现版本，目前已替换为ReacTran::advection.1D()。plot_type参数支持两种可视化模式：「pipe」展示水管内温度的时间序列，「sensors」展示加热器前后传感器的温度序列。

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()原本多用于模拟流动水中溶质的浓度，在我的水温模拟场景中使用该函数是否合理？ #

问题解答 #

ReacTran::advection.1D()用在水温模拟场景里完全合理。
平流过程本质上描述的是随流体运动的标量输运规律，不管这个标量是溶质浓度还是温度，控制方程的形式是完全一致的，区别只是物理量的单位和参数的物理意义，核心的离散求解逻辑没有任何差异。你模型里设置的环形边界（C.up = v[L], C.down = v[1]）也完全匹配环形管道的物理场景，实现是正确的。
另外针对你最开始的时间步长疑问补充：在导数函数里直接叠加扩散项和平流项的导数是正确做法，deSolve的求解器会自动处理时间积分的步长和精度，不需要手动拆分两个过程分步计算，拆分反而可能引入额外的算子分裂误差。 #

内容的提问来源于stack exchange，提问作者Bakaburg