关于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对比,此情况下的二次修正就显得有意义。请问,这种想法对吗?



  • 矩阵A跟向量H应该是一样的,但是我记得书上说第二次更正的时候需要用到第一次更正的值(速度跟压力都是)



  • 你说得有道理,压力方程的源在第二次修的时候看起来跟第一次修正时候的源是一样的。

    不清楚这个东西volVectorField HbyA(“HbyA”, U);是不是跟上一次修正最后的更新速度有关系。


  • OpenFOAM教授

    OpenFOAM认为在低库朗数情况下,动量方程中的非线性项 $\nabla \cdot (\textbf u \textbf u)$ 对系统的影响远比速度压力耦合对系统的影响小,因此在循环时不更新速度。



  • 你有比较过带fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)跟不带它的结果么?



  • 这个问题我应该很早就告诉你答案了,以下是$\nabla$-气固大神的回复:

    那个HbyA在第二个piso loop的时候的值不同了,但是A是不变的,所以H的值使用了上一个piso loop最后的U。

    在2~3次修正之后,最后的U才认为是守恒的。


  • OpenFOAM讲师

    同问,这个非正交修正到底是个什么意思


  • 管理员

    @yhdthu
    在离散梯度的时候,面法相梯度只有正交网格的时候计算才是精准的。对于非正交网格,openfoam引入了修正,就是非正交修正。不过这个修正会降低对角占优引起不稳定,因此多试试看,limited一下。

    参考:Jasak these, page 84


  • OpenFOAM讲师

    @cfd-china 谢谢回复,是不是在求解
    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
    的时候,Ap和An中包含了梯度离散项的原因,所以在反复迭代中产生值的变化?


  • 网格教授 OpenFOAM教授 管理员

    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也就不变啊,看不出来非正交修正体现在哪个地方,请问这里怎么理解?


  • 网格教授 OpenFOAM教授 管理员

    @hongfu2233

    Jasah建议一般情况下不进行压力非矩形修正,因为本身SIMPLE,PISO算法就可以进行速度压力的修正计算。如果要进行这一部分修正,隐形的每一次新的计算都调用了新的非正交修正。并且是迭代进行的,see:

    if
    (
      corr == nCorr-1 && nonOrth == nNonOrthCorr
    )
    {
      pEqn.solve(mesh.solver("pFinal"));
    }
    else
    {
      pEqn.solve();
    }
    

    这个代码实际上做的就是potentiaoFoam求解,然而细致的解析需要首先对非矩形修正做分析,目前并没有足够的添加内容。你可以运行potentialFoam算例看压力非矩形修正的效果。

    OpenFOAM编程指南里面有详细的对比。



  • @李东岳 多谢啊!我去看看



  • @李东岳

    @李东岳 在 关于piso循环的二次修正的一个疑问 中说:

    piso循环的二次修正是建立在第一次修正的基础上的。第二次修正的连续性方程是涉及到u**。上面的代码里,UEqn并没有做更新。如果忽略函数fvc::ddtCorr(U, phi),pEqn也是没有变动的。
    问题1: 这是不是就意味着对同一个pEqn重复一次迭代运算?

    虽然构建的速度矩阵Ueqn没有更新,但是HbyA有变化,因此求解的pEqn不是同一个pEqn

    请教一下,HbyA内循环里是通过哪个函数更新的呀?实在没有头绪



  • @李东岳

    重新看了一下icoFoam,恍然大悟。因为通过压力的修正,速度也得到更新。H 操作相当于也更新了HbyA。
    但是貌似这种求解思路与原始的piso算法还是有些区别的。最明显的是,没有了二次修正。我对此的理解是,因为icoFoam里面没有忽略掉式7-39(computational methods for fluid dynamics P175)右边的第二项,
    0_1489643730610_upload-27534b17-accf-4624-8809-3469da030b9e
    所以这个二次修正式可以省略。

    请问这样的理解对吗?



  • @yuan_neu 这个是压力修正方程,和icoFoam完全对不上啊。



  • @李东岳 在 关于piso循环的二次修正的一个疑问 中说:

    Jasah建议一般情况下不进行压力非矩形修正,

    http://www.cfd-china.com/topic/840/运行出错/5
    好像这个帖子里面进行非正交修正解决问题了?



  • 请教一下,大家想要了解piso循环的算法是可以理解的,但是大多数情况下piso循环部分的代码好像并不需要修改,我理解只要写出正确的UEqn就可以了,piso循环对于一般的UEqn都是适用的,是这样吗?什么情况下需要修改piso循环中的代码呢?谢谢!



  • 我觉得UEqn并不是重要的,更重要的是pEqn,有的时候UEqn并没有求解…


  • OpenFOAM讲师

    @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格式的差值也是这么搞的。

    这种半隐式离散理论上必须达到最终收敛才是和全隐式离散是一样的解,但是只要显式部分远小于隐式部分,问题也不大,甚至收敛阶数和速度也不太受影响。不过对于压力方程还是解得准确点儿好。