简介：这是面向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(x+L) = 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 ≠ 0，x 方向的两个对应面之间在 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, tol=1e-3, rpPrefix='RP_E'): pbc3d_kernel.buildPeriodicConstraints( instanceName, tol=tol, rpPrefix=rpPrefix )这部分代码的核心是分离 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, rtol=1e-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 是匹配容差，单位跟模型一致。如果模型单位是 mm，rtol 取 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，人工抽查几十对，能发现很多网格层面的隐藏问题。
容差参数的实际影响很大。同样是 rtol=1e-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 承载 ε11，RP_E22 承载 ε22，RP_E33 承载 ε33，RP_E12 承载 ε12，RP_E13 承载 ε13，RP_E23 承载 ε23。约束方程把节点相对位移和参考点位移连起来，最终只对参考点施加边界条件就能控制整个单胞的变形。
def build_equations(pairs, shift_map, Lvec, rp_prefix='RP_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(terms=tuple(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( u1=0.01, u2=0.0, u3=0.0)这个最小模型也是后续做等效性能提取的基底。得到边界节点的反力后，用体积平均应力除以施加的应变，就能得到刚度矩阵的一列。把六个应变分量分别算一遍，等效弹性矩阵就齐了。
再往后扩展，可以考虑把 GUI 对话框加上，把实例名、容差、参考点前缀都暴露到界面上，做成真正意义上的工具栏插件。也可以把温度场的周期性约束加进去，自由度 11 对应温度，方程形式完全一样，只是把位移换成温度。大变形超弹性材料的 PBC 是另一个话题，线性方程约束在小变形前提下成立，大变形要引入更复杂的位移梯度控制，这个方向值得单独写一套实现。
踩过网格不周期的坑之后，我养成了一个习惯：在划分网格阶段就把对面种子数拍好照片存进文档，每次做新单胞都先核对这一步。插件写得再好，网格不周期一样白费。希望这套实现思路能帮你在做 3D 周期单胞分析时少走一段弯路。
本文还有配套的精品资源，点击获取
