不同半徑圓柱誘導(dǎo)CH4/空氣預(yù)混燃燒的計算研究

余亦澤, 徐勝利, 張夢萍, 張竹鶴

(1. 中國科學(xué)技術(shù)大學(xué) 數(shù)學(xué)學(xué)院, 安徽 合肥 230026; 2. 清華大學(xué) 航空航天學(xué)院, 北京 100084)

0 引 言

為捕捉化學(xué)非平衡湍流流動(燃燒等)的精細結(jié)構(gòu)，Boger等開展了直接模擬(DNS)[1-2]和大渦模擬(LES)[3-4]等計算研究。LES在亞網(wǎng)格尺度采用和非定常雷諾平均(URANS)湍流計算相似的?；枷耄档土擞嬎憬Y(jié)果對湍流模型的依賴性，即湍流模型只用于亞網(wǎng)格尺度。已有湍流模型主要來自不可壓縮流動思想，湍流模型及其常數(shù)和工況相關(guān)，這給可壓縮湍流流場計算帶來不確定因素。DNS計算不需要湍流模型，但要求網(wǎng)格和時間尺度均小于對應(yīng)的Kolmogorov尺度。和URANS相比，LES和DNS計算都采用細網(wǎng)格，計算量增大。

URANS高精度計算是開展LES和DNS計算的基礎(chǔ)，URANS采用高精度數(shù)值格式可減小計算量。Balsara和Shu的LES計算結(jié)果表明[5]：要得到近似相同的計算結(jié)果，和9階格式相比，2階格式在各坐標軸方向的網(wǎng)格數(shù)約為前者4倍。當網(wǎng)格數(shù)相同，9階格式計算量是2階格式的3倍。對三維問題，要得到近似相同的計算結(jié)果，2階格式計算量是9階格式44/3(≈85)倍。

可壓縮化學(xué)非平衡流動高精度差分計算遇到的問題是：(1)如何處理復(fù)雜幾何域?(2)高精度數(shù)值格式如何滿足保自由流要求?貼體變換是處理復(fù)雜幾何域常用方法之一。在變換后的曲線坐標系中，如果格式不能滿足自由流條件，會導(dǎo)致計算出現(xiàn)數(shù)值誤差，抹平細微湍流結(jié)構(gòu)或者氣動聲學(xué)波等微弱的物理擾動。同時，不保自由流條件產(chǎn)生的計算誤差還會導(dǎo)致高階格式數(shù)值不穩(wěn)定[6-7]。

針對兩類高精度中心型緊致格式，在三維廣義坐標系中，Visbal和Gaitonde[7]仔細研究了Thomas和Lombard[8]提出的守恒形式網(wǎng)格導(dǎo)數(shù)差分誤差，發(fā)現(xiàn)網(wǎng)格導(dǎo)數(shù)項和對流項采用相同差分格式，緊致格式才能保持自由流條件。但標準高階WENO格式數(shù)值通量是直接重構(gòu)分解后的通量，通量包含了網(wǎng)格導(dǎo)數(shù)。當非線性重構(gòu)通量時，很難將Visbal和Gaitonde想法[7]應(yīng)用到標準WENO格式，從而無法滿足流表面的網(wǎng)格導(dǎo)數(shù)表面守恒(surface conservation law)和體守恒(volume conservation law)。因此，標準WENO格式在曲線網(wǎng)格的多維流場計算難以保持精確的自由流解[5]。JIANG等人[9]提出基于對解進行WENO插值的另一種WENO格式，隨后將Visbal和Gaitonde[7]想法應(yīng)用到新WENO格式[10]。針對多種不規(guī)則和運動網(wǎng)格以及超聲速圓柱繞流，JIANG等研究新WENO格式保自由流和保渦(vortex-preserving)性質(zhì)。結(jié)果表明：不論是固定網(wǎng)格還是動網(wǎng)格，相對于標準WENO 格式，新WENO格式(簡稱保自由流WENO格式)都可保證自由流條件和渦條件。

在本文條件下，預(yù)混燃燒主要受化學(xué)動力學(xué)機理影響。對圓柱誘導(dǎo)預(yù)混燃燒問題，圓柱誘導(dǎo)的激波和火焰駐點距離可檢驗不同基元反應(yīng)模型的點火延遲[11-14]。采用9 組分/12步基元反應(yīng)模型，Soetrisno[15]計算了CH4/空氣的激波和火焰耦合流場。采用由33組分/149步基元反應(yīng)模型得到的簡化模型(19組分/52步反應(yīng))，Yungster和Robinowitz[11]計算了CH4/空氣解耦的爆轟波流場。采用14組分/19步基元反應(yīng)模型，Toshimitsu[12]較好預(yù)測了CH4/空氣預(yù)混燃燒的點火延遲。沖壓加速器[16]和斜爆轟發(fā)動機[17]推進原理也是基于超聲速來流預(yù)混燃燒。

為開展超聲速化學(xué)反應(yīng)流場高精度計算，針對圓柱誘導(dǎo)CH4/空氣超聲速來流預(yù)混燃燒問題，本文在曲線坐標系應(yīng)用保自由流WENO格式進行求解，研究圓柱不同半徑對計算結(jié)果的影響。

1 控制方程和定解條件

1.1 曲線坐標系下的控制方程

本文計算采用Euler方程形式為:

(1)

其中，坐標變換取t=τ，x=x(ξ,η,ζ)，y=y(ξ,η,ζ)，z=z(ξ,η,ζ)。

Q=J-1Q，Q=[ρ1,…,ρns,ρu,ρv,ρw,E]T，

E(Q)=[ρ1u,…,ρnsu,ρu2+p,ρuv,ρuw,u(E+p)]T，

F(Q)=[ρ1v,…,ρnsv,ρuv,ρv2+p,ρvw,v(E+p)]T，

G(Q)=[ρ1v,…,ρnsv,ρuw,ρvw,ρw2+p,w(E+p)]T，

(2)

其中，Ru是普適氣體常數(shù)，T是混合物溫度，Wi是第i組元摩爾質(zhì)量。

1.2 熱力學(xué)參數(shù)

對ns組元量熱完全氣體構(gòu)成的混合物，第i組元定壓比熱cpi和hi由多項式擬合得到，即:

(3)

(4)

其中，a1i,a2i,…,a6i取自JANAF表[18]。

1.3 基元反應(yīng)模型

包含ns組元、NR反應(yīng)步的基元反應(yīng)統(tǒng)一表示為:

(5)

由質(zhì)量作用定律得到第i組元的質(zhì)量生成速率為:

(6)

其中，kf和kb分別為正向和逆向反應(yīng)速率常數(shù)，通常為Arrhenius定律形式，具體為:

(7)

其中，Aj、βj、Eaj分別為第j步反應(yīng)的指前因子、溫度指數(shù)和活化能。平衡常數(shù)Kcj由Gibbs自由焓確定[19]，逆反應(yīng)速率常數(shù)kbj由kfi和Kcj確定，即:

(8)

本文采用的基元反應(yīng)模型[9]見表1所列。

2 數(shù)值方法簡介

2.1 保自由流WENO格式

采用保自由流5階WENO格式求解方程(1)，方程(1)半離散格式為:

(9)

表1 CH4/空氣基元反應(yīng)模型Table 1 Chemistry model for methane air mixture

保自由流WENO和標準WENO格式的差異在于數(shù)值通量重構(gòu)方法不同，前者數(shù)值通量由解和解的導(dǎo)數(shù)插值后得到[10]。以ξ方向為例，有

(10)

2.2 三階TVD Runge-Kutta方法

采用三階TVD Runge-Kutta方法離散方程(1)時間導(dǎo)數(shù)項。對:

(11)

具體形式為[20]:

(12)

其中，Δt為時間步長，Lh為空間離散算子。

3 計算結(jié)果分析和討論

3.1 數(shù)值方法驗證

從格式精度和保自由流性質(zhì)兩方面驗證本文選用的數(shù)值方法。參考文獻[9]，取如下變換:

x=ξ+0.01sinξcosη

y=η-0.02cosξsinη

對方程(1)，初值取為:

ρ=1.0+0.2sin(x+2y)

p=1.0,u=0.5,v=-0.5

在(x，y)采用以2π為周期的周期性邊界條件。表2給出t=2密度計算值的L1和L∞誤差以及對應(yīng)的階。從表2看出，本文采用的WENO格式具有5階精度。

表2 ρ的L1和L∞誤差以及對應(yīng)的階(t=2)Table 2 L1 errors and L∞ errors and related order of ρ at t=2

本文選用Jiang[9]的算例驗證新WENO格式保自由流條件。取如下計算網(wǎng)格:

xi,j,k=-2+(i-1)Δx0+

0.2sin [π(j-1)Δy0]sin[π(j-1)Δy0]

yi,j,k=-2+(j-1)Δy0+

0.2sin[π(k-1)Δz0]sin[π(i-1)Δx0]

zi,j,k=-2+(k-1)Δx0+

0.2sin[π(i-1)Δx0]sin[π(j-1)Δy0]

其中，i=1,2,…,Imax，j=1,2,…,Jmax，k=1,2,…,Kmax，Imax=Jmax=Kmax=21，Δx0=Δy0=Δz0=0.2，Δt取為0.05，最大計算時間取為10，ρ和聲速設(shè)為常數(shù)，自由流沿x方向且馬赫數(shù)為0.5。表3給出了標準WENO和保自由流WENO格式計算得到v和w的L2誤差。表3表明：保自由流WENO格式的L2誤差接近機器誤差，標準WENO的L2誤差約為10-3，這驗證了保自由流WENO的保自由流性質(zhì)，也校驗了本文計算程序。要說明的是：保自由流WENO格式更多精度分析和算例驗證見文獻[9, 10]。

表3 v和w的L2誤差Table 3 L2 errors of v and w

3.2 圓柱誘導(dǎo)CH4/空氣超聲速預(yù)混燃燒

本文選擇和非平衡化學(xué)反應(yīng)流計算的的標模算例結(jié)果[11]進行對照。算例1給出圓柱誘導(dǎo)超聲速預(yù)混燃燒示意圖(圖1(a))，圓柱半徑R為0.5 mm，圓心位于x=-1。計算參數(shù)和邊界條件簡述如下：

自左至右來流靜溫T∞、靜壓p∞和速度U∞分別為295 K、51 kPa和2330 m/s，對應(yīng)的馬赫數(shù)Ma∞為6.61，預(yù)混氣為化學(xué)計量比CH4/空氣混合物，計算網(wǎng)格取40×30，見圖1(b)。入口邊界AD指定超聲速來流參數(shù)，出口邊界AB和CD取外推條件，固壁BC采用絕熱滑移條件。

圖1 圓柱誘導(dǎo)超聲速預(yù)混燃燒和計算網(wǎng)格示意圖Fig.1 Schematic of a cylinder induced combustion and grid distribution

圖2和圖3分別給出壓力和溫度等值線與云圖分布。圖2清楚地顯示圓柱上游產(chǎn)生的弓形激波位置。圖3也顯示了激波和火焰陣面位置，對應(yīng)著兩道溫度間斷。自左至右的第一個溫度間斷是弓形激波產(chǎn)生的波后氣流溫度上升，第二個溫度間斷對應(yīng)化學(xué)反應(yīng)(燃燒)放熱導(dǎo)致的溫度上升。圖2和圖3均表明：激波和火焰陣面的厚度很薄，特別是火焰陣面，這表明預(yù)混燃燒的化學(xué)反應(yīng)速率很大。仔細觀察還可看出，沿著圖3火焰陣面遠離對稱線，火焰厚度隨著距離增大而略有增大，相當于火焰厚度略微變厚。原因是：弓形激波后的當?shù)貧饬鱉a數(shù)逐漸變大且溫度降低。根據(jù)方程(7)和(8)，當?shù)厝紵齥fj和kbj略有減小，化學(xué)反應(yīng)速率略有降低，從而導(dǎo)致預(yù)混火焰陣面略微變厚。

圖2 壓力分布等值線和云圖(算例1)Fig.2 Pressure contours and snapshots (case1)

圖3 溫度分布等值線和云圖(算例1)Fig.3 Temperature contours and snapshots (case1)

圖4給出壓力和溫度沿過駐點水平線的分布。圖4(a)顯示壓力只有一次上升，這表明火焰面兩側(cè)壓力近似相等，近似為等壓燃燒。圖4(b)顯示了溫度兩次上升，這和圖3也是對應(yīng)的。為了定量考察計算結(jié)果的網(wǎng)格無關(guān)性，圖4還給出了網(wǎng)格加密一倍(80×60)后的壓力和溫度計算結(jié)果。圖4表明：當計算網(wǎng)格加密一倍(80×60)，粗細網(wǎng)格計算得到的壓力和溫度沿駐點線分布近似重合。這表明本文粗網(wǎng)格計算得到的計算結(jié)果也是可靠的。圖4(c)還給出典型反應(yīng)物(CH4、O2)和生成物(H2O、CO2)質(zhì)量百分數(shù)沿過駐點線分布。圖4(c)表明：化學(xué)反應(yīng)區(qū)約占據(jù)9 網(wǎng)格點。顯然，本文采用網(wǎng)格量是可以分辨化學(xué)反應(yīng)區(qū)的。圖5給出了Yungster[11](網(wǎng)格數(shù)91×91)和本文計算得到的激波和火焰陣面位置。這些數(shù)據(jù)分別從Yungster[11]溫度等值線和本文圖3提取的。圖5表明：盡管本文計算采用網(wǎng)格數(shù)較少(40×30)，但是，兩者的激波形狀近似重合，火焰面位置也大致重合，只是在遠離過駐點線的重合度略差，可能原因是本文和Yungster采用不同基元反應(yīng)模型。

在達到同樣收斂解前提下，采用5階WENO格式、較少網(wǎng)格數(shù)(40×30)和2階TVD格式、較多網(wǎng)格數(shù)(91×91)，所用CPU運算時間分別約為6 h和10 h，這表明：在本文條件下，采用高精度格式和較少網(wǎng)格數(shù)可提高計算效率。

(a) 壓力

(b) 溫度

(c) 組元質(zhì)量百分數(shù)

圖5 本文和Yungster[11]激波和火焰面位置(算例1)Fig.5 Shapes of a shock and a flame front in contrast to those in Ref[11](case 1)

3.3 不同圓柱半徑的影響

以算例1為參考，本節(jié)研究不同半徑圓柱誘導(dǎo)CH4/空氣超聲速預(yù)混燃燒的影響。來流參數(shù)和邊界條件同算例1，圖6給出了R分別為1.5 mm和3.5 mm對應(yīng)的溫度等值線與云圖分布。

圖6表明：和圖3類似，不同半徑圓柱上游流場均產(chǎn)生兩道溫度間斷，分別對應(yīng)弓形激波(自左至右第一個間斷)和火焰面(自左至右第二個間斷)的氣流溫度上升。R分別為1.5 mm和3.5 mm的圓柱，對應(yīng)的激波駐點距離分別為0.18mm、0.59mm和1.82mm，駐點距離隨著圓柱半徑增大而增大。即使不考慮來流預(yù)混氣的化學(xué)反應(yīng)，不同半徑圓柱對上游來流的阻礙作用不同，所產(chǎn)生的弓形激波形狀或強度也不相同，即弓形激波后的氣流溫度不相同。根據(jù)方程(7、8)，不同溫度對應(yīng)不同的kfj和kbj，從而影響當?shù)鼗瘜W(xué)反應(yīng)速率，造成預(yù)混燃燒流場溫度和壓力也不相同。

R分別為0.5 mm、1.5 mm和3.5 mm圓柱，對應(yīng)的火焰面駐點距離分別為0.05 mm、0.51 mm和1.75 mm。和激波駐點距離類似，火焰面駐點距離隨圓柱半徑增大而增大。

(a) R=1.5 mm(算例2)

(b) R=3.5 mm(算例3)

從圖6還可看出，弓形火焰面(化學(xué)反應(yīng)區(qū))幾何形狀和激波形狀類似。盡管來流條件相同，但是，不同半徑圓柱對應(yīng)的激波和火焰面之間距離是不同的，兩者沿駐點線的距離是最小的。原因是：來流經(jīng)過弓形激波壓縮，弓形激波波后氣流沿激波陣面法線的Ma數(shù)是不同的，偏離駐點線越遠，沿激波陣面法線Ma數(shù)就越小。根據(jù)激波關(guān)系，對應(yīng)的氣流溫度和壓力也就越低。由基元反應(yīng)質(zhì)量作用定律知，組元密度時間變化率和溫度、分壓力(對應(yīng)摩爾濃度)是直接相關(guān)的。因此，弓形激波陣面不同位置對應(yīng)的預(yù)混氣流反應(yīng)速率也是不同的。沿激波陣面離開駐點線，弓形激波后氣流化學(xué)反應(yīng)由當?shù)貋喡曀俎D(zhuǎn)變?yōu)槌曀?。不同的波后氣流壓力和溫度對?yīng)著不同的反應(yīng)速率，宏觀上表現(xiàn)為圖6弓形激波和弓形火焰面之間的不同距離，該距離通常稱為點火延時距離(ignition induction length)。和當?shù)貧饬魉俣汝P(guān)聯(lián)，由點火延時距離可得到點火延時(ignition delay)。測量過駐點線的圓柱上游激波和火焰之間距離，換算為基元反應(yīng)動力學(xué)機理對應(yīng)的點火延時，表4給出本文不同半徑圓柱對應(yīng)的點火延時計算值。為和實驗數(shù)據(jù)對比，根據(jù)激波管甲烷點火延時數(shù)據(jù)整理的擬合公式[22]:

(16)

其中：τ為點火延時，μs；T為溫度，K；[O2]和[CH4]分別為O2和CH4摩爾濃度，mol/cm3。從表4看出，擬合公式[22]給出的點火延時和本文計算值較接近。點火延時隨圓柱半徑增大而增大。

表4 不同半徑圓柱對應(yīng)的點火延時計算值和擬合值Table 4 Computed and fitted ignition delay at different cylindrical radiuses

對應(yīng)圖2、圖3和圖6，圖7給出壓強和溫度沿過駐點線分布。作為比較，圖7還給出R為0.5 mm圓柱無化學(xué)反應(yīng)流場激波駐點距離。為清楚地顯示不同算例結(jié)果的差別，圖7只顯示圓柱表面左側(cè)附近的區(qū)域。圖7(a)表明：不同半徑圓柱對應(yīng)的激波駐點位置是不同的，圓柱半徑越大，對應(yīng)的駐點距離也越大。圖7(b)表明：來流經(jīng)過弓形激波和火焰面到達圓柱表面，不同半徑圓柱對應(yīng)不同的火焰面位置，但燃燒產(chǎn)物在圓柱表面的最大溫度相同，這表明圓柱上游的預(yù)混氣燃燒是充分的。對R為3.5 mm圓柱，激波和火焰面沿過駐點線的距離最短。這表明：圓柱半徑越大，對預(yù)混燃燒流場的影響也越大。原因是：弓形激波下游氣流是局部亞聲速流動，圓柱對流場擾動(阻擋作用)向上游傳播，從而影響激波和火焰駐點距離。這和圖6與圖7(a)給出的結(jié)論也是一致的。

(a) 壓強

(b) 溫度

圖8給出流場Ma數(shù)云圖分布。圖8表明：弓形激波和火焰面下游氣流Ma數(shù)低于來流。經(jīng)過激波后，超聲速氣流速度降低，但溫度和壓力升高，氣流動能轉(zhuǎn)變?yōu)閮?nèi)能，有利于可燃預(yù)混氣燃燒。圖8還給出過駐點線兩側(cè)附近的亞聲速流場分布，不同半徑圓柱對應(yīng)的亞聲速區(qū)域大小也不相同。其中，R=0.5 mm圓柱對應(yīng)的亞聲速區(qū)域最小，R=3.5 mm圓柱對應(yīng)的亞聲速區(qū)域最大。在弓形激波波后流場，Mach數(shù)隨著離開駐點線距離增大而逐漸增大，即由亞聲速區(qū)過渡到超聲速區(qū)，分別對應(yīng)著亞聲速燃燒和超聲速燃燒。圖8還可解釋不同半徑圓柱在駐點線附近流場存在的壁面、火焰面和激波之間相互干擾。圓柱壁面對來流擾動通過亞聲速區(qū)域向上游傳播，直至影響到激波和火焰陣面。

對R為1.5 mm和3.5 mm圓柱預(yù)混燃燒流場，圖9給出本文和Yungster[11]計算得到的激波和火焰面位置。和圖5類似，圖9的激波和火焰位置數(shù)據(jù)也是從溫度等值線提取的。圖9表明：當R=1.5 mm，Yungster[11]與本文得到的激波和火焰面位置基本重合，只是遠離過駐點線的火焰面位置略有差異。當R=3.5 mm，Yungster和本文計算得到的激波位置大致重合，但在遠離駐點線的區(qū)域，激波和火焰面位置存在差異，其中，火焰面位置的差異尤其明顯。主要原因是Yungster采用了不同的基元反應(yīng)模型和網(wǎng)格分布。兩者在過駐點線附近區(qū)域的激波和火焰面位置完全重合。本文和Yungster得到的弓形激波無量綱駐點距離均為xshock/d=0.3。其中，xshock和d分別是激波駐點距離和圓柱直徑。這和文獻[21]激光全息得到xshock/d=0.29值非常接近。

圖8 馬赫數(shù)分布云圖Fig.8 Mach number counter map

(a)R=1.5 mm

(b)R=3.5 mm

圖9 本文和Yungster[11]激波和火焰面位置
(算例2和算例3)
Fig.9 Shapes of a shock and a flame front in contrast
to those in Ref[11](case2 and case3)

4 結(jié) 論

1) 在本文條件下，采用保自由流5階WENO格式計算圓柱誘導(dǎo)CH4/空氣預(yù)混燃燒流場，消除了曲線網(wǎng)格下不保自由流產(chǎn)生的計算誤差，表現(xiàn)了較好計算穩(wěn)定性，這和文獻[9]結(jié)論是相同的。

2) 本文得到的激波和火焰面形狀與文獻計算結(jié)果相符，特別是過駐點線的激波和火焰位置完全相同。和文獻[11]二階TVD格式相比，采用保自由流5階WENO格式和近似四分之一網(wǎng)格數(shù)，得到近似相同的計算結(jié)果，提高了計算效率。

3) 圓柱不同半徑影響著超聲速來流預(yù)混燃燒的激波和火焰面駐點距離。圓柱半徑越大，對應(yīng)的火焰面和激波之間距離也越小，相應(yīng)的點火延時也越小。圓柱表面對來流阻礙作用或擾動，通過圓柱上游的亞聲速區(qū)傳播，改變著過駐點線附近的弓形激波和火焰面形狀和駐點距離。