技术求助：用Python绘制Abaqus INP文件的10节点四面体网格并映射应力应变

针对你的大规模Abaqus模型（近百万C3D10单元），以下是用Python实现网格绘制及应力/应变映射的可行方案，适配多组仿真结果对比需求：

核心依赖库 #

优先选择PyVista（基于VTK，适合大规模3D网格可视化），搭配numpy处理数据：

pip install pyvista numpy

步骤1：解析INP文件关键数据 #

INP文件主要需提取三类信息：节点坐标、单元节点索引、单元应力/应变数据。以下是简化的解析函数（需根据你的INP格式调整，比如场数据的关键字段）：

import numpy as np
import pyvista as pv

def parse_inp(file_path):
    nodes = {}
    elements = []
    el_stress = {}
    current_section = None

    with open(file_path, 'r') as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith('**'):
                continue
            # 识别INP段
            if line.startswith('*NODE'):
                current_section = 'NODE'
                continue
            elif line.startswith('*ELEMENT, TYPE=C3D10'):
                current_section = 'ELEMENT'
                continue
            elif line.startswith('*EL FILE,'):
                current_section = 'EL_STRESS'
                # 跳过表头行（比如S11,S22等）
                next(f)
                continue

            # 读取节点数据
            if current_section == 'NODE':
                parts = line.split(',')
                node_id = int(parts[0].strip())
                x, y, z = map(float, [p.strip() for p in parts[1:4]])
                nodes[node_id] = (x, y, z)
            # 读取单元数据（C3D10是10个节点）
            elif current_section == 'ELEMENT':
                parts = line.split(',')
                el_id = int(parts[0].strip())
                node_ids = [int(p.strip()) for p in parts[1:11]]
                elements.append((el_id, node_ids))
            # 读取单元应力数据（示例为S11，按需修改）
            elif current_section == 'EL_STRESS':
                parts = line.split(',')
                el_id = int(parts[0].strip())
                s11 = float(parts[1].strip())
                el_stress[el_id] = s11

    # 转换为numpy数组适配PyVista
    node_coords = np.array([nodes[pid] for pid in sorted(nodes.keys())])
    # PyVista的单元格式：先指定单元类型（VTK_QUADRATIC_TETRA=24），再按节点索引（注意要减1，因为VTK是0基）
    cells = []
    for el_id, node_ids in elements:
        cells.append(10)  # 单元节点数
        cells.extend([pid-1 for pid in node_ids])  # 转换为0基索引
    cells = np.array(cells)

    # 创建PyVista网格对象
    grid = pv.UnstructuredGrid(cells, [pv.CellType.QUADRATIC_TETRA]*len(elements), node_coords)
    # 添加应力数据到网格
    stress_array = np.array([el_stress[el_id] for el_id, _ in elements])
    grid['S11'] = stress_array

    return grid

步骤2：单组结果可视化 #

解析完成后，可直接绘制带应力映射的3D网格：

# 读取单个INP文件
grid = parse_inp('your_model.inp')

# 创建绘图窗口
plotter = pv.Plotter()
# 添加网格，按S11着色
plotter.add_mesh(grid, scalars='S11', cmap='viridis', show_edges=False)
# 添加颜色条
plotter.add_scalar_bar(title='Stress S11')
# 显示
plotter.show()

步骤3：多组结果对比（迭代仿真） #

针对多组载荷下的结果，可创建子图对比：

# 读取多组仿真结果
grid_list = [parse_inp(f'model_load_{i}.inp') for i in range(3)]
titles = ['Load Case 1', 'Load Case 2', 'Load Case 3']

# 创建多子图窗口
plotter = pv.Plotter(shape=(1,3))
for i in range(3):
    plotter.subplot(0, i)
    plotter.add_mesh(grid_list[i], scalars='S11', cmap='viridis', clim=[grid['S11'].min(), grid['S11'].max()])
    plotter.add_scalar_bar(title='S11')
    plotter.add_title(titles[i])
plotter.link_views()  # 同步视角
plotter.show()

性能优化建议 #

• 对于百万级单元，建议关闭边显示（show_edges=False），减少渲染压力；
• 若内存不足，可使用PyVista的grid.extract_subset()提取部分单元预览，或使用pyvista.global_theme.load_theme('document')降低渲染精度；
• 解析INP时可采用分块读取，避免一次性加载全部数据到内存。

内容的提问来源于stack exchange，提问作者siddarth swaminathan