数值海洋与大气模式（五）：湍流模式与参数化方案

zyw · 发表于海洋技术专业 2021-11-26 13:00:32

1.参数化方案　　上文提到，POM模式对湍流项采取的是Mellor-Yamada方案。为了讲明白湍流参数化方案的物理意义，首先从参数化方案是什么来讲起。
　　目前为止，前文所讨论的针对模式方程求解所使用的方法均为有限差分法，在使用有限差分法时，首先要做的就是设置时间和空间步长，将时间空间坐标网格化。步长取得越小，网格越密集，分辨率也就越高。然而，不管网格划分的多密，网格点都是离散的，只能无限逼近于连续。这么做就带来一个问题，任何一个离散网格都有一个最小可分辨尺度。在有限差分的模式中，运动的最小可分辨尺度为二倍格距的波长。例如，在使用50-100km分辨率的全球天气预报模式进行模拟时，将无法分辨大气中波长为几百米的波动。因此，可以将运动分为网格尺度和次网格尺度。
　　次网格尺度的运动虽然是小尺度运动，但如果在模式计算中将其直接忽略，结果会产生显著偏差。因其运动尺度远远小于网格的间距，也无法通过其物理机理去精确计算。网格尺度的运动和次网格尺度的运动必然存在一定程度的相互作用，可以通过引入一些参数，来近似建立网格尺度的物理量与次网格尺度的物理量之间的关系，而确定这些参数的方案，被称为参数化方案。在大气与海洋模式中，湍流是次网格尺度的运动。本文将从湍流是什么开始，介绍几种经典的湍流模式。这几个湍流模式十分经典，尽管很少被直接应用于模式之中，但弄明白其物理意义是理解成熟的湍流参数化方案的基础。
2.湍流和雷诺应力　　说起湍流，就不得不提到雷诺数。雷诺数所表征的是惯性力比粘性力，由下式可以看出，分母是粘性系数，分子则包含流动速度项。也就是说雷诺数越大，惯性力就越占主导，越小则粘性力越占主导。著名的雷诺实验，就是通过调整雷诺数，发现了粘性流体的两种流动状态，层流和湍流。

$R_{e}=\frac{\rho U D}{\mu} \\$

登录/注册后可看大图

　　实验发现，当雷诺数大于临界雷诺数时，流体是湍流状态，流场中出现了很多小旋涡。而雷诺数小于临界雷诺数时，流体沿着管道平行的方向做比较平滑的直线运动，动量交换主要靠分子粘性力。不严谨的讲，分子粘性力使得周围流体趋于朝同一方向运动，小雷诺数时，粘性力影响大，而粘性力又会使得速度较大的那一层流体带着速度较慢的那一层流体一起朝同一个方向运动。而雷诺数较大时，相邻层间就会出现滑动，而不是统一步调，不同层间也出现混合现象。
　　为什么存在湍流目前还没有一致的说法，因此不做赘述。但需要理解湍流的几个特点。湍流具有随机性，运动几乎无法预测，但是其平均态又具有确定性。因此，为了更好的研究湍流的特点，把变量分解成平均项和脉动项。下式f可以指代各个变量，f可以是速度u，也可以是v。

$f=\bar{f}+f^{\prime} \\$

登录/注册后可看大图

　　而

$\bar{f}$

登录/注册后可看大图

是时间平均后的结果。

$\bar{f}(x, y, z, t)=\frac{1}{T} \int_{t-T / 2}^{t+T / 2} f(x, y, z, t) \mathrm{d} t\\$

登录/注册后可看大图

　　理解了这两项的物理意义，就理解了其运算规则。

$\begin{array}{l} f+g=\bar{f}+\bar{g} \\ a f=a \bar{f} \\ f \cdot g=\bar{f} \cdot \bar{g}+f^{\prime} \cdot g^{\prime} \\ \partial f / \partial s=\partial \bar{f} / \partial s \\ \bar{f}=\bar{f} \\ f^{\prime}=0 \\ f^{\prime} \cdot \bar{g}=0 \\ \bar{f} \cdot g=\bar{f} \cdot \bar{g} \end{array}\\$

登录/注册后可看大图

　　上述运算规则会在后面推导雷诺应力时用到，若对哪项不理解可以代到时间平均的积分中就能很好的推出来。下面这一项可能略微难理解，在此稍微解释下。因为

$f^{\prime}$

登录/注册后可看大图

是脉动项，从数学的角度来理解的话，此项平均项的波动项，在求时间平均后肯定是为0的，因为这一项所表达的是对平均的波动。

$\bar{f^{\prime}}=0\\$

登录/注册后可看大图

　　湍流的运动尺度虽小，但大于分子尺度。因此控制方程的核心前提，连续介质假设是成立的。所以湍流运动是满足连续方程和NS方程的，然而湍流的随机性导致方程不可解，因为每一个瞬时流动都是不同的。必须找到统计平均态的控制方程，才可以避免随机性。接下来，开始推导雷诺方程。下面用张量的形式来表达方程。

$\begin{array}{c} \frac{\partial u_{i}}{\partial x_{i}}=0 \ \rho \frac{\partial u_{j}}{\partial t}+\rho u_{i} \frac{\partial u_{j}}{\partial x_{i}}=\frac{\partial \tau_{i j}}{\partial x_{i}} \end{array}\\$

登录/注册后可看大图

　　因为脉动项的引入，方程的项数过于多，在湍流方程的推导时经常会使用如下的表达方式，称之为爱因斯坦约定求和法则。即当出现哑标时，则意为求和。

$A_{i} B_{i}=A_{j} B_{j}=A_{1} B_{1}+A_{2} B_{2}+A_{3} B_{3}\\$

登录/注册后可看大图

　　上式的i和j都为哑标，下式的i是自由标，而j是哑标。

$A_{i j} B_{j}=A_{i 1} B_{1}+A_{i 2} B_{2}+A_{i 3} B_{3}\\$

登录/注册后可看大图

　　结合上面的方程组的第一个方程来说，展开就很容易看出来是描述的连续方程，其中

$u_{2}$

登录/注册后可看大图

就是速度

$v$

登录/注册后可看大图

，

$x_{2}$

登录/注册后可看大图

就是

$y$

登录/注册后可看大图

，以此类推。

$\frac{\partial u_{i}}{\partial x_{i}}=\frac{\partial u_{1}}{\partial x_{1}}+\frac{\partial u_{2}}{\partial x_{2}}+\frac{\partial u_{3}}{\partial x_{3}}\\$

登录/注册后可看大图

　　这两个方程左右取平均再联立，即可得到雷诺方程。

$\rho \frac{\partial \bar{u}{j}}{\partial t}+\rho \bar{u}{i} \frac{\partial \bar{u}{j}}{\partial x{i}}=-\frac{\partial \bar{p}}{\partial x_{j}}+\mu \frac{\partial^{2} \bar{u}{j}}{\partial x{i} \partial x_{i}}+\frac{\partial}{\partial x_{i}}\left(-\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}\right)\\$

登录/注册后可看大图

　　中间的推导过程略多，暂时跳过，如有需要可以查阅相关书籍。但最重要的是，需要知道雷诺应力项（下式），是来源于NS方程中的非线性项，也就是等号左边偏导数后面那一项，通常写成

$u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z}$

登录/注册后可看大图

（以x方向运动方程的非线性项为例）。

$\frac{\partial}{\partial x_{i}}\left(-\rho \overline{u_{i}^{\prime} u_{j}^{\prime}}\right)\\$

登录/注册后可看大图

　　至此，考虑湍流的流动的控制方程已经得到了。但是这么做又引入了一个新的问题，就是方程不封闭了。通俗的讲，就是又引入了一大堆新的变量，而方程数少于未知量的数量。为了使得方程组封闭，就需要通过更多的关系来建立方程，确立待求解变量直接的关系。
3.混合长度模式（代数涡粘模式）　　在力学中，应该都接触过牛顿黏性定律。假设如下图所示的场景。

登录/注册后可看大图

　　摩擦力可以表达成下式，该理论认为摩擦力和速度的剪切是正比关系，而这一比例系数为

$\mu$

登录/注册后可看大图

。

$F=\mu A \frac{\mathrm{du}}{\mathrm{d} y} \\$

登录/注册后可看大图

　　接下来会略微有些抽象，还是以上述平行剪切流为例。考虑

$y$

登录/注册后可看大图

层和

$y+l$

登录/注册后可看大图

层这两层，因为脉动，一个流体质点从

$y+l$

登录/注册后可看大图

层运动到

$y$

登录/注册后可看大图

层，在这个过程中流体质点没有任何碰撞（在分子间隙以内），即速度保持不变。为了保证该脉动是在分子间隙之内，也就是为了保证流体速度不变，这个

$l$

登录/注册后可看大图

就是分子平均自由程的长度尺度。该质点到达

$y$

登录/注册后可看大图

层时，引起脉动。通过这一思想，可以分析平均脉动量的数量级。

$\overline{\left|u^{\prime}\right|}=\bar{u}(y+l)-\bar{u}(y) \sim l \frac{d \bar{u}}{d y}\\$

登录/注册后可看大图

　　可以很合理的假设

$u^{\prime}$

登录/注册后可看大图

和

$v^{\prime}$

登录/注册后可看大图

为相同数量级。因此很容易推导下式。

$-\rho \overline{u^{\prime} v^{\prime}} \sim \rho\left(l \frac{d \bar{u}}{d y}\right)^{2}\\$

登录/注册后可看大图

　　这样就可以得到下式。其中

$l_{m}$

登录/注册后可看大图

被被称为混合长度，与湍流结构有关，具体数值可以通过实验等方式测得，是个经验性质的参数。

$\mu_{t}$

登录/注册后可看大图

被称为涡粘系数。

$\begin{array}{c} -\rho \overline{u^{\prime} v^{\prime}}=\mu_{t} \frac{d \bar{u}}{d y} \\ \mu_{t}=\rho l_{m}^{2}\left|\frac{d \bar{u}}{d y}\right| \end{array} \\$

登录/注册后可看大图

　　以上是二维平板剪切流的表达形式。在三维时，雷诺应力就变成了如下形式。

$-\rho u_{i}^{\prime} u_{j}^{\prime}=\mu_{t}\left(\frac{\partial \bar{u}{i}}{\partial x{j}}+\frac{\partial \bar{u}{j}}{\partial x{i}}\right)-\frac{2}{3} k \delta_{i j}\\$

登录/注册后可看大图

　　对下式雷诺应力这九项都通过这样的方式去建立关系，方程组就封闭了。然而，混合长度的经验性，给模式的求解带来了不确定性。具体在实验中，对三维流动效果不如二维好。

$\tau^{t}=\left[\begin{array}{ccc} -\rho \overline{u^{\prime 2}} & -\rho \overline{u^{\prime} v^{\prime}} & -\rho \overline{u^{\prime} w^{\prime}} \\ -\rho \overline{u^{\prime} v^{\prime}} & -\rho \overline{v^{\prime 2}} & -\rho \overline{v^{\prime} w^{\prime}} \\ -\rho \overline{u^{\prime} w^{\prime}} & -\rho \overline{v^{\prime} w^{\prime}} & -\rho \overline{w^{\prime 2}} \end{array}\right]\\$

登录/注册后可看大图

4.k-ε模型　　k-ε模型对涡粘系数的处理方式如下，是由上述模型发展而来的，并且考虑了历史效应。其中

$k$

登录/注册后可看大图

为湍动能，

$\mathcal{E}$

登录/注册后可看大图

为湍动能耗散率。因为涡粘系数用这种方式表达，又引出了两个变量，因此为了保证方程组的封闭，还需要两个方程。

$\begin{aligned} -\overline{u_{i}^{\prime} u_{j}^{\prime}} &=v_{t}\left(\frac{\partial \bar{u}{i}}{\partial x{j}}+\frac{\partial \bar{u}{j}}{\partial x{i}}\right)-\frac{2}{3} k \delta_{i j} \\ v_{t} &=C_{\mu} \frac{k^{2}}{\varepsilon} \end{aligned} \\$

登录/注册后可看大图

　　为了封闭

$k$

登录/注册后可看大图

，需要湍动能方程。为了封闭

$\mathcal{E}$

登录/注册后可看大图

，需要湍动能耗散率的输运方程，分别如下所示。若想了解方程中这些系数的确定，可以查阅文末参考资料。在此不过多赘述。

$\begin{array}{c} \frac{\partial k}{\partial t}+\bar{u}_{i} \frac{\partial k}{\partial x_{i}}=2 v_{t} \bar{s}_{i j} \frac{\partial \bar{u}_{i}}{\partial x_{j}}-\frac{\partial}{\partial x_{i}}\left[\left(v+\frac{v_{t}}{\sigma_{k}}\right) \frac{\partial k}{\partial x_{i}}\right]-\varepsilon \\ \frac{\partial \varepsilon}{\partial t}+\bar{u}_{i} \frac{\partial \varepsilon}{\partial x_{i}}=C_{\varepsilon 1} \frac{\varepsilon}{k}\left(2 v_{t} \bar{s}_{i j} \frac{\partial \bar{u}_{i}}{\partial x_{j}}\right)-\frac{\partial}{\partial x_{i}}\left[\left(v+\frac{v_{t}}{\sigma_{\varepsilon}}\right) \frac{\partial \varepsilon}{\partial x_{i}}\right]-C_{\varepsilon 2} \frac{\varepsilon^{2}}{k} \end{array}\\$

登录/注册后可看大图

　　以上为标准的k-ε模型，除此之外，k-ε模型还有种种变种。想要进一步了解可以查阅文末参考资料。
5.POM模式的Mellor and Yamada方案　　现在，再来回看POM的几个基本控制方程，可以发现运动方程和温盐方程都带着湍流项。因为湍流具有混合性，倾向于使流体变得更加混合。因此，这些项在运动方程中，体现的更多的是粘性作用。在温盐方程中，起的是扩散作用。

$\begin{array}{l} \frac{\partial U}{\partial t}+U \frac{\partial U}{\partial x}+V \frac{\partial U}{\partial y}+W \frac{\partial U}{\partial z}-f V=-\frac{1}{\rho_{0}} \frac{\partial p}{\partial x}+\frac{\partial}{\partial z}\left(K_{M} \frac{\partial U}{\partial z}\right)+F_{x} \\ \frac{\partial V}{\partial t}+U \frac{\partial V}{\partial x}+V \frac{\partial V}{\partial y}+W \frac{\partial V}{\partial z}+f U=-\frac{1}{\rho_{0}} \frac{\partial p}{\partial y}+\frac{\partial}{\partial z}\left(K_{M} \frac{\partial V}{\partial z}\right)+F_{y} \\ \frac{\partial T}{\partial t}+U \frac{\partial T}{\partial x}+V \frac{\partial T}{\partial y}+W \frac{\partial T}{\partial z}=\frac{\partial}{\partial z}\left(K_{H} \frac{\partial T}{\partial z}\right)+F_{T} \\ \frac{\partial S}{\partial t}+U \frac{\partial S}{\partial x}+V \frac{\partial S}{\partial y}+W \frac{\partial S}{\partial z}=\frac{\partial}{\partial z}\left(K_{H} \frac{\partial S}{\partial z}\right)+F_{S} \end{array}\\$

登录/注册后可看大图

　　而模式中对湍流的参数化方案，就是要确定这些湍流项前面的系数。其主要方程如下，具体系数的推导较为复杂，若想要深入了解POM模式所用的湍流闭合方案可查阅论文Development of a Turbulence Closure Model for Geophysical Fluid Problems。

$\begin{array}{l} \frac{\partial q^{2}}{\partial t}+U \frac{\partial q^{2}}{\partial x}+V \frac{\partial q^{2}}{\partial y}+W \frac{\partial q^{2}}{\partial z}=\frac{\partial}{\partial z}\left(K_{q} \frac{\partial q^{2}}{\partial z}\right)+2 K_{M}\left[\left(\frac{\partial U}{\partial Z}\right)^{2}+\left(\frac{\partial V}{\partial Z}\right)^{2}\right]+\frac{2 g}{\rho_{0}} K_{H} \frac{\partial \rho}{\partial z}-\frac{2 q^{3}}{B_{1} l}+F_{q} \\ \frac{\partial q^{2} l}{\partial t}+U \frac{\partial q^{2} l}{\partial x}+V \frac{\partial q^{2} l}{\partial y}+W \frac{\partial q^{2} l}{\partial z}=\frac{\partial}{\partial z}\left(K_{q} \frac{\partial q^{2} l}{\partial z}\right)+l E_{1} K_{M}\left[\left(\frac{\partial U}{\partial Z}\right)^{2}+\left(\frac{\partial V}{\partial Z}\right)^{2}\right]+\frac{l E_{1} g}{\rho_{0}} K_{H} \frac{\partial \rho}{\partial z}-\frac{q^{3}}{B_{1}} \widetilde{W}+F_{l} \\ K_{M} \equiv l q S_{M} \\ K_{H} \equiv l q S_{H} \\ K_{q} \equiv l q S_{q} \end{array}\\$

登录/注册后可看大图

6.总结　　以上所讲的湍流模式都有一个共同点，就是试图用统计平均的方式封闭雷诺平均方程，使得平均运动可解。如何对方程进行闭合，通常被称为湍流闭合方案。而如何确定湍流闭合所带来的系数，通常被称为湍流参数化方案。
　　除了建立湍流统计模式，研究湍流时还会使用模拟的方式，比如湍流的直接数值模拟（DNS）和大涡模拟（LES）。
版权声明　　本文创作的初衷是用于帮助数值模式的学习者。欢迎转载，转载请私信并注明作者和出处，请勿用于任何商业用途。
参考书目湍流理论与模拟（第2版）. 张兆顺、崔桂⾹、许春晓编著. 清华⼤学出版社，2005.9.
A Manual for POM and GOMO. Xiaomeng Huang, Xing Huang.
Development of a turbulence closure model for geophysical fluid problems. Mellor G L, Yamada T. Reviews of Geophysics, 1982.

数值海洋与大气模式（五）：湍流模式与参数化方案

相关帖子

您看了很久哦，登陆下吧！