0

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

我正在通过以下方式对热扩散进行建模：

热扩散

它使用局部温度曲率来定义时间梯度。

正如这里所建议的deSolve：具有两个连续动力学的微分方程，我ReacTran用来解决水的运动，但我对参数化感到困惑。
作为第一次尝试，我尝试了：

ReacTran::advection.1D(v, v = vel, dx = 1, C.up = v[L], C.down = v[1], adv.method = 'Up')$dC

模拟温度，因为它是溶质的浓度。这里v表示将管道细分为L带温度的段v[i]。vel以每个时间步长定义的速度扩散dx（我对此不确定100％，所以如果有人可以对此发表评论，请+1）。我将返回的dC（浓度增量）等同于每个时间步长的温度变化dT 。

最终，在定义沿每个管段的总温度变化时，将由于热扩散引起的dT与热运动的 dT 相加。

我想知道我的方法是否正确，特别是：

  • advection.1D() 可以用来模拟温度流动，因为它是一种浓度吗？
  • 简单地将由于热扩散的dT添加到与水运动相关的 dT 是否正确，即两者是否相加？

型号详情：

我的模型有多种温度变化源：自发对称热扩散；水流；如果传感器（放置在热源之前）附近的温度低于下限阈值时打开，如果高于上限阈值则关闭热源；由周期性外部温度确定的恒定热扩散。

该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

            trigger_root <- 0
        }

        if (y[sensorP] > thresh[2] & heat == TRUE & t > 1) { # if heat is TRUE -> T was below the threshold

            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
        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 based on cycling ambient temperature
        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°')
    }
}

这是 50000 个时间点的模型图。每个子图的最后一点显示加热是打开还是关闭以及它的温度；两条水平线定义了打开/关闭加热的温度阈值。时间/空间的趋势似乎很现实，但我不确定它是否具有物理意义。

管道图

4

0 回答 0