参考资料：

东岳博文：http://dyfluid.com/icoFoam.html
牛奕博文：http://wap.sciencenet.cn/blog-3410930-1175782.html?mobile=1
cloudbird07博文：https://blog.csdn.net/cloudbird07/article/details/105441128
OpenFOAMwiki： https://openfoamwiki.net/index.php/IcoFoam
OpenFOAM C++ Source Code Guide：https://cpp.openfoam.org/v8/icoFoam_8C_source.html

OpenFOAM所提供的icoFoam求解器，所求解的是非定常不可压缩层流流动问题（连续方程+不可压缩Navier-Stokes方程），其使用PISO算法来解耦速度与压力。

1 控制方程

∇ ⋅ u = 0 ∂ u ∂ t + ∇ ⋅ ( u u ) = − ∇ ( p ρ ) + ∇ ⋅ ( ν ∇ u ) \nabla \cdot {\bold{u}}=0 \\ \frac{\partial{\bold{u}}}{\partial{t}} + \nabla \cdot ({\bold{u}\bold{u}})=-\nabla\left( \frac{p}{\rho} \right) + \nabla \cdot ({\nu\nabla{\bold{u}}}) ∇⋅u=0∂t∂u+∇⋅(uu)=−∇(ρp)+∇⋅(ν∇u)

由于不可压缩问题中密度 ρ \rho ρ为常数，可以直接将其从方程中去除，即认为 ρ = 1 \rho=1 ρ=1，将 p / ρ p/\rho p/ρ计作 p p p。运动粘性 ν = μ / ρ \nu=\mu/\rho ν=μ/ρ，同样不包含密度 ρ \rho ρ。对流项 ∇ ⋅ ( u u ) \nabla \cdot ({\bold{u}\bold{u}}) ∇⋅(uu)为非线性项，通常将头一个 u \bold{u} u固化为上个时刻的 u n \bold{u}^n un。

2 算法流程

PISO算法流程如下：

a. 离散和求解动量方程来获取新速度场 u ⋆ \bold{u}^\star u⋆，注意压力梯度项已经从源项中剥离出来；

a C u u C + ∑ F ∼ N B ( C ) a F u u F = b C u − ∇ p C ( n ) a_C^\bold{u}\bold u_C+\sum_{F \sim NB(C)}a_F^\bold u \bold u_F=\bold b_C^{\bold{u}}-\nabla p_C^{(n)} aCuuC+F∼NB(C)∑aFuuF=bCu−∇pC(n)

b. 基于新速度场 u C ∗ \bold u_C^* uC∗，构建 H b y A C ⋆ \bold{HbyA}_C^\star HbyAC⋆；

H b y A C ⋆ = 1 a C u ( b C u − ∑ F ∼ N B ( C ) a F u u F ⋆ ) \bold{HbyA}_C^\star=\frac{1}{a_C^\bold u}\left(\bold b_C^{\bold u}-\sum_{F \sim NB(C)}a_F^\bold u \bold u_F^\star\right) HbyAC⋆=aCu1⎝⎛bCu−F∼NB(C)∑aFuuF⋆⎠⎞

c. 求解Possion方程得到新的压力场

∇ ⋅ ( H b y A ⋆ ) = ∇ ⋅ ( 1 a C u ∇ p ⋆ ) \nabla \cdot (\bold {HbyA}^\star) = \nabla \cdot \left( \frac{1}{a_C^\bold u} \nabla p^\star \right) ∇⋅(HbyA⋆)=∇⋅(aCu1∇p⋆)

d. 根据新的压力场 p ⋆ p^\star p⋆，更新界面通量 u f ⋆ \bold u_f^\star uf⋆和速度场 u C ⋆ \bold u_C^\star uC⋆；

u f ⋆ = H b y A f ⋆ − 1 a f u ∇ p f ⋆ u C ⋆ = H b y A C ⋆ − 1 a C u ∇ p C ⋆ \bold u_f^{\star}=\bold{HbyA}^\star_f - \frac{1}{a_f^{\bold u}} \nabla p^\star_f \\ \bold u_C^{\star}=\bold{HbyA}^\star_C - \frac{1}{a_C^{\bold u}} \nabla p^\star_C uf⋆=HbyAf⋆−afu1∇pf⋆uC⋆=HbyAC⋆−aCu1∇pC⋆

e. 重复b-d直至收敛，然后进行下一时间步的推进直至终止时间。

3 代码分析

/*---------------------------------------------------------------------------*\=========                 |\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox\\    /   O peration     | Website:  https://openfoam.org\\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation\\/     M anipulation  |-------------------------------------------------------------------------------LicenseThis file is part of OpenFOAM.OpenFOAM is free software: you can redistribute it and/or modify itunder the terms of the GNU General Public License as published bythe 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 WITHOUTANY WARRANTY; without even the implied warranty of MERCHANTABILITY orFITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public Licensefor more details.You should have received a copy of the GNU General Public Licensealong with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.ApplicationicoFoamDescriptionTransient solver for incompressible, laminar flow of Newtonian fluids.\*---------------------------------------------------------------------------*/// fvCFD.h头文件中包含了许多计算所需要的头文件，里面定义了大量类和函数// pisoControl.h头文件中则定义了pisoControl类，用于控制piso算法流程，// 有time-loop和piso-loop信息。#include "fvCFD.H"#include "pisoControl.H"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //int main(int argc, char *argv[]){// setRootCaseLists.H头文件中设置算例case的根root目录信息// createTime.h中定义了时间类Time的对象runTime// createMesh.h中定义了网格类fvMesh的对象mesh#include "setRootCaseLists.H"#include "createTime.H"#include "createMesh.H"// 定义piso算法控制类pisoControl的对象piso，用mesh构造pisoControl piso(mesh);// createFields.H中定义粘性标量nu、体心标量场压力p、体心矢量场速度U、面心标量场通量phi，// 以及指标型参考压力单元pRefCell和标量型参考压力值pRefValue等。#include "createFields.H"// initContinuityErrs.H中定义标量型累计连续误差cumulativeContErr，初始化为0。#include "initContinuityErrs.H"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //// 输出“开始时间循环了！”Info<< "\nStarting time loop\n" << endl;// 时间步逐个推进while (runTime.loop()){// 输出当前时刻Info<< "Time = " << runTime.timeName() << nl << endl;// 计算Courant数#include "CourantNo.H"// Momentum predictor// a. 动量预测，u_C^star// 离散动量方程，获得矢量方程UEqn，注意压力梯度项并未包含其中。// 这里的fvm表明均为隐式处理，即离散到线性代数方程的左端项系数矩阵中，以及右端项中。fvVectorMatrix UEqn(fvm::ddt(U)+ fvm::div(phi, U)- fvm::laplacian(nu, U));// 如果有动量预测步骤的话，求解动量方程，用所得速度来作为预测速度，// 不然的话直接用上一时刻的速度作为预测速度（即省略了步骤a，最开始不求动量方程了）// 注意，在求解动量方程的时候把压力梯度项-fvc::grad(p))显式拿了进来if (piso.momentumPredictor()){solve(UEqn == -fvc::grad(p));}// --- PISO loop// PISO循环的压力速度解耦迭代过程while (piso.correct()){// b. 构建HbyA_C^star// 体心标量场rAU，存储着动量方程对角线系数的倒数，即1/a_C^UvolScalarField rAU(1.0/UEqn.A());// 体心矢量场HbyA，UEqn.H()为动量方程离散后的右端项（不含压力梯度项）// 减去 左端非对角元系数与预测速度的乘积，rAU*UEqn.H()则再除以对角元系数// constrainHbyA(..., U, p)则再对HbyA的边界值进行一定的修正。volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));// 面形心标量场phiHbyA，为HbyA的面通量，// ddtPhiCorr项通过提取插值速度与面通量的差值，引入了面速度场的散度，// 第二项可能是做的某些修正，不是很清楚……surfaceScalarField phiHbyA("phiHbyA",fvc::flux(HbyA)+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi));// 调整进口和出口通量使其符合连续性，这对压力解的存在性，很有帮助。adjustPhi(phiHbyA, U, p);// Update the pressure BCs to ensure flux consistency// 更新压力边界条件来确保通量一致性。constrainPressure(p, U, phiHbyA, rAU);// Non-orthogonal pressure corrector loop// 非正交压力修正循环// Laplacian的非正交部分是由最近求得的压力值来计算的，使用延迟修正方法引入。// 这里多说两句：// 对于非正交网格而言，扩散项离散的时候会有交叉扩散项的引入，其一般采用当次迭代的// 变量值计算，并当做源项处理，即延迟修正引入，再进行2-3次迭代求解以使得该非正交// 修正充分引入，同时注意不要引入过多的非正交修正，其既不会收敛又耗费计算时间。// 还有就是即便是正交网格，下面的循环体部分也会执行一次，即至少求解一次该压力方程。while (piso.correctNonOrthogonal()){// Pressure corrector// c. 压力修正// 组装压力修正方程，获得标量方程pEqnfvScalarMatrix pEqn(fvm::laplacian(rAU, p) == fvc::div(phiHbyA));// 在不可压缩流动中，只有压力的相对值才是有效的，所以需要设置某个单元为// 压力参考单元，其压力值为参考压力值，来确保解的唯一性。pEqn.setReference(pRefCell, pRefValue);// 求解压力修正方程pEqn.solve();// d. 更新界面通量// 如果是最后一步非正交修正，那么修正界面通量phi值，即u_f^star修正。if (piso.finalNonOrthogonalIter()){phi = phiHbyA - pEqn.flux();}}#include "continuityErrs.H"// d. 更新速度，修正速度u_C^starU = HbyA - rAU*fvc::grad(p);// 修正速度场u的边界条件U.correctBoundaryConditions();}// 视情况输出runTime.write();// 显示计算时间、计算耗时Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"<< "  ClockTime = " << runTime.elapsedClockTime() << " s"<< nl << endl;}// 显示终止Info<< "End\n" << endl;return 0;}// ************************************************************************* //

OpenFOAM学习笔记_01_icoFoam理解相关推荐

OpenFOAM学习笔记案例1之Cavity（1）
OpenFOAM照比商业软件如Fluent Star-CCM 较难,是因为其操作需要较高的编程能力.而我不具备这种能力,只能通过一个一个的案例来学习,也希望对其他OpenFOAM新手提供一些帮助. 在 ...
openFOAM学习笔记（一）—— C++基础
一.C++基础很久不用C++,学习openFOAM之前复习一下C和C++相关的基础知识 1.1常用操作符 ++,– 整型变量的自加自减,用法很多.最简单的有i++,i–.运算速度会比i=i+1更快. ...
openFOAM学习笔记(五)——chemFoam的运行过程
在前面的帖子中已经大概给出了chemFoam主程序的结构,这里给出一个比较全面的总结首先程序结构如下: 添加头文件//*****************************// int main ...
openFOAM学习笔记（二）—— openFOAM的安装和网络资料汇总
openFOAM的安装安装参考了这三篇帖子: http://dyfluid.com/docs/install.html https://blog.csdn.net/u011786352/articl ...
TCP/IP学习笔记-如何理解
任何技术的掌握都需要做到应用技能的熟练掌握,比如让你写一个实现亮灯的程序,你本能的知道加载头文件,写main函数,这就是熟练掌握的应用技能,让一个刚学C的人,肯定就会为为什么家在头文件,为什么要写ma ...
tipi 深入理解php内核 pdf_大牛的学习笔记-深入理解Linux内核(完整版)
第一章.绪论 1.Unix文件可以是下列类型之一: a.正规文件(regular file) b.目录(directroy) c.符号链(symbolic link) d.块设备文件(block-or ...
JavaScript --- [学习笔记]观察者模式理解对象工厂模式构造函数模式
说明本系列(JS基础梳理)为后面TCP的模拟实现做准备本篇的主要内容: 观察者模式.工厂模式.构造函数模式和对对象的理解 1. 观察者模式参考JavaScript设计模式 1.1 消息注册方 ...
SpringMVC:学习笔记(1)——理解MVC及快速入门
SprigMVC-理解MVC及快速入门说明: 传统MVC-->JSPModel2-->Front Controller + Application Controller + Page C ...
OpenFOAM 学习笔记
目录目录 OpenFOAM在Centos7系统集群上OpenMPI并行的配置 1)环境说明 2)建立节点间SSH无密码连接 3)NFS的配置(已安装NFS) (1)管理节点配置 (2)计算节点配置 ...

OpenFOAM学习笔记_01_icoFoam理解

1 控制方程

2 算法流程

3 代码分析

OpenFOAM学习笔记_01_icoFoam理解相关推荐

最新文章

热门文章