可壓縮湍流邊界層壁面函數(shù)方法綜述

毛枚良，閔耀兵，王新光，*，陳 琦，葉 濤

(1. 中國空氣動(dòng)力研究與發(fā)展中心，綿陽 621000；2. 西南科技大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，綿陽 621010)

0 引 言

盡管近幾十年來計(jì)算機(jī)資源飛速增長，但是巨大的計(jì)算資源消耗仍然是限制直接數(shù)值模擬和大渦模擬在工程復(fù)雜湍流問題上應(yīng)用的一個(gè)根本因素。即使是采用計(jì)算資源消耗最少的RANS模擬方法，要獲得準(zhǔn)確的壁面摩阻和熱流，通常仍然需要在黏性底層中布置一定數(shù)目的網(wǎng)格點(diǎn)，從而使得壁面附近的網(wǎng)格十分細(xì)密。這樣，不僅會(huì)極大地增加迭代收斂的計(jì)算步數(shù)，還會(huì)帶來較為嚴(yán)重的數(shù)值剛性問題，導(dǎo)致迭代計(jì)算的穩(wěn)定性下降。湍流邊界層模擬采用壁面函數(shù)方法，可以大幅放寬壁面第一層網(wǎng)格的尺度，從而達(dá)到加速計(jì)算過程收斂的效果。另外，壁面函數(shù)方法不僅可以應(yīng)用于RANS模擬，也可以用于大渦模擬[1]。因此，目前主流商業(yè)軟件在工程湍流模擬應(yīng)用中默認(rèn)采用壁面函數(shù)方法。

圖1給出了湍流邊界層沿壁面法向的典型速度分布，在緊挨壁面的黏性底層中，由于受到壁面的限制，速度脈動(dòng)通常小到可以忽略，流動(dòng)黏性中分子黏性占主導(dǎo)，湍流渦黏性可以忽略；在對(duì)數(shù)律層和外層，湍流脈動(dòng)得到充分發(fā)展，流動(dòng)黏性由湍流渦黏性主導(dǎo)，分子黏性效應(yīng)可以忽略；在黏性底層和對(duì)數(shù)律層之間的過渡層，分子黏性和湍流黏性同等重要，兩者共同影響速度分布。壁面函數(shù)基于如下的事實(shí)：對(duì)于不可壓縮流動(dòng)，從壁面到對(duì)數(shù)律層外緣之間在局部平衡的邊界層中存在相似解。由于最初使用“壁面律（law of the wall）”僅指普適的對(duì)數(shù)律層解，而“混合壁面律（hybrid law of the wall）”或“自適應(yīng)壁面律（adaptive law of the wall）”則主要強(qiáng)調(diào)在整個(gè)近壁區(qū)域內(nèi)一致有效。對(duì)于不可壓縮平板湍流邊界層而言，在Prandtl混合長度理論和若干假設(shè)下，利用若干實(shí)驗(yàn)觀察/測量結(jié)果，可以證明平均速度的壁面函數(shù)就是邊界層方程（RANS方程在Prandtl邊界層假設(shè)下的近似）的近似解[2]。

圖1 湍流邊界層示意圖Fig. 1 A schematic diagram of turbulent boundary layer

對(duì)于管道和平板湍流邊界層的分析，大多數(shù)早期的壁面函數(shù)并不考慮壓力梯度和傳熱的影響。1938年Millikan[3]給出了剪切力占主導(dǎo)作用的不可壓縮較大Reynolds數(shù)流動(dòng)的對(duì)數(shù)律層的漸近解，然而在流動(dòng)分離點(diǎn)和再附點(diǎn)附近壁面剪切力為零，故Millikan的漸近解在分離點(diǎn)和再附點(diǎn)附近不再成立。Tennekes和Lumley[4]推導(dǎo)了考慮壓力梯度和零剪應(yīng)力情況下邊界層的漸近解，但Tennekes和Lumley的解不能適用于壓力梯度為零的情況。Nichols[5]、Viegas[6]和Smith[7]等都相繼發(fā)展了包含壁面?zhèn)鳠嵝?yīng)的壁面函數(shù)，但這些壁面函數(shù)都不是統(tǒng)一的速度分布函數(shù)，在黏性底層與對(duì)數(shù)律層之間的過渡區(qū)域也是不連續(xù)的，Wilcox[8]引入一種壁面匹配方法?；赟palding[9]統(tǒng)一壁面函數(shù)理論，F(xiàn)rink[10]發(fā)展了絕熱壁的壁面函數(shù)，Shih等[11]基于Spalding的工作也給出了一個(gè)包含壓力梯度效應(yīng)的統(tǒng)一壁面函數(shù),但仍然沒有將流動(dòng)可壓縮性和壁面?zhèn)鳠嵝?yīng)耦合進(jìn)來。Nichols和Nelson[12]在White和Christoph[13]工作的基礎(chǔ)上，忽略壓力梯度的影響，將速度和溫度函數(shù)耦合起來，發(fā)展了適用于不可壓縮流動(dòng)、可壓縮流動(dòng)、絕熱壁和等溫壁湍流邊界層的壁面函數(shù)。由于壁面律在壓力梯度不可忽略的區(qū)域內(nèi)是否仍然有效尚無定論[14-15]，Nichols和Nelson忽略壓力梯度影響的壁面函數(shù)對(duì)于分離點(diǎn)和再附點(diǎn)附近流動(dòng)的模擬會(huì)存在困難。

實(shí)際上，考慮到壁面函數(shù)的使用能有效減少湍流邊界層數(shù)值模擬結(jié)果對(duì)壁面網(wǎng)格尺度的關(guān)聯(lián)性，RANS模擬中湍流模型在具有壓力梯度的流動(dòng)中表現(xiàn)非常重要。Knopp等[16]指出：盡管壁面函數(shù)的研究工作已經(jīng)持續(xù)開展了50多年，其中大多數(shù)研究成果的應(yīng)用工作并不樂觀，如壁面函數(shù)與用于全局流動(dòng)問題模擬的湍流模型不相容的問題[10,17-18]，容易引起數(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)物理特性恰當(dāng)?shù)姆直?，特別是對(duì)非平衡效應(yīng)顯著的流動(dòng)駐點(diǎn)區(qū)域、逆壓梯度導(dǎo)致的流動(dòng)分離/再附點(diǎn)附近區(qū)域，以獲得比較準(zhǔn)確的氣動(dòng)力/熱特性。此外，Kalitzin等[19]也針對(duì)壁面函數(shù)在RANS模擬中的應(yīng)用問題開展了相關(guān)研究工作。

從公開的文獻(xiàn)來看，關(guān)于湍流邊界層壁面函數(shù)的研究，原創(chuàng)性工作幾乎都來自國外，國內(nèi)相關(guān)工作較少，比較典型的研究工作有：賀旭照等[20-21]將Nichols的考慮流動(dòng)可壓縮性和熱傳導(dǎo)的壁面函數(shù)耦合到了k-ω兩方程模型中；吳曉軍[22]和肖紅林[23]分別在RANS模擬和大渦模擬中采用壁面函數(shù)；Gao[24]和Tao等[25]分別通過數(shù)值實(shí)驗(yàn)和風(fēng)洞實(shí)驗(yàn)對(duì)壁面函數(shù)的系數(shù)進(jìn)行了修正研究；Wu等[26]對(duì)Nichols的壁面函數(shù)進(jìn)行了修正研究。目前可壓縮湍流特別是強(qiáng)壓縮和顯著傳熱的湍流理論研究，幾乎都是對(duì)不可壓湍流研究成果的修正、延拓和推廣。國內(nèi)開展高超聲速湍流邊界層壁面函數(shù)的研究工作相對(duì)較少，缺乏系統(tǒng)性，也缺乏原創(chuàng)性。

本文的目的是對(duì)現(xiàn)有研究工作進(jìn)行梳理，從理論上比較系統(tǒng)地把握湍流邊界層壁面函數(shù)理論，為研制自主可控的高超聲速工程實(shí)用的CFD軟件提供理論支撐。

1 標(biāo)準(zhǔn)壁面函數(shù)

二維可壓縮定常平板湍流邊界層問題，通過量級(jí)分析，忽略壓力在邊界層內(nèi)沿壁面法向的變化，可以得到湍流邊界層內(nèi)的流動(dòng)控制方程為：

其中，x為來流方向，y為壁面法向。由x方向的動(dòng)量方程（2）可知總剪切力在邊界層內(nèi)沿壁面法向保持不變。由于在黏性底層速度脈動(dòng)被壁面抑制，湍流渦黏性可以忽略不計(jì)，記：

在遠(yuǎn)離壁面一定距離后，湍流脈動(dòng)充分發(fā)展，湍流渦黏性遠(yuǎn)大于分子黏性，忽略分子黏性的影響可以得到：

對(duì)于不可壓流動(dòng)，密度ρ為常數(shù)，對(duì)壁面剪切力定義式（4）在壁面附近積分，并考慮到速度在壁面處的無滑移邊界條件可以得到：

其中，u+=u/uτ為無量綱速度，y+=ρyuτ/μ為無量綱壁面距離，為摩擦速度。為黏性底層外緣處的無量綱常數(shù)，主要由實(shí)驗(yàn)觀察結(jié)果確定，不同的文獻(xiàn)略有差異，大致取值范圍為∈(5,15)。

在湍流脈動(dòng)充分發(fā)展的對(duì)數(shù)律層內(nèi)，依據(jù)Prandtl混合長度理論給出：

其中，lm=κy為對(duì)數(shù)律層內(nèi)的Prandtl混合長度，κ≈0.41，為von Kármán常數(shù)。考慮到湍流脈動(dòng)充分發(fā)展時(shí)可忽略分子黏性的影響（式（5）），容易得到：

對(duì)式（8）進(jìn)行積分可以得到對(duì)數(shù)律表達(dá)式為：

其中常數(shù)B∈[5,5.5]。

由上述分析可知，壁面函數(shù)實(shí)質(zhì)上就是邊界層方程的近似解。需要特別指出的是：τw不僅是壁面處的分子黏性應(yīng)力，也等于邊界層剖面中的總黏性應(yīng)力，即湍流渦黏性應(yīng)力與分子黏性應(yīng)力之和），故τw既可以通過壁面的分子黏性應(yīng)力公式（4）計(jì)算得到，也可以通過對(duì)數(shù)律層中的湍流渦黏性應(yīng)力計(jì)算得到。

在靠近壁面的剪切流區(qū)域，其流動(dòng)特性與Couette流動(dòng)相似，在湍動(dòng)能生成與耗散處于局部平衡的條件下，壁面剪切應(yīng)力τw與湍動(dòng)能k之間滿足關(guān)系：

其中Cμ=0.09為常數(shù)，于是平均速度剖面可以改寫為：

對(duì)于可壓縮流動(dòng)，特別是高超聲速流動(dòng)，溫度和密度在邊界層內(nèi)存在明顯的變化，導(dǎo)致分子黏性系數(shù)在邊界層內(nèi)的變化不可忽略。Fernholz和Finley[27]指出，若平均速度采用變換：

則可壓縮流動(dòng)的對(duì)數(shù)律可以表述為[28]：

或者表示為：

對(duì)于黏性底層，有不少工作忽略可壓縮性的影響，仍然采用不可壓縮壁面函數(shù)的表達(dá)式。實(shí)際上，對(duì)于高超聲速流動(dòng)，黏性底層中密度變化可能是不可忽略的，圖2給出了不同馬赫數(shù)情況下平板邊界層[29-31]密度分布，可以看出，隨著馬赫數(shù)的增加，黏性底層內(nèi)的密度變化愈加明顯。當(dāng)來流馬赫數(shù)為11.1時(shí)，黏性底層頂部附近，密度變化接近40%。根據(jù)邊界層方程，類似地引入變換：

則可壓縮流動(dòng)的黏性底層速度分布可表述為：

圖2 高超聲速平板邊界層密度分布Fig. 2 The density distributions of hypersonic boundary layers

對(duì)于空氣而言，μ/μw是T/Tw的函數(shù)，在邊界層內(nèi)壓力沿壁面法向保持不變，由狀態(tài)方程可知ρ/ρw也僅是T/Tw相關(guān)，只要獲得了溫度在邊界層內(nèi)的表達(dá)式，就可以唯一確定可壓縮流動(dòng)條件下的平均速度壁面函數(shù)。

溫度在邊界層內(nèi)的分布可由邊界層能量方程得到，簡略推導(dǎo)如下：將邊界層能量方程（3）沿邊界層法向積分，得到：

其中用到了總剪切應(yīng)力在湍流邊界層內(nèi)層沿壁面法向不變的條件，qw=-(k?T/?y)y=0為壁面熱流，k為傳熱系數(shù)。在黏性底層，分子黏性起主要作用，而在對(duì)數(shù)律層，則湍流黏性起主要作用，故有：

在黏性底層（0＜y+＜） 考慮傳熱系數(shù)（18）和速度分布（6），由能量積分方程（17）可得黏性底層（0＜y+＜） 中的溫度分布為：

在對(duì)數(shù)律層湍流脈動(dòng)充分發(fā)展，湍流渦黏性起主要作用，可忽略分子黏性的影響即：

考慮傳熱系數(shù)式（18），由能量積分方程（17）可得對(duì)數(shù)律層的溫度關(guān)系為：

從對(duì)數(shù)律層下緣開始對(duì)式（21）積分得到：

湍流邊界層沿壁面法向的剖面包括黏性底層、過渡層、對(duì)數(shù)律層以及邊界層外層等多層結(jié)構(gòu)，對(duì)應(yīng)地，湍流邊界層的壁面函數(shù)通常描述從壁面到完全湍流區(qū)域的剖面。因此，嚴(yán)格來說應(yīng)該采用三層模型，但有不少研究者忽略對(duì)過渡層的模擬，而采用兩層模型。對(duì)于速度壁面函數(shù)，當(dāng)y+＜11.225時(shí)，采用黏性底層公式（6）計(jì)算，否則，采用對(duì)數(shù)律層公式（9）描述。而對(duì)于溫度壁面函數(shù)，由于其分布還同壁面溫度等相關(guān)。實(shí)際上式（19）和式（22）共同構(gòu)成了一個(gè)兩層的溫度壁面函數(shù)，在交界點(diǎn)y+=處，采用黏性底層公式（19）和對(duì)數(shù)律層公式（22）計(jì)算得到的無量綱溫度應(yīng)相同，則有：

Huang等[32]假定黏性底層的等效Prandtl數(shù)與湍流Prandtl數(shù)相同，則式（19）和式（22）相同。從已有的數(shù)值試驗(yàn)結(jié)果來看，對(duì)速度分布進(jìn)行可壓縮修正時(shí)，這種假定對(duì)壁面摩擦應(yīng)力造成的誤差可以忽略，且簡化了計(jì)算，但對(duì)等溫壁壁面熱流的影響如何，則需要進(jìn)一步研究。

2 統(tǒng)一壁面函數(shù)

即使采用兩層模型，在實(shí)際使用過程中也不太方便，且易造成平均速度分布不連續(xù)。Splading[9]提出了在黏性底層和對(duì)數(shù)律層內(nèi)一致有效的絕熱不可壓流動(dòng)的湍流邊界層壁面函數(shù)：

Nichols[12]在考慮傳熱和可壓縮性修正時(shí)，采用White[13]的對(duì)數(shù)層壁面率公式（14）來替換Splading公式（23）中的絕熱不可壓壁面率，得到：

其中：

一個(gè)比較自然的想法，將式（14）和式（16）類似于式（23）進(jìn)行合并，可以得到：

圖3比較了統(tǒng)一壁面函數(shù)的不同表達(dá)式（23）、式（24）和式（25）的速度型，當(dāng)來流馬赫數(shù)為11.1時(shí)，在黏性底層內(nèi)的溫度、密度變化較大，整體考慮可壓縮效應(yīng)的統(tǒng)一壁面函數(shù)公式（25）更接近于數(shù)值解，而在來流馬赫數(shù)為5時(shí)三種表達(dá)形式差異不大。如果忽略黏性底層溫度變化對(duì)速度分布的影響，則公式（25）退化到不可壓縮絕熱邊界層的表達(dá)式。在不可壓縮絕熱流動(dòng)時(shí)，式（25）也和式（24）一樣，能夠自動(dòng)退化到式（23），同時(shí)還比較全面地考慮了可壓縮性和傳熱對(duì)黏性底層和對(duì)數(shù)律層的影響。

圖3 典型高超聲速邊界層平均速度分布 = 3.7×107 m-1, T∞ = 68.79 K, Tw = 300 KFig. 3 The mean velocity profiles of typical hypersonic boundary layers ReL = 3.7×107 m-1,T∞ = 68.79 K, Tw = 300 K

在可壓縮邊界層中，溫度分布函數(shù)一般采用Crocco-Busemann公式：

式（19）、式（22）和式（26）的比較如圖4所示，當(dāng)Ma∞= 5時(shí)，黏性底層內(nèi)大體相同，僅略有差異，而當(dāng)Ma∞= 11.1時(shí)差異較大，整體來看，公式（19）和式（22）在黏性底層和對(duì)數(shù)律層分別取層流和湍流普朗特?cái)?shù)時(shí)，更接近于數(shù)值模擬結(jié)果。

圖4 典型高超聲速邊界層溫度分布 = 3.7×107 m-1, T∞ = 68.79 K, Tw = 300 KFig. 4 The mean temperature profiles of typical hypersonic boundary layers ReL = 3.7×107 m-1,T∞ = 68.79 K, Tw = 300 K

此外，Kader[33]也提出了一個(gè)混合方法，將黏性底層與對(duì)數(shù)律層的分布函數(shù)整合為一個(gè)表達(dá)式，以便統(tǒng)一描述近壁區(qū)域（黏性底層、過渡層、完全湍流區(qū)）的物理量分布。

3 增強(qiáng)型壁面函數(shù)

上述壁面函數(shù)的推導(dǎo)，均未涉及壓力梯度。實(shí)際飛行器的邊界層流動(dòng)，由于壁面通常是曲面，并可能出現(xiàn)凸起和凹坑，從而使得流動(dòng)具有壓力梯度。眾所周知，流場中的逆壓梯度和順壓梯度對(duì)壁面流動(dòng)的影響是非常重要的。逆壓梯度不僅會(huì)使流動(dòng)減速，甚至可能造成流動(dòng)從壁面分離，使得壁面附近的流動(dòng)失去了經(jīng)典意義下的邊界層特性；而順壓梯度則會(huì)加速邊界層中的流動(dòng)，也可能會(huì)壓制湍流脈動(dòng)，使湍流邊界層流動(dòng)層流化。增強(qiáng)型壁面函數(shù)設(shè)計(jì)的一個(gè)重要考量就是壓力梯度的影響，使之能夠適應(yīng)流動(dòng)分離和再附點(diǎn)附近壁面流動(dòng)的模擬。對(duì)于逆壓梯度和零剪應(yīng)力邊界層，Tennekes[4]推導(dǎo)了一個(gè)漸近解：

其中：

常數(shù)α≈5.0,β≈8由Stratford[34]的實(shí)驗(yàn)數(shù)據(jù)確定。顯然，Tennekes的漸近解(式27)在壓力梯度趨于零時(shí)不再成立，因?yàn)榇藭r(shí)速度up也趨于零。而在流動(dòng)分離和再附點(diǎn)附近，剪切應(yīng)力和摩擦速度趨于零，基于剪切應(yīng)力特性的壁面函數(shù)（式（7））亦不再成立。

Shih等[11]基于Tennekes的分解方法，將不可壓縮流動(dòng)邊界層方程的積分形式：

分裂為：

其中：

利用量綱分析方法，容易得到漸近解：

其中：

在黏性底層中，湍流應(yīng)力可以忽略，由式（28）得到速度剖面：

由此可知，邊界層流動(dòng)由壁面剪切應(yīng)力和壁面壓力梯度決定。不再會(huì)出現(xiàn)uτ+p= 0的情況，因?yàn)閡τ和up不會(huì)同時(shí)為零。對(duì)于壓力梯度起主要作用的流動(dòng)，沿邊界層法向黏性剪應(yīng)力不再是常數(shù)；對(duì)于零壓力梯度的情況，流動(dòng)由壁面剪應(yīng)力確定，邊界層中黏性應(yīng)力沿壁面法向?yàn)槌?shù)，且等同于壁面剪切應(yīng)力。

式(31)和式(32)分別為對(duì)數(shù)律層和黏性底層的漸近解，兩層之間的區(qū)域（5≤uτ+py/v≤30）稱為過渡層（buffer layer），在過渡層中，黏性應(yīng)力和湍流應(yīng)力處于同一量級(jí)，尚未有文獻(xiàn)在理論上給出過渡層的漸近解。在過渡層中采用簡單的湍流應(yīng)力模型和基于量級(jí)分析的匹配過程，能夠得到唯一的解析函數(shù)，即適用于整個(gè)湍流邊界層（包括黏性底層、過渡層、對(duì)數(shù)律層）的統(tǒng)一壁面函數(shù)。Spalding[9]首先構(gòu)造了一個(gè)統(tǒng)一壁面函數(shù)（23）（沒有考慮壓力梯度的影響）。為以示區(qū)別，將Spalding的統(tǒng)一壁面函數(shù)（23）記為：

其逆函數(shù)記為：

當(dāng)壓力梯度占主導(dǎo)作用時(shí)的壁面函數(shù)為：

其逆函數(shù)記為：

結(jié)合函數(shù)（34）和式（36）可以得到滿足零壁面摩擦應(yīng)力和零壓力梯度的漸近解，即能夠應(yīng)用于流動(dòng)加速、減速、分離和回流區(qū)域的速度壁面函數(shù)為：

其中y+=uτ,py/v。在黏性底層和過渡層函數(shù)f1和f2分別為關(guān)于和的分片多項(xiàng)式函數(shù)，而在對(duì)數(shù)律層均為對(duì)數(shù)函數(shù)，其具體表達(dá)式由Shih等[11]給出。

對(duì)于可壓縮流動(dòng)，White[35]給出的對(duì)數(shù)律層平均速度分布函數(shù)為：

其中：

需要特別指出的是：由于White等在推導(dǎo)上式時(shí)采用的混合長表達(dá)式lm=κy，不能退化到絕熱不可壓對(duì)數(shù)律（31），在使用過程中一定要牢記這一點(diǎn)。如果壓力梯度為零，則式（38）可退化為：

此式等同于上文的式（13）。在黏性底層，忽略湍流黏性的影響由邊界層方程可以得到：

也可進(jìn)一步改寫為：

在傳熱、可壓縮性和壓力梯度的共同作用下，黏性底層的高度與絕熱不可壓流動(dòng)之間的差異，可能是一個(gè)需要重視的問題。

4 壁面函數(shù)的應(yīng)用

RANS方法在模擬湍流流動(dòng)時(shí)，采用壁面函數(shù)的初衷是在固壁處提供一個(gè)邊界條件，降低數(shù)值解（特別是摩擦阻力和熱流等）對(duì)臨近壁面第一層網(wǎng)格節(jié)點(diǎn)位置的依賴性。以表示第一個(gè)網(wǎng)格節(jié)點(diǎn)到固壁面的距離（采用黏性長度尺度uτ+p/v 進(jìn)行無量綱化），當(dāng)網(wǎng)格雷諾數(shù)較低時(shí)（≈1），RANS方程可以直接積分到壁面，在壁面處的邊界條件可以采用速度無滑移和相應(yīng)的溫度邊界條件；而當(dāng)網(wǎng)格雷諾數(shù)較高時(shí)30），此時(shí)第一層網(wǎng)格節(jié)點(diǎn)位于湍流充分發(fā)展的對(duì)數(shù)律層，只能在對(duì)數(shù)律層內(nèi)求解RANS方程，在壁面處必須采用壁面函數(shù)才能求得比較準(zhǔn)確的壁面剪切應(yīng)力、熱流和溫度（具體由壁面溫度邊界條件確定）等物理量；對(duì)于中等網(wǎng)格雷諾數(shù)問題，第一層網(wǎng)格節(jié)點(diǎn)位于黏性底層和對(duì)數(shù)律層之間的過渡層，壁面邊界的處理與高網(wǎng)格雷諾數(shù)情況類似。

在實(shí)際計(jì)算中，網(wǎng)格生成是無法預(yù)先確定壁面第一層網(wǎng)格中每個(gè)網(wǎng)格點(diǎn)位于湍流邊界層的具體位置，只有在計(jì)算過程中，根據(jù)其速度、溫度和壁面邊界條件，通過壁面函數(shù)才能確定。對(duì)于不可壓縮流動(dòng)，邊界層的具體位置可由第一層網(wǎng)格對(duì)應(yīng)的網(wǎng)格雷諾數(shù)唯一確定。因?yàn)楹愠闪?，因此，在黏性底層中y+=恒成立。據(jù)此可以依據(jù)Rey的 大小來判斷臨近壁面第一層網(wǎng)格節(jié)點(diǎn)處于邊界層的哪個(gè)子層中，同時(shí)壁面摩擦速度、剪切應(yīng)力等參數(shù)均可由Rey唯一確定?；诖苏J(rèn)識(shí)，有文獻(xiàn)給出了不可壓縮流動(dòng)的壁面摩擦速度、剪切應(yīng)力等參數(shù)與Rey的擬合函數(shù)，在具體計(jì)算中，無需對(duì)非線性方程組進(jìn)行迭代求解來得到壁面摩擦應(yīng)力等相關(guān)參數(shù)，以進(jìn)一步節(jié)省計(jì)算時(shí)間，提高計(jì)算效率。對(duì)于壓力梯度影響不可忽略的不可壓縮流動(dòng)，原則上，也可以構(gòu)造雙參數(shù)的擬合函數(shù)，但目前尚未看到相關(guān)的研究工作。對(duì)于可壓縮流動(dòng)，邊界層內(nèi)速度、溫度的分布所依賴的參數(shù)，除了第一層網(wǎng)格對(duì)應(yīng)的網(wǎng)格雷諾數(shù)和壓力梯度外，還有壁面溫度邊界條件參數(shù)、來流馬赫數(shù)和溫度等，由于相關(guān)參數(shù)較多，構(gòu)造擬合函數(shù)的工作量和難度會(huì)急劇增大。

湍流邊界層壁面函數(shù)的理論和實(shí)驗(yàn)研究工作，主要基于平板邊界層，特別是理論研究。實(shí)際飛行器的表面，普遍是三維曲面。將二維平板湍流邊界層的理論近似結(jié)果，推廣到三維曲面邊界層流動(dòng)，其可靠性需通過數(shù)值模擬試驗(yàn)同風(fēng)洞和實(shí)際飛行試驗(yàn)結(jié)果的一致性來評(píng)估。

從平面推廣到曲面，一個(gè)自然的作法，就是把邊界層方程及其相關(guān)理論表達(dá)式轉(zhuǎn)換到貼體坐標(biāo)系中。以當(dāng)?shù)厍邢蛩俣茸鳛榱鲃?dòng)方向，將三維流動(dòng)進(jìn)行局部二維化處理，在此基礎(chǔ)上采用湍流邊界層的壁面函數(shù)，可以得到壁面剪切應(yīng)力、熱流（或壁面溫度）等。對(duì)于壓力梯度影響可以忽略的流動(dòng)，這種做法似乎是一種比較自然的推廣，但對(duì)于壓力梯度與流動(dòng)方向存在明顯差異的流動(dòng)情況，如何推廣應(yīng)用已有的壁面函數(shù)研究成果還需要進(jìn)一步研究。

壁面函數(shù)在CFD計(jì)算過程中的具體實(shí)施方法，Tidriri[36]從計(jì)算區(qū)域分解和疊加求解的角度，給出了一個(gè)十分清晰的闡述，如圖5所示。

圖5 壁面附近計(jì)算區(qū)域分解[16]Fig. 5 The near-wall domain decomposition with full overlap[16]

記Ωδ?Ω為計(jì)算區(qū)域Ω的近壁區(qū)域，內(nèi)邊界位于對(duì)數(shù)層內(nèi)，那么壁面函數(shù)在CFD中的應(yīng)用問題可以分解為兩個(gè)問題：1）全局流動(dòng)計(jì)算問題，在整個(gè)計(jì)算區(qū)域Ω內(nèi)求解RANS方程，利用壁面函數(shù)改進(jìn)邊界條件（在固壁面強(qiáng)加剪應(yīng)力條件，而不是速度無滑移條件，等溫固壁面強(qiáng)加熱流條件而不僅僅給定壁面溫度），以彌補(bǔ)由于網(wǎng)格無法分辨湍流邊界層中的近壁流動(dòng)結(jié)構(gòu)產(chǎn)生的誤差；2）近壁區(qū)域Ωδ內(nèi)湍流邊界層的求解問題。顯然，解的具體形式就是湍流邊界層方程的近似解即壁面函數(shù)，在固壁面Γw上始終自動(dòng)滿足原始邊界條件（速度無滑移條件、等溫或絕熱壁面條件），在近壁區(qū)域Ωδ的外緣Γδ上與全局的RANS數(shù)值解匹配一致。該方法對(duì)Γδ處于任意位置都是適定的，所得的解與Γδ的具體位置無關(guān)（只要Γδ位于湍流邊界層的對(duì)數(shù)律層以內(nèi)的近壁流動(dòng)區(qū)域中），這意味著壁面壓力、摩阻和熱流等物理量與網(wǎng)格分布密度無關(guān)。目前實(shí)際計(jì)算中幾乎都設(shè)定Γδ為離開壁面的第一層網(wǎng)格形成的曲面。由此可見，在RANS模擬中采用壁面函數(shù)方法，其實(shí)質(zhì)是避免在近壁區(qū)域數(shù)值求解完全可壓縮RANS方程和湍流模型方程。

由上所述易知，不僅求解RANS方程需要用到壁面函數(shù)邊界條件，求解湍流模型方程也同樣需要。在實(shí)際工程應(yīng)用的航空航天飛行器湍流流動(dòng)的模擬中，最常用的湍流模型為以SA模型為代表的一方程模型和以SSTk-ω模型為代表的兩方程模型等。Nichols和Nelson[12]在忽略壓力梯度的情況下，基于剪切應(yīng)力沿邊界層法向基本不變的假設(shè)，利用壁面函數(shù)公式、臨近壁面第一層網(wǎng)格單元上的速度、溫度（或絕熱邊界條件）和壁面距離等物理流場和幾何信息，采用迭代方法獲得壁面摩擦應(yīng)力和熱流（或壁面溫度），并以此為基礎(chǔ)，修正臨近壁面網(wǎng)格單元上的湍流渦黏性系數(shù)和湍流模型方程中的求解變量，其具體實(shí)現(xiàn)方法如下：

1）不可壓縮流動(dòng)的湍流黏性系數(shù)：

2）可壓縮流動(dòng)的湍流黏性系數(shù)：

其中：

3）SA一方程湍流模型[37]：

其中常數(shù)cv1=7.1。

4）k-ω兩方程湍流模型[38-39]：

其中常數(shù)Cμ=0.09。

Knopp等[16]研究了在壁面處湍流變量的修正同湍流模型方程的相容性問題。不考慮壓力梯度影響時(shí)，湍流模型相容性條件意味著需要特定的壁面函數(shù)模型模型。不同版本的SA模型和k-ω模型的近壁剖面幾乎重疊，充分說明了壁面函數(shù)與對(duì)應(yīng)的湍流模型是相容的，其具體表達(dá)式如下：

1）SA湍流模型：

2）k-ω湍流模型：

圖6給出了SA和k-ω湍流模型壁面函數(shù)的曲線圖，Knopp等[16]采用Reichardt的壁面律：

與經(jīng)典對(duì)數(shù)律函數(shù)Flog=lny+/κ+B相融合，采用如下的混合函數(shù)：

以及帶參數(shù)的Spaldings[9]壁面律可以表示為：

為了避免在過渡層中ω偏離其低雷諾數(shù)的解，將標(biāo)準(zhǔn)混合函數(shù)（式（45））替換為：

其中：

上述壁面函數(shù)均基于臨近壁面計(jì)算單元上的速度（可壓縮情況下，還包括溫度）構(gòu)造，可以適用于基于一方程和二方程湍流模型的RANS模擬，正如方程(11)，基于湍動(dòng)能的速度壁面函數(shù)，比較方便用于包含湍動(dòng)能的湍流模型RANS模擬中，但以上壁面函數(shù)都只適用于不可壓縮流動(dòng)，而沒有考慮可壓縮性和壓力梯度的影響。

圖6 SA和k-ω模型壁面函數(shù)曲線圖[16]Fig. 6 Wall Functions for SA and k-ω models[16]

國內(nèi)Gao等[24]在RANS模擬中采用壁面函數(shù)時(shí)，著重研究了壁面函數(shù)與RANS數(shù)值解之間的一致性問題，計(jì)算結(jié)果表明：如果壁面函數(shù)與RANS數(shù)值解一致性好，則相應(yīng)的近壁面網(wǎng)格限制條件可由y+≤100放寬到y(tǒng)+≤400。

最近提出的“子網(wǎng)格壁面函數(shù)”[40-42]，將第一層網(wǎng)格進(jìn)一步細(xì)分成一個(gè)子網(wǎng)格，在第一層網(wǎng)格區(qū)域的子網(wǎng)格上數(shù)值求解邊界層方程，而不再追求獲得其解析解，文獻(xiàn)作者聲稱與同等密度的單一網(wǎng)格求解方法相比，計(jì)算效率可以提高一個(gè)量級(jí)。顯然，這種做法，與Tidriri[36]的區(qū)域分解思想一致，可以避免由于簡化得到的解析壁面函數(shù)與實(shí)際復(fù)雜情況下的邊界層流動(dòng)解一致性很差的問題，從發(fā)展趨勢(shì)上講，是值得進(jìn)一步研究的技術(shù)途徑。

5 結(jié) 論

圍繞發(fā)展自主可控CFD軟件的目標(biāo)，以建立高超聲速湍流邊界層工程實(shí)用模擬方法為目的，從湍流壁面函數(shù)是湍流邊界層方程近似解的角度，對(duì)相關(guān)公開文獻(xiàn)的研究工作進(jìn)行了梳理，得到以下認(rèn)識(shí)：

1）解析形式的壁面函數(shù)，在復(fù)雜工程問題求解中已經(jīng)得到了廣泛的應(yīng)用。但對(duì)于以高超聲速流動(dòng)為代表的具有強(qiáng)壓縮性、顯著傳熱和大壓力梯度的流動(dòng)，現(xiàn)有的表達(dá)形式同微分方程系統(tǒng)細(xì)密網(wǎng)格上解的一致性問題沒有得到很好的解決，可能會(huì)導(dǎo)致計(jì)算結(jié)果存在明顯的誤差，特別是壓力梯度比較大而壁面剪切應(yīng)力比較小的區(qū)域，如分離、再附點(diǎn)附近的流動(dòng)。

2）“子網(wǎng)格”壁面函數(shù)，其實(shí)質(zhì)是通過數(shù)值求解邊界層方程得到的數(shù)值形式的壁面函數(shù)，能夠比較充分地反映流動(dòng)可壓縮、傳熱和壓力梯度以及物面幾何特征等的影響，對(duì)復(fù)雜飛行器的壁面湍流流動(dòng)模擬應(yīng)該是一種能夠同時(shí)兼顧計(jì)算精度和效率的解決方案，但如何高效求解是“子網(wǎng)格”壁面函數(shù)方法在復(fù)雜問題中發(fā)揮作用的關(guān)鍵。

壁面函數(shù)的研究已有將近60年的歷史，本文僅僅是對(duì)RANS湍流邊界層壁面函數(shù)的部分研究成果的歸納總結(jié)。在目前的計(jì)算機(jī)資源條件下，系統(tǒng)研究“子網(wǎng)格”壁面函數(shù)，并與數(shù)據(jù)挖掘技術(shù)相結(jié)合，有可能發(fā)展出更加普適的解析壁面函數(shù)，以滿足對(duì)包含強(qiáng)可壓縮性、顯著傳熱和強(qiáng)壓梯度以及化學(xué)反應(yīng)等復(fù)雜物理流動(dòng)現(xiàn)象和復(fù)雜幾何構(gòu)型的壁面湍流流動(dòng)的高效高精度模擬需求。

另外，需要特別指出的是：近年來大渦模擬在科學(xué)研究和工程設(shè)計(jì)中應(yīng)用越來越廣泛，對(duì)于簡單機(jī)翼繞流算例，當(dāng)雷諾數(shù)Re= 1×109時(shí)，壁模型大渦模擬的網(wǎng)格量相較于近壁區(qū)分辨（wall-resolved）大渦模擬可以減少3個(gè)數(shù)量級(jí)[43]，也涌現(xiàn)了大量關(guān)于壁模型的研究，因此，壁模型的研究對(duì)于大渦模擬在工程應(yīng)用中的推廣具有積極意義，相關(guān)的綜述性文章可參考文獻(xiàn)[44-46]。我們也將會(huì)持續(xù)關(guān)注該方向的研究進(jìn)展。