孟江山，路川藤，羅小峰，丁 偉

（1.中交廣州航道局有限公司，廣東 廣州 510000；2.南京水利科學(xué)研究院，江蘇 南京 210029；3.河海大學(xué)，江蘇 南京 210024）

1 研究背景

廈門周邊海域島嶼眾多，潮汐汊道交錯(cuò)，水流運(yùn)動(dòng)較為復(fù)雜。廈門海域潮汐屬于正規(guī)半日潮［1］，潮波主要受來自臺(tái)灣海峽的潮波傳播影響，水動(dòng)力強(qiáng)度相對較弱，且主要由潮汐作用引起［2］。由于廈門附近海域近年來圍墾較多，因此眾多學(xué)者研究廈門海域潮流泥沙運(yùn)動(dòng)較多，而系統(tǒng)性的潮波傳播研究成果相對較少，如謝森揚(yáng)［3］、陸榮華［4］均系統(tǒng)研究了廈門歷史圍填海的水動(dòng)力、水環(huán)境問題，許婷［5］、陳純等［6］還針對具體工程問題研究了局部海域水沙環(huán)境問題。潮波傳播理論體系目前已相對成熟［7］，研究方法以實(shí)測資料分析和數(shù)學(xué)模型模擬兩種手段為主，實(shí)測資料分析以調(diào)和分析［8-9］、神經(jīng)網(wǎng)絡(luò)［10］等方法為主，數(shù)學(xué)模型模擬可以彌補(bǔ)實(shí)測資料時(shí)段短及區(qū)域覆蓋面局限等缺點(diǎn)，因此被廣泛應(yīng)用［11-12］，近年來，隨著計(jì)算機(jī)科學(xué)的發(fā)展，數(shù)學(xué)模型計(jì)算正朝著高分辨率、高計(jì)算效率方向發(fā)展［13］，特別是GPU（Graphics Processing Unit）并行計(jì)算，大幅度節(jié)省了數(shù)學(xué)模型的計(jì)算時(shí)間，使得數(shù)學(xué)模型可以快速實(shí)現(xiàn)精細(xì)化、實(shí)時(shí)化模擬預(yù)報(bào)，如龔佳輝等［14］在城市雨洪模擬中，認(rèn)為GPU并行可提速23.88～158.72倍。本文在前人研究的基礎(chǔ)上，采用GPU并行算法模擬廈門周邊海域潮波傳播，深入研究廈門附近海域潮波傳播特征，為相關(guān)工程建設(shè)和規(guī)劃提供技術(shù)支撐。

2 GPU并行算法

2.1 離散方法 數(shù)學(xué)模型采用具有自主知識(shí)產(chǎn)權(quán)的數(shù)值模擬系統(tǒng)CJK3D-WEM計(jì)算，該系統(tǒng)于2014年取得國家軟件著作權(quán)登記，適用于江河湖泊、河口海岸等涉水工程中的水動(dòng)力、泥沙、水質(zhì)、溫排和溢油模擬預(yù)測研究。

對于水平尺度遠(yuǎn)大于垂直尺度的河口海岸地區(qū)，通過一定的假設(shè)條件（均質(zhì)不可壓、靜壓、剛蓋、Boussinesq等假定），沿水深方向積分取平均，可得到平面二維水動(dòng)力方程：

式中：H為總水深，m；z為水位，m；u、v流速矢量V沿x、y方向的速度分量，m/s；t為時(shí)間，s；f為科氏力，s-1；g為重力加速度，m/s2；Nx、Ny為x、y向水流紊動(dòng)黏性系數(shù)，m2/s；C為謝才系數(shù)。

控制方程式（1）向量形式可表示為：

其中：

其中:

紊動(dòng)擴(kuò)散項(xiàng)：

其中:

源項(xiàng)M表示為：

式中：Mox、Moy分別為x、y方向的河床底部高程變化；Mfx、Mfy分別為x、y方向的底摩擦項(xiàng)。

采用三角形網(wǎng)格作為控制單元，水深布置在網(wǎng)格頂點(diǎn)，其他物理變量配置在每個(gè)單元的中心。將第i號(hào)控制元記為Ωi，在Ωi上對向量式的基本方程組（2）進(jìn)行積分，并利用Green公式將面積分化為線積分，得：

式（3）求解主要分三部分，一為對流項(xiàng)求解，二為紊動(dòng)項(xiàng)求解，三為底坡項(xiàng)處理。對流項(xiàng)的數(shù)值通量采用Roe格式的近似Riemann解，紊動(dòng)項(xiàng)采用單元交界面的平均值，計(jì)算通過該界面紊動(dòng)黏性項(xiàng)的數(shù)值通量，底坡項(xiàng)采用斜底方法處理。

式（3）中對流項(xiàng)可表示為：

(E?n)R和(E?n)L分別為交界面兩側(cè)的數(shù)值通量，可表示為：

其中：

這里λk和ek（k=1，2，3）分別為對流項(xiàng)特征矩陣A的特征值和特征向量，滿足：ΔE?=AΔU ，特征矩陣形式如下：

矩陣中的u?、v?和c?分別采用質(zhì)量加權(quán)平均計(jì)算，表達(dá)如下：

矩陣A的特征值和特征向量均可計(jì)算得出，見文獻(xiàn)［15］。

水平擴(kuò)散項(xiàng)含有二次項(xiàng)，是離散淺水方程的難點(diǎn)之一。本文采用單元交界的平均值計(jì)算通過該界面擴(kuò)散項(xiàng)的數(shù)值通量，這樣處理算法簡單、效率高。式（3）中的水平擴(kuò)散項(xiàng)可表示為：

對于上式的偏導(dǎo)數(shù)可用如下方式求解：

底坡項(xiàng)采用斜底模型求解，記某個(gè)控制體三個(gè)頂點(diǎn)坐標(biāo)為xi，yi，zi（i=1，2，3），三個(gè)頂點(diǎn)按逆時(shí)針方向排序，由這三個(gè)頂點(diǎn)確定的平面方程為：

其中：

其中：當(dāng)Dz≠0時(shí)，

2.2 GPU并行計(jì)算 并行計(jì)算是指同時(shí)使用多種計(jì)算資源解決計(jì)算問題的過程，是提高計(jì)算機(jī)系統(tǒng)計(jì)算速度和處理能力的一種有效手段。目前，流體力學(xué)并行計(jì)算一般基于CPU和GPU兩種手段，CPU單核計(jì)算能力強(qiáng)，但核心數(shù)量較少，從而限制了并行計(jì)算能力，GPU雖單核計(jì)算能力差，但核心數(shù)量多，因此并行計(jì)算效率較CPU強(qiáng)。

GPU并行計(jì)算程序基于CUDA架構(gòu)下的C語言編寫，采用單精度算法。GPU并行計(jì)算流程見圖1，當(dāng)CPU計(jì)算參數(shù)通過cudaMemcpy（dev_input，input，sizeof（int），cudaMemcpyHostToDevice）函數(shù)傳遞給GPU后，GPU通過核函數(shù)kernel實(shí)現(xiàn)大量并行執(zhí)行線程，從而達(dá)到并行計(jì)算的效果，GPU計(jì)算結(jié)果通過cudaMemcpy（result，dev_result，sizeof（int），cudaMemcpyDeviceToHost）函數(shù)反饋給CPU，需要說明的是CPU與GPU之間的數(shù)據(jù)交換耗時(shí)較長，因此實(shí)際計(jì)算過程中應(yīng)盡量降低二者之間的數(shù)據(jù)交換頻率，以過丘流算例為例進(jìn)行說明，計(jì)算采用的CPU型號(hào)為E5-1630V4CPU（8核心），GPU型號(hào)為NVIDIA Ge?force RTX2080（2944核心）。

圖1 GPU并行計(jì)算流程圖

建立了長60 m，寬6 m的數(shù)值水槽，以三角形單元?jiǎng)澐志W(wǎng)格，網(wǎng)格總數(shù)約為20 000個(gè)，空間步長最小約為0.13 m，模型范圍與地形等見圖2（a）和圖2（b），不考慮紊動(dòng)黏性系數(shù)，模擬有激波條件下的過丘流，水槽上游徑流取0.18 m3/s，下游水位-0.07 m，模擬時(shí)間1 h，由于控制方程離散格式為顯式，時(shí)間步長受Courant數(shù)約束，因此時(shí)間步長取值0.003s，需運(yùn)行1 200 000次循環(huán)方可完成模擬。圖2（c）為CPU單核（單精度）與GPU并行（單精度）計(jì)算結(jié)果與解析解的對比，計(jì)算精度一致。圖3為GPU并行計(jì)算過程中，GPU與CPU之間的數(shù)據(jù)交換頻次對耗時(shí)的統(tǒng)計(jì)，當(dāng)采用CPU單核進(jìn)行計(jì)算時(shí)，計(jì)算用時(shí)約36 h，若GPU并行計(jì)算過程中，GPU與CPU之間實(shí)時(shí)（即每步）進(jìn)行數(shù)據(jù)交換，則計(jì)算效率嚴(yán)重降低，耗時(shí)超過72 h，隨著GPU與CPU數(shù)據(jù)交換頻次的降低，計(jì)算時(shí)間呈指數(shù)下降趨勢，當(dāng)GPU每模擬15 s（5000步）與CPU數(shù)據(jù)交換1次后，計(jì)算耗時(shí)趨于穩(wěn)定，約為11～12 min，與CPU單核計(jì)算相比，計(jì)算效率提高約180倍。

圖2 有激波的過丘流數(shù)值模擬

圖3 GPU并行計(jì)算效率統(tǒng)計(jì)

3 廈門海域數(shù)學(xué)模型

3.1 研究范圍與模型參數(shù) 廈門海域數(shù)學(xué)模型范圍如圖4所示，包括同安灣、九龍江口、圍頭灣等廈門島及周邊水域，模型北側(cè)邊界至泉州市崇武鎮(zhèn)，南邊界至漳州市東山縣北側(cè)，模型總長約178 km，寬約105 km。外海邊界采用潮位控制，由全球潮汐預(yù)報(bào)模型Nao-test計(jì)算。模型網(wǎng)格總數(shù)為157 818個(gè)，工程水域網(wǎng)格加密，見圖4，最大網(wǎng)格邊長1825 m，最小網(wǎng)格邊長30 m，模型水深采用2019年最新電子海圖數(shù)據(jù)。數(shù)學(xué)模型時(shí)間步長取1 s，糙率采用附加糙率0.013+0.012/h（h表示水深，當(dāng)h小于1 m時(shí)，糙率取0.015），紊動(dòng)黏性系數(shù)采用0.1hU*計(jì)算（U*為摩阻流速）。

圖4 數(shù)學(xué)模型范圍及網(wǎng)格示意圖

3.2 模型驗(yàn)證 模型采用2019年3月22—23日同安灣大潮同步水文測驗(yàn)資料和2019年國家海洋信息中心廈門站潮汐表資料作為模型驗(yàn)證資料，驗(yàn)證點(diǎn)位置見圖4。圖5和圖6為2019年3月大潮潮位驗(yàn)證，圖7和圖8為對應(yīng)的潮流驗(yàn)證。由圖5—圖8知，潮位計(jì)算值與實(shí)測值吻合較好，高低潮位偏差基本在0.1 m之內(nèi)，流速和流向?qū)崪y值和計(jì)算值吻合亦較好。

圖5 T1站潮位驗(yàn)證

圖6 T2站潮位驗(yàn)證

圖7 A3站潮流驗(yàn)證

圖8 A6站潮流驗(yàn)證

圖9為廈門站2019年全年潮汐表潮位驗(yàn)證，采用GPU并行計(jì)算2019年全年的潮波傳播過程共耗時(shí)約30 h，GPU每模擬30 min與CPU數(shù)據(jù)交換1次，采用CPU單核計(jì)算，耗時(shí)約5475 h，計(jì)算效率提高約180倍。對于長時(shí)序驗(yàn)證，可用Skill值表達(dá)模型驗(yàn)證質(zhì)量，Skill值計(jì)算見式（11），當(dāng)Skill值為1時(shí)，模型驗(yàn)證完美，Skill值越小，計(jì)算偏差越大，本文廈門站潮位統(tǒng)計(jì)Skill值為0.926，說明數(shù)學(xué)模型計(jì)算的廈門附近海域潮波傳播相似性較好。

圖9 廈門站潮汐表潮位驗(yàn)證

式中：Xmod為模型值；Xobs為觀測值（本文為潮汐表預(yù)報(bào)值）；為觀測平均值。

4 潮波傳播空間分布特征

4.1 汊道潮流特征 廈門周邊海域島嶼眾多，地形邊界相當(dāng)復(fù)雜，漲落潮流在不同區(qū)域存在時(shí)間或者空間的不平衡。漲潮時(shí)（圖11（a）），外海水流自圍頭灣和廈門灣兩個(gè)口門向廈門島方向運(yùn)動(dòng)，廈門灣口門處的漲潮流運(yùn)動(dòng)至金門島西側(cè)時(shí)，在廈門島的阻水作用下，水流分為向九龍江口和同安灣兩個(gè)方向的水流，向同安灣方向的水流與圍頭灣的漲潮流在廈門島東北側(cè)匯合，進(jìn)而向同安灣運(yùn)動(dòng)。落潮時(shí)（圖11（b）），同安灣的水體自廈門島東西兩側(cè)向外海運(yùn)動(dòng)，金門島北側(cè)由于水深相對較淺，大部分落潮流向圍頭灣運(yùn)動(dòng)，小部分向烈?guī)Z鄉(xiāng)方向運(yùn)動(dòng)，九龍江口的落潮流在廈門島西南側(cè)與同安灣的落潮流匯合，并向外海運(yùn)動(dòng)。

圖11 廈門附近海域漲落急流場分布圖

表1統(tǒng)計(jì)了廈門島周邊的年平均漲落潮流量，斷面位置見圖10。由表1可知，廈門島南側(cè)（斷面A）與東側(cè)漲潮流量（斷面B+C）基本相當(dāng)，圍頭灣處的漲潮流量（斷面D）相對較小，不足前述斷面的三分之一，落潮流量分布趨勢與漲潮基本一致，由于金門島北側(cè)部分落潮流向烈?guī)Z鄉(xiāng)方向，因此廈門島東側(cè)落潮流量大于南側(cè)。由此可見，廈門島周邊的潮量主要來自廈門灣口門。

表1 廈門灣各汊道年平均流量統(tǒng)計(jì) （單位：m3/s）

圖10 采樣點(diǎn)位置圖

4.2 汊道潮位特征 圖12為廈門附近海域汊道沿程年平均高、低潮位分布特征。潮波自外海向口門內(nèi)傳播，在徑流、地形及岸線邊界的影響下，高潮位總體呈逐漸升高的趨勢，圍頭灣及九龍江口沿線高潮位變化相對較平緩，同安灣沿線高潮位在廈門島東側(cè)波動(dòng)較大，主要因?yàn)樽詮B門灣口門和圍頭灣口門進(jìn)入的潮波在廈門島東北側(cè)（5#點(diǎn)附近）匯合，而金門島北側(cè)（6#點(diǎn)）高潮位明顯高于廈門島東側(cè)（19#點(diǎn)），因此在水流匯合處（5#點(diǎn)）高潮位升高幅度較大，約為0.11 m，且同安灣高潮位明顯高于九龍江口和圍頭灣。低潮位自外海至口門內(nèi)，圍頭灣及九龍江口沿線呈逐漸降低的趨勢，潮差逐漸增大；同安灣沿線低潮位變化較為復(fù)雜，對于灣內(nèi)的1#—5#點(diǎn)影響，受徑流、地形及圍頭灣低潮位低的影響，越往灣頂?shù)统蔽辉礁撸?9#以外各點(diǎn)低潮位受外海低潮位的影響，低潮位升高約0.10 m。由于廈門島東北側(cè)是廈門灣和圍頭灣兩股潮波傳播的交匯點(diǎn)，因此該處高低潮位波動(dòng)較大，但總體來看交匯處5#采樣點(diǎn)的高低潮位與圍頭灣處的6#點(diǎn)更為接近，與廈門島東側(cè)的19#點(diǎn)差異較大，因此同安灣的潮波特性與圍頭灣潮波特性更為相似。

圖12 廈門附近海域沿程年平均高、低潮位分布

5 潮汐調(diào)和分析

潮波從外海進(jìn)入河口區(qū)后，口門位置和形態(tài)的差異導(dǎo)致潮波在傳播過程中的變形程度也不同，本節(jié)采用具有自主知識(shí)產(chǎn)權(quán)的CJK-TIDE軟件，對數(shù)學(xué)模型計(jì)算結(jié)果進(jìn)行調(diào)和分析計(jì)算，以此分析潮波在傳播過程中的特征和差異，分潮數(shù)采用63個(gè)分潮。圖13為廈門島附近18#點(diǎn)（位置見圖10）部分分潮振幅的分布，由圖13知，半日分潮M2分潮占絕對主導(dǎo)地位，其振幅占所有分潮振幅和的41.2%，S2分潮振幅次之，二者占所有分潮振幅和的53.3%，全日分潮K1分潮振幅值最大。圖14和圖15為廈門附近海域半日分潮M2、全日分潮K1分潮振幅與初相角沿程分布，由圖14和圖15知，M2分潮自外海向同安灣及九龍江口方向，振幅和初相角總體呈增大趨勢，同安灣與外海的振幅比約為1.13，九龍江口與外海的振幅比約為1.24，潮波自圍頭灣向同安灣方向傳播時(shí)，至5#、6#點(diǎn)時(shí)，振幅降低，這主要因?yàn)樵撎幨軓B門灣和圍頭灣兩股潮波傳播影響，而廈門灣的振幅明顯低于圍頭灣（圖14）；K1分潮自圍頭灣向同安灣傳播過程中，振幅先增大后減小，在圍頭灣7#位置達(dá)到最大，K1分潮自外海向九龍江口傳播過程中，振幅持續(xù)增大，由圖15可明顯看出，九龍江口全日振幅明顯大于同安灣。圖16為廈門附近海域淺水影響系數(shù)沿程分布變化，淺水影響系數(shù)一般用HM4/HM2表示，表征潮波變形的程度，由圖16知，淺水影響系數(shù)自外海向同安灣和九龍江口方向均呈逐漸增大的趨勢，說明潮波向岸邊傳播過程中，潮波變形逐漸加劇，同安灣潮波變形程度大于九龍江口。

圖13 廈門島附近18#采樣點(diǎn)部分分潮振幅分布

圖14 廈門附近海域M2分潮振幅與初相角沿程分布特征

圖15 廈門附近海域K1分潮振幅與初相角沿程分布特征

圖16 廈門附近海域淺水影響系數(shù)沿程分布特征

6 結(jié)論

（1）構(gòu)建了廈門及其附近海域高分辨率二維潮流數(shù)學(xué)模型，研究了GPU并行計(jì)算效率，結(jié)果表明GPU并行較CPU單核計(jì)算效率可提高約180倍。

（2）廈門島南側(cè)（斷面A）與東側(cè)漲落潮流量（斷面B+C）較為接近，圍頭灣處的漲落潮流量（斷面D）相對較小，說明廈門島附近潮量主要來自廈門灣。潮波自外海向九龍江口、圍頭灣方向傳播，高潮位逐漸升高，低潮位逐漸降低，潮差增大。同安灣方向潮波受廈門灣和圍頭灣兩股潮波影響，其潮波特性與圍頭灣更為相似。

（3）潮汐調(diào)和分析結(jié)果表明，廈門海域M2分潮占絕對主導(dǎo)地位，其振幅占所有分潮的41.2%，S2分潮振幅次之，二者占所有分潮振幅和的53.3%，全日分潮K1分潮振幅值最大。同安灣與九龍江口M2分潮振幅接近，K1分潮九龍江口振幅較大，潮波變形程度同安灣較大。