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

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

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

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

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

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

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

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

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

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

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

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

1.1?控制方程

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

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

1.4?基于局部時間步的高效多重時間法

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

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

2?結(jié)果與討論

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

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

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

圖2?火焰中心剖面局步時間步分布

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

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

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

圖5?中心軸線上時均值比較

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

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

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

圖6?時間步分布

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

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

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

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

3?結(jié)?語

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

一般而言，局部時間步法針對穩(wěn)態(tài)流動或統(tǒng)計上呈現(xiàn)穩(wěn)態(tài)的準(zhǔn)穩(wěn)態(tài)流動問題，其得到的結(jié)果是比較可靠的，但對于強(qiáng)瞬態(tài)的流場演化問題，其方法的適用性則需要進(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］ 杜永樂，閻?超. DSMC方法中的自適應(yīng)當(dāng)?shù)貢r間步長法[J]. 北京航空航天大學(xué)學(xué)報，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］ 姜婷婷，錢?坤，陳偉芳. 高超聲速混合流動的粒子模擬耦合算法[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］ 譚志軍，黃云清. 對一維守恒律的一種局部時間步長自適應(yīng)網(wǎng)格方法[J]. 湘潭大學(xué)自然科學(xué)學(xué)報，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］ 吳?迪，蔚喜軍，徐?云. 局部時間步長間斷有限元方法求解三維歐拉方程[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］ 胡?鵬，韓健健，雷云龍. 基于局部分級時間步長方法的水沙耦合數(shù)學(xué)模擬[J]. 浙江大學(xué)學(xué)報(工學(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)過程調(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燃煤鍋爐污泥摻燒試驗與數(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．

國家自然科學(xué)基金資助項目(91641110)；國家重點研發(fā)計劃資助項目(2019YFB1704200).

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

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

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