2. 多面体骨料与纤维混合从零搭建ABAQUS参数化插件搞混凝土细观模拟的朋友应该都有体会在ABAQUS里手动建立随机骨料模型简直就是一场灾难。每次想生成一批随机分布的多面体骨料和乱向纤维都要写一堆Python脚本调参调到怀疑人生换个骨料体积分数又得从头再来。我前前后后被这个问题折磨了快两个月最后咬咬牙写了一个参数化插件把多面体骨料生成、纤维随机分布、碰撞检测、嵌入约束全部封装成可视化界面现在点几下鼠标就能生成一个完整的细观模型。这篇文章就把插件的核心思路、源代码解析、以及调试过程中踩过的坑一次性说清楚希望对正在做混凝土、岩石或者复合材料细观模拟的朋友有价值。2.1 插件能做什么解决了什么问题这套插件解决的核心痛点是随机几何模型的可复现性和参数化控制。传统方式下每改一次骨料体积分数就得重新跑一遍随机撒点算法而且一旦生成结果不满意很难定向调整某一个参数。插件把所有几何参数都变成GUI里的输入框基体尺寸、骨料粒径范围、骨料体积分数、纤维体积分数、纤维长径比、随机种子改任何一个参数点一下就能重新生成模型。配合固定随机种子同一套参数可以生成完全一致的模型这在做正交试验、参数敏感性分析时特别重要。比如研究纤维掺量对强度的影响理论上应该保证骨料分布完全一致只改变纤维掺量这一个变量固定种子号就能实现。插件适合三类用户一是做混凝土细观力学研究的科研人员二是做纤维增强复合材料仿真分析的工程师三是对ABAQUS二次开发感兴趣的Python开发者。如果只是想快速得到一批随机骨料模型不需要了解底层算法直接用GUI就能满足需求。3. 插件整体架构与核心设计思路3.1 大框架从输入参数到CAE模型的数据流整个插件的架构非常清晰核心就是参数解析 → 几何生成 → 碰撞检测 → 装配嵌入四个环节。源文件包含三个模块param_utils.py负责参数校验和默认值管理geometries.py负责骨料和纤维的几何实体生成placement.py负责随机分布和碰撞检测逻辑。GUI层通过ABAQUS的RSG对话框构造器生成收集到的参数会封装成一个字典传到主程序。主程序的执行流程是读取GUI参数校验合法性根据基体尺寸创建基体部件一个长方形体循环生成骨料每次生成后做碰撞检测失败则重新生成循环生成纤维同样做碰撞检测全部几何体生成完毕后执行装配和布尔合并创建Embedded Region约束把骨料和纤维嵌入基体赋材料属性、划分网格、设置分析步关键设计决策是把几何生成和碰撞检测完全分离。geometries.py只负责生成几何数据纯粹返回坐标信息不关心模型是否自相交。placement.py只负责判断几何体之间是否冲突不关心几何体长什么样。这样的好处是以后想换成椭球骨料、矩形截面纤维只需要改geometries.py放置算法可以原封不动地复用。3.2 多面体骨料凸包算法生成近似真实形态真实骨料不是球形的而是不规则的凸多面体。插件用随机点凸包法近似生成在目标粒径范围内随机撒一组空间点求这些点的凸包得到的凸多面体就是一颗骨料。凸包算法我直接用了scipy.spatial.ConvexHull优点是准确可靠但有一个麻烦ABAQUS自带的Python环境不一定预装SciPy。如果目标机器没有SciPy插件会运行失败。稳妥的做法是自己实现一个简化版的凸包算法或者把SciPy和NumPy一起打包部署。我自己在写插件时选择打包部署因为NumPy在ABAQUS里基本都有但SciPy版本兼容性不好说打包能一劳永逸。凸包生成的具体流程是先在半径为R的球面上随机生成若干点模拟骨料表面凸起再在球内部随机生成少量点模拟内部核心全部点丢进ConvexHull最后把顶点坐标缩放和平移到目标位置。点数量控制在10-20个之间比较合适太少骨料形状太简单太多计算量增大且形状过度扭曲。注意骨料凹凸程度由表面点到球心的距离波动范围控制。波动范围占粒径的5%-20%超过30%会出现非常畸形的梭形骨料在网格划分时容易产生极小单元。3.3 纤维模型圆柱加半球帽的组合体纤维相对简单用圆柱体加两端半球帽近似保证纤维端部平滑避免尖锐棱角给网格带来麻烦。半径可控长度可控方向由两个随机角度确定。纤维的核心变量是长径比即长度与直径的比值。钢纤维混凝土中长径比通常在40-80之间短切纤维增强复合材料中可以到几百甚至上千这个参数直接决定纤维数量。因为纤维体积分数固定时长径比越大单根纤维体积越大总根数越少。生成纤维的代码逻辑是随机一个起点坐标确保整根纤维都在基体包围盒内随机生成两个角度theta和phi确定空间方向根据纤维长度和方向向量计算终点坐标在起点和终点位置创建两个半球帽在中间创建圆柱体需要注意纤维不能与基体边界相交否则嵌入约束会出错。我处理的办法是生成纤维前先判断起点和终点是否都在基体内部不在则重新生成。还有一种更精细的处理方式是把基体边界外的部分裁剪掉但裁剪操作会引入布尔运算的开销不如直接拒绝越界纤维简单高效。4. 碰撞检测与随机分布的核心算法4.1 多面体与多面体的碰撞分离轴定理实现骨料之间不能重叠这是细观模型最基本的要求。多面体碰撞检测我采用的是分离轴定理SAT逻辑很直观两个凸多面体不相交时一定存在一个分离平面使得两个多面体分别落在这个平面两侧。实际实现中只需要检测两类候选轴所有面的法向量和所有边方向两两叉积的结果。把所有顶点投影到候选轴上检查两段投影区间是否重叠只要某条轴上不重叠就说明两个多面体未碰撞。计算量分析一个骨料平均12个面两两检测就是12条面法向量加12×12条叉积轴总共156条轴每条轴要对所有顶点投影计算。如果生成50颗骨料最坏情况是50×49/2≈1225次两两检测每颗骨料顶点按20个算每次检测计算量约156×406240次操作整体量级在千万次Python跑起来可能要几秒钟。我在插件里加了一个阶段式检测的优化先用球体包围盒做粗判距离大于两倍最大粒径的直接跳过只有包围盒相交时才做精确SAT检测这样计算量能下降一个量级。提示凸包算法保证骨料是凸多面体是实现SAT的前提。凹多面体需要先做凸分解算法复杂度会大幅提升。实际骨料大多是近凸的凸包近似在工程上足够。4.2 纤维与多面体的碰撞线段与面相交测试纤维是细长圆柱体与骨料的碰撞检测其实可以简化为判断线段纤维轴线到多面体表面的距离是否小于纤维半径。如果纤维轴线离多面体表面足够远圆柱体必然不会穿透。具体做法是遍历骨料的每个三角面计算纤维线段到三角形的最短距离。最短距离如果小于纤维半径加一个安全余量就判定为碰撞。三角形三维空间中任意朝向的三角形到线段的最短距离算法在图形学里有标准实现思路是把线段参数化点在两个端点之间运动三角形区域用重心坐标描述构造一个6维优化问题用迭代法求最小值如果最小距离点落在线段端部或三角形边界上采用点到线段、点到三角形的距离公式这里要注意一个坑直接用ABAQUS自带的getDistance函数并不适用它只能算两个点之间的距离而我们需要的是线段到面的最小距离。我当时查了一下午资料最终决定自己实现这个几何工具函数并且在geometries.py里封装好。还有个实际工程问题当纤维穿过骨料时从几何上是可以做布尔减运算把纤维切成两段然后只保留基体外的部分。这个方案看起来很完美但ABAQUS布尔运算处理几十上百根纤维时极易失败而且切出来的小碎块网格质量极差。我在实际测试中放弃了这种做法转而采用禁止纤维与骨料相交的策略。这样做的代价是纤维分布更加稀疏一些但胜在模型稳定、网格质量可控。4.3 纤维与纤维的碰撞空间距离判断纤维之间碰撞检测比上面两个简单得多因为圆柱之间不相交等价于两线段之间的最短距离大于两倍半径。三维空间两条线段最近距离有经典的解析算法先计算了两条线段的参数化表达式通过求解4个线性方程得到最近点参数再根据参数是否落在[0,1]区间内裁剪修正。纤维数量通常比骨料多一个量级如果做两两全检测复杂度是O(n²)当纤维数量超过500时会显著拖慢生成速度。我的策略是空间哈希网格把基体立方体划分为若干小格子每根纤维只检测所在格子及相邻26个格子里的其他纤维。这样复杂度降到O(n)空间换时间很划算。格子大小选择有个技巧我把它设为直径的2到3倍。太大则格子内纤维数量多检测效率降低太小则一根纤维跨很多格子重复检测次数增加。实际测试下来格子边长取2.5倍直径时性能最好。4.4 随机分布策略与边界处理的细节随机分布算法本身不复杂难的是在高体积分数下保持高效的生成效率。我的做法是随机放置加最大尝试次数控制每颗骨料或纤维生成后检测碰撞失败则重新生成同一颗粒重新生成超过200次仍然失败就认为达到最大堆积密度跳出循环。这带来一个现实问题目标骨料体积分数50%以上时小粒径骨料很容易塞进空隙里但大粒径骨料尝试上百次也难以找到合适位置。遇到这种情况我在GUI里加了一个警告提示同时在代码里做了粒径分级投放先投放大粒径骨料再投放小粒径骨料填补空隙。这个策略模仿了实际混凝土的级配过程能有效提高体积分数的上限。最大堆积密度受很多因素影响包括粒径分布范围、形状规则程度、尝试次数上限、随机数质量等。根据我的测试经验采用级配投放配合100次尝试上限六面体基体内多面体骨料的体积分数上限大概在55%-65%之间。想突破这个上限需要引入周期性边界条件让骨料可以跨边界生成这在placement.py里预留了接口但目前还没有正式启用。边界处理还有一个细节骨料离边界太近时嵌入约束之后的网格过渡很不自然表面出现畸形单元。我在生成骨料后加了一个边界收缩操作把基体边界向内偏移一个安全距离在这个收缩后的区域内生成骨料保证骨料不会紧贴基体表面。5. 装配、嵌入约束与网格划分实操5.1 装配策略单一Part内布尔合并还是多Part装配生成几何体只是第一步ABAQUS里怎么把这些几何体组织成能分析的模型才是关键问题。两种方案方案A所有几何体在同一个Part内做布尔合并。好处是材料赋给一个Part就够了Embedded Region约束简单坏处是布尔合并操作极其消耗资源几百个实体合并经常失败或者卡死。方案B基体一个Part骨料和纤维分别生成Part然后装配在一起用Embedded Region嵌入。好处是几何生成和装配分开失败时定位容易坏处是嵌入约束需要逐个指定嵌入体和主区域。我最终选了方案B因为稳定性和错误定位对调试更重要。对于骨料我会先建立所有骨料Part再装配到装配体里。对于纤维如果全部当成独立Part数量太多上百根我采用了一个折中方案把同一批次生成的所有纤维合并成一个Part这个Part内部包含多个实体但对外是一整个Part。这样装配体里的Part数量很少同时嵌入约束的指定也简化了。5.2 Embedded Region约束的正确打开方式Embedded Region约束是细观模型的核心它让骨料和纤维作为嵌入体嵌在基体这个主区域里。ABAQUS会在分析时自动把嵌入体的节点约束到主区域的单元上相当于一种多点约束。操作步骤是在Interaction模块里创建Constraint类型选Embedded Region然后分别指定嵌入体区域和主区域。这里有几个容易踩的坑坑一嵌入体选择的是整个Part几何体还是某个Set要选Set不要选整个几何体。最好在Mesh模块里提前创建好Set比如All_AGG_Embedded和All_FIBER_Embedded分别包含所有骨料和纤维单元。如果直接选整个Part几何体ABAQUS有时候会把基体自身也嵌进去约束混乱。坑二嵌入体的单元类型必须正确。骨料用实体单元C3D4或C3D10纤维如果用实体建模也用实体单元如果用梁建模应该用T3D2桁架单元。混合使用实体和桁架单元在同一嵌入约束里完全可行但桁架单元嵌入后要确保节点被正确约束到主区域节点上检查信息文件里有没有warning。重要创建Embedded Region后一定要在Job模块提交分析前检查一下约束是否生效。方法是在Mesh模块里选择Display Group单独显示嵌入体单元看这些单元的位置是否正确有没有偏离到基体外。5.3 网格过渡与质量控制细观模型的网格划分是整个流程里最影响分析成败的环节。骨料和基体的网格尺寸差异过大时嵌入约束处的节点不匹配会严重影响计算精度。我的经验是基体网格尺寸取骨料最小粒径的1/3到1/4。比如骨料粒径范围5-10mm基体网格尺寸取1.5-2mm左右。纤维半径通常很小0.5mm左右纤维周围的基体网格必须足够密才能约束住纤维节点否则纤维只能被“空泡”约束受力完全传不进去。网格划分顺序建议是先划基体再划骨料和纤维。基体用结构化或扫掠网格骨料由于是不规则多面体只能使用自由网格Tet单元类型选C3D4线性四面体或者C3D10M修正二次四面体。C3D10M在处理接触和嵌入问题时比C3D4稳定很多代价是计算量大一些。纤维如果是实体建模的圆柱自由网格也会生成四面体单元数量会爆炸。所以我更推荐纤维用T3D2桁架单元建模几何上不再需要圆柱实体直接用线几何创建Part再赋予梁截面能节省大量单元。不过如果模拟目标是准确捕捉纤维与基体界面的应力传递还是应该用实体纤维配合Cohesive界面或Embedded约束做精细分析。网格质量检查要在Mesh模块里过一遍不要有负体积单元最小角度不要小于10度翘曲度别太大。碰到畸形单元时不要臆想通过加密网格解决先检查几何是不是有问题最常见的原因是骨料凸包上存在极薄的面或者极尖的角。这时候回到geometries.py里调整点集分布比在后处理修网格要高效得多。6. 材料属性定义与单元选择6.1 骨料、基体、纤维的材料本构细观模型里三种材料角色完全不同基体水泥砂浆或混凝土基体一般用塑性损伤本构Concrete Damaged PlasticityCDP描述拉压不对称行为。CDP需要输入受拉和受压的应力应变曲线、膨胀角、偏心率、双轴抗压强度比这些参数。如果没有试验数据可以先从ABAQUS官方手册的混凝土默认参数表取初值之后通过反演标定。骨料实际骨料强度远高于基体在普通荷载水平下骨料基本处于线弹性状态。直接用弹性本构定义弹性模量和泊松比就行除非模拟的是骨料破碎这种极端情况才需要引入脆性材料本构。纤维使用T3D2桁架单元配合弹塑性本构输入弹性模量、屈服强度和硬化参数。钢纤维的屈服强度通常在400MPa以上弹性模量200GPa左右注意单位统一。材料参数的单位一致性是个重量级陷阱。ABAQUS没有内置单位系统如果你用mm建模弹性模量应该用MPa密度用t/mm³力用N时间用s。很多人在这步出错得到的分析结果莫名其妙回头检查才发现是单位换算错误。我的插件在GUI里默认全部采用mm-N-MPa单位制尽可能减少转换工序。6.2 接触设定与界面行为嵌入约束和接触约束这两种方式处理界面时有着完全不同的力学假设。嵌入约束本质上是把纤维节点刚性绑定到基体单元不滑移、不脱粘适合宏观上模拟纤维的增强效果。接触约束Surface-to-Surface Contact允许界面的法向压力和切向摩擦能够模拟纤维滑移和脱粘但收敛难度剧增。我的插件目前默认使用Embedded Region原因很简单接触分析在几百个骨料界面存在时收敛极其困难。如果研究的主要目标是纤维掺量对整体强度的影响趋势嵌约束足矣要精确模拟纤维拔出和桥接作用建议后续改用Cohesive Behavior界面代价是建模和分析成本大幅提升。提示如果未来加入了Cohesive界面模型必须在网格划分前创建界面几何通常用骨料表面偏移一层零厚度Cohesive单元实现。这个工作在ABAQUS里用常规GUI做很繁琐也是插件后续扩展的重点方向。7. 插件GUI设计与参数校验7.1 界面布局与参数分组插件的GUI通过ABAQUS的RSG对话框构造器生成不需要手工写控件代码。界面分为三组基体参数、骨料参数、纤维参数外加一个全局参数区。基体参数就是三个几何尺寸长度、宽度、高度。骨料参数包括粒径最小值、粒径最大值、骨料体积分数、随机种子、凸包随机点数。纤维参数包括纤维直径、纤维长度、纤维体积分数、随机种子。这里有个经验之谈骨料和纤维分别设定独立的随机种子。这样做的用意很大调整纤维分布时不会影响骨料分布。比如做参数化研究时你先固定骨料种子只改纤维种子就能得到保持骨料完全相同、仅纤维乱向不同的系列模型非常有利于后期的统计分析。7.2 用户输入校验规则插件要做参数校验不能什么值都接收。我在param_utils.py里定义了一套校验规则粒径最小值必须是正数最大值必须大于最小值骨料体积分数必须大于0且小于一个最大允许值默认0.6纤维体积分数不能超过5%因为超过5%后随机分布的纤维极难达到体积分数目标大量纤维会聚集或者穿透纤维长度不能超过基体最小尺寸的0.8倍否则大量纤维无法找到合适位置随机种子必须是正整数参数不合法时插件会弹出明确错误消息指出哪个参数有问题而不是等生成到一半才崩溃。这些校验逻辑虽然不起眼但在被人高频使用时能省下大量调试时间。8. 实操流程从打开插件到提交分析8.1 一步步生成模型整个操作流程可以概括为七个步骤在ABAQUS主界面Plug-ins菜单里启动插件填写基体尺寸比如100×100×100 mm设定骨料粒径4-8mm体积分数0.4随机种子42设定纤维直径0.5mm长度20mm体积分数0.02随机种子2024点击Generate按钮等待进度条跑完在CAE里查看生成的部件和装配体赋材料属性划分网格创建Job提交计算步骤间的一个小坑插件运行过程中ABAQUS的GUI可能处于假死状态。这是正常现象因为几何生成和装配操作量大尤其是骨料数量超过100时这个过程可能需要几十秒到几分钟。我建议在插件里加上一个简单的进度提示虽然用纯Python实现进度条不太优雅但至少让用户知道程序没挂。8.2 实际生成的模型形态验证生成完模型后先自检再提交分析。我的自检方法是在Assembly模块里查看装配体用Render Style切换到Shaded显示看骨料和纤维是否均匀分布在基体内部用Query功能检查最密区域的骨料间隙最小间隙不应小于0.05mm随机选择几根纤维确认两端都在基体内部数值统计层面用Python脚本读取所有骨料和纤维的体积加和后除以基体体积与目标体积分数对比偏差控制在5%以内这些自检看起来不起眼但能挡住90%的无效计算。我曾经有一次忘了检查等Job算了两小时才发现纤维全部偏置到基体一侧原因只是随机数生成器的初始化种子设置错误从头再来损失惨重。9. 常见问题与排查技巧实录9.1 高体积分数下生成死循环症状插件长时间没反应CPU占用率100%但进度条不动。 原因碰撞检测退化为死循环尝试生成一个新骨料无论怎么随机都碰撞反复重试超过200次仍然失败。 排查思路第一时间检查目标骨料体积分数是否过高尝试减小最大粒径缩小粒径分布范围检查随机种子是否导致分布偏向某一侧。 解决办法插件内置的级配投放策略能改善但不完美最直接的办法是降低体积分数目标值。如果必须达到高填充率建议考虑改变骨料形状参数增加每颗粒的表面积让形状更接近球形即减少凸包点的数量从而提高空间利用率。经验分享我遇到过体积分数0.55时生成过程稳定0.58时突然开始死循环。用行列式统计了每次定位失败的骨料粒径发现几乎所有定位失败都发生在最大粒径级别。于是我把投放策略调整为按粒径从大到小分批投放大粒径先占位置小粒径填充缝隙解决问题。9.2 Embedded Region约束失败症状提交Job后报错“The host region for embedded elements is not defined”或者“Embedded elements are not contained within the host region”。原因分析最常见的是嵌入体和主区域的几何位置没有真实重合比如骨料由于碰撞检测精度不足实际上跟基体存在微小穿透或者间隙或者嵌入体的部分单元落到了主区域单元的外部。排查步骤检查几何生成阶段有没有骨料凸包面极度变形的情况导致表面局部穿透检查装配模块里基体Part和骨料Part的坐标系是否一致有没有发生意外偏移——我遇到过装配时不小心用了非原点坐标系所有骨料偏移了一个向量嵌入约束自然失败最后尝试重新划分主区域网格加密外部边界层。经验分享ABAQUS对嵌入约束的几何检查很严格嵌入体节点的容差设置太小也会失败。在Embedded Region约束编辑框中有一个Tolerance参数默认值很小必要时可以手工调大。但要注意调得过大会把本该在外面的节点也吸进来导致约束失真所以要谨慎操作。9.3 网格划分时单元畸变症状骨料网格划分成功但单元质量极差包含大量尖角退化单元有的单元甚至出现负Jacobian。原因骨料凸包表面存在极尖的角或者极薄的三角形面自由网格在这种情况下生成的四面体质量非常差。改进方案在SAT碰撞检测阶段加一个表面质量筛查剔除凸包表面三角形角度过小的骨料。具体判断标准是三角形最小内角小于5度的直接丢弃。这会稍微降低骨料生成效率但大幅提升网格质量。最有效的方法其实是优化凸包点集分布在生成表面点时增加一项随机扰动让表面趋于光滑同时避免生成非常接近共面的点。9.4 纤维只嵌入了一半症状渲染视图里纤维的一端嵌在基体内另一端悬在空气中像插在蛋糕上的牙签。原因纤维在生成时起点在基体内部但终点因为方向向量计算误差跑到了边界外。碰撞检测只检测了纤维与骨料碰撞没有检测纤维与基体边界的交互。修复方案在placement.py中增加边界检测函数判断纤维线段是否完全落在基体包围盒内部。我当时测试时发现单纯判断起点终点还不够因为纤维可能横跨基体角落部分跑出去。正确的做法是用线段与包围盒的裁剪算法Slab Method检测整根线段是否在包围盒内裁剪后如果线段长度变短了说明有部分越界拒绝这根纤维。碰到这种问题要记住纤维越界不会导致分析直接失败但会导致结果完全失真因为悬空部分的纤维节点没有被约束到任何基体单元上承载能力等于零而你以为它还在工作。9.5 参数微调后结果发生跳变症状只改了一个小参数比如纤维体积分数从2%改成2.1%算出来的强度反而大幅下降不规律跳变。原因随机种子没固定每次生成都是新的随机分布模型差异远大于参数本身的差异。解决办法固定随机种子是细观模拟的铁律。我在插件GUI里把随机种子放在最显眼的位置而且在参数校验时强制种子必须是正整数。提醒固定了种子也未必保证分布完全一致因为ABAQUS自身网格划分以及几何引擎在版本更新后可能产生细微差别这是无法控制的。10. Python源码结构解析10.1 模块划分与核心类的设计源码文件组织如下abaqus_agg_fiber/ ├── __init__.py # 插件入口注册GUI菜单 ├── plugin_gui.py # RSG对话框控件定义与参数回调 ├── generator.py # 主流程控制参数解析与模型组装 ├── geometries.py # 多面体(point cloud convex hull)、纤维的几何算法 ├── placement.py # 空间数据结构和碰撞检测系列算法 ├── embedder.py # 装配、Embedded Region约束创建封装 └── param_utils.py # 参数校验、默认值管理、单位转换geometries.py里最重要的类是Polyhedronvertices顶点坐标列表faces三角形面的顶点索引列表volume()计算体积translate()平移scale_random()随机缩放模拟粒径分布placement.py里最重要的类是SpatialHashGrid负责空间索引加速碰撞检测。所有几何体的包围盒信息会注册到这个哈希网格里检测碰撞时只遍历邻域格子内的候选体大大减少计算量。10.2 核心函数代码逐段分解举一个最核心的函数——generate_polyhedron_aggregate的伪代码分解def generate_polyhedron_aggregate(target_radius, surface_points12, irregularity0.2, seedNone): 生成一个随机多面体骨料。 target_radius: 目标粒径对应半径 surface_points: 表面随机点数 irregularity: 表面起伏幅度相对半径的比例 rng np.random.default_rng(seed) # 第一步在球面上生成随机点模拟骨料表面 # 生成方式随机单位向量乘以半径再乘以(1irregularity*随机偏移) directions rng.normal(size(surface_points, 3)) directions / np.linalg.norm(directions, axis1)[:, None] radii target_radius * (1.0 irregularity * rng.uniform(-1, 1, surface_points)) surface_pts directions * radii[:, None] # 第二步在球内部生成少量点模拟内部核心 inner_pts target_radius * 0.3 * rng.normal(size(4, 3)) pts np.vstack([surface_pts, inner_pts]) # 第三步凸包 hull ConvexHull(pts) # 第四步为了数值稳定性把顶点坐标稍微落在一个安全范围内 # 凸包可能产生非常接近于0面积的三角形这里做一次简化去重、去退化面 vertices, faces _simplify_hull(hull) return Polyhedron(vertices, faces)这个函数只用了NuMPy和SciPy不依赖ABAQUS API。可以独立测试。注意_simplify_hull这个函数特别关键它负责删除退化面避免凸包出现面积几乎为零的三角形这是我踩过无数坑后总结出来的经验。纤维生成函数相对简单def generate_fiber(length, radius, start_point, direction): 返回一个包含圆柱体和两个半球帽的Fibre对象。 # 圆柱部分两个底面圆心分别在start_point和end_point end_point start_point direction * length cylinder create_cylinder(start_point, end_point, radius) # 半球帽两个半球的圆心分别在两个端点 cap1 create_hemisphere(start_point, radius, direction) cap2 create_hemisphere(end_point, radius, -direction) return Fiber(cylinder, cap1, cap2, start_point, end_point)圆柱和半球帽在ABAQUS里最终都以Circle截面和Sweep生成几何体在布尔合并时会被处理成一个整体。要注意的是ABAQUS创建几何体用的是自家的内核通常基于Parasolid或ACIS直接用Python调用mdb.models[].Part和part.Cylinder这类命令生成基础几何再通过part.BooleanCut或part.BooleanMerge组合。10.3 碰撞检测代码的关键实现多面体碰撞用SAT核心代码结构def sat_collision(poly1, poly2): # 收集所有候选轴面法向量 边方向叉积 axes [] for face in poly1.faces: axes.append(face_normal(poly1, face)) for face in poly2.faces: axes.append(face_normal(poly2, face)) for edge1 in poly1.edges: for edge2 in poly2.edges: axes.append(np.cross(edge1, edge2)) # 去重方向相同或相反的轴合并 axes unique_axes(axes) # 每条轴上投影检测 for axis in axes: proj1 project_polyhedron(poly1, axis) proj2 project_polyhedron(poly2, axis) if not overlap(proj1, proj2): return False # 找到分离轴说明不相交 return True # 全部轴上都有重叠说明相交这里有一个非常影响性能的细节候选轴去重。如果不做去重两个各12面的多面体会有1212144168条轴很多轴方向相同或相反但向量数值不完全一样导致重复计算投影。去重后通常只剩30-50条独立轴计算量能再降一个量级。SAT准确测量两个凸多面体的距离还不行它只能返回是否碰撞。如果需要最小安全距离比如想让骨料之间保留0.05mm间隙需要在SAT基础上加一步每一次投影重叠区间宽度记录取所有轴上重叠宽度的最小值就是两个多面体的嵌入深度。嵌入深度小于间隙阈值就通过检测这个技巧我用来实现均匀间隙控制。10.4 随机数生成策略与可复现性保障插件使用numpy.random.default_rng(seed)作为核心随机数生成器。每个骨料或纤维生成时会从全局RNG读取一组随机值。这种做法的好处是可复现坏处是添加新逻辑会导致历史种子号失效。比如你固定种子42源程序生成50颗骨料。后来你在生成逻辑里多加了一个参数如骨料旋转角RNG的消耗序列变化种子42生成的骨料分布就完全不一样了。这在前向兼容性上很麻烦。我的做法是每个实体使用独立的子流base_rng np.random.default_rng(seed) for i in range(target_num): entity_seed int(base_rng.integers(1, 2**31 - 1)) entity_rng np.random.default_rng(entity_seed) # 用entity_rng生成几何这样即便修改了生成逻辑只要子流的分配顺序不变旧模型依然可以复现。这个设计值得实现它让我在做模型回溯时省了很多力气。11. 扩展能力与实际应用效果11.1 更换几何类型从多面体到椭球、球形geometries.py和placement.py分离设计的最大好处是扩展类型极其容易。想改用椭球骨料只需要在geometries.py里新增一个generate_ellipsoid_aggregate函数返回一个Polyhedron的近似版本——用密集网格化的椭球面三角剖分来近似椭球然后placement.py里的SAT检测依然适用。想改用球骨料更简单Polyhedron退化为一个正二十面体就能很好的近似球体而且计算量骤然降低因为候选轴数量大幅减少。11.2 周期性边界条件的预留细观模拟里周期性边界条件非常重要它能让模型边界处的骨料和纤维分布不再受制于边界效应。我在placement.py里预留了一个use_periodic参数当开启后随机生成时可以把实体复制到对侧边界思想是在空间上“卷绕”。触发条件是把一个实体越出包围盒的部分镜像到对侧等效于在无限周期域中放置实体。这个功能目前在生产代码里还没完全启用因为镜像后的实体与基体之间的布尔运算极其复杂ABAQUS的几何引擎处理周期边界镜像特征时效率很低。如果后续用单位胞方法做均匀化分析这个功能会变得非常必要。11.3 插件的实际应用案例我用手头这套插件做过几个方向的研究一是混凝土单轴压缩细观模拟基体100×100×100mm骨料体积分数0.4粒径4-8mm纤维掺量1%-3%不等固定种子号做系列对比。嵌入约束条件下纤维掺量从0到2%时峰值强度提升约25%-40%和文献趋势一致。二是氯离子扩散模拟这个场景稍微不同因为不需要力学分析只需要几何模型我把插件生成的装配体导出为DXF或STL再导入到COMSOL做扩散分析。省掉了在COMSOL里手动生成随机骨料的步骤效率提升明显。三是界面过渡区ITZ建模骨料和基体之间的ITZ层厚度通常50微米左右直接几何建模比较困难。我用插件的碰撞检测间隙控制功能让骨料之间保持固定间距再在骨料表面偏移一层薄壳模拟ITZ这个思路目前正在试验。11.4 性能数据在标准工作站Intel i7-1270032GB内存上的实际耗时基体100³mm骨料40颗纤维50根总生成时间约8秒基体300³mm骨料350颗纤维200根总生成时间约3分钟网格划分耗时很长通常比几何生成慢3-5倍C3D10M网格50万单元约耗时5分钟提交单个Job在隐式分析下单轴压缩模型约3万迭代步耗时数小时到一天取决于接触和收敛情况12. 个人经验总结与后续想法写这个插件的过程让我体会到一件事ABAQUS二次开发里最耗时的不是算法设计而是把算法和ABAQUS的几何内核、装配机制、网格划分完美好地衔接起来。最初我用一个O(n²)的两两检测来验证正确性一天之内跑完1000颗粒的模型就很吃力。后来意识到直接用空间哈希一定能提速但为了省事没尽早实现反而因为单次生成耗时过久影响了大量调试工作。从那以后我的原则变成了能一次验证正确性的逻辑不拖到后面再去简化性能。目前插件还是以“生成模型”为主材料属性和网格参数仍然需要自己在CAE里设定。下一步我想把材料库也集成进去让CDP本构的默认参数自动填入网格参数也从GUI传入真正做到一键建模、一键出网格。对于想做类似功能的人我的建议是先明确目标你需要的到底是能看的三维几何模型还是能直接提交算力分析的完整CAE模型。这两者的工作量和坑的密度差一个数量级前者半天写完后者我折腾了一个月。如果你也想写类似插件强烈建议优先固定随机种子这是所有可复现性工作的基石。先把一个简单的球形骨料生成器跑通再做多面体化、碰撞检测、嵌入约束每步都验证过再往后走。代码量不大但每一步都有暗坑。