滑坡涌浪高性能有限元计算软件研发

(群马大学大学院理工学府， 日本 桐生 376-8515)

滑坡、泥石流及冰川等高速运动体入水冲击水库、湖泊、河道及海湾等封闭水域内的水体时会诱发涌浪，并由此产生大规模的灾害。其代表性的实例为1963年意大利瓦依昂水库库岸发生大面积整体岩质滑坡，滑坡体长2 km，宽约1.6 km，滑坡体积达2.4亿m3，下滑运动速度达15～30 m/s的滑坡体突入水库时产生高达90～130 m的涌浪，给库岸及大坝下游集镇造成毁灭性破坏[1]。在日本2011年7月6号台风时奈半利川平锅大坝上游流入大量的泥石流，引发超过5 m高的涌浪，给大坝闸门造成损伤并产生漫顶，同时损坏了大坝上游的平锅吊桥。紫坪铺大坝为高158 m，库容量达11.1亿m3的混泥土面板堆石坝，它与2008年5月12日发生的汶川地震震中距离只有17 km。汶川地震引发了紫坪铺滑坡，滑体中心到库水面的高差达700 m，约45万m3的滑坡岩体高速冲入水库，产生高达25 m的涌浪，大约70名湖边垂钓人被涌浪卷入水库造成死亡，10台汽车也被卷入湖中[2]。表1列出了一些典型的涌浪及其造成的人员伤亡情况。

表1 历史性的滑坡或部分水下滑坡产生的涌浪

为防止涌浪灾害需要能确立涌浪特性及其规模等的预测方法。浅水长波方程常用作涌浪的控制方程。因为浅水长波方程中包含有移流项，数值计算的稳定性很差，有必要解决数值计算的稳定问题。另外，有限元离散得到的总系数矩阵是非对称的，因此大规模涌浪分析时为缩短计算时间，需要开发线性方程组的多核多线程并行解法。本研究中，利用Intel MKL提供的共享内存机器上实现的稀疏线性方程求解器PARDISO，开发了涌浪高性能有限元分析程序，并用一些理论解和实验结果验证了开发的程序。

1 涌浪有限元分析

1.1 控制方程

用竖直方向平均流速表示的浅水长波方程包括连续方程(1)、x和y方向的运动方程(2)和(3)[3-4]。

(1)

(2)

(3)

式中：η为水面标高；z为地表面标高(见图1)；水深h=η-z；x方向流束U=hu；y方向流束V=hv；g为重力加速度；n为曼宁系数；ε为水深平均涡动黏性系数。为考虑滑坡运动引起的地表面标高的变化，与一般的浅水长波方程有所不同，式(1)中同时引入了水面和地表面标高。

图1 涌浪示意图

1.2 稳定化有限元离散

用有限元对上述3个控制方程进行有限元离散，并引入SUPG(Streamline upwind Petrov Galerkin)项[5-6](其物理意义见图2)和Shock capturing项，则连续方程可表示为

(4)

x方向的运动方程可以表示为

(5)

(6)

图2 权重函数(左)Galerkin有限元法和(右)SUPG有限元法

同理可得y方向的运动方程。

上述方程如用矩阵表示，则连续方程表示为

(7)

x和y方向的运动方程可分别表示为

(8)

(9)

式中，η、U、V分别为各节点水面标高x和y方向流束组成的向量。另外，

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

另外，τSUPG为一维稳定化参数[7]扩展到多维条件下的稳定化参数[8-9]。

(26)

(27)

(28)

(29)

单元尺寸可按下式计算[10]

(30)

Shock capturing系数可按下式定义[11]

νSHOC=τSHOC(uint)2

(31)

(32)

1.3 时间的直接积分

考虑时间间隔为Δt(Δt=tn+1-tn)的两个状态tn+1和tn时刻的水面标高、x和y方向流束向量分别为ηn+1、Un+1、Vn+1和ηn、Un、Vn，则水面标高、x和y方向流束向量的时间微分可表示为：

(33)

(34)

(35)

tn+θh时刻的水面标高、tn+θu时刻x和y方向流束向量的时间微分可表示为：

ηn+θh=θhηn+1+(1-θh)ηn

(36)

Un+θu=θuUn+1+(1-θu)Un

(37)

Vn+θu=θuVn+1+(1-θu)Vn

(38)

式中，θh和θu为时间积分参数，θh和θu≥1/2时时间积分是无条件稳定的。本文取θh=θu=1/2，即用Crank-Nicolson法进行时间积分。将式(33)─(38)代入式(7)─(9)，整理可得：

(39)

式中：

(40)

A12=θu(-B1+B1S)

(41)

A13=θu(-B2+B2S)

(42)

A21=θh(G1+G1S)

(43)

(44)

A31=θh(G2+G2S)

(45)

A33=A22

(46)

(47)

(48)

(49)

(50)

式中，Δtn-1和Δtn分别为时步n-1和n的时步长。当这两个时步长相等时，式(50)可简化成等时步Adams-Bashforth法。

2 多核多线程并行求解器PARDISO

2.1 求解器PARDISO简介

空间与时间离散得到的大规模稀疏矩阵作为系数的一次联立线性方程组Ax=b的求解为涌浪有限元分析的中心任务之一，因此采用适合于稀疏线性方程组的高速鲁棒的求解器就非常重要。直接法对迭代法不收敛的问题也可能求解，特别是矩阵对称正定时，如果没有数值误差则一定可以求解。因为这些特长，直接法得到了广泛的应用。通过考虑矩阵变带宽一维存储内部的零元素，可进一步减少计算量和节约内存的稀疏线性方程组直接求解法现在已成为主流。包括PARDISO[12]在内的稀疏线性方程组直接法求解器的计算步骤(图3)一般包括重排序、符号分解、LU分解、前进后退代入等4步的顺次计算。以下简单介绍一下各步的计算。

Ⅰ：矩阵A的非零元素位置发生变化时 Ⅱ：矩阵A的非零元素位置不变但值发生变化时 Ⅲ：只有b发生变化时

重排序就是用合适的置换矩阵P，使LU分解PAPT时填入元素尽可能的少。这里填入元素是指LU分解前为零但分解后为非零的元素。重排序的方法有等幅缩小(最小次数法，Reverse Cuthill-McKee法)，三角化(Markowitz法，Tewarson法)，分块化(Stewart法，Nested Dissection法)等许多算法[13]。PARDISO用METIS算法包中最小次数算法或者Nested Dissection算法进行重排序。

符号分解并不进行具体元素的分解计算，而只着眼于矩阵A的LU分解时非零元素的分布形式，找出LU分解后的非零元素的位置。由此算出LU分解所需要的内存量和计算量，确保LU分解后保存非零元素所需的内存，并记录非零元素的位置。为高效进行符号分解，利用列消去树的概念，把问题归为有效图的路径探索问题，从而使高速计算成为可能。

用符号分解时确保的内存量进行实际的LU分解。作为主要的LU分解方法，参照更新列的右侧的right-looking算法，左侧的left-looking算法广为人知。本程序所用的大规模稀疏线性方程组的多核多线程并行求解器PARDISO组合利用了left-looking和right-looking算法以实现高效并行计算。

用LU分解后的下三角形矩阵L进行前进代入，用上三角形矩阵U进行后退代入可求得解x。如果反复求解系数矩阵A相同的一次联立线性方程组(例如，一个加载步用修正Newton-Raphson法进行计算时)，则只要迭代进行前进后退代入计算就可求解。系数矩阵的非零元素的构造相同时则必须返回到LU分解进行计算。边界条件或分析范围等变化时则必须重新合成系数矩阵并从重排序开始求解过程。

2.2 超级节点

所谓超级节点为上三角矩阵L中全为非零，各列有相同的非零构造的列的集合。图4所示{1, 2}，{3, 4}，{5}，{6, 7, 8}分别为阶数2，2，1，3的超级节点。left-looking算法引入超级节点后，可促进分块化，也就是数据访问的局部化，可大幅提高阶层构造内存计算机的计算速度。把不同的非零构造的零元素假定成非零元素，生成超级节点，在有些情况下可以得到更有效的超级节点。

多核多线程并行计算稀疏线性方程组求解器PARDISO在消去树的并行化、节点水平的并行化、数据通道处理的并行化3个水平上实现了并行化。

图4 稀疏矩阵A(左)、L和U的非零构造及超级节点(右)

PARDISO用于涌浪有限元计算时需要生成线性方程组系数矩阵A的非零元素的位置列表。PARDISO使用以行为主的存储方法，即变形CSR(compressed sparse row)形式，对称矩阵只存储上半三角元素。该方法以行为单位存储每个非零数据。PARDISO对一个稀疏矩阵A的存储包括了3个数组：

1)values-矩阵A的非零元素。矩阵A的非零数据通过下面的columns与rowindex映射到values数组中。

2)columns-values中每个元素所在矩阵的列。

3)rowindex-给出每一行的元素在values中的位置。

具体可参照参考文献[14]。

3 计算实例

用研发的软件首先计算了图5所示长20 m的水槽中水位差0.8 m的涌浪传播。这个计算实例为静止的水坝瞬时坍塌放流的现象。为与完全流体的理论解进行比较，计算中假定水深平均涡动黏性系数ε=0，另外水槽侧壁假定为slip条件。图5所示为水坝坍塌1 s后的水深和流速的数值计算结果和理论解的比较。图5表明计算结果和理论解非常一致。

注：(左)1 s后的水深，(右)1 s后的流速

为确认移动边界问题的计算精度，计算了水槽长10 m，中心左侧水位高0.2 m，右侧为无水干床水槽中的水坝瞬时坍塌问题。同前例，计算中假定水深平均涡动黏性系数ε=0，水槽侧壁为slip条件。图6所示水坝坍塌1 s后水深计算结果和理论解的比较。由图6可知干床水槽条件下计算结果和理论解也非常一致。

图7所示为长7.2 m水深0.2 m的水槽内长2.44 m的底面活塞向上运动形成涌浪的实验[15]。活塞向上运动速度ζ(t)由下式表示

ζ(t)=ζ0(1-e-αt)

(51)

式中：ζ为活塞运动位置；ζ0与α为活塞运动控制参数，模型试验中ζ0=0.2h，α=0.235 8。

图6 干床水槽中水坝坍塌1 s后水深的计算结果与理论解的比较

图7 水槽及造波活塞

计算条件为水深平均涡动黏性系数ε=0，曼宁系数n=0，时步长Δt=0.014 29(无量纲表示)，采用0.2 m均匀网格。计算结果和实验结果的比较见图8，图8(上)为活塞左端，(下)为活塞右端的水位随时间的变化过程。图8表明计算结果和实验结果相当一致。

图8 活塞造波实验结果(左图实线)与计算结果(右图)的比较

4 结 论

高速运动的滑坡、泥石流及冰川等入水冲击水库、湖泊、河道及海湾等封闭水域内的水体时产生涌浪的预测经常是个大规模的计算问题。为缩短大规模涌浪分析的计算时间，本研究利用Intel MKL提供的共享内存机器上实现的稀疏线性方程组求解器PARDISO，开发了涌浪高性能有限元分析程序，并用一些理论解和实验结果验证了开发的程序。垂直水坝溃坝的理论解和数值计算结果，以及室内实验结果和数值计算结果的比较表明研发的涌浪高性能有限元软件可以快速得到可靠的计算结果。