简介这是一套基于ABAQUS的二维随机颗粒生成插件面向从事颗粒材料有限元模拟的工程师与科研人员。核心功能通过Python脚本调用ABAQUS内建函数自动完成颗粒尺寸、形状、材料属性及分布区域的参数化设置并借助随机数算法生成位置方向均无序的颗粒模型免去手动逐个建模的繁琐流程。压缩包共1个文件类型为py脚本大小仅2KB体量轻巧、即下即用适合在颗粒填充、堆积、接触行为及应力传递等仿真场景中快速搭建初始模型。目前已有657人学习使用者只需掌握基础的ABAQUS操作和Python语法即可上手。借助该工具可大幅缩短前处理时间同时通过ABAQUS后处理验证颗粒排列与边界条件提升建模效率与仿真精度。1. 为什么二维随机颗粒模型不能靠手摆在ABAQUS里做颗粒增强复合材料、散体介质或土体细观力学分析时最耗时间的步骤通常不是求解而是建模。手动放置十几个圆形颗粒还能接受一旦目标数量到几百颗半径还服从分布再逐个画圆、装配、定义接触光几何前处理就能耗掉半天。suijishengchengkeli.py 这个插件在 ABAQUS 的 Python 接口里实现了一个二维随机颗粒生成器输入区域尺寸、粒径范围、颗粒间距和目标数量脚本自动生成互不重叠的圆心与半径并转成可分析的 Part 和 Instance。它适合做细观仿真、写论文需要控制随机种子复现结果的工程师和研究生也适合做 CAE 二次开发、想把建模流程脚本化的分析人员。用脚本生成随机颗粒表面上是“撒点”真正的难点是让颗粒不重叠、体积分数可预期、随机序列可复现。接下来的内容按“算法原理 → 脚本实现 → 运行调试 → 验证技巧”的顺序展开。2. 颗粒生成的底层逻辑随机、干涉和几何约束2.1 随机数种子与粒径分布常见的做法是把随机过程分成两个独立部分颗粒几何随机和位置随机。粒径可以通过random.uniform(r_min, r_max)获得连续均匀分布如果需要模拟级配材料也可以换成random.gauss(mu, sigma)或按对数正态分布。这里的关键是random.seed(seed)固定种子后每次运行生成的颗粒完全一致这一特性对参数敏感性分析和论文复现非常重要。如果不传种子每次生成结果都不相同适合做批量工况的统计学分析。ABAQUS 的内核脚本运行的是标准 Python 解释器random模块直接可用不需要额外安装第三方包。颗粒位置通常按区域矩形范围均匀采样。位置采样区域建议比颗粒半径留出一个安全边距这样生成的颗粒不会跨出模型边界。如果你需要模拟特定级配可以先用筛分曲线拟合出累积分布函数再用逆变换采样得到每个颗粒的半径但插件默认使用均匀分布因为它在粒径信息不充分时不会引入额外偏差。2.2 颗粒重叠检测圆心距判据随机生成的颗粒如果只是简单撒点两个圆有概率重叠尤其在体积分数较高时。对于圆颗粒判断两个圆是否重叠只需看圆心距是否小于半径之和加最小间隙。这个判据在二维和三维中同样适用只不过三维中要判断球心距。代码上我习惯把“是否与已有颗粒重叠”抽成独立函数def is_overlap(cx, cy, r, existing_particles, min_gap): for px, py, pr in existing_particles: if (cx - px) ** 2 (cy - py) ** 2 (r pr min_gap) ** 2: return True return False这段代码的意思很直接遍历所有已生成的颗粒计算新的圆心距离的平方再与半径之和加上最小间隙后的平方比较。用平方比较是为了省去开平方操作在颗粒数量上千时能省不少时间。min_gap通常取一个大于 0 的数值比如比模型最小网格尺寸还大一点这样后续在 ABAQUS 中划分网格时不会因为颗粒间距离过近而产生退化单元。如果颗粒数量需求很高找不到足够位置这个函数会一直重试。因此生成逻辑里必须设置最大尝试次数否则脚本可能陷入死循环。实际项目中目标数量超过 500 时我会把尝试次数上限放宽到 100 万次同时打印进度避免用户误以为程序卡死。2.3 从随机几何到 ABAQUS 部件生成坐标和半径之后需要在 ABAQUS 中建立几何。二维颗粒最常见的做法是把每个颗粒单独建为一个可变形壳部件再由rootAssembly生成 Instance。单独建 Part 的原因有两个一是不同颗粒可以赋予不同材料属性或接触属性二是便于在后处理中单独观察某个颗粒的应力和位移。如果只是做纯统计研究所有颗粒材料相同也可以把整个区域建成一个 Part 再画多个圆但那样颗粒间接触定义会比较麻烦所以不推荐。在这个阶段ABAQUS 的ConstrainedSketch和Part命令被反复调用。为了避免脚本执行时在界面里不断弹出绘图窗口建议把草图名称设为临时前缀或者使用 noGUI 方式运行脚本。产品中集成这个插件时也尽量不要在绘制草图时启动显示否则上千个颗粒会显著拖慢界面响应。2.4 体积分数与目标数量换算实际使用插件时用户输入的是颗粒数量而不是体积分数。我一般先根据目标体积分数估算数量n floor(vf * region_area / (pi * r_mean^2))。但直接这样写有个隐患如果半径分布较宽平均半径计算出的数量可能导致实际体积分数偏高。所以更稳妥的做法是生成后计算真实体积分数再微调min_gap或num_particles。这个换算关系可以帮助你不盲目试参数。比如区域是 100×60目标体积分数 25%平均半径 2.0那么先设颗粒数大约为0.25 * 6000 / (3.1416 * 4) ≈ 119。脚本生成后如果实际体积分数偏低说明颗粒间空隙偏大可以适当降低min_gap或增加颗粒数量。3. suijishengchengkeli.py 脚本实现与参数设置3.1 脚本整体流程与参数表suijishengchengkeli.py 的工作流程可以归纳为读取参数 → 初始化随机数 → 循环生成候选颗粒 → 重叠检查 → 生成 Part 和 Instance → 返回实例列表。我通常不把材料属性写进生成函数因为随机几何和物理属性是两件事。材料参数留给后面的截面和材料模块处理这样脚本可以复用到不同分析类型中从二维平面应力到平面应变都能用。生成阶段需要用户输入的参数如下参数名类型说明我的建议region_width / region_heightfloat生成区域宽高应大于 2 倍最大粒径r_min / r_maxfloat颗粒半径范围按筛分数据或占地比推算min_gapfloat颗粒间最小间隙建议 ≥ 一个网格尺寸num_particlesint目标颗粒个数根据体积分数和粒径反算seedint随机种子填固定值可复现填 None 真随机shapestr颗粒形状当前版本主要为圆形可扩展到多边形这张表里最容易被忽略的是区域尺寸与粒径的关系。如果region_width比2 * r_max还小random.uniform生成坐标时会出现下界大于上界脚本会抛出ValueError。所以生成函数里要先做合法性检查。3.2 核心生成函数下面的函数是整个插件的核心它把“随机撒点 拒绝采样”完整串起来import random def generate_random_particles(region_width, region_height, r_min, r_max, min_gap, num_particles, seedNone): if seed is not None: random.seed(seed) if region_width 2 * r_max or region_height 2 * r_max: raise ValueError(区域尺寸必须大于2倍最大半径) particles [] max_attempts num_particles * 2000 attempts 0 while len(particles) num_particles and attempts max_attempts: attempts 1 r random.uniform(r_min, r_max) x random.uniform(r, region_width - r) y random.uniform(r, region_height - r) overlap False for px, py, pr in particles: if (x - px) ** 2 (y - py) ** 2 (r pr min_gap) ** 2: overlap True break if not overlap: particles.append((x, y, r)) if len(particles) num_particles: raise RuntimeError(无法在给定区域内生成足够的非重叠颗粒) return particles代码逻辑先固定随机种子再检查区域尺寸然后用while循环不断尝试生成候选颗粒。每次尝试都只随机一个半径和一个圆心位置再和已有颗粒逐个比较距离。如果完全不重叠就追加到列表否则丢弃这次尝试继续下一轮。max_attempts是防御性写法避免区域已经塞满时陷入死循环。我把上限设为目标颗粒数乘以 2000这个倍数在多数二维圆颗粒场景下够用如果体积分数很低实际尝试次数会远小于上限。random.uniform(r, region_width - r)这一步保证了颗粒整体都在区域内。注意这里的边界收缩量是当前颗粒自己的半径而不是r_max所以小颗粒可以在更靠近边界的位置出现这更合理。还有一个可选优化把半径序列先按从大到小排序再逐个放置因为大颗粒更碍事先放大的能显著减少后续尝试次数。但为了保持统计均匀性默认函数没有做排序。3.3 把颗粒列表转为 ABAQUS 部件得到颗粒列表后还需把它变成 ABAQUS 可用的 Part。下面的函数基于 abaqus Python API 实现from abaqus import mdb from abaqusConstants import * def create_particle_instances(model_name, region_width, region_height, particles): model mdb.models[model_name] assembly model.rootAssembly instances [] for idx, (cx, cy, r) in enumerate(particles): sketch model.ConstrainedSketch( name__profile__%d % idx, sheetSizemax(region_width, region_height)) sketch.CircleByCenterPerimeter(center(cx, cy), point(cx r, cy)) part model.Part(nameParticle-%03d % idx, dimensionalityTWO_D_PLANAR, typeDEFORMABLE_BODY) part.BaseShell(sketchsketch) inst assembly.Instance(nameInst-%03d % idx, partpart, dependentON) instances.append(inst) return instances这个函数遍历所有颗粒坐标为每个颗粒创建独立的草图圆和二维可变形部件。CircleByCenterPerimeter需要两个点圆心和圆上任意一点这里用(cx r, cy)作为半径方向的参考点。BaseShell把草图转换成一个可赋截面属性的二维壳几何。如果后续要赋予面内厚度并画网格这种二维壳可以搭配平面应力单元使用。装配时使用dependentON表示实例依赖部件几何之后对 Part 的网格和截面积分控制会自动同步到所有实例。对于需要单独给某个颗粒设置材料的情况可以在每个 Part 上先做好截面赋值再装配。使用dependent比OFF更节省装配内存接触定义也更稳定所以在颗粒数量较多时我一般保持默认依赖。3.4 扩展颗粒形状和方向除了圆形颗粒形状可以扩展到椭圆或多边形。对椭圆而言除了圆心和长短轴还需要随机旋转角度。生成时重叠检测从“圆心距比较半径之和”变成“椭圆距离函数”实现成本会高不少。我一般用外接圆做快速丢弃再精算真实重叠。这个插件保留 shape 参数就是给这类扩展留接口。如果你只需要圆形直接复用上面两个函数即可。4. 在 ABAQUS 中运行插件与调试4.1 三种运行方式脚本写好之后有三种常见运行方式在 CAE 交互界面点击 File → Run Script选择 suijishengchengkeli.py。这种方式可以在脚本执行后直接在视口中看到生成的颗粒实例适合少量颗粒。使用命令行 noGUI 运行abaqus cae -noGUI suijishengchengkeli.py。这种方式不启动图形窗口适合批量生成和服务器环境。把脚本注册为插件菜单项通过 Abaqus Plugin Builder 打包成.rsg文件。这需要额外编写 GUI 定义适合交付给不熟悉 Python 的用户。我推荐先用第 2 种方式调试因为它不会受到图形界面重绘的影响。生成上千个颗粒时交互界面里显示每个 Part 草图会非常慢noGUI 模式只打印输出信息速度明显更快。如果使用第 1 种方式脚本运行前最好先打开一个新的模型数据库保证活动模型名是Model-1否则mdb.models[Model-1]会报 KeyError。4.2 常见运行错误与处理以下是这个插件在运行中常遇到的问题和我的处理思路错误现象可能原因处理方式mdb.models[Model-1]报 KeyError脚本执行时没有活动模型脚本开头建立新模型或从已有 CAE 文件导入ValueError: uniform(a,b)下界大于上界区域尺寸小于 2 倍r_max检查region_width和r_max的关系生成颗粒数少于预期min_gap过大或num_particles过高减小min_gap或加大区域重新计算体积分数界面卡成假死状态草图对象过多CAE 重绘太久改用 noGUI 运行或降低颗粒数量屏幕上出现 libpng error显示缓存或图形驱动异常关闭 CAE 图形窗口用 noGUI 模式重跑其中“界面卡死”在颗粒数量超过 300 时比较容易出现。ABAQUS 每次新建 Part 都会触发图形渲染即使把草图名隐藏也不能完全避免。碰到这种情况不要反复点鼠标不让 CAE 继续重绘而是直接启动任务管理器结束进程。还有一个常被问到的“中断不了”问题脚本在 while 循环里卡住时CAE 菜单中的停止按钮并不会中断正在执行的内核脚本我通常直接关闭 CAE 重启或在提交 job 前用mdb.jobs[Job-1].kill()控制任务生命周期。4.3 与网格和材料赋值的衔接几何生成完成后下一步是赋值网格和材料。在主模型里我会把粒子实例按粒径分组例如半径大于某个阈值的颗粒作为“粗颗粒”单独设置一种材料其余作为“细颗粒”。这样能减少后续接触定义时需要手动选取的对象的数量。一个常用的衔接片段是在生成实例后统一给所有 part 创建壳截面material_name Grain if material_name not in model.materials: model.Material(namematerial_name) model.materials[material_name].Elastic(table((210e3, 0.3),)) for name, part in model.parts.items(): if name.startswith(Particle-): sec_name Section_ name model.HomogeneousShellSection(namesec_name, materialmaterial_name, thickness1.0) faces part.faces part.Set(nameAll, facesfaces) part.SectionAssignment(regionpart.sets[All], sectionNamesec_name)这段脚本给所有以Particle-开头的部件创建弹性材料并赋予壳截面。thickness1.0对二维平面应力问题只是名义厚度不影响面内应力结果如果是平面应变问题厚度值同样可以设为 1。如果遇到“截面类型不支持”的错误可以先检查 Part 的维度是否真的是TWO_D_PLANAR以及是否存在孤立网格没有几何。网格种子尺寸建议设置在min_gap的十分之一以内避免颗粒之间因为距离过近生成高畸变单元。颗粒与基体之间的接触定义属于另一层工作这里跳过。5. 进阶用统计量验证随机模型质量5.1 体积分数与最小间距核验生成几何后我建议先用一段小脚本做数学验证而不是直接进求解器。体积分数是二维圆颗粒最重要的统计量它等于全部圆面积除以区域面积。最小间距则用来判断是否满足接触参数假设。下面的脚本可以打印这两项import math def check_model(particles, region_width, region_height): area_region region_width * region_height area_grain sum(math.pi * r * r for _, _, r in particles) volume_fraction area_grain / area_region min_dist float(inf) for i in range(len(particles)): for j in range(i 1, len(particles)): dx particles[i][0] - particles[j][0] dy particles[i][1] - particles[j][1] d math.hypot(dx, dy) - particles[i][2] - particles[j][2] if d min_dist: min_dist d print(体积分数: %.4f % volume_fraction) print(最小间隙: %.6f % min_dist)体积分数和你在 ABAQUS 里看到的面积占比应该一致最小间隙应大于等于min_gap的设定值。如果最小间隙接近 0说明颗粒几乎相切网格划分阶段可能无法生成高质量的单元。5.2 边界和可视化检查清单在 CAE 里我习惯按这个清单检查生成结果颗粒边缘是否超出区域边界尤其注意小颗粒靠近边界的情况最小间隙是否全部满足通过测量两个相邻颗粒的最近点随机种子是否被固定如果两次运行结果不同说明 seed 没有生效颗粒数是否与体积分数匹配避免出现“区域看起来稀疏、但统计值偏高”的错觉。对于只关心统计规律的模型还可以用最近邻距离分布来描述随机性。颗粒数量足够多时最近邻距离的均值会落在自由体积理论预测区间附近。如果偏差太大就需要怀疑min_gap是不是塞得太满导致颗粒位置偏向区域中心。5.3 保存和复用生成参数最后一个小技巧每次生成颗粒后把颗粒列表输出成 CSV 文件下次需要同一模型时直接读取 CSV不需要重新随机。这个习惯可以大量节省调试时间。写一个简单的导出函数就是几行代码import csv def export_particles(particles, filename): with open(filename, w, newline) as f: writer csv.writer(f) writer.writerow([x, y, radius]) writer.writerows(particles)在脚本末尾调用export_particles就可以把随机几何保存为可追溯的“颗粒布置文件”。之后在 ABAQUS 之外用 Plot 或 Matlab 复核也能直接读同一个文件。恢复时用csv.reader读入并调用create_particle_instances即可。本文还有配套的精品资源点击获取
简介这是一套基于ABAQUS的二维随机颗粒生成插件面向从事颗粒材料有限元模拟的工程师与科研人员。核心功能通过Python脚本调用ABAQUS内建函数自动完成颗粒尺寸、形状、材料属性及分布区域的参数化设置并借助随机数算法生成位置方向均无序的颗粒模型免去手动逐个建模的繁琐流程。压缩包共1个文件类型为py脚本大小仅2KB体量轻巧、即下即用适合在颗粒填充、堆积、接触行为及应力传递等仿真场景中快速搭建初始模型。目前已有657人学习使用者只需掌握基础的ABAQUS操作和Python语法即可上手。借助该工具可大幅缩短前处理时间同时通过ABAQUS后处理验证颗粒排列与边界条件提升建模效率与仿真精度。1. 为什么二维随机颗粒模型不能靠手摆在ABAQUS里做颗粒增强复合材料、散体介质或土体细观力学分析时最耗时间的步骤通常不是求解而是建模。手动放置十几个圆形颗粒还能接受一旦目标数量到几百颗半径还服从分布再逐个画圆、装配、定义接触光几何前处理就能耗掉半天。suijishengchengkeli.py 这个插件在 ABAQUS 的 Python 接口里实现了一个二维随机颗粒生成器输入区域尺寸、粒径范围、颗粒间距和目标数量脚本自动生成互不重叠的圆心与半径并转成可分析的 Part 和 Instance。它适合做细观仿真、写论文需要控制随机种子复现结果的工程师和研究生也适合做 CAE 二次开发、想把建模流程脚本化的分析人员。用脚本生成随机颗粒表面上是“撒点”真正的难点是让颗粒不重叠、体积分数可预期、随机序列可复现。接下来的内容按“算法原理 → 脚本实现 → 运行调试 → 验证技巧”的顺序展开。2. 颗粒生成的底层逻辑随机、干涉和几何约束2.1 随机数种子与粒径分布常见的做法是把随机过程分成两个独立部分颗粒几何随机和位置随机。粒径可以通过random.uniform(r_min, r_max)获得连续均匀分布如果需要模拟级配材料也可以换成random.gauss(mu, sigma)或按对数正态分布。这里的关键是random.seed(seed)固定种子后每次运行生成的颗粒完全一致这一特性对参数敏感性分析和论文复现非常重要。如果不传种子每次生成结果都不相同适合做批量工况的统计学分析。ABAQUS 的内核脚本运行的是标准 Python 解释器random模块直接可用不需要额外安装第三方包。颗粒位置通常按区域矩形范围均匀采样。位置采样区域建议比颗粒半径留出一个安全边距这样生成的颗粒不会跨出模型边界。如果你需要模拟特定级配可以先用筛分曲线拟合出累积分布函数再用逆变换采样得到每个颗粒的半径但插件默认使用均匀分布因为它在粒径信息不充分时不会引入额外偏差。2.2 颗粒重叠检测圆心距判据随机生成的颗粒如果只是简单撒点两个圆有概率重叠尤其在体积分数较高时。对于圆颗粒判断两个圆是否重叠只需看圆心距是否小于半径之和加最小间隙。这个判据在二维和三维中同样适用只不过三维中要判断球心距。代码上我习惯把“是否与已有颗粒重叠”抽成独立函数def is_overlap(cx, cy, r, existing_particles, min_gap): for px, py, pr in existing_particles: if (cx - px) ** 2 (cy - py) ** 2 (r pr min_gap) ** 2: return True return False这段代码的意思很直接遍历所有已生成的颗粒计算新的圆心距离的平方再与半径之和加上最小间隙后的平方比较。用平方比较是为了省去开平方操作在颗粒数量上千时能省不少时间。min_gap通常取一个大于 0 的数值比如比模型最小网格尺寸还大一点这样后续在 ABAQUS 中划分网格时不会因为颗粒间距离过近而产生退化单元。如果颗粒数量需求很高找不到足够位置这个函数会一直重试。因此生成逻辑里必须设置最大尝试次数否则脚本可能陷入死循环。实际项目中目标数量超过 500 时我会把尝试次数上限放宽到 100 万次同时打印进度避免用户误以为程序卡死。2.3 从随机几何到 ABAQUS 部件生成坐标和半径之后需要在 ABAQUS 中建立几何。二维颗粒最常见的做法是把每个颗粒单独建为一个可变形壳部件再由rootAssembly生成 Instance。单独建 Part 的原因有两个一是不同颗粒可以赋予不同材料属性或接触属性二是便于在后处理中单独观察某个颗粒的应力和位移。如果只是做纯统计研究所有颗粒材料相同也可以把整个区域建成一个 Part 再画多个圆但那样颗粒间接触定义会比较麻烦所以不推荐。在这个阶段ABAQUS 的ConstrainedSketch和Part命令被反复调用。为了避免脚本执行时在界面里不断弹出绘图窗口建议把草图名称设为临时前缀或者使用 noGUI 方式运行脚本。产品中集成这个插件时也尽量不要在绘制草图时启动显示否则上千个颗粒会显著拖慢界面响应。2.4 体积分数与目标数量换算实际使用插件时用户输入的是颗粒数量而不是体积分数。我一般先根据目标体积分数估算数量n floor(vf * region_area / (pi * r_mean^2))。但直接这样写有个隐患如果半径分布较宽平均半径计算出的数量可能导致实际体积分数偏高。所以更稳妥的做法是生成后计算真实体积分数再微调min_gap或num_particles。这个换算关系可以帮助你不盲目试参数。比如区域是 100×60目标体积分数 25%平均半径 2.0那么先设颗粒数大约为0.25 * 6000 / (3.1416 * 4) ≈ 119。脚本生成后如果实际体积分数偏低说明颗粒间空隙偏大可以适当降低min_gap或增加颗粒数量。3. suijishengchengkeli.py 脚本实现与参数设置3.1 脚本整体流程与参数表suijishengchengkeli.py 的工作流程可以归纳为读取参数 → 初始化随机数 → 循环生成候选颗粒 → 重叠检查 → 生成 Part 和 Instance → 返回实例列表。我通常不把材料属性写进生成函数因为随机几何和物理属性是两件事。材料参数留给后面的截面和材料模块处理这样脚本可以复用到不同分析类型中从二维平面应力到平面应变都能用。生成阶段需要用户输入的参数如下参数名类型说明我的建议region_width / region_heightfloat生成区域宽高应大于 2 倍最大粒径r_min / r_maxfloat颗粒半径范围按筛分数据或占地比推算min_gapfloat颗粒间最小间隙建议 ≥ 一个网格尺寸num_particlesint目标颗粒个数根据体积分数和粒径反算seedint随机种子填固定值可复现填 None 真随机shapestr颗粒形状当前版本主要为圆形可扩展到多边形这张表里最容易被忽略的是区域尺寸与粒径的关系。如果region_width比2 * r_max还小random.uniform生成坐标时会出现下界大于上界脚本会抛出ValueError。所以生成函数里要先做合法性检查。3.2 核心生成函数下面的函数是整个插件的核心它把“随机撒点 拒绝采样”完整串起来import random def generate_random_particles(region_width, region_height, r_min, r_max, min_gap, num_particles, seedNone): if seed is not None: random.seed(seed) if region_width 2 * r_max or region_height 2 * r_max: raise ValueError(区域尺寸必须大于2倍最大半径) particles [] max_attempts num_particles * 2000 attempts 0 while len(particles) num_particles and attempts max_attempts: attempts 1 r random.uniform(r_min, r_max) x random.uniform(r, region_width - r) y random.uniform(r, region_height - r) overlap False for px, py, pr in particles: if (x - px) ** 2 (y - py) ** 2 (r pr min_gap) ** 2: overlap True break if not overlap: particles.append((x, y, r)) if len(particles) num_particles: raise RuntimeError(无法在给定区域内生成足够的非重叠颗粒) return particles代码逻辑先固定随机种子再检查区域尺寸然后用while循环不断尝试生成候选颗粒。每次尝试都只随机一个半径和一个圆心位置再和已有颗粒逐个比较距离。如果完全不重叠就追加到列表否则丢弃这次尝试继续下一轮。max_attempts是防御性写法避免区域已经塞满时陷入死循环。我把上限设为目标颗粒数乘以 2000这个倍数在多数二维圆颗粒场景下够用如果体积分数很低实际尝试次数会远小于上限。random.uniform(r, region_width - r)这一步保证了颗粒整体都在区域内。注意这里的边界收缩量是当前颗粒自己的半径而不是r_max所以小颗粒可以在更靠近边界的位置出现这更合理。还有一个可选优化把半径序列先按从大到小排序再逐个放置因为大颗粒更碍事先放大的能显著减少后续尝试次数。但为了保持统计均匀性默认函数没有做排序。3.3 把颗粒列表转为 ABAQUS 部件得到颗粒列表后还需把它变成 ABAQUS 可用的 Part。下面的函数基于 abaqus Python API 实现from abaqus import mdb from abaqusConstants import * def create_particle_instances(model_name, region_width, region_height, particles): model mdb.models[model_name] assembly model.rootAssembly instances [] for idx, (cx, cy, r) in enumerate(particles): sketch model.ConstrainedSketch( name__profile__%d % idx, sheetSizemax(region_width, region_height)) sketch.CircleByCenterPerimeter(center(cx, cy), point(cx r, cy)) part model.Part(nameParticle-%03d % idx, dimensionalityTWO_D_PLANAR, typeDEFORMABLE_BODY) part.BaseShell(sketchsketch) inst assembly.Instance(nameInst-%03d % idx, partpart, dependentON) instances.append(inst) return instances这个函数遍历所有颗粒坐标为每个颗粒创建独立的草图圆和二维可变形部件。CircleByCenterPerimeter需要两个点圆心和圆上任意一点这里用(cx r, cy)作为半径方向的参考点。BaseShell把草图转换成一个可赋截面属性的二维壳几何。如果后续要赋予面内厚度并画网格这种二维壳可以搭配平面应力单元使用。装配时使用dependentON表示实例依赖部件几何之后对 Part 的网格和截面积分控制会自动同步到所有实例。对于需要单独给某个颗粒设置材料的情况可以在每个 Part 上先做好截面赋值再装配。使用dependent比OFF更节省装配内存接触定义也更稳定所以在颗粒数量较多时我一般保持默认依赖。3.4 扩展颗粒形状和方向除了圆形颗粒形状可以扩展到椭圆或多边形。对椭圆而言除了圆心和长短轴还需要随机旋转角度。生成时重叠检测从“圆心距比较半径之和”变成“椭圆距离函数”实现成本会高不少。我一般用外接圆做快速丢弃再精算真实重叠。这个插件保留 shape 参数就是给这类扩展留接口。如果你只需要圆形直接复用上面两个函数即可。4. 在 ABAQUS 中运行插件与调试4.1 三种运行方式脚本写好之后有三种常见运行方式在 CAE 交互界面点击 File → Run Script选择 suijishengchengkeli.py。这种方式可以在脚本执行后直接在视口中看到生成的颗粒实例适合少量颗粒。使用命令行 noGUI 运行abaqus cae -noGUI suijishengchengkeli.py。这种方式不启动图形窗口适合批量生成和服务器环境。把脚本注册为插件菜单项通过 Abaqus Plugin Builder 打包成.rsg文件。这需要额外编写 GUI 定义适合交付给不熟悉 Python 的用户。我推荐先用第 2 种方式调试因为它不会受到图形界面重绘的影响。生成上千个颗粒时交互界面里显示每个 Part 草图会非常慢noGUI 模式只打印输出信息速度明显更快。如果使用第 1 种方式脚本运行前最好先打开一个新的模型数据库保证活动模型名是Model-1否则mdb.models[Model-1]会报 KeyError。4.2 常见运行错误与处理以下是这个插件在运行中常遇到的问题和我的处理思路错误现象可能原因处理方式mdb.models[Model-1]报 KeyError脚本执行时没有活动模型脚本开头建立新模型或从已有 CAE 文件导入ValueError: uniform(a,b)下界大于上界区域尺寸小于 2 倍r_max检查region_width和r_max的关系生成颗粒数少于预期min_gap过大或num_particles过高减小min_gap或加大区域重新计算体积分数界面卡成假死状态草图对象过多CAE 重绘太久改用 noGUI 运行或降低颗粒数量屏幕上出现 libpng error显示缓存或图形驱动异常关闭 CAE 图形窗口用 noGUI 模式重跑其中“界面卡死”在颗粒数量超过 300 时比较容易出现。ABAQUS 每次新建 Part 都会触发图形渲染即使把草图名隐藏也不能完全避免。碰到这种情况不要反复点鼠标不让 CAE 继续重绘而是直接启动任务管理器结束进程。还有一个常被问到的“中断不了”问题脚本在 while 循环里卡住时CAE 菜单中的停止按钮并不会中断正在执行的内核脚本我通常直接关闭 CAE 重启或在提交 job 前用mdb.jobs[Job-1].kill()控制任务生命周期。4.3 与网格和材料赋值的衔接几何生成完成后下一步是赋值网格和材料。在主模型里我会把粒子实例按粒径分组例如半径大于某个阈值的颗粒作为“粗颗粒”单独设置一种材料其余作为“细颗粒”。这样能减少后续接触定义时需要手动选取的对象的数量。一个常用的衔接片段是在生成实例后统一给所有 part 创建壳截面material_name Grain if material_name not in model.materials: model.Material(namematerial_name) model.materials[material_name].Elastic(table((210e3, 0.3),)) for name, part in model.parts.items(): if name.startswith(Particle-): sec_name Section_ name model.HomogeneousShellSection(namesec_name, materialmaterial_name, thickness1.0) faces part.faces part.Set(nameAll, facesfaces) part.SectionAssignment(regionpart.sets[All], sectionNamesec_name)这段脚本给所有以Particle-开头的部件创建弹性材料并赋予壳截面。thickness1.0对二维平面应力问题只是名义厚度不影响面内应力结果如果是平面应变问题厚度值同样可以设为 1。如果遇到“截面类型不支持”的错误可以先检查 Part 的维度是否真的是TWO_D_PLANAR以及是否存在孤立网格没有几何。网格种子尺寸建议设置在min_gap的十分之一以内避免颗粒之间因为距离过近生成高畸变单元。颗粒与基体之间的接触定义属于另一层工作这里跳过。5. 进阶用统计量验证随机模型质量5.1 体积分数与最小间距核验生成几何后我建议先用一段小脚本做数学验证而不是直接进求解器。体积分数是二维圆颗粒最重要的统计量它等于全部圆面积除以区域面积。最小间距则用来判断是否满足接触参数假设。下面的脚本可以打印这两项import math def check_model(particles, region_width, region_height): area_region region_width * region_height area_grain sum(math.pi * r * r for _, _, r in particles) volume_fraction area_grain / area_region min_dist float(inf) for i in range(len(particles)): for j in range(i 1, len(particles)): dx particles[i][0] - particles[j][0] dy particles[i][1] - particles[j][1] d math.hypot(dx, dy) - particles[i][2] - particles[j][2] if d min_dist: min_dist d print(体积分数: %.4f % volume_fraction) print(最小间隙: %.6f % min_dist)体积分数和你在 ABAQUS 里看到的面积占比应该一致最小间隙应大于等于min_gap的设定值。如果最小间隙接近 0说明颗粒几乎相切网格划分阶段可能无法生成高质量的单元。5.2 边界和可视化检查清单在 CAE 里我习惯按这个清单检查生成结果颗粒边缘是否超出区域边界尤其注意小颗粒靠近边界的情况最小间隙是否全部满足通过测量两个相邻颗粒的最近点随机种子是否被固定如果两次运行结果不同说明 seed 没有生效颗粒数是否与体积分数匹配避免出现“区域看起来稀疏、但统计值偏高”的错觉。对于只关心统计规律的模型还可以用最近邻距离分布来描述随机性。颗粒数量足够多时最近邻距离的均值会落在自由体积理论预测区间附近。如果偏差太大就需要怀疑min_gap是不是塞得太满导致颗粒位置偏向区域中心。5.3 保存和复用生成参数最后一个小技巧每次生成颗粒后把颗粒列表输出成 CSV 文件下次需要同一模型时直接读取 CSV不需要重新随机。这个习惯可以大量节省调试时间。写一个简单的导出函数就是几行代码import csv def export_particles(particles, filename): with open(filename, w, newline) as f: writer csv.writer(f) writer.writerow([x, y, radius]) writer.writerows(particles)在脚本末尾调用export_particles就可以把随机几何保存为可追溯的“颗粒布置文件”。之后在 ABAQUS 之外用 Plot 或 Matlab 复核也能直接读同一个文件。恢复时用csv.reader读入并调用create_particle_instances即可。本文还有配套的精品资源点击获取
)
)
