湍流燃燒反應(yīng)數(shù)值模擬的研究與實(shí)現(xiàn)

孫兆國，劉志勤，劉 濤，劉敏賢，吳穎川

（1.西南科技大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，四川 綿陽621010；2.中國空氣動(dòng)力研究與發(fā)展中心，四川綿陽621000）

0 引 言

燃燒流動(dòng)的數(shù)值模擬是指應(yīng)用計(jì)算機(jī)為工具，將流體力學(xué)、傳熱學(xué)、化學(xué)反應(yīng)動(dòng)力學(xué)和數(shù)值計(jì)算方法相結(jié)合所得到的求解化學(xué)流體力學(xué)基本方程的理論和方法。數(shù)值模擬的基本思想是：把空間與時(shí)間坐標(biāo)中連續(xù)的物理量的場（如速度場，溫度場，濃度場等），用一系列有限個(gè)離散點(diǎn)上的值的集合來代替，通過一定的原則建立起這些離散點(diǎn)上變量值之間關(guān)系的代數(shù)方程，并根據(jù)所建立起來的代數(shù)方程求解這些變量的近似值。數(shù)值模擬的研究方法對(duì)燃燒室性能評(píng)估有一定的參考價(jià)值。

本文的研究重點(diǎn)是基于RANS［1］求解的定長可壓縮湍流燃燒流動(dòng)的數(shù)值模擬。湍流是工程技術(shù)領(lǐng)域中常見的流動(dòng)形式，在湍流流動(dòng)中流體的各種物理參數(shù)，如速度、壓力、溫度等都隨時(shí)間與空間發(fā)生隨機(jī)變化。流體做湍流燃燒時(shí)涉及到流體流動(dòng)、對(duì)流換熱以及化學(xué)反應(yīng)等過程，數(shù)值模擬比較繁瑣。在燃燒流動(dòng)數(shù)值模擬中需要借助計(jì)算流體力學(xué)平臺(tái)。當(dāng)前計(jì)算流體力學(xué)數(shù)值模擬領(lǐng)域多是商業(yè)軟件，如Fluent、CFX等。由于商業(yè)機(jī)密等原因，其源代碼難以開放，當(dāng)計(jì)算結(jié)果不理想時(shí)，針對(duì)特定問題的判斷和解決較困難。本文燃燒流動(dòng)數(shù)值模擬采用C＋＋語言開發(fā)的開源計(jì)算流體力學(xué)軟件OpenFOAM［2］，這對(duì)于開發(fā)新的針對(duì)特定問題的解算器有重要的幫助作用。

在OpenFOAM平臺(tái)上做燃燒數(shù)值模擬的研究尚處于初級(jí)階段，Aalborg University的 Christian Andersen等［3］使用OpenFOAM中的dieselFoam解算器對(duì)基于Burner Flow Reactor（BFR）的柴油噴射進(jìn)行數(shù)值模擬，中國科學(xué)技術(shù)大學(xué)的陸陽、武文等［4－6］改進(jìn)了OpenFOAM平臺(tái)自帶的基于RANS求解的非定常燃燒流動(dòng)解算器reactingFoam對(duì)湍流燃燒進(jìn)行數(shù)值模擬，這些都是對(duì)OpenFOAM中自帶的解算器進(jìn)行改進(jìn)且是對(duì)非定?？蓧嚎s燃燒過程做數(shù)值模擬。

本文根據(jù)定?？蓧嚎s湍流燃燒流動(dòng)機(jī)理和化學(xué)動(dòng)力反應(yīng)模型庫Cantera設(shè)計(jì)出湍流燃燒流動(dòng)的解算器，然后使用該解算器對(duì)Sydney鈍體駐定火焰［7－8］進(jìn)行數(shù)值模擬。模擬得到了燃燒流動(dòng)組分濃度分布圖和溫度曲線變化圖等。通過將模擬的溫度變化與Sydney鈍體駐定火焰的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比驗(yàn)證表明：模擬效果滿足燃燒組分變化的研究要求，流場溫度較好的符合Sydney鈍體駐定火焰HM1的實(shí)驗(yàn)數(shù)據(jù)，這驗(yàn)證了解算器的合理性和可行性。

1 數(shù)學(xué)模型

在湍流燃燒反應(yīng)中，反應(yīng)物湍流流動(dòng)過程和化學(xué)反應(yīng)過程相互關(guān)聯(lián)又相互影響，湍流流動(dòng)組分濃度及溫度的脈動(dòng)可以強(qiáng)化組分的混合與傳熱，因而影響時(shí)平均化學(xué)反應(yīng)速率，化學(xué)反應(yīng)放熱過程迅速放熱而引起密度變化，同時(shí)使流體輸運(yùn)系數(shù)變化，進(jìn)而影響著湍流流動(dòng)過程。這種相互關(guān)聯(lián)和相互影響使得湍流燃燒反應(yīng)極為復(fù)雜。湍流燃燒流動(dòng)數(shù)值模擬的基本思想是分別獨(dú)立描述湍流流動(dòng)和化學(xué)反應(yīng)過程，然后考慮湍流流動(dòng)和化學(xué)反應(yīng)相互作用，基于這個(gè)思想，在數(shù)值模擬過程中使用的燃燒模型包括湍流模型，湍流燃燒相互作用模型和化學(xué)反應(yīng)模型。

湍流模型主要包括層流（laminar），直接數(shù)值模擬（direct numerical simulation，DNS），大渦模擬（large eddy simulation，LES）和雷諾時(shí)均數(shù)值模擬（RANS）［9］。雷諾時(shí)均數(shù)值模擬（RANS）是將瞬態(tài)N—S方程的瞬時(shí)量分解為時(shí)均值和脈動(dòng)值之和（雷諾分解），再取時(shí)間平均，得到雷諾時(shí)均方程，然后利用某些模擬假設(shè)，將方程中的高階的未知關(guān)聯(lián)項(xiàng)用低階項(xiàng)或時(shí)均量來表達(dá)，從而使雷諾（reynolds）時(shí)均方程封閉。

本文中使用的模型湍流是基于RANS的k－ε雙方程模型和Cantera化學(xué)動(dòng)力反應(yīng)模型庫。

1.1 湍流模型

在反應(yīng)物進(jìn)行湍流混合時(shí)，使用k－ε雙方程模型。kε雙方程模型對(duì)于一般的流動(dòng)情況，如平壁邊界層、無力平面射流、管流等流動(dòng)雙方程模型能給出相當(dāng)滿意的計(jì)算結(jié)果，而且計(jì)算工作量小。k－ε雙方程模型［10］，模型如下

在上述方程中，Gk表示由于平均速度梯度引起的湍動(dòng)能產(chǎn)生，Gb表示用于浮力影響引起的湍動(dòng)能產(chǎn)生，YM表示可壓速湍流脈動(dòng)膨脹對(duì)總的耗散率的影響。

k－ε雙方程模型的湍動(dòng)能k和耗散率ε方程求解式為式中：Ux′，Uy′，Uz′——x，y，z這3個(gè)方向的速度矢量，l——特征長度，Cμ設(shè)定——0.09。湍流粘性系數(shù)求解式為

1.2 化學(xué)反應(yīng)模型

化學(xué)反應(yīng)模型使用Cantera化學(xué)動(dòng)力反應(yīng)模型庫［11］。Cantera是跨平臺(tái)的面向?qū)ο蟮拈_源軟件，主要應(yīng)用于化學(xué)反應(yīng)，熱力學(xué)和傳輸過程等問題的求解。Cantera化學(xué)反應(yīng)模型是基于美國氣體研究所和加州大學(xué)的研究成果GRIMech3.0詳細(xì)化學(xué)反應(yīng)機(jī)理［12］。GRI－Mech3.0詳細(xì)化學(xué)反應(yīng)機(jī)理總共涉及到53種組分和325個(gè)基元反應(yīng)，其中包括C1反應(yīng)、C2反應(yīng)、C3反應(yīng)和NOX的生成機(jī)理。

對(duì)于燃料CmHn，其和氧氣的總包反應(yīng)可以寫成［13］

對(duì)于CO，其和氧氣的總包反應(yīng)可以寫成

2 仿真平臺(tái)的建立

2.1 OpenFOAM平臺(tái)

本文解算器設(shè)計(jì)使用近幾年流行起來的計(jì)算流體力學(xué)軟件OpenFOAM作為平臺(tái)。OpenFOAM的全稱是Open Field Operation and Manipulation，它的核心是使用 C＋＋編寫的面向?qū)ο蟮挠?jì)算流體力學(xué)類庫，有較高效率的偏微分方程求解模塊。OpenFOAM是一款開源軟件，其源代碼是對(duì)外公開的，這就可以方便用戶編寫解算器。Open－FOAM軟件可以模擬復(fù)雜的流體流動(dòng)、化學(xué)反應(yīng)、傳熱傳質(zhì)等問題，還可以進(jìn)行結(jié)構(gòu)動(dòng)力學(xué)分析、電磁場分析等的數(shù)值仿真。

從程序?qū)崿F(xiàn)功能的角度來看，OpenFOAM軟件包括核心解算器、前處理和后處理3大模塊。前處理主要用網(wǎng)格設(shè)計(jì)的工具來進(jìn)行，如OpenFOAM本身自帶的BlockMesh工具，也可以使用其他商業(yè)網(wǎng)格設(shè)計(jì)軟件，如GridGen等。在本文中，燃燒室網(wǎng)格設(shè)計(jì)采用GridGen，然后轉(zhuǎn)換成OpenFOAM可以識(shí)別的網(wǎng)格文件。解算器主要是求解計(jì)算流體力學(xué)問題常用的一些解算器，OpenFOAM的解算器支持的流動(dòng)模型有層流（laminar），基于雷諾時(shí)均湍流模型（RANS），大渦模擬（LES）以及直接數(shù)值模擬（DNS）等。后處理主要是將計(jì)算結(jié)果轉(zhuǎn)化成可視化圖像，如云圖，燃燒等值面圖等。OpenFOAM的組件構(gòu)成如圖1所示。

圖1 OpenFOAM 的組件構(gòu)成［14－15］

2.2 解算器

解算器位于/opt/application/solvers/combustion/turbulentFoam目錄下。解算器包括Make目錄，頭文件和主程序3部分。其中Make目錄包括options文件和files文件，options文件定義編譯選項(xiàng)，用于指定編譯用到的頭文件位置及其動(dòng)態(tài)庫，files文件用于指定當(dāng)前要編譯的文件。頭文件是主程序編譯過程中調(diào)用用到的程序，解算器主程序是turbulentFoam.C。解算器結(jié)構(gòu)如圖2所示。

圖2 解算器結(jié)構(gòu)

3 算例

3.1 物理模型

計(jì)算區(qū)域如圖3所示，鈍體半徑為50mm，燃料入口半徑為3.6mm。鈍體網(wǎng)格結(jié)構(gòu)簡圖如圖4所示，使用結(jié)構(gòu)化網(wǎng)格離散計(jì)算區(qū)域，網(wǎng)格數(shù)為5580，并在軸向劃分5°，以便在OpenFOAM中模擬軸對(duì)稱問題。

網(wǎng)格有上下兩個(gè)進(jìn)氣口，上口為伴隨流，注入成分為空氣，為了簡化，使用3組分空氣構(gòu)成，分別為N2，O2和Ar，伴隨流的溫度為300K。射流入口注入燃料煤氣，煤氣中包含CH4，CO，H2，CO2，H2O和N2。反應(yīng)工況如表1所示。反應(yīng)開始時(shí)，鈍體內(nèi)充滿空氣。燃燒反應(yīng)的邊界條件見表1。

表1 燃燒反應(yīng)的邊界條件和初始條件

3.2 離散方法

對(duì)控制方程采用有限體積法進(jìn)行離散，對(duì)控制方程中的瞬態(tài)項(xiàng)采用LU－SGS隱式迭代法求解，對(duì)流項(xiàng)差分采用二階迎風(fēng)的高斯差分格式，擴(kuò)散項(xiàng)差分采用二階線性修正的高斯差分格式。由于是定?？蓧嚎s流動(dòng)計(jì)算，本文采用SIMPLE（semi－implicit method for pressure－linked equations）［16］半隱式壓力校正方程算法，SIMPLE算法是流體流場數(shù)值模擬中最重要的模擬方法。SIMPLE算法的求解流程如圖5所示。

圖5 SIMPLE算法求解流程

3.3 求解過程及分析

計(jì)算平臺(tái)使用基于MPI的并行計(jì)算系統(tǒng)［17］。并行計(jì)算原理即在計(jì)算前將燃燒室網(wǎng)格剖分為4個(gè)子區(qū)域分配給各CPU內(nèi)核，把子區(qū)域的初始相關(guān)信息分別裝載入各子區(qū)域?qū)?yīng)的CPU內(nèi)存，每一個(gè)CPU內(nèi)核分別啟動(dòng)一個(gè)計(jì)算進(jìn)程，這樣同時(shí)有四個(gè)進(jìn)程并行運(yùn)算，可以加快計(jì)算時(shí)間。

在數(shù)值模擬階段，將Coalgas和Air分別從fuelInlet和airInlet以118m/s和40m/s的速度注入燃燒室。反應(yīng)物注入后在燃燒室內(nèi)進(jìn)行燃燒反應(yīng)，經(jīng)過21600s的計(jì)算，燃燒反應(yīng)穩(wěn)定。圖6，圖7和圖8分別為燃燒反應(yīng)穩(wěn)定后產(chǎn)物二氧化碳，水，速度分布圖。從圖中可以看出燃燒反應(yīng)從噴口到狂漲段都富含燃燒反應(yīng)產(chǎn)物二氧化碳和水，并在輸出口富集。分布云圖較好的反映了反應(yīng)流場在反應(yīng)穩(wěn)定時(shí)的一些基本特征。

圖8 速度分布（m/s）

燃燒反應(yīng)的參數(shù)溫度是與實(shí)驗(yàn)結(jié)果對(duì)比的重要參數(shù)。圖9是燃燒反應(yīng)穩(wěn)定時(shí)燃燒室截面中取一條平行于X軸的端點(diǎn)坐標(biāo)為（0，0）和（0，0.4）的直線并根據(jù)此直線上坐標(biāo)對(duì)應(yīng)的溫度得到的溫度分布曲線圖，從曲線圖可以得出最高溫度超過2200K。

Sydney鈍體駐定火焰是國際湍流擴(kuò)散火焰測量與計(jì)算研討會(huì)（TNF workshop）的標(biāo)準(zhǔn)系列火焰之一，常用于與仿真結(jié)果做對(duì)比。實(shí)驗(yàn)溫度分布與Sydney鈍體駐定火焰的HM1實(shí)驗(yàn)結(jié)果基本吻合，尤其是燃燒的最高溫度等變化情況的計(jì)算結(jié)果很好。

圖9 燃燒反應(yīng)穩(wěn)定時(shí)溫度分布（T/K）

從圖6，圖7和圖8可以看出，火焰上游沒有出現(xiàn)由燃料射流擴(kuò)散而形成的回流區(qū)，這使得模擬出的效果偏向?qū)恿魅紵?，因此在?shù)值模擬時(shí)，k－ε雙方程模型中k和ε的設(shè)置還得進(jìn)一步探討，以便得到更好的模擬效果。另外，由于在數(shù)值模擬時(shí)考慮涉及詳細(xì)化學(xué)反應(yīng)機(jī)理，求解的規(guī)模巨大，為了加快計(jì)算時(shí)間，將網(wǎng)格規(guī)模設(shè)計(jì)較小，這使得有些區(qū)域網(wǎng)格不夠密集，尤其是火焰上游區(qū)域。這對(duì)數(shù)值模擬精度也有一定的影響。

4 結(jié)束語

本文在開源計(jì)算流體力學(xué)平臺(tái)OpenFOAM和Cantera化學(xué)動(dòng)力反應(yīng)模型庫的基礎(chǔ)上，設(shè)計(jì)出基于定常可壓縮的湍流燃燒流動(dòng)反應(yīng)解算器，并通過該解算器對(duì)Sydney鈍體駐定湍流擴(kuò)散火焰HM1進(jìn)行數(shù)值模擬，在Coalas/Air詳細(xì)反應(yīng)機(jī)理的條件下得到的計(jì)算結(jié)果并得出了燃燒流場中各組分濃度和溫度變化曲線圖，仿真效果較好的符合燃燒組分變化的研究要求，尤其是溫度曲線變化效果比較理想。模擬效果的具體情況基本驗(yàn)證了湍流燃燒流動(dòng)解算器設(shè)計(jì)的合理性，同時(shí)，對(duì)OpenFOAM平臺(tái)和Cantera反應(yīng)庫的研究對(duì)于將開源軟件應(yīng)用在燃燒流動(dòng)數(shù)值模擬領(lǐng)域有重要的開拓意義。

［1］ZHANG Zhaoshun，CUI Guixiang，XU Chunxiao.Theory and modeling of turbulence ［M］.Beijing：Tsinghua University Press，2005（in Chinese）.［張兆順，崔桂香，許春曉.湍流理論與模擬 ［M］.北京：清華大學(xué)出版社，2005.］

［2］Hrvoje Jasak.OpenFOAM：Open source CFD in research and industry ［J］.International Journal of Naval Architecture and Ocean Engineering，，2009，1（2）：89－94.

［3］Niels E L Nielsen.Numerical investigation of a BFR using OpenFOAM ［R］.Fluids and Combustion Engineering，2008.

［4］LU Yang.A study on flamelet model in combustion simulation［D］.Hefei：University of Science and Technology of China，2009：93－117（in Chinese）.［陸陽.燃燒計(jì)算中火焰面模型的研究［D］.合肥：中國科學(xué)技術(shù)大學(xué)博士論文，2009：93－117.］

［5］WU Wen，HUANG Wei，ZHAO Pinghui，et al.Numerical simulation of turbulent CH4/Air of jet diffusion ［C］.Hefei：The Chinese Society of Engineering Thermophysics Conference，2009（in Chinese）.［武文，黃威，趙平輝，等.湍流CH4/Air射流擴(kuò)散燃燒數(shù)值模擬 ［C］.合肥：中國工程熱物理學(xué)會(huì)學(xué)術(shù)會(huì)議論文，2009.］

［6］WU Wen，HUANG Wei，ZHAO Pinghui，et al.Study of RANS based flamelet/progress variable turbulent combustion model ［J］.Journal of University of Science and Technology of China，2010，40（10）：1016－1022（in Chinese）.［武文，黃威，趙平輝，等.基于RANS求解的火焰面/反應(yīng)進(jìn)度變量湍流燃燒模型研究 ［J］.中國科學(xué)技術(shù)大學(xué)學(xué)報(bào)，2010，40（10）：1016－1022.］

［7］Dally B B，Masri A R，Barlow R S，et al.Instantaneous and mean compositional steucture of bluff body stabilized non－premixed flames［J］.Combustion F1ame，1998，114（1－2）：119－148.

［8］HUANG Qing，ZHU Minming，YE Taohong，et al.PDF simulation of bluff－body stabilized turbulent non－premixed flame［J］.Chinese Journal of Comutational Physics，2008，25（6）：733－743（in Chinese）.［黃慶，朱旻明，葉桃紅，等.PDF方法模擬鈍體駐定的湍流擴(kuò)散火焰 ［J］.計(jì)算物理，2008，25（6）：733－743.］

［9］LU Yu.Kinetic reduction for combustion chemistry and investigation of methane turbulent jet flame using direct numerical simulation［D］.Hangzhou：Zhejiang University，2011（in Chinese）.［呂鈺.燃燒化學(xué)機(jī)理簡化及甲烷湍流射流火焰的直接數(shù)值模擬研究［D］.杭州：浙江大學(xué)碩士論文，2011.］

［10］SONG Baojun.On tow－equation turbulence model for complex turbulent flow by using DI approximation ［C］.Kyoto，Japan：The 9th International Conference of Turbulent Shear Flow，1993.

［11］Moffat H K.CADS：Cantera aerosol dynamics simulator ［M］.Albuquerque：Sandia National Laboratories，2007：7－9.

［12］Michael Frenklach.GRI－Mech30 ［EB/OL］.http：//www.me.berkeley.edu/gri＿mech/version30/text30.html，2011.

［13］WANG Hui，HOU Lingyun.Numerical simulation of hydrocarbon gas combustion and supersonic shear flow ［J］.Journal of Aerospace Power，2007，22（4）：559－564（in Chinese）.［王慧，侯凌云.碳?xì)淙細(xì)獬曀偌羟辛鲃?dòng)燃燒數(shù)值模擬［J］.航空動(dòng)力學(xué)報(bào)，2007，22（4）：559－564.］

［14］OpenFOAM Company.OpenFOAM user guide ［M］.Berkshire：OpenCFD Company，2010.

［15］Mattijs janssens andy heather sergio ferraris graham macpherson helene blanchnet henry weller CFD toolbox－user guide［M］.Berkshire：OpenCFD Company，2010.

［16］WANG Fujun.The analysis of computational fluid dynamics－principles and applications of CFD software ［M］.Beijing：Tsinghua University Press，2005：74－75（in Chinese）.［王福軍.計(jì)算流體動(dòng)力學(xué)分析－CFD軟件原理與應(yīng)用 ［M］.北京：清華大學(xué)出版社，2005：74－75.］

［17］Edgar Gabriel.OpenMPI：Open source high performance computing ［EB/OL］.http：//www.open－mpi.org/，2011.