浅水方程大时间步长离散方法研究进展

许仁义，赵 磊，王远见，徐 波

（1.扬州大学 水利科学与工程学院，江苏 扬州 225009；2.黄河水利科学研究院，河南 郑州 450003）

在实际的工程项目中，对数值模拟计算效率的需求是永无止境的。比如上游发生溃坝，洪峰何时到达下游城市，迅速划出人口撤离范围；或者上游发生化工污染事件，污染物何时会到达下游取水口。这需要流域管理单位根据模拟软件的模拟结果迅速做出反应。这时，模拟速度就显得尤为重要。模拟计算速度的提升方法有两个方面，一方面是提升硬件水平，在科技日新月异的今天，计算机的主频提升速度相对放缓，硬件方面可以通过多核来提升计算效率；另一方面是追求好的高效算法，大时间步长格式（Large Time Step，LTS）［1］就是其中一种，它的主要思想是在一个时间步长内完成传统格式在多个时间步长内完成的工作，这种算法被证明能大幅提升计算效率。本文第一部分从线性标量双曲线方程求解开始，介绍了大时间步长格式的基本思想和方法，然后将其扩展到非线性双曲方程组求解，以及二维非结构网格上的大时间步长格式更新方法；第二部分提出了目前大时间步长格式存在的问题以及解决方法，最后对其未来进行了展望。

1 离散方法

1.1 大时间步长格式的引出——标量方程的大时间步长格式

针对线性波动方程（双曲型偏微分方程）：

式中：u为振幅；a为波的传播速率，大约为330 m／s。

如果采用如下一阶迎风型显格式离散方法：

那么，需要满足CFL定理，柯朗数小于等于1：

如图1所示，在i单元与i－1单元的界面xi－1／2处产生的波经过一个时间步长Δt后到达ε点。采用式（2）离散时，这个ε点需要落在i单元内。

图1 传统格式单个时间步长内单波传输只限于本单元内

波在时间步长Δt内所行走的距离必然不能超过一个单元的长度，也就是空间步长Δx。

而大时间步长格式是指考虑波在一个时间步长Δt内穿过多个单元的情况，如图2所示。

图2 LTS格式单个时间步长内单波传输穿过多个单元

此时式（2）的更新方法不再适用。把式（2）改变一个写法：

式中：ζ为大时间步长格式的更新系数。

从式（5）可以看出，从n时刻到n＋1时刻，在这个i－1／2处发出的单波的影响下，变化量是），这个变化量由两部分组成，其中一个是：

Δu为波强，在传播过程中是不会变化的。波速a乘以时间Δt得到单波所行进的路程，除以空间步长，结果就是行进路程占整个空间步长的比重，在使用传统格式时，这个就是柯朗数，需要满足小于1。当扩展到大时间步长格式时（见图2），在界面xi－1／2处产生的单波经过一个时间步长Δt后，完整地穿过了i、i＋1单元，最终到达的ε点落在i＋2单元内。由前面的定义，当单波完整地穿过某个单元时，波在单元内所行进的路程与单元的空间步长是相等的。

floor是舍去法求整函数。比如Cr是2.6，那么floor（Cr）就是2，代入上式得到ζi＋2＝0.6。因此下一时刻的更新方法就变成：

这就是由Leveque［1］提出的大时间步长格式。这种方法的优点是：把原来几个时间步长才能实现的计算工作，在一个时间步长内实现，极大地提升了计算效率。

1.2 非线性双曲守恒律的大时间步长格式

对于非线性的双曲守恒律，以浅水方程为例：

式中：h为水深，m；u为流速，m／s；g为重力加速度，取9.8 m／s2。

由标量方程实现大时间步长格式的计算，需要知道波速来求得单波在一个时间步长内穿过了多少单元，所经过的单元都要进行更新；需要知道波强以求得这些单元需要怎样的更新。而这两个量的求得也有两种方法，第一种方法是黎曼近似解，就是通常所说的Roe方法［2］。对于相邻两个单元先做Roe平均：

由平均值求得雅可比矩阵：

再对雅可比矩阵进行特征化：

求得两个单元间的界面所发出波的速度：

求波强：

有了波强和波速，就可以按照式（9）对每个单元进行更新。

第二种方法是黎曼精确解，这种方法需要迭代。对于相邻的两个单元，有着不同的初始变量，可以理解为一个局部间断，那么这就是一个黎曼问题。非线性双曲守恒系统的黎曼问题的解有4种形式［3］，它们是稀疏波与激波的4种排列组合：①左右两个都是激波；②左稀疏波右激波；③左激波右稀疏波；④两侧都为稀疏波。问题求解的关键是求得中间状态。如果是激波，则满足条件：

式中：S为激波的速度。F和U见（10）式，下标S表示中间状态，下标L表示左边，下标R表示右边。

如果左边是稀疏波，则：

如果右边是稀疏波，则：

通过式（17）～式（19）可以求出中间状态，与两侧已知的左右状态的差值即为波强，而波速也可以通过上面的关系式求出。

1.3 含底坡源项的大时间步长格式

含底坡源项的浅水方程，若采用简单的中心格式则会产生虚假流动，传统格式为解决这个问题采用的是和谐（well⁃balanced）格式［4－6］。大时间步长格式对此也有两种方法，第一种方法是黎曼近似解，对底坡源项进行特征分解［7］，把方程（10）变成：

式中：z为河底高程。

此时雅可比矩阵的右矩阵为

此时虽然由前面的两个单波变成3个，但是其中一个波速为0，另外两个波速不变：

波强为

那么：

第二种方法同样是黎曼精确解。采用Godonov方法离散时，两个相邻单元的物理状态是不连续的，同样在非齐次的浅水方程中，底坡也是不连续的，像楼梯的台阶一样，我们称之为阶梯黎曼问题（Step Riemann Problem，SRP）。相比于齐次黎曼问题，阶梯黎曼问题复杂许多，因为在阶梯位置出现了一个静态的激波，所以待求的中间状态变成了两个，遵循的关系依然是式（18）～式（20）。

1.4 二维大时间步长格式

将一维格式拓展到结构网格的二维大时间步长格式，采用Strane维分裂的方法［8］：当时间步长n是奇数时，x方向上前进半个时间步长，然后y方向上前进一个时间步长，最后回到x方向上再前进半个时间步长；当时间步长n是偶数时，y方向上前进半个时间步长，然后x方向上前进一个时间步长，最后回到y方向上再前进半个时间步长。

如果是二维非结构网格，问题就会变得复杂［9－10］，还是按照前面的方法，先求出波速和波强。对于二维浅水方程：

式中：v为横向流速；＝icosθ＋jsinθ为单元边的外法线方向。

两个相邻的单元做Roe平均得到：

于是得到雅可比矩阵：

对其进行特征化：

左矩阵为

右矩阵为

波速为

波强为

于是得到：

得到波速和波强的表达式后，下面来看非结构网格上的大时间步长格式是如何实现的。

如图3所示的非结构三角形网格，在三角形单元Cw上一边w产生的单波，如果是传统格式，只会影响Cw单元本身，而在大时间步长格式里还要扩展到Cw单元的相邻单元C1和C2。与一维大时间步长格式不同的是，在波的传播方向上，一维格式的相邻单元只有一个，而二维三角形网格有两个，这就需要对波强进行分配。

图3 二维非结构网格上的LTS

2 存在的问题

2.1 非线性问题

大时间步长格式从线性标量方程扩展到非线性双曲守恒律遇到的第一个问题是非线性的问题，在大时间步长格式中，每个单元需要经历不同界面的波，如图4所示，如果波是线性的或者非线性比较弱，那么只需要将这些波的波强简单累加就可以；如果非线性比较强，那么不同波相遇后产生新波的波速会发生变化，如图5所示。空气动力学的欧拉方程中，是需要考虑这种非线性叠加的［11－12］，而在浅水方程中，考虑这种非线性带来的变化几乎可以忽略［13］。

图4 多个单波同时穿过一个单元

图5 非线性单波叠加

2.2 稀疏波的分解

采用Godonov型格式进行离散时，每两个单元间的界面上是一个局部间断，就是一个黎曼问题。对非线性双曲系统，这个黎曼问题会产生稀疏波和激波。我们通常使用的Roe［2］、Osher［14］、HLL［15］等格式是指对这个黎曼问题的近似解。在这些传统格式中，忽略了稀疏波的存在，将其视为激波。因为时间步长不够大，所以稀疏波的宽度扩张不是很严重，如图6中的红色部分。但是在大时间步长格式中，这个问题就凸显出来。如果仍然采用传统方法不做任何处理，就会出现稀疏波分裂的情况（如图7所示）。

图6 时间步长增大稀疏波波头波尾距离变大

图7 稀疏波断裂

Morales－Hernández等［16］在黎曼近似解中强制对单波进行分解，也取得了不错的效果。另一种方法就是采用黎曼精确解。这种方法在传统的格式中没有被广泛使用的原因是在每个局部间断需要做Newton迭代来取得精确解，而且迭代的步数无法控制，可能造成迭代步数过多甚至发散。而这种方法最大的好处就是可以完美地求解出稀疏波，得到波头和波尾的速度，使得稀疏波的分解变得简单。实践证明，大时间步长格式中迭代步数一般为3步左右，而对于柯朗数动辄几十的大时间步长格式来说，这个计算成本完全可以抵消［17］。

2.3 底坡源项

浅水方程有底坡和糙率两个源项，使得原来的齐次双曲偏微分方程组变成了非齐次，而大时间步长格式解决这个问题有两大类方法。第一大类是黎曼近似解，因为和谐格式是用物理通量平衡源项，所以这种方法是不能采用传统和谐格式的，而此方法中，物理通量没有被计算，那只有对底坡源项进行特征分解，将其融合到波强中，随着单波一起穿过多个单元［18］。

第二大类是黎曼精确解，这类方法比较复杂。首先求解的理论就分为两种，一种是基于动量和质量守恒的方法［19］，一种是基于能量质量守恒的求解方法［20］；其次是对底坡源项间断的处理方法，一种认为这是一个垂直间断［21－22］，一种认为这是一个单调连续的变量［23－24］（宽度趋向于0）。阶梯黎曼问题在很多情况下不能得到唯一解，只有通过附加其他条件才能使得解唯一［25－26］。因为求解的出发点不同，所以这些方法求得的结果不能统一。

2.4 振 荡

振荡问题可以说是大时间步长格式的顽疾，采用随机选取法（Random Choice Method，RCM）只能在一定程度上抑制振荡［27］，随着柯朗数的增大，这种方法也逐渐失去了效用。在计算机程序中，随机数不是真正的随机数，而是固定的。Wang等［28］用范德科皮特（Van der Corput）序列来制造随机数，效果没有什么变化。如果采用滤波法来抑制振荡，这显然需要预先知道振荡发生的位置和幅度，在工程中不太可能实现。所以，到目前为止对这个问题还没有一个很好的解决方法。

3 总结与展望

大时间步长格式从提出至今已有40多a时间，最近几年仍然有很多成果，展现了其强大的生命力。从另一方面来看，它仍然有问题没有得到解决。其基本思想就是通过单个时间步长内单波穿过多个单元，在一个时间步长内完成了传统格式多个时间步长的计算，通过这种方法可极大地提高计算效率，这在线性波动方程里是显而易见的，但是到了非线性的方程组里就会出现一系列问题，其中最重要的问题就是振荡，这种振荡与传统格式的振荡在现象上有所区别，不是发生在激波附近，而是分布在整个平台。其解决方法也不同于传统格式，这是个顽疾，如果这个问题得到解决，大时间步长格式将会得到更广阔的应用。