研究锂金属电池的人十有八九都被枝晶搞过心态。我做COMSOL仿真这几年踩过最多的坑就是界面移动和应力耦合叠在一起之后疯狂不收敛。最近这个项目正好把泰森多边形、粉末锂金属负极和应力模型放在一起做了一遍出图效果和物理过程都很满意借这个机会把整个建模思路、物理场设定和调试过程完整记录下来给正在折腾枝晶仿真的朋友一个能直接上手的参考。这个项目解决的问题其实很具体锂金属负极表面在充电过程中会长出枝晶枝晶一旦刺穿隔膜电池就直接短路报废。传统的仿真模型大多把电极表面当成一块平整的金属这当然能定性看出趋势但真实负极尤其是粉末锂金属负极表面是无数颗粒堆起来的颗粒之间还有孔隙、晶界和应力集中点。想回答“应力到底怎么影响枝晶在哪个位置长、长得有多快”就必须把几何形貌和多物理场耦合一起做进去。这也是我选择泰森多边形做电极微观结构、再叠加应力模型的核心原因。整套模型跑下来TSATaylor series analysis不太合适了……好吧不开玩笑直接说项目。所用软件是COMSOL Multiphysics几何上用泰森多边形Voronoi近似负极粉末颗粒的排布电化学部分用Butler-Volmer方程描述锂沉积动力学固体力学部分耦合扩散诱导应力界面移动用ALE移动网格实现最后用类似“海拔图”的方式输出界面形貌随时间的变化也就是大家常说的horizon图。1. 项目整体设计泰森多边形、粉末锂金属负极与应力模型为什么放在一起1.1 锂枝晶模拟要回答的核心问题锂金属负极的理论比容量高达3860 mAh/g是石墨负极的十倍以上但它至今没法大规模商业化核心障碍就是枝晶。模拟锂枝晶生长并不是单纯为了复现“树枝状形貌”好看而是为了回答三个工程问题枝晶在什么位置最容易形核是颗粒尖角、晶界交界处还是孔隙附近沉积过程中内部产生的应力会不会反过来影响局部电流密度从而加剧或抑制枝晶生长如果负极粉末颗粒的初始形貌改变比如粒径缩小、孔隙率提高枝晶风险是否显著下降这三个问题单靠“平整表面均匀电场”的二维模型回答不了。你至少需要一个能体现颗粒间几何约束和应力集中的微观结构模型。这时候泰森多边形就派上用场了。1.2 为什么是泰森多边形颗粒形貌的“抽象但足够真实”泰森多边形是一种空间剖分方式给定一组种子点每个点周围划分出一个区域使得区域内的任意位置到该种子点的距离都小于到其他种子点距离。用大白话说就是把空间按“最近邻原则”切成一个个小房子每个种子点管一块地盘。在材料微观结构模拟里泰森多边形被广泛用来生成多晶和颗粒堆积结构。因为真实粉末锂金属负极中的锂颗粒在压制和电化学循环后会形成类似多晶金属的微观组织颗粒大小不一、边界不规则、颗粒之间有一定孔隙。如果完全按照SEM照片重构几何建模成本极高而且随机性太大很难做参数化研究。泰森多边形的好处是只要调整种子点的密度和位置就能控制颗粒尺寸分布和孔隙率几何边界天然具有“颗粒边界”的物理含义容易与周期性边界条件配对模拟结果可以近似代表宏观电极行为。所以在这个项目里我用二维泰森多边形表示负极截面每个多边形代表一颗锂金属颗粒颗粒之间的窄缝代表孔隙和电解质通道。对于粉末锂金属负极来说这种多颗粒几何比“一整块锂片”更接近真实微观环境。1.3 应力模型不能省沉积诱导应力如何反过来支配枝晶锂沉积过程中锂离子在电极表面还原成锂原子新沉积的锂占据体积周围颗粒会挤压变形。粉末负极看起来全是孔隙似乎有空间容纳膨胀但颗粒接触点附近仍然会产生显著应力。更麻烦的是应力改变了局部的电化学平衡电位导致电流密度分布不均匀电流密度大的地方沉积更快应力进一步集中形成正反馈。应力模型在COMSOL中通过固体力学接口实现核心方程是应力平衡∇·σ Fv 0其中σ是柯西应力张量。对于小变形弹塑性分析应变张量可以分解为弹性应变、塑性应变和浓度膨胀应变而浓度应变项跟锂离子浓度c直接相关通常写成ε_c β (c - c_ref) I这里的β是锂化引起的体积膨胀系数c_ref是初始浓度。这个式子把电化学浓度场和力学场串在一起了是整个耦合模型的关键纽带。如果你不做应力模型相当于默认颗粒在沉积过程中不发生任何变形和力反馈。但真实情况下应力会显著影响等效平衡电位公式近似为E_eq E_0 (σ_n V_m) / (zF)其中σ_n是界面处的法向应力V_m是锂的摩尔体积z1F是法拉第常数。颗粒接触点附近应力大E_eq会改变局部过电位η也跟着变最终影响Butler-Volmer电流密度。这就是“应力模型不能省”的根本原因。2. COMSOL里的模型搭建从几何到物理场2.1 泰森多边形的生成与导入我走过的三条路COMSOL没有内置的“一键生成泰森多边形”功能至少到目前的主流版本都没有。我试过三种方案逐个说下效果第一种MATLAB LiveLink for COMSOL。在MATLAB里用内置的polyshape或voronoi函数生成泰森多边形的顶点和边然后通过LiveLink把几何数据直接传给COMSOL。这个方案最灵活缺点是必须额外安装LiveLink而且MATLAB版本和COMSOL版本匹配容易出幺蛾子。第二种外部工具生成DXF再导入。用Neper或者是Python的scipy.spatial.Voronoi生成多边形顶点坐标再用dxfwrite或matplotlib输出成DXF文件最后在COMSOL中“导入CAD文件”。这个方案不依赖LiveLink是绝大多数情况下的好选择。第三种COMSOL内置几何操作“硬拼”。在COMSOL中手工输入几个圆的圆心坐标再取交线效率极低只适合两三个颗粒的示意图不适合做真实粉末电极模型。我最终推荐第二种。用scipy.spatial.Voronoi生成周期性泰森多边形大概流程是在矩形区域内随机撒N个种子点坐标范围从0到1调用Voronoi函数得到顶点和边裁剪掉超出边界的线段保留闭合多边形导出DXF文件每个多边形闭合在COMSOL中导入再执行“转换为实体”。如果颗粒数量比较多比如50个以上建议在Python里就处理好拓扑关系不要让COMSOL靠CAD内核去自动缝合否则很容易出现“转换为CAD内核时不支持的拓扑”这类报错。2.2 几何检查与拓扑修复别让CAD内核卡住导入泰森多边形之后第一步不是设置物理场而是做几何清理。COMSOL底层CAD内核有时候会把一些共享边识别成两条重复边导致布尔操作失败。我的经验是使用“修复”节点把小于某个绝对容差的短边合并掉把颗粒之间的狭缝统一处理成“孔隙域”不要让某个三角形状的微小孔隙孤立存在如果出现“转换为CAD内核时不支持的拓扑”错误优先改用“从网格生成实体”的办法也就是先用网格划分器生成表面网格再用“网格到几何”重建实体模型。这一步看起来不起眼但非常关键。几何拓扑不干净后面划分网格要么卡死要么出来的网格在尖角处极度畸形直接拉低求解收敛性。2.3 电化学与应力物理场的耦合设定在COMSOL中我用的是“三次电流分布”接口来描述电解液中的离子传输和电势分布用“稀物质传递”接口描述锂离子浓度场用“固体力学”接口描述负极颗粒的应力和变形。三个接口通过边界条件互相咬合。电解液区域遵循质量守恒和电荷守恒核心方程是Nernst-Planck方程简化后的形式∂c/∂t ∇·(-D ∇c - z u F c ∇φ_l) 0其中D是锂离子扩散系数u是离子迁移率φ_l是电解液电位。界面处的锂沉积通量由Butler-Volmer方程控制i i0 [ exp(αa F η / (RT)) - exp(-αc F η / (RT)) ]这里i0是交换电流密度αa和αc分别是阳极和阴极传递系数η是局部过电位η φ_s - φ_l - E_eqφ_s是电极电位φ_l是电解液电位E_eq是平衡电位。涉及应力时需要把上节提到的应力修正项写进E_eq里。电流密度和界面移动速度的关系是v_n -i M / (F ρ)M是锂的摩尔质量ρ是锂的密度。这个表达式直接赋给ALE移动网格的边界法向速度。2.4 枝晶生长的界面表达ALE移动网格模拟枝晶生长有几条路线水平集法、相场法、锐界面追踪法和ALE移动网格法。相场法物理信息最丰富但计算量巨大水平集法适合大拓扑变化但界面宽度和参数选择比较玄学。考虑到我要耦合应力模型用ALE移动网格是更稳妥的方案。ALE的思路是网格的节点随物质界面一起运动但内部网格节点可以做平滑处理以避免过度扭曲。好处是界面位置始终清晰固体域和电解液域的边界保持追踪应力计算比较自然。ALE移动网格的关键设置有三点界面速度必须是法向速度不能是简单的x或y方向速度否则界面会在运动中“脱皮”网格平滑方式选“Laplace平滑”或“Winslow平滑”比“超弹性平滑”更容易收敛界面位移太大时比如枝晶长度超过网格尺寸好几倍必须定期重新剖分网格否则会得到严重畸形的网格单元直接导致雅可比行列式变负。3. 实操过程与出图从求解器设置到horizon图3.1 参数与边界条件清单先给一个可以直接抄作业的参数表。这套参数对应的是二维模型电解质为液态负极材料为锂金属温度设为室温300K。参数数值说明负极初始锂浓度 c_076.6 mol/L锂金属的本体浓度电解液锂离子浓度 c_e1 mol/L典型商业电解液扩散系数 D1e-10 m²/s电解液中锂离子扩散系数交换电流密度 i010 A/m²典型液态电解液数值阳极传递系数 αa0.5对称传递阴极传递系数 αc0.5对称传递锂摩尔体积 V_m1.3e-5 m³/mol锂的摩尔体积锂摩尔质量 M6.94e-3 kg/mol锂的相对原子质量锂密度 ρ534 kg/m³锂金属密度杨氏模量 E7.8e9 Pa锂金属模量泊松比 ν0.36锂金属泊松比初始形核过电位0.05 V界面拉入/拉出量因为沉积层只有几百纳米到几微米初始几何里在颗粒界面上预设一个“小凸起”作为形核种子。这个小凸起会造成初始电流密度集中后面枝晶的发育方向就会从这里开始。如果没有预设形核点界面演化会完全被网格噪声主导不容易复现实验结果。边界条件方面模型最左侧是集流体边界设置恒定电位φ_s 0最右侧是对称边界或电解液本体边界设置φ_l 0上下边界用周期性边界条件模拟无限大电极。泰森多边形颗粒内部是固体力学域颗粒与电解液接触界面是电化学活性边界赋Butler-Volmer电流密度和法向移动速度。3.2 网格划分与求解器收敛调试网格是这套模型最容易翻车的地方。我的策略是全域用自由三角形网格但在颗粒/电解液界面上添加边界层网格保证浓度梯度和界面移动的解析精度。界面附近的网格尺寸设为0.02 μm而颗粒内部的网格尺寸放宽到0.2 μm这样可以显著减少计算量。求解器设置上我直接忽略瞬态求解器的自适应时间步长默认值自己手动设置最大时间步长不超过0.01 s。因为枝晶生长过程中局部电流密度可能在几微秒内剧烈变化自动时间步长往往为了满足整体误差把步长放得太大导致界面速度在时间轴上“跳步”最后得到锯齿形枝晶轮廓。另外应力场和浓度场需要同时迭代。COMSOL的默认分离式求解器有时能跑通但经常会因为浓度场剧烈变化导致固体力学场振荡。我建议改成全耦合求解器打开“恒定牛顿”或“自动牛顿”选项。如果内存吃紧就把“阻尼因子”从1降到0.7左右收敛稳定性会好很多。3.3 “horizon图”怎么做时变形貌图和应力分布图“horizon图”这个词我刚接触时也懵了一下其实它对应的就是把电池界面看成一个地形表面把沉积厚度或者界面位移当“海拔”画出来的投影图。你可以把它理解成一张“界面海拔图”能直观看到枝晶在哪个位置长得高、在哪个位置塌陷。在COMSOL后处理里我的做法是在“二维绘图组”里新建一个“表面图”表达式选择move变量里的displacement的y方向分量或者直接选择ALE移动网格的边界位移勾选“高度表达式”把这个表面图变成带海拔起伏的3D图再新建一个“颜色表达式”用锂沉积厚度或von Mises应力做颜色映射。出来的图形就是一块起伏的“地形”颜色代表应力或浓度z方向高度代表界面位移这就是我要的horizon图。这种图比单纯看浓度云图高级在哪儿它能同时呈现形貌演化和应力分布两个信息。比如颗粒接触点附近应力很高同时界面形貌显示那里的沉积厚度反而偏薄那说明应力抑制了沉积这个现象在锂金属负极里非常重要会直接决定枝晶是否会穿透颗粒边界。3.4 结果怎么解读从浓度云图到应力演化我跑完一个典型的37颗粒泰森多边形模型后第一感受是应力集中在颗粒接触点和尖角处而不一定在电流密度最大处。这就引出一个很有意思的结论粉末锂金属负极虽然整体表面积大单位面积电流密度小但颗粒尖角处的应力集中可能反而成为枝晶形核的“帮凶”。从horizon图上看枝晶并不是从所有颗粒表面均匀长出来的而是优先在几个应力相对较低的孔隙通道中发育。这说明应力对枝晶的路径选择有很强的调控作用。换句话说不能只靠“降低电流密度”来抑制枝晶还需要从颗粒形貌和力学约束角度去设计负极结构。4. 常见问题与避坑实录4.1 “转换为CAD内核时不支持的拓扑”错误这个报错在导入泰森多边形后非常常见。根本原因通常是DXF文件里的多边形边界在共享边处没有精确重合。CAD内核在转换时无法判断哪些线段应该缝合从而报错。处理办法有三个层次在Python导出DXF时把每个顶点的坐标四舍五入到小数点后8位确保共享点完全一致在COMSOL导入后打开“修复”节点删除短于1e-9 m的边把细小裂缝忽略掉最暴力但有效的办法隐藏几何建模改用“导入网格”方式先生成表面网格再重建几何。4.2 弹塑性应变变量在迭代未收敛如果你在COMSOL里用了塑性模型可能在迭代过程中碰到“用于查找弹塑性应变变量在迭代未收敛”这类提示。这通常不是模型设置错误而是塑性本构更新算法在局部单元上发散。我的解决办法是降低最大牛顿迭代步数或使用阻尼牛顿法把塑性参数中的硬化模量稍微提高一点这是因为纯理想塑性模型在接触点容易导致局部奇异性检查是否其实不需要塑性模型如果应力远低于屈服强度只做弹性就够了别为了“完整”硬加塑性。4.3 泰森多边形的尖角导致网格畸变泰森多边形经常会出现角度极小的尖角尤其是种子点距离太近时。这些尖角对几何求解和网格划分都是噩梦。处理策略是在生成种子点时加一个“最小距离限制”有点像分子动力学里的排斥势让任意两个种子点之间距离不小于某个阈值比如颗粒平均粒径的20%。如果你不想重跑Python脚本也可以在COMSOL里对尖角做“倒角”或“圆角”特征。不过对大量颗粒来说手动处理不现实。建议还是在几何生成阶段就把分布控制好。4.4 电脑配置低怎么跑大型模型泰森多边形颗粒数量超过30个以后三维模型基本就不是普通电脑能跑的了。我通常做两层近似一是把三维模型压成二维截面虽然会损失一部分厚度方向信息但枝晶生长的主趋势依然成立二是只选取局部代表性区域比如3×3个颗粒的中心区域做了周期边界来捕捉颗粒接触点的局部应力。如果必须要跑三维建议减小模型尺寸为单个颗粒加半个孔隙的代表体积单元同时关闭“极大网格”选项使用“较细化”。我实测过三维模型颗粒数超过50个后内存消耗会成倍增长求解时间长达数天已经不适合参数扫描。5. 写在最后的个人建议与扩展方向5.1 这个模型还能怎么玩泰森多边形加应力耦合这套框架不只是能模拟锂枝晶稍微改改物理场就能迁移到多个方向如果研究固态电解质可以把孔隙域替换成固态电解质域并在界面处加上接触阻抗项观察晶界处的锂镍渗透路径如果研究SEI膜可以在颗粒外加一层薄壳域给壳域单独赋杨氏模量和断裂准则看应力积累到多少会让SEI破裂如果研究负极体积膨胀对全电池的影响可以把泰森多边形颗粒放在一个更大的电极单元里通过均匀化获得应力-应变响应。5.2 我踩过坑沉淀下来的几条守则最后分享几条个人经验算不上什么高深理论但确实能救命第一泰森多边形模型的灵魂在几何清理不在物理场设置。几何不干净后面所有操作都会受到牵连。我一般会花一半精力在处理几何上。第二应力模型要胆大心细。别一上来就全耦合先跑通浓度场和电化学场再把应力场加进来。分步加物理场至少能让你定位是哪一步引入的不收敛。第三horizon图是很好的展示工具但不要被它骗了。颜色范围如果自动变化两张不同时刻的图看起来会差很多。出图时务必固定颜色标尺的范围否则视觉上的“枝晶长高”可能只是标尺缩放了。第四不要吝啬边界层网格。枝晶生长的核心驱动力来自界面附近的浓度梯度和电位梯度界面解析精度不够后续应力场再准也白搭。这个项目做完之后我对粉末锂金属负极的理解比只看文献深刻了很多。泰森多边形并不完美但它给了我们一个把“颗粒形貌”“孔隙结构”“应力集中”和“枝晶路径”统一进来的最小模型框架。如果你正在做类似方向希望这篇记录能帮你少走一些弯路。