高阶 BDF 方法在高压雷诺数下不可压缩 Navier-Stokes 方程的稳健收敛性

《Communications in Nonlinear Science and Numerical Simulation》：Robust convergence of high-order BDF methods for the incompressible Navier-Stokes equations at high Reynolds numbers

　　

摘要

针对不可压缩流体的雷诺数鲁棒高阶后向差分格式（BDF）离散化方案的发展仍然具有挑战性，其主要原因在于二阶以上的BDF方法不具备A稳定性，因而不适用于标准的能量方法。本文针对高雷诺数下的不可压缩Navier-Stokes方程，构建了一组最高达五阶的隐式-显式BDF（IMEX-BDF）格式。我们发展了一个统一的全球能量框架，为这些高阶格式建立了鲁棒的L2范数误差估计。其核心要素是一类离散正交卷积核，它们将非A稳定的BDF格式重构为等价的低阶三角卷积系统，并实现了一个全局离散能量论证。我们证明了相关的收敛常数独立于雷诺数，从而使所提出的IMEX-BDF方法特别适用于高雷诺数不可压缩流。统一的分析还进一步实现了稳定且高效的自适应变阶策略的构造，该策略在高阶BDF离散化中平衡了稳定性与精度。针对高达$Re=10^8$雷诺数的广泛数值实验，包括圆柱绕流问题，证实了所提出格式的鲁棒性、精度和效率。

引言

考虑区域 $\Omega \subset \mathbb{R}^2$ 上的二维不可压缩Navier-Stokes方程（NSEs）：

$$\begin{cases} \partial_t \mathbf{u} - \nu \Delta \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} + \nabla p = 0, & \mathbf{x} \in \Omega \times (0, T], \\ \nabla \cdot \mathbf{u} = 0, & \mathbf{x} \in \Omega \times (0, T], \end{cases}$$

其中施加初始条件 $\mathbf{u}(\mathbf{x}, 0) = \mathbf{u}_0(\mathbf{x})$，$\mathbf{x} \in \Omega$，以及在 $\partial\Omega$ 上的无滑移边界条件。这里，$T > 0$ 表示终止时间，$\mathbf{u} = (u_1, u_2)^\top$ 表示流体速度，$p$ 为其压力。运动粘度系数为 $\nu = 1/Re$，其中雷诺数 $Re \gg 1$。设动能为 $E[\mathbf{u}] = \int_\Omega \frac{1}{2}|\mathbf{u}|^2 \, d\mathbf{x}$，线性动量为 $M[\mathbf{u}] = \int_\Omega \mathbf{u} \, d\mathbf{x}$，角动量为 $M_{\mathbf{x}}[\mathbf{u}] = \int_\Omega \mathbf{u} \times \mathbf{x} \, d\mathbf{x}$。在边界贡献消失的边界条件下，NSEs (1.1) 具有以下耗散和守恒性质：

- 动能耗散律：$\frac{d}{dt} E[\mathbf{u}] + \nu \|\nabla \mathbf{u}\|^2 = 0$，对 $t > 0$ 成立。
- 线性动量守恒律：$\frac{d}{dt} M[\mathbf{u}] = 0$，对 $t > 0$ 成立。
- 角动量守恒律：$\frac{d}{dt} M_{\mathbf{x}}[\mathbf{u}] = 0$，对 $t > 0$ 成立。

NSEs在流体力学中起着核心作用，并为准确高效的数值方法设计持续带来重大挑战，特别是在高雷诺数区域。一个关键困难源于不可压缩性、非线性对流和长时间稳定性之间的相互作用。标准的Galerkin离散化方法，如Taylor-Hood有限元，在离散水平上并不能精确保持散度自由约束，这可能导致物理不变量（如能量耗散和动量守恒）在长时间模拟中退化 [1]。为应对这些问题，结构保持格式已被广泛研究。特别地，EMAC（能量、动量和角动量守恒）格式已被证明即使仅施加弱散度自由约束（即对所有 $(\mathbf{u}, q) \in X \times Q$，满足 $(\nabla \cdot \mathbf{u}, q) = 0$，其中 $X := H^1_0(\Omega)$，$Q := L^2_0(\Omega)$），也能保持基本的物理定律。这一特性使得基于EMAC的离散化方法在速度仅为弱散度自由时，对不可压缩流的长时间模拟具有吸引力，参见 [1], [2], [3]。

设修正压力定义为 $p_{\text{emac}} := p - \frac{1}{2}|\mathbf{u}|^2$，令 $\nabla \mathbf{u}$ 的对称部分定义为 $\mathbb{D}(\mathbf{u}) := \frac{1}{2}(\nabla \mathbf{u} + (\nabla \mathbf{u})^\top)$。动量方程 (1.1) 的EMAC形式 [1] 如下：

$$\begin{cases} \partial_t \mathbf{u} - \nu \Delta \mathbf{u} + \nabla p_{\text{emac}} + 2\mathbb{D}(\mathbf{u})\mathbf{u} + (\nabla \cdot \mathbf{u})\mathbf{u} = 0, & \mathbf{x} \in \Omega \times (0, T], \\ \nabla \cdot \mathbf{u} = 0, & \mathbf{x} \in \Omega \times (0, T]. \end{cases}$$

为简便起见，当不引起混淆时，本文后续中仍用符号 $p$ 表示 $p_{\text{emac}}$。设双线性形式为 $a_\nu(\mathbf{u}, \mathbf{v}) := \nu(\nabla \mathbf{u}, \nabla \mathbf{v})$ 和 $b(q, \mathbf{v}) := (\nabla \cdot \mathbf{v}, q)$，其中 $\mathbf{u}, \mathbf{v} \in X$，$q \in Q$。设三线性形式为 $c(\mathbf{u}, \mathbf{v}, \mathbf{w}) := (2\mathbb{D}(\mathbf{u})\mathbf{v} + (\nabla \cdot \mathbf{u})\mathbf{v}, \mathbf{w})$，其中 $\mathbf{u}, \mathbf{v}, \mathbf{w} \in X$。方程 (1.2) 的有限元半离散格式为：寻找 $\mathbf{u}_h \in X_h \subset X$ 和 $p_h \in Q_h \subset Q$，使得

$$\begin{cases} (\partial_t \mathbf{u}_h, \mathbf{v}) + a_\nu(\mathbf{u}_h, \mathbf{v}) - b(p_h, \mathbf{v}) + c(\mathbf{u}_h, \mathbf{u}_h, \mathbf{v}) = 0, \\ (\nabla \cdot \mathbf{u}_h, q) = 0, \end{cases} \quad \forall \mathbf{v} \in X_h, \forall q \in Q_h.$$

在下文中，我们假设不可压缩约束仅以弱形式施加。Charnyi等人 [1] 证明了基于后向Euler（BE）、BDF2和Crank-Nicolson（CN）方法的完全离散隐式时间步进格式，在应用于弱形式 (1.3) 时，能保持所有三个基本物理性质，即能量耗散、动量守恒和角动量守恒。对于CN格式，这些结构保持性质的详细验证可参见 [1, 子节3.3]。对线性化EMAC型离散化的进一步扩展在 [2] 中进行了研究，文中证明此类格式保留了一部分不变量，通常只能保持能量耗散或线性动量与角动量守恒，但不能同时保持全部三个。对于空间半离散格式 (1.3)，Olshanskii和Rebholz [3, 子节3.1] 证明其结构保持特性具有长时间数值稳定性。此外，利用EMAC对流项的特定抵消结构，他们推导了一个Re半鲁棒的误差估计，其中Gronwall常数独立于雷诺数，尽管未实现完全的压力鲁棒性；参见 [3, 定理3.2, 注3.3]。然而，这些Re半鲁棒结果本质上依赖于对EMAC非线性项的完全隐式处理。因此，该论证并不直接适用于EMAC对流项的显式或半显式处理。

对于EMAC非线性的显式或半显式处理，通常需要额外的稳定化或重构。Han等人 [4] 提出了一种基于标量辅助变量（SAV）重构的一阶后向Euler格式。他们的方法保持了一个修正的离散能量律，同时保持了动量守恒和角动量守恒。然而，所得的误差估计中收敛常数对雷诺数表现出强依赖性；参见 [4, 定理4.1]。在相关方向上，Demir等人 [5] 研究了EMAC型时间滤波方法，并推导了一个二阶精度的格式，具有类似类型的误差界，但仍继承了依赖于雷诺数的常数。为克服这一限制，grad-div稳定化已被广泛用于在不可压缩流模拟中获得Re鲁棒稳定性，特别是在高雷诺数区域；参见例如 [6], [7], [8], [9]。最近，Lan等人 [10] 将grad-div稳定化项 $-\kappa \nabla(\nabla \cdot \mathbf{u})$（其中 $\kappa > 0$）引入EMAC格式，并建立了基于SAV的BE和BDF2格式的Re鲁棒收敛结果，从而消除了 [4] 中观测到的雷诺数依赖性。关于Navier-Stokes方程的动态正则化Lagrange乘子方法的最新进展，我们亦参见 [11], [12], [13]。然而，现有的Re鲁棒分析本质上仍为低阶时间精度，由于高阶BDF方法固有的A稳定性缺失，无法直接推广到高阶格式。因此，本工作的主要目标是为一组高阶隐式-显式BDFk（$k = 3, 4, 5$）格式开发针对Navier-Stokes方程的 $Re$ 鲁棒L2范数误差估计的统一框架。

综上所述，本文的主要贡献如下：

(i) 我们构造了一组IMEX-BDFk格式（$k = 3, 4, 5$），并证明它们在完全离散水平上保持离散线性动量和离散角动量。

(ii) 我们引入了一个全局离散能量框架，用于非A稳定高阶BDF离散化的 $Re$ 鲁棒分析。更精确地说，存在常数 $C$，独立于逆粘度 $\nu^{-1}$，使得对所有 $k = 3, 4, 5$ 有：

$$\|\mathbf{u}(t_n) - \mathbf{u}_h^n\| + \|\nabla \cdot (\mathbf{u}(t_n) - \mathbf{u}_h^n)\| \leq C \exp(t_n)(\tau^k + h^r).$$

(iii) 在上述误差分析的指导下，我们设计了一种自适应变阶策略，根据解的局部时间行为调整BDF阶数，从而在不可压缩流模拟中平衡稳定性与精度。

本文其余部分组织如下。第2节介绍了最高达五阶的IMEX-BDF方法和离散正交卷积（DOC）重构，为贡献(i)提供了简洁的证明。第3节发展了全局离散能量框架并建立了贡献(ii)中的结果。第4节展示了数值实验，证明所提出的IMEX-BDFk格式结合(iii)中的自适应策略，在多尺度流问题上优于固定阶方法。最后，第5节总结了全文。

各节概要

**高阶IMEX-BDF格式**

对于 $s \geq 0$ 和 $1 \leq m \leq \infty$，设 $W^{s,m}(\Omega)$ 表示带有范数 $\|\cdot\|_{W^{s,m}}$ 的标准Sobolev空间；设 $H^s(\Omega) = W^{s,2}(\Omega)$，并用 $|\cdot|_{H^s}$ 和 $\|\cdot\|_{H^s}$ 分别表示相应的半范数和范数。对于 $s = 0$，简记 $L^m(\Omega) = W^{0,m}(\Omega)$，范数写为 $\|\cdot\|_{L^m}$；特别地，在 $L^2(\Omega)$ 上，范数和内积分别记为 $\|\cdot\|$ 和 $(\cdot, \cdot)$。向量值空间定义为 $W^{s,m}(\Omega) = [W^{s,m}(\Omega)]^2$ 和 $H^s(\Omega) = W^{s,2}(\Omega)$，遵循相同约定；对于 $L^2(\Omega)$，为简洁起见，范数和内积同样记为 $\|\cdot\|$ 和 $(\cdot, \cdot)$。设

**鲁棒L2范数误差估计**

本节致力于对所提出的IMEX-BDFk格式 (2.6) 进行 $Re$ 鲁棒的L2范数误差分析。子节3.1引入了一个Re鲁棒的Stokes投影并建立了其正则性性质。子节3.2分析了离散卷积结构引起的一致性误差，并推导了EMAC型对流项的估计。子节3.3为速度误差分析发展了一个统一的Re鲁棒框架。最后，子节3.4建立了相应的Re鲁棒误差估计。

**数值结果**

本节首先评估所提出的IMEX-BDFk格式的精度。然后设计自适应变阶算法，以展示将IMEX-BDFk方法与该策略相结合在多尺度时间相关问题中的优势。最后，将该变阶算法应用于圆柱绕流模拟。

**结论性评述**

本文为最高五阶（$k = 3, 4, 5$）的隐式-显式BDF方法应用于高雷诺数不可压缩Navier-Stokes方程发展了一个统一的Re鲁棒收敛框架。主要的分析困难源于高阶BDF方法缺乏A稳定性，这阻碍了标准能量论证的直接使用。为克服这一困难，我们采用DOC核技术，将高阶BDF方法重构为等价的低三角…

**致谢**

作者感谢编辑和匿名审稿人提出的宝贵意见和建议，这些意见改进了本文的质量。

**CRediT作者贡献声明**

Ji Bingquan：撰写——原稿。Wang Yushun：概念构思。Gong Yuezheng：项目管理。

**利益冲突声明**

作者声明不存在已知的可能影响本文所报告工作的竞争性经济利益或个人关系。

Ji Bingquan | Wang Yushun | Gong Yuezheng

教育部NSLSCS重点实验室，南京师范大学数学科学学院，南京，210023，中国