馬斌,谷翠軍,徐先偉,黃艷麗
(中車青島四方機車車輛研究所 重慶研發中心,重慶 401133)
近年來,隨著國民經濟的發展,鐵路車輛的速度越來越快、運量越來越多,這對列車結構提出更高要求。列車結構的損壞大多在焊接結構上,其破壞形式一般為疲勞斷裂,故迫切需要找到一種準確計算焊縫疲勞壽命的方法。
計算焊縫疲勞壽命時,初始輸入值一般為靜強度計算中焊縫結構的節點力。有限元計算分析軟件Ansys具有強大的強度分析能力,其輸出的靜強度計算的節點力信息包含在rst文件中。目前,對rst文件數據結構的解析較少。本文對rst文件的數據結構進行解讀,指出關鍵焊縫節點力信息的讀取方法,為分析焊縫疲勞強度提供鋪墊。由于MATLAB的編程環境簡單、完善,編程語言可操作性強,數據處理能力強大,因此通過MATLAB讀取rst文件。
rst文件是Ansys對有限元模型進行強度分析計算后自動生成的擴展名為.rst的數據文件,亦稱結果文件。rst文件包含有限元模型的數據,還包含靜強度分析后各單元節點的應力、應變等信息。本文用MATLAB讀取rst文件,從中提取出單元節點的節點力數據,方便后續采用主-曲線理論計算焊縫的疲勞壽命以及焊縫各節點位置處的應力損傷比。
rst文件采用Fortran語言編寫,屬于二進制文件。要讀取rst文件中包含的節點力數據,首先需要知道整個rst文件的數據結構形式。rst文件數據結構的特殊性決定rst文件的特殊性,須將整個rst文件劃分成多個數據塊讀取。
rst文件包含很多信息,但組成每條完整信息的字節數不一定相同,字節類型也不盡相同。其中,字節類型包括int型、double型、char型、long int型等,long int型字節可以用2個int型字節表示。因此,用MATLAB讀取rst文件時,必須以字節為單位進行讀取,一次可以同時讀取幾個單位字節,但需要確保一次讀取的字節所表達的意思是完整的。
用MATLAB讀取rst文件時,如果命令框中顯示的結果無意義,那么需要改變一次性讀取字節的長度或者類型,有時需要多次改變才能使讀出的數據意義完整。Ansys自帶二進制文件翻譯,可將rst文件的數據塊結構以文本文檔的形式顯示出來,驗證rst文件的讀取進度。圖1為本文模型的rst文件翻譯后部分文檔的截圖,分析圖1可知rst文件的文件頭包含的字節數及其代表的意義等。

圖1 rst文件翻譯后部分文檔截圖
rst文件的數據結構形式并非一成不變,不同的有限元模型使用的單元類型不同、材料屬性不同、加載條件不同,都會導致rst文件的數據結構模塊發生變化。參照Ansys自帶的幫助文檔,本文介紹一般情況下rst文件的數據結構形式。
初始模型為1塊長方形板,在不考慮溫度影響下對其進行有限元網格劃分,單元類型為三維線性4節點殼單元SHELL 181。利用Ansys進行靜態強度分析,形成rst文件。材料屬性為單一線性材料,板厚單一。
每個rst文件包括1個標準的文件頭,文件頭包含100個整形數據的字節,每個字節對應相應的Iterm。這個文件頭包含一些rst文件信息,在利用MATLAB讀取rst文件時可以知道字節所表達的意思是否完整。當然,一些Iterm的內容是字符或者字符串形式,而MATLAB在讀取二進制文件時不能顯示類似數據,讀取比較麻煩,需要切換讀取類型,由整形int型變為char型或其他類型。
文件頭中的信息如下:
Iterm1表示文件數,MATLAB讀取值為12,默認為12;
Iterm2表示文件形式,Iterm2有1和-1共2個值可選,MATLAB讀取值為-1,表示此模型為大數形式;
Iterm3表示時間,MATLAB讀取值為171659,對照翻譯文檔中的時間信息,可知讀取進度正確;
Iterm4表示日期,MATLAB讀取值為1212013,可知讀取字節的意思表達完整;
Iterm5表示單位,MATLAB讀取值為0,代表選擇單位為用戶自定義,符合模型單位;
Iterm6~9在表中無顯示,MATLAB讀取的信息在翻譯文檔中找不到對應信息;
Iterm10表示軟件的版本號,MATLAB不能顯示;
Iterm11為版本發布日期,MATLAB讀取值為20111024;
Iterm12~14、15~16、17~18、19、20~22、23~25分別為機器識別符、工作名稱、Ansys產品名稱、用戶名以及機器識別符,這些都是字符串,MATLAB不能顯示;
Iterm26為系統記錄尺寸,MATLAB讀取值為16384,即此rst文件中包含16384條記錄;
Iterm27為最大文件長度,MATLAB讀取值為312525;
Iterm28為最大記錄號,MATLAB讀取值為20;
Iterm31~38為工作名稱、字符串類型;
Iterm41~60為分析主題、字符串類型;
Iterm61~80為第一子標題、字符串類型;
Iterm95為文件切割點,查看翻譯文檔為未知;
Iterm97~98為長整形文件尺寸,MATLAB讀取值為3125250,占2個整形字節;
Iterm99~100無明顯意義,MATLAB讀取值為0654321;
至此,一個標準的100字節的文件頭表達完整。通過讀取此文件頭,可以知道此文件計算的時間、記錄長度、軟件版本號等相關信息。這些信息是整個rst文件的概括特征,并沒有各節點的節點力等類似數據,要想知道類似信息,還需要繼續往下讀取。讀取方法與上述類似,不再詳述。
在后續的rst文件讀取中,可以知道單元號、節點號、節點坐標、材料屬性、實常數等有限元模型的基本信息,還可以知道各單元節點的應力、應變以及節點自由度等相關信息,所求的節點力信息就包含在其中。讀取完整個rst文件,就可以知道節點力信息具體在哪里。注意:這些節點力信息屬于所有節點,后文會介紹怎樣提取焊縫節點的節點力。
本文的有限元模型基于有限元前處理軟件HyperMesh建立,其有限元網格劃分前處理功能強大,且軟件用戶界面良好,便于操作。HyperMesh支持多種求解器對應文件的輸入、輸出格式,在使用HyperMesh劃分完三維幾何模型的有限元網格后,可以直接把建好的有限元模型轉化為不同求解器對應的文件格式,利用相應的求解器進行后續的分析計算。HyperMesh中求解器為Ansys的cdb文件和Abaqus的inp文件,其中包含焊縫信息。由于讀取inp文件為焊縫信息,不涉及節點力,因此此處不展開。
不同的rst文件具有不同的數據結構,為便于解析,本文取1塊100 mm×100 mm×10 mm的鋼板,其彈性模量=2.1×10MPa,泊松比=0.3,密度=7.85×10g/mm。其幾何模型見圖2。

圖2 幾何模型
劃分有限元模型,網格尺寸為5 mm。有限元模型見圖3,單元數為400個,節點數為441個,采用三維4節點殼單元SHELL 181劃分。

圖3 有限元網格模型
疲勞加載工況為:鋼板在長度方向上受最大拉力為2 100 N的脈動循環載荷,循環次數為1×10次。工況說明:設置最大拉力為2 100 N主要是為加載方便,鋼板沿寬度方向正好有21個節點,每個節點上可加載整100 N的拉力。
脈動循環是一種典型循環方式,在一定的條件下,脈動循環方式可以與對稱循環方式相互轉化。在后續的程序準確性驗證中,可以把每個節點上的載荷調整為50 N,循環方式改為對稱循環,循環次數不變。查看得到的焊縫應力損傷比結果,是否與100 N脈動循環情況下的計算結果一致。建立完整的有限元模型后,將cdb文件導入到Ansys中進行計算,即可得到rst文件,此文件包含程序計算所需的節點力信息數據。
1條完整的焊縫包含1排焊縫單元和1排焊縫節點。焊縫節點的位置即是真實模型中焊縫的焊趾或焊跟部位的位置,而焊縫單元的位置即是真實焊縫中緊鄰焊縫焊趾或焊跟部位的那一部分母材位置。用1排焊縫單元和1排焊縫節點才能表達1條完整焊縫的原因如下:
本文使用SHELL 181單元建立有限元模型,SHELL 181單元是三維4節點單元,每個單元包含4個節點,單元中的節點會被相鄰單元共用。同一個節點對相鄰不同單元的節點力不一定相同,其與單元所處的位置、模型的形狀、載荷和約束有一定關系。僅憑節點編號提取出的節點力并不準確,單用1排節點表示1條焊縫并不能準確提供疲勞計算所需的節點力信息。從單元坐標系的角度講,每個單元都有其對應的單元坐標系,而單元坐標系的建立需要單元的3個節點才能確定,并且單元坐標系可決定焊縫的走向以及計算時的先后順序,因此需要建立1排包含焊縫節點的焊縫單元,并且這些焊縫單元需要在對應焊縫節點的同側,確保計算和單元坐標系的有序性。
在有限元模型中建立焊縫。為使實驗數據更具有參考性,可建立2條焊縫。鋼板的左、右2邊各1條,2條焊縫走向一致,焊縫的節點數和焊縫單元數也相同。焊縫1的焊縫節點數為21,單元數為20;焊縫2的節點數為21,單元數為20。2條焊縫的節點編號和單元編號見表1。建立焊縫時,單元和節點都按順序選擇,而單元編號和節點編號并非按照編號大小依次排列。

表 1 焊縫的單元編號和節點編號
在2條焊縫的建立過程中,焊縫單元和焊縫節點都按順序沿同一方向選取,便于后續計算,更省去有限元模型單元坐標系建立過程中的麻煩。建立的焊縫模型見圖4。

圖4 有限元模型中焊縫
在建立焊縫時,需要注意焊縫單元集名稱和焊縫節點集名稱的寫法。由于在讀取inp文件控件的回調函數中,對焊縫信息的提取是以程序需要識別特定的焊縫名稱為前提的,因此有限元模型中的焊縫名稱需要特定的寫法:焊縫單元需要寫成weld_1_e,焊縫節點需要寫成weld_1_n。名稱中間的數字可以隨著焊縫條數的增加而增大。1排焊縫單元對應1排焊縫節點,焊縫單元與焊縫節點依次建立。建立好焊縫模型后,導出對應求解器為Abaqus的inp文件。
以上完成有限元模型的建立,在Ansys模板下生成對應的cdb文件,在Abaqus模板下生成對應的inp文件。先把cdb文件導入Ansys中進行靜強度計算,得到的結果是模型在疲勞工況作用下所受的應力范圍,或稱應力幅值,其生成的rst文件中包含節點力信息,也稱節點力幅值。得到rst文件后,即可進行文件的讀入和計算。
為掌握rst文件的路徑和文件名,在讀取rst文件時還需要用到uigetfile命令和fopen命令,且以只讀方式打開rst文件。rst文件屬于二進制數據文件,不同的模型對應的rst文件的字節長度不一樣,在讀取rst文件時需按順序存儲每個字節的信息。rst文件模型單元的節點力信息并不是完全按照節點編號從小到大依次排列的,其排列方式與單元的節點順序有一定關系。每個單元包含4個節點,這4個節點按照一定的順序(或順時針或逆時針)排列,每個節點可以有1種或多種節點力,若其相鄰的單元數是4個,則此節點會有4種不同的節點力,見圖5。表2為節點的不同節點力,表中、和表示沿整體坐標系、和這3個方向的節點力,、和表示3個方向的節點力矩。

圖5 與4個單元相鄰的節點示意

表 2 節點的不同節點力
因此,在讀取節點力時需要甄別哪些是焊縫節點的節點力信息。單純按照節點編號的順序讀取焊縫的節點力信息,不能得到想要的節點力信息。焊縫信息包含單元編號和節點編號,其中單元編號尤為重要,其決定所要讀取的節點力信息。
知道節點力信息的排列方式后,就不難提取出焊縫節點對應的節點力信息。為便于后續使用數據方便,按照焊縫信息將節點力信息也放進元胞矩陣中,并且節點力信息的排列順序與相應焊縫的節點編號順序對應。在提取并存放完所有的焊縫節點力數據后,需要對節點力進行轉化。
后續結構應力計算所使用的節點力是單元坐標系(見圖6)下的節點力,而存放好的節點力是整體坐標系下的節點力,需要將整體坐標系下的節點力轉化為單元坐標系下的節點力。使用for循環語句對節點力數據逐個進行轉化,轉化時使用方向余弦向量。單元坐標系需要3個點才能確定,而焊縫單元包含2個焊縫節點,需要再提取出此單元包含的另外一個節點的單元坐標系,同時確定單元節點之間的距離,即單元寬度。這些功能可以在讀取節點力文件時完成,從而得到單元坐標系下的節點力。

圖6 單元坐標系示意
整個讀取過程可按照以下5個步驟進行:
(1)通過HyperMesh建模及Ansys有限元計算,得到rst文件。
(2)建模時規定命名規則,即在Abaqus模板下生成對應的inp文件。采用HyperMesh在Abaqus模板下生成對應的inp文件,其中記錄單元編號及節點編號。
(3)使用MATLAB讀取文件,使用uigetfile命令和fopen命令打開文件。
(4)根據文件頭釋義,結合數據類型和數據長度讀取文件,并將整體坐標系轉化為單元坐標系。
(5)根據焊縫節點及焊縫單元的排列規則,進一步讀取相應的節點力信息。
將上述方法得到的節點力信息轉化為等效結構應力,并選擇合適的-曲線,即可計算出焊縫的應力損傷比。將該結果與疲勞仿真軟件FE-WELD的計算結果進行對比,見表3。用本文所述方法得到的節點力有一定的準確性,并且2種計算結果具有一定的相似性,見圖7。

表 3 焊縫1的應力損傷比結果對比

圖7 2種計算結果對比
本文采用三維4節點殼單元建立典型方形鋼板的有限元模型,加載并計算后得到包含單元節點力信息的rst文件。重點介紹rst文件的數據結構及節點力信息的排序規則,以及利用MATLAB讀取節點力信息的方法。將讀取的節點力轉化后得到的結果與專業疲勞仿真軟件的計算結果進行對比,驗證讀取信息的準確性以及讀取方法的有效性,為焊縫節點力的獲得提供一種方法。