基于格子Boltzmann方法的起落架阻尼油孔參數研究

2023-01-31 07:47:24甘盛勇魏小輝房興波

振動與沖擊 2023年2期

甘盛勇，魏小輝，房興波

(1.南京航空航天大學航空學院，南京 210016； 2.南京航空航天大學機械結構力學及控制國家重點實驗室, 南京 210016； 3.南京航空航天大學飛行器先進設計技術國防重點學科實驗室, 南京 210016)

油氣式緩沖器有著緩沖效率高、結構穩定性好的優點，被廣泛應用在各種型式的飛機起落架中。油孔阻尼是緩沖器吸收著陸沖擊能量和減小地面滑跑垂向擾動的關鍵部件，在緩沖過程中將飛機降落的機械能轉化為內能，油孔阻尼產生的載荷則直接影響飛機垂向過載。為了滿足機載設備工作環境和乘客舒適性的設計要求，需要對油孔阻尼性能進行準確分析和設計。

目前起落架緩沖動力學分析基于二自由度質量模型，緩沖器運動特性通過兩質量運動差異得出。緩沖器運動特性、油孔結構和初始充填參數共同決定緩沖載荷的大小。其中緩沖器壓縮速度特性和油孔截面積不僅直接被用于對阻尼力的計算公式中，還對流動縮流系數的大小產生影響，進一步影響阻尼力變化趨勢。目前對于縮流系數的確定方法可以分為工程估算法、經驗公式法和計算流體動力學(computational fluid dynamics ,CFD)分析法。工程估算法計算油孔阻尼力時，在整個緩沖過程中將縮流系數確定為常數[1-2]，對于直角入口圓柱孔建議取值區間為0.632～0.707[3]。經驗公式法是對特定結構油孔在有限參數區間內進行試驗，根據試驗結果，得到縮流系數關于部分關鍵參數的計算公式，利用該公式計算試驗參數區間內的縮流系數[4-6]。CFD分析法基于數值計算方法，利用特定參數下的試驗數據驗證計算模型結果的準確性后，對一般結構參數的油孔進行流動特性分析[7-8]。工程估算法在設計完成后利用起落架落震試驗結果對流孔結構進行修正，需要較多的時間成本和經濟成本。經驗公式法得到的公式只適用于試驗參數區間內，其不可拓展性限制了該方法的應用范圍。CFD分析法具備時間和經濟成本上的優勢，在驗證計算準確性的前提下，對于阻尼油孔分析具有一定有效性。

格子Boltzmann方法(lattice Boltzmann method,LBM)是一種基于粒子的CFD方法，通過模擬粒子的傳播和碰撞過程來模擬流體的宏觀運動，在處理復雜邊界、介觀-微觀相互作用和算法并行化方面具備優勢[9]，被廣泛應用在醫學[10]、化學[11]和工程[12]中。基于LBM的大渦模擬(large eddy simulation,LES)分析高雷諾數流動的精度與傳統的有限體積法(finite volume method，FVM)的LES相當，且LBM的計算效率更高[13]。浸入邊界-LBM在計算流固耦合問題的有效性[14]可以應對更為復雜及可變的油孔結構阻尼特性分析，例如對動結構邊界的處理方法[15]。

現有文獻中關于緩沖器油孔阻尼的CFD分析更多關注油孔幾何參數在起落架落震性能中的影響[16-18]，以及油孔流動靜壓力和流速分布情況[19]。油孔幾何參數對落震性能的影響規律會隨著起落架其他參數變動而失去代表性，油孔流動靜壓力和流速分布也不能直接指導油孔結構設計。為了得到直角入口圓柱孔縮流系數的一般計算公式，基于格子Boltzmann方法對油孔在不同幾何參數和工作環境參數下的阻尼性能進行仿真計算，在單參數分析結果的基礎上，利用多元線性回歸方法得到縮流系數在設計參數區間內的擬合方程。本研究對不同參數的油氣式起落架定截面油孔阻尼設計具有一定指導意義。

1 油孔結構及流動特性分析

1.1 油孔結構

油氣式緩沖器內部被油孔分割成兩個腔體，其中一側腔體充填滿油液，稱為油腔，油液為不可壓縮流體，另一側腔體充填部分液體和部分氣體，稱為氣腔。飛機降落或滑跑過程中，緩沖器被壓縮后，油腔內部靜壓力首先增大，油液在靜壓力作用下，通過油孔流入氣腔減小氣體體積，氣體壓力增加為起落架提供氣體彈簧載荷。緩沖器內部結構如圖1所示。

圖1 緩沖器及流孔流動示意圖Fig.1 Diagram of shock absorber and oil hole flow

1.2 油孔理論特性分析

油液流動截面從壓油面積突然減小至流孔面積，流動截面收縮產生的局部阻力是吸能的主要因素。油液流動截面收縮后，在油孔入口后一定距離處產生一個射流直徑最小的收縮截面c-c，收縮截面位置根據油孔長度可能分布在流孔內部或者流孔出口后端，收縮截面位置對油孔流動特性有一定影響。因為來流速度遠小于流孔內部截面平均流速，忽略來流速度的動能，在不考慮油孔產生的局部阻力和沿程阻力條件下，根據伯努利方程，油孔出口理論平均流速應該為

(1)

式中：vth為油孔出口理論平均流速；ΔP為油孔兩側靜壓力差；ρ為油液密度。根據縮流系數的物理意義，即油孔流動實際質量流量和理論質量流量的比值，則縮流系數的計算形式可以表述為

(2)

(3)

2 格子Boltzmann方法及仿真設置

2.1 格子Boltzmann方法

通過對Boltzmann-BGK方程在速度、時間和空間上的離散得到格子Boltzmann-BGK方程[20]為

(4)

(5)

式中：wi為權系數；u為流體宏觀平均速度；cs為格子特征聲速。平衡態分布函數的計算取決于流體局部密度ρ和速度u，這兩個量可以通過局部粒子分布函數統計值計算得來，而流體宏觀速度則可以通過動量密度和質量密度的比值得到，具體表達形式為

(6)

格子Boltzmann-BGK方程描述粒子運動可以分為兩個過程：碰撞項和傳播項，表述形式分別為

(7)

圖2 格子Boltzmann方法計算過程Fig.2 Lattice Boltzmann method calculation process

2.2 仿真設置

緩沖器定截面油孔為直角入口圓柱孔，影響流孔阻尼性能的主要參數有：油孔長度和直徑、流孔出口處靜壓力。采用控制變量法，在其他參數保持不變的情況下，改變研究變量的大小，設置仿真工況進行計算，各參數在仿真中的基礎值如表1所示。

表1 油孔參數基礎值Tab.1 Oil hole parameters base value

仿真中邊界條件的設定盡量模擬真實的緩沖過程。邊界條件設置如表2所示。

表2 仿真邊界條件設置Tab.2 Simulation boundary condition setting

根據在參數為基礎值工況下進行的計算區域格子尺寸收斂性分析結果，將油孔處格子大小設置為0.18 mm，格子尺寸沿遠離油孔方向依次加倍，共四個尺寸層級，格子設置如圖3所示。流體介質為15號航空液壓油，其在溫度為40 ℃時的密度為837 kg/m3，運動黏度為13.85 mm2/s。

圖3 多層級流體格子設置Fig.3 Multi-level fluid lattice setup

3 計算結果與分析

3.1 模型驗證

利用Cheng等[21]總結的縮流系數經驗公式、Jiang等發表的直角圓柱孔流動試驗數據和基于雷諾平均方程(Reynolds-averaged navier-stokes，RANS)的有限體積法數值結果來分析LBM計算模型的準確性。模型中流孔幾何參數、邊界條件和流體屬性與Jiang等研究中試驗參數設置保持一致。圖4中對比了試驗數據、RANS數值結果、經驗公式和格子Boltzmann方法計算模型結果中縮流系數隨靜壓力差變化規律。

圖4 縮流系數隨靜壓力差變化曲線Fig.4 Variation curve of discharge coefficient with static pressure

根據圖4可得，縮流系數在壓差計算區間內呈逐漸增大規律，增大速度逐漸放緩。格子Boltzmann方法計算結果與經驗公式結果的相對誤差較小，最大誤差為1.51%，與試驗結果的相對誤差較大，最大誤差為5.48%。RANS數值結果與試驗結果的最大相對誤差為6.56%。通過四組縮流系數的對比，驗證了格子Boltzmann方法在分析直角圓柱孔流動特性中的準確性，可以用于對起落架油孔阻尼性能分析中。

3.2 流孔長度

通過對各長度流孔的流動特性計算，得出阻尼力和縮流系數變化趨勢。考慮一般流孔厚度，分別選取長度為5 mm，10 mm，20 mm，30 mm和40 mm的流孔進行計算，其余參數設定為表 1中的基礎值。得到阻尼力與壓縮速度曲線和縮流系數與壓縮速度曲線如圖5所示。

圖5 各流孔長度阻尼性能隨壓縮速度變化曲線Fig.5 Variation curve of damping performance with compression velocity for each length of oil hole

從圖5(a)可以看出，流孔長度從5～30 mm變化時，阻尼力依次減小，流孔長度從30～40 mm變化時，阻尼力增大。當流孔長度為5 mm時，流孔長徑比為1.3，屬于薄壁孔的范疇，流動經收縮后擴散沒有附壁，流動特性只與壓力差有關。而隨著流孔長度的增加，流動收縮后擴散開始附壁，收縮界面處流速最大，靜壓力則最小，形成的負壓空間有助于液體通過流孔，所以阻尼力會依次減小。流動從收縮截面擴散后附壁，沿程阻力增加，又增大阻尼力，解釋了流孔長度為40 mm時阻尼力略有增大的情況。圖5(b)中縮流系數在壓縮速度為0.1～1.5 m/s時，隨著壓縮速度的增大而增加，壓縮速度在1.5～2.0 m/s時，縮流系數隨著壓縮速度的增加而減小，縮流系數隨流孔長度增加呈先增大后減小趨勢，內在原因和阻尼力變化的因素一樣。

圖6為壓縮速度為2.0 m/s時的各流孔長度下水平流速分布情況，因為流動質量流量一致，所以各長度流速大小基本一致。從圖6(a)中可以看出，流孔長度為5 mm時，液體收縮截面位于流孔內部，但液體收縮后擴散沒有附壁，總體依然呈薄壁孔流動特性。圖6(b)流孔長度增加到10 mm，流動擴散后已經能夠附壁，當流孔長度繼續增大，流動附壁長度加長，增加了沿程阻力，并且形成的收縮截面直徑變小，因此兩側形成的回流空間也更大。在圖6(d)中能夠觀察到明顯的液體收縮和擴散的流線變化，而在圖6(a)中流動收縮后幾乎沿水平方向流動，沒有發生明顯擴散。

圖6 各流孔長度液體水平流速分布Fig.6 Horizontal flow velocity distribution for each length of oil hole

圖7對比了2 m/s壓縮速度下流動附壁前后流孔內部靜壓力分布變化。圖7(a)中是流孔長度為5 mm時的靜壓力分布，由于收縮截面在流孔內部，流孔內形成了負壓空間，集中在流孔入口處兩側區域。圖7(b)是流孔長度為10 mm時的靜壓力分布，此結構下流動已經附壁，形成的負壓空間在流孔入口處的連續區域內。因此流動附壁前后流孔內部都出現了一定區域的負壓空間，這是收縮截面處流速最大導致的，然而流動附壁后形成的負壓區域分布在流孔入口后的連續截面上，形成的區域比附壁前更大。較大和較連續的負壓區域加強了液體的通過能力，進一步驗證了流孔長度對縮流系數的影響趨勢。

圖7 流動附壁前后靜壓力分布Fig.7 Static pressure distribution before and after flow attached wall

3.3 流孔直徑

流孔直徑是唯一直接體現在阻尼力計算公式中的流孔幾何參數，所以也是唯一在縮流系數工程估算方法中使用到的幾何參數，因而流孔直徑對阻尼的影響分為直接作用和間接作用，間接作用指的是流孔直徑對縮流系數的影響。考慮緩沖器流孔直徑尺寸設置和對阻尼力大小的變動區間，分別選取流孔直徑為3.4 mm，3.6 mm，3.8 mm，4.0 mm和4.2 mm的流孔進行計算。得到阻尼力與壓縮速度曲線和縮流系數與壓縮速度曲線如圖8所示。

圖8 各流孔直徑阻尼性能隨壓縮速度變化曲線Fig.8 Variation curve of damping performance with compression velocity for each diameter of oil hole

圖8(a)對比了各流孔直徑阻尼力計算結果，阻尼力隨著流孔直徑的增大而減小。流孔直徑減小0.2 mm，阻尼力最大增加了1.96 kN，并且隨著流孔直徑的減小，相同直徑差產生的阻尼力差值逐漸增大。圖8(b)中對比了各流孔直徑縮流系數計算結果，縮流系數隨流速增加先增大后減小，隨流孔直徑變化在三個區域內呈現不同的規律。在圖中左側區域，縮流系數隨流孔直徑增大總體呈增大趨勢，說明壓縮速度較低時，孔內流動速度相應較小，流動沒有發生分離，所以流動的有效面積隨著流孔的增大相應的增加。中間區域縮流系數隨流孔直徑增大呈減小趨勢。在圖中右側區域，壓縮速度繼續增加，縮流系數隨流孔直徑增大呈增大趨勢，同左側區域變化趨勢一致。

壓縮速度增加使得流孔內部流速相應增大，流孔內部流動由層流逐漸過渡為湍流，轉變發生的臨界壓縮速度與孔內流速相關。孔內截面理論平均流速隨緩沖器壓縮速度變化趨勢如圖9所示，孔內流速隨著壓縮速度呈線性變化，直徑較小的油孔孔內流速更大。孔內流動在更低的壓縮速度下就轉變為湍流，縮流系數也隨著流動特性產生變化。故在圖8(b)中部區域小直徑流孔縮流系數首先增大，在右側區域縮流系數首先減小。

圖9 理論平均流速分布Fig.9 Theoretical average flow velocity distribution

圖10對比了各流孔直徑下縮流系數隨流孔截面理論平均速度變化趨勢，除了在低流速區域中的不同，各流孔的變化趨勢基本一致。當理論流速低于40 m/s時，縮流系數隨著流孔直徑增大而減小。當理論流速大于40 m/s時，各流孔縮流系數曲線基本重合，說明在孔內流速較大時，流孔直徑對縮流系數的影響僅限于改變了同壓縮速度下的孔內流速，而不同直徑在有限區間內的變化對縮流系數幾乎沒有影響。

圖10 縮流系數隨理論平均流速變化曲線Fig.10 Variation curve of discharge coefficient with theoretical average flow velocity

3.4 后端壓力

工程估算方法中沒有考慮后端壓力，且緩沖器行程變化過程中，氣腔壓力持續變化，對流孔流動特性產生影響。考慮起落架氣腔初始充填壓力和工作最大壓力，分別選取后端壓力0.5 MPa，1.5 MPa，2.5 MPa和3.5 MPa進行計算，得到阻尼力與壓縮速度曲線和縮流系數與壓縮速度曲線如圖11所示。

圖11 各后端壓力阻尼性能隨壓縮速度變化曲線Fig.11 Variation curve of damping performance with compression velocity for each backpressure

不同后端壓力下阻尼力在低壓縮速度時的差異很小，在較高壓縮速度時，開始出現明顯的差異，阻尼力隨著后端壓力增大而減小。縮流系數隨壓縮速度增大先增大后減小，且在各個壓縮速度下隨后端壓力增大均呈增加趨勢，隨壓縮速度增加縮流系數變化規律基本一致。在低壓縮速度時，各后端壓力下縮流系數間的差距小于在較高壓縮速度下的間距。后端壓力主要影響流動分離的臨界條件，后端壓力越大流動越不容易發生分離，則流動特性更有利于流體通過流孔，導致更小的阻尼力和更大的縮流系數。流動分離產生的情況隨著流速的不斷增加而變得嚴重，所以在高流速情況下，由于后端壓力增加產生的改善流動特性的效果則更明顯，體現在縮流系數增長的更多。

4 回歸方程及驗證

上述通過對LBM計算結果的分析，得到了各參數對流孔縮流系數的影響結果，然而單參數的敏感性分析結果不能直接用于緩沖器流孔設計，需要將各參數對阻尼性能的影響表示成統一的形式，用于阻尼力計算中。

選擇流孔長度、流孔直徑、后端壓力和壓縮速度作為設計變量，選擇流孔處阻尼力作為目標值，考慮參數個數和研究水平個數基本處于較低的數值，使用全因子試驗方法，利用LBM計算樣本點處阻尼力數值。各設計變量取值范圍根據上述分析中使用的結構參數基礎值進行確定，選取的各參數設計區間如表3所示。

表3 阻尼力變量設計區間Tab.3 Interval of design variables for damping force

樣本點數據計算在搭載了兩個CPU的計算機上執行，CPU型號為Intel Xeon E5-2650 v4，每個CPU含有12核心24線程，378個樣本點數據的計算時間總共為68 h。基于最小二乘法擬合得到了阻尼力回歸方程，均方根誤差為41.3，R2為0.997,說明回歸方程的擬合精度可靠。圖12中樣本點數據分布與45°對角線相對位置也說明了回歸方程的準確度。上述對比結果說明了回歸方程在參數區間內可以被用于阻尼力的計算，結合數值優化方法，可以得到使得起落架緩沖性能最優的油孔結構參數。