減速器內部金屬屑運動仿真分析

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

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

0 引言

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

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

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

1 VOF-DPM 耦合模型驗證

1.1 顆粒受力方程

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

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

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

式中：a1、a2、a3為常數，在不同雷諾數范圍內有不同的值[14]

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

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

對于球形顆粒，其轉動慣量Ip為

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

1.2 編程計算

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

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

圖1 計算流程

1.3 CFD仿真

由于3 維減速器中顆粒運動軌跡復雜，本文以簡化的2 維減速器進行分析。建立齒輪箱平面如圖2 所示。齒輪箱長434 mm，高300 mm，壓力出口寬20 mm，對其劃分非結構化網格以適應復雜流動狀態，以壓力出口的質量流量為監測量，經過網格無關性驗證，網格數量最終確定為18742，網格質量最大偏斜度為0.5812，符合CFD仿真計算標準。

圖2 齒輪箱平面

齒輪模數、齒數、轉速、浸油深度、顆粒物性參數等齒輪及滑油的主要參數設置見表1。

表1 齒輪及滑油的主要參數

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

圖3 流場計算模型初始狀態

1.4 對比驗證

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

圖4 編程計算方法與CFD方法計算的顆粒運動狀態

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

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

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

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

2.1 虛質量力

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

式中：Cvm為虛質量系數，默認值為0.5。

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

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

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

圖6 有、無附加虛質量力的顆粒運動狀態及油液分布對比

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

2.2 壓力梯度力

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

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

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

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

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

2.3 Saffman升力

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

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

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

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

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

2.4 Magnus升力

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

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

對于旋轉升力系數CRL，常用以下2種計算方法。

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

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

自旋參數S定義為

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

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

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

2.5 熱泳力

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

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

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

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

2.6 Brownian力

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

式中：δij為Kronecker δ函數，且有

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

式中：ξ為與單位方差無關的高斯隨機數，均值為零。

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

3 結論

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

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

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

（4）本文的研究方法方法可以拓展到3 維減速器仿真，以實現更準確的預測減速器內部金屬屑運動狀態，避免金屬屑在減速器內部沉積，提升故障監控效果。