计算钝体涡激振动时发散(pimpleDyMFoam)
-
各位老师好,我使用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的时间延长了,但速度场和压力场还是一开始就不正确。两种网格质量都没问题。
我附加了我的算例设置文件,和流场图以供参考。
-
粗网格及流场
-
细网格及流场
-
边界条件(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; } }
- 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); }
- 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; } } }
请各位老师帮忙看看是什么问题,困扰好久了。。。谢谢!
-
-
@wying 没有人看嘛。。。救救孩子。。。
-
@wying 动网格这面最大的问题就是网格变形过大导致的库朗数特别大,需要大幅度降低时间步长保证收敛。你需要保证时间步长充分的小。
debug的时候要检查边界条件是否有问题可以先忽略网格变形,先直接算看能不能算出大体的结果。如果正确的话说明不是边界条件的问题。
-
@李东岳 感谢李老师的回复!我按照您的建议,先计算了钝体绕流(不加动网格)的情况,发现能够得到合理的流场和大致结果,这应该说明了边界条件设置没有问题。那么就是动网格的设置出错了?我计算的是圆柱+不同长度分流板的流致振动,奇怪的是,大部分算例可以正常计算,只有某些长度的分流板会导致结果发散。所有算例的网格划分方式和参数设置都一样的,真是百思不得其解。。。
-
@wying 可能某些长度的板子网格变形太大了
-
@李东岳 我也想过这个可能性。如果是因为网格变形增大导致的发散,那么至少在开始一段时间内变形较小的时候计算能正常进行。但是这些发散的case从一开始的流场就有问题,而不是说随着网格变形增大,逐渐发散。。
-
@wying 如果发散的算例网格比较少、也不涉密、可以发上来我给你看看
-
-
@wying 不好意思明天我就要去杭州了... 国庆后回来 只能看看有没有其他大佬能给你看看了,
-
-
好的 我给你瞅瞅最新的这个case.zip
-
@wying 试试看对p不要进行松弛
-
@cresendo 谢谢您的建议!去掉p松弛之后,算例确实可以计算下去,但是CFL保持在20左右,且速度场和压力场也明显有问题。所以还是不太对,不过还是感谢您的帮助!
-
@wying 那应该就是动网格的问题了,试试在dynamicMeshDict里将松弛因子调成0.4试试
-
@cresendo 请问您说的dynamicMeshDict里的松弛因子是 accelerationRelaxation 0.9 和 accelerationDamping 0.95 吗?
-
@wying 前者
-
@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”。因此我担心设置较小的加速度松弛会不会对结果产生太大影响?因为没有实验结果作为参照,所以无法判断结果的准确性。不知道您一般是如何处理这种误差的?谢谢!
-
@wying 你说的这种松弛方式是不正确的,应该是上一时间步的0.6加上当前时间步的0.4,这样松弛以后获得最终的加速度。至于wikidict上说的这个建议,松弛应该影响是的最终稳定的速度,但是可能的确会导致响应运动不是那么准确。因此,还是建议你寻找一个标模实验验证一下(这种实验很多)
-
@cresendo 我重新查看了source code,aRel的松弛方式确实如您所说,感谢指出错误!我对比了不同aRel的结果,发现Rel = 0.4-0.8都可以正常运行且结果几乎相同,但是aRel >0.8之后会发散。我不确定为什么wikidict上说太小的aRel会造成不准确的响应,至少我的结果并不受aRel取值的影响。相反,恰当的aRel可以使本来发散的case收敛。我之前没有指定过aRel,也就是取默认值1.0,因此比较好奇为什么您会设定aRel=0.4而不是其他值,是根据模拟的经验吗?谢谢!
-
@wying 0.4是比较保守的设置,我也是根据自带的DTCHull算例参考的。但是对于涡激运动问题,以我个人2D的尝试,0.9是可以跑的。我本人是做涡激运动的,这是我的songtao.chen@sjtu.edu.cn,有问题可以一起交流~
-