泥沙局部沖淤二維數(shù)值模擬仿真

張芝永，拾兵

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

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

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

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

1 水動力模型

1.1 控制方程

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

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

1.2 湍流模型

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

ε方程:

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

2 泥沙輸運模型

2.1 推移質(zhì)輸沙模型

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

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

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

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

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

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

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

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

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

2.2 床面變形模型

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

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

3 動邊界處理方法

3.1 動邊界處理方法

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

圖1 動網(wǎng)格更新效果Fig.1 Mesh deformation

3.2 沙滑模型

由于泥沙傳輸?shù)膹碗s性和與流場相互作用的非線性，在海床高程更新過程中會發(fā)生反射，導致網(wǎng)格畸變，造成不存在的沖淤形態(tài)和數(shù)值的不穩(wěn)定.為解決此問題，采用Liang提出的沙滑模型［10］:當局部沖淤斜坡角度超過泥沙休止角時，對相應的海床節(jié)點進行調(diào)節(jié)，使沖淤斜坡角等于休止角.

4 數(shù)值解法及計算流程

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

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

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

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

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

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

3)采用動網(wǎng)格技術(shù)，改變床面節(jié)點位置，同時重新調(diào)整計算區(qū)域內(nèi)部網(wǎng)格;

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

5 模型驗證

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

5.1 水渠沖淤

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

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

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

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

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

圖4 不同時刻網(wǎng)格更新情況Fig.4 Mesh deformation in different time

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

5.2 海底管道局部沖刷

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

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

圖6 計算區(qū)域Fig.6 Computational domain

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

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

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

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

6 結(jié)論

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

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

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

［2］彭靜，河原能久.丁壩群近體流動結(jié)構(gòu)的可視化實驗研究［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，F(xiàn)REDS?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］韋雁機，葉銀燦.床面上短圓柱體局部沖刷三維數(shù)值模擬［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］李昌良.泥沙運動與底床變形的數(shù)值模擬［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.