做混凝土细观数值仿真的人应该都有这种体会:模型图怎么画都好看,真正要把它建出来,骨料的几何尺寸、空间位置、体积分数、界面层厚度,每一项都是在给自己挖坑。我最早在ANSYS LS-DYNA里尝试三维混凝土细观骨料模型的时候,就是很老实的一个人工建模路线,球体一个个复制平移,尺寸靠眼睛调,放进去之后连“骨料体积到底占了多少”都只能估个大概。后来被审稿人一句话问住——“ITZ层怎么处理的,网格有没有畸变?”我才意识到,这个建模过程必须参数化,必须能用命令流重复生成,才能称得上是一个可用于研究的建模方案。
这篇文章就围绕这个需求展开,目标很明确:在LS-DYNA对应的关键字体系下,把三维细观混凝土模型的“三组级配可调、ITZ精细化处理、骨料体积分数自动计算”这三件事讲透。适合正在做混凝土冲击破坏、爆炸响应、局部剥落等显式动力学仿真的研究生和工程师,尤其是已经对LS-PrePost和关键字文件有一定基础,但还没找到一套可靠细观建模流程的人。
1. 细观模型到底在解决什么问题
1.1 宏观模型和细观模型的本质差异
LS-DYNA里常见的混凝土模型,比如MAT_072R3(CSCM)、MAT_084/085(Winfrith)、MAT_111等,都是宏观均质模型。它们把混凝土看作一个连续均匀材料,用一套本构和损伤参数来描写“平均意义下”的力学响应。计算效率非常高,单元尺寸到了10毫米量级可以用,工程整体响应和破口尺寸大致能对得上。
但一旦关心的是裂纹路径、骨料从砂浆中拔出、骨料自身断裂这种局部破坏机制,宏观假设就跟实验现象对不上了。真实混凝土的破坏路径几乎全部发生在骨料与砂浆之间,也就是ITZ区域附近。裂纹遇到骨料会绕过去,遇到软弱界面会沿着界面扩展,只有应力水平足够高的时候才出现骨料劈裂。宏观模型永远无法回答“裂纹为什么会走这条路”的问题,因为细观结构信息根本没有进到模型里。
细观模型就是为了把这种局部机制还原出来的。把混凝土拆成三类或四类组成部分:粗骨料、砂浆基体、ITZ层,以及可能存在的初始孔隙。骨料用弹性或高强弹塑材料描述,砂浆用带软化的率相关材料描述,ITZ则用更弱的性能参数或专门的粘聚失效模型去模拟。然后让裂纹自己选择路径,而不是人为指定破坏面。
1.2 LS-DYNA里建细观模型的特殊性
很多人在细观建模时第一反应是去用ABAQUS或者COMSOL。这两款软件做二维细观模型很方便,但放到LS-DYNA里,事情立刻变得不一样:LS-DYNA是显式求解器,对网格质量、最小单元尺寸、沙漏控制都极其敏感。一套细观模型往往有几十万甚至上百万单元,骨料与骨料之间还夹杂着ITZ薄层,单元尺寸差距可能达到一到两个数量级。
这意味着,在LS-DYNA里建细观模型,不能只考虑“几何上放得下”就行,还必须考虑“求解器能不能扛得住”。ITZ层太薄,最薄单元会直接把你整个模型的时间步拖到微秒级以下,计算时长从几小时变成几周。骨料之间贴得太近,四面体网格的极瘦单元会产生负体积,求解中途直接终止。这些都不是画图阶段能看出来的,要在建模命令里提前留出约束。
所以,这篇文章里我讲到的“建模命令”,不只是一段生成骨料位置和节点坐标的脚本,更像是一套完整的建模约束逻辑:级配参数控制骨料尺寸,投放逻辑控制空间分布,ITZ厚度控制界面网格层的尺度,体积分数作为全局校验指标。每一环都跟LS-DYNA的求解稳定性挂钩。
2. 整体设计思路:三组级配、ITZ、体积分数怎么协同配合
2.1 三组级配可调的核心逻辑
所谓“三组级配”,对应的是工程混凝土粗骨料通常按照粒径范围分成三级或四级的做法。常见做法是把4.75毫米到26.5毫米范围内的骨料分成三档:
- 细档:4.75~9.5毫米
- 中档:9.5~16毫米
- 粗档:16~26.5毫米
每一档通过最小直径和最大直径两个参数定义。级配可调的意思,就是这两个直径边界,加上每一档的体积占比,全部做成参数。想模拟细骨料多的混凝土,就把细档占比调高;想模拟大石料密集的碾压混凝土,就把粗档占比拉起来。整个过程不需要重新画图,只需要修改参数后重新运行生成命令。
这里其实有一个比例尺问题需要想清楚。骨料直径区间和体积占比都必须对应同一个参考级配曲线。比较常用的是富勒级配或者简化富勒级配。富勒级配给出的累计通过率与粒径的平方根成正比,对应到各档骨料数量上,并不简单地是“大骨料体积占30%”。如果直接用体积比例分配样本空间,得到的中小骨料数量会少得离谱。
我自己更建议的做法是:目标总骨料体积分数是固定的,比如0.35或0.40,各档骨料按“体积贡献比例”分配到总骨料体积中去。也就是如果细档目标体积占总骨料的30%,那就用0.4×0.3=0.12作为细档在总模型体积中的目标体积,再根据每个骨料球的体积去倒推需要生成多少个骨料。这个换算过程是体积分数自动计算的基础,也是三组级配能真正“可调”的保证。
2.2 ITZ精细化处理的三种可选方案
ITZ,也就是界面过渡区,是骨料表面附近厚度只有几十微米的薄层区域。在真实混凝土中,这个区域水灰比偏高、氢氧化钙晶体富集、孔隙率高于砂浆基体,因此强度和刚度都比砂浆低,通常被认为是混凝土中最薄弱的环节。
在LS-DYNA的细观建模里,ITZ的处理方案大体上有三种。我用一张表列出来:
| 处理方案 | 建模方式 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|---|
| 独立薄层单元 | 每个骨料表面向外扩展一个壳层,单独划分网格并赋予弱化材料 | 直接模拟ITZ厚度,力学机理最清楚 | 网格量巨大,最小单元尺寸受ITZ厚度控制,计算步长严重下降 | 小尺寸试件、单骨料验证模型 |
| Tiebreak接触 | 不在骨料与砂浆之间建物理薄层,而是在界面接触定义拉伸/剪切软化失效 | 网格尺寸不因ITZ而受限,计算效率高 | ITZ厚度不具有几何物理意义,只能通过接触参数等效 | 大尺寸细观模型、动力破坏分析 |
| 等效属性折减 | 将骨料外层相邻一圈砂浆单元的强度、模量按ITZ比例折减 | 无需特殊接触,单元连续性最好 | ITZ厚度受网格尺寸影响,不够精细 | 轴对称或二维简化模型 |
从我自己的经验来看,纯显式大变形破坏分析里,方案一“独立薄层单元”是理想化程度最高、也最容易把求解器搞崩的方案。方案二“Tiebreak接触”和方案三“等效属性折减”更实用。标题里说的“ITZ精细化处理”,我建议理解成三层意思:厚度可设、材料强度可调、界面失效行为可控。不一定非要在几何上建出那几十微米的壳层。
2.3 体积分数自动计算在建模流程里的角色
体积分数不是建模结束后的统计报告,它是整个生成过程的约束条件。如果你投放到最后发现骨料体积分数只有0.25,而目标值是0.40,那说明投放算法拒绝率太高或者级配参数换算错误,模型不能带着这个误差交给LS-DYNA求解。
我这里说的“自动计算”,有两个层面。第一个层面是生成过程中实时累加:每成功投放一个骨料,就把当前骨料体积加到总骨料体积变量里面,投放完成后立刻输出实际体积分数。第二个层面是生成网格后,在LS-PrePost里按单元材料ID统计所有骨料单元的总体积,与模型总体积相除。两个结果彼此校验。
把体积分数自动计算做成流程的一部分,还有一个更大的好处:当三组级配的参数变化时,你可以快速检查不同配比下体积分数是否维持稳定。比如调整粗档最大粒径从26.5毫米改成31.5毫米,粗骨料数量会明显减少,如果体积分数保持不变,说明抽样逻辑是自洽的。
3. 实操:从参数设定到LS-DYNA关键字生成
3.1 Step 1:设定几何参数与级配参数
建模第一步,先把所有全局参数写到脚本头部。以最常见的100毫米立方体试件为例,我通常会这样设定:
模型边长: 100 mm 骨料总体积分数目标: 0.35 级配档数: 3 第1档: 4.75 ~ 9.5 mm,体积占比 35% 第2档: 9.5 ~ 16 mm,体积占比 45% 第3档: 16 ~ 26.5 mm,体积占比 20% ITZ厚度: 0.05 mm 骨料之间最小净距: 0.10 mm 随机种子: 20250620(可改)这里有三个细节需要注意。
第一,骨料总体积分数选择0.35左右是常见工程范围,高强混凝土和自密实混凝土可以到0.45,但显式模型里体积分数越高,骨料投放算法越难收敛,模型中的浆体润滑层就越薄。
第二,三档的体积占比加起来必须等于100%,每个档的骨料目标体积就是0.35乘以对应百分比,也就是0.1225、0.1575、0.07。这个换算关系在代码里要固定,否则容易出比例错乱。
第三,随机种子必须保留。细观模型的随机性很强,同一个级配参数下换一个种子,骨料位置完全不一样。如果后处理发现某组结果有异常,需要回到完全相同的模型,就必须锁住这个种子。
3.2 Step 2:骨料随机投放算法与碰撞控制
骨料投放是整个流程中最容易卡住的地方。核心思路是“先大后小”的拒绝-重试方法:先生成所有大骨料,然后依次放入中档和小档骨料;每放一个骨料,跟已放置的所有骨料做距离判断,如果重叠就重新随机生成位置并重新判断。
用伪代码来说,核心逻辑是这样的:
for size_group in groups: for i in range(num_particles[size_group]): attempt = 0 while attempt < max_attempts: x, y, z = random_position_in_cube() r = random_radius(d_min, d_max) / 2 if not overlap_with_existing(x, y, z, r): place_particle(x, y, z, r) total_volume += 4/3 * PI * r**3 break attempt += 1这里有个很多人会忽略的细节:碰撞判断用的是什么距离。如果ITZ层厚度是0.05毫米,骨料球半径是R1和R2的两个骨料,它们能允许的最小中心距离应该是:
distance >= R1 + R2 + ITZ厚度 x 2 + 最小净距ITZ厚度为什么要乘以2?因为每个骨料两侧都有自己的ITZ层,两个骨料各自外扩0.05毫米后,两者界面之间还需要留出至少0.10毫米的砂浆净距。如果是用独立薄层单元法建模,这个约束直接决定ITZ单元不会在骨料间交叉重叠;如果只用Tiebreak接触,这个约束也能保证砂浆单元有足够空间,不会生成极薄且畸变的四面体。
在实际实现中还有一个优化策略:不要每放一个骨料就扫描所有已经放置的骨料列表。大模型可能有几千个骨料,两两判断的计算量是O(n²)。可以先把立方体空间划分成若干子网格,碰撞检查只遍历当前骨料所在子网格及其相邻子网格内的骨料,速度能快出一个数量级。
3.3 Step 3:ITZ层的生成与材料赋予方式
如果选择了独立薄层单元路线,ITZ层的几何生成有一个比较省事的技巧:以骨料球心为球心,用半径r + ITZ_thickness再生成一个同一尺寸、同球心的“外扩球”。夹在原始骨料球表面和外扩球表面之间的壳状空间,就是ITZ层的几何区域。
但直接拿球壳去划分网格并不方便。LS-DYNA普通的前处理流程里没有一个命令能一键生成球壳单元。比较现实的做法是,在骨料球面和外扩球面之间人工布置若干层径向单元,层数不要太多,1到2层就够。这个操作可以在LS-PrePost里通过对单个骨料手动划分,或者通过脚本控制网格划分软件批量生成。
如果选择Tiebreak接触方案,建模流程里其实不生成ITZ几何,而是在骨料外表面和砂浆内表面之间定义接触:
*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE_TIEBREAK关键参数是法向拉伸失效强度FS和切向剪切失效强度FD。Tiebreak接触的作用是让两片网格在未失效前像一体一样共同变形,一旦界面应力超过设定强度,接触关系释放,骨料和砂浆就分离了。这种方法本质是用接触参数等效ITZ层的弱化行为,不算严格的厚度建模,但工程上已经很实用。
这里有一个明确建议:ITZ层千万不能做成一个网格内极小厚度。比如ITZ厚度设成0.02毫米,在20毫米骨料旁边,这套模型的最小单元尺寸就成了0.02毫米量级。LS-DYNA显式时间步长与最小单元尺寸成正比,这个薄层会让整个模型的时间步下降到可以忽略的程度,算力全浪费在一个细小到对宏观破坏贡献不大的壳层上。
3.4 Step 4:生成LS-DYNA关键字文件的核心要点
骨料位置算完,ITZ方案选定,接下来就是把这些几何和材料信息写成LS-DYNA能运行的.k文件。一个最简可用的细观模型关键字文件,至少需要以下几块:
*NODE *ELEMENT_SOLID *PART *SECTION_SOLID *MAT_ELASTIC *CONTROL_TIMESTEP *CONTROL_TERMINATION你的建模命令脚本要做的,就是把节点、单元、PART和材料卡自动写出来。骨料单元统一赋PART ID 1,砂浆赋PART ID 2,ITZ层如果有独立几何则赋PART ID 3。所有单元的类型、节点编号、材料ID必须一一对应。这里建议在生成脚本里做一个单元核对环节:每个单元的质心坐标落到哪个骨料球里面,就分配哪个PART。这个“归属校验”能有效避免网格划分阶段边界节点归属错乱。
材料卡部分,最简模型可以这样给:
*MAT_ELASTIC $ MID, RO, E, PR 1, 2.7e-9, 50000, 0.2这是标准的units单位制,密度、弹性模量、泊松比。注意LS-DYNA里没有强制单位系统,关键是全模型自洽,质量、长度、时间单位的换算不能混。
3.5 Step 5:体积分数自动统计与校验
这一步很好写。生成完所有骨料单元后,遍历整个模型的所有体单元,按材料ID分类累加体积。用LS-PrePost的*DATABASE_EXTENT或者后处理Model Check功能都可以做。但最可靠的是直接在生成脚本里自己算:每个体单元的几何体积积分求和。
我常用的校验方式是分两道:
第1道:几何级配校验——统计每个档位的球体数量是否与预期接近 第2道:网格级配校验——统计模型内骨料单元总体积 / 模型总体积如果两道校验的偏差都小于2%,模型就可以放心提交求解。如果偏差很大,优先检查骨料投放阶段被拒绝的数量。一般来说,最大重试次数达到50次还没有成功放入的骨料,不是随机种子的问题,而是“目标体积分数已经超过了该填充算法的上限”,需要调低骨料体积分数或者缩小骨料尺寸范围。
4. 材料模型与接触设置经验
4.1 骨料、砂浆、ITZ分别用什么材料
骨料在LS-DYNA里用一个弹性材料通常够了,因为骨料的破坏不是细观模型研究的主要对象,除非加载水平非常高。弹性材料参数按照花岗岩或石灰岩来给,密度2600~2800 kg/m³,弹性模量40~70 GPa,泊松比0.18~0.22。
砂浆基体则是破坏主体,建议用带拉压不对称和损伤退化能力的材料。MAT_072R3(CSCM)或者MAT_111都可以。但要警惕一个问题:LS-DYNA这些混凝土材料模型内部有自己的软化长度和能量释放参数,它们默认针对宏观单元尺寸标定。细观模型的砂浆单元尺寸往往只有1~2毫米,远小于宏观标定范围,材料参数必须重新标定或至少做一个单单元验证。
ITZ材料如果采用独立薄层单元,最简单的办法是给一个比砂浆低40%~60%的强度参数,弹性模量降20%~30%。这里我给一组参考值思路:砂浆抗拉强度4.5 MPa,ITZ抗拉强度给2.5 MPa左右;砂浆弹性模量26000 MPa,ITZ给18000 MPa。实际数值取决于你的试验数据和标定目标,但这一比例关系反映了ITZ的弱化本质。
4.2 网格共节点与接触分离怎么选
细观模型里,骨料和砂浆网格如果完全共节点,那么在拉伸破坏时,裂纹不能沿骨料表面滑移延伸,只能“穿透”界面,这就是最典型的强度高估。如果不共节点,界面之间就是接触对,接触算法会引入接触搜索、摩擦和滑移,物理上更接近骨料与砂浆之间的界面行为。
但这不代表“不共节点就更好”。LS-DYNA里所有不共节点界面都需要定义接触。界面数量一多,接触搜索的开销很大,而且接触初始穿透问题会频繁出现,尤其是骨料表面经过网格划分后曲率突变,会在局部产生微小穿透。设置*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE_TIEBREAK时必须预留IGNORE和MORTAR等参数组合来忽略初始微穿透。
我自己更推荐的计算路线是:砂浆和骨料网格在几何上完全连续,但给骨料外表面周围一层砂浆单元一个弱化材料赋值,模拟ITZ;如果需要完整裂纹分支,再把骨料与砂浆之间的接触改成Tiebreak型。这种方案网格质量最好,时间步最稳定,也不容易出现接触穿透导致的求解发散。
4.3 显式时间步、沙漏与质量缩放
细观模型单元小,必然导致计算时间成本高。LS-DYNA的时间步受最小单元尺寸控制,而细观模型最小单元通常出现在骨料与骨料之间的砂浆缝隙或者ITZ薄壳上面。如果你发现时间步小到离谱,优先要做的不是无限质量缩放,而是检查是不是有极端瘦长单元存在。
沙漏是另一个LS-DYNA老问题。单点积分六面体单元默认最容易出现沙漏模式,表现为单元变形呈锯齿状但应变能很低。细观模型里骨料大多是刚体或弹性材料,沙漏能被放大到整个模型。一个比较稳的配置是:
*CONTROL_HOURGLASS IHQ=4,QH=0.05同时给骨料单元用全积分单元公式(ELFORM=2),砂浆单元保留单点积分但开启增强沙漏控制。骨料全积分算得慢一点,但至少不会出现骨料被压成“眼睛”形状的诡异画面。
5. 常见问题与排查技巧实录
5.1 骨料投放率低、体积分数达不到目标
这个问题我在第一次做模型时栽得很惨。目标体积分数0.4,生成的骨料球怎么也放不进去,程序反复重试,最后实际体积分数只有0.26。排查后发现,我设置的最小骨料间距等于ITZ厚度的两倍,但骨料半径分布是从4.75毫米到26.5毫米,大量中档骨料彼此之间几乎没有空隙。
解决思路有三条:第一,把最大重试次数从20提高到100;第二,对每个档位单独设置随机种子,避免低效的伪随机聚集;第三,允许骨料半径在给定档位范围内做轻微缩径,缩径幅度不超过5%。这个缩径操作对体积分数影响很小,但对填充率提升明显。另外还可以先做一次“空投测试”,仅用几何碰撞检查统计该级配理论最大填充率,如果连理论值都低于目标,就真的需要调整级配比例了。
5.2 网格生成后出现负体积或极差单元
细观模型的网格质量是LS-DYNA显式计算的生死线。常见的负体积出现在两种位置:骨料之间间距太窄形成的压扁四面体,以及ITZ薄壳层与相邻骨料表面相切导致的重叠单元。
解决办法是给骨料间距增加一个下界。比如最小净距设置0.15毫米,同时保证ITZ厚度小于这个净距的一半,这样所有砂浆间隙里至少能塞进一个完整单元。如果已经生成了坏网格,不要直接在求解器里硬算,回LS-PrePost里对缺口区域做局部重划分,或者干脆丢弃该区域的一层网格,重新局部建模。
5.3 裂纹沿界面扩展不明显,模型破坏形态失真
这种情况多半是ITZ弱化没有真正进入计算。如果你用等效属性折减方案,但骨料表面第一层砂浆单元网格尺寸已经是1毫米,这个网格单元里的平均强度受到周围整体区域影响,实际削弱效果被模糊掉了。Tiebreak接触方案里,则要检查接触失效参数有没有真的在变形中被激活。
建议在后处理时单独输出骨料-砂浆界面处的单元失效应变或接触面法向应力。如果应力早已超过设定强度但裂纹没有打开,很可能是Tiebreak接触连接没生效,检查接触卡里的OPT参数是否设成了1(法向分离+切向滑移失效模式)。还有一点:扁平骨料表面在LS-DYNA中很难自然生成Tiebreak接触对,因为接触搜索默认是主从面法向,曲率大的表面识别并不好,尽量在骨料几何阶段就使用球体或缓变多面体。
5.4 求解中途终止:时间步无限变小和负体积
负体积在LS-DYNA错误信息里几乎成了细观模型的代名词。通常发生在高应变率冲击阶段,某个砂浆单元被压缩到翻转。解决办法按优先级排列:
- 把ITZ厚度稍微调大,或者把骨料之间缝隙放宽;
- 把砂浆单元公式从默认的单点积分改成2号全积分;
- 对目标区域做局部网格细化,而不是全模型统一细分;
- 最后才考虑适当质量缩放。质量缩放系数过大时,模型惯性效应会掩盖真实破坏行为。
我在实际项目中还发现一个规律:负体积大概率发生在粗骨料的尖锐角部附近。原因是尖锐骨料几何在网格划分时必然留下短边元素。所以我后来生成骨料几何时都会对骨料顶点做圆角化处理,牺牲一点点形状真实度,换取单元质量的大幅提升。
5.5 后处理裂纹云图不连续
LS-DYNA后处理里,细观模型的裂纹显示很容易出现“断断续续”的效果。这不是你的模型不对,而是损伤变量输出频率太低。默认的*DATABASE_BINARY_D3PLOT输出间隔可能在0.01毫秒级别,对于骨料断裂这种瞬间事件,帧与帧之间裂纹已经错位。把DT改小到0.002毫秒,同时增加*DATABASE_EXTENT_BINARY里的失效变量输出选项,裂纹演化就会连贯起来。
6. 这套流程还能怎么扩展
最后说一点个人经验。三组级配、ITZ、体积分数自动计算这套流程,熟练之后最大的价值不是“能建出模型”,而是“能批量生成对比算例”。比如研究骨料体积分数对混凝土抗侵彻性能的影响,直接改β从0.3到0.45,脚本跑四个模型出来,三个小时就能完成全部建模。看裂纹路径差异,看能量吸收曲线差异,看破片分布差异,这就是细观模型区别于宏观模型的真正输出。
我自己后续的扩展方向有两个。一是引入更多样的骨料几何,不只是球体,加入随机多面体和椭球体,这样骨料咬合作用更明显,骨料体积分数也能做得更接近真实混凝土结构。二是把ITZ方案从等效接触升级为连续级配界面模型,在骨料和砂浆之间插入一层可弯曲的粘聚单元,尝试在不同加载率下再现界面破坏的应变率敏感性。这两个方向都建立在现有参数化命令流之上,只改几何生成模块和材料卡,结构不需要大改。
上面这些内容,如果你也正在被细观混凝土模型的网格质量折磨,建议先不要急着追求更高精度的骨料形状,先把三组级配可调、ITZ处理、体积分数校验这条链跑通。模型在LS-DYNA里能稳定求解、能复现裂纹路径,远比你画一个漂亮的骨料球重要得多。