杜磊 孫波 代春良

摘要：基于開(kāi)源計(jì)算流體力學(xué)軟件OpenFOAM中的密度基求解器rhoCentralFoam植入典型混合格式AUSM+。求解器通過(guò)兩組算例進(jìn)行驗(yàn)證：二維圓柱繞流數(shù)值模擬和二維高超聲速進(jìn)氣道數(shù)值模擬。與圓柱繞流模型相似的高超聲速進(jìn)氣道鈍化前緣，是高超聲速進(jìn)氣道典型流動(dòng)特征。數(shù)值計(jì)算結(jié)果分別與商業(yè)軟件CFD++計(jì)算結(jié)果和試驗(yàn)結(jié)果進(jìn)行對(duì)比，結(jié)果表明，基于OpenFOAM植入的AUSM+格式可以精確捕捉脫體激波、準(zhǔn)確構(gòu)建流場(chǎng)結(jié)構(gòu)，并得到可靠的流場(chǎng)參數(shù)，與商業(yè)軟件CFD++計(jì)算結(jié)果和試驗(yàn)結(jié)果均具有良好的一致性，滿足高超聲速進(jìn)氣道流場(chǎng)計(jì)算需求。

關(guān)鍵詞：高超聲速流動(dòng)；OpenFOAM；AUSM+；通量分裂；數(shù)值模擬

中圖分類號(hào)：V211.3文獻(xiàn)標(biāo)識(shí)碼：ADOI：10.19452/j.issn1007-5453.2020.11.015

對(duì)于高超聲速流動(dòng)數(shù)值計(jì)算而言，一般認(rèn)為求解器的數(shù)值格式應(yīng)具備激波捕捉、高黏性分辨率以及避免出現(xiàn)“粉刺現(xiàn)象”等特性[1]。通量差分分裂（FDS）格式，對(duì)接觸間斷和邊界層都有很高的分辨率，但其魯棒性較差。矢通量分裂（FVS）格式在激波捕捉方面有較好的魯棒性，但在接觸間斷和邊界層內(nèi)有較大的數(shù)值耗散[2]。AUSM類格式兼有FDS格式的高分辨率特性和FVS格式的強(qiáng)魯棒特性。

高超聲速進(jìn)氣道作為吸氣式推進(jìn)系統(tǒng)的重要組成部分，一般位于高超聲速飛行器流場(chǎng)上游區(qū)域，其流動(dòng)特性對(duì)下游流場(chǎng)具有重要影響。高超聲速進(jìn)氣道的鈍化前緣，與圓柱繞流形成的脫體激波流場(chǎng)結(jié)構(gòu)具有一定相似性。鈍化前緣在引起進(jìn)氣道流場(chǎng)結(jié)構(gòu)變化的同時(shí)進(jìn)一步影響進(jìn)氣道的工作特性，準(zhǔn)確預(yù)測(cè)前緣鈍化后的流場(chǎng)特性，對(duì)進(jìn)氣道的修正設(shè)計(jì)提供指導(dǎo)和依據(jù)[3]。

OpenFOAM中的求解器大多數(shù)為壓力基求解器，主要用于求解不可壓?jiǎn)栴}；為了準(zhǔn)確預(yù)測(cè)高速流動(dòng)的流場(chǎng)結(jié)構(gòu)及參數(shù)，OpenFOAM發(fā)展了密度基求解器rhoCentralFoam，并應(yīng)用了Kurganov和Tadmor[4]的中心迎風(fēng)格式。Borm[5]等開(kāi)發(fā)了用于葉輪機(jī)械數(shù)值模擬的densityBasedTurbo求解器，對(duì)流通量的計(jì)算使用了以Roe格式為代表的Godunov類格式。Shen[6]等對(duì)densityBasedTurbo求解器中的數(shù)值格式進(jìn)行了較為詳細(xì)的驗(yàn)證，并在此基礎(chǔ)上利用時(shí)間導(dǎo)數(shù)預(yù)處理方法將密度基求解器的應(yīng)用推廣至低速流。本文擬基于開(kāi)源軟件OpenFOAM植入典型的混合格式AUSM+格式，基于rhoCentralFoam求解器采用算子分裂的方法進(jìn)行解算，以滿足高超聲速進(jìn)氣道流場(chǎng)的計(jì)算需求。通過(guò)二維圓柱繞流的數(shù)值模擬驗(yàn)證AUSM+格式計(jì)算所得的主要流場(chǎng)參數(shù)，如壓力分布、激波距離等，驗(yàn)證AUSM+格式激波捕捉能力。選取二維高超聲速進(jìn)氣道，對(duì)其流場(chǎng)結(jié)構(gòu)及參數(shù)進(jìn)行計(jì)算，驗(yàn)證所植入的數(shù)值格式可信度。

1數(shù)值方法及其實(shí)現(xiàn)

1.1控制方程

1.3代碼實(shí)現(xiàn)

開(kāi)源求解器gFoam含有Roe等Godunov類格式，此處參考gFoam求解器實(shí)現(xiàn)AUSMPLUSFlux（…）通量計(jì)算函數(shù)。求解器程序結(jié)構(gòu)框圖如圖1所示，基于rhoCentralFoam求解器，通過(guò)在其代碼主循環(huán)中調(diào)用通量計(jì)算函數(shù)進(jìn)行AUSM+格式在OpenFOAM框架下的實(shí)現(xiàn)。通量計(jì)算函數(shù)的實(shí)現(xiàn)利用了OpenFOAM中提供的類，如const scalar alpha = 3/16語(yǔ)句，其中scalar為OpenFOAM中提供的C++張量類，表示零階張量（標(biāo)量）。標(biāo)量alpha對(duì)應(yīng)著分裂壓強(qiáng)定義式中的α；類似地，相應(yīng)的代碼實(shí)現(xiàn)中phi、phiUp、phiEp即為AUSM+理論方法中的通量項(xiàng)。求解器通過(guò)調(diào)用AUSMPLUSFlux（…）函數(shù)分別給出phi、phiUp、phiEp在求解域內(nèi)部和邊界上的對(duì)流通量值。

2計(jì)算結(jié)果與分析

本節(jié)采用與商業(yè)軟件CFD++和試驗(yàn)結(jié)果分別進(jìn)行對(duì)比的方法來(lái)驗(yàn)證基于OpenFOAM植入的數(shù)值格式可信度。驗(yàn)證算例分別為二維圓柱繞流和二維高超聲速進(jìn)氣道，出口均采用零梯度邊界條件且壁面為無(wú)滑移等溫壁面，來(lái)流條件見(jiàn)算例。

2.1圓柱算例

選取Holden[11]開(kāi)展的低密度層流圓柱繞流試驗(yàn)?zāi)Ｐ?，試?yàn)過(guò)程中整個(gè)流場(chǎng)為低焓值層流狀態(tài)，且未發(fā)生化學(xué)反應(yīng)。計(jì)算所采用的幾何模型為半徑R=0.0381m的圓柱，簡(jiǎn)化求解域在二維對(duì)稱條件下進(jìn)行數(shù)值模擬，且求解域在來(lái)流方向與垂直來(lái)流方向尺寸分別為0.75R、2.5R。來(lái)流條件馬赫數(shù)為16.01，速度為2111.045m/s，溫度為43.214K，壓強(qiáng)為21.974Pa，壁面溫度為300.333K。

2.1.1與試驗(yàn)結(jié)果對(duì)比及分析

本節(jié)通過(guò)三種近壁面首層網(wǎng)格高度不同的網(wǎng)格進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證，Grid 1、Grid 2以及Grid 3的近壁面首層網(wǎng)格高度分別為1μm、20μm、50μm，且三種網(wǎng)格之間均有一定繼承性。圖2為圓柱壁面沿程壓力分布，三種網(wǎng)格計(jì)算結(jié)果一致性較高。Grid 3計(jì)算結(jié)果與Grid 1計(jì)算結(jié)果重合度較高，但是分布曲線振蕩較為明顯。Grid 2計(jì)算結(jié)果分布曲線雖然沒(méi)有同Grid 1一樣產(chǎn)生振蕩，并且與Grid 1計(jì)算結(jié)果也較為一致，但是其計(jì)算結(jié)果在駐點(diǎn)附近產(chǎn)生了較大的振蕩?？紤]到計(jì)算準(zhǔn)確度，本文選用Grid 1（近壁面首層網(wǎng)格高度1μm）的計(jì)算結(jié)果進(jìn)行下文分析與驗(yàn)證。

對(duì)高超聲速二維圓柱繞流試驗(yàn)所采集的圓柱壁面沿程壓力系數(shù)分布試驗(yàn)數(shù)據(jù)與OpenFOAM_AUSM+格式數(shù)值計(jì)算結(jié)果相對(duì)比，如圖3所示[12]。結(jié)果表明，AUSM+格式計(jì)算所得圓柱壁面沿程壓強(qiáng)分布與試驗(yàn)基本吻合，在駐點(diǎn)附近有一定差異且曲線變得不光滑，在壁面沿程下游重合度較好且沿程曲線差異呈逐漸減小趨勢(shì)，計(jì)算結(jié)果在駐點(diǎn)附近相比于壁面沿程下游更為符合試驗(yàn)數(shù)據(jù)，其在壁面沿程下游均高于試驗(yàn)數(shù)據(jù)，所采用的數(shù)值格式能精確計(jì)算壓力等流場(chǎng)參數(shù)。

2.1.2與CFD++結(jié)果對(duì)比及分析

OpenFOAM與商業(yè)軟件CFD++的計(jì)算結(jié)果對(duì)比如圖4所示，分別為流場(chǎng)壓力云圖（見(jiàn)圖4（a））和溫度云圖（見(jiàn)圖4（b）），云圖上下部分分別為CFD++和OpenFOAM計(jì)算結(jié)果，其中CFD++為默認(rèn)的HLLC格式，OpenFOAM則為典型的混合格式AUSM+格式。

計(jì)算結(jié)果表明，OpenFOAM_AUSM+格式計(jì)算所得激波距離與CFD++_HLLC格式計(jì)算所得激波距離基本一致，各云圖中對(duì)稱位置處等值線均能一一對(duì)應(yīng)。微小區(qū)別出現(xiàn)在OpenFOAM云圖一側(cè)靠近激波位置處，其壓力云圖、溫度云圖中等值線均不夠光滑，但等值線分布均能與CFD++計(jì)算結(jié)果云圖相匹配。

綜合數(shù)值格式計(jì)算結(jié)果的對(duì)比表現(xiàn)，OpenFOAM能夠準(zhǔn)確預(yù)測(cè)出鈍化前緣處的脫體激波以及波后流場(chǎng)，是準(zhǔn)確預(yù)測(cè)高超聲速進(jìn)氣道內(nèi)部流動(dòng)的前提。

2.2進(jìn)氣道算例

計(jì)算對(duì)象采用德國(guó)航空航天中心（DLR）在高超聲速風(fēng)洞H2K中進(jìn)行試驗(yàn)的二維高超聲速進(jìn)氣道模型，如圖5所示[13]。來(lái)流條件馬赫數(shù)為7，進(jìn)氣道處于亞額定狀態(tài)，總溫為500K，總壓為7×105Pa，密度為0.0123kg/m3，分別對(duì)應(yīng)來(lái)流靜溫為46K，靜壓為170Pa。壁面邊界條件采用無(wú)滑移等溫邊界，壁面溫度取300K。在保證壁面網(wǎng)格正交性的前提下，對(duì)近壁面網(wǎng)格進(jìn)行局部加密，二維網(wǎng)格總數(shù)達(dá)到10萬(wàn)量級(jí)。采用k -ωSST湍流模型，首層網(wǎng)格高度滿足增強(qiáng)壁面函數(shù)所需y+[14-17]。

2.2.1與試驗(yàn)結(jié)果對(duì)比及分析

高超聲速進(jìn)氣道入口分離區(qū)較大，激波系復(fù)雜，隔離段內(nèi)部則表現(xiàn)為單激波反射，并且激波系沿下游方向逐漸減弱，如圖6所示為進(jìn)氣道數(shù)值紋影圖。唇罩激波（a）與自前體發(fā)展而來(lái)的邊界層干擾形成了分離區(qū)（b），氣流在進(jìn)氣道肩部先膨脹后壓縮，此處壓縮波（c）即為分離區(qū)誘導(dǎo)而來(lái)的分離激波，并與唇罩激波（a）相交，形成異側(cè)激波相交結(jié)構(gòu)。分離區(qū)（b）的再附激波（d）在上壁面處與邊界層干擾，形成了與分離區(qū)（b）結(jié)構(gòu)相似的分離區(qū)（f），此分離區(qū)相對(duì)于分離區(qū)（b）較小，原因則是激波強(qiáng)度的減弱。同樣的，隨著沿流動(dòng)方向激波系強(qiáng)度的減弱，再附激波（g）、反射激波（h、i）并沒(méi)有導(dǎo)致壁面處出現(xiàn)分離區(qū)。

進(jìn)氣道內(nèi)部上、下壁面的壓力系數(shù)分布如圖7所示。下壁面的分離激波（c）、再附激波（d）、反射激波（h）分別導(dǎo)致在上壁面x=0.41m、x=0.475、x=0.56m位置附近出現(xiàn)波后壁面壓力峰值，并且峰值大小逐步減小。同樣的，進(jìn)氣道下壁面的壁面壓力峰值出現(xiàn)原因與上壁面相似，唇口激波（a）、再附激波（g）分別導(dǎo)致在下壁面x=0.438m、x=0.515位置附近出現(xiàn)波后壁面壓力峰值，變化趨勢(shì)與上壁面相似。隔離段下游的壓力分布曲線相對(duì)上游較為平滑，單激波反射致使壁面出現(xiàn)了典型的隔離段壓力峰值交變現(xiàn)象。進(jìn)氣道上、下壁面的壓力分布曲線變化趨勢(shì)與試驗(yàn)數(shù)據(jù)較為符合，在隔離段中壓力分布曲線與試驗(yàn)數(shù)據(jù)有一定的差異，體現(xiàn)在壓力峰值的大小與峰值點(diǎn)出現(xiàn)位置不同?？紤]到數(shù)值計(jì)算部分邊界條件的選取可能與試驗(yàn)條件有所區(qū)別，如壁面溫度的選取，且誤差范圍較小，所以計(jì)算結(jié)果較為符合實(shí)際流動(dòng)[18]。

2.2.2與CFD++結(jié)果對(duì)比及分析

圖8為進(jìn)氣道馬赫數(shù)分布云圖，上、下部分分別為OpenFOAM計(jì)算結(jié)果與CFD++計(jì)算結(jié)果，圖中使用紅色線條標(biāo)注出分離區(qū)流線，二者計(jì)算所得分離區(qū)（b）大小基本一致，并且分離點(diǎn)位置基本相同，均稍靠后于唇罩前緣。云圖中的激波系位置一一對(duì)應(yīng)，自前體發(fā)展而來(lái)的邊界層厚度也基本一致，且OpenFOAM計(jì)算結(jié)果在隔離段中的反射激波較為明顯，形成的分離區(qū)（f）較為明顯。二者計(jì)算結(jié)果差別較小，僅在進(jìn)氣道隔離段出口處等值線形狀出現(xiàn)微小差別，除此之外等值線形狀與位置基本一致。圖9為進(jìn)氣道溫度分布云圖，兩者對(duì)比結(jié)果與馬赫數(shù)云圖類似，計(jì)算結(jié)果一致性較高。

為了進(jìn)一步對(duì)比OpenFOAM計(jì)算結(jié)果與CFD++計(jì)算結(jié)果在隔離段出口的微小差別，高超聲速進(jìn)氣道隔離段出口馬赫數(shù)分布如圖10所示。OpenFOAM計(jì)算結(jié)果與CFD++計(jì)算結(jié)果變化趨勢(shì)與試驗(yàn)測(cè)量值的變化趨勢(shì)一致，數(shù)值模擬結(jié)果一致性較高，二者預(yù)測(cè)出的隔離段出口峰值馬赫數(shù)相比于試驗(yàn)測(cè)量值均偏大，可能是由于試驗(yàn)值并非連續(xù)測(cè)量的結(jié)果，試驗(yàn)過(guò)程中未能準(zhǔn)確測(cè)得峰值大小及其位置所在。

綜合流場(chǎng)結(jié)構(gòu)、流場(chǎng)參數(shù)的對(duì)比結(jié)果，OpenFOAM_ AUSM+格式計(jì)算結(jié)果與試驗(yàn)值和商業(yè)軟件CFD++_HLLC格式計(jì)算結(jié)果一致性較高，能夠準(zhǔn)確預(yù)測(cè)流場(chǎng)，計(jì)算結(jié)果滿足高超聲速進(jìn)氣道流場(chǎng)預(yù)測(cè)需求。

3結(jié)論

本文基于開(kāi)源計(jì)算流體力學(xué)軟件OpenFOAM植入典型混合格式AUSM+格式，對(duì)二維圓柱繞流、二維高超聲速進(jìn)氣道進(jìn)行數(shù)值模擬，并與CFD++計(jì)算結(jié)果和試驗(yàn)結(jié)果進(jìn)行對(duì)比，得出以下結(jié)論：

（1）高超聲速可壓縮流動(dòng)控制方程離散關(guān)鍵在于對(duì)流項(xiàng)數(shù)值格式的選擇，基于OpenFOAM植入的AUSM+格式可以準(zhǔn)確預(yù)測(cè)流場(chǎng)結(jié)構(gòu)、流場(chǎng)參數(shù)，本文所選取的二維圓柱繞流、二維高超聲速進(jìn)氣道驗(yàn)證算例計(jì)算結(jié)果與商業(yè)軟件CFD++和試驗(yàn)數(shù)據(jù)均具有良好的一致性，滿足高超聲速進(jìn)氣道流場(chǎng)預(yù)測(cè)需求。

（2）本文選用的計(jì)算對(duì)象均為二維試驗(yàn)?zāi)Ｐ退憷瑢?duì)三維模型的復(fù)雜流場(chǎng)結(jié)構(gòu)準(zhǔn)確計(jì)算需要進(jìn)一步驗(yàn)證。對(duì)于AUSM+的改進(jìn)格式，如界面處聲速計(jì)算方法，仍有待進(jìn)一步研究。

（3）OpenFOAM中已植入的模塊均可與求解器進(jìn)行耦合調(diào)用，如湍流模型、化學(xué)反應(yīng)模型等。高超聲速流動(dòng)帶來(lái)的高溫氣體效應(yīng)可基于此求解器耦合相應(yīng)模塊后進(jìn)行深入研究。

參考文獻(xiàn)

[1]錢戰(zhàn)森，楊希明，李椿萱.高超聲速鈍頭體繞流氣動(dòng)熱計(jì)算的數(shù)值格式研究中存在的問(wèn)題[C]//探索創(chuàng)新交流——中國(guó)航空學(xué)會(huì)青年科技論壇， 2014. Qian Zhansen， Yang Ximing， Lee Chunhian. Problems existing in the numerical scheme research of hypersonic blunt nose body flow [C]// Exploration and Innovation Exchange-Youth scienceandTechnologyForumofChina Aeronautical Association，2014.（in Chinese）

[2]Shen C，F(xiàn)engxian S，Xinlin X. Analysis on capabilities of density-basedsolverswithinOpenFOAMtodistinguish aerothermal variables in difusion boundary layer[J]. Chinese Journal ofAeronautics，2013（6）：21-30.

[3]楊波，白菡塵，范孝華，等.前緣鈍度對(duì)馬赫數(shù)6平面壓縮進(jìn)氣道流場(chǎng)的影響分析[J].推進(jìn)技術(shù)，2018，39（3）：501-509. Yang Bo， Bai Hanshen， Fan Xiaohua， et al. Influence of leading edge bluntness on flow field of Mach 6 plane compression inlet [J]. Propulsion Technology， 2018，39 （3）： 501-509.（in Chinese）

[4]Kurganov A，Tadmor E. New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[J]. Journal of Computational Physics，2000，160（1）：241-282.

[5]Borm O，Jemcov A，Kau H P. Density based Navier Stokes solver for transonicflows[C]// 6th OpenFOAM Workshop，Penn State University，USA，2011.

[6]Shen C，Sun F，Xia X. Implementation of density-based solver for all speeds in the framework of OpenFOAM[J]. Computer Physics Communications，2014，185（10）：2730-2741.

[7]Greenshields C J，Weller H G，Gasparini L，et al. Implementation of semi-discrete，non-staggered central schemes in a colocated，polyhedral，finite volume framework，for high-speed viscous flows[J]. International Journal for Numerical Methods in Fluids，2010，63（1）：1-21.

[8]Liou M S，Jr Steffen C J. A new flux splitting scheme[J]. Journal of Computational Physics，1993，107（1）：23-39.

[9]Liou M S. A sequel to AUSM：AUSM～ + [J]. Journal of Computational Physics，1996，129（2）：364-382.

[10]Liou M S. Ten years in the making AUSM-family[C]// 15th AIAAComputational Fluid Dynamics Conference，2001.

[11]Holden M，Kolly J，Martin S. Shock/shock interaction heating in laminar and low-density hypersonic flows[C]// Thermo physics Conference，1996.

[12]Lee J，Rho O. Accuracy of AUSM+ scheme in hypersonic blunt body flow calculations[C]// 8th AIAA International Space Planes and Hypersonic Systems and Technologies Conference，1998.

[13]Hohn O，Gülhan A. Experimental investigation on the influence of yaw angle on the inlet performance at Mach 7[C]// 48th AIAAAerospace Sciences Meeting，2010.

[14]Menter F R，Kuntz M，Langtry R. Ten years of industrial experience with the SST turbulence model[J]. Heat and Mass Transfer，2003，4（1）：625-632.

[15]Menter F R. Zonal two equation k-ωmodels for aerodynamic flows[R]. 1993.

[16]Menter F R. Two-equation eddy-viscosity turbulence models for engineering applications[J]. AIAA Journal，1994，32（8）：1598-1605.

[17]Chan W Y K，Jacobs P A，Mee D J. Suitability of the k-ωturbulencemodelforscramjetflowfieldsimulations[J]. International Journal for Numerical Methods in Fluids，2012，70（4）：493-514.

[18]Haberle J，Gülhan A. Investigat-ion of two-dimensional scramjet inlet flowfield at Mach 7[J]. Journal of Propulsion and Power，2008，24（3）：446-459.（責(zé)任編輯陳東曉）

作者簡(jiǎn)介

杜磊（1997-）男，碩士研究生。主要研究方向：高超聲速氣體動(dòng)力學(xué)。

E-mail：dul@njust.edu.cn

孫波（1980-）男，博士，副教授。主要研究方向：高超聲速氣體動(dòng)力學(xué)。

Tel：13912999104E-mail：hypersun@126.com

代春良（1995-）男，博士研究生。主要研究方向：高超聲速氣體動(dòng)力學(xué)。

E-mail：DCL3839@126.com

Application of AUSM+ Scheme in Numerical Simulation of Hypersonic Inlet

Du Lei，Sun Bo*，Dai Chunliang

Nanjing University of Science and Technology，Nanjing 210094，China

Abstract： Based on the density-base solver rhoCentralFoam of open source computational fluid dynamics（CFD） software OpenFOAM， a typical hybrid scheme AUSM+ is implemented. The solver is verified by two sets of numerical examples： numerical simulation of the flow around the two-dimensional cylinder； numerical simulation of the twodimensional hypersonic inlet. The blunted leading edge of the hypersonic inlet， which is similar to the cylindrical flow model， is a typical flow characteristic of the hypersonic inlet. The numerical calculation results are compared with the commercial software CFD++ calculation results and experiment data respectively. The results show that the AUSM+ scheme based on OpenFOAM can capture the detached shock wave accurately， construct the flow field structure accurately and obtain the reliable flow field parameters. The results are in agreement with the commercial software CFD++ calculation results and experimental results， and satisfy the caculation requirements of hypersonic inlet flow field.

Key Words： hypersonic flow； OpenFOAM； AUSM+； flux splitting； numerical simulation