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

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

0 引 言

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

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

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

對于管道和平板湍流邊界層的分析，大多數(shù)早期的壁面函數(shù)并不考慮壓力梯度和傳熱的影響。1938年Millikan[3]給出了剪切力占主導作用的不可壓縮較大Reynolds數(shù)流動的對數(shù)律層的漸近解，然而在流動分離點和再附點附近壁面剪切力為零，故Millikan的漸近解在分離點和再附點附近不再成立。Tennekes和Lumley[4]推導了考慮壓力梯度和零剪應力情況下邊界層的漸近解，但Tennekes和Lumley的解不能適用于壓力梯度為零的情況。Nichols[5]、Viegas[6]和Smith[7]等都相繼發(fā)展了包含壁面?zhèn)鳠嵝谋诿婧瘮?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的工作也給出了一個包含壓力梯度效應的統(tǒng)一壁面函數(shù),但仍然沒有將流動可壓縮性和壁面?zhèn)鳠嵝詈线M來。Nichols和Nelson[12]在White和Christoph[13]工作的基礎上，忽略壓力梯度的影響，將速度和溫度函數(shù)耦合起來，發(fā)展了適用于不可壓縮流動、可壓縮流動、絕熱壁和等溫壁湍流邊界層的壁面函數(shù)。由于壁面律在壓力梯度不可忽略的區(qū)域內是否仍然有效尚無定論[14-15]，Nichols和Nelson忽略壓力梯度影響的壁面函數(shù)對于分離點和再附點附近流動的模擬會存在困難。

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

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

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

1 標準壁面函數(shù)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

或者表示為：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中：

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

圖3比較了統(tǒng)一壁面函數(shù)的不同表達式（23）、式（24）和式（25）的速度型，當來流馬赫數(shù)為11.1時，在黏性底層內的溫度、密度變化較大，整體考慮可壓縮效應的統(tǒng)一壁面函數(shù)公式（25）更接近于數(shù)值解，而在來流馬赫數(shù)為5時三種表達形式差異不大。如果忽略黏性底層溫度變化對速度分布的影響，則公式（25）退化到不可壓縮絕熱邊界層的表達式。在不可壓縮絕熱流動時，式（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所示，當Ma∞= 5時，黏性底層內大體相同，僅略有差異，而當Ma∞= 11.1時差異較大，整體來看，公式（19）和式（22）在黏性底層和對數(shù)律層分別取層流和湍流普朗特數(shù)時，更接近于數(shù)值模擬結果。

圖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ù)整合為一個表達式，以便統(tǒng)一描述近壁區(qū)域（黏性底層、過渡層、完全湍流區(qū)）的物理量分布。

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

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

其中：

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

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

分裂為：

其中：

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

其中：

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

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

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

其逆函數(shù)記為：

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

其逆函數(shù)記為：

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

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

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

其中：

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

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

也可進一步改寫為：

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

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

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

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

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

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

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

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

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

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

1）SA湍流模型：

2）k-ω湍流模型：

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

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

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

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

其中：

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

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

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

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

5 結 論

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

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

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

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

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