基于有限體積法的三維局部沖刷數(shù)值模型研究及應(yīng)用

李登松，戴光清，金 鑫，曾 赟，漆力健

(1.四川農(nóng)業(yè)大學 水利水電學院，雅安 625014；2.四川大學 水力學與山區(qū)河流開發(fā)保護國家重點實驗室，成都 610065；3.成都理工大學 能源學院，成都610059)

1 引 言

水流在流經(jīng)橋梁墩臺、丁壩和河海岸線等水工建筑物時，極易產(chǎn)生局部沖刷，導致水工建筑物的毀壞。因此，局部沖刷的破壞機制和程度一直備受關(guān)注，目前主要借助物理模型試驗進行研究。大量的試驗結(jié)果表明，沖刷破壞的主要原因是漩渦和下潛水流[1,2]，其發(fā)展大致分為形成初始沖坑、沖坑發(fā)展和平衡沖深三階段[3]。破壞程度主要根據(jù)經(jīng)驗或半經(jīng)驗性沖深公式的計算結(jié)果進行預測，如Garde等[4]的丁壩沖深公式和Melville等[5，6]的橋臺公式。但是試驗存在比尺效應(yīng)，而且難以準確全面地捕捉?jīng)_刷坑發(fā)展過程和坑內(nèi)流場信息等。

與物理模型試驗相比，數(shù)值計算彌補了上述缺點，并且可以預測沖坑地形的演變。目前，數(shù)值模型的開發(fā)主要基于有限體積法(FVM)、有限差分法(FDM)和光滑粒子法(SPH)等。SPH法根據(jù)拉格朗日法求解，將其應(yīng)用在局部沖刷模擬的過程中則為追蹤泥沙顆粒的運動軌跡[7]。FVM和FDM法均是根據(jù)歐拉法求解N-S方程，不同之處在于方程的求解形式[8]。在以歐拉法為基礎(chǔ)建立局部沖刷仿真平臺時，核心問題是水沙(流固)模型的耦合計算，以及水沙運動交界面的耦合處理。水沙模型的耦合包括單向和雙向耦合，大多采用單向耦合，如文獻[9，10]在模擬沖刷演變的過程中數(shù)據(jù)傳遞按照流場(N-S方程)→地形更新(泥沙連續(xù)性方程)的方向進行。水沙運動交界面的耦合處理包括動邊界法(MBM)和浸入邊界法(IBM)。動邊界法直接處理運動邊界和水沙耦合面[11]，文獻[12，13]通過該方法實現(xiàn)了圓柱橋墩和丁壩周圍局部沖刷時運動邊界的網(wǎng)格重構(gòu)。

本文以有限體積法和非結(jié)構(gòu)化的計算網(wǎng)格為基礎(chǔ)，水沙交界面數(shù)據(jù)采用單向傳遞，運動交界面應(yīng)用動邊界法，最終構(gòu)建了局部沖刷的三維數(shù)值模型。應(yīng)用該數(shù)值平臺模擬橋墩和丁壩周圍的局部沖刷，預測了沖深和沖刷過程中地形的變化，并通過相關(guān)研究者的試驗結(jié)果驗證了三維數(shù)值模型的可靠性和準確度。

2 控制方程

2.1 三維非恒定水流控制方程

局部沖刷水流為三維非恒定水流，連續(xù)性方程和納維-斯托克斯(Navier-Stokes)方程能較好地進行模擬，控制方程的積分形式如式(1,2)。

連續(xù)性方程的積分形式為

(1)

動量方程的積分形式

(2)

N-S方程(2)為瞬時運動方程，按時間平均后為雷諾方程(RANS)。雷諾方程中的切應(yīng)力項為

(3)

在SSTk-ω紊流模型中，運動渦粘系數(shù)νT、紊動能k和耗散率ω分別由式(4～6)計算。

(4)

(5)

(6)

式中F1,F2,S,σk,σω 2,Pk,β*,α,β和σk為閉合系數(shù)。

2.2 河床變形方程

在局部沖刷的過程中，水力學條件變化較為迅速，并不會立刻或較大地改變懸移質(zhì)輸沙率。而短時間內(nèi)河床地形的改變主因為推移質(zhì)通量變化。因此，忽略泥沙顆粒沿水深方向輸移擴散的通量守恒法適用于水工建筑物周圍的泥沙輸運，滿足通量守恒法的河床變形控制方程為

(7)

qB L=qbuB L|ds||δb|

(8)

式中qb為單寬體積推移質(zhì)輸沙率；uB L為平行于床面的速度矢量；ds為單元體面的邊緣長度矢量；δb為沿Z坐標軸的矢量，大小為推移質(zhì)層厚度。

在計算推移質(zhì)通量時，常用的單寬體積推移質(zhì)輸沙率公式有Meyer-Peter-Muller公式和Van Rijn公式。本文選用Van Rijn公式,即

(9)

T=(τb-τb,c)/τb,c

(10)

式中τb為床面切應(yīng)力，本文通過數(shù)值計算獲取瞬時床面切應(yīng)力；τb,c為泥沙顆粒起動時的臨界床面切應(yīng)力，該切應(yīng)力由Shields曲線獲取。

2.3 邊界條件

邊界條件分為水流進出口和壁面。水流進口邊界采用速度進口，出口邊界采用壓力出口。當水流流過固壁邊界時，主要涉及壁面上的剪切條件和粗糙度參數(shù)。剪切條件是對固壁剪切力的設(shè)定，文中移動壁面設(shè)定為無滑移條件。粗糙度常數(shù)Cs主要由顆粒類型決定，當顆粒類型與均勻沙粒相距甚遠時，結(jié)果偏差很大，在 0.5～1.0 間調(diào)整該值。但是，不存在通用的參數(shù)準則適用于任意顆粒類型。由于泥沙沖刷模擬需要泥沙顆粒粒徑，因此應(yīng)用有效粗糙率ks和糙率常數(shù)值Cs模擬均勻沙和非均勻沙屬性，取值列入表1。

表1 有效粗糙率ks取值

Tab.1 Value of effective roughnessks

壁面類型有效粗糙率ks光滑壁面0非均勻沙粒d50/dm均勻沙粒沙粒高度

此外，臨近壁面區(qū)域的網(wǎng)格處理精度直接影響附近流場求解結(jié)果的精度。因此臨近壁面的第一個網(wǎng)格幾何中心的尺寸應(yīng)大于有效粗糙率ks。

k-ω模型的近壁面處理依據(jù)以下方法,

?k/?n=0, ?w/?n=0

(11)

式中n為壁面邊界法向量。

在臨近壁面邊界的第一個網(wǎng)格中心處的值定義為

(12)

式中Cμ為調(diào)節(jié)湍動粘性幅度的常數(shù),K=0.4為卡門常數(shù),uτ為摩阻流速,yp為第一層網(wǎng)格節(jié)點與壁面的距離。

3 空間離散與求解

3.1 河床變形方程的空間離散

根據(jù)散度定理離散基于有限體積法的河床變形方程(7)后，該方程變?yōu)?/p>

(13)

式中AH為網(wǎng)格面投影在水平面上的面積,n為網(wǎng)格面上邊緣線的單位外法向量,ie為單元體面的四條邊。

由于數(shù)據(jù)存儲于非結(jié)構(gòu)化體網(wǎng)格或面網(wǎng)格的中心，在重構(gòu)網(wǎng)格面上的通量qB L時，定義中心網(wǎng)格Cell_0和鄰近網(wǎng)格Cell_1如圖1所示。

(14)

式中qb,e s為通過面es的單寬體積推移質(zhì)輸沙率,de s為面es的面積(三維)或邊緣線長度(二維),uB L為面es平行于床面的速度矢量,n01為網(wǎng)格中心和鄰近網(wǎng)格單元的法向量，方向始終由單元Cell_0指向單元Cell_1。

圖1 網(wǎng)格面通量重構(gòu)

Fig.1 Reconstructed flux of grid face

(15)

式中qb,d為中心網(wǎng)格下游的單元體，如圖2所示，采用一階迎風格式判斷單元體是上游或下游單元體；(qb)c為中心網(wǎng)格的梯度，根據(jù)高斯定理獲取，如式(17)；d為中心網(wǎng)格指向下游網(wǎng)格的矢量，大小為兩個網(wǎng)格中心的距離。

(16)

式中Ω為單元體的三維體積或二維面積,dS為單元面的面積或二維邊緣線的長度。

3.2 沙滑模型

天然河床或者沖坑絕大多數(shù)都具有一定坡度，斜坡上的泥沙顆粒存在泥沙休止角，大于泥沙休止角時，失去穩(wěn)定狀態(tài)，下滑直至達到新的穩(wěn)定狀態(tài)。因此,在更新網(wǎng)格節(jié)點位移值時，需要沙滑模型修正床面變形位移值，實現(xiàn)斜坡床面坍塌。

以圖3的網(wǎng)格單元C19作為分析對象，闡述沙滑模型。當網(wǎng)格面的坡度超過泥沙顆粒休止角后，調(diào)整節(jié)點位置，使坡度滿足泥沙休止角要求。在C19網(wǎng)格上的節(jié)點由鄰近四個網(wǎng)格單元共用，節(jié)點調(diào)整距離ΔZ受到共用網(wǎng)格調(diào)整值的影響，調(diào)整值由式(17)計算。

(17)

式(17)中，ΔZ19和ΔZ均為未知量，添加泥沙守恒方程(18)使其封閉[14]。

圖2 一維非結(jié)構(gòu)中心網(wǎng)格和上下游網(wǎng)格的位置

Fig.2 Location of up and down grids of one-dimensional unstructured central grid and central grid

(18)

式中i為共用鄰近網(wǎng)格，A19和Ai為單元面水平投影面積。

3.3 水沙模型的耦合與求解

在水沙模型耦合時，將計算域單獨分為流體域、泥沙固體域和水沙交界面，在流體域內(nèi)基于ANSYS-FLUENT平臺計算水動力控制方程，在固體域內(nèi)基于用戶定義函數(shù)(UDF)實現(xiàn)泥沙輸移模型和底床變形方程的數(shù)值求解，然后將位移變化量賦給水沙交界面。

水沙模型的耦合計算以上一時間步結(jié)果作為水沙交界面當前步的計算數(shù)據(jù)，即采用單向弱耦合，具體步驟如下。

(1)根據(jù)上一時間n-1步的水流模塊和泥沙模塊計算數(shù)據(jù)，執(zhí)行當前時間n步底床變形方程的顯式數(shù)值計算，通過反距離權(quán)重系數(shù)得到交界單元面上節(jié)點的位移值。

(2)根據(jù)當前時間n步的位移值，預估交界面任何位置坡度，若超過泥沙休止角，執(zhí)行沙滑模型，進而調(diào)整n步位移值給定最終更新值。

(3)根據(jù)最終位移更新值修改交界面地形，重構(gòu)計算域內(nèi)網(wǎng)格。若沖淤平衡則停止計算，否則重構(gòu)后的網(wǎng)格將為下一時間n+1步的水流模塊提供邊界條件，執(zhí)行步驟(1)和(2)，直至達到平衡。

圖3 斜坡上泥沙顆粒滑動

Fig.3 Sediment particles sliding on the slope

應(yīng)用ANSYS-FLUENT用戶自定義函數(shù)接口(UDF)實現(xiàn)泥沙模塊和底床地形變化模塊的數(shù)值求解后，得到了床面位移的更新值，通過動網(wǎng)格模型中DEFFINE _GRID_MOTION 更新和定義邊界運動。

4 模型驗證與應(yīng)用

4.1 Sumer動床圓柱沖刷

4.1.1 模型及工況

模型的可靠性和準確性驗證基于Roulund等[15]的動床圓柱沖刷試驗。在動床試驗中，橋墩的直徑為D=10 cm，距離穩(wěn)定后的水流進口為6.6 m，前方長為3.85 m的底床范圍均為動床。試驗明渠長為9.9 m，寬為3.6 m。水流行進流速為V= 0.46 m/s，水深H=0.4 m，水流邊界層厚度為0.2 m，行進流速與泥沙顆粒起動臨界流速比值為1.25，底床泥沙顆粒粒徑為d50=0.26 mm。圓柱周圍的水流特性和沖刷過程記錄采用水下微型攝像機，流速采用DANTEC 2.5 mm。

數(shù)值計算域和計算工況與文獻[14]的數(shù)值試驗保持一致，列入表2。

4.1.2 平衡沖深時地形與試驗對比

圖4為平衡沖深時沖坑形狀和Sumer試驗的對比?？梢钥闯?，通過本文的三維局部沖刷數(shù)值模擬得到最終沖坑形狀投影至平面為圓形，與試驗測量結(jié)果吻合良好。

表2 計算工況

Tab.2 Computational conditions

泥沙粒徑/mm橋墩直徑/m行進流速/m·s-1水深/m0.260.1000.460.2

圖4 平衡沖深時沖坑形狀和試驗對比

Fig.4 Comparison of the shape of the scour hole with the experiment in equilibrium depth

圖5為平衡沖深時縱軸線本文數(shù)值計算的沖深與Sumer試驗測量值和數(shù)值計算值的比較。從圖5的沖深值比較可以看出，本文數(shù)值的圓柱迎水段沖深值大于Sumer試驗和數(shù)模值，最大誤差值分別達到20%和30%，但變化趨勢一致，其原因在于本文初始給定的水流穩(wěn)定段直接影響了迎水段泥沙輸移。圓柱尾流后的沖深值與Sumer數(shù)值計算值間的誤差為5%左右。

4.1.3 沖刷坑內(nèi)水流特性與試驗對比

由于Sumer的動床試驗未測量沖坑內(nèi)水流特性，與Unger等[16]的試驗測量結(jié)果定性比較分析軸對稱面圓柱迎水端沖坑內(nèi)不同時刻的流場特征，如 圖6 所示。可以看出，沖坑發(fā)展時坑內(nèi)流場表現(xiàn)的特征基本與試驗測量結(jié)果吻合。在初始時刻，如圖6(a)所示，近底床附近的渦旋和下潛水流以及近自由水面的上升趨勢均有捕捉；在沖刷發(fā)展一段時間后，如圖6(b)所示，沖坑從圓柱邊界逐漸向后發(fā)展，形成了一定的坡面，整體坡面角度均小于泥沙休止角，且坑內(nèi)水流為明顯的下潛水流，與Unger試驗測量結(jié)果一致，Unger分析后指出此階段的下潛水流偏轉(zhuǎn)至圓柱橫向兩側(cè)，加劇此處泥沙的沖刷；當沖刷繼續(xù)發(fā)展時，如圖6(c)所示，沖坑內(nèi)形成了明顯的渦旋，近床面附近的水流下潛趨勢也明顯。此外，從自由水面在整個沖刷演變過程中的發(fā)展可以看出，靠近圓柱附近的水面上升趨勢未 改變。

圖5 縱軸線數(shù)值計算的沖深與Sumer試驗測量值比較

Fig.5 Comparison of the numerical depth of the longitudinal axis with the experimented by Sumer

圖6 軸對稱面數(shù)值計算的沖坑內(nèi)流場(左)與Unger試驗(右)比較

Fig.6 Comparison of numerical calculation(left)and Unger’s experiment(right)of the flow field of axisymmetric plane in the scour hole

根據(jù)沖坑內(nèi)水流特性與試驗的定性比較，可以認為本文建立的三維數(shù)值局部沖刷模型能準確模擬沖坑內(nèi)水流發(fā)展。

4.1.4 三維沖刷坑的演變結(jié)果分析

圓柱局部沖刷在不同時刻的三維沖坑地形如圖7所示?？梢钥闯?，地形演變特征如下。在初始發(fā)展階段，即T≤50 s時，圓柱上游的沖坑形狀接近半圓形，斜坡斜率小于或等于泥沙休止角，在圓柱下游的泥沙淤積成類似坎的形狀；隨著沖刷持續(xù)發(fā)展至T≤1230 s，沖坑內(nèi)水流表現(xiàn)出圖6(a)的特性后逐漸演變?yōu)閳D6(b,c)所示的特性，由于下潛水流偏移至圓柱兩側(cè)導致沖坑橫向擴展和縱向逐漸沖深，最終形成了圓形沖坑，坑內(nèi)的泥沙帶走后大面積淤積為舌狀，當水流演變?yōu)橄聺撍骱蜏u旋共同作用后，水流持續(xù)沖刷沖坑底部，使得沖坑橫向和縱向范圍擴大，且伴隨斜坡坍塌形成倒錐形坡度，坑后淤積泥沙向圓柱橫向和下游擴散；當沖刷發(fā)展至T=2415 s時，即將達到穩(wěn)定狀態(tài)，此時沖坑范圍擴大趨勢減小；當沖刷發(fā)展至T=3600 s時，達到穩(wěn)定沖深階段，深度增加值和范圍基本保持穩(wěn)定，床面泥沙沖淤達到平衡狀態(tài)。此外，在整個過程中能較為清晰地觀察到初始時下游和水流入口段形成的堆積泥沙(沙丘)逐漸遷移發(fā)展的演變。以上總結(jié)的演變特征與Sumer試驗過程中觀察到的演變特征一致。

圖7 不同時刻的沖深地形

Fig.7 Topography of scour hole at different times

經(jīng)過上述動床沖坑演變、沖深變化和坑內(nèi)水流特性可以看出,本文建立的三維局部沖深數(shù)值模型的準確度能滿足要求。

4.2 Michiue動床丁壩局部沖刷試驗

4.2.1 模型及工況

采用Michiue等[17]在明渠水槽中的非淹沒正交丁壩局部沖刷試驗進一步驗證數(shù)值模型準確度。試驗水槽長為7 m、寬為0.4 m，坡度為1/300；流量為0.004 m3/s，平均水深為2.85 cm，行進流速為 0.351 m/s；泥沙為均勻沙，中值粒徑為0.6 mm。丁壩尺寸長為0.1 m，寬為0.01 m。數(shù)值計算的區(qū)域長為8.0 m、寬為0.4 m、高為0.3 m，如圖8所示。

4.2.2 平衡沖深時地形與試驗對比

圖9比較了丁壩沖刷數(shù)值計算地形與Michiue的試驗測量值?？梢钥闯觯诙伪沉鱾?cè)，數(shù)值計算的沖刷坑形態(tài)與試驗測量基本一致，計算的最大沖刷深度平均達到了5 cm；在丁壩壩頭位置，沖刷坑的外圍形態(tài)均表現(xiàn)為長舌狀；但數(shù)值計算的最大沖刷深度僅位于壩頭，試驗測量值則延伸至迎流側(cè)的上游約15 cm處，其原因是由于數(shù)值計算模型中添加的沙滑模型導致斜坡面滿足泥沙休止角后就處于穩(wěn)定狀態(tài)，不再受到?jīng)_刷，使得此處與實測值存在一定偏差。圖9中，丁壩的計算最大沖刷深度約為8 cm，試驗測量值為9 cm，兩者的誤差為11%。此外，數(shù)值計算結(jié)果沿y方向的擴展比測量值大，這種趨勢與文獻[14，18]數(shù)值計算的丁壩沖刷時得到的沖刷范圍明顯比實測范圍大的結(jié)果 相似。

圖8 數(shù)值計算域示意

Fig.8 Schematic diagram of numerical calculation area

4.2.3 三維沖刷坑的演變結(jié)果分析

圖10為丁壩周圍局部沖刷演變示意。可以看出，丁壩沖刷至時間T=50 s時的時間段內(nèi)為初始沖刷發(fā)展階段，且初始位置始于迎流側(cè)和壩頭，帶走的泥沙分別向壩頭的沖坑外圍和壩后淤積成由兩道梁匯成的舌狀體；隨著水流持續(xù)沖刷至 550 s 時，沖坑范圍持續(xù)性的擴展變大，壩后舌狀堆積體逐漸往邊墻附近移動，丁壩迎流和背流側(cè)的沖坑均往外后退擴展至邊墻；當沖刷至1050 s時，壩后的淤積形態(tài)和壩頭沖坑地形形態(tài)基本不再發(fā)生變化，此時迎流側(cè)靠近邊墻的沖坑斜坡面繼續(xù)沖深；當沖刷至2300 s和3600 s時，丁壩沖刷的最大沖刷深度誤差接近0.2%，沖坑整體形態(tài)保持穩(wěn)定，即沖刷達到穩(wěn)定沖深階段。

圖9 沖刷數(shù)值計算地形與Michiue試驗實測值比較(單位:cm)

Fig.9 Comparison of numerical calculation and Michiue’s experiment of the topography of scour hole(unit:cm)

圖10 丁壩周圍不同時刻局部沖刷坑的地形

Fig.10 Topography of local scour hole around dike at different times

5 結(jié) 論

基于有限體積法和非結(jié)構(gòu)化的計算網(wǎng)格建立了包括三維非恒定水流、泥沙輸移和河床變形方程等控制方程的局部沖刷三維數(shù)學模型。應(yīng)用ANSYS-FLUENT軟件模擬三維非恒定水流，利用用戶自定義函數(shù)接口(UDF)二次開發(fā)實現(xiàn)泥沙輸移和底床變形方程的數(shù)值計算，并結(jié)合動網(wǎng)格技術(shù)重構(gòu)變形網(wǎng)格，最終實現(xiàn)了局部沖刷的三維數(shù)值模擬。主要結(jié)論如下。

與文獻[14]的圓柱動床局部沖刷試驗以及文獻[16]進行的丁壩局部沖刷試驗結(jié)果對比后驗證了該三維模型的可靠性，并成功模擬了局部沖刷過程中地形的變化、最大沖深位置和平衡沖深時的最終地形。

以切應(yīng)力觀點為基礎(chǔ)建立的三維局部沖刷數(shù)值模型，因數(shù)值計算過程中切應(yīng)力較易獲取，具有簡易性和較強的可靠性；在利用有限體積法開發(fā)類似的局部沖刷數(shù)值模型時，單元體泥沙通量重構(gòu)和坡度影響下,泥沙顆粒起動、輸沙率計算以及網(wǎng)格重構(gòu)過程中坡度修正等均影響數(shù)值模型的精度。

在地形變化導致網(wǎng)格更新時，ANSYS-FLUENT提供的動網(wǎng)格模型基本能滿足體型較為簡單的工程(如本文圓柱橋墩沖刷和丁壩沖刷等)局部沖涮大變形區(qū)域網(wǎng)格更新需求。但對于體型復雜的工程問題已不再適應(yīng)，應(yīng)開發(fā)相應(yīng)的網(wǎng)格重構(gòu)替代模型。