二級(jí)沙波水流特性的大渦模擬

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

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

二級(jí)沙波水流特性的大渦模擬

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

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

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

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

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

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

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

1.1基本方程

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

1.2邊界條件

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

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

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

1.3網(wǎng)格生成

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

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

1.4計(jì)算模型

本研究共建4個(gè)具有不同沙波形態(tài)的數(shù)值模型，其概念模型如圖2所示。一級(jí)沙波采用Mielo和Ruiter［11］試驗(yàn)時(shí)的沙波尺度；二級(jí)沙波疊加在一級(jí)沙波迎流面上，波長(zhǎng)0.5 m，波高0.05 m，這與Fernandez［2］試驗(yàn)組合相似。試驗(yàn)過(guò)程中，二級(jí)沙波處于圖2所示的3種不同水平位置（模型b、c和d），二級(jí)沙波波谷距一級(jí)沙波波峰的長(zhǎng)度lS分別為0.5 m、0.16 m和0 m，其中模型a、c和d與Fernandez［2］和Frias［10］試驗(yàn)組合一致，模型b主要針對(duì)二級(jí)沙波離一級(jí)沙波波峰較遠(yuǎn)時(shí)的水流特性進(jìn)行研究。試驗(yàn)過(guò)程中，平均流速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 計(jì)算結(jié)果與分析

2.1模型驗(yàn)證

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

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

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

2.2縱向流速

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

2.3垂向流速

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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