【OpenFOAM入门】How to run your first simulation

视频参考链接：

https://www.bilibili.com/video/BV1ga41157pp?p=2

版本：OpenFOAM 10
示例：$FOAM_TUTORIALS/incompressible/icoFoam/elbow

icoFoam

Transient solver for incompressible, laminar flow of Newtonian fluids (isothermal; PISO-loop)
pressure: $p=\frac{1}{\rho }\bar{p}$

Case setup

.
├── 0
│   ├── U
│   └── p
├── constant
│   └── physicalProperties
├── elbow.msh
└── system
    ├── controlDict
    ├── foamDataToFluentDict
    ├── fvSchemes
    └── fvSolution

elbow_tri/0/p：

dimensions      [0 2 -2 0 0 0 0]; // pressure divided by density, m^2/s^2

internalField   uniform 0;

elbow_tri/system/controlDict，改结束时间75s

endTime         75;

elbow_tri/system/fvSolution：

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 2;
}

PISO参数理解：cd /opt/openfoam10/applications/solvers/incompressible/icoFoam看源码：两层循环的迭代次数

#include "fvCFD.H"
#include "pisoControl.H"

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

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

    pisoControl piso(mesh);

    #include "createFields.H"
    #include "initContinuityErrs.H"

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

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

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.userTimeName() << nl << endl;

        #include "CourantNo.H"

        // Momentum predictor

        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

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

        // --- PISO loop
        while (piso.correct())
        {
            volScalarField rAU(1.0/UEqn.A());
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
            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
            while (piso.correctNonOrthogonal())
            {
                // Pressure corrector

                fvScalarMatrix pEqn // pressure equation of the PISO-loop
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );

                pEqn.setReference(pRefCell, pRefValue);

                pEqn.solve();

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

            #include "continuityErrs.H"
			
            U = HbyA - rAU*fvc::grad(p); // update all velocities with the new pressure
            U.correctBoundaryConditions();
        }

        runTime.write();

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

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

    return 0;
}

Info<< "Reading physicalProperties\n" << endl;

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

dimensionedScalar nu
(
    "nu",
    dimViscosity,
    physicalProperties.lookup("nu")
);

Info<< "Reading field p\n" << endl;
volScalarField p // pressure divided by density
(
    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
);


#include "createPhi.H" // creating the flux


label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solution().dict().subDict("PISO"), pRefCell, pRefValue);
mesh.schemes().setFluxRequired(p.name());

Initial values

elbow_tri/0

Mesh

三角形网格

导入网格

$ fluentMeshtoFoam elbow.msh
$ cd constant/polyMesh/
$ ls
boundary  cellZones  faceZones  faces  neighbour  owner  pointZones  points

• 查看几何外形：启动paraview，Open，选择All files，选择controlDict，选择OpenFOAMReader，Apply
  • solid color-surface with edges：查看网格
  • 选择mesh regions查看不同区域网格

四边形网格

下载网格文件elbow_quad.msh，其他步骤一致

网格细化

$ refineMesh -overwirte

OpenFOAM 10会有警告不知道如何处理，但simlulation可以继续。

--> FOAM Warning :
    From function void Foam::cellCuts::setFromCellCutter(const Foam::cellLooper&, const Foam::List<Foam::refineCell>&)
    in file meshCut/cellCuts/cellCuts.C at line 2393
    Found loop on cell 3319 that resulted in an unexpected bad cut.
    Suggestions:
      - Turn on the debug switch for 'cellCuts' to get geometry files that identify this cell.
      - Also keep in mind to check the defined reference directions, as these are most likely the origin of the problem.

网格细化了，为保持Courant Number小，时间间隔要相应减小。
system/controlDict：

stopAt          endTime;

endTime         75;

deltaT          0.025;

writeControl    timeStep;

writeInterval   40;

Simulation

$ icoFoam
Courant Number mean: 0.0810349 max: 0.501578
smoothSolver:  Solving for Ux, Initial residual = 1.29438e-05, Final residual = 3.2696e-08, No Iterations 1
smoothSolver:  Solving for Uy, Initial residual = 8.85889e-06, Final residual = 8.85889e-06, No Iterations 0
DICPCG:  Solving for p, Initial residual = 0.000425241, Final residual = 1.97738e-05, No Iterations 3
DICPCG:  Solving for p, Initial residual = 3.732e-05, Final residual = 1.53043e-06, No Iterations 47
DICPCG:  Solving for p, Initial residual = 1.25666e-05, Final residual = 9.64473e-07, No Iterations 4
time step continuity errors : sum local = 2.43042e-10, global = -3.45032e-12, cumulative = 4.61643e-07
DICPCG:  Solving for p, Initial residual = 8.68874e-05, Final residual = 4.14519e-06, No Iterations 8
DICPCG:  Solving for p, Initial residual = 1.40088e-05, Final residual = 8.68304e-07, No Iterations 15
DICPCG:  Solving for p, Initial residual = 4.2701e-06, Final residual = 7.87489e-07, No Iterations 1
time step continuity errors : sum local = 1.9844e-10, global = -4.48494e-11, cumulative = 4.61598e-07
ExecutionTime = 2.08762 s  ClockTime = 2 s

• Courant Number
• smoothSolver: 源码中求解动量方程：solve(UEqn == -fvc::grad(p));
• DICPCG: PISO-loop，外循环nCorrectors次，内循环nNonOrthogonalCorrectors+1次

Postprocessing

查看结果数据

$ pwd
/path/to/tutorials/incompressible/icoFoam/elbow_tri
$ cd 75
$ ls
U  p  phi  uniform
$ vim p # 查看p
dimensions      [0 2 -2 0 0 0 0];

internalField   nonuniform List<scalar>
918
(
......
)
;

boundaryField 
{
......
}

p不再是初值时的internalField uniform 0而是internalField nonuniform List<scalar>。每个值对应一个cell，共918个cell。最后boundaryField是边界条件。查看U的值则为向量。

paraview

• 启动paraview，Open，选择All files，选择controlDict，选择OpenFOAMReader，Apply；或者cp controlDict controlDict.foam，打开controlDict.foam
• legend
• timesteps

不同网格对比

• 左下角Transforming菜单移动位置
• 三角网格出口的速度更发散
• 四边形网格，能算出回流

结果格式

• foamToVTK将结果转为VTK格式，会多一个VTK/文件夹

结果导出

• 设置legend格式：标题、字体、范围、colormap、位置
• 保存截图：File->save screenshot
• 导出曲线：Filters->Data Analysis->Plot over line
  • 注意：曲线点的坐标不能根据contour移动后的坐标给，画了很久出不来线，最后关闭paraview重新打开controlDict.foam就可以了
  • 导出数据：File->Save Data