李春良，倪曉雯，梅國永

（1.山東省交通規(guī)劃設(shè)計院，山東 濟(jì)南 250031；2.山東省建筑科學(xué)研究院，山東 濟(jì)南 250031）

?

FVCOM在長江口水動力數(shù)值模擬中的應(yīng)用

李春良1，倪曉雯2，梅國永2

（1.山東省交通規(guī)劃設(shè)計院，山東 濟(jì)南 250031；2.山東省建筑科學(xué)研究院，山東 濟(jì)南 250031）

本文以“大江河口濕地演變退化的評價體系”項目為背景，應(yīng)用FVCOM潮流及形態(tài)動力學(xué)模型建立長江口三維潮流數(shù)值計算模型，建立了包括長江口、杭州灣及鄰近海域大范圍的三維潮流數(shù)值模型，基于Linux平臺下的并行計算使得變尺度大范圍河口地區(qū)的模擬效率得到了很大的提高。運用實測潮位、流速、流向?qū)δＰ偷南嗨菩赃M(jìn)行驗證，計算驗證結(jié)果與實測值比較吻合，模擬流場能夠比較好的反映長江口地區(qū)往復(fù)流場和口外區(qū)域順時針旋轉(zhuǎn)流特征，可以用于長江口潮流的進(jìn)一步研究。

長江河口區(qū)；潮流界；三維；數(shù)值模擬

1 長江口概況

1.1 概況

長江口是長江在東海入?？诘囊欢嗡颍瑢儆谳^為典型的潮汐河口，潮區(qū)界位于安徽省銅陵市和蕪湖市之間，距離長江河口約640 km，潮流界在江蘇省江陰市以下，長度約240 km。按照河口地區(qū)潮流潮汐特征，通常把上自安徽省銅陵市大通鎮(zhèn)，下至水下三角洲地區(qū)前緣（約東經(jīng)123°）的河段稱為河口區(qū)［1］。在徐六經(jīng)以下，由于科氏力的作用，落潮流偏向南、漲潮流偏向北的現(xiàn)象較為明顯，長江在崇明島西側(cè)被分為南支和北支，在吳淞口以下南支又被長興島和橫沙島分為兩支，即南港和北港，南港又因為九段沙的分割分為南槽和北槽，河槽自西往東呈現(xiàn)較有規(guī)律的分叉。最終南、北支，南、北港，南、北槽呈三級分汊、四口入海的格局［2］（圖1）。

圖1　長江口河勢

1.2 潮汐與潮流

長江河口地區(qū)屬于不正規(guī)淺海半日潮，潮汐現(xiàn)象主要受外海的潮流潮波影響，潮汐日不等現(xiàn)象較為顯著［2］。潮波進(jìn)入長江河口地區(qū)后，受到岸灘、河底河床抬高和上游徑流等因素的影響，潮波在上游的變形要大于下游，上游潮位要高于下游，上游潮差要小于下游，越往上游漲潮歷時逐步縮短，落潮歷時相應(yīng)逐步延長［3］。

長江河口以外海域的潮流多為順時針旋轉(zhuǎn)流，潮位和流速在時間節(jié)點上基本能夠保持一致。潮流在河口進(jìn)入各個河槽以后，因受到河底和岸灘約束多呈往復(fù)流形態(tài)，水位和流速隨時間發(fā)生變化，不再同步且相位發(fā)生偏差。長江口的南支和北支相對而言，南支無論從漲潮量還是落潮量上都大于北支，南支的落潮流占主導(dǎo)地位，但是北支從 1958年以后，徑流量分配比逐步減少，漲潮量大于落潮量成為趨勢，演變?yōu)闈q潮流占主導(dǎo)地位［2］。

2 FVCOM海洋數(shù)值模型

FVCOM（Finite Volume Coast and Ocean Model）是美國麻州大學(xué)陳長勝教授研究課題組編程開發(fā)的，是基于無結(jié)構(gòu)三角形網(wǎng)格架構(gòu)、有限體積、三維預(yù)報原始方程近岸海洋數(shù)值模型。其計算原理是通過在水平方向上建立非結(jié)構(gòu)化三角形網(wǎng)格，運用有限體積法解控制方程，通過積分的運算方式計算三角形控制體單元中的通量［4］。

2.1 FVCOM主要特征

1）垂向混合系數(shù)

FVCOM的垂向混合系數(shù)是基于二階湍流閉合模型來明確的，其目的是為了盡量擺脫人為干擾因素的影響。由湍流假設(shè)理論推導(dǎo)形成閉合湍流模型后移植到FVCOM模型中來［4］。

該模式在垂直方向上采用了s坐標(biāo)系，對于處理地形的顯著變化具有良好的適應(yīng)性，這對于研究受潮汐影響的河口海岸水動力學(xué)非常有必要。

3）差分格式和模態(tài)分離

模型采用顯式水平時間差分和隱式垂向時間差分，從而在水平方向和垂直方向都有比較高的分辨率。其算法采用的是時間分裂算法，其中，外模應(yīng)用二維方程求解，選取的時間步長較短；內(nèi)模應(yīng)用三維方程來求解，所選取的時間步長較長，采用內(nèi)外模式能夠節(jié)省大量的計算時間［4］。

4）干濕判別和并行計算

干濕處理法目前已經(jīng)廣泛應(yīng)用于河口潮間帶的水動力模型數(shù)值模擬中，通過移動邊界在不同時間步長時的干濕來判斷，干是不過水狀態(tài)，濕是過水狀態(tài)。模型運用MPI平行計算模塊，在方便調(diào)試和運行的同時極大的節(jié)約了運算時間。

2.2 原始方程組

模式控制方程組由不可壓連續(xù)、動量、溫度、密度和鹽度方程［4～6］等組成聯(lián)合方程組：

2.3非結(jié)構(gòu)網(wǎng)格設(shè)計

本模型在計算區(qū)域內(nèi)采用的非結(jié)構(gòu)化三角形網(wǎng)格，計算范圍內(nèi)由三角形單元組成，分別有節(jié)點、中心和三邊組成，并且每個單元都不重合。把變量和等設(shè)在三角單元節(jié)點上，u、v設(shè)在三角單元中心，這樣的設(shè)定將高度、速度、溫度和研通量的計算準(zhǔn)確性大大提高了。通過選取所有相鄰三角單元中心所圍成切面的凈通量來計算節(jié)點上的變量值，通過選取進(jìn)出三角單元三邊的凈通量來計算中心點的變量值［7］。

圖2　非結(jié)構(gòu)網(wǎng)格中參量位置示意

2.4 Linux與MPI

計算程序以 Fortran77為藍(lán)本編譯，在河海大學(xué)聯(lián)想機群網(wǎng)絡(luò)系統(tǒng)上運行計算。機群操作系統(tǒng)采用局域網(wǎng)連接和星型拓?fù)浣Y(jié)構(gòu)［8］，集群所有節(jié)點上的操作系統(tǒng)均采用RedHat Linux 9.0。

圖3　MPI的執(zhí)行過程

MPI是可以被 Fortran77/C/Fortran90/C++調(diào)用的函數(shù)庫系統(tǒng)，其具有可移植性、高效性、靈活性和易用性的特點，目前MPICH是使用最廣泛的版本，它既可以在純 Linux環(huán)境下運行，也可以在Windows環(huán)境下運行，并且使用比較方便。本文的計算就是基于 Linux平臺，運用 MPICH進(jìn)行的FVCOM并行計算。

3 數(shù)值模型的計算與驗證

3.1 模擬范圍及網(wǎng)格

長江口徑流和潮汐等水動力作用比較強，模型的模擬計算范圍應(yīng)注意潮流界在河口地區(qū)的變動范圍和潮汐作用對河口地區(qū)的綜合影響，減小選取的邊界條件對研究區(qū)域的干擾。模型選定范圍自西往東從江陰開始，東到外海50 m等深線，南北邊界介于呂四港和舟山島南側(cè)之間，包括長江口河口區(qū)域、杭州灣及相鄰部分外海區(qū)域［9］，鑒于模型計算范圍較大，為提高運算速度和效率，對部分復(fù)雜岸邊界進(jìn)行了簡化，從大范圍的角度分析長江口的三維潮流動力作用。模型的計算范圍示意以及模型范圍的地形見圖4、圖5。

圖4　模型計算范圍示意

圖5　模型計算范圍地形高程填充

圖6　計算網(wǎng)格

圖7　深水航道網(wǎng)格加密示意

模型采用 SMS軟件生成全三角形網(wǎng)格，共有52 085個單元，27 099個節(jié)點，將水深自動插值到網(wǎng)格節(jié)點上，網(wǎng)格尺度250～10 000 m。計算網(wǎng)格見圖 6。模型在上邊界、崇明島、長興島、深水航道和舟山群島附近都進(jìn)行了加密，加密區(qū)域的網(wǎng)格尺度為250～500 m，深水航道段加密網(wǎng)格見圖7。由于FVCOM模式采用的是干濕判別法，所以為了節(jié)省計算時間，我們將計算范圍中的邊界和島嶼進(jìn)行了整合，利于模型計算更加穩(wěn)定。

3.2 地形資料及模型參數(shù)

模型的地形采用 2005年實測水下地形圖，計算基準(zhǔn)面采用 1985年國家高程基面。采用冷啟動方式進(jìn)行初始啟動，設(shè)定初始為0的水位和流速，上游邊界采用實地測量的流量數(shù)值進(jìn)行控制輸入，外海邊界采用東中國海潮波數(shù)值模型提供的潮位控制輸入，在時間序列上連續(xù)模擬潮流場。結(jié)合網(wǎng)格尺度，我們估算時間步長為2.0 s。

FVCOM為三維大范圍海洋數(shù)值模式，底部摩阻由 brough_ud.F子代碼控制，在計算過程中，可以根據(jù)自己需要修改代碼分塊設(shè)定糙率。

3.3 模型驗證及結(jié)果分析

模型采用實測水文站統(tǒng)計的數(shù)據(jù)進(jìn)行對比驗證，通過實地測量的流速流向值、潮位值和垂直方向的平均流速值與模型預(yù)測數(shù)據(jù)進(jìn)行對比分析。

1）潮位驗證分析

對比各潮位觀測站驗證曲線，除位于北支中段的青龍港站外，絕大多數(shù)觀測點模型推算潮位值與實地觀測值比較接近，誤差低于10 cm，能有效控制在誤差范圍內(nèi)。共選擇 11個潮位對比站點，整個模擬區(qū)域站點布置見圖8。

圖8　潮位驗證站

經(jīng)過對驗證曲線進(jìn)行分析可以得出：堡鎮(zhèn)、橫沙、連興港和高橋站模型推算出的潮位最高值比實地測量的略小，潮位最低值比實地測量的略大，推算的潮差比實地測量數(shù)值略小。經(jīng)分析，認(rèn)為主要因為南北支地形和自然條件比較復(fù)雜，個別河段河底高程起伏較大，而且FVCOM模式采用的插值地形與實測地形也必然存在一定差異，所以造成個別站點的部分?jǐn)?shù)據(jù)出現(xiàn)誤差是正常的。此外，還有一個比較重要的原因是在模型的底摩阻系數(shù)選擇上，我們?yōu)榱丝紤]大范圍的潮波傳播過程，在大河段采用統(tǒng)一底摩阻系數(shù)，這無疑也影響了模擬的精度。由于外海邊界是由東中國海潮波數(shù)值模型提供，其提供邊界與實際的真實數(shù)據(jù)必然存在的一些誤差，以上幾個原因?qū)е聜€別站點小潮計算值偏高（徐六經(jīng)站、南門站和石洞口站）。除此之外，其余觀測點的驗證結(jié)果比較理想。我們可以得出，F(xiàn)VCOM數(shù)值模型推算的結(jié)果比較可靠，與選取模擬范圍的實際潮波傳播過程相吻合［10］。

2）流速驗證分析

我們選取9個流速驗證觀測點（位置見圖9），對觀測點的表、底層以及垂直方向平均流速模擬值進(jìn)行提取，與實地觀測值進(jìn)行對比分析，并將流向圖與實際進(jìn)行對照。

圖9　流速驗證點示意

對比分析流速驗證圖，F(xiàn)VCOM數(shù)值模型模擬結(jié)果中提取的站點的表層、底層及垂直方向平均流速流向與實地測量結(jié)果驗證較好，對流速和流向的變化能夠比較準(zhǔn)確的反映。河口地區(qū)因受漲潮落潮和徑流相互作用表現(xiàn)為往復(fù)流形態(tài)，模擬得出的落、漲潮流速，落、漲潮歷時均與實際情況相符。北支相對于南支來說，因為北支地形地貌條件比較復(fù)雜，水量較小，模擬的精度也必然受到影響，個別流速流向與實際出現(xiàn)偏差［10］。

3）流場驗證分析

最后根據(jù)模型視頻輸出模式形成漲落潮流場（見圖10），通過觀察可以得出，長江河口地區(qū)南、北支，南、北港，南、北槽漲落潮的水流流向和流態(tài)變化均能清楚的呈現(xiàn)出來。模擬范圍河段內(nèi)潮流的往復(fù)流形態(tài)和口門外順時針旋轉(zhuǎn)流形態(tài)也能夠較好的呈現(xiàn)。

圖10　長江口落、漲急流場

4 結(jié) 語

本文對FVCOM海洋數(shù)值模式基本理論以及模型在Linux操作系統(tǒng)環(huán)境下的構(gòu)建和應(yīng)用進(jìn)行了簡要介紹，并對利用MPI并行在校園網(wǎng)絡(luò)平臺上進(jìn)行計算的原理進(jìn)行了較為詳細(xì)的闡述。

利用FVCOM數(shù)值模型對長江河口大范圍的三維潮流進(jìn)行模擬，并根據(jù)提取的節(jié)點表、底以及垂向平均流速與實地測量的數(shù)據(jù)進(jìn)行對比分析，模擬的結(jié)果吻合度較高，往復(fù)流形態(tài)和長江口門地區(qū)順時針旋轉(zhuǎn)流的特征明顯，說明該模型在潮流水動力特性方面驗證合理，可以用來模擬長江口潮流水動力狀況。

［1]禔李來， 李誼純. 長江口整治工程對鹽水入侵影響研究［J］. 海洋工程， 2005.

［2]丁遠(yuǎn)靜. 長江口擬建工程對枯季水動力環(huán)境的影響研究［D］. 河海大學(xué)， 2007.

［3]鄭金海， 諸裕良. 長江河口鹽淡水混合的數(shù)值模擬計算［J］. 海洋通報， 2001， 20（4）: 1-10.

［4]朱學(xué)明. 中國近海潮汐潮流的數(shù)值模擬與研究［D］. 中國海洋大學(xué)， 2009.

［5]Haosheng Huang， Changsheng Chen， Jackson O Blanton，et al. A numerical study of tidal asymmetry in Okatee Creek， South Carolina［J］. Estuarine， Coastal and Shelf Science， 2008.

［6]Changsheng Chen， Robert C Beardsley， Geoffrey Cowles. An unstructured grid， finite-volume coastal ocean model FVCOM user manual［J］.

［7]汪守東. 基于 Lagrange追蹤的海上溢油預(yù)報模型研究［D］.大連理工大學(xué)， 2008.

［8]陳林. 基于 Linux機群的大型結(jié)構(gòu)并行有限元方法研究［D］. 河海大學(xué)， 2006.

［9]謝銳. EFDC模型在長江口及相鄰海域三維水流模擬中的開發(fā)應(yīng)用［J］. 水動力學(xué)研究與進(jìn)展: A輯， 2010.

［10]王振華. 海平面上升對長江口水動力及污染物運動軌跡影響的數(shù)值模擬研究［D］. 大連理工大學(xué)， 2011.

Application of FVCOM in Numerical Simulation of Hydrodynamic Force at Yangtze River Estuary

Li Chunliang1，Ni Xiaowen2，Mei Guoyong2
（1.Shandong Provincial Communications Planning and Design Institute， Jinan Shandong 250031， China；2.Shandong Academy of Building Research， Jinan Shandong 250031， China）

Based on the evaluation system for the evolution and degradation of great river estuarine wetlands， 3D numerical model is built to calculate the tidal current at Yangtze River Estuary by using FVCOM tidal flow and form dynamics model. By now， 3D tidal current models have been established for Yangtze River estuary， Hangzhou Bay and the surrounding sea area within a large range. The parallel computation on Linux platform can improve the variable-metric simulation of river estuarine area within a large range. The above models are tested from the aspects of the measured tidal level， velocity and direction of current. The calculation results agree very well with the measured values. The simulated flow field shows the reversing flow filed at Yangtze River delta and the clockwise rotation characteristics outside the river estuary well， which applies to the further study of the tidal current at Yangtze River estuary.

Yangtze River estuary area； tidal current limit； three-dimensional； numerical simulation

TV148+.5

A

1004-9592（2016）03-0001-05

10.16403/j.cnki.ggjs20160301

2016-01-14

李春良（1983-），男，工程師，主要從事港口及航道工程設(shè)計和海岸動力數(shù)值模擬工作。