甘盛勇， 魏小輝 ， 房興波

(1.南京航空航天大學(xué) 航空學(xué)院，南京 210016； 2.南京航空航天大學(xué) 機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點(diǎn)實(shí)驗(yàn)室, 南京 210016； 3.南京航空航天大學(xué) 飛行器先進(jìn)設(shè)計技術(shù)國防重點(diǎn)學(xué)科實(shí)驗(yàn)室, 南京 210016)

油氣式緩沖器有著緩沖效率高、結(jié)構(gòu)穩(wěn)定性好的優(yōu)點(diǎn)，被廣泛應(yīng)用在各種型式的飛機(jī)起落架中。油孔阻尼是緩沖器吸收著陸沖擊能量和減小地面滑跑垂向擾動的關(guān)鍵部件，在緩沖過程中將飛機(jī)降落的機(jī)械能轉(zhuǎn)化為內(nèi)能，油孔阻尼產(chǎn)生的載荷則直接影響飛機(jī)垂向過載。為了滿足機(jī)載設(shè)備工作環(huán)境和乘客舒適性的設(shè)計要求，需要對油孔阻尼性能進(jìn)行準(zhǔn)確分析和設(shè)計。

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

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

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

1 油孔結(jié)構(gòu)及流動特性分析

1.1 油孔結(jié)構(gòu)

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

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

1.2 油孔理論特性分析

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

(1)

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

(2)

(3)

2 格子Boltzmann方法及仿真設(shè)置

2.1 格子Boltzmann方法

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

(4)

(5)

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

(6)

格子Boltzmann-BGK方程描述粒子運(yùn)動可以分為兩個過程：碰撞項(xiàng)和傳播項(xiàng)，表述形式分別為

(7)

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

2.2 仿真設(shè)置

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

表1 油孔參數(shù)基礎(chǔ)值Tab.1 Oil hole parameters base value

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

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

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

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

3 計算結(jié)果與分析

3.1 模型驗(yàn)證

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

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

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

3.2 流孔長度

通過對各長度流孔的流動特性計算，得出阻尼力和縮流系數(shù)變化趨勢。考慮一般流孔厚度，分別選取長度為5 mm，10 mm，20 mm，30 mm和40 mm的流孔進(jìn)行計算，其余參數(shù)設(shè)定為表 1中的基礎(chǔ)值。得到阻尼力與壓縮速度曲線和縮流系數(shù)與壓縮速度曲線如圖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變化時，阻尼力增大。當(dāng)流孔長度為5 mm時，流孔長徑比為1.3，屬于薄壁孔的范疇，流動經(jīng)收縮后擴(kuò)散沒有附壁，流動特性只與壓力差有關(guān)。而隨著流孔長度的增加，流動收縮后擴(kuò)散開始附壁，收縮界面處流速最大，靜壓力則最小，形成的負(fù)壓空間有助于液體通過流孔，所以阻尼力會依次減小。流動從收縮截面擴(kuò)散后附壁，沿程阻力增加，又增大阻尼力，解釋了流孔長度為40 mm時阻尼力略有增大的情況。圖5(b)中縮流系數(shù)在壓縮速度為0.1～1.5 m/s時，隨著壓縮速度的增大而增加，壓縮速度在1.5～2.0 m/s時，縮流系數(shù)隨著壓縮速度的增加而減小，縮流系數(shù)隨流孔長度增加呈先增大后減小趨勢，內(nèi)在原因和阻尼力變化的因素一樣。

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

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

圖7對比了2 m/s壓縮速度下流動附壁前后流孔內(nèi)部靜壓力分布變化。圖7(a)中是流孔長度為5 mm時的靜壓力分布，由于收縮截面在流孔內(nèi)部，流孔內(nèi)形成了負(fù)壓空間，集中在流孔入口處兩側(cè)區(qū)域。圖7(b)是流孔長度為10 mm時的靜壓力分布，此結(jié)構(gòu)下流動已經(jīng)附壁，形成的負(fù)壓空間在流孔入口處的連續(xù)區(qū)域內(nèi)。因此流動附壁前后流孔內(nèi)部都出現(xiàn)了一定區(qū)域的負(fù)壓空間，這是收縮截面處流速最大導(dǎo)致的，然而流動附壁后形成的負(fù)壓區(qū)域分布在流孔入口后的連續(xù)截面上，形成的區(qū)域比附壁前更大。較大和較連續(xù)的負(fù)壓區(qū)域加強(qiáng)了液體的通過能力，進(jìn)一步驗(yàn)證了流孔長度對縮流系數(shù)的影響趨勢。

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

3.3 流孔直徑

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

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

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

壓縮速度增加使得流孔內(nèi)部流速相應(yīng)增大，流孔內(nèi)部流動由層流逐漸過渡為湍流，轉(zhuǎn)變發(fā)生的臨界壓縮速度與孔內(nèi)流速相關(guān)?？變?nèi)截面理論平均流速隨緩沖器壓縮速度變化趨勢如圖9所示，孔內(nèi)流速隨著壓縮速度呈線性變化，直徑較小的油孔孔內(nèi)流速更大?？變?nèi)流動在更低的壓縮速度下就轉(zhuǎn)變?yōu)橥牧?，縮流系數(shù)也隨著流動特性產(chǎn)生變化。故在圖8(b)中部區(qū)域小直徑流孔縮流系數(shù)首先增大，在右側(cè)區(qū)域縮流系數(shù)首先減小。

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

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

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

3.4 后端壓力

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

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

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

4 回歸方程及驗(yàn)證

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

選擇流孔長度、流孔直徑、后端壓力和壓縮速度作為設(shè)計變量，選擇流孔處阻尼力作為目標(biāo)值，考慮參數(shù)個數(shù)和研究水平個數(shù)基本處于較低的數(shù)值，使用全因子試驗(yàn)方法，利用LBM計算樣本點(diǎn)處阻尼力數(shù)值。各設(shè)計變量取值范圍根據(jù)上述分析中使用的結(jié)構(gòu)參數(shù)基礎(chǔ)值進(jìn)行確定，選取的各參數(shù)設(shè)計區(qū)間如表3所示。

表3 阻尼力變量設(shè)計區(qū)間Tab.3 Interval of design variables for damping force

樣本點(diǎn)數(shù)據(jù)計算在搭載了兩個CPU的計算機(jī)上執(zhí)行，CPU型號為Intel Xeon E5-2650 v4，每個CPU含有12核心24線程，378個樣本點(diǎn)數(shù)據(jù)的計算時間總共為68 h?；谧钚《朔〝M合得到了阻尼力回歸方程，均方根誤差為41.3，R2為0.997,說明回歸方程的擬合精度可靠。圖12中樣本點(diǎn)數(shù)據(jù)分布與45°對角線相對位置也說明了回歸方程的準(zhǔn)確度。上述對比結(jié)果說明了回歸方程在參數(shù)區(qū)間內(nèi)可以被用于阻尼力的計算，結(jié)合數(shù)值優(yōu)化方法，可以得到使得起落架緩沖性能最優(yōu)的油孔結(jié)構(gòu)參數(shù)。

圖12 樣本點(diǎn)數(shù)據(jù)分布Fig.12 Sample points data distribution

5 結(jié) 論

主要結(jié)論如下：

(1)格子Boltzmann方法計算得到的阻尼力與試驗(yàn)值最大相對誤差不超過5.48%，說明格子Boltzmann方法在分析直角圓柱孔流動特性中的準(zhǔn)確性，可以用于對起落架油孔阻尼性能的分析中。

(2)油孔長度會改變流動最小收縮截面與油孔的相對位置和流動沿程阻力，從而對阻尼性能產(chǎn)生影響；油孔直徑對阻尼力影響最大，對縮流系數(shù)的影響僅限于改變了同壓縮速度下的孔內(nèi)流速，而油孔直徑在有限區(qū)間內(nèi)的變化對流動特性幾乎沒有影響；后端壓力增大會減緩油孔內(nèi)部流動分離，增大各壓縮速度下的縮流系數(shù)。

(3)阻尼力關(guān)于油孔長度、油孔直徑、后端壓力和質(zhì)量流量的回歸方程具有較高的擬合度，能夠用于參數(shù)區(qū)間內(nèi)的阻尼力計算，說明格子Boltzmann方法對阻尼性能的分析能夠被應(yīng)用在一般性的起落架緩沖器油孔設(shè)計中。