ABSTRACT

In recent years, magnetic hyperthermia as a new tumor treatment method has been studied. In this paper, the Non-Fourier heat transfer equation with magnetic nanoparticles as the heat source is analyzed. The heat transfer model of single-phase-lag equation is established by introducing relaxation time based on the Fourier heat transfer model, and the analytical solution of the equation is obtained by using the method of separating variables. The effects of external magnetic field and properties of magnetic nanoparticles on temperature distribution are investigated.

Keywords:Non-Fourier Heat Transfer Model, Relaxation Time, Magnetic Hyperthermia, Magnetic Nanoparticles

磁纳米粒子为热源的非傅里叶热传输

李秋月，菅永军^*

内蒙古大学数学科学学院，内蒙古 呼和浩特

收稿日期：2020年6月8日；录用日期：2020年6月23日；发布日期：2020年6月30日

摘 要

近年来，磁感应热疗作为一种新兴的肿瘤治疗方法得到较多的关注。本文分析了以磁纳米粒子作为热源的傅里叶传热方程，在傅里叶热传导模型的基础上通过引入弛豫时间建立了单相方程的热传递模型，并利用分离变量法求得了该方程的解析解。分析了外部磁场性质和磁性纳米粒子特性对温度分布的影响。

关键词 :非傅里叶传热模型，弛豫时间，磁感应热疗，磁纳米粒子

Copyright © 2020 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

现代热疗作为一种肿瘤治疗技术，其研究热点主要集中在靶区温度的合理性和均匀性等方面上。其中磁感应热疗作为一种有效的热疗方法，具有加热均匀性好、对正常细胞组织损伤较小和联合治疗效果佳等优点。磁热疗技术是通过注射方式，将磁纳米粒子注射进肿瘤区域附近使其聚集，再施加一外部交变磁场，使磁纳米粒子因磁滞、弛豫或感应涡流等现象被加热，最终破坏肿瘤细胞 [1] [2] [3] [4]。2014年，Bohara [5] 等人对乳腺癌细胞进行了体外热疗试验，发现癌细胞死亡率可达到80%，证实了磁流体热疗治疗乳腺癌的潜在可能性。

进行磁热疗过程中最重要的工作是考虑加热时人体中温度的分布情况，精准的对病患位置进行加热且不损伤正常组织是第一需求。因此在以磁纳米粒子为热源的磁热疗中考虑温度分布是必要的。对于常规条件下稳态热传导方程，人们通常用Fourier定律来描述热流密度与温度梯度之间的关系。表述为介质内任意一点的热流密度矢量与该处的温度降度成正比 [6]。自傅里叶热传导定律提出后，一直被作为各种传热现象的基础理论。但在其描述瞬态导热过程时，与能量方程的结合所得到的抛物型方程会出现热扰动传播速度无限大的情况，即达到一个平衡状态的弛豫时间是无限趋近于零，显然是不符合物理实际的。因此，在傅里叶方程的基础上，考虑热流矢量的传播和温度梯度形成之间的时间延迟现象是必要的 [7] [8] [9] [10]。Jou [11] 提出了基于熵产理论的热波模型，解决了热力学转变过程中非平衡态超速传热的问题。Tzou [12] 从另一个角度出发，考虑热流矢量的传播和温度梯度的形成之间有一个时间延迟，提出了一个描述非平衡热力学转变的宏观方程，即单相延迟的双曲热传导方程。Bai [13] 等人的研究表明典型的双曲型热传导方程在某些条件下会低于绝对零度，此时连续介质假说不再成立。Coleman [14] 修正了其内能的表示方程，得到了修正双曲型热传导方程。

在傅里叶传热方程中加入了热流量对时间的导数，使方程变成双曲型的波动方程，表示如下：

$q+{\tau }_q\frac{\partial q}{\partial t}=-k\frac{\partial T}{\partial x}$ (1)

其中q为热流密度，是单位时间通过单位面积的热流量，T为温度，k为导热率，又称导热系数，衡量物体导热能力的参数，τ_q为弛豫时间，表示热流量与温度梯度之间的相位差。傅里叶定律作为导热的本构方程，在一个多世纪以来用来解决实际工程问题，所得结果能较好地与实际相吻合，且适用于各种情况下的导热问题。但实际情况下，由于反应时间不够，理论上在时间t形成的温度梯度其实是t + τ_q时热流矢量到达x点处才建立的。

本文拟研究以磁纳米粒子为热源的修正双曲型热传导方程。利用分离变量法，求得考虑弛豫时间的非傅里叶传热方程的解析解，得到磁场参数和磁纳米粒子性质对磁热疗中组织温度分布的影响，为实验和实际应用提供理论基础。

2. 模型描述

将式(1)热流矢量方程关于t进行泰勒级数展开后，再代入Pennes方程 [15] (因为相位延迟时间τ_q很小，故忽略其二次项和高阶项，并假设热导率k为常数)，得到一个含自变量x和t的双曲型热传导方程：

$\rho c\left(1+{\tau }_q\frac{\partial }{\partial t}\right)\frac{\partial T}{\partial t}+w_b{c}_b\left(1+{\tau }_q\frac{\partial }{\partial t}\right)\left(T-T_b\right)=k\frac{{\partial }^{2}T}{\partial x^2}+Q.$ (2)

上述模型满足以下边界条件和初始条件：

${T|}_{x=0}=T_0,{\frac{\partial T}{\partial x}|}_{x=l}=0,{T|}_{t=0}=0,{\frac{\partial T}{\partial t}|}_{t=0}=0.$ (3)

式(2)就是包含弛豫时间的双曲型传热方程，方程左边w_b和ρ分别是血液灌注率和组织和血液密度，τ_q是弛豫时间，c_b和c分别表示血液和组织的比热，Q表示磁纳米粒子产热项，T₀是x = 0处的温度，取值37℃，l是组织的厚度为0.02 m。

对上述方程进行无量纲化，其无量纲量为：

$\begin{array}{l}{x}^{*}=\frac{x}{\left[x\right]},L=\frac{l}{\left[x\right]},\left[x\right]=\sqrt{\frac{k}{w_b{c}_b}},t^{*}=\frac{t}{\left[t\right]},\left[t\right]=\frac{\rho c}{w_b{c}_b},{\tau }_q^{*}=\frac{{\tau }_q}{\left[{\tau }_q\right]},\left[{\tau }_q\right]=\frac{\rho c}{w_b{c}_b},\\ T^{*}=\frac{\left(T-T_b\right)}{\left[T\right]},\left[T\right]=T_0-T_b,Q^{*}=\frac{Q}{\left[Q\right]},\left[Q\right]=w_b{c}_b\left(T_0-T_b\right).\end{array}$ (4)

将以上无量纲量代入式(2)，(3)中，化简得到：

$\left(1+{\tau }_q^{*}\frac{\partial }{\partial t^{*}}\right)\frac{\partial T^{*}}{\partial t^{*}}+\left(1+{\tau }_q^{*}\frac{\partial }{\partial t^{*}}\right)T^{*}=\frac{\partial T^{*}}{\partial x^{*2}}+Q^{*},$ (5a)

${T^{*}|}_{x=0}=1\text{,}{\frac{\partial T^{*}}{\partial x^{*}}|}_{x^{*}=L}=0,{T^{*}|}_{t^{*}=0}=0,{\frac{\partial T^{*}}{\partial t^{*}}|}_{t^{*}=0}=0.$ (5b)

3. 方程求解

由于热源稳定发热，与温度场的位置和时间都无关。我们令： $T^{**}=T^{*}-Q^{*}$，方程(4)，(5)可化为：

$\left(1+{\tau }_q^{*}\frac{\partial }{\partial t^{*}}\right)\frac{\partial T^{**}}{\partial t^{*}}+\left(1+{\tau }_q^{*}\frac{\partial }{\partial t^{*}}\right)T^{**}=\frac{{\partial }^2{T}^{\text{**}}}{\partial x^{*}{}^2},$ (6)

${T^{\text{**}}|}_{x^{*}=0}=1-Q^{*}\text{,}{\frac{\partial T^{**}}{\partial x^{*}}|}_{x^{*}=L}=0,{T^{**}|}_{t^{*}=0}=-Q^{*},{\frac{\partial T^{\text{**}}}{\partial t^{\text{*}}}|}_{t^{*}=0}=0.$ (7)

现得到一个温度的微分方程和非齐次的边界条件。为了求解上述方程，可由叠加原理将其分解成一个齐次边界条件的偏微分方程以及一个常微分方程。为了表示方便，去掉后面所有的表达式中表示无量纲的*号。设：

$T\left(x,t\right)={\theta }_1\left(x,t\right)+{\theta }_2\left(x\right).$ (8)

其中θ₂满足方程：

$\frac{{\text{d}}^2{\theta }_2\left(x\right)}{\text{d}{x}^2}-{\theta }_2\left(x\right)=0,$ (9)

${{\theta }_2\left(x\right)|}_{x=0}=1-Q\text{,}{\frac{\partial {\theta }_2\left(x\right)}{\partial x}|}_{x=L}=0.$ (10)

解得：

(11)

θ₁(x, t)满足：

$\left(1+{\tau }_q\frac{\partial }{\partial t}\right)\frac{\partial {\theta }_1}{\partial t}+\left(1+{\tau }_q\frac{\partial }{\partial t}\right){\theta }_1=\frac{{\partial }^2{\theta }_1}{\partial x^2},$ (12)

${{\theta }_1|}_{x=0}=0,{\frac{\partial {\theta }_1}{\partial t}|}_{x=L}=0,{{\theta }_1|}_{t=0}=-Q-{\theta }_2,{\frac{\partial {\theta }_1}{\partial t}|}_{t=0}=0.$ (13)

使用分离变量法，设： ${\theta }_{\text{1}}=X\left(x\right)T\left(t\right)$，则有： $\left(1+{\tau }_q\right)X{T}^{\prime }+{\tau }_{q}X{T}^{″}+XT=X^{″}T$，

可以通过整理得到：

$\left(1+{\tau }_q\right)\frac{T^{\prime }}{T}+{\tau }_q\frac{T^{″}}{T}+1=\frac{X^{″}}{X}=-\lambda .$ (14)

设 $\frac{X^{″}\left(x\right)}{X\left(x\right)}=-\lambda$ 和 $\left(1+{\tau }_q\right)\frac{T^{\prime }}{T}+{\tau }_q\frac{T^{″}}{T}+1=-\lambda$。

先求解二阶线性常微分方程组：

$X^{″}\left(x\right)+\lambda X\left(x\right)=0,$ (15a)

$X\left(0\right)=0，X^{\prime }\left(L\right)=0.$ (15b)

利用特征值法求解X(x)，解得：

$X_n\left(x\right)=\sin \sqrt{{\lambda }_n}x.$ (16)

其中 ${\lambda }_n={\left[\frac{\left(n+1/2\right)\pi }{L}\right]}^2,n=0,1,2,\cdots$。

再对时间函数T(t)进行求解

$T^{″}+\frac{1+{\tau }_q}{{\tau }_q}{T}^{\prime }+\frac{1+{\lambda }_n}{{\tau }_q}T=0.$ (17)

若 ${\left(1-{\tau }_q\right)}^2-4{\tau }_q{\lambda }_n<0$ 则 $w_n^2=\frac{{\lambda }_n}{{\tau }_q}-{\left[\frac{\left(1-{\tau }_q\right)}{2{\tau }_q}\right]}^2$，

$T_n\left(t\right)={\text{e}}^{-\frac{\alpha t}{2}}\left(c_1\cos w_{n}t+c_2\sin w_{n}t\right),$ (18a)

若 ${\left(1-{\tau }_q\right)}^2-4{\tau }_q{\lambda }_n\ge 0$ 则 $r_n^2=-\frac{{\lambda }_n}{{\tau }_q}+{\left[\frac{\left(1-{\tau }_q\right)}{2{\tau }_q}\right]}^2$，

$T_n\left(t\right)={\text{e}}^{-\frac{\alpha t}{2}}\left(c_3{\text{e}}^{r_{n}t}+c_4{\text{e}}^{-r_{n}t}\right).$ (18b)

其中 $\alpha =\frac{\left(1+{\tau }_q\right)}{{\tau }_q}$。

解得：

${\left(1-{\tau }_q\right)}^2-4{\tau }_q{\lambda }_n<0$ 时 ${\theta }_1\left(x,t\right)={\text{e}}^{-\frac{\alpha t}{2}}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left(c_1\cos w_{n}t+c_2\sin w_{n}t\right)}\left(\sin \sqrt{{\lambda }_n}x\right),$ (19a)

${\left(1-{\tau }_q\right)}^2-4{\tau }_q{\lambda }_n\ge 0$ 时 ${\theta }_1\left(x,t\right)={\text{e}}^{-\frac{\alpha t}{2}}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left(c_3{\text{e}}^{r_{n}t}+c_4{\text{e}}^{-r_{n}t}\right)}\left(\sin \sqrt{{\lambda }_n}x\right).$ (19b)

将式(19)代入边界条件(13)求解系数可得：

$c_1=-\frac{4Q}{\left(2n-1\right)\pi }-\frac{\left(2n-1\right)\left(1-Q\right)\pi }{L^2+{\left(\frac{2n-1}{2}\pi \right)}^2},c_2=\frac{\left(1+{\tau }_q\right)}{2w_n{\tau }_q}{c}_1,c_3=\frac{\alpha +2r_n}{4r_n}{c}_1,c_4=\frac{-\alpha +2r_n}{4r_n}{c}_1.$ (20)

将式(11)和式(19)带入式(8)和得到：

${\left(1-{\tau }_q\right)}^2-4{\tau }_q{\lambda }_n<0$ 时，

$T=Q+{\text{e}}^{-\frac{\alpha t}{2}}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left(c_1\cos w_{n}t+c_2\sin w_{n}t\right)}\left(\sin \sqrt{{\lambda }_n}x\right)+\left(1-Q\right)\frac{\cosh \left(x-L\right)}{\cosh L},$ (21a)

${\left(1-{\tau }_q\right)}^2-4{\tau }_q{\lambda }_n\ge 0$ 时，

$T=Q+{\text{e}}^{-\frac{\alpha t}{2}}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left(c_3{\text{e}}^{r_{n}t}+c_4{\text{e}}^{-r_{n}t}\right)}\left(\sin \sqrt{{\lambda }_n}x\right)+\left(1-Q\right)\frac{\cosh \left(x-L\right)}{\cosh L}.$ (21b)

4. 分析和讨论

上述研究分析了存在时间滞后特性的傅里叶热传导模型，利用分离变量法求得该模型在Dirchlet和Neumann混合边界条件下的解析解，得出了温度分布的函数。在实际问题中，典型的参数值选取为：组织导热系数k = 0.5 W/(m.℃)，血液灌注率w_b = 0.52 kg/(m³.s)，动脉血温T_b = 36℃，血液比热c_b = 4186 J/(kg.℃)，组织比热c = 4186 J/(kg.℃)，组织和血液的密度为ρ = 1000 kg/m [16]。

磁疗中的热疗效应实质上是由在中频范围内磁滞损耗和弛豫τ引起。Rosensweig在文献 [17] 中给出了关于弛豫生热、磁性粒子理化性质与外场的关系的公式：

$Q=\pi {\mu }_0{\chi }_0{H}^{2}f\frac{2\pi f\tau }{1+{\left(2\pi f\tau \right)}^2},$ (22)

其中μ₀ = 4π × 10⁻⁷ N/A²是真空的磁导率，χ₀为磁化率，χ₀ = M/H，M是磁化强度描述磁性体磁性强弱的物理量，取425 kA∙m^-1，H为外加磁场强度，f为交变磁场频率。τ为磁纳米粒子的有效弛豫时间。由式(22)，影响磁纳米粒子产热的参数主要分为两个方面：外磁场的性质和磁纳米粒子自身的性质，外磁场性质中包括磁场强度和交变磁场频率。另一方面，金属粒子的弛豫时间τ会随颗粒大小改变，本文为简单起见τ取定值3.378 × 10⁻⁶ [17]。在外加磁场的选取方面，由于通常人体可耐受的磁场强度最大为10 kA/m，所以对H的选择一般小于10 kA/m。在0~10 kA/m之间选取3个不同的磁场强度，比较其温度的变化情况。以同样方式选取了3个不同交变磁场频率的数值，分别为55 kHz，65 kHz，75 kHz。

4.1. 外磁场性质改变对温度变化的影响

图1给出了不同磁场强度情况下，无量纲温度曲线随时间和距离的变化。类似的，图2是不同交变磁场频率情况下，无量纲温度曲线随时间和距离的变化：从图1(a)和图2(a)可以看出，在固定点处，温度曲线总体趋势是随着时间变化温度逐渐升高的。在初始阶段有一段微弱的波动。这段波动体现了考虑弛豫时间的传热方程的特征，即温度变化和热量传播之间存在一定的延迟。在波动阶段后温度快速升高，并在到达一定温度后逐渐趋于平缓。图1(b)和图2(b)可以看出，在固定时刻，离磁纳米粒子热源的位置距离越远，温度逐渐下降。另外，提高磁场强度或交变磁场频率，都会使得温度升高。

图1. 不同磁场强度下的无量纲温度曲线。(a) 在0.02 m处随时间变化；(b) 400 s时随距离变化

图2. 不同交变磁场频率下的无量纲温度曲线。(a) 0.02 m处随时间变化；(b) 400 s时随距离变化

4.2. 不同传热弛豫时间温度曲线的变化

图3描述了不同弛豫时间下温度随时间的变化，选取了两个位置(r= 0.02 m和r= 0.002 m)，并截取局部进行放大。可以看出，在初始阶段，弛豫时间越小的温度场温度越高，震荡的时间越短，较早的趋于稳定。这是因为弛豫时间描述了热流矢量的传播和温度梯度的形成之间存在的时间差，而弛豫时间越短意味着人体内的这个时间差越小。但这个时间差并不影响温度的大小，所以随着加热时间的增加，最终三种弛豫时间下达到的温度会趋于一致，从瞬态过渡到了稳态。越靠近热源，达到稳定的时间越短。

图3. 不同弛豫时间在(a) r = 0.02 m和(b) r = 0.002 m处随时间变化的无量纲温度曲线(H = 7 kAm⁻¹, f = 65 kHz)

图4展示了不同弛豫时间下在某一时刻随热源距离变化的温度曲线。在远离热源处，弛豫时间较长的温度曲线震荡幅度较大。在考虑弛豫时间的情况下，某点处的温度梯度的形成过程是震荡的，弛豫时间越长，震荡幅度越大，震荡时间越长，达到稳态温度越慢。

图4. 不同弛豫时间在(a) t = 100 s和(b) t= 400 s处随距离变化的无量纲温度曲线(H = 7 kAm⁻¹, f = 65 kHz)

5. 结论

为了研究磁纳米粒子的生热效应，本文建立了以磁纳米粒子为热源的考虑滞后时间修正的单项延迟Pennes热传导方程，通过分离变量法求出上述方程的解析解。分析热源项，从外磁场性质和磁性纳米粒子特性两个方面改变参数得到不同的升温曲线。证实了提高磁场强度，交变磁场频率或磁纳米粒子粒径都可以有效地提高温度场可达到的最高温度。升温曲线由于存在弛豫时间，所以在初始阶段有一定程度的震荡，但外磁场性质的改变对震荡幅度影响不大。改变传热过程中的弛豫时间并不影响温度场最终达到的温度，但是会影响初期的波动幅度。

基金项目

本研究由国家级大学生创新创业训练计划“磁纳米粒子热源的生物热传输分析(201910126021)”项目资助。

文章引用

李秋月,菅永军. 磁纳米粒子为热源的非傅里叶热传输
Non-Fourier Heat Transfer with Heat Source of Magnetic Nanoparticles[J]. 应用数学进展, 2020, 09(06): 972-979. https://doi.org/10.12677/AAM.2020.96115

参考文献