icoFoam的一些细节问题

sejabs

东岳大神写了关于icoFoam solver的一些解释，不过还有一些细节问题不太清楚，特来请教（OpenFOAM 3.0.1）：

这段代码：

while (piso.correct())
        {
            volScalarField rAU(1.0/UEqn.A());
            .....
        },

这里的UEqn.A和后面的UEqn.H是自动根据U的新值进行刷新的吗？
如果是固定不变的，那就没有必要写在while循环里；如果是自动更新的，那么不考虑什么非正交循环或特殊的算法，只是最简单的不断根据旧值求新值，solver是不是可以写成下面这种？

    fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );
        while (......) {
        solve(UEqn == .....)
        }

0_1463646525604_捕1111111111获.JPG

但是在代码里写的是

            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                (fvc::interpolate(HbyA) & mesh.Sf())
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );

            adjustPhi(phiHbyA, U, p);

为什么要加上 fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) 这一项？以及在后面adjustPhi(phiHbyA, U, p) ？因为只按照图片中的方程显示，只要 fvc::interpolate(HbyA) & mesh.Sf()) 就够了。

CFD中文网

修正代码显示问题。

@李东岳

:sunglasses:

sejabs

谢谢了~~~，有没有发帖格式的说明？

CFD中文网

@sejabs

有关adjustphi()看这里：adjustPhi的作用是检查边界条件？

代码高亮在回帖的时候，文本框右上方有个很小的“编写？”点击可以看规则。

sejabs

看了给的帖子，以及 https://openfoamwiki.net/index.php/IcoFoam 的说明，虽然细节还不懂，但是有了粗糙的理解。
由于phiHbyA 对应 0_1463648362011_upload-bdb99637-9e62-4dc4-b290-510f5e3d1f1a 在face上的插值，对内部的面而言无所谓，但是在边界上就需要调整，保持一致性。那么phi 也是在面上的插值，由于它对应是 U 的插值，U 已经定义了边界条件，所以就不用调整了。或者 adjustPhi(phi, U, p) 也是可以写的，但是只是在边界上重复一遍没什么用吧。
不知道这样的理解对不对。

CFD中文网

@sejabs

不是很理解你的意思，
adjustphi()是用于对压力全部为Neumann边界条件的情况下，给定一个附加条件使得压力有解，这个附加条件通过adjustphi()进行。

// adjusts the inlet and outlet fluxes to obey continuity, which is necessary for creating a well-posed
// problem where a solution for pressure exists.

这句话也就是要使得压力有解，有必要进行adjustphi()

李东岳

这里的UEqn.A和后面的UEqn.H是自动根据U的新值进行刷新的吗？

UEqn.A()对应的为$A_p$在当前时内迭代中其为不变的，在进入下一个时间步更新；
UEqn.H()对应的是icoFoam解析公式(15)右侧括号内的项，其中$A_n$为不变的但是UEqn.H()是变化的；

 fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );
        while (......) {
        solve(UEqn == .....)
        }

这段代码在UEqn后跟随的solve()，就是把动量方程求解。但是并没有用到连续性方程。还是得需要求解压力修正方程和压力泊松方程。上述代码只是动量预测用。可有可无。

P.S. 关于ddtcorr，需要码大量公式，择日补充。

sejabs

好的，谢谢两位的指点！

李东岳

有关ddtphicorr，简单来说是一种使得计算更精确的方法。在下面这篇文献中对增加ddtphicorr和省略ddtphicorr做了比较。在误差放大100倍的时候才可以对比。

看这篇论文应该会有所启发：NOTE ON THE USE OF MOMENTUM INTERPOLATION METHOD FOR UNSTEADY FLOWS

在ville的论文里直接对ddtphicorr做了对比。误差均非常小。

cfd-online一个讨论帖：http://www.cfd-online.com/Forums/openfoam-solving/60096-ddtphicorr.html
然而cfd-online的讨论大部分还是很具有挑战性。对比文献后我比较赞成eugene的看法。初步看了下文献，和代码还没具体的对上，但是确实是这种技术。以后有时间把公式和方程对应一下。

李东岳

最近看到了一些资料，目前可以确定的是ddtcorr()的引入是由于：“非稳态问题的时间步长较小的时候，采用Rhie-chow插值仍然会导致压力波”。因此“求解和时间相关（尤其在小时间步长的时候），建议采用Choi提出的格式”。且“Rhie-chow插值求解的结果和松弛因子有关”

参考《数值传热学》250页以及《流动与传热数值计算—若干问题的研究和探讨》，我会跟进Choi格式并和ddtcorr()对应一下。

Note on the use of momentum interpolation method for unsteady flows, SK Choi - Numerical Heat Transfer: Part A: Applications, 1999 - Taylor & Francis

dyj19901127

@李东岳按照《流动与传热数值计算—若干问题的研究和探讨》第2.8节和Discussion on momentum interpolation method for collocated grids of incompressible flow. Numerical Heat Transfer: Part B: Fundamentals，choi格式与亚松弛因子无关，仍然与时间步长相关。

李东岳

@dyj19901127

参考choi的论文，其引用Majumdar的文章表明"the resulting converged solution is relaxation factor dependent". 具体原理，还需研究。

我刚看了Yu bo的Discussion on momentum interpolation method for collocated grids of incompressible flow. Numerical Heat Transfer: Part B: Fundamentals，他也是说相关的。他们引用的Majumdar和Miller的文章。确认下？

Anyway，多研究研究，看看choi原理是怎么实现的。Choi那个论文很短。Yo bo师从陶院士？不一般啊。

队长别开枪

@李东岳
我觉得关键点在于进行当前时间步Rhie-chow动量插值时上一时间步界面速度$\mathbf{U}_{f}^{n}$的计算，如果没有ddtcorr()，$\mathbf{U}_{f}^{n}$只是简单通过相邻单元基于动量方程系数矩阵主对角元素加权平均得到，如果加了修正项，那么$\mathbf{U}_{f}^{n}$使用的就是上一时间步里基于动量插值得到的结果，所以在phiHbyA的计算里需要修正界面的体积通量，下面的代码是EulerDdtScheme.C527~545行（OpenFOAM 2.4.0）中的一部分

fluxFieldType phiCorr
(
      phi.oldTime() - (mesh().Sf() & fvc::interpolate(U.oldTime()))
);

- (mesh().Sf() & fvc::interpolate(U.oldTime()))用于消去计算HbyA时加上的体中心处的旧值。修正的原因是即便体中心处的速度场满足连续方程，基于动量方程系数矩阵主对角元素加权平均往面中心插值，时间步长越小，越倾向于简单算术平均，再加上非结构不规则网格等因素，很难再保证插值后的面速度场仍然保持0散度条件。
个人浅见，欢迎批评指正。

yhdthu

@李东岳前辈好，我有个问题，这个UEqn.H()分为两个部分，一个是AnUn，另一个是源项S，在同一个时间步内，随着修正次数变化的只是AnUn，但是S是不变的对不？

李东岳

S不一定，如果S和求解的变量有关，那肯定变啊。

yhdthu

@李东岳前辈，关于动量方程中的扩散项，由于其处理方式为半隐半显，即：

tmp<fvVectorMatrix> laminar::divDevRhoReff(volVectorField& U) const 
{
  return
        (
          - fvm::laplacian(muEff(), U)
          - fvc::div(muEff()*dev2(T(fvc::grad(U))))
        );
}

我目前认为这个显式部分是进入了.H()算子中，我想问的就是这部分是不是在压力修正时（即同一个时间步内）的计算不变的？

李东岳

我目前认为这个显式部分是进入了.H()算子中

这个目前我不确定，你可以看一看。:expressionless:

我是河滩

@东岳对流项是不是对UEqn.A()没有贡献？因为单元面上通量phi求和为0，满足连续性方程。

我是河滩

@CFD中文网压力全为零梯度边界，压力泊松方程离散后的系数矩阵奇异，adjustPhi修改了矩阵，让矩阵不再奇异？

李东岳

@我是河滩在 icoFoam的一些细节问题中说：

对流项是不是对UEqn.A()没有贡献？

取决于你的离散格式，大部分格式有贡献

我是河滩

@东岳是的，还和网格有关，均匀网格和非均匀网格产生的结果也不同。