求助: SonicFoam求解二维叶栅时无法收敛(求解文件可下载)



  • 大家好,最近一直在利用OpenFOAM的SonicFoam求解二维叶栅的跨声速流动,由于缺乏计算经验,计算一直无法收敛,现在通过这个帖子向大家求教,恳求大家帮忙指点。原文件如下:https://pan.baidu.com/s/1i_9ZtTNnw_uyRIB9WqPHeg 提取码:3vx0

    我选用的计算对象是rotor67的70%截面,预计进口Ma数为1.2,进出口边界条件如图所示,进口给定总温总压以及速度方向的进口条件,出口给定背压。因为一开始怕计算发散,采用先给定小进口总压,然后逐步抬高,进口总压为1.2E5时的进口速度大约是180m/s。
    c2a0aee9-84a6-4446-aad8-fca41215625f-图片.png 进出口边界示意图

    网格采用ICEM生成,yplus是在来流速度400左右确定为30,,第一层网格大约为0.03mm。

    该网格可以在rhoCentralFoam求解器中采用速度入口条件可以收敛,所以网格和边界条件这一块应该没有什么问题。原先选用其结果作为初场,可是计算也是发散。

    我这边一直怀疑是自己的边界和求解控制那一块有问题,但是一直没能找到问题解决办法,恳请大家帮忙提些建议。

    这里给出我的压力、温度、速度条件

    pOut            1e5;
    pAll            1.2e5;
        INLET
        {
            type            totalPressure;
            psi             thermo:psi;  
            gamma           1.4;      
            p0              uniform $pAll;
            value           uniform $pOut;
        }
        OUTLET
        {
            type            waveTransmissive;
            field           p;
            psi             thermo:psi;
            gamma           1.4;
            fieldInf        100000;
            lInf            0.3;
            value           uniform 100000;
        }
        PER1
        {
            type            cyclicAMI;
        }
        PER2
        {
            type            cyclicAMI;
    
        }
        BLADE
        {
            type             zeroGradient;
        }
    
    Tinlet
    INLET
        {
            type       totalTemperature;
            gamma      1.4;
            psi        thermo:psi;
            T0         uniform 304.975;
            value      uniform 288.15;
        }
        OUTLET
        {
            type            inletOutlet;
            inletValue      uniform $Tinlet;
            value           $inletValue;
        }
        PER1
        {
            type            cyclicAMI;
        }
        PER2
        {
            type            cyclicAMI;
    
        }
        BLADE
        {
            type            zeroGradient;
        }
    
    boundaryField
    {
        INLET
        {
            type            pressureDirectedInletVelocity;
            inletDirection  uniform (0.5 0.86605 0);    //angle=60
            value           uniform (90 159.35 0);
        }
        OUTLET
        {
            type            inletOutlet;
            inletValue      uniform $Uinlet;
            value           uniform $Uinlet;
        }
        PER1
        {
            type            cyclicAMI;
        }
        PER2
        {
            type            cyclicAMI;
        }
        BLADE
        {
            type            noSlip;
        }
    
    

    这里给出fvSolution和fvSchemes。

    solvers
    {
        "rho.*"
        {
            solver          diagonal;
        }
    
        "p.*"
        {
            solver          smoothSolver;
            smoother        symGaussSeidel;
            tolerance       1e-08;
            relTol          0;
        }
    
        "(U|e|R).*"
        {
            $p;
            tolerance       1e-05;
        }
    
        "(k|omega).*"
        {
            $p;
            tolerance       1e-08;
        }
    }
    
    PIMPLE
    {
        nOuterCorrectors 2;
        nCorrectors      1;
        nNonOrthogonalCorrectors 2;
    }
    fieldBounds
    {
    
        p      1000     1e6;
        T      100    3000;
        U      1000;
    }
    
    relaxationFactors
    {
        fields
        {
            p 0.7;
            rho 0.01;
        }
        equations
        {
            p 0.7;
            U 0.3;
            "(e|h|k|epsilon|omega)" 0.2;
        }
    }
    
    ddtSchemes
    {
        default         Euler;
    }
    
    gradSchemes
    {
        default         Gauss linear;
    }
    
    divSchemes
    {
        default         none;
        div(phi,U)      Gauss limitedLinearV 1;
        div(phi,e)      Gauss limitedLinear 1;
        div(phid,p)     Gauss limitedLinear 1;
        div(phi,K)      Gauss limitedLinear 1;
        div(phiv,p)     Gauss limitedLinear 1;
        div(phi,k)      Gauss upwind;
        //div(phi,epsilon) Gauss upwind;
        //div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
        div(phi,omega)  Gauss upwind;
        div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
        div((muEff*dev2(grad(U).T())))  Gauss linear; 
    }
    
    laplacianSchemes
    {
        default         Gauss linear corrected;
    }
    
    interpolationSchemes
    {
        default         linear;
    }
    
    snGradSchemes
    {
        default        corrected;
    }
    
    wallDist
    {
        method meshWave;
        //nRequired       false;
    }
    


  • 不收敛可能性太多了,周一给你调试下你的算例,现在手头没有openfoam

    既然rhoCentralFoam能收敛,为何不用rhoCentralFoam



  • @东岳 真的是太感谢您了,rhoCentralFoam的总温总压边界再配上进口速度方向,会使得计算发散,我现在正在想办法解决。在压气机内流叶栅的计算里面,一般都会给进口总温总压边界条件,和实验对比的比较好



  • 源文件里面的松弛因子给的不合理,修改以后依然没什么效果
    修改后

    relaxationFactors
    {
        fields
        {
            p 0.7;
            rho 0.8;
        }
        equations
        {
            p 0.7;
            U 0.3;
            "(e|h|k|epsilon|omega)" 0.2;
        }
    


  • 你的算例我下来之后,发现lowRe这个计算没有问题?



  • @东岳 计算是要过一段时间才会显示发散,看云计算流场应该是乱的



  • @东岳 这是5e-5时刻的流场云图
    b4b29278-89ba-421f-83e6-acef2c2935d0-图片.png



  • 现在我能确定的是我的整个边界条件在低速下是可以用的,用pisofoam,湍流模型选用kepsilon收敛的很好,现在正在尝试kOmega。但是在高速求解器如我现在要用的sonicfoam 和原先尝试的rhoCentralFoam都会出现问题。在sonicfoam 里面的问题就是会出现内部流场十分紊乱,所以我觉得是不是我求解和算法控制那块是不是有问题



  • @宝丁 这个看起来出现了震荡,我在测试,就是算起来太慢了,



  • @宝丁 预计进口流速为60m/s,还是总温总压入口



  • @东岳 我在这里找到了一个您的邮箱,给您发了一个这个网格量一般的网格,您能收到吗?
    振荡是从一开始就振荡。时间步长为1e-8或者5e-9时,计算到0.0001就可以看到不稳定的流场。



  • 我把你的算例,把速度进口替换成zeroGradient计算就没问题了,你试试

    但你这个是不是应该是cyclic边界?另外,即使用cyclicAMI,是不是应该是rotational的?
    另外,网格正交性是不是可以调整下。个人感觉你调试这个算例,不用满足y+,用3000 5000调试就行, 要不算起来太慢了



  • 我也尝试过,确实是可以的,但是我的工况要求我定义进口的方向,尤其是这种跨声速流场,对方向很敏感



  • 我的边界条件是平移周期边界,cyclicAMI和cyclic相比只是不需要交界面两边网格一一对应

    ed798635-7613-4f0f-a894-34b0618270f0-1490441868405-cyclic-resized.jpg



  • 把你的tin和blk发来看看,我调整下网格,



  • @东岳 当时试的时候一开始确实没什么问题,可是算到后面zeroGradient 边界条件也会发散,但是会好很多,不知道是不是几何不合适,就没有往下尝试



  • @东岳 好的 文件不大,就直接放到这里了
    R67_70.zip



  • 速度角度和网格不成直线,对流项格式不好设置。不太清楚你这个结果应该是什么样,最后我这个稳定到这样了。因为我重新画了网格,你用你的网格试试下面的fvScheme,我感觉可能是格式的问题,因为看起来是震荡。

    还得进一步研究一下,不能用这么迎风的格式,你用稳态求解器试过没

    捕获.JPG

    ddtSchemes
    {
        default         Euler;
    }
    
    gradSchemes
    {
        default         cellLimited leastSquares 01;
    }
    
    divSchemes
    {
        default         none;
    
        div(phi,U)      Gauss upwind;
        div(phi,e)      Gauss upwind; 
        div(phid,p)     Gauss  upwind; 
        div(phi,K)      Gauss  upwind; 
        div(phiv,p)     Gauss  upwind; 
    
    
        //div(phi,U)      Gauss limitedLinearV 1;
        //div(phi,e)      Gauss limitedLinear 1;
        //div(phid,p)     Gauss limitedLinear 1;
        //div(phi,K)      Gauss limitedLinear 1;
        //div(phiv,p)     Gauss limitedLinear 1;
        div(phi,k)      Gauss upwind;
        div(phi,omega) Gauss upwind;
        div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
    }
    
    laplacianSchemes
    {
        default         Gauss linear limited corrected 0.5;
    }
    
    interpolationSchemes
    {
        default         linear;
    }
    
    snGradSchemes
    {
        default         corrected;
    }
    
    wallDist
    { 
        method          meshWave;
    }
    


  • @东岳 好的,非常非常感谢,我马上就去尝试



  • @东岳 我用稳态算过,也是发散



  • @东岳 这个是用rhoCentralFoam 采用速度入口算出来的结果,我感觉挺好的

    87c00dc6-c017-4e3b-b02f-ef2dba50ac4e-5646029ed87f78ea52a2af3fe421329.png



  • @东岳 这个是用您的对流项格用比较密的网格算出来的结果,稳定了,但是求解的流场还有一些不太理想,我再好好调一下,感谢东岳老师。
    b27179dd-fd4d-4a28-9564-95ac0a98237f-图片.png
    57e5ad47-6366-4e0d-9ffe-220de175ccc4-图片.png



  • @东岳 东岳老师,感觉问题并没有解决,我用upwind的算法求解稀疏一些的网格算的久一些就会出现发散的情况,网格量大的却能收敛,收敛后的流场如上所示也会有一些不好的地方,现在尝试提高离散项精度,但是很难找到合适的离散项,有些不知所措了,产生这个的根本原因是什么。



  • cascade.zip

    这是我最后得到的结果,不知道准不准,从下面那个图看起来差不多,你可以顺着我的这个算例,逐个调整参数直到准确位置

    不要运行blockMesh,忘记删了

    Screenshot from 2019-04-17 09-32-16.png



  • @东岳 好的,谢谢老师



  • @东岳 老师,这个算例运行的时候残差很高,然后到0.01s 左右就会发散了,不知道是不是网格的问题



  • @宝丁 降低时间步长试试



  • @东岳 好的



  • @东岳 @东岳 好像也没有作用,我也用我的网格试了一下,只改网格,出口那个地方也是出现了不稳定现象,但是还能算,感觉收敛不了



  • @东岳 可以尝试把速度 压力 温度的离散格式全改成二阶迎风,目前看效果挺好的



  • @宝丁 多谢分享 楼主好人 :xiexie:



  • @东岳 额 还没有算收敛,这个好难啊



  • 用rhoPimpleFoam试试?



  • @东岳 好的



  • @东岳 老师,不过我在用rhoSimpleFoam中遇到了一个很奇怪的现象,就是进口的压力会出现异常,使得整个计算发散,就像这样。这又会出现数值上面的异常。可以用层流模型试一下。这个问题也是困扰了我好久

    e06b455e-eafa-4128-ac62-adf23fe1b0ad-图片.png



  • @东岳 不对,是rhoCentralFoam



  • @东岳 老师,我已经调通了,非常感谢老师的知道。具体做法就是只修改

        default         none;
       Gauss limitedUpwind limited
        div(phi,e)      Gauss limitedUpwind limited;
        div(phid,p)     Gauss limitedUpwind limited;
        div(phi,K)      Gauss limitedUpwind limited;
        div(phiv,p)     Gauss limitedUpwind limited;
        div(phi,k)      Gauss upwind;
         div(phi,omega)  Gauss upwind;
        div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
        div((muEff*dev2(grad(U).T())))  Gauss linear; 
    

    注意,现在计算出来的结果并不是很理想,感觉粘性太大了,所以想问一下东岳老师

    Gauss limitedUpwind limited
    
     wallDist  {method Poisson ;}
    

    这两项的含义
    附一张计算的结果
    eadabeaa-0e5f-4865-8b4a-e3785c391d20-cda3d1260df247c4f73e55f20ce2f7b.png



  • 非常非常非常感谢东岳老师的帮助



  • Gauss linearUpwind limited:

    • Gauss: 高斯积分计算散度

    • linearUpwind: 线性迎风格式

    • limited: 在线性迎风格式需要指定U的梯度怎么计算,limited用来调用U的梯度计算格式

    wallDist  {method Poisson ;}
    

    是用来计算壁面距离的,Poisson是一种计算避免距离的方法,比较适用于复杂几何,meshWave方法适用于简单几何。kOmegaSST模型会用到壁面距离。



  • @东岳 哦哦,不知道怎么复制错了,是那个midmod:wocao:



  • @东岳 还有,老师,您这些都是在哪里看的,感觉从来没见过的用法



  • 老铁,Minmod是一个特别出名的TVD格式啊。这些都是在书上看到的

    :qichuang:

    捕获.JPG



  • @东岳 好的,谢谢老师,做这个的时间太短了,后面一定加强总结



  • @东岳 真心求教这是什么书:shangxue:



  • Finite Volume Methods for Hyperbolic Problems

    • 对于刚学CFD的,我不推荐这本书,这本书单纯讲求解$\frac{\p \phi}{\p t}+\nabla\cdot(\bfU\phi)=0$,并且略针对可压缩空气动力学

    • 学CFD几年的可以看看,这本书讲解的很细致思路很独特,扩展性很强



  • @东岳 收到



  • @东岳 老师,国内有没有比较好的讲这种算法的书,我的英语水平如果要看这样的书花费的时间太长了



  • 国内有没有比较好的讲这种算法的书

    @宝丁 那本书即使英语水平好,看起来也费尽,哈哈 :zoule: 目前中文版的我只看过《数值传热学》,但《数值传热学》没有介绍这部分的内容。目前没有其他中文书推荐了。博士的话得加强英语啊。现在好多学生写SCI,我要是改的话基本等于重写。

    莫非是硕士?硕士毕业就无所谓了,怎么快怎么来



  • @东岳 我是硕士,入门不久,但是也想多学点东西。现在对老师真是钦佩的不行。



  • @宝丁 要读博士的话,趁硕士期间多看CFD书,等到博士中后期,好多时间都要用在发SCI写论文上了。现在打基础完全来得及。我当时硕士期间研究方向没定,但比较喜欢推方程,自己硬看了2年CFD教材和C++。CFD这面东西太多了,算法太多,感觉不看个几十本大厚书搞不出个全貌..硕士时间太紧了

    不读博士的话,硕士这样直接算毕业完全也没问题。要多学点东西的话,看《数值传热学》,CFD的基本流程也能有点感觉


Log in to reply