本發(fā)明屬于計(jì)算流體力學(xué)(cfd)領(lǐng)域，具體涉及一種考慮邊界層燃燒放熱效應(yīng)的可壓縮壁函數(shù)計(jì)算方法，可有效提高超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道摩阻/熱流計(jì)算效率。

背景技術(shù)：

超燃沖壓發(fā)動(dòng)機(jī)被認(rèn)為是未來吸氣式高超聲速飛行器的理想動(dòng)力裝置，其技術(shù)發(fā)展具有重大戰(zhàn)略意義，因此受到了各國的高度重視。由于超聲速燃燒相關(guān)實(shí)驗(yàn)實(shí)施困難且某些參數(shù)限于測量技術(shù)而無法直接獲取，計(jì)算流體力學(xué)(cfd)在超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道的流場分析中得到了廣泛的應(yīng)用，并逐漸成為超燃沖壓發(fā)動(dòng)機(jī)設(shè)計(jì)中非常重要的技術(shù)手段。

超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道壁面摩阻和熱流的數(shù)值預(yù)測對計(jì)算流體力學(xué)(cfd)技術(shù)構(gòu)成了很大挑戰(zhàn)。超聲速可壓縮湍流邊界層內(nèi)的速度和溫度在近壁區(qū)呈現(xiàn)劇烈變化，使得摩阻和熱流的計(jì)算結(jié)果高度依賴于壁面法向網(wǎng)格的劃分方式，特別是壁面法向第一層網(wǎng)格的間距。很多研究均指出，要準(zhǔn)確計(jì)算高速可壓縮湍流邊界層的摩阻和熱流，壁面法向第一層網(wǎng)格間距一般要滿足y⁺＝ρwuτδy/μw<1的要求(其中ρw為壁面密度，uτ為壁面剪切速度，δy為壁面第一層網(wǎng)格間距，μw為壁面分子粘性系數(shù))。這要求近壁面法向網(wǎng)格需非常之密，一般要求在10^-6m甚至更小。對于實(shí)際的復(fù)雜超燃沖壓發(fā)動(dòng)機(jī)構(gòu)型,由于網(wǎng)格過密而大大限制了積分時(shí)間步長，再加上內(nèi)部復(fù)雜化學(xué)反應(yīng)模擬的高計(jì)算量，導(dǎo)致計(jì)算周期過長甚至是無法接受的。因此,對整個(gè)湍流邊界層進(jìn)行詳細(xì)模擬而得到壁面摩阻和熱流的計(jì)算方式在實(shí)際工程應(yīng)用時(shí)存在困難，因此很多工程應(yīng)用中均采用了壁函數(shù)邊界條件的方法來得到壁面摩阻和熱流，即對近壁區(qū)的湍流邊界層不進(jìn)行模擬而直接應(yīng)用速度和溫度壁面律來模擬速度和溫度的分布，從而間接得到壁面摩擦和熱流。因此，壁函數(shù)方法的使用可有效放寬壁面摩阻和熱流模擬對網(wǎng)格密度的限制，從而降低計(jì)算量。超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流動(dòng)的復(fù)雜性給壁函數(shù)方法的應(yīng)用帶來很大困難：燃燒化學(xué)反應(yīng)會(huì)在可壓縮湍流邊界層內(nèi)發(fā)生，國內(nèi)外大量的實(shí)驗(yàn)和數(shù)值模擬研究均顯示湍流邊界層內(nèi)的燃燒放熱效應(yīng)將大幅改變壁面摩阻及熱流，因此壁函數(shù)模型中必須要考慮邊界層內(nèi)的燃燒放熱效應(yīng)。然而，國內(nèi)外現(xiàn)有的所有壁函數(shù)方法均是建立在邊界層內(nèi)無化學(xué)反應(yīng)的假設(shè)基礎(chǔ)上，目前尚無考慮邊界層內(nèi)化學(xué)反應(yīng)放熱效應(yīng)的壁函數(shù)模型，這也導(dǎo)致壁函數(shù)方法不能應(yīng)用于超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道的摩阻/熱流模擬中。

技術(shù)實(shí)現(xiàn)要素：

為了能在超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道壁面摩阻和熱流的數(shù)值模擬中應(yīng)用壁函數(shù)方法，從而大大提高計(jì)算效率，本發(fā)明提出了一種考慮邊界層燃燒放熱效應(yīng)的新型可壓縮壁函數(shù)計(jì)算方法。通過理論推導(dǎo)提出了燃燒化學(xué)反應(yīng)條件下的速度壁面律和溫度壁面律，基于此構(gòu)造了可同時(shí)兼容有/無化學(xué)反應(yīng)可壓縮湍流邊界層的統(tǒng)一壁函數(shù)計(jì)算方法，并設(shè)計(jì)了與cfd程序耦合求解的詳細(xì)流程(見參考文獻(xiàn)[1]：gaoz.x.,jiangc.w.,pans.w.,leec.h.,combustionheatreleaseeffectonsupersoniccompressibleturbulentboundarylayers,aiaajournal53(7)(2015)1949-1968)。

本發(fā)明提出的一種考慮邊界層燃燒放熱效應(yīng)的新型可壓縮壁函數(shù)計(jì)算方法，可大幅提高超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道摩阻/熱流的數(shù)值模擬效率。壁函數(shù)方法中最核心的部分即為速度壁面律和溫度壁面律，因此要建立考慮邊界層內(nèi)燃燒放熱效應(yīng)的壁函數(shù)方法，則需要考慮燃燒的放熱效應(yīng)對可壓縮速度壁面律和溫度壁面律的影響?，F(xiàn)有壁函數(shù)邊界條件中的速度壁面律和溫度壁面律建立均是針對充分發(fā)展的二維平板湍流邊界層建立的，而要建立燃燒放熱效應(yīng)下的速度壁面律和溫度壁面律需要構(gòu)造包含邊界層燃燒的標(biāo)準(zhǔn)流動(dòng)模型。

具體的，本發(fā)明的考慮邊界層燃燒放熱效應(yīng)的可壓縮壁函數(shù)計(jì)算方法包括如下步驟：

步驟1：燃燒放熱效應(yīng)下的速度壁面律；

在推導(dǎo)二維充分發(fā)展湍流邊界層的速度壁面律過程中，需要邊界層內(nèi)沿法向速度與密度的關(guān)系式，邊界層內(nèi)的燃燒放熱會(huì)改變這一關(guān)系式。為了得到燃燒條件下的邊界層內(nèi)速度和密度關(guān)系式，引入燃燒理論中的shvab-zeldovich耦合參數(shù)zho，zot和zht：

其中，和分別為邊界層內(nèi)氧氣和氫氣組分的質(zhì)量分?jǐn)?shù)濃度，h為氣體的敏感焓，而常數(shù)f和δq分別為反應(yīng)當(dāng)量比以及單位質(zhì)量氫氣燃燒放出的熱量。令z＝zho，zot或zht，則z的輸運(yùn)方程與速度u的輸運(yùn)方程相似，且在簡化模型下兩者的邊界條件也相似，因此有：

z＝d1+d2u(2)

其中，d1和d2為常數(shù)，對于每個(gè)耦合參數(shù)由相應(yīng)的邊界條件確定。

所述的簡化模型是指對于包含邊界層燃燒的二維平板湍流邊界層流動(dòng)模型，在快速化學(xué)反應(yīng)假設(shè)下，火焰前鋒上部氫氣濃度為0，而下部氧氣濃度為0，由此得到簡化模型。

接下來，利用zho，zot和zht三個(gè)耦合參數(shù)與速度u的關(guān)系可以推導(dǎo)邊界層內(nèi)溫度與速度的關(guān)系，在邊界層內(nèi)法向壓強(qiáng)梯度為零的假設(shè)下溫度與速度的關(guān)系可轉(zhuǎn)化為邊界層內(nèi)燃燒混合氣體的密度和速度的關(guān)系。結(jié)果由下式給出：

其中，

上式中ρ和u分別代表燃燒混合氣體的密度和速度，ρw、hw和chw則表示壁面上混合氣體的密度、敏感焓和氫氣質(zhì)量分?jǐn)?shù)，而ue、he、coe則代表邊界層外緣混合氣體的速度、敏感焓和氧氣質(zhì)量分?jǐn)?shù)。另外，r為常數(shù)恢復(fù)因子。α、b1、b2和a均為中間變量。

以上公式中下標(biāo)中含有e代表為邊界層外緣參數(shù)，下標(biāo)含有w代表壁面參數(shù)。

在得到以上密度與速度的關(guān)系式(3)之后，根據(jù)混合長度以及等切應(yīng)力假設(shè)有：

其中，κ為卡門常數(shù)，取值0.41。u⁺和y⁺分別為以壁面剪切速度uτ等壁面參數(shù)無量綱化的速度和壁面距離。接下來將式(3)中的ρw/ρ表達(dá)式代入式(4)，進(jìn)行整理后得到：

積分上式(5)可得到：

其中，

q＝(β²+4γ)^1/2,

qw、tw、μw、kw分別為壁面上的熱流、溫度、分子粘性系數(shù)、熱傳導(dǎo)系數(shù)。此外，b為速度壁面律中的常數(shù)，取值5.0。為常數(shù)。式(7)即為包含燃燒放熱效應(yīng)的速度對數(shù)律，此式中并未顯示包含燃燒放熱相關(guān)的變量，燃燒放熱效應(yīng)是通過改變式(7)中參數(shù)γ和β中的壁面剪切速度uτ和熱流qw的值而改變速度型的。另外，γ和β中的參數(shù)hw、分子粘性系數(shù)μw和熱傳導(dǎo)系數(shù)kw在包含邊界層燃燒的條件下也不再只與壁面溫度相關(guān)，還會(huì)受壁面組分濃度的影響，這一效應(yīng)也很重要。

式(7)只在湍流邊界層的對數(shù)層成立，在實(shí)際應(yīng)用時(shí)需要構(gòu)造在湍流邊界層內(nèi)層一致成立的速度壁面律。本發(fā)明中引入了spalding構(gòu)造湍流邊界層內(nèi)層一致有效不可壓速度壁面律的方法，最后得到的考慮邊界層燃燒放熱效應(yīng)的速度壁面律由下式給出：

該速度壁面律的一個(gè)重要特征是在形式上對于邊界層內(nèi)有/無燃燒放熱效應(yīng)的情況均成立，這對于實(shí)際壁函數(shù)方法在超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流場的應(yīng)用非常重要，可以避免在應(yīng)用時(shí)要判斷流場當(dāng)?shù)剡吔鐚觾?nèi)是否有燃燒化學(xué)反應(yīng)發(fā)生而實(shí)現(xiàn)兩種速度壁面律的切換。

步驟2：燃燒放熱效應(yīng)下的溫度壁面律；

對于壁函數(shù)中需要的溫度壁面律，傳統(tǒng)方法中應(yīng)用crocco-busemann速度-溫度關(guān)系式來獲得，但在燃燒放熱效應(yīng)下crocco-busemann速度-溫度關(guān)系式必然不再成立。本發(fā)明中提出采用另外的策略，即不采用速度-溫度的關(guān)系，而采用速度與靜焓h的關(guān)系。此時(shí)的靜焓為所有組分敏感焓和生成焓的和，因此靜焓并不會(huì)因?yàn)槿紵瘜W(xué)反應(yīng)的發(fā)生而變化。由二維條件下總焓(靜焓加動(dòng)能)與速度輸運(yùn)方程的相似性以及相應(yīng)的邊界條件，可得到速度u-靜焓h的關(guān)系式如下：

其中，hw和he分別為壁面上和邊界層外緣混合氣體的靜焓，r為常數(shù)恢復(fù)因子。通過數(shù)學(xué)變換可去掉式(10)中邊界層外緣的參數(shù)并將速度u的形式變換為無量綱的速度u⁺的形式，最后得到速度-靜焓的關(guān)系式如下：

其中，cpw為壁面上混合氣體的定壓比熱容。

此速度-靜焓關(guān)系式在無邊界層燃燒的情況下會(huì)自動(dòng)退化為crocco-busemann速度-溫度關(guān)系式，因此其也可兼容邊界層內(nèi)有或無燃燒放熱效應(yīng)的不同情況。

步驟3：考慮邊界層燃燒放熱效應(yīng)的壁函數(shù)方法；

基于以上所發(fā)展的速度壁面律式(9)和速度-靜焓關(guān)系式(11)可以構(gòu)造考慮邊界層燃燒放熱效應(yīng)的壁函數(shù)方法。該壁函數(shù)需要與cfd計(jì)算程序耦合求解，具體的求解流程如下：

首先，根據(jù)速度壁面律式(9)以及cfd程序計(jì)算得到的壁面第一層網(wǎng)格上的速度u1，迭代求解壁面剪切速度uτ，進(jìn)而求得壁面剪切應(yīng)力τw＝ρwuτ²，即作為當(dāng)?shù)乇诿娴哪Σ翍?yīng)力。特別地，式(9)中γ和β的表達(dá)式中μw和kw一定要基于wilke的混合律獲得。此外，還需要將求得的當(dāng)?shù)乇诿娴哪Σ翍?yīng)力τw值加入到cfd粘性子程序中替換掉動(dòng)量方程中對應(yīng)的壁面切應(yīng)力值τ1/2。

其次,根據(jù)以上速度壁面律中得到的u1⁺以及cfd程序計(jì)算得到的壁面第一層網(wǎng)格上的靜焓h1,由速度-靜焓關(guān)系式(11)式得到熱流qw,所求得的熱流qw即作為當(dāng)?shù)乇诿娴臒崃鳌Ｍ瑯?需要將求得的熱流qw值加入到cfd粘性子程序中替換掉能量方程中對應(yīng)的壁面熱傳導(dǎo)項(xiàng)q1/2。

本發(fā)明的優(yōu)點(diǎn)在于：

(1)利用本發(fā)明中的壁函數(shù)方法可大幅提高超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道摩阻/熱流的模擬效率。如果不采用壁函數(shù)方法，計(jì)算時(shí)壁面第一層網(wǎng)格的無量綱化的壁面距離y⁺需要在1左右。而采用壁函數(shù)之后，實(shí)際計(jì)算結(jié)果表明壁面第一層網(wǎng)格的無量綱化的壁面距離y⁺可在100左右仍可得到可靠的摩阻和熱流，因此可有效提高數(shù)值模擬的積分時(shí)間步長，大幅縮短計(jì)算周期。

(2)技術(shù)實(shí)現(xiàn)難度低。本發(fā)明中的壁函數(shù)方法在使用時(shí)可以作為原有cfd程序的一種壁面邊界條件，只需作為子程序加入，而與cfd程序耦合計(jì)算時(shí)對原有程序的修改很少。

(3)所提出的壁函數(shù)方法同時(shí)適用于有或無邊界層燃燒放熱效應(yīng)的情況。本壁函數(shù)中的速度壁面律和溫度壁面律的形式在邊界層內(nèi)有/無燃燒化學(xué)反應(yīng)的情況下均成立，因此此壁函數(shù)方法對單組分無化學(xué)反應(yīng)、多組分混合、多組分有化學(xué)反應(yīng)的可壓縮湍流邊界層均適用，無需在計(jì)算時(shí)判斷流場邊界層內(nèi)有無化學(xué)反應(yīng)。

附圖說明

圖1為包含邊界層燃燒的二維流動(dòng)模型；

圖2為邊界層燃燒的二維簡化流動(dòng)模型；

圖3為壁函數(shù)子程序與cfd主程序耦合求解方法：(a)粘性項(xiàng)；(b)熱傳導(dǎo)項(xiàng)；

圖4為包含邊界層燃燒條件下摩阻/熱流計(jì)算的實(shí)驗(yàn)算例；

圖5為數(shù)值模擬得到的溫度分布云圖；

圖6為無壁函數(shù)密網(wǎng)格(y⁺約為1)條件下摩阻/熱流計(jì)算結(jié)果的驗(yàn)證；

圖7為網(wǎng)格稀疏后壁面第一層網(wǎng)格的y⁺沿流向的分布；

圖8為采用壁函數(shù)在稀網(wǎng)格條件下計(jì)算的摩阻/熱流結(jié)果對比(其中cwfbc代表本專利中的壁函數(shù)計(jì)算方法，nicholswfbc代表現(xiàn)有未考慮燃燒放熱效應(yīng)的壁函數(shù)方法)。

具體實(shí)施方式

下面結(jié)合附圖和實(shí)施例對本發(fā)明進(jìn)行詳細(xì)說明。

本發(fā)明提供一種考慮邊界層燃燒放熱效應(yīng)的可壓縮壁函數(shù)計(jì)算方法，基于圖1所示的包含邊界層燃燒的二維平板湍流邊界層流動(dòng)模型進(jìn)行理論推導(dǎo)，以獲得燃燒放熱效應(yīng)下的速度壁面律和溫度壁面律。在圖1中的二維流動(dòng)模型，氫氣由狹槽平行于壁面噴射進(jìn)入壁面附近空氣流動(dòng)形成的可壓縮湍流邊界層內(nèi)，此后氫氣/空氣擴(kuò)散火焰在邊界層內(nèi)發(fā)生并放出熱量，在快速化學(xué)反應(yīng)假設(shè)下，火焰前鋒上部氫氣濃度ch為0，而火焰下部氧氣濃度co為0，壁面上的溫度tw為常數(shù)，表示壁面上混合氣體中氫氣質(zhì)量分?jǐn)?shù)chw為常數(shù)，由此得到如圖2所示的二維簡化流動(dòng)模型。

本發(fā)明提出了一種考慮邊界層燃燒放熱效應(yīng)的新型可壓縮壁函數(shù)計(jì)算方法，可大幅提高超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道摩阻/熱流的數(shù)值模擬效率。該壁函數(shù)方法可作為一種壁面邊界條件與cfd程序耦合求解，可在較稀疏的壁面法向網(wǎng)格條件下得到可靠的壁面摩阻和熱流。該壁函數(shù)計(jì)算方法具體的實(shí)施步驟如下：

步驟1：編寫壁函數(shù)計(jì)算方法子程序，作為一種壁面邊界條件加入現(xiàn)有cfd計(jì)算程序中。該子程序要實(shí)現(xiàn)基于以上速度壁面律式(9)和溫度-靜焓關(guān)系式(11)迭代求解其中的未知量uτ和qw。該子程序需要由cfd主程序輸入壁面上的參數(shù)值包括ρw、tw、hw、μw、kw、cpw，以及壁面法向第一層網(wǎng)格上計(jì)算得到的速度值u1以及靜焓值h1；

步驟2：按以下方法修改原有cfd主程序：由壁函數(shù)子程序得到的uτ可得到壁面摩擦應(yīng)力τw，將此壁面摩擦應(yīng)力τw值加入到cfd粘性子程序中替換掉動(dòng)量方程中對應(yīng)的壁面切應(yīng)力值,如圖3(a)中所示的τ1/2。此外，將壁函數(shù)子程序求得的qw值加入到cfd粘性子程序中替換掉能量方程中對應(yīng)的壁面熱傳導(dǎo)項(xiàng),如圖3(b)中所示的q1/2。

步驟3：在實(shí)際的超燃沖壓發(fā)動(dòng)機(jī)內(nèi)流道數(shù)值模擬時(shí)，將壁面第一層網(wǎng)格的寬度保持y⁺在100左右。計(jì)算時(shí)實(shí)現(xiàn)cfd主程序和壁函數(shù)子程序的迭代求解，收斂后壁面的摩擦應(yīng)力和熱流值以壁函數(shù)子程序求解得到的τw和qw為準(zhǔn)。

實(shí)施例：基于本發(fā)明中提出的考慮邊界層燃燒放熱效應(yīng)的可壓縮壁函數(shù)方法，對圖4所示的一個(gè)包含邊界層燃燒的實(shí)驗(yàn)算例進(jìn)行了模擬。空氣來流條件為：馬赫數(shù)4.42，溫度1120k，壓強(qiáng)83kpa。氫氣在上壁面通過一個(gè)高度為3mm的臺(tái)階以音速噴射入空氣形成的超聲速可壓縮湍流邊界層內(nèi)，圖5中給出的計(jì)算得到的溫度分布云圖可以看到邊界層內(nèi)形成的擴(kuò)散燃燒火焰。圖6首先展示了在密網(wǎng)格(壁面第一層網(wǎng)格尺度5×10^-7m，對應(yīng)的y⁺在1左右)條件下不采用壁函數(shù)計(jì)算得到的有/無邊界層燃燒的摩阻和熱流之比與實(shí)驗(yàn)數(shù)據(jù)的對比，可見與實(shí)驗(yàn)數(shù)據(jù)相當(dāng)吻合。接下來，生成兩套稀網(wǎng)格，壁面第一層網(wǎng)格尺度分別為5×10^-5m和1×10^-4m，其所對應(yīng)的y⁺在100左右，如圖7所示。以經(jīng)過實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證的5×10^-7m密網(wǎng)格摩阻和熱流結(jié)果為準(zhǔn)，圖8將基于兩套稀網(wǎng)格并采用本發(fā)明提出的壁函數(shù)計(jì)算方法所得到的摩阻和熱流與密網(wǎng)格結(jié)果進(jìn)行了對比，可以看到本發(fā)明的壁函數(shù)計(jì)算方法可保證在兩套稀網(wǎng)格上均可獲得與密網(wǎng)格一致的摩阻和熱流分布。顯然，稀網(wǎng)格可降低計(jì)算周期，本算例中1×10^-4m稀網(wǎng)格的計(jì)算時(shí)間只有5×10^-7m密網(wǎng)格計(jì)算時(shí)間的20％，因此本發(fā)明中的壁函數(shù)方法大幅提高了計(jì)算效率。此外，圖8中的雙點(diǎn)劃線同時(shí)還給出了本算例使用現(xiàn)有未考慮邊界層燃燒放熱效應(yīng)的壁函數(shù)在稀網(wǎng)格上的計(jì)算結(jié)果，可以看到熱流結(jié)果嚴(yán)重偏離了準(zhǔn)確值，表明本發(fā)明中壁函數(shù)方法可有效考慮邊界層內(nèi)燃燒對摩阻/熱流的影響。