关于correctPhi.H这个函数
-
@金石为开
fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divUdivU在没初始化的时候就是0。即GeometricZeroField还有我认为这是求解之前强制的做了一个压力求解和通量修正，我不清楚这是出于什么的数值上的考虑
我觉得是数值上的。因为correctPhi有的时候可以关闭。并不是必须要开。不过从教材上来看，物理上也必须符合空间守恒。所以。我也不是很明确。
另外，pimpleDyMFoam有个方程写错了，已经修改，参考Peric书方程12.13。不过代码还对不上，从方程12.13来看代码应该为：
fvm::laplacian(rAU, pcorr) == fvc::div(mesh.phi()) - divU感兴趣你可以找个简单测试一下？要把correctPhi设置为on。
-
理论上是为了满足压力方程的相容性条件，实际上调用了adjustPhi.H。
参考
adjustPhi和pRef针对的是压力方程为全Neuman边界的情况，这会有两个问题:
- 方程系数矩阵M为奇异，解有无穷多，相差一个常数，所以需要pRef
- 方程相容性问题，需要adjustPhi
pRef
当方程为：
$$
\Delta p = f \text{ on }\Omega
$$
$$
\nabla p \cdot \mathbf n = g \text{ on }\partial \Omega
$$
首先，如果$p_1$是一个解，那么$p_2=p_1+const$也是方程的解，方程解不唯一，所以离散方程会有奇异性，无法求解。
解决方法是固定一个点的压力为给定值，就是设置pRef。adjustPhi
p变量前面是直接加了微分的，所以需要对于全Neumann条件的情形还需要满足相容性条件：
$$
\int_{\Omega}{f dV} =\int_{\Omega}{\Delta p dV} = \int_{ \partial\Omega}{\nabla p \cdot\mathbf n dS}= \int_{ \partial\Omega}{g dS}
$$这就是AdjustPhi去强制满足这个条件。
带入f和g的表达式$$
f= \nabla\cdot\mathbf u^*
$$
$$
g=0
$$可得：
$$
\int_{\Omega}{\nabla \cdot \mathbf u^* dV} =\int_{\partial \Omega}{\mathbf u^*dS}= 0
$$你再看看adjustPhi出现的位置，是在解压力方程之前，对phiHbyA进行的调整。其实也就是这里的$\mathbf u^*$项。必须使其满足以上积分的相容性条件。
-
@程迪 多谢分享，不过上面的帖子说的不是correctPhi.H吗？
你这里提到的是adjustPhi(phiHbyA, U, p)，我看icoFoam里面是在求解Poisson方程前，在pimpleoam和pimpleDyMFoam都是在pEqn.H里面，你提到的pRef也是，所以是为了求解Poisson方程作的准备。我觉得我们说的不是一个东西，这里是ajustPhi的帖子adjustPhi的作用是检查边界条件？ -
@Haining-LUO
correctPhi调用的就是adjustPhi。不过adjustPhi通常针对p，correctPhi是针对的pcorr。道理都一样的。保证微分方程相容性。 -
@程迪 明白你的意思了，在pimpleFoam和icoFoam里面都是在解p方程前调用了adjustPhi以保证相容性（通过修正流量保证连续方程的方式：其中涉及了两种边界条件fixedValue和inletOutletFvPatchField，修改的是inletOutletFvPatchField上的phi）。
但仔细看pimpleDyMFoam里面的correctPhi.H
if (mesh.changing()) // 如果是动网格，肯定会执行的部分 { forAll(U.boundaryField(), patchI) { if (U.boundaryField()[patchI].fixesValue()) { U.boundaryField()[patchI].initEvaluate(); } } forAll(U.boundaryField(), patchI) { if (U.boundaryField()[patchI].fixesValue()) { U.boundaryField()[patchI].evaluate(); phi.boundaryField()[patchI] = U.boundaryField()[patchI] & mesh.Sf().boundaryField()[patchI]; } } } { // 这里我想说妈蛋……这个大括号是什么鬼 volScalarField pcorr ( IOobject ( "pcorr", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensionedScalar("pcorr", p.dimensions(), 0.0), pcorrTypes ); dimensionedScalar rAUf("rAUf", dimTime, 1.0); // 构建了一个单位且均匀的有时间量纲的场 while (pimple.correctNonOrthogonal()) // 这里其实还没有进入真正的pimpleFoam循环（while (pimle.loop())），但corrNonOrtho_已经存在。但是这里没有adjustPhi !! 为什么 ?? 这里虽然求解的不是物理场p，但是pcorr也得有个边界条件吧? (于是我在createPcorrTypes.H里面找到pcorrTypes几乎是zeroGradient，如果p设定的条件中有fixedValue，对应的pcorr也是fixedValue，所以如果没有fixedValue的p边界这里也是需要adjustPhi的!) { fvScalarMatrix pcorrEqn ( fvm::laplacian(rAUf, pcorr) == fvc::div(phi) ); pcorrEqn.setReference(pRefCell, pRefValue); pcorrEqn.solve(); if (pimple.finalNonOrthogonalIter()) { phi -= pcorrEqn.flux(); } } #include "continuityErrs.H" }在这一段代码之后才是pimple.loop()，所以我觉得这是因为
mesh.update(); // Calculate absolute flux from the mapped surface velocity phi = mesh.Sf() & Uf;用新网格算了新phi，而新的phi不是无散度的，因此用以上来修正之后再进入N-S方程求解。不过@李东岳 指出的代码和我这里OpenFOAM/2.3.1-foss-2016a不大一样，我觉得这里求解的并不是Peric书方程12.13，你觉得呢？
-
Peric 12.13是空间守恒律吧。网格的东西，和这玩意儿有啥关系？
据我所知，volScalarField存的是体心的物理量，而不是体心物理量×体积，网格变形相当于物理量分布变化了，新phi和老phi值是一样，分布不同而已。
