金纳米颗粒的光热效应这几年在纳米医学、等离激元传感、光热抗菌、光声成像领域被反复提到但真要把“光照一颗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测量值或者用自由电子气模型自行构建。实验数据的好处是准坏处是版本和数据来源差异会导致结果有百分之几的差别不同文献使用的数据点也不一样只要保证峰位偏移在可接受范围内就可以。自由电子气模型的好处是参数少、形式清楚可以方便地加入尺寸修正。对孤立纳米颗粒当颗粒半径接近电子平均自由程时电子与颗粒边界的碰撞会显著增加阻尼等效于在费米速度基础上增加一个散射频率这会导致共振峰展宽、吸收效率降低。在仿真中可以通过修改介电函数的虚部来实现这个修正。无论用哪种方式有一条必须反复检查光学手册上的复折射率写成nikk为正代表介质有吸收换算成介电常数时金在可见光区的实部是负的虚部是正的。如果符号搞反模型里的金就会从吸收体变成增益介质仿真结果立刻变成荒谬的负吸收或发散温度这是新手最容易踩的坑之一。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μmPML外部厚度d_pml设为300nm。几何单位建议直接使用nmCOMSOL内部会自动换算成国际单位制但在写自定义表达式和积分时要注意吸收功率密度是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微瓦时颗粒表面r40nm处温升理论值大概是0.6K量级如果仿真结果比这个数量级差出十倍以上说明热源设置或者边界条件有问题。4.4 光热转换效率的工程意义对光热应用来说最后往往要看的是光热转换效率定义是吸收效率占消光效率的比例也就是吸收截面/(吸收截面散射截面)。10~50nm的金颗粒散射占比小光热效率可以到90%以上颗粒长大到100nm以上散射比例显著增加光热效率下降。这个指标在COMSOL里也很容易算散射截面可以用远场探针算出或者用消光截面减去吸收截面。做光热治疗剂量设计时真正起作用的就是吸收截面而不是消光截面这个区别要特别注意。5. 实际跑仿真中的坑和排查实录5.1 吸收峰位置不对先检查网格和数据来源我最常见的情况是学生在群里发来一条吸收谱线峰位比文献值偏蓝了二三十纳米。排查顺序固定三件事材料介电常数是不是用的金而不是银等其它金属水的折射率是否设置成了1真空而不是1.33颗粒表面的网格是不是太粗。网格太粗时颗粒内部的局域场分辨率不够会人为压低吸收效率、压低红移让峰位看起来蓝移。如果这三项都没问题再检查PML是否有效吸收边界反射会形成驻波导致吸收谱出现异常波纹。5.2 温度全为零或者温度低得离谱十有八九是热源没接上传热解出来一点温升都没有最直接的原因是“电磁热”多物理场耦合没有启用。很多新手先单独建了电磁模型又单独建了传热模型物理场之间没有通过多物理场节点连接热源Q在传热方程里就是0。另一个更隐蔽的原因是背景场用了E01V/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水介质域半径不含PMLd_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操作是充分条件但真正决定结果可不可信的永远是“解析校验——吸收守恒——网格无关性”这三层验收关卡别急着改参数先把基准算例摸透后面才会越跑越顺。
金纳米颗粒的光热效应这几年在纳米医学、等离激元传感、光热抗菌、光声成像领域被反复提到但真要把“光照一颗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测量值或者用自由电子气模型自行构建。实验数据的好处是准坏处是版本和数据来源差异会导致结果有百分之几的差别不同文献使用的数据点也不一样只要保证峰位偏移在可接受范围内就可以。自由电子气模型的好处是参数少、形式清楚可以方便地加入尺寸修正。对孤立纳米颗粒当颗粒半径接近电子平均自由程时电子与颗粒边界的碰撞会显著增加阻尼等效于在费米速度基础上增加一个散射频率这会导致共振峰展宽、吸收效率降低。在仿真中可以通过修改介电函数的虚部来实现这个修正。无论用哪种方式有一条必须反复检查光学手册上的复折射率写成nikk为正代表介质有吸收换算成介电常数时金在可见光区的实部是负的虚部是正的。如果符号搞反模型里的金就会从吸收体变成增益介质仿真结果立刻变成荒谬的负吸收或发散温度这是新手最容易踩的坑之一。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μmPML外部厚度d_pml设为300nm。几何单位建议直接使用nmCOMSOL内部会自动换算成国际单位制但在写自定义表达式和积分时要注意吸收功率密度是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微瓦时颗粒表面r40nm处温升理论值大概是0.6K量级如果仿真结果比这个数量级差出十倍以上说明热源设置或者边界条件有问题。4.4 光热转换效率的工程意义对光热应用来说最后往往要看的是光热转换效率定义是吸收效率占消光效率的比例也就是吸收截面/(吸收截面散射截面)。10~50nm的金颗粒散射占比小光热效率可以到90%以上颗粒长大到100nm以上散射比例显著增加光热效率下降。这个指标在COMSOL里也很容易算散射截面可以用远场探针算出或者用消光截面减去吸收截面。做光热治疗剂量设计时真正起作用的就是吸收截面而不是消光截面这个区别要特别注意。5. 实际跑仿真中的坑和排查实录5.1 吸收峰位置不对先检查网格和数据来源我最常见的情况是学生在群里发来一条吸收谱线峰位比文献值偏蓝了二三十纳米。排查顺序固定三件事材料介电常数是不是用的金而不是银等其它金属水的折射率是否设置成了1真空而不是1.33颗粒表面的网格是不是太粗。网格太粗时颗粒内部的局域场分辨率不够会人为压低吸收效率、压低红移让峰位看起来蓝移。如果这三项都没问题再检查PML是否有效吸收边界反射会形成驻波导致吸收谱出现异常波纹。5.2 温度全为零或者温度低得离谱十有八九是热源没接上传热解出来一点温升都没有最直接的原因是“电磁热”多物理场耦合没有启用。很多新手先单独建了电磁模型又单独建了传热模型物理场之间没有通过多物理场节点连接热源Q在传热方程里就是0。另一个更隐蔽的原因是背景场用了E01V/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水介质域半径不含PMLd_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操作是充分条件但真正决定结果可不可信的永远是“解析校验——吸收守恒——网格无关性”这三层验收关卡别急着改参数先把基准算例摸透后面才会越跑越顺。
