中性大气环境湍流动能的自保持 | 附有算例下载
-
我最近参考了杨老师的这一篇行业标杆文章
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(右图中的红线,跟黑线有偏移)。
我个人怀疑是nut的壁面函数的因素。但是杨老师文章里面说用粗糙壁面函数就可以。想问问相关大佬,是不是我还忽略了什么?
更新:算例下载在这里 https://www.cfd-china.com/topic/6231/中性大气环境湍流动能的自保持/8
-
-
@李东岳 李老师您好,我们这边之前在fluent里面有验证过该入流条件,具体结果可参见刘安彬的硕士论文《平衡大气边界层的CFD数值模拟》。但是在模拟时上部流场的湍动能并无明显衰减,只是近壁面有衰减,并在施加壁面剪应力之后可有效改善。您可以考虑下顶面边界条件、壁面剪应力这两方面的设置。
下图模拟结果截取自上述硕士论文,未加任何修正
下图是施加壁面剪应力后的模拟结果
-
@Hope 我看你这面的常数跟杨老师那面不一样。我在主贴中更新了杨老师的这几个参数。是说这几个参数要调整么?你们这面的这几个参数是哪里得来的?
您可以考虑下顶面边界条件、壁面剪应力这两方面的设置。
yang2009这篇文章引用了近400,按道理来说应该可以通过这一篇sci复现相关的结果。如果考虑这些影响就可以正确的话,sci里面应该有提及。我在想是不是某些因素我没关注到,忽略掉了。比如这两个,在sci里面有隐含的提及么
-
@李东岳 李老师,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保持住。 -
@李东岳 李老师你好,请问你这个保持住是指多大的距离,你在图中标的inlet、middle、outlet是多大的计算域?
-
@疏影横斜水清浅 进口处,2500米处,出口处。一共5000米
-
@李东岳 好的,了解了,谢谢李老师
-
@李东岳 在 中性大气环境湍流动能的自保持 中说:
@疏影横斜水清浅 进口处,2500米处,出口处。一共5000米
李老师,看杨老师论文中计算域尺寸比较小,与同济风洞实验室尺寸一致。有无可能是算的全尺寸,网格太粗糙导致结果不同
-
能不能自保持,是要看壁面函数和入口能否协调的。我是用LES做大气边界层的。LES里跑预前模拟生成的流场就能让湍动能处处保持平衡。LES这边还有用入口生成方法的,但一般自保持不了,流体跑着跑着湍动能就会衰减,是因为壁面用的wall model和入口人工合成的流场不协调。
-
@coolhhh 我Yang2009算的是Yang2009的算例,HW07算的是HW07的算例。我暂时还没有尝试用hw07的方法算yang2009,以及用yang2009的方法算hw07,这个过几天试一下...
-
@Vortex
感谢分享,目前还没到LES那一阶段 :-) 下一阶段要上kOmega以及SST。LES估计要几个月后了。
-
@李东岳 介绍个专门做ABL的code,叫LESGO (https://lesgo.me.jhu.edu/),Meneveau大佬以前是拿它研究LES模型的。不稳定、中性、稳定的ABL都可以做,还能考虑科氏力。
-
@Vortex 在 中性大气环境湍流动能的自保持 中说:
@李东岳 介绍个专门做ABL的code,叫LESGO (https://lesgo.me.jhu.edu/),Meneveau大佬以前是拿它研究LES模型的。不稳定、中性、稳定的ABL都可以做,还能考虑科氏力。
想问一下LESGO有现成的考虑浮力驱动的求解器吗?
-
@chszkc 有,解温度输运方程,求解器有植入。
-
李东岳
-
-
@李东岳 收到,谢谢老师
-
@李东岳 李老师,算了您给的算例,计算结果跟您的结果一样,还没找到原因。但有以下问题:
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 在 中性大气环境湍流动能的自保持 | 附有算例下载 中说:
尝试设置为
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记事本过去。
-
@李东岳 李老师,我尝试按照自己的理解,用fluent尝试算了Yang2009算例,参照OpenFOAM的设置,结果能对上,但和原论文结果有点区别:
- 按照Yang2009设置,顶面设置为
slip wall
,结果基本能自保持,但接近顶部的k
无法自保持,猜测是顶面边界条件设置问题。
fluent中slip wall
设置如下,把壁面剪切应力设置为0。
- 把顶面设置为
symmetry
,结果就基本一样。实在不知道为何OpenFOAM算出来k
衰减那么多
- 按照Yang2009设置,顶面设置为
-
似乎找到原因了,是要把算例0文件夹中的
U
、k
、epsilon
中用codedFixedValue
编写的边界name
命名为不同名称,基于李老师算例文件修改后如下,用OpenFOAM-9直接Allrun
Yang2009OFCase.tarinlet { type codedFixedValue; value uniform 0.1; //default value name linearTBC1; //name of new BC type
修改后的算例结果如下,近壁面有点偏差,但总体是保持住了
-
完美解决!
太感谢了!这个小bug让我卡了一个月。太无语了
。。。。。
-
@coolhhh 在 中性大气环境湍流动能的自保持 | 附有算例下载 中说:
把顶面设置为symmetry,结果就基本一样
另外,symmetry与slip的差异极小。在openfoam这面,symmetry是一种几何性的限定,所有的场各种场都会严格的对称。但是slip只是速度会这样,其他场不会,比如
HbyA
。但是这个是fluent的结果,那fluent为什么这样,就不知道了。 -
@李东岳 或许是我这边fluent有某个设置错了,杨老师论文中顶面用的
slip
结果不会出现k
衰减。我再试了fluent中顶面设置no-slip
,结果是比较符合预期,接近顶面壁面平均速度减小,湍动能增大
-
李东岳
-
@李东岳 李老师,我最近尝试复现Yang2017的kOmegaSST,按照文献的设置和参数设置,k都无法保持稳定。不知道老师您是否复现出来了。
-
@SHUKK kOmegaSST没有尝试。你可以把你的算例打包发上来我给你看看
-
@李东岳 好的,谢谢李老师!y2017kOmegaSST.zip