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

李东岳

我最近参考了杨老师的这一篇行业标杆文章New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer in computational wind engineering，杨老师在这篇文章里面写的很明确，并且看起来结果很好。我在openfoam里面弄进去了，但是结果对不上，仔细核对了几遍，也没看见植入的跟文章有什么差异。在此问问大佬，求个经验。

杨老师这篇文章里面，我总结了一下，

在进口，植入3个方程，分别是文章中的方程11、21、22，把这三个作为计算域的进口。其中$C_1=-0.17,C_2=0.62$
同时要把kEpsilon模型系数改一下：$C_1=1.5,C_2=1.92,C_\mu=0.028,\sigma_k=1.67,\sigma_\epsilon=2.51,u^*=0.511,z_0=0.000225$，参考3.1节最后
参考表一，nut使用粗糙壁面函数，Cs=0.75，Ks=0.0025

其他没啥了。我按照上面的设置，没出来理想的结果，下图中出口的k（右图中的红线，跟黑线有偏移）。

捕获.PNG

我个人怀疑是nut的壁面函数的因素。但是杨老师文章里面说用粗糙壁面函数就可以。想问问相关大佬，是不是我还忽略了什么？

更新：算例下载在这里 https://www.cfd-china.com/topic/6231/中性大气环境湍流动能的自保持/8

Hope

@李东岳李老师您好，我们这边之前在fluent里面有验证过该入流条件，具体结果可参见刘安彬的硕士论文《平衡大气边界层的CFD数值模拟》。但是在模拟时上部流场的湍动能并无明显衰减，只是近壁面有衰减，并在施加壁面剪应力之后可有效改善。您可以考虑下顶面边界条件、壁面剪应力这两方面的设置。
下图模拟结果截取自上述硕士论文，未加任何修正

下图是施加壁面剪应力后的模拟结果

李东岳

@Hope 我看你这面的常数跟杨老师那面不一样。我在主贴中更新了杨老师的这几个参数。是说这几个参数要调整么？你们这面的这几个参数是哪里得来的？

您可以考虑下顶面边界条件、壁面剪应力这两方面的设置。

yang2009这篇文章引用了近400，按道理来说应该可以通过这一篇sci复现相关的结果。如果考虑这些影响就可以正确的话，sci里面应该有提及。我在想是不是某些因素我没关注到，忽略掉了。比如这两个，在sci里面有隐含的提及么

Hope

@李东岳李老师，yang2009这篇sci论文里面我没看到有壁面剪应力的相关设置。由于这部分内容是我师兄做的，我先依据硕士论文解释您提出的几个问题，具体的我还需要再查证之后再给您回复。
硕士论文中引用的入流条件也是出自这篇sci文献，有出入的是壁面函数的粗糙高度取值，硕士论文中的粗糙高度是依据文献（Blocken et al. 2007, CFD simulation of the atmospheric boundary layer: wall function problems）计算所得：

壁面剪应力的取值是依据硕士论文中的一些假定提出的，而非yang2009中相关参数。

李东岳

@Hope 多谢多谢！我也研究研究。

李东岳

我还是在考虑这个问题。文章我没有时间太细的去查看细节了，我让一个学生去查看一下是否哪方面我还是设置的不到位。

李东岳

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

捕获.PNG

疏影横斜水清浅

@李东岳李老师你好，请问你这个保持住是指多大的距离，你在图中标的inlet、middle、outlet是多大的计算域？

李东岳

@疏影横斜水清浅进口处，2500米处，出口处。一共5000米

疏影横斜水清浅

@李东岳好的，了解了，谢谢李老师

coolhhh

@李东岳在中性大气环境湍流动能的自保持中说：

@疏影横斜水清浅进口处，2500米处，出口处。一共5000米

李老师，看杨老师论文中计算域尺寸比较小，与同济风洞实验室尺寸一致。有无可能是算的全尺寸，网格太粗糙导致结果不同

Vortex

能不能自保持，是要看壁面函数和入口能否协调的。我是用LES做大气边界层的。LES里跑预前模拟生成的流场就能让湍动能处处保持平衡。LES这边还有用入口生成方法的，但一般自保持不了，流体跑着跑着湍动能就会衰减，是因为壁面用的wall model和入口人工合成的流场不协调。

李东岳

@coolhhh 我Yang2009算的是Yang2009的算例，HW07算的是HW07的算例。我暂时还没有尝试用hw07的方法算yang2009，以及用yang2009的方法算hw07，这个过几天试一下...

李东岳

@Vortex 感谢分享，目前还没到LES那一阶段 :-) 下一阶段要上kOmega以及SST。LES估计要几个月后了。

Vortex

@李东岳介绍个专门做ABL的code，叫LESGO (https://lesgo.me.jhu.edu/)，Meneveau大佬以前是拿它研究LES模型的。不稳定、中性、稳定的ABL都可以做，还能考虑科氏力。

chszkc

@Vortex 在中性大气环境湍流动能的自保持中说：

@李东岳介绍个专门做ABL的code，叫LESGO (https://lesgo.me.jhu.edu/)，Meneveau大佬以前是拿它研究LES模型的。不稳定、中性、稳定的ABL都可以做，还能考虑科氏力。

想问一下LESGO有现成的考虑浮力驱动的求解器吗？

Vortex

@chszkc 有，解温度输运方程，求解器有植入。

李东岳

Yang2009.tar.xz

算例在这里

openfoam-9可以直接./Allrun运行。我是看不出什么bug了。但是结果对不上。希望各位大佬给看看

@coolhhh 直接给你发这了

coolhhh

@李东岳收到，谢谢老师

coolhhh

@李东岳李老师，算了您给的算例，计算结果跟您的结果一样，还没找到原因。但有以下问题：

blockMeshDict中把inlet和outlet设置成了wall，但都改为patch后，结果还是跟上面的计算结果基本一样，k偏差仍较大。

算例中blockMeshDict初始设置

convertToMeters 1;

vertices
(
    (0  0   0)
    (12 0   0)
    (12 1.8 0)
    (0  1.8 0)
    (0  0   1.8)
    (12 0   1.8)
    (12 1.8 1.8)
    (0  1.8 1.8)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (120 42 18) simpleGrading (1 4 1)
);

edges
(
);

boundary
(
    top
    {
        type patch;
        faces
        (
            (3 7 6 2)
        );
    }
    inlet
    {
        type wall;
        faces
        (
            (0 4 7 3)
        );
    }
    outlet
    {
        type wall;
        faces
        (
            (2 6 5 1)
        );
    }
    bottom
    {
        type wall;
        faces
        (
            (1 5 4 0)
        );
    }
    sides
    {
        type wall;
        faces
        (
            (0 3 2 1)
            (4 5 6 7)
        );
    }
);

Yang2009中网格，第一层网格高度0.01m，增长率是1.06。simpleGrading (1 4 1)这个设置反算对应第一层网格高度为0.0197，网格增长率为1.034，计算Y+结果大约在300：

# y+ ()       
# Time        	patch         	min           	max           	average       
0             	bottom	3.602087e+02	3.602087e+02	3.602087e+02
0             	sides	0.000000e+00	0.000000e+00	0.000000e+00
100           	bottom	3.026569e+02	4.241725e+02	3.196787e+02
100           	sides	2.660222e-03	1.234693e+02	2.231677e+01
200           	bottom	2.978864e+02	4.241725e+02	3.151323e+02
200           	sides	2.712471e-03	1.696112e+01	7.879751e-01
300           	bottom	2.978864e+02	4.241725e+02	3.151323e+02
300           	sides	3.569842e-03	1.696111e+01	5.413341e-01
400           	bottom	2.978864e+02	4.241725e+02	3.151323e+02
400           	sides	1.932649e-03	2.113821e+01	6.643026e-01
500           	bottom	2.978864e+02	4.241725e+02	3.151323e+02
500           	sides	1.392717e-03	1.696111e+01	5.073163e-01

尝试设置为simpleGrading (1 10.9 1)，网格总数不变，对应的网格和Yang2009中网格接近，第一层网格高度0.01m，增长率是1.06。计算的Y+数值在80左右，但计算结果还是跟上面的计算结果基本一样，k偏差仍较大。


# Time        	patch         	min           	max           	average       
0             	bottom	1.868099e+02	1.868099e+02	1.868099e+02
0             	sides	0.000000e+00	0.000000e+00	0.000000e+00
100           	bottom	1.598811e+02	2.174476e+02	1.701706e+02
100           	sides	2.192814e-02	1.108888e+02	2.684494e+01
200           	bottom	1.539062e+02	2.174475e+02	1.625827e+02
200           	sides	8.207370e-04	2.373603e+01	1.893065e+00
300           	bottom	1.537618e+02	2.174475e+02	1.625742e+02
300           	sides	4.922483e-03	8.914920e+00	7.251943e-01
400           	bottom	1.537618e+02	2.174475e+02	1.625742e+02
400           	sides	1.176434e-03	3.476801e+01	9.170413e-01
500           	bottom	1.537618e+02	2.174475e+02	1.625742e+02
500           	sides	3.049010e-03	8.458892e+00	7.410303e-01

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

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

coolhhh

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

尝试设置为simpleGrading (1 10.9 1)，网格总数不变，对应的网格和Yang2009中网格接近，第一层网格高度0.01m，增长率是1.06。计算的Y+数值在80左右，但计算结果还是跟上面的计算结果基本一样，k偏差仍较大。

这里更正，计算的Y+数值大概在160左右

李东岳

@coolhhh 非常感谢非常感谢！

目前看起来无解了。Yang2017的kOmegaSST我也试了，基本也是类似的结果。抓不到头脑。

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

是的，类似的我也尝试过。Yang2017里面给出了两种进口，一种log，一种幂指数。俩种进口型线一样。但是我画出来也不一样，区别很大。

真是想不出来为啥了。你在国内么，留个地址，给你弄两个CFD记事本过去。

coolhhh

@李东岳李老师，我尝试按照自己的理解，用fluent尝试算了Yang2009算例，参照OpenFOAM的设置，结果能对上，但和原论文结果有点区别：

按照Yang2009设置，顶面设置为slip wall，结果基本能自保持，但接近顶部的k无法自保持，猜测是顶面边界条件设置问题。

fluent中slip wall设置如下，把壁面剪切应力设置为0。

把顶面设置为symmetry，结果就基本一样。实在不知道为何OpenFOAM算出来k衰减那么多

coolhhh

似乎找到原因了，是要把算例0文件夹中的U、k、epsilon中用codedFixedValue编写的边界name命名为不同名称，基于李老师算例文件修改后如下，用OpenFOAM-9直接Allrun
Yang2009OFCase.tar

    inlet
    {
        type            codedFixedValue;
        value           uniform 0.1; //default value
        name            linearTBC1; //name of new BC type

修改后的算例结果如下，近壁面有点偏差，但总体是保持住了