动态 SGS 模型

输运方程

可压缩流动的N-S方程如下：

$\begin{array}{l} \frac{\partial \rho }{\partial t} + \frac{\partial \rho {u_i}}{\partial {x_i}} = 0\\ \frac{\partial \rho u_i} {\partial t} + \frac{\partial \rho u_iu_j} {\partial x_j} = - \frac{\partial p}{\partial x_i} + \frac{\partial \sigma _{ij}} {\partial x_j} + \rho f_i \end{array}$

其中，$\sigma _{ij} = 2\mu S_{ij} - \frac{2}{3}\mu (\frac{\partial u_k} {\partial x_k}) \delta _{ij}\ ,S_{ij} = \frac{1}{2}(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i})$。

对N-S方程进行雷诺平均，并引入Favre平均 ：$\widetilde {f_i} = \frac {\overline {\rho f_i} } {\bar \rho }$，可以得到：

$\begin{array}{l} \frac{\partial \bar \rho }{\partial t} + \frac{\partial \bar \rho \widetilde {u_i}}{\partial x_i} = 0\\ \frac{\partial }{\partial t}(\bar \rho \widetilde {u_i}) + \frac{\partial \bar \rho \widetilde {u_i}\widetilde {u_j}}{\partial x_j} = - \frac{\partial \bar p}{\partial x_i} + \frac{\partial \widetilde {\sigma _{ij}}}{\partial x_j} + \bar \rho \widetilde {f_i} + \frac{\partial }{\partial x_j}( - \bar \rho \widetilde {u'_iu'_j}) \end{array}$

其中，

$\widetilde {\sigma _{ij}} = 2\mu \widetilde {S_{ij}} - \frac{2}{3}\mu \widetilde {S_{kk}}{\delta _{ij}}{\rm{ , }}\widetilde {S_{ij}} = \frac{1}{2}(\frac{\partial \widetilde {u_i}}{\partial {x_j}} + \frac{\partial \widetilde {u_j}}{\partial x_i})$

$\widetilde {u_iu_j} = \widetilde {u_i}\widetilde {u_j} - (\widetilde {u_i}\widetilde {u_j} - \widetilde {u_iu_j})$

LES的动量方程为：

$\frac{\partial }{\partial t}(\bar \rho \widetilde {u_i}) + \frac{\partial \bar \rho \widetilde {u_i}\widetilde {u_j}}{\partial {x_j}} = - \frac{\partial \bar p}{\partial x_i} + \frac{\partial \widetilde {\sigma _{ij}}}{\partial x_j} + \bar \rho \widetilde {f_i} - \frac{\partial }{\partial x_j}\bar \rho (\widetilde {u_iu_j} - \widetilde {u_i}\widetilde {u_j})$

定义亚格子应力为$\tau _{ij} = \bar \rho (\widetilde {u_iu_j} - \widetilde {u_i}\widetilde {u_j})$。

这样，动量方程就变为：

$\frac{\partial }{\partial t}(\bar \rho \widetilde {u_i}) + \frac{\partial \bar \rho \widetilde {u_i}\widetilde {u_j}}{\partial x_j} = - \frac{\partial \bar p}{\partial x_i} + \frac{\partial \widetilde {\sigma _{ij}}}{\partial x_j} + \bar \rho \widetilde {f_i} - \frac{\partial \tau _{ij}}{\partial x_j}$

Smagorinsky模型：

$\tau _{ij}^{Smag} - \frac{1}{3}\tau _{kk}^{Smag}{\delta _{ij}} = - 2{\mu _{SGS}}\widetilde {S_{ij}}$

其中，${\mu _{SGS}} = \bar \rho {({C_s}\bar \Delta )^2}|\widetilde S|$,$|\widetilde S| = \sqrt {2\widetilde {S_{ij}}\widetilde {S_{ij}}}$。

Yoshizawa.1986对其进行了可压缩修正，我们称之为可压 Smagorinsky 模型：

$\begin{array}{l} \tau _{ij} - \frac{1}{3}\tau _{kk}\delta _{ij} = - 2\textcolor{red}{\bar \rho (C_s{\bar \Delta })^2 |\widetilde S|}(\widetilde {S_{ij}} - \frac{1}{3}\widetilde {S_{kk}}\delta _{ij}) = C_s^2\alpha _{ij}\\ \tau _{kk} = C_I 2\bar \rho {\bar \Delta }^2|\widetilde S{|^2} = {C_I}\alpha \end{array}$

其中，${C_s} = 0.16{\rm{ }},{\rm{ }}{C_I} = 0.09{\rm{ }},{\rm{ }}|\widetilde S| = \sqrt {2\widetilde {S_{ij}}\widetilde {S_{ij}}}$。
1991 年，Germano 提出了动态 Smagorinsky 模型，该模型通过局部自相似假设，将两个不同尺度的亚格子应力模型联系起来，得到在计算过程中动态决定模型系数的方法。动态大涡模型的过程中，引入了二次滤波，称为检验滤波（test filter），检验滤波的滤波宽度为$\widehat \Delta$，比一次滤波的滤波宽度要大，一般为$\widehat \Delta = 2\bar \Delta$ 。但是 Germano 提出的动态模型是基于不可压流体的，这里我们采用Moin.1991提出的可压动态 SGS 模型。对一次滤波之后的动量方程进行检验滤波，得到如下方程：

$\frac{\partial }{\partial t}(\widehat {\bar \rho \widetilde {u_i}}) + \frac{\partial \widehat {\bar \rho \widetilde {u_i}\widetilde {u_j}}}{\partial x_j} = - \frac{\partial \widehat {\bar p}}{\partial x_i} + \frac{\partial \widehat {\widetilde {\sigma _{ij}}}}{\partial x_j} + \widehat {\bar \rho \widetilde {f_i}} - \frac{\partial }{\partial x_j}\widehat {\bar \rho (\widetilde {u_iu_j} - \widetilde {u_i}\widetilde {u_j})}$

根据$\bar \rho \widetilde {u_i} = \overline {\rho {u_i}}$ 和$\bar \rho \widetilde {u_iu_j} = \overline {\rho {u_i}{u_j}}$ ，我们得到：

$\frac{\partial } {\partial t} (\widehat{\overline{\rho {u_i}}}) +\frac{\partial } {\partial {x_j}} (\frac{\widehat{\overline{\rho u_i}}\widehat{\overline{\rho u_i}}}{\widehat{\bar\rho}}) = - \frac{\partial \widehat{\bar p}} {\partial x_i} + \frac{\partial \widehat{\widetilde{\sigma_{ij}}}} {\partial x_j} + \widehat{\bar\rho \widetilde{f_i}} - \frac{\partial } {\partial x_j} \widehat{\overline{\rho u_iu_j}} \bcancel{+\frac{\partial }{\partial {x_j}}\widehat{\bar\rho \widetilde{u_i}\widetilde{u_j}}} + \frac{\partial } {\partial x_j} (\frac{\widehat{\overline{\rho u_i}\widehat{\overline{\rho u_j}}}} {\widehat{\bar\rho }}) \bcancel{-\frac{\partial }{\partial x_j}\widehat{\bar{\rho }\widetilde{u_i}\widetilde{u_j}}}$

$\frac{\partial }{\partial t} (\widehat{\overline{\rho {u_i}}})+\frac{\partial }{\partial {x_j}}(\frac{\widehat{\overline{\rho u_i}}\widehat{\overline{\rho {u_i}}}}{\widehat{\bar{\rho }}})=-\frac{\partial \widehat{\bar{p}}}{\partial x_i}+\frac{\partial \widehat{\widetilde{\sigma_{ij}}}}{\partial x_j}+\widehat{\bar{\rho }\widetilde{f_i}}-\frac{\partial T_{ij}}{\partial x_j}$

$T_{ij}= \widehat{\overline{\rho u_iu_j}} - \frac{\widehat{\overline{\rho u_i}}\widehat{\overline{\rho u_j}}} {\widehat{\bar\rho }}$

$T_{ij}$ 是两次滤波之后的亚格子应力。

$\begin{array}{l} L_{ij}=T_{ij}-\widehat {\tau _{ij}}=\widehat{\overline{\rho {u_i}{u_j}}}-\frac{\widehat{\overline{\rho {u_i}}} \widehat{\overline{\rho u_j}}}{\widehat{\bar{\rho }}} -(\widehat{\bar\rho }\widetilde{u_iu_j} -\widehat{\bar\rho } \tilde{u_i}\tilde{u_j})\\ \text{ }=\bcancel{\widehat{\overline{\rho {u_i}u_j}}}-\frac{\widehat{\overline{\rho {u_i}}}\widehat{\overline{\rho u_j}}}{\widehat{\bar{\rho }}}\bcancel{-\widehat{\overline{\rho {u_{i}}{u_j}}}}+\frac{\widehat{\overline{\rho u_i}\overline{\rho {u_j}}}}{\hat{\bar{\rho }}}\\\text{ }=\frac{\widehat{\overline{\rho {u_{i}}}\overline{\rho u_j}}}{\hat{\bar{\rho }}}-\frac{\widehat{\overline{\rho {u_{i}}}}\widehat{\overline{\rho {u_{j}}}}}{\widehat{\bar{\rho }}}\\\text{ }=\widehat{\bar{\rho }\widetilde{u_{i}}\widetilde{u_{j}}}-\frac{\widehat{\bar{\rho }\widetilde{u_i}}\widehat{\bar{\rho }\widetilde{u_j}}}{\widehat{\bar{\rho }}} \end{array}$

这就是所谓的 Leonard 应力，上次称为 Germano 恒等式？

分别对$\widehat{\tau_{ij}}$ 和$T_{ij}$ 采用可压 Smagorinsky 模型：

$\widehat {\tau_{ij}}-\frac{1}{3}\widehat{\tau_{kk}}=C_s^2\widehat{\alpha_{ij}} \\ T_{ij}-\frac{1}{3} {T_{kk}}\delta_{ij}=-C_s^22\widehat{\bar\rho} \widehat \Delta^2|\ddot S|(\ddot S_{ij}-\frac{1}{3}\ddot S_{kk}\delta_{ij})=C_s^2\beta_{ij}\\ where: \ddot f=\frac{\widehat{\overline{\rho f}}}{\widehat{\bar\rho}}\\ \alpha_{ij} =- 2\textcolor{red}{\bar \rho {\bar \Delta }^2 |\widetilde S|}(\widetilde {S_{ij}} - \frac{1}{3}\widetilde {S_{kk}}\delta _{ij})\\ \beta_{ij}=-2\widehat{\bar\rho} \widehat \Delta^2|\ddot S|(\ddot S_{ij}-\frac{1}{3}\ddot S_{kk}\delta_{ij})$

由$L_{ij} = T_{ij} - \widehat {\tau _{ij}}$，我们可以得到：

${L_{ij}} - \frac{1}{3}\widehat {L_{kk}}{\delta _{ij}} = C_s^2{\beta _{ij}} - C_s^2\widehat {\alpha _{ij}} = C_s^2{M_{ij}}\\ where: M_{ij}=\beta_{ij}-\widehat {\alpha _{ij}}$

方程组是超定的，因此上式应该看作是在某种平均意义上成立。

Martín.2000采用Lilly.1992提出的最小二乘法，使式（15）的误差最小化，同时采用对均匀方向平均的形式：

$C_s^2 = \frac{\left\langle {L_{ij}{M_{ij}}} \right\rangle }{\left\langle {M_{kl}{M_{kl}}} \right\rangle }{\text{ }},{\text{ }}{C_I} = \frac{\left\langle {L_{kk}} \right\rangle }{\left\langle {\beta - \widehat \alpha } \right\rangle }{\text{ }},{\text{ }}\beta = 2{\widehat \Delta ^2}\widehat {\overline \rho }|\ddot S|^2$

这就是可压 Smagorinsky 模型中需要的。

$\begin{array}{l} \tau _{ij} - \frac{1}{3}\tau _{kk}\delta _{ij} = - 2\textcolor{red}{\bar \rho (C_s{\bar \Delta })^2 |\widetilde S|}(\widetilde {S_{ij}} - \frac{1}{3}\widetilde {S_{kk}}\delta _{ij}) = C_s^2\alpha _{ij}\\ \tau _{kk} = C_I 2\bar \rho {\bar \Delta }^2|\widetilde S{|^2} = {C_I}\alpha \end{array}$

滤波宽度比$\widehat \Delta /\overline \Delta$ 是动态模型中唯一需要给定的参数。