湍流模型（2）——雷诺平均方程

  本文内容主要是不可压的雷诺平均方程（Reynolds-averaged Navier-Stokes equations，RANS）的推导。在推导雷诺方程之前先来总结一些关于时间平均的运算规律。

1. 关于时间平均的运算规律证明

  对于瞬时量$\varphi =\Phi +{\varphi }^{\prime }$和 $\psi =\Psi +{\psi }^{\prime }$，$\Phi$和$\Psi$已经是与时间无关的量，所以根据时间平均的定义有
$$\overline{{\varphi }^{\prime }}=\overline{{\psi }^{\prime }}=0\text{\hspace{1em}(1.1)}$$
$$\overline{\Phi }=\frac{1}{\Delta t}{\int }_0^{\Delta t}\Phi dt=\Phi \frac{1}{\Delta t}{\int }_0^{\Delta t}dt=\Phi \text{\hspace{1em}(1.2)}$$
$$\begin{array}{rl}\overline{\frac{\partial \varphi }{\partial s}}& =\frac{1}{\Delta t}{\int }_0^{\Delta t}\frac{\partial \varphi }{\partial s}dt=\frac{\partial }{\partial s}\left(\frac{1}{\Delta t}{\int }_0^{\Delta t}\varphi dt\right)=\frac{\partial \bar{\varphi }}{\partial s}=\frac{\partial \Phi }{\partial s}\end{array}\text{\hspace{1em}(1.3)}$$
$$\begin{array}{rl}\overline{\int \varphi ds}& =\frac{1}{\Delta t}{\int }_0^{\Delta t}\left(\int \varphi ds\right)dt=\int \left(\frac{1}{\Delta t}{\int }_0^{\Delta t}\varphi dt\right)ds=\int \bar{\varphi }ds=\int \Phi ds\end{array}\text{\hspace{1em}(1.4)}$$
$$\begin{array}{rl}\overline{\varphi +\psi }& =\frac{1}{\Delta t}{\int }_0^{\Delta t}(\varphi +\psi )dt=\frac{1}{\Delta t}{\int }_0^{\Delta t}\varphi dt+\frac{1}{\Delta t}{\int }_0^{\Delta t}\psi dt=\Phi +\Psi \end{array}\text{\hspace{1em}(1.5)}$$
$$\overline{\varphi \Psi }=\frac{1}{\Delta t}{\int }_0^{\Delta t}(\varphi \Psi )dt=\Psi \frac{1}{\Delta t}{\int }_0^{\Delta t}\varphi dt=\Phi \Psi \text{\hspace{1em}(1.6)}$$
$$\overline{{\varphi }^{\prime }\Psi }=\frac{1}{\Delta t}{\int }_0^{\Delta t}({\varphi }^{\prime }\Psi )dt=\Psi \frac{1}{\Delta t}{\int }_0^{\Delta t}{\varphi }^{\prime }dt=\Psi \overline{{\varphi }^{\prime }}=0\text{\hspace{1em}(1.7)}$$
$$\begin{array}{rl}\overline{\varphi \psi }=\overline{(\Phi +{\varphi }^{\prime })(\Psi +{\psi }^{\prime })}& =\overline{\Phi \Psi +\Phi {\psi }^{\prime }+{\varphi }^{\prime }\Psi +{\varphi }^{\prime }{\psi }^{\prime }}\\ & =\overline{\Phi \Psi }+\overline{\Phi {\psi }^{\prime }}+\overline{{\varphi }^{\prime }\Psi }+\overline{{\varphi }^{\prime }{\psi }^{\prime }}\\ & =\overline{\Phi \Psi }+0+0+\overline{{\varphi }^{\prime }{\psi }^{\prime }}\\ & =\Phi \Psi +\overline{{\varphi }^{\prime }{\psi }^{\prime }}\end{array}\text{\hspace{1em}(1.8)}$$
  以上均为标量的运算，矢量的散度和梯度计算可以通过上述公式扩展而来。
  本文中均以粗体字母代表矢量，如速度矢量$\mathbf{u}$和其他任意矢量$\mathbf{a}$等，以$\nabla \cdot \mathbf{u}$表示矢量$\mathbf{u}$的散度，以$\nabla \varphi$表示标量$\varphi$的梯度。根据雷诺分解，瞬时速度矢量$\mathbf{u}=\mathbf{U}+{\mathbf{u}}^{\prime }$和瞬时标量$\varphi =\Phi +{\varphi }^{\prime }$有如下运算规律：
$$\mathbf{u}=(u,v,w{)}^T$$

$$\bar{\mathbf{u}}=\mathbf{U}=(\bar{u},\bar{v},\bar{w},)=(U,V,W{)}^T$$

$${\mathbf{u}}^{\prime }=(u^{\prime },v^{\prime },w^{\prime }{)}^T$$
$$\begin{array}{rl}\overline{\nabla \cdot \mathbf{u}}& =\overline{\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}}=\frac{\partial \bar{u}}{\partial x}+\frac{\partial \bar{v}}{\partial y}+\frac{\partial \bar{w}}{\partial z}\\ & =\frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}+\frac{\partial W}{\partial z}=\nabla \cdot \mathbf{U}=\nabla \cdot \bar{\mathbf{u}}\end{array}\text{\hspace{1em}(1.9)}$$
$$\overline{\varphi \mathbf{u}}=(\overline{\varphi u},\overline{\varphi v},\overline{\varphi w}{)}^T=\left[\begin{array}{c}\Phi U+\overline{{\varphi }^{\prime }{u}^{\prime }}\\ \\ \Phi V+\overline{{\varphi }^{\prime }{v}^{\prime }}\\ \\ \Phi W+\overline{{\varphi }^{\prime }{w}^{\prime }}\end{array}\right]=\Phi \mathbf{U}+\overline{{\varphi }^{\prime }{\mathbf{u}}^{\prime }}\text{\hspace{1em}(1.10)}$$
$$\overline{\nabla \cdot (\varphi \mathbf{u})}=\nabla \cdot \overline{\varphi \mathbf{u}}=\nabla \cdot (\Phi \mathbf{U})+\nabla \cdot (\overline{{\varphi }^{\prime }{\mathbf{u}}^{\prime }})\text{\hspace{1em}(1.11)}$$
$$\overline{\nabla \varphi }=\overline{{\left(\frac{\partial \varphi }{\partial x},\frac{\partial \varphi }{\partial y},\frac{\partial \varphi }{\partial z}\right)}^T}={\left(\frac{\partial \bar{\varphi }}{\partial x},\frac{\partial \bar{\varphi }}{\partial x},\frac{\partial \bar{\varphi }}{\partial x}\right)}^T={\left(\frac{\partial \Phi }{\partial x},\frac{\partial \Phi }{\partial x},\frac{\partial \Phi }{\partial x}\right)}^T=\nabla \Phi \text{\hspace{1em}(1.12)}$$
$$\overline{\nabla \cdot \nabla \varphi }=\nabla \cdot \overline{\nabla \varphi }=\nabla \cdot \nabla \Phi \text{\hspace{1em}(1.13)}$$

  $\nabla \cdot \nabla$也可写作${\nabla }^2$,即拉普拉斯算子。

2. 雷诺平均方程的推导过程

  有了上述准备，推导雷诺平均方程就容易多了，这里考虑三维笛卡尔坐标系下瞬态不可压的连续性方程和N-S方程，忽略体积力，方程如下：
$$\nabla \cdot \mathbf{u}=0\text{\hspace{1em}(2.1)}$$

$$\frac{\partial u}{\partial t}+\nabla \cdot (u\mathbf{u})=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu {\nabla }^{2}u\text{\hspace{1em}(2.2a)}$$

$$\frac{\partial v}{\partial t}+\nabla \cdot (v\mathbf{u})=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\nu {\nabla }^{2}v\text{\hspace{1em}(2.2b)}$$

$$\frac{\partial w}{\partial t}+\nabla \cdot (w\mathbf{u})=-\frac{1}{\rho }\frac{\partial p}{\partial z}+\nu {\nabla }^{2}w\text{\hspace{1em}(2.2c)}$$
  雷诺平均方程的思想就是先将方程中的流动变量分解成平均量和脉动量之和，然后对方程两边同时取时间平均。根据雷诺分解，方程中的各瞬时量分解如下：
$$\mathbf{u}=\mathbf{U}+{\mathbf{u}}^{\prime }u=U+u^{\prime }v=V+v^{\prime }w=W+w^{\prime }p=P+p^{\prime }\text{\hspace{1em}(2.3)}$$
将$(2.3)$各式代入到式$(2.1)$~$(2.2c)$中，然后对方程取时间平均，再根据第1节中的运算规律对方程进行整理。

先来推导方程$(2.1)$，根据式$(1.9)$
$$\overline{\nabla \cdot \mathbf{u}}=\nabla \cdot \overline{\mathbf{u}}=\nabla \cdot \mathbf{U}$$
  于是得到平均流的连续性方程：
$$\nabla \cdot \mathbf{U}=0\text{\hspace{1em}(2.4)}$$
  其中$\mathbf{U}$是平均流速度矢量，写成分量形式为
$$\frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}+\frac{\partial W}{\partial z}=0$$

推导方程$(2.2a)$：
  方程左边第一项，根据式$(1.3)$
$$\overline{\frac{\partial u}{\partial t}}=\frac{\partial U}{\partial t}$$
  左边第二项，根据式$(1.11)$
$$\overline{\nabla \cdot (u\mathbf{u})}=\nabla \cdot (U\mathbf{U})+\nabla \cdot (\overline{u^{\prime }{\mathbf{u}}^{\prime }})$$
  方程右边第一项，
$$\overline{-\frac{1}{\rho }\frac{\partial p}{\partial x}}=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x}=-\frac{1}{p}\frac{\partial P}{\partial x}$$
  右边第二项，根据式$(1.13)$
$$\overline{\nu {\nabla }^{2}u}=\nu {\nabla }^{2}U$$
  综上，对方程$(2.2a)$取时间平均得到$x$分量的时均动量方程：

$$\begin{array}{r}\frac{\partial U}{\partial t}+\nabla \cdot (U\mathbf{U})+\nabla \cdot (\overline{u^{\prime }{\mathbf{u}}^{\prime }})=-\frac{1}{\rho }\frac{\partial P}{\partial x}+\nu {\nabla }^{2}U\\ (1)(2)(3)(4)(5)\end{array}\text{\hspace{1em}(2.5a)}$$

同样的推导过程可以得到$y$和$z$分量的时均动量方程：
$$\frac{\partial V}{\partial t}+\nabla \cdot (V\mathbf{U})+\nabla \cdot (\overline{v^{\prime }{\mathbf{u}}^{\prime }})=-\frac{1}{\rho }\frac{\partial P}{\partial y}+\nu {\nabla }^{2}V\text{\hspace{1em}(2.5b)}$$

$$\frac{\partial W}{\partial t}+\nabla \cdot (W\mathbf{U})+\nabla \cdot (\overline{w^{\prime }{\mathbf{u}}^{\prime }})=-\frac{1}{\rho }\frac{\partial P}{\partial z}+\nu {\nabla }^{2}W\text{\hspace{1em}(2.5c)}$$
以方程$(2.5a)$为例，可以看到第（2）、（3）项都是来自原方程的对流项，第（3）项是湍流脉动量的乘积，是时间平均后多出来的项，其代表的物理含义是湍流涡的对流所引起动量输运，这就类似于流体在粘性剪应力的作用下又附加了湍流剪应力。所以习惯上，这一项会移到方程的右边和粘性项放在一起，以反应其物理属性。为了展示湍流应力的结构，把第（3）项中的脉动速度矢量写出分量形式，并移到等号右边，则有
$$\begin{array}{rl}\frac{\partial U}{\partial t}+\nabla \cdot (U\mathbf{U})=& -\frac{1}{\rho }\frac{\partial P}{\partial x}+\nu {\nabla }^{2}U\\ & +\frac{1}{\rho }\left[\frac{\partial (-\rho \overline{u^{\prime 2}})}{\partial x}+\frac{\partial (-\rho \overline{u^{\prime }{v}^{\prime }})}{\partial y}+\frac{\partial (-\rho \overline{u^{\prime }{w}^{\prime }})}{\partial z}\right]\end{array}\text{\hspace{1em}(2.6a)}$$
$$\begin{array}{rl}\frac{\partial V}{\partial t}+\nabla \cdot (V\mathbf{U})=& -\frac{1}{\rho }\frac{\partial P}{\partial y}+\nu {\nabla }^{2}V\\ & +\frac{1}{\rho }\left[\frac{\partial (-\rho \overline{u^{\prime }{v}^{\prime }})}{\partial x}+\frac{\partial (-\rho \overline{v^{\prime 2}})}{\partial y}+\frac{\partial (-\rho \overline{v^{\prime }{w}^{\prime }})}{\partial z}\right]\end{array}\text{\hspace{1em}(2.6b)}$$
$$\begin{array}{rl}\frac{\partial W}{\partial t}+\nabla \cdot (W\mathbf{U})=& -\frac{1}{\rho }\frac{\partial P}{\partial z}+\nu {\nabla }^{2}W\\ & +\frac{1}{\rho }\left[\frac{\partial (-\rho \overline{u^{\prime }{w}^{\prime }})}{\partial x}+\frac{\partial (-\rho \overline{v^{\prime }{w}^{\prime }})}{\partial y}+\frac{\partial (-\rho \overline{w^{\prime 2}})}{\partial z}\right]\end{array}\text{\hspace{1em}(2.6c)}$$
湍流应力一共有6个独立分量，包括3个法向应力分量：
$${\tau }_{xx}=-\rho \overline{u^{\prime 2}}{\tau }_{yy}=-\rho \overline{v^{\prime 2}}{\tau }_{zz}=-\rho \overline{w^{\prime 2}}$$
和3个切应力分量：
$${\tau }_{xy}={\tau }_{yx}=-\rho \overline{u^{\prime }{v}^{\prime }}{\tau }_{xz}={\tau }_{zx}=-\rho \overline{u^{\prime }{w}^{\prime }}{\tau }_{zy}={\tau }_{yz}=-\rho \overline{w^{\prime }{v}^{\prime }}$$
这些湍流应力称为雷诺应力，方程式$(2.4)$和$(2.6a)$~$(2.6c)$称为雷诺平均N-S方程(RANS)。

参考资料：

1. Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.