湍流燃燒反應數值模擬的研究與實現

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

（1.西南科技大學 計算機科學與技術學院，四川 綿陽621010；2.中國空氣動力研究與發展中心，四川綿陽621000）

0 引 言

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

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

在OpenFOAM平臺上做燃燒數值模擬的研究尚處于初級階段，Aalborg University的 Christian Andersen等［3］使用OpenFOAM中的dieselFoam解算器對基于Burner Flow Reactor（BFR）的柴油噴射進行數值模擬，中國科學技術大學的陸陽、武文等［4－6］改進了OpenFOAM平臺自帶的基于RANS求解的非定常燃燒流動解算器reactingFoam對湍流燃燒進行數值模擬，這些都是對OpenFOAM中自帶的解算器進行改進且是對非定常可壓縮燃燒過程做數值模擬。

本文根據定常可壓縮湍流燃燒流動機理和化學動力反應模型庫Cantera設計出湍流燃燒流動的解算器，然后使用該解算器對Sydney鈍體駐定火焰［7－8］進行數值模擬。模擬得到了燃燒流動組分濃度分布圖和溫度曲線變化圖等。通過將模擬的溫度變化與Sydney鈍體駐定火焰的實驗數據進行對比驗證表明：模擬效果滿足燃燒組分變化的研究要求，流場溫度較好的符合Sydney鈍體駐定火焰HM1的實驗數據，這驗證了解算器的合理性和可行性。

1 數學模型

在湍流燃燒反應中，反應物湍流流動過程和化學反應過程相互關聯又相互影響，湍流流動組分濃度及溫度的脈動可以強化組分的混合與傳熱，因而影響時平均化學反應速率，化學反應放熱過程迅速放熱而引起密度變化，同時使流體輸運系數變化，進而影響著湍流流動過程。這種相互關聯和相互影響使得湍流燃燒反應極為復雜。湍流燃燒流動數值模擬的基本思想是分別獨立描述湍流流動和化學反應過程，然后考慮湍流流動和化學反應相互作用，基于這個思想，在數值模擬過程中使用的燃燒模型包括湍流模型，湍流燃燒相互作用模型和化學反應模型。

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

本文中使用的模型湍流是基于RANS的k－ε雙方程模型和Cantera化學動力反應模型庫。

1.1 湍流模型

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

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

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

1.2 化學反應模型

化學反應模型使用Cantera化學動力反應模型庫［11］。Cantera是跨平臺的面向對象的開源軟件，主要應用于化學反應，熱力學和傳輸過程等問題的求解。Cantera化學反應模型是基于美國氣體研究所和加州大學的研究成果GRIMech3.0詳細化學反應機理［12］。GRI－Mech3.0詳細化學反應機理總共涉及到53種組分和325個基元反應，其中包括C1反應、C2反應、C3反應和NOX的生成機理。

對于燃料CmHn，其和氧氣的總包反應可以寫成［13］

對于CO，其和氧氣的總包反應可以寫成

2 仿真平臺的建立

2.1 OpenFOAM平臺

本文解算器設計使用近幾年流行起來的計算流體力學軟件OpenFOAM作為平臺。OpenFOAM的全稱是Open Field Operation and Manipulation，它的核心是使用 C＋＋編寫的面向對象的計算流體力學類庫，有較高效率的偏微分方程求解模塊。OpenFOAM是一款開源軟件，其源代碼是對外公開的，這就可以方便用戶編寫解算器。Open－FOAM軟件可以模擬復雜的流體流動、化學反應、傳熱傳質等問題，還可以進行結構動力學分析、電磁場分析等的數值仿真。

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

圖1 OpenFOAM 的組件構成［14－15］

2.2 解算器

解算器位于/opt/application/solvers/combustion/turbulentFoam目錄下。解算器包括Make目錄，頭文件和主程序3部分。其中Make目錄包括options文件和files文件，options文件定義編譯選項，用于指定編譯用到的頭文件位置及其動態庫，files文件用于指定當前要編譯的文件。頭文件是主程序編譯過程中調用用到的程序，解算器主程序是turbulentFoam.C。解算器結構如圖2所示。

圖2 解算器結構

3 算例

3.1 物理模型

計算區域如圖3所示，鈍體半徑為50mm，燃料入口半徑為3.6mm。鈍體網格結構簡圖如圖4所示，使用結構化網格離散計算區域，網格數為5580，并在軸向劃分5°，以便在OpenFOAM中模擬軸對稱問題。

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

表1 燃燒反應的邊界條件和初始條件

3.2 離散方法

對控制方程采用有限體積法進行離散，對控制方程中的瞬態項采用LU－SGS隱式迭代法求解，對流項差分采用二階迎風的高斯差分格式，擴散項差分采用二階線性修正的高斯差分格式。由于是定常可壓縮流動計算，本文采用SIMPLE（semi－implicit method for pressure－linked equations）［16］半隱式壓力校正方程算法，SIMPLE算法是流體流場數值模擬中最重要的模擬方法。SIMPLE算法的求解流程如圖5所示。

圖5 SIMPLE算法求解流程

3.3 求解過程及分析

計算平臺使用基于MPI的并行計算系統［17］。并行計算原理即在計算前將燃燒室網格剖分為4個子區域分配給各CPU內核，把子區域的初始相關信息分別裝載入各子區域對應的CPU內存，每一個CPU內核分別啟動一個計算進程，這樣同時有四個進程并行運算，可以加快計算時間。

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

圖8 速度分布（m/s）

燃燒反應的參數溫度是與實驗結果對比的重要參數。圖9是燃燒反應穩定時燃燒室截面中取一條平行于X軸的端點坐標為（0，0）和（0，0.4）的直線并根據此直線上坐標對應的溫度得到的溫度分布曲線圖，從曲線圖可以得出最高溫度超過2200K。

Sydney鈍體駐定火焰是國際湍流擴散火焰測量與計算研討會（TNF workshop）的標準系列火焰之一，常用于與仿真結果做對比。實驗溫度分布與Sydney鈍體駐定火焰的HM1實驗結果基本吻合，尤其是燃燒的最高溫度等變化情況的計算結果很好。

圖9 燃燒反應穩定時溫度分布（T/K）

從圖6，圖7和圖8可以看出，火焰上游沒有出現由燃料射流擴散而形成的回流區，這使得模擬出的效果偏向層流燃燒，因此在數值模擬時，k－ε雙方程模型中k和ε的設置還得進一步探討，以便得到更好的模擬效果。另外，由于在數值模擬時考慮涉及詳細化學反應機理，求解的規模巨大，為了加快計算時間，將網格規模設計較小，這使得有些區域網格不夠密集，尤其是火焰上游區域。這對數值模擬精度也有一定的影響。

4 結束語

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

［1］ZHANG Zhaoshun，CUI Guixiang，XU Chunxiao.Theory and modeling of turbulence ［M］.Beijing：Tsinghua University Press，2005（in Chinese）.［張兆順，崔桂香，許春曉.湍流理論與模擬 ［M］.北京：清華大學出版社，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）.［陸陽.燃燒計算中火焰面模型的研究［D］.合肥：中國科學技術大學博士論文，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射流擴散燃燒數值模擬 ［C］.合肥：中國工程熱物理學會學術會議論文，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求解的火焰面/反應進度變量湍流燃燒模型研究 ［J］.中國科學技術大學學報，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方法模擬鈍體駐定的湍流擴散火焰 ［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）.［呂鈺.燃燒化學機理簡化及甲烷湍流射流火焰的直接數值模擬研究［D］.杭州：浙江大學碩士論文，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）.［王慧，侯凌云.碳氫燃氣超聲速剪切流動燃燒數值模擬［J］.航空動力學報，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）.［王福軍.計算流體動力學分析－CFD軟件原理與應用 ［M］.北京：清華大學出版社，2005：74－75.］

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