做细观力学模拟这些年Voronoi镶嵌一直是我工具箱里离不开的一个几何工具。尤其在做混凝土、陶瓷基复合材料这类多相材料时用二维Comsol把Voronoi边界设置好就能快速生成接近真实骨料分布的随机多边形骨料再搭配纤维骨料模型几乎可以覆盖绝大多数的细观尺度分析需求。这篇内容我就把自己常用的整套流程摊开讲包括Voronoi边界怎么设、多边形骨料怎么提取、纤维怎么嵌进去以及我实际踩过的坑和处理办法。1. Voronoi在骨料建模中的核心思路与设计选择1.1 为什么二维问题选Voronoi从真实细观结构到数值模型真实混凝土或岩石的骨料分布并不是人工排列的而是由地质过程、搅拌过程随机堆积形成的。直接对真实截面做CT扫描再矢量化的确能还原细观结构但成本高、普适性差而且每换一组级配就要重新扫描一次。Voronoi图从数学本质上就对应了“在空间里按距离竞争划分区域”的过程——每个点控制一块领地领地边界是相邻两点连线的垂直平分线。这与真实骨料在浆体中相互“抢占空间”的堆积逻辑非常相似所以用它来生成二维骨料模型既有随机性又有物理合理性。从计算角度看Voronoi镶嵌还有一个额外的好处它天然把区域剖分成不重叠的凸多边形集合。凸多边形的几何质量普遍比凹多边形好后续划分网格时少出现尖角和细长面元求解器收敛也更省心。二维问题里Voronoi多边形的边数通常在5到8之间这与真实细观截面上骨料颗粒的等效多边形轮廓也比较接近。尤其是当你要研究裂纹扩展路径时Voronoi骨料的棱角分布会迫使裂纹沿着界面或穿过骨料产生自然的分支和偏转这比画一堆圆形骨料要真实得多。在软件选型上二维问题选Comsol而不是其他有限元工具主要看中它的几何序列非常灵活可以直接用坐标表或外部脚本生成Voronoi点源然后通过几何编辑构造多边形材料分配、网格划分和后续后处理全在一个序列里完成。对侧重机理分析而非大规模统计的研究场景这种工作流明显更顺手。1.2 多边形骨料与纤维骨料的分工两种细观单元的建模差异很多刚开始接触细观建模的人会有一个误区以为把多边形骨料做好就够了纤维只是可选项。实际上这两类骨料在物理机制里承担的角色完全不同。多边形骨料主要提供刚度骨架和体积填充。在混凝土里粗骨料的体积率能到40%到50%它的弹性模量通常高于水泥砂浆受力时会把大部分载荷通过骨料骨架传递同时在骨料与砂浆的界面过渡区产生应力集中成为裂纹萌生的温床。所以多边形骨料建模的重点是形状随机、尺寸符合级配、体积率准确、界面处理得当。纤维骨料则是增韧和拉结的关键。钢纤维、PVA纤维或碳纤维在基体里体积占比往往只有0.5%到2%但引入后能将材料的峰值后韧性和抗裂能力提升一个量级。纤维的长度一般在6到30毫米直径只有几十到几百微米长径比从几十到上百如果用实体单元去划分网格纤维与基体之间的尺寸跨度会造成灾难性的网格数量膨胀。因此纤维建模的核心思路是降维处理用线单元或梁单元代表纤维通过嵌入约束或内聚力模型把它与基体耦合起来只保留纤维轴向力学行为的关键信息。这两类骨料在一个模型里共存时几何构建的难度会显著上升。多边形骨料是面域纤维是线域线要穿过面、落在面边界上或者完全在基体里每一种情况对边界设置的要求都不一样。后面我会专门讲怎么把这两类单元在同一几何序列里编排避免出现“骨料盖住纤维”或者“纤维悬空”的低级错误。2. Voronoi边界设置完整流程基于Comsol2.1 第一步生成Voronoi点源与几何拓扑Voronoi生成的前提是有一组离散点点的分布直接决定了骨料形态和空间均匀性。我在实际工作里常用的点源生成方式有三种。第一种是最简单的随机均匀点在矩形或圆形区域里用Matlab、Python或Comsol自带的随机函数生成坐标点数量由目标骨料数量决定。这种方式快但会出现局部点过密或过疏生成的Voronoi多边形尺寸差异很大如果是为了做级配研究需要后续筛选和调整。第二种是泊松圆盘采样点与点之间的距离不小于给定最小值这样生成的Voronoi多边形尺寸更均匀骨料之间不会出现极端的狭窄间隙。在Comsol里可以用“全局常微分和微分代数方程”配合随机数辅助生成或者直接在外部生成坐标后以表格形式导入。第三种是分层级配方式先确定目标级配曲线比如Fuller级配或者等效面积级配把不同粒径范围的骨料颗粒数和面积分数算出来再按面积权重生成点。这种方式最接近真实混凝土但前置计算量稍大适合要做定量体积率匹配的研究场景。点源坐标生成后我自己习惯在Comsol里头用“几何 → 转换 → 多边形”的方式快速手动描边但如果骨料数量多手描不现实。更高效的做法是用外部脚本先算好Voronoi图的顶点坐标和边连接关系然后以DXF或文本坐标表的方式导入Comsol再通过“转换为曲线→转换为实体”构建多边形。注意COMSOL版本不同“曲线转实体”菜单路径可能有差异但原理不变——先用样条或多段线把Voronoi边轮廓画出来再把这些闭合轮廓转换成面域。建议在导入前先清理掉重复点和自相交边否则后续布尔操作会报“几何不可修复”的错误。2.2 第二步构建多边形边界与布尔操作拿到Voronoi轮廓后第二步是把骨料从完整镶嵌里分出来。常规做法是把整个Voronoi镶嵌当作一个覆盖全域的多边形集合然后用“布尔操作 → 分区”把区域拆成两大部分骨料域和多边形之外的基体域。具体操作上我会先把所有Voronoi单元按粒径大小分组。例如目标粗骨料最大粒径是20毫米那么面积等效直径大于10毫米的保留为骨料小于10毫米的填充到基体区域里当作细骨料或砂浆的一部分。这一步直接用“几何 → 选择 → 按面积筛选”是做不到的我通常用外部脚本先导出每个多边形的面积和边长筛选好之后再导入Comsol。布尔操作的常见坑是多边形之间如果存在共享边或者说完全邻接Comsol的“并集”“差集”操作有时会出现边界重叠警告。解决方法是引入极小容差偏移比如把骨料多边形略做0.001毫米的内缩或外扩保证域边界不会完全贴合。2.3 第三步材料边界与物理场边界的一致性几何边界设置好后接下来要把材料边界和物理场边界统一起来。一个常见的错误是只给域分配材料却不检查边界是否封闭导致后续求解时部分边界被当成开放边界结果出现异常泄漏。在静力学分析中基体与骨料的界面通常用“连续性”默认绑定如果你是研究界面脱粘就需要把共享边界单独建立成“接触边界对”或“内聚力边界”。这要求在几何序列里为每一条共享边定义命名选择不能只用域选择。传热或扩散问题里Voronoi骨料的边界还可能涉及界面热阻或渗透率的跳跃这时候要在“边界 → 边界条件”里指定不连续条件而不能简单用连续条件带过。我在做电池电极的多孔结构模拟时界面参数的影响非常大这一块设置不对温度场和浓度场就能差出几个百分比。2.4 边界处理3个关键参数与隐藏坑参数名推荐范围或设置隐藏坑最小多边形边长大于网格最小尺寸的2到3倍太短会生成畸形面元求解器卡死相邻骨料间隙大于最小单元尺寸间隙过小自由网格划分时直接报“边界层重叠”Voronoi多边形圆滑角0.05到0.2毫米过渡圆角尖角处应力奇异但圆角过大又会失真这里重点说第三个圆角的处理。Voronoi多边形天然是尖锐转角在真实骨料中棱角会有一定的圆化。如果直接保留锐角有限元分析时尖角处应力会趋于奇异数值上表现为“应力集中随网格加密而不断增大”这并不是真实材料行为。我一般用Comsol里的“圆角”或“倒角”功能给每个多边形顶点加0.05到0.1毫米的圆弧。注意圆角半径不能太大否则大面积改变了骨料的形状和体积率也不能太小否则形同虚设。3. 多边形骨料建模与物性参数分析方法3.1 从Voronoi单元到骨料颗粒如何“保形”及粒径调整Voronoi图直接生成的单元形态不完全等同于真实骨料的级配分布。真实骨料的级配指的是不同粒径区间颗粒的数量比例和面积比例而Voronoi图由点源密度决定两者之间需要一个转换过程。我的处理方法是先定义目标级配区间例如筛孔尺寸2.36、4.75、9.5、16毫米对应四个粒径区间。然后按面积等效原则将每个Voronoi多边形的面积换算成等效圆直径落入对应粒径区间。对于不满足级配要求的单元有两种调整方式一是修改点源密度和分布重新生成二是对已有Voronoi单元做局部缩放或合并。“局部缩放”在实践中很容易引入几何重叠两个相邻多边形几乎共享边界其中一个缩放变大必然挤压到相邻单元。所以我的建议是优先采用“重新生成筛选”策略而不是后处理缩放。生成500个候选Voronoi单元从中按级配和体积率筛选出200到300个剩余区域自动成为基体这个方法稳定可靠。3.2 界面过渡区与接触边界处理多边形骨料与基体之间的薄层区域在水工和混凝土细观力学中常被称为界面过渡区ITZ。ITZ的厚度一般在10到50微米范围远小于骨料尺寸直接按实体建模会带来极大的网格规模。我常用的简化手段有三种第一种是“等效界面属性法”不用单独画出ITZ层而是将骨料与基体界面的粘接刚度折算成等效薄层属性在Comsol的“薄层”边界特征里输入厚度和弹性模量。这种方法最省事适合只关心宏观响应的分析。第二种是“内聚力边界法”沿骨料边界设置内聚力模型定义最大拉应力和断裂能在载荷作用下允许界面发生脱粘和滑移。这个方法能模拟裂纹沿界面扩展的过程但需要同时使用统计分布函数来设置界面强度的随机性否则所有骨料会在同一时刻一起脱粘不符合实际。第三种是“有限厚度环法”对关键区域比如局部裂纹要穿越的路径单独生成ITZ薄环实体其他区域用等效边界代替。混合方法兼顾精度和成本适合做裂纹路径的详细研究。经验之谈如果你做的是单轴压缩或拉伸下的混凝土细观破坏界面强度的Weibull分布参数尤其重要。我常用平均拉应力2到4兆帕、形状参数6到10来标定界面强度效果比均匀强度模型好很多。3.3 多边形骨料的后处理与接口提取模型算完不等于工作结束。多边形骨料模型的后处理重点通常集中在三个位置骨料内部最大主应力、骨料与基体界面上的最大剪切应力、基体中的拉应力集中区。在Comsol的后处理中我会先建立“探测点组”在每个骨料的形心处布置探测点提取应力和应变时程。对于随机骨料不同骨料的应力水平差异很大简单看云图容易忽略局部最不利骨料。用“派生值 → 全局计算”把所有骨料的最大主应力排序就能准确找到起裂关键骨料位置。如果研究的是裂纹扩展我还会输出“损伤变量”或“塑性应变”的等值线把骨料边界处的损伤分布单独用“体素图”或“切片图”表达方便后续生成论文插图或工程报告。这里提醒一点后处理里最好固定颜色标尺不然不同载荷步之间的云图颜色不对应看不出裂纹扩展过程。4. 纤维骨料建模与耦合分析4.1 纤维的几何构建线实体与实体折中纤维骨料的建模和多边形骨料很不一样。多边形骨料是面域纤维是长条细杆如果严格按真实几何建模一根直径0.2毫米、长20毫米的钢纤维与尺寸为几十毫米的基体模型放在一起网格尺寸比能达到几十倍甚至上百倍计算效率直线下降。因此我日常建模会遵守一个准则如果纤维的直径方向不是研究重点一律用线实体边界或边代替如果必须考虑纤维的弯曲刚度或拔出效应再用梁单元或三维实体但这时模型规模必须严格控制。二维线实体纤维的做法是在基体域内画一条线段或样条曲线线段两端嵌入基体边界或穿过边界延伸到外部。在Comsol的“几何”里直接添加“贝塞尔多边形”或“线段”即可。对于随机分布的纤维我常用Matlab或Python生成起点、终点和偏转角坐标再一次性导入Comsol。纤维数量多时建议在外部脚本里直接生成几何对象然后批量导入不要在GUI里一条一条画。4.2 纤维-基体界面的力学/热学耦合纤维建模的下一个核心问题是纤维和基体怎么耦合。最简单的方式是“共用节点”把纤维线实体嵌入基体网格时通过“一致网格”保证纤维节点与基体节点重合。这种方式适合纤维和基体完全粘接、不发生滑移的情况但会带来网格限制纤维线上节点的位置必须与基体网格对齐。更灵活的方式是“嵌入约束”用Comsol的“装配”功能让纤维作为独立几何组件嵌入基体然后在两者之间设置“对”的接触或内聚力。这种方式允许网格各自划分先做基体网格再做纤维线网格再用一致对关联。计算上需要多花一点时间做映射但模型几何独立后期调整纤维数量或分布非常方便。热学分析中的纤维耦合类似纤维可以作为高导热通道在基体中形成热桥此时需要在纤维与基体之间设置界面热阻否则热量会在界面上出现不连续。这里同样可以通过“薄层”边界特征输入界面热导率不需要单独画实体。4.3 纤维多边形联合模型的分析策略当多边形骨料和纤维共存时建模顺序非常关键。我的标准顺序是先生成多边形骨料并把基体域建立好在基体域内生成纤维线注意纤维不能穿过多边形骨料内部否则会造成几何自交或布尔失败纤维线生成后用“分割”命令把所有重叠边界切成多个分支保证每个分支归属于单一材料域最后统一分配材料属性。纤维穿骨料的处理是实操里最常见的失败点。真实混凝土中纤维很难穿过粗骨料绝大多数都绕开骨料分布在砂浆中所以建模时应该对纤维的生成范围进行限制。方法是在生成纤维的脚本里先排除骨料对应的多边形区域再随机生成线段。联合模型的分析参数上纤维体积分数的控制公式是纤维体积分数 纤维总面积 / 基体总面积。用线实体代替时纤维的“面积”需要按等效直径折算例如一根长度L、等效直径d的钢纤维在二维模型中贡献面积近似为 L × d。因此在设置纤维数量时要先把目标体积分数换算成目标总长度。5. 常见报错、网格问题与效率优化速查5.1 网格划分失败的典型原因Voronoi骨料模型的网格失败一半以上都是几何问题而不是网格算法问题。最常见的报错包括“无法生成网格检测到输入几何的自相交”、“边界过度狭窄”、“网格尺寸小于最小单元尺寸”。多边形骨料模型的自相交通常来自布尔操作残留的微小重叠区域。解决方案是在做布尔前先检查每条边是否有微小间隙或者用“合并”特征把所有共享边界统一起来。另一个根源是骨料多边形本身有凹角和非常小的内角尤其是相邻Voronoi单元靠得太近时共享边和内角会变得特别尖锐。网格划分时最好先设置“最小单元质量”检查把边长比大于0.1或最小角小于10度的单元过滤掉。如果发现过滤后网格数量依然很大优先回来检查几何而不是硬调网格参数。5.2 计算效率优化与参数化扫描建议骨料数量在200到500个之间、纤维数量在50到200根之间的二维模型网格规模大概在20万到100万自由度。这个规模在桌面工作站上可以求解但如果你要做多次参数化扫描效率问题就会很突出。我的优化方法有四个一是把纤维模型中的梁单元单独归类在主求解前先做一个小规模的“代表性验证”确认纤维的嵌入约束设置正确避免大规模计算到一半才发现纤维和基体根本没耦合。二是对骨料内部网格做“粗分”基体网格做“细分”。因为骨料内部的应力梯度相对平滑可以适当放宽界面和基体的局部区域才需要加密。三是在扫描不同纤维体积分数时尽量复用几何网格只改材料属性或边界条件。Comsol里“研究”和“几何”分离的设计正好支持这一套流程。这样每次扫描只需要更新求解器不用重新剖分网格。四是如果分析类型是稳态线性问题尽量选择直接求解器加多重网格预处理器的组合。对二维问题这种组合的收敛速度快于迭代求解器在内存占用上也能接受。5.3 我踩过的三个坑实战经验帖第一个坑是“纤维穿过骨料”。我在做钢纤维混凝土试件模拟时随机生成的纤维有一部分直接穿过了粗骨料多边形。添加物理场后纤维节点和骨料节点在交叉点出现了几何体积重叠网格划分直接失败。后来我改成先检测骨料范围再在骨料外区域生成纤维问题就再没出现过。第二个坑是“Voronoi点数过多但筛掉的太少”。做级配研究时我随手生成了1000个点结果筛完之后骨料数量还是太多模型规模膨胀单次求解就要半小时。后来我养成了先把目标级配曲线和面积分数算清楚反推点源数量生成500个点左右筛出200到300个计算量一下就下来了。第三个坑是“周期性边界条件的共享边没设置好”。做代表性体积单元分析时我希望模型左右边界是周期性对应的但Voronoi多边形在左右边界上的切割方式不同导致周期边界上的节点无法一一对应。后来我直接用Comsol的“周期网格”功能先让几何边界匹配再指定周期节点映射这个问题就解决了。如果你要在Voronoi模型上做RVE分析这一步最好提前规划而不是等网格划分完再补救。做细观模拟这些年我越来越体会到Voronoi边界设置和骨料建模本身并不难难的是每一步都保持几何一致性和物理合理性。多花一点时间在几何准备上把边界条件、材料分区和网格策略提前想清楚后续的计算和后处理会顺畅很多。希望这篇分享能让你少走几个弯路。
