祝世興，李佩窈，祝恒佳，魏戩

(1.中國民航大學航空工程學院，天津 300300；2.昆明理工大學機電工程學院，云南昆明 650000)

0 前言

緩沖器是所有現代飛機不可或缺的部分，目前油氣緩沖器在現代飛機上應用最為廣泛[1-2]。雙腔式油氣緩沖器包含低壓腔和高壓腔2個空氣腔，擁有單腔緩沖器優(yōu)點的同時在應對粗糙路況時的吸能特性也得到明顯改善。雙腔式油氣緩沖器剛度特性對改善飛機在沖擊振動、高速滑跑時的運行品質，提升乘客的舒適度，保護飛機零部件以及儀表設備等方面影響較大。目前，緩沖器的剛度特性一般采用激振試驗的方法獲得，但是由于試驗自身的滯后性、周期長以及成本高等缺點，采用仿真模擬的方法精確地模擬其動態(tài)特性成為研究趨勢[3-5]。雙腔式油氣緩沖器內部結構復雜，一般的仿真模擬并沒有考慮到流體在工作過程中的流動損失對緩沖器整體剛度特性的影響。因此，為更接近緩沖器實際工作情況，采用流固耦合仿真方法對雙腔式油氣緩沖器整體的剛度特性進行研究是非常重要的。

早期的雙腔式緩沖器多應用于直升機上，并且大部分是對起落架整體進行研究，而針對雙腔式緩沖器的單獨研究較少。文獻[6]中研究了符合直升機地面共振穩(wěn)定性要求的起落架剛度及阻尼的優(yōu)化設計方法。文獻[7]中利用動力學分析軟件ADAMS對起落架進行了落震動力響應的模擬與分析。文獻[8]中對大柔性飛機起落架緩沖器參數進行了優(yōu)化設計。文獻[9]中利用ADAMS建立了柔性活塞桿的緩沖器動力學模型并進行了落震仿真模擬，探討活塞桿剛度對緩沖器吸能特性的影響。大多數研究中均將起落架簡化為二自由度彈簧質量塊模型，將緩沖器簡化為所需力的形式加載到質量塊上，但這并未考慮到起落架結構的彈性特性和緩沖器的空氣彈簧和阻尼的非線性特性。因此，文獻[10]中利用哈密頓變分原理，建立了彈性結構的起落架動力學模型，以彌補上述缺點。對于直升機起落架的防“地面共振”設計和抗墜毀能力，文獻[11]中設計了一種單級雙腔緩沖器，其中高壓腔與低壓腔分離，并且保留了單級雙腔緩沖器的優(yōu)點，串聯雙腔緩沖器目的是將飛機的起飛著陸和耐墜毀緩沖行程進行分離。文獻[12]中利用等效線性化方法，分析了串聯雙腔緩沖器在多種頻率下的緩沖器等效剛度和當量阻尼分析值。文獻[13]中利用流固耦合及動網格方法，對緩沖器的阻尼特性進行了仿真。但是到目前為止，仍缺少考慮緩沖器復雜的流固耦合問題尤其是緩沖器工作過程中流體的流動損失有限元模擬問題對緩沖器剛度特性影響的相關算例與研究。

本文作者以某大型民航客機主起落架為例，建立適用于ABAQUS的雙腔油氣緩沖器仿真模型。利用工程估算方法模擬緩沖器沖擊過程并與ABAQUS中的仿真結果進行比對，驗證ABAQUS中模型的可用性；采用正弦激勵，分析緩沖器隨幅值和頻率變化時的剛度特性；采用隨機激勵，模擬實際工況，分析緩沖器的動態(tài)特性，從而為雙腔式緩沖器結構性能的優(yōu)化提供依據，并為進一步利用ABAQUS軟件開發(fā)緩沖器樣機設計奠定基礎。

1 雙腔緩沖器模型描述

如圖1所示，雙腔式緩沖器低壓腔與高壓腔通常充入一定初始氣壓的干燥氮氣或空氣，工作時活塞運動使油腔中的油液流向低壓腔壓縮氣體，同時吸收能量。油液在外力的作用下，流經油孔耗散掉起落架一部分能量。當低壓腔的空氣壓力與油液阻尼力增加到等于高壓腔初始壓力時，位于高低壓腔之間的活塞開始移動，高壓腔開始被壓縮[14]。壓縮緩沖器使得活塞上移稱為“正行程”。經過初始沖擊后，空氣腔氣體作為彈性體開始膨脹釋放能量，迫使油液又反流回油腔，活塞向下移動，此過程稱為“反行程”。經過一正一反2個行程完成一個循環(huán)，經多個循環(huán)，可耗散掉絕大部分飛機降落時產生的能量，進一步達到減震的效果。飛機的著陸撞擊時間在0.6～0.9 s之間[15]。正反行程總時間不超過0.8 s。

圖1 緩沖器結構

2 緩沖器數學建模

緩沖器軸向力Fs是由空氣彈簧力Fa、油液阻尼力Fh、緩沖器結構限制力Fl以及內部摩擦力Ff構成[16]。

Fs=Fa+Fh+Fl+Ff

(1)

結合氣體多變方程式可以得到空氣彈簧力的力學方程為

(2)

利用MATLAB擬合緩沖行程和活塞行程為

s=8.582×104t5-1.064×105t4+4.769×

104t3-1.025×104t2+1 180t

(3)

s′=1.276×104t4-1.203×104t3+3 180t2-124t

(4)

油液阻尼力是油液流過油孔時，油孔兩側存在壓力差而產生的[17]。假設油液為不可壓縮流體，由Bernoulli方程和連續(xù)性方程[18]可以得到油液阻尼力為

(5)

其中：Ah是緩沖器有效壓油面積；Aoil是油孔截面積；Cd是油孔縮流系數；ρ是油液密度。

緩沖器內部庫侖摩擦力假定為一個由緩沖器彎曲位移引起的內部摩擦力和內部壓力引起的內部摩擦力的合力，表達式為

Ff=Ff1+Ff2

(6)

緩沖器彎曲位移引起的內部摩擦力為

(7)

其中：μb是庫侖摩擦因數；Nu和NL分別是活塞對缸體的作用力與柱塞對缸體的作用力。

緩沖器內部壓力引起的內部摩擦力為

(8)

其中：Db是軸套直徑；Hb是軸套高度。

緩沖器的設計要求指出:緩沖器系統應適當保留某些行程余量，吸收給定過載下的使用功量[1]。只有在極端條件下，才允許緩沖系統吸收最大的功量。正常使用情況下,緩沖器結構限制力為0，只有在過度伸長與極限壓縮的情況下才考慮結構限制力。文中仿真緩沖器正常的使用情況，所以結構限制力為0。

緩沖器的模型參數如表1所示。

表1 緩沖器模型參數

3 有限元模型建立

緩沖器在工作過程中，大致分為兩類接觸問題：第一類是缸體與活塞和柱塞的摩擦接觸，ABAQUS中對于此類接觸問題常采用拉格朗日法進行計算；第二類是氮氣和油液與缸體、活塞和柱塞的接觸，應采用流固耦合的方法進行計算。對于流固耦合問題，ABAQUS中通常采用歐拉-拉格朗日(CEL)接觸算法來模擬流體與固體邊界處高速的相互作用。CEL方法中固體采用拉格朗日單元、流體采用歐拉單元，歐拉單元可以克服流體大變形下網格嚴重扭曲的現象，并且與拉格朗日單元不同之處在于:拉格朗日單元隨材料移動而移動，歐拉單元是材料在單元內流動。如圖2所示。

圖2 連續(xù)體變形

在默認情況下，歐拉體內部是沒有任何材料的，因此采用CEL方法的關鍵在于材料的填充，應在Load模塊下對相應的歐拉體區(qū)域進行操作。該模型材料的屬性定義如表2所示，其中Eos: US-UP為Mie-Grüneisen狀態(tài)方程中定義流體的材料參數。

表2 油氣材料屬性定義

在通常情況下，建立油氣緩沖器有限元模型時需要考慮氮氣腔、油腔、柱塞、活塞、缸體等，為提升計算速度，應將模型進行必要的簡化。假定油液密封性良好，將密封圈簡化為柱塞的一部分，不予以考慮,去除不必要的轉角與溝槽,忽略緩沖器自身重力的影響等。

該模型中，緩沖器柱塞與缸體采用C3D8R六面體單元，氣體與液體采用EC3D8R單元，進行材料的填充。利用ABAQUS/CAE建立緩沖器三維實體模型，而后建立離散單元模型，離散剛體單元有41 339個，其中C3D8R有26 438個、EC3D8R有14 901個，如圖3所示。

圖3 網格劃分

應用ABAQUS/Explicit模塊進行動態(tài)分析，緩沖器有限元模型計算過程分為2個部分：(1)在特定的頻率與幅值下給予上端一個沖擊位移，模擬緩沖器的沖擊過程，輸出隨時間變化的緩沖器軸向力曲線；(2)給予緩沖器不同的動態(tài)激勵，模擬它的運動狀態(tài)與工作特性。

4 模型仿真與結果分析

有限元模型建立后，利用軟件按照設定的要求進行仿真分析，直到設定的仿真時間為止。緩沖器的正反行程總時間不超過0.8 s，模型設定的仿真時間到0.5 s為止。

4.1 沖擊過程

給予緩沖器上端240 mm的垂向位移，模擬緩沖器沖擊過程。圖4和圖5分別為壓縮過程中緩沖器的內部流速與內部液體壓力。

圖4 緩沖器的內部流速

圖5 緩沖器內部液體壓力

基于理論基礎建立數學模型，采用MATLAB軟件編寫相應的計算程序，得到空氣彈簧力、油液阻尼力、緩沖器軸向力隨時間變化的曲線，并將ABAQUS仿真結果與MATLAB結果進行比對，分析二者之間的誤差。

緩沖器正行程時間約為0.23 s，圖6所示是空氣彈簧力的理論值與仿真值隨時間變化曲線?？芍涸趖=0.13 s之前只有低壓腔被壓縮，t=0.13 s后高壓腔處于壓縮階段;在t=0.19 s后，由于空氣腔的氣體在高壓狀態(tài)下液化導致氣體溫度升高，內能增加并伴隨了一定的動能損失，所以在t=0.19 s之后理論值與仿真值產生了較大誤差。

圖6 空氣彈簧力隨時間變化曲線 圖7 油液阻尼力隨時間變化曲線

圖7所示為油液阻尼力理論值與仿真值隨時間變化曲線?？芍河鸵鹤枘崃﹄S時間的增加先增加后減小，這是因為緩沖器內部回油孔前后壓差處于先增加后減小的趨勢，使得油液流速處于先增加后減慢的情況。

數值仿真和有限元仿真的誤差可以通過均方根值來定義：

(9)

式中：m表示離散力信號的數據個數；FNUM和FSIM分別表示數值仿真和有限元仿真的輸出力。

表3所示為空氣彈簧力和油液阻尼力的數值仿真和有限元仿真之間的誤差大小，可以看出：空氣彈簧力的誤差較大，但由于空氣彈簧力的平均值為17 609 N，因此誤差是可以接受的，驗證了文中建立的緩沖器有限元模型的可用性。

表3 不同仿真模型之間的誤差 單位：N

圖8所示為緩沖器軸向力理論值與仿真值隨時間變化的曲線，可得：理論值與仿真值相對誤差為0.96%，從而驗證了仿真模型的準確性。

圖8 緩沖器軸向力隨時間變化曲線

經比對，空氣彈簧力、油液阻尼力和軸向力仿真值與理論值誤差均小于1%，證明了該模型的可用性，為后續(xù)研究緩沖器特性提供了參考。

4.2 不同激振對緩沖器等效動剛度影響

在緩沖系統運行過程中，外界環(huán)境溫度與充氣壓力的改變會直接影響緩沖系統的性能。初始擾動的大小會影響起落架的振動幅值，飛機裝載不同會引起機體轉動慣量的變化，從而導致起落架固有頻率的改變[19]。因此，本文作者給予緩沖器不同位移激勵方式，改變頻率與幅值，觀察其變化情況。

圖9所示為沖擊激勵下緩沖器內部流體壓力，顯示了緩沖器在工作過程中油液流動時的噴濺現象。由于空氣壓縮的速度很快，以至于在壓縮過程中空氣所產生的熱量來不及向外界傳遞，因此空氣壓縮的過程可近似看成絕熱過程。但油液高速流動時會發(fā)生噴濺現象，與低壓腔的氣體進行強烈的熱交換，因此空氣壓縮的過程為多變過程。圖9也從側面說明了在壓縮過程中低壓腔的多變過程。

圖9 沖擊激勵下緩沖器內部流體壓力

設對緩沖器施加的位移激勵為簡諧激勵，表達式為

x(t)=H0+Hsinωt

(10)

其中：ω為激振力的圓頻率，ω=2πf；f為激勵頻率；H為激勵幅值；H0為初始壓縮量。

根據雙腔式磁流變減振器的軸向力與位移功量圖，結合式(11)求得緩沖器的等效動剛度K為

(11)

假定緩沖器的振動幅值保持12 mm不變，緩沖器等效動剛度隨振動頻率的變化如圖10所示?？梢钥吹剑弘S頻率的增加等效動剛度降低，當頻率超過8 Hz時，等效動剛度降低速度減慢。

圖10 頻率對等效動剛度的影響

假定緩沖器的振動頻率保持2 Hz不變，以不同振幅的正弦位移信號作為輸入，圖11所示為緩沖器在激振幅值為6 mm時的力-位移遲滯曲線，其中位移為正時表示壓縮，為正行程；位移為負時表示拉伸，為負行程。圖12所示為幅值對等效動剛度的影響,可以看出：幅值與等效動剛度呈現嚴重的非線性特性，隨振動幅值增加等效動剛度減??；與頻率對等效動剛度的影響相比，幅值的變化對等效動剛度的影響較大。

圖11 激振幅值為6 mm時的力-位移曲線 圖12 幅值對等效動剛度的影響

然而，正弦激勵在實際工作過程中很少遇到，作為隔振元件的油氣緩沖器被廣泛應用于隨機振動的環(huán)境中。隨機信號可以通過離散傅里葉變換分解為一系列具有不同頻率和幅值的簡諧波。諧波疊加法基本原理是將路面的不平度表示為大量具有隨機相位的正弦或余弦諧波之和[20]。

在ABAQUS中，可以考慮采用周期形幅值曲線來近似代替隨機激勵，以產生偽隨機激勵曲線。周期形幅值曲線用傅里葉級數表示，表達式為

(12)

式中：N為傅里葉級數項的個數；t0為起始時刻；A0為初始幅值；An為cos項的系數；Bn為sin項的系數。

取A0=0.01、t0=0，頻率與幅值按照表4取值。

表4 偽隨機激勵曲線幅值與頻率部分取值

圖13所示為偽隨機位移激勵時域曲線，圖14所示為0.2～1 s緩沖器輸出的軸向力。選取此時間段是整個仿真過程中的隨機選擇，以方便觀察軸向力變化趨勢。由圖14可以看出：有限元方法輸出緩沖器軸向力略小于數值方法輸出的緩沖器軸向力，原因是油液與筒壁之間產生摩擦從而引起沿程損失，但總體趨勢一致。

圖13 偽隨機位移激勵時域曲線 圖14 0.2～1 s內的軸向力

5 結論

通過建立起落架雙腔緩沖器的數學模型與有限元模型，模擬不同工況下的緩沖器工作特性，得到以下結論：

(1)與數學模型相比，基于流固耦合有限元模型可以較好地處理油液在流動過程中產生的壓力損失，與緩沖器的實際工況更為接近；有限元仿真分析結果與數學模型計算結果在大位移沖擊作用下相對誤差均小于1%，進一步說明了有限元模型的可行性；

(2)在正弦激勵工況下，隨激振頻率與幅值的改變，緩沖器等效動剛度表現出分段非線性動力學特性，隨頻率與幅值的增加等效動剛度減小，并且幅值的改變對等效動剛度產生的影響較大；

(3)在隨機激勵工況下，由于油液的流動損失導致了仿真結果小于理論結果，但仿真結果的總體趨勢與理論結果基本保持一致，吻合度較高。