趙 瑞，閻 超，李新亮，蔣海軍

（1.北京航空航天大學(xué)國家計算流體力學(xué)實驗室，北京 100191；2.中國科學(xué)院力學(xué)所高溫氣體動力學(xué)重點實驗室，北京 100190；3.北京臨近空間飛行器系統(tǒng)工程研究所，北京 100076）

0 引言

熱力學(xué)第二定律，即熵增原理，被廣泛地應(yīng)用在計算傳熱學(xué)以及計算流體力學(xué)領(lǐng)域。熵是有效能量耗散程度的指標(biāo)，表征自然界中能量演化的不可逆過程。Naterer和Camberos［1］總結(jié)了熵概念的具體應(yīng)用：1）基本的熱力學(xué)準(zhǔn)則；2）解的唯一性，即附加熵條件來排除多余的廣義解；3）提高數(shù)值穩(wěn)定性；4）熵函數(shù)的凹性以及其他特有的性質(zhì)；5）時間步長限制；6）反設(shè)計方法；7）數(shù)值誤差以及收斂性判斷準(zhǔn)則；8）上風(fēng)格式精度影響；9）流體機(jī)械的優(yōu)化設(shè)計。另外，Roe等人采用網(wǎng)格單元的熵通量作為判斷指標(biāo)，進(jìn)行網(wǎng)格自適應(yīng)加 密［2－3］。最近，趙瑞等人［4］通過對超聲速壓縮拐角進(jìn)行直接數(shù)值模擬（DNS），發(fā)現(xiàn)無論在附著流動區(qū)域還是在拐角處的分離區(qū)，熵函數(shù)沿壁面法向都能夠保持很好的單調(diào)性，并利用該性質(zhì)重構(gòu)Baldwin－Lomax（BL）湍流模型的長度尺度，顯著提高了BL模型的魯棒性。

本文將繼續(xù)文獻(xiàn)［4］的研究工作，通過對比不同馬赫數(shù)下（從低速到高超聲速）沿壁面法向的熵增變化規(guī)律，提出比熵增概念。比熵增代表單位機(jī)械能的耗散程度，研究發(fā)現(xiàn)該物理量在數(shù)值上既能保持很好的單調(diào)性，同時具有馬赫數(shù)無關(guān)性，進(jìn)而增強(qiáng)了采用熵概念分辨邊界層的能力。文章結(jié)構(gòu)安排如下：第一節(jié)定義熵函數(shù)并通過對不同馬赫數(shù)下平板邊界層的DNS模擬提出比熵增概念，并構(gòu)造熵函數(shù)fs。第二節(jié)描述本文采用的數(shù)值方法。指出原始的BL 模型存在的缺陷以及如何利用熵函數(shù)進(jìn)行修正。第三節(jié)通過對后臺階流動（Ma＝0.128）以及高超聲速錐柱裙組合體（Ma＝7.05）進(jìn)行數(shù)值模擬，考核該修正方式的性能。最后一節(jié)對熵概念的應(yīng)用做綜合評價，并對該概念在CFD 領(lǐng)域的發(fā)展方向進(jìn)行展望。

1 比熵增及熵函數(shù)fs

在自然界和工程實際應(yīng)用中，最經(jīng)常發(fā)生的流動狀態(tài)是湍流。湍流是流體微團(tuán)的不規(guī)則運動，湍流運動產(chǎn)生的質(zhì)量、動量和能量的輸運將遠(yuǎn)遠(yuǎn)大于分子熱運動產(chǎn)生的宏觀輸運，同時湍流脈動也會導(dǎo)致額外的能量耗散，引起氣動加熱和摩阻的增加?？梢灶A(yù)見到，在壁面附近由于劇烈的湍流脈動以及壁面摩擦作用，當(dāng)?shù)氐臋C(jī)械能不可逆地轉(zhuǎn)化為內(nèi)能。而熵，作為有效能量耗散程度的指標(biāo)，可以在數(shù)值上很好地表征這一現(xiàn)象。

對于完全氣體，熵增與能量變化有如下關(guān)系［5］：

其中，s為單位質(zhì)量的熵，q為單位質(zhì)量的熱量，T代表當(dāng)?shù)販囟取?/p>

將能量守恒定律，即熱力學(xué)第一定律：

代入方程（1），整理得：

其中，u為單位質(zhì)量內(nèi)能，u＝cVT；p、ρ分別為當(dāng)?shù)貕毫?，密度；cV為質(zhì)量定容熱容；r為比熱比。由式（3）可以看出，s正比于p／ρr，本文中定義熵函數(shù)：

可以看出，熵s只與當(dāng)?shù)氐膲毫?、密度有關(guān)，是個點函數(shù)，和物系的具體變化過程無關(guān)，只表征當(dāng)?shù)啬芰康暮纳⒊潭取４送?，?dāng)外部勢流流經(jīng)物體時，由于壁面阻滯作用，速度逐步減小為0，可以預(yù)見到，在不同馬赫數(shù)下，所耗散掉的機(jī)械能不同，因此壁面的熵值也會大不一樣。為減小這種對不同來流情況的依賴性，本文提出比熵增的概念（用ˉs表示），用來表征單位機(jī)械能的耗散程度。

李新亮等人［6－8］采用DNS方法對不同馬赫數(shù)下的平板邊界層進(jìn)行模擬，實驗條件如表1所示。

表1 平板邊界層DNS計算參數(shù)Table 1 Simulation conditions

從中分別提取沿壁面法向的熵和比熵增，如圖1所示?？梢钥吹剑煌R赫數(shù)下，壁面熵值相差一個量級（圖1（a）），而壁面處比熵增基本在0.1～0.5之間，并且比熵增ˉs＞0 的范圍與邊界層重合（圖1（b））。

圖1 不同馬赫數(shù)下平板邊界層沿壁面法向熵和比熵增Fig.1 Entropy and entropy increment ratio profiles normal to the wall under different Mach numbers

此外，通常在網(wǎng)格生成時需要在壁面進(jìn)行適度加密，本文引入壁面距離來緩解高超聲速流動中激波引起的局部的熵增，構(gòu)造以下過渡函數(shù)：

其中dw為網(wǎng)格點到壁面最小距離，δs為沿壁面法向，fs＝0.99到最近壁面距離。

熵函數(shù)fs有以下性質(zhì)：當(dāng)fs＜1.0時，所指示的區(qū)域為邊界層，當(dāng)fs＝1.0時，所指示區(qū)域為外部勢流，如圖2所示。

圖2 不同馬赫數(shù)下平板邊界層沿壁面法向fs 分布Fig.2 Function fsprofiles normal to the wall with different Mach numbers

下文中，將以BL 模型的長度尺度修正為例，詳細(xì)描述熵函數(shù)fs的應(yīng)用。

2 計算方法

2.1 控制方程

數(shù)值計算采用有限體積法求解可壓縮雷諾平均方程，采用量熱完全氣體假設(shè)，粘性通量采用中心差分格式進(jìn)行離散，無粘通量選擇Roe的FDS格式，時間推進(jìn)采用LU－SGS隱式方法。

2.2 Baldwin－Lomax湍流模型

BL模型［9］對湍流邊界層的內(nèi)層與外層采用不同的混合長假設(shè)，其渦粘性如下所示：

這里yc是（μt）inner＝（μt）outer時y到壁面最小距離。

對于內(nèi)層，即y≤yc，有：

其中ρ為密度，Ω為渦量，K＝0.4為Karman常數(shù)，Van Driest衰減因子D如下：

其中τw為壁面摩擦力。對于外層，即y＞yc，有：

尾跡函數(shù)Fwake如下式所示：

其中Fmax＝max（yΩD），ymax是函數(shù)F＝y(tǒng)ΩD達(dá)到最大值Fmax的位置；Cwk＝1.0，Udiff是平均速度分布中最大值與最小值之差。

Klebanoff間歇函數(shù)Fkleb如下式所示：

其中Ckleb＝0.3，邊界層厚度δ＝y(tǒng)max／Ckleb。

大量研究表明，原始的BL模型只能準(zhǔn)確地預(yù)測附體流動以及弱分離流動，對于激波／邊界層干擾以及分離、再附等復(fù)雜流動卻無能為力。究其根本原因，一方面BL 模型中的經(jīng)驗常數(shù)Ccp＝1.6、Ckleb＝0.3是基于跨聲速的平衡流動得出，并不適用于其他馬赫數(shù)下的復(fù)雜流動；另一方面外層渦粘性的捕捉嚴(yán)重依賴于Fmax、ymax的確定，原始的BL 模型并沒有限制尋找Fmax的范圍，在復(fù)雜流動中會導(dǎo)致長度尺度ymax的不確定性［5，10］。

類似文獻(xiàn)［5］采用的修正方法，以BL 模型在超聲速壓縮拐角x／δ＝0處的長度尺度修正為例，整個搜尋過程分兩步完成：1）沿壁面法向由內(nèi)向外搜索直到fs＞0.95停止，儲存停止位置。2）函數(shù)F由停止位置反向搜索，確定Fmax與ymax，如圖3所示。另外，根據(jù)文獻(xiàn)［11］的建議，修正后的模型在低速流動中Ccp＝1.6、Ckleb＝0.65，在超聲速及高超聲速流動中采用Ccp＝2.08、Ckleb＝0.3。

圖3 搜索過程示意圖Fig.3 Searching steps

3 算例驗證

3.1 低速后臺階流動（Ma＝0.128）

后臺階流動［11］由于外形簡單，并且流動分離點確定，是考核湍流模型性能的標(biāo)準(zhǔn)算例。來流條件如表2所示，低速流動在臺階處發(fā)生大分離，并在下游逐漸恢復(fù)邊界層流動，流動示意圖如圖4（a）所示。壁面為無滑移的絕熱壁，網(wǎng)格在壁面以及臺階分離區(qū)進(jìn)行適度加密，第一層壁面網(wǎng)格高度保證y＋＜1（如圖4（b）所示）。

表2 后臺階流動來流條件（x／H＝－4，H：臺階高度）Table 2 Inflow conditions for backward－facing step（at x／H＝－4，H：the step height）

低速流體流經(jīng)臺階時發(fā)生分離再附，原始的BL模型（BL－origin）在分離區(qū)處由于沒有限制長度尺度的搜索范圍，會導(dǎo)致在該處渦粘性預(yù)測過小，出現(xiàn)明顯的渦粘性“斷層”，如圖6（a）所示。修正后的BL 模型（BL－entropy）采用的熵函數(shù)fs能夠合理的指示出邊界層的范圍（圖5），進(jìn)而限制長度尺度的搜索范圍，獲得較為合理的渦粘性分布（圖6（b））。但由于BL模型的平衡態(tài)假設(shè)，即渦粘性只由當(dāng)?shù)氐钠骄鲃訁?shù)的代數(shù)關(guān)系式所決定，因此對于后臺階這種大分離的復(fù)雜流動，對湍流場的預(yù)測能力要明顯低于SA 一方程湍流模型，如圖6（c）所示，SA 模型預(yù)測的渦粘性分布要比BL模型更為連續(xù)，而且在數(shù)值上明顯大于BL模型。

圖4 后臺階流場結(jié)構(gòu)及網(wǎng)格示意圖Fig.4 The flow structures and grid construction for the backward－facing step

圖5 后臺階流動中熵函數(shù)fs 分布云圖（只顯示fs＜0.95）Fig.5 Function fsdistributions in the flowfield（only depict fs＜0.95for convenience）

圖6 三種湍流模型預(yù)測的渦粘性分布云圖Fig.6 The distributions of eddy－viscosity around the step

后臺階下壁面壓力和摩擦力系數(shù)分布如圖7所示，可以看到，BL－origin預(yù)測的分離區(qū)明顯過大，修正后的BL－entropy和SA 模型預(yù)測的結(jié)果與實驗值較為吻合，但由于BL－entropy計算的渦粘性偏低，分離區(qū)預(yù)測依舊偏大。

圖7 壁面壓力系數(shù)及摩擦力系數(shù)分布Fig.7 Pressure and skin－friction distributions along the step－side wall

3.2 高超聲速錐柱裙組合體流動（Ma＝7.05）

Kussoy和Horstman［12］對高超聲速錐柱裙組合體流動進(jìn)行了實驗研究，具體實驗條件如表3所示。高超聲速氣流流至拐角前段時，在逆壓梯度的作用下，邊界層變厚，從而使主流提前產(chǎn)生一系列的壓縮波，并在離拐角一定距離處匯聚成一道斜激波，激波后邊界層往往不能承受強(qiáng)逆壓梯度而發(fā)生分離，形成一個“凸包”，高超聲速氣流流經(jīng)該“凸包”后緊接著在斜面的壓縮作用下，形成新的壓縮波區(qū)，流動示意圖如圖8（a）所示。圖8（b）為該算例的網(wǎng)格示意圖，在壁面以及分離區(qū)附近進(jìn)行加密，保證壁面第一層網(wǎng)格高度y＋＜1。由于拐角處強(qiáng)烈的激波／邊界層干擾，BL－origin在該算例中計算發(fā)散。

表3 錐柱裙組合體流動入口條件（x＝－6cm）Table 3 The inflow conditions for cylinder with conical flare（at x＝－6cm）

圖8 35°錐柱裙組合體流動結(jié)構(gòu)及網(wǎng)格示意圖Fig.8 The flow structures and grid construction for cylinder with conical flare at 35°

圖9為該算例中熵函數(shù)fs的分布云圖，可以看到fs準(zhǔn)確地捕捉到邊界層在再附區(qū)域的發(fā)展歷程。然而對于高超聲速復(fù)雜流動，BL－entropy和SA 模型預(yù)測能力都明顯下降。如圖10所示，在分離區(qū)以及再附區(qū)域，雖然fs將BL－entropy預(yù)測的渦粘性很好地限制在壁面附近，但沿流向出現(xiàn)明顯間斷，而SA模型在再附點附近預(yù)測的渦粘性明顯不足。

圖9 熵函數(shù)fs 分布云圖（只顯示fs＜0.95）Fig.9 Function fsdistributions in the flowfield（only depict fs＜0.95for convenience）

圖10 渦粘性分布云圖Fig.10 The distributions of eddy－viscosity

計算所得的壁面壓力與熱流分布如圖11所示，BL－entropy能夠很好的地預(yù)測處分離位置，但在再附點處壓力出現(xiàn)明顯的過沖，而SA 模型預(yù)測的分離區(qū)過小，壓力的恢復(fù)速度也要更快。在高超聲速流動中，即使在附著流動中（x≤0cm），SA 模型預(yù)測的熱流值也明顯過大見圖11（b），整個再附區(qū)域的熱流值也要遠(yuǎn)遠(yuǎn)大于實驗值。與壓力分布類似，BL－entropy在再附點處出現(xiàn)熱流過沖，但拐角區(qū)域的熱流預(yù)測與實驗值十分吻合。說明本文提出的BL－entropy對于高超聲速復(fù)雜流動也能夠進(jìn)行合理的預(yù)測。

圖11 壁面壓力及熱流分布Fig.11 Pressure and heat flux distributions along the wall

4 結(jié)論

比熵增表征單位機(jī)械能的耗散，由此構(gòu)造的熵函數(shù)fs同時具備單調(diào)性和馬赫數(shù)無關(guān)性，因此可以從數(shù)值上保證對邊界層的分辨能力。采用熵函數(shù)重構(gòu)BL湍流模型的長度尺度，并通過對各種工況下的（從低速流動到高超聲速流動）分離再附流動、激波／邊界層干擾流動進(jìn)行數(shù)值模擬，充分驗證了該修正方法的可靠性。但同時介于BL模型本身的平衡態(tài)假設(shè)，渦粘性的預(yù)測不如SA 模型的光滑，在低速流動中預(yù)測性能也要低于SA 湍流模型。

諸多湍流模型、轉(zhuǎn)捩模型以及可壓縮修正方法都是借助邊界層流動進(jìn)行模型參數(shù)的標(biāo)定，因此通過對邊界層范圍的合理預(yù)測來約束以及修正這些參數(shù)，將會提升這些模型對不同流動狀況的預(yù)測能力。熵概念從能量耗散的角度上預(yù)測邊界層的范圍，更具有普適性。因此如何通過合理的構(gòu)造熵函數(shù)改進(jìn)現(xiàn)有模型的性能，也將是我們下一步工作。

［1］NATERER G F，CAMBEROS J A.The role of entropy and the second law in computational thermofluids［R］.AIAA 2001－2758.

［2］FIDKOWSKI K J，ROE P L.Entropy－based mesh refinement，I：the entropy adjoint approach［R］.AIAA 2009－3790.

［3］ZAIDE D W，ROE P L.Entropy－based mesh refinement，II：a new approach to mesh movement［R］.AIAA 2009－3791.

［4］趙瑞，閻超，于劍.一種基于熵限制的Baldwin－Lomax湍流模型［J］.北京航空航天大學(xué)學(xué)報，2012，38（2）：175－179.

［5］GREINER W，NEISE L，ST?CKER H.Thermodynamcis and statistical mechanics［M］.Springer－Verlag ：New York Inc.1995.

［6］ZHOU Y，LI X L，F(xiàn)U D X，MA Y W.Coherent structures in the transition of a flat－plate boundary layer atMa＝0.7［J］.ChinesePhysicsLetter，2007，24（1）：147－150.

［7］LI X L，F(xiàn)U D X，MA Y W，et al.Acoustic calculation for super－sonic turbulent boundary layer flow［J］.ChinesePhysicsLetter，2009，26（9）：181－184.

［8］LI X L，F(xiàn)U D X，MA Y W.Direct numerical simulation of a spa－tially evolving supersonic turbulent boudary layer atMa＝6［J］.ChinesePhysicsLetter，2006，23（6）：1519－1522.

［9］BALDWIN B S，LOMAX H.Thin layer approximation and alge－braic model for separated turbulent flows［R］.AIAA Paper 78－0257，1978.

［10］VISBAL M，KNIGHT D.The Baldwin－Lomax turbulence model for two－dimensional shock－wave／boundarylayer interactions［J］.AIAAJournal，1984，22（7）：921－928.

［10］DRIVER D M.Features of a reattaching turbulent shear layer in divergent channel flow［J］.AIAAJournal，1985，23（1）：163－171.

［11］SETTLES G，DODSON L.Hypersonic shock／boundarylayer interaction database［R］.AIAA Paper 91－1763，1991.