旋轉(zhuǎn)水槽顆粒流試驗(yàn)的流體形態(tài)

倪龍良，梁 龍，鄧檢良

（1.上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海 200240；2.浙江省工程勘察設(shè)計(jì)院集團(tuán)有限公司，浙江寧波 315012）

0 引 言

泥石流對(duì)生命財(cái)產(chǎn)和生態(tài)環(huán)境等造成了嚴(yán)重危害，主要形式為破壞河道，也包括侵襲公路、鐵路等交通線路，例如2022 年貴州榕江站動(dòng)車脫線事故就由泥石流侵襲動(dòng)車路軌引發(fā)。因此，國(guó)內(nèi)外學(xué)者對(duì)泥石流的運(yùn)動(dòng)機(jī)理展開了系統(tǒng)的研究，取得了一系列成果，文獻(xiàn)［1-2］中以Savage-Hunter（SH）模型［1］為理論基礎(chǔ)，以大型泥石流試驗(yàn)結(jié)果為實(shí)驗(yàn)基礎(chǔ)，建立了基于深度平均法的泥石流理論體系。在SH 模型中，假設(shè)溝底與顆粒流之間的摩擦服從庫(kù)侖摩擦定律，同時(shí)假設(shè)橫截面上的縱向壓力可以采用底壓與土壓力系數(shù)的乘積確定。SH 模型及其擴(kuò)展形式已被許多水槽試驗(yàn)［2-4］和野外觀測(cè)記錄證實(shí)［5］。Iverson［6］對(duì)沿粗糙河床穩(wěn)定移動(dòng)的理想化泥石流混合物進(jìn)行分析，證實(shí)該模型可以預(yù)測(cè)實(shí)驗(yàn)?zāi)嗍鞯乃俣群蜕疃取Ａ硪环矫?，很多研究者提出?duì)土壓力系數(shù)進(jìn)行了改進(jìn)，例如采用Rankine 公式的改進(jìn)方法［7-8］；或者將土體壓力作為各向同性處理，從而直接將土壓力系數(shù)取為1［9-10］。

試驗(yàn)研究方面，近年采用了旋轉(zhuǎn)水槽的方法再現(xiàn)泥石流取得了進(jìn)展［11-13］，該方法的特征是用水槽的旋轉(zhuǎn)模擬無(wú)限長(zhǎng)斜面，并使流體形成穩(wěn)定的流動(dòng)［14］。由于顆粒流（包括泥石流）流動(dòng)問題的復(fù)雜性，目前尚無(wú)出現(xiàn)采用SH模型及其擴(kuò)展形式對(duì)旋轉(zhuǎn)水槽顆粒流試驗(yàn)穩(wěn)定流動(dòng)的模擬計(jì)算，這導(dǎo)致人們對(duì)旋轉(zhuǎn)水槽試驗(yàn)結(jié)果的有效性產(chǎn)生懷疑。

本文以Savage-Hunter（SH）模型為基礎(chǔ)，通過修正土壓力系數(shù)及考慮溝底傾斜角變化和離心力影響建立旋轉(zhuǎn)水槽顆粒流流體的二維流動(dòng)模型；同時(shí)用數(shù)值計(jì)算和試驗(yàn)研究對(duì)比，重點(diǎn)分析旋轉(zhuǎn)水槽試驗(yàn)的顆粒流深度、長(zhǎng)度及龍頭龍尾位置變化等規(guī)律，驗(yàn)證該模型和數(shù)值計(jì)算用于分析旋轉(zhuǎn)水槽顆粒流試驗(yàn)的流體形態(tài)變化規(guī)律的有效性。

1 改進(jìn)的Savage-Hunter模型

在旋轉(zhuǎn)水槽試驗(yàn)中的顆粒流流體形態(tài)有特殊性：在穩(wěn)定流動(dòng)狀態(tài)下，其平均絕對(duì)速度為零；顆粒流長(zhǎng)度與旋轉(zhuǎn)水槽轉(zhuǎn)速有關(guān)。為研究這一流動(dòng)特性，采用SH模型并假設(shè)：11 忽略垂直于床層的速度分量，且運(yùn)動(dòng)過程中的質(zhì)量和密度保持恒定，不考慮動(dòng)量擴(kuò)散效應(yīng)；22 以深度平均流速代表實(shí)際速度，且溝底與顆粒流之間的摩擦角（簡(jiǎn)稱界面摩擦角）與速率無(wú)關(guān)。

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

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

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

圖1 顆粒流沿靜止的圓弧形溝底運(yùn)動(dòng)

根據(jù)SH 模型，考慮溝底的傾角ζ 是隨位置發(fā)生變化的，并以圓弧形溝底可產(chǎn)生離心力導(dǎo)致溝底正應(yīng)力加大［5］，可得圓弧形溝底的動(dòng)量守恒方程和質(zhì)量守恒方程：

式中：Kact為主動(dòng)土壓力系數(shù)；Kpass為被動(dòng)土壓力系數(shù)；θ為顆粒的內(nèi)摩擦角。

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

1.2 旋轉(zhuǎn)的圓弧形溝底的SH模型

在旋轉(zhuǎn)水槽試驗(yàn)中，圓弧形溝底處于旋轉(zhuǎn)運(yùn)動(dòng)，其動(dòng)量守恒方程與前述靜止?fàn)顟B(tài)下的圓弧水槽對(duì)應(yīng)的動(dòng)量守恒方程有兩點(diǎn)不同：

式中，uflume為水槽槽底處運(yùn)動(dòng)速度。需要說明的是，式（4）中的是平均化的絕對(duì)速度，并不是溝底處的絕對(duì)速度，而式（4）僅僅用于判斷摩擦力方向，對(duì)穩(wěn)定流動(dòng)狀態(tài)下的計(jì)算影響不大。

（2）土壓力系數(shù)的取值不同。旋轉(zhuǎn)水槽試驗(yàn)結(jié)果表明，在水槽旋轉(zhuǎn)速度不高的條件下，顆粒流流動(dòng)狀態(tài)以層流狀態(tài)為主，層面與溝底底面大致平行［15］。因此，可認(rèn)為顆粒流內(nèi)部發(fā)生塑性流動(dòng)時(shí)的剪切破壞面與溝底底面大致平行。溝底底面方向上的正應(yīng)力為pyy，該應(yīng)力大致等于上述剪切破壞面上的應(yīng)力；相對(duì)應(yīng)地，pxx為與剪切破壞面垂直方向上的正應(yīng)力，如圖2所示。對(duì)應(yīng)的土壓力系數(shù)

圖2 沿著內(nèi)部破壞的莫爾應(yīng)力圓和破壞包絡(luò)線

式（5）和式（3）的主要區(qū)別在于對(duì)顆粒體內(nèi)部的破壞面的判斷不同。式（3）的依據(jù)是固體顆粒在床層發(fā)生滑移時(shí)的莫爾圓分析結(jié)果，并認(rèn)為這種破壞形式也存在于顆粒流內(nèi)部。而式（5）是直接依據(jù)旋轉(zhuǎn)水槽試驗(yàn)觀測(cè)結(jié)果：無(wú)論溝底的界面摩擦角多大，顆粒流內(nèi)的大部分的剪切破壞面都與溝底底面大致平行，因此采用圖2 所示莫爾圓，這個(gè)莫爾圓與界面摩擦角δ 無(wú)關(guān)，且沒有主動(dòng)土壓系數(shù)和被動(dòng)土壓系數(shù)的區(qū)別。

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

2 數(shù)值計(jì)算及結(jié)果

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

圖3 拉格朗日法網(wǎng)格單元符號(hào)的定義

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

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

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

圖4 視摩擦角定義

圖5 uflume ＝0.08 m／s時(shí)顆粒流數(shù)值計(jì)算和實(shí)測(cè)結(jié)果

3 旋轉(zhuǎn)水槽試驗(yàn)

3.1 試驗(yàn)方法

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

試驗(yàn)顆粒材料是硅砂（密度2.65 g／cm3；50 μm孔徑硅砂）。采用動(dòng)態(tài)試驗(yàn)方法測(cè)量［16］，測(cè)得θ ＝39°，硅砂與有機(jī)玻璃之間δ ＝31°。

試驗(yàn)初始狀態(tài)及旋轉(zhuǎn)水槽初始速度與數(shù)值計(jì)算的輸入初始狀態(tài)一致。水槽旋轉(zhuǎn)一段時(shí)間后，顆粒流會(huì)達(dá)到穩(wěn)定狀態(tài)。在后處理階段，采用AE 軟件根據(jù)攝影記錄，提取顆粒流的位置坐標(biāo)，每秒取5 點(diǎn)。需要說明的是，為保證流動(dòng)的穩(wěn)定，旋轉(zhuǎn)水槽速度uflume設(shè)定在0.06～0.1 m／s之間。

3.2 試驗(yàn)結(jié)果與數(shù)值計(jì)算結(jié)果的比較

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

圖6 顆粒流數(shù)值計(jì)算與實(shí)測(cè)變化曲線

圖7 旋轉(zhuǎn)水槽試驗(yàn)

4 結(jié) 語(yǔ)

本文在Savage-Hunter（SH）模型基礎(chǔ)上研究旋轉(zhuǎn)水槽試驗(yàn)中的顆粒流流動(dòng)特性，建立了流體形態(tài)計(jì)算模型及相應(yīng)數(shù)值計(jì)算方法，并進(jìn)行數(shù)值計(jì)算與試驗(yàn)結(jié)果對(duì)比分析，結(jié)果表明：

（1）考慮到顆粒流的流動(dòng)特性，應(yīng)對(duì)SH 模型中的土壓力系數(shù)進(jìn)行修正。水槽旋轉(zhuǎn)速度不高的條件下，顆粒流流動(dòng)狀態(tài)以層流狀態(tài)為主，層面與溝底底面大致平行，可以認(rèn)為顆粒體（硅砂體）發(fā)生塑性流動(dòng)時(shí)的剪切破壞面與溝底底面大致平行。因此土壓力系數(shù)須做相應(yīng)修正。

（2）計(jì)算與試驗(yàn)結(jié)果一致，顆粒流最終進(jìn)入穩(wěn)定流動(dòng)狀態(tài)。穩(wěn)定流動(dòng)狀態(tài)下，視摩擦角的實(shí)測(cè)值和計(jì)算值分別為30°和28°；流體直線長(zhǎng)度的實(shí)測(cè)值和計(jì)算值分別為150 mm 和165 mm。實(shí)測(cè)值和計(jì)算值相差不大。修正后的計(jì)算模型［式（6）］和相關(guān)的數(shù)值計(jì)算方法可用于研究旋轉(zhuǎn)水槽試驗(yàn)中的硅砂顆粒流流體形態(tài)。

本文的研究為旋轉(zhuǎn)水槽試驗(yàn)在泥石流試驗(yàn)研究領(lǐng)域的進(jìn)一步發(fā)展提供理論依據(jù)和數(shù)據(jù)支撐。但對(duì)于水槽旋轉(zhuǎn)速度對(duì)顆粒流運(yùn)動(dòng)特征的影響，缺少更多的試驗(yàn)結(jié)果，尚需進(jìn)一步研究積累水槽旋轉(zhuǎn)速度對(duì)顆粒流流體形態(tài)的影響。