不同湍流模式下錢(qián)塘江涌潮水流三維模擬

汪求順，潘存鴻

(1. 浙江省水利河口研究院，浙江 杭州 310020; 2. 浙江省海洋規(guī)劃設(shè)計(jì)研究院，浙江 杭州 310020; 3. 浙江省河口海岸重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310016)

不同湍流模式下錢(qián)塘江涌潮水流三維模擬

汪求順1, 2, 3，潘存鴻1, 2, 3

(1. 浙江省水利河口研究院，浙江 杭州 310020; 2. 浙江省海洋規(guī)劃設(shè)計(jì)研究院，浙江 杭州 310020; 3. 浙江省河口海岸重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 310016)

錢(qián)塘江涌潮具有動(dòng)力作用強(qiáng)和流速變化快等特點(diǎn)。涌潮水流紊動(dòng)復(fù)雜，流速的垂向分布和紊動(dòng)強(qiáng)度息息相關(guān)。通過(guò)涌潮水流實(shí)測(cè)資料的分析可以發(fā)現(xiàn)，涌潮作用下流速垂向分布在底部和上層存在差異。為研究涌潮作用下流速垂向分布的特征，應(yīng)用基于非結(jié)構(gòu)網(wǎng)格下有限體積法模型FVCOM對(duì)錢(qián)塘江涌潮河段水流運(yùn)動(dòng)進(jìn)行三維數(shù)值模擬?？紤]到涌潮紊動(dòng)作用復(fù)雜且對(duì)流速的垂向分布起著重要影響，采用不同的湍流模式對(duì)涌潮傳播過(guò)程中水流的運(yùn)動(dòng)特征開(kāi)展研究。通過(guò)與涌潮河段實(shí)測(cè)資料的驗(yàn)證，復(fù)演涌潮到達(dá)前后水流運(yùn)動(dòng)特征，給出涌潮水流湍動(dòng)能的變化過(guò)程。研究成果有助于深入認(rèn)識(shí)涌潮水流紊動(dòng)特征和流速的分布規(guī)律，為涌潮作用下物質(zhì)輸運(yùn)的研究提供基礎(chǔ)。

涌潮；三維模擬；湍動(dòng)能；錢(qián)塘江河口；FVCOM

Abstract: The tidal bore in the Qiantang estuary has the characteristics of intense hydrodynamics and rapid-variation velocity. The turbulence of tidal bores is complicated, and turbulence intensity is closely relevant to the vertical distribution of velocity. It is found that the distribution of the velocity along the water depth is different at the bottom and the top layers under the tidal bore by the analysis of the field data. This study is to simulate the three-dimensional feature of the tidal bore in the Qiantang estuary based on the unstructured grids finite-volume method FVCOM model so as to investigate the distribution of velocity along the water depth. As the turbulence of the tidal bore is complicated and it plays an important role in the distribution of velocity, the tidal flows are studied during the propagation of tidal bores by different turbulence models. The movement of the flow is reproduced before and after the arrival of tidal bores by the calibration with the measured data, the time-varying process of turbulent kinetic energy is exhibited. The present results are helpful to deepen the understanding of turbulence characteristics and the distribution of velocity, and it will provide a basis for the study of mass transport.

Keywords: tidal bore; 3D simulation; turbulent kinetic energy; Qiantang estuary; FVCOM

目前，對(duì)于涌潮作用下二維水流運(yùn)動(dòng)的數(shù)值研究成果較多。潘存鴻等[1]通過(guò)守恒型淺水方程采用間斷捕捉方法對(duì)錢(qián)塘江河口涌潮的平面運(yùn)動(dòng)特征進(jìn)行了模擬。在涌潮水流現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)的分析中可以發(fā)現(xiàn)，涌潮到達(dá)前后流速不僅在平面上存在突變，而且水深方向的流速在底部和上層存在差異。涌潮水流具有明顯的時(shí)空梯度分布特征。為完整地反映涌潮水流的結(jié)構(gòu)，需通過(guò)三維模型來(lái)復(fù)演水流的運(yùn)動(dòng)特征。李紹武等[2]通過(guò)建立三維潮流數(shù)學(xué)模型研究概化河道中涌潮水流運(yùn)動(dòng)。王燦星等[3]采用FLUENT軟件對(duì)涌潮影響下的錢(qián)塘江三維水流結(jié)構(gòu)進(jìn)行了模擬。Liu等[4]基于無(wú)網(wǎng)格光滑粒子方法對(duì)涌潮的三維流場(chǎng)進(jìn)行了模擬。謝東風(fēng)等[5]基于FVCOM模型對(duì)錢(qián)塘江河口涌潮的三維運(yùn)動(dòng)進(jìn)行了模擬，并指出底部糙率對(duì)流速垂向分布的準(zhǔn)確模擬起著主要作用。陳永平等[6]在潮流模擬過(guò)程中選取不同渦黏系數(shù)類(lèi)型研究三維水流的結(jié)構(gòu)，結(jié)果表明垂向渦黏系數(shù)對(duì)水平流速的垂向分布起著重要作用。可以得出，涌潮作用下水流垂向分布不僅和底部糙率有關(guān)，還與湍流垂向紊動(dòng)有著緊密聯(lián)系。已有的涌潮作用下水流平面運(yùn)動(dòng)的數(shù)值研究表明，錢(qián)塘江河床糙率偏小[1]。當(dāng)前對(duì)錢(qián)塘江河口涌潮作用下水流紊動(dòng)特征的研究成果很少，對(duì)其進(jìn)行深入研究有助于探討水流的垂向分布機(jī)理。

涌潮在傳播過(guò)程中會(huì)產(chǎn)生強(qiáng)烈的湍流紊動(dòng)，國(guó)外學(xué)者通過(guò)物理模型試驗(yàn)、大渦模擬技術(shù)對(duì)涌潮作用下湍流過(guò)程進(jìn)行了研究[7-9]。Xie和Pan[10]通過(guò)鹽官站點(diǎn)涌潮水流的實(shí)測(cè)數(shù)據(jù)，對(duì)涌潮作用下水流的湍流紊動(dòng)進(jìn)行了初步研究，給出了雷諾應(yīng)力的大小。湍流紊動(dòng)和水流垂向分布有著緊密的聯(lián)系。湍流紊動(dòng)會(huì)引起水體湍動(dòng)能的變化，進(jìn)而影響水流的結(jié)構(gòu)分布。熊偉等[11]在水動(dòng)力方程中結(jié)合不同湍流模型研究了淺灘海域風(fēng)暴作用下潮流的垂向結(jié)構(gòu)。對(duì)于錢(qián)塘江涌潮河段，在強(qiáng)涌潮作用下，湍流紊動(dòng)過(guò)程復(fù)雜，這對(duì)水體的對(duì)流擴(kuò)散起著重要影響。為說(shuō)明涌潮作用下水流在垂向上的分布差異，對(duì)涌潮水流紊動(dòng)中湍動(dòng)能的定量分析有助于認(rèn)識(shí)涌潮水流的流速分布機(jī)理。

現(xiàn)有的研究準(zhǔn)確地模擬了涌潮的平面運(yùn)動(dòng)特征，對(duì)錢(qián)塘江涌潮的形態(tài)進(jìn)行了很好地復(fù)演，但反映涌潮水流垂向分布特征的數(shù)值研究很少。為準(zhǔn)確地復(fù)演涌潮水流的垂向分布，應(yīng)用基于非結(jié)構(gòu)網(wǎng)格下有限體積法的開(kāi)源模型FVCOM，通過(guò)結(jié)合不同的湍流模式，對(duì)錢(qián)塘江涌潮達(dá)到前后流速的分布進(jìn)行研究，給出涌潮初始段和強(qiáng)涌潮段湍動(dòng)能的變化過(guò)程。

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

1.1控制方程

控制方程采用FVCOM模型中經(jīng)雷諾時(shí)均的三維σ坐標(biāo)下Navier-Stokes方程。為提高三維模型在實(shí)際河口區(qū)域的計(jì)算效率，考慮到垂向采用靜壓假定的三維模型在模擬涌潮時(shí)具有很好的精度[12]，因此垂向動(dòng)量方程采用靜壓假定。

式中：t為時(shí)間；x，y和σ分別為水平和垂向坐標(biāo)；ζ為水面的高度；u，v，w分別為x，y，σ方向流速；D為總水深；H為靜水深；f為科氏常數(shù)；g為重力加速度；Am為水平紊動(dòng)黏性系數(shù)；Km為垂向渦黏系數(shù)。考慮到水平紊動(dòng)系數(shù)對(duì)流速分布的影響小[6]，采用Smagorinsky方程進(jìn)行計(jì)算。湍流的垂向紊動(dòng)對(duì)流速沿水深分布影響大，將通過(guò)不同湍流模式進(jìn)行分析。

1.2湍流模式

1.2.1 M-Y模式

FVCOM模型中垂向紊動(dòng)的默認(rèn)湍流模式為Mellor-Yamada方程。該湍流模式的湍動(dòng)能和混合長(zhǎng)度的方程如下

式中：q2/2為湍動(dòng)能；l為湍流混合長(zhǎng)度；B1，E1，E2為模型閉合常數(shù)；Kq為湍流垂向擴(kuò)散系數(shù)，即Kq=0.2lq；垂向渦黏系數(shù)Km=lqSm，Sm為穩(wěn)定函數(shù)；壁面函數(shù)中的L=(ζ-z)-1+(H+z)-1；κ為卡門(mén)常數(shù)，即κ=0.4；Fq，F(xiàn)l分別為湍動(dòng)能和混合長(zhǎng)度的水平擴(kuò)散項(xiàng)。為減小水平擴(kuò)散影響，兩者取值均為0。

水面和底部的邊界條件：

式中：uτs,uτb分別為水面和底部的摩阻流速。

1.2.2 k-ε模式

k-ε湍流模式適合高雷諾數(shù)下水流紊動(dòng)的求解。為分析涌潮水流的強(qiáng)紊動(dòng)特征，采用該湍流模式進(jìn)行水體紊動(dòng)強(qiáng)度的分析。對(duì)于k-ε湍流模型可表示為

式中：k為湍流動(dòng)能；ε為湍流耗散率；垂向渦黏系數(shù)Km=cμk2/ε；cμ，c1，c2，k和ε分別為0.09、1.44、1.92、1.00和1.30。

水面條件：

底部的邊界條件：

1.3水流邊界

對(duì)于涌潮水流的表面，水平速度的垂向梯度為0。垂向速度按下式給出：

對(duì)于水流的底邊界，水平流速的分布和水深、底部切應(yīng)力以及垂向渦黏系數(shù)有關(guān)，按式(13)確定。垂向速度按式(14)計(jì)算。

式中：τbx,τby分別為x,y方向底部切應(yīng)力。

底部切應(yīng)力按如下關(guān)系進(jìn)行計(jì)算：

式中：Cd為底部摩阻系數(shù)；ρ0為水的密度；ub,vb分別為x,y方向床面的流速。

在FVCOM模型中，底部摩阻系數(shù)默認(rèn)采用式(16)。

式中：z0b為床面粗糙度，取為2.5d50，d50為床面泥沙的中值粒徑；z1為近底層網(wǎng)格中心到底部的距離。

根據(jù)已有錢(qián)塘江河口涌潮水流的數(shù)值研究結(jié)果，錢(qián)塘江河床糙率偏小。在涌潮平面運(yùn)動(dòng)的數(shù)值模擬中，一般采用較小糙率系數(shù)的曼寧公式即式(17)進(jìn)行底部摩阻系數(shù)的計(jì)算。

式中：n為曼寧糙率系數(shù)，反映床面粗糙程度。

1.4離散處理

考慮到有限體積法能保證計(jì)算區(qū)域物理量的守恒，模型控制方程采用有限體積法進(jìn)行數(shù)值離散。對(duì)流項(xiàng)的空間離散采用二階精度的計(jì)算格式，時(shí)間步進(jìn)采用二階精度的龍格庫(kù)塔法，垂向速度采用隱式格式進(jìn)行求解。模型能很好地計(jì)算水躍等強(qiáng)間斷水動(dòng)力問(wèn)題[13]。該模型在模擬涌潮水流的突變等性質(zhì)方面具有一定優(yōu)勢(shì)，將利用該開(kāi)源模型研究不同湍流模式下錢(qián)塘江河口水流的三維運(yùn)動(dòng)特征。

2 模型驗(yàn)證

2.1模型計(jì)算范圍

為模擬涌潮從生成到衰減過(guò)程的水流變化特征，模型計(jì)算區(qū)域從涌潮生成段的澉浦附近開(kāi)始到衰減段的閘口。模型計(jì)算范圍見(jiàn)圖1，其中東西向長(zhǎng)為89 120 m，南北向?qū)挒?8 840 m。涌潮傳播速度和網(wǎng)格劃分的分辨率有關(guān)，高分辨率網(wǎng)格能細(xì)化局部地形[5]。另外，受河口平面形態(tài)和沙坎影響，涌潮在傳播到鹽官時(shí)強(qiáng)度達(dá)到最大。為俘獲涌潮傳播過(guò)程中的最大流速，水平網(wǎng)格的分辨率從鹽官區(qū)域的50 m向上下游邊界逐漸增大到100～400 m。模型計(jì)算區(qū)域網(wǎng)格劃分的單元為45 010個(gè)，節(jié)點(diǎn)為23 306個(gè)。

圖1 模型計(jì)算范圍和觀測(cè)站點(diǎn)的位置Fig. 1 Computational domain and locations of measured stations

2007年10月25日12:00～30日12:00在錢(qián)塘江涌潮河段進(jìn)行了12個(gè)站點(diǎn)潮流數(shù)據(jù)的測(cè)量，并利用潮位觀測(cè)站點(diǎn)進(jìn)行同步潮位的測(cè)量，潮位和潮流觀測(cè)站點(diǎn)布置如圖1所示。在無(wú)涌潮時(shí)每小時(shí)記錄一次潮位數(shù)據(jù)，涌潮到達(dá)后每1～2分鐘記錄一次數(shù)據(jù)。對(duì)于潮流的觀測(cè)，每個(gè)站位在水深大于4 m時(shí)使用垂向6點(diǎn)法進(jìn)行測(cè)量，即分別位于表層(距離水面0.5 m)、0.2倍水深、0.4倍水深、0.6倍水深、0.8倍水深和底層(距離床面0.5 m)。無(wú)涌潮時(shí)，每小時(shí)記錄一次流速數(shù)據(jù)，涌潮到達(dá)前后每分鐘記錄一次流速數(shù)據(jù)。

上、下游開(kāi)邊界根據(jù)實(shí)測(cè)潮位給定。模型計(jì)算開(kāi)始時(shí)刻，區(qū)域內(nèi)的潮位和流速均設(shè)為0。垂向分12層。采用模型中已有的內(nèi)外模分裂算法，外模計(jì)算時(shí)間步長(zhǎng)為0.1 s，內(nèi)模為1.0 s。模型計(jì)算時(shí)段從2007年10月25日00:00開(kāi)始到10月31日00:00結(jié)束，模型驗(yàn)證從計(jì)算后12小時(shí)即10月25日12:00開(kāi)始以消除初始啟動(dòng)影響。考慮到涌潮水流模擬計(jì)算的耗時(shí)性，采用刀片服務(wù)器進(jìn)行并行計(jì)算，對(duì)錢(qián)塘江從澉浦附近到閘口的涌潮河段進(jìn)行三維水流模擬。

2.2數(shù)值驗(yàn)證

2.2.1 強(qiáng)涌潮生成

為復(fù)演鹽官站強(qiáng)涌潮到達(dá)時(shí)間和水位的猛增過(guò)程，利用不同底部摩阻系數(shù)類(lèi)型公式進(jìn)行模型的率定。圖2給出了不同類(lèi)型阻力系數(shù)計(jì)算公式(16)和(17)中底摩阻系數(shù)在模型區(qū)域的取值分布范圍。可以看出，在曼寧糙率系數(shù)n=0.009 5時(shí)，底摩阻系數(shù)分布在0.000 37～0.001 9之間。

考慮到模型的計(jì)算采用分裂算法，即三維模型通過(guò)二維模型得出的潮位ζ進(jìn)行求解。而在二維模型中，底部摩阻的大小直接影響涌潮高度和到達(dá)時(shí)間的模擬結(jié)果。圖3給出了利用不同底部切應(yīng)力公式進(jìn)行數(shù)值計(jì)算后鹽官站點(diǎn)10月28日涌潮到達(dá)時(shí)刻水位的比較。從圖3可以看出，在較大的底部摩阻系數(shù)下，鹽官站涌潮到達(dá)時(shí)間落后于較小底摩阻系數(shù)下的結(jié)果。同時(shí)，生成的涌潮高度低于較小底摩阻系數(shù)下的潮位值。在三維模型的控制方程中，底部切應(yīng)力作為二階項(xiàng)對(duì)流速垂向分布影響大，但對(duì)水位和平面流速的影響小。結(jié)合M-Y湍流模式對(duì)計(jì)算時(shí)段中模型區(qū)域測(cè)量站點(diǎn)的潮位和各層潮流進(jìn)行了驗(yàn)證。

圖2 不同類(lèi)型底部摩阻系數(shù)的分布Fig. 2 Distribution of different types of bottom friction coefficient

圖3 涌潮到達(dá)時(shí)刻水位比較Fig. 3 Comparison of water level at the arrival time of tidal bore

2.2.2 潮位和潮流

為驗(yàn)證模型在計(jì)算涌潮傳播過(guò)程中水流結(jié)果的可靠性，結(jié)合M-Y湍流模式對(duì)五天漲落潮過(guò)程中的全部測(cè)量站點(diǎn)的計(jì)算結(jié)果與實(shí)測(cè)資料進(jìn)行了對(duì)比。圖4為計(jì)算區(qū)域中強(qiáng)涌潮段鹽官和大缺口站點(diǎn)的模型計(jì)算潮位與實(shí)測(cè)數(shù)據(jù)的驗(yàn)證結(jié)果，其中曼寧系數(shù)取為0.009 5。圖中空心框點(diǎn)為實(shí)測(cè)潮位，實(shí)線(xiàn)為模型計(jì)算結(jié)果。

將模型計(jì)算的潮流流速和流向分別與觀測(cè)數(shù)據(jù)中的表層、中層(0.6倍水深)和底層進(jìn)行驗(yàn)證。圖5和圖6給出了曼寧系數(shù)取為0.025時(shí)強(qiáng)涌潮站點(diǎn)703和涌潮初始段站點(diǎn)712模型計(jì)算的流速和流向與實(shí)測(cè)數(shù)據(jù)的對(duì)比，其中空心圓點(diǎn)為實(shí)測(cè)流速，空心方點(diǎn)為實(shí)測(cè)流向，實(shí)線(xiàn)為模型計(jì)算結(jié)果。

圖4 站點(diǎn)的潮位驗(yàn)證Fig. 4 Validation of tidal level at stations

圖5 703站點(diǎn)的表、中和底層潮流驗(yàn)證Fig. 5 Validation of tidal current at station 703

圖6 712站點(diǎn)的表、中和底層潮流驗(yàn)證Fig. 6 Validation of tidal current at station 712

從站點(diǎn)的潮位、潮流流速和流向的驗(yàn)證結(jié)果可以看出，各站點(diǎn)模型計(jì)算值和實(shí)測(cè)值的平均絕對(duì)誤差較小。703站點(diǎn)表、中、底層流速的相對(duì)誤差分別為13%、29%、15%；流向的相對(duì)誤差分別為12%、18%、14%。712站點(diǎn)各層流速的相對(duì)誤差分別為8%、6%和13%；流向的相對(duì)誤差分別為13%、14%和16%。模型計(jì)算的表、中、底層流速和流向值與實(shí)測(cè)值基本一致。相對(duì)于先前學(xué)者的研究結(jié)果[5]，本文采用的三維模型能夠更好地復(fù)演流速的垂向分布特征，較為準(zhǔn)確地模擬錢(qián)塘江涌潮河段潮流各層變化。

3 數(shù)值結(jié)果分析

3.1不同湍流模式下涌潮水流的垂向分布

3.1.1 沿程分布

為說(shuō)明不同湍流模式下水流運(yùn)動(dòng)的三維變化特征，提取了圖1中模型計(jì)算區(qū)域鹽官縱剖面1涌潮水流流速。圖7給出了M-Y和k-ε湍流模式下鹽官縱剖面在涌潮到達(dá)時(shí)流速的垂向分布(圖中左側(cè)為錢(qián)塘江河口的上游，右側(cè)為下游)。從圖7(a)可以看出，在涌潮到達(dá)前，上游河道的水位為2.1 m，水流以較小的速度向下游移動(dòng)。在漲潮過(guò)程中隨著下游水流向上游區(qū)域涌動(dòng)，下游水流到達(dá)鹽官段形成強(qiáng)涌潮，水位猛增到4.0 m。涌潮水流的上層流速達(dá)到3.5 m/s，并呈現(xiàn)由上層到底層減小的分布特征。從斷面的沿程流速分布可以看出，涌潮形成的后方水流以更大的流速向上游運(yùn)動(dòng)，其后上部水體的最大流速可達(dá)到5.0 m/s，并在水深方向呈現(xiàn)梯度分布。從圖7(b)可以看出，k-ε湍流模式在涌潮到達(dá)時(shí)水位猛增高度比M-Y湍流模式的計(jì)算結(jié)果偏小0.2 m。另外，k-ε湍流模式在計(jì)算涌潮到達(dá)同一點(diǎn)時(shí)比M-Y湍流模式慢1 000 m。同時(shí)，涌潮表面的水流速度低于M-Y湍流模式得出的結(jié)果。在強(qiáng)涌潮生成時(shí)，潮頭后方存在流速較大的水流。由涌潮觀測(cè)可知，涌潮過(guò)后，漲潮流速仍迅速增加，此類(lèi)為快水現(xiàn)象。從流速的沿程分布可以看出，本文采用的模型能復(fù)演在涌潮作用下河段區(qū)域的快水。

圖7 涌潮到達(dá)時(shí)刻鹽官河段流速垂向分布Fig. 7 Profile velocity of Yanguan along the water depth at the arrival of tidal bores

3.1.2 不同時(shí)刻涌潮水流的垂向分布

圖8給出了兩種湍流模式下涌潮到達(dá)前后鹽官站在不同時(shí)刻流速沿水深的分布，其中圖8(a)為M-Y湍流模式計(jì)算結(jié)果，圖8(b)為k-ε湍流模式計(jì)算結(jié)果。從圖中可以看出，在13:30-13:50時(shí)刻，兩種湍流模式得出的流速分布基本一致。在涌潮達(dá)到后的時(shí)刻14:00，M-Y湍流模式計(jì)算得出的潮位比k-ε湍流模式得出的結(jié)果偏高0.2 m。M-Y湍流模式計(jì)算得出的流速在14:00時(shí)刻比k-ε湍流模式得出的流速偏大。在1小時(shí)過(guò)后的最高潮位上，M-Y湍流模式計(jì)算得出的潮位比k-ε湍流模式得出的結(jié)果稍高0.1 m，流速的分布基本一致，沿水深分布存在梯度。在涌潮到達(dá)的一段時(shí)間內(nèi)，表層水流以5.0 m/s的速度沿水深呈遞減分布，并向上游運(yùn)動(dòng)。在15:30時(shí)刻以后，涌潮水流強(qiáng)度減弱，呈現(xiàn)一般漲潮流運(yùn)動(dòng)特征。潘存鴻等[1]指出在鹽官段大于5.0 m/s快水現(xiàn)象出現(xiàn)在涌潮后15 min，持續(xù)時(shí)間能達(dá)到33 min。從涌潮過(guò)后不同時(shí)刻斷面流速的分布可以看出，數(shù)學(xué)模型很好地模擬了強(qiáng)涌潮區(qū)域的快水，給出的鹽官?gòu)?qiáng)涌潮區(qū)域的快水在水深方向存在分布梯度。本文采用的三維模型結(jié)合不同的湍流模式能很好地反映出涌潮水流沿水深分布的運(yùn)動(dòng)特征。

圖8 不同時(shí)刻兩種湍流模式下鹽官站流速垂向分布Fig. 8 Distribution of time-series velocity at Yanguan under two types of turbulence models

圖9 不同湍流模式下潮位過(guò)程和實(shí)測(cè)值的對(duì)比Fig. 9 Comparison of tidal level under two types of turbulence models

3.1.3 潮位和流速過(guò)程與實(shí)測(cè)值的對(duì)比

圖9給出了兩種湍流模式下鹽官站潮位過(guò)程和實(shí)測(cè)值的對(duì)比。從圖中可以看出，兩種湍流模式均反映了涌潮到達(dá)時(shí)刻水位的猛增過(guò)程，但k-ε湍流模式得出涌潮到達(dá)時(shí)刻的結(jié)果比M-Y模式慢6 min。同時(shí)，在得到的水位高度上稍低于M-Y模式的結(jié)果。圖10給出了兩種湍流模式下鹽官站流速過(guò)程和實(shí)測(cè)值的對(duì)比。兩種湍流模式得出的流速變化基本一致，但在相位上k-ε湍流模式落后于M-Y湍流模式計(jì)算結(jié)果。這和涌潮作用下產(chǎn)生的湍流紊動(dòng)有著緊密的關(guān)系。

圖10 不同湍流模式下流速過(guò)程和實(shí)測(cè)值的對(duì)比Fig. 10 Comparison of velocity under two types of turbulence models

3.2湍動(dòng)能

為反映不同湍流模式中涌潮作用下水流運(yùn)動(dòng)產(chǎn)生差異的內(nèi)部機(jī)理，給出兩種湍流模式下在涌潮開(kāi)始生成和強(qiáng)涌潮處湍動(dòng)能的變化過(guò)程。提取了圖1中模型區(qū)域涌潮開(kāi)始生成處的代表站點(diǎn)A和強(qiáng)涌潮處的代表站點(diǎn)鹽官的湍動(dòng)能。圖11給出了涌潮開(kāi)始生成處的代表點(diǎn)A在兩種湍流模式下湍動(dòng)能的變化。圖12給出了強(qiáng)涌潮處鹽官站點(diǎn)在兩種湍流模式下湍動(dòng)能的變化。

圖11 涌潮初始生成處兩種湍流模式下湍動(dòng)能Fig. 11 Turbulent kinetic energy of two types of turbulence models at the initial formation period of tidal bores

從圖11可以看出，在涌潮初始生成處，湍動(dòng)能較小，呈現(xiàn)表層到底層增大的趨勢(shì)。M-Y湍流模式計(jì)算得出的湍動(dòng)能稍大于k-ε湍流模式的結(jié)果。從圖12可以看出，在強(qiáng)涌潮生成處，湍動(dòng)能增大到一個(gè)量級(jí)以上。M-Y湍流模式得出的湍動(dòng)能呈現(xiàn)由表、中、底層增大的分布特征。k-ε湍流模式計(jì)算得出的湍動(dòng)能表現(xiàn)為底層和中層較一致，表層小。兩種湍流模式給出的湍動(dòng)能在量級(jí)上差異不明顯。在強(qiáng)涌潮段，k-ε湍流模式得出的湍動(dòng)能較大，反映涌潮水流的脈動(dòng)作用較為強(qiáng)烈。這說(shuō)明了k-ε湍流模式在鹽官段計(jì)算得出的涌潮到達(dá)時(shí)間稍慢于M-Y湍流模式結(jié)果的原因。另外，涌潮的強(qiáng)紊動(dòng)特征使得k-ε湍流模式計(jì)算得出的涌潮高度稍低于M-Y湍流模式的計(jì)算結(jié)果。

圖12 強(qiáng)涌潮處兩種湍流模式下湍動(dòng)能Fig. 12 Turbulent kinetic energy of two types of turbulence models under the strong tidal bore

為定量地表示兩種湍流模式下湍動(dòng)能的大小，表1和表2分別給出了M-Y湍流模式和k-ε湍流模式計(jì)算得出的最小、最大和平均湍動(dòng)能。從兩表中可以看出，在涌潮初始生成處，M-Y湍流模式計(jì)算得出的最大湍動(dòng)能在中層和底層分別為0.011 7 m2/s2和0.022 2 m2/s2，大于k-ε湍流模式得出的結(jié)果。k-ε湍流模式在表層最大湍動(dòng)能為0.003 5 m2/s2，大于M-Y湍流模式得出的結(jié)果。在強(qiáng)涌潮生成處，k-ε湍流模式在表、中和底層的湍動(dòng)能均大于M-Y湍流模式計(jì)算值?？梢缘贸?，在涌潮水流的三維模擬中，k-ε湍流模式給出了較強(qiáng)的水流紊動(dòng)特征，使得在涌潮水位和流速上稍小于M-Y湍流模式的計(jì)算結(jié)果。

表1M-Y模式計(jì)算湍動(dòng)能

Tab.1TurbulentkineticenergybyturbulencemodelofMellor-Yamada(m2·s-2)

代表點(diǎn)表層中層底層最小最大平均值最小最大平均值最小最大平均值涌潮開(kāi)始處A點(diǎn)2．18′10-60．00240．00066．78′10-60．01170．00266．51′10-60．02220．0047強(qiáng)涌潮處鹽官站2．21′10-50．10760．01295．19′10-50．25110．03016．97′10-50．35490．0425

表2k-ε模式計(jì)算湍動(dòng)能(m2·s-2)

Tab.2Turbulentkineticenergybyturbulencemodelofk-ε(m2·s-2)

代表點(diǎn)表層中層底層最小最大平均值最小最大平均值最小最大平均值涌潮開(kāi)始處A點(diǎn)2．57′10-60．00350．00083．55′10-60．00720．00165．54′10-60．01300．0029強(qiáng)涌潮處鹽官站9．49′10-50．33440．04141．06′10-40．47220．05848．40′10-50．49810．0617

4 結(jié) 語(yǔ)

基于有限體積法FVCOM模型，并結(jié)合湍流模式對(duì)錢(qián)塘江涌潮水流的三維運(yùn)動(dòng)進(jìn)行了準(zhǔn)確的數(shù)值模擬。從三維數(shù)值模擬結(jié)果的分析可以看出，采用的模型能很好地模擬涌潮到達(dá)前后流速的垂向變化特征。通過(guò)兩種湍流模式下涌潮到達(dá)時(shí)間和涌潮高度等進(jìn)行的數(shù)值研究分析結(jié)果，可以發(fā)現(xiàn)，k-ε湍流模式計(jì)算得出的涌潮達(dá)到時(shí)間稍遲、潮位稍低，兩種模式計(jì)算的涌潮流速分布結(jié)果偏差較小，均能模擬涌潮水流三維運(yùn)動(dòng)特征。不同湍流模式給出的湍動(dòng)能在量級(jí)上基本一致，在強(qiáng)涌潮作用下水流湍動(dòng)能比涌潮初始生成段大一個(gè)量級(jí)以上。通過(guò)對(duì)錢(qián)塘江涌潮水流的三維模擬和湍動(dòng)能的分析，能為深入認(rèn)識(shí)涌潮水流運(yùn)動(dòng)提供技術(shù)手段，并為錢(qián)塘江河口物質(zhì)輸運(yùn)提供研究基礎(chǔ)。

[1] 潘存鴻, 林炳堯, 毛獻(xiàn)忠. 錢(qián)塘江涌潮二維數(shù)值模擬[J]. 海洋工程，2007, 25(1): 50-56. (PAN Cunhong, LIN Bingyao, MAO Xianzhong. 2D numerical simulation of tidal bore in Qiantang River[J]. The Ocean Engineering, 2007, 25(1), 50-56. (in Chinese))

[2] 李紹武, 盧麗鋒, 時(shí)鐘. 河口準(zhǔn)三維涌潮數(shù)學(xué)模型研究[J]. 水動(dòng)力學(xué)研究與進(jìn)展：A輯. 2004, 19(4): 407-415. (LI Shaowu, LU Lifeng, SHI Zhong. A quasi 3D numerical model for estuarine tidal bore[J]. Chinese Journal of Hydrodynamics, 2004, 19(4): 407-415. (in Chinese))

[3] 王燦星, 陳菊芳, 金晗輝, 等. 涌潮對(duì)錢(qián)塘江河道流場(chǎng)影響的三維數(shù)值模擬研究[J]. 水動(dòng)力學(xué)研究與進(jìn)展：A輯. 2012, 27(4): 367-375. (WANG Canxing, CHEN Jufang, JIN Hanhui, et al. Three-dimensional numerical simulation of tidal bore of Qiantang River[J].Chinese Journal of Hydrodynamics, 2012, 27(4): 367-375. (in Chinese))

[4] LIU H, LI J, SHAO S, et al. SPH modeling of tidal bore scenarios [J]. Natural Hazards, 2015, 75(2): 1247-1270.

[5] 謝東風(fēng), 潘存鴻, 吳修廣. 基于FVCOM模式錢(qián)塘江河口涌潮三維數(shù)值模擬研究[J]. 海洋工程, 2011, 29(1): 47-52. (XIE Dongfen, PAN Cunhong, WU Xiuguang. Three-dimensional numerical modeling of tidal bore in Qiantang based on FVCOM[J]. The Ocean Engineering, 2011, 29(1), 47-52. (in Chinese))

[6] 陳永平, 劉家駒, 喻國(guó)華. 潮流數(shù)值模擬中紊動(dòng)黏性系數(shù)的研究[J]. 河海大學(xué)學(xué)報(bào)：自然科學(xué)版, 2002, 30(1): 39-43. (CHEN Yongping, LIU Jiaju ,YU Guohua. A study on eddy viscosity coefficient in numerical tidal simulation[J]. Jounal of Hohai University, 2002, 30(1): 39-43. (in Chinese))

[7] CHANSON H, TAN K. Turbulent mixing of particles under tidal bores: An experimental analysis [J]. Journal of Hydraulic Research, 2010, 48(5): 641-649.

[8] DOCHERTY N J, CHANSON H. Physical modeling of unsteady turbulence in breaking tidal bores [J]. Journal of Hydraulic Engineering-ASCE, 2012, 138(5): 412-419.

[9] LU B P, CHANSON H, GLOCKNER S. Large eddy simulation of turbulence generated by a weak breaking tidal bore [J]. Environmental Fluid Mechanics, 2010, 10(5): 587-602.

[10] XIE Dongfeng, PAN Cunhong. A preliminary study of the turbulence features of the tidal bore in the Qiantang River, China [J]. Journal of Hydrodynamics, 2013, 25(6): 903-911.

[11] 熊偉, 朱志夏, 董佳, 等. 湍黏系數(shù)對(duì)淺灘海域三維風(fēng)暴潮的影響[J]. 水利水運(yùn)工程學(xué)報(bào), 2014(6): 45-51. (XIONG Wei, ZHU Zhixia, DONG Jia, et al. Effects of turbulent viscosity coefficient on 3-D storm surge within shallow seas[J] Hydro-Science and Engineering, 2014(6): 45-51. (in Chinese))

[12] GUSEV A V, LYAPIDEVSKII V Y. Turbulent bore in a supercritical flow over an irregular bed [J]. Fluid Dynamics, 2005, 40(1): 54-61.

[13] HUANG H, CHEN C S, COWLES G W, et al. FVCOM validation experiments: Comparisons with ROMS for three idealized barotropic test problems [J]. Journal of Geophysical Research-Oceans, 2008, 113(C07): 1-14.

Three-dimensional simulation of tidal bore in the Qiantang estuary under different turbulence models

WANG Qiushun1, 2, 3, PAN Cunhong1, 2, 3

(1. Zhejiang Institute of Hydraulics and Estuary, Hangzhou 310020, China; 2. Zhejiang Institute of Marine Planning and Design, Hangzhou 310020, China; 3. Key Laboratory of Estuarine and Coast of Zhejiang Province, Hangzhou 310016, China)

TV148

A

10.16483/j.issn.1005-9865.2017.01.009

1005-9865(2017)01-0080-10

2016-04-21

國(guó)家自然科學(xué)基金資助項(xiàng)目(51379190；41306082)

汪求順(1984-)，男，安徽桐城人，博士，主要從事河口海岸動(dòng)力學(xué)研究。E-mail:wangqiushun57@163.com