張嶸釗 ，熊文 ?，劉川渟

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

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

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

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

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

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

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

1 SPH基本理論

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

1.1 SPH法的積分插值理論

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

核函數(shù)W是關(guān)于光滑長(zhǎng)度h的函數(shù)，其作用是將連續(xù)微分方程基于插值函數(shù)在特定點(diǎn)處進(jìn)行積分，任意物理量連續(xù)形式用F可表示為

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

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

式中：q為粒子間距與平滑長(zhǎng)度的比值.

1.2 控制方程的離散形式

N-S 方程中質(zhì)量守恒方程與動(dòng)量守恒方程分別表示為式（4）與式（5）：

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

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

式中：下標(biāo)a表示中心粒子；下標(biāo)b表示鄰域粒子；下標(biāo)ab表示粒子a與b的差；?a為針對(duì)粒子a的哈密頓算子；m為粒子質(zhì)量；ρa(bǔ)與ρb分別表示粒子a與粒子b的密度.

通過(guò)Favre 平均法對(duì)密度進(jìn)行加權(quán)平均，可將SPS 項(xiàng)引入SPH 方法中，動(dòng)量守恒方程改寫為SPH離散形式［16］：

式中：μ0為運(yùn)動(dòng)黏度；τij為SPS應(yīng)力矢量.

2 基于SPH的兩相流數(shù)值模型

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

2.1 泥沙相非牛頓流體模型

HBP模型計(jì)算流體表觀黏度如下：

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

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

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

2.2 泥沙侵蝕模型

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

2.2.1 泥沙屈服強(qiáng)度模型

Drucker-Prager（DP）屈服準(zhǔn)則一般形式為：

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

式中：φ為內(nèi)摩擦角；c為土體黏聚力.

根據(jù)Fourtakas 等［13］的研究，J2可由土粒子所受水體剪切應(yīng)力ταβ計(jì)算：

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

聯(lián)立式（11）、式（12）可以建立J2與ⅡD間聯(lián)系：

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

當(dāng)沉積物粒子所受切應(yīng)力大小不超過(guò)其屈服強(qiáng)度|τy|時(shí)，土體粒子由于高黏度特性幾乎不發(fā)生剪切變形，其表觀黏度μapp由|τy|代入式（8）進(jìn)行計(jì)算.

對(duì)于河床表面與水接觸的土粒子，若滿足由沉積物轉(zhuǎn)化為推移質(zhì)的條件，土體屈服強(qiáng)度采用Shields準(zhǔn)則計(jì)算并替換，詳見(jiàn)2.2.2節(jié).

2.2.2 泥沙起動(dòng)判定

初始條件下，將泥沙層置于水層下方，當(dāng)泥沙粒子所受切應(yīng)力不超過(guò)臨界值時(shí)，泥沙粒子固定不動(dòng)；達(dá)到臨界值時(shí)，泥沙由沉積物轉(zhuǎn)化為推移質(zhì)，隨水流沿床面輸運(yùn).

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

轉(zhuǎn)化為推移質(zhì)的泥沙粒子須位于河床表面，同時(shí)考慮到推移質(zhì)粒子為飽和土的特性，其轉(zhuǎn)化條件設(shè)置以下兩個(gè)：

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

2）待轉(zhuǎn)化泥沙粒子實(shí)際質(zhì)量應(yīng)小于閾值質(zhì)量.

當(dāng)泥沙粒子同時(shí)滿足以上條件時(shí)，基于Shields侵蝕準(zhǔn)則，計(jì)算水平床面臨界轉(zhuǎn)化剪切應(yīng)力：

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

θcr求解基于雷諾數(shù)Re［12］：

式中：R*=，v為水體運(yùn)動(dòng)黏度，u*為摩阻流速，按式（17）計(jì)算：

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

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

式中：ks為當(dāng)量砂徑糙率，按經(jīng)驗(yàn)取泥沙中值粒徑.

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

水流沖刷的起動(dòng)力τb由Einstein 對(duì)數(shù)流速分布公式［20］計(jì)算：

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

若泥沙粒子所受剪切應(yīng)力τb超過(guò)臨界剪切應(yīng)力τbcr，則判斷泥沙粒子由沉積物轉(zhuǎn)變?yōu)橥埔瀑|(zhì).轉(zhuǎn)化后粒子表觀黏度需結(jié)合HBP模型進(jìn)行更新，將τbcr代入式（8）中代替屈服應(yīng)力力τc，計(jì)算得到黏度μori.

2.2.3 泥沙上升與沉降判定

已屈服推移質(zhì)粒子在滿足條件時(shí)會(huì)轉(zhuǎn)化為懸移質(zhì)，被水流挾帶而發(fā)生輸運(yùn).當(dāng)流速減小、大量懸移質(zhì)粒子聚集時(shí)，懸移質(zhì)再次沉降，轉(zhuǎn)化為推移質(zhì).

懸移質(zhì)周圍泥沙粒子體積濃度Cv，i不超過(guò)臨界濃度Cvcr：

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

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

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

若u≤uset，則判斷懸疑質(zhì)沉降，轉(zhuǎn)化為推移質(zhì)，黏度按推移質(zhì)計(jì)算.

2.2.4 臨界起動(dòng)切應(yīng)力斜坡修正

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

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

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

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

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

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

抵抗起動(dòng)力FC由庫(kù)倫摩擦力產(chǎn)生，與F方向相反，在極限狀態(tài)下，

式中：φ表示泥沙內(nèi)摩擦角.當(dāng)泥沙粒子達(dá)到臨界起動(dòng)狀態(tài)，即滿足

由式（28），推導(dǎo)出|FD|表達(dá)式為：

由式（29）得到斜坡床面上的臨界起動(dòng)切應(yīng)力大小.對(duì)于水平床面，通過(guò)受力分析易得：

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

式中：k為臨界起動(dòng)切應(yīng)力斜坡修正系數(shù)：

水平面的臨界起動(dòng)切應(yīng)力經(jīng)過(guò)斜坡修正為

2.3 泥沙模塊運(yùn)行流程

本文沖刷模型中，初始狀態(tài)下泥沙粒子皆為沉積物，經(jīng)由每個(gè)時(shí)間步進(jìn)行處理，將依次完成沉積物向推移質(zhì)轉(zhuǎn)化、推移質(zhì)向懸移質(zhì)轉(zhuǎn)化以及懸移質(zhì)向推移質(zhì)轉(zhuǎn)化.轉(zhuǎn)化完成后，分別存儲(chǔ)3 種狀態(tài)的泥沙粒子，進(jìn)入SPH 求解步驟，直至計(jì)算結(jié)束，計(jì)算流程示意圖如圖3所示.

2.3.1 預(yù)處理模塊

獲取上一時(shí)間步中泥沙粒子狀態(tài)信息；迭代計(jì)算摩阻流速并應(yīng)用Shields侵蝕準(zhǔn)則得到臨界Shields數(shù)；按Mastbergen公式計(jì)算上升與沉降流速.

2.3.2 推移質(zhì)模塊

沉積物粒子轉(zhuǎn)化為推移質(zhì)粒子須滿足2.2.2 節(jié)所述起動(dòng)條件.因此，需判斷其是否處于沙床表面，進(jìn)而利用Einstein 對(duì)數(shù)流速分布公式計(jì)算該粒子所受流體切應(yīng)力.若切應(yīng)力超過(guò)屈服應(yīng)力，則儲(chǔ)存為新的推移質(zhì)粒子，并結(jié)合HBP 模型利用Shield 準(zhǔn)則計(jì)算屈服強(qiáng)度，獲得推移質(zhì)粒子表觀黏度.

2.3.3 懸移質(zhì)模塊

推移質(zhì)粒子轉(zhuǎn)化為懸移質(zhì)粒子，須滿足2.2.3 節(jié)所述條件.由于水流裹挾懸移質(zhì)輸送存在濃度閾值，首先判斷粒子鄰域內(nèi)泥沙體積濃度，當(dāng)濃度小于閾值時(shí)，判斷其速度大小，若速度大于上升流速，則判斷為懸移質(zhì)粒子并依據(jù)Vand模型更新粒子黏度.

對(duì)于懸移質(zhì)粒子，其速度若小于沉降流速，則又轉(zhuǎn)化為推移質(zhì)粒子，返回推移質(zhì)模塊.

2.3.4 沉積物模塊

泥沙粒子不起動(dòng)，基于DP準(zhǔn)則結(jié)合HBP模型更新其屈服強(qiáng)度.

3 潰壩沖刷數(shù)值仿真

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

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

3.1 基于DualSPHysics的沖刷仿真計(jì)算流程

本文在DualSPHysics 開(kāi)源代碼的基礎(chǔ)上，實(shí)現(xiàn)沉積物粒子起動(dòng)判斷與類型轉(zhuǎn)化.二次開(kāi)發(fā)模塊分為三部分執(zhí)行：①初始化模塊；②判斷模塊；③黏度計(jì)算模塊.程序的運(yùn)行流程如圖4所示.

圖4 程序運(yùn)行流程Fig.4 Flow of program executing

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

3.2 數(shù)值模型參數(shù)設(shè)置

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

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

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

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

圖6 流體粒子間距結(jié)果對(duì)比Fig.6 Comparison of the distance between particles

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

3.3 潰壩沖刷形態(tài)及誤差分析

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

圖7 網(wǎng)格模型示意圖Fig.7 Diagram of mesh-model

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

圖8 數(shù)值模型示意圖Fig.8 Diagram of the numerical model

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

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

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

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

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

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

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

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

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

綜上，相比于Fourtakas模型，本文建立的改進(jìn)模型能夠?qū)υ囼?yàn)河床整體沖刷形態(tài)、演變及推進(jìn)過(guò)程進(jìn)行更為精確地模擬；兩種模型最大沖刷深度與真實(shí)結(jié)果均存在一定誤差.總體上看，使用本文提出的兩相流SPH改進(jìn)算法，仿真計(jì)算結(jié)果較為可信.

4 結(jié)論與展望

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

2）考慮水流升力與水流拖曳力，推導(dǎo)了在三維斜坡床面上泥沙起動(dòng)臨界切應(yīng)力修正系數(shù)，作為水平床面泥沙侵蝕準(zhǔn)則的補(bǔ)充.

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

4）當(dāng)前模型最大沖深與試驗(yàn)結(jié)果仍存在一定誤差，模型中亦未能考慮水體滲透影響，后續(xù)工作將針對(duì)上述問(wèn)題開(kāi)展進(jìn)一步研究.