國立交通大學
土木工程學系
碩士論文
泥流型和礫石型土石流有限元素分析初步探討
The Preliminary Study of Finite Element Analysis of
Mudflow and Stony Debris Flow
研 究 生：王靖
指導教授：潘以文 博士
泥流型和礫石型土石流有限元素分析初步探討
The Preliminary Study of Finite Element Analysis ofMudflow and Stony Debris Flow
研 究 生：王靖 Student：Ching Wang 指導教授：潘以文 博士 Advisor：Dr. Yii-Wen Pan
國立交通大學 土木工程學系
碩士論文
A Thesis
Submitted to Department of Department of Civil Engineering College of Engineering
National Chiao Tung University In partial Fulfillment of the Requirements
For the Degree of Master
In
Civil Engineering
July 2013
Hsinchu, Taiwan, Republic of China
中華民國一 O 二年七月
泥流型和礫石型土石流有限元素分析初步探討 學生：王靖 指導教授：潘以文 博士 國立交通大學土木工程學系碩士班
中文摘要
本研究嘗試運用流體動力學之二相流分析以模擬土石流之運移 與 堆 積 行 為 。 數 值 模 擬 係 以 COMSOL Multiphysics （ 以 下 簡 稱 COMSOL）有限元素模擬軟體為工具，將土石流視為黏性流體，以 等位函數(Level Set) 法分辨空氣與土石流之自由介面，以便探討土石 流運移的特性。COMSOL 可同時輸入多重物理量，其中計算流體力 學(CFD)模組可以同時模擬二相流之流體行為。土石流模擬所考慮之 控制方程式包括連續方程式與動量方程式(Navior-Stokes 方程式)，土 石流的流變模式為非牛頓流，可包含降伏剪應力、黏滯力、顆粒碰撞 產生的離散剪應力及土石流流床的粗糙度所產生的紊流力等項。 回顧檢討前人對土石流流變模式和參數的研究，可知土石流流變 性質會隨顆粒粒徑分布、礦物性質、土石流的泥砂體積濃度等物理量 的不同而有所差異，應合理考慮土石流的流變模式以及其流變性質。 本研究整理文獻收集的土石流流變模式參數範圍，以供為隨後模擬的 依據。本研究繼而以 COMSOL 為分析平台，建立土石流運移的模型， 允許兼而考量適合描述泥流型與礫石型土石流的複合流變模式。針對 泥流型土石流的流變模式，採用等效賓漢流變行為的黏滯係數，令原 本用於模擬牛頓流體的 Navior-Stokes 方程式可以模擬賓漢流體的流 變行為。針對礫石型土石流的流變模式，考慮顆粒碰撞阻力/紊流力 與剪應變率平方正相關，將顆粒碰撞阻力/紊流力以體積力（Volume force）加入原本用於模擬牛頓流體之 Navior-Stokes 方程式。流變模 式中降伏剪應力為常數、黏滯力與流速流度梯度成正比，顆粒碰撞阻 力/紊流力與流速流度梯度平方成正比。當忽略最後一項，此複合流 變模式即簡化至賓漢流變模式。 本研究檢討上述之模擬方法是否可合理地模擬土石流啟動後之 運移與推積行為，並探討不同因子、流變模式成分、與模式參數值如 何影響土石流啟動後之運移與堆積行為。本研究中所使用等位函數法 (Level set Method)旨在區別土石流和空氣來做介面，但卻會發生土石 流質量不守恆之問題，網格大小和二相流介面厚度可影響質量不守恆 之程度。
關鍵字：土石流、有限元素、等位函數法(Level set Method) 、Comsol 軟體、流變模式
The Preliminary Study of Finite Element Analysis of Mudflow and Stony Debris Flow
Student：Ching Wang Advisors：Dr. Yii-Wen Pan
Abstract
This study aims to investigate the behavior of migration and deposition for debris flow using two-phase flow analysis of fluid dynamics. This work made use of “COMSOL Multiphysics” (herein referred to as "COMSOL") as the numerical simulation tool to simulate fluid behavior of two-phase flow. The rheological models of debris flow are assumed to be Non-Newtonian fluid, including yield shear stress, viscous stress, dispersive shear stress from particles collision, and turbulent stress from bed roughness, and so on.
Reviewing rheological models and parameters for debris flows from previous studies, debris flow rheological properties vary with the particle size distribution, mineral properties, debris flow’s sediment concentration, and other physical quantities. It should be reasonable to consider the debris flow rheological model and its corresponding rheological properties together. This study deliberately validates the ranges of parameters which are collected in order, and view them as basis for simulation.
This study regards COMSOL as an analysis platform to build debris flow migration model and to design debris flow rheological model as well. With the rheological model for mudflow, this study uses effective viscosity from the effective Bingham rheological behavior so that the Navior-Stokes equation originally simulating Newtonian fluid is able to
simulate rheological behavior of Bingham fluid. For rheological model of stony debris flow, considering that particles collision force and turbulent force are square proportional to shear rate, this study adds particles collision force and turbulent force as volume forces to the Navior-Stokes equation originally simulating Newtonian fluid. Assuming that the yield shear stress is constant, viscous flow is proportional to the velocity flow gradient; the shear stress due to particle collision is proportional to the square of flow rate of gradient. When the shear stress is neglected, this rheological model is reduced to the Bingham model.
This study evaluates whether the behavior of the migration and deposition of debris flow is reasonably simulated. The major role of the Level Set method is to distinguish the interface between debris flow and air. However, the problem of poor mass conservation for debris flow appears. It appears the mesh size and the interface thickness of two phase flow have an influence on the quantity of poor mass conservation.
Key Words：Debris flow; Finite Element Method; Comsol Multiphysics; Rheological Parameters; Level set Method
致謝
這段時間相當感謝恩師 潘以文老師，由於我的題目是用到實驗 室以前沒人用過的模擬工具，在模擬的建構和調整上，也多虧潘以文 老師願意每週找我討論模擬對策，使研究更加完善。在土石流文獻的 部分，感謝廖志中和潘以文老師對文獻部分做篩選，使模擬所用的文 獻能更貼近本研究所要的。即使有時研究做到後面失去方向或遇到瓶 頸而茫然導致進度延後，潘老師還是有耐心的方式拉回研究的正確方 向。感謝口試委員吳建宏老師、王承德老師、吳博凱老師及吳柏林博 士，費心審閱學生論文提出寶貴意見。課堂上承蒙方永壽老師、黃安 斌老師、單信瑜老師及林志平老師於學識的指導。做研究期間也相當 辛苦，非常感謝李國維(貓博)學長在這段時間對我的照顧。我在這很 感謝偉欽學長給我的研究立下起點，也感謝柏林學長和皮托科技的工 程師前輩們提供給我 COMSOL 的使用對策。也很感謝明萬、阿康、 慧蓉、里泉、麵包、阿蒲、小昱學長姐提供學習上的幫助和做研究的 指引。也感謝同門的婉容、阿村、包包、輝凱和芋頭同學一起為各自 的研究拼命。同時也感謝小短、咪咪、宥達、孝存、建文、俊廷在口 試場地幫忙和提供電腦供我跑模擬。在這碩士生涯中，我感謝我父親 教導我做研究的方法和該要有的態度和提供我畢業後的規劃，還有耐心聆聽我問題的母親。以及禪學社的師兄姐、聖研社的們的弟兄們、 我妹和知己好友小白學弟的陪伴。感恩！
目錄
中文摘要 ... i Abstract ... iii 致謝 ... v 目錄 ... vii 圖目錄... xi 表目錄... xix 符號表... xxi 第一章 緒論 ... 1 1.1. 研究動機與目的 ... 1 1.2. 研究目的 ... 3 1.3. 研究內容和流程 ... 4 第二章 文獻回顧 ... 6 2.1. 土石流種類 ... 6 2.2. 土石流的發展過程 ... 15 2.3. 土石流發生特性和運作 ... 17 2.4. 土石流流變模式 ... 212.4.1. 牛頓流體流變模式(Newtonian fluid model) ... 21
2.4.3. 賓漢流變模式(Bingham model) ... 25 2.4.4. 赫許布雷流變模式(Herschel-Bulkley model) ... 35 2.4.5. 碰撞/膨脹流變模式(Collision/Dilatant model) ... 39 2.4.6. 廣義黏塑性流變模式 ... 45 2.4.7. 二項式流變模式(Quadratic model) ... 46 2.4.8. Voellmy 流變模型 ... 47 2.5. 土石流流速 ... 48 2.6. 土石流規模與流量 ... 53 2.6.1. 土石流體積估計 ... 53 2.6.2. 土石流流量與洪水流量 ... 55 2.6.3. 土石流尖峰流量和規模 ... 56 2.7. FLO-2D 簡介 ... 57 2.7.1. 控制方程式 ... 58 2.7.2. FLO-2D 流變模式和組成律 ... 60 2.7.3. FLO-2D 輸入參數和案例 ... 60 2.8. 土石流數值模擬之文獻討論... 64 2.8.1. 土石流現地估計和流變計實驗 ... 64 2.8.2. 土石流模擬數值方法 ... 70 2.9. COMSOL 雪崩案例回顧 ... 83
第三章 研究方法 ... 91 3.1. COMSOL Multiphysics 簡介和處理程序 ... 93 3.1.1. COMSOL Multiphysics 簡介 ... 93 3.1.2. COMSOL Multiphysics 的處理程序 ... 93 3.1.3. 有限元素法簡介 ... 95 3.2. 數值模型... 97 3.2.1. Navior-Stoke 方程式和連續方程式 ... 97 3.2.2. 從流變模型角度修改 Navior-Stoke 方程式 ... 99
3.2.3. 等位函數法(Level set Method)方程式 ... 104
3.3. 模擬參數的設定... 106
3.4. 幾何模型和邊界條件的選擇 ... 107
3.5. 網格設定... 115
3.6. 求解器和收斂標準 ... 117
3.6.1. 求解器簡介 ... 117
3.6.2. COMSOL 的收斂準則(Convergence criteria) ... 120
3.7. 後處理 ... 121
3.8. 虛擬案例研究規劃 ... 122
第四章 方法驗證與參數探討 ... 130
4.1.1. 牛頓流情況 ... 133 4.1.2. 非牛頓流情況 ... 136 4.2. 密度的探討 ... 137 4.3. 黏滯係數的探討... 139 4.3.1. 黏滯係數為常數之條件 ... 139 4.3.2. 賓漢雙黏度流變之條件 ... 142 4.3.3. Cross 等效賓漢流之條件 ... 145 4.4. 碰撞-紊流係數探討 ... 148 4.5. 等位函數法質量守恆問題探討... 151 4.5.1. 網格測試 ... 152 4.5.2. 介面厚度大小測試 ... 160 第五章 虛擬流槽土石流之運動模擬 ... 161 5.1. 坡度大小探討 ... 162 5.2. 泥流型土石流體積濃度影響 ... 171 5.3. 礫石型土石流之碰撞-紊流係數影響 ... 182 第六章 結論與建議 ... 194 6.1. 結論 ... 194 6.2. 建議 ... 196 參考文獻 ... 198
圖目錄
圖 1-1 研究流程圖 ... 5 圖 2-1 水砂含量和運動關係圖(Coussot et al,1996)... 9 圖 2-2 天然壩破壞引發的土石流(kuang et al.，1977) ... 10 圖 2-3 土石流發展示意圖(Johnson et al.，1984) ... 15 圖 2-4 土體受侵蝕而形成土石流(Tognacca et al.，1997) ... 16 圖 2-5 山崩誘發型土石流的的流動過程(Takahashi，2001) ... 17 圖 2-6 土石流流動剖面圖(詹錢登，1998) ... 19 圖 2-7 土石流區域分布圖(詹錢登，1998) ... 19 圖 2-8 土石流運移過程圖(謝正倫，1991) ... 20 圖 2-9 常用的土石流變模型(Tatahashi,2007) ... 21 圖 2-10 將半橢圓截面的流道簡化成矩形截面的流道 ... 27 圖 2-11 含粗顆粒漿體流變特性之分析概念 ... 30 圖 2-12 以賓漢流體為基底，前人實驗賓漢黏滯係數和體積濃度的關 係 (間接摘自 Takahashi(2007)) ... 31 圖 2-13 賓漢模型(實線部分)及賓漢雙黏性模型(虛線部分)的比較 (Dentand et al.，1986) ... 32 圖 2-14 賓漢模型和 Cross 等效賓漢模型的比較 ... 34 圖 2-15 在不同剪切率的 Bingham 流體直線和 Herschel-Bulkley 流體曲線比較(摘自王志賢，2007) ... 38 圖 2-16 高速顆粒流顆粒之間的應力關係 (Tsujimoto et al.，1997) . 41 圖 2-17 顆粒體積濃度分別對應碰撞應力和顆粒動能應力的大小分 布(圖中的σ為𝜌𝑠，摘自於 Tsujimoto et al. (1997))... 43 圖 2-18 賓漢流體流速分布 ... 48 圖 2-19 山洪的流域面積和土石流流出土方量關係 (Mizuhara,1990) ... 54 圖 2-20 FLO-2D 所使用的坐標系統 ... 58 圖 2-21 土石流潛勢溪流（花蓮 069）土石流運動模擬 ... 63 圖 2-22 沖積扇的模擬和實驗的結果比較(Takahashi,1987) ... 71 圖 2-23 山崩土塊的數學模組(間接摘自於 Takahashi(2007)) ... 72 圖 2-24 山崩誘發型土石流的數值模型 (Satofuka et al. 2003) ... 72 圖 2-25 數值分析網格和六邊形角柱土塊邊界剪應力 ... 73 圖 2-26 水俣市寶川內地區某邊坡破壞前後的形狀(Satofuka,2004) . 74 圖 2-27 依邊坡破壞深度，所建模出來的土塊初始深度分布圖 (Satofuka,2004) ... 74 圖 2-28 深度和土塊高度分布之模擬結果(Satofuka,2004) ... 75
圖 2-29 Rickenmann et al. (2006a)用 HB 流變模式模擬瑞士 Varunasch 集水區土石流的堆積厚度。 ... 76
圖 2-31 以特徵線法計算土石流停止後之前緣位置與模擬結果之比較 ... 77 圖 2-32 以特徵線法計算中央斷面之土石流流深分佈並與模擬結果 作比較(黃名村等人，2002) ... 78 圖 2-33 用 DFEM-2D 模擬 Wartschenbach 發生土石流時堆積厚度。 ... 79 圖 2-34 一維斜面上的潰壩型土石流(Martinez et al.，2008) ... 80 圖 2-35 土石流運移距離解析解、實驗數據和數值解的比較 ... 80 圖 2-36 土石流對直立壩的模擬畫面 (聶志成，2007) ... 81 圖 2-37 土石流對直立壩的流場模擬畫面 (聶志成，2007) ... 81 圖 2-38 直立壩所受衝擊力的分布 (聶志成，2007) ... 82 圖 2-39 直立壩在模擬時間內的的衝擊力的變化 ... 82 圖 2-40 雪崩模組勢函數分布圖示(Bovet et al., 2007) ... 84 圖 2-41 雪崩模組模擬所使用的網格(Bovet et al.,2007) ... 87 圖 2-42 牛頓流雪崩模擬結果 (Bovet et al.,2007) ... 87 圖 2-43 牛頓流雪崩模擬在 t=1.2 秒時流速的分布 ... 88 圖 2-44 非牛頓流雪崩模擬在 t=1.2 秒時雪崩黏滯係數大小的分布 (Bovet et al.,2007) ... 88 圖 3-1 數值求解流程圖 ... 95 圖 3-2 流槽模組在 t=0 時 Level set 平滑函數φ的大小分布 ... 106
圖 3-4 流槽模組壁邊界條件(Wall Boundary Condition)的分布 ... 109
圖 3-5 壁邊界條件依滑移距離 b 來分類 ... 109
圖 3-6 底部的邊界條件設為 No Slip Wall 容易使空氣參雜在底部 110 圖 3-7 No slip wall 邊界條件的位置 ... 112
圖 3-8 Partial Slip Wall 邊界條件的位置 ... 113
圖 3-9 土石流和空氣的兩相介面位置 ... 113 圖 3-10 流槽頂部邊界用對稱邊界 ... 114 圖 3-11 右端邊界設成出口邊界 ... 115 圖 3-12 坡道為 20 度的建模以 CFD 模組規格正常網格大小所建的網 格，網格數約為 5071 個... 117 圖 3-13 虛擬案例模擬規劃的流程圖 ... 129 圖 4-1 COMSOL 3.5a 原檔所使用的幾何模型 ... 130 圖 4-2 COMSOL 4.3b 所重建的雪崩模型 ... 131 圖 4-3 COMSOL 4.3b 所重建的雪崩模型所使用的網格 ... 132 圖 4-4 以 COMSOL 4.3b 重建雪崩模型 t=0(s)的模擬結果 ... 133 圖 4-5 以 COMSOL 4.3b 重建雪崩模型 t=1(s)的模擬結果 ... 134 圖 4-6 以 COMSOL 4.3b 重建雪崩模型 t=2(s)的模擬結果 ... 134 圖 4-7 以 COMSOL 4.3b 重建雪崩模型 t=4(s)的模擬結果 ... 135 圖 4-8 以 COMSOL 4.3b 重建雪崩模型 t=12(s)的模擬結果 ... 135
圖 4-9 以 COMSOL 4.3b 重建雪崩模型 t=1.2(s)的流速模擬結果 .... 136 圖 4-10 非牛頓流雪崩黏滯係數分布結果 ... 137 圖 4-11 密度為300 𝑘𝑔 𝑚⁄ 3的流體在 t=20(s)的流動情形 ... 138 圖 4-12 密度為1907.5 𝑘𝑔 𝑚⁄ 3的流體在 t=20(s)的流動情形 ... 139 圖 4-13 黏滯係數為 20Pa*s 在 t=5.5(s)的流速分布圖 ... 140 圖 4-14 黏滯係數為 400Pa*s 在 t=5.5(s)的流速分布圖 ... 140 圖 4-15 黏滯係數為 20Pa*s 在 t=60(s)的流體流況圖 ... 141 圖 4-16 黏滯係數為 400Pa*s 在 t=60(s)的流體流況圖 ... 142 圖 4-17 降伏應力 5000Pa、黏滯係數為 50Pa*s 在 t=0(s)的黏滯係數大 小分布 ... 143 圖 4-18 降伏應力 5000Pa、黏滯係數為 50Pa*s 在 t=0.1(s)的黏滯係數 大小分布 ... 144 圖 4-19 降伏應力 5000Pa、黏滯係數為 50Pa*s 在 t=0.7(s)的黏滯係數 大小分布 ... 144 圖 4-20 降伏應力 5000Pa、黏滯係數為 50Pa*s 在 t=0.8(s)的黏滯係數 大小分布 ... 145 圖 4-21 降伏應力 5000Pa、黏滯係數為 50Pa*s 在 t=0(s)的等效黏滯係 數大小分布(白線為土石流和空氣的介面形狀) ... 146 圖 4-22 降伏應力 5000Pa、黏滯係數為 50Pa*s 在 t=0.5(s)的等效黏滯
係數大小分布(白線為土石流和空氣的介面形狀) ... 146 圖 4-23 降伏應力 5000Pa、黏滯係數為 50Pa*s 在 t=1(s)的等效黏滯係 數大小分布(白線為土石流和空氣的介面形狀) ... 147 圖 4-24 降伏應力 5000Pa、黏滯係數為 50Pa*s 在 t=2(s)的等效黏滯係 數大小分布(白線為土石流和空氣的介面形狀) ... 147 圖 4-25 P2+P1 離散選相因誤差放大而產生的高密度為例 ... 149 圖 4-26 t=0 秒時碰撞紊流的阻抗分布 ... 150 圖 4-27 t=3 秒時碰撞紊流的阻抗分布 ... 150 圖 4-28 t=12 秒時碰撞紊流的阻抗分布 ... 151 圖 4-29 模擬測試所使用正常大小網格 ... 152 圖 4-30 模擬測試所使用粗化大小網格 ... 153 圖 4-31 模擬測試所使用細化大小網格 ... 154 圖 4-32 模擬測試所使用正常大小加上底部邊界加密的網格 ... 155 圖 4-33 使用適性網格約 0~6 秒所使用的網格 ... 156 圖 4-34 模擬 t=6(s)時，使用圖 4-33 網格的情形 ... 156 圖 4-35 使用適性網格約 6~12 秒所使用的網格 ... 156 圖 4-36 模擬 t=6(s)時，改用圖 4-35 網格的情形 ... 157 圖 4-37 依網格大小所造成的面積流失測試 ... 158 圖 4-38 固定介面厚度，網格對面積損失的影響 ... 159
圖 4-39 介面厚度大小設定對體積變化率的影響 ... 160 圖 5-1 模擬面積守恆情形 ... 163 圖 5-2 模擬最大流速測量結果(模擬時間 60 秒) ... 164 圖 5-3 模擬最大流速結果(取前面 20 秒測量結果) ... 164 圖 5-4 模擬平均流速測量結果(模擬時間 60 秒) ... 165 圖 5-5 模擬堆積區流動距離測量結果(模擬時間 60 秒) ... 166 圖 5-6 [(tpf.Vf1>=0.5)*(y)]的表面最大值的測量結果 ... 167 圖 5-7 坡度為 15 度的土石流模擬至 24.7 秒的堆積情形 ... 169 圖 5-8 將圖 5-7 的等位線輸出的座標值 ... 169 圖 5-9 坡度為 18 度的土石流模擬至 21.8 秒的堆積情形 ... 170 圖 5-10 將圖 5-9 的等位線輸出的座標值 ... 170 圖 5-11 模擬之面積守恆情形 ... 172 圖 5-12 模擬最大流速測量結果(模擬時間 60 秒) ... 173 圖 5-13 模擬最大流速結果(取前面 20 秒測量結果) ... 174 圖 5-14 模擬平均流速測量結果(模擬時間 60 秒) ... 174 圖 5-15 模擬堆積區流動距離測量結果(前 7 秒的模擬結果) ... 175 圖 5-16 模擬堆積區流動距離測量結果(模擬時間 60 秒) ... 175 圖 5-17 模擬[(tpf.Vf1>=0.5)*(y)]的表面最大值的測量結果 ... 176 圖 5-18 體積密度為 45%土石流模擬至 26.6 秒的堆積情形 ... 177
圖 5-19 將圖 5-18 的等位線輸出的座標值 ... 177 圖 5-20 體積密度為 55%土石流模擬至 20.4 秒的堆積情形 ... 178 圖 5-21 體積密度為 65%土石流模擬至 18.9 秒的堆積情形 ... 179 圖 5-22 模擬至 13 秒，體積濃度為 45%的泥流型土石流 ... 183 圖 5-23 模擬至 13 秒，體積濃度為 45%、𝜇𝑐𝑡為 0.5 的礫石型土石流 ... 183 圖 5-24 模擬的面積守恆情形 ... 184 圖 5-25 模擬最大流速測量結果(模擬時間 60 秒) ... 185 圖 5-26 模擬最大流速結果(取前面 20 秒測量結果) ... 185 圖 5-27 模擬平均流速測量結果(模擬時間 60 秒) ... 185 圖 5-28 模擬堆積區流動距離測量結果(前 7 秒的模擬結果) ... 186 圖 5-29 模擬堆積區流動距離測量結果 ... 187 圖 5-30 模擬[(tpf.Vf1>=0.5)*(y)]的表面最大值的測量結果 ... 187 圖 5-31 𝜇𝑐𝑡為 0.5 的土石流在 25 秒時的堆積情形 ... 188 圖 5-32 𝜇𝑐𝑡為 1 的土石流在 25 秒時的堆積情形 ... 189 圖 5-33 𝜇𝑐𝑡為 1.5 的土石流在 25.5 秒時的堆積情形 ... 189
表目錄
表 2-1 不同土石流組成分類(熊剛，1996) ... 8 表 2-2 由無因次量範圍區分土石流(Julien et al.，1991) ... 12 表 2-3 無因次參數大小參考值（Iverson，1997） ... 13 表 2-4 土石流在運移過程中因坡度而所產生的流動特性 ... 20 表 2-5 1998 年義大利 Acquabona 溪流流域土石流監測結果... 29 表 2-6 各研究者所提出之 Herschel-Bulkley 模式𝜗值範圍 ... 37 表 2-7 土石流平均流速整理表，摘自 Rickenmann (1999)。 ... 52 表 2-8 台灣土石流土砂流出量迴歸公式表 (行政院農委會,2001) .... 54 表 2-9 FLO-2D 前人使用與流變參數決定方式 ... 62 表 2-10 花蓮土石流潛勢溪流 FLO-2D 各參數輸入值 ... 63 表 2-11 現地估計方法所得之降伏應力與黏滯係數整理表 ... 66 表 2-12 現地取樣後室內試驗所得之賓漢降伏應力 ... 67 表 2-13 NS 模組參數表格(Bovet et al.,2007) ... 89 表 2-14 CD 模組參數表格(Bovet et al.,2007) ... 90 表 3-1 各種土體飽和狀態下的內摩擦角(周，1991) ... 103 表 3-2 坡度參數參數探討之相關參數組合 ... 127 表 3-3 體積濃度參數探討之相關參數組合 ... 127 表 3-4 礫石型模型參數探討之參數組合 ... 128表 4-1 質量守恆性模擬測試所使用的網格的資料 ... 158 表 5-1 體積濃度對應其相關物理量 ... 181
符號表
2.1~2.8 的部分 (x 表示長度座標、y 表示寬度座標、z 表示高度或深度的座標) 𝑇𝑑∗ 無因次紊流力 𝐷𝑣∗ 無因次碰撞力 𝜏∗ 無因次剪力 𝑁𝑆𝑎𝑣 Savage 無因次指標數 𝑁𝐵𝑎𝑔 Bagnold 無因次指標數 𝑁𝑓𝑟𝑖𝑐 摩擦力指標數 𝜌𝑠 土石固態顆粒密度(粗粒料為主) 𝜌𝑓 流體密度 𝑑𝑝 顆粒平均直徑 𝛾̇ 剪(切)應變率(此處定義為𝑑𝑢 𝑑𝑧⁄ ) μ 流體黏滯係數(對流體狀態黏塑性物質也稱為動黏滯係數) μ0 水的黏滯係數 CV 體積濃度 𝜌𝑚 泥流土石流密度(只含粒徑為 2mm 以下的細粒料) 𝜌𝑤 水密度 𝐺𝑠 土石固體比重𝑤𝑐 含水量 𝐶𝑣𝑚 顆粒之最大極限濃度 σ𝑛 顆粒間相互接觸在剪力面產生的正向力 c 凝聚力 φ𝑐 顆粒內摩擦角 τB 賓漢降伏應力 μB 賓漢黏滯係數 Tc 臨界流動深度 Θ 堆積傾斜坡度 ρb 特殊大石塊單位重 𝑙𝑚𝑖𝑛 特殊大石塊在 x 軸、y 軸和 z 軸中最短的長度 𝑛𝑣 特殊大石塊被水淹沒的體積佔所有體積的比例 W𝑝 塞流(弱剪區)寬度 h 流動深度 𝑊 流道寬度 𝛼1 賓漢降伏應力迴歸參數 𝛽1 賓漢降伏應力指數迴歸參數 𝛼2 賓漢黏滯係數迴歸參數 𝛽2 賓漢黏滯係數指數迴歸參數
Cvg 純粗粒料部分的體積濃度 Cvf 純細粒料部分的體積濃度 𝛾𝑙𝑠𝑟̇ 較小的剪(切)應變率 𝜇𝑙𝑠𝑟 較小的剪(切)應變率狀態下的黏滯係數 𝛾ℎ𝑠𝑟̇ 較大的剪(切)應變率 𝜇ℎ𝑠𝑟 較大的剪(切)應變率狀態下的黏滯係數 𝜇𝑒𝑓𝑓 Cross 流變模型的等效黏滯係數 τHB 赫許布雷流變模式的降伏應力 μHB 赫許布雷流變模式的流體行為指標 𝜗 流體稠度指標或流體行為指標 𝜇𝑑 膨脹/碰撞流體行為參數 𝑃𝑑 顆粒離散應力 λ 線性濃度 𝐶𝑏 顆粒靜止接觸時的最大體積濃度 𝜏𝑑 離散剪應力 𝛼𝑖 顆粒離散應力的方向和垂直線(z 軸)之間的夾角 𝑎𝑖 顆粒運動係數 𝜏resist 土石流流床對土石流所產生的阻抗剪應力 𝜏𝑐𝑜𝑙𝑙𝑖 土石顆粒碰撞所產生的剪應力
h 流動深度 𝜌𝑑 土石流整體密度(皆含粗細粒料) g 重力加速度 𝑈𝑎𝑣 土石流截面的平均流速 𝑓𝑠 顆粒碰撞對流速所產生的摩擦係數 𝜏𝑐 碰撞剪應力 𝑔0 (線性濃度+1) 𝑒 碰撞恢復係數 𝜏 顆粒動能應力 𝑓𝑝 顆粒運動的速度的高斯分布函數 𝐶𝑡 平均紊流速度 T 溫度 b𝑝𝑟 顆粒運移的物理量 Q𝑝𝑟 顆粒運移的物理量的平均值 𝑙′ 顆粒運移的距離 𝑡𝑓 顆粒運移的時間 μgvn GVN 模式的流體行為指標 𝜏𝑐𝑦 凝聚性降伏應力 𝜏𝑚𝑐 莫爾-庫倫剪應力
𝜏𝑣 黏滯剪應力 𝜏𝑡 紊流剪應力 𝜏𝑐 碰撞離散剪應力 𝜏OB 二項式流變的降伏應力 𝑙𝑚 Prandtl 混摻長度係數 𝑓𝑣 岩石間摩擦力 n 曼寧係數 H 最大流動深度 S 河床坡降 𝑧′ 從表面往下算，其剪力總和等於賓漢降伏剪應力的深度。 τB ρ ⁄ 𝑢∗ 剪力速度 √𝑔 𝑠 𝑛 (=𝑧 ⁄ ) 𝐶∞ 土石流的平衡濃度 ́ ́ 常數 ∗ (=𝑢∗ 𝜇⁄ ) 𝑜 (= 2 𝑎 𝑖𝑠 𝑛𝛼𝑖⁄ 𝜌2 𝑠⁄ 𝑑ρ𝑑 𝑝⁄ 2) 0 (=𝑎𝑜𝜇𝑜⁄𝑢∗ ) A𝑑 山洪的流域面積 𝑀𝑜𝑢𝑡 土石流流出土方量 M
qs 單位寬度土石流量 qw 水流流量 𝐶∗ 通常為堆積層土體濃度(其定義為(1-土體孔隙率)*100%) 𝑆𝑏 堆積層含水飽和度 Qp 土石流尖峰流量 M* 歷史平均土石量體 u x 軸方向流速 v y 軸方向(寬度)流速 i 降雨強度 t 時間 S𝑓 摩擦坡降 𝑆𝑓𝑥 x 軸方向摩擦坡降 𝑆𝑓𝑦 y 軸方向摩擦坡降 𝑆𝑏𝑥 x 軸方向底床坡降 𝑆𝑏𝑦 y 軸方向底床坡降 𝑆𝑦 降伏坡降 𝑆𝑡𝑐 紊流-碰撞坡降 K 層流阻滯係數 Cw 重量濃度
𝜏𝑖 剪應力張量 𝜏𝑦 廣義降伏應力 𝛾𝑖 剪應變張量 𝑞𝑥 x 方向單位寬度的流量 𝑑 下游堆積區的傾斜度 𝐶𝑝𝑐 pseudo-Che´zy 常數 2.9、第三章、第四章和第五章的部分 (x 表示長度座標、y 表示高度或深度的座標) 𝜇1 兩相流流體 1 的黏滯係數(指雪崩) 𝜌1 兩相流流體 1 的密度(指雪崩) 𝜇2 兩相流流體 2 的黏滯係數(指空氣) 𝜌2 兩相流流體 2 的密度(指空氣) 𝑠 勁度矩陣 𝑛 結點上的未知數矩陣 聯立方程組的常數項 𝜌 流體密度 p 孔隙水壓 Uij 流速向量( u, v T) u x 軸方向流速
v y 軸方向(高度)流速 μ 黏滯係數 p 孔隙水壓 Fvol 體積力 𝛾̇ 剪應變率(=𝑑𝑢 𝑑𝑦⁄ ) 𝛬𝑖 流體變形應變張量 D̅ 流體在某單一方向的剪應變率張量 γ𝑖𝑗̇ 剪應變率張量 φ 等位函數法的是函數或平滑函數 τv,ij 黏滯應力張量 單位矩陣 𝜏𝑐𝑡,𝑖 碰撞-紊流應力張量 𝑣𝑜𝑙,𝑥 x 方向體積力 𝑣𝑜𝑙,𝑦 y 方向體積力 𝜌𝑑 土石流總密度 g 重力加速度 h 土石流的流動深度 流動坡度角 𝑐 土石的 Coulomb 摩擦角
γr 重新初始化參數 ls 參數控制介面厚度，其數值取決於網格大小。 ls 兩相介面單位法向量 𝐶 體積濃度 𝜌𝑤 水的密度 𝐺𝑠 土石比重 b 滑移距離 法向量 0 初始壓力 L 上三角矩陣 U 下三角矩陣 M 預測矩陣 相對總誤差 𝑖 第 i 個參數的誤差評估 𝑎𝑡 絕對全誤差 N 自由度數量 𝑐,𝑡 某時間點土石流的質心位置(x 方向) 𝑦𝑐,𝑡 某時間點土石流的質心位置(y 方向) Uc,t→t+1 某時間間隔的質心流速
1.
第一章 緒論
1.1. 研究動機與目的
台灣本島位在歐亞板塊和菲律賓海板塊之交接處，地形上較為險 峻且地質脆弱、地震頻繁，又因台灣屬於海島型氣候，春季時有梅雨 季、夏和秋季時因為處於西太平洋颱風路徑上時有颱風和豪雨，處在 在這種強度降雨和高累積雨量的情況下常會造成表層土壤的流失或 崩塌地滑的危險，使得下游地區的生命財產受到極大威脅，並伴隨著 重大災害的發生，土石流即台灣最常見的災害之一。「土石流」在中 國大陸又稱泥石流，英文常用”debris flow”，”mud flow”，” debris avalanche” 等 ， 日 本 稱 謂 ” 山 津 波 ( や ま つ な み ， 羅 馬 拼 音 為”Yama-tsunami”) ”。 土石流基本上是由土石粒料和水混合並受重 力作用沿著傾斜的地表流動的流動體。 不論就人的居住地、農業生產、工商業活動及遊憩休閒等等，河 口三角洲常是人類活動的重要地點。人類往往利用溪流的谷口、或山 坡地居住。台灣地狹人稠，山高坡陡，河短流急，加上人口與工業經 濟的快速成長，本島的平地資源幾乎用到極限，被迫往山坡地開發， 一旦土石流發生時，會挾帶著大量土砂向下游運移，潛在對受影響區 域的生命財產造成極大的危害。民國 88 年發生了「九二一大地震」之後，因地震造成的地質鬆 動現象，讓土石流發生的機會更多且發生地區也增加不少。民國 79 年的歐菲颱風間花蓮縣銅門村土石流造成土石流，民國 93 年的敏督 利颱風間台中縣松鶴地區多處造成土石流，到民國 98 年的莫拉克颱 風間高雄縣小林村產生土石流而村子被淹沒，都是土石流威脅人類生 命財產的例子。 現今土石流之災害防治對策，可分為兩方面：一方面是針對土石 潛勢流建構攔砂壩、護岸、堤防、山坡地水土保持工程等；另一方面 是透過土石流的預警係統、土石流扇狀地之危險區域畫定以及相關的 法令訂定等。然而只單純靠護岸、攔砂壩等工法來防制土石流的發生， 除了對於河川整體生態平衡造成破壞，也使得上游泥砂受到阻礙使得 下游地區失去了主要的泥砂來源，對於中、下游河床的變化以及河口 地區的地形變化會造成嚴重的影響。因此建立土石流預警系統，可說 是將土石流對人類的生命財產的傷害降到最低且又可兼顧生態的好 方法。而建立預警系統，需要能預測土石流流動行為，需探討土石流 的流動速度和性質以及會在何處堆積，本文因此嘗試藉 CAE 電腦輔 助數值分析軟體模擬土石流流體和流床的流變行為，以數值分析方法 模擬土石流的運移和堆積行為，潛在可供土石流預警之參考。
1.2. 研究目的
本研究採用 CAE 電腦輔助數值分析軟體 COMSOL Multiphysics 4.3b 的數值流體動力學模組 (CFD model)針對土石流進行模擬。 COMSOL Multiphysics 所用的數值模擬方法，運用有限元素法(Finite
elements method)可耦合模擬多重物理（可同時求解滿足數種不同物理 問題之控制方程式）之問題。基於文獻探討，本研究嘗試分別或兼而 考慮黏性泥型態與礫石流型態土石流的流變行為與特性，並探討不同 土石流流變性質對土石流流動特性和流動相關參數(堆積距離和流速 等)的影響。本研究為簡化問題，並不考慮土石流運移過程中流床之 侵蝕和土石流體體積濃度的變化，密度分布均勻且不隨運移過程而改 變。
1.3. 研究內容和流程
本論文分為六章，本章介紹研究的動機和目的。第二章介紹土石 流特性與種類、探討各種土石流的流變模式以研究土石流的流變行為， 蒐集以往土石流文獻所記錄的土石流物理和流變參數對土石流的流 動行為以探討參數的合理性，並收集土石流數值模擬的案例和模擬兩 相流雪崩的案例。第三章說明本研究所採取的方法，介紹兩相流土石 流流動機制的數值方法和工具，說明模型建立、模擬分析的流程和研 究規劃。第四章先驗證 Bovet et al.(2007)雪崩模型的可行性，再探討 土石流物理、流變參數對模擬結果的影響。第五章是虛擬案例模擬結 果展示，依研究規劃的案例，探討土石流物理參數對象流速和堆積距 離等流動特性參數的影響。最後第六章為本研究的結論和建議。研究 流程圖如圖 1-1 所示。圖 1-1 研究流程圖 研究開始 文獻回顧 土石流的流變模型和前 人土石流參數探討 模擬方法建立 土石流數值模型和分析 方法建立 模擬觀察和驗證 驗證模擬方法的合理性和觀 察參數對模擬的影響 結論與建議 本研究建構方法的檢討 和未來展望 修正 研究方法調整 模擬結果分析 模擬虛擬案例結果 和流動特性討論 結果不合理 結果合理
2.
第二章 文獻回顧
土石流是指泥、砂、礫石以及巨石和水混合並且和水混合而成的 混合物，並且受重力作用而在有傾斜地形下流動的混合流體。純液體 的水可以用牛頓流體來分析，但土石流卻未必適合牛頓流體的範疇， 其力學特性也比純水液體複雜。本章先回顧土石流的特性和流變模式， 並蒐集的土石流案例，亦參考回顧 Bovet(2007)等學者用 COMSOL 所 做的雪崩模擬案例。 (註：在本章節內除 2.9 節外，其他章節的座標系統是借用 FLO-2D 的 座標系統，x 表示長度座標、y 表示寬度座標、z 表示高度或深度的 座標)2.1. 土石流種類
不同的學者針對土石流也有不同的分法，以下是不同學者針對土 石流的分類方式： 1. 周必凡(1991) 依土石流流動特性及其組成將土石流分為兩 類： a. 黏性土石流(又稱結構型土石流)：土石流內含超過 3% 的細粒料，體積濃度 Cv介於 40%~80%。土石流全為層流流動， 且流動過程中土石流消散情形比較不會發生。運移停止開始堆 積時，堆積物會保持原狀。b. 稀性土石流(又稱紊流型土石流) ：土石流內細粒料含 量小於 3%，體積濃度平均 Cv在 40%以下，最低為 10%。這種 土石流的流動會形成高速的紊流流動，當流動漸減到停下來堆 積時，堆積物的分布會隨著粒徑而變。 2. 水土保持手冊(2005)依土石流內粒徑大小而分為三類： a. 礫石型土石流(Granular Flow)：細粒料粒徑 0.1mm 以 下含量小於總粒料含量的十分之一，同時粒徑 2.0mm 以上的粗 粒徑含量超過一半。 b. 泥流型土石流(Mud Flow)：細粒料(粒徑 0.1mm 以下) 含量超過一半，同時粒徑 2.0mm 以上含量低於總粒料含量的十 分之一。 c. 一般型土石流(Debris flow)：粒徑分布非礫石型和泥 流型的分類結果及稱呼之。 3. 熊剛 (1996)對中國大陸各地方所採集土石樣本進行粒徑篩 選後，從體積濃度、粒徑大小分布比例和粒徑大小來針對土 石流會出現的流動特性進行分類。分類如表 2-1。
表 2-1 不同土石流組成分類(熊剛，1996) 4. Coussot et al. (1996)藉由圖 2-1橢圓水砂含量循環橢圓形對於 水體和土石的相互組合有詳細的說明。由圖 2-1的上半圓來 看，幾乎都是顆粒較小土石材料為主，如果從土石塊體(solid) 往逆時針來看，隨著含水量增加，波地面的運動行為也從塊 體運動漸漸變成地滑(landslide)、泥流(mud flow)或土石流 (debris flow)、高含砂水流(hyperconcentrated flow )、河川水流
(stream flow)到一般水流(water flow)。若以圖 2-1的下半圓來
看的話，是以粗顆粒的土石材料為主，從水(water)往逆時針 來看，因為土石顆粒的增加，運動行為從原本的水流轉變成 二項流(two-phase flow)、顆粒流(granular flow)到土石崩落 (debris avalanche)。
圖 2-1 水砂含量和運動關係圖(Coussot et al,1996) 5. Takahashi (1977)依土石流發生機制整理成三種破壞類型： a. 滑動型土石流：起初，土石開始滑動時，除滑動面附 近的位移量較大之外，滑動體的變形量不會很大。然而如果土 石還是持續移動到一定程度，塊體內的變形量逐漸增大，使得 整體結構逐漸瓦解，如果在這此時土壤的含水量多到使土壤飽 和液化的話，就會形成土石流。 b. 潰壩型土石流：因堰塞湖形成的天然壩或人工所建造 的土石壩因偉分使強度減弱，當土石壩強度小到一定程度，壩 體開始崩潰土石和水混和在斜坡流動。 c. 沖蝕型土石流：假設一無限寬度邊坡，土層為均勻飽 和的土壤，受到降雨強度增加和逕流水位上升造成地表水深變 高，使地表某些土石隨者雨水沖蝕，沖蝕出的泥沙一旦量多和 水結合，如果受重力作用的話就會形成土石流。
圖 2-2 天然壩破壞引發的土石流(kuang et al.，1977)
6. 如圖 2-2 所示， kuang et al.(1988)針對天然壩破壞而引起的 土石流分出三種破壞機制：
a. 溢流沖蝕破壞(Erosion destruction due to overtopping)：
在強度大、均勻且透水性較低土壤的天然壩，一旦天然壩滲流 速度跟不上上游湖面水位上昇速率，天然壩發生溢流破壞，水 對天然壩表面的侵蝕力的增加，水和土石混和，而引發土石 流。
b. 滑動潰決(Sudden sliding collapse)：如果天然壩的土
壤結構很鬆散，且透水性強，而且有很明顯地滑動面的話，一 旦天然壩上游水位高於天然壩弱面，水就會進入弱面使摩擦力 降低，最後天然壩驟然潰決，天然壩土體汗水混和，而造成土 石流。 c. 漸進破壞(Retrogressive failure)：強度偏低、透水性高 的天然壩，天然壩後面的水會滲入天然壩，然後會滲流到下游 去。在滲流的同時，天然壩面向下游那一面的強度會減弱，並
個天然壩強度無法承受上游的水壓的話，就發生潰壩而產生土 石流。
7. Julien et al.(1991)引述 O’Brien et al.(1985)流變模式定義(此定 義會在 2.4.7 節會再作介紹)： τ 𝜏OB 𝜇B(𝑑𝑢𝑑 ) (𝑎𝑖𝜌𝑠 2𝑑 𝑝2 𝜌𝑚𝑙𝑚2 ) (𝑑𝑢𝑑 ) 2 (2-1) 將式(2-1)無因次化後： 𝜏∗ 1 1 𝑇 𝑑∗ 𝑎𝑖𝐷𝑣∗ (2-2) 𝜏∗ (2-3) 𝜏∗ 1 為賓漢流體(此會在 2.4.3 節會作介紹) 𝐷𝑣∗ 𝑑 (𝑑𝑢 𝑑 ) 顆粒碰撞力 黏滯力 (2-4) 𝑇𝑑∗ 𝑙 𝑎 𝑑 紊流影響 顆粒碰撞力 (2-5) 然後由式(2-4)和式(2-5)兩個無因次量的範圍量土石流分為三類， 如表 2-2。
表 2-2 由無因次量範圍區分土石流(Julien et al.，1991) 8. 以土石流內顆粒流運動過程所產生的力學機制來分類，這過 程中碰撞機制為主要影響因子（Savage(1984)；Savage and Hutter(1989)），因顆粒體之間或與底床互相碰撞、接觸，而 土石流內顆粒流內的顆粒動量出現消散和傳遞的現象，Iverson （1997）從大型試驗觀察得到的明顯影響影響土石流流動因子 （如碰撞、摩擦、黏滯性等）提出一系列無因次參數： 𝑁𝑆𝑎𝑣 ̇ 𝑑 𝑔ℎ collision orce riction orce (2-6) 𝑁𝐵𝑎𝑔 ̇𝑑 collision orce visco s orce (2-7) 𝑁𝑓𝑟𝑖𝑐 riction orce visco s orce (2-8) 式(2-6)到式(2-8)中： 𝜌𝑠、𝜌𝑓：分別為顆粒密度和流體密度 𝑑𝑝 ：顆粒平均直徑 𝛾̇：剪應變率(𝑑𝑢 𝑑𝑧⁄ ) λ： ( )1⁄3 1 1)顆粒體積分率(此在 2.4.5 節會介紹) 土石流種類 礫石流 (granular flow) 泥流 (mud flow) 含砂洪流 (mud flood) 無因次值範圍 𝐷𝑣 ∗> 400 和𝑇𝑑∗< 1 𝐷𝑣∗< 30 𝐷𝑣∗> 400 和𝑇𝑑∗> 1 主控力 顆粒碰撞力 降伏剪應力 和黏滯力 剪應力
從這三個無因次參數得知坡度變陡的話，土石滑動速率變動，剪 切速率、𝑁𝑆𝑎𝑣和𝑁𝐵𝑎𝑔也跟著變大。對於顆粒直徑大小，雖和摩擦數無 關，但顆粒太小則會被視為泥流的一部份(Iverson（1997）)。這三個 無因次參數的大小範圍，會主控土石流的力學因素，其因數大小的對 照表如表 2-3所示： 表 2-3 無因次參數大小參考值（Iverson，1997） 9.林美聆等人(2002)整理 Takahashi(1991)、五居有恆(1980)和水山 高久(1979)的土石流資料，並歸納出六種土石流型態： a. 邊坡滑動型土石流-邊坡上的土石不穩定，土體滑落到溪谷並 和水混和而形成土石流。 b. 位在溪谷兩側的堆積土，因暴雨或融雪(冰)而補充土體的水分 使得土體不穩定，最後水分足夠且坡地傾斜度夠大而形成土 石流。 c. 堰塞湖的天然壩體潰決而觸發的土石流 無因次參數 參考值 主控力學因素 參考文獻 𝑁𝑆𝑎𝑣 >0.1 碰撞控制 Savage et al. (1989) <0.1 摩擦控制 𝑁𝐵𝑎𝑔 >200 碰撞控制 Bagnold (1954), Iverson (1997) <40 黏性控制 𝑁𝑓𝑟𝑖𝑐 >1000 摩擦控制 Iverson (1997) <1000 黏性控制
d. 溪谷兩側邊破山腹經雨水沖刷，而不穩定坡址上的土石崩壞 且滑落至溪谷，而加大土石流的流量。 e. 火山噴發型土石流(有時稱為 lahar)-火山噴發通常伴隨著豪雨 的發生，其雨水將鬆軟的火山灰和噴發物混和成稠狀物，向 火山坡腳沖刷而下的土石流。 f. 在北歐、加州和阿拉斯加等地會出現特殊且粒徑相當小的流 黏土(Quick clay，又稱快泥，高緯度地區較常見)。流黏土遇 水、受擾動或震動，就容易出現液化，且流速很快，而形成 土石流。
2.2. 土石流的發展過程
圖 2-3 土石流發展示意圖(Johnson et al.，1984)
Johnson & Rodine(1984)把土石流的發展分六個階段：
一、在上游的堆積土體處於不穩定狀態，當土體開始發生破壞時， 會先從土體底部的剪動區進行進行破壞，如圖 2-3a。 二、剪動區漸近破壞滑動的同時，小的滑動土體也開始和水混合， 如圖 2-3b。 三、當土體和水混和並且滑動，土體軟化，後續土體越過前起 (snout)，流動將要發生，如圖 2-3c。
四、大部份的堆積土體材料和充份的水混合，且土體在運移過程 中受到翻滾、旋轉和擠壓等作用，土體內的凝聚力不斷減小， 並開始流動，如圖 2-3d。 五、土石流在流動路途中，會把流床的沿岸周圍土體也帶過來， 如圖 2-3e。 六、當土石流動速度開始減速甚至到停止時，土石流內的水分會 開始減少，但此時若後續土石粒料再加入，又和水再混合後， 土石流可能又開始流動，如圖 2-3f。 圖 2-4 土體受侵蝕而形成土石流(Tognacca et al.，1997) Tognacca et al. (1997)指出水先浸潤並滲入土體使土體強度減弱， 而後來流入的水會侵蝕因滲入水而強度弱化的土體，水混合成土石流 並流動，直到土石流因坡度、地面粗糙程度而減速到停下堆積，如圖
圖 2-5 山崩誘發型土石流的的流動過程(Takahashi，2001) Takahashi(2001)對於因山崩而誘發的土石流分為以下步驟： 一、在邊坡頂端部分，如果滲水發達，又地層蘊含很多地下水， 也因此出現滑動面，使土石開始滑落。 二、在土體滑動的過程中，土體內的水分集中在邊坡面較近的底 部，使底部含水量過大而出現液化，而變成土石流並在坡面 上流動。上層未液化的土塊的流速比下層已液化的土石漿體 快。 三、當土石流流道平緩的堆積區後，土石流減速停下，未液化的 土塊沉入已液化土石漿體內。
2.3. 土石流發生特性和運作
謝正倫(1991)和詹錢登(1998)針對土石流的研究和觀察，發現並 整理土石流有以下特性：1. 土石流發生地點區域內的崩積物厚度、水文特性、地質成 分和地形特性會和土石流開始流動的時間有相關性。 2. 土石流發生突然，很難預測發生的時間，而且土石流歷時 很短，最多只會歷經幾小時而已。 3. 土石流的流速剖面來看，底床的流速恆小於於表面流速。 4. 在土石流流動的橫切面來看，粗粒料通常集中在土石流前 端，後端含細粒料多於粗粒料，使土石流前端形狀隆起、後段形 狀偏向凹陷。如圖 2-6 所示。 5. 如果土石流的粒徑非均勻分布的話，流動會開始不穩定且 出現陣流波動的現象。就算前端土石流停下堆積，換段土石流會 因慣性作用而推向前端部分，使前端土石流又開始流動。 6. 當土石流流道坡度傾角很大且泥位又量體暴增，就算流道 有轉彎和障礙物，高流速土石流的強大的直進和衝擊力道，使土 石流動衝擊力大到可以跨越或沖毀障礙物和衝破流道的彎道。 7. 土石流的發展過程中依坡角的大小範圍而區分為發生區、 流動區、堆積區和停止區。發生區會在土石流流道波度 15~30 度 時，此區域易發生土石流。流動區在土石流流道坡度 6~15 度，在 這段土石流持續流動。堆積區發生於坡度 3~6 度的坡面，土石流
漸漸減速並開始堆積。最後流到坡度 3 度以下的坡面，即停止區，
流速漸漸歸零而停止並開始堆積堆積。如圖 2-7 和圖 2-8所示。
土石流一旦流到坡度緩且寬度寬的溪谷出口，會形成沖積扇。
圖 2-6 土石流流動剖面圖(詹錢登，1998)
圖 2-8 土石流運移過程圖(謝正倫，1991) Takahashi(1997)針對土石流運移過程中，依流道的坡度不同而出 現流動特性，如表 2-4所示： 表 2-4 土石流在運移過程中因坡度而所產生的流動特性 坡度 土石流流動特徵 𝛉 > 𝟐𝟎° 土石流開始啟動 𝟏𝟓°< 𝛉 < 𝟐𝟎° 土石流運移流動 𝟏𝟎°< 𝛉 < 𝟏𝟓° 流速開始遞減，並下坡運移 𝟑°< 𝛉 < 𝟏𝟎° 流速持續遞減，流動產生的動力波前開始靜止 𝟎°< 𝛉 < 𝟑° 堆集並趨於靜止
2.4. 土石流流變模式
土石流為固態的土石和液體的水混合的流體，固體和液體的交互 作用動量的傳遞過程中，包含漿體間的黏滯力與紊流強度及顆粒間的 摩擦力與碰撞力，而這些因素的強弱除了土石流體的組成成分有關， 而且與其流動速度有關，使得土石流流動情況較為複雜。要說明其流 動的情況，需要有適當的流變模式來探討土石流流體內部之剪應力𝜏 與剪應變率𝛾̇(=du/dz)間的關係，找出土石流的流變特性。到目前為止， 有很多學者提出很多流變模式。以下就介紹常用在描述土石流的流變 模式(先假設這些流變模型不隨時間改變)： 圖 2-9 常用的土石流變模型(Tatahashi,2007)2.4.1. 牛頓流體流變模式(Newtonian fluid model)
此模式是最早用在土石流的流變模式，如圖 2-9 的(a)直線，牛 頓流體是指流體的黏滯係數在恆溫下為一定值，且黏滯係數也不受濃 度和流速等因子影響。其流變方程式如下：
τ μ𝛾̇ (2-9) 式中： μ：黏滯係數， τ：剪應力， 𝛾̇ 𝑑𝑢 𝑑𝑧⁄ ：剪應變率。 此模式比較適合用於含泥砂量較低(Cv<20%)的含砂水體，分析泥 砂量較高的土石流較不適合，因無法描述土石流內顆粒間的凝聚力和 摩擦力。Einstein 在 1906 年首先提出非黏性球體在牛頓流體算出黏滯 係數的方法，這計算方法經證實在體積濃度(Cv)小於 2%的流體適用 牛頓流體表示。其濃度和體積濃度的關係式： μ μ0 1 .5CV (2-10) 式中： μ0：水的黏滯係數，約為10 3 (Pa*s) CV：體積濃度，其物理定義為： CV 1+ 1𝑤 (2-11) 𝜌𝑚：土石(泥)流密度 𝜌𝑤：水密度 𝜌𝑠：土石固態密度 𝐺𝑠：土石固態比重 𝑤𝑐(%)：含水量
Roscoe 在 1952 年提出兩個方程式，分別為水體含均勻顆粒(式 (2-12))和非均勻顆粒(式(2-13))：
μ μ0 1 1.35 CV 2. (2-12)
μ μ0 1 CV 2. (2-13)
Krieger and Dougherty 於 1959 年針對式(2-13)加入極限濃度的概 念，提出另一黏滯係數評估式： μ μ0 1 CV Cv 𝐵 (2-14) 式中： 𝐶𝑣𝑚：顆粒之最大極限濃度，在泥沙顆粒為球形下，𝐶𝑣𝑚約為 0.74， ：B 值約 2.5。 Bagnold(1954)依實驗結果，假設流體呈現黏性的時候，顆粒的存 在對黏滯性的影響如下式所示： μ . 5μ0 1. (2-15) 其中 為線性濃度，其和含砂濃度的關係為： λ 1 1 (2-16) 式中𝐶𝑣𝑚為顆粒之最大極限濃度。
2.4.2. 摩擦流變模式(Friction model) 顆粒間摩擦效性強的顆粒流之流域比較接近於摩擦模式。假設顆 粒流內孔隙水含量相當少，在此流域中，顆粒流裡面顆粒運動速度較 為緩慢而顆粒間緊緊相距，使得顆粒和顆粒之間所產生滑動摩擦成為 顆粒流間動量轉換的方式，因此這流域則稱為摩擦域(friction regime)。 依據 Coulomb 摩擦理論，在摩擦流域的任一斷面內的剪應力和正向 力用下面的關係式表示： τ σ𝑛 φ𝑐 (2-17) 式中： τ：剪應力， σ𝑛：顆粒間相互接觸在剪力面產生的內正向力(internal normal force)， c：凝聚力(cohesive strength)， φ𝑐：顆粒內摩擦角，一般會隨體積濃度 Cv而變，但通常視為定 值
2.4.3. 賓漢流變模式(Bingham model)
賓漢模式是由 Bingham and Green 在 1920 所提出的流變關係式， 如圖 2-9 的(b)直線，又稱黏塑性模式（Visco-Plastic model），此模 式中含有一賓漢降伏應力(Bingham Yield stress)，流體所受的剪應力 需要大於賓漢降伏應力才會出現剪應變率，其方程式如下： τ τB μB𝛾̇ (2-18) 式中τB為賓漢降伏應力，μB為賓漢黏滯係數，此兩參數會受到土石流 的含砂量、粒徑分布、礦物性質、溫度、酸鹼度等影響(Johnson, 1970、 O’Brien et al., 1988)。 賓漢降伏應力τB的形成主要源自黏性細顆粒泥砂的絮凝結構和 泥砂顆粒間的摩擦阻力。以泥砂漿體來說，細顆粒泥砂的絮凝結構是 形成降伏應力τB主要原因。Johnson（1965）經由實驗和現場觀測土 石流將降伏應力用摩擦模式表示，摩擦模式凝聚力和內正向力，其中 摩 擦 力 正 比 於 作 用 在 剪 力 面 的 正 向 應 力 。 Johnson 建 議 之 Coulomb-Viscous model 如下： τ σ𝑛 φ𝑐 μB𝛾̇ (2-19) 等號右邊前兩項為降伏應力的部份，如果土石流開始流動，須先 超過此降伏應力。τB及μB兩參數並非定值，會隨泥砂體積濃度、粒子 形狀、黏土性質、粒徑分布與溫度等情況而改變。
Johnson et al. (1984)提出下列在現地估計方法，以估計賓漢降伏 應力： a. 臨界坡度法： 假設邊坡的寬度式無限大，如所示，用土石流停止流動時的 深度來求得賓漢降伏應力： τ𝐵 Tcρ (2-20) Tc：臨界流動深度(m) ：堆積傾斜坡度( °) b. 特殊大石塊尺寸(於 1917 年在 surprise canyon 土石流案例求 得的)： τ𝐵 0. 19𝑙𝑚𝑖𝑛 ρb 𝑛𝑣ρ (2-21) ρb：特殊大石塊單位重( ⁄ 3 ) 𝑙𝑚𝑖𝑛：特殊大石塊在 x 軸、y 軸和 z 軸中最短的長度 𝑛𝑣：特殊大石塊被水淹沒的體積 特殊大石塊總體積
c. 半橢圓形截面土石流流道估算： 圖 2-10 將半橢圓截面的流道簡化成矩形截面的流道 (Johnson et al，1984)(圖中的 d 為流動深度 h) Johnson(1970)認為，土石流流道比較常是半橢圓形截面，通常土 石流流動時離流道壁面最近的土石流會先降伏成流體流動；離較遠的 流道壁面較遠的，因流動剪力為大於降伏剪應力而不會產生流動，形 成塞流(plug flow)，會被下面流動區域搬運運移。圖 2-10 表示半橢圓 形截面流道(圖 2-10(a))可以簡化成矩形截面流道(圖 2-10(b))，其原 因是因土石流降伏強度，而矩形截面流道離尖角部分較近的部分不會 流動，而形成”未流動死水區(dead regions)” (Johnson et al, 1984)。 其兩者的截面，土石流流速分布幾乎差不多。假設土石流流道截面為
線對稱的半橢圓截面，流道寬度 W 約為深度 h 的兩倍，塞流寬度W𝑝， 其賓漢降伏剪力強度為τ𝐵為： τ𝐵 W𝑝⁄ 𝜌4 𝑚𝑠 𝑛 (2-22) 𝜌𝑚：泥漿密度 ：流道傾斜坡度( °) 然後再假設 𝑝,𝑚𝑎𝑥為土石流中最大表面流速(同時也為塞流的表面流 速)，其黏滯係數μB為： μB ( 𝑣 , ) ( ) 1 2 (2-23) 在 Pierson(1986)在 3 號監測站監測土石流流道，其流道寬度 W 不等於深度 h 的兩倍，泥漿密度𝜌𝑚為 .0 𝑚3，流道傾斜坡度 為 7 度，其賓漢降伏剪力強度為τ𝐵為： τ𝐵 (W𝑝⁄ )𝜌𝑚𝑠 𝑛 𝑊 ⁄ ⁄ 2 1 (2-24)
在 Pia R Tecca et al.(2003)中有提到在 1998 年，有用到此現地估 測方法來監測義大利 Eastern Dolomites 地區內 Acquabona 溪流流域的 土石流，也因此求得賓漢流變參數，其結果如表 2-5 所示：
表 2-5 1998 年義大利 Acquabona 溪流流域土石流監測結果 (Pia R Tecca et al.,2003)
賓漢黏滯係數(μB)和牛頓流體的黏滯係數相當，賓漢黏滯係數可
用牛頓流體的黏滯係數的計算方式。然而大部分的賓漢黏滯係數是從 室內試驗方法所得到的。一般取土石流漿體部份，依試驗儀器限制決 定取樣的最大粒徑。較常用的為旋轉式流變計：試驗時量測轉矩 (moment)與轉速，求取管壁剪力與外視剪應變率(apparent shear strain rate)，經由迴歸分析求得賓漢降伏應力與賓漢黏滯係數。 雖然賓漢模型可描述泥流的流動行為，但對於粗粒料的流動行為 並無法完整描述，因賓漢模型不包含粗粒料碰撞而產生的離散力和紊 流力。然而，如果假設粗粒料碰撞而產生的離散力小到可忽略，也可 以用賓漢模型來逼近。王志賢(2000)用兩個總體積濃度(即包含所有粒 徑粒料的體積濃度，也就是 Cv)相同實驗樣本。一個是純細粒料泥漿 體，另一個是含粗、細粒料混和的泥漿體。兩者都用 Brookfield DV-III
水平旋轉流變計，測量出各自的賓漢降伏應力和賓漢黏滯係數。然後 針對含粗、細粒料混和的泥漿體，以粒徑 2mm 為界限，分成純粗粒 料的體積濃度(Cvg)和純細粒料的體積濃度(Cvf)。然後將含粗、細粒料 混和的泥漿體實驗樣品所測出的降伏應力和黏滯係數減去純細粒料 泥漿體的部分而得到差值。減出來的差值用純粗粒料的體積濃度、粗 粒料粒徑和粗細粒料粒徑界限(即 2mm)進行迴歸出參數式。最後把原 本純細粒料的降伏應力和黏滯係數表示式加上粗粒料的部分由迴歸 得到的參數式，可得到含粗粒料賓漢流變模式表示式。整個實驗的概 念圖如圖 2-11 所示。 圖 2-11 含粗顆粒漿體流變特性之分析概念 註：粗顆粒濃度 Cvg= 粗顆粒體積 細顆粒體積+水體積+粗顆粒體積 (2-25) 細顆粒濃度 Cvf = 粗顆粒體積 細顆粒體積+水體積+粗顆粒體積 (2-26) 由王志賢(2000)所得到的實驗結果，粗粒料加入細粒料漿體皆可 略提高賓漢降伏應力和賓漢黏滯係數，但賓漢降伏應力提高的程度比 賓漢黏滯係數還多。
通常做土石流流槽實驗和數值模擬，如果要使用賓漢流變模式， 需針對賓漢降伏應力τB和賓漢黏滯係數μB對體積濃度 Cv的迴歸關係 式來表示土石流的材料性質。賓漢降伏應力τB和賓漢黏滯係數μB的迴 歸式會和體積濃度 Cv呈現自然對數關係，其迴歸式如下： 𝜏𝑦 𝛼1 (2-27) 𝜇B 𝛼2 (2-28) 對於式(2-28)，在圖 2-12由前人用賓漢模式為主作流槽實驗作出 體積濃度𝐶 對應賓漢黏滯係數𝜇B的實驗結果圖，將圖 2-12上面各𝐶 值對應的𝜇B作迴歸，其迴歸式會與(2-28)差不多。此迴歸計算方式在 之後 O’Brien et al.(1985)所提出的二項式流變模式也用此方法估算體 積濃度對降伏應力和黏滯係數的關係。 圖 2-12 以賓漢流體為基底，前人實驗賓漢黏滯係數和體積濃度的關係 (間 接摘自 Takahashi(2007))
Dent and Lang (1983)模擬雪崩時針對藉由雙折線法(雙黏性)來考 慮賓漢流變模型中的賓漢降伏應力τB： 當賓漢黏塑性物體所受的剪力未達τB時，且剪應變率非常小(假 設其剪應變率為𝛾𝑙𝑠𝑟，假設𝛾𝑙𝑠𝑟̇ < 10 1 𝑠 1 )時，其處在固體狀態的黏 塑性物體內的黏滯係數非常大(假設其黏滯係數為𝜇𝑙𝑠𝑟，假設 𝜇𝑙𝑠𝑟 > 10 𝑃𝑎 ∗ 𝑠 )，大到黏塑性物體無法流動。到達一定的剪應變 率時，所受的剪力約大於降伏應力τB，其原本在固體狀態的黏塑性物 體開始液化。其說法可以解釋泥流型土石流受剪力才開始流動的物理 現象，其流變模型稱為"賓漢雙黏性流變模型(Bi-viscous Bingham model)”，如圖 2-13。陳孟志(2011)和柯昱明(2012)分別模擬海嘯沖刷 坑的發展和高濃度泥沙沖刷行為，也採用此賓漢雙黏性流變模型。 圖 2-13 賓漢模型(實線部分)及賓漢雙黏性模型(虛線部分)的比較(Dentand et al.，1986)
Shao & Lo(2003)曾用過 Barnes H.A. et al. (1989) 所提出的 Cross 流變模式，以等效黏滯係數，將賓漢流體簡化成牛頓流體的形式。 Cross 流變模式的式子如下： 𝑘𝛾̇ 𝑚 (2-29) 𝜇lsr和𝜇 sr分別則表示在極小和極大剪應變率下的黏滯係數。𝜇𝑒𝑓𝑓 為等效黏滯係數，其定義等於 總剪應力 總剪應變率 ̇。k 和 m 皆為常數。 在式(2-18)的賓漢模式同除於剪應變率𝛾̇所得到的等效黏滯係數： 𝜇𝑒𝑓𝑓 ̇ μB (2-30) 賓漢流體屬於含降伏應力的線性流體，則令式(2-29)的 m 值為 1， 由式(2-29)所得到的等效黏滯係數為： 𝜇𝑒𝑓𝑓 + ̇ 1+ ̇ (2-31) 為了使式(2-30)和式 (2-31)相近，則假設式 (2-31)的參數： ：常數 𝑘 𝜇 sr：極高剪應變率的黏滯係數，在賓漢模式，𝜇 sr等於μB
𝜇lsr：極低剪應變率的黏滯係數。Hammad and Vradis(1994)數值
模擬賓漢流，當𝜇lsr等於103μ
B，整個 Cross 等效流變模型連續性較佳，
賓漢模型和 Cross 等效賓漢模型的比較如圖 2-14 所示：
圖 2-14 賓漢模型和 Cross 等效賓漢模型的比較 (Shao et al.，2003)
2.4.4. 赫許布雷流變模式(Herschel-Bulkley model)
Herschel-Bulkley (1925) 提出此流變模式，如圖 2-9 的(c)曲線(𝜗 <1 的情形)，屬於具降伏應力之冪次方模式(yield power law model)的 一種。 τ τHB μHB 𝛾̇ (2-32) 由上式可知，當剪力大於屈服應力τHB且𝜗不為 1 時，流型曲線不 再是直線而是曲線。而當𝜗等於 1 時，就是賓漢模式，此一模式融合 了賓漢模式與冪定理模式，因此應用上相當廣泛。在土石流研究上， 亦有許多研究者採用此一流變模式描述土石流體的流變特性。 此模式具有三個未知參數而分析起來較為複雜，而優點是應用較 為廣泛，且適合具有降伏應力且流變為非線性的土石流漿體，為了解 決參數多不易分析的問題，法國學者 Coussot et al.(1995)認為泥漿體 中的𝜗值是一定值，經實驗該值約為 1/3，但其它的研究者所得的𝜗值 卻不一定與該值接近，王志賢 (2007)在表 2-6 整理出學者們針對 Herschel-Bulkly 模式𝜗值範圍及相關條件。此外，Herschel-Bulkly 模 式與冪定理模式有一樣的問題，即當𝜗值不等於「1」時，μHB不可用 黏滯係數(Viscosity)來定義，須改用流體稠度指標(Consistency index) 或流體行為指標(Fluid behavior index)來定義，且單位也不能用黏滯係 數的單位，因其單位也會帶有冪次方。其流變模式用於土石流最主要 的目的，是在於用賓漢模式來解釋土石流的流變關係時，容易出現高
估降伏應力的情形。赫許布雷模式在𝜗小於 1 時，剪應變速率較小的
情況下(<10𝑠 1)，剪應變率下降到為 0，剪應力下降比賓漢流體還多。
使赫許布雷模式的降伏應力τHB比賓漢模式的降伏應力τB還要小，如
圖 2-15所示。Coussot(1994)指出，土石流粒徑小於 40μm 的粒料佔粒
表 2-6 各研究者所提出之 Herschel-Bulkley 模式𝜗值範圍 (王志賢，2007)
研究者 𝜗值範圍 泥砂體積濃度 剪切率範圍
Coussot and Piau (1994) 0.22~0.363 15% < Cv < 41% 0 < 𝑑𝑢
𝑑𝑧 < 50𝑠 1 Coussot and Piau
(1995a) 0.33 38% < Cv < 43% 0 < 𝑑𝑢
𝑑𝑧 < 100𝑠 1 Coussot and Piau
(1995b) 0.15~0.38 61% < Cv < 80% 0 < 𝑑𝑢 𝑑𝑧 < 100𝑠 1 Coussot et al. (1998) 0.33 30% < Cv < 84% 0 < 𝑑𝑢 𝑑𝑧 < 50𝑠 1 Martino (2003) 1.66~1.85 61% < Cv < 80% 0 < 𝑑𝑢 𝑑𝑧 < 100𝑠 1 Malet et al. (2003) 0.25 ~ 0.4 30% < Cv < 60% 0 < 𝑑𝑢 𝑑𝑧 < 50𝑠 1 Schatzmann et al. (2003) 0.33 ~ 0.5 22.5% < Cv < 30% 0 < 𝑑𝑢 𝑑𝑧 < 100𝑠 1
圖 2-15 在不同剪切率的 Bingham 流體直線和 Herschel-Bulkley 流體曲線比較(摘 自王志賢，2007)
2.4.5. 碰撞/膨脹流變模式(Collision/Dilatant model)
碰撞模式主要適用於顆粒間碰撞效應較為明顯之情況，為冪次方 模式( power law model) 𝜗值為 2 的情形，如圖 2-9 的(d)曲線。在此種 流域顆粒流動時保持一定的距離，顆粒變形速度也較快，顆粒和顆粒 間接觸的時間很短，因此顆粒間相互碰撞成為動量的轉換主要來源， 和摩擦模式的滑動摩擦不同。本流域稱為碰撞流域、完全動力流域 （fully dynamic regime）或慣性顆粒流域（grain-inertia regime），此 模式最早為 Bagnold（1954）所提出，又稱為膨脹模式（Dilatant model） 或離散模式（Dispersive model），形式如下： τ 𝑃𝑑 𝜇𝑑 𝛾̇ 2 (2-33) 其中𝜇𝑑參數和顆粒濃度及粒徑大小的分布有關，此參數則定義為 稠度指標。由 Bagnold(1954)的實驗結果顯示𝜇𝑑 𝑎𝑖 𝑠𝛼𝑖 𝜌𝑠𝑑𝑠 2， 其中λ 為線性濃度，其定義為： λ ( )1⁄3 1 1 (2-34) 𝐶 ：體積濃度 𝐶𝑏：顆粒靜止接觸時的最大體積濃度 Bagnold(1954)也提出顆粒離散應力(Dispersive stress)理論，其求 得的顆粒離散應力和其剪應力為： 𝑃𝑑 𝑎𝑖 𝑠𝛼𝑖 𝜌𝑠𝑑𝑝 2 𝛾̇ 2 (2-35) 𝜏𝑑 𝑃𝑑𝑡𝑎𝑛𝛼𝑖 𝑎𝑖𝑠 𝑛𝛼𝑖 𝜌𝑠𝑑𝑝 2 𝛾̇ 2 (2-36)
上兩式中，𝜌𝑠為泥漿中顆粒的密度，𝑑𝑝為顆粒的粒徑，𝛼𝑖為顆粒
離散應力的方向和垂直線(z 軸)之間的夾角，𝑎𝑖為顆粒運動係數。至
於顆粒運動係數𝑎𝑖，通常是用實驗結果所用的值，Bagnold 認為該值
為 0.042， Takahashi 認為該值為 0.3 ~ 0.5 之間，Julien and Lan (1991) 則認為該值為 0.0087。
Takahashi(1991)根據 Bagnold(1954)和 Takahashi(1978)定義礫石 型土石流(stony debris flow)，底床剪應力假設為阻抗應力(resistance
stress,指𝜏resist)，Bagnold(1954)定義顆粒碰撞所產生之剪應力(collision
stress,指𝜏𝑐𝑜𝑙𝑙𝑖)之和： τ 𝜏resist 𝜏𝑐𝑜𝑙𝑙𝑖 𝜌𝑠 𝜌𝑚 𝑔 𝐶 ∗ 𝑠 ∗ 𝑡𝑎𝑛 𝑐 𝜌𝑑𝑓𝑠𝑈 (2-37) 𝐶 ：體積濃度 h：流動深度 𝜌𝑚：土石流漿體部分的密度 𝜌𝑠：土石流顆粒固體部分的密度 ：流道坡度 𝜌𝑑：土石流整體密度(𝜌𝑑 𝐶 𝑔 ∗ 𝜌𝑠 (1 𝐶 𝑔) ∗ 𝜌𝑚) 𝑐：土石的摩擦角 𝑈：截面的平均流速(m/sec) 𝑓𝑠：顆粒碰撞對流速所產生的摩擦係數( 𝑠 2 𝑎 𝑠𝑖𝑛 𝑑 𝑑ℎ 2)
Tsujimoto et al. (1997)對於礫石型土石流內顆粒的碰撞、運動所
產生的碰撞應力(Collision stress, 𝜏𝑠)、顆粒之間相斥所產生的壓力
(particle repulsive pressure, 𝑐)和顆粒動能應力(Kinetic stress, 𝜏 )有進
行推導。其顆粒之間碰撞的示意圖如圖 2-16所示。 圖 2-16 高速顆粒流顆粒之間的應力關係 (Tsujimoto et al.，1997) Tsujimoto et al. (1997)利用兩顆粒相撞，以單一顆粒碰撞前後表 面每一點的動量差值做體積分後所產生的應力張量而得到碰撞剪應 力𝜏𝑐和顆粒之間相斥所產生的壓力 𝑐： 𝑐 𝜌𝑠 𝐶𝑣 2𝑔0 1 𝑒 (2-38) 𝜏𝑐 𝜌𝑠 𝐶𝑣 2𝑔 0 1 𝑒 𝑑𝑝√ 𝛾̇ (2-39)
