一維數(shù)學(xué)模型在堤前灘面短期演變計(jì)算中的應(yīng)用

張景新,王 穎,王 偉,寇 軍,王佳飛,劉 樺

(1.上海交通大學(xué)船舶海洋與建筑工程學(xué)院,上海 200240;2.上?；瘜W(xué)工業(yè)區(qū)物業(yè)管理有限公司,上海 201507)

沿海地區(qū)夏季頻繁的臺(tái)風(fēng)對(duì)海岸防護(hù)工程有諸多危害,涉及越浪、護(hù)岸結(jié)構(gòu)的破壞及堤身的穩(wěn)定等。臺(tái)風(fēng)引起的風(fēng)浪在灘面上傳播,伴隨著波浪的破碎,引起灘面上大量的床沙起動(dòng)、懸浮及輸移。灘面沖淤變化不僅影響著波浪的傳播、爬高及越浪,同時(shí)也影響著堤身結(jié)構(gòu)的穩(wěn)定,甚至引起海岸災(zāi)害。

鑒于海岸防護(hù)工程的重要性,波浪作用下灘面沖淤演變的研究一直以來(lái)備受關(guān)注,國(guó)內(nèi)外近年來(lái)已有大量的研究成果[1-4],涉及試驗(yàn)研究及數(shù)值模擬。試驗(yàn)研究借助于大型水槽設(shè)備的開(kāi)發(fā)及測(cè)量手段的不斷完善,已取得了豐碩的研究成果[4]。數(shù)值模擬手段在海灘輸沙的計(jì)算中也發(fā)揮著重要的作用。近年來(lái),復(fù)雜流體運(yùn)動(dòng)模型或水沙兩相流模型在海岸帶潮流場(chǎng)、波浪場(chǎng)及泥沙運(yùn)動(dòng)的研究中逐漸得以發(fā)展,在水沙運(yùn)動(dòng)的機(jī)理性研究中發(fā)揮著重要的作用[3,5]。但對(duì)于實(shí)際海岸工程而言,大型數(shù)值模型的開(kāi)發(fā)及應(yīng)用還難以達(dá)到工程實(shí)際應(yīng)用的程度。早期發(fā)展的海岸帶橫向輸沙的一維數(shù)學(xué)模型形式相對(duì)簡(jiǎn)單,計(jì)算工作量較小,借助于實(shí)測(cè)資料及合理的經(jīng)驗(yàn)系數(shù)的選取,該類(lèi)模型在實(shí)際工程應(yīng)用中具有一定的適用性。一維數(shù)學(xué)模型應(yīng)用中的首要問(wèn)題是確定橫向輸沙率,目前已有眾多輸沙率計(jì)算方法[6-7]。依據(jù)輸沙的動(dòng)力因素不同,研究者將泥沙在空間上分為波浪上涌帶和波浪淺化帶及破碎帶分別加以考慮[7]。泥沙輸沙率計(jì)算公式繁多,如何選擇或發(fā)展適合特定水動(dòng)力條件下的輸沙率公式非常重要。一維數(shù)學(xué)模型涉及的輸沙率計(jì)算公式多為經(jīng)驗(yàn)型公式,但借助于適當(dāng)?shù)膶?shí)測(cè)資料及經(jīng)驗(yàn)系數(shù)的選取,往往可獲得較滿(mǎn)意的計(jì)算結(jié)果。

筆者建立了波浪作用下灘面演變的一維數(shù)學(xué)模型,波浪場(chǎng)模擬采用波浪能量守恒方程,該方程計(jì)入了波浪破碎效應(yīng)。波浪模型給出波參數(shù)的空間分布,應(yīng)用于輸沙率的計(jì)算。橫向輸沙率計(jì)算公式參考了眾多研究成果,計(jì)入了波浪上涌帶泥沙輸沙率、波浪淺化帶及破碎帶輸沙率及波生流對(duì)泥沙輸運(yùn)的作用。將該模型應(yīng)用于杭州灣某海防工程堤前灘面的沖淤演變計(jì)算,模擬了某次臺(tái)風(fēng)過(guò)境期間風(fēng)浪所引起的灘面沖淤變化。通過(guò)數(shù)值模擬驗(yàn)證了該模型的精度,分析了若干影響因素及該模型的適用性。

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

1.1 波浪場(chǎng)計(jì)算

文中波浪場(chǎng)的計(jì)算基于規(guī)則波,通過(guò)波浪能量(或波作用量)守恒方程的求解,獲得波浪要素的空間分布。近岸水域波浪運(yùn)動(dòng)受地形等因素的影響,運(yùn)動(dòng)形式復(fù)雜,淺化、折射及繞射現(xiàn)象并存,筆者僅考察入射波浪在橫向的傳播演化及其所引起的橫向輸沙,將控制方程簡(jiǎn)化為橫向的一維模型。近岸帶水體運(yùn)動(dòng)受到諸多因素的影響,如床面摩阻作用隨水深的減小逐漸顯著。而波浪破碎是該水域波浪運(yùn)動(dòng)的顯著特征,本文模型僅概化了破波帶內(nèi)波浪運(yùn)動(dòng)的破碎效應(yīng),控制方程如下:

式中:Cg(x)為波浪群速度;E(x)為單位波長(zhǎng)的波能;Hw(x)為波高的空間分布;Dw(x)為由波浪破碎引起的波能耗散項(xiàng);W為破波因子函數(shù)(W建立了空間波高與當(dāng)?shù)厮畹年P(guān)系,通過(guò)設(shè)定破波指標(biāo)可將波浪破碎這一物理過(guò)程模式化,大量應(yīng)用成果表明該處理方法具有相當(dāng)廣泛的適用性),本文破波因子函數(shù)計(jì)算式來(lái)源于文獻(xiàn)[8];K為計(jì)算系數(shù);d為當(dāng)?shù)厮?γ為破波指標(biāo),γ=0.73。

利用有限差分法求解上述波浪能量守恒方程,求得波浪參數(shù)的空間分布,進(jìn)而可計(jì)算床面應(yīng)力等物理量,并進(jìn)一步用于計(jì)算橫向輸沙率。

1.2 床面演化控制方程

床面橫向演化的一線(xiàn)方程給出了局部床面高程的時(shí)間變化與輸沙率橫向變化的關(guān)系,表達(dá)式如下:

式中:z為床面高程;ε為床沙孔隙率;ρs和 ρ分別為泥沙的密度和水的密度;g為重力加速度;q為輸沙率;qR為波浪上涌帶的輸沙率;qw為波浪淺化帶及破碎帶內(nèi)的輸沙率。

海岸帶橫向輸沙率的空間分布如圖1所示。圖中區(qū)域1,2,3分別代表破波帶外、破波帶和波浪上涌帶;xB,xR,xM分別為破波點(diǎn)的橫坐標(biāo)、上涌帶的輸沙率峰值點(diǎn)的橫坐標(biāo)和波浪最大爬升高度處的橫坐標(biāo);dB,dR分別為xB,xR處的水深。觀察圖1可知,破波帶(區(qū)域2)內(nèi)的輸沙強(qiáng)度最大,波浪上涌帶的輸沙強(qiáng)度隨著波浪的爬升逐漸減弱。利用迎風(fēng)格式的有限差分法數(shù)值求解式(2),而式(2)的求解首先需要給出橫向輸沙率的計(jì)算表達(dá)式。

圖1 近岸帶橫向輸沙強(qiáng)度空間分布[7]

1.3 海岸帶橫向輸沙率

海岸帶橫向輸沙率按空間分布可分為2部分,即波浪上涌帶的輸沙率qR和波浪淺化帶及破碎帶內(nèi)的輸沙率qw。參考已有的研究成果[7-8],分別給出各輸沙率的計(jì)算表達(dá)式。海浪上涌帶的qR計(jì)算公式如下:

式中:βeq為床沙平衡角;β為床面坡度;c1為系數(shù),c1=10-3～10-1;R為波浪爬升高度;Hs0為入射波的有效波高;L0為深水波波長(zhǎng)。

上涌帶外側(cè)qR值通過(guò)式(4)計(jì)算得到:

式中:c3為系數(shù),c3=0.2～0.3。

淺化帶及破波帶內(nèi)的波浪輸沙率qw可分為推移質(zhì)輸沙率qb和懸移質(zhì)輸沙率qs,參考已有的研究成果[3,7],分別給出各自的計(jì)算表達(dá)式:

式中:εb為推移質(zhì)輸沙效率 ,εb=0.1;εs為懸移質(zhì)輸沙效率,εs=0.01;Φ為床沙顆粒的內(nèi)摩擦角,tan Φ=0.6;fw為波浪摩阻系數(shù);Ud為水深平均的離岸流速;F和B分別為床面摩阻和波浪破碎引起的動(dòng)能耗散率;T為波周期;Uw為波浪邊界層處水質(zhì)點(diǎn)運(yùn)動(dòng)速度幅值;z0為床面粗糙尺度;ωs為泥沙沉速;d50為床沙中值粒徑。

離岸流體運(yùn)動(dòng)主要分布于下層水體,該流動(dòng)與向岸的波浪場(chǎng)質(zhì)量輸移運(yùn)動(dòng)和破碎波引起的上層水體的向岸運(yùn)動(dòng)相平衡,Ud可由式(7)計(jì)算:

式中:c為波浪相速度;E為波能。

懸移質(zhì)輸沙率計(jì)算表達(dá)式中F及B的計(jì)算參考相關(guān)研究成果[7]:

2 模型應(yīng)用及討論

將模型應(yīng)用于杭州灣某海防工程堤前灘面沖淤演變的數(shù)值模擬研究中,該海防工程位于杭州灣北岸,工程位置及地形測(cè)量斷面位置見(jiàn)圖2(圖中符號(hào)G代表斷面)。2005年“麥莎”臺(tái)風(fēng)前后,分別開(kāi)展了堤前水下地形的測(cè)量工作,2次測(cè)量前后間隔大約3個(gè)月。2次局部地形的變化主要由此次臺(tái)風(fēng)過(guò)境期間的風(fēng)浪所引起。首先根據(jù)現(xiàn)場(chǎng)實(shí)測(cè)資料確定計(jì)算水文條件,臺(tái)風(fēng)風(fēng)向ESE,取入射波的有效波高Hs0=2.7m,相應(yīng)周期T=6.5s,堤前平均潮位 η0=2.0m,波向由實(shí)測(cè)資料確定為正向入射。模型計(jì)算波參數(shù)取相應(yīng)的有效波參數(shù),計(jì)算歷時(shí)24h(實(shí)際臺(tái)風(fēng)過(guò)境歷時(shí)約72h)。堤前灘面床沙中值粒徑d50=0.03mm,將2005年5月地形作為初始地形,采用所設(shè)定的水文參數(shù),模型給出了此次臺(tái)風(fēng)過(guò)境期間風(fēng)浪所引起的沖淤變化量。該海防工程堤前建有混凝土護(hù)岸的緩坡,故在輸沙率的計(jì)算中未考慮波浪上涌帶的泥沙輸沙率。

測(cè)量資料給出了沿岸4個(gè)固定斷面的地形資料,將各斷面的模擬結(jié)果與實(shí)測(cè)資料比較,如圖3所示。圖3分別給出了4個(gè)固定斷面2005年5月和8月的實(shí)測(cè)地形值和臺(tái)風(fēng)后灘面地形的模擬值,x軸負(fù)向?yàn)橥夂７较?正向?yàn)榈躺矸较?。?jì)算結(jié)果顯示在堤腳的鄰近區(qū)域各斷面普遍出現(xiàn)灘面沖刷,而離岸方向普遍出現(xiàn)淤積。實(shí)測(cè)結(jié)果顯示在鄰近堤腳的區(qū)域內(nèi)灘面沖刷明顯,離岸方向沖淤幅度減弱。計(jì)算值與實(shí)測(cè)值在鄰近堤腳的區(qū)域吻合較好。分析流場(chǎng)動(dòng)力因素,隨著水深的減小,入射波波高逐漸增大,直至破碎。破波帶內(nèi)波浪對(duì)泥沙的輸運(yùn)起著非常重要的作用,而在遠(yuǎn)離堤腳的水域,潮流作為泥沙輸運(yùn)的動(dòng)力因素逐漸增強(qiáng)。分析該水域水文實(shí)測(cè)資料,潮流運(yùn)動(dòng)表現(xiàn)為沿岸往復(fù)流動(dòng),與波向近似正交,其輸沙方向?yàn)檠匕丁１疚哪Ｐ蛢H關(guān)注了波浪作用下的泥沙橫向輸運(yùn),未計(jì)及潮流對(duì)泥沙輸運(yùn)的影響。以G1斷面為例,模擬結(jié)果顯示該斷面內(nèi)床沙的沖淤總量平衡,而實(shí)測(cè)值顯示該斷面的總沖刷量大于淤積量,意味著部分起懸泥沙可能為潮流所輸運(yùn),表現(xiàn)為沿岸輸沙。除潮流影響外,破碎波引起的沿岸流動(dòng)也是重要的動(dòng)力因素之一。沿岸流在淺水近岸處受地形等因素影響較弱,泥沙輸沙方向基本維持在斷面內(nèi)。遠(yuǎn)離堤腳,隨著水深的增加,地形等影響作用減弱,沿岸流動(dòng)增強(qiáng),沿岸輸沙強(qiáng)度漸增,這一因素未計(jì)入本模型,也是影響遠(yuǎn)離堤腳處計(jì)算偏差漸增的可能原因之一。

圖2 工程位置圖及地形測(cè)量斷面位置

圖3 灘面地形沖淤變化

3 結(jié) 語(yǔ)

海岸帶橫向泥沙輸運(yùn)的一維數(shù)學(xué)模型以其形式簡(jiǎn)單、計(jì)算高效等特點(diǎn),具有較高的實(shí)用性。模型成功應(yīng)用的前提是橫向輸沙率的準(zhǔn)確計(jì)算。海岸帶流動(dòng)條件復(fù)雜,波浪的淺化、破碎及上涌都對(duì)泥沙運(yùn)動(dòng)起著重要的作用,不同輸沙形式是由不同的動(dòng)力因素決定的。面對(duì)種類(lèi)繁多的輸沙率公式,還有待分析總結(jié)及發(fā)展更適合的計(jì)算模型。深入研究該水域的流體運(yùn)動(dòng),量化各流動(dòng)物理量是改進(jìn)模型的前提。此外,如文中提到的這類(lèi)工程設(shè)施,不僅需要考慮風(fēng)浪,還需考慮潮流作用,這方面的工作還有待深入。并且,文中模型應(yīng)用于實(shí)際海岸的床面沖淤,針對(duì)規(guī)則波,而實(shí)際風(fēng)浪為隨機(jī)波列。雖然對(duì)于近岸床面的沖淤模擬結(jié)果較好,但仍需進(jìn)一步研究隨機(jī)波作用下的灘面沖刷。

:

[1]RAKHA K A,DEIGAARD R,BROKER I.A phase-resolving cross shore sediment transport model forbeach profile evolution[J].Coastal Engineering,1997,31(4):231-261.

[2]SPIELMANN K,ASTRUC D,THUAL O.Analysis of some key parametrizations in a beach profile morphodynamical model[J].Coastal Engineering,2004,51(10):1021-1049.

[3]LEONT'YEV I O.Short-term shoreline changes due to crossshore structures:a one-line numerical model[J].Coastal Engineering,1997,31(1):59-75.

[4]van RIJN L C.Prediction of dune erosion due to storms[J].Coastal Engineering,2008,56(4):441-457.

[5]CHRISTENSEN E D,WALSTRA D J,EMERAT N.Vertical variation of the flow across the surf zone[J].Coastal Engineering,2002,45(3):169-198.

[6]錢(qián)寧,萬(wàn)兆惠.泥沙運(yùn)動(dòng)力學(xué)[M].北京:科學(xué)出版社,1983:270-359.

[7]LEONT'YEV I O.Numericalmodeling of beach erosion during storm event[J].Coastal Engineering,1996,29(2):187-200.

[8]KIRBY J T,DALRYMPLE R A.Modeling waves in surf-zones and around islands[J].Journal of Waterway,Port,Coastal and Ocean Engineering,1986,112(1):78-93.