Abaqus 3D周期单胞周期性边界条件插件开发实战

简介这是面向Abaqus CAE 2018用户的周期边界条件定义插件基于Python开发主要解决三维几何模型中周期性边界的快速施加问题适合从事复合材料、超材料及晶胞类结构有限元分析的工程师与研究人员。通过该插件可在两个网格化曲面间自动添加约束并支持对特殊节点或区域实现豁免避免手工逐个节点操作。资源包共13个文件、约2.04MB包含三个Py脚本、一个INP示例模型、一个CAE模型、使用说明PDF、界面预览图及LICENSE文件其中脚本分模块控制插件交互与核心边界计算便于使用者二次开发或迁移到其他版本。配套的CAE与INP文件可直接导入查看帮助理解周期性边界条件在实际模型中的设定方式。该资源已有886人学习适合希望提高Abaqus建模效率、需要结合晶胞模型开展周期性分析的仿真用户也可作为Abqus插件开发与Python脚本编写的参考案例。1. 为什么 3D 周期单胞的边界条件这么难写以及这个插件到底帮你做了什么做复合材料、多孔材料或超材料的等效性能分析时几乎绕不开周期性边界条件取一个代表体积元跑一次静力分析通过均匀化得到刚度矩阵。Abaqus CAE 2018 里没有现成的“施加周期性边界条件”按钮很多人就自己写 Python 脚本去生成约束方程。问题在于 3D 单胞和 2D 完全不同三对边界面上各有几百个节点棱上的节点同时属于两个面角点同时属于三个面稍不注意就会把同一个节点写进多个方程直接过约束。这个标题描述的就是这类插件用 Python 在 Abaqus CAE 2018 里读取三维网格自动完成边界节点的周期配对生成周期性位移约束方程写回当前模型。适合做线性小变形均质化分析的工程师尤其是用 C3D8、C3D20 这类三维实体单元建模、又被边界条件折腾过的人。它也值得放进你自己的脚本库里——下一批单胞分析能省掉一半体力活。2. 周期边界约束的数学与软件实现先厘清方程再动代码2.1 单胞位移场的数学形式u(xL) u(x) ε̄·L周期性边界条件的本质是把单胞看成无限周期结构中的一个胞元相邻胞元共享同一边界面位移连续、应力连续。用数学语言写对于沿第 j 个方向的周期长度 L_j两个周期等价点 x 和 x L_j·e_j 之间的位移满足u_i(x L_j·e_j) - u_i(x) ε̄_ij·L_j这里的 ε̄_ij 是宏观均匀应变张量下标 i 是位移分量方向j 是周期方向。注意右边是位移差不是零。很多人一开始理所当然地认为“对应面位移相同”那样只适用于无剪切变形的特殊情况。一旦宏观应变里有剪切分量比如 ε̄_12 ≠ 0x 方向的两个对应面之间在 U2 方向上必须存在一个固定的位移差否则单胞之间会撕开或者重叠。这个方程要分别对三个位移分量 U1、U2、U3 写所以每个节点对要生成三条独立的标量方程。3D 单胞里有三对边界面对应三个平移方向但棱和角点上的节点需要多个方向的叠加补偿。比如一条沿 x 方向的棱边它既是 y 方向面的边界又是 z 方向面的边界它的位移差要同时满足 y 方向和 z 方向两个周期条件的叠加。对称边界条件看起来简单很多人想用“约束面法向位移为零”代替 PBC这是个明显的坑。对称边界强制边界保持平面等价于对单胞施加了额外的刚度约束剪切变形被完全锁死算出来的等效弹性模量会偏高。PBC 允许边界自由翘曲只约束相对位移这才是均匀化分析的正确边界条件。2.2 在 Abaqus 里选 *EQUATION 而不是 *MPC约束格式与自由度Abaqus 实现这类线性多节点约束最常见的是 *EQUATION 和 *MPC。用户手册里对两者的定位很清楚MPC 是主从节点形式的线性约束适合梁、壳的刚性连接*EQUATION 是通用线性约束方程允许自由指定多个节点和自由度的加权和等于零。PBC 的形式是“主节点位移 - 从节点位移 - 周期位移 0”属于典型的线性约束用 *EQUATION 最直接。*EQUATION 在输入文件里的格式是*EQUATION 3 node_m, 1, 1.0, node_s, 1, -1.0, node_rp, 1, -Lx第一行的 3 表示方程项数每行按“节点编号、自由度、系数”为一组。Abaqus 把所有这些项的加权和等于零所以如果想让 U_m - U_s - Lx·ε 0 成立第三个项的系数就是 -Lx并且需要有一个承载 ε 的自由度。这就是参考点驱动法的由来为每个宏观应变分量建立一个参考点自由度边界条件最终只施加在这些参考点上。在 CAE 里用 Python API 创建约束不写输入文件用的是 Model.Equation 方法。每个方程是一个 terms 元组每一项是 (nodeLabel, dof, coefficient) 三元组。这个 API 的好处是约束创建后可以在 Interaction 模块的 Equation Constraints 管理器里直接看到和修改排查问题比对着 inp 文件方便得多。2.3 网格必须周期对应划分阶段就要满足的三个条件代码写得再好网格不是周期的就白搭。周期性边界条件要求每个边界节点在对面都有一个“完全对应的伙伴节点”节点坐标允许有数值误差但必须在容差范围内精确匹配。如果对面网格种子数不一致或者节点位置错位脚本就找不到配对数。网格周期化要在划分阶段解决不能靠约束脚本补救。常见做法是三条第一用 Partition 把三个方向的对应面切割成相同的拓扑结构保证面分区一一对应第二在对边上设置相同的种子数用 structured 或 sweep 算法划分第三检查节点坐标精度导入的网格经常有 1e-6 到 1e-4 量级的坐标误差匹配容差至少要覆盖到这个量级。还有一类专门针对二阶单元的坑。C3D20 这类二次六面体单元的边中点在边界面上也是一等公民节点必须参与配对。只看单元面会漏掉它们正确的做法是直接从网格节点集合里取全体节点坐标按坐标找周期等价关系而不是从几何面集合里选点。后面第 4 章的算法就是按坐标匹配设计的能同时覆盖角点、棱点和面内节点。3. 插件架构与数据准备从 CAE 里把边界节点找出来3.1 插件的目录结构与注册方式一个能在 Abaqus CAE 2018 里用的插件通常是一个写在 abaqus_plugins 目录下的 Python 包。常见结构是一个主插件脚本负责注册菜单按钮一个 kernel 脚本放真正的计算逻辑。kernel 脚本不依赖 GUI可以直接用 File → Run Script 执行这个特性在调试阶段非常有用。# pbc3d_plugin.py # -*- coding: utf-8 -*- # 在 Abaqus CAE 2018 (Python 2.7) 下运行 from abaqus import * from abaqusConstants import * import pbc3d_kernel # 插件入口函数纯计算逻辑交给 kernel def createPBCForInstance(instanceName, tol1e-3, rpPrefixRP_E): pbc3d_kernel.buildPeriodicConstraints( instanceName, toltol, rpPrefixrpPrefix )这部分代码的核心是分离 GUI 和计算逻辑。我把插件入口定义成普通函数而不是类因为调试时可以直接在 CAE 的命令行里调用createPBCForInstance(PartInstance-1)看到结果之后再考虑注册菜单按钮。把容差和参考点前缀作为参数暴露出来是为了方便在命令行里快速试不同容差。菜单注册不是重点RSG 对话框和 GUI Toolkit 的细节随版本变化很大工作中我一般先做纯 kernel 版本验证算法没问题后再套 GUI 外壳。Abaqus 2018 的内置 Python 是 2.7这一点要特别注意文件开头写# -*- coding: utf-8 -*-提示信息尽量用英文避免中文字符串在打印时报 Unicode 错误。3.2 从装配体中提取网格节点与几何包围盒周期配对的第一步是把装配体里某个 Part Instance 的全部节点坐标取出来。Abaqus Python API 里MeshNode 对象的 coordinates 属性是三维元组label 是节点编号。注意要用装配体实例而不是原始 Part 的节点因为装配体的坐标才是最终计算用的坐标。# pbc3d_kernel.py 片段 from abaqus import mdb from abaqusConstants import * import math def _get_instance_nodes(instanceName): # 从装配体里取出指定的实例 inst mdb.models[Model-1].rootAssembly.instances[instanceName] nodes inst.nodes # 组装 (label, 坐标) 列表 node_list [(n.label, n.coordinates[0], n.coordinates[1], n.coordinates[2]) for n in nodes] return node_list这段代码的逻辑是先把节点 label 和坐标打包成列表后面所有处理都基于这个列表避免反复调用 API 查询。这里有两个容易忽略的点第一nodes是 MeshNodeArray直接迭代即可但单个节点的 coordinates 是浮点元组不是 numpy 数组第二如果实例被移动过坐标已经是装配坐标后面求最小值和周期长度直接对这个坐标做运算。拿到节点列表后下一步是算包围盒。遍历所有坐标分别求三个方向的最小值和最大值周期长度就是最大值减最小值。我一般会打印出来核对一下这个值直接决定方程里的系数错了整套约束都会偏。def _get_bounding(node_list): xs [c[1] for c in node_list] ys [c[2] for c in node_list] zs [c[3] for c in node_list] min_x, max_x min(xs), max(xs) min_y, max_y min(ys), max(ys) min_z, max_z min(zs), max(zs) Lx max_x - min_x Ly max_y - min_y Lz max_z - min_z return (min_x, min_y, min_z), (Lx, Ly, Lz)参数说明返回值里 min_x 是整体坐标的偏移量后续取模运算必须先减去它否则坐标从 100 开始算模 L 的结果全是碎的根本没法匹配。3.3 周期等价类的建立按坐标模长分桶配对算法最稳的思路不是按“哪两个面配对”来写而是把所有边界节点投影到局部周期坐标里把“坐标对周期长度取模”后相等的节点归为同一个周期等价类。这样面节点、棱节点、角点统一处理从数据结构上避免了重复约束。def classify_periodic(node_list, origin, Lvec, rtol1e-3): ox, oy, oz origin Lx, Ly, Lz Lvec groups {} for label, x, y, z in node_list: # 平移到单胞局部坐标后再取模 nx (x - ox) % Lx ny (y - oy) % Ly nz (z - oz) % Lz # 量化到整数区间消除浮点误差 key (round(nx / rtol), round(ny / rtol), round(nz / rtol)) groups.setdefault(key, []).append(label) return groups原理说明%运算把坐标折叠到 [0, L) 区间周期等价点的局部坐标理论上完全相同但浮点误差会导致 3.999999 和 4.000001 对不上所以再用 rtol 量化。rtol 是匹配容差单位跟模型一致。如果模型单位是 mmrtol 取 1e-3 意味着匹配精度是微米级一般够用如果网络是从其他软件导入的先查一下坐标打印精度再调这个值。这一步跑完groups 的每个 value 就是一个等价节点集合。理想单胞里每个集合应该有两个节点一个在 min 面一个在 max 面。只有一个节点的情况通常是三周期方向上自等价点比如正中心节点需要单独标记。超过两个节点的情况多半是模型不是单胞或者坐标有严重偏离此时应该停下来查网格而不是继续生成约束。4. 核心实现节点配对、棱角去重与 *EQUATION 生成4.1 周期等价类一对一的配对算法与容差设置分桶之后的配对变得非常简单每个 group 内部就是一个周期等价类选其中一个节点作 master其余作 slave。3D 单胞的边界节点理论上每个等价类只有两个成员所以配对就是二选一。我习惯选 label 较小的那个作为 master原因纯粹是排查约束时看 inp 文件更容易定位。def build_pairs(groups): pairs [] # (master, slave) 节点编号 singles [] # 只有自己的等价类需要告警 overlapped [] # 多于两个节点的等价类 for key, labels in groups.items(): if len(labels) 1: singles.append(labels[0]) elif len(labels) 2: m min(labels) s max(labels) pairs.append((m, s)) else: overlapped.append(labels) return pairs, singles, overlapped逻辑说明这个函数把上一步的分组结果转成后续生成方程用的 master-slave 对。singles 要排查overlapped 要直接报错原因留到第 5 章的避坑部分细说。配对阶段我还建议把结果写成一个 CSV 文件备份记录了 master、slave 和分组 key人工抽查几十对能发现很多网格层面的隐藏问题。容差参数的实际影响很大。同样是 rtol1e-3在 mm 单位模型里对应微米级偏差在 m 单位模型里对应毫米级偏差后者可能把两个本该不同的节点错误合并。所以我的习惯是先看模型单位再定 rtol计算完成后在日志里输出配对总数如果配对数量远小于理论边界节点数的一半就要回头查容差。4.2 棱节点和角节点的去重让每个节点只出现在一个方程里很多人写 PBC 脚本的翻车点是按面处理先找 x 方向的一对面生成一批约束再找 y 方向再找 z 方向。这样一条棱上的节点会同时出现在两个方向的循环里角点出现在三个方向的循环里最终同一对节点被反复约束两到三次。Abaqus 求解器遇到冗余约束不会直接拒绝但会在 Data check 阶段报一堆节点自由度奇异结果完全没法用。我的做法是彻底绕开“按面循环”只依赖周期等价类。每个节点对在等价类里天然只有一个映射从数据结构上就不会出现重复约束。生成方程时还要计算这个节点对跨了几个周期方向也就是平移向量。这个向量决定方程右侧要叠加几个参考点位移项。def _shift_vector(node_m, node_s, origin, Lvec): # 返回三个方向的整数平移分量取值通常为 -1, 0, 1 ox, oy, oz origin Lx, Ly, Lz Lvec dx round((node_m[1] - node_s[1]) / Lx) dy round((node_m[2] - node_s[2]) / Ly) dz round((node_m[3] - node_s[3]) / Lz) return (int(dx), int(dy), int(dz))这里 node_m 和 node_s 是包含 label 和坐标的元组。dx 为 0 表示该节点对在 x 方向上没有周期跨越dx 为 1 或 -1 表示跨了一个周期。方程右侧叠加时每个有值的分量都要加上对应的 L_j·ε̄_ij 项。角点的 dx、dy、dz 全部非零这是它比棱点复杂的地方但代码上只是同一个循环多几个分支。去重还有一个细节同一个节点如果同时是某条棱的 master 和另一个等价类的 slave会不会冲突不会。因为 PBC 约束的是“周期等价节点之间的相对位移”一个节点只需和一个伙伴建立方程它在哪个 group 里由局部坐标决定不会同时作为 master 出现在两个 group 里。等价类分组天然保证了一对一关系这就是我坚持用坐标取模而不是面集合的原因。4.3 用参考点承载应变分量生成方程的完整代码现在到了核心输出环节。我给每个宏观应变分量准备一个参考点自由度RP_E11 承载 ε11RP_E22 承载 ε22RP_E33 承载 ε33RP_E12 承载 ε12RP_E13 承载 ε13RP_E23 承载 ε23。约束方程把节点相对位移和参考点位移连起来最终只对参考点施加边界条件就能控制整个单胞的变形。def build_equations(pairs, shift_map, Lvec, rp_prefixRP_E): Lx, Ly, Lz Lvec equations [] # 每个节点对生成三个方向的位移方程 for m_label, s_label in pairs: dx, dy, dz shift_map[(m_label, s_label)] # 自由度 1,2,3 对应 U1, U2, U3 for dof in (1, 2, 3): terms [(m_label, dof, 1.0), (s_label, dof, -1.0)] # x 方向周期贡献 if dx ! 0: # 对应 ε11, ε12, ε13 由自由度 1,2,3 承担 rp_label _get_rp_label(rp_prefix, dof, 1) terms.append((rp_label, 1, -dx * Lx)) # y 方向周期贡献 if dy ! 0: rp_label _get_rp_label(rp_prefix, dof, 2) terms.append((rp_label, 1, -dy * Ly)) # z 方向周期贡献 if dz ! 0: rp_label _get_rp_label(rp_prefix, dof, 3) terms.append((rp_label, 1, -dz * Lz)) equations.append(terms) return equations参数说明系数 -dx·Lx 里的负号对应方程U_m - U_s - L·ε 0。dx 的正负由 master 和 slave 的位置决定如果 master 在 min 面、slave 在 max 面dx 为 -1-dx·Lx 为正方程右侧等价的参考点位移取正值。映射函数 _get_rp_label 负责把 (dof, 方向) 组合到正确的参考点自由度 1 配正向 ε11自由度 2 配 ε12自由度 3 配 ε13以此类推。完整的写回模型动作是调用m.Equation(termstuple(terms))。每个 terms 元组的第一个元素是主项系数为正后面的都是带负号的从项。方程创建完成后在 CAE 的 Interaction 模块里能看到一排 Equation constraints名字类似 Eqn-1、Eqn-2 这样点开任意一个都能看到具体的节点和系数这一步建议抽查两三个。参考点本身要预先创建在装配体的几何中心或任意位置数量固定为 6 个名字由 rp_prefix 控制。创建参考点后要对它们施加边界条件否则方程里引用了它们的自由度却没有给它约束求解器依然会报奇异。通常的做法是目标应变分量对应的参考点施加指定的位移值其余参考点自由度固定为零。5. 排错与避坑跑通这个插件前先看这 5 个高频翻车点5.1 过约束Data check 报 zero pivot 或 redundant constraints现象提交 Job 后 Data check 失败报错信息里有大量节点出现 zero pivot、negative eigenvalue或者直接提示 redundant constraints。模型明明很常规网格也没有问题。原因同一对节点被写进了多个约束方程。最常见的来源是按“面循环”生成约束棱节点和角节点被重复处理还有一个隐蔽来源是等价类分组时容差太小导致同一节点被分到两个不同的 key产生了两组 master-slave 关系。解决改用周期等价类分组保证每个节点只出现在一个等价类里代码里加一个断言检查 pairs 列表里 master 和 slave 是否有重复出现重复立刻中止。这个断言值得留成永久检查因为网格一变边界情况就可能复现。5.2 节点配对为空容差、网格周期性与坐标精度现象运行脚本后日志打印 “Found 0 periodic pairs”或者配对数量远小于理论边界节点数的一半。偶发现象是部分面配上了部分面一条都没配上生成的约束少得可怜。原因三选一。网格本身不周期对面种子数不一致匹配容差 rtol 设得比网格坐标误差还小导入的网格坐标被舍入比如从 inp 读取后只有 6 位小数0.000001 级偏差在取模后被放大。解决先在 CAE 里用 Query 功能检查对应面的节点坐标肉眼对比一下是否一一对应再打印局部坐标的直方图看看 key 的分桶是不是都集中在整数附近。容差从 1e-3 起步确认模型单位后逐步调大但不要超过网格尺寸的十分之一否则会误并节点。5.3 剪切应变结果差一倍工程应变与张量应变的系数现象施加 ε12 0.01 的剪切变形后后处理里看到的剪应力数值和手算差了一倍或者等效刚度矩阵里剪切分量明显不对。原因工程剪应变 γ12 和张量剪应变 ε12 之间差一个 2 的因子。方程里写 Lx·ε12 时用的是张量应变但边界条件给参考点赋位移时如果按工程应变的习惯给了 γ12 0.01实际输入的 ε12 就变成 0.005最后应力自然差一倍。解决统一约定。我习惯在插件文档里明确写“参考点位移值等于张量应变分量”并且生成一个 README 式的注释块把六个参考点和六个应变分量的对应关系列成表每次建模前看一眼。如果喜欢工程应变那就把方程系数整体乘 2但没必要张量应变才是 Abaqus 内部使用的定义。5.4 print 语法报错Abaqus 2018 内嵌 Python 2.7现象脚本运行到某一行直接语法错误报错指向 print 语句或者提示某个 Python 3 的语法不支持。原因Abaqus CAE 2018 内置解释器是 Python 2.7不是外界常见的 Python 3。很多人从 2019 版本转过来或者平时写 Python 3 写习惯了把 f-string、print(xxx) 这类语法直接带进插件。解决文件头加上from __future__ import print_function把 print 统一成函数形式字符串格式化用 format 方法不用 f-string。另一个习惯是把调试用的 print 集中到一个函数里替换成日志输出这样将来迁移到更高版本的 Abaqus 时改动集中在一处。5.5 刚体位移奇异PBC 只约束相对位移不约束刚体模态现象求解器报大量 DOF singularity或者位移云图出现整体平移量级的巨大数值应力却正常。有时候加上一个看似多余的边界条件后一切恢复正常让人百思不解。原因周期性边界条件约束的是节点之间的相对位移对整个单胞的刚体平动和刚体转动没有任何约束。方程再完备刚度矩阵里依然存在零特征值求解器自然报奇异。解决加一个最小约束把某个角点或者参考点的三个平动自由度固定。注意不要把多个节点同时固定否则会引入额外的内部约束破坏 PBC 的本意。我一般直接固定第一个参考点的 U1、U2、U3反正各应变分量都由参考点控制在参考点上加约束最容易被检查出来。6. 验证从最小单胞开始三分钟检查约束是否正确再谈扩展拿到插件生成的一堆约束后不要直接跑复合材料模型。先做一个最小验证创建一个 1×1×1 的立方体单胞用 C3D8 单元划分网格保证三对边界的种子数一致然后调用插件生成约束。这一步能在一分钟内暴露 90% 的网格和配对问题。验证分三步走。第一步施加应变 ε11 0.01把其他参考点固定为零提交分析第二步在后处理里检查 x 方向两个对应面上节点的 U1 差值理论上应该等于 0.01第三步检查应力云图是否均匀单胞内部的应力分布应该在边界附近没有任何突变。如果剪切应变方向不对反过来检查参考点编号映射。# 验证用例的参数设置 # 参考点边界条件U_RP1 0.01其余全部固定为 0 step mdb.models[Model-1].steps[Step-1] mdb.models[Model-1].boundaryConditions[BC_RP1].setValues( u10.01, u20.0, u30.0)这个最小模型也是后续做等效性能提取的基底。得到边界节点的反力后用体积平均应力除以施加的应变就能得到刚度矩阵的一列。把六个应变分量分别算一遍等效弹性矩阵就齐了。再往后扩展可以考虑把 GUI 对话框加上把实例名、容差、参考点前缀都暴露到界面上做成真正意义上的工具栏插件。也可以把温度场的周期性约束加进去自由度 11 对应温度方程形式完全一样只是把位移换成温度。大变形超弹性材料的 PBC 是另一个话题线性方程约束在小变形前提下成立大变形要引入更复杂的位移梯度控制这个方向值得单独写一套实现。踩过网格不周期的坑之后我养成了一个习惯在划分网格阶段就把对面种子数拍好照片存进文档每次做新单胞都先核对这一步。插件写得再好网格不周期一样白费。希望这套实现思路能帮你在做 3D 周期单胞分析时少走一段弯路。本文还有配套的精品资源点击获取