CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

    计算钝体涡激振动时发散(pimpleDyMFoam)

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

      各位老师好,我使用pimpleDyMFoam计算钝体涡激振动时,算例的CFL越来越大并计算一段时间后发散。通过查看速度和压力场,发现一开始就不太正常,然后随着时间推移,误差逐渐积累最后发散。报错如下:

      #0  Foam::error::printStack(Foam::Ostream&) at ??:?
      #1  Foam::sigFpe::sigHandler(int) at ??:?
      #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
      #3  Foam::GAMGSolver::scale(Foam::Field<double>&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const at ??:?
      #4  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
      #5  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
      #6  Foam::fvMatrix<double>::solveSegregated(Foam::dictionary const&) at ??:?
      #7  Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/home/ying/OpenFOAM41/OpenFOAM-4.1/platforms/linux64GccDPInt32Opt/bin/pimpleDyMFoam"
      #8  ? in "/home/ying/OpenFOAM41/OpenFOAM-4.1/platforms/linux64GccDPInt32Opt/bin/pimpleDyMFoam"
      #9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
      #10  ? in "/home/ying/OpenFOAM41/OpenFOAM-4.1/platforms/linux64GccDPInt32Opt/bin/pimpleDyMFoam"
      Floating point exception (core dumped)
      

      我查了一下有人说很可能是网格或边界条件的问题。我的边界条件是比较基础的速度入口/压力出口,感觉应该没问题;我也使用了refine的网格重新计算,发现crash的时间延长了,但速度场和压力场还是一开始就不正确。两种网格质量都没问题。

      我附加了我的算例设置文件,和流场图以供参考。

      1. 粗网格及流场
        coarse.png
        coarse_U.png
        coarse_p.png

      2. 细网格及流场
        refine.png
        refine_U.png
        refine_p.png

      3. 边界条件(U/p/pointDisplacement)

      FoamFile
      {
          version     2.0;
          format      ascii;
          class       volVectorField;
          location    "0";
          object      U;
      }
      // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      dimensions      [0 1 -1 0 0 0 0];
      
      internalField   uniform (0 0 0);
      
      boundaryField
      {
          INLET
          {
              type            fixedValue;
              value           uniform (1 0 0);
          }
          OUTLET
          {
              type            zeroGradient;
          }
          TOP
          {
              type            symmetry;
          }
          BOTTOM
          {
              type            symmetry;
          }
          CYLINDER
          {
              type            movingWallVelocity;
              value           uniform (0 0 0);
          }
          frontAndBackPlanes
          {
              type            empty;
          }
      }
      
      FoamFile
      {
          version     2.0;
          format      ascii;
          class       volScalarField;
          location    "0";
          object      p;
      }
      // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      
      dimensions      [0 2 -2 0 0 0 0];
      
      internalField   uniform 0;
      
      boundaryField
      {
          INLET
          {
              type            zeroGradient;
          }
          OUTLET
          {
              type            fixedValue;
              value           uniform 0;
          }
          TOP
          {
              type            symmetry;
          }
          BOTTOM
          {
              type            symmetry;
          }
          CYLINDER
          {
              type            zeroGradient;
          }
          frontAndBackPlanes
          {
              type            empty;
          }
      }
      
      FoamFile
      {
          version     2.0;
          format      ascii;
          class       pointVectorField;
          location    "0";
          object      pointDisplacement;
      }
      // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      
      dimensions      [0 1 0 0 0 0 0];
      
      internalField   uniform (0 0 0);
      
      boundaryField
      {
          INLET
          {
              type            fixedValue;
              value           uniform (0 0 0);
          }
          OUTLET
          {
              type            fixedValue;
              value           uniform (0 0 0);
          }
          TOP
          {
              type            symmetry;
          }
          BOTTOM
          {
              type            symmetry;
          }
          CYLINDER
          {
              type            calculated;
          }
          frontAndBackPlanes
          {
              type            empty;
          }
      }
      
      1. fvScheme & fvSolution
      FoamFile
      {
          version     2.0;
          format      ascii;
          class       dictionary;
          object      fvSchemes;
      }
      // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      ddtSchemes
      {
          default         backward;
      }
      
      gradSchemes
      {
          default         Gauss linear;
          grad(p)         Gauss linear;
          grad(U)         Gauss linear;
      }
      
      divSchemes
      {
          default         none;
          div(phi,U)      Gauss linearUpwind grad(U);
          div(div(phi,U)) Gauss linear;
          div(phi,k)      Gauss limitedLinear 1;
          div(phi,omega)  Gauss limitedLinear 1;
          div((nuEff*dev2(T(grad(U))))) Gauss linear;
      }
      
      laplacianSchemes
      {
          default         Gauss linear corrected;
      }
      
      interpolationSchemes
      {
          default         linear;
      }
      
      snGradSchemes
      {
          default         corrected;
      }
      
      wallDist
      {
          method meshWave;
      }
      
      FoamFile
      {
          version     2.0;
          format      ascii;
          class       dictionary;
          object      fvSolution;
      }
      // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
      solvers
      {
          "pcorr.*"
          {
              solver           GAMG;
              tolerance        0.02;
              relTol           0;
              smoother         GaussSeidel;
          }
      
          "(p|Phi)"
          {
              $pcorr;
              tolerance        1e-7;
              relTol           0.01;
          }
      
          pFinal
          {
              $p;
              tolerance        1e-7;
              relTol           0;
          }
      
          "(U|k|omega)"
          {
              solver          smoothSolver;
              smoother        symGaussSeidel;
              tolerance       1e-06;
              relTol          0.1;
          }
      
          "(U|k|omega)Final"
          {
              $U;
              tolerance       1e-06;
              relTol          0;
          }
      
          cellDisplacement
          {
              solver          GAMG;
              tolerance       1e-5;
              relTol          0;
              smoother        GaussSeidel;
          }
      }
      
      PIMPLE
      {
          correctPhi          yes;
          nOuterCorrectors    2;
          nCorrectors         2;
          nNonOrthogonalCorrectors 1;
      }
      
      relaxationFactors
      {
          fields
          {
              p                  0.3;
          }
          equations
          {
              "(U|k|omega)"      0.7;
              "(U|k|omega)Final" 1.0;
          }
      }
      
      potentialFlow
      {
          nNonOrthogonalCorrectors 15;
      }
      
      cache
      {
          grad(U);
      }
      
      1. dynamicMeshDict
      FoamFile
      {
          version         2;
          format          ascii;
          class           dictionary;
          object          dynamicMeshDict;
      }
      
      dynamicFvMesh   dynamicMotionSolverFvMesh;
      
      motionSolverLibs ( "libsixDoFRigidBodyMotion.so" );
      
      solver          sixDoFRigidBodyMotion;
      
      sixDoFRigidBodyMotionCoeffs
      {
          patches         ( CYLINDER );
          innerDistance   2.5;
          outerDistance   20;
          mass            11.4586;
          centreOfMass    ( 0.280485 0 0 );
          momentOfInertia ( 0.6414 0.6414 0.24 );
          g               ( 0 0 0 );
          rho             rhoInf;
          rhoInf          1.225;
          report          on;
          solver
          {
              type            Newmark;
              gamma           0.5;
              beta            0.25;
          }
          constraints
          {
              yLine
              {
                  sixDoFRigidBodyMotionConstraint line;
                  centreOfRotation ( 0.280485 0 0 );
                  direction       ( 0 1 0 );
              }
              noRotation
              {
                  sixDoFRigidBodyMotionConstraint orientation;
              }
          }
          restraints
          {
              verticalSpring
              {
                  sixDoFRigidBodyMotionRestraint linearSpring;
                  anchor          ( 0.280485 0 0 );
                  refAttachmentPt ( 0.280485 0 0 );
                  stiffness       10.70694631;
                  damping         0;
                  restLength      0;
              }
          }
      }
      

      请各位老师帮忙看看是什么问题,困扰好久了。。。谢谢!

      W C 2 条回复 最后回复 回复 引用
      • W
        WYing @WYing 最后由 编辑

        @wying 没有人看嘛。。。救救孩子。。。

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

          @wying 动网格这面最大的问题就是网格变形过大导致的库朗数特别大,需要大幅度降低时间步长保证收敛。你需要保证时间步长充分的小。

          debug的时候要检查边界条件是否有问题可以先忽略网格变形,先直接算看能不能算出大体的结果。如果正确的话说明不是边界条件的问题。

          线上CFD课程开始报名:http://www.dyfluid.com/class.html

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

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

            @李东岳 感谢李老师的回复!我按照您的建议,先计算了钝体绕流(不加动网格)的情况,发现能够得到合理的流场和大致结果,这应该说明了边界条件设置没有问题。那么就是动网格的设置出错了?我计算的是圆柱+不同长度分流板的流致振动,奇怪的是,大部分算例可以正常计算,只有某些长度的分流板会导致结果发散。所有算例的网格划分方式和参数设置都一样的,真是百思不得其解。。。

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

              @wying 可能某些长度的板子网格变形太大了

              线上CFD课程开始报名:http://www.dyfluid.com/class.html

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

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

                @李东岳 我也想过这个可能性。如果是因为网格变形增大导致的发散,那么至少在开始一段时间内变形较小的时候计算能正常进行。但是这些发散的case从一开始的流场就有问题,而不是说随着网格变形增大,逐渐发散。。

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

                  @wying 如果发散的算例网格比较少、也不涉密、可以发上来我给你看看

                  线上CFD课程开始报名:http://www.dyfluid.com/class.html

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

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

                    @李东岳 感谢您的回复! crash.zip 这是我的case文件,由于大小限制,文件夹内没有包含最后一个时间步,但是大概在time=1.8所有停止计算。欢迎您提出的任何意见或建议!

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

                      @wying 不好意思明天我就要去杭州了... 国庆后回来 只能看看有没有其他大佬能给你看看了,:135:

                      线上CFD课程开始报名:http://www.dyfluid.com/class.html

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

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

                        @李东岳 谢谢李老师的回复!我尝试使用高版本的OF计算,发现相同设置下算例可以正常计算,不确定这是不是低版本的OF涉及动网格时存在bug。不过有个case在所有版本OF下都发散。。。能麻烦李老师有空的话帮忙看一下吗?算例文件case.zip (大约在time=1.43发散)十分感谢!!

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

                          好的 我给你瞅瞅最新的这个case.zip

                          线上CFD课程开始报名:http://www.dyfluid.com/class.html

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

                          1 条回复 最后回复 回复 引用
                          • C
                            cresendo @WYing 最后由 编辑

                            @wying 试试看对p不要进行松弛

                            W 1 条回复 最后回复 回复 引用
                            • W
                              WYing @cresendo 最后由 编辑

                              @cresendo 谢谢您的建议!去掉p松弛之后,算例确实可以计算下去,但是CFL保持在20左右,且速度场和压力场也明显有问题。所以还是不太对,不过还是感谢您的帮助!

                              C 1 条回复 最后回复 回复 引用
                              • C
                                cresendo @WYing 最后由 编辑

                                @wying 那应该就是动网格的问题了,试试在dynamicMeshDict里将松弛因子调成0.4试试

                                W 1 条回复 最后回复 回复 引用
                                • W
                                  WYing @cresendo 最后由 编辑

                                  @cresendo 请问您说的dynamicMeshDict里的松弛因子是 accelerationRelaxation 0.9 和 accelerationDamping 0.95 吗?

                                  C 1 条回复 最后回复 回复 引用
                                  • C
                                    cresendo @WYing 最后由 编辑

                                    @wying 前者

                                    W 1 条回复 最后回复 回复 引用
                                    • W
                                      WYing @cresendo 最后由 编辑

                                      @cresendo 我按照您的建议,调整了accelerationRelaxation 0.4,计算可以正常进行,CFL数小于1, 流场看上去也正常。

                                      但有一个小疑问,通过查询dynamicMeshDict文件的官方描述,得知accelerationRelaxation直接减小刚体的加速度,比如在某时间步solver求出的加速度为10m/s^2,如果设置accelerationRelaxation=0.4,那么实际应用的加速为4m/s^2,这是一个很大的松弛。而官方描述中建议的accelerationRelaxation范围是0.9-1,并且提到 “Be careful with this accelerationRelaxation. Too low of a value will mean that the Body does not respond to the fluid forces correctly”。因此我担心设置较小的加速度松弛会不会对结果产生太大影响?因为没有实验结果作为参照,所以无法判断结果的准确性。不知道您一般是如何处理这种误差的?谢谢!

                                      C 1 条回复 最后回复 回复 引用
                                      • C
                                        cresendo @WYing 最后由 编辑

                                        @wying 你说的这种松弛方式是不正确的,应该是上一时间步的0.6加上当前时间步的0.4,这样松弛以后获得最终的加速度。至于wikidict上说的这个建议,松弛应该影响是的最终稳定的速度,但是可能的确会导致响应运动不是那么准确。因此,还是建议你寻找一个标模实验验证一下(这种实验很多):146:

                                        W 1 条回复 最后回复 回复 引用
                                        • W
                                          WYing @cresendo 最后由 编辑

                                          @cresendo 我重新查看了source code,aRel的松弛方式确实如您所说,感谢指出错误!我对比了不同aRel的结果,发现Rel = 0.4-0.8都可以正常运行且结果几乎相同,但是aRel >0.8之后会发散。我不确定为什么wikidict上说太小的aRel会造成不准确的响应,至少我的结果并不受aRel取值的影响。相反,恰当的aRel可以使本来发散的case收敛。我之前没有指定过aRel,也就是取默认值1.0,因此比较好奇为什么您会设定aRel=0.4而不是其他值,是根据模拟的经验吗?谢谢!

                                          C 1 条回复 最后回复 回复 引用
                                          • C
                                            cresendo @WYing 最后由 编辑

                                            @wying 0.4是比较保守的设置,我也是根据自带的DTCHull算例参考的。但是对于涡激运动问题,以我个人2D的尝试,0.9是可以跑的。我本人是做涡激运动的,这是我的songtao.chen@sjtu.edu.cn,有问题可以一起交流~

                                            W 1 条回复 最后回复 回复 引用
                                            • W
                                              WYing @cresendo 最后由 编辑

                                              @cresendo 在 计算钝体涡激振动时发散(pimpleDyMFoam) 中说:

                                              songtao.chen@sjtu.edu.cn

                                              谢谢!非常感谢您的建议!

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