icoFoam

基本介绍

1. 适用范围：用于不可压，层流，牛顿流体的瞬态求解器

2. 控制方程：

$$
\nabla\cdot\mathbf u=0\tag{1}
$$

$$
\rho\frac{\partial \mathbf u}{\partial t}+\rho\nabla\cdot{\partial (\mathbf u\mathbf u)}=-\nabla p+\mu\nabla\cdot(\nabla\mathbf u)\tag{2}
$$

在不可压缩流体的计算中，通常将压力定义为$p/\rho$

之后便是采用之前提到的PISO算法计算，此处不再赘述

代码分析

Application
    icoFoam

Description
    Transient solver for incompressible, laminar flow of Newtonian fluids.

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

#include "fvCFD.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"	//设置算例的根目录

    #include "createTime.H"		//创建时间对象
    #include "createMesh.H"		//创建网格对象
    
---------------------------------------createFields.H-------------------------    
    //#include "createFields.H"
    Info<< "Reading transportProperties\n" << endl;

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );

    //typedef dimensioned<scalar> dimensionedScalar;
    dimensionedScalar nu	
    (
        transportProperties.lookup("nu")
    );

    Info<< "Reading field p\n" << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );


    Info<< "Reading field U\n" << endl;
    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

-------------------------------createPhi.H-----------------------------------
//#   include "createPhi.H"'
    Info<< "Reading/calculating face flux field phi\n" << endl;
    
    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        linearInterpolate(U) & mesh.Sf()
        //Sf(): Return cell face area vectors
    );
--------------------------------createPhi.H-----------------------------------    


    label pRefCell = 0;
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
---------------------------------createFields.H------------------------------
   
---------------------------------initContinuityErrs.H-------------------------
    //#include "initContinuityErrs.H"
    scalar cumulativeContErr = 0;
---------------------------------initContinuityErrs.H-------------------------
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;
----------------------------------readPISOControls.H--------------------------
        //#include "readPISOControls.H"
        const dictionary& pisoDict = mesh.solutionDict().subDict("PISO");
        //fvSolution中的PISO字段
    
        const int nOuterCorr =
            pisoDict.lookupOrDefault<int>("nOuterCorrectors", 1);
    
        const int nCorr =
            pisoDict.lookupOrDefault<int>("nCorrectors", 1);
    
        const int nNonOrthCorr =
            pisoDict.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
    
        const bool momentumPredictor =
            pisoDict.lookupOrDefault("momentumPredictor", true);
    
        const bool transonic =
            pisoDict.lookupOrDefault("transonic", false);
----------------------------------readPISOControls.H--------------------------   
        
----------------------------------CourantNo.H--------------------------------
        //#include "CourantNo.H"
    	//库朗数的计算：参见scalarTransportFoam
        scalar CoNum = 0.0;
        scalar meanCoNum = 0.0;
        
        if (mesh.nInternalFaces())	//网格存在内部面
        {
            scalarField sumPhi
            (
                fvc::surfaceSum(mag(phi))().internalField()
            );
        
            CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
        
            meanCoNum =
                0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
        }
        
        Info<< "Courant Number mean: " << meanCoNum
            << " max: " << CoNum << endl;
-----------------------------------CourantNo.H--------------------------------
        fvVectorMatrix UEqn	//fvm：隐式求解
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        solve(UEqn == -fvc::grad(p));	//对应方程(5)

        // --- PISO loop

        for (int corr=0; corr<nCorr; corr++)
        {
            volScalarField rAU(1.0/UEqn.A());	//1/A_D

            volVectorField HbyA("HbyA", U);
            HbyA = rAU*UEqn.H();
            //H()代表方程（5）中A_{OD}和Q部分
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf())
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );
            //ddtCorr:用来考虑差值速度场计算的流律和实际流律的差别

            //- 调整边值，以便满足连续性条件
			//  For cases which do not have a pressure boundary.
			//  Return true if the domain is closed.
            adjustPhi(phiHbyA, U, p);

            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
            //压力修正步骤
            {
                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );
                //求解泊松方程（10）

                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve();

                if (nonOrth == nNonOrthCorr)
                {
                    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;
}


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

参考文献

[1] http://dyfluid.com/icoFoam.html

[2] http://blog.sina.com.cn/s/blog_5fdfa7e60100d6le.html