王宇帆， 劉 虓， 陳超核

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

1 前言

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

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

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

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

2 慣性釋放原理及實現(xiàn)

假設結(jié)構(gòu)有限元模型由p個塊單元、q個節(jié)點組成，其慣性釋放算法在有限元中的實現(xiàn)步驟如下：

（1）計算結(jié)構(gòu)體在不平衡外載荷下產(chǎn)生的加速度

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

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

Mt和Mr為有限元模型的平動質(zhì)量矩陣及轉(zhuǎn)動質(zhì)量矩陣；

at和ar為結(jié)構(gòu)體在外載荷作用下產(chǎn)生的平移加速度及繞質(zhì)心的旋轉(zhuǎn)加速度。

根據(jù)公式（1），計算出平移加速度at和旋轉(zhuǎn)加速度ar。

（2）根據(jù)加速度計算節(jié)點慣性載荷

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

平移加速度等效節(jié)點慣性力：

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

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

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）所示，同樣尺寸和材質(zhì)的梁，去掉跨中載荷，僅保留兩端集中力F=ql，在考慮慣性釋放的情況下，該問題與圖2 a）等效。

圖2 實體梁

3.1 理論解

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

3.2 有限元解

有限元網(wǎng)格，如圖3 所示：對于三維塊單元模型，需要對6 個自由度方向的位移進行限制以防止剛體運動。邊界條件如下：約束A 點的UX、UY、UZ 自由度，B 點的UY、UZ 自由度；C 點的UY 自由度。需要注意的是，計算完畢后要檢查這6 個自由度方向?qū)闹Х戳κ欠駷榱?。如果不為零，則說明邊界條件錯誤，因為它干擾了結(jié)構(gòu)的變形。

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

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

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

3.3 結(jié)果對比

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

表1 實體梁算例結(jié)果對比

4 算例二

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

4.1 有限元模型建立

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

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

表2 鋼管參數(shù)

4.2 載荷施加

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

圖6 吊裝力載荷模擬

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

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

將法向力分解為豎直向上和水平向內(nèi)的集中力：

將F1和F2施加到接觸區(qū)域的節(jié)點上，如圖7 所示。

圖7 鋼管起吊載荷

4.3 結(jié)果評估

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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