关于piso循环的二次修正的一个疑问
-
比如说，在icoFoam的piso循环内
for (int corr=0; corr<nCorr; corr++) { volScalarField rAU(1.0/UEqn.A()); volVectorField HbyA("HbyA", U); HbyA = rAU*UEqn.H();//公式(15) surfaceScalarField phiHbyA ( "phiHbyA", (fvc::interpolate(HbyA) & mesh.Sf())//此处依据Rhie-chow插值原理，HbyA使用线性插值得到， 即需要在算例中设定interpolate(HbyA)的格式 + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi) );//phiHbyA即为以HbyA来计算的通量，第一次内循环的时候，此处phiHbyA并不守恒，not divergence free adjustPhi(phiHbyA, U, p);//此函数功能请参考：更新4 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::laplacian(rAU, p) == fvc::div(phiHbyA)//公式(26)，有关div(phiHbyA)请参考：更新3 ); pEqn.setReference(pRefCell, pRefValue); pEqn.solve();//公式(26) if (nonOrth == nNonOrthCorr) { phi = phiHbyA - pEqn.flux();//使用求解的压力场修正通量场，在最后一次修正的时候通量守恒，Issa指出，大约需要2-3次内循环步。 对应公式(25)，pEqn.flux()返回公式(25)方程右边第二项，也为fvc::interpolate(rUA)*fvc::snGrad(p)*mag(mesh.Sf())。 某些可压缩求解器其中的pEqn.flux()可能为+号，即为phi = phiHbyA + pEqn.flux()。这是因为pEqn中的laplacian项为−−号 } } #include "continuityErrs.H" U = HbyA - rAU*fvc::grad(p);//公式(17)，这里并没有采用重组守恒通量phi的方法来计算速度，即U=fvc::reconstruct(phi)； Henry表示采用fvc::reconstruct(phi)会引起比较大的守恒误差。目前这种误差有多大是未知的。更多相关的讨论请参考：更新1，更新2 U.correctBoundaryConditions(); }piso循环的二次修正是建立在第一次修正的基础上的。第二次修正的连续性方程是涉及到u**。上面的代码里，UEqn并没有做更新。如果忽略函数fvc::ddtCorr(U, phi)，pEqn也是没有变动的。
问题1： 这是不是就意味着对同一个pEqn重复一次迭代运算？问题2：我对问题1的理解是，在icoFoam的算法里，速度u通过连续性方程被隐性的耦合到了pEqn里面。第二次的修正是没有必要的。这与原始的PISO算法有区别。请问这样的理解对吗？
问题3：另外，在sonicLiquidFoam里面的pimple循环里：
while (pimple.correct()) { volScalarField rAU("rAU", 1.0/UEqn.A()); surfaceScalarField rhorAUf ( "rhorAUf", fvc::interpolate(rho*rAU) ); U = rAU*UEqn.H(); surfaceScalarField phid ( "phid", psi *( (fvc::interpolate(U) & mesh.Sf()) + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho) ) );// phi = (rhoO/psi)*phid;//important fvScalarMatrix pEqn//for solution of pEqn, construction with respect to p0 and rho0 ( fvm::ddt(psi, p)//=0 if psi=0 + fvc::div(phi)//// + fvm::div(phid, p)//=0, if psi=0 - fvm::laplacian(rhorAUf, p)// );// pEqn.solve(); phi += pEqn.flux();// solve(fvm::ddt(rho) + fvc::div(phi));//solve for rho #include "compressibleContinuityErrs.H" U -= rAU*fvc::grad(p); U.correctBoundaryConditions(); }pEqn里面的phid的构建是与U相关的。与icoFoam对比，此情况下的二次修正就显得有意义。请问，这种想法对吗？
-
-
piso循环的二次修正是建立在第一次修正的基础上的。第二次修正的连续性方程是涉及到u**。上面的代码里，UEqn并没有做更新。如果忽略函数fvc::ddtCorr(U, phi)，pEqn也是没有变动的。
问题1： 这是不是就意味着对同一个pEqn重复一次迭代运算？虽然构建的速度矩阵
Ueqn没有更新，但是HbyA有变化，因此求解的pEqn不是同一个pEqn问题2：我对问题1的理解是，在icoFoam的算法里，速度u通过连续性方程被隐性的耦合到了pEqn里面。第二次的修正是没有必要的。这与原始的PISO算法有区别。请问这样的理解对吗？
不是很清楚你说的哪个二次修正，如果你说的是楼上指的非正交修正，一般情况下不需要进行，除非在那个
potentialFoam中。因为这个求解器没有进行速度-压力耦合迭代。 -
@cfd-china
在非正交修正的循环内即，fvScalarMatrix pEqn ( fvm::laplacian(rAU, p) == fvc::div(phiHbyA)//公式(26)，有关div(phiHbyA)请参考：更新3 ); pEqn.setReference(pRefCell, pRefValue); pEqn.solve();//公式(26)反复计算pEqn，rAU与phiHbyA 都不变，每次更新的p也就不变啊，看不出来非正交修正体现在哪个地方，请问这里怎么理解？
-
Jasah建议一般情况下不进行压力非矩形修正，因为本身SIMPLE，PISO算法就可以进行速度压力的修正计算。如果要进行这一部分修正，隐形的每一次新的计算都调用了新的非正交修正。并且是迭代进行的，see:
if ( corr == nCorr-1 && nonOrth == nNonOrthCorr ) { pEqn.solve(mesh.solver("pFinal")); } else { pEqn.solve(); }这个代码实际上做的就是potentiaoFoam求解，然而细致的解析需要首先对非矩形修正做分析，目前并没有足够的添加内容。你可以运行
potentialFoam算例看压力非矩形修正的效果。OpenFOAM编程指南里面有详细的对比。
-
@李东岳 在 关于piso循环的二次修正的一个疑问 中说：
Jasah建议一般情况下不进行压力非矩形修正，
http://www.cfd-china.com/topic/840/运行出错/5
好像这个帖子里面进行非正交修正解决问题了？ -
@hongfu2233
非正交修正是这样的：- fvm::laplacian返回的是fvMatrix，fvMatrix既包含矩阵M，又包含源项b。
- fvc::laplacian返回的volScalarField，相当于只有b，没有M。
非正交修正是把laplacian算子分成正交部分和非正交部分
- 对于fvc而言，反正两部分最后要加在一起，无所谓啦，都一样；
- 而fvm只能对正交部分进行隐式离散，相关的系数进入矩阵M（并非不能对非正交部分进行隐式离散，而是这样的话M矩阵就不是和mesh对应的lduMatrix相同的稀疏结构了），非正交部分显式离散，进入源项b中。
- 这种伎俩叫延迟修正(deferred correction)，我看来其实就是部分隐式部分显式的离散。
部分显式部分隐式的离散会带来信息滞后的问题，比如说$\Delta p^{(n)} = s$方程整成了$L_I(p^{(n)})=s - L_E(p^{(n-1)})$，虽然$\Delta p^{(n)} = L_I(p^{(n)})+L_E(p^{(n)})$，但是$\Delta p^{(n)} \ne L_I(p^{(n)})+L_E(p^{(n-1)})$，除非$p^{(n)}=p^{(n-1)}$，也就是$p$已经收敛
所以说对于fvc而言，这没有任何问题，但对于fvm而言，你一次求解只能求得：$p^{(n)}=L_I^{-1}(s - L_E(p^{(n-1)})) \ne \Delta^{-1}(s)$，这个误差只有靠不断地求解同样的方程（比如pEqn），也就是非正交修正(nonOrthogonalCorrection)中看似没有变化的那个pEqn来消除。在非正交修正的迭代过程中，其实每次解的方程不是不变的，只是其中的fvm::laplacian生成的隐式项系数没有变化，显式项是和你给的$p$有关的，是有变化的，所以实际上实现的是上述的$p^{(n-1)}$到$p^{(n)}$的迭代修正循环。
其实OpenFOAM中很多的延迟修正都是类似的部分隐式部分显式离散，比如High Resolution格式和High Order格式的差值也是这么搞的。
这种半隐式离散理论上必须达到最终收敛才是和全隐式离散是一样的解，但是只要显式部分远小于隐式部分，问题也不大，甚至收敛阶数和速度也不太受影响。不过对于压力方程还是解得准确点儿好。
