王新光，陳 琦，萬 釗，張愛婧，陳堅強

(1. 空氣動力學國家重點實驗室，綿陽 621010；2. 中國空氣動力研究與發(fā)展中心計算所，綿陽 621000)

0 引 言

近年來，國家對于研發(fā)具有自主知識產權的自主品牌計算流體力學軟件，給予了前所未有的重視。眾所周知，湍流流動廣泛存在于航空、航天、航海以及人們的日常生活中，湍流模擬方法是計算流體力學界長期以來面臨的研究難題和研究熱點。商業(yè)軟件中通常使用RANS方法模擬湍流，常常需要在黏性底層中布置一定數(shù)目的網(wǎng)格點(通常y+≈1)，從而使得壁面附近的網(wǎng)格十分細密，不僅會極大地增加迭代收斂的計算步數(shù)，還會帶來較為嚴重的數(shù)值剛性問題，從而導致迭代計算的穩(wěn)定性下降。在復雜外形的數(shù)值計算中，網(wǎng)格量驟增，導致90%的計算時間都消耗在壁面附近1%計算域內。因此，目前主流商業(yè)軟件在工程湍流模擬應用中默認采用壁面函數(shù)方法對湍流邊界層進行模擬[1-2]，壁面函數(shù)方法可以大幅放寬壁面第一層網(wǎng)格的尺度(通常30

基于對數(shù)定律的標準壁面函數(shù)已經(jīng)被大量不可壓縮數(shù)值實驗確認，其中包括DNS模擬零壓力梯度邊界層和槽道流等結果。隨著近年DNS模擬能力的增強，涌現(xiàn)了大量可壓縮湍流邊界層DNS結果,其中文獻[4]使用DNS模擬了可壓縮低雷諾數(shù)槽道湍流，當馬赫數(shù)較小時標準壁面函數(shù)和DNS結果較為吻合，但是隨著馬赫數(shù)的增加基于不可壓流動的標準壁面函數(shù)和DNS之間的偏差逐漸增大。對于Ma2.48的槽道湍流，在黏性底層邊緣附近y+=10，相較于壁面參數(shù)，密度減小了50%，溫度增高了60%，表明低雷諾數(shù)超聲速槽道湍流存在強可壓縮性。

盡管壁面函數(shù)的研究工作已經(jīng)持續(xù)開展了50多年，但其中大多數(shù)研究成果的應用工作并不樂觀，如壁面函數(shù)與用于全局流動問題模擬的湍流模型不相容的問題[5-6]，容易引起數(shù)值模擬結果對壁面附近第一層網(wǎng)格節(jié)點位置的強烈依賴性，導致分離流動的數(shù)值模擬結果精度很差，因而提出了網(wǎng)格和流動自適應的壁面函數(shù)方法，保證在壁面函數(shù)有效的情況下對近壁流動物理特性進行恰當?shù)姆直?，特別是針對非平衡效應顯著的流動區(qū)域、逆壓梯度導致的流動分離和再附點附近區(qū)域，以獲得比較準確的氣動力/熱特性。此外，文獻[7]也針對壁面函數(shù)在RANS模擬中的應用問題開展了相關研究工作。

近年來可壓縮壁面函數(shù)的代表性工作，文獻[8]使用Crocco-Busemann方程，忽略壓力梯度的影響，將速度和溫度函數(shù)耦合起來，發(fā)展了適用于可壓縮湍流邊界層的壁面函數(shù)。國內壁面函數(shù)研究起步較晚，其中文獻[9-10]將文獻[8]中考慮流動可壓縮性和熱傳導的壁面函數(shù)耦合到了k-ω兩方程模型中；文獻[11]和文獻[12]分別通過數(shù)值實驗和風洞實驗對壁面函數(shù)的系數(shù)進行了修正研究；文獻[13]對文獻[8]中的壁面函數(shù)進行了修正研究。由于壁面律在存在逆壓梯度的復雜流動中是否適用無法確定[14-15]，文獻[8]忽略對流項和壓力梯度影響的壁面函數(shù)對于分離點和再附點附近流動的模擬會存在困難。

除了標準壁面函數(shù)，由于近年來解析壁面函數(shù)[16]在壁面處不涉及太多假設，保留了邊界層方程中的對流項和壓力梯度項，因此在不可壓縮流體中已經(jīng)得到應用。由于解析壁面函數(shù)在粗網(wǎng)格上的預測精度可以接近低雷諾數(shù)模型的結果，而計算時間比低雷諾數(shù)低一到兩個數(shù)量級。同時由于解析壁面函數(shù)在實際編程，計算中的魯棒性和可操作性等優(yōu)點，都使得其工程應用得到廣泛關注。而將解析壁面函數(shù)用于可壓縮流體中時，由于需要考慮流體可壓縮性，其方法研究對于先進壁面函數(shù)在復雜超聲速和高超聲速流動中的應用也具有實際意義。本文基于OpenFoamV 5.0軟件平臺將目前應用于不可壓縮流動的解析壁面函數(shù)，通過可壓縮湍流邊界層控制方程的簡化，考慮密度變化、對流項變化和黏性耗散項對解析壁面函數(shù)的影響，探索其在可壓縮流動中的應用，并通過二維超聲速激波邊界層干擾算例進行驗證。

1 物理模型

本文解算器使用Favre平均NS方程，無黏通量由Roe-Pike格式計算，黏性通量采用二階中心差分格式，時間推進采用一階歐拉格式。根據(jù)完全氣體假設，狀態(tài)方程為p=ρRT，空氣黏性系數(shù)使用Sutherland方程計算。壁面函數(shù)方法基于高雷諾數(shù)RANS模型，即標準k-ε模型。對于存在壁面的流動問題，局部雷諾數(shù)Rt=k2/ν/ε很小，黏性效應不能忽略，因此文獻[17]采用基于Rt的阻尼函數(shù)來計算近壁面湍流尺度，并在ε方程中添加Yap修正[18]，來降低在分離流動中ε方程得到的近壁面湍流，即低雷諾Launder-Sharmak-ε模型。

2 壁面函數(shù)

2.1 標準壁面函數(shù)(Standard Wall Function)

在局部平衡湍流邊界層內，壁面定律為：

(1)

其中:U+和y+定義為：

(2)

(3)

(4)

(5)

其中:常數(shù)cl=2.55，下標n表示壁面第一次網(wǎng)格面,如圖1所示。近壁面網(wǎng)格點P的壁面剪切應力為：

(6)

2.2 解析壁面函數(shù)(Analytical Wall Function)

解析壁面函數(shù)[16]在壁面附近通過可壓縮湍流邊界層假設簡化輸運方程，考慮對流項和壓力梯度的影響，壁面附近x-y平面上簡化的動量和能量方程為：

(7)

(8)

其中:P表示壓力;T表示溫度;μ和μt分別表示層流黏性系數(shù)和湍流黏性系數(shù);Pr和Prt分別表示層流普朗特數(shù)和湍流普朗特數(shù)。湍流黏性系數(shù)為：

μt=0 當y

(9)

(10)

(11)

其中:C表示動量簡化式(7)中的壓力梯度和對流項。

(12)

(13)

(14)

(15)

(16)

湍動能方程中的在壁面附近生成項為：

(17)

對于壁面附近的耗散項，在全湍流區(qū)為:

(18)

在黏性底層內耗散項為：

(19)

為了消除在黏性底層交界面處的間斷(圖1(a))，重新選取一個位置(使用下標d表示)，令耗散項在壁面網(wǎng)格內連續(xù)分布(圖1(b))，需滿足條件：

圖1 耗散項在壁面網(wǎng)格內分布示意圖Fig.1 Schech map of dissipation term in the wall cell

(20)

因此，在壁面網(wǎng)格內的平均耗散項為：

(21)

對于壁面溫度的解析表達式可以通過類似于壁面速度的兩次積分得到，對于絕熱壁有：

(22)

其中，Cth下標1表示黏性底層的值，2表示全湍流區(qū)的值：

等溫壁為：

(23)

原始解析壁面函數(shù)詳情可見文獻[16]。

2.3 可壓縮修正(Modification of AWF)

可壓縮流動的能量方程通常可表示為：

(24)

式中:E表示總能,σ表示流體應力張量,q表示熱流。

為了增強壁面函數(shù)的魯棒性，本文發(fā)展的解析壁面函數(shù)和主網(wǎng)格的能量方程統(tǒng)一使用式(24)，對于可壓縮湍流邊界黏性耗散項不能忽略[19]，因此在解析壁面函數(shù)能量方程中需包含粘性耗散項。

湍流壁面函數(shù)方法使用粗網(wǎng)格進行數(shù)值計算，第一層網(wǎng)格壁面距離通常滿足20

(25)

(26)

式(26)中右端第二項表示黏性耗散項，使用式(25)積分得到的解析速度計算，方程中的對流項和壓力梯度分別表示為：

上式D、C、Dth中下標1表示黏性底層的值，2表示全湍流區(qū)的值,在黏性底層內和全湍流區(qū)對簡化動量式(25)和能量式(26)積分，可得到解析速度和溫度表達式，具體形式如下：

(27)

(28)

其中：

系數(shù)N表示：

黏性底層溫度解析表達式為：

(29)

其中：

全湍流區(qū)溫度解析表達式為：

(30)

其中：

通過壁面第一層網(wǎng)格內的解析速度和溫度分布可以得到壁面剪切應力和壁面熱流。同時考慮到在使用粗網(wǎng)格求解時，為了更加準確描述主網(wǎng)格控制方程的黏性耗散項，即主方程中壁面網(wǎng)格內的平均粘性耗散項為：

在程序中通過解析速度表達式(27)、(28)數(shù)值求解。

3 數(shù)值結果

本文使用斜射激波-湍流邊界層干擾來驗證發(fā)展的可壓縮修正解析壁面函數(shù)(MAWF)，三種來流馬赫數(shù)實驗來流條件列于表1中[20-22]。來流邊界為充分發(fā)展的湍流邊界層，其邊界層厚度δ0和動量厚度θ0如表1所示，具有1.5%的自由來流湍流強度和黏性系數(shù)比率μt/μ=10。

表1 斜激波邊界層干擾來流條件Table 1 Inflow conditions for the impinging shock interactions

計算中使用結構網(wǎng)格，對于低雷諾數(shù)湍流模型，壁面第一層網(wǎng)格較密，第一層網(wǎng)格的y+≈1,對于高雷諾數(shù)湍流模型，第一層網(wǎng)格y+大約為30。為保證計算結果與網(wǎng)格無關，在開展結果對比之前，分別對表1中算例開展了網(wǎng)格無關性研究，計算結果表明所有算例均取得網(wǎng)格無關性結果，詳情可參考文獻[23]，其中Ma=2.9算例標準壁面函數(shù)法粗網(wǎng)格結果如圖2所示,其中xI表示無黏性情況下斜激波入射位置坐標。

圖2 Ma=2.9斜激波邊界層干擾的網(wǎng)格收斂性研究Fig.2 Grid independency study for the Ma=2.9 impinging shock interaction

3.1 Ma=2.9斜激波邊界層干擾[20]

計算網(wǎng)格入口邊界的下部使用充分發(fā)展的湍流邊界層，入口邊界的上部使用無粘斜激波關系式得到的斜激波后參數(shù)，在壁面處，使用等溫無滑壁邊界條件，壁溫為Tw=271K。在頂部和出口邊界處，使用Neumann邊界條件。本算例密網(wǎng)格為200×90，粗網(wǎng)格為150×60。

圖3比較了10°和13°斜激波的壁面壓力、摩擦系數(shù)和壁面熱流，圖中湍流模擬方法分別為標準壁面函數(shù)(SWF)、原始解析壁面函數(shù)(AWF)、可壓縮修正解析壁面函數(shù)(MAWF)和低雷諾數(shù)Launder-Sharmak-ε模型(LS)。對于β=10°，所有模型均返回相似的壁面壓力，和實驗數(shù)據(jù)[20]吻合。從摩擦系數(shù)結果來看，在無粘激波反射點附近存在一個小的分離區(qū)。LS模型的結果與實驗結果非常接近，壁面函數(shù)不能復現(xiàn)分離區(qū)，其中SWF傾向于低估分離區(qū)上下游的摩擦系數(shù)。壁面熱流不同模型的結果差異較大，其中MAWF結果更接近于密網(wǎng)格LS模型。對于β=13°，所有模型預測的壁面壓力和摩擦系數(shù)結果與實驗值較為吻合，MAWF預測的壁面熱流和密網(wǎng)格LS模型較為接近，而SWF和AWF均不能準確預測壁面熱流，在分離區(qū)前后甚至出現(xiàn)熱流符號相反的現(xiàn)象。

圖3 Ma=2.9算例不同入射角壁面壓力(上)、摩擦系數(shù)(中)和壁面熱流(下)Fig.3 Wall pressure(top)，skin-friction(middle) and wall heat-flux(bottom) distribution for the Ma=2.9 case with different impinging angles

3.2 Ma=5.0斜激波邊界層干擾[21]

本算例網(wǎng)格類似于Ma=3.0算例，其入口邊界使用充分發(fā)展的可壓縮湍流邊界層，其中動量厚度如表1所示，在壁面處使用等溫無滑壁邊界條件，壁溫為Tw=300 K。本算例密網(wǎng)格為240×80，粗網(wǎng)格為120×45。在相同設置的情況下，不同湍流模擬方法消耗的計算機時間對比如圖4所示，通過對比表明使用粗網(wǎng)格的壁面函數(shù)將大幅減少計算時間，僅為密網(wǎng)格低雷諾數(shù)模型計算時間的5%左右，一是由于網(wǎng)格量大幅減少，二是由于壁面距離的增加，在相同CFL數(shù)時計算時間步長的增加，收斂速度加快。

圖4 Ma=5.0算例不同模型計算消耗時間對比Fig.4 The CPU time comparison using different turbulence models for the Ma=5.0 case

圖5比較了不同入射角的壁面壓力、摩擦系數(shù)和壁面熱流，從圖中可知四種數(shù)值模擬方法得到的分離區(qū)較試驗結果[21]較小。當激波邊界層干擾加強時，AWF得到的壁面摩擦系數(shù)和Stanton數(shù)在分離起始位置出現(xiàn)非物理振蕩，主要原因是AWF使用壁面第一個網(wǎng)格點信息計算對流項，如式(11)，由于分離區(qū)附近流向速度變化較大，對流項在整個壁面網(wǎng)格內取常數(shù)，將導致AWF預測能力降低。而MAWF認為對流項在壁面處為零，以無量綱壁面距離的平方項遞增，從物理上更符合邊界層內流動規(guī)律，從而提高預測精度。對于三種不同角度斜激波，MAWF相較于AWF，壁面熱流的預測能力大幅上升，主要得益于能量方程簡化中保留了黏性耗散項，使得粗網(wǎng)格預測的壁面熱流精度大幅改善，接近密網(wǎng)格LS的結果。

3.3 Ma=7.2斜激波邊界層干擾[22]

本算例由于激波生成器拐角處的膨脹波會對下表面產生影響，因此在數(shù)值模擬的過程中加入了激波產生器的固壁邊界，激波產生器夾角β=7.5°和15°，且由于本算例風洞試驗模型是軸對稱體，因此z方向按照試驗外形旋轉2°，其余邊界設置和Ma=2.9和Ma=5.0算例相同，壁面溫度為300K。本算例密網(wǎng)格為420×105，粗網(wǎng)格為280×40。

圖6給出了四種不同湍流模型的數(shù)值結果與試驗值[22]的對比，其中MAWF大幅度提升了壁面熱流的預測精度，在分離區(qū)前后得到和密網(wǎng)格LS模型相同的結果，且避免了AWF在分離區(qū)起始位置附近的非物理振蕩。和實驗值相比，數(shù)值結果高估了分離區(qū)的大小和分離區(qū)內的壁面壓力、摩擦系數(shù)和壁面熱流，這主要是RANS模型的局限性導致，文獻[24]總結了不同RANS模型預測的β=15°壁面熱流結果，所有模型均出現(xiàn)高估分離區(qū)熱流的情況。

圖6 Ma=7.2算例不同入射角的壁面壓力(上)、摩擦系數(shù)(中)和壁面熱流(下)Fig.6 Wall pressure(top)， skin-friction(centre) and Stanton number (bottom) distribution for the Ma=7.2 case with different impinging angles

4 結 論

本文通過可壓縮湍流邊界層特點，考慮邊界層內密度的變化，將基于不可壓縮流動的解析壁面函數(shù)推廣到可壓縮流動原始解析壁面函數(shù)(AWF)，并考慮邊界層內對流項變化和粘性耗散項的影響，構造了適用于可壓縮修正的解析壁面函數(shù)(MAWF)，通過三個不同來流馬赫數(shù)的斜激波邊界層干擾算例進行測試，并與LS和SWF進行了比較。主要結論如下：

1)整體來看四種數(shù)值模擬方法預測的壁面壓力差異不大，且與實驗吻合良好。

2)整體來看四種數(shù)值模擬方法預測的壁面摩擦系數(shù)差異不大，標準壁面函數(shù)預測的數(shù)值偏小，可壓縮解析壁面函數(shù)由于考慮了對流項在邊界層內分布，消除了原始解析壁面函數(shù)在分離區(qū)起始位置的非物理振蕩。

3)壁面熱流的數(shù)值結果差異較大，整體來看本文構造的可壓縮解析壁面函數(shù)得到的結果和密網(wǎng)格Launder-Sharmak-ε模型結果最為接近，而原始解析壁面函數(shù)和標準壁面函數(shù)，不能準確預測壁面熱流。

4)壁面函數(shù)模型使用粗網(wǎng)格，相較于低雷諾數(shù)模型可以大幅減少計算時間，對于Ma=5算例，計算時間約為Launder-Sharmak-ε模型的5.3%，而相較與標準壁面函數(shù)僅增加了1%左右的時間。

總的來說，四種方法都顯示出良好的壁面壓力和摩擦系數(shù)預測能力。然而，標準壁面函數(shù)和原始解析壁面函數(shù)嚴重低估壁面熱流，即使在分離區(qū)前后其結果也和密網(wǎng)格Launder-Sharmak-ε模型結果相差甚遠。本文構造的解析壁面函數(shù)充分考慮可壓縮湍流邊界層流動的特征，彌補了原始解析壁面的缺點，從而大幅提升了壁面熱流的預測精度，且不需要在近壁面區(qū)域內的精細網(wǎng)格。