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

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

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

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

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

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

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

二級沙波；大渦模擬；水平位置；Q準則

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

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

1 數學模型

1.1基本方程

大渦模擬的思想是在流動區域內對N?S方程進行網格過濾，從而得到較大尺度漩渦的基本方程組。

1.2邊界條件

當水流流經一系列沙波時，各個沙波上水流結構存在周期性變化規律，單個沙波周期進出口水流特性基本相同。在保證計算準確性的基礎上，為了減少模型網格數和計算時間，從而取為周期性邊界條件，以此可以保證在單位周期沙波長度的模擬情況下，沙波水流充分發展為湍流，其中周期性條件為

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

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

1.3網格生成

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

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

1.4計算模型

本研究共建4個具有不同沙波形態的數值模型，其概念模型如圖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，雷諾數Re為171 000，弗汝德數Fr為0.29。

圖2 4種床面形態示意圖Fig.2Four morphology categories of bed surface

2 計算結果與分析

2.1模型驗證

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

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

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

2.2縱向流速

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

2.3垂向流速

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

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

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

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

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

2.4渦旋結構

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

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

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

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

圖7 某一時刻渦結構分布圖（Q=9）Fig.7Distribution of vortex structure at a specific time（Q=9）

3 結論

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

［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，Falconer 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］.水科學進展，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

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

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

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