泥沙局部沖淤二維數值模擬仿真

張芝永，拾兵

(中國海洋大學工程學院，山東青島266100)

泥沙沖淤是河流、海洋中比較普遍的一種自然現象，在近床面位置存在結構物情況下，床面泥沙運動劇烈，局部地形變化強烈，局部沖淤問題嚴重.對于此類問題，國內外許多學者主要是通過模型實驗手段來進行研究分析，其研究對象主要包括橋墩［1］、丁壩［2-3］、海底管道［4-5］等結構物周圍的局部沖淤問題.

隨著計算機技術的發展，局部沖淤的數值模擬方法越來越引起研究者的重視.Ouillon等［6］采用標準的三維k-ε紊流模型模擬了丁壩周圍的流場，并對平整床面時的床面切應力沿河床的分布與最終沖刷坑的床面形態進行了定性比較.彭靜等［7］將線性與非線性三維k-ε模型應用于丁壩繞流模擬，并比較了各種模型的模擬效果.Brφrs［8］提出了新的二維沖淤模型采用標準的k-ε紊流模型和同時考慮懸移質和推移質輸運，對海底管道局部沖刷過程進行了模擬.Liang和Cheng［9-10］基于有限差分法發展了一個可以精確地預測海底管道沖刷時間歷程的二維模型.Liu［11］利用VOF和動網格技術對水平射流情況下的底床沖淤進行了數值模擬，該模型考慮了水面的變化情況，沖坑變化過程模擬值與試驗值吻合較好.韋燕機［12］利用OpenFOAM和動網格技術對樁周局部沖刷進行了三維數值模擬.Zhi-wen Zhu［13］利用網格技術對橋墩沖刷進行了數值模擬，驗證了網格較好的地形邊界擬合效果.綜上所述，以上大部分模型均是通過動網格技術來模擬泥沙的局部沖刷，這說明動網格技術要比其他如兩相流技術來模擬泥沙更為準確.但上述模型也存在一些不足，如大部分局部沖淤模擬均是通過非通用計算程序來完成的，這限制了其數值模擬方法的推廣.

FLUENT是一款目前較為流行的商用CFD軟件，它在航空航天，能源、水利、環境等領域得到了廣泛的應用，是目前全球功能最強大的計算流體力學商用軟件.該軟件在水動力計算模塊有著廣泛的通用性，但該軟件并沒有利用動網格來計算泥沙沖淤的模塊.有鑒于此，作者充分利用FLUENT的二次開發接口，即自定義函數功能，通過C++編寫了泥沙輸運模塊，并將其嵌入到FLUENT中，實現了泥沙局部沖淤的二維數值模擬.

1 水動力模型

1.1 控制方程

水動力模塊中，其控制方程為連續性方程和動量方程:

式中:ui為流體速度在i方向上的分量;ρ是流體的密度;μ是流體粘度;p是作用在流體微元上的壓力;v是流體的運動粘滯系數，Sij是平均應變速率張量.ui'是i方向速度脈動值為雷諾應力張量，vT紊流運動粘滯系數，δij是克羅內克爾符號.

1.2 湍流模型

對于不可壓縮流體且不考慮用戶自定義的源項時，標準k-ε模型方程為k方程:

ε方程:

式中:Gk是由于平均速度梯度引起的湍動能k的產生項;C1ε和C2ε為經驗常數，取C1ε=1.44，C2ε=1.92;σk和σε分別是與湍動能k和耗散率ε對應的 Prandtl數，取 σk=1.0，σε=1.3.對于床面近底層，采用壁面函數法將壁面上的物理量與湍流核心區求得物理量聯系起來.

2 泥沙輸運模型

2.1 推移質輸沙模型

床面泥沙的起動可以通過臨界希爾茲數來確定，其希爾茲數的表達式為

式中:τ為床面切應力，ρ為水的密度，g為重力加速度，s為泥沙比重，d50為泥沙中值粒徑.

水平床面上泥沙臨界起動希爾茲數為

其中，D*是無量綱泥沙顆粒尺寸，定義為

在有坡度的床面上泥沙臨界起動希爾茲數 按式(10)進行修正［14］:

式中:α是砂床面和水平面的坡角，規定上坡為正角，下坡為負角;φ是泥沙的休止角.

平面單寬體積輸沙率q0由式(11)求得

推移質單寬體積輸沙率可用下式表示:

式中:q0為平面單寬體積輸沙率，τ為床面剪應力，τx為床面剪應力x方向分量，h為床面高程，經驗系數C取值范圍為 1.5 ～2.3［8］.

2.2 床面變形模型

根據輸沙平衡，床面高程h的變化可以用沖淤方程表示為

式中:p0為床沙孔隙率，qbx為x方向推移質單寬體積輸沙率.

3 動邊界處理方法

3.1 動邊界處理方法

動網格模型可以用來模擬流場形狀由于邊界運動而隨時間改變的問題.其廣泛應用于活塞、閥門及柔性體等有運動邊界工況的模擬.床面的沖淤變形形式類似于柔性體變形，需要對各個節點的運動情況進行描述.因此本文采用FLUENT的DEFINE_GRID_MOTION宏命令來給出邊界節點的運動方式.在床面節點位置更新后，需要對區域內部網格進行調整.FLUENT提供了3種網格更新方式:1)光順網格法，2)動態層網格法，3)局部網格重畫法.非結構化網格比結構化網格更加適用于不規則邊界的變形問題，因此本模型采用非結構化網格，其網格更新方法為光順網格法和局部網格重畫法相結合的方式，這2種方法的網格更新效果可見圖1.圖1為下文算例中水流下坡面處局部網格變化情況.從圖中可以看出，在底部邊界發生變化后，區域內網格進行重新劃分，各區域網格尺寸尺度基本保持不變，靠近下部變形邊界的網格仍然比較密集，這對于確保近壁面流場計算的準確性至關重要.

圖1 動網格更新效果Fig.1 Mesh deformation

3.2 沙滑模型

由于泥沙傳輸的復雜性和與流場相互作用的非線性，在海床高程更新過程中會發生反射，導致網格畸變，造成不存在的沖淤形態和數值的不穩定.為解決此問題，采用Liang提出的沙滑模型［10］:當局部沖淤斜坡角度超過泥沙休止角時，對相應的海床節點進行調節，使沖淤斜坡角等于休止角.

4 數值解法及計算流程

水動力控制方程及紊流方程使用基于單元中心的有限體積法離散，利用非穩態求解器進行求解，瞬態項采用二階隱格式，對流項采用二階迎風差分格式，擴散項采用中心差分格式;壓力與速度的耦合使用SIMPLE算法.

床面變形方程式(13)通過有限差分法進行離散，其離散后形式為

式中:p代表床面節點序號;n為當前時間步;n+1為下一時間步;Δtb為床面更新時間步長.

為節省計算時間，流場計算的時間步長和床面更新的時間步長采用不同值，床面更新時間步長Δtb要遠大于流場計算時間步長Δtflow.其數值模擬過程可概括為以下幾個過程:

1)計算多個步長的速度、壓力、湍流數及切應力等;

2)利用最后一步得到的床面切應力數據，調用泥沙輸運模塊，計算推移質輸沙率和床面高程變化情況;

3)采用動網格技術，改變床面節點位置，同時重新調整計算區域內部網格;

4)返回步驟1)，重新上面的計算，直至到達指定時間.

5 模型驗證

為驗證本文所建模型的可靠性，通過2個算例對其進行了驗證.

5.1 水渠沖淤

Van Rijn［15］對有坡度水渠的水流情況進行了試驗研究，李昌良［16］在其基礎上又對其沖淤情況進行了進一步的研究.其試驗布置方案如圖2所示.左邊進口平均流速0.51 m/s.水深0.39 m，沙床泥沙平均粒徑d50=0.16 mm.在數值模擬中，其計算區域與試驗布置相同，左邊進口設置為速度入口，右邊為自由出流邊界，水面邊界設置為對稱邊界.底面邊界設置為wall邊界，其粗糙度為2.5d50.整個區域采用非結構化網格進行離散.

圖2 試驗布置Fig.2 Layout of experiment

圖3列出了在初始地形情況下4個垂直斷面的流速分布情況，從圖中可以看出數值模擬結果與試驗結果是較為一致的，這說明了數值模擬中水面采用剛蓋假設的合理性，同時為進一步的泥沙沖淤過程研究奠定了基礎.

圖4為不同時段的網格情況及床面地形變化情況，可以發現在床面地形發生變化后，采用光順網格法和局部網格重畫法更新后的非結構化網格很好地擬合了床面地形的變化，在變化過程中，網格的疏密程度也得到了較好的保持.圖5給出了3 h后床面沖淤情況.從圖中可以看出，在溝內由于水深的增加，水流流速減小，進而導致床面切應力減小，泥沙在此淤積.在水流下坡區域，床面切應力減小，泥沙淤積較為嚴重，而在水流上坡區域，下游流速加大，床面切應力增大，從而導致該區域沖刷嚴重.數值結果與前人結果比較一致，準確地反映了溝內床面的變化規律.

圖3 各斷面流速分布情況Fig.3 Velocity distributions at different crossvsections

圖4 不同時刻網格更新情況Fig.4 Mesh deformation in different time

圖5 3 h后床面地形Fig.5 Bed profile at t=3.0 h

5.2 海底管道局部沖刷

本算例對近床面的海底管道局部沖刷的沖刷歷程進行詳細的分析.

Liang在文獻［17］中對有間隙海底管道的局部沖刷進行了詳細的數值模擬研究.在這里采用其中一計算算例進行了數值計算，算例中管道直徑D為10 cm，管道與床面之間垂直距離G為0.5D，來流流速為0.5 m/s，泥沙中值粒徑為 0.36 mm，水深 3.5D.選擇距管道中心上下游各20D距離的區域為計算區域.左邊界為速度進口邊界，右邊界為自由出流邊界，上部邊界采用對稱邊界，下部邊界為wall邊界(如圖6所示).整個區域采用非結構化網格進行離散，在近壁面處和管道周圍進行加密.

圖6 計算區域Fig.6 Computational domain

圖7 不同時刻x方向流速等值線(單位:m/s)Fig.7 Contours of x veolicity in different time(unit:m/s)

由于動網格計算耗時較長，本文只模擬了80 min內的床面沖刷過程.圖7列出了沖刷過程中，不同時刻的x方向流速分布圖，該圖同樣直觀地描述了床面的沖刷變形情況.在t=5 min時，由于初始沖刷階段間隙內流速較大，管道局部沖刷程度較為嚴重，管道正下方出現沙坑，下游區域出現沙丘.在t=15 min時，上游沖坑深度加大，下游區沙丘逐漸向下游移動并并且沙丘高度逐漸減小.在t=80 min時，下游沙丘消失，整個過程與實際的沖刷過程規律是較為一致的.從圖中還可以發現，管道下游并未出現渦旋脫落現象，這可能是由網格尺度的限制以及啟動動網格模型而導致的計算精度下降等原因造成的.而在實際試驗中，管后區域是存在渦旋脫落現象的.

Liang［17］的研究結果表明，有渦旋脫落時，管道后尾流區范圍比較大而且變化比較劇烈，大范圍的尾流區使得更多水流通過管道下方孔道流向下游，進而導致近床面流速加大，管下和管后床面切應力也相應增大，同時，渦旋的不斷脫落又使得管道下方及后方床面的切應力波動變化，波動幅值有可能會是平均值的數倍.因此，有渦旋釋脫落時的沖刷深度要大于無渦旋脫落情況.圖8列出了t=80 min時，Liang的模擬結果與本文模擬結果的對比情況.從圖中可以看出，管道上游的床面沖刷模擬結果較為接近，管道下游區域床面沖刷模擬結果與Liang的無渦旋脫落時的沖刷結果比較接近，而與有渦旋脫落時的模擬結果則有較大偏差.其原因在于本文模擬過程中并未出現渦旋脫落現象，從而造成管道下方及后方床面切應力計算值偏小，沖刷程度較小.

圖8 t=80 min時床面地形Fig.8 Bed profile at t=80 min

6 結論

通過對商用CFD軟件FLUENT的二次開發，建立了局部沖淤二維數值模型.該模型的泥沙輸運模型通過FLUENT的DEFINE_GRID_MOTION宏命令嵌入，床面邊界的沖淤變化利用動網格技術來實現.

數值計算的結果表明，泥沙沖淤模型可以比較準確的預測不同條件下泥沙沖淤過程，可用于模擬以推移質輸沙為主的二維泥沙沖淤問題，具有廣闊的應用前景.但該模型同樣存在一些問題，如對于有渦旋脫落時的泥沙沖淤模擬偏差較大以及計算時間過長等問題.這都有待于進一步優化和研究.

［1］MELVILLE B W.Pier and abutment scour:integrated Approach［J］.Journal of Hydraulic Engineering，1997，123(2):125-136.

［2］彭靜，河原能久.丁壩群近體流動結構的可視化實驗研究［J］.水利學報，2000，31(3):44-47.PENG Jing，YOSHIHISA Kawahara.Visualization of flow structure around submerged spur dikes［J］.Journal of Hydraulic Engineering，2000，31(3):44-47.

［3］ROGER A，KUHNLE C V，ALONSO F，et al.Local scour associated with angled spur dikes［J］.Journal of Hydraulic Engineering，2002，128(12):1087-1093.

［4］SUMER B M，JENSEN H R，FREDS?E J.Effect of leewake on scour below pipelines in current［J］.J Waterway，Port，Coastal and Ocean Engineering，1988，114(5):599-614.

［5］CHIEW Y M.Prediction of maximum scour depth at submarine pipelines［J］.Journal of Hydraulic Engineering，1991，117(4):452-466.

［6］OUILLON S，DARTUS D.Three-dimensional computation of flow around groyne［J］.Journal of Hydraulic Engineering，1997，123(11):962-970.

［7］彭靜，河源能久.線性與非線性紊流模型及其在丁壩繞流中的應用［J］.水動力學研究與進展:A輯.2003，18(5):589-594.PENG Jing ，KAWAHARA Yoshihisa.Application of linear and non-linear turbulent models in spur dike flow［J］.Chinese Journal of Hydrodynamics，2003，18(5):589-594.

［8］BR?RS B.Numerical modeling of flow and scour at pipelines［J］.Journal of Hydraulic Engineering，1999，125(5):511-523.

［9］LIANG D F，CHENG L，LI F.Numerical modelling of scour below a pipeline in currents.part I:flow simulation［J］.Coastal Engineering，2004，52(1):25-42.

［10］LIANG D F，CHENG L，LI F.Numerical modelling of scour below a pipeline in currents.partⅡ:Scour simulation［J］.Coastal Engineering，2004，52(1):43-62.

［11］LIU Xiaofeng，GARCíA M H.Three-dimensional numerical model with free water surface and mesh deformation for local sediment scour［J］.J Waterway，Port，Coastal and Ocean Engineering，2008，134(4):203-217.

［12］韋雁機，葉銀燦.床面上短圓柱體局部沖刷三維數值模擬［J］.水動力學研究與進展:A輯，2008，23(6):655-661.WEI Yanji，YE Yincan.3D numerical modeling of flow and scour around short cylinder［J］.Chinese Journal of Hydrodynamics，2008，23(6):655-661.

［13］ZHU Zhiwen，LIU Zhenqing.CFD prediction of local scour hole around bridge piers［J］.Journal of Central South University of Technology，2012，19(1):273-281.

［14］ALLEN J R L.Simple models for the shape and symmetry of tidal sand waves:statically stable equilibrium forms［J］.Marine Geology，1982，48(12):31-49.

［15］VAN RIJN L C.Principles of sediment transport in rivers，estuaries and coastal seas［M］.Amsterdam:Aqua Publications，1993:1245-1247.

［16］李昌良.泥沙運動與底床變形的數值模擬［D］.青島:中國海洋大學，2008:58-61.LI Changliang.The numerical simulation of sediment transport and bed evolution［D］.Qingdao:Ocean University of China，2008:58-61.

［17］LIANG Dongfang，CHENG Liang，YEOW K.Numerical study of the Reynolds-number dependence of two-dimensional scour beneath offshore pipelines in steady currents［J］.Ocean Engineering，2005，32(4):1590-1607.