基于Lattice Boltzmann方法的洪水演进模拟

刘宁宁，时 勇，孙华林

（山东省海河淮河小清河流域水利管理服务中心，济南 250100）

0 引 言

洪水演进是指根据河段上断面的洪水过程推求河段下断面未来的洪水过程，其理论依据是洪水波在河道中的传播规律。天然河道中的洪水波属于浅水长波，它的运动规律可以采用圣维南方程组进行描述，但圣维南方程组属于双曲型拟线性的偏微分方程，直到目前还不能给出这类方程的解析解，只能用数值方法近似求解，不仅计算复杂，而且需要详细的河道地形资料。因此，常在实践中采用近似的计算方法，即采用圣维南方程组的简化形式。圣维南方程组的简化形式多采用运动波方程和扩散波方程。经研究发现，扩散波不考虑惯性项但考虑附加比降项，能反映回水影响，能用于平原河网地区，其适用范围相较运动波更广［1-3］。

Lattice Boltzmann 方法是20世纪80年代末期提出的一种模拟流体流动［4］的数值计算方法，是建立在介观动理学模型基础上，以粒子速度分布函数为研究对象，建立微观粒子与宏观量之间的联系。它打破了传统数值计算方法的理论框架，不再受连续介质条件的束缚，并在求解非线性偏微分方程组时，将非线性偏微分方程组转化为Lattice Boltzmann 演进方程，从而变为简单的线性方程，使得计算效率得到大大提高。这些特点促使Lattice Boltzmann 方法在模拟大规模流动和求解非线性方程中得到了广泛的应用［5-8］。在水力学方面：程永光等［9］采用Lattice Boltzmann 方法模拟了一维、二维明渠非恒定流，并模拟了实际工程中二维圆形溃坝波的传播；张东辉等［10］通过采用Lattice Boltzmann 方法求解Richards 方程，对土壤下渗过程进行了模拟；吴家阳［11］提出了基于Lattice Boltzmann 方法的河库瞬变流模拟方法，并与GPU 结合进行并行计算，使得计算效率大大提升。

在水文学中，所谓“物质特性”的传输，表面上是流速、流量、浓度的改变，但其物理本质却是质量、动量、能量的传输，从微观角度来理解，实质上就是组成宏观量的粒子的碰撞和迁移，因此，Lattice Boltzmann 方法应用于水文学中，具有较强的物理基础。扩散波方程是修正的Burgers 方程参数取特定值时的特例，是其在水文学中的应用，近年来，Lattice Boltzmann 方法在求解Burgers 方程方面有较多的研究［12-14］。本文将Lattice Boltzmann方法运用到扩散波方程中，通过算例进行对比分析。

1 扩散波方程

河道中的洪水波运动可以采用扩散波方程进行描述，当河段下边界为不受回水顶托影响的自由下边界，初始流量为零，且不考虑旁侧入流时，可用下列定解问题描述：

式中：Q为流量；Cd为扩散波的波速；μ为扩散波的扩散系数；x为河长；t为时间；I（t）为河段上断面的入流过程；L为河段长度。

2 Lattice Boltzmann 模型

2.1 Lattice Boltzmann方程

Lattice Boltzmann 方法是将空间划分为均匀的网格，并在每个节点上设置粒子分布函数，根据粒子分布函数在网格上的运动特性，碰撞项采用BGK 近似［15］，得到Lattice BGK方程：

式中：fα为流体粒子的分布函数；Ωα为碰撞项；τ为弛豫时间；为平衡态分布函数；α为方向，与选取的离散速度模型有关。

Lattice Boltzmann 方法认为在Δt时间内粒子能以恒定的速度e→沿α方向从一个节点迁移到下一个相邻节点，并在这个节点上与原有粒子发生碰撞，取代原有分布函数，其碰撞和迁移过程分别由式（3）和（4）确定：

碰撞过程：

在碰撞过程中，碰撞项需要满足质量和动量守恒条件：

2.2 扩散波方程的Lattice Boltzmann 解法

构造Lattice Boltzmann 模型的关键在于选择合适的平衡态分布函数，而平衡态分布函数的具体表达式又与离散速度模型的构造有关［16］。Qian 等人提出的DnQb模型是最具有代表性的离散速度模型，其中，Dn表示空间维数，n可取1、2、3；Qb表示离散的速度数，b可取3、5、7、9、15、19 等［6］。本文尝试采用D1Q5速度模型，以期提高计算精度。D1Q5 为一维五速模型，即在每个网格点上存在5 类不同速度的粒子，粒子速度=c[0，-1，-2，1，2］，α=0，1，2，3，4，c为网格步长Δx与时间步长Δt之比，即：c=ΔxΔt。同时，采用对空间坐标不进行多尺度处理，而对时间坐标引入五个时间尺度的多尺度展开方式［10］，得到：

式中：ε为Knudsen 数，表示粒子平均自由路径与宏观流动中特征长度的比值。

并将Chapman-Enskog展开式应用于分布函数，得到：

对式（2）的第一项进行二阶Taylor展开，并将式（6）和（7）代入，最终可得到五个不同时间尺度的离散Boltzmann方程：

由于在Chapman-Enskog 展开中假设局部范围内已接近平衡态，故对式（7）两边求和，得到：

为建立D1Q5 模型的平衡态分布函数与扩散波方程中流量这个宏观量之间的联系，并检验多尺度展开方式选取的正确性，需将五个不同时间尺度的离散Boltzmann 方程通过求解各阶矩的方式还原到式（1）。具体步骤详见文献［17］可得：

将式(9)+ε×式(10)+ε2×式(11)+ε3×式(12)五个方向求和得：

与宏观方程式（1）相比，当k=时，式（15）只增加了一项o(ε4)，可以通过τ来控制。这样就得到了具有四阶截断误差的扩散波方程。D1Q5 模型对应的扩散波方程的平衡态分布函数由式（14）求解，得到：

根据宏观初始流量，由式（16）计算各节点平衡态分布函数，通过粒子的碰撞和迁移过程，更新每个节点的分布函数，并通过式（14）即可得到不同时刻各节点的宏观流量。

3 应用实例

此次采用两个算例，分别取自文献［18］和［19］。算例一：长江上游龙街—乔家河段，河段长247 km，平均河床坡度0.001 12。1967年6月18日至20日，在河段上游和下游两端的两个站点测量了一次洪水事件。根据龙街和巧家两水文站的有关资料得到波速Cd=5.92 m/s，扩散系数μ=19 456 m2/s。算例二：湖南资水柘溪水电站下游江南—夫溪河段，河段长20.1 km，1967年5月13日至14日，发生一次洪水过程。根据江南和夫溪两水文站的有关资料得到波速Cd=2.792 m/s，扩散系数μ=8 376 m2/s。

3.1 模拟结果

采用Muskingum 法、Lattice Boltzmann 方法的D1Q5 模型以及解析解法对所给算例进行求解。算例一：对于Muskingum法，为了保证长河段条件下的线性条件，需将长河段分成4 个子河段即N=4，作分段连续演算，子河段xe=0.45，Ke=Δt=3 h。对于解析解法，为保证计算精度，划分矩形条块时，取计算时段Δt=3 h。对于Lattice Boltzmann 方法的D1Q5 模型，取弛豫时间τ=1.5，空间步长Δx=1 000 m，时间步长Δt=10 s。算例二：对于Muskingum 法，将长河段分成4 个子河段即N=4，作分段连续演算，子河段xe=0，Ke=Δt=0.5 h。对于解析解法，取计算时段Δt=0.5 h。对于Lattice Boltzmann 方法的D1Q5 模型，取弛豫时间τ=1.5，空间步长Δx=402 m，时间步长Δt=5 s。上述3 种解法所得计算结果及精度分析见表1和图1。

由表1和图1可以看出，这3种方法计算得到的洪峰相对误差、峰现时间、确定性系数基本一致，这表明Lattice Boltzmann方法可以应用于河道洪水波的演算。

图1 3种方法计算与实测流量过程比较Fig.1 The comparisons of observed and computed hydrograph

表1 3种方法计算精度统计表Tab.1 The comparisons of observed and computed values

3.2 影响因素分析

Lattice Boltzmann 方法的计算精度主要受到弛豫时间、步长选取、速度模型、边界格式等因素的制约［10］。在水文预报中，判断预报精度的指标有4 个：①洪峰误差；②径流量误差；③峰现时间；④确定性系数。洪峰、径流量的误差在20%以内，峰现时间在一个计算时长之内，确定性系数在0.7 以上，即认为是合格的预报，对于误差的分析，4 个指标应综合考虑。经计算，步长和弛豫时间的选取对峰现时间没有影响，径流量误差可以由确定性系数反应，故以洪峰误差和确定性系数作为分析对象，将Lattice Boltzmann方法的计算结果与实测流量进行比较。

3.2.1 步长对计算精度的影响

首先分析空间步长对计算精度的影响。选择弛豫时间τ=2，时间步长Δt=5s，综合考虑河长、计算效率、计算稳定性等因素，对于算例一，分别取空间步长Δx=500 m、Δx=1 000 m、Δx=2 470 m、Δx=4 940 m、Δx=12 350 m 进行计算；对于算例二，分别取空间步长Δx=402 m、Δx=804 m、Δx=1 005 m、Δx=1 340 m、Δx=2 010 m进行计算，计算结果详见表2。

表2 不同空间步长对计算精度的影响Tab.2 The influence of different space step on calculation accuracy

其次分析时间步长对计算精度的影响。选择弛豫时间τ=2。综合考虑河长、计算效率、计算稳定性等因素，对于算例一，在空间步长Δx=2 470 m 的情况下，分别取时间步长Δt=1 s、Δt=2 s、Δt=5 s、Δt=10 s、Δt=20 s进行计算；对于算例二，在空间步长Δx=1 005 m 的情况下，分别取时间步长Δt=1 s、Δt=2 s、Δt=5 s、Δt=10 s、Δt=20 s进行计算，计算结果详见表3。

由表2～表3可知，时间步长比空间步长对计算精度的影响要小一些。这是因为，空间步长的选取直接影响到相同河段长情况下，网格划分的多少，进而影响粒子的碰撞迁移次数。因此，在采用Lattice Boltzmann 方法进行洪水演算时，应结合上下游断面长度，选取合适的空间步长。同时，在计算效率允许的情况下，也应综合考虑粒子速度（Δx/Δt）的取值，因为粒子速度（Δx/Δt）取值太小，会导致计算不稳定。

表3 不同时间步长对计算精度的影响Tab.3 The influence of different time step on calculation accuracy

3.2.2 弛豫时间对计算精度的影响

弛豫时间τ是一个无纲量参数，表示粒子分布函数趋于平衡态的时间。在理论上，为了满足Lattice Boltzmann 方程的线性稳定性，弛豫时间τ的取值应大于0.5［20］。但从计算精度的角度出发，对于扩散波方程的Lattice Boltzmann 方法D1Q5 速度模型，应分析弛豫时间τ在哪个取值范围，可以使计算误差最小。为讨论弛豫时间τ的取值对计算精度的影响，应在特定的空间和时间步长情况下，选取不同的弛豫时间τ进行计算分析。对于算例一，选取空间步长Δx=1 000 m，时间步长Δt=5 s；对于算例二，选取空间步长Δx=402 m，时间步长Δt=5 s。计算结果详见表4。

表4 不同弛豫时间对计算精度的影响Tab.4 The influence of different relaxation time on calculation accuracy

由表4可知，弛豫时间τ的取值过大，会导致误差增大，取值过小，在特定的粒子速度（Δx/Δt）下，会导致计算不稳定。因此对于洪水演算，建议弛豫时间τ的取值取在［1.5，3］范围内为宜。

4 结 论

本文选用扩散波方程来描述河道洪水波运动，采用Lattice Boltzmann 方法的D1Q5 模型来求解该偏微分方程，并给出了详细的多尺度处理和分布函数确定步骤。采用两个算例来说明Lattice Boltzmann 方法的应用，同时采用Muskingum 法和解析解法来验证Lattice Boltzmann 方法的准确性，证明了Lattice Boltzmann 方法可以很好地求解扩散波方程和预测洪水过程。此外，分析了步长和弛豫时间的选取对计算精度的影响，结果表明，空间步长的取值比时间步长的取值对精度的影响大，建议弛豫时间τ的取值取在［1.5，3］范围内为宜。