最近做一批厚壁筒和隧道围岩的算例时，我又遇到了那个绕不开的需求：把Abaqus结果里的积分点应力取出来，换算成径向应力，再和对应的位移一起沿着半径方向画成曲线。问这个问题的后台消息也挺多，简单回几句讲不清楚，干脆写一篇把这套流程完整拆开。本文围绕积分点的读写、从全局直角坐标到柱坐标的应力张量旋转变换、节点位移与积分点位移的区别、以及最后如何导出成工程可用的数据文件这几个环节展开，目标是用Python脚本把"提取积分点径向应力与位移"这件事做成一条能直接复用的流水线，减少来回翻后处理界面、手工抄数据的时间。
先给结论：如果模型是轴对称单元（CAX系列），提取径向应力可以直接读S11，径向位移直接读U1，几乎不用动脑子；但如果你建的是三维实体模型或平面模型，默认输出场在全局直角坐标系下，这时必须自己做张量旋转，否则拿到的"径向应力"是错的。这两种情况我都会讲到，并且给出核心代码和一套能直接跑的完整脚本，涵盖从ODB文件读取、单元信息判断到CSV落盘的全过程，适合刚接触Abaqus二次开发的人直接抄作业，也适合已经写过后处理脚本但被坐标变换、单位问题坑过的人对照排错。
1. 为什么默认后处理拿不到"径向应力"：先弄清应力存在哪里
1.1 应力场输出规则：积分点才是数据的老家
Abaqus在求解完成后，默认会把一系列场变量写入ODB文件。其中应力张量S、应变E这类跟本构积分有关的量，绝大多数情况下是存在**单元积分点（integration point）**上的，而不是节点上。也就是说，你在后处理界面看到的应力云图，实际上是程序用积分点上的数值经过外推（extrapolation）和平均后，映射到节点再渲染出来的"光滑结果"。云图颜色很直观，但它不是原始数据。
脚本读取则不同，我们能直接接触到最原始的一层。通过fieldOutputs['S'].values遍历到的每个value对象，携带了elementLabel（单元编号）、integrationPoint（积分点序号）、data（应力分量元组）、position（位置类型）等属性。当position显示为INTEGRATION_POINT时，读取到的就是求解器真实计算并存储的积分点应力，没有经过任何后处理平均，是最"干净"的数据源。
所以第一条经验：想要提取积分点径向应力，读的必须是积分点场输出；如果你不小心读的是节点场（比如某些只有节点输出的自定义变量），那已经是加工过的数据了。查看Pythonscript里打印val.position，养成这个习惯能帮你规避很多莫名其妙的数值偏差。
1.2 默认坐标系与"径向"的矛盾
Abaqus默认输出的应力张量，是针对**全局直角坐标系（X,Y,Z）**的六个分量，三维模型对应S11、S22、S33、S12、S13、S23，平面模型对应S11、S22、S33、S12。这个坐标系在建模过程里一般不会特意调整，圆柱孔的中心落在哪里，坐标原点就在哪里。
问题随之而来：所谓"径向应力"，在柱坐标系下的定义是沿着半径方向的正应力分量σ_r。直角坐标系的σ_x、σ_y，只是全局轴向的正应力，跟"径向"没有直接对应关系。打个比方，你站在圆柱外壁某一点，你面朝圆心方向的"挤压"程度，才是径向应力；而默认输出的是"沿着全局X轴的拉压"和"沿着全局Y轴的拉压"，这两个量只有在恰好位于坐标轴上的点时才有明确的物理解释。只要待提取的位置不在坐标轴上，就必须做坐标旋转。
这是整篇内容的核心痛点，也是为什么很多人在后处理里尝试PlotOnPath画路径曲线时，发现"沿径向的应力曲线"要么不存在、要么需要自己定义坐标系却不知道对不对的原因。默认后处理无法直接给你一条正确的σ_r沿半径分布曲线，除非你用脚本拿到原始积分点分量再自己算。
1.3 "提取"真正要产出的东西是什么
工程里最常见的需求是把结果整理成一条曲线或一张表格。比如厚壁筒内压分析，想看到沿壁厚方向径向应力从内壁的-10MPa逐渐过渡到外壁的0MPa，同时径向位移从内壁的最大值衰减到外壁的较小值。这时候你要输出的就不是单一应力值，而是一组对应关系：
- 每个积分点或节点的位置（用半径r表示）
- 该位置上的径向应力σ_r
- 该位置上的径向位移u_r
所以脚本的工作流可以概括为三步：先拿到积分点/节点的原始应力与位移分量；再根据该位置相对轴心的方位角，做坐标变换；最后把结果按半径排序并导出。搞清楚这个流程，后面的代码就只是具体执行细节而已。
2. 从ODB读数据的第一条路：轴对称模型直接读分量
2.1 轴对称单元的输出约定：S11就是径向
Abaqus对轴对称实体单元（CAX4R、CAX6M等）有一套固定的局部方向约定：1方向是径向（r），2方向是轴向（z），3方向是环向（θ）。这意味着结果文件里：
- 应力张量S的分量中，S11就是径向应力σ_r，S22是轴向应力σ_z，S33是环向应力σ_θ，S12是径向-轴向剪切应力τ_rz。
- 位移场U中，U1就是径向位移u_r，U2是轴向位移u_z。
这个细节很多教程不会单独拎出来讲，但它极其关键。也就是说，轴对称模型的"径向应力提取"方案非常简单：写出一个小脚本遍历积分点S的data[0]即可，不需要任何应力张量旋转。
我猜很多人的模型其实就是用轴对称单元建的，结果看到三维模型的变换公式以为很复杂，反而绕了远路。所以脚本里第一步要做的，不是闷头写旋转矩阵，而是先判断模型到底是不是CAX单元。
2.2 判断单元类型并提取分量的最小脚本
下面这段是我常用的"开局三板斧"：打开ODB、定位帧、判断单元类型。
# -*- coding: utf-8 -*- from odbAccess import openOdb odb = openOdb('C:/sim/cylinder.odb', readOnly=True) step = odb.steps['Step-1'] frame = step.frames[-1] # 最后一帧 # 遍历装配体中的实例，看第一个单元类型 for inst in odb.rootAssembly.instances.values(): etypes = {e.type for e in inst.elements} if all(t.startswith('CAX') for t in etypes): print(f'实例 {inst.name} 是轴对称模型，S11->径向应力，U1->径向位移') # 直接读应力场 stress_field = frame.fieldOutputs['S'] for val in stress_field.values: if val.instance.name == inst.name: sr = val.data[0] # 径向应力 sz = val.data[1] # 轴向应力 st = val.data[2] # 环向应力 srz = val.data[3] # 剪切应力 # 这里可以继续保存到列表注意val.data在不同模型类型下长度不一样：轴对称和平面模型通常有4个分量（S11、S22、S33、S12），三维模型有6个分量。索引别写错了。
2.3 积分点位置没有现成值时，用平均节点坐标近似
读到了应力，还得知道这个积分点到底在哪。麻烦的是，积分点坐标并不是S场输出里的标准属性，val.data只有应力分量，没有XYZ坐标。一个最简单的可行方案是用该单元所有节点的平均坐标来代表积分点位置，对线性单元来说误差很小，因为一阶单元的积分点基本位于单元中心附近。
nodes = inst.getElementFromLabel(val.elementLabel).connectivity coords = [inst.getNodeFromLabel(n).coordinates for n in nodes] x_avg = sum(c[0] for c in coords) / len(coords) y_avg = sum(c[1] for c in coords) / len(coords) # 轴对称模型里，坐标的1方向就是半径方向，2方向是轴向 radius = x_avg如果精度要求更高，可以用单元形函数在积分点局部坐标处做等参插值，得到精确积分点坐标。线性四边形单元（quad4）的积分点在局部坐标(±1/√3, ±1/√3)，用四节点双线性形函数插值即可。但对大多数径向应力统计场景，平均节点坐标已经够用，不会影响曲线的趋势判断。另外，Abaqus还支持在Step模块里额外输出COORD场，这样ODB中会直接存放积分点坐标，脚本里可以优先检查frame.fieldOutputs.get('COORD')，有就直接取val.data，没有再用上面的近似。工程上我建议建模时顺手勾上COORD输出，省去后面所有估算。
3. 三维模型的正经做法：旋转应力张量到柱坐标系
3.1 旋转矩阵推导：不要背公式，理解一次就够
三维模型默认输出的是全局直角坐标下的应力分量，设一点的坐标为(x, y, z)，该点相对柱轴（假设沿Z方向）的方位角为θ = atan2(y, x)。径向方向单位向量是(cosθ, sinθ, 0)，环向方向单位向量是(-sinθ, cosθ, 0)。应力张量要发生坐标旋转，本质上就是把局部的两个正交基向量和全局基向量的方向余弦关系代入二阶张量变换式。
落到具体分量上，绕Z轴旋转角度θ后，有关系式：
- σ_r = σ_x·cos²θ + σ_y·sin²θ + 2·τ_xy·sinθ·cosθ
- σ_θ = σ_x·sin²θ + σ_y·cos²θ - 2·τ_xy·sinθ·cosθ
- τ_rθ = (σ_y - σ_x)·sinθ·cosθ + τ_xy·(cos²θ - sin²θ)
- σ_z 保持不变
- τ_rz = τ_xz·cosθ + τ_yz·sinθ
- τ_θz = -τ_xz·sinθ + τ_yz·cosθ
Abaqus的S场数据顺序是(S11, S22, S33, S12, S13, S23)，对应σ_x, σ_y, σ_z, τ_xy, τ_xz, τ_yz。直接从val.data索引得到六个分量，套用公式即可。
3.2 代码落地：一个可以直接复制改写的函数
import math def s_to_cylindrical(data, theta): """ 将全局直角坐标应力分量转换到柱坐标系。 data顺序: S11, S22, S33, S12, S13, S23 theta: 该点相对旋转中心的方位角，单位弧度 """ s11, s22, s33 = data[0], data[1], data[2] s12, s13, s23 = data[3], data[4], data[5] c, s = math.cos(theta), math.sin(theta) c2, s2 = c * c, s * s cs = c * s sr = s11 * c2 + s22 * s2 + 2.0 * s12 * cs st = s11 * s2 + s22 * c2 - 2.0 * s12 * cs srt = (s22 - s11) * cs + s12 * (c2 - s2) sz = s33 srz = s13 * c + s23 * s stz = -s13 * s + s23 * c return sr, st, sz, srt, srz, stz调用时注意θ的计算要基于相对旋转中心，而不是盲目用全局坐标：
# 假设旋转中心（圆柱轴心）在 (0.5, 0.0) cx, cy = 0.5, 0.0 theta = math.atan2(y - cy, x - cx) sr = s_to_cylindrical(val.data, theta)[0]这里有个常见误区：如果圆柱的中心不在全局坐标原点，直接用atan2(y, x)，那所有点的方位角都算错了，后面的变换自然全部失真。正确做法是先做坐标平移，再求角度，所有变换公式保持一致，但几何参考点必须显式地传给函数。
3.3 本构方向被设置成柱坐标系时，小心重复旋转
还有个容易忽略的坑：如果你的模型在材料属性里通过*ORIENTATION或截面指派，把单元材料方向定义成了圆柱坐标系，那么应力输出可能已经是在该局部柱坐标系下了，而不是全局直角坐标。这时候fieldValue.localCoordSystem不会为空，Abaqus在向ODB写结果时，会把应力旋转到材料方向坐标系。
脚本里处理这种情形时，要加一道判断：
val = stress_field.values[0] if val.localCoordSystem is not None: print('应力已经输出在局部坐标系中，注意不要重复旋转！')对这类输出，你要么直接采用已有的S11作为径向应力，要么先把它变换回全局坐标系再按自己的需求重算。我的建议是统一强制模型输出全局坐标结果，然后在后处理脚本里统一做变换，避免一会儿局部一会儿全局导致数据不对。想让输出到ODB的应力用全局坐标，可以在Step输出请求里不要勾选"Use Local CSYS"之类的选项，或在*ORIENTATION指派处检查一下，这点在批量处理多个模型时尤其重要。
4. 位移提取没那么简单：节点解、单位换算与径向分量
4.1 位移场是节点量，不是积分点量
跟应力不同，Abaqus里的位移场U存储在节点上，遍历frame.fieldOutputs['U'].values时，每个value带的是nodeLabel而不是integrationPoint。所以提取径向位移跟你提取的"应力点"位置可能并不严格一致——应力在积分点上，位移在节点上，两者天然不重合。
很多实际工程关心的是某条特定径向路径上的应力与位移曲线，而路径常常定义为从内壁到外壁的一段节点序列。我常用的处理方式有两种：
- 如果模型是轴对称的，直接提取目标节点的U1作为径向位移，同时取落在这些节点附近的单元积分点应力（按单元中心或最近原则匹配），然后按半径排序画在一张图里。
- 如果确实需要积分点处的位移，可以用相邻节点位移做线性插值，或者根本不需要这么精确时，直接用节点位移。对绝大多数工程判断，节点位移足够。
位移场的value还有position=NODE的属性标记，写代码时建议确认一下，避免误把某些经过外推的单元变量当位移来读。
4.2 单位换算：SI制与mm制的巨大差异
Abaqus本身没有内置单位制，你建模时用毫米、千克、秒、牛顿的mm制，结果文件里的应力就是MPa，位移就是mm；用SI制的米、千克、秒、牛顿，应力就是Pa，位移就是m。有的脚本提取出来的径向应力数量级是十几，还有人提取出的是上千万，其实都是同一个模型，只是单位不同。
我在实际处理时会在脚本开头定义一个单位转换字典：
unit = 1.0 # 若是SI制，保持Pa和m # unit_mm = 1e3? 如果模型用mm，应力按MPa输出，想转成Pa需乘1e6，位移乘1e-3 stress_scale = 1e6 if abaqus_units == 'mm' else 1.0 disp_scale = 1e-3 if abaqus_units == 'mm' else 1.0如何知道模型单位？看几何坐标数量级最直观。如果模型坐标在几十到几百的范围，多半是mm制；如果在0.01到1之间，多半是m制。提交计算前在CAE左侧模型树里看下part尺寸基本能判断。这个细节不处理，理论解校验、多模型数据对比时会被坑得很惨。
4.3 三维模型位移的径向分量
三维模型的位移默认是全局直角坐标分量(U1, U2, U3)，径向位移同样需要旋转：u_r = U1·cosθ + U2·sinθ。乘以方位角的cos和sin即可，非常简单，但别忘了它针对的是该节点相对圆柱轴心的方位角。
如果模型是轴对称模型，U1就是径向位移，直接取用。只有三维或平面模型的U才需要做这个变换。
5. 把提取逻辑做成流水线：完整脚本与CSV落盘
5.1 一条可复现的提取主线
上面的逻辑拆开讲了不少，实际落地时我更倾向把它们收拢进一个脚本，输入ODB路径和旋转中心，输出一张CSV表格。表格第一列是半径，第二列是径向应力，第三列是径向位移，方便导入Excel或直接用matplotlib、Origin画图。
以下是我整理后的核心脚本框架，把轴对称模型和三维模型统一处理：
# -*- coding: utf-8 -*- import math import csv from odbAccess import openOdb def extract_r_stress_disp(odb_path, step_name, center=(0.0,0.0), out='extract.csv'): odb = openOdb(odb_path, readOnly=True) frame = odb.steps[step_name].frames[-1] S = frame.fieldOutputs['S'] U = frame.fieldOutputs['U'] rows = [] cx, cy = center for inst in odb.rootAssembly.instances.values(): types = {e.type for e in inst.elements} is_axisym = all(t.startswith('CAX') for t in types) label_to_node = {} if is_axisym: # 轴对称：S11直接是径向应力，U1直接是径向位移 node_ur = {} for uv in U.values: if uv.instance.name == inst.name: node_ur[uv.nodeLabel] = uv.data[0] for sv in S.values: if sv.instance.name != inst.name: continue elements = inst.elements nlabels = inst.getElementFromLabel(sv.elementLabel).connectivity coords = [inst.getNodeFromLabel(n).coordinates for n in nlabels] r = sum(c[0] for c in coords) / len(coords) sr = sv.data[0] ur = sum(node_ur.get(n, 0.0) for n in nlabels) / len(nlabels) rows.append([r, sr, ur]) else: # 三维：需要按柱坐标旋转 for sv in S.values: if sv.instance.name != inst.name: continue nlabels = inst.getElementFromLabel(sv.elementLabel).connectivity coords = [inst.getNodeFromLabel(n).coordinates for n in nlabels] x = sum(c[0] for c in coords) / len(coords) y = sum(c[1] for c in coords) / len(coords) theta = math.atan2(y - cy, x - cx) r = math.hypot(x - cx, y - cy) sr = s_to_cylindrical(sv.data, theta)[0] # 节点位移径向分量的均值近似 ur_list = [] for n in nlabels: uv = U.values # 实际应建立nodeLabel到U值的映射 # 简写：实际代码可先构建字典 rows.append([r, sr, ur]) rows.sort(key=lambda row: row[0]) with open(out, 'w', newline='') as f: writer = csv.writer(f) writer.writerow(['radius', 'radial_stress', 'radial_disp']) writer.writerows(rows) print('已输出', out)实际使用三维模型时，位移部分建议先构建一个节点位移字典再按节点查询，避免在循环里频繁遍历U.values。
5.2 多帧提取：从收敛结果到加载历史的完整曲线
有时候我们需要的不只是最后一帧，而是整个加载过程中某个积分点或节点位置上的径向应力、径向位移演化曲线。比如风偏、土压力等工况下，载荷是逐步施加的，最后一步最大值固然重要，但中间历史也常常需要回看。
多帧提取的做法很简单，在外面套一层帧循环：
for i, frame in enumerate(odb.steps[step_name].frames): if frame.fieldOutputs.get('S') is None: continue # 单帧内的提取逻辑同上，额外在行数据里加上时间 frame.frameValueframe.frameValue返回该帧对应的分析步时间（增量步结束时的时间）或频率值。将所有帧的结果都放进一个列表，写CSV时多出一列time，就能画成"径向应力随加载时间变化"的曲线。
5.3 封装成可复用函数：换模型也能直接用
强烈建议把上述代码封装成三个函数：get_integration_point_r_stress负责应力部分，get_node_r_disp负责位移部分，write_csv_rows负责落盘。这样你换一个新的ODB算例时，只需改朋友路径、分析步名称和旋转中心，其他逻辑完全复用。
我自己的模板里还会额外输出一个summary.txt，记录模型类型、单元数量、输出帧数、所用坐标系、中心点坐标，方便后来复盘。调试阶段全靠这个文件确认本次提取的条件，避免数据出处含糊。
6. 实战中容易翻车的六个地方：从场输出陷阱到坐标原点
6.1 场输出列表里根本没有S：变量被压缩或未请求
新手最容易踩的坑：ODB打开后直接fieldOutputs['S']，结果报KeyError。原因是计算前在Step模块的场输出请求里没有勾选"Stress"；或者为了减小ODB体积，只输出了某些里程对应的结果。更隐蔽的是开启*OUTPUT, FREQUENCY后，某些增量步之间没有输出数据，会导致fieldOutputs['S'].values为空列表，但又不报错。
解决办法是脚本里先做个健壮性判断：
if 'S' not in frame.fieldOutputs: print(f'帧 {frame.frameId} 没有应力输出，跳过')如果整个分析都找不到S，回到Step模块检查场输出请求，把S、U、COORD全部勾上重新算一遍。靠后处理脚本是没法凭空造出没写入ODB的数据的。
6.2 输出频率设置导致结果稀疏：增量步多但帧很少
有人算了一个很长的动力分析，结果ODB里只有几帧数据，画出来的应力曲线稀稀拉拉。这往往是场输出频率设置成了N=100或更稀疏。对提取径向应力的需求来说，如果只关心准静态峰值，问题不大；但如果要分析加载历史，就得在Step模块把场输出频率调高，或者在inp里用*OUTPUT, FIELD, FREQUENCY=1之类的设置。不然脚本写得再好也是巧妇难为无米之炊。
6.3 变换公式里的θ方向反了：正负号影响应力数量级
张量旋转公式里的sin/cos平方项，决定了旋转方向。如果你用的θ是-atan2(y,x)，会导致环向应力和剪应力的正负号翻转，径向正应力本身因为cos²/sin²的关系影响不到，但一旦模型对称性不是很好，或者你后续要做强度校核，这个符号问题就是致命的。
统一建议：以全局X轴为0度方向，以X-Y平面内逆时针为正，这样θ =atan2(y - cy, x - cx)。脚本里不要多个地方各自算θ，而是统一写一个get_theta(x, y, cx, cy)函数，全工程调用同一套约定。
6.4 单元类型混杂导致积分点数量不确定
有些模型包含六面体单元、楔形单元、甚至少量四面体单元。不同单元类型的积分点数量不一样，积分点的局部位置也不一样。如果你用单元平均节点坐标代表积分点位置，对楔形和四面体单元的误差会明显偏大。
处理这类混杂网格时，我一般先把单元按类型分组，分别统计。对六面体、四面体分别采用不同的插值方案，或者干脆只提取主要单元类型（比如六面体网格）的目标路径，避免把不同单元类型的应力混杂进同一条曲线里。脚本里可以按val.baseElementType判断：
if sv.baseElementType not in ('C3D8R', 'C3D8'): continue # 只取六面体单元6.5 旋转中心选错：全模型应力分布看似正常但路径曲线全乱
这是最隐蔽也最要命的一个坑。三维厚壁筒如果建模时圆心在(0,0)，但你随口把center设成了(0.5,0)，那么所有点的θ角都错了，径向应力曲线自然面目全非。更糟糕的是，由于模型本身有对称性，如果你提取的是环向一圈的数据，平均值也许还能看，沿径向逐点取的曲线就完全不能用了。
我在每个脚本里都强制打印中心点坐标，并且用一段校验逻辑：找到模型最内层的节点，计算r，然后打印该节点的径向应力，和理论解对比一下。对比通过才继续大批量提取，这个习惯救过我很多次。
6.6 不要忘记用理论解或ABAQUS后处理交叉验证
脚本跑完输出一大堆数据，如果直接拿去写报告，一旦变换公式有偏差，很难发现。我的建议是每次设计一个新算例的提取脚本时，先找一个能解析求解的简单模型（厚壁筒、带孔平板等）做交叉验证。
比如厚壁筒内压问题，理论解可以直接算出内壁处的径向应力等于内压的负值（-p），环向应力等于p*(ri²+ro²)/(ro²-ri²)之类的表达式，把脚本提取结果往上一比对，数量级、正负号、随半径的变化趋势都能一目了然。哪怕只是验证一个点，也能说明整条提取链路是通的。之后再拿这个脚本处理复杂模型，心里就有底了。
实际做过几次之后，我的体会是：提取径向应力和位移这件事，难点不在Abaqus的API调用本身，而在于你对自己模型的理解——到底是轴对称还是三维，应力存在积分点还是节点上，坐标系旋转中心在哪里，输出单位是什么。这四个问题想清楚，脚本半小时就能写完；想不清楚，哪怕抄一段现成代码，出来也是错误结果。最后再分享一个小技巧：给新算例写提取脚本时，先打印几个关键点的原始分量值和方位角，肉眼确认过坐标转换逻辑无误，再运行全量提取，会省下很多排查时间。
