interFoam计算气泡槽道流时的压力问题

学流体的小明

各位老师同学好，我现在在算一个槽道，里面有很多气泡。基于interFoam自己添加了一些需要的功能，所以水和气都是不可压的。现在的问题是流场的压力结果总是不好。
算例是这样的：

x方向为流动方向，周期边界条件；y方向是壁面，上下都是壁面；z方向是展向，也是周期边界条件。流场通过meanVelocityForce驱动，流场的平均速度会基本保持在一个量上。没有重力。
之前先使用pimpleFoam算过一个单相的湍流槽道流，达到稳定后，把湍流槽道流的速度场mapFields到现在这个算例当中，然后用setFields把气泡初始化，具体来说是alpha.water、p_rgh、p都初始化。
使用了自适应网格加密，更好地捕捉界面，运行起来网格大约在7e6 ~ 8e6。
由于压力边界条件没有在任何一个边界上给定一个压力值，所以interFoam求解器（其它pimple算法求解器同样也）会启用“压力参考点”这一套东西。就是指定一个pRefCell，给一个pRefValue，在每一次压力计算结束之后，流场会整体地加一个数值，使得这个pRefCell的压力总等于pRefValue。具体代码如下：

p == p_rgh + rho*gh;

    if (p_rgh.needReference())
    {
        p += dimensionedScalar
        (
            "p",
            p.dimensions(),
            pRefValue - getRefCellValue(p, pRefCell)
        );
        p_rgh = p - rho*gh;
    }

当槽道中仅有一个气泡时，每次修正的值并不大，还可以接受，流场中某点的压力随时间的变化（由于没有重力p=p_rgh）如下图，是可以接受的，由于有气泡经过，该点p会有几次峰。

但是在我计算下面这个很多气泡的流场时

流场中一些监测点的压力会变成下面这样，可以看到很多点的压力基本是一样的，而且振荡很厉害，完全不连续。这是因为每次修正的值，也就是上面公式中的pRefValue - getRefCellValue(p, pRefCell)很大，淹没了该点本身有意义的值。

从上面的压力场中我无法进行进一步的分析，平均值和脉动平方均值都感觉没什么意义——因为修正操作造成的波动太大了。

自己分析有以下几种解决方案，也尝试了一下

是因为压力计算没有收敛吗？使用合适的fvSolution和fvScheme可以解决？fvSolution和fvScheme我放在下一层。
关掉压力修正操作，流场的压力可能会一直增大或者减小，最后可能是1e7Pa这样很大的值，但是压力梯度本身还是有用的。
关键原因是求解方程，求解出来pRefCell的值可能就是500，然后流场就会整体减去500让它等于0。

while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

        if (pimple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - p_rghEqn.flux();

            p_rgh.relax();

            U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
            U.correctBoundaryConditions();
            fvOptions.correct(U);
        }
    }

各位大佬有遇到过相似的问题吗？求指点🙏🙏🙏

学流体的小明

fvSolution

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    "alpha.water.*"
    {
        nAlphaCorr      2;
        nAlphaSubCycles 3;
        cAlpha          1.0;

        MULESCorr       yes;
        nLimiterIter    5;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-07;
        relTol          0.01;
        smoother        GaussSeidel;
    }

    "pcorr.*"
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-8;
        relTol          0;
    }

    ".*(rho|rhoFinal)"
    {
        solver          diagonal;
    }   

    p_rghFinal
    {
        $p_rgh;
        tolerance       1e-07;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-08;
        relTol          0;
        nSweeps         1;
    }

    "(T|k|B|nuTilda).*"
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-8;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor no;
    nOuterCorrectors 1;
    nCorrectors     3;
    nNonOrthogonalCorrectors 0;
    correctPhi     yes;

    pRefPoint      (0.0201 0.0101 0.0101);
    pRefValue      0;
}


// ************************************************************************* //

fvScheme

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss linear;

    div(rhoPhi,U)  Gauss limitedLinearV 1;
    div(phi,p)      Gauss linearUpwind grad(p);//limitedLinear 1;
    div(phi,k)      Gauss linear;
    
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

// ************************************************************************* //

dynamicMeshDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      dynamicMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dynamicFvMesh   dynamicRefineFvMesh;

// How often to refine
refineInterval  5;

// Field to be refinement on
field           alpha.water;

// Refine field inbetween lower..upper
lowerRefineLevel 0.001;
upperRefineLevel 0.999;

// If value < unrefineLevel unrefine
unrefineLevel   10;

// Have slower than 2:1 refinement
nBufferLayers   1;

// Refine cells only up to maxRefinement levels
maxRefinement   1;

// Stop refinement if maxCells reached
maxCells        36000000;

// Flux field and corresponding velocity field. Fluxes on changed
// faces get recalculated by interpolating the velocity. Use 'none'
// on surfaceScalarFields that do not need to be reinterpolated.
correctFluxes
(
    (phi none)
    (nHatf none)
    (rhoPhi none)
    (alphaPhi0.water none)
    (ghf none)
    (alphaPhiUn none)
);

// Write the refinement level as a volScalarField
dumpLevel       true;


// ************************************************************************* //

coolhhh

@学流体的小明我没算过多相流，按照算单向流设置压力参考点有些疑问：

我在模拟建筑绕流时候，用的速度入口。此时出口可以设置fixedValue=0或者zeroGradient。设置fixedValue=0，此时压力参考点没有调用。设置zeroGradient，此时压力参考点对结果影响就很大，需要设置在建筑物位置上方一定距离且要在建筑物影响范围外，此时建筑物位置附近就不受非物理压力脉动影响。
对不可压缩流体流体，速度入口不满足每个瞬时的质量通量平衡，会导致流场发生非物理压力脉动。直观理解就是整个流场每个时刻输入的流量不同，时多时少，就会造成额外的人为挤压或抽空，压力脉动就大。这个时候可以对入口做个简单的质量通量修正，来抑制非物理压力脉动。参考文章：《Kim, Y., I.P. Castro and Z. Xie, Divergence-free turbulence inflow conditions for large-eddy simulations with incompressible flow solvers. Computers & Fluids, 2013. 84: p. 56-68.》
当速度入口满足质量通量平衡，无论出口采用fixedValue=0或者zeroGradient，流场里都不会有压力脉动。当入口质量不平衡，采用fixedValue=0在出口边界压力恒为0，但流场中有压力脉动；采用zeroGradient，在参考点附近的压力接近与0，此处的建筑模拟结果就会合理，但流场其他地方仍然是存在压力脉动。
采用mapFields,可能会出现入口质量通量不平衡，但是用的周期边界似乎不存在这个问题。
流场中是否有气泡完全不经过的地方？压力参考点的设置的不受气泡影响的位置处
尝试入口如果只加平均速度剖面，出口设置fixedValue=0，这个时候流场肯定没有压力脉动。再同样测试下气泡槽道中的同样位置的监测点压力，看看是否也是压力脉动过大。如果是的话，那就说明气泡产生的影响不可忽略，就类似于流场中放了个建筑绕流一样，建筑周边就有大的压力波动。
原来的单向流槽道用meanVelocityForce驱动，因为边界没有指定压力值，应该也调用了压力参考点作用？查看算例是指定了pRefCell 1001; pRefValue 0;。整个流场应该是压力梯度驱动流动，但查看of自带的Channel395算例结果的压力均值，好像没有明显看到沿着X向的压力梯度？

学流体的小明

@coolhhh 感谢你的回复和讨论。

关于质量通量平衡这一点，确实就像你说的，我用周期边界应该不受这个影响。

coolhhh 在 interFoam计算气泡槽道流时的压力问题中说：
需要设置在建筑物位置上方一定距离且要在建筑物影响范围外，此时建筑物位置附近就不受非物理压力脉动影响。

参考点的位置确实会极大地影响到它附近的压力和流场周围的压力，比如下面这张图片，参考点附近的pPrime2Mean是很小的。

Fig. t=0.45s，展向中间截面，pPrime2Mean

coolhhh 在 interFoam计算气泡槽道流时的压力问题中说：
流场中是否有气泡完全不经过的地方？压力参考点的设置的不受气泡影响的位置处

有的，从当前算出来的结果来看，气泡不会附着到壁面上，也很少靠近壁面，将参考点设置在壁面附近是一种选择，但这会带来另一个问题
——由于我的算例上下都是壁面，所以参考点设置在靠近壁面的地方会导致两个壁面上的结果不一致。下面图片是单相槽道流的壁面pPrime2Mean结果，参考点设置是pRefCell=0，pRefValue=0，可以看到下壁面的四个角附近的pPrime2Mean是蓝色的很小的值。由于我是要对壁面上的压力进行分析，所以发现这个问题之后，我重新算了参考点设置在计算域中心的的算例。

coolhhh 在 interFoam计算气泡槽道流时的压力问题中说：
尝试入口如果只加平均速度剖面，出口设置fixedValue=0，这个时候流场肯定没有压力脉动。再同样测试下气泡槽道中的同样位置的监测点压力，看看是否也是压力脉动过大。如果是的话，那就说明气泡产生的影响不可忽略，就类似于流场中放了个建筑绕流一样，建筑周边就有大的压力波动。

感觉气泡产生的影响确实是不可忽略的。

coolhhh 在 interFoam计算气泡槽道流时的压力问题中说：
原来的单向流槽道用meanVelocityForce驱动，因为边界没有指定压力值，应该也调用了压力参考点作用？查看算例是指定了pRefCell 1001; pRefValue 0;。整个流场应该是压力梯度驱动流动，但查看of自带的Channel395算例结果的压力均值，好像没有明显看到沿着X向的压力梯度？

我计算的单相槽道流也是基于channel395的各种设置的，是调用了压力参考点，但是单相槽道流中，每次修正的值都非常小，压力数据是可以用的。下面的图是壁面上沿流向的一条线(y=0,z=0.01,x变化)的时间历程，pPrime=p-pMean都是直接从OpenFOAM的后处理结果处理出来的，这张图就很nice。

因为是槽道流，设置的压力梯度是和壁面剪切力平衡了的，确实不会看到沿X向的压力梯度。meanVelocityForce相当于加了源项，和重力的表现不一样。

coolhhh

@学流体的小明那如果流场中都是充满气泡，难以有合适的地方设置压力参考点，感觉就只能用速度入口加湍流的方式模拟。或者能否把单向槽道流某一个面的风速时程提取出来作为气泡槽道的速度入口加载。

学流体的小明

@coolhhh 但这样会导致很多问题，入口速度设置好之后，流场还要向下游发展，如何保证湍流没有被耗散掉——这个时候应该不能用meanVelocityForce了。
我是要得到一个稳定的气泡输运过程，如果按照您的思路，气泡在上游初始化，之后得和湍流一起运动发展一段时间，到达计算域中游下游时候，我可以采集我需要的信息，这会造成我的计算域非常长。而且也很难保证最终是一个遍历态。
还是先想想从求解压力的方面入手吧😂

李东岳

流场中一些监测点的压力会变成下面这样，可以看到很多点的压力基本是一样的

为什么不同的网格点的压力是一样的？

当槽道中仅有一个气泡时，每次修正的值并不大，还可以接受，流场中某点的压力随时间的变化（由于没有重力p=p_rgh）如下图，是可以接受的，由于有气泡经过，该点p会有几次峰。

单气泡算例，p符合预期么（跟表面张力的关系）

学流体的小明

李东岳在 interFoam计算气泡槽道流时的压力问题中说：
为什么不同的网格点的压力是一样的？

A点的值可能是100+0.1，B点的值可能是100+0.2，就是修正的量太大了，把它们本身的空间起伏都淹没了

李东岳在 interFoam计算气泡槽道流时的压力问题中说：
单气泡算例，p符合预期么（跟表面张力的关系）

符合，多气泡算例的表面张力关系也是符合的。下面这个截面的图，p_water = -289 Pa, p_gas = -165 Pa，虽然气泡半径0.001 m，应该是差144 Pa，但考虑到变形，算出来的结果差124 Pa也还是符合预期的。

下面这张图，中间那个红点是压力参考点，p = 0 Pa，两侧的橘色的斑块是气泡。就可以看到在参考点之外，流场的液体的压力迅速就降低到-289 Pa了。这算是没收敛么

coolhhh

@学流体的小明因为这算的是多相流，压力参考点能够指定在water还是gas成分上吗？整个water的压力基本一致，如果压力参考点是以water为基准，那么整个water的压力应该都为0。会不会是默认的压力参考点以gas为准，而gas的压力很小接近于0，然后以这个值为参考点，water的压力值就变为-289？

李东岳

这个看起来是算法的问题。

单核跑是否有问题？
correctPhi如果变成no是否有问题？

学流体的小明

@coolhhh
从pEqn.H的代码来看，并没有这种water / gas之分
@李东岳

没试过单核跑，单核算得会非常慢，而且估计得算几步才能看结果，我试一下。
correctPhi是必须要是yes的，改成no的话会造成自适应加密出问题，流场所有地方的alpha.water会变得不严格是1，比1小一些，然后流场的所有地方都会加密。

学流体的小明

@李东岳东岳老师，最近试了很多参数，发现了点规律。

应该是自适应网格的问题。如果使用自适应网格而不显式地指定correctPhi，程序会自动调用correctPhi。并且自适应网格必须correctPhi，否则会导致alpha.water出问题，全场都加密。

bool correctPhi
(
    pimple.dict().getOrDefault("correctPhi", mesh.dynamic())
);

不使用自适应网格，但是correctPhi yes，流场的压力结果也正常。
单核跑，即使是自适应网格，自动调用correctPhi，流场也计算正常。

应该是并行 + 自适应网格的问题？是不是要在dynamicMeshDict中进行correctFluxes操作呀？但我看Tutorial里面都没correctFluxes。

李东岳

我看到Correctphi代码里面：

    while (pcorrControl.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 (pcorrControl.finalNonOrthogonalIter())
        {
            phi -= pcorrEqn.flux();
        }
    }

第10行是硬植入的，你看看把pcorrEqn.setReference(0, 0)改成p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell))是否有作用。

学流体的小明

我按照您说的方法修改了代码，没有什么大的改变。虽然这行代码确实不太对。
还有就是不太明白setRefrence()这个函数

template<class Type>
void Foam::fvMatrix<Type>::setReference
(
    const label celli,
    const Type& value,
    const bool forceReference
)
{
    if ((forceReference || psi_.needReference()) && celli >= 0)
    {
        source()[celli] += diag()[celli]*value;
        diag()[celli] += diag()[celli];
    }
}

把方程中，RefCell位置的源项加上一个值，然后系数矩阵对角线上RefCell位置翻倍？

还发现一点是自适应加密的间隔可以调小，改成1。每次都运行自适应这个模块，但不是每一步都会导致网格变化，所以调小一些也没事。可能需要nBufferLayers调大到2。
自适应加密间隔调小之后，流场也算得比较好了。自适应加密间隔再从1调整到5，压力马上就出问题了。下面图中0.09s之后就是自适应加密间隔调大的后果。

学流体的小明

学流体的小明在 interFoam计算气泡槽道流时的压力问题中说：

我按照您说的方法修改了代码，没有什么大的改变

甚至导致结果更差了。

不进行这样的修改的话，相同的设置下，相同测点的压力如下图，可以看到紫色线有明显的台阶上升，这是气泡经过了这个测点，测到了气泡内部的压强。

李东岳

把方程中，RefCell位置的源项加上一个值，然后系数矩阵对角线上RefCell位置翻倍？

信息量很大。我稍后回复你这个问题。不过圣诞节要放假了。明天没回复你就得月底了。

这个看起来是个bug。不过要让官方debug测试，他们需要你提供一个小网格的算例他们复现一下。你这个原本的3D-LES他们debug不动。

学流体的小明

我搞了一个125000网格的小算例，也可以复现流场整体异常压力跳跃的问题，不用跑完，跑几分钟就可以用里面的gnuplot脚本画图看到流场的压力。
算例放上来了，李老师@李东岳您有空帮我看看，非常感谢！
singleBubbleChannel_debug.zip

李东岳

把方程中，RefCell位置的源项加上一个值，然后系数矩阵对角线上RefCell位置翻倍？

我在无痛苦ns方程笔记里面详细的写了这个，你看一下setReference()这一节。

我下载看一下能不能跑动。感觉有点够呛，我这是单核虚拟机。

学流体的小明

单核跑算出来的结果非常好，图中所有测点的压强降到-140左右是因为气泡经过了参考点，参考点（此时在气泡内部）的压强为0，外部的水的压强就应该是-144。这个算例的参考点、测点的设置和我上传的那个不一样。

学流体的小明

学流体的小明在 interFoam计算气泡槽道流时的压力问题中说：

应该是并行 + 自适应网格的问题？

更正一下，是并行+自适应网格+压力参考点三个要素的共同作用，去掉任何一个都可以算出来好的结果。

一个新的发现是壁面上会有奇异点，如下图，着色是p_rgh，这些奇异点都处于并行分区的各个界面上。更奇怪的是，上壁面并没有这样的问题。Note:参考点放置在流场中间。

现在暂时放弃自适应加密了，直接全场加密算，压力就是正确的。
😂

CFD中文网

interFoam计算气泡槽道流时的压力问题