基于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的模型轉化程序,其設計思想過于簡單?，F有轉化程序大多以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程序中存在的不足會隨著各分析軟件的不斷更新以及開發者的修復改進,而在將來的研究中予以完善。