摘要:不可压Navier-Stokes方程是流体力学的基本控制方程, 其高精度数值模拟具有重要的科学意义. 本综述性文章回顾了求解Navier-Stokes方程的投影方法, 重点介绍了时空一致四阶精度的GePUP方法. 该方法用一个广义投影算子对不可压Navier-Stokes方程进行了重新表述, 使得投影流速的散度由一个热方程控制, 保持了UPPE方法的优点. 与UPPE方法不同的是, GePUP方法的推导不依赖于Leray-Helmholtz投影算子的各种性质, 并且GePUP表述中的演化变量无需满足散度为零的条件, 因此数值近似Leray-Helmholtz投影算子的误差对精度和稳定性的影响非常透明. 在GePUP方法中, 时间积分和空间离散是完全解耦的, 因此对这两个模块都能以“黑匣子”的方式自由替换. 时间积分模块的灵活性实现了时间上的高阶精度, 并使得GePUP算法能同时适用于低雷诺数流体和高雷诺数流体. 空间离散模块的灵活性使得GePUP算法能很好地适应不规则边界. 理论分析和数值测试结果都显示, 相对于二阶投影方法, GePUP方法无论在精度上还是效率上都具有巨大优势.
关键词: 不可压Navier-Stokes方程/
投影方法/
时空一致四阶精度/
无滑移边界条件/
广义投影算子/
压力Poisson方程

English Abstract

全文HTML

--> --> -->

1900年, 著名数学家希尔伯特在国际数学联盟IMU (International Mathematical Union)巴黎会议上提出了23个数学问题, 这个列表极大地影响了20世纪的数学发展. 1998年, 菲尔兹奖得主Smale^[1]应IMU的邀请, 给出了18个他认为在21世纪非常重要的数学问题. 2000年5月, 克雷数学研究所发布了七个千禧年数学问题^[2], 由于这些问题的重要性, 解决其中任何一个都能获得100万美元的奖金. 在千禧年数学问题和Smale的18个问题中, 有一个问题同时存在, 它就是三维空间中Navier-Stokes方程全局解的存在性和正则性.
数学家们对Navier-Stokes方程的重视源于该方程极其广泛的应用, 尤其是该方程的解对于理解湍流很可能具有根本性的重要意义. 如果一个牛顿流体是性质均匀且各向同性的, 其不可压Navier-Stokes方程 (INSE)的基本形式为

$ \frac{\partial {u}}{\partial t} = -{u}\cdot\nabla{u} +{g} -\nabla p + \nu\Delta {u},~~ \text{in} \;\varOmega\times (0, +\infty),\tag{1a} $

$ \nabla \cdot {u} = 0, \quad \text{in}\; \varOmega\times (0, +\infty),\tag{1b} $

其中t是时间变量, $ {x}\in{\mathbb{R}}^ { {\rm{D}}} $是空间变量($ { {\mathrm{D}}} = 2, 3 $), $ {g} $是外力项, $ \nu $是动力黏性系数, $ p({x}, t) $是压强, 而$ {u}({x}, t) $是流速. 假定流速场的初始条件${u}_0\in {\cal C}^{\infty}$光滑无源, 并且其空间导数满足一定的有界条件, 对于外力项的空间时间偏导数也做类似假定. 在千禧年问题中, 空间域$ \varOmega $是$ {\mathbb{R}}^{3} $或3-环面($ {\mathbb{R}}^3/{\mathbb{Z}}^3 $), 因此$ \varOmega $没有边界. 但是在实际物理问题中, 往往需要指定(1)式中的边界条件. 除周期边界条件外, 本文主要讨论最普遍也是最重要的无滑移边界条件:

$ {u} = {0}, \quad \text{on }\; \partial \varOmega. \tag{2} $

在无滑移边界条件下, 中高雷诺数不可压流体往往含有多个时间尺度和多个空间尺度的流动结构. 因此在进行数值模拟时, 数值方法也必须拥有足够的精度来捕捉所有相关尺度的特征, 这对计算效率提出了极大的挑战. 一阶和二阶方法的计算结果虽然可以收敛, 但是为了达到给定的计算精度, 有限的计算资源可能会很快地消耗殆尽. 另外, 至关重要的湍流特征往往取决于流速场的梯度或高阶导数. 例如, 涡量是流速的一阶导数, 而转捩点是切向流速的二阶法向导数为零的位置. 低阶方法的流速场计算结果虽然可以收敛, 但是对于流速导数决定的重要物理量往往存在较大的误差. 综上所述, 我们希望发展一个时空一致高阶方法, 达到以下五个目标:
(Goal-1) 在时间上和空间上同时达到四次收敛阶;
(Goal-2) 在每个时间步内只需求解常数个关于压强或关于流速的线性方程组, 并且方程组求解的复杂度为$ O(N) $, 这里N是控制体个数;
(Goal-3) 时间积分方法的具体细节不出现在整体算法中;
(Goal-4) 时间积分与空间离散完全解耦, 也就是说用户可以在一定范围内自由选择空间离散和时间积分模块组合成一个整体;
(Goal-5) 同时适用于高雷诺数和低雷诺数流体.
众所周知, 当动力黏性系数很大时, INSE的半离散系统具有刚性, 需要对(1a)式中的扩散项采用隐式时间积分, 否则时间步长的选取将受到数值稳定性非常苛刻的限制. 而当动力黏性系数很小时, 对流占优使得INSE的半离散系统不具有刚性, 此时用显式时间积分方法也可以使时间步长免受稳定性的限制. 如果不能满足(Goal-3), 那么针对不同情况更换时间积分方法将是一件非常繁复的任务. 反之, 如果能以“黑匣子”的方式调用高阶时间积分模块, 将会实现(Goal-4), 同时能为整体算法提供最大的灵活性, 从而轻松实现(Goal-5). 因此, (Goal-5)可以看成是(Goal-3)和(Goal-4)的副产品.
在第2节预备知识和第3节经典投影方法的基础上, 本文重点介绍了一个满足上述所有目标((Goal-1)—(Goal-5))的四阶投影方法GePUP^[3—5], 该方法达到核心目标((Goal-3)和(Goal-4))的关键是策略性地利用投影算子处理不可压约束条件(1b)式. 第4节推导INSE的另一种表述, 作为GePUP方法的控制方程, 并简要总结了在有限体积框架下的数值算法. 由于(Goal-4)带来的灵活性, 用针对不规则区域的离散格式替换经典有限体积离散格式就得到了一个时空一致四阶精度, 适用于不规则区域的投影方法. 第5节中的数值实验结果验证了GePUP算法的时空一致四阶精度及其捕捉真实物理过程的能力, 其中5.4节系统地论述了GePUP算法相对二阶方法^[6,7]在精度和效率上的优越性.

2

-->

2.1.Helmholtz分解定理

与流速场不同, 动量方程(1a)中的压强p既不需要初始条件也不需要边界条件. 这是因为Helmholtz分解定理指出, 有界区域$\varOmega\subseteq{\mathbb{R}}^{ { {{D}}}}$上任何充分光滑的向量场${v}^*$都可以被唯一地分解为一个无源场${v}$和一个无旋场$ \nabla \phi $之和:

$ {v}^* = {v} + \nabla \phi, \tag{3a}$

$ \nabla\cdot {v} = 0, \qquad \nabla\times\nabla \phi = {0}. \tag{3b} $

该分解可以通过求解Neumann边界条件下的Poisson方程实现:

$ \Delta \phi = \nabla\cdot{v}^* \qquad \quad \quad \ \text{ in }\; \varOmega,\tag{4a} $

$ {n}\cdot \nabla\phi = {n}\cdot({v}^* - {v}) \qquad \text{ on }\; \partial\varOmega, \tag{4b}$

其中$ {n} $表示区域边界$ \partial \varOmega $的单位外法向量.
如果定义欧拉加速度

$ {a}: = \frac{\partial {u}}{\partial t},\qquad {a}^*: = -{u}\cdot\nabla{u} +{g} + \nu\Delta {u}, $

则可将动量方程(1a)写成(3a)式的形式:

$ {a}^* = {a} + \nabla p. $

因此压强可由欧拉加速度$ {a}^* $在相差一个常数的意义下由Helmholtz定理唯一确定, 并且在这个确定过程中无需任何p的边界条件, 或者说p的边界条件也是由欧拉加速度$ {a}^* $和$ {a} $在边界上的法向分量确定的. 另外一个从流速场$ {u} $得到压强的方法是对动量方程(1a)求散度, 应用(1b)式后可以得到关于压强的Poisson方程, 从而求解出压强及其梯度.
2 -->

2.2.Leray-Helmholtz投影算子

下面的定理在很多偏微分方程的教科书上都可以找到, 如文献[8]:
定理1　如果f和g足够光滑, 具有Neumann边界条件的Poisson方程

$ \Delta \phi = f \; \text{ in }\;\varOmega,\qquad {n}\cdot\nabla \phi = g \;\text{ on }\;\partial\varOmega $

解存在的充要条件是

$ \int_{\varOmega} f = \oint_{\partial \varOmega} g. $

并且$ \phi $在相差一个常数的意义下是唯一确定的.
容易验证(4)式满足可解条件(8)式且$ \nabla \phi $唯一确定. 因此对于任意给定的向量场$ {v}^* $都可以唯一确定(3a)式中的无源场$ {v} $. 由这个唯一性可以定义一个线性算子${\mathscr{P}}$, 称为Leray-Helmholtz投影算子, 从而把Helmholtz定理中的分解步骤整合其中,

$ {\mathscr{P}} {v}^*: = {v} = {v}^* -\nabla\phi. $

由(3b)式可知, 在任意时刻 $ {\mathscr{P}} $都满足

$ {\mathscr{P}}^2 = {\mathscr{P}}, \quad \nabla\cdot {\mathscr{P}}{v}^* = 0,\quad {\mathscr{P}}\nabla\phi = {0}. $

将Leray-Helmholtz投影算子作用在欧拉加速度$ {a}^* $上可以得到$ {a} $, 从而得到压强梯度$ \nabla p = ({\cal I}- {\mathscr{P}}){a}^* $, 这个压强梯度的提取方式与时间无关.
在周期边界条件下, Leray-Helmholtz投影算子${\mathscr{P}}$和Laplace算子$ \Delta $是可交换的, 将${\mathscr{P}}$作用于(1a)式, 再结合(10)式中最后一个性质可得

$ \frac{\partial {u}}{\partial t} = {\mathscr{P}}\left(-{u}\cdot\nabla{u} +{g}\right) + \nu\Delta {u}. $

在应用Method-of-lines (MOL)方法数值求解INSE时, 由于(1b)式不是一个动力过程而是一个限制条件, 直接对(1)式进行离散得到的常微分方程组很难在求解过程中保证收敛阶^[9]. 相比来说, (11)式是关于流速场的一个动力过程, 可以直接应用MOL方法. 这是(11)式相对于(1)式的一个显著优势.
2 -->

2.3.Laplace-Leray交换子

在无滑移条件下, Leray-Helmholtz投影算子${\mathscr{P}}$和Laplace算子$ \Delta $是不可交换的. 为了度量这两者的不可交换性并推导出不带约束条件的关于流速场的动力过程, 定义Laplace-Leray交换子为

$ [\Delta, {\mathscr{P}}]: = \Delta {\mathscr{P}} - {\mathscr{P}} \Delta. $

为了推导其表达式, 先定义梯度-散度交换子

$ {\cal B} = [\nabla, \nabla\cdot]: = \Delta - \nabla \nabla\cdot, $

满足

$ \nabla\cdot{\cal B} = 0,\tag{14a} $

$ \Delta {\mathscr{P}} = {\cal B},\tag{14b} $

其中(14a)式是由(13)式以及$ \Delta $和$ \nabla\cdot $的交换性得到的, 而(14b)式是由(9)式, $ \Delta $和$ \nabla $的交换性, (4a)式以及(13)式得到的:

$\begin{split} \Delta {\mathscr{P}} {v}^* =\;& \Delta({v}^* - \nabla\phi) = \Delta {v}^* - \nabla\Delta\phi \\=\;& \Delta {v}^* - \nabla \nabla\cdot {v}^* = {\cal B} {v}^*. \end{split} $

应用${\mathscr{P}}$于(15)式的第一和第三项, 结合(10)式中最后一个性质, 可得

$\begin{array}{l} {\mathscr{P}} \Delta {\mathscr{P}} {v}^* = {\mathscr{P}} \Delta {v}^*. \end{array} $

因此Laplace-Leray交换子满足

$\begin{array}{l} [\Delta, {{\mathscr{P}}}]: = \Delta {{\mathscr{P}}} - {{\mathscr{P}}} \Delta = ({\cal I}- {{\mathscr{P}}})\Delta {{\mathscr{P}}} = ({\cal I}- {{\mathscr{P}}}) {\cal B}, \end{array} $

其中$ {\cal I} $表示恒等算子. 由(9)式知 $({\cal I}- {\mathscr{P}}) {v}^* = \nabla\phi$, 因此$[\Delta, {\mathscr{P}}]{v}^*$是某标量场的梯度, 若令$ {v}^* $为(1)式中的流速$ {u} $, 得到的标量场称为 Stokes压强, 其梯度满足

$\begin{array}{l} \nabla p_{\rm s}: = (\Delta {\mathscr{P}} - {\mathscr{P}} \Delta) {u}. \end{array} $

由(17)式和(18)式可得

$\begin{array}{l} {\cal B}{u} = {\mathscr{P}} {\cal B}{u} + \nabla p_{\rm s}. \end{array} $

结合(10)式中第二个性质以及(14a)式可得$ \Delta p_{\rm s} = 0 $, 即Stokes压强既无旋也无源, 是调和的.
定义标量$p_{\rm c}$满足

$\begin{array}{l} \nabla p_{\rm c}: = ({\cal I} - {\mathscr{P}})\left({g} -{u}\cdot\nabla{u}\right). \end{array} $

应用Leray-Helmholtz投影算子于(1a)式, 再由(6)式, (17)式, (18)式及(20)式可得

$\begin{array}{l} \nabla p = \nabla p_{\rm c} + \nu \nabla p_{\rm s}. \end{array} $

因此INSE中的压强梯度由两部分构成, 其中$ \nabla p_{\rm c} $平衡外力项和非线性对流项的散度, 而$ \nabla p_{\rm s} $度量了 Laplace算子和Leray-Helmholtz投影算子的非交换性. 在$ \nu\rightarrow 0 $和$ \nu\rightarrow +\infty $这两种极限情况下, 压强梯度相应地由$ \nabla p_{\rm c} $和$ \nu \nabla p_{\rm s} $所主导, 因而INSE的性质在这两种情况下会有很大的差别.

2

-->

3.1.鞍点法

鞍点法^[10]将 (1)式中的流速和压强耦合起来同时求解. 例如, 由Crank-Nicolson方法可得

$ \begin{split}&\frac{{u}^{n+1}-{u}^n}{\Delta t}+ {{G}} p^{n+\tfrac{1}{2}} \\=\;& -{N}\big({u}^{n+\tfrac{1}{2}}\big)+{g}^{n+\tfrac{1}{2}}+\nu {{L}}{u}^{n+\tfrac{1}{2}}, \end{split}\tag{22a}$

$ {{D}}{u}^{n+1} = 0, \tag{22b}$

其中, ${u}^{n+\tfrac{1}{2}} = \dfrac{1}{2}\left({u}^{n+1}+{u}^n\right)$, ${N}\big({u}^{n+\tfrac{1}{2}}\big)\approx $$ \big[\big({u}^{n+\tfrac{1}{2}}\cdot\nabla \big){u}^{n+\tfrac{1}{2}}\big]$是对非线性对流项的二阶近似. 将(22)式写成一个线性方程组:

$ { A} \begin{bmatrix} {u}^{n+1} \\ p^{n+\tfrac{1}{2}} \end{bmatrix} = {F}. $

由于矩阵A具有鞍点结构, 这个方法被称为“鞍点法”. 同时求解流速和压强使得线性方程组中的未知量个数最大化, 因此这种方法对方程组的求解效率提出了很大的挑战^[11]. 另外, (23)式中矩阵A的推导依赖于空间离散和具体时间积分方法的细节. 对于较为复杂的时间积分方法, 如高阶Additive Runge-Kutta (ARK)方法, 矩阵A的复杂性往往是难以接受的. 因此, 这个方法难以实现(Goal-3)和(Goal-4), 从而也难以做到(Goal-1)和(Goal-2).
2 -->

3.2.Chorin的一阶投影方法

1968年, Chorin^[12]提出了以下一阶投影方法:

$ \frac{ {u}^*-{u}^n}{\Delta t} = -{C}({u}^*, {u}^n)+{g}^n+\nu {{L}}{u}^*, \tag{24a}$

$ {u}^{n+1} = {{P}}{u}^*, \tag{24b}$

其中$ {C}({u}^*, {u}^n) $是对对流项的近似, ${{P}} = {I}- $$ {{G}} {{L}}^{-1} {{D}}$是对Leray-Helmholtz投影算子的一个近似. 该方法首先在忽略压强梯度项的情况下计算一个辅助流速场$ {u}^* $, 再将$ {u}^* $投影到散度为零的空间上来实现不可压约束(1b)式. 由于流速和压强在计算域内解耦, 在每个时间步内, 都只需要求解一系列关于流速或压强的边值问题, 因此线性方程组的规模较鞍点法会小很多, 这种效率优势可能是投影方法得到广泛应用的重要原因之一.
2 -->

3.3.分步二阶投影方法

在Chorin的基础上, ****们提出了许多二阶投影方法^[13,6,14-17], 出发点是先对INSE进行时间上的二阶隐式离散:

$ \begin{split}&\frac{{u}^{n+1}-{u}^n}{\Delta t}+\nabla p^{n+1/2}= \\ \;& -[({u}\cdot\nabla){u}]^{n+ \tfrac12} \!+\! {g}^{n+ \tfrac{1}{2}} \!+\! \frac{\nu}{2}\Delta({u}^{n+1}\!+\! {u}^n),\qquad \end{split} \tag{25a} $

$ \nabla\cdot {u}^{n+1} = 0, \tag{25b}$

这里$ [({u}\cdot\nabla){u}]^{n+1/2} $表示对流项在时刻$ t^{n+1/2} $的显式二阶近似. 如果再进行空间离散就得到了鞍点法的(22)式. 与Chorin的投影方法类似, 可以引入辅助变量将压强和流速解耦. 不同的是, 达到二阶精度需要引入两个辅助变量.
定义2　分步二阶投影方法由以下三个步骤构成^[17]:
(Proj-1) 选择一个边界条件$ {B}({u}^*) = {0} $以求解辅助变量$ {u}^* $:

$\begin{split} \frac{{u}^*-{u}^n}{\Delta t}+\nabla q =\;& -[({u}\cdot\nabla){u}]^{n+1/2}+{g}^{n+{1}/{2}}\\&+\frac{\nu}{2}\Delta({u}^*+{u}^n), \end{split} \tag{26} $

(Proj-2) 利用投影算子实现不可压约束(25b)式:

$ {u}^* = {u}^{n+1}+\Delta t\nabla\phi^{n+1}, \tag{27a}$

$ \nabla\cdot {u}^{n+1} = 0. \tag{27b}$

注意$ \phi^{n+1} $的边界条件要与$ {u}^{n+1} $的物理边界条件以及$ {u}^* $的边界条件$ {B}({u}^*) = {0} $相容.
(Proj-3) 更新压强$ p^{n+1/2} $:

$ p^{n+1/2} = q+{U}(\phi^{n+1}). $

一个具体的分步二阶投影方法由以上步骤中的三个要素确定, 即压强近似函数q、辅助变量$ {u}^* $的边界条件$ {B}({u}^*) $以及压强更新函数$ {U} $. 其中投影步骤(Proj-2)通过求解Neumann边界条件下的Poisson方程实现:

$ \Delta t\Delta\phi^{n+1} = \nabla\cdot{u}^*, \tag{29a} $

$ \Delta t\frac{\partial\phi^{n+1}}{\partial{n}} = {n}\cdot{u}^{*}|_{\partial\varOmega}. \tag{29b} $

将(27a)式代入(26)式, 结合(25a)式知(Proj-3)步骤中的压强更新公式(28)式可写为

$ p^{n+1/2} = q+\phi^{n+1}-\frac{\nu\Delta t}{2}\Delta\phi^{n+1}. $

以下简要回顾在定义2的框架下的两个经典的二阶投影方法.
3 -->

3.3.1.Bell, Colella & Glaz (BCG)

Bell, Colella和Glaz^[6]指定

$ \begin{split} &q = p^{n-1/2}, \quad {B}({u}^*) = {u}^*|_{\partial\varOmega} = {0}, \quad \\ &\nabla p^{n+1/2} = \nabla p^{n-1/2}+\nabla\phi^{n+1}. \end{split} $

在文献[17]中, 作者运用正则模式分析和数值实验验证了在$ L_{\infty} $范数度量下, 该方法对于流速的求解是二阶收敛的, 而压强的计算结果则仅有一阶收敛. 造成压强降阶的原因是(31)式中对压强梯度的更新公式与(30)式是不匹配的, 忽略(30)式中的高阶项会导致压强p的精度损失(通常会出现边界层)以及压强梯度在边界上的不准确. 的确, (29b)式和(31)式意味着

${n}\cdot\nabla\phi^{n+1}|_{\partial\varOmega} = 0, \quad {n}\cdot\nabla p^{n+1/2} = {n}\cdot\nabla p^{n-1/2}.$

而这一般是不对的. 但若修改(31)式中的压强梯度更新公式使得其与(30)式一致, 即

$ \begin{split} &q = p^{n-1/2}, \quad {B}({u}^*) = {u}^*|_{\partial\varOmega} = {0}, \quad \\ &\nabla p^{n+1/2} = \nabla p^{n-1/2}+\nabla\phi^{n+1}-\frac{\nu\Delta t}{2}\nabla\Delta\phi^{n+1}. \end{split} $

则对流速和压强的求解都能达到二阶精度^[17].
Martin, Colella和Graves提出了BCG方法的改进版, 称为MCG方法, 见5.4节.
3 -->

3.3.2.无压投影

在Kim和Moin^[13]提出的无压投影方法中, 压强近似变量q不出现在辅助变量$ {u}^* $的求解过程(Proj-1)中:

$ q = 0, \tag{33a} $

$ {n}\cdot{u}^*|_{\partial\varOmega} = 0, \quad {\tau}\cdot{u}^*|_{\partial\varOmega} = \Delta t {\tau}\cdot\nabla\phi^n|_{\partial\varOmega}, \tag{33b} $

$ \nabla p^{n+1/2} = \nabla\phi^{n+1} - \frac{\nu\Delta t}{2}\nabla\Delta\phi^{n+1}, \tag{33c} $

其中$ {\tau} $为$ \partial\varOmega $的单位切向量, 而$ q = 0 $这一选择导致$ {u}^* $与$ {u}^{n+1} $之间差别较大; 如果达到二阶精度, 其边界条件不能简单地提为$ {u}^*|_{\partial\varOmega} ={0} $. 由于$ {n}\cdot{u}^{n+1} $的正确性总能由投影步骤(Proj-2)保证, 因此从理论上说$ {u}^* $的法向分量$ {n}\cdot{u}^* $的选择具有很大的灵活性. 但从数值计算角度看, $ {u}^* $的边界条件选取方式会影响$ {u}^* $在靠近计算域边界处的性质, 也会对压强p在边界处的性质产生影响. 的确, (29a)式和(33c)式意味着

$ p^{n+1/2} = \phi^{n+1}-\frac{\nu}{2}\nabla\cdot{u}^*.$

在文献[17]中, Brown等用数值实验展示了除非选择$ {u}^* $的边界条件使得其在靠近边界处光滑, 否则该方法求解压强时不能达到二阶精度.
分步二阶投影方法虽然取得了巨大的成功, 但要直接将其推广到更高阶精度是非常困难的, 主要表现在定义2中三个不确定因素的耦合: 若改变时间积分方法, 则三个因素的选取方式全都需要重新推导, 这无疑是极其繁琐的. 其次, 辅助变量$ {u}^* $不具有物理意义, 其边界条件必须通过具体时间积分方法的细节来进行推导. 另外, 压强和流速尽管在计算域内部是解耦的, 但它们在边界上仍然通过辅助变量$ {u}^* $耦合在一起. 综上所述, 用分步二阶投影方法很难达到第1节中的五个目标. 也正是因为这些原因, 之前很多****把投影方法称为时间分步法, 并认为它是不可能达到四阶精度的.
2 -->

3.4.Gauge方法

在二阶投影方法中, 对辅助变量$ {u}^* $提恰当的边界条件通常比较困难. Weinan和Liu^[15] 提出了Gauge方法, 该格式的出发点是用规范变量$ \chi $来替代压强p, 他们还引入辅助变量

$ {m}: = {u}+\nabla\chi, $

并策略性地选择$ \chi $和$ {m} $使得由此导出的流速$ {u} $和压强p满足INSE. 考虑

$ {m}_t+({u}\cdot\nabla){u} = {g}+\nu\Delta{m}\quad \text{ in } \;\varOmega, \tag{35a}$

$ {u} = {0}\quad \text{ on } \;\partial\varOmega.\tag{35b} $

由(34)式, (35a)式和(1a)式可得

$ p = \chi_t-\nu\Delta\chi. $

(34)式和(35)式的优势在于可以用规范变量的自由度来对$ {m} $和$ \chi $提明确的边界条件, 例如

$ {n}\cdot\nabla\chi = 0, \;\; {n}\cdot{m} = 0, \;\; {\tau}\cdot{m} = {\tau}\cdot\nabla\chi \quad \text{ on }\; \partial\varOmega, $

或

$ \chi = 0,\; \; {n}\cdot{m} = {n}\cdot\nabla\chi,\; \;{\tau}\cdot{m} = 0 \quad \text{ on }\; \partial\varOmega. $

对控制方程(34)式, (35)式及边界条件(37)式离散后得到以下定义.
定义3　Gauge方法包括以下步骤:
(GM-1) 求解规范变量$ {m}^{n+1} $:

$\begin{split} \frac{{m}^{n+1}-{m}^n}{\Delta t} =\;& -[({u}\cdot\nabla){u}]^{n+1/2}+{g}^{n+1/2}\\&+\frac{\nu}{2}\Delta({m}^{n+1}+{m}^n), \end{split} $

注意, $ {m}^{n+1} $的边界条件要与$({m}^{n+1}\!-\!\nabla\chi^{n+1})|_{\partial\varOmega} = $$ {0}$相容,

$ {n}\cdot{m}^{n+1} = 0, \;\; {\tau}\cdot{m}^{n+1} = {\tau}\cdot\nabla(2\chi^n-\chi^{n-1}) \text{ on } \partial\varOmega, $

其中第二个等式是由外插公式($\chi^{n+1}\approx 2\chi^n- \chi^{n-1}$)近似后得到的.
(GM-2) 求解流速$ {u}^{n+1} $:

$\begin{array}{l} {u}^{n+1} = {m}^{n+1}-\nabla\chi^{n+1}, \end{array} $

其中$ \chi^{n+1} $满足

$ \Delta\chi^{n+1} = \nabla\cdot{m}^{n+1} \;\text{ in }\; \varOmega,\tag{42a} $

$ {n}\cdot\nabla\chi^{n+1} = 0 \;\text{ on }\; \partial\varOmega.\tag{42b} $

(GM-3) 如果需要的话, 压强p可由(36)式得到:

$ p^{n+1/2} = \frac{\chi^{n+1}-\chi^n}{\Delta t}-\frac{\nu}{2}\Delta(\chi^{n+1}+\chi^n). $

该方法对流速$ {u} $和压强p都能达到二阶收敛, 并且(40)式中使用的外插公式是必要的, 若仅使用滞后值$ \chi^n $则只有一阶精度^[15]. 由边界条件(40)式可知, 该方法与前面提到的方法具有同样的问题, 即具体算法依赖于时间积分方法的细节.
2 -->

3.5.PPE方法

Pressure Poisson Equation (PPE)方法^[18-24]的核心理念是用一个压力Poisson方程替换流速的不可压条件:

$ \frac{\partial {u}}{\partial t} = -{u}\cdot\nabla{u}+{g}-\nabla p+\nu\Delta{u} \quad \text{ in }\;\varOmega, \tag{44a}$

$ {u} = {0} \qquad \text{ on } \;\partial\varOmega, \tag{44b} $

$ \Delta p = \nabla\cdot({g}-{u}\cdot\nabla{u}+\nu\Delta{u}) \quad \text{ in }\; \varOmega,\tag{44c} $

$ {n}\cdot\nabla p = {n}\cdot({g}+\nu\Delta{u}) \quad \text{ on }\; \partial\varOmega, \tag{44d}$

其中(44c)式是对动量方程(1a)取散度后结合不可压约束(1b)式得到的, 而(44d)式是由(44a)式的法向分量及无滑移边界条件得到的. 若将$ \nabla\cdot{u} = 0 $作为(44)式的边界条件, 则(44)式等价于INSE^[19].
相对于前面四小节中的方法, PPE方法最大的优势在于可将时间积分模块完全看成“黑匣子”: 流速$ {u} $是(44)式中唯一的演化变量, 而压强p可视为$ {u} $的隐函数. 在任一时刻, 压强都可以按照(44c)式从流速中提取出来, 并且p和$ {u} $的边界条件不再与时间积分方法的细节相耦合. 当应用一个时间积分方法的时候, 只要给出了(44a)式右端项在该方法指定时刻的“样本”值, 就可以根据该方法的“契约”得到流速在时间上的高阶解. 在这个过程中, 无需知道该时间积分方法的任何细节.
PPE方法(44)式也存在一些不足之处. 首先, 其推导过程假设了流速的不可压条件总是成立的; 由于(44c)式的推导过程依赖于(1b)式, 要分析近似(1b)式产生的误差对INSE解的影响是很困难的. 更重要的是, (44a)式和(44c)式意味着 $\dfrac{\partial \nabla\cdot{u}}{\partial t} = 0$, 即流速散度的耗散是退化的, 也就是说在应用该方法时对流速散度是没有控制的. 数值实验结果显示, 直接对PPE方法的控制方程进行MOL离散得到的四阶方法是不稳定的.
2 -->

3.6.UPPE方法

利用2.3节中的Laplace-Leray交换子估计, Liu等^[22]提出了以下UPPE (Unconstrained PPE)方法:

$ \frac{\partial {u}}{\partial t} = {g}-{u}\cdot\nabla{u}-\nabla p+\nu\Delta{u}\;\text{ in }\; \varOmega,\tag{45a} $

$ {u} = {0}\; \text{ on } \;\partial\varOmega, \tag{45b} $

$ \Delta p = \nabla\cdot({g}-{u}\cdot\nabla{u}) \;\text{ in }\; \varOmega, \tag{45c} $

$ {n}\cdot\nabla p = {n}\cdot({g}+\nu\Delta{u}-\nu\nabla\nabla\cdot{u})\; \text{ on }\; \partial\varOmega. \tag{45d} $

该方法相对PPE方法的优势在于对(45a)式取散度后结合(45c)式可得

$ \frac{\partial(\nabla\cdot{u})}{\partial t} = \nu\Delta(\nabla\cdot{u}). $

因此流速的散度由一个热方程所控制, 根据热方程的极大值原理, UPPE格式(45)式对流速散度$ \nabla\cdot{u} $的控制要比原始PPE格式(44)式可靠得多. 在有限元离散的框架下, 数值实验结果^[23]印证了 UPPE方法的稳定性.

在把(45)式中的UPPE表述应用于基于MOL的有限差分和有限体积方法的时候我们遇到了困难. 首先, 直接对控制方程(45)进行MOL离散再求解常微分方程组得到的算法是不稳定的. 虽然在时间积分的过程中对流速的中间变量进行投影有时可以得到稳定的收敛结果, 但是由于控制方程(45)中并没有出现任何投影算子, 这个对流速中间变量进行投影的额外步骤是直接违背原方程的. 另外, 即使可以把控制方程(45)中的一些$ {u} $用$ {{\mathscr{P}}}{u} $进行替换, 也很难找到一个离散的四阶投影算子$ {{P}} $满足Leray-Helmholtz投影性质 (10)式. 换言之, 离散四阶投影$ {{P}} $难以满足以下任何条件:

$ {{P}}^2 = {{P}},\quad {{D}} {{P}}\left\langle {{u}} \right\rangle = 0,\quad {{P}} {{G}}\left\langle {\phi} \right\rangle = {0}, $

其中$ \left\langle {{u}} \right\rangle $和$ \left\langle {\phi} \right\rangle $是定义在控制体上的平均值, $ {{D}} $和$ {{G}} $分别是离散散度算子和离散梯度算子. 由于计算域的非周期性和四阶精度要求, 这一困难在交错网格上仍然存在.
UPPE表述(45)式的推导过程依赖于Leray-Helmholtz投影的性质(10)式. 在这种情况下, 如果仍然把不可压流速场作为动力系统的演化变量, 那么用一个不满足(47)式的离散投影算子$ \tilde{ {{P}}} $来近似Leray-Helmholtz投影算子产生的离散误差将对数值稳定性造成怎样的影响呢? 这个问题在UPPE表述下是难以回答的, 我们认为这也是导致UPPE表述的MOL离散格式不稳定的根本原因.
为了在理论上刻画不满足(47)式的离散投影算子, 4.1节定义了一个广义投影算子(generic projection). 在4.2节中, 我们首先推导Laplace算子和广义投影算子的交换子, 然后在4.3节中将它作用在INSE上得到GePUP表述. 这三小节理论推导背后的主要思想是, 既然满足(47)式的离散投影算子难以找到, 投影以后的流速场也不一定能完全满足不可压条件, 那么索性将散度不为零的流速场作为动力系统的主变量. 之后四小节内容讨论在GePUP表述基础上的离散算法. 其中4.4节介绍经典的有限体积离散, 4.5节列举了一些四阶离散算子的形式, 4.6节简述了如何把4.4节和4.5节中对于规则区域的空间离散推广到不规则区域上, 4.7节引入时间积分, 总结了GePUP-IMEX算法的主要步骤.
2

-->

4.1.广义投影算子

广义投影算子是向量空间上的一个线性算子, 其输入向量与输出向量之差为一个梯度场:

$ {\cal{P}} {v}^* = {v}: = {v}^* -\nabla\psi \ \;\text{ in } \;\varOmega, $

其中$ \nabla\cdot{v} = 0 $不一定成立. 对于非周期边界条件的有界区域, 还要求广义投影算子${\cal{P}}$保持输入流场的法向边界条件, 即

$ {n}\cdot {\cal{P}} {v}^* = {n}\cdot{v}^*\\ \text{ on } \;\partial \varOmega. $

无穿透边界上的Leray-Helmholtz投影算子显然是一个广义投影算子.
对一个给定的广义投影算子${\cal{P}}$, 由(48)式和(49)式知以下Poisson方程

$ \Delta \psi = \nabla\cdot({v}^*- {\cal{P}} {v}^*) \ \;\text{ in }\; \varOmega, \quad {n}\cdot\nabla\psi = 0 \;\text{ on } \partial\varOmega $

满足定理1中的可解条件, 从而可唯一确定$ \nabla\psi $. 但是通常

$\nabla\psi\ne {v}^*- {\mathscr{P}} {v}^* = \nabla \phi. $

例如$ \psi-\phi $可以是2.3节中提及的Stokes压强.
因此, (48)式中的${\cal{P}} {v}^* = {v}^*-\nabla \psi$并不是对$ {v}^* $的唯一分解, 而是由给定的$ {\cal{P}} $直接决定了梯度场$ \nabla \psi $. 由于许多向量场无法写成梯度场的形式, 因此(48)式和(49)式定义的广义投影算子是一种特殊的线性算子, 其特殊性在于总是为输入的$ {v}^* $选择一个输出${\cal{P}} {v}^*$, 使得两者之差为一个梯度场.
为什么可以用广义投影算子来刻画那些不满足(47)式的离散投影呢? 这是因为离散投影通常是Leray-Helmholtz投影的一个高阶($p\geqslant 4$)近似, 从而有

${v}^* = {\mathscr{P}} {v}^* + \nabla \phi = {{P}} {v}^* + O(h^p) + \nabla \phi.$

和$ \nabla \phi $相比, 高阶小量$ O(h^p) $往往可以忽略.
2 -->

4.2.Laplace算子和广义投影算子的交换子$[{{\nabla}}, {{{\cal{P}}}}]$

首先定义梯度-散度-投影交换子为

$ [\nabla\nabla\cdot, {\cal{P}}]: = \nabla\nabla\cdot {\cal{P}} - {\cal{P}}\nabla\nabla\cdot. $

由(48)式、$ \Delta $和$ \nabla $的交换性以及(13)式可得

$\begin{split}\;& \Delta {\cal{P}} {v}^* = \Delta({v}^* - \nabla\phi) = \Delta {v}^* - \nabla\nabla\cdot\nabla\phi \\=\;& \Delta {v}^* - \nabla\nabla\cdot ({v}^*-{v}) = {\cal B} {v}^* +\nabla \nabla\cdot {\cal{P}}{v}^*. \end{split} $

因此

$\begin{array}{l} \Delta {\cal{P}} = {\cal B} + \nabla \nabla\cdot {\cal{P}}. \end{array} $

对(52)式中第三项和最后一项应用$ {\cal{P}} $, 有

$ {\cal{P}} \Delta - {\cal{P}} \nabla\nabla\cdot ({\cal I}- {\cal{P}}) = {\cal{P}} {\cal B} + {\cal{P}} \nabla\nabla\cdot {\cal{P}}, $

整理后得到

$ {\cal{P}} \Delta = {\cal{P}} {\cal B} + {\cal{P}} \nabla\nabla\cdot. $

由(53)式, (54)式和 (51)式可得 Laplace算子和广义投影算子的交换子为

$ [\Delta, {\cal{P}}]: = \Delta {\cal{P}} - {\cal{P}} \Delta = ({\cal I}- {\cal{P}}) {\cal B} + [\nabla\nabla\cdot, {\cal{P}}]. $

若$ {\cal{P}} $满足(10)式, 则$ [\nabla\nabla\cdot, {\cal{P}}]\equiv {0} $, 此时(55)式简化为(17)式.
2 -->

4.3.INSE的GePUP表述

由4.1节中的讨论可知, 当广义投影算子$ {\cal{P}} $作用在一个散度为零的流场$ {u} $上时, 可以看成是对$ {u} $加一个梯度扰动项使其散度不一定为零. 在文献[5]中, 我们对动量方程(1a)应用广义投影算子$ {\cal{P}} $并结合 (55)式和(51)式推导了 INSE的GePUP表述. 本文的推导过程在形式上略有不同.
由(6)式, (53)式, (5)式以及(15)式可得

$ {\cal{P}}\frac{\partial{u}}{\partial t} = {\cal{P}}\left( {a}^* - \nabla p \right) -\nu{\cal B}{u} + \nu\Delta {\cal{P}}{u} - \nu\nabla\nabla\cdot {\cal{P}}{u} $

$ \Rightarrow\ \frac{\partial {\cal{P}}{u}}{\partial t} = {g} -{u}\cdot\nabla{u} -\nabla q + \nu\Delta {\cal{P}}{u} -\nu \nabla\nabla\cdot {\cal{P}}{u}, $

其中梯度场$ \nabla q $定义为

$ \nabla q: = {a}^* - {\cal{P}} {a} - \nu\nabla\nabla\cdot{u}. $

由(6)式及广义投影算子的定义可知 (58)式的右端确实是个梯度场. 标量q可由如下Neumann边值问题确定:

$ \Delta q = \nabla\cdot \left( {g} -{u}\cdot\nabla{u} \right) - \nabla\cdot {\cal{P}} {a}\ \; \text{ in }\; \varOmega,\tag{59a} $

$ {n}\cdot\nabla q = {n}\cdot ({a}^* - \nu\nabla\nabla\cdot{u}) -{n}\cdot {\cal{P}} {a}\ \;\text{ on } \;\partial\varOmega, \tag{59b}$

其中(59a)式是对(58)式取散度后由$ \nabla\cdot $和$ \Delta $的交换性得到的, 而(59b)式是 (58)式的法向分量. 在无穿透边界条件下, 由(5)式和(49)式可得

$\begin{split} {n}\cdot {\cal{P}}{a} =\;& {n}\cdot {\cal{P}}\left(\frac{\partial{u}}{\partial t}\right) = {n}\cdot\frac{\partial {\cal{P}}{u}}{\partial t}\\=\;& {n}\cdot \frac{\partial{u}}{\partial t} = 0 \;\text{ on }\; \partial\varOmega. \end{split} $

另外, 无滑移条件(2)式意味着对流项在边界上为零:

$ {u}\cdot\nabla{u} = {0} \;\text{ on }\; \partial\varOmega. $

再结合(5)式知可用$ {g}+\nu\Delta {u} $替换 (59b)式中的$ {a}^* $.
定义投影流速

$\begin{array}{l} {w}: = {\cal{P}} {u}, \end{array} $

并假设其散度满足

$ \nabla\nabla\cdot {w} = {0},\tag{63a} $

$ \frac{\partial \nabla\cdot{w}}{\partial t} = 0,\tag{63b} $

可以得到无滑移边界条件下INSE的GePUP表述:

$ \frac{\partial{ w} }{\partial t} = { g } -{ u}\cdot\nabla{ u} -\nabla q + \nu\Delta{ w}\quad {\rm{in }}\; \varOmega,{\tag{64a}} $

$ { w} = {0}\quad {\rm{on}}\ \; \partial \varOmega,\tag{64b} $

$ { u } = {\mathscr{P}} { w}\quad {\rm{in }}\ \;\varOmega, \tag{64c}$

$ {{u}}\cdot {{n}} = 0\quad {\rm{on }}\ \;\partial \varOmega ,\tag{64d} $

$ \Delta q = \nabla\cdot \left( {{g}} -{{u}}\cdot\nabla{{u}} \right) \quad{\rm{in }}\ \;\varOmega, \tag{64e}$

$ {n}\cdot\nabla q = {n}\cdot \left( {g} + \nu\Delta{u}-\nu\nabla\nabla\cdot{u} \right)\quad {\rm{on }}\ \;\partial \varOmega, \tag{64f}$

其中(64a)式是由 (57)式和假定(63a)式得到的, (64e)式是由(59a)式和假定(63b)式得到的, 边界条件(64f)式来自于(59b)式, (60)式以及(61)式. 由于广义投影算子对不可压流场$ {u} $的扰动是增加了一个梯度项, 而任何梯度都在Leray-Helmholtz投影算子的零空间中, 因此(64c)式显然成立.
在(64)式这个三变量系统中, 投影流速$ {w} $是唯一演化变量, 无源流速$ {u} $由$ {w} $唯一确定, 而压强梯度$ \nabla q $又由$ {u} $唯一确定. 因此(64)式可以看成是关于$ {w} $的动力过程, 对其进行空间MOL离散之后得到一个常微分方程组, 而求解这个常微分方程组的时间积分方法可作为“黑匣子”使用, 这样就可以实现第1节中的(Goal-3).
上述GePUP表述保持了UPPE表述的核心优势, 即投影流速的散度由一个热方程所控制. 确实, 对(64a)式取散度并结合(64e)式可得

$ \frac{\partial (\nabla\cdot{w})}{\partial t} = \nu\Delta (\nabla\cdot{w}), $

若在$ \partial\varOmega $上指定边界条件$ \nabla\cdot{w} = 0 $, 则由热方程的极大值原理可知$ \nabla\cdot{w} $在$ \varOmega $上随时间呈指数衰减^①.
GePUP表述相对UPPE表述的优势有以下几点:
(I) 推导过程不依赖于任何Leray-Helmholtz投影算子${\mathscr{P}}$的性质(10)式;
(II) 演化变量$ {w} $没有不可压限制条件;
(III) 投影算子${\mathscr{P}}$出现在(64)式的右端, 因此用离散投影算子来近似${\mathscr{P}}$时对数值精度和算法稳定性的影响非常清晰.
对一个满足不可压条件的$ {w} $的初始值来说, 如果INSE的解存在且唯一, 热方程(65)将使得假定(63)式成立. 但如果INSE的解本身将爆破的话, 其散度必先爆破. 因此只有当INSE在理论上确实有解的时候我们的计算才有意义, 换言之, 我们无法通过重新表述INSE的方式回避其解存在性的千禧年问题. 这是假定(63)式背后的主要思想.
2 -->

4.4.矩形区域上的空间离散

用正方形网格划分计算域$ \varOmega $, 并以多元索引$ {{i}}\in{\mathbb{Z}}^ { {{D}}} $代表第$ {{i}} $个控制体$ {\mathcal C}_{ {{i}}} $, 其区域为

$ {\mathcal C}_{ {{i}}}: = \bigl[ {x}_O + {{i}} h,\ {x}_O + \left( {{i}}+ {\mathds{1}}\right)h \bigr], $

控制体$ {{i}} $在维度d上的高侧界面表示为

$ {\mathcal F}_{ {{i}}+\frac{{1}}{2} {{e}}^d}: = \left[ {x}_O + \left( {{i}} + {{e}}^d\right)h, \ {x}_O + \left( {{i}}+ {\mathds{1}}\right)h \right], $

其中$ {x}_O\in {\mathbb{R}}^{ { {{D}}}} $是固定的坐标原点, h是均匀网格间距, ${\mathds{1}}\in {\mathbb{Z}}^{ { {{D}}}}$是分量全为1的多元索引, $ {{e}}^d\in {\mathbb{Z}}^{ { {{D}}}} $表示第d个分量为1, 其余分量全为0的多元索引.
如图1所示, $ \phi $在控制体$ {\mathcal C}_{ {{i}}} $上的平均值定义为
图 1 面平均值和控制体平均值. 面平均值用带$ \left\langle {\cdot} \right\rangle $且下标为分数的符号表示, 而控制体平均值对应下标为整数. 因此 $ \left\langle {\phi} \right\rangle_{ {{i}}+\frac{1}{2} {{e}}^0} $, $ \left\langle {\phi} \right\rangle_{ {{i}}- {{e}}^0+\frac{1}{2} {{e}}^1} $和$ \left\langle {u_0} \right\rangle_{ {{i}}+\frac{3}{2} {{e}}^0 + {{e}}^1} $表示面平均值, 而$ \left\langle {\phi} \right\rangle_{ {{i}}} $是控制体平均值^[5]
Figure1. Notation of face-averaged and cell-averaged values. A symbol with angled brackets denotes either a cell-averaged value if the subscript is an integer multi-index, or a face-averaged value if the subscript is a fractional multi-index. Hence $ \left\langle {\phi} \right\rangle_{ {{i}}+\frac{1}{2} {{e}}^0} $, $ \left\langle {\phi} \right\rangle_{ {{i}}- {{e}}^0+\frac{1}{2} {{e}}^1} $, and $ \left\langle {u_0} \right\rangle_{ {{i}}+\frac{3}{2} {{e}}^0 + {{e}}^1} $ are face-averaged values while $ \left\langle {\phi} \right\rangle_{ {{i}}} $ is a cell-averaged value. Horizontal and vertical hatches represent the averaging processes over a vertical cell face and a horizontal cell face, respectively. Light gray area represents averaging over a cell^[5].

$ \langle {\phi} \rangle_{ {{i}}}: = \frac{1}{h^{ { {{D}}}}}\int_{{\mathcal C}_{ {{i}}}} \phi\left({x}\right) {\rm{d}} {{x}}, $

而$ \phi $在控制体边界面$ {\mathcal F}_{ {{i}}+\tfrac{{1}}{2} {{e}}^d} $上的平均值定义为

$ \left\langle {\phi} \right\rangle_{ {{i}}+\frac{{1}}{2} {{e}}^d}: = \frac{1}{h^{ { {{D}}}-1}} \int_{{\mathcal F}_{ {{i}}+\frac{{1}}{2} {{e}}^d}} \phi\left({x}\right) {\rm{d}} {{x}}. $

这两种平均值之间的关系可通过微积分基本定理得到^[25]:

$ \begin{split} \left\langle {\phi} \right\rangle_{ {{i}}+\tfrac{{1}}{2} {{e}}^d} =\;& \frac{1}{12}\Big( - \left\langle {\phi} \right\rangle_{ {{i}}+2 {{e}}^d} + 7 \left\langle {\phi} \right\rangle_{ {{i}}+ {{e}}^{d}} \\ &+ 7 \langle {\phi} \rangle_{ {{i}}} \!-\! \left\langle {\phi} \right\rangle_{ {{i}}- {{e}}^d} \Big) \!+\! O(h^4), \end{split}\tag{70a} $

$ \begin{split} \left\langle {\frac{\partial \phi}{\partial x_d}} \right\rangle_{ {{i}}+\frac{{1}}{2} {{e}}^d} = \;&\frac{1}{12 h}\Big( - \left\langle {\phi} \right\rangle_{ {{i}}+2 {{e}}^d} \!+\! 15 \left\langle {\phi} \right\rangle_{ {{i}}+ {{e}}^{d}} \\ & \!\!\! - 15 \langle {\phi} \rangle_{ {{i}}} \!+\! \left\langle {\phi} \right\rangle_{ {{i}}- {{e}}^d} \Big) \!+\! O(h^4).\;\;\;\;\;\; \end{split}\tag{70b} $

由于控制体中心的点值和控制体平均值之间相差$ O(h^2) $, 二阶有限体积方法不需要对它们进行区分, 而在四阶方法中则必须考虑这一差异. 此外, 作用于平均值上的有限体积算子与作用于点值上的有限差分算子是不同的.
2 -->

4.5.四阶离散算子

离散梯度、离散散度以及离散Laplace算子分别定义如下:

$ \begin{split} {{G}}_d \left\langle {\phi} \right\rangle_{ {{i}}}: =\;& \frac{1}{12 h} \Big( - \left\langle {\phi} \right\rangle_{ {{i}}+2 {{e}}^d} + 8 \left\langle {\phi} \right\rangle_{ {{i}}+ {{e}}^d} \\ &- 8 \left\langle {\phi} \right\rangle_{ {{i}}- {{e}}^d} + \left\langle {\phi} \right\rangle_{ {{i}}-2 {{e}}^d} \Big), \end{split} $

$ \begin{split} {{D}} \left\langle {{u}} \right\rangle_{ {{i}}}: =\;& \frac{1}{12 h} \sum_d\Big( - \left\langle {u_d} \right\rangle_{ {{i}}+2 {{e}}^d} + 8 \left\langle {u_d} \right\rangle_{ {{i}}+ {{e}}^d} \\ &- 8 \left\langle {u_d} \right\rangle_{ {{i}}- {{e}}^d} + \left\langle {u_d} \right\rangle_{ {{i}}-2 {{e}}^d} \Big), \end{split} $

$ \begin{split} {{L}} \left\langle {\phi} \right\rangle_{ {{i}}}: =\;& \frac{1}{12h^2} \sum_{d} \Big( -\left\langle {\phi} \right\rangle_{ {{i}}+ 2 {{e}}^{d}} + 16 \left\langle {\phi} \right\rangle_{ {{i}}+ {{e}}^{d}} \\ &- 30 \left\langle {\phi} \right\rangle_{ {{i}}} + 16 \left\langle {\phi} \right\rangle_{ {{i}}- {{e}}^{d}} - \left\langle {\phi} \right\rangle_{ {{i}} - 2 {{e}}^{d}} \Big). \end{split} $

离散散度算子$ {{D}} $还可作用在张量平均值上

$ {{D}} \left\langle {{u}{u}} \right\rangle_{ {{i}}}: \!=\! \frac{1}{h} \sum_d \left( {F}\left\langle {u_d,{u}} \right\rangle_{ {{{i}}+\tfrac{1}{2}{{e}}^d}} -{F}\left\langle {u_d,{u}} \right\rangle_{ {{{i}}-\tfrac{1}{2}{{e}}^d}} \right), $

其中两标量乘积的面平均值可近似为

$ \begin{split} &{F}\left\langle {\varphi,\psi} \right\rangle_{ {{i}}+\tfrac{{1}}{2} {{e}}^d}: =\left\langle {\varphi} \right\rangle_{ {{i}}+\tfrac{{1}}{2} {{e}}^d} \left\langle {\psi} \right\rangle_{ {{i}}+\tfrac{{1}}{2} {{e}}^d} \\ & \qquad +\frac{h^2}{12} \sum \limits_{d'\ne d} \left( {{G}}_{d'}^{\perp}\varphi \right)_{ {{i}}+\tfrac{{1}}{2} {{e}}^d} \left( {{G}}_{d'}^{\perp}\psi \right)_{ {{i}}+\tfrac{{1}}{2} {{e}}^d}, \end{split} $

而$ {{G}}_{d'}^{\perp} $表示横向离散梯度算子

$ \left( {{G}}_{d'}^{\perp} \varphi \right)_{ {{i}}+\tfrac{{1}}{2} {{e}}^d}: = \frac{1}{2h}\left( \left\langle {\varphi} \right\rangle_{ {{i}}+\tfrac{{1}}{2} {{e}}^d+ {{e}}^{d'}} - \left\langle {\varphi} \right\rangle_{ {{i}}+\tfrac{{1}}{2} {{e}}^d- {{e}}^{d'}} \right). $

命题4　(71)式—(75)式中的离散算子都具有四阶精度:

$ {{G}}_d\left\langle {\phi} \right\rangle_{ {{i}}} = \frac{1}{h^{ { {{D}}}}}\int_{\mathcal C_{{i}}}\frac{\partial\phi}{\partial x_d} {\rm{d}}{x} + O(h^4),\tag{77a} $

$ {{D}}\left\langle {{u}} \right\rangle_{{i}} = \frac{1}{h^{ { {{D}}}}}\int_{\mathcal C_{{i}}}\nabla\cdot{u} {\rm{d}}{x} + O(h^4),\tag{77b} $

$ {{L}}\left\langle {\phi} \right\rangle_{{i}} = \frac{1}{h^{ { {{D}}}}}\int_{\mathcal C_{{i}}}\Delta\phi {\rm{d}}{x}+O(h^4),\tag{77c} $

$\begin{split}& {{D}}\left\langle {{u}} \right\rangle_{{i}} = O(h^4) \Rightarrow {{D}}\left\langle {{u}{u}} \right\rangle_{{i}} \\=& \frac{1}{h^{ { {{D}}}}}\int_{\mathcal C_{{i}}}({u}\cdot\nabla){u} {\rm{d}}{x} + O(h^4),\end{split}\tag{77d} $

$ {F}\left\langle {\varphi, \psi} \right\rangle_{ {{i}}+\frac{{1}}{2} {{e}}^d} = \frac{1}{h^{ { {{D}}}-1}}\int_{\mathcal F_{ {{i}}+\frac{{1}}{2} {{e}}^d}}\varphi\psi {\rm{d}}{x} + O(h^4). \tag{77e} $

证明见文献[4].
如图2所示, 我们在计算域边界之外再设置两层虚拟单元. 如果是周期性边界条件则直接把计算域内对应控制体的值复制到虚拟单元中, 否则由给定边界条件外插得到虚拟单元的平均值. 例如, 对于Dirichlet边界条件$ \psi = 0 $, 用以下五阶公式来填充虚拟单元$ {{i}}+ {{e}}^d $及$ {{i}}+2 {{e}}^d $:
图 2 计算域边界上的两层虚拟单元$ (d = 1) $. 粗实线代表域边界, 细实线表示内部控制体, 而虚线表示虚拟单元^[5]
Figure2. An example ($ d = 1 $) of ghost filling for the domain boundary. The thick line represents a physical boundary, the thin solid lines interior cells, and the dotted lines ghost cells^[5].

$ \begin{split} \left\langle {\psi} \right\rangle_{ {{i}}+ {{e}}^d} =\;& \frac{1}{12} ( -77 \left\langle {\psi} \right\rangle_{ {{i}}} + 43 \left\langle {\psi} \right\rangle_{ {{i}}- {{e}}^d}\\ &-17 \left\langle {\psi} \right\rangle_{ {{i}}-2 {{e}}^d} + 3 \left\langle {\psi} \right\rangle_{ {{i}}-3 {{e}}^d} ) + O(h^5), \\ \left\langle {\psi} \right\rangle_{ {{i}}+ 2 {{e}}^d} =\;& \frac{1}{12} ( -505 \left\langle {\psi} \right\rangle_{ {{i}}} + 335 \left\langle {\psi} \right\rangle_{ {{i}}- {{e}}^d} \\ &-145 \left\langle {\psi} \right\rangle_{ {{i}}-2 {{e}}^d} + 27 \left\langle {\psi} \right\rangle_{ {{i}}-3 {{e}}^d}) + O(h^5). \end{split} $

在确定了所有虚拟单元上的平均值之后, 再用格式(71)式—(74)式计算离散算子. 这样做的好处是将边界条件和连续算子的离散形式完全解耦.
用(71)式—(73)式可定义离散投影算子

$\begin{array}{l} {{P}}: = {{I}} - {{G}} {{L}}^{-1} {{D}}. \end{array} $

由$ {{L}}\ne {{D}} {{G}} $可知$ {{P}}^2\ne {{P}} $. 在周期边界条件下, $ {{P}} $是对Leray-Helmholtz投影算子${\mathscr{P}}$的四阶近似, 并且可以证明其谱半径和2-范数都是1^[4]. 在无滑移边界下, $ {{P}} $不是对称矩阵, 但是数值实验验证了$ {{P}} $的谱半径不大于1^[5]. 总之, $ {{P}} $可以看成是有限维线性空间上的广义投影算子, 作为对GePUP表述(64)式中Leray-Helmholtz投影算子的近似.
2 -->

4.6.不规则区域上的空间离散

如果计算域$ \varOmega $是不规则的, 我们仍然使用正方形控制体离散一个包含了$ \varOmega $的矩形区域, 并定义控制体流相$ {\cal V}_{ {{i}}}: = {\cal C}_{ {{i}}}\cap\varOmega $. 当$ {\cal V}_{ {{i}}} = {\cal C}_{ {{i}}}, \emptyset $时, 分别称$ {\cal C}_{ {{i}}} $为纯流相和无流相控制体, 这两种以外的情况称为边界控制体. 记$ \|{\cal V}\| $为控制体流相$ {\cal V} $的体积, 定义$ \phi $在$ {\cal V} $上的平均值为

$\begin{array}{l} \left\langle {\phi} \right\rangle_{{\cal V}}: = \frac{1}{\|{\cal V}\|} \int_{{\cal V}}\phi({x}) {\rm{d}} {x}. \end{array} $

对于远离边界的纯流相控制体, 空间算子$ {\cal L} $在其上的平均值$ \left\langle {{\cal L} \phi} \right\rangle_{ {{i}}} $可以用4.5节中的离散算子近似到四阶精度. 在图3(a)的例子中, $ {\cal L} = \varDelta $. 对于边界控制体$ {\cal C}_{ {{i}}} $或其附近的纯流相控制体, 我们也希望用一些控制体平均值的线性组合来近似$ \left\langle {{\cal L}\phi} \right\rangle_{{\cal C}_{ {{i}}}} $. 如图3(b)所示, 主要思路是用选定控制体的平均值以及边界条件拟合一个多维完全多项式. 其核心难点是如何选择$ {\cal V}_{ {{i}}} $附近的一组控制体流相使得多项式拟合的方程组有唯一解, 这个问题最终归结为在$ {\mathbb{Z}}^{ { {{D}}}} $中找到 $ {n+ { {{D}}} \choose { {{D}}}} $个点满足一些适定条件. 我们提出了一个基于对称群作用的回溯算法^[27], 见图4.
图 3 在不规则计算域上用正交网格离散Laplace算子$ \varDelta $^[26]
Figure3. Discretizing the Laplacian on irregular domains with rectangular grids^[26].

图 4 从给定起始点开始选择一组格点以拟合多维n次完全多项式. 用拟合的完全多项式可以离散空间算子, 并且能很方便地满足边界条件^[26]
Figure4. For a given starting point, we choose a set of nodes to fit a high-order, multi-variate, complete polynomial, which facilitates the discretization of spatial operators and the enforcement of boundary conditions^[26].

2 -->

4.7.时间积分

对(64)式在控制体上积分, 并用离散算子近似连续算子的积分, 可得常微分方程组:

$ \frac{ \mathrm{d} \left\langle {{w}} \right\rangle}{ \mathrm{d} t} = \left\langle {{g}} \right\rangle- {{D}}\left\langle {{u}{u}} \right\rangle - {{G}}\left\langle {q} \right\rangle + \nu {{L}} \left\langle {{w}} \right\rangle, \tag{81a}$

$ {{L}} \left\langle {q} \right\rangle = {{D}}\bigl( \left\langle {{g}} \right\rangle - {{D}}\left\langle {{u}{u}} \right\rangle \bigr), \tag{81b}$

$ \left\langle {{u}} \right\rangle = {{P}} \left\langle {{w}} \right\rangle,\tag{81c} $

其中外力项$ \left\langle {{g}} \right\rangle $已知, 未知量是控制体平均值向量$ \left\langle {{w}} \right\rangle $, 而$ \left\langle {{u}} \right\rangle $和$ \left\langle {q} \right\rangle $都可以通过$ \left\langle {{w}} \right\rangle $得到. (81)式的初始条件为$ \left\langle {{w}} \right\rangle(t_0) = \left\langle {{u}} \right\rangle(t_0) $, 边界条件可由(64)式得到.
由于$ {{P}} $是对Leray-Helmholtz投影算子$ {\mathscr{P}} $的四阶近似, 结合命题4可知常微分方程组(81)是对 GePUP表述(64)式的四阶离散.
在求解半离散系统(81)式时, 时间积分方法可以完全作为黑匣子使用. 在文献[5]中我们把GePUP和IMEX方法^[28]耦合起来得到了一组求解INSE的半隐式时间推进算法, 称为GePUP-IMEX. 这里以Kennedy与Carpenter^[29] 提出的ERK-ESDIRK方法作为一个例子来阐述GePUP-IMEX算法.
对于

$ \frac{ \mathrm{d} \psi}{ \mathrm{d} t} = {{X}}^{[\text{E}]}(\psi, t) + {{X}}^{[\text{I}]}(\psi), $

ERK-ESDIRK方法在每个时间步的求解步骤如下:

$ \psi^{(1)} = \psi^n \approx \psi(t^n), \tag{83a} $

$ \begin{split} &\text{for } s = 2,\;3,\;\cdots, \;n_{ \textsf{s}}, \\ &\bigl(I - {\Delta t}\gamma {{X}}^{[{\rm {I}}]} \bigr) \psi^{(\rm s)} = \psi^{n} + {\Delta t} \sum\limits^{s-1}_{j = 1} a^{[\rm E]}_{s,j} {{X}}^{[\rm{E}]} \left(\psi^{(j)}, t^{(j)}\right) \\ &\qquad \qquad \qquad\qquad + {\Delta t} \sum\limits^{s-1}_{j = 1} a^{[\rm I]}_{s,j} {{X}}^{[\rm{I}]}\psi^{(j)}, \end{split} \tag{83b} $

$ \psi^{n+1} = \psi^{( n_{ \textsf{s}})} + {\Delta t} \sum\limits^{ n_{\textsf{s}}}_{j = 1} \left(b_j - a^{[\rm E]}_{ n_{ \textsf{s}},j}\right) {{X}}^{[\rm{E}]} \left(\psi^{(j)}, t^{(j)}\right), \tag{83c} $

其中上标$ ^{(\rm s)} $代表中间阶段, $ t^{(\rm s)} = t^n + c_{s} {\Delta t} $是该阶段对应时刻; $ n_{ {s}} $是总阶段数; $ A, {b}, {c} $是Butcher表中的系数. $ A^{[\rm E]}, {b} $和$ {c} $表示一个ERK方法, 而$ A^{[\rm I]}, {b} $和$ {c} $表示一个ESDIRK方法. 在每个中间阶段, $ \psi^{(\rm s)} $通过求解线性方程组(83b)得到, 其中右端项是结合前面步骤值计算得到的, 故(83)式确实属于半隐式格式. 另外(83c)式是由 $b^{[\rm E]}_{j} = b^{[I]}_{j} = b_j = a^{[\rm I]}_{ n_{ \textsf{s}}, j}$得到的. 有关该方法的更多细节请参考文献[29].
对GePUP半离散系统(81)定义

$ {{X}}^{[\rm{E}]}: = \left\langle {{g}} \right\rangle- {{D}}\left\langle {{u}{u}} \right\rangle- {{G}}\left\langle {q} \right\rangle; \quad {{X}}^{[\rm{I}]}: = \nu {{L}}\left\langle {{w}} \right\rangle, $

并套用ERK-ESDIRK方法(83)式的契约, 我们得到如下的GePUP-IMEX算法:

$ \left\langle {{w}} \right\rangle^{(1)} = \left\langle {{w}} \right\rangle^n,\tag{85a} $

$ \left\{ \begin{aligned} &~~\text{for } s = 2,\; 3, \;\cdots,\; n_{ \textsf{s}}, \\ &\bigl( {{I}} - {\Delta t}\nu\gamma {{L}}\bigr) \left\langle {{w}} \right\rangle^{(s)} \\ =& \left\langle {{w}} \right\rangle^n + {\Delta t} \sum^{s-1}_{j = 1} a^{[\rm{E}]}_{s,j} {{X}}^{[\rm{E}]} \left(\left\langle {{u}} \right\rangle^{(j)}, t^{(j)}\right) \\ &~~ + {\Delta t}\nu \sum^{s-1}_{j = 1} a^{[\rm{I}]}_{s,j} {{L}}\left\langle {{w}} \right\rangle^{(j)}, \\ &~~\left\langle {{u}} \right\rangle^{(s)} = {{P}} \left\langle {{w}} \right\rangle^{(s)}, \end{aligned} \right.\tag{85b} $

$ \left\{ \begin{aligned} &\left\langle {{w}} \right\rangle^{*} = \left\langle {{w}} \right\rangle^{( n_{ \textsf{s}})} + {\Delta t} \sum^{ n_{\textsf{s}}}_{j = 1} \left(b_j - a^{[\rm E]}_{ n_{ \textsf{s}},j}\right) \\ & ~~\qquad\times {{X}}^{[\rm{E}]} \left(\left\langle {{u}} \right\rangle^{(j)}, t^{(j)}\right), \\ &\left\langle {{u}} \right\rangle^{n+1} = {{P}} \left\langle {{w}} \right\rangle^{*}, \\ &\left\langle {{w}} \right\rangle^{n+1} = \left\langle {{u}} \right\rangle^{n+1}, \end{aligned} \right. \tag{85c} $

在初始时刻$ t = t_0 $, 有$ \left\langle {{w}} \right\rangle^0 = \left\langle {{u}} \right\rangle^0 \approx \left\langle {{u}(t_0)} \right\rangle $.
GePUP-IMEX算法(85)式的具体步骤直接来自于将ERK-ESDIRK方法应用到GePUP的半离散系统(81)式上. 其中(85b)式的第一步对应于(81a)式和(81b)式, 而(85b)式的第二步对应于(81c)式, 即用离散投影$ {{P}} $将$ \left\langle {{w}} \right\rangle^{(s)} $映成$ \left\langle {{u}} \right\rangle^{(s)} $, 以便为下一阶段计算${{X}}^{[\rm {E}]}(\left\langle {{u}} \right\rangle^{(j)}, $$ t^{(j)})$做好准备. 在(85c)式中, 算子$ {{X}}^{[\text{E}]} $中的非线性对流项导致$ \left\langle {{w}} \right\rangle^* $的散度要远大于中间阶段$ \left\langle {{w}} \right\rangle^{(s)} $的散度, 因此我们将该时间步的最终解$ \left\langle {{w}} \right\rangle^{n+1} $设置为$ \left\langle {{u}} \right\rangle^{n+1} $.
根据(85b)式和(81)式, 在每个中间阶段, 需要求解三个线性方程组.
(HLS-1) Neumann边界条件下的Poisson方程(81b): 用来从$\left\langle {{u}} \right\rangle$中提取$ \left\langle {q} \right\rangle $;
(HLS-2)无滑移边界条件下的Helmholtz方程(85b): 用来在时间上推进$ \left\langle {{w}} \right\rangle $;
(HLS-3) Neumann边界条件下的Poisson方程(81c): 用来对$\left\langle {{w}} \right\rangle^*$进行投影.
求解这些方程组的具体算法见文献[5]. 在(85c)式中, 需要求解(HLS-1)和(HLS-3), 因此在每个时间步内, GePUP-IMEX算法($ n_{ \textsf{s}} = 6 $)总共需要求解17个线性方程组. 我们使用的方法是标准的多重网格V循环, 其中光滑迭代为加权Jacobi方法($ \omega = {2}/{3} $), 网格转移算子为单射算子和常值限制算子^[30]. 由于$ - {{L}} $属于本性正型^{[31, 32]}, 具有M-矩阵的良好性质, 例如正定性等. (HLS-2)中的算子$ \bigl( {{I}}- {\Delta t}\nu\gamma {{L}}\bigr) $的特征值的实部全为正, 因此经典几何多重网格方法对三个线性方程组的求解都相当高效.
为了估算稳定的库朗数范围, 我们将对流项线性化, 再用与处理对流扩散方程^[25]同样的方法可估算稳定的库朗数范围为

$ {{\textsf{Cr}}} \leqslant {2.91}/{ { {{D}}}}. $

综上所述, GePUP-IMEX算法是将IMEX时间积分方法直接应用到GePUP半离散系统上得到的. 这个简单的思路对于其他时间积分方法如显式Runge-Kutta方法也是适用的, 具体步骤往往会更简单. 例如, 我们将经典四阶Runge-Kutta方法应用到GePUP半离散系统上得到的算法称为GePUP-ERK算法^[5].

2

-->

5.1.不规则区域上求解INSE

在不规则区域(倾斜正方形)

$\begin{split} { {\varOmega}} =\;& \left\{ \begin{bmatrix} \cos\dfrac{\pi}{6} & -\sin\dfrac{\pi}{6} \\ \sin\dfrac{\pi}{6} & \cos\dfrac{\pi}{6} \end{bmatrix} \begin{pmatrix} \xi \\ \eta \end{pmatrix} \right.\\&\left.+ \begin{pmatrix} \dfrac{1}{2} \\ 0 \end{pmatrix}: (\xi, \eta)\in [0,1]^2 \right\} \end{split} $

上给定流速和压强的精确解

$ {u}(x,y,t) = \left( \begin{array}{c} \sin^2(\pi x) \sin(2 \pi y) \\ -\sin(2 \pi x) \sin^2(\pi y) \end{array} \right)\pi \cos{t}, $

$ p(x,y,t) = -\cos(\pi x) \sin(\pi y) \cos {t}, $

可以根据动量方程(1a)推导出本节测试中使用的外力项${g}$.
在计算域(87)上, 用(88)式设置流速场的初始条件和无滑移边界条件, 用GePUP-ERK算法将INSE的解从$ t_0 = 0 $推进到$ t_{\rm e} = 0.1 $. 表1列出了最终时刻结果的误差和收敛阶, 显然流速$ {u} $和压强p都能达到四阶收敛.

h	${1}/{32}$	Rate	${1}/{64}$	Rate	$1/{128}$	Rate	$1/{256}$
${u}$ $ L_{\infty} $	2.63 × 10–4	3.94	1.71 × 10–5	3.98	1.08 × 10–6	4.00	6.78 × 10–8
$ {u} $ $ L_1 $	1.63 × 10–4	3.87	1.12 × 10–5	3.97	7.14 × 10–7	3.99	4.48 × 10–8
$ {u} $ $ L_2 $	1.49 × 10–4	3.89	1.01 × 10–5	3.98	6.38 × 10–7	4.00	4.00 × 10–8
p $ L_{\infty} $	1.63 × 10–3	3.99	1.02 × 10–4	4.00	6.39 × 10–6	4.01	3.97 × 10–7
p $ L_1 $	3.26 × 10–4	4.00	2.04 × 10–5	4.03	1.25 × 10–6	3.97	7.97 × 10–8
p $ L_2 $	4.60 × 10–4	3.99	2.89 × 10–5	4.01	1.80 × 10–6	3.99	1.13 × 10–7

表1用GePUP-ERK算法在不规则计算域(87)上求解以(88)式和(89)式为解析解的INSE得到的误差和收敛阶. Re = 1000, $ t_0 = 0 $, $ t_{\rm e} = 0.1 $, $ {{\textsf{Cr}}} = 0.2 $^[5]
Table1.Errors and convergence rates of GePUP-ERK for solving the INSE on the irregular domain (87) with Eq. (88) and Eq. (89) as the analytic solutions. Re = 1000, $ t_0 = 0 $, $ t_{\rm e} = 0.1 $, $ {{\textsf{Cr}}} = 0.2 $^[5].

2 -->

5.2.三维顶盖驱动方腔流

三维顶盖驱动方腔内的不可压缩流动几乎表现出流体力学中所有常见的现象: 漩涡、二次流、混沌粒子运动、不稳定性、过渡、湍流和如展向涡的复杂三维流动结构^[33], 因此其数值模拟具有重要的科学意义. 另外由于几何设定非常简单, 顶盖驱动方腔流是检验INSE数值求解方法最常使用的基准测试之一.
在规则计算区域$ [0, 1]\times[0, 1]\times[0, 2] $上, 用GePUP-IMEX算法(85)式对三维顶盖驱动方腔流进行数值模拟. 选择的具体参数为: 流速场$ {u} = (u, v, w) $的边界条件为在$ x = 1 $处, $ v = 1 $, 其余边界上$ {u} = {0} $; 而在初始时刻, $ {u} $在计算域内全为0, $ t_0 = 0 $及$ t_{\rm e} = 12 $. 雷诺数Re = 1000, 库朗数$ {{\textsf{Cr}}} = 0.5 $, 均匀网格间距$h = {1}/{128}$. 下面的数值结果可以在文献[5]中找到.
图5给出了在时刻$ t = 2, 4, 6, 8, 10, 12 $, 对称平面$ z = 1 $内流速场数值解的流线. 数值模拟刚开始时, 盖子的移动会沿腔盖产生很大的剪切应力, 并在$ t = 2 $之前形成一个漩涡: 位于漩涡中心控制体的流速要比附近所有控制体流速小50倍左右. 随着模拟继续进行, 逆时针旋转漩涡附近的高速流体逐渐穿透最初静止的流体, 在$ t = 4 $之后不久, 又有另外一个顺时针旋转的漩涡开始形成, 并在$ t = 6 $变得显著. 然后二次涡开始增强并移动到计算域左上角. 最终在$ t = 12 $时刻, 主涡、二次涡和流速场达到稳定状态. 这些特征与文献[34]中的实验结果完全符合. 三维与二维顶盖驱动方腔流动的一个主要区别是沿翼展方向壁端涡旋的演化和对称面附近的Taylor-G?rtler型涡旋, 图6清晰地显示了时刻$ t = 8, 10, 12 $的这些突出特征.
图 5 三维顶盖驱动方腔流测试中均匀网格上($h = {1}/{128}$)对称平面$ z = 1 $内的流线快照. 不同颜色代表不同的$ \|{u}\|_2 $^[5]
Figure5. Snapshots of streamlines for the 3D lid-driven cavity test inside the symmetry plane $ z = 1 $ on a uniform grid with $h = {1}/{128}$. Different colors represent different values of $ \|{u}\|_2 $^[5].

图 6 三维顶盖驱动方腔流测试中均匀网格上($h = {1}/{128}$)涡量的x-分量等值面. 红色, 橙色, 蓝色和青色对应的涡量x-分量值分别为–0.50, –0.25, 0.25和0.50^[5]
Figure6. Isosurfaces of the x component of the vorticity vector for the 3D lid-driven cavity test on a uniform grid with $h = {1}/{128}$. The values of the vorticity component for the red, orange, blue, cyan surfaces are –0.50, –0.25, 0.25, and 0.50, respectively^[5].

本文还对GePUP-IMEX算法结果和 Guermond等^[34]的实验和计算结果进行了定量比较: 图7和图8分别给出了主漩涡的中心位置和在平面$ z = 1, 0.5 $内的流速剖面图. 可以看到, GePUP-IMEX算法的计算结果和Guermond等^[34]的实验和计算结果都非常符合. 此外, GePUP-IMEX算法的某些计算结果, 例如流速的v分量, 要比Guermond等^[34]的计算结果更加符合实验数据.
图 7 时刻$ t = 1, 2, 4, 6, 8, 10, 12 $对称平面$ z = 0 $上主旋涡的中心位置比较. “$ \circ $”表示本文的数值结果(由最小化$ \|{u}\|_2 $得到), 而“+”和“$ \square $”分别是Guermond等^[34]得到的实验和计算结果^[5]
Figure7. A comparison of the center locations of the primary eddy in the symmetry plane $ z = 0 $ at time instances $ t = 1, 2, 4, 6, 8, 10, 12 $. The circles represent numerical results of this work, which are determined by minimizing the speed $ \|{u}\|_2 $. The experimental and computational results of Guermond et. al.^[34] are represented by the crosses and the squares, respectively^[5].

图 8 三维顶盖驱动方腔流在时刻$ t = 12 $的流速剖面图比较($h = {1}/{128}$). 实线表示本文的数值结果, 即$ -\frac{1}{2}u $和$ \frac{1}{2}v $分别作为 $ \frac{1}{2}-y $和$ \frac{1}{2}-x $的函数. “$ \times $”和“$ \cdot $”分别表示Guermond等^[34]得到的实验和计算结果^[5]
Figure8. A comparison of velocity profiles for the three dimensional lid-driven cavity test ($h = {1}/{128}$) at $ t = 12 $. Solid lines represent computational results of this work, i.e. $ -\frac{1}{2}u $ as a function of $ \frac{1}{2}-y $ and $ \frac{1}{2}v $ as a function of $ \frac{1}{2}-x $. The computational and experimental results of Guermond et. al.^[34] are represented by the solid dots and crosses, respectively^[5].

5.1节验证了基于GePUP表述的INSE投影方法能够达到时空一致四阶精度. 在这个基础上, 本节的结果进一步验证了GePUP方法精确捕捉实际物理过程的能力.
2 -->

5.3.二维顶盖驱动方腔流

本节在倾斜正方形((87)式)上模拟二维顶盖驱动方腔流, 并将得到的数值结果与对应规则区域$ [0, 1]^2 $上的数值结果进行比较. 在由点$P = $$ \left(\dfrac{\sqrt{3}}{2}, \dfrac{\sqrt{3}}{2}+\dfrac{1}{2}\right)$和点$Q = \left(\dfrac{\sqrt{3}}{2}+\dfrac{1}{2}, \dfrac{1}{2}\right)$确定的域边界线段上, 流速场$ {u} $的边界条件为

$ \begin{array}{l} {u}({x}, t) = \left( \begin{array}{c} - {1}/{2} \\ {\sqrt{3}}/{2} \end{array} \right) \eta(r)(1-{\rm e}^{-\omega t}){\rm e}, \end{array} $

其中

$\begin{split} \eta(\alpha) =\;& \left \{ \begin{aligned} &{\rm e}^{-\tfrac{1}{1-\alpha^2}} & \text{ if } |\alpha| < 1, \\ &0 & \text{ otherwise,} \end{aligned} \right. \\&\;\; r = 2\|{x}-P\|_2-1, \quad \omega = 2. \end{split} $

其余计算域边界上$ {u} = {0} $, 而在初始时刻, $ {u} $在计算域内全为0, $ t_0 = 0 $, $ t_{\rm e} = 60 $, $ {{\textsf{Cr}}} = 0.1 $. 雷诺数Re = 1000, 均匀网格间距$h = {1}/{128}$.
图9和图10分别给出了在时刻$ t = 4, 8, 32 $时的流速场和涡量快照. 在规则区域和不规则区域上, 这两个物理量的特征完全一致.
图 9 二维顶盖驱动方腔流测试中均匀网格上$(h = {1}/{128})$的流速场
Figure9. Numerical results of the velocity field in the two dimensional lid-driven cavity test on the unit box and the rotated box with $h = {1}/{128}$.

图 10 二维顶盖驱动方腔流测试中均匀网格上$(h = {1}/{128})$的涡量快照
Figure10. Numerical results of the vorticity in the two dimensional lid-driven cavity test on the unit box and the rotated box with $h = {1}/{128}$.

2 -->

5.4.GePUP与一个经典二阶方法的性能比较

尽管在许多应用中^[35-37], 高阶方法都被认为要优于二阶方法, 但有****也报告了某些四阶方法的效率事实上是不尽如人意的^[38-40], 原因是这些“四阶方法”^[38-40]仅能在空间上达到四阶精度, 在时间上的精度仍为二阶. 由于求解误差主要是由二阶时间误差主导的, 这些方法在求解非定常流动时精度和效率都相对较低. 因此即使一个方法在空间上具有高精度甚至谱精度, 如果时间上只有二阶精度, 我们仍然认为该方法是一个二阶方法. 另一方面, 这表明要想获得一个真正的四阶方法并充分发挥其优势, 必须保证时间上的四阶精度. 现有文献中对四阶方法与二阶方法的比较缺乏系统的论述, 本节回顾文献[5]中给出的一个性能比较公式.
在比较具有不同收敛阶的INSE数值方法的效率时, 作以下两个假定:
(ECA-1) 求解线性方程组的时间远大于显式计算离散算子所耗费的时间;
(ECA-2) 求解单个线性方程组消耗的CPU时间与其未知量个数成线性关系.
基于这两个假定, 有:
命题5　 定义一个$ \xi $阶方法相对于一个$ \eta $阶方法的性能加速比为

$ S_{\xi|\eta}(\epsilon): = \frac{T_{\eta}(\epsilon)}{T_{\xi}(\epsilon)}, $

其中$ \xi, \eta\in{\mathbb{R}}^+ $, $\xi\geqslant\eta$, $ T_{\xi}(\epsilon) $和$ T_{\eta}(\epsilon) $分别是$ \xi $阶方法和$ \eta $阶方法达到给定精度$ \epsilon $需要的CPU时间. 则对于任一给定问题有^[5]:

$ S_{\xi|\eta}(\epsilon) = \frac{1}{n_{\rm T}} \left(\frac{E_{\eta}}{E_{\xi}}\right)^{\tfrac{ { {{D}}}+1}{\eta}} \left(\frac{E_{\xi}}{\epsilon}\right)^{\tfrac{(\xi-\eta)( { {{D}}}+1)}{\xi\eta}}, $

其中$n_{\rm T}$是在同一网格下求解该问题时 $ \xi $阶方法与$ \eta $阶方法所需CPU时间的比值, 而$ E_{\xi} $和$ E_{\eta} $分别是该固定网格下相应的误差, $E_{\xi}\leqslant E_{\eta}$.
证明　 设$ h_{\xi}(\epsilon) $和$ h_{\eta}(\epsilon) $分别是该$ \xi $阶方法和$ \eta $阶方法为达到给定精度$ \epsilon $所取的网格间距, 则

$\begin{array}{l} h_{\xi}(\epsilon) = \varTheta\left(\epsilon^{1/{\xi}}\right), \qquad h_{\eta}(\epsilon) = \varTheta\left(\epsilon^{1/{\eta}}\right), \end{array} $

其中$ f = \varTheta(g) $是指存在常数$ c_l, c_u > 0 $使得 $c_lg\leqslant $$ f \leqslant c_ug$.
根据假定(ECA-1)和(ECA-2), 可得

$ T_{\xi}(\epsilon) = c_\xi\left(\frac{1}{h_{\xi}(\epsilon)}\right)^{ { {{D}}}+1}, \quad T_{\eta}(\epsilon) = c_{\eta}\left(\frac{1}{h_{\eta}(\epsilon)}\right)^{ { {{D}}}+1}, $

其中$ c_{\xi} $和$ c_{\eta} $是两个正常数, $\left(\dfrac{1}{h_{\xi}(\epsilon)}\right)$和$\left(\dfrac{1}{h_{\eta}(\epsilon)}\right)$的幂次中的“+1”是由于初始值需要被推进$\varTheta\left({1}/{h}\right)$个时间步才能得到最终解. 由性能比$ S_{\xi|\eta}(\epsilon) $的定义(92)式, (94)式和(95)式可得

$ S_{\xi|\eta}(\epsilon) = \frac{T_{\eta}(\epsilon)}{T_{\xi}(\epsilon)} = c_{\xi|\eta}\epsilon^{-\tfrac{(\xi-\eta)( { {{D}}}+1)}{\xi\eta}}, $

其中$ c_{\xi|\eta} $对于该给定问题是个常数. 在相同网格上, $ E_{\xi} $和$ E_{\eta} $是运行$ \xi $阶方法和$ \eta $阶方法得到的精度, 而$n_{\rm T}$是这两个方法消耗CPU时间的比值. 假设$ \eta $阶方法的精度由$ E_{\eta} $提升到$ E_{\xi} $需要对网格加密$ n_{\eta\rightarrow\xi} $次, 且相邻两网格间距比率为r, 则

$ n_{\eta\rightarrow\xi} = \log_{r^{\eta}}\frac{E_{\eta}}{E_{\xi}},\tag{97a} $

$ \frac{T_{\eta}(E_{\xi})}{T_{\xi}(E_{\xi})} = \frac{r^{( { {{D}}}+1)n_{\eta\rightarrow\xi}}}{n_{\rm T}},\tag{97b} $

其中(97a)式中$ r^{\eta} $的幂次$ \eta $是 $ \eta $阶方法的收敛阶, (97b)式是由(95)式得到的. 在(96)式中取$ \epsilon $为$ E_{\xi} $, 再结合(97b)式可得

$ c_{\xi|\eta} = \frac{1}{n_{\rm T}} \bigl(E_{\xi}\bigr)^{\tfrac{(\xi-\eta)( { {{D}}}+1)}{\xi\eta}} \left(\frac{E_{\eta}}{E_{\xi}}\right)^{\tfrac{D+1}{\eta}}. $

将(98)式代入 (96)式即可得到(92)式.
利用(93)式, 可以对四阶投影方法GePUP-IMEX和Martin等^[7]提出的经典二阶投影方法MCG^②的效率作比较.
表2给出了在不同的测试中为了达到同样的$ L_2 $流速计算精度$ \epsilon $, GePUP-IMEX相对于MCG的性能加速比$ S_{4 |2} $. 对于较低的相对精度, 如$\epsilon = $$ 10^{-4}$, GePUP-IMEX相对于MCG在效率上的优势还不太明显, 而要达到适中的相对精度, 如$\epsilon = $$ 10^{-6},\; 10^{-8}$, GePUP-IMEX的性能优势是相当明显的. 尽管对于其他问题, $ S_{4 |2} $的具体值可能与表2中的结果不同, 但是$ S_{4 |2}(\epsilon) $是$ 1/\epsilon $的幂函数这一增长模式是始终保持不变的. 此外, 即使数值精度不是首要考虑的因素, GePUP-IMEX也允许使用更粗糙的网格来获得与MCG相同的精度, 因此在实际计算时GePUP-IMEX依然具有巨大的效率优势.

$ S_{4 |2} $	$ \epsilon = 10^{-4} $	$ \epsilon = 10^{-6} $	$ \epsilon = 10^{-8} $	$ \epsilon = 10^{-10} $
单涡量测试, Re = $ 2\times 10^4 $	1.73	5.50 × 101	1.73 × 103	5.47 × 104
二维粘盒测试, Re = $ 10^4 $	5.03	1.59 × 102	5.03 × 103	1.59 × 105
二维粘盒测试, Re = 100	1.26 × 101	4.00 × 102	1.26 × 104	3.98 × 105
三维粘盒测试, Re = $ 10^4 $	8.65 × 101	8.65 × 103	8.65 × 105	8.65 × 107
三维粘盒测试, Re = 100	4.24 × 101	4.24 × 103	4.24 × 105	4.24 × 107

表2不同的测试中为了达到同样的$ L_2 $流速计算精度$ \epsilon $, GePUP-IMEX相对于MCG^[7]的性能加速比$ S_{4 |2} $^[5]
Table2.The speedup $ S_{4 |2} $ of GePUP-IMEX over MCG^[7] for achieving the same $ L_2 $ accuracy $ \epsilon $ of the velocity^[5].

由于涡量$ \nabla\times{u} $是由流速场的一阶导数得到的, GePUP-IMEX和MCG对涡量的计算在$ L_{\infty} $范数下分别具有三阶和一阶精度, 表3给出了在不同测试中为了达到同样的$ L_{\infty} $涡量计算精度$ \epsilon $, GePUP-IMEX相对于MCG的性能加速比$ S_{3 |1} $. 即使对于很低的相对精度, 如$ \epsilon = 10^{-2} $, GePUP-IMEX相对于MCG在效率上的优势也是非常明显的.

$ S_{3 |1} $	$ \epsilon = 10^{-2} $	$ \epsilon = 10^{-3} $	$ \epsilon = 10^{-4} $	$ \epsilon = 10^{-5} $
单涡量测试, Re = $ 2\times 10^4 $	1.89 × 104	1.89 × 106	1.89 × 108	1.89 × 1010
二维粘盒测试, Re = $ 10^4 $	2.30 × 105	2.30 × 107	2.30 × 109	2.30 × 1011
三维粘盒测试, Re = $ 10^4 $	3.66 × 107	1.70 × 1010	7.88 × 1012	3.66 × 1015

表3不同的测试中为了达到同样的$ L_{\infty} $涡量计算精度$ \epsilon $, GePUP-IMEX相对于MCG^[7]的性能加速比$ S_{3 |1} $^[5]
Table3.The speedup $ S_{3 |1} $ of GePUP-IMEX over MCG^[7] for achieving the same $ L_{\infty} $ accuracy $ \epsilon $ of the vorticity^[5].

以上展示的测试均是在作者的个人电脑(Intel^? i7-3930, 138 GFlop/s)上进行的, 截至2020年11月,世界上最快的超级计算机是日本的 Fugaku (7630848 核, peak performance 537212 TFlop/s)^[43], 因此由个人电脑切换到超级计算机上进行计算所带来的性能提升最大为$ 4\times 10^6 $. 而注意到表2和表3中许多结果是大于$ 4\times 10^6 $的, 这说明为了在一些测试问题中达到给定精度要求, 在作者的个人电脑上运行四阶GePUP-IMEX算法比在国际上最快的超级计算机上运行二阶MCG方法都要快!
在数值模拟高雷诺数不可压流体的物理过程时, 往往要用到依赖于流速的一阶和二阶导数的物理量, 例如涡量和曲率等. 表3中的结果说明了四阶GePUP-IMEX投影方法对精确捕捉高雷诺数流体涡量场的演化过程具有重要的科学意义.
随着流体力学的飞速发展, 我们越来越需要研究具有多长度尺度和多时间尺度的复杂物理问题, 这对数值模拟在精度上和效率上都提出了越来越高的要求. 四阶和更高阶精度的算法极大地提高了计算精度和计算效率, 若进一步将这些算法与并行计算和自适应网格结合, 我们将能解决更多的实际物理问题.

本文首先回顾了一些求解不可压Navier-Stokes方程的经典数值方法, 然后基于广义投影算子和UPPE表述推导了GePUP表述和相应的数值方法. GePUP方法最突出的优势在于时间积分方法和空间离散格式是完全解耦的, 因此这两个算法模块都能够作为“黑匣子”进行处理. 时间积分方法的灵活自由选择便于实现INSE解在时间上的高阶精度, 并使得GePUP格式同时适用于高低雷诺数流体. 空间离散的灵活自由选择使得算法可以很好地适应不规则边界. 数值测试结果印证了该方法的四阶精度及其捕捉实际物理过程的能力, 再结合一个简单的效率比较公式就可以看出四阶GePUP方法相对于现有二阶投影方法在效率上和精度上都具有巨大优势.
目前我们已经将GePUP方法推广到静止不规则区域上, 我们计划将其与高阶界面追踪方法^[44,45]耦合, 针对动边界不可压流体提出一个时空一致四阶精度的投影算法.