0 引言

1959年，华罗庚先生在《大哉数学之为用》(数与量)^[1]一文中关于“时空”有以下论述：四维空间听来好象有些神秘，其实早已有之，即以“宇宙”二字来说，“往古来今谓之宙，四方上下谓之宇”(《淮南子·齐俗训》)，就是宇是东西、南北、上下三维扩展的空间，而宙是一维的时间。

接着他引用相对论和生活的常识叙述了时空既独立又统一的性质：爱因斯坦不再把“宇”、“宙”分开来看，也就是时间也在进行着。每一瞬间三维空间中的物质都占有它一定的位置。他根据麦克斯韦——洛伦兹的光速不变假定，并继承牛顿的相对性原理，提出了狭义相对论。狭义相对论中的洛伦兹变换把时空联系在一起，但并不消灭时空特点。如向东走三里再向西走三里就回到原处，但时间则不然，共用了走六里的时间。时间是一去不复返地流逝着。

华先生意在说明数学在描述“宇宙之大”中的作用。本文的引用，旨在强调时空之间的关系深刻地植根于计算流体力学的模型理解和算法构造。离开了时空演化，一切归于“稳态”。正因为时空演化，流体的世界才色彩斑澜、波澜壮阔。

定性甚至定量描述流体世界的时空演化并不容易，数值模拟同样是貌似简单，实则不易。下面举一个众所周知的例子:

$ {\partial _t}u(x,t) + a{\partial _x}u(x,t) = 0 $

(1)

其中a是一个常数，x表示空间坐标，t是时间变量。此模型表示了变量u(x, t)的空间变化(扰动、波)与时间变化(传播)的关系，描述了一个时空关联的性质。

$ u(x,t + \Delta t) = u(x - a\Delta t,t) $

(2)

式(2)事实上比式(1)更加根本和直接，并且不牵涉函数u(x, t)本身更多性质, 比如正则性。即使u(x, t)是一个方波或脉冲, 式(2)仍然有效。用平移算子e^-aΔt∂_x来表示式(2)，则有：

$ u(x,t + \Delta t) = u(x - a\Delta t,t) = {{\rm{e}}^{ - a\Delta t{\partial _x}}}u(x,t) $

(3)

从这里出发，如果想精确表达式(1)的输运性质，则需要进行展开(如泰勒展开):

$ \begin{array}{c} u(x, t+\Delta t)=u(x-a \Delta t, t)=\mathrm{e}^{-a \Delta t \partial_x} u(x, t) \\ =u(x, t)-a \Delta t \partial_{x} u(x, t)+\frac{\left(-a \Delta t \partial_{x}\right)^{2}}{2} u(x, t)+\cdots \end{array} $

(4)

为了设计实际工程中的计算方法(简称算法)，需进行必要操作：一是截断处理，涉及算法与控制方程的相容性，以及精度的概念；二是对空间变化∂_xu(x, t)的处理(如有限差分)，不同的处理得到不同的算法。这些都与函数u(x, t)的正则性和模型(1)的内在性质息息相关：正则性越好，精度越高; 模型(1)的时空关联性需要通过式(2)精确描述。这些处理充满技巧和无限可能，怎样“令人信服”地设计时空耦合算法且精准描述相应模型的内在时空关联性，是一个需要深入思考的本质问题。

近年来，计算流体力学蓬勃发展，各种算法和相应的应用软件目不暇接。有些算法具有坚实的理论基础，如基于变分原理的有限元方法等。而流体力学中关于可压缩流动的研究，由于流场结构非常复杂(激波、物质界面、涡团和湍流等), 人们对流动(解的)性质缺乏很好理解，模型的适定性和相关数值模拟中算法的探讨相当工程化，形式化的表达往往导致似是而非的结果，常常是模型具有很好的性质，而相应数值模拟的离散模型失去了该保持的性质，比如模型的时空关联性和算法的时空耦合性等。不同的出发点导出的数值算法往往截然不同。本文试图在这方面做一点尝试: 基于连续介质假定探讨流体力学模型时空变化的内在关联性，阐述相应算法应保持的时空耦合性质。

从流体力学微团法的直接建模开始思考模型的时空关联性。通过对通量的研究，发现经典的瞬时Cauchy通量无法反映时空关联性，提出了某一时间段内时间区间通量，在任意特定时间段内反映流体的动力学过程^[2-3]，进而提出有限体积格式的基本原理：在连续介质假定性下，任一时间段内流体通过控制体边界的时间区间通量关于边界的扰动是Lipschitz连续的。这一Lipschitz连续性恰好体现流场时空关联性的本质特征，实现了数学上流体力学方程组弱解和物理上积分型平衡律之间的统一，后者其实就是全离散有限体积格式的物理起源。通过对流体力学方程组时空关联性质的研究，实现有限体积框架下的时空耦合。在算法的实现过程中，物理量(质量、动量、能量等)以及它的变化率这一对量起着重要的作用。一个量的变化率通过它的空间变化关系来反映。变化率不是一类常规的数学演算，而是反映了内在的物理规律。例如，加速度是速度的变化率，它等于控制体表面受到的力，即应力，这是牛顿第二定律。在流体力学计算方法中使用这样的量是自然的也是必须的，这反映了时空变化的联系。在文献[4]中，笔者倡导的时空耦合高精度算法正是这一思想的具体表现。本文再次描述实现时空耦合算法的一种基本流程，并通过算例展示时空耦合算法的数值表现。

从偏微分方程理论上的Cauchy-Kowalevski方法到数值解中的Lax-Wendroff方法，都是某种意义下的时空耦合方法。现代计算流体力学的算法时空观常有体现，如迎风格式、广义黎曼解法器^[5-7]、守恒元/解元(CE/SE)方法^[8]、气体动理学解法器^[9-10]等。最近几年，笔者致力于思考高精度时空耦合算法的必要性与理论基础，但远远不够深入。随着计算技术的提高和工程需求的精细化，这方面的研究应该得到重视。应当建立相应的时空观，这对算法的设计与实现很有帮助。故作此文，抛砖引玉。

1 流体力学方程组的时空关联性

众所周知，流体力学建模过程常用微团法，即在某一时间段(t, t+δt)，研究控制体Ω(t)=(x, x+δx)内流体质团的运动。时间间隔δt非常重要，给出了微团运动的动力学过程响应时间。为了描述方便，基于连续介质假定，流体力学方程组写成如下形式：

给定一个控制体Ω(t), 流体运动满足下列的积分平衡律(略去外力场等效应)：

$ \frac{\mathrm{d}}{\mathrm{d} t} \int_{\varOmega(t)} \boldsymbol{u}(\boldsymbol{x}, t) \mathrm{d} \boldsymbol{x}=-\int_{\partial \varOmega(t)} \boldsymbol{f} \cdot \boldsymbol{n} \mathrm{d} \varGamma $

(5)

这里u是密度函数，n是∂Ω(t)的单位外法向。相应地，称:

$ \boldsymbol{m}(t)=\int_{\varOmega(t)} \boldsymbol{u}(\boldsymbol{x}, t) \mathrm{d} \boldsymbol{x} $

(6)

为控制体Ω(t)上的质量，右边:

$ \boldsymbol{\mathcal{C}} (\partial \varOmega ; t)=-\int_{\partial \varOmega(t)} \boldsymbol{f} \cdot \boldsymbol{n} \mathrm{d} \varGamma $

(7)

为Cauchy通量，f为通量密度函数，简称为通量函数，它是u以及空间变化量∇u等的函数f=f(u, ∇u, …)。方程(5)反映了质量时间变化率与通过控制体边界∂Ω(t)的通量之间的瞬时关系。

Cauchy在特定连续介质假定下，论述了Cauchy通量的性质^[11-12]，这是他对连续介质力学最重要的贡献^[11]。要深入理解这一关系并不容易，这与Gauss-Green定理密切相关。基于这一定理，可以将式(5)写成散度型偏微分方程(PDE)形式：

$ {\partial _t}\mathit{\boldsymbol{u}} + \nabla \cdot \mathit{\boldsymbol{F}}(\mathit{\boldsymbol{u}},\nabla \mathit{\boldsymbol{u}}, \cdots ) = 0. $

(8)

这里，略去雷诺输运定理等细致讨论，并设Ω不随流场运动，即欧拉型控制体。相应地，式(1)中的控制体为拉格朗日型控制体。

方程组(8)对光滑流动是有效的，可由偏微分方程的知识建立时空的关联性。注意这里的通量函数F(u, ∇u, …)与式(7)略有不同，它涉及雷诺输运定理等有关变换，这里不做详细讨论。以后欧拉控制体边界上的通量用大写字母表示，否则用小写字母。对于实际流动来说，这一转化是非常粗糙的，微妙之处在于Cauchy通量$\boldsymbol{\mathcal{C}}$(∂Ω; t)的数学性质，即正则性。可压缩流场蕴涵丰富的现象，如冲击波、物质界面、湍流等效应，流场(或称方程的解)的正则性非常弱，人们对它的认识很少，故而该定理的应用并不显然。历史上有很多研究^[13-16]，直到最近，这方面的研究仍是如火如荼^[17]，遗憾的是所取得的结果与实际的期望仍然相差甚远。正如最近的专著^[18]的1.3节有这样一段论述, “Regarding the right-hand side of (1.3) one needs to keep in mind the following comment concerning the identification of the boundary flux: the drawback of this functional analytic, demonstration is that it does not provide any clues on how the q_D may be computed from A”，其中q_D对应于这里的f(u)·n，A对应于一个给定的应力张量。也就是说：人们还不知道如何从应力张量中计算出通量密度函数。

仔细审视可以看出Cauchy通量$\boldsymbol{\mathcal{C}}$(∂Ω; t)是表示一个特定时刻t流动的瞬时行为，方程(5)只是表示一个非常弱的“时-空”关联性质。现在常用的一个观点是下面的“弱解”概念，它来自于偏微分方程(8)，利用了分布理论。

定义1.1  设Ω是$\mathbb{R}$ⁿ中的一个有界区域，n=1, 2, 3, 0≤t₁≤t₂。对任意试验函数φ(x, t)∈C₀^∞(Ω×[t₁, t₂]), 如果函数u(x, t)满足:

$ \int_{{t_1}}^{{t_2}} {\int_\mathit{\Omega } {\left[ {\mathit{\boldsymbol{u}}(\mathit{\boldsymbol{x}},t){\partial _t}\phi + \mathit{\boldsymbol{F}} \cdot \nabla \phi } \right]} } {\rm{d}}x{\rm{d}}t = 0 $

(9)

则称u(x, t)是方程(8)的弱解。

这个表达式比式(5)前进了一大步：如果解u(x, t)是光滑的，式(8)和式(9)是等价的。近年来偏微分方程理论得到极大的发展，人们对于式(8)的认识比直接对式(5)的认识要深刻得多，且偏微分方程(8)的时空关联性一目了然：流场的空间摄动可以通过式(8)反映到时间的变化中, 线性波的传播式(1)正是这种时空关联性的特例。另外，当一个间断面在无限薄的假定下，可以导出Rankine-Hugoniot间断关系，即间断面(线)两侧解的“迹”的关系。

定义1.1并没有对函数u(x, t)做过多假设，只要式(9)成立即可。一旦流场出现复杂间断或其它流场结构时，这个定义并不能给出太多的信息。我们回避不了的事实是：在可压缩流动中，流场即使初始性质非常“好”，但在不久的将来由于复杂的非线性相互作用，可能变得非常“糟糕”。这可以从最简单的Burgers方程看到^[19-20]。因此，人们不得不面对复杂的情形，深入思考为什么应用Gauss-Green定理理解Cauchy通量非常困难。

任何流动都有一个动力学过程，这也许是个常识, 但仍需要严格论证。最近的研究表明^[3-4]，流动的时空关系既相互独立又彼此关联。下面的原理深刻地反映这一事实。

有限体积方法基本原理。设u(x, t)是定义1.1下方程(8)的弱解，Ω是$\mathbb{R}$ⁿ中任一紧集。假定:

(ⅰ) u(x, t)是局部有界；

(ⅱ) 质量m(t)是时间t的连续函数。

那么有下列结论：

(ⅰ) 对任意δ>0, 函数H(x; t, t+δt)=$ \int_t^{t + \delta t} $F(x, t)dt满足:

$ \nabla_{\boldsymbol{x}} \cdot \boldsymbol{H}(\boldsymbol{x} ; t, t+\delta t) \in L_{\mathrm{loc}}^{\infty}\left(\mathbb{R}^{n}\right) $

(10)

(ⅱ) 对任意δ>0，时间段(t, t+δt)内流过Ω边界∂Ω的时间区间通量:

$ \mathcal{F}(\partial \varOmega ; t, t+\delta t)=\int_{t}^{t+\delta t} \int_{\partial \varOmega} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma \mathrm{d} t $

(11)

关于Ω的边界∂Ω扰动是Lipschitz连续的。

时间区间通量$\mathcal{F}$(∂Ω; t, t+δt)和Cauchy通量$\boldsymbol{\mathcal{C}}$(∂Ω；t)明显不同，特别是Lipschitz连续性，进一步反映了流场的时空关联性质；而Cauchy通量只是一个瞬时行为，不能完全反映动力学的过程。直观地来说，这里通过时间效应，换取空间的正则性，这是和Cauchy通量数学性质上的根本差异。

由此，我们给出方程(8)的弱解另一种表述，它恰恰是物理上常见的积分型平衡律。

定义1.2  设Ω是$\mathbb{R}$ⁿ中任一紧集, δt>0。如果它满足下列条件：

(ⅰ) u(x, t)局部有界，且质量m(t)是时间t的连续函数；

(ⅱ) 整体通量(4)有意义且关于Ω的边界∂Ω扰动是连续的；

(ⅲ) 积分平衡律成立:

$ \int_{\varOmega} \boldsymbol{u}(\boldsymbol{x}, t+\delta t) \mathrm{d} \boldsymbol{x}-\int_{\varOmega} \boldsymbol{u}(\boldsymbol{x}, t) \mathrm{d} \boldsymbol{x}=-\int_{t}^{t+\delta t} \int_{\partial \Omega} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma \mathrm{d} t $

(12)

其中n是∂Ω的单位外法向, 则称函数u(x, t)是偏微分方程(8)的弱解。

事实上，已经证明定义1.1和定义1.2是等价的^[3-4]，并且给出了通量(11)的正则性质。这个定义说明后面讨论的有限体积方法就是计算流体力学的基本格式，和物理最原始的建模一致。流体力学方程组(8)可理解为：如果流场光滑，用偏微分方程组(8)刻画；但当流场失去正则性时，解需要满足积分平衡律式(12)。流体力学有限体积格式的构造过程就是基于式(12)的直接建模过程：在偏微分方程(8)诱导出的时空关联解辅助下，构造出和式(12)相容的离散模型，精准反映真实的物理过程。

2 有限体积方法及其时空耦合性

前面论述到：有限体积方法构造过程就是在偏微分方程(8)辅助下的一个直接建模过程。通常的有限体积方法从偏微分方程组(8)出发，在每一个计算单元(控制体)上应用Gauss-Green定理，得到有限体积公式。具体地，将计算区域Ω剖分为N个计算单元(控制体)${\mathit{\Omega} _j}, \mathit{\Omega} = \cup _{j = 1}^N{\mathit{\Omega} _j}, \partial {\mathit{\Omega} _j} = { \cup _\ell }{\mathit{\Gamma} _{jl}}$。有下列两种形式：半离散和全离散有限体积方法。通过上一节的讨论看出，后者具有时空耦合性。

2.1 半离散有限体积方法

在Ω_j上，对式(8)的空间散度项应用Gauss-Green公式, 可得半离散有限体积方程:

$ \frac{\mathrm{d}}{\mathrm{d} t} \int_{\varOmega_{j}} \boldsymbol{u}(\boldsymbol{x}, t) \mathrm{d} \boldsymbol{x}=-\int_{\partial \varOmega_{j}} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma, \quad j=1, \cdots, N $

(13)

它直接对应于流体力学方程组(5)。在初始时刻t=0, 给定u(x, t)的分布:

$ \mathit{\boldsymbol{u}}(\mathit{\boldsymbol{x}},0) = \mathit{\boldsymbol{P}}_0^k(x) \in {V^k} $

(14)

其中V^k是一个可持续函数空间(通常为分片多项式函数空间，用高阶多项式等光滑函数去逼近正则性很弱的函数, 如间断函数，常常会导致振荡现象，数据重构/数据投影仍是有限体积格式中未彻底解决的难题，意味着在每个时间层t=t_n>0, u(x, t_n)经过投影后所得P_n^k∈V^k中)。由于守恒性的基本要求:

$ \int_{\varOmega_{j}} \boldsymbol{u}\left(\boldsymbol{x}, t_{n}\right) \mathrm{d} \boldsymbol{x}=\int_{\varOmega_{j}} \boldsymbol{P}_{n}^{k}(\boldsymbol{x}) \mathrm{d} \boldsymbol{x} $

(15)

P_n^k(x)一般是沿着每个控制体Ω_j的边界Γ_jl是间断的。为了利用式(13)更新解，需要沿着每个控制体Ω_j的边界Γ_jl求解式(8)，得到瞬时解。数值上，使用后面所描述的黎曼解法器或近似黎曼解法器，构造时刻t的数值通量:

$ \left| {{\mathit{\Gamma }_{jl}}} \right|{\mathit{\boldsymbol{F}}_{jl}}(t) \cdot {\mathit{\boldsymbol{n}}_{jl}} \approx \int_{\mathit{\Gamma }_{jl}} \mathit{\boldsymbol{F}} \cdot \mathit{\boldsymbol{n}}{\rm{d}}\mathit{\Gamma } $

(16)

将之带入(13)即得:

$ \begin{array}{*{20}{c}} {\quad \frac{{\rm{d}}}{{{\rm{d}}t}}\int_{{\mathit{\Omega }_j}} \mathit{\boldsymbol{u}} (\mathit{\boldsymbol{x}},t){\rm{d}}\mathit{\boldsymbol{x}} = - \int_{\partial {\mathit{\Omega }_j}} \mathit{\boldsymbol{F}} \cdot \mathit{\boldsymbol{n}}{\rm{d}}\mathit{\Gamma }}\\ { \approx - \sum\limits_\ell {\left| {{\mathit{\Gamma }_{jl}}} \right|} {\mathit{\boldsymbol{F}}_{jl}}(t) \cdot {\mathit{\boldsymbol{n}}_{jl}},\quad j = 1, \cdots ,N} \end{array} $

(17)

接着使用常微分方程求解器，如典型的Runge-Kutta方法，求解这个时间依赖(常微分)方程。需要指出的是：这个方程需要在每个时间层t=t_n求解数值通量，并把式(13)在下一时间层的解投影到V^k中，得到P_n+1^k∈V^k。如果用每个时间步用多级Runge-Kutta方法，则需要多次求解数值通量和多次投影。这个投影过程又称为数据重构。用算子形式，半离散方法可以表示为：

$ \boldsymbol{P}_{n+1}^{k}(x)=\mathcal{P} \cdot \mathcal{R} \mathcal{K} \cdot \mathcal{A}_{r i e}\left[\boldsymbol{P}_{n}^{k}(\boldsymbol{x})\right], $

(18)

其中$\mathcal{A}$_rie表示用黎曼解法器构造瞬时数值通量，$\mathcal{RK}$代表Runge-Kutta型的时间推进过程，$\mathcal{P}$表示数据重构。

我们注意到，由于对Cauchy通量∫_ΓjlF·ndΓ还没有清晰的认识，不得不借助局部黎曼解法器(见下节)进行数值通量(16)的逼近。一般来说，所得到的局部误差为:

$ \left|\varGamma_{j l}\right| \boldsymbol{F}_{j l}(t) \cdot \boldsymbol{n}_{j l}-\int_{\varGamma_{j l}} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma=\mathrm{O}\left(\left|(\Delta \boldsymbol{u})_{j l}\right|\right)\left|\varGamma_{j l}\right| $

(19)

这里|(Δu)_jl|表示数据P_n^k(x)跨过控制体边界Γ_jl的局部跳跃。

半离散有限体积方法属于线方法，从某种意义上来说是时空解耦方法。在光滑流场中，流体力学方程组的偏微分方程关系可以直接反映时空关联性。但对间断解来说，式(19)中的局部误差的累积会给数值模拟带来巨大伤害。

2.2 全离散有限体积方法

本文将把式(12)应用到时空控制体Ω_j×[t_n, t_n+1], t_n+1=t_n+Δt_n, Δt_n>0为时间步长，得到有限体积格式的全离散形式:

$ \int_{\varOmega_{j}} \boldsymbol{u}\left(\boldsymbol{x}, t_{n+1}\right) \mathrm{d} \boldsymbol{x}-\int_{\varOmega_{j}} \boldsymbol{u}\left(\boldsymbol{x}, t_{n}\right) \mathrm{d} \boldsymbol{x}=-\int_{t_{n}}^{t_{n+1}} \int_{\partial \varOmega_{j}} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma \mathrm{d} t $

(20)

也可以把式(20)看作是式(13)在时间段(t_n, t_n+1)上的积分。与半离散情形下求解瞬时通量(16)不同，这里需要利用偏微分方程(8)的时空关联性质逼近[t_n, t_n+1]上时间区间通量:

$ \sum\limits_{\ell}\left|\varGamma_{j l}\right| \boldsymbol{F}_{j l}\left(t_{n}, t_{n+1}\right) \cdot \boldsymbol{n}_{j l} \approx \sum\limits_{\ell} \int_{t_{n}}^{t_{n+1}} \int_{\varGamma_{j l}} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma \mathrm{d} t $

(21)

将之代入(20)，得到有限体积公式:

$ \overline{\boldsymbol{u}}_{j}^{n+1}=\overline{\boldsymbol{u}}_{j}^{n}-\frac{1}{\left|\varOmega_{j}\right|} \sum\limits_{\ell}\left|\varGamma_{j l}\right| \boldsymbol{F}_{j l}\left(t_{n}, t_{n+1}\right) \cdot \boldsymbol{n}_{j l} $

(22)

这里记逼近解在控制体Ω_j上的平均值为:

$ \bar{u}_{j}^{n}=\frac{1}{\left|\varOmega_{j}\right|} \int_{\varOmega_{j}} \boldsymbol{u}\left(\boldsymbol{x}, t_{n}\right) \mathrm{d} \boldsymbol{x} $

(23)

由有限体积方法基本原理，可以得到:

$ \begin{array}{c} \left|\varGamma_{j l}\right| \boldsymbol{F}_{j l}\left(t_{n}, t_{n+1}\right) \boldsymbol{n}_{j l}-\int_{t_{n}}^{t_{n+1}} \int_{\varGamma_{j l}} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma \mathrm{d} t \\ =\mathrm{O}\left(\Delta t_{n}\right)^{q+1}\left|\varGamma_{j l}\right| \end{array} $

(24)

其中q>0。值得注意的是，式(24)中的误差和式(19)中的误差非常不同, 式(19)中的误差是用解的局部跳跃(变差)来测量，而式(24)中的误差直接用时间步长(等价于网格尺度)来测量。当流场相对光滑时，这两者是等价的, 但当流场中有剧烈间断、脉动时，两者差别是明显的。即使网格加密，式(19)也得不到收敛解(在标量方程的研究中，网格加密是可以得到收敛解的，原因在于全局变差有界性质。到目前为止还没有例证可证明该性质在流体力学方程组成立^[7])。后文算例4.1可以看出数值通量的重要性。

同半离散方法相似，在每个时间层t=t_n上投影可得到在可持续空间V^k中的逼近解P_n^k。全离散有限体积方法可以用符号表示为：

$ \boldsymbol{P}_{n+1}^{k}(x)=\mathcal{P} \cdot \mathcal{E} \cdot \mathcal{A}\left[\boldsymbol{P}_{n}^{k}(\boldsymbol{x})\right] $

(25)

这里$\mathcal{A}$表示通量的逼近，$\mathcal{E}$表示解的有限体积发展过程，$\mathcal{P}$代表投影过程(数据重构)。一旦$\mathcal{P}$被确定，$\mathcal{E}$是自然的且没有任何误差产生。数据重构部分是重要的，将在后面评注。

2.3 有限体积方法的相容性

所谓相容性，描述的是离散格式与背景方程之间的关系。传统上常常使用Taylor展开的方式研究数值格式与偏微分方程(如式(8))之间的相容性，这样的运算只对光滑流场成立。当流场比较复杂时，目前大部分的数值分析某种意义上只是启发性的，比如以标量方程为模型建立相关的理论^[21]。

现在基于积分平衡律(12)来建立有限体积格式的相容性，即离散格式(18)或(25)须与式(12)进行直接比较。由式(18)或式(25)可知, 误差来源于数据投影步$\mathcal{P}$以及通量的逼近步$\mathcal{A}$_rie或$\mathcal{A}$。对于前者，虽然有大量的文献存在，数据重构步主要依赖相关的空间信息。文献[22]只是做了初步的尝试，把未来时间的数据重构与过去的信息进行了关联。本文假定投影步满足所有的精度要求，讨论通量的逼近误差。

对半离散格式(18)来说，在每个固定时刻给出的通量误差至多为:

$ \begin{aligned} \sum\limits_{\ell}\left|\varGamma_{j l}\right| & \boldsymbol{F}_{j l}(t) \cdot \boldsymbol{n}_{j l}-\sum\limits_{\ell} \int_{\varGamma_{j l}} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma \\ &=\mathrm{O}\left((\Delta \boldsymbol{u})_{j l}\right)\left|\varGamma_{j l}\right| \\ &=\mathrm{O}\left((\Delta \boldsymbol{u})_{j}\right)\left|\varGamma_{j}\right| d_{j} \end{aligned} $

(26)

其中(Δu)_j表示在控制体附近解的局部跳跃，|Γ_j|是Ω_j边的最大长度，d_j=diam(Ω_j)。并且在每个时间层，随着Runge-Kutta步的增加，误差也会累积。在式(26)中，d_j来源于不同方向上通量差的获利，在3.5节可以进一步看到。

对于全离散方法(25)，我们讨论每个时间层数值通量的逼近，并期望有下列的估计:

$ \begin{array}{c} \sum\limits_{\ell}\left|\varGamma_{j l}\right| \boldsymbol{F}_{j l}\left(t_{n}, t_{n+1}\right) \cdot \boldsymbol{n}_{j l}-\sum\limits_{\ell} \int_{t_{n}}^{t_{n+1}} \int_{\varGamma_{j l}} \boldsymbol{F} \cdot \boldsymbol{n} \mathrm{d} \varGamma \mathrm{d} t \\ =\mathrm{O}\left(\Delta t_{n}\right)^{q+2}\left| \Gamma_{j}\right|. \end{array} $

(27)

其中q>0。比较式(26)和式(27)，它们有本质的差别，即(Δu)_j与Δt_n的差别。对光滑解来说它们是等价的。这样我们给出下列的定义。

定义2.1 (有限体积格式的相容性)。设Ω_j是任一控制体，j=1, …, N, 如果对于某个q>0, 式(27)成立，则有限体积格式(25)相容于平衡律(12)，并具有q阶相容性。

如上所述，相容性概念常常相对与偏微分方程所言。如果式(8)是双曲守恒律，Lax和Wendroff给出了有限差分方法的相容性，常称为Lax相容性^[23]。由于Lax相容性概念影响深远，有必要进行回顾。

定义2.2 (Lax相容性^[23])。考虑双曲守恒律方程组:

$ \partial_{t} \boldsymbol{u}+\partial_{x} \boldsymbol{f}(\boldsymbol{u})=0 $

(28)

其2p+1点守恒性差分格式为:

$ \boldsymbol{u}_{j}^{n+1}=\boldsymbol{u}_{j}^{n}-\lambda\left[\boldsymbol{g}_{j+\frac{1}{2}}^{n}-\boldsymbol{g}_{j-\frac{1}{2}}^{n}\right], \quad \lambda=\frac{\Delta t}{\Delta x} $

(29)

其中Δx是网格尺寸, Δt是时间步长, $\mathit{\boldsymbol{g}}_{j + \frac{1}{2}}^n = \mathit{\boldsymbol{g}}(\mathit{\boldsymbol{u}}_{j - p + 1}^n, \ldots, \mathit{\boldsymbol{u}}_{j + p}^n)$。如果成立:

$ \mathit{\boldsymbol{g}}(\mathit{\boldsymbol{u}}, \cdots ,\mathit{\boldsymbol{u}}) = \mathit{\boldsymbol{f}}(\mathit{\boldsymbol{u}}) $

(30)

则称差分方法(29)与双曲守恒律(28)相容。

基于这个定义，建立了影响深远的Lax-Wendroff收敛定理，直至今日，该定理仍被广泛应用。随着高阶精度格式的发展以及工程应用的需求，这个定理常常被误用，典型表现为：

(ⅰ) 相容性定义(30)只适用于传统一阶精度数值方法，对于目前广泛使用的高阶精度方法并不适用，如高阶有限体积方法等。在定义2.2中，数值通量$\mathit{\boldsymbol{g}}_{j + \frac{1}{2}}^n$并不包含P_n^k(x)的空间变化信息。

(ⅱ) 在此基础上建立的Lax-Wendroff定理只适用于一致网格或结构网格，对非结构网格并不适用^[24-25]。数值验证过程中使用的网格加密方法测试收敛阶，需要倍加小心。

实际上，定义2.1第一次给出高精度有限体积数值格式与积分平衡律之间的相容性^[3]。有限体积方法与网格的结构无关，因此适用于无结构网格。特别地，收敛阶测试时一般针对光滑流进行，数据重构部分能够达到应有的精度，通量逼近阶事实上就是收敛阶。但是，当流场含有间断或其他复杂结构时，数据重构仍存在诸多争议，通量的逼近效果往往被忽略。通过以上分析可以看到，如果通量相容性(27)不成立，逼近解不可能收敛。也就是说，式(27)是有限体积格式相容性的一个必要条件。

2.4 再访Godunov方法

谈到流体力学中的有限体积方法，需要回顾Godunov方法^[26]。经过半个多世纪的发展和检验，它已成为现代计算流体力学基石。考虑一维可压缩欧拉方程组:

$ {\partial _t}\mathit{\boldsymbol{u}} + {\partial _\mathit{\boldsymbol{x}}}\mathit{\boldsymbol{F}}(\mathit{\boldsymbol{u}}) = 0 $

(31)

其中守恒量u=(ρ, ρu, ρE)^T, F(u)=(ρu, ρu²+p, u(ρE+p))^T。这个方程组由Gibbs关系式封闭，给出热力学量之间的关系，这里不再赘述。记${I_j} = \left({{x_{i - \frac{1}{2}, }}, {x_{j + \frac{1}{2}}}} \right), {x_j} = \frac{1}{2}\left({{x_{j - \frac{1}{2}}} + {x_{j + \frac{1}{2}}}} \right), $$ \Delta {x_j} = {x_{j + \frac{1}{2}}} - {x_{j - \frac{1}{2}}}, {t_{n + 1}} = {t_n} + \Delta {t_n} $，时空控制体C_jⁿ=I_j×[t_n, t_n+1]。那么式(31)可以写成:

$ \begin{array}{c} \int_{I_{j}} \boldsymbol{u}\left(x, t_{n+1}\right) \mathrm{d} x-\int_{I_{j}} \boldsymbol{u}\left(x, t_{n}\right) \mathrm{d} x \\ +\int_{t_{n}}^{t_{n+1}} \boldsymbol{F}\left(\boldsymbol{u}\left(x_{j+\frac{1}{2}}, t\right)\right) \mathrm{d} t \\ -\int_{t_{n}}^{t_{n+1}} \boldsymbol{F}\left(\boldsymbol{u}\left(x_{j-\frac{1}{2}}, t\right)\right) \mathrm{d} t=0 \end{array} $

(32)

即积分平衡律(12)的形式。对应于式(22)有限体积格式为:

$ \overline{\boldsymbol{u}}_{j}^{n+1}=\overline{\boldsymbol{u}}_{j}^{n}-\frac{\Delta t_{n}}{\Delta x_{j}}\left[\boldsymbol{F}_{j+\frac{1}{2}}^{n}-\boldsymbol{F}_{j-\frac{1}{2}}^{n}\right] $

(33)

其中F_{${j + \frac{1}{2}}$}ⁿ是数值通量。给定初始的数据分布u(x, t_n)=P_n⁰∈V⁰, Godunov格式的一个中心思想是在每个网格边界点(x_{${j + \frac{1}{2}}$}, t_n)处求解相应黎曼问题:

$ \begin{array}{c} {\partial _t}\mathit{\boldsymbol{u}} + {\partial _x}\mathit{\boldsymbol{F}}(u) = 0,x \in \mathbb{R},t > {t_n},\\ \mathit{\boldsymbol{u}}\left( {x,{t_n}} \right) = \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{\bar u}}_j^n,}&{x < {x_{j + \frac{1}{2}}}}\\ {\bar u_{j + 1}^n,}&{x > {x_{j + \frac{1}{2}}}} \end{array}} \right. \end{array} $

(34)

得到网格边界上解的值u_{${j + \frac{1}{2}}$}ⁿ。对于简单的状态方程，如理想气体状态方程，欧拉方程组(31)的黎曼解可以清晰给出，记为u(x, t)=R(ξ, u_jⁿ, u_j+1ⁿ), $\xi = \frac{{x - {x_{j + \frac{1}{2}}}}}{{t - {t_n}}}, \mathit{\boldsymbol{u}}_{j + \frac{1}{2}}^n = \mathit{\boldsymbol{R}}(0;\mathit{\boldsymbol{\bar u}}_j^n, \mathit{\boldsymbol{\bar u}}_{j + 1}^n)$。这样F_{${j + \frac{1}{2}}$}ⁿ=F(u_{${j + \frac{1}{2}}$}ⁿ), 从而得到Godunov格式:

$ \overline{\boldsymbol{u}}_{j}^{n+1}=\overline{\boldsymbol{u}}_{j}^{n}-\frac{\Delta t_{n}}{\Delta x_{j}}\left[\boldsymbol{F}\left(\boldsymbol{u}_{j+\frac{1}{2}}^{n}\right)-\boldsymbol{F}\left(\boldsymbol{u}_{j-\frac{1}{2}}^{n}\right)\right] $

(35)

由于

$ \frac{1}{\Delta t_{n}} \int_{t_{n}}^{t_{n+1}} \boldsymbol{F}\left(\boldsymbol{u}\left(x_{j+\frac{1}{2}}, t\right)\right) \mathrm{d} t \equiv \boldsymbol{F}\left(\boldsymbol{u}_{j+\frac{1}{2}}^{n}\right) $

可知式(35)和式(12)完全一致。因此，从积分平衡律的逼近来说，基于分片常数的初值逼近，式(35)与式(12)完全相容或通量有无穷阶相容性。从整体逼近式(25)来说，所有的误差来自于投影过程。因此习惯上称Godunov格式为一阶格式。

如果用逼近的黎曼解法解，界面上的值一般可以写为:

$ \widetilde{\boldsymbol{u}}_{j+\frac{1}{2}}^{n}=\frac{1}{2}\left(\overline{\boldsymbol{u}}_{j}^{n}+\overline{\boldsymbol{u}}_{j+1}^{n}\right)-\frac{\alpha\left(\overline{\boldsymbol{u}}_{j}^{n}, \overline{\boldsymbol{u}}_{j+1}^{n}\right)}{2}\left[\boldsymbol{F}\left(\overline{\boldsymbol{u}}_{j+1}^{n}\right)-\boldsymbol{F}\left(\overline{\boldsymbol{u}}_{j}^{n}\right)\right] $

(36)

这里α(u_jⁿ, uⁿ_j+1)称为黏性系数。一般来说，通量的逼近误差为:

$ \begin{array}{c} \frac{1}{{\Delta {t_n}}}\int_{{t_n}}^{{t_{n + 1}}} \mathit{\boldsymbol{F}} \left( {\mathit{\boldsymbol{u}}\left( {{x_{j + \frac{1}{2}}},t} \right)} \right){\rm{d}}t - \mathit{\boldsymbol{F}}\left( {{\mathit{\boldsymbol{\tilde u}}^n}_{j + \frac{1}{2}}} \right)\\ = {\rm{O}}\left( {\left| {\mathit{\boldsymbol{\bar u}}_{j + 1}^n - \mathit{\boldsymbol{\bar u}}_j^n} \right|} \right) \end{array} $

(37)

这个误差在强间断的情形下对计算结果有巨大的伤害。如式(19)所述，这个误差并不随着网格加密而减小。

在多维情形下，例如二维双曲守恒律情形:

$ {\partial _t}\mathit{\boldsymbol{u}} + {\partial _x}\mathit{\boldsymbol{F}}(\mathit{\boldsymbol{u}}) + {\partial _y}\mathit{\boldsymbol{G}}(\mathit{\boldsymbol{u}}) = 0 $

(38)

在每个时间层t=t_n的数据为分片常数，即每个控制体Ω_j内为常数u_jⁿ, 通过求解相应的黎曼问题, 构造数值通量。这时相应的黎曼问题可分为两类：

(ⅰ) 相对于界面的法向黎曼问题

由于流体力学方程组的伽利略不变性，通过旋转变换，总可以设x-方向为Γ_jl的外法向，Γ_jl两侧单元的值记为u_L和u_R，法向黎曼问题定义为:

$ \begin{array}{c} \partial_{t} \boldsymbol{u}+\partial_{x} \boldsymbol{F}(\boldsymbol{u})=0, \quad(x, y) \in \mathbb{R}^{2}, t>0 \\ \boldsymbol{u}(x, y, 0)=\left\{\begin{array}{ll} \boldsymbol{u}_{L}, & x<0 \\ \boldsymbol{u}_{R}, & x>0 \end{array}\right. \end{array} $

(39)

它的求解和上面的一维黎曼问题(34)没有本质差异。显然，此解只反映了沿Γ_jl法向流场的变化情况。

(ⅱ) 顶点处二维黎曼问题

这是真正的多维问题，可以写成如下形式^[27]:

$ \begin{array}{c} \partial_{t} \boldsymbol{u}+\partial_{x} \boldsymbol{F}(\boldsymbol{u})+\partial_{y} \boldsymbol{G}(\boldsymbol{u})=0, \\ (x, y) \in \mathbb{R}^{2}, t>0 \\ \boldsymbol{u}(x, y, 0)=\boldsymbol{u}_{\ell}, \quad(x, y) \in \varTheta_{\ell} \end{array} $

(40)

这里Θ_ℓ，ℓ=1, …, K, 是从原点出发的角形区域, 见图 1。由于有限传播性质，从中心点发出的复杂波结构只会影响有限的区域。仅从通量的构造来说，顶点处解对界面通量的贡献占比很小, 可以适当限制Courant数，减少顶点周边流场对数值通量的影响，从而可以忽略此问题的求解。在移动网格方法中，特别是拉格朗日方法中，需要用顶点解确定顶点运动速度，于是顶点处二维黎曼问题(40)的解非常重要。但由于解的结构过于复杂，倒不如通过顶点近似黎曼解法器来处理更为方便可靠^[28]。

黎曼问题的(逼近)解被广泛用到双曲守恒律，并延展到更一般的流体力学方程组半离散数值方法中。与下面广义黎曼问题的解进行比较可以发现，这个解不能充分反映出流体的动力学过程；即使从微观角度，欧拉方程反映了粒子无穷碰撞的结果。再者，在每个控制单元上及初始时间层，流场处于常状态，空间的变化依赖相邻单元之间的间断来实现。因此，Godunov格式的解不能充分反映瞬时行为。对长时间、多物理以及多尺度现象，数值模拟结果常常出现似是而非的现象。

到目前为止黎曼问题只限于对双曲守恒律进行研究, 相关的拓展可以从下面的广义黎曼问题研究中看出。

2.5 高阶数值通量与广义黎曼问题

有限体积方法中高阶数值通量的构造实质上等价于相应广义黎曼问题的数值求解，即广义黎曼问题解法器以及高阶数值实现(算法构造)。即在t=t_n时刻，给定初始数据P_n^k∈V^k，求解方程组(8)。与黎曼问题类似，广义黎曼问题包括界面法向广义黎曼问题以及顶点广义黎曼问题。后者只是在依赖网格顶点速度的移动网格方法中起着重要作用，有兴趣的读者可参阅文献[28]。下面只讨论前者。

一般地，方程组(8)的广义黎曼问题可以写成如下形式:

$ \begin{array}{c} \partial_{t} \boldsymbol{u}+\nabla \cdot \boldsymbol{F}(\boldsymbol{u}, \nabla \boldsymbol{u}, \cdots)=0,\quad \boldsymbol{x} \in \mathbb{R}^{n}, \quad t>0 \\ \boldsymbol{u}(\boldsymbol{x}, \boldsymbol{t}=0)=\left\{\begin{array}{l} \boldsymbol{u}_{L}(\boldsymbol{x}), & \varPhi(\boldsymbol{x})<0 \\ \boldsymbol{u}_{R}(\boldsymbol{x}), & \varPhi(\boldsymbol{x})>0 \end{array} \right. \end{array} $

(41)

其中Φ(x)=0是定向的n-1维光滑曲面，经过适当的坐标变换，可记该曲面为x=0 (记空间坐标为x=(x, y, z))，即式(41)可化为如下形式:

$ \begin{array}{c} \partial_{t} \boldsymbol{u}+\nabla \cdot \boldsymbol{\tilde F}(\boldsymbol{u}, \nabla \boldsymbol{u}, \cdots) =\boldsymbol{\tilde H}(\boldsymbol{x}, \boldsymbol{u}), \boldsymbol{x} \in \mathbb{R}^{n}, t>0 \\ \boldsymbol{u}(\boldsymbol{x}, \boldsymbol{t}=0) =\left\{\begin{array}{l} \boldsymbol{u}_{L}(\boldsymbol{x}), & \boldsymbol{x}<0 \\ \boldsymbol{u}_{R}(\boldsymbol{x}), & \boldsymbol{x}>0 \end{array} \right. \end{array} $

(42)

其中x=0称为界面。这里未知量仍采用同样的记号，由变换后所得形式$ \mathit{\boldsymbol{\tilde H}} $(x, u)来自于曲面Φ(x)=0的几何效果，下面讨论中将省略记号“~”。可以看出，广义黎曼问题与相应黎曼问题至少有三方面不同：(ⅰ) 初始数据不同，u_L(x)与u_R(x)是非常状态的两个光滑函数。除了精度以外，数据的非一致性包含更多的物理信息，比如熵变化。(ⅱ) 黎曼问题只对双曲守恒律而言，而广义黎曼问题的控制方程可以针对任意的偏微分方程，比如具有黏性力、外力等效应的流体力学方程组, 或更一般的Boltzmann型方程。(ⅲ) 正如在式(41)所表达那样，当Φ(x)=0不是一个平面时，它的拓扑变化也嵌入到解的变化中。

根据有限体积方法基本原理，需要直接逼近:

$ \mathcal{F}\left(A_{\varepsilon} ; 0, \delta\right)=\int_{A_{\varepsilon}} \int_{0}^{\delta} \boldsymbol{F}(\boldsymbol{u}(0, y, z ; t)) \mathrm{d} y \mathrm{~d} z \mathrm{~d} t $

(43)

这里A_ε={(0, y, z); y∈(y₀-ε, y₀+ε), z∈(z₀-ε, z₀+ε)}是面x=0上的一部分。在一维情况下，(43)即为:

$ \mathcal{F}(0 ; 0, \delta t)=\int_{0}^{\delta t} \boldsymbol{F}(\boldsymbol{u}(0 ; t)) \mathrm{d} t $

(44)

根据已有流体力学方程组相关的偏微分方程理论与应用研究，可以有效地逼近或求解式(42)，满足与式(27)对应的相容性要求。具体地，式(44)有:

$ \begin{array}{c} \mathcal{F}(0 ; 0, \delta t)= \\ \delta t\left[\boldsymbol{F}\left(\boldsymbol{u}\left(0 ; 0^{+}\right)\right)+\frac{\delta t}{2} \partial_{t} \boldsymbol{F}\left(\boldsymbol{u}\left(0 ; 0^{+}\right)\right)\right]+\mathrm{O}\left(\delta t^{3}\right) \end{array} $

(45)

这里需要被积函数F(u(0;t))满足关于时间t的正则性。对于式(42)中的初始数据，这一点可以得到保障，从而可以对式(44)进行相容性逼近。

定义2.3 (有限体积方法的时空耦合性)。考虑时空关联模型(8)的有限体积方法(25)。如果方程(8)的解可以用来逼近数值通量，并使式(25)在定义2.1的意义下相容, 数据重构也充分利用式(8)的时空关联信息，则算法(25)是时空耦合的。

注意文献[22]给出的Hermite型数据重构方法是时空耦合的，该方法已用在两步四阶方法^{[4, 29]}中。文献中数据重构的时空耦合性研究还不充分，上述定义中“数据重构也充分利用式(8)的时空关联信息”这句话比较模糊，需要做更深入的研究。

2.6 几个时空耦合算法的例子

下面给出几个时空耦合算法的例子。

2.6.1 线性输运方程

对于线性方程(1), 其有限体积格式可以写成:

$ \begin{aligned} \bar{u}_{j}^{n+1}=& \bar{u}_{j}^{n}-\lambda_{j}^{n}\left[\frac{1}{\Delta t} \int_{t_{n}}^{t_{n+1}} u\left(x_{j+\frac{1}{2}}, t\right) \mathrm{d} t-\right.\\ &\left.\frac{1}{\Delta t} \int_{t_{n}}^{t_{n+1}} u\left(x_{j-\frac{1}{2}}, t\right) \mathrm{d} t\right], \lambda=\frac{a \Delta t_{n}}{\Delta x_{j}} \end{aligned} $

(46)

假设a>0, λ_jⁿ < 1。如果初始数据是分片常数:

$ u\left(x, t_{n}\right)=\bar{u}_{j}^{n}, x \in I_{j} $

(47)

则解在(t_n, t_n+1)时间内有u(x_{${j + \frac{1}{2}}$}, t)=u_jⁿ，及:

$ \frac{1}{\Delta t} \int_{t_{n}}^{t_{n+1}} u\left(x_{j+\frac{1}{2}}, t\right) \mathrm{d} t=\bar{u}_{j}^{n} $

(48)

从而由(46)可得:

$ \bar{u}_{j}^{n+1}=\bar{u}_{j}^{n}-\lambda_{j}^{n}\left(\bar{u}_{j}^{n}-\bar{u}_{j-1}^{n}\right). $

(49)

这就是迎风格式。数值通量的构造完全使用了解的时空关联信息，因此迎风格式(49)是唯一的一阶时空耦合格式。众所周知，迎风格式在所有一阶稳定格式中具有“最优”性质。

如果初始数据是分片多项式，特别是分片线性函数时:

$ u\left(x, t_{n}\right)=\bar{u}_{j}^{n}+\sigma_{j}^{n}\left(x-x_{j}\right), x \in I_{j}, $

(50)

则式(1)的解式(2)可写为:

$ \begin{array}{c} u\left(x_{j+\frac{1}{2}}, t\right)=\bar{u}_{j}^{n}+\sigma_{j}^{n}\left[\frac{\Delta x}{2}-a\left(t-t_{n}\right)\right], \\ t \in\left(t_{n}, t_{n+1}\right) \end{array} $

(51)

直接计算可得:

$ \frac{1}{\Delta t_{n}} \int_{t_{n}}^{t_{n+1}} u\left(x_{j+\frac{1}{2}}, t\right) \mathrm{d} t=u\left(x_{j+\frac{1}{2}}, t_{n}+\frac{\Delta t}{2}\right)=: u_{j+\frac{1}{2}}^{n+\frac{1}{2}}, $

(52)

以及

$ \left\{\begin{array}{l} \bar{u}_{j}^{n+1}=\bar{u}_{j}^{n}-\lambda_{j}^{n}\left(u_{j+\frac{1}{2}}^{n+\frac{1}{2}}-u_{j-\frac{1}{2}}^{n+\frac{1}{2}}\right), \quad \lambda_{j}^{n}=\frac{a \Delta t_{n}}{\Delta x_{j}} \\ \sigma_{j}^{n+1}=\frac{1}{\Delta x_{j}}\left[u\left(x_{j+\frac{1}{2}}, t_{n+1}\right)-u\left(x_{j-\frac{1}{2}}, t_{n+1}\right)\right] \end{array}\right. $

(53)

注意梯度更新σ_jⁿ⁺¹也使用了解(51)的时空关联信息，所以格式(53)是时空耦合的。这里梯度σ_jⁿ⁺¹更新后得到的分片线性函数仍可能有振荡现象，并不能有效逼近间断函数，所以数据重构仍然值得深入探讨。

更一般地，我们可以利用文献[29, 22]的方式来构造高阶时空耦合方法，或者Lax-Wendroff型的单步高阶方法，只是推广到实际工程问题时，真正非线性和多维性会给单步方法带来实质性的困难。

2.6.2 线性对流扩散方程

考虑下列方程:

$ \partial_{t} u+a \partial_{x} u=\mu \partial_{x}^{2} u, \quad t>0 $

(54)

这里物理黏性系数μ>0。把它写成散度形式有:

$ \partial_{t} u+\partial_{x}\left(a u-\mu \partial_{x} u\right)=0, \quad \mu>0 $

(55)

在时空控制体上，把它表示为式(32)的形式，进行通量逼近。文献[8]中的守恒元/解元方法展示其时空耦合的属性。由于其构造需要一点篇幅展示，有兴趣的读者可以查阅文献[8]及其后来的一系列文章。

相对来说, Navier-Stokes方程组及其相关时空关联模型的时空耦合算法研究较少. 尽管形式上可以进行简单的时空对换处理，但对耗散项等高阶空间导数项的处理仍需要严格数学论证。当然，GKS方法间接提供了Navier-Stokes方程的数值算法，具有时空耦合性^{[9, 30]}。

2.6.3 线性松弛系统

考虑下列简单的松弛系统:

$ \begin{array}{c} \partial_{t} u+a \partial_{x} u=\frac{v-u}{\tau}=: \frac{Q(u, v)}{\tau}, t>0 \\ u(x, 0)=u_{0}(x) \end{array} $

(56)

其中τ>0是松弛参数，a为常数，v也是常数。这个方程的有限体积形式为：

$ \begin{array}{c} \bar{u}_{j}^{n+1}=\bar{u}_{j}^{n}-\lambda\left(u_{j+\frac{1}{2}}^{n+\frac{1}{2}}-u_{j-\frac{1}{2}}^{n+\frac{1}{2}}\right)+ \\ \frac{1}{\Delta x} \int_{t_{n}}^{t_{n+1}} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} \frac{Q(v, u)}{\tau} \mathrm{d} x \mathrm{~d} t \end{array} $

(57)

为了构造时空耦合算法，可以用式(56)的解表达式:

$ u(x,t) = [{u_0}(x - at) - v]{{\rm{e}}^{ - \frac{t}{\tau }}} + v $

(58)

来计算数值积分u_{${j + \frac{1}{2}}$}^{$n + \frac{1}{2}$}以及$\frac{{Q\left({v, u} \right)}}{\tau }$的积分。为了更一般的应用，这里使用文献[30]中的方法，即DUGKS方法，得到:

$ \bar{u}_{j}^{n+1}=\bar{u}_{j}^{n}-\lambda\left(u_{j+\frac{1}{2}}^{n+\frac{1}{2}}-u_{j-\frac{1}{2}}^{n+\frac{1}{2}}\right)+\frac{\Delta t}{2 \tau}\left(Q_{j}^{n}+Q_{j}^{n+1}\right) $

(59)

其中u_{${j + \frac{1}{2}}$}^{$n + \frac{1}{2}$}是由下列形式给出:

$ \begin{array}{c} u_{j+\frac{1}{2}}^{n+\frac{1}{2}}=u\left(x_{j+\frac{1}{2}}-\frac{a}{2} \Delta t, t_{n}\right)+ \\ \frac{1}{2 \tau}\left[Q\left(x_{j+\frac{1}{2}}-\frac{a}{2} \Delta t, t_{n}\right)+Q_{j+\frac{1}{2}}^{n+\frac{1}{2}}\right] \end{array} $

(60)

进而可以利用式(58)或使用式(60)类似的想法，得到u(x_{${j + \frac{1}{2}}$}, t_n+1)，并用来更新梯度

$ \sigma _j^{n + 1} = \frac{1}{{\Delta {x_j}}}[u({x_{j + \frac{1}{2}}},{t_{n + 1}}) - u({x_{j - \frac{1}{2}}},{t_{n + 1}})] $

(61)

很显然，格式(59)~式(61)是时空耦合的，并具有时空二阶精度。

3 基于广义黎曼解法器(GRP solver)的时空耦合算法

数值通量的构造是执行有限体积格式的两个核心步骤之一。对于非线性问题我们一般无法像上述线性方程那样, 得到解的表达式。下面将要利用广义黎曼问题(41)或(42)解的局部正则性质，发展广义黎曼解法器。其主要思想可以概括为：

为了方便叙述，这一节统一把界面放在x=0, 把广义黎曼问题的求解点平移到坐标原点(0, 0), 初始数据表示成式(42)的形式。记:

$ {\mathit{\boldsymbol{u}}^*}: = \mathop {\lim }\limits_{t \to 0 + } \mathit{\boldsymbol{u}}({\bf{0}},t),{({\partial _t}\mathit{\boldsymbol{u}})^*}: = \mathop {\lim }\limits_{t \to 0 + } {\partial _t}\mathit{\boldsymbol{u}}({\bf{0}},t). $

(62)

广义黎曼问题(GRP)解法器：给定控制方程以及形如式(42)的分片光滑函数作为初始数据，求解式(42) 并得到形如式(62)的极限值。

一旦有了极限值(62)，受益于解u(x, t)在界面上关于时间t的连续性，我们有:

$ \mathit{\boldsymbol{u}}(0,\delta t) = {\mathit{\boldsymbol{u}}^*} + {({\partial _t}\mathit{\boldsymbol{u}})^*}\delta t + {\rm{O}}(\delta {t^2}) $

(63)

由于u^*只表示了一种瞬时行为，其解由相应的黎曼解给定:

$ {\mathit{\boldsymbol{u}}^*} = R(0;{\mathit{\boldsymbol{u}}_L}(0),{\mathit{\boldsymbol{u}}_R}(0)) $

(64)

接下来的关键任务是求值(∂_tu)^*。“非常不严格”地由控制方程(42)直接得到:

$ {\partial _t}\mathit{\boldsymbol{u}} = - \nabla \cdot \mathit{\boldsymbol{\tilde F}}(\mathit{\boldsymbol{u}},\nabla \mathit{\boldsymbol{u}}, \cdots ) + \mathit{\boldsymbol{\tilde H}}(\mathit{\boldsymbol{x}},\mathit{\boldsymbol{u}}) $

(65)

然后求得极限值。这样解u(x, t)的变化率可以通过空间的变化反映，就是经典的Lax-Wendroff方法^[23]，或Cauchy-Kowalevski方法^[31]。通过这种途径把模型的时空关联性质嵌入到数值格式中，所得的算法具有时空耦合性质。

需要指出的是, 一般的单步Lax-Wendroff方法需要计算解的高阶时间导数∂_t^ku，但对流体力学方程组来说，需要避免这样的操作，原因是: (ⅰ) 由于解中常有间断，这样的运算失去了数学与物理意义；(ⅱ) 即使允许如此运算，所得出的通量估算过于复杂，失去了实际工程意义；(ⅲ) 最近发展的两步四阶方法^{[4, 29]}表明，这一对值具有明确物理意义，可以用来作为构造高阶方法的基石，避免了高阶时间导数的运算。

3.1 线性GRP解法器

对于线性系统:

$ \begin{aligned} \partial_{t} \boldsymbol{u}+\boldsymbol{A} \partial_{x} \boldsymbol{u} &=\boldsymbol{C} \boldsymbol{u}+\boldsymbol{D},\ \ \ \ \ \ \ \ t>0 \\ \boldsymbol{u}(x, 0) &=\left\{\begin{array}{ll} \boldsymbol{u}_{L}(x), & & x<0 \\ \boldsymbol{u}_{R}(x), & & x>0 \end{array}\right. \end{aligned} $

(66)

其中A、C、D是常数矩阵，A可以对角化，Λ=R^-1AR, 其中Λ是由A的特征值λ_i组成的矩阵，Λ=diag{λ_i}。记w=R^-1u, 有:

$ \begin{array}{l} \partial_{t} \boldsymbol{w}+\boldsymbol{\varLambda}\partial_{x} \boldsymbol{w}=\boldsymbol{R}^{-1} \boldsymbol{C} \boldsymbol{u}+\boldsymbol{R}^{-1} \boldsymbol{D}, \ \ t>0 \\ \boldsymbol{w}(\boldsymbol{x}, 0)=\left\{\begin{array}{ll} \boldsymbol{w}_{L}(x), \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ & x<0 \\ \boldsymbol{w}_{R}(\boldsymbol{x}), \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ & x>0 \end{array}\right. \end{array} $

(67)

进而记Λ⁺=diag{max(λ_i, 0)}，Λ^-=diag{min(λ_i, 0)}，I⁺=$\frac{1}{2}$diag{1+sign(λ_i)}, I^-=-$\frac{1}{2}$diag{1-sign(λ_i)}。显然我们可以得到：

$ \begin{array}{l} ({\partial _t}{\mathit{\boldsymbol{w}}^*} = - ({\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^ + }{\kern 1pt} \mathit{\boldsymbol{w}}_{^L}^\prime + {\mathit{\boldsymbol{ \boldsymbol{\varLambda} }}^ - }{\kern 1pt} \mathit{\boldsymbol{w}}_{^R}^\prime ) + \\ \left( {{\mathit{\boldsymbol{I}}^ + }{\kern 1pt} {\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{u}}_L} + {\mathit{\boldsymbol{I}}^ - }{\kern 1pt} {\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{u}}_R}} \right) + {\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{D}} \end{array} $

(68)

其中u_L=u_L(0), u_R=u_R(0), w′_L=w′_L(0), w′_R=w′_R(0)。再记A^±=RΛ^±R^-1, 回到原始变量u, 我们有：

$ \begin{array}{L} \left(\partial_{t} \boldsymbol{u}\right)^{*}=-\left(\boldsymbol{A}^{+}\ \boldsymbol{u}_{L}^{\prime}+\boldsymbol{A}^{-}\ \boldsymbol{u}_{R}^{\prime}\right)+ \\ \ \ \ \ \left(\boldsymbol{R} \boldsymbol{I}^{+}\ \boldsymbol{R}^{-1} \boldsymbol{C} \boldsymbol{u}_{L}+\boldsymbol{R} \boldsymbol{I}^{-}\ \boldsymbol{R}^{-1} \boldsymbol{C} \boldsymbol{u}_{R}\right)+\boldsymbol{D} \end{array} $

(69)

特别当u_L=u_R=u^*时有:

$ \left(\partial_{t} \boldsymbol{u}\right)^{*}=-\left(\boldsymbol{A}^{+}\ \boldsymbol{u}_{L}^{\prime}+\boldsymbol{A}^{-}\ \boldsymbol{u}_{R}^{\prime}\right)+\boldsymbol{C} \boldsymbol{u}^{*}+\boldsymbol{D} $

(70)

从中可以看出如何应用Lax-Wendroff解法器和间断追踪计算(∂_tu)^*。

多维情形是类似的。例如考虑二维线性方程组：

$ \begin{array}{ll} \partial_{t} \boldsymbol{u}+\boldsymbol{A} \partial_{x} \boldsymbol{u}+\boldsymbol{B}\partial_{y} \boldsymbol{u}=\boldsymbol{C} \boldsymbol{u}+\boldsymbol{D}, \ \ t>0 \\ \boldsymbol{u}(x, y, \boldsymbol{0})=\left\{\begin{array}{ll} \boldsymbol{u}_{L}(x, y), \ \ \ \ \ \ \ \ \ \ \ \ & x<0 \\ \boldsymbol{u}_{R}(x, y), \ \ \ \ \ \ \ \ \ \ \ \ & x>0 \end{array}\right. \end{array} $

(71)

这里A、B、C、D是常矩阵，且A、B可对角化。切向变化量B∂_yu可看作相对法向的一个扰动, 将线性方程组写成如下形式：

$ \begin{array}{ll} \partial_{t} \boldsymbol{u}+\boldsymbol{A} \partial_{x} \boldsymbol{u}=-\boldsymbol{B}\partial_{y} \boldsymbol{u}+\boldsymbol{C} \boldsymbol{u}+\boldsymbol{D}, \ \ t>0 \\ \boldsymbol{u}(x, y, 0)=\left\{\begin{array}{ll} \boldsymbol{u}_{L}(x, y), \ \ \ \ \ \ \ \ \ \ \ \ & x<0 \\ \boldsymbol{u}_{R}(x, y), \ \ \ \ \ \ \ \ \ \ \ \ & x>0 \end{array}\right. \end{array} $

(72)

与上面类似方法，我们可以得到:

$ \begin{array}{c} \left(\partial_{t} \boldsymbol{u}\right)^{*}=-\left[\boldsymbol{A}^{+}\ \left(\partial_{x} \boldsymbol{u}\right)_{L}+\boldsymbol{A} ^{-}\ \left(\partial_{x} \boldsymbol{u}\right)_{R}\right]- \\ \left.\boldsymbol{R} \boldsymbol{I}^{+}\ \boldsymbol{R}^{-1}\ \boldsymbol{B}\left(\partial_{y} \boldsymbol{u}\right)_{L}+\boldsymbol{R} \boldsymbol{I}^{-}\ \boldsymbol{R}^{-1} \boldsymbol{B}\left(\partial_{y} \boldsymbol{u}\right)_{R}\right]+ \\ \left(\boldsymbol{R} \boldsymbol{I}^{+} \ \boldsymbol{R}^{-1} \boldsymbol{C} \boldsymbol{u}_{L}+\boldsymbol{R} \boldsymbol{I}^{-}\ \boldsymbol{R}^{-1} \boldsymbol{C} \boldsymbol{u}_{R}\right)+\boldsymbol{D} \end{array} $

(73)

这里记号与上述意思相同，并有(∂_xu)_L=$\mathop {\lim }\limits_{x \to {0^ - }}$ ∂_xu(x, y), (∂_xu)_R=$\mathop {\lim }\limits_{x \to {0^ + }} $∂_xu(x, y)等。从式(70)和式(73)可以看出，所有的信息都是时空关联的，最终的算法是迎风且时空耦合的。

特别强调，通过GRP解法器把切向变化自然地蕴含到数值通量中，这是法向(一维)黎曼解法器做不到的^[38]。换句话说，如果只使用黎曼问题解法器，导致格式失去多维性，该格式应用到多维问题模拟时自然会有数值缺陷，不得不用其他方式进行补偿。再者，GRP解法器是保证数值格式真正多维性的有效手段，在合理的数据投影基础上，算法具有涡量保持等性质。

3.2 声波近似——线性化GRP解法器

对于非线性系统来说，当流场是光滑的或者波的强度不大时，使用近似解法器进行通量的逼近就可以取得不错的效果，即声波近似或线性化GRP解法器。Toro等的ADER解法器就是GRP的声波近似版本^[35]。

先看下面的一维双曲平衡律:

$ \begin{array}{c} \partial_{t} \boldsymbol{u}+\partial_{x} \boldsymbol{F}(\boldsymbol{u})=\boldsymbol{H}(x, \boldsymbol{u}), \quad x \in \mathbb{R}, t>0 \\ \mathit{\boldsymbol{u}}(x,0) = \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{u}}_L}(x),}&{x < 0}\\ {{\mathit{\boldsymbol{u}}_R}(x),}&{x > 0} \end{array}} \right. \end{array} $

(74)

可以把它看作描写相对于网格界面法向效应的方程。如果u_L(0-0)=u_R(0+0), 而u′_L(0-0)≠u′_R(0+0), 只有线性波(声波)从(0, 0)点发出。记A(u)=∂_uF(u), u^*=u_L(0-0)=u_R(0+0), u=u^*+v, v可看成u在u^*附近的扰动。然后在u=u^*处线性化，由式(73)可得:

$ \partial_{t} \boldsymbol{v}+\boldsymbol{A}\left(\boldsymbol{u}^{*}\right) \partial_{x} \boldsymbol{v}=\boldsymbol{H}\left(x, \boldsymbol{u}^{*}\right), \quad \boldsymbol{u}=\boldsymbol{u}^{*}+\boldsymbol{v} $

(75)

这是一个线性系统。类似于式(66)，从中可以计算(∂_tu)=(∂_tv)^*,

$ \left(\partial_{t} \boldsymbol{v}\right)^{*}=-\boldsymbol{A}^{+}\ \left(\boldsymbol{u}^{*}\right) \boldsymbol{v}_{L}^{\prime}(0)-\boldsymbol{A}^{-}\ \left(\boldsymbol{u}^{*}\right) \boldsymbol{v}_{R}^{\prime}(0)+\boldsymbol{H}\left(0, \boldsymbol{u}^{*}\right) $

(76)

这里A^±(u^*)=R(u^*)Λ^±(u^*)R^-1(u^*), Λ^±(u^*)是由Λ(u^*)的“±”部分组成Λ⁺=diag{max(λ_i), 0)}，Λ^-=diag{min(λ_i), 0)}。从而

$ \left(\partial_{t} \boldsymbol{u}\right)^{*}=-\boldsymbol{A}^{+}\ \left(\boldsymbol{u}^{*}\right) \boldsymbol{u}_{L}^{\prime}(0)-\boldsymbol{A}^{-}\ \left(\boldsymbol{u}^{*}\right) \boldsymbol{u}_{R}^{\prime}(0)+\boldsymbol{H}\left(x, \boldsymbol{u}^{*}\right) $

(77)

这就是线性化GRP解法器(当源项H(x, u)存在时，即使u′_L=u′_R≡0, (∂_tu)^*≠0，意味着GRP解法器仍然起着重要作用)。

当u_L(0-0)≠u_R(0+0)时，声波近似策略如下：以局部黎曼解u^*=R(0;u_L(0), u_R(0))为背景解, 线性化式(74)得到式(75)，从而可得线性化GRP解法器式(77), 具体见文献[32]。

对于多维系统(仍然以二维为例):

$ \begin{array}{c} \partial_{t} \boldsymbol{u}+\partial_{x} \boldsymbol{F}(\boldsymbol{u})+\partial_{y} \boldsymbol{G}(\boldsymbol{u})=\boldsymbol{H}(x, y, \boldsymbol{u}), \quad t>0 \\ \boldsymbol{u}(x, y, 0)=\left\{\begin{array}{ll} \boldsymbol{u}_{L}(x, y), & x<0 \\ \boldsymbol{u}_{R}(x, y), & x>0 \end{array}\right. \end{array} $

(78)

采用声波近似的策略，线性化这个方程组得到:

$ \begin{array}{c} \partial_{t} \boldsymbol{v}+\boldsymbol{A}\left(\boldsymbol{u}^{*}\right) \partial_{x} \boldsymbol{v}=-\boldsymbol{B}\left(\boldsymbol{u}^{*}\right) \partial_{y} \boldsymbol{v}+\boldsymbol{H}\left(x, y, \boldsymbol{u}^{*}\right), t>0 \\ \boldsymbol{v}(x, y, 0)=\left\{\begin{array}{ll} \boldsymbol{v}_{L}(x, y), & x<0 \\ \boldsymbol{v}_{R}(x, y), & x>0 \end{array}\right. \end{array} $

(79)

从而得到类于式(73)的解法器, 这里B(u^*)=∂_uG(u)。

3.3 真正非线性GRP解法器

在式(74)中，如果初始数据在原点(相对网格尺寸)有“很大”跳跃，就会出现强非线性波。线性化GRP解法器不足以分辨强非线性波，需要使用真正非线性GRP解法器，否则就会出现明显的数值缺陷^{[7, 34-35]}。一般的界定标准为：

$ \left\| {{\mathit{\boldsymbol{u}}_L}(0) - {\mathit{\boldsymbol{u}}_R}(0)} \right\|{\rm{ }} \gg \alpha \Delta x $

(80)

参数α>0是一个重要的经验参数。求解GRP解法器的基本思想是：解析强稀疏波，追踪强间断。特别需要指出，在强间断的情形下，热力学效应起着重要作用，真正非线性的GRP解法器反映了热力学相容的性质^[7]。

本文不具体给出真正非线性GRP解法器。对于可压缩欧拉方程组，可参阅文献[5, 36]。后期发展的不依赖于坐标系统的GRP解法器，详见文献[6]。对一般的双曲方程组，可参见文献[37]。

对于多维情形，仍然可以把切向的变化和间断面的变化看作扰动，考虑有下列形式问题:

$ \begin{array}{c} \partial_{t} \boldsymbol{u}+\partial_{x} \boldsymbol{F}(\boldsymbol{u})=-\partial_{y} \boldsymbol{G}(\boldsymbol{u})+\boldsymbol{H}(x, y, \boldsymbol{u}), t>0 \\ \boldsymbol{u}(x, y, 0)=\left\{\begin{array}{ll} \boldsymbol{u}_{L}(x, y), & x<0 \\ \boldsymbol{u}_{R}(x, y), & x>0 \end{array}\right. \end{array} $

(81)

从而利用3.1和3.2节中法向解法器求解式(72)和式(73)，这主要源自于一个关键的事实：切向扰动在法向的传播是连续的。由此切向的变化在通量中得到充分反映，得到的数值方法具有真正的多维性^[38]。

正如前面所讨论的那样，除非涉及(依赖于网格移动)中心拉格朗日型数值方法，这里不关注顶点多维黎曼解法器。由于真正多维黎曼问题的理论基础很不成熟^{[27, 39-41]}，真正多维顶点黎曼解法器常常带来不稳定的数值结果^[42]。在实践中，人们更倾向使用健壮的逼近GRP解法器, 例如，Maire等利用新的框架并结合拉格朗日型GRP解法器^[43]，发展了稳定的中心型高阶拉氏方法。

3.4 一般时空关联模型的GRP解法器

对于一般的时空关联模型，如多相流、多介质模型^[44]和Navier-Stokes方程等，广义黎曼问题解法器可以进行类似的构造，简述如下。

对于多相流、多介质模型，广义黎曼问题解法器可以归结为一般的非线性平衡律的范畴，模型的时空关联性质在GRP解法器中可以得到充分体现^[45-47]。困难在于介质组份以及混合引起的数值振荡、物质分数为负等非物理解，涉及物理建模本身的缺陷以及有限体积格式的投影(平均)过程。由于热力学关系的高度非线性，投影过程不能反映精确的物理过程, 例如内能的平均过程导致物质界面附近的压力振荡. 因此，在投影步实现时空耦合，充分利用解的信息也许可以提高相关数值质量。

对于Navier-Stokes方程组等宏观模型，基本的想法是类似的。对于对流项(欧拉部分)，使用上述标准的广义黎曼解法器。需要指出的是：初始数据式(41)需要进行适当的匹配。也就是说该初始数据至少是二阶以上分片多项式。另外可使用的策略为，对于对流占优问题，采用松弛策略，把相应模型双曲化^[48-49]。这样的策略是基于能量原理——在对流占优情形下，能量集中于初始能量传播影响区域内(双曲性质)。

从Boltzmann方程直接出发的动理学方法是设计流体力学数值方法的一条有效途径。一般来说，很难写出Boltzmann方程的解，但在特定的近似假定下，如BGK模型^[50]，可以类似于式(58)那样隐式得出解的表达式，并将之用于数值通量的构造，得出的数值方法如气体动理学格式(GKS)^[9]、统一气体动理学格式UGKS^[10]等。特别地，GKS可以用来模拟Navier-Stokes方程描述的流动，可以作为Navier-Stokes的解法器来使用。按照这样的定义，在通量的构造部分，GKS和UGKS解法器显然是时空耦合的。

3.5 高精度方法中黎曼解法器和GRP解法器的比较

在高精度数值方法中，黎曼解法器和GRP解法器都可用来构造数值通量。前者在高阶线方法中常用，如Runge-Kutta型高阶方法，后者用在两步四阶方法中。下面仍以一维双曲守恒律方程组(31)来比较两种解法器的差别。

考察控制体$\left[{{x_{j- \frac{1}{2}}}, {x_{j + \frac{1}{2}}}} \right] \times [{t_n}, {t_{n + 1}})$边界上通量。在x=x_{${j + \frac{1}{2}}$}处, 记:

$ \begin{array}{l} \boldsymbol{u}_{j+\frac{1}{2}}^{n}=\lim\limits _{t \rightarrow t_{n}+} \boldsymbol{u}\left(x_{j+\frac{1}{2}}, t\right) \\ \left(\partial_{t} \boldsymbol{u}\right)_{j+\frac{1}{2}}^{n}=\lim\limits _{t \rightarrow t_{n}+} \partial_{t} \boldsymbol{u}\left(x_{j+\frac{1}{2}}, t\right) \end{array} $

(82)

由解关于时间t的正则性得到:

$\begin{array}{c} \boldsymbol{F}\left(\boldsymbol{u}\left(x_{j+\frac{1}{2}}, t\right)\right)=\boldsymbol{F}\left(\boldsymbol{u}_{j+\frac{1}{2}}^{n}\right)+\\ \boldsymbol{F}^{\prime}\left(\boldsymbol{u}_{j+\frac{1}{2}}^{n}\right)\left(\partial_{t} \boldsymbol{u}\right)_{j+\frac{1}{2}}^{n}\left(t-t_{n}\right)+\mathrm{O}\left(\left(t-t_{n}\right)^{2}\right) \end{array} $

(83)

(ⅰ) 黎曼解法器。选取F_{${j + \frac{1}{2}}$}^Godu=F(u_{${j + \frac{1}{2}}$}ⁿ), 从而有:

$ \begin{array}{c} \int_{t_{n}}^{t_{n+1}}\left[\boldsymbol{F}\left(\boldsymbol{u}\left(x_{j+\frac{1}{2}}, t\right)\right)-\boldsymbol{F}\left(\boldsymbol{u}\left(x_{j-\frac{1}{2}}, t\right)\right)\right] \mathrm{d} t \\ \quad-\Delta t_{n}\left[\boldsymbol{F}_{j+\frac{1}{2}}^{\text {Godu }}-\boldsymbol{F}_{j-\frac{1}{2}}^{\text {Godu }}\right] \\ =\frac{\left(\Delta t_{n}\right)^{2}}{2}\left(\boldsymbol{F}^{\prime}\left(\boldsymbol{u}_{j+\frac{1}{2}}^{n}\right)\left(\partial_{t} \boldsymbol{u}\right)_{j+\frac{1}{2}}^{n}-\right. \\ \left.\boldsymbol{F}^{\prime}\left(\boldsymbol{u}_{j-\frac{1}{2}}^{n}\right)\left(\partial_{t} \boldsymbol{u}\right)_{j-\frac{1}{2}}^{n}\right)+\mathrm{O}\left(\left(\Delta t_{n}\right)^{3}\right) \end{array} $

(84)

所以对于间断解来说，通量的截断误差为:

$ E_j^n = {\rm{O}}({(\Delta {t_n})^2}) $

(85)

误差阶q=0。不过，对于光滑解来说，式(84)中的差赋予了一个“阶数”，

$ E_j^n = {\rm{O}}\left( {{{\left( {\Delta {t_n}} \right)}^3}} \right) $

(86)

从而误差阶为q=1。

(ⅱ) GRP解法器。选取$\mathit{\boldsymbol{F}}_{j + \frac{1}{2}}^\text{GRP} = \mathit{\boldsymbol{F}}\left({\mathit{\boldsymbol{u}}_{j + \frac{1}{2}}^n} \right) + \frac{{{{(\Delta {t_n})}^2}}}{2}\mathit{\boldsymbol{F'}}\left({\mathit{\boldsymbol{u}}_{j + \frac{1}{2}}^n} \right)({\partial _t}\mathit{\boldsymbol{u}})_{j + \frac{1}{2}}^n$, 从而有:

$ \begin{array}{c} \int_{t_{n}}^{t_{n+1}}\left[\boldsymbol{F}\left(\boldsymbol{u}\left(x_{j+\frac{1}{2}}, t\right)\right)-\boldsymbol{F}\left(\boldsymbol{u}\left(x_{j-\frac{1}{2}}, t\right)\right)\right] \mathrm{d} t \\ \qquad=\Delta t_{n}\left[\boldsymbol{F}_{j+\frac{1}{2}}^{\mathrm{GRP}}-\boldsymbol{F}_{j-\frac{1}{2}}^{\mathrm{GRP}}\right] \\ =\frac{\left(\Delta t_{n}\right)^{3}}{6}\left(\partial_{t}^{2} \boldsymbol{F}\left(\boldsymbol{u}_{j+\frac{1}{2}}^{n}\right)-\partial_{t}^{2} \boldsymbol{F}\left(\boldsymbol{u}_{j-\frac{1}{2}}^{n}\right)\right)+O\left(\left(\Delta t_{n}\right)^{4}\right) \end{array} $

(87)

与上面讨论类似，对于间断解，GRP解法器有一阶精度q=1，而对于光滑解有二阶精度q=2。

3.6 时空耦合数值边界条件

边界条件的数值处理是计算流体力学中的一个核心问题，大部分的数值处理基于对虚拟区域的延拓(外插技巧)或特征传播理论^[51]。最近逆Lax-Wendroff方法^[52]用来数值处理边界条件，这种方法蕴含了时空耦合性。我们认识到，流体力学方程组的数值边界条件等价于单边广义黎曼问题的数值求解^[53]。事实上，先前的工作已经认识到单边黎曼问题在数值边界条件处理上的重要性^[54-55]。这些工作与直接的逆Lax-Wendroff方法相比，有以下特点：

(ⅰ) 在有限体积框架下，计算区域的边界自然为控制体的边界，所以数值边界条件处理等价于边界上的通量数值逼近，它归结为单边黎曼解法器的构造。

(ⅱ) 单边黎曼解法器中的(∂_tu)^*其实反映了边界条件上的受力情况，即牛顿第二定律的数学表现，这是时空耦合算法的体现。在文献[53]中，这一事实是单边黎曼解法器的基石。

(ⅲ) 在技术上，如此处理数值边界条件可以尽量少使用虚拟网格。在时空二阶格式中无需使用虚拟网格；相比较其它更高阶方法，至多使用一半的虚拟网格。

(ⅳ) 无需使用更高阶导数，避免了人为的复杂性和虚假物理性质等^[52]。

单边黎曼解法器的使用将极大简化边界条件的数值处理，特别是解决无结构网格下虚拟区域的插值问题^[56]。

3.7 高阶时空耦合算法

高阶精度数值方法对小尺度问题的数值模拟非常重要，如湍流等。基于有限体积格式的框架式(25)，高阶方法需要在数值通量和数据重构两部分具有高阶逼近性质。在通量的逼近部分，可以采用Lax-Wendroff方法，但流体力学方程组(8)的非线性性质导致高阶展开过于复杂，缺乏实际工程价值。更糟糕的是，解的间断性质意味着高阶展开的运算没有意义, 因此需要寻求其他方式。在过去的几十年中，各种空间重构技术以及Runge-Kutta型时间推进方法取得了巨大的进步，可以容易地查阅到相关文献。

最近，文献[29]发展了两步四阶时间推进方法，成功地融合了Runge-Kutta型方法简单性优点和基于Lar-Wendroff型解法器的时空耦合性质。文献[4]详述了该类方法的特性：(ⅰ) 计算效率。同等条件下其计算效率是Runge-Kutta方法的一半(二维情形为例)^[57]。(ⅱ) 紧致性。由于时间推进步骤减少一半，模版就会减少一半；时空耦合的HWENO型重构^{[22, 58]}可以使计算格式更加紧致。(ⅲ) 稳定性。已经证明其稳定区域比Runge-Kutta大^[59]。(ⅳ) 兼容性。该方法很容易和其它方法相结合，如GKS解法器^[60]、多矩方法^[61]以及CRP重构技术^[62]等。

目前这类方法还局限在“显式”框架内，相关的研究处于起步阶段。为了应用的需要，亟须发展“隐式”两步四阶方法，以适应“强刚性”等多尺度问题的求解。

3.8 投影过程时空耦合性的一点注释

高阶精度数值方法的投影过程常称为数据重构，原因在于有限体积的发展过程$\mathcal{E}$一般只给出守恒量的平均值u_jⁿ⁺¹，通过它们得到V^k中的函数P_n+1^k，WENO型重构过程一般采用这种策略^[63]。撇开很多细节，WENO型重构模版太大，基于非结构网格的重构面临很多困难，后来的HWENO重构有效改善了模版太大的问题^[64]。

本文更愿意称V^k中函数P_n+1^k的构造过程为投影过程，它可通过时空关联解u(x, t_n+1)来实现，从而所得到的算法具有“完全的”时空耦合性。GRP解法器提供了这种可能性，它提供的边界点处的极限值式(62)不仅可以用来构造通量，还可以用在投影过程中^{[5-6, 37]}:

$ \boldsymbol{u}_{j+\frac{1}{2}}^{n+1}=\boldsymbol{u}_{j+\frac{1}{2}}^{n}+\Delta t\left(\partial_{t} \boldsymbol{u}\right)_{j+\frac{1}{2}}^{n} $

(88)

从而有:

$ \boldsymbol{\sigma}_{j}^{n+1}=\frac{1}{\Delta x_{j}}\left(\boldsymbol{u}_{j+\frac{1}{2}}^{n+1}-\boldsymbol{u}_{j-\frac{1}{2}}^{n+1}\right) $

(89)

再次强调，注意极限值(62)已用在时间区间通量的构造中，无需额外的处理。由于(∂_tu)_{${j + \frac{1}{2}}$}ⁿ被赋予了时空关联信息，如此构造的梯度具有时空耦合特性，这与前文2.6节中几个线性模型的例子是一致的。

在最近的两步四阶方法中，我们也使用了这一策略^[22]，得到的数据更加紧致，从而格式也具有紧致性。

4 数值算例

在可压缩流体的算法设计及数值模拟中，有大量基准算例。随着算法进步和工程需求的提高, 基准算例应该与时俱进，以面对相当极端的环境。在文献[65]中，一系列基准算例值得用新的数值算法进行测试。下面几个算例，可以看出时空耦合性的重要性。

4.1 大压力比(大密度比)问题

这个问题又称为LeBlanc问题，它是经典的激波管问题的极端情形，从中可以看出数值方法对强稀疏波的分辨能力以及它对激波产生的影响。控制方程为一维欧拉方程(31)，初始数据具有如下形式：

$ (\rho, u, p)(x, 0)=\left\{\begin{array}{ll} \left(10^{4}, 0,10^{4}\right), & 0 \leqslant x<0.3 \\ (1,0,1), & 3 \leqslant x \leqslant 1.0 \end{array}\right. $

(90)

多方指数取1.4，计算的终止时间t=0.12。文献[66]利用这个例子检测了多个数值格式的性能。这里我们同样用这个算例针对性讨论本文涉及的一些观点，数值结果来源于文献[7, 29]。

首先在图 2中展现使用不同解法器的一阶格式。可以看出一阶Godunov方法随着网格加密，可以收敛到精确解；使用逼近解法器，收敛相对较慢，但仍然具有收敛性。图 3使用了空间二阶重构，而时间离散使用一阶向前欧拉方法，解法器分别使用黎曼解法器和逼近的黎曼解法器(HLLC, Roe), 数值表现很差，不具有收敛性。正如在3.5节所述, 当空间具有高阶精度，使用一阶黎曼解法器构造数值通量等，算法不具有相容性；类似的情形反映在图 4, 即使时间离散使用二阶Runge-Kutta方法。

图 5考察了在这种极端情况下使用线性化方法的数值表现，黎曼解使用声波近似解法器。其实文献[34-35]已经指出线性化解法器的数值缺陷。图 6使用了非线性GRP解法器，数值结果立即改善，说明了强非线性出现后真正非线性GRP解法器的重要性。图 7展示的结果是基于GRP解法器的两步四阶方法^[29]，从中看到用100网格点可以很好分辨间断，300网格点可以得到完美效果。这是许多数值方法很难达到的，从而说明了时空耦合算法的必要性。

4.2 激波和熵波相互作用的问题

这个算例是Shu-Osher算例的推广，考虑了更极端的情形，又称为Titarev-Toro问题^[67]。控制方程依然为欧拉方程，初始数据为:

$ (\rho, u, p)(x, 0)=\left\{\begin{array}{l} (1.515696,0.52336,1.805), \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -5 \leqslant x<-4.5 \\ (1+\sin (20 \pi x), 0,1), \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ -4.5 \leqslant x \leqslant 5 \end{array}\right. $

(91)

这里使用了GKS解法器^[9]得到图 8的数值结果，详细说明见文献[60]。

4.3 激波和悬浮刚体的相互作用

这个数值结果取自文献[54]，模拟了隧道[0, 1]×[0, 0.2]内激波撞击一个长a=0.06、宽b=0.03矩形刚体的流场情况。控制方程为二维欧拉方程组，我们使用了时空耦合二阶GRP方法，并用单边广义黎曼解法器刻画方块移动时的数值边界条件。初始时刻一个马赫数3的激波位于x=0.08处，波前状态为(ρ₀, u₀, v₀, p₀)=(1.3, 0, 0, 0.1)；该方块的密度为ρ_Ω=10ρ₀, 质量为M_Ω=abρ_Ω。初始时刻该方块以倾斜度$\frac{{\rm{ \mathsf{ π} }}}{4}$被放置于隧道中，重心为(0.15，0.06)，惯性力为${A_\mathit{\Omega} } = \frac{{{M_\mathit{\Omega} }}}{3}({a^2} + {b^2})$。图 9展示了两个时刻的流场压力分布情况。

4.4 多介质激波和气泡的相互作用

这个例子展示了激波和气泡相互作用的情况(图 10)，该问题已经成为多相流领域一个标准算例^[68]。下面的结果取自文献[45]，分别用能量分裂的Godunov方法和GRP方法模拟所得结果(图 11)。显然，GRP很好地提高了流场的分辨率。

4.5 各向同性可压缩湍流的模拟

下面的算例采用了两步四阶GKS方法^[60]模拟的超音速各向同性可压缩湍流的情况^[69](图 12)，从中可查阅具体的计算参数，其中计算区域是[－π, π]³。该文展示了两步四阶方法和二阶方法同样健壮。对许多高阶方法来说，分辨激波碎片(shocklet)是个巨大的挑战^[70]，该方法中GKS解法器的时刻耦合性质继承了物理本身的属性。

5 讨论与展望

流体力学的时空关联模型刻画了物理量空间变化在时间方向上的演化的传播，如各种波的传播等。当进行数值模拟时，这种性质理应得到继承，相应地所用算法应该具备时空耦合特性。实践过程中这个看似简单课题理应得到关注。遗憾的是，相对于时空解耦方法，时空耦合算法需要更好理解模型的物理性质或数学理论，人们自然“避难就易”。一个直白的问题是——即使物理问题的数学模型十分完美，当相应的算法缺乏相应完美性质时，其数值结果怎么令人信服呢？

本文首先严格论述有限体积方法。正如引言所解释的原因：(ⅰ) 无论从物理定律的形成还是数值算法的构造，有限体积方法是解流体力学问题最自然的框架，本文总结了这个框架形成的数学理论和内在涵义；(ⅱ) 流体现象的复杂性(如间断等)客观地阻碍了有限体积方法严格数学理论的形成，对于这一框架的认识常出现诸多似是而非的争论; (ⅲ) 对于许多实际工程问题，基于无结构网格的算法设计是一个自然要求，在此之上许多方法相互冲突^[24]，有必要从最基本的地方出发建立一个根本准则。幸运地是，许多工程软件，如Fluent, 自然选择在有限体积框架下生成；许多工程人员当遇到难以跨越的困难时，往往借用有限体积框架度过难关。从论证过程可以看到时空耦合是算法的根本属性。

有限体积格式的步骤很简单：投影过程和数值通量的构造。目前CFD算法的大部分研究者将注意力集中于前者，基本上在空间逼近论的范畴内进行探索。尽管黎曼不变量等物理量以及其他技术用于重构过程，但是在随机选取方法中重要的时空耦合性质相当缺乏。黎曼问题的解被用在随机选取中，在某种意义下，它可以看成是一种时空耦合投影方法。本文侧重于后者的讨论，给出了有限体积格式与积分平衡律(12)之间的相容性准则。特别对于数值通量介绍了黎曼解法器和GRP解法器，并在3.5节中进行了对比。简单地总结如下，具体内容可以到文献[2]中查阅。

(ⅰ) 关于Riemann解法器。对于双曲守恒律(欧拉方程)，当初始数据是分片常数时，由Riemann解法器产生的数值通量具有无穷阶相容性，即Godunov格式就是积分平衡律；当初始数据是非常数的分片光滑函数时，Riemann解法器产生的数值通量对光滑解是一阶相容的，但对间断解是不相容的(见3.5节)。也就是说对于高阶空间重构，如果不能有匹配的数值通量，产生的数值格式是不相容的。值得注意的是，如果控制方程包含外力项时，Godunov格式永远是不相容的。

(ⅱ) 关于GRP解法器。GRP解法器给出时空耦合的数值通量。对于光滑流场，GRP解法器得到二阶相容的数值通量；但对于间断解只有一阶时间精度。GRP解法器是保证有限体积格式的逼近解收敛的一个必要条件。

这里之所以再次强调这点并细致给出第算例4.1节的, 原因在于为了看清它与传统理解上的差异。事实上，算法时空耦合性也体现在其他方面，比如最近研究气体动理学的渐近性质时，提出了统一保持性质(unified preserving property, UP) 的概念^[71]。对于一个动理学方法，不仅应该具有欧拉层面的渐近保持性质(Asymptotic preserving property, AP)，还应该根据需要具有Navier-Stokes或更深层面的UP性质，这个过程实际上是通过时空耦合的思想实现的。该文分别用时空耦合DUGKS方法以及时空解耦CLR格式，对二维不可压缩Taylor涡进行数值模拟，见图 13，具体见文献[71]。特别指出的是，这一问题本质上是多尺度问题。在多尺度问题及其相关的研究中，时空耦合策略是一条有效途径。

由于篇幅限制，本文只给出了有限体积方法基本原理的相容性准则以及通过GRP解法器实现相容性的技术细节，对很多根本性质还没有涉及，如时空相容格式的熵稳定性等。对于可压缩流体力学来说，熵稳定性对应于热力学相容性^[7]，这部分留在将来的文章中探讨。

客观地说，时空耦合算法的研究非常不成熟，目前只在相对比较“纯粹”的模型上进行分析和验证，但是这一思想应该加以推广与应用，比如如何将时空耦合算法的思想应用于一般的时空关联湍流模型^[72]、粒子建模和时空耦合算法等。就算法本身来说，时空耦合的隐式格式研究还没有开展，这一领域的研究应该是下一阶段的重点。

在工程应用中时空耦合的程序显得更少，一方面是受限算法模块的历史延续性制约，另一方面是工程领域内的“理性”力学越来越少。加强工程算法的科学性是需要目前亟待解决的观念问题。

后记。作为一个数学工作者，很忐忑地接受邀请在力学的专业期刊《空气动力学学报》奉献此类稿件，毕竟隔行如隔山。但想到数学、力学本就是“一家”，向力学家“学习”本就是“数学”的一大研究动机，心里稍微坦然点。

致谢: 此文的想法是在与下列合作者合作过程中逐步形成的，在此一并致谢，他们是Matania Ben-Artzi、齐进、汪玥、杜知方、雷昕、徐昆、郭照立、李启兵，等。特别感谢曹贵瑜博士提供的关于各向同性可压缩湍流模拟算例4.5, 感谢徐昆教授、郭照立教授和汪玥副研究员提出许多根本性改进意见，感谢陈海波给予了文字上的润色。