基于參變量變分原理的直升機(jī)系留載荷高性能計(jì)算方法

張燕輝， 高 強(qiáng)

(1.中國(guó)直升機(jī)設(shè)計(jì)研究所，景德鎮(zhèn) 333001；2.大連理工大學(xué) 工程力學(xué)系 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024)

1 引 言

當(dāng)直升機(jī)在艦船上停留或者在機(jī)庫(kù)長(zhǎng)期存放時(shí)，通常要承受因艦船運(yùn)動(dòng)產(chǎn)生的慣性力、直升機(jī)重力作用和風(fēng)載等載荷作用，這些載荷可能會(huì)引起直升機(jī)在艦船上發(fā)生側(cè)翻或滑移，從而使直升機(jī)受損。為了保證直升機(jī)安全地停放在艦船上，就必須使用系留裝置將直升機(jī)系留在艦船上，這就需要對(duì)直升機(jī)的系留索具和系留接頭上的載荷進(jìn)行分析，從而給出合理的系留方案，使直升機(jī)和系留設(shè)備上承受的載荷比較均勻和合理。

分析艦載直升機(jī)的系留系統(tǒng)，不僅要判別索具的拉伸和松弛狀態(tài)，而且還需考慮起落架的變形。當(dāng)索具受拉伸荷載作用時(shí)，其與普通桿件沒有區(qū)別；而當(dāng)索具受壓縮荷載作用時(shí)，因無(wú)法承受壓力作用出現(xiàn)松弛而失效。對(duì)于起落架，當(dāng)其受到壓縮載荷作用時(shí)，與普通桿件沒有區(qū)別；而當(dāng)起落架受到拉伸荷載作用時(shí)，起落架離開地面不承力而失效。因此，索具和起落架分別在受拉和受壓時(shí)表現(xiàn)出的強(qiáng)非線性特征，是直升機(jī)系留系統(tǒng)的一個(gè)顯著特點(diǎn)。同時(shí)，由于索具和起落架的變形，將導(dǎo)致整個(gè)系統(tǒng)發(fā)生較大變形和轉(zhuǎn)動(dòng)，這意味著直升機(jī)系留問題還必須考慮結(jié)構(gòu)大變形。一般情況下，系留索具的數(shù)量往往要多于直升機(jī)機(jī)體的運(yùn)動(dòng)自由度，直升機(jī)系留問題本質(zhì)上屬于超靜定問題。對(duì)于直升機(jī)系留問題表現(xiàn)為材料和幾何雙重非線性特征，需要針對(duì)這類問題發(fā)展有效的求解算法。目前關(guān)于系留載荷的計(jì)算方法主要有矩陣力法、有限元法和能量法等。早期關(guān)于系留載荷的計(jì)算通常采用矩陣力法[1]求解，但需要先判斷索具拉壓狀態(tài)，采用矩陣力法求解較繁瑣，且不便于實(shí)現(xiàn)程序化。隨著有限元法在飛機(jī)設(shè)計(jì)過(guò)程中的廣泛應(yīng)用，利用有限元法計(jì)算系留載荷已成為重要手段，而且通常借助成熟的商業(yè)軟件對(duì)系留問題進(jìn)行模擬，如NASTRAN[2]，ABAQUS[3]和ANSYS[4]等。主要方法是，機(jī)身采用剛體單元模擬，索具采用拉壓剛度不同的桿單元模擬，而起落架采用間隙單元模擬。盡管利用有限元法計(jì)算系留載荷能夠獲得相對(duì)較準(zhǔn)確的計(jì)算結(jié)果，但在一些特殊情況下會(huì)發(fā)生不收斂現(xiàn)象。另一類求解系留載荷的重要方法是能量法，該方法主要思想是利用能量原理建立平衡方程后，采用迭代方法進(jìn)行求解。文獻(xiàn)[5]利用能量原理對(duì)系留索具載荷進(jìn)行求解，但并沒有考慮起落架的彈性變形和機(jī)身的轉(zhuǎn)動(dòng)。在此基礎(chǔ)上，文獻(xiàn)[6]對(duì)其進(jìn)行了改進(jìn)，考慮了起落架的彈性變形，而且增加了機(jī)身轉(zhuǎn)動(dòng)自由度。文獻(xiàn)[7]基于虛功原理推導(dǎo)了機(jī)體剛體位移方程組，并采用牛頓迭代法求解非線性方程組?；谀芰吭淼挠?jì)算方法在求解直升機(jī)系留載荷方面有兩個(gè)難點(diǎn)。一方面，不容易判斷索具和起落架的拉壓狀態(tài)；另一方面，由于該問題具有材料和幾何雙重非線性，計(jì)算需要不斷迭代，容易發(fā)生不收斂。

目前關(guān)于我國(guó)艦載直升機(jī)系留問題，大部分采用ABAQUS和NASTRAN軟件進(jìn)行分析，存在較多困難。首先，直升機(jī)在艦船停放時(shí)，會(huì)受到艦船運(yùn)動(dòng)的慣性力、直升機(jī)重力和風(fēng)載等載荷作用，因此不同的系留方案、直升機(jī)重量狀態(tài)、航速、浪向和風(fēng)載方向，組合后的工況數(shù)量較多，可達(dá)幾萬(wàn)甚至幾十萬(wàn)種工況，當(dāng)采用有限元軟件批處理分析時(shí)，需要不斷讀取軟件產(chǎn)生的結(jié)果數(shù)據(jù)文件，會(huì)耗費(fèi)較多的計(jì)算時(shí)間。每更新一次系留方案，幾天甚至幾十天才能給出計(jì)算結(jié)果，嚴(yán)重影響任務(wù)進(jìn)度。其次，ABAQUS和NASTRAN軟件在計(jì)算非線性問題時(shí)，在一些工況下會(huì)發(fā)生不收斂，因此在進(jìn)行批處理時(shí)，對(duì)于不收斂工況，需要挑選出來(lái)后將間隙單元改成桿單元或彈簧單元計(jì)算，從而導(dǎo)致計(jì)算結(jié)果與實(shí)際偏離較大，而且由于該過(guò)程需要挑選不收斂工況，處理過(guò)程復(fù)雜繁瑣，計(jì)算效率低，無(wú)法滿足目前的型號(hào)需求。最后，由于需要不斷讀取軟件計(jì)算結(jié)果文件，效率較低，因而無(wú)法對(duì)系留方案進(jìn)行優(yōu)化設(shè)計(jì)。本文針對(duì)以上難點(diǎn)，基于參變量變分原理提出了計(jì)算系留載荷的新方法。該方法具有高效、收斂性好和方便實(shí)現(xiàn)程序化等優(yōu)點(diǎn)。

參變量變分原理是由鐘萬(wàn)勰等[8]提出的，并廣泛應(yīng)用于彈塑性[9]、接觸[10]和摩擦[11]等非線性問題分析，在求解拉壓剛度不同的桿件系統(tǒng)時(shí)非常有效[12]。艦載直升機(jī)的系留問題可簡(jiǎn)化為由機(jī)身、索具和起落架組成的拉壓剛度不同的桿件系統(tǒng)。本文將機(jī)身簡(jiǎn)化為剛體，并利用參變量變分原理建立了機(jī)體重心處的平衡方程組，同時(shí)考慮了索具和起落架的材料非線性以及結(jié)構(gòu)大變形。利用參變量變分原理引入?yún)⒆兞亢退沙谧兞?，從而?zhǔn)確地判斷索具和起落架的拉壓狀態(tài)，將索具和起落架拉壓剛度不同的材料非線性問題轉(zhuǎn)換為線性互補(bǔ)問題求解，從而極大地提高了結(jié)果的收斂性。

2 直升機(jī)系留模型簡(jiǎn)介

艦載直升機(jī)系留方式如圖1所示，系留系統(tǒng)由機(jī)身、索具和起落架組成。索具一端與機(jī)身系留點(diǎn)連接，另一端與艦船甲板上的系留窩連接；起落架一端與機(jī)身連接，另一端與艦船甲板接觸。機(jī)身簡(jiǎn)化為剛體，索具和起落架可通過(guò)拉壓剛度不同的桿單元模擬。

圖1 直升機(jī)系留

索具只能承受拉伸荷載作用，當(dāng)索具處于拉伸狀態(tài)時(shí)，其受軸向拉力；當(dāng)索具松弛或未變形時(shí)，索具不受力。起落架和索具的重量相對(duì)于直升機(jī)重量非常微小，可忽略不計(jì)。同時(shí)假設(shè)起落架可以沿著水平方向運(yùn)動(dòng)，且不考慮摩擦力影響。設(shè)索具的數(shù)量為N，則第i(i=1,2,…,N)根索具的本構(gòu)關(guān)系如圖2(a)所示，可表示為

(1)

(2)

圖2 索具和起落架的本構(gòu)關(guān)系

在風(fēng)浪較大時(shí)，直升機(jī)系留裝置可能會(huì)受到較大載荷，導(dǎo)致索具和起落架的變形較大。因此對(duì)于直升機(jī)系留問題，不僅要考慮結(jié)構(gòu)本身的材料非線性，而且還要考慮因結(jié)構(gòu)發(fā)生大變形而導(dǎo)致的幾何非線性，這表明直升機(jī)的系留問題表現(xiàn)出較強(qiáng)的非線性特點(diǎn)，結(jié)果很容易發(fā)生不收斂現(xiàn)象。本文在考慮結(jié)構(gòu)大變形的情況下，利用參變量變分原理，將拉壓剛度不同桿件的材料和幾何雙重非線性問題轉(zhuǎn)換為僅考慮幾何非線性的互補(bǔ)問題，從而極大地提高了結(jié)果的收斂性。

3 基于參變量變分原理的系留載荷分析

3.1 拉壓剛度不同桿件的參變量變分原理

對(duì)于任意一根拉壓剛度不同的桿件(不考慮索具預(yù)緊力和初始變形)，假設(shè)桿件的拉伸和壓縮剛度分別為K(+)和K(-)，且K(+)≠K(-)，則桿件的本構(gòu)關(guān)系為

F=KΔl

(3)

式中F為桿件軸向力，Δl為桿件變形量，且

(4)

當(dāng)K(+)>K(-)時(shí)，本構(gòu)關(guān)系(3)可表示為

F=K(+)(Δl+λ)

(5)

式中λ為參變量，

(6)

式(6)與互補(bǔ)關(guān)系(7)等價(jià)，即

(K(-)-K(+))Δl-K(+)λ+ν=0

(λ≥0,ν≥0,λν=0) (7)

而當(dāng)K(+)

F=K(-)(Δl-λ)

(8)

則要求

(9)

式(9)與互補(bǔ)關(guān)系(10)等價(jià)，即

(K(-)-K(+))Δl-K(-)λ+ν=0

(λ≥0,ν≥0,λν=0)

(10)

上述兩種情況下，對(duì)應(yīng)的方程(5,7,8,10)可統(tǒng)一表示為

F=Kmax(Δl-sλ)

sKmax-KminΔl-Kmaxλ+ν=0

(λ≥0,ν≥0,λν=0)

(11)

式中s=sign (K(-)-K(+))，Kmax=max (K(-),K(+))，Kmin=min (K(-),K(+))。符號(hào)sign表示取符號(hào)，其定義為

(12)

容易證明，式(11)給出的參變量變分原理與式(3,4)給出的平衡方程等價(jià)，具體證明過(guò)程參見文獻(xiàn)[13]。

3.2 直升機(jī)重心平衡方程推導(dǎo)

(13)

式中δX=[δx,δy,δz]T，且

(14)

(15)

(16)

(17)

通過(guò)參變量λi可判斷索具的拉壓狀態(tài)，當(dāng)λi=0時(shí)，索具處于拉伸狀態(tài)；當(dāng)λi>0時(shí)，索具處于松弛狀態(tài)。Fi沿著x，y和z方向的分量分別為

(18)

(19)

式中符號(hào)?表示兩個(gè)向量的叉積運(yùn)算。

(20)

(21)

(22)

(23)

起落架壓力對(duì)重心點(diǎn)的力矩為

Mj=Fj?(X′j-X0)

(24)

整個(gè)結(jié)構(gòu)的平衡方程可寫為

(25)

再考慮約束條件

(26)

(27)

式(26,27)分別對(duì)應(yīng)于索具和起落架的約束方程。

根據(jù)以上推導(dǎo)，直升機(jī)系留載荷計(jì)算最終歸結(jié)為求解帶有約束條件(26,27)的關(guān)于方程組(25)的二次規(guī)劃問題。由于式(25)為非線性方程組，需要通過(guò)迭代方法進(jìn)行求解。本文在求解時(shí)運(yùn)用MATLAB軟件的fmincon函數(shù)，計(jì)算時(shí)選擇序列二次規(guī)劃算法(SQP)進(jìn)行求解。通過(guò)帶有約束條件(25～27)即可求參變量λ，然后根據(jù)方程(16,21)可分別計(jì)算出索具的拉力和輪胎的壓力。

4 數(shù)值算例

以某艦載直升機(jī)在艦船上的系留問題為例，證實(shí)本文方法的性能。直升機(jī)質(zhì)量為6815 kg，利用機(jī)身上的10個(gè)系留點(diǎn)(A′,B′,C′,D′,E′，每側(cè)各一點(diǎn))布置14根系留索具，其中點(diǎn)B′和點(diǎn)D′，每個(gè)點(diǎn)2根；點(diǎn)A′、點(diǎn)C′ 和點(diǎn)E′每個(gè)點(diǎn)1根，艦船甲板上布置14個(gè)系留窩，機(jī)身左右兩側(cè)分別布置7個(gè)系留窩；起落架采用后三點(diǎn)布置，由左右兩個(gè)主起落架和一個(gè)尾起落架組成，系留方案如圖3所示，機(jī)身系留點(diǎn)、艦船甲板系留窩和起落架坐標(biāo)列入表1。直升機(jī)重心坐標(biāo)為(8512,1,3417)。索具的彈簧剛度為900 N/mm，主起落架的剛度為705.9 N/mm，尾起落架的剛度為437.5 N/mm。分別考慮3種不同的載荷工況，其對(duì)應(yīng)的載荷數(shù)據(jù)列入表2，載荷作用于重心處。

圖3 某型艦載直升機(jī)系留方案

為驗(yàn)證本文方法的性能，分別使用通用有限元軟件NASTRAN和ABAQUS建立全機(jī)系留模型進(jìn)行計(jì)算(不考慮預(yù)緊力和起落架初始?jí)嚎s)，并將計(jì)算結(jié)果與本文方法的計(jì)算結(jié)果進(jìn)行比較，三種方法在計(jì)算時(shí)直升機(jī)機(jī)身均簡(jiǎn)化為剛體。本文算例均是在中央處理器為12核、32 G內(nèi)存和主頻3.6 GHz的計(jì)算機(jī)上進(jìn)行。運(yùn)用NASTRAN軟件計(jì)算時(shí)，索具采用只承拉不承壓的桿單元來(lái)模擬，而起落架采用Gap單元模擬；運(yùn)用ABAQUS軟件計(jì)算時(shí)，索具和起落架均采用非線性連接單元(Connector單元)模擬索具和起落架非線性剛度特性。三種工況下的計(jì)算結(jié)果分別列入表3～表5。

表1 機(jī)身系留點(diǎn)、甲板系留窩和起落架坐標(biāo)(單位：mm)

表2 重心處作用三種載荷工況

表3 工況1：系留載荷計(jì)算結(jié)果(單位：N)

表4 工況2：系留載荷計(jì)算結(jié)果(單位：N)

表5 工況3：系留載荷計(jì)算結(jié)果(單位：N)

對(duì)于工況1，由表3可知，本文方法的計(jì)算結(jié)果與NASTRAN和ABAQUS軟件的計(jì)算結(jié)果符合良好，所有索具和起落架的拉壓狀態(tài)均一致，這證實(shí)了本文方法的正確性。同時(shí)，表3還表明本文方法與ABAQUS軟件的計(jì)算結(jié)果更接近，因?yàn)锳BAQUS軟件在求解非線性接觸問題時(shí)精度更高。對(duì)于工況2，NASTRAN軟件不收斂，無(wú)法進(jìn)行計(jì)算，而本文方法和ABAQUS軟件能夠很好地收斂，并在表4中給出了對(duì)應(yīng)的計(jì)算結(jié)果，可以看出結(jié)果符合良好；而對(duì)于工況3，ABAQUS軟件不收斂，無(wú)法進(jìn)行計(jì)算，而本文方法和NASTRAN軟件能夠很好地收斂，并在表5中給出了本文方法和NASTRAN軟件的計(jì)算結(jié)果。以上分析表明，ABAQUS和NASTRAN軟件在某些工況下可能不收斂，而本文方法具有較好的收斂性。

為了說(shuō)明本文方法的高效性，在均采用并行計(jì)算的情況下，隨著工況數(shù)量的增加，將本文方法與NASTRAN和ABAQUS軟件批處理計(jì)算效率進(jìn)一步比較。圖4給出了NASTRAN軟件批處理、ABAQUS軟件批處理和本文方法隨工況數(shù)量增加的CPU時(shí)間?？梢钥闯?，與NASTRAN軟件和ABAQUS軟件相比，本文方法具有非常高的計(jì)算效率，而且由于運(yùn)用NASTRAN和ABAQUS軟件批處理計(jì)算時(shí)需要讀取結(jié)果命令文件，隨著工況數(shù)量的增加，NASTRAN和ABAQUS軟件的計(jì)算時(shí)間急劇增加，而本文方法方便實(shí)現(xiàn)程序化處理，且不需要讀取其他軟件的數(shù)據(jù)文件。因此，本文方法在計(jì)算多工況的直升機(jī)系留載荷時(shí)具有非常高的效率。

圖4 NASTRAN軟件、ABAQUS軟件和本文方法的CPU時(shí)間

5 結(jié) 論

基于參變量變分原理建立了直升機(jī)系留載荷求解的有效方法，將由索具和起落架組成的拉壓不同剛度的材料和幾何雙重非線性問題，轉(zhuǎn)換為僅幾何非線性互補(bǔ)問題求解，極大地提高了結(jié)果的收斂性。數(shù)值算例中，通過(guò)與NASTRAN和ABAQUS軟件比較，結(jié)果表明,本文方法對(duì)于計(jì)算艦載直升機(jī)的系留載荷具有較好收斂性和較高計(jì)算效率。該方法對(duì)提高系留載荷計(jì)算效率、快速優(yōu)化系留方案和設(shè)計(jì)直升機(jī)機(jī)身系留接頭提供了更可靠的計(jì)算依據(jù)。