找到vof中interface的位置



  • 大家好,我目前在用VOF-DEM耦合做气泡(VOF)和颗粒(DEM)相互作用的模拟。现在想添加一个颗粒和气泡间的作用力。这个作用力由颗粒与气液interface间的距离决定。目前想用alpha1=0.5作为interface,想请教一下大家,如何在代码层面尽量快速高效的找到alpha1=0.5的所有网格的位置,要遍历所有网格去找吗?


  • 自由表面模型副教授 OpenFOAM教授

    求解过程中的0.5等值面提取没搞过啊,只有后处理的时候提取过iso-surface。你说的这个需要在求解过程中确定interface的位置,使用OpenFoam最新的iso-Advector或者PLIC-VOF吧,这两个都属于几何重构方法,求解过程中可以得到interface的位置信息。



  • 这个作用力由颗粒与气液interface间的距离决定

    你要处理颗粒一半在水里、一半在空气中的情况?
    是的。你这个在VOF里面确切的获得界面不容易,主要是VOF目前很难有严格的alpha=0.5的情况,除非采用其他算法几何性的重构出来。



  • @队长别开枪 没用过iso-Advector,主要我还改了相方程(加了颗粒相体积分数的影响),不知道iso-Advector相方程好修改吗?



  • @东岳 不一定是一半一半的情况,就是颗粒距离气泡表面有一定的距离的时候,会有力的作用,类似于MPPICinterFOAM里面在界面处施加的里,只不过我是跟界面距离有关,MPPICinterFOAM是相分数梯度的函数。


  • 自由表面模型副教授 OpenFOAM教授

    @linhan-ge 不知道你的相方程怎么修改的,iso-Advector 修改起来不难,我自己在做多面体网格的PLIC解析算法,advection部分因为不是我的研究范围就用的iso-Advector提供的算法。



  • @队长别开枪 其实我的修改很简单,就是加入一个颗粒相体积分数的影响,0_1528067861730_fc6c44b9-e94e-4d5f-a7c1-7102a25b4c1d-image.png 这是对MULES方法的修改,针对isoadvector,应该就需要前两项。不知道以大神的经验,可行性如何,代码实现难度大不大。我目前对isoadvector的了解非常浅薄。


  • 自由表面模型副教授 OpenFOAM教授

    @linhan-ge $\alpha_f$ 是指颗粒相体积分数吗?因为一般下标$f$指代有限控制体的边界(face)



  • @队长别开枪 $\alpha_f$是所有流体的体积分数,包括空气和水。$\alpha_f$=$\alpha_1$+$\alpha_2$



  • @队长别开枪 你好,因为研究内容界面张力主导,需要两相界面更精确一点,想要在OF中实现结构化网格的PLIC界面重构,找了一些资料,包括Fortran版的实现,几何计算实在看得不是很明白:qichuang:
    请教一下有讲得好一点的资料推荐吗?计算几何的知识可以怎么补一下?感觉很多论文里网格层面操作很多啊:mihu:



  • @wallong 为什么不用iso-advector 啊?有什么特别的考虑吗?


  • 自由表面模型副教授 OpenFOAM教授

    @wallong 结构网格中的PLIC算法相对容易,PLIC说白了就是求解interface近似平面的位置,数学上一般使用$\vec{n} \cdot \vec{X} + D = 0$表示,需要求解interface orientation vector $\vec{n}$ 和 signed distance $D$,$\vec{n}$一般使用主相体积分数的梯度,即$\vec{n} = - \frac{\nabla \alpha}{| \nabla \alpha |}$,也有使用其他辅助函数的,例如CLSVOF方法,使用的是level-set函数的梯度计算$\vec{n}$,优点是光滑。总的来说,$\vec{n}$的计算主要是梯度的计算,在非结构网格中无非高斯定理或者最小二乘,而$D$的计算就和网格单元类型相关了,也更复杂一些,六面体单元算是比较简单的一类了,可以参考 http://iccfd9.itu.edu.tr/assets/pdf/papers/ICCFD9-2016-288.pdf 。所有的网格操作都是为了计算$D$,结构网格最简单。
    iso-advector使用的是 iso-surface的基本思想,先假定isoValue = 0.5,将interface cell中的isoValue 等值面作为初始interface近似,然后调整isoValue 的数值使得等值面的位置满足主相体积分数$alpha$的值,isoValue 相当于PLIC中的$D$。与PLIC相比,iso-advector重构出来的interface近似不一定是平面,会有翘曲。


  • 自由表面模型副教授 OpenFOAM教授

    @linhan-ge 这样的话没有颗粒的话$\alpha_f = 1$,考虑颗粒的影响的话$\alpha_f = \alpha_1 + \alpha_2$就小于1了,我的理解对吗?如果我的理解对的话,那情况就复杂了。。。



  • @队长别开枪 抱歉,这个地方搞错了,$\alpha_1+\alpha_2=1$是始终成立的,是液体或气体占流体的体积分数,不是占三相的,$\alpha_f$实际上是流体占流固两相的体积分数。



  • @linhan-ge iso-advector很棒,不过考虑到自己添加自适应功能的难度,观望中...



  • @linhan-ge 作者一直在这方面努力,好像最近有新的进展了,上周Yoube更新了新的AMR视频,之前说六月可能会发布 https://www.youtube.com/channel/UCt6Idpv4C8TTgz1iUX0prAA



  • @队长别开枪 非常感谢,很有帮助,抱歉回复得有点晚了
    我起初参考论文[1],通过PLIC重构,进而得到一个每个面上AOF(area of fraction),来作为alphaEqn.H中的phiAlpha相可以显著提高界面精度。此外论文中实现了非结构化网格的PLIC,依据的是Maric T[2]的迭代方法,这部分可能和你的工作相关。
    isoAdvector的实现确实很有参考意义,学习中,最近作者好像在测试AMR了,可能不久会发布。
    之前Maric T说过等完成了会公开Code,见https://www.cfd-online.com/Forums/openfoam-programming-development/89713-plic-interfoam.html,后来似乎没有。猜测可能出书,办培训班的原因,看了一眼去年他的博士论文,里面涉及了PLIC和其他界面追踪方法,个人底子差,确实难懂...
    关于我的课题,除了精确界面,还需要关注的是界面张力项,CSF模型有更大影响,所以在试着实现一个新的张力模型,在你给的论文中作者反复提到了“ without smearing, dispersing or wrinkling.”,大家都知道OpenFOAM原来的方法在smearing方面表现不好,其余wrinking方面是不是也表现不好,在2-D vortex测试中会出现不明锯齿状界面,那dispersing方面呢?
    [1]Dianat M, Skarysz M, Garmory A. A Coupled Level Set and Volume of Fluid method for automotive exterior water management applications[J]. International Journal of Multiphase Flow, 2017, 91: 19-38.
    [2]Maric T, Marschall H, Bothe D. voFoam-a geometrical volume of fluid algorithm on arbitrary unstructured meshes with local dynamic adaptive mesh refinement using OpenFOAM[J]. arXiv preprint arXiv:1305.3417, 2013.



  • @wallong 看各位高手讨论学到很多,对于我们这种主要侧重于应用这些方法的人来说,跟不上你们的节奏:zoule: 。大神对于我上面提出的问题,即在相方程中加入颗粒相的影响,有什么看法吗?代码层面的实现难度怎么样?由于有毕业压力,目前不敢轻易的换方法。



  • @linhan-ge 我顶多算刚入门,聊八卦有点多,没怎么聊正事:mianmo:
    VOF-DEM不了解,不好说,有什么影响,难度如何,建议先试一试吧,试下就知道了。下次可以把参考论文post出来



  • 就是颗粒距离气泡表面有一定的距离的时候,会有力的作用,

    建议直观的,用公式或者图的方式说明一下



  • @东岳 0_1528431512437_7454f070-5718-44ac-9088-6072e19772e0-image.png
    这是我们组之前用DEM做的,中间蓝色的是气泡,灰色的是颗粒,都用DEM 来模拟。当颗粒在离气泡表面一定距离时,会有一些表面力,如范德华力,疏水力等,lubrication force 等等,他们的大小都是颗粒与气泡表面距离的函数。举个例子,疏水力的公式是
    0_1528432104731_67f6afbf-84b6-4eec-896d-de5702b8d429-image.png
    其中k,$\lamda$,Rp都是已知的参数,H就是颗粒表面与气泡表面的距离,n是气泡径中心与颗粒中心相连的方向向量。



  • @东岳 0_1528432472586_e0564bbf-15f0-483b-a40f-ba221dc4dcad-image.png
    这个图可能好点,也是我们组做的,但是是油滴和颗粒,还是DEM模拟。


  • 自由表面模型副教授 OpenFOAM教授

    @linhan-ge emmm...这个说实话超出我的研究范围了,可以考虑将$\alpha_f$代入方程(8),然后进行简化,难点在对流项的计算,使用有限体积法转化为面积分后就是控制体表面的流量的计算,无论几何重构还是代数方法都是为这个服务的,但是我没有你这个方向上的经验,因此不好给出任何建议,代码上难度还可以,我在考虑写一个iso-advector的代码解析,不过有毕业压力,可能近期不会完成了。。。


  • 自由表面模型副教授 OpenFOAM教授

    @wallong dispersing我认为这是VOF方法天生的缺陷,无论对于几何重构还是代数方法,只能增加网格分辨率,减少时间步长,别的没啥好的解决办法。几何重构相对代数方法一定程度上克服了代数方法的numerical diffusion,根子还是VOF函数在interface处是个不连续函数。
    我们组有一部分工作是PLIC-VOF中任意多面体单元界面重构的解析算法,文献[2]中用的迭代法求解效率不适合用在实际应用中,研究成果发表后我们将会公布源代码。
    我们在研究中发现现有的CSF模型在capillary force主导的流动中存在计算的表面张力偏大的问题。



  • @队长别开枪 多谢解答。修改MULES时非常简单,直接在相应位置把通量或体积分数替换掉即可。不知道在isoadvector中是不是也是类似的操作。:mihu:



  • @linhan-ge找到vof中interface的位置 中说:

    H就是颗粒表面与气泡表面的距离,n是气泡径中心与颗粒中心相连的方向向量。

    以疏水力举例,主要就是定义H和n。提醒几点:

    • 如要计算颗粒表面与气泡表面的距离H,最好先求n
    • 在求n的时候,颗粒中心有了。你的气泡中心,目前OpenFOAM里面是VOF,这个是重构出来的。你要精准的获得气泡中心,你需要非常高的网格分辨率,这样才能比较准。如果有了比较高的分辨率,可以通过alpha值提取出气泡边界,通过这个边界,依据你们的算法计算出中心位置。这样n向量就有了
    • 有了n之后,H要顺着n求,这个比较好求

    个人感觉,你们这个网格分辨率要很大。VOF是直接模拟。DEM是介尺度。假定一个气泡4mm,那网格如果是0.1mm的话,你的颗粒直径最好小于0.01mm左右。DEM这面求解精度高度要求网格大于颗粒直径,VOF那面要求网格远小于气泡直径。可以衡量下。

    比如你第二个图,液滴和颗粒这个图,这个图液滴和颗粒基本同量级。网格方面会存在冲突。网格太密DEM那面不好用,网格太粗液滴界面出不来。如果模拟单气泡,应该比较好处理。

    PS 图片很直观,经典。

    如何在代码层面尽量快速高效的找到alpha1=0.5的所有网格的位置,要遍历所有网格去找吗?

    回到你最开始的问题,这个需要有非常好的网格分辨率。比如网格尺度是气泡尺度的百分之一或千分之一。才能有这个界面。如果网格糙,界面就存在厚度。你可以寻找$0.2<\alpha <0.8$,类似这样给两个用户自定义阈值。

    目前,能想到的只有遍历网格去找。另一种方法是参考DEM那面,DEM那面通过粒子穿过网格面的算法去定位粒子所在的网格位置。But,非常复杂且有点小bug,有时候粒子沿着网格定点对角线穿过的时候粒子就丢了。



  • @东岳 非常感谢如此详尽的回复,高手就是高手,直戳问题的要点,信息量较大,我需要消化消化。由于DEM的问题,目前我肯定无法使用非常细的网格,在另一篇帖子里,我想尝试用dual grid来解决这个问题,但可能毕业前没有时间尝试了。在力的处理上,我感觉可能需要用简化的模型代替,来回避对精确界面的依赖。我再仔细琢磨琢磨您的回复。



  • @linhan-ge 有什么进展么最近



  • @东岳 最终还是通过简化力的模型来解决这个问题。正如你之前说的,通过找到界面进行几何分析来计算相应的力需要网格分别率很高, 但是再小的网格也比颗粒要大,这就无法提供足够精确的几何关系。



  • @队长别开枪 您是做带有几何重构的VOF算法的,我想请问一下,现有的算法都是在不可压缩流体中使用,对于可压缩流体,这些几何重构算法(比如PLIC)还能适用吗,为什么现在没有相应的文章呢,是没有人做,还是说在算法上有难以解决的问题呢,这个问题困扰我好久了。


  • 自由表面模型副教授 OpenFOAM教授

    @litong189456 目前我的工作还没考虑气相的可压缩性,考虑可压缩性后,控制体内的密度就不再是VOF的函数了,一些控制方程也要改,应该不容易搞





  • @队长别开枪 最近想钻研下isoadvector,想请教大神一些基本问题。
    相方程中的时间项在isoadvector代码的哪个部分求解的呢?找了一圈只找到advection那项。
    @队长别开枪找到vof中interface的位置 中说:

    求解过程中的0.5等值面提取没搞过啊,只有后处理的时候提取过iso-surface。你说的这个需要在求解过程中确定interface的位置,使用OpenFoam最新的iso-Advector或者PLIC-VOF吧,这两个都属于几何重构方法,求解过程中可以得到interface的位置信息。

    这里,您说的interface的位置信息具体是怎么得到的呢,代码上如何处理?


  • 自由表面模型副教授 OpenFOAM教授

    @linhan-ge 绝大部分几何重构方法都是使用一阶欧拉显式方法处理VOF方程中的时间项的,具体到代码就是

        // Initialising dVf with upwind values
        // i.e. phi[facei]*alpha1[upwindCell[facei]]*dt
        dVf_ = upwind<scalar>(mesh_, phi_).flux(alpha1_)*mesh_.time().deltaT();
    
        // Do the isoAdvection on surface cells
        timeIntegratedFlux();
    
        // Synchronize processor patches
        syncProcPatches(dVf_, phi_);
    
        // Adjust dVf for unbounded cells
        limitFluxes();
    
        // Advect the free surface
        alpha1_ -= fvc::surfaceIntegrate(dVf_);
    

    fvc::surfaceIntegrate(dVf_)就是用来计算当前时间步流过控制体表面的主相体积占控制体体积的比值。

    isoAdvector获取重构好的interface比较困难。在PLIC方法中则会定义交界面为$\vec{n} \cdot \vec{X} + D_0 = 0$,通过求解$\vec{n}$和$D_0$得到interface的近似平面。我们组最近会开源一个基于PLIC方法的二维求解器,在里面可以很方便地得到这个interface的空间位置,三维的预计年底开源,文章正在审。参考 https://doi.org/10.1002/fld.4664



  • 这个蛮好的嘛

    0_1531965646653_捕获.JPG


  • 自由表面模型副教授 OpenFOAM教授

    @东岳 在前人的基础上做的研究,只求老板给毕业 :xinlei:



  • @队长别开枪 非常有用的信息,期待你们求解器的开源。



  • @队长别开枪
    你主要做界面这面的哪个模型?PLIC?


  • 自由表面模型副教授 OpenFOAM教授

    @东岳 做的比较杂,主要是PLIC,提高多边形多面体网格中PLIC的求解效率,PLIC + overset mesh,使用PLIC研究气蚀,就这些



  • @队长别开枪 队教授:chigua: 我目前在使用VOF代数重构的方法在做相关研究,对于几何重构的优点不是很明白,PLIC算法在of中的应用能否提供简单的例子学习一下啊,非常感谢


  • 自由表面模型副教授 OpenFOAM教授

    @fireztw 给个图吧 344e0754-2b10-484f-8441-cfd15db9d26e-image.png
    左边是isoAdvector,属于几何重构方法里的一种,右边代数方法。纠正一下说法,代数方法没有重构这一步的,优点在于求解效率和网格单元形状的适应性。



  • @队长别开枪 大神您好 刚读了上面您的PLIC-VOF的论文 请问这套代码开源了吗 谢谢


  • 自由表面模型副教授 OpenFOAM教授

    @nbyjn 留个邮箱。



  • @队长别开枪 992906307@qq.com 万分感谢! 另外想请教您一个问题 plic界面重构方法和动网格(没有拓扑变形的动网格)会有冲突吗 因为上文提到的isoadvector方法只能用于一部分动网格 不知道plic方法会不会有类似问题


  • 自由表面模型副教授 OpenFOAM教授

    @nbyjn PLIC可以用于重叠网格,其他动网格技术还没试过,不好评价。



  • @队长别开枪 @nbyjn 最近在计算液滴撞击固体壁面的飞溅现象,使用Fluent可以产生飞溅,interFoam没有该现象(修改了壁面接触角模型),而且isoAdvector表面张力和压力平衡又不太好。考虑应该是VOF界面处理引起的,可以求一份PLIC的代码吗?1611689516@qq.com 不胜感激。


  • 自由表面模型副教授 OpenFOAM教授

    @wallong找到vof中interface的位置 中说:

    1611689516@qq.com

    已发。不过advection部分我们暂时使用的iso-advector的时间积分方案,可能对你的问题帮助有限。



  • @队长别开枪 好的,我研究一下,十分感谢!



  • @wallong找到vof中interface的位置 中说:

    而且isoAdvector表面张力和压力平衡又不太好。考虑应该是VOF界面处理引起的

    这个结论是您测试的结果还是有相关文献呢?



  • @linhan-ge 测试过,另外可以参看一下Johan Roenby的说法:"For you surface tension people, I believe you will also need an improved estimate of the curvature, right? This information is somehow contained in the difference in surface normal direction between adjacent surface cells. Hmm... still thinking about how to cook this up in a consistent and robust way for general meshes..." 来自https://www.cfd-online.com/Forums/openfoam-news-announcements-other/180439-isoadvector-release.html


Log in to reply
 

京ICP备15017992号-2