湍流是流体动力学中与层流相对应的一种运动形式, 其特征是在平均运动的基础上, 又叠加了一种以流体微团的形式随机运动。研究湍流运动对减小能量耗散和提高能量传递速度, 以及加速化学反应率和提高热交换率等具有重要意义, 且湍流混合对于海洋中氧、盐和热分布至关重要^[1]。湍流早期研究多使用雷诺平均模型, 为了减少计算成本和忽略最小的长度尺度, 人们提出大涡模拟模型^[2-3], 并运用到海洋垂直混合研究^[4-5]和海洋底边界层研究^[6-8]。湍流动能收支往往很难进行直接观测, 一般用湍流动能平衡方程中的每一项作为特征量, 如速度输运、剪切生成、压力输运、湍流输运和耗散等, 这些特征量在湍流研究中有特别关键的作用。Inall等^[9]讨论地形上分层流的湍流动能产生和耗散的演化和分布。刘欢等^[10]指出可以通过一些特征量来描述海底边界层的特性, 认为控制湍流能量收支平衡的因素之一为湍流耗散率。在湍流动能收支项中, 湍流能量耗散的研究一直受到广泛关注^[11-13]。

Calhoun和Street^[14]使用大涡模拟研究波形面上的层流, 比较三种不同海底地形高度的影响, 提出在许多地球物理流动中, 波形面上的湍流显著影响动量通量、混合、和输运等量。Grigoriadis等^[15]用大涡模拟分析了波纹床面上湍流边界层, 提出波纹陡度对流动分离、涡度动力学、湍流、壁应力和阻力的影响很强, 对相对强度的影响则较弱。最近, Jalali等^[16]采用三维高分辨率的大涡模拟, 模拟了作用在一个孤立的超临界障碍上的内波和湍流, 设置一个平滑三角形形状、具有超临界斜率的障碍, 被视作是一个实验室模型的二维海岭, 湍流主要产生于射流中的剪切、反射流的对流不稳定性及瞬态背风波的破碎, 能量收支分析表明随着其定义的无量纲数的增加, 能量转换大幅下降以及局部能量损失大幅增加。湍流是底边界层中显著的海洋动力过程, 是浅海底边界层垂直混合的重要动力因子, 研究分析海底地形对湍流, 尤其是湍流动能收支的影响, 对深入了解海洋混合过程有着重要的意义。为此, 本文通过数值模拟探讨地形对湍流能量收支的影响及湍流耗散随地形改变的垂向结构变化。

1 模式简介

并行大涡模拟模型(The Parallelized Large-Eddy Simulation Model, PALM)是由德国莱布尼茨大学汉诺威气象研究所开发的一个用于大气和海洋流动的大涡模拟模式, 是特别设计用于执行大规模并行计算机体系结构^[17]。PALM是基于Fortran 95以及Fortran 2003的代码且已经被应用于多种大气和海洋边界层的模拟。大涡模拟是计算流体力学中湍流的数学模型。大涡模拟主要思想是减少计算成本和忽略最小的长度尺度, 通过低通滤波求解Navier-Stokes方程。

本文使用模型基于Boussinesq近似下的非流体静力学, 求解经过滤波的不可压缩的Navier-Stokes方程。在大涡模拟, 速度和密度场通过一个空间滤波操作, 分离为大尺度和次网格小尺度范围。在下面的方程组中, 上划线表示滤波后除去次网格项的值。下标0表示表面值。方程中的变量被离散化隐含地滤波, 但是为了方便起见, 这里使用连续形式的等式。双撇号表示次网格尺度(subgrid-scale, SGS)变量。在笛卡尔网格上, 得到通过网格体积滤波的质量, 能量, 位温和盐度的守恒方程。忽略科式力项的大涡模拟模型的基本方程为:

$ \frac{{\partial {u_i}}}{{\partial t}} = - \frac{{\partial {u_i}{u_j}}}{{\partial {x_j}}} - \frac{1}{{{\rho _0}}}\frac{{\partial {{\rm{ \mathsf{ π} }}^*}}}{{\partial {x_i}}} + \frac{\mu }{\rho }\left[{\frac{{{\partial ^2}{u_i}}}{{\partial x_k^2}} + \frac{1}{3}\frac{\partial }{{\partial {x_i}}}\left( {\frac{{\partial {u_k}}}{{\partial {x_k}}}} \right)} \right], $

(1)

$ \frac{{\partial {u_j}}}{{\partial {x_j}}} = 0, $

(2)

$ \frac{{\partial \theta }}{{\partial t}} = - \frac{{\partial {u_j}\theta }}{{\partial {x_j}}} - \frac{\partial }{{\partial {x_j}}}\left( {\overline {{{u''}_j}\theta ''} } \right) - \frac{{{L_V}}}{{{c_p}\varPi }}{\varPsi _{{q_v}}}, $

(3)

$ \frac{{\partial {q_v}}}{{\partial t}} = - \frac{{\partial {u_j}{q_v}}}{{\partial {x_j}}} - \frac{\partial }{{\partial {x_j}}}\left( {\overline {{{u''}_j}{{q''}_v}} } \right) + {\varPsi _{{q_v}}}, $

(4)

$ \frac{{\partial S}}{{\partial t}} = - \frac{{\partial {u_j}S}}{{\partial {x_j}}} - \frac{\partial }{{\partial {x_j}}}\left( {\overline {{{u''}_j}S''} } \right) + {\varPsi _S}. $

(5)

这里, 模型使用笛卡儿坐标系, $i, j, k \in \left\{ {1, 2, 3} \right\}$。${u_i}$是速度分量(${u_1} = u, {u_2} = v, {u_3} = w$), 坐标${x_i}\left( {{x_1} = x, {x_2} = y, {x_3} = z} \right)$, t是时间, $\rho $是海水密度, p是液体静压力, $e = \frac{1}{2}\overline {{{u''}_i}{{u''}_i}} $是次网格湍流动能, g是重力加速度。$\frac{\mu }{\rho } = \upsilon $是运动黏度, $\mu $是动态黏度, $\theta $是位温, ${q_v}$是比湿度, S是盐度。$\varepsilon $是列维-齐维他符号。${\varPsi _{{q_v}}}$和${\varPsi _{{\rm{Sa}}}}$分别是${q_v}$和Sa的源和汇。双撇号代表SGS, 上横杠代表平均值。

本文数值区域设计为100.0×100.0×100.0 m³(x, y, z方向长度), 空间分辨率为1 m, 计算时间为7 200 s, 时间步长为1 s, 海表位温为275 K, –50 m≤z≤0时垂直梯度为–0.015 K/m, –100 m≤z≤–50 m时垂直梯度为–0.01 K/m。海水盐度为32.0, –50 m≤z≤0时垂直梯度为0.015 /m, –100 m≤z≤–5 m时垂直梯度为0.01 /m. x方向的背景流为1 m/s, y方向的背景流为0, 科式力被忽略, 粗糙度为0.1。

本文使用三角地形来简化模拟海底地形情况, 地形结构是完全密闭、不可透过的, 符合海底山丘的简化特征。在空间直角坐标系中, 设计三维地形模型, 地形底边半长用r表示, 地形高度用h表示, 水深沿z方向, 地形表面光滑。采用三角地形斜坡的坡陡作为地形参数, 记为δ, δ=h/r, 得到的参数δ是无量纲的。在本文中, 设置三种地形: δ=0.5, 1, 2。定义δ值为1时为临界地形, 小于1时为亚临界地形, 大于1为超临界地形, 设置的三种地形恰好代表了三种地形状态。三种地形设置的具体参数值见表 1所示。三种地形设置在x-z方向上的剖面如图 1所示。

在这里主要介绍对模拟区域x轴和z轴的无量纲化。对于x轴坐标, 本文采取的措施是除以x方向上地形的特征量半底边长r; 对于z轴坐标, 采取的是除以z方向上水深H。

根据公式(1), 通过张量的收缩, 可以得到水平平均的SGS-TKE方程^[18], 如下:

$ \begin{array}{l} \underbrace {\frac{{\partial e}}{{\partial t}}}_C = \underbrace {- w\frac{{\partial e}}{{\partial z}}}_{{T_m}}\underbrace {- \left( {\overline {u''w''} } \right)\frac{{\partial \bar u}}{{\partial z}}}_S + \\ \underbrace {\frac{{\rm{g}}}{{{\rho _{\theta, 0}}}}\overline {w''\rho _\theta ^{''}} }_B\underbrace {- \frac{\partial }{{\partial z}}\left[{\overline {w''\left( {e + \frac{{p''}}{{{\rho _0}}}} \right)} } \right]}_{P, T} + \underbrace {\upsilon \frac{{{\partial ^2}k}}{{\partial w\partial w}}}_M\underbrace { -\epsilon }_D, \end{array} $

(6)

这里, 双撇号代表SGS, 上划线代表平均值。e为湍流动能。对于${\rho _\theta }$, 它可由海水状态方程计算。${\rho _\theta }$的值取决于盐度Sa、温度T, 和压力p。为简化, 在模拟中${\rho _\theta }$值是常数1.025 7×10³ kg/m³。在公式(6)右边, 第一项代表平均速度输运, 记为${T_{\rm{m}}}$。第二项代表剪切生成, 记为S。第三项代表浮力生成, 记为B。第四项和第五项代表压力输运和湍流输运, 记为P和T。第六项代表分子黏性扩散, 记为M。最后一项代表耗散, 记为D。方程等式左边项记为C, 代表SGS-TKE时间变化率。在本文中, 由于给定恒定密度${\rho _\theta }$, 浮力生成(B)为零。方程右边第六项(M)比右边其他项小若干个数量级, 所以它可以被忽略。整个方程显示SGS- TKE随时间的变化率(C)等于平均速度输运(${T_{\rm{m}}}$)、剪切生成(S)、输运(P, T)和耗散(D)之和。

本文对湍流动能平衡方程进行了无量纲化处理, 分析方程每一项的物理意义及量纲, 对此, 采用的方法是方程两边同除$\frac{{v_0^2}}{{{t_0}}}$, ${v _0}$=1 m/s, 为背景流速, ${t_0}$=3 600 s, 为模拟时长7 200 s(2小时)的一半。对湍流动能收支每一项进行这样的无量纲处理, 可以减弱背景流速初始条件设置以及模拟时长对结果的影响。

2 数值模拟结果 2.1 湍流动能时间序列

在分析湍流动能收支之前, 有必要对其中涉及的变量进行单个分析。通过提取模型输出的全场平均时间序列量, 对湍流动能e进行绘图, 横坐标为时间, 纵坐标为湍流动能e, 结果如图 2a所示。考虑到模式模拟有一个启动过程, 开始的运行结果不太稳定, 通过时间序列图可以看到在模拟三种地形情况时虽然在模式启动时都出现了不稳定的情况, 但持续时间并不长, 所有案例在1 000 s之后都趋于稳定状态, 没有大的波动, 基本符合理论情况。因此在结果呈现中, 各物理量均为模式输出的3 600s和7 200s瞬时值的平均。

在湍流动能平衡方程(6)中, 在等式右边的项有一个变量常常用到, 即动量通量w″u″, w″u″为x方向次网格垂直动量通量, 对应雷诺模型中的雷诺应力, 是脉动运动的平均动量输运。通过提取模型输出的时间序列量, 对动量通量w″u″进行绘图, 横坐标为时间, 纵坐标为动量通量w″u″, 结果如图 2b所示。

2.2 流体速度

x-z剖面x方向速度分布情况如图 3所示, x-z剖面y方向速度分布情况如图 4所示, x-z剖面z方向速度分布情况如图 5所示。第一行、第二行、第三行分别表示在地形坡陡$\delta = 0.5$, 1, 2条件下速度分布情况。本文选取了三维空间中y=5, 20, 35 m三个区域截面, y=20 m为模拟区域在y方向的中心面, 然后对称选取了y=5 m和y=35 m这两个截面作为代表展示。第一列、第二列、第三列分别表示在y=5, 20, 35 m时x-z截面速度分布情况。可以直观对比不同地形坡陡和不同y方向剖面的速度分布情况。

由图 3可以看出整个区域可以分为三个部分:海表附近部分, 地形顶部至海底部分, 以及地形顶部至海表附近部分。在地形顶部至海底区域, 地形的作用使得地形左右两边均产生了回流, 回流最大值可达到–2 m/s, 地形坡陡越大, 回流速度也越大。在地形顶部至海表附近这部分, 速度变大, 尤其是在超临界地形情况最为明显, 最大速度可达5 m/s。海表附近部分速度梯度线分布较为密集。

由图 4可以看出, 在速度y方向上, 速度变化梯度较为密集, 变化也特别紊乱没有次序, 速度梯度在地形顶部梯度达到最大, 由于速度的黏滞性也会产生能量的耗散。在亚临界地形($\delta = 0.5$), 迎风坡和背风坡的流速情况基本呈对称情况, 速度波动不大, 最大值和最小值的绝对值均未超过1 m/s, 分布特别紊乱。在临界地形δ = 1, 在y方向剖面的左右两侧, 地形迎风坡速度分别是负的和正的, 也就是流体速度到达迎风坡后歧分两支, 向两侧流出。在超临界地形($\delta = 2$), 也出现了速度y分量在y轴剖面于迎风坡一正一负情况, 不过情况更为紊乱。

对于速度z分量, 初始背景流速为0 m/s, 这一方向的速度剪切主要是由地形设置产生。在亚临界地形($\delta = 0.5$)和临界地形($\delta = 1$), 迎风坡和背风坡两侧出现略微的速度反向, 在迎风坡水流遇到地形后向下走, 速度z分量为负, 不过速度不是很大。对于超临界地形, 速度变化特别明显, 且速度等值线不像其他两种地形那么密集, 出现成块特征, 可以看到受地形影响较为强烈, 迎风坡前区域速度变为负值, 可达到-2 m/s, 其影响可一直延伸至海表区域。

2.3 流场分布

与速度分布情况对应, 流场分布情况也从x-z剖面呈现, 将与这两个坐标对应的速度x、z分量进行矢量合成, 分析流场分布。本文选取了y方向上中间的剖面y=20 m进行呈现。

在亚临界地形($\delta = 0.5$)下的速度流场分布如图 6a所示, 在水深–0.5以上影响不是特别大, 有一些流动扰动。在水深–0.5以下, 尤其是在地形顶部到海底, 速度大小和方向都有特别明显的影响, 湍流特征突出。在地形底部迎风坡左边, 出现了一个小的涡旋。在地形背风坡右侧区域, 速度较小, 速度方向出现了一个回流涡旋。在峰顶区域有一股比较明显的流。在临界地形($\delta = 1$)下的速度流场分布如图 6b所示, 可以看出由于地形坡陡增加, 对水流影响更为强烈, 海洋上层流体紊乱。在水深–0.6以下, 水流遇到地形底部迎风坡后向下流, 左侧出现了涡旋, 涡旋大小比亚临界地形大许多, 涡旋位置也可以通过旋度看出, 旋转程度较大。在地形背风坡底部, 回流的流体向上涌, 速度较小。同样地, 地形顶部对流场影响很大。在超临界地形($\delta = 2$)下的速度流场分布如图 6c所示, 随着坡陡增大, 地形高度已经设置特别大, 占据水深的近80%, 很接近海表区域, 对于整个水域的影响比较大。在地形迎风坡左侧中上部流体向下运动, 回流形成一个涡旋。在三种地形情况下, 亚临界、临界、超临界地形迎风坡左侧均出现了涡旋, 而且涡旋的位置随着坡陡的增大从海底部逐渐上升, 形状大小也越来越大。

2.4 湍流动能收支

图 7展示湍流动能收支情况, 红线表示速度剪切, 玫红色点画线表示输运, 黑色虚线表示速度输运, 蓝色虚线表示耗散。右图均为左图在地形顶部附近的局部放大图。

在图 7a, 对于亚临界地形下水平湍流动能收支, 有两个特别的区域需要关注, 当z接近于海表和z=–0.815(底地形顶所在高度)。在这两个区域, 剪切生成和耗散的值突然变化, 而压力和湍流输运改变较小。该接近海面区域的变化是由于在海表面复杂的运动。同时, 在底地形顶部的变化是由于地形的剪切等产生剪应力和耗散。在底部山顶的区域, 变化率最大, 与山顶点的位置完全相同。当生成和耗散的运动发生, 必须有压力和湍流输运的变化来传递能量。对于临界地形和超临界地形下水平湍流动能收支, 和亚临界地形类似, 当z接近于海表和底地形顶所在高度, 速度剪切和耗散值突然变化, 而压力和湍流输运改变较小。

3 讨论 3.1 地形顶部耗散与地形陡度关系

在地形顶部区域出现了局部强耗散, 为此, 我们选取耗散项进行讨论。三种地形顶部附近耗散值对比情况如图 8a所示, 蓝色点画线、玫红色虚线、红色实线分别表示亚临界地形($\delta = 0.5$)、临界地形($\delta = 1$)和超临界地形($\delta = 2$)下的耗散值。地形和层结流体间的相互作用, 流体的结构受到地形摩擦的影响, 强烈的剪切导致湍流的生成和发展, 导致底边界附近强耗散。在底地形的山脊处能量耗散的变化, 是由于地形的剪切和摩擦, TKE产生剪应力和耗散, 在底部山顶的区域, 变化率最大, 与山顶点的位置完全相同。

探究地形顶点处耗散值(D₁)与地形坡陡的关系, 对这两个参数进行了回归分析, 得出呈指数分布趋势, 为了是结论更具有代表性, 又设置了四种不同坡陡地形进行模拟计算, 结果见表 2。地形顶部耗散与地形坡陡关系回归图如图 8b所示, 得出的指数关系式为:

$ D =-0.371{e^{2.704\delta }}. $

(7)

3.2 海表耗散与地形陡度关系

除了在地形顶部区域的湍流动能收支情况明显增大, 在海表处的剧增也值得关注, 该接近海面区域的变化是由于在海表面复杂的运动, 速度梯度较密集, 速度黏性作用下TKE产生剪切应力和耗散。三种地形下海表处耗散对比情况如图 9a所示, 蓝色点画线、玫红色虚线、红色实线分别表示亚临界地形($\delta = 0.5$)、临界地形($\delta = 1$)和超临界地形($\delta = 2$)下的耗散值。

探究海表处耗散值(D₂)与地形坡陡的关系, 对这两个参数进行了回归分析, 得出呈指数分布趋势, 同样地, 又增加了四种不同坡陡地形进行模拟, 结果见表 3。海表处耗散与地形坡陡关系回归图如图 9b所示, 得出的指数关系式为:

$ D =-6.446{e^{2.514\delta }}. $

(8)

4 结论

本文使用大涡模拟模型, 以坡陡作为无量纲地形参数($\delta $), 设置了亚临界、临界和超临界三种地形($\delta = 0.5$, 1和2)三种地形, 模拟了海底地形对湍流动能收支的影响。对速度和湍流动能收支的分布, 以及湍流动能平衡方程各参量进行了分析。

(1) 速度在x、y和z方向分布的分析结果都展示了在地形顶点附近的流有大的速度梯度, 所以它意味着在地形顶点附近区域有很强的剪切作用, 导致大的湍流动能生成。地形迎风坡和背风坡附近有较大的速度变化, 通过进一步分析速度流场分布, 三种地形在迎风坡左侧均回流形成涡旋, 而且涡旋的位置随着坡陡的增大从海底部逐渐上升, 形状也越来越大。

(2) 湍流动能收支的分析结果表明除了靠近海洋表面的区域由于复杂的海表运动会产生速度底部山顶区域也出现湍流。正如结果表明, 当浮力生成(B)被忽略, 剪切生成(S)支配湍流动能。然而, 平均速度输运(T_m)是不重要的。此外, 剪切生成在山顶处有最大值, 并且在背风坡比在迎风坡更大。在每个水平层, 湍流动能收支各项S, T_m, (P, T)和D之和接近于0。

(3) 通过回归分析分别探究了地形顶部和海表处耗散, 与地形坡陡的关系, 得出呈指数形式变化, 地形顶部耗散与地形坡陡关系式为: $D =-0.371{e^{2.704\delta }}$; 海表处耗散与地形坡陡关系式为: $D =-6.446{e^{2.514\delta }}$。可为工程实践提供理论依据。