Effect of Non-Fourier Heat-Flux Boundary Conditions on Heat Conduction Behavior of Laser-Irradiated Biological Tissues

Xu Guangying; Xue Dawen; Wang Jinbao

doi:10.3788/CJL202047.1207003

Abstract

The dual phase lagging (DPL) non-Fourier heat transfer model can reflect the transient interaction process between pulsed laser and biological tissues. However, in many literatures where the DPL model is used to study the heat conduction mechanism of biological tissues, there exist both Fourier and non-Fourier boundary conditions and thus many induced conclusions are contradictory. In this paper, the control equation based on the DPL non-Fourier model is adopted and the Fourier and non-Fourier boundary conditions are derived. Meanwhile, the analytical solutions under the above conditions are obtained by integral transformation and Laplace transformation. The biological tissues are taken as an example and the calculation results show that as for the non-Fourier control equation, the predicted temperature distribution in tissues based on the non-Fourier boundary condition is in accordance with the energy conservation law, while the result based on the Fourier boundary condition is not. The conclusions on the temperature rising amplitude and temperature rising rate are opposite for the two kinds of boundary conditions. Moreover, the thermal damage predicted under the Fourier boundary conditions is overly conservative and obviously lower than that under the non-Fourier boundary conditions. Finally, from the point view of energy conservation, the DPL non-Fourier boundary condition is just the DPL energy conservation equation of boundary, while the Fourier boundary condition is the energy conservation equation of the Fourier model. The Fourier control equation of bio-heat conduction should be matched with the Fourier boundary conditions, while the non-Fourier control equation of bio-heat conduction should be matched with the non-Fourier boundary conditions.

Keywords

biotechnology DPL non-Fourier boundary conditions Laplace transformation no-Fourier bio-heat conduction thermal wave

OCIS code 170.6935; 080.2720; 140.6810

1 引言

激光技术已广泛应用于现代医学治疗^[1-2]。在激光治疗中,将适量的热量输送到指定部位且不对周围其他组织造成热损伤,这依赖于对生物组织热传导、温度响应及温度分布的正确预测。

显然,建立能够准确描述激光辐照下生物组织传热的数学模型对激光热疗具有重要的指导意义。迄今为止,已有多种生物传热模型被提出,包括基于傅里叶定律的生物传热模型^[3]和基于单相滞(single phase lag, SPL)和双相滞(dual phase lag, DPL)定律的非傅里叶模型^[4-7]。在傅里叶定律中,热流密度与温度梯度成正比。

傅里叶模型隐含着热扰动的无限传播速度^[8-9],难以描述激光脉冲加热引起的组织传热非平衡状态。单相滞定律也称作C-V定律,其描述了热流及热流相延迟项与温度梯度的正比关系。热流相延迟项表征热流的速度。基于单相滞定律建立的热传导方程能够描述热扰动的快速瞬态过程,但不能反映生物体组织内部血液和组织之间的传热作用。双相滞定律表述了热流及热流相延迟与温度梯度及温度梯度相延迟之间的正比关系。温度梯度相延迟项表征内部组织和血液之间的热传递。双相滞生物传热模型既能描述热扰动的快速瞬态过程,又能表征生物微观结构相互作用效应,被广泛应用于研究生物热现象,特别是在烧伤、冷冻、激光和肿瘤热疗等领域。相关研究^[10-20]均指出,DPL模型和傅里叶传热模型对温度的预测存在很大的差异。DPL模型基于热流相延迟时间 $τ_{q}$ 和温度梯度相延迟时间 $τ_{T},$ 能够描述傅里叶模型无法描述的波动效应,同时也能反映组织和血液之间的相互热传递。目前研究者对DPL生物传热模型的研究主要集中在控制方程上,但边界条件的运用很混乱。一些研究使用非傅里叶模型时,出现热传导控制方程和边界条件不一致的问题;傅里叶模型、SPL模型和DPL模型存在同时使用非傅里叶边界条件和傅里叶边界条件的问题。杨洪钦等^[10]基于DPL传热模型,研究了皮肤组织受高斯激光辐照的热响应和热损伤问题,其中激光源作为内热源处理,考虑了边界的对流冷却,给出的是傅里叶第三类边界条件。Kumar等^[11]利用傅里叶模型、热波(C-V)模型和DPL模型,研究了皮肤组织受不同表面热源激励时组织内部的温度响应和热损伤,但给出的边界条件都是傅里叶边界条件。Ziaei等^[12]和Ahmadikia等^[13]分别基于C-V模型和DPL模型,研究了皮肤受常热流和周期热流作用时的温度响应和热损伤变化,但给出的均是傅里叶Neumman边界条件。当任意处的温度发生变化时,经典傅里叶边界条件预测的温度变化的速度最快,升温或降温幅度最大。C-V模型预测的温度是有限速度热波,其变化速度最慢,热损伤最小。DPL传热模型预测的结果介于二者之间。从传热模型分析,热波模型描述的是有限速度的热传播,这意味着微观传热载流子交换能量是需要时间的,传播速度较慢,这意味着能量在当地易积累;经典傅里叶模型描述的是具有无限传播速度的热波,这意味着传热载流子是瞬时完成能量交换的,传播速度较快,能量在当地不易累积。有累积效应的当地温度应该高于没有累积效应的当地温度,因此非傅里叶模型预测的温度应该高于傅里叶模型预测的温度。温度越高热损伤越大,因此非傅里叶模型预测的热损伤也应该高于傅里叶模型预测的热损伤。1973年,Maurer等^[14]处理高通量热流边界条件时,通过引入热流的松弛效应,成功解决了瞬时表面温度突然升高的问题,预测结果比傅里叶边界条件预测的结果更符合现实。2002年,Al-Nimr等^[15]首次采用微元表面单元热平衡法,建立了C-V模型的非傅里叶边界条件。Antaki^[6]利用非傅里叶模型成功解释了冷藏处理肉中的温度波变化。Antaki^[16]基于热波模型,研究了第三类边界条件下的热波传播特征,并给出了非傅里叶第三类边界条件,研究结果表明,表面温度变化存在突然跳跃的特点,而傅里叶模型预测的结果没有这个特点。Shirmohammadi等^[17]直接将激光热流代入C-V方程,并将其作为非傅里叶边界条件。Biswas等^[18]对傅里叶定律方程进行了延迟化处理,并将其作为边界条件,然后将其泰勒级数一阶展开作为边界条件,这里给出的非傅里叶边界条件在形式上与文献[ 15]是一致的。

事实上,对于生物组织热疗领域的非傅里叶热传导方程,经典的Neumann边界条件并不适用于非傅里叶传热方程。研究者对皮肤组织温度响应的非傅里叶边界条件进行了相关研究。Furghani等^[19]对比了傅里叶Neumman边界条件和非傅里叶Neumman边界条件对传热的影响,但只比较了同样条件下的温度变化特点,没有对不同边界条件造成的不同传热机制给予物理上的解释,也没有明确说明如何应用非傅里叶边界条件。Kumar等^[20]虽然也给出了通用非傅里叶DPL边界条件,但给出的通用非傅里叶边界条件与文献[ 13]中给出的第三类边界条件是一致的,且与文献[ 15]中给出的第二类边界条件相比,热流没有松弛时间项,这种形式是值得商榷的。王晋宝等^[21]研究了受周期热流作用的皮肤组织温度和热损伤变化的特点,通过直接采用C-V非傅里叶边界条件,得出了组织内部温度和热损伤变化幅度均大于傅里叶模型下的结果,这种变化特点符合非傅里叶模型的基本特征。

对于已经刊出的关于傅里叶边界条件和非傅里边界条件的研究,很少有实验加以佐证,文献[ 10]、[12]与[21]同样是采用生物组织的非傅里叶模型,为什么在热流作用下预测的温度和热损伤规律不一样?在文献[ 19]中,同一非傅里叶传热方程分别施加傅里叶边界条件和非傅里叶边界条件,得出的传热机理不同。目前,这方面存在的问题有:1)在物理上,非傅里叶边界条件和傅里叶边界条件的物理意义区别在哪里?2)在数学上,傅里叶边界条件和非傅里叶边界条件如何配合对应的热传导方程?针对这些问题,本文从能量守恒方程出发,从理论上导出了非傅里叶Neumann边界条件。利用DPL生物传热方程,得到了基于经典傅里叶Neumann边界条件和非傅里叶Neumann边界条件的解析解,通过比较两种不同边界条件下的温度分布、温度响应和热损伤变化,详细讨论了非傅里叶Neumann边界条件的有效性和合理性。

2 数学模型

2.1 控制方程

本文重点讨论双相滞生物传热模型。考虑图1所示的生物组织,其厚度为 $L 。$ 假设组织均匀,各向同性,初始温度为灌注血流温度 $T_{a}, 在 t > 0^{+}$ 时刻(0⁺为大于0的无限小时间),左侧突然受到激光 $J_{0} f (t)$ 照射,其中, $J_{0}$ 为激光强度, $f (t)$ 为时变函数,且组织表面对光源表现出强烈的吸收性,右侧按绝热处理。

Figure 1.Schematic of laser heating of biological tissues

生物组织内部的传热服从双相滞非傅里叶本构关系:

- k [\frac{\partial T}{\partial x} + τ_{T} \frac{\partial^{2} T}{\partial x \partial t}] = q (x, t) + τ_{q} \frac{\partial q}{\partial t}, (1)

式中: $k$ 为组织的导热系数; $T$ 为组织的当地温度; $q$ 为组织的当地热流密度; $τ_{T}$ 为温度梯度相延迟时间,表征组织内部的相互作用; $τ_{q}$ 为热流相延迟时间,表征组织热载体之间的相互作用。

结合引入血流灌注的Pennes生物能量守恒方程:

ρ_{t} c_{t} \frac{\partial T}{\partial t} = - \frac{\partial q}{\partial x} + ρ_{b} c_{b} {\tilde{w}}_{b} (T_{a} - T) + q_{met}, (2)

得到双相滞后生物传热方程为

ρ_{t} c_{t} τ_{q} \frac{\partial^{2} T}{\partial t^{2}} + (τ_{q} {\tilde{w}}_{b} ρ_{b} c_{b} + ρ_{t} c_{t}) \frac{\partial T}{\partial t} = k (\frac{\partial^{2} T}{\partial x^{2}} + τ_{T} \frac{\partial}{\partial t} \frac{\partial^{2} T}{\partial x^{2}}) - {\tilde{w}}_{b} ρ_{b} c_{b} (T - T_{a}) + (q_{met} + τ_{q} \frac{\partial q_{met}}{\partial t}), (3)

式中: $T_{a}$ 为动脉血液温度; $ρ_{t} 和 c_{t}$ 分别为组织密度和比热容; $ρ_{b} 和 c_{b}$ 分别为血液密度和比热容; ${\tilde{w}}_{b}$ 为血液灌注率; $q_{met}$ 为新陈代谢热。

当忽略内部相互作用即 $τ_{T} = 0$ 时,得到单相滞生物传热方程

\begin{array}{l} ρ_{t} c_{t} τ_{q} \frac{\partial^{2} T}{\partial t^{2}} + (τ_{q} {\tilde{w}}_{b} ρ_{b} c_{b} + ρ_{t} c_{t}) \frac{\partial T}{\partial t} = \\ k \frac{\partial^{2} T}{\partial x^{2}} - {\tilde{w}}_{b} ρ_{b} c_{b} (T - T_{a}) + (q_{met} + τ_{q} \frac{\partial q_{met}}{\partial t}) \end{array} 。 (4)

进一步忽略热流快速加热效应,得到经典的傅里叶型生物传热方程

ρ_{t} c_{t} \frac{\partial T}{\partial t} = k \frac{\partial^{2} T}{\partial x^{2}} - {\tilde{w}}_{b} ρ_{b} c_{b} (T - T_{a}) + q_{met} 。 (5)

2.2 边界条件和初始条件

为了建立基于非傅里叶传热机制的Neumman边界条件,在图1所示的有限厚度组织中,考虑左侧Δx区域的能量平衡,其左侧和右侧的热流密度分别为 $J_{0} f (t) 和 q (Δ x, t),$ 且对于生物组织,还要兼顾血流灌注项和新陈代谢。即

\begin{array}{l} ρ_{t} c_{t} Δ x \frac{\partial T}{\partial t} = J_{0} f (t) - q (Δ x, t) + \\ ρ_{b} c_{b} {\tilde{w}}_{b} (T_{a} - T) Δ x + q_{met} Δx \end{array} 。 (6)

考虑DPL模型,则 $q (Δ x, t)$ 满足

- k [\frac{\partial T (Δ x, t)}{\partial x} + τ_{T} \frac{\partial^{2} T (Δ x, t)}{\partial^{2} x}] = q (Δ x, t) + τ_{q} \frac{\partial}{\partial t} q (Δ x, t) 。 (7)

当 $Δ x \to 0,$ (6)式和(7)式分别为

q (x, t) = J_{0} f (t), x = 0, (8)

\begin{array}{l} - k [\frac{\partial T (x, t)}{\partial x} + τ_{T} \frac{\partial^{2} T (x, t)}{\partial x \partial t}] = \\ q (x, t) + τ_{q} \frac{\partial}{\partial t} q (x, t), x = 0 \end{array} 。 (9)

将(8)式代入(9)式,获得

\begin{array}{l} - k [\frac{\partial T (x, t)}{\partial x} + τ_{T} \frac{\partial^{2} T (x, t)}{\partial x \partial t}] = \\ J_{0} [f (t) + τ_{q} \frac{\partial f}{\partial t}], x = 0 \end{array} 。 (10)

(10)式就是DPL非傅里叶模型对应的热流Neumman边界条件,当忽略内部相互作用即 $τ_{T} = 0$ 时,得到SPL非傅里叶模型对应的热流Neumman边界条件:

- k \frac{\partial T (x, t)}{\partial x} = J_{0} [f (t) + τ_{q} \frac{\partial f}{\partial t}], x = 0 。 (11)

当考虑傅里叶导热机制时,第二类热流Neumman边界条件为

- k \frac{\partial T}{\partial x} = J_{0} f (t), t > 0, x = 0 。 (12)

可以看出,(10)式和(11)式右边含有热流的松弛效应,(10)式左边还含有温度梯度相延迟项,它们都是基于能量方程和非傅里叶定律获得的,因此(10)式和(11)式是DPL和SPL非傅里叶边界能量守恒方程。(12)式是经典傅里叶边界能量守恒方程,该方程是基于微元体能量方程和傅里叶定律获得的,两者在数学表达式上是有差别的。

对应图1的右侧边界条件,采用的绝热边界条件为

\frac{\partial T (x, t)}{\partial x} = 0, t > 0, x = L 。 (13)

初始时刻认为组织完全被血流灌注,动脉血处于平衡态,且组织本身处于准静态。因此初始条件可以表示为

\{\begin{array}{l} T (x, t) = T_{a}, t = 0, 0 < x < L \\ \frac{\partial T (x, t)}{\partial t} = 0, t = 0, 0 < x < L \end{array}, (14)

3 温度场的求解

设 $q = J_{0} f (t), 其中 f (t) = U (t) - U (t - t_{0}), t_{0}$ 为脉冲时间, $U (t)$ 为单位阶跃函数。

为了求解方便,将SPL、DPL和Pennes生物热传导方程及其对应的边界条件进行无量纲处理。引入无量纲量

\{\begin{array}{l} x^{*} = \frac{x}{L} \\ t^{*} = \frac{at}{L^{2}} \\ T^{*} = \frac{k (T - T_{a})}{L J_{0}} \\ q_{met}^{*} = \frac{q_{met} L}{J_{0}} \end{array}, (15)

\{\begin{array}{l} {(τ_{T}, τ_{q}, t_{0})}^{*} = \frac{a}{L^{2}} (τ_{T}, τ_{q}, t_{0}) \\ {\tilde{w}}^{*}_{b} = \frac{{\tilde{w}}_{b} L^{2}}{a} \end{array}, (16)

式中: $a = k / (ρ_{t} c_{t}) 。$

3.1 非傅里叶和傅里叶边界条件DPL模型的解

考虑定常 $q_{met},$ 通过无量纲化,将上标星号略去,则得到DPL生物热传导方程为

\begin{array}{l} (1 + τ_{q} {\tilde{w}}_{b} β) \frac{\partial T}{\partial t} + τ_{q} \frac{\partial^{2} T}{\partial t^{2}} = \\ (1 + τ_{T} \frac{\partial}{\partial t}) \frac{\partial^{2} T}{\partial x^{2}} - {\tilde{w}}_{b} βT + q_{met} \end{array}, (17)

式中: $β = ρ_{b} c_{b} / ρ_{t} c_{t} 。$

加热侧非傅里叶边界条件为

- [\frac{\partial T}{\partial x} + τ_{T} \frac{\partial^{2} T}{\partial x \partial t}] = U (t) - U (t - t_{0}) + τ_{q} [δ (t) - δ (t - t_{0})], t > 0, x = 0, (18)

式中: $δ (\cdot)$ 为狄拉克函数。

加热侧傅里叶边界条件为

- \frac{\partial T}{\partial x} = U (t) - U (t - t_{0}), t > 0, x = 0 。 (19)

另一侧按绝热边界对待,有

\frac{\partial T (x, t)}{\partial x} = 0, t > 0, x = 1 。 (20)

初始条件

\{\begin{array}{l} T (x, t) = 0, t = 0,0 < x < 1 \\ \frac{\partial T (x, t)}{\partial t} = 0, t = 0,0 < x < 1 \end{array} 。 (21)

从(17)式和边界条件(18)~(20)式,运用有限余弦傅里叶变换和拉普拉斯变换进行求解,其最后的解为

T (x, t) = θ_{0} (t) + 2 \overset{\infty}{\sum_{n = 1}} θ_{n} (t) \cos (n π x) 。 (22)

非傅里叶边界条件时

\begin{array}{l} θ_{0} (t) = φ_{10} (t) - φ_{10} (t - t_{0}) + \\ \frac{τ_{q} q_{met}}{τ_{q} β {\tilde{w}}_{b} - 1} [1 - \exp (- \frac{t}{τ_{q}})] \end{array}, (23)

式中: $φ_{10} (t) = \frac{1}{β {\tilde{w}}_{b} - 1} [1 - \exp (- β {\tilde{w}}_{b} t)] 。$

傅里叶边界条件时

θ_{0} (t) = φ_{20} (t) - φ_{20} (t - t_{0}) + \frac{q_{met} τ_{q}}{1 - {\tilde{w}}_{b} β τ_{q}} \{\frac{1}{{\tilde{w}}_{b} β - 1} [1 - \exp (- β {\tilde{w}}_{b} t)] - \frac{τ_{q}}{1 - τ_{q}} [1 - \exp (- \frac{t}{τ_{q}})]\}, (24)

式中: $φ_{20} (t) = φ (t) + τ_{T} \overset{t}{\int_{0}} φ (ξ) d ξ$ ,其中

\begin{array}{l} φ (t) = \frac{1}{{\tilde{w}}_{b} β τ_{q} - 1 `} \{\frac{1}{1 - τ_{q}} [1 - \exp (- \frac{t}{τ_{q}})] + \frac{1}{τ_{T} {\tilde{w}}_{b} β - 1} [\exp (- {\tilde{w}}_{b} βt) - \exp (- \frac{t}{τ_{q}})]\} \\ θ_{n} (t) = φ_{1 n} (t) - φ_{1 n} (t - t_{0}) 。 \end{array} 。 (25)

当 ${(1 + β {\tilde{w}}_{b} τ_{q} + n^{2} π^{2} τ_{T})}^{2} > 4 τ_{q} (n^{2} π^{2} + {\tilde{w}}_{b} β)$ 时,有

\begin{matrix} φ_{1 n} (t) = \frac{2}{ε_{1}} [α_{i} \exp (- \frac{1 + β τ_{q} {\tilde{w}}_{b} + n^{2} π^{2} τ_{T}}{2 τ_{q}} t) \sinh (\frac{ε_{1}}{2 τ_{q}} t) + \overset{t}{\int_{0}} \exp (- \frac{1 + β τ_{q} {\tilde{w}}_{b} + n^{2} π^{2} τ_{T}}{2 τ_{q}} ξ) \sinh (\frac{ε_{1}}{2 τ_{q}} ξ) d ξ] \end{matrix}, (26)

当 ${(1 + β {\tilde{w}}_{b} τ_{q} + n^{2} π^{2} τ_{T})}^{2} < 4 τ_{q} (n^{2} π^{2} + {\tilde{w}}_{b} β)$ 时,有

\begin{matrix} φ_{1 n} (t) = \frac{2}{ε_{2}} [α_{i} \exp (- \frac{1 + β τ_{q} {\tilde{w}}_{b} + n^{2} π^{2} τ_{T}}{2 τ_{q}} t) \sin (\frac{ε_{2}}{2 τ_{q}} t) + \overset{t}{\int_{0}} \exp (- \frac{1 + β τ_{q} {\tilde{w}}_{b} + n^{2} π^{2} τ_{T}}{2 τ_{q}} ξ) \sin (\frac{ε_{2}}{2 τ_{q}} ξ) d ξ] \end{matrix}, (27)

式中:ε₁= $\sqrt[]{{(1 + {\tilde{w}}_{b} β τ_{q} + n^{2} π^{2} τ_{T})}^{2} - 4 τ_{q} (n^{2} π^{2} + {\tilde{w}}_{b} β)}$ ;ε₂= $\sqrt[]{4 τ_{q} (n^{2} π^{2} + \tilde{w} β) - {(1 + {\tilde{w}}_{b} β τ_{q} + n^{2} π^{2} τ_{T})}^{2}}$ ;i=1对应非傅里叶边界且 $α_{1} = τ_{q}; i = 2$ 对应傅里叶边界且 $α_{2} = τ_{T} 。$

3.2 Pennes模型的解

令(17)~(18)式中 $τ_{q} = 0$ 以及 $τ_{T} = 0,$ 则获得Pennes生物传热方程及相应的初始条件和边界条件,其解的形式与(22)式相同,其中

\begin{array}{l} θ_{0} (t) = φ_{30} (t) - φ_{30} (t - t_{0}) + \\ \frac{q_{met}}{β {\tilde{w}}_{b} - 1} [1 - \exp (- β {\tilde{w}}_{b} t)] \end{array}, (28)

θ_{n} (t) = φ_{3 n} (t) - φ_{3 n} (t - t_{0}), (29)

式中: $φ_{30} (t) = 1 - \exp (- β {\tilde{w}}_{b} t)$ ; $φ_{3 n} (t) = 1 - \exp [- (β {\tilde{w}}_{b} + n^{2} π^{2}) t] 。$

4 算例和讨论

为了讨论非傅里叶Neumann边界条件对生物组织内传热的影响,选取文献[ 21]中的皮肤作为算例对象,皮肤近似看作单层模型,厚度 $L$ =0.006m, $q_{met} = 368.1 W \cdot m^{- 3},$ 其热物性如表1所示。以下计算中取 $J_{0} = 1000 W \cdot m^{- 2}, t_{0} = 20 s, τ_{q} = 16 s, τ_{T} = 0,4, 8 s 。$

图2显示了傅里叶和非傅里叶Neumman边界条件下双相滞模型的温度分布对比,其中Pennes模型对应的方程是傅里叶热传导方程,边界条件是傅里叶边界条件。在图2(a)中,Pennes模型对应于最高温度和最大影响区域。 $τ_{q} = 16 s, τ_{T} = 0$ 时DPL模型退化为SPL热波模型,其温度最低,影响区域最小。而对于处于二者之间的DPL热传导方程模型,温度分布区域和高低则与两个弛豫时间有关,这里统一取 $τ_{q} = 16 s, τ_{T} = 0,4, 8 s 。$ 图2(a)中松弛时间 $τ_{T}$ 越长,温度越高,影响区域越大。然而,在非傅里叶边界条件对应的图2(b)中, $τ_{q} 和 τ_{T}$ 的取值与图2(a)中相同,但温度分布较复杂。对于 $0 < x < 0.6$ 区域,Pennes模型的温度分布是最低的,热波模型预测的温度却是最高的,DPL模型的温度分布处于二者之间。对于DPL模型, $τ_{T}$ 越大,温度越低,这与图2(a)相反。在 $x > 0.6$ 区域,Pennes模型的温度最高,热波模型的温度幅值最低,DPL非傅里叶模型居于两者之间。

Figure 2.Temperature distributions under different boundary conditions when t=10s and $τ_{q}$ =16s. (a) Fourier boundary conditions; (b) non-Fourier boundary conditions

Parameter	Tissuedensity /(kg·m^-3)	Blooddensity /(kg·m^-3)	Blood specificcapacity /(J·kg^-1·℃)	Tissue specificcapacity /( J·kg^-1·℃)	Conductivecoefficient /(W·m^-1·℃)	Bloodperfusionrate /s^-1
Value	1190	3600	1060	3770	0.235	1.33×10^-4

Table 1. Thermo-physical parameters of skin

View all Tables

可以看出,傅里叶边界条件和非傅里叶边界条件对生物体内的温度分布的影响明显不一样,究竟哪种温度分布更合理呢?从热波传播速度来看,热传播速度越快,温度应该越低;相反,热传播速度越慢,温度应该越高。其次,从能量守恒定律分析,在0~10s内,同一时刻接收到的能量应该是一样的,包括Pennes模型在内的各条曲线与横轴围成的面积应该是一样的。但从图2(a)可以发现,各条温度曲线围成的面积显然不同。图2(b)却符合能量守恒定律。图2(a)、(b)中Pennes模型对应的面积是一样的。在1.2节导出边界条件时,DPL模型对应的是DPL非傅里叶边界条件,傅里叶模型对应的是傅里叶边界条件。傅里叶边界条件使得DPL模型对应的定解方程不满足能量守恒定律,由此可见,傅里叶边界条件用于非傅里叶方程时部分能量被傅里叶边界条件过滤了,这也是造成图2(a)所示现象的原因。图2(b)是非傅里叶边界条件用于DPL非傅里叶方程时得出的分布曲线,DPL非傅里叶边界条件使得DPL模型对应的定解方程满足能量守恒定律。

为了进一步说明非傅里叶边界条件与经典傅里叶边界条件所带来的预测温度的差别,图3给出了表面温度随时间的变化情况。图3(a)是傅里叶边界条件下计算出的温度响应,在 $0 < t < 20 s$ 区间,Pennes热传导方程( $τ_{q} = 0, τ_{T} = 0$ )预测的温度幅值最高,温度变化速度最快,热波方程( $τ_{q} = 16 s, τ_{T} = 0$ )预测的温度幅值最小,温度变化速度最慢。这里得出的结论与文献[ 10]~[13]的结论几乎是一致的,与[14]和[16]中的结论是相反的。显然其中所给出的边界条件是有缺陷的,也是不合理的。另外,从热传播速度分析,热波方程含有最慢的速度,能量积累应该较多,温度应该较高。而图3(a)所示的结论相反,这是傅里叶边界条件用于双相滞后非傅里叶热传导方程引起的。图3(b)是非傅里叶边界条件下计算出的温度响应。在 $0 < t < 20 s$ 区间,热波方程对应的温度变化速度最快,升温幅度也最大。Pennes方程预测的温度最低,温度变化速度最慢,刚好与图3(a)中的结论相反,这是由于非傅里叶延迟效应减小了能量的传递速度,能量出现累积,不同延迟时间表征不同的累积效应,因此温升幅度不一样。这里傅里叶边界条件下没有出现能量累积,从而表面的温升幅度最小,温升速率最慢。这个结论再次揭示了非傅里叶传热方程必须结合非傅里叶边界条件。

Figure 3.Temperature variations at x=0 under different boundary conditions when τ_q=16s. (a) Fourier boundary conditions; (b) non-Fourier boundary conditions

图4进一步给出了两类边界条件下x=1.2mm处的温度响应对比。图4(a)是傅里叶边界条件下多种传热模型的温度响应。可以看出,各种模型温度开始变化的时间起点是不同的。热波模型变化的滞后时间最大,含有无穷传播速度的Pennes模型几乎没有滞后时间, $τ_{q} = 16 s, τ_{T} = 4 s$ 时的温度变化曲线和 $τ_{q} = 16 s, τ_{T} = 8 s$ 时的温度变化曲线处于热波模型和Pennes模型之间,但 $τ_{T} = 8 s$ 时变化的滞后时间要小于 $τ_{T} = 4 s 。$ 图4(b)是非傅里叶边界条件下的温度响应。从变化的延迟时间来看,图4(a)与图4(b)是一样的。但从温度的变化速度和温升幅值来看,图4(a)与图4(b)明显不同。图4(b)所示的温度变化幅值要大于图4(a),图4(b)的温升速度也要高于图4(a)。这个结论与Mitra等^[23]的结论吻合,且文献[ 24]利用分数阶微分进一步描述了实验数据,即本文提出的边界条件得到了实验验证。

Figure 4.Temperature variations at x=1.2mm under different boundary conditions when $τ_{q}$ =16s. (a) Fourier boundary conditions; (b) non-Fourier boundary conditions

在图2(b)、图3(b)和图4(b)中, $τ_{T} = 0$ 对应的是SPL热波模型,其温度分布和温度变化速率在波阵面上呈现陡峭特点,该特点是修正的非傅里叶边界条件带来的,而 $τ_{T} \neq 0$ 的DPL热传导模型则没有该特点。通过控制(17)式中扩散项的时间空间混合导数项 $τ_{T} \frac{\partial}{\partial t} (\frac{\partial^{2} T}{\partial x^{2}}),$ 可直接破坏波动方程的结构,也就使热传导机制复杂化,当 $τ_{T} < τ_{q}$ 时,热传导具有扩散和波动的特性,当 $τ_{T} > τ_{q}$ 时热传导仅具有扩散特性^[22]。

在生物医学工程中,热损伤是激光手术的重要参数。一般认为,热损伤始于基底层温度,表皮与真皮的界面温度上升到44℃时热损伤就开始发生^[25]。Henriques等^[26]提出了评估热损伤最常用的数学公式:

Ω = \overset{t}{\int_{0}} A \exp (- \frac{E_{a}}{RT}) d t, (30)

式中:频率因子 $A = 3.1 \times 10^{9}; E_{a} 和 R$ 分别为活化能和理想气体常数,且 $E_{a}$ =211493324+2688.367ln A。一般 $Ω = 0.53$ 定义为一度损伤, $Ω = 1$ 是二度损伤, $Ω = 10^{4}$ 是三度损伤。活化能 $E_{a}$ 和频率因子 $A$ 是不随温度变化的常数。在实际中观测到,40~41℃时细胞会失去活性,但这是暂时的,一段时间后又会恢复活性且细胞还具有耐热性。如果稳态温度保持在42 ℃阈值处,热损伤与激光照射时间成正比。

图5给出了两类不同边界条件下x=0和x=1.2mm处的热损伤变化对比。在图5(a)中,热波模型的表面热损伤幅是最小的,达到一度、二度和三度损伤的速度是最慢的,而Pennes模型达到一度、二度和三度损伤的速度是最快的。这一点与图3(a)的预测结果是一致的,也是不合理的结果。这里的规律与文献[ 10]和[13]中描述的规律也是一致的,原因在于文献[ 10]和[13]采用了傅里叶边界条件。在图5(b)中,最早达到一度、二度和三度损伤的是热波模型,最慢达到的是Pennes模型,这个规律与文献[ 18]中描述的规律是一致的。由此可知,非傅里叶模型在傅里叶边界条件下不但预测的组织温度与物理原理是矛盾的,而且热损伤预测结论也是错误的,这给实际热疗带来错误的认识。图5(c)与图5(d)展示了内部位置(x=1.2mm)在两类边界条件下的热损伤变化。相比图5(a)、(b),Pennes模型的热损伤只是幅度衰减,而热波模型不但在幅度上有衰减,而且在时间上也有延迟,DPL模型预测的结果介于二者之间。但非傅里叶边界条件下预测的结果和傅里叶边界条件下预测的结果是明显不同的,甚至是相反的。

Figure 5.Thermal damage variations at x=0 and x=1.2mm under different boundary conditions when τ_q=16s. (a) Fourier boundary conditions at x=0; (b) non-Fourier boundary conditions at x=0; (c) Fourier boundary conditions at x=1.2mm; (d) non-Fourier boundary conditions at x=1.2mm

5 结论

针对激光加热中生物传热的非傅里叶效应,推导了非傅里叶生物传热方程,建立了相应的非傅里叶Neumann边界条件。通过傅里叶变换和拉普拉斯变换得到解析解,对傅里叶和非傅里叶Neumann边界条件下的非傅里叶生物热方程的解析解进行了比较和讨论,得出以下结论:

1)边界条件本质上是描述边界上的热平衡方程,非傅里叶传热方程对应的边界条件应该是非傅里叶边界条件,傅里叶传热方程对应的边界条件应该是傅里叶边界条件。

2)边界条件和控制方程应是相容的。对于非傅里叶传热控制方程,采用傅里叶边界条件时,温度分布和温度变化会出现与物理学原理相违背的情况,而非傅里叶传热方程与非傅里叶边界条件相结合,所得结论与物理学原理一致。

3)利用非傅里叶Neumann边界条件,非傅里叶传热方程预测的温度和热损伤值均分别高于Pennes预测的温度和热损伤值,且所得结论与傅里叶边界条件下所得结论是相反的,这一点是对一些现有文献所得结论的纠正,在以后的使用中应予以重视。

References

[1] Chen X, Lu J L, Li P. Viscoelasticity measurement of biological tissues using laser speckle techniques: a review[J]. Chinese Journal of Lasers, 45, 0207005(2018).

[2] Fu Z L, Li R Z, Li H Y et al. Progress in biomedical imaging based on terahertz quantum cascade lasers[J]. Chinese Journal of Lasers, 47, 0207014(2020).

[3] Pennes H H. Analysis of tissue and arterial blood temperatures in the resting human forearm[J]. Journal of Applied Physiology, 85, 5-34(1948).

[4] Liu J, Ren Z P, Wang C C[J]. Interpretation of living tissue's temperature oscillations by thermal wave theory Chinese Science Bulletin, 1995, 1493-1495.

[5] Xu F, Lu T J. Introduction to skin biothermomechanics and thermal pain[M]. Heidelberg: Springer(2011).

[6] Antaki P J. New interpretation of non-Fourier heat conduction in processed meat[J]. Journal of Heat Transfer, 127, 189-193(2005).

[7] Xu F, Lu T J, Seffen K A. Dual-phase-lag model of skin bioheat transfer. [C]//2008 International Conference on BioMedical Engineering and Informatics, May 27-30, 2008, Sanya, China. New York: IEEE, 10029091(2008).

[8] Xu G Y, Wang J B, Han Z. Study on the transient temperature field based on the fractional heat conduction equation for laser heating[J]. Applied Mathematics and Mechanics, 36, 844-854(2015).

[9] Wang Y Z, Wang Q, Liu D et al. Generalized thermoelastic solutions to the problems of thermal shock on elastic half space[J]. Applied Mathematics and Mechanics, 35, 640-651(2014).

[10] Yang H Q, Chen J L, Wang Y H et al. Fourier and non-Fourier heat conduction effects in biological tissue under laser irradiation[J]. Chinese Journal of Lasers, 36, 2582-2586(2009).

[11] Kumar D, Singh S, Rai K N. Analysis of classical Fourier, SPL and DPL heat transfer model in biological tissues in presence of metabolic and external heat source[J]. Heat and Mass Transfer, 52, 1089-1107(2016).

[12] Ziaei P, Moosavi H, Moradi A. Analysis of the dual phase lag bio-heat transfer equation with constant and time-dependent heat flux conditions on skin surface[J]. Thermal Science, 20, 1457-1472(2016).

[13] Ahmadikia H, Fazlali R, Moradi A. Analytical solution of the parabolic and hyperbolic heat transfer equations with constant and transient heat flux conditions on skin tissue[J]. International Communications in Heat and Mass Transfer, 39, 121-130(2012).

[14] Maurer M J, Thompson H A. Non-Fourier effects at high heat flux[J]. Journal of Heat Transfer, 95, 284-286(1973).

[15] Al-Nimr M A, Khadrawi A F, Hammad M. A generalized thermal boundary condition for the hyperbolic heat conduction model[J]. Heat and Mass Transfer, 39, 69-79(2002).

[16] Antaki P J. Analysis of hyperbolic heat conduction in a semi-infinite slab with surface convection[J]. International Journal of Heat and Mass Transfer, 40, 3247-3250(1997).

[17] Shirmohammadi R, Moosaie A. Non-Fourier heat conduction in a hollow sphere with periodic surface heat flux[J]. International Communications in Heat and Mass Transfer, 36, 827-833(2009).

[18] Biswas P, Singh S, Bindra H. A closed form solution of dual-phase lag heat conduction problem with time periodic boundary conditions[J]. Journal of Heat Transfer, 141, 031302(2019). http://www.researchgate.net/publication/330346344_A_Closed_Form_Solution_of_Dual-Phase_Lag_Heat_Conduction_Problem_With_Time_Periodic_Boundary_Conditions

[19] Forghani P, Ahmadikia H, Karimipour A. Non‐Fourier boundary conditions effects on the skin tissue temperature response[J]. Heat Transfer Research, 46, 29-48(2017).

[20] Kumar P, Kumar D, Rai K N. Numerical simulation of dual-phase-lag bioheat transfer model during thermal therapy[J]. Mathematical Biosciences, 281, 82-91(2016).

[21] Wang J B, Wang J, Tian M L et al. Study on multi-pulse heat flow heating biological tissue based on thermal wave model[J]. Journal of Engineering Thermophysics, 40, 1904-1912(2019).

[22] Tzou D Y. Macro- to microscale heat transfer: the lagging behavior[M]. 2nd ed. Chichester :John Wiley & Sons, Ltd.(2015).

[23] Mitra K, Kumar S, Vedevarz A et al. Experimental evidence of hyperbolic heat conduction in processed meat[J]. Journal of Heat Transfer, 117, 568-573(1995). http://scitation.aip.org/getabs/servlet/GetabsServlet?prog=normal&id=JHTRAO000117000003000568000001&idtype=cvips&gifs=Yes

[24] Xu H Y. Research on the application of fractional calculus in anomalous transportation[D]. Jinan: Shandong University(2017).

[25] Torvi D A, Dale J D. A finite element model of skin subjected to a flash fire[J]. Journal of Biomechanical Engineering, 116, 250-255(1994).

[26] Henriques F, Montz A. Studies of thermal injury: I. The conduction of heat to and through skin and the temperature attained therein. A theoretical and experimental investigation[J]. American Journal of Pathology, 23, 530-549(1947).