CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

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

    OpenFOAM
    7
    30
    885
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • 李东岳
      李东岳 管理员 最后由 李东岳 编辑

      我最近参考了杨老师的这一篇行业标杆文章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

      CFD高性能服务器 http://dyfluid.com/servers.html

      H 1 条回复 最后回复 回复 引用
      • H
        Hope @李东岳 最后由 编辑

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

        博士生涯收敛

        1 条回复 最后回复 回复 引用
        • 李东岳
          李东岳 管理员 最后由 编辑

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

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

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

          CFD高性能服务器 http://dyfluid.com/servers.html

          H 1 条回复 最后回复 回复 引用
          • H
            Hope @李东岳 最后由 编辑

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

            博士生涯收敛

            李东岳 1 条回复 最后回复 回复 引用
            • 李东岳
              李东岳 管理员 @Hope 最后由 编辑

              @Hope 多谢多谢!我也研究研究。:duang:

              CFD高性能服务器 http://dyfluid.com/servers.html

              1 条回复 最后回复 回复 引用
              • 李东岳
                李东岳 管理员 最后由 编辑

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

                CFD高性能服务器 http://dyfluid.com/servers.html

                1 条回复 最后回复 回复 引用
                • 李东岳
                  李东岳 管理员 最后由 编辑

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

                  捕获.PNG

                  CFD高性能服务器 http://dyfluid.com/servers.html

                  疏影横斜水清浅 1 条回复 最后回复 回复 引用
                  • 疏影横斜水清浅
                    疏影横斜水清浅 @李东岳 最后由 编辑

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

                    李东岳 1 条回复 最后回复 回复 引用
                    • 李东岳
                      李东岳 管理员 @疏影横斜水清浅 最后由 编辑

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

                      CFD高性能服务器 http://dyfluid.com/servers.html

                      疏影横斜水清浅 C 2 条回复 最后回复 回复 引用
                      • 疏影横斜水清浅
                        疏影横斜水清浅 @李东岳 最后由 编辑

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

                        1 条回复 最后回复 回复 引用
                        • C
                          coolhhh @李东岳 最后由 编辑

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

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

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

                          李东岳 1 条回复 最后回复 回复 引用
                          • V
                            Vortex 最后由 编辑

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

                            李东岳 1 条回复 最后回复 回复 引用
                            • 李东岳
                              李东岳 管理员 @coolhhh 最后由 编辑

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

                              CFD高性能服务器 http://dyfluid.com/servers.html

                              1 条回复 最后回复 回复 引用
                              • 李东岳
                                李东岳 管理员 @Vortex 最后由 编辑

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

                                CFD高性能服务器 http://dyfluid.com/servers.html

                                V 1 条回复 最后回复 回复 引用
                                • V
                                  Vortex @李东岳 最后由 编辑

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

                                  C 1 条回复 最后回复 回复 引用
                                  • C
                                    chszkc @Vortex 最后由 编辑

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

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

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

                                    V 1 条回复 最后回复 回复 引用
                                    • V
                                      Vortex @chszkc 最后由 编辑

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

                                      1 条回复 最后回复 回复 引用
                                      • Referenced by  李东岳 李东岳 
                                      • 李东岳
                                        李东岳 管理员 最后由 李东岳 编辑

                                        Yang2009.tar.xz

                                        算例在这里

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

                                        @coolhhh 直接给你发这了

                                        CFD高性能服务器 http://dyfluid.com/servers.html

                                        C 2 条回复 最后回复 回复 引用
                                        • C
                                          coolhhh @李东岳 最后由 编辑

                                          @李东岳 收到,谢谢老师

                                          1 条回复 最后回复 回复 引用
                                          • C
                                            coolhhh @李东岳 最后由 编辑

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

                                            1. 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)
                                                    );
                                                }
                                            );
                                            
                                            

                                            1. 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
                                            

                                            1. 尝试设置为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
                                            

                                            1. Yang2009中,提到了SKE1和SKE2两种模型,并模拟对比。
                                              c59352fa-c8e2-4ea4-af17-011af48def7f-image.png
                                              9d478c4c-9e46-4643-b4c2-bd5c9ec76d3d-image.png
                                              尝试根据式(25)和(26)把这两个k曲线画出对比,这两条曲线还是有挺大区别,但在原论文中两个工况的inlet结果看起来像是一样的
                                              k剖面对比.png
                                            C 1 条回复 最后回复 回复 引用
                                            • C
                                              coolhhh @coolhhh 最后由 编辑

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

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

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

                                              李东岳 1 条回复 最后回复 回复 引用
                                              • 李东岳
                                                李东岳 管理员 @coolhhh 最后由 李东岳 编辑

                                                @coolhhh 非常感谢非常感谢!

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

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

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

                                                :zoule:

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

                                                CFD高性能服务器 http://dyfluid.com/servers.html

                                                C S 2 条回复 最后回复 回复 引用
                                                • C
                                                  coolhhh @李东岳 最后由 编辑

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


                                                  1. 按照Yang2009设置,顶面设置为slip wall,结果基本能自保持,但接近顶部的k无法自保持,猜测是顶面边界条件设置问题。
                                                    Ux剖面.png TKE剖面.png
                                                    fluent中slip wall设置如下,把壁面剪切应力设置为0。
                                                    43329f4b-3ec9-42fb-bf45-7b1e48e79db2-image.png

                                                  1. 把顶面设置为symmetry,结果就基本一样。实在不知道为何OpenFOAM算出来k衰减那么多
                                                    Ux剖面.png TKE剖面.png
                                                  C 1 条回复 最后回复 回复 引用
                                                  • C
                                                    coolhhh @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
                                                    

                                                    修改后的算例结果如下,近壁面有点偏差,但总体是保持住了
                                                    95463c69-25ce-4003-89e4-7e5b8db5db1e-image.png
                                                    0bf33617-11ed-4328-ac7a-9995380adc5f-image.png

                                                    1 条回复 最后回复 回复 引用
                                                    • 李东岳
                                                      李东岳 管理员 最后由 编辑

                                                      完美解决!

                                                      太感谢了!这个小bug让我卡了一个月。太无语了

                                                      。。。。。

                                                      CFD高性能服务器 http://dyfluid.com/servers.html

                                                      1 条回复 最后回复 回复 引用
                                                      • 李东岳
                                                        李东岳 管理员 最后由 编辑

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

                                                        把顶面设置为symmetry,结果就基本一样

                                                        另外,symmetry与slip的差异极小。在openfoam这面,symmetry是一种几何性的限定,所有的场各种场都会严格的对称。但是slip只是速度会这样,其他场不会,比如HbyA。但是这个是fluent的结果,那fluent为什么这样,就不知道了。

                                                        CFD高性能服务器 http://dyfluid.com/servers.html

                                                        C 1 条回复 最后回复 回复 引用
                                                        • C
                                                          coolhhh @李东岳 最后由 编辑

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

                                                          1 条回复 最后回复 回复 引用
                                                          • Referenced by  李东岳 李东岳 
                                                          • S
                                                            SHUKK @李东岳 最后由 编辑

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

                                                            李东岳 1 条回复 最后回复 回复 引用
                                                            • 李东岳
                                                              李东岳 管理员 @SHUKK 最后由 编辑

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

                                                              CFD高性能服务器 http://dyfluid.com/servers.html

                                                              S 1 条回复 最后回复 回复 引用
                                                              • S
                                                                SHUKK @李东岳 最后由 编辑

                                                                @李东岳 好的,谢谢李老师!y2017kOmegaSST.zip

                                                                1 条回复 最后回复 回复 引用
                                                                • First post
                                                                  Last post