AUSM+格式在高超聲速進氣道數(shù)值模擬中的應(yīng)用

杜磊 孫波 代春良

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

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

中圖分類號：V211.3文獻標識碼：ADOI：10.19452/j.issn1007-5453.2020.11.015

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

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

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

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

1.1控制方程

1.3代碼實現(xiàn)

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

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

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

2.1圓柱算例

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

2.1.1與試驗結(jié)果對比及分析

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

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

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

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

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

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

2.2進氣道算例

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

2.2.1與試驗結(jié)果對比及分析

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

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

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

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

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

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

3結(jié)論

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

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

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

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

參考文獻

[1]錢戰(zhàn)森，楊希明，李椿萱.高超聲速鈍頭體繞流氣動熱計算的數(shù)值格式研究中存在的問題[C]//探索創(chuàng)新交流——中國航空學(xué)會青年科技論壇， 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]楊波，白菡塵，范孝華，等.前緣鈍度對馬赫數(shù)6平面壓縮進氣道流場的影響分析[J].推進技術(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.（責任編輯陳東曉）

作者簡介

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

E-mail：dul@njust.edu.cn

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

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

代春良（1995-）男，博士研究生。主要研究方向：高超聲速氣體動力學(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