做混凝土数值模拟的朋友应该都有同感：宏观模型算起来快，但很多宏观现象解释不了。裂缝从哪来、怎么穿过骨料还是沿界面走、砂浆和骨料的协同变形怎么分配——这些问题都得往细观尺度去看。我一直在做这方面的研究，最近完成了一个含界面过渡区的多面体骨料密堆积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命令调高布尔容差，但注意不要设太大，否则会把本不该合并的边界合并掉。我一般先设BTOL=1e-5，不行再微调往上走。另外，导入几何后要养成习惯用VLIST、VPLOT逐个查看体编号和位置，确认没有缺失或错位。
4.3 材料参数设置与单元选择
三个区域的材料参数按典型值设置。骨料：弹性模量50GPa、泊松比0.2、密度2700kg/m³；砂浆：弹性模量25GPa、泊松比0.2、密度2200kg/m³；ITZ：弹性模量取砂浆的60%，即15GPa、泊松比0.22。如果后续要扩展成非线性损伤模型，再给每个区域配塑性参数或内聚力参数，本次按线弹性展开。
单元类型上，骨料和砂浆都用SOLID185，ITZ因为几何很薄，也用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高阶单元在细观模型里有时会让计算时间翻好几倍，精度提升却非常有限。
