Smagorinsky模型系数问题
-
-
李 李东岳 被引用 于这个主题
-
对，因为湍流模型不太懂，我猜测
dynSmagorinsky差个$\sqrt 2$可能是这个原因。Smagorinsky模型不是直接植入这个公式，是通过把系数代入$\nu_{sgs}=C_k\Delta \sqrt{k_{sgs}}$后得到等价的$\nu_{sgs}=C_k \sqrt{\frac{2C_k}{C_e} } \Delta ^2 \sqrt{\mathbf{D} :\mathbf{D} }
=C_k \sqrt{\frac{C_k}{C_e} } \Delta ^2 \sqrt{2\mathbf{D} :\mathbf{D} }=\left ( C_s \Delta\right )^2 \sqrt{2\mathbf{D} :\mathbf{D} } $ -
@李东岳 在 Smagorinsky模型系数问题 中说：
你的意思就是，OpenFOAM里面植入的应该是$\nu_{sgs} = (C_s^{of})^2 \Delta^2 \sqrt{ \mathbf{D} :\mathbf{D} } $，其中$(C_s^{of})^2=C_s^2\sqrt{2}$
而不是$\nu_{sgs} = C_s^2 \Delta^2 \sqrt{ \mathbf{D} :\mathbf{D} } $
是这样么
https://bugs.openfoam.org/view.php?id=816
这个确实是个bug，不过OpenFOAM新版已经把动态Smagorinsky删掉了，所以无关紧要了。
-
我按照$\nu_t=(C_s\Delta)^2\sqrt{2\bfS:\bfS}$的方式植入进入了。从云图来看没啥区别。希望有大佬来详细验证一下是否均衡。下面这个文件，可以在OpenFOAM-10下进行编译，步骤：
wmake cd pitzDaily ./Allrun注意算例中的$C_k$的值忘记改了，应该改为0.065
个人疑问是植入：$ \nu _{sgs} =\left ( C_s \Delta\right )^2 \sqrt{2\mathbf{D} :\mathbf{D} }$ 那么式 (1d) 中的k是怎么计算的？这会不会就是两个方法的主要差别？
@coolhhh 应该是从$\nu_t$反推出来的，$k=\nu_t^2/(C_k\Delta)^2$
-
@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的取值范围类似？
