殷亞明,張金善,宋志堯

(1.江蘇省交通規(guī)劃設(shè)計(jì)院,江蘇 南京 210005;2.南京水利科學(xué)研究院,江蘇 南京 210029;3.南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210097)

河口海岸二維懸沙數(shù)學(xué)模型挾沙力公式修正及應(yīng)用

殷亞明1,張金善2,宋志堯3

(1.江蘇省交通規(guī)劃設(shè)計(jì)院,江蘇 南京 210005;2.南京水利科學(xué)研究院,江蘇 南京 210029;3.南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室,江蘇 南京 210097)

在河口海岸二維懸沙數(shù)學(xué)模型中,半潮或全潮平均的潮流挾沙能力并不能完全反映其逐時(shí)的實(shí)際情況,為此進(jìn)行了分析研究。研究結(jié)果表明:1)應(yīng)用半潮或全潮平均的潮流挾沙力公式進(jìn)行懸沙數(shù)值模擬時(shí),應(yīng)引入一個(gè)修正系數(shù);2)該修正系數(shù)與公式中的指數(shù)成反比,與相對(duì)水深成正比,且與相位角也存在一定聯(lián)系;3)該修正系數(shù)理論正確、實(shí)用方便、驗(yàn)證成果良好,實(shí)例計(jì)算表明能夠較好地復(fù)演天然懸沙場(chǎng)分布。

河口海岸;挾沙力;懸沙模型;修正系數(shù)

Abstract:The paper iswritten to dealwith the problem that the integralmean ormean of sediment transport capacity couldn′t reflect the instantaneous capacity of two-dimensional suspended sedimentmathematicalmodel in estuarine and coastalwaters.The results show that:1)A correction coefficient should be introduced into the integralmean ormean of sediment transport capacity formula for the numerical simulation of suspended sediment.2)The correction coefficient is proportional to the index of the formula,inverse proportional to relativewater depth,and also has some connectionwith the phasic angle.3)The correction coefficient isprecise,practical,well-grounded and convenient.Themethod that this paper has introduced could effectively reflect the natural distribution of sediment concentration by practical calculation.

Key words:estuarine and coast;sediment carrying capacity;numericalmodelof suspended sediment;correction coefficient

水流挾沙力是反映底床沖淤平衡狀態(tài)下水流挾帶泥沙能力的綜合性指標(biāo)[1]。天然水體(如湖泊、河流、河口、海岸及海洋)的流動(dòng)都或多或少地挾帶著泥沙,通常將一定水流條件下單位水體能夠挾帶的泥沙含量稱為水流的挾沙力。在涉水工程中,無論是水工建筑物的規(guī)劃設(shè)計(jì)、底床演變的沖淤分析,還是確定物理模型的相似比尺或數(shù)學(xué)模型的沖淤函數(shù)等,都必須考慮水流挾沙能力[2]。

基于挾沙力概念的二維懸沙輸運(yùn)方程一般可寫:式中:s為懸移質(zhì)含沙量,D為總水深,t為時(shí)間,u、v分別為x、y方向的水深平均流速,εx、εy分別為x、y方向懸沙的水平擴(kuò)散系數(shù),α為沉降機(jī)率或恢復(fù)飽和系數(shù),ωf為懸沙沉速,s＊為水流挾沙力公式。

國(guó)內(nèi)海岸河口泥沙運(yùn)移和底床沖淤變化預(yù)測(cè)的數(shù)值模擬中,潮流挾沙能力更是關(guān)鍵問題之一[3]。由于潮流運(yùn)動(dòng)的非恒定性,現(xiàn)有的潮流挾沙力公式大都依據(jù)恒定水流移植的,為了獲得較高的相關(guān)系數(shù),一般通過水深、流速和含沙量的現(xiàn)場(chǎng)數(shù)據(jù)的半潮或全潮平均值確定其系數(shù)[4]。常用的挾沙能力公式有以下兩種[5]:

式中:S＊為半潮或全潮平均的潮流挾沙能力,K為包含量綱的系數(shù),U為半潮或全潮平均流速,h為半潮或全潮平均水深,g為重力加速度,m為指數(shù)(一般在0.3～2.0之間)。

在實(shí)際模擬計(jì)算時(shí),通常認(rèn)為式(2)或(3)中的K、m不變,并以瞬時(shí)流速u來替代U,總水深D來替代h,以此來構(gòu)建潮流挾沙力公式s＊,對(duì)式(3)就有

如此簡(jiǎn)單處理,盡管通過調(diào)試模型參數(shù)α可得到定性上正確、定量上合理的模擬結(jié)果,但其科學(xué)合理性迄今并未進(jìn)行相關(guān)研究。我們認(rèn)為半潮或全潮平均的潮流挾沙能力并不能完全反映其瞬時(shí)的實(shí)際情況,這樣處理可能引起一個(gè)潮周期T內(nèi)的潮流挾沙能力不守恒,即

為此,通過在s＊和S＊之間引入比例系數(shù)β來修正河口海岸二維懸沙數(shù)學(xué)模型,從理論上討論并得出取值的影響因素之間的關(guān)系,為懸沙數(shù)值模擬提供科學(xué)合理的依據(jù)。

1 理論分析

考慮到在河口海岸水域,式(3)的應(yīng)用更為普遍,因此對(duì)式(3)進(jìn)行深入理論分析是不失一般性的。由式(3)可知,潮流挾沙能力的影響因素主要是潮位和潮流,其次是式中的系數(shù)[6]K和m。

1.1 潮汐及潮流

作為理論分析,可假定潮位與流速為一單頻過程,并由于摩阻等作用而存在一個(gè)相位差φ,理論研究及實(shí)測(cè)資料分析可知取值范圍在0～π/2之間,基此可知流速u、潮位η和總水深D的表達(dá)式:

式中:um為流速振幅;ηm為潮位振幅;ω=2π/T為角頻率;T為潮周期。

由此可知,半潮平均潮流速為

1.2 修正系數(shù)

為了保證一個(gè)潮周期T內(nèi)的潮流挾沙能力的守恒性,有必要在瞬時(shí)挾沙力s＊和半潮平均挾沙力S＊之間引入修正系數(shù)β,使得

將式(5)、(6)和(7)代入式(10),可得

1.3 修正系數(shù)β的數(shù)值計(jì)算

由式(11)可以看出,修正系數(shù)β與變量cosφ、ηm/h、m有關(guān),綜合分析了這三個(gè)變量隨β的變化關(guān)系,可數(shù)值積分的方法求出一系列β的數(shù)值解。圖1給出了當(dāng)φ=0、π/4、π/2時(shí),不同相對(duì)水深ηm/h情況下,β隨m的變化關(guān)系曲線。

圖1 不同相對(duì)水深情況下修正系數(shù)β與m的變化關(guān)系曲線Fig.1 The relation curve between correction coefficient andmat differentwater depth

圖1(a)中20條曲線分別代表相位角φ=0時(shí),相對(duì)振幅ηm/h從0.1變化到2.0時(shí)的曲線。從圖中可以看出,β與相對(duì)振幅ηm/h成正比;當(dāng)相對(duì)振幅較小時(shí),β值隨m的變化先增后減,相對(duì)振幅較大時(shí),β值隨m的增大迅速增加。相位差為0時(shí)的β值變化曲線均呈現(xiàn)拋物線型圖形。

圖1(b)中的20條曲線明顯呈現(xiàn)出兩種形式:一種是拋物線型;另一種是折線型。當(dāng)相對(duì)振幅ηm/h較小時(shí),β值隨m的變化先增后減,為拋物線型曲線;相對(duì)振幅ηm/h較大時(shí),β值在m=0.02左右時(shí)發(fā)生一個(gè)跳躍,然后隨m的增大迅速變大,曲線呈現(xiàn)折線形態(tài)。其中ηm/h=1.4與1.5的曲線在m=0.6左右時(shí)有所交叉,是因?yàn)棣莔/h=1.5的曲線在m=0.02左右時(shí)已經(jīng)有一個(gè)不明顯跳躍,為不明顯折線型曲線。因此ηm/h=1.5曲線是拋物線型與折線型曲線的分水嶺。

從圖1(c)可以看出,拋物線型與折線型曲線區(qū)分相當(dāng)明顯。當(dāng)ηm/h>1時(shí)為折線型曲線;ηm/h<1時(shí)為拋物線型曲線。但ηm/h=1的曲線與ηm/h=0.9曲線有所交叉,曲線形式不是很明顯,處于過渡狀態(tài)。

綜合比較φ=0、π/4與π/2的圖形可以看出,β值隨φ從0～π/2的變化是先減后增的。φ=0時(shí),20條曲線均為拋物線型;φ=π/4時(shí),ηm/h>1.5的曲線變?yōu)檎劬€型,且β值比φ=0時(shí)整體減小;φ=π/2時(shí),ηm/h>1的曲線變?yōu)檎劬€型,β值又整體增大。

2 修正系數(shù)的擬合

采用曲面擬合的最小二乘方法,可以推導(dǎo)出β的通用公式。曲面的擬合是借用曲線擬合的最小二乘法思想并推廣到多元的情形,再求解高斯矩陣的解可得擬合系數(shù)。由β的解析解關(guān)系式(12)可以看出,β的取值與相對(duì)振幅ηm/h、相對(duì)水深φ及指數(shù)m有關(guān)。依據(jù)解析解中β的形式,假定m=0.5、1.0、2.0時(shí)的擬合關(guān)系式。

擬合結(jié)果比較滿意,相關(guān)系數(shù)高達(dá)0.996以上,系數(shù)取值如表1所示。

表1ηm/h<1擬合結(jié)果Tab.1 Thefittedresultswhileηm/h<1

由于挾沙力修正公式的推導(dǎo)過程中采用了ηm/h<1時(shí)的假定,故ηm/h>1時(shí)的擬合關(guān)系式不能采用式(13)的形式。由式(13)所求的β值還應(yīng)當(dāng)與m有關(guān),因此假設(shè)擬合關(guān)系式為式(14)形,并推廣到ηm/h>1的情形。

曲面擬合結(jié)果較為滿意(見表2),基本能反映β值隨三個(gè)自變量ηm/h、φ和m的變化情形,擬合公式的最終形式:

表2 擬合結(jié)果Tab.2 Thefittingresults

3 應(yīng)用實(shí)例

采用2006年1月冬季大潮期間5個(gè)潮位站和16個(gè)流速及含沙量觀測(cè)站的實(shí)測(cè)資料進(jìn)行驗(yàn)證,在平面二維水流泥沙數(shù)學(xué)模型中,代入式(15)中修正系數(shù)的擬合結(jié)果進(jìn)行計(jì)算,并對(duì)比分析了傳統(tǒng)數(shù)值模型與本模型的差異來證明β值的準(zhǔn)確性。

3.1 流場(chǎng)驗(yàn)證

連云港海區(qū)小模型的邊界由東中國(guó)海潮波大模型提供,小模型采用非結(jié)構(gòu)三角形網(wǎng)格建立二維潮流數(shù)值模型(圖2),在模型驗(yàn)證過程中,對(duì)邊界處的潮位和流速進(jìn)行微調(diào)以達(dá)到驗(yàn)證相似。邊界條件采用流量控制,即采用水位和流速聯(lián)合控制。潮流數(shù)學(xué)模型計(jì)算時(shí)糙率n=0.015～0.030,糙率在連云港外航道區(qū)適當(dāng)減小。模型網(wǎng)格單元數(shù)26 679,節(jié)點(diǎn)數(shù)14 042。時(shí)間步長(zhǎng)取為1 s,紊動(dòng)粘性系數(shù)取為15～25,動(dòng)邊界控制水深為0.01 m。

1#點(diǎn)布置在海頭附近與2#、3#點(diǎn)在同一斷面上,4#在秦山島附近與5#、16#號(hào)測(cè)站同一斷面上,8#點(diǎn)布置在敦子口附近與9#、10#點(diǎn)在同一斷面上,16#點(diǎn)位于灌河口外海-10m等深線處,為典型的旋轉(zhuǎn)流,驗(yàn)證結(jié)果相當(dāng)吻合。由于各方面復(fù)雜因素影響,除少數(shù)點(diǎn)的驗(yàn)證情況稍差,其余絕大部分驗(yàn)證點(diǎn)驗(yàn)證結(jié)果均較好,說明模型確定的計(jì)算范圍合理,東中國(guó)海潮波數(shù)學(xué)模型提供的邊界條件相對(duì)精確,連云港小模型能夠正確反映研究海域的潮流運(yùn)動(dòng)特征。

3.2 含沙量驗(yàn)證

基于上述考慮,河口、海岸二維懸沙數(shù)學(xué)模型修正形式為

式中:s為垂線平均含沙量;εx、εy為紊動(dòng)擴(kuò)散系數(shù);α為沉降幾率,α=0.2;ωf為泥沙沉速,取0.000 5 m/s;s＊為半潮或全潮平均的潮流挾沙能力;β為修正系數(shù)。圖3為潮流計(jì)算結(jié)果推算出的β分布場(chǎng)。

圖2 連云港區(qū)泥沙數(shù)值模型的測(cè)點(diǎn)布置Fig.2 Arrangement of observation stations for numerical simulation of suspended sediment in Lian Yungang

圖3 修正系數(shù)場(chǎng)Fig.3 Correction coefficient field

分析連云港區(qū)2006年冬季大潮16個(gè)水文觀測(cè)站的實(shí)測(cè)潮流及含沙量報(bào)表,率定出半潮平均水流挾沙力式(3)中的系數(shù)K=19.68、m=0.838 5。整個(gè)連云港區(qū)潮位流速的相位差φ=π/2,由潮流數(shù)值模型計(jì)算結(jié)果得出,連云港區(qū)潮差分布大體上由西南向東北遞減,在灌河口處局部地區(qū)明顯增大。由修正系數(shù)場(chǎng)分布圖3可以看出,修正系數(shù)值基本與等深線平行,且與懸沙濃度場(chǎng)很相似,從另一側(cè)面反映出懸沙濃度近岸大、遠(yuǎn)岸小,沿等深線呈梯度變化的特征。

利用實(shí)測(cè)期間的潮汐動(dòng)力條件,不考慮風(fēng)浪的影響(由于實(shí)測(cè)期間為無風(fēng)或小風(fēng)),分別使用傳統(tǒng)挾沙力公式和修正后的挾沙力公式數(shù)值模型做計(jì)算,對(duì)現(xiàn)場(chǎng)16個(gè)測(cè)點(diǎn)大潮的含沙量過程進(jìn)行了驗(yàn)證(圖4)。根據(jù)歷史資料及近期的實(shí)測(cè)資料分析,正常天氣下除灌河口外,整個(gè)海域的垂線平均含沙量較小,其中海州灣北部含沙量最小。從含沙量分布圖5中可以看出,灌河口的11#點(diǎn)位于灌河入?？谔?實(shí)測(cè)含沙量達(dá)到3.0 kg/m3左右,由于水深較淺且潮差大,該測(cè)站附近平均β值達(dá)到1.2左右,經(jīng)數(shù)模修正放大后計(jì)算結(jié)果較為滿意。

圖4 2006年冬季大潮含沙量驗(yàn)證Fig.4 Verification of sediment concentration during flood period in 2006

從驗(yàn)證結(jié)果可以看出,除少數(shù)點(diǎn)外,含沙量計(jì)算過程線與實(shí)測(cè)過程線基本吻合,大部分測(cè)點(diǎn)修正后更接近實(shí)測(cè)值。其中1、4、6、7、11號(hào)測(cè)點(diǎn)計(jì)算值比實(shí)測(cè)值偏小,經(jīng)修正計(jì)算放大后更接近真值。第3、13號(hào)測(cè)點(diǎn)計(jì)算值比實(shí)測(cè)值偏大,經(jīng)修正縮小后結(jié)果更為準(zhǔn)確。除了2、12號(hào)測(cè)點(diǎn)不太理想外,其他測(cè)點(diǎn)均與實(shí)測(cè)值符合良好。表3為16個(gè)水文觀測(cè)站分別用本方法和傳統(tǒng)方法所計(jì)算出的潮段平均含沙量平均值標(biāo)準(zhǔn)差δ;表4為16個(gè)水文觀測(cè)站分別用本方法和傳統(tǒng)方法所計(jì)算出的潮段平均最大峰值差,除了個(gè)別號(hào)點(diǎn)使用本方法的標(biāo)準(zhǔn)差比傳統(tǒng)方法略大,其他測(cè)點(diǎn)的標(biāo)準(zhǔn)差都比傳統(tǒng)方法計(jì)算所得的要小很多。因此,所建立的模型能正確反映連云港區(qū)的含沙量分布,并有力地證明了挾沙力公式經(jīng)修正后更為準(zhǔn)確。

圖5 5大潮含沙量等值線Fig.5 The isoline of sediment concentration during flood period

表3 觀測(cè)站平均值標(biāo)準(zhǔn)差Tab.3 The average standard deviation of observation stations

表4 觀測(cè)站最大峰值差Tab.4 The maximum difference of observation stations

4 結(jié) 語

在建立二維懸沙擴(kuò)散方程時(shí)引入了一個(gè)修正系數(shù),并結(jié)合半潮平均的挾沙力公式給出了其計(jì)算表達(dá)式。

結(jié)果表明該修正系數(shù)與相對(duì)振幅、潮流與潮差的相位差及挾沙力公式中的指數(shù)有著一定關(guān)系,由此修正的二維懸沙輸運(yùn)擴(kuò)散方程簡(jiǎn)單實(shí)用,且在理論上很完善地解決了半個(gè)或一個(gè)潮周期內(nèi)的潮流挾沙能力不守恒問題。本修正模型作為水動(dòng)力懸沙數(shù)值模擬的一個(gè)成功嘗試,為進(jìn)一步開展港區(qū)、航道和環(huán)境泥沙問題奠定了基礎(chǔ),對(duì)以后用經(jīng)驗(yàn)挾沙力公式計(jì)算泥沙數(shù)值模型有一定的實(shí)用價(jià)值。

[1] 張瑞瑾.河流泥沙動(dòng)力學(xué)[M].北京:水利電力出版社,1989.

[2] 錢 寧,萬兆惠.泥沙運(yùn)動(dòng)力學(xué)[M].北京:科學(xué)出版社,1983.

[3] 孔 俊,宋志堯,章衛(wèi)勝.挾沙能力公式系數(shù)的最佳確定[J].海洋工程,2005,23(1):93-96.

[4] 邢 云,宋志堯,孔 俊,等.長(zhǎng)江口水流挾沙力公式初步研究[J].水文,2008,28(1):64-66.

[5] 李瑞杰,羅 鋒,周華民.水流挾沙力分析與探討[J].海洋湖沼通報(bào),2009(1):88-94.

[6] 舒安平.水流挾沙力公式的驗(yàn)證與評(píng)述[J].人民黃河,1993(1):7-10.

[7] 李義天,趙明登,曹志芳.河道平面二維水沙數(shù)學(xué)模型[M].北京:中國(guó)水利水電出版社,2001.

Modification and application of numericalmodel of suspended sediment transport in estuarine and coastalwaters

YIN Ya-ming1,ZHANGJin-shan2,SONG Zhi-yao3
(1.Jiangsu Provincial Communications Planning and Design Institute Co.,Ltd.,Nanjing 210005,China;2.Nanjing Hydraulic Research Institute,Nanjing 210029,China;3.Key Laboratory of Virtual Geographic Environment(Ministry of Education),Nanjing Normal University,Nanjing 210097,China)

TV148

A

1005-9865(2011)01-0082-07

2010-05-28

“十一五”交通重大專項(xiàng)資助項(xiàng)目

殷亞明(1984-),男,江蘇鎮(zhèn)江人,碩士,主要從事港口、航道與海岸工程研究。E-mail:283665360@qq.com