Cheng-Yu Ku

Pei-Jia Hsiao

Ping-Cheng Hsieh

蕭沛佳

謝平城

顧承宇

National Chung Hsing University Institutional Repository

 341土石流運動數值模式之建立 －結合水力機制與不連續變形分析法 蕭沛佳（1）  顧承宇（2）  謝平城（3）＊ 摘 要 土石流為固態顆粒與液體混合體之流動現象，對於礫石型土石流，土石流材料多以礫石等較大之土石顆粒為主體，使以連續體為主之分析方法有其限制。本研究以非連續體力學理論中之不連續變形分析法(DDA)為主體進行土石流運動機制之研究，考慮土石流係以各個大小不同之顆粒所組成，以探討土石顆粒在孔隙水壓增加後所發生之土石流化機制，同時進一步在非連續體力學機制下結合水力機制以分析土石流之問題。 （關鍵詞：不連續變形分析法、土石流、水力機制） Development of a Numerical Model for Debris Flow with the Coupling of DDA and hydro-mechanism Pei-Jia Hsiao（1） ,Ping-Cheng Hsieh（3） Graduate Student（1） and Associate Professor（3）, Department of Soil and Water Conservation National Chung Hsing University, Taichung, 402 Taiwan, R.O.C. Cheng-Yu Ku（2） Researcher（2）, Geotechnical Research Center, Sintotech Engineering Consultant. Inc. ABSTRACT A debris flow is a mixture of massive fluid and solids moving downhill under the force of gravity. Key solids in a debris flow consist of materials varying in size from clay to rocks up to 100m in length and/or width. Debris flow models based on the continuum mechanics are known to be lack of the mechanics of adequately describing the behavior of the stony debris flow and the granular flow. Hence, the Discontinuous Deformation Analysis (DDA) based on the discontinuum mechanics was employed to model a debris flow in this study.  （1）國立中興大學水土保持學系碩士班研究生 （2）中興工程顧問社大地工程研究中心研究員兼組長 （3）國立中興大學水土保持學系副教授       *通訊作者  水土保持學報 37(4)：341 - 348 (2005) Journal of Soil and Water Conservation, 37(4)：341 - 348 (2005)   342The results of this study showed that the DDA model had provided more accurate simulation of debris flows and was indeed a better simulation tool with emphasis on the particle mechanics. (Keywords：Discontinuous Deformation Analysis, Debris Flow, Hydro-mechanism) 一、前言  土石流為固態顆粒與液體混合體之流動現象，土石流體中由於含有大量之泥砂、巨石等固態物質，使其流動機制較一般流體複雜。為瞭解土石流發生之運動模式與力學機制，有從流體力學或泥砂運動力學觀點視土石流為連續之黏滯性流體，利用傳統牛頓流體模式或是土石流流變模式(如：賓漢、膨脹性及擬塑性流體)等。唯土石流流變模式研究多侷限於土石流體中之細泥成分，對於礫石型土石流，前述泥流型土石流之方法也不甚適用。目前研究與諸多災害案例顯示，土石流材料多以礫石等較大之土石顆粒為主體，使以連續體為主之分析方法有其限制。本研究將以非連續體力學理論中之不連續變形分析法(DDA) (見 Shi and Goodman, 1984)為主體進行土石流運動機制之研究，同時進一步在非連續體力學機制下結合水力機制以分析土石流之問題。(參見顧和陳, 2002; Ku and Chern, 1998; Tsesarsky and Hatzor, 2005; Wu et. al, 2004 ) 目前國際間有關非連續體之分析方法可概分為兩類，分別為分離元素法 (Distinct Element Method, DEM)與不連續變形分析法(Discontinuous Deformation Analysis, DDA) (如：Zhang and Li, 2005)。兩方法基礎理論各不相同，然而分離元素法多已被包裝成商用程式，無法對其程式碼進行修改，且其對於顆粒材料之模擬僅限於圓形或是橢圓形，對於模擬不規則形狀之土石流顆粒有限制。因此本研究之分析方法將採用 DDA 法，以顆粒力學之觀念探討土石顆粒在孔隙水壓增加時所發生之土石流化機制。(參見 Lin et. al, 1996; Jing et. al, 2001) DDA 法最先係由反算一塊體系統內之位移與變形之關係所發展(Shi and Goodman, 1984)，之後 Shi (1989)更進一步將 DDA 發展為可分析完整塊體系統之數值模式。由於DDA 法係非連續體之力學模式，其以考慮固體力學模擬為主，故本研究以 DDA 法之非連續體力學理論架構下進一步結合水力機制使其擴大以分析土石流之問題分析。 二、不連續變形分析 2.1 適用環境 Shi and Goodman(1984)提出不連續變形分析法(DDA)，具有完整的塊體運動架構，能處理離散塊體間之接觸碰撞行為及大變位，對於處理不連續體有極佳的能力，可應用於邊坡穩定、落石問題、岩盤隧道開挖與卵礫石層明挖等工程問題上。 2.2 假設與限制 對於單一塊體而言，塊體變形使用一階線性變形函數，塊體內部的應力與應變不論在任何狀態下均固定，但在許多的情況下，應力的分佈常是不規則的，在此由於塊體以一階線性變形函數計算，使得塊體的應力與應變能力受到限制。當塊體與其他塊體接觸時，塊體假設為彈性，其邊界間的剪阻抗力遵守莫爾-庫倫準則(Mohr-Coulomb criteria)。 2.3 控制方程式 2.3.1 不連續變形分析(DDA)公式     Shi(1989)以詳細的矩陣形式來呈現不連續變形分析，對於一個二維的 DDA 數值方法，根據以上的假設，採用一階近似，其不 蕭沛佳、顧承宇、謝平城：土石流運動數值模式之建立－結合水力機制與不連續變形分析法  343連續塊體結構的大位移變形函數如下： ( ) ( ) ( )( ) ( ) ( )⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−−−−−−−=⎥⎦⎤⎢⎣⎡xyyxrvuxxyyxxyyxxyyvuγεε00000000020102001  或 [ ] [ ] 166212×××=⎥⎦⎤⎢⎣⎡ DTvu  (1) 其中 [ ] [ ] [ ] txyyxtiiiiiii rvuddddddD γεε000654321 ==   (2) [ ] ( ) ( ) ( )( ) ( ) ( ) ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−−−−−−−=20102001000000xxyyxxyyxxyyT(3)       [ ]D 為位移矩陣，包含六個與位移相關的變數 ，[ ]T 為轉換矩陣，( )00 , yx 為塊體質心座標，( )yx, 為塊體內部任一點之座標，( )00 ,vu 為點( )00 , yx 在 X 與 Y 方向上的位移， 0r 為點( )00 , yx 的轉動量， ( )yx εε , 為塊體的正向應變， xyγ 為塊體的剪應變，( )vu, 為塊體內部任一點在 X 與 Y 方向上的位移。 2.3.2 最小能量法     由熱力學第二定律，在有載重(內部載重或外部載重)存在情況下，系統必須沿產生系統最小能量的方向運動或變形，對於包含 N塊塊體的系統的總勢能而言，其為所有勢能來源(即各種力和應力)的總和，包括：初始應力、外部載重、塊體變形、系統限制、塊體動力勢能等，其總勢能函數如下： ∑∫∫=+=∏Niit CdA1 21 εσ                 (4) 其中∑∫∫=NiitdA1 21 εσ 為變形勢能，C為其他勢能的總和。而對於系統中第 i 塊體而言，個別計算每個力或應力的勢能及根據最小能量原理，將其勢能對位移作微分，微分如下： N , ... , 1      ; 6 , ... , 1     ,  ==∂∏∂ irdri          (5) 以上諸項係由 6×1 子矩陣產生，所以將其加入 [ ]F 的子矩陣中。 Njisrdd sjri , ... , 1 ,      ; 6, ... , 1 ,      ,  2==∂∂∏∂     (6) 以上諸項係由 6×6 子矩陣產生，所以將其加入 [ ]K 的子矩陣中。 最後可得到第 i 塊體的 [ ]K 矩陣和 [ ]F 矩陣，系統中獨立塊體經由塊體接觸與單一塊體位移限制，利用塊體運動接觸彈簧的作用，將獨立塊體組合起來，考慮 N 個塊體的系統，其總平衡方程式如下：  FFFF DDDDN321N321321333323122322211131211⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡MMLMOMMMLLLNNNNNNNNKKKKKKKKKKKKKKKK 或 [ ][ ] [ ]FDK =                       (7) 其中 ijK 為一個 6×6 矩陣是系統中第 i 塊與第j 塊互相作用的子矩陣， iD 為 6×1 的位移子矩陣， iF 為 6×1 荷重子矩陣，而在總平衡式中 [ ][ ] [ ]FDK = 中， [ ]K 為一個 6N×6N 勁度矩陣， [ ]D 和 [ ]F 分別為 6N×1 的位移和載重矩陣。 塊體運動(block kinematics)(即塊體間無穿透和無張力作用)與沿塊體介面滑動的庫倫摩擦力(Coulomb friction)相關的多變性系 水土保持學報 37(4)：341 - 348 (2005) Journal of Soil and Water Conservation, 37(4)：341 - 348 (2005)   344統限制了式(7)的解，而其解是為了求得位移的變數，系統的最後解遵守以下的過程：首先，檢查解是否滿足限制，假如發現張力或穿透沿著任何接觸面，重新選擇新的地點和限制位置調整限制條件，並且從新的解中得到修正後的 [ ]K 和 [ ]F ，重複此過程直到沿著所有的塊體接觸面沒有張力產生和貫入發生為止。因此，藉由疊代過程(iterative process)確實地得到所給定的時間間距(time step)的位移變數。 2.4 塊體運動(block kinematics)     在塊體系統中，塊體間的接觸大致上可以分成三種形式(見圖 1)：邊對邊、角對角、角對邊，其中角對邊是基本的形式，其他兩種組合可以分解成此基本形式。當塊體一角貫入參考線(reference line)(Shi, 1989)時，貫入作用(penetration)可能會發生，不過，塊體接觸時，塊體系統必須遵守兩個條件，即塊體間無貫入發生也沒有張力作用存在，所以當塊體系統遇到上述穿透情況時，在接觸點可以增加或移除其彈力。當塊體接觸發生時，其接觸力以接觸彈簧來描述，而接觸彈簧在法線和切線方向上彈簧勁度分別是 kn、ks，如圖 2 所示(參考 Jing, 1998)。  圖 1. 塊體間接觸的形式：(a) 邊對邊接觸 (b) 角對角接觸 (c) 角對邊接觸 Figure 1. Different types of contacts between two blocks: (a) edge-to-edge; (b) vertex-to-vertex and (c) vertex-to-edge.   圖 2. 接觸間的行為 (a)正向接觸 (b)切向接觸 Figure 2. DDA mechanical behavior during contact: (a) normal contact and (b) shear contact. 2.5 靜力與動力分析     DDA 方法可以處理靜力和動力分析，對靜力分析而言，DDA 係以擬靜力 (Pseudo static)方式進行，即塊體系統中每個塊體速度在每個時間間隔的開始，假設為零。另一方面，對動力分析而言，在目前的時間間隔內塊體系統的速度為前一時刻速度的累加。 2.6 不連續變形分析的優缺點 優點： 1. 允許個別之塊體進行大範圍的塊體運動。         2. 可處理靜力分析和動力分析。 3. 在沒有任何外力作用的計算過程中，塊體破壞依據所假設的破壞準則。 缺點：         1. 需要輸入精確的塊體幾何形狀。 2. 即使是計算簡單問題的結果，也是需要長時間的計算。         3. 當實際使用來分析和設計時，要求更多的經驗與其他的方法結合。     由於 DDA 是一個相當新的方法，所以它的應用範圍依然遭到限制，這篇文章中提出 蕭沛佳、顧承宇、謝平城：土石流運動數值模式之建立－結合水力機制與不連續變形分析法  345DDA 與水力機制之結合，可使其更廣泛的應用到土石流的問題上。 三、結合水力機制分析  3.1 塊體邊界的流體壓力 施加於塊體邊界的流體壓力只是額外作用於岩石塊體的邊界條件。就剛性塊體而言來自流體壓力的作用被視為額外的淨力和力矩增量來計算。假定代表性的塊體 i 有N個邊界(表面)和N個交叉點(節點)，即表示定義第 i 塊邊界的潮濕部分的兩個點座標( 1x , 1y )和( 2x , 2y )，並且擁有水力頭 1h 和 2h  (如圖 3)。  圖 3. 塊體邊界的流體壓力 Figure 3. Hydraulic pressure on the block boundaries.  淨力值 F∆ 及其在邊緣 i 上的作用點 e 座標如下： ( ) ( ) ( )212212212 yyxxhhgF f −+−+=∆ ρ   (8) ( )( )( )( )( )( )( )( )26262121121221211212yyhhyyhhyxxhhxxhhxee+++−−=+++−−=        (9) 其中 fρ 為流體密度，g為重力加速度。而邊界 i 的方向餘弦( il , im )可以由下式計算之： ( ) ( )( ) ( )2122121221221212yyxxyymyyxxxxlii−+−−=−+−−=           (10) 並且ΔF的方向餘弦( Fl , Fm ) 可藉由解下列方程式而得： 0=+ iFiF mmll , ( ) ( ) 122 =+ FF ml       (11) 因為ΔF垂直邊界 i ，應小心確保ΔF的方向指向邊界。由流體壓力所造成塊體的力量( )yx FF , 和扭矩(τ )的增量可藉由下式計算出 ( )∑∑==∆=∆=NiiFNixx FlFF11            (12) ( )∑∑==∆=∆=NiiFNiyy FmFF11           (13) [ ]( )( ) ( )( )[ ]∑∑==∆−−∆−=∆−∆=NiFiceFiceNixyFlyyFmxxFyFx11   τ                    (14) 式中，( icx , icy )是塊體質心的座標。 3. 2 不連續塊體分析的流體壓力  在DDA裡流體壓力項使用相似的能量最小化方法予以闡述。在此常見的慣例是壓力呈線性分佈於塊體邊界兩端點( 1x , 1y ) 和( 2x , 2y ) 所定義的長度L上(見圖 4)。 ( )( )⎩⎨⎧+−=+−=112112yyyyxxxxωω                (15a)  ( ) ( )212212 yyxxL −+−=          (15b)   水土保持學報 37(4)：341 - 348 (2005) Journal of Soil and Water Conservation, 37(4)：341 - 348 (2005)   346 圖 4. 流體壓力作用在 (a) 剛體 (b) 三角形元素 (c) 四角形元素 Figure 4. Fluid pressure on (a) a rigid block (b) triangular element and (c) a quadrilateral element. 10 ≤≤ ω 代表沿著塊體邊界的相關距離。假設向量 ( ){ } ( ) ( ){ }tyx ppp ωωω ,= 沿著塊體邊界潮濕部分呈線性變化，則可以表示如下： ( ) ( )( ) ( )⎩⎨⎧+−=+−=112112yyyyxxxxppppppppωωωω       (16) 因為流體壓力所造成的剛體或元素的勢能表示如下： ( ){ } { }( ) ( ){ }( ) ( ){ } ( )[ ] [ ]{ }dXdTppLLduuppLdupyxyxyxtw⎟⎠⎞⎜⎝⎛−=⎭⎬⎫⎩⎨⎧−=−=∏∫∫∫101010,       ,       ωωωωωωωωω{ } [ ]{ }        dXFL tp−=                (17) 具有封閉形式積分，在此向量{d}是元件節點的位移向量(或剛性塊體在其形心的位移向量)，矩陣[X]是塊體或元件的形狀矩陣，向量{ PF }隨其元件的類型而變化。在式(17)裡函數 w∏ 的最小化及對節點位移向量 { }d 進行微分： { } { } [ ]XFLd tpw =∂∏∂  →  [ ]F           (18)  然後將其增加至相關元件所屬的塊體載重向量 [ ]F (load vector)，之後，將塊體的載重向量加入總體載重向量內，如總體運動方程式(即式(7))。 式(18)可用來模擬非線性壓力的分佈，在處理自由水面與元件邊界相交的自由水面位置時更有彈性，在這樣的狀態下，只有部分元件(或塊體)邊界可能被浸沒在自由水面下並且部分是保持乾燥的。而對所有浸沒的元件來說，只需設定( 1x , 1y )和( 2x , 2y )的座標等於交界的座標，即可求得浸沒在水中的部分。結合水力機制後之土石流運動如圖5所示  圖 5. 結合水力機制後以DDA進行土石流模擬示意圖 Figure 5. Modeling debris flow using DDA with the combination of fluid pressure. 四、結語  在原本的不連續變形分析法(DDA)中，尚缺乏水理方面的機制，如何考慮結合水力機制便成為一個重要的問題。本研究即在利用塊體間的水壓力，結合流體與土體間的力學機制，修正不連續變形分析(DDA)法。不過，此法目前僅考慮水體流動後所產生之水壓作 蕭沛佳、顧承宇、謝平城：土石流運動數值模式之建立－結合水力機制與不連續變形分析法  347用，未來尚需考量水動力機制以提高與水力機制結合的真實性。 五、參考文獻  1. 顧承宇、陳錦清(2002)「以流形方法分析地下水對節理岩體之影響」 ，第五屆地下水資源及水質保護研討會，雲林。 2. Jing, L. (1998), “Formulations of discontinuous deformation analysis for block systems,” Eng. Geol., vol. 49, pp. 371–381. 3. Jing, L., Y. Ma and Z. Fang (2001), “Modeling of fluid flow and solid deformation for fractured rocks with discontinuous deformation analysis (DDA) method,” International Journal of Rock Mechanics and Mining Sciences., Vol. 38, pp.343-355. 4. Ku, C.Y. and J.C. Chern (1998), “Modification of the DDA method for rigid block problems, ” International Journal of Rock Mechanics and Mining Science, Vol. 35, No.6, pp. 683-693. 5. Lin, C.T., B. Amadei, J. Jung and J. Dwyer (1996), “Extensions of discontinuous deformation analysis for jointed rock masses,” International Journal of Rock Mechanics and Mining Science & Geomechanics., Vol. 33, pp. 671-694. 6. Shi, G.H. (1989), “Block System Modeling by Discontinuous Deformation Analysis,” Ph.D. Thesis. 7. Shi, G.H. and R.E. Goodman (1984), “Discontinuous Deformation Analysis,”  8. Proceedings of the 25th U.S. Symposium  9. on Rock Mechanics, pp.269-277. 10. Tsesarsky, M. and Y.H. Hatzor (2005), “Tunnel roof deflection in blocky rock masses as a function of joint spacing and friction – A parametric study using discontinuous deformation analysis (DDA),” Tunnelling and Underground Space Technology, In Press. 11. Wu, J.H., Y. Ohnishi and S. Nishiyama (2004), “Simulation of the mechanical behavior of inclined jointed rock masses during tunnel construction using Discontinuous Deformation Analysis (DDA),” International Journal of Rock Mechanics and Mining Sciences., Vol. 41, pp. 731-743. 12. Zhang, M., H. Yang and Z. Li (2005), “A Coupling Model of the Discontinuous Deformation Analysis Method and the Finite Element Method,” Tsinghua Science & Technology., Vol. 10, pp. 221-226.   94 年 10 月 04 日 收稿 94 年 11 月 05 日 修改 94 年 11 月 20 日 接受         水土保持學報 37(4)：341 - 348 (2005) Journal of Soil and Water Conservation, 37(4)：341 - 348 (2005)   348 

Development of a Numerical Model for Debris Flow with the Coupling of DDA and hydro-mechanism

http://ir.lib.nchu.edu.tw/bitstream/11455/78739/1/37-4-3.pdf

Development of a Numerical Model for Debris Flow with the Coupling of DDA and hydro-mechanism

Abstract

Similar works

Full text

Available Versions

National Chung Hsing University Institutional Repository