基于OpenFOAM的燃燒數(shù)值模擬研究

王耀彬+蔡德琛

摘要：現(xiàn)有商業(yè)計算軟件對燃燒室內(nèi)的燃燒流場進行計算時存在精度不高、非開源、不可擴展等問題，難以滿足工程設(shè)計需要。在前期完成冷態(tài)霧化模擬基礎(chǔ)上，開發(fā)了基于開源計算平臺OpenFOAM的燃油霧化燃燒求解系統(tǒng)，構(gòu)建大渦求解器，對加州大學(xué)伯克利分校測量的液霧兩相燃燒實驗進行了數(shù)值模擬。通過計算得出液滴溫度維持在350K左右，表明模擬效果較符合燃燒流場溫度變化，系統(tǒng)能較好地展示湍流燃燒流動的反應(yīng)特性，符合工程研究需要。

關(guān)鍵詞：燃燒數(shù)值模擬；并行計算；OpenFOAM

DOIDOI：10.11907/rjdk.171423

中圖分類號：TP319

文獻標(biāo)識碼：A 文章編號：1672-7800（2017）006-0117-03

0 引言

化工、航空航天等工程應(yīng)用領(lǐng)域存在液霧-空氣兩相流動現(xiàn)象。發(fā)動機燃燒室內(nèi)多為復(fù)雜湍流流動，傳統(tǒng)的液體燃料直接與空氣混合往往導(dǎo)致燃燒不充分，耗能高、效用低。預(yù)先將燃料經(jīng)燃油噴嘴霧化，然后與空氣混合燃燒，理論上可明顯提高燃燒效率，減少發(fā)動機燃油殘渣殘留[1-2]。傳統(tǒng)的液霧燃燒模擬多使用雷諾時均（RANS）或者直接模擬（DNS）。雷諾平均方法特點是將流動方程組統(tǒng)計平均后建模，只能提供湍流的平均信息。大渦模擬是預(yù)先對湍流劃分尺度，大尺度渦采用直接求解，小尺度湍流脈動建立亞網(wǎng)格燃燒模型，這樣做的優(yōu)點是對空間分辨率要求小，計算速度較快，并且能獲得比RANS更多的湍流信息。但現(xiàn)有商業(yè)計算軟件對燃燒流場計算存在精度不足、非開源、不可擴展問題，不能很好地滿足工程設(shè)計的需要[3-4]。

針對上述問題，本文采用大渦方法對N-S方程進行過濾封閉處理，結(jié)合相應(yīng)燃燒湍流模型，采用PIMPLE算法計算湍流動方程，設(shè)計開發(fā)了基于開源流體力學(xué)平臺OpenFOAM的大渦求解器，并通過構(gòu)建算例完成液體燃料液霧兩相燃燒模擬。對比分析實驗驗證了該燃燒流場模擬系統(tǒng)的可靠性與穩(wěn)定性。

1 OpenFOAM簡介

OpenFOAM中文譯為開源的場運算和處理軟件，遵守GNU通用公共許可證，用戶可以依據(jù)自身需求對軟件進行針對性拓展。其前身FOAM是Hrvoje Jasak所寫，后來將其更名為OpenFOAM。代碼是一款基于C++編寫的面向?qū)ο蟮拈_源流體力學(xué)軟件包，采用有限容積方法求解，在模擬湍流動、物理化學(xué)反應(yīng)、復(fù)雜流體流動等方面效果非常好。它不僅提供了許多編譯好的模型庫、輔助工具和求解器等，而且自帶偏微分求解模塊以及網(wǎng)格生成工具，可以在PolyMesh目錄里生成多面體網(wǎng)格，支持多種流動模型。

2 燃燒流場求解模塊設(shè)計

2.1 物理模型

加州大學(xué)伯克利分校在燃燒實驗領(lǐng)域一直居于領(lǐng)先地位，實驗過程中采用基于激光成像原理的設(shè)備記錄瞬間數(shù)據(jù)，所以實驗數(shù)據(jù)全面、準(zhǔn)確，可以為數(shù)值模擬出的數(shù)據(jù)進行對比分析[5-6]，本文基于OpenFOAM平臺對液霧燃燒實驗進行數(shù)值模擬工作。

燃燒室物理模型如圖1所示，此圖為剖面圖，裝置高200mm，以乙醇進口管為中心左右軸對稱，燃燒室內(nèi)橫置一層巢狀網(wǎng)格，網(wǎng)格與底部平板填充大量的玻璃小珠，乙醇燃油從底部中心管以小液滴的形式噴入燃燒室，霧化蒸發(fā)與兩側(cè)的空氣伴隨流充分混合，網(wǎng)格底端中心點為Inlet，頂端中心點為outlet。

反應(yīng)開始時，燃油經(jīng)Inlet噴入燃燒室霧化，然后與兩側(cè)空氣混合燃燒。實驗以乙醇為燃料，使用空氣作為伴隨流，得到穩(wěn)定態(tài)的氣體溫度、液體溫度及軸向速度等數(shù)據(jù)。實驗噴嘴直徑為10mm，燃料壓力在1.4-2.6bar之間變化，液體流量在0.39-0.54g/s之間變化，空氣伴隨流的速度則設(shè)定在0-0.64m/s之間。

2.2 幾何模型

系統(tǒng)根據(jù)流體力學(xué)數(shù)值求解流程，在OpenFOAM平臺下構(gòu)建計算區(qū)域幾何模型 [9]。

本文使用Gambit軟件來劃分網(wǎng)格文件。生成的網(wǎng)格文件默認為fluent（CFD商業(yè)軟件）格式，Gambit默認的長度為毫米，OpenFOAM中默認為米。為了將Gambit軟件繪制的網(wǎng)格文件導(dǎo)入OpenFOAM平臺下，必須轉(zhuǎn)化兩者單位。將網(wǎng)格文件轉(zhuǎn)化為OpenFOAM可操作的結(jié)構(gòu)化文件，導(dǎo)入到OpenFOAM里，使用paraView工具查看網(wǎng)格形狀，如圖2所示。三維計算區(qū)域網(wǎng)格半徑為100mm，高度為160mm，網(wǎng)格總數(shù)為8.5萬個。

2.3 燃燒流場數(shù)值算法

本文采用PIMPLE算法作為燃燒流場數(shù)值算法。PIMPLE算法實際是根據(jù)SIMPLE和PISO算法各自特點，取長補短融合的一種新型算法，算法步驟如圖3所示。其將每個時間步長看成穩(wěn)態(tài)流動，循環(huán)求解，特點如下：

（1）PIMPLE算法和瞬態(tài)PISO算法相似，不同之處在于前者在時間步長內(nèi)增加了速度壓力耦合循環(huán)過程。

（2）PIMPLE算法在時間步進內(nèi)使用亞松弛技術(shù)解決兩個時間段物理量變化大的問題，其實質(zhì)是將每個時間步看成穩(wěn)態(tài)流動，這一點采用了SIMPLE算法思想。

2.4 求解模塊設(shè)計

燃燒室中反應(yīng)物的流動過程大都是湍流流動，所用的物理量都是時間和空間的隨機變量，其瞬時量可以用流體力學(xué)經(jīng)典的N-S（納維-斯托克斯）方程描述，而燃燒模型采用部分攪拌模型（Partially Stirred Reactor Combustion Model，簡稱PaSR模型）。

大渦模擬方法求解燃燒流場復(fù)雜湍流動問題，關(guān)鍵在于對湍流動尺度進行劃分，采用直接求解法計算大尺度湍流動，采用亞網(wǎng)格模型模擬小尺度渦。大渦模擬中常用的亞網(wǎng)格模型有Smagorinsky模型、尺度相似模型、LRR差分應(yīng)力模型等，其中Smagorinsky模型是目前應(yīng)用較廣泛的亞網(wǎng)格模型，其優(yōu)點是可通過添加渦黏系數(shù)模塊高效地對控制方程進行數(shù)值計算，OpenFOAM對渦黏系數(shù)提供支持[7-8]。

本文采用Smagorinsky模型作為亞網(wǎng)格模型，在第一次動量預(yù)測時選取湍流模型。求解器由Make目錄、頭文件和主程序myLesFoam.C組成。頭文件包括一些方程求解文件、場文件等，作用是聲明變量，從文件中讀入初值等。myLesFoam.C為求解器主程序，由PIMPLE算法程序及湍流動求解程序組成。

3 燃油霧化算例設(shè)計

液體燃料燃燒過程大致包含霧化、蒸發(fā)、摻混及燃燒4個部分，燃油霧化為第一步，步驟如下：①燃油通過噴油嘴流出形成液滴柱；②燃油射流初始噴出產(chǎn)生湍流，在周圍氣體作用下，液滴產(chǎn)生波動、最終分離出液滴小碎片；③在液滴表面張力作用下，液滴碎片壓縮成球形小液滴；④在氣動力作用下，小液滴進一步破碎。

模擬采用KHRT（Kelvin-Helmholtz-Rayleigh-Taylor）模型作為液滴破碎模型，該破碎模型與TAB模型是當(dāng)前使用最為廣泛的破碎模型。KHRT模型分為兩次破碎，KH破碎是由于氣體與液體之間的速度差異引起的，RT破碎則考慮到液滴-氣體界面的加速度引起的液滴表面波增長。

算例分為time目錄、constant目錄和system目錄3部分。time目錄下對燃燒反應(yīng)的邊界條件及初始條件進行設(shè)定，包括O2、k（湍動能產(chǎn)生）、ε（耗散率）以及速度、壓力和溫度的設(shè)定；Chemkin目錄為OpenFOAM，耦合了化學(xué)反應(yīng)動力學(xué)數(shù)據(jù)庫Chemkin的化學(xué)反應(yīng)機理，在constant目錄下的thermophysicalProperties調(diào)用化學(xué)反應(yīng)機理文件。turbulenceProperties設(shè)定使用的模擬方式為LES，LESProperties設(shè)定模型。polyMesh為網(wǎng)格文件，其中的核心文件為blockMeshDict，燃燒室網(wǎng)格結(jié)構(gòu)是為解讀此文件而構(gòu)造的；system目錄包含3個文件：fvSchemes、fvSolution 和 controlDict文件。具體部署步驟：①選擇代數(shù)方程求解器。在fvSolution文件里對壓力采用預(yù)條件共軛梯度法，對速度則采用預(yù)條件雙共軛梯度法，二者區(qū)別在于前者是求解對稱矩陣，后者用來求解反對稱矩陣；

②設(shè)置k-ε湍流模型。采用預(yù)條件雙共軛梯度法求解k，ε；

③選擇方程求解方法后，進行離散格式設(shè)置。首先采用歐拉格式離散，然后對梯度進行離散，本文采用高斯線性插值方法；散度離散中對流項采用高斯理論，TVD格式中的限制型線性差分，擴散項采用迎風(fēng)差分格式；

④設(shè)置拉普拉斯項離散，采用高斯方法，二階線性守恒格式。

設(shè)置OpenFOAM時間控制參數(shù)字典controlDict文件步驟如下：

①在算例開始求解時，程序要控制計算的初始時間，設(shè)定值為startTime，含義是從startTime指定的時間開始計算；

②設(shè)置計算程序結(jié)束時間。值為endTime，從endTime指定時間結(jié)束計算；

③設(shè)置計算時間步長。為了精確求解，本文采用5x10-6s作為一個時間步長。在數(shù)據(jù)輸出過程中會產(chǎn)生時間文件夾，這里按照0.001s物理時間寫入一次直到計算結(jié)束；

④設(shè)置寫入過程不覆蓋，數(shù)據(jù)的輸出精度為6位，計算程序運行時參數(shù)允許改變，并將改變立即反饋給程序。

4 實驗結(jié)果

對算例求解收斂后得到燃油濃度分布云圖，從圖中可以看出，燃油噴入燃燒室后，首先形成一錐形空心油膜，隨之迅速變薄破碎成小液珠，在湍流作用下，液珠蒸發(fā)并與燃料出口空氣充分混合，形成新的燃燒場，如圖4所示。

對圖5和圖6得到的模擬溫度數(shù)據(jù)與實驗數(shù)據(jù)進行對比：系統(tǒng)模擬的液體溫度與伯克利分校液霧試驗數(shù)據(jù)相差較小，反應(yīng)剛開始時要加快反應(yīng)速度，提高燃燒室溫度；當(dāng)燃料破碎形成液滴后，液體溫度維持在350K左右，模擬得到的液滴溫度能夠較好地與實驗數(shù)據(jù)吻合。

5 結(jié)語

本文在開源流體力學(xué)平臺OpenFOAM上設(shè)計大渦求解器，通過構(gòu)建算例完成液體燃料液霧兩相燃燒模擬，得到的燃油濃度分布云圖表明液體與空氣混合較好，火焰的溫度場分布總體趨勢符合實驗數(shù)據(jù)，驗證了本燃燒流場模擬系統(tǒng)的可靠性與穩(wěn)定性。

參考文獻：

[1]王海鷗，羅坤，白云，等. 不同湍流強度下液霧燃燒的直接數(shù)值模擬[J]. 燃燒科學(xué)與技術(shù)，2014 （3）：222-226.

[2]周力行，李科，王方，等. 液霧燃燒大渦模擬的應(yīng)用研究進展[J]. 化工學(xué)報，2010（11）：2769-2775.

[3]HAKAN NILSSON.Evaluation of OpenFOAM for CFD of turbulent flow in water turbines[C]. Proceedings of the 23rd IAHR Symposium in Yokohama， 2006.

[4]李科，周力行. 乙醇-空氣液霧兩相流動和燃燒的大渦模擬[J]. 工程熱物理學(xué)報，2011（6）：965-968.

[5]陳靖，朱旻明，劉明侯，等. 乙醇-空氣稀液霧值班火焰的大渦模擬研究[J]. 推進技術(shù)，2015（2）：276-284.

[6]楊建山，羅坤，邵長孝，等. 氣液兩相旋流燃燒的時均統(tǒng)計特性[J]. 工程熱物理學(xué)報，2015（2）：335-337.

[7]陸陽. 燃燒計算中火焰面模型的研究[D]. 合肥：中國科學(xué)技術(shù)大學(xué)， 2009.

[8]武文，黃威，趙平輝，等.基于RANS求解的火焰面/反應(yīng)進度變量湍流燃燒模型研究[J]. 中國科學(xué)技術(shù)大學(xué)學(xué)報，2010， 40（10）：1016-1022.

[9]王華僑. 結(jié)構(gòu)有限元分析中的網(wǎng)格劃分技術(shù)及其應(yīng)用實例[J]. CAD/CAM與制造業(yè)信息化，2005（1）：42-47.

（責(zé)任編輯：杜能鋼）