三維塊單元慣性釋放方法研究

王宇帆， 劉 虓， 陳超核

（華南理工大學土木與交通學院，廣州 510640）

1 前言

航行中的船舶、飛行中的飛機、行駛中的汽車均屬于自由無約束結構，采用有限元法對此類結構進行靜力分析時，仍需要施加約束來消除剛度矩陣的奇異。當此類結構受到的外載荷滿足靜力平衡條件時，求解后約束點處的支反力為零；由于在載荷的計算和施加過程中存在諸多近似，很難得到絕對平衡力系，約束點處產生的支反力會改變結構的實際受力狀態，引起局部應力集中、變形失真、傳力路徑錯誤等異常現象。為解決這一問題，研究人員將慣性釋放技術引入此類結構的有限元強度計算中，這是一種基于達朗貝爾原理的慣性平衡技術，它首先計算結構體在不平衡外載荷下產生的加速度，再根據加速度計算出慣性載荷，最后將慣性載荷與外載荷疊加從而使結構體達到靜態平衡。

商用有限元軟件MSC、ANSYS、ABAQUS 中均具備該功能，被廣泛應用于船舶、航空、汽車領域。在船舶領域，張少雄[1]等使用慣性釋放功能對某油船艙段結構進行強度計算，消除了邊界條件對應力計算結果的影響；趙振宇[2]等應用慣性釋放功能對某浮船塢結構進行強度計算，得到了精度較高的應力響應；在航空領域，路林華[3]等使用慣性釋放功能對某無人直升機進行了整機強度和剛度分析；Boni[4]等使用慣性釋放方法研究了處于自由飛行中的太陽帆在不同輻射角下的結構行為；在汽車領域，孫輝[5]等使用慣性釋放功能對某車架進行了強度分析，得到了較符合實際的車架應力分布；Vdovin[6]等對某車架分別使用慣性釋放法和傳統約束法進行強度分析，結果表明慣性釋放法得到的結果更加真實；Banshan[7]等針對某汽車平衡軸架使用慣性釋放和子模型法得到了結構的應力分布，并通過與實際結構的疲勞部位對比證實了該方法的有效性。

以上研究借助國外軟件的慣性釋放功能，取得了理想的效果。但國內極少有人研究慣性釋放算法原理并開發對應的程序。慣性釋放算法主要有兩個關鍵點：（1）計算結構體在外載荷作用下產生的加速度aI；（2）計算aI引起的節點慣性載荷fI。ANSYS 幫助文件[8]介紹了aI的計算方法，但可能出于保護商業機密的考慮，對fI并沒有具體介紹；國產自主有限元系統SiPESC.FEM 具備慣性釋放功能，但開發者路林華[9]也沒有給出節點慣性載荷fI的具體算法。

本文根據能量原理，獨立推導三維塊單元由平移和旋轉加速度引起的等效節點慣性力，以該算法為基礎編寫了慣性釋放程序，并通過兩個算例驗證了算法的精度、可靠性和工程適用性。程序和算例文件，可從http://www.huagongchuanhai.cn/fem/獲取。

2 慣性釋放原理及實現

假設結構有限元模型由p個塊單元、q個節點組成，其慣性釋放算法在有限元中的實現步驟如下：

（1）計算結構體在不平衡外載荷下產生的加速度

由文獻[8]可知，根據達朗貝爾原理，外載荷和慣性載荷組成平衡方程：

式中：Ft和Fr為質心處所受合外力及合外力矩；

Mt和Mr為有限元模型的平動質量矩陣及轉動質量矩陣；

at和ar為結構體在外載荷作用下產生的平移加速度及繞質心的旋轉加速度。

根據公式（1），計算出平移加速度at和旋轉加速度ar。

（2）根據加速度計算節點慣性載荷

這是慣性釋放方法的一大難點，目前國內尚沒有學者提供計算單元節點慣性力的具體算法，在國外公開的文獻中也查不到。本文根據能量方法推導出等效節點慣性力，分別考慮了平移加速度和旋轉加速度引起的等效節點慣性力。

平移加速度等效節點慣性力：

將節點慣性力施加在外載荷向量上，得到新的平衡外載荷向量:

在以上算法的基礎上，編寫了具有慣性釋放功能的三維實體八節點等參單元計算程序，下面通過兩個算例來驗證其可靠性和工程適用性。

3 算例一

如圖2 a）所示：一條長、寬、高分別為L=2 000 mm、B=100 mm、H=120 mm 的 梁，密 度ρ=7.85 e-5 N/mm3，受均布載荷q=30 N/mm，由兩端的反力ql維持平衡，求跨中上緣x 方向應力；如圖2 b）所示，同樣尺寸和材質的梁，去掉跨中載荷，僅保留兩端集中力F=ql，在考慮慣性釋放的情況下，該問題與圖2 a）等效。

圖2 實體梁

3.1 理論解

由彈性力學[10]可知，該梁跨中x 方向應力為：

3.2 有限元解

有限元網格，如圖3 所示：對于三維塊單元模型，需要對6 個自由度方向的位移進行限制以防止剛體運動。邊界條件如下：約束A 點的UX、UY、UZ 自由度，B 點的UY、UZ 自由度；C 點的UY 自由度。需要注意的是，計算完畢后要檢查這6 個自由度方向對應的支反力是否為零。如果不為零，則說明邊界條件錯誤，因為它干擾了結構的變形。

圖3 實體梁有限元模型及邊界條件

設置材料密度后，分別使用本文程序和ANSYS（使用SOLID185 單元）進行了慣性釋放計算，計算結果如圖4所示。本文程序的積分方式采用6×6×6高斯積分，ANSYS 的單元技術采用默認的全積分法。

圖4 實體梁x 方向應力云圖

3.3 結果對比

值得注意的是，在ANSYS 求解中，通過設置不同的KY（2）值將采用不同的單元技術，計算結果也會不同。對比計算結果，如表1 所示。

表1 實體梁算例結果對比

4 算例二

本算例以某鋼管（取自某導管架平臺組裝件）為對象，分別采用本文程序和ANSYS 對鋼管吊裝過程進行慣性釋放計算，并對比計算結果。

4.1 有限元模型建立

鋼管結構參數見表2，有限元模型見圖5。通過約束A 點的UX、UY、UZ 自由度，B 點的UY、UZ 自由度，C 點的UY 自由度，限制鋼管的剛體運動。

圖5 導管架平臺組裝件——鋼管

表2 鋼管參數

4.2 載荷施加

將起吊力以集中力的形式施加在吊繩與鋼管接觸的半圓柱面節點上，如圖6 a）所示：起吊時采用兩條吊繩，分別安裝在距離鋼管兩端500 mm 處；吊繩選取直徑為80 mm 的鋼絲繩，起吊后鋼絲繩與鋼管的接觸寬度為鋼絲繩直徑的一半，即接觸寬度w=40 mm。

圖6 吊裝力載荷模擬

如圖6 b）所示，吊繩與半徑為r 的圓柱接觸處單位長度所受法向支持力[11]：

假設沿圓周每隔dθ布置一個節點，沿接觸寬度方向布置m個節點，則接觸區域每個節點受指向圓心的法向力：

將法向力分解為豎直向上和水平向內的集中力：

將F1和F2施加到接觸區域的節點上，如圖7 所示。

圖7 鋼管起吊載荷

4.3 結果評估

約束點的反力為零，可以作為慣性釋放解有效性的判別準則[12]。程序計算完成后，首先檢查程序輸出的節點支反力文件，獲得約束節點的支反力均為零，表明邊界條件對結構變形并沒有影響。

使用可視化后處理軟件Tecplot 打開程序計算輸出的節點應力和位移文本，可獲得鋼管的VonMises 應力云圖及變形圖，如圖8 所示：采用慣性釋放計算后得到的結構變形是相對于約束節點的，可以看出主要變形部位位于接觸區域，應力也集中在吊繩與鋼管接觸區域，尤其在接觸區域底部的鋼管內側應力達到最大；在約束點處沒有出現類似支座附近的應力集中、變形失真等異常情況，比較合理地反映了鋼管在起吊過程中的應力分布規律。

圖8 鋼管VonMises 應力云圖及變形圖（程序解）

采用相同的數據模型，使用ANSYS 進行了同樣的計算（單元技術采用默認的全積分法），結果如圖9所示。

對比圖8 和圖9 可以發現：程序與ANSYS 計算結果的應力和變形分布一致，二者的最大應力略有差別，程序最大應力為0.866 0 Mpa，ANSYS最大應力為0.937 9 Mpa，相對偏差為7.67%。

圖9 鋼管VonMises 應力云圖及變形圖（ANSYS 解）

5 結論

本文在國內首次提出了三維塊單元節點慣性力的計算方法，在此基礎上編制了相應程序，并通過兩個算例進行對比驗證，結論如下：

（1）理論算例證明本文提出的三維塊單元節點慣性載荷計算方法是可靠的；

（2）約束節點的支反力為零，驗證了本文提出的慣性力計算方法能夠消除邊界節點對變形的干擾；

（3）本文算法獲得了符合實際的應力分布，合理地反映了鋼管在起吊過程中的應力分布規律：主要變形和應力集中在吊繩與鋼管接觸區域；約束節點處沒有出現應力集中、變形失真等異常情況；計算結果與ANSYS 相當，驗證了本文算法的工程適用性。

本文從算法設計到程序實現乃至工程應用，既參考了國外軟件和國內外同行的成果，也獨立研究了核心算法。文中提出的慣性力計算方法除了三維塊單元，還可擴展到平面應力單元、梁單元、薄板彎曲單元等，對國產有限元的自主開發具有一定的借鑒意義。