胡 鵬，姬奧飛，陶俊余

(浙江大學(xué) 海洋學(xué)院，浙江 舟山 316021)

淺水方程被廣泛用于模擬潮流過程[1-5]。在20世紀七八十年代，普萊斯曼(Preissmann)隱格式[6]因結(jié)構(gòu)簡單、無條件穩(wěn)定以及對時間步長沒有限制等優(yōu)點，被廣泛用于求解淺水方程[7]。夏云峰[8]采用SIMPLE、SIMPLEC和SIMPLER算法對動量方程和連續(xù)方程耦合求解。近年，基于壓力分裂模式和θ半隱算法，并采用歐拉-拉格朗日方法(ELM)求解對流項的淺水模型，因其“無條件穩(wěn)定”的特征而得到了關(guān)注[9]。淺水方程的解可能包含間斷/激波[10-13]，表現(xiàn)為水跌或水躍、風暴潮、涌潮等。為捕捉這些現(xiàn)象，基于黎曼近似解通量的有限體積法被廣泛應(yīng)用[14-15]，此類模型時間步長受到約束，如CFL穩(wěn)定條件[16]。當網(wǎng)格尺度或水流不均勻時，網(wǎng)格之間的最大允許時間步長可能相差很大。為程序編寫方便，這類模型往往采用整體最小時間步長，將單元時間步長取為所有網(wǎng)格允許最大時間步長的全局最小值，這簡化了程序編寫，但限制了計算效率。為提高計算效率，局部時間步長(LTS)方法受到關(guān)注。相對于傳統(tǒng)方法，LTS法使每個網(wǎng)格采用滿足穩(wěn)定性的盡可能大的時間步長[17-19]。然而，當最大和最小時間步長相差倍數(shù)太大時，二維淺水模型的LTS法可能存在不穩(wěn)定性[19]。胡鵬等[16, 20]通過限制激波/間斷處的時間步長空間梯度，得到了穩(wěn)定的局部時間步長方法。在此基礎(chǔ)上，建立高性能的平面二維潮流數(shù)值模型，以長江口、杭州灣和渤海、黃海、東海的潮流為例，驗證其計算效率和精度。

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

1.1 基于有限體積法的淺水模型

平面二維淺水方程可寫成向量形式如下：

(1)

(2)

非結(jié)構(gòu)三角形網(wǎng)格示意如圖1所示。

圖1 非結(jié)構(gòu)三角形網(wǎng)格示意Fig. 1 Sketches of the unstructured triangular meshes

圖1(a)為計算區(qū)域內(nèi)部某單元(編號為i，i=1,2,3,……,Nc；Nc為總單元數(shù))以及三個相鄰單元(編號分別記為：i1、i2、i3)；圖1(b)為某個界面(編號為j，j=1,2,3,……,Nf；Nf為總界面數(shù))及其左右兩側(cè)的單元(編號分別記為：jL、jR)。

針對單元i，采用有限體積法對控制方程(1)進行離散可得[21]：

(3)

Enj=FHLLC(UjL,UjR)

(4)

其中，UjL和UjR分別為界面j兩側(cè)的黎曼狀態(tài)；FHLLC為HLLC黎曼算子。調(diào)用HLLC黎曼算子之前，采用非負水深重構(gòu)界面兩側(cè)的黎曼狀態(tài)[24-25]。HLLC黎曼算子的具體表達式可參考文獻[21]。

1.2 LTS方法及其實現(xiàn)過程

首先通過下式計算各單元允許的最大時間步長：

(5)

式中：Δtami為第i單元的最大允許時間步長；Cr為克朗數(shù)，設(shè)定為0.9；uij，vij為第i單元第j界面法向水流流速；hi為第i個單元的水深；Rij為單元中心到第j界面的距離；hthr=10-6m為臨界水深。其次，計算整體最小時間步長Δtmin和各單元時間分級指數(shù)mi：

(6)

(7)

(8)

1) 計算界面數(shù)值通量時，如果特征值(np-1)/2mfj為整數(shù)，則計算，否則不更新。

(9)

2 LTS方法對模擬結(jié)果的影響

采用L(f)量化LTS方法與傳統(tǒng)整體最小時間步長方法之間的相對誤差，這里L(f)描述的相對誤差是指采用LTS法計算的某一變量值與采用整體最小時間步長方法計算出的這一變量值之間的相對差值，并對計算區(qū)域中所有節(jié)點的相對差值進行統(tǒng)計，最終得到一個面上的統(tǒng)計值。采用ε量化全局水體質(zhì)量相對誤差。L(f)及ε計算式為：

(10)

(11)

式中：參數(shù)f表示水動力變量，如水位、流速等；N為網(wǎng)格節(jié)點總數(shù)；fLTS和fGMi分別表示LTS方法和傳統(tǒng)整體最小時間步長方法的計算結(jié)果；V(t1)和V(t2)表示計算區(qū)域在t1和t2兩個時刻包含的水體總體積；δV表示在這兩個時刻之間從邊界凈流入或流出的水體體積。

計算區(qū)域如圖2所示。三角形網(wǎng)格總數(shù)為117 323，最小邊長約180 m，最大邊長約38 000 m。地形采用2016年2月實測地形資料(徐六涇—口門附近)和杭州灣、渤海、黃海、東海海圖。長江徑流邊界取在三江營，采用大通流量過程作為邊界條件。長江口潮流界在江陰以下，潮區(qū)界位于銅陵和蕪湖間[26]。因此，三江營不受潮流影響。錢塘江給流量800 m3/s。外海開邊界位于大陸架內(nèi)側(cè)，為水位驅(qū)動，考慮M2等13個分潮，通過海潮模型TPXO計算。初始水位、流速設(shè)為0；各單元的糙率采用公式0.01+0.01/h計算[27]。

圖2 計算區(qū)域、非結(jié)構(gòu)網(wǎng)格和實測數(shù)據(jù)站點位置Fig. 2 Computational domain, unstructured grids and stations for measuring data

模擬2016年7月10日至7月13日的潮流過程，考慮muser取值從1到7共7個LTS方法工況，與傳統(tǒng)整體最小時間步長方法對比(即取muser=0)。表1是兩種方法計算效率和相對誤差的對比。由表1可知：1) LTS方法能大幅度提高模型計算效率。模型計算耗時隨著muser的增加而降低。傳統(tǒng)方法計算耗時12.75 h，而LTS法(取muser=7)僅耗時1.77 h，LTS方法的計算效率提高約7.2倍。2) LTS方法帶來的相對誤差可以忽略，遠小于傳統(tǒng)方法與實測數(shù)據(jù)之間的相對誤差。雖然LTS方法帶來的相對誤差隨著muser增大略有增加，但最大僅為12 mm(水位)、0.66 mm/s(流速)，遠小于傳統(tǒng)方法與實測數(shù)據(jù)之間的相對誤差：179 mm(水位)，186 mm/s(流速)。3) 基于LTS方法的潮流模型具有滿意的守恒性。采用式(11)計算全局水體質(zhì)量相對誤差ε，其數(shù)值越接近0說明計算格式守恒性越高。結(jié)果顯示，基于LTS法的有限體積法(muser=7)的全局水體質(zhì)量相對誤差ε僅為10-12，表明模型的守恒性令人滿意。

表1 LTS方法與傳統(tǒng)方法之間的對比Tab. 1 Comparison between LTS approach and traditional method

3 模式在長江口、杭州灣及渤海、黃海、東海的應(yīng)用

3.1 渤海、黃海、東海潮汐特征

圖3為計算所得M2分潮在渤海、黃海、東海的同潮圖，包括潮差分布(圖3(a))和遲角分布(圖3(b))。從圖可以看出，潮波從外海邊界傳入后，向西北和西南兩個方向傳播。向西北方向傳播的潮波，首先在渤海、黃海分別形成兩個逆時針旋轉(zhuǎn)的潮波系統(tǒng)，兩個中心點分別在(34°50′N，120°35′E)和(37°59′N，122°50′E)附近；之后潮波繼續(xù)向北，經(jīng)渤海海峽后再次分成兩個方向，部分沿渤海海峽方向向西傳播，另一部分右轉(zhuǎn)向北傳播，再次形成兩個逆時針旋轉(zhuǎn)的渤海半日潮波系統(tǒng)，其中心點分別在(38°26′N，119°5′E)和(40°15′N，120°48′ E)附近。上述結(jié)果與張衡等[28]所述的潮波傳播特征基本一致。如圖3(a)，M2分潮最大振幅約為3.14 m，位置在朝鮮半島西部的江華灣附近(圖3(a)中黑色圓點所示)。朱學(xué)明等[29]基于FVCOM海洋數(shù)值模式模擬了渤海、黃海、東海的潮汐、潮流，認為江華灣頂M2分潮最大振幅達3.2 m，支持了本文模擬結(jié)果。黃海區(qū)域和渤海區(qū)域分別存在兩個振幅接近于0的點，即無潮點，如圖3(a)黑色三角形所示，與沈育疆[30]所得無潮點位置(圖3(a)黑色五角星)接近。對比圖3(a)和圖3(b)還可看出，無潮點周圍對應(yīng)一個逆時針旋轉(zhuǎn)的潮波系統(tǒng)；無潮點周圍分潮振幅等值線與遲角等值線大致垂直。

圖3 M2分潮同潮圖Fig. 3 The distribution of coamplitude and cophase of the M2 tide constituent

3.2 潮位與潮流驗證

根據(jù)可得實測數(shù)據(jù)，驗證長江口四個測站(石洞口、雞骨礁、南槽東、北槽中)的潮位(圖4，2016年7月10日至8月10日)，兩個測站的流速和流向(圖5，2016年7月21日至月22日)；杭州灣四個測站(洋山港、北侖、岱山、綠華)的潮位(圖6，2015年6月29日至7月9日)。測站位置如圖2所示。從圖4～圖6可以看出，潮位、潮流(流速和流向)的模擬值與實測值吻合良好。

圖4 長江口計算水位與觀測值比較Fig. 4 Comparison of calculated and observed water levels in the Yangtze River estuary

圖5 垂向平均流速和流向驗證Fig. 5 Validation of vertical average current speed and direction

圖6 杭州灣計算水位與觀測值比較Fig. 6 Comparison of calculated and observed water levels in Hangzhou Bay

采用均方根誤差RMSE、相關(guān)系數(shù)CC和技術(shù)評分SS等參數(shù)進一步量化誤差，計算式如下

(12)

(13)

(14)

誤差統(tǒng)計結(jié)果如表2所示。潮位計算值與實測值之間的均方根誤差范圍在0.179 m到0.362 m之間(綠華站最小，石洞口站最大)。考慮到最大潮差約為4.5 m(洋山港)，這樣的潮位誤差可接受。流速、流向的均方根誤差范圍在0.186 m/s到0.287 m/s之間(NGN4S站最小，CS9S站最大)和13°37′到24°4′之間(CS9S站最小，NGN4S站最大)。考慮到最大流速約為2.2 m/s(CS9S站)，流向最大變化范圍約為180°，模型計算的流速、流向誤差也可接受。相關(guān)系數(shù)CC和技術(shù)評分SS，分別在0.9以上和0.76以上，表明模型計算值與實測值吻合程度很好[31]。

表2 誤差分析Tab. 2 Error analysis

4 結(jié) 語

建立了基于LTS方法的平面二維潮流數(shù)值模型。模型采用非結(jié)構(gòu)三角形網(wǎng)格，對局部區(qū)域加密，通過LTS方法提高模型計算效率。對比表明，LTS方法和傳統(tǒng)整體最小時間步長方法兩者之間的水位、流速、流向相對誤差均較小，但LTS法可使模型計算效率大幅度提高(本文工況計算效率提高了約7.2倍)。模型全局水體質(zhì)量相對誤差，在muser=7時僅為10-12，表明模型計算格式的守恒性也較高。最后，采用所建立的高效率平面二維模型，成功模擬了長江口、杭州灣、渤海、黃海、東海的潮流過程，計算所得潮流傳播特征、振幅分布、無潮點均與前人結(jié)果相符，與長江口、杭州灣十余個站點實測數(shù)據(jù)吻合良好。需要說明的是：LTS方法提高模型效率的具體倍數(shù)會隨網(wǎng)格和流速非均勻性變化，在非常理想的情況下，即網(wǎng)格和流速都很均勻時，LTS的加速效果將減弱。實際上，網(wǎng)格局部加密非常普遍，實際水流往往也是非均勻的，文中模型在大多具體實例中的加速效果將是顯著的。