“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 calculatedFvPatchFieldon 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 momentumsolve(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