基于Python的SAP2000向ANSYS模型转化技术及软件开发研究

马小平

(中铁第一勘察设计院集团有限公司,西安 710043)

引言

随着我国经济实力与科学技术的迅猛发展,一系列新颖复杂的建筑结构不断涌现,以满足人民生活的需求。为保证此类复杂建筑结构的可靠性,工程人员常利用一些成熟的商业分析软件对结构进行设计和校核。然而,结构设计是一个不断反复调整以获得最优结果的过程,特别是对于复杂的大跨度或高层建筑,除常规的弹性分析之外,非线性分析也必不可少,例如:弹塑性时程分析、静力稳定分析、多尺度受力分析等[1-3]。因此,仅依靠单一的设计软件或结构分析软件,难以实现对复杂结构多类型力学响应的综合分析。

工程人员为验证设计结构的正确性和可靠性,常采用多种类型的软件进行比较分析。除常规的设计软件(如:PKPM和YJK)之外,SAP2000、ETABS等结构专业设计分析软件,因其良好的三维结构分析能力、便捷的空间建模能力以及完善的荷载计算功能,受到了业内专业人士的广泛应用[4-6]。同时,ANSYS、ABAQUS等大型通用有限元分析软件,因其超强的拓展性和灵活性,加上内嵌精确的数值计算方法,故在特殊构件、复杂结构的分析验算中得到使用[7-10]。然而,仅依靠ANSYS、ABAQUS自带的前处理功能对较为复杂的结构模型进行重构,不仅耗时费力,且意义不大。考虑各类型软件模型建立方法上的相似性,故可编译软件之间数据转换的接口程序,将人工建模的工作量转移到程序完成,极大地提高建模效率和准确性[11-13]。

结构设计分析软件与通用有限元软件之间并不能直接进行数据交换共享,因此,众多企业和科研团队基于各类开发平台推出了一系列模型转化程序。北京盈建科软件股份有限公司已相继开发了YJK和ETABS接口软件(YJK-ETABS)、YJK和MIDAS接口软件(YJK-MIDAS)、YJK和SAP2000接口软件(YJK-SAP2000)、YJK和ABAQUS接口软件(YJK-ABAQUS)等,为YJK设计模型向各类型结构分析软件的转化提供了便捷。葛金刚等[14]基于Python语言,开发了将结构模型由SAP2000及MIDAS向ABAQUS转换的接口程序,并以天津某展览馆屋盖结构为例,对比验证程序转化结果的正确性。孟仲永[15]以AutoCAD作为图形处理平台,利用ObjectArx及VS2008和C++语言开发了从SAP2000到ABAQUS的模型转换程序,并应用到了两个实际大跨度复杂空间结构的模型转化中,验证了转化结果可靠性。在此基础上,王杰[16]通过类似的方法研发了MIDAS/GEN到ABAQUS的模型转换程序,实现了结构楼层信息、节点及单元信息、截面及材料信息、荷载及约束信息的高效转化。基于SAP2000和ABAQUS软件之间的内在转换逻辑,张月强等[17]编制了用以模型自动转换的MTR1.0程序,实现了SAP2000几何模型向ABAQUS中的精确转化。基于此,祝辉庆等[18]也利用MATLAB软件编制SAP2000到ABAQUS的模型数据转换程序,实现了节点、单元、材料、截面、荷载、约束等信息的自动转化。曹伟良等[19]采用VB6.0结合ACCESS数据库程序技术,将SAP2000导出的MDB文件转化为ANSYS命令流文件,实现了SAP2000常规构件及简单荷载类型的转化。

通过对现有模型转化技术及软件开发研究的深入调研,可以看出,众学者分享了不同软件之间的转化思想以及相应程序平台的开发,为相关研究的进一步发展提供了有力的技术支撑。然而,针对结构设计分析软件SAP2000与通用有限元软件ANSYS之间的模型转化,其相关研究整体上仍较为匮乏,且存在以下重要问题亟待改进。

(1)现有SAP2000向ANSYS的模型转化程序,其转化内容过于简单。部分自编程序仅支持几何模型的转化,相对成熟的程序虽然转化内容较为全面,但也仅支持简单的低阶单元(如Beam44梁单元、Shell63壳单元等)、常规截面类型(如矩形、工字形等)、常规荷载(均布荷载、集中荷载等)。当需要依赖高阶单元(Beam188、Shell181等)确保非线性分析的准确性时,或结构包含特殊截面(如变截面、自定义截面等)和荷载类型(三角形分布荷载、梯形荷载等),上述转化程序就会失效。

(2)现有SAP2000向ANSYS的模型转化程序,其功能性严重不足。对于实际结构分析中常用到连接件设置、节点坐标系变换、梁单元坐标系变换、梁端自由度释放、刚性域等功能,现有转化程序仍然不够完善。此外,由于未考虑控制网格划分数量的功能,对计算精度或效率有要求的模型,现有转化程序难以提供有效帮助。

(3)现有SAP2000向ANSYS的模型转化程序,其设计思想过于简单。现有转化程序大多以ANSYS中的节点(Node)代替SAP2000中结点(Joint),直接建立构件或结构的有限元模型。该方法虽然能够快速实现模型的创建,但对于有网格划分要求或荷载施加形式复杂的模型,以上述思想设计的程序无法实现模型的准确转化。

鉴于此,开发团队确定了以关键点作为模型转化基础,并考虑了多种单元类型、材料本构、荷载形式的转化需求,研究了SAP2000模型结构信息在ANSYS中的实现条件及方法,进而开发了“SAP2000 To ANSYS模型转化软件 V1.0”(简称“STAMT V1.0”),以期实现SAP2000结构分析模型向ANSYS有限元仿真模型的全方面转化,为结构设计的从业人员提供必要的技术支撑。

1 SAP2000 To ANSYS模型转化程序

1.1 STAMT程序的目标及功能

STAMT程序的核心目标是依据SAP2000导出的.s2k文件,自动读取文件内数据,按照指定的网格划分数,转化模型的全部信息,生成ANSYS有限元软件可读取的参数化命令流,并以文档的形式导出,实现SAP2000结构分析模型向ANSYS有限元分析模型的精准转化。

STAMT程序具备了以下几个主要功能。

(1)SAP2000文件的导入和数据读取功能,如:材料信息读取、几何坐标读取、构件截面读取、荷载分布读取等。

(2)自定义框架梁柱单元的网格划分数量,即用户可根据需求指定梁单元网格划分的数量。

(3)多种结构单元类型及基本材料属性的转化,如:梁单元、板壳单元、弹簧单元、质量单元、线弹性材料、弹塑性材料等。

(4)多种结构构件功能性需求的转化,如:节点坐标系变换、梁单元坐标系变换、梁端自由度释放、创建刚性域等。

(5)多种荷载工况的转化,如:节点力、节点位移、梁柱分布载荷、梁柱集中载荷、板壳均布载荷、传递到梁的板壳载荷等。

STAMT程序面向的用户主要是土木工程领域相关的结构设计及分析人员。因此,程序应具有良好的用户交互界面,方便不熟悉ANSYS参数化设计语言(ANSYS Parametric Design Language,APDL)的人员进行系统应用,且转化结果清晰直观,可供相关人员快速、高效地实现模型的转化,节省模型重建所耗费的时间,大大提高工程人员的结构验核效率。

1.2 STAMT程序的开发环境

为满足上述目标及功能,STAMT程序需具备良好的用户交互界面,快速实现文件的导入导出,并具备模型数据的批量转化和数值计算功能。Python作为一款解释性、编译性、面向对象的脚本语言,不仅具有易于学习、阅读、维护等特点,且包含了各类的标准库,可实现不同文本数据的读入和写出,以及各类函数的数值运算[20-21]。同时,Python中的PyQt5库涵盖了丰富的功能函数用于交互界面设计,借助Eric6和Qt Designer可更加快速、便捷地实现软件操作界面的创建。因此,本软件采用Python 3.9作为开发环境,借助Eric6集成开发软件,采用用户界面与业务逻辑分离思想进行软件的整体架构设计。

1.3 STAMT程序的关键技术及流程

STAMT程序内部的模型转化过程主要包括读取SAP2000文件,生成并导出相关的APDL代码,如定义单元类型、材料属性、实常数、截面属性等,创建关键点、几何对象、单元等,变换单元坐标轴及施加约束,施加荷载及质量等。具体转化流程如图1所示。

图1 STAMT程序内部模型转化流程Fig.1 Model transformation process of the STAMT program

1.3.1 读取SAP2000文件

SAP2000软件可导出.s2k格式文件,该文件涵盖了SAP2000模型的所有数据库表格,包含了所有在交互界面中设置的模型信息。STAMT程序会首先读取该文件,将数据分类并以列表或字典格式存储。

1.3.2 单元类型的转化

为满足更高的计算需求,STAMT程序提供了梁单元(Beam188)、壳单元(Shell181)、弹簧单元(Combin14及Combin39)、质量单元(Mass21)、网格单元(Mesh200)的转化条件。其中,壳单元通过关键项KEYOPT(1)分为两类,即考虑薄膜及弯曲刚度的壳单元,其对应SAP2000中的壳单元,以及仅考虑薄膜刚度的壳单元,其对应SAP2000中的膜单元。弹簧单元包含线性弹簧单元(Combin14)和非线性弹簧单元(Combin39)。通过Combin14的关键项KEYOPT(2),分别定义了6个自由度方向上的线性弹簧单元。通过Combin39的关键项KEYOPT(1)和KEYOPT(2),在6个自由度方向上分别定义了弹性非线性弹簧单元以及弹塑性非线性弹簧单元。以弹簧单元为例,首先,建立转化弹簧单元的子函数。然后,依次编写一维线性弹簧、一维弹性非线性弹簧、一维弹塑性非线性弹簧的转化代码。最后,将列表返回并调用子函数。

def Spr_Elem_Type():

APDL_Spr_Elem_Type+=['et,2,combin14'+' '+'keyopt,2,1,0'+' '+'keyopt,2,2,1'+' ']

APDL_Spr_Elem_Type+=['et,8,combin39'+' '+'keyopt,8,1,0'+' '+'keyopt,8,2,0'+' '+'keyopt,8,3,1'+' ']

APDL_Spr_Elem_Type+=['et,14,combin39'+' '+'keyopt,14,1,1'+' '+'keyopt,14,2,0'+' '+'keyopt,14,3,1'+' ']

return APDL_Spr_Elem_Type

1.3.3 材料属性的转化

为使转化后模型的本构关系更具通用性,STAMT程序提供了多种基本材料属性的转化条件,包括:弹性、理想弹塑性、双线性弹塑性、多线性弹塑性等。

1.3.4 实常数的转化

STAMT程序提供了梁单元、壳单元、弹簧单元、质量单元的实常数转化条件。梁单元、壳单元仅创建了实常数编号,以满足网格划分的需要,无其他实际意义。弹簧单元的实常数包含了线性弹簧单元实常数、非线性弹性弹簧单元实常数、非线性塑性弹簧单元实常数三部分,具体设置参数和对应关系如表1所示。

表1 弹簧单元的对应关系及实常数设置Table 1 Correspondence of spring units and real constants

质量单元的实常数可设置3个平动方向的质量以及3个转动方向的质量惯性矩。确定了模型信息中实常数的对应关系,以及整体的转化逻辑,可编译相关代码实现实常数信息的转化。

1.3.5 截面属性的转化

STAMT程序提供了梁壳单元截面属性定义的转化条件。壳单元截面仅需设置板壳的厚度。梁单元截面可转化的类型有矩形、工字形、箱形、圆管、自定义截面和变截面。其中,自定义截面和变截面的转换较为特殊,现有转化程序都未涉及,而这两种类型在实际工程中应用广泛。为此,结合ANSYS中有关自定义梁截面和变截面梁的APDL命令,编译可转化此类截面的相关代码。

ANSYS中自定义截面的APDL命令如下:

APDL_Beam_Sec+=['sectype,'+Frame_Sec_Num+',beam,asec,'+Frame_Sec_Name+' '+'secdata,'+ASEC_Area[1]+','+ASEC_Iyy[1]+','+ASEC_Iyz[1]+','+ASEC_Izz[1]+','+ASEC_Iw+','+ASEC_Jt[1]+','+ASEC_CGy+','+ASEC_CGz+','+ASEC_SHy+','+ASEC_SHz+','+ASEC_TKz+','+ASEC_TKy+','+ASEC_TSxz+','+ASEC_TSxy]

上述代码中,Frame_Sec_Num为截面编号;Frame_Sec_Name为截面名称;ASEC_Area为自定义截面面积;ASEC_Iyy为自定义截面对y轴的惯性矩;ASEC_Iyz为自定义截面的惯性积;ASEC_Izz为自定义截面对z轴的惯性矩;ASEC_Jt为自定义截面的扭转惯性矩。上述参数在SAP2000导出的数据文件中均有对应信息,可直接转化得到。然而,自定义截面的翘曲惯性矩ASEC_Iw、截面重心y坐标ASEC_CGy、截面重心z坐标ASEC_CGz、截面剪切中心y坐标ASEC_SHy、截面剪切中心z坐标ASEC_SHz、截面沿z轴厚度ASEC_TKz、截面沿y轴厚度ASEC_TKy、截面xz剪切修正系数ASEC_TSxz、截面xy剪切修正系数ASEC_TSxy,无法从SAP2000文件中直接获取并转化得到。因此,程序在输出的文件中保留有参数符号及说明,使用者可利用其他截面计算软件求得上述参数后,填写入命令流中进行仿真计算。

ANSYS中变截面的APDL命令如下:

APDL_Frame_Sec+=['sectype,'+Frame_Sec_Num+',taper,,'+Frame_Sec_Name+' '+'secdata,'+Sec_Num_Start+','+Joint_I_X+','+Joint_I_Y+','+Joint_I_Z+' '+'secdata,'+Sec_Num_End+','+Joint_J_X+','+Joint_J_Y+','+Joint_J_Z]

如上述代码所示,首先,需要确定变截面梁的I节点(Frame_Joint_I)和J节点(Frame_Joint_J)的编号,并根据节点坐标数据(Joint_Coor)分别遍历得到I、J节点的坐标。然后,在截面数据中,确定梁始端的截面编号(Sec_Num_Start)及其对应I节点坐标,梁末端的截面编号(Sec_Num_End)及其对应J节点坐标。通过以上步骤,即可转化得到定义梁变截面属性的APDL命令流。

1.3.6 关键点、几何对象、单元划分的转化

利用SAP2000文件中的“连接数据”,通过编译转化代码,可实现几何模型(包括点、线、面)的创建。此处代码无特殊难点,故不再赘述。STAMT程序提供了梁单元划分、壳单元划分、弹簧单元创建的转化能力。其中,梁单元划分的数量可根据需要输入,壳单元按默认方式划分。弹簧单元创建的代码较为复杂,此处将详细阐述。

首先,在连接(Link)位置上创建几何线,选择线并赋予空的材料属性和实常数,以及网格单元类型Mesh200,并划分一个单元。

APDL_Spr_Elem+=['l,'+Link_Joint_I+','+Link_Joint_J]

APDL_Spr_Elem+=['lsel,s,line,,Line_Num_Max+1']

APDL_Spr_Elem+=['latt,Mat_Num_Max+1,Real_Num_Max+1,23']

APDL_Spr_Elem+=['lesize,all,,,1']

APDL_Spr_Elem+=['lmesh,all']

以x方向一维线性弹簧为例,确定单元类型和实常数后,即可根据连接(Link)节点I和J的坐标创建弹簧单元。其他单元类似。

APDL_Spr_Elem+=['type,2']

APDL_Spr_Elem+=['real,'+Real_Link_Linear]

APDL_Spr_Elem+=['en,Elem_Num_Max+1,node(kx('+Link_Joint_I+')'+','+'ky('+Link_Joint_I+')'+','+'kz('+Link_Joint_I+')'+')'+',node('+'kx('+Link_Joint_J+')'+','+'ky('+Link_Joint_J+')'+','+'kz('+Link_Joint_J+')'+')'+' ']

1.3.7 坐标轴变化及约束的转化

STAMT程序提供了节点坐标轴变化、梁端自由度释放、梁单元坐标轴变化、创建刚性域、施加支座约束等功能的转化条件。依据SAP2000文件中“节点局部轴指定1-标准”“框架释放指定1-通用”“节点约束指定”“节点刚性约束指定”“框架局部轴指定1-标准”等模块的数据,即可编译代码实现相应模型信息的转化。节点坐标轴变换、梁端自由度释放、施加支座约束分别通过APDL中的nmodif命令、endrelease命令、d命令等实现,其代码较为简单,此处不再赘述。此处,重点介绍梁单元坐标轴变化及刚性域创建。

梁单元坐标轴变化的部分代码如下。

首先,确定梁单元坐标系x轴的方向向量,并计算向量的模和归一化后的单位向量。

Frame_x_Axis_Vec=np.array([(Joint_J_X-Joint_I_X),(Joint_J_Y-Joint_I_Y),(Joint_J_Z-Joint_I_Z)])

Frame_x_Axis_Vec_Len=np.linalg.norm(Frame_x_Axis_Vec)

Frame_x_Axis_Vec_Norm=Frame_x_Axis_Vec/Frame_x_Axis_Vec_Len

然后,将x轴的单位方向向量与整体坐标系下Z轴的方向向量作比较,若二者平行,则单元坐标系的y轴方向与整体坐标系下Y轴的方向一致,若不平行,则通过二者叉积确定。单元坐标系z轴的方向则通过单元坐标系的y轴单位方向向量与x轴单位方向向量的叉积得到。

Global_Z_Axis=np.array([0,0,1])

if Frame_x_Axis_Vec_Norm.dot(Global_Z_Axis)/(np.sqrt(Frame_x_Axis_Vec_Norm.dot(Frame_x_Axis_Vec_Norm))*np.sqrt(Global_Z_Axis.dot(Global_Z_Axis)))==1 or Frame_x_Axis_Vec_Norm.dot(Global_Z_Axis)/(np.sqrt(Frame_x_Axis_Vec_Norm.dot(Frame_x_Axis_Vec_Norm))*np.sqrt(Global_Z_Axis.dot(Global_Z_Axis)))==-1:

Frame_y_Axis_Vec=np.array([0,1,0])

Frame_z_Axis_Vec=np.cross(Frame_x_Axis_Vec_Norm,Frame_y_Axis_Vec)

else:

Frame_y_Axis_Vec=np.cross(Global_Z_Axis,Frame_x_Axis_Vec_Norm)

Frame_z_Axis_Vec=np.cross(Frame_x_Axis_Vec_Norm,Frame_y_Axis_Vec)

接着,确定旋转中心轴(Rot_Axis)为单元坐标轴x,旋转角度为Rot_Ang,起始向量(Vec_Origin)为单元坐标轴z,则可计算得到目标向量(Vec_Target),进而将目标向量的元素与梁单元起始节点(I节点)的坐标对应相加,即可得到旋转后梁单元的方向节点。

Rot_Axis=Frame_x_Axis_Vec_Norm

Vec_Origin=Frame_z_Axis_Vec

Vec_Target=Vec_Origin*math.cos(Rot_Ang)+np.cross(Rot_Axis,Vec_Origin)*math.sin(Rot_Ang)+Rot_Axis*(np.dot(Rot_Axis,Vec_Origin))*(1-math.cos(Rot_Ang))

New_Frame_Z_Axis_X_Coor=Vec_Target[0]+Joint_I_X

New_Frame_Z_Axis_Y_Coor=Vec_Target[1]+Joint_I_Y

New_Frame_Z_Axis_Z_Coor=Vec_Target[2]+Joint_I_Z

最后,通过创建新的方向控制节点,并利用emodif命令修改梁单元的属性,实现梁单元坐标轴的变化。

APDL_Beam_Axis+=['allsel,all']

APDL_Beam_Axis+=['n,'+'Node_Num_Max+'+str(i+1)+','+str(New_Frame_Z_Axis_X_Coor)+','+str(New_Frame_Z_Axis_Y_Coor)+','+str(New_Frame_Z_Axis_Z_Coor)]

APDL_Beam_Axis+=['lsel,s,line,,'+Frame_Line_Name]

APDL_Beam_Axis+=['esll,s']

APDL_Beam_Axis+=['emodif,all,-3,Node_Num_Max+'+str(i+1)+' ']

刚性域创建的部分代码如下。首先,根据节点约束信息(Constraint_Joint),明确需要创建刚性体约束的节点,将其筛选出来并定义为一个组件(Node_Rigid_Region)。

for k in range(0,len(Constraint_Joint)):

if k==0:

APDL_Rigid_Region+=['ksel,s,kp,,'+Constraint_Joint[k]]

elif k>0:

APDL_Rigid_Region+=['ksel,a,kp,,'+Constraint_Joint[k]]

APDL_Rigid_Region+=['nslk,s']

APDL_Rigid_Region+=['cm,Node_Rigid_Region'+str(i+1)+',node']

接着,从每个刚性体约束节点组中选第一个节点作为主节点(Node_Main),并利用cerig命令在主结构与其他节点之间建立刚性约束。

APDL_Rigid_Region+=['Node_Main='+'node('+'kx('+Constraint_Joint[0]+'),ky('+Constraint_Joint[0]+'),kz('+Constraint_Joint[0]+')'+')']

APDL_Rigid_Region+=['cmsel,s,Node_Rigid_Region'+str(i+1)+',node']

APDL_Rigid_Region+=['cerig,Node_Main,all,'+'ux']

1.3.8 荷载及质量的转化

STAMT程序提供了节点质量、节点荷载、支座位移、梁柱分布荷载、梁柱集中荷载、楼面均布荷载、楼面传递至梁的荷载等转化条件。节点质量、节点荷载、支座位移的转化较为简单,此处不再详述。下面以梁柱荷载和楼面荷载的转化作为关键技术点,进行展开介绍。

以梁柱分布荷载的转化为例,由于对模型进行了网格划分,因此框架梁上的分布荷载并不能直接通过sfbeam命令进行施加。鉴于此,分4个步骤编译了梁单元分布荷载施加的转化代码。第一步,选取需要施加分布荷载的框架梁(Frame),确定其起始节点的几何坐标以及梁长。第二步,基于Python语言,编译APDL命令流,命令流内容为选取ANSYS模型中梁(Beam)上所有节点,将其录入至所创建列表中。第三步,与第二步方法一致,在ANSYS中遍历每个梁上的节点,并将其定义为节点组件。以上步骤转化代码不再详述。第四步,读取分布荷载数据,判断荷载方向,确定其与单元坐标轴的夹角,最后施加于梁单元上。此处,着重阐述第四步的实现方法。

以单元坐标轴x和荷载沿重力方向为例,其他方向类似。首先,判断框架分布荷载的方向,定义荷载的方向向量。

if Frame_Distr_Load_Dir=='Gravity’ or Frame_Distr_Load_Dir=='Gravity Proiected':

Load_Vec=np.array([0,0,-1])

然后,判断梁单元坐标方向与荷载方向的夹角。单元坐标轴的方向向量(Frame_x_Axis_Vec)已在前文说明。

Frame_x_Axis_Vec_Len=np.sqrt(Frame_x_Axis_Vec.dot(Frame_x_Axis_Vec))

Cos_Angle_x=Frame_x_Axis_Vec.dot(Load_Vec)/(Frame_x_Axis_Vec_Len*Load_Vec_Len)

Angle_x=np.arccos(Cos_Angle_x)*180/np.pi

进而,判断框架分布荷载是否为投影荷载,若为投影荷载,则根据SAP2000程序的规定,将框架分布荷载乘以荷载与坐标轴夹角的正弦值。

if Frame_Distr_Load_Dir=='Gravity Proiected':

Frame_Distr_Load_A=str(float(Frame_Distr_Load_A)*math.sin(Angle_x))

Frame_Distr_Load_B=str(float(Frame_Distr_Load_B)*math.sin(Angle_x))

最后,将框架分布荷载在单元坐标轴x上投影,并根据分布荷载在框架梁上的作用距离,筛选梁上节点,分段将荷载施加到梁单元上。

Load_a_x=float(Frame_Distr_Load_A)*math.cos(Angle_x)

Load_b_x=float(Frame_Distr_Load_B)*math.cos(Angle_x)

APDL_Beam_Distr_Load+=['Dis_I='+str(Len_as)+'-(i-1)*'+str(Elem_Size)]

APDL_Beam_Distr_Load+=['Dis_J=(i)*'+str(Elem_Size)+'-'+str(Len_bs)]

APDL_Beam_Distr_Load+=["node_1=strcat('Node_',chrval(i))"]

APDL_Beam_Distr_Load+=["node_2=strcat('Node_',chrval(i+1))"]

APDL_Beam_Distr_Load+=['cmsel,s,node_1,node']

APDL_Beam_Distr_Load+=['cmsel,a,node_2,node']

APDL_Beam_Distr_Load+=['esln,s,1,corner']

APDL_Beam_Distr_Load+=['sfbeam,all,3,pres,'+str(Load_a_x)+','+str(Load_b_x)+',,,Dis_I,Dis_J']

壳单元均布荷载的施加方法与梁单元相似,重点在于单元坐标系z轴的确定。图2为默认情况下,壳单元的单元坐标系。e1、e2、e3分别为单元坐标1轴、2轴、3轴;S1和S2为第一和第二参考方向;关键项KEYOPT(11)为单元x轴的默认方向,若KEYOPT(11)=0,则单元x轴的默认方向与质心的第一个参数方向S1一致。

图2 ANSYS中Shell181单元的默认单元坐标系Fig.2 Element coordinate system of Shell181 in ANSYS

参考方向S1和S2与单元的形函数有关。根据文献[22],参考方向S1的计算方法如下

(1)

(2)

其中,{x}I、{x}J、{x}K、{x}L为总体坐标系下节点的坐标。

默认单元坐标轴e3由下式计算得到

e3=S1×S2/|S1×S2|

(3)

面划分网格后,壳单元四个节点(IJKL)仅能在ANSYS程序中提取分析。因此,参考方向与单元坐标的确定通过APDL命令流实现,其转化代码如下。首先,确定施加荷载的面(Area_Name),挑选其上编号最大的壳单元,定义四个节点的数组以分别存放节点坐标,此处以I节点为例。

APDL_Shell_Unif_Load+=['*dim,Area'+Area_Name+'_NodeI_Vec,array,3,1']

APDL_Shell_Unif_Load+=['Area'+Area_Name+'_NodeI_Vec(1,1)=nx(Area'+Area_Name+'_NodeI)']

APDL_Shell_Unif_Load+=['Area'+Area_Name+'_NodeI_Vec(2,1)=ny(Area'+Area_Name+'_NodeI)']

APDL_Shell_Unif_Load+=['Area'+Area_Name+'_NodeI_Vec(3,1)=nz(Area'+Area_Name+'_NodeI)']

然后,根据式(1)和式(2),计算得到壳单元的参考方向S1,S2与之相似。

APDL_Shell_Unif_Load+=['*voper,Area'+Area_Name+'_Parr1,Area'+Area_Name+'_NodeJ_Vec,sub,Area'+Area_Name+'_NodeI_Vec']

APDL_Shell_Unif_Load+=['*voper,Area'+Area_Name+'_Parr2,Area'+Area_Name+'_NodeK_Vec,sub,Area'+Area_Name+'_NodeL_Vec']

APDL_Shell_Unif_Load+=['*voper,Area'+Area_Name+'_Parr3,Area'+Area_Name+'_Parr1,add,Area'+Area_Name+'_Parr2']

APDL_Shell_Unif_Load+=['*voper,Area'+Area_Name+'_S1_A,Area'+Area_Name+'_Parr3,mult,1/4']

APDL_Shell_Unif_Load+=['Area'+Area_Name+'_S1_B=(Area'+Area_Name+'_S1_A(1,1)**2+Area'+Area_Name+'_S1_A(2,1)**2+Area'+Area_Name+'_S1_A(3,1)**2)**(1/2)']

APDL_Shell_Unif_Load+=['*voper,Area'+Area_Name+'_S1,Area'+Area_Name+'_S1_A,div,Area'+Area_Name+'_S1_B']

最后,根据式(3),计算得到壳单元坐标轴e3。得到壳单元的面外法向坐标轴,即可通过与梁单元相似的方法,编译施加壳单元均布荷载的代码。

APDL_Shell_Unif_Load+=['*voper,Area'+Area_Name+'_Parr7,Area'+Area_Name+'_S1,cross,Area'+Area_Name+'_S2']

APDL_Shell_Unif_Load+=['Area'+Area_Name+'_Parr8=(Area'+Area_Name+'_Parr7(1,1)**2+Area'+Area_Name+'_Parr7(1,2)**2+Area'+Area_Name+'_Parr7(1,3)**2)**(1/2)']

APDL_Shell_Unif_Load+=['*voper,Area'+Area_Name+'_Parr9,Area'+Area_Name+'_Parr9,div,Area'+Area_Name+'_Parr8']

APDL_Shell_Unif_Load+=['*dim,Area'+Area_Name+'_esys3,array,1,3']

APDL_Shell_Unif_Load+=['Area'+Area_Name+'_esys3(1,1)=Area'+Area_Name+'_Parr9(1,1)']

APDL_Shell_Unif_Load+=['Area'+Area_Name+'_esys3(1,2)=Area'+Area_Name+'_Parr9(2,1)']

APDL_Shell_Unif_Load+=['Area'+Area_Name+'_esys3(1,3)=Area'+Area_Name+'_Parr9(3,1)']

2 STAMT程序算例验证

2.1 STAMT程序操作界面

STAMT程序的操作界面包括转化程序(图3(a))和使用说明(图3(b))两个基本模块。

图3 STAMT程序操作界面Fig.3 STAMT program operation interface

转化程序模块提供了梁单元划分数的输入、SAP2000文件录入、APDL文件输出、转化进度显示等功能。使用时,根据需要输入梁单元网格划分个数,并点击“浏览”按钮选取转化和导出的文件,或在输入框中自定义工作目录。完成上述步骤后,点击“转化”按钮,就会出现转化进度条,出现“转化完成”时,代表模型转化结束。

2.2 算例一:单层工业厂房

厂房基本布局为:纵向3跨,单跨4 m;横向1跨,跨度12 m;屋架下弦高度6 m;檐口高度7.5 m;屋脊高度8.5 m。柱采用箱形钢截面,梁采用工字形钢截面,桁架采用工字形钢桁架,屋面板为混凝土薄壳板。该结构的SAP2000模型如图4(a)所示,转换后的ANSYS模型如图4(b)。

图4 单层工业厂房模型Fig.4 Single-story industrial plant model

通过计算,对比两个软件的模态变形,如图5、图6所示。结果表明,SAP2000计算的前3阶模态与ANSYS的计算结果基本一致。

图5 未释放自由度单层工业厂房模型的模态对比Fig.5 Modal comparison of single-storey industrial models with unreleased degree of freedom

图6 释放自由度后单层工业厂房模型的模态对比Fig.6 Modal comparison of single-storey industrial models with released degree of freedom

为更加准确地验证STAMT程序的转化效果,分别按释放桁架自由度和不释放桁架自由度两种情况,对比SAP2000模型和ANSYS模型的质量和周期,如表2所示。可以看出,SAP2000模型的总质量和ANSYS模型一致。对于不释放自由度的模型,SAP2000的计算结果与ANSYS程序相近,误差平均值为2.44%,最大误差不超过5%。对于释放自由度后的模型,SAP2000计算的前3阶周期与ANSYS相近。

表2 单层工业厂房周期对比Table 2 Period comparison of single-storey industrial plant

但在较大周期(>4阶)时,二者的计算误差增加,这是由于SAP2000软件使用的是节点集中质量,而ANSYS软件在释放梁端自由度之后,模态计算无法打开集中质量开关(lumpm, on),从而造成二者的计算结果出现差异,且误差在高阶模态下不断增加。

整体而言,通过STAMT程序可以实现SAP2000模型向ANSYS模型的准确转化,特别是低阶模态的计算误差不超过5%。然而,对于释放自由度后,集中质量与一致质量的转化问题,仍需在今后的软件优化中进一步研究。

2.3 算例二:多层框架结构

框架结构长30 m,宽18 m,高18 m,结构由下部钢筋混凝土构件组成。梁柱构件截面种类为矩形或方形。整体框架结构的SAP2000模型如图7(a)所示。通过STAMT程序转化SAP2000模型,生成APDL命令流,将命令流导入ANSYS中,创建的ANSYS有限元模型如图7(b)所示。

图7 多层框架结构模型Fig.7 Multi-storey framework structural model

对比两个模型的模态变形,如图8所示。可以看出,SAP2000计算的模态变形图与ANSYS的计算结果一致,表明转化后的模型满足要求。

图8 多层框架结构模型的模态对比Fig.8 Modal comparison of multi-storey frame structures

进一步比较二者的质量以及计算的前5阶周期,如表3所示。结果表明,两个模型的结构总质量一致,而周期的误差在2%以内,且平均误差为1.32%,说明转化模型的计算精度满足要求。

表3 多层框架结构质量及周期对比Table 3 Comparison of weight and period for multi-storey frame structures

通过以上的算例分析,证明了本文所编译开发的STAMT程序可以快速、准确地将SAP2000模型转化为ANSYS有限元分析模型,可满足相关技术人员对模型转化的需求。

3 结论与展望

基于Python语言,编译并开发了STAMT程序,实现了SAP2000结构设计分析模型向ANSYS有限元模型的转化,并通过两个算例的计算验证了该程序的可行性及准确性。相比于现有的模型转化程序,本文开发的STAMT程序具有以下优点。

(1)丰富了单元类型、材料属性、截面形式、荷载形式的转化类型,提高了模型转化的精确性。

(2)涵盖了节点坐标系变换、梁单元坐标系变换、梁端自由度释放、创建刚性域等必备功能,进一步满足了模型转化过程中的功能性需求。

(3)以ANSYS的关键点(Keypoint)对应SAP2000的结点(Joint),使用先建立几何模型后生成有限元模型的思路,实现了梁单元网格数量的自定义设置以及板壳荷载向梁的传递。

开发的STAMT程序可为相关从业人员更加快速、准确地实现模型信息的转化,为结构分析结果的对比验证提供便捷。但由于结构设计分析软件与通用有限元软件之间的差异性,难以实现SAP2000模型向ANSYS模型的完全转化,STAMT程序中存在的不足会随着各分析软件的不断更新以及开发者的修复改进,而在将来的研究中予以完善。