局部時(shí)間步法在低馬赫燃燒模擬中的適用性研究

姚?衛(wèi)，孫?超，劉?杭，吳?梅，楊少波，曹順利

局部時(shí)間步法在低馬赫燃燒模擬中的適用性研究

姚?衛(wèi)1, 2，孫?超3，劉?杭1，吳?梅4，楊少波3，曹順利3

(1. 中國(guó)科學(xué)院力學(xué)研究所高溫氣體動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100190；2. 中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院，北京 100049；3.中國(guó)船舶工業(yè)系統(tǒng)工程研究院，北京 100094；4.合肥工業(yè)大學(xué)土木與水利工程學(xué)院，合肥 230009)

在復(fù)雜燃燒模擬中由于整場(chǎng)流速的不均勻性和局部網(wǎng)格尺寸的差異，各局部流場(chǎng)區(qū)域的CFL數(shù)差異較大．傳統(tǒng)的基于整場(chǎng)最大CFL數(shù)定義的整體時(shí)間步法嚴(yán)重制約計(jì)算效率．本文首次考察了基于當(dāng)?shù)谻FL數(shù)限制的局部時(shí)間步法在低馬赫數(shù)湍流燃燒模擬中的適用性．對(duì)開(kāi)放空間中甲烷池火(1065萬(wàn)網(wǎng)格)和封閉空間建筑火災(zāi)(320萬(wàn)網(wǎng)格)的大渦模擬表明，采用局部時(shí)間步法相比于整體時(shí)間步法分別實(shí)現(xiàn)了6倍和8倍的加速比．加速比隨網(wǎng)格尺度減小呈增加趨勢(shì)．研究進(jìn)一步從兩個(gè)方面驗(yàn)證了局部時(shí)間步法在低馬赫數(shù)燃燒模擬中的準(zhǔn)確性：①與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比表明，由于低馬赫數(shù)燃燒的準(zhǔn)穩(wěn)態(tài)特性，局部和整體時(shí)間步法均較為準(zhǔn)確地預(yù)測(cè)了溫度的時(shí)間變化特性；②對(duì)時(shí)均流場(chǎng)的比較表明，除微量痕跡物(質(zhì)量分?jǐn)?shù)小于0.1%)以外，兩種方法對(duì)時(shí)均溫度、時(shí)均速度和氧氣體積分?jǐn)?shù)的預(yù)測(cè)差異均較?。芯恐羞€對(duì)現(xiàn)有的PaSR湍流燃燒模型和壓力求解算法進(jìn)行了改進(jìn)和優(yōu)化，以分別提高其物理準(zhǔn)確性和魯棒性．

湍流燃燒；局部時(shí)間步法；大渦模擬；CFL數(shù)；OpenFOAM

低馬赫數(shù)燃燒是日常生活中最常見(jiàn)的一種燃燒形式，例如各類火災(zāi)和非動(dòng)力裝置類鍋爐內(nèi)部燃燒．在各類跨聲速甚至高超聲速航空航天發(fā)動(dòng)機(jī)燃燒室中也存在較大區(qū)域的局部低馬赫數(shù)燃燒．低馬赫數(shù)燃燒的特點(diǎn)是：①壓力波動(dòng)較小近似恒定，因此可按照不可壓縮流體處理；②燃燒化學(xué)反應(yīng)特征時(shí)間相對(duì)流動(dòng)特征時(shí)間較小，因此快速反應(yīng)假設(shè)成立且流動(dòng)化學(xué)反應(yīng)可以近似解耦計(jì)算；③流速遠(yuǎn)小于聲速，壓力波近似以無(wú)限快速度傳播，無(wú)黏流動(dòng)部分呈現(xiàn)橢圓型流動(dòng)特征．針對(duì)上述流動(dòng)特點(diǎn)，低馬赫燃燒求解可以采取基于壓力驅(qū)動(dòng)速度原理的壓力求解器，并且基于氣體狀態(tài)方程耦合燃燒化學(xué)反應(yīng)效應(yīng)實(shí)現(xiàn)流動(dòng)和燃燒的隱式求解以降低對(duì)時(shí)間步的CFL(Courant-Friedrichs-Lewy)約束(CFL數(shù)＝速度×?xí)r間步/網(wǎng)格尺度，一般要求此數(shù)小于某一定值以保證數(shù)值格式的穩(wěn)定性)．在復(fù)雜燃燒流場(chǎng)中由于整場(chǎng)流速的不均勻性和局部網(wǎng)格尺寸的差異，各局部流場(chǎng)區(qū)域的CFL數(shù)差異較大．特別是在模擬復(fù)雜幾何結(jié)構(gòu)或需要捕捉局部精細(xì)流場(chǎng)結(jié)構(gòu)(如近壁邊界層或火焰反應(yīng)區(qū))時(shí)通常需要使用非均勻或局部自適應(yīng)網(wǎng)格，根據(jù)CFL約束時(shí)間步正比于網(wǎng)格尺度．傳統(tǒng)方法中整場(chǎng)使用單一時(shí)間步，通常使得時(shí)間步受限于整場(chǎng)最小網(wǎng)格尺度，因而顯著降低了計(jì)算效率．而如果能夠針對(duì)局部流場(chǎng)采用多重時(shí)間步，將可在滿足CFL約束的前提下顯著提高計(jì)算效率．

基于局部時(shí)間步(local time-stepping)的局部時(shí)間步法由Osher與Sanders[1]在求解一維標(biāo)量守恒律方程時(shí)提出，其主要原理為在滿足局部穩(wěn)定性條件下，每個(gè)網(wǎng)格單元擁有獨(dú)立的局部時(shí)間步長(zhǎng)并可在計(jì)算過(guò)程中自適應(yīng)更新．相比于所有網(wǎng)格單元按統(tǒng)一時(shí)間層更新的方法，一方面提高了計(jì)算效率，另一方面也加速了流場(chǎng)收斂．空間差異的多重時(shí)間法等效于對(duì)控制方程的空間預(yù)處理，可以一定程度上消除空間梯度較大區(qū)域的剛性[2]．多重時(shí)間法成功應(yīng)用于基于笛卡爾網(wǎng)格的自適應(yīng)網(wǎng)格加密(AMR)算例[3-4]，其通過(guò)解耦非加密區(qū)域與局部加密區(qū)域的時(shí)間步有效提高了計(jì)算效率并且實(shí)現(xiàn)了對(duì)加密區(qū)域流場(chǎng)結(jié)構(gòu)(如激波)的準(zhǔn)確捕捉．Dumbser等[5]和Taube等[6]在涉及網(wǎng)格尺度劇烈變化的對(duì)彈性波和Maxwell方程的任意高階間斷有限元(DG-FEM)模擬中耦合了顯式局部時(shí)間步法以提高計(jì)算效率和收斂速率．Qi等[7]通過(guò)耦合間斷有限元和單步高階泰勒式時(shí)間積分格式實(shí)現(xiàn)了對(duì)局部時(shí)間步法在基于多尺度網(wǎng)格的瞬態(tài)電磁波模擬中的支持．Grote等[8]推導(dǎo)了局部時(shí)間步法與任意高階顯式Runge-Kutta方法的耦合形式并成功應(yīng)用于多尺度網(wǎng)格波方程求解. Lilia[9]在基于局部時(shí)間步的二階Runge-Kutta格式中構(gòu)造了更加準(zhǔn)確的二次多項(xiàng)式以計(jì)算時(shí)間步相異的相鄰網(wǎng)格界面值．Ashbourne[10]進(jìn)一步將Lilia[9]的方法拓展到三階和四階Runge-Kutta格式并改進(jìn)了網(wǎng)格界面插值的多項(xiàng)式逼近方法．杜永樂(lè)等[11]和姜婷婷等[12]在稀薄氣體的DSMC模擬中采用了類似的局部自適應(yīng)時(shí)間步方法表明其可以在保真的前提下顯著提高流場(chǎng)收斂速率．譚志軍等[13]基于局部時(shí)間步法和自適應(yīng)網(wǎng)格加密方法對(duì)一維對(duì)流擴(kuò)散方程所做的測(cè)試表明實(shí)現(xiàn)了2～3倍的加速，同時(shí)未見(jiàn)明顯的精度降低．吳迪等[14]基于局部時(shí)間步法在超聲速后臺(tái)階流的三維歐拉方程模擬中相比于整體時(shí)間步長(zhǎng)法最高取得了5倍的加速比．胡鵬等[15]基于局部時(shí)間步長(zhǎng)的水沙耦合模擬表明，該方法能夠在保證精度的前提下節(jié)省高達(dá)92%的計(jì)算時(shí)間．Espinoza等[16]在基于局部時(shí)間步法的可壓縮高速和超聲速穩(wěn)態(tài)流動(dòng)模擬中分別觀測(cè)到了2.56和8.96倍的迭代收斂速度．最近，除了基于當(dāng)?shù)谻FL數(shù)控制局部時(shí)間步以外，Kalkote等[17]提出了一種基于當(dāng)?shù)亟財(cái)嗾`差控制的局部自適應(yīng)時(shí)間步方法．Jeanmasson等[18]則對(duì)局部時(shí)間步法進(jìn)一步發(fā)展了通量守恒修正．

目前基于局部時(shí)間步的局部時(shí)間步法在湍流燃燒中的應(yīng)用還不多見(jiàn)．而湍流燃燒通常需要對(duì)主要釋熱的反應(yīng)區(qū)域重點(diǎn)求解，意味著需要采用較密的局部網(wǎng)格，而對(duì)其余的無(wú)反應(yīng)區(qū)域可以采用相對(duì)較粗的網(wǎng)格．特別是在火災(zāi)模擬中，由于模擬場(chǎng)景相對(duì)火源區(qū)域一般較大，為了提高計(jì)算效率一般僅對(duì)火源區(qū)域局部加密．另外在涉及復(fù)雜幾何結(jié)構(gòu)的燃燒室或建筑火災(zāi)模擬中，對(duì)于曲率較大的結(jié)合拐角需要局部加密以逼近幾何外形．燃燒和火災(zāi)一般發(fā)生于受限空間內(nèi)，如發(fā)動(dòng)機(jī)和室內(nèi)，為了準(zhǔn)確模擬近壁效應(yīng)一般需要對(duì)壁面邊界層區(qū)域采用附加膨脹層的方法局部加密．因此針對(duì)湍流燃燒模擬，如果采用整體時(shí)間層統(tǒng)一的時(shí)間步，將使得時(shí)間推進(jìn)受限于火源、拐角和近壁等區(qū)域處的最小時(shí)間步．而通常局部加密區(qū)域僅占整體計(jì)算域較小的體積分?jǐn)?shù)．按照通常火災(zāi)模擬中火源、拐角和近壁等區(qū)域局部加密設(shè)置，最大最小網(wǎng)格尺度一般相差3～5倍來(lái)估算，采用局部時(shí)間步法預(yù)計(jì)至少將可使計(jì)算效率提高3～5倍．在并行環(huán)境下，位于火源局部加密區(qū)域的分區(qū)節(jié)點(diǎn)因需求解化學(xué)反應(yīng)計(jì)算載荷通常較大，而其他對(duì)應(yīng)無(wú)反應(yīng)區(qū)的節(jié)點(diǎn)計(jì)算負(fù)擔(dān)較小，從而導(dǎo)致計(jì)算載荷的嚴(yán)重不均衡．而采用局部時(shí)間步法，對(duì)應(yīng)無(wú)反應(yīng)區(qū)的計(jì)算節(jié)點(diǎn)可將額外的計(jì)算資源用于時(shí)間上的加速推進(jìn)(大時(shí)間步一般意味著更多的迭代次數(shù)和收斂時(shí)間)，因而也可以在一定程度上緩解計(jì)算載荷的不均衡性．

火災(zāi)和內(nèi)燃等面向?qū)嶋H應(yīng)用的低馬赫數(shù)燃燒模擬的兩大主要特點(diǎn)為：(1)模擬工況尺寸較大且內(nèi)部幾何結(jié)構(gòu)復(fù)雜，計(jì)算區(qū)域尺寸可跨越數(shù)十米甚至上百米；(2)完整的燃燒起始、發(fā)展和熄火過(guò)程模擬需要求解較長(zhǎng)的時(shí)間歷程，對(duì)應(yīng)物理時(shí)間可達(dá)數(shù)分鐘(約600s)甚至1h．例如大型建筑火災(zāi)[19-20]、隧道火??災(zāi)[21-22]、艦船火災(zāi)[23]、工業(yè)燃燒鍋爐[24-26]模擬等．這些低馬赫數(shù)燃燒場(chǎng)景模擬一方面對(duì)計(jì)算效率有特殊的實(shí)際需求，另一方面又存在需要著重求解的區(qū)域(如火源、疏散口、重點(diǎn)熱防護(hù)部件等)．模擬中需要兼顧計(jì)算效率與求解精細(xì)度，而對(duì)最終結(jié)果一般僅關(guān)注重點(diǎn)區(qū)域的穩(wěn)態(tài)流場(chǎng)．傳統(tǒng)的燃燒模擬程序側(cè)重于對(duì)小尺寸(如射流火焰)或局部燃燒場(chǎng)的高精度、高保真模擬，所采用的時(shí)空高階格式一方面顯著增加了計(jì)算代價(jià)，另一方面也削弱了求解復(fù)雜流場(chǎng)的魯棒性，在工程燃燒模擬中的適用性較弱．

本研究針對(duì)火災(zāi)等低馬赫數(shù)燃燒模擬需求開(kāi)發(fā)了基于局部時(shí)間步法的壓力校正隱式燃燒求解器，并將其應(yīng)用于開(kāi)放空間池火和封閉空間建筑火災(zāi)模擬中以進(jìn)一步驗(yàn)證其準(zhǔn)確性和計(jì)算效率．本文將首先介紹局部時(shí)間步法的實(shí)現(xiàn)算法、對(duì)PaSR湍流燃燒模型的改進(jìn)以及對(duì)壓力求解算法的穩(wěn)定性優(yōu)化等技術(shù)細(xì)節(jié)，進(jìn)而結(jié)合具體算例分析局部時(shí)間步法與整體時(shí)間步法對(duì)湍流燃燒時(shí)均流場(chǎng)的量化影響，以考察局部時(shí)間步法在低馬赫數(shù)湍流燃燒中的適用性．

1?數(shù)學(xué)和物理模型

1.1?控制方程

1.2?改進(jìn)型PaSR湍流燃燒模型

1.3?基于魯棒增強(qiáng)型PIMPLE算法的隱式壓力求解

1.4?基于局部時(shí)間步的高效多重時(shí)間法

1.5?基于OpenFOAM的求解器與算例設(shè)置

本求解器開(kāi)發(fā)的主要目的是為火災(zāi)等低馬赫數(shù)準(zhǔn)穩(wěn)態(tài)燃燒模擬提供高效、準(zhǔn)確且魯棒性較好的計(jì)算平臺(tái)．為提高算例的普適性，擬分別針對(duì)開(kāi)放和封閉空間兩種典型燃燒情形進(jìn)行計(jì)算驗(yàn)證．為提高算例的標(biāo)準(zhǔn)性，算例選取為FireFOAM自帶的標(biāo)準(zhǔn)開(kāi)放空間池火和封閉空間單室建筑火災(zāi)算例．其中開(kāi)放池火為理想算例，研究主要比較兩種方法預(yù)測(cè)的差異；而單室建筑火災(zāi)為實(shí)驗(yàn)測(cè)試算例，研究結(jié)合實(shí)驗(yàn)數(shù)據(jù)對(duì)兩種方法預(yù)測(cè)的準(zhǔn)確性進(jìn)行了驗(yàn)證．在以下計(jì)算中，整體時(shí)間步算例采用原FireFOAM求解器(不支持多重時(shí)間步功能)計(jì)算，多重時(shí)間步則采用改進(jìn)的求解器(命名為Amber)進(jìn)行．本研究主要目的為驗(yàn)證局部時(shí)間步法與傳統(tǒng)的整體時(shí)間步法在預(yù)測(cè)關(guān)鍵燃燒場(chǎng)特性(如溫度、速度、組分濃度等)的準(zhǔn)確性與計(jì)算效率．物理模型和數(shù)值方法對(duì)燃燒物理過(guò)程的影響分析不是本研究的重點(diǎn)，擬在后續(xù)研究中開(kāi)展．

2?結(jié)果與討論

2.1?開(kāi)放空間內(nèi)池火模擬

該算例模擬了一個(gè)典型的開(kāi)放空間內(nèi)池火工況．計(jì)算域?yàn)?m×1m×1m的立方形空間，在底部中央有一0.2m×0.2m的方形火源區(qū)域．燃料為氣態(tài)甲烷，常溫(300K)常壓(101325Pa)下以法向速度0.01m/s噴入并發(fā)生燃燒．模擬中考慮重力效應(yīng)，火焰受浮力誘導(dǎo)的自然對(duì)流驅(qū)動(dòng)．由于計(jì)算區(qū)域的規(guī)則性，網(wǎng)格多為正六面體正交網(wǎng)格，整體區(qū)域網(wǎng)格尺度為0.5mm，由于整場(chǎng)網(wǎng)格尺度已經(jīng)滿足大渦模擬要求(由網(wǎng)格尺度決定的截?cái)嗖〝?shù)位于慣性子區(qū))，火源區(qū)域不再局部加密，整體網(wǎng)格數(shù)目為800萬(wàn)．該算例總模擬時(shí)間為100s，以特征速度1m/s計(jì)算約為100個(gè)整場(chǎng)流通時(shí)間(FTT)，其單FTT時(shí)間為0＝1s．圖1(a)比較了不同網(wǎng)格條件下的中心軸線速度分布，可見(jiàn)隨著網(wǎng)格密度增加，預(yù)測(cè)結(jié)果趨于相同，并且最終與整體時(shí)間步法的預(yù)測(cè)接近．這里定義誤差為與最密網(wǎng)格(1065萬(wàn))預(yù)測(cè)的相對(duì)誤差．由圖1(b)可見(jiàn)，平均誤差與最大誤差均隨網(wǎng)格尺度減小而減小，平均誤差隨網(wǎng)格尺度近似呈指數(shù)衰減(對(duì)數(shù)空間表現(xiàn)為線性)．在網(wǎng)格密度大于800萬(wàn)的情況下，平均誤差小于0.5%，最大誤差小于2%，滿足網(wǎng)格獨(dú)立性要求．以下分析如不特殊指明，均基于1065萬(wàn)網(wǎng)格結(jié)果．

圖1?局部時(shí)間步法的網(wǎng)格收斂性分析

圖2?火焰中心剖面局步時(shí)間步分布

圖3 無(wú)量綱計(jì)算周期和加速比對(duì)網(wǎng)格尺寸的依賴性關(guān)系

圖5定量比較了中心軸線的時(shí)均值．火焰高度(＝0.6m)以下的時(shí)均值分布幾乎相同，之后差異越往下游越大，但整體差異仍較?。捎镁植繒r(shí)間步法，溫度和速度峰值與采用整體時(shí)間步法預(yù)測(cè)基本一致，但平均氧氣消耗率略低，火焰溫度和最終浮力驅(qū)動(dòng)速度也略低于整體時(shí)間步法的預(yù)測(cè)．在中心軸線上各關(guān)鍵變量的平均差異分別為1.3%(溫度)，0.7%(速度)，3.9%(CH4)，3.7%(O2)，均小于5%，可以近似認(rèn)為采用局部時(shí)間法和整體時(shí)間法對(duì)時(shí)均結(jié)果的預(yù)測(cè)一致．

多重時(shí)間步法等效于對(duì)時(shí)間項(xiàng)的權(quán)重預(yù)處理，有限體積方法中，相鄰網(wǎng)格單元交界面的質(zhì)量通量計(jì)算與時(shí)間無(wú)關(guān)，根據(jù)單元內(nèi)質(zhì)量、動(dòng)量、組分和能量等的守恒原則時(shí)間權(quán)重影響各網(wǎng)格單元值更新的快慢．時(shí)間步較大的區(qū)域快速演化迅速趨于準(zhǔn)穩(wěn)態(tài)，而時(shí)間步較小的區(qū)域各變量相當(dāng)于處于相對(duì)“凍結(jié)”狀態(tài)．不同區(qū)域交界處互為邊界條件，意味著時(shí)間演化較快的區(qū)域通過(guò)提供一個(gè)快速趨于準(zhǔn)穩(wěn)態(tài)的邊界加速了其相鄰時(shí)間演化較慢區(qū)域的收斂．本算例中，以射流邊界為區(qū)分，可將火焰流場(chǎng)分為演化較慢的火焰中心區(qū)域和快速趨于穩(wěn)態(tài)的伴流．伴流通過(guò)為中心射流提供了更趨于穩(wěn)態(tài)的射流邊界條件，加速了中心火焰區(qū)的收斂，類似于Espinoza等[16]在可壓縮高速流動(dòng)中的觀察．因而采用局部時(shí)間步法從兩個(gè)方面促進(jìn)了準(zhǔn)穩(wěn)態(tài)流場(chǎng)的建立：一方面通過(guò)增加時(shí)間步加速了無(wú)反應(yīng)區(qū)域的演化，另一方面通過(guò)邊界效應(yīng)促進(jìn)了相對(duì)演化較慢的火焰區(qū)域往準(zhǔn)穩(wěn)態(tài)收斂的進(jìn)度.

圖5?中心軸線上時(shí)均值比較

2.2?封閉空間內(nèi)池火模擬

計(jì)算區(qū)域?yàn)榈湫偷膬墒医ㄖ?，其中左邊房間尺寸為0.4m(長(zhǎng))×0.4m(寬)×0.6m(高)，右邊緊鄰房間的尺寸為0.4m(長(zhǎng))×0.4m(寬)×0.4m(高)．直徑0.095m的圓形火源位于右側(cè)房間底部中央．在兩個(gè)相鄰房間由一堵厚度為25.4cm的墻隔開(kāi)，墻的底部和頂部各留有30cm的間隙分別用于通風(fēng)和排煙．整體網(wǎng)格數(shù)目為320萬(wàn)，其中圍繞火源區(qū)域的0.2m(長(zhǎng))×0.2m(寬)×0.3m(高)的區(qū)域局部加密網(wǎng)格．左側(cè)房間頂部和側(cè)邊為開(kāi)放邊界條件，固定壓力為環(huán)境大氣壓(101325Pa)．右側(cè)房間為完全封閉空間，僅通過(guò)通風(fēng)口和排煙孔和外界平衡壓力．

在火焰的發(fā)展演化過(guò)程中，新鮮空氣從左側(cè)吹入帶來(lái)的動(dòng)量使火焰偏向右方．同時(shí)由于左側(cè)頂部的抽吸，熱煙氣層向左上方移動(dòng)．相比較而言，右側(cè)房間中靠近右壁面一側(cè)的氣流速度較大，相應(yīng)地，圖6中的局部時(shí)間步右側(cè)較?。髠?cè)房間通過(guò)快速時(shí)間演化相當(dāng)于為右側(cè)房間的流場(chǎng)發(fā)展提供了一個(gè)穩(wěn)定的邊界條件，而右側(cè)房間內(nèi)的計(jì)算域時(shí)間步較小，數(shù)值耗散更低，有利于捕捉到更加豐富的火焰瞬態(tài)效應(yīng).這也意味著，通過(guò)采用多重時(shí)間步，對(duì)著重關(guān)注的瞬態(tài)流動(dòng)區(qū)域指定較小的時(shí)間步，在對(duì)物理過(guò)程精細(xì)捕捉的同時(shí)兼顧計(jì)算效率．

圖6?時(shí)間步分布

圖7比較了采用局部時(shí)間步法和整體時(shí)間步法的瞬態(tài)溫度和O2體積分?jǐn)?shù)預(yù)測(cè)結(jié)果．總體上兩種方法的預(yù)測(cè)與實(shí)驗(yàn)測(cè)量值的趨勢(shì)均吻合較好．局部時(shí)間步法可以比整體時(shí)間步法更快地收斂到穩(wěn)態(tài)解[16]，因此局部時(shí)間步法更適用于穩(wěn)態(tài)流場(chǎng)模擬．對(duì)于低馬赫數(shù)燃燒問(wèn)題流動(dòng)速度遠(yuǎn)小于聲速，壓力波以近乎無(wú)限快的速度在整場(chǎng)傳播，驅(qū)動(dòng)速度場(chǎng)更快地達(dá)到平衡，并且流動(dòng)特征時(shí)間遠(yuǎn)大于化學(xué)反應(yīng)特征時(shí)間，流動(dòng)和燃燒均趨于準(zhǔn)穩(wěn)態(tài)，因此局部時(shí)間步法較為適用．局部時(shí)間步方法更適用于最終達(dá)到穩(wěn)態(tài)或準(zhǔn)穩(wěn)態(tài)的流場(chǎng)模擬，通過(guò)快速演化區(qū)域與慢演化區(qū)域的互為邊界原理實(shí)現(xiàn)了整體的加速收斂．由于時(shí)間步較大區(qū)域的單步時(shí)間內(nèi)變量增幅較大，相應(yīng)地會(huì)表現(xiàn)出較大的脈動(dòng)值，如圖中所示．在火源上方0.14m和0.26m高度位置整體時(shí)間步法的預(yù)測(cè)更貼近實(shí)驗(yàn)值，0.38m處前半段整體時(shí)間步法預(yù)測(cè)較好而后半段局部時(shí)間步預(yù)測(cè)更貼近測(cè)量值．在初始點(diǎn)火(0～40s)和熄火階段(310～350s)，由于釋熱速率的急劇變化導(dǎo)致瞬態(tài)效應(yīng)較強(qiáng)，兩種方法的預(yù)測(cè)均存在一定偏差．靠近火源根部(＝0.14m高度處)燃料和空氣摻混后反應(yīng)劇烈，表現(xiàn)在測(cè)量和預(yù)測(cè)上溫度空間波動(dòng)均較大，而高處羽流區(qū)反應(yīng)微弱因而溫度波動(dòng)更為平穩(wěn)．圖7(d)中兩種方法對(duì)氧氣濃度的預(yù)測(cè)在中間準(zhǔn)穩(wěn)態(tài)階段均與實(shí)驗(yàn)測(cè)量值吻合較好，其中由于時(shí)間步的增加，局部時(shí)間步法的預(yù)測(cè)表現(xiàn)出較大的脈動(dòng).

圖8比較了穩(wěn)定燃燒階段(?。?00s處燃燒速率)基于局部時(shí)間步法和整體時(shí)間步法預(yù)測(cè)的時(shí)均值分布．為方便對(duì)比，在圖例范圍內(nèi)標(biāo)識(shí)了5條等間距輪廓線．整體而言，兩種方法預(yù)測(cè)的時(shí)均流場(chǎng)相似度較高．對(duì)于時(shí)均溫度，采用局部時(shí)間法預(yù)測(cè)的火焰高溫核心區(qū)域略寬而高度相似．時(shí)均速度未見(jiàn)明顯差異．從溫度和速度分布可以明顯看出，由于左下部空氣的通入，使得火焰向右側(cè)傾斜明顯．又由于左上部的抽吸作用，熱羽流向左側(cè)移動(dòng)．整體火羽流路徑呈現(xiàn)反C字型．從產(chǎn)物H2O和CO2的分布看，熱燃燒產(chǎn)物在浮力作用下上升到頂棚以后，左側(cè)的部分被排煙通道抽吸走，右側(cè)的部分被右側(cè)倉(cāng)壁阻擋后沿著壁面向下堆積．與溫度的分布類似，采用局部時(shí)間步法預(yù)測(cè)的核心羽流區(qū)(高溫高產(chǎn)物濃度輪廓線標(biāo)識(shí))的寬度稍窄．在火災(zāi)發(fā)展過(guò)程中，左房間為半開(kāi)放空間(上、左、前、后均為自由開(kāi)放邊界條件)，因而左側(cè)房間的流場(chǎng)趨于穩(wěn)態(tài)，兩種方法的預(yù)測(cè)更為相似．

圖8?火源中心剖面時(shí)均溫度、速度、CO2和H2O質(zhì)量分?jǐn)?shù)分布

圖9量化比較了圖8中從火源到頂棚(＝0.1～0.39m)不同高度處橫向的時(shí)均量分布．對(duì)于燃料C7H16，在＝0.1m靠近火源底部濃度較大兩種方法預(yù)測(cè)的偏差較小．隨著高度增加，燃料濃度逐漸降低至千分之一量級(jí)，微量痕跡物的預(yù)測(cè)相對(duì)誤差增加．對(duì)于氧氣質(zhì)量分?jǐn)?shù)，兩種方法的預(yù)測(cè)在不同高度處量級(jí)和趨勢(shì)均較為相似．對(duì)于溫度分布二者的預(yù)測(cè)幾乎沒(méi)有差別，尤其是對(duì)于火源區(qū)域的溫度峰值和近壁區(qū)域的溫度陡降均吻合較好．在反C型流場(chǎng)的中部，由于底部進(jìn)氣道向右氣流動(dòng)量和火羽流向上浮力的驅(qū)動(dòng)會(huì)形成一回流區(qū)．速度分布顯示該回流區(qū)底部(0.1m)進(jìn)氣口位置在采用局部時(shí)間法時(shí)略低于整體時(shí)間步法．0.3m高度處局部時(shí)間步法預(yù)測(cè)的速度峰值略低，這也對(duì)應(yīng)圖8中較窄的高溫區(qū)域(浮力較弱)．各時(shí)均量在4個(gè)高度處的平均相對(duì)誤差分別為1.1%(溫度，0.1～0.39m)、6.1%(速度，0.1～0.39m)、1.0%(氧氣質(zhì)量分?jǐn)?shù)，0.1～0.39m)、5.0% (C7H16濃度，0.1～0.2m高度)，＞100%(C7H16濃度，0.3～0.39m高度)，可見(jiàn)除了燃料因濃度衰減較快導(dǎo)致遠(yuǎn)離火焰區(qū)相對(duì)誤差較大以外，其他量的預(yù)測(cè)誤差均較?。捎镁植繒r(shí)間步法，增加了弱流動(dòng)區(qū)域的時(shí)間步，因而一定程度上增加了數(shù)值耗散性，這有可能是引起上述差異的原因，然而整體上二者的時(shí)均結(jié)果十分相近，表明了局部時(shí)間步法在低馬赫數(shù)準(zhǔn)穩(wěn)態(tài)燃燒問(wèn)題中的適用性，以及犧牲非主要燃燒反應(yīng)區(qū)微量痕跡物的預(yù)測(cè)精度以顯著提高整體計(jì)算效率是可行的．

3?結(jié)?語(yǔ)

本研究基于局部時(shí)間步法、改進(jìn)PaSR湍流燃燒模型以及優(yōu)化壓力求解算法開(kāi)發(fā)了面向低馬赫數(shù)的壓力求解器，并選取了開(kāi)放空間和封閉空間兩種典型的低馬赫數(shù)燃燒算例進(jìn)行了驗(yàn)證．一般可認(rèn)為，局部時(shí)間步法的加速比近似與網(wǎng)格差異化程度成正比．按照通常的火源、拐角、近壁等區(qū)域局部加密設(shè)置，最大/最小網(wǎng)格尺度一般至少相差3～5倍，采用局部時(shí)間步法預(yù)計(jì)至少可使計(jì)算效率提高3～5倍．本文中的網(wǎng)格設(shè)置使得局部時(shí)間步法相比于整體時(shí)間步法分別實(shí)現(xiàn)了6倍和8倍的實(shí)際加速比．實(shí)際加速比通常略小于由平均時(shí)間步之比定義的理想值，這是因?yàn)槌藭r(shí)間步計(jì)算本身的額外代價(jià)外，輻射和剛性化學(xué)反應(yīng)求解一般耗時(shí)較多且與流動(dòng)時(shí)間步關(guān)系較?。傮w上計(jì)算周期隨網(wǎng)格尺寸減小近似呈線性增加，其中整體時(shí)間步法計(jì)算周期的上升梯度更大，表明密網(wǎng)格條件下局部時(shí)間步法有更好的加速效果．與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比表明，局部和整體時(shí)間步法均較為準(zhǔn)確地預(yù)測(cè)了溫度的時(shí)間變化特性和除微量痕跡物以外的時(shí)均流場(chǎng)特性，表明了局部時(shí)間步法在低馬赫數(shù)準(zhǔn)穩(wěn)態(tài)燃燒問(wèn)題中的適用性，以及犧牲非主要燃燒反應(yīng)區(qū)微量痕跡物的預(yù)測(cè)精度以顯著提高整體計(jì)算效率是可行的．

一般而言，局部時(shí)間步法針對(duì)穩(wěn)態(tài)流動(dòng)或統(tǒng)計(jì)上呈現(xiàn)穩(wěn)態(tài)的準(zhǔn)穩(wěn)態(tài)流動(dòng)問(wèn)題，其得到的結(jié)果是比較可靠的，但對(duì)于強(qiáng)瞬態(tài)的流場(chǎng)演化問(wèn)題，其方法的適用性則需要進(jìn)一步的研究探討．

［1］ Osher S，Sanders R. Numerical approximations to nonlinear conservation laws with locally varying time and space grids[J].，1983，41(164)：321-336.

［2］ Powell K G，van Leer B. A genuinely multi-dimensional upwind cell-vertex scheme for the Euler equations[C]//. Reno，USA，1989.

［3］ Berger M J，Oliger J. Adaptive mesh refinement for hyperbolic partial differential equations[J].，1984，53(3)：484-512.

［4］ Berger M J，Colella P. Local adaptive mesh refinement for shock hydrodynamics[J].，1989，82(1)：64-84.

［5］ Dumbser M，K?ser M，Toro E F. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes—V. Local time stepping and-adaptivity[J].，2007，171(2)：695-717.

［6］ Taube A，Dumbser M，Munz C D，et al. A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations[J].：，，2009，22(1)：77-103.

［7］ Qi H，Wang X，Zhang J，et al. An ADER discontinuous Galerkin method with local time-stepping for transient electromagnetics[J].，2018，229：106-115.

［8］ Grote M J，Mehlin M，Mitkova T. Runge-Kutta-based explicit local time-stepping methods for wave propagation[J].，2015，37(2)：A747-A775.

［9］ Lilia Krivodonova. An efficient local time-stepping scheme for solution of nonlinear conservation laws[J].，2010，229(22)：8537-8551.

［10］ Ashbourne A. Efficient Runge-Kutta Based Local Time-Stepping Methods[D]. Waterloo，Canada：University of Waterloo，2016.

［11］ 杜永樂(lè)，閻?超. DSMC方法中的自適應(yīng)當(dāng)?shù)貢r(shí)間步長(zhǎng)法[J]. 北京航空航天大學(xué)學(xué)報(bào)，2006，32(4)：388-390.

Du Yongle，Yan Chao. Adaptive local time step method for DSMC code [J].，2006，32(4)：388-390(in Chinese).

［12］ 姜婷婷，錢(qián)?坤，陳偉芳. 高超聲速混合流動(dòng)的粒子模擬耦合算法[J]. 固體火箭技術(shù)，2017，40(3)：283-288.

Jiang Tingting，Qian Kun，Chen Weifang. Study on hybrid particle simulation method for hypersonic mixed flows[J].，2017，40(3)：283-288(in Chinese).

［13］ 譚志軍，黃云清. 對(duì)一維守恒律的一種局部時(shí)間步長(zhǎng)自適應(yīng)網(wǎng)格方法[J]. 湘潭大學(xué)自然科學(xué)學(xué)報(bào)，2003，25(2)：110-116.

Tan Zhijun，Huang Yunqing. An adaptive grid method with local time stepping for one-dimensional conservation laws[J].，2003，25(2)：110-116(in Chinese).

［14］ 吳?迪，蔚喜軍，徐?云. 局部時(shí)間步長(zhǎng)間斷有限元方法求解三維歐拉方程[J]. 計(jì)算物理，2011，28(1)：1-9.

Wu Di，Wei Xijun，Xu Yun. A discontinuous Galerkin method with local time stepping for Euler equations[J].，2011，28(1)：1-9(in Chinese).

［15］ 胡?鵬，韓健健，雷云龍. 基于局部分級(jí)時(shí)間步長(zhǎng)方法的水沙耦合數(shù)學(xué)模擬[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版)，2019，53(4)：743-752.

Hu Peng，Han Jianjian，Lei Yunlong. Coupled modeling of sediment-laden flows based on local-time-step approach[J].()，2019，53(4)：743-752(in Chinese).

［16］ Espinoza D E R，Scanlon T J，Brown R E. Validation of tools to accelerate high-speed CFD simulations using OpenFOAM[C]// 20. Glasgow，UK，2015.

［17］ Kalkote N，Assam A，Eswaran V. Towards developing an adaptive time stepping for compressible unsteady flows[J].，2019，29(2)：487-503.

［18］ Jeanmasson G，Mary I，Mieussens L. On some explicit local time stepping finite volume schemes for CFD[J].，2019，397：108818.

［19］ Cant R S. Modelling and simulation of accidental combustion in complex geometries[C]//. Edinburgh，UK，2007.

［20］ 姜?蓬，邱?榕，蔣?勇. 基于數(shù)值模擬的某大廈特大火災(zāi)過(guò)程調(diào)查[J]. 燃燒科學(xué)與技術(shù)，2007，13(1)：76-80.

Jiang Peng，Qiu Rong，Jiang Yong. Investigation of a mansion fire based on numerical simulation[J].，2007，13(1)：76-80(in Chinese).

［21］ Wang H Y. Prediction of soot and carbon monoxide production in a ventilated tunnel fire by using a computer simulation[J].，2009，44(3)：394-406.

［22］ Ji J，Guo F，Gao Z，et al. Numerical investigation on the effect of ambient pressure on smoke movement and temperature distribution in tunnel fires[J].，2017，118：663-669.

［23］ Su S，Wang L. Three dimensional reconstruction of the fire in a ship engine room with multilayer structures[J].，2013，70：201-207.

［24］ 王秋紅，王?超，張小桃，等. 四角切圓鍋爐爐內(nèi)燃燒數(shù)值模擬[J]. 熱力發(fā)電，2016，45(9)：61-66.

Wang Qiuhong，Wang Chao，Zhang Xiaotao，et al. Numerical simulation of combustion process in a tangentially pulverized-coal-fired boiler under variable load operation[J].，2016，45(9)：61-66(in Chinese).

［25］ 張?成，朱天宇，殷立寶，等. 100MW燃煤鍋爐污泥摻燒試驗(yàn)與數(shù)值模擬[J]. 燃燒科學(xué)與技術(shù)，2015，21(2)：114-123.

Zhang Cheng，Zhu Tianyu，Yin Libao，et al. Field test and numerical simulation for co-combustion of sludge in a 100 MW coal fired boiler[J].，2015，21(2)：114-123(in Chinese).

［26］ 馬?侖，方慶艷，張?成，等. 600MW W型火焰鍋爐拱上二次風(fēng)低NO燃燒特性的數(shù)值模擬及優(yōu)化[J]. 燃燒科學(xué)與技術(shù)，2016，22(1)：64-70.

Ma Lun，F(xiàn)ang Qingyan，Zhang Cheng，et al. Numerical simulation and optimization for the influence of the arch secondary airon low-NOcombustion characteristics of a 600MW down-fired boiler[J].，2016，22(1)：64-70(in Chinese).

［27］ Yoshizawa A. Statistical theory for compressible turbulent shear flows，with the application to subgrid modeling[J].，1986，29(7)：2152-2164.

［28］ Golovitchev V I，Nordin N，Jarnicki R，et al. 3-D diesel spray simulations using a new detailed chemistry turbulent combustion model[C]//. Paris，F(xiàn)rance，2000，2000-01-1891.

［29］ Kjaldman L，Brink A，Hupa M. Micro mixing time in the eddy dissipation concept[J].，2000，154(1)：207-227.

［30］ Chomiak J，Karlsson A. Physical and chemical effects in diesel spray ignition[C]//21. Interlaken，Swizerland，1995.

［31］ Fureby C. LES for supersonic combustion[C]//. Tours，F(xiàn)rance，2012，2012-5979.

［32］ Chomiak J，Karlsson A. Flame liftoff in diesel sprays[J].()，1996，26(2)：2557-2564.

［33］ Hallaji M，Mazaheri K. Numerical simulation of turbulent non-premixed combustion in diluted hot coflow using PaSr combustion model[C]//. Sardinia，Italy，2011.

［34］ Yao W，Lu Y，Wu K，et al. Modeling analysis of an actively cooled scramjet combustor under different kerosene/air ratios[J].，2018，34(4)：975-991.

［35］ John D A J. Inviscid high-temperature nonequilibrium flows.[M]. 2nd Edition. Virginia，USA：AIAA，2006.

［36］ Weller H G，Tabor G，Jasak H，et al. A tensorial approach to computational continuum mechanics using object oriented techniques[J].，1997，12(6)：620-631.

［37］ Murthy J Y，Mathur S R. Finite volume method for radiative heat transfer using unstructured meshes[J].，1998，12(3)：313-321.

Study on Applicability of Local Time-Stepping Method in Low-Mach Combustion Simulation

Yao Wei1, 2，Sun Chao3，Liu Hang1，Wu Mei4，Yang Shaobo3，Cao Shunli3

(1. Key Laboratory of High-Temperature Gas Dynamics，Institute of Mechanics，Chinese Academy of Sciences，Beijing 100190，China；2. School of Engineering Science，University of Chinese Academy of Sciences，Beijing 100049，China；3. China Shipbuilding Industry Systems Engineering Research Institute，Beijing 100094，China；4. College of Civil Engineering，Hefei University of Technology，Hefei 230009，China)

Due to the non-uniformity of flow speed in the whole field and the difference in local grid size in complicated combustion simulations，the local Courant-Friedrichs-Lewy(CFL)number varies significantly among different local flow field regions. The computational efficiency is also severely restricted by the conventional global time-stepping method，which is on the basis of the definition of maximum CFL number. In this paper，the applicability of the local time-stepping method restricted by the local CFL number to low-Mach turbulent combustion simulation is investigated for the first time. Through the large-eddy simulation of an open methane pool fire(with 10.65 million cells)and an enclosed building fire(with 3.2 million cells)，it is shown that the local time-stepping method can achieve speedup ratios of 6 and 8 compared with the global time-stepping method，respectively. The speedup ratio rises with the reduction of grid size. Moreover，the accuracy of the local time-stepping method applied in the low-Mach combustion simulation is validated from two aspects. First，in comparison with the experimental data，the time-varying characteristics of temperature are accurately predicted by both the local and global methods，which is probably due to the quasi-steadiness of low-Mach combustion. Second，in comparison with the time-averaged flow field，the differences in time-averaged temperature，time-averaged velocity，and the volume fraction of oxygen are relatively smaller except for some minor trace species(with a mass fraction of less than 0.1%). In addition，the existing PaSR turbulent combustion model and the pressure-momentum coupling algorithm are also improved and optimized，thereby enhancing their physical accuracy and robustness.

turbulent combustion；local time-stepping method；large eddy simulation；CFL number；OpenFOAM

TK11

A

1006-8740(2021)03-0321-13

10.11715/rskxjs.R2003003

2020-06-20．

國(guó)家自然科學(xué)基金資助項(xiàng)目(91641110)；國(guó)家重點(diǎn)研發(fā)計(jì)劃資助項(xiàng)目(2019YFB1704200).

姚?衛(wèi)（1983—??），男，博士，副研究員.

姚?衛(wèi)，weiyao@imech.ac.cn.

(責(zé)任編輯：武立有)