OpenFOAM学习笔记_01_icoFoam理解

最新推荐文章于 2024-11-21 10:52:54 发布

梅冠华

最新推荐文章于 2024-11-21 10:52:54 发布

阅读量3.4k

点赞数 11

分类专栏： OpenFOAM学习笔记文章标签： OpenFOAM icoFoam piso

本文链接：https://blog.csdn.net/meiguanhua/article/details/113678619

版权

OpenFOAM学习笔记专栏收录该内容

1 篇文章

订阅专栏

本文介绍了OpenFOAM中的icoFoam求解器，用于解决非定常不可压缩层流问题，采用PISO算法。控制方程包括连续方程和不可压缩Navier-Stokes方程，通过离散和求解动量方程、迭代更新压力和速度场来实现。详细步骤包括：离散动量方程、构建HbyA矩阵、求解Poisson方程、更新界面通量和速度场。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

参考资料：

东岳博文：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 控制方程

$\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}}})$

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

2 算法流程

PISO算法流程如下：

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

$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)}$

b. 基于新速度场 $\bold u_C^*$ ，构建 $\bold{HbyA}_C^\star$ ；

$\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)$

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

$\nabla \cdot (\bold {HbyA}^\star) = \nabla \cdot \left( \frac{1}{a_C^\bold u} \nabla p^\star \right)$

d. 根据新的压力场 $p^\star$ ，更新界面通量 $\bold u_f^\star$ 和速度场 $\bold u_C^\star$ ；

$\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$

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  |
 -------------------------------------------------------------------------------
 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.
 
 \*---------------------------------------------------------------------------*/
 
 // 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^U
             volScalarField 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. 压力修正
 
                 // 组装压力修正方程，获得标量方程pEqn
                 fvScalarMatrix 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^star
             U = 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;
 }
 
 
 // ************************************************************************* //