葉健，陶和平，陳錦雄，陳曉清

(中國(guó)科學(xué)院 成都山地災(zāi)害與環(huán)境研究所，四川 成都，610041)

我國(guó)幅員遼闊，地貌差異顯著，地質(zhì)環(huán)境復(fù)雜，災(zāi)害性氣象及氣候條件出現(xiàn)頻繁，再加上不斷擴(kuò)大的人類活動(dòng)因素，使得山地災(zāi)害頻繁發(fā)生，嚴(yán)重威脅山區(qū)建設(shè)和人民生命財(cái)產(chǎn)安全。其中，高速遠(yuǎn)程滑坡-碎屑流是一種特殊的地質(zhì)災(zāi)害，它是指高速遠(yuǎn)程滑坡或崩塌在運(yùn)動(dòng)過(guò)程中轉(zhuǎn)化而成的碎屑流體[1]。目前，很多礦山、采石場(chǎng)等將廢棄物放于附近的山坡上，當(dāng)遇到地震或附近的爆破作業(yè)等情況時(shí)，或者當(dāng)滑坡體與山體碰撞粉碎后，散體堆積物就會(huì)沿坡滑下[2-4]，形成碎屑流。形成后的碎屑流以高速度和大的位移引發(fā)災(zāi)難性事故，往往會(huì)造成嚴(yán)重的生命財(cái)產(chǎn)損失[1]。2008年5月12日汶川地震后，文家溝頂部距離汶川地震發(fā)震斷層僅3.6 km處，山體在長(zhǎng)時(shí)間強(qiáng)烈地震震動(dòng)作用下失穩(wěn)而形成滑坡?；略诟咚傧禄^(guò)程中轉(zhuǎn)化為碎屑流[5]。2009年5月發(fā)生在我國(guó)重慶武隆的雞尾山滑坡-碎屑流，運(yùn)行了1.5 km，等值摩擦因數(shù)僅為0.27，高速運(yùn)行的碎屑流掩埋了下游的鐵礦，造成 10人死亡，64人失蹤[6]。由于碎屑流災(zāi)害近年來(lái)的頻繁發(fā)生，相關(guān)部門已著手研究碎屑流防災(zāi)減災(zāi)對(duì)策。在碎屑流防災(zāi)減災(zāi)的措施中，工程治理是防治巖石碎屑流災(zāi)害的主要措施之一，應(yīng)用極其廣泛，對(duì)消除和減輕地質(zhì)災(zāi)害對(duì)人民生命財(cái)產(chǎn)的危害發(fā)揮了重要的作用。攔砂壩是工程防治中較常采用的手段，利用攔砂壩可抑制碎屑流的發(fā)生，消減碎屑流流量，延緩碎屑流到達(dá)時(shí)間，減少碎屑流土砂總流出量，阻滯碎屑流顆粒輸移、減勢(shì)、穩(wěn)定溝床及減輕兩岸侵蝕[7]。在進(jìn)行碎屑流工程防治的過(guò)程中，研究碎屑流運(yùn)動(dòng)規(guī)律，逐步建立和完善碎屑流模型，并根據(jù)碎屑流災(zāi)害防治工程設(shè)計(jì)規(guī)范建立攔沙壩模型，使其能夠以三維可視化的方式再現(xiàn)碎屑流運(yùn)動(dòng)和堆積的過(guò)程，將有助于人們加深對(duì)碎屑流災(zāi)害的認(rèn)識(shí)，并有助于進(jìn)行災(zāi)害的防災(zāi)減災(zāi)措施。盡管目前在地質(zhì)學(xué)領(lǐng)域里，已經(jīng)有不少學(xué)者利用離散元法對(duì)碎屑流進(jìn)行數(shù)值模擬，但把離散元法運(yùn)用到三維仿真領(lǐng)域并進(jìn)行碎屑流與攔砂壩、排導(dǎo)槽交互的可視化工作卻較少。因此，本研究旨在將三維可視化理論與技術(shù)引入到碎屑流現(xiàn)象及其防災(zāi)減災(zāi)的三維建模與可視化中，建立碎屑流模型，不僅可以利用傳統(tǒng)的數(shù)值模擬技術(shù)，同時(shí)可以借助動(dòng)態(tài)、三維交互的三維仿真技術(shù)來(lái)模擬碎屑流現(xiàn)象及其碎屑流的防災(zāi)減災(zāi)過(guò)程。在利用計(jì)算機(jī)模擬離散物體現(xiàn)象的方法中，基于散體理論建立的離散元(discrete element method, DEM)法具有不依賴網(wǎng)格和描述破壞問題便利的特點(diǎn)，能夠更好地反映一些自然現(xiàn)象發(fā)生過(guò)程的本質(zhì)，它是散體力學(xué)分析的一種有效的工具[8-9]，該算法非常適合于進(jìn)行碎屑流現(xiàn)象的研究。而另一方面，隨著圖形技術(shù)的飛速發(fā)展,近年來(lái)傳統(tǒng)的圖形硬件已經(jīng)演化為可編程圖形處理單元(GPU)，其強(qiáng)大的計(jì)算性能、靈活的可編程性以及面向個(gè)人消費(fèi)市場(chǎng)的低廉價(jià)格，吸引了越來(lái)越多的研究者將 GPU 應(yīng)用于解決圖形渲染之外的通用計(jì)算(GPGPU)[10-11]，并利用GPU高效的并行計(jì)算來(lái)進(jìn)行復(fù)雜的物理計(jì)算。本研究基于散體理論的離散元法建立的碎屑流模型，呈現(xiàn)出高度的數(shù)據(jù)并行性，因此，可以充分利用GPU的高度并行性來(lái)進(jìn)行碎屑流的建模與三維可視化的實(shí)時(shí)渲染。可視化場(chǎng)景中攔砂壩和排導(dǎo)槽的建模及兩者與碎屑流之間邊界問題的正確處理在本研究中較重要，直接影響到巖石碎屑流與攔砂壩和排導(dǎo)槽交互場(chǎng)景的可視化表達(dá)。為此，本文作者首先分析碎屑流在場(chǎng)景中運(yùn)動(dòng)的整個(gè)過(guò)程；然后，構(gòu)建攔砂壩和排導(dǎo)槽模型，并確立它們的邊界；最后，根據(jù)碎屑流顆粒與攔砂壩和排導(dǎo)槽邊界之間的受力分析，建立有效的碎屑流顆粒與攔砂壩和排導(dǎo)槽邊界問題的處理方法，并充分利用GPU的高度并行性和可編程性，實(shí)現(xiàn)巖石碎屑流與攔砂壩和排導(dǎo)槽交互場(chǎng)景的三維建模與可視化。

1 基于離散元法的巖石碎屑流建模

1.1 巖石碎屑流運(yùn)動(dòng)特點(diǎn)

巖石碎屑流常常是由高速滑坡的發(fā)生而引起的，但巖石碎屑流與滑坡不同，不是沿滑動(dòng)面的整體，而是呈顆粒介質(zhì)的運(yùn)動(dòng)，液體對(duì)運(yùn)動(dòng)不起重要作用[12]。從啟動(dòng)到靜止，高速遠(yuǎn)程滑坡或崩滑在運(yùn)動(dòng)過(guò)程中逐漸轉(zhuǎn)化為巖石碎屑流，最終以高速遠(yuǎn)程滑坡-碎屑流的復(fù)合形式出現(xiàn)[13]。

巖石碎屑流發(fā)生的整個(gè)過(guò)程通常可歸納為啟程階段、加速階段和堆積階段。

第一階段(啟程階段)：當(dāng)遇到地震或附近的爆破作業(yè)等情況時(shí)，或者當(dāng)滑坡體與山體碰撞粉碎后，散體堆積物就會(huì)沿坡滑下。

第二階段(加速階段)：高速巖石崩滑體在經(jīng)歷了啟程階段后，與運(yùn)動(dòng)路徑兩側(cè)或前方突出的阻擋山體發(fā)生強(qiáng)烈的沖擊碰撞崩解形成巖石碎屑流，并以這種形態(tài)進(jìn)入加速階段，此時(shí)，巖石碎屑流勢(shì)能不斷釋放，一部分轉(zhuǎn)化為巖石碎屑流動(dòng)能，一部分轉(zhuǎn)化為摩擦熱能。

第三階段(堆積階段)：巖石碎屑流的運(yùn)動(dòng)隨著山體的坡角減小，阻力變大，這時(shí)隨著動(dòng)能慢慢變小，巖石碎屑流運(yùn)動(dòng)速度也越來(lái)越小，并開始堆積，直至所有碎屑物質(zhì)堆積完畢為止。

巖石碎屑流在運(yùn)動(dòng)過(guò)程中常常呈現(xiàn)“流態(tài)化”的特征，這是它能高速運(yùn)動(dòng)的主要原因之一[14-18]。在巖石碎屑流的啟動(dòng)、加速和堆積的整個(gè)過(guò)程中，巖石碎屑流不會(huì)產(chǎn)生大體積的“飛行”運(yùn)動(dòng)，且?guī)r石碎屑流空隙度較大，具有離散性、流動(dòng)性和高運(yùn)動(dòng)性的特點(diǎn)。

1.2 基于離散元法的巖石碎屑流模型

離散元法最早是由Cundall于1971年提出的一種不連續(xù)的數(shù)值方法。其基本思想是把不連續(xù)體分離為剛性元素的集合，使各個(gè)剛性元素滿足運(yùn)動(dòng)方程，用時(shí)步迭代的方法求解各剛性元素的運(yùn)動(dòng)方程，繼而求得不連續(xù)體的整體運(yùn)動(dòng)形態(tài)[19-20]。離散元法自問世以來(lái)，在巖土工程和顆粒散體工程這兩大傳統(tǒng)的應(yīng)用領(lǐng)域中發(fā)揮了其他數(shù)值算法不可替代的作用，特別是在計(jì)算和模擬材料具體的破壞形式和破壞的整個(gè)過(guò)程方面，離散元法顯示出明顯的優(yōu)勢(shì)[21-22]。離散元法能模擬巖石顆粒之間存在的滑動(dòng)、平移和轉(zhuǎn)動(dòng)等復(fù)雜過(guò)程，具有宏觀上的不連續(xù)性，可以較真實(shí)地模擬碎屑流顆粒的運(yùn)動(dòng)，也可以模擬高速大位移運(yùn)動(dòng)。由于高速遠(yuǎn)程滑坡后所形成的巖石碎屑流具有離散度大、流動(dòng)顯著且高速大位移運(yùn)動(dòng)等特點(diǎn)，因此，適合于利用離散元法進(jìn)行模擬。

離散元法是基于牛頓第二運(yùn)動(dòng)定理和顆粒在接觸點(diǎn)上的力與位移之間的本構(gòu)關(guān)系所建立的數(shù)值模擬方法[23-24]。在顆粒運(yùn)動(dòng)的過(guò)程中，首先由外邊界條件決定顆粒所受外力，然后根據(jù)顆粒之間的位置關(guān)系來(lái)確定顆粒間的相互作用，最后求出顆粒所受的合力，根據(jù)這個(gè)合力求出顆粒的加速度和速度，進(jìn)而求出各個(gè)顆粒的位移。

當(dāng)2個(gè)顆粒之間發(fā)生碰撞，顆粒在前進(jìn)方向受到阻力，此時(shí)顆粒之間會(huì)產(chǎn)生接觸力Fc，并且顆粒在接觸點(diǎn)處發(fā)生變形，當(dāng)顆粒達(dá)到最大位移變形時(shí)，顆粒處于靜止?fàn)顟B(tài)，然后沿運(yùn)動(dòng)的方向反彈。若此時(shí)顆粒間的碰撞為偏心碰撞，則相撞點(diǎn)處的接觸力可以分解為垂直于接觸面的法向接觸力Fcn和平行于接觸面的切向接觸力Fcs。顆粒之間的接觸模型是離散元法的核心，本模型在顆粒之間施加了變形元件彈簧和阻尼器，采用的是彈簧-阻尼器模型[25]，模型中顆粒j和顆粒i之間的法向接觸力Fcn可以表示為：

其中：fi,s為顆粒i和顆粒j碰撞的排斥力，顆粒在碰撞過(guò)程中，顆粒阻尼力fi,d使用阻尼器模擬；κ，η，d，ri,j和vi,j分別為彈簧彈性系數(shù)、阻尼系數(shù)、顆粒直徑，相對(duì)位置和顆粒j對(duì)顆粒i的相對(duì)速度。

碰撞后，顆粒的動(dòng)能會(huì)產(chǎn)生一定的損失，損失的大小與顆粒的彈性阻尼系數(shù)及顆粒間的相對(duì)速度有關(guān)。而切向接觸力Fcs與相對(duì)切向速度vij,t可描述為：

相對(duì)切向速度vij,t的計(jì)算公式為：

作用于剛體的力Fc和力矩Tc為施加于所用剛體顆粒之上的力和力矩的和：

其中：ri為當(dāng)前顆粒i質(zhì)心的相對(duì)位置。

1.3 離散元法在GPU通用計(jì)算平臺(tái)上的實(shí)現(xiàn)

為了更好地模擬巖石碎屑流的細(xì)節(jié)，碎屑流顆粒的數(shù)目往往達(dá)到數(shù)十萬(wàn)個(gè)，若此時(shí)將碎屑流顆粒間相互作用的計(jì)算完全放在CPU上，碎屑流運(yùn)動(dòng)的模擬幾乎不能滿足實(shí)時(shí)性和交互性的要求，而利用離散元法模擬碎屑流可以實(shí)現(xiàn)碎屑流顆粒之間的交互計(jì)算，非常適合在GPU上實(shí)現(xiàn)。在GPU通用計(jì)算平臺(tái)上，由于將每個(gè)顆粒的計(jì)算都分配到每個(gè)GPU線程上執(zhí)行，這樣就可以充分利用GPU的并行性，可極大地加快算法的計(jì)算速度，從而滿足數(shù)十萬(wàn)個(gè)碎屑流顆粒相互作用的模擬場(chǎng)景實(shí)時(shí)性和交互性的要求。

在本研究所構(gòu)建的碎屑流模型中，碎屑流顆粒間相互作用力的計(jì)算需要建立高效的顆粒查找和排序機(jī)制以及顆粒間碰撞檢測(cè)機(jī)制，因此，首先需要查詢鄰近顆粒并計(jì)算它們之間的碰撞。為了避免對(duì)鄰近顆粒進(jìn)行n2次搜索，采用鄰近顆粒搜索算法，該算法與文獻(xiàn)[26]中的相似，可以概括為：

(1) 劃分仿真區(qū)域?yàn)榻y(tǒng)一網(wǎng)格；

(2) 利用每一顆粒的空間位置去查找它所屬的單元；

(3) 把顆粒單元位置作為哈希函數(shù)的輸入；

(4) 根據(jù)空間哈希排序顆粒；

(5) 根據(jù)它們的哈希值在線性緩沖區(qū)中重新排序顆粒。

因此，把三維空間劃分為統(tǒng)一均勻的網(wǎng)格單元，每一顆粒被賦予一個(gè)基于中心點(diǎn)的網(wǎng)格單元，顆粒的索引被存儲(chǔ)在這個(gè)網(wǎng)格單元中，并以該網(wǎng)格單元索引作為該顆粒的哈希值。一旦建立了網(wǎng)格結(jié)構(gòu)，就能使用它來(lái)加速顆粒和顆粒之間的交互計(jì)算。

當(dāng)計(jì)算顆粒所在的網(wǎng)格單元時(shí)，可以循環(huán)遍歷相鄰的26個(gè)網(wǎng)格單元(3×3×3個(gè)單元)并檢測(cè)是否與這26個(gè)網(wǎng)格單元中的顆粒之間發(fā)生碰撞：若發(fā)生碰撞，則通過(guò)計(jì)算便可得到該顆粒所在的網(wǎng)格單元，繼而寫入該網(wǎng)格單元索引到正確位置的網(wǎng)格數(shù)組中，而此時(shí)顆粒和網(wǎng)格的對(duì)應(yīng)關(guān)系已經(jīng)發(fā)生了變化。由于更新顆粒位置后的顆粒-網(wǎng)格對(duì)應(yīng)關(guān)系是亂序的，必須通過(guò)排序算法以網(wǎng)格為序排列起來(lái)才能在計(jì)算中快速地取得鄰近顆粒的信息。顆粒的排序是基于哈希值來(lái)排序顆粒的，采用的排序算法為并行基數(shù)排序[27]。利用并行基數(shù)算法可以得到顆粒和網(wǎng)格單元的正確對(duì)應(yīng)關(guān)系，從而保證在計(jì)算中可以快速地取得鄰近顆粒的信息。由于本研究將離散元法的全部模擬計(jì)算都分配到GPU中處理，因此，充分利用了GPU的高并行性和可編程性，可利用離散元法實(shí)現(xiàn)碎屑流的實(shí)時(shí)計(jì)算和模擬。

2 攔砂壩和排導(dǎo)槽的建模

工程治理是防治巖石碎屑流災(zāi)害的主要措施之一，應(yīng)用極其廣泛，在消除和減輕地質(zhì)災(zāi)害對(duì)人民生命財(cái)產(chǎn)的危害中發(fā)揮了重要的作用。修建于溝道中的攔砂壩是工程治理的重要手段，通過(guò)在巖石碎屑流運(yùn)動(dòng)的通道中構(gòu)建層層攔截的攔砂壩，以攔砂壩的規(guī)模效應(yīng)形成了較強(qiáng)的攔沙、削峰、滯洪能力和上攔下保的作用。不僅如此，攔砂壩還具有固定河床，抬高河底高程，穩(wěn)定邊坡，遏制溝岸擴(kuò)張、溝底下切和溝頭前進(jìn)，減輕溝道侵蝕，大大減少泥沙下泄量的作用[28]。而排導(dǎo)槽由于結(jié)構(gòu)簡(jiǎn)單、就地取材、施工方便和效果顯著等優(yōu)點(diǎn)得到廣泛應(yīng)用[29]。

由于攔砂壩和排導(dǎo)槽在山地災(zāi)害工程治理中發(fā)揮了重要的作用，本研究將以巖石碎屑流作為研究對(duì)象，通過(guò)利用離散元法來(lái)建立巖石碎屑流模型，然后構(gòu)建攔砂壩和排導(dǎo)槽兩大巖石碎屑流防護(hù)措施模型，重點(diǎn)模擬巖石碎屑流形成以后,巖石碎屑流蓄滿在第一級(jí)攔砂壩后被第二級(jí)攔砂壩攔截的防災(zāi)減災(zāi)過(guò)程。

為了構(gòu)建碎屑流與攔沙壩、排導(dǎo)槽交互的虛擬場(chǎng)景，需要構(gòu)建攔砂壩和排導(dǎo)槽的防護(hù)措施模型。要構(gòu)建該防護(hù)措施模型，其關(guān)鍵在于如何確定和處理攔砂壩及其排導(dǎo)槽的邊界。為了更好地再現(xiàn)虛擬場(chǎng)景中攔砂壩的防災(zāi)減災(zāi)過(guò)程，在虛擬場(chǎng)景中構(gòu)建了二級(jí)巖石碎屑流攔砂壩，當(dāng)巖石碎屑流發(fā)生時(shí)，攔砂壩能將巖石碎屑流攔截并堆積在第一級(jí)攔砂壩庫(kù)內(nèi)，直到該級(jí)攔砂壩蓄滿。若大量巖石碎屑流從壩的上方溢出，則溢出的巖石碎屑流開始在第二級(jí)攔砂壩內(nèi)蓄積，由于溢出的巖石碎屑流經(jīng)過(guò)上一級(jí)攔砂壩的攔擋，其流量和規(guī)模將大幅度減少，流速也減小。經(jīng)過(guò)二級(jí)攔砂壩攔蓄與消能后，其危害可以大大減小。

2.1 邊界的確定

巖石碎屑流模型建立后，需要確定攔砂壩和排導(dǎo)槽的邊界。其邊界的俯視圖如圖1所示。

圖1 攔沙壩和排導(dǎo)槽的邊界平面圖Fig.1 Boundary planform of debris dam and drainage canal

圖2 攔沙壩與排導(dǎo)槽的剖面示意圖Fig.2 The diagrammatic crosssection of debris dam and drainage canal

圖2中排導(dǎo)槽的斜坡段L1的坡度α為30°，L2段為排導(dǎo)槽的水平段，巖石碎屑流自排導(dǎo)槽頂端在重力的作用下沖向排導(dǎo)槽底端，巖石碎屑流首先要觸碰排導(dǎo)槽L1段的底端和兩側(cè)，然后被攔砂壩攔截并堆積在第一級(jí)攔砂壩A的庫(kù)內(nèi)，直到該級(jí)攔砂壩蓄滿為止；然后，巖石碎屑流開始從攔砂壩A溢出并蓄積在攔砂壩B內(nèi)，經(jīng)過(guò)兩級(jí)攔砂壩的攔截，巖石碎屑流基本被攔截在攔砂壩內(nèi)，只有少量巖石碎屑流顆粒會(huì)落入排導(dǎo)槽L2段。

結(jié)合圖1和圖2可知：首先需要確定L1和L2段排導(dǎo)槽底端及其兩側(cè)的邊界，然后，必須確定攔砂壩A和攔砂壩B的邊界。根據(jù)以上分析，在場(chǎng)景中構(gòu)建了攔砂壩和排導(dǎo)槽的模型，其效果圖如圖3所示。

圖3 攔砂壩和排導(dǎo)槽的三維建模效果Fig.3 Effect of three-dimensional modeling of debris dam and drainage canal

2.2 邊界的處理

在確定攔砂壩和排導(dǎo)槽的邊界后，需要處理碎屑流顆粒與邊界的碰撞。碎屑流顆粒與邊界的碰撞過(guò)程如圖4所示，顆粒首先從邊界面的某一方向撞擊邊界面，在碰撞的過(guò)程中，當(dāng)顆粒碰撞邊界面的那一時(shí)刻，顆粒發(fā)生形變，在回彈的那一瞬間，由于能量損失，反彈后的速度減小。

假設(shè)顆粒與邊界在碰撞前一瞬間的速度為v，這時(shí)，可以把顆粒沿著邊界面的法向和切向分解為vn和vt。其顆粒與邊界碰撞前和碰撞后速度分解示意圖如圖5所示，碰撞后的法向和切向速度將變?yōu)椋?/p>

圖4 巖石碎屑流顆粒與邊界的碰撞過(guò)程示意圖Fig.4 Sketch map of collision process of rock avalanche granule and boundary

圖5 巖石碎屑流顆粒與邊界碰撞前后Fig.5 Rock avalanche granule before and after colliding with boundary

其中：α為法向速度修改系數(shù)，α=1時(shí)，顆粒的碰撞無(wú)能量損失；α=0時(shí)，顆粒的法向分量速度為 0。本實(shí)驗(yàn)α取為0.6，碰撞后法向速度反向。β為表面摩擦因數(shù)，控制障礙物表面的滑移特性。當(dāng)β=1 時(shí)，顆粒在切向碰撞時(shí)，無(wú)摩擦；當(dāng)β=0 時(shí)，顆粒沿邊界面的切向速度將變?yōu)?。本實(shí)驗(yàn)β取為0.75。為了進(jìn)行顆粒與邊界的碰撞檢測(cè)及顆粒碰撞后速度的計(jì)算，需在程序初始化之前保存攔砂壩和排導(dǎo)槽邊界距離場(chǎng)的信息，當(dāng)顆粒與邊界的距離d小于邊界碰撞閾值T時(shí)，可以確定顆粒與邊界發(fā)生碰撞，然后進(jìn)行受力和速度的計(jì)算。由于本研究?jī)H處理靜止的邊界，因此，當(dāng)檢測(cè)到碰撞且顆粒已進(jìn)入邊界的內(nèi)部時(shí)，邊界不會(huì)發(fā)生任何的移動(dòng)或變化，僅僅需要處理顆粒的速度。

3 實(shí)驗(yàn)結(jié)果

為了驗(yàn)證攔砂壩攔截巖石碎屑流的效果，把實(shí)驗(yàn)分為有攔砂壩攔截和無(wú)攔砂壩攔截2組。

本實(shí)驗(yàn)采用Intel(R) Core (TM)2 Duo 3.00 GHz，2.0 GB內(nèi)存，NVIDIA GeForce 8800GT 顯卡的PC機(jī)。實(shí)驗(yàn)軟件平臺(tái)為 Windows XP和 VS2005，并使用NVIDIA CUDA作為GPU計(jì)算架構(gòu)。

第1組實(shí)驗(yàn)如圖6所示，在重力作用下，巖石碎屑流從排導(dǎo)槽的頂部滑下，巖石碎屑流的勢(shì)能逐漸轉(zhuǎn)換為動(dòng)能，在沒有阻擋物阻擋的情況下，巖石碎屑流的速度越來(lái)越快。從實(shí)驗(yàn)結(jié)果可以看出：巖石碎屑流的運(yùn)動(dòng)將對(duì)排導(dǎo)槽的水平段造成巨大的沖擊。

圖6 巖石碎屑流在沒有攔砂壩攔截下的全過(guò)程Fig.6 Overall process of rock avalanche granule without debris dam

在第2組實(shí)驗(yàn)中，添加了攔砂壩，實(shí)驗(yàn)如圖7所示。盡管巖石碎屑流在滑落的過(guò)程中速度越來(lái)越快，但由于受到第一級(jí)攔砂壩的消能和攔截，部分巖石碎屑流將蓄積在第一級(jí)攔砂壩內(nèi)，而另一部分巖石碎屑流將繼續(xù)沖擊排導(dǎo)槽的水平段，在經(jīng)過(guò)第二級(jí)攔砂壩的消能和攔截后，巖石碎屑流將基本蓄積在攔砂壩內(nèi)，只有部分巖石碎屑流會(huì)飛濺到攔砂壩的水平段。從圖7(b)可以看出：巖石碎屑流與攔沙壩的碰撞伴隨著撞擊和爬高，其可視化效果能夠較真實(shí)地再現(xiàn)自然界中以下特點(diǎn)：當(dāng)巖石碎屑流勢(shì)能基本轉(zhuǎn)化為動(dòng)能后，巖石碎屑流運(yùn)動(dòng)速度快，慣性力大；當(dāng)前進(jìn)遇障礙阻擋時(shí)，可沖擊爬高、翻越障礙等。并再現(xiàn)了巖石碎屑流堆積的過(guò)程。另外，為了模擬出較豐富的細(xì)節(jié)效果，當(dāng)本模型初始化的粒子數(shù)達(dá)到122.2萬(wàn)個(gè)時(shí)，實(shí)時(shí)模擬速度仍然能達(dá)到10幀/s以上。

圖7 巖石碎屑流被攔砂壩攔截的全過(guò)程Fig.6 Overall process of rock avalanche granule intercepted by debris flow

4 結(jié)論

(1) 利用離散元法構(gòu)建了巖石碎屑流模型，并充分利用GPU的高度并行性特點(diǎn)，將巖石碎屑流模擬的大部分模擬放在GPU中計(jì)算，使得巖石碎屑流顆粒數(shù)為數(shù)百萬(wàn)時(shí)仍能滿足實(shí)時(shí)的要求。

(2) 在繪制巖石碎屑流與攔砂壩和排導(dǎo)槽交互的場(chǎng)景中，通過(guò)有效地確定和處理巖石碎屑流與攔砂壩和排導(dǎo)槽的邊界條件，模擬出巖石碎屑流蓄滿在第一級(jí)攔砂壩后被第二級(jí)攔砂壩攔截的過(guò)程，同時(shí)也展示了巖石碎屑流與攔沙壩和排導(dǎo)槽交互過(guò)程中撞擊和爬高的細(xì)節(jié)，并實(shí)現(xiàn)了巖石碎屑流堆積的過(guò)程。

(3) 為巖石碎屑流防災(zāi)工程的處理設(shè)計(jì)與施工提供了直觀的可視化分析平臺(tái)，此平臺(tái)可結(jié)合巖石碎屑流災(zāi)害防治工程設(shè)計(jì)規(guī)范和巖石碎屑流運(yùn)動(dòng)參數(shù)的計(jì)算結(jié)果設(shè)計(jì)出的巖石碎屑流防災(zāi)減災(zāi)工程，再現(xiàn)巖石碎屑流防災(zāi)減災(zāi)的過(guò)程。

在巖石碎屑流的動(dòng)態(tài)繪制中，尚未考慮巖石碎屑流顆粒級(jí)配的問題，所繪制的巖石碎屑流顆粒粒徑是均勻的，這與自然現(xiàn)象中的巖石碎屑流有所不同。為了提高巖石碎屑流模擬的效果，應(yīng)該考慮模擬不同顆粒半徑的巖石碎屑流。另外，應(yīng)充分發(fā)掘GPU更加廣闊的通用計(jì)算性能，利用能夠表達(dá)連續(xù)流體的物理模型繪制流體的動(dòng)態(tài)場(chǎng)景，并將模型應(yīng)用于工程設(shè)計(jì)和災(zāi)害評(píng)估。

[1]張明, 殷躍平, 吳樹仁, 等.高速遠(yuǎn)程滑坡-碎屑流運(yùn)動(dòng)機(jī)理研究發(fā)展現(xiàn)狀與展望[J].工程地質(zhì)學(xué)報(bào), 2010, 18(6)∶ 805-812.ZHANG Ming, YIN Yueping, WU Shuren, et al.Development status and prospects of studies on kinematics of long runout rock avalanches[J].Journal of Engineering Geology, 2010, 18(6)∶805-812.

[2]賀可強(qiáng), 安振遠(yuǎn).崩滑碎屑流的形成條件與形成類型[J].河北地質(zhì)學(xué)院學(xué)報(bào), 1996, 19(3/4)∶ 344-351.HE Keqiang, AN Zhenyuan.A preliminary analysis on the forming conditions and the Forming types of debris flow caused by landslides and falls[J].Journal of Hebei College of Geology,1996, 19(3/4)∶ 344-351.

[3]陳自生, 張曉剛.1994-04-30 四川省武隆雞冠嶺滑坡→崩塌→碎屑流→堵江災(zāi)害鏈[J].山地研究, 1994, 12(4)∶ 225-229.CHEN Zisheng, ZHANG Xiaogang.A hazard-chain of landslide→collapse→debris flow→river stoppage in Wulong county, Sichuan Province on April 30, 1994[J].Mountain Research, 1994, 12(4)∶ 225-229.

[4]朱平一, 王成華, 唐邦興.西藏特大規(guī)模碎屑流堆積特征[J].山地學(xué)報(bào), 2000, 18(5)∶ 453-456.ZHU Pingyi, WANG Chenghua, TANG Bangxing.The deposition characteristic of supper debris flow in Tibet[J].Journal of Mountain Science, 2000, 18(5)∶ 453-456.

[5]黃河清, 趙其華.汶川地震誘發(fā)文家溝巨型滑坡-碎屑流基本特征及成因機(jī)制初步分析[J].工程地質(zhì)學(xué)報(bào)，2010, 18(2)∶168-177.HUANG Heqing, ZHAO Qihua.Basic characteristics and preliminary mechanlism analysis of large scale rockslide-sturzstrom at WenJiagou triggered by Wenchuan Earthquake[J].Journal of Engineering Geology, 2010, 18(2)∶168-177.

[6]殷躍平.斜傾厚層山體滑坡視向滑動(dòng)機(jī)理研究∶ 以重慶武隆雞尾山滑坡為例[J].巖石力學(xué)與工程學(xué)報(bào), 2010, 29(2)∶217-226.YIN Yueping.Mechanism of apparent dip slide of inclined bedding rockslide∶ A case study of Jiweishan rockslide in Wulong, Chongqing[J].Chinese Journal of Rock Mechanics and Engineering, 2010, 29(2)∶ 217-226.

[7]熊道錕, 徐世民.泥石流攔擋壩之虞[J].中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào), 2010, 21(4)∶ 136-138.XIONG Daokun, XU Shiming.The fear for debris dam[J].The Chinese Journal of Geological Hazard and Control, 2010, 21(4)∶136-138.

[8]成名, 劉維甫, 劉凱欣.高速?zèng)_擊問題的離散元法數(shù)值模擬[J].計(jì)算力學(xué)學(xué)報(bào), 2009, 26(4)∶ 591-594.CHENG Ming, LIU Weipu, LIU Kaixin.Numerical simulation for high-speed impact problems by discrete element method[J].Chinese Journal of Computational Mechanics, 2009, 26(4)∶591-594.

[9]楊洋, 唐壽高.顆粒流的離散元法模擬及其進(jìn)展[J].中國(guó)粉體技術(shù), 2006, 38(5)∶ 38-43.YANG Yang, TANG Gaoshou.Discrete element method simulation of granula flow and Its advances[J].China Powder Science and Technology, 2006, 38(5)∶ 38-43.

[10]Owens J D, Luebke D, Govindaraju N, et al.A survey of general-purpose computation on graphics hardware[J].Computer Graphics Forum, 2007, 26(1)∶ 802113.

[11]吳恩華, 柳有權(quán).基于圖形處理器(GPU)的通用計(jì)算[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2004, 16(5)∶ 6012612.WU Enhua, LIU Youquan.General purpose computation on GPU[J].Journal of Computer-Aided Design & Computer Graphics, 2004, 16(5)∶ 6012612.

[12]魯曉兵, 張旭輝, 崔鵬.碎屑流沿坡面運(yùn)動(dòng)的數(shù)值模擬[J].巖土力學(xué), 2009, 30(增刊 2)∶ 524-527.LU Xiaobing, ZHANG Xuhui, CUI Peng.Numerical simulation of clastic grain flow along a slope[J].Rock and Soil Mechanics,2009, 30(Suppl 2)∶ 524-527.

[13]孫萍, 張永雙, 殷躍平, 等.東河口滑坡-碎屑流高速遠(yuǎn)程運(yùn)移機(jī)制探討[J].工程地質(zhì)學(xué)報(bào), 2009, 17(6)∶ 737-744.SUN Ping, ZHANG Yongshuang, YIN Yueping, et al.Discussion on long run-out sliding mechanism of Donghekou landslide-debris flow[J].Journal of Engineering Geology, 2009,17(6)∶ 737-744.

[14]Hsu K J.Catastrophic debris streams(sturzstroms)generated by rockfalls[J].Geological Society of America Bulletin, 1975, 86∶129-140.

[15]Heim A, von Elm der Bergsturz.Deutsch[J].Geol Zeitschr, 1882,34∶ 74-115.

[16]Buss E, Heim A.von Elm der Bergsturz[J].Zurich, Wurster &Cie, 1881∶ 163.

[17]Bagnold R A.Deposition in the process of hydraulic transport[J].Sedimentolgy, 1968, 10∶ 45-56.

[18]Bagnold R A.The nature of saltation and of bed-load transport in water[C]//Proceedings, Royal Society of London, Series A , 1973,332∶ 473-504.

[19]CU NDALL P A.A computer model for simulating progressive large scale movement in block rock system[J].Symposium ISRM, 1971, 2∶ 129-136.

[20]劉凱欣, 高凌天.離散元法研究的評(píng)述[J].力學(xué)進(jìn)展, 2003,33(4)∶ 483-490.LIU Kaixin, GAO Lingtian.A review on the discrete element method[J].Advances in Mechanics, 2003, 33(4)∶ 483-490.

[21]Liu K, Zheng W, Gao L, et al.A numerical analysis for stress wave propagation of anisotropic solids by discrete element method[C]//Chiba A, Tanimura s, Hokamoto K,eds.Proceedings of the 4th International Symposium on Impact Engineering.Kumamoto, Japan.UK∶ Elsevier Science Ltd, 2001∶ 589-594.

[22]Sawamoto Y, Tsubota H, Kasai Y, et al.Analytical studies on local damage to reinforced concrete structures under impact loading by discrete element method[J].Nuclear Engineering and Design, 1998, 179∶ 157-177.

[23]Zhu H P, Yu A B.A theoretical analysis of the force models in discrete element method[J].Powder Technol, 2006, 161∶122-129.

[24]Nezami E G, Hashash Y M.A fast contact detection algorithm for 3-D discrete element method[J].Compute Geotechnol, 2004, 31∶575-587.

[25]Mishra B K.A review of computer simulation of tumbling mills by the discrete element method[J].Part 1—Contact Mechanics.International Journal of Mineral Processing, 2003, 71(1)∶ 73-93.

[26]Hubert Nguyen.GPU Gems[M].3ed, Boston∶ Addison Wesley Professional, 2007∶ 531-549.

[27]Satish N, Harris M, Garland M.Designing efficient sorting algorithms for manycore gpus[R].Santa Clara∶ NVIDIA Technical Report NVR-2008-001, NVIDIA Corporation, 2008.

[28]李發(fā)斌, 涂建軍, 胡凱衡.基于GIS的泥石流攔砂壩優(yōu)化設(shè)計(jì)模型[J].中國(guó)科學(xué) E 輯∶ 技術(shù)科學(xué), 2003, 33(增刊)∶ 146-156.LI Fabing, TU Jianjun, HU Kaiheng.Design model of debris dam Based on GIS[J].Science in China Series E∶ Technological Sciences, 2003, 33(suppl)∶ 146-156.

[29]李德基.泥石流減災(zāi)理論與實(shí)踐[M].北京∶ 科學(xué)出版社,1997∶ 178-186.LI Deji.Theory and practice on debris flow hazard reduction[M].Beijing∶ Science Press, 1997∶ 178-186.