pimpleFoam源码解析

写柴火的小火柴

于 2024-05-27 20:35:29 发布

阅读量263

点赞数 8

文章标签：算法 linux

本文链接：https://blog.csdn.net/BeanChai/article/details/139246903

版权

本文对OpenFOAM-7中的pimpleFoam求解器进行了注释，鉴于其与icoFoam相似，仅对pimple算法与ico算法不同的地方进行了详细注释

若不懂icoFoam请先行阅读这篇推文：icoFoam源码解析（附测试用的源代码文件）-CSDN博客

pimple算法理论参考：不可压稳态SIMPLE算法 — OpenFOAM|CFD

这里我认为pimple算法就是PISO和SIMPLEC算法的结合，再PISO算法的基础上，对省去的临点项进行修正。

推荐参考博客：

OpenFOAM学习笔记_openfoam欧拉离散方案-CSDN博客

pimpleFoam 代码浅析 - 哔哩哔哩

主函数：

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2019 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
    pimpleFoam

Description
    Transient solver for incompressible, turbulent flow of Newtonian fluids,
    with optional mesh motion and mesh topology changes.

    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "CorrectPhi.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "postProcess.H"

    #include "setRootCaseLists.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"
    #include "initContinuityErrs.H"
    #include "createDyMControls.H"
    #include "createFields.H"//与icofoam不同的是，这里会读取湍流模型参数，以及"constant/fvOptions"中应用的函数
    #include "createUfIfPresent.H"

    
    turbulence->validate();

    if (!LTS)
    {
        #include "CourantNo.H"
        //Info<< "**************************" << nl << endl;
        #include "setInitialDeltaT.H"//这里会读取controlDict中的函数,但看源代码好像没有关于读取函数的语句？？？？？？
        
    }

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
        #include "readDyMControls.H"

        if (LTS)
        {
            #include "setRDeltaT.H"
        }
        else
        {
            #include "CourantNo.H"
            #include "setDeltaT.H"
        }

        runTime++;

        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
            {
                mesh.update();
                if (mesh.changing())
                {
                    MRF.update();

                    if (correctPhi)
                    {
                        // Calculate absolute flux
                        // from the mapped surface velocity
                        phi = mesh.Sf() & Uf();

                        #include "correctPhi.H"

                        // Make the flux relative to the mesh motion
                        fvc::makeRelative(phi, U);
                    }

                    if (checkMeshCourantNo)
                    {
                        #include "meshCourantNo.H"
                    }
                }
            }

            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                //Info<< "==========================" << nl << endl;
                laminarTransport.correct();
                turbulence->correct();
            }
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
        
        //return 0;
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

速度方程

// Solve the Momentum equation

MRF.correctBoundaryVelocity(U);

tmp<fvVectorMatrix> tUEqn
(
    //这里对流项fvm::div(phi, U)为隐式方法，如果改为fvc则为显式，显式需要库朗数或CFL数特别低
    //否则会出现严重的速度场振荡，这一点在槽道流中特别明显
    fvm::ddt(U) + fvm::div(phi, U)
  + MRF.DDt(U)//关于旋转坐标系的科里奥利力
  + turbulence->divDevReff(U)
 ==
    fvOptions(U)//源项
);
fvVectorMatrix& UEqn = tUEqn.ref();

UEqn.relax();

fvOptions.constrain(UEqn);

if (pimple.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p));

    fvOptions.correct(U);
}

压力方程求解 *********核心***********

//rAU 为矩阵对角元素
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf))
);

MRF.makeRelative(phiHbyA);

if (p.needReference())
{
    fvc::makeRelative(phiHbyA, U);
    adjustPhi(phiHbyA, U, p);
    fvc::makeAbsolute(phiHbyA, U);
}

tmp<volScalarField> rAtU(rAU);

//这里用的是SIMPLEC算法，参考东岳流体的SIMPLE算法一节http://www.dyfluid.com/simple.html
if (pimple.consistent())
{
    //rAtU就是方程37柏松方程中左侧系数，UEqn.H1()=-(矩阵当前行的非对角元素合)
    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
    //phiHbyA就是方程38的HbyA_simplec 是面通量(面上的值*面积)
    phiHbyA +=
        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
    //phiHbyA就是方程38的HbyA_simplec 是体心值
    HbyA -= (rAU - rAtU())*fvc::grad(p);
}

if (pimple.nCorrPiso() <= 1)
{
    tUEqn.clear();
}

//后面的内容与icoFoam一致
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAtU(), MRF);

// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
    fvScalarMatrix pEqn
    (
        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
    );

    pEqn.setReference(pRefCell, pRefValue);

    pEqn.solve();

    if (pimple.finalNonOrthogonalIter())
    {
        phi = phiHbyA - pEqn.flux();
    }
}

#include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector
p.relax();

U = HbyA - rAtU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);

// Correct Uf if the mesh is moving
fvc::correctUf(Uf, U, phi);

// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);

学习用的pimpleFoam求解器pimpleFoam_learn下载链接：

https://download.csdn.net/download/BeanChai/89365394?spm=1001.2014.3001.5503