前言

如題，本章主要準(zhǔn)備講解的是：基于蒙特卡洛(MC)和pennes模型的三維生物熱傳導(dǎo)的具體理論原理和實現(xiàn)普监。預(yù)計這章寫完可以開始寫基于CUDA的GPU運算加速(大概)贵试。

附上前文鏈接：

生物熱仿真(1)：無生命材料的熱傳導(dǎo)簡析

生物熱仿真(2)：蒙特卡洛一維仿真

理論

從傅里葉到pennes

在生物熱仿真(1)：無生命材料的熱傳導(dǎo)簡析中提到的一維傅里葉熱傳導(dǎo)經(jīng)典公式如下：

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

式 $(1)$ 中假設(shè)了 $\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)產(chǎn)熱項 $Q_m$ 豌蟋，得到pennes模型^[1]的傳熱式：

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

pennes模型屬于生物傳熱的最經(jīng)典模型之一，在許多相關(guān)研究中都有用到桑滩，例如低溫外科手術(shù)方面的Lisa X^[2]梧疲，腫瘤加熱治療方面的李和杰^[3]，在微波溫?zé)嶂委煼矫娴?em>Thamire^[4]等人运准，均使用pennes模型作為傳熱仿真的理論基礎(chǔ)幌氮。

pennes模型具體說明

在pennes模型中，有：
$Q_{\mathrm胁澳}=W_{\mathrm该互} C_{p, \mathrm}\left(T_{\mathrm{a}}-T\right), \tag{4}$
其中 $W_{\mathrm韭畸}$ 為體積血液灌注率（ $kg/m^{3} \cdot s$ ）宇智，對于此參數(shù)的選擇尚無特別好的結(jié)論； $C_{p, \mathrm陆盘}$ 為血液比熱容普筹； $T_{\mathrm{a}}$ 為動脈血液溫度，一般可設(shè)為體核溫度（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）進(jìn)行計算；對于單純的傳導(dǎo)式加熱則可使用固定的代謝值進(jìn)行處理,例如Deng^[5]使用式 $(5)$ 作為膚下高度血管化腫瘤的代謝值酸员。

三維均勻熱傳導(dǎo)率物體的顯式差分迭代

將公式 $(3)$ 中的溫度場設(shè)定為三維溫度場蜒车，并認(rèn)為導(dǎo)熱率不變，可得到如下形式的公式：

$\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}$

根據(jù)公式 $(4)$ 和 $(5)$ 幔嗦，上式可化為：
$\frac{\partial T}{\partial t}= \alpha \nabla ^2 T+ \frac{W_{\mathrm酿愧} C_{p, \mathrm}}{\rho c } \left(T_{\mathrm{a}}-T\right) + \frac{Q_m}{\rho c } , \tag{7}$

簡單整理邀泉，將常數(shù)項移到最右邊：

$\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嬉挡} C_{p, \mathrm}}{\rho c }T + Q, \tag{7'}$

其中有 $Q = \frac{Q_m}{\rho c } + \frac{W_{\mathrm汇恤} C_{p, \mathrm庞钢}}{\rho c }T_{\mathrm{a}}$ ，將上式如生物熱仿真(1)：無生命材料的熱傳導(dǎo)簡析中進(jìn)行一樣的處理因谎，時間采用前向差分基括，空間采用中間差分，三維空間步長一致（ $\Delta x = \Delta y = \Delta z$ ）财岔，并對其中無導(dǎo)數(shù)的 $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}$
進(jìn)行松弛风皿，則對第 $s$ 個時刻河爹，坐標(biāo)為 $(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} C_{p, \mathrm桐款}}{\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} C_{p, \mathrm鲁僚}}{\rho c } \tag{9.2}$

然后可得簡化式(m為維度炊苫，現(xiàn)在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)$ 項的系數(shù)之和為1，從而符合MC方式的要求冰沙，對右邊提取公因式侨艾，首先設(shè) $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$ 來代表代謝產(chǎn)熱和血流散熱的影響盖灸，更加適用于生物傳熱場景蚁鳖。

當(dāng)令 $W=0$ 時，有 $F_1=F_2, Q=0$ 赁炎，上式會回歸為工程傳熱的三維公式醉箕，可見pennes模型和傳統(tǒng)模型是完全兼容的。

顯式差分的問題

使用MC方式計算式 $(13)$ 徙垫，首先需要計算概率讥裤，各個位置的概率如下：

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

松弛因子為 $\beta=0.5$ ，熱擴散率 $\alpha=0.432$ 姻报，則有：

$W = \frac{W_{\mathrm己英} C_{p, \mathrm}}{\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$

為了維持差分穩(wěn)定性吴旋，需要滿足 $0 < \frac{F}{F_1} < \frac{1}{6}$ 损肛，由于 $\Delta t$ 常常小于1，所以只需要滿足 $F_1-6F>0$ 即可荣瑟。

所以穩(wěn)定性條件為： $[W(1-\beta) + \frac{6\alpha}{(\Delta x)^2}] \Delta t < 1$ 荧关，當(dāng) $\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}$

隱式差分的迭代式

由于顯式已經(jīng)推導(dǎo)過仙辟，推導(dǎo)過程類似同波，迭代式為：
$\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*}$
然而，雖然隱式差分能克服一定的穩(wěn)定性問題叠国，但由于本身需要聯(lián)立求解未檩，對于多線程加速十分不友好。

P.S. 補充一個穩(wěn)定性分析

只要所有系數(shù)都大于0即可滿足MC方法下的穩(wěn)定性粟焊，即

$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$ 恒滿足

所以穩(wěn)定性條件可化為 $\Delta t < \frac{1}{W(1-\beta)}$

在 $\beta = 0.5, W=0.0005$ 時冤狡，有 $\Delta t < 4000$

可見由于不依賴于 $\Delta x$ ， $\Delta t$ 的取值空間得到了很大的改善

基于蒙特卡洛的求解思路-隨機游走

在調(diào)研了許多論文之后项棠，發(fā)現(xiàn)基于MC思路求解的方法往往基于隨機游走(Random Walk)進(jìn)行實現(xiàn)悲雳，其核心思路在于基于統(tǒng)計方法獲取內(nèi)部點和邊界點/初始點的溫度關(guān)系，從而實現(xiàn)：邊界點迭代->內(nèi)部點跟隨進(jìn)行變化這樣的過程香追。優(yōu)勢在于便于并行計算合瓢，計算量相對小等^[6]。但是透典，其概率的計算仍然基于有限差分晴楔，所以仍然需要滿足有限差分的穩(wěn)定性條件。

代碼實現(xiàn)

隨機游走算法描述

為了提高差分穩(wěn)定性峭咒，此處使用隱式形式的差分税弃。

//計算(x,y,z)位置, n次隨機游走的結(jié)果
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));

    //判斷狀態(tài)，流體不更新溫度
    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;
        //判斷鄰居凑队，如果為流體的鄰居则果，則使用牛頓冷卻定律先進(jì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;
}

引用

劉靜，王存誠顽决，等．生物傳熱學(xué)[M]．北京：科學(xué)出版社．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 ?
李和杰赋访，劉靜可都，張學(xué)學(xué)，等．激光汽化活體生物組織的傳熱過程分析[J]．航天醫(yī)學(xué)與醫(yī)學(xué)工程蚓耽，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. ?
王輝, 吳建國. 生物傳熱學(xué)及其醫(yī)學(xué)應(yīng)用[J]. 同濟大學(xué)學(xué)報:醫(yī)學(xué)版, 2004, 025(002):159-162. ?

生物熱仿真(3)：由一維到三維鲤孵，由一般到特殊

生物熱仿真(3)：由一維到三維壶栋，由一般到特殊

前言

理論