各向异性NMC颗粒脱锂浓度、应力与晶界损伤有限元模拟
文献【Shaofeng Li, Zhisen Jiang, Jiaxiu Han，et al. Mutual modulation between surface chemistry and bulk microstructure within secondary particles of nickel-rich layered oxides. Nature Communications, 2020, 11: 4433】通过将软X射线纳米探针（表面Ni价态表征，~5 nm探测深度）与硬X射线全视场透射显微镜（体相三维微结构重构，~30 nm分辨率）相结合，在单个NMC811颗粒上实现了表面化学与体相微结构的关联表征。
图1 成像装置与原理示意图：TXM全视场成像（探测体相）与软X射线纳米探针扫描（探测表面，~5 nm深度）
同时，文献还建立了二维模型模拟NMC颗粒脱锂过程中Li浓度、应力和晶间断裂损伤的协同演化过程。
图2 有限元模拟锂浓度分布与晶界损伤在20 s、400 s、720 s的演化
文献模拟结果动画展示
文献中具体的模拟方法为：将NMC 二次颗粒简化为一个圆，内部采用 Voronoi方法生成大量随机尺寸和取向的多边形一次颗粒，以模拟真实多晶 NMC 的微观结构。采用各向异性的 Fick 定律模拟锂的扩散过程，设定ab 面层内扩散系数 是 c 轴方向的10 倍。初始 Li 浓度设为 NMC 的最大浓度 Cmax（对应完全锂化态）。
为模拟颗粒表面不同区域与电解液/导电网络接触程度的差异，在颗粒表面施加非均匀的恒定浓度边界条件：将颗粒表面划分为四个象限，各象限分别固定为：0.2Cmax、0.3Cmax、0.5Cmax、0.8Cmax。这种设置是对颗粒表面电化学条件不均匀性的简化模拟，例如部分颗粒表面被粘结剂覆盖、或与导电碳接触不好，导致局部脱锂速率不同。
采用线弹性材料本构模型计算弹性和脱锂诱导应变过程。力学边界条件设定为：颗粒中心施加零位移约束；颗粒外表面为自由边界。
采用内聚力模型（Cohesive Zone Model, CZM） 模拟晶界处的裂纹萌生与扩展。定义损伤函数 d ：d=0为弹性响应阶段，界面完好；d=1为应变能释放率超过断裂能，界面完全失效（开裂）；0<d<1 是线性插值的过渡损伤状态。
参照本文的模拟方法，模拟设置和原始文献不完全一致。本人采用COMSOL软件复现了文献中的脱锂浓度、应力、晶界损伤的有限元模型。首先展示模拟结果：
本人模拟结果动画展示：锂浓度cs/csmax，扩散系数各向异性不考虑各晶向力学异性
本人模拟结果动画展示：应力和晶界损伤，扩散系数各向异性不考虑各晶向力学异性
本人模型中，以整个颗粒的容量为基准，颗粒表面以一定倍率的电流发生脱锂反应，模型考虑了颗粒表面电化学反应速率不均匀的现象；表面不断脱锂，颗粒内部的锂根据菲克定律向表面扩散，颗粒内的一次晶粒具有不同的取向，模型考虑了锂扩散的各向异性，ab方向层内扩散系数比层间c轴方向扩散更大。
通过COMSOL中插层应变模型模拟颗粒脱离过程中的应变，由于锂浓度不均匀产生不均匀应力。考虑弹性应变和脱锂本征应变，考虑各晶向不同的弹性模量。采用内聚力模型模拟晶界处的损伤。
下面是该模型详细的建模过程，如果感兴趣可继续付费阅读（可提供模型comsol源文件）。如果不想付费阅读可退出，或者也可以根据文献的模拟方法来自己建模。
第一步，采用Voronoi算法生成圆形颗粒内的多边形一次颗粒。我采用在线工具(websvg.github.io/voronoi/)生成了初步的几何模型，半径5μm的圆形颗粒内包含13个多边形一次颗粒。然后，自己设置了如图所示的各个晶粒的晶向，二维几何中晶向轴分别为a和c，a轴表示层内方向，c轴层间方向。后面的扩散系数和弹性模量各向异性根据此晶向来设定。
Voronoi生成在线工具界面
NMC圆形颗粒及晶粒的晶向
