徑流和風(fēng)對涌潮影響的數(shù)值模擬

張舒羽，潘存鴻

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

涌潮往往發(fā)生在潮差較大的喇叭形河口或海灣，是漲潮波前鋒淺水變形的結(jié)果。涌潮來臨時，其勢兇猛，轟轟作響，潮頭陡立，猶如一道直立的水墻高速前進，壯觀無比。全世界大約有450 個河口或海灣存在涌潮現(xiàn)象［1］，但以錢塘江為最甚。錢塘江涌潮是獨特的、寶貴的自然景觀和旅游資源，每年秋季大潮，吸引數(shù)十萬中外游客前往觀潮。另一方面涌潮現(xiàn)象對兩岸海塘堤防的破壞很嚴(yán)重，危及兩岸廣大人民的生命財產(chǎn)安全。同時，涌潮也給航運安全造成很大威脅。研究涌潮規(guī)律，對于涌潮保護、涌潮防災(zāi)具有實際意義。

影響涌潮特征的因素錯綜復(fù)雜，大致可歸納為潮汐、徑流、江道地形和氣象因素四個方面［2-3］。林炳堯等［4］根據(jù)2000年錢塘江實測資料分析了涌潮傳播規(guī)律，得到了涌潮高度與當(dāng)?shù)爻辈畹恼染€性關(guān)系;趙雪華［5］應(yīng)用一維涌潮數(shù)值模型研究了徑流對涌潮高度的影響，認(rèn)為山水流量增大，對漲潮流有頂托作用，可使涌潮高度稍有增加，但山水流量增大到某一臨界值后，涌潮高度反而減小;曾劍等［6］應(yīng)用神經(jīng)網(wǎng)格模型分析了鹽官站涌潮高度與當(dāng)?shù)爻辈睢⒌匦魏蛷搅鞯亩筷P(guān)系，研究結(jié)果表明，當(dāng)上游徑流量小于7 000 m3/s 時，鹽官涌潮高度隨徑流量增大而增加;當(dāng)上游徑流量超過此值時，涌潮高度隨徑流量增大而減小。氣象因素對涌潮高度也產(chǎn)生一定影響，尤其是在臺風(fēng)暴潮期間影響較大。其影響主要表現(xiàn)在直接影響和間接影響兩個方面:直接影響是氣象因素直接增大或減小涌潮高度，當(dāng)東風(fēng)持續(xù)勁吹時，涌潮高度將增大;間接影響是氣象因素通過增減水，從而增大或減小涌潮河段下游潮差來實現(xiàn)的，氣象因素對涌潮高度的影響還沒有定量的結(jié)論［3］。由于缺乏涌潮特征與影響因素的長系列觀測資料，并且在實際江道中，影響涌潮特征的因素眾多，各影響因素之間并非都是獨立變量，不易得到某些因素(如徑流、風(fēng)向風(fēng)速)對涌潮特征的影響規(guī)律。為此，以錢塘江河口為基礎(chǔ)，建立概化強潮河口的水動力數(shù)學(xué)模型，使復(fù)雜問題相對簡單化，通過數(shù)值模擬探討徑流、風(fēng)向風(fēng)速對涌潮的影響。

1 涌潮模擬的水動力模型

尋找一個合適的數(shù)值模型是模擬涌潮的首要條件。涌潮數(shù)值模擬有兩個關(guān)鍵問題，一是強間斷的模擬。涌潮前后，水位、流速(流量)存在突變，非線性效應(yīng)很強，常規(guī)的數(shù)值方法要么數(shù)值耗散太大，計算結(jié)果將突變的物理量抹平;要么產(chǎn)生數(shù)值振蕩，甚至失穩(wěn)。二是方程源項的處理。為準(zhǔn)確求解涌潮傳播速度，控制方程必須采用守恒型方程，從而帶來了方程左端壓力項與方程右端底坡源項的“和諧”問題［7］。

采用基于無碰撞二維Boltzmann 方程的KFVS(kinetic flux vector splitting)格式求解二維淺水方程，該格式能模擬間斷流動。無碰撞二維Boltzmann 方程:

式中:q 為平衡態(tài)時分子速度分布函數(shù);cx、cy分別為分子在x、y 方向的分子速度;φ 為外力作用項，這里考慮非平底引起的重力、阻力、風(fēng)應(yīng)力、柯氏力和鹽度密度引起的壓力等外力。

式中:g 為重力加速度;Sfx、Sfy分別為x、y 方向的阻力項;S0x、S0y分別為x、y 方向的底坡項;ρ 為鹽水密度;f 為柯氏系數(shù);Wx、Wy分別為x、y 方向的風(fēng)應(yīng)力;Qs為取(排)水口的水流流量;β 為取(排)水方向與x 軸的夾角。

式(1)不考慮碰撞項的結(jié)果是，其推得的對應(yīng)控制方程沒有計及二階粘性項。將式(1)乘以(1，cx，cy)T，并對分子速度空間積分，利用水流宏觀變量水深h，x 和y 方向的流速u、v 與分子微觀變量f、cx、cy的關(guān)系:

可得控制方程

式中:

采用無結(jié)構(gòu)三角形單元離散，并采用網(wǎng)格中心格式。設(shè)Ωi為第i 個三角形單元域，Γ 為其邊界，對式(7)應(yīng)用有限體積法離散，經(jīng)推導(dǎo)可得基本數(shù)值解公式［8］

式中:Ai為三角形單元Ωi的面積;Fn=Fcosθ +Gsinθ，(cosθ，sinθ)為Γ 外法向單位向量;Δt 為時間步長;下標(biāo)j 表示i 單元第j 邊;lj為三角形邊長;上標(biāo)n 為時間步。

求解式(12)的核心是法向數(shù)值通量的計算、底坡源項的處理，以及動邊界處理［3］。這里應(yīng)用能模擬間斷流的干底Riemann 解處理動邊界問題。上述模型已經(jīng)經(jīng)典算例和實際問題的應(yīng)用檢驗［3，8-10］。

2 模型概化

以錢塘江富春江電站至蘆潮港的江段為原型，按照電站、聞家堰、閘口、七堡、倉前、鹽官、澉浦、乍浦、蘆潮港斷面的平均高程和河寬對概化模型進行控制，忽略斷面上河床地形的變化，構(gòu)造一個軸對稱的概化模型，見圖1，沿河道縱向的地形見圖2。模型構(gòu)造需考慮錢塘江涌潮形成的2 個基本條件［2］:一是錢塘江河口在平面上的急劇收縮使潮波能量積聚，潮差增大;二是乍浦以上的水下沙坎堆積，使水深迅速減小，潮波淺水效應(yīng)增強。概化模型電站到蘆潮港全長278 km，其中聞家堰至澉浦長116 km。各控制斷面的河寬采用實際河寬，從下往上平面上急劇收縮。自乍浦以上，河床高程抬高，至倉前達到最高，而后至聞家堰高程降低。由此可見，構(gòu)造的概化模型滿足涌潮發(fā)生的兩個條件。概化模型下邊界蘆潮港采用正弦函數(shù)來模擬潮汐過程，周期12 h，潮差3 m。模型計算域采用三角形網(wǎng)格離散，三角形單元194 432個，節(jié)點99 625 個，最小空間步長50 m，位于關(guān)注區(qū)域聞家堰～鹽官河段，計算時間步長0.5 s。漲潮糙率0.002 ～0.02，落潮糙率0.003 ～0.03。概化模型在鹽官位置涌潮高度可達2.4 m左右。

圖1 概化模型網(wǎng)格示意Fig.1 Grid of generalized model

3 涌潮傳播規(guī)律

當(dāng)上游電站流量為500 m3/s，下游蘆潮港潮差為3 m，無風(fēng)情況下時，繪出反映涌潮傳播規(guī)律的特征圖并闡述規(guī)律如下:

圖3 繪出了高、低潮位、潮差和涌潮高度的沿程變化。鹽官至倉前，河底坡度大，潮波在傳播過程中，高、低水位均向上游迅速抬升，低潮位的抬升值大于高潮位，潮差沿程遞減。倉前至聞家堰，高潮位變化不大，低潮位繼續(xù)緩慢抬升，潮差減小。鹽官至七堡，涌潮高度沿程減小，七堡至聞家堰，涌潮高度略有增加。

圖4 繪出了代表站位潮位隨時間的變化過程。涌潮到達時，水位迅速抬升，如鹽官站水位迅速上升2.42 m，七堡站水位迅速上升1.76 m。

圖2 縱向地形示意Fig.2 Schematic diagram of longitudinal terrain

圖3 高低潮位、潮差和涌潮高度沿程變化Fig.3 Variation of high and low water level，tidal range and tidal bore height along the river

圖4 代表站位潮位過程線Fig.4 Tidal level of representative station

圖5 繪出了每間隔1 h 沿程水面線的變化。從圖中清晰可見每個時刻涌潮到達的位置。涌潮到達時，水面迅速抬升。在圖中所示的3 h 里，涌潮行進了74 km。

圖5 不同時刻水面線Fig.5 Water surface profile at different time

4 不同徑流下涌潮特征變化規(guī)律

上游徑流量采用500、1 000、5 000、7 000、10 000 m3/s 進行計算，各特征值的統(tǒng)計結(jié)果:根據(jù)計算，當(dāng)流量超過7 000 m3/s，概念模型里七堡上游無漲潮流;當(dāng)流量超過10 000 m3/s，倉前上游無漲潮流。各流量下漲落潮差、漲落歷時、涌潮高度特征值見表1。各流量下沿程各站涌潮高度的變化見圖6。一般情況下，上游流量越小，沿程各站潮差越大。當(dāng)流量小于1 000 m3/s 時，漲潮歷時的變化差異不大;當(dāng)流量繼續(xù)增加，漲潮歷時一般隨著流量的增加而減小，但當(dāng)流量超過7 000 m3/s，七堡以上只有水位的抬升和下降，而無漲潮流，當(dāng)上游流量繼續(xù)增加，七堡以上水位抬升的時間隨著流量的增加而增加。涌潮高度的最大值上下游河段有所不同，在下游鹽官、倉前站，涌潮高度在5 000 m3/s時達到最大，在上游七堡、閘口、聞家堰，涌潮高度在1 000 m3/s 達到最大。說明在徑流不大時，涌潮高度隨著徑流的增加而增加，到達一定臨界流量后，水深加大，受下泄流量頂托，潮汐減弱，涌潮高度反而降低，并且降低幅度較大。這與文獻［6］的結(jié)論在定性上是一致的;越往上游，臨界流量值越小。以鹽官、七堡為代表，畫出各流量下潮位隨時間的變化見圖7。涌潮到達時，潮位迅速抬高;流量越大，漲潮時間越滯后。

圖6 各流量下沿程涌潮高度Fig.6 Tidal bore height along the river under different discharges

表1 各流量下涌潮特征值Tab.1 Characteristic values of tidal bore under different runoffs

涌潮傳播速度指涌潮潮頭的行進速度。根據(jù)計算結(jié)果，統(tǒng)計各流量下各區(qū)段的涌潮傳播速度見表2。同一流量下，下游的潮波傳播速度小于上游的速度。同一區(qū)段，隨著流量的增加，大部分情況潮波傳播速度隨著流量的增加而減小，但倉前～七堡在5 000 m3/s 時潮波傳播速度達到最大。根據(jù)文獻［3］分析，涌潮傳播速度c 為:

式中:h 為水深;Δh 為涌潮高度;u 為流速;下標(biāo)d 和u 分別標(biāo)識下游側(cè)和上游側(cè)。

由式(13)可見，涌潮傳播速度與涌潮前落潮流速uu、涌潮高度Δh、潮前水深hu和潮后水深hd等因素有關(guān)。隨著上游徑流量的增大，uu絕對值增大，hu增大，hd變化相對較小，Δh 變化有增有減，如圖6 所示，表2為各流量下涌潮傳播速度的數(shù)模計算結(jié)果，由表可知，鹽官以上河段c 隨著徑流量的增大而減小;而澉浦～鹽官河段，由于徑流量從500 m3/s 增大到5 000 m3/s 時Δh 增大較多，導(dǎo)致在這一流量區(qū)間c 變化不大，而當(dāng)徑流量增大到7 000 m3/s 以后，因uu、hu和Δh 三個因子均促使c 減小，從而c 急劇下降。因此，數(shù)模計算結(jié)果與解析式(13)的分析結(jié)果是一致的。

圖7 代表站位各流量下潮位過程線Fig.7 Tidal level hydrograph of representative stations under different discharge

表2 各流量下涌潮傳播速度Tab.2 Propagation velocity of tidal bore under different discharges

5 不同風(fēng)況下涌潮特征變化規(guī)律

5.1 不同風(fēng)向下涌潮特征變化規(guī)律

加風(fēng)時間為鹽官低潮位前2.5 h，風(fēng)速15 m/s，上游徑流500 m3/s。計算了東風(fēng)、西風(fēng)、無風(fēng)三種工況。加風(fēng)后潮汐特征見表3 和圖8。

由圖8 和表3 可見，風(fēng)向?qū)Τ辈?、涌潮高度、漲潮歷時和涌潮傳播速度影響很大，無論是潮差、涌潮高度還是漲潮歷時和涌潮傳播速度，均是“西風(fēng)(逆風(fēng))”＜“無風(fēng)”＜“東風(fēng)(順風(fēng))”。各站影響程度有所不同，如鹽官站“東風(fēng)”方案潮差比“無風(fēng)”方案大0.83 m，增大17%;“西風(fēng)”方案潮差比“無風(fēng)”方案小0.90 m，減小18%;“東風(fēng)”方案比“西風(fēng)”方案潮差大1.73 m，增大44%。鹽官站“東風(fēng)”方案漲潮歷時比“無風(fēng)”方案延長0∶43，增大46%;“西風(fēng)”方案漲潮歷時比“無風(fēng)”方案縮短0∶15，減小16%;“東風(fēng)”方案比“西風(fēng)”方案漲潮歷時長0∶58，增大73%。鹽官站“東風(fēng)”方案涌潮高度比“無風(fēng)”方案大0.02 m，僅增大1%;“西風(fēng)”方案涌潮高度比“無風(fēng)”方案小0.21 m，減小9%;“東風(fēng)”方案比“西風(fēng)”方案涌潮高度增大10%。

一般而言，因上游吹程遠，風(fēng)的作用增大，以閘口站為例，“東風(fēng)”方案潮差比“無風(fēng)”方案大1.20 m，增大59%;“西風(fēng)”方案潮差比“無風(fēng)”方案小0.95 m，減小47%;“東風(fēng)”方案比“西風(fēng)”方案潮差大2.15 m，增大197%?！皷|風(fēng)”方案漲潮歷時比“無風(fēng)”方案延長0∶01，增大1%;“西風(fēng)”方案漲潮歷時比“無風(fēng)”方案縮短1∶21，減小96%;“東風(fēng)”方案比“西風(fēng)”方案漲潮歷時長1∶22，增大2 800%?！皷|風(fēng)”方案涌潮高度比“無風(fēng)”方案大0.50 m，增大26%;“西風(fēng)”方案涌潮高度比“無風(fēng)”方案小0.83 m，減小43%;“東風(fēng)”方案比“西風(fēng)”方案涌潮高度增大122%。

表3 各風(fēng)向下涌潮特征值Tab.3 Characteristic values of tidal bore under different wind directions

在風(fēng)速15 m/s 條件下，“西風(fēng)”方案沿程涌潮傳播速度比“無風(fēng)”方案減小6% ～13%，“東風(fēng)”方案沿程涌潮傳播速度增加5% ～9%。為檢驗這一結(jié)果的合理性，對實際發(fā)生的臺風(fēng)期及臺風(fēng)前的情況進行統(tǒng)計對比，如表4。錢塘江臺風(fēng)期間風(fēng)向為東偏北，實測資料表明，臺風(fēng)期涌潮傳播速度普遍快于臺風(fēng)前，由于實際臺風(fēng)的風(fēng)速風(fēng)向均隨時間而變，變化復(fù)雜，因此臺風(fēng)期涌潮傳播速度的增加幅度各次臺風(fēng)相差較大。但在定性上，概化模型計算得到的涌潮傳播規(guī)律與實際一致。另外，在臺風(fēng)期間，經(jīng)常形成較高的高潮位、大潮差和強涌潮，計算結(jié)果也解釋了這一現(xiàn)象。就涌潮傳播速度而言，概化模型和實際情況相比，主要有兩點不同。一是河口進行對稱概化后，潮流上溯的阻力減少，潮波上溯速度變快，動力變強;二是概化的模型在橫斷面上沒有深槽和淺灘的變化，阻力減少，潮波上溯速度變快。如表2 和表4，概化模型中，流量500 m3/s 時，鹽官到倉前的涌潮傳播速度為6.59 m/s;實際臺風(fēng)作用下，表4 所列的實際江道中發(fā)生的涌潮最大傳播速度為5 m/s。

表4 臺風(fēng)前后潮波傳播速度Tab.4 Propagation velocity of tidal bore before or after typhoon

5.2 不同風(fēng)速下涌潮特征變化規(guī)律

加風(fēng)時間為鹽官低潮位前2.5 h，上游徑流500 m3/s 下，計算了無風(fēng)、東風(fēng)風(fēng)速7.5 m/s，東風(fēng)風(fēng)速15 m/s三種工況。加風(fēng)后潮汐特征見表5 和圖9。

由圖9 和表5 可見，無論是潮差、涌潮高度還是漲潮歷時和涌潮傳播速度，均是“無風(fēng)”＜“東風(fēng)風(fēng)速7.5 m/s”＜“東風(fēng)風(fēng)速15 m/s”。風(fēng)速在7.5 m/s 以下，風(fēng)速對以上涌潮特征值的影響相對小一些，風(fēng)速在15 m/s以上，對涌潮特征值的影響較大。如鹽官站“東風(fēng)7.5 m/s”方案潮差比“無風(fēng)”方案大0.14 m，增大3%;“東風(fēng)15 m/s”方案潮差比“無風(fēng)”方案大0.83 m，增大17%。鹽官站“東風(fēng)7.5 m/s”方案漲潮歷時比“無風(fēng)”方案延長0∶09，增大10%;“東風(fēng)15 m/s”方案漲潮歷時比“無風(fēng)”方案延長0:43，增大46%。鹽官站“東風(fēng)7.5 m/s”方案涌潮高度比“無風(fēng)”方案大0.01 m，“東風(fēng)15 m/s”方案涌潮高度比“無風(fēng)”方案大0.02 m。

“東風(fēng)7.5 m/s”方案沿程涌潮傳播速度比“無風(fēng)”方案增加0.7% ～4%，“東風(fēng)15 m/s”方案沿程涌潮傳播速度增加6% ～9%。

圖8 各風(fēng)向下潮波傳播速度Fig.8 Propagation velocity of tidal bore under different wind directions

圖9 各風(fēng)速下潮波傳播速度Fig.9 Propagation velocity of tidal bore under different wind speeds

表5 各風(fēng)速下涌潮特征值Tab.5 Characteristic values of tidal bore under different wind speed

6 結(jié) 語

以錢塘江河口平面尺度和河床地形為基礎(chǔ)，建立了概化河口的二維水動力數(shù)學(xué)模型，通過數(shù)值試驗探討了徑流、風(fēng)向?qū)τ砍钡挠绊憽Ｖ饕Y(jié)論如下:

1)徑流量對涌潮的影響。上游流量越小，沿程各站潮差越大。當(dāng)流量小于1 000 m3/s 時，漲潮歷時的變化差異不大;當(dāng)流量繼續(xù)增加，漲潮歷時一般隨著流量的增加而減小。徑流量對涌潮高度的影響比較復(fù)雜，不同河段存在涌潮高度最大值的相應(yīng)臨界徑流量，越往下游，臨界徑流量越大，在下游鹽官、倉前站，涌潮高度在5 000 m3/s 時達到最大，在上游七堡、閘口、聞家堰，涌潮高度在1 000 m3/s 達到最大。涌潮傳播速度與潮前、潮后水深、涌潮高度和潮前落潮流速有關(guān)，大多情況下涌潮傳播速度隨著徑流量的增大而減小。

2)風(fēng)況對涌潮的影響。計算結(jié)果表明，潮差、涌潮高度、漲潮歷時和涌潮傳播速度，均是“西風(fēng)(逆風(fēng))”＜“無風(fēng)”＜“東風(fēng)(順風(fēng))”。根據(jù)臺風(fēng)期間調(diào)研和實測資料分析，計算結(jié)果與實際情況符合。順風(fēng)條件下，風(fēng)速越大，涌潮高度越大。

3)涌潮傳播規(guī)律研究成果是基于定床下的概化模型，在各流量下，地形對流量有不同的響應(yīng)，如洪水期上游江道沖深，下游淤積，致涌潮傳播規(guī)律也會有所差異，這個影響有待于下一階段進行研究。

［1］Chanson H.Current knowledge in tidal bores and their environmental，ecological and cultural impacts［J］.Environmental Fluid Mechanics，2011，11:77-98.

［2］韓曾萃，戴澤蘅，李光炳.錢塘江河口治理開發(fā)［M］.北京:中國水利水電出版社，2003.

［3］潘存鴻，魯海燕，曾 劍.錢塘江涌潮特性及其數(shù)值模擬［J］.水利水運工程學(xué)報，2008(2):1-9

［4］林炳堯，黃世昌，毛獻忠，等.錢塘江河口潮波變化過程［J］.水動力學(xué)研究與進展:A 輯，2002，17 (6):665-675.

［5］趙雪華.幾項主要因素對錢塘江涌潮形成和發(fā)展的影響［R］.杭州:浙江省河口海岸研究所，1982.

［6］曾 劍，熊紹隆，潘存鴻，等.運用神經(jīng)網(wǎng)絡(luò)理論研究錢塘江涌潮的影響因素［J］.長江科學(xué)院院報，2006，23(5):14-16，20。

［7］潘存鴻.淺水間斷流動數(shù)值模擬研究進展［J］.水利水電科技進展，2010，30(5):78-84.

［8］潘存鴻，徐 昆.三角形網(wǎng)格下求解二維淺水方程的KFVS 格式［J］. 水利學(xué)報，2006，37(7):858-864.

［9］潘存鴻，于普兵，魯海燕.淺水動邊界的干底Riemann 解模擬［J］.水動力學(xué)研究與進展:A 輯，2009，24(3):305-312.

［10］潘存鴻，魯海燕.二維淺水間斷流數(shù)值模型在涌潮模擬中的應(yīng)用［J］.浙江大學(xué)學(xué)報:工學(xué)版，2009，43(11):2107-2113.