[转载]solver>>basic>>potentialFoam的说明

su_junwei

//creatField.H

//提示创建压力场

Info<< "Reading field p\n" << endl;

//压力场为标量场，网格中心存储变量（OpenFOAM用的是非结构化同位网格），下面为创建标量场压力，两个参数IOobject对象和网格对象mesh，IOobject主要从事输入输出控制
    volScalarField p
    (
        IOobject
        (
            "p",          //压力场初始文件的名字。
            runTime.timeName(), //位置，该位置由case中的system/controlDict中的startTime控制
            mesh,               //网格对象，主要从事对象注册，以便由runTime.write()控制输出
            IOobject::MUST_READ,//该对象由文件读出创建，因此，需要READ
            IOobject::NO_WRITE //不输出压力场
        ),
        mesh  //压力场所用的网格对象，在createMesh.H创建
    );

//压力场初始化为0，单位为上面文件中的单位。dimensionedScalar 为带单位的标量，初始化三个

//参数，名字，单位，数值。也可采用如下形式

//dimensionedScalar("zero",dimensionSet(1,-1,-2,0,0 ,0, 0),0.0);

//应当注意，OpenFOAM中的大部分case对动量的求解都是求解的速度场，压力单位初始化时候，

//一般为除去密度后的值.

//dimensionSet中有7个参数，他们依次为质量、长度、时间、温度、摩尔数、安培、坎德拉。



p = dimensionedScalar("zero", p.dimensions(), 0.0);   



//提示读入速度场

Info<< "Reading field U\n" << endl;
    //创建速度场，向量场，体心存储变量。

volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE  //自动写，根据runTime.write()或者U.write();
        ),
        mesh
    );



//初始化速度场。这里初始化速度为0 ，如果初始化Ux=1，Uy=2，Uz=3 （单位为国际单位）可采用

//U = dimensionedVector("0", U.dimensions(), vector(1,2,3));

U = dimensionedVector("0", U.dimensions(), vector::zero);



//表面场，phi，界面流率，存储在体之间的交接面上。表面场(surface...)不能和体积场(vol...)

//直接计算，因为他们存储在不同地方，大小不同。

//可以将体积场转化为表面场（运用fvc::interpolate())

//或者由表面场转化为体积场（运用fvc::reconstruct())进行计算。

surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
    //U是体积场，运用插值转为表面场和表面积场进行相乘来初始化流率

//mesh.Sf()返回网格交接面面积矢量。

    fvc::interpolate(U) & mesh.Sf()

);

//压力参考cell
    label pRefCell = 0;

//压力参考值
    scalar pRefValue = 0.0;
    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
   //只有求解区域所有的压力边界都为第二类边界条件是，上面的值才会用到。如果有第一类边界条件，

//压力参考值为这点边界值。对于不可压缩流动压力值为相对值，上面的参考值的大小对结果无影响。





potentialFoam.C

#include "fvCFD.H"


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

//主程序入口

int main(int argc, char *argv[])
{

//增加新的有效输入参数。

argList::validOptions.insert("writep", "");

//设置case目录

#   include "setRootCase.H"

//创建时间，对计算流程进行控制，主要是名为runTime的对象

#   include "createTime.H"

//创建网格对象mesh
#   include "createMesh.H"

//创建场变量，前面已经说过
#   include "createFields.H"

//读simple控制参数
#   include "readSIMPLEControls.H"

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

//计算势流提示

    Info<< nl << "Calculating potential flow" << endl;

//对流率进行调整，以保证方程的连续性

    adjustPhi(phi, U, p);

//网格非正交性循环，如果网格是正交的，可以设定nNonOrthCorr=1

  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {

//创建压力方程，该方程为标量稀疏矩阵类
        fvScalarMatrix pEqn
        (
            fvm::laplacian  //拉普拉斯离散，隐式
            (
                dimensionedScalar
                (
                    "1",
                    dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
                    1
                ),
                p
            )
         ==
            fvc::div(phi)//速度离散，显示。因为是压力方程，其他变量只能显示
        );

//设定压力参考值

        pEqn.setReference(pRefCell, pRefValue);

//求解压力方程，调用的fvMatrix成员函数
        pEqn.solve();

//流率修正，应当注意OpenFOAM中对压力本身求解，而非压力变化值。

//关于simple算法和PISO算法的OpenFOAM实现，会在simpleFoam及其icoFoam中详细说明。

if (nonOrth == nNonOrthCorr)
        {
            phi -= pEqn.flux();
        }
    }

//提示连续性方程求解误差

    Info<< "continuity error = "
        << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
        << endl;

//根据表面场重建速度场

    U = fvc::reconstruct(phi);

//对速度场边界进行修正，主要针对第二类边界条件下的边界场
    U.correctBoundaryConditions();

    Info<< "Interpolated U error = "
        << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi)))
          /sum(mesh.magSf())).value()
        << endl;

    // Force the write

//直接对速度进行输出
    U.write();

//界面流率输出。注意当前场存储在界面上，phi的大小（个数）和U的大小（个数）不相同的
    phi.write();

//如果用户计算是，让输出p，即输入了writep，则输出p

    if (args.options().found("writep"))
    {
        p.write();
    }

//提示执行时间，cpu时间

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

//提示程序结束

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

//返回0   

return(0);
}
转自OpenFOAM研究：http://blog.sina.com.cn/openfoamresearch

[ 本帖最后由 su_junwei 于 2009-5-2 18:29 编辑 ]