转载自苏军伟微博“影响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。