楊 飛，楊文俊，楊 森

（1.長江科學院水利部江湖治理與防洪重點實驗室，武漢 430010；2.重慶市武隆縣水務局，重慶 408500）

ELCIRC源程序代碼分析

楊 飛1，楊文俊1，楊 森2

（1.長江科學院水利部江湖治理與防洪重點實驗室，武漢 430010；2.重慶市武隆縣水務局，重慶 408500）

開放源程序代碼ELCIRC是采用基于水平向無結(jié)構(gòu)網(wǎng)格、垂向z坐標體系和半隱格式的歐拉－拉格朗日有限體積／有限差分方法來解淺水方程。研究了ELCIRC中除控制方程組求解以外的部分，分析了源程序在插值計算、拓撲結(jié)構(gòu)、定解條件、分層信息和干濕法等技術(shù)上的具體操作，尤其是對歐拉拉格朗日模型特有的逆向追蹤算法做了詳細闡述。同時討論了一些細節(jié)問題，在不恰當?shù)牡胤浇o出一些參考和建議，避免小數(shù)做分母、相近數(shù)相減引起的較大誤差。然而，這些經(jīng)過實際應用的算法原理，可以為我們開發(fā)水動力學模型提供有價值的參考。

ELCIRC；逆向追蹤；源代碼；三維水動力學模型；算法

1 研究背景

ELCIRC［1］模型（Eulerian-Lagrangian CIRCulation）是美國俄勒岡健康與科學大學海岸陸地邊緣陸緣研究中心（Oregon Health＆Science University OHSU’s Center for Coastal and Land-Margin Research）基于對哥倫比亞河的研究，開發(fā)的哥倫比亞環(huán)流模型CORIE的一部分，可以單獨用于模擬計算。模型已進行基準測試，并應用到哥倫比亞河的物理過程的模擬。國內(nèi)外應用逐漸增加，Myers和Aikman2003年研究St.Johns河流域洪泛區(qū)淹沒過程［2］。楊金艷等2006年利用ELCIRC計算了長江口的潮流場［3］，該模型的中文介紹可以參考楊金艷［3］碩士畢業(yè)論文。Gong等2008年將ELCIRC應用到海南新村瀉湖研究中［4］，Gong等2009年使用ELCIRC檢驗了Chesapeake海灣水位對選定的3組東北風響應及正壓模式下潮下帶海灣河口水交換［5］。Dias等2009年利用水動力ELCIRC和輸移VELA，VELApart模型在葡萄牙福爾摩挲河口瀉湖的應用，研究入口位置改變的影響［6］。Picado等2010年以葡萄牙阿威羅瀉湖為例研究了局部地形調(diào)整對潮汐可能產(chǎn)生的影響，采用ELCIRC模型研究在過水面積增加、瀉湖的水動力條件減弱下可能的潮汐變化［7］。

ELCIRC采用了基于水平向無結(jié)構(gòu)網(wǎng)格、垂向z坐標體系下，半隱格式的歐拉－拉格朗日有限體積／有限差分方法來解淺水方程。算法上考慮了多種紊流閉合模式，也包括了潮汐勢，大氣壓梯度項及水氣交換，可以模擬多種物理過程［1］。圖1為以3D斜壓哥倫比亞河環(huán)流模擬為例的ELCIRC求解過程結(jié)構(gòu)圖，網(wǎng)格水平節(jié)點33 634個，水平單元50 389個，垂向分層62個，對應總的活動柱面約為2.2×1012個［1］。本文主要涉及的是動量方程組求解以外的參數(shù)處理等，文中大部分程序語言來自文獻［8］中的ELCIRC5＿4c。

2 拓撲關(guān)系

程序首先使用平面的單元、邊、頂點來表示實際區(qū)域，然后進行垂向劃分。實際的計算三維網(wǎng)格為柱體，實際中用到的變量定義位置如圖2。字母t和s（或tem和sal）表示溫度和鹽度，q2和xl表示紊流模型中的紊動能和混合長，字母p（或nd）表示節(jié)點處，字母e（或el）和c表示單元中心，字母s（或sd）表示側(cè)面中心，字母z表示垂向分層，0表示初始，diff表示擴散，tmp表示過程量，half表示相鄰垂向分層平均，tot表示垂向累積，rho表示密度，bt表示逆向追蹤。對于流速，下標n和t表示法向和切向，下標u，v和w表示節(jié)點處流速，后面的數(shù)字1和2對應迭代計算過程不同時段。圖中d表示厚度，但程序中多數(shù)變量會將d作為首字母，表示雙精度數(shù)值，將n作為首字母表示整數(shù)值，如單元、節(jié)點和邊的序號。

圖1 ELCIRC計算流程圖［1］Fig.1 Flow chart of ELCIRC［1］

圖2 程序常用變量定義位置Fig.2 Position of program variables in grid

2.1 平面位置關(guān)系數(shù)組

nx（4，4，3），考慮三角形和四邊形2種單元，見式

（1）。數(shù)組實現(xiàn)循環(huán)三角形1－2－3循環(huán)和四邊形的

1－2－3－4循環(huán)。注意，在平面網(wǎng)格生成文件hgrid.in中，單元對應的節(jié)點就是按照空間位置逆時針的順序表述，數(shù)組可以循環(huán)尋址，起點任意。實際中調(diào)用非零元素來得出邊與節(jié)點在單元中的相對位置關(guān)系。

2.2 邊正向判定

離散網(wǎng)格的每條邊其實是矢量，有（正）方向，端點定義為起止節(jié)點，邊的正向和單元編號關(guān)系緊密。一般情況下單元編號是按照空間位置由下到上、由左到右順序編排的，這給邊的定義帶來極大的方便。對于內(nèi)部邊，對應2個相鄰單元，邊的法線方向為由編號低的垂直邊指向編號高的；對于組成區(qū)域邊界的邊，其只對應一個單元，即邊界單元，邊的法線方向為由單元內(nèi)部垂直邊指向外。邊的切向均為法線方向逆時針旋轉(zhuǎn)90度。假定邊為j，存儲時采用x軸正向到邊夾角的正弦值snx（j），余弦值sny（j），邊的法向單位向量表示為（snx（j），sny（j）），切向單位向量為（－sny（j），snx（j））。

在計算單元i的水位時，其中j邊的法向正負符號Ssign（i，j）表達式具體為

3 逆向追蹤

逆向追蹤從邊的中點（側(cè)面中心x0，y0，z0）出發(fā)，根據(jù)該點t時刻的速度，顯式倒推判斷Δt以前該質(zhì)點位置（xt，yt，zt）。網(wǎng)格無水和流速小于10－6m／s，不需要逆向追蹤計算，垂向流速為側(cè)面四頂點垂向流速平均值。由于方程建立在邊上，每條邊都要考慮逆向追蹤方程，上面是依據(jù)邊來搜尋逆向追蹤的需要網(wǎng)格，邊界以外的邊都對應2個單元，一般優(yōu)先考慮編號小的單元。

3.1 定位（quick-search）

（x0，y0）在nel單元內(nèi)，倒推得到的（xt，yt）也在nel內(nèi)時，確定zt所在層jlev即可，但不能溢出邊界；（xt，yt）不在nel內(nèi)時，找出nel與（x0，y0）、（xt，yt）連線（這時（x0，y0）可以先微調(diào)內(nèi)移，使其不在邊上）相交的一邊及交點（用到下面的intersect2函數(shù)），求出交點的高程，并用交點代替（x0，y0）作為起點，共用該邊的單元作為新的nel單元，繼續(xù)迭代直到（x0，y0）、（xt，yt）在同一單元內(nèi)。（x0，y0）的速度值是取所屬邊速度（兩端點速度均值，求解動量方程不考慮時間因素，流場過程知道時求解物質(zhì)輸移方程時要考慮速度隨時間變化），這樣簡化默認相鄰節(jié)點間的速度差異在誤差允許范圍內(nèi)。筆者認為可以根據(jù)端點線性插值交點速度，保證一階精度。每次確定新起點后，（xt，yt）根據(jù)剩余時間和起點速度重新確立。整個過程如圖3所示，該圖時間步長dtb不同于文獻1中圖4所示整個時間步長dt。從t到t－dtb的迭代次數(shù)不能預先知道，因此要限定以避免進入死循環(huán)。

圖3 一次逆向追蹤迭代示意圖Fig.3 Back-track iteration for one sub step

整個過程中并沒有考慮z向速度的逆向追蹤，這與動量積分方程建立在水平法向和切向是一致的。逆向追蹤遇到邊界時，只取切向速度。時間步長Δt劃分越小，誤差就越小，因此把原來的時間步長在計算有效性（整個程序的計算時間長短和計算量的大?。┣疤嵯逻m當細分ndelt個子步長dtb，進行迭代ndelt次逆向追蹤迭代（dtb＝Δt／ndelt），每次迭代得到的（xt，yt）作為下次的起點。

3.2 交叉點（intersect2）

判斷端點為（x1，y1）、（x2，y2）和（x3，y3）、（x4，y4）兩線段有無交點，若有交點，則存在0－1之間的2個數(shù)tt1，tt2，滿足：

判斷tt1，tt2是否在0和1之間，即可確定有無交點，若有，則可求出交點。

3.3 逆向追蹤（back-track）

利用quick-search函數(shù)找到（x0，y0）對應的（xt，yt）以及nel后，確定（xt，yt）所在的平面，垂向插值nel單元在該平面頂點值，然后根據(jù)在單元的相對位置，平面插值（xt，yt）。

單元為三角形時，（xt，yt）與頂點連線將其分為3個小三角形，（xt，yt）的值為3頂點值的加權(quán)平均，權(quán)重為頂點所對應的小三角形與單元面積之比。單元為四邊形時，需要雙線性投影變換，中心對應新的坐標原點，原來的點1，2，3，4或點A，B，C，D映射坐標為（－1，－1），（1，－1），（1，1），（－1，1），如圖4。分別討論平行四邊形、梯形、不規(guī)則四邊形，任意一點P的投影變換后的坐標（ξ，η）。

圖4 四邊形單元雙線性投影Fig.4 Bilinear projection of quadrilateral unit

平行四邊形時，如圖5（a）所示。

式中：X，E分別是2條邊的中點，O為平行四邊形的中心，（xt，yt）即為P點。

梯形時，以其中一組對邊AB和DC平行為例，另一種情況一樣。對邊中點連線將梯形分為4個區(qū)域。如圖5（b），M，N為對角線BD和AC的中點，E，F(xiàn)，G，H為AB，BC，CD，DA的中點，O為HF和EG的交點。SP平行底邊，延長線交BC于I，ξ大小為SP／SI，η大小為OS／OG。

圖5 圖形變換Fig.5 Graphic transformation

一般對不規(guī)則四邊形，如圖5（c）所示，其算法為

式中：α＝2 H→F×2M→N；β＝2 H→F×2→EG－4 O→P×2 M→N；γ＝－4 O→P×2→EG。求解一元二次方程式（9）得到ξ的2個根，選擇絕對值小于1的根，代入式（10）中求解η。

3.4 算法誤差

程序原有的計算不規(guī)則四邊形（圖5（c））使用到向量M→N和→AD×C→B（記為d）。數(shù)值計算要求數(shù)

eta值穩(wěn)定，在設計算法時還應盡量避免誤差危害，防止有效數(shù)字損失，要避免兩相近數(shù)相減和用絕對值很小的數(shù)做除數(shù)，還要注意運算次序和減少運算次數(shù)。在程序判斷語句中使用了M→N的分量數(shù)值和d是否

eta大于小量small3＝10－5作為判斷依據(jù)，對于動量方程求解器計算精度為7位有效數(shù)字的程序來說，繼續(xù)將其作為計算乘數(shù)和除數(shù)，筆者認為會帶來不可忽略的誤差。考慮到在網(wǎng)格剖分基本優(yōu)化為正交網(wǎng)格，將會有更多的單元格處在臨界條件附近，盡管采用雙倍字長計算，相近數(shù)相減仍會帶來有效位數(shù)的減少。因此需要對算法進行修改。

先利用式（7）求出近似值ξ0，令ξ＝ξ0＋Δξ，Δξ是一個未知校正量，然后將ξ0＋Δξ代入式（9）得到式（11），參照式（9）表示成式（12）。

這樣由（13）或迭代計算式（14）求解會使得計算的精度有所改善。對求得的小量Δξ加以判斷，Δξ太大（如大于0.1）則不采用。

2

2

2

程序中包含判斷ξ＞1.1，η＞2.0是否成立的語句，顯示出了求解結(jié)果的誤差，當然，由于速度梯度不大，這種求解誤差給結(jié)果帶來的影響不是很明顯，但是較大的計算量和不可忽略的誤差確實暴露出這種插值方法的缺陷，現(xiàn)實中存在一些較為成熟的插值算法。為了追求更高的計算精度和求解可靠性，有必要詳細根據(jù)流體特征分析速度的插值算法。

最后，由所得的ξ，η來確定加權(quán)系數(shù)，進行插值。

對于物質(zhì)輸移的逆向追蹤，要將單元細分為4個子區(qū)，四邊形是采用對邊中點連線劃分成4個子四邊形，三角形則是3條中位線劃分成4個子三角形，判斷逆向追蹤后質(zhì)點所屬子區(qū)，在子區(qū)內(nèi)部進行插值計算（圖6、圖7）。圖6，圖7中areas（1），areas（2），areas（3），areas（4）表示所在子三角形面積，作為插值計算的權(quán)重。

圖6 物質(zhì)輸移的四邊形插值Fig.6 Interpolation on quadrilateral unit for transport equation

圖7 物質(zhì)輸移的三角形插值Fig.7 Interpolation on triangle unit for transport equation

4 定解條件

初始場分布條件是在節(jié)點處給定的，溫度和鹽度都要進行邊、單元初始化賦值。陸地邊界和開邊界在計算中有很大的差別，陸地邊界實際限定了流場的范圍，干濕判別標準是根據(jù)水深是否大于某一臨界值確定。在高出水面的陸地區(qū)域，也要參與網(wǎng)格剖分，所有剖分的網(wǎng)格都會參與計算；垂向網(wǎng)格剖分時，水體以下部分的柱體也參與計算，相同水平位置的無水柱體，均賦予同一水平位置上最底部有水單元的值，這樣處理是為了網(wǎng)格整齊，方便調(diào)用求解。

對于水位，一般會有全區(qū)域的單一初始值，在計算過程中會逐漸調(diào)整為合理值。恒定水流求解問題水流只需給出流量的上邊界（進口）條件和水位的下邊界（出口）條件，鹽度和溫度場則根據(jù)需要給出進口處的變化過程，非恒定水流求解問題要給出流量的上邊界過程和水位的下邊界過程。要將給出的進口處流量轉(zhuǎn)化為流速，只是簡單的利用開邊界寬度來除流量，得到單寬流量，再除以邊界水深，因此研究區(qū)域和進口邊界有相當距離，避免失真流速的影響。

5 網(wǎng)格前處理

在預先分層的網(wǎng)格上確定實際水體的頂?shù)讓樱ㄔ闯绦驑俗⒅蠱表示頂層號，m表示底層號），初始化單元、節(jié)點、邊的分層信息。每次迭代求解出水位值后，都要重新計算單元、節(jié)點、邊的分層信息。以單元為例進行說明，實際水深小于某一特定值時，單元為無水單元，頂層號kfe＝0；否則計算單元新的底層號kbe，頂層號kfe，各層厚度dc和上下鄰層平均厚度dchalf。如圖8，ztot指的是層面的垂直坐標，nvrt個層對應nvrt＋1個層面，jlev表示層號。逆向追蹤中用到的zt表示某質(zhì)點的實際高度，zup指該點到所在層層頂?shù)木嚯x。頂?shù)淄瑢訒r，dc＝dchalf＝水深；不同層時頂?shù)讓雍馾c和dchalf，如圖9所示。

圖9 頂?shù)讓邮疽鈭DFig.9 Schem atic diagram of top and bottom elements

圖8 垂向分層示意圖Fig.8 Schematic diagram of verticalmeshes

dzhalf，dz要在動量方程做分母（dzhalf（i，kbs－1）、dzhalf（i，kfs）不做分母，kbs，kfs表示邊的底層號和頂層號），所以要檢查兩者是否小于某限定值denom（略小于判定網(wǎng)格干濕時的臨界水深），初始水位線應盡量靠近分層的中上部，作為簡化，程序中的初始水位設定要求水位線不能靠近層面線。迭代過程中，這種檢查只針對頂?shù)淄瑢?。ug，vg表示j邊上的第k層面上的流速可通過式（15）、（16）線性插值得到

圖10 四邊形單元節(jié)點插值Fig.10 Interpolation w ith nodes on quadrilateral unit

其中utmp（k），vtmp（k）在柱體第k層側(cè)面（平面網(wǎng)格中仍對應j邊）面心上的橫向和縱向流速。

如圖10所示根據(jù)距離倒數(shù)加權(quán)平均求得4條邊nd1，nd2，nd3，nd4的公共節(jié)點i的x方向流速uu2（i，k）為

式中u1，u2，u3，u4和d1，d2，d3，d4分別對應nd1，nd2，nd3，nd4四邊中點x方向流速和四邊長度。Y方向流速vv2（i，k）為

v

g（nd1，k），vg（nd2，k），vg（nd3，k），vg（nd4，k）分別對應nd1，nd2，nd3，nd4四邊中點y方向流速，distj（nd1），distj（nd2），distj（nd3），distj（nd4），即式（17）中的d1，d2，d3，d4，表示所在邊的長度。節(jié)點垂向流速ww2（i，k）則是按周圍4個單元垂向流速面積加權(quán)平均求得。

其中areas（1），areas（2），areas（3），areas（4）為4單元面積，we1，we2，we3，we4為4單元垂向流速。節(jié)點垂向流速則是按周圍4個單元垂向流速面積加權(quán)平均求得。

6 結(jié) 語

ELCIRC作為水動力學開放源代碼，在河口海岸數(shù)值模擬中取得一定的成功。文中給出了部分源代碼的具體分析，可以為水動力學模型開發(fā)人員提供方法實現(xiàn)的參考。源代碼凝結(jié)了美國俄勒岡健康與科學大學海岸陸地邊緣陸緣研究中心科研人員的辛勤勞動，在此表示感謝，筆者只是對其優(yōu)秀的程序進行部分粗略的學習，重點著眼于其核心部分的歐拉拉格朗日模型逆向追蹤過程，發(fā)現(xiàn)了其中的不完善的地方，希望能給源程序在改進過程中盡一點努力。ELCIRC的開發(fā)人員已經(jīng)在此代碼的基礎上進行改進，開發(fā)出SELFE程序，本文的大部分都可以為讀者學習SELFE提供幫助。

致謝：非常感謝清華大學研究生李健在ELCIRC源代碼學習上給予悉心指導和莫大的幫助！

［1］ ZHANG Y L，BAPTISTA A M，MYERSE P.A Crossscale Model for 3D Baroclinic Circulation in Estuary-Plume-Shelf Systems：I.Formulation and Skill Assessment［J］.Continental Shelf Research，2004，24（18）：2174－2214.

［2］ MYERSE P，AIKMAN F，ZHANG A J.A Forecast Circulation Model of the St.Johns River，F(xiàn)lorida［C］∥Proceedings of the Eighth International Conference on Estuarine and Coastal Modeling，November 3－5，2003， Monterey，California，USA：144－156.

［3］ 楊金艷.ELCIRC模型在長江口的應用［D］.南京：河海大學，2006.（YANG Jin-yan.Application of ELCRIC Model in Changjiang Estuary［D］.Nanjing：Hohai University，2006.（in Chinese））

［4］ GONG W P，SHEN J，WANG D R.Mean Water Level Setup／Setdown in the Inlet-Lagoon System Induced by Tidal Action：A Case Study of Xincun Inlet，Hainan Island in China［J］.Acta Oceanologica Sinica，2008，27（5）：63－80

［5］ GONG W P，SHEN J，CHO K H，et al.A Numerical Model Study of Barotropic Subtidal Water Exchange Between Estuary and Subestuaries（Tributaries）in the Chesapeake Bay During Northeaster Events［J］.Ocean Modelling，2009，26：170－189

［6］ DIAS JM，SOUSA M C，BERTIN X，et al.Numerical Modeling of the Impactof the Anc~ao Inlet Relocation（Ria Formosa，Portugal）［J］.Environmental Modelling＆Software，2009，24（6）：711－725.

［7］ PICADO A，DIAS JM，F(xiàn)ORTUNATO A B.Tidal Changes in Estuarine Systems Induced by LocalGeomorphologic Modifications［J］.Continental Shelf Research，2010，30（17）：1854－1864.

［8］ OHSU’s Center for Coastal and Land-Margin Research.ELCIRC Description：3D Baroclinic Circulation Model［EB／OL］.（2006－05－01）［2011－05－14］，http：／／www.stccmop.org／CORIE／software／elcirc／r ecent.html

（編輯：劉運飛）

Algorithm of ELCIRC Source Code

YANG Fei1，YANGWen-jun1，YANG Sen2
（1.Key Laboratory of River Regulation and Flood Control of MWR，Yangtze River Scientific Research Institute，Wuhan 430010，China；2.Water Affairs Bureau ofWulong County，Chongqing City，Chongqing 408500，China）

Open-source code ELCIRC（Eulerian-Lagrangian CIRCulation）solves the shallow water equations using a semi-implicit Eulerian-Lagrangian finite volume／finite difference method with horizontally unstructured grids and vertically unstretched z-coordinates.ELCIRC source code aside from the governing equation solution is analyzed in this paper.The operation of interpolation，topological structure，definite condition，hierarchy information，and wetting and dryingmethod are described in detail.Back-track，as the key technology of ELCIRC，is expounded comprehensively.Moreover，some detailed problems are discussed，and a few references and suggestions are given to avoid the errors caused by similar number subtraction and by employing small number as the denominator.ELCIRC has been applied in practice and could serve as a valuable reference for developing hydrodynamic models.

ELCIRC；back-track；source code；3-D hydrodynamic model；algorithm

TV148

A

1001－5485（2013）05－0097－06

2013，30（05）：97－102

10.3969／j.issn.1001－5485.2013.05.021

2012－04－25；

2012－06－13

國家自然科學基金面上項目（51079008）

楊 飛（1985－），男，山東泰安人，碩士研究生，從事水力學及河流動力學研究，（電話）13810609717（電子信箱）yangfeihaoyun＠163.com。