Smagorinsky模型系数问题
-
@cfdem小白 可以看3#的帖子，这是《Moin, P. and J. Kim, Numerical investigation of turbulent channel flow. Journal of fluid mechanics, 1982. 118: p. 341-377.》建议取
Cs=0.065，然后根据of2206中planeChannel算例的设置，LES采用Smagorinsky，提到设置Ce=1.048，Ck=0.0265463553，能够Updated to give Cs = 0.065。我也用这个系数计算过Channel Flow，结果跟实验吻合还比较好 -
- 李老师，重新看了这篇文章《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.》，文中对这两种植入方法有个阐述：
（1）第一种方法：
$$\nu_{sgs}=(C_s\Delta)^2\sqrt{2\bfS:\bfS} \tag{1}$$
from an assumed local equilibrium balance between shear production and dissipation in the SGS turbulent kinetic energy (TKE) equation。（2）第二种方法：
$$\nu_{sgs}=C_k\Delta\sqrt{k_{sgs}} \tag{2}$$
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.
①看文章意思像是直接对文献中的方程(8)求解，当假设$P=D$的时候，就能与第一种方法得到等价的$C_s$:
$$C_s=\left ( C_k \sqrt{\frac{C_k}{C_e} } \right ) ^{1/2} \tag{3}$$
但根据Smagorinsky SGS model in OpenFOAM | CFD WITH A MISSION 的 推导，of中求解的方程是通过Local equilibrium得到的，即假设了$P=D$，才得到了$k_{sgs}$的求解方程，与文献中的直接求解TKE方程不需要平衡假设的说法不太一样。
②并且根据论文阐述，Moeng(1984)文章中提到$k_{sgs}$的值是用于计算压力
- 在23#提到第一种方法的$k_{sgs}$是反推出来的，根据式(2)可反推得到$k_{sgs}=\nu_{sgs}^2/(C_k\Delta)^2$。疑问是这里的$C_k$，对于第一种方法是否取值就是$C_s$？
如果$C_k$与$C_s$是不同的取值，那么对对第一种方法仍然还是有未知数$C_k$，无法求$k_{sgs}$；如果$C_k=C_s$，产生的疑问见第3个问题。
- 第一种方法如果$C_k=C_s$，根据式(3)，可知：
（1）第一种方法：取$C_s=0.065$，$\nu_{sgs1}$由式（1）计算，那么$k_{sgs1}=\nu_{sgs1}^2/(C_s\Delta)^2$
（2）第二种方法：$C_e=1.048，C_k=0.0265463553$，可得到等价的$C_s=0.065$，此时得到的$\nu_{sgs2}=\nu_{sgs1}$。那么根据式（2）得到$k_{sgs2}=\nu_{sgs2}^2/(C_k\Delta)^2$
（3）若要使得两种方法计算的$k_{sgs1}=k_{sgs2}$，由于$\nu_{sgs2}=\nu_{sgs1}$，因此要求
$$C_s=C_k \tag{4}$$
再联立式（3），可解得两个方法完全等价时需要满足条件：$C_s=C_k, C_e=1/C_k$。
- Moeng(1984)文章中提到$k_{sgs}$的值是用于计算压力。比较两个方法的不同，构造不同系数组合$C_k, C_e$, 使得两种方法有相同的$C_s=0.065$，但计算的$k_{sgs}$差异较大的工况，对比压力是否可能能对比出结果差异？
- 李老师，重新看了这篇文章《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.》，文中对这两种植入方法有个阐述：
-
@李东岳 李老师，我把您的程序跑了对比下
1. case1
用第一种方法
mySmagorinsky
设置Ck=0.065,也就是Cs=0.065，此时$k_{sgs1}=\nu_{sgs1}^2/(C_s\Delta)^2$中的系数$1/(C_s^2)=236.68$2. case2
用第二种方法：of10中的
Smagorinsky
设置Ce=1.048，Ck=0.0265463553, 满足Cs=0.065。
$k_{sgs2}=\nu_{sgs2}^2/(C_k\Delta)^2$中的系数$1/(C_k^2)=1419$3. case3
用第二种方法：of10中的
Smagorinsky
根据30#的猜测，设置了一组可能不合理的组合：Ce=56020, Ck=1, 但能满足Cs=0.065。
此时$k_{sgs2}=\nu_{sgs2}^2/(C_k\Delta)^2$中的系数$1/(C_k^2)=1$，与前面两个算例相差特别大
根据$C_e=C_k^3/C_s^4$，可画出下面曲线。
4. 结果对比
根据模拟结果，似乎结果都差别不大。对于of10中的
Smagorinsky，感觉只要构造出等价的Cs=0.065，结果基本不变，也就是说$k_{sgs}$似乎没有参与压力计算。（1）瞬时速度
（2）平均速度
（3）脉动速度
（4）瞬时压力
（5）平均压力
（6）脉动压力
-
1. 补充说明
因为$k_{sgs}$的值与系数$C_k/C_e$有关。上面33#中的case2和case3，虽然$k_{sgs}$差个1419倍，但注意到这时$C_k/C_e$的数值已经很小了，所以计算结果差别不大。
2. 增加case4
（1）of10中的
Smagorinsky，设置Ce=0.001，Ck=0.002613472，同样能够满足Cs=0.065。
$k_{sgs2}=\nu_{sgs2}^2/(C_k\Delta)^2$中的系数$1/(C_k^2)=146407$，这个时候计算的$k_{sgs}$数值相较于前面的工况大非常多。
（2）case4模拟结果
case4的瞬时速度剖面和平均压力剖面，和cae1（第一种方法）非常相近。这是否可以推测：第一种方法与第二种方法效果要基本一样时，此时的Ce和Ck取值要与case4的取值范围类似？
-
针对你的case2与case3，你看看nut以及k的结果，看看有区别么？
把下面的放在controlDict里面可以输出$k_{sgs}$
cacheTemporaryObjects ( k ); functions { #includeFunc writeObjects("k") }我觉得这些东西，可以写个英文让大家知道一下了。投个一般的期刊当做一个note都可以。目前我是没发现有人研究这个。并且我认为研究很有意义。
主要是看这东西之前是不是确实没人玩过。如果没人搞过，那研究一下绝对是发现。写成英文的目的不是为了发文章，而是真正的让大家知道一下这个事情：1）两种Smagorinsky的植入区别与对比，2）第二种植入的模型系数问题。我是没发现有人研究这个。 -
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的结果，与前面工况的速度场有点不同，还有计算的数值误差等其他原因导致。
-
- 李老师，我也没具体看到哪个文章用了这两个默认系数，但搜到一个帖子有个观点，帖子中的超链接无法查看：OpenFOAM大涡模拟湍流模型之Smagorinsky模型代码详解
- 无意中看到一个日本CFD网站，看起来$k_{sgs}$与压力计算值有关。
然后看Smagorinsky的推导，这个推导过程也类似的出现$Ck$和$C_e$参数。疑问依然还是第一种方法是如何计算$k_{sgs}$，如果这个问题知道了，那应该能推断出第一种方法与第二种方法的系数$Ck$和$C_e$的关系。
-
@李东岳 李老师，最近把
Smagorinsky湍流模型两种方式植入方式再梳理了下，下面只讨论不可压缩湍流情况，重新完整总结如下（由于有的公式latex显示异常，就采用截图形式）：1. OpenFOAM中LES求解方程形式
根据 【OpenFOAMv2012】Smagorinskyの解説 | mtk_birdman's blog
，对于不可压缩湍流，OpenFOAM中求解的方程如下
其中，
\begin{equation}
\nu_{eff}=\nu+\nu_{sgs} \tag{3}
\end{equation}\begin{equation}
\overline{P}=\overline{p}+\frac{2}{3}\rho k_{sgs} \tag{4}
\end{equation}在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{equation}
\mathbf{\bar{D}}=\frac12\Big(\nabla\mathbf{\bar{U}}+\nabla\mathbf{\bar{U}}^\mathrm{T}\Big) \tag{5}
\end{equation}$$亚格子应力张量（subgrid scale stress tensor）：
其中 $k_{sgs}$定义为
\begin{equation}
k_{sgs}=\frac12\operatorname{tr}\left(\boldsymbol{\tau}_{sgs}\right) \tag{7}
\end{equation}
根据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{equation}
P=\varepsilon\Leftrightarrow\overline{\mathbf{D}}:\mathbf{\tau_{sgs}}+C_\varepsilon\frac{k_{\mathrm{sgs}}^{1.5}}\Delta=0 \tag{9}
\end{equation}
其中局部平衡假设公式参考资料：
- 《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{equation}
\nu_{\mathrm{sgs}}=C_k\Delta\sqrt{k_{\mathrm{sgs}}} \tag{10}
\end{equation}$$事实上，式(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{equation}
\left|\bar{\mathbf{D}}\right|=\sqrt{2\bar{\mathbf{D}}:\bar{\mathbf{D}}} \tag{14}
\end{equation}把式(13)代入式(11)，得到
\begin{equation}
k_{\mathrm{sgs}}=\frac ca=\frac{C_k\Delta^2\left|\overline{\mathbf{D}}\right|^2}{C_\varepsilon} \tag{15}
\end{equation}\begin{equation}
\nu_{\mathrm{sgs}}=C_k\Delta\sqrt{k_{\mathrm{sgs}}}=C_k\sqrt{\frac{C_k}{C_\varepsilon}}\Delta^2\left|\overline{\mathbf{D}}\right| \tag{16}
\end{equation}式(16)对比第1种方法式(18)，可得到关系式：
\begin{equation}
C_s=\left ( C_k \sqrt{\frac{C_k}{C_\varepsilon} } \right ) ^{1/2} \tag{17}
\end{equation}OpenFOAM中，默认系数设置为 $C_\varepsilon=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{div}(\bar{\mathbf{U} })=0$得到的，但根据 pisoFoam计算的{U}，求div(U)结果为何不是严格等于0， $\mathrm{div}(\bar{\mathbf{U}})$实际上并不严格等于0，化简的式(15)与原公式(11)不是严格等价的。
【讨论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 energywhich is needed to construct the actual pressure.-
Sullivan(1994)对第2种方法评价，意思是根据式(4)，真实的压力应减去 $k_{sgs}$，即 $\overline{p}=\overline{P}-\frac23\rho k_{\mathrm{sgs}}$。但OpenFOAM求解的是$\overline{P}=\overline{p}+\frac23\rho k_{\mathrm{sgs}}$，所以无论 $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{sgs}}=\left(\frac{\nu_{sgs}}{C_k\Delta}\right)^2$，此时只要 $C_k$取值很小，$k_{sgs}$的值就会很大，但计算的解析$k_{resolved}$基本不变，因此解析湍流能占比就会很小。但如果$C_k$取值很大，$k_{sgs}$的值就会很小，解析湍流能占比就会较大。所以检查LES的解析湍动能占比，Smagorinsky模型经验系数选择有很大影响。
-
为什么第2种方法结果是满足局部平衡假设，还会出现上述问题？原因是 $\nu_{\mathrm{sgs}}=C_k\Delta\sqrt{k_{\mathrm{sgs}}}$这个关系式是一个假设关系，对于不同湍流，$C_k$应有一个合理取值范围，而不能任意取，需要通过实验进行校准。
4. 第1种方法讨论
第1种方法的植入如式(18)所示。
$$\begin{equation}
\nu_{\mathrm{sgs}}=\left(C_s\Delta\right)^2\sqrt{2\bar{\mathbf{D}}:\bar{\mathbf{D}}}=\left(C_s\Delta\right)^2\left|\bar{\mathbf{D}}\right| \tag{18}
\end{equation}$$第1种方法只能设置$C_s$一个参数得到合理的$\nu_{sgs}$，但无法直接得到$k_{sgs}$。但OpenFOAM最后求解得到的压力是$\overline{P}=\overline{p}+\frac23\rho k_{\mathrm{sgs}}$，所以没有求解$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_\varepsilon$
因此，第1种方法要求解$k_{sgs}$，还需要输入系数$C_\varepsilon$。
或者输入系数$C_k$，根据 $\nu_{\mathrm{sgs}}=C_k\Delta\sqrt{k_{\mathrm{sgs}}}$反算$k_{sgs}$。注意这里的$C_k$与$C_s$并不直接等价。
【讨论5】
第1种方法：已知$\nu_{sgs}$
- 输入系数$C_\varepsilon$，根据式(19)就能求解$k_{sgs}$。这个时候计算结果自然满足Local equilibrium假设。
- 或者输入系数$C_k$，根据 $\nu_{\mathrm{sgs}}=C_k\Delta\sqrt{k_{\mathrm{sgs}}}$反算$k_{sgs}$。再将$\nu_{sgs}$,$k_{sgs}$代入式(19)就能求解满足Local equilibrium假设的系数$C_\varepsilon$
- 因此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_\varepsilon$
- 根据式(11)求解$k_{sgs}$ ，结果自然满足Local equilibrium假设。
综上，Smagorinsky湍流模型的2种植入方式都可以满足Local equilibrium假设。
5. 总结
- 第1种方法只设置1个参数$C_s$只能得到$\nu_{sgs}$，要求解$k_{sgs}$还需要额外指定系数$C_k$或者$C_\varepsilon$。
- 第2种方法通过设置2个参数$C_k,C_\varepsilon$，可以得到$\nu_{sgs}$,$k_{sgs}$，但同样的等价$C_s=\left ( C_k \sqrt{\frac{C_k}{C_\varepsilon} } \right ) ^{1/2} $前提下，$C_k$系数的不同取值对$k_{sgs}$影响很大。
- OpenFOAM求解的是$\overline{P}=\overline{p}+\frac23\rho k_{\mathrm{sgs}}$，所以只要 $\nu_{sgs}$是相同的，那么得到压力结果也会是相同的。因此第2种方法通过设置 $C_k,C_\varepsilon$得到等价的$C_s=\left ( C_k \sqrt{\frac{C_k}{C_\varepsilon} } \right ) ^{1/2} $，第1个方法设置同样的$C_s$值，此时两种方法的结果是近似等价的。
- 两个方法设置等价的$C_s$值，结果仍有点小差异，主要是因为速度压力耦合求解过程中，$\mathrm{div}(\bar{\mathbf{U}})$实际并不严格等于0，导致两个方法计算的$\nu_{sgs}$有略微差别，但整体结果是很接近的。
- Smagorinsky湍流模型要得到合理的结果，需要确定$C_s,C_k,C_\varepsilon$三者中的任意2个即可。但因为OpenFOAM求解的是$\overline{P}=\overline{p}+\frac23\rho k_{\mathrm{sgs}}$，可以不需要求解$k_{sgs}$的值，因此获得合理结果只要确定合理的$C_s$即可。
- 目前模拟
channel flow，用的比较多取值是$C_s=0.065$。对应的，OpenFOAM自带算例planeChannel建议系数取值为$C_k=0.0265463553,C_\varepsilon=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/planeChannelSmagorinskyCoeffs { Ce 1.048; Ck 0.0265463553; // Updated to give Cs = 0.065 }
