滑坡涌浪高性能有限元計算軟件研發

(群馬大學大學院理工學府， 日本 桐生 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，開發了涌浪高性能有限元分析程序，并用一些理論解和實驗結果驗證了開發的程序。垂直水壩潰壩的理論解和數值計算結果，以及室內實驗結果和數值計算結果的比較表明研發的涌浪高性能有限元軟件可以快速得到可靠的計算結果。