基于光滑粒子元的水-沙兩相流沖刷數值仿真

張嶸釗 ，熊文 ?，劉川渟

（1.東南大學 交通學院，江蘇 南京 211189；2.同濟大學 土木工程學院，上海 200092）

水流侵蝕導致河床沖刷是水中建筑結構物常見現象，頻繁發生于橋梁基礎、防波堤、海底管線與堤壩下游等.2018 年，四川省堰塞湖泄洪導致下游河床劇烈沖刷，造成下游竹巴龍金沙江大橋沖毀，多座大橋下部基礎被掏空而發生傾斜［1］；余文疇等［2］研究發現，大通站長江干流1950—2002 年年平均徑流量中，每立方米水體平均含沙量達到47.2%，可見河床侵蝕作用十分顯著.因此，有必要針對河床沖刷演進過程進行研究，以指導水工建筑物抗水設計，妥善預防及處置沖刷災害，減少損失.

床底泥沙受到水流沖擊，會在表面產生切應力，當泥沙表面屈服時，細顆粒泥沙將以推移質和懸移質的形式向下游推移、輸送，并逐漸形成沖刷坑.沖刷誘因多種多樣：暴雨、泄洪、潰壩等短時間內產生巨大峰值流量皆會造成河床劇烈沖刷.沖刷本質皆為河床自由表面發生高度非線性變形、破碎以及水和泥沙相互夾帶.傳統數值模擬方法基于歐拉坐標系建立流域模型，通過數值求解輸沙方程和河床變形方程，獲得床面各處高程變化，最后采用動網格技術以及體積分數法（VOF）實現河床形態可視化［3］.然而，動網格技術對于網格質量具有較高要求，沖刷伴隨的河床大變形通常會造成網格畸變、計算不收斂等情況［4-5］；VOF 法通過捕捉水-沙交界面網格中材料體積分數分布情況實現河床形態可視化，因此精度受到交界面網格劃分密度限制，對河床劇烈變形的區域，往往不能精準捕捉［6］.可見，傳統網格法沖刷數值模型仍存在技術缺陷.

基于拉格朗日坐標系的光滑粒子流體動力學（Smoothed Particle Hydrodynamics，SPH）法，克服了傳統網格法的固有局限，與歐拉法相比，SPH 法粒子本身具有質量，無需額外計算就能保證質量守恒；模擬復雜自由表面流時，無需追蹤流體邊界及不同流體交界面，對于復雜液面具有優秀的求解能力，可適用于模擬泥沙沖刷、輸運等問題［7-11］.現階段，SPH 法模擬河床沖刷主要有兩種思路：其一，不進行土體粒子建模，使滿足條件的邊界粒子轉化為土粒子進入流域參與計算.轉化后土粒子視為離散元（DEM），不賦予材料本構，其運動方程受牛頓第二定律控制；土體與水體間考慮流固耦合作用力，土體粒子間考慮接觸、碰撞作用［8，11］.然而，此類模型計算相對復雜，與真實沖刷關聯較弱.其二，基于SPH多相流理論進行水體與土體建模［9-10］.土體粒子作為流體相可賦予材料本構，其變形與輸運具有真實物理意義.

目前，基于SPH 多相流理論的泥沙沖刷研究多集中于河床潰壩沖刷，土體多為非黏性泥沙顆粒（以下簡稱為泥沙），泥沙粒子流變特性通常采用非牛頓流體進行模擬［8-14］.Manenti等［12］在沖淤研究中，基于非牛頓流體理論構建了土體模型，并分別對Mohr-Coulomb 屈服理論和Shields 理論進行了比較分析，推薦后者作為泥沙模型的材料屈服強度；Fourtakas等［13］引入了廣義Herschel-Bulkley-Papanastasiou（HBP）模型進行土體建模，結合DP 強度準則確定材料屈服強度，并應用于水-沙兩相流沖刷模擬中，獲得了良好效果.

然而，現階段基于SPH 多相流的沖刷數值仿真方法仍處于發展階段：泥沙本構模型不統一，多數研究采用單一的土體強度準則，并集中于二維模型，對于三維泥沙沖刷模型及組合強度準則研究較少，對泥沙狀態轉變考慮不夠詳細［8-14］.由于SPH方法對于計算資源要求較高，三維模型需要耗費大量計算時間；同時沖刷過程中土體狀態變化急促、劇烈，單一準則難以完整表征泥沙由靜止到起動乃至輸運的全過程特性變化，最終造成河床沖刷形態與實際仍有差距.

針對上述技術難題，本文提出了基于SPH 多相流模型的河床沖刷數值仿真改進算法，主要包含三個步驟：①構建泥沙粒子；②根據侵蝕模型和起動準則判斷泥沙粒子狀態并分類轉化；③計算不同狀態下泥沙粒子流變特性.具體來說，基于開源軟件DualSPHysics 進行二次開發，編寫出泥沙沖刷計算模塊，利用開源CUDA 框架，進行GPU 計算加速；基于SPH 多相流理論與HBP 非牛頓流體模型，實現河床潰壩沖刷數值建模；引入Drucker-Prager 與Shields應力模型作為泥沙強度計算準則，考慮泥沙沖刷全過程中由沉積物到推移質再到懸移質的狀態轉變；最后，通過潰壩沖刷數值模擬，結合水槽試驗驗證模型準確性.

1 SPH基本理論

光滑粒子流體動力學（SPH）方法將連續流體離散為具有各種物理量的粒子，通過核函數實現粒子間相互作用，并按Navier-Stokes 方程進行運動控制.本節介紹了SPH 法控制方程及其離散形式，結合層流黏性應力與亞粒子（Sub-Particle Scale）紊流應力模型封閉N-S方程.

1.1 SPH法的積分插值理論

SPH 法根據相鄰粒子物理性質在每個粒子周圍進行局部積分，從而更新其物理量.相鄰粒子由核函數檢索，用W表示.

核函數W是關于光滑長度h的函數，其作用是將連續微分方程基于插值函數在特定點處進行積分，任意物理量連續形式用F可表示為

式中：F代表粒子的任意物理量；r和r'分別表示中心粒子和鄰域粒子位置矢量；h為光滑長度.將F在A點插值并離散化，則其離散形式表達為

式中：下標a表示中心粒子；下標b表示鄰域粒子；Δvb表示鄰域粒子體積.核函數選取參照Fourtakas等［13］的研究采用五次核函數［15］：

式中：q為粒子間距與平滑長度的比值.

1.2 控制方程的離散形式

N-S 方程中質量守恒方程與動量守恒方程分別表示為式（4）與式（5）：

式中：u為流速矢量；g為重力加速度矢量；P為壓力項；ρ為流體密度；Γ表示流體黏性力耗散項.

將質量守恒方程改寫為SPH離散形式：

式中：下標a表示中心粒子；下標b表示鄰域粒子；下標ab表示粒子a與b的差；?a為針對粒子a的哈密頓算子；m為粒子質量；ρa與ρb分別表示粒子a與粒子b的密度.

通過Favre 平均法對密度進行加權平均，可將SPS 項引入SPH 方法中，動量守恒方程改寫為SPH離散形式［16］：

式中：μ0為運動黏度；τij為SPS應力矢量.

2 基于SPH的兩相流數值模型

基于SPH 多相流模型的沖刷模擬，通常將泥沙視為非牛頓流體相，常見模型有Bingham 模型［17］、HB（Herschel-Bulkley）［18］模型、HBP（Herschel-Bulkley-Papanastasiou）模型［13］等.此類非牛頓模型在達到屈服強度前后，其表觀黏度變化顯著：材料所受剪切應力達到屈服強度前，黏度通常較大，流體表現為“難流通”，而在材料屈服后，表觀黏度通常急劇下降，逐漸轉變為“易流通”，因而適宜模擬沖刷起動前后的泥沙流變行為.然而，Bingham模型及HB模型在低剪切速率時其應力非連續，容易造成數值計算不收斂.本文采用HBP模型建立水-沙兩相流模型，將沉積相分為沉積物、推移質、懸移質三種狀態，賦予不同流變特性進行計算.

2.1 泥沙相非牛頓流體模型

HBP模型計算流體表觀黏度如下：

式中：μori為泥沙表觀黏度；μ為黏性系數；n為Herschel-Bulkley 冪指數參數，與剪切應力相關；m為Papanastasiou 參數，控制應力指數增長速率，使應力在未屈服區域保持小剪切速率，在已屈服區域保持線性增長；ⅡD為二階不變剪切應變率張量；τc為材料屈服強度.

Herschel-Bulkley 模型是廣義Bingham 模型，屈服后具有Bingham 模型流變特性，而Papanastasiou 模型適合描述剪切速率趨于0 時材料較大的表觀黏度，兩模型結合使得HBP 模型在模擬沖刷時穩定性更高，可連續表達土體從低剪切速率發展到高剪切速率的全過程［13］（圖1）.此外，當m=0，n=1時，該模型降階為牛頓流體模型，使得多相流計算中不同流相具有相同廣義模型.因而，在DualSPHysics 多相流耦合計算中，對于任意粒子核函數半徑內的被檢索粒子，僅需基于HBP 模型對所有被檢索粒子進行黏性力求解，從而規避了對土體粒子的分類檢索計算，簡化了計算框架，加快了計算效率［13］.

圖1 HBP模型應力曲線隨m和n的變化Fig.1 Diagram of the shear stress of HBP model developed with the m and n

2.2 泥沙侵蝕模型

為了模擬泥沙侵蝕與輸運過程，本文將泥沙劃分為三種狀態：沉積物、推移質、懸移質［19］.泥沙粒子狀態轉化與鄰域粒子相關，受其體積濃度、速度及所受剪切應力等因素影響，根據不同狀態，賦予不同黏度參與SPH計算.

2.2.1 泥沙屈服強度模型

Drucker-Prager（DP）屈服準則一般形式為：

式中：J2為二階不變剪切應力張量；p為飽和泥沙粒子上作用的靜水壓力；α與β由Mohr-Coulomb 屈服準則參數給出：

式中：φ為內摩擦角；c為土體黏聚力.

根據Fourtakas 等［13］的研究，J2可由土粒子所受水體剪切應力ταβ計算：

剪切應力ταβ可由二階不變剪切應變率張量ⅡD表述：

聯立式（11）、式（12）可以建立J2與ⅡD間聯系：

進而，|J2|=|αp-β|可作為HBP 模型中材料屈服強度限定條件.計算床面下靜止土體屈服強度時，引入Drucker-Prager（DP）屈服準則：

當沉積物粒子所受切應力大小不超過其屈服強度|τy|時，土體粒子由于高黏度特性幾乎不發生剪切變形，其表觀黏度μapp由|τy|代入式（8）進行計算.

對于河床表面與水接觸的土粒子，若滿足由沉積物轉化為推移質的條件，土體屈服強度采用Shields準則計算并替換，詳見2.2.2節.

2.2.2 泥沙起動判定

初始條件下，將泥沙層置于水層下方，當泥沙粒子所受切應力不超過臨界值時，泥沙粒子固定不動；達到臨界值時，泥沙由沉積物轉化為推移質，隨水流沿床面輸運.

泥沙臨界切應力由侵蝕準則計算.Manenti等［12］在沖淤研究中，對Mohr-Coulomb 屈服理論和Shields理論進行了比較分析.研究表明，Shields理論對于模型參數均具有較高敏感性；同時Shields 公式中物理量更易由試驗測得，相較Mohr-Coulomb 準則中應變率等參數，計算靈敏度更高，本文最終選擇Shields理論作為泥沙狀態轉化判斷準則.

轉化為推移質的泥沙粒子須位于河床表面，同時考慮到推移質粒子為飽和土的特性，其轉化條件設置以下兩個：

1）在泥沙粒子的鄰域粒子中，至少包括一個水粒子.

2）待轉化泥沙粒子實際質量應小于閾值質量.

當泥沙粒子同時滿足以上條件時，基于Shields侵蝕準則，計算水平床面臨界轉化剪切應力：

式中：τbcr，0為水平床面上的臨界剪切應力；d為特征粒徑，對于非均勻材料取中值粒徑d50；ρs為泥沙飽和密度；θcr為臨界Shields數，是泥沙粒子所受臨界起動力與最大抗力之比；ρw為水的密度；g為重力加速度.

θcr求解基于雷諾數Re［12］：

式中：R*=，v為水體運動黏度，u*為摩阻流速，按式（17）計算：

根據Prandtl 混合長度理論，假定沉積物表面存在厚度為δ的層流亞層，且流速u(z)按線性分布；在其上紊流層中，流速u(z)按對數分布：

式中：κ為馮·卡曼常數，取0.41；z為水粒子到泥沙-水交界處的垂直距離；z0為床底粗糙度，計算如下：

式中：ks為當量砂徑糙率，按經驗取泥沙中值粒徑.

由式（17）、式（19）可知，摩阻流速u*與δ和z0有關，因此需迭代求解.設定u*初始值u*0從ksu*v＜5的穩態流開始迭代，直至流速曲線u(z)在邊界層厚度δ收斂，即u(δ-)=u(δ+)，此時可求得摩阻流速u*，繼而求出水平床面的臨界剪切應力τbcr，0，以及斜坡修正臨界剪切應力τbcr.

水流沖刷的起動力τb由Einstein 對數流速分布公式［20］計算：

式中：d為粒子特征粒徑；Δu為泥沙粒子與附近水粒子的速度差，取中心泥沙粒子鄰域內所有水粒子進行積分插值：

若泥沙粒子所受剪切應力τb超過臨界剪切應力τbcr，則判斷泥沙粒子由沉積物轉變為推移質.轉化后粒子表觀黏度需結合HBP模型進行更新，將τbcr代入式（8）中代替屈服應力力τc，計算得到黏度μori.

2.2.3 泥沙上升與沉降判定

已屈服推移質粒子在滿足條件時會轉化為懸移質，被水流挾帶而發生輸運.當流速減小、大量懸移質粒子聚集時，懸移質再次沉降，轉化為推移質.

懸移質周圍泥沙粒子體積濃度Cv，i不超過臨界濃度Cvcr：

式中：i表示中心泥沙粒子；j表示鄰域粒子；臨界體積濃度Cvcr取0.3［12］，當滿足式（22）時，泥沙粒子可以作為牛頓流體處理，可以進一步考察流速條件.

引入Mastbergen 公式［21］，計算泥沙臨界上升流速與臨界沉降流速，如式（23）：

若u≥ulift，則判斷推移質轉化為懸移質，作為牛頓流體參與計算.其黏度采用Vand 膠體［22］方程計算：

若u≤uset，則判斷懸疑質沉降，轉化為推移質，黏度按推移質計算.

2.2.4 臨界起動切應力斜坡修正

對于任意坡度下泥沙起動分析，需對水平床面臨界剪切應力τbcr，0進行修正，得到適用于斜坡的臨界剪切應力τbcr.Van Rijn［23］將修正系數分解為縱向坡和橫向坡單獨影響；Chen 等［24］根據受力平衡推導出任意三維斜坡臨界起動切應力修正系數.本文在現有工作基礎上，考慮水流升力作用影響，進行臨界起動切應力斜坡修正.

假設泥沙粒子為一質點，處于斜坡床面上，其受力包括水下重力W（包含重力與浮力）、水流拖曳力FD、水流升力FL、抵抗起動力FC（本文中即土體的庫倫摩擦力）.斜坡、水流方向粒子受力如圖2所示.

圖2 斜坡床面向量示意圖Fig.2 Diagram of the vector of bed surface on the slope

圖2中，i、j、k為正交直角坐標系的單位向量；l、t所在平面為斜坡床面，其中l處于XOZ平面，與X軸夾角為α，t處于YOZ平面，與Y軸夾角為β；n為傾斜坡面單位法向量，與Z軸正方向夾角為γ，由幾何關系可知cosγ=水流拖曳力FD方向與流速u方向相同.根據Dey［25］的研究，假定水流升力FL方向與坡面法向量n一致，大小與拖曳力比值η為0.85.

水流拖曳力與水下重力沿坡面的分量合成泥沙起動力F：

式中：θx、θy、θz分別為單位向量e與X、Y、Z軸正方向的夾角.

抵抗起動力FC由庫倫摩擦力產生，與F方向相反，在極限狀態下，

式中：φ表示泥沙內摩擦角.當泥沙粒子達到臨界起動狀態，即滿足

由式（28），推導出|FD|表達式為：

由式（29）得到斜坡床面上的臨界起動切應力大小.對于水平床面，通過受力分析易得：

由式（29）、式（30）可得：

式中：k為臨界起動切應力斜坡修正系數：

水平面的臨界起動切應力經過斜坡修正為

2.3 泥沙模塊運行流程

本文沖刷模型中，初始狀態下泥沙粒子皆為沉積物，經由每個時間步進行處理，將依次完成沉積物向推移質轉化、推移質向懸移質轉化以及懸移質向推移質轉化.轉化完成后，分別存儲3 種狀態的泥沙粒子，進入SPH 求解步驟，直至計算結束，計算流程示意圖如圖3所示.

2.3.1 預處理模塊

獲取上一時間步中泥沙粒子狀態信息；迭代計算摩阻流速并應用Shields侵蝕準則得到臨界Shields數；按Mastbergen公式計算上升與沉降流速.

2.3.2 推移質模塊

沉積物粒子轉化為推移質粒子須滿足2.2.2 節所述起動條件.因此，需判斷其是否處于沙床表面，進而利用Einstein 對數流速分布公式計算該粒子所受流體切應力.若切應力超過屈服應力，則儲存為新的推移質粒子，并結合HBP 模型利用Shield 準則計算屈服強度，獲得推移質粒子表觀黏度.

2.3.3 懸移質模塊

推移質粒子轉化為懸移質粒子，須滿足2.2.3 節所述條件.由于水流裹挾懸移質輸送存在濃度閾值，首先判斷粒子鄰域內泥沙體積濃度，當濃度小于閾值時，判斷其速度大小，若速度大于上升流速，則判斷為懸移質粒子并依據Vand模型更新粒子黏度.

對于懸移質粒子，其速度若小于沉降流速，則又轉化為推移質粒子，返回推移質模塊.

2.3.4 沉積物模塊

泥沙粒子不起動，基于DP準則結合HBP模型更新其屈服強度.

3 潰壩沖刷數值仿真

SPH 法具有可以模擬流體大變形、復雜液面、不同液相間相互裹挾作用的優勢.在兩相流沖刷模型中，水-沙相互作用將根據控制方程自動計算，無需單獨求解輸運方程.本文數值仿真分為兩步：①基于DualSPHysics 開源代碼進行二次開發，完成泥沙狀態轉化功能；②進行潰壩水流沖刷模擬，結合水槽試驗結果驗證數值模型.

需要注意，本節重點對提出的泥沙模塊進行驗證，對于DualSPHysics 水體計算速度場、壓力場和湍流場精確性驗證，不在本文的討論范圍內.Sato 等［26］對多種場景下DualSPHysics 流域計算精度進行了驗證，限于篇幅，不贅述.

3.1 基于DualSPHysics的沖刷仿真計算流程

本文在DualSPHysics 開源代碼的基礎上，實現沉積物粒子起動判斷與類型轉化.二次開發模塊分為三部分執行：①初始化模塊；②判斷模塊；③黏度計算模塊.程序的運行流程如圖4所示.

圖4 程序運行流程Fig.4 Flow of program executing

針對不同問題與工程條件，本文二次開發模塊具有高度可定制性，可根據研究目標靈活地進行土層參數、泥沙起動理論變換.精確捕捉兩相流體交界面處粒子黏性與屈服特性，需盡可能提高粒子分辨率，即縮小粒子間距、增加粒子數量.為利用有限計算機資源提高SPH 模擬效率，本文基于CUDA 框架利用GPU 實現泥沙模塊加速計算，以解決三維模型對于計算資源要求高的問題.

3.2 數值模型參數設置

Dey［25］在研究黎曼波理論時引用了Louvain 潰壩侵蝕試驗，該試驗已被多名學者［9，11，27-28］用作沖刷模型評估標準，可為本文提供驗證結果.試驗中采用長度2.5 m、寬度0.1 m、高度0.35 m的透明水槽，在底部鋪滿厚度5～6 cm、等效直徑3.5 mm、密度1 540 kg/m3的PVC顆粒.試驗開始時，通過高速抬升水閘釋放高度為10 cm的水柱，并采用高速攝像機記錄沖刷過程.

本文三維數值模型水槽寬度設置為0.1 m，在厚度為6 cm、長度為200 cm 的泥沙層上建立深度為10 cm、長度為100 cm的水體.當模擬開始時，水體將迅速跌落，水-沙耦合作用將引起泥沙起動與輸送，從而發生潰壩侵蝕.

SPH 法中粒子流速與其設置間距緊密關聯，參照Zubeldia 等［19］的研究，通過設置4 組不同間距（0.016 m、0.008 m、0.004 m、0.002 m）模型組別與水槽試驗數據分別進行對比，從而確定適宜粒子間距，模型設置如圖5 所示.水槽空間布置同潰壩試驗相同，長2.5 m，水箱右側與水槽連接位置（X=0 m）設有0.02 m 開口，水槽底部鋪設泥沙，為避免泥沙隨水流發生輸運，其密度、屈服強度及黏度預設值偏大，可視為固定床面.水粒子動力黏度取1×10-6m2·s，密度為1 000 kg/m3，模型求解時長為8 s.各組別SPH 模型與試驗水頭在t=8 s 時前進位置繪制如圖6 所示.由圖6 可知，隨著粒子間距不斷減小，流體運動模擬結果逐漸逼近真實值.據此建立數值模型，粒子間距取dp=0.002 m，共產生3 087 630個粒子.

圖5 水槽模型示意圖Fig.5 Diagram of tank model

圖6 流體粒子間距結果對比Fig.6 Comparison of the distance between particles

泥沙與水的密度分別設置為ρs=1 540 kg/m3和ρw=1 000 kg/m3.泥沙相作為非牛頓流體，參照Fourtakas 等［13］的研究，HBP模型參數取m=100，n=1.8，動力黏度取0.001 m2/s，庫倫黏聚力c=0，內摩擦角φ=38°，以模擬物理試驗中無黏聚力PVC 材料特性；水作為牛頓流體，取m=0，n=1，動力黏度取1 × 10-6m2/s.沉積物轉化推移質質量閾值參照文獻［29］取1 250 kg/m3.

3.3 潰壩沖刷形態及誤差分析

為驗證SPH 法模擬河床沖刷等大變形及強非線性問題的精度優勢，本文采用Flow-3D建立相同空間布置網格計算模型，泥沙及水體材料參數設置與SPH模型相同，頂面為壓力邊界，其余為壁面，采用LES湍流模型，網格尺寸為0.002 m，模型如圖7所示.

圖7 網格模型示意圖Fig.7 Diagram of mesh-model

本次SPH 計算平臺為Intel i5-12600KF+NIVIDA 3080Ti，模擬時長為1 s，初始狀態模型如圖8 所示，求解時間為2 h 22 min.圖9、圖10 展示了SPH 模型、網格模型與Louvain 潰壩試驗在0.25 s、0.50 s、0.75 s、1.00 s 時剖面比較結果，其中：圖9（a）為Fraccarollo 物理試驗剖面圖；圖9（b）為本文模擬結果，其中藍色表示水體，紅色表示沉積物；圖9（c）為Flow-3D 網格模型模擬結果，其中藍色為水體，灰色為沉積物.由于篇幅限制，本文僅展示0.25 s 時刻試驗及各模型結果.圖10 展示了Louvain 試驗結果、Fourtakas 模擬結果、本文模擬結果以及網格模型模擬結果的水、沙表面高程對比圖，通過計算同一時刻試驗與模擬結果液面曲線豎直方向高程均方根誤差（ERMS），量化比較兩者輪廓線吻合程度.需要注意，在試驗與各數值模擬中，水與沉積物均置于三維水槽中，圖中展示僅為水槽側面視圖；同時，由于粒子總數量龐大造成其分辨率較高，部分懸疑質粒子緊鄰水-沙交界面，影響到真實河床表面高程觀測，因此實際河床沖刷后高程應以圖10為準.

圖8 數值模型示意圖Fig.8 Diagram of the numerical model

圖10 水面及沖刷地形輪廓線對比Fig.10 Comparison of free surface and scour morphology profile

由圖9 可知，在t=0.25 s 時刻，本文模型和Fourtakas 模型均能較好地再現Louvain 試驗中水頭行進位置與形狀；水-沙交界面處均較為準確地再現了最大沖刷坑深度與位置.從圖9（b）局部放大圖中可以看出，交界面處漂移泥沙粒子較多，且受到水流沖擊和裹挾作用于近床處發生懸移，河床真實沖刷位置應位于該類粒子下方.本文模型較好地模擬出沉積物泥沙受水體沖蝕發生狀態轉化，并在水流裹挾下推移、輸送，最終沉降、堆積.圖9（c）網格模型中河床在水流侵蝕下形成沖刷坑，未再現交界面處水-沙裹挾及水頭后方位置泥沙堆積行為，其沖刷坑位置更為靠后，水頭模擬誤差較大.圖10（a）反映出本文模型與Fourtakas 模型總體均方根誤差相近，但本文水頭行進位置與試驗結果更為相符，而網格模型水面線及沖刷地形均方根誤差均為SPH模型數倍以上.

在t=0.50 s 時刻［圖10（b）］，本文與Fourtakas 模型水頭與推移質行進長度均大于試驗結果，而本文所得最大沖刷深度相較于Fourtakas模型更為貼近試驗結果；網格模型中仍未出現泥沙堆積，僅有沖蝕作用；由高程對比可知，針對自由液面與床沙表面整體形態模擬，本文模型均方根誤差均優于Fourtakas 模型，網格模型模擬誤差最大.

在t=0.75 s 時刻［圖10（c）］和t=1.00 s 時刻［圖10（d）］，本文模型自由液面均方根誤差略大于Fourtakas 模型，而水沙交界面顯著優于Fourtakas 模型，網格模型誤差最大.本文模型所得最大沖刷深度略小于試驗值.總體上看，相較于SPH 模型，網格模型難以精準捕捉水沙交界面處大變形、強非線性行為，水面及沖刷地形模擬誤差均較大.

為進一步對SPH 數值仿真結果的局部誤差進行分析，選取壩內斷面（S1）、最大沖深位置斷面（S2）、Louvain 試驗水頭處斷面（S3）共3處斷面對本文模型與Fourtakas 模型模擬誤差進行比較，各斷面分布匯總如圖11 所示.分別提取各時刻兩種模型所有斷面的自由液面高程及泥沙界面高程，并計算局部誤差率繪制成圖12.由圖12（a）可知，兩種模型在水頭處（S3）液面模擬誤差相比于其他斷面均偏大，這與水頭推移位置模擬略遠于試驗結果相關.Fourtakas 模型液面高程整體上略低于本文模型，此種現象是由其河床沖蝕深度整體大于試驗結果造成的，河床沖蝕越深使得水位高程相應地沉降越多.本文模型河床沖蝕趨勢與試驗結果更為貼近，但液面整體上略高于試驗結果，推測原因可能為：Louvain 潰壩侵蝕試驗采用PVC 塑料顆粒模擬泥沙，而顆粒間存在間隙，隨著試驗的進行，水流會發生滲流，導致水位下降；Fourtakas 模型與本文模型均采用SPH 方法模擬潰壩沖刷，未能考慮水流向下滲透作用，最終造成本文模型沖刷結果相近而水位偏高的情況，此亦能較好地解釋Fourtakas模型沖蝕大于試驗結果而水位幾乎與試驗持平的現象.

圖11 對比斷面分布示意圖Fig.11 Configuration of the compared sections

圖12 誤差率變化圖Fig.12 Development of percentage of error

圖12（b）中S1、S3 斷面處泥沙高程模擬，本文模型明顯優于Fourtakas 模型，Fourtakas 模型在最大沖深位置前后的模擬誤差較大，本文模型改進了這一不足；而各個時刻最大沖深結果，兩模型各有優劣.圖13 給出了Louvain 試驗、Fourtakas 模型以及本文模型最大沖深處高程隨時間變化趨勢圖.由圖13 可知，在0.50 s 前，本文模型可以較為精準地預測最大沖刷深度.本文模型和Fourtakas 模型模擬的最大沖刷深度均出現在0.50 s 時刻，而實際最大沖深出現在 0.75 s 時刻；自0.50 s 后，本文模型和Fourtakas 模型最大沖刷深度均出現回升，推測原因為沖刷坑周圍泥沙隨水流回填至最大沖深處，造成高程回升，而此現象亦在原試驗0.75 s 之后發生.整體上看，本文模型所得全過程最大沖深略小于試驗值，而Fourtakas 模型略大于試驗值.另外，泥沙回填與其流動性相關，最終由其黏度決定，本文泥沙采用HBP 模型，其黏度受m、n值控制，為統一變量從而進行沖刷精度驗證，本文參數取值參照了Fourtakas 模型；然而，在本文模型中由于引入了Shields準則作為屈服泥沙強度計算依據，原有模型參數取值可能造成黏度計算誤差，從而加劇了沖刷坑泥沙回填現象.事實上，Louvain 潰壩侵蝕試驗對照研究中［9，11，13，28］，即使采用相同泥沙本構，也未出現模型參數取值統一，基于SPH 法潰壩沖刷最大深度的精確模擬，仍需在后續研究工作中對泥沙模型及參其數取值深入研究.

圖13 最大沖深位置高程對比Fig.13 Comparison of the elevation of maximum scour position

綜上，相比于Fourtakas模型，本文建立的改進模型能夠對試驗河床整體沖刷形態、演變及推進過程進行更為精確地模擬；兩種模型最大沖刷深度與真實結果均存在一定誤差.總體上看，使用本文提出的兩相流SPH改進算法，仿真計算結果較為可信.

4 結論與展望

1）提出一種基于光滑粒子元的水-沙兩相流沖刷數值仿真改進方法.將泥沙視為一種非牛頓流體，引入HBP 模型描述其非牛頓流體行為.根據泥沙運動狀態將其劃分為沉積物、推移質、懸移質三種類型，引入上升與沉降流速模型作為三種不同狀態泥沙轉化判斷依據.應用Drucker-Prager 與Shields 侵蝕準則計算泥沙起動臨界切應力，并為不同狀態泥沙粒子賦予不同流變特性.

2）考慮水流升力與水流拖曳力，推導了在三維斜坡床面上泥沙起動臨界切應力修正系數，作為水平床面泥沙侵蝕準則的補充.

3）基于所提出的SPH 多相流模型河床沖刷數值仿真改進算法，采用三維潰壩沖刷算例與試驗結果以及其他同類仿真結果對比驗證模型精度.結果表明，相比于其他同類模型，本文建立的改進模型能高精度模擬試驗河床整體沖刷形態、演變及推進過程.

4）當前模型最大沖深與試驗結果仍存在一定誤差，模型中亦未能考慮水體滲透影響，后續工作將針對上述問題開展進一步研究.