滑坡涌浪高性能有限元計(jì)算軟件研發(fā)

(群馬大學(xué)大學(xué)院理工學(xué)府， 日本 桐生 376-8515)

滑坡、泥石流及冰川等高速運(yùn)動(dòng)體入水沖擊水庫(kù)、湖泊、河道及海灣等封閉水域內(nèi)的水體時(shí)會(huì)誘發(fā)涌浪，并由此產(chǎn)生大規(guī)模的災(zāi)害。其代表性的實(shí)例為1963年意大利瓦依昂水庫(kù)庫(kù)岸發(fā)生大面積整體巖質(zhì)滑坡，滑坡體長(zhǎng)2 km，寬約1.6 km，滑坡體積達(dá)2.4億m3，下滑運(yùn)動(dòng)速度達(dá)15～30 m/s的滑坡體突入水庫(kù)時(shí)產(chǎn)生高達(dá)90～130 m的涌浪，給庫(kù)岸及大壩下游集鎮(zhèn)造成毀滅性破壞[1]。在日本2011年7月6號(hào)臺(tái)風(fēng)時(shí)奈半利川平鍋大壩上游流入大量的泥石流，引發(fā)超過(guò)5 m高的涌浪，給大壩閘門(mén)造成損傷并產(chǎn)生漫頂，同時(shí)損壞了大壩上游的平鍋吊橋。紫坪鋪大壩為高158 m，庫(kù)容量達(dá)11.1億m3的混泥土面板堆石壩，它與2008年5月12日發(fā)生的汶川地震震中距離只有17 km。汶川地震引發(fā)了紫坪鋪滑坡，滑體中心到庫(kù)水面的高差達(dá)700 m，約45萬(wàn)m3的滑坡巖體高速?zèng)_入水庫(kù)，產(chǎn)生高達(dá)25 m的涌浪，大約70名湖邊垂釣人被涌浪卷入水庫(kù)造成死亡，10臺(tái)汽車(chē)也被卷入湖中[2]。表1列出了一些典型的涌浪及其造成的人員傷亡情況。

表1 歷史性的滑坡或部分水下滑坡產(chǎn)生的涌浪

為防止涌浪災(zāi)害需要能確立涌浪特性及其規(guī)模等的預(yù)測(cè)方法。淺水長(zhǎng)波方程常用作涌浪的控制方程。因?yàn)闇\水長(zhǎng)波方程中包含有移流項(xiàng)，數(shù)值計(jì)算的穩(wěn)定性很差，有必要解決數(shù)值計(jì)算的穩(wěn)定問(wèn)題。另外，有限元離散得到的總系數(shù)矩陣是非對(duì)稱(chēng)的，因此大規(guī)模涌浪分析時(shí)為縮短計(jì)算時(shí)間，需要開(kāi)發(fā)線性方程組的多核多線程并行解法。本研究中，利用Intel MKL提供的共享內(nèi)存機(jī)器上實(shí)現(xiàn)的稀疏線性方程求解器PARDISO，開(kāi)發(fā)了涌浪高性能有限元分析程序，并用一些理論解和實(shí)驗(yàn)結(jié)果驗(yàn)證了開(kāi)發(fā)的程序。

1 涌浪有限元分析

1.1 控制方程

用豎直方向平均流速表示的淺水長(zhǎng)波方程包括連續(xù)方程(1)、x和y方向的運(yùn)動(dòng)方程(2)和(3)[3-4]。

(1)

(2)

(3)

式中：η為水面標(biāo)高；z為地表面標(biāo)高(見(jiàn)圖1)；水深h=η-z；x方向流束U=hu；y方向流束V=hv；g為重力加速度；n為曼寧系數(shù)；ε為水深平均渦動(dòng)黏性系數(shù)。為考慮滑坡運(yùn)動(dòng)引起的地表面標(biāo)高的變化，與一般的淺水長(zhǎng)波方程有所不同，式(1)中同時(shí)引入了水面和地表面標(biāo)高。

圖1 涌浪示意圖

1.2 穩(wěn)定化有限元離散

用有限元對(duì)上述3個(gè)控制方程進(jìn)行有限元離散，并引入SUPG(Streamline upwind Petrov Galerkin)項(xiàng)[5-6](其物理意義見(jiàn)圖2)和Shock capturing項(xiàng)，則連續(xù)方程可表示為

(4)

x方向的運(yùn)動(dòng)方程可以表示為

(5)

(6)

圖2 權(quán)重函數(shù)(左)Galerkin有限元法和(右)SUPG有限元法

同理可得y方向的運(yùn)動(dòng)方程。

上述方程如用矩陣表示，則連續(xù)方程表示為

(7)

x和y方向的運(yùn)動(dòng)方程可分別表示為

(8)

(9)

式中，η、U、V分別為各節(jié)點(diǎn)水面標(biāo)高x和y方向流束組成的向量。另外，

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(24)

(25)

另外，τSUPG為一維穩(wěn)定化參數(shù)[7]擴(kuò)展到多維條件下的穩(wěn)定化參數(shù)[8-9]。

(26)

(27)

(28)

(29)

單元尺寸可按下式計(jì)算[10]

(30)

Shock capturing系數(shù)可按下式定義[11]

νSHOC=τSHOC(uint)2

(31)

(32)

1.3 時(shí)間的直接積分

考慮時(shí)間間隔為Δt(Δt=tn+1-tn)的兩個(gè)狀態(tài)tn+1和tn時(shí)刻的水面標(biāo)高、x和y方向流束向量分別為ηn+1、Un+1、Vn+1和ηn、Un、Vn，則水面標(biāo)高、x和y方向流束向量的時(shí)間微分可表示為：

(33)

(34)

(35)

tn+θh時(shí)刻的水面標(biāo)高、tn+θu時(shí)刻x和y方向流束向量的時(shí)間微分可表示為：

η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為時(shí)間積分參數(shù)，θh和θu≥1/2時(shí)時(shí)間積分是無(wú)條件穩(wěn)定的。本文取θh=θu=1/2，即用Crank-Nicolson法進(jìn)行時(shí)間積分。將式(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分別為時(shí)步n-1和n的時(shí)步長(zhǎng)。當(dāng)這兩個(gè)時(shí)步長(zhǎng)相等時(shí)，式(50)可簡(jiǎn)化成等時(shí)步Adams-Bashforth法。

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

2.1 求解器PARDISO簡(jiǎn)介

空間與時(shí)間離散得到的大規(guī)模稀疏矩陣作為系數(shù)的一次聯(lián)立線性方程組Ax=b的求解為涌浪有限元分析的中心任務(wù)之一，因此采用適合于稀疏線性方程組的高速魯棒的求解器就非常重要。直接法對(duì)迭代法不收斂的問(wèn)題也可能求解，特別是矩陣對(duì)稱(chēng)正定時(shí)，如果沒(méi)有數(shù)值誤差則一定可以求解。因?yàn)檫@些特長(zhǎng)，直接法得到了廣泛的應(yīng)用。通過(guò)考慮矩陣變帶寬一維存儲(chǔ)內(nèi)部的零元素，可進(jìn)一步減少計(jì)算量和節(jié)約內(nèi)存的稀疏線性方程組直接求解法現(xiàn)在已成為主流。包括PARDISO[12]在內(nèi)的稀疏線性方程組直接法求解器的計(jì)算步驟(圖3)一般包括重排序、符號(hào)分解、LU分解、前進(jìn)后退代入等4步的順次計(jì)算。以下簡(jiǎn)單介紹一下各步的計(jì)算。

Ⅰ：矩陣A的非零元素位置發(fā)生變化時(shí) Ⅱ：矩陣A的非零元素位置不變但值發(fā)生變化時(shí) Ⅲ：只有b發(fā)生變化時(shí)

重排序就是用合適的置換矩陣P，使LU分解PAPT時(shí)填入元素盡可能的少。這里填入元素是指LU分解前為零但分解后為非零的元素。重排序的方法有等幅縮小(最小次數(shù)法，Reverse Cuthill-McKee法)，三角化(Markowitz法，Tewarson法)，分塊化(Stewart法，Nested Dissection法)等許多算法[13]。PARDISO用METIS算法包中最小次數(shù)算法或者Nested Dissection算法進(jìn)行重排序。

符號(hào)分解并不進(jìn)行具體元素的分解計(jì)算，而只著眼于矩陣A的LU分解時(shí)非零元素的分布形式，找出LU分解后的非零元素的位置。由此算出LU分解所需要的內(nèi)存量和計(jì)算量，確保LU分解后保存非零元素所需的內(nèi)存，并記錄非零元素的位置。為高效進(jìn)行符號(hào)分解，利用列消去樹(shù)的概念，把問(wèn)題歸為有效圖的路徑探索問(wèn)題，從而使高速計(jì)算成為可能。

用符號(hào)分解時(shí)確保的內(nèi)存量進(jìn)行實(shí)際的LU分解。作為主要的LU分解方法，參照更新列的右側(cè)的right-looking算法，左側(cè)的left-looking算法廣為人知。本程序所用的大規(guī)模稀疏線性方程組的多核多線程并行求解器PARDISO組合利用了left-looking和right-looking算法以實(shí)現(xiàn)高效并行計(jì)算。

用LU分解后的下三角形矩陣L進(jìn)行前進(jìn)代入，用上三角形矩陣U進(jìn)行后退代入可求得解x。如果反復(fù)求解系數(shù)矩陣A相同的一次聯(lián)立線性方程組(例如，一個(gè)加載步用修正Newton-Raphson法進(jìn)行計(jì)算時(shí))，則只要迭代進(jìn)行前進(jìn)后退代入計(jì)算就可求解。系數(shù)矩陣的非零元素的構(gòu)造相同時(shí)則必須返回到LU分解進(jìn)行計(jì)算。邊界條件或分析范圍等變化時(shí)則必須重新合成系數(shù)矩陣并從重排序開(kāi)始求解過(guò)程。

2.2 超級(jí)節(jié)點(diǎn)

所謂超級(jí)節(jié)點(diǎn)為上三角矩陣L中全為非零，各列有相同的非零構(gòu)造的列的集合。圖4所示{1, 2}，{3, 4}，{5}，{6, 7, 8}分別為階數(shù)2，2，1，3的超級(jí)節(jié)點(diǎn)。left-looking算法引入超級(jí)節(jié)點(diǎn)后，可促進(jìn)分塊化，也就是數(shù)據(jù)訪問(wèn)的局部化，可大幅提高階層構(gòu)造內(nèi)存計(jì)算機(jī)的計(jì)算速度。把不同的非零構(gòu)造的零元素假定成非零元素，生成超級(jí)節(jié)點(diǎn)，在有些情況下可以得到更有效的超級(jí)節(jié)點(diǎn)。

多核多線程并行計(jì)算稀疏線性方程組求解器PARDISO在消去樹(shù)的并行化、節(jié)點(diǎn)水平的并行化、數(shù)據(jù)通道處理的并行化3個(gè)水平上實(shí)現(xiàn)了并行化。

圖4 稀疏矩陣A(左)、L和U的非零構(gòu)造及超級(jí)節(jié)點(diǎn)(右)

PARDISO用于涌浪有限元計(jì)算時(shí)需要生成線性方程組系數(shù)矩陣A的非零元素的位置列表。PARDISO使用以行為主的存儲(chǔ)方法，即變形CSR(compressed sparse row)形式，對(duì)稱(chēng)矩陣只存儲(chǔ)上半三角元素。該方法以行為單位存儲(chǔ)每個(gè)非零數(shù)據(jù)。PARDISO對(duì)一個(gè)稀疏矩陣A的存儲(chǔ)包括了3個(gè)數(shù)組：

1)values-矩陣A的非零元素。矩陣A的非零數(shù)據(jù)通過(guò)下面的columns與rowindex映射到values數(shù)組中。

2)columns-values中每個(gè)元素所在矩陣的列。

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

具體可參照參考文獻(xiàn)[14]。

3 計(jì)算實(shí)例

用研發(fā)的軟件首先計(jì)算了圖5所示長(zhǎng)20 m的水槽中水位差0.8 m的涌浪傳播。這個(gè)計(jì)算實(shí)例為靜止的水壩瞬時(shí)坍塌放流的現(xiàn)象。為與完全流體的理論解進(jìn)行比較，計(jì)算中假定水深平均渦動(dòng)黏性系數(shù)ε=0，另外水槽側(cè)壁假定為slip條件。圖5所示為水壩坍塌1 s后的水深和流速的數(shù)值計(jì)算結(jié)果和理論解的比較。圖5表明計(jì)算結(jié)果和理論解非常一致。

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

為確認(rèn)移動(dòng)邊界問(wèn)題的計(jì)算精度，計(jì)算了水槽長(zhǎng)10 m，中心左側(cè)水位高0.2 m，右側(cè)為無(wú)水干床水槽中的水壩瞬時(shí)坍塌問(wèn)題。同前例，計(jì)算中假定水深平均渦動(dòng)黏性系數(shù)ε=0，水槽側(cè)壁為slip條件。圖6所示水壩坍塌1 s后水深計(jì)算結(jié)果和理論解的比較。由圖6可知干床水槽條件下計(jì)算結(jié)果和理論解也非常一致。

圖7所示為長(zhǎng)7.2 m水深0.2 m的水槽內(nèi)長(zhǎng)2.44 m的底面活塞向上運(yùn)動(dòng)形成涌浪的實(shí)驗(yàn)[15]。活塞向上運(yùn)動(dòng)速度ζ(t)由下式表示

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

(51)

式中：ζ為活塞運(yùn)動(dòng)位置；ζ0與α為活塞運(yùn)動(dòng)控制參數(shù)，模型試驗(yàn)中ζ0=0.2h，α=0.235 8。

圖6 干床水槽中水壩坍塌1 s后水深的計(jì)算結(jié)果與理論解的比較

圖7 水槽及造波活塞

計(jì)算條件為水深平均渦動(dòng)黏性系數(shù)ε=0，曼寧系數(shù)n=0，時(shí)步長(zhǎng)Δt=0.014 29(無(wú)量綱表示)，采用0.2 m均勻網(wǎng)格。計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果的比較見(jiàn)圖8，圖8(上)為活塞左端，(下)為活塞右端的水位隨時(shí)間的變化過(guò)程。圖8表明計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果相當(dāng)一致。

圖8 活塞造波實(shí)驗(yàn)結(jié)果(左圖實(shí)線)與計(jì)算結(jié)果(右圖)的比較

4 結(jié) 論

高速運(yùn)動(dòng)的滑坡、泥石流及冰川等入水沖擊水庫(kù)、湖泊、河道及海灣等封閉水域內(nèi)的水體時(shí)產(chǎn)生涌浪的預(yù)測(cè)經(jīng)常是個(gè)大規(guī)模的計(jì)算問(wèn)題。為縮短大規(guī)模涌浪分析的計(jì)算時(shí)間，本研究利用Intel MKL提供的共享內(nèi)存機(jī)器上實(shí)現(xiàn)的稀疏線性方程組求解器PARDISO，開(kāi)發(fā)了涌浪高性能有限元分析程序，并用一些理論解和實(shí)驗(yàn)結(jié)果驗(yàn)證了開(kāi)發(fā)的程序。垂直水壩潰壩的理論解和數(shù)值計(jì)算結(jié)果，以及室內(nèi)實(shí)驗(yàn)結(jié)果和數(shù)值計(jì)算結(jié)果的比較表明研發(fā)的涌浪高性能有限元軟件可以快速得到可靠的計(jì)算結(jié)果。