劉 彬,周武洲,趙朋程,李 瑞
(1.燕山大學 信息科學與工程學院,河北 秦皇島 066004;2.燕山大學 電氣工程學院,河北 秦皇島 066004)
隨著近些年水泥行業(yè)的高速發(fā)展,熟料燒結過程中的高能耗和高排放的問題,引起越來越多的關注[1-2]。每生產一噸水泥,一般要消耗3 GJ能量,同時排放出比所生產水泥總量更多的溫室氣體[3]。研究回轉窯工況對窯內流場、溫度分布以及煤粉燃燒的影響,可以為水泥廠的節(jié)能減排、改進水泥生產工藝提供理論指導,同時也可以為燃燒器的優(yōu)化提供理論依據。然而由于回轉窯內傳熱過程復雜、測量手段不完善等原因[4],無法準確得到回轉窯內部的流場、溫度分布信息,故而通過數值模擬的方法分析回轉窯內的流場和溫度場的分布情況具有很大研究價值。
針對回轉窯內部環(huán)境的數值模擬,國內外已經做過很多研究,如Elattar等[5]建立回轉窯的二維模型,研究了回轉窯主要操作對不同燃料下回轉窯內火焰行為的影響;Mujumdar等[6]提出將熟料看作一種偽流體,分別對回轉窯內的空氣和熟料進行三維建模,以偽均勻近似法來模擬固體-固態(tài)反應,建立了水泥回轉窯內的綜合CFD模型;Manju等[7]利用兩相流方法對回轉窯氣動煤注入和燃燒過程進行了建模,研究了回轉窯操作參數對注入煤粉顆粒的分布及其提供熱負荷的影響;Kangwanpongpan等[8]通過優(yōu)化灰度氣體模型的加權和,并結合離散坐標法(Discrete ordinate,DO)模擬了富氧燃燒環(huán)境下的輻射熱傳遞,氣體火焰溫度的預測效果得到了改進;Gaikwad[9]等建立了簡化CFD模型,該模型采用DO模型與加權和灰度氣體方法耦合的輻射熱傳遞模型,以該模型研究了空氣、氧蒸汽等燃燒環(huán)境對溫度和NO濃度分布的影響;廖斌[10]等針對一種低NOX四通道燃燒器,利用FLUENT軟件研究了影響回轉窯內NOX產生的主要因素以及降低NOX排放的方法;馬愛純等[11]利用商業(yè)CFD軟件,對裝備四通道旋流燃燒器的回轉窯進行了數值模擬,分析了不同工況下窯內流場分布、溫度場分布、NOX生成和火焰的特點。Liedmann等[12]建立了固定斜度的3D回轉窯模型,利用商業(yè)CFD軟件研究了煤粉摻雜垃圾衍生燃料(refuse derived fuels,RDF)后,燃燒器操作參數對回轉窯內溫度場的影響。
雖然上述學者所建立的模型能較好地模擬窯內燃料燃燒的特性,然而由于未能綜合考慮回轉窯斜度的影響,所以窯內溫度場與流場的數值模擬結果會和實際情況有較大的差別。
本文在以上研究的基礎上,建立回轉窯內湍流流動、輻射換熱和煤粉燃燒的數學模型。為了研究回轉窯斜度對窯內流場以及煤粉燃燒等的影響,保持燃燒器水平方向不變,對所建立的4個不同斜度的回轉窯-燃燒器流場三維幾何模型進行全結構化網格劃分,采用數值模擬的方法,利用FLUENT軟件研究斜度、旋流風速和軸流風速對窯內流場、窯內空氣溫度場與熟料表面空氣溫度場的影響。
回轉窯內的湍流流動、煤粉燃燒和輻射換熱數值模擬,涉及湍流模型、煤粉燃燒模型以及輻射模型等一系列物理化學過程子模型。為了使仿真結果中回轉窯內部的流場和溫度場更符合實際情況,需要根據回轉窯內的流場特點以及四通道旋流燃燒器的燃燒特點進行子模型的選擇,下面對各子模型進行分析。
在文獻[11]中發(fā)現,標準k-ε模型的計算結果與回轉窯實驗測試結果更吻合,因此氣相湍流模型選擇標準k-ε模型。標準k-ε模型是基于湍流動能k及其耗散率ε的模型傳輸方程的模型。k的模型傳輸方程是從精確方程得出的,而ε的模型傳輸方程是使用物理推理獲得的,與其數學上的精確方程非常相似,標準k-ε模型的傳輸方程為
Gk+Gb-ρε-YM+SK,
(1)
(2)
(3)
式中,μt為湍流(或渦流)粘度,Gk表示平均速度梯度產生的湍流動能,Gb表示由浮力的影響產生的浮力動能,YM表示在可壓縮湍流中的波動膨脹對總體耗散率的貢獻,常量C1ε=1.44、C2ε=1.92、σk=1、σε=1.3。
煤粉燃燒過程大致可以分為熱解和煤焦燃盡兩個階段。對于熱解反應,本文采用目前工程模擬上廣泛使用的雙競爭模型[13],該模型將煤粉的熱解看成兩個熱解反應,這兩個熱解反應同時進行但反應速率不同,表示如下:
(4)
(5)
式中,k1,k2為反應速率常數,服從Arrhenius公式,α1表示煤粉的低溫析出性能,α2為煤粉的高溫析出性能。在低溫段k1>k2,而高溫段k1 煤焦的燃盡過程比煤粉的熱解反應要慢得多,因此煤粉的燃燒時間主要受煤焦燃盡的影響,現在主流的研究觀點認為煤焦的燃燒速率同時受化學反應速率和氧氣擴散速率影響[14],煤焦燃燒速率為 (6) 擴散反應速率常數為 (7) 化學反應速率為 R=C2e-(E/RTP), (8) 式中,A為煤粉顆粒表面積,Pox為煤粉顆粒周圍氧氣分壓,Do為擴散系數,R為化學反應速率常數。 DO輻射模型通過有限數量的離散立體角建立輻射傳遞方程(Radiation transfer equation,RTE),每個離散立體角與全局笛卡爾系統(tǒng)中固定的向量方向相關聯。因為DO模型是通過使用像素化實現的[15],故DO模型最大的特點是空間方向的離散化和相點的處理。在DO模型中,任何位置處的4π角空間的象限被離散為Nθ×Nη個立體角ψ,稱其為控制角,其中角度θ和η分別是極坐標和方位角。 DO模型將RTE視為場方程,因此RTE方程可以寫成 ▽·(I(r,s)s)+(α+σs)I(r,s)= (9) 該模型可以使用灰度模型對非灰度輻射進行建模。光譜強度Iλ(r,s)[27]的RTE為 ▽·(Iλ(r,s)s)+(αλ+σs)Iλ(r,s)=αλn2Ibλ+ (10) 式中,λ為波長,αλ是光譜吸收系數,Ibλ是由普朗克函數給出的黑體強度,并假設散射系數,散射相位函數和折射率n與波長無關。 為了模擬回轉窯實際生產環(huán)境,本文回轉窯幾何模型采用Ф4.0 m×70 m尺寸,但考慮到回轉窯內耐火磚的厚度,實際回轉窯內徑設置為3.6 m。 本仿真所采用的燃燒器為目前廣泛使用的四通道旋流燃燒器,在按實際燃燒器和回轉窯1∶1比例建立幾何模型后,抽取其內部流體域模型,燃燒器內部流體域結構剖面圖如圖1所示。燃燒器4個風道由外向內依次為軸流分道、煤粉風道、旋流風道和中心風道。 熟料平面的位置可以由動態(tài)休止角β和物料中心角ω確定,如圖2所示。物料中心角可以根據填充率φ計算: 2πφ=ω-sinω。 (11) 設置計算工況物料填充率為10%,回轉窯轉速為2 r/min,計算可得ω=93.16°。根據文獻[16],水泥物料的動態(tài)休止角β在25°~33°之間,取β=27°。根據上述計算,建立去除熟料的回轉窯內部流體域模型。 圖1 燃燒器內部流體域幾何模型結構圖 圖2 回轉窯剖面熟料區(qū)域示意圖 實際生產中,回轉窯斜度i一般控制在0.03~0.04之間,故本文建立的4個回轉窯模型的斜度i分別取0、0.03、0.035、0.04。 本文采用工程上廣泛使用的網格劃分軟件ICEM對燃燒器內部流場進行全結構化網格劃分。與非結構化網格相比,結構化網格大大減小了網格數量,因此極大降低了仿真時的計算量,同時增加了計算結果的準確性并提高了計算的收斂性。此外,結構化網格的計算效率也更高,因此,本文網格劃分采用全結構化網格劃分方法。燃燒器網格如圖3所示。 回轉窯及燃燒器入口網格輪廓如圖4所示。由于回轉窯長度過長,故以斜度為0.03的回轉窯為例,將燃燒器和該回轉窯組成的整體網格的剖面圖局部顯示如圖5所示。 圖3 燃燒器網格 圖4 回轉窯及燃燒器入口網格輪廓 圖5 斜度為0.03的回轉窯網格局部剖面圖 本仿真參考實際回轉窯生產工況進行仿真邊界條件設置,設置的標準工況邊界條件參數如表1所示。實際中煤粉顆粒直徑大小為離散分布,假設煤粉顆粒直徑服從Rosin-Rammler分布。為了使仿真更符合實際情況,將煤粉顆粒設置為9種不同直徑的單元,按照工況煤流量總量計算各個煤粉單元所占煤流量大小,煤粉流量分布表格如圖6所示,總煤流量設置為2.46 kg/s,回轉窯標準工況斜度設置為0.03, 氧含量設置為0.2315。燃煤的元素分析和工業(yè)分析如表2所示?;剞D窯壁面設置為絕熱無滑移邊界。 表1 標準工況邊界條件參數表Tab.1 Boundary condition parameter table under standard operating conditions 圖6 不同煤粉顆粒直徑流量分布 表2 燃煤的元素分析和工業(yè)元素分析Tab.2 Ultimate and proximate analyses of the coal burned 為了確保所劃分的網格能夠滿足仿真計算的需求,需要進行網格無關性驗證,以確保仿真結果不會因網格數量過少而影響計算精度。邊界條件均按表1進行設置,劃分3組不同密度的網格模型,對3組網格模型分別進行仿真計算,計算結果如表3所示。 根據表3,綜合考慮計算精度和計算量的要求,選擇網格數量為453 114的網格步長進行不同斜度模型的網格劃分。 表3 網格無關性驗證結果Tab.3 Simulation result of grid independent verification 為了研究不同回轉窯斜度以及不同回轉窯工況對回轉窯內流場、煤粉燃燒等的影響,本文在2.2節(jié)邊界條件的基礎上設置兩組實驗:斜度對流場和煤粉燃燒的影響實驗,將回轉窯運行參數按表2設置,用不同斜度的回轉窯模型進行仿真模擬計算;燃燒器工況對流場和煤粉燃燒的影響實驗,固定回轉窯斜度為0.03,分別改變旋流風道風速和軸流風道風速,其他參數按表2設定,進行仿真模擬,將仿真結果與標準工況的進行對比分析。 基于第一部分建立的數學模型,利用FLUENT軟件對每組實驗進行數值模擬。先開啟1.1節(jié)描述的標準k-ε模型,進行450次冷態(tài)無反應流仿真迭代,以使回轉窯內流場趨于穩(wěn)定,然后開啟1.2節(jié)的煤燃燒模型來點火,進行1 350次熱態(tài)反應流仿真迭代,再開啟1.3節(jié)中的DO輻射模型,進行300次帶輻射的反應流仿真,最后開啟粒子輻射相互作用,進行300次帶粒子輻射相互作用的反應流仿真。 根據上節(jié)邊界條件進行仿真實驗,由于回轉窯長度過長,這里僅顯示窯頭15 m范圍溫度場的仿真結果,窯頭部分垂直方向軸線位置處的溫度場仿真結果如圖7所示。 圖7 不同斜度下窯頭部分垂直方向溫度分布圖 從圖7可以看出,相較于i=0時回轉窯內的火焰,i≠0時回轉窯內火焰存在黑火頭,垂直方向上火焰形狀不規(guī)則情況明顯,然而在i≠0的3組實驗結果中,i=0.035時的火焰形狀明顯比i=0.03及i=0.04的更細長且更規(guī)則。 圖8給出了i=0時軸線方向YZ平面上的流場速度分布,由圖8可以看出,中心風道附近存在回流區(qū),該回流區(qū)是由旋流風的強旋流效應和旋流風與中心風之間的速度差共同引起的。而圖8中的回流區(qū)方向明顯偏下,這是由于熟料區(qū)的存在使回轉窯內的流場在靠近熟料區(qū)附近時風速更高,進而吸引中心回流區(qū)向下偏斜,而斜度的存在會加劇這種偏斜,造成回流區(qū)的減弱甚至消失。 圖9給出了i=0.035時的流場速度分布,在圖9中,可以發(fā)現中心回流區(qū)消失,這正是回轉窯內下方流場流速過高,與上方流場形成較大速差射流引起的。 中心回流區(qū)可以卷吸下游高溫煙氣,有助于煤粉的點火,中心回流區(qū)的減弱或消失使得煤粉的著火點后移,形成黑火頭。黑火頭的存在可以防止燃燒器溫度過高,延長燃燒器壽命,但過長的黑火頭會使煤粉燃燒效率下降、煤耗增加等不利于熟料燒結的因素產生。 圖8i=0時軸線位置垂直方向流場速度分布 圖9i=0.035時軸線位置垂直方向流場速度分布 圖10給出了不同斜度下回轉窯軸線上溫度分布,由此可以看出,i=0.035時回轉窯內的溫度分布要明顯高于i=0.03和i=0.04的。圖11給出了不同斜度下熟料平面軸線上溫度分布,由圖11也可以看出,i=0.035時,熟料平面上溫度高于1 600 K燒成帶范圍要明顯大于其余兩種斜度的。所以i=0.035更適合熟料的燒結。 圖10 不同斜度下回轉窯軸線上溫度分布 圖11給出了不同斜度下熟料平面軸線上溫度分布。 由圖11可以發(fā)現i=0.03和i=0.04的溫度曲線幾乎重合,但圖10中i=0.04的軸線溫度曲線峰值溫度明顯滯后于i=0.03的,由此可以看出,提高回轉窯斜度會使燒成帶區(qū)域靠近窯頭,進而使冷卻區(qū)范圍縮小。圖11中回轉窯出口位置熟料平面溫度升高,這是由出口處高溫煙氣回流造成的。 圖11 不同斜度下熟料平面軸線上溫度分布 3.2.1旋流風速的影響 設置旋流風速分別為70 m/s、80 m/s、100 m/s和110 m/s,其余邊界條件按標準工況進行設置,對比標準工況下的仿真結果。不同旋流風速下的溫度分布對比結果如圖12和圖13所示。 圖12 不同旋流風道風速下回轉窯軸線溫度分布圖 由圖12可以看出,s=80 m/s時的回轉窯軸線峰值溫度最高,但是在圖13中,s=80 m/s時的熟料溫度分布整體并非最高的。造成這種現象的原因主要是:旋流風速的降低使入口切向速度降低,入口切向速度的降低有利于煤粉的聚集,這使得s=80 m/s時火焰有較高的溫度,同時過低的旋流風速大大降低了窯內氣相流場的旋流強度,進而減弱了旋流風對火焰的約束,降低了火焰的穩(wěn)定性,使得火焰熱量不能有效地用于熟料的燒結。圖14為不同旋流風速下窯頭部分水平方向溫度分布圖,圖14(a)和(b)的火焰形狀明顯要比(c)和(d)的更不規(guī)則,也說明了減弱旋流風降低了火焰的穩(wěn)定性。 圖13 不同旋流風速下熟料平面軸線上溫度分布 圖14 不同旋流風速下窯頭部分水平方向溫度分布圖 在圖12中,s=100 m/s時溫度在初始階段的升高速度明顯要低于s=80 m/s和s=90 m/s的,而圖14(e)中火焰的黑火頭明顯要比其余狀態(tài)下的長,這是由于旋流風速的提高增加了入口流場的軸向速度和切向速度,切向速度的增加使煤粉可以與空氣充分混合,同時增強對下游高溫煙氣的卷吸,但是由于斜度的存在使中心回流區(qū)消失,旋流風對下游高溫煙氣的卷吸能力被削弱,而入口流場軸向速度的增加結合旋流風卷吸高溫煙氣能力的減弱共同導致了著火點的后移。 在實際回轉窯運行中,旋流風速的增加,會導致在燃燒器出口區(qū)域的煤粉顆粒擴散,這會使高旋流風速的火焰變粗,而低旋流風速的火焰更細長,對比圖14 中各圖可以發(fā)現,圖14(c)中的火焰明顯要更細長,這與實際運行情況相符,而圖14(a)和(b)中的火焰由于旋流風速過低,火焰穩(wěn)定性不足,造成火焰形狀不規(guī)則。對比s=90 m/s和s=100 m/s在圖12和圖13中的溫度曲線,可以發(fā)現二者十分相似,旋流風速的提高并沒有明顯提升回轉窯內溫度分布水平。結合以上分析,可以認為90 m/s是比較好的旋流風速選擇。 3.2.2軸流風速的影響 分別設置回轉窯軸流風速a=150 m/s、a=160 m/s、a=170 m/s和a=200 m/s,其余邊界條件按標準工況進行設置,對比標準工況下a=180 m/s的仿真結果。軸流風速為150、160、170、180和200 m/s的溫度分布對比結果如圖15和圖16所示。 由圖15可以看出,隨著軸流風速的提高,回轉窯內火焰峰值溫度明顯下降,且峰值溫度隨軸流風速的提高而向窯尾偏移。軸流風速a由150 m/s增加到180 m/s,回轉窯火焰峰值溫度降低了近150 K。在實際回轉窯運行過程中,過高的軸流風速會形成大速差射流,這種強烈的速差射流會在燃燒器軸線位置形成較大的負壓,進而阻礙中心回流區(qū)的形成,使中心回流區(qū)變小,甚至消失。同時,提高軸流風速雖然可以卷吸更多的高溫二次風,但是風速的增加也會使這些高溫二次風與煤粉的混合位置后移,使煤粉的預熱延遲。中心回流區(qū)的減小以及煤粉預熱的延遲共同引起了回轉窯內火焰峰值溫度的后移。并且提高軸流風速,會加速煤粉的運動速度,這將分散煤粉的燃燒,使回轉窯火焰溫度降低。 圖15 不同軸流風速下回轉窯軸線溫度分布圖 圖16 不同軸流風速下熟料平面軸線上溫度分布 圖17為不同軸流風速下窯頭溫度輪廓圖,對比圖17(a)~(e),可以看出,當軸流風速增加到200 m/s時,回轉窯內的火焰嚴重變形,這是由于強烈的軸流風速破壞了旋流風對煤粉的卷吸,使煤粉嚴重擴散,難以聚集燃燒引起的。而對比圖17(a)~(d)可以發(fā)現,隨著軸流風速的增加,黑火頭長度逐漸縮短,這也驗證了提高軸流風會卷吸更多的高溫二次風對煤粉進行預熱的理論。 通過建立的回轉窯內部CFD模型,以及對不同斜度下回轉窯-燃燒器幾何模型進行的全結構化網格劃分,利用FLUENT軟件,對不同斜度及燃燒器工況下的回轉窯進行數值模擬實驗研究,得出以下結論。 1) 回轉窯斜度的存在會使窯內上下區(qū)域流場形成射流速差,且斜度越大,射流速差越大,該射流速差不利于中心回流區(qū)的形成。但斜度的存在同時也有利于燒成帶的擴大。所以回轉窯斜度大小的選擇要適中,而不同燃燒器所造成的射流速差也不同,所以最佳斜度的確定,要依據具體燃燒器的數值模擬實驗來進行選擇。 2) 旋流風速的提高可以卷吸更多下游高溫煙氣,縮短黑火頭,穩(wěn)定火焰燃燒,但是旋流風速的提高對提升窯內溫度場的溫度分布作用不大,而且由于斜度的存在,中心回流區(qū)消失,旋流風對下游高溫煙氣的卷吸作用被弱化,過高的旋流風速不僅不能縮短黑火頭,還會加劇煤粉顆粒擴散。所以旋流風速大小的選擇要適當。 3) 軸流風可以卷吸大量的高溫二次風,對提升燃燒初始階段的溫度有很大作用。提高軸流風速可以縮短黑火頭的長度,但會導致回轉窯火焰溫度降低,且過高的軸流風速會破壞火焰燃燒的穩(wěn)定。 圖17 不同軸流風速下窯頭部分水平方向溫度分布圖1.3 DO輻射模型
2 幾何模型與邊界條件
2.1 回轉窯-燃燒器幾何模型的網格劃分
Fig.1 The structure of the geometric model of the fluid domain inside the burner
Fig.2 A schematic diagram of the clinker area on the rotary kiln profile
Fig.3 The mesh of burner
Fig.4 Rotary kiln and burner entrance mesh profile
Fig.5 Partial profile of rotary kiln mesh with gradient of 0.032.2 邊界條件
Fig.6 Flow distribution of different pulverized coal particle diameter3 仿真結果與分析
3.1 斜度對流場和煤粉燃燒的影響
Fig.7 The temperature distribution in the vertical direction of the kiln portion at different slopes
Fig.8 The flow velocity distribution at the axis position in the vertical direction wheni=0
Fig.9 The flow velocity distribution at the axis position in the vertical direction wheni=0.035
Fig.10 Temperature distribution on the rotary kiln axis at different slopes
Fig.11 Temperature distribution on the plane axis of clinker at different slopes3.2 燃燒器工況對流場和煤粉燃燒的影響
Fig.12 Temperature distribution on rotary kiln axis at different swirling wind speeds
Fig.13 Temperature distribution on clinker plane at different swirling wind speeds
Fig.14 Temperature distribution in horizontal direction at different swirling wind speeds
Fig.15 Temperature distribution of rotary kiln under different axial flow rates
Fig.16 Temperature distribution of clinker plane axis under different axial flow rates4 結論
Fig.17 Temperature distribution in horizontal direction at different axial wind speeds