基于节点-连接关系图的3D管系有限元快速仿真系统的开发

2022-08-08颜孙挺陈永贵

压力容器 2022年6期

颜孙挺，唐萍，陈永贵

(1.浙江省特种设备科学研究院，杭州 310020;2.浙江省特种设备安全检测技术研究重点实验室，杭州 310020)

0 引言

压力管道是用于输送高压气体或液体的特种设备，其在初始设计、生产制造和安装运行过程中的疏忽均可引发爆炸或泄漏中毒等事故[1-3]。根据统计，在工业管道严重损坏事故中有60%～65%是由生产制造安装过程中的原初缺陷引起的；35%～40%是由使用中的缺陷和力学损伤引起的[4]，因此，对压力管道进行准确可靠的安全评定十分重要。

GB/T 19624—2019《在用含缺陷压力容器安全评定》为压力容器的安全评定提供了可靠的技术手段，但其评价过程较为复杂，须完成应力线性化后，获取薄膜应力、弯曲应力、峰值应力等参数。参数一般只能通过有限元分析等手段获取，不仅需要在前处理阶段建立3D几何模型，而且需要划分高质量网格，并进行多个管件网格的匹配和拼接。由于工程实际中的管系往往十分复杂，分析过程繁杂，且对人员素质要求高，因此，开展管道系统的有限元仿真系统性工具的研究开发具有重要意义。

国内较多学者对压力管道仿真技术开展了一系列研究。林远龙[5]以GB/T 19624—2004《在用含缺陷压力容器安全评定》为基础，利用Visual Basic.net编程语言开发了含缺陷压力容器、压力管道的安全评定系统；陈江等[6]结合API-1104开发了在役压力管道缺陷安全评定专家系统；马业华等[7]分析了油气管道并综合考虑管道表面体积型缺陷等影响因素，利用Ansys有限元分析软件的二次开发功能，开发了含缺陷管道安全评定系统;闫怀磊等[8-9]采用Ansys的APDL语言开发了含未焊透缺陷承压管道的剩余强度评价软件;何剑坤等[10]就地质灾害中的管道系统安全性问题，同样利用Ansys有限元APDL语言，基于VC++和Ansys的接口技术，开发了地质灾害下管道的安全评定有限元分析系统，该系统还实现了对地质灾害造成的悬空和漂移管道的参数化建模、有限元分析以及安全评定功能。

国外商业软件CAESAR 2采用一维管单元模型，能获取管道受力情况和变形结果，检查管道是否满足规范要求，用户可以自定义支吊架形式使其应力符合标准[11]，该软件不仅适用于工业复杂管系，还广泛应用于埋地管道分析[12]。但CAESAR 2 是完全采用管单元的一维分析软件，无法获取3D模型复杂的应力分布和开展应力线性化分析，因此,不能直接利用CAESAR 2结合GB/T 19624—2019完成含缺陷管道评价。

鉴于国内暂无用于分析复杂管系的仿真工具和国外软件CAESAR 2无法开展3D实体有限元分析的现状，本文提出一种先对管系整体布局；再局部建模；最后整体拼接的开发思路。先将整体复杂管系布局抽象成为一系列节点及其连接关系组成的图，表征拓扑连接关系和空间变换；然后，针对每一个组件提出参数化建模方案和网格划分方法，即将管系整体拼接和单个管道的建模完全分离；最后，利用空间变换和节点耦合方法组成3D管系，完成后续加载和计算。基于此思路开发的软件系统(命名为ABAP)实现全参数化建模，只需要修改几何参数即可实现整体管系的更新并开展后续参数分析，对复杂管系的安全评价具有较大的应用价值。

1 管道节点-连接关系图

以图1(a)中3D管系为例，该模型由直管、弯管、三通组成，不在一个平面内，点1,2之间为直管,点2,3之间为弯管，点3,4,5和点6，7，8对应三通，点5,6之间为直管，点8,9对应弯管。虽然3D管系空间结构看似复杂，但其拓扑关系和空间布置的信息可抽象为图1(b)所示的节点-连接关系图。

(a)

(b)

在图1(b)中，点1，2可表示图1(a)中直管1-2 的两个端部位置，点2,3和中间的虚拟节点(1)可确定弯管2-3的两个端部位置、弯管折角及空间朝向，三通管3-4-5可抽象为图1(b)中的点3,4,5和中间虚拟节点(2)，点5,6可表示直管5-6，以此类推。节点和虚拟节点的引入可唯一确定各个类型管件间的拓扑连接关系，若节点a和节点b属于同一个管件，则节点a，b用线段连接。根据图论(Graph theory)，抽象的节点-连接关系图是一个无向图，图节点表示各个管件关键节点的坐标和序号，线段表示单个管件内部的连接关系。节点-连接关系图可将3D管系的拓扑连接关系和具体单个管件的几何参数分离开，并分别加以处理，这一思路和主流3D建模软件的边界表示法[13]类似。

2 空间变换和网格划分及匹配

2.1 空间变换和管系组装

(1)

弯管节点2的外法线矢量为：

(2)

(3)

(4)

式中，·表示内积;×表示外积;‖·‖2表示欧式空间范数。

(a)三通建模坐标系

(b)坐标变换示意

2.2 管件几何建模和网格划分

2.2.1 直管、变径管

ABAP软件将参数分为实质参数和缺省参数。缺省参数可从节点-连接关系图中计算，如式(1)的弯管折角可由3个节点坐标唯一确定；实质参数表示不能从节点-连接关系图中获取几何参数，需要指定，如壁厚、外直径等。因此，直管在节点-连接关系图中以1号节点和2号节点两个节点表征，实质参数为外直径和壁厚，缺省参数为长度；变径管以端部1号节点和2号节点两个节点表征，可划分为开始直管段、中间过渡段和末尾直管段三部分，实质参数为开始段外直径、末尾段外直径、开始段壁厚、末尾段壁厚、开始段长度和中间过渡段长度，缺省参数为末尾段长度。

图3为变径管几何参数图。图4示出直管和变径管的网格划分方案，壁厚方向有6层网格，以保证弯曲应变的精度，同时在轴向和环向上设置均匀分布网格。

图3 变径管几何结构参数

图4 直管、变径管网格划分

2.2.2 弯管

弯管由两段直管和圆弧过渡段构成，图5示出弯管的几何参数，其中，弯管需要3个节点，即1,2,3号节点来表征，其中1,3号节点位于两个端部中心，2号节点对应两个直管中心线的交点。弯管实质参数为弯管外直径、壁厚和倒圆半径，缺省参数为直管段长度1、直管段长度2和折角。实际上，倒圆半径规定了过渡段的范围，由相切管系固定了过渡段的长度。图6示出了弯管网格划分方案。壁厚方向划分为6层网格，在直管段，环向和轴向均匀划分网格；而在过渡段，环向均分网格，轴向以扫掠方式生成网格。

图5 弯管几何结构参数

图6 弯管网格划分

2.2.3 三通

三通管是几何参数较为复杂的管件，也是工业管道系统中至关重要的管件。三通管的节点包括主管端部1号和2号节点，支管4号节点和虚拟3号节点。通常三通结构存在等壁厚过渡和等圆角过渡两种形式[14]，等壁厚过渡中，外圆角半径等于内圆角半径和壁厚之和；等圆角过渡中，内外圆角半径大小相等，因此有必要考虑过渡圆角的参数化。三通管实质参数包括主管外直径、直管外直径、主管壁厚、支管壁厚、过渡段外部倒圆半径和过渡段内部倒圆半径；缺省参数包括主管道长度1、主管道长度2和直管长度。内外倒圆表征主管和支管过渡区域，倒圆半径参数确定了过渡段的大小和壁厚分布(见图7)。在过渡段上，壁厚从主管壁厚逐渐连续过渡到支管壁厚，而越大的倒圆半径表征越大的过渡区，在此情况下，壁厚变化较为缓慢，而当倒圆半径很小时，壁厚剧烈变化。

图7 三通管几何结构参数

过渡区的存在使三通结构化有限元网格的划分尤为困难，根据支管外径与主管外径的比值，将三通分为两类：(1)主管和支管外直径相差不大(比值大于0.95)；(2)支管外直径明显小于主管外直径(比值小于0.95)。对于(1)类，由于主管和支管外直径相差不大(见图8)，因此，相贯线接近于45°平面和管道的交线，通过平移45°平面，生成带45°倾角的多个网格过渡区，对过渡区分别划分网格即可；对于(2)类，其网格生成方法比较复杂，先以内外倒圆对应的空间曲线构造两个直纹面，通过两个直纹面切割实体，再过三通中心点，作和直纹面最外侧点偏转5°而成的主管切分面，在主管上形成倒三角形区域，相较于传统划分方法，该方法保证了主管和支管环向网格数一致，使管件的匹配变得容易。由图9可以看出，三通结构化网格质量较好。

图8 斜切型三通结构化网格(适用于主管和接管直径接近的情况)

(a)模型分块

(b)结构性网格生成

2.3 网格匹配性问题

输入误差或结构固有特性均会引起网格不匹配的问题，即组装后有限元网格节点不一定完全匹配，ABAP软件在设计之初就考虑了此问题。如图10所示，两根直管存在微小的不共线偏折，若不加处理，两个端面上网格不能连续过渡。有些网格不匹配问题是管系本身性质造成的，如图11中管系由1根三通管和1根弯管组成，由于三通管和弯管端部均为均分网格，当弯管平面和三通管平面的夹角非0°和非90°时，网格一般无法匹配(扭转不一致型不匹配)。为此，本文提出了匹配网格的概念，匹配算法如下。

(1)设沿环向有m个节点，沿厚度方向有n个节点，数据结构上用2个矩阵Mij,Sij表示第i个环向节点和第j个厚度方向节点，其中1≤i≤m,1≤j≤n。

(3)沿着中心线方向，将原管道收缩微小的长度，即2节点收缩到4节点和5节点，储存4,5节点对应截面的节点坐标，中间空间用于匹配网格生成，具体如图12所示。

图10 两直管非共线型网格不匹配示意

(4)按照步骤(2)中所得的i0，连接所有的4,5节点对应截面上的最相邻节点，即将节点4截面的环向第1个节点和节点5对应截面的环向第i0节点直线连接，依次类推，一般地，将节点4截面环向第j个节点和节点5对应截面环向第[i0+j-1]个节点直线连接(其中，[i0+j-1]表示i0+j-1除以m的余数)。按此方法，对截面所有节点重复上述操作，连接全部生成的直线段，生成六面体网格。

图11 管道扭转不一致型网格不匹配示意

图12 匹配网格生成过程示意

上述算法对于微小不匹配的适用性很强，图13示出了图11中弯管和三通的匹配网格：生成的匹配网格通过扭转方式连续地过渡了两个管道的网格。

图13 三通管、弯管的匹配网格

3 截面载荷处理

ABAP软件可以设置任意端部截面上的总力(X,Y,Z方向)和总力矩(X,Y,Z方向)。总力矢量可分解为垂直于截面的分量和平行于截面的分量，并采用均匀分布的切应力和正应力分别等效平行分量和垂直分量。力矩载荷分解为弯矩载荷和扭矩载荷，采用垂直于截面的线性分布正应力等效弯矩载荷，采用沿环向的均匀切应力等效扭矩载荷。等效方法保证了截面总力和总力矩的正确性(见图14)。

图14 端部载荷等效方法示意

4 ABAP软件架构简介

ABAP软件模块架构见图15。

图15 ABAP软件模块架构

ABAP软件由人机交互主界面程序和Abaqus PDE二次开发模块组成，在主界面部分，节点编辑器用于输入节点-连接关系图的节点坐标和编号，连接关系编辑器用于指定单个管件中从属节点的编号，管道属性编辑器用于指定管件的实质参数。另外，利用OpenGL编写了模型可视化窗口，可用于即时预览管系模型，工作文件生成器生成传递到Abaqus PDE的自定义格式文件。当主界面生成工作文件后，将采用多线程模式调用Abaqus Python开发环境中的工作文件解释模块、3D管件建模模块、3D结构化网格划分和匹配模块、管道系统组装模块、管道载荷加载模块5个Python模块。其中，工作文件解释模块将主界面传来的工作文件做进一步语法解读和数据整理；3D管件建模模块根据实质参数和节点-连接关系图，生成三维几何模型；3D结构化网格划分和匹配模块将管件三维几何模型划分为结构化有限元网格，并在必要时，生成匹配网格；最后，由管道系统组装模块将管件和相关匹配网格按照空间变换拼接成为管道系统，并由管道载荷加载模块完成添加载荷、边界条件设置和作业提交的任务。

图16示出ABAP软件生成的由6个管件(1个三通、3个弯管、1个直管和1个变径管)组成管道系统的全过程，用户仅参与定义节点-连接关系图和指定实质参数，最终查看计算结果云图即可，整个过程仅需要几分钟即可完成，可明显提高仿真计算的效率。

图16 ABAP软件中模型生成和计算全过程(长度单位：m)

5 管系外压垮塌仿真案例

现有文献中多涉及单一直管或变径的外压垮塌分析[15-19]，缺乏对管系整体的外压垮塌分析研究。以图16中的管系为例，借助ABAP软件，研究壁厚对管系特征垮塌压力的影响，其中，三通支管、变径管小头和两者之间的连接弯管的壁厚t1、外径(D1=0.8 m)分别相等，三通主管、直管、右侧两个弯头、变径管大头的壁厚t2=0.1 m、外径D2=1.0 m分别相等。参数分析中，材料弹性模量为200 GPa，泊松比0.3，其他参数不变，仅改变t1数值，依次计算了t1=0.05,0.06,0.07,0.08 m四种情况下，管系在外压下的特征值屈曲的模态和垮塌压力，计算结果如图17所示。垮塌压力随着t1增大而增大，分别为197.2，314.0，475.4，625.5 MPa。当t1=0.05 m时，屈曲模态中，变形集中在三通支管，其余管件未发生明显变形;当t1=0.08 m时，直管、变径的变形较大，这说明管系的外压垮塌压力取决于较为薄弱的区域，管系一般不会同时发生垮塌。采用圆柱壳外压垮塌公式DTMB[20](见式(5))、von Mises公式[20](见式(6))、平面应变公式[16](见式(7))验证FEA垮塌压力(pcr)的合理性。根据节点-连接关系图(见图16)，壁厚t=t1的管道长度约为2.5 m(公式计算中薄弱区域等效为圆柱壳，但实际薄弱区域包含三通接管和弯头，因此2.5 m为粗略近似值)，取L=2.5 m，t=t1，a=D1/2-t1，E=200 GPa，υ=0.3，计算结果见图18。可以看出，平面应变公式计算的数值明显小于FEA数值，而考虑了长度L影响的DTMB和von Mises公式的计算结果与FEA接近，验证了FEA计算结果的合理性。

(5)