本案例演示利用OpenFOAM创建自定义插值格式。

1 文件结构

利用下面的命令创建文件结构：

run
foamNewApp demo15
cd demo15
touch demo15.H

文件结构如下图所示。

• 修改Make/files文件

demo15.C
LIB = $(FOAM_USER_LIBBIN)/libcustomScheme

• 修改Make/options文件

EXE_INC = 
    -I$(LIB_SRC)/finiteVolume/lnInclude 
    -I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = 
    -lOpenFOAM 
    -lfiniteVolume 
    -lmeshTools

2 程序代码

• 头文件demo15.H

#ifndef demo15_H
#define demo15_H
#include "surfaceInterpolationScheme.H"
#include "volFields.H" 
 
namespace Foam
{
    // 创建一个继承自surfaceInterpolationScheme的模板类
    template <typename T>
    class customScheme : public surfaceInterpolationScheme
    {
        const surfaceScalarField &faceFlux_;
    public:
        // 运行时类的信息
        TypeName("customScheme");
        //----------------构造函数-----------------------------------------
        // 从面通量进行初始化
        customScheme(const fvMesh &mesh, const surfaceScalarField &faceFlux):
        surfaceInterpolationScheme(mesh),
        faceFlux_(faceFlux)
        {
        }
 
 
        customScheme(const fvMesh &mesh, Istream &is) :
        surfaceInterpolationScheme(mesh),
        faceFlux_(mesh.lookupObject(word(is)))
        {
        }
 
 
        customScheme(const fvMesh &mesh, const surfaceScalarField &faceFlux, Istream &):
        surfaceInterpolationScheme(mesh),
        faceFlux_(faceFlux)
        {
        }
 
 
        // ---------------成员函数--------------------------------------
        // 实现基类的虚函数
        /*
        - 返回插值权重因子。
        注：关键函数，该函数定义相邻网格对该面上的值的贡献。它继承自SurfaceInterpolationScheme基类
        (在Surface eInterpolationScheme.H中)。由于基类中存在纯虚函数weights，
        因此每个派生类都必须自己实现此函数。
        在当前示例中，该插值格式用于处理scalarTransportFoam.C中定义的标量传输方程的div项：
                fvScalarMatrix TEqn
                (
                    fvm::ddt(T)
                  + fvm::div(phi, T)      <=== 本案例针对这一项进行处理
                  - fvm::laplacian(DT, T)
                 ==
                    fvOptions(T)
                );
        计算出的权重最终用于Surface InterpolationScheme.C
        使用插值方法，利用owner（P）和neighbour（N）的值来确定面值：
        sfi[fi] = lambda[fi]*vfi[P[fi]] + y[fi]*vfi[N[fi]];
        */
        virtual tmp weights(const GeometricField &) const
        {
            // 创建权重。注意使用临时场变量来匹配函数返回类型，这里采用线性插
            tmp weights_linear = this->mesh().surfaceInterpolation::weights();
            // 还可以通过获取位于正通量方向的网格值来创建迎风权重
            tmp weights_upwind = pos0(this->faceFlux_);
 
 
            // 存储以备后用(由于tmp结构的工作方式，一旦使用它们就会立即解除分配)。
            surfaceScalarField wl = weights_linear();
            surfaceScalarField wu = weights_upwind();
 
 
            // 将权重相加以形成混合格式
            const scalar fBlend = 0.2;
            tmp retWeights = weights_linear * (1.0 - fBlend) + weights_upwind * fBlend;
 
            //-----------------下面的代码只是演示，完全可以不要------------------
            // 实例化一个实际的标量场
            surfaceScalarField w = retWeights();
            scalarField phi = this->faceFlux_();
 
            label iFace = 10;
            const labelUList &owner = this->mesh().owner();
            const labelUList &neighbour = this->mesh().neighbour();
            const vectorField &Cf = this->mesh().faceCentres();
            const vectorField &C = this->mesh().cellCentres();
            const vectorField &Sf = this->mesh().faceAreas();
 
            Info << "Face " << iFace << " x=" << Cf[iFace] << " phi=" << phi[iFace]
                 << " (inferred velocity phi/A=" << (phi[iFace] / mag(Sf[iFace])) << ")"
                 << nl
                 << "x_owner=" << C[owner[iFace]] << " x_neighbour=" << C[neighbour[iFace]]
                 << nl
                 << "w_upwind=" << wu[iFace] << " w_linear=" << wl[iFace] << " w_total=" << w[iFace]
                 << nl << endl;
            //-----------------上面的代码没啥用，只是演示，完全可以删掉--------------------
 
            return retWeights;
        }
        void operator=(const customScheme &) = delete;
    };
}
#endif
 
 

• 源文件demo15.C

#include "demo15.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * //
namespace Foam
{
    makeSurfaceInterpolationScheme(customScheme);
}

编译生成动态库libcustomScheme.so。

3 测试案例

这里主要说明如何使用上面编译完成的动态库。

• 在system/controlDict文件中加载动态库

libs ("libcustomScheme.so");

• 在system/fvSchemes文件中修改散度项离散格式

divSchemes
{
    default         none;
    div(phi,T)      Gauss customSchmeme;
}

运行计算测试修改了离散方式后的的计算结果。这里使用的是一种混合格式，简单的问题里面感觉不出来有什么明显的差异，这里就懒得贴结果图了。不过道友们可以采用相同的套路构建更加复杂的离散格式。

总体上来说，将离散格式转化为代码还是很简单的，真正难的是构造出快准稳的离散格式。

（完毕）