icoFoam代码解析
icoFoam是OpenFOAM最基础的求解器之一,想要以OpenFOAM为基础进行二次开发的同学必须对icoFOAM的代码有所了解。关于此部分内容,由于篇幅限制,更多内容可以参考本博文对应的资源"OpenFOAM编程日志,单相不可压缩流动"或去本人的CSDN资源空间搜索下载。
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
icoFoam
Description
Transient solver for incompressible, laminar flow of Newtonian fluids.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "pisoControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
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)
);
//如果进行速度预测,则求解完整的动量方程得到预测速度
//如果不进行速度预测,则预测速度直接取当前已知时间步的速度
if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}
// --- PISO loop
while (piso.correct())
{
//UEqn.A 取出的是式12中的Ap
//下面的UEqn.H()对应解析中的H算符
//rAU的结果就是解析中的A的逆
//注意,进行H运算和A运算的时候,针对的都是最早组建的动量方程
//也就是前面所说的更新HbyA时三个系数并不更新。
//关于这部分代码的具体实现,可以参考附录7
volScalarField rAU(1.0/UEqn.A());
//可参考constrainHbyA的源码,这部分的程序主要考虑的是代码上的考虑
//就结果上来说,返回的结果就是rAU*UEqn.H()
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
//ddtPhiCorr对于在非结构网格,运用压力-速度耦合算法时,可以减少压力,速度以及速度通量的解耦
//(解耦就是用数学方法将两种运动分离开来处理问题,常用解耦方法就是忽略或简化对所研究问题影响较小的一种运动,
//只分析主要的运动的发生)。
//但新版本中使用了更具一般性的ddtCorr,替代并增强了前者的作用。
//可参见附录2
//据说fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)的作用是保证速度通量的全局守恒,以确保压力方程有解
//此部分代码的具体意义还有待商榷。
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);
//泊松方程如果给定Neumann边界条件,需要给定一个附加条件。
//在NS方程中,这个附加条件就是就是要保证速度边界条件是守恒的。
//adjustPhi即是做此修正。
adjustPhi(phiHbyA, U, p);
//修改压力边界条件确保通量守恒
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);
// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector
//构建压力方程并进行求解
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}
#include "continuityErrs.H"
//更新速度
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //