許雪峰,羊天柱,孫志林,施偉勇,聶 源

(1.浙江大學 建筑工程學院,浙江 杭州 310012;2.國家海洋局 第二海洋研究所 工程海洋學研究中心,浙江 杭州 310012)

中國東部沿海灘涂分布廣闊,如蘇北淺灘、上海蘆潮港東部淺灘、杭州灣慈溪淺灘、溫州淺灘、溫州南部蒼南沿海潮灘以及東部沿海港灣內(nèi)均有大量潮灘。它們少則十幾平方公里,多則上百平方公里。此外,中國東部沿海潮差較大,除舟山群島以外,一般潮差達到4 m以上,于是落潮時大片潮灘露灘,漲潮時漫灘。在這一類海岸的潮流數(shù)值模擬中,漫灘潮流的模擬是重點,但也是難點,主要原因是：(1)由于觀測方法和技術(shù)的原因,對漫灘潮流觀測難度較大,因此難以對漫灘潮流的數(shù)值模擬進行驗證。但是隨著觀測技術(shù)發(fā)展,漫灘潮流的觀測成為可能,如2004年5月國家海洋局第二海洋研究所在杭州灣慈溪淺灘的東南側(cè)進行了漫灘和露灘流過程的觀測,為本研究提供了基礎資料。(2)目前對漫灘動力邊界條件的處理一般采用干濕網(wǎng)格法[1],邊界上水位一般取臨近網(wǎng)格節(jié)點的平均值,這可能導致邊界上計算水深失真。本研究介紹一種漫灘邊界上向上取水位值的方法,使得邊界上地形突變時模型穩(wěn)定性增強。

本研究基于實測漫灘潮流資料,建立平面二維數(shù)值模型[2],在模型的漫灘邊界上采用不同的水位處理方法,并驗證方法的適宜性,最后模擬漫灘潮流的全過程。

1 研究區(qū)域的實測漫灘潮流

本文選擇杭州灣海域為大范圍研究區(qū)域,選擇慈溪淺灘為漫灘流重點研究區(qū)域。杭州灣概況如圖1所示,圖中標注了本文用來驗證和參考資料的各水文站位置：共有 9個潮位站;6個潮流站,其中主槽有4個潮流站,漫灘區(qū)有2個潮流站。

圖1 杭州灣的概況、以及各水文測站的位置Fig.1 Overview of the Hangzhou Bay,and locations of hydrological stations

A(30°17′54.5"N,121°25′58.9"E),B(30°18′23.9"N,121°26′21"E)為兩個漫灘流觀測站,位于慈溪淺灘上(見圖1)。由于低水位時露灘,高水位時漫灘,故A,B點布放測站時采用底座式觀測方式,儀器采用1.2MK的ADCP,從2005年5月20～23日連續(xù)觀測3 d,觀測間隔10 min,觀測層次設0.5 m一層。底座式ADCP的固定方法和工作方式見圖2。

圖2 ADCP固定及工作方式示意圖Fig.2 Scheme of the fixation and mechanism of ADCP

取72 h有水時段的平均流速,并取出兩個潮周期的流速資料,分別統(tǒng)計漲落潮流速,流向。統(tǒng)計KP50.5,KP49.5兩個處于灘涂上的潮流測站最大流速、平均流速統(tǒng)計見表1、表2。

表1 灘涂上實測最大潮流Tab.1 The measured maximum flow on tidal flat

表2 灘涂上實測平均流速Tab.2 The measured average flow on tidal flat

取KP49.5,KP50.5兩點的垂向平均的潮流資料,作出流速、流向和水深的過程曲線見圖3。

從測量數(shù)據(jù),本文得到了漫灘區(qū)潮流的以下幾個特性：(1)漫灘流的流速明顯小于前沿深水區(qū)域的潮流。從漲落潮流速來看,落潮流速略占優(yōu)勢。A測站最大漲潮流流速為116cm/s,最大落潮流流速145cm/s,平均漲潮流流速為 47.4cm/s,平均落潮流流速為51.7cm/s;B測站最大漲潮流流速為100cm/s,最大落潮流流速99cm/s,平均漲潮流流速為32.7cm/s,平均落潮流流速為 38.4cm/s。而與此同時,潮灘前沿的主槽上的 KP40站的最大漲潮流速為 391cm/s,最大落潮流速為 296cm/s。(2)漫灘區(qū)潮流漫灘過程(漲潮)和露灘過程(落潮)的時間并不對稱,但主槽上的漲落潮時間較為對稱。本次現(xiàn)場觀測得到以下特征數(shù)據(jù)：A點漲落潮時間比為1.38∶1,B點漲落潮時間比為 1.30∶1;主槽上 KP40站的漲落潮時間比為1.08∶1。(3)潮灘上最大漲落潮流出現(xiàn)的時刻與主槽上最大漲落潮流出現(xiàn)的時刻不一致。潮灘上A,B站最大漲潮流速一般出現(xiàn)在最高水位時刻之前的1.5～2.5 h,最大落潮流速一般出現(xiàn)在最高水位時刻之后的0.5～1 h,最大漲落潮流速出現(xiàn)的時刻間隔僅2～3 h;而主槽上的 KP40站漲落潮的最大流速一般都出現(xiàn)在中水位時刻,最大漲落潮流速出現(xiàn)的時刻間隔一般在6 h左右。因此在漲落潮過程中有可能出現(xiàn)主槽上的潮流方向和潮灘上的漫灘流方向相反的現(xiàn)象。(4)從漲落潮流向來看,研究區(qū)域的漫灘潮流的流向與潮灘的岸線并不垂直,漲潮流向表現(xiàn)為攏流特征,落潮流向表現(xiàn)為開流特征。KP49.5站漲潮流主流向一般在 285°～305°之間,在 115°～145°之間;KP50.5站漲潮流主流向一般在 275°～295°之間,落潮流向約在95°～135°左右。

圖3 KP49.5,KP50.5流速、流向、水位過程曲線Fig.3 The velocity,flow,and depth curves at observations KP 49.5 and 50.5

2 潮灘海域的數(shù)值模擬方法

2.1 模型方程及初始邊界條件

在正交曲線網(wǎng)格下采用 N-S控制方程[1],并引用傳統(tǒng)的 ADI法[2]對方程進行離散求解,初始條件取水位、流速均為零[3]。

模型開邊界采用水位邊界,采用歷史資料的調(diào)和分析得到[4]。作者搜集了杭州灣附近海域多個潮位站的資料,其中選取了鎮(zhèn)海、金塘瀝港、馬目、大漁山、小洋山、蘆潮港潮位站(圖1)的潮位資料作為模型邊界條件輸入,而乍浦、金山潮位站的潮位資料用于模型的驗證。

2.2 漫灘邊界處理方法

在固定閉邊界上,取法向流速為 0,切線方向上的邊壁流速取對數(shù)分布[5]。

漫灘邊界采用干濕網(wǎng)格法處理[5-6]。其中在漫灘邊界上水位ζ又分別采用了傳統(tǒng)處理方法和本文介紹的“向上取水位值”的方法。傳統(tǒng)方法對漫灘邊界上水位的處理一般取臨近網(wǎng)格節(jié)點的水位平均值,因此在灘涂上的地形突變區(qū)域,當插值所得的水位低于海底高程時,總水深將為負值,邊界附近節(jié)點的計算將無法進行,但此時灘涂上節(jié)點的水位卻仍然是高于外側(cè),是明顯不合理的,如圖4所示。在模型計算中,若多個連續(xù)的計算節(jié)點同時出現(xiàn)上述的情況,可能導致灘涂上水位和流速計算值不合理,最終可能導致模型計算不穩(wěn)定并終止。

圖4 在潮灘上不同水位處理方法的效果Fig.4 The effects of different water treatment methods on tidal flat

“向上取水位值”的方法：由于上文提到的常用的方法導致的模型不穩(wěn)定均發(fā)生在露灘過程中(圖4),因此首先要判斷該計算節(jié)點是否為露灘過程的干濕邊界點,模型中先用水深來判斷是否可能為正在露灘的計算點,當時可能正在露灘,再采用計算節(jié)點的水位、相鄰點的水位和流速來綜合判別是否需要對該點 “向上取水位值”以保證計算水深值不失真,如圖4所示。具體判斷方法如式(1),水位ζ和水深的處理方法如式(1)～(5),這里Hmin為模型中預先設定的極小水深,本文經(jīng)過多次試算認為Hmin=0.1 m較為合適。

在曲線網(wǎng)格ξ方向上,當

在正交網(wǎng)格中,計算節(jié)點的總水深Hm,n由η方向的水深dη和水位ξζ相加得到最終參與計算的節(jié)點水深值：

在普通的計算節(jié)點上：

在曲線網(wǎng)格η方向上,的取值與ξ方向上類似。

3 模型實例

本文選用杭州灣慈溪潮灘為研究區(qū)域,采用曲線坐標網(wǎng)格對計算域進行剖分,網(wǎng)格范圍東西跨度124 km,南北最大跨度 96 km。計算域內(nèi)剖分成500×600,總共有 300 000個網(wǎng)格,最大的網(wǎng)格邊長取 1 000 m左右,主要研究區(qū)域附近的網(wǎng)格邊長為150 m左右,計算時間步長為1 min,計算水深取自1：10 000的杭州灣海圖,模型計算網(wǎng)格見圖5。

圖5 計算網(wǎng)格以及邊界位置Fig.5 Grid and boundary location

3.1 模型驗證

潮位采用金山潮位站、乍浦站 11d的實測資料進行驗證;漫灘潮流采用A,B兩站的實測漫灘潮流資料進行驗證,以此證明模型的可靠性。

3.1.1 潮位驗證

潮位驗證誤差見表3。由表可見,大、小潮期間實測潮位與模擬計算的潮位之間擬合得較好,最大潮差的模擬誤差一般在10～15cm左右,僅小潮個別誤差在30cm左右,大潮模擬潮差相對誤差均在3%以內(nèi),模擬和實測總體擬合較好。

表3 乍浦、金山潮位驗證誤差Tab.3 Tide validation error of Zhapu and Jinshan

3.1.2 漫灘上潮流驗證

漫灘潮流驗證誤差見表4。由于該研究區(qū)域為寬廣的淺灘,潮流漫灘、露灘過程時間較長。從模擬結(jié)果來看,漫灘流的各項特征得到了很好的模擬,如落潮強于漲潮,漫灘和露灘過程的時間不對稱性等。從驗證誤差來看,最大流速誤差一般在 10cm/s左右,而平均流速的驗證誤差一般都為 5cm/s左右,流向誤差也均控制在15°以內(nèi)?？偟膩碚f,漫灘潮流的模擬驗證情況較好。

表4 漫灘潮流驗證誤差Tab.4 Validation error of floodplain flow

實測站點KP49.5和KP50.5的漫灘潮流驗證過程曲線見圖6。

圖6 KP49.5,KP50.5潮流站流速、流向驗證Fig.6 Validation error of floodplain flow of KP49.5 and 50.5

3.2 漫灘區(qū)域的潮流場

為了更清晰地演示漫灘潮流的變化過程,本文給出了漫灘過程中的1 h間隔的潮流場,見圖7。

圖7 研究區(qū)漫灘潮流全過程流場圖Fig.7 The whole process of floodplain flow of the study area

由圖7可見：(1)漫灘區(qū)域的漲潮時間要略早于主航道上,當主航道上仍有在落潮或轉(zhuǎn)流的時候,漫灘上的流已經(jīng)呈現(xiàn)漲潮的方向,這也與實測資料的分析結(jié)果一致。(2)漫灘區(qū)域的落潮時間也要略早于主航道上,當主航道上仍在漲潮或轉(zhuǎn)流的時候,漫灘流已經(jīng)呈現(xiàn)落潮的方向。(3)灘涂上的漫灘流在空間上可分為兩個區(qū)域,一為干濕相交的水陸分界面區(qū)域,二為分界面向水一側(cè)的淹水區(qū)域。在漲潮過程中,漫灘岸界面的漲潮流流向偏向岸線的法向,淹水區(qū)域漫灘潮流將偏向主航道方向(岸線切線方向);落潮流與漲潮流情況類似。

4 結(jié)論

通過實測資料分析和潮流數(shù)值模擬及驗證,可以看出,采用漫灘邊界處向上取水位值是漫灘潮流模擬中一種簡便可行的邊界處理方法。此方法有效地解決了灘涂上地形突變導致的模型不收斂問題。從模擬結(jié)果和實測資料的驗證分析來看,誤差較小,模擬能夠真實地反應漫灘潮流的各項特征。

[1]張越美,孫英蘭.渤海灣三維變動邊界潮流數(shù)值模擬[J].青島海洋大學學報(自然科學版),2002,32(3)：337-344.

[2]華秀箐,呂玉麟.二維淺水域潮流數(shù)值模擬的ADI-Quick格式[J].水動力學研究與進展,A輯,1996,11(1)：77-92.

[3]湯立群.正交曲線坐標下河口二維潮流過程計算[J].水動力學研究與進展,2003,18：16-25.

[4]陳宗鏞.潮汐學[M].北京：科學出版社,1980.

[5]毛獻忠,潘存鴻.移動邊界淺水問題的數(shù)值研究[J].水動力學研究與進展,2002,17(4)：33-43.

[6]史峰巖,孫文心.極坐標變換變邊界模型及其應用[J].海洋與湖沼,1995,26(4)：369-376.