中性大气环境湍流动能的自保持 | 附有算例下载
-
我最近参考了杨老师的这一篇行业标杆文章
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的壁面函数的因素。但是杨老师文章里面说用粗糙壁面函数就可以。想问问相关大佬,是不是我还忽略了什么?
更新:算例下载在第24楼
-
-
@李东岳 李老师,yang2009这篇sci论文里面我没看到有壁面剪应力的相关设置。由于这部分内容是我师兄做的,我先依据硕士论文解释您提出的几个问题,具体的我还需要再查证之后再给您回复。
硕士论文中引用的入流条件也是出自这篇sci文献,有出入的是壁面函数的粗糙高度取值,硕士论文中的粗糙高度是依据文献(Blocken et al. 2007, CFD simulation of the atmospheric boundary layer: wall function problems)计算所得:
壁面剪应力的取值是依据硕士论文中的一些假定提出的,而非yang2009中相关参数。 -
@Vortex 在 中性大气环境湍流动能的自保持 中说:
@李东岳 介绍个专门做ABL的code,叫LESGO (https://lesgo.me.jhu.edu/),Meneveau大佬以前是拿它研究LES模型的。不稳定、中性、稳定的ABL都可以做,还能考虑科氏力。
想问一下LESGO有现成的考虑浮力驱动的求解器吗?
-
-
@李东岳 李老师,算了您给的算例,计算结果跟您的结果一样,还没找到原因。但有以下问题:
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结果看起来像是一样的
-
-
-
@李东岳 测试的网格如下,网格尺寸比较大
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 ( );
-
@李东岳 李老师,我重新看了下原论文,做了些测试:
1. Yang2009论文中
k
剖面是缩尺的风洞拟合结果,不能用于全尺寸。Yang2009论文中k剖面,随着高度在增加如下所示。
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的方法,能否直接用于二维网格还没有验证,并且目前对比的结果都是取计算域中间的结果,不考虑边界的影响。
@疏影横斜水清浅 在 中性大气环境湍流动能的自保持 | 附有算例下载 中说:
@李东岳 测试的网格如下,网格尺寸比较大
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 ( );
还是有挺多疑问,请各位大佬解惑