GREAT = 1.0e+6;
VGREAT = 1.0e+37;
ROOTVGREAT = 1.0e+18;
SMALL = 1.0e-6;
VSMALL = 1.0e-37;
ROOTVSMALL = 1.0e-18

@李东岳 这个速度波动的问题，在fluent理论指南中找到如下说法：

20.4.1.1. 压力插值格式

注意：

对于一个包含多孔介质的算例，体积力加权格式只应用在非多孔面上，它考虑了显式体积力（例如，重力，旋转，科里奥利力）的不连续性和密度迅速变化流动（例如，自然对流，VOF）压力梯度的不连续性。所有内部和外部的多孔表面都使用一种特殊的格式处理，尽管阻力是不连续的，这种格式可以保持单元体表面上法向速度的连续性。

但是我没有找到fluent具体采用的到底是哪种特殊的压力插值格式。

另外OpenFOAM中需要用压力插值格式吗？暂时对不上号。

在 bug report 中，henry 说

Using

U = fvc::reconstruct(phi);

causes a HUGE momentum conservation error but if momentum conservation is not important in your cases then yes it does remove spurious "wiggles" of many kinds.

为啥会有动量守恒的问题呢？

还有如何解决上述 bug report 中报告的多孔介质与自由流界面处的速度波动问题？ 在相关的 帖子 中，说

This is not a problem.

Just change under-relaxation factor to 0.0001, and you will not have this oscillations.

OpenFOAM software is the best software among all CFD, and being unexperienced in CFD, does not meam that OpenFOAM has bugs, that should be reported like that...

这到底是bug 吗？解决方案就是小松弛因子？

@cfd-china 厉害了我的哥

move to Gas-dispersion-by-OpenFOAM

另外，
参考
fvc::div(u)=fvc::div(phi)。

对应的代码如下：

template<class Type>
tmp
<
    GeometricField
    <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>
>
gaussDivScheme<Type>::fvcDiv
(
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp
    <
        GeometricField
        <typename innerProduct<vector, Type>::type, fvPatchField, volMesh>
    > tDiv
    (
        fvc::surfaceIntegrate
        (
            this->mesh_.Sf() & this->tinterpScheme_().interpolate(vf)
        )
    );

    tDiv().rename("div(" + vf.name() + ')');

    return tDiv;
}

即div(U)的计算方法也是先插值到单元体表面上，然后进行求和。

correctBoundaryConditions()
{
 this->setUpToDate();
 storeOldTimes();
 boundaryField_.evaluate();
}

然后就是具体调用各个边界条件中的evaluate()函数了。

可参考http://www.cfd-online.com/Forums/openfoam-programming-development/129271-how-boundary-conditions-called-openfoam-solvers.html

但是，如果选了calculated边界，边界值就是被算出来的。但是这个值是怎么给定的我暂时不清楚，比如U=d*P，这都是针对体心值的操作，不清楚选了calculated时边界值是如何更新的

@cfd-china 数理方程 这门课讲了

mark

@hurricane007 试试，说不定能发个science

http://wiki.jikexueyuan.com/project/tensorflow-zh/tutorials/pdes.html

流不流弊

@程迪 厉害。确实有一些计算机方面的方法可以搞，不过用vtk解决了，多谢！

@李东岳 生成的随机坐标与前一个点的值有关

准备试试vtk，好像没有信息损失。。

发现paraview导出csv然后处理比较简单，虽然这样没有网格信息了，只能基于网格点上的值对任意点的T进行插值

目前做一个后处理算法，需要按照随机的位置对计算结果进行采样。后处理中需要用到复杂的算法，目前是基于python的类库实现。
最简单的想法是在python中调用probe，但是这样每采样一次就需要重新从硬盘读取整个网格的结果，太慢。
另一个想法是写一个of的后处理程序，首先读入整个网格的结果进内存，然后调用python代码，有个问题是此时python怎么访问读进来的volScalarField？即：

int main(int argc, char *argv[])
{
    #include "setRootCase.H"

    #include "createTime.H"
    #include "createMesh.H"

    Info<< "Reading field T\n" << endl;
    volScalarField T
    (
        IOobject
        (
            "T",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    调用一段python，里面有一个函数的功能是：对于python随机生成的坐标(x,y,z)，能够返回T[mesh.findCell(point(x,y,z))]
}

如果能在python中实现读入mesh和field，以及mesh.findCell功能，问题就简单了。

各位有何想法？多谢！

@piteqiu

    Info<< "Reading field S\n" << endl;
    volScalarField S  //1.48e-5 kg/s
    (
        IOobject
        (
            "S",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("S",dimensionSet(1 ,-3, -1 ,0 ,0 ,0 ,0),0.0)
    );
    forAll(S,i)
    {
        S[i]=0;
    }
    S[mesh.findCell(point(0,0,0))]=Q.value()/mesh.V()[mesh.findCell(point(0,0,0))];

@piteqiu 把这些物性参数都弄成场，然后随便咋赋值。。