基于Flow－3D 的東安工程樁壩沖刷數(shù)值模擬研究

栗銘陽(yáng)，張寶森，方國(guó)華，汪自力

（1.水利部堤防安全與病害防治工程技術(shù)研究中心，河南 鄭州 450003； 2.黃河水利委員會(huì)黃河水利科學(xué)研究院，河南 鄭州 450003； 3.河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098）

透水樁壩施工快速、方便，物料好選、經(jīng)濟(jì)實(shí)用，且防沖作用顯著［1］，20 世紀(jì)70 年代已開(kāi)始用于河道整治工程，在河南省內(nèi)目前已應(yīng)用在韋灘工程、孤柏嘴工程、張王莊工程、東安工程等黃河河道整治工程中［2］。透水樁壩主要利用樁壩的滯流減速作用達(dá)到控導(dǎo)河勢(shì)、固灘保堤的目的，設(shè)計(jì)時(shí)可作為不搶險(xiǎn)建筑物或少搶險(xiǎn)建筑物。 但在實(shí)際工程運(yùn)行中，在復(fù)雜水流條件下，透水樁壩會(huì)發(fā)生倒塌、斷裂，危及黃河堤防的安全［3］。 因此，研究透水樁壩在不同水流條件下的沖刷情況很有必要。

透水樁壩的研究主要圍繞類似樁式水工建筑物的繞流特性、樁基礎(chǔ)附近沖刷坑的演變過(guò)程和沖刷坑的形態(tài)、深度進(jìn)行。 在透水樁結(jié)構(gòu)的影響下，樁結(jié)構(gòu)附近的水流形態(tài)發(fā)生改變，水流流過(guò)樁壩后產(chǎn)生回流。 此外，樁結(jié)構(gòu)的阻水作用導(dǎo)致產(chǎn)生下潛水流，下潛水流與繞壩水流互相產(chǎn)生影響，形成馬蹄形旋渦［4］。 下潛水流的垂直射流能力較強(qiáng)，對(duì)河床床沙產(chǎn)生直接沖擊，馬蹄形旋渦會(huì)產(chǎn)生較大剪切應(yīng)力，提升水流輸沙能力，在垂直射流和馬蹄形旋渦的作用下，樁基礎(chǔ)附近發(fā)生局部沖刷。

透水樁壩局部沖刷強(qiáng)度的影響因素可以歸結(jié)為結(jié)構(gòu)參數(shù)、水力參數(shù)和河床泥沙參數(shù)三類。 高先剛等［5］初步探討了透水樁壩前沖刷深度的影響因素，認(rèn)為單寬流量和入流角度對(duì)沖刷深度影響最大。 Masoud Gh?odsian 等［6］通過(guò)物理試驗(yàn)研究了90°彎道丁壩長(zhǎng)度與沖刷深度之間的關(guān)系，表明沖刷深度與丁壩長(zhǎng)度正相關(guān)。 劉燕等［7］基于東安工程、韋灘工程進(jìn)行長(zhǎng)系列試驗(yàn)，對(duì)透水樁壩的導(dǎo)流效果進(jìn)行分析。 考慮到物理試驗(yàn)的復(fù)雜性，采用數(shù)值模擬方法研究樁基礎(chǔ)的泥沙沖刷已經(jīng)成為近些年的熱點(diǎn)。 鈄錦周等［8］利用Flow－3D對(duì)海上風(fēng)電單樁進(jìn)行數(shù)值模擬，分析了流速和樁徑對(duì)局部沖刷結(jié)果的影響。 樊新建等［9］進(jìn)行了水力插板透視丁壩結(jié)構(gòu)數(shù)值模擬研究，得到了不同體形參數(shù)下水力插板透水丁壩水流結(jié)構(gòu)的變化規(guī)律。 黃佳麗等［10］針對(duì)Flow－3D 泥沙模型，確定了影響射流沖刷過(guò)程的關(guān)鍵參數(shù)。 戚藍(lán)等［11］采用RNGk－ε模型，用VOF 方法進(jìn)行自由表面的追蹤，并對(duì)丁壩局部沖刷坑形態(tài)及沖刷深度進(jìn)行了數(shù)值模擬。

透水樁壩的倒塌主要原因在于樁壩附近床面的泥沙受到水流的長(zhǎng)期作用，導(dǎo)致河床泥沙發(fā)生了液化、流動(dòng)，并最終影響了樁壩的穩(wěn)定。 然而，該方面的研究比較缺乏。 考慮到物理試驗(yàn)的復(fù)雜性，筆者利用Flow－3D 軟件對(duì)東安透水樁壩工程（簡(jiǎn)稱東安工程）進(jìn)行了不同入流角度下的樁基礎(chǔ)沖刷模擬，利用現(xiàn)場(chǎng)測(cè)量得到的水深、流速等確定模擬參數(shù)，建立東安透水樁壩沖刷模型，對(duì)其沖刷機(jī)理進(jìn)行研究，以期為透水樁壩的設(shè)計(jì)和后期安全運(yùn)行提供依據(jù)。

1 工程概況

東安工程位于河南省焦作市武陟縣，地處黃河、沁河交匯處。 東安工程所處河段從20 世紀(jì)80 年代開(kāi)始有計(jì)劃地進(jìn)行了整治。 按照控導(dǎo)主流、護(hù)灘保堤的方針，采用因勢(shì)利導(dǎo)、以壩護(hù)灣、以灣導(dǎo)流的方法整治中水河槽。 在棗樹(shù)溝控導(dǎo)工程、東安透水樁壩工程束縛下，該河段已成為人工控制的彎曲型河段［12］。

東安工程采用單排鋼筋混凝土灌注樁結(jié)構(gòu)，灌注樁直徑0.8 m、樁長(zhǎng)29 m，相鄰樁中心間距1.2 m。 樁壩上部設(shè)有工作便橋，工作便橋上布置欄桿、扶手，欄桿采用立方體樁式結(jié)構(gòu)，扶手采用鏈接式，建造時(shí)對(duì)兩兩相鄰的樁壩進(jìn)行連續(xù)梁澆筑，形成連續(xù)的樁排。 工程始建于2000 年，運(yùn)行使用過(guò)程中在原基礎(chǔ)上延伸，直到2020 年整體完工，最后一次擴(kuò)建由原始長(zhǎng)度5 748 m擴(kuò)至6 714 m，長(zhǎng)度增加了966 m。 東安樁壩于2021 年4 月11 日前后出險(xiǎn)，出險(xiǎn)時(shí)根據(jù)花園口水文站記錄，流量為995 m3／s，在小流量作用下透水樁壩連續(xù)兩次發(fā)生樁排倒塌。 據(jù)現(xiàn)場(chǎng)拍攝視頻顯示，出險(xiǎn)時(shí)透水樁壩本身鋼筋混凝土結(jié)構(gòu)沒(méi)有發(fā)生破壞，向河岸側(cè)整體倒塌。 此外，在多年運(yùn)行下東安工程樁壩混凝土結(jié)構(gòu)磨損嚴(yán)重，相鄰樁間泥沙淤積，堵塞了樁間的過(guò)流通道，降低了樁壩的緩流落淤效果。 樁壩的中軸線也在水流的長(zhǎng)期作用下發(fā)生偏移［13］。

2 模擬計(jì)算方程

Flow－3D 模擬計(jì)算的連續(xù)性方程［14］為

動(dòng)量方程為

式中：u、v、w分別為x、y、z方向的瞬時(shí)速度分量；t為時(shí)間；VF為可流動(dòng)的體積分?jǐn)?shù)；Ax、Ay、Az為x、y、z方向的可流動(dòng)面積分?jǐn)?shù)；ρ為流體密度；p為瞬時(shí)壓強(qiáng)；Gx、Gy、Gz為x、y、z方向的重力加速度；fx、fy、fz為x、y、z方向黏性加速度。

湍流模型采用RNGk－ε，具體計(jì)算控制方程見(jiàn)參考文獻(xiàn)［14］，這里重點(diǎn)介紹Flow－3D 模擬計(jì)算所采用的泥沙控制方程。

無(wú)量綱臨界希爾茲數(shù)θcr由基于Soulsby 公式［6］計(jì)算得到的曲線確定，當(dāng)局部希爾茲數(shù)θi大于臨界希爾茲數(shù)θcr時(shí)，河床表面泥沙起動(dòng)，臨界希爾茲數(shù)計(jì)算公式為

其中

式中：ρs為泥沙密度；d50為泥沙中值粒徑，本文為0.385 mm；μ為動(dòng)力黏滯系數(shù)；g為重力加速度向量。

斜坡上泥沙顆粒受到重力作用，在順斜坡方向存在重力的分力更容易讓泥沙顆粒在水流作用下產(chǎn)生移動(dòng)，為了減弱這種影響，將式（1）修正為

式中：θ′cr為修正臨界希爾茲數(shù)；ˉω為河床坡度；ψ為水流方向與坡面夾角；φ為泥沙休止角。

水流挾帶床沙的上升速度ulift基于Mastbergen 和Berg 理論［15］進(jìn)行計(jì)算，公式為

式中：ai為泥沙挾帶系數(shù)；ns為法向向量；τ為局部切應(yīng)力。

河床泥沙顆粒的沉降與起動(dòng)在理論上是兩種相反的過(guò)程，但實(shí)際上這兩種過(guò)程往往同時(shí)發(fā)生。 本文采用Soulsby 提出的沉降速度公式計(jì)算床沙的沉降速度：

推移質(zhì)泥沙顆粒運(yùn)動(dòng)往往是顆粒本身在河床表面發(fā)生滾動(dòng)、跳躍等，采用Van Rijn 輸沙率公式［16］計(jì)算河床表面上泥沙的單位寬度輸沙率：

其中

式中：φ為泥沙休止角；β為推移質(zhì)系數(shù)，取0.053；d為泥沙顆粒直徑。

推移質(zhì)層的厚度δ采用Van Rijn 公式計(jì)算：

3 數(shù)值模型布置及驗(yàn)證

3.1 模型布置

根據(jù)現(xiàn)場(chǎng)情況，以“口門”處倒塌的樁壩作為原型，根據(jù)工程實(shí)際尺寸，選取倒塌位置的16 根樁壩，以現(xiàn)場(chǎng)采用多波束測(cè)深系統(tǒng)和ADCP 測(cè)量所得水深等作為參數(shù)，分別進(jìn)行45°、60°、90°入流角度的沖刷模擬。將計(jì)算域高度設(shè)置為29 m（其中河床厚度10 m），寬度設(shè)置為30 m，為了使水流平順，進(jìn)出口距離為60 m。模型布置及模擬區(qū)域位置見(jiàn)圖1。

圖1 模型布置及模擬區(qū)域位置

計(jì)算域采用均質(zhì)泥沙，粒徑選用黃河泥沙平均粒徑0.1 mm，密度為2 650 kg／m3，臨界希爾茲數(shù)為0.025，泥沙休止角為32°，表面粗糙度選取2.5d50， 泥沙挾帶系數(shù)為0.018。 模型中進(jìn)出口邊界均為壓力水深邊界，進(jìn)口水深設(shè)置為出險(xiǎn)時(shí)現(xiàn)場(chǎng)測(cè)量的平均水深7 m，出口水深設(shè)置為5 m，底部采用無(wú)滑移邊界條件，其余為對(duì)稱邊界。 計(jì)算域采用結(jié)構(gòu)化正交網(wǎng)格，網(wǎng)格總數(shù)為120 萬(wàn)個(gè)。 采用被廣泛應(yīng)用的修正湍流黏度的RNGk － ε模型進(jìn)行計(jì)算，模擬時(shí)間為20 min。

3.2 模型驗(yàn)證

模型計(jì)算結(jié)果采用多波束測(cè)深儀測(cè)量得到的水下地形數(shù)據(jù)進(jìn)行驗(yàn)證，圖2 為利用sufer 軟件處理得到的東安工程出險(xiǎn)位置水下地形等高線分布，在“口門”附近切割4 個(gè)斷面，從水下地形斷面的分布情況定量得到?jīng)_刷坑的深度，與模擬結(jié)果進(jìn)行對(duì)比，從而驗(yàn)證模擬結(jié)果。

圖2 出險(xiǎn)位置水下地形等高線分布（m）

從提取的4 個(gè)斷面地形圖可以大致判斷沖刷坑的深度為2 m 左右，寬度約40 m，沖刷坑分布在河道中間位置，沖刷坑的形態(tài)中間部位有輕微的凸起，各個(gè)斷面最大沖刷深度見(jiàn)表1，實(shí)測(cè)斷面地形見(jiàn)圖3。 模擬區(qū)域不同入流角度的沖刷深度與現(xiàn)場(chǎng)實(shí)測(cè)的沖刷深度相差不大，60°、90°入流角度下的沖刷深度與實(shí)測(cè)結(jié)果擬合較好，最大誤差僅為9.4%。

圖3 實(shí)測(cè)斷面地形

表1 實(shí)測(cè)與模擬最大沖刷深度對(duì)比

此外，從45°、60°、90°入流角度下沖刷坑模擬結(jié)果與實(shí)測(cè)三維地形圖比較（見(jiàn)圖4）來(lái)看，沖刷坑的形態(tài)與實(shí)測(cè)形態(tài)基本一致，即樁基附近形成類似于V 形的沖刷坑，沖刷坑后產(chǎn)生淤積。

圖4 模擬結(jié)果與實(shí)測(cè)結(jié)果對(duì)比

綜上，可以認(rèn)為通過(guò)數(shù)值模擬可研究區(qū)域的沖刷情況。

4 數(shù)值模擬結(jié)果

4.1 不同入流角度下沖刷湍流強(qiáng)度分析

不同入流角度下樁壩附近湍流強(qiáng)度云圖及沖刷深度等值線見(jiàn)圖5～圖7（圖中負(fù)號(hào)代表沖刷）。 隨著入流角度的增大，樁壩后的湍流強(qiáng)度峰值逐漸增大，且湍流強(qiáng)度的峰值區(qū)逐漸向下移動(dòng)。 45°入流角度下樁壩后湍流強(qiáng)度普遍不大，60°入流角度下壩后湍流強(qiáng)度峰值區(qū)分布在上半部分，90°入流角度下樁壩后湍流強(qiáng)度峰值區(qū)分布在中下部。 結(jié)合沖刷深度的模擬結(jié)果發(fā)現(xiàn)，樁壩附近的沖刷深度與水體的湍流有關(guān)，水流流過(guò)樁壩產(chǎn)生湍流，湍流促進(jìn)了沖刷坑的發(fā)展，湍流強(qiáng)度大的地方?jīng)_刷坑的深度大，沖刷坑的形態(tài)與湍流強(qiáng)度云圖的峰值區(qū)形態(tài)類似［17］，45°、60°、90°入流角度下樁基礎(chǔ)周圍最大湍流強(qiáng)度分別為39.2%、40.0%、56.7%。

圖5 45°入流條件下湍流強(qiáng)度云圖與沖刷深度等值線

圖6 60°入流條件下湍流強(qiáng)度云圖與沖刷深度等值線

圖7 90°入流條件下湍流強(qiáng)度云圖與沖刷深度等值線

4.2 不同入流角度下的剪切應(yīng)力與剪切速率分析

入流角度為45°、90°時(shí)樁壩附近最大沖刷深度點(diǎn)的剪切應(yīng)力、剪切速率歷時(shí)曲線見(jiàn)圖8。 可以看出，剪切應(yīng)力曲線在初始階段有一個(gè)陡升的過(guò)程，而后趨于穩(wěn)定，入流角度為90°時(shí)剪切應(yīng)力曲線陡升后還存在陡降的過(guò)程。 考慮到剪切應(yīng)力與流速正相關(guān)，初始時(shí)刻流速由零增大到某一數(shù)值，因此剪切應(yīng)力在初始時(shí)刻較大；之后受樁壩的緩流作用，流速逐漸減小，對(duì)應(yīng)的剪切應(yīng)力減小，隨著流場(chǎng)的穩(wěn)定，剪切應(yīng)力穩(wěn)定在某一范圍內(nèi)。 剪切速率是表征水流中各層之間流速變化的量［18］，初始階段流場(chǎng)未穩(wěn)定時(shí)，各層水流間存在流速差，剪切速率增大；流場(chǎng)穩(wěn)定時(shí)，相鄰層水體間的流速差較小，因而比較穩(wěn)定。

圖8 不同入流角度下最大沖刷深度點(diǎn)的剪切應(yīng)力、剪切速率歷時(shí)曲線

入流角度為90°對(duì)應(yīng)的剪切應(yīng)力顯著大于入流角度為45°時(shí)的，入流角度為90°對(duì)應(yīng)的剪切應(yīng)力穩(wěn)定在10.4 N／m2，而入流角度為45°時(shí)穩(wěn)定在2.8 N／m2；入流角度為90°對(duì)應(yīng)的剪切速率穩(wěn)定在0.10 s－1附近，而入流角度為45°時(shí)為0.04 s－1左右。

4.3 不同入流角度下沖刷深度分析

不同入流角度下樁壩基礎(chǔ)附近沖刷深度歷時(shí)曲線見(jiàn)圖9。 可以看出，隨著入流角度的增大，樁壩周圍最大沖刷深度也相應(yīng)增大。 沖刷歷時(shí)1 200 s 時(shí)，入流角度為45°、60°、90°對(duì)應(yīng)的沖刷深度分別為1.6、2.3、2.4 m。 入流角度越大，沖刷過(guò)程的回淤現(xiàn)象越不明顯，原因是入流角度增大時(shí)，湍流強(qiáng)度增大，沖刷作用增強(qiáng)。 45°、90°入流角度下的沖刷深度歷時(shí)曲線中出現(xiàn)水平段，原因是沖刷到某一限度時(shí)，進(jìn)一步?jīng)_刷需要一定的緩沖時(shí)間。 入流角度為45°時(shí)可能因沖刷強(qiáng)度不足而導(dǎo)致回淤，入流角度為90°時(shí)湍流強(qiáng)度過(guò)大導(dǎo)致沖刷過(guò)程存在停滯階段。 此外，入流角度為60°時(shí)最大沖刷深度略小于入流角度為90°時(shí)的，但沖刷過(guò)程有所區(qū)別，由此可見(jiàn)，工程實(shí)際中不能僅關(guān)注主流頂沖這一工況，對(duì)入流角度較大的工況也應(yīng)關(guān)注。

圖9 沖刷深度歷時(shí)曲線

此外，對(duì)比現(xiàn)場(chǎng)測(cè)量情況來(lái)看，60°、90°入流角度下的模擬結(jié)果與實(shí)際沖刷深度較為接近，可以初步認(rèn)為本次出險(xiǎn)是入流角度過(guò)大導(dǎo)致樁排倒塌。

5 結(jié)論

以黃河下游東安混凝土透水樁壩工程出險(xiǎn)為背景，基于工程現(xiàn)場(chǎng)多波束測(cè)深系統(tǒng)、ADCP 等設(shè)備測(cè)量所得數(shù)據(jù)，通過(guò)Flow－3D 軟件對(duì)不同入流角度條件的樁壩沖刷進(jìn)行了數(shù)值模擬，得到以下主要結(jié)論：

（1）透水樁壩受上游來(lái)水螺旋環(huán)流淘刷，壩前形成近乎平行于樁壩的沖刷坑，穿過(guò)透水樁壩的含沙水流會(huì)對(duì)壩后產(chǎn)生沖刷。 模擬結(jié)果顯示，湍流的高強(qiáng)度區(qū)域隨著入流角度的增大逐漸向下移動(dòng)。 沖刷坑形態(tài)與湍流高強(qiáng)度區(qū)的分布有關(guān)，沖刷較深的區(qū)域往往是湍流高強(qiáng)度區(qū)分布的地方，模擬結(jié)果中45°、60°、90°入流角度下樁基礎(chǔ)附近最大湍流強(qiáng)度分別為39.2%、40.0%、56.7%。

（2）90°入流角度對(duì)應(yīng)的剪切應(yīng)力、剪切速率穩(wěn)定在10.4 N／m2、0.10 s－1左右，45°入流角度對(duì)應(yīng)的剪切應(yīng)力、剪切速率穩(wěn)定在2.8 N／m2、0.04 s－1左右。 90°入流角度對(duì)應(yīng)的剪切應(yīng)力約為45°入流角度時(shí)的4 倍，剪切速率約為2 倍。 模擬結(jié)果中45°入流角度對(duì)應(yīng)的最大沖刷深度為1.6 m，而90°入流角度對(duì)應(yīng)的最大沖刷深度為2.4 m，對(duì)比45°、90°入流條件下的剪切應(yīng)力和沖刷深度，沖刷深度隨著剪切應(yīng)力的增大而增大，相較剪切應(yīng)力的變化，沖刷深度的增幅要小得多。

（3）入流角度越大，沖刷深度越大，沖刷過(guò)程中回淤現(xiàn)象越不明顯。 但入流角度小時(shí)，回淤導(dǎo)致沖刷過(guò)程不平穩(wěn)；入流角度過(guò)大時(shí)，湍流強(qiáng)度過(guò)大導(dǎo)致沖刷過(guò)程不平穩(wěn)；60°入流角度對(duì)應(yīng)的最大沖刷深度略小于90°入流角度的。 從工程的安全性考慮，不僅需考慮主流頂沖時(shí)樁壩的穩(wěn)定性問(wèn)題，也應(yīng)對(duì)小角度、小流量的沖刷穩(wěn)定性進(jìn)行關(guān)注。 對(duì)比現(xiàn)場(chǎng)測(cè)量情況，出險(xiǎn)位置最大沖刷深度為2 m 左右，與60°、90°入流角度下模擬結(jié)果接近，初步認(rèn)為入流角度過(guò)大是造成本次樁排倒塌的原因。