基于局部時(shí)間步長(zhǎng)方法的潮流數(shù)值模型研究及應(yīng)用

胡 鵬，姬奧飛，陶俊余

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

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

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

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

(3)

Enj=FHLLC(UjL,UjR)

(4)

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

1.2 LTS方法及其實(shí)現(xiàn)過(guò)程

首先通過(guò)下式計(jì)算各單元允許的最大時(shí)間步長(zhǎng)：

(5)

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

(6)

(7)

(8)

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

(9)

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

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

(10)

(11)

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

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

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

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

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

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

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

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

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

3.2 潮位與潮流驗(yàn)證

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

圖4 長(zhǎng)江口計(jì)算水位與觀測(cè)值比較Fig. 4 Comparison of calculated and observed water levels in the Yangtze River estuary

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

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

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

(12)

(13)

(14)

誤差統(tǒng)計(jì)結(jié)果如表2所示。潮位計(jì)算值與實(shí)測(cè)值之間的均方根誤差范圍在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°，模型計(jì)算的流速、流向誤差也可接受。相關(guān)系數(shù)CC和技術(shù)評(píng)分SS，分別在0.9以上和0.76以上，表明模型計(jì)算值與實(shí)測(cè)值吻合程度很好[31]。

表2 誤差分析Tab. 2 Error analysis

4 結(jié) 語(yǔ)

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