gauss upwind和gauss linearUpwind grad(U)结果差异大



  • 在计算建筑物外流场时,使用的是simpleFoam稳态求解器,我首先使用了gauss upwind作为div项的离散格式,因为我发现这个格式计算的快一些。
    在计算稳定后,提取了屋面的摩擦速度,(此时计算已经稳定了),然后我把gauss upwind格式改为gauss linearUpwind grad(U),其他设置不动,接着这个时间步进行往下计算,算了一段时间后,再提取屋面摩擦速度,发现gauss linearUpwind grad(U)的摩擦速度特别怪异,请大家帮忙分析一下原因。。
    屏幕截图 2021-01-29 143832.jpg
    麻烦点在于,我看论文一般都用的二阶格式,但是我的这个二阶格式算出来太怪异了,效果还有没有一阶的gauss upwind好,论文里不敢说自己用的一阶的格式。。


  • 教授

    @Samuel-Tu
    绝大部分情况下,动量方程对流项离散用 gauss upwind 会导致很大的误差,不建议使用。
    gauss linearUpwind grad(U) 算出来有震荡,可以试试对 grad(U) 使用 gradient limiter,应该会有改善。不知道你现在 gradSchemes 是怎么设置的?



  • @xpqiu 感谢,又学到一个一阶迎风的知识。我的gradSchemes是采用的Gauss linear,laplacianSchemes是用的Gauss linear corrected。我去查一下怎么加gradient limiter



  • 太奇怪了。。。fluent里面用二阶迎风格式:
    0fbba8ab-31a5-4acf-a1de-c5583adc9420-image.png
    算出来没问题。但是OpenFOAM里面对div(phi,U)项只能用Gauss upwind才能和fluent里的结果吻合,如果对div(phi,U)项采用高阶项,例如Gauss linearUpwind grad(U),对grad(U)进行限制、leastSquares,cellLimited,结果全都对不上。。太奇怪了。。。。
    湍流的div项我是用的二阶迎风,改了格式影响也很小。。这不知道怎么解释了。。。论文里确实没看见过对div(phi,U)用一阶迎风的。。。


  • 教授

    @Samuel-Tu fvSchemes 和 fvSolution 贴上来看看



  • @xpqiu
    fvSchemes:

    ddtSchemes
    {
        default         steadyState;
    }
    
    gradSchemes
    {
        default         Gauss linear;
    
        grad(U)         cellLimited leastSquares 1;//还用了cellLimited gauss linear和gauss linear
        grad(k)         cellLimited leastSquares 1;
        grad(epsilon)     cellLimited leastSquares 1;
    }
    
    divSchemes
    {
        default         none;
    
        div(phi,U)      Gauss limitedLinear 1;//还用了gauss upwind, gauss linearUpwind grad(U), bounded Gauss linearUpwind grad(U), 发现就gauss upwind和计算结果最吻合,gauss limitedLinear 1吻合的还行,其他的就完全有问题了。。
    
        div(phi,k)      bounded Gauss linearUpwind limited;//k和epsilon换了几个格式,发现影响不大
        div(phi,epsilon) bounded Gauss linearUpwind limited;
    
        div((nuEff*dev2(T(grad(U))))) Gauss linear;
    }
    
    laplacianSchemes
    {
        default         Gauss linear corrected;
    }
    
    interpolationSchemes
    {
        default         linear;
    }
    
    snGradSchemes
    {
        default         corrected;
    }
    
    wallDist
    {
        method meshWave;
    }
    

    fvSolution:

    solvers
    {
        p
        {
            //solver          PCG;
            //preconditioner  DIC;
            solver          GAMG;
            tolerance       1e-7;
            relTol          0.1;
            smoother        DICGaussSeidel;
        }
    
        "(U|k|epsilon|omega)"
        {
            solver          PBiCGStab;
            preconditioner  DILU;
            //solver          smoothSolver;
            //smoother        symGaussSeidel;
            tolerance       1e-7;
            relTol          0.001;
        }
    }
    
    SIMPLE
    {
        nNonOrthogonalCorrectors 2;
        consistent      yes;
    
        residualControl
        {
            p               1e-7;
            U               1e-7;
            "(k|epsilon|omega|f|v2)" 1e-7;
        }
    }
    
    /*relaxationFactors
    {
        equations
        {
            U               0.9;
            ".*"            0.8;
        }
    }*/
    relaxationFactors
    {
        fields
        {
            p               0.3;
        }
        equations
        {
            U               0.7;
            "(k|omega|epsilon).*" 0.7;
        }
    }
    

    我换了另一个小模型,算的快一些,结果如下(红线的数据不见了,但我记得形状):
    cca97ee0-5054-4553-95e5-98340b8608e8-image.png



  • wholeModel.png



  • 0f0345fa-d59b-4030-9578-0dcb60ab0d6f-image.png


  • 教授

    @Samuel-Tu
    算例如果可以分享的话,发给我玩玩?



  • @xpqiu 好滴,稍等我下载一下



  • @xpqiu 不知道咋回事,网站上传不了,分卷压缩也不行,我用奶牛快传了。
    这是地址,取件码:118173



  • @xpqiu
    链接:https://pan.baidu.com/s/1Ts5W8eHIITgI4bjOOHyS6g
    提取码:kk6h
    或者用这个百度网盘的,我用的simpleFoam,这里面多加了一个50时间步的计算结果,这是用gauss linear计算比较稳定后的结果,方便改了div schemes再计算,能更快的得到稳定后的结果。。


  • 教授

    @Samuel-Tu
    结果怎么对比呢,你有后处理脚本吗?



  • @xpqiu 我是用postProcess -func wallShearStress 和 writeCellCentres把roof这个patch 输出出来,再在paraview里面用calculater: 摩擦速度=sqrt(mag(wallShearStress)),最后在excel里画的roof里的每个面单元Cx(横坐标)和对应的摩擦速度的曲线。。我研究一下怎么整合出一个后处理脚本。。主要是postProcess似乎没有sqrt这个命令,如果有就可以不用paraview了。。


  • 教授

    @Samuel-Tu 明白了,我算几个看看



  • @xpqiu 麻烦了,我研究一下后处理能不能总结一下,,减少一下工作量。


  • 教授

    @Samuel-Tu
    8198a8ae-f0dd-4498-bb4b-6cb3e9b6732b-image.png

    这个是 div(phi,U) 我用 Gauss linearUpwind 的结果,看起来跟 fluent 的比较接近。



  • @xpqiu 能不能分享一下fvSchemes。。为什么我算出来不是这样的。。。我看看设置


  • 教授

    @Samuel-Tu
    你这个算例设置得不太合理的地方有好几个,我没有一个一个去考察每个因素的影响,只是直接改成了我认为合理的设置。我把我改过的文件打包了,如下:
    OFPlate0.zip

    另外,你的网格也不太好
    U_grid.png

    这个网格咋一看,直觉就是边界层附近网格太粗了,看了一下最后一步的 yPlus分布,只看 roof 这个面,

    yPlus.png

    果然,最小都300 多了。你用的湍流模型是 rke,这个是合理的选择,但是这个模型最好把 yPlus 控制在 30-100 范围内。
    最后,我把我这边 gauss upwind 和 gauss linearUpwind 的结果对比放上来,
    756323d2-5df0-414e-8ba2-61bd03b263f2-image.png

    虽然看起来可能还是 upwind 的结果在中间这一段跟文献的更接近,但是左边这一段实际上这两个结果都与文献值有明显的偏差。我觉得你还是先把网格改得合理之后再算算看,应该会有改进的。



  • @xpqiu 感谢了。。我先学习一下这个设置。。。网格这个我也很困惑,因为这个按照论文的最小网格尺寸设置的,我也发现了按论文的设置yPlus会很大。。但是作者做过网格无关性验证,尺寸改小后摩擦速度没有超过5%的变化,然后作者就这么用了,,后来 又来了两篇类似文献,全都是参考这个论文的设置,。。


  • 教授

    @Samuel-Tu
    论文描述的最小网格尺寸,可能是除了 prism layer 以外区域的最小尺寸啊,他有没有提他的prism layer是怎么生成的,first layer thickness,expansion ratio,number of layers 这些?



  • @xpqiu 论文只提了minimum size,他们用fluent做的,所以肯定没有用snappyHexMesh了,但是他们课题组发的其他文章差不多都是控制最小尺寸为0.1m左右。我估计因为他们是土木方向的,所以可能对壁面函数对yPlus的限制不太了解。。。


  • 教授

    @Samuel-Tu
    那么,你上面发的paper结果只是 fluent仿真结果是吧,有没有实验结果参照?



  • @xpqiu 确实是只有仿真结果,没有直接的试验结果参考。。但是这个仿真是为风吹雪这个物理过程做参考的,他们是做了风吹雪的风洞试验的,按照这个仿真结果算风吹雪的现象和风洞试验对上了。。。
    另外我看您的设置有两个问题:
    1.我看k文件里把壁面kLowREWallFucntion改成了kqrWallFunction。我记得kLowREWallFunction是适合低、高雷诺的,kqrWallFunction只适合高雷诺。我就是按比较保险的设置,所以用的kLowREWallFunction(当然现在看来yPLus那么大,没有必要用这个低、高雷诺都可以的壁面)。您这里为什么改成kqrWallFunction呢
    2.fvSolution里面把nNonOrthogonalCorrectors从2改成0了。我记得这个是好像用来修正网格非正交的,我的网格有非正交的网格,这里改动的目的是啥呢。。



  • cfdOnline里关于这个非正交修正次数的讨论
    Non-orthogonal correctorsare here to save you if your code is blowing up because the mesh is so non-orthogonal that the first solution is driving the velocity to be stupid. If your velocity is OK, you just keep doing "normal" correctors, without special need for non-orthogonal ones.

    I use them on bad meshes (some people call them "industrial") when the solver is giving me trouble. Usually, 1 is enough, and I never used more than 3.

    Hope this helps,

    Hrvoje
    按照hjasak的说法,nNonOrthogonalCorrectors实际不是修正网格非正交的,而是使得计算不容易发散。又学到了一点


  • 教授

    @Samuel-Tu
    kLowReWallFunction 这个是针对特定的湍流模型来使用的,不能用于 rke 模型。



  • @xpqiu 感谢,学到了


  • 副教授

    @xpqiu 想请教一下,OpenFOAM的壁面函数有没有很系统的文献啥的,我一直参考的是源代码和这个。有时候我会碰到同样的网格和同样的湍流模型,用某个壁面函数就可以,但是用另外一个壁面函数就会发散;或者是同样的网格和同样的壁面函数,一种湍流模型(kEpsilon和RNGkEpsilon)可以算,另一种(realizableKE)就会发散。感觉难点在于湍流量(k, epsilon, omega)和湍流粘度(nut)都需要壁面函数,nut按定义是根据湍流量算出来的,为什么还要给壁面函数呢?


  • 教授

    @cccrrryyy
    这个链接里面解释了为什么nut需要设置壁函数:https://www.slideshare.net/fumiyanozaki96/openfoam-36426892
    这几个视频里面解释了壁函数以及k,epsilon的壁面处理的原则:
    https://www.youtube.com/watch?v=fJDYtEGMgzs&t=4s
    https://www.youtube.com/watch?v=zSIgKsQSa9k
    https://www.youtube.com/watch?v=RKoXFpwi2go

    其他资料:
    On the Wall Boundary Condition for Turbulence Models , By JONAS BREDBERG

    Kalitzin, G., Medic, G., Iaccarino, G., & Durbin, P. (2005). Near-wall behavior of RANS turbulence models and implications for wall functions. Journal of Computational Physics, 204(1), 265–291.

    Popovac, M., & Hanjalic, K. (2007). Compound wall treatment for RANS computation of complex turbulent flows and heat transfer. Flow, Turbulence and Combustion, 78(2), 177–202.

    OpenFOAM里面的代码实现通常跟文献的描述有点差异,所以也只能是尽可能去了解原理,结合代码来理解,然后实践。也不是所有的出错都能找到合理的解释:xinlei:


  • 副教授

    好的,感谢,我先看看这些!

    确实,有时候我感觉参考理论指导就可以,但有时候光靠这个不行,和代码的植入可能也有一定关系。我以前用过很长时间FLUENT,它里面基本都是成套的,比如kEpsilon下面就有好像是4个壁面函数可以选择,并不是每个量都要给一个。另外它也没有提供nut的边界条件,不知道是没用还是用了些什么东西但是它不明说。



  • @cccrrryyy fluent里面壁面函数只用设置一下,而OF里所有湍流变量都要设置。。我都搞不明白什么时候该搭配什么。只是知道根据y+简单判断一下该用哪种。。


Log in to reply
 


CFD中文网 | 东岳流体学术 | 东岳流体商业 | 吉ICP备20003622号-1