OpenFOAM 11（2020）：提高颗粒检测的VOF模拟算法改进

软件X 11（2020）100382原始软件出版物使用OpenFOAMAndy Spitzenberger，Sebastian Neumann，Martin Heinrich，Rüdiger Schwarze弗赖贝格工业大学力学与流体动力学研究所，Lampadiusstr.4，09599 Freiberg，Germanyar t i cl e i nf o文章历史记录：收到2019年收到修订版2019年11月8日接受2019年12月9日保留字：VOF欧拉粒子检测OpenFOAM数值验证a b st ra ct这项工作的目的是提高OpenFOAM中包含的流体体积（VOF）模拟相场中的颗粒检测。颗粒特性的可靠评估是一个先决条件用于将初级破碎的VOF方法与次级破碎的欧拉-拉格朗日模拟相耦合。为了实现正确的检测结果的代码的相关变化。四个数值算例表明，改进后的算法能够准确地检测出评价层上的粒子及其大小、位置和速度©2019作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本V1此代码版本使用的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2019_260Code Ocean compute capsule无法律代码许可证GNU GPL V3使用git的代码版本控制系统使用C++的软件代码语言、工具和服务编译要求，操作环境依赖OpenFOAM v1906如果可用，链接到开发人员文档/手册问题支持电子邮件Andy. imfd.tu-freiberg.de1. 动机和意义雾化在许多工业应用中是重要的，例如内燃机或金属粉末的生产。特别令人感兴趣的是可控制的崩解过程，其导致特定的粒度分布。特别是在金属粉末的生产中，熔体的高温使得实验难以进行因此，计算流体动力学（CFD）是一个很有前途的替代方案，以获得深入了解复杂的物理雾化。在雾化过程中，液体射流通过喷嘴出现在最初完整的液核的界面上，由于气动力或湍流而出现不稳定性。这些不稳定性增长，并最终导致液体核心分裂成韧带。雾化过程的这一部分称为初级分裂。*通讯作者。电子邮件地址：Andy. imfd.tu-freiberg.de（A. Spitzenberger）。https://doi.org/10.1016/j.softx.2019.100382在二次破碎期间，由一次破碎产生的韧带在空气动力的持续影响二次破碎发生，直到表面张力变得占主导地位，并提供稳定和接近球形的液滴[1，2]。通常使用界面捕获方法（ICM）（如标记物和细胞[3]、水平集[4，5]或流体体积（VOF）[6]）对初次破裂进行建模。所有这些方法都能够重新求解任意形状的界面，从而适用于稠密喷雾区的液核和非球形韧带由于液滴必须通过网格充分良好地分辨，因此这些方法与高计算成本相关联[2，7]。由于同样的原因，ICM不适合应用于稀喷雾区，在稀喷雾区，液滴很小，接近球形且分散。解决大量的小液滴将大大增加计算成本。由于拉格朗日框架内存在函数二次破碎模型[8，9]，因此2352-7110/©2019作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx2A. Spitzenberger，S.诺伊曼，M。Heinrich等人粤公网安备44010802000011号+−阿 （）+·=ρρu1拉格朗日公式要求单元尺寸比颗粒直径大数倍，这与ICM相比显著降低了计算成本。但由于这些模型仅限于特定的液滴形状或大小特征，因此违反这些限制可能导致严重的不准确性。为了有效地检查整个原子化过程的数值，除其他外，欧拉和拉格朗日方法之间的过渡是必要的。过渡必须基于欧拉VOF模拟的相场内分离的流体粒子的检测。结合检测，评估关于颗粒的尺寸、位置和速度的信息。本文将这一过程称为欧拉粒子检测。然后，所获得的信息可以作为初始值转移到拉格朗日框架中的点粒子。Hebert等人[10]和Herrmann [11]已经做出了与欧拉粒子检测一致的第一次努力。 Hebert 等人[10]介绍了一种基于连接组件标记（CCL）的方法，用于在VOF模拟中自动识别液滴。欧拉粒子的检测是在模拟域的三维部分中进行的。几个研究小组也应用了类似的3D方法[12Grosshans等人[21]报告了一种略有不同的方法，他们在预定义的耦合层（2D）上检测颗粒。然后，可以将初次破碎后获得的液滴特性转移到随后的二次破碎的欧拉-拉格朗日模拟中本工作的目的是使用开源CFD库OpenFOAM v1906 [22]改进和数值验证VOF模拟中的欧拉粒子检测。它最初是在OpenFOAMv1612中引入的，基于2D方法。由于发现OpenFOAM提供的颗粒检测产生不正确的结果，因此对源代码进行了更改。提出了改进的检测算法，并给出了改进的结果为了进行数值验证，本文介绍了四个测试用例。这些测试用例在要检测的颗粒的形状和数量2. 软件描述2.1. 液体体积法本文采用流体体积法（VOF）[6在欧拉框架下对两种不混溶流体的界面追踪。因此，引入相分数α来区分两相：2.2. 粒子检测算法增强的2D检测算法基于OpenFOAM v1906的正式实现的extractEulerianParticles函数对象[22]。它使用任意放置的faceZone 作为评估层，并完全并行化。为了识别VOF模拟中的欧拉粒子及其大小、位置和速度，采用了一种基于连通分量标记（CCL）方法的CCL起源于图论，例如用于计算机视觉[24，25]。为了应用空间CCL方法，CFD网格可以被视为由边连接的节点组成的图CCL将图划分为子集（连接结构），其中子集内的所有节点都可以通过边彼此访问。在使用评估层进行颗粒检测的情况下，对于2D空间CCL，仅考虑满足α > αt的细胞面，其中αt是相分数阈值（图11）。①的人。这些单元面表示图形的节点。连接这些单元面的面边形成图形的边。 由于空间CCL的过程已经在许多不同的出版物[10，12，15，21]中描述，因此省略了详细的说明，并在图1中示意性地给出。1.一、在每个时间步长处执行空间CCL以识别与评估层相交的颗粒。然后，通过额外的时间CCL将不同时间步长的切片组合成粒子，这是源代码更改的主要主题。修改的目的是正确处理和保存两个连续时间步长的切片之间的所有连接，这是不可能的。这是最初实施的情况。如第3节所示，这导致了对欧拉粒子的检测的显著改进。关于增强的时间CCL的每个时间步长处的过程如下（图2）。（2）：(1) 生成当前时间步（ti）和前一时间步（ti1）的切片之间的所有连接。如果先前和当前时间步长的切片至少在一个面处重叠，则切片之间存在时间连接。不同时间的切片之间的这种连接的一个例子如图所示。二、因此，图随着每个时间步长不断扩展。(2) 找到形成粒子的切片的连接结构这是通过从当前时间步的切片开始并搜索与前一时间步的切片的连接来完成的。从潜在连接的先前切片，该算法再次搜索连接的当前切片etα={00<α1 1phase2，接口，阶段1（跟踪阶段）。（一）等等。如果没有更多的可访问切片，则该过程从尚未是任何相干结构的一部分如果所有当前切片均为它描述了某一相所占据的细胞体积，并指示了两相的空间分布。通过等温和不可压缩的质量和动量守恒方程（Navier-Stokes方程）以及相分数传输方程描述流动，n·u=0，（2）处理后，程序结束。(3) 相干结构的处理取决于它们包括的当前和先前切片的数量：(a) 一个当前切片，没有旧切片：生成一个新的欧拉粒子并与当前切片链接。(b) 一个或多个当前切片和一个先前切片：asso-1+将当前切片与属于旧切片的粒子相关联。(c)一 或 更 电流和 两 或多个先前ααu0，（4）阿勒特其中u是速度，t是时间，ρ是密度，p是压力，ν是运动粘度。符号g和f分别表示重力加速度和表面张力。根据Brackbill [23]引入的连续表面力（CSF）模型，在动量方程中表示表面张力σ。切片：将可能的多个部分组合在一起，属于旧切片的粒子，并将该组合粒子与相干结构的当前切片相关联。这方面的一个例子是图中的蓝色粒子。二、两个先前分离的相干结构直到ti才结合成一个粒子，因为当前切片1（ti）连接到先前切片1和2（ti-1）。A. Spitzenberger，S.诺伊曼，M。Heinrich等人粤公网安备44010802000011号3⏐⏐=··Fig. 1. 2D空间CCL：相分数场α（左）代表黄色和蓝色粒子。通过相分数阈值αt，作为颗粒的一部分的面被识别（中间）。通过CCL将标记的细胞分组为代表颗粒切片的相干结构（右）。图二. 时间CCL：两个粒子和不同时间的评估层的位置（左）。这些粒子包含一个图形，粒子切片作为节点，它们的时间连接作为边。简化的1D评估层表示从2D空间CCL获得的切片标签（右）。如果连续时间步的两个切片至少在一个面上重叠，则它们之间存在连接(d)无当前切片和一个或多个先前切片：与先前切片相关联的粒子被完全跟踪，并计算最终粒子属性。跟踪的欧拉粒子的粒子属性（指数p）计算如下：Vn，tαn，t Φn，tVp=∑ ∑<$Vn，t，（6）表1水和空气的流体性质水空调相分数α1 0密度ρ[kg/ m3]1000 1运动粘度ν [m2/s]1·10−61. 48·10−6表面张力σ[kg/m] 0.07具有特定的直径d0和位置x0。它的移动速度速度u就像周围的空气一样。表1总结了所有相关dp=1036Vp/π，（7）0两个阶段的流体性质和情况概述Xp=（∑∑n，t·xn，t）/Vp，（8）设置见表2。本出版物中的所有模拟均使用interFoam求解器[26，27]执行，该求解器基于t n上的VOF方法，是开源CFD库的一部分up=（∑∑Vn，t·un，t）/V p.（九）OpenFOAM v1906.使用OpenFOAM的frozenFlow选项，t n动量方程不求解以保持速度场固定指数n和t表示面指数和时间。每个时间步长通过表面的体积是Φn，t，Φn， t是通过表面的体积通量。与原算法相比，改进算法在计算相位分数αn，t时考虑了相位分数αn， t。通过对若干面和时间步长上的相应面值求和并在适用时用ΔVn，t对其加权来计算所识别颗粒的体积Vp、位置xp和速度up等效直径dp等于体积为Vp的球体的直径。此外，改进的颗粒检测现在允许同时进行多个评估层。代码的示意图可以在图中找到。3.第三章。3. 说明性实例创建了一个带有空气中水滴的数值测试案例，以验证检测算法（图1）。 4 a）。初始化液滴在模拟过程中。因此，表面张力被忽略，这导致相分数的纯平流。域因此，液滴通过出口离开域，并将通过入口重新进入一旦液滴通过评估层30次，模拟停止。其余边界定义为滑动墙。除了时间离散，这是一阶精度，所有的数值方案是二阶。此外，还引入了三个测试用例来测试检测算法的准确性（表2）。在其中一种情况下，液滴被等效直径为d0的H形颗粒（H形）取代（图4b）。由于frozenFlow选项而保留的这种不寻常的粒子形状考虑了识别复杂粒子的能力。H形颗粒产生了一个具有挑战性的情况，其中评估层不n4A. Spitzenberger，S.诺伊曼，M。Heinrich等人粤公网安备44010802000011号∇ ·∇=、∼、、、d0速度u0[m/s]（0 0−10）（0 0−10）（0 0−10）（0 10−10）表3数值方案。Term方案欧拉方程高斯线（ρφu）高斯线性迎风梯度（U）（φα）Gauss vanLeer（φrbα）高斯界面压缩（τij）高斯线性ηu高斯线性校正线性插值颗粒形状由于表面张力。最后，分析了一个测试案例（多个），其中共有12个不同尺寸和形状的颗粒（图1）。4c）。由于平行于评估层的附加速度分量，入口和出口不再是循环的，并且每个颗粒仅被跟踪一次。所有的模拟都是并行运行的，以测试其适用性的检测算法的雾化模拟，其中并行化是不可避免的，由于高计算成本。计算层由处理器边界细分，因此并行化直接影响粒子检测过程。所用的数值方案见表3。为了评估，跟踪的粒子属性xp和Vp（等式2）被设置为：（6），（8））通过初始值归一化如下：VV p/V0，10x=（xx p−xx0）2+（xy p−xy0）2 1。（十一）图三. 粒子检测算法的示意图。 程序的修改部分用灰色方框突出显示。的虚线框并不是本工作中所介绍的代码的一部分临时看到两个不同的切片，这两个切片必须由算法指定给同一粒子。另一个带有H形粒子的测试用例（H形 *）考虑了质量和动量方程（无frozenFlow），从而产生动态变化的由于评估层的z坐标是固定的，因此通过在平行于评估层的平面内各种测试用例的结果如图所示。 五、请注意，原始算法（左侧）和修改后的算法（右侧）的比例有很大不同。如图2的箱线图所示。 5 a表明，Open-FOAMv1906的原始实现无法正确检测球形液滴。方框的顶部和底部表示下四分位数和上四分位数，而条带表示中位数。须线覆盖其余数据，但其长度限制为四分位距（IQR）的两个须以外的数据跟踪的体积和位置与预期结果有很大差异。 就颗粒体积而言，盒子非常狭窄，围绕V0 ，离群值广泛分布（0 1，这在理论上是不可能的。此外，粒子位置也显示出非常低的差异。在使用H形颗粒结合表面张力影响的情况下（图1）。5d），跟踪的体积V不像纯平流那样散射。由于H形颗粒在评估层的几个通道内被转移成近似球形，从而减小了其表面积，因此数值扩散的影响较小与球形液滴相比，跟踪的体积不太准确，但其偏差再次小于1%，这揭示了增强的检测算法的良好性能。关于标准化位置x，在该测试用例中可以观察到增加的偏差。然而，这不是修改后的算法的故障，而是由收缩到液滴引起的实际移位总共有12个粒子的测试情况也产生了良好的结果（图5e）。由于在计算域中存在多个粒子，因此必须使用基于粒子生成程序的较不精确的V0测定。 因此，出现了相对较大的V值的分散。然而，跟踪的体积仍然在可接受的范围内，最大偏差约为1.5%。由于横向速度分量，只有x-分量用来计算x的值。同样，通过以下方式实现了良好的精度：x= 0。01.此外，颗粒速度和检测到的颗粒数量与所有测试用例的正确值一致4. 影响和结论在这项工作中，粒子检测VOF模拟与开源CFD库OpenFOAMv1906的改进和数值验证。由于最初的实现产生了6A. Spitzenberger，S.诺伊曼，M。Heinrich等人粤公网安备44010802000011号图五. 箱形图（a）使用OpenFOAM v1906原始检测算法的液滴测试用例，左侧为刻度，以及（b）-（e）所有测试用例，增强的检测算法，在右手边缩放不正确的结果，提出了修改版本。为了验证的目的，介绍了四个具有不同颗粒形状、颗粒数量和基础物理的测试用例。增强的检测算法能够在评估层正确地检测颗粒及其尺寸、位置和速度，所有测试案例与实际颗粒体积的偏差通常小于1%。改进后的算法为今后高效模拟整个雾化过程奠定了基础。它提供了一个耦合的VOF模拟的初级破碎与欧拉-拉格朗日模拟的二次破碎的可能性耦合可以基于直接[12，14，15]或统计[21]方法。当然，颗粒检测也可用于分析初级破碎的所得颗粒尺寸分布。为了改善雾化模拟中的检测结果，可以实施两种不同相分数阈值的组合。较低的阈值可以用于尽可能正确地确定粒子体积，而较高的阈值指示粒子的正确数量。因此，可以实现改进的颗粒检测，特别是在液滴不能很好地分辨或彼此靠近的困难条件下。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作致谢这项工作得到了Deutsche Forschungsgemein- schaft（DFG，德国研究基金会）- Projektnummer的支持第 799 章 . 这 一 支 持 得 到 了作 者 。 所 有 作 者 也 感 谢 SebastianBorrmann的协助。附录A. 补充数据与本文相关的补充材料可以在https://doi.org/10.1016/j.softx.2019.100382上找到。引用[1]Marmottant P，Villermaux E.关于喷雾形成。流体机械杂志2004;498：73-111.http://dx.doi.org/10.1017/S0022112003006529网站。[2]放大图片创作者：Gorokhovski M. 一次雾化建模。 流体机械年鉴2008;40：343-66。http://dx.doi.org/10.1146/annurev.fluid的网站。40.111406.102200。[3]Harlow FH，韦尔奇JE。 含自由表面粘性不可压缩流体时变流动的数值计算物理流体1965;8（12）：2182-9. http://dx.doi.org/10.1063/1.1761178网站。[4]Osher S，Sethian JA.以曲率相关速度传播的锋面：基于汉密尔顿-雅可比公式的 算 法 。 J Comput Phys 1988;79 ： 12-49. http ： //dx.doi.org/10.1016/0021-9991（88）90002-2.[5] Sussman M，Smereka P，Osher S.计算不可压缩两相流解的水平集方法JComput Phys 1994;114 ： 146-59. http://dx.doi.org/10.1006/jcph.1994.1155网站。[6] 希特CW，尼科尔斯BD。自由边界动力学的流体体积法。J Comput Phys1981;39：201-25. http://dx.doi.org/10.1016/0021-9991（81）90145-5.[7]Warncke K，Gepperth S，Sauer B，Sadiki A，Janicka J，Koch R，et al.空气喷射液膜初始破裂的实验和数值研究。Int J Multiph Flow 2017;91：208-24.http://dx.doi的网站。org/10.1016/j.ijmultiphaseflow.2016.12.010.[8]O'Rourke PJ，Amsden AA.喷雾液滴破碎数值计算的调整片法。SAE技术文件872089，1987，http://dx.doi.org/10。4271/872089。A. Spitzenberger，S.诺伊曼，M。Heinrich等人粤公网安备44010802000011号7[9]Tanner FX.非蒸发柴油喷雾的液体射流雾化和液滴破碎模拟。SAE技术文件970050，1997，http://dx.doi。org/10.4271/970050。[10]Hebert DA，Schmidt DP，Knaus DA，Phillips S，Magari PJ.非结构网格中的并行vof喷雾液滴识别。在ILASS美国，第21届年度会 议 在喷雾系统。佛 罗里 达 州奥兰多，2008年。[11]赫尔曼·M界面追踪的精细水平集网格方法。 In：Annu res brief. Stanford，CA：Center for Turbulence Research; 2005，p. 3比18[12]赫尔曼·M并行欧拉界面跟踪/拉格朗日点粒子多尺度耦合程序。J Comput Phys2010;229（3）：745-59. http://dx.doi.org/10.1016/j.jcp.2009.10.009网站。[13] Ling Y，Zaleski S，Scardovelli R.用拉格朗日点粒子模型表示的小液滴雾化的多尺度模拟IntJMul-tiphFlow2015;76：122http://dx.doi.org/10.1016/j.ijmultiphaseflow的网站。2015.07.002。[14] Tomar G， Fuster D， Zaleski S，Popinet S.初级雾化 的多尺度模拟。ComputFluids2010;39（10）：1864-74.http://dx.doi.org/10.1016/j.compfluid.2010.06.018。[15]瓦利埃河空化模拟-从大蒸汽结构到小气泡动力学 [论文]，隆德大学; 2013年。[16] 王文，王文，等.复杂几何条件下一次雾化的直接数值模拟和大涡模拟. AtSpray2013;23（11）：1001-48.http://dx.doi.org/10.1615/AtomizSpr.2013007679网站。[17]周L，夏J，申乔J，凯恩斯A，克鲁夫L，布莱克希尔H。喷雾过程混合多尺度模 拟 方 法 的 发 展 。 工 艺 仪 表 机 械 工 程 师D 2016;230 （ 3 ） ： 322-31 。http://dx.doi.org/10.1177/0954407015585687网站。[18]金D，Desjardins O，Herrmann M，Moin P.同轴气体流对圆形液体射流的初始破碎。在ILASS Am，2007年在芝加哥[19] Zuzio D，Estivalèzes J-L，DiPierro B.一种改进的多尺度欧拉-拉格朗日方法 用 于 雾 化 过 程 的 模 拟 。 ComputFluids2018;176 ：285http://dx.doi.org/10.1016/j.compfluid.2016.12.018网站。[20]Edelbauer W，Birkhold F，Rankel T，Pavlović Z，Kolar P.用流体体积法模拟adblue注射器的液体破裂，然后进行离线耦合拉格朗日粒子跟踪。ComputFluids 2017;157：294-311. http://dx.doi.org/10.1016/j.compfluid.2017.09.003网站。[21]Grosshans H，Szász R-Z，Fuchs L.一种有效的流体统计体积-拉格朗日粒子跟踪耦合方法的发展。Int J Numer Methods Fluids 2014;74（12）：898-918.http://dx.doi.org/10.1002/fld.3879网站。[22]Openfoam：开源CFD工具箱。2019年，http://www.openfoam。com/.[23]作者简介：王建民，王建民.模拟表面张力的连续介质方法。J Comput Phys1992;100（2）：335-54. http://dx.doi.org/10的网站。1016/0021-9991（92）90240-Y。[24]Alnuweiti HM，Prasanna VK.平行架构和算法用于图像组件标记。IEEETransPatternAnalMachIntell1992;14 （ 10 ） ： 1014-34.http://dx.doi.org/10.1109/34.159904网站。[25]Choudhary A，Thakur R.在共享与分布记忆体多处理机上的连通元件标示演算法评估。在Proc int conf parallel processBeverly Hills，CA中。1992年，页。362-5[26]Deshpande SS，Anumolu L，Trujillo MF. 评价两相流解算器interfoam性能 。 Comput Sci Discov 2012;5. 014016. http://dx.doi.org/10.1088/1749-4699/5/1/014016网站。[27] 放大图片作者：Klostermann J，Schaake K，Schwarze R.表面压缩VOF法对单个上升气泡的数值模拟。Int J Numer Methods Fluids 2013;71（8）：960-82. http://dx.doi.org/10.1002/fld.3692网站。