搞混凝土抗冲击仿真的同仁多多少少都跟LS-DYNA打过交道宏观均质模型算整板、算梁柱都很方便可一旦要回答“裂缝从哪儿绕过去、为什么偏偏走某个界面”这类细观问题模型里如果没有真实的骨料和界面很多机理根本说不清楚。我最近在Ansys LS-DYNA环境里搭了一套三维混凝土细观骨料模型三组级配可调、界面过渡区ITZ精细化壳层处理骨料体积分数按目标含量反算修正整体可以像参数化模板一样反复使用。这篇把建模思路、关键命令流程和踩过的坑都整理出来适合已经有ANSYS APDL基础、想往细观损伤方向深入的工程师也适合做混凝土冲击、爆破、劈裂仿真的研究生参考。1. 项目拆解与三相体系设计1.1 细观模型为什么必须区分三相材料混凝土在宏观有限元里通常被当成均质材料方便当然方便但有个老问题始终解决不了裂缝到底沿着哪里扩展。实际切开试件或者看CT扫描混凝土内部粗骨料随机分布骨料周围有一圈疏松多孔的界面过渡区也就是ITZ。裂纹进入这个区域后极易偏转绕开强度高的骨料沿界面和砂浆走。把模型拆成骨料、砂浆、ITZ三相以后这类破坏路径才可能被模拟出来。我刚开始做细观模型时也怀疑过直接用带随机强度的均质单元不也能模拟出裂缝形态吗试过几轮就发现不行。均质模型的裂纹基本上是沿着最大主应力方向单刀直入而真实混凝土裂纹会绕着骨料走形成明显的迂回路径。只有把骨料强度、砂浆强度、ITZ强度分开赋予单元失效顺序才会因为局部应力集中而改变宏观上才能看到锯齿状、分叉状的裂缝形态。所以三相体系不是可选项是细观建模的基本盘。1.2 三组级配、体积分数与ITZ的耦合关系细观模型不是一个单纯的几何拼图。骨料粒径分布决定整体级配级配又直接影响堆积密度和孔隙率而ITZ的体积占比又跟骨料总表面积挂钩。颗粒越小越细总表面积越大ITZ占比越高模型整体就越“软”。所以在设计参数时三组级配的粒径范围、目标骨料体积分数、ITZ厚度必须放在一起权衡不能分开拍脑袋。有一个很容易踩的坑如果ITZ厚度取得过大而级配下限又卡得比较小细骨料周围的两层界面会互相搭接甚至把砂浆区域割裂成孤岛。结果就是模型里砂浆体积被严重压缩强度整体拉低。我一般把ITZ厚度控制在0.5到1毫米再根据级配下限做一次快速校核只有保证任意两个相邻骨料表面的ITZ外缘距离大于0模型才算几何自洽。体积分数方面骨料体积分数一般按试件总体积的40%到55%来定三组级配分别承担对应的体积比例投放完成后要回头统计实际体积分数误差控制在1%以内。2. 参数化设计与核心参数选型2.1 三组级配怎么划分最实用“三组级配可调”里的三组我习惯按粗、中、细三段划分。比如做标准混凝土试件粗骨料16到20毫米中骨料10到16毫米细骨料5到10毫米。如果你做的是小尺寸试件或者砂浆试块粒径要整体缩放但分组逻辑不变。级配的分布规律可以按富勒级配曲线来取核心公式是累计通过率P(D)与最大粒径Dmax的平方根关系P(D)等于(D/Dmax)的0.5次方。用这个曲线反推各组粒径范围内的骨料数量占比再结合体积分数换算成每一组到底要投放多少个球。实际操作里我一般先用Python或者MATLAB把坐标和半径清单生成好再输出成文本文件给APDL读入。为什么不在APDL内部写完整随机投放逻辑因为APDL的数组操作和循环效率不高几百个球做重叠检测时速度很慢而且一旦随机种子没控制好后处理阶段想重现结果很麻烦。外部脚本生成坐标一是可以精确控制随机种子二是方便对级配分布做统计分析算完体积分数再投放出错率低很多。每一组骨料数量不是拍脑袋定的。我先把试件总体积乘以目标骨料体积分数得到骨料总体积再按富勒曲线算出三组粒径段的体积占比通常粗骨料约占一半中骨料和细骨料各占四分之一左右最后用球的体积公式逐组反算数量。这里有个细节随机粒径在每组范围内按均匀分布取但为了逼近富勒连续级配每组的下限不要卡太死可以留一点重叠模拟自然筛分时的过渡状态。2.2 ITZ厚度与精细化网格尺度的权衡ITZ的“精细化”到底细到什么程度取决于你的计算诉求。如果只是为了在宏观破坏形态上看到界面脱粘ITZ做一层薄实体单元就够厚度取0.5到1毫米如果研究的内容是ITZ内微裂纹起裂、扩展的细观机制厚度可能要更薄甚至引入内聚力单元。但越薄的几何体在网格划分阶段越容易出问题单元尺寸一旦大于ITZ厚度这一层要么被网格吞掉要么畸变严重算出来全是负体积。我常用的做法是让ITZ厚度与最小网格尺寸之间保持2到3倍的关系。比如ITZ取0.6毫米那ITZ内部至少划分两层单元单元尺寸控制在0.2到0.3毫米。骨料和砂浆区域的单元可以适当放大到1毫米左右通过过渡网格连接。这样既保证了ITZ区域的应力梯度能被捕捉到又不至于让整体单元数量爆炸。顺便说一句ITZ材料参数很重要弹性模量通常取砂浆的50%到70%抗拉强度比砂浆低20%到40%密度可以按相同值处理主要影响的是局部刚度分配和损伤起始位置。2.3 体积分数的反向校验体积分数不是投完球就算结束。球体之间相互重叠会被判定无效而跳过实际投放进去的骨料体积大概率小于目标值模型整体可能偏软。解决办法是在投放循环里加一个实时累加器每成功投放一个球就把球体积加进去投放结束后计算实际体积分数再看与目标值的偏差。偏差超过1%时不要盲目增加粒径上限而是先看看是哪一组贡献不足。多数情况下是因为大粒径骨料投放困难可容纳位置少重叠率高。这时可以把大粒径组的数量先放足再补小粒径因为小颗粒可以钻到缝隙里空间利用率高。另外边界位置也要注意球体如果有一部分伸出试件边界在布尔减运算时会整体损失实际留在基体内体积变小。我习惯在做体积分数校验时直接按被截断后的体积计算而不是简单累加完整球体积这样和后续有限元模型里的真实体积更接近。3. 基于APDL命令的建模实现全流程3.1 坐标文件准备与数组VREAD读入APDL命令流的起点是把外部生成的骨料坐标文件读进来。下面是一段常用的VREAD读取方式! 假设坐标文件为agg_coord.txt每行编号, x, y, z, 半径 *DIM, AGG, ARRAY, 500, 5 *VREAD, AGG(1,1), agg_coord, txt, , 5, 500, 5 (F8.0, 4F12.6)读完数组后用DO循环遍历骨料。这里的5列分别存编号、x、y、z、半径。用APDL的好处是后续建模命令和参数设置都在同一个环境里不需要来回切换几个软件。文件路径要写成绝对路径或者放在当前工作目录否则VREAD容易读不出数据这个错误很隐蔽因为我第一次用的时候卡了快半小时。坐标文件里的数据预处理也很关键。生成骨料球时我建议先对所有骨料按半径从大到小排序再写入文件。大颗粒先投放能占据更合理的空间小颗粒后投放去填缝整体体积分数更容易达标。这一步在外部脚本里做APDL里只负责按顺序读入建模。3.2 球骨料、ITZ壳层与砂浆基体布尔运算在APDL里生成球骨料的命令组合我们可以用SPHERE结合坐标系抓取的思路。APDL的SPHERE命令并不像CAD软件那么直观它需要给六个参数分别是x1、x2、y1、y2、z1、z2相对当前坐标系的值。我通常先创建一个局部坐标系原点放在骨料中心再在这个局部坐标系里用SPHERE生成长方体包围球然后通过VSBV等布尔运算切出球体。不过说实话用局部坐标系逐个建球体几百个球的时候命令流会非常冗长。更高效的做法是先建一个“中心在原点、半径为1”的球体再用VGEN或平移命令移动到目标位置。代码大致是这样! 先建立单位球 SPHERE, -1, 1, -1, 1, -1, 1 ! 循环读取数组并平移 *DO, I, 1, N_AGG, 1 LOCAL, 11, 0, AGG(I,2), AGG(I,3), AGG(I,4), , , , , , VGEN, , 1, , , , , AGG(I,2), AGG(I,3), AGG(I,4), , , , 0 ! 缩放半径如果半径不为1 ! 这里可以配合 VLSCALE 对体做半径缩放 *ENDDO用VGEN平移体和用VLSCALE缩放的组合对圆形骨料来说逻辑非常直接。但注意VLSCALE是按体自身的坐标系原点缩放的所以要先平移后缩放或者先局部坐标建球再整体处理顺序不能反。缩放后需要再次把体放到正确位置建议直接在外部脚本里生成半径一致的球心坐标或者干脆把统一球的几何体进行缩放后重新定位。ITZ壳层建模我一般用“两球相减”的办法。核心骨料半径R加上ITZ厚度t形成半径Rt的外球。先建外球再减去内球剩下的空心球壳就是ITZ实体。代码如下! 假设已经用半径Rt建好外球体编号为V1用半径R建好内球编号为V2 VSEL, S, VOLU, , V1 VSEL, A, VOLU, , V2 VSBV, V1, V2 ! 从外球中减去内球最后处理砂浆基体。基体是长方体把所有骨料球和ITZ壳都从长方体里减掉剩下的区域就是砂浆实体。这里注意布尔运算顺序里一定要先减ITZ壳层再减骨料球否则ITZ和骨料之间的共享面可能不完整。减完之后用VGLUE把三相区域粘在一起保证后续网格划分时节点连续、单元共享这样在LS-DYNA里就不需要额外设置骨料与砂浆之间的接触省了很多麻烦。布尔运算在骨料数量多的时候非常耗时尤其是几百个球在同一个基体里反复做布尔减可能出现无法求解或者是体退化问题。我的经验是分批次处理一次读入50个左右骨料做完这一批的布尔减再继续下一批。如果需要中间停一下可以用SAVE保存模型状态调试时很方便。3.3 网格划分策略与LS-DYNA导出配置三相几何体划分网格时单元类型必须在进入求解器之前设定好。做LS-DYNA显式分析模型最终要转成SOLID164单元。在APDL前处理阶段可以直接建立SOLID164单元类型然后和普通实体单元一样划分网格。不过SOLID164是八节点六面体单元网格划分会受几何形状影响较大。对于球形骨料六面体网格很难保证质量我更推荐先用SOLID185划分出高质量的网格再在后处理和导出阶段转换成LS-DYNA可用格式或者使用ANSYS的LS-DYNA求解器输出时自动转换。实际流程是在前处理器中把ET设置为SOLID164划分网格时其节点数和拓扑和六面体一致但APDL对不同区域的网格算法会有自适应调整中心部分可能退化成四面体。如果对单元质量敏感可以先做四面体网格再通过MOPT转换但会牺牲一点求解精度。我的建议是骨料和ITZ用六面体映射网格砂浆基体用自由网格在布尔运算生成的交界面上插入过渡层保证界面位置单元对齐。网格尺寸的全局控制通过LESIZE、KESIZE或者ESIZE设置。骨料表面是曲率变化最大的位置加密程度要高于内部砂浆。实际操作中我先把骨料表面的单元尺寸定为0.5毫米ITZ层径向细分三层砂浆区域用1毫米到1.5毫米总单元数量控制在80万以内这样在单机上进行显式冲击模拟时计算时间还可以接受。完成网格后进入LS-DYNA输出配置。如果是在ANSYS Workbench的LS-DYNA环境中直接切换分析类型并生成输入文件即可。经典APDL环境里可以通过EDWRITE命令把模型写成LS-DYNA的k文件EDWRITE, LSDYNA, concrete_meso, k写文件之前要把每个体的材料号、单元类型、Part编号定义清楚。材料参数如果还没赋导出后k文件里只有几何和网格材料卡片必须在LS-DYNA的keyword文件里补充或者在前处理中定义好。我见过不少人导出k文件后直接送求解器结果报错Pid没有材料就是因为漏了材料定义。4. 材料参数映射与关键字文件配置4.1 三相材料本构怎么选LS-DYNA里材料卡片非常多细观混凝土建模没有必要每个相都选复杂的本构关键是强度和失效行为要区分开。骨料是脆性材料用MAT_BRITTLE_DAMAGE或者MAT_ELASTIC加失效都可以弹性模量高、抗拉强度高模拟压碎和穿透时能看到骨料本身的断裂。砂浆是基体相用MAT_CSCM_CONCRETE或者MAT_PLASTIC_KINEMATIC都能得到不错的拉压不等强效果损伤软化也能捕捉到。ITZ是薄弱面一般用*MAT_PLASTIC_KINEMATIC把屈服强度和失效应变调低或者用带粘结的接触面模型模拟脱粘。我的经验是ITZ不要单独用非常复杂的材料模型舍简取繁容易出问题。它的真实厚度在毫米量级单元本就不多用简单弹塑性加失效应变配合较低的失效应变参数就能复现“最先开裂、沿着界面扩展”的宏观现象。真正决定裂纹走向的是三相材料强度差值和失效应变的搭配而不是本构模型的复杂度。4.2 Part定义与接触策略在LS-DYNA里每个Part对应一种几何体材料组合。布尔运算后用VGLUE生成的是共享节点模型骨料、砂浆、ITZ之间不需要定义接触天然就是连续的传递应力不会产生间隙。这种建模方式适合弹性阶段到损伤初期如果模拟到骨料和砂浆大面积脱粘、甚至颗粒飞溅建议还是把三相区域分开网格用接触来定义界面比如*CONTACT_TIED_SURFACE_TO_SURFACE_FAILURE允许界面在达到强度之后失效脱开。接触做失效的时候要注意失效准则参数的标定。我不是直接把抗拉强度填进去而是先算单元应力看破坏发生的位置是否集中在ITZ附近。如果ITZ单元全部先失效再逐步扩展到砂浆说明参数合理如果一开始就在骨料内部失效可能是骨料的强度设得太低或者是接触刚度出了问题。边界条件方面LS-DYNA里常用关键字是BOUNDARY_PRESCRIBED_MOTION_RIGID或者LOAD_BODY_Z。如果是做冲击模拟可以把试件底部固定或者给冲击体施加初速度。对细观模型来说边界效应是绕不开的靠近边界的骨料数量偏少、基体体积偏大如果不关心边界附近破坏建模时建议在目标区域旁边垫一层同材料保护区减少试件边界对核心区域应力场的影响。5. 实操过程中踩过的坑与排查方法5.1 骨料随机投放卡死与重叠问题随机投放骨料最烦人的问题就是重叠。几百个球放下去后放的球只要和前面任何一个球距离小于半径之和就必须重新生成位置。最简单直接的方法是每放一个球都遍历所有已放球的距离一旦重叠就重新抽样。但规模一大这个计算量是O(N^2)APDL循环吃不消。我的方案是引入空间网格哈希把试件划分成若干小立方体网格每个球只跟所在格子以及相邻格子的球做重叠检测计算量可以从平方级降到接近线性。重叠判据不能只考虑球心距还要留出最小净距。因为骨料膨胀、网格过渡和布尔运算的容差需要额外空间如果两个球刚好相切后续布尔运算经常产生退化面网格划分就会失败。我留的净距一般不小于0.2毫米球心距大于等于两半径和再加0.2。另外在圆柱形或长方体边界附近球心离边界太近也容易导致模型不对称我通常把边界净距设为直径的0.1倍防止骨料被边界切得只剩一点点。5.2 布尔运算失败与单元退化布尔运算是APDL建模里最大的不稳定因素。球壳减核心球、还有基体减所有球壳每一步都依赖几何容差判断。骨料数量多尤其两个球之间净距很小时减运算容易出现重合面或者生成非常薄的片体最后网格划分时产生大量畸形单元。遇到特殊情况我的做法是把异常骨料的半径稍微调小0.01毫米重新执行布尔减大多数情况下就能正常出结果。如果仍然失败就单独删掉这个骨料后续用一组更小的骨料补充体积分数。网格划分阶段常见的负体积报错多半是ITZ壳层单元高度方向只有一层而厚度方向应力梯度又很陡。这时候不要急着缩小全局单元尺寸那样会让单元数量爆炸更合理的做法是单独对ITZ层的线设置分段数强制径向划分三层以上与此同时骨料内部仍然保持稀疏网格。5.3 计算成本控制与质量缩放细观模型单元数量大显式分析的时间步长受最小单元尺寸控制。如果裂纹扩展区集中在骨料界面附近而那边单元尺寸只有0.3毫米时间步长会被拉得非常小计算可能几周都算不完。常用的缓解手段是质量缩放LS-DYNA里通过*CONTROL_TIMESTEP的参数DT2MS来统一限制时间步长把过小单元对应的质量增加从而放大稳定时间步长。这个参数不是随便设的通常让增加的质量占比控制在5%以内否则算出来的惯性效应会扭曲冲击响应。我最终跑完一遍的模型总单元数约70万质量缩放后时间步长是原始步长的两倍多整体计算时间从原本预计的10天降低到4天得到的裂纹形貌和未缩放模型差别很小。这里提醒一下质量缩放对静力或者准静态分析影响非常大但混凝土冲击模拟本身就是高速动态过程只要增重比例控制严格仍然在可接受范围。5.4 体积分数反算与级配修正投放完成后统计模型里骨料单元体积与总体积的比值经常发现比目标值低2%到3%。原因有几方面随机投放产生空洞边界截断导致部分骨料体积丢失布尔运算从基体里减去球壳时如果球壳厚度接近单元尺寸一部分体积在网格划分阶段被膨胀或者压缩掉。我的修正方法是给每一组骨料设置一个体积补偿系数在下一轮循环里把每一组的投放数量乘以这个系数再重新生成坐标。一般两轮迭代之后体积分数就能稳定在目标值正负0.5%的范围内。还有个容易忽略的点三组级配对体积分数的影响并不均匀。粗骨料单个体积大数量少随机性波动很大细骨料数量多体积占比相对稳定。所以在做级配修正时要先调粗骨料组再微调中骨料和细骨料不然反复调三轮都不收敛。6. 这套模型的扩展空间模型跑通之后可以做的事情其实很多。比如把球骨料替换成随机凸多面体骨料更接近真实碎石形状破坏路径会有些差别再比如把ITZ从实体壳换成内聚力单元模拟界面脱粘时会更有物理意义或者把三组级配增加成四组五组配合更精细的连续级配曲线可以复现更多真实配比混凝土的力学响应。我自己的体会是细观模型的价值不在于算得特别豪华而在于能把三两个关键问题回答清楚。建模命令和参数体系一旦搭好后续做参数敏感性分析、级配对比、ITZ强度影响分析都只是换一套数据重新跑一遍的问题。如果你们也在做类似的方向强烈建议把整个建模流程脚本化尽量不依赖GUI手工点选这样调一组粒径、改一个ITZ厚度代价只是重新生成一次坐标文件而不是重新画一遍模型。