本文基于多篇icoFoam解析的文章，做了一些补充以供参考。
本文给出的源代码中包含了大量测试命令，便于大家边运行边学习，提高学习效率。
为不污染OpenFOAM源代码，单独将icoFoam源码提取出来并重命名为icoFoam_learn，方便大家使用。
参考网站：
总体算法参考：东岳流体不可压瞬态PISO算法 — OpenFOAM|CFD
其他细节参考：
无痛苦NS方程http://www.dyfluid.com/theory.pdf
OpenFOAM>>solver>>incompressible>>icoFoam的说明_苏军伟_新浪博客
学习顺序建议：
1、参考李东岳翻译的东岳流体运行方腔驱动流
2、参考 C++ 教程 | 菜鸟教程 学习C++语法，重点看明白继承，多态，模板（类模板比较抽象看可以上B站看视频课）
3、 看明白基本算法 东岳流体不可压瞬态PISO算法 — OpenFOAM|CFD
4、可以下载我注释好了的 icoFoam源代码，边运行边学，下载链接：
https://download.csdn.net/download/BeanChai/89351775?spm=1001.2014.3001.5503
/*---------------------------------------------------------------------------*\
========= |
\\ / 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"//创建对象runTime，成员函数runTime.timeName()返回当前时间
#include "createMesh.H"//读取网格
pisoControl piso(mesh);
#include "createFields.H"//创建场对象U，p，phi=flux(U)
#include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//新建测试用的中间标量场
Info<< "setting field test_out_scar\n" << endl;
volScalarField test_out_scar
(
IOobject
(
"test_out_scar",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("test_out_scar", dimensionSet(0,2,-2,0,0,0,0), 0.0),
"zeroGradient"
);
//新建测试用的中间向量场
Info<< "setting field test_out_vect\n" << endl;
volVectorField test_out_vect
(
IOobject
(
"test_out_vect",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
//dimensionedVector("test_out_vect", dimensionSet(0,1,-1,0,0,0,0), vector::zero)
//"zeroGradient"
);
//将场数据输出直接用其内置的write函数，会将此时的场输出到该时刻的文件夹中
//test_out_vect.write();
//return 0;
Info<< "\nStarting time loop\n" << endl;//显示进入时间循环
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;//显示当前运行时间
#include "CourantNo.H"//计算库朗数，并输出到终端，最大值计算公式为0.5*max(sumPhi/mesh.v*dt),sumPhi为网格通量（速度*面积）绝对值的加和
// Momentum predictor
//test_out_vect = U;//测试用 测试动量方程推进后是否更新速度场U
//这里会构造线性方程组Ax=b，存储A和b
fvVectorMatrix UEqn
(
//这里每一项返回一个系数矩阵A或右端向量b(LHS AND RHS)然后累加起来
//如果用显式fvc离散，则会计算右端项b
fvm::ddt(U)//隐式离散生成矩阵类，该项会被加到UEqu的A上
+ fvm::div(phi, U)//隐式离散生成矩阵类，该项会被加到UEqu的A上
- fvm::laplacian(nu, U)//隐式离散生成矩阵类，该项会被加到UEqu的A上
);
//从if可看出momentum prediction这一步并不是一定要求解的，那么如何控制下面momentumPredictor是否求解呢？
//我们在fvSolution里面可看到这一项。那里面momentumPredictor选择yes即求解。如果求解，会得到一个新的U，
//语句下面用到的任何关于U的函数都是momentumPredictor中求解出来的U。如果不求解，那下面一步不会被执行，之后的U用的是上一个步长的
//作者：流年非留年 https://www.bilibili.com/read/cv14971258/?from=search&spm_id_from=333.337.0.0 出处：bilibili
//求解动量方程，如果前面用fvm隐式命名空间中的函数，这里迭代多次求解，如果前面用fvc显式命名空间中的函数，这里迭代0次
if (piso.momentumPredictor())
{
//这里把系数矩阵组装起来，然后求解线性方程组
//这里求解完线性方程组后，计算得到的是u*，存储在U中，之前的结果存储在U.oldTime中
//Info<< "---------------------------"<< endl;//测试用
solve(UEqn == -fvc::grad(p));//依据初始条件求解预测速度U*，http://www.dyfluid.com/piso.html,在终端输出残差和迭代信息
}
// test_out_vect = test_out_vect - U;//测试用
// test_out_vect.write();//测试用
// U.write();//测试用
// return 0 ;//测试用
// --- PISO loop
//
//fvVectorMatrix.A()表示返回系数矩阵的对角线元素，返回volScalarField
//fvVectorMatrix.H()表示去除对角元素的系数矩阵乘以速度矩阵再加上右端项b，返回volVectorField
//参考李东岳无痛苦NS方程
while (piso.correct())
{
volScalarField rAU(1.0/UEqn.A());//标量场，因为系数矩阵A固定，即使piso循环多次，在当前时间步不会被更新
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));//向量场，因为.H命令会使用U速度场，每次piso循环会发生变动，可以用test_out_vect测试
//因为预测的速度需要满足边界条件，因此这里需要让rAU*UEqn.H()也满足速度边界条件，调用constrainHbyA施加边界条件到HbyA
surfaceScalarField phiHbyA//通过方程(32)组建HbyA，方程出处：http://www.dyfluid.com/piso.html
(
"phiHbyA",
fvc::flux(HbyA)//显式插值，有限体积法中.flux()是面通量，一般是面上的值*面的大小
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)//这是一个修正项，利用上一时间步通量和速度向面上的插值之差来修正本时间步的值，
//会增加额外的耗散，好处是计算鲁棒性好，//出自：CFD中文网-->OpenFOAM-->pEqn.H中phiHbyA计算修正项的问题
);
//Info<< "============================="<< endl;//测试用
adjustPhi(phiHbyA, U, p); //对压力全为Neumann边界条件下，给定一个附加条件使得压力有解，作者：流年非留年 出处：bilibili
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);//更正压力边界条件，保证通量守恒。作者：流年非留年 出处：bilibili
// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())//网格非正交时循环
//若在fvSolution字典文件中nNonOrthogonalCorrectors设置为0,就只求解控制方程一次,如果设置为n,则求解控制方程n+1次。OpenFOAM 采用的是非结构化的同位网格，
//非结构化网格可能出现非正交网格，需要进行非正交修正。而同位网格为了消除压力振荡，需要进行动量插值 //作者：流年非留年 出处：bilibili
//经验上，除非使用四面体网格，一般都不用非正交循环
{
// Pressure corrector
//Info<< "*****************************"<< endl;//测试用
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)//构建压力泊松方程，求解方程(33)求解压力 方程出处：http://www.dyfluid.com/piso.html
);
pEqn.setReference(pRefCell, pRefValue);//不可压缩流动中，因为只有压力梯度项，因此需要设置压力参考点
//test_out_scar = p;//测试用
//test_out_scar.write();//测试用
pEqn.solve();//求解压力p'的泊松方程,这里求解完成后p被更新为 p* = p_t + p‘，在终端输出残差和迭代信息
//p.write();//测试用
//return 0;//测试用
//Info<< "+++++++++++++++++++++++++++++"<< endl;//测试用
if (piso.finalNonOrthogonalIter())
{
//Info<< "&&&&&&&&&&&&&&&&&&&&&&&&&"<< endl;//测试用
phi = phiHbyA - pEqn.flux();//根据方程(32)修正面上的通量，注意面通量 phi=U_f & S.f（面速度*面的面积）方程出处：http://www.dyfluid.com/piso.html
//pEqn.flux() = fvc::interpolate(rUA)*fvc::snGrad(p)*mag(mesh.Sf()) //来源https://zhuanlan.zhihu.com/p/432907267
}
}
#include "continuityErrs.H"
U = HbyA - rAU*fvc::grad(p);//通过方程(30)依据压力p*以及HbyA*更新速度有U**。出处：http://www.dyfluid.com/piso.html 网站里是根据方程(32)修正
U.correctBoundaryConditions();//使用边界条件修正速度
}
runTime.write();//输出控制
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
//break;//测试用
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
扫一扫
3139
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
