影响CFD计算量的因素分析及在OpenFOAM中的参数调整

Posted 2021-06-21 希望先生

转载自苏军伟微博“影响CFD计算量的因素分析及在OpenFOAM中的参数调整”

影响 CFD 计算量的因素很多，大概可以归为一下几个部分
1）物理问题本身
物理问题本身的复杂程度直接关系到计算量。 一般而言，非线性模型的计算量要高于线性模型，多相流计算量大于单相流动。 如果单纯从求解方程个数及其方程类型而言， 方程个数越多计算量越大，比如提供例子中 square 需要解 3 个标量方程（标量 p 和向量 U （2d） ），而 dambreak 需要求解 6 个方程（标量 p、向量 U（2d）、标量体积分率 alpha，标量湍流强度 k 和标量湍流强度耗散率 epsilon）。因此 dambreak 的计算量要高于方块绕流 。湍流模拟而言，大涡模拟的计算量要高于雷诺时均。
2）计算网格单元数目和维度
计算网格的单元数直接关系到最终代数方程组的个数（每个单元求解一个代数方程）。计算单元的个数越多，代数方程组越难求解，计算量越大，因此在满足工程需求的情况下，应尽量减少网格数目，以减低计算量。 当网格数目相同时，计算区域的维度越大，得到的代数方程组越难求解。也就是说，网格数目相同的情况下也就是说 3d 的网格较 2d 的网格难求解。
3）计算网格的相对大小
对于显式或者半隐式算法（SIMPLE 或者 PISO，cfd 中较常采用），计算的稳定性受库朗数（Co=u*dt/dx<1）的限制，物理上认为，计算过程中流体在设定的时间步长内流动距离不能超过一个网格。因此，在剖分网格时，应该避免出现体积过小的网格。
4）代数方程求解器
代数方程求解器是关系计算量的关键因素之一，当计算的网格单元数比较多时， 代数方程求解器的选择尤其重要。 OpenFOAM 中的代数方程求解器有 3 类，多重网格求解器，共轭梯度求解器，光滑求解器。一般而言，计算速度的排列可以为共轭梯度求解器+多重网格预条件器  >  多重网格求解器  >  共轭梯度求解器+一般与预
条件器>光滑求解器。代数方程的求解难度和方程类型有关，对于速度方程和一般的标量方程（k 或者 epsilon），相对容易求解， 因为这些方程通常具有对角占优特性。而对于，泊松类型的方程（比如压力p 的方程），由于对角占优特性较差，则很难求解。故，对于网格数目较多的问题，通常对压力方程采用“共轭梯度求解器+多重网格预条件器”或者“多重网格求解器”来进行求解，而对于普通的标量方程则使用“共轭梯度求解器”。求解器的选择可以通过 case 文件夹下的 system 下的 fvSolution 来完成，典型的设置方法如下：
“共轭梯度求解器+多重网格求解器”
p //求解变量
{
    solver PCG; //共轭梯度求解器（使用对称矩阵，否则 PBICG）
    preconditioner
    {
        preconditioner GAMG; //代数网格预条件器
        tolerance 1e-05; //预条件器绝对残差
        relTol 0; //预条件器相对残差
        smoother DICGaussSeidel;
        nPreSweeps 0;
        nPostSweeps 2;
        nFinestSweeps 2;
        cacheAgglomeration false; //可以改为 true，用内存换时间
        nCellsInCoarsestLevel 10;
        agglomerator faceAreaPair;
        mergeLevels 1;
    }
    tolerance 1e-05; //绝对残差
    relTol 0; //相对残差
    maxIter 100; //最大迭代数目，一般而言，该求解器可以保证在100步以内收敛
}
“多重网格求解器”
p //求解变量
{
    solver GAMG; //代数多重网格求解器
    tolerance 1e-07; //绝对误差
    relTol 0.01; //相对误差
    smoother DIC; //光滑器
    nPreSweeps 0;
    nPostSweeps 2;
    nFinestSweeps 2;
    cacheAgglomeration true;
    nCellsInCoarsestLevel 10;
    agglomerator faceAreaPair;
    mergeLevels 1;
}
“共轭梯度求解器” （非对称矩阵）
"(U|k)Final"
{
    solver  PBiCG;
    preconditioner DILU;
    tolerance 1e-08;
    relTol 0;
}
“共轭梯度求解器”（对称矩阵）
p
{
    solver PCG;
    preconditioner DIC;
    tolerance 1e-8;
    relTol 0.01;
}
5）代数方程残差要求
第 4）部分给出了常见代数方程器的设置。应当注意到设置过程中有两个残差，他们分别为 relTol（相对残差）和 tolerance（绝对残差），方程组迭代过程中，哪个残差最先达到，则方程组停止迭代，这两个值设置越小，则计算量越大，但精度会越高。在实际计算过程中可以将其中的一个设置为 0，则达到另外一个残差，则停止迭代，当然不能两个都设置为 0（因为残差为 0 永远达不到）。在 OpenFOAM 求解器中压力方程常见有 p(或者 p_rgh)及其pFinal(p_rghFinal)，p(或者 p_rgh)  用来求解压力修正方程的中间步骤，而 pFinal(p_rghFinal)为最后一次修正过程的代数方程求解的残差。 一般而言，将中间步骤的残差设置较大以降低计算量，而将最后一次修正过程的残差设置较小，以提交求解精度。
6）代数方程类型
代数方程类型有两方面可能影响程序的计算量：（1）微分方程离散后得到的代数方程组的对角占优特性，对角占优的方程组更容易收敛；（2） 微分方程离散后得到的代数方程组的是否主对角对称 （目前代数方程求解器大都针对对称矩阵方程组和反对称方程组，不对称系数矩阵方程组可以转化为对称矩阵和反对称矩阵之和。因此，一般而言，非对称矩阵的计算量是对称矩阵的两倍） 。常见的求解变量的对上述两个条件满足状况如下压力方程 p（不可压缩非稳态流动） ：对称矩阵，但对角相等，非对角占优，因此计算量较大。速度方程 U（不可压缩非稳态流动） ：非对称矩阵，对角占优，计算量通常较 p 小。一般标量方程 k，epsilon，T 等（非稳态流动）：非对称矩阵，对角占优，计算量较小。非稳态纯扩散方程（比如纯导热问题）：对称矩阵，对角占优，计算量最小；因此，在计算过程中，可以对压力方程 p 选择“共轭梯度求解器+代数多重网格预条件器”而速度 U 和其他标量方程选择一般的共轭梯度求解器。选择的方法见第 4）部分。
7）离散格式
对流项的离散格式在一定程度上会影响计算量，比如一阶迎风（upwind） ，数值粘性比较大，计算稳定，可以适当增加时间步长。而对于高阶格式如 Gamma 格式，数值粘性相对较较小，计算稳定性相对较差，时间步长可能要小一点，从而影响计算量。同时，高阶格式求解插值系数计算量也较大。离散格式的调整可以通过 system 下的 fvScheme 文件下的 divScheme 字典下进行填写。当不知道选择什么的时候，可以随便填写一个字符串，然后根据屏幕提示进行填写。比如：
div(phi, U) Gauss aaa;
时间项的离散也会影响到计算量，比如欧拉格式（Euler）计算量较 backward 和 CN 格式小。但精度也想应降低。 时间格式在 system 下的 fvScheme 文件中 ddtScheme 字典下进行填写。同样，当不知道填写什么的时候，可以随便填写一个字符串， 然后根据提示进行填写。比如：
ddt(U) aaa;
8）时间步长
时间步长越小在给定的时间内计算所需要的计算时间越长，精度也就越高。 但是，一般而言，时间步长的选择通常根据当前计算的库朗数进行选择，也就是dt <(dx/U)。OpenFOAM 中通过在 system 下的 controlDict 中的 deltaT（以秒为单位）来设置。OpenFOAM 中有的求解器是可以自动调节时间步长的，比如给出的算例中 dambreak 和hotRoom，如果要自动调节时间步长需要将 controlDict 文件中的 adjustTimeStep  后面设置为on，并通过 maxCo（最大库朗数）来自动条件时间步长。maxCo 应该小于 1，为了改善计算精度可以将 maxCo 调节的小一些，但会增加计算量。   应当指出， Co 和 deltaT 呈线性关系，如果你设置了 maxCo  为 0.5 可以比设置为 0.25 降低一倍的计算量。一般将 maxCo 设置为0.5 即可。对于 damBreak 这个例子还有一个 maxAlphaCo，他的意义和 maxCo 类似，决定了求解体积分率方程 alpha1 的步进步长，可以将 maxAlphaCo 设置和 maxCo 一样的值。
9）最终结束时间
最终结束时间不会影响单步计算时间，但会影响总体运行时间。可以通过 system 下的controlDict 文件中的 endTime 来调节。
10）并行分块方法
并行过程中分块方法可能影响到计算量，一般而言，应该让各个处理器之间的交接面最小，以减少通信量，具体方法和参数请查看附件中 decomposeParDict。