“cannot be called for a calculatedFvPatchField“——OpenFOAM中实现ρU(rhoU)边界条件的源码学习

因研究需要,特写一篇非单一变量(ρU, rhoU)边界条件的实现过程。在解可压NS方程时,rhoCentralFoam(解析)对动量方程的ρU进行直接的插值求解。换句话说,以ρU作为一个守恒变量,在方程中先进行求解,再分别更新ρ和U来。
我们知道, 有限体积法的边界条件就是计算边界面上的通量。在OpenFOAM中单一变量ρ和U的边界条件都必须要分别设置好,这意味着它们边界面上的通量如何计算将由我们来进行决定,其中最常用的有fixedValue(第一类边界条件)和zeroGradient(第二类边界条件)。那么问题来了,以ρU为方程中的独立变量进行求解,它的边界条件是怎么被设定好的?
这似乎是一个自动的过程,毕竟在0文件中也没有让我们设置ρU的边界条件。然而自己试着求解类似的方程时,却报错了,信息如下:

--> FOAM FATAL ERROR: 
cannot be called for a calculatedFvPatchField
    on patch CYLINDER of field f in file "/home/dyfluid/Desktop/cylinder/0.0002/f"
    You are probably trying to solve for a field with a default boundary condition.

翻译过来就是“您可能正在尝试求解具有默认边界条件的流场”。在“CYLINDER”这个壁面边界上,ρ为零梯度,U为无滑移,那ρU呢?我知道这个应该按无滑移处理,不过怎么在代码中实现呢?这意味着,类似于这种守恒变量的边界需要额外的处理,而在OpenFOAM中(如rhoCentralFoam)这些是替我们设置好了的,而我们并不清楚其中的细节。所以为了解决上面的问题,我们有必要学习一下源码中的实现方式。

rhoCentralFoam的实现方式

rhoCentralFoam求解的是无粘的NS方程(欧拉方程),使用显式离散进行求解。在rhoCentralFoam中,我们并不需要在0文件里面创建rho的文档并设置边界和初值,因为它在求解器里面通过psiThermo(psi是ψ,compressibility,即为ρ/p)进行创建的:

volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo.rho()
);

其中thermo就是psiThermo类型变量。在rhoCentralFoam中,和边界有关的代码:

 // --- Solve density
solve(fvm::ddt(rho) + fvc::div(phi));
 // --- Solve momentum
 solve(fvm::ddt(rhoU) + fvc::div(phiUp));

U.ref() =rhoU()/rho();
U.correctBoundaryConditions();
rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField();

rhoCentralFoam在求解完连续性方程(解ρ)动量方程(以ρU为守恒变量)之后,通过ρ导出了U。以correctBoundaryConditions()修正了U的边界条件,再对ρU的边界进行修正。
而自己仿照其形式写的却报错了(前面提到),那原因一定出在更深一层的代码上。先在网上看看有没有相似的报错信息,也许会找到有用的线索。在CFD中文网上有相关的帖子,帖子中的同学在计算p场时使用了calculated边界条件而报的一样错,而这是解方程时不被允许的(那什么时候能用calculated?)cfd-online上也有相似的问题,回答也是说用了calculated边界条件。那么,calculated边界条件是什么?我在0文件中给变量配置的边界条件里并没有设置calculated(只有fixedValue和zeroGradient),为什么也会报错?

calculated边界条件

查看calculatedFvPatchField.H,找到关于这个边界条件的描述:

This boundary condition is not designed to be evaluated; it is assmued that the value is assigned via field assignment, and not via a call to e.g. \\c updateCoeffs or \\c evaluate.
该边界条件的设计目的不是进行评估;而是假定该值是通过字段赋值而不是通过调用来赋值的,例如\\c updateCoeffs或\\c evaluate。

再看看fixedValue的:

This boundary condition supplies a fixed value constraint, and is the base class for a number of other boundary conditions.
此边界条件提供固定值约束,是许多其他边界条件的基类。

和zeroGradient:

This boundary condition applies a zero-gradient condition from the patch internal field onto the patch faces.
此边界条件以零梯度的方式将边界单元上的值赋值到边界面上。

至少从描述上来理解,我所采用的fixedValue和zeroGradient边界条件应该是不会调用到calculated的。
再看看FOAM aborting后面几行的信息:

FOAM aborting
#0 Foam::error::printStack(Foam::Ostream&) at ??:?
#1 Foam::error::abort() at ??:?
#2 Foam::calculatedFvPatchField<Foam::Vector>::valueInternalCoeffs(Foam::tmp<Foam::Field > const&) const in"/home/dyfluid/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt32Opt/bin/shenFoam"
#3 Foam::fv::gaussConvectionScheme<Foam::Vector>::fvmDiv(Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<Foam::Vector, Foam::fvPatchField, Foam::volMesh> const&) const at ??:?

注意到报错的函数是calculatedFvPatchField::valueInternalCoeffs, 且是在离散对流项格式时,可以看到#3有fvmDiv函数。也许问题不是出现在边界格式上,而是我对ρU这个场的创建过程。

创建ρU场

在rhoCentralFoam中ρU场(rhoU)是在createFields.h中直接以rho*U定义好的:

volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    thermo.rho()
);

volVectorField rhoU
(
    IOobject
    (
        "rhoU",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    rho*U
);

在这里,注意跟IO有关的字段“IOobject::NO_READ”。rho场和rhoU场都是不读取文件,这意味着,它们是凭空建立的场,初值、边界等都是继承于其他场(或自己设置)的。rho场由thermo.rho()返回得到,而rhoU场则是由rho*U得到,这里面必然有涉及到边界的继承。寻根索源,我们看一看thermo.rho()。在psiThermo.C里找到代码:

Foam::tmp<Foam::volScalarField> Foam::psiThermo::rho() const
{
    return p_*psi_;
}

在这里,psi_是psiThermo类中的成员变量,p_是basicThermo类中的成员变量(psiThermo->fluidThermo->basicThermo)。也就是说,rho是由压强和可压性决定的。p_在一开始读取时就安排好了,我们看到basicThermo.C的构造函数:

Foam::basicThermo::basicThermo
 (
    const fvMesh& mesh,
    const word& phaseName
)
 :
     IOdictionary
     (
        IOobject
         (
            phasePropertyName(dictName, phaseName),
             mesh.time().constant(),
             mesh,
             IOobject::MUST_READ_IF_MODIFIED,
             IOobject::NO_WRITE
         )
     ),
 
    phaseName_(phaseName),
  	p_(lookupOrConstruct(mesh, "p")),

而psi_在psiThermo中:

    psi_
    (
        IOobject
        (
            phasePropertyName("thermo:psi"),
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        dimensionSet(0, -2, 2, 0, 0)
    ),

一方面,p_是通过basicThermo读取了0文件中的p,另一方面,psi_由psiThermo创建。p是有边界和初始条件的(从文件中读取),但psi_看起来似乎是没有安排初值和边界的,这似乎是由thermophysicalProperties中设置thermoType之后由系统默认的。问题仍然没有解决,为什么p* psi* U(ρ*U)组成的场就没有问题了呢?

NO_READ!

除了必要了0文件的p和U场以外,前面几个场的建立(psi,rho,rhoU)等等都设置了不要读取文件。这似乎成为了某种误导。在某博文的列表403中提到了这个问题,并称这是因为“由于使用NO_READ标志创建该场,相当于没有边界条件被提供。”要想改正这个问题,不是把NO_READ改成MUST_READ(作者试过),而是需要删除创建该场时所传递的参数,在计算之前再进行赋值更新。例如:

volScalarField T 
(
IOobject 
(
"T",
runTime.timeName (), 
mesh ,
IOobject ::NO_READ , 
IOobject :: AUTO_WRITE
), 
mesh , 
dimensionedScalar("zero", dimensionSet (0, 0, 0, 0, 0), 0.0) 
);

第12行“dimension…”这个传递参数的一行应当删去,具体的原理作者还不能确定,但这与psi一样是没有自定义初始化的。在cfd-online另一篇帖子也提到这个问题,且有相似的解决方案:
在这里插入图片描述
所以我在算例的0文件夹中创建了f,并在求解器中创建场:

volScalarField rho
(
	IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh
);

volVectorField f
(
    IOobject
    (
        "f",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::NO_WRITE
    ),
    mesh
);

计算之前再赋值:

rho==rho0+rhop;
f==rhop*U+rho*up;

就不再报“cannot be called for a calculatedFvPatchField”的问题了!(还有其他的bug)。

总结

“cannot be called for a calculatedFvPatchField”一方面可能是因为边界使用了calculated的边界条件而报错,另一方面可能是自己创建的场在定义时,传递了别的参数,使得这与其中的某些设置起到了冲突。因此在求解新的变量场时(如ρU这种形式),需要注意定义时不初始化,计算时再赋值。

参考网址

  1. http://dyfluid.com/rhoCentralFoam.html)
  2. https://www.cfd-online.com/Forums/main/233758-error-cannot-called-calculatedfvpatchfield.html
  3. https://www.cfd-online.com/Forums/main/233758-error-cannot-called-calculatedfvpatchfield.html
  4. https://gitee.com/poplee/openFoamUserManual/blob/master/57.Under_the_hood_of_OpenFOAM.md
  5. https://www.cfd-online.com/Forums/openfoam-programming-development/107479-gradientinternalcoeffs-cannot-called-calculatedfvpatchfield-print.html
  6. OpenFoam c++ sourcs code guide

© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发
头像
欢迎您留下宝贵的见解!
提交
头像

昵称

取消
昵称表情代码图片