王新光，陳 琦，萬(wàn) 釗，張愛婧，陳堅(jiān)強(qiáng)

(1. 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621010；2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心計(jì)算所，綿陽(yáng) 621000)

0 引 言

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

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

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

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

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

1 物理模型

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

2 壁面函數(shù)

2.1 標(biāo)準(zhǔn)壁面函數(shù)(Standard Wall Function)

在局部平衡湍流邊界層內(nèi)，壁面定律為：

(1)

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

(2)

(3)

(4)

(5)

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

(6)

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

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

(7)

(8)

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

μt=0 當(dāng)y

(9)

(10)

(11)

其中:C表示動(dòng)量簡(jiǎn)化式(7)中的壓力梯度和對(duì)流項(xiàng)。

(12)

(13)

(14)

(15)

(16)

湍動(dòng)能方程中的在壁面附近生成項(xiàng)為：

(17)

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

(18)

在黏性底層內(nèi)耗散項(xiàng)為：

(19)

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

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

(20)

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

(21)

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

(22)

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

等溫壁為：

(23)

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

2.3 可壓縮修正(Modification of AWF)

可壓縮流動(dòng)的能量方程通?？杀硎緸椋?/p>

(24)

式中:E表示總能,σ表示流體應(yīng)力張量,q表示熱流。

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

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

(25)

(26)

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

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

(27)

(28)

其中：

系數(shù)N表示：

黏性底層溫度解析表達(dá)式為：

(29)

其中：

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

(30)

其中：

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

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

3 數(shù)值結(jié)果

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

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

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

圖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]

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

圖3比較了10°和13°斜激波的壁面壓力、摩擦系數(shù)和壁面熱流，圖中湍流模擬方法分別為標(biāo)準(zhǔn)壁面函數(shù)(SWF)、原始解析壁面函數(shù)(AWF)、可壓縮修正解析壁面函數(shù)(MAWF)和低雷諾數(shù)Launder-Sharmak-ε模型(LS)。對(duì)于β=10°，所有模型均返回相似的壁面壓力，和實(shí)驗(yàn)數(shù)據(jù)[20]吻合。從摩擦系數(shù)結(jié)果來(lái)看，在無(wú)粘激波反射點(diǎn)附近存在一個(gè)小的分離區(qū)。LS模型的結(jié)果與實(shí)驗(yàn)結(jié)果非常接近，壁面函數(shù)不能復(fù)現(xiàn)分離區(qū)，其中SWF傾向于低估分離區(qū)上下游的摩擦系數(shù)。壁面熱流不同模型的結(jié)果差異較大，其中MAWF結(jié)果更接近于密網(wǎng)格LS模型。對(duì)于β=13°，所有模型預(yù)測(cè)的壁面壓力和摩擦系數(shù)結(jié)果與實(shí)驗(yàn)值較為吻合，MAWF預(yù)測(cè)的壁面熱流和密網(wǎng)格LS模型較為接近，而SWF和AWF均不能準(zhǔn)確預(yù)測(cè)壁面熱流，在分離區(qū)前后甚至出現(xiàn)熱流符號(hào)相反的現(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ā)展的可壓縮湍流邊界層，其中動(dòng)量厚度如表1所示，在壁面處使用等溫?zé)o滑壁邊界條件，壁溫為Tw=300 K。本算例密網(wǎng)格為240×80，粗網(wǎng)格為120×45。在相同設(shè)置的情況下，不同湍流模擬方法消耗的計(jì)算機(jī)時(shí)間對(duì)比如圖4所示，通過對(duì)比表明使用粗網(wǎng)格的壁面函數(shù)將大幅減少計(jì)算時(shí)間，僅為密網(wǎng)格低雷諾數(shù)模型計(jì)算時(shí)間的5%左右，一是由于網(wǎng)格量大幅減少，二是由于壁面距離的增加，在相同CFL數(shù)時(shí)計(jì)算時(shí)間步長(zhǎng)的增加，收斂速度加快。

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

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

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

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

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

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

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

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

3)壁面熱流的數(shù)值結(jié)果差異較大，整體來(lái)看本文構(gòu)造的可壓縮解析壁面函數(shù)得到的結(jié)果和密網(wǎng)格Launder-Sharmak-ε模型結(jié)果最為接近，而原始解析壁面函數(shù)和標(biāo)準(zhǔn)壁面函數(shù)，不能準(zhǔn)確預(yù)測(cè)壁面熱流。

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

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