離散元素法(Discrete Element Method)是在大地工程領域應用很多的一種數值
分析方法，其將連續的時間分割為多個時間步驟，並計算個別時間步驟內元素的
物理行為。由於精準的離散元素法需要計算大量元素之間的運動與碰撞，亦需要
微小的時間分割，因此在使用離散元素法時往往會有運算時間過長的問題。本研
究撰寫一套應用繪圖處理器 (GPU) 執行離散元素法之 GPGPU 程式，以降低離
散元素法之分析時間。
本研究 1)鑽研商業軟體 PFC3D實作離散元素法之演算法、流程與應用；2)利
用物件導向程式設計的概念撰寫一執行離散元素法之程式，並使用 PFC3D驗證運
算結果；3)比較碰撞偵測之演算法，並選擇出適合之演算法實作於程式之中；4)
開發新的運算引擎，透過 OpenCL 將離散元素分析流程實作於 GPU。
本研究並深入探討 GPU 與 CPU 執行離散元素法之結果，分別比較不同精度
之浮點數，以及分析不同離散元素個數時之效能與運算結果。結果顯示，GPU 的
效能在大量元素時可超過 CPU 效能 100 倍。而以目前所開發的程式而言，其效
能主要取決於記憶體頻寬 (memory bound) 而非計算效能 (compute bound)。

Discrete Element Method (DEM) is a numerical method that has many applications
in geotechnical engineering. It divides the continuous time into many time steps and
calculate the physical behavior of discrete elements in each time step. Accurate DEM
analyses needs to calculate many elements and collisions between them with tiny time
discretization and thus many time steps. Therefore, simulations using DEM takes too
long to be practical. This study develops a program that utilizes graphical processing
unit (GPU) to conduct DEM analyses with reduced computing time.
This study consists of the following four main steps. The commercial software
PFC3D
.is firstly studied to understand its inner work. A GPGPU code is then
developed using object-oriented programming to allow flexibility and extensibility.
The calculation results from the developed code are compared and validated with
PFC3D. Afterwards, several collision detection algorithm are studied and one is then
selected to be implemented in the program. Finally, the code is extended with a new
computing engine that uses OpenCL to drive the GPU device to execute the full discrete
elements analysis.
This study further compares the performance of the DEM calculated by GPU and
CPU under different problem sizes and different floating-point precisions. It is
concluded that the DEM performance using GPU can be 100 times better than CPU
when the number of discrete elements is large. It is also concluded that the
performance of the developed code is memory bound (by the memory bandwidth)
rather than bound by the compute capability (CPU bound).

[1] P. A.Cundall, “A computer model for simulating progressive, large-scale movements in blocky rock systems,” Proc. Symp. Int. Soc. Rock Mech. Soc. Rock Mech. (ISRM), Fr. II-8., pp. 129–136, 1971.
[2] Itasca Consulting Group, Particle Flow Code in3 Dimensions,Version 5.0. Minneapolis, USA: Itasca Consulting Group Inc., 2015.
[3] A. J.Walsh, “Discrete element modeling of the influences of density and confining pressure on the compression and shear behavior of granular materials,” 1998.
[4] 魏佳韻, “粒狀土壤抗剪強度異向性之數值模擬 TT - Numerical simulation of the anisotropic shear strength for granular soils.,” 國立中興大學, 台中市, 2001.
[5] 徐薇, “以三維個別元素法模擬卵礫石層的物理特性 TT - Modeling of Physical Properties in Gravel Formation Using Three-Dimensional Distinct Element Method,” 朝陽科技大學, 台中市, 2009.
[6] 林宏哲, “PFC3D探討互層岩體之力學行為 TT - none,” 國立中央大學, 桃園縣, 2013.
[7] 林郁修, “分離元素法於岩石貫切破壞試驗之模擬分析 TT - Distinct Element Approach on Rock Indentation test,” 國立臺北科技大學, 台北市, 2007.
[8] 王紹宇, “分離元素法於接觸破壞之刀刃磨耗與雙刀效應之模擬暨耦合有限差分法之數值初探 TT - Distinct Element Approach to the Effect of Wear Flat and Doubled- indenters on Contact Fracture and the Study of Coupled Numerical Method Using Distinct/Finite Different Element,” 國立臺北科技大學, 台北市, 2008.
[9] 江明軒, “FLAC耦合PFC2D應用於雙楔刀貫切破壞之初探 TT - Coupled FLAC/PFC2D Numerical Simulation to Indentation Fracture in Rock by Double Wedge Indenters,” 國立臺北科技大學, 台北市, 2008.
[10] 林耕白, “大規模崩塌地形衍育之數值模擬－甲仙案例 TT - Simulation for the Evolution of a Large Scale landslide－The Jiasian Case as an Example.,” 國立交通大學, 新竹市, 2016.
[11] 韋廷樺, “大規模崩塌地形衍育概念模型-鹿場案例 TT - A Conceptual Model for the Previous Evolution of a Large Scale landslide－The Luchang Case as an Example.,” 國立交通大學, 新竹市, 2016.
[12] R.A andZ.R, “Rockfall and Mitigation Evaluation with 3-D Discrete Element Modeling.,” pp. 1082–9, 2011.
[13] Christer Ericson, Real-Time Collision Detection. CRC Press; HAR/CDR edition, 2004.
[14] L.Lu, Z.Gu, K.Lei, S.Wang, andK.Kase, “An efficient algorithm for detecting particle contact in non-uniform size particulate system,” Particuology, vol. 8, no. 2, pp. 127–132, 2010.
[15] L. E.Walizer andJ. F.Peters, “A bounding box search algorithm for DEM simulation,” Comput. Phys. Commun., vol. 182, no. 2, pp. 281–288, 2011.
[16] E.Perkins andJ. R.Williams, “A fast contact detection algorithm insensitive to object sizes,” Eng. Comput., vol. 18, pp. 48–62, 2001.
[17] K.He, S.Dong, andZ.Zhou, “Multigrid contact detection method,” Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., vol. 75, no. 3, pp. 1–7, 2007.
[18] K.Han, Y. T.Feng, andD. R. J.Owen, “Performance comparisons of tree-based and cell-based contact detection algorithms,” Eng. Comput., vol. 24, no. 2, pp. 165–181, 2007.
[19] 林壽佑, “應用繪圖處理器於離散元素法計算之初探 TT - A Preliminary Study of Using Graphic Processors on Discrete Element Method Computation,” 國立臺灣科技大學, 台北市, 2008.
[20] Z.Langbert andM.Lewis, “Processing Hard Sphere Collisions on a GPU Using OpenCL,” World-Comp.Org, pp. 35–41.
[21] M.Durand, P.Marin, F.Faure, andB.Raffin, “DEM-based simulation of concrete structures on GPU,” Eur. J. Environ. Civ. Eng., vol. 16, no. 9, pp. 1102–1114, 2012.
[22] K.Opencl, W.Group, L.Howes, andB.Sochacki, OpenCL C Specification. Khronos OpenCL Working Group, 2016.
[23] E.Gamma, R.Helm, R.Johnson, J.Vlissides, andG.Booch, Design Patterns: Elements of Reusable Object-Oriented Software. 1994.