一种适用于可压缩流体-固体多相扩散界面模型的稳健能量-对流迎风格式和分裂压力方法

《International Journal of Multiphase Flow》：A robust energy-convective upwind and split pressure method for compressible fluid–solid multiphase diffuse interface model

　　•开发了一种E-CUSP方案，适用于可压缩流体-固体扩散界面模型。
•对流项的分割完全考虑了拉伸张量。
•自适应耗散可以抑制振荡并保持材料界面的完整性。
•能够清晰地捕捉固体和流体中的激波和膨胀波。
•该方案能够稳健地解决复杂的流体-固体相互作用问题。

引言
可压缩流体-固体多相问题的数值模拟是一个极具挑战性的课题。这些问题涉及可压缩流体中的激波和膨胀波、固体变形以及固体中的纵向和横向激波，同时还涉及材料界面。传统上，流体和固体状态通常采用两种不同的方法（欧拉描述和拉格朗日描述）来描述，因为它们的行为有所不同，这也会增加模拟的复杂性。
基于连续介质力学理论，Godunov和Romenskii在1972年及2003年建立了一个热力学上一致的保护性欧拉非线性弹性模型，在该模型中，粒子的变形和旋转可以通过一个非对称畸变张量来表征。这种模型具有双曲特性，并促进了各种流体-固体相互作用数值方法的开发，这些方法大致可以分为尖锐界面方法（SIM）和扩散界面方法（DIM）。虽然SIM方法能够保持尖锐的界面，但它们往往不具备守恒性，且难以扩展到三维问题。相比之下，DIM方法允许界面在多个网格单元上扩散，从而确保守恒性，并便于将用于单相流的激波捕捉方案扩展到三维问题中。
基于Godunov和Romenskii（1972, 2003）以及Favrie等人（2009）的工作，推导出了一个适用于极端变形下弹性固体-可压缩流体相互作用的欧拉扩散界面模型。在小变形极限情况下，Favrie等人给出了雅可比矩阵的特征值（五种波），并通过应用三波近似改进了Harten-Lax-van Leer接触（HLLC）求解器（Toro, 1999）来求解该模型。考虑到许多固体材料可以用特定的状态方程框架来描述，在零强度极限下该框架可简化为Mie-Gruneisen状态方程，Barton（2019）提出了一种新的扩散界面模型，该模型适用于涉及N个组分的问题，但只需要演化一个运动学张量。Barton（2019）讨论了该模型的特征结构（包括七种不同的波），然后扩展了HLLD求解器来求解它。在该求解器中，通过将两对横向波近似为单个横向激波来简化黎曼问题。HLLD求解器最初由Miyoshi和Kusano（2005）提出，用于求解理想磁流体动力学问题，他们考虑了由一个熵波和两个阿尔芬波分隔的四态黎曼扇形结构。随后，Ortega等人（2014）通过改进的幽灵流体方法将HLLD求解器扩展到了弹塑性流动问题，假设在快波和界面之间只存在一个慢波。对于用于求解弹塑性模型的HLLD求解器（Barton, 2019），除了黎曼扇形的近似（即五波近似）外，还估算了声速而不是直接计算声学张量，并假设剪切波速是恒定的。正如Barton所指出的，根据初始状态的选择，非线性弹性的黎曼问题解可能包括由七种不同波分隔的八个常数状态。这些简化有助于数值计算，但可能会产生一些不希望出现的数值问题，如振荡和超调。
与HLL型方案（如HLL、HLLC、HLLD）不同，CUSP型（对流迎风和压力分离）方案（Zha等人，2011；Liou和Steffen，1993；Liou，1996）首先将通量分为对流通量和压力通量（以下简称压力通量），然后对流项采用迎风格式处理，而压力通量则通过分离处理来实现。根据对流通量的组成成分，CUSP方案可以分为H-CUSP、E-CUSP和K-CUSP（H、E、K分别表示在能量方程中分离的对流通量包含焓、总能量和动能）。例如，AUSM（对流迎风分离方法）系列方案属于H-CUSP组。对于可压缩流体，E-CUSP方法中的压力分离项的特征值仅与声速有关，使其更适用于描述亚音速区域的声波传播，并具有更清晰的物理意义。Shen等人（2012）将E-CUSP方案扩展到了理想磁流体动力学方程中，这些方程具有与固体-流体欧拉模型相似的波成分（都包括七种不同的波）。Rodriguez等人（2019）为线性粘弹性（假设变形可以忽略不计，因此适用线性粘弹性）多相流开发了AUSM算法。
在这项工作中，我们基于Barton的可压缩流体-固体扩散界面模型开发了一种E-CUSP方案。在新方案中，通过将拉伸通量的第二项（Bijk=?(V?ijeuk?V?kjeui)/?xk）放在方程的右侧或作为另一个通量，并采用不同的离散化方法进行离散化，从而保持了拉伸张量的严格对流。因此，在对流通量中不会出现零通量Bkjk。此外，在对流项中自适应添加数值耗散和界面切换函数，可以帮助抑制固体激波附近的非物理振荡，并防止材料界面附近的过度数值耗散。文章的结构如下：第2节简要介绍了Barton模型的控制方程；第3节描述了HLLD求解器和我们E-CUSP方案的数值算法；第4节通过几个数值示例验证了新方法的准确性和稳健性；第5节给出了结论。

控制方程
Barton（2019）的流体-固体模型包括各组分的质量连续性守恒定律、动量守恒定律、能量守恒定律、除一个组分外所有组分体积分数的对流方程以及拉伸张量的非守恒方程。本文以二维模型为例来说明数值方法。
?q?t + ?f?x + ?g?y = ?se，
其中 q = ρ(l)?(l)ρu1ρu2ρu3ρE?(l)
V?ijef = u1ρ(l)?(l)u1ρu1?σ11u1ρu2?σ21u1ρu3?σ31u1ρE?uk
σk1u1?(l)u1V?ije?ui
V?1jeg = u2ρ(l)?(l)u2ρu1?σ12u2ρu2?

数值算法
二维控制方程（1）的半离散形式可以表示为：
?q?t + 1Δxfi + 1/2,j ? fi + 1/2,j + 1Δygi,j + 1/2 ? gi,j ? 1/2 = ?S
其中 fi+1/2,j 和 gi,j+1/2 分别是界面 (i+1/2,j) 和 (i,j+1/2) 处的数值通量。

数值示例
为了评估新方案的性能，解决了几个一维和二维的多相初始值问题。如Barton（2019）所述，使用带有van Leer限制器（van Leer, 1979；Sweby, 1984）的二阶MUSCL方法来计算原始变量的重构。采用三阶总变差减少（Total Variation Diminishing, TVD）Runge-Kutta时间积分方法（Shu和Osher, 1989），CFL（Courant-Friedrichs-Lewy）数为0.4。对于一维情况，时间步长为...

结论
构建了一种E-CUSP方案来求解可压缩流体-固体扩散界面模型。首先，将通量分为三项：第一项（对流通量）具有相同的特征对流速度u，第二项（广义压力通量）包含总应力，与拉伸张量相关的其他项放入第三项中。然后，为了控制数值耗散，引入了一个附加项的对流迎风方法。

作者贡献声明
Yi Cheng：撰写——审稿与编辑、初稿撰写、软件开发、调研、形式分析、数据整理、概念化。
Yiqing Shen：撰写——审稿与编辑、资源获取、方法论研究、资金争取。
Li Li：撰写——审稿与编辑、软件开发、方法论研究。

利益冲突声明
作者声明他们没有已知的财务利益冲突或个人关系可能影响本文所述的工作。

致谢
本研究工作得到了国家自然科学基金（NSFC）的支持，项目编号分别为12472299和12172364。
Yi Cheng | Yiqing Shen | Li Li