大顆粒氣力輸運(yùn)的CFD與EDEM耦合計(jì)算方法

王歡， 孫科， 任博揚(yáng)， 宋江濤

(中國飛行試驗(yàn)研究院發(fā)動機(jī)所， 西安 710089)

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

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

1 冰雹拋射機(jī)構(gòu)

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

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

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

拋雹裝置在運(yùn)行時，首先打開壓縮空氣供應(yīng)系統(tǒng)的調(diào)節(jié)閥，使得在拋射管內(nèi)部形成穩(wěn)定的流場。后將冰雹投放機(jī)構(gòu)的閥門打開，將冰雹持續(xù)不斷的投放至發(fā)射管入口處，冰雹受到氣流作用產(chǎn)生加速。

2 計(jì)算方法

2.1 流體相計(jì)算方法

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

由于冰雹顆粒與流體相之間存在相互作用，在流體相的守恒方程中額外添加一個體積分?jǐn)?shù)項(xiàng)來修正。修正后的動量方程為

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

(1)

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

動量匯S計(jì)算公式為

(2)

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

流體項(xiàng)體積分?jǐn)?shù)計(jì)算公式為

(3)

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

2.2 冰雹曳力計(jì)算方法

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

(4)

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

(5)

自由流阻力模型中阻力系數(shù)計(jì)算式為

(6)

2.3 耦合計(jì)算方法

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

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

3 拋射原理分析

3.1 計(jì)算模型

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

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

3.2 冰雹加速過程分析

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

表1 基本工況參數(shù)設(shè)置

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

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

4 影響因素分析

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

4.1 氣源壓力影響

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

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

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

4.2 拋射管徑影響

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

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

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

4.3 拋射管長影響

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

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

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

4.4 冰雹流量影響

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

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

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

5 結(jié)論

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

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

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

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