VOF添加斥力模型

[学流体的小明]

更新一下，顺便就新的问题求助。
我现在全场加密之后算气泡，然后添加了气泡间斥力，是添加在了phig当中。

pEqn.H
surfaceScalarField phig
    (
        (
            mixture.surfaceTensionForce()
          + repulsion // 防止气泡融合的斥力
          - ghf*fvc::snGrad(rho)
        )*rAUf*mesh.magSf()
    );

然后发现当斥力生效的时候，监测的压强出现了我上面帖子中同样的情况，流场中的压力同步地增大或减小，如下图3.5s、4.5s附近的压力曲线。

Fig1. 多个probe的p_rgh随时间的变化

流场如下图，这是一个二维的算例（三维的一样有这个问题，所以拿二维的debug），三个泡在槽道中流动，无重力，上下边界是周期边界，左右是slip的wall，有压力参考点，位置就在左边图的左侧边界附近的“小气泡”的旁边。时间是3.5s，对应Fig1.压力曲线剧烈波动的那一阶段。
左边的图是alpha.water，中间图是p_rgh，右边图是我添加的斥力。
斥力仅在两个气泡接近到一定程度时，加在气泡边界上，就是右图中的两个小竖条。问题在于斥力生效的时候，压力参考点附近很奇怪的冒出了气体，这些气体似乎是从1号气泡上抽出来的，由于压力参考点附近有了气体，相当于有了一个气泡，所以由于表面张力，流场中的液体的压强就会下降，可以看到中间图中流场大部分p_rgh是负的几十帕，压力参考点附近有个红点。
这是Fig1.中流场压力同步地减小的直接原因。

Fig2. 流场图像。左，alpha.water；中，p_rgh;右，我添加的斥力源项。

除了上述的“抽离气体”行为之外，还有可能凭空生成一些气体，表现为全场积分得到的气体体积超过初始生成的气体体积。
还可能是压力参考点附近出现alpha.water>1的网格，也会导致有grad(alpha.water)，然后导致流场其它地方的压力下降。我之前采用自适应加密时候，看到的监测点压力曲线同步增大减小的直接原因就是这种alpha.water>1的现象。

所以直接原因是压力参考点附近出现了气体，更深层次的原因我感觉是和流场的连续性有关。看计算的log文件，发现出现问题的这一个时间步，time step continuity errors明显大一些。下面是log文件，第一个时间步是没有斥力生效，第二个时间步时候斥力生效了，time step continuity errors也比第一个时间步大了两个数量级。

Courant Number mean: 0.062496661 max: 0.49675831
Interface Courant Number mean: 0.00012625052 max: 0.3551864
deltaT = 5.565631e-06
Time = 0.02073116103

PIMPLE: iteration 1
smoothSolver:  Solving for alpha.water, Initial residual = 0.0006698861, Final residual = 2.4349078e-13, No Iterations 3
Phase-1 volume fraction = 0.99579311  Min(alpha.water) = 1.2926026e-11  Max(alpha.water) = 1.0000346
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.99579311  Min(alpha.water) = -0.0043784079  Max(alpha.water) = 1.0000346
smoothSolver:  Solving for alpha.water, Initial residual = 0.00066896527, Final residual = 2.4380054e-13, No Iterations 3
Phase-1 volume fraction = 0.99579311  Min(alpha.water) = 1.2926035e-11  Max(alpha.water) = 1.0000345
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.99579311  Min(alpha.water) = -0.0043784079  Max(alpha.water) = 1.0000345
GAMG:  Solving for p_rgh, Initial residual = 0.055540806, Final residual = 7.6812173e-05, No Iterations 3
Pressure gradient source: uncorrected Ubar = 1.0363001, pressure gradient = 333.76278
~~time step continuity errors : sum local = 5.1237013e-08, global = -3.751065e-20, cumulative = -4.8212224e-09
The pressure at the reference cell is corrected from 0.13807379 to 0~~
GAMG:  Solving for p_rgh, Initial residual = 0.0012770611, Final residual = 1.733568e-06, No Iterations 5
Pressure gradient source: uncorrected Ubar = 1.0363001, pressure gradient = 335.51285
time step continuity errors : sum local = 1.1477203e-09, global = -3.6709073e-20, cumulative = -4.8212224e-09
The pressure at the reference cell is corrected from 0.0099254384 to 0
GAMG:  Solving for p_rgh, Initial residual = 0.00012244276, Final residual = 9.5463703e-09, No Iterations 43
Pressure gradient source: uncorrected Ubar = 1.0363001, pressure gradient = 334.6004
time step continuity errors : sum local = 6.3417109e-12, global = -3.7236763e-20, cumulative = -4.8212224e-09
The pressure at the reference cell is corrected from 0.018604574 to 0
ExecutionTime = 112420.57 s  ClockTime = 112564 s

fieldAverage fieldAverage write:
    Calculating averages

volFieldValue gas_in_domain write:
    volIntegrate(region0) of alpha.air = 6.7310313e-08

Courant Number mean: 0.062903353 max: 0.50213272
Interface Courant Number mean: 0.00012710096 max: 0.37880418
deltaT = 5.5439467e-06
Time = 0.02073670498

PIMPLE: iteration 1
smoothSolver:  Solving for alpha.water, Initial residual = 0.00066738165, Final residual = 2.5095471e-13, No Iterations 3
Phase-1 volume fraction = 0.99579311  Min(alpha.water) = 1.292604e-11  Max(alpha.water) = 1.0000345
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.99579311  Min(alpha.water) = -0.0043784079  Max(alpha.water) = 1.0000345
smoothSolver:  Solving for alpha.water, Initial residual = 0.00066646251, Final residual = 2.5131431e-13, No Iterations 3
Phase-1 volume fraction = 0.99579311  Min(alpha.water) = 1.2926051e-11  Max(alpha.water) = 1.0000345
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.99579311  Min(alpha.water) = -0.0043784079  Max(alpha.water) = 1.0000345
GAMG:  Solving for p_rgh, Initial residual = 0.75668049, Final residual = 0.00077440161, No Iterations 4
Pressure gradient source: uncorrected Ubar = 1.0362943, pressure gradient = 669.91266
~~time step continuity errors : sum local = 1.9935763e-06, global = -6.2120888e-09, cumulative = -1.1033311e-08
The pressure at the reference cell is corrected from 5523.05 to 0~~
GAMG:  Solving for p_rgh, Initial residual = 0.0010213942, Final residual = 2.0013854e-06, No Iterations 69
Pressure gradient source: uncorrected Ubar = 1.036315, pressure gradient = -548.48592
time step continuity errors : sum local = 1.5198456e-08, global = -6.2120888e-09, cumulative = -1.72454e-08
The pressure at the reference cell is corrected from 5341.8877 to 0
GAMG:  Solving for p_rgh, Initial residual = 7.1997514e-05, Final residual = 8.5095615e-09, No Iterations 86
Pressure gradient source: uncorrected Ubar = 1.036315, pressure gradient = -549.38251
time step continuity errors : sum local = 6.2499269e-09, global = -6.2120888e-09, cumulative = -2.3457489e-08
The pressure at the reference cell is corrected from 5299.7144 to 0
ExecutionTime = 112459.83 s  ClockTime = 112604 s

fieldAverage fieldAverage write:
    Calculating averages

所以有大佬遇到过相同的问题吗？求指点，救救孩子吧，卡在压力这个问题上大半年了🥺🥺🥺🙏🙏🙏

[李东岳]

如果确定是斥力的问题的话，为何不这样：找个1万网格的算例，弄两个气泡进去，单核跑，专门debug你添加的斥力。

[学流体的小明]

在用二维算例debug了，应该是斥力加得太突兀的原因，也许
继续求指点

[李东岳]

斥力方程写一下，或者sci发一下

[学流体的小明]

Short-range repulsive force model for near-contact interactions of bubbles
Lingxin Zhang, Kai Peng, Xueming Shao, and Jian Deng
Phys. Rev. E 106, 045110 – Published 25 October 2022
https://doi.org/10.1103/PhysRevE.106.045110

[李东岳]

1. 遍历所有相分数在0-1直接的网格，并做标记；
2. 在标记的网格上，通过方程10判断是要遍历的网格，取第i个网格，算${\mathbf{F}}_{rep,i}=\sum ^N\frac{K}{h_{ij}^3}\mathbf{n}$；
3. 对i遍历，有做标记网格上的所有的${\mathbf{F}}_{rep,i}$，加到动量方程或压力方程

你这样debug一下：

用你的2D算例，网格搞到几千，然后算一下，

1. 把你的相分数、斥力的云图弄出来，
2. 显示斥力矢量图，类似下面这个图（显示的是网格上的表面张力）

首先要把气泡往外推开的过程复现

---

他们挺有意思，之前还用多相VOF来算这个，一个气泡一哥相，暴力！

Direct numerical simulation of deformable rising bubbles at low Reynolds numbers

[学流体的小明]

斥力矢量图我晚些再搞一下。

---

气泡往外推开是复现了的。
现在发现压力参考点附近的异常是结果，原因还是斥力有些突兀。不使用压力参考点也会出现这种压力振荡。
下面是不使用压力参考点的一次模拟。压力振荡是存在的，而且是和斥力的出现是同步的。

最开始的时候气泡弹开了，如果不加斥力，上面两个气泡是会融合的。
3bbsNoReferencevideo1.avi

压力振荡最厉害的一段时间：3-4s。
3bbsNoReferencevideo2.avi

[李东岳]

一楼从你的描述来看有2个问题：1）空泡发生，2）压力不稳。视频算的可以啊，没看到空泡的发生，第一个问题解决了是吧？

[学流体的小明]

也不算是空泡，其实第二个视频中也可以看到有斥力的那对泡，左边的泡变小了，反正又是一个很奇怪的问题，只是在有压力参考点的时候，缺少的那部分气体跑到了压力参考点附近。
压力不稳是最关键的问题，感觉解决它之后，第一个问题应该也能解决。

[李东岳]

压力那个不太好搞。就像你上次发那个帖子一样。我记得你说3者选2就没问题，3者都开启就不行。我初步只能感觉是那个方向的问题，至于咋处理现在也没啥头绪。

至于相分数变小，倒是可以继续debug一下。

[学流体的小明]

网格上的表面张力是一个surfaceScalarField，要显示的话，看其他人的方法是使用foamToVTK先转化，再用paraview看，感觉有点麻烦。李老师您有什么快速便捷的方法吗？

[李东岳]

surfaceScalarField是看不了，不过他那个表面张力是从体场插值过去的，你把体场输出就好了

[学流体的小明]

现在大致上解决了，确实是添加的斥力太大的问题，之前代码写得不太对，写成了$K/h^3\mathrm{\nabla }\alpha$，而正确的表达式应该是$\frac{K}{h^3}\frac{\mathrm{\nabla }\alpha }{\left|\mathrm{\nabla }\alpha \right|}$，导致某些情况下两个气泡接近时候，由于界面上某些网格的梯度特别大，导致实际添加的斥力会非常非常大。
现在修改了代码，之前的问题很大程度上解决了，还有一点问题是还是会出现添加斥力比较大的时候，因为总会有气泡不断接近到比较近的时候，从而产生一个比较大的斥力，也会导致全局的压力振荡。就比如下图中这些上下跳的曲线。

---

表面张力是$\sigma \kappa \mathrm{\nabla }\alpha$，和斥力的表达式还不一样，我再都算一下volScalarField，不往面上插值什么的，比较一下它俩。

[李东岳]

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

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

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

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

http://dyfluid.com/reactingTwoPhaseEulerFoam.html

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

[学流体的小明]

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

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

我在输出两个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的最大值。

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

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

[李东岳]

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

[学流体的小明]

之前试算了很多，发现压力参考点不是最根本的原因，最根本的原因还是添加了这个斥力源项导致的流场压力不稳定。不用压力参考点也算出来了这种压力振荡。
应该是和具体的代码有关系，和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();
        #};
    }
}