☰
基于Comsol的相场法多晶介电击穿电树枝动态模拟
2026/10/7 22:15:36 网站建设 项目流程

搞介电材料和绝缘诊断的人,估计都见过那张经典的“电树枝”照片:一束银白色的枝杈从针电极尖端朝介质深处疯狂生长,像闪电劈进冰面留下的裂纹。但真正轮到自己上手去模拟这玩意儿,才发现实验切片只能看静态结果,想要复现动态生长路径,尤其是多晶结构里的那种“贴着晶界走”的诡异倾向,靠经验公式根本推不出来。我在这个坑里折腾了大半年,最后选定用Comsol做多晶介电击穿的相场模拟,总算把电树枝的路径、分叉和晶界偏转都跑出来了。这篇就从头讲一遍:为什么相场能行、多晶结构怎么搭、方程怎么写进Comsol,以及那些文档里从来不会写的坑。

1. 为什么要用相场方法去算电树枝

1.1 传统模型的局限

先说清楚背景。电树枝的本质是介质内部局部电场超过材料承受极限后,发生不可逆的损伤并逐步扩展成树枝状通道。早期模拟常用蒙特卡洛随机网络模型——把区域离散成格子,按局部场强超过阈值的概率随机击穿格点,一步步生长。这类模型能长出漂亮的树形图案,统计分形维数也能和实验对上一些,但有一个致命伤:格子概率是人为拍出来的,路径随机性大,很难系统地研究材料微观结构的影响,比如晶界、气孔、第二相颗粒。你就算跑一万次,也回答不了“为什么这批试样总是在晶界处分叉”这种问题。

后来也见过用有限元加生死单元的方式模拟,电场每步迭代、超过阈值就“杀死”单元。逻辑上很直观,但问题在于暴露了界面依赖网格的毛病:树枝走哪条路、晶体里的取向如何影响扩展,都要先预设裂纹路径,而且每次单元变形后电场重新求解,容易振荡。说白了,它还是在“猜路径”,而不是让物理自行决定路径。

1.2 相场法的核心优势

相场法的思路完全不同,它不再去追踪一条明确的尖裂纹或导电通道边界,而是用一个连续变量s来区分材料和损伤状态,s=0代表绝缘完好,s=1代表已经完全击穿导电(或者退化成缺陷通道)。这个s在空间上有一个梯度过渡区域,叫扩散界面,界面本身会随着自由能驱动自发演化,不需要任何先验路径假设。

通俗点理解:传统模型就像用笔在地图上画一条登山路线,每走一步人工判断哪里好走;相场模型则像在水面上撒一层染色剂,不同区域的能量差会让染色剂自然扩散,哪里低势垒就往哪里流,路径完全由势能地形自己决定。电场集中、晶界弱化、材料各向异性都能直接体现在自由能泛函的构成里,这样一个框架就能统一处理多晶复杂微结构和电树枝的动态扩展。

1.3 多晶场景为什么非它不可

多晶材料的结构特点是:晶粒内部有序、晶界处原子排布混乱,还常常富集缺陷和杂质。实验中观察到的电树枝大量沿晶界扩展,原因就是晶界区域的介电强度低、电导率高,在电场下更容易热失稳。要把这种差异写进模拟,传统随机模型只能用“晶界格子击穿概率高”这种后验修正,而相场模型可以在自由能项里直接把晶界定义为一个独立的物理区域,赋予不同的势垒高度、界面能和介电响应。换句话说,晶界不是一个需要人为指定的“路”,而是材料内部的一个真实能量景观,电树枝走过去只是因为那边的能量门槛更低。

我当时看到文献里有人用相场模拟铁电薄膜的击穿路径,晶界对树枝的导向作用体现得非常清晰,才意识到这套方法才是研究多晶结构的正路。

2. 多晶结构与相场方程的整体设计

2.1 自由能泛函怎么搭

相场模型的起点是写下一个总自由能泛函,通常是三块:

F = ∫Ω [ W(s) + (κ/2)|∇s|² + f_e(s, E) ] dΩ

第一项是双阱势,形状像M型的双坑曲线,最典型的写法是:

W(s) = 4γ s²(1-s)²

γ决定了s=0和s=1两个稳态之间的势垒高度。物理上,这个势垒对应击穿所需越过的能量门槛。第二项是梯度界面能,κ控制s从0过渡到1的界面宽度和能量,κ越大界面越宽,树枝越“钝”,分叉倾向也会变化。第三项是电场的能量贡献,f_e耦合电场变量与相场变量。

电场能具体形式取决于你采用哪种击穿假设。我做的是介电主导的退化模型,写法是:

f_e = -1/2 ε0 εr(s) |E|²

其中介电响应随s变化,最简单的线性插值是:

εr(s) = εr完好 + (εr击穿 - εr完好) · s

回到最开始的问题:击穿区域到底该是介电常数升高还是降低?文献里两种都有。如果模拟的是“局部热致退化”,通常让εr下降,导致该处电场应力集中加剧,形成正反馈自加速。如果模拟的是“已击穿形成导电通道”,那么更合理的做法是引入电导率项而不是单纯改变介电常数。我在实际建模时推荐第二种:在击穿相引入高电导率σbr,在完好相保留极低绝缘电导率σins,这样击穿通道在电位分布上更接近等势体,树枝尖端的电场集中也更真实。

2.2 演化方程与电场方程

相场变量的松弛演化用Allen-Cahn方程描述:

∂s/∂t = -M ( δF/δs )

展开来看,δF/δs = W'(s) - κ∇²s + ∂f_e/∂s,其中W'(s) = 8γs(1-s)(1-2s)。M是迁移率,相当于树枝扩展速度的动力学系数,它的数值标定直接影响模拟时间与实际时间的对应关系。

电场方程我是在Comsol的“电流(ec)”物理场里实现的,而不是纯静电场:

∇·[ σ(s) ∇V ] = 0

σ(s) = σins + (σbr - σins) · s

当s一点点变成1时,局部电导率从绝缘体量级跳到导体量级,电场重新分布。未击穿区域的电场被拉高,已击穿通道内部电场被拉到很低,这正是实验里电树枝通道发黑、尖端发亮的物理对应。

整个耦合系统的逻辑就是:电场分布驱动s演化(电场能进入自由能),s演化反过来改变电导率和介电常数,进而改变电场分布。如此反复迭代,树枝就长出来了。

2.3 无量纲化的必要性

这里必须提醒一句:相场各参数如果直接用国际单位,收敛性和数值尺度会把人折磨疯。演化和电场方程耦合时,界面宽度通常是微米级甚至纳米级,而击穿阈值电场是kV/mm级,两者量纲差八九个数量级,直接求解非线性方程几乎不可能稳定收敛。

我在建模时把长度、时间、电场都做了归一化。参考长度选l0 = 1 μm(特征微结构尺度),参考电场选E0 = 100 kV/mm(材料介电强度量级),时间尺度由迁移率M和特征长度共同决定。归一化之后,迁移率M和界面宽度κ都落在0.1到10之间,数值上温和得多,求解器也好伺候。具体比例因子等于给整个方程除了一圈标度,Comsol里直接用“缩放”功能就行。

3. Comsol实操:一步步把模型搭出来

3.1 物理场选择

Comsol里我实际用到的物理场只有两个:一个是“电流(ec)”或“静电(es)”,另一个是“系数型偏微分方程(c)”。如果你用最新版的“PDE通用形式”,可以直接输入Allen-Cahn方程,但要记得勾选“瞬态求解”。通行的做法是在Global Definitions里定义变量s和V,然后在PDE物理场里把s设成因变量,在ec物理场里把V设成因变量,两者通过源项耦合起来。

我个人的偏好是用“电流(ec)”而非“静电(es)”,原因上面说了:击穿通道的等电位行为靠电导率突变实现更稳定,纯静电的介电常数突变容易在界面上造成电荷奇点。

3.2 生成多晶几何的两种路径

多晶结构建模是整个模拟里最花时间的前处理。最简单的路径是直接在Comsol里用“生成Voronoi图”功能:在二维平面随机撒点,每个点作为晶粒种子,生成Voronoi镶嵌,再给每个多边形赋予不同的晶粒编号。操作上,Comsol的“模型向导”里有“晶粒(Grain)”几何特征,配合“随机种子”就能生成多晶结构。

但如果你对晶界厚度有定量要求(比如想要晶界厚度100 nm vs 1 μm的对比),我更推荐用MATLAB生成Voronoi图后导入。做法是:在MATLAB里用内置的voronoi函数生成顶点和边,再通过Comsol的CAD导入接口生成面域;或者直接编程计算晶界区域,把晶界和晶粒分成不同材料域。纯Comsol几何也能做,只是控制晶界厚度比较绕。我在本机用python的scipy.spatial.Voronoi生成顶点坐标,然后通过Comsol LiveLink for MATLAB把几何数据传进去,这套流程我稳定复现了十几次。

3.3 晶界怎么赋参数

多晶Voronoi图生成后,晶粒之间天然共享边界,但座标上并没有厚度概念。要让晶界成为一个独立的弱化区,有两种处理方式。一种是直接在几何上把晶界建造成一条细带,比如给每个共享边界偏置出一个宽度为1 μm的矩形区域,单独给它分配低介电强度、高电导率的材料属性。另一种更“相场”一点,不显式建晶界带,而是在相场模型里额外加入一个晶界位势场G(x,y),当点落在晶界附近时,G=1,此时势垒γ降为0.3倍甚至更低,电场能耦合项也相应放大。

我强烈推荐第二种,原因有两个:一是不需要额外处理几何布尔运算,建细带容易产生大量狭窄网格单元;二是晶界势垒的宽度可以通过高斯函数平滑控制,不会引入非物理的锐利界面。

3.4 边界条件设置

电树枝模拟的经典“针-板”结构:上方是高压针电极,尖端的曲率半径非常重要(我用的是5 μm),下方是接地平板电极,中间是厚度约100 μm的多晶介质。左右边界设为电绝缘,即法向电流密度为零。

针尖附近的电场集中是树枝起始的关键,如果针尖曲率太钝,起始场强不够,相场势垒过不去,树枝根本不会发芽。我另外做了一个小姑子实验:把针尖改成半径1 μm和10 μm两组对比,前者起始时间明显提前,树枝也更细密。网格上针尖附近必须做局部加密,元素尺寸控制在界面宽度的1/3以内,否则前期不发芽、后期突然爆长。

3.5 求解器配置

方程是强非线性的瞬态问题,默认的求解器大概率不收敛。我给的是手动设置:

  • 时间步进用BDF隐式方法,初始步长取归一化时间的1e-4,最大步长不超过总模拟时长的1/100
  • 非线性求解器用牛顿法,阻尼因子初始设为0.1,打开“自动稳定”
  • 如果使用电流物理场,注意电导率跨越多个数量级(σins=1e-14,σbr=1e-4),必须把对数标度切换到变量或采用分段线性插值,否则Jacobian矩阵条件数爆掉

网格上,除了针尖和晶界附近加密,其余的网格粗一点就行。多晶结构本身网格数量不小,100 μm见方的模型跑下来大概3~5万个自由度,单次计算时间在十几分钟到一小时,完全可以接受。

4. 关键参数标定与稳定性技巧

4.1 参数对照表

为了方便直接套用,我把一套能跑通多晶电树枝的参数列出来。单位都是归一化后的,具体物理量级表后说明。

参数符号数值说明
势垒高度γ1.0晶粒内部稳态势差,越大约难击穿
晶界势垒因子λgb0.3晶界处γ乘以该因子,弱化程度
界面能系数κ0.02控制树枝的界面宽度,默认0.02
迁移率M0.05相变速率,决定生长速度
完好相电导率σins1e-6归一化量级,对应绝缘
击穿相电导率σbr10对应导电通道,比绝缘高7个数量级
参考电场E01归一化阈值,实际对应材料介电强度
施加电压Vapp2.5归一化电压,针极需明显超阈值

这里每个参数都牵一发动全身。势垒γ是“多难触发”,迁移率M是“多快扩展”,界面能κ是“树枝长多细”——κ偏大时界面厚,树尖钝,分叉少;κ偏小则长出大量细密侧枝,接近实验上那类“灌木状”电树枝。如果你想突出晶界导引效果,就把λgb调低,同时保持晶粒内部γ不变,你会发现树枝自发地沿晶界绕行,而不是穿切进去。

4.2 时间步与稳定性

这个模型最容易翻车的地方就是时间步。Allen-Cahn方程对界面宽度的要求严格:数值解必须在扩散界面内至少布置3个单元,否则会出现空间振荡,树枝变成锯齿状。可以粗略估算:如果界面宽度l = 0.5(归一化尺度),网格最大尺寸不能超过0.15,不然电场捕捉不到界面处能量的急剧变化。

时间步方面,我习惯开自适应步长,但限制最大增长因子到1.5。初始阶段电场在针尖聚集,s其实变化很小,但如果时间步太大,一个步子里界面会跳过多个网格单元,形成伪拓扑变化;这时树枝看起来像瞬移,而不是连续生长。实测中,初始dt=1e-4到1e-5能稳定跑过前几十步,后面树枝生长进入稳定期,dt可以涨到1e-2。这块建议做一次步长收敛性测试,把dt减半看路径是否变化,如果不变化再放心用大步长。

4.3 电场重分布的逻辑

相场模拟和普通电-力耦合最不同的点在“电场重分布”的反馈方式。当局部s升高时,电导率上升,电流会更快地绕过未损伤区域,于是电场被重新分配:树枝内部的E迅速下降,而尖端前方的E被抬高。这个前置电场增强区域就是下一次击穿跃迁的位置,电树枝的分叉也是因为尖端周围出现两个不等价的增强点,哪个先突破势垒哪个方向就优先生长。

用多晶结构来看,晶界本来就是电导率偏高、介电强度偏低的区域,树枝尖端推进到晶界附近时,前锋场强超过晶界的局部阈值,于是沿晶界方向扩展比穿晶方向更容易,宏观上就是贴界走。这个机制用相场模型展现得非常直观,你只要后处理时把晶界位置和s=0.5等值面叠在一起,就能看到树枝和晶界高度重合。

5. 常见问题与排查技巧实录

5.1 树枝死活不发芽

这是新手上手最容易撞的问题,表现是算了几百步s全场还是0,不是模型错,而是电压或势垒没匹配。排查顺序:

  • 先检查针尖最大电场是否超过了击穿阈值。用后处理画一根穿过针尖的电场线,归一化电场超过1.2左右再说。
  • 再看势垒γ。如果γ设得过高(比如5.0),相变驱动力根本不够,可暂时把γ降到0.3验一下能否发芽,能长再调回去。
  • 最后看时间尺度。前期迁移率M太小,时间步下演化缓慢,需要让总模拟时间足够长,而不是把时间步硬调小。

5.2 树枝长成一片“糊状”而非清晰枝状

典型原因是界面能系数κ太大。κ大则界面厚,界面移动更像平滑扩展,电树枝细枝分叉全部被抹平。把κ从0.05降到0.005,再加密过渡区网格,形态立刻改善。另外,网格太粗也是元凶,检查s=0.1和s=0.9之间的单元层数,少于4层就说明网格分辨率不够。

5.3 树枝沿晶界和不沿晶界的判据

很多人碰到的问题是,明明晶界弱化了,树枝还是直穿晶粒,完全没有偏转。原因多半是晶界势垒因子不够低,或者晶界区域在几何上没有真正形成连续介质。我排查时发现,Voronoi导入后晶界带有时是断开的,或者晶界带宽度远大于电场集中区,弱化效果被稀释。一个快速验证方法是单独缩小模型范围,只放一条晶界,通电看树枝是否沿它偏转,能更快定位问题。

5.4 求解不收敛和网格依赖

再稳的模型也会偶尔发散,最常见的报错是“找不到初始解”或“时间步长降到最小”。基本应对套路是:把电流物理场的电导率从指数突变改成平滑过渡,比如σ(s) = σins + (σbr-σins) * s^3,cuadratic变三次有助于Jacobian稳定;或者把PDE里的梯度项系数κ调大一点让界面更宽。如果还是不收敛,把迁移率M调小一个数量级,发散了再处理。

6. 一些个人心得和后续扩展方向

最后聊点实操层面的体会。相场模拟电树枝这件事,技术门槛确实不低,但最大的难点其实不是写方程,而是把材料物理合理映射进参数。不同聚合物、陶瓷的晶界电导率差异可能是数量级的,模拟前最好先通过阻抗谱或者击穿切片实验拿到一些参考数据,再反过来标定γ、σbr和λgb。我自己的经验是先复现一套文献里“没有晶界”的均匀介质模型,确保电树枝生长形态和理想情况一致,再加多晶结构,一次只增加一个复杂度,这样出错时定位容易得多。

如果你后续想把模型引向实际工程应用,有几个方向值得考虑:一是加入热-电耦合,模拟焦耳热对相变势垒的修正,这会让击穿路径更接近真实放电过程;二是把机械应力场也耦合进来,电致伸缩或热应力在多晶界面产生的应力集中会显著影响树枝路径;三是和优化算法结合,比如用相场模拟作为评估函数,反向搜索晶界结构设计,让材料获得更高的抗电树枝能力。我现在就在往热-电-力三耦合方向推进,头很疼但很上头,等有进展了再回来更新,把这套Comsol模型继续往下盘。

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

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

立即咨询