濾波白噪聲法的單輪起落架滑跑模型

雷繼超， 石鑫剛， 蔡良才， 王觀虎， 梁 磊

(空軍工程大學航空工程學院， 西安， 710038)

飛機對道面動載作用的影響因素主要有飛機自身的結構參數、道面平整度狀況、飛機滑跑速度等，對于飛機動載的研究主要是實測法。NASA[1]在飛機起落架安裝傳感器測量不平度激勵下起落架的荷載響應；許金余[2]對多種機型進行了動載測試，統計了不平度激勵下動載分布規律。理論建模法。研究者通過對功率譜密度與不平度的關系生成仿真道面，在對起落架模型適度簡化的基礎上求解動力學方程，進而求得飛機動載。孫璐建立了國際平整度指數(International Roughness Index，IRI)與路面功率譜密度的轉換關系[3]，并在假設機場道面平整度服從均值為0的高斯隨機過程基礎上，推導了功率譜密度與隨機動載的表達式[4-5]。張獻民[6]對文獻[3]提出的動載計算模型進行了空氣動力修正，定義了升力影響因子，并對各影響參數做了敏感度分析，但采用的升力影響因子僅是經驗公式，起落架和輪胎的仿真采用的是線彈性阻尼。虛擬樣機法。朱立國[7]采用虛擬樣機技術，對飛機滑跑過程進行了動力學分析，提供了一種道面動載計算的新思路。

1 PSD與IRI的轉換關系

機場道面關于不平度的研究主要參考公路平整度的評價方法，前期公路平整度主要采用功率譜密度(Power Spectral Density，PSD)Gq(n)表示：

Gq(n)=Gq(n0)(n/n0)-ω

(1)

式中：n為不平度的空間頻率，表示每米含有幾個波長；n0為空間參考頻率，一般取值為0.1 m-1；Gq(n0)為對應等級路面下的功率譜密度，詳細數值見文獻[8]；ω為頻率指數，一般取值為2。

綜合路面構造分類、路面平整度波長及其空間頻率范圍，確定了IRI(簡記I)與Gq(n0)關系式[9]：

(2)

代入式(1)，可得I與不平度換算關系，見表1。

表1 國際平整度指數分級標準與對應的PSD值

2 濾波白噪聲法仿真道面不平度的參數確定

根據隨機振動理論和線性系統響應規律，推導得到濾波白噪聲路面不平度時域模型[10]：

(3)

對上述方程進行仿真，得到圖1模型。

圖1 濾波白噪聲Simulink模型

在band-limited White noise中生成白噪聲原始數據，需要確定采樣時間間隔。文獻[11～12]均采用固定的simple time，但由于采樣頻率固定，白噪聲產生的數據相同，不同速度下經過Gian增益模塊和Integrated模塊積分變換得到的路面譜幅值不同，圖2為v=10 m/s和v=20 m/s的路面不平度曲線，2組路面在仿真過程中僅僅是速度不同，轉化到空間頻域內不平度的幅值應該保持不變，但2組路面的幅值存在差異，采用固定采樣頻率確定的白噪聲并不符合實際。

圖2 固定采樣時間間隔時濾波白噪聲仿真

王永生[13]對白噪聲2種生成方式White noise和band-limited White noise中的參數取值問題進行了詳細的討論，確定了仿真步長與數據間隔的關系。采用變采樣頻率，采樣頻率與速度成反比時，仿真得到的數據轉換到空域范圍滿足空域不變原則。這里根據最小采樣定理確定的變采樣頻率fs=2nmaxv，仿真步長Δt取值為1/2nmaxv，對A級路面不同速度下路面譜分別進行仿真，得到v=10 m/s和v=20 m/s的路面不平度曲線圖3，可以看到，當速度增大1倍，響應的功率譜幅值一定，幅頻增加1倍，符合道面不平度空間不變規律，進而可以根據位移時間換算關系x=vt，將時間頻域不平度轉化到空間頻域，見圖4。

圖3 變采樣時間間隔時濾波白噪聲仿真

圖4 濾波白噪聲法生成的隨機路面不平度

濾波白噪聲法采用Simulink生成白噪聲，通過函數積分器很方便求解不平度，圖5為Welch周期法計算的隨機路面其功率密度函數與標準路面功率譜的對比，擬合效果良好。

圖5 濾波白噪聲法功率譜密度及其與標準譜對比

3 單輪起落架力學響應模型

3.1 模型的建立

類比于國際平整度指數的四分之一車輛力學響應模型，考慮起落架的緩沖系統受力、輪胎線性阻尼和滑行狀態升力的影響，建立圖6所示單輪起落架動力響應計算模型。

圖6 單輪起落架力學響應模型

圖6中：m2、m1分別為彈性支撐質量(包括機體、機翼及緩沖器的外筒)和非彈性支撐質量(包括起落架支柱及緩沖器、機輪)，kg；Foil為緩沖系統油液阻力；Ffri為緩沖系統摩擦力；Flim為緩沖系統限制力；Fair為緩沖系統空氣彈簧力；L為飛機升力；k1為輪胎剛度，N/m；c1為輪胎阻尼系數；μ為皮碗摩擦系數；q為道面激勵豎向位移，即道面不平度，m。

由達朗貝爾原理可知，飛機滑行時平衡位置處(重力減升力)機體動力學模型如下：

(4)

機輪動載為平衡狀態輪胎受力與不平度激勵輪胎受力的合力，則機輪動載系數D由下式進行計算：

(5)

有關升力影響下的動載系數的研究，張獻民[6]將動載系數修正為不考慮升力時的動載系數與升力影響因子的乘積，物理意義不明顯。也有很多學者[14-16]在采用達朗貝爾原理進行飛機起落架滑跑受力分析時，考慮了系統的重力，實際上飛機的滑跑過程初始條件應為系統的平衡狀態，系統的2組彈簧均處于受壓狀態，彈簧力與重力大小相等方向相反，飛機的簧載質量和非簧載質量均在各自的平衡位置上下振動。但在計算輪胎受力時，需要對這部分受力進行疊加，因此采用式(5)計算的動載系數更為合理。

3.2 相關力的計算

3.2.1 飛機升力

L=CLρairSwinv2/2

(6)

式中：CL為升力系數；ρair為空氣密度；Swin為機翼有效面積；v為空速。

飛機起飛離地速度時，升力等于重力：

(7)

式中：vq為起飛離地速度。

飛機以速度v在光滑的道面勻速滑行時，對道面的作用力F為：

(8)

3.2.2 空氣彈簧力Fair[17]

Fair=[P0(V0/(V0-SAa))γ-Patm]Aa

(9)

式中：P0、V0分別為氣腔的初始壓力和初始體積；Patm為標準大氣壓；S為緩沖行程；Aa為氣腔有效壓氣面積；γ為氣體的多變指數，取值為1.0～1.4。

式(9)為非線性方程，且方程的參數較多，求解過程中很容易出現不收斂的情況，因此，在研究道面不平度時，為簡化計算，轉化為線性微分方程，常假設采用停機下線性剛度進行計算。具體為[18]：

(10)

3.2.3 油液阻尼力Foil[19]

(11)

3.2.4 摩擦力

緩沖器的摩擦力包括皮碗摩擦力和彎曲摩擦力，因為彎曲摩擦力為緩沖支柱彎曲上下支撐點產生的摩擦力，所以這里不考慮水平方向的力。

為簡化計算，僅考慮皮碗摩擦力，如式(12)：

Ffrc=μmFair

(12)

3.2.5 結構限制力

緩沖器伸長或壓縮可能超過限制行程，此時，外筒結構會限制緩沖器進一步變形，從而產生結構限制力，很多模型忽略了結構限制力的影響，對于研究不平度激勵下的位移影響不大，但對于著陸沖擊階段的位移往往不符合實際，因此采用式(13)計算：

(13)

式中：Ks為緩沖器剛度；Smax為緩沖器的極限行程。

簡化后的模型(見圖7)和方程為：

(14)

式中：k2為停機時空氣彈簧當量線性剛度，N/m；c2為油液當量線性阻尼，N/m；sign為符號函數。

(15)

Matlab中Simulink計算功能很強，對于耦合微分方程的求解直觀、高效，本文采用Simulink求解飛機著陸滑跑過程的動力學響應。

圖7 簡化后的單輪起落架力學響應模型

4 仿真計算

4.1 仿真參數的確定

仿真采用的相關參數取值見表2[20]。建立的Simulink模型見圖8。

表2 仿真參數及其取值

圖8 求解機輪動載的Simulink模型

由于道面不平度在濾波白噪聲法中看作為均值為0的高斯隨機過程，由道面不平度產生的機輪對道面的激勵荷載同樣為高斯隨機場。本文采用3σ原則(99.73%的保證率)確定最大動載系數，采用式(16)進行計算：

(16)

式中：N為仿真樣本數量；Di為第i個樣本點的動載系數。

這里除了考慮這些因素，還根據IRI的計算公式，定義了起落架彈性支撐質量m2和非彈性支撐質量m1的累積相對位移，用Za表示，即：

(17)

4.2 仿真初始條件確定

該仿真過程有6個初始條件，分別為道面功率譜密度、滑跑速度、簧載質量和非簧載質量的下落速度以及輪胎和緩沖器的初始豎向位移。上述模型并未考慮起落架著陸過程的起轉回彈，僅能模擬飛機的滑跑狀態，簧載質量和非簧載質量的下落速度均取值為0。簧載質量和非簧載質量均在各自的平衡位置上下振動，輪胎和緩沖器初始豎向位移為0時，在生產隨機路面時功率譜密度取值為0的條件下，動載系數保持為常數。但實際彈簧會有壓縮，壓縮量等于系統縮減重力與彈簧剛度的比值。

4.3 仿真結果分析

將初始條件代入模型中，對不同速度和道面平整度耦合作用下的動載系數進行仿真，對動載系數均值、動載系數標準差、3σ準則確定的動載系數最大值和緩沖器的累計相對位移進行統計分析。

仿真計算發現，不同道面平整度，相同滑跑速度下的動載系數的均值保持不變，分析圖9可知，動載系數均值隨著滑跑速度的增加逐漸減小，二次關系擬合良好，方程可由式(18)表達：

(18)

式(18)表明，道面不平度引起的動載系數均值與道面不平度無關，僅與特定滑跑速度下的重力與升力的差值有關。這是因為一方面在進行起落架模型簡化時假設系統為線彈性體，另一方面采用的濾波白噪聲法仿真產生的道面不平度激勵均假設為均值為0的高斯隨機過程，動載系數均值等于既定速度下的系統平衡狀態動載系數。

圖9 動載系數均值隨滑跑速度的變化規律

由圖10可知，滑跑速度相同時，道面不平度越差，飛機顛簸越劇烈，動載標準差越大。圖11顯示，同一道面不平度下，當滑跑速度較小時，動載系數隨滑跑速度的增大而持續增大，當滑跑速度增大到60～70 m/s時，各平整度下的平均動載系數減小2.61%，原因可能與系統的固有振動頻率有關。由于模型假設起落架在平衡位置上下振動，飛機剛好達到起飛離地速度仍與道面接觸時，還會在不平度激勵下產生較大波動脫離跑道，動載系數突變為0。

圖10 動載系數標準差隨IRI等級變化規律

圖11 動載系數標準差隨滑跑速度變化規律

由圖12可知，當滑跑速度為0時，飛機為靜止狀態，動載系數恒等于1。滑跑速度一定時，動載系數最大值隨著道面不平度惡化而增大，近似成線性關系。圖13顯示，當道面平整度一定時，隨著滑跑速度的增大，動載系數最大值先增大后減小。飛機在低速滑跑狀態下，對道面的動載作用更大，動載系數的最大值均出現在滑跑速度20 m/s以下，IRI=6時的動載系數最大值可達1.23，比靜載狀態下的機輪荷載提高23%。隨著滑跑速度的進一步增大，由于飛機升力的作用，動載系數迅速減小。

圖12 動載系數最大值隨IRI等級變化規律

圖13 動載系數最大值隨滑跑速度變化規律

為了確定不同平整度下飛機滑跑的最不利滑跑速度和對應的動載系數最大值，要對上述仿真數據進行擬合。文獻[6]將起落架和輪胎均假設為帶一次阻尼的線性彈簧，通過頻響函數推導了隨機動載標準差的計算公式：

(19)

H1(ω)=[ω4-iα2(1+h)ω3-β2(1+h)ω2]/

[ω4-i(α1+hα2+α2)ω3-(α1α2+β1+β2+hβ2)ω2+

i(α1β2+α2β1)ω+β1β2]

(20)

式中：h=m2/m1；α1=c2/m1；α2=c2/m2；β1=k1/m1；β2=k2/m2。

根據3σ準則，計算的最大動載系數為：

(21)

化簡可得：

(22)

(23)

式中：c0=10-3s0.5m-0.5

將不同等級的IRI代入可以得到對應等級的DCL表達式，分別對式(23)進行求導，可得各道面平整度下飛機最不利滑跑速度和對應的最大動載系數，見表3。

表3 飛機最不利滑跑速度和最大動載系數

與最大動載系數的變化規律相似，在一定的速度下，累計相對位移與道面IRI等級近似成線性關系，道面不平度越差，累計相對位移越大。通過比較折線的斜率，可以看出隨著速度的增大，飛機升力增大，累計相對位移隨不平度的變化逐漸減小(見圖14、圖15)。

圖14 累計相對位移隨IRI等級的變化規律

圖15 累計相對位移隨滑跑速度的變化規律

當道面不平度等級保持不變時，累計相對位移隨速度的增大先增大后減小。低速時，累計相對位移的變化較大，速度在10 m/s附近時普遍達到最大值。此后隨速度的增大，累計相對位移的變化逐漸減小，如飛機滑跑速度從10 m/s增大到70 m/s，按照10 m/s的間隔，各不平度下的平均累計相對位移分別減小23.64%、17.53%、13.45%、10.56%和11.60%，究其原因是IRI指標采用的是車輛模型，滑跑過程中未考慮升力的影響。

四分之一車輛模型是在標準車身懸架剛度、輪胎剛度和懸架阻尼條件下，測定80 km/h速度懸掛系統的累計相對位移。此處為起落架模型，也可以計算不同IRI對應的實際累計相對位移，見表4。

表4 各道面平整度下飛機起落架累計相對位移

相對于四分之一標準車輛模型，單輪起落架模型累計相對位移較大。隨著平整度等級均勻增大，從1增大到6，以1為間隔，累計相對位移分別增大70.41%、35.41%、23.33%、17.46%和14.51%，增大速率逐漸減小，原因是本文的起落架模型采用的是非線性油液阻尼，顛簸越劇烈，緩沖系統受到的阻尼力呈二次關系增大，增速越來越大，導致累計相對位移逐漸減小。

5 結論

考慮起落架緩沖器系統、輪胎受力及不同速度下升力的影響，本文建立了單輪起落架力學響應計算模型，研究了道面不平度和速度耦合作用下機輪動載的統計變化規律。研究發現：道面不平度引起的動載系數均值與道面不平度無關，僅與特定滑跑速度下的重力與升力的差值有關；滑跑速度相同時，道面不平度越差，動載系數標準差越大，同一道面不平度下，滑跑速度較小時，動載系數標準差隨滑跑速度的增大而持續增大，但存在與系統固有頻率相關的最敏感速度；根據3σ準則確定的動載系數最大值，隨滑跑速度的增大先增大后減小。飛機在低速滑跑時對道面的動載作用更大，計算了對應平整度等級的最不利滑跑速度；由于緩沖系統二次油液阻尼的影響，隨著IRI數值均勻增大，起落架模型的累計相對位移增速越來越小。