泥沙沖刷模擬的FLOW-3D數值方法研究

黃佳麗 ，倪福生 ，顧磊

（1.河海大學疏浚技術教育部工程研究中心，江蘇 常州 213022；2.河海大學機電工程學院，江蘇 常州 213022）

0 引言

泥沙沖刷是很多水利工程中常見的現象，常見于樁柱沖刷、管道局部沖刷及水射流沖刷等。黃雄合等[1]利用動床模型水槽試驗，研究了樁柱周圍局部床面由平整發展到沖刷坑的過程，獲得了各工況下沖刷坑的形態參數；姜萌等[2]采用物理模型試驗方法，探索了泥沙中值粒徑、立管-樁傾斜角度及水深等對海洋平臺立管系統底部沖刷坑的影響。張浩等[3]在粗沙沙床條件下對二維垂向淹沒射流沖坑深度進行了實驗研究，研究了射流速度對動態坑深的影響；顧磊等[4]探究了雙股射流沖刷時噴嘴間距對沖刷坑形、沖坑深度和沖坑截面積的影響。

通過對前人的研究以及相關試驗過程的分析發現試驗研究存在很多缺點：在試驗前對沖刷的泥沙的準備有很高的要求，每次都需要耗費大量的時間；試驗過程中無法觀測試驗現象，尤其是泥沙粒徑很小時，射流沖刷流場很渾濁，無法對其內部機理進行捕捉；沖刷結束后需要緩慢將水排放，不僅耗時而且還對沖坑有一定的影響；最后在對沖坑進行測量時，沖坑表面會有很多的泥沙堆積，對其靜態坑深的結果有影響，同時由于觀測不便，動態沖坑深度的變化無法準確的進行測量。

近年來，在與泥沙沖刷的有關數值模擬中，FLOW-3D軟件應用較為廣泛。與試驗相比，數值模擬方法具有省錢、省力和靈活性較大等優勢，由于試驗參數的改變方便，對于大量重復試驗能夠輕易完成[5]。周麗丹等[6]采用泥沙沖刷模型和湍流模型k-ε相結合，對均勻來流海底管道沖刷過程進行分析，發現了沖刷坑的形成原因及沖刷坑最大深度發生位置；王建軍等[7]采用FLOW-3D軟件基于泥沙沖刷模型進行了數值模擬，發現數值模擬結果與試驗結果較為吻合；顧磊等[8]應用FLOW-3D軟件研究了雙噴嘴射流沖刷砂床時沖坑形態的發展過程以及沖坑形態、深度和面積隨噴嘴間距的變化規律，并通過模型試驗對其進行了驗證。

本文在前人研究的基礎上發現：數值模擬軟件FLOW-3D在泥沙沖刷的模擬中有較完善的模型和方法，但是在泥沙沖刷領域對其泥沙模型關鍵參數設置的研究較少，大多采用經驗值或相互參考，較多文章中參數幾乎一致，沒有將模型參數與試驗可以獲得參數進行較好的關聯。所以本文采用FLOW-3D數值模擬的手段從射流沖刷泥沙出發，對泥沙模型中的關鍵參數進行了探究，以期為該軟件在泥沙沖刷模擬上的應用提供一定的參考。

1 理論模型

1.1 控制方程

連續性方程

式中：VF為流體流動部分的體積分數；ρ為流體密度；RDIF為湍流擴散項；RSOR為質量源項；u、v、w分別為在坐標方向（x、y、z）上的速度分量；Ax、Ay、Az分別表示流體在x、y、z方向的投影面積；R和ξ為與坐標系統的選擇相關的系數。

動量方程

式中：Gx、Gy、Gz為流體重力加速度；fx、fy、fz為黏滯加速度；bx、by、bz為流體通過多孔介質或導板的能量損失。

1.2 泥沙模型

在泥沙模型中，泥沙在水中以懸移質和推移質兩種狀態存在。由于單個泥沙顆粒周圍的流體動力和交界面處的邊界層難以計算，所以須使用經驗模型。

坡面法線與重力方向的夾角懸移質計算模型中，泥沙的上升速度由式（3）計算：

式中：ulift為用來計算由底沙轉化為懸沙的泥沙量；α為挾帶系數；ns為沙床面的外法線方向；d*為泥沙顆粒的無量綱粒徑；θ為床面希爾茲數；θcr為臨界希爾茲數；g為重力加速度；ds為泥沙顆粒的直徑；ρs為泥沙的密度。

床面希爾茲數θ由床面剪切應力計算：

式中：τ為床面剪切應力。

臨界希爾茲數θcr基于Mastbergen和Von den Ber[9]提出的經驗模型，以及Shields-Rouse方程[10]來預測。

推移質計算基于Meyer-Peter和Muller[11]的計算模型，床面的單寬體積推移質輸沙率由該模型預測可得，公式如下：

式中：β是推移質系數；Φ是無量綱的推移質輸沙率；Φ與單寬體積推移質輸沙率qb的關系是：

由式（6）可以計算出每個單元中各個時刻的單寬體積輸沙率。

2 網格質量

網格是數值模擬計算的基礎，網格劃分決定了計算域的劃分和控制方程的離散程度，其尺度決定了計算的穩定性和準確度。因此在數值模擬中對于網格的探索十分必要。在使用FLOW-3D進行水射流沖刷泥沙的數值模擬時，對網格的探索如下。

由于建立模型的簡單，可以直接在軟件內進行繪制，網格覆蓋整個模型區域，為保證所有結構都能被識別，在關鍵部位增加網格線，進行局部加密。

在所建模型中，噴嘴的結構是圓柱體，而FLOW-3D中的網格形式為結構式有限差分網格，單元格呈現為矩形棱柱。不同于CFD常見的貼體網格，FLOW-3D的網格對復雜幾何的描述精度較差，必要的局部加密才能實現較好的描述，網格覆蓋整個模型但因其全結構式，有很高的計算效率，累計的數值誤差較小。

將尺寸相對于整體結構較小的噴嘴進行加密處理。圖1為3種相鄰網格尺寸比情況下，劃分的網格經過favor檢查之后的噴嘴內徑與外徑的識別情況，深色區域表示所識別的固體材料。

圖1 不同網格尺寸比的噴嘴識別情形Fig.1 Nozzle recognition under different mesh size ratios

由圖1中可以看出相鄰網格尺寸比越小，識別的噴嘴結構越圓。在相鄰網格尺寸比較大時，噴嘴只在網格的節點以方形結構顯示，所以在進行網格加密時，同一坐標方向上相鄰兩個結構的網格寬度之比，應控制在一定的范圍之內，一般要求不超過1.25。但考慮計算效率等因素，經過對各種網格尺寸比的探索，發現比值在1.0～2.0以內時，計算結果的數值沒有明顯的誤差。此外單個網格在3個坐標方向上的寬度比例也應該保持在合適的范圍以內，一般不超過3.0，盡可能接近一致。

因此在對實際問題進行數值計算時，其網格比例通常是不確定的，需要通過不斷的調試和比較，才能獲得吻合具體問題的網格。

本文研究模型較為簡單，所得出的加密尺寸僅供參考，其他不同情況下的模型網格還是需要大量探索。

3 泥沙模型關鍵參數

在水射流過程中，沖刷為一種隨機現象，泥沙會呈現出懸浮、擴散、遷移和沉降等狀態，FLOW-3D軟件中的泥沙模型可較好地預測泥沙的各種狀態，模擬出與試驗吻合的泥沙運動狀態。泥沙模型參數設置對泥沙的懸浮、擴散以及堆積情況有很大影響，其中影響最大的是泥沙臨界體積分數Cv、泥沙的水下休止角φ、泥沙挾帶系數α和推移質系數β。由于研究的泥沙粒徑不同，各種參數的具體設置也不同。采用單因素分析法，以0.09 mm粒徑的沙為例，每次只改變其中1個參數，觀察該參數對模擬過程的影響。

3.1 泥沙臨界體積分數Cv

泥沙的臨界體積分數Cv對泥沙沖刷形態有很大影響，其設置范圍為（0，1）。分別設置體積分數為0.1，0.3和0.9，其他參數保持默認，模擬結果如圖2所示。

圖2 不同泥沙體積分數的射流沖刷濃度分布Fig.2 Distribution of jet scouring concentration under different maximum packing fractions

由圖2中可看出，泥沙臨界體積分數值不同時，射流沖坑深度相差不大，但是泥沙懸浮量有很大不同，臨界體積分數為0.1時懸浮泥沙擴散最高，當其值增大時懸浮泥沙擴散高度降低，且沿著沙床方向移動，運動過程呈現為懸移質與推移質的交換過程。即臨界體積分數越大，懸浮的泥沙量越少?？偨Y前人關于FLOW-3D的研究[11]以及在對射流沖刷泥沙試驗中需要研究的實際參數的比較，發現泥沙臨界體積分數與試驗中泥沙的孔隙率有關，泥沙臨界體積分數=1-泥沙孔隙率。

當研究的泥沙粒徑不同時，泥沙臨界體積分數不同。對于不均勻泥沙，其值會較大；而對于均勻泥沙或泥沙顆粒形狀不規則，其值會較小。

3.2 泥沙的水下休止角φ

在射流沖刷完成后，沖坑呈現為有一定角度的斜坑，在模擬時，可以通過設置泥沙模型中水下休止角的參數模擬沖坑角度。在對該參數進行設置時發現，休止角設置為某角度最終呈現的模擬效果就是該角度，十分直觀。圖3為分別將參數值設置為32°和45°時的沖坑截面圖。

圖3 不同水下休止角的射流沖坑形狀Fig.3 The shape of jet scouring pit at different underwater angles of repose

從圖3的沖坑角度來看，圖3（a）設置為32°，模擬結果顯示休止角為31°；圖3（b）設置為45°，模擬結果顯示休止角為45°，可見最終沖坑的休止角與所設置的參數值一致，但如此并不符合實際情況，試驗時得到的沖坑不是一個確定的角度。所以在設置泥沙水下休止角時要先對所試驗的沙粒進行水下休止角測量，同時考慮泥沙粒徑和水下斜坡效果的影響，如此才能更好地貼合實際情況。

3.3 泥沙挾帶系數α

泥沙挾帶系數α控制該沉積物受力大于臨界剪切應力時侵蝕的速率。挾帶系數可以用于縮放沖刷率。零值完全關閉挾帶模型。圖4為挾帶系數分別為0.001、0.018和0.4時的射流沖刷濃度分布情況。

在進行模擬時，只改變挾帶系數的值，其他參數保持為默認值，從圖4可以看出，圖4（a）中沖刷坑深以及沖刷寬度都很小，懸沙沒有沉降而是隨著水流的上升全部上升，此情況不符合實際情況；圖4（b）中，沖刷坑深與沖坑寬度比圖4（a）有明顯增加，懸沙在渦流挾帶下在兩側交替增多，兩側懸沙形態呈現一定的對稱性，有沉積現象且泥沙向四周擴散；圖4（c）中沖坑深度與沖坑寬度遠大于前兩種情況，泥沙起動之后直接向四周擴散，沒有出現明顯的泥沙隨渦流懸浮運動的狀況。綜合3種情況可以看出，挾帶參數的選擇對于泥沙運動模擬情形有很大的影響。在模擬過程中，挾帶模型的速度方向是垂直于沙面向外的，該數值越大對底床沖刷的速度越快，隨沖刷起動的沙更多。

圖4 不同挾帶系數的射流沖刷濃度分布Fig.4 Distribution of jet scouring concentration under different entrainment coefficients

3.4 推移質系數β

推移質系數β在床載荷傳輸速率方程中使用，并控制該床載荷傳輸發生在剪切應力比臨界剪切應力更高時的速率。不同推移質系數下射流沖刷濃度分布如圖5所示。

圖5 不同推移質系數的射流沖刷濃度分布Fig.5 Distribution of jet scouring concentration under different bed load coefficients

由圖5可以看出推移質系數小時，泥沙懸浮量很小，在沙坑兩邊有輕微堆積泥沙現象，當推移質系數較大時，沖刷沙坑中泥沙懸浮明顯，且沙坑兩邊泥沙堆積較為明顯?？梢钥闯觯和埔瀑|模型中泥沙的運動方向是沙面附近的水流方向，推移質系數值越大，沉積物沿河床的傳輸量越大，研究發現[11]從低運輸的5.0到非常高的沙運輸的13.0的值，沙子和礫石的典型報告值為5.7，在對推移質系數進行設置時可以此為參考。

4 結語

本文應用FLOW－3D軟件中的泥沙沖刷模型，對垂直水射流沖刷泥沙進行了數值模擬。本文的研究可以對FLOW-3D在泥沙沖刷的應用提供一定的參考。

1）FLOW-3D軟件在用于水射流沖刷的數值模擬可以與試驗結果有較好的吻合。

2）網格加密時，需要將相鄰網格比以及網格本身3個坐標方向上的寬度比保持在合適的比例，且需要經過不斷的調試和比較才能獲得計算時間短且符合試驗的結果。

3）在泥沙參數設置中，泥沙的臨界體積分數越大，懸浮的泥沙量越少；水下休止角則需要進行試驗驗證；挾帶系數越大，隨沖刷起動的沙更多；推移質系數越大，沉積物沿河床的傳輸量越大。