謝東風(fēng),潘存鴻,吳修廣

(浙江省水利河口研究院,浙江 杭州 310020)

基于FVCOM模式錢塘江河口涌潮三維數(shù)值模擬研究

謝東風(fēng),潘存鴻,吳修廣

(浙江省水利河口研究院,浙江 杭州 310020)

應(yīng)用FVCOM模型進(jìn)行錢塘江涌潮的二維和三維數(shù)值模擬,并使用2007年10月大潮期間的現(xiàn)場(chǎng)觀測(cè)結(jié)果進(jìn)行模型驗(yàn)證。結(jié)果表明,平面二維模擬結(jié)果與實(shí)測(cè)值吻合較好,與基于Godunov算法求解錢塘江涌潮二維特征的結(jié)果基本一致。三維計(jì)算結(jié)果再現(xiàn)了流速在垂向上的差異,以及在落潮轉(zhuǎn)向漲潮時(shí)刻,漲潮先從底部開始,垂向上從底層向表層逐步從落潮流轉(zhuǎn)為漲潮流。敏感性分析顯示模型網(wǎng)格大小和曼寧系數(shù)的選取對(duì)模擬結(jié)果有較大的影響。

錢塘江;涌潮;數(shù)值模擬;FVCOM

Abstract:In this contribution,the FVCOM model has been used to simulate the tidal bore in the Qiantang River.The depth-averaged 2-D model results agreewellwith the field data,which is consistentwith the research of Pan et al.(2007),who used a depth-averaged 2-D numericalmodel based on a Godunov-type scheme.The 3-Dmodel reproduced the vertical distribution characteristicsof current velocitiesand the physical processes,and the phenomenon thatwhen the ebb currents turn to flood currents,velocitiesat bottom increase first.Numerical experiments revealed that the cell size andmanning coefficient could influence themodel results to a large extent.

Key words:Qiantang River;tidal bore;numericalmodeling;FVCOM

涌潮是一種特殊的自由面重力流動(dòng),水流穿越涌波,水位、流速急劇變化。強(qiáng)涌潮是世界罕見的自然景觀之一,世界上許多強(qiáng)潮河口都存在涌潮,如英國的塞文河、巴西的亞馬遜河、長江口北支等。錢塘江涌潮以兇猛、多變、驚險(xiǎn)而堪稱一絕,從古至今吸引不計(jì)其數(shù)的觀潮者。涌潮是一種罕見的自然景觀,其巨大能量對(duì)河口過程有重要的影響,也影響到航運(yùn)和河口生態(tài)系統(tǒng)。近年來在錢塘江河口治江圍縮窄、建橋等人類活動(dòng)對(duì)涌潮的影響以及涌潮的水力學(xué)結(jié)構(gòu)方面進(jìn)行了大量研究[1-2]。在其他河口,如長江口北支涌潮的形成及其影響因素亦取得了較大進(jìn)展[3]。

東海潮波進(jìn)入杭州灣后,潮差急速增大,傳至灣頂澉浦,潮差與水深的比值已在2倍以上。受淺水效應(yīng)影響,潮波非線性變形加劇,在澉浦上游形成水位驟然升高的漲潮波前鋒線,即為涌潮。涌潮形成之初,潮后的潮位仍在上漲,因此強(qiáng)度逐漸增大,直到涌潮成為整個(gè)漲潮面,潮后潮位隨即下降,涌潮較強(qiáng)時(shí),分裂成前、后兩股,逐漸減弱,直到湮滅[4]。與普通的潮波傳播相比,涌潮前后存在水位、流速的突變,是典型的淺水間斷流動(dòng),比潮波尺度要小好幾個(gè)量級(jí),因而給物理模擬和數(shù)值模擬都造成相當(dāng)大的困難。從20世紀(jì)60年代開始,就有學(xué)者使用激波裝配法和激波捕捉法探索錢塘江河口涌潮的一維或平面二維數(shù)值模擬[5-9]。2001年以來,潘存鴻等[10-15]提出了水位床底法同時(shí)結(jié)合底坡源項(xiàng)離散技術(shù),建立了四邊形網(wǎng)格和三角形網(wǎng)格下二維涌潮數(shù)值模型的Godunov格式和KFVS格式(kinetic flux vector splitting)。數(shù)值模擬給出精細(xì)的流速場(chǎng),彌補(bǔ)了因測(cè)量困難,對(duì)流速分布狀況知之不多的缺憾。然而,到目前為止,還沒有對(duì)錢塘江涌潮進(jìn)行三維數(shù)值模擬研究的報(bào)道。嘗試使用現(xiàn)今國際上最新的河口、陸架和海洋模型FVCOM進(jìn)行錢塘江河口涌潮的三維數(shù)值模擬,加深對(duì)錢塘江河口涌潮傳播特性的認(rèn)識(shí)。

1 模型建立與資料搜集

FVCOM模型為美國MassachusettsDartmouth州立大學(xué)陳長勝所領(lǐng)導(dǎo)的研究小組開發(fā)[16],其控制方程類似于POM模型,包括自由表面、非線性平流項(xiàng)、耦合的密度和速度場(chǎng)、徑流、垂直混合的2.5階湍流閉合模型等,數(shù)值模型采用有限體積法,可應(yīng)用于各種河口、海灣、陸架和海洋問題。FVCOM模型采用σ垂向坐標(biāo)和水平三角形非結(jié)構(gòu)網(wǎng)格坐標(biāo)的結(jié)合,可以使模型應(yīng)用對(duì)感興趣區(qū)域處的網(wǎng)格進(jìn)行加密處理,既可控制計(jì)算量,又不犧牲笛卡爾坐標(biāo)的特性。FVCOM采用類似POM的外模、內(nèi)模分裂的模型求解。二維外模數(shù)值格式為基于三角形網(wǎng)格的有限體積法,將連續(xù)方程、動(dòng)量方程在三角形區(qū)域積分后,通過改進(jìn)的四階龍格-庫塔法求解。三維內(nèi)模的動(dòng)量方程求解采用簡單的顯式和隱式相結(jié)合的差分格式求解,其中流速的局部變換采用一階精度的迎風(fēng)格式,對(duì)流項(xiàng)采用二階精度的改進(jìn)龍格-庫塔時(shí)間推進(jìn)格式,垂向擴(kuò)散則用隱式求解。此外,對(duì)于潮間帶的處理,FVCOM引入了干/濕網(wǎng)格技術(shù)。

2007年10月大潮期間,組織了一次大規(guī)模錢塘江涌潮觀測(cè),上游富春江電站、下游澉浦?jǐn)嗝妗Ｓ^測(cè)期間設(shè)置了20個(gè)潮位站,本研究范圍內(nèi)有11個(gè)(見圖1),無涌潮時(shí)每小時(shí)記錄一次潮位數(shù)據(jù),涌潮到達(dá)后每1～2分鐘記錄一次數(shù)據(jù)。使用ADP進(jìn)行流速觀測(cè),設(shè)置7個(gè)觀測(cè)站位,本研究范圍內(nèi)有3個(gè)測(cè)站,它們是703站、705站、709站(圖1),每個(gè)站位使用6點(diǎn)法進(jìn)行觀測(cè),即分別位于水面、0.2倍水深、0.4倍水深、0.6倍水深、0.8倍水深和海底。無涌潮時(shí),每小時(shí)記錄一次流速數(shù)據(jù),涌潮到達(dá)前后每分鐘記錄一次流速數(shù)據(jù)。觀測(cè)期間最高潮位6.05 m,最高平均高潮位5.41 m,出現(xiàn)在倉前測(cè)站;最低潮位-3.50 m,最低平均低潮位-3.07m,出現(xiàn)于蓋北中沙臨時(shí)站。漲潮歷時(shí)向西逐漸縮短,相應(yīng)落潮歷時(shí)延長。在聞家堰以下河段,平均漲、落潮歷時(shí)差約3小時(shí),至杭州段,平均漲潮歷時(shí)只有約1小時(shí),而平均落潮歷時(shí)長達(dá)11小時(shí)以上。實(shí)測(cè)最大測(cè)點(diǎn)漲潮流速5.83m/s,出現(xiàn)在709測(cè)點(diǎn)表層,最大落潮流速3.84m/s,出現(xiàn)在709測(cè)點(diǎn)0.2倍水深層。最大漲、落潮流速在垂向上一般自表層向底層逐漸減弱,即表層最大,中層次之,底層最弱。通過觀測(cè),對(duì)錢塘江涌潮傳播、發(fā)展規(guī)律有了進(jìn)一步認(rèn)識(shí),同時(shí)也為數(shù)值模擬研究提供了詳盡的資料。

圖1 錢塘江河口Fig.1 Qiantang Estuary

計(jì)算范圍從澉浦至倉前的72 km范圍,河寬從澉浦站20 km縮窄到上邊界倉前2 km(圖1),進(jìn)行岸邊界平滑處理以后,使用無結(jié)構(gòu)網(wǎng)格離散計(jì)算域,模型由17 409個(gè)節(jié)點(diǎn)和32 298個(gè)三角形單元組成,平均網(wǎng)格大小約300m,單元最小邊長60 m,內(nèi)外模的時(shí)間步長分別取10 s和2 s。地形數(shù)據(jù)使用2007年10月地形數(shù)值化以后插值得到,垂向上分7個(gè)等距的σ層。曼寧系數(shù)漲潮期間取值范圍為0.004～0.01,落潮期間取值范圍為0.007～0.02。模型海面邊界對(duì)應(yīng)風(fēng)應(yīng)力,海底邊界條件對(duì)動(dòng)量方程為引入底摩擦應(yīng)力,對(duì)物質(zhì)輸運(yùn)方程為垂向無通量。上下游開邊界采用實(shí)測(cè)潮位數(shù)據(jù)。在模型驗(yàn)證與分析的基礎(chǔ)上,通過數(shù)值試驗(yàn)對(duì)網(wǎng)格大小和曼寧系數(shù)對(duì)涌潮計(jì)算的影響進(jìn)行敏感性分析。

2 模型驗(yàn)證與應(yīng)用

2.1 潮位驗(yàn)證與平面二維模型計(jì)算結(jié)果

驗(yàn)證計(jì)算模擬了錢塘江涌潮形成、發(fā)展、及其衰減的過程,限于篇幅,文中繪出了2007年10月16～17日鹽官、大缺口、曹娥江口三站的潮位驗(yàn)證過程,示于圖2。由圖可知,潮位的驗(yàn)證情況較好,高低潮位、潮差和位相均與實(shí)際吻合。當(dāng)涌潮到達(dá)時(shí)水位上漲很快,這一跳躍現(xiàn)象在模型中被復(fù)演(圖3),說明FVCOM模型能夠捕捉到涌潮前后的水位間斷特征。在計(jì)算域內(nèi),涌潮逐漸形成,并向上游逐漸增強(qiáng),這一現(xiàn)象在模型中清楚地被重現(xiàn)。在下游曹娥江口段,當(dāng)漲潮鋒面到達(dá)時(shí)水位升幅數(shù)厘米。當(dāng)涌潮到達(dá)鹽官站時(shí),涌潮高度達(dá)到最大,之后逐漸衰減,涌潮高度逐漸降低。

圖2 潮位驗(yàn)證Fig.2 Tidal level validations

平面二維模型模擬結(jié)果與實(shí)際情況較為吻合,703、705、709三站漲落潮的最大計(jì)算流速及流速過程線與實(shí)測(cè)接近,受篇幅限制,這里不展示流速驗(yàn)證結(jié)果。落潮垂線平均流速大小在1～2 m/s之間,但是在涌潮到達(dá)以后漲潮垂向平均流速可以達(dá)到4～5 m/s?？傮w上看,垂向平均的最大流速向上游逐漸減小,三個(gè)站位的最大漲潮垂向平均流速分別為 5.0、4.0、3.7 m/s。涌潮傳播過程中在平面上呈現(xiàn)一些特殊形態(tài),習(xí)稱“潮景”,數(shù)值計(jì)算再現(xiàn)了部分“潮景”,如交叉潮、一線潮、回頭潮?？傮w上,潮位、流速、潮景的二維計(jì)算結(jié)果與基于Godnov格式和KFVS格式的二維數(shù)值模型[11-12]的計(jì)算結(jié)果基本一致。

2.2 三維模型計(jì)算結(jié)果

圖4為709站表、中、底三層流速流向計(jì)算驗(yàn)證情況。結(jié)果表明,中層流速與實(shí)際情況吻合較好,表層流速與實(shí)際相比略小,而底層流速比實(shí)測(cè)值略大。盡管三維計(jì)算結(jié)果能夠反映表層流速最大,中層次之,底層流速最小的流速向下遞減的分布特征,但是各層差異不如實(shí)際情況明顯。表層漲落潮最大流速實(shí)測(cè)值分別為5.51m/s和2.89m/s,計(jì)算值則分別為4.98m/s和1.97m/s,底層漲落潮最大流速實(shí)測(cè)值分別為3.56m/s和1.74 m/s,計(jì)算值則分別為4.40 m/s和1.37 m/s。

涌潮的三維特性與二維特性的主要差異表現(xiàn)為流速在垂線上瞬時(shí)分布上的不同。二維在垂線上均勻,三維在垂線上是變化的。在轉(zhuǎn)流過程中先是底部增大,垂線上逐步過渡為滿足對(duì)數(shù)分布。數(shù)值模擬結(jié)果復(fù)演了這一過程。圖5展示709站在從落潮到漲潮的轉(zhuǎn)流過程各σ層垂向流速分布。轉(zhuǎn)流過程的垂線變化一方面是由水質(zhì)點(diǎn)運(yùn)動(dòng)的慣性引起的,即落潮時(shí)期表層流速大,因而克服慣性慢,底層流速小,易于克服慣性,另一方面是由于下游鹽度密度較大引起的。

為分析垂向流速分布的影響因素,從網(wǎng)格粗細(xì)和糙率系數(shù)兩個(gè)方面進(jìn)行了敏感性數(shù)值試驗(yàn)。通常,計(jì)算區(qū)域網(wǎng)格大小的選取取決于模型計(jì)算工作量和對(duì)計(jì)算結(jié)果精度的要求。由于涌潮流速時(shí)空分布的不均勻性,大流速可能只出現(xiàn)某一局部區(qū)域,且大流速出現(xiàn)的時(shí)間有先后。因此,較粗的計(jì)算網(wǎng)格難以捕捉到大的涌潮流速。上述模型的網(wǎng)格大小約300m,在此基礎(chǔ)上將709站附近的局部區(qū)域網(wǎng)格加密到70m左右,其他參數(shù)保持不變。結(jié)果表明,網(wǎng)格加密以后,709站的漲落潮最大流速分別增加了2 m/s和0.5 m/s左右(圖 6)。

圖3 鹽官涌潮到達(dá)前后潮位的實(shí)測(cè)值與計(jì)算值Fig.3 Observed and calculated tidal levels at Yanguan when the bore arrives

底床糙率系數(shù)是模擬潮流運(yùn)動(dòng)的一個(gè)重要參數(shù),但是在現(xiàn)場(chǎng)很難對(duì)糙率系數(shù)進(jìn)行直接觀測(cè),因此糙率系數(shù)是在數(shù)值模擬中初始條件和邊界條件都較為準(zhǔn)確的情況下進(jìn)行模型校正的一個(gè)重要方式[17]。本模型使用曼寧系數(shù)表示糙率,漲潮和落潮階段的曼寧系數(shù)分別取0.004～0.01和0.007～0.02,這種情況下,流速的垂向差異并不明顯。在不改變其他參數(shù)的前提下,將漲落潮曼寧系數(shù)增大至0.01～0.02。結(jié)果顯示。隨著曼寧系數(shù)的增大,各層的垂向流速都有所減小,漲潮流速從4.40～4.98m/s減小為1.40～2.30 m/s,落潮流速從1.37～1.97 m/s減小為0.96～1.55 m/s。各層流速大小的差異則有所增加,表底層漲落潮流速的差異從12%和30%分別增大到39%和38%(表1)。說明曼寧系數(shù)能夠在較大程度上影響流速的垂向分布特征。

圖4 709站流速流向驗(yàn)證Fig.4 Comparison of computed andmeasured three-dimensional velocities at 709

表1 不同曼寧系數(shù)取值下709站表、中、底層漲落潮流速最大值Tab.1 Maximal flood and ebb velocitiesat surface,m iddle and bottom of 709 with differentmanning coefficients

圖5 709站落潮到漲潮的轉(zhuǎn)流過程各層流速分布Fig.5 Current velocities atσlayerswhen ebb currents turn to flood currents at 709

圖6 709站模型加密前后流速-時(shí)間曲線Fig.6 Depth-averaged velocity vs.time at 709

3 結(jié) 語

使用FVCOM模型進(jìn)行了錢塘江河口涌潮的三維數(shù)值模擬,結(jié)果表明FVCOM可以捕捉到涌潮到達(dá)時(shí)的流速和水位突變過程。當(dāng)選用較小的曼寧系數(shù)時(shí)可以使平面二維模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)吻合良好,典型的潮景被再現(xiàn),這與前人基于Godunov和KFVS格式的平面二維數(shù)值模擬結(jié)果基本一致。三維計(jì)算再現(xiàn)了流速在垂向上的差異,即表層最大、中層次之,底層最小。在落潮轉(zhuǎn)向漲潮的過程中,先是底部流速增大,垂向上逐步過渡為滿足對(duì)數(shù)分布。然而模型結(jié)果在各層流速的差異上不如實(shí)測(cè)數(shù)據(jù)明顯。數(shù)值試驗(yàn)顯示,模型網(wǎng)格粗細(xì)可以在較大程度上影響計(jì)算流速大小,當(dāng)網(wǎng)格大小從300m加密到70 m,局部流速增大約1 m/s。在小的曼寧系數(shù)下,流速的垂向差異并不明顯,僅在高潮位憩流期間各層流速差異較大。如果選用較大的曼寧系數(shù),流速在垂向上的差異增大,但是流速減小。此外,為準(zhǔn)確求解涌潮的傳播速度,必須應(yīng)用守恒型計(jì)算方法。數(shù)值實(shí)驗(yàn)表明[19],由于FVCOM的控制為非守恒型,涌潮的傳播速度與時(shí)間步長有關(guān)。

涌潮是一種極為特殊的水流運(yùn)動(dòng),流速大、水流特性復(fù)雜。錢塘江涌潮到達(dá)時(shí),流速一般為6～7m/s,實(shí)測(cè)最大測(cè)點(diǎn)流速達(dá)12 m/s。前人計(jì)算結(jié)果表明,往往需要選用較小的糙率系數(shù)以使計(jì)算流速與實(shí)測(cè)流速吻合。然而,三維數(shù)值模擬結(jié)果顯示,小的曼寧系數(shù)使得垂向流速分布不顯著。其可能的原因是:水流運(yùn)動(dòng)方程中的阻力項(xiàng)采用謝才公式或修正形式,這種粗糙的處理使得原本定義明確的糙率系數(shù)變得包羅萬象,相應(yīng)地物理意義也變得極其含糊。因此,要較好地模擬涌潮作用下垂向流速的差異,很大程度上還有賴于糙率研究的進(jìn)展。

[1] 魯海燕,潘存鴻,盧祥興,等.錢塘江海寧三期治江圍涂工程對(duì)涌潮影響的數(shù)值模擬研究[J].水動(dòng)力學(xué)研究與進(jìn)展:A輯,2008,23(5):484-491.

[2] 楊火其,潘存鴻,周建炯,等.涌潮水力學(xué)特性試驗(yàn)研究[J].水電能源科學(xué),2008,26(4):136-138.

[3] 陳沈良,谷國傳,劉勇勝.長江口北支涌潮的形成條件及其初生地探討[J].水利學(xué)報(bào),2003(11):30-36.

[4] 林炳堯,黃世昌,毛獻(xiàn)忠,等.錢塘江河口潮波變化過程[J].水動(dòng)力學(xué)研究與進(jìn)展:A輯,2002,17(6):665-675.

[5] 金旦華,劉國俊,周本華.一維涌潮計(jì)算[J].應(yīng)用數(shù)學(xué)與計(jì)算數(shù)學(xué),1965,3(2):183-195.

[6] 趙雪華.錢塘江涌潮的一維數(shù)學(xué)模型[J].水利學(xué)報(bào),1985(1):50-54.

[7] 譚維炎,胡四一,韓曾萃,等.錢塘江涌潮的二維數(shù)值模擬[J].水科學(xué)進(jìn)展,1995,6(2):83-93.

[8] 蘇銘德,徐 昕,朱錦林.數(shù)值模擬在錢塘江涌潮分析中的應(yīng)用-Ⅰ數(shù)值計(jì)算方法[J].力學(xué)學(xué)報(bào),1999,31(5):521-533.

[9] 蘇銘德,徐昕,朱錦林.數(shù)值模擬在錢塘江涌潮分析中的應(yīng)用-Ⅱ計(jì)算結(jié)果分析[J].力學(xué)學(xué)報(bào),1999,31(6):700-16.

[10] HuiW H,Pan CH.Water level-bottom topography for the shallow-water flow with application to the tidal boreson the Qiantang River[J].Computational and Dynamics Journal,2003,12(3):549-55.

[11] 潘存鴻,林炳堯,毛獻(xiàn)忠.求解二維淺水流動(dòng)方程的Godunov格式[J].水動(dòng)力學(xué)研究與進(jìn)展:A輯,2003,l8(1):16-23.

[12] 潘存鴻,徐 昆.三角形網(wǎng)格下求解二維淺水方程的 KFVS格式[J].水利學(xué)報(bào),2006,37(7):858-864.

[13] 潘存鴻,林炳堯,毛獻(xiàn)忠.錢塘江涌潮二維數(shù)值模擬[J].海洋工程,2007,25(1):50-56.

[14] 潘存鴻,魯海燕,曾 劍.錢塘江涌潮特性及其數(shù)值模擬[J].水利水運(yùn)工程學(xué)報(bào),2008(2):1-9.

[15] Pan C H,Lin B Y,Mao X Z.Case study:Numericalmodeling of the tidal bore on the Qiantang River,China[J].Journal of Hydraulic Engineering,2007,133(2):130-138.

[16] Chen C,Liu H,Beardsley R.An unstructured grid,finite-volume,three-dimensional,primitive equationsoceanmodel:Application to coastal ocean and estuaries[J].Journal of Atmospheric and Oceanic Technology,2003,20(1):159-186.

[17] You ZJ.Estimation of bed roughness from mean velocitiesmeasured at two levels near the seabed[J].Continental Shelf Research,2005,25(9):1043-1051.

[18] 倪晉仁,薛 安,秦華鵬.影響水流運(yùn)動(dòng)研究進(jìn)展的若干關(guān)鍵問題I.糙率研究[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào),1997,5(4):371-380.

[19] 潘存鴻.淺水間斷流動(dòng)數(shù)值模擬研究進(jìn)展[J].水利水電科技進(jìn)展,2010,30(5):77-84.

Three-dimensional numericalmodeling of tidal bore in Qiantang based on FVCOM

XIEDong-feng,PAN Cun-hong,WU Xiu-guang
(Zhejiang Institute of Hydraulics and Estuary,Hangzhou 310020,China)

TV142

A

1005-9865(2011)01-0047-06

2010-04-08

國家自然科學(xué)基金資助項(xiàng)目(41006053,40806037,50809062);國家水體污染控制與治理科技重大專項(xiàng)(2009ZX07424-001);浙江省創(chuàng)新團(tuán)隊(duì)建設(shè)資助項(xiàng)目(2009F20024);水利部公益性行業(yè)科研專項(xiàng)經(jīng)營資助項(xiàng)目(200801072)

謝東風(fēng)(1981-),男,四川安岳人,博士,主要從事海洋沉積動(dòng)力學(xué)研究。E-mail:eastwind111@hotmail.com