在分析这一部分代码之前,需要先了解该求解器的控制方程以及其基本离散格式,具体可参考东岳老师的网站:
在理解了icoFoam求解器所要解决的问题以及相应方法后,需要对主程序代码做进一步分析:本文参考了[openfoamwiki.net]和CloudBird07的文章做了进一步整理。
// 有限体积库头文件的集合
#include "fvCFD.H"
//定义PISO循环,使用PISO循环必备头文件
#include "pisoControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//进入主程序
int main(int argc, char *argv[])
{
#include "setRootCase.H" //根据输入参数argc和argv设置算例根目录rootcase,必备头文件
#include "createTime.H" //创建时间对象,由runTime控制,非定常求解器必备头文件
#include "createMesh.H" //由 $case/constant/polyMesh,创建网格对象
//从网格mesh对象构造类pisoControl(基类:solutionControl)的对象piso
pisoControl piso(mesh);
#include "createFields.H" //创建场对象
#include "initContinuityErrs.H" //初始化连续误差
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl; //显示进入时间循环
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;//显示当前运行时间
#include "CourantNo.H" //计算库朗数并输出至屏幕
// Momentum predictor,去除压力项,构造动量方程算子,隐式离散
fvVectorMatrix UEqn //构造体积矢量场矩阵
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);
//如果进行速度预测(on),则求解完整的动量方程得到预测速度
//如果不进行速度预测(off),则预测速度直接取当前已知时间步的速度
//momentumPredictor:动量预测求解开关,对多相流以及低雷诺数一般设置为off;
if (piso.momentumPredictor())
{
//求解动量方程,由上一时刻的压力p进行计算,因此又称速度预测
solve(UEqn == -fvc::grad(p));//solve是FOAM名称空间的全局函数,参数为矩阵(fvMatrix)
}
// --- PISO loop
while (piso.correct())
{
// rAU:在速度的最后一个解中,从矩阵中提取对角项并存储倒数
//由于对流的非线性,矩阵系数是U的函数
volScalarField rAU(1.0/UEqn.A());
//UEqn.H是UEqn减去非对角项和U的乘积
//Jocobi 迭代
//可参考constrainHbyA的源码,这部分主要是代码上的考虑
//返回结果:rAU*UEqn.H()
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
//将上述体积矢量场转化为面心标量场(参考有限体积法)
//fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)的作用是保证速度通量的全局守恒,以确保压力方程有解
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA) //(fvc::interpolate(U) & mesh.Sf())
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
//调整边界值,保证速度边界条件守恒
adjustPhi(phiHbyA, U, p);
//更正压力边界条件,保证通量守恒
constrainPressure(p, U, phiHbyA, rAU);
// 网格非正交压力修正循环
//如果在fvSolution字典文件中设置为0,就只求解控制方程一次,如果设置为n,则求解控制方程n-1次
// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector
// 构造压力方程进行求解
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);
//在不可压缩流动中,只有相对压力才起作用。除非存在边界压力
//一个网格的压力可以任意设定,以产生独特的压力解决方案
pEqn.setReference(pRefCell, pRefValue);
//通过查询system/fvSolution下的PISO更新压力参考网格并重新设定参考值
pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
//在最后一次非正交校正中,使用最新压力校正通量
if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}
//计算连续性方程误差
#include "continuityErrs.H"
//校正速度,满足边界条件(主要针对第二类边界条件)
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
//输出计算结果,写入时间文件夹
runTime.write();
//显示求解器执行时间及CPU耗时
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
//提示运行结束。
Info<< "End\n" << endl;
//返回 0
return 0;
}