水平著陸飛行器著陸架系統耦合壓縮觸地過程仿真分析

劉 剛，鄭宏濤，李 洋

(中國運載火箭技術研究院研究發展中心， 北京 100076)

0 引言

水平著陸飛行器著陸架觸地是一個短時間、高動態的復雜運動過程。為了實現飛行器安全著陸和觸地，保證飛行器在此過程中不后翻、尾部不擦地、輪胎和減震器壓縮量不超限，需要對觸地過程進行精確的仿真分析，才能確定或優化相關的設計參數。著陸架觸地過程中輪胎和減震器的壓縮運動存在耦合因素，因此輪胎壓縮量及其受力的求解成為問題的關鍵環節。高澤迥[1]給出了著陸架減震器和輪胎的力學模型，但是未給出壓縮量的求解方法。羅琳胤等[2]采用商業軟件對起落架多體動力學進行建模與分析，但是軟件中的壓縮量求解方法未知。陳麗城等[3]采用Matlab/Simulink實現了無人機地面動力學建模及分析，但是其模型中不包含減震器和輪胎壓縮量的耦合因素。針對水平著陸飛行器著陸架系統耦合壓縮觸地過程仿真分析問題，本文采用一種基于瞬時平衡條件的求解方法，先求解減震器壓縮量，再求解輪胎壓縮量，進而求得輪胎法向力，實現水平著陸飛行器觸地過程的仿真分析。

1 著陸架系統耦合壓縮模型

飛行器搖臂式著陸架系統示意圖如圖1所示。其中AB為懸臂梁，與飛行器機體固連。ACD為搖臂，在A與懸臂梁鉸接，D為輪軸。BC為減震器，在B與懸臂梁鉸接，在C與搖臂鉸接。

圖1 搖臂式著陸架系統示意圖Fig.1 Schematic figure of a rock arm undercarriage system

fa=fa(S)

(1)

(2)

(3)

輪胎法向力N是輪胎壓縮量δ的函數,滾動摩擦力f與法向力成線性關系，如式(4)、式(5)所示，式中cf為輪胎的滾動摩擦系數。

N=N(δ)

(4)

f=cfN

(5)

在輪胎下邊沿與著陸場地面接觸后，飛行器開始觸地過程。由于飛行器在觸地時具有一定的垂直方向的下沉速度，會造成輪胎迅速壓縮，輪胎壓縮后輪軸受到支持力和滾動摩擦力。這兩個力對搖臂產生一個對A點的逆時針力矩，當這個力矩增大到一定程度時，搖臂會逆時針轉動。搖臂的逆時針轉動一方面會造成減震器壓縮；另一方面會改變著陸架系統的幾何構型，使輪軸D的高度增大。輪軸D的高度增大會造成輪胎壓縮量減小。因此著陸架的減震器和輪胎的壓縮量是耦合的。

2 基于瞬時平衡條件考慮著陸架系統耦合壓縮的觸地過程仿真分析方法

2.1 瞬時平衡條件推導

搖臂的運動過程滿足繞A點轉動的微分方程：

(6)

其中，JACD為搖臂的轉動慣量，ωACD為其轉動角速度，∑MA為搖臂受到的繞A點的力矩之和。

通過直接求解這個微分方程對著陸架系統進行仿真具有一定的困難。一方面，壓縮量的耦合因素使得每一步仿真中的受力和力矩難以確定；另一方面，如果直接求解微分方程，由第k-1步遞推到第k步時，可能出現減震器受摩擦力和阻尼力與其運動趨勢不符合的情況，即使仿真步長比飛行器六自由度仿真步長低一個數量級，求解結果也會發生振蕩或跳變。為解決此問題，本文采用一種基于瞬時平衡條件的著陸架系統耦合壓縮量和受力的仿真分析方法。

對上述搖臂的運動方程進行分析。一般情況下，搖臂的轉動慣量與角加速度的乘積相對于來自輪胎和減震器的力矩是小量，其數值一般低3個以上數量級。因此可忽略此小量，即得到搖臂的瞬時平衡公式：

∑MA=0

(7)

2.2 著陸架減震器運動趨勢判斷方法

搖臂所受的對A點的力矩包括5部分：輪胎法向力N對應的力矩MN，輪胎滾動摩擦力f對應的力矩Mf，減震器彈簧力fa對應的力矩Mfa，減震器摩擦力ff對應的力矩Mff，減震器阻尼力fd對應的力矩Mfd。由于減震器摩擦力和阻尼力的方向與減震器運動方向相反，因此需要先判斷減震器的運動方向，再分不同的情況以不同的方法計算合力矩。

在第k步仿真中，來自飛行器六自由度運動方程的相關參數(包括飛行器質心高度hk，俯仰角φk，滾轉角φk，質心在飛行器的坐標Xgk、Ygk、Zgk)被更新。此時，假設減震器壓縮量暫時維持與第k-1步的值不變，則根據第k-1步仿真中輪軸D的坐標Xdk-1、Ydk-1、Zdk-1和減震器壓縮量Sk-1以及相關幾何關系可計算出一個輪胎壓縮量δtemp：

δtemp=δtemp(hk,φk,φk,Xgk,Ygk,Zgk,

Xdk-1,Ydk-1,Zdk-1,Sk-1)

(8)

3)減震器壓縮量不變情況：若不滿足上述兩種情況，則減震器壓縮量維持第k-1步的值不變。

2.3 減震器和輪胎的壓縮量耦合求解方法

分別針對上述3種情況，給出第k步仿真中減震器和輪胎的壓縮量的耦合求解方法。

(1)減震器壓縮量不變情況

此種情況最簡單，壓縮量按式(9)、式(10)計算：

Sk=Sk-1

(9)

δk=δtemp

(10)

(2)減震器壓縮量增大情況

(11)

其中,step為仿真步長，與飛行器六自由度仿真步長相同。

如圖1所示，減震器壓縮量改變后，著陸架系統的幾何關系將發生改變，C點變為C1，D點變為D1。因此需要根據S*k的值重新計算幾何關系，計算出輪軸D1的坐標Xd*k、Yd*k、Zd*k，重新計算輪胎壓縮量δ*k：

δ*k=δ*k(hk,φk,φk,Xgk,Ygk,Zgk,

Xd*k,Yd*k,Zd*k,S*k)

(12)

對A點的和力矩為：

∑M*A=M*N+M*f-(M*fa+M*ff+M*fd)

(13)

注意到，每對應一個假定的減震器壓縮量ΔS*,可得到一個∑M*A=∑M*A(ΔS*)，這是一個以ΔS*為自變量的一元函數。按前文推導的瞬時平衡條件可得：

∑M*A(ΔS*)=0

(14)

式(14)是一個以ΔS*為自變量的一元非線性方程，可用割線法[4]等數值法求解，解方程時需要注意限制自變量的取值在合理范圍。

其解即為第k步的減震器壓縮量增量ΔSk，則第k步減震器壓縮量為：

Sk=Sk-1+ΔSk

(15)

根據Sk重新更新幾何關系，包括輪軸的坐標Xdk、Ydk、Zdk，各桿的長度和各角的角度等，進而可計算出第k步的輪胎壓縮量δk。

δk=δk(hk,φk,φk,Xgk,Ygk,Zgk,

Xdk,Ydk,Zdk,Sk)

(16)

(3)減震器壓縮量減小情況

此種情況求解方法與減震器壓縮量增大情況類似，唯一不同的是求解的瞬時平衡方程如式(17)：

∑MA(ΔS) =Mfa-(Mff+Mfd+MN+Mf)=0

(17)

2.4 飛行器觸地過程仿真分析方法

根據2.3節的方法分3種情況可求得第k步的減震器壓縮量Sk和輪胎壓縮量δk。根據輪胎法向力與輪胎壓縮量的函數關系或插值表可求得輪胎法向力Nk。由滾動摩擦系數可求解滾動摩擦力fk。在仿真第k步，把飛行器視為一個整體，將Nk、fk及其對飛行器質心的力矩加入飛行器六自由度運動方程并進行求解，即可實現考慮著陸架系統耦合壓縮的觸地過程仿真分析。

3 考慮著陸架系統耦合壓縮的觸地過程仿真分析算例

采用本文給出的方法，選取觸地下沉速度-1.51m/s、觸地俯仰角11.2°、觸地滾轉角-0.02°的典型工況，對水平著陸飛行器的觸地過程進行仿真分析，仿真步長取1ms。仿真結果如圖2～圖9。

圖2 左減震器壓縮量Fig.2 Compressing amount of the left shock absorber

圖3 右減震器壓縮量Fig.3 Compressing amount of the right shock absorber

圖4 左主輪壓縮量Fig.4 Compressing amount of the left main tire

圖5 右主輪壓縮量Fig.5 Compressing amount of the right main tire

圖6 飛行器俯仰角Fig.6 Pitch angle of the vehicle

圖9 飛行器尾尖與地面距離Fig.9 Distance between the vehicle tail and ground

由如圖2～圖9可見，左、右主輪都被彈起一次，并進行了二次觸地。減震器壓縮量和輪胎壓縮量都經歷了從零到最大再回落并穩定的過程。飛行器的俯仰角、滾轉角和下沉率都逐漸收斂于零附近，并保持穩定。在觸地過程中飛行器未發生后翻、尾部擦地、輪胎或減震器壓縮量超限現象。

4 結論

針對水平著陸飛行器觸地過程六自由度運動與著陸架減震器和輪胎的耦合壓縮運動聯合仿真問題，將復雜的著陸架的機構多體動力學進行簡化，推導了搖臂式著陸架的力矩瞬時平衡條件，給出了著陸架減震器運動趨勢判斷方法，通過求解以減震器壓縮量增量為自變量的一元非線性方程求得減震器壓縮量，進而求得輪胎壓縮量和輪胎受力，實現飛行器觸地過程仿真分析。該方法可在飛行器六自由度仿真模型基礎上直接擴展，無需改變其仿真步長，著陸架模型求解無需依賴商業軟件，可對水平著陸飛行器觸地過程相關參數和性能進行有效的仿真評估和驗證。

參考文獻

[1] 高澤迥.飛機設計手冊第14冊：起飛著陸系統設計[M].北京：航空工業出版社，2002.

[2] 羅琳胤，邊寶龍.飛機起落架緩沖性能仿真分析[J].機械設計，2012,29(4)：56-58+62.

[3] 陳麗城, 李春濤, 張孝偉, 等. 無人機地面動力學建模及分析[J]. 計算機仿真, 2016 (6): 13-18.

[4] 顏慶津.數值分析[M].北京：北京航空航天大學出版社，2000.