Smagorinsky模型系数问题

[李东岳]

@coolhhh thats a good questions.

不过目前我就没见过详细讨论OpenFOAM植入smagorinsky的方法，更别提可压缩适配了。大佬可以搜搜sci，你这整的挺细啊，这些讨论可以做sci的appendix了

[李东岳]

https://github.com/Unofficial-Extend-Project-Mirror/openfoam-extend-foam-extend-3.1/blob/master/src/turbulenceModels/incompressible/LES/dynSmagorinsky/dynSmagorinsky.C

我最近看到一个老版本的动态模型植入。植入的不是${\nu }_{sgs}={\left(C_s\mathrm{\Delta }\right)}^2\sqrt{2\mathbf{D}:\mathbf{D}}$，而是${\nu }_{sgs}={\left(C_s\mathrm{\Delta }\right)}^2\sqrt{\mathbf{D}:\mathbf{D}}$，大佬怎么看

[coolhhh]

@李东岳 李老师， 这两个链接在讨论of中的Smagorinsky模型时，也类似讨论of中为何植入的是$\sqrt{\mathbf{D}:\mathbf{D}}$，相比理论公式看似缺少$\sqrt{2}$系数。但其实是一致的，这是因为of中把$\sqrt{2}$加到了前面的系数中。下面截图中的理论公式符号代表：$\left|\overline{S}\right|=\sqrt{2\mathbf{D}:\mathbf{D}}$，但OF公式代表：$\left|{\overline{S}}_{of}\right|=\sqrt{\mathbf{D}:\mathbf{D}}$，在CFD-online上的帖子因为把符号写成一样导致了混淆。
Smagorinsky model details -- CFD Online Discussion Forums
[转载]有关OpenFoam中LES smagorinsky模型的推导_馨予_新浪博客

---

湍流模型我不太懂，在想dynSmagorinsky是否类似于CFD-online的讨论，是否也是把系数$\sqrt{2}$加到了前面的系数中？也就是cD(D)可能比dynSmagorinsky理论公式多了$\sqrt{2}$系数？

void dynSmagorinsky::updateSubGridScaleFields(const volSymmTensorField& D)
{
    nuSgs_ = cD(D)*sqr(delta())*sqrt(magSqr(D));
    nuSgs_.correctBoundaryConditions();
}

[李东岳]

是不是这个公式里面，也应该有一个$\sqrt{2}$？ 然后就可以了？

[coolhhh]

@李东岳 这个公式结论是对的，不用加$\sqrt{2}$。CFD-online讨论帖子初始问题混淆点是把OF计算的$\sqrt{\mathbf{D}:\mathbf{D}}$当做$\left|\overline{S}\right|$，而理论公式的$\left|\overline{S}\right|$表示的是$\sqrt{2\mathbf{D}:\mathbf{D}}$。OF是单独计算了$\sqrt{\mathbf{D}:\mathbf{D}}$，为与理论对应，$\sqrt{2}$在前面的系数里计算再相乘。

Smagorinsky SGS model in OpenFOAM | CFD WITH A MISSION这个帖子推导就没有混淆符号，统一用$\sqrt{2\mathbf{D}:\mathbf{D}}$表示$\left|\overline{S}\right|$

[李东岳]

你的意思就是，OpenFOAM里面植入的应该是${\nu }_{sgs}=(C_s^{of}{)}^2{\mathrm{\Delta }}^2\sqrt{\mathbf{D}:\mathbf{D}}$，其中$(C_s^{of}{)}^2=C_s^2\sqrt{2}$

而不是${\nu }_{sgs}=C_s^2{\mathrm{\Delta }}^2\sqrt{\mathbf{D}:\mathbf{D}}$

是这样么

[coolhhh]

对，因为湍流模型不太懂，我猜测dynSmagorinsky差个$\sqrt{2}$可能是这个原因。Smagorinsky模型不是直接植入这个公式，是通过把系数代入${\nu }_{sgs}=C_k\mathrm{\Delta }\sqrt{k_{sgs}}$后得到等价的${\nu }_{sgs}=C_k\sqrt{\frac{2C_k}{C_e}}{\mathrm{\Delta }}^2\sqrt{\mathbf{D}:\mathbf{D}}=C_k\sqrt{\frac{C_k}{C_e}}{\mathrm{\Delta }}^2\sqrt{2\mathbf{D}:\mathbf{D}}={\left(C_s\mathrm{\Delta }\right)}^2\sqrt{2\mathbf{D}:\mathbf{D}}$

[李东岳]

那我知道了，我感觉那个应该是个bug。毕竟这个代码已经过时了。我跟作者确认一下。

[李东岳]

LES那面为啥所有文章非要定义$|\mathbf{S}|=\sqrt{2\mathbf{S}:\mathbf{S}}$呢，正常张量这玩意就应该是$|\mathbf{S}|=\sqrt{\mathbf{S}:\mathbf{S}}$。不知道最开始那个大佬咋想的。

被这玩意彻底搞蒙了

[dxl]

是不是这个

[李东岳]

@dxl 没太看明白 老铁

[李东岳]

@李东岳 在 Smagorinsky模型系数问题 中说：

你的意思就是，OpenFOAM里面植入的应该是${\nu }_{sgs}=(C_s^{of}{)}^2{\mathrm{\Delta }}^2\sqrt{\mathbf{D}:\mathbf{D}}$，其中$(C_s^{of}{)}^2=C_s^2\sqrt{2}$

而不是${\nu }_{sgs}=C_s^2{\mathrm{\Delta }}^2\sqrt{\mathbf{D}:\mathbf{D}}$

是这样么

https://bugs.openfoam.org/view.php?id=816

这个确实是个bug，不过OpenFOAM新版已经把动态Smagorinsky删掉了，所以无关紧要了。

[李东岳]

我按照${\nu }_t=(C_s\mathrm{\Delta }{)}^2\sqrt{2\mathbf{S}:\mathbf{S}}$的方式植入进入了。从云图来看没啥区别。希望有大佬来详细验证一下是否均衡。下面这个文件，可以在OpenFOAM-10下进行编译，步骤：

wmake
cd pitzDaily
./Allrun

Smagorinsky.tar.xz

注意算例中的$C_k$的值忘记改了，应该改为0.065

---

个人疑问是植入：${\nu }_{sgs}={\left(C_s\mathrm{\Delta }\right)}^2\sqrt{2\mathbf{D}:\mathbf{D}}$ 那么式 (1d) 中的k是怎么计算的？这会不会就是两个方法的主要差别？

@coolhhh 应该是从${\nu }_t$反推出来的，$k={\nu }_t^2/(C_k\mathrm{\Delta }{)}^2$

[coolhhh]

@李东岳 李老师，这个反推公式$k={\nu }_t^2/(C_k\mathrm{\Delta }{)}^2$是如何得到的，有更具体的推导吗？

[李东岳]

我刚才更新了下笔记，你看下4.2.5 RANS-LES 混合模型、DES 模型

还有这个：10.1007/s10494-005-6916-y

[cfdem小白]

@李东岳 李老师你好，这边Smagorinsky模型中Ce=1.048，Ck=0.065系数是比较好的吗

[李东岳]

这二者都可以得到cs=0.065，Ce=1.048，Ck=0.065是标准系数。然后$k_{sgs}$也跟这两个系数的比有关。不过取别的值也能达到这个效果。 @coolhhh 在这方面应该有更深的见解。

[coolhhh]

@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，结果跟实验吻合还比较好

[cfdem小白]

@coolhhh 这边李老师提到了Ce=1.048，Ck=0.065是标准系数，代入上述的Cs公式算出来Cs在0.127左右。

[coolhhh]

@李东岳

1. 李老师，重新看了这篇文章《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）第一种方法：
$$\begin{array}{}\text{(1)}& {\nu }_{sgs}=(C_s\mathrm{\Delta }{)}^2\sqrt{2\mathbf{S}:\mathbf{S}}\end{array}$$
from an assumed local equilibrium balance between shear production and dissipation in the SGS turbulent kinetic energy (TKE) equation。

（2）第二种方法：
$$\begin{array}{}\text{(2)}& {\nu }_{sgs}=C_k\mathrm{\Delta }\sqrt{k_{sgs}}\end{array}$$
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$:
$$\begin{array}{}\text{(3)}& C_s={\left(C_k\sqrt{\frac{C_k}{C_e}}\right)}^{1/2}\end{array}$$
但根据Smagorinsky SGS model in OpenFOAM | CFD WITH A MISSION 的 推导，of中求解的方程是通过Local equilibrium得到的，即假设了$P=D$，才得到了$k_{sgs}$的求解方程，与文献中的直接求解TKE方程不需要平衡假设的说法不太一样。
②并且根据论文阐述，Moeng(1984)文章中提到$k_{sgs}$的值是用于计算压力

---

2. 在23#提到第一种方法的$k_{sgs}$是反推出来的，根据式(2)可反推得到$k_{sgs}={\nu }_{sgs}^2/(C_k\mathrm{\Delta }{)}^2$。疑问是这里的$C_k$，对于第一种方法是否取值就是$C_s$？
   如果$C_k$与$C_s$是不同的取值，那么对对第一种方法仍然还是有未知数$C_k$，无法求$k_{sgs}$；如果$C_k=C_s$，产生的疑问见第3个问题。

---

3. 第一种方法如果$C_k=C_s$，根据式(3)，可知：
   （1）第一种方法：取$C_s=0.065$，${\nu }_{sgs1}$由式（1）计算，那么$k_{sgs1}={\nu }_{sgs1}^2/(C_s\mathrm{\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\mathrm{\Delta }{)}^2$
   （3）若要使得两种方法计算的$k_{sgs1}=k_{sgs2}$，由于${\nu }_{sgs2}={\nu }_{sgs1}$，因此要求
   $$\begin{array}{}\text{(4)}& C_s=C_k\end{array}$$
   再联立式（3），可解得两个方法完全等价时需要满足条件：$C_s=C_k,C_e=1/C_k$。

---

4. Moeng(1984)文章中提到$k_{sgs}$的值是用于计算压力。比较两个方法的不同，构造不同系数组合$C_k,C_e$, 使得两种方法有相同的$C_s=0.065$，但计算的$k_{sgs}$差异较大的工况，对比压力是否可能能对比出结果差异？