找到vof中interface的位置
-
@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近似不一定是平面,会有翘曲。 -
@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 dispersing我认为这是VOF方法天生的缺陷,无论对于几何重构还是代数方法,只能增加网格分辨率,减少时间步长,别的没啥好的解决办法。几何重构相对代数方法一定程度上克服了代数方法的numerical diffusion,根子还是VOF函数在interface处是个不连续函数。
我们组有一部分工作是PLIC-VOF中任意多面体单元界面重构的解析算法,文献[2]中用的迭代法求解效率不适合用在实际应用中,研究成果发表后我们将会公布源代码。
我们在研究中发现现有的CSF模型在capillary force主导的流动中存在计算的表面张力偏大的问题。 -
@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,有时候粒子沿着网格定点对角线穿过的时候粒子就丢了。
-
@litong189456 目前我的工作还没考虑气相的可压缩性,考虑可压缩性后,控制体内的密度就不再是VOF的函数了,一些控制方程也要改,应该不容易搞
-
@队长别开枪 最近想钻研下isoadvector,想请教大神一些基本问题。
相方程中的时间项在isoadvector代码的哪个部分求解的呢?找了一圈只找到advection那项。
@队长别开枪 在 找到vof中interface的位置 中说:求解过程中的0.5等值面提取没搞过啊,只有后处理的时候提取过iso-surface。你说的这个需要在求解过程中确定interface的位置,使用OpenFoam最新的iso-Advector或者PLIC-VOF吧,这两个都属于几何重构方法,求解过程中可以得到interface的位置信息。
这里,您说的interface的位置信息具体是怎么得到的呢,代码上如何处理?
-
@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 -
-
@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
CFD中文网2016-2023 | 京ICP备15017992号-2