关于Open Foam9的RheoTool流变工具包模拟毫米量级的几何模型不收敛问题
-
@李东岳 在 关于Open Foam9的RheoTool流变工具包模拟毫米量级的几何模型不收敛问题 中说：
你看一下你的求解器源代码是瞬态算法还是稳态算法。从你的fvSolution以及fvScheme来看，两者混到一起了。
我查阅了一下，在rheoTool算例文件
中入口速度是和时间有关的
这应该是一个瞬态的算例，但是在fvSolution中，也有SImpleSIMPLE { nInIter 1; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; residualControl { } } relaxationFactors { fields { p 1; } equations { U 1; tau 1; theta 1; } }此外我将Simple改为PISO，运行的时候会报错，提示fvsolution文件中无SIMPLE。我查了一下用户文档，发现关于fvSolution有这样一段描述
这又表明当SIMPLE存在时，确实是稳态的情况，尽管rheoFoam是瞬态求解器。现在我搞不清楚到底该怎么设置。从所有的算例文件夹下，无论是瞬态还是稳态都是有SIMPLE这个字段的。
另外，在用户指导最后一页看到这样一句话
说是ηp>ηs是，尽可能将耦合设置为 none，按照他说的设置了，仍然不起作用。
总之，试了很多方法仍然不收敛，期待大佬帮忙解决 -
@tiankai @李东岳
补充一下：和李老师讲的一样，初步判断是对流项的数值格式不稳定导致的。打开FENE-P的算例文件夹中，模拟了一个宽度为1m的槽道流，入口速度指定为1m/s。并且在该算例中，忽略掉了对流项
divSchemes { default none; div(tau) Gauss linear; div(grad(U)) Gauss linear; div(phi,U) GaussDefCmpw none; div(phi,theta) GaussDefCmpw cubista; div(phi,1|f) GaussDefCmpw cubista; div(phi,A) GaussDefCmpw cubista; div(phi,tau) GaussDefCmpw cubista; div(phi,C) GaussDefCmpw cubista; }执行之后，发现是收敛的。
随后我将入口速度改为了100m/s，1000m/s，执行之后不收敛。
然后修改对流项格式 div(phi,U) GaussDefCmpw cubista; 运行之后收敛。不加对流项不收敛，添加上对流项反而收敛，有点意思？最后，我将槽道的宽度设置为1mm，入口速度设为1000mm/s（也就是将通道尺度缩小1000倍，入口速度缩小1000倍），同样不忽略对流项，运行之后，发现不收敛。在对流项几种格式都进行了尝试后，还是不收敛。难道是对流项的数值格式有bug？
希望同感兴趣的同行一块讨论：tiankai0721@163.com -
@李东岳 李老师，这是对上面的内容重新梳理，谢谢您啦
我使用rheoTool工具包，模拟了几何长度为36mm，宽度为4mm，在通道的两侧有4mm×4mm方腔的二维流动；模型为FENE-P非牛顿模型。我设置入口流量为3.62e-7m/s（或入口速度0.008m/s），使用了自适应时间步长，最大Co数设置为1。运行后发现不收敛。求解器为rheoFoam瞬态求解器。
针对不收敛的情况进行了尝试
1、在fvschemes文件中，时间使用了Euler，fvsolution中算法是SIMPLE。这里需要说明的是在rheoTool工具包中，SIMPLE是瞬态的。因为在rheoTool中，所有的算例的fvsolution文件中，均有SIMPLE字段。因为在算例文件中，大量的入口条件是和时间相关的，这种流动问题显然不可能是稳态流动。（我在将SIMPLE修改为PISO报错了，提示不存在SIMPLE字段）下面是用户手册的一段关于fvsolution的说明
手册中说，SIMPLE确实是稳态的算法，但是rheoFoam默认是瞬态求解器，只要在ddtSchemes中设置为空便是稳态求解器。另外手册中是推荐使用瞬态算法的，因为这样会使矩阵尽可能的对角占优，数值格式更稳定一点。
2、另外我将松弛因子改小也不能保证收敛。
3、在用户指导最后一页看到这样一句话，说是ηp>ηs是，尽可能将耦合设置为 none，按照他说的设置了，仍然不起作用。
4、当我将模型改为Oldroyd-B时，发现算例可以收敛，所以肯定是数值格式的问题。
5、所以我又运行了FENE-P的算例文件：模拟了一个宽度为1m的槽道流，入口速度指定为1m/s。运行之后发现可以收敛。于是我将入口速度改为100m/s，1000m/s，运行之后就发现不收敛了。通过查看fvSchemes发现，对流项的导数设置为空，于是我修改对流项的格式为div(phi,U) GaussDefCmpw cubista; ，运行之后便收敛了。（这里忽略对流项发散，添加对流项收敛，有点意思？）
最后，我将槽道的宽度设置为1mm，入口速度设为1000mm/s（也就是将通道尺度缩小1000倍，入口速度缩小1000倍），同样不忽略对流项，运行之后，发现不收敛。下图是对流项几种可供下选择的数值格式。在对流项几种格式都进行了尝试后，还是不收敛。难道是对流项的数值格式有bug？
总结：通过一系列的尝试，终于找到了原因。发现就是FENE-P的数值格式不稳定导致的，大概率是对流项的原因。也就是说FENE-P在毫米尺度上对流项会导致数值不稳定，从而发散。但是我看到一篇关于FENE-P模型的文章（10.1103/PhysRevFluids.6.033304），论文中使用的也是mm量级的几何模型，却可以达到收敛。请求大佬帮忙解决，或者感兴趣的同行一块交流：tiankai0721@163.com.
下面是我的各个文件的设置情况
0文件夹设置如下
\*---------------------------------------------------------------------------*/ FoamFile { format ascii; class volVectorField; location "0"; object U; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 1 -1 0 0 0 0]; internalField uniform (0 0 0); boundaryField { INLET { type flowRateInletVelocity; volumetricFlowRate 3.36e-7; value uniform (0 0 0); //猜测初始速度 } OUTLET { type zeroGradient; } LOWERWALL { type fixedValue; value uniform (0 0 0); } UPPERWALL { type fixedValue; value uniform (0 0 0); } frontAndBack { type empty; } } // ************************************************************************* // \*---------------------------------------------------------------------------*/ FoamFile { format ascii; class volSymmTensorField; object tau; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [1 -1 -2 0 0 0 0]; internalField uniform (0 0 0 0 0 0); boundaryField { INLET { type fixedValue; value uniform (0 0 0 0 0 0); } OUTLET { type zeroGradient; } LOWERWALL { type linearExtrapolation; value uniform (0 0 0 0 0 0); } UPPERWALL { type linearExtrapolation; value uniform (0 0 0 0 0 0); } frontAndBack { type empty; } } // ************************************************************************* // FoamFile { format ascii; class volScalarField; object p; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dimensions [0 2 -2 0 0 0 0]; internalField uniform 0; boundaryField { INLET { type zeroGradient; } OUTLET { type fixedValue; value uniform 0; } LOWERWALL { type zeroGradient; } UPPERWALL { type zeroGradient; } frontAndBack { type empty; } } 物性参数设置如下 FoamFile { format ascii; class dictionary; object constitutiveProperties; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // parameters { type FENE-P; rho rho [1 -3 0 0 0 0 0] 1.; etaS etaS [1 -1 -1 0 0 0 0] .01; etaP etaP [1 -1 -1 0 0 0 0] .99; lambda lambda [0 0 1 0 0 0 0] 0.5; L2 L2 [0 0 0 0 0 0 0] 100; solveInTau false; stabilization coupling; } passiveScalarProperties { solvePassiveScalar no; D D [ 0 2 -1 0 0 0 0 ] 1e-9; } // ************************************************************************* // fvSchemes文件如下 \*---------------------------------------------------------------------------*/ FoamFile { format ascii; class dictionary; object fvSchemes; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ddtSchemes { default Euler; } gradSchemes { default Gauss linear; grad(p) Gauss linear; grad(U) Gauss linear; linExtrapGrad Gauss linear; } divSchemes { default none; div(tau) Gauss linear; div(grad(U)) Gauss linear; div(phi,U) GaussDefCmpw cubista; div(phi,1|f) GaussDefCmpw cubista; //FENE-P本构方程中的f div(phi,A) GaussDefCmpw cubista; //A div(phi,theta) GaussDefCmpw cubista; div(phi,tau) GaussDefCmpw cubista; div(phi,C) GaussDefCmpw cubista; // } laplacianSchemes { default none; laplacian(eta,U) Gauss linear orthogonal; laplacian(p|(ap-H1)) Gauss linear orthogonal; } interpolationSchemes { default linear; } snGradSchemes { default orthogonal; } fluxRequired { default no; p; } // ************************************************************************* // fvSolution 如下 \*---------------------------------------------------------------------------*/ FoamFile { format ascii; class dictionary; object fvSolution; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // solvers { p { solver PCG; preconditioner DIC; tolerance 1e-10; relTol 0; minIter 0; maxIter 800; } U { solver PBiCG; preconditioner { preconditioner DILU; } tolerance 1e-8; relTol 0; minIter 1; maxIter 1000; } "(theta|tau|C|A)" { solver PBiCG; preconditioner { preconditioner DILU; } tolerance 1e-12; relTol 0; minIter 0; maxIter 1000; } } SIMPLE { nInIter 1; nNonOrthogonalCorrectors 0; pRefCell 0; pRefValue 0; residualControl { } } relaxationFactors { fields { p 0.9; } equations { U 0.5; tau 0.5; theta 0.5; } } // ************************************************************************* // ControlDict文件如下 FoamFile { format ascii; class dictionary; object controlDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // application rheoFoam; startFrom latestTime; startTime 0; stopAt endTime; endTime 1; deltaT 1e-4; writeControl runTime; writeInterval 0.0001; purgeWrite 0; writeFormat ascii; writePrecision 12; writeCompression compressed; timeFormat general; timePrecision 10; graphFormat raw; runTimeModifiable yes; adjustTimeStep on; maxCo 0.8; maxDeltaT 0.001;
