基于GPU的巖石碎屑流與攔砂壩交互場(chǎng)景的三維建模與可視化

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

(中國(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)格和描述破壞問(wèn)題便利的特點(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)槽的建模及兩者與碎屑流之間邊界問(wèn)題的正確處理在本研究中較重要，直接影響到巖石碎屑流與攔砂壩和排導(dǎo)槽交互場(chǎng)景的可視化表達(dá)。為此，本文作者首先分析碎屑流在場(chǎng)景中運(yùn)動(dòng)的整個(gè)過(guò)程；然后，構(gòu)建攔砂壩和排導(dǎo)槽模型，并確立它們的邊界；最后，根據(jù)碎屑流顆粒與攔砂壩和排導(dǎo)槽邊界之間的受力分析，建立有效的碎屑流顆粒與攔砂壩和排導(dǎo)槽邊界問(wèn)題的處理方法，并充分利用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ò)程通?？蓺w納為啟程階段、加速階段和堆積階段。

第一階段(啟程階段)：當(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)越小，并開(kā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]。離散元法自問(wèn)世以來(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í)攔砂壩蓄滿。若大量巖石碎屑流從壩的上方溢出，則溢出的巖石碎屑流開(kāi)始在第二級(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í)攔砂壩蓄滿為止；然后，巖石碎屑流開(kāi)始從攔砂壩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)能，在沒(méi)有阻擋物阻擋的情況下，巖石碎屑流的速度越來(lái)越快。從實(shí)驗(yàn)結(jié)果可以看出：巖石碎屑流的運(yùn)動(dòng)將對(duì)排導(dǎo)槽的水平段造成巨大的沖擊。

圖6 巖石碎屑流在沒(méi)有攔砂壩攔截下的全過(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í)配的問(wèn)題，所繪制的巖石碎屑流顆粒粒徑是均勻的，這與自然現(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]張明, 殷躍平, 吳樹(shù)仁, 等.高速遠(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)_擊問(wèn)題的離散元法數(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.