Smagorinsky模型系数问题

[李东岳]

@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}$差异较大的工况，对比压力是否可能能对比出结果差异？

[coolhhh]

@cfdem小白 应该是笔误，of的默认系数见这个链接

[coolhhh]

@coolhhh 更正下符号，当满足Local equilibrium时，$P=ϵ$

[coolhhh]

@李东岳 李老师，我把您的程序跑了对比下

1. case1

用第一种方法mySmagorinsky
设置Ck=0.065,也就是Cs=0.065，此时$k_{sgs1}={\nu }_{sgs1}^2/(C_s\mathrm{\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\mathrm{\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\mathrm{\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）脉动压力

[coolhhh]

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\mathrm{\Delta }{)}^2$中的系数$1/(C_k^2)=146407$，这个时候计算的$k_{sgs}$数值相较于前面的工况大非常多。
（2）case4模拟结果
case4的瞬时速度剖面和平均压力剖面，和cae1（第一种方法）非常相近。这是否可以推测：第一种方法与第二种方法效果要基本一样时，此时的Ce和Ck取值要与case4的取值范围类似？

[cfdem小白]

@coolhhh 大佬，这边问你一个问题。你图中的脉动速度Uprime2Mean应该是基于一段时间计算且在某一时刻显示的，这段时间的平均值吧？如果想要知道某时刻的瞬时脉动速度，是否是用该时刻的瞬时速度U减去该时刻的Umean得到？

[coolhhh]

Uprime2Mean是fieldAverage的计算结果，不是瞬时的脉动速度。风工程领域经常把速度或压力方差，称作脉动速度或脉动压力。如果要求瞬时的脉动速度，那应该按照你说的这么算

@cfdem小白 在 Smagorinsky模型系数问题 中说：

如果想要知道某时刻的瞬时脉动速度，是否是用该时刻的瞬时速度U减去该时刻的Umean得到？

[李东岳]

针对你的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$的关系。