A spacetime outlook on CFD: Spacetime correlated models and spacetime coupled algorithms
IAPCM ,北京 年 月 日 计算流体力学的时空观 : 模型的时空关联性与算法的时空耦合性 计算流体力学的时空观:模型的时空关联性及算法的时空耦合性 李杰权 北京应用物理与计算数学研究所计算物理重点实验室; 北京大学应用物理与技术中心 摘要: 流体力学中波的有限传播、粒子的碰撞、各种力之间相互作用,无不体现时空关联效应。本文以计算方法的视角探讨计算流体力学的时空观,即流体力学模型的时空关联性,计算方法的时空耦合性。从流体力学微团法建模出发,明确模型时空关联性的涵义,建立有限体积格式的基本原理,阐述算法时空耦合的必要性,实现流体力学基本控制方程物理建模与有限体积格式数学原理的统一。在实践中,给出时空耦合高精度数值方法设计思路,通过算例比较它与时空解耦方法的差别。需要指出,本文的大部分内容适用于连续介质假定下的计算流体力学研究,部分仅适用于可压缩流动。 ⼀ 引⾔ 四维空间听来好象有些神秘,其实早已有之,即以“宇宙“二字来说,‘往古来今谓之宙,四方上下谓之宇”(《淮南于·齐俗训》)就是宇是东西、南北、上下三维扩展的空间,而宙是一维的时间。 接着他引用相对论和生活的常识叙述了时空既独立又统一的性质: 爱因斯坦不再把“宇”、“宙”分开来看,也就是时间也在进行着。每一瞬间三维空间中的物质在占有它一定的位置。他根据麦克斯韦一洛伦兹的光速不变假定,并继承了牛顿的相对性原理而提出了狭义相对论。狭义相对论中的洛伦兹变换把时空联系在一起,当然并不是消灭了时空特点。如向东走三里,再向西走三里,就回到原处,但时间则不然,共用了走六里的时间。时间是一去不复返地流逝着。 华先生意在说明数学在描述“宇宙之大”的作用,本文的引用旨在强调 时空之间的关系 深刻地植根于计算流体力学的模型理解和算法构造。离开了时空演化,一切归于“稳态”。正因为时空演化,流体的世界才是色彩斑澜,波澜壮阔。 定性甚至定量描述流体世界的时空演化并不容易,数值模拟同样是貌似简单,实则不易。下面是众所周知的例子, 𝜕 ! 𝑢(𝑥, 𝑡) + 𝑎 𝜕 " 𝑢(𝑥, 𝑡) = 0, (1) 基金项目 :国家自然科学基金(11771054, 91852207),国家重大项目 GJXM92579, 计算物理重点实验室基金 作者简介 :李杰权( ),男,安徽,研究员,研究方向:计算流体力学、偏微分方程,数值分析 通信作者 :李杰权, Email : [email protected]. 其中 𝑎 是一个常数, 𝑥 表示空间坐标, 𝑡 是时间变量。此模型表示了变量 𝑢(𝑥, 𝑡) 的空间变化(扰动、波)与时间变化(传播)的关系,描述了一个时空关联的性质 𝑢(𝑥, 𝑡 + Δ𝑡) = 𝑢(𝑥 − 𝑎Δ𝑡, 𝑡). (2) 此表达式事实上比(1)更加根本和直接,并不牵涉函数 𝑢(𝑥, 𝑡) 本身更多性质,比如正则性. 即使 𝑢(𝑥, 𝑡) 是一个方波或脉冲, (2)仍然有效。用平移算子 𝑒 ! 来表示(2),则有 𝑢(𝑥, 𝑡 + Δ𝑡) = 𝑢(𝑥 − 𝑎Δ𝑡, 𝑡) = 𝑒 ! 𝑢(𝑥, 𝑡). (3) 从这里出发,如果想精确表达(1)的输运性质,则需要进行展开(如泰勒展开) 𝑢(𝑥, 𝑡 + Δ𝑡) = 𝑢(𝑥 − 𝑎Δ𝑡, 𝑡) = 𝑒 ! 𝑢(𝑥, 𝑡) = 𝑢(𝑥, 𝑡) − 𝑎Δ𝑡 𝜕 " 𝑢(𝑥, 𝑡) + (−𝑎Δ𝑡 𝜕 " ) ( (4) 为了设计实际工程中的计算方法(简称算法),需进行必要操作:一是截断处理,涉及算法与控制方程的相容性,涉及精度的概念;二是对空间变化 𝜕 " 𝑢(𝑥, 𝑡) 的(如有限差分)处理,不同的处理得到不同的算法。这些都与函数 𝑢(𝑥, 𝑡) 的正则性和模型(1)的内在性质息息相关:正则性越好,精度越高; 模型(1)的时空关联性需要通过(2)精确描述. 这些处理充满技巧和无限可能,怎样“令人信服”地设计时空耦合算法且精准描述相应模型的内在时空关联性,是一个自然的、但要深入思考的本质问题。 近年来,计算流体力学得到蓬勃的发展,各种算法和相应的应用软件目不暇接。有些算法具有坚实的理论基础,如基于变分原理的有限元方法等。而流体力学中关于可压缩流动的研究,由于流场结构非常复杂(激波、物质界面、涡团和湍流等),人们对流动(解的)性质缺乏很好理解,模型的适定性和相关数值模拟中算法的探讨相当工程化,形式化的表达往往导致似是而非的结果,常常是模型具有很好的性质,而相应数值模拟的离散模型失去了该保持的性质,比如模型的时空关联性和算法的时空耦合性等。不同的出发点导出的数值算法往往截然不同。本文试图在这方面做一点尝试: 基 于 连 续 介 质 假 定 探 讨 流 体 力 学 模 型 时 空 变 化 的 内 在 关 联 性 ( inherent spacetime correlation ),阐述相应算法应保持的时空耦合( spacetime coupling )性质 。 本文从流体力学微团法的直接建模开始思考模型的时空关联性。通过对通量的研究,发现经典的瞬时 Cauchy 通量无法反映时空关联性,提出了某一时间段内总通量( lump flux ),在任意特定时间段内反映流体的动力学过程【2,3】,进而提出 有限体积格式的基本原理:在连续介质假定性下,任一时间段内流体通过控制体边界的总通量关于边界的扰动是 Lipschitz 连续的 。这一 Lipschitz 连续性恰好体现流场时空关联性的本质特征,实现了数学上流体力学方程组弱解和物理上积分型平衡律( integral balance laws) 之间的统一,后者其实就是全离散有限体积格式的物理起源。通过对流体力学方程组时空关联性质的研究,实现有限体积框架下的时空耦合。在算法的实现过程中, 物理量(质量、动量、能量等)以及它的变化率 这一对量( pair )起着重要的作 用,一个量的变化率通过它的空间变化关系来反映。变化率不是一类常规的数学演算,而是反映了内在的物理规律。例如,加速度是速度的变化率,它等于控制体表面受到的力,即应力,这是牛顿第二定律。在流体力学计算方法中使用这样的量是自然的也是必须的,反映了时空变化的联系。在文【4】,笔者倡导的时空耦合高精度算法正是这一思想的具体表现。本文再次描述实现时空耦合算法的一种基本流程,并通过算例展示时空耦合算法的数值表现。 这里倡导的计算流体力学时空观不是新的,从偏微分方程的
Cauchy-Kowalevski 方法到
Lax-Wendroff 方法就是某种意义下的时空耦合方法. 现代计算流体力学的算法时空观常有体现,如对流方程的迎风格式、广义黎曼解法器【5,6,7】、守恒元/解元(
CE/SE )方法【8】、气体动理学解法器【9,10】等。最近几年,笔者致力于思考高精度时空耦合算法的必要性与理论基础,但远远不够深入. 随着计算技术的提高和工程需求的精细化,这方面的研究应该得到重视。故作此文,抛砖引玉。 二 流体力学方程组的时空关联性 众所周知,流体力学建模过程常用微团法,即在某一时间段 (𝑡, 𝑡 + 𝛿𝑡) , 研究控制体 (𝒙, 𝒙 + 𝛿𝒙) 内流体质团的运动。时间间隔 𝛿𝑡 非常重要,给出了微团运动的动力学过程响应时间。为了数学的描述方便,基于连续介质假定,流体力学方程组写成如下形式:给定一个控制体 Ω(𝑡) ,流体运动满足下列的关系 𝑑𝑑𝑡 7 𝒖(𝒙, 𝑡) 𝑑𝒙 = − 7 𝒇 ⋅ 𝒏 𝑑Γ ’)(!))(,) (5) 这里 𝒖 是密度函数, 𝒏 是 𝜕Ω(𝑡) 的单位外法向。相应地,称 𝒎(𝑡) = 7 𝒖(𝒙, 𝑡) 𝑑𝒙 )(!) (6) 为控制体 Ω(𝑡) 上的质量,右边
𝒞(𝜕Ω; 𝑡) = − 7 𝒇 ⋅ 𝒏 𝑑Γ ’)(!) (7) 为
Cauchy 通量( flux ), 𝒇 为通量密度函数,简称为通量函数,它是 𝒖 以及空间变化量 𝛁𝒖 的函数 𝑓 = 𝒇(𝒖, 𝛁𝒖, ⋯ ) 。方程组(5)反映了质量时间变化率与通过控制体边界 𝜕Ω(𝑡) 的通量之间的瞬时关系。 Cauchy 在特定连续假定下,论述了
Cauchy 通量的性质【11,12】,这是他对连续介质力学(
Continuum Mechanics )最重要的贡献【11】。要深入理解这一关系并不容易,这与
Gauss-Green 定理密切相关。基于这一定理,可以将(5)写成散度型偏微分方程(
PDE )形式 略去外力场等效应. 略去雷诺输运定理等细致讨论,并设 Ω 不随流场运动,即欧拉型控制体。相应地,(1)中的控制体为拉格朗日控制体。 𝜕 ! 𝒖 + ∇ ⋅ 𝑭(𝒖, ∇𝒖, ⋯ ) = 0. (8) 这个方程对光滑流动是有效的,可由偏微分方程的知识建立时空的关联性。对于实际流动来说,这一转化是非常粗糙的,难点在于 Cauchy 通量 𝒞(𝜕Ω; t) 的数学性质,即正则性。可压缩流场蕴涵丰富的现象,如冲击波、物质界面、湍流等效应,流场(或称方程的解)的正则性非常弱,人们对它的认识很少,故而该定理的应用并不显然。历史上有很多研究【13,14,15,16】,直到最近这方面的研究工作还是热点【17】,遗憾的是所取得的结果与实际的期望仍然相差甚远。正如最近的专著【18 , Section “ 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 𝑞 - may be computed from A ” , 其中 𝑞 - 对应于这里的 𝒇(𝒖) ⋅ 𝒏 , 𝐴 对应于一个给定的应力张量( stress tensor) 。也就是说: 人们还不知道如何从应力张量中计算出通量函数。 仔细审视可以看出 Cauchy 通量 𝓒(𝝏𝛀; 𝒕) 是表示一个特定瞬时 𝑡 的行为,方程(5)只是表示一个非常弱的“时-空”关联性质。现在常用的一个观点是下面的“弱解”概念,它来自于偏微分方程(8): 定 义 1.1. 设 Ω 是 ℝ . 中 的 一 个 有 界 区 域 , 𝑛 = 1,2,3, / ≤ 𝑡 ( . 对 任 意 试 验 函 数 𝜙(𝑥, 𝑡) ∈ 𝐶 (Ω × [𝑡 / , 𝑡 ( ]) , 如果函数 𝑢(𝑥, 𝑡) 满足 ! 𝜙 + 𝑭 ⋅ ∇ϕ]𝑑𝒙𝑑𝑡 = 0 )! " ! , (9) 则称 𝒖(𝒙, 𝑡) 是(8)的弱解。 这个表达式比(5)前进了一大步:如果解 𝒖(𝒙, 𝑡) 是光滑的,(8)和(9)是等价的。近年来偏微分方程理论得到极大的发展,人们对于(8)的认识比直接对(5)的认识要深刻得多,且偏微分方程(8)的时空关联性一目了然: 流场的空间摄动可以通过(8)反映到时间的变化中, 线性波的传播(1)正是这种时空关联性的例子。另外,当一个间断面在无限薄的假定下,可以导出 Rankine-Hugoniot 间断关系,即间断面(线)两侧解的“ 迹” 的关系。 定义1.1并没有对函数 𝒖(𝑥, 𝑡) 做过多假设,只要(9)成立即可。一旦流场出现复杂间断或其它流场结构时,这个定义并不能给出我们太多的信息。我们回避不了的事实是:在可压缩流动中,流场即使初始性质非常“好”,但在不久的将来由于复杂的非线 性 相 互 作 用 , 可 能 变 得 非 常 “ 糟 糕 ” 。 这 可 以 从 最 简 单 的 Burgers 方 程 看 到【19,20】。因此,人们不得不面对复杂的情形,深入思考为什么应用
Gauss-Green 定理理解
Cauchy 通量非常困难。 任何流动都有一个动力学过程 (dynamic process) , 这也许是个常识, 但需要严格论证。最近的研究表明【3,4】,流动的时空关系既是相互独立又是彼此关联。下面的原理深刻地反映这一事实。 有限体积方法基本原理. 设 𝒖(𝐱, t) 是定义1.1意义下方程(8)的弱解, Ω 是 ℝ 中任一紧集。假定 (i) 𝐮(𝐱, t) 是局部有界; (ii) 质量 𝐦(t) 是时间 t 的连续函数。 那么有下列结论: (i) 对任意 δ > 0, 函数 𝑯(𝒙; 𝑡, 𝑡 + 𝛿𝑡) = ∫ 𝑭(𝒙, t) !34,, dt 满足 ∇ 𝒙 ⋅ 𝑯(𝒙; 𝑡, 𝑡 + 𝛿𝑡) ∈ 𝐿 (ℝ . ). (10) (ii) 对任意 δ > 0 ,时间段( 𝑡, 𝑡 + 𝛿𝑡) 内流过 Ω 边界 𝜕Ω 的总通量( lump flux ) 𝓕(𝛛𝛀; 𝑡, 𝑡 + 𝛿𝑡) = 7 7 𝑭 ⋅ 𝒏 𝒅𝜞𝒅𝒕 𝛛𝛀!3;!! (11) 关于 Ω 的边界 𝜕Ω 扰动是 Lipschitz 连续的。 总通量
𝓕(𝛛𝛀; 𝑡, 𝑡 + 𝛿𝑡) 和 Cauchy 通量 𝓒d𝝏𝛀 ; 𝒕e 明显不同,特别是 Lipschitz 连续性,进一步反映了流场的时空关联性质 。 而 Cauchy 通量只是一个瞬时行为,不能完全反映动力学的过程。这些反映了它们数学性质的根本差异。直观地来说, 这里通过时间效应,换取空间的正则性。 由此给出 (8) 的弱解另一种表述,它恰恰是物理上常见的 积分型平衡律。 定义1.2. 设 Ω 是 ℝ . 中任一紧集, 𝛿𝑡 > 0 . 如果它满足下列条件: (i) 𝒖(𝒙, 𝑡) 局部有界,且质量 𝒎(𝑡) 是时间 𝑡 的连续函数; (ii) 整体通量(4)有意义且关于 Ω 的边界 𝜕Ω 扰动是连续的; (iii) 积分平衡律成立 𝛀 = − 7 7 𝑭 ⋅ 𝒏 𝑑Γ 𝑑𝑡 ’)!3;!!) , (12) 其中 𝒏 是 𝜕Ω 的单位外法向, 则称函数 𝒖(𝒙, 𝑡) 是偏微分方程(8)的弱解。 事实上,已经证明定义1.1和1.2是等价的, 见【3,4】,并且给出了通量(11)的正则性质。这个定义说明将要讨论的有限体积方法就是计算流体力学的基本格式,和物理最原始的建模一致。流体力学方程组(8)可理解为 :如果流场光滑,用偏微分方程组(8)刻画;但当流场失去正则性时,解需要满足积分平衡律(12)。 流体力学有限体积格式的构造过程就是基于(12)的直接建模过程: 在偏微分方程(8)诱导出的解辅助下,构造出和(12)相容的离散模型,精准反映真实的物理过程。 三 有限体积方法及其时空耦合性 前面论述到:有限体积方法构造过程就是在偏微分方程(8)辅助下的一个直接建模过程。通常的有限体积方法从偏微分方程组(8)出发,在每一个计算单元(控制体)上应用 Gauss-Green 定理,得到有限体积公式。具体地,将计算区域 Ω 剖分为 𝑁 个计算单 元(控制体) Ω < , Ω = ⋃ Ω <=<>/ , 𝜕Ω ? = ⋃ Γ <ℓℓ . 有下列两种形式: 半离散和全离散有限体积方法。 在 Ω ? 上,对 (8) 的空间散度项应用 Gauss-Green 公式可得 半离散有限体积方程 𝑑𝑑𝑡 7 𝒖(𝒙, 𝑡) 𝑑𝒙 = − 7 𝑭 ⋅ 𝒏 𝑑Γ , 𝑗 = 1, ⋯ , 𝑁. ’) $ ) $ (13) 它直接对应于流体力学方程组 (5) 。在初始时刻 𝑡 = 0 , 给定 𝒖(𝒙, 𝑡) 的分布, 𝒖(𝒙, 0) = 𝑷 (𝑥) ∈ 𝒱 C , (14) 其中 𝒱 C 是一个可持续函数空间(通常为分片多项式函数空间 ),意味着在每个时间层 𝑡 = 𝑡 . > 0 , 𝒖(𝒙, 𝑡 . ) 经过投影后所得 𝑃 .C ∈ 𝒱 C 中。由于守恒性的基本要求 𝒏 ) 𝛀 𝒋 𝑑𝒙 = 7 𝑷 𝒏𝒌 (𝒙)𝒅𝒙 𝛀 𝒋 , (15) 𝑷 𝒏𝒌 (𝑥) 一般是沿着每个控制体 Ω ? 的边界 Γ <ℓ 是间断的。为了利用(13)更新解,需要在每个控制体 Ω ? 的边界 Γ <ℓ 附近求解(8), 数值上使用后面所描述的黎曼解法器或近似黎曼解法器, 得到在时刻 𝑡 的数值通量 l𝚪 𝒋𝓵 l 𝑭 𝒋𝓵 (𝑡) ⋅ 𝒏 𝒋𝓵 ≈ 7 𝑭 ⋅ 𝒏 𝒅𝚪 𝚪 𝐣𝓵 , (16) 带入(13)即得 𝑑𝑑𝑡 7 𝒖(𝒙, 𝑡) 𝑑𝒙 = − 7 𝑭 ⋅ 𝒏 𝑑Γ ≈ − ol𝚪 𝒋𝓵 l 𝑭 𝒋𝓵 (𝑡) ⋅ 𝒏 𝒋𝓵ℓ , 𝑗 = 1, ⋯ , 𝑁. ’) $ ) $ (17) 接着使用常微分方程求解器,如典型的 Runge-Kutta 方法,求解这个时间依赖(常微分)方程。需要指出的是:这个方程需要在每个时间层 𝑡 = 𝑡 . 求解数值通量,并把(13)在下一时间层的解投影到 𝒱 C 中,得到 𝑃 .3/C ∈ 𝒱 C 。如果用每个时间步用多级 Runge-Kutta 方法,则需要多次求解数值通量和多次投影。这个投影过程又称为 数据重构 。用算子形式,半离散方法可以表示为: 𝑷 .3/C (𝒙) = 𝒫 ∘ ℛ𝒦 ∘ 𝒜 HIJ [𝑷 .C (𝒙)], (18) 其中 𝒜 HIJ 表示用黎曼解法器构造瞬时数值通量, ℛ𝒦 代表 Runge-Kutta 型的时间推进过程, 𝒫 表示数据重构。 用高阶多项式等光滑函数去逼近正则性很弱的函数(如间断函数),常常会导致振荡现象(Gibbs 现象). 数据重构(数据投影)仍是有限体积格式中未彻底解决的难题。 我们注意到,由于对 Cauchy 通量 ∫ 𝑭 ⋅ 𝒏 𝒅𝚪 𝚪 𝐣𝓵 还没有清晰的认识,不得不借助局部黎曼解法器(见下节)进行数值通量(16)的逼近。一般来说,所得到的局部误差为 l𝚪 𝒋𝓵 l 𝑭 𝒋𝓵 (𝑡) ⋅ 𝒏 𝒋𝓵 − 7 𝑭 ⋅ 𝒏 𝒅𝚪 = 𝓞(|(𝚫𝒖) 𝒋𝒍 |) 𝚪 𝐣𝓵 l𝚪 𝒋𝓵 l, (19) 这里 |(𝚫𝒖) 𝒋𝒍 | 表示数据 𝑷 .C (𝒙) 跨过控制体边界 Γ <6 的局部跳跃。 半离散有限体积方法属于 线方法 ( Method of Line ),从某种意义上来说是 时空解耦方法。 在光滑流场中,流体力学方程组的
PDE 关系可以直接反映时空关联性;但对间断解来说,(19)中的局部误差的累积给数值模拟带来巨大困难。 本文将把(12)应用到时空控制体 Ω < × [𝑡 . , 𝑡 .3/ ] , 𝑡 .3/ = 𝑡 . + Δ𝑡 . , Δ𝑡 . > 0 为时间步长,得到有限体积格式的全离散形式 .3/ )𝑑𝒙 − 7 𝒖(𝒙, 𝒕 𝒏 )𝒅𝒙 𝛀 𝐣 = − 7 7 𝑭 ⋅ 𝒏 𝑑Γ 𝑑𝑡. ’) ( ! )* ! ) ) ( (20) 也可以把(20)看作是(13)在时间段 (𝑡 . , 𝑡 .3/ ) 上的积分。与半离散情形下求解瞬时通量(16)不同,这里需要利用偏微分方程(8)的时空关联性质逼近 [𝑡 . , 𝑡 .3/ ] 上总通量, ol𝚪 𝒋𝓵 l 𝑭 𝒋ℓ (𝑡 . , 𝑡 .3/ ) ⋅ 𝐧 ?ℓ𝓵 ≈ o 7 7 𝑭 ⋅ 𝒏 𝒅𝚪𝒅𝒕 𝚪 𝒋𝓵 𝒕 𝒏*𝟏 𝒕 𝒏 𝓵 . (21) 将之代入(20),得到有限体积的公式 𝒖y <.3/ = 𝒖y <. − 1lΩ < l ol𝚪 𝒋𝓵 l 𝑭 𝒋𝓵 (𝑡 . , 𝑡 .3/ ) ℓ ⋅ 𝒏 <6 , (22) 这里记逼近解在控制体 Ω < 上的平均值为 𝑢z <. = 1lΩ < l 7 ) ( 𝒖(𝒙, 𝑡 . )𝑑𝒙. (23) 由有限体积方法基本原理,可以得到 l𝚪 𝒋𝓵 l 𝑭 𝒋𝓵 (𝑡 . , 𝑡 .3/ ) − 7 7 𝑭 ⋅ 𝒏 𝒅𝚪𝒅𝒕 𝚪 𝒋𝓵 𝒕 𝒏*𝟏 𝒕 𝒏 = 𝓞(𝚫𝒕 𝒏 ) 𝒒3𝟏 l𝚪 𝒋𝓵 l , (24) 其中 𝑞 > 0 。值得注意的是,(24)中的误差和(19)中的误差非常不同:(19)中的误差是用解的局部跳跃(变差)来测量,而(24)中的误差直接用时间步长(等价于网格尺度)来测量。当流场相对光滑时,这两者是等价的;但当流场中有剧烈间断、脉动时,二 者差别是明显的.即使网格加密,(19)也得不到收敛解 【7】。后文算例 6.1 可以看出数值通量的重要性。 同半离散方法相似,在每个时间层 𝑡 = 𝑡 . 上投影可得到在可持续空间 𝒱 C 中的逼近解 𝑃 .C 。全离散有限体积方法可以用符号表示为: 𝑷 .3/C (𝒙) = 𝒫 ∘ ℰ ∘ 𝒜[𝑷 .C (𝒙)], (25) 这里 𝒜 表示通量的逼近, ℰ 表示解的有限体积发展过程, 𝒫 代表投影过程(数据重构)。一旦 𝒜 被确定, ℰ 是自然的且没有任何误差产生。数据重构部分是重要的,将在后面评注。 所谓相容性,描述的是离散格式与背景方程之间的关系。传统上常常使用 Taylor 展开的方式研究数值格式与微分方程(如(8))之间的相容性,这样的运算只对光滑流场是成立的。当流场比较复杂时,目前大部分的数值分析某种意义上只是启发性的,比如以标量方程为模型建立相关的理论【21】。 现在基于积分平衡律(12)来建立有限体积格式的相容性,即离散格式(18)或(25)须与(12)进行直接比较。由(18)或(25)可知,误差来源于数据投影步 𝒫 以及通量的逼近步 𝒜 HIJ 或 𝒜 。对于前者,虽然有大量的文献存在,数据重构步主要依赖相关的空间信息。文【22】只是做了初步的尝试,把未来时间的数据重构与过去的信息进行了关联。本文假定投影步满足所有的精度要求,只讨论通量的逼近误差。 对半离散格式(18)来说,在每个固定时刻给出的通量误差至多为 ol𝚪 𝒋𝓵 l 𝓵 𝑭 𝒋𝓵 (𝑡) ⋅ 𝒏 𝒋𝓵 − o 7 𝑭 ⋅ 𝒏 𝒅𝚪 = 𝓞d(𝚫𝒖) 𝒋𝒍 e 𝚪 𝐣𝓵 l𝚪 𝒋𝓵 l 𝓵 = 𝒪d(𝚫𝒖) 𝒋 el𝚪 𝒋 l𝑑 𝒋 , (26) 其中 (𝚫𝒖) 𝒋 表示在控制体附近解的局部跳跃, |𝚪 𝒋 | 是 Ω < 边的最大长度, 𝑑 < = 𝑑𝑖𝑎𝑚(Ω < ) 。并且在每个时间层,随着 Runge-Kutta 步的增加,误差也会累积。在(26)中, 𝑑 < 来源 于不同方向上通量差的获利,在 4.5 节可以进一步看到。 对于全离散方法(25),我们讨论每个时间层的数值通量的逼近,并期望有下列的估计 ol𝚪 𝒋𝓵 l 𝓵 𝑭 𝒋𝓵 (𝑡 . , 𝑡 .3/ ) − o 7 7 𝑭 ⋅ 𝒏 𝒅𝚪𝒅𝒕 𝚪 𝒋𝓵 𝒕 𝒏*𝟏 𝒕 𝒏 𝓵 = 𝓞(𝚫𝒕 𝒏 ) 𝒒3𝟐 l𝚪 𝒋 l . (27) 其中 𝑞 > 0 。比较 (26) 和 (27) ,它们有本质的差别,即 (𝚫𝒖) 𝒋 与 𝚫𝒕 𝒏 的差别 。 对光滑解来说它们是等价的。这样我们给出下列的定义: 在标量方程的研究中,网格加密是可以得到收敛解的,原因在于全局变差有界性质 . 到目前为止还没有 例证可证明该性质在流体力学方程组成立。 定义 3.1(有限体积格式的相容性). 设 Ω < 是任一控制体, 𝑗 = 1, ⋯ , 𝑁 , 如果对于某个 𝑞 > 0 , (27)成立,则有限体积格式(25)相容于平衡律(12),并具有 𝑞 阶相容性。 如上所述,相容性概念常常相对与微分方程所言. 如果(8)是双曲守恒律, Lax 和 Wendroff 给出了有限差分方法的相容性,常称为
Lax 相容性【23】。由于
Lax 相容性概念影响深远,有必要进行回顾。 定义 3.2(Lax 相容性【23】)。 考虑双曲守恒律方程组 𝜕 ! 𝒖 + 𝜕 " 𝒇(𝒖) = 0, (28) 其
2𝑝 + 1 点守恒性差分格式为 𝒖 <.3/ = 𝒖 <. − 𝜆 (cid:129)𝒈 <3/(. − 𝒈 < (cid:131) , 𝜆 = 𝛥𝑡𝛥𝑥, (29) 其中 𝛥𝑥 是网格尺寸, 𝛥𝑡 是时间步长, 𝒈 <3 . = 𝒈(𝒖 < , ⋯ , 𝒖 <3P. ) 。如果成立 𝒈(𝒖, ⋯ , 𝒖) = 𝒇(𝒖) , (30) 则称差分方法(29)与双曲守恒律(28)相容。 基于这个定义,建立了影响深远的 Lax-Wendroff 收敛定理,直至今日,该定理仍被广泛应用。但是随着高阶精度格式的发展以及工程应用的需求,这个定理常常被误用,典型表现为: (i) 相容性定义(30)只适用于传统一阶精度数值方法,对于目前广泛使用的高阶精度方法并不使用,特别是高阶有限体积方法等。 (ii) 在此基础上建立的 Lax-Wendroff 定理只适用于一致网格或结构网格,对非结构网格并不适用【24,25】. 数值验证过程中使用的网格加密方法测试收敛阶,需要倍加小心。 定义 3.1 实际上第一次给出高精度有限体积数值格式与积分平衡律之间的相容性【3】。有限体积方法与网格的结构无关,因此适用于无结构网格。特别地,收敛阶测试时一般针对光滑流进行,数据重构部分能够达到应有的精度,通量逼近阶事实上就是收敛阶。但是,当流场含有间断或其它复杂结构时,数据重构仍存在诸多争议,通量的逼近效果往往被忽略。通过以上分析可以看到,如果通量相容性(27)不成立,逼近解不可能收敛。也就是说,(27)是有限体积格式相容性的一个必要条件。 谈到流体力学中的有限体积方法,需要回顾
Godunov 方法【26】。 经过半个多世纪的发展和检验,它已成为现代计算流体力学基石。考虑一维可压缩欧拉方程组 𝜕 ! 𝒖 + 𝜕 " 𝑭(𝒖) = 0, (31) 其中守恒量 𝒖 = (𝜌, 𝜌𝑢, 𝜌𝐸) Q , 𝑭(𝒖) = d𝜌𝑢, 𝜌𝑢 ( + 𝑝, 𝑢(𝜌𝐸 + 𝑝)e Q . 这个方程组由 Gibbs 关 系 式 封 闭 , 这 里 不 再 赘 述 。 记 时 空 控 制 体 𝐶 <. = 𝐼 < × [𝑡 . , 𝑡 .3/ ] , 𝐼 < = (𝑥 I , , 𝑥 <3 ) , 𝑥 < = /( (cid:136)𝑥 < + 𝑥 <3 (cid:137) , Δ𝑥 < = 𝑥 <3 − 𝑥 < , 𝑡 .3/ = 𝑡 . + Δ𝑡 . 。 那么(31)可以写成 " 𝒖(𝑥, 𝑡 ./0 )𝑑𝑥 − " 𝒖(𝑥, 𝑡 . )𝑑𝑥 + " 𝑭 -𝒖 .𝑥 , 𝑡/0 𝑑𝑡 − !" ! $ $ " 𝑭 -𝒖 .𝑥 , 𝑡/0 𝑑𝑡 = 0 !" ! , (32) 即积分平衡律(12)的形式。对应于(22)有限体积格式为 𝒖y <.3/ = 𝒖y <. − Δ𝑡 . Δ𝑥 < (cid:129)𝑭 <3/(𝒏 − 𝑭 𝒋 (cid:131) , (33) 其中 𝑭 <3 𝒏 是数值通量。给定初始的数据分布 𝒖(𝒙, 𝒕 𝒏 ) = 𝑷 𝒏 𝟎 ∈ 𝓥 , Godunov 格式的一个中心思想是在每个网格边界点 (cid:136)𝑥 <3 , 𝑡 . (cid:137) 处求解相应黎曼问题 𝜕 ! 𝒖 + 𝜕 " 𝑭(𝒖) = 0, 𝑥 ∈ ℝ, 𝑡 > 𝑡 . , 𝒖(𝑥, 𝑡 . ) = (cid:139) 𝒖y <. , 𝑥 < 𝑥 <3//( ,𝒖y <3/. , 𝑥 > 𝑥 <3//( . (34) 求得网格边界上解的值 𝒖 𝒋3 𝒏 。欧拉方程组(31)的黎曼解可以清晰给出,记为 𝒖(𝑥, 𝑡) =𝑅d𝜉, 𝒖 𝒋𝒏 , 𝒖 𝒋3𝟏𝒏 e , 𝜉 = " (* ! ) , 𝒖 𝒋3 𝒏 = 𝑹 (𝟎; 𝒖y 𝒋𝒏 , 𝒖y 𝒋3𝟏𝒏 ) 。 这样 𝑭 <3 𝒏 = 𝑭(𝒖 <3 . ) , 从而得到 Godunov 格式, 𝒖y <.3/ = 𝒖y <. − Δ𝑡 . Δ𝑥 < (cid:129)𝑭 (cid:145)𝒖 <3/(. (cid:146) − 𝑭 (cid:145)𝒖 < (cid:146)(cid:131). (35) 这是和(12)完全一致的 . " 𝑭 -𝒖 .𝑥 , 𝑡/0 𝑑𝑡 = 𝑭 - 𝒖 𝑗+12𝑛 . !" ! 因此,从积分平衡律的逼近来说,基于分片常数的初值逼近,(35)与(12)完全相容或有无穷阶相容性。从整体逼近(25)来说,所有的误差来自于投影运算, 因此习惯上称 Godunov 格式为一阶格式。 如果用逼近的黎曼解法解,界面上的值一般可以写为 𝒖(cid:147) .<3/( = 12 d𝒖y <. + 𝒖y <3/. e − 𝛼d𝒖y <. , 𝒖y <3/. e2 (cid:149)𝑭d𝒖y <3/. e − 𝑭d𝒖y <. e(cid:150), (36) 这里 𝛼d𝒖y <. , 𝒖y .<3/ e 称为粘性系数。一般来说,通量的逼近误差为 . " 𝑭 -𝒖 .𝑥 , 𝑡/0 𝑑𝑡 − 𝑭 - 𝒖 𝑛𝑗+12 = 𝒪(|𝒖 𝑗+1𝑛 − 𝒖 𝑗𝑛 |) !" ! . (37) 这个误差在强间断的情形下对计算结果有巨大的伤害。如(19)所述,这个误差不随着网格加密而减小。 在多维情形下,例如二维双曲守恒律情形 𝜕 ! 𝒖 + 𝜕 " 𝐹(𝒖) + 𝜕 Y 𝑮(𝒖) = 0, (38) 在每个时间层 𝑡 = 𝑡 . 的数据为分片常数,即每个控制体 𝛀 𝒋 内为常数 𝒖y 𝒋𝒏 , 通过求解相应的黎曼问题,构造数值通量。这时相应的黎曼问题可分为两类: (i) 相对于界面的法向黎曼问题 由于流体力学方程组的伽利略( Galilean) 不变性,通过旋转变换,总可以设 𝑥 -方向为 Γ <ℓ 的外法向, Γ <ℓ 两侧单元的值记为 𝒖 Z 和 𝒖 [ ,法向黎曼问题定义为 𝝏 𝒕 𝒖 + 𝜕 " 𝑭(𝒖) = 0, (𝑥, 𝑦) ∈ ℝ ( , 𝑡 > 0,𝒖(𝑥, 𝑦, 0) = (cid:154)𝒖 Z , 𝑥 < 0,𝒖 [ , 𝑥 > 0. (39) 它的求解和上面的一维黎曼问题(34)没有本质差异。显然,此解只反映了沿 Γ <ℓ 法向流场的变化情况。 (ii) 顶点处二维黎曼问题 这是真正的多维问题,可以写成如下形式【27】 𝜕 ! 𝒖 + 𝜕 " 𝐹(𝒖) + 𝜕 Y 𝑮(𝒖) = 0, (𝑥, 𝑦) ∈ ℝ ( , 𝑡 > 0,𝒖(𝑥, 𝑦, 0) = 𝒖 ℓ , (𝑥, 𝑦) ∈ Θ ℓ , (40) 这里 Θ ℓ , ℓ = 1, ⋯ , 𝐾 ,是从原点出发的角形区域,见图 3.1。由于有限传播性质,从中心点发出的复杂的波结构只会影响有限的区域。仅从通量的构造来说,顶点处解对界面通量的贡献占比很小,可以适当限制 Courant 数,减少顶点周边流场对数值通量的影响。在移动网格方法中,特别是拉格朗日方法中,需要用顶点解确定顶点运动速度,顶点处二维黎曼问题的解非常重要。但由于解的结构过于复杂,往往通过顶点近似黎曼解法器来处理【28】。 Θ ( Θ \ Θ / Θ ] 图 3.1. 二维黎曼问题的初始数据分布 黎曼问题的(逼近)解被广泛用到双曲守恒律及延展到更一般的流体力学方程组半离散数值方法中。与下面广义黎曼问题的解进行比较可以发现,这个解不能充分反映出流体的动力学过程;即使从微观角度,欧拉方程反映了粒子无穷碰撞的结果。再 者,在每个控制单元上及初始时间层,流场处于常状态,空间的变化依赖相邻单元之间的间断来实现。因此, Godunov 格式的解不能充分反映瞬时行为。对长时间、多物理以及多尺度现象,数值模拟结果常常出现似是而非的现象。 到目前为止黎曼问题只限于对双曲守恒律进行研究, 相关的拓展可以从下面的广义黎曼问题研究中看出。 高阶数值通量的构造实质上等价于相应广义黎曼问题的数值求解,即广义黎曼问题解法器以及高阶数值实现 ( 算法构造 ) 。即在 𝑡 = 𝑡 . 时刻,给定初始数据 𝑃 .C ∈ 𝒱 C ,求解方程组 (8) 。与黎曼问题类似, 广义黎曼问题包括界面法向广义黎曼问题以及顶点广义黎曼问题。后者只是在依赖网格顶点速度的移动网格方法中起着重要作用,有兴趣的读者可参阅【 】。下面只讨论前者。 一般地,方程组 (8) 广义黎曼问题可以写成如下形式 𝜕 ! 𝒖 + ∇ ⋅ 𝑭(𝒖, ∇𝒖, ⋯ ) = 0 , 𝒙 ∈ ℝ . , 𝑡 > 0,𝒖(𝒙, 𝒕 = 𝟎) = (cid:154)𝒖 𝑳 (𝒙), 𝚽(𝒙) < 𝟎,𝒖 𝑹 (𝒙), 𝚽(𝒙) > 𝟎, (41) 其中 𝚽(𝒙) = 𝟎 是定向的 𝑛 − 1 维光滑曲面 , 经过适当的坐标变换,可记该曲面为 𝑥 = 0 ( 记空间坐标为 𝒙 = (𝑥, 𝑦, 𝑧) ) ,即 (41) 可化为如下形式, 𝜕 ! 𝒖 + ∇ ⋅ 𝑭(cid:160)(𝒖, ∇𝒖, ⋯ ) = 𝑯(cid:160) (𝒙, 𝒖) , 𝒙 ∈ ℝ . , 𝑡 > 0,𝒖(𝒙, 𝒕 = 𝟎) = (cid:154)𝒖 𝑳 (𝒙), 𝑥 < 𝟎,𝒖 𝑹 (𝒙), 𝑥 > 𝟎, (42) 其中 𝑥 = 0 称为界面. 这里未知量仍采用同样的记号,由变换后所得形式 𝑯(cid:160) (𝒙, 𝒖) 来自于曲面
𝚽(𝒙) = 𝟎 的几何效果,下面讨论中将省略记号“~”。可以看出,广义黎曼问题与相应黎曼问题至少有三方面不同: (i) 初始数据不同, 𝒖 𝑳 (𝒙) 与 𝒖 𝑹 (𝒙) 是非常状态的两个光滑函数. 除了精度以外,数据的非一致性包含更多的物理信息,比如熵变化。(ii) 黎曼问题只对双曲守恒律而言,而广义黎曼问题的控制方程可以任意的偏微分方程,比如具有粘性力、外力等效应的流体力学方程组, 或更一般的 Boltzmann 型方程。(iii) 正如在(41)所表达那样,当
𝚽(𝒙) = 𝟎 不是一个平面时,它的拓扑变化也嵌入到解的变化中。 根据有限体积方法基本原理,需要直接逼近
𝓕(𝐴 ‘ ; 0, 𝛿) = 7 7 𝑭d𝒖(0, 𝑦, 𝑧; 𝑡)e𝑑𝑦𝑑𝑧 𝑑𝑡, ;0a (43) 这里 𝐴 ‘ = {(0, 𝑦, 𝑧); 𝑦 ∈ (𝑦 − 𝜖, 𝑦 + 𝜖), 𝑧 ∈ (𝑧 − 𝜖, 𝑧 + 𝜖)} 是面 𝑥 = 0 上的一部分。在一维情况下,(43)即为 𝓕(0; 0, 𝛿𝑡) = 7 𝑭d𝒖(0; 𝑡)e𝑑𝑡. ;!0 (44) 根据已有流体力学方程组相关的偏微分方程理论与应用研究,可以有效地逼近或求解(42),满足与(27)对应的相容性要求。具体地,(44)有
𝓕(0; 0, 𝛿𝑡) = 𝛿𝑡 ⁄𝑭d𝒖(0; 0 )e + 𝛿𝑡2 𝜕 ! 𝑭d𝒖(0; 0 )e¥ + 𝒪(𝛿𝑡 \ ). (45) 这里需要被积函数 𝑭d𝒖(0; 𝑡)e 满足关于时间 𝑡 的正则性. 对于(42)中的初始数据,这一点可以得到保障,从而可以对(44)进行相容性逼近。 定义3.3 (有限体积方法的时空耦合性). 考虑时空关联模型(8)的有限体积方法(25). 如果方程(8)的解可以用来逼近数值通量,并使(25)在定义3.1的意义下相容, 数据重构也充分利用(8)的时空关联信息,则算法(25)是时空耦合的。 注3.4. 文【22】给出的Hermite型数据重构方法是时空耦合的,该方法已用在两步四阶方法【4,29】中。文献中数据重构的时空耦合性研究还不充分,上述定义中“数据重构也充分利用(8)的时空关联信息”这句话比较模糊,需要做更深入的研究。 下面给出几个时空耦合算法的例子。 (i)线性输运方程 对于第一节中引出的线性方程(1), 其有限体积格式可以写成 𝑢z <.3/ = 𝑢z <. − 𝜆 <. [ 1Δt 7 𝑢 (cid:136)𝑥 <3/( , 𝑡(cid:137) 𝑑𝑡 − 1Δt 7 𝑢 (cid:136)𝑥 < , 𝑡(cid:137) 𝑑𝑡], 𝜆 = 𝑎Δ𝑡 . Δ𝑥 < . ! )* ! ) ! )* ! ) (46) 假设 𝑎 > 0 , 𝜆 <. < 1 。如果初始数据是分片常数 𝑢(𝑥, 𝑡 . ) = 𝑢z <. , 𝑥 ∈ 𝐼 < , (47) 则解在 (𝑡 . , 𝑡 .3/ ) 时间内有 𝑢 (cid:136)𝑥 <3 , 𝑡(cid:137) = 𝑢z <. ,及 <3/( , 𝑡(cid:137) 𝑑𝑡 = ! )* ! ) 𝑢z <. . (48) 从而由(46)可得 𝑢z <.3/ = 𝑢z <. − 𝜆 <. d𝑢z <. − 𝑢z < e. (49) 这就是迎风格式。数值通量的构造完全使用了解的信息,因此迎风格式(49)是时空耦合的。众所周知,迎风格式在所有一阶稳定格式中具有“最优”性质。 如果初始数据是分片多项式,特别是分片线性函数时 𝑢(𝑥, 𝑡 . ) = 𝑢z <. + 𝜎 <. d𝑥 − 𝑥 < e, 𝑥 ∈ 𝐼 < , (50) 则(1)的解(2)可写为 𝑢 (cid:136)𝑥 <3/( , 𝑡(cid:137) = 𝑢z <. + 𝜎 <. (cid:136)Δ𝑥2 − 𝑎(𝑡 − 𝑡 . )(cid:137) , 𝑡 ∈ (𝑡 . , 𝑡 .3/ ). (51) 直接计算可得 . <3/( , 𝑡(cid:137) 𝑑𝑡 = ! )* ! ) 𝑢 (cid:136)𝑥 <3/( , 𝑡 . + Δ𝑡2 (cid:137) =: 𝑢 <3/(.3/( , (52) 以及 ⎩⎪⎨⎪⎧ 𝑢z <.3/ = 𝑢z <. − 𝜆 <. (cid:145)𝑢 <3/(.3/( − 𝑢 < (cid:146), 𝜆 <. = 𝑎Δ𝑡 . Δ𝑥 < , 𝜎 <.3/ = 1Δ𝑥 < (cid:145)𝑢 (cid:136)𝑥 <3/( , 𝑡 .3/ (cid:137) − 𝑢 (cid:136)𝑥 < , 𝑡 .3/ (cid:137)(cid:146) . (53) 注意梯度更新 𝜎 <.3/ 也使用了解(51)的信息。这样构造的格式完全利用了解(51)的信息,所以格式(53)是时空耦合的。这里梯度 𝜎 <.3/ 更新后得到的分片线性函数仍可能有震荡现象,并不能有效逼近间断函数。从这一简单例子可以看出,数据重构仍然是值得探讨的。 更一般地,我们可以利用【29,22】的方式来构造高阶时空耦合方法,或者 Lax-Wendroff 型的单步高阶方法,只是推广到实际工程问题时,真正非线性和多维性会给单步方法带来实质性的困难。 (ii) 线性对流扩散方程 考虑下列方程 𝜕 ! 𝑢 + 𝑎𝜕 " 𝑢 = 𝜇𝜕 "( 𝑢, 𝑡 > 0, (54) 这里物理粘性系数 µ > 0 . 把它写成散度形式有 𝜕 ! 𝑢 + 𝜕 " (𝑎𝑢 − 𝜇𝜕 " 𝑢) = 0, 𝜇 > 0. (55) 然后在时空控制体上,把它表示为(32)的形式,进行通量逼近。文【8】中的守恒元/解元( Conservation Element/Solution Element )方法展示其时空耦合的属性。由于其构造需要一点篇幅展示,有兴趣的读者可以查阅【8】及其后来的一系列文章。 相对来说 , Navier-Stokes 方程组及其相关时空关联模型的时空耦合算法研究较少. 尽管形式上可以进行简单的时空对换处理,但对耗散项等高阶空间导数项的处理仍需要严格数学论证。当然,GKS 方法间接提供的
Navier-Stokes 方程的数值算法具有时空耦合性【9,30】。 (iii)线性松弛系统 考虑下列简单的松弛系统 𝜕 ! 𝑢 + 𝑎𝜕 " 𝑢 = 𝑣 − 𝑢𝜏 =: 𝑄(𝑢, 𝑣)𝜏 , 𝑡 > 0,𝑢(𝑥, 0) = 𝑢 (𝑥), (56) 其中 τ > 0 是松弛参数, 𝑎 为常数, 𝑣 也是常数。这个方程的有限体积形式为: 𝑢z <.3/ = 𝑢z <. − 𝜆 (cid:145)𝑢 <3/(.3/( − 𝑢 < (cid:146) + 1Δ𝑥 7 7 𝑄(𝑣, 𝑢)𝜏 𝑑𝑥𝑑𝑡 " (* " (8 ! )* ! ) . (57) 为了构造时空耦合算法,可以用(56)的解表达式, 𝑢(𝑥, 𝑡) = (𝑢 (𝑥 − 𝑎𝑡) − 𝑣)𝑒 + 𝑣, (58) 来计算数值积分 𝑢 <3 .3 以及 c(d,e)b 的积分。为了更一般的应用,这里使用【30】中的方法,即 DUGKS 方法,得到 𝑢z <.3/ = 𝑢z <. − 𝜆 (cid:145)𝑢 <3/(.3/( − 𝑢 < (cid:146) + Δ𝑡2τ d𝑄 <. + 𝑄 <.3/ e, (59) 其中 𝑢 <3 .3 是由下列形式给出 𝑢 <3/(.3/( = 𝑢 (cid:136)𝑥 <3/( − 𝑎2 Δ𝑡, 𝑡 . (cid:137) + 12𝜏 (cid:145)𝑄 (cid:136)𝑥 <3/( − 𝑎2 Δ𝑡, 𝑡 . (cid:137) + 𝑄 <3/(.3/( (cid:146) . (60) 进而可以利用(58)或使用(60)类似的想法,得到 𝑢 (cid:136)𝑥 <3 , 𝑡 .3/ (cid:137) ,并用来更新梯度 𝜎 <.3/ = 1Δ𝑥 < (cid:145)𝑢 (cid:136)𝑥 <3/( , 𝑡 .3/ (cid:137) − 𝑢 (cid:136)𝑥 < , 𝑡 .3/ (cid:137)(cid:146). (61) 很显然,格式(59)-(61)是时空耦合的,并具有时空二阶精度。 四 基于广义黎曼解法器 (GRP solver)的时空耦合算法 数值通量的构造是执行有限体积格式的两个核心步骤之一。对于非线性问题我们一般无法像上述线性方程那样,得到解的表达式。下面将要利用广义黎曼问题(41)或(42)解的局部正则性质,发展广义黎曼解法器。其主要思想可以概括为: Lax-Wendroff 型解法器 + 奇异性分辨(间断追踪) 为了方便叙述,这节统一把界面放在 𝑥 = 0 , 把广义黎曼问题的求解点平移到坐标原点 (0,0) , 初始数据表示成(42)的形式。记 𝒖 ∗ : = lim !→0 * 𝒖(𝟎, 𝑡), (𝜕 ! 𝒖) ∗ : = lim !→0 * 𝜕 ! 𝒖(𝟎, 𝑡). (62) 广义黎曼问题(GRP)解法器: 给定控制方程以及形如(42)的分片光滑函数作为初始数据,求解 (42) 并得到形如(61)的极限值。 一旦有了极限值(62),受益于解 𝒖(𝒙, 𝑡) 在界面上关于时间 𝑡 的连续性,我们有 𝒖(0, 𝛿𝑡) = 𝒖 ∗ + (𝜕 ! 𝒖) ∗ 𝛿𝑡 + 𝒪(𝛿𝑡 ( ). (63) 由于 𝒖 ∗ 只表示了一种瞬时行为,其解由相应的黎曼解给定 𝒖 ∗ = 𝑅d0; 𝒖 𝑳 (0), 𝒖 𝑹 (0)e. (64) 接下来的关键任务是求值 (𝜕 ! 𝒖) ∗ 。“非常不严格”地由控制方程(42)直接得到 𝜕 ! 𝒖 = −∇ ⋅ 𝑭(cid:160)(𝒖, ∇𝒖, ⋯ ) + 𝑯(cid:160) (𝒙, 𝒖) , (65) 然后取极限。这样解 𝑢(𝑥, 𝑡) 的变化率可以通过空间的变化反映,就是经典的 Lax-Wendroff 方法【23】,或
Cauchy-Kowalevski 方法【31】。从中可以看出,通过这种途径把模型的时空关联性质嵌入到数值格式中,从而所得的算法具有时空耦合性质。 注 : 一般的单步 Lax-Wendroff 方法需要计算解的高阶时间导数 ∂ ,h 𝐮 ,但对流体力学方程组来说,需要避免这样的操作,原因是:(i)由于解中常有间断,这样的运算失去了数学与物理意义;(ii)即使允许如此运算,所得出的通量估算过于复杂,失去了实际工程意义;(iii)最近发展的两步四阶方法【4, 29】表明,这一对值具有明确物理意义,可以用来作为构造高阶方法的基石,并避免高阶时间导数的运算。 对于线性系统 𝝏 𝒕 𝒖 + 𝑨𝝏 𝒙 𝒖 = 𝑪𝒖 + 𝑫, 𝒕 > 𝟎,𝒖(𝒙, 𝟎) = (cid:154)𝒖 𝑳 (𝒙), 𝒙 < 𝟎,𝒖 𝑹 (𝒙), 𝒙 > 𝟎, (66) 其中 𝑨, , 𝑪, 𝑫 是常数矩阵 , 𝑨 可以对角化 , 𝚲 = 𝑹 𝐀𝑹 , 其中 𝚲 是由 𝑨 的特征值 𝜆 I 组成的矩阵, 𝚲 = 𝑑𝑖𝑎𝑔{𝜆 𝒊 } 。记 𝒘 = 𝑹 𝒖 , 有 𝝏 𝒕 𝒘 + 𝚲𝝏 𝒙 𝒘 = 𝑹 𝑪𝒖 + 𝑹
𝑫, 𝒕 > 𝟎,𝒘(𝒙, 𝟎) = (cid:154)𝒘 𝑳 (𝒙), 𝒙 < 𝟎,𝒘 𝑹 (𝒙), 𝒙 > 𝟎, (67) 进 而 记 𝚲 = 𝑑𝑖𝑎𝑔(max(𝜆 I , 0)} , 𝚲 = 𝑑𝑖𝑎𝑔(min(𝜆 I , 0)} , 𝑰 = /( 𝑑𝑖𝑎𝑔{1 + 𝑠𝑖𝑔𝑛(𝜆 I )} , 𝑰 = − /( 𝑑𝑖𝑎𝑔{1 − 𝑠𝑖𝑔𝑛(𝜆 I )} . 显然我们可到 (𝝏 𝒕 𝒘) ∗ = −(𝚲 𝒘 𝑳j + 𝚲 𝒘 𝑹j ) + d𝑰 𝑹 𝑪𝒖 𝑳 + 𝑰 𝑹 𝑪𝒖 𝑹 e + 𝑹 𝑫, (68) 其中 𝒖 𝑳 = 𝒖 𝑳 (𝟎) , 𝒖 𝑹 = 𝒖 𝑹 (𝟎) , 𝒘 𝑳j = 𝒘 𝑳 ′(𝟎) , 𝒘 𝑹j = 𝒘 𝑹 ′(𝟎) . 再记 𝑨 ± = 𝑹𝚲 ± 𝑹 , 回到原始变量 𝒖 , 我们有 (𝝏 𝒕 𝒖) ∗ = −(𝑨 𝒖 𝑳j + 𝑨 𝒖 𝑹j ) + d𝑹𝑰 𝑹 𝑪𝒖 𝑳 + 𝑹𝑰 𝑹 𝑪𝒖 𝑹 e + 𝑫. (69) 特别当 𝒖 𝑳 = 𝒖 𝑹 = 𝒖 ∗ 时有 (𝝏 𝒕 𝒖) ∗ = −(𝑨 𝒖 𝑳j + 𝑨 𝒖 𝑹j ) + 𝑪𝒖 ∗ + 𝑫. (70) 从中可以看出如何应用 Lax-Wendroff 解法器和间断追踪计算 (𝝏 𝒕 𝒖) ∗ 。 多维情形是类似的.例如考虑二维线性方程组 𝝏 𝒕 𝒖 + 𝑨𝝏 𝒙 𝒖 + 𝑩𝝏 𝒚 𝒖 = 𝑪𝒖 + 𝑫, 𝒕 > 𝟎,𝒖(𝑥, 𝑦, 𝟎) = (cid:154)𝒖 𝑳 (𝑥, 𝑦), 𝑥 < 𝟎,𝒖 𝑹 (𝑥, 𝑦), 𝑥 > 𝟎, (71) 这里 𝑨, 𝑩, 𝑪, 𝑫 是常矩阵 , 且 𝑨, 𝑩 可对角化. 切向变化量 𝑩𝝏 𝒚 𝒖 可看作相对法向的一个扰动, 将线性方程组写成如下形式 𝝏 𝒕 𝒖 + 𝑨𝝏 𝒙 𝒖 = −𝑩𝝏 𝒚 𝒖 + 𝑪𝒖 + 𝑫, 𝑡 > 0,𝒖(𝑥, 𝑦, 0) = (cid:154)𝒖 Z (𝑥, 𝑦), 𝑥 < 0,𝒖 [ (𝑥, 𝑦), 𝑥 > 0. (72) 与上面类似方法,我们可以得到 (𝝏 𝒕 𝒖) ∗ = −[𝑨 (𝝏 𝒙 𝒖) 𝑳 + 𝑨 (𝝏 𝒙 𝒖) 𝑹 ] − [𝑹𝑰 𝑹 𝑩(𝝏 𝒚 𝒖) 𝑳 + 𝑹𝑰 𝑹 𝑩(𝝏 𝒚 𝒖) 𝑹 ] +(𝑹𝑰 𝑹 𝑪𝒖 𝑳 + 𝑹𝑰 𝑹 𝑪𝒖 𝑹 ) + 𝑫, (73) 这里记号与上述意思相同,并有 (𝜕 " 𝒖) Z = lim "→0 𝜕 " 𝒖 (𝑥, 𝑦) , (𝜕 " 𝒖) [ = lim "→03 𝜕 " 𝒖 (𝑥, 𝑦) 等。 从(70)和(73)可以看出,所有的信息都是迎风且时空耦合的。 特别强调,通过 GRP 解法器把切向变化自然地蕴含到数值通量中,这是法向(一维)黎曼解法器做不到的【38】。换句话说,如果只使用黎曼问题解法器,导致格式失去多维性,该格式应用到多维问题模拟时自然会有数值缺陷,不得不用其它方式进行补偿。 对于非线性系统来说,当流场是光滑的或者波的强度不大时,使用近似解法器进行通量的逼近就可以取得不错的效果,即声波近似或线性化 GRP 解法器。Toro 等的ADER 解法器就是 GRP 的声波近似版本【35】。 先看下面的一维双曲平衡律( hyperbolic balance laws ) 𝝏 𝒕 𝒖 + 𝜕 " 𝑭(𝒖) = 𝑯(𝑥, 𝒖), 𝑥 ∈ ℝ, 𝑡 > 0,𝒖(𝑥, 0) = (cid:154)𝒖 Z (𝑥), 𝑥 < 0,𝒖 [ (𝑥), 𝑥 > 0. (74) 可以把它看作描写相对于网格界面法向效应的方程。如果 𝒖 Z (0 − 0) = 𝒖 [ (0 + 0) , 而 𝒖 Zj (0 − 0) ≠ 𝒖 [j (0 + 0) , 只有线性波(声波)从 (0,0) 点发出.记 𝑨(𝒖) = 𝜕 𝒖 𝑭(𝒖) , 𝒖 ∗ =𝒖 Z (0 − 0) = 𝒖 [ (0 + 0) , 𝒖 = 𝒖 ∗ + 𝒗, 𝒗 可看成 𝒖 在 𝒖 ∗ 附近的扰动。然后在 𝒖 ∗ 处线性化(73)可得 𝝏 𝒕 𝒗 + 𝑨(𝒖 ∗ )𝝏 𝒙 𝒗 = 𝑯(𝑥, 𝒖 ∗ ), 𝒖 = 𝒖 ∗ + 𝒗, (75) 这是一个线性系统。类似于(66),从中可以计算 (𝜕 ! 𝒖) ∗ = (𝜕 ! 𝒗) ∗ , (𝝏 𝒕 𝒗) ∗ = −𝑨 (𝒖 ∗ )𝒗 𝑳j (0) − 𝑨 (𝒖 ∗ )𝒗 𝑹j (0) + 𝑯(0, 𝒖 ∗ ), (76) 这 里 𝑨 ± (𝒖 ∗ ) = 𝑹(𝒖 ∗ )𝚲 ± (𝒖 ∗ )𝑹 (𝒖 ∗ ) , 𝚲 ± (𝒖 ∗ ) 是 由 𝚲(𝒖 ∗ ) 的 “ ± ” 部 分 组 成 𝜦 =𝑑𝑖𝑎𝑔{max(𝜆 I , 0)} , 𝜦 = 𝑑𝑖𝑎𝑔{min(𝜆 I , 0)} 。 从而 (𝜕 ! 𝒖) ∗ = −𝑨 (𝒖 ∗ )𝒖 Zj (0) − 𝑨 (𝒖 ∗ )𝒖 [j (0) + 𝑯(𝑥, 𝒖 ∗ ). (77) 这就是线性化 GRP 解法器 . 当 𝒖 Z (0 − 0) ≠ 𝒖 [ (0 + 0) 时,声波近似策略如下:以局部黎曼解 𝒖 ∗ =𝑅(0; 𝒖 Z (0), 𝒖 𝑹 (0)) 为背景解, 线性化(74)得到(75),从而可得线性化 GRP 解法器(77),具体见【32】。 对于多维系统(仍然以二维为例), 𝝏 𝒕 𝒖 + 𝜕 " 𝑭(𝒖) + 𝜕 Y 𝑮(𝒖) = 𝑯(𝑥, 𝑦, 𝒖), 𝑡 > 0,𝒖(𝑥, 𝑦, 0) = (cid:154)𝒖 Z (𝑥, 𝑦), 𝑥 < 0,𝒖 [ (𝑥, 𝑦), 𝑥 > 0, (78) 采用声波近似的策略,线性化这个方程组得到 𝝏 𝒕 𝒗 + 𝑨(𝒖 ∗ )𝜕 " 𝒗 = −𝑩(𝒖 ∗ )𝜕 Y 𝒗 + 𝑯(𝑥, 𝑦, 𝒖 ∗ ), 𝑡 > 0,𝒗(𝑥, 𝑦, 0) = (cid:154)𝒗 Z (𝑥, 𝑦), 𝑥 < 0,𝒗 [ (𝑥, 𝑦), 𝑥 > 0, (79) 当源项 𝑯(𝑥, 𝒖) 存在时,即使 𝒖 𝑳: = 𝒖 𝑹: ≡ 𝟎, (𝜕 𝑢) ∗ ≠ 0 ,意味着 GRP 解法器仍然起着重要作用。 从而得到类于(73)的解法器, 这里 𝑩(𝒖 ∗ ) = 𝜕 𝒖 𝑮(𝒖) 。 在(74)中,如果初始数据在原点(相对网格尺寸)有“很大”跳跃,就会出现强非线性波,线性化 GRP 解法器不足以分辨强非线性波,这时候需要使用真正非线性 GRP解法器,否则就会出现明显的数值缺陷【7,34,35】. 一般的界定标准为 ∥ 𝒖 𝑳 (𝟎) − 𝒖 𝑹 (𝟎) ∥≫ 𝜶𝚫𝒙 , (80) 参数 𝛼 > 0 是一个重要的经验参数。求解 GRP 解法器的基本思想是: 解析强稀疏波,追踪强间断。 特别需要指出,在强间断的情形下,热力学效应起着重要作用,真正非线性的 GRP 解法器反映了 热力学相容的性质 【7】。 本文不具体给出真正非线性
GRP 解法器.对于可压缩欧拉方程组,可参阅【5,36】,后期发展的不依赖于坐标系统的
GRP 解法器,详见【6】;对一般的双曲方程组,可参见【37】。 对于多维情形,仍然可以把切向的变化和间断面的变化看作扰动,考虑有下列形式问题 𝝏 𝒕 𝒖 + 𝜕 " 𝑭(𝒖) = −𝜕 Y 𝑮(𝒖) + 𝑯(𝑥, 𝑦, 𝒖), 𝑡 > 0,𝒖(𝑥, 𝑦, 0) = (cid:154)𝒖 Z (𝑥, 𝑦), 𝑥 < 0,𝒖 [ (𝑥, 𝑦), 𝑥 > 0. (81) 从而利用 4.1 和 4.2 节中法向解法器求解(72)和(73),这主要源自于一个关键的事实: 切向扰动在法向的传播是连续的。 由此切向的变化在通量中得到充分反映,得到的数值方法具有真正的多维性【38】。 正如前面所讨论的那样,除非涉及(依赖于网格移动)中心拉格朗日型数值方法,这 里 不 关 注 顶 点 多 维 黎 曼 解 法 器。 由 于 真 正 多 维 黎 曼 问 题 的 理 论 基 础 很 不 成 熟【27,39,40,41】,真正多维顶点黎曼解法器常常带来不稳定的数值结果【42】。 在实践中,人们更倾向使用健壮的逼近 GRP 解法器, 例如,Maire 等利用新的框架并结合拉格朗日型 GRP 解法器【43】,发展了稳定的中心型高阶拉氏方法。 对于一般的时空关联模型,如多相流、多介质模型【44】和 Navier-Stokes 方程等,广义黎曼问题解法器可以进行类似的构造,现进行简述. 对于多相流、多介质模型,广义黎曼问题解法器可以归结为一般的非线性平衡律的范畴,模拟的时空关联性质在 GRP 解法器中可以得到充分体现【45,46,47】。困难在于由介质组份以及混合引起的数值振荡、物质分数为负等非物理解,涉及物理建模本身的缺陷以及有限体积格式的投影(平均)过程。由于热力学关系的高度非线性,投 影过程不能反映精确的物理过程, 例如内能的平均过程导致物质界面附近的压力振荡. 因此,在投影步实现时空耦合,充分利用解的信息也许可以提高相关数值质量。 对于
Navier-Stokes 方程组等宏观模型,基本的想法是类似的。对于对流项(欧拉部分),使用上述标准的广义黎曼解法器。需要指出的是: 初始数据(41)需要进行适当的匹配。 也就是说该初始数据至少是二阶以上分片多项式。另外可使用的策略为: 对于对流占优问题,采用松弛策略,把相应模型双曲化 【48,49】。这样的策略是基于能量原理: 在对流占优情形下,能量集中于初始能量传播影响区域内(双曲性质)。 从 Boltzmann 方程直接出发的动理学方法 (Kinetic methods) ,是设计流体力学数值方法的一条有效途径. 一般来说,很难写出
Boltzmann 方程的解,但在特定的近似假定下,如 BGK 模型【50】,可以类似于(58)那样隐式得出解的表达式,并将之用于数值通量的构造,得出的数值方法如气体动理学格式 (GKS) 【9】、统一气体动理学格式
UGKS 【10】等。特别地,
GKS 可以用来模拟
Navier-Stokes(N-S) 方程描述的流动,可以作为
Navier-Stokes 的解法器来使用。按照这样的定义,在通量的构造部分,
GKS 和 UGKS 解法器显然是时空耦合的。 在高精度数值方法中,黎曼解法器和 GRP 解法器都可用来构造数值通量。前者在高阶线方法( high order methods of line )中常用,如
Runge-Kutta 型高阶方法,后者用在两步四阶方法中。下面仍以一维双曲守恒律方程组(31)来比较两种解法器的差别。考察控制体 ⁄𝑥 < , 𝑥 <3 ¥ × [𝑡 . , 𝑡 .3/ ) 边界上通量。在 𝑥 = 𝑥 <3 , 记 𝒖 <3/(. = lim !→! ) 𝒖 (cid:136)𝑥 <3/( , 𝑡(cid:137) , (𝜕 ! 𝒖) <3/(. = lim !→! ) 𝜕 ! 𝒖 (cid:136)𝑥 <3/( , 𝑡(cid:137). (82) 由解关于时间 𝑥 的正则性得到 𝑭 (cid:145)𝒖 (cid:136)𝑥 <3/( , 𝑡(cid:137)(cid:146) = 𝑭 (cid:145)𝒖 <3/(. (cid:146) + 𝑭′ (cid:145)𝒖 <3/(. (cid:146) (𝜕 ! 𝒖) <3/(. (𝑡 − 𝑡 . ) + 𝒪((𝑡 − 𝑡 . ) ( ). (83) (i) 黎曼解法器. 选取 𝑭 <3 n7oe = 𝑭 (cid:136)𝒖 <3 . (cid:137) , 从而有 <3/( , 𝑡(cid:137)(cid:146) − 𝑭 (cid:145)𝒖 (cid:136)𝑥 < , 𝑡(cid:137)(cid:146)(cid:131) 𝒕 𝒏*𝟏 𝒕 𝒏 𝑑𝑡 − Δ𝑡 . [𝑭 <3/(n7oe − 𝐹 < ]= (Δ𝑡 . ) j (cid:145)𝒖 <3/(. (cid:146) (𝜕 ! 𝒖) <3/(. − 𝑭 j (cid:145)𝒖 < (cid:146) (𝜕 ! 𝒖) < (cid:146) + 𝒪((Δ𝑡 . ) \ ) (84) 所以对于间断解来说,通量的截断误差为 𝐸 <. = 𝒪((Δ𝑡 . ) ( ) (85) 误差阶 𝑞 = 0 。不过,对于解的光滑区域来说,(84)中的差赋予了一个“阶数”, 𝐸 <. = 𝒪((Δ𝑡 . ) \ ) (86) 从而误差阶为 𝑞 = 1 。 (ii) GRP 解法器. 选取 𝑭 <3 n[p = 𝑭 (cid:136)𝒖 <3 . (cid:137) + (%! ) ) " 𝑭′ (cid:136)𝒖 <3 . (cid:137) (𝜕 ! 𝒖) <3 . , 从而有 <3/( , 𝑡(cid:137)(cid:146) − 𝑭 (cid:145)𝒖 (cid:136)𝑥 < , 𝑡(cid:137)(cid:146)(cid:131) 𝒕 𝒏*𝟏 𝒕 𝒏 𝑑𝑡 − Δ𝑡 . [𝑭 <3/(n[p − 𝐹 < ]= (Δ𝑡 . ) !( 𝑭 (cid:145)𝒖 <3/(. (cid:146) − 𝜕 !( 𝑭 (cid:145)𝒖 < (cid:146)˛ + 𝒪((Δ𝑡 . ) r ). (87) 与上面讨论类似,对于间断解, GRP 解法器有一阶精度 𝑞 = 1 ,而对于光滑解有二阶精度 𝑞 = 2 。 边界条件的数值处理是计算流体力学中的一个核心问题,大部分的数值处理基于对虚拟区域的延拓或特征传播理论【51】。最近逆 Lax-Wendroff 方法【52】用来数值处理边界条件,这种方法蕴含了时空耦合性。我们认识到,流体力学方程组的数值边界条件等价于单边广义黎曼问题问题的数值求解( one-sided GRP solver )【53】。事实上,先前的工作已经认识到单边黎曼问题在数值边界条件处理上的重要性【54,55】。这些工作与直接的逆
Lax-Wendroff 方法相比,有以下特点: (i) 在有限体积框架下,计算区域的边界自然为控制体的边界,所以数值边界条件处理等价于边界上的通量数值逼近,它归结为单边黎曼解法器的构造。 (ii) 单边黎曼解法器中的 (𝜕 ! 𝑢) ∗ 其实反映了边界条件上的受力情况,即牛顿第二定律的数学表现,这是时空耦合算法的体现。在【53】中,这一事实是单边黎曼解法器的基石。 (iii) 在技术上,如此处理数值边界条件可以尽量少使用虚拟网格。在时空二阶格式中无需使用虚拟网格; 相比较其它更高阶方法,至多使用一半的虚拟网格。 (iv) 无需使用更高阶导数,避免了人为的复杂性和虚假物理性质等【52】。 单边黎曼解法器的使用将极大简化了边界条件的数值处理,特别是解决了无结构网格下虚拟区域的插值问题【56】。 高阶精度数值方法对小尺度问题的数值模拟非常重要,如湍流等。基于有限体积格式的框架(25),高阶方法需要在数值通量和数据重映两部分具有高阶逼近性质.在通 量的逼近部分,可以采用 Lax-Wendroff 方法,但流体力学方程组(8)的非线性性质导致高阶展开过于复杂,缺乏实际工程价值。更糟糕的是,解的间断性质意味着高阶展开的运算没有意义, 因此需要寻求其它方式。在过去的几十年中,各种空间重构技术以及
Runge-Kutta 型时间推进方法取得了巨大的进步,读者可以容易查阅,这里不再赘述。下面简述一类时空耦合方法。 文【29】发展了两步四阶时间推进方法,成功地融合了
Runge-Kutta 型方法简单性优点和基于 GRP 解法器的时空耦合性质。文【4】详述了该类方法的特性: (i)计算效率:同等条件下其计算效率是 Runge-Kutta 方法的一半(二维情形为例)【57】;(ii) 紧致性:由于时间推进步骤减少一半,模版就会减少一半;时空耦合的 HWENO 型重构【22, 58】 可以使计算格式更加紧致; (iii) 稳定性:已经证明其稳定区域比 Runge-Kutta 大【59】;(iv) 兼容性:该方法很容易和其它方法相结合,如 GKS 解法器【60】、多矩方法【61】以及 CRP 重构技术【62】等。 目前,该类方法还局限在“显式”框架内,相关的研究处于起步阶段。为了应用的需要,亟须发展“隐式”高阶方法,以适应“强刚性”等多尺度问题的求解。 五 数值算例 在可压缩流体的算法设计及其数值模拟中,有大量基准算例(
Benchmarks )。随着算法进步和工程需求的提高,基准算例应该与时俱进,面对相当极端的环境. 在【63】中,一系列基准算例值得新的数值算法测试。下面几个算例,可以看出时空耦合性的重要性。 算例 5.1. 大压力比(大密度比)问题 这个问题又称为 LeBlanc 问题,它是经典的激波管问题的极端情形,从中可以看出数值方法对强稀疏波的分辨以及对激波的影响。控制方程为一维欧拉方程(31),初始数据具有如下形式 (𝜌, 𝑢, 𝑝)(𝑥, 0) = (cid:154)(10 r , 0,10 r ), 0 ≤ 𝑥 < 0.3,(1,0,1), 3 ≤ 𝑥 ≤ 1.0. (88) 多方指数取 ,计算的终止时间 𝑡 = 0.12 。文【64】利用这个例子检测了多个数值格式的性能。这里我们同样用这个算例针对性讨论本文涉及的一些观点,数值结果来源于【7,29】。 首先在图 5.1(A)中展现使用不同解法器的一阶格式。可以看出一阶 Godunov 方法随着网格加密,可以收敛到精确解;使用逼近解法器,收敛相对较慢,但仍然具有收敛性. 图 5.1(B)使用了空间二阶重构,而时间离散使用一阶向前欧拉方法,解法器分别使用黎曼解法器和逼近的黎曼解法器(HLLC,Roe), 数值表现很差,不具有收敛性。正如在 4.5 节所述,当空间具有高阶精度,使用一阶黎曼解法器构造数值通量等,算法不具有相容性;类似的情形反映在图 5.1(C), 即使时间离散使用二阶 Runge-Kutta 方法。 图 5.1(A) 一阶格式的数值结果. 图标 Godunov, HLLC, Roe 分别表示不同的解法器 图 5.1(B) 空间二阶时间一阶的数值结果. 间断处格式不具有收敛性 图 5.1(D)考察了在这种极端情况下使用线性化方法的数值表现,黎曼解使用声波近似解法器。其实【34,35】已经指出线性化解法器的数值缺陷。图 5.1(E)使用了非线性 GRP 解法器,数值结果立即改善,说明了强非线性出现后真正非线性 GRP 解法器的重要性。 图 5.1(F)展示的结果是基于 GRP 解法器的两步四阶方法【29】,从中看到用 100网格点可以很好分辨间断,300 网格点可以得到完美效果。这是许多数值方法很难达到的,从而说明了时空耦合算法的必要性。 图 5.1(C) Runge-Kutta 型二阶方法,使用精确和逼近黎曼解法器构造数值通量 图 5.1(D) 单步声波近似的 GRP 方法 图 5.1(E) 二阶 GRP 格式,使用了真正非线性 GRP 解法器 算例 5.2. 激波和熵波相互作用的问题 这个算例是 Shu-Osher 算例的推广,考虑了更极端的情形,又称为
Titarev-Toro 问题【62】。控制方程依然为欧拉方程,初始数据为 (𝜌, 𝑢, 𝑝)(𝑥, 0) = (cid:154)(1.515696,0.52336,1.805) , − 5 ≤ 𝑥 < −4.5,(1 + sin (20𝜋𝑥),0,1) , − 4.5 ≤ 𝑥 ≤ 5. (89) 这里使用了 GKS 解法器【9】得到图 5.2 的数值结果,详细说明见【60】。 图 5.1(F) 高阶方法. 这里 m 表示所用的网格点数,RK4-WENO5 表示使用空间 5阶 WENO 重构,4 阶 Runge-Kutta 时间推进方法;GRP4-HWENO5 表示使用空间 5阶 HWENO 重构,两步四阶时间推进方法 图 5.2.
Titarev-Toro 问题 .这里使用 1000 网格点数,空间重构采用了不同的重构策略,计算终止时间为 𝑡 = 5 算例 5.3. 激波和悬浮刚体的相互作用 这个数值结果取自【54】,模拟了隧道 [0,1] × [0,0.2] 内激波撞击一个长 𝑎 = 0.06 、宽 𝑏 = 0.03 矩形刚体的流场情况。控制方程为两维欧拉方程方程组,我们使用了时空耦合二阶 GRP 方法,并用单边广义黎曼解法器刻画方块移动时的数值边界条件。初始时刻一个马赫 3 的激波位于 𝑥 = 0.08 处,波前状态为 (𝜌 , 𝑢 , 𝑣 , 𝑝 ) = (1.3,0,0,0.1) ;该 方块的密度为 𝜌 ) = 10𝜌 ,质量为 𝑀 ) = 𝑎𝑏𝜌 ) 。初始时刻该方块以倾斜度 sr 被放置于隧道中,重心为 (0.15 , ,惯性力为 𝐴 ) = t = \ (𝑎 ( + 𝑏 ( ) 。图 5.3 展示了两个时刻的流场压力分布情况。 图 5.3. 激波和悬浮刚体的相互作用后压力的分布情况。这里使用
800 × 160 个网格,上图计算终止时刻为 𝑡 = 0.6 ,下图为 𝑡 = 1.0 算例 5.4. 多介质激波和起泡的相互作用 这个例子展示了激波和气泡相互作用的情况,该问题已经成为多相流领域一个标准算例【66】。下面的结果取自【45】,分别用能量分裂的 Godunov 方法和 GRP 方法模拟所得结果。显然,GRP 很好提高了流场的分辨率。 图 5.4(A).激波撞击氦气泡的装置示意图。初始时刻激波的马赫数是 𝑀 u = 1.22 ,计算网格为 , 其它参数详见【45】的算例 D 图 5.4(B).激波撞击气泡不同时刻的阴影图,每个子图的左边是 Godunov 格式的结果,右边是 GRP 的结果. (a) 𝑡 = 32𝜇𝑠 ; (b) 𝑡 = 62𝜇𝑠; (c) 𝑡 = 72𝜇𝑠; (d) 𝑡 = 102𝜇𝑠; (e) 𝑡 = 427𝜇𝑠; (f) 𝑡 = 674𝜇𝑠. 算例 5.5. 各向同性可压缩湍流的模拟 下面的算例采用了两步四阶 GKS 方法【60】模拟的超音速各向同性可压缩湍流的情况【67】,从中可查阅具体的计算参数,其中计算区域是 [−𝜋, 𝜋] \ 。该文展示了两步四阶方法和二阶方法同样健壮。对许多高阶方法来说,分辨激波碎片(shocklet)是个巨大的挑战【68】,该方法中 GKS 解法器的时刻耦合性质继承了物理本身的属性。 图 5.5. 各向同性可压缩湍流的模拟【68】。其中湍流马赫数 𝑀𝑎 ! = 1.2 , 初始泰勒微尺度雷诺数 𝑅𝑒 v = 72 , 速度梯度张量第二不变量的等值面 𝑄 = 六 讨论与展望 流体力学的时空关联模型刻画了物理量空间变化在时间方向上的演化的传播,如各种波的传播等。当进行数值模拟时,这种性质理应得到继承,相应地所用算法应该具备时空耦合特性. 实践过程中这个看似简单课题理应得到关注。遗憾的是,相对于时空解耦方法,时空耦合算法需要更好理解模型的物理性质或数学理论,人们自然“避难就易”。一个直白的问题是: “即使物理问题的数学模型十分完美,当相应的算法缺乏相应完美性质时,其数值结果怎么令人信服呢 ?” 本文首先严格论述有限体积方法. 正如引言所解释的原因:(1)无论从物理定律的形成还是数值算法的构造,有限体积方法是解流体力学问题最自然的框架,本文总结了这个框架形成的数学理论和内在涵义;(2)流体现象的复杂性(如间断等)客观地阻碍了有限体积方法严格数学理论的形成,对于这一框架的认识常出现诸多似是而非的争论; (3)对于许多实际工程问题,基于无结构网格的算法设计是一个自然要求,在此之上许多方法相互冲突【24】,有必要从最基本的地方出发建立一个根本准则。幸运地是,许多工程软件,如 Fluent , 自动地在有限体积框架下生成;许多工程人员当遇到难以跨越的困难时,往往借用有限体积框架度过难关。从论证过程可以看到时空耦合是算法的根本属性。 有限体积格式的步骤很简单: 数据重构(投影)和数值通量的构造 。目前 CFD 算法的大部分研究者将注意力集中于前者,基本上在空间逼近论的范畴内进行探索。尽管黎曼不变量(
Riemann invariants )等物理量以及其它技术用于重构过程,但是重要的时空耦合性质相当缺乏 。本文侧重于后者的讨论,给出了有限体积格式与积分平衡律(12)之间的相容性准则。特别对于数值通量介绍了黎曼解法器和 GRP 解法器,并在4.5 节中进行了对比。简单地总结如下,具体内容可以到【2】中查阅。 (i) 关于 Riemann 解法器 . 对于双曲守恒律(欧拉方程),当初始数据是 分片 常 数 时 。 由 Riemann 解 法 器 产 生 的 数 值 通 量 是 无 穷 阶 相 容 , 即Godunov 格式就是积分平衡律;当初始数据是 非常数的分片光滑函数 时,Riemann 解法器产生的数值通量对光滑解是一阶相容的,但对间断解是不相容的, 见 4.5 节。也就是说对于高阶空间重构,如果不能有匹配的数值通量,产生的数值格式是不相容的。值得注意的是,如果控制方程包含外力项时,Godunov 格式永远是不相容的。 (ii) 关于 GRP 解法器 . GRP 解法器给出时空耦合的数值通量.对于光滑流场,GRP 解法器得到二阶相容的数值通量;但对于间断解只有一阶时间精度。GRP 解法器是保证有限体积格式的逼近解收敛的一个必要条件。 这里之所以再次强调这点并细致给出算例 5.1, 原因在于为了看清它与传统理解的差异。事实上,算法时空耦合性也体现在其它方面,比如最近研究气体动理学的渐近性质时,提出了统一保持性质 ( unified preserving property, UP ) 的概念【69】。对于一个动理学方法,不仅应该具有欧拉层面的渐近保持性质( Asymptotic preserving property, 在随机选取( Random Choice )方法中, 黎曼问题问题的解被用在随机选取中. 在某种意义下,它可以看成是一种时空耦合投影方法。 AP ),还应该根据需要具有 Navier-Stokes 或更深层面的 UP 性质,这个过程实际上是通过时空耦合的思想实现的。该文分别用时空耦合 DUGKS 方法以及时空解耦 CLR 格式,对两维不可压缩
Taylor 涡的进行数值模拟,见图 6.1,具体见【69】。 图 6.1. 关于二维不可压缩 Taylor 涡的 DUGKS(时空耦合)和 CLR(时空解耦)模拟的比较 由于篇幅限制,本文只给出了有限体积方法基本原理的相容性准则以及通过 GRP解法器实现相容性的技术细节,对很多根本性质还没有涉及,如时空相容格式的熵稳定性等。对于可压缩流体力学来说,熵稳定性对应于热力学相容性【7】,这部分留在将来的文章中探讨。 客观地说,时空耦合算法的研究非常不成熟,目前只在相对比较“纯粹”的模型上进行分析和验证,但是这一思想应该加以推广与应用,比如如何将时空耦合算法的思想应用于一般的时空关联湍流模型【70】,粒子建模和时空耦合算法等。就算法本身来说,时空耦合的隐式格式研究还没有开展,这一领域的研究应该是下一阶段的重点。 在工程应用中时空耦合的程序显得更少,一方面是受限算法模块的历史延续性制约;另一方面是工程领域内的“理性”力学越来越少。加强工程算法的科学性是需要目前亟待解决的观念问题。 后记. 作为一个数学工作者 ,很忐忑地接受邀请在力学的专业刊物《空气动力学学报》奉献此类稿件,毕竟隔行如隔山;但想到数学、力学本就是“一家”,向力学家“学习”本就是“数学”的一大研究动机,心里稍微坦然点。 致谢: 此文是和我的下列合作者合作过程中完成的,在此一并致谢:他们是
Matania Ben-Artzi, 齐进、汪 玥 、杜知方、雷昕、徐昆、郭照立、李启兵等。特别感谢曹贵瑜博士提供的关于各向同性可压缩湍流模拟算例 5.5, 感谢徐昆教授、郭照立教授和汪 玥 副研究员提出许多根本性改进意见,感谢陈海波给予了文字上的润色。 A spacetime outlook on CFD: Spacetime correlated models and spacetime coupled algorithms
Jiequan Li Institute of Applied Physics and Computational Mathematics, Beijing; Center for Applied Physics and Technology, Peking University
Abstract:
A spacetime outlook on Computational Fluid Dynamics is advocated: models in fluid mechanics often have the spacetime correlation property, which should be inherited and preserved in the corresponding numerical algorithms. Starting from the fundamental formulation of fluid mechanics under continuum hypothesis, this paper defines the meaning of spacetime correlation of the models, establishes the fundamental principle of finite volume schemes, expounds the necessity of spacetime coupling of algorithms, as well as realizes the physical and mathematical unification of basic governing equations of fluid mechanics and finite volume schemes. In practice, the design methodology of spacetime coupling high order numerical algorithms is presented, and the difference from spacetime decoupling method is compared. It should be pointed out that most of the contents in this paper are suitable for computational fluid dynamics under the assumption of continuous medium, and some are only suitable for compressible flow. 参考文献 【1】 华罗庚,大哉数学之为用,人民日报, 年 月 日 . 【2】 M. Ben-Artzi and J. Q. Li, Consistency of finite volume approximations to nonlinear hyperbolic balance laws, Math. Comp. 2020, doi:10.1090/mcom/3569. 【3】
M. Ben-Artzi and J. Q. Li, Regularity of fluxes in nonlinear hyperbolic balance laws, 2020, arXiv:2006.09253. 【4】
J. Q. Li, Two-stage fourth order: Temporal-spatial coupling in computational fluid dynamics (CFD),
Advances in Aerodynamics, 【5】
M. Ben-Artzi and J. Falcovitz,
A second-order Godunov-type scheme for compressible fluid dynamics,
J. Comput. Phys.
55 (1984), no. 1, 1–32. 【6】
M. Ben-Artzi, J. Q. Li and G. Warnecke, A direct Eulerian GRP scheme for compressible fluid flows,
J. Comput. Phys.
218 (2006), no. 1, 19–43. 【7】
J. Q. Li and Y. Wang, Thermodynamical effects and high resolution methods for compressible fluid flows,
J. Comput. Phys.
343 (2017), 340–354. 【8】
S. -C. Chang, The method of spacetime conservation element and solution element—A new approach for solving the Navier-Stokes and Euler equations,
J. Comput. Phys., 119(1995), 295-324. 【9】
K. Xu, A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method,
J. Comput. Phys.
171 (2001), no. 1, 289–335. 【10】
K. Xu and J. -C. Huang, A unified gas-kinetic scheme for continuum and rarefied flows,
J. Comput. Phys.
229 (2010), no. 20, 7747–7764. 【11】
A. L. Cauchy, Recherches sur l’équilibre et le mouvement intérieur des corps solides ou fluides, élastiques or non élastiques,
Bull. Soc. Philomathique , 10(2) (1823), 9–13. 【12】
A. L. Cauchy, Da la pression ou tension dans un corps solide , Exercises de Matehématiques , 2(2)(1827), 42–56 【13】
H. Federer, The Gauss–Green theorem,
Trans.Am. Math. Soc.
58 (1945), 44–76. 【14】
H. Federer,
Geometric Measure Theory , Springer-Verlag , 1969. 【15】
M. E. Gurtin and L. C .Martins, Cauchy’s theorem in classical physics,
Arch. Rat. Mech. Anal., 【16】
M. Šilhavý, The existence of the flux vector and the divergence theorem for general Cauchy fluxes,
Arch. Rat. Mech.Anal.,
90 (1985), 195–212. 【17】
G. -Q. Chen, G.E. Comi and M. Torres, Cauchy fluxes and Gauss-Green formulas for divergence measure fields over general open sets,
Arch. Rat. Mech.Anal.,
233 (2019), 87–166. 【18】
C. M. Dafermos, “Hyperbolic Conservation Laws in Continuum Physics, Fourth Edition”,
Grundlehren der Mathematischen Wissenschaften, vol. 325,
Springer, 【19】
E. Hopf, The partial differential equation 𝑢 ! + 𝑢𝑢 " = 𝜇𝑢 "" , Comm. Pure Appl. Math., 【20】
J. Bec and K. Ihanin, Burgers Turbulence,
Physics Reports , 447 (2007), 1-66. 【21】
T. Tang and Z. H. Teng, The sharpness of Kuznetsov's
𝒪d√Δ𝑥e 𝐿 / -error estimate for monotone difference schemes, Math. Comp . 64 (1995), 581-589. 【22】
Z. F. Du and J. Q. Li, A Hermite WENO reconstruction for fourth order temporal accurate schemes based on the GRP solver for hyperbolic conservation laws,
J. Comput. Phys.,
355 (2018), 385–396. 【23】
P.D. Lax and B. Wendroff, Systems of conservation laws.
Comm. Pure Appl. Math.
13 (1960), 217–237. 【24】
R. Eymard, T. Gallouet and R. Herbin,
Finite Volume Methods , in
Handbook of Numerical Analysis , Vol. VII, Eds. P.G. Ciarlet and J.-L. Lions, North-Holland, 2000, 713--1020. 【25】
V. Elling, A Lax-Wendroff type theorem for unstructured quasiuniform grids,
Math. Comp. , 76 (2007), 251--272. 【26】
S. K. Godunov, A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics,
Mat. Sb. (N.S.) , :3 (1959), 271–306. 【27】 J. Q. Li, T. Zhang and S. L. Yang, The two-dimensional Riemann problem in gas dynamics. Pitman Monographs and Surveys in Pure and Applied Mathematics, 98.
Longman, Harlow, 【28】
P. H. Maire, R. Abgrall, J. Breil, J. Ovadia, A cell-centered Lagrangian scheme for two-dimensional compressible flow problems,
SIAM Journal on Scientific Computing
29 (2007), 1781-1824. 【29】
J. Q. Li and Z. F. Du, A two-stage fourth order time-accurate discretization for Lax-Wendroff type flow solvers I. Hyperbolic conservation laws,
SIAM J. Sci. Comput.
38 (2016), no. 5, A3046–A3069. 【30】
Z. L. Guo, K. Xu, and R. J. Wang, Discrete unified gas kinetic scheme for all Knudsen number flows: Low-speed isothermal case,
Phys. Rev. E , 88 (2013), 033305. 【31】
L. Evans, Partial differential equations. Second edition. Graduate Studies in Mathematics, 19.
American Mathematical Society, Providence, RI, 【32】
M. Ben-Artzi and J. Q. Li, Hyperbolic balance laws: Riemann invariants and the generalized Riemann problem,
Numer. Math.
106 (2007), no. 3, 369–425. 【33】
E. Toro and V. Titarev, Derivative Riemann solvers for systems of conservation laws and ADER methods,
J. Comput. Phys.
212 (2006), no. 1, 150–165. 【34】
J. Z. Qian, J. Q. Li and S. H. Wang, The generalized Riemann problems for compressible fluid flows: towards high order,
J. Comput. Phys.
259 (2014), 358–389. 【35】
G. Montecinos, C. E. Castro, M. Dumbser and E. Toro, Comparison of solvers for the generalized Riemann problem for hyperbolic systems with source terms,
J. Comput. Phys.
231 (2012), no. 19, 6472–6494. 【36】
M. Ben-Artzi and J. Falcovitz, Generalized Riemann problems in computational fluid dynamics. Cambridge Monographs on Applied and Computational Mathematics, 11.
Cambridge University Press, Cambridge, 【37】
M. Ben-Artzi and J. Q. Li, Hyperbolic balance laws: Riemann invariants and the generalized Riemann problem,
Numer. Math.
106 (2007), no. 3, 369–425. 【38】
X. Lei and J. Q. Li, Transversal effects of high order numerical schemes for compressible fluid flows,
Applied Mathematics and Mechanics (English Edition) , 40(2019), 343–354. 【39】
J. Q. Li, On the two--dimensional gas expansion for compressible Euler equations,
SIAM Journal on Applied Mathematics , 62 (2001/2002), 831--852. 【40】
J. Q. Li and Y. X. Zheng, Interaction of rarefaction waves of the two-dimensional self-similar Euler equations,
Archive Rational Mechanics and Analysis , 193 (2009), 623--657. 【41】
J. Q. Li and Y.X. Zheng, Interaction of four rarefaction waves in the bi-symmetric class of the two-dimensional Euler equations,
Communications in Mathematical Physics,
296 (2010), 303--321. 【42】
D. Balsara and B. Nkonga, Multidimensional Riemann problem with self-similar internal structure—part III—a multidimensional analogue of the HLLI Riemann solver for conservative hyperbolic systems,
J. Comput. Phys.
346 (2017), 25–48. 【43】
J. Q. Li and Z. F. Sun, Remark on the generalized Riemann problem method for compressible fluid flows,
J. Comput. Phys.,
222 (2007), 796--808. 【44】
M. R. Baer, J. W. Nunziato, A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials,
Int. J. Multiphase Flow
12 (6) (1986) 861–889. 【45】 X . Lei and J. Q. Li, A non-oscillatory energy-splitting method for the computation of compressible multi-fluid flows,
Physics of Fluids , 30 (2018), 006891. 【46】
X. Lei, Z. F. Du and J. Q. Li, The simulation of compressible multi-fluid flows by a GRP-based energy-splitting method,
Computers and Fluids, 【47】
X. Lei and J. Q. Li, A staggered-projection Godunov-type method for the Baer-Nunziato two-phase model, arXiv:2007.14355v1,2020. 【48】
Y. Wang, Two relaxation methods for fluid dynamical systems, Ph.D Thesis, Beiijng Normal University, 2012. 【49】
G. Monecinos and E. Toro, Reformulations for general advection-diffusion-reaction equations and locally implicit ADER schemes,
J. Comput. Phys.
275 (2014), 415–442. 【50】
P. L. Bhatnagar, E. P. Gross, and M. Krook, A model for collision processes in gases I: Small amplitude processes in charged and neutral one-component systems,
Phys. Rev . 94, 511 (1954). 【51】
T. J. Poinsot and S. K. Lee, Boundary conditions for direct simulations of compressible viscous flows,
J. Comput. Phys.
101 (1992), 104–129. 【52】
J. Lu, J. Fang, S. Tan, C. Shu and M. Zhang, Inverse Lax-Wendroff procedure for numerical boundary conditions of convection-diffusion equations,
J. Comput. Phys.
317 (2016), 276–300. 【53】
J. Q. Li and Q.L. Zhang, One-sided GRP solver and high order numerical boundary conditions for compressible fluid flows, preprint, 2020. 【54】
Z. F. Du and J. Q. Li, Accelerated piston problem and high order moving boundary tracking method for compressible fluid flows,
SIAM J. Scientific Computing , 42(3) (2020), A1558-A1581. 【55】
Z. F. Du and J. Q. Li, A two-stage fourth order time-accurate discretization for Lax-Wendroff type flow solver, II High order numerical boundary conditions,
Journal of Computational Physics , 369 (2018), 125-147. 【56】
J. Q. Li and Y.J. Zhang, Numerical boundary treatment over unstructured meshes for compressible fluid flows, work in preparation, 2020. 【57】
J. Cheng, Z. F. Du, X. Lei, Y. Wang and J. Q. Li, A two-stage fourth-order discontinuous Galerkin method based on the GRP solver for the compressible Euler equations,
Computers and Fluids, 【58】
F. X. Zhao, X. Ji, W. Shyy and K. Xu, Compact higher-order gas-kinetic schemes with spectral-like resolution for compressible flow simulations,
Adv. Aerodyn.
1, 13 (2019). https://doi.org/10.1186/s42774-019-0015-6. 【59】
Y. H. Yuan and H. Z. Tang, On the explicit two-stage fourth-order accurate time discretizations, 2020 , arXiv:2007.02488. 【60】 L. Pan, K. Xu, Q. B. Li and J. Q. Li, An efficient and accurate two-stage fourth-order gas-kinetic scheme for the Euler and Navier-Stokes equations , Journal of Computational Physics , Vol 326 (2016), 197-221. 【61】
Y. Chen, C. Chen, F. Xiao, X. Li and X. Shen, A two-stage fourth-order multi-moment global shallow water model on cubed sphere, Mon. Wea. Rev. 1-40 (2020), https://doi.org/10.1175/MWR-D-20-0004.1. 【62】
C. Zhang, Q. B. Li, Z. J. Wang, J. Q. Li and S. Fu, A two-stage fourth-order gas-kinetic CPR method for Navier-Stokes equations on triangular meshes, preprint, 【63】
L. Pan, J. Q. Li and K. Xu, A Few b enchmark test cases for higher-order Euler solvers, Numerical Mathematics: Theory, Methods and Applications, Vol.10(2017), 489-514. 【64】
H. Z. Tang and T. G. Liu, A note on the conservative schemes for the Euler equations,
J. Comput. Phys.,
218 (2006), pp. 451–459. 【65】
V. A. Titarev and E. F. Toro, Finite volume WENO schemes for three-dimensional conservation laws,
J. Comput. Phys.
201 (2014), 238–260. 【66】
J. J. Quirk and S. Karni, On the dynamics of a shock-bubble interaction,
J. Fluid Mech. , 318(1996), 129-163. 【67】
G. Y. Cao, L. Pan, K. Xu, Three dimensional high-order gas-kinetic scheme for supersonic isotropic turbulence I: criterion for direct numerical simulation,
Computers & Fluids,
192 (2019), 104273. 【68】
R. Samtaney, D. I. Pullin, B. Kosovic, Direct numerical simulation of decaying compressible turbulence and shocklet statistics,
Physics of Fluids , 13 (2001),1415-1430. 【69】
Z. L. Guo, J. Q. Li and K. Xu, On unified preserving properties of kinetic schemes, preprint, 【70】
G. W. He, G. D. Jin, and Y. Yang, Spacetime correlations and dynamic coupling in turbulent flows,