ansys土木工程应用实例_BKM ? 35期 — Python解析ANSYS文件

2020年10月3日作者注释：

官方API已更新，本文部分函数已无法调用，请参考官方API，本文仅提供路径参考，感谢。

——————————————————————

ansys是融结构、流体、电场、磁场、声场分析于一体的大型通用有限元分析软件。在核工业、铁道、石油化工、航空航天、机械制造、能源、汽车交通、国防军工、电子、土木工程、造船、生物医学、轻工、地矿、水利、日用家电等领域有着广泛的应用。

最近需要解析一下ansys生成的文件，包括cdb和rst文件，但是没有学习过ansys方面的开发，意外中找到了Python关于解析ansys文件的包pyansys。国内目前没有关于此包的使用方法，借助此机会介绍Python通过pyansys解析ansys文件。

这是官方链接：https://akaszynski.github.io/pyansys/ansys_control.html

一、准备材料：

1. 安装pyansys

pip install pyansys 

2. cdb文件，rst文件，full文件（有需要文末获取）

二、ANSYS APDL交互式控件

1. 非交互式绘图，需要安装ansys软件（目前来看需要15及其以上版本）。

在生成几何体和模数时绘制几何体和等点通常很有用，并且为了调试（或脚本）目的，在 pyansys 中绘制也很有用。要启用交互式绘图，在启动 ANSYS 时设置interactive_plotting=True。绘制命令（如 APLOT、EPLOT 和 KPLOT）将打开一个 matploblib。

import 

运行完成后

2. 互动断点（需安装ansys软件）

在大多数情况下，尤其是在生成几何图形时，不打开APDL GUI就是不可能的。使用内联绘图很难识别几何图形项，因此pyansys有一种open_gui方法可以无缝打开GUI，而又不会丢失工作或重新启动会话，可以使用以下代码。

import pyansys

# run ansys with interactive plotting enabled
ansys = pyansys.Mapdl()

# create a square area using keypoints
ansys.prep7()
ansys.k(1, 0, 0, 0)
ansys.k(2, 1, 0, 0)
ansys.k(3, 1, 1, 0)
ansys.k(4, 0, 1, 0)
ansys.l(1, 2)
ansys.l(2, 3)
ansys.l(3, 4)
ansys.l(4, 1)
ansys.al(1, 2, 3, 4)

# open up the gui
ansys.open_gui()

# it resumes where you left off...
ansys.et(1, 'MESH200', 6)
ansys.amesh('all')
ansys.eplot()

这种方法避免了必须在交互式会话和脚本会话之间来回切换的麻烦。相反，可以进行一个脚本会话，并从脚本会话中打开GUI，而不会丢失工作或进度。此外，在GUI中所做的任何更改都不会影响脚本。可以在GUI中进行实验，python脚本不会受到影响。

3. 运行批处理（需要安装ansys软件）

可以定义运行ansys的函数，而不是通过使用输入文件调用ANSYS来运行ANSYS批处理。本示例基于具有扭转载荷的圆柱体的最大应力进行网格收敛研究。

import numpy as np
import pyansys

def cylinder_batch(elemsize, plot=False):
    """ Report the maximum von Mises stress of a Cantilever supported cylinder"""

    # clear
    ansys.finish()
    ansys.clear()

    # cylinder parameters
    radius = 2
    h_tip = 2
    height = 20
    force = 100/radius
    pressure = force/(h_tip*2*np.pi*radius)

    ansys.prep7()
    ansys.et(1, 186)
    ansys.et(2, 154)
    ansys.r(1)
    ansys.r(2)

    # Aluminum properties (or something)
    ansys.mp('ex', 1, 10e6)
    ansys.mp('nuxy', 1, 0.3)
    ansys.mp('dens', 1, 0.1/386.1)
    ansys.mp('dens', 2, 0)

    # Simple cylinder
    for i in range(4):
        ansys.cylind(radius, '', '', height, 90*(i-1), 90*i)

    ansys.nummrg('kp')

    # mesh cylinder
    ansys.lsel('s', 'loc', 'x', 0)
    ansys.lsel('r', 'loc', 'y', 0)
    ansys.lsel('r', 'loc', 'z', 0, height - h_tip)
    # ansys.lesize('all', elemsize*2)
    ansys.mshape(0)
    ansys.mshkey(1)
    ansys.esize(elemsize)
    ansys.allsel('all')
    ansys.vsweep('ALL')
    ansys.csys(1)
    ansys.asel('s', 'loc', 'z', '', height - h_tip + 0.0001)
    ansys.asel('r', 'loc', 'x', radius)
    ansys.local(11, 1)
    ansys.csys(0)
    ansys.aatt(2, 2, 2, 11)
    ansys.amesh('all')
    ansys.finish()

    if plot:
        ansys.view(1, 1, 1, 1)
        ansys.eplot()

    # new solution
    ansys.slashsolu()
    ansys.antype('static', 'new')
    ansys.eqslv('pcg', 1e-8)

    # Apply tangential pressure
    ansys.esel('s', 'type', '', 2)
    ansys.sfe('all', 2, 'pres', '', pressure)

    # Constrain bottom of cylinder/rod
    ansys.asel('s', 'loc', 'z', 0)
    ansys.nsla('s', 1)

    ansys.d('all', 'all')
    ansys.allsel()
    ansys.psf('pres', '', 2)
    ansys.pbc('u', 1)
    ansys.solve()
    ansys.finish()

    # access results using ANSYS object
    result = ansys.result

    # to access the results you could have run:
    # resultfile = os.path.join(ansys.path, '%s.rst' % ansys.jobname)
    # result = pyansys.read_binary(result file)

    # Get maximum von Mises stress at result 1
    nodenum, stress = result.principal_nodal_stress(0)  # 0 as it's zero based indexing

    # von Mises stress is the last column
    # must be nanmax as the shell element stress is not recorded
    maxstress = np.nanmax(stress[:, -1])

    # return number of nodes and max stress
    return nodenum.size, maxstress


# initialize ANSYS
ansys = pyansys.Mapdl(override=True, loglevel='error')

result_summ = []
for elemsize in np.linspace(0.6, 0.15, 15):
    # run the batch and report the results
    nnode, maxstress = cylinder_batch(elemsize, plot=False)
    result_summ.append([nnode, maxstress])
    print('Element size %f: %6d nodes and maximum vom Mises stress %f'
          % (elemsize, nnode, maxstress))

# Exit ANSYS
ansys.exit()

运行结果

Element size 0.600000:   9657 nodes and maximum vom Mises stress 142.623505
Element size 0.567857:  10213 nodes and maximum vom Mises stress 142.697800
Element size 0.535714:  10769 nodes and maximum vom Mises stress 142.766510
Element size 0.503571:  14177 nodes and maximum vom Mises stress 142.585388
Element size 0.471429:  18371 nodes and maximum vom Mises stress 142.825684
Element size 0.439286:  19724 nodes and maximum vom Mises stress 142.841202
Element size 0.407143:  21412 nodes and maximum vom Mises stress 142.945984
Element size 0.375000:  33502 nodes and maximum vom Mises stress 142.913437
Element size 0.342857:  37877 nodes and maximum vom Mises stress 143.033401
Element size 0.310714:  59432 nodes and maximum vom Mises stress 143.328842
Element size 0.278571:  69106 nodes and maximum vom Mises stress 143.176086
Element size 0.246429: 110547 nodes and maximum vom Mises stress 143.499329
Element size 0.214286: 142496 nodes and maximum vom Mises stress 143.559128
Element size 0.182143: 211966 nodes and maximum vom Mises stress 143.953430
Element size 0.150000: 412324 nodes and maximum vom Mises stress 144.275406

4. ANSYS对象方法

classpyansys.Mapdl(exec_file=None, run_location=None, jobname='file', nproc=2, override=False, loglevel='INFO', additional_switches='', start_timeout=120, interactive_plotting=False, log_broadcast=False, check_version=True, prefer_pexpect=True, log_apdl='w')

该类在后台打开ANSYS，并允许将命令传递给持久会话。

exec_file (str, optional)– ANSYS可执行文件的位置。保留为默认值None时将使用缓存的位置。
run_location (str, optional) – ANSYS工作目录。默认为临时工作目录。
jobname (str, optional) – ANSYS作业名。默认为'file'。
nproc（int ，optional）–处理器数量。默认为2。
override (bool, optional)–尝试删除run_location处的锁定文件。在先前的ANSYS会话过早退出并且未删除锁定文件时很有用。
wait（bool ，optional）–为True时，等待初始化ANSYS之后再初始化python ansys对象。将此设置为False进行调试。
loglevel（str ，optional）–设置将哪些消息打印到控制台。默认的“ INFO”打印所有ANSYS消息，“ WARNING”仅打印包含ANSYS警告的消息，“ ERROR”仅打印错误消息。
Additional_switches（str ，optional）–用于ANSYS的其他开关，例如aa_r和学术研究许可证，将添加：
Additional_switches =”-aa_r”避免添加-i -o或-b之类的开关，因为它们已经包含在启动ANSYS MAPDL服务器中。
start_timeout（float ，optional）–引发ANSYS无法启动的错误之前要等待的时间。
interactive_plotting (bool, optional)–使用启用交互式绘图matplotlib。默认为False。
log_broadcast（bool ，optional）– ansys解决方案进度的其他日志记录。默认值True，在日志级别“ INFO”下可见。
check_version（bool ，optional）–检查二进制文件的版本，并在无效时引发异常。
preferred_pexpect（bool ，optional）–启用后，将避免在CORBA服务器模式下使用ansys APDL，并将生成一个进程并使用pexpect控制它。默认为False。
log_apdl（str ，optional）–在当前ANSYS工作目录中打开一个APDL日志文件。默认为“ w”。设置为“ a”以追加到现有日志。

三、读写ANSYS档案（不需要安装ansys软件）

1. 阅读ANSYS档案

可以使用Archive加载包含实体元素（传统元素和现代元素）的ANSYS存档文件，然后将其转换为vtk对象：

import pyansys
from pyansys import examples

# Sample *.cdb
filename = 'examples/sector.cdb'

# Read ansys archive file
archive = pyansys.Archive(filename)

# Print raw data from cdb
for key in archive.raw:
   print("%s : %s" % (key, archive.raw[key]))

# Create a vtk unstructured grid from the raw data and plot it
grid = archive.parse_vtk(force_linear=True)
grid.plot(color='w', show_edges=True)

输出运行结果

此外还可以通过启用read_parameters参数来有选择地读入存档文件中的任何已存储参数。

import pyansys archive = pyansys.Archive('examples/HexBeam.cdb') # parameters are stored as a dictionary print(archive.raw.keys())

import pyansys
archive = pyansys.Archive('examples/HexBeam.cdb')

# parameters are stored as a dictionary
print(archive.raw.keys())

所有的参数以字典形式存在，上述代码中输出了所有的key值，对应可以获取value。

2. 编写ANSYS档案

使用VTK生成的非结构化网格可以转换为ANSYS APDL存档文件，并使用加载到任何版本的ANSYS中pyansys.save_as_archive。以下使用内置存档文件的示例演示了此功能。

import pyvista as pv
from pyvista import examples
import pyansys

grid = pv.UnstructuredGrid(examples.hexbeamfile)
script_filename = '/tmp/grid.cdb'
pyansys.save_as_archive(script_filename, grid)

# read in archive in ANSYS and generate cell shape quality report
ansys = pyansys.ANSYS()
ansys.cdread('db', script_filename)
ansys.prep7()
ansys.shpp('SUMM')

3. 使用ANSYS结果文件（rst）

ANSYS结果文件是FORTRAN格式的二进制文件，其中包含从ANSYS分析写入的结果。结果至少包含所分析模型的几何形状以及节点和单元结果。根据分析，这些结果可能是从模态位移到节点温度的任何结果。目前，仅支持以下结果：

• 节点自由度来自静态分析或模态分析。
• 节点自由度来自循环静态或模态分析。
• 节点平均分量应力（即x，y，z，xy，xz，yz）
• 节点主应力（即S1，S2，S3，SEQV，SINT）

# Load the reader from pyansys
import pyansys
from pyansys import examples

# Sample result file
rstfile = 'examples/file.rst'

# Create result object by loading the result file
result = pyansys.read_binary(rstfile)

# Beam natural frequencies
print(result.resultheader)
result.plot_nodal_solution(0, 'x', label='Displacement')

其中result可以获取以下属性。

result.resultheader #结果标头 result.time_values #分析的时间和频率值 result.nnum #排序节点 result.enum #元素编号 result.geometry #几何内容

其中result.geometry包含以下键：

'nnum' （已排序的节点编号）

• 'nodes' （节点位置）
• 'etype' （元素类型）
• 'enum' （与elem数组关联的未排序元素编号）
• 'elem' （numpy数组显示与每个元素关联的节点，-1表示未使用的条目）
• 'ekey' （2xn元素类型引用数组）
• 'coord systems' （坐标系）

本期内容和BIM没有太大关系，但是也需要思考怎样使用Python将解析出来的结果附着到BIM模型中去，比如IFC或者常用的Revit软件。