当前位置：服务支持 > 软件文章 > 用Abaqus算颈椎节段活动度，基于节点向量变化的方法挺好用

用Abaqus算颈椎节段活动度，基于节点向量变化的方法挺好用

阅读数 8

引言

在脊髓损伤的生物力学仿真中，颈椎模型在屈伸、侧弯或旋转运动中的节段活动度（Range of Motion, ROM）是一个至关重要的参数。它不仅是验证模型生物逼真度的关键指标，也是不同工况下评估脊髓应力应变状态的输入条件。

在我学生时代，曾使用Abaqus进行颈椎-脊髓系统的有限元计算仿真。其中一个基础但棘手的问题就是：如何从仿真结果中准确、可靠地提取出各椎体节段的角度随时间变化曲线？ 本文旨在回顾并分享当时所采用的一种方法，其核心是通过计算刚体椎体上特征向量夹角的变化来反推节段角度，该方法有效避免了常用简化方法带来的误差，并与实验测量原理保持一致。

一、问题背景与技术路线的选择

在Abaqus中获取各节段相对角度，课题组内主要有两种技术路线：

连接器单元法：在椎体参考点之间建立连接器（Connector），并定义铰链关节等行为(具体关节名称我已经遗忘)，直接输出连接器上的相对旋转角度。这种方法在单节段分析中非常直观有效。但其局限性在于：连接器依赖一个固定的参考坐标系，而在多节段、大范围的复杂运动中，坐标系不随运动移动会导致角度计算累积误差，使得角度计算复杂化。
节点坐标计算法：直接读取Abaqus输出数据库（.odb文件）中特定节点的坐标随时间变化的历史，通过数学计算得到角度。这正是我所采用的方案。它的优势在于不依赖任何可能产生歧义的中间单元或参考系，直接处理最原始的节点位移数据，原理清晰，结果稳健。

二、为什么不是“点点连线”？

很多人的第一直觉是：分别在上下两个椎体上选取一个点，连接这两点形成一个矢量，然后计算相邻节段矢量之间的夹角。如下图所示：

图.1：A.直觉做法：计算相邻矢量V1 (A->B) 和 V2 (B->C) 的夹角 θ'；B.特征向量法；C.传统研究测角方法图示。

这种方法将复杂的椎间盘简化为一个理想的铰链（hinge），它隐含了一个假设：椎体只发生旋转，没有平移。然而，真实的颈椎屈伸运动是耦合的，即同时包含旋转和平移。这种简化会忽略椎体的平移分量，导致计算出的角度θ'大于真实的节段活动角度。

值得注意的是，这种直接连线计算的逻辑其实与采用颈椎活动度计或体表标志物进行测量的实验研究方法一致。然而，此类方法通常只能提供整个颈椎的总活动角度，而无法准确反映各节段的具体活动情况，因此往往无法在验证研究中使用。

三、我们的方法：基于椎体自身向量的变化

我们的模型并不发生椎体破坏，因此椎体被简化为了刚体以节约了计算资源，也符合我们研究聚焦于脊髓生物力学的需求。刚体的一个重要特性是：其内部任意两点间的向量方向与整个刚体的旋转直接相关。

因此，我们方法的思路在于：

为每一个椎体定义一个特征向量：在每个椎体上选取两个不重合的节点，由它们构成一个向量 V。这个向量可以代表该椎体的方向。
追踪向量方向的变化：在仿真过程中，通过读取这两个节点的位移场，计算出每一时刻该向量 V(t) 的新方向。
计算相对角度：对于相邻的两个椎体（例如C4和C5），我们计算C4椎体的向量方向和C5椎体的向量方向之间的夹角。这个夹角的变化量就是C4-C5节段的相对活动度。

这种方法的核心在于：每个椎体的运动（旋转+平移）信息都被其自身向量的方向变化所捕获。当我们计算两个椎体向量之间的夹角时，自然就剔除掉了它们共同的平移运动分量，得到的正是我们想要的纯旋转角度。

这与实验生物力学中测量节段活动度的经典方法——叠加X光片法——在原理上高度一致(图1C)。实验人员通过叠印对齐同一椎体来消除平移，再测量角度的变化。我们的计算方法可视为这一实验原理在数字仿真中的完美实现。

四、代码实现与解析（Abaqus Script）

以下是我们为实现上述计算所编写的Abaqus Script脚本。它通过odbAccess模块读取计算结果，并执行上述向量运算。

Abaqus Script 本质上是基于 Python 的派生脚本语言，随 Abaqus 一同安装部署。在调用脚本处理 ODB 文件时，同样遵循“高版本兼容低版本”的规律。因此，建议生成 ODB 的 Abaqus 版本与分析所使用的版本保持一致，以避免潜在的兼容性问题。需要注意的是，我们当时选用的Abaqus(6.14-5)求解器所附带的Abaqus Script其Python版本较低，不支持中文注释。示例代码中的中文注释仅为便于理解，在实际应用时需要删除。

代码块

Python

自动换行

复制代码

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224

# -*- coding: utf-8 -*-

"""
ABAQUS ODB 后处理脚本
=================================

目的:
本脚本用于自动计算 Abaqus 输出数据库 (ODB) 文件中颈椎（C1-C7）相邻节段间的
相对转动角度。它会处理指定分析步中的每一帧，计算出预先定义的椎体节段的角度变化。

如何运行:
在 Abaquse CAE 的命令行窗口或系统命令行中，使用 Abaqus Python 解释器执行此脚本。
abaqus python <脚本文件名>.py <odb文件名>.odb

示例:
abaqus python calculate_angles.py spine_simulation.odb

日期: 2022-09-12
"""

# =============================================================================
# 1. 导入库与初始设置
# =============================================================================
# 将所有需要用到的库统一在脚本开头导入，便于管理。
import sys
from odbAccess import openOdb
import numpy as np
from math import acos, degrees

# =============================================================================
# 2. 用户自定义参数 (配置区域)
# =============================================================================
# 将所有可能需要修改的变量集中放在这里，方便用户使用，无需深入代码内部。

# 指定需要处理的分析步名称。
# 这个名称必须与您在 Abaqus 模型中定义的分析步名称完全一致。
STEP_NAME = 'Step-1'

# --- 选择计算角度的解剖平面 ---
# 使用描述性的字符串来定义平面，比使用数字 (0, 1, 2) 更直观，不易出错。
# 可选值: 'SAGITTAL' (矢状面), 'CORONAL' (冠状面), 'HORIZONTAL' (水平面)
PLANE_OF_INTEREST = 'SAGITTAL'

# --- 定义节段节点对 ---
# 定义代表每个椎体节段的节点对。脚本将基于这两个节点构建一个方向向量。
# 格式: ['实例名', 节点1标签, 节点2标签]
SEGMENT_NODE_PAIRS = [
    ['C1-1', 28436,  28462],   # C1 节段
    ['C2-1', 415019, 415069],  # C2 节段
    ['C3-1', 454031, 453996],  # C3 节段
    ['C4-1', 159373, 159342],  # C4 节段
    ['C5-1', 303573, 303550],  # C5 节段
    ['C6-1', 218572, 252133],  # C6 节段
    ['C7-1', 263059, 262972]   # C7 节段 (或作为参考的胸椎)
]

# =============================================================================
# 3. 脚本初始化与验证
# =============================================================================

# --- 检查命令行参数 ---
# 确保用户在运行脚本时，提供了 ODB 文件的路径作为参数。
if len(sys.argv) != 2:
    print("错误：您必须在命令行中提供 ODB 文件的名称。")
    print("用法: abaqus python <脚本文件名>.py <odb文件名>.odb")
    sys.exit(1) # 使用非零代码退出，表示程序因错误而终止。

odb_path = sys.argv[1]

# --- 平面选择逻辑 ---
# 使用字典将人类可读的平面名称映射到计算所需的索引和描述性名称。
# 这种方法结构清晰，易于扩展和维护。
plane_map = {
    # 水平面内的转动，是 Z 轴的旋转，因此关注 X 和 Y 坐标分量。
    'HORIZONTAL': {'index': 0, 'name': 'Z 平面 (水平面)'},
    # 矢状面内的转动（屈伸），是 X 轴的旋转，因此关注 Y 和 Z 坐标分量。
    'SAGITTAL':   {'index': 1, 'name': 'X 平面 (矢状面)'},
    # 冠状面内的转动（侧屈），是 Y 轴的旋转，因此关注 Z 和 X 坐标分量。
    'CORONAL':    {'index': 2, 'name': 'Y 平面 (冠状面)'}
}

if PLANE_OF_INTEREST not in plane_map:
    print(f"错误: 无效的平面名称 '{PLANE_OF_INTEREST}'")
    print("请从 'SAGITTAL', 'CORONAL', 'HORIZONTAL' 中选择。")
    sys.exit(1)

# 根据用户选择，获取对应的索引和描述性名称。
plane_index = plane_map[PLANE_OF_INTEREST]['index']
plane_name = plane_map[PLANE_OF_INTEREST]['name']

# --- 打印脚本运行信息 ---
print(f"正在处理 ODB 文件: {odb_path}")
print(f"计算分析步: [{STEP_NAME}] | 计算平面: [{plane_name}]")

# =============================================================================
# 4. 数据准备与 ODB 访问
# =============================================================================

# --- 打开 ODB 文件 ---
# 使用 try-except 结构来捕获可能发生的错误，例如文件不存在或已损坏。
try:
    # readOnly=True 表示只读模式，速度更快且更安全。
    odb = openOdb(path=odb_path, readOnly=True)
except Exception as e:
    print(f"打开 ODB 文件时出错: {e}")
    sys.exit(1)

# 从 ODB 获取根装配体 (rootAssembly) 对象。
assembly = odb.rootAssembly

# --- 初始化数据结构 ---
# 使用更具描述性的变量名，并利用 NumPy 创建数组以提高计算效率。
num_segments = len(SEGMENT_NODE_PAIRS)
# 该数组用于存储每个节段在模型初始状态（未变形）时的方向向量。
initial_vectors = np.zeros((num_segments, 3))
# 该数组用于存储在每一帧中，每个节段的三个平面角度（弧度制）。
segment_angles_rad = np.zeros((num_segments, 3))

# --- 计算初始方向向量 (第 0 帧) ---
print("\n正在计算初始方向向量...")
for i, (instance_name, node1_label, node2_label) in enumerate(SEGMENT_NODE_PAIRS):
    # 从 Abaqus 装配体结构中获取节点对象。
    node1 = assembly.instances[instance_name].getNodeFromLabel(node1_label)
    node2 = assembly.instances[instance_name].getNodeFromLabel(node2_label)
    
    # 计算初始向量 (从节点2指向节点1)。
    # 这个初始状态是后续所有角度计算的参考基准。
    initial_vector = np.array(node1.coordinates) - np.array(node2.coordinates)
    initial_vectors[i] = initial_vector

# =============================================================================
# 5. 主处理循环：逐帧分析
# =============================================================================

# --- 准备输出内容的表头 ---
print('帧号       C1-C2      C2-C3      C3-C4      C4-C5      C5-C6      C6-C7')
print('-' * 72)

# 访问指定分析步中的所有帧。
try:
    step_frames = odb.steps[STEP_NAME].frames
except KeyError:
    print(f"错误: 在 ODB 文件中未找到名为 '{STEP_NAME}' 的分析步。")
    odb.close()
    sys.exit(1)
    
# 遍历分析步中的每一帧。
# 使用 enumerate 可以同时获得帧的索引 (frame_index) 和帧对象本身 (frame)。
for frame_index, frame in enumerate(step_frames):
    # 获取当前帧的位移场输出 ('U')。
    displacement_field = frame.fieldOutputs['U']
    
    # 创建一个列表，用于存储当前帧计算出的所有相对角度，最后统一打印。
    relative_angles_deg_for_frame = []

    # --- 计算每个节段在当前帧的变形向量和角度 ---
    for i, (instance_name, n1_label, n2_label) in enumerate(SEGMENT_NODE_PAIRS):
        # 再次获取节点对象，以便查询它们的位移。
        node1_instance = assembly.instances[instance_name].getNodeFromLabel(n1_label)
        node2_instance = assembly.instances[instance_name].getNodeFromLabel(n2_label)

        # 从位移场中提取这两个节点的位移数据。
        disp1 = displacement_field.getSubset(region=node1_instance).values[0].data
        disp2 = displacement_field.getSubset(region=node2_instance).values[0].data

        # 计算变形后的新坐标：新坐标 = 初始坐标 + 位移。
        deformed_coord1 = np.array(node1_instance.coordinates) + np.array(disp1)
        deformed_coord2 = np.array(node2_instance.coordinates) + np.array(disp2)

        # 计算变形后的新方向向量。
        deformed_vector = deformed_coord1 - deformed_coord2
        
        # --- 使用向量投影法计算角度 ---
        # 我们通过计算初始向量和变形后向量在各个解剖平面上的投影之间的夹角，来确定转动角度。
        
        # 为了代码清晰，定义向量分量
        v_init = initial_vectors[i]
        v_deformed = deformed_vector
        
        # 投影到 XY 平面 (用于计算水平面/Z轴转角)
        v_init_xy = v_init[:2]      # [vx, vy]
        v_deformed_xy = v_deformed[:2]
        
        # 投影到 YZ 平面 (用于计算矢状面/X轴转角)
        v_init_yz = v_init[1:]      # [vy, vz]
        v_deformed_yz = v_deformed[1:]
        
        # 投影到 ZX 平面 (用于计算冠状面/Y轴转角)
        v_init_zx = np.array([v_init[2], v_init[0]])       # [vz, vx]
        v_deformed_zx = np.array([v_deformed[2], v_deformed[0]])

        # 使用点积公式计算角度的余弦值: cos(theta) = (A . B) / (|A| * |B|)
        # 注意: np.clip() 函数至关重要！由于浮点数计算可能存在微小的精度误差，
        # 导致计算出的余弦值略微超出[-1.0, 1.0]的范围 (例如 1.00000001)。
        # 这会使 acos 函数报错。np.clip 函数将数值强制限制在此范围内，确保代码的鲁棒性。
        cos_angle_xy = np.clip(np.dot(v_init_xy, v_deformed_xy) / (np.linalg.norm(v_init_xy) * np.linalg.norm(v_deformed_xy)), -1.0, 1.0)
        cos_angle_yz = np.clip(np.dot(v_init_yz, v_deformed_yz) / (np.linalg.norm(v_init_yz) * np.linalg.norm(v_deformed_yz)), -1.0, 1.0)
        cos_angle_zx = np.clip(np.dot(v_init_zx, v_deformed_zx) / (np.linalg.norm(v_init_zx) * np.linalg.norm(v_deformed_zx)), -1.0, 1.0)
        
        # 将计算出的角度（弧度制）存储起来。
        segment_angles_rad[i][0] = acos(cos_angle_xy)
        segment_angles_rad[i][1] = acos(cos_angle_yz)
        segment_angles_rad[i][2] = acos(cos_angle_zx)

    # --- 计算相邻节段间的相对角度 ---
    for i in range(num_segments - 1):
        # 相对角度 = 上一个节段的绝对角度 - 下一个节段的绝对角度。
        relative_angle_rad = segment_angles_rad[i][plane_index] - segment_angles_rad[i+1][plane_index]
        relative_angles_deg_for_frame.append(degrees(relative_angle_rad))

    # --- 打印当前帧的计算结果 ---
    # 使用现代化的 f-string 格式化方法，使输出代码更整洁、易读。
    # map 函数可以高效地将格式化操作应用到列表中的每一个角度值。
    formatted_angles = "   ".join(map(lambda angle: f"{angle:9.2f}", relative_angles_deg_for_frame))
    print(f"{frame_index:<5}    {formatted_angles}")

# =============================================================================
# 6. 清理工作
# =============================================================================
print('-' * 72)
print('--------------------------- 计算结束 ---------------------------')
# 关闭 ODB 文件至关重要！这会释放文件锁和内存。
# 如果忘记关闭，可能导致 Abaqus 在后续操作中变慢、卡死或无法访问该文件。
odb.close()

      复制成功

脚本目标

本脚本旨在自动化处理颈椎模型的 Abaqus 仿真结果。其核心功能是量化每一个仿真时间步（帧）下，各椎间节段（如 C1-C2, C2-C3）在特定解剖平面（矢状面、冠状面或水平面）的相对转动角度。

核心知识点

1、Abaqus 脚本接口：Abaqus 提供了一个强大的 Python API (odbAccess)，允许我们通过编程方式访问和操作输出数据库 (.odb) 文件中的海量数据。这极大地提高了后处理效率，避免了通过图形界面手动提取数据的繁琐工作。

2、向量数学：本次分析的核心是向量计算。我们用两个节点定义一个代表椎体节段的方向向量。该向量与其初始位置的夹角，代表了这个节段的绝对转动。而我们更关心的相对转动（即椎间盘的变形），则是相邻两个节段绝对转动角度的差值。

3、解剖平面：为了分析复杂的3D运动，我们通常将其分解到2D的解剖平面上进行研究。

矢状面 (Sagittal Plane)：观察屈伸运动 (flexion/extension)。我们将3D向量投影到 YZ 平面进行分析。
冠状面 (Coronal Plane)：观察侧屈运动 (lateral bending)。我们将3D向量投影到 ZX 平面进行分析。
水平面 (Horizontal Plane)：观察轴向转动 (axial rotation)。我们将3D向量投影到 XY 平面进行分析。

4、数值计算的鲁棒性：计算机使用浮点数进行计算，这可能引入微小的精度误差。例如，在计算一个角度的余弦值时，结果可能是 1.000000001 或 -1.000000001。虽然误差极小，但这会导致 acos() 函数数学计算失败。因此，我们必须用 clip (钳位) 操作将计算结果严格限制在 [-1.0, 1.0] 的有效区间内，这是编写科学计算程序的重要技巧。

脚本各部分详解

第 1 & 2 部分: 配置区域

目的: 将所有用户可能需要修改的参数（如分析步名、节点编号）集中放置在脚本顶部。这是一个优秀的编程习惯，用户无需阅读和修改核心代码，就能轻松使用此脚本。
脚本使用 'SAGITTAL' 这样的描述性字符串来代替 plane = 1 这样的“魔法数字”。这使得代码本身就具有说明性，降低了出错的概率。我们用一个字典 plane_map 来清晰地管理这些字符串与实际计算逻辑的对应关系。

第 3 部分: 初始化与验证

目的: 该部分负责执行“准备工作”。它首先会检查用户是否在命令行正确地输入了 ODB 文件名。一个健壮的程序应当总是先验证输入，并向用户提供清晰的错误提示。
如果用户运行命令的格式不正确，脚本会打印出明确的“用法”说明。同时，它还会验证用户设置的 PLANE_OF_INTEREST 是否为合法值。

第 4 部分: 数据准备

目的: 在进入主循环之前，我们需要一个计算的基准。此部分打开 ODB 文件，并计算出模型在未变形状态（第0帧）下，每个椎体节段的方向向量。这些初始向量是衡量后续所有运动的“尺子”。
我们使用 NumPy 的 np.zeros() 来预先分配数组内存，这比在循环中动态创建 Python 列表效率更高。变量名 initial_vectors (初始向量) 比 coord0 这样的命名要清晰得多。

第 5 部分: 主处理循环

这是脚本的核心，它会遍历指定分析步中的每一帧。

1、获取位移: 对每一帧，首先从场输出中抓取位移数据 'U'。

2、计算变形状态: 对每个节段，通过“初始坐标 + 位移”计算出节点在当前帧的新坐标，并由此得到新的方向向量。

3、计算角度:

将3D的初始向量和新向量分别投影到我们关心的2D平面上。
利用点积公式计算这两个2D投影向量之间的夹角 θ：
我们使用 np.linalg.norm() 计算向量的模 (即)，用 np.dot() 计算点积。
在送入 acos() 函数之前，必须用 np.clip() 来保证数值的有效性。

4、计算相对角度: 临床上更有意义的是椎体间的相对运动。我们通过计算相邻节段绝对角度的差值（例如：C2的角度 - C3的角度）来得到该值。

5、打印输出: 将当前帧的所有计算结果格式化后，打印到控制台。

脚本使用 NumPy 内置的函数 (np.dot, np.linalg.norm) 进行向量运算，比手动计算高效、简洁且不易出错。打印输出时，使用了 f-string 格式化，这是 Python 现代、推荐的字符串处理方式。

第 6 部分: 清理工作

目的: 这仍然是脚本中重要的一步。计算完成后，必须调用 odb.close() 来关闭 ODB 文件。
如果不关闭文件，Abaqus 可能会在后台一直占用该文件，导致内存泄漏。这会拖慢系统，甚至可能使你无法在图形界面中打开该 ODB 文件，或者在下次运行脚本时失败。

这一Abaqus Script脚本是通过“abaqus python 脚本名.py ODB名称.odb”来调用的。相邻椎体在同一平面上的角度之差，即为该节段（如C1-2, C2-3）的活动角度。结果被逐帧打印出来。

五、总结与启示

这种方法虽然需要额外的后处理脚本，但它提供了极高的精度和可靠性。它严格遵循了刚体运动学原理，并复现了实验测量的逻辑，避免了Abaqus内置功能在某些复杂场景下的局限性。

毕业多年后我早已不再从事相关领域，相比直接进入故纸堆，希望这份回顾和代码分享能对从事相关研究的朋友有所帮助。

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删

返回上级列表

联系我们

，获取更多内容

Abaqus后处理自动化脚本，多区域最大Mises应力和主应变一键追踪

矩阵奇异问题有多致命？数值计算里这个坑必须避开

Abaqus仿真中计算颈椎节段活动度：一种基于节点向量变化的方法