本文简单描述OpenFOAM中线性方程组的求解。

注：本文内容整理自Wolf Dynamics公司的公开培训教材“AOpenFOAM® Introductory Training Online session – 2020 Edition”。建议英文过得去的道友自行查看英文原版。

”

在计算区域的每一个控制体内进行空间离散和时间离散后，控制方程可以转化为下面的形式：

对计算域内的所有网格单元，离散方程可以组装成矩阵方程：

此矩阵方程可使用迭代法或直接法进行求解。

1 fvSolution字典

fvSolution字典文件中包含了线性方程组求解过程中的一些控制参数（如求解方法、残差等）。如下面的fvSolution字典示例：

solvers
{
    p
    {
        // 指定求解器为PCG
        solver           PCG;
        // 选择预处理器为DIC
        preconditioner   DIC;
        // 指定计算残差为1e-6
        tolerance        1e-06;
        // 指定相对残差为0，表示禁用
        relTol           0;
    }
    // pFinal指的是最终压力校正
    pFinal
    {
        // 采用与p方程相同的求解参数
        $p;
        relTol    0;
    }
 
    U
    {
        // 采用PBiCGStab求解U方程
        solver           PBiCGStab;
        // 预处理器选用DILU
        preconditioner      DILU;
        // 残差设置为1e-8
        tolerance        1e-08;
        relTol           0;
        // 设置最小迭代次数
        minIter            3;
        // 设置最大迭代次数
        maxIter            1000;
    }
} 
 
PISO
{
    nCorrectors    2;
    nNonOrthogonalCorrectors    1;
}

在上面的示例文件中：

• 求解器子词典中指定了用于每个要求解的方程所采用的线性求解器。如示例中p方程使用PCG求解器，U方程使用PBiCGStab求解器
• 线性求解器区分对称矩阵和非对称矩阵
• 如果忘记定义线性方程求解器或使用了错误的求解器，OpenFOAM在求解计算时会以错误提示的方式进行通知
• 采用迭代法进行线性方程组求解，迭代残差是一个非常重要的衡量计算结果的指标，残差越小意味着解的精度越高
• 在求解特定物理场方程之前，会根据该物理场的当前值评估初始残差。每次求解迭代后，将重新计算残差。
• 如果达到以下任一条件，求解器将停止：残差低于所设定的求解器残差、当前残差与初始残差之比低于求解器相对容差relTol、迭代次数超过最大迭代次数maxIter。
• 关键字maxIter是可选的，默认值为1000，还可以使用关键字miniter定义最小迭代次数，默认值为0。

2 线性求解器

以下是OpenFOAM提供的线性求解器：

• GAMG：多重网格求解器
• PCG：Newton-Krylov求解器
• PBiCG：Newton-Krylov求解器
• smoothSolver：光顺求解器
• PBiCGStab：Newton-Krylov求解器
• diagonalSolver

求解器源代码位于文件夹$WM_PROJECT_DIR/src/OpenFOAM/matrices/lduMatrix/solvers中，当使用Newton-Krylov求解器时，需要定义预条件。

OpenFOAM中的预条件包括：DIC、FDIC、diagonal、DILU、GAMG、noPreconditioner。预条件的源代码路径为$WM_PROJECT_DIR/src/OpenFOAM/matrices/lduMatrix/preconditioners。

smoothSolver求解器需要指定光顺器，OpenFOAM提供的光顺器包括：DIC、DICGaussSeidel、DILU、DILUGaussSeidel、FDIC、GaussSeidel、nonBlockingGaussSeidel、symGaussSeidel。光顺器的源代码位于$WM_PROJECT_DIR/src/OpenFOAM/matrices/lduMatrix/smoothers。

3 线性求解器的残差

关于OpenFOAM中的线性方程求解器：

• 线性求解器的选择没有理论建议，其与具体的问题及计算机硬件有关(网格类型、涉及的物理模型、处理器缓存、网络连接、分区方法等)。
• 大多数情况下，可以使用GAMG方法(geometric algebraic multigrid，几何代数多重网格方法)，该方法对于对称矩阵(例如压力)来说是最佳选择。
• GAMG方法收敛速度快(通常迭代次数小于20次)。如果想要迭代更多次数，可尝试更改光顺器。
• 如果计算时间太长或计算不稳定，可以使用PCG求解器。
• 当运行多核(1000个以上)时，使用PCG可能是更好的选择。
• 对于非对称矩阵，带DILU预条件的PBiCGStab方法是一个很好的选择。
• 带有GaussSeider光顺器的smoothSolver求解器的性能也非常好。
• 如果带有DILU预处理器的PBiCGStab方法神秘地崩溃，并出现与预条件器相关的错误，可以使用smoothSolver解算器或更换预条件器，但通常情况下PBiCGStab求解器比smoothSolver求解器更快。
• 非对称矩阵由速度场(U)和输运标量场(k、ω、ε、T等)组合而成。通常情况下，速度场和输运标量场的计算开销小且计算很快，因此可以对这些物理场使用更严格的残差(如1e-8)。
• diagonal求解器用于反向替换，例如，在使用状态方程计算密度时(我们知道p和T)。

对于线性方程求解器的残差：

• 残差并不能直接表明正在向正确的求解方向收敛。
• 第一个时间步内可能计算不收敛，这通常是可以接受的。
• 可能需要在第一次迭代期间使用较小的时间步长来保证求解计算的稳定性。
• 如果求解在一段时间后不收敛，可尝试减小时间步长。

那么残差该如何设置呢？一般情况下可以参考：

• 压力方程非常重要，应该被准确地求解
• 压力方程求解是整个迭代过程中计算开销最大的部分
• 对于压力方程(对称矩阵)，可以在设置初始tolerance为1e-6、relTol等于0.01，计算一段时间后，将这些值分别更改为1e-6和0.0。
• 如果线性解算器占用的时间太长，可以将收敛残差更改为1e-4，relTol修改0.05
• 对于速度场(U)和输运标量场(非对称矩阵)，求解这些变量的成本相对较低，因此您可以一开始以严格的残差开始。

一个比较松的残差设置：

p
{
    solver            PCG;
    preconditioner DIC;
    tolerance         1e-6;
    relTol           0.01;
}
 
U
{
    solver            PBiCGStab;
    preconditioner DILU;
    tolerance         1e-8;
    relTol           0.001;
}
 

一个比较严格的残差设置示例：

p
{
    solver            PCG;
    preconditioner DIC;
    tolerance         1e-6;
    relTol           0.0;
}
 
U
{
    solver            PBiCGStab;
    preconditioner DILU;
    tolerance         1e-8;
    relTol           0.0;
}

PISO或PIMPLE方法与MomentumPredictor选项(默认情况下处于启用状态)一起使用时，可以选择设置最终压力校正步(PFinal)的残差。采用这种方式可以将所有计算工作仅放在最后一个校正步中(在本例中为pFinal)。而对于所有中间校正步(P)可以使用更宽松的收敛标准。如果按照这种方式操作，建议至少执行2个校正步骤(nCorrectors)。如下面pFinal残差定义：

pFinal
{
    solver            PCG;
    preconditioner DIC;
    tolerance         1e-6;
    relTol           0.0;
}

4 矩阵重排序

系数矩阵求解时，矩阵对角线越大，收敛速度越快，因此强烈建议在进行计算之前利用renumberMesh工具对系数矩阵进行重排，命令调用方式为：

renumberMesh -overwrite

renumberMesh可以显著提高线性求解器的速度，特别是在第一次迭代期间。重新排序背后的思路是使系数矩阵更加对角占优，从而加快迭代求解器的速度。

5 多重网格求解器

多重网格解算器(OpenFOAM®中的GAMG)的发展，以及高分辨率TVD格式和并行计算的发展，是CFD历史上最显著的成就之一。在大多数情况下都可以使用GAMG线性解算器。但是如果看到GAMG线性求解器收敛时间过长或收敛的迭代次数超过100次，则最好改用PCG线性解算器。OpenFOAM中的GAMG线性求解器在将计算扩展到500个以上的处理器时性能表现不是很好，此外对于某些多相流情况，PCG方法的性能优于GAMG方法。

对于对称矩阵（如压力场）采用GAMG线性求解器，大多数情况下可以使用下面的残差设置。

p
{
    solver            GAMG;
    tolerance         1e-6;
    relTol           0.01;
    smoother          GaussSeidel;
    nPreSweeps        0;
    nPostSweeps      2;
    cacheAgglomeration  on;
    agglomerator      faceAreaPair;
    nCellsInCoarsestLevel  100;
    mergeLevels       1;
    minIter 3;
}
 
pFinal
{
    solver            GAMG;
    tolerance         1e-6;
    relTol           0;
    smoother          GaussSeidel;
    nPreSweeps        0;
    nPostSweeps      2;
    cacheAgglomeration  on;
    agglomerator      faceAreaPair;
    nCellsInCoarsestLevel  100;
    mergeLevels       1;
    minIter 3;
}

6 稳态模拟

前面的残差设置都是针对瞬态计算的，对于稳态模拟，可以采用类似的残差控制方式，如下面的SIMPLE方法残差设置示例：

SIMPLE
{
    nNonOrthogonalCorrectors    2;
    residualControl
    {
        p    1e-4;
        U    1e-4;
    }
}

这里使用了一个子字典residualControl进行残差设置。