屈思博，張?緯，由長(zhǎng)福,2

顆粒群在燃燒環(huán)境中的曳力研究

屈思博1，張?緯1，由長(zhǎng)福1,2

(1. 清華大學(xué)能源與動(dòng)力工程系，熱科學(xué)與動(dòng)力工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084；2. 清華大學(xué)山西清潔能源研究院，太原 030000)

對(duì)顆粒半徑為100～500mm、體積分?jǐn)?shù)為0.075～0.35的固定隨機(jī)分布的顆粒群在303K與1173K環(huán)境溫度下的受力及相應(yīng)的流場(chǎng)進(jìn)行模擬．模擬工具為虛擬區(qū)域法．它采用統(tǒng)一形式的結(jié)構(gòu)化網(wǎng)格描述顆粒與氣體，并實(shí)現(xiàn)兩相的相互作用，兼有較高的準(zhǔn)確度與效率．相關(guān)結(jié)果顯示，顆粒群在燃燒環(huán)境中形成的流場(chǎng)更穩(wěn)定，這主要是由于氣體在燃燒環(huán)境中的黏性更大，導(dǎo)致雷諾數(shù)更?。瑫r(shí)，顆粒群在燃燒環(huán)境中所受黏性力更大．在較低與較高的雷諾數(shù)下，顆粒群在燃燒環(huán)境中的曳力分別更大與更?。S著顆粒尺寸增大，顆粒群的燃燒強(qiáng)度與曳力均增大；隨著顆粒體積分?jǐn)?shù)增大，顆粒群的燃燒強(qiáng)度減小但是曳力增大．最后，本研究對(duì)于顆粒群的終端速度、曳力系數(shù)以及氣-固動(dòng)量交換系數(shù)提出了在燃燒環(huán)境中的修正．

多相流；煤燃燒；顆粒群；曳力；虛擬區(qū)域法

顆粒群是氣-固兩相流在顆粒材料輸送、煤粉燃燒等工業(yè)應(yīng)用中重要的顆粒形式．尤其在稠密的顆粒流中，不僅顆粒與氣體存在相互作用，顆粒之間的相互作用也不可忽視．以流化床[1]為代表的燃燒環(huán)境中存在氣-固相間動(dòng)量交換、傳質(zhì)與傳熱，氣體物性變化與化學(xué)反應(yīng)，導(dǎo)致多物理場(chǎng)的復(fù)雜變化，因此，受到了較多的關(guān)注[2-6]．目前，傳統(tǒng)的顆粒曳力模型及其相應(yīng)的修正[7-12]一般是基于冷態(tài)環(huán)境所提出，而基于燃燒環(huán)境的研究仍然較少．近年來，相關(guān)領(lǐng)域內(nèi)關(guān)注斯蒂芬流或孔隙對(duì)顆粒曳力的影響的研究較多.在斯蒂芬流方面，Kurose等[13]模擬了一個(gè)具有斯蒂芬流的固定顆粒，發(fā)現(xiàn)斯蒂芬流傾向于減少曳力. Jayawickrama等[14]基于類似的工作，進(jìn)一步指出斯蒂芬流主要影響曳力的黏性力部分．Chen等[15]模擬了具有斯蒂芬流的固定隨機(jī)分布顆粒群，發(fā)現(xiàn)斯蒂芬流對(duì)曳力的削弱程度隨著顆粒體積分?jǐn)?shù)的增大而減小.在孔隙方面，Bhattacharyya等[16]基于滲透率模型，發(fā)現(xiàn)曳力隨著達(dá)西數(shù)的增大而減小. Wittig等[17-18]采用顆粒群描述多孔顆粒，發(fā)現(xiàn)在低與高雷諾數(shù)下，曳力隨著孔隙率的增大而分別減小與增大．這些研究針對(duì)性地考查了特定顆粒形式對(duì)曳力的影響，為后續(xù)相關(guān)工作的展開奠定了基礎(chǔ)．但是，不足在于基于冷態(tài)環(huán)境的模擬并不完全符合一般的物理事實(shí)，即：顆粒在高溫或燃燒環(huán)境中才產(chǎn)生明顯的出流，或產(chǎn)生大量的孔隙．Watanabe等[19]通過加熱煤顆粒實(shí)驗(yàn)觀察到明顯的熱解氣體從顆粒表面流出，并通過模擬在高溫環(huán)境中的顆粒發(fā)現(xiàn)，顆粒在氣體逆流方向產(chǎn)生表面出流的情況下，曳力趨于減小．Zhang等[20]模擬了處于燃燒狀態(tài)且具有斯蒂芬流的單個(gè)顆粒，發(fā)現(xiàn)在燃燒環(huán)境中斯蒂芬流反而增大了曳力．本研究之前通過一系列在燃燒環(huán)境中單個(gè)顆粒的模擬，發(fā)現(xiàn)在燃燒環(huán)境中由于氣體黏性更大，顆粒雷諾數(shù)更小，生成的流場(chǎng)更穩(wěn)定，總體上呈現(xiàn)出所受黏性力更大、壓力更小的趨勢(shì)．因此，加熱或燃燒環(huán)境中顆粒曳力的模擬不可直接套用基于冷態(tài)環(huán)境所提出的傳統(tǒng)曳力模型，前者涉及的系統(tǒng)內(nèi)物質(zhì)物性變化的影響是不可忽視的.

如Watanabe等[19]的實(shí)驗(yàn)所確認(rèn)，顆粒出流主要發(fā)生在揮發(fā)分釋放階段，而作為煤燃燒主體的焦炭燃燒階段，相對(duì)于燃燒環(huán)境，顆粒出流對(duì)曳力的影響較?。瑫r(shí)，受Wittig等[17-18]的工作所啟發(fā)，多孔顆?？梢圆捎米鳛橐粋€(gè)整體的具有一定體積分?jǐn)?shù)的顆粒群來描述．因此，研究燃燒環(huán)境中顆粒群的曳力特性對(duì)于進(jìn)一步揭示氣-固兩相間的復(fù)雜相互作用機(jī)制具有重要意義．為了詳細(xì)地描述顆粒(群)，采用小于顆粒尺寸的單位網(wǎng)格的完全解析模擬可以滿足最高的精度的要求．例如浸入式邊界法[21]在以往實(shí)踐[14-15]中表現(xiàn)良好．近年來，從浸入式邊界法發(fā)展而來的虛擬區(qū)域法[22-24]由于其對(duì)顆粒更嚴(yán)格的剛體約束，在對(duì)攜帶大量相間密度差大的顆粒(群)的氣-固(反應(yīng))流中更表現(xiàn)出極大的優(yōu)勢(shì)．因此，本研究采用虛擬區(qū)域法對(duì)顆粒半徑為100～500mm、體積分?jǐn)?shù)為0.075～0.35的固定隨機(jī)分布的顆粒群在303K與1173K環(huán)境溫度下的受力及相應(yīng)的流場(chǎng)進(jìn)行模擬，并對(duì)顆粒群所形成的物理場(chǎng)與曳力的組成(黏性力與壓力)進(jìn)行分析．最后，以常溫環(huán)境為基準(zhǔn)，本研究對(duì)燃燒環(huán)境中顆粒群的顆粒終端速度、曳力系數(shù)與氣-固動(dòng)量交換系數(shù)提出相應(yīng)的修正，以期為指導(dǎo)清潔高效的燃燒技術(shù)設(shè)計(jì)提供參考．

1?數(shù)值方法

對(duì)于氣-固反應(yīng)流多物理場(chǎng)的描述，本研究采用納維-斯托克斯方程來求解速度與壓力場(chǎng)、傳質(zhì)與傳熱方程來求解組分與溫度場(chǎng)．類似地，本研究基于經(jīng)典的氣-固傳熱理論來描述顆粒的溫度變化．如圖1所示，虛擬區(qū)域法采用結(jié)構(gòu)化網(wǎng)格統(tǒng)一描述顆粒與氣體．被藍(lán)色圓點(diǎn)(顆粒)包圍的區(qū)域與被紅色三角點(diǎn)(氣體)包圍的區(qū)域分別適用于顆粒與氣體控制方程. 一系列黑色圓點(diǎn)(相界面)被用于描述顆粒與氣體的相界面. 本研究通過在相界面所在網(wǎng)格的控制方程中設(shè)置相應(yīng)的源相來描述氣-固兩相的相互作用．虛擬區(qū)域法通過在氣體控制方程內(nèi)設(shè)置基于顆粒控制方程得到的“虛擬源項(xiàng)”，使得氣體的控制方程(組)適用于全計(jì)算域，即通過添加了“虛擬源相”的氣體控制方程求解的顆粒物理量與從顆?？刂品匠糖蠼獾玫降念w粒物理量在數(shù)值上一致．氣體各控制方程“虛擬源項(xiàng)”的具體形式、方程離散與數(shù)值求解算法，以及虛擬區(qū)域法對(duì)氣-固反應(yīng)流的動(dòng)力學(xué)與化學(xué)反應(yīng)描述的準(zhǔn)確性與網(wǎng)格無關(guān)性驗(yàn)證可參考文獻(xiàn)[23-24]．

圖1?虛擬區(qū)域法計(jì)算域[11]

值得說明的是，本研究基于理想氣體假設(shè)描述氣體的密度、壓強(qiáng)與溫度關(guān)系，如式(1)所示．

為簡(jiǎn)化問題，本研究采用表1[11-12]所列適用范圍較廣的3個(gè)單步不可逆反應(yīng)來描述由焦炭顆粒(群)燃燒導(dǎo)致的系統(tǒng)內(nèi)的物質(zhì)變化．

表1?化學(xué)反應(yīng)動(dòng)力學(xué)參數(shù)

Tab.1?Kinetic parameters for chemical reactions

2?模擬設(shè)置與流場(chǎng)分析

2.1?顆粒位置與邊界條件

圖2?模擬設(shè)置

2.2?流場(chǎng)分析

圖3?顆粒群所形成的速度場(chǎng)(e＝0.15)

圖4?顆粒群所形成的壓力場(chǎng)(e＝0.15)

表2?顆粒雷諾數(shù)(＝0.15)

Tab.2?Particle Reynolds number(e＝0.15)

圖5?顆粒群所形成的溫度場(chǎng)(e＝0.15)

圖6?顆粒群所形成的O2質(zhì)量分?jǐn)?shù)場(chǎng)(e＝0.15)

表3?顆粒平均溫度與燃燒速率(＝0.15)

Tab.3 Average particle temperature and combustion rate (e＝0.15)

此外，根據(jù)壓力場(chǎng)的變化可以發(fā)現(xiàn)，隨著增大，無論是在常溫還是燃燒環(huán)境中，流場(chǎng)的壓降均增大．因此，可知更加密集的顆粒群形成了更大的氣流阻力，或換言之，受到了更大的曳力，從而導(dǎo)致了更大的流場(chǎng)壓降．

圖7?顆粒群所形成的速度場(chǎng)(rp＝300mm)

圖8?顆粒群所形成的壓力場(chǎng)(rp＝300mm)

表4?顆粒雷諾數(shù)(p＝300mm)

Tab.4?Particle Reynolds number(rp＝300mm)

圖9?顆粒群所形成的溫度場(chǎng)(rp＝300mm)

圖10?顆粒群所形成的O2質(zhì)量分?jǐn)?shù)場(chǎng)(rp＝300mm)

表5?顆粒平均溫度與燃燒速率(p＝300mm)

Tab.5 Average particle temperature and combustion rate (rp＝300mm)

3?顆粒曳力組成分析及經(jīng)驗(yàn)修正

3.1?顆粒曳力組成分析

如式(2)所示，本研究通過分別積分顆粒表面的黏性應(yīng)力和壓強(qiáng)，并將該兩者求和以計(jì)算顆?？偟囊妨Γ艿筋w粒尺寸和體積分?jǐn)?shù)的影響，在不同的氣流速度和環(huán)境溫度條件下，氣體施加于顆粒群的黏性力與壓力對(duì)其整體的曳力貢獻(xiàn)程度是不同的．

（a）p＝100mm

（b）p＝200mm

（c）p＝400mm

（d）p＝500mm

圖11?顆粒平均曳力及其組成(e＝0.15)

總的來說，本研究的定量結(jié)果符合前文2.2節(jié)的推測(cè)，且對(duì)于顆粒群在燃燒環(huán)境中曳力特性變化及顆粒尺寸影響的發(fā)現(xiàn)與之前單顆粒的情況是類似的，即在燃燒環(huán)境中，較大的氣體黏性保持了顆粒(群)流較高的穩(wěn)定性，且使得顆粒(群)所受曳力的氣體黏性力部分增大．

（a）＝0.075

（b）＝0.15

（c）＝0.25

（d）＝0.35

圖12?顆粒平均曳力及其組成(rp＝300mm)

表6?顆粒平均曳力及其組成(g為0.5m/s，2.5m/s)

Tab.6 Average particle drag force composition(ug＝0.5m/s，2.5m/s)

3.2?顆粒曳力修正

（a）t,c/t,0

（b）d,c/d,0

（c）c/0

圖13?顆粒群動(dòng)力學(xué)特性變化

Fig.13?Particle group dynamic characteristic variations

4?結(jié)論與展望

(1)相比于常溫環(huán)境，顆粒在燃燒環(huán)境中形成的流場(chǎng)更加穩(wěn)定．這是由于后者的氣體黏性更大，導(dǎo)致在相同的顆粒尺寸、氣流速度與體積分?jǐn)?shù)條件下p更小．

(2)通過顆粒流的形態(tài)與顆粒曳力組成的分析可知，在p較小的情況下，常溫與燃燒環(huán)境中的流場(chǎng)均較為穩(wěn)定，兩種環(huán)境中顆粒所受的壓力相差不大，但是燃燒環(huán)境中的顆粒受到更大的黏性力與曳力；在p的較大的情況下，盡管燃燒環(huán)境中的顆粒所受黏性力更大，但是更加穩(wěn)定的流場(chǎng)使得顆粒受到的更小的壓力與曳力．

(3)隨著顆粒尺寸增大，顆粒群的燃燒強(qiáng)度與曳力均增大；隨著增大，顆粒群的燃燒強(qiáng)度有所降低，但是曳力增大．

本研究分別將上述若干算例中顆粒群在兩種環(huán)境中的終端速度、曳力系數(shù)與氣-固動(dòng)量交換系數(shù)之比作為修正值，提出了相應(yīng)的修正公式，以期為指導(dǎo)清潔高效的燃燒技術(shù)設(shè)計(jì)提供參考.

值得說明的是，本研究對(duì)單個(gè)顆粒采用簡(jiǎn)單的二維實(shí)心純碳球體假設(shè)，未考慮實(shí)際工業(yè)應(yīng)用中顆粒的三維狀態(tài)、多組分、孔隙結(jié)構(gòu)以及非球體形態(tài)的影響；同時(shí)，對(duì)燃燒環(huán)境溫度的設(shè)置局限于平均溫度較低、溫度范圍較窄的循環(huán)流化床環(huán)境．因此，對(duì)于具有更加復(fù)雜組分、結(jié)構(gòu)與形態(tài)的單個(gè)顆?；蝾w粒群，以煤粉鍋爐為代表的具有更高溫度的燃燒環(huán)境以及本研究所得到的模型在不同工況中的適用性，相關(guān)的研究仍然需要在未來得到進(jìn)一步探索與完善．

［1］ Alobaid F，Almohammed N，F(xiàn)arid M，et al. Progress in CFD simulations of fluidized beds for chemical and energy process engineering[J].，2022，92：100930.

［2］ Stroh A，Alobaid F，von Bohnstein M，et al. Numerical CFD simulation of 1MWth circulating fluidized bed using the coarse grain discrete element method with homogenous drag models and particle size distribution[J].，2018，169：84-93.

［3］ Xie J，Zhong W，Shao Y，et al. Coupling of CFD-DEM and reaction model for 3D fluidized beds[J].，2019，353：72-83.

［4］ Mu L，Buist K A，Kuipers J A M，et al. Hydrodynamic and heat transfer study of a fluidized bed by discrete particle simulations[J].，2020，8(4)：463.

［5］ Wu Y，Liu D，Ma J，et al. Three-dimensional Eulerian–Eulerian simulation of coal combustion under air atmosphere in a circulating fluidized bed combustor[J].，2017，31(8)：7952-7966.

［6］ Wu Y，Liu D，Ma J，et al. Effects of gas-solid drag model on Eulerian-Eulerian CFD simulation of coal combustion in a circulating fluidized bed[J].，2018，324：48-61.

［7］ Schiller L，Naumann A. A drag coefficient correlation[J].，1935，77：318-320.

［8］ Wen C，Yu Y. Mechanics of fluidization [J].，1966，62：100-111.

［9］ Gibilaro L，Di Felice R，Waldram S，et al. Generalized friction factor and drag coefficient correlations for fluid-particle interactions[J].，1985，40(10)：1817-1823.

［10］ Syamlal M，O'Brien T. Derivation of a Drag Coefficient from Velocity-Voidage Correlation[R]. United States：National Energy Technology Laboratory，1987.

［11］ Gidaspow D. Multiphase flow and fluidization：Continuum and kinetic theory description[J].，1994，55(2)：207-208.

［12］ Beetstra R，van der Hoef M，Kuipers J. Drag force of intermediate Reynolds number flow past mono‐and bidisperse arrays of spheres[J].，2007，53(2)：489-501.

［13］ Kurose R，Makino H，Komori S，et al. Effects of outflow from the surface of a sphere on drag，shear lift，and scalar diffusion[J].，2003，15(8)：2338-2351.

［14］ Jayawickrama T，Haugen N，Babler M，et al. The effect of Stefan flow on the drag coefficient of spherical particles in a gas flow[J].，2019，117：130-137.

［15］ Chen X，Du S，Zhao L，et al. Effect of Stefan flow on the drag force in flow past random arrays of spheres [J].，2021，412(7)：128691.

［16］ Bhattacharyya S，Dhinakaran S，Khalili A. Fluid motion around and through a porous cylinder[J].，2006，61(13)：4451-4461.

［17］ Wittig K，Golia A，Nikrityuk P. 3D numerical study on the influence of particle porosity on heat and fluid flow[J].，2012，12(2/3)：207-219.

［18］ Wittig K，Nikrityuk P，Richter A. Drag coefficient and Nusselt number for porous particles under laminar flow conditions[J].，2017，112：1005-1016.

［19］ Watanabe M，Yahagi J. Effects of nonuniform outflow and buoyancy on drag coefficient acting on a spherical particle[J].，2017，5(4)：99-110.

［20］ Zhang H，Luo K，Haugen N，et al. Drag force for a burning particle[J].，2020，217：188-199.

［21］ Uhlmann M. An immersed boundary method with direct forcing for the simulation of particulate flows[J].，2005，209(2)：448-476.

［22］ Glowinski R，Pan T，Hesla T，et al. A distributed Lagrange multiplier/fictitious domain method for particulate flows[J].，1999，25(5)：755-794.

［23］ Zhang L，Liu K，You C. Fictitious domain method for fully resolved reacting gas-solid flow simulation[J].，2015，299：215-228.

［24］ Qu S，You C. Direct numerical simulation(DNS)of alkali metals released during char combustion[J].，2019，255：115763.

［25］ Linstrom P，Mallard W. NIST Chemistry WebBook，NIST Standard Reference Database[EB/OL]. http:// webbook.nist.gov，2022.

Drag Force on Particle Group Under Combustion Condition

Qu Sibo1，Zhang Wei1，You Changfu1, 2

(1. Key Laboratory for Thermal Science and Power Engineering of Ministry of Education，Department of Energy and Power Engineering，Tsinghua University，Beijing 100084，China；2. Shanxi Research Institute for Clean Energy，Tsinghua University，Taiyuan 030000，China)

In this study，a set of gas-solid flow simulations with fixed randomly distributed particles was implemented. The particle radius was set as 100—500mm，the particle volume fraction was set as 0.075—0.35，and the temperature was set as 303K and 1173K，representing the normal and combustion conditions，respectively. The fictitious domain method was adopted as the numerical tool，which uses a unified structured mesh to describe particle and gas and realizes two-phase interaction，with high accuracy and efficiency. Relevant results show that the flow field formed by particle group in the combustion environment is more stable，which is mainly due to the smaller particle Reynolds number resulting from larger gas viscosity. The viscous force of particle group is greater in the combustion environment. The drag force of particle group in the combustion environment is large at low Reynolds numbers and small at high Reynolds numbers. With the increase of particle size，both the combustion intensity and drag force of particle group increase；with the increase of particle volume fraction，the combustion intensity of particle group decreases but the drag force increases. Finally，the corrections of particle terminal velocity，drag force coefficient and gas-solid momentum exchange coefficient in the combustion environment were proposed.

multi-phase flow；coal combustion；particle group；drag force；fictitious domain method

TK11

A

1006-8740(2023)06-0624-11

10.11715/rskxjs.R202209004

2022-09-02.

國家自然科學(xué)基金資助項(xiàng)目(51761125012)；華能集團(tuán)總部科技項(xiàng)目——基礎(chǔ)能源科技研究專項(xiàng)資助項(xiàng)目(HNKJ20-H50).

屈思博(1995—??)，男，博士研究生，qusibo2013@163.com.

張?緯，女，博士，助理研究員，zw_sim@163.com.

(責(zé)任編輯：隋韶穎)