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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

或者表示為：

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

則可壓縮流動的黏性底層速度分布可表述為：

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

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

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

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

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

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

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

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

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

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

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

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

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

其中：

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

圖3比較了統(tǒng)一壁面函數(shù)的不同表達(dá)式（23）、式（24）和式（25）的速度型，當(dāng)來流馬赫數(shù)為11.1時，在黏性底層內(nèi)的溫度、密度變化較大，整體考慮可壓縮效應(yīng)的統(tǒng)一壁面函數(shù)公式（25）更接近于數(shù)值解，而在來流馬赫數(shù)為5時三種表達(dá)形式差異不大。如果忽略黏性底層溫度變化對速度分布的影響，則公式（25）退化到不可壓縮絕熱邊界層的表達(dá)式。在不可壓縮絕熱流動時，式（25）也和式（24）一樣，能夠自動退化到式（23），同時還比較全面地考慮了可壓縮性和傳熱對黏性底層和對數(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時，黏性底層內(nèi)大體相同，僅略有差異，而當(dāng)Ma∞= 11.1時差異較大，整體來看，公式（19）和式（22）在黏性底層和對數(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]也提出了一個混合方法，將黏性底層與對數(shù)律層的分布函數(shù)整合為一個表達(dá)式，以便統(tǒng)一描述近壁區(qū)域（黏性底層、過渡層、完全湍流區(qū)）的物理量分布。

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

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

其中：

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

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

分裂為：

其中：

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

其中：

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

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

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

其逆函數(shù)記為：

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

其逆函數(shù)記為：

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

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

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

其中：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中：

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

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

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

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

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

1）SA湍流模型：

2）k-ω湍流模型：

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

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

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

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

其中：

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

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

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

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

5 結(jié) 論

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

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

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

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

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