關(guān)于混合氣體輸運系數(shù)計算方案的影響研究

王曉棟,盧德勇

(中國空氣動力研究與發(fā)展中心 計算空氣動力研究所,綿陽 621000)

0 引 言

吸氣式超燃沖壓發(fā)動機允許來流以超聲速狀態(tài)進入燃燒室,是一種新型吸氣式動力裝置,在未來航空、航天領(lǐng)域中具有廣闊的應(yīng)用前景。在過去的幾十年中,人們對超燃沖壓發(fā)動機技術(shù)進行了大量試驗、理論及數(shù)值計算研究[1]。對超燃沖壓發(fā)動機設(shè)計研究來講,理解超燃沖壓發(fā)動機燃燒室內(nèi)部的流場結(jié)構(gòu)、流動特性是非常重要的。單純依靠風(fēng)洞試驗難以給出詳細的流場細節(jié),而與風(fēng)洞試驗緊密結(jié)合的數(shù)值模擬在提供流場細節(jié)的同時,也有助于減少風(fēng)洞試驗內(nèi)容,提高風(fēng)洞試驗效率。超燃沖壓發(fā)動機內(nèi)部流場數(shù)值模擬增加了對組分方程的求解過程,這使得流場模擬的計算量急劇增加。為了降低計算量,提高流場求解效率,除了進一步發(fā)展高效算法、采用并行技術(shù)外,通常還需要通過簡約組分來降低組分守恒方程的規(guī)模。此外,在組分方程的求解過程中,為了獲得動量、組分以及能量的輸運強度,頻繁計算組分及混合氣體的輸運系數(shù),這是造成計算量劇增的主要因素之一。在混合氣體輸運系數(shù)計算方面,混合氣體分子黏性系數(shù)由各組分的黏性系數(shù)表達,混合氣體導(dǎo)熱系數(shù)和擴散系數(shù)通常有兩種計算方案[2]:1)假設(shè)普朗特數(shù)、施密特數(shù)為常數(shù),由混合氣體的分子黏性系數(shù)關(guān)聯(lián)相應(yīng)的導(dǎo)熱系數(shù)和擴散系數(shù);2)在計算組分的導(dǎo)熱系數(shù)和擴散系數(shù)基礎(chǔ)上,利用半經(jīng)驗公式計算混合氣體的導(dǎo)熱系數(shù)和擴散系數(shù)。受單個組分的基本輸運系數(shù)數(shù)據(jù)匱乏和計算機能力的限制,早期多采用第一類方法計算混合氣體導(dǎo)熱系數(shù)和擴散系數(shù)[3],現(xiàn)在針對空氣的流場計算仍主要采用該方案[4]。目前,組分基本輸運系數(shù)有了很大的擴充,計算機能力也有了大幅度的提高,第二類計算方案已成為主流。

第二類計算方案雖然提高了混合氣體輸運系數(shù)計算的準(zhǔn)確度,但是卻增加了流場的計算量,由此降低了流場求解效率。在大規(guī)模、多組分復(fù)雜流場計算時,為了提高流場的求解效率,根據(jù)流場實際靈活選擇導(dǎo)熱系數(shù)和擴散系數(shù)計算方案,這就需要開展混合氣體輸運系數(shù)求解方案對流場計算影響的對比研究。當(dāng)前,國內(nèi)外超聲速燃燒流場數(shù)值技術(shù)主要關(guān)注數(shù)值格式和計算方法[5],均直接采用其中一種方案計算混合氣體輸運系數(shù),并沒有對比研究輸運系數(shù)各計算方案對流場數(shù)值模擬的影響。

因此,為有效運用混合氣體輸運系數(shù)近似計算方案,降低多組分流場數(shù)值模擬中對混合氣體輸運過程的計算量,進一步提高流場求解效率,有必要開展混合氣體輸運系數(shù)近似計算方案對流場計算的影響特點研究。本文基于Burrow和Kurkov實驗構(gòu)型及其實驗數(shù)據(jù)[6],進行了多組分混合氣體輸運系數(shù)計算方案影響數(shù)值模擬,對比研究了混合氣體分子導(dǎo)熱系數(shù)、組分?jǐn)U散系數(shù)的不同近似計算方案對邊界層及射流剪切層流場計算的影響,以期為在超燃沖壓發(fā)動機燃燒室流場的數(shù)值模擬中有效運用混合氣體輸運系數(shù)求解方案提供技術(shù)基礎(chǔ)。

1 計算模型及參數(shù)定義

1.1 計算模型及來流條件

Burrow和Kurkov實驗研究了二維后臺階燃燒室中的順噴氫燃料與來流的混合、燃燒特性,給出了較為詳細的燃燒室出口溫度、壓力及組分分布,該實驗已成為常用的驗證算例[7-9]。

圖1為實驗(計算)模型構(gòu)型示意圖。該實驗?zāi)Ｐ驮谂_階之前為等截面構(gòu)型,在臺階之后為擴張截面構(gòu)型,用下壁面的擴張抵消燃燒室下游邊界層的影響。實驗過程中,模型與燃燒加熱風(fēng)洞的噴管相連,氫氣燃料由實驗?zāi)Ｐ拖卤诿娴暮笈_階以聲速注入,噴注壓強為一個大氣壓,氫氣燃料出流方向與污染空氣來流平行、同向(圖中從左至右)。

圖1 實驗?zāi)Ｐ蜆?gòu)型示意[6](單位：cm)Fig.1 Geometry of test section[6](unit：cm)

表1 計算及實驗的來流參數(shù)Table 1 Inflow parameters for the CFD and test

本文針對純混合工況開展了數(shù)值模擬,其來流條件由表1給出(純混合實驗時,污染空氣來流中沒有O2組分,以確保不發(fā)生化學(xué)反應(yīng))。

1.2 計算區(qū)域及邊界條件

計算區(qū)域如圖1所示,為保證流向后臺階位置的邊界層厚度為0.01 m,自后臺階(圖中0 m處)向上游延伸0.42 m。網(wǎng)格結(jié)點總數(shù)為461×161,其中臺階上游部分的網(wǎng)格結(jié)點數(shù)為201×101,近壁網(wǎng)格的法向距離為1×10-5m,以保證k-ωSST類湍流模型所需的近壁網(wǎng)格雷諾數(shù)(y＋＜2)。本文僅開展相同網(wǎng)格條件下不同近似計算方案對比研究,沒有特別考慮網(wǎng)格無關(guān)性的影響。同時,由于計算網(wǎng)格同向相鄰網(wǎng)格間距比小于1.25,并且受近壁網(wǎng)格雷諾數(shù)的限制,因此網(wǎng)格無關(guān)性影響相對較小。

在計算中,① 壁面邊界:根據(jù)實驗條件,設(shè)定模型壁面為等溫壁(壁溫為298 K),并采用無滑移條件;② 來流邊界:來流為超聲速流動,固定入口截面的來流參數(shù);③ 噴口邊界:固定噴嘴出口截面的噴流參數(shù);④ 出流邊界:核心區(qū)的出流為超聲速流動,對原始變量進行一階外推。

1.3 摩阻系數(shù)、傳熱系數(shù)定義

壁面摩阻系數(shù)、傳熱系數(shù)的計算式為[10]:

式中,ρ∞、V∞、Cp∞、Tt∞分別為來流核心區(qū)流場的密度、速度、定壓比熱及總溫,τw、qw分別為壁面切應(yīng)力和熱流量,定義為:

式中,μ、μt、κ、κt分別為流場中混合氣體的分子黏性系數(shù)、渦黏性系數(shù)、分子導(dǎo)熱系數(shù)、湍流導(dǎo)熱系數(shù),u、T分別為流場的速度、溫度。

2 控制方程和數(shù)值方法

2.1 控制方程和計算方法簡介

計算坐標(biāo)系下,含k-ωSST湍流模型[11]、組分方程的三維密度平均完全Navier-Stokes方程的守恒形式為:

這里,Q＝(ρ,ρu,ρv,ρw,ρE,ρk,ρω,ρfi)T為 守 恒變量,其中的E為單位質(zhì)量總內(nèi)能;E、F、G、Ev、Fv以及Gv分別為控制方程的對流通量項和擴散通量項;S為湍流及化學(xué)源項;J為兩坐標(biāo)系之間的體積變換系數(shù)。

計算中,采用Liou[12]的AUSM＋-up格式離散控制方程中的無黏項,采用二階中心差分離散黏性項,采用Jameson[13]的LU-SSOR方法進行隱式求解。不對S進行隱式化處理時,LU-SSOR方法可表示為(去掉了方程中的上標(biāo)“~”):

式中,σA＝max[|λ(A)|],λ(A)為A＝?E/?Q的特征值,A±[A±σAI],其它從略。

為了適用于計算超聲速自由剪切流,需對初始k-ωSST湍流模型進行可壓縮性修正,本文采用了Sarkar類型修正[14-15]。計算源項S時,對湍流源項采用隱式處理,對化學(xué)源項則采用了預(yù)處理[16]。

2.2 輸運模型

在本文所關(guān)心的多組分流場中,氣體輸運過程主要引起以下三種現(xiàn)象,即:由速度梯度引起動量輸運的黏性現(xiàn)象、由溫度梯度和組分?jǐn)U散引起能量輸運的導(dǎo)熱現(xiàn)象及由濃度梯度引起組分輸運的擴散現(xiàn)象。為了獲得流場中的動量、組分以及能量輸運強度,需要給出組分的分子輸運系數(shù),并以此求出混合氣體的分子輸運系數(shù)。在多組分流場中,精細計算混合氣體的分子輸運系數(shù)耗時巨大,通常采用近似計算方法,其常用近似計算公式[17]如下文所述。

2.2.1 混合氣體輸運系數(shù)的近似計算

計算混合氣體分子黏性系數(shù)μ的常用公式為Wilke的半經(jīng)驗關(guān)系式,即:

其中,Xi、Wi和μi分別為組分的摩爾分?jǐn)?shù)、摩爾質(zhì)量(分子量)及分子黏性系數(shù)。

計算混合氣體的分子導(dǎo)熱系數(shù)κ的常用近似公式類同Wilke的半經(jīng)驗公式,而較為精確的近似公式則為Brokaw的半經(jīng)驗公式,即:

其中,κ′i為組分i的導(dǎo)熱系數(shù)中的平動部分為由轉(zhuǎn)動、振動等自由度產(chǎn)生的附加項。

與流體力學(xué)控制方程中組分的平均速度、擴散速度和擴散通量相容,組分i穿過其余混合氣體的分子擴散系數(shù)Dim的半經(jīng)驗公式為[18]:

式中,Yi為組分i的質(zhì)量分?jǐn)?shù),Dij為各組分之間的雙組元擴散系數(shù)。實際計算中,普遍采用下述兩個半經(jīng)驗公式計算混合氣體分子擴散系數(shù),即:

公式(18)導(dǎo)自基于摩爾分?jǐn)?shù)梯度的質(zhì)量平均速度和擴散通量,公式(19)則導(dǎo)自基于摩爾分?jǐn)?shù)梯度的摩爾平均速度和擴散通量,雖然理論上均與流體力學(xué)控制方程存在相容性問題,但因其計算量相對較小,得到了普遍使用,公式(17)雖與流體力學(xué)控制方程相容,但是目前鮮見實際應(yīng)用。

2.2.2 導(dǎo)熱系數(shù)、擴散系數(shù)的關(guān)聯(lián)計算

多組分流場中,即使采用了半經(jīng)驗公式,求解混合氣體分子輸運系數(shù)仍是一個耗時運算,還需要在計算效率和計算精度間進行權(quán)衡。關(guān)聯(lián)計算方案預(yù)設(shè)常值普朗特數(shù)和施米特數(shù),根據(jù)混合氣體分子黏性系數(shù)計算導(dǎo)熱系數(shù)和擴散系數(shù),即:

該方案計算精度稍遜,但計算效率很高,為提高復(fù)雜流場求解效率,仍需研究其適用空間。

3 結(jié)果分析

3.1 基于Burrow&Kurkov實驗的驗證計算

驗證計算的來流條件仍由表1給出,混合氣體的輸運系數(shù)由半經(jīng)驗關(guān)系式求解。其中,混合氣體分子擴散系數(shù)由公式(17)計算。圖2為數(shù)值模擬給出的全流場壓強等值線,并且給出了基于壓強的局部流場波系結(jié)構(gòu)。等直段的流場簡單,屬于受限平板流。由于氫氣燃料的平行、同向噴注,且噴注靜壓強為一個大氣壓,后臺階區(qū)域的流場結(jié)構(gòu)主要表現(xiàn)為噴口唇口構(gòu)型所引起的底部流場結(jié)構(gòu),其中剪切層上緣的左傳波系與上壁面作用,其反射波系在與下壁面作用前傳出計算域,上緣波系與噴口的噴流波系疊加,經(jīng)過下壁面反射、底部及射流自由剪切層的透射、反射,形成了噴口近場波系。

圖2 壓強等值線Fig.2 Pressure contour

圖3 溫度等值線Fig.3 Temperature contour

圖3為流場等溫度線。由圖可見,近壁區(qū)域的溫度梯度大,掩蓋了核心區(qū)溫度場變化細節(jié),沒有刻畫出類似壓強場給出的波系結(jié)構(gòu)。射流溫度接近壁溫,射流上緣的溫度梯度大,在一定程度上顯示了臺階下游區(qū)域射流剪切層的發(fā)展情況。

圖4給出了臺階及出口位置的皮托壓強、總溫和組分剖面。圖中Experiment指實驗,CFD SST指采用的湍流模型為k-ωSST模型,CFD SST-CC指采用的湍流模型僅對k-ωSST模型的k-ε模型部分進行可壓縮性修正。湍流可壓縮性效應(yīng)主要作用于邊界層外緣、自由剪切層,對出口流場參數(shù)剖面影響較為明顯。整體上,計算的皮托壓強、總溫及組分剖面雖然與實驗值有偏差,但是反映了變化規(guī)律,表明該數(shù)值技術(shù)可以用于本文的數(shù)值模擬研究。

圖4 皮托壓強、總溫及組分剖面Fig.4 Profile of P p,T t and f i

3.2 輸運系數(shù)近似計算方案對比分析

根據(jù)表1的來流條件,針對輸運系數(shù)近似計算方案,開展了湍流/層流狀態(tài)下的數(shù)值模擬,具體的算例由表2給出。計算中,混合氣體分子黏性系數(shù)均由公式(12)和公式(13)計算給出,并且作為近似處理其它輸運系數(shù)的基礎(chǔ)。

表2 輸運系數(shù)近似計算方案Table 2 Approximate formula for K and D im

根據(jù)圖2、圖3所示,臺階上游流場屬于受限平板邊界層流場,顯示了邊界層的流動,下游流場包含射流剪切層流動,顯示了臺階下游射流剪切層發(fā)展情況。由此,本文通過對一個狀態(tài)的計算可以同時對比輸運系數(shù)的近似計算方案對邊界層、射流剪切層流動的影響特點。

3.2.1 近似計算方案對壁面邊界層流動的影響

圖5為根據(jù)導(dǎo)熱系數(shù)、黏性系數(shù)以及定壓比熱計算給出的混合氣體普朗特數(shù)等值線,顯示了混合氣體普朗特數(shù)在流場中的變化特點。計算中,臺階上游沒有燃料注入,污染空氣來流的實際普朗特數(shù)變化較為平緩,且與本文預(yù)先給定的普朗特數(shù)相差不大(本文給定Pr＝0.7)。由于燃料的噴入,臺階下游近壁區(qū)域的氫氣組分占優(yōu),使得該區(qū)域的混合氣體普朗特數(shù)與預(yù)設(shè)普朗特數(shù)相差較大,同時射流上緣區(qū)域的普朗特數(shù)場梯度大,也可在一定程度上顯示臺階下游區(qū)域射流剪切層的發(fā)展情況。

圖5 普朗特數(shù)等值線Fig.5 Prandtl number contour

圖6為下壁面在臺階上游10 mm位置的邊界層速度和溫度剖面。這里,“Lam”表示層流條件下的計算值。因計算混合氣體分子黏性系數(shù)的公式相同,速度剖面幾近相同,算例1、算例2和算例3采用同一公式計算分子導(dǎo)熱系數(shù),因此溫度剖面相同,算例4采用關(guān)聯(lián)法獲得分子導(dǎo)熱系數(shù),其溫度剖面與其余算例的結(jié)果出現(xiàn)差異。

圖6 x＝-10 mm位置的速度、靜溫剖面Fig.6 Profile of V x and T s at x＝-10 mm

圖7為模型下壁面的壁面摩阻系數(shù)、傳熱系數(shù)分布。算例1、算例2和算例3同樣有相同的壁面摩阻系數(shù)和壁面?zhèn)鳠嵯禂?shù),算例4的摩阻系數(shù)與其余算例值接近,傳熱系數(shù)則與其余算例的結(jié)果有一定的差異。結(jié)合圖6所給出的速度、溫度剖面可見,分子導(dǎo)熱系數(shù)的計算方案對邊界層溫度剖面及壁面?zhèn)鳠嵯禂?shù)影響較大,通過能量方程間接影響了邊界層速度剖面及壁面摩阻系數(shù)分布,但是其影響程度相對較小。對比數(shù)據(jù)的差異表明,雖然混合氣體分子導(dǎo)熱系數(shù)計算方案對邊界層速度及溫度剖面、壁面摩阻系數(shù)及傳熱系數(shù)分布計算產(chǎn)生了影響,但是其偏差仍在可接受的范圍。由此,通過預(yù)先設(shè)定適當(dāng)普朗特數(shù),由分子黏性系數(shù)關(guān)聯(lián)混合氣體分子導(dǎo)熱系數(shù),在降低計算量的基礎(chǔ)上可獲得與采用公式計算混合氣體分子導(dǎo)熱系數(shù)的直接計算方案相近的壁面摩阻系數(shù)和壁面?zhèn)鳠嵯禂?shù)。

3.2.2 近似計算方案對射流剪切層流動的影響

由圖3、圖5的參數(shù)等值線可見,臺階下游的射流剪切層在發(fā)展過程中逐步與壁面邊界層產(chǎn)生了流動干擾。為排除其干擾因素的影響,本文在臺階下游選擇10、100 mm位置對比輸運系數(shù)近似計算方案對射流剪切層流動的影響,并同時給出燃燒室出口位置(356 mm)的參數(shù)剖面,以對比顯示輸運系數(shù)近似計算方案對綜合流場的影響。

圖7 摩阻系數(shù)、傳熱系數(shù)分布Fig.7 Distribution of C f and C h

圖8為臺階下游10、100、356 mm位置的流場氫氣組分剖面,顯示了近似計算方案對射流剪切層區(qū)域內(nèi)氫氣組分剖面的影響特點。由圖可見,層流狀態(tài)下,算例1、算例2和算例3的流場組分剖面近乎相同,算例4的組分剖面則與其余算例的結(jié)果出現(xiàn)差異;湍流狀態(tài)下,混合氣體“渦黏性”輸運效應(yīng)掩蓋了對應(yīng)的分子輸運效應(yīng),使各算例(輸運系數(shù)近似計算方案)給出了相同的組分剖面。

圖9為臺階下游10、100、356 mm位置的流場速度、溫度剖面。湍流狀態(tài)下,各算例給出相同的流場速度、溫度剖面;層流狀態(tài)下,算例1、算例2和算例3的速度、溫度剖面非常接近,而算例4的速度、溫度剖面則與算例1、算例2以及算例3的對應(yīng)結(jié)果有相對較大的差異。結(jié)合圖8可見,隨著流場向下游發(fā)展,層流流場中組分剖面的差異主要出現(xiàn)在組分?jǐn)U散波陣面的上升沿,且差異程度逐漸減弱。與此相呼應(yīng),層流流場中對應(yīng)的速度、溫度剖面差異程度也隨流場向下游的發(fā)展而減弱。

圖8 不同x向位置的H2組分剖面Fig.8 Profile of H2 composition at different x position

圖9 不同x向位置的速度、靜溫剖面Fig.9 Profile of v x and T s at different x position

圖10 不同x向位置的H2組分剖面Fig.10 Profile of H2 composition at different x position

在算例1的基礎(chǔ)上,開展了兩種可壓縮性修正方案下的流場數(shù)值模擬,進一步對比顯示混合氣體的“渦黏性”輸運對組分剖面的影響。圖10給出了臺階下游10、100、356 mm位置的流場氫氣組分的摩爾分?jǐn)?shù)剖面。圖中CFD SST-CC-1指采用的湍流模型對k-ωSST模型中的k-ε和k-ω模型部分進行可壓縮性修正。結(jié)合圖8、圖9可見,湍流流場中各湍流模型的“渦黏性”輸運系數(shù)是導(dǎo)致流場組分剖面(包括速度剖面、溫度剖面)出現(xiàn)差異的主要因素,此時混合氣體分子輸運系數(shù)近似計算方案對流場數(shù)值模擬的影響則退居于次要地位。

4 結(jié) 論

利用Burrow和Kurkov的實驗構(gòu)型及相關(guān)研究數(shù)據(jù),本文開展了多組分混合氣體導(dǎo)熱系數(shù)、組分?jǐn)U散系數(shù)近似計算方案對流場參數(shù)影響的數(shù)值模擬研究,獲得以下主要結(jié)論:

(1)與采用由相關(guān)公式近似計算混合氣體導(dǎo)熱系數(shù)、擴散系數(shù)的直接近似計算方案相比,采用由普朗特數(shù)、施密特數(shù)關(guān)聯(lián)混合氣體導(dǎo)熱系數(shù)、擴散系數(shù)的關(guān)聯(lián)近似計算方案,流場的數(shù)值模擬在獲得合理計算結(jié)果的同時,可大為降低計算量。

(2)拋開相容性問題,混合氣體擴散系數(shù)直接近似計算方案對流場參數(shù)計算影響不明顯。

(3)相對層流計算,湍流計算時邊界層/剪切層區(qū)域的混合氣體湍流輸運系數(shù)大于相應(yīng)的分子輸運系數(shù),由此弱化了混合氣體分子輸運系數(shù)不同近似計算方案對流場數(shù)值模擬差異的影響。

(4)相對輸運系數(shù)近似計算方案,湍流可壓縮修正模型對組分剖面的影響更為明顯。

本文主要對比研究了冷態(tài)流場下的混合氣體輸運系數(shù)計算方案對流場參數(shù)計算的影響,獲得了有一定參考價值的結(jié)論。同時,基于分子運動論的氣體輸運系數(shù)理論計算公式表明,施密特數(shù)受流場壓強、溫度和組成等因素的影響較大,因此下一步需要開展化學(xué)非平衡狀態(tài)下的計算研究,以期更加全面地數(shù)值模擬對比、分析研究輸運系數(shù)計算方案對流場參數(shù)計算的影響特點。