关于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。处理格式比较费时。

zhanghan

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

金石为开

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

Haining LUO

@金石为开我觉得这里是因为网格动到了新的位置，而当前的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^*$项。必须使其满足以上积分的相容性条件。

Haining LUO

@程迪多谢分享，不过上面的帖子说的不是correctPhi.H吗？
你这里提到的是adjustPhi(phiHbyA, U, p)，我看icoFoam里面是在求解Poisson方程前，在pimpleoam和pimpleDyMFoam都是在pEqn.H里面，你提到的pRef也是，所以是为了求解Poisson方程作的准备。我觉得我们说的不是一个东西，这里是ajustPhi的帖子adjustPhi的作用是检查边界条件？

程迪

@Haining-LUO
correctPhi调用的就是adjustPhi。不过adjustPhi通常针对p，correctPhi是针对的pcorr。道理都一样的。保证微分方程相容性。

Haining LUO

@程迪明白你的意思了，在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值是一样，分布不同而已。

Haining LUO

@程迪

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

李东岳

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

这个是OpenFOAM的代码风格。

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

我是河滩

@程迪压力都是零梯度边界条件，压力泊松方程的系数矩阵是奇异的，adjustphi让矩阵不奇异了吗？奇异矩阵怎么解的？