基于子单元流体体积界面重建的三维多材料任意拉格朗日-欧拉方法

《Computers & Fluids》：A three-dimensional multimaterial arbitrary Lagrangian-Eulerian method based on the subcell volume-of-fluid interface reconstruction

　　摘要：本文将Sha和Jia在《Comput. Math. Appl.》160卷第86-107页（2024年）提出的二维流体体积法扩展到三维情形。在此基础上，我们为六面体网格开发了一种分层多材料任意拉格朗日-欧拉方法。所提出的三维流体体积法通过将每个界面单元划分为八个子单元来运作，进而通过最小二乘法最小化这些子单元中参考流体的重构体积分数与给定体积分数之间的差异，从而确定线性界面的法向量。同时，通过确保整个界面单元中参考流体的体积分数符合给定值来确定平面常数。该三维流体体积法保留了二维版本的一个关键特性，即无需依赖相邻单元的信息。不过，与使用不重叠三角形子单元的二维方法不同，我们的三维方法中的子单元是相互重叠的。所开发的任意拉格朗日-欧拉方法采用拉格朗日加重映射的架构：在拉格朗日阶段，运用兼容的有限体积离散化方法，并结合基于边缘的人工粘性机制来捕捉冲击波，以及Tipton压力松弛模型来实现界面单元的闭合；在重映射阶段，则通过多面体细分与相交算法来完成处理。为了全面评估所提方法，我们进行了包括静态界面重建、动态平流以及多材料ALE模拟在内的多项数值试验。结果表明，三维流体体积法在静态测试中能达到二阶精度，其性能优于某些传统流体体积法，精度上可与流体矩方法相媲美。动态测试进一步证明了该方法能够准确捕捉发生大幅变形的界面。最后，多材料ALE测试显示，该流体体积法能够生成更平滑的界面，且产生的碎片数量少于流体矩方法。

引言：涉及多种材料的可压缩流体的数值建模主要基于两种不同的空间框架：拉格朗日框架和欧拉框架。拉格朗日方法采用随流体运动而变形的网格，其优势在于能够精确呈现不同材料之间的界面。但其缺点是在面对极端材料变形时稳定性较差。欧拉方法则使用固定网格，虽然在大变形情况下具有更好的稳定性，但由于固有的平流误差会导致材料界面模糊。为寻求一种更为优化的解决方案，Hirt等人率先提出了任意拉格朗日-欧拉方法，该方法允许用户自定义网格的运动方式，从而将拉格朗日方法的精确跟踪能力与欧拉网格的强稳定性有机结合。尽管经过数十年的发展，任意拉格朗日-欧拉方法已成为工业领域多材料模拟的核心技术，但其传统形式在面临强剪切作用时仍存在诸多问题。正是这一缺陷促使Peery等人开发出了多材料任意拉格朗日-欧拉方法，该方法最初是为解决惯性约束聚变过程中的严重剪切变形问题而设计的。该方法的创新之处在于将界面捕捉理念融入任意拉格朗日-欧拉框架之中，这种混合方法使得该方法在模拟网格严重缠绕，以及强剪切作用下界面的撕裂与合并等传统任意拉格朗日-欧拉方法难以处理的场景时具有出色的表现。关于这些发展的详细历史和技术背景，Barlow等人的综述文章是一份极具参考价值的现代文献，同时也与其他基础研究著作形成了互补关系。

界面捕捉方法是多材料任意拉格朗日-欧拉方法中的核心组成部分，它指的是在流体动力学模拟中通过数值方法追踪并捕捉流体界面的位置与运动轨迹。这一方法对于理解流体界面的动态行为、预测流体的混合与分离过程具有重要意义，主要包括流体体积法、流体矩方法、水平集方法以及界面元素追踪方法等多种类型。流体体积法通过体积分数函数来描述流体界面，该函数表示每个计算单元中参考流体的体积分数。大多数流体体积法会将实际流体界面近似为线性界面，而要确定特定界面单元中线性界面的法向量，就需要利用相邻单元的体积分数信息，同时该界面的具体位置则由该单元中参考流体的给定体积分数决定。流体矩方法则同时利用参考流体的体积分数和质心信息来重建界面，其目标是在满足界面单元中参考流体体积分数给定值的前提下，尽可能减小重构质心与给定质心之间的差异。与流体体积法不同，流体矩方法仅依赖当前计算单元的信息，无需借助相邻单元的数据，且其精度已被证明可达二阶。水平集方法则通过连续的水平集函数来表示流体界面，该函数在界面处值为零，在流体内部为正，外部为负，通过求解该函数的变化方程即可追踪界面的运动与变形情况。界面元素追踪方法则直接追踪界面元素如拉格朗日点或界面网格的运动与变形，该方法能够高精度地模拟界面运动，但计算成本通常较高，同时还需处理界面元素的再生与合并等问题。

相比而言，流体体积法具有较好的守恒性，能够处理复杂的界面问题，且计算效率较高，因此被广泛应用于众多领域。流体体积法可分为几何型、代数型以及几何代数混合型三类。几何型流体体积算法通常精度较高，但实现起来相对复杂，尤其是在非结构化网格上；代数型流体体积算法则相对简单，能够适应任意拓扑结构的单元，但精度较低且存在较大的数值耗散现象。双曲切线界面捕捉算法属于几何代数混合型流体体积算法，它无需进行显式的界面重建或额外的压缩处理，只需设定合适的陡度参数，即可消除界面处的数值扩散现象，并使过渡层的厚度得到良好控制。在后续研究中，该算法得到了进一步改进，能够用于多维模拟以及非结构化网格的计算。如前所述，传统流体体积法需要利用相邻单元的体积分数来估算界面单元中界面的法向量，在采用六面体网格的三维情况下，甚至需要依赖多达26个相邻单元的数据。近期，Sha和Jia在相关研究中提出了一种名为SVOF的新流体体积方法，该方法仅需利用界面单元内若干子单元中的参考流体体积分数信息即可重建界面，无需依赖相邻单元的数据。具体而言，线性界面的法向量是通过求解一个非线性优化问题来确定的，该问题的目标是最小化这些子单元中给定体积分数与重构体积分数之间差异的平方和。在相关研究中，还采用了改进的二次拟合黄金分割算法来寻找最优的近似界面，其中界面的截距计算方式与传统流体体积法类似，即通过确保每个界面单元中参考流体的体积分数符合给定值来确定。实际上，SVOF方法可以视为流体体积法与流体矩方法的结合体——它像流体体积法一样利用流体的零阶矩信息，又像流体矩方法一样不依赖相邻单元的数据。此外，该方法的精度与流体矩方法相当，值得进一步推广与应用。本文正是将SVOF方法扩展到三维场景，并基于此方法开发了分层多材料任意拉格朗日-欧拉方法。该三维流体体积法沿用了二维版本中确定线性界面法向量与平面常数的基本策略，依然具备无需依赖相邻单元信息的特性，即通过最小二乘法最小化界面单元内若干子单元中重构体积分数与给定体积分数之间的差异来确定法向量，同时通过整个界面单元中参考流体的给定体积分数来确定线性界面的精确位置。不过，通过一些新的改进措施，该方法在精度与稳定性方面都表现出更优异的性能。一方面，我们采用了名为动态插值BFGS算法的高效优化算法来寻找最优的近似界面，该算法通过显式的动态插值函数来近似隐式的目标函数，并在优化过程中用插值函数的解析偏导数替代原始隐式目标函数的数值偏导数，从而显著提升了计算效率与精度。另一方面，我们对界面单元内的子单元定义进行了改进，允许子单元之间存在重叠，而二维版本中的子单元则是不重叠的三角形结构。这样的定义方式使得三维流体体积法更具稳定性，因为子单元之间的关联性有助于避免出现一些非物理现象。基于这种三维流体体积法进行界面单元的重建，我们进一步开发出了分层多材料任意拉格朗日-欧拉方法，该方法属于拉格朗日加重映射类型。在拉格朗日阶段，我们采用了兼容的拉格朗日有限体积方法，同时运用边缘粘性机制来捕捉冲击波，再结合Tipton压力松弛模型来实现界面单元的热力学闭合。此外，还采用了相关研究中的多面体细分与相交算法来进行重映射操作。数值实验结果表明，SVOF方法能够达到二阶精度，而所开发的分层多材料任意拉格朗日-欧拉方法也表现出良好的稳定性。本文的其余部分结构如下：第二节简要回顾流体体积法界面重建技术，第三节详细介绍三维流体体积法，第四节基于所提出的SVOF方法阐述分层多材料任意拉格朗日-欧拉算法，第五节报告各项数值试验结果，包括两项静态测试、两项动态平流测试以及两项多材料ALE模拟试验，最后第六节给出总结意见。

片段摘录：
流体体积法概述：考虑一个界面多边形单元Ωj，其体积记为|Ωj|。在Ωj中指定一种材料作为参考材料，假设该参考材料占据了一个非平凡的开集ω，其体积可通过|ω|=∫ωdx来计算。则单元Ωj中参考材料的体积分数为fj(ω)=|ω||Ωj|。分隔参考材料与其他材料的线性界面可以用以下平面方程表示：n?x+d?=0，其中n为法向单位向量。
子单元的定义：本节介绍适用于一般六面体网格的SVOF方法的三维扩展版本。该方法基于界面单元内若干子单元中参考流体的体积分数来重建材料界面。考虑一个面可以是非平面的通用六面体单元Ωj，设子单元的总数为Ns，将这些子单元记为Ωj,k，k=1,…,Ns。假设整个单元及其各子单元中的参考流体体积分数已知。
分层多材料任意拉格朗日-欧拉方法：基于第三节中提出的新SVOF方法，我们开发了一种分层多材料任意拉格朗日-欧拉方法，属于拉格朗日加重映射类型。该方法是对三维多材料任意拉格朗日-欧拉方法的扩展，采用方向上不分裂的算法，在拉格朗日阶段，首先更新每个界面单元及其子单元中的物理变量与体积分数，从tn变为tn+1，然后再将该数值从tn+1时刻的拉格朗日网格重映射到tn+1时刻的新网格上。图8展示了该多材料任意拉格朗日-欧拉方法的流程图。
数值试验：本节通过两项静态测试、两项动态平流测试以及两项多材料ALE模拟测试，来评估所提出的SVOF方法以及基于该方法的多材料任意拉格朗日-欧拉方法的精度与稳定性。同时，这些测试中也提供了使用流体矩方法进行界面重建的数值结果作为对比。SVOF方法和MOF方法均采用[32]中提出的DIBFGS优化算法来确定法向量。我们在研究中使用了两种类型的网格。

结论
本文构建了一个用于三维多材料流动模拟的完整框架，其贡献主要体现在两方面：一种3D子单元VOF（SVOF）方法——该方法在2D方法的基础上进行扩展，通过八个子单元在每个界面单元内实现界面重构。通过最小二乘法结合八个子单元中的参考体积分数来优化界面法向量，同时通过确保整个计算域的体积守恒来确定平面常数。

作者贡献声明
邓旺：写作——审稿与编辑、写作——初稿撰写、验证、软件、方法论、研究、形式分析。贾祖鹏：写作——审稿与编辑、监督、软件、研究、资金获取、形式分析、概念设计。

未引用的参考文献
[57]、[58]

利益冲突声明
作者声明没有已知的可能影响本文研究成果的财务利益或个人关系。

致谢
本研究得到了中国国家自然科学基金（项目编号12272060）的支持。同时，也非常感谢匿名审稿人的宝贵意见。

邓旺|贾祖鹏
中国北京100094，应用物理与计算数学研究所