二級沙波水流特性的大渦模擬

董偉良，諸裕良，馬殿光，徐俊鋒

（1.河海大學(xué)港口海岸與近海工程學(xué)院，南京210098；2.交通運輸部天津水運工程科學(xué)研究所，天津300456）

二級沙波水流特性的大渦模擬

董偉良1，諸裕良1，馬殿光2，徐俊鋒2

（1.河海大學(xué)港口海岸與近海工程學(xué)院，南京210098；2.交通運輸部天津水運工程科學(xué)研究所，天津300456）

采用FLuent軟件平臺大渦模擬方法（LES）對二級沙波水流特性進行數(shù)值模擬，進出口邊界采用周期性條件以模擬一系列沙波。通過與水槽試驗數(shù)據(jù)對比發(fā)現(xiàn)，模擬計算結(jié)果與試驗數(shù)據(jù)吻合較好，在此基礎(chǔ)上對3種不同水平位置處的二級沙波進行模擬。計算結(jié)果表明：當(dāng)二級沙波距一級沙波波峰較遠時，二級沙波對一級沙波起遮掩作用，使得一級沙波回流區(qū)范圍減小，垂向負流速區(qū)更加接近一級沙波波峰，同時二級沙波背流面的水流會對一級沙波迎流面形成侵蝕，進而減小一級沙波尺度；當(dāng)一、二級沙波接近融合或融合時，兩者之間的非線性相互作用使得回流區(qū)范圍擴大，紊動增強。Q準(zhǔn)則等值面顯示，二次沙波的疊加增強了沙波渦旋結(jié)構(gòu)的尺度和強度。

二級沙波；大渦模擬；水平位置；Q準(zhǔn)則

在河流和海洋環(huán)境中，因底床形態(tài)不斷發(fā)展演變，沙波地形也處于不斷疊加、融合的過程中，有效地保存了河床底部沉積物的形態(tài)［1］。Fernandez［2］試驗結(jié)果表明：二級沙波沿一級沙波迎流面從波谷移向波峰時，波谷產(chǎn)生的非線性作用使得沙波水流流態(tài)處于不斷的演變之中。數(shù)值模擬方法是研究沙波水流條件及泥沙運動規(guī)律的主要手段之一，目前常采用k?ε紊流模型，如：Mendoza［3］、Johns［4］、Yoon［5］等人成功模擬了不同條件下沙波水流特性，但模擬結(jié)果與實測數(shù)據(jù)在紊動較強的波谷差距較大；為了克服k?ε紊流模型只能模擬時均值及偏差較大的缺點，Yue［6］、Stoesser［7］、Grigoriadis［8］、Xie［9］等選擇利用大渦數(shù)值模型（LES）來進行研究，相對而言，LES模型可以更為精細地刻畫沙波水流條件。

然而，目前沙波水流數(shù)值模擬研究主要針對一級沙波形態(tài)，較少對二級沙波影響下的沙波水流特性進行研究。Frias等［10］曾對Fernandez［2］的試驗過程進行模擬，主要對沙波接近融合和融合2種狀態(tài)進行研究，而對于二級沙波遠離一級沙波波峰的情況暫未涉及。本文擬采用LES模型對二級沙波上的水流特性進行模擬，通過與已有的物理試驗資料進行對比，分析水平位置對沙波水流流態(tài)的影響，以此研究沙波疊加、融合的演變機制。

1 數(shù)學(xué)模型

1.1基本方程

大渦模擬的思想是在流動區(qū)域內(nèi)對N?S方程進行網(wǎng)格過濾，從而得到較大尺度漩渦的基本方程組。

1.2邊界條件

當(dāng)水流流經(jīng)一系列沙波時，各個沙波上水流結(jié)構(gòu)存在周期性變化規(guī)律，單個沙波周期進出口水流特性基本相同。在保證計算準(zhǔn)確性的基礎(chǔ)上，為了減少模型網(wǎng)格數(shù)和計算時間，從而取為周期性邊界條件，以此可以保證在單位周期沙波長度的模擬情況下，沙波水流充分發(fā)展為湍流，其中周期性條件為

式中：Φ為縱向流速u、橫向流速v、垂線流速w和穩(wěn)定動能k等水力參數(shù)；λ為沙波波長。

由于近壁湍流的脈動能量主要集中在小尺度渦，從而給粗網(wǎng)格的大渦模擬帶來困難，為使Smagorinsky?Lilly亞網(wǎng)格模型很好地適應(yīng)近壁面的湍流邊界層，采用壁面函數(shù)法處理近壁湍流，即用半經(jīng)驗公式將自由流中的紊流與壁面附件的流動鏈接起來。在固壁上應(yīng)用無滑移條件，計算域的上邊界采用對稱邊界條件：?/?z=0和w=0。

1.3網(wǎng)格生成

計算域由二級沙波疊加在一級沙波迎流面上組成。使用FLUENT軟件前處理程序Gambit生成計算區(qū)域幾何體，再進行網(wǎng)格劃分，得到如圖1所示的四面體網(wǎng)格單元。模型計算域長為1.6 m，寬為1.0 m，高為0.372 m，采用結(jié)構(gòu)網(wǎng)格進行劃分，其中沿x、y和z軸3個方向分別為160、50和70個網(wǎng)格（如圖1所示，其中網(wǎng)格每隔3個顯示1個），相鄰網(wǎng)格比不超過1.05。

圖1 模型計算范圍及網(wǎng)格剖分Fig.1Calculation scope and computational grid

1.4計算模型

本研究共建4個具有不同沙波形態(tài)的數(shù)值模型，其概念模型如圖2所示。一級沙波采用Mielo和Ruiter［11］試驗時的沙波尺度；二級沙波疊加在一級沙波迎流面上，波長0.5 m，波高0.05 m，這與Fernandez［2］試驗組合相似。試驗過程中，二級沙波處于圖2所示的3種不同水平位置（模型b、c和d），二級沙波波谷距一級沙波波峰的長度lS分別為0.5 m、0.16 m和0 m，其中模型a、c和d與Fernandez［2］和Frias［10］試驗組合一致，模型b主要針對二級沙波離一級沙波波峰較遠時的水流特性進行研究。試驗過程中，平均流速uˉ為0.52 m/s，雷諾數(shù)Re為171 000，弗汝德數(shù)Fr為0.29。

圖2 4種床面形態(tài)示意圖Fig.2Four morphology categories of bed surface

2 計算結(jié)果與分析

2.1模型驗證

為了驗證模型的準(zhǔn)確性，選取Mielo和Ruiter［11］的試驗數(shù)據(jù)進行驗證。沙波波長1.6 m，波高0.08 m，如圖1所示，沙波表面泥沙中值粒徑D50=1.6 mm。本文選取T6組次試驗數(shù)據(jù)來對模型進行驗證，其中流量Q=0.257 m3/s，水深d=0.334 m，弗汝德數(shù)Fr=0.29。

圖3顯示了不同位置處縱向流速垂線分布情況，其中各測點距前一個波峰分別為0.13 m、0.48 m、0.82 m和1.27 m。水流從波峰流向波谷，再流向波峰，縱向流速剖面因背流面回流區(qū)影響垂線分布不均勻，后因迎流面坡度作用水流加速而產(chǎn)生次生流，使得縱向流速垂線分布又趨于均勻［12］?？傮w而言，LES計算結(jié)果和實測值吻合較好，只是在位置0.48 m和0.82 m處，底部流速稍大于實測值。

圖3 縱向平均流速垂線分布計算結(jié)果與實測值對比Fig.3Comparison of calculation results with measured values of mean longitudinal velocity distribution along vertical

2.2縱向流速

圖4為縱向平均流速等值線分布圖，流速隨過流斷面減小而增大，最大流速值出現(xiàn)在二級沙波波峰處。當(dāng)二級沙波從波谷移向波峰，u=0.6 m/s等值面不斷擴大，甚至出現(xiàn)u=0.7 m/s等值面。值得注意的是：當(dāng)水流流經(jīng)沙波時，因沙波地形影響，波谷縱向流速u＜0 m/s。以u=0 m/s作為u＜0 m/s的界限，由此可以觀察到，當(dāng)二級沙波距一級沙波波峰較遠時（模型b），縱向平均流速u＜0 m/s范圍小于模型a；而當(dāng)一、二級沙波接近融合和融合時（模型c、d），隨著二級沙波越接近一級波峰u＜0 m/s范圍越大。分析認為：當(dāng)二級沙波達到一定尺寸后，在其背流面會形成較為明顯的回流區(qū)，此時因二級沙波距一級波峰較遠，一、二級沙波回流沒有產(chǎn)生相互作用，二級沙波對一級沙波起到遮掩作用，減小了一級沙波回流區(qū)的強度和范圍。模型b與Schwammle等［13］提出新月形沙波孤立波現(xiàn)象時的水流條件相似，當(dāng)二級沙波尺度較大時，其背流面的水流會對一級沙波形成侵蝕，同時二級沙波阻礙了一級沙波的來沙，致使一級沙波尺寸不斷減少，二級沙波尺寸不斷增加，從而形成孤立波現(xiàn)象。

2.3垂向流速

圖4 縱向平均流速等值線分布圖Fig.4Distribution of mean longitudinal velocity isoline

圖5給出了4種模型的垂向流速等值線分布圖，由圖5可見，盡管二級沙波的疊加改變了一級沙波迎流面上的分布情況，但一般分布規(guī)律并沒有改變：因沙波的迎、背流面地形影響，在迎流面會出現(xiàn)垂向流速w＞0 m/s區(qū)域（標(biāo)記為U），背流面垂向流速w＜0 m/s區(qū)域（標(biāo)記為D）。垂向負流速區(qū)D因u＞0 m/s、w＜0 m/s，水流沖擊床面，泥沙運動較為劇烈，床面侵蝕也剖為嚴(yán)重。

圖5 垂向平均流速等值線分布圖Fig.5Distribution of mean vertical velocity isoline

因二級沙波位置影響，一級沙波背流面垂向負流速區(qū)D發(fā)生明顯改變。模型b中，一、二級沙波的垂向負流速區(qū)D沒有融合，其中以w=0 m/s和0.04 m/s等值線作為參考，從圖5中可以看出一級沙波背流面后的垂向負流速區(qū)D相對于模型a明顯減小。研究認為，垂向負流速區(qū)D的減小主要由以下2個因素引起：（1）二級沙波對一級沙波起遮掩作用；（2）二級沙波疊加在一級沙波迎流面上時，二級沙波上形成的垂向正流速區(qū)U在一定程度上限制了垂向負流速區(qū)D的發(fā)展。與Fernandez［2］試驗的試驗結(jié)果相同，在二級沙波愈加靠近一級沙波波峰的情形下，一、二級沙波背流面垂向負流速區(qū)D相融合，同時由圖5可以看出，D范圍也愈大，垂向流速w=-0.06 m/s等值面也離一級沙波波峰愈遠，甚至還出現(xiàn)了w=-0.08 m/s的垂向流速。

因一級沙波波陡λ=1/20小于二級沙波波陡λ=1/10，將垂向流速w=0.02 m/s作為垂向正流速區(qū)U的界限，可以看出：模型a的垂向正流速區(qū)U相對比較扁平，模型b因空間有限，受上游垂向負流速區(qū)D的擠壓，垂向正流速區(qū)U最小；隨著二級沙波愈加靠近一級沙波波峰，和垂向負流速區(qū)D一樣，垂向正流速區(qū)U范圍也逐漸增加。

2.4渦旋結(jié)構(gòu)

剪切層Kelvin?Helmholtz不穩(wěn)定性是沙波產(chǎn)生渦旋的主要原因［7、9］，圖6為某一時刻4種不同床面形態(tài)下橫向位置y=0.5 m截面上橫向渦量ωy的分布圖，其中在一、二級沙波背流面均形成了較明顯的大尺度渦旋。從圖6中可以看出，模型a中橫向渦量ωy主要貼附底床，上層水體中ωy值普遍偏??；模型b中，一、二級沙波產(chǎn)生的橫向渦量ωy并沒有相互作用，在二級沙波背流面可以明顯觀測到二級沙波產(chǎn)生的橫向渦量ωy，較強的渦旋會對床面造成侵蝕；模型c、d中，一、二級沙波產(chǎn)生的橫向渦量ωy相互作用，其中模型d中，在上層水體中也可以觀察到較大的橫向渦量ωy。

圖6 某一時刻橫向渦量ωy分布圖Fig.6Distribution of transverse vorticityωyat a specific time

為了更清楚地顯示流場中的渦旋特征及破碎過程，采用Q準(zhǔn)則［14］方法進行識別。圖7顯示了Q=9對應(yīng)的渦旋結(jié)構(gòu)。從圖7中可知，渦旋結(jié)構(gòu)主要集中在波谷附近出現(xiàn)，根據(jù)前人研究分類，沙波水體中渦旋結(jié)構(gòu)一般可分為：橫向渦、發(fā)卡（馬蹄）渦和縱向渦3類［8-9］。由于剪切層Kelvin?Helmholtz不穩(wěn)定性，在一、二級沙波波峰處首先會出現(xiàn)橫向渦，隨著向下游發(fā)展，橫向渦逐漸演變成發(fā)卡渦。隨著流場進一步向下游發(fā)展，發(fā)卡渦逐漸向斜上方演變，由于上層水體流速大，下層流速小，渦旋結(jié)構(gòu)被不斷拉長扭曲，其中一部分到達水面形成水面翻滾現(xiàn)象，一部分演變成縱向渦進而破碎成隨機分布的小蝸，并逐漸消失。

當(dāng)只有一級沙波時（模型a），與橫向渦量ωy顯示結(jié)果相同，渦旋結(jié)構(gòu)主要集中在波峰線以下水體。當(dāng)二級沙波疊加時，背流面渦旋結(jié)構(gòu)尺度變大，其中在二級波峰頂處可以明顯觀察到渦旋結(jié)構(gòu)B，這主要是由于二級沙波產(chǎn)生的渦旋結(jié)構(gòu)在上部相對較為平靜的水體中，隨著流場的發(fā)展，能夠較好地保持結(jié)構(gòu)形態(tài)而不破碎。圖7顯示，當(dāng)二級沙波離一級沙波波峰較遠時（模型b）與接近融合和融合（模型c、d）狀態(tài)不同，二級沙波產(chǎn)生的渦旋結(jié)構(gòu)沒有能夠繼續(xù)發(fā)展，部分已經(jīng)破碎，同時在上層水體中也沒出現(xiàn)模型c、d中大尺度的渦旋結(jié)構(gòu)K。

圖7 某一時刻渦結(jié)構(gòu)分布圖（Q=9）Fig.7Distribution of vortex structure at a specific time（Q=9）

3 結(jié)論

本文利用LES模型對不同水平位置處的二級沙波進行數(shù)值模擬。比較本文的計算結(jié)果和實測數(shù)據(jù)，可以看出LES可以準(zhǔn)確刻畫沙波上的水流特性。二級沙波從波谷移向波峰時，其對一級沙波的作用經(jīng)歷了從遮掩侵蝕到非線性相互作用的演變：當(dāng)二級沙波離一級波峰較遠時，二級沙波形態(tài)對下游一級沙波起遮掩作用，使得一級沙波回流區(qū)面積縮小，垂向負流速區(qū)域更加靠近一級波峰，同時二級沙波背流面的水流會對一級沙波迎流面形成侵蝕，甚至?xí)p小一級沙波尺度；當(dāng)一、二級沙波接近融合或融合時，兩者非線性作用使得背流面回流區(qū)面積擴大，紊動增強。Q準(zhǔn)則等值面顯示，二次沙波的疊加使得沙波水體中渦旋結(jié)構(gòu)尺度和強度均勻增加，因上層水體紊動較小，二級沙波產(chǎn)生的渦旋結(jié)構(gòu)更易保持形態(tài)而不破碎。

［1］Best J，Blois G，Barros J，et al.The dynamics of bedform amalgamation：new insights from a very thin flume［C］//Fourth Internation?al Conference on Marine and River Dune Dynamics.Belgium：Bruges，2013：29-34.

［2］Fernandez R，Best J，López F.Mean flow，turbulence structure，and bed form superimposition across the ripple?dune transition［J］. Water Resources Research，2006，42（5）：72-88.

［3］Mendoza C，Wen S H.Investigation of turbulent flow over dunes［J］.Journal of Hydraulic Engineering，1990，116（4）：459-477.

［4］Johns B，Soulsby R L，Xing J X.A comparison of numerical model experiments of free surface flow over topography with flume and field observations［J］.Journal of Hydraulic Research，1993，31（2）：215-228.

［5］Yoon J Y，Patel V C.Numerical Model of Turbulent Flow over Sand Dune［J］.Journal of Hydraulic Engineering，1996，122（1）：10-18.

［6］Yue W，Lin C L，Patel V C.Numerical investigations of turbulent free surface flows using level set method and large eddy simula?tion［R］.Iowa City：Iowa Inst.of Hydraul.Res，2003.

［7］Stoesser T，Braun C，Garcia-Villalba M，et al.Turbulence structures in flow over two?dimensional dunes［J］.Journal of Hydraulic Engineering，2008，134（1）：42-55.

［8］Grigoriadis D G E，Balaras E，Dimas A A.Large?eddy simulations of unidirectional water flow over dunes［J］.Journal of Geophysical Research：Earth Surface，2009，114（F2）：91-100.

［9］Xie Z，Lin B，F(xiàn)alconer R A.Turbulence characteristics in free?surface flow over two?dimensional dunes［J］.Journal of Hydro?envi?ronment Research，2014（8）：200-209.

［10］Frias C E，Abad J D.Mean and turbulent flow structure during the amalgamation process in fluvial bed forms［J］.Water Resources Research，2013，49（10）：6 548-6 560.

［11］Mierlo M，De Ruiter J C C.Turbulence measurements above artificial dunes［R］.Delft：Delft Hydraulics Lab，1988.

［12］馬殿光，董偉良，徐俊鋒.沙波迎流面流速分布公式［J］.水科學(xué)進展，2015，26（3）：396-403. MA D G，DONG W L，XUN J F.Velocity distribution of nonuniform flow on the stoss side over dunes［J］.Advances in Water Sci?ence，2015，26（3）：396-403.

［13］Schwammle V，Hermann H J.Solitary wave behaviour of sand dunes［J］.Nature，2003，426（11）：619-620.

［14］Jeong J，Hussain F.On the identification of a vortex［J］.Journal of Fluid Mechanics，1995，285（4）：69-94.

Large?eddy simulation of flow characteristics over secondary dunes

DONG Wei?liang1,ZHU Yu?liang1,MA Dian?guang2,XU Jun?feng2
（1.College of Harbor,Coastal and Offshore Engineering，Hohai University，Nanjing 210098，China;2.Tianjin Research Institute for Water Transport Engineering,M.O.T.,Tianjin 300456,China）

Based on the Fluent software,the large eddy simulation(LES)technique was used to model the flow characteristics over secondary dunes.Periodic boundary condition was imposed at the inlet and outlet to simulate a series of dunes.Good agreement was obtained between the simulation and experimental data,and secondary dunes at three different positions were simulated.The results show that the secondary dune has a shelter effect on the pri?mary dune when far from the crest of the primary dune,which reduces the range of recirculation zone and makes ver?tical negative velocity zone more close to the primary dune crest.Additionally,the current at the lee side of second?ary dunes can cause erosion at the upstream side of the primary dune,thus diminish the dimension of the primary dunes.However,recirculation zone becomes large and turbulence intensities increase due to nonlinear interaction between current and dunes when secondary dunes approach and merge with the primary dune.The size and strength of vortex structures,visualized by the isosurface of Q?criterion,were strengthened by the secondary dune.

secondary dune;large eddy simulation;horizontal position;Q?criterion

TV 135；O 35

A

1005-8443（2016）02-0154-05

2015-07-29；

2015-11-30

中央級公益性科研院所基本科研業(yè)務(wù)費專項資金項目（TKS150102，TKS130105）

董偉良(1990-)，男，江蘇省人，碩士研究生，主要從事河流水力學(xué)、港口航道研究。

Biography：DONG Wei?liang（1990-），male，master student.