伸縮臂結構超級單元建模及失穩荷載快速搜索方法

卓 英 鵬， 齊 朝 暉*， 王 剛， 徐 金 帥， 趙 天 驕

( 1.大連理工大學 工業裝備結構分析國家重點實驗室， 遼寧 大連 116024；2.大連理工大學 海洋科學與技術學院， 遼寧 盤錦 124221 )

0 引 言

流動式起重機因其靈活機動性強、伸展范圍大，廣泛應用在工程機械各領域，其中伸縮臂結構作為其主要的承載部件，起著不可替代的作用[1-3]．起重機實際工作時的荷載高度可達到幾十米甚至上百米，伸縮臂展開后的結構通常都是細長形的，具有明顯的幾何非線性效應，此外，結構的穩定性也是限制起重性能的主要因素．

伸縮臂結構分析方面已經積累了大量的相關研究[4-7]，這些研究在一定程度上解決了很多理論和工程實際問題，其中相當部分的工作是直接采用商業軟件中殼單元完成整個結構的有限元建模，而龐大的自由度和復雜的接觸約束往往嚴重影響計算效率和收斂性．現有的很多文獻、規范[8-13]中將伸縮臂視為變截面階梯柱，采用能量法、瑞利-利茲法、精確有限元法等求解其變形、失穩等情況，實質上仍停留在截面變化的單個構件問題上，由多個截面厚度變化的空心單體組成的多結構復雜伸縮臂問題鮮有研究．而且，這種等效階梯梁模型采用的是工程上常見的重疊部分抗彎抗扭模量疊加的方法，無疑會在重疊部分產生虛假的約束反力．

工程中承載后的細長結構整體表現為大變形，但局部區域內仍處于小變形范疇，結合這一特點，很多專家學者提出相應的幾何非線性計算方法[14-15]，共旋坐標法[16]相對于早期參考初始構型的TL列式、參考當前構型的UL列式[17-19]等呈現出簡單、高效的優勢：可以應用現有的線性單元，省去復雜非線性單元的構造，深受工程設計人員的青睞[20-22]．

傳統起重機臂架結構穩定性分析常采用的特征值屈曲分析方法適用于發生屈曲時仍處于線性變形狀態的工況．對于重載條件下大變形失穩問題，追蹤結構位移-荷載曲線平衡路徑成為近代非線性穩定性理論發展的主要方向[23]，包括增量法、弧長法及改進的弧長法等各種方法[24-27]應運而生．這些方法在工程實踐中被廣泛應用，針對不同的具體問題，體現出不同的優缺點，但是均面臨在保證精度前提下如何調整步長的困惑．

伸縮臂結構具有組合規則、計算工況多的特點，往往需要大批量性能計算，因此，很值得研究的問題是：能否在準確描述變形幾何非線性效應及臂節間嵌套約束關系的前提下，盡可能減少自由度規模，快速地搜索結構的失穩荷載．鑒于此，本文提出一種可以滿足上述要求的伸縮臂結構的超級單元建模及快速搜索失穩荷載的方法：采用梁單元對各臂節進行離散，以相鄰臂節間接觸截面節點、油缸鉸接面節點作為邊界點建立子結構；子結構內部節點自由度可通過內部平衡方程由邊界節點自由度表示，進而轉換為形式上兩節點超級單元；以超級單元節點自由度作為系統變量，建立含荷載參數的結構非線性平衡方程及相應的切線剛度陣；最后采用平衡方程的微分形式快速求解結構的變形平衡路徑及失穩荷載．

1 變截面臂節的梁單元拼接方法

組成伸縮臂結構的臂節長度尺寸遠大于截面尺寸，是典型的空心細長結構，因此，臂節可采用常見的梁單元建模．梁單元的變形虛功率

(1)

其中，單元節點參數可表示為

(2)

(3)

圖1 同一臂節中臂厚不同的截面Fig.1 The change of cross-section in the same boom segment with different thickness

(4)

其中

(5)

節點不在形心處梁單元變形虛功率和重力虛功率

(6)

其中

(7)

(8)

通過選擇合適的截面特征點作為節點，建立截面形心和節點間的位移協調條件，實現變截面臂節各梁單元的拼接．

2 子結構單元位移分解及自由度凝聚

針對細長結構的幾何非線性問題，最早由Wempner[28]和Belytschko等[29]分別提出的共旋坐標法將單元的位移場分解為隨單元坐標系的剛體轉動以及相對于單元的小位移．這類單元在具體應用中也體現出明顯的優勢[15]．但另一方面它與傳統方法是一樣的，如果要做非線性分析，結構中所有單元都要按非線性單元處理，最終導致結構平衡方程是一組維數很高的非線性方程．事實上，根據伸縮臂結構的特點，每個臂節可以合理地劃分若干子結構，幾何非線性效應主要體現在子結構連體坐標系大位移大轉動上，子結構內部各節點相對于連體坐標系的位移轉動可按小位移小轉動處理．

圖2所示的子結構由n個梁單元組成，將梁截面坐標系原點作為廣義梁單元的節點，其相對于總體坐標系的矢徑分別為r0，r1，…，rn，總體轉動參數(卡爾丹角)為θ0，θ1，…，θn．R0，R1，…，Rn分別為截面瞬時方位矩陣，Ri(i=0，1，…，n)的列向量可由總體轉動參數θi中各元素αi、βi、γi表示為

圖2 子結構內梁單元Fig.2 Beam elements in the substructures

(9)

其中

(10)

定義子結構最左側截面坐標系{t0，b0，s0}為子結構連體坐標系，子結構內任意一點的矢徑

(11)

及其變化率

(12)

(13)

(14)

子結構內部任意一點的總體自由度虛速度分解為子結構連體坐標系的平動轉動自由度和當前連體坐標系的局部自由度虛速度．由于連體坐標系中的局部自由度和總體坐標系下的總體自由度不一致，為將虛功率方程轉化為代數方程，需要給出總體自由度和局部自由度之間的函數關系．梁單元節點在子結構連體坐標系下的位移可表示為

(15)

它的變化率

(16)

(17)

它的變化率

(18)

其中ωi是截面坐標系轉動的角速度，它與總體轉動參數θi中各元素αi、βi、γi轉換關系為

(19)

最左側截面節點位移轉角相對于子結構連體坐標系為零，在子結構連體坐標系內描述的所有節點位移轉角均為小位移小轉角，滿足線性梁單元的特征，且每個子結構中的節點可分為兩類：(1)與其他子結構連接的邊界節點集nb；(2)內部節點集ni．因此，可將子結構內部自由度凝聚到邊界自由度上．

子結構坐標系下，對組成子結構的每個梁單元進行剛度陣組裝并按內部節點自由度和邊界節點自由度進行分塊：

(20)

對應的變形虛功率

(21)

(22)

體力虛功率

(23)

式中：fg為子結構內任意一點的重力線密度；F0、T0是體力等效節點力相對于連體坐標系原點的合力和合力矩，表示為

(24)

(25)

從中可得

(26)

臂節被劃分為多個子結構后，臂節大位移大轉動可以描述為子結構坐標系的大位移大轉動，子結構內部的梁單元相對于子結構坐標系的變形屬于小變形，因此可以認為單元節點在子結構坐標系下的位移和轉動為小量．在此條件下，子結構剛度陣中元素為常量，因而

(27)

將式(26)、(27)代入式(21)中可得

(28)

體力虛功率

(29)

其中關聯矩陣

(30)

等效剛度陣、等效重力影響系數分別表示為

(31)

這樣，由n個梁單元拼接組裝的子結構就縮減成了用兩端節點位移和轉角表示所有節點自由度的超級廣義梁單元．

3 伸縮臂子結構劃分及并聯約束處理

常見多物體結構由受約束的典型的局部剛性區域依次串聯而成，可稱之為串聯連接；伸縮臂結構是由多個空心結構通過層層環繞嵌套而成的系統，伸縮臂結構中臂節與臂節的連接方式區別于傳統多物體結構，工作狀態時臂節沿軸向的滑移運動被約束，整個臂節均為可發生大變形的柔體，這種連接可稱為并聯式連接．如圖3所示，伸縮臂各臂節間通過銷軸銷孔連接，它們雖然存在縫隙，但與臂節的截面相比十分微小，且各臂節抗彎模量在同一數量級，接觸區域都局限在內層臂左端和外層臂右端很窄的區域內；每個臂節最多有4個截面受其他臂節的作用：連接銷軸的左端面、與其外層臂右端面接觸的截面、與其內層臂銷接的銷孔所在截面、與內層梁接觸的右端面，這4個特殊截面節點可分別稱之為左端節點、外接觸點、內接觸點、右端節點．根據伸縮臂的這一特點，以這4個節點作為邊界點對每個臂節劃分子結構，每一個子結構內部可繼續劃分多個梁單元．這樣，具有m個臂節的復雜伸縮臂結構可由共含有6(4m-2)個節點自由度的3m-2個子結構進行幾何非線性分析．

圖3 伸縮臂系統子結構劃分Fig.3 Substructure division of the telescopic boom system

圖4 伸縮臂結構臂節間的約束關系Fig.4 Constraints between boom segments of telescopic boom structure

(32)

第二類約束關系中截面②、④間的滑塊限制二者沿截面內主軸方向的相對平移以及繞面外法向的相對轉動，具有多體系統理論中棱柱鉸的特點，對應的約束方程可寫為

(33)

依據式(33)，選擇r2、α2、β2、γ2、r4為獨立變量，α4、β4、γ4為非獨立變量．由此，式(32)、(33)共同組成伸縮臂結構兩兩相鄰內外層臂間約束方程．

4 附加結構約束關系及等效節點力

傳統結構問題的邊界條件通常是直接約束位移，對這類約束的處理方法已經十分成熟．對于流動式起重機的伸縮臂結構，最外層臂(基本臂)與轉臺間為銷軸鉸接，存在一個機構自由度；該自由度則由臂節變形和工作狀態下固定長度的變幅油缸支撐共同決定．此外，臂節所受到的外力不僅包括自重以及荷載，還包括臂頭、起升繩施加的附加外力．這些附加結構向臂節施加的外力不僅與結構上的外力作用點位移相關，而且和約束結構本身的運動有關．

4.1 基本臂邊界約束條件

轉臺與基本臂鉸接部分如圖5所示．由于銷軸位置受到很大集中力，該處會進行加強處理，與臂節左端面焊接．

圖5 基本臂與轉臺連接部分Fig.5 Connection between basic boom and turntable

左端截面A0的節點矢徑r1與轉動參數α1、β1、γ1存在約束關系：

(34)

其中r0為銷軸中心點的矢徑，選擇初始狀態轉軸方向與銷軸軸向重合的角度β1作為獨立變量，剩余變量r1、α1、γ1為非獨立變量，則R0可表示為

(35)

變幅油缸下鉸點與轉臺鉸接，上鉸點與基本臂鉸接．工作狀態下，假設油缸長度保持恒定，本質上對上下鉸點間形成固定距離約束．上鉸點可由其所接觸的截面A2的節點矢徑r2和轉動參數α2、β2、γ2表示為

(36)

其中R2為截面坐標系基矢量組成的矩陣．變幅油缸對基本臂的約束關系為

(37)

式中：d為油缸的長度，選擇節點r2第一個分量作為獨立變量，截面A2的其余描述參數為非獨立變量．

4.2 臂頭、荷載、起升繩附加節點力

流動式起重機伸縮臂臂頭焊接安裝在最內層臂的右端面上，如圖6所示．起升繩一端連接于轉臺的卷揚上，另一端相切于導向輪表面．

圖6 臂頭與起升繩連接部分Fig.6 Connection between boom head and rope

臂頭所連的最內層臂截面節點用D表示，將臂頭視為質量為m的剛體，臂頭的質心、起升繩作用點和荷載作用點在與臂頭相連的截面坐標系Rn中的矢徑分別為rm、rg、rp．則臂頭質量和荷載對節點D處附加等效節點力可分別表示為

Tm=(mg(Rnrm)×mg)

(38)

(39)

起升繩繞過導向輪后，與n倍率起升滑輪組相連，繩索拉力大小為荷載G的1/n；工作狀態下的方向與節點位移有關，可表示為

(40)

轉換為節點D處附加等效節點力

Tn=(Gn/nRnrp×Gn/n)

(41)

5 結構平衡方程及切線剛度陣

以各臂節凝聚后的子結構邊界節點自由度作為系統變量，將其按平動和轉動劃分，伸縮臂結構整體虛功率方程可以表達為

(42)

(43)

(44)

(45)

(46)

其中

(47)

(48)

子結構k的等效節點力和力矩對應的虛功率方程可改寫為

(49)

定義系統描述變量q，結合式(42)和(49)，伸縮臂結構的整體虛功率方程可表示為

(50)

不失一般性，結構的約束方程可表示為

g(q)=0

(51)

根據相鄰臂節間約束關系和邊界條件，系統變量劃分為非獨立和獨立兩部分：

q=(qdpqip)

(52)

式(51)所派生的變分約束方程可表示為

(53)

非獨立部分虛變分可由獨立部分表示，存在關系

(54)

因而，由虛功率方程(50)可得用于求解的代數方程

(55)

此時，方程的個數與獨立變量個數相同，需要添加與非獨立變量相同個數的獨立方程才可滿足求解完備性，將平衡方程(55)與約束方程(51)相結合，可得到求解節點整體參數所需的方程：

(56)

這是一組高度非線性的方程，給出相應的切線剛度陣可大幅度提高求解效率．

為獲得平衡方程的切線剛度陣，需要給出方程(55)的微分形式：

(57)

式子第1項中

Jfq=?f/?q

(58)

(59)

(60)

(61)

節點局部參數與總體參數的微分之間的關系可由式(45)、(46)得到，即

(62)

(63)

式子第2項需要將廣義力f代入整體分析：

(64)

其中fdp為廣義力陣f中獨立自由度對應的子矩陣．補充的約束方程對應的微分形式可由式(53)得到，即

dg=Gdpdqdp+Gipdqip?Gdq

(65)

將式(58)～(64)代入式(57)，結合式(65)，可得所需的切線剛度陣．

6 結構非線性平衡方程的微分形式

伸縮臂結構的變形求解方程是一組高度非線性方程，初值的選取嚴重影響迭代結果的收斂性．荷載是造成節點位移變化的主要參量，通過跟蹤伸縮臂結構在不同荷載下的平衡路徑可以快速求得任意荷載對應的變形以及最終的失穩荷載．整體伸縮臂結構的平衡方程可寫作以下形式：

R(q)+λG0=0

(66)

式中：q是結構的節點位移和轉角組成的列陣；荷載G按照某種比例施加在結構上，可表示為一個單位荷載增量G0和一個荷載控制參數λ的形式，即G=λG0．

傳統荷載增量法將施加荷載分多個荷載步逐級增加，以前一步的計算結果作為當前步的方程初值，大大改善收斂程度，進而來追蹤平衡路徑．但是這種方法需要人為設置一個固定的荷載增量，該固定增量設置太小增加計算量，降低計算效率，太大又可能收斂失敗，而且由于事先無法準確判斷失穩荷載所在區間，太大的荷載增量很有可能直接越過極值點，導致搜索失穩荷載失敗．

鑒于此，本文考慮將伸縮臂結構平衡方程對荷載控制參數求導，轉化為微分方程，借助常規微分方程求解器中自動調節步長的優勢，實現根據系統當前荷載狀態對應的非線性程度自動調節荷載步大小的功能，在保證每一步求解收斂的前提下，快速追蹤平衡路徑和搜索失穩荷載．其主要操作如下：

(1)結構平衡方程(66)對荷載控制參數λ求導，得到微分方程

(67)

式中：(?R/?q)為平衡方程的切線剛度陣，?q/?λ反映當前荷載下位移-荷載曲線的非線性程度．

(2)依據平衡方程(66)，計算空載(λ=0)對應的位移轉角q0，作為微分方程求解的初值．

(3)給定荷載控制系數上限λmax，在區間λ∈(0，λmax)積分微分方程(67)，積分過程中：(i)每一步用平衡方程(66)修正積分變量；(ii)時刻檢查終止條件(5)．

(4)觸發終止條件(5)，計算終止，得到失穩荷載以及加載至失穩荷載的平衡路徑；沒有觸發終止條件(5)，計算至λmax，得到荷載上限Gmax=λmaxG0以及加載至荷載上限的平衡路徑．

(5)如果系統遇到極值點失穩，求解微分方程的方法無法越過極值點，數值表現為位移對荷載系數的導數趨于無窮大，因此可設置

(68)

來判斷結構是否失穩，其中ε為給定的閾值，代表位移-荷載控制參數值變化的倍數，可根據具體情況設定，當微分求解過程中位移對荷載控制參數的導數值變化倍數達到給定閾值時，此時對應的荷載就是結構的失穩荷載．

上述方法借助常規微分求解器自動調節步長的優勢，可以在初始荷載增加階段(線性變形階段)以較大步長快速求解；同時在接近失穩階段，迅速自動調節以較小步長逐步逼近失穩荷載．伸縮臂的綜合性能由失穩荷載和強度破壞荷載共同決定，但失穩荷載的搜索為破壞荷載計算的前提，前者為后者提供了大量數據，因此在此基礎上也可進一步完成對破壞荷載的計算搜索．

7 數值算例

圖7所示為由彈性模量2.1×1011Pa、泊松比0.3、密度7 850 kg/m3的7節U形截面臂節相互嵌套組成的流動式起重機伸縮臂結構；臂節長度為L，每個臂節含1個銷軸、4個銷孔(從左至右，銷孔號依次為1、2、3、4)，它們沿臂節長度位置分別記為L0、L1、L2、L3、L4；下蓋板在長度位置為Ls處厚度發生變化，變化量為Δt；截面形式采用本文提到的由規則形狀的圓弧組成的非規則復雜截面；各臂節具體相關尺寸見表1、2．基本臂左端面固連的轉軸銷接于回轉平臺，通過變幅油缸的支撐實現變幅．

圖7 流動式起重機伸縮臂結構Fig.7 The telescopic boom structure for mobile crane

表1 臂節長度、銷軸銷孔位置、下蓋板厚度變化位置和變化量(單位：mm)Tab.1 Boom segment length, position of pin shaft and hole, position and change of thickness of lower cover plate (units: mm)

表2 臂節的截面尺寸(單位：mm)Tab.2 Cross-sectional size of boom segment (units: mm)

7.1 臂節子結構劃分

以臂節間滑塊相互接觸截面節點、變幅油缸鉸接截面節點作為邊界點，沿長度方向將各臂節劃分為如圖8所示的子結構，其中3號節點為油缸連接截面節點，最內層臂劃分為2個子結構，其余臂均為3個子結構，7節臂共需27個節點、20個子結構超級單元即可劃分完畢．

圖8 伸縮臂子結構單元劃分Fig.8 Element division of the telescopic boom substructures

7.2 不同工況下伸縮臂結構的變形

針對圖7所示的伸縮臂結構，應用本文方法得到了不同工況下的變形結果．具體的工況選擇如表3所示，其中θ為初始狀態下變幅角，k1、k2、k3、k4、k5、k6分別為由外層至內層1～6節臂相鄰內層臂左端銷軸依次插入的銷孔號．

表3 計算工況Tab.3 Calculated working condition

(a) 工況1

由圖9可看出：相同工況下，隨著荷載的增大，伸縮臂的幾何非線性效應逐漸增大；變幅角增大后，結構的承壓效果明顯，相同荷載作用時，彎曲變形會明顯減小．

方便對比，采用ANSYS商業軟件的Shell181單元建立工況2、4、5對應的有限元模型，變幅油缸使用Link180，相鄰臂節間滑塊接觸部分的自由度局部耦合．

表4、5分別列出了工況2、4、5在不同荷載作用時，利用本文方法計算得到的臂頭導向輪軸中心在豎直方向的位移和計算時間，以及與ANSYS 殼單元建模分析結果的對比．由表中數據可見，本文方法得到的位移解與ANSYS解差別很小．同時，采用ANSYS殼單元建模時劃分82 806個節點，82 279個單元，需要求解496 836個非線性方程；本文方法劃分了27個節點，20個超級單元，需求解162個非線性方程，求解過程變步長的實現會節約線性變形階段的求解時間．此外，只需表1、2所示的尺寸參數即可實現程序化建模，無須重復ANSYS中每次建模的過程，為流動式起重機起重性能大批量快速計算提供了可靠的模型．

表4 工況2、4、5不同荷載作用下的豎直位移Tab.4 The vertical displacements under different loads of conditions 2, 4, 5

表5 工況2、4、5不同荷載作用下的計算時間Tab.5 The calculation time under different loads of conditions 2, 4, 5

7.3 全伸狀態下伸縮臂結構的失穩荷載

圖10給出了參考點位移相對于荷載的變化曲線，其中u表示第7節臂最右端節點位移，v表示第4節臂最右端節點位移．從曲線可看出，當荷載較小時，豎直位移基本保持線性；當荷載接近失穩時，位移曲線明顯變陡，且表現出明顯的非線性；事實上追蹤平衡路徑時，當荷載增加到接近極值點時，步長增量會逐漸趨于無窮小，當荷載位移曲線的斜率達到某個值時，已經非常接近失穩荷載，繼續往下計算會浪費大量的計算時間，實際工程中可綜合計算效率和精度要求，調整ε．

圖10 伸縮臂參考點位移Fig.10 The displacements of reference points of telescopic boom

圖11顯示了搜索失穩荷載平衡路徑過程中荷載為209.33 kN時結構的變形圖，為顯示清晰，位移放大5倍．可以看出，此時伸縮臂結構的變形十分明顯，具有非線性大位移的特征．

圖11 失穩荷載下伸縮臂的變形Fig.11 Deformation of telescopic boom under critical load

8 結 論

(1)考慮腹板厚度變化，選擇合適的邊界節點建立子結構，通過內部自由度凝聚，獲得可用于描述幾何非線性效應的兩節點超級單元．

(2)給出了相鄰臂節間并聯式約束關系方程、邊界條件及附加節點力，推導了含荷載控制參數的結構非線性平衡方程及相應切線剛度陣．

(3)將平衡方程轉化為其微分形式，結合現有微分方程的數值方法，得到一種快速計算伸縮臂結構變形路徑及失穩荷載的方式．通過數值算例，驗證了方法的準確性以及在程式化和計算效率方面的優勢．