旋轉水槽顆粒流試驗的流體形態

倪龍良，梁 龍，鄧檢良

（1.上海交通大學船舶海洋與建筑工程學院，上海 200240；2.浙江省工程勘察設計院集團有限公司，浙江寧波 315012）

0 引 言

泥石流對生命財產和生態環境等造成了嚴重危害，主要形式為破壞河道，也包括侵襲公路、鐵路等交通線路，例如2022 年貴州榕江站動車脫線事故就由泥石流侵襲動車路軌引發。因此，國內外學者對泥石流的運動機理展開了系統的研究，取得了一系列成果，文獻［1-2］中以Savage-Hunter（SH）模型［1］為理論基礎，以大型泥石流試驗結果為實驗基礎，建立了基于深度平均法的泥石流理論體系。在SH 模型中，假設溝底與顆粒流之間的摩擦服從庫侖摩擦定律，同時假設橫截面上的縱向壓力可以采用底壓與土壓力系數的乘積確定。SH 模型及其擴展形式已被許多水槽試驗［2-4］和野外觀測記錄證實［5］。Iverson［6］對沿粗糙河床穩定移動的理想化泥石流混合物進行分析，證實該模型可以預測實驗泥石流的速度和深度。另一方面，很多研究者提出對土壓力系數進行了改進，例如采用Rankine 公式的改進方法［7-8］；或者將土體壓力作為各向同性處理，從而直接將土壓力系數取為1［9-10］。

試驗研究方面，近年采用了旋轉水槽的方法再現泥石流取得了進展［11-13］，該方法的特征是用水槽的旋轉模擬無限長斜面，并使流體形成穩定的流動［14］。由于顆粒流（包括泥石流）流動問題的復雜性，目前尚無出現采用SH模型及其擴展形式對旋轉水槽顆粒流試驗穩定流動的模擬計算，這導致人們對旋轉水槽試驗結果的有效性產生懷疑。

本文以Savage-Hunter（SH）模型為基礎，通過修正土壓力系數及考慮溝底傾斜角變化和離心力影響建立旋轉水槽顆粒流流體的二維流動模型；同時用數值計算和試驗研究對比，重點分析旋轉水槽試驗的顆粒流深度、長度及龍頭龍尾位置變化等規律，驗證該模型和數值計算用于分析旋轉水槽顆粒流試驗的流體形態變化規律的有效性。

1 改進的Savage-Hunter模型

在旋轉水槽試驗中的顆粒流流體形態有特殊性：在穩定流動狀態下，其平均絕對速度為零；顆粒流長度與旋轉水槽轉速有關。為研究這一流動特性，采用SH模型并假設：11 忽略垂直于床層的速度分量，且運動過程中的質量和密度保持恒定，不考慮動量擴散效應；22 以深度平均流速代表實際速度，且溝底與顆粒流之間的摩擦角（簡稱界面摩擦角）與速率無關。

SH模型的基本方程是針對平直的溝底情況而提出的，本文首先將SH 模型推廣到靜止的圓弧形溝底情況，然后進一步推廣到旋轉的圓弧形溝底情況。

1.1 靜止的圓弧形溝底的SH模型

考慮固體顆粒狀材料沿旋轉水槽裝置的圓弧形底部輪廓的表面流動（見圖1），假設顆粒流不可壓縮流體，且水槽的旋轉速度為零，即圓弧形溝底靜止不動。

圖1 顆粒流沿靜止的圓弧形溝底運動

根據SH 模型，考慮溝底的傾角ζ 是隨位置發生變化的，并以圓弧形溝底可產生離心力導致溝底正應力加大［5］，可得圓弧形溝底的動量守恒方程和質量守恒方程：

式中：Kact為主動土壓力系數；Kpass為被動土壓力系數；θ為顆粒的內摩擦角。

由式（1）和（2）組成的偏微分方程組包含未知數深度＾h（x，t）和橫向平均速度。如果已知界面摩擦角δ、內摩擦角θ和溝底幾何形狀b（x），以及初始輪廓和速度分布，則就能確定顆粒流的深度＾h（x，t）和平均速度（x，t）。

1.2 旋轉的圓弧形溝底的SH模型

在旋轉水槽試驗中，圓弧形溝底處于旋轉運動，其動量守恒方程與前述靜止狀態下的圓弧水槽對應的動量守恒方程有兩點不同：

式中，uflume為水槽槽底處運動速度。需要說明的是，式（4）中的是平均化的絕對速度，并不是溝底處的絕對速度，而式（4）僅僅用于判斷摩擦力方向，對穩定流動狀態下的計算影響不大。

（2）土壓力系數的取值不同。旋轉水槽試驗結果表明，在水槽旋轉速度不高的條件下，顆粒流流動狀態以層流狀態為主，層面與溝底底面大致平行［15］。因此，可認為顆粒流內部發生塑性流動時的剪切破壞面與溝底底面大致平行。溝底底面方向上的正應力為pyy，該應力大致等于上述剪切破壞面上的應力；相對應地，pxx為與剪切破壞面垂直方向上的正應力，如圖2所示。對應的土壓力系數

圖2 沿著內部破壞的莫爾應力圓和破壞包絡線

式（5）和式（3）的主要區別在于對顆粒體內部的破壞面的判斷不同。式（3）的依據是固體顆粒在床層發生滑移時的莫爾圓分析結果，并認為這種破壞形式也存在于顆粒流內部。而式（5）是直接依據旋轉水槽試驗觀測結果：無論溝底的界面摩擦角多大，顆粒流內的大部分的剪切破壞面都與溝底底面大致平行，因此采用圖2 所示莫爾圓，這個莫爾圓與界面摩擦角δ 無關，且沒有主動土壓系數和被動土壓系數的區別。

依據式（4）和（5），且將控制方程式（1）和（2）的所有變量轉換為有量綱形式。得如下控制方程：

2 數值計算及結果

對式（6）采用拉格朗日法［1］求解，將顆粒流切分成一系列的“四邊形＋弓形”單元體（見圖3），網格單元的邊界旋轉水槽底邊切線垂直。在時間步長為n－1 時的網格邊界點定義為xj，n－1，邊界點處的速度定義為uj，n－1，單元塊的速度ˉui，n－1，邊界點對應的顆粒流深度定義為＾hi，n－1，其中下標i對應的是網格單元中心，j對應的是靠近左側網格i的網格單元邊界點j。

圖3 拉格朗日法網格單元符號的定義

假定在初始時刻已知邊界點速度uj，n－1、邊界點位置xj，n－1、hj，n－1的初始數值，將已知的顆粒塊的速度用插值法轉換為邊界點的速度uj，n－1后，計算經過時間步長dt后得到邊界單元點xj，n新位置（式（7）），然后確定單元中心點i上單元體的深度＾hi，n，最后，根據式（6）求解單元塊i在n時刻的速度：

偏微分的計算方法同文獻［1］。本文數值計算方法與SH模型計算方法最大的不同是在離散單元體中考慮傾角ζ的變化；另外，采用了式（4）和（5）判斷摩擦力方向和計算土壓力系數，即：

在圖3 坐標系下，旋轉水槽槽底表達為式（9），設初始時刻顆粒流的上表面為式（10），內摩擦角θ ＝39°，界面摩擦角δ ＝31°，時間步長dt＝2 ms，uflume＝0.06～0.10 m／s。采用數學工具軟件Mathcad 15 M050 版本編寫程序，求解h、L、u、β 及龍頭位置等參數，其中視摩擦角、龍頭等的定義如圖4 所示。數值模擬與試驗測試結果如圖5 和6 所示。

圖4 視摩擦角定義

圖5 uflume ＝0.08 m／s時顆粒流數值計算和實測結果

3 旋轉水槽試驗

3.1 試驗方法

如圖7所示為旋轉水槽裝置（或稱為滾筒試驗裝置［14］）制備顆粒流。系統由旋轉水槽（有機玻璃，內徑φ29 cm，槽底寬6 cm）、動力系統（步進電動機和傳動裝置）、攝影系統（50 f／s）三部分組成。試驗裝置的特征是水槽的轉速由計算機驅動模塊精確控制。

試驗顆粒材料是硅砂（密度2.65 g／cm3；50 μm孔徑硅砂）。采用動態試驗方法測量［16］，測得θ ＝39°，硅砂與有機玻璃之間δ ＝31°。

試驗初始狀態及旋轉水槽初始速度與數值計算的輸入初始狀態一致。水槽旋轉一段時間后，顆粒流會達到穩定狀態。在后處理階段，采用AE 軟件根據攝影記錄，提取顆粒流的位置坐標，每秒取5 點。需要說明的是，為保證流動的穩定，旋轉水槽速度uflume設定在0.06～0.1 m／s之間。

3.2 試驗結果與數值計算結果的比較

由于旋轉水槽試驗自身的特殊性，在水槽旋轉一段時間后，可以觀測到顆粒流穩定的流動狀態（見圖7（b））。在此狀態下，ˉu（x，t）≈0（見圖5（b）），這一現象表現為穩定流動狀態下的β 穩定不變。由圖6（a）可知，在不同的旋轉速度條件下，試驗測定的β與旋轉速度無關，為30°；計算得到的β 也與旋轉速度無關，為28°。計算結果與試驗結果接近，結果表明：＾h（x）在不同時刻存在一定差異［見圖5（a）］，但是差異很小，最大差異僅僅3 mm；而顆粒流上表面局部的小波動就超過1 mm［見圖7（b）］；顆粒流龍頭、中端和龍尾位置隨時間變化明顯［見圖5（c）］，這導致對深度＾h（x，t）的直接的測量和試驗驗證存在一定的困難。另一方面，在顆粒流體積不變的條件下，L（t）的變化可以大致反應＾h（x）最大值＾hmax的變化，因此，由圖6（b）中L（t）的實測值和計算值表明，L（t）在穩定流動狀態下，實測值和計算值分別為150 和165 mm，二者相差10%，屬于可接受的誤差范圍內。

圖6 顆粒流數值計算與實測變化曲線

圖7 旋轉水槽試驗

4 結 語

本文在Savage-Hunter（SH）模型基礎上研究旋轉水槽試驗中的顆粒流流動特性，建立了流體形態計算模型及相應數值計算方法，并進行數值計算與試驗結果對比分析，結果表明：

（1）考慮到顆粒流的流動特性，應對SH 模型中的土壓力系數進行修正。水槽旋轉速度不高的條件下，顆粒流流動狀態以層流狀態為主，層面與溝底底面大致平行，可以認為顆粒體（硅砂體）發生塑性流動時的剪切破壞面與溝底底面大致平行。因此土壓力系數須做相應修正。

（2）計算與試驗結果一致，顆粒流最終進入穩定流動狀態。穩定流動狀態下，視摩擦角的實測值和計算值分別為30°和28°；流體直線長度的實測值和計算值分別為150 mm 和165 mm。實測值和計算值相差不大。修正后的計算模型［式（6）］和相關的數值計算方法可用于研究旋轉水槽試驗中的硅砂顆粒流流體形態。

本文的研究為旋轉水槽試驗在泥石流試驗研究領域的進一步發展提供理論依據和數據支撐。但對于水槽旋轉速度對顆粒流運動特征的影響，缺少更多的試驗結果，尚需進一步研究積累水槽旋轉速度對顆粒流流體形態的影響。