• 返回主页
  • potentialFoam解析

    Reviewed by:
    求解器类型:速度势求解器,可用于粗略预测速度场。

    potentialFoam求解速度势方程(无粘无旋)后写入速度场。其在OpenFOAM主要用于粗略设置初始速度场,或使得速度通量守恒,为后续的使用其他求解器的计算加快收敛速度。其植入的控制方程如下:

    \begin{equation} \nabla \cdot u =0 \end{equation} \begin{equation} \nabla^2 \phi =0 \end{equation} 其中$\phi$为速度势,在OpenFOAM的代码中用p来表示,注意不要当做压力。下面进入代码分析:

    首先是createFields.H

        Info<< "Reading field p\n" << endl;
        volScalarField p
        (
            IOobject
            (
                "p",
                runTime.timeName(),
                mesh,
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            ),
            mesh
        );
    
        p = dimensionedScalar("zero", p.dimensions(), 0.0);//初始化速度势为0
    
    
        Info<< "Reading field U\n" << endl;
        volVectorField U
        (
            IOobject
            (
                "U",
                runTime.timeName(),
                mesh,
                IOobject::MUST_READ,
                IOobject::AUTO_WRITE
            ),
            mesh
        );
    
        U = dimensionedVector("0", U.dimensions(), vector::zero);//初始化速度为0
    
        surfaceScalarField phi
        (
            IOobject
            (
                "phi",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            fvc::interpolate(U) & mesh.Sf()
        );//计算phi,由于速度初始化为0,因此phi的内部场均为0
    
        if (args.optionFound("initialiseUBCs"))
        {
            U.correctBoundaryConditions();
            phi = fvc::interpolate(U) & mesh.Sf();
        }
    
        label pRefCell = 0//备注(1);
        scalar pRefValue = 0.0//设定参考压力为0;
        setRefCell
        (
            p,
            potentialFlow,
            pRefCell,
            pRefValue
        );//备注(2)
    	
    


    备注(1):此处设置参考速度势网格单元为0。对于某些算例并没有提供速度势固定值边界条件,这对于求解速度势泊松方程来讲会导致无解。例如:如果我们考虑一维速度势泊松方程: \begin{equation} \frac{\partial}{\partial x} (\frac{\partial p}{\partial x}) = 0 \end{equation} 其解为$p=ax+b$。如果定义速度势固定值边界条件$p(0)=0$,$p(1)=1$,其中$a$和$b$的值就可以通过他们求出。如果我们制定两个零梯度边界条件我们有:$\frac{\partial p}{\partial x} = 0$,通过俩个零梯度速度势边界条件可以求出$a=0$,但是我们并求不出来$b$的值。可以看出来如果速度势的边界条件均为零梯度边界条件,速度势泊松方程无解。因此我们需要设定参考速度势。

    备注(2):在OpenFOAM早期的版本,对于没有提供固定值边界条件的速度势条件,在求解器中强制执行参考速度势单元为0以及参考速度势值为0(或者类似的值)。因此并不需要setRefCell()函数,在OpenFOAM目前的版本中,用户具有更大的灵活,可以通过上述代码自行设置参考速度势单元以及参考速度势值。另外serRefCell()函数的一个作用为:假如用户设定的参考速度势网格单元碰巧是cyclic、sysmetry、processor边界条件,则取最近的非上述边界条件的网格单元。

    下面进入主程序:

    #include "fvCFD.H"
    #include "fvIOoptionList.H"
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    int main(int argc, char *argv[])
    {
        argList::addBoolOption("writep", "write the final pressure field");//类似的addBoolOption可以在运行程序的时候添加附加命令参数例如:
    potentialFoam –writep,这将写入速度势场
    argList::addBoolOption ( "initialiseUBCs", "initialise U boundary conditions" ); #include "setRootCase.H" #include "createTime.H" #include "createMesh.H" #include "readControls.H" #include "createFields.H" #include "createFvOptions.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< nl << "Calculating potential flow" << endl; // Since solver contains no time loop it would never execute // function objects so do it ourselves runTime.functionObjects().start(); fvOptions.relativeFlux(phi); adjustPhi(phi, U, p);


    上部分代码中adjustPhi(phi,U,p)遍历边界条件并进行通量加和以确保通量守恒。对于不守恒的算例,修改压力边界条件为zeroGradient边界条件的通量使其守恒(例如出口通量等于进口通量)。如果不守恒,压力方程无解并会提示:

    Continuity error cannot be removed by adjusting the"
    
    " outflow.\nPlease check the velocity boundary conditions"
                               
    " and/or run potentialFoam to initialise the outflow."
    


    对于某些算例(例如封闭体系空腔流),此函数自动返回为真,继续代码分析:

        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)//此处进入迭代求解,备注(3)
        {
            fvScalarMatrix pEqn
            (
                fvm::laplacian
                (
                    dimensionedScalar
                    (
                        "1",
                        dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
                        1
                    ),
                    p//此处dimTime的单位为(0, 0, 1, 0, 0), 首先除以从文件中读取的速度势单位(0, 2, -2, 0, 0),
    再乘以dimensionSet(0, 2, -2, 0, 0),再乘从文件中读取的速度势单位(0, 2, -2, 0, 0),最后进行梯度操作后单位变为(0, 0, -1, 0, 0)
    ) == fvc::div(phi)//此处div(phi)的单位为(0, 0, -1, 0, 0),因为在div()函数执行的时候,在函数内除以了体积的单位,且由于通量的单位为 $m^3 s^{-1}$, );//因此div(phi)的单位为(0, 0, -1, 0, 0),方程左右量纲一致,此处代码对应公式(2) pEqn.setReference(pRefCell, pRefValue);//在构建速度势矩阵后设定参考速度势值 pEqn.solve();//求解速度势泊松方程,同时提升速度势边界条件,因此我们在代码中不需要p.correctBoundayConditions() if (nonOrth == nNonOrthCorr) { phi -= pEqn.flux();//从最后收敛的速度势场组建守恒通量场,公式为  $\phi=\nabla p \cdot S_f$ } } fvOptions.absoluteFlux(phi); Info<< "continuity error = " << mag(fvc::div(phi))().weightedAverage(mesh.V()).value() << endl; U = fvc::reconstruct(phi);//从守恒的通量场组建速度场 U.correctBoundaryConditions();//修正边界条件函数一般用于下述情况:对于求解的方程组,如果调用了solve()函数,则不需要调用此函数,
    因为solve()函数内部已经调用修正边界条件函数,如果场不是通过solve()进行求解,比如本例中的U为通过重构phi而来,则需要调用“修正边界
    条件函数”来修正边界条件,这可以确保边界条件是一致的 Info<< "Interpolated U error = " << (sqrt(sum(sqr((fvc::interpolate(U) & mesh.Sf()) - phi))) /sum(mesh.magSf())).value() << endl; // Force the write U.write(); phi.write(); if (args.optionFound("writep")) { p.write(); } runTime.functionObjects().end(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; Info<< "End\n" << endl; return 0; }

    如何引用:

    李东岳. potentialFoam解析.[OL]. http://dyfluid.com/potentialFoam.html, 2016.3.10

    东岳流体
    info@dyfluid.com