做混凝土数值模拟的朋友应该都有同感宏观模型算起来快但很多宏观现象解释不了。裂缝从哪来、怎么穿过骨料还是沿界面走、砂浆和骨料的协同变形怎么分配——这些问题都得往细观尺度去看。我一直在做这方面的研究最近完成了一个含界面过渡区的多面体骨料密堆积3D模型整体跑在ANSYS里效果比较理想。这篇文章就把整个建模思路、算法原理和实操流程完整梳理一遍给想做混凝土细观模拟的同行做个参考尤其是刚入手细观建模、打算用ANSYS落地三相模型的朋友。这个模型解决的核心问题是把混凝土从“均质材料”还原成“骨料砂浆界面过渡区ITZ”的三相复合材料并且在三维空间里把不规则多面体骨料尽量密实地堆积起来。相比常见的球形骨料模型它更贴近真实混凝土的细观结构能更好地还原应力分布、裂缝萌生与扩展路径。下面我从方案选型、算法实现、ITZ建模、ANSYS实操到问题排查一步步展开。1. 混凝土细观建模的整体思路与方案选型1.1 为什么要把混凝土建成细观模型混凝土在宏观层面被近似为均质材料一上细观尺度就完全不是那么回事。粗骨料、砂浆、界面过渡区三者的力学性能差别很大骨料弹性模量普遍在50GPa以上砂浆通常只有20多GPa界面过渡区比砂浆还要弱。这种性能差异直接决定了混凝土受力时的应力分配路径——裂缝容易在骨料附近萌生然后沿着界面或穿过砂浆扩展。宏观模型把材料等效成一个均匀体天然丢掉了这层信息。当你研究尺寸效应、裂缝扩展路径、氯离子在界面区域的传输行为时宏观模型就力不从心了。细观模型把混凝土真实地离散成“骨料砂浆ITZ”三相每相分配独立的几何形状和材料参数。好处非常直接可以看到应力在骨料之间的传递路径可以看到ITZ中的应力集中可以为内聚力模型提供开裂界面甚至可以统计不同骨料含量下的宏观响应差异。做单轴压缩、劈裂拉伸、断裂模拟时细观模型的预测结果通常比宏观模型更贴近试验数据。这也是近几年混凝土细观模拟在论文和工程分析中越来越普遍的原因。1.2 为什么选择多面体骨料而不是球形骨料我刚接触细观建模时也走过弯路最早用的是球形骨料因为程序好写投放算法简单。但球形骨料的缺点太明显球与球之间是点接触骨料在砂浆中的嵌锁效应完全体现不出来真实混凝土里的骨料几乎没有球形的基本都是破碎或浑圆的随机多面体。用球形骨料算出来的裂缝路径往往过于平直应力集中也不够准确和试验数据对不上。多面体骨料的好处在于几何上更接近真实尤其是拉拔和受压模拟中多面体的棱角会造成明显的应力集中和界面剪切错动这恰恰是混凝土破坏机理的关键部分。此外多面体的堆积体积分数可以做得比较高颗粒之间的嵌锁和咬合效应更真实。代价是程序复杂度和计算成本上去了但以现在的计算机性能这个代价完全可以接受。1.3 为什么选ANSYS做落地平台细观建模常用的平台有ANSYS、ABAQUS、COMSOL也有些人用自编程序。我选择ANSYS主要基于三点考虑。第一ANSYS的APDL完全可脚本化几何创建、网格划分、材料赋值、求解设置全部可以用代码控制这对细观模型这种“批量生成、批量装配”的场景特别友好。第二ANSYS的实体单元库比较全SOLID185这种低阶单元在细观模拟中兼具精度和效率非线性求解时比高阶单元更不容易发散。第三Workbench环境里的SpaceClaim和Meshing模块能处理复杂装配体的几何拓扑配合APDL做单元级别的精细控制流程很顺畅。另外ITZ这种薄层结构在ANSYS里的处理方式很灵活。你可以用实体薄层单元模拟也可以用内聚力单元做界面开裂分析选择空间大。这套思路在ABAQUS里同样能做但APDL脚本写起来更直接对纯几何建模需求来说上手更快。2. 多面体骨料生成与密堆积的算法解析2.1 骨料几何的生成方法从球到随机多面体在ANSYS里直接建模生成随机多面体并不方便推荐路径是先用外部程序我用Python生成骨料几何数据再导入ANSYS。numpy处理凸包和分离轴检测非常顺手整个生成流程可以完全参数化。生成单个随机多面体的思路是这样的先定义一颗球半径设为目标等效粒径的一半然后在球面上随机生成一组平面每个平面由球面上的一个随机点作为法向穿球点用这个平面去切割球体保留包含球心的一侧多次切割之后球体就变成了一个随机凸多面体切割次数决定了多面体的面数和棱角程度。切6到10刀得到的是介于八面体和二十面体之间的形状看起来已经有天然的破碎感切太多次就接近球了失去多面体的意义。得到多面体之后要做两个后处理一是按体积等效直径缩放到目标粒径二是做归一化处理让所有多面体在投放时有一个统一的参考尺寸。我这里的归一化基准是“等效球”即该多面体的体积等于一个直径为D的球的体积后续按粒径级配投放时就拿这个等效直径来卡。2.2 级配曲线与骨料数量的确定真实混凝土的骨料级配一般按富勒级配或ASTM标准曲线来配细观模型里常用富勒级配的简化形式P(D) (D / Dmax)^0.5这个式子的含义是累计通过筛孔D的骨料体积占总骨料体积的比例。有了这个分布就可以把设计粒径区间分成若干档逐档计算骨料个数。比如模型粒径范围是5mm到20mm可以分成5-8、8-12、12-16、16-20四档每档取中心粒径对应的体积再按目标体积分数和级配比例算出每档需要投放多少颗骨料。一个容易忽略的细节是三维模型中骨料级配是按体积比例计算的实际投放却按个数进行所以从体积级配换算到个数级配时一定要除以单颗骨料的体积。这会导致小粒径骨料的数量远多于大粒径骨料。假设目标体积分数是50%粒径跨度从5mm到16mm最后统计出来小粒径骨料可能占了几百颗大粒径骨料只有几十颗这是正常现象和真实混凝土的颗粒分布规律一致。2.3 随机投放与碰撞检测的工程实现密堆积的核心算法是随机投放也叫蒙特卡洛砖石法。基本流程是按粒径从大到小的顺序逐个投放骨料每个骨料先在立方体区域内随机生成球心坐标再随机生成旋转矩阵然后检查它与所有已投放骨料是否重叠不重叠就落定重叠就重试。达到最大重试次数后仍然失败就适当减小粒径或调整投放范围后再尝试。碰撞检测是性能瓶颈环节。对两个凸多面体做精确相交判定工程上最常用的是分离轴定理SAT。原理一句话两个凸体不相交当且仅当在某个轴向上它们的投影区间是分离的。对多面体A和B候选分离轴包括A的全部面法向量、B的全部面法向量以及A的每条边与B的每条边叉积得到的向量。逐个轴做投影判断只要发现某个轴上的投影区间不相交就能立即判定为两者不重叠。实际写代码时还需要配合加速技巧。先用球心距离做快速排除如果两个骨料包围球半径之和小于球心距则肯定不相交再用精确SAT检测。对已投放骨料用网格分桶存储空间索引只检测邻近区域内的骨料时间复杂度能从O(n^2)降到接近O(n)。我实测下来这套组合能做到几百颗骨料的投放在一分钟内完成完全够用。3. 界面过渡区ITZ的建模方法与实践选择3.1 ITZ是什么为什么影响宏观力学ITZ的全称是Interfacial Transition Zone也就是界面过渡区。混凝土搅拌后骨料表面会形成一层水灰比偏高的水膜水化产物在这里的堆积密度低于砂浆本体导致ITZ孔隙率较高、微裂缝较多。宏观上它的强度和弹性模量都比砂浆低是混凝土里的“薄弱环节”。试验数据显示ITZ的弹性模量通常只有砂浆的50%-70%拉伸强度更低很多宏观裂缝就是从ITZ区域开始萌生的。细观模型中单独建出ITZ核心意义就是让模型具备描述这种弱化行为的能力。如果不单独建ITZ骨料和砂浆之间直接绑定或接触算出来的界面强度会偏高裂缝更倾向于穿过砂浆而不是沿界面发展这与实际破坏模式不符。特别是做混凝土拉伸、劈裂、断裂模拟时ITZ的几何位置和力学参数直接决定了裂缝是“绕骨料走”还是“穿骨料走”这是整个细观模型成败的关键之一。3.2 ITZ几何建模的三种常用方案ITZ在真实混凝土中的厚度只有20到50微米直接按这个尺寸建模网格规模将不可接受所以细观模型中通常采用“等效厚度放大”的处理方式。常见方案有三种我简单对比一下。第一种是薄层实体单元法。在骨料表面沿法向外扩一个等效厚度壳层壳层与骨料、砂浆之间共节点连接壳层赋予ITZ材料参数。这个方法实现直观能在模型里看到实实在在的ITZ区域也能输出该区域的应力应变。缺点是薄层区域网格容易畸变厚度方向至少要保证两层单元否则容易算歪。第二种是零厚度单元/接触内聚力法。骨料和砂浆之间不建实体层直接插入一层零厚度内聚力单元用双线性牵引-分离本构模拟界面开裂。这个方法更灵活不需要放大厚度天然适合模拟ITZ开裂过程。缺点是内聚力参数标定麻烦而且ANSYS里这类接触型内聚力单元的收敛性需要调很多参数。第三种是等效均质法。把ITZ的弱化效应折算到骨料表面一层单元或砂浆单元上比如给骨料表面单元赋降低后的模量。这个方法省事但严格来说不算真正的三相模型我做损伤分析时基本不用。3.3 本次模型的ITZ处理方式与理由我这次做的模型采用第一种薄层实体单元法同时把ITZ厚度设置为0.5mm。真实ITZ只有几十微米但在细观模拟中ITZ的作用是“在正确的位置提供一个弱化层”而不是精确复制它的几何厚度。把厚度控制在骨料最小粒径的10%以内不会对骨料形态和堆积密度造成太大干扰同时又能保证计算可行性。在ANSYS里实现时有一个小技巧先用APDL的VSBV命令把骨料体从砂浆母体中“减掉”让骨料和砂浆共用边界但互不交叉然后在骨料外表面用VOFFST命令向外偏移一个面得到ITZ壳层几何再把这个壳层与砂浆体做布尔运算。这样做的好处是三个区域的几何严格贴合没有缝隙后期网格划分不会产生意外的孤立面。4. ANSYS中完整建模实操流程4.1 建模方案对比APDL脚本还是Workbench这个模型我建议以APDL为主线。原因很现实细观模型的核心是大量重复的几何创建与材料赋值在Workbench里用鼠标点几百次不现实APDL则可以用循环批量完成。Workbench在这个场景里的定位是“检查几何和看结果”。具体流程是Python生成的多面体几何以IGES格式导出APDL里用VGLUE、VSBV做布尔装配材料赋值用MAT命令加循环网格用ESIZE控制尺寸后VMESH批量划分求解和边界条件全部写入APDL脚本。这样从几何到求解一轮计算可以很快跑通对参数化研究极其重要。顺带提一句ANSYS的安装和许可配置在这些自动化流程里是基础前提建议把工作目录和内存工作空间提前设好避免批量跑脚本时因为路径或内存问题中断。4.2 几何创建与装配的关键细节用IGES导入几十个多面体后第一步是检查几何完整性。Python生成的多面体面片和ANSYS的实体建模体系之间存在容差差异容易出现小裂缝或重叠面。我的经验是Python导出时把所有面都用闭合三角片构建导入后先用VGLUE把所有体粘接起来再用VSBV把骨料区域从砂浆中减去。这两个布尔操作是装配的核心前提是所有体都是“干净的”——没有退化面、没有零面积面、没有重复体。如果布尔操作报错多数情况是几何容差问题。ANSYS里可以用BTOL命令调高布尔容差但注意不要设太大否则会把本不该合并的边界合并掉。我一般先设BTOL1e-5不行再微调往上走。另外导入几何后要养成习惯用VLIST、VPLOT逐个查看体编号和位置确认没有缺失或错位。4.3 材料参数设置与单元选择三个区域的材料参数按典型值设置。骨料弹性模量50GPa、泊松比0.2、密度2700kg/m³砂浆弹性模量25GPa、泊松比0.2、密度2200kg/m³ITZ弹性模量取砂浆的60%即15GPa、泊松比0.22。如果后续要扩展成非线性损伤模型再给每个区域配塑性参数或内聚力参数本次按线弹性展开。单元类型上骨料和砂浆都用SOLID185ITZ因为几何很薄也用SOLID185但厚度方向至少两层单元。这里有个重要细节SOLID185默认全积分薄层单元在高纵横比下容易锁死建议打开增强应变选项也就是设置KEYOPT(2)2或者换用SOLID186高阶单元。我实测SOLID185加增强应变在ITZ区域表现稳定求解效率也高。4.4 网格划分要点与质量控制网格划分是整个流程中最容易翻车的地方。骨料和砂浆是不规则多面体直接自由网格划分会得到大量四面体能用但精度和收敛性不如六面体。由于多面体形状不规则全六面体网格几乎不可能折中方案是骨料内部用四面体砂浆区域用四面体ITZ薄层用扫描网格从骨料表面向砂浆方向生成至少两层。ITZ薄层网格要特别注意厚度太薄时扫描网格容易产生高纵横比单元导致求解矩阵病态。我控制ITZ厚度方向网格尺寸为厚度的0.5倍也就是至少两层面内尺寸控制在骨料面尺寸的三分之一左右。这样网格数量会大一些但单元质量有保障。划完网格后用CHECK命令检查负雅可比把畸形单元修正掉再进求解。4.5 边界条件、载荷与求解设置边界条件要跟试验工况对应。做单轴压缩模拟就在模型底面固定Z向位移顶面施加位移载荷做单轴拉伸底面固定顶面施加指定位移。关键点是细观模型不像宏观模型那样适合直接在模型表面施加载荷因为载荷直接加在骨料或砂浆单元上会造成局部应力集中。正确做法是在顶面和底面各加一块刚性垫板可以用SOLID185建一个弹性模量极高的区域载荷加在垫板上应力分布更均匀也更符合试验条件。求解设置方面先做静力分析打开大变形效应NLGEOM,ON打开自动时间步初始步长设为总时长的十分之一。非线性计算时打开线性搜索LNSRCH,ON和预测器PRED,ON能显著提高收敛概率。如果还是不收敛重点检查接触状态和网格畸变逐步细化。5. 常见问题与排查技巧实录5.1 骨料投放卡死或重叠率过高投放程序卡死常见原因有两个一是投放区域太小目标体积分数定得太高后续骨料找不到位置二是碰撞检测算法有遗漏SAT检测里候选轴构建不完整导致明明重叠的骨料被判为不重叠。解决思路给目标体积分数留余量从理论值往下调3%-5%投放顺序严格按粒径从大到小每次随机尝试时适当放宽旋转自由度检测时先用包围球快速排除再精确检测。我自己的经验是粒径跨度越大、骨料数量越多体积分数上限越低。设计模型时要根据计算资源和目标体积分数动态调整投放次数和重试次数别一口气追求60%以上的体积分数很容易卡死。5.2 ITZ薄层单元网格畸变薄层单元畸变是细观模型最容易出问题的地方典型表现是计算到一半突然发散或者结果出现异常尖峰。原因往往不是材料参数而是ITZ单元纵横比太大。比如厚度方向只有一层单元面内尺寸又大这种单元在荷载下很容易产生沙漏或锁死。我的处理流程第一步检查单元质量把纵横比大于10的单元挑出来看集中在哪些区域第二步如果畸形单元过于集中回到建模阶段把ITZ厚度加厚或调整网格尺寸第三步打开SOLID185的增强应变选项能极大缓解锁死问题。如果做了这三步仍发散那就检查ITZ与骨料、砂浆之间是否真正共节点有时候微小的几何间隙会在求解时造成局部应变异常。5.3 模型规模过大与内存不足细观模型的节点和单元数量很容易突破百万级。我记得最早做一批骨料600颗的模型时网格划分完就直接把机器卡死了。后来的经验是骨料数量控制在模型可承载范围内一般3D细观模型控制在100到300颗比较合理网格尺寸不要盲目加密ITZ区域加密骨料内部和砂浆远处用较粗网格ANSYS里通过设置WORKSPACE限制内存增长同时降低结果输出频率只保留关键载荷步的结果。还有一个实用技巧把模型切成四分之一或八分之一做对称简化。很多压缩和拉伸工况都有对称性能直接把计算规模缩小数倍而且对结果精度影响很小。对称面上加相应的对称边界条件即可。5.4 计算不收敛的排查思路细观模型不收敛时先不要急着改材料参数。第一步查网格把网格畸变单元查一遍用CHECK命令看有没有负雅可比第二步查接触或绑定设置骨料-ITZ-砂浆之间如果不是共节点而是用了接触对要检查接触状态是否稳定第三步逐步缩小时间步配合自动时间步和线性搜索通常就能收敛。如果特定载荷步反复不收敛我会回去检查几何装配骨料之间是不是出现了意外的小间隙布尔运算后是不是残留了极微小单元。这些问题在几何阶段就消除比后期调求解参数省事得多。记住一个原则细观模型的收敛性80%取决于网格质量剩下的才是求解参数的事。我个人在这套流程里反复跑了很多轮最深的体会是细观建模的瓶颈不在软件操作而在几何生成和网格质量这两个前置环节。前期把骨料生成和投放算法写扎实了ANSYS里的装配、网格和求解反而都是顺水推舟的事。另外ITZ的厚度和材料参数对结果非常敏感做参数敏感性分析时最好先把基准工况定下来再逐项扫参数。最后再分享一个小技巧机器内存有限时骨料数量超过200颗尽量用SOLID185低阶单元SOLID186高阶单元在细观模型里有时会让计算时间翻好几倍精度提升却非常有限。