简介这是一份面向ABAQUS用户的Voronoi图生成代码资源适用于多晶体材料、颗粒复合材料、裂纹扩展与颗粒分布等有限元分析场景。Voronoi图本质上是基于种子点的空间划分又称泰森多边形可在ABAQUS中用于构建具有随机几何特征的非均质模型。压缩包共8个文件包含4个ipynb交互式脚本、3个txt说明文档和1个py核心脚本包体仅159KB轻量易用其中ipynb适合逐步演示生成流程py脚本可直接运行或修改txt文档用于记录依赖与参数说明。目前已有522人浏览学习。读者借助这些脚本可完成种子点定义、Voronoi区域计算、有限元网格划分、边界条件设置及后处理分析形成一套从几何生成到力学计算的完整链路。代码结构清晰便于有Python和ABAQUS脚本基础的用户直接复用、二次开发或嵌入自有模型对初学者而言也可通过跟随交互式脚本理解Voronoi图与有限元结合的基本原理是兼顾实用与学习价值的轻量工具包。1. 拿到 Abaqus-Voronoi 工程包之前先弄清楚这套东西解决什么问题拿到 Abaqus-Voronoi--master.zip 这类工程包最常见的动作是解压、找 .inp、提交作业然后被网格或材料报错劝退。Voronoi 剖分把随机分布的种子点变成凸多边形或多面体每个单元代指一个晶粒于是 Abaqus 里就能做晶粒尺度的力学、热学与断裂分析金属多晶各向异性、沿晶开裂、焊接热影响区组织都从这个几何基础出发。标题拆开看Abaqus 是求解器Voronoi 是晶粒几何算法master.zip 意味着一套可二次开发的模板工程。适合做过常规各向同性有限元、接下来要处理晶粒取向、晶界本构和微结构统计的工程师。下面按几何生成、网格与材料、晶界与边界、报错与验证的顺序走通链路每一步都给可复现的命令和参数。2. Voronoi 几何生成种子点规则决定晶粒分布两条路线对照2.1 剖分算法保证什么、不保证什么Voronoi 剖分的定义很简洁给定 n 个种子点空间里每个点被划分给离它最近的那个种子。二维情形下每个晶粒是一个凸多边形三维是凸多面体。它保证晶粒之间无缝隙、无重叠晶界由相邻种子点的中垂线构成几何上是完备的但不保证晶粒尺寸的均匀性。用完全均匀随机分布的种子点生成晶粒面积会呈近似指数分布小晶粒偏多这与真实金属里常见的对数正态晶粒尺寸分布相差很远。所以建模第一步不是写网格代码而是先设定目标晶粒尺寸分布。常见做法是先按对数正态分布采样等效直径 diameq再反推种子点密度和分布范围。遇到大晶粒包围小晶粒的组织可以在生成剖分后做一次 Lloyd 松弛迭代把每个种子移到自己单元的质心再重新剖分让结构趋近 centroidal Voronoi。这样晶粒尺寸更均匀、单元形状也更接近正多边形后续网格质量会明显改善Abaqus 里出现畸形单元的几率大幅下降。2.2 路线一Neper 一键生成带周期性的多晶网格如果目标是几百个晶粒以上的代表性体积单元 RVE我一般直接用开源工具 Neper而不是自己写剖分和网格代码。Neper 分两步先用 -T 模块生成剖分再用 -M 模块加密网格并输出 Abaqus inp。# 生成 200 个晶粒、等效直径服从对数正态分布的周期剖分 neper -T -n 200 -morpho diameq:lognormal(15,0.3) -periodic all -o poly # 在上述剖分基础上生成网格输出 Abaqus 输入文件 neper -M -n 200 -periodic all -format inp -o poly第一条命令里 -n 指定晶粒个数-morpho 里的 diameq:lognormal(15,0.3) 让等效直径的均值和形状参数按对数正态设定-periodic all 表示三个方向都做周期剖分。这是做 RVE 时最值钱的一个选项它保证相对面的网格节点一一对应后面加周期性位移边界条件时不用自己配对节点。第二条命令的 -format inp 直接把 mesh 写成 poly.inp里面已经包含节点、单元、晶粒 elset 和周期节点对应关系。Neper 输出的 inp 里单元按晶粒分好了 elset省去手工归类步骤。需要注意两个地方一是它默认生成四面体单元想换六面体要显式指定单元类型参数二是如果后续要插入 cohesive 晶界Neper 的共享面信息不会自动写进 inp需要在后处理时自己重建晶粒邻居关系这一点放到第四章展开。另外输出版本较老的 Neper 生成的 inp 里关键字缩进风格和 Abaqus 新版本有细微差异提交前用 abaqus datacheck 跑一遍最稳。2.3 路线二Python scipy.spatial.Voronoi 自己造几何不用 Neper、只想在现有 Python 流程里快速出模型时scipy.spatial.Voronoi 是最短路径。但它有两个坑边界上的单元会退化到无穷远必须裁剪完全随机的种子会让晶粒大小起伏过大。下面这段是带 3×3 镜像扩展的最小实现镜像的目的是消除非周期边界让中心区域的晶粒在边界处自然闭合。import numpy as np from scipy.spatial import Voronoi def voronoi_tiled(seeds, L(100.0, 100.0)): 3x3 平铺种子点做 Voronoi中心区域得到周期连续的多边形。 tiles np.array([(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)]) offset tiles * np.array(L) sup (seeds[None, :, :] offset[:, None, :]).reshape(-1, 2) return Voronoi(sup) rng np.random.default_rng(7) seeds rng.uniform(0, 100, size(80, 2)) vor voronoi_tiled(seeds) # vor.regions 中的多边形顶点在 vor.vertices 中 # 取所有顶点都落在 [0,100]x[0,100] 内的区域即为可用晶粒这段代码把种子点复制到周围 8 个相邻区块对 9 倍密度的点集做一次 Voronoi然后只保留落在中心区块内的有限多边形。做出来的晶界在边界处与对侧晶界自然连续等效实现了周期剖分。参数上seeds 的个数决定平均晶粒面积L 是 RVE 边长两者配合决定晶粒数密度如果做晶粒尺寸敏感性分析建议固定 L只改种子的对数正态采样参数否则几何变化和尺寸变化耦合在一起后处理里分不清是谁引起的力学响应差异。拿到多边形之后网格生成这一步 scipy 帮不上忙。常见做法是把多边形集合导出成 Gmsh 几何.geo或者直接在 Abaqus 的 Part 模块用 Partition Face 逐个分区再自由网格。我习惯对每个晶粒单独铺三角形再合并确保网格不跨晶界这个细节直接关系到第四章 cohesive 单元的生成。2.4 两条路线怎么选看晶粒数、周期性和后续修改频率Neper 和 Python 自建路线的差别主要在周期支持、网格生成和中间文件透明度。按我的使用经验对比没有绝对优劣看场景选即可。对比项NeperPython scipy周期剖分-periodic all 原生支持需 3×3 平铺自行实现晶粒尺寸分布diameq:lognormal(μ,σ) 直接指定需自己采样种子并做松弛网格划分内置网格器直接出 inp只出几何网格借 Gmsh 或 Abaqus二次开发CLI 参数控制脚本友好全 Python可嵌入现有数据流典型适用大型 RVE、晶体塑性参数反演教学演示、快速原型、对接实验数据选型时我的判断标准是晶粒数少于 300、并且最终要自己加 cohesive 或自己控制晶粒取样的场景Python 路线更顺手晶粒数上千、还要保证周期性和网格均匀性Neper 性价比高得多。还有一个容易被忽略的点Neper 输出的晶粒边界是理想化平面实验重构的 EBSD 数据往往是任意多面体如果你的模型要逐晶粒校准取向两条路线都得额外写数据映射层把 EBSD 取向表按晶粒号关联到 elset。3. 把 Voronoi 装进 AbaqusPart 拓扑、inp 拼装与晶粒取向3.1 先定模型拓扑单 Part 多 Section 和多 Part 装配的区别Voronoi 多晶模型在 Abaqus 里最常见的失败原因是拓扑选错。多 Part 装配方案里每个晶粒是独立 Part界面网格必然不共节点后期要靠 *TIE 绑定这会引入主从面选择和接触搜索的麻烦晶界 cohesive 也难做。我一般不用这个方案而是把整个 RVE 建成一个 Part所有晶粒共享连续网格晶粒之间的边界只是 elset 划分的产物。在 CAE 里实现这套拓扑的标准做法是先建一个矩形或立方体 Part用 Sketcher 导入 Voronoi 多边形线框再在 Part 模块用 Partition Face 按多边形逐个分区最后统一划分网格。网格跨晶粒自动共节点这是后续一切晶粒级分析的前提。划分时不要让单元穿过晶界2D 用 CPS3 或 CPS6M 三角形时边界天然贴合3D 建议用四面体并按晶粒设置局部种子密度界面两侧尺寸突变过大会在晶界处产生畸形单元求解时先报错的地方往往就在这里。3.2 从 Voronoi 网格直接写出可提交的 .inp不想在 CAE 里手工分区时把上一章 Python 生成的多边形直接写成 inp 是最可控的路径。下面这段把节点、单元、晶粒 elset 和取向一起写进文件以 2D 三角形单元 CPS3 为例。def write_voronoi_inp(path, polys): polys: [(verts_array, euler_angle), ...]输出单 Part 的 Abaqus inp。 node_id, elem_id 1, 1 lines, nid_map [*NODE], {} grain_elems, grain_angles [], [] for verts, angle in polys: ids [] for x, y in verts: key (round(float(x), 6), round(float(y), 6)) if key not in nid_map: nid_map[key] node_id lines.append(f{node_id}, {x:.6f}, {y:.6f}) node_id 1 ids.append(nid_map[key]) # 扇形式三角化以第一个顶点为公共点拆分凸多边形 tris [(elem_id k, ids[0], ids[k], ids[k 1]) for k in range(1, len(ids) - 1)] elem_id len(tris) grain_elems.append(tris) grain_angles.append(angle) lines.append(*ELEMENT, TYPECPS3, ELSETALLELEM) for tris in grain_elems: for eid, a, b, c in tris: lines.append(f{eid}, {a}, {b}, {c}) for gi, (tris, angle) in enumerate(zip(grain_elems, grain_angles), 1): lines.append(f*ELSET, ELSETGRAIN_{gi}) lines.append(, .join(str(t[0]) for t in tris)) lines.append(f*SOLID SECTION, ELSETGRAIN_{gi}, MATERIALCU_SINGLE) lines.append(f*ORIENTATION, NAMEORI_{gi}, SYSTEMRECTANGULAR) lines.append(f{angle:.4f}, 0.0, 0.0, 0.0, 1.0, 0.0) with open(path, w) as f: f.write(\n.join(lines) \n)写 inp 时最容易错的是 elset 与单元号的对应。这段脚本在三角化时直接按晶粒记录单元号列表再拼 *ELSET 行单元归属和几何生成在同一次循环里完成不会错位。节点去重用了六位小数的坐标作为字典键浮点噪声会被 round 吸收。*ORIENTATION 的两个方向点定义局部坐标系的 X 轴和 XY 平面晶粒的欧拉角通过这个局部坐标系施加到各向异性材料上如果做晶体塑性还要在这里配 *DEPVAR 和 UMAT先用单晶弹性矩阵验证取向是否生效再叠加晶体塑性本构排错面会小很多。3.3 晶粒取向、单晶弹性矩阵和 Schmidt 因子检查材料卡片上多晶模型的起点是单晶立方弹性。铜、铝这类面心立方金属用 *ELASTIC, TYPEANISO 时给 C11、C12、C44 三个独立常数即可。晶粒取向的正交性检查是新手最容易忽略的*ORIENTATION 定义的局部坐标必须保证三轴正交否则 Abaqus 会警告并强制正交化取向就偏离你输入的欧拉角。各晶粒的 inp 数据块关系可以用一张表理清。数据块作用常见错误*NODE节点坐标浮点重复导致共节点失败*ELEMENT单元连接晶界穿越、单元反转*ELSET晶粒归属单元号错位取向张冠李戴*ORIENTATION晶粒局部坐标三轴不正交被强制修正提交前用便宜的办法验证取向是否正确把材料暂时设成各向同性弹性给模型单向拉伸位移边界后处理里看每个晶粒的应力分布。因为各晶粒弹性刚度相同但方向不同晶粒内应力会呈现与取向相关的明暗差异。若所有晶粒应力场完全相同基本可以断定 *ORIENTATION 没被 *SOLID SECTION 引用到或者欧拉角转换时多了 90 度偏差。这种检查在焊接仿真里同样适用热力耦合计算前先跑一个纯力学小模型确认晶粒应力响应符合预期再叠加温度场排错成本会低很多。4. 焊接仿真里的 Abaqus-Voronoi 组合cohesive 晶界、热源与周期边界4.1 焊接热影响区组织演化的 Voronoi 近似abaqus 焊接仿真里Voronoi 最多被用在焊缝金属和热影响区 HAZ 的晶粒组织重建。真实焊缝凝固组织是柱状晶HAZ 是粗化的等轴晶两者可以用不同种子参数的两套 Voronoi 拼合近似焊缝区种子密度低、晶粒沿焊接方向拉长HAZ 晶粒尺寸用峰值温度决定的对数正态分布。热源方面常用 Goldak 双椭球热源写成 DFLUX 用户子程序配合温度相关的导热系数和比热做热力耦合分析。SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT, COORDS,JLTYP,TEMP,PRESS,SNAME) C 双椭球热源前半段热流密度按高斯分布衰减 INCLUDE ABA_PARAM.INC DIMENSION COORDS(3), FLUX(2), TIME(2) REAL Q, A, B, C, X0, Y0, Z0, V Q 0.85 * 20000.0 ! 热效率乘功率单位 W A 0.004 ! 熔池半宽m B 0.004 ! 熔池深度m C 0.006 ! 前半椭球长度m V 0.005 ! 焊接速度m/s X0 V * TIME(2) ! 热源中心沿焊缝移动 FLUX(1) (6.0*SQRT(3.0)*Q/(A*B*C*3.1415927*SQRT(3.1415927))) * EXP(-3.0*(COORDS(1)-X0)**2/A**2) * EXP(-3.0*(COORDS(2)-0.0)**2/B**2) * EXP(-3.0*(COORDS(3)-0.0)**2/C**2) RETURN END子程序里 FLUX(1) 是单位面积热流三个指数项分别控制沿焊接方向、板宽方向和深度方向的热流衰减。焊接速度 V 和热源功率 Q 是标定重点Q 偏大熔池过宽A、B、C 三轴参数要和焊缝截面金相照片对照校核。Voronoi 晶粒模型在这个阶段的作用是给 HAZ 提供温度相关的各向异性热膨胀和高温流变行为晶粒取向不同导致的热应力分布差异直接决定后续晶界开裂的起裂位置。4.2 晶界 cohesive零厚度单元插入顺序cohesive 和 voronoi 结合的场景主要是沿晶开裂晶界比晶粒内部弱裂纹沿 Voronoi 边扩展。做法是在共享边上插入零厚度 COH2D4。操作顺序固定先找出所有被两个晶粒共用的边然后把这条边上的节点复制一份原节点归左侧晶粒、新节点归右侧晶粒用两对节点构造一个 COH2D4 单元单元初始厚度为零由牵引-分离本构决定开裂行为。from collections import defaultdict def find_shared_edges(elems): elems: {grain_id: [(n1,n2,n3), ...]}返回被两个晶粒共用的边。 edge_owner defaultdict(set) for gid, tris in elems.items(): for tri in tris: for a, b in [(0, 1), (1, 2), (2, 0)]: edge_owner[frozenset((tri[a], tri[b]))].add(gid) return {e: sorted(gs) for e, gs in edge_owner.items() if len(gs) 2}返回字典的 key 是共享边两端节点value 是两个晶粒号。后续为新节点重新编号、构造 cohesive 单元时必须以 value 里的归属性决定节点归属否则会出现晶粒 A 的单元引用了属于晶粒 B 的节点拓扑直接乱掉。插入完成后原单元集合要整体重编号cohesive 单元单独放一个 *ELSET方便后处理里输出 SDEG、CSDMG 等损伤变量。注意共节点多的晶界三叉点三个晶粒交会处不能简单插入 cohesive需要特殊的节点分裂规则常见的工程做法是给三叉点留一小段无 cohesive 的钝化区避免单元退化。cohesive 材料参数是焊接沿晶开裂仿真的核心金属晶界典型量级参考下表具体值要按实验或文献标定。参数含义金属晶界典型量级Enn / Ett法向/切向牵引刚度1e5 – 1e7 N/mm³Tmax最大牵引应力100 – 1000 MPaGIC / GIIC法向/切向断裂能10 – 100 J/m²BK 指数混合模式耦合1.0 – 2.0对应字符块如下QUADS 准则控制损伤起始ENERGY 方式控制损伤演化MIXED MODEBK 处理法向和切向同时加载的情况。*MATERIAL, NAMEGB_COHESIVE *ELASTIC, TYPETRACTION 100000., 100000., 0. *DAMAGE INITIATION, CRITERIONQUADS 50., 80., 80. *DAMAGE EVOLUTION, TYPEENERGY, MIXED MODEBK, POWER1.5 10., 20., 0.牵引刚度 Enn 不是随便取的物理量它决定 cohesive 单元在开裂前贡献的额外柔度。取值过小整体模型刚度被晶界显著软化取值过大显式求解的稳定时间步急剧缩小。经验做法是把 Enn 设为晶粒弹性模量除以一个特征晶界宽度焊接近熔点温度附近材料刚度和强度都下降Tmax 和 GIC 要按高温段折减否则裂纹全部起于低温区结果与实测不符。4.3 周期边界与 RVE 约束的取舍RVE 的周期边界在 Abaqus 里用 *EQUATION 配对上相对面的节点自由度。Neper 输出的 inp 里周期节点列表现成可用Python 路线则在 3×3 平铺时就能记录边界配对关系写成下面的方程约束让对面节点位移差等于零宏观应变载荷通过额外参考点施加。*EQUATION 2 node_right, 1, 1.0, node_left, 1, -1.0每行一条方程前导数字 2 表示该方程涉及两个节点后面依次是节点号、自由度、系数。right 和 left 是对应面上的节点对这样的三维三个方向要写若干组。焊接仿真如果只取焊缝局部区域不建议加周期边界移动热源引起的温度梯度在边界上不满足周期性强制周期假设会人为约束热变形产生虚假应力。这种场景的常见做法是取足够大的板材模型远场边界用热-力对称条件焊缝区域局部细化网格。Voronoi 多晶模型加周期边界是均匀化分析的标准配置一旦叠加移动热源就要从 RVE 切回带真实热流边界的几何这两个模型不要混用。5. Abaqus 的 libpng error 与 Voronoi 模型交付前的三个验证5.1 libpng error 的触发场景与处理顺序Abaqus 的 libpng error 最常见于向 Sketch 导入 PNG 图片描 Voronoi 线框或者后处理导出位图时出现报错形如 libpng error: IDAT: invalid distance too far back。根源是 Abaqus 加载的 libpng 版本较旧读不了用高版本 zlib 压缩、16 位深度或带异常色块的 PNG。不要试图改 PNG 文件头去骗解析器最稳妥的做法是重编码# 8 位 RGB 重编码去除索引色和异常色块 python -c from PIL import Image; Image.open(voronoi.png).convert(RGB).save(voronoi_8bit.png) # ImageMagick 方案顺手把长边压到 4096 像素以内 magick voronoi.png -depth 8 -strip -resize 4096x4096 voronoi_fixed.png参数上 -depth 8 强制 8 位位深-strip 去掉 PNG 附带的颜色配置信息-resize 降低异常大图触发解压问题的概率。如果重编码后仍报错再排查运行时环境Linux 下检查 LD_LIBRARY_PATH 是否指向旧 libpng用 ldd 确认求解器进程实际加载的动态库版本Windows 下则是核对安装目录里的 PNG 解析动态库是否损坏。处理顺序永远是先重编码图片、再排查动态库直接替换系统动态库会影响其他组件不作为第一步推荐。5.2 晶粒尺寸分布与网格质量的双重自检模型几何对不对用统计量说话。生成完所有晶粒多边形后用等效直径拟合对数正态分布和设计值对比from scipy import stats def check_grain_size(polys): areas [abs(poly_area(p)) for p in polys] # 多边形面积 d_eq 2.0 * np.sqrt(np.array(areas) / np.pi) # 等效圆直径 mu, sigma stats.lognorm.fit(d_eq) print(f拟合期望{np.exp(mu):.2f}, 形状参数{sigma:.2f}) # 与设计参数偏差超过 15% 时回头调种子点采样参数拟合偏差超过 15% 说明种子点分布与目标组织不符常见原因是种子点数量太少导致统计波动大或者对数正态采样时直接用了正态分布。网格质量检查通常在提交前对模型做 aspect ratio 过滤ABAQUS 后处理里按单元质量指标着色重点看晶界两侧单元尺寸比超过 5 倍的区域要加密或调整种子。5.3 最小验收流程先跑 5 晶粒模型最省事的验收办法是拿 5 个晶粒的模型先跑一遍材料设成单晶弹性后处理里按晶粒显示应力云图确认每个晶粒的颜色对应你指定的欧拉角方向再给一个晶界组合插入单排 cohesive 单元做单轴拉伸看 SDEG 是否沿晶界连续演化。两步都过了这套 Voronoi 几何和 inp 组织方式才算合格再换回完整晶粒数。本文还有配套的精品资源点击获取 Abaqus-Voronoi 多晶建模：从几何生成到焊接仿真全流程解析
简介这是一份面向ABAQUS用户的Voronoi图生成代码资源适用于多晶体材料、颗粒复合材料、裂纹扩展与颗粒分布等有限元分析场景。Voronoi图本质上是基于种子点的空间划分又称泰森多边形可在ABAQUS中用于构建具有随机几何特征的非均质模型。压缩包共8个文件包含4个ipynb交互式脚本、3个txt说明文档和1个py核心脚本包体仅159KB轻量易用其中ipynb适合逐步演示生成流程py脚本可直接运行或修改txt文档用于记录依赖与参数说明。目前已有522人浏览学习。读者借助这些脚本可完成种子点定义、Voronoi区域计算、有限元网格划分、边界条件设置及后处理分析形成一套从几何生成到力学计算的完整链路。代码结构清晰便于有Python和ABAQUS脚本基础的用户直接复用、二次开发或嵌入自有模型对初学者而言也可通过跟随交互式脚本理解Voronoi图与有限元结合的基本原理是兼顾实用与学习价值的轻量工具包。1. 拿到 Abaqus-Voronoi 工程包之前先弄清楚这套东西解决什么问题拿到 Abaqus-Voronoi--master.zip 这类工程包最常见的动作是解压、找 .inp、提交作业然后被网格或材料报错劝退。Voronoi 剖分把随机分布的种子点变成凸多边形或多面体每个单元代指一个晶粒于是 Abaqus 里就能做晶粒尺度的力学、热学与断裂分析金属多晶各向异性、沿晶开裂、焊接热影响区组织都从这个几何基础出发。标题拆开看Abaqus 是求解器Voronoi 是晶粒几何算法master.zip 意味着一套可二次开发的模板工程。适合做过常规各向同性有限元、接下来要处理晶粒取向、晶界本构和微结构统计的工程师。下面按几何生成、网格与材料、晶界与边界、报错与验证的顺序走通链路每一步都给可复现的命令和参数。2. Voronoi 几何生成种子点规则决定晶粒分布两条路线对照2.1 剖分算法保证什么、不保证什么Voronoi 剖分的定义很简洁给定 n 个种子点空间里每个点被划分给离它最近的那个种子。二维情形下每个晶粒是一个凸多边形三维是凸多面体。它保证晶粒之间无缝隙、无重叠晶界由相邻种子点的中垂线构成几何上是完备的但不保证晶粒尺寸的均匀性。用完全均匀随机分布的种子点生成晶粒面积会呈近似指数分布小晶粒偏多这与真实金属里常见的对数正态晶粒尺寸分布相差很远。所以建模第一步不是写网格代码而是先设定目标晶粒尺寸分布。常见做法是先按对数正态分布采样等效直径 diameq再反推种子点密度和分布范围。遇到大晶粒包围小晶粒的组织可以在生成剖分后做一次 Lloyd 松弛迭代把每个种子移到自己单元的质心再重新剖分让结构趋近 centroidal Voronoi。这样晶粒尺寸更均匀、单元形状也更接近正多边形后续网格质量会明显改善Abaqus 里出现畸形单元的几率大幅下降。2.2 路线一Neper 一键生成带周期性的多晶网格如果目标是几百个晶粒以上的代表性体积单元 RVE我一般直接用开源工具 Neper而不是自己写剖分和网格代码。Neper 分两步先用 -T 模块生成剖分再用 -M 模块加密网格并输出 Abaqus inp。# 生成 200 个晶粒、等效直径服从对数正态分布的周期剖分 neper -T -n 200 -morpho diameq:lognormal(15,0.3) -periodic all -o poly # 在上述剖分基础上生成网格输出 Abaqus 输入文件 neper -M -n 200 -periodic all -format inp -o poly第一条命令里 -n 指定晶粒个数-morpho 里的 diameq:lognormal(15,0.3) 让等效直径的均值和形状参数按对数正态设定-periodic all 表示三个方向都做周期剖分。这是做 RVE 时最值钱的一个选项它保证相对面的网格节点一一对应后面加周期性位移边界条件时不用自己配对节点。第二条命令的 -format inp 直接把 mesh 写成 poly.inp里面已经包含节点、单元、晶粒 elset 和周期节点对应关系。Neper 输出的 inp 里单元按晶粒分好了 elset省去手工归类步骤。需要注意两个地方一是它默认生成四面体单元想换六面体要显式指定单元类型参数二是如果后续要插入 cohesive 晶界Neper 的共享面信息不会自动写进 inp需要在后处理时自己重建晶粒邻居关系这一点放到第四章展开。另外输出版本较老的 Neper 生成的 inp 里关键字缩进风格和 Abaqus 新版本有细微差异提交前用 abaqus datacheck 跑一遍最稳。2.3 路线二Python scipy.spatial.Voronoi 自己造几何不用 Neper、只想在现有 Python 流程里快速出模型时scipy.spatial.Voronoi 是最短路径。但它有两个坑边界上的单元会退化到无穷远必须裁剪完全随机的种子会让晶粒大小起伏过大。下面这段是带 3×3 镜像扩展的最小实现镜像的目的是消除非周期边界让中心区域的晶粒在边界处自然闭合。import numpy as np from scipy.spatial import Voronoi def voronoi_tiled(seeds, L(100.0, 100.0)): 3x3 平铺种子点做 Voronoi中心区域得到周期连续的多边形。 tiles np.array([(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)]) offset tiles * np.array(L) sup (seeds[None, :, :] offset[:, None, :]).reshape(-1, 2) return Voronoi(sup) rng np.random.default_rng(7) seeds rng.uniform(0, 100, size(80, 2)) vor voronoi_tiled(seeds) # vor.regions 中的多边形顶点在 vor.vertices 中 # 取所有顶点都落在 [0,100]x[0,100] 内的区域即为可用晶粒这段代码把种子点复制到周围 8 个相邻区块对 9 倍密度的点集做一次 Voronoi然后只保留落在中心区块内的有限多边形。做出来的晶界在边界处与对侧晶界自然连续等效实现了周期剖分。参数上seeds 的个数决定平均晶粒面积L 是 RVE 边长两者配合决定晶粒数密度如果做晶粒尺寸敏感性分析建议固定 L只改种子的对数正态采样参数否则几何变化和尺寸变化耦合在一起后处理里分不清是谁引起的力学响应差异。拿到多边形之后网格生成这一步 scipy 帮不上忙。常见做法是把多边形集合导出成 Gmsh 几何.geo或者直接在 Abaqus 的 Part 模块用 Partition Face 逐个分区再自由网格。我习惯对每个晶粒单独铺三角形再合并确保网格不跨晶界这个细节直接关系到第四章 cohesive 单元的生成。2.4 两条路线怎么选看晶粒数、周期性和后续修改频率Neper 和 Python 自建路线的差别主要在周期支持、网格生成和中间文件透明度。按我的使用经验对比没有绝对优劣看场景选即可。对比项NeperPython scipy周期剖分-periodic all 原生支持需 3×3 平铺自行实现晶粒尺寸分布diameq:lognormal(μ,σ) 直接指定需自己采样种子并做松弛网格划分内置网格器直接出 inp只出几何网格借 Gmsh 或 Abaqus二次开发CLI 参数控制脚本友好全 Python可嵌入现有数据流典型适用大型 RVE、晶体塑性参数反演教学演示、快速原型、对接实验数据选型时我的判断标准是晶粒数少于 300、并且最终要自己加 cohesive 或自己控制晶粒取样的场景Python 路线更顺手晶粒数上千、还要保证周期性和网格均匀性Neper 性价比高得多。还有一个容易被忽略的点Neper 输出的晶粒边界是理想化平面实验重构的 EBSD 数据往往是任意多面体如果你的模型要逐晶粒校准取向两条路线都得额外写数据映射层把 EBSD 取向表按晶粒号关联到 elset。3. 把 Voronoi 装进 AbaqusPart 拓扑、inp 拼装与晶粒取向3.1 先定模型拓扑单 Part 多 Section 和多 Part 装配的区别Voronoi 多晶模型在 Abaqus 里最常见的失败原因是拓扑选错。多 Part 装配方案里每个晶粒是独立 Part界面网格必然不共节点后期要靠 *TIE 绑定这会引入主从面选择和接触搜索的麻烦晶界 cohesive 也难做。我一般不用这个方案而是把整个 RVE 建成一个 Part所有晶粒共享连续网格晶粒之间的边界只是 elset 划分的产物。在 CAE 里实现这套拓扑的标准做法是先建一个矩形或立方体 Part用 Sketcher 导入 Voronoi 多边形线框再在 Part 模块用 Partition Face 按多边形逐个分区最后统一划分网格。网格跨晶粒自动共节点这是后续一切晶粒级分析的前提。划分时不要让单元穿过晶界2D 用 CPS3 或 CPS6M 三角形时边界天然贴合3D 建议用四面体并按晶粒设置局部种子密度界面两侧尺寸突变过大会在晶界处产生畸形单元求解时先报错的地方往往就在这里。3.2 从 Voronoi 网格直接写出可提交的 .inp不想在 CAE 里手工分区时把上一章 Python 生成的多边形直接写成 inp 是最可控的路径。下面这段把节点、单元、晶粒 elset 和取向一起写进文件以 2D 三角形单元 CPS3 为例。def write_voronoi_inp(path, polys): polys: [(verts_array, euler_angle), ...]输出单 Part 的 Abaqus inp。 node_id, elem_id 1, 1 lines, nid_map [*NODE], {} grain_elems, grain_angles [], [] for verts, angle in polys: ids [] for x, y in verts: key (round(float(x), 6), round(float(y), 6)) if key not in nid_map: nid_map[key] node_id lines.append(f{node_id}, {x:.6f}, {y:.6f}) node_id 1 ids.append(nid_map[key]) # 扇形式三角化以第一个顶点为公共点拆分凸多边形 tris [(elem_id k, ids[0], ids[k], ids[k 1]) for k in range(1, len(ids) - 1)] elem_id len(tris) grain_elems.append(tris) grain_angles.append(angle) lines.append(*ELEMENT, TYPECPS3, ELSETALLELEM) for tris in grain_elems: for eid, a, b, c in tris: lines.append(f{eid}, {a}, {b}, {c}) for gi, (tris, angle) in enumerate(zip(grain_elems, grain_angles), 1): lines.append(f*ELSET, ELSETGRAIN_{gi}) lines.append(, .join(str(t[0]) for t in tris)) lines.append(f*SOLID SECTION, ELSETGRAIN_{gi}, MATERIALCU_SINGLE) lines.append(f*ORIENTATION, NAMEORI_{gi}, SYSTEMRECTANGULAR) lines.append(f{angle:.4f}, 0.0, 0.0, 0.0, 1.0, 0.0) with open(path, w) as f: f.write(\n.join(lines) \n)写 inp 时最容易错的是 elset 与单元号的对应。这段脚本在三角化时直接按晶粒记录单元号列表再拼 *ELSET 行单元归属和几何生成在同一次循环里完成不会错位。节点去重用了六位小数的坐标作为字典键浮点噪声会被 round 吸收。*ORIENTATION 的两个方向点定义局部坐标系的 X 轴和 XY 平面晶粒的欧拉角通过这个局部坐标系施加到各向异性材料上如果做晶体塑性还要在这里配 *DEPVAR 和 UMAT先用单晶弹性矩阵验证取向是否生效再叠加晶体塑性本构排错面会小很多。3.3 晶粒取向、单晶弹性矩阵和 Schmidt 因子检查材料卡片上多晶模型的起点是单晶立方弹性。铜、铝这类面心立方金属用 *ELASTIC, TYPEANISO 时给 C11、C12、C44 三个独立常数即可。晶粒取向的正交性检查是新手最容易忽略的*ORIENTATION 定义的局部坐标必须保证三轴正交否则 Abaqus 会警告并强制正交化取向就偏离你输入的欧拉角。各晶粒的 inp 数据块关系可以用一张表理清。数据块作用常见错误*NODE节点坐标浮点重复导致共节点失败*ELEMENT单元连接晶界穿越、单元反转*ELSET晶粒归属单元号错位取向张冠李戴*ORIENTATION晶粒局部坐标三轴不正交被强制修正提交前用便宜的办法验证取向是否正确把材料暂时设成各向同性弹性给模型单向拉伸位移边界后处理里看每个晶粒的应力分布。因为各晶粒弹性刚度相同但方向不同晶粒内应力会呈现与取向相关的明暗差异。若所有晶粒应力场完全相同基本可以断定 *ORIENTATION 没被 *SOLID SECTION 引用到或者欧拉角转换时多了 90 度偏差。这种检查在焊接仿真里同样适用热力耦合计算前先跑一个纯力学小模型确认晶粒应力响应符合预期再叠加温度场排错成本会低很多。4. 焊接仿真里的 Abaqus-Voronoi 组合cohesive 晶界、热源与周期边界4.1 焊接热影响区组织演化的 Voronoi 近似abaqus 焊接仿真里Voronoi 最多被用在焊缝金属和热影响区 HAZ 的晶粒组织重建。真实焊缝凝固组织是柱状晶HAZ 是粗化的等轴晶两者可以用不同种子参数的两套 Voronoi 拼合近似焊缝区种子密度低、晶粒沿焊接方向拉长HAZ 晶粒尺寸用峰值温度决定的对数正态分布。热源方面常用 Goldak 双椭球热源写成 DFLUX 用户子程序配合温度相关的导热系数和比热做热力耦合分析。SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT, COORDS,JLTYP,TEMP,PRESS,SNAME) C 双椭球热源前半段热流密度按高斯分布衰减 INCLUDE ABA_PARAM.INC DIMENSION COORDS(3), FLUX(2), TIME(2) REAL Q, A, B, C, X0, Y0, Z0, V Q 0.85 * 20000.0 ! 热效率乘功率单位 W A 0.004 ! 熔池半宽m B 0.004 ! 熔池深度m C 0.006 ! 前半椭球长度m V 0.005 ! 焊接速度m/s X0 V * TIME(2) ! 热源中心沿焊缝移动 FLUX(1) (6.0*SQRT(3.0)*Q/(A*B*C*3.1415927*SQRT(3.1415927))) * EXP(-3.0*(COORDS(1)-X0)**2/A**2) * EXP(-3.0*(COORDS(2)-0.0)**2/B**2) * EXP(-3.0*(COORDS(3)-0.0)**2/C**2) RETURN END子程序里 FLUX(1) 是单位面积热流三个指数项分别控制沿焊接方向、板宽方向和深度方向的热流衰减。焊接速度 V 和热源功率 Q 是标定重点Q 偏大熔池过宽A、B、C 三轴参数要和焊缝截面金相照片对照校核。Voronoi 晶粒模型在这个阶段的作用是给 HAZ 提供温度相关的各向异性热膨胀和高温流变行为晶粒取向不同导致的热应力分布差异直接决定后续晶界开裂的起裂位置。4.2 晶界 cohesive零厚度单元插入顺序cohesive 和 voronoi 结合的场景主要是沿晶开裂晶界比晶粒内部弱裂纹沿 Voronoi 边扩展。做法是在共享边上插入零厚度 COH2D4。操作顺序固定先找出所有被两个晶粒共用的边然后把这条边上的节点复制一份原节点归左侧晶粒、新节点归右侧晶粒用两对节点构造一个 COH2D4 单元单元初始厚度为零由牵引-分离本构决定开裂行为。from collections import defaultdict def find_shared_edges(elems): elems: {grain_id: [(n1,n2,n3), ...]}返回被两个晶粒共用的边。 edge_owner defaultdict(set) for gid, tris in elems.items(): for tri in tris: for a, b in [(0, 1), (1, 2), (2, 0)]: edge_owner[frozenset((tri[a], tri[b]))].add(gid) return {e: sorted(gs) for e, gs in edge_owner.items() if len(gs) 2}返回字典的 key 是共享边两端节点value 是两个晶粒号。后续为新节点重新编号、构造 cohesive 单元时必须以 value 里的归属性决定节点归属否则会出现晶粒 A 的单元引用了属于晶粒 B 的节点拓扑直接乱掉。插入完成后原单元集合要整体重编号cohesive 单元单独放一个 *ELSET方便后处理里输出 SDEG、CSDMG 等损伤变量。注意共节点多的晶界三叉点三个晶粒交会处不能简单插入 cohesive需要特殊的节点分裂规则常见的工程做法是给三叉点留一小段无 cohesive 的钝化区避免单元退化。cohesive 材料参数是焊接沿晶开裂仿真的核心金属晶界典型量级参考下表具体值要按实验或文献标定。参数含义金属晶界典型量级Enn / Ett法向/切向牵引刚度1e5 – 1e7 N/mm³Tmax最大牵引应力100 – 1000 MPaGIC / GIIC法向/切向断裂能10 – 100 J/m²BK 指数混合模式耦合1.0 – 2.0对应字符块如下QUADS 准则控制损伤起始ENERGY 方式控制损伤演化MIXED MODEBK 处理法向和切向同时加载的情况。*MATERIAL, NAMEGB_COHESIVE *ELASTIC, TYPETRACTION 100000., 100000., 0. *DAMAGE INITIATION, CRITERIONQUADS 50., 80., 80. *DAMAGE EVOLUTION, TYPEENERGY, MIXED MODEBK, POWER1.5 10., 20., 0.牵引刚度 Enn 不是随便取的物理量它决定 cohesive 单元在开裂前贡献的额外柔度。取值过小整体模型刚度被晶界显著软化取值过大显式求解的稳定时间步急剧缩小。经验做法是把 Enn 设为晶粒弹性模量除以一个特征晶界宽度焊接近熔点温度附近材料刚度和强度都下降Tmax 和 GIC 要按高温段折减否则裂纹全部起于低温区结果与实测不符。4.3 周期边界与 RVE 约束的取舍RVE 的周期边界在 Abaqus 里用 *EQUATION 配对上相对面的节点自由度。Neper 输出的 inp 里周期节点列表现成可用Python 路线则在 3×3 平铺时就能记录边界配对关系写成下面的方程约束让对面节点位移差等于零宏观应变载荷通过额外参考点施加。*EQUATION 2 node_right, 1, 1.0, node_left, 1, -1.0每行一条方程前导数字 2 表示该方程涉及两个节点后面依次是节点号、自由度、系数。right 和 left 是对应面上的节点对这样的三维三个方向要写若干组。焊接仿真如果只取焊缝局部区域不建议加周期边界移动热源引起的温度梯度在边界上不满足周期性强制周期假设会人为约束热变形产生虚假应力。这种场景的常见做法是取足够大的板材模型远场边界用热-力对称条件焊缝区域局部细化网格。Voronoi 多晶模型加周期边界是均匀化分析的标准配置一旦叠加移动热源就要从 RVE 切回带真实热流边界的几何这两个模型不要混用。5. Abaqus 的 libpng error 与 Voronoi 模型交付前的三个验证5.1 libpng error 的触发场景与处理顺序Abaqus 的 libpng error 最常见于向 Sketch 导入 PNG 图片描 Voronoi 线框或者后处理导出位图时出现报错形如 libpng error: IDAT: invalid distance too far back。根源是 Abaqus 加载的 libpng 版本较旧读不了用高版本 zlib 压缩、16 位深度或带异常色块的 PNG。不要试图改 PNG 文件头去骗解析器最稳妥的做法是重编码# 8 位 RGB 重编码去除索引色和异常色块 python -c from PIL import Image; Image.open(voronoi.png).convert(RGB).save(voronoi_8bit.png) # ImageMagick 方案顺手把长边压到 4096 像素以内 magick voronoi.png -depth 8 -strip -resize 4096x4096 voronoi_fixed.png参数上 -depth 8 强制 8 位位深-strip 去掉 PNG 附带的颜色配置信息-resize 降低异常大图触发解压问题的概率。如果重编码后仍报错再排查运行时环境Linux 下检查 LD_LIBRARY_PATH 是否指向旧 libpng用 ldd 确认求解器进程实际加载的动态库版本Windows 下则是核对安装目录里的 PNG 解析动态库是否损坏。处理顺序永远是先重编码图片、再排查动态库直接替换系统动态库会影响其他组件不作为第一步推荐。5.2 晶粒尺寸分布与网格质量的双重自检模型几何对不对用统计量说话。生成完所有晶粒多边形后用等效直径拟合对数正态分布和设计值对比from scipy import stats def check_grain_size(polys): areas [abs(poly_area(p)) for p in polys] # 多边形面积 d_eq 2.0 * np.sqrt(np.array(areas) / np.pi) # 等效圆直径 mu, sigma stats.lognorm.fit(d_eq) print(f拟合期望{np.exp(mu):.2f}, 形状参数{sigma:.2f}) # 与设计参数偏差超过 15% 时回头调种子点采样参数拟合偏差超过 15% 说明种子点分布与目标组织不符常见原因是种子点数量太少导致统计波动大或者对数正态采样时直接用了正态分布。网格质量检查通常在提交前对模型做 aspect ratio 过滤ABAQUS 后处理里按单元质量指标着色重点看晶界两侧单元尺寸比超过 5 倍的区域要加密或调整种子。5.3 最小验收流程先跑 5 晶粒模型最省事的验收办法是拿 5 个晶粒的模型先跑一遍材料设成单晶弹性后处理里按晶粒显示应力云图确认每个晶粒的颜色对应你指定的欧拉角方向再给一个晶界组合插入单排 cohesive 单元做单轴拉伸看 SDEG 是否沿晶界连续演化。两步都过了这套 Voronoi 几何和 inp 组织方式才算合格再换回完整晶粒数。本文还有配套的精品资源点击获取