HbyA,phiHbyA,fvc::div(phiHbyA)计算错误问题
-
最近做pEqn.H中的演算，输出了HbyA，phiHbyA和fvc::div(phiHbyA)，通过手动计算进行验证
首先测试了一个简单的331的立方体模型，由HbyA计算phiHbyA，再进行求和并除以体积求div，计算结果和输出的结果一致，没问题。
换了个模型，吊诡的事情出现了，由HbyA计算phiHbyA没问题，但是由phiHbyA进行求和算div居然是错的？！我很奇怪难道连加法也能算错？
以下是我的一个结果，测试的cell编号是47794：
手动计算的div：-6.3169e-3
输出的div：1.90343e-6且不说正负号，cell体积的量级就是10^-8，手动求和phiHbyA的量级是10^-10，check了一下代码，也没有其它干扰了，实在想不通内部是怎么算的
-
-
@李东岳 前辈好，是这样，我打算测试压力泊松方程源项对计算结果的影响
fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)以上是压力方程，我想看看**源项fvc::div(phiHbyA)**的计算结果和我自己手动计算的是否一样
所以下一步我需要求phiHbyA，我看了fvc::surfaceIntegrate的源码，fvc::div(phiHbyA)的计算就是将一个cell中不同face的通量phiHbyA相加，之后除以体积mesh.V()
关于phiHbyA的计算，查看代码如下，
surfaceScalarField phiHbyA ( "phiHbyA", fvc::flux(HbyA) + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi) );有两个部分组成，一个是通过volVectorField类型的HbyA中心插值到face上再与面矢量做点积，第二个修正项是通过上一时间步对本时间的修正，具体算法代码中有
目前采用的都是2D的算例，一个cell有4个neighbour，具体细节如下：
对于简单算例，几何模型如下，是一个3 * 3 * 1的一个立方体，选择4号单元进行验证，具体数据如图
由HbyA手动计算得得phiHbyA与输出的结果一致，再由phiHbyA计算fvc::div(phiHbyA)也没问题；
对于网格量比较大的复杂问题
由HbyA手动计算得得phiHbyA与输出的结果一致但是，由phiHbyA求和除以mesh.V()所得的结果，与输出的fvc::div(phiHbyA)相差了好几个量级！也就是之前列出来的结果了。
求前辈指教，我现在不懂为啥这个求和还能求错了？或者是求这个fvc::div(phiHbyA)时根本用的就不是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就是对速度做体积分，通过加和的方式计算。 -
-
@李东岳 前辈 您说的这个算例我没异议，这个就和我的那个简单例子一样的，我目前已经掌握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)了？？
-
-
@李东岳 前辈 具体的计算过程如下，这里只说有错误的第二个算例（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
