使用 COMSOL Multiphysics® 模拟超快传热

2026年 1月 22日

超快脉冲激光技术的进步使得对极短时间尺度上的热传递的理解变得尤为重要。 在亚纳秒尺度下,传统的传热方程失效,必须采用更精细的模型,例如双曲线传热方程和双温度模型。本文,我们将简要回顾最常用的非傅里叶传热模型的相关理论,并探讨如何在 COMSOL Multiphysics® 软件中实现这些模型。

重新审视传热方程

首先,我们将简要回顾科学工程领域最常见的偏微分方程之一:固体传热方程。其起点是能量方程,或者说是温度连续性方程 T

C \frac{\partial T}{\partial t} = -\nabla \cdot \mathbf{q} + S

 

其中 C=\rho C_p 为体积热容,\mathbf{q} 为热通量, S 为体积热源。

该方程本质上是能量守恒表达,因此无论时空尺度如何变化都必须成立。为获得 T的闭合形式偏微分方程,需先确定\mathbf{q} 的表达式。最简单直观的选择是傅里叶定律,该定律指出热通量与温度梯度成正比:

\mathbf{q} = -k\nabla T

 

其中比例常数 k 称为热导率。将其代入连续性方程即可得到经典的热传导方程:

C \frac{\partial T}{\partial t} = \nabla \cdot(k\nabla T) + S

 

这是一个抛物型偏微分方程,故称为抛物型单步(POS)模型。(”单步”的含义将在本文后续内容中阐明。)

通过考察半无限介质域(x \geq 0 )可深入了解该方程的特性:当x=0 表面温度保持为 T_0 ,而域内其余区域初始温度为 T_i 时,该问题存在如下解析解:

T(x,t)=T_0 + (T_i-T_0)\, \mathrm{erf} \frac{x}{2\sqrt{\alpha t}}

 

其中,\alpha = k/C 为介质的热扩散率,\mathrm{erf} 表示误差函数。该解的动画演示见下图。由上述表达式可知,边界条件 x=0 的影响会立即传递至整个区域,这违反了狭义相对论。此现象源于推导该方程时采用的强假设——傅里叶定律。

尽管名为傅里叶定律,其实质是一种唯象关系,其从第一性原理的严谨推导涉及多重近似。但需注意的是:在常规的时间与长度尺度下,能量的空间传输速度足够快,可视为瞬时过程,因此对于大多数工程应用而言,热方程的这一缺陷在多数工程应用中可安全地忽略不计。

热方程解在简单半无限域中的时域演化动画。

超越傅里叶定律

当我们转向更短的时间尺度时,无法再像前文所述那样假设能量的空间传输是瞬时的。此外,还必须考虑电子与晶格之间能量传递的有限速度。下文将阐明前者如何导致双曲型热传导,后者又如何催生双温模型(参考文献1)。

双曲型传热

我们不再假设热通量能瞬时响应温度梯度,而是提出如下”滞后型”傅里叶定律:

\mathbf{q}(\mathbf{r}, t+\tau) = -k\nabla T(\mathbf{r}, t)

 

滞后时间 \tau 应与电子弛豫时间(即电子平均自由时间)相对应。如果滞后时间足够短时,上述表达式可展开为泰勒级数:

\mathbf{q}(\mathbf{r}, t) +\tau \frac{\partial}{\partial t}\mathbf{q}(\mathbf{r}, t) = -k\nabla T(\mathbf{r}, t)

 

该表达式被称为 Cattaneo–Vernotte 模型(参考文献2)。将其代入连续性方程(在对后者进行时间微分后)可得到双曲型一步(HOS)模型:

C \frac{\partial T}{\partial t} + \tau C \frac{\partial ^2 T}{\partial t^2} = \nabla \cdot (k \nabla T )+ S +\tau \frac{\partial S}{\partial t}

 

该双曲型偏微分方程与波动方程相似,而且正如我们将在下文进行阐述的,热传导的有限速度确实会导致”热波”现象的出现。

双温度模型

若摒弃电子与晶格始终处于局部热平衡的这一假设,则需引入两个独立的温度场(晶格温度 T_l 与电子温度 T_e )及其耦合关系项 G 。初始阶段,激光能量被电子吸收导致 T_e 急剧上升(步骤1),随后与 T_l 发生耦合(步骤2)。 在亚纳秒时间尺度下,热量通过晶格的扩散相对于电子输运可忽略不计,因此最终得到以下控制方程:

\begin{align*}
C_e \frac{\partial T_e}{\partial t} &= \nabla \cdot (k_e \nabla T_e ) -G(T_e-T_l) + S \\
C_l \frac{\partial T_l}{\partial t} &= G(T_e-T_l)
\end{align*}

 

该模型称为抛物型双步(PTS)模型。虽然可消去 T_e 以获得 T_l 的闭合形式偏微分方程,但会产生繁琐的交叉项。保留这种耦合的偏微分方程与域常微分方程的形式更为便捷。电子-晶格耦合参数 G 的估计,可通过对所有可能跃迁的电子-声子能量交换率进行积分获得(参考文献3)。

若考虑电子热能传播的有限速度,可进一步建立双曲型双步(HTS)模型。与上述模型类似,当忽略晶格扩散时,该模型可表示为:

\begin{align*}
C_e \frac{\partial T_e}{\partial t} + \tau_e C_e \frac{\partial ^2 T_e}{\partial t^2} &= \nabla \cdot (k_e \nabla T_e ) -G(T_e-T_l) -\tau_e \frac{\partial }{\partial t}[G(T_e-T_l)]+ S +\tau_e \frac{\partial S}{\partial t}\\
\quad C_l \frac{\partial T_l}{\partial t} &= G(T_e-T_l)
\end{align*}

 

进一步推广

为获得精确结果,我们可能还需要考虑以下若干附加效应,例如:

  • 除热通量外,引入温度滞后时间,由此形成双相位滞后(DPL)模型
  • 在上述双温度模型中纳入晶格扩散效应,由此衍生出双抛物型双步模型(DPTS)与双曲型双步(DHTS)模型
  • 采用双温度漂移-扩散型公式显式建模载流子动力学,即半经典双步(STS)模型
  • 考虑辐射传热效应

虽然这超出了本文的讨论范围,但对于许多应用而言,考虑超快加热过程中可能发生的更复杂的过程(如烧蚀和熔化)至关重要。

在 COMSOL Multiphysics® 中实现非傅里叶传热

默认情况下, COMSOL Multiphysics® 的传热接口仅支持 POS 模型,但借助偏微分方程(PDE)和常微分方程(ODE)接口,可在软件中轻松实现双曲型热传导与双温度模型。

以示例问题为例:考虑一层厚度为200 nm 的金膜,受到一个激光脉冲的照射,其吸收的激光脉冲通量 \Phi _p = 1\,\mathrm{J}/\mathrm{m}^2 (此处假设与其电子温度和晶格温度无关),脉冲具有高斯时间分布,脉宽 t_p = 1 \, \mathrm{ns}, \; 100 \, \mathrm{ps} , \; 100 \, \mathrm{fs},10 \, \mathrm{fs}\Phi _p 包含激光波长信息,因其同时取决于薄膜反射率。 我们还假设光束半径远大于其穿透深度 d ,从而使问题简化为一维模型。因此,由该脉冲产生的体积热源可表示为(参考文献1):

S=I/d=\sqrt{\frac{4\ln 2}{\pi}}\frac{\Phi _p}{t_pd}e^{-x/d}e^{-4\ln 2 \left(\frac{t}{t_p}\right)^2}

 

前置因子确保了总吸收脉冲能量的正确归一化。该脉冲中心位于 t=0,故模拟起始点设为 t=-2t_p 以捕捉完整的脉冲。我们采用参考文献3中报道的金的材料参数,并假设其保持恒定。

双曲型传热

为在 HOS 模型中包含时间的二阶导数项,需要使用 一般形式偏微分方程 接口。相关设置如下所示。

软件中 HOS 模型设置界面的截图。 在 COMSOL Multiphysics® 中实现 HOS 模型。为了包含二阶时间导数项,需使用一般形式偏微分方程接口。

双温度模型

为了实现 PTS 模型,可以使用预定义的传热 接口模拟电子温度。由于该模型忽略晶格扩散,仅需一个域常微分方程即可求解晶格温度。所需设置如下图所示。

在 COMSOL Multiphysics® 中实现 PTS 模型。对于电子温度(左侧),可使用默认的传热接口,并添加因耦合引起的额外源项;对于晶格温度(右侧),则使用分布式常微分方程接口,并将耦合项指定为源项。

与上述设置类似,我们也可以通过替换电子温度的传热 接口为一般形式偏微分方程 接口来实现HTS模型。需要注意的是:本模型中热导率被视为恒定标量,但传热 接口与一般形式偏微分方程 接口均支持例如温度依赖的系数,以及用于模拟各向异性传热的张量形式。

结果与讨论

让我们首先来看下 t_p = 1 \, \mathrm{ns} 时的计算结果。下图展示了辐照后晶格温度的空间分布。可以看到,在该脉冲持续时间下,正如预期的那样,不同建模方案的结果基本一致。

四种模型的晶格温度分布曲线图分别用浅蓝色、深蓝色、绿色和红色表示,在该时间尺度下,它们的结果几乎相同。 脉冲持续时间为 1 纳秒时晶格温度随深度变化的曲线。在此时间尺度下,四种模型的预测结果高度一致。

接下来,我们观察 t_p = 100 \, \mathrm{ps} 时的同类图谱(见下图)。 单步模型(POS 和 HOS)与双步模型(PTS 和 HTS)之间出现了明显差异:由于忽略了相对较快的电子热传递,单步模型预测的辐照表面附近晶格温度更高。然而,抛物型模型(POS 和 PTS)与双曲型模型(HOS 和 HTS)之间尚无明显差异。

脉冲持续时间为 100 皮秒时的晶格温度分布曲线图,其中四种模型分别用浅蓝色、深蓝色、绿色和红色表示,结果各不相同。 脉冲持续时间为 100 皮秒时晶格温度随深度变化的曲线。单步模型(POS 和 HOS)预测在辐照表面附近的晶格温度更高。

当脉冲进入飞秒尺度( t_p = 100 \, \mathrm{fs})时,这种差异变得更为显著,如下图所示。

脉冲持续时间为 100 飞秒时的晶格温度分布曲线图,显示不同颜色的模型得出了不同的结果。 脉冲持续时间为 100 飞秒时晶格温度随深度变化曲线。双曲型模型(HOS 和 HTS)中包含的二阶时间导数的影响开始显现。

此时,单步模型的预测结果要高出数百度,而双曲型模型因电子热传导存在有限速度,在辐照表面附近呈现出高于对应抛物型模型的温度值。

到目前为止,我们只考虑了晶格温度。若将电子温度绘制在同一坐标系(采用不同刻度),可直观展现电子温度与晶格温度的差异。结果如下方动画所示。

动画展示了 100 飞秒脉冲作用下电子温度(左)与晶格温度(右)随时间演变的分布。

接下来,为了进一步阐明包含二阶时间导数的影响,让我们比较双步模型对 t_p = 10 \, \mathrm{fs} 条件下预测的电子温度分布。现在可以清晰地观察到热波在薄膜中传播的现象,如下图所示。需要注意的是,在该脉冲长度下,电子温度达到 了足以需要考虑辐射传热的数值,而当前模型中忽略了这一效应。

一张电子温度图,其中五种不同的模型用不同的颜色表示。 脉冲持续时间为 10 飞秒时,电子温度随深度变化的曲线图。可以观察到由HTS模型预测的热波传播现象。

最后,为了了解薄膜中热传递的动态过程,可绘制前表面(受辐照表面)与后表面温度随时间的变化曲线。但因尺度差异较大,必须将每组数据按其对应的最高温度进行归一化处理。

下图展示了脉冲持续时间为 100 飞秒时的结果。在此时间尺度下,抛物型模型与双曲型模型的差异尚可接受,但单步模型无法预测后表面温度的定性行为。 单步模型显示温度单调递增,而双步模型预测后表面温度会达到峰值——这是由于电子-晶格弛豫的过程比电子传热更缓慢所致。这些结果与参考文献 3 中的图 3 和图 4 高度吻合。

归一化电子温度的曲线图,其中以蓝色和绿色展示了各种不同的结果。 脉冲持续时间为100飞秒时,薄膜前表面与后表面归一化电子温度(双步模型)及晶格温度(单步模型)的时间演化曲线。

结论

在本篇博客中,我们通过实现四种简单而重要的非傅里叶传热模型——抛物型单步模型(POS)、双曲型单步模型(HOS)、抛物型双步模型(PTS)和双曲型双步模型(HTS)——展示了如何利用 COMSOL Multiphysics® 模拟超快传热过程。

下一步

下载本文展示的示例模型文件:

参考文献

  1. V.E. Alexopoulou and A.P. Markopoulos, “A Critical Assessment Regarding Two-Temperature Models: An Investigation of the Different Forms of Two-Temperature Models, the Various Ultrashort Pulsed Laser Models and Computational Methods,” Arch Computat Methods Eng, vol. 31, pp. 93–123, 2024. https://doi.org/10.1007/s11831-023-09974-1
  2. Y. Zhang, D.Y. Tzou, and J.K. Chen, “Micro- and Nanoscale Heat Transfer in Femtosecond Laser Processing of Metals”, High-Power and Femtosecond Lasers: Properties, Materials and Applications, eds. P.H. Barret and M. Palmerm, Nova Science Publishers, Inc., Hauppauge, NY, ch. 5, pp. 159–206, 2009. https://arxiv.org/abs/1511.03566
  3. T.Q. Qiu and C.L. Tien, “Heat Transfer Mechanisms During Short-Pulse Laser Heating of Metals.” ASME J. Heat Transfer, vol. 115, iss. 4, pp. 835–841, 1993. https://doi.org/10.1115/1.2911377

评论 (0)

正在加载...
浏览 COMSOL 博客