VOF添加斥力模型

李东岳

还有一点问题是还是会出现添加斥力比较大的时候，因为总会有气泡不断接近到比较近的时候，从而产生一个比较大的斥力

是的，在最开始最开始你在讨论这个问题的时候，我就想到了可能会有这个问题。所以就跟你讨论了一下。因为我最近从OpenFOAM学习到一个方法来处理这个问题。

你这个添加斥力的模型，从另外一种角度，可以近似的看做添加一种分散力。同理的一个现象，在气固流动那面，存在一个相压力模型，在相分数超过设定值的时候，添加一个分散力，将alpha打散。正常的算法就是在压力方程中，比如你得处理方式，把这个分散力加进去。但是OpenFOAM那面把这个问题在alpha方程中也处理了。我测试过，如果仅仅是压力方程处理，可能会导致发散（alpha越界）。如果压力方程+相方程同时处理，就非常稳定。

我觉得你可以借鉴这个方法。

http://dyfluid.com/reactingTwoPhaseEulerFoam.html

PS 你这个斥力是$1/h^3$，那面的力是exp指数级别的函数（最终是一个非常陡暴增，力会突然变得高几十个数量级）

学流体的小明

谢谢李老师提供的思路。我的理解是有两种方法：

处理相方程，粗略看了一下，感觉这个方法挺复杂。
可以改一改斥力的表达式，让它不那么大。

我在输出两个volVectorField，用来比较斥力和表面张力的大小，具体编程上也想请教一下，控制体体心的alpha的单位梯度怎么写。我目前是这样的，单位梯度必须加一个1e-300防止发散。

// Cell gradient of alpha
const volVectorField gradAlpha(fvc::grad(alpha1));
repulsion_Compare = repulsion_Sum*gradAlpha/(mag(gradAlpha)+dimensionedScalar(dimensionSet(0,-1,0,0,0,0,0), 1e-300));
surfaceTensionMy=mixture.sigmaK()*gradAlpha;

OpenFOAM在算表面张力的时候，算曲率 $\kappa$ 的时候用到了面心的alpha单位梯度，也分享一下代码。

const surfaceVectorField& Sf = mesh.Sf();

    // Cell gradient of alpha
    const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));

    // Interpolated face-gradient of alpha
    surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));


    // Face unit interface normal
    surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
  
    correctContactAngle(nHatfv.boundaryFieldRef(), gradAlphaf.boundaryField());

    // Face unit interface normal flux
    nHatf_ = nHatfv & Sf;

   // Simple expression for curvature
    K_ = -fvc::div(nHatf_);

李东岳

可以改一改斥力的表达式，让它不那么大

你直接调节那个K的值不就好了? 反正也是个经验参数。

学流体的小明

K调太小的话，气泡间斥力无法弹开气泡，两个气泡会融合。
现在还挺混乱下面这张图，左侧轴是压力测点的压力，右侧轴是volScalarField repulsion_Sum的最大值。

有的时候没有斥力，也有全局的压力下降，比如下面的椭圆圈出来的。
有的时候斥力生效，但是压力没问题。方框中基本水平的蓝色和橙色线。
总体来说斥力生效会极大地增加压力出问题的风险。

红色箭头是斥力，白色箭头是表面张力，斥力还是比表面张力大一些。

李东岳

才想起来，你这个还是周期性网格，所以需要给定压力参考点，我自己简单算了个纯粹的气泡上升，压力还没啥问题。

学流体的小明

之前试算了很多，发现压力参考点不是最根本的原因，最根本的原因还是添加了这个斥力源项导致的流场压力不稳定。不用压力参考点也算出来了这种压力振荡。
应该是和具体的代码有关系，和SIMPLE算法里面预测方程有关系。
我没有往UEqn中加这个斥力项，仅加在了phig里面，也就是在压力泊松方程里面加了，后面速度修正那一行代码中，斥力也是体现在phig当中。
李老师您怎么写的代码？

李东岳

我没有写全，就在function里面测试了一下这个力的大小（很多参数都是拍脑袋硬植入）：


functions
{
    setDeltaT
    {
        type coded;
        libs            ("libutilityFunctionObjects.so");
        writeControl    writeTime;
        executeControl  timeStep;

        code
        #{
        #};

        codeExecute
        #{
            const volScalarField& alpha
            (
                mesh().lookupObject<volScalarField>("alpha.air")
            ); 

            const surfaceVectorField& Sf = mesh().Sf();

            const volVectorField gradAlpha("gradAlpha", fvc::grad(alpha));

            surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));

            const dimensionedScalar deltaN_
            (
                "deltaN",
                1e-8/pow(average(mesh().V()), 1.0/3.0)
            );

            surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
            volVectorField nHatv("nHatv", gradAlpha/(mag(gradAlpha) + deltaN_));

            surfaceScalarField nHatf = nHatfv & Sf;

            volScalarField K("K", -fvc::div(nHatf));
            volScalarField interface("interface", alpha);
            volVectorField force("force", nHatv);

            List<label> markedCell;

            forAll(mesh().C(), i)
            {
                force[i] = Zero;
                if (alpha[i] < 0.99 && alpha[i] > 0.01)
                {
                    interface[i] = 1.0;
                    markedCell.append(i);
                }
                else
                {   
                    interface[i] = 0.0;
                }
            }

            forAll(markedCell, i)
            {
                label celli = markedCell[i];
                forAll(markedCell, j)
                {
                    label cellj = markedCell[j];
                    scalar dij = mag(mesh().C()[celli] - mesh().C()[cellj]);
                    scalar posneg = nHatv[celli] & nHatv[cellj];

                    if (posneg < 0.0)
                    {
                        if (dij < 0.15)
                        {
                            scalar K = 1e-2;
                            force[celli] += K*nHatv[celli]/pow3(dij);
                        }
                    } 
                }
            }

            markedCell.clear();
            K.write();
            nHatv.write();
            gradAlpha.write();
            force.write();
            interface.write();
        #};
    }
}

学流体的小明

李老师，想请教一下源项的添加的问题。
比如我现在在有一个repulsion源项，要想植入，应该是有三个地方要添加。

一，速度预测方程，即UEqn.H中的fvMatrix UEqn。但我不清楚这里是否要加上源项。
加的理由是：fvOptions这种人工添加的源项都在。
不加的理由是：写的那个icoFoam解析，以及这里实际上只保留了对流项和粘度项，压力、表面张力和重力都没有加，这个斥力在形式上是和表面张力有关的，就像表面张力那样处理。

fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
      + MRF.DDt(rho, U)
      + turbulence->divDevRhoReff(rho, U)
     ==
        fvOptions(rho, U)
      + repulsion_Vector // add
    );

二，真实速度方程，一般的想法都是加上这个源项，就是最后的速度修正项是

但是在reactingTwoPhaseEulerFoam中，我看到这些源项都是直接加到最后的速度上的，不知道这里是您省略它们的具体形式了，还是这些$M$确实就不除以$ A_{d,p} $

就是这两种方式应该是哪一种？还是两种都可以？

三，上面两种不同方式的速度方程，对应不同的压力泊松方程，也就差个$ A_{d,p} $

李东岳

速度预测方程，即UEqn.H中的fvMatrix UEqn。但我不清楚这里是否要加上源项。

不需要加

源项你直接1）加到phi上面组建压力方程，2）最后的速度更新要加上源项

我那个M都忘记除掉A了，是要除掉，我有遗漏，我周末更新过去
http://dyfluid.com/reactingTwoPhaseEulerFoam.html 这个里面41-46之前写的一直没更新 29之前是我最近更新的没问题

学流体的小明

不需要加

源项你直接1）加到phi上面组建压力方程，2）最后的速度更新要加上源项

这种方法算得的结果就是那些有压力振荡不稳定的。

另一种方法：仅添加到UEqn当中，压力方程不做处理，斥力的贡献依靠UEqn.H参与到后续的计算当中。结果就没有压力振荡了。
目前初步的结果是这样的，后续可能还有其它问题。

非常感谢李老师这么久这么多的指点，帮助太大了。

李东岳

另一种方法：仅添加到UEqn当中，压力方程不做处理，斥力的贡献依靠UEqn.H参与到后续的计算当中。结果就没有压力振荡了。

Interesting！

李东岳

@学流体的小明大佬后来你们多相流直接模拟的工作发出来没

学流体的小明

@李东岳谢谢李老师关心呀。算例模拟完了，文章😂我还在写。
算出来的结果不是很好，气泡槽道流收敛性还比较差，几个参数之间的规律不是很明显，故事不太好讲。

李东岳

不好讲故事是比较困境了，哎..

CFD中文网

VOF添加斥力模型