0 引　言

尽管近几十年来计算机资源飞速增长，但是巨大的计算资源消耗仍然是限制直接数值模拟和大涡模拟在工程复杂湍流问题上应用的一个根本因素。即使是采用计算资源消耗最少的RANS模拟方法，要获得准确的壁面摩阻和热流，通常仍然需要在黏性底层中布置一定数目的网格点，从而使得壁面附近的网格十分细密。这样，不仅会极大地增加迭代收敛的计算步数，还会带来较为严重的数值刚性问题，导致迭代计算的稳定性下降。湍流边界层模拟采用壁面函数方法，可以大幅放宽壁面第一层网格的尺度，从而达到加速计算过程收敛的效果。另外，壁面函数方法不仅可以应用于RANS模拟，也可以用于大涡模拟^[1]。因此，目前主流商业软件在工程湍流模拟应用中默认采用壁面函数方法。

图1给出了湍流边界层沿壁面法向的典型速度分布，在紧挨壁面的黏性底层中，由于受到壁面的限制，速度脉动通常小到可以忽略，流动黏性中分子黏性占主导，湍流涡黏性可以忽略；在对数律层和外层，湍流脉动得到充分发展，流动黏性由湍流涡黏性主导，分子黏性效应可以忽略；在黏性底层和对数律层之间的过渡层，分子黏性和湍流黏性同等重要，两者共同影响速度分布。壁面函数基于如下的事实：对于不可压缩流动，从壁面到对数律层外缘之间在局部平衡的边界层中存在相似解。由于最初使用“壁面律（law of the wall）”仅指普适的对数律层解，而“混合壁面律（hybrid law of the wall）”或“自适应壁面律（adaptive law of the wall）”则主要强调在整个近壁区域内一致有效。对于不可压缩平板湍流边界层而言，在Prandtl混合长度理论和若干假设下，利用若干实验观察/测量结果，可以证明平均速度的壁面函数就是边界层方程（RANS方程在Prandtl边界层假设下的近似）的近似解^[2]。

对于管道和平板湍流边界层的分析，大多数早期的壁面函数并不考虑压力梯度和传热的影响。1938年Millikan^[3]给出了剪切力占主导作用的不可压缩较大Reynolds数流动的对数律层的渐近解，然而在流动分离点和再附点附近壁面剪切力为零，故Millikan的渐近解在分离点和再附点附近不再成立。Tennekes和Lumley^[4]推导了考虑压力梯度和零剪应力情况下边界层的渐近解，但Tennekes和Lumley的解不能适用于压力梯度为零的情况。Nichols^[5]、Viegas^[6]和Smith^[7]等都相继发展了包含壁面传热效应的壁面函数，但这些壁面函数都不是统一的速度分布函数，在黏性底层与对数律层之间的过渡区域也是不连续的，Wilcox^[8]引入一种壁面匹配方法。基于Spalding^[9]统一壁面函数理论，Frink^[10]发展了绝热壁的壁面函数，Shih等^[11]基于Spalding的工作也给出了一个包含压力梯度效应的统一壁面函数,但仍然没有将流动可压缩性和壁面传热效应耦合进来。Nichols和Nelson^[12]在White和Christoph^[13]工作的基础上，忽略压力梯度的影响，将速度和温度函数耦合起来，发展了适用于不可压缩流动、可压缩流动、绝热壁和等温壁湍流边界层的壁面函数。由于壁面律在压力梯度不可忽略的区域内是否仍然有效尚无定论^[14-15]，Nichols和Nelson忽略压力梯度影响的壁面函数对于分离点和再附点附近流动的模拟会存在困难。

实际上，考虑到壁面函数的使用能有效减少湍流边界层数值模拟结果对壁面网格尺度的关联性，RANS模拟中湍流模型在具有压力梯度的流动中表现非常重要。Knopp等^[16]指出：尽管壁面函数的研究工作已经持续开展了50多年，其中大多数研究成果的应用工作并不乐观，如壁面函数与用于全局流动问题模拟的湍流模型不相容的问题^[10,17-18]，容易引起数值模拟结果对壁面附近第一层网格节点位置的强烈依赖性，导致分离流动的数值模拟结果精度很差，因而提出了网格和流动自适应的壁面函数方法，保证在壁面函数有效的情况下对近壁流动物理特性恰当的分辨，特别是对非平衡效应显著的流动驻点区域、逆压梯度导致的流动分离/再附点附近区域，以获得比较准确的气动力/热特性。此外，Kalitzin等^[19]也针对壁面函数在RANS模拟中的应用问题开展了相关研究工作。

从公开的文献来看，关于湍流边界层壁面函数的研究，原创性工作几乎都来自国外，国内相关工作较少，比较典型的研究工作有：贺旭照等^[20-21]将Nichols的考虑流动可压缩性和热传导的壁面函数耦合到了k-ω两方程模型中；吴晓军^[22]和肖红林^[23]分别在RANS模拟和大涡模拟中采用壁面函数；Gao^[24]和Tao等^[25]分别通过数值实验和风洞实验对壁面函数的系数进行了修正研究；Wu等^[26]对Nichols的壁面函数进行了修正研究。目前可压缩湍流特别是强压缩和显著传热的湍流理论研究，几乎都是对不可压湍流研究成果的修正、延拓和推广。国内开展高超声速湍流边界层壁面函数的研究工作相对较少，缺乏系统性，也缺乏原创性。

本文的目的是对现有研究工作进行梳理，从理论上比较系统地把握湍流边界层壁面函数理论，为研制自主可控的高超声速工程实用的CFD软件提供理论支撑。

1 标准壁面函数

二维可压缩定常平板湍流边界层问题，通过量级分析，忽略压力在边界层内沿壁面法向的变化，可以得到湍流边界层内的流动控制方程为：

$\frac{{\partial \rho u}}{{\partial x}} + \frac{{\partial \rho v}}{{\partial y}} = 0$

(1)

$\frac{\partial }{{\partial y}}\left( {{\tau _{xy}} - \rho \widetilde {u''v''}} \right) = 0$

(2)

$\frac{\partial }{{\partial y}}\left( {{\tau _{xy}}u + k\frac{{\partial T}}{{\partial y}}} \right) = 0$

(3)

其中，x为来流方向，y为壁面法向。由x方向的动量方程（2）可知总剪切力在边界层内沿壁面法向保持不变。由于在黏性底层速度脉动被壁面抑制，湍流涡黏性可以忽略不计，记：

${\tau _w} \triangleq {\left( {{\tau _{xy}} - \rho \widetilde {u''v''}} \right)_{y = 0}}{\rm{ = }}{\left. {\mu \frac{{\partial u}}{{\partial y}}} \right|_{y = 0}}$

(4)

在远离壁面一定距离后，湍流脉动充分发展，湍流涡黏性远大于分子黏性，忽略分子黏性的影响可以得到：

$ - \rho \widetilde {u''v''} = {\tau _w}$

(5)

对于不可压流动，密度ρ为常数，对壁面剪切力定义式（4）在壁面附近积分，并考虑到速度在壁面处的无滑移边界条件可以得到：

${u^ + } = {y^ + },\;0 < {y^ + } < y_L^ + $

(6)

其中， ${u^ + } = u/{u_\tau }$ 为无量纲速度， ${y^ + } = \rho y{u_\tau }/\mu $ 为无量纲壁面距离， ${u_\tau } = \sqrt {{\tau _w}/\rho } $ 为摩擦速度。 $y_L^ + $ 为黏性底层外缘处的无量纲常数，主要由实验观察结果确定，不同的文献略有差异，大致取值范围为 $y_L^ + \in \left( {5,15} \right)$ 。

在湍流脉动充分发展的对数律层内，依据Prandtl混合长度理论给出：

$ - \rho \widetilde {u''v''} = \rho {l_m}{l_m}\left| {\frac{{\partial u}}{{\partial y}}} \right|\frac{{\partial u}}{{\partial y}}$

(7)

其中， ${l_m} = \kappa y$ 为对数律层内的Prandtl混合长度， $\kappa \approx $ $ 0.41$ ，为von Kármán常数。考虑到湍流脉动充分发展时可忽略分子黏性的影响（式（5）），容易得到：

$\frac{{\partial u}}{{\partial y}} = \frac{{{u_\tau }}}{{\kappa y}}\;$

(8)

对式（8）进行积分可以得到对数律表达式为：

${u^ + } = (1/\kappa )\ln {y^ + } + B,\;30 < {y^ + } < 1\;000$

(9)

其中常数 $B \in \left[ {5,5.5} \right]$ 。

由上述分析可知，壁面函数实质上就是边界层方程的近似解。需要特别指出的是： ${\tau _w}$ 不仅是壁面处的分子黏性应力，也等于边界层剖面中的总黏性应力（ ${\tau _w} = {\tau _{xy}} - \rho \widetilde {u''v''}$ ，即湍流涡黏性应力与分子黏性应力之和），故 ${\tau _w}$ 既可以通过壁面的分子黏性应力公式（4）计算得到，也可以通过对数律层中的湍流涡黏性应力计算得到。

在靠近壁面的剪切流区域，其流动特性与Couette流动相似，在湍动能生成与耗散处于局部平衡的条件下，壁面剪切应力 ${\tau _w}$ 与湍动能k之间满足关系：

$\frac{{{\tau _w}}}{\rho } = \sqrt {{C_\mu }} k$

(10)

其中 ${C_\mu } = 0.09$ 为常数，于是平均速度剖面可以改写为：

$\begin{split}& u^* = y^*,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 < y^* < y_L^* \\& u^* = \frac{1}{\kappa }\ln y^* + B,\;\;\;\;30 < y^* < 1\;000\; \end{split} $

(11)

其中， $u^* = {{uC_\mu ^{1/4}k_{}^{1/2}} / {\left( {{\tau _w}/\rho } \right)}}$ , $y^* = {{\rho yC_\mu ^{1/4}k_{}^{1/2}} / \mu }$ 。

对于可压缩流动，特别是高超声速流动，温度和密度在边界层内存在明显的变化，导致分子黏性系数在边界层内的变化不可忽略。Fernholz和Finley^[27]指出，若平均速度采用变换：

${u_{c,T}} = \int {{{\left( {{\rho / {{\rho _w}}}} \right)}^{1/2}}{\rm{d}}u} $

(12)

则可压缩流动的对数律可以表述为^[28]：

$u_{c,T}^ + \buildrel \Delta \over = \frac{{{u_{c,T}}}}{{{u_\tau }}} = \frac{1}{\kappa }\ln {y^ + } + B$

(13)

或者表示为：

${y^ + } = {{\rm{e}}^{ - \kappa B}}{{\rm{e}}^{\kappa u_{c,T}^ + }}$

(14)

对于黏性底层，有不少工作忽略可压缩性的影响，仍然采用不可压缩壁面函数的表达式。实际上，对于高超声速流动，黏性底层中密度变化可能是不可忽略的，图2给出了不同马赫数情况下平板边界层^[29-31]密度分布，可以看出，随着马赫数的增加，黏性底层内的密度变化愈加明显。当来流马赫数为11.1时，黏性底层顶部附近，密度变化接近40%。根据边界层方程，类似地引入变换：

${u_{c,L}} = \int_0^u {\left( {{\mu / {{\mu _w}}}} \right){\rm{d}}u} $

(15)

则可压缩流动的黏性底层速度分布可表述为：

$u_{c,L}^ + \buildrel \Delta \over = \frac{{{u_{c,L}}}}{{{u_\tau }}} = {y^ + }$

(16)

对于空气而言， ${\mu / {{\mu _w}}}$ 是 ${T/ {{T_w}}}$ 的函数，在边界层内压力沿壁面法向保持不变，由状态方程可知 ${\rho / {{\rho _w}}}$ 也仅是 ${T / {{T_w}}}$ 相关，只要获得了温度在边界层内的表达式，就可以唯一确定可压缩流动条件下的平均速度壁面函数。

温度在边界层内的分布可由边界层能量方程得到，简略推导如下：将边界层能量方程（3）沿边界层法向积分，得到：

${\tau _w}u + k\frac{{\partial T}}{{\partial y}}{\rm{ = }}{\left( {{\tau _w}u + k\frac{{\partial T}}{{\partial y}}} \right)_{y = 0}}{\rm{ = - }}{q_w}$

(17)

其中用到了总剪切应力在湍流边界层内层沿壁面法向不变的条件， ${q_w} = - {\left( {k{{\partial T} / {\partial y}}} \right)_{y = 0}}$ 为壁面热流，k为传热系数。在黏性底层，分子黏性起主要作用，而在对数律层，则湍流黏性起主要作用，故有：

$\begin{split}& k \approx \frac{{{\mu _t}{C_p}}}{{P{r_t}}},\;\;\;\;30 < {y^ + } < 1\;000\; \\& k \approx \frac{{{\mu _L}{C_p}}}{{P{r_L}}},\;\;\;\;0 < {y^ + } < y_L^ + \end{split} $

(18)

在黏性底层（ $0 < {y^ + } < y_L^ + $ ）考虑传热系数（18）和速度分布（6），由能量积分方程（17）可得黏性底层（ $0 < {y^ + } < y_L^ + $ ）中的温度分布为：

$\frac{T}{{{T_w}}} = 1 - P{r_L}\frac{{{T_\tau }}}{{{T_w}}}{u^ + } - P{r_L}\frac{{\gamma - 1}}{2}M_\tau ^2{ {{u^ + }} ^2}$

(19)

其中， ${T_\tau } = {q_w}/\left({\rho _w}{u_\tau }{C_p}\right)$ ， ${M_\tau } = {u_\tau }/\sqrt {\gamma R{T_w}} $ 。

在对数律层湍流脉动充分发展，湍流涡黏性起主要作用，可忽略分子黏性的影响即：

${\tau _w} = - \rho \widetilde {u''v''} = {\mu _T}\frac{{\partial u}}{{\partial y}}$

(20)

考虑传热系数式（18），由能量积分方程（17）可得对数律层的温度关系为：

$ - \frac{{{C_p}{\tau _w}}}{{P{r_T}}}\frac{{\partial T}}{{\partial y}} = {q_w}\frac{{\partial u}}{{\partial y}} + u{\tau _w}\frac{{\partial u}}{{\partial y}}$

(21)

从对数律层下缘开始对式（21）积分得到：

$\frac{T}{{{T_w}}} = {c_1} - P{r_t}\frac{{{T_\tau }}}{{{T_w}}}{u^ + } - P{r_t}\frac{{\gamma - 1}}{2}M_\tau ^2{ {{u^ + }} ^2}$

(22)

其中 ${c_1} = \dfrac{{{T_1}}}{{{T_w}}} + P{r_t}\dfrac{{{T_\tau }}}{{{T_w}}}u_1^ + + P{r_t}\dfrac{{\gamma - 1}}{2}M_\tau ^2{ {u_1^ + } ^2}$ 。

湍流边界层沿壁面法向的剖面包括黏性底层、过渡层、对数律层以及边界层外层等多层结构，对应地，湍流边界层的壁面函数通常描述从壁面到完全湍流区域的剖面。因此，严格来说应该采用三层模型，但有不少研究者忽略对过渡层的模拟，而采用两层模型。对于速度壁面函数，当 ${y^ + } < 11.225$ 时，采用黏性底层公式（6）计算，否则，采用对数律层公式（9）描述。而对于温度壁面函数，由于其分布还同壁面温度等相关。实际上式（19）和式（22）共同构成了一个两层的温度壁面函数，在交界点 ${y^ + } = y_{{j}}^ +$ 处，采用黏性底层公式（19）和对数律层公式（22）计算得到的无量纲温度应相同，则有：

${c_1} = 1 + \left( {P{r_t} - P{r_L}} \right)\left( {\frac{{{T_\tau }}}{{{T_w}}}u_j^ + + \frac{{\gamma - 1}}{2}M_\tau ^2{{ {u_{j}^ + } }^2}} \right)$

Huang等^[32]假定黏性底层的等效Prandtl数与湍流Prandtl数相同，则式（19）和式（22）相同。从已有的数值试验结果来看，对速度分布进行可压缩修正时，这种假定对壁面摩擦应力造成的误差可以忽略，且简化了计算，但对等温壁壁面热流的影响如何，则需要进一步研究。

2 统一壁面函数

即使采用两层模型，在实际使用过程中也不太方便，且易造成平均速度分布不连续。Splading^[9]提出了在黏性底层和对数律层内一致有效的绝热不可压流动的湍流边界层壁面函数：

${y^ + } = {u^ + } + {{\rm{e}}^{ - \kappa B}}\left[ {{{\rm{e}}^{\kappa {u^ + }}} - \sum\limits_{n = 0}^3 {\frac{{{{\left( {\kappa {u^ + }} \right)}^n}}}{{n!}}} } \right]$

(23)

Nichols^[12]在考虑传热和可压缩性修正时，采用White^[13]的对数层壁面率公式（14）来替换Splading公式（23）中的绝热不可压壁面率 ${{\rm{e}}^{ - \kappa B}}{{\rm{e}}^{\kappa {u^ + }}}$ ，得到：

${y^ + } = {u^ + } + y_{{\rm{white}}}^ + - {{\rm{e}}^{ - \kappa B}}\sum\limits_{n = 0}^3 {\frac{{{{\left( {\kappa {u^ + }} \right)}^n}}}{{n!}}} $

(24)

其中：

$ y_{{\rm{white}}}^ + = \exp \left\{ { {\frac{\kappa }{{\sqrt {{\rm{ }}\varGamma } }}} \left[ {{{\sin }^{ - 1}}\left( {\frac{{ {2\varGamma {u^ + } - \beta } }}{Q}} \right) - \phi } \right]} \right\} \times {{\rm{e}}^{ - \kappa B}} $

$ \varGamma \!=\! \frac{{ru_\tau ^2}}{{2{C_p}{T_w}}},\;\beta \!=\! - \frac{{{q_w}{\mu _w}}}{{{\rho _w}{T_w}{k_w}{u_\tau }}},\;\phi \!=\! {\sin ^{ - 1}}\left( {\frac{\beta }{Q}} \right),\;Q \!=\! \sqrt {{\beta ^2} + 4\varGamma } $

一个比较自然的想法，将式（14）和式（16）类似于式（23）进行合并，可以得到：

${y^ + } = u_{c,L}^ + + y_{{\rm{white}}}^ + - {{\rm{e}}^{ - \kappa B}}\sum\limits_{n = 0}^3 {\frac{{{{\left( {\kappa u_{c,T}^ + } \right)}^n}}}{{n!}}} $

(25)

图3比较了统一壁面函数的不同表达式（23）、式（24）和式（25）的速度型，当来流马赫数为11.1时，在黏性底层内的温度、密度变化较大，整体考虑可压缩效应的统一壁面函数公式（25）更接近于数值解，而在来流马赫数为5时三种表达形式差异不大。如果忽略黏性底层温度变化对速度分布的影响，则公式（25）退化到不可压缩绝热边界层的表达式。在不可压缩绝热流动时，式（25）也和式（24）一样，能够自动退化到式（23），同时还比较全面地考虑了可压缩性和传热对黏性底层和对数律层的影响。

在可压缩边界层中，温度分布函数一般采用Crocco-Busemann公式：

$T = {T_w}(1 - \beta {u^ + } - \varGamma {{u^ + }^2})$

(26)

式（19）、式（22）和式（26）的比较如图4所示，当Ma_∞ = 5时，黏性底层内大体相同，仅略有差异，而当Ma_∞ = 11.1时差异较大，整体来看，公式（19）和式（22）在黏性底层和对数律层分别取层流和湍流普朗特数时，更接近于数值模拟结果。

此外，Kader^[33]也提出了一个混合方法，将黏性底层与对数律层的分布函数整合为一个表达式，以便统一描述近壁区域（黏性底层、过渡层、完全湍流区）的物理量分布。

3 增强型壁面函数

上述壁面函数的推导，均未涉及压力梯度。实际飞行器的边界层流动，由于壁面通常是曲面，并可能出现凸起和凹坑，从而使得流动具有压力梯度。众所周知，流场中的逆压梯度和顺压梯度对壁面流动的影响是非常重要的。逆压梯度不仅会使流动减速，甚至可能造成流动从壁面分离，使得壁面附近的流动失去了经典意义下的边界层特性；而顺压梯度则会加速边界层中的流动，也可能会压制湍流脉动，使湍流边界层流动层流化。增强型壁面函数设计的一个重要考量就是压力梯度的影响，使之能够适应流动分离和再附点附近壁面流动的模拟。对于逆压梯度和零剪应力边界层，Tennekes^[4]推导了一个渐近解：

$\frac{u}{{{u_p}}} = \alpha \ln \left( {\frac{{{u_p}y}}{\nu}} \right) + \beta $

(27)

其中：

${u_p} = {\left( {\frac{\nu }{\rho }\left| {\frac{{{\rm{d}}{p_w}}}{{{\rm{d}}x}}} \right|} \right)^{\frac{1}{3}}}$

常数 $\alpha \approx 5.0,\;\beta \approx 8$ 由Stratford^[34]的实验数据确定。显然，Tennekes的渐近解(式27)在压力梯度趋于零时不再成立，因为此时速度 ${u_p}$ 也趋于零。而在流动分离和再附点附近，剪切应力和摩擦速度趋于零，基于剪切应力特性的壁面函数（式（7））亦不再成立。

Shih等^[11]基于Tennekes的分解方法，将不可压缩流动边界层方程的积分形式：

$ - \overline {u'v'} + \nu \frac{{\partial u}}{{\partial y}} = \frac{{{\tau _w}}}{\rho } + \frac{y}{\rho }\frac{{\partial {p_w}}}{{\partial x}}$

(28)

分裂为：

$ - {\left( {\overline {u'v'} } \right)_1} + \nu \frac{{\partial {u_1}}}{{\partial y}} = \frac{{{\tau _w}}}{\rho }$

(29)

$ - {\left( {\overline {u'v'} } \right)_2} + \nu \frac{{\partial {u_2}}}{{\partial y}} = \frac{y}{\rho }\frac{{\partial {p_w}}}{{\partial x}}$

(30)

其中：

$ u = {u_1} + {u_2},\; - \overline {uv} = {\left( { - \overline {uv} } \right)_1} + {\left( { - \overline {uv} } \right)_2} $

利用量纲分析方法，容易得到渐近解：

$\begin{split} \frac{u}{{{u_{\tau + p}}}} = &\frac{{{\tau _w}}}{{\rho {u_\tau }{u_\tau }}}\frac{{{u_\tau }}}{{{u_{\tau + p}}}}\left[ {\frac{1}{\kappa }\ln \left( {\frac{{{u_\tau }y}}{\nu }} \right) + B} \right] +\\& \frac{\nu }{\rho }\frac{{\partial {p_w}}}{{u_p^3\partial x}}\frac{{{u_p}}}{{{u_{\tau + p}}}}\left[ {\alpha \ln \left( {\frac{{{u_p}y}}{\nu }} \right) + \beta } \right] \end{split}$

(31)

其中：

${u_{\tau + p}} = {u_\tau } + {u_p} = \sqrt {\frac{{\left| {{\tau _w}} \right|}}{\rho }} + {\left( {\frac{\nu }{\rho }\left| {\frac{{{\rm{d}}{p_w}}}{{{\rm{d}}x}}} \right|} \right)^{\frac{1}{3}}}$

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

$\begin{split} u =& \frac{{{\tau _w}}}{\rho }\frac{y}{\nu } + \frac{1}{2}\frac{\nu }{\rho }\frac{{{\rm{d}}{p_w}}}{{{\rm{d}}x}}{\left( {\frac{y}{\nu }} \right)^2} \\ {\tau _w} =& {\left. {\mu \frac{{\partial u}}{{\partial y}}} \right|_{y = 0}} \end{split} $

(32)

由此可知，边界层流动由壁面剪切应力和壁面压力梯度决定。不再会出现 ${u_{\tau + p}}$ = 0的情况，因为 ${u_\tau }$ 和 ${u_p}$ 不会同时为零。对于压力梯度起主要作用的流动，沿边界层法向黏性剪应力不再是常数；对于零压力梯度的情况，流动由壁面剪应力确定，边界层中黏性应力沿壁面法向为常数，且等同于壁面剪切应力。

式(31)和式(32)分别为对数律层和黏性底层的渐近解，两层之间的区域（ $5 \leqslant {u_{\tau + p}}y/\nu \leqslant 30$ ）称为过渡层（buffer layer），在过渡层中，黏性应力和湍流应力处于同一量级，尚未有文献在理论上给出过渡层的渐近解。在过渡层中采用简单的湍流应力模型和基于量级分析的匹配过程，能够得到唯一的解析函数，即适用于整个湍流边界层（包括黏性底层、过渡层、对数律层）的统一壁面函数。Spalding^[9]首先构造了一个统一壁面函数（23）（没有考虑压力梯度的影响）。为以示区别，将Spalding的统一壁面函数（23）记为：

$Y_\tau ^ + \approx u_1^ + + \exp \left( { - \kappa B} \right)\left\{ {\exp \left( {\kappa u_1^ + } \right) - \left[ {1 + \sum\limits_{n = 1}^3 {\frac{{{{\left( {\kappa u_1^ + } \right)}^n}}}{{n!}}} } \right]} \right\}$

(33)

其逆函数记为：

$\frac{{{u_1}}}{{{u_\tau }}} = \frac{{{\tau _w}}}{{\rho u_\tau ^2}}{f_1}\left( {Y_\tau ^ + } \right),Y_\tau ^ + = \frac{{{u_\tau }y}}{\nu }$

(34)

当压力梯度占主导作用时的壁面函数为：

${\left( {Y_p^ + } \right)^2} \approx u_2^ + + \exp \left( { - \frac{\beta }{\alpha} } \right)\left\{ {\exp \left( {\frac{u_2^ +} {\alpha }} \right) - \left[ {1 + \left( {\frac{{u_2^ + }}{\alpha }} \right)} \right]} \right\}$

(35)

其逆函数记为：

$\dfrac{{{u_2}}}{{{u_p}}} = \dfrac{{\left| {\dfrac{{d{p_w}}}{{dx}}} \right|}}{{\dfrac{{d{p_w}}}{{dx}}}}{f_2}\left( {Y_p^ + } \right),\;\;Y_p^ + = \frac{{{u_p}y}}{\nu }$

(36)

结合函数（34）和式（36）可以得到满足零壁面摩擦应力和零压力梯度的渐近解，即能够应用于流动加速、减速、分离和回流区域的速度壁面函数为：

$\dfrac{u}{{{u_{\tau ,p}}}} = \frac{{{u_\tau }}}{{{u_{\tau ,p}}}}\dfrac{{{\tau _w}}}{{\rho u_\tau ^2}}{f_1}\left( {\dfrac{{{u_\tau }}}{{{u_{\tau ,p}}}}{y^ + }} \right) + \dfrac{{{u_p}}}{{{u_{\tau ,p}}}}\dfrac{{\left| {\dfrac{{d{p_w}}}{{dx}}} \right|}}{{\dfrac{{d{p_w}}}{{dx}}}}{f_2}\left( {\frac{{{u_p}}}{{{u_{\tau ,p}}}}{y^ + }} \right)$

(37)

其中 ${y^ + } = {u_{\tau ,p}}y/\nu $ 。在黏性底层和过渡层函数 ${f_1}$ 和 ${f_2}$ 分别为关于 $Y_\tau ^ + $ 和 $Y_p^ + $ 的分片多项式函数，而在对数律层均为对数函数，其具体表达式由Shih等^[11]给出。

对于可压缩流动，White^[35]给出的对数律层平均速度分布函数为：

$\begin{split} {u^ + } = &\frac{1}{{2\gamma }}\left\{ {\beta + Q\sin \left[ {\phi + 2\frac{{\sqrt \gamma }}{\kappa }\left( {S - {S_0}} \right) + } \right.} \right.\\& \left. {\left. {\frac{{\sqrt \gamma }}{\kappa }\ln \left( {\frac{{S - 1}}{{S + 1}}\frac{{{S_0} + 1}}{{{S_0} - 1}}} \right)} \right]} \right\} \end{split}$

(38)

其中：

$ \begin{split}& \;{y^ + }\! \!=\!\! \frac{{{\rho _w}y{u_\tau }}}{{{\mu _w}}},\;\phi \!\!=\!\! {\sin ^{ - 1}}\left( {\frac{{ - \beta }}{Q}} \right),\;Q \!\!=\!\! \sqrt {{\beta ^2} + 4\gamma } ,\;S \!\!=\!\! \sqrt {1 + \alpha {y^ + }} ,\; \\& {S_0} = \sqrt {1 + \alpha y_0^ + } ,\;\alpha = \frac{{{\nu _w}}}{{{\tau _w}{u_\tau }}}{\left( {\frac{{{\rm{d}}p}}{{{\rm{d}}x}}} \right)_w},\;\beta = \frac{{{q_w}{\mu _w}}}{{{\rho _w}{T_w}{k_w}{u_\tau }}}{\rm{, }} \\& \gamma = P{r^{\frac{1}{3}}}\frac{{{u_\tau }{u_\tau }}}{{2{C_p}{T_w}}},\;y_0^ + = {{\rm{e}}^{ - \kappa B}} \approx 0.110\;8\;{\rm{.}} \end{split} $

需要特别指出的是：由于White等在推导上式时采用的混合长表达式 ${l_m} = \kappa y$ ，不能退化到绝热不可压对数律（31），在使用过程中一定要牢记这一点。如果压力梯度为零，则式（38）可退化为：

${u^ + } = \frac{1}{{2\gamma }}\left\{ {\beta + Q\sin \left[ {\phi + \frac{{\sqrt \gamma }}{\kappa }\ln \left( {\frac{{{y^ + }}}{{y_0^ + }}} \right)} \right]} \right\}$

(39)

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

${u_{c,L}} = \frac{{{\tau _w}}}{{{\mu _w}}}y + \frac{{{y^2}}}{{2{\mu _w}}}\frac{{\partial {p_w}}}{{\partial x}}$

(40)

也可进一步改写为：

$u_{c,L}^ + = {y^ + }\left( {1 + \frac{1}{2}\alpha {y^ + }} \right)$

(41)

在传热、可压缩性和压力梯度的共同作用下，黏性底层的高度与绝热不可压流动之间的差异，可能是一个需要重视的问题。

4 壁面函数的应用

RANS方法在模拟湍流流动时，采用壁面函数的初衷是在固壁处提供一个边界条件，降低数值解（特别是摩擦阻力和热流等）对临近壁面第一层网格节点位置的依赖性。以 $y_1^ + $ 表示第一个网格节点到固壁面的距离（采用黏性长度尺度 ${{{u_{\tau + p}}} / \nu }$ 进行无量纲化），当网格雷诺数较低时（ $y_1^ + \approx 1$ ），RANS方程可以直接积分到壁面，在壁面处的边界条件可以采用速度无滑移和相应的温度边界条件；而当网格雷诺数较高时（ $y_1^ + > $ $ 30$ ），此时第一层网格节点位于湍流充分发展的对数律层，只能在对数律层内求解RANS方程，在壁面处必须采用壁面函数才能求得比较准确的壁面剪切应力、热流和温度（具体由壁面温度边界条件确定）等物理量；对于中等网格雷诺数问题，第一层网格节点位于黏性底层和对数律层之间的过渡层，壁面边界的处理与高网格雷诺数情况类似。

在实际计算中，网格生成是无法预先确定壁面第一层网格中每个网格点位于湍流边界层的具体位置，只有在计算过程中，根据其速度、温度和壁面边界条件，通过壁面函数才能确定。对于不可压缩流动，边界层的具体位置可由第一层网格对应的网格雷诺数唯一确定。因为 ${u^ + }{y^ + } = R{e_y} = {{\rho u{y_1}} / \mu }$ 恒成立，因此，在黏性底层中 ${y^ + } = \sqrt {R{e_y}} $ 恒成立。据此可以依据 $R{e_y}$ 的大小来判断临近壁面第一层网格节点处于边界层的哪个子层中，同时壁面摩擦速度、剪切应力等参数均可由 $R{e_y}$ 唯一确定。基于此认识，有文献给出了不可压缩流动的壁面摩擦速度、剪切应力等参数与 $R{e_y}$ 的拟合函数，在具体计算中，无需对非线性方程组进行迭代求解来得到壁面摩擦应力等相关参数，以进一步节省计算时间，提高计算效率。对于压力梯度影响不可忽略的不可压缩流动，原则上，也可以构造双参数的拟合函数，但目前尚未看到相关的研究工作。对于可压缩流动，边界层内速度、温度的分布所依赖的参数，除了第一层网格对应的网格雷诺数和压力梯度外，还有壁面温度边界条件参数、来流马赫数和温度等，由于相关参数较多，构造拟合函数的工作量和难度会急剧增大。

湍流边界层壁面函数的理论和实验研究工作，主要基于平板边界层，特别是理论研究。实际飞行器的表面，普遍是三维曲面。将二维平板湍流边界层的理论近似结果，推广到三维曲面边界层流动，其可靠性需通过数值模拟试验同风洞和实际飞行试验结果的一致性来评估。

从平面推广到曲面，一个自然的作法，就是把边界层方程及其相关理论表达式转换到贴体坐标系中。以当地切向速度作为流动方向，将三维流动进行局部二维化处理，在此基础上采用湍流边界层的壁面函数，可以得到壁面剪切应力、热流（或壁面温度）等。对于压力梯度影响可以忽略的流动，这种做法似乎是一种比较自然的推广，但对于压力梯度与流动方向存在明显差异的流动情况，如何推广应用已有的壁面函数研究成果还需要进一步研究。

壁面函数在CFD计算过程中的具体实施方法，Tidriri^[36]从计算区域分解和叠加求解的角度，给出了一个十分清晰的阐述，如图5所示。

记 ${{\varOmega }_{\rm{\delta }}} \subset {\varOmega }$ 为计算区域 ${\varOmega }$ 的近壁区域，内边界位于对数层内，那么壁面函数在CFD中的应用问题可以分解为两个问题：1）全局流动计算问题，在整个计算区域 ${\varOmega }$ 内求解RANS方程，利用壁面函数改进边界条件（在固壁面强加剪应力条件，而不是速度无滑移条件，等温固壁面强加热流条件而不仅仅给定壁面温度），以弥补由于网格无法分辨湍流边界层中的近壁流动结构产生的误差；2）近壁区域 ${{\varOmega }_{\rm{\delta }}}$ 内湍流边界层的求解问题。显然，解的具体形式就是湍流边界层方程的近似解即壁面函数，在固壁面 ${{\varGamma }_{w}}$ 上始终自动满足原始边界条件（速度无滑移条件、等温或绝热壁面条件），在近壁区域 ${{\varOmega }_{\rm{\delta }}}$ 的外缘 ${{\varGamma }_{\rm{\delta }}}$ 上与全局的RANS数值解匹配一致。该方法对 ${{\varGamma }_{\rm{\delta }}}$ 处于任意位置都是适定的，所得的解与 ${{\varGamma }_{\rm{\delta }}}$ 的具体位置无关（只要 ${{\varGamma }_{\rm{\delta }}}$ 位于湍流边界层的对数律层以内的近壁流动区域中），这意味着壁面压力、摩阻和热流等物理量与网格分布密度无关。目前实际计算中几乎都设定 ${{\varGamma }_{\rm{\delta }}}$ 为离开壁面的第一层网格形成的曲面。由此可见，在RANS模拟中采用壁面函数方法，其实质是避免在近壁区域数值求解完全可压缩RANS方程和湍流模型方程。

由上所述易知，不仅求解RANS方程需要用到壁面函数边界条件，求解湍流模型方程也同样需要。在实际工程应用的航空航天飞行器湍流流动的模拟中，最常用的湍流模型为以SA模型为代表的一方程模型和以SST $k - \omega $ 模型为代表的两方程模型等。Nichols和Nelson^[12]在忽略压力梯度的情况下，基于剪切应力沿边界层法向基本不变的假设，利用壁面函数公式、临近壁面第一层网格单元上的速度、温度（或绝热边界条件）和壁面距离等物理流场和几何信息，采用迭代方法获得壁面摩擦应力和热流（或壁面温度），并以此为基础，修正临近壁面网格单元上的湍流涡黏性系数和湍流模型方程中的求解变量，其具体实现方法如下：

1）不可压缩流动的湍流黏性系数：

$\frac{{{\mu _t}}}{{{\mu _w}}} = \kappa {{\rm{e}}^{ - \kappa B}}\left[ {{{\rm{e}}^{\kappa u_1^ + }} - \sum\limits_{n = 0}^2 {\frac{{{{\left( {\kappa u_1^ + } \right)}^n}}}{{n!}}} } \right]$

(42)

2）可压缩流动的湍流黏性系数：

$\frac{{{\mu _t}}}{{{\mu _w}}} = 1 + \frac{{\partial y_{{\rm{white}}}^{\rm{ + }}}}{{\partial y}} - \kappa {{\rm{e}}^{ - \kappa B}}\sum\limits_{n = 0}^2 {\frac{{{{\left( {\kappa u_1^ + } \right)}^n}}}{{n!}}} - \frac{{{\mu _{w + 1}}}}{{{\mu _w}}}$

(43)

其中：

$\frac{{\partial y_{{\rm{white}}}^{\rm{ + }}}}{{\partial y}} = 2y_{{\rm{white}}}^{\rm{ + }}\frac{{\kappa \sqrt \varGamma }}{Q}{\left[ {1 - {{\left( {\frac{{2\varGamma {u^ + } - \beta }}{Q}} \right)}^2}} \right]^{ - 1/2}}$

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

${ {{{\tilde \mu }_T}} ^4} = {\mu _T}\left[ {{{ {{{\tilde \mu }_T}} }^3} + {{\left( {{c_{v1}}{\mu _L}} \right)}^3}} \right]$

(44)

其中常数 ${c_{v1}} = 7.1$ 。

4） $k - \omega $ 两方程湍流模型^[38-39]：

$\begin{split} {\mu _t} =& \frac{{\rho k}}{\omega }\;,\;\omega = \sqrt {{\omega _i}{\omega _i} + {\omega _{\rm{o}}}{\omega _{\rm{o}}}}\;, \\ {\omega _i} = &\frac{{80{\mu _w}}}{{{\rho _w}{y^2}}}\;,\;{\omega _{\rm{o}}} = \frac{{{u_\tau }}}{{\sqrt {{C_\mu }} ky}} \end{split} $

(45)

其中常数 ${C_\mu }=0.09$ 。

Knopp等^[16]研究了在壁面处湍流变量的修正同湍流模型方程的相容性问题。不考虑压力梯度影响时，湍流模型相容性条件 $\nu _t^{{\rm{bl}}} = {\nu _t}$ 意味着需要特定的壁面函数模型模型。不同版本的SA模型和 $k - \omega $ 模型的近壁剖面几乎重叠，充分说明了壁面函数与对应的湍流模型是相容的，其具体表达式如下：

1）SA湍流模型：

$\begin{split} {F_{{\rm{SA}}}} =& \left( {1 - {\phi _{{\rm{SA}}}}} \right){F_{{\rm{Sp,5}}}} + {\phi _{{\rm{SA}}}}{F_{{\rm{Rei,}}m}} \\ {\phi _{{\rm{SA}}}} =& \tanh \left( {{{\arg }^3}} \right),\;\arg = {y^ + }/24 \end{split}$

(46)

2） $k - \omega $ 湍流模型：

$\begin{split} {F_{k\omega }} = &\left( {1 - {\phi _{k\omega }}} \right){F_{{\rm{Sp,3}}}} + {\phi _{k\omega }}{F_{{\rm{Rei}},m}} \\ {\phi _{k\omega }} =& \tanh \left( {{{\arg }^2}} \right),\;\arg = {y^ + }/50 \end{split} $

(47)

图6给出了SA和 $k - \omega $ 湍流模型壁面函数的曲线图，Knopp等^[16]采用Reichardt的壁面律：

$ {u^ + } = {F_{{\rm{Rei}}}}\left( {{y^ + }} \right) = \dfrac{{\ln \left( {1 + 0.4{y^ + }} \right)}}{\kappa } + 7.8\left( {1 - {{\rm{e}}^{ - \tfrac{{{y^ + }}}{{11}}}} - \dfrac{{{y^ + }}}{{11}}{{\rm{e}}^{ - \tfrac{{{y^ + }}}{3}}}} \right) $

(48)

与经典对数律函数 ${F_{\log }} = {{\ln {y^ + }} / \kappa } + B$ 相融合，采用如下的混合函数：

$\begin{split} {F_{{\rm{Rei,}}m}}\left( {{y^ + }} \right) =& \left( {1 - {\phi _{{\rm{bl}}}}} \right){F_{{\rm{Rei}}}} + {\phi _{{\rm{bl}}}}{F_{\log }} \\ {\phi _{{\rm{bl}}}} = &\tanh \left( {{{\arg }^4}} \right),\;\arg = {y^ + }/27 \end{split} $

(49)

以及带参数的Spaldings^[9]壁面律可以表示为：

${y^ + } = F_{{\rm{Sp,}}N}^{ - 1}\left( {{u^ + }} \right) = {u^ + } + {{\rm{e}}^{ - \kappa B}}\left[ {{{\rm{e}}^{\kappa {u^ + }}} - \sum\limits_{n = 0}^N {\frac{{{{\left( {\kappa {u^ + }} \right)}^n}}}{{n!}}} } \right]$

(50)

为了避免在过渡层中ω偏离其低雷诺数的解，将标准混合函数（式（45））替换为：

$\begin{split} \omega = &{\phi {\omega _{b1}} + \left( {1 - \phi } \right){\omega _{b2}}} \\ \phi =& {\tanh \left( {{{\arg }^4}} \right),\;\arg = {y^ + }/10} \end{split}$

(51)

其中：

${\omega _{b1}} = {\omega _{{\rm{vis}}}} + {\omega _{\log }},\;{\omega _{b2}} = {\left( {\omega _{{\rm{vis}}}^{6/5} + \omega _{\log }^{6/5}} \right)^{5/6}}$

${\omega _{{\rm{vis}}}} = \frac{{6\mu }}{{{\beta _\omega }\rho {y^2}}},\;{\omega _{\log }} = \frac{{{u_\tau }}}{{\sqrt {{\beta _k}} \kappa y}}$

上述壁面函数均基于临近壁面计算单元上的速度（可压缩情况下，还包括温度）构造，可以适用于基于一方程和二方程湍流模型的RANS模拟，正如方程(11)，基于湍动能的速度壁面函数，比较方便用于包含湍动能的湍流模型RANS模拟中，但以上壁面函数都只适用于不可压缩流动，而没有考虑可压缩性和压力梯度的影响。

国内Gao等^[24]在RANS模拟中采用壁面函数时，着重研究了壁面函数与RANS数值解之间的一致性问题，计算结果表明：如果壁面函数与RANS数值解一致性好，则相应的近壁面网格限制条件可由 ${y^ + } \leqslant 100$ 放宽到 ${y^ + } \leqslant 400$ 。

最近提出的“子网格壁面函数” ^[40-42]，将第一层网格进一步细分成一个子网格，在第一层网格区域的子网格上数值求解边界层方程，而不再追求获得其解析解，文献作者声称与同等密度的单一网格求解方法相比，计算效率可以提高一个量级。显然，这种做法，与Tidriri^[36]的区域分解思想一致，可以避免由于简化得到的解析壁面函数与实际复杂情况下的边界层流动解一致性很差的问题，从发展趋势上讲，是值得进一步研究的技术途径。

5 结　论

围绕发展自主可控CFD软件的目标，以建立高超声速湍流边界层工程实用模拟方法为目的，从湍流壁面函数是湍流边界层方程近似解的角度，对相关公开文献的研究工作进行了梳理，得到以下认识：

1）解析形式的壁面函数，在复杂工程问题求解中已经得到了广泛的应用。但对于以高超声速流动为代表的具有强压缩性、显著传热和大压力梯度的流动，现有的表达形式同微分方程系统细密网格上解的一致性问题没有得到很好的解决，可能会导致计算结果存在明显的误差，特别是压力梯度比较大而壁面剪切应力比较小的区域，如分离、再附点附近的流动。

2）“子网格”壁面函数，其实质是通过数值求解边界层方程得到的数值形式的壁面函数，能够比较充分地反映流动可压缩、传热和压力梯度以及物面几何特征等的影响，对复杂飞行器的壁面湍流流动模拟应该是一种能够同时兼顾计算精度和效率的解决方案，但如何高效求解是“子网格”壁面函数方法在复杂问题中发挥作用的关键。

壁面函数的研究已有将近60年的历史，本文仅仅是对RANS湍流边界层壁面函数的部分研究成果的归纳总结。在目前的计算机资源条件下，系统研究“子网格”壁面函数，并与数据挖掘技术相结合，有可能发展出更加普适的解析壁面函数，以满足对包含强可压缩性、显著传热和强压梯度以及化学反应等复杂物理流动现象和复杂几何构型的壁面湍流流动的高效高精度模拟需求。

另外，需要特别指出的是：近年来大涡模拟在科学研究和工程设计中应用越来越广泛，对于简单机翼绕流算例，当雷诺数Re = 1×10⁹时，壁模型大涡模拟的网格量相较于近壁区分辨（wall-resolved）大涡模拟可以减少3个数量级^[43]，也涌现了大量关于壁模型的研究，因此，壁模型的研究对于大涡模拟在工程应用中的推广具有积极意义，相关的综述性文章可参考文献[44-46]。我们也将会持续关注该方向的研究进展。