張芝永，拾兵

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

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

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

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

1 水動(dòng)力模型

1.1 控制方程

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

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

1.2 湍流模型

對(duì)于不可壓縮流體且不考慮用戶自定義的源項(xiàng)時(shí)，標(biāo)準(zhǔn)k-ε模型方程為k方程:

ε方程:

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

2 泥沙輸運(yùn)模型

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

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

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

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

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

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

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

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

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

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

2.2 床面變形模型

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

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

3 動(dòng)邊界處理方法

3.1 動(dòng)邊界處理方法

動(dòng)網(wǎng)格模型可以用來模擬流場(chǎng)形狀由于邊界運(yùn)動(dòng)而隨時(shí)間改變的問題.其廣泛應(yīng)用于活塞、閥門及柔性體等有運(yùn)動(dòng)邊界工況的模擬.床面的沖淤變形形式類似于柔性體變形，需要對(duì)各個(gè)節(jié)點(diǎn)的運(yùn)動(dòng)情況進(jìn)行描述.因此本文采用FLUENT的DEFINE_GRID_MOTION宏命令來給出邊界節(jié)點(diǎn)的運(yùn)動(dòng)方式.在床面節(jié)點(diǎn)位置更新后，需要對(duì)區(qū)域內(nèi)部網(wǎng)格進(jìn)行調(diào)整.FLUENT提供了3種網(wǎng)格更新方式:1)光順網(wǎng)格法，2)動(dòng)態(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)格進(jìn)行重新劃分，各區(qū)域網(wǎng)格尺寸尺度基本保持不變，靠近下部變形邊界的網(wǎng)格仍然比較密集，這對(duì)于確保近壁面流場(chǎng)計(jì)算的準(zhǔn)確性至關(guān)重要.

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

3.2 沙滑模型

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

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

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

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

式中:p代表床面節(jié)點(diǎn)序號(hào);n為當(dāng)前時(shí)間步;n+1為下一時(shí)間步;Δtb為床面更新時(shí)間步長.

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

1)計(jì)算多個(gè)步長的速度、壓力、湍流數(shù)及切應(yīng)力等;

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

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

4)返回步驟1)，重新上面的計(jì)算，直至到達(dá)指定時(shí)間.

5 模型驗(yàn)證

為驗(yàn)證本文所建模型的可靠性，通過2個(gè)算例對(duì)其進(jìn)行了驗(yàn)證.

5.1 水渠沖淤

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

圖2 試驗(yàn)布置Fig.2 Layout of experiment

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

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

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

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

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

5.2 海底管道局部沖刷

本算例對(duì)近床面的海底管道局部沖刷的沖刷歷程進(jìn)行詳細(xì)的分析.

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

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

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

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

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

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

6 結(jié)論

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

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

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

［2］彭靜，河原能久.丁壩群近體流動(dòng)結(jié)構(gòu)的可視化實(shí)驗(yàn)研究［J］.水利學(xué)報(bào)，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］彭靜，河源能久.線性與非線性紊流模型及其在丁壩繞流中的應(yīng)用［J］.水動(dòng)力學(xué)研究與進(jìn)展: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ī)，葉銀燦.床面上短圓柱體局部沖刷三維數(shù)值模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展: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］李昌良.泥沙運(yùn)動(dòng)與底床變形的數(shù)值模擬［D］.青島:中國海洋大學(xué)，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.