不同半徑圓柱誘導CH4/空氣預混燃燒的計算研究

2020-03-11 08:06:42余亦澤徐勝利張夢萍張竹鶴

空氣動力學學報 2020年1期

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

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

0 引言

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

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

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

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

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

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

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是普適氣體常數，T是混合物溫度，Wi是第i組元摩爾質量。

1.2 熱力學參數

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

(3)

(4)

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

1.3 基元反應模型

包含ns組元、NR反應步的基元反應統一表示為:

(5)

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

(6)

其中，kf和kb分別為正向和逆向反應速率常數，通常為Arrhenius定律形式，具體為:

(7)

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

(8)

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

2 數值方法簡介

2.1 保自由流WENO格式

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

(9)

表1 CH4/空氣基元反應模型Table 1 Chemistry model for methane air mixture

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

(10)

2.2 三階TVD Runge-Kutta方法

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

(11)

具體形式為[20]:

(12)

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

3 計算結果分析和討論

3.1 數值方法驗證

從格式精度和保自由流性質兩方面驗證本文選用的數值方法。參考文獻[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∞誤差以及對應的階。從表2看出，本文采用的WENO格式具有5階精度。

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

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

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，ρ和聲速設為常數，自由流沿x方向且馬赫數為0.5。表3給出了標準WENO和保自由流WENO格式計算得到v和w的L2誤差。表3表明：保自由流WENO格式的L2誤差接近機器誤差，標準WENO的L2誤差約為10-3，這驗證了保自由流WENO的保自由流性質，也校驗了本文計算程序。要說明的是：保自由流WENO格式更多精度分析和算例驗證見文獻[9, 10]。

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

3.2 圓柱誘導CH4/空氣超聲速預混燃燒

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

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

圖1 圓柱誘導超聲速預混燃燒和計算網格示意圖Fig.1 Schematic of a cylinder induced combustion and grid distribution

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

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

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

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

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

(a) 壓力

(b) 溫度

圖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為參考，本節研究不同半徑圓柱誘導CH4/空氣超聲速預混燃燒的影響。來流參數和邊界條件同算例1，圖6給出了R分別為1.5 mm和3.5 mm對應的溫度等值線與云圖分布。

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

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

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

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

從圖6還可看出，弓形火焰面(化學反應區)幾何形狀和激波形狀類似。盡管來流條件相同，但是，不同半徑圓柱對應的激波和火焰面之間距離是不同的，兩者沿駐點線的距離是最小的。原因是：來流經過弓形激波壓縮，弓形激波波后氣流沿激波陣面法線的Ma數是不同的，偏離駐點線越遠，沿激波陣面法線Ma數就越小。根據激波關系，對應的氣流溫度和壓力也就越低。由基元反應質量作用定律知，組元密度時間變化率和溫度、分壓力(對應摩爾濃度)是直接相關的。因此，弓形激波陣面不同位置對應的預混氣流反應速率也是不同的。沿激波陣面離開駐點線，弓形激波后氣流化學反應由當地亞聲速轉變為超聲速。不同的波后氣流壓力和溫度對應著不同的反應速率，宏觀上表現為圖6弓形激波和弓形火焰面之間的不同距離，該距離通常稱為點火延時距離(ignition induction length)。和當地氣流速度關聯，由點火延時距離可得到點火延時(ignition delay)。測量過駐點線的圓柱上游激波和火焰之間距離，換算為基元反應動力學機理對應的點火延時，表4給出本文不同半徑圓柱對應的點火延時計算值。為和實驗數據對比，根據激波管甲烷點火延時數據整理的擬合公式[22]:

(16)

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

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

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

(a) 壓強

(b) 溫度

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

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