無結(jié)構(gòu)網(wǎng)格二維海洋模式的正壓梯度力算法改進(jìn)

陳睿

（國家海洋局東海預(yù)報中心,上海201206）

?

無結(jié)構(gòu)網(wǎng)格二維海洋模式的正壓梯度力算法改進(jìn)

（國家海洋局東海預(yù)報中心,上海201206）

摘要：基于一套自主研制的無結(jié)構(gòu)網(wǎng)格二維河口海洋數(shù)值模式A2D，在大圓湖理想模型下，通過與解析解進(jìn)行比較分析，采用不同架構(gòu)配置，改進(jìn)設(shè)計正壓梯度力計算方法。改進(jìn)后的算法中引入了從算架構(gòu)的配置，以配合主算架構(gòu)，得到更佳的穩(wěn)定性。通過水位場平面分布與單點過程線可以發(fā)現(xiàn)，三組試驗的算法均獲得了較好的精度和比原算法更好的穩(wěn)定性，其中TSNS配置算法（中心點計算水位、邊中點計算流速的主算架構(gòu)，配合節(jié)點計算水位、邊中點計算流速的從算架構(gòu)）由于其主算架構(gòu)更接近結(jié)構(gòu)網(wǎng)格下的C網(wǎng)格，在守恒性、移動潮灘邊界處理等方面具有一定優(yōu)勢和便利性，有利于在實際海洋中的計算。將TSNS配置算法在江浙沿海進(jìn)行試算，水位驗證結(jié)果與實測基本符合，與原A2D模式計算水位之間無顯著差異。TSNS算法在穩(wěn)定性方面的改進(jìn)，有助于提升模式升級為三維后的穩(wěn)定性，為今后模式成功升級為三維打下基礎(chǔ)。

關(guān)鍵詞：無結(jié)構(gòu)網(wǎng)格；海洋模式；解析解；正壓梯度力；穩(wěn)定性

1　引言

無結(jié)構(gòu)網(wǎng)格憑借其易擬合岸界、可局部加密等優(yōu)勢，正越來越多地被用于海洋水動力數(shù)值模式。當(dāng)前國際上知名的模式多來自國外[1-12]，其中采用無結(jié)構(gòu)網(wǎng)格的模式有FVCOM[10-11]、ADCIRC[12]等，這些模式在國內(nèi)被廣泛使用，而國內(nèi)自主開發(fā)的模式[13-14]則非常缺乏。

鑒于此，2006—2011年期間，筆者所在研究組自主研制了一套無結(jié)構(gòu)網(wǎng)格二維河口海洋數(shù)值模式[15-16]，其中2011年最終版本后文稱為A2D[16]。A2D模式采用無結(jié)構(gòu)三角形網(wǎng)格，基于有限體積法求解。水位在網(wǎng)格中心（重心）求解，流速在網(wǎng)格邊中點求解（同時求解x和y方向流速）。這種配置結(jié)合了各種已有的無結(jié)構(gòu)網(wǎng)格模式的長處。一方面，它與國際上較為通行的C網(wǎng)格[17-18]配置相近，具有在網(wǎng)格中心求解水位之利于干濕判別的優(yōu)點；另一方面，直接求解x-方向和y-方向流速使得離散簡潔高效，無須水平坐標(biāo)轉(zhuǎn)換，且能更方便地利用有限體積法，提升守恒性。A2D模式在時間上顯式求解，采用預(yù)估修正法[19-22]。通過理想試驗、黃浦江和長江口等環(huán)境下的試驗計算，有效地驗證了A2D模式的精度[16]。

但作為新生事物，A2D模式尚不成熟，存在一些不足，如模式還只是二維、在理想試驗下有微小數(shù)值波動等。筆者所在研究組曾將A2D模式以相同架構(gòu)和算法升級至三維版本（稱之為A3D）[23]，但在對長江口海域的計算中發(fā)現(xiàn)A3D的水動力不如A2D穩(wěn)定，所以三維版本暫未獲得成功。A3D的不穩(wěn)定性，很有可能源自于A2D的架構(gòu)和算法，其中求解流速時正壓梯度力項為主要項，是關(guān)注的重點。本文對A2D的正壓梯度力算法進(jìn)行分析與改進(jìn)，提升模式在一個理想的大圓湖解析解模型下的穩(wěn)定性表現(xiàn)，并投入實際海洋中試算。

2　大圓湖理想模型簡介

Csanady提出了一個理想的大圓湖模型[24]，計算了該圓湖在風(fēng)的作用下的表明波動。Birchfield指出了其中一些錯誤并給出了正確的解析解表達(dá)式[25]。該模型為一個半徑67.5 km，水深75 m的平底大圓湖（見圖1），初始時刻水位靜止為0，施加以恒定西風(fēng)3 m/s，科氏力系數(shù)f取常數(shù)0.0001，忽略底摩擦。在風(fēng)與科氏力的共同作用下，湖表將會產(chǎn)生水位波動，并且可以得到該波動的解析解表達(dá)式。通過這個模型，可以將數(shù)值模式的數(shù)值解與解析解進(jìn)行比較，從而探討模式的精度。

圖1　大圓湖網(wǎng)格及輸出點A位置

由于這個大圓湖的空間尺度很大，根據(jù)量綱分析的結(jié)果，流速平流項和擴(kuò)散項對水位和流速變化的貢獻(xiàn)極其微小，可忽略不計[24-25]，所以影響模式計算精度的主要項為正壓梯度力項，因此A3D在此模型下的穩(wěn)定性問題，應(yīng)是緣于正壓梯度力項。重新設(shè)計模式的算法，使得正壓梯度力項得到更合理地求解，并在大圓湖理想模型下取得較好的結(jié)果，是本文的主要目標(biāo)。

圖2　A2D模式的數(shù)值解（紅色）與解析解（藍(lán)色）的過程線比較（U是東向流速，D是水深）

圖3　A3D模式的數(shù)值解（黑色）與解析解（紅色）的水位過程線比較

圖4　解析解在1 h、1 d、5 d時的水位場分布

圖5　A2D模式在1 h、1 d、5 d時的水位場分布

3　原模式的穩(wěn)定性分析

采用A2D和A3D分別對大圓湖進(jìn)行模擬（網(wǎng)格見圖1），時間步長均為5 s，并在A點輸出站位過程線作比較。模式A2D與解析解的過程線比較結(jié)果如圖2所示?？梢钥吹?，模式的計算結(jié)果與解析解十分吻合，且總體比較光滑，沒有出現(xiàn)明顯的不穩(wěn)定，說明模式的水動力達(dá)到了較高的精度。但在4.89 d左右，水位的數(shù)值解（紅色）存在微小的波動，即存在微小的不穩(wěn)定跡象。而在相同架構(gòu)和算法的三維模式A3D下，水位的不穩(wěn)定性表現(xiàn)得更為明顯（圖3）。對比解析解（見圖4）和A2D（見圖5）的水位場分布也可以發(fā)現(xiàn)，在模式5d時，A2D的計算結(jié)果的等值線略顯抖動，這也反映出模式在穩(wěn)定性上略有欠缺。

A2D在理想試驗下的微小不穩(wěn)定，在實際海洋中計算時并不會發(fā)生，因為無頻散的流速平流項會使得模式保持穩(wěn)定。但當(dāng)三維開發(fā)版本A3D模式在實際海洋中計算時，由于不穩(wěn)定性增大，所以難免出現(xiàn)個別區(qū)域流速或流向異常，影響計算結(jié)果的正確性。所以，本文改進(jìn)A2D正壓梯度力項算法，使得模式更為穩(wěn)定，是模式升級為穩(wěn)定的三維版本的基礎(chǔ)。

4　算法改進(jìn)和探討

4.1原模式控制方程組

A2D模式包含水動力計算和鹽度計算，而溫度暫時不作計算，取為常數(shù)T=10°C。本文不牽涉溫鹽算法的改進(jìn)，主要對水動力控制方程組的相關(guān)部分作簡介。對流體不可壓縮、Boussinesq和靜力近似下的海洋動力學(xué)原始方程組作垂向積分，可得到垂向平均的二維控制方程組：

式中：t為時間，ζ表示海表水位，D代表總水深（總水深D=H+ζ，H為固定不變的基準(zhǔn)水深），分別表示水平x方向和y方向的垂向平均流速（其中u和v分別為空間中某點的水平x方向和y方向的局地流速），F(xiàn)x和Fy為x方向和y方向的流速水平擴(kuò)散，AM為水平湍流粘滯系數(shù)，f為柯氏力系數(shù)，g為重力加速度，和為海表應(yīng)力，和為海底摩擦應(yīng)力，ρ為密度，ρ0為參考密度。其中密度ρ=ρ（S,T,p）為鹽度S、溫度T（°C）和海壓p（Pa）的函數(shù)，采用1980年國際海水狀態(tài)方程（EOS80）計算。

海表應(yīng)力和為風(fēng)應(yīng)力在x和y方向上的分量：

式中：ρa(bǔ)為空氣密度；Ua和Va分別為x和y方向的風(fēng)速，為風(fēng)矢量的絕對值大??；CD為海水對風(fēng)的拖曳系數(shù)，它根據(jù)Large和Pond改進(jìn)的穩(wěn)定狀態(tài)拖曳系數(shù)計算[26]：

海底摩擦應(yīng)力和根據(jù)謝才公式得到：

式中：nb為曼寧系數(shù)，一般取值為0.015到0.018之間。

4.2原模式正壓梯度力算法分析

在已有的A2D模式中，正壓梯度力項的計算方法如圖6所示。構(gòu)造至多4個控制體A1、A2、A3和A4（未必一定有4個，網(wǎng)格邊緣和干點附近會有缺失），每個控制體的3個頂點均有計算得到的水位，于是可以通過格林公式計算控制體內(nèi)的水位梯度。再將各控制體的水位梯度平均，得到邊j上的水位梯度，從而計算得到邊j上的正壓梯度力結(jié)果。

觀察此算法，較為顯著的問題為未能利用到節(jié)點水位進(jìn)行正壓梯度力計算，控制體遠(yuǎn)端的網(wǎng)格點距離較遠(yuǎn)，不利于數(shù)值穩(wěn)定。而如果對水位進(jìn)行插值，則由于水位梯度計算本來就對空間位置敏感，會導(dǎo)致精度不高，這在A2D研發(fā)階段已做過嘗試，其效果不如A2D最終方案理想[16]。

圖6　A2D正壓梯度力算法示意圖，淺藍(lán)色區(qū)域為控制體A1，至多有4個控制體

在當(dāng)前配置架構(gòu)下暫時無法找到更好算法的情況下，嘗試更多不同的配置架構(gòu)是較好的途徑。簡便起見，將不同配置架構(gòu)的試驗采用其配置特征命名，分別用N、S、T 3個字母代表三角形網(wǎng)格的節(jié)點、邊中點、網(wǎng)格中心，試驗的前兩位字母代表水位計算點和流速計算點的位置。如原A2D模式的命名為：TS。有些試驗采用一種配置架構(gòu)作為主算架構(gòu)，另一種配置架構(gòu)作為從算架構(gòu)，從算架構(gòu)作為第2套同時計算的輔助配置架構(gòu)為主算架構(gòu)提升穩(wěn)定性，這些試驗的命名中第1—2位字母代表了主算架構(gòu)的配置，第3—4位字母代表了從算架構(gòu)的配置，如TSNS即代表網(wǎng)格中心水位、邊中點流速的主算架構(gòu)結(jié)合節(jié)點水位、邊中點流速的從算架構(gòu)的試驗。在經(jīng)過對多種不同配置算法的嘗試后發(fā)現(xiàn)，部分配置架構(gòu)下得到了較穩(wěn)定的結(jié)果，其中包括TSNS。TSNS在節(jié)點計算輔助的水位，有助于提升模式的穩(wěn)定性。這種配置的主算架構(gòu)仍然是中心水位、邊中點流速，但增加計算節(jié)點水位、邊中點流速的從算架構(gòu)，使得中心與節(jié)點同時具有水位，更有利于邊中點位置的正壓梯度力求解。

圖7　求解水位的控制體模型

4.3 TSNS配置算法和計算結(jié)果

在TSNS配置的算法中，通過連續(xù)方程，利用邊中點的流速同時求出網(wǎng)格中心點和網(wǎng)格節(jié)點的水位，并通過動量方程，利用網(wǎng)格中心點和網(wǎng)格節(jié)點的水位求出邊中點的流速。

求網(wǎng)格中心點i的水位時，根據(jù)連續(xù)方程（1），仍然沿用原A2D的求法，在三角形控制體A（圖7a）中根據(jù)格林公式求得：

式中：l為繞A一周的正向曲線，即l的方向為逆時針，A始終在其左邊。

求網(wǎng)格節(jié)點m的水位時，也以格林公式（13）求解，不過控制體A變?yōu)榘鼑?jié)點m的多邊形（圖7b）?？刂企wA有可能完全包圍節(jié)點m，也有可能如圖7b一般缺失部分角度，無論哪種情形，l均為繞A一周的正向曲線。在圖7b的情形下，僅需計算j4-i3-j3-i2-j2-i1-j1上的線積分，因為j1-m-j4上線積分等于0無通量進(jìn)出。每一小段的流速為該小段所連接邊中點的流速。

求邊中點j的流速時，對正壓梯度力項的算法作了修改。在大圓湖模型中，無斜壓梯度力項，平流項也非常小可忽略不計，故對精度影響最大的項正是正壓梯度力項。TSNS配置算法中，動量方程僅正壓梯度力項較原A2D模式有改變，其它項均維持不變。TSNS配置算法在計算邊中點j的流速時，通過格林公式

圖8　求解正壓梯度力項的控制體

由于在以上計算方法中，同時計算了網(wǎng)格中心與網(wǎng)格節(jié)點兩套水位，在時間積分久后兩套的數(shù)值解可能存在分離的現(xiàn)象，所以將從算架構(gòu)中的節(jié)點水位按照一定速度向主算架構(gòu)中的中心點水位回歸。這里取每時間步0.05的回歸比例進(jìn)行趨近。

圖9　TSNS配置結(jié)果分析（U是東向流速，D是水深，下同）

將上述TSNS算法在大圓湖模型中試驗，時間步長仍取5 s，發(fā)現(xiàn)TSNS的精度與原A2D接近，穩(wěn)定性更好（見圖9），第4.89 d未出現(xiàn)不穩(wěn)定抖動，水位場等值線也更平整。

4.4其它配置的計算結(jié)果

除了TSNS配置外，本文還測試了一些其它配置，其中NSTS配置和ST配置這兩組試驗也得到了不錯的結(jié)果。在NSTS配置中，主算與從算架構(gòu)與TSNS對換，NS主算、TS從算，回歸系數(shù)仍取0.05。其結(jié)果與TSNS略有不同（見圖10），精度接近TSNS，穩(wěn)定性較好。而ST配置下也得到了穩(wěn)定性尚可的結(jié)果（見圖11）。通過統(tǒng)計點A水位時間序列的平均誤差和均方根誤差來比較這3組試驗與原A2D模式的精度，可以發(fā)現(xiàn)TSNS配置的精度在所有4種配置中最好（見表1），同時誤差在模式4-5d時較之前時間段有所增大。由于模式最終應(yīng)用時，以TS作為主算架構(gòu)的TSNS配置更接近結(jié)構(gòu)網(wǎng)格下的C網(wǎng)格，在守恒性、移動潮灘邊界處理等方面具有一定優(yōu)勢和便利性，而NSTS配置和ST配置在各種邊界條件設(shè)計中存在一定的困難，故NSTS配置和ST配置不作重點介紹，具體算法細(xì)節(jié)在此省略，僅展示理想試驗下的結(jié)果，后續(xù)在實際海洋中的試算僅基于TSNS配置下進(jìn)行。

圖10　NSTS配置結(jié)果分析

圖11　ST配置結(jié)果分析

5　江浙沿海海域試算

5.1模式網(wǎng)格和基本設(shè)置

基于改進(jìn)了正壓梯度力算法的TSNS配置模式和原A2D模式，對東海區(qū)范圍內(nèi)包含呂四測站、嵊山測站和定海測站的江浙沿海海域進(jìn)行試算，以觀察模擬效果，并進(jìn)行水位驗證，其中TSNS配置模式因配置變化對邊界條件進(jìn)行了部分完善。模式的網(wǎng)格范圍見圖12，基于54坐標(biāo)，包括了整個長江口、杭州灣和鄰近海區(qū)。東邊至124.5°E附近，北邊到34.3°N左右，南邊到28.4°N左右，長江上游邊界取在大通。長江口內(nèi)、深水航道附近和島嶼附近的網(wǎng)格作了局部加密，最小網(wǎng)格分辨率可達(dá)100 m，而口外網(wǎng)格則被放大，最大超過10 000 m。模式中深水航道的導(dǎo)堤和丁壩漲潮時淹沒、落潮時露出，是作為動邊界處理的。

模式的時間步長統(tǒng)一取為1 s。外海開邊界處的水位利用16個分潮（M2，S2，N2，K2，K1，O1，P1，Q1，MU2，NU2，T2，L2，2N2，J1，M1，OO1）的調(diào)和常數(shù)計算得出，而在長江口上游大通處則利用實測徑流量資料給出通量邊界條件。海表面的風(fēng)場以每6 h為1組、分辨率為0.5°×0.5°經(jīng)緯度的氣象預(yù)報后處理結(jié)果給出（可從網(wǎng)址http ://dss.ucar.edu/datasets/ ds744.4/data/處下載）。底摩擦曼寧系數(shù)在模式中設(shè)為0.015。模式從2008年11月5日起計算，共計算30 d。通過搜集的呂四站、嵊山站和定海站的實測水位資料，對原A2D模式和改進(jìn)后的TSNS配置模式進(jìn)行比對驗證。

表1　大圓湖模型下各模式在點A的水位誤差統(tǒng)計表（單位：mm）

5.2水位驗證結(jié)果

模式共運行30 d，輸出第10—30 d的水位，與呂四站、嵊山站和定海站的實測資料進(jìn)行比對（圖13—15，圖中上子圖的藍(lán)實線為原A2D計算水位，黑虛線為TSNS配置計算水位，紅點為實測水位；下子圖的黑實線為TSNS計算水位減去A2D計算水位的差，藍(lán)實線為0軸）?？梢园l(fā)現(xiàn)，TSNS配置模式的水位計算結(jié)果與原A2D模式（TS配置）的水位計算結(jié)果非常接近，兩者差異在3站均不到0.1 m，同時模式計算結(jié)果與實測基本符合，基本反映了3個測站的水位變化規(guī)律，其中嵊山站和定海站存在較明顯的潮汐日不等現(xiàn)象，而定海站的振幅略偏大。

圖12　江浙沿海海域網(wǎng)格

圖13　呂四站TSNS配置與原A2D結(jié)果對比

由于此番改進(jìn)主要目的是探討是模式本身的性能，故在江浙沿海海域進(jìn)行試算并驗證時，底摩擦曼寧系數(shù)統(tǒng)一取了0.015，未針對不同海域進(jìn)行分區(qū)設(shè)置，對驗證效果略有影響。而TSNS配置模式與原A2D模式雖然在算法上存在差異，但由于案例設(shè)定的物理環(huán)境相同，兩者受相同的邊界條件驅(qū)動，故兩者的水位計算結(jié)果與變化特征并沒有出現(xiàn)顯著差異。

圖14　嵊山站TSNS配置與原A2D結(jié)果對比

圖15　定海站TSNS配置與原A2D結(jié)果對比

6　結(jié)論

本文通過對模式配置架構(gòu)的改變，重新設(shè)計正壓梯度力項的算法，得到了3種比原二維A2D模式穩(wěn)定性更高的算法，在大圓湖理想模型下模擬得到了較穩(wěn)定的結(jié)果，同時精度與原A2D接近，均與解析解較為符合。其中TSNS配置算法中由于其主算架構(gòu)TS配置更接近結(jié)構(gòu)網(wǎng)格下的C網(wǎng)格，在守恒性、移動潮灘邊界處理等方面較好，故在完善了相應(yīng)的邊界條件后，在江浙沿海海域進(jìn)行了試算并與實測資料進(jìn)行了對比驗證。水位驗證結(jié)果與實測基本符合，同時與原A2D模式計算水位之間無顯著差異。

原A2D的算法雖然與TSNS算法一樣也能在二維情形的實際海洋中穩(wěn)定，且計算結(jié)果相近，但A2D的三維開發(fā)版本A3D在實際海洋計算中卻遭遇穩(wěn)定性不佳的問題。大圓湖理想模型下更穩(wěn)定的表現(xiàn)，有助于提升模式升級為三維后的穩(wěn)定性，所以TSNS算法在穩(wěn)定性方面的改進(jìn)，為今后模式成功升級為三維打下基礎(chǔ)。經(jīng)過改進(jìn)后的TSNS配置模式可在風(fēng)暴潮模擬等場合進(jìn)行應(yīng)用，同時由于其網(wǎng)格配置特性，可進(jìn)一步優(yōu)化邊界條件設(shè)定，根據(jù)需要靈活設(shè)置堤壩等特殊情形，具有一定的發(fā)展前景。

參考文獻(xiàn)：

[1] Blumberg A F. A Primer for ECOM-si[R]. Technical Report of HydroQual, Mahwah, New Jersey, 1994.

[2] Blumberg A F, Mellor G L. A Description of a Three-Dimensional Coastal Ocean Circulation Model[M]. Washington, DC: American Geophysical Union, 1987.

[3] Casulli V. Semi-implicit Finite Difference Methods for the Two-dimensional Shallow Water Equations[J]. Journal of Computational Physics, 1990, 86(1): 56-74.

[4] Casulli V. A Semi-implicit Finite Difference Method for Non-hydrostatic, Free-surface Flows[J]. International Journal for Numerical Methods in Fluids, 1999, 30(4): 425-440.

[5] Casulli V, Walters R A. An Unstructured Grid, Three-dimensional Model based on the Shallow Water Equations[J]. International Journal for Numerical Methods in Fluids, 2000, 32(3): 331-348.

[6] Ham D A, Kramer S C, Stelling G S, et al. The Symmetry and Stability of Unstructured Mesh C-grid Shallow Water Models Under the Influence of Coriolis[J]. Ocean Modelling, 2007, 16 (1-2): 47-60.

[7] Zhang Y L, Baptista A M, Myers III E P. A Cross-scale Model for 3D Baroclinic Circulation in Estuary-plume-shelf Systems: I. Formulation and Skill Assessment[J]. Continental Shelf Research, 2004, 24(18): 2187-2214.

[8] Fringer O B, Gerritsen M, Street R L. An Unstructured-grid, Finite-volume, Nonhydrostatic, Parallel Coastal Ocean Simulator [J]. Ocean Modelling, 2006, 14(3-4): 139-173.

[9] DHI. MIKE 21 & MIKE 3 Flow Model FM Hydrodynamic and Transport Module Scientific Documentation[R]. Copenhagen: DHI, 2011.

[10] Chen C S, Liu H D, Beardsley R C. An Unstructured Grid, Finite-volume, Three-dimensional, Primitive Equations Ocean Model: Application to Coastal Ocean and Estuaries[J]. Journal of Atmospheric and Oceanic Technology, 2003, 20(1): 159-186.

[11] Chen C S, Beardsley R C, Cowles G. An Unstructured Grid, Finite-volume Coastal Ocean Model: FVCOM User Manual[M]. New Bedford, Mass: SMAST/UMASSD, 2006.

[12] Luettich R, Westerink J. Formulation and Numerical Implementation of the 2D/3D ADCIRC Finite Element Model Version 44.XX[R]. 2004.

[13]顧峰峰,倪漢根,戚定滿,等. AUSM格式在二維淺水方程求解中的應(yīng)用[J].水科學(xué)進(jìn)展, 2008, 19(5): 624-629.

[14]趙棣華,戚晨,庾維德,等.平面二維水流-水質(zhì)有限體積法及黎曼近似解模型[J].水科學(xué)進(jìn)展, 2000, 11(4): 368-374.

[15]陳昞睿,朱建榮,吳輝,等.無結(jié)構(gòu)網(wǎng)格二維河口海岸水動力數(shù)值模式的建立及其應(yīng)用[J].海洋學(xué)報, 2010, 32(2): 31-39.

[16]陳昞睿.無結(jié)構(gòu)網(wǎng)格二維河口海洋數(shù)值模式的研制[D].上海:華東師范大學(xué), 2011.

[17] Arakawa A, Lamb V R. Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model[J]. Methods of Computational Physics, 1977, 17: 173-265.

[18] Kleptsova O, Pietrzak J D, Stelling G S. On the Accurate and Stable Reconstruction of Tangential Velocities in C-grid Ocean Models[J]. Ocean Modelling, 2009, 28(1-3): 118-126.

[19]馮士筰,李鳳岐,李少菁.海洋科學(xué)導(dǎo)論[M].北京:高等教育出版社, 1999.

[20]朱建榮.海洋數(shù)值計算方法和數(shù)值模式[M].北京:海洋出版社, 2003.

[21]朱建榮,楊隴慧,朱首賢.預(yù)估修正法對河口海岸海洋模式穩(wěn)定性的提高[J].海洋與湖沼, 2002, 33(1): 15-22.

[22]朱建榮,朱首賢. ECOM模式的改進(jìn)及在長江河口、杭州灣及鄰近海區(qū)的應(yīng)用[J].海洋與湖沼, 2003, 34(4): 364-374.

[23] Chen B R, Zhu J R, Li L. Accelerating 3D Ocean Model Development by using GPU Computing[M]//Deng W. Future Control and Automation: Proceedings of the 2nd International Conference on Future Control and Automation. Berlin Heidelberg: Springer, 2012, 1: 37-43.

[24] Csanady G T. Motions in a Model Great Lake Due to a Suddenly Imposed Wind[J]. Journal of Geophysical Research, 1968, 73(20): 6435-6447.

[25] Birchfield G E. Response of A Circular Model Great Lake to A Suddenly Imposed Wind Stress[J]. Journal of Geophysical Research, 1969, 74(23): 5547-5554.

[26] Large W G, Pond S. Open Ocean Momentum Flux Measurements in Moderate to Strong Winds [J]. Journal of Physical Oceanography, 1981, 11(3): 324-406.

Algorithm Improvement of Barotropic Force in an Unstructured Grid Two-dimensional Ocean Model

CHEN Bing-rui
（East China Sea Marine Forecasting Center, State Oceanic Administration, Shanghai 201206 China）

Abstract：Based on A2D, an independently developed unstructured grid two-dimensional ocean model, and by comparison with analytical solution under an ideal Model Great Lake, the algorithm of barotropic force was improved via employing different computational designs. In the improved algorithm, an assistant design was introduced to cooperate with the major design in order to acquire better stability. Results of elevation field and site time series showed that algorithms of three experiments got satisfactory accuracy and better stability. The TSNS algorithm was among the three, in which the major design solves elevation at centroid and velocity at mid-point of side, and the assistant design solves elevation at node and velocity at mid-point of side. Due to the similarity of its major design to C-grid design in a structured grid, the TSNS algorithm had advantages in conservation and movable tide-flat boundary treatments, which made the algorithm easier to apply for real ocean simulations. The TSNS algorithm was applied to the simulation in real sea near Jiangsu and Zhejiang. The simulated elevation had a good agreement with observed data, and was similar with results from the original A2D model. The improvement in stability will help TSNS algorithm get better stability in three-dimensional upgraded version, and will be the foundation of a successful three-dimensional version in future.

Key words：unstructured grid; ocean model; analytical solution; barotropic force; stability

作者簡介：陳昞睿（1980-），男，工程師，博士，從事海洋數(shù)值預(yù)報研究。E-mail：cu238@163.com

基金項目：國家海洋局青年海洋科學(xué)基金（2013202）。

收稿日期：2015-07-20

中圖分類號：P731

文獻(xiàn)標(biāo)識碼：A

文章編號：1003-0239（2016）01-0027-10