王宇帆， 劉 虓， 陳超核

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

1 前言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

在以上算法的基礎(chǔ)上，編寫了具有慣性釋放功能的三維實(shí)體八節(jié)點(diǎn)等參單元計(jì)算程序，下面通過(guò)兩個(gè)算例來(lái)驗(yàn)證其可靠性和工程適用性。

3 算例一

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

圖2 實(shí)體梁

3.1 理論解

由彈性力學(xué)[10]可知，該梁跨中x 方向應(yīng)力為：

3.2 有限元解

有限元網(wǎng)格，如圖3 所示：對(duì)于三維塊單元模型，需要對(duì)6 個(gè)自由度方向的位移進(jìn)行限制以防止剛體運(yùn)動(dòng)。邊界條件如下：約束A 點(diǎn)的UX、UY、UZ 自由度，B 點(diǎn)的UY、UZ 自由度；C 點(diǎn)的UY 自由度。需要注意的是，計(jì)算完畢后要檢查這6 個(gè)自由度方向?qū)?yīng)的支反力是否為零。如果不為零，則說(shuō)明邊界條件錯(cuò)誤，因?yàn)樗蓴_了結(jié)構(gòu)的變形。

圖3 實(shí)體梁有限元模型及邊界條件

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

圖4 實(shí)體梁x 方向應(yīng)力云圖

3.3 結(jié)果對(duì)比

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

表1 實(shí)體梁算例結(jié)果對(duì)比

4 算例二

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

4.1 有限元模型建立

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

圖5 導(dǎo)管架平臺(tái)組裝件——鋼管

表2 鋼管參數(shù)

4.2 載荷施加

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

圖6 吊裝力載荷模擬

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

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

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

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

圖7 鋼管起吊載荷

4.3 結(jié)果評(píng)估

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

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

圖8 鋼管VonMises 應(yīng)力云圖及變形圖（程序解）

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

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

圖9 鋼管VonMises 應(yīng)力云圖及變形圖（ANSYS 解）

5 結(jié)論

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

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

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

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

本文從算法設(shè)計(jì)到程序?qū)崿F(xiàn)乃至工程應(yīng)用，既參考了國(guó)外軟件和國(guó)內(nèi)外同行的成果，也獨(dú)立研究了核心算法。文中提出的慣性力計(jì)算方法除了三維塊單元，還可擴(kuò)展到平面應(yīng)力單元、梁?jiǎn)卧⒈“鍙澢鷨卧?，?duì)國(guó)產(chǎn)有限元的自主開發(fā)具有一定的借鑒意義。