土石流模擬數值方法

FLO-2D 是常用來模擬洪水和土石流常用的工具，除前述 FLO-2D 之數值軟體以外，還有其他學者以控制土石流行為的控制方程組(一 組偏微分方程式)為基礎，建立數學模型，再以數值方法取求解。以 下列出一些 FLO-2D 以外之土石流數值模擬方法的例子。

(1).Takahashi(1987)以一維的土石流模型以有限差分法來模擬土 石流流道平面堆積區形成沖積扇。並且以相同的坡度、粒徑、密度和 體積濃度等物理量來做流槽實驗。數值模擬和實驗所得到的堆積面積 和深度分布在土石流開始流動後 10 秒、20 秒和堆積完畢後的比較如

圖 2-22所示。

比較數值模擬結果和流槽實驗結果，發現模擬出來沖積扇的深度 等位線較圓滑、實驗出來的沖積扇較不圓滑且較瘦長。可表示數值模 擬的沖積扇是偏向理想，堆積深度平均分布。

圖 2-22 沖積扇的模擬和實驗的結果比較(Takahashi,1987)

(2).對於因山崩而誘發形成的土石流，簡化為一維問題將斜坡上 山崩形成的土塊(Earth Block)分為三層：未飽和土塊層、飽和土塊層 和液化土石流層，如圖 2-23 所示：

圖 2-23 山崩土塊的數學模組(間接摘自於 Takahashi(2007))

Satofuka et al.(2003)針對三維山崩誘發型土石流，坡面上含液化 的土石流和土塊，將液化的土石流上的土塊切割成六邊形角柱，每塊 切成六邊形角柱的土塊分成不飽和區和飽和區。其土塊體動速度會受 到重力、六邊形角柱土塊體底部對液化區土石流和坡面之間的摩擦力、

六邊形角柱土塊體之間在滑動過程中的互(引)斥力和六邊形角柱土 塊體表面的剪力等影響，如圖 2-24 所示。

圖 2-24 山崩誘發型土石流的數值模型 (Satofuka et al. 2003)

圖 2-25 數值分析網格和六邊形角柱土塊邊界剪應力 的佔有比率 (Satofuka et al. 2003)

此模型所使用的數值計算方法為有限差分法 (finite difference method)，利用網格的節點範圍內來模擬土塊和液化後土石流的移 動。

此數值模型可以模擬土石流堆積時每一六角型區塊的在土石流 啟動、運移和堆積的分布，以預測土石流堆積時的堆積深度如何。但 此方法僅適用於體積濃度較大且是由山崩觸發而成的且含有未液化 土塊的土石流。Satofuka(2004)針對 2003 年日本熊本縣水俣市寶川內 地區因豪雨使集川河水氾濫而造成的土石流，運用此數值模擬的方法 來重現其災害。

圖 2-26 水俣市寶川內地區某邊坡破壞前後的形狀(Satofuka,2004)

圖 2-27 依邊坡破壞深度，所建模出來的土塊初始深度分布圖(Satofuka,2004)

圖 2-28 深度和土塊高度分布之模擬結果(Satofuka,2004)

(3).Laigle(1997)採用赫許布雷模式(Herschel-Bulkly model，簡稱 HB 模式)模擬土石流的黏塑性行為。其控制方程式用陡坡淺水方程式 (shallow water equation)，計算方法是用正方形網格用有限體積法 (Finite volume method)來做計算，輸入參數為 和^μ

(即 ^降伏應力

土石流密度和

黏滯係數

降伏應力)。用赫許布雷模式比較適用於細粒料較多、體積濃度較高且 不含粗粒料的土石流流體。Rickenmann et al. (2006a)使用 HB 模式模 擬 1987 年在瑞士東南部格勞賓登州(Canton of Grisons)的 Varunasch 集水區所發生的土石流案例，其集水區的面積為 6.5𝑘𝑚²，其流變參 數是用現地調察，如圖 2-29 所示。在圖 2-29 上，輸入參數分別為 τ_HB⁄𝜌_𝑚 1. 𝑚²⁄𝑠²和μ_HB⁄τ_HB 0.3𝑠^{1 3⁄} ， 左 半 部 (a) 圖 表 示 為

=600𝑚³⁄ 的模擬結果、右半部(b)圖表示為 𝑠 =900𝑚³⁄ 的模擬結果，𝑠 紫色線段則表示真實土石流堆積擴散區域。

圖 2-29 Rickenmann et al. (2006a)用 HB 流變模式模擬瑞士 Varunasch 集水區土石 流的堆積厚度。

(4).黃名村等人(2002)對神木村現地做數值分析使用 Julien &

Lan(1991)的二項式流變模式，假設土石流流道深度很大，強剪層可 因此忽略(含底床充實效應，即剪應力張量元素的絕對值大於降伏應 力，其算式為|𝜏_𝑖 | > 𝜏_𝑦，又名邊界層)，取土石流弱剪層(不含底床充 實效應，即剪應變張量元素為零，其算式為|𝛾_𝑖 | > 0，又稱為塞流層)，

依動量連續方程式和運動邊界條件寫出控制方程式。再以有限差分法 將控制方程式離散成差分式，最後分別以黃名村(1996)所提出的一維

水平渠道解析和特徵線法解出一維數值解(圖 2-30)和二維數值解(圖 2-31 和圖 2-32)。所求出的數值解和現地情況的誤差在 0.5%以內。

圖 2-30 一維土石流之流動深度數值解與解析解之比較 (黃名村等人，2002)

圖 2-31 以特徵線法計算土石流停止後之前緣位置與模擬結果之比較 (黃名村等人，2002)

圖 2-32 以特徵線法計算中央斷面之土石流流深分佈並與模擬結果作比較(黃名 村等人，2002)

(5). Rutschmann(1993) 在 瑞 士 聯 邦 科 技 研 究 所 (Swiss Federal Research Institute)的水利實驗室發展出針對流體力學的有限元素工具 箱(Finite Element Method TOOLbox，簡稱 FEMTOOL)，後來加入土 石流的流變模型和會隨著流速和深度而不同的乾濕過渡區域，就變成 土石流有限元素法模型(簡稱 DFEM model)，但此模式很少被用來做 土石流模擬。此模型是由隨意三角形網格組成，運算方式則採用上風 法(upwind scheme)來提升數值運算的穩定性，通常是在模擬開放流道，

此數值模型所使用的流變模型為 Voellmy 流變模型。而 Voellmy 流變 模型在 x 方向的摩擦坡降形式如下：

𝑆_𝑓𝑥 ^√

ℎ _𝑑 φ_𝑐 (2-74) 𝑞_𝑥：x 方向單位寬度的流量(𝑞_𝑥 𝑢)

𝑑：下游堆積區的傾斜度 φ_𝑐：顆粒內摩擦角

𝐶

Rickenmann et al. (2006a)用 DFEM-2D 模擬奧地利的地羅爾州 (Tyrol)東部的 Wartschenbach 集水區，其面積為 2.5𝑘𝑚²。模擬所使用 的網格 1576 個網格，其模擬結果如圖 2-33所示。在圖 2-33中，紅線 所包圍的區塊為含粗粒料所堆積的區域。綠線為土石流初始堆積時的 區塊。內部藍線為土石流再出現在堆積時，細粒料所向外堆積的區塊。

外部藍線為土石流堆積完後，土石流堆積和洪水靜止的區塊。

圖 2-33 用 DFEM-2D 模擬 Wartschenbach 發生土石流時堆積厚度。

輸入參數為尖端流量Q_𝑝=100𝑚³⁄ 、𝐶𝑠 _𝑝𝑐 11 𝑚^0. ⁄ 𝑠 (Rickenmann et al.，2006a)

(6). Martinez et al.(2008)以有限元素法做一維潰壩型土石流數值 分析，其數值模擬所使用的流變模型為賓漢 Cross 等效流變模型，其 模擬示意圖如下圖所示：

圖 2-34 一維斜面上的潰壩型土石流(Martinez et al.，2008)

在圖 2-34，傾斜角 為 11 度、壩體所使用的材料為高嶺土

(Kaolinte)、體積濃度𝐶 是 13.05%、壩體面積大小 A 為 24.7 𝑚²、壩 體在斜坡長度 ₀為 0.16m。其數值模型在模擬土石流運移的距離和實 驗值吻合。和解析解(其解析解並不精確)相比，數值解小於解析解，

但當模擬約到 0.7 秒，數值解開始大於解析解，Martinez et al.研判可 能和乾濕介面的數值處理有關，如圖 2-35所示。

圖 2-35 土石流運移距離解析解、實驗數據和數值解的比較 (Martinez et al.，2008)

(7). 聶 志 成 (2007) 使 用 有 限 體 積 法 計 算 流 體 力 學 模 擬 軟 體 FLUENT 6.2 選用 Realizable 紊流物理模組、標準壁面模式和以 流體體積法(Volume of Fluid，簡稱 VOF 模式)之兩項流模擬土石 流對壩體的衝擊力。依壩體形狀、土石流密度和黏滯係數的更換，比 較出壩體衝擊力的不同。

圖 2-36 土石流對直立壩的模擬畫面 (聶志成，2007)

圖 2-37 土石流對直立壩的流場模擬畫面 (聶志成，2007)

圖 2-38 直立壩所受衝擊力的分布 (聶志成，2007)

圖 2-39 直立壩在模擬時間內的的衝擊力的變化 (聶志成，2007)