关于correctPhi.H这个函数



  • 对于多相流VOF方法要先运行一下correctPhi.H这个文件,我见功能是修改通量,Required during start-up, restart, mesh-motion etc.在开始的地方运行,这个怎么理解?采用什么公式方法,请教大家~~~



  • correctPhi.H

    Flux correction functions to ensure continuity.
    
    Required during start-up, restart, mesh-motion etc. when non-conservative fluxes may adversely affect the prediction-part of the solution algorithm (the part before the first pressure solution which would ensure continuity). This is particularly important for VoF and other multi-phase solver in which non-conservative fluxes cause unboundedness of the phase-fraction.
    

    interFoam有一个动网格求解器叫interDyMFoam,两个求解器可以互相调用。从描述来看:correctPhi强制连续性方程。因此通过下述代码:

    while (pimple.correctNonOrthogonal())
        {
            // Solve for pcorr such that the divergence of the corrected flux
            // matches the divU provided (from previous iteration, time-step...)
            fvScalarMatrix pcorrEqn
            (
                fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
            );
    
            pcorrEqn.setReference(0, 0);
            pcorrEqn.solve();
    
            if (pimple.finalNonOrthogonalIter())
            {
                phi -= pcorrEqn.flux();
            }
        }
    

    求解由连续性方程构造的泊松方程后,重组守恒的phi场:phi -= pcorrEqn.flux();

    这是一个数值上的考虑,可有可无。



  • @李东岳 东岳老师,我看过你写的pimpleDyFoam里面的讲解,对于动网格,变动之后有体积守恒定律,网格动了之后要进行通量修正。代码中的divU其实就是0场(公式里面对1做时间导数为0),不知道理解的对不对?
    还有我认为这是求解之前强制的做了一个压力求解和通量修正,我不清楚这是出于什么的数值上的考虑或者是物理上的考虑?发散?稳定?精度还是其他谢谢老师。另外您的fluent理论指南好了没有呀?很期待看看有什么学习参考的,谢谢老师



  • @金石为开 fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU divU在没初始化的时候就是0。即GeometricZeroField

    还有我认为这是求解之前强制的做了一个压力求解和通量修正,我不清楚这是出于什么的数值上的考虑

    我觉得是数值上的。因为correctPhi有的时候可以关闭。并不是必须要开。不过从教材上来看,物理上也必须符合空间守恒。所以。我也不是很明确。

    另外,pimpleDyMFoam有个方程写错了,已经修改,参考Peric书方程12.13。不过代码还对不上,从方程12.13来看代码应该为:

    fvm::laplacian(rAU, pcorr) == fvc::div(mesh.phi()) - divU
    

    感兴趣你可以找个简单测试一下?要把correctPhi设置为on。



  • 哦对了还有fluent理论指南。最近琐事太多了:sad: 。都已经翻译完了。就是格式太难统一了。因为是各自翻译的。用的是word。处理格式比较费时。



  • @李东岳 李老师,您有研究过把ansys的计算结果导到openfoam中,或者是其他软件中的数据导到openfoam中吗??望您给一些指点!



  • @李东岳 coreectPhi.C出错了??phi变成mesh.phi()??如果这里错了,那,,所有用到correctPhi的地方都错了呀?。。是不是呢?。。。



  • @金石为开 我觉得这里是因为网格动到了新的位置,而当前的phi相对于当前的网格就不再divergence free,网格更新了而phi场也需要更新,在当前网格下对absolute phi进行修正,而不是对应的方程12.13。假设了一个pcorrEqn场的拉普拉斯等于当前phi的不为零散度,然后用phi -= pcorrEqn.flux();修正。



  • 理论上是为了满足压力方程的相容性条件,实际上调用了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,你觉得呢?



  • @Haining-LUO

    Peric 12.13是空间守恒律吧。网格的东西,和这玩意儿有啥关系?

    据我所知,volScalarField存的是体心的物理量,而不是体心物理量×体积,网格变形相当于物理量分布变化了,新phi和老phi值是一样,分布不同而已。



  • @程迪

    phi是surfaceScalarField,在phi = mesh.Sf() & Uf;之前,phi还是定义在老网格的位置,通过Uf=fvc::interpolate(U)得到Uf后(我觉得是相对新网格的插值),这样就定义了一个新phi:phi = mesh.Sf() & Uf;,于是新的phi定义在新的网格面上,不再满足连续方程(旧phi是满足的),所以在进入pimple循环前通过pcorr来修正。



  • { // 这里我想说妈蛋……这个大括号是什么鬼

    这个是OpenFOAM的代码风格。

    其他内容,我需要周末研究下才能跟你们讨论 :cheeky:


登录后回复
 

与 CFD 中国 的连接断开,我们正在尝试重连,请耐心等待