Smagorinsky模型系数问题

[李东岳]

针对你的case2与case3，你看看nut以及k的结果，看看有区别么？

把下面的放在controlDict里面可以输出$k_{sgs}$

cacheTemporaryObjects
(
    k
);

functions
{
    #includeFunc writeObjects("k")
}

我觉得这些东西，可以写个英文让大家知道一下了。投个一般的期刊当做一个note都可以。目前我是没发现有人研究这个。并且我认为研究很有意义。 主要是看这东西之前是不是确实没人玩过。如果没人搞过，那研究一下绝对是发现。写成英文的目的不是为了发文章，而是真正的让大家知道一下这个事情：1）两种Smagorinsky的植入区别与对比，2）第二种植入的模型系数问题。我是没发现有人研究这个。

[coolhhh]

@李东岳

1. nut结果对比

nut结果中，case1~case3数值范围都比较接近，但是case4理论上应该与前几个算例结果近似才对，但由于比例系数$1/(C_k^2)=146407$过大，nut计算结果实际是偏大，但数量级都是1e-4

2. tmp<k>结果对比

case3的系数$1/(C_k^2)=1$最小，以这个工况为基准。case3的tmp<k>最大值为0.0035319。
case2的系数$1/(C_k^2)=1419$，理论上case2的最大值应该为0.0035319*1419=5.0117，实际计算的最大值5.0779，与理论结果基本一致。
case4的系数$1/(C_k^2)=146407$，数值特别大。理论上case4的最大值应该为0.0035319*146407=517.096，但实际计算的最大值是42185，比理论结果偏大，这可能是结果显示的最终时刻0.1s的结果，与前面工况的速度场有点不同，还有计算的数值误差等其他原因导致。

[李东岳]

case4的nut特别高，case1-3都可以预测类似的nut所以流动行为差不多。因为给了不同的Ck，所有k区别非常大。这是不是说明Smagorinsky预测的nut是合理的，但是预测的k存在任意性，也侧方面的表示了OpenFOAM没有输出k的原因？

这些k肯定有一种是趋向于合理的。但具体哪个合理目前还不清楚。两种方法的产生与耗散是否相等也未知。我觉得第二种方法的产生与耗散应该是一样的，因为本身那个二次方程就是从P=D推出来的。但是第一种可能够呛。

OpenFOAM原始的Ck、Ce的值是哪里搞出来的你知道么？ @coolhhh

[coolhhh]

@李东岳

1. 李老师，我也没具体看到哪个文章用了这两个默认系数，但搜到一个帖子有个观点，帖子中的超链接无法查看：OpenFOAM大涡模拟湍流模型之Smagorinsky模型代码详解

2. 无意中看到一个日本CFD网站，看起来$k_{sgs}$与压力计算值有关。
   

然后看Smagorinsky的推导，这个推导过程也类似的出现$Ck$和$C_e$参数。疑问依然还是第一种方法是如何计算$k_{sgs}$，如果这个问题知道了，那应该能推断出第一种方法与第二种方法的系数$Ck$和$C_e$的关系。

[李东岳]

OpenFOAM原始的Ck、Ce的值是哪里搞出来的你知道么？

我找了几个原始文献（主要是Lilly1967，DearDorff1970，以及Moeng and Wyngaard1988），里面讨论了这个Ck Ce的东西，我没太细看。不过看起来这个值是推导出来的，要满足。不过Cs的值都不一样。我再看看。

[李东岳]

最近没细看这个玩别的去了。我参考Spectral analysis of large eddy simulations of the convective boundary layer昨天晚上推了个八九不离十。这个ck ce是符合湍流能谱的一个常数。这是前人已经做过的工作。没法写英文了。我以为别人没玩过这个。哈哈

[lyc]

@coolhhh 太细致了！！！

[coolhhh]

@李东岳 李老师，最近把Smagorinsky湍流模型两种方式植入方式再梳理了下，下面只讨论不可压缩湍流情况，重新完整总结如下（由于有的公式latex显示异常，就采用截图形式）：

1. OpenFOAM中LES求解方程形式

根据 【OpenFOAMv2012】Smagorinskyの解説 | mtk_birdman's blog
，对于不可压缩湍流，OpenFOAM中求解的方程如下

其中，
$$\begin{array}{}\text{(3)}& {\nu }_{eff}=\nu +{\nu }_{sgs}\end{array}$$

$$\begin{array}{}\text{(4)}& \overline{P}=\overline{p}+\frac{2}{3}\rho k_{sgs}\end{array}$$

在OpenFOAM中，不可压缩湍流的 $\rho =1$

---

【讨论1】

• OpenFOAM中看起来是将 $\overline{P}=\overline{p}+\frac{2}{3}\rho k_{sgs}$作为整体计算，最后的压力计算结果并没有减去 $\frac{2}{3}\rho k_{sgs}$。因此OpenFOAM中LES计算得到的压力并不是真实的压力，其偏差程度取决于 $k_{sgs}$的大小。
• 因此，OpenFOAM中只要不同亚格子模型得到的 ${\nu }_{sgs}$的值相同，那么其压力结果理论上也应相同，因为 $k_{sgs}$不参与最终的压力计算。

类似的问题在RANS模型也有，OpenFOAM也是把k忽略掉了，详见下面两个帖子的讨论
OpenFOAM 不可压缩湍流模型的 divDevReff 函数 | Giskard's CFD Learning Tricks
Calculating divDevReff - Page 2 -- CFD Online Discussion Forums
上面两个帖子针对OpenFOAM中RANS模型没有减去k有2个观点：
第1个观点是k对于压力的计算结果影响很小因此可以省去。
第2个观点是将k放到了压力项中一起计算，但OpenFOAM最后没有把压力减去k，究竟k对于压力结果影响大小是未知的。

---

2. 符号定义与2个假设

解析应变率张量（resolved-scale strain rate tensor）：
$$\begin{array}{}\text{(5)}& \overline{\mathbf{D}}=\frac{1}{2}(\mathrm{\nabla }\overline{\mathbf{U}}+\mathrm{\nabla }{\overline{\mathbf{U}}}^{\mathrm{T}})\end{array}$$

亚格子应力张量（subgrid scale stress tensor）：

其中 $k_{sgs}$定义为
$$\begin{array}{}\text{(7)}& k_{sgs}=\frac{1}{2}\mathrm{tr}\left({\mathit{\tau }}_{sgs}\right)\end{array}$$

---

根据Smagorinsky SGS model in OpenFOAM | CFD WITH A MISSION，有两个假设：

• Assumption 1: Eddy viscosity approximation，并根据式(7)，可得到
  

• Assumption 2: Local equilibrium（assumption of the balance between the subgrid scale energy production and dissipation）
  $$\begin{array}{}\text{(9)}& P=\epsilon \iff \overline{\mathbf{D}}:{\tau }_{\mathbf{s}\mathbf{g}\mathbf{s}}+C_{\epsilon }\frac{k_{\mathrm{s}\mathrm{g}\mathrm{s}}^{1.5}}{\mathrm{\Delta }}=0\end{array}$$

其中局部平衡假设公式参考资料：

• 《Sullivan, P.P., J.C. McWilliams and C. Moeng, A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows. Boundary-Layer Meteorology, 1994. 71: p. 247-276.》
  
• 【OpenFOAMv2012】Smagorinskyの解説 | mtk_birdman's blog
  

3. 第2种方法（OpenFOAM的植入方式）讨论

若直接根据 Local equilibrium 假设求解式(9)，存在2个未知数: ${\nu }_{sgs}$,$k_{sgs}$，因此还需要构造 ${\nu }_{sgs}$和 $k_{sgs}$的关系式：

$$\begin{array}{}\text{(10)}& {\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}=C_k\mathrm{\Delta }\sqrt{k_{\mathrm{s}\mathrm{g}\mathrm{s}}}\end{array}$$

事实上，式(10)额外增加了未知数$C_k$，但这最后是作为经验参数直接输入，需要根据不同的模拟流场类型改变。

式(10)来源可参考《Sullivan, P.P., J.C. McWilliams and C. Moeng, A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows. Boundary-Layer Meteorology, 1994. 71: p. 247-276.》

根据Smagorinsky SGS model in OpenFOAM | CFD WITH A MISSION，接着将式(8)和(10)代入式(9)，就可求解 $k_{sgs}$，再将 $k_{sgs}$代入式(10)可求得 ${\nu }_{sgs}$。

式中，

对于不可压缩湍流，式(12)化简为：

式中，
$$\begin{array}{}\text{(14)}& \left|\overline{\mathbf{D}}\right|=\sqrt{2\overline{\mathbf{D}}:\overline{\mathbf{D}}}\end{array}$$

把式(13)代入式(11)，得到
$$\begin{array}{}\text{(15)}& k_{\mathrm{s}\mathrm{g}\mathrm{s}}=\frac{c}{a}=\frac{C_k{\mathrm{\Delta }}^2{\left|\overline{\mathbf{D}}\right|}^2}{C_{\epsilon }}\end{array}$$

$$\begin{array}{}\text{(16)}& {\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}=C_k\mathrm{\Delta }\sqrt{k_{\mathrm{s}\mathrm{g}\mathrm{s}}}=C_k\sqrt{\frac{C_k}{C_{\epsilon }}}{\mathrm{\Delta }}^2\left|\overline{\mathbf{D}}\right|\end{array}$$

式(16)对比第1种方法式(18)，可得到关系式：
$$\begin{array}{}\text{(17)}& C_s={\left(C_k\sqrt{\frac{C_k}{C_{\epsilon }}}\right)}^{1/2}\end{array}$$

OpenFOAM中，默认系数设置为 $C_{\epsilon }=1.048,C_k=0.094\Rightarrow C_s=0.1678$， $C_s=0.1678$这个数值是从均匀各向同性湍流推导出来的。

https://www.openfoam.com/documentation/guides/latest/doc/guide-turbulence-les-smagorinsky.html

---

【讨论2】

• 对于不可压缩湍流，式(13)的化简是基于 $\mathrm{d}\mathrm{i}\mathrm{v}(\overline{\mathbf{U}})=0$得到的，但根据 pisoFoam计算的{U}，求div(U)结果为何不是严格等于0， $\mathrm{d}\mathrm{i}\mathrm{v}(\overline{\mathbf{U}})$实际上并不严格等于0，化简的式(15)与原公式(11)不是严格等价的。

【OpenFOAMv2012】Smagorinskyの解説 | mtk_birdman's blog

---

【讨论3】

《Sullivan, P.P., J.C. McWilliams and C. Moeng, A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows. Boundary-Layer Meteorology, 1994. 71: p. 247-276.》对第2种方法评价：
An alternative SGS model is to solve the TKE equation explicitly. Advantages of such models are that no equilibrium assumption is required, and the prognostic equation provides a direct means of calculating the SGS kinetic energy which is needed to construct the actual pressure.

• Sullivan(1994)对第2种方法评价，意思是根据式(4)，真实的压力应减去 $k_{sgs}$，即 $\overline{p}=\overline{P}-\frac{2}{3}\rho k_{\mathrm{s}\mathrm{g}\mathrm{s}}$。但OpenFOAM求解的是$\overline{P}=\overline{p}+\frac{2}{3}\rho k_{\mathrm{s}\mathrm{g}\mathrm{s}}$，所以无论 $k_{sgs}$实际值是多少，只要 ${\nu }_{sgs}$是相同的，那么得到压力结果也会是相同的。

• 这也正是为何我们只要控制$C_s={\left(C_k\sqrt{\frac{C_k}{C_e}}\right)}^{1/2}$为固定值0.065，即使改变 $C_k$值使得 $k_{sgs}$值很大，得到的压力结果也差别不大。但此时我们如果求解析$k_{resolved}$占比，就会对结论产生很大的影响。

• 当$C_s={\left(C_k\sqrt{\frac{C_k}{C_e}}\right)}^{1/2}$设置为固定值，根据式(16)，${\nu }_{sgs}$也为固定值。根据$k_{\mathrm{s}\mathrm{g}\mathrm{s}}={\left(\frac{{\nu }_{sgs}}{C_k\mathrm{\Delta }}\right)}^2$，此时只要 $C_k$取值很小，$k_{sgs}$的值就会很大，但计算的解析$k_{resolved}$基本不变，因此解析湍流能占比就会很小。但如果$C_k$取值很大，$k_{sgs}$的值就会很小，解析湍流能占比就会较大。所以检查LES的解析湍动能占比，Smagorinsky模型经验系数选择有很大影响。

• 为什么第2种方法结果是满足局部平衡假设，还会出现上述问题？原因是 ${\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}=C_k\mathrm{\Delta }\sqrt{k_{\mathrm{s}\mathrm{g}\mathrm{s}}}$这个关系式是一个假设关系，对于不同湍流，$C_k$应有一个合理取值范围，而不能任意取，需要通过实验进行校准。

---

4. 第1种方法讨论

第1种方法的植入如式(18)所示。

$$\begin{array}{}\text{(18)}& {\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}={\left(C_s\mathrm{\Delta }\right)}^2\sqrt{2\overline{\mathbf{D}}:\overline{\mathbf{D}}}={\left(C_s\mathrm{\Delta }\right)}^2\left|\overline{\mathbf{D}}\right|\end{array}$$

第1种方法只能设置$C_s$一个参数得到合理的${\nu }_{sgs}$，但无法直接得到$k_{sgs}$。但OpenFOAM最后求解得到的压力是$\overline{P}=\overline{p}+\frac{2}{3}\rho k_{\mathrm{s}\mathrm{g}\mathrm{s}}$，所以没有求解$k_{sgs}$也对结果无影响。

---

【讨论4】

$k_{sgs}$的定义如式(7)所示，因此要求解$k_{sgs}$得从求解 ${\tau }_{sgs}$ 入手。根据式(8)，目前未知数有2个未知数：$k_{sgs}$,${\tau }_{sgs}$。与第2种方法不同的地方，我们现在是已知${\nu }_{sgs}$，求解$k_{sgs}$。再根据 Local equilibrium，式(11)的前2行：

结果发现式(19)中仍然还有2个未知数：$k_{sgs}$, $C_{\epsilon }$
因此，第1种方法要求解$k_{sgs}$，还需要输入系数$C_{\epsilon }$。
或者输入系数$C_k$，根据 ${\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}=C_k\mathrm{\Delta }\sqrt{k_{\mathrm{s}\mathrm{g}\mathrm{s}}}$反算$k_{sgs}$。注意这里的$C_k$与$C_s$并不直接等价。

---

【讨论5】

第1种方法：已知${\nu }_{sgs}$

• 输入系数$C_{\epsilon }$，根据式(19)就能求解$k_{sgs}$。这个时候计算结果自然满足Local equilibrium假设。
• 或者输入系数$C_k$，根据 ${\nu }_{\mathrm{s}\mathrm{g}\mathrm{s}}=C_k\mathrm{\Delta }\sqrt{k_{\mathrm{s}\mathrm{g}\mathrm{s}}}$反算$k_{sgs}$。再将${\nu }_{sgs}$,$k_{sgs}$代入式(19)就能求解满足Local equilibrium假设的系数$C_{\epsilon }$
• 因此Sullivan(1994)对第1种方法的评价是：which results from an assumed local equilibrium balance between shear production and dissipation in the SGS turbulent kinetic energy (TKE) equation.

第2种方法：已知 $C_k,C_{\epsilon }$

• 根据式(11)求解$k_{sgs}$ ，结果自然满足Local equilibrium假设。

综上，Smagorinsky湍流模型的2种植入方式都可以满足Local equilibrium假设。

---

5. 总结

• 第1种方法只设置1个参数$C_s$只能得到${\nu }_{sgs}$，要求解$k_{sgs}$还需要额外指定系数$C_k$或者$C_{\epsilon }$。
• 第2种方法通过设置2个参数$C_k,C_{\epsilon }$，可以得到${\nu }_{sgs}$,$k_{sgs}$，但同样的等价$C_s={\left(C_k\sqrt{\frac{C_k}{C_{\epsilon }}}\right)}^{1/2}$前提下，$C_k$系数的不同取值对$k_{sgs}$影响很大。
• OpenFOAM求解的是$\overline{P}=\overline{p}+\frac{2}{3}\rho k_{\mathrm{s}\mathrm{g}\mathrm{s}}$，所以只要 ${\nu }_{sgs}$是相同的，那么得到压力结果也会是相同的。因此第2种方法通过设置 $C_k,C_{\epsilon }$得到等价的$C_s={\left(C_k\sqrt{\frac{C_k}{C_{\epsilon }}}\right)}^{1/2}$，第1个方法设置同样的$C_s$值，此时两种方法的结果是近似等价的。
• 两个方法设置等价的$C_s$值，结果仍有点小差异，主要是因为速度压力耦合求解过程中，$\mathrm{d}\mathrm{i}\mathrm{v}(\overline{\mathbf{U}})$实际并不严格等于0，导致两个方法计算的${\nu }_{sgs}$有略微差别，但整体结果是很接近的。
• Smagorinsky湍流模型要得到合理的结果，需要确定$C_s,C_k,C_{\epsilon }$三者中的任意2个即可。但因为OpenFOAM求解的是$\overline{P}=\overline{p}+\frac{2}{3}\rho k_{\mathrm{s}\mathrm{g}\mathrm{s}}$，可以不需要求解$k_{sgs}$的值，因此获得合理结果只要确定合理的$C_s$即可。
• 目前模拟channel flow，用的比较多取值是$C_s=0.065$。对应的，OpenFOAM自带算例planeChannel建议系数取值为$C_k=0.0265463553,C_{\epsilon }=1.048$

根据《Moin, P. and J. Kim, Numerical investigation of turbulent channel flow. Journal of fluid mechanics, 1982. 118: p. 341-377.》，建议取Cs=0.065

/OpenFOAM-v2206/tutorials/incompressible/pimpleFoam/LES/planeChannel

    SmagorinskyCoeffs
    {
        Ce              1.048;
        Ck              0.0265463553; // Updated to give Cs = 0.065
    }

[李东岳]

@coolhhh 哎我去，你这搞完了，可以写个Smagorinsky的review了...

[Voynich]

@coolhhh 太牛逼了！膜拜！另外关于压力没有减去k的问题，Fluent文档（2022R2版本, Ansys Fluent Theory Guide, p118）里有这样的论述：

文献是《G. Erlebacher, M. Y. Hussaini, C. G. Speziale, and T. A. Zang. "Toward the Large-Eddy Simulation
of Compressible Turbulent Flows". J. Fluid Mech. 238. 155–185. 1992.》，是否可以认为k就是忽略了

[coolhhh]

@Voynich 目前看到的处理都是把k忽略掉，这个帖子OpenFOAM 不可压缩湍流模型的 divDevReff 函数 | Giskard's CFD Learning Tricks中的2个观点，其实本质都是k被忽略掉，只是不知清楚k忽略对结果影响具体是多大，我也没具体搜过有无人做过这方面的比较

[coolhhh]

最近看《Holzmann T. Mathematics, numerics, derivations and OpenFOAM》，里边也提到大部分情况忽略k，对结果影响不大：