李春良，倪曉雯，梅國(guó)永

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

?

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

李春良1，倪曉雯2，梅國(guó)永2

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

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

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

1 長(zhǎng)江口概況

1.1 概況

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

圖1　長(zhǎng)江口河勢(shì)

1.2 潮汐與潮流

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

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

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

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

2.1 FVCOM主要特征

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

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

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

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

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

4）干濕判別和并行計(jì)算

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

2.2 原始方程組

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

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

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

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

2.4 Linux與MPI

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

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

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

3 數(shù)值模型的計(jì)算與驗(yàn)證

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

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

圖4　模型計(jì)算范圍示意

圖5　模型計(jì)算范圍地形高程填充

圖6　計(jì)算網(wǎng)格

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

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

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

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

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

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

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

1）潮位驗(yàn)證分析

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

圖8　潮位驗(yàn)證站

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

2）流速驗(yàn)證分析

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

圖9　流速驗(yàn)證點(diǎn)示意

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

3）流場(chǎng)驗(yàn)證分析

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

圖10　長(zhǎng)江口落、漲急流場(chǎng)

4 結(jié) 語(yǔ)

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

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

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

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

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

［4]朱學(xué)明. 中國(guó)近海潮汐潮流的數(shù)值模擬與研究［D］. 中國(guó)海洋大學(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ù)報(bào)模型研究［D］.大連理工大學(xué)， 2008.

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

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

［10]王振華. 海平面上升對(duì)長(zhǎng)江口水動(dòng)力及污染物運(yùn)動(dòng)軌跡影響的數(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è)計(jì)和海岸動(dòng)力數(shù)值模擬工作。