企业官网建设流程全解析

第一次用COMSOL算石墨烯太赫兹透射时，我干过一件很天真的蠢事：在几何里画一个0.34nm厚的矩形块当石墨烯，然后让网格器自己去划分。结果模型直接卡死，网格质量掉到0.001量级，求解器一步都迈不出去。后来我才彻底想明白，石墨烯这种单原子层二维材料，在COMSOL里根本不该用三维体单元去建，而是压缩成一个边界，用表面电导率去描述它的电磁响应。表面电导率怎么选，又牵出德鲁得和Kubo两个经典模型。

这篇文章把我从“画一个薄块”到“能在同一个模型里切换太赫兹德鲁得模型和近红外Kubo模型”的完整思路、公式落地、边界条件挂接和调试经验整理出来。如果你正在做石墨烯超表面、太赫兹调制器、近红外光电探测器或等离激元器件仿真，这里面的内容应该能帮你少踩不少坑。我用的COMSOL版本是6.x，5.x的界面差异不大，也能照方抓药。

1. 为什么石墨烯不能按三维体材料建：尺度灾难与二维边界方案

1.1 0.34nm与波长之间隔着的六个数量级

石墨烯单层厚度0.34nm，这个数字看着不起眼，放到电磁仿真里就是一场灾难。太赫兹频段对应波长约30μm到3mm，近红外频段对应波长约0.7μm到1.6μm。你算一下，0.34nm和700nm这段最“短”的近红外波长之间也差了差不多三个数量级，和太赫兹波长相比差距更大，接近六个数量级。

如果把石墨烯当成三维薄板来建模，网格就必须同时解析0.34nm的厚度和微米级波长。这意味着最小网格与最大网格尺寸之间要跨几个数量级，生成的单元几乎全是超高纵横比的畸形单元，网格质量经常掉到0.01以下。MUMPS求解器对这种网格会极度敏感，要么算不动，要么算完也毫无可信度。

有人可能会想，那我把石墨烯的厚度人为加厚到1nm甚至5nm，然后在材料参数里相应调低体电导率，不就行了吗？这种等效思路在静态场里勉强能用，但在光频和太赫兹频段的电磁波仿真里并不可靠。因为你改变厚度后，层内阻抗、表面电流分布和相位积累都会变化，尤其涉及石墨烯等离激元或透射相位的时候，结果会明显偏离真实物理。

记住一句话：石墨烯的电磁响应本质上是二维的，COMSOL里最自然的处理方式是把它当作一个零厚度边界，而不是一个三维体材料。这样网格只需要解析波长尺度，不需要去管0.34nm这种原子级厚度。

1.2 两个边界入口：表面电流密度与过渡边界条件

把石墨烯压成边界之后，物理上真正进入电磁场方程的是它的表面电导率。对一个位于介质分界面上的导电薄膜，切向磁场跨越边界时会产生跳变，跳变的大小就是面板电流密度：

J_s = σ_s E_t

这里σ_s就是石墨烯的表面电导率，E_t是边界切向电场，单位是V/m，σ_s单位是西门子S，J_s单位是A/m。在COMSOL的“电磁波，频域”接口中，实现这个关系通常有两个入口。

第一个是“表面电流密度”边界条件。它直接在选定的边界上施加一个面电流密度表达式，最贴近上述方程，也是我最常用的一种方式。比如J_sx = sigmaEx，J_sy = sigmaEy，只要电场切向分量的提取方式正确，结果就很干净。

第二个是“过渡边界条件”TBC。它把薄层当成一个传递阻抗来处理，需要设定厚度、相对介电常数、电导率等参数。TBC的优势是能描述层内非均匀场和垂直方向响应，但对石墨烯这种厚度远小于趋肤深度的材料，TBC和表面电流密度在大多数频段下结果几乎一致。真正需要二选一的场景，我会在后面第5章里展开。

网格策略上，石墨烯边界本身需要细划，尤其是把它当作超表面或天线结构一部分时，边界网格尺寸要能分辨特征金属结构。而周围的空气或衬底网格按波长控制即可，太赫兹频段所有结构电尺寸都很小，网格压力小，近红外频段则需要更细。

2. 太赫兹频段的德鲁得模型：从物理推导到变量落地

2.1 为什么太赫兹频段是德鲁得的主场

太赫兹的光子能量非常小。1THz对应的光子能量只有约4.1meV，而石墨烯的费米能级通常通过掺杂或栅极电压调到0.2eV到1eV之间。也就是说，太赫兹频段光子能量远小于费米能级，电子在能带内做带内跃迁是绝对主导的物理过程，带间跃迁因为泡利阻塞被压制掉。

这种情况下，石墨烯的表面电导率可以用经典的德鲁得形式描述。常用的强掺杂近似式是：

σ_Drude(ω) = (e²E_F)/(πħ²) × i/(ω + i/τ)

其中e是电子电荷，E_F是相对狄拉克点的费米能级，ħ是约化普朗克常数，τ是动量弛豫时间。这个公式在零频极限下回到直流电导率σ_DC = e²E_Fτ/(πħ²)，物理含义很直观：载流子越多、弛豫时间越长，导电性越好。

这里必须提一个容易踩的坑：i/(ω + i/τ)这个写法里的虚部符号，与你使用的时谐因子约定强相关。如果你从某篇使用e^{-iωt}约定的论文里抄公式，填进COMSOL时很可能发现虚部符号和预期相反。COMSOL本身的时谐约定在不同物理场接口下也有说明，我建议你不要硬背“COMSOL是正还是负”，而是做一个解析验证，这个习惯能省很多脑细胞。

2.2 COMSOL里怎样定义德鲁得表面电导率

以COMSOL 6.x为例。先在“全局定义→参数”里把材料相关量都放好：

然后在“全局定义→变量”里定义角频率和电导率：

omega = 2*pi*freq[1/s]
sigma_drude = e_const^2*EF_energy/(pi*hbar_const^2)*i/(omega+i/tau)

注意变量单位要设为S，COMSOL的单位检查如果报警告，可以在变量表达式末尾手动补上[S]，或者忽略警告直接使用，但最好保持单位正确，后面做结果后处理时才不会乱。

物理场部分，选中代表石墨烯的那个边界，添加“表面电流密度”边界条件。如果石墨烯平面正好落在xy平面内，可以直接写：

J_sx = sigma_drude * ewfd.Ex
J_sy = sigma_drude * ewfd.Ey

如果石墨烯边界是斜的或曲面，就不能直接用全局的Ex、Ey，而要把电场投影到边界局部坐标系的t1、t2方向上，这个细节会让很多人卡住。

2.3 太赫兹扫描设置与参数敏感性

研究部分选“频域”，把频率范围设成0.1THz到10THz，步长0.1THz或0.2THz。太赫兹频段下，如果几何尺寸在微米级，整个模型比波长小很多，MUMPS直接求解器最稳，基本不会出现迭代不收敛的问题。

参数敏感性方面，τ是第一个要调的。高迁移率的悬浮石墨烯τ可以达到皮秒量级，CVD石墨烯因为缺陷和声子散射，τ通常在几十飞秒到几百飞秒。同一个模型，把τ从50fs改成500fs，太赫兹段的表面电导率实部变化可能不太离谱，但虚部会变很多，进而影响反射相位和吸收峰宽度。

另外要注意，低频端1/τ可能和ω同量级甚至更大，这时德鲁得模型的实部几乎与频率无关，虚部随频率升高而下降，这是典型的自由载流子响应。如果你仿出来透射率在很宽频带内是平的、没有吸收峰，往往不是模型错了，而是自由载流子本来就没有共振，你要找的共振来自结构而不是材料本身。

3. 近红外频段的Kubo模型：带内与带间的完整展开

3.1 近红外为什么必须换Kubo模型

到了近红外频段，光子能量变成0.7eV到1.7eV，和石墨烯费米能级已经可比。一旦光子能量超过两倍费米能级，也就是ħω > 2E_F，电子从价带直接激发到导带的带间吸收通道就被打开，这部分吸收比带内自由载流子吸收大得多。

如果近红外仿真继续用德鲁得模型，等于丢掉了带间项，算出来的吸收率会明显偏低，尤其在石墨烯费米能级较低的时候，偏差会非常夸张。这就是Kubo模型登场的地方，它在德鲁得的基础上，把带内和带间两项都包含进去，是一个既适合低频也适合高频的完整表面电导率表达式。

3.2 Kubo公式的有限温度形式与COMSOL实现

我实际使用的Kubo模型通常分成两项。带内项的形式是：

σ_intra(ω) = (2e²k_BT)/(πħ²) × i/(ω + i/τ) × ln[2cosh(E_F/(2k_BT))]

带间项的形式是：

σ_inter(ω) = e²/(4ħ) × { 1/2 + 1/π·atan[(ħω - 2E_F)/(2k_BT)] - i/(2π)·ln[((ħω + 2E_F)² + (2k_BT)²) / ((ħω - 2E_F)² + (2k_BT)²)] }

总电导率就是两项相加：σ_Kubo = σ_intra + σ_inter。

上面的式子看起来长，但在COMSOL里只是几个内置函数的嵌套。在“全局定义→变量”里继续写：

kT = k_B_const*T0
sigma_intra = 2*e_const^2*kT/(pi*hbar_const^2)*(i/(omega+i/tau))*log(2*cosh(EF_energy/(2*kT)))
sigma_inter = e_const^2/(4*hbar_const)*(0.5+1/pi*atan((hbar_const*omega-2*EF_energy)/(2*kT))-i/(2*pi)*log(((hbar_const*omega+2*EF_energy)^2+(2*kT)^2)/((hbar_const*omega-2*EF_energy)^2+(2*kT)^2)))
sigma_kubo = sigma_intra + sigma_inter

这里有几个细节容易出错。COMSOL里log默认是自然对数，不是常用对数。atan和cosh都是内置函数，直接用即可。所有能量项必须统一单位，EF_energy和kT都要是焦耳，hbar_const*omega也是焦耳，这样括号内无量纲，单位检查才能通过。

3.3 化学势、温度和散射率对近红外结果的影响

Kubo模型里最敏感的参数是EF。低掺杂时，比如0.1eV，近红外光子很容易超过2E_F，带间吸收会很明显，单层石墨烯吸收率可以逼近2.3%。高掺杂时，比如0.6eV，带间吸收边被推到1.2eV以上，低于这条线的光几乎无法引起带间跃迁，透射率反而升高。所以石墨烯太赫兹调制器的原理本质上就是通过栅压改变EF，让带间吸收边扫过目标激光波长。

温度项看似不起眼，实则影响很大。300K对应k_BT约为25.9meV，它决定了带间吸收边的展宽。如果用零温近似公式里的阶跃函数，得到的透射光谱会在ħω=2E_F处突然跳变，和实际实验曲线对不上。COMSOL仿真应当使用有限温度版本，也就是带atan和log的完整形式，这样曲线在吸收边附近有平滑的过渡。

τ的影响在高频段变得微妙。当ħω远大于2E_F时，带间项几乎不依赖τ，所以近红外下如果主要看石墨烯本征吸收，τ的误差影响相对小。但如果做的是金属-石墨烯混合等离激元纳米结构，共振品质因子Q受τ影响很大，τ太小会把共振峰直接抹平。这时候τ就一定要仔细校准。

4. 太赫兹到近红外：德鲁得与Kubo的切换问题

4.1 德鲁得其实是Kubo带内项的特例

很多人会有疑问：既然Kubo模型全频段都能用，为什么不直接丢掉德鲁得模型？原因有两层。第一，德鲁得形式简单，拟合参数少，在太赫兹段做参数反演时更直观。第二，理解两者关系能帮你判断仿真结果的合理性。

当E_F远大于k_BT时，带内项中的ln[2cosh(E_F/(2k_BT))]可以近似为E_F/(k_BT)，代入σ_intra后正好退化成德鲁得形式。所以德鲁得不是另一个独立模型，它是Kubo带内项在强掺杂低温近似下的简写。在太赫兹频段，带间项被2E_F远远隔离，于是德鲁得和Kubo算出来几乎重合。

在实际模型中，我建议把sigma_drude和sigma_kubo都定义在全局变量里，然后用一个开关参数切换。比如设置use_kubo这个参数，表面电流密度边界条件表达式写成：

sigma_graphene = if(use_kubo>0.5, sigma_kubo, sigma_drude)

这样同一个模型，改一个参数就能横扫太赫兹和近红外，不用每次改动边界条件。

4.2 跨度三个数量级的频率扫描设置

太赫兹到近红外的频率跨度非常大，从0.1THz到400THz，超过三个数量级。我不建议在一个频域研究里从低扫到高，因为PML厚度、网格尺度和求解器策略在两端差异太大。分别做两个研究是最省事的。

太赫兹研究用频率扫描，范围0.1THz到10THz；近红外研究可以直接扫波长，从1600nm扫到700nm，波长扫出来更直观。COMSOL研究设置里的频率参数可以写range(0.1[THz],0.2[THz],10[THz])或者range(1600[nm],-50[nm],700[nm])，注意近红外扫波长的单位写法。

还有一个细节：如果研究中同时有参数扫描和频率扫描，COMSOL会先执行参数扫描再执行频率扫描，生成的解序列非常大。建议参数扫描的EF、τ值控制在一二十个以内，否则输出文件会庞大无比。

4.3 网格、PML和求解器的频段差异化处理

太赫兹频段一切结构都比波长小很多，PML厚度不是特别敏感，网格也只需要保证石墨烯边界形状准确，周围空气网格可以放到粗糙档。但近红外频段不同，PML至少要取四分之一波长，也就是几百纳米；模型内部最大网格尺寸不能超过150nm，如果有金属纳米结构或纳米间隙，局部加密到10nm到20nm也不稀奇。

求解器方面，太赫兹频段的模型自由度通常小，MUMPS直接求解器最稳。近红外模型如果比较复杂，自由度很容易超过几百万，考虑换成PARDISO或迭代求解器GMRES。不管用哪种，我都建议先单独算一个频点，确认收敛后再做频率扫描，不然一个不收敛的频率点会在参数扫描列表中拖垮整个研究。

5. 仿真结果对不上文献时的五个排查方向

5.1 表面电导率和体电导率的单位陷阱

这是最常发生的低级错误。文献里石墨烯表面电导率σ_s用的是西门子S，如果你在COMSOL材料节点里直接填入σ_s，实际上体电导率的单位是S/m，量纲差出米分之一，结果必然离谱。反之，如果要用TBC，需要设置等效体电导率，这时要做除法σ_bulk = σ_s / d，d是假设厚度，单位用米。

排查方法很简单：先建一个自由空间中的单层石墨烯薄层模型，不做任何复杂结构。平面波垂直入射，石墨烯面用表面电导率边界条件，理论上功率透射率为：

T = |2/(2 + η0·σ)|²

η0约376.73Ω。把COMSOL算出来的透射率和这个公式对比，偏差超过千分之几就说明边界条件设置有问题。这个测试模型我强烈建议保留在模型库里，以后每个新版本COMSOL都能用来做基准验证。

5.2 散射率τ和衰减率γ的2倍关系

不同文献对Kubo模型里散射项的写法并不统一。有的用弛豫时间τ，有的用散射率Γ，还有的写γ。部分公式里γ等于1/τ，部分则定义γ为1/(2τ)。如果你从两篇文献里各抄一段公式拼在一起，低频极限的直流电导率会差出2倍，太赫兹段实部虚部比例也会对不上。

解决思路是不要依赖记忆，而是从直流极限校准。把omega设成0，检查σ_intra的低频表达式是否还原成你期望的σ_DC数值。比如GO天线? 石墨烯典型σ_DC在EF=0.3eV、τ=100fs时约为1mS量级，如果算出来明显大了或小了，就检查散射率的系数。

5.3 切向电场分量取错导致极化方向异常

在平面且与坐标轴对齐的情况下，把J_sx写成sigma*ewfd.Ex、J_sy写成sigma*ewfd.Ey是没问题的。但一旦石墨烯边界是斜面、曲面或任意走向，这种写法就错了。因为全局Ex、Ey包含了法向分量，而后者的贡献不应该进入表面电流公式。

正确的做法是用边界局部坐标系的切向方向。在COMSOL边界上，电场可以投影到边界的两个切线方向，通常可以用边界坐标系变量或使用投影算子。如果实在找不到对应变量，最简单的办法是手动把几何坐标系旋转到与石墨烯平面一致，再用全局Ex、Ey。这个方法虽然笨，但在大量规则结构里足够稳定。

5.4 Kubo带间项出现阶跃或负吸收时查温度项

用零温Kubo公式时，带间项里的阶跃函数会把σ_inter的实部直接切成0和e²/(4ħ)两段，虚部在边界附近会出现对数奇异性。在COMSOL扫频结果里表现为透射率曲线出现生硬的台阶，甚至吸收率在某些频点为负值，这是典型的数值非物理现象。

解决办法很简单：把T设成实际测量温度，或者在最极端情况下至少给kT设一个10meV到25.9meV的热展宽。有限温度Kubo公式的两个项都是解析光滑的，不会出现除零或负吸收。这个坑我踩过一次，排查了整整两天才发现是公式退化问题，不是模型结构问题。

5.5 高频下TBC与表面电流密度结果不一致

做太赫兹大尺寸结构时，TBC和表面电流密度边界条件差异可以忽略。但到了近红外，尤其涉及石墨烯等离激元共振的纳米结构时，TBC的薄层近似可能引入额外误差，现象是共振峰位置偏移几十纳米，品质因子也偏高。

我个人的经验是：一旦模型特征尺寸进入纳米量级，且你关心共振波长，优先使用表面电流密度边界条件。它直接定义面板电流，和实验中等效电路的理解更一致。TBC更适合关注薄层内部场分布的情况，比如研究双层石墨烯或异质结垂直电场时，它才有不可替代的价值。做超表面和传感器的人，默认先从表面电流密度开始，大概率不会走弯路。

最后分享一点自己的习惯：我总会在模型的全局变量里同时保留德鲁得和Kubo两套表达式，用use_kubo这个开关参数做切换。这样同一个几何结构既可以跑到太赫兹波段验证带内输运，也可以切到近红外观察带间吸收。石墨烯这个材料最迷人的地方就在这里——同样是碳原子单层，在不同的频率窗口里展现出的物理过程完全不同。德鲁得和Kubo背后其实是同一套微观跃迁，只是不同频段下谁占主导的问题。理解了这层关系，你在COMSOL里调参数时就不容易迷失方向。