HbyA,phiHbyA,fvc::div(phiHbyA)计算错误问题

李东岳

#include "fvCFD.H"

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

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

    volVectorField U
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

	surfaceScalarField phi
	(
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
       	fvc::interpolate(U) & mesh.Sf()
    );

	Info << "fvc::interpolate(U) " << fvc::interpolate(U) << nl;

	Info << "phi is " << phi << nl;

	Info << "mesh surface vector is " << mesh.Sf() << nl;

	Info << "fvc::div(phi) is " << fvc::div(phi) << nl;

    return 0;
}

请运行上述代码，并把U初始化为：

dimensions      [0 1 -1 0 0 0 0];

internalField   nonuniform List<vector> 
9
(
(0 1 0)
(0 2 0)
(0 3 0)
(0 4 0)
(0 5 0)
(0 6 0)
(0 7 0)
(0 8 0)
(0 9 0)
)
;


//internalField   uniform (0 1 0);

boundaryField
{
    wall
    {
        type            fixedValue;
        value           uniform (0 1 0);
    }
    defaultFaces
    {
        type            empty;
    }
}

我没看出计算错误的问题。fvc::div就是对速度做体积分，通过加和的方式计算。

李东岳

Yes, 我觉得这个算例非常好，因为OKSS2上来就是高斯积分计算，也算是我对你这个感兴趣的动力，我会把这个放在课上分析。

由phiHbyA求和除以mesh.V()所得的结果，与输出的fvc::div(phiHbyA)相差了好几个量级！

用上面那个简单试试，可能你那个太复杂了，小数点很多，数也不工整。客观来讲，不会出现这个现象。

yhdthu

@李东岳前辈您说的这个算例我没异议，这个就和我的那个简单例子一样的，我目前已经掌握HbyA到fvc::div(phiHbyA)的计算方法。

但是，问题来了，我用相同的计算方法去计算我的第二个例子的某个cell，发现求和的时候出现了错误。

即，我手动计算的fvc::div(phiHbyA)为：-6.3169e-3
但是输出的结果却是1.90343e-6

如果只看phiHbyA求和，其结果是10^-10量级，cell的体积是10^-7量级，两个相除就是10^-3次方量级了；

输出的fvc::div(phiHbyA)的结果居然是10^-6量级，这可是除完体积之后的值，也就是说之前的phiHbyA的求和量级是10^-13量级

前辈，这就是我想表达的意思，为什么从phiHbyA就算不到fvc::div(phiHbyA)了？？

李东岳

所以就是说：你实际上已经知道手动计算的方式，但是你算的和代码结果不一样。如果可以的话，请把4号网格4个面的phiHbyA输出，我看看你是怎样 手动计算的fvc::div(phiHbyA) 的，就是说：把你计算-6.3169e-3的过程写出来。

输出的fvc::div(phiHbyA)的结果居然是10^-6量级

这个非常正常，匹配你的压力方程残差。
你可以看出来，这只不过是加减法问题，可能你手动计算存在一个非常小的bug。

yhdthu

@李东岳前辈，我的意思是：

我计算那个方形的4号cell全都没问题，以上第一个方形网格的算例给出了我所查找的所有参数；

有问题的是，用相同的计算方法，去计算第二个算例的47794号cell，从phiHbyA到fvc::div(phiHbyA)的计算出了问题。

我明白量级什么的无所谓，最重要的是数得算的对，但是现在手动算的和程序算的根本就对不上了，而且差了n个量级，这就说不过去了:crying:

yhdthu

@李东岳前辈具体的计算过程如下，这里只说有错误的第二个算例（cell_47794）

需要进行计算的场输出为：
HbyAVol = HbyA;
phiHbyASuf = phiHbyA;
sourcePoissionU = fvc::div(phiHbyA);

根据输出的phiHbyA，4个面的值如下：

E：0.000125002
W：0.000124669
N：-1.94136e-05
S：-1.908e-05
F 和 B 为empty，0

加和：E - W + N - S = 6e-10

mesh.V() = 9.49835e-8

fvc:div(phiHbyA) = (E - W + N - S) / mesh.V() = -0.0063169

但是，输出的fvc::div(phiHbyA) = 1.90343e-6

李东岳

E：0.000125002
W：0.000124669
N：-1.94136e-05
S：-1.908e-05

这个找的正确么

yhdthu

@李东岳前辈，这个百分百没问题，我是用paraview一个一个面点出来看的

李东岳

你是怎么在paraview查看face的值的？

yhdthu

@李东岳通过paraview找到cell编号，查找polymesh里面的owner和neighbor文件定位行号，打开surfacefield对应找就好啦～

李东岳

方法是对的，不过为什么错了？我也很奇怪。

yhdthu

@李东岳更奇怪的是，程序前半段计算phiHbyA都是对的，但是之后求div的时候（加法）却错了，难道加法也能算错:upset:

yhdthu

@李东岳前辈，解决了，原来是输出精度不够导致的，我的锅:lol:

CFD中文网

HbyA,phiHbyA,fvc::div(phiHbyA)计算错误问题