大顆粒氣力輸運的CFD與EDEM耦合計算方法

王歡， 孫科， 任博揚， 宋江濤

(中國飛行試驗研究院發動機所， 西安 710089)

為保證航空發動機運行過程中的安全可靠，國內外均制定了相應的適航標準[1-2]。中國民用航空局制定的《航空發動機適航規定》[3]明確提出，民用航空發動機在突然遭遇濃度達到CCAR 33部附錄B中定義的冰雹時，在其整個規定的工作包線范圍內仍具有可接受的工作能力。發動機可接受的工作能力是指在任何連續30 s的降雹周期內發動機不熄火、不降轉、不發生持續或不可恢復的喘振或失速，或不失去加速或減速的能力。還必須證明吸入之后沒有不可接受的機械損壞、不可接受的功率或推力損失或其他不利的發動機異常情況。大氣冰雹濃度為3 660～4 570 m時最高為10%。為對發動機安全性進行充分考核，需進行發動機吸入冰雹的驗證試驗。有學者對航空發動機吸雨吸雹進行了適航審定符合性研究[4]。根據發動機進氣流量大小及冰雹水含量分析得到冰雹流量，在進行吸雹試驗時，需在短時間內向發動機入口拋射上以萬記的冰雹。由于在試驗時，無法真實模擬真實大氣中的冰雹尺寸分布，可用CCAR 33部附錄B中的平均直徑代替。國外發動機進行連續吸雹試驗所用的冰雹直徑為12.7 mm左右。連續冰雹的速度及打擊位置不是關鍵試驗條件，但應使得連續冰雹能夠順利全部進入發動機中。使用壓縮空氣氣力輸運大量冰雹，在拋射管內部冰雹加速運動，離開拋射管進入發動機進口。類似壓縮氣體輸運顆粒問題廣泛應用于化工、冶金、制藥、熱力發電等行業，國內外學者針對這種問題進行了大量的仿真研究。文獻[5]中應用計算流體力學-離散單元法(computational fluid dynamics-discrete element method，CFD-DEM)耦合方法進行了離心泵內的固液兩相流動模擬。文獻[6]中應用CFD-DEM耦合方法進行氣力輸送彎管內懸浮氣固兩相流動的數值模擬。文獻[7]中使用CFD-DEM方法對一種氣力輸送裝置進行了研究，文獻[8]中使用CFD-DEM方法進行了連續鼓泡流化床速度及壓降的模擬研究。文獻[9-10]針對煤層氣井大粒徑巖屑運動規律進行了研究。文獻[11]中應用CFD-DEM方法研究了粗糙壁面顆粒趨壁沉積過程，上述研究中針對顆粒直徑僅為2～3 mm的小顆粒，且其中流體相速度較低，均視為不可壓縮流動。文獻[12]中研究了單個大顆粒在流場中的沉降。連續冰雹拋射是一個可壓縮流氣力輸運大量大尺寸顆粒的問題，對連續冰雹這類大尺寸顆粒的氣力輸運問題的研究還未見報道。

針對較大尺寸冰雹的連續拋射問題，分析冰雹拋射裝置的氣力輸運過程，其中氣相為可壓縮流動，固相冰雹顆粒尺寸及流量都較大，應用CFD-EDEM耦合計算方法對冰雹拋射原理進行計算研究。

1 冰雹拋射機構

拋雹裝置采用壓縮空氣進行冰雹顆粒的氣力輸運，采用如圖1所示的拋射管對冰雹顆粒進行運輸。如圖2所示為拋雹炮管示意圖。在拋射管入口處投入冰雹，由于拋射管內部形成了較為穩定的摩擦管流，管內較快速度的氣流與顆粒間存在速度差，因此存在動量交換，顆粒受到高速氣流的作用，在管道內部加速運動，直至運動至拋射管出口。

圖1 拋雹炮管Fig.1 Hail ejection tube

圖2 拋雹炮管示意圖Fig.2 Diagram of the hail ejection tube

拋雹裝置在運行時，首先打開壓縮空氣供應系統的調節閥，使得在拋射管內部形成穩定的流場。后將冰雹投放機構的閥門打開，將冰雹持續不斷的投放至發射管入口處，冰雹受到氣流作用產生加速。

2 計算方法

2.1 流體相計算方法

由于冰雹顆粒大，為實現冰雹顆粒的較大速度，所用的壓縮空氣壓力需較大。因此流體相計算采用可壓流計算方法?？諝饷芏认嗖捎美硐霘怏w狀態方程進行求解。采用Fluent軟件求解流量方程、動量方程、能量方程。采用標準k-ε湍流模型。

由于冰雹顆粒與流體相之間存在相互作用，在流體相的守恒方程中額外添加一個體積分數項來修正。修正后的動量方程為

-?p+ρfεfg+μfεf?2uf-S

(1)

式(1)中：εf為流體項的體積分數，取1；ρf為流體密度，kg/m3；uf為流體速度，m/s；t為時間，s；μf為流體動力黏度，N·s/m2；p為氣流靜壓，Pa；g為重力加速度，m/s2；S為動量匯，即顆粒和流體間的體積作用力。氣固兩相之間的耦合由相對運動產生阻力的動量匯大小來實現。

動量匯S計算公式為

(2)

式(2)中：Fs,i為固體顆粒i的曳力，N；V為CFD網格單元的體積，m3；n為控制體中所含的顆粒數。

流體項體積分數計算公式為

(3)

式(3)中：Vs,i為第i個固體顆粒的體積。

2.2 冰雹曳力計算方法

使用自由流阻力模型對顆粒的受力情況進行分析，受力和顆粒與流體間速度差、流體密度及阻力系數等有關。顆粒受力Fs計算公式[13-14]為

(4)

式(4)中：CD為阻力系數；ds為冰雹顆粒直徑，m；us為冰雹顆粒速度，m/s；阻力系數取決于顆粒流動雷諾數Res，其定義為

(5)

自由流阻力模型中阻力系數計算式為

(6)

2.3 耦合計算方法

考慮冰雹拋射時，冰雹體積分數較大，從而對流場分布產生影響，因此需要進行雙向耦合計算，迭代求解冰雹顆粒在拋射管內部的運動過程。如圖3所示為耦合計算流程圖，首先使用計算流體力學軟件對拋射管內部的流場分布進行求解，將流場信息傳遞至離散元分析軟件中，在拋射管入口處等間隔時刻放置冰雹，從而冰雹在高速氣流作用下運動，使用離散元方法對顆粒受力及速度、位置進行計算；再將顆粒造成的體積分數變化及動量匯等信息傳遞至流場計算中，以此形成計算流體力學與離散單元法計算的雙向耦合。待冰雹拋射速度計算結果穩定時，即可停止迭代耦合計算。

圖3 CFD與離散單元法耦合計算流程示意圖Fig.3 CFD and discrete element method coupling calculation diagram

3 拋射原理分析

3.1 計算模型

針對冰雹拋射裝置，進行CFD與EDEM耦合的冰雹顆粒氣力輸運計算。如圖4所示為冰雹拋射管的流體域計算模型。管道進口處設置穩定的壓力入口條件，拋射管的壁面為絕熱壁面，在拋射管出口處劃分一個大圓柱形的膨脹區域，計算模型出口為壓力出口，模擬壓縮空氣在拋射管出口處的膨脹過程。

圖4 冰雹拋射過程計算模型Fig.4 Computational model of hail projectile process

3.2 冰雹加速過程分析

通過基本工況條件下的計算，對冰雹在拋射管內的運動過程及加速過程進行分析。該工況的具體參數如表1所示。拋射管選擇為5 m，在拋射管出口處劃分一個直徑為2 m、長為2 m的圓柱，對拋射管出口處的射流進行計算。

表1 基本工況參數設置

如圖5所示為基本算例對應的計算結果，其中對單一流體相、流體/顆粒兩相耦合這兩種條件的計算結果進行了對比。圖5所示為單一流體相計算速度分布，及流體/顆粒兩相耦合計算中流體相與顆粒相在拋射管內部和其延長段上的速度分布。對于流體單相計算，在等截面摩擦管流作用下，氣流速度隨著管長發展有明顯的增長。在拋射管出口，氣流噴射至空氣中，速度逐漸衰減。在兩相耦合計算中，由于顆粒的存在對氣流速度的發展形成阻礙，因此管道中與單一流體相情況相比，氣流流速較小，在拋射管出口處，氣流速度迅速降低。對于顆粒而言，由于在拋射管入口處，氣流速度與顆粒速度差較大，因此輸運曳力較大，顆粒有較大的速度增長趨勢，隨著顆粒速度增大，與氣流速度差逐漸減小，從而曳力減小，因此顆粒速度基本沒有繼續增大的趨勢。

圖5 流體相及顆粒相速度分布Fig.5 velocity distribution of fluid phase and particle phase

4 影響因素分析

冰雹拋射規律受到多個因素影響。拋射氣源壓力、拋射管徑、拋射管長、冰雹拋射流量等均會影響冰雹拋射規律及拋射速度。采用拋射管外射流段處的顆粒平均速度，作為評價冰雹拋射速度大小的指標。

4.1 氣源壓力影響

針對不同氣源壓力進行拋射規律計算，對基本工況中的拋射氣源壓力進行更改，其余參數均與基本算例保持一致。分別在拋射壓力為160、180、200、220、240 kPa條件下進行計算，從而得到不同拋射壓力對冰雹拋射速度的影響規律。如圖6所示為不同氣源壓力條件下的冰雹拋射速度變化規律。由于氣源壓力大時，氣流輸運冰雹的動力較強，因此隨著氣源壓力的增大，冰雹拋射速度逐漸增大。如圖7所示為不同氣源壓力條件下的空氣流量，對不帶顆粒氣流單相流動與帶顆粒的氣-固兩相流動情況下的空氣流量進行了對比。由于固體顆粒的存在對氣流相形成了阻力，因此氣-固兩相流動情況下的空氣流量與單相流動情況下相比明顯偏小。隨著氣源壓力的增大，空氣流量基本線性增大。

圖6 不同氣源壓力下冰雹速度Fig.6 Hail velocity under different air source pressure

圖7 不同氣源壓力下空氣流量Fig.7 Air flowrate under different air source pressure

4.2 拋射管徑影響

為分析拋射管徑對顆粒輸運的影響，分別進行了拋射管內徑為50、60、70、80、100 mm條件下的顆粒氣力輸運計算模擬，其余參數設置均與基本算例保持一致。如圖8所示為不同管徑條件下冰雹顆粒速度的計算結果對比，隨著管徑的增大，冰雹速度有較大的增長。這是由于單位時間內冰雹顆粒數量保持一致，管徑較大的條件下，冰雹顆粒相體積分數較小，氣流提供給冰雹的動力較強。但隨著拋射管徑的增大，空氣流量平方級增長，如圖9所示為不同管徑條件下的氣流流量計算結果。如保持長時間持續冰雹拋射，則需要有大量的壓縮空氣，造成較大的資源消耗。因此進行冰雹拋射管徑的選型時，無須一味增大管徑。

圖8 不同管徑冰雹顆粒速度Fig.8 Hail particle velocity with different pipe diameters

圖9 不同管徑對應空氣流量Fig.9 Air flowrate with different pipe diameters

4.3 拋射管長影響

進行不同拋射管長條件下的冰雹顆粒輸運計算，對管長造成的拋射運動規律變化進行分析。分別進行了拋射管長度為3、4.5、5、5.5、6 m的計算，其余參數設置均與基本算例保持一致，與基本算例進行對比。如圖10所示為不同管長條件下冰雹拋射速度計算結果，隨著管長的變化，冰雹拋射速度無明顯的變化。如圖11所示為不同管長條件下氣流流量的結果，隨著管長的增大，空氣流量逐漸減小，這是由于管長增大時，氣流所受的壁面摩阻作用更大。另外由于拋射管長度增加，冰雹顆粒加速運動路徑增長，兩個因素綜合作用下，冰雹顆粒速度隨著管長的增長并無明顯變化。冰雹拋射管較長時，對拋射管的加工要求也較高，因此綜合考慮，拋射管長度無須太長。

圖10 不同管長條件下冰雹顆粒速度Fig.10 Hail particle velocity under different pipe lengths

圖11 不同管長條件下氣流流量Fig.11 Air flowrate under different pipe lengths

4.4 冰雹流量影響

通過進行不同拋雹流量條件下的計算，對冰雹流量因素進行分析。分別進行冰雹拋射流量為200、400、600、800 個/s，其余參數設置均與基本算例保持一致，與基本算例進行對比。如圖12所示為不同冰雹拋射流量條件下的計算結果對比，隨著冰雹拋射流量的增大，冰雹速度逐漸減小，這是由于隨著冰雹顆粒流量的增大，冰雹顆粒相體積分數增大，在入口壓力同樣的條件下，單個冰雹所分的動力變小。如圖13為不同冰雹拋射流量條件下氣流流量的計算結果，由于冰雹顆粒的增多，造成對氣流相的阻力增大，因此隨著冰雹顆粒的增多，氣流流量逐漸減小。

圖12 不同拋雹量下冰雹顆粒速度Fig.12 Hail particle velocity with different hail flowrate

圖13 不同拋雹量下氣流流量Fig.13 Air flowrate with different hail flowrate

5 結論

針對壓縮空氣輸運大尺寸冰雹顆粒的氣力輸運問題，建立耦合迭代計算方法，并針對不同工況條件下進行計算，為冰雹持續拋射裝置的設計提供支撐。主要得到以下結論。

(1)建立了適用于可壓縮流輸運大尺寸顆粒的雙向耦合迭代計算方法，通過計算流體力學與離散單元法的迭代計算，得到冰雹在拋射管內的運動過程。

(2)通過進行不同氣源壓力、不同拋射管徑、不同拋射管長、不同冰雹拋射流量條件下的計算。形成了以下結果：氣源壓力的增長會使得冰雹速度增大，所消耗氣流流量逐漸增長；拋射管徑的增大會使得冰雹速度增大，所消耗氣流流量平方級增長；拋射管長的增大對冰雹速度的影響微乎其微，所消耗氣流流量有小幅度減??；冰雹流量的增大會使得冰雹速度減小，所消耗氣流流量減小。

(3)通過分析上述計算結果，可初步選擇管徑為45 mm、管長為4 m的拋射管進行冰雹的氣力輸運。