前言

如题，本章主要准备讲解的是：基于蒙特卡洛(MC)和pennes模型的三维生物热传导的具体理论原理和实现。预计这章写完可以开始写基于CUDA的GPU运算加速(大概)。

附上前文链接：

生物热仿真(1)：无生命材料的热传导简析

生物热仿真(2)：蒙特卡洛一维仿真

理论

从傅里叶到pennes

在生物热仿真(1)：无生命材料的热传导简析中提到的一维傅里叶热传导经典公式如下：

$\frac{\partial T}{\partial t}=\alpha \frac{\partial^{2} T}{\partial x^{2}} \tag{1}$

式 $(1)$ 中假设了 $\alpha$ 是均匀的，将此式拓展到三维可得到：

$\frac{\partial T}{\partial t}=\alpha \nabla^2 T \tag{2}$

其中有 $\alpha = k/\rho c$ ，并且 $\alpha$ 不一定均匀，所以上式也可写为：

$\rho c \frac {\partial T}{\partial t}=\nabla (k \nabla T) \tag{2*}$

为式 $(2^*)$ 添加血液(blood)传热项 $Q_b$ 和代谢(metabolism)产热项 $Q_m$ ，得到pennes模型^[1]的传热式：

$\rho c \frac{\partial T}{\partial t}=\nabla (k \nabla T) + Q_b + Q_m \tag{3}$

pennes模型属于生物传热的最经典模型之一，在许多相关研究中都有用到，例如低温外科手术方面的Lisa X^[2]，肿瘤加热治疗方面的李和杰^[3]，在微波温热治疗方面的Thamire^[4]等人，均使用pennes模型作为传热仿真的理论基础。

pennes模型具体说明

在pennes模型中，有：
$Q_{\mathrm{b}}=W_{\mathrm{b}} C_{p, \mathrm{b}}\left(T_{\mathrm{a}}-T\right), \tag{4}$
其中 $W_{\mathrm{b}}$ 为体积血液灌注率（ $kg/m^{3} \cdot s$ ），对于此参数的选择尚无特别好的结论； $C_{p, \mathrm{b}}$ 为血液比热容； $T_{\mathrm{a}}$ 为动脉血液温度，一般可设为体核温度（37度）。

$Q_{m}=\left\{\begin{array}{ll}4,200 \mathrm{W} / \mathrm{m}^{3}, & x, y, z \notin \Omega_{t} \\ 42,000 \mathrm{W} / \mathrm{m}^{3}, & x, y, z \in \Omega_{t}\end{array}\right. \tag{5}$

$Q_m$ 的计算依赖于具体的加热方式，对于RF加热常使用比吸收率（SAR）进行计算；对于单纯的传导式加热则可使用固定的代谢值进行处理,例如Deng^[5]使用式 $(5)$ 作为肤下高度血管化肿瘤的代谢值。

三维均匀热传导率物体的显式差分迭代

将公式 $(3)$ 中的温度场设定为三维温度场，并认为导热率不变，可得到如下形式的公式：

$\frac{\partial T}{\partial t} = \alpha (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} + \frac{\partial^{2} T}{\partial z^{2}}) + \frac{Q_b}{\rho c } + \frac{Q_m}{\rho c }, \tag{6}$

根据公式 $(4)$ 和 $(5)$ ，上式可化为：
$\frac{\partial T}{\partial t}= \alpha \nabla ^2 T+ \frac{W_{\mathrm{b}} C_{p, \mathrm{b}}}{\rho c } \left(T_{\mathrm{a}}-T\right) + \frac{Q_m}{\rho c } , \tag{7}$

简单整理，将常数项移到最右边：

$\frac{\partial T}{\partial t}= \alpha (\frac {\partial ^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} + \frac{\partial^{2} T}{\partial z^{2}}) -\frac{W_{\mathrm{b}} C_{p, \mathrm{b}}}{\rho c }T + Q, \tag{7'}$

其中有 $Q = \frac{Q_m}{\rho c } + \frac{W_{\mathrm{b}} C_{p, \mathrm{b}}}{\rho c }T_{\mathrm{a}}$ ，将上式如生物热仿真(1)：无生命材料的热传导简析中进行一样的处理，时间采用前向差分，空间采用中间差分，三维空间步长一致（ $\Delta x = \Delta y = \Delta z$ ），并对其中无导数的 $T=T(x,y,z,t)$ 使用
$T = \beta T(x_{i},y_{j},z_{k}, t_{s+1}) + (1-\beta) T(x_{i},y_{j},z_{k}, t_{s}) \tag{8}$
进行松弛，则对第 $s$ 个时刻，坐标为 $(i,j,k)$ 的位置的温度迭代式为：
$\begin{array}{c} \frac{T\left(x_{i}, y_{j}, z_{k}, t_{s+1}\right)-T\left(x_{i}, y_{j}, z_{k}, t_{s}\right)}{\Delta t}= \\ \frac{\alpha}{(\Delta x)^{2}}\left[T\left(x_{i+1}, y_{j}, z_{k}, t_{s}\right)+T\left(x_{i-1}, y_{j}, z_{k}, t_{s}\right)-2 T\left(x_{i}, y_{j}, z_{k}, t_{s}\right)\right] \\ +\frac{\alpha}{(\Delta x)^{2}}\left[T\left(x_{i}, y_{j+1}, z_{k}, t_{s}\right)+T\left(x_{i}, y_{j-1}, z_{k}, t_{s}\right)-2 T\left(x_{i}, y_{j}, z_{k}, t_{s}\right)\right] \\ +\frac{\alpha}{(\Delta x)^{2}}\left[T\left(x_{i}, y_{j}, z_{k+1}, t_{s}\right)+T\left(x_{i-1}, y_{j}, z_{k-1}, t_{s}\right)-2 T\left(x_{i}, y_{j}, z_{k}, t_{s}\right)\right] \\ -\frac{W_{\mathrm{b}} C_{p, \mathrm{b}}}{\rho c}\left[\beta T\left(x_{i}, y_{j}, z_{k}, t_{s+1}\right)+(1-\beta) T\left(x_{i}, y_{j}, z_{k}, t_{s}\right)\right]+Q \end{array} \tag{9}$
为了简化上式，使用如下定义：
$X = (x_{i},y_{j},z_{k}), X_{1}^{r} = (x_{i+r},y_{j},z_{k}),\\ X_{2}^{r} = (x_{i},y_{j+r},z_{k}), X_{3}^{r} = (x_{i},y_{j},z_{k+r}), \tag{9.1}$

$F = \frac {\alpha \Delta t}{(\Delta x)^2}, W = \frac {W_{\mathrm{b}} C_{p, \mathrm{b}}}{\rho c } \tag{9.2}$

然后可得简化式(m为维度，现在m=3)：
$T(X, t_{s+1})-T(X, t_{s}) = F\sum _{i=1}^{m}{ [ T(X_{i}^{+1},t_s)+T(X_{i}^{-1},t_s)-2T(X,t_s) ]}\\ -W\Delta t[\beta T(X,t_{s+1})+(1-\beta)T(X,t_{s})] + Q\Delta t \tag{10}$
移项可得：
$\begin {array}{l} (1+W \beta \Delta t) T\left(X, t_{s+1}\right)=[1-W(1-\beta) \Delta t-2 m F] T\left(X, t_{s}\right) \\ \quad+F \sum_{i=1}^{m}\left[T\left(X_{i}^{+1}, t_{s}\right)+T\left(X_{i}^{-1}, t_{s}\right)\right]+Q \Delta t \end {array} \tag{11}$
最后可得迭代式：
$\begin{array}{c} T\left(X, t_{s+1}\right)=\frac{1-W(1-\beta) \Delta t-2 m F}{1+W \beta \Delta t} T\left(X, t_{s}\right) \\ +\frac{F}{1+W \beta \Delta t} \sum_{i=1}^{m}\left[T\left(X_{i}^{+1}, t_{s}\right)+T\left(X_{i}^{-1}, t_{s}\right)\right]+\frac{Q \Delta t}{1+W \beta \Delta t} \end{array} \tag{12}$
为了使所有的 $T(X,t)$ 项的系数之和为1，从而符合MC方式的要求，对右边提取公因式，首先设 $F_{1} = 1-W(1-\beta)\Delta t, F_{2}=1+W\beta \Delta t$ ，可将上式化为:

$\begin{array}{c} T\left(X, t_{s+1}\right)=\frac{F_{1}}{F_{2}}\left[\frac{F_{1}-2 m F}{F_{1}} T\left(X, t_{s}\right)\right. \\ \left.+\frac{F}{F_{1}} \sum_{i=1}^{m} T\left(X_{i}^{+1}, t_{s}\right)+\frac{F}{F_{1}} \sum_{i=1}^{m} T\left(X_{i}^{-1}, t_{s}\right)+\frac{Q \Delta t}{F_{1}}\right] \end{array} \tag{13}$

从式 $(13)$ 中可以看到，三维的显式差分方程和一维差别不算很大，由s时刻的该位置和周围位置的温度，可以计算s+1时刻的温度；另外，基于pennes模型，多考虑了一项 $Q$ 来代表代谢产热和血流散热的影响，更加适用于生物传热场景。

当令 $W=0$ 时，有 $F_1=F_2, Q=0$ ，上式会回归为工程传热的三维公式，可见pennes模型和传统模型是完全兼容的。

显式差分的问题

使用MC方式计算式 $(13)$ ，首先需要计算概率，各个位置的概率如下：

原位置： $1-\frac{6F}{F_1}$
邻居位置： $\frac{F}{F_1}$

松弛因子为 $\beta=0.5$ ，热扩散率 $\alpha=0.432$ ，则有：

$W = \frac{W_{\mathrm{b}} C_{p, \mathrm{b}}}{\rho c } = 0.0005 \mathrm{W} / \mathrm{m}^{3} \cdot K$
$F_1=1-0.0005\times0.5\times\Delta t=1-0.00025\Delta t$
$F=\frac{\alpha \Delta t}{(\Delta x)^2}=0.432\Delta t/(\Delta x)^2$

为了维持差分稳定性，需要满足 $0 < \frac{F}{F_1} < \frac{1}{6}$ ，由于 $\Delta t$ 常常小于1，所以只需要满足 $F_1-6F>0$ 即可。

所以稳定性条件为： $[W(1-\beta) + \frac{6\alpha}{(\Delta x)^2}] \Delta t < 1$ ，当 $\Delta x$ 较小时， $\Delta t$ 会被次方级的限制至一个非常小的空间。

例如，若我们使用 $0.1mm=0.0001m$ 的空间步长，则有
$\Delta t < \frac{(\Delta x)^2}{2.592+0.00025(\Delta x)^2} \approx 3.8580e{-9}$

隐式差分的迭代式

由于显式已经推导过，推导过程类似，迭代式为：
$\begin{array}{l} T\left(X, t_{s+1}\right)=\frac{1-W(1-\beta) \Delta t}{1+2 m F+W \beta \Delta t} T\left(X, t_{s}\right)+\frac{Q \Delta t}{1+2 m F+W \beta \Delta t} \\ \quad+\frac{F}{1+2 m F+W \beta \Delta t} \sum_{i=1}^{m}\left[T\left(X_{i}^{+1}, t_{s+1}\right)+T\left(X_{i}^{-1}, t_{s+1}\right)\right] \end{array} \tag{12*}$
然而，虽然隐式差分能克服一定的稳定性问题，但由于本身需要联立求解，对于多线程加速十分不友好。

P.S. 补充一个稳定性分析

只要所有系数都大于0即可满足MC方法下的稳定性，即

$1-W(1-\beta)\Delta t > 0$

$1+2mF+W\beta \Delta t > 0$ 在 $\beta >0, W>0$ 时,恒满足

$F = \frac{\alpha \Delta t}{(\Delta x)^2} > 0$ 恒满足

所以稳定性条件可化为 $\Delta t < \frac{1}{W(1-\beta)}$

在 $\beta = 0.5, W=0.0005$ 时，有 $\Delta t < 4000$

可见由于不依赖于 $\Delta x$ ， $\Delta t$ 的取值空间得到了很大的改善

基于蒙特卡洛的求解思路-随机游走

在调研了许多论文之后，发现基于MC思路求解的方法往往基于随机游走(Random Walk)进行实现，其核心思路在于基于统计方法获取内部点和边界点/初始点的温度关系，从而实现：边界点迭代->内部点跟随进行变化这样的过程。优势在于便于并行计算，计算量相对小等^[6]。但是，其概率的计算仍然基于有限差分，所以仍然需要满足有限差分的稳定性条件。

代码实现

随机游走算法描述

为了提高差分稳定性，此处使用隐式形式的差分。

//计算(x,y,z)位置, n次随机游走的结果
double MonteCarlo::computeTempWithRandomWalk(int x, int y, int z, int n)
{
    // 计算各个方向的概率
    std::vector<double> probs = this->computeProbs(x,y,z);
    std::vector<double> interval = GlobalFun::getProbsInterval(probs);
    std::vector<std::vector<int>> neighs;
    GlobalFun::getNeighPos({x,y,z}, neighs, this->shape->getSize().toVectorInt());

    // 开始迭代
    int N = (n < 1) ? 1 : n;
    double result = 0.0;
    std::vector<int> cur_pos = {x, y, z};
    srand((unsigned)time(NULL));

    //判断状态，流体不更新温度
    if(this->shape->getState(x,y,z)==ShapeState::FLUENT){
        result = this->temp_3d[x][y][z];
    } else {
        for (int i = 0; i < N; ++i) {
            for (int s=0; s<this->curr_iterate_time; ++s) { 
                double rand_prob = rand() / double(RAND_MAX);
                for(int neigh_idx=0; neigh_idx<neighs.size(); ++neigh_idx){
                    if(rand_prob>=interval[neigh_idx] && rand_prob<=interval[neigh_idx+1]){
                        cur_pos = neighs[neigh_idx];
                    }
                }
            }
            result += this->temp_3d[cur_pos[0]][cur_pos[1]][cur_pos[2]];        
        }
        result /= N;
        //判断邻居，如果为流体的邻居，则使用牛顿冷却定律先进行低精度仿真
        for(std::vector<int> neigh : neighs){
            if(this->shape->getState(neigh[0],neigh[1],neigh[2])==ShapeState::FLUENT){
                double delta_temp = this->temp_3d[x][y][z] - this->fluent_value;
                if(delta_temp > 0.0){
                    result -= this->k*delta_temp*this->step;
                }
            }
        }
    }
    
    return result;
}

引用

刘静，王存诚，等．生物传热学[M]．北京：科学出版社．1997 ↩
Xu LX，Zhang At．，Sandison GA，et al．A microscale model for prediction of the breast cancer cell damage during cryosurgery[C]．ASME Summer Heat Transfer Conference，July 20—23，2003 ↩
李和杰，刘静，张学学，等．激光汽化活体生物组织的传热过程分析[J]．航天医学与医学工程，2002，1(2)：103—107. ↩
Thamire CS．Bellary S．A numerical,study of microwave thermotherapy for benign prostatic hyperplasia [C]． ASME Summer Heat Transfer Co nference，July 20—23，2003 ↩
Deng Z, Liu J. MONTE CARLO METHOD TO SOLVE MULTIDIMENSIONAL BIOHEAT TRANSFER PROBLEM[J]. Numerical Heat Transfer Part B-fundamentals, 2002, 42(6): 543-567. ↩
王辉, 吴建国. 生物传热学及其医学应用[J]. 同济大学学报:医学版, 2004, 025(002):159-162. ↩

生物热仿真(3)：由一维到三维，由一般到特殊