中性大气环境湍流动能的自保持 | 附有算例下载

coolhhh

@李东岳或许是我这边fluent有某个设置错了，杨老师论文中顶面用的slip结果不会出现k衰减。我再试了fluent中顶面设置no-slip，结果是比较符合预期，接近顶面壁面平均速度减小，湍动能增大
Ux剖面.png TKE剖面.png

SHUKK

@李东岳李老师，我最近尝试复现Yang2017的kOmegaSST，按照文献的设置和参数设置，k都无法保持稳定。不知道老师您是否复现出来了。

李东岳

@SHUKK kOmegaSST没有尝试。你可以把你的算例打包发上来我给你看看

SHUKK

@李东岳好的，谢谢李老师！y2017kOmegaSST.zip

chszkc

@Vortex 前辈您好，请问可以向您请教一下LESGO温度分层算例的相关设置吗？

疏影横斜水清浅

请问有没有人做大涡模拟的中性大气环境湍动能自保持，大涡模拟的有没有参考文献？

疏影横斜水清浅

@李东岳老师您好，请问Yang2009中的垂向网格设置，第一层网格必须很小吗，和y+差不多。如果将垂向网格的尺寸变大，算例会发散，有人遇到过这种情况吗？

李东岳

@疏影横斜水清浅

如果将垂向网格的尺寸变大，算例会发散，有人遇到过这种情况吗？

我提供的这个算例么

疏影横斜水清浅

@李东岳是的

李东岳

@疏影横斜水清浅明天我试一下我记得我放大过很大的网格也能算

疏影横斜水清浅

@李东岳测试的网格如下，网格尺寸比较大

convertToMeters 1;

vertices
(
    (0 0 0)//0
    (5000 0 0)//1
    (5000 500 0)//2
    (0 500 0)//3
    (0 0 1)//4
    (5000 0 1)//5
    (5000 500 1)//6
    (0 500 1)//7
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (500 50 1) simpleGrading (1 10 1)
);

edges
(
);

boundary
(
    inlet
    {
        type       patch;
        faces
        (
             (3 0 4 7)
        );
    }
    outlet
    {
        type       patch;
        faces
        (
             (1 2 6 5)
        );
    }
    topWall
    {
        type       patch;
        faces
        (
             (2 3 7 6)
        );
    }
    flatWall1
    {
        type       wall;
        faces
        (
             (0 1 5 4)
        );
    }
    front
    {
        type empty;
        faces
        (
             (4 5 6 7)
        );
    }
    back
    {
        type empty;
        faces
        (
             (3 2 1 0)
        );
    }
);

mergePatchPairs
(
);

李东岳

你这个不是网格尺寸变大，是网格尺度变大。网格尺寸放大，是计算域不变，网格变粗。

yang本身是小尺度，你变成了大气环境大尺度，类似工业放大。或许会有问题。

疏影横斜水清浅

@李东岳这样呀，那我明白了，HW07那个可以放大计算域，但是那个的缺点是湍动能不随高度变化，所以我想着Yang09的也可以放大的，试了试，每次都发散

李东岳

@疏影横斜水清浅是，HW07本身计算雨很大，可以放大计算域，Yang2009这个我记得我尝试也不行，现在有点记不清了。

疏影横斜水清浅

@李东岳我明白了，HW07的湍动能和耗散率都是定值，不随高度变化，而Yang09的乘了一个系数，这个系数是个开平方数，而里面的数随着高度的增加会减小，当高度超过一定高度（与粗糙度系数有关）时，这个数就变成了负数，那么开方运算就不可以了，程序报错。

李东岳

@疏影横斜水清浅厉害厉害！非常厉害！从算法上解决了问题。有没有什么解决办法呢。

疏影横斜水清浅

@李东岳如果想要改进，就需要改公式里的C1和C2这两个参数，使其可以适应大尺度，Yang2009中的应该是根据实验数据或者其他的进行曲线拟合、调整得到的，我理解的，如果要调整就需要一些野外的实验观测数据什么的，根据实验得到一些经验参数。我暂时解决不了，这些经验参数不能拍脑袋确定，需要一定的依据。再有方法就是对后面乘的那个系数公式进行改进了，这都需要风工程的大佬们来试验解决了。我暂时也就想到这些，抛个砖

李东岳

@疏影横斜水清浅屌屌屌，厉害厉害

这位老铁是计算风工程专业的， @coolhhh 怎么看这个问题

coolhhh

@coolhhh 在中性大气环境湍流动能的自保持 | 附有算例下载中说：

Yang2009中，提到了SKE1和SKE2两种模型，并模拟对比。

尝试根据式(25)和(26)把这两个k曲线画出对比，这两条曲线还是有挺大区别，但在原论文中两个工况的inlet结果看起来像是一样的

更正之前的错误

再仔细看了论文，SKE1模型用的是 $C_\mu=0.09$，SKE2模型用的是 $C_\mu=0.028$，两条曲线结果是基本一致的。
k剖面对比.png

coolhhh

@李东岳李老师，我重新看了下原论文，做了些测试：

1. Yang2009论文中`k`剖面是缩尺的风洞拟合结果，不能用于全尺寸。

Yang2009论文中k剖面，随着高度在增加如下所示。
k剖面对比-10m.png

2. 缩尺风洞目标值换算为全尺寸

（1）缩尺风洞目标值换算为全尺寸，要考虑缩尺比，包括：几何缩尺比（$S_L=L_{model}/L_{full}$），速度缩尺比（$S_V=V_{model}/V_{full}$），时间缩尺比（$S_T=T_{model}/T_{full}=S_L/S_V$）。
（2）对数率平均风、湍动能、耗散率剖面：
$$u=\frac{u_*}{\kappa}\rm{ln}(\frac{z+z_0}{z_0}) \tag{1}$$

$$k=\frac{u_*^2}{\sqrt{C_\mu}} \sqrt{C_1\rm{ln}(\frac{z+z_0}{z_0})+C_2} \tag{2}$$

$$\epsilon=\frac{u_*^3}{\kappa (z+z_0)} \sqrt{C_1\rm{ln}(\frac{z+z_0}{z_0})+C_2} \tag{3}$$

从缩尺换算到全尺寸，需换算的参数：$u_* /S_V$，$z_0 /S_L$

（3）Yang2009原论文中的风洞试验拟合结果如下，其中$k/u_* ^2$剖面沿着高度的均值等于2.86（SKE1）和2.91（SKE2）。

（4）假如取$S_L=1/400$，$S_V=1/6$（风洞试验10m/s，全尺寸取60m/s），对Yang2009论文的结果放大为全尺寸结果如下，其中$k/u_* ^2$剖面沿着高度的均值等于3.17（SKE1）和3.22（SKE2）

3. 对上述全尺寸的参数，进行RANS模拟

（1）对原算例网格每个方向都放大400倍（of9写法与of8不一样）

# transform to full scale
runApplication transformPoints "scale=(400 400 400)"

（2）原算例设置算全尺寸会发散，修改fvSchemes，增加 grad(U) cellLimited Gauss linear 1;：

gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;
}

（3）0文件夹的U、k、epsilon对应修改为全尺寸参数
算例下载链接为：Yang2009FullScale_20230610.zip ，用of9执行Allrun或Allrun-parallel

（4）模拟结果如下，平均风剖面不变，但k基本衰减为0。

4.分析与疑问

（1）Yang2009原论文提到，参数 $k/u_* ^2$是根据风洞试验测量得到，具体是如何通过$k/u_* ^2$计算得到$C_\mu$？

（2）SKE1和SKE2模型，参数$C_\mu, C_1, C_2$组合都能很好的拟合实验结果，但RANS模拟自保持结果差别很大，其中$C_\mu$感觉像是重要影响参数。问题同（1），如何根据不同目标值，合理确定$C_\mu$？

（3）全尺寸模拟的k完全无法自保持，猜测原因：
①直接用缩尺目标值放大为全尺寸，参数组合$C_\mu, C_1, C_2$可能不一定适用。
②数值格式对全尺寸网格耗散更大，耗散程度对于网格尺寸增加成非线性增加的关系。但也尝试过把每个方向网格加密一倍，145万网格的算例结果也是跟上面的工况一致，k基本耗散为0

（4）OpenFOAM中提供的壁面函数，是否可用于全尺寸模拟？目前壁面函数的经验公式基本都是实验拟合出来的，虽然是无量纲化的。但能否用于全尺寸仍有待确定。本工况按原网格直接放大400倍，$y^+$值达到了2000~70000，这个数值对于全尺寸是否合理？

# y+ ()       
# Time        	patch         	min           	max           	average       
1000          	top	2.315308e+01	2.312809e+05	4.448700e+03
1000          	bottom	4.092494e+04	4.498426e+05	7.201201e+04
1000          	sides	4.473468e-01	2.178080e+05	2.860156e+03

（5）HW07模拟结果，k剖面的数值只有1.3左右，对于全尺寸而言这个数值是否太小不合理？如果HW07也模拟同样的目标值，结果会如何？

@李东岳在中性大气环境湍流动能的自保持 | 附有算例下载中说：

我在OpenFOAM里面植入了一个新的算法。我也把相关的内容更新在了《无痛苦NS方程笔记》。这个方法简称HW07，参考的是D.M. Hargreaves and N.G. Wright 2007，目前这一套方法计算出来的结果可以把k保持住。

（6）这个网格是二维平面网格（宽度只有一层），对于Yang2009的方法，能否直接用于二维网格还没有验证，并且目前对比的结果都是取计算域中间的结果，不考虑边界的影响。

@疏影横斜水清浅在中性大气环境湍流动能的自保持 | 附有算例下载中说：