基于κ_sgs与θ_sgs²输运的中性大气边界层大涡模拟

Large-Eddy Simulation of Neutral Atmospheric Boundary Layer Based on κ_sgs and θ_sgs² Transports

摘要: 大气边界层(ABL)中标量的仿真对于掌握大气中温度、湿度等的分布，污染物的扩散规律及预测雨水的生成等均具有重要意义。大涡模拟(LES)是当前ABL湍流仿真的主要方法，其关键是如何利用解析尺度的速度场和标量场信息构建亚网格(SGS)应力模型和通量模型。对于ABL湍流这一典型的高雷诺数条件下的非各向同性湍流，以Boussinesq假设为出发点构建的粘性SGS模型存在先验误差较大、后验耗散过强等许多问题。本研究引入一种新的非线性SGS通量模型，该模拟不预先假设能量的传输方向，基于速度和标量的梯度计算SGS通量的结构(矢量分量的相对大小)，且摈弃局部平衡假设，采用动力学方程考虑能量反传，预测SGS通量强度的动态演化过程。研究采用一个中性ABL基准工况，通过与已有的理论预测和各种流动统计的对比进行系统的考察来对模型进行评估，同时也与传统粘性SGS模型的表现进行比较。具体地，当LES达到统计稳定后，我们重点关注了在不同网格条件下模型对无量纲速度梯度、无量纲标量梯度、能谱及流场结构等的预测结果。结果表明，除了可以获取可靠的流场结构，相比传统粘性SGS模型，新模型对无量纲梯度预测的准确度有明显提升，且对湍流能谱的预测有显著改善。此外，我们讨论了新模型预测效果提升的原因，相较于传统粘性SGS模型存在耗散较强的问题，新模型采用动态非线性的建模方法，可以预测ABL湍流中能量的逆向输运，并更好地捕捉小尺度的涡旋。

Abstract: The simulation of scalars in the atmospheric boundary layer (ABL) is of great significance for understanding the distribution of temperature, humidity, and other factors in the atmosphere, as well as the diffusion patterns of pollutants and the prediction of rain generation. Large-Eddy Simulation (LES) is currently the main method for ABL turbulence simulation, and the key is how to construct sub-grid-scale (SGS) stress model and flux model by using the velocity and scalar field information of analytic scale. For ABL turbulence, a typical non-isotropic turbulence under high Reynolds number conditions, the viscous SGS model based on the Boussinesq hypothesis has many problems, such as large prior errors and strong posterior dissipation. This study introduces a new nonlinear SGS flux model, which simulates the structure of SGS flux (relative magnitude of vector components) based on velocity and scalar gradients without assuming the direction of energy transfer. In addition, the local equilibrium hypothesis is discarded, and kinetic equations are adopted to consider the reverse energy transfer and predict the dynamic evolution of SGS flux intensity. The study used a neutral ABL benchmark operating condition and systematically evaluated the model by comparing it with existing theoretical predictions and various flow statistics, as well as comparing its performance with traditional viscous SGS models. Specifically, when LES reaches statistical stability, we focus on the prediction results of the model for dimensionless velocity gradient, dimensionless scalar gradient, energy spectrum, and flow field structure under different grid conditions. The results show that in addition to obtaining reliable flow field structure, compared with the traditional viscous SGS model, the new model significantly improves the accuracy of dimensionless gradient prediction and turbulence energy spectrum prediction. In addition, we discussed the reasons for the improved prediction performance of the new model. Compared with the traditional viscous SGS model, which has a strong dissipation problem, the new model adopts a dynamic nonlinear modeling method, which can predict the reverse energy transport in ABL turbulence and better capture small-scale vortices.

1. 引言

大气边界层(Atmospheric Boundary Layer, ABL)是指靠近地球表面的大气分层，其厚度通常情况下在300米到3公里之间，动植物生存、人类活动以及风电场运行等等都处在ABL中，因此对ABL的仿真对于掌握大气中温度、湿度等的分布，雨水的生成，污染物的扩散规律，以及风能的利用等均具有重要意义。因为ABL湍流的雷诺数基本上都在10⁸以上，对其仿真求解仍依赖于湍流模型。大涡模拟(Large-Eddy Simulation, LES)是当前ABL仿真的主要方法之一，其关键是如何利用解析尺度的速度场和标量场信息构建亚网格(Sub-Grid-Scale, SGS)应力模型和通量模型[1]，

${\tau }_{ij}=\widetilde{u_i{u}_j}-\widetilde{u_i}\widetilde{u_j}$ (1)

$q_i=\widetilde{u_i\theta }-\widetilde{u_i}\tilde{\theta }$ (2)

其中，三维空间过滤操作在文中表示为波浪线(~)，${\tau }_{ij}$ 为SGS应力，$q_i$ 为SGS通量。对于均匀各向同性湍流这种简单的理想化湍流来说，不同SGS模型给出的差异有限，但对于ABL这种典型的高雷诺数高剪切的非各向同性湍流来说，SGS模型对结果的影响较大[2] [3]。已有的研究表明，涡黏性SGS模型在ABL模拟中暴露出不少问题，包括对小尺度能量耗散过快，无量纲梯度过大等等[4] [5]。这都与ABL的特殊性质与涡黏性模型基本建模假设的局限有关，包括ABL在内的边界层湍流作为非各向同性湍流的典例，存在像发卡涡、大尺度低速条带等等大尺度的相干结构。而要在仿真中捕捉这些结构，务必要求湍流模型能够提供从小尺度到大尺度的逆向传输，且不假设湍动能量的局部平衡。

Smagorinsky模型[6]以及对应的通量模型是最早提出的SGS模型，

${\tau }_{ij}-\frac{1}{3}{\delta }_{ij}{\tau }_{kk}=-2v_{sgs}{\tilde{S}}_{ij}$ (3)

$q_i=-\frac{{\nu }_{sgs}}{P{r}_{sgs}}\frac{\partial \tilde{\theta }}{\partial x_i}$ (4)

其中，${\tilde{S}}_{ij}$ 是应力张量，${\nu }_{sgs}={\left(C_S\tilde{\Delta }\right)}^2\left|\tilde{S}\right|$ 为各向同性的SGS黏性项，$\left|\tilde{S}\right|=2{\tilde{S}}_{ij}{\tilde{S}}_{ij}$ 为应力强度，$C_S$ 为常系数，$P{r}_{sgs}$ 是SGS普朗特数(对应温度场，或者$S{c}_{sgs}$ 对应其他标量场)。以该模型为代表的黏性SGS模型在湍流建模理论上需要依赖的基础假设包括：1) 亚网格尺度流动的影响主要表现在湍流能量上，因此，采用能量传递的平衡关系来描述尺度之间的相互作用；2) 从解析大尺度到亚网格尺度的能量传递机制沿用以扩散表示的分子黏性机制；3) 尺度分离假说假设亚网格尺度与解析尺度完全分离，因此它们之间的相互作用只能通过从较大的解析尺度到较小的亚网格尺度的能量传输来模拟；4) 局部平衡假说假设湍流场中任何区域的湍动能耗散都与大尺度到小尺度的能量传输平衡，这常常和尺度相似性假设一起，被用于确定黏性模型的系数，它与第二和第三个假设一起，常常在被修正之后用于对黏性模型的改进研究；5) 模型广泛采用各向同性假设，在不同方向上，采用同样大小的湍流黏性。务必要提及的是，前两个假设继承了雷诺平均(Reynolds-Averaged Navier-Stokes, RANS)方法的主要建模思想，即Boussinesq假设，用黏性项来建模。黏性假设使得模型项与精确项之间存在较大误差，并且模型项与精确项的坐标变换不一致性[7]，再加上各向同性假设和对于能量传递机制的假设无法满足非各向同性湍流物理特性，使得在包括ABL的许多湍流仿真中黏性模型给出的结果与实际测量有较大的偏差，例如在ABL表层，速度和温度的无量纲垂直梯度可能被高估20%以上[8] [9]。

非线性模型(也称为结构型模型)是SGS建模的另外一个方向[1]。非线性模型对不同湍流尺度之间的相互作用不做人为预设的假设，而主要关注如何精确地模拟SGS项。早期的非线性模型包括相似性模型[10]和梯度模型[11] [12]等都表现出较强的不稳定性，主要是因为它们对湍动能大小的预测不够准确，因此错估跨尺度的能量传输快慢强弱[13]。经过改进之后，许多高级模型可以在许多复杂条件下的湍流模拟中给出高精度的大涡模拟[2]，尤其是在ABL研究当中，结果表明它们能够给出比黏性SGS模型更好的预测[7] [14] [15]。

本研究将采用中性ABL基准工况，考察SGS湍动能和SGS标量方差的输运方程耦合梯度型结构SGS模型在ABL仿真中的效果。采用中性ABL工况可以有效地解耦速度场与标量场之间的耦合关系，标量场以被动标量的形式出现在ABL工况中，因此也扩展了标量的通用性，可以代表温度、湿度、污染扩散研究中的化学组分等等。虽然SGS湍动能输运方程与梯度型结构模型的耦合在各向同性湍流、旋转湍流[7] [16]，以及复杂的湍流燃烧、喷雾和喷雾燃烧模拟中已被证明能够给出低网格敏感且高精度的LES结果，但在以ABL为代表的边界层湍流中还没有得到有效的论证，这是本研究的首要任务。更重要的是，本研究将通过引入SGS标量方差的输运方程这一方法来对大气边界层条件下的SGS通量模型进行建模，并通过与经典的SGS黏性模型进行比较等方式对新模型进行考察。

2. 基于$k_{sgs}$ 与${\theta }_{sgs}^2$ 输运的梯度型结构模型

在LES中，结果的准确性与选用的亚网格模型有很大关系。本研究提出采用非线性梯度型结构亚网格应力模型和亚网格通量模型来求解ABL，

${\tau }_{ij}=2k_{sgs}\left(\frac{{\tilde{G}}_{ij}}{{\tilde{G}}_{mm}}\right)+{\nu }_{sgs}^u{\nabla }^2{\tilde{S}}_{ij}$ (5)

$q_i=\left|\text{q}\right|\left(\frac{{\tilde{G}}_{\theta ,i}}{\left|{\tilde{G}}_{\theta }\right|}\right)+\frac{{\nu }_{sgs}^u}{S{c}_{sgs}}{\nabla }^2\frac{\partial \tilde{\theta }}{\partial x_i}$ (6)

其中，$k_{sgs}=\frac{1}{2}\left(\widetilde{u_i{u}_j}-\widetilde{u_i}\widetilde{u_j}\right)$ 为亚网格动能，$\left|\text{q}\right|$ 为亚网格通量的大小，$\tilde{\theta }$ 为滤波后的任意标量，${\tilde{G}}_{ij}=\frac{{\tilde{\Delta }}^2}{12}\frac{\partial {\tilde{u}}_i}{\partial x_m}\frac{\partial {\tilde{u}}_j}{\partial x_m}$ 和${\tilde{G}}_{\theta ,i}=\frac{{\tilde{\Delta }}^2}{12}\frac{\partial \tilde{u}{}_i}{\partial x_m}\frac{\partial \tilde{\theta }}{\partial x_m}$ 为梯度项(考虑到网格的各向异性，本研究中梯度张量和梯度向量由$G_{ij}=\frac{\tilde{\Delta }{x}^2}{12}\frac{\partial {\tilde{u}}_i}{\partial x}\frac{\partial {\tilde{u}}_j}{\partial x}+\frac{\tilde{\Delta }{y}^2}{12}\frac{\partial {\tilde{u}}_i}{\partial y}\frac{\partial {\tilde{u}}_j}{\partial y}+\frac{\tilde{\Delta }{z}^2}{12}\frac{\partial {\tilde{u}}_i}{\partial z}\frac{\partial {\tilde{u}}_j}{\partial z}$ 和${\tilde{G}}_{\theta ,i}=\frac{\tilde{\Delta }{x}^2}{12}\frac{\partial {\tilde{u}}_i}{\partial x}\frac{\partial \tilde{\theta }}{\partial x}+\frac{\tilde{\Delta }{y}^2}{12}\frac{\partial {\tilde{u}}_i}{\partial y}\frac{\partial \tilde{\theta }}{\partial y}+\frac{\tilde{\Delta }{z}^2}{12}\frac{\partial {\tilde{u}}_i}{\partial z}\frac{\partial \tilde{\theta }}{\partial z}$ 来进行计算)，$\left|{\tilde{G}}_{\theta }\right|=\sqrt{{\tilde{G}}_{\theta ,1}^2+{\tilde{G}}_{\theta ,2}^2+{\tilde{G}}_{\theta ,3}^2}$ 。同时，研究采用亚网格高黏性${\nu }_{sgs}^u={C^{\prime }}_k{\tilde{\Delta }}^3\sqrt{k_{sgs}}$ ，替代了传统的亚网格黏性，以减小对大尺度湍流的耗散。根据Lu等人[17]的研究，选择经验常数${C^{\prime }}_k=0.008$ 。这一模型将建模分成两部分：归一化梯度向量用于模拟亚网应力张量和通量矢量的结构(每个方向的相对量级)；并且需要一个单独的模型来计算亚网格应力和通量的大小。

进一步地，本研究摒弃了ABL在解析尺度与亚网格尺度之间的局部平衡假设，而采用输运方程来计算亚网格湍动能，

$\frac{\partial k_{sgs}}{\partial t}+{\tilde{u}}_i\frac{\partial k_{sgs}}{\partial x_i}=-\frac{{\tau }_{ij}}{\overline{\rho }}\frac{\partial {\tilde{u}}_i}{\partial x_j}-C_{\epsilon }\frac{k_{sgs}^{3/2}}{\tilde{\Delta }}+\frac{\partial }{\partial x_i}\left({\nu }_{sgs}\frac{\partial k_{sgs}}{\partial x_i}\right)$ (7)

上式从左到右分别表示的是亚格子湍动能的瞬态项、对流项、源项、耗散项和扩散项。其中，${\nu }_{sgs}=C_k\sqrt{k_{sgs}}\tilde{\Delta }$ ，根据Menon等人[18]的研究，经验常数取$C_k=0.05$ ，$C_{\epsilon }=1.0$ 。亚网格通量的大小用下面的方程计算，

$\left|\text{q}\right|=u_{sgs}{\theta }_{sgs}$ (8)

其中，亚网格速度标度与亚网格动能的平方根成正比，即$u_{sgs}=\sqrt{2k_{sgs}}$ ，用于刻画标量脉动程度的标量标准差${\theta }_{sgs}$ 的计算，依赖于亚网格标量方差${\theta }_{sgs}^2=\widetilde{\theta \theta }-\tilde{\theta }\tilde{\theta }$ 求解。对于标量通量的模化，本研究也摒弃局部平衡假设，而是求解其动态的输运方程，

$\frac{\partial {\theta }_{sgs}^2}{\partial t}+{\tilde{u}}_i\frac{\partial {\theta }_{sgs}^2}{\partial x_i}=\frac{\partial }{\partial x_i}\left(\frac{{\nu }_{sgs}}{S{c}_{sgs}}\frac{\partial {\theta }_{sgs}^2}{\partial x_i}\right)-q_i\frac{\partial \tilde{\theta }}{\partial x_i}-\sqrt{2}{C}_{\epsilon \theta }\frac{{\theta }_{sgs}^2\sqrt{k_{sgs}}}{\tilde{\Delta }}$ (9)

本研究将这节展示的基于$k_{sgs}$ 与${\theta }_{sgs}^2$ 输运的梯度型结构模型用于ABL大涡模拟仿真，该方法是一种求解$k_{sgs}$ 与${\theta }_{sgs}^2$ 动力学过程的湍流模化方法，因此缩写为“GDSM (Gradient-type Dynamic Structure Modeling)”。

3. 数值模拟

本研究采用成熟应用于许多ABL数值模拟研究中的LES方法[5] [19]-[24]求解基本守恒方程，包括滤波后的连续性方程、动量守恒方程和标量输运方程。由于ABL的雷诺数基本上都在108以上，湍流的黏性占据了绝对的主导(往往在数量级上大于分子黏性千万倍以上)，因此空气的分子黏性可以忽略不计。目前研究中常用的标准黏性亚网格模型(公式(3))及以此为基础的扩散型通量模型(公式(4)) [6]。最早的涡流黏性模型假设$C_S$ 为常数，即Smagorinsky系数。该模型已被广泛应用，但在高雷诺数大气边界层的大涡模拟中，标准涡黏性模型很难预测表层的平均风速、温度和剪切力，进而会产生相应错误的剖面，其亚网格项与实际项的相关性较低，且该模型在所有方向上使用相同的涡黏度/扩散率[25]。考虑到所提出的模型的新颖性主要来自于其没有使用常用的涡黏性/扩散模型，为了进行比较，本文还将使用具有相似复杂度(零方程、非动态)和计算成本的涡黏性/扩散模型获得的结果。然而，由于流动的各向异性，特别是高雷诺数大气边界层中表面附近存在强平均剪切，使这些系数的最佳值与各向同性的系数不同[14] [26] [27]。一种常见的做法是以特定的方式指定系数。Mason and Thomson [4]提出的特设阻尼函数可以改写成：$C_S={\left(C_0^{-n}+{\left(\kappa \left(\frac{z}{\tilde{\Delta }}+\frac{z_0}{\tilde{\Delta }}\right)\right)}^{-n}\right)}^{-1/n}$ ，其中n为可调参数，且研究已经表明，与使用常系数的Smagorinsky模型相比，当$C_0$ 的值在0.1至0.3这一范围且n = 1, 2或3时，该公式可以在表层提供更真实的对数速度分布，且研究还发现$S{c}_{sgs}$ 值的范围在0.33到0.7之间[4] [5] [25]。根据Port-Agel等人的研究[5]，本次研究中该模型采用的阻尼系数为$C_0=0.17$ 及n = 1 (修正的Smagorinsky模型缩写为“SM”)，并在研究中采用$S{c}_{sgs}=0.5$ 及$S{c}_{sgs}=0.7$ 这两个常数值。

本次模拟的大气边界层水平均匀，其水平方向采用伪谱离散，垂直导数用二阶中心差分近似。计算域的高度为H = 1000 [m]，模拟体积的水平尺寸为L_x = L_y = 2πH。计算域被划分为N_x，N_y和N_z个均匀间隔的网格点，本研究以N_x × N_y × N_z = 32 × 32 × 32, 48 × 48 × 48和64 × 64 × 64的网格分辨率进行模拟。网格平面在垂直方向上交错，其中第一个垂直速度平面距离表面为$z=\frac{H}{N_z-1}$ ，而第一个水平速度平面距离表面为$\frac{\Delta z}{2}$ 。在底部通过应用Monin-Obukhov相似理论计算瞬时壁应力[5] [19]：${{\tau }_{i3}|}_{\varpi }=-u_{\ast }^2\frac{{\tilde{u}}_i}{U\left(z\right)}\varpi =-{\left(\frac{U\left(z\right)\kappa }{\ln \left(z/z_0\right)-{\Psi }_M}\right)}^2\frac{{\tilde{u}}_i}{U\left(z\right)}$ ，其中$\kappa$ 为von 常数，$u_{\ast }$ 为摩擦速度，$z_0$ 为粗糙长度，${\Psi }_M$ 为动量的稳定性校正，而$U\left(z\right)$ 是平面平均分解水平速度。本文使用常见的公式$\tilde{\Delta }=\sqrt[3]{{\Delta }_x{\Delta }_y{\Delta }_z}$ 来计算过滤器大小，其中${\Delta }_x=\frac{L_x}{N_x}$ ，且${\Delta }_y=\frac{L_y}{N_y}$ 。这一研究根据3/2规则在非线性项中校正相应的混叠误差，且时间推进使用二阶精确的Adams-Bashforth方案来执行[28]。

本次模拟采用经典的数值设置[5] [19] [24]。模拟流动由x方向上的恒定压力梯度$-u_{\ast }^2/H$ 进行驱动，并参照之前的研究设置$u_{\ast }=0.45$ [m/s]及$z_0=0.1$ [m] [5] [24] [25]。上边界条件设置为$\partial \tilde{u}/\partial z=0$ ，$\partial \tilde{v}/\partial z=0$ ，$\tilde{w}=0$ 及$\partial \tilde{\theta }/\partial z=0$ 。在底部，中性稳定性结果为${\Psi }_M=0$ 。通过施加恒定的向下表面通量${q_3|}_w=-u_{\ast }{\theta }_{\ast }$ ，将类似于先前研究的被动标量场引入模拟中[19] [25] [29]，并设置${\theta }_{\ast }=0.9$ [K]。

4. 模拟结果

图1显示的是在高度为z/H = 0.1的水平面上LES获得的瞬态流向速度云图，可以看出相对较为狭长的低速条带，高速带则较宽，这一观测和其他模拟[24]的结果是一致的。在数值模拟达到统计学稳定状态后，本研究对两种模型的湍流统计数据进行了收集。在本文中，水平面和时间的综合平均值表示为$\langle \cdot \rangle$ ，而相应的分解变量的波动则表示为${\tilde{f}}^{\prime }=\tilde{f}-\langle \tilde{f}\rangle$ ；本次研究以32³、48³和64³这三种空间解析率下获得的模拟结果为基础来进行检验。

Figure 1. Transient flow velocity cloud map obtained from 64³ LES on a horizontal plane with a height of z/H = 0.1

图1. 高度z/H = 0.1水平面上64³ LES获得的瞬态流向速度云图

4.1. 速度与标量

Figure 2. Non-dimensional velocity profile

图2. 无量纲速度剖面

长期以来，用LES方法求解ABL湍流存在平均风剖面和温度剖面不准确的问题。因此本小节主要将新方法的模拟预测结果与相似性理论结果进行比较，以更好了解新模型的性能。冯卡门[30]于1931首次发表的对数律剖面是一种半经验关系，用于描述湍流边界层表面上方水平风速的垂直分布。该研究指出，湍流边界层中某一点的平均流向速度与该点到表面距离的对数成正比。后面建立的包括热效应的Monin-Obukhov相似理论[8] [9]已经在很多测量中得到了证实，因此新建立的亚网格模型应该与之进行比较。中性ABL条件下的风速剖面可以写成已知的对数定律公式：$\langle \tilde{u}\rangle =\frac{u_{\ast }}{\kappa }\ln \left(\frac{z}{z_0}\right)$ ，其中空气动力学粗糙度$z_0$ 必然非零。对数定律能较好近似ABL湍流在表层的速度分布，其占据大气边界层中最低处z/H = 10%~15%以下的流动。图2比较了在64³网格条件下使用SM和GDSM这两个模型获得的平均流向速度剖面。可以看到，使用SM获得的结果与之前的研究结果一致[5]，并且很明显使用SM获得的平均速度大于对数定律预测的平均速度；而在表层，GDSM能够更有效获得对数剖面，并且随着模拟使用的网格分辨率的增加，风速剖面与对数定律的偏离度随之减少。

为了更严格地评估模型性能，可以检查作为垂直位置函数的解析流向速度的无量纲垂直梯度值。平均解析流向速度的无量纲垂直梯度定义为：

${\Phi }_M=\frac{\kappa z}{u_{\ast }}\frac{\partial \langle \tilde{u}\rangle }{\partial z}$ (13)

根据测量结果和量纲分析[8] [9] [30]，在表层${\Phi }_M=1$ ，因此与对数层之间的未匹配可以清楚显示，且有助于定量评估模型性能。图3显示了两种模型所得的无量纲垂直速度梯度值。与其他研究一致[4] [5]，SM的模拟在z/H = 10%以下流动处产生的${\Phi }_M$ 明显大于理论值1 (最大相对误差约为30%)，说明这一模型仍然过于耗散，在解析场去除了较多的动能，这将导致表层中的过度剪切，进而产生较大${\Phi }_M$ 值。可以注意到，尽管SM在更高的水平高度(z/H > 0.15)可以更好地预测梯度值，但如图2所示，其模拟所得速度的大小也随之增加。考虑到相似性曲线仅适用于大气边界层中最低处z/H = 10%~15%以下的流动。相比之下，在表层GDSM这一新模型产生的${\Phi }_M$ 更接近于理论值(最大相对误差约18%)，因而可以更好地表示预期的对数关系。

Figure 3. Non-dimensional vertical velocity gradient map

图3. 无量纲垂直速度梯度图

对于标量，本次研究主要通过作为垂直位置函数的平均分解标量浓度的无量纲垂直梯度值来进行检验。无量纲标量梯度的定义为：

${\Phi }_{\theta }=\frac{\kappa z}{{\theta }_{\ast }}\frac{\partial \langle \tilde{\theta }\rangle }{\partial z}$ (14)

许多文献表明[8] [9]，在中性大气边界层条件下，ABL表层有${\Phi }_{\theta }=0.74$ 。使用两个不同的$S{c}_{sgs}$ 值进行SM模拟获得的无量纲标量梯度如图4(a)所示，可以看到，标量梯度具有与速度场对应相同的一般特征。与其他研究一致[4]，当前研究还表明，模拟中${\Phi }_{\theta }$ 仍然存在与0.74这一值的偏差，且采用其他涡流扩散模型时也会产生类似的${\Phi }_{\theta }$ 过冲[19] [31]。图4右显示了使用新模型在不同网格条件下进行模拟获得的无量纲标量梯度。由图可知，新模型的模拟显然产生了更好的${\Phi }_{\theta }$ 值，并且结果显示新模型对网格分辨率的敏感性较低，且32³网格条件下的模拟获得的${\Phi }_{\theta }$ 值略接近于预期值，这可能是由于在较粗糙的网格分辨率下数值扩散增加的影响。此外，该模拟在低网格点处略低于标准${\Phi }_{\theta }$ 值，而在高网格点处将大于标准值。Port-Agel [19]之前的研究已经检验了动态涡流粘度模型和动态涡流扩散模型的组合，得出的${\Phi }_{\theta }$ 值在地面附近很小(约0.6)，在表层急剧增加(至约0.9) [32] [33]。

(a) SM (b) GDSM

Figure 4. Non-dimensional scalar gradient simulated by two models: (a) SM; (b) GDSM

图4. 两种模型所得的无量纲标量梯度图：(a) SM；(b) GDSM

4.2. 能量谱

标量场的功率谱具有惯性子区间和耗散区间。在惯性子区间内，能量谱表现为遵循−5/3幂律标度[13]，特别是在中性ABL流动中，惯性子区间应扩展到相对较小的尺度范围，即$k_1\ge z^{-1}$ ，其中z为高度，$k_1$ 为流向波数。需要注意到，耗散区间在高雷诺数湍流的大涡模拟中没有进行解析，因此不需要考虑这个范围。

图5、图6和图7分别是在64³网格条件下使用SM和GDSM进行模拟获得的一维流向速度能谱图、垂直速度能谱图及标量能谱图。可以看到，两种模型的所有模拟结果均遵循−5/3幂律标度。图5(a)、图6(a)和图7(a)分别为64³网格条件下用SM进行数值模拟所获得的无量纲一维流向速度、垂直速度和标量能谱图。通过一维傅里叶变换计算频谱，然后在展向方向和时间上对其进行平均。为了检查惯性子区域内曲线可能的下降，需要将其通过${\theta }_{\ast }^{2}z$ 来进行归一化，并根据$k_{1}z$ 来对能量谱进行绘制。对于相对较小的标度($k_{1}z\ge 1$ )，−5/3幂律标度关系预计会保持不变，能量谱表现出良好的下降情况。在表层，大多数网格分辨率尺度都在惯性子区间之外，能量谱会迅速下降，这表明这种湍流扩散方法会产生过度耗散。相关研究表明[19] [32] [33]，在表面附近，标准动态闭合产生的能谱斜率太平，且在最小解析尺度下具有不符合实际的标量方差堆积。此外也有研究表明[25]，测试的涡流粘度/扩散模型可以在小尺度上产生过度耗散或者明显的功率密度累积。使用涡流扩散模型获得的结果在很大程度上取决于模型系数的值，尤其是在表层，且目前的数值结果也支持这一说法。

图5(b)、图6(b)和图7(b)分别为64³网格条件下用GDSM进行数值模拟所获得的无量纲一维流向速度、垂直速度和标量能谱图。GDSM的模拟显然能够在惯性子区间内实现−5/3幂律标度。在表层，将这些结果与SM的模拟结果进行比较，揭示了新模型两个主要的改进特征：GDSM在功率包含区间(低模式)中产生更大的值，且在滤波器尺度上的耗散明显更小，因此GDSM可以更好地再现速度和标量方差向亚网格尺度的传递速率。

(a) SM (b) GDSM

Figure 5. One-dimensional flow velocity energy spectrum simulated by two models: (a) SM; (b) GDSM

图5. 两种模型所得的一维流向速度能谱图：(a) SM；(b) GDSM

(a) SM (b) GDSM

Figure 6. One-dimensional vertical velocity energy spectrum simulated by two models: (a) SM; (b) GDSM

图6. 两种模型所得一维垂直速度能谱图：(a) SM；(b) GDSM

(a) SM (b) GDSM

Figure 7. One-dimensional scalar energy spectrum simulated by two models: (a) SM; (b) GDSM

图7. 两种模型所得一维标量能谱图：(a) SM；(b) GDSM

4.3. 二阶统计

图8显示了从64³网格的模拟中获得的标准化总剪应力及已分部剪应力(包括已解析的剪应力和亚网格剪应力)的垂直分布，标准化壁面总法向通量及已分部法向通量(包括已解析的法向通量和亚网格法向通量)的垂直分布，以及从两个较粗糙的网格分辨率(32³网格及48³网格)条件下获得的标准化亚网格剪应力和亚网格法向通量。正如预期的那样，较低分辨率条件下模拟产生的亚网格应力和亚网格通量大于较高分辨率条件下对应的模拟结果。模拟的流动由恒定的压力梯度驱动，并在流动中施加恒定的表面通量，因此在没有粘性效应的情况下，标准化总湍流应力和总湍流通量的大小分别从表面的±1近似线性下降到顶部的0。直接数值模拟的研究说明了总湍流应力和总湍流通量特征之间的相似性[34]，表明标量波动也与速度波动一样间歇性产生。此外，总湍流通量的近线性特征与对数区域的直接数值模拟结果和中性ABL流动的大涡模拟结果都能很好地吻合[19] [29] [34]。这些结果也证实了该模型方案的平稳性及动量守恒。

Figure 8. GDSM simulation results: (a) shear stress diagram; (b) wall-normal flux diagram

图8. GDSM模拟所得的：(a) 剪应力图；(b) 壁面法向通量图

5. 结论与展望

准确的ABL仿真对于精准掌握大气中温度、湿度等的分布，污染物的扩散规律及预测雨水的生成等均具有重要意义。传统的基于Smagorinsky模型的涡黏性/扩散SGS模型需要通过假设SGS脉动的产生与消耗速率达成局部平衡来模化SGS应力和通量，这一假设会给以高雷诺数高剪切为基本特征的ABL湍流的仿真带来较强的限制，导致预测结果存在较大偏差。本研究提出不采用局部平衡这一假设，通过引入输运方程这一方式来计算亚网格动能和标量方差，并结合梯度方法计算亚网格应力和通量，不需要额外的测试滤波，使得该方法在计算上十分高效。

本文通过对中性大气边界层进行模拟来检验所提出的模型。众所周知，在亚网格运动占总湍流很大一部分的ABL表层，ABL的大涡模拟对亚网格模型相当敏感。传统的涡黏性模型与表层的Monin-Obukhov相似性形式存在偏差，这种偏差很容易在速度和标量的廓线中观察到，在更大程度上在其无量纲垂直导数中也可以观察到。通过比较可以看到，新方法相对于简单的涡黏性模型有了显著的改善，除了能够获取可靠的流场结构，进一步的结果还表明，该方法可以获得更符合对数定律的风速剖面，新模型无量纲垂直速度梯度预测的相对误差提升了12%，而无量纲标量梯度预测的相对误差则提升了将近20%，且新模型可以更好地再现速度和标量方差向亚网格尺度的传递速率。此外，我们讨论了新模型预测效果提升的原因，相较于传统粘性SGS模型存在耗散较强的问题，新模型采用动态非线性的建模方法，可以预测ABL湍流中能量的逆向输运，并更好地捕捉小尺度的涡旋。

致 谢

数值计算在合肥先进计算中心完成。