用于双曲守恒律的解耦通量限制高阶有限差分方法

《Journal of Computational Physics》：Decoupled flux limiting high order finite difference methods for hyperbolic conservation laws

　　摘要：为了解决求解标量守恒律的高阶数值方法中保持数值解上下界的问题，将基于通量校正的传输方法推广到一种保守型高阶有限差分方法，从而在保持高阶精度的同时确保数值解的上下界得到约束。该方法的优势在于推导概念简单，且可广泛应用于高阶差分法、有限体积法或不连续Galerkin方法中，同时仅需较弱的CFL约束。然而，这一优势也带来不便，即限制参数以不等式组的形式相互关联，使得其在多维问题及系统中的应用变得复杂。本文提出了一种新型的解耦通量限制方法，将高阶有限差分框架中的每个限制参数独立处理。同时，我们证明了该方式在较弱的CFL约束下既能保持全局界限，又能维持高阶精度。通过对标量问题的数值模拟，结果显示该方法在保持全局最大值和最小值的同时，仍能实现预设的精度阶数。此外，该解耦通量限制方法在实现上简单，且易于推广，应用于计算可压缩Euler方程解的高阶有限差分方法时，还能以较低的计算成本保持正密度和正压力。

引言：保持界限的高阶数值方法及其在处理具有自然界限的物理量问题上的应用已引起广泛关注。满足物理界限不仅能提升数值结果的准确性，还能在某些情况下增强模拟的稳定性，例如通过确保正密度和正压力来可靠地计算可压缩Euler方程。对于求解标量守恒律的数值方法，其保持界限的一般思路是通过多项式或通量限制，使更新后的数值满足ujm≤ujn+1≤ujM，其中ujm和ujM为合适的局部上下界，用以抑制高阶近似带来的异常振荡。即便使用二阶精度的方法来求解守恒律问题（1.1），这一要求依然存在。不过，保持界限并不一定意味着完全变差稳定性，这曾是早期通量限制方法的主要目标，相关研究可参考[1]、[2]及其他文献。作为数值通量限制技术之一，基于通量校正的传输型方法试图寻找一阶单调通量与高阶通量的凸组合，期望以此在较弱的时间步长约束下同时实现保持界限与高阶精度。这类方法的基础可见于[3]，并由[4]进一步改进。例如，对于标量守恒律，可行的局部界限可定义为ujm=mink∈{j?1,j,j+1}ukn，ujM=maxk∈{j?1,j,j+1}ukn；而采用一阶单调格式时，则定义为ujm=mink∈{j?1,j,j+1}u?kn+1，ujM=maxk∈{j?1,j,j+1}u?kn+1；对于双曲系统，甚至可从u?n+1中提取界限，如[4]中所述。这类方法在平滑极值点附近的精度最多为二阶。[5]、[6]提出的通量限制方法可视为将通量校正传输方法推广到标量守恒律及可压缩Euler系统的高阶格式中。对于标量守恒律，[6]中仅通过ujm=minxu0(x)，ujM=maxxu0(x)来保证界限，从而实现高阶精度。确定局部界限是高阶差分方法在不牺牲精度阶数的前提下保持界限所面临的技术挑战之一。另一个造成技术困难的因素是大多数显式格式采用基于Euler前进的时间近似。以基本平流问题为例，设f(u)=u，初始条件为u(x,0)=sinx。在均匀网格上且无空间近似误差的情况下（这正是最高阶近似的核心理念），有uj1=sin(xj)?Δtcos(xj)=1+Δt2sin(xj+α)，这意味着经过一次Euler前进迭代后，数值至少会超出最大值O(Δt2)，具体取决于初始值的选取方式。因此，若要在一个Euler前进步骤内以rth阶精度保持界限，就必须满足O(Δxr)≥O(Δt2)，也就是说空间近似误差占主导地位。r表示空间近似的预设精度阶数。如果不考虑这一点，无论空间近似有多精确，最终都会导致精度下降。因此，[6]中提出了放宽的中间界限，用于多级TVD Runge-Kutta时间积分，或者[7]中提出了对积分后数值通量的通量限制，以此应对这一技术挑战。当没有空间近似误差时，有限体积格式的类似代数计算结果为uˉj1=1+Δt2sin(Δx/2)Δx/2cos(xj+α)，不过其在最小值或最大值超出的程度上不如有限差分格式明确。但从泰勒展开来看，其幅值的首项形式为1+Δt22?Δx224。若要在单个Euler前进步骤内同时保持数值界限和高阶精度，就需要满足Δt22≤Δx224这一约束，这等价于无论采用何种限制技术，都需满足CFL≤112，这一结论可通过[6]中的完整截断误差分析得到验证。从这个角度看，基于积分通量限制器[7]的通量限制方法具有明显优势，因为它允许使用更大的CFL值进行计算。本文的工作正是基于这一优势展开的。当将这类保持界限的通量限制方法应用于多维问题或系统时，其缺点也十分明显[5]、[8]，相比[9]提出的用于高阶有限体积法和不连续Galerkin方法的多项式限制技术，它显得较为粗糙且不够精细。当保持界限的约束以一组非线性不等式的形式出现时，这一问题更为突出，例如在计算可压缩Euler方程时需要保持正压力。其主要原因在于，用于更新单个网格点数值的数值通量限制参数是相互关联的。本文提出并研究了这种解耦的通量限制方法，其目的是充分利用保持界限的通量限制方法的优势——该方法概念简单，且允许使用更大的CFL值以提高效率。同时，通过引入成本较低的“中间”通量来实现通量限制，无需在相互关联的线性或非线性不等式中确定限制参数。此前已有学者尝试过类似的方法，可作为临时解决方案，相关研究可见[10]、[11]、[12]及其他相关文献。其中关键的是选择具有数学上合理精度保持能力的中间通量或辅助值，而目前高阶有限差分框架中还缺乏这一机制。此外还需指出，[12]中关于有限体积框架的研究提供了解析层面的支持。本文结构如下：第2节简要回顾了应用于积分后数值通量的参数化保持界限通量限制方法[7]，该方法在空间和时间维度上都具有高阶特性。我们介绍了用于一维及多维标量守恒律的高阶有限差分方法的解耦型保持界限通量限制方法。第3节通过截断误差分析给出了主要的精度结果。第4节将解耦通量限制技术扩展到可压缩Euler问题的正密度和正压力保持上。第5节通过大量数值实验验证了带有解耦保持界限通量限制器的高阶有限差分方法的精度，同时还通过计算一些极端可压缩Euler问题给出了数值结果，以证明所提方法的稳定性。最后，第6节对全文进行总结，并阐述未来的研究计划。

章节片段：
高阶有限差分方法的保持界限通量限制器：为说明[6]、[7]中参数化通量限制技术的基本思想，我们考虑在均匀空间离散化条件下对(1.1)式进行保守型高阶有限差分近似。一维标量守恒律(1.1)的完全离散化保守型高阶近似通常可表示为ujn+1=ujn?λ(H^j+12?H^j?12)，其中的数值通量H^j±12是通过高阶空间重构和时间积分得到的。λ表示ΔtΔx。例如：

截断误差分析：以三阶有限差分方法为例，我们希望证明对于一维标量问题，有H^j+12?H?j+12=O(Δx3+Δt3)。如果我们能够证明f^j?12(H^j+12+H^j?12)=O(Δx3+Δt3)，这一结论就会很显然，因为ujn?2λ(H^j+12?f^j)的超出最大值的幅度最多为O(Δx3+Δt3)。因此，只需将积分后的数值通量H^j+12修正到O(Δx3+Δt3)的精度，就能确保整个近似(2.7)满足保持界限的要求。

可压缩Euler方程的正性保持通量限制：在本节中，我们探讨将解耦通量限制方法应用于计算可压缩Euler方程ut+f(u)x=0，其中u=(ρ,ρu,E)T，f(u)=(ρu,ρu2+p,(E+p)u)T，ρ为密度，u为速度，p为压力，E=12ρu2+pγ?1来自状态方程(EOS)，γ为比热容比（空气的γ值为1.4）。将标量守恒律的保持界限方法推广到可压缩Euler方程，以实现正密度和正压力的保持。

数值实验：在本节中，我们对新开发的通量限制方法在标量守恒律及双曲系统上的应用进行了测试，主要目的是验证所设计的精度阶数是否得到保持。虽然也进行了保持界限性的验证，但这一环节并非必需，因为该方法本身就为数值方法满足局部或全局界限提供了充分条件。本文仅关注全局界限的保持，因为局部界限更为……

结论：本文提出并研究了一种用于求解双曲守恒律的高阶有限差分方法的新型解耦保持界限通量限制技术。通过数值实验和截断误差分析，我们可以确定这种新的通量限制方法在保持数值解的上下界方面十分有效。将其推广到可压缩Euler方程后，不仅证明了其有效性，也显示了其实现的简便性。

CRediT作者贡献声明：郑福徐：方法设计、研究实施、资金获取、正式分析、数据整理、概念构思。

利益冲突声明：作者声明不存在任何可能影响本文研究的已知财务利益或个人关系。郑福徐