pisoFoam计算的U，求div(U)结果为何不是严格等于0

李东岳

@coolhhh 在 pisoFoam计算的U，求div(U)结果为何不是严格等于0 中说：

为何of中最后输出phi不用fvc::flux(U)，而是用了phi = phiHbyA - pEqn.flux()？

                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );
                pEqn.setReference(pRefCell, pRefValue);
                pEqn.solve();
                phi = phiHbyA - pEqn.flux();

上面这几行代码里面，pEqn实际上是建立在sum(phi)= sum(phiHbyA - pEqn.flux())=0，也就是说，在sum(phi)= sum(phiHbyA - pEqn.flux())=0的时候，才有fvm::laplacian(rAU, p) == fvc::div(phiHbyA)这个形式。所以求解fvm::laplacian(rAU, p) == fvc::div(phiHbyA)，sum(phi)= sum(phiHbyA - pEqn.flux())=0。

coolhhh

@李东岳谢谢李老师，这样能理解，还有个小疑问。我们可以从U计算phi=fvc::flux(U)，那么如果已知phi，能不能反算出U？这样的话U与phi就能完全对应上，此时计算就能满足div(U)=div(phi)。但of中后面接着用U = HbyA - rAU*fvc::grad(p)计算U，才导致了这个不对应问题，而U才是我们最后要得到的结果，但这个结果却没能满足零散度。

李东岳

那么如果已知phi，能不能反算出U？

就是因为这个搞不出来，所以就只是phi守恒了。各种从phi到U的算法，都涉及到各种插值之类。所以你从守恒的phi，到U（这里面有一次插值），再从U到phi（又来一次插值），误差就出来了。

有一个算法是reconstruct(phi)，不过这个里面也存在数学假设。

按道理来说守恒的场就是phi场。我们不能说速度是守恒的。我们只能说通量是守恒的。类似只能说能量是守恒的，不能说温度是守恒的。

coolhhh

理解了，谢谢李老师

coolhhh

李老师，我做了简单的对比：

下面计算的div(U) 与 div(phi) 结果基本相等，区别是小数点后五六位
①直接求div(U)
②已知U -> phi = fvc::flux(U) -> div(phi)
接着对比reconstruct(phi)。编译mypisoFoam，把pEqn.H中计算U修改为

U = fvc::reconstruct(phi);
// U = HbyA - rAU*fvc::grad(p);

同样计算计算of2206自带一个算例：Decay of homogeneous isotropic turbulence，模拟结果：
（1）对比div(U) 与 div(phi)的均值和标准差。采用reconstruct(phi)，得到的div(U) 的均值和标准差要小一点。这里有个注意点是，of计算的div(phi) ，基本数值都很小，但会有一个比较大的数，下面的计算结果是剔除了这个最大数的计算结果。

（2）对比湍流特性。可以看出reconstruct(phi)结果要更差，还是原本of的U = HbyA - rAU*fvc::grad(p)计算准确

李东岳

对比湍流特性。可以看出reconstruct(phi)结果要更差，还是原本of的U = HbyA - rAU*fvc::grad(p)计算准确

对比你下面那两个图，我觉得你可以把这部分总结一下写在文章里面，当做appendix。这是2种计算方法，从方程上来说后者是一致的。从物理上来说前者的U看起来更加守恒。不过从图上来看，reconstruct耗散掉了一部分没有后者好。

你这个做的很细，值得一写。

coolhhh

谢谢李老师，文章已经发表了，只是好奇做个对比。在想的问题是，一直以来各种湍流生成方法都在追求实现生成零散度的湍流场U，我们直接构造的是0时刻的U，而不是phi。

根据测试结果看，如果严格让div(U)=0, 此时U -> phi = fvc::flux(U) -> div(phi)得到的div(phi)也基本等于0。但这种情况下0时刻的湍流特性会变化的比较大，结果会变差。

of在后面时间步计算输出的U实际上不满足div(U)=0，所以想法是0时刻提供的初始U也不需要强制要求div(U)=0，这样得到的结果可能才是更合理的。

李东岳

of在后面时间步计算输出的U实际上不满足div(U)=0，所以想法是0时刻提供的初始U也不需要强制要求div(U)=0，这样得到的结果可能才是更合理的。

你目前0时间步的U与phi，是前者散度为0，还是后者。

不过我感觉0时间步的问题不大，稍微步进一个时间步，后面就都是div(phi)=0了。

厉害啊，小伙子，我记得你好像做建筑领域的，现在在做DNS么

coolhhh

0时间步的U，是用湍流生成方法直接生成的，可以计算得到div(U)。然后就可以计算div(phi)： U -> phi = fvc::flux(U) -> div(phi)。测试结果是，这种方式计算的div(U)和div(phi)结果基本一样的，区别是小数点后五六位。此时div(U)=div(phi)=0

稍微步进一个时间步，of输出的U和phi，就会出现div(U)不等于0，div(phi)=0

是做建筑领域的，DNS搞不了，还是做的LES

李东岳

0时间步的U，是用湍流生成方法直接生成的，可以计算得到div(U)。然后就可以计算div(phi)： U -> phi = fvc::flux(U) -> div(phi)。测试结果是，这种方式计算的div(U)和div(phi)结果基本一样的，区别是小数点后五六位。此时div(U)=div(phi)=0

那这个算法很屌。不错。

CFD中文网

pisoFoam计算的U，求div(U)结果为何不是严格等于0