CFD中文网

@Albert-Lee 老师您好，请问您对这个问题有新的理解了吗？谢谢!

博主您好，我用v1812版本得wedge网格旋转一周生成360°的3D网格似乎很顺利，没有遇到博主您说的问题，可能差别在于axis这里：

nLayers             36;

expansionRatio      1.0;

sectorCoeffs
{
    axisPt      (0 0 0.0);
    axis        (0 0 -1);
    angle       360;  // For nLayers=1 assume symmetry so angle/2 on each side
}

但是仿真的时候就有点问题了。请教一下博主，我后续的计算是用wedge网格生成的3D圆管状网格进行两相流动计算，在对称轴位置处，界面波动非常剧烈，（举例来说就像包子的褶皱），请问博主您是否有解决类似问题的经验？

感谢。

@李东岳 东岳老师，我网格160w。
想请教您是否有相关的调整或优化建议？谢谢！
因为我用了第三方求解器（interFlow）

我的环境及求解器信息：
openfoam: v2312
solver: interFlow
AMR: BlastAMR

@李东岳 东岳老师好，感谢您的关注和回复。
这里是我之前有点表述得太混乱。

简单来说，上述第二个图里梯形的加密是我用SHM加进去的，相当于：第一个图是AMR加密时的结果，第二个图是SHM加密时的结果（没用AMR）。我在左上角新加了标注。

但不管是哪种加密，都有一个问题：在网格加密等级出现变化的地方，会出现不连续的解场。
并且我对比发现要么是（p连续，则U不连续），要么是（U连续，则p不连续）。

由于SHM的结果也是这样，所以我猜应该和AMR没有关系，可能是我的求解设置哪里有点问题，所以发帖咨询一下大家。

@李东岳东岳老师，各位老师，大家好。

最近尝试使用AMR和SHM算同一个算例，但是输出的时间步结果中，p或者U总有一个在加密网格周围不连续，显示出非物理的结果。

我最开始怀疑这可能是因为解场在网格细化和重新平衡后需要多个时间步才能完全调整和平滑（这个想法很呆），还向blastAMR库的作者留言咨询了一下，他说“Generally, we cannot run a CFD simulation for even a single step with a degrading field - it will mess up the accuracy from that point onward; so, this is not allowed to happen.”

从气泡的运动和计算结果来看，AMR和SHM的结果都能和实验结果对比得很好。
所以这会不会是结果输出或者通量校正哪里没设置好，比如correcPhi, correctFluxes之类的，不知道大家有没有遇到类似的问题，想请大家指导，给点经验，谢谢!

💡 AMR结果（展示的是后半段结果）：

💡 SHM结果（展示的是前半段结果）：

以下是fvSchemes文件信息

/·························分割线·························/
ddtSchemes
{
    default         CrankNicolson 0.9;
}

gradSchemes
{
    default         pointCellsLeastSquares;//Gauss linear;//
}

divSchemes
{
    default             none;
    div(rhoPhi,U)	Gauss limitedLinearV 1;
    div(phi,alpha)	Gauss vanLeer;
    div(phirb,alpha)	Gauss interfaceCompression;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default	Gauss linear corrected;
}

interpolationSchemes
{
    default	linear;
}

snGradSchemes
{
    default	corrected;
}

fluxRequired
{
    default         no;
    p_rgh;
    pcorr;
    alpha.water;
}

以下是fvSolution文件信息

solvers
{
    "alpha.water.*"
    {
        advectionScheme isoAdvection;//geoAdvection;//
        reconstructionScheme plicRDF; //RDFadvect isoInverseDistance

        // advectionScheme MULESScheme;//geoAdvection;//
        // reconstructionScheme isoAlpha; //plicRDF; //RDFadvect isoInverseDistance

        // isoAdvector parameters
        vof2IsoTol      1e-6;
        surfCellTol     1e-6;
        nAlphaBounds    3;
        snapTol         1e-12;
        clip            true;
        gradAlphaNormal false;
        writeVTK        false;

        // MULES parameters
        nAlphaCorr      2;
        cAlpha          1;
        icAlpha         0;
        MULESCorr       yes;
        nLimiterIter    3;
        alphaApplyPrevCorr yes;
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-8;
        relTol          0;
        minIter         1;

        // parameters read by both isoAdvector and MULES
        nAlphaSubCycles 1;
    }

    pcorr
    {
        solver          GAMG;
        smoother        DIC;
        tolerance       1e-10;
        relTol          0;
    }

    pcorrFinal
    {
        $pcorr;
        relTol          0;
    }

    p_rgh
    {
        solver          GAMG;
        smoother        DIC;
        tolerance       1e-07;
        relTol          0.001;
    }

    p_rghFinal
    {
        $p_rgh;
        tolerance       1e-8; /* Should be 1 to 2 orders of 
                                 magnitude smaller than 
                                 alpha.phase1.*.surfCellTol */
        relTol          0;
    }

    "(U|k|epsilon).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-07;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor   yes;
    nOuterCorrectors    1;
    nCorrectors         3;
    nNonOrthogonalCorrectors 1;
    correctPhi true;
    pRefCell            0;
    pRefValue           0;
}

relaxationFactors
{
    equations
    {
        ".*"            1;
    }
}

补充上我的环境及求解器信息：

openfoam: v2312
solver: interFlow
AMR: BlastAMR

@capillaryFix 好的老师，感谢！另外，如果可以，您能方便描述一下此前植入contact line pinning效果不好的具体表现吗？比如是接触线固定不住，还是像这篇帖子提到的这种“界面褶皱”，想听听前辈您的经验。

@李东岳 东岳老师，这两天我系统调研了一下动态接触角模型相关的文献，后续我整理一下发出来。
这里先回复您：Model B （Stick-slip model）的公式是这个：

刚好这个Stick-slip model的公式好像是动态接触角模型中形式最为简单的，不知道这个公式的植入好不好弄？想请老师您帮忙看看，膜拜巨蟹！

@capillaryFix 谢谢老师您的回复。

如果可以的话，您植入过的contact line pinning可以分享一份吗？巨蟹！
763400067@qq.com

此外，我发现把气泡底部的网格SHM多加密一道（就是涉及接触角的周围网格），影响就会小很多，您需要的话可以试试。

@李东岳 东岳老师好，感谢您的关注和回复。
这个是需要植入一个新的接触角条件，具体的是（Yuming Chen, 2009IJMF）那张截图中框起来的那个。

另外，我还没太关注原生 interFoam 是怎么处理这个问题的，我尽快再调研下，然后补充回复您。

东岳老师@李东岳 各位大佬@capillaryFix ，大家好。

我对这个问题调研有一段时间了。最开始发帖咨询这个问题，是源于在使用第三方求解器仿真浸没孔口注气气泡生长的时候，发现孔口周围的接触线那里有点奇怪。这两天在站里又逛了下，觉得这个帖子的问题也有点类似。结合最近的调研，我觉得原因可能在接触角边界条件。
现在简单总结一下自己的思考和想法，可能有错，新看到的朋友一定要慎重参考。

/·························分割线·························/

先汇报一下和Henning博士(TwoPhaseFlow求解器作者)在github上的沟通情况。

• 插入一段Henning博士的信息，他和IsoAdvector作者之一的Roenby合作很密切，但拿不准他们之间的团队关系，总之Henning博士水平肯定很顶。

Henning博士的回复其实可以总结为两点：

1. interIsoFoam 不考虑接触角，当接触角不等于 90 度时，底部不会根据接触角来作修正，所以结果也肯定不对（这也是为什么interIsoFoam仿真注气气泡时最底部不会出现薄片褶皱。
   
   2.Henning对于他的求解器中接触角的实现是很有信心的，但为什么出现这个问题，他也觉得奇怪，但他没时间特地dubug。

/·························分割线·························/

现在整体调研下来后，我的观点是Henning博士的求解器没有问题，问题出在边界条件上。解决的办法是要新写一个边界条件，以实现接触线“钉扎效应”。

💡先介绍一下接触线钉扎效应，flow3d公司的这个仿真案例用液滴做了个很好的介绍：


这篇文献专门对液滴的钉扎效应做了边界条件处理，其所提供的附加材料更直观。

Cai, Z., Song, Y., 2021. Implementing Contact Angle Hysteresis in Moving Mesh-Based Two-Phase Flow Numerical Simulations. ACS Omega 6, 35711–35717.

简单来说，我理解的钉扎效应就是允许接触角在指定的范围内角度值可以任意变，但是接触线的位置不能动。

💡那么对于气泡，当注气针管的管壁无限薄，或者润湿性特别高的浸没孔口注气（许多文献称之为Mode A气泡，如下图a），也属于接触线钉扎效应的问题，因为其接触线固定在孔口，接触角可以随意改变，因此不该用constantAlphaContactAngle静态接触角边界条件仿真这个问题。

更细节的，用constantAlphaContactAngle静态接触角边界条件，仿真注气气泡生长会导致的问题可以借助(Albadawi et al., 2013CES)的水中注气气泡生长的实验参数进行的OpenFOAM仿真结果来说明，在注气流量非常低，即表面张力作用影响显著的工况下：

1.气泡生长阶段：可以看到使用constantAlphaContactAngle边界条件，因为interFlow求解器几何VOF方法要更正界面接触角，但是这种更正相当于强行按压界面回到指定接触角度，这就导致了气泡生长过程中界面不断在波动，并且速度分布也特别奇怪。强调一下这里展示的还是楔形计算的结果，三维结果会更加明显。

2.气泡颈缩脱离阶段，就更奇怪了，颈缩期间出现了十分明显的褶皱，并且在这些褶皱的地方，由于曲率变化太大，很多点出现压力激增，最终导致气泡界面变得扭曲，没法看。

所以我认为问题出在边界条件上。对于这种Mode A 气泡的生长，解决的办法是要新写一个边界条件，以实现接触线“钉扎效应”。

/·························分割线·························/

OK，紧接着从文献里看看几个学术团队对于这种Mode A 气泡的接触线边界条件是怎么处理的：

①💡 Yuming Chen，University of Stuttgart 
Chen, Y., Mertz, R., Kulenovic, R., 2009. Numerical simulation of bubble formation on orifice plates with a moving contact line. Int. J. Multiphase Flow 35, 66–77.

（注意框起来的是模型B，这里是在介绍模型，和Mode A气泡没有必然联系，即：并不是Mode A气泡就对应Model A接触线模型）

②💡 Kuipers团队，Eindhoven University of Technology
Mirsandi, H., Rajkotwala, A.H., Baltussen, M.W., Peters, E.A.J.F., Kuipers, J.A.M., 2018. Numerical simulation of bubble formation with a moving contact line using Local Front Reconstruction Method. Chem. Eng. Sci. 187, 415–431.

然后...然后就开始有点吃力了，之前还编译了这个帖子中的边界条件dynamicAlphaContactAngleMikeWorth，结果他那个好像说得不是一回事。
现在自己想基于OpenFOAM自带的两个基础的接触角边界条件来改，但对constantAlphaContactAngle 和dynamicAlphaContactAngle的代码有点看不懂，请东岳老师，以及各位大佬出手指点，巨蟹！

此外，贴子码了一大篇，有很多自己的不严谨观点，请大家批评指正，这样能尽快学到东西。

@capillaryFix 下午好，铁子。您那边有新的解决方案和想法吗？

@浪迹天大 老师您好，按照这个方法，我需要更改controlDict中的maxDeltaT。但是会发现文件更新不起作用。

/*--------------------------------*- 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          controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     interFlowVb1;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         0.3;

deltaT          0.00005;

writeControl    adjustableRunTime;

writeInterval   0.05;

purgeWrite      0;

writeFormat     binary;

writePrecision  8;

writeCompression off;

timeFormat      general;

timePrecision   8;

runTimeModifiable yes;

adjustTimeStep  yes;

maxCo           0.5;
maxAlphaCo      0.5;
maxCapillaryNum 10; // stability criteria is 1

 maxDeltaT       2e-5;//controlDict.a;
// maxDeltaT       1e-5;//controlDict.b;
// maxDeltaT       5e-6;//controlDict.c;


functions
{
    fileUpdate1
    {
        type              timeActivatedFileUpdate;
        libs              ("libutilityFunctionObjects.so");
        writeControl      timeStep;
        writeInterval     1;
        fileToUpdate      "$FOAM_CASE/system/controlDict";
        timeVsFile
        (
            (-1 "$FOAM_CASE/system/controlDict.a")
            (1e-4 "$FOAM_CASE/system/controlDict.b")
            (2e-4 "$FOAM_CASE/system/controlDict.c")
        );
    }
}

三个文件只有maxDeltaT有差别，请老师捞捞@李东岳 @浪迹天大

运行后结果提示：

timeActivatedFileUpdate: copying file
"~/orifice3Da/system/controlDict.c"
to:
"~/orifice3Da/system/controlDict"

Courant Number mean: 2.0171174e-08 max: 0.3369988
Interface Courant Number mean: 2.978204e-10 max: 0.22052861
Capillary Number: 2.4526254
deltaT = 2e-05
Time = 0.00022

文件有更新，并且maxDeltaT也更新为了1e-5; 但是它没有发挥作用，可以看到文件更新后的deltaT 依然是 2e-05，如果是我自己运行过程中手动改小maxDeltaT，它会马上从deltaT = 2e-05，便为deltaT = 1e-05。

@好运来辣 我以前测试的时候，记得paraview不支持中文路径，注意检查这个影响因素。
包括paraview的安装路径（盘符有中文名称不影响），以及仿真文件的所在路径，并且文件名中的特殊符号最好也避免一下。

@好运来辣
直接去paraview的官方网页下载
paraview官网
最好下载绿色版的直接解压就能用，win服务器没试过，需要你自己试一试。可以的话回来跟帖和大家分享一下
具体的见图

没看懂你具体问的是啥问题，是不是问windows下的wsl2的具体地址？

你在wsl2命令行输这个命令试试：

explorer.exe .

@李东岳 谢谢东岳老师，我用的就是v2106，然后我看了下我按照v2106的211215的补丁，用wmake -build-info看了下补丁版本是0，自己尝试手动打补丁（替换文件），但是最终没起效果。

最后干脆用了简单粗暴的办法，把openfoam-OpenFOAM-v2106.211215这个版本下载下来之后直接编译。
依然感谢东岳老师。

各位老师，最近在GitHub发现一位作者开源了它的求解器，其中明确提到需要用OpenFOAM v2106
我尝试用东岳老师提供的OpenFOAM全系列虚拟机下载中的v2106编译，一直报错，通过不了。

后来查询到作者有这么一句提示：
If using a pre-compiled version of OpenFOAM v2106, use pre-compiled patch 211215.

我对应去看了这个v2106的211215的补丁

请问一下各位老师，怎么打这个补丁，甚至是如何查看自己的版本是否已经打过这个补丁了？

谢谢各位老师 这个打补丁 完全不知道怎么下手
感谢东岳老师的全系列全集 非常好用@李东岳

@chengan-wang 投影不清楚怎么弄，缩放行不行？用transform把x轴的缩放因子改成0.00001，你看看效果？

@lyc 你好，我很少遇到。即便遇到都是重启就能解决的那种。

@lwjetmann 老师您好，再请教一下您。

1. 用wallshearstress后处理算壁面粘性剪切力，然后在paraview中计算出在壁面各点法向和切向的投影；
2. 基于wallshearstress functionObject进行修改，输出投影后的结果

按照您提到的上述方案1，对于圆球绕流（球体界面是“noSlip”边界条件），我顺利的输出了粘性剪切力，结果也验证的上。

想请您再帮忙看看另外一种情况：绕流气泡，即将边界条件改为Symmetry：
因为wallshearstress只对wall类型的边界计算粘性剪应力，因此不能再通过【方案1】输出界面上的粘性剪应力结果。
关于设置Symmetry边界条件的原因是源于这里：Stokes流中的圆球绕流的阻力系数，怎么设置才能算准。

对此，针对绕流气泡的结果（边界条件为Symmetry）：

我尝试手动在结果中将Symmetry边界条件改为Wall，然后执行一次wallshearstress；
通过对比文献中气泡的压差阻力系数、法向粘性应力和切向粘性应力结果，输出的结果并不对

除此以外，是否只有像您提到的【方案二】基于wallshearstress functionObject进行修改，才能输出“绕流气泡”的粘性剪切力结果。
如果您方便的话，能否请您指点一些修改思路，谢谢老师。