一种用于不可压缩纳维-斯托克斯方程的不连续伽辽金一致分裂方法

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：A discontinuous Galerkin consistent splitting method for the incompressible Navier–Stokes equations

　　多米尼克·T·斯蒂尔（Dominik T. Still）|娜塔莉亚·内布利什维利（Natalia Nebulishvili）|理查德·舒斯尼格（Richard Schussnig）|卡塔琳娜·科尔曼（Katharina Kormann）|马丁·克龙比希勒（Martin Kronbichler）
慕尼黑工业大学计算力学研究所，波尔兹曼街15号（Boltzmannstra?e 15），加兴（Garching b. München），D-85748，巴伐利亚州，德国

**摘要**
本文提出了刘（Liu）在2009年提出的连贯分裂方案（consistent splitting scheme）的不连续伽辽金离散化方法。该方法隐式地满足了无散度约束，消除了速度-压力兼容性条件，并摒弃了压力边界层。即使对于具有开放边界和牵引边界的情形，也施加了连贯的边界条件。因此，时间的准确性不再受分裂误差的限制。对于二阶空间导数，采用了对称内部惩罚伽辽金方法（symmetric interior penalty Galerkin method）。对流项以线性隐式方式处理，这既放宽了显式方案（explicit schemes）的CFL限制，又避免了完全隐式公式（fully implicit formulations）所需的非线性方程求解。为了改进质量守恒，将Leray投影（Leray projection）与散度和法向连续性惩罚项结合使用。通过为速度场的散度和对流算子的散度选择适当的通量，可以证明连贯的压力边界条件仅由加速度和粘性项贡献得出，适用于L2离散化。

每个时间步长内，由于该方案在速度场和压力场方面的解耦特性，导致一个压力泊松方程（pressure Poisson equation）和一个矢量值对流-扩散-反应方程（vector-valued convection–diffusion–reaction equation）。我们验证了该方法在空间和时间上的最优收敛率，并证明了其与高阶时间积分方案（higher-order time integration schemes）的兼容性。一系列数值实验，包括围绕圆柱体的二维流动基准问题和三维泰勒-格林涡旋问题（three-dimensional Taylor–Green vortex problem），验证了其适用于实际相关流动问题的能力。

**1. 引言**
不可压缩流体流动在科学和工程中起着重要作用，具有巨大的社会价值。计算流体动力学（Computational Fluid Dynamics）旨在通过计算机模拟提供近似解，已成为一个基础学科，并迅速应用于其他研究和技术领域。尽管在过去几十年取得了巨大进展，但目前可用的高性能高保真算法实现仍然对当今最大的计算集群构成挑战。数十年来，相关数值算法及其实现的准确性和吞吐量一直在不断提高。然而，由于对准确性、求解时间、数值稳健性或守恒特性和物理一致性需求的增加，仍迫切需要进一步改进。这些创新通常由能源、航空或燃烧领域中的高分辨率应用驱动，或者受到多查询应用（如不确定性量化、敏感性分析、模型降阶、逆问题分析和机器学习）需求的推动。

瞬态不可压缩流动的纳维-斯托克斯方程（Navier–Stokes equations）通常通过隐式积分器（implicit integrators）随时间推进，这是由于无散度或不可压缩性约束。通过无散度条件，压力充当拉格朗日乘数（Lagrange multiplier），强制质量守恒。这种结构要求对速度和压力的近似空间进行特定选择，以满足所谓的inf-sup（Ladyzhenskaya–Babu?ka–Brezzi）条件[1]、[2]。自Hood和Taylor[3]首次使用inf-sup稳定的速度-压力对以来，有限元方法（finite element method）已成为有限体积方法（finite volume method）的可行替代方案。关于inf-sup稳定性和合适的稳定化方法，参见[4]、[5]、[6]、[7]、[8]。

标准的时间离散化方法将控制速度场和压力场的混合问题视为一个大型微分-代数系统（large differential–algebraic system），导致具有鞍点结构（saddle-point structure）的线性系统。根据流动状态和问题参数，压力Schur补数（pressure Schur complement）可以通过压力质量算子和泊松算子（pressure mass and Poisson operators）来近似[9]；而对于对流占主导的情况，可以通过例如基于增广拉格朗日的方法（augmented–Lagrangian-based methods）[10]、[11]、[12]、[13]来处理。然而，这些方法可能导致速度-速度块（velocity-velocity block）条件不良，从而需要专门的多网格方法（multigrid methods），参见[14]。其他通常对雷诺数（Reynolds number）不太敏感的替代方法包括压力对流-扩散（pressure convection–diffusion）[15]、[16]和最小二乘交换子（least-squares commutator）[17]、[18]预处理器，这些方法在合理的时间步长和网格质量下表现良好。其中一些概念也已扩展到更复杂的物理现象，最简单的情况是包含非线性粘性（nonlinear viscosity），例如用于地幔对流[19]、多相流动[20]或广义牛顿流体[8]。最近的研究试图重新表述鞍点问题背后的偏微分方程，以简化线性代数方面[21]。

为了（表面上）避免与鞍点结构相关的复杂性，一类称为分裂（splitting）、投影（projection）或分数步长（fractional-step）的方案通过解耦速度和压力的控制方程来避免在每个时间步长求解整体系统。然而，请注意，时间分裂或投影方案可能仍然受inf-sup条件的限制，从而限制了允许的有限元配对数量，参见[6]、[22]。将耦合问题分解为单独的步骤显著简化了预处理器的设计，尽管代价是压力控制的泊松问题边界条件变得更加复杂。相关的例子包括[23]、[24]中描述的增量压力修正方法（incremental pressure-correction methods）的旋转形式，Karniadakis等人提出的双重分裂方案[25]（也参见[26]），以及Guermond和Shen[27]以及Liu[28]提出的连贯分裂方案（也参见Huang和Shen开发的相关方法[29]）。最近的发展方向包括向复杂流变定律（complex rheological laws）的非平凡扩展[30]、[31]、[32]，遵循Huang和Shen[29]思想的无条件稳定性超越二阶精确方案，以及多相流动[33]、[34]和流体-结构相互作用[35]。

本文重点讨论了刘[28]提出的连贯分裂方案的高阶精确不连续伽辽金（higher-order accurate discontinuous Galerkin, DG）空间离散化方法及其框架的全面测试。选择连贯分裂方案而非传统投影方案的动机是避免压力边界层，并在时间上实现更高阶的准确性，同时避免解决混合公式中出现的整体系统。与Liu[28]采用连续有限元的方法不同，DG离散化提供了高阶精度和对流主导情况下的改进稳健性[36]。接下来的内容将简要评论最接近的竞争方案。

首先，与Karniadakis等人提出的双重分裂方案[25]及其由Krank等人[37]和Fehn等人[22]、[38]提出的稳健DG离散化相比，本方案对压力泊松方程采用连贯边界条件，并根据数值证据实现了等阶C0-连续插值[7]、[39]。此外，它是迈向高阶时间步进方案无条件稳定性的重要一步，因为该方法可以采用Huang和Shen[29]的思想进行扩展。其次，Rosales等人[40]采用了连贯分裂方案，并使用了高阶隐式-显式龙格-库塔（higher-order implicit–explicit Runge–Kutta）方案。然而，他们引入了涡度ω=?×u作为额外的未知数来自然地强制所谓的电边界条件（electric boundary conditions），但由于我们希望强制不同的边界条件，因此可以避免这种复杂性，这将在下文讨论。第三，Li等人[41]提出了一种基于连贯分裂思想的高阶方法，用于变密度纳维-斯托克斯方程（variable-density Navier–Stokes equations）。他们采用DG方法对密度演化进行有界限制（bound-preserving limiter），并对压力和动量方程使用连续伽辽金方案（continuous Galerkin schemes），结合了完全显式和隐式-显式龙格-库塔时间步进器（fully explicit and implicit–explicit Runge–Kutta time steppers）。与他们的工作相比，我们不考虑可变粘性，但旨在通过DG离散化动量方程来提高对流效应的稳健性。

本文的结构如下：第2节介绍了连贯分裂方案背后的基本思想，引出了第3节中提出的时间离散化方法。第4节讨论了L2一致性的DG离散化。第5节批判性地研究了所设计数值方法的性质，包括第5.1节中对精确解的空间和时间收敛性研究，以及第5.2节中围绕圆柱体的实际相关流动基准问题和第5.3节中的三维泰勒-格林涡旋问题[43]、[44]。第6.2节给出了结论性意见。

**2.1. 不可压缩流动的纳维-斯托克斯方程**
我们考虑在空间域Ω?Rd（d=2或3）内控制不可压缩流体流动的纳维-斯托克斯方程。它们包括动量平衡方程（1）
?u?t + (u??)u ? νΔu + ?p = f
在Ω×(0,T)，
其中ν是运动粘度（kinematic viscosity），以及连续性方程（2）
??u = 0
在Ω×(0,T)。

系统用无散度初始条件（3）
u(x,t) = u0(x)
在Ω中初始化。

边界Γ=?Ω被划分为不重叠的狄利克雷边界（Dirichlet boundary）和诺伊曼边界（Neumann boundary），分别用ΓD和ΓN表示，使得Γ=ΓD∪ΓN且ΓD∩ΓN=0?。在狄利克雷边界和诺伊曼边界上，满足以下条件：
(4)
u = g?onΓD×[0,T],
(5)
(ν?u ? pI)n = h?onΓN×[0,T]。
这里，n表示单位外法向量（unit outward normal vector）。

**2.2. 连贯分裂**
采用混合有限元方法（mixed finite element method）对纳维-斯托克斯问题（1）–(5)进行标准隐式离散化会导致具有鞍点形式（saddle-point form）的非线性系统。连贯分裂方案提供了一个有吸引力的替代方案，它将耦合的速度-压力系统分离，从而得到压力的泊松问题（Poisson problem）和动量方程的对流-扩散问题（convection–diffusion problem）。正如刘[28]所提出的，连贯分裂方法用压力泊松方程（PPE）
?Δp = ??(u??u) ? ??f
在Ω×[0,T]
来替代不可压缩性约束，通过取方程（1）的散度并在适当的规则性假设下使用时间导数和粘性项的散度。

为了推导狄利克雷部分ΓD上PPE的边界条件，方程（1）与法向量n点乘，得到
(7)
?p?n = ??g?t ? (g??)u + f?n ? ν?×(?×u)?n?onΓD×[0,T]。
这里，我们使用向量恒等式Δu≡?(??u)??×(?×u)来以旋度形式表示粘性项，这对分裂方法有利[45]。

在诺伊曼边界上，压力gp可以通过将其投影到外法方向上来获得，参见[7]（8）
?gp?p = ν?u:n?n ? I ? h?onΓN×[0,T]。
或者，如果已知压力值，可以直接在边界上施加压力gp。请注意，这种选择也与最实际相关的场景一致，因为通常需要在开放边界（如出流）上规定压力值。有关此分裂方案的详细讨论，请参见[46]。

**3. 时间离散化**
在这项工作中，采用后向差分公式（BDF）进行时间积分。感兴趣的时间区间（0,T]被分解为大小为Δt的时间步长tn到tn+1。方案按照刘[28]的方法编写，而各个步骤的呈现类似于[31]。

**3.1. 连贯分裂方法**
作为该方法的第一种变体，我们首先考虑先求解动量方程，然后在右侧使用外推的压力项，之后使用压力泊松方程隐式计算更新后的压力。

对于动量方程中的速度时间导数，选择J阶的BDF方案。对于右侧的压力项，采用Jc阶的外推，其中Jc=J或Jc=J?1，以提高数值稳定性——这对于高阶情况尤为重要。时间tn+1处的速度un+1通过求解
(9)
γ0Δtun+1 + u???un+1 ? ν???un+1 = fn+1 + ∑i=1JαiΔtun+1?i ? ∑i=1Jcβi??pn+1?i
来计算，其中γ0、αi和βi（i=1,…,J）以及βi?（i=1,…,Jc）是合适的时间积分和外推常数（参见Hairer等人[47]）。外推的速度定义为
?u??∑i=1Jβiun+1?i。这种线性隐式的对流项形式放宽了严格的CFL条件，同时避免了完全显式处理对流时产生的非线性方程组，但所有选项都能产生有效的方案，更多关于隐式-显式方案的讨论参见[31]、[48]。

**3.2. 压力步骤**
利用速度un+1，通过（10）求解压力pn+1的PPE
?Δpn+1 = ??un+1??un+1 ? ??fn+1，
使用连贯边界条件（11）
?pn+1?n = ?gn+1??un+1?n ? ν?×?×un+1?n ? ?tgn+1?n+fn+1?n?onΓD。
如果边界条件?tgn+1的时间导数无法解析得到，采用J阶的BDF方案进行近似[22]。在诺伊曼边界ΓN上，可以直接应用方程（8）的分解得到
(12)
gpn+1 = ν?un+1:n?n ? I ? hn+1?n?onΓN，
或者如果已知可以直接应用gpn+1。

**3.3. Leray投影**
只有在使用刘[28]提出的连贯分裂方案时，才隐式地强制执行无散度条件，该方案基于PPE替换了连续性方程。为了进一步抑制速度场中的散度，我们采用Helmholtz–Leray分解
(13)
u = u? ? ??，
参见Liu[28]和Pacheco等人[31]的更详细讨论。在这种表示中，u?是无散度的（在每一点都无散度），而??是无旋的（irrotational）。这种分裂可以用来推导所谓的Leray投影，通过取方程（13）的散度
(14)
??u = ??u? ? ???????u = ?????。方程（14）在齐次狄利克雷和诺伊曼边界条件（15）下得到完整表述：n???=0在ΓD上，以及?=0在ΓN上。由此得到的速度场u?=u+??是完全无散度的。然而，请注意，u?并不满足原始问题的切向边界条件。这一点可以通过将亥姆霍兹-勒雷分解（13）与切向单位向量τ进行点积来轻易证明，即τ?u?=τ?u+τ???。同时，由于?在ΓD上的法向导数的边界条件（见方程（15）的存在，我们得到了正确的法向边界条件un?u?=n?u+n???=n?udue。这种困境迫使我们必须在以下两种情况中做出选择：（a）一个不严格无散度的速度场，但满足初始条件；或者（b）一个无散度的速度场u?，但不满足切向边界条件。遵循Liu [28]和Pacheco等人[31]的方法，我们仅对旧时间步的速度场应用勒雷投影，这类似于在经典投影方案中省略投影步骤的做法（参见Guermond等人[45]的方程（3.17）–（3.18））。无散度的速度场u?被用于方程（9）右侧的加速度项中，该方程可以重写为（16）γ0Δtun+1+u???un+1?ν???un+1=fn+1+∑i=1JαiΔtu?n+1?i?∑i=1Jcβi??pn+1?i=fn+1+∑i=1JαiΔt(un+1?i+??n+1?i)?∑i=1Jcβi??pn+1?i。首先求解包含勒雷修正的动量方程。利用新的速度场，PPE（10）可以得到压力pn+1。勒雷变量?通过（17）?Δ?n+1=??un+1来确定，其边界条件由方程（15）给出。这导致每个时间步都需要求解一个额外的泊松方程，其中?n+1将在下一个时间步中使用。

3.3. 带有修正压力变量的方案
带有勒雷修正的一致分裂方案包括方程（16）–（17）和PPE（10）。通过与标准方案相比增加的成本，可以通过定义一个修正的压力变量p?来消除。因此，我们从动量方程（16）开始，其中修正的加速度项为（18）γ0Δtun+1+u???un+1?ν???un+1=fn+1+∑i=1JαiΔtun+1?i+∑i=1JαiΔt??n+1?i?∑i=1Jcβi??pn+1?i=fn+1+∑i=1JαiΔtun+1?i??p?n+1。修正后的压力p?n+1不是单独计算的，而是通过（19）p?n+1=∑i=1Jcβi?pn+1?i?∑i=1JαiΔt?n+1?i来定义的。请注意，动量方程（18）与方程（16）相比没有变化。相反，勒雷修正项和外推压力项被合并为一个变量。为了推进一个时间步，首先通过求解一个修正的压力泊松方程来计算修正后的压力p?n+1，该方程是通过线性组合之前时间步的（瞬时）PPE（10）和方程（17）得到的。修正后的PPE表示为（20）?Δp?n+1=?∑i=1Jcβi?Δpn+1?i+∑i=1JαiΔtΔ?n+1?i=∑i=1Jcβi???(un+1?i??un+1?i?fn+1?i)?∑i=1JαiΔt??un+1?i，其中所有涉及的项都来自之前的时间步。因此，一个时间步包括以下两个步骤：
1. 求解修正后的PPE（20）以得到修正后的压力p?n+1。
2. 求解经过勒雷修正的动量方程（18）以得到新的速度un+1，这等同于方程（16）。
步骤的顺序与上述推导中的顺序相反，但每个时间步仍然保持PPE的隐式形式。修正后的压力变量p?n+1是一种一致的压力近似（参见[28]、[31]），因此可以直接使用。这消除了勒雷投影所带来的计算开销：在计算速度un+1之后，不需要分别求解两个泊松问题——一个用于压力pn+1，另一个用于勒雷变量?n+1——而是在求解动量步骤之前就计算了修正后的压力p?n+1，结合了两个泊松问题的贡献。方程（20）中的项∑i=1Jcβi?fn+1?i?fn+1用于外推强迫项，使得修正后的PPE最终可以写为（21）?Δp?n+1=∑i=1Jcβi???(un+1?i??un+1?i)???fn+1?∑i=1JαiΔt??un+1?i。在一致边界条件下，对流项再次以Jc阶进行外推，而对于粘性项，可以选择Jp阶的外推。合适的Jp值可以是J或J?1（参见[49]、[50]）。然后，修正后的PPE受到边界条件（22）?p?n+1?n=∑i=1Jcβi??gn+1?i??un+1?i+∑i=1Jpβi??ν?×?×un+1?i?n??tgn+1?n+fn+1?nonΓD的约束，其中合适的外推系数βi?，i=1,…,Jp（参见例如Hairer等人[47]）。选择不同的J、Jc和Jp阶数的效果将在表7中的数值示例中展示。由于方程（15）中给出的?的边界条件，方程（22）不受影响，保持不变。在诺伊曼边界上，压力g?pn+1由（23）g?pn+1=∑i=1Jhβi?ν?un+1?i:n?n?I?hn+1?nonΓN给出，相应的系数βi?，i=1,…,Jh同样不受方程（15）的影响。这种方法遵循Liu等人[49]的压力近似形式，他们证明了该方法对于J阶的BDF方案可以达到O(ΔtJ)的收敛阶数。与[49]相比的唯一区别在于，修正后的PPE和边界条件中对流项的处理采用了Jc阶，而不是J阶，从而提供了更多的灵活性。

3.4. 迭代方案
为了在强非线性情况下提高方法的准确性和鲁棒性，可以对方案进行轻微调整，通过迭代PPE和动量方程直到收敛。迭代通过求解方程（21）、（18）来获得速度u?的起始值来初始化。一次迭代包括求解压力步骤（24）?Δp?=??(u???)u????fn+1?∑i=1JαiΔt??un+1?i（边界条件来自方程（11）–（12）和动量步骤，将対流速度设置为上次迭代的速度u?=u?），然后求解（25）γ0Δtu?+u???u??ν???u?=fn+1+∑i=1JαiΔtun+1?i??p?。一旦达到收敛，压力和速度分别设置为pn+1=p?和un+1=u?。需要注意的是，这种提高的鲁棒性和准确性是以需要求解多个线性系统为代价的。每次迭代的这个额外成本是否值得取决于具体的应用场景，但大多数流行的时间分裂或投影方法通常为了性能原因都是无迭代的。因此，我们在第5.3节中主要介绍迭代变体，以便进一步了解这些方法，而不是因为它每次迭代的成本较高。不过，与标准的混合求解方法相比，无论如何都可以避免鞍点问题。

3.5. 概述
表1总结了通过在单个时间步内应用动量和压力方程及其相关边界条件的一致分裂方案。表1. 用于推进一个时间步的一致分裂方案版本：相应的动量和压力步骤以及相关的边界条件（BCs）。
| 方案 | 步骤 |
|------|----------------------------------------|
| 一致分裂 | 1. 动量方程（9） | 2. 带有边界条件的PPE（10）：（11）–（12） |
| | 勒雷投影 | |
| | 校正后的动量方程（16） | 2. 带有边界条件的PPE（10）：（11）–（12） |
| | 勒雷投影（17） | |
| 修正压力 | 1. 带有边界条件的修正PPE（21）：（22）–（23） | 2. 校正后的动量方程（18） |
| | 迭代方案 | |
| | 1. 带有边界条件的修正PPE（21）：（22）–（23） | 2. 校正后的动量方程（18） |
| | 迭代：（a） 带有边界条件的修正PPE（24）：（11）–（12） | （b） 校正后的动量方程（25） |

3.6. 与具有周期性边界条件的高阶对偶分裂的相似性
下面展示了带有修正压力变量的方案与Karniadakis等人[25]提出的高阶对偶分裂之间的相似性。对偶分裂方案由以下方程给出：（26）γ0u??∑i=0J?1(αiun?i)Δt=?∑i=0J?1βi(un?i??)un?i+fn+1，（27）?Δpn+1=?γ0Δt??u?，（28）u??=u??Δtγ0?pn+1，（29）γ0Δtun+1?νΔun+1=γ0Δtu??。将方程（26）代入（27）得到带有勒雷修正的PPE（方程（21）。进一步地，如果将方程（26）和（28）代入方程（29），它类似于动量方程（9）。不同之处在于对流项的处理方式，对偶分裂方案中是对流的显式处理，而这里提出的一致分裂方法是线性隐式的。然而，请注意，在对偶分裂方案中也可以考虑线性隐式处理（参见[51]及相关工作[26]、[32]、[45]、[52]）。对于在提出的不一致分裂方法中对对流算子的显式处理，两种方案仅在边界条件上有所不同，见附录7。在周期性边界条件下，两种方案是相同的，假设只考虑时间离散化，即问题在空间上保持连续。对于提出的L2-一致的空间离散化，方程（26）的右侧可以使用经典的通量边界条件，而PPE的右侧则包含了针对测试函数梯度的对流项（见方程（35）。

这种等价性引发了对这里提出的一致分裂方法中速度和压力函数空间选择的质疑。Karniadakis等人[25]的对偶分裂方法已知受到inf–sup条件的限制，而一致分裂方案（不使用勒雷投影）据报道不受此类条件的限制[28]、[49]、[53]。然而，为了获得最佳收敛率，采用勒雷投影可能会在某种意义上重新引入inf–sup条件。因此，使用等阶插值处理粘性流动时可能导致的潜在精度损失并不是选择ku=kp+1的唯一原因。最后，我们指出，所提出的一致分裂方案也适用于H1-一致或H(div)-一致的离散化，因为它避免了对方程（26）中质量矩阵系统的求解。第5.3节提供了当前分裂方案与Karniadakis等人[25]的对偶分裂方案的比较。

4.1. 压力泊松步骤
对于PPE（21），选择了对称内部惩罚伽辽金（SIPG）离散化[54]、[55]、[56]以避免引入辅助变量。原始公式为：寻找p?n+1∈Vhp，使得（31）lpp?n+1,q=rpui=1,…,Jn+1?i,fn+1,gi=1,…,Jn+1?i,q，对于所有q∈Vhp和所有元素Ωe∈Ωh。这里，左侧lp定义为（32）?lp(p?n+1,q)??p?n+1,?qΩe??q,12p?n+1n?Ωe?Γh?q,{{?p?n+1}}?n?Ωe?Γh+q,τp?n+1?Ωe?Γh??q,p??n+1n?Ωe∩ΓN?q,?p??n+1?n?Ωe∩ΓN+q,2τp??n+1?Ωe∩ΓN。右侧rp包括来自体力f的贡献、对流算子的项、时间导数、粘性项、勒雷投影和SIPG离散化的非齐次部分：（33）?rpui=1,…,Jn+1?i,fn+1,gi=1,…,Jn+1?i,q?rf(fn+1,q)+∑i=1Jcβi?rconv(un+1?i,gn+1?i,q)?+rSIPG(gpn+1,un+1?i,gn+1,q)?∑i=1JαiΔtrLeray(un+1?i,gn+1?i,q)。通过选择中心通量并将强迫项、对流算子和勒雷投影项分别积分一次，方程（22）中给出的一致边界条件简化为仅来自加速度项和粘性项的贡献。所有其余项在边界上相互抵消。强迫项写为（34）?rf(f,q)?(f,?q)Ωe?(f?n,q)?Ωe?Γh?f?n,q?Ωe∩ΓN。在狄利克雷边界上，强迫项?f?n,q?Ωe∩ΓD与一致边界条件中的强迫项f?n,q?Ωe∩ΓD相互抵消。对流项以对流公式表示，并且也分别积分一次。在狄利克雷边界上，平均算子的贡献是g??u，方程（22）中的对流贡献是?g??u，两者相互抵消，得到（35）?rconv(u,g,q)??u??u,?qΩe+q,{{u??u}}?n?Ωe?Γh+q,u?n?Ωe∩ΓN。勒雷投影项的散度算子得到边界条件（36）?rLeray(u,g,q)???q,uΩe+q,{{u}}?n?Ωe?Γh+q,u?n?Ωe∩ΓN。SIPG离散化的非齐次部分包括一致边界条件的其余部分（方程（22）），由（37）?rSIPG(gpn+1,un+1?i,gn+1,q)?q,2τgpn+1?Ωe∩ΓN??q,gpn+1n?Ωe∩ΓN?q,γ0Δtgn+1?n?Ωe∩ΓD?∑i=1Jpβi?q,ν?×?×un+1?i?n?Ωe∩ΓD给出，其中gpn+1来自方程（23）。注意，由于方程（15）中给出的?的边界条件，方程（22）不受影响，保持不变。在诺伊曼边界上，压力g?pn+1由（23）g?pn+1=∑i=1Jhβi?ν?un+1?i:n?n?I?hn+1?nonΓN给出，相应的系数βi?，i=1,…,Jh同样不受方程（15）的影响。这种方法遵循Liu等人[49]的压力近似形式，他们证明了该方法对于J阶的BDF方案可以达到O(ΔtJ)的收敛阶数。与[49]的唯一区别在于，修正后的PPE和边界条件中对流项的处理采用了Jc阶，而不是J阶，提供了更多的灵活性。

3.4. 迭代相反，我们遵循Krank等人[37]的方法，首先通过标准的逐元素L2投影计算涡度?ω??×u，然后再次计算旋度，以获得所需的?×ω≡?×?×u。4.2. 动量步骤 动量步骤包括找到un+1∈Vhu，使得对于所有v∈Vhu和所有元素Ωe∈Ωh，都有(38)luun+1,u?,gn+1,v=ruui=1,…,Jn+1?i,u?,p?n+1,gn+1,gpn+1,hun+1,v成立。左侧lu(un+1,u?,v)包括时间导数项m(u,v)、对流项c(u,u?,v)和粘性部分v(u,v)。为了稳定性，考虑了散度惩罚项adiv(u,u?,v)和强制法线连续性的惩罚项acont(u,u?,g,v)，从而得到(39)?luun+1,u?,gn+1,v?m(un+1,v)+c(un+1,u?,v)+v(un+1,v)+adiv(un+1,u?,v)+acont(un+1,u?,gn+1,v)。时间导数项表示为(40)?m(u,v)?γ0Δtu,vΩe。对流项c(u,u?,v)可以以对流形式cconv(u,u?,v)或散度形式cdiv(u,u?,v)来表示。对流项的对流公式使用上风通量(41)?cconv(u,u?,v)?u???u,vΩe?{{u?}}?nu,v?Ωe?Γh+{{u?}}?n{{u}}+12|{{u?}}?n|u,v?Ωe?Γh+v,|u??n|?u??nu?Ωe∩ΓD。在对流项的散度公式中采用了局部Lax–Friedrichs通量(42)?cdiv(u,u?,v)??u?u?,?vΩe+({{u?u?}}n)+12Λu,v?Ωe?Γh+v,|u??n|u?Ωe∩ΓD+v,(u?u?)n)?Ωe∩ΓN，其中Λ=ζLFmax(|u+??n|,|u???n|)，ζLF可以设置为0.5以允许更大的时间步长。注释：以斜对称形式表示对流算子并不像在连续Galerkin设置中那样具有相同的优势（例如，在能量稳定性方面）。这种差异源于在对流和散度公式中都出现的面贡献，即使选择了相同的通量。与压力算子类似，SIPG方法也用于粘性部分(43)?v(u,v)?ν?u,?vΩe?n??v,12νu?Ωe?Γh?v,ν{{n??u}}?τu?Ωe?Γh?n??v,νu?Ωe∩ΓD?v,νn??u?2ντu?Ωe∩ΓD。张量积有限元的稳定性参数τ再次根据Hillewaert [57]来选择。散度惩罚项定义为(44)?adiv(u,u?,v)?ζDhe‖u?ˉ‖ku+1??u,??vΩe，其中惩罚参数ζD=1，‖u?ˉ‖是逐元素体积平均的外推速度，he=Ve1/3是特征元素长度，Ve是元素体积[38]。注意，速度散度惩罚与连续元素的grad–div稳定性有关[59]，[60]。连续性惩罚项按照[37]，[38]，[58]在内部面上和Dirichlet边界上沿法线方向应用，定义为(45)?acont(u,u?,g,v)?v?n,ζC‖u?ˉ‖(u?n)?Ωe?Γh+v?n,2ζC‖u?ˉ‖(u?g)?n?Ωe∩ΓD，其中惩罚参数ζC=1。由于分裂方法具有统一的边界条件，惩罚项可以直接纳入动量方程。也可以像[38]中那样在后处理（或投影）步骤中包含这些项。注释：散度和连续性惩罚可以看作是强制解决方案空间类似于H(div)元素，前者调节速度场的无散度性，后者调节法线方向的连续性[59]。右侧包括强迫项f(f,v)、时间导数项mr(ui=1,…,Jn+1?i,v)、压力项bp(p,gp,v)、对流项cr(u?,g,v)和粘性项vr(g,hu,v)，表示为(46)?ru?f(f,v)+mr(ui=1,…,Jn+1?i,v)+bp(p?n+1,gpn+1,v)+cr(u?,gn+1,v)+vr(gn+1,hun+1,v)，其中(47)?f(f,v)?(fn+1,v)Ωe和BDF积分器的先前时间速度，(48)?mr(ui=1,…,Jn+1?i,v)?∑i=1JαiΔtun+1?i,vΩe。压力梯度项的公式是通过两次分部积分得到的(49)?bp(p,gp,v)???p,vΩe+12pn,v?Ωe?Γh+(p?gp)n,v?Ωe∩ΓN。尽管压力梯度项仅出现在动量平衡方程的右侧，但Fehn等人[22]的数值结果表明，为了实现预期的收敛率，确实需要分部积分。对流公式中对流项的非均匀贡献表示为(50)?cr,conv(u?,g,v)??v,u??n?|u??n|g?Ωe∩ΓD。对于散度公式，右侧为(51)?cr,div(u?,g,v)?v,|u??n|g?(g?u?)n?Ωe∩ΓD。粘性部分的贡献是由于SIPG离散化(52)?vr(g,hu,v)??n??v,νg?Ωe∩ΓD+v,2ντg?Ωe∩ΓD+hu,v?Ωe∩ΓN，其中hu=ν(?u)n。注意，结合方程(49)和(52)可以恢复完整的牵引力h=ν(?u)n?gpn，因此只需要在Neumann边界指定h。4.3. 时间步进算法 本节总结了算法1中修改后的压力方案的时间步进过程。为了简化介绍，这里省略了启动过程，该过程包括在所需时间点插值给定的精确解，使用逐渐增加阶数的低阶BDF方案或其他适当阶数的时间积分器，直到收集到必要的数据un+1?i。下载：下载高分辨率图像（164KB）下载：下载全尺寸图像5. 数值结果 本节检验了所提出方法的数值特性。首先，在第5.1节中展示了具有已知解析解的示例在空间和时间上的收敛性，随后评估了该方案在实际相关场景中准确预测流体力学量的能力，经典基准测试包括围绕圆柱体的二维流动[42]和Wang等人[44]提出的三维设置中的Taylor–Green涡旋[43]。基于矩阵和基于矩阵的自由实现都可通过BDFConsistentSplitting方案在ExaDG2 [61]中获得，该方案基于deal.II [62]。PPE使用通过多重网格算法预处理的共轭梯度求解器来解决，而动量方程使用GMRES求解器解决，因为得到的线性系统不是对称的。除非另有说明，否则相对求解器容忍度设置为10?6，绝对求解器容忍度设置为10?12，以未预处理的残差来衡量。除非另有说明，否则结果都是针对修改后的压力方案。5.1. 人造解 首先，我们关注二维Taylor–Green涡旋测试案例[43]，也参见Hesthaven和Warburton [36]。不可压缩Navier–Stokes问题的解由(53)u(x,t)=(?sin(2πx2)+sin(2πx1))exp(?4νπ2t)，p(x,t)=?cos(2πx1)cos(2πx2)exp(?8νπ2t)给出，右侧为零，f=0。考虑的域?Ω?[?0.5,0.5]2，边长为1的正方形，并在时间间隔0≤t≤T=1内进行模拟，粘度ν=0.025。Dirichlet和Neumann边界条件及初始条件是从精确解中得出的。速度和压力的相对L2误差计算为(54)?p=‖p?ph‖L2(Ω)‖p‖L2(Ω)|t=T和?u=‖u?uh‖L2(Ω)‖u‖L2(Ω)|t=T。5.1.1. 时间收敛性 为了验证预期的时间收敛率，网格均匀细化了四次，得到一个均匀的笛卡尔网格，he=1/32。对于速度使用多项式度数ku=5的四边形元素，对于压力使用kp=4的四边形元素。因此，空间离散化的误差比时间误差小。对于BDF-1和BDF-2，Jc和Jp选择等于J，而对于其他情况，使用Jp=Jc=J?1。图1显示了BDF-1到BDF-4的最佳收敛性，即O(ΔtJ)，对于速度和压力都是如此，即使对于J>2也是如此。在当前场景中，达到了空间误差的水平，但没有观察到J>2时收敛率的降低。然而，通常不能排除收敛率降低到Jc=J?1的情况。这种潜在行为的彻底分析超出了当前研究的范围。下载：下载高分辨率图像（320KB）下载：下载全尺寸图像图1. 通过一致分裂方案获得的2D Taylor–Green涡旋基准测试中速度（左）和压力（右）的时间收敛性，对于BDF-1到BDF-4。5.1.2. 空间收敛性 为了计算空间误差，将时间步长设置为6.1035×10?5，应用BDF-2方案，使得时间误差低于空间误差。图2显示了网格细化下的相对L2误差。对于所有多项式度数ku=2,…,5和kp=1,…,4，可以观察到hk+1的最佳收敛率。对于ku=2，线性求解器的绝对容忍度设置为10?8，以实现两个最细网格的适当收敛率。注释：添加散度和连续性惩罚项不会改变收敛阶数，也不会影响此特定示例中误差的大小。然而，忽略Leray投影——在当前的分裂方案形式中，它仅导致压力变量被修改——会降低方法的收敛阶数：压力只观察到hkp+1，而不是最佳收敛率hkp+1，而速度则观察到hku?1，而不是最佳hku+1。同样，没有Leray投影时的时间收敛率也受到影响。这表明添加Leray投影对于所提出方案的性能至关重要。实施工作和算子复杂性受到的影响很小，因为在时间步长期间不需要解决额外的线性系统。下载：下载高分辨率图像（308KB）下载：下载全尺寸图像图2. 通过一致分裂方案获得的2D Taylor–Green涡旋基准测试中速度（左）和压力（右）的空间收敛性，对于不同的多项式度数。5.1.3. 一致分裂变体的比较 表2和表8列出了第3.5节中提出的方案变体的相对L2误差和收敛率。所有变体都实现了最佳收敛率，除了Leray投影方案。在这种变体中，首先用外推压力解决动量方程，然后在加速度项中应用Leray校正（参见表1）。使用J=3或J=4会导致一些时间不稳定性，主要是对于后者。这略微降低了实现的时间收敛率。图3提供了BDF-3的这些数据的可视化，显示迭代版本的结果误差最低，因为它解决了选定的容忍度内的非线性。我们还提供了采用全牵引边界条件（23）的结果，表示为“Traction”。这需要在PPE的Neumann边界项中对速度进行外推，与修改后的压力方案相比，增加了误差。此外，当考虑BDF-3及更高阶数时的稳定性原因时，Jh=J减少到Jh=2。根据这些结果，修改后的压力版本似乎是所提出变体中的最佳选择，因为它结合了高计算效率和准确的结果。表2. BDF阶数1到4和各种一致分裂方案的相对速度L2误差和观察到的收敛率（括号内）。对于相应的压力误差和收敛率，请参见表8。方案Δt/T×10?2BDF-1BDF-2BDF-3BDF-4Leray3.133.46e?2 (1.15)1.48e?3 (2.33)9.15e?4 (2.43)5.24e?5 (3.42)projection1.651.58e?2 (1.13)3.01e?4 (2.29)1.64e?4 (2.48)4.88e?6 (3.42)0.787.30e?3 (1.11)6.37e?5 (2.24)2.90e?5 (2.50)4.31e?7 (3.50)Modified3.122.35e?2 (0.98)1.00e?3 (2.03)4.77e?5 (3.08)2.43e?6 (4.11)pressure1.561.18e?2 (0.99)2.47e?4 (2.02)5.79e?6 (3.04)1.54e?7 (3.98)0.785.90e?3 (1.00)6.11e?5 (2.01)7.18e?7 (3.01)1.31e?8 (3.55)Iterated3.121.18e?2 (0.99)2.47e?4 (2.02)5.78e?6 (3.04)1.46e?7 (4.05)1.565.90e?3 (1.00)6.11e?5 (2.01)7.10e?7 (3.03)9.40e?9 (3.95)0.782.94e?3 (1.00)1.52e?5 (2.01)8.76e?8 (3.02)9.42e?10 (3.32)Traction3.122.16e?2 (1.20)6.94e?4 (2.50)3.22e?4 (3.03)1.84e?5 (4.27)1.569.96e?3 (1.12)1.41e?4 (2.30)4.00e?5 (3.01)1.11e?6 (4.05)0.784.75e?3 (1.07)3.15e?5 (2.17)4.99e?6 (3.00)6.81e?8 (4.03)下载：下载高分辨率图像（328KB）下载：下载全尺寸图像图3. 使用BDF-3总结的一致分裂方案的不同变体的比较。由于需要减少外推阶数以确保稳定性，Leray校正方案显示出次优的收敛率。对于具有修改后压力变量的方案，没有观察到这种减少。5.2. 绕圆柱体的流动 使用Sch?fer等人[42]提出的2D-3案例研究了绕圆柱体的流动，目的是为了展示在实际相关设置中使用更高阶多项式度数的效果。考虑的几何形状是一个长度为L=2.2、高度为H=0.41的通道。在通道中放置了一个直径为D=0.1的圆柱体，见图4。在x1=0处的入口边界条件为(55)gu(x1=0,x2,t)=(4Umx2(H?x2)H2sin(πt/T)。平均速度为Um=1.5，粘度为ν=10?3，因此得到的最大雷诺数为Remax=100。时间间隔选为0≤t≤T=8。在x1=L处，给出一个开放的牵引（出口）边界，gp=0，而在所有其他边界处规定无滑移条件，u=0。5.2.1. 与参考值的比较 使用不同多项式度数的分裂方案获得的结果与文献[22]，[42]，[63]中的值进行了比较。文献[22]中的结果是使用ku=10和大约106个自由度最大阻力值cD,max（参考值2.95091839，来源于[22]，见表3）针对测试案例2D-3，使用显式、线性隐式和隐式对流项版本，在不同的CFL数值下进行模拟，采用了BDF-2时间积分方案。未能收敛的模拟结果用破折号标记。CFL数值：0.4 1.0 2.0 4.0 8.0 16.0

显式方法：2.95091536 2.95091496
线性隐式方法：2.95091557 2.95091625 2.95092445 2.95096155 2.95113375 2.95232103
隐式方法：2.95091523 2.95091335 2.95091206 2.95091085 2.95091590

下载：高分辨率图像（398KB）
下载：全尺寸图像

图4. 圆柱体周围流动的基准测试：经过2次均匀细化后的计算网格（上图）以及产生涡脱落现象时的速度矢量的水平分量ux（下图）。

5.2.2. 更高阶时间积分器和时间稳定性
为了展示更高阶时间积分与一致分裂方案结合的潜力，我们验证了各种对流项线性化变体的时间稳定性。表5展示了在不同目标CFL数值下，采用显式、线性隐式或完全隐式对流项公式时获得的结果，这些CFL数值用于控制时间步长。对于一致分裂方案，选择了对流项的显式公式，因为在这种示例中它允许更大的时间步长。结果显示，在具有800个元素的网格中，当ku=5时可以获得最大阻力系数。这种空间离散化在增加目标CFL数值时保持不变。显式方法仅在小CFL数值下保持稳定，而线性隐式公式至少在CFL数值为16时仍能保持稳健性。令人惊讶的是，完全隐式方法并未达到与线性隐式方案相同的CFL极限。然而，对于较大的CFL数值，隐式版本得到的结果更接近参考解，尽管需要求解非线性系统。

为了与双重分裂方案进行比较，表6展示了在不同CFL数值下，采用一致分裂方案和双重分裂方案（及其线性隐式公式）时获得的结果。对于双重分裂方案，连续性和散度惩罚项是在后处理步骤中添加的（参见[64]）。为了更好地比较，一致分裂方案也采用了相同的方法。这两种方案在使用BDF方案时都保持无条件稳定，直到包括二阶（即J=Jc=Jp=2），而对于BDF-3及以上阶数则会出现差异。对于一致分裂方案，PPE右侧对流项的外推阶数设置为Jc=J?1，这允许更大的CFL数值，而在双重分裂方案中该项的外推阶数为Jc=J。注意，在这种情况下，两种方案都使用Jp=J?1。此外，当CFL>1.0时，BDF-3相对于BDF-2的精度提升变得明显。

表6. 使用线性隐式一致分裂（CS）和双重分裂（DS）方案以及不同的BDF时间积分方案，针对测试案例2D-3获得的最大阻力值cD,max（参考值2.95091839，见表3）。未能收敛的模拟结果用破折号标记。

CFL数值：0.4 1.0 2.0 4.0 8.0 16.0
CS
BDF-2：2.95091557 2.95091625 2.95092445 2.95096155 2.95113375 2.95232103
BDF-3：2.95091538 2.95091365 2.95092369 2.95092490 2.95093743 2.95106423
BDF-4：2.95092587 2.95092424

下载：高分辨率图像（396KB）
下载：全尺寸图像

图5. 一致分裂方法在2100万个速度自由度（ku=7）下的动能（左图）、分子耗散（中图）和数值耗散（右图）与文献[38][65]中的结果进行比较。文献中的数据仅适用于t∈[0,10]的时间范围。较低的数值耗散意味着更精确的结果，但这受到固定时空分辨率的限制。

(i) 双重分裂方案：完全显式，CFL = 0.4，对流项采用[22]中的对流公式；
(ii) 耦合解方法：完全隐式，CFL = 0.4，对流项采用[22][32]等文献中的对流公式；
(iii) 一致分裂方案：线性隐式，CFL = 0.4，对流项采用本工作中的对流公式。

结果表明，所有这些方法都产生了可比较的结果，仅在使用相同空间分辨率的情况下显示出微小差异。作为参考，Fehn等人的参考结果[38]再次包含在图6中。这三种方案，即双重分裂方案、耦合解方法以及分裂方案，都在ExaDG [61][64]中实现。

下载：高分辨率图像（418KB）
下载：全尺寸图像

图6. 耦合解方法、双重分裂方案和具有786432个速度自由度的一致分裂方法在动能（左图）、分子耗散（中图）和数值耗散（右图）方面的比较。所有方案的结果都与Fehn等人的参考结果[38]相当，后者使用了32.21亿个速度自由度。较低的数值耗散意味着更精确的结果，但这受到固定时空分辨率的限制。

5.3.3. 更高阶时间积分器和时间稳定性
为了展示更高阶时间积分器对数值耗散的影响，我们固定空间离散化和目标CFL数值为1.6（BDF-4的稳定性极限），并使用线性隐式变体增加时间积分阶数。同样，我们选择Jc=Jp=J for BDF-1和BDF-2，而对于BDF-3和BDF-4，则选择Jc=Jp=J?1。图7中的结果显示，BDF-1方案与其他方案相比表现出更高的数值耗散。BDF-2方案接近更高阶方案，但在耗散率?上仍存在差异。BDF-3和BDF-4方案在t>14.0时的耗散率?有所不同，否则它们的结果相当接近，表明空间误差占主导地位。BDF-2和更高阶方案在初始时间步长期间观察到的增加的数值耗散是由于初始化过程引起的，该过程涉及执行较低阶方法的单个时间步长，直到收集到所需的BDF启动数据。

最后，为了使用线性隐式对流项变体在不同CFL数值下比较各种方案，我们定义了总数值耗散能量为（58）∫t=0^20 ?dEdt??dt=?(E(t=20)?E(t=0))?∫t=0^20?dt。

为了评估对流项的外推阶数Jc的影响，我们比较了Jc=J?1和Jc=J情况下的一致分裂方法。为了完整性，也包括了使用双重分裂方案获得的结果。由于BDF-1显示出较大的耗散且BDF-2是无条件稳定的，因此不再进一步分析BDF-1。表7显示了BDF-2、BDF-3和BDF-4的总数值耗散能量，其中后者采用了迭代版本的一致分裂方案。

下载：高分辨率图像（449KB）
下载：全尺寸图像

图7. 在CFL数值为1.6时，使用不同BDF方案的一致分裂方法（CS）的总数值耗散能量（左图）、动能（中图）和分子耗散（右图）。较低的数值耗散意味着更精确的结果，但这受到固定时空分辨率的限制。再次强调，所有方案都在后处理步骤中应用了惩罚项。选择Jc=Jp=J?1时，使用BDF-3的一致分裂方法在更高的CFL数值下保持稳定，与Jc=J和Jp=J?1的一致分裂方案和双重分裂方案相比，产生了相似的质量结果。结果表明，选择Jc=J?1而不是Jc=J可以提高一致性方案的稳健性，而不影响精度。BDF-2与BDF-3的比较再次证明了更高阶方法的准确性改进。对于迭代方法（CS?，Jc=Jp=J），与无迭代版本CS?和双重分裂方案相比，可以观察到更高的稳定性。BDF-4的CFL数值可以增加到2.0，而精度仍与无迭代版本CS相当。请注意，在这里无迭代一致分裂方案中压力的外推阶数降低了，这进一步证明了使用更高效的无迭代方法的必要性。

总之，本节的结果表明：(i) 该方案在空间和时间上具有更高阶精度；(ii) 对流项的隐式和线性隐式变体允许在实际示例中使用远大于1的目标CFL数值；(iii) 使用BDF-3和更高阶方法时，即使对于非平滑解，也可以通过使用更大的时间步长来提高精度，如在3D Taylor–Green涡旋基准测试中所展示的；(iv) 使用Jc=Jp=J?1的一致分裂方案与双重分裂方案相比，显示出更高的稳定性，并且与耦合解方法具有相似的精度。

表7. 使用线性隐式一致分裂方案（CS，当使用Jc=J?1时）和双重分裂（DS）方法以及不同的BDF时间积分方案获得的总数值耗散能量。较低的数值耗散意味着更精确的结果，但这受到固定时空分辨率的限制。未能收敛的模拟结果用破折号标记。

CFL数值：0.4 1.0 2.0 4.0
CS
BDF-2：0.0205 10.0235 10.0287 30.0346
BDF-3：0.0198 30.0215 20.0193 30.01255
CS?：BDF-3：0.0198 50.02179
––
CS?：BDF-3：0.0216 40.0246 50.02624
––
BDF-4：0.0213 50.0244 90.01949

6. 结论性评论
本研究提出了刘[28]提出的一致分裂方案的不连续Galerkin空间离散化方法，重点在于L2一致设置下的通量项公式以及一致边界条件的整合。对流项以线性隐式方式公式化，克服了显式公式严格的CFL条件，同时避免了完全隐式公式会导致求解非线性系统的问题。通过为对流算子的散度选择适当的通量，以及对PPE的速度散度进行积分，一致压力边界条件简化为仅来自加速度和粘性项的贡献。与Liu等人[49]相比，虽然降低了项的外推阶数，但不影响时间收敛率（如5.1节所示），但使用更高阶BDF方案时稳定性得到改善，Karniadakis等人[25]的相关方案也观察到了这一点。为了稳定性，添加了连续性惩罚项和散度惩罚项，同时使用Leray投影来改进质量守恒。与连续Galerkin离散化相比，使用高斯定理（对每个元素单独应用）重写curl–curl项，并用一阶导数替换二阶导数（类似于[28]或[21]），可能会导致观察到的收敛率降低。因此，对于高粘性问题，如果curl–curl项无法充分解决，使用当前方案的最低等阶离散化可能会导致精度损失。由于Leray投影要求使用inf–sup稳定对，这些缺点似乎可以忽略不计。从性能角度来看，压力Poisson方程通常是计算瓶颈，因此通常更倾向于使用inf–sup稳定方法中的更丰富的速度近似。

尽管存在这些小限制，所提出的工作提供了两个主要优势：(i) 与具有鞍点线性系统的耦合解方法相比，速度和压力是解耦的，导致问题更小，可以使用标准工具进行有效预处理；(ii) 与标准投影方法相比，避免了分裂误差。分裂误差常常导致虚假的压力边界层的产生，这些边界层是由于人为设定的边界条件而引入的，目的是为了使压力泊松问题得以求解。此外，这也使得BDF-3和BDF-4成为提高精度的可行选择，尽管这些时间积分器由于缺乏A稳定性而存在一定的局限性。本研究探讨的变体在具有已知解析解的示例中展示了更高的空间和时间收敛率，从而从数值上验证了提出的主张。进一步地，我们在实际相关的基准测试中（如圆柱体绕流问题和泰勒-格林涡流问题，参见Sch?fer等人[42]、Taylor和Green[43]以及Wang等人[44]）显示出了与最先进不可压缩流求解器的出色一致性。在后一个示例中，该方案与其最接近的竞争方案——Karniadakis等人[25]提出的双分裂方案进行了比较，两者性能非常相似，但重要的区别在于：（i）当前方案不受分裂误差的影响；（ii）最近的研究表明，BDF-3甚至更高阶的分裂方案在时间上也是无条件稳定的[29]——这一性质此前仅在其他分裂或投影方案中得到了二阶证明，但可能需要更高的误差常数。因此，一致分裂方案是求解不可压缩流问题的理想选择。从这个意义上说，本研究是朝着成功实现这一方法族的L2-一致不连续Galerkin离散化方向迈出的一步。

为了避免需要额外的稳定化项（尤其是连续性稳定化项），一致分裂方案可以扩展到适用于速度变量的H(div)一致元素，同时保持L2离散化作为PPE的自然假设空间。这就引出了关于这些方法性能的问题，需要考虑近似质量和求解时间，而在对流主导场景中对动量方程进行预处理的研究进展在此过程中起着关键作用。

**作者贡献声明：**
Dominik T. Still：撰写 – 审稿与编辑，撰写 – 原始草稿，方法论，概念化。
Natalia Nebulishvili：撰写 – 审稿与编辑，撰写 – 原始草稿，方法论，概念化。
Richard Schussnig：撰写 – 审稿与编辑，撰写 – 原始草稿，方法论，概念化。
Katharina Kormann：撰写 – 审稿与编辑，方法论，资金获取，概念化。
Martin Kronbichler：撰写 – 审稿与编辑，方法论，资金获取，概念化。