面向環境力學的離散元分析軟件研發和工程應用

季順迎　趙金鳳　狄少丞　孫珊珊

摘要：離散元計算分析軟件對解決環境力學中的離散介質問題有至關重要的作用.針對環境災害中的非規則顆粒單元，以圓球顆粒為基本單元，分別構造鑲嵌組合單元、黏結組合單元、擴展圓盤單元和擴展多面體單元等，并在此基礎上開發基于球形顆粒離散元方法的分析軟件（Software of Spherical Particlebased Discrete Element Method，SDEM）.該軟件可模擬碎冰、巖石和道砟等顆粒材料的力學行為，能直觀展現這些力學過程的發生、發展和演化；基于GPU的并行計算實現離散元大規模計算的高效性.對SDEM軟件在地質災害、工程海冰和鐵路道床等領域的應用進行介紹.

關鍵詞：離散元方法； 環境力學； 寒區工程； 地質災害； 有砟鐵路道床

中圖分類號： O347.7； TP319

文獻標志碼： B

0前言

美國學者CUNDALL等[1]在20世紀70年代提出的離散元法（Discrete Element Method，DEM） 發展迅速.特別是近年來，隨著計算機技術的發展，該方法已廣泛地應用于巖土、化工、生物和環境等諸多領域.[23]在DEM的工程應用中，一直伴隨著相應軟件的不斷發展.CUNDALL和STRACK最早研發二維圓盤程序BALL和三維圓球程序TRUBAL.[4]1989年英國Aston大學THORNTON將TRUBAL改進為TRUBALAston版，后定名GRANULE.目前，基于DEM的商業軟件主要有美國Itasca公司的PFC2D和PFC3D，英國DEMSolution公司的EDEM.此外，澳大利亞的相關學者也開展出色的工作.

在具有獨立知識產權的國產軟件研發中，吉林大學基于CAD成功解決機械部件與球體顆粒接觸作用的仿真過程.[5]北京大學基于DEM的工程計算軟件的前后處理系統，以人機交互方式實現建模、數據傳輸和動畫演示等全部過程，并主要應用于沖擊、碰撞和斷裂等動力學問題中.[6]最近，大連理工大學面向環境力學問題，自行研發出基于球形顆粒離散元方法的分析軟件（Software of Spherical Particlebased Discrete Element Method，SDEM）.該軟件通過對環境力學領域中相關工程問題的發生機理、演化過程進行數值模擬和動態顯示，為理論分析和工程應用提供數值方法和后處理分析技術.此外，基于GPU并行計算的高性能數值算法，實現離散元方法的大規模計算.目前，SDEM已成功應用到工程海冰、地質災害和鐵路道床等環境力學領域.

1SDEM中離散單元的構造

在自然條件下，鐵路道砟、浮冰和塊石等散體材料都以非規則形態存在.為此，在SDEM中通過黏結、鑲嵌和擴展等不同方式進行非規則顆粒構造，以精確描述顆粒幾何形態.

1.1鑲嵌式組合顆粒單元

針對巖土材料中塊石顆粒的非規則性，可采用具有鑲嵌模式的組合顆粒單元進行構造.該組合顆粒由不同粒徑的球形顆粒按一定組合方式構成，見圖1.對于該組合顆粒的質量、質心和轉動慣量采用有限分割法確定，并通過四元數方法進行組合顆粒的動力學分量在整體坐標和局部坐標下的轉換.[7]采用鑲嵌式組合單元可以根據設定的顆粒形態進行任意組合，應用簡便.

1.2黏結式組合顆粒單元

對于可破碎的非規則顆粒單元，可以通過球形顆粒黏結的模式進行構造，見圖2.該非規則單元可在一定作用力下按一定的破壞準則發生破碎.單元間的黏結方式可采用接觸黏結或平行黏結.[4]在模擬巖石材料的破碎特性時，一般采用平行黏結方式，通過黏結球體單元間的作用力和力矩，確定其法向應力和切向應力.采用該平行黏結模式不僅可以構造非規則顆粒單元，而且可以對連續體材料的破壞過程進行模擬.[89]

1.3擴展多面體單元

擴展多面體單元具有計算簡便、精度高的特點，可精確描述自然環境下顆粒材料的幾何形態.[1011]在一個基礎多面體模型的表面設置若干個等粒徑的球體顆粒，所有球體顆粒的球心位于基礎多面體的表面.采用該擴展模型可將空間上的一個點擴展為一個球體，將一條線段擴展為一個圓柱及柱端的兩個半球體，將任一平面擴展為由圓柱體包絡的具有一定厚度的板，將任一多面體擴展為具有光滑棱邊及棱角的多面體.采用擴展多面體模型構造的塊石顆粒見圖3.

類似于擴展多面體單元，還可構造擴展圓盤單元：首先對所需要構造的圓形單元進行函數定義，然后采用無限多個球型顆粒對其填充，以構造出所設定單元形態，見圖4.由此，圓盤單元間的相互作用轉變為球體單元的接觸判斷和作用力計算.

為便于離散元方法在環境力學領域的應用，需要研發操作簡易、應用靈活的前后處理系統.SDEM可顯示散體材料內部力鏈的分布，直觀地展現結構體與顆粒以及顆粒之間的接觸狀態.

2SDEM在環境力學中的應用

目前，SDEM已應用于工程海冰、鐵路道床和地質災害等不同研究領域.

2.1SDEM在工程海冰中的應用

2.1.1船舶在碎冰區航行的離散元模擬

近年來，隨著極地科學考察和油氣資源開發以及北極夏季通航可行性研究的開展，冰區船舶結構設計和航運性能分析引起人們更多的重視.針對碎冰區海冰的結構特點，采用三維圓盤單元描述碎冰塊，對其與船體結構的相互作用過程進行離散元分析.當海冰單元直徑在4.0～6.0 m隨機分布時，初始密集度為50%，船舶以4.0 m/s的速度行進，此時海冰與船體的相互作用的整體狀況見圖5，海冰單元與船體的接觸碰撞見圖6.[12]在不同冰速、冰厚、冰塊尺寸和冰塊密集度等海冰條件下，對海冰與船體動力作用過程進行數值分析，可為冰區船舶的安全運行和結構設計提供一定借鑒.

2.1.2海冰與海洋結構相互作用的離散元模擬

在冰區油氣開發中，海冰會導致海洋結構的強烈振動，并對結構累積疲勞損傷、上部油氣管線和設備以及操作人員健康構成很大威脅.目前，我國渤海冰區的導管架海洋平臺主要包括直立、錐體和自升式結構等.采用離散元模型可對海冰與不同類型海洋結構相互作用中的隨機破碎過程、冰載荷時程和冰激結構振動響應進行數值計算.

在海冰與直立腿導管架海洋平臺結構相互作用的離散元模擬中，將海冰離散為若干個規則排列且具有黏結破碎功能的顆粒單元，并通過海冰單軸壓縮試驗對單元間的黏結強度進行確定.在此基礎上，對直立海洋平臺結構的隨機作用過程進行數值計算，獲得不同樁徑下的冰載荷和結構冰振響應.海冰與直立海洋平臺結構相互作用過程見圖7.由此可以看出顆粒的平均作用力，也反映冰內應力場分布的大致規律[8].

在冰區油氣開發中，海洋平臺上的錐體結構可以有效降低冰力，避免強烈的冰激振動，是目前渤海油氣平臺的主要結構形式.采用離散元模型可對海冰與錐體結構的相互作用過程進行數值模擬，對比分析上錐體和下錐體冰力的大小和特征.[9]在海冰與錐體結構的相互作用過程中，海冰主要以彎曲破碎為主，顆粒間發生拉伸和剪切破壞.將錐體平臺簡化為具有等效的質量剛度阻尼的計算模型，并且不考慮錐體的變形.海冰與錐體結構相互作用過程見圖8.計算結果表明，DEM可以有效模擬海冰與錐體結構的相互作用過程，有效確定不同冰厚、冰速和冰強度條件下的冰載荷，并可分析不同錐角和錐徑下的冰載荷.

在海冰與海洋結構相互作用的離散元模擬中，結構尺度的增大對計算規模和計算效率提出迫切的需求.為此，采用CUDA編程平臺實現基于GPU的大規模離散元并行計算.為模擬大面積海冰分布，

取顆粒單元數目為16萬個，其與自升式海洋平臺的作用過程見圖9.由此可確定海冰與多樁腿結構的冰載荷，并分析海冰與海洋結構相互作用過程中的破碎、重疊和堆積特性.

2.2SDEM在地質災害中應用

2.2.1山體滑坡過程的離散元模擬

在地質災害中，邊坡的破壞形式有很多種.這里采用SDEM對重慶武隆雞尾山的垮塌過程進行離散元分析.該垮塌在空間分布上具有明顯的崩滑區、鏟刮區、滑坡主堆積區和碎屑流堆積區等四部分.[13]針對以上滑坡特性，采用具有破碎特性的黏結顆粒模型模擬巖體材料.

首先對地質體進行有限元網格劃分，并提取地質體表面的三角形節點坐標構造地質體模型，然后對關鍵塊體進行顆粒填充，見圖10.通過設置顆粒間的不同黏結強度以模擬地質體的力學性質.該山體垮塌過程的離散元模擬見圖11.由此可對該山體垮塌的形成機理和演化過程進行數值分析.

2.2.2巖石破壞中聲發射現象的離散元模擬

對于巖質邊坡，聲發射特征是典型的內部破壞狀態與宏觀可測物理量間的指標.這里采用SDEM軟件對巖石試樣在單軸壓縮破壞過程及其伴隨的聲發射現象進行數值分析.[14]巖石試樣在單軸壓縮下的離散元模擬見圖12.巖石試樣在不同壓縮應變下的力鏈結構和破壞狀態見圖13.通過離散元模擬得到的應力應變曲線及顆粒間的破碎規律可確定聲發射特征.

2.3SDEM在風沙區鐵路有砟道床中的應用

為模擬鐵路有砟道床的動力特性，分別采用組合顆粒單元和擴展多面體單元模擬道砟碎石顆粒.在此基礎上，分別對風沙貫入下的道砟力學特性和往復載荷作用下的動力特性進行離散元分析.

2.3.1含沙鐵路道砟變形模量的離散元模擬

針對風沙區有砟鐵路道床的結構特性，研究其在細沙貫入下的有效變形模量，有助于理解風沙影

響下的道床動力特性.采用離散單元模型對不同含沙率道砟材料的有效變形模量進行數值分析.

在離散元數值模擬中，細沙和道砟顆粒隨機放置于剛性圓筒內，見圖14.計算得到的不同含沙率道砟材料變形模量見圖15.可知，在低含沙率時，沙石混合體的變形模量基本保持不變；在高含沙率時，變形模量隨含沙率的增加呈線性降低.基于離散單元模型的數值模擬，還可在細觀尺度上對沙石混合體的力鏈強度、空間分布及配位數進行分析，揭示沙石混合體有效變形模量隨含沙量變化的內在機理.[15]

2.3.2鐵路道砟在往復載荷下的離散元模擬

對有砟鐵路道床在往復載荷作用下的動力特性研究有助于道床結構的合理設計.采用擴展多面體單元可構造非規則的道砟顆粒，并對道砟材料在往復載荷作用下的動力特性進行離散元分析，研究其沉降量和形變剛度的演變過程.

選用12種不同幾何形態的擴展多面體道砟顆粒，對道砟材料在往復載荷下的動力特性進行離散元分析，見圖16.道砟顆粒粒徑在36～63 mm間隨機分布.

在90個加載周期內，軌枕沉降量的變化情況見圖17.結果表明，軌枕位移在隨往復載荷上下振蕩過程中，其有效剛度和累積沉降量在初始階段波動較大，然后隨道砟材料的密實過程不斷增加并趨于穩定.該擴展多面體單元能夠合理地模擬道砟材料的動力行為，由其細觀作用過程揭示相應的宏觀演化規律.

3結束語

針對環境力學中的顆粒材料計算力學問題，研制出SDEM軟件.該軟件可通過構造鑲嵌組合單元、黏結組合單元以及擴展的圓盤和多面體單元對不同形態的顆粒單元進行相對精確的描述.此外，還開展基于GPU的大規模離散元數值計算，以更好地解決環境力學問題的計算規模和計算效率.采用SDEM對工程海冰、地質災害和風沙區有砟道床的動力特性進行離散元數值分析，有效解決相應的工程問題.

在目前SDEM基礎上，將針對環境力學中的顆粒材料計算力學問題，進一步開展顆粒滲流、復雜邊界條件和顆粒形態等方面的研究.

致謝：SDEM的研發得到美國Clarkson大學Hayley H SHEN教授的指導；美國陸軍寒區研究與工程實驗室（CREEL） Mark A HOPKINS博士對擴展多面體離散元模型的程序研發給予大力支持，在此深表謝意.

參考文獻：

[1]CUNDALL P A， STRACK O D L. A discrete numerical model for granular assembles[J]. Géotechnique， 1979， 29（1）： 4765.

[2]CLEARY P W. DEM prediction of industrial and geophysical particle flows[J]. Particuology Sci & Technol Particles， 2010， 8（2）： 106118.

[3]徐泳， 孫其誠， 張凌， 等. 顆粒離散元法研究進展[J]. 力學進展， 2003， 33（2）： 251260.

[4]POTYONDYA D O， CUNDALL P A. A bondedparticle model for rock[J]. Int J Rock Mech & Mining Sci， 2004， 41（8）： 13291364.

[5]付宏， 商慧， 于建群. 三維離散元法軟件開發研究[J]. 華中科技大學學報： 自然科學版， 2005， 33（S1）： 345348.

[6]鄭文剛， 劉凱欣. 離散元法工程計算軟件的前后處理系統[J]. 計算機工程與科學， 2000， 22（6）： 1416.

[7]YAN Ying， JI Shunying. Discrete element modeling of direct shear tests for a granular material[J]. Int J Numer & Anal Methods Geomechanics， 2010， 34（9）： 978990.

[8]季順迎， 狄少丞， 李正， 等. 海冰與直立結構相互作用的離散單元數值模擬[J]. 工程力學， 2013， 30（1）： 463469.

[9]DI Shaocheng， JI Shunying， YUE Qianjin， et al. Ice loads on conical offshore structures based on discrete element simulation[C]// Proc 21st IAHR Int Symp Ice， Dalian， 2012.

[10]HOPKINS M A. Discrete element modeling with dilated particles[J]. Eng Comput， 2004， 21（24）： 422430.

[11]HOPKINS M A， SHEN H H. Simulation of pancakeice dynamics in a wave field[J]. Annals of Glaciology， 2001， 33（1）： 355360.

[12]季順迎， 李紫麟， 李春花， 等. 碎冰區海冰與船舶結構相互作用的離散元分析[J]. 應用力學學報， 2013， 30（4）： 520526.

[13]馮振， 殷躍平， 李濱， 等. 重慶武隆雞尾山滑坡視向滑動機制分析[J]. 巖土力學， 2012， 33（9）： 27042712.

[14]JI Shunying， DI Shaocheng. Discrete element modeling of acoustic emission in rock fracture[J]. Theor & Appl Mech Lett， 2013， 3（2）： 4852.

[15]嚴穎， 狄少丞， 蘇勇， 等. 風沙影響下鐵路道砟變形模量的離散元數值分析[J]. 計算力學學報， 2012， 29（3）： 439445.

（編輯于杰）

參考文獻：

[1]CUNDALL P A， STRACK O D L. A discrete numerical model for granular assembles[J]. Géotechnique， 1979， 29（1）： 4765.

[2]CLEARY P W. DEM prediction of industrial and geophysical particle flows[J]. Particuology Sci & Technol Particles， 2010， 8（2）： 106118.

[3]徐泳， 孫其誠， 張凌， 等. 顆粒離散元法研究進展[J]. 力學進展， 2003， 33（2）： 251260.

[4]POTYONDYA D O， CUNDALL P A. A bondedparticle model for rock[J]. Int J Rock Mech & Mining Sci， 2004， 41（8）： 13291364.

[5]付宏， 商慧， 于建群. 三維離散元法軟件開發研究[J]. 華中科技大學學報： 自然科學版， 2005， 33（S1）： 345348.

[6]鄭文剛， 劉凱欣. 離散元法工程計算軟件的前后處理系統[J]. 計算機工程與科學， 2000， 22（6）： 1416.

[7]YAN Ying， JI Shunying. Discrete element modeling of direct shear tests for a granular material[J]. Int J Numer & Anal Methods Geomechanics， 2010， 34（9）： 978990.

[8]季順迎， 狄少丞， 李正， 等. 海冰與直立結構相互作用的離散單元數值模擬[J]. 工程力學， 2013， 30（1）： 463469.

[9]DI Shaocheng， JI Shunying， YUE Qianjin， et al. Ice loads on conical offshore structures based on discrete element simulation[C]// Proc 21st IAHR Int Symp Ice， Dalian， 2012.

[10]HOPKINS M A. Discrete element modeling with dilated particles[J]. Eng Comput， 2004， 21（24）： 422430.

[11]HOPKINS M A， SHEN H H. Simulation of pancakeice dynamics in a wave field[J]. Annals of Glaciology， 2001， 33（1）： 355360.

[12]季順迎， 李紫麟， 李春花， 等. 碎冰區海冰與船舶結構相互作用的離散元分析[J]. 應用力學學報， 2013， 30（4）： 520526.

[13]馮振， 殷躍平， 李濱， 等. 重慶武隆雞尾山滑坡視向滑動機制分析[J]. 巖土力學， 2012， 33（9）： 27042712.

[14]JI Shunying， DI Shaocheng. Discrete element modeling of acoustic emission in rock fracture[J]. Theor & Appl Mech Lett， 2013， 3（2）： 4852.

[15]嚴穎， 狄少丞， 蘇勇， 等. 風沙影響下鐵路道砟變形模量的離散元數值分析[J]. 計算力學學報， 2012， 29（3）： 439445.

（編輯于杰）

參考文獻：

[1]CUNDALL P A， STRACK O D L. A discrete numerical model for granular assembles[J]. Géotechnique， 1979， 29（1）： 4765.

[2]CLEARY P W. DEM prediction of industrial and geophysical particle flows[J]. Particuology Sci & Technol Particles， 2010， 8（2）： 106118.

[3]徐泳， 孫其誠， 張凌， 等. 顆粒離散元法研究進展[J]. 力學進展， 2003， 33（2）： 251260.

[4]POTYONDYA D O， CUNDALL P A. A bondedparticle model for rock[J]. Int J Rock Mech & Mining Sci， 2004， 41（8）： 13291364.

[5]付宏， 商慧， 于建群. 三維離散元法軟件開發研究[J]. 華中科技大學學報： 自然科學版， 2005， 33（S1）： 345348.

[6]鄭文剛， 劉凱欣. 離散元法工程計算軟件的前后處理系統[J]. 計算機工程與科學， 2000， 22（6）： 1416.

[7]YAN Ying， JI Shunying. Discrete element modeling of direct shear tests for a granular material[J]. Int J Numer & Anal Methods Geomechanics， 2010， 34（9）： 978990.

[8]季順迎， 狄少丞， 李正， 等. 海冰與直立結構相互作用的離散單元數值模擬[J]. 工程力學， 2013， 30（1）： 463469.

[9]DI Shaocheng， JI Shunying， YUE Qianjin， et al. Ice loads on conical offshore structures based on discrete element simulation[C]// Proc 21st IAHR Int Symp Ice， Dalian， 2012.

[10]HOPKINS M A. Discrete element modeling with dilated particles[J]. Eng Comput， 2004， 21（24）： 422430.

[11]HOPKINS M A， SHEN H H. Simulation of pancakeice dynamics in a wave field[J]. Annals of Glaciology， 2001， 33（1）： 355360.

[12]季順迎， 李紫麟， 李春花， 等. 碎冰區海冰與船舶結構相互作用的離散元分析[J]. 應用力學學報， 2013， 30（4）： 520526.

[13]馮振， 殷躍平， 李濱， 等. 重慶武隆雞尾山滑坡視向滑動機制分析[J]. 巖土力學， 2012， 33（9）： 27042712.

[14]JI Shunying， DI Shaocheng. Discrete element modeling of acoustic emission in rock fracture[J]. Theor & Appl Mech Lett， 2013， 3（2）： 4852.

[15]嚴穎， 狄少丞， 蘇勇， 等. 風沙影響下鐵路道砟變形模量的離散元數值分析[J]. 計算力學學報， 2012， 29（3）： 439445.

（編輯于杰）