potentialFoam

基本介绍

1. 基本概念：势流描述了速度场是一个标量函数的梯度，它称为速度势，势流常表征为无旋流动。
2. 在不可压缩流动中，势流速度遵循拉普拉斯方程
3. 应用范围：机翼的外部流场、水波、电渗流和地下水流动。对于强涡度效应的流动(或其部分)，势流近似是不适用的。

描述

$$
\mathbf v=\nabla \varphi\tag{1}
$$

由定义知梯度的旋度为0

东岳流体，张量公式（18）

$$
\nabla\times\mathbf v=\nabla\times\nabla\varphi=0\tag{2}
$$

不可压缩流

在不可压缩流动中，速度的散度为0，即
$$
\nabla\cdot \mathbf v=0\tag{3}
$$
代入(1)得
$$
\nabla^2\varphi=0\tag{4}
$$
$\nabla^2=\nabla\cdot\nabla=\Delta$：称为拉普朗普算子

二维流动分析

令$z=x+iy$和$w=\varphi+i\psi$

我们定义投影f

$$
f(x+iy)=\varphi+i\psi\quad or \quad f(z)=w\tag{5}
$$
（5）满足**Cauchy–Riemann equations**
$$
\frac{\partial \varphi}{\partial x}=\frac{\partial \psi}{\partial y},\quad \frac{\partial \varphi}{\partial y}=-\frac{\partial \psi}{\partial x}\tag{6}
$$

$$
\frac{df}{dz}=u-iv
$$

$$
u=\frac{\partial \varphi}{\partial x}=\frac{\partial \psi}{\partial y},\quad v=\frac{\partial \varphi}{\partial y}=-\frac{\partial \psi}{\partial x}\tag{7}
$$

$\varphi$和$\psi$满足拉普拉斯方程
$$
\Delta \varphi=\frac{\partial^2\varphi}{\partial x^2}+\frac{\partial^2\varphi}{\partial y^2}=0\tag{8}
$$

$$
\Delta \psi=\frac{\partial^2\psi}{\partial x^2}+\frac{\partial^2\psi}{\partial y^2}=0\tag{9}
$$

又
$$
\nabla \varphi\cdot \nabla \psi=\frac{\partial \varphi}{\partial x}\frac{\partial \psi}{\partial x}+\frac{\partial \varphi}{\partial y}\frac{\partial\psi}{\partial y}=\frac{\partial \psi}{\partial y}\frac{\partial \psi}{\partial x}-\frac{\partial \psi}{\partial x}\frac{\partial\psi}{\partial y}=0\tag{10}
$$
由此可知当势函数线与流函数线相垂直

求解过程（参考东岳流体分析）

定义速度势函数为：
$$
-\nabla\varphi=\mathbf U\tag{11}
$$

如此定义可保证速度势和压力的作用同向

速度势满足的方程为（4），将求解方程演变为
$$
\nabla \cdot (\nabla P)=-\nabla \cdot \mathbf U\tag{12}
$$
同时定义通量场
$$
\phi=\mathbf U_f\cdot \mathbf S_f\tag{13}
$$
并将其内部场初始化为0，同时通过$\mathbf U$给定$\mathbf \phi$的边界条件。这种方法可以将边界条件的影响通过源项来植入方程(12)。因为在方程离散的过程中，边界条件控制也是通过进入方程的源项影响待求变量的结果。

对于动量方程，假设稳态，忽略扩散项
$$
\nabla\cdot\mathbf U=-\nabla \frac{p}\rho\tag{14}
$$

代码分析

代码知识点

1. addOption, addBoolOption

//- Add to an option to validOptions with usage information
//  An option with an empty param is a bool option
static void addOption
(
    const word& opt,
    const string& param = "",
    const string& usage = ""
);

//- Add to a bool option to validOptions with usage information
static void addBoolOption
(
    const word& opt,
    const string& usage = ""
)
{
    addOption(opt, "", usage);
}

2. $\phi$
   //可以将体积场转化为表面场（运用fvc::interpolate())
   //或者由表面场转化为体积场（运用fvc::reconstruct())进行计算。

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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
    potentialFoam

Description
    Potential flow solver which solves for the velocity potential
    from which the flux-field is obtained and velocity field by reconstructing
    the flux.

    This application is particularly useful to generate starting fields for
    Navier-Stokes codes.

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

#include "fvCFD.H"				//有限体积法头文件
#include "fvIOoptionList.H"		//IOoptionList

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

int main(int argc, char *argv[])
{
    //增加新的有效输入参数。
    argList::addOption
    (
        "pName",
        "pName",
        "Name of the pressure field"
    );

    argList::addBoolOption
    (
        "initialiseUBCs",
        "Initialise U boundary conditions"
    );

    argList::addBoolOption
    (
        "writePhi",
        "Write the velocity potential field"
    );

    argList::addBoolOption
    (
        "writep",
        "Calculate and write the pressure field"
    );

    argList::addBoolOption
    (
        "withFunctionObjects",
        "execute functionObjects"
    );

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    //#include "readControls.H"
    const dictionary& potentialFlow =
        mesh.solutionDict().subDict("potentialFlow");	
    //寻找fvSolution中的potentialFlow关键字部分并设为字典

    const int nNonOrthCorr =
        potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
    //读取potentialFlow中的非正交修正的值，如没有，默认为0

    //#include "createFields.H"
    Info<< "Reading velocity 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);
    //Construct given a name, a value and its dimensionSet.
    //dimensioned(const word&, const dimensionSet&, const Type);
    //初始化速度场。这里初始化速度为0 ，如果初始化Ux=1，Uy=2，Uz=3 （单位为国际单位）可采用 
	//U = dimensionedVector("0", U.dimensions(), vector(1,2,3));

    //表面场，phi，界面流率，存储在体之间的交接面上。表面场(surface...)不能和体积场(vol...)
	//直接计算，因为他们存储在不同地方，大小不同。
	//可以将体积场转化为表面场（运用fvc::interpolate())
	//或者由表面场转化为体积场（运用fvc::reconstruct())进行计算。
    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::interpolate(U) & mesh.Sf()
    );

    if (args.optionFound("initialiseUBCs"))
    {
        U.correctBoundaryConditions();
        phi = fvc::interpolate(U) & mesh.Sf();
    }


    // Default name for the pressure field
    word pName("p");

    // Update name of the pressure field from the command-line option
    args.optionReadIfPresent("pName", pName);

    // Infer the pressure BCs from the velocity BCs
    wordList pBCTypes
    (
        U.boundaryField().size(),
        fixedValueFvPatchScalarField::typeName
    );

    forAll(U.boundaryField(), patchi)
    {
        if (U.boundaryField()[patchi].fixesValue())
        {
            pBCTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
        }
    }

    Info<< "Constructing pressure field " << pName << nl << endl;
    volScalarField p
    (
        IOobject
        (
            pName,
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar(pName, sqr(dimVelocity), 0),	//m^2/s^-2
        //对于不可压缩里面的压力通常为p/rho
        pBCTypes
    );
/*
template<class Type, template<class> class PatchField, class GeoMesh>
Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
(
    const IOobject& io,
    const Mesh& mesh,
    const dimensionSet& ds,
    const word& patchFieldType
)
:
    DimensionedField<Type, GeoMesh>(io, mesh, ds, false),
    timeIndex_(this->time().timeIndex()),
    field0Ptr_(NULL),
    fieldPrevIterPtr_(NULL),
    boundaryField_(mesh.boundary(), *this, patchFieldType)
{
    if (debug)
    {
        Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
               "creating temporary"
            << endl << this->info() << endl;
    }

    readIfPresent();
}
*/            

    Info<< "Constructing velocity potential field Phi\n" << endl;
    volScalarField Phi
    (
        IOobject
        (
            "Phi",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionedScalar("Phi", dimLength*dimVelocity, 0),	//m^2*s^-1
        p.boundaryField().types()
    );

    label PhiRefCell = 0;
    scalar PhiRefValue = 0;
     //只有求解区域所有的压力边界都为第二类边界条件时，上面的值才会用到。如果有第一类边界条件，
	//压力参考值为这点边界值。对于不可压缩流动压力值为相对值，上面的参考值的大小对结果无影响。
    setRefCell
    (
        Phi,
        potentialFlow,
        PhiRefCell,
        PhiRefValue
    );
/*
//- If the field needs referencing find the reference cell nearest
//  (in index) to the given cell looked-up for field, but which is not on a
//  cyclic, symmetry or processor patch.
void setRefCell
(
    const volScalarField& field,
    const dictionary& dict,
    label& refCelli,
    scalar& refValue,
    const bool forceReference = false
);
*/

    //#include "createFvOptions.H"
    fv::IOoptionList fvOptions(mesh);    

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

    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.makeRelative(phi);

    adjustPhi(phi, U, p);

    // Non-orthogonal velocity potential corrector loop
    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {
        fvScalarMatrix PhiEqn
        (
            fvm::laplacian(dimensionedScalar("1", dimless, 1), Phi)
         ==
            fvc::div(phi)//速度离散，显示。因为是压力方程，其他变量只能显示
        );

        PhiEqn.setReference(PhiRefCell, PhiRefValue);
        //设置压力参考值
        PhiEqn.solve();
        //求解压力方程，调用的fvMatrix成员函数

        if (nonOrth == nNonOrthCorr)
        {
            phi -= PhiEqn.flux();	//流量修正
        }
    }

    fvOptions.makeAbsolute(phi);

    //提示连续性方程求解误差
    Info<< "Continuity error = "
        << mag(fvc::div(phi))().weightedAverage(mesh.V()).value()
        << endl;

    //根据表面场重建速度场
    U = fvc::reconstruct(phi);
    //对速度场边界进行修正，主要针对第二类边界条件下的边界场
    U.correctBoundaryConditions();

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

    // Write U and phi
    U.write();
    phi.write();

    // Optionally write Phi
    if (args.optionFound("writePhi"))
    {
        Phi.write();
    }

    // Calculate the pressure field
    if (args.optionFound("writep"))
    {
        Info<< nl << "Calculating approximate pressure field" << endl;

        label pRefCell = 0;
        scalar pRefValue = 0.0;
        setRefCell
        (
            p,
            potentialFlow,
            pRefCell,
            pRefValue
        );

        // Calculate the flow-direction filter tensor
        volScalarField magSqrU(magSqr(U));	//U*U
        volSymmTensorField F(sqr(U)/(magSqrU + SMALL*average(magSqrU)));

        // Calculate the divergence of the flow-direction filtered div(U*U)
        // Filtering with the flow-direction generates a more reasonable
        // pressure distribution in regions of high velocity gradient in the
        // direction of the flow
        volScalarField divDivUU
        (
            fvc::div
            (
                F & fvc::div(phi, U),
                "div(div(phi,U))"
            )
        );

        // Solve a Poisson equation for the approximate pressure
        //由于divDicUU是其他变量，所以为fvc显式
        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
        {
            fvScalarMatrix pEqn
            (
                fvm::laplacian(p) + divDivUU
            );

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

        p.write();
    }

    runTime.functionObjects().end();

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

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

    return 0;
}


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

势流绕流圆柱的解析解（待验证）

$$
\begin{cases}
u=U_\infin [1-(\frac r d)^2cons(2\theta)] \
v=U_\infin(\frac r d )^2sin(2\theta)\
\end{cases}
$$

参考文献

1. http://blog.sina.com.cn/s/blog_5fdfa7e60100d5ic.html
2. http://dyfluid.com/potentialFoam.html
3. https://en.wikipedia.org/wiki/Potential_flow