简介面向混凝土细观模拟、随机骨料建模与三维图形学研究这份Python代码基于vor3d库生成3D Voronoi多面体用于模拟混凝土中骨料颗粒的随机分布与不规则形态适合材料科学、建筑工程、数值仿真及图形学方向的开发者参考。资源包内容高度聚焦仅含1个py脚本文件压缩包体积约1KB便于快速下载与阅读源码。目前已有528人学习可用于快速掌握Voronoi分区的核心实现思路。脚本完整覆盖了种子点随机生成、3D Voronoi结构构建、多面体边界处理与可视化等关键步骤并对无限面剔除、骨料尺寸与间隙控制等实际问题给出可扩展的代码框架读者可在此基础上调整参数将其接入有限元分析或混凝土性能预测流程为后续结构模拟提供几何建模基础。1. 混凝土随机骨料建模为什么 3D Voronoi 是第一选择做混凝土细观数值试验的人都会卡在同一道坎试件里那几百颗棱角分明的骨料怎么放进代码。画球太假CT 扫描太贵手工拼多面体太慢而 3D Voronoi 生成随机骨料恰好把三个问题一起绕开。先随机布点再做 Voronoi 剖分每个种子自动长成一颗凸多面体天然铺满试件、天然不重叠、尺寸可控收缩一步就得到带浆体间隙的骨料形貌。下文把布点、剖分、裁剪、收缩、导出整条链路拆开讲代码和参数都能直接抄适合做混凝土细观有限元、离散元模拟也适合批量生成骨料样本喂给强度预测模型的工程师。2. 3D Voronoi 多面体为什么像真实骨料从 Delaunay 对偶到体积分数2.1 Voronoi 剖分的几何定义每颗骨料都是「最近种子点的领地」给定一组种子 s₁…sₙVoronoi 胞的定义是空间中离某个种子比其他所有种子都近的点集V(sᵢ) {x | ‖x−sᵢ‖ ≤ ‖x−sⱼ‖, ∀j ≠ i}。每个胞是若干半空间的交集所以一定是凸多面体相邻两个胞的公共面是种子连线的垂直平分面。这也是为什么 Voronoi 骨料看起来像破碎石料而不是鹅卵石——破碎、断裂形成的骨料表面本来就是平面占主导而 Voronoi 胞的每一个面都是平面。三维泊松随机布点均匀随机种子产生的 Voronoi 胞平均面数在 15 个左右、平均棱数约 40正好落在破碎石灰岩和花岗岩骨料的棱面数统计区间内。这个统计上的接近不是巧合Voronoi 胞的顶点由三邻种子构成、棱由两邻种子构成面数和棱长分布完全由种子密度与随机性决定天然具备「碎块感」。相比之下在球面上随机取点再取凸包得到的多面体面数少、形态偏圆而且没有天然的尺度分布拿来当骨料几何既不像、也不好控制级配。2.2 铺满空间是杀手锏Voronoi 把「堆积问题」变成「抽样问题」用球、椭球或任意凸多面体做随机骨料核心难题一直是堆积要在正方体内放进体积分数 40% 以上的互不重叠颗粒得跑 DEM 模拟或写专门的随机顺序堆积算法速度慢、收敛差、边界效应难处理。Voronoi 剖分不存在这个问题——剖分完成后整块域被胞完全铺满覆盖率 100%没有任何重叠也没留任何孔隙。于是骨料体积分数不再靠几何试错而是一次精确算术先铺满全部胞占 100% 体积再按目标体积分数保留部分胞或对所有胞做统一收缩把缩掉的体积留给砂浆。我一般用收缩这条路因为每个胞都保留下来骨料空间分布均匀不会出现局部过密或过稀后续界面层厚度的控制也直接落在收缩系数上。2.3 等效粒径与富勒级配Voronoi 胞的尺寸怎么映射到筛分曲线Voronoi 胞不是球得定义等效粒径才能跟筛分级配对上。常用的是等体积球径 d_eq (6V/π)^(1/3)把每个胞的体积折算成直径再排序。混凝土细观建模的级配通常直接引富勒曲线累计通过率 P(d) (d/D_max)^nn 取 0.45~0.70D_max 是最大骨料粒径。混凝土类型典型骨料体积分数建模备注C30 普通0.40–0.43多峰级配细骨料占比高C50 中等0.45–0.47细观模型默认值C80 高强0.50–0.55骨料偏细界面层更薄砂浆对照组0用于分离界面效应一个经常踩的换算坑实验报告给的是质量分数或面积分数转体积分数要乘密度比或对面积分数做 3/2 次幂修正。CT 切片上看到 30% 的骨料面密度对应体积分数约 45%别直接拿切片数值当目标。3. vor3d 最小可运行流程硬核布点、周期剖分、收缩导出不同版本的 vor3d 工具接口可能完全不同有的是命令行、有的是 Python 库但剖分链路是一致的布点 → 周期剖分 → 跨界处理 → 收缩 → 导出。下面用最通用的 scipy 实现把链路走通换到任何封装都只需要替换函数签名。3.1 布点别用裸 np.random硬核拒绝采样控制最小粒径均匀随机种子在三维空间里会随机扎堆扎堆处会生成极薄的多面体胞三角面退化成细条后面网格划分直接报错。要给种子之间设一个最小间距 min_gap做拒绝采样import numpy as np from scipy.spatial import Voronoi, ConvexHull def hardcore_seeds(n, box, gap, max_try200_000): 在 [0,box]^3 内生成 n 个最小间距为 gap 的种子点 pts [] tries 0 while len(pts) n and tries max_try: p np.random.uniform(0, box, 3) if all(np.linalg.norm(p - q) gap for q in pts): pts.append(p) tries 1 if len(pts) n: raise RuntimeError(f只布了 {len(pts)}/{n} 个点减小 gap 或增大 box) return np.array(pts)min_gap 直接决定最小粒径的下限gap 越小允许越细的骨料gap 越大细骨料直接被拒掉级配整体变粗。n 与 gap 要一起试先跑 50 颗种子看 d_eq 分布再定最终参数。拒绝采样在 n 到几千时没问题到几万层就要换成 Bridson 泊松盘采样或六面体预划分否则 while 循环退化成超长尾。3.2 周期延拓剖分27 个镜像种子封住开边界直接对试件里的种子做 scipy Voronoi边界上的胞会向无穷远开口。常见做法是把种子沿三个方向各平移 −1/0/1 倍边长复制出 26 组镜像凑成 3×3×327 组一起剖分这样每个胞都有完整邻居边界被镜像种子封住def periodic_voronoi_polys(seeds, box): 对中心装箱内的种子生成周期 Voronoi返回(顶点,三角面,体积)列表 shifts np.array([[i, j, k] for i in (-1, 0, 1) for j in (-1, 0, 1) for k in (-1, 0, 1)]) all_pts np.concatenate([seeds g * box for g in shifts]) # 27 倍点集 vor Voronoi(all_pts, qhull_optionsQbb Qc Qz) n len(seeds) polys [] for i in range(n): reg vor.regions[vor.point_region[13 * n i]] # 中心域块 if -1 in reg: continue v vor.vertices[reg] tol 1e-9 inside (v -tol).all(1) (v box tol).all(1) if not inside.all(): continue # 跨界胞直接丢弃表面留砂浆层见 4.4 if len(v) 4: continue hull ConvexHull(v) polys.append((v, hull.simplices.copy(), hull.volume)) return polysshifts 列表里 (0,0,0) 是第 13 个元素所以中心域种子在 all_pts 里的索引是 13ni这个偏移写错最常见的症状是胞错乱、导出后骨料互相切割。qhull_options 的 Qbb 把输入归一化到单位包围盒Qc 保留共面点Qz 处理近共面退化这三项对三维 Voronoi 基本是标配。27 倍镜像的内存开销是硬伤n600 约 1.6 万点没问题n5000 就是 13.5 万点scipy 的时间与内存都开始吃紧这个规模我一般换 voro 或 neper 这类 C 周期剖分工具算法逻辑不变。3.3 收缩与体积分数标定V′ s³V先统计保留胞的体积分数再按目标值反推收缩系数。由于绕任意点均匀缩放体积都按 s³ 变化这一步是精确的box 100.0 # 试件边长单位 mm seeds hardcore_seeds(600, box, 4.0) polys periodic_voronoi_polys(seeds, box) tot sum(p[2] for p in polys) vf_keep tot / box ** 3 # 去掉跨界胞后的保留率一般 0.5~0.8 Vt 0.45 # 目标骨料体积分数 s (Vt / vf_keep) ** (1.0 / 3) # 收缩系数例如 0.6^(1/3)≈0.84 cells [] for vtx, fcs, _ in polys: c vtx.mean(0) vs c (vtx - c) * s # 绕质心近似点均匀收缩体积精确乘 s³ hull ConvexHull(vs) cells.append((vs, hull.simplices.copy(), hull.volume))不要担心vtx.mean(0)不是真质心会导致体积公式失效——均匀缩放绕任何点做体积都乘 s³标定不受影响。真正要算的是丢胞补偿跨界胞被丢弃后 vf_keep 低于 1直接用目标 Vt 去比会整体偏小所以 s 必须用 vf_keep 反推这正是上面代码做的事。参数常用范围主要影响box50–150 mm试件尺寸建议 ≥ 5×D_max 控制边界效应n300–3000骨料总颗数与平均粒径min_gap2–8 mm最小粒径与最差网格单元质量s0.70–0.95浆体/界面层厚度与最终体积分数缩完立刻做一道自检sum(c.volume)/box³应等于 Vt误差 1% 以内对不上就是丢胞比例估错了或缩放写错。3.4 导出 OBJ 与可视化直接进 ParaView 和 gmshdef write_obj(path, cells, header#): with open(path, w) as f: f.write(f{header}\n) base 1 for vtx, fcs, _ in cells: for x, y, z in vtx: f.write(fv {x:.6f} {y:.6f} {z:.6f}\n) for tri in fcs: f.write(ff {base tri[0]} {base tri[1]} {base tri[2]}\n) base len(vtx) write_obj(aggregates.obj, cells, header# seed600 gap4.0 vf0.45)OBJ 里 v 是顶点、f 是三角面索引base 变量负责把每颗胞的顶点索引平移到全局偏移。ParaView 可视化可以走 pyvista注意面索引要跟着顶点拼接一起重映射import pyvista as pv v_all, f_all, base [], [], 0 for vtx, fcs, _ in cells: v_all.append(vtx) f_all.append(fcs base) base len(vtx) mesh pv.PolyData(np.vstack(v_all), np.hstack([[3, *t] for t in np.vstack(f_all)])) mesh.plot()画完发现破洞优先怀疑共面顶点产生的重复三角形如果骨料表面出现非流形边检查顺序是「去重 → 面积阈值过滤 → 统计边共用次数」。进有限元前 OBJ 还要做四面体化gmsh 里按逐胞导入再网格化或者直接走 Abaqus 的部件装配流程。4. 随机骨料生成的 3 个核心参数与 4 个典型坑min_gap、qhull、收缩系数4.1 min_gap 与最小粒径的换算先试算再定标无法从 gap 一步解析出 d_eq_min但统计上 gap 与最小胞尺寸同量级gap4mm 时 d_eq 分布基本从 5~8mm 起跳600 颗种子在 100mm 试件里等效粒径大致落在 5~25mm。正确做法是先小规模试算deq np.sort([(6 * p[2] / np.pi) ** (1/3) for p in polys]) print(fd_eq 范围 {deq[0]:.2f} ~ {deq[-1]:.2f} mm, 中位 {np.median(deq):.2f} mm)若需要明确的级配区间例如 5~20mm 双档用「生成后筛留」更直接按 d_eq 排序保留落在档位区间内的胞再重新标定体积分数。靠反复调 gap 去精确命中级配曲线是浪费时间的搜索直接按目标曲线抽样种子或筛留胞都更可控。提示种子间距分布直接决定最终胞形。想要长条状骨料在布点阶段对某一轴做坐标缩放想要更圆的骨料收缩前对顶点做 Laplace 光顺。这两个操作都放在剖分之后、收缩之前。4.2 qhull 数值问题共面、退化与重复面ConvexHull 输出三角剖分面当胞内出现四个以上共面顶点时同一平面被拆成多个三角形导出后部分三角形面积趋近零。处理办法是去重加面积阈值过滤clean [] for vtx, fcs, _ in cells: fcs np.unique(np.sort(fcs, axis1), axis0) tri vtx[fcs] area np.linalg.norm(np.cross(tri[:, 1] - tri[:, 0], tri[:, 2] - tri[:, 0]), axis1) / 2 fcs fcs[area 1e-6] clean.append((vtx, fcs, ConvexHull(vtx[fcs].reshape(-1, 3)).volume))退化面清完再做流形检查统计每条边被几个三角形共用封闭曲面应为 2。from collections import Counter edge_cnt Counter() for _, fcs, _ in clean: for a, b, c in fcs: edge_cnt[tuple(sorted((a, b)))] 1 edge_cnt[tuple(sorted((b, c)))] 1 edge_cnt[tuple(sorted((a, c)))] 1 bad [e for e, k in edge_cnt.items() if k ! 2] print(非流形边数量:, len(bad))出现 1 或 3 次共用的边说明该处网格有裂缝或翻转多半是共面顶点没清干净回到上一步而不是直接改收缩系数。4.3 收缩系数与 ITZ 厚度的换算真实混凝土界面过渡区 ITZ 厚度约 20~50 μm细观模型通常放大到 0.5~1.0 mm 以便网格捕捉。收缩产生的间隙 g 近似为 g ≈ (1−s)·R_effR_eff 是该胞等效半径 d_eq/2。给定目标 g0.6mm、平均粒径 12mm则 s ≈ 0.9反过来先定 s 再统计间隙分布两组数对不上就要查是不是有胞收缩后仍相碰。相碰检查用表面采样每颗胞表面取 100 个随机点做 cKDTree 最近邻查询最小距离应大于 0 且接近设计的 g。4.4 四类典型翻车症状、原因、改法症状根因修改动作网格剖分内存爆掉27 倍镜像加种子过多减少 n或换 voro/neper 周期剖分骨料两两刺穿s 太大或 min_gap 太小减小 s或增大 min_gap 重布点级配偏细没有大骨料泊松 Voronoi 尺寸分布天然偏细按 d_eq 筛留大胞或按级配曲线布点表面缺骨料、整体 Vf 偏低跨界胞被丢弃且未补偿丢胞后按 vf_keep 重算 s循环 2~3 次OBJ 破面、重复面共面顶点退化三角形去重 面积阈值过滤 流形边检查最后一行最容易被忽略先定体积分数、再生成几何、最后算网格。踢回重来的顺序是「删跨界胞 → 重算 vf_keep → 重算 s → 重收缩」循环两三次就收敛别手工硬凑。5. 用 CT 扫描代码校准 vor3d批量产出抗压强度数据集5.1 拿 CT 分割结果校核统计特征合成骨料不能自说自话真实样本的统计特征是验收基准。CT 扫描代码的标准链路是滤波去噪 → 阈值分割 → 连通域标记 → marching cubes 提面 → 算球度。import numpy as np from skimage import filters, measure from scipy import ndimage vol load_ct_volume() # 自己的 CT 体数据 mask vol filters.threshold_otsu(vol) mask ndimage.binary_opening(mask, iterations2) lab measure.label(mask) for p in measure.regionprops(lab): if p.volume 500: # 体素个数滤噪 continue single (lab p.label) verts, faces, _, _ measure.marching_cubes(single, 0.5) A measure.mesh_surface_area(verts, faces) V p.volume * voxel_vol # voxel_vol 单个体素体积 psi np.pi ** (1/3) * (6 * V) ** (2/3) / A # Wadell 球度每一颗分割骨料的球度、纵横比、凸度直方图就是 vor3d 合成骨料的对照基准。绝大多数情况下 Voronoi 胞的球度会偏高说明棱角不够锐利把 Laplace 光顺步数减小或干脆不做球度自然回落。若 CT 骨料有明显扁平主轴则在布点阶段做各向异性缩放。CT 的职责不是给几何而是给统计矩。5.2 批量生成训练样本把第三章脚本包成gen_mesostructure(random_seed, Vt, gap)一个函数循环里只改随机种子其余全固定——体积分数、级配、网格尺寸、材料本构、加载速率。每个样本跑一次单轴压缩有限元把峰值应力、峰值应变、断裂能、破坏模式标签和对应几何统计特征d_eq 均值、球度、体积分数写进同一行 CSV。这就是「混凝土抗压强度数据集」的一种自产路径真实数据不够用细观模拟批量补关键是让输入特征和标签来自同一套可控几何而不是把不同来源的数据拼在一起。5.3 两个五分钟验收技巧批量生成的返工成本最高所以跑循环前固定两道自检第一道体积守恒sum(vol)/box³与目标 Vt 偏差小于 1%第二道最小间隙cKDTree 表面采样最近邻统计最小值不为负。两道都过才进有限元。最后再安利一个小习惯把随机种子、gap、s、Vt 全部写进 OBJ 文件头注释就是 3.4 里的 header 参数。样本一旦多起来文件名先撑不住注释里留参数任何一颗骨料都能回溯到生成现场训练数据出异常时靠这一行注释能省下半天对账时间。本文还有配套的精品资源点击获取