COMSOL 雪崩案例回顧

本研究嘗試以 COMSOL 模擬土石流，參考 Bovet 等(2007)以該 軟體模擬瑞士阿爾卑斯山上的雪崩流動行為的文獻作為土石流數值 模擬比較基準。該模擬運用 COMSOL Multiphysics 的 Navier-Stokes 不可壓縮流模組(Incompressible Navier-Stokes Model，簡稱 NS 模組) 和等位函數法(Level Set Method)的勢函數守恆耦合而成的模組。(註：

在本節內座標系統，x 表示長度座標、y 表示高度或深度的座標，其 流速分別為 u 和 v)

不可壓縮流的控制方程式(governing equation)為：

ρ _t [μ ( U_ij ( U_ij^T))] ρ(U_ij )U_ij p F_vol (2-75)

Bovet 等(2007)使用等位函數法(Level Set Method)來決定雪崩和 空氣的流向並且以雪崩和空氣之間的交介面(interface)，定義方程式 φ(X,t)=0，最後而推導出兩項介質隨時間改變的方程式為：

t U_ij φ 0 (2-76)

上式中，φ為勢函數，將密度和黏滯係數大的雪崩的勢函數設為 φ > 0；將密度和黏滯係數大的雪崩的勢函數設為φ < 0，如圖 2-40 所示。對於隨著時間點所對應出的勢函數的傳遞，是用對流與擴散模 組(Convection and Diffusion Model，簡稱 CD 模組)對勢函數作暫態分 析，而得到每個時間點雪崩的流動情況。對流與擴散的控制方程式如 下：

𝑡𝑠

t D_𝑖𝑠𝑜 φ U_ij φ 0 (2-77) 式(2-77)中，令 _𝑡𝑠為 1、D_𝑖𝑠𝑜和 R 為 0 的話，式(2-77)會與式(2-76) 一樣。

圖 2-40 雪崩模組勢函數分布圖示(Bovet et al., 2007)

對於雪崩模擬用勢函數來決定雪崩和空氣的流向，不論是密度和 黏滯係數，都可以寫成一個勢函數的方程式來表示兩項流中兩種物質 的分布和流動情況。勢函數(φ)為 0 就是位於兩項流體的交介面，然 對於兩項流的流體分布方程式會出現不連續的現象。所以在雪崩和空 氣的交介面需要在用平滑步階方程式(smooth step function)做平滑效 果，使流體分布方程式在兩相流交介面上(勢函數φ 0時)，得以連 (Heaviside step function)，如式(2-80)表示：

{0, < 0 (Coulomb force)和黏滯力(viscous force)。當雪崩流到平緩的地面時，

摩擦所產生的庫倫力和水平流速 u 呈現(1 𝑒 ^𝑢)的倍率，也會對地面 產生摩擦的黏滯力會和水平流速 u 成正比。而由庫倫力和黏滯力所產

生的摩擦力和庫倫力和黏滯力的總和差在多成了一個摩擦係數𝜇_𝑐。在 這裡𝜇_𝑐 10 0.17 3。其中角度為 10 度的角度為雪崩的安息 角（angle of repose），也就是坡面的傾斜角度小於 10 度的話，坡面 上的雪崩的流速會漸慢到停止。

Bovet 等模擬的雪崩案例共有四個案例，本研究只探討前兩個案 例。第一個案例是假設雪崩為牛頓流，其黏滯係數為定值；第二個案 例是假設雪崩為會出現剪應變率遞減(shear-thinning)的現象的非牛頓 流，其黏滯係數會隨著雪崩流速而減少。圖 2-42 為第一個案例模擬 在 t=0 秒、1 秒、2 秒 4 秒和 12 秒時雪崩所流到的位置，t=12 秒顯示 雪崩流速幾乎停止且開始堆積。紅色為雪崩，藍色為空氣。圖 2-43

為第一個案例模擬在 t=1.2 秒時雪崩流速的分布。圖 2-44為第二個案 例模擬在 t=2 秒時雪崩的流速的黏滯係數大小分布。其模擬結果是用 COMSOL 3.5a 的 UMFPACK 求解器來求解(此求解器在 COMSOL 4.0 以後，因太占記憶體，而被刪除)。

圖 2-41 雪崩模組模擬所使用的網格(Bovet et al.,2007)

圖 2-42 牛頓流雪崩模擬結果 (Bovet et al.,2007)

圖 2-43 牛頓流雪崩模擬在 t=1.2 秒時流速的分布 (Bovet et al.,2007)

圖 2-44 非牛頓流雪崩模擬在 t=1.2 秒時雪崩黏滯係數大小的分布(Bovet et al.,2007)

綜合雪崩模擬的需求， 表 2-13 和表 2-14 分別列出 NS 模組和 CD 模組在雪崩模擬中所用到的參數：

表 2-13 NS 模組參數表格(Bovet et al.,2007)

註：1. abs(數值或參數)為括弧內數值的絕對值

2.flc2hs(數值或參數,步階解析尺度)為二維的海維賽階梯函數(2-D Heaviside step function)

3. Fx 內的 flc2hs(0.5-y,h_scale)，是指水平方向的摩擦力會集中在 y=0.5 軸的上下，

使摩擦力集中出現在地勢平緩的地方(也就是整個雪崩模組 y 值以 0.5 為中心而 y 值介於 0 到 1 之間的區塊)。

4.等向擴散係數的目的是為了增加 NS 模組計算流場擴散上的穩定性。

5.文獻下面附屬的表格物理參數定義有誤，eta1 和 rho1 分別為空氣的黏滯係數和 密度(數值較低)；eta2 和 rho2 分別為雪崩的黏滯係數和密度(數值較高)。需這樣

子區域(Subdomain):N-S 不可壓縮流 (Incompressible Navier-Stokes,簡稱 NS 模組)

參數符號 參數定義 COMSOL

eta eta1+H*(eta2-eta1) Pa*s

Fy 垂直體積力(

重力)

Fy gy*rho*(1-H) ³

Fx 水平體積力

(摩擦力)

Fx flc2hs(0.5-y,h_scale)*gy*

rho*0.1763*(1-exp(-u)+3*u)

定義，COMSOL 參數 eta(式(2-78))和 rho(式(2-79))才能成立。

6. COMSOL 3.5a 進行幾何建模是可用輸入幾何方程式的方式，ph_0 表示半圓形 方程式來表示雪崩的初始形狀和位置。但在 COMSOL4.0 以後，就無法以輸入幾 何方程式輸入幾何模型。

表 2-14 CD 模組參數表格(Bovet et al.,2007)

註：流線常數式增加 CD 模組勢函數傳遞和流場計算上的穩定性 子區域(Subdomain):熱傳擴散(Convection and Diffusion,簡稱 CD 模組)

參數定義 參數符號 COMSOL 參數或數值 單位

時間尺度 δ_𝑠𝑑 1 無

擴散常數 D_𝑖𝑠𝑜 0 無

反應速率 R 0 ³

水平速度 u u m/s

垂直速度 v v m/s

流線常數 δ_𝑠𝑑 0.25 無

勢函數 ph 無