☰
金纳米颗粒光热效应的COMSOL仿真:从局域等离激元到稳态温度场
2026/10/14 7:26:20 网站建设 项目流程

金纳米颗粒的光热效应,这几年在纳米医学、等离激元传感、光热抗菌、光声成像领域被反复提到,但真要把“光照一颗40纳米的金球,它周围的水到底热了多少度”这件事算清楚,很多同学第一反应是查文献,第二反应是打开COMSOL,然后卡在“介电函数怎么给、边界条件怎么设、网格怎么剖、吸收截面怎么提”这一连串环节里。这篇博文就是来打通这条路的:从局域表面等离激元共振的物理机制,到麦克斯韦方程和热传导方程的控制方程,再到COMSOL Multiphysics里一步一步建模、求解、后处理,最后给出一套可以直接照抄的基准算例参数和脚本化复现思路。适合正在做纳米光热计算、准备把实验和仿真对照、或者刚接触COMSOL电磁-热耦合的同学,读完之后你能独立把单个金纳米颗粒的光热仿真完整跑通,并且看懂每一步在干什么。

1. 先弄清金纳米颗粒为什么会发热

1.1 局域等离激元共振:自由电子在纳米尺度内的“集体振荡”

金纳米颗粒之所以对光热转换特别敏感,根源在于局域表面等离激元共振。把一颗直径几十纳米的球形金颗粒放在水里,入射光打到它表面时,颗粒内部的自由电子气会在电场驱动下做集体振荡。由于颗粒尺寸远小于入射光波长,颗粒的曲率边界会把这种振荡“憋”在纳米尺度范围内,形成与颗粒几何、周围介质折射率强相关的共振模式。用个不太严谨但好记的类比:一盆水端在手里,手一抖,水会沿着盆壁周期性晃荡;颗粒边界就相当于盆壁,入射光就相当于手抖的激励,震荡幅度最大的那个频率就是共振频率。

在共振波长处,金颗粒对光的吸收和散射都会急剧增强,尤其是吸收,因为金的介电常数虚部在可见光到近红外区恰好比较大。这就是为什么金纳米颗粒被称为“光热转换的微型锅炉”:它把光能先变成电子气振荡的动能,再通过电子-声子碰撞把能量交给晶格,最后传导成周围介质的宏观温升。这个机制本身并不复杂,但要在仿真里精确复现,核心就是抓住两个量:颗粒在入射波长处的介电函数,以及颗粒周围的电磁场分布。

1.2 从光到热的三个时间尺度,决定我们用稳态还是瞬态

光热过程不是瞬时就完成的热扩散,它跨越了三个明显不同的时间尺度。第一阶段是自由电子的受激振荡和电子气内部热化,发生在飞秒到皮秒量级,电子被“加热”后的能量分布迅速通过电子-电子散射达到费米分布;第二阶段是电子-声子耦合,电子把能量交给晶格,典型时间是几皮秒到几十皮秒,也就是颗粒自身“体温”上升的阶段;第三阶段才是颗粒通过热传导把多余热量传递给周围水、形成微米尺度温度场的过程,这个时间取决于颗粒尺寸和介质热扩散率,通常在纳秒到微秒量级。

阶段时间尺度物理机制在仿真中的处理
电子热化0.1~1 ps电子-电子散射不建模,统计平均即可
电子-声子耦合1~100 ps晶格升温不建模,共振吸收已包含能量转换
热扩散到介质ns~μs颗粒向水导热的堆积连续光时可以直接算稳态温度场

对连续激光照射的典型金纳米颗粒光热实验,前两个阶段持续极短,对宏观温升没有持久影响,所以我们在仿真里可以忽略时间演化细节,直接解频域电磁场、再把吸收功率作为热源交给稳态热传导。只有研究脉冲激光激发、飞秒时间分辨光谱时才需要进入瞬态传热求解,且还要引入电子-声子耦合的双温度模型。这篇文章覆盖的连续波场景,稳态模型已经足够。

1.3 单颗粒仿真为什么不能省

实验上测单颗金纳米颗粒的绝对温升非常麻烦:远场光学方法给的是系综平均,近场测温手段(比如掺杂荧光分子的温度探针)又受空间分辨率和荧光淬灭限制。仿真则可以直接给出颗粒表面温度、周围介质的热流密度、温度梯度等实验难以直接测得的细节。这些细节才是光热治疗剂量控制、光热流控、纳米颗粒热泳等应用的底层数据。对单个颗粒的模拟也是一切复杂系统的基础:你连单颗粒的吸收截面和温升都算不准,后面做颗粒二聚体、颗粒阵列、细胞包裹颗粒的模型,误差只会被放大。

2. 控制方程和模型参数:把物理问题变成可算的数学问题

2.1 金的介电函数是整个仿真的“心脏”

光热仿真最容易出错、也最决定结果质量的地方,就是金的介电函数。在COMSOL里,光学仿真用到的金属材料一般不是简单地填一个常数电导率,而是需要给出随波长(或频率)变化的复介电常数。常见方式有几种:直接使用光学手册的实验数据插值表(比如对金、银、铝从紫外到近红外的复折射率n和k测量值),或者用自由电子气模型自行构建。实验数据的好处是准,坏处是版本和数据来源差异会导致结果有百分之几的差别,不同文献使用的数据点也不一样,只要保证峰位偏移在可接受范围内就可以。

自由电子气模型的好处是参数少、形式清楚,可以方便地加入尺寸修正。对孤立纳米颗粒,当颗粒半径接近电子平均自由程时,电子与颗粒边界的碰撞会显著增加阻尼,等效于在费米速度基础上增加一个散射频率,这会导致共振峰展宽、吸收效率降低。在仿真中可以通过修改介电函数的虚部来实现这个修正。无论用哪种方式,有一条必须反复检查:光学手册上的复折射率写成n+ik,k为正代表介质有吸收;换算成介电常数时,金在可见光区的实部是负的,虚部是正的。如果符号搞反,模型里的金就会从吸收体变成增益介质,仿真结果立刻变成荒谬的负吸收或发散温度,这是新手最容易踩的坑之一。

2.2 频域电磁场方程与吸收功率密度的关系

电磁场部分,我们在COMSOL中使用“电磁波,频域”接口求解频域麦克斯韦方程,本质是在给定入射波长下,把电场分布解出来。对于线性介质、连续波入射,不需要做时域模拟,频域求解效率高得多。金颗粒内部的场分布完全由复介电常数、颗粒几何和周围介质的折射率决定。颗粒内外的电场计算出之后,吸收功率密度可以直接写成时间平均的焦耳热形式:

Q = 0.5 * Re(J · E*)

如果使用介电常数形式,还可以写成与介电常数虚部直接相关的表达式,它本质上是电磁场在损耗介质中的功率耗散。在COMSOL的“电磁热”多物理场耦合中,这个耗散密度会被自动识别并作为热源加到传热方程里,不需要手动写积分。这也是用COMSOL做这类问题的最大优势:电磁场和热源之间的耦合是内置的,你只要把物理场接口选对、耦合节点启用,后面热分析自然会看到正确的热源分布。

2.3 传热方程和纳米尺度边界条件的细节

温度场部分,我们求解的是稳态热传导方程:

∇·(-k∇T) = Q

其中k是导热系数,Q就是上一步算出的电磁功率损耗密度。对水和金颗粒这个体系,金的热导率很高但体积很小,颗粒内部温度梯度几乎可以忽略,主要的温度梯度集中在周围水介质里。因此,决定稳态温升的主要物理量其实是水的导热系数和颗粒半径。

这里有一个经常被问到的问题:颗粒和水之间是否存在界面热阻?在超快时间分辨实验中,金-水界面的热导确实会影响热量从颗粒进入液体的速率,典型数值在几十到几百MW/(m²·K)量级。但对连续激光照射的稳态模型,只要吸收功率不是极高,界面热阻导致的温度阶跃通常很小。一个数量级估算:吸收功率0.1微瓦、颗粒半径40nm时,通过颗粒表面的热流密度大约只有几万W/m²,除以界面热导后温降不足1K。因此在基准算例里忽略界面热阻是合理的,算完再根据实际功率量级判断是否需要补充。

2.4 为什么单向耦合就够,不用做循环迭代

电磁场和温度场之间理论上存在双向耦合:温度升高会改变金的介电常数,而介电常数改变又会影响电场分布和吸收功率。但在单颗粒、连续光、温升只有几K到几十K的范围内,金的介电常数随温度的变化非常小,对吸收截面的影响微乎其微。所以完全没必要做电磁场-温度场的全耦合迭代,单向耦合就够了:先解电磁场,把热源确定,再解传热。这也是COMSOL里“电磁热”多物理场节点的默认逻辑,它会把电磁损耗作为固定热源输入,不进行反向反馈。只有在极高光强、颗粒接近相变温度或者研究热光非线性效应时,才需要升级成全耦合。

3. COMSOL 实操:从几何到求解一步步搭起来

3.1 建模方案:首选二维轴对称,别一上来就搞三维

单个球形金颗粒在均匀水环境中的光热问题,具有严格的旋转对称性,只要入射平面波的波矢沿着旋转对称轴,整个物理系统绕轴对称轴完全一致。这时候用二维轴对称模型是数学上严格等价的,但计算量比三维满模型低一到两个数量级。我建议所有第一次做这个仿真的同学,都不要直接建三维球,先用二维轴对称把流程跑通。三维模型留给研究颗粒朝向、阵列排列、非对称入射等问题时再用。

二维轴对称模型的几何其实很简单:一个代表水的矩形域,中间嵌一个半圆形金颗粒,最外层包一圈PML区。矩形域的高度和宽度取1.2微米左右(对40nm颗粒已经是颗粒尺寸的30倍),PML厚度取0.5个入射波长附近。模型外边界用PML吸收层模拟无限大空间。这样设置能把边界反射对颗粒吸收截面的干扰压到极低。

3.2 几何、单位与参数化设计

打开COMSOL后,建议把全局参数一次性定义好,之后所有设置都引用参数名,这样后面做波长扫描、粒径扫描时不需要改几何或重新建模。基准算例的参数表我放在文末,这里先说明几个关键几何参数:颗粒半径r0设为40nm,水域半径R_w设定为1.2μm,PML外部厚度d_pml设为300nm。几何单位建议直接使用nm,COMSOL内部会自动换算成国际单位制,但在写自定义表达式和积分时要注意,吸收功率密度是W/m³,体积分结果才是瓦特。

对二维轴对称模型,几何绘制时用一个半径为r0的半圆弧表示金球剖面,放在原点附近;用水域矩形包裹;最外层再画一个PML矩形层。旋转轴的设置让COMSOL自动处理,不需要额外指定。

3.3 物理场选择和材料参数的正确导入

物理场选择上,电磁部分选“电磁波,频域”(ewfd),热部分选“固体传热”(ht)。这里容易有个误解:水虽然是液体,但在纳米尺度静水中我们只关心热传导,而且纳米颗粒周围不存在能形成自然对流的尺度和温差条件,所以完全可以用固体传热接口模拟水的导热,没必要引入流体流动的共轭传热。

材料参数是重头戏。水的参数很简单:相对介电常数取1.33的平方(相当于折射率1.33),导热系数0.6W/(m·K),恒压热容4179J/(kg·K),密度1000kg/m³。金的导热系数取318W/(m·K),密度19300kg/m³,恒压热容129J/(kg·K)。更重要的是金的介电常数:通过插值表读入实验光学常数,波长范围覆盖你要扫描的400~700nm。金的数据文件是一个二维表,第一列是波长,第二列是相对介电常数实部,第三列是虚部。注意有些材料库默认给的是折射率n、k,转换的时候保留复数的虚部为正。这个文件在COMSOL里通过“外部材料-插值”创建,曲线平滑性不要求,但插值范围必须覆盖扫描波长区间,靠外推的数据不可信。

3.4 边界条件:PML、背景场和恒温边界

电磁边界条件是这个模型的核心。在二维轴对称模型里,外边界选择PML,并在PML设置中定义其厚度和内部层数。不要把PML域直接加在电磁波接口的几何域内部就行,需要在大几何选择里指定PML的物理位置。

入射光的实现方式,我推荐用散射场公式:总电场=背景场+散射场。背景场设定为沿旋转轴方向传播的平面波,电场方向垂直于旋转轴。这样做的好处一是计算域边界上散射场可以自然衰减,二是提取颗粒的散射场信息很方便。要注意电场的幅值直接决定热源的绝对强度。如果后续要对比某个激光功率的实验,就把E0设置成与该强度对应的值;如果只是要谱线形状,可以先设成1V/m,等找到共振波长后再按强度缩放。

热边界条件就简单很多:水域最外层设为恒定温度293.15K,模拟无穷远处环境温度。颗粒表面不需要特殊边界,默认连续传热。PML层在传热里不应该参与,可以把PML域排除在传热接口之外,或者在PML边界也设成恒温,反正它对颗粒附近温度场的影响可忽略。

3.5 网格划分:既要收敛,也别把内存撑爆

网格是光热仿真里最容易翻车的环节。电磁场必须在颗粒附近精细求解,因为共振模式下颗粒内部的场增强和快速衰减都发生在几十纳米尺度内。我的推荐网格策略分三层:金颗粒内部自由三角形最大单元不超过10nm,颗粒外部近水区域(从颗粒表面向外0.3微米)最大单元不超过20~30nm,远场水区域最大单元不超过波长/10即可。在金颗粒表面再局部加密并添加边界层网格,边界层可以更准确地捕捉近场热流。

网格无关性验证是必须做的。把颗粒内部最大单元尺寸从10nm改成5nm重新计算一次吸收效率,如果两次结果差别超过1%,说明网格还需要加密。如果差别在1%以内,10nm的网格方案就收敛了。在二维轴对称模型里,这套网格规模通常只有几万自由度,求解非常快,你可以放心加密。如果是三维模型,同样的网格策略自由度会轻易到几百万,内存低一点的机器直接卡死,这也是我反复劝大家先用二维轴对称的原因。

3.6 求解器设置与波长扫描策略

求解顺序上,建议先做电磁场波长扫描,确定共振峰位置,再在共振波长下做传热稳态求解。这样有两个好处:一是传热不需要扫描所有波长,省时间;二是只用共振波长的吸收功率作为热源,出来的温度场信号最强。COMSOL中可以建两个研究:研究1是频域扫描,研究2是稳态传热。在研究2里,热源来自研究1的解,通过“电磁热”多物理场自动耦合。如果只想要一个波长,可以不用扫描,直接固定λ0。

波长扫描时,先粗扫(比如450到650nm,步长10nm)看曲线轮廓,找到峰附近后再细化到2nm甚至1nm扫一圈。金颗粒在可见区的吸收峰比较尖锐,粗扫很容易漏掉峰位,导致后续温度计算用了一个偏离共振的波长,结果明显偏低,这个操作顺序一定要养成习惯。求解器用默认的频域直接求解器即可,线性问题不需要迭代,PARDISO在自由度较大时比较稳定。

4. 后处理与指标提取:算完要怎么解读

4.1 吸收效率和吸收截面的标准提取方法

电磁场解完之后,最关心的是颗粒到底吸收了多少功率。在“派生值-体积积分”里,选择金颗粒域,对电磁功率损耗密度做体积分,得到P_abs,单位是W。然后吸收截面σ_abs = P_abs / I0,其中I0是入射光强,等于0.5nε0cE0²。吸收效率Q_abs是吸收截面除以几何截面πr0²。这个Q_abs是等离激元纳米颗粒最重要的无量纲指标,相当于“有效吸收面积是其物理尺寸的多少倍”。对一颗直径80nm的金球在水里,共振波长附近的Q_abs通常在3~6之间,也就是说它吸收光的有效面积比几何截面大好几倍,这正是等离激元光热转换功率高的原因。

4.2 共振峰位的判断与尺寸效应

扫描完成后,绘制Q_abs随波长的变化曲线。对金纳米颗粒在纯水中的情况,峰位一般在520~560nm附近,具体位置取决于颗粒尺寸。粒径增加时,共振峰会发生红移,同时吸收和散射的相对占比也会改变:小颗粒以吸收为主,大颗粒散射比例上升。这个规律在做仿真时可以用来检查模型的合理性:如果你算的40nm金球吸收峰出现在480nm或600nm,多半是材料参数、背景折射率或者网格出了问题,而不是物理上真的有这么奇怪的峰。

4.3 温度场和热影响范围的解读

固定共振波长求解稳态温度后,可以绘制温度云图。典型形态是颗粒中心温度最高,向外迅速衰减,到几百纳米外基本回到环境温度。因为金颗粒热导率高,颗粒内部温度几乎均匀;温度梯度的主要部分落在水介质中。沿着径向画一条温度剖面线,你可以看到一个大致按1/r衰减的趋势,这种趋势可以用球壳导热的解析公式验证:

ΔT(r) = P_abs / (4πk_w·r)

其中k_w是水的导热系数,r是离球心的距离。这个公式在远场区非常准,可以用它来检验COMSOL的稳态温度结果是否合理。比如40nm颗粒,吸收功率0.2微瓦时,颗粒表面r=40nm处温升理论值大概是0.6K量级,如果仿真结果比这个数量级差出十倍以上,说明热源设置或者边界条件有问题。

4.4 光热转换效率的工程意义

对光热应用来说,最后往往要看的是光热转换效率,定义是吸收效率占消光效率的比例,也就是吸收截面/(吸收截面+散射截面)。10~50nm的金颗粒散射占比小,光热效率可以到90%以上;颗粒长大到100nm以上,散射比例显著增加,光热效率下降。这个指标在COMSOL里也很容易算:散射截面可以用远场探针算出,或者用消光截面减去吸收截面。做光热治疗剂量设计时,真正起作用的就是吸收截面,而不是消光截面,这个区别要特别注意。

5. 实际跑仿真中的坑和排查实录

5.1 吸收峰位置不对,先检查网格和数据来源

我最常见的情况是学生在群里发来一条吸收谱线,峰位比文献值偏蓝了二三十纳米。排查顺序固定三件事:材料介电常数是不是用的金而不是银等其它金属;水的折射率是否设置成了1(真空)而不是1.33;颗粒表面的网格是不是太粗。网格太粗时,颗粒内部的局域场分辨率不够,会人为压低吸收效率、压低红移,让峰位看起来蓝移。如果这三项都没问题,再检查PML是否有效吸收,边界反射会形成驻波,导致吸收谱出现异常波纹。

5.2 温度全为零或者温度低得离谱,十有八九是热源没接上

传热解出来一点温升都没有,最直接的原因是“电磁热”多物理场耦合没有启用。很多新手先单独建了电磁模型,又单独建了传热模型,物理场之间没有通过多物理场节点连接,热源Q在传热方程里就是0。另一个更隐蔽的原因是背景场用了E0=1V/m做归一化扫描,找到共振波长后没有把E0改成实际光强对应的值就直接算温度。因为热源正比于E0的平方,E0从1V/m改成75000V/m,温度要差几十亿倍,这种错误会让结果看起来像“颗粒根本没发热”。我建议从一开始就给E0设为实际光强对应的场强值,例如想模拟1kW/cm²(1×10⁷W/m²)的入射强度,E0大约为75000V/m,这样热源和温度都是真实物理量,后处理不需要任何缩放。

5.3 模型不收敛或直接被内存击穿

二维轴对称模型几乎不会爆内存,但如果有人在三维模型里同时用了加密网格、PML和宽波长扫描,自由度数上千万是常有的事,消费级电脑直接卡到不能操作。这种情况下我的建议是先砍需求再逐步加复杂度:先做二维轴对称,把物理搞明白;确定必须三维时,至少把扫描分成几段跑,每段跑完再合并数据;或者用粗网格先找峰,细网格只算峰值附近的几个波长。求解器上,三维模型优先用PARDISO,不要用默认迭代器硬扛。

5.4 能量不守恒检查:校验收敛的最硬标准

判断一个电磁仿真做得好不好,最好用的办法是能量守恒核算。在金颗粒内部积分得到吸收功率P_abs,再通过远场积分算出散射功率P_sca,两者相加应该等于消光截面乘以入射光强。如果算出来的消光效率和文献里同尺寸颗粒的对不上,差出10%以上,基本可以断定网格、材料参数或者PML有问题。这个守恒校验比单纯盯着温度曲线可靠得多,因为它暴露的是物理层面的错误,不是后处理显示层面的错误。

5.5 二维轴对称和三维结果不一致,通常是背景场设置有问题

有的同学按教程跑了二维轴对称,又拿三维模型验证,发现吸收峰对不上。最常见的原因是三维模型里平面波入射方向设置成了沿某个轴,但电场极化方向没有明确指定,导致等效背景场在颗粒处产生不正确的相位分布。球颗粒散射理论本身与入射波的极化方向无关,但如果背景场在计算域边界上没有被正确描述,散射场和总场的叠加就会出错。解决办法是把三维背景场的波矢方向、极化方向都显式写清楚,并确保计算域边界距离颗粒足够远(一般建议至少一个波长)。

6. 可复现的基准算例与参数速查

6.1 完整参数清单

下面这张表可以直接作为COMSOL全局参数输入,跑通后再按自己需求改。

参数值说明
r040 nm金颗粒半径
R_w1.2 μm水介质域半径(不含PML)
d_pml300 nmPML厚度
lambda0530 nm入射波长(可扫描)
n_w1.33水的折射率
E075000 V/m对应入射强度约1kW/cm²
T0293.15 K环境温度
k_w0.6 W/(m·K)水的导热系数
k_gold318 W/(m·K)金的导热系数

入射波沿旋转对称轴方向,背景场电场方向垂直于旋转轴,这是二维轴对称模型的配套设置。金介电常数插值表覆盖400~700nm,数据来源采用实验光学常数,虚部保留正值。

6.2 快速复现的建模步骤速查

步骤很简单:第一步新建二维轴对称模型,加入“电磁波,频域”和“固体传热”两个物理场;第二步用参数画几何,半圆颗粒、水域矩形、PML层;第三步材料设置,水区用自定义折射率和热物性,金区用插值介电常数和热物性;第四步电磁波接口开启背景场,设置E0和波矢方向,外边界选择PML;第五步传热接口设置恒温边界,把PML域排除;第六步添加“电磁热”多物理场耦合节点;第七步网格配置按前面策略;第八步研究1做波长扫描,确定峰位;第九步固定峰位做稳态传热;第十步后处理积分吸收功率、绘制温升曲线。整个过程在一台普通办公电脑上几分钟内就能跑完。

6.3 COMSOL with MATLAB 脚本化修改关键参数

如果你有多个尺寸或多个材料要批量算,完全没必要每次在GUI里手动改参数。COMSOL with MATLAB 或者App开发器可以做到参数化复现。下面是一个轮廓脚本,核心是用模型对象设置参数、运行研究、导出结果:

model = mphopen('AuNP_photoheat.mph'); model.param.set('r0', '40[nm]'); model.param.set('lambda0', '530[nm]'); model.param.set('E0', '75000[V/m]'); model.study('std1').run(); % 读取吸收功率(体积分) Pabs = model.result.numerical('intop_pabs').getReal(); fprintf('Absorbed power: %.3e W\n', Pabs); % 修改粒径并重新计算 model.param.set('r0', '60[nm]'); model.study('std1').run();

这段代码只是个骨架,实际使用时需要根据你模型里的组件标签和派生值的名称做微调。核心思路是把几何尺寸、入射波长、入射强度这些参数全部集中到全局参数里,脚本只负责改参数-求解-提数据,能够把几十个粒径的扫描工作从手动点改中解放出来。

这个内容后续想继续扩展的话,可以把单颗粒改成二聚体甚至周期性阵列,研究颗粒间距对热场的影响;也可以在颗粒表面加一层二氧化硅壳或生物分子层,模拟实际功能化颗粒的光热响应。但所有这些扩展都建立在今天这个单颗粒基准算例跑得准确、收敛、可重复的基础上。我个人踩过几次坑之后的体会是:做这类多物理场仿真,理论公式是必要条件,COMSOL操作是充分条件,但真正决定结果可不可信的,永远是“解析校验——吸收守恒——网格无关性”这三层验收关卡,别急着改参数,先把基准算例摸透,后面才会越跑越顺。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询