張祺祥 ，韋 坤 ，王泓林 ，陸鳳霞

（1.中國航空發(fā)動機集團有限公司，北京 100097；2.南京航空航天大學(xué)直升機傳動技術(shù)重點實驗室，南京 210016）

0 引言

直升機減速器內(nèi)軸承發(fā)生異常磨損時會產(chǎn)生金屬屑，通過高速旋轉(zhuǎn)的齒輪攪動油氣使其向底部流動，若其能被機匣底部的信號器吸附及報警，可實現(xiàn)減速器故障監(jiān)控，因此亟需研究金屬屑在中減內(nèi)部固-液-氣的運動狀態(tài)及附加力對其的影響規(guī)律。

Ran 等[1]分析了粒子的邊界層特征和對粒子的作用力，應(yīng)用拉格朗日方法提出了氣固旋轉(zhuǎn)流中粒子運動的3 維數(shù)學(xué)模型，表明在旋轉(zhuǎn)邊界層中離心力和薩夫曼升力對顆粒與氣固旋轉(zhuǎn)流的分離過程中起著重要作用；Mithun 等[2]提出了一種三相流完全可壓縮模型和浸入邊界模型，用于預(yù)測不可凝氣體（Non-Condensable Gas,NCG）存在下2 維齒輪泵發(fā)生的空化現(xiàn)象；Zhang 等[3]提出了一種多尺度的離散相模型（Discrete Particle Mode , DPM）和直接數(shù)值模擬法（Direct Numerical Simulation,DNS）耦合方法DPM-DNS，表明采用DPM-DNS 計算得到的合成相力分布和濃度分布與采用得到的DNS 基本相同，計算速度較DNS 高30 倍以上，但遠(yuǎn)低于DPM；喻澤雄[4]提出了由高斯光誘導(dǎo)的溫度梯度導(dǎo)致的表面張力系數(shù)梯度，研究了表面張力梯度對微粒的驅(qū)動特性；賴科等[5]以毫米級球形顆粒為例，獲得顆粒在流場中的運動軌跡，并分析了顆粒在旋轉(zhuǎn)流場中的受力情況，表明離心力和科氏力是顆粒在旋轉(zhuǎn)流場中運動的主要影響因素；夏鋮[6]基于計算流體力學(xué)-離散元法（Computational Fluid Dynamics-Discrete Element Method, CFD-DEM）耦合方法，選取歐拉-拉格朗日法計算液固兩相流運動，對比分析了不同形態(tài)和體積分?jǐn)?shù)的顆粒在2 級混流泵內(nèi)部的運動規(guī)律來得到顆粒的流動特性；Drew 等[7]推導(dǎo)了虛擬質(zhì)量加速度的一般形式，并通過試驗確定了虛質(zhì)量力的計算公式；黃社華等[8-10]針對任意流場內(nèi)稀疏顆粒非線性運動方程，采用積分變換方法對其進行解析，并提出了以P（EC）k多步法差分格式為基礎(chǔ)的數(shù)值計算方法；李振中等[11]應(yīng)用理論和數(shù)值方法計算稀疏固氣兩相湍流中的顆粒所受各種相間力，同時考慮剪切升力和旋轉(zhuǎn)升力，表明除曳力和重力外，在主流和翼展方向上Basset 力較重要，而在壁面法向上拖曳力、剪切升力和Basset 力均重要；Yang 等[12]結(jié)合動網(wǎng)格法（dynamic mesh，DM）和DPM 模型計算顆粒運動軌跡和分布，結(jié)果表明顆粒均呈離心運動趨勢，粒子的平均運動速度與齒輪轉(zhuǎn)速和中心距呈非線性正、負(fù)關(guān)系；陳有斌等[13]研究了旋流場作用下固相顆粒的動力學(xué)特性，并進行粒子圖像測速（Particle Image Velocimetry, PIV）流場測試分析，研究固相顆粒的運動軌跡、停留時間及受力特點。

目前針對顆粒的運動分析多以顆粒的運動特性仿真分析為主，缺少通過顆粒受力的運動方程揭示運動機理研究。本文根據(jù)顆粒受力平衡方程，基于單向數(shù)值計算方法，獲得減速器內(nèi)顆粒的運動軌跡，采用多相流模型（Volume of Fluid, VOF）和DPM 相結(jié)合的CFD 仿真方法，計算雙向耦合作用下的顆粒運動軌跡；對比驗證VOF-DPM 耦合模型的正確性，并利用CFD 仿真對比分析不同附加力對顆粒運動軌跡的影響。

1 VOF-DPM 耦合模型驗證

1.1 顆粒受力方程

根據(jù)牛頓第二定律得到顆粒的載荷平衡方程

式中：mp為顆粒質(zhì)量，kg；u為連續(xù)相的速度，m/s；up為顆粒速度，m/s；ρ為連續(xù)相的密度，kg/m3；ρp為顆粒的密度，kg/m3；F為附加力，N；（u-up）/τr為顆粒受到的流體阻力（曳力），N；τr為顆粒的弛豫時間，s

式中：μ為連續(xù)相的分子粘度，Pa·s；dp為顆粒直徑，m；Cd為阻力系數(shù)；Re為相對雷諾數(shù)，其定義為

式中：a1、a2、a3為常數(shù)，在不同雷諾數(shù)范圍內(nèi)有不同的值[14]

此外，顆粒旋轉(zhuǎn)對其在流體中的運動狀態(tài)也有著一定影響，尤其對于高轉(zhuǎn)動慣量的顆粒，影響更為明顯。為考慮顆粒的旋轉(zhuǎn)，需求解關(guān)于角動量的常微分方程

式中：Ip為顆粒轉(zhuǎn)動慣量，kg·m2；ωp為顆粒角速度，rad/s；ρf為流體密度，kg/m3；Cω為旋轉(zhuǎn)阻力系數(shù)；T為作用在顆粒上的力矩，N·m；|Ω|為相對顆粒-流體角速度，rad/s

對于球形顆粒，其轉(zhuǎn)動慣量Ip為

由于DPM 模型基于拉格朗日坐標(biāo)系，將顆粒視為質(zhì)點，忽略顆粒間的相互作用力、碰撞、摩擦、變形等，通過式（1）從流場中的受力即可得出其加速度，根據(jù)加速計算速度，最后得出顆粒在流場中的運動軌跡。

1.2 編程計算

根據(jù)顆粒的受力方程，編程計算顆粒的運動軌跡及運動速度，為節(jié)省計算成本，未添加任何附加力，則式（1）簡化為

每次迭代時需判斷Cd大小，尤其是在空氣與潤滑油的交界面，此時連續(xù)相的粘度與密度均會發(fā)生突變。此外還需判定當(dāng)前時刻顆粒的空間位置是否撞擊到箱體或齒輪壁面，如果撞擊到則需變換顆粒的速度矢量繼續(xù)求解。計算流程如圖1所示。

圖1 計算流程

1.3 CFD仿真

由于3 維減速器中顆粒運動軌跡復(fù)雜，本文以簡化的2 維減速器進行分析。建立齒輪箱平面如圖2 所示。齒輪箱長434 mm，高300 mm，壓力出口寬20 mm，對其劃分非結(jié)構(gòu)化網(wǎng)格以適應(yīng)復(fù)雜流動狀態(tài)，以壓力出口的質(zhì)量流量為監(jiān)測量，經(jīng)過網(wǎng)格無關(guān)性驗證，網(wǎng)格數(shù)量最終確定為18742，網(wǎng)格質(zhì)量最大偏斜度為0.5812，符合CFD仿真計算標(biāo)準(zhǔn)。

圖2 齒輪箱平面

齒輪模數(shù)、齒數(shù)、轉(zhuǎn)速、浸油深度、顆粒物性參數(shù)等齒輪及滑油的主要參數(shù)設(shè)置見表1。

表1 齒輪及滑油的主要參數(shù)

為更好研究顆粒的運動狀態(tài)，在齒輪嚙合處預(yù)置單個顆粒，流場計算模型初始狀態(tài)如圖3 所示。開啟雙向耦合模型，時間步長設(shè)為5×10-6s，計算0.1 s 即主動輪轉(zhuǎn)動約10圈。

圖3 流場計算模型初始狀態(tài)

1.4 對比驗證

應(yīng)用編程計算方法與CFD 方法分別計算顆粒運動軌跡，如圖4（a）所示。從圖中可見，在初始時刻顆粒從齒輪嚙合處至進入油液中，2 種方法計算的顆粒軌跡極為接近；在顆粒射入油液中后，在浮力和流體曳力的綜合作用下，顆粒軌跡有2 處上浮的位置，而CFD 方法僅有1 處（圖中黑色方框處），而其余位置軌跡誤差很小；待顆粒隨齒輪轉(zhuǎn)動再次飛出油液后，編程計算方法計算的顆粒飛出角度略大，導(dǎo)致其后續(xù)運動軌跡與CFD方法的略有差異。

圖4 編程計算方法與CFD方法計算的顆粒運動狀態(tài)

編程計算方法與CFD 方法計算的顆粒運動速度如圖4（b）所示。從圖中可見，與CFD 仿真分析方法相比，編程計算方法中顆粒剛從齒輪嚙合處飛出時速度較小，其余時刻均較大，導(dǎo)致計算所得運動軌跡與CFD仿真分析方法的略有差異，其原因在于CFD方法考慮了連續(xù)相與離散相的耦合作用。與編程計算方法中的單向耦合相比，CFD方法計算的顆粒載荷略有差別，二者計算的運動軌跡誤差約為5%，速度誤差除剛開始時刻外其余約為3%，充分驗證了VOF-DPM耦合模型的正確性，同時也證明了在本文顆粒質(zhì)量輕、數(shù)量少的參數(shù)特點下，連續(xù)相與離散相的耦合作用對顆粒軌跡的影響較小，可以忽略。

2 附加力對金屬屑運動的影響

顆粒在流場中除受到流體浮力、曳力、顆粒自身重力外，由于流體與顆粒的相對運動和流體自身的粘性，還受到其他的附加表面力作用。根據(jù)顆粒的承載類型，可將作用力分為如下類型：重力與浮力、氣動阻力（流體曳力）、虛質(zhì)量力、壓力梯度力、熱泳力、Brownian力、Saffman升力、Magnus升力等。

本章應(yīng)用CFD 方法仿真計算顆粒在流場中的運動狀態(tài)，并分析不同附加力對金屬屑運動狀態(tài)的影響。

2.1 虛質(zhì)量力

顆粒周圍流體加速運動時作用在顆粒上的附加力為虛質(zhì)量力（Virtual Force）

式中：Cvm為虛質(zhì)量系數(shù)，默認(rèn)值為0.5。

從式（10）中可見，當(dāng)流體的密度遠(yuǎn)小于顆粒的密度時，虛質(zhì)量力數(shù)值會很小，例如本文中空氣和顆粒的密度比僅為0.00014，但是油液與顆粒的密度比達(dá)到了0.11，此時虛質(zhì)量力不可忽略，因此需要考慮虛質(zhì)量力對顆粒軌跡的影響。

有/無附加虛質(zhì)量力的顆粒運動軌跡對比如圖5所示；主動齒輪分別旋轉(zhuǎn)1、10、20、30、40、50 圈時，顆粒運動狀態(tài)和油液體積分布如圖6所示，其中左側(cè)無附加力，右側(cè)添加虛質(zhì)量力，顆粒的運動位置見圖中黃色方框處（下同）。

圖5 有、無附加虛質(zhì)量力的顆粒運動軌跡對比

圖6 有、無附加虛質(zhì)量力的顆粒運動狀態(tài)及油液分布對比

結(jié)合式（10）可見，當(dāng)顆粒在空氣中運動時，顆粒的密度遠(yuǎn)大于空氣密度，此時空氣對顆粒運動軌跡的影響很小，如圖6（a）、（b）所示；但當(dāng)顆粒進入油液之后，連續(xù)相流體的密度突變，此時連續(xù)相對顆粒附加的虛質(zhì)量力不可忽略，從圖6（c）、（d）中可見，主動輪轉(zhuǎn)動10 圈后，顆粒已經(jīng)脫離油液，且由于虛質(zhì)量力的影響，顆粒脫離油液的射出角度不同，導(dǎo)致受到虛質(zhì)量力的顆粒已經(jīng)達(dá)到箱體的上半部分，撞擊壁面反彈的入射角也不同，顆粒軌跡開始發(fā)生較大的變化，但運動趨勢依舊是由齒輪的旋轉(zhuǎn)和油液來主導(dǎo)，最終隨著油液落至油池底部。

2.2 壓力梯度力

流體的壓力梯度作用在顆粒上的附加力為壓力梯度力（Pressure Gradient Force）

壓力梯度力與虛質(zhì)量力相同，取決于連續(xù)相與離散相的密度比，因此同樣需要考慮壓力梯度力對顆粒運動軌跡的影響。

根據(jù)式（11），流體的密度對顆粒的壓力梯度力有很大影響，有/無附加壓力梯度力的顆粒運動軌跡對比如圖7 所示，有/無附加壓力梯度力的顆粒運動狀態(tài)及油液分布對比如圖8 所示。從圖中可見，主動齒輪轉(zhuǎn)動10 圈時，顆粒尚未擺脫油液，且在后續(xù)運動中始終處于油液當(dāng)中；右側(cè)添加壓力梯度力相較于左側(cè)無附加力的顆粒，齒輪轉(zhuǎn)動50 圈內(nèi)并沒有圍繞著齒輪運動，反而是受油液的壓力梯度力影響一直位于箱體右半部分，可見油液中壓力梯度力對顆粒運動狀態(tài)影響極為明顯。

圖7 有、無附加壓力梯度力的顆粒運動軌跡對比

圖8 有、無附加壓力梯度力的顆粒運動狀態(tài)及油液分布對比

2.3 Saffman升力

顆粒在有速度梯度的流場中運動時，由于顆粒兩側(cè)的流速不一樣，會產(chǎn)生一項由低速指向高速方向的升力，這就是Saffman升力

式中：K=2.594；dij為變形張量。

從式（12）中可見，Saffman 升力考慮的是剪切變形，這僅在小的顆粒雷諾數(shù)的流體運動中較為明顯，因此Saffman 升力對亞微尺度的顆粒影響較大。有/無附加Saffman 升力的顆粒運動軌跡對比如圖9 所示，有/無附加Saffman 升力的顆粒運動狀態(tài)及油液分布對比如圖10 所示。受到Saffman 升力的顆粒始終未能擺脫油液，除去左側(cè)無附加力的顆粒撞擊齒輪箱壁面反彈后的運動軌跡外，左右兩側(cè)有無Saffman升力的顆粒運動軌跡并無明顯區(qū)別，因此在本文中可以選擇忽略Saffman升力對顆粒運動軌跡的影響。

圖9 有、無附加Saffman升力的顆粒運動軌跡對比

圖10 有、無附加Saffman升力的顆粒運動狀態(tài)及油液分布對比

2.4 Magnus升力

顆粒在流體中旋轉(zhuǎn)時會產(chǎn)生Magnus 升力，這是由沿粒子表面的壓力差引起的，高雷諾數(shù)情況下，Magnus升力FRL為

式中：Ap為顆粒投影面積，m2；V為流體-顆粒相對速度，m/s；Ω為流體-顆粒相對角速度，rad/s。

對于旋轉(zhuǎn)升力系數(shù)CRL，常用以下2種計算方法。

（1）Oesterle and Bui Dinh[15]。旋轉(zhuǎn)升力系數(shù)依賴于旋轉(zhuǎn)雷諾數(shù)Reω與顆粒雷諾數(shù)Rep，該公式在Rep＜2000時與試驗值吻合較好

（2）Tsuji et al[16]。旋轉(zhuǎn)升力系數(shù)定義為自旋參數(shù)S的函數(shù)，該公式在Rep＜1600時廣泛采用

自旋參數(shù)S定義為

有/無附加Magnus 升力的顆粒運動軌跡、運動狀態(tài)及油液分布對比如圖11、12 所示。從圖中可見，在齒輪轉(zhuǎn)動低于20 圈時，有/無附加Magnus 升力的顆粒運動軌跡較為接近，但當(dāng)顆粒第1 次撞擊到齒輪箱上壁面時（圖12（c）），右側(cè)受到Magnus 升力的顆粒被油液裹挾向后運動，由此開始二者的運動軌跡產(chǎn)生了較大的差異，因此Magnus升力對顆粒的運動影響不可忽略。

圖11 有、無附加Magnus升力的顆粒運動軌跡對比

圖12 有、無附加Magnus升力的顆粒運動狀態(tài)及油液分布對比

2.5 熱泳力

存在溫度梯度的流體中懸浮的小顆粒受到的力與溫度梯度方向相反，該現(xiàn)象稱為熱泳，可以在附加力中加入對粒子的熱泳效應(yīng)

式中：Dt，p為熱泳系數(shù)，可定義為常數(shù)、多項式或用戶自定義函數(shù)，也可采用Talbot[17]提出的形式

式中：Kn=2λ/dp，為Knudsen 數(shù)，λ為流體的平均自由程，m；K=k/kp，k為流體熱導(dǎo)率，W/（m·K）；kp為顆粒熱導(dǎo)率，W/（m·K）；Cs=1.17；Ct=2.18；T為流體溫度，K；μ為流體粘度，Pa·s。

本文中減速器內(nèi)部流場馬赫數(shù)小于0.3，空氣為不可壓縮相，根據(jù)以往的研究成果，可不考慮溫度變化和能量方程，因此忽略熱泳力的影響。

2.6 Brownian力

對于亞微尺度的粒子，需要考慮布朗運動的影響。將布朗力的分量建模為譜強度為Sn，ij的高斯白噪聲過程[18]

式中：δij為Kronecker δ函數(shù)，且有

式中：γ為流體的運動粘度，m2/s；Cc為Cunningham 修正系數(shù)；kB為Boltzmann 常數(shù)。布朗力分量的幅值定義為

式中：ξ為與單位方差無關(guān)的高斯隨機數(shù)，均值為零。

每個時間步長計算布朗力分量的幅值，布朗力只適用于層流模擬。由于減速器內(nèi)齒輪的旋轉(zhuǎn)，箱體內(nèi)時刻有湍渦流產(chǎn)生，不適用Brownian 力，因此可不考慮Brownian力對顆粒運動的影響。

3 結(jié)論

（1）采用編程計算與CFD 仿真對比的方式，驗證了仿真模型的準(zhǔn)確性。與單純的連續(xù)相作用相比，連續(xù)相和離散相耦合作用對金屬屑顆粒運動的影響較小，在2 種計算條件下的軌跡誤差約5%、速度誤差約3%，可忽略不計。

（2）基于CFD 方法對2 維減速器的仿真結(jié)果表明，虛質(zhì)量力、壓力梯度力和Magnus 升力對顆粒運動軌跡的影響較大；Saffman 升力對尺寸較大的顆粒運動軌跡影響較小，但對亞微尺寸顆粒運動有影響；熱泳力與布朗力不適用于本文研究對象。在實際分析時應(yīng)當(dāng)結(jié)合實際情況施加相應(yīng)的附加力。

（3）通過2 維減速器模型揭示了顆粒在減速器內(nèi)部固、氣、液多相耦合作用下的運動機理，可為預(yù)測減速器內(nèi)部金屬屑運動狀態(tài)及運動軌跡提供技術(shù)方法。

（4）本文的研究方法方法可以拓展到3 維減速器仿真，以實現(xiàn)更準(zhǔn)確的預(yù)測減速器內(nèi)部金屬屑運動狀態(tài)，避免金屬屑在減速器內(nèi)部沉積，提升故障監(jiān)控效果。