本站小編為你精心準備了設計正壓梯度力的計算參考范文,愿這些范文能點燃您思維的火花,激發您的寫作靈感。歡迎深入閱讀并收藏。
《海洋預報雜志》2016年第3期
摘要:
基于一套自主研制的無結構網格二維河口海洋數值模式A2D,在大圓湖理想模型下,通過與解析解進行比較分析,采用不同架構配置,改進設計正壓梯度力計算方法。改進后的算法中引入了從算架構的配置,以配合主算架構,得到更佳的穩定性。通過水位場平面分布與單點過程線可以發現,三組試驗的算法均獲得了較好的精度和比原算法更好的穩定性,其中TSNS配置算法(中心點計算水位、邊中點計算流速的主算架構,配合節點計算水位、邊中點計算流速的從算架構)由于其主算架構更接近結構網格下的C網格,在守恒性、移動潮灘邊界處理等方面具有一定優勢和便利性,有利于在實際海洋中的計算。將TSNS配置算法在江浙沿海進行試算,水位驗證結果與實測基本符合,與原A2D模式計算水位之間無顯著差異。TSNS算法在穩定性方面的改進,有助于提升模式升級為三維后的穩定性,為今后模式成功升級為三維打下基礎。
關鍵詞:
無結構網格;海洋模式;解析解;正壓梯度力;穩定性
1引言
無結構網格憑借其易擬合岸界、可局部加密等優勢,正越來越多地被用于海洋水動力數值模式。當前國際上知名的模式多來自國外[1-12],其中采用無結構網格的模式有FVCOM[10-11]、ADCIRC[12]等,這些模式在國內被廣泛使用,而國內自主開發的模式[13-14]則非常缺乏。鑒于此,2006—2011年期間,筆者所在研究組自主研制了一套無結構網格二維河口海洋數值模式[15-16],其中2011年最終版本后文稱為A2D[16]。A2D模式采用無結構三角形網格,基于有限體積法求解。水位在網格中心(重心)求解,流速在網格邊中點求解(同時求解x和y方向流速)。這種配置結合了各種已有的無結構網格模式的長處。一方面,它與國際上較為通行的C網格[17-18]配置相近,具有在網格中心求解水位之利于干濕判別的優點;另一方面,直接求解x-方向和y-方向流速使得離散簡潔高效,無須水平坐標轉換,且能更方便地利用有限體積法,提升守恒性。A2D模式在時間上顯式求解,采用預估修正法[19-22]。通過理想試驗、黃浦江和長江口等環境下的試驗計算,有效地驗證了A2D模式的精度[16]。但作為新生事物,A2D模式尚不成熟,存在一些不足,如模式還只是二維、在理想試驗下有微小數值波動等。筆者所在研究組曾將A2D模式以相同架構和算法升級至三維版本(稱之為A3D)[23],但在對長江口海域的計算中發現A3D的水動力不如A2D穩定,所以三維版本暫未獲得成功。A3D的不穩定性,很有可能源自于A2D的架構和算法,其中求解流速時正壓梯度力項為主要項,是關注的重點。本文對A2D的正壓梯度力算法進行分析與改進,提升模式在一個理想的大圓湖解析解模型下的穩定性表現,并投入實際海洋中試算。
2大圓湖理想模型簡介
Csanady提出了一個理想的大圓湖模型[24],計算了該圓湖在風的作用下的表明波動。Birchfield指出了其中一些錯誤并給出了正確的解析解表達式[25]。該模型為一個半徑67.5km,水深75m的平底大圓湖(見圖1),初始時刻水位靜止為0,施加以恒定西風3m/s,科氏力系數f取常數0.0001,忽略底摩擦。在風與科氏力的共同作用下,湖表將會產生水位波動,并且可以得到該波動的解析解表達式。通過這個模型,可以將數值模式的數值解與解析解進行比較,從而探討模式的精度。由于這個大圓湖的空間尺度很大,根據量綱分析的結果,流速平流項和擴散項對水位和流速變化的貢獻極其微小,可忽略不計[24-25],所以影響模式計算精度的主要項為正壓梯度力項,因此A3D在此模型下的穩定性問題,應是緣于正壓梯度力項。重新設計模式的算法,使得正壓梯度力項得到更合理地求解,并在大圓湖理想模型下取得較好的結果,是本文的主要目標。
3原模式的穩定性分析
采用A2D和A3D分別對大圓湖進行模擬(網格見圖1),時間步長均為5s,并在A點輸出站位過程線作比較。模式A2D與解析解的過程線比較結果如圖2所示。可以看到,模式的計算結果與解析解十分吻合,且總體比較光滑,沒有出現明顯的不穩定,說明模式的水動力達到了較高的精度。但在4.89d左右,水位的數值解(紅色)存在微小的波動,即存在微小的不穩定跡象。而在相同架構和算法的三維模式A3D下,水位的不穩定性表現得更為明顯(圖3)。對比解析解(見圖4)和A2D(見圖5)的水位場分布也可以發現,在模式5d時,A2D的計算結果的等值線略顯抖動,這也反映出模式在穩定性上略有欠缺。A2D在理想試驗下的微小不穩定,在實際海洋中計算時并不會發生,因為無頻散的流速平流項會使得模式保持穩定。但當三維開發版本A3D模式在實際海洋中計算時,由于不穩定性增大,所以難免出現個別區域流速或流向異常,影響計算結果的正確性。所以,本文改進A2D正壓梯度力項算法,使得模式更為穩定,是模式升級為穩定的三維版本的基礎。
4算法改進和探討
4.1原模式控制方程組A2D模式包含水動力計算和鹽度計算,而溫度暫時不作計算,取為常數T=10°C。本文不牽涉溫鹽算法的改進,主要對水動力控制方程組的相關部分作簡介。對流體不可壓縮、Boussinesq和靜力近似下的海洋動力學原始方程組作垂向積分,可得到垂向平均的二維控制方程組:式中:t為時間,ζ表示海表水位,D代表總水深(總水深D=H+ζ,H為固定不變的基準水深),U=1D∫-Hζudz和V=1D∫-Hζvdz分別表示水平x方向和y方向的垂向平均流速(其中u和v分別為空間中某點的水平x方向和y方向的局地流速),Fx和Fy為x方向和y方向的流速水平擴散,AM為水平湍流粘滯系數,f為柯氏力系數,g為重力加速度,<wu(0)>和<wv(0)>為海表應力,<wu(-1)>和<wv(-1)>為海底摩擦應力,ρ為密度,ρ0為參考密度。其中密度ρ=ρ(S,T,p)為鹽度S、溫度T(°C)和海壓p(Pa)的函數,采用1980年國際海水狀態方程(EOS80)計算。海表應力<wu(0)>和<wv(0)>為風應力在x和y方向上的分量:式中:ρa為空氣密度;Ua和Va分別為x和y方向的風速,||Va=U2a+V2a為風矢量Va的絕對值大小;CD為海水對風的拖曳系數,它根據Large和Pond改進的穩定狀態拖曳系數計算[26]:海底摩擦應力<wu(-1)>和<wv(-1)>根據謝才公式得到:式中:nb為曼寧系數,一般取值為0.015到0.018之間。
4.2原模式正壓梯度力算法分析在已有的A2D模式中,正壓梯度力項的計算方法如圖6所示。構造至多4個控制體A1、A2、A3和A4(未必一定有4個,網格邊緣和干點附近會有缺失),每個控制體的3個頂點均有計算得到的水位,于是可以通過格林公式計算控制體內的水位梯度。再將各控制體的水位梯度平均,得到邊j上的水位梯度,從而計算得到邊j上的正壓梯度力結果。觀察此算法,較為顯著的問題為未能利用到節點水位進行正壓梯度力計算,控制體遠端的網格點距離較遠,不利于數值穩定。而如果對水位進行插值,則由于水位梯度計算本來就對空間位置敏感,會導致精度不高,這在A2D研發階段已做過嘗試,其效果不如A2D最終方案理想[16]。在當前配置架構下暫時無法找到更好算法的情況下,嘗試更多不同的配置架構是較好的途徑。簡便起見,將不同配置架構的試驗采用其配置特征命名,分別用N、S、T3個字母代表三角形網格的節點、邊中點、網格中心,試驗的前兩位字母代表水位計算點和流速計算點的位置。如原A2D模式的命名為:TS。有些試驗采用一種配置架構作為主算架構,另一種配置架構作為從算架構,從算架構作為第2套同時計算的輔助配置架構為主算架構提升穩定性,這些試驗的命名中第1—2位字母代表了主算架構的配置,第3—4位字母代表了從算架構的配置,如TSNS即代表網格中心水位、邊中點流速的主算架構結合節點水位、邊中點流速的從算架構的試驗。在經過對多種不同配置算法的嘗試后發現,部分配置架構下得到了較穩定的結果,其中包括TSNS。TSNS在節點計算輔助的水位,有助于提升模式的穩定性。這種配置的主算架構仍然是中心 位、邊中點流速,但增加計算節點水位、邊中點流速的從算架構,使得中心與節點同時具有水位,更有利于邊中點位置的正壓梯度力求解。
4.3TSNS配置算法和計算結果在TSNS配置的算法中,通過連續方程,利用邊中點的流速同時求出網格中心點和網格節點的水位,并通過動量方程,利用網格中心點和網格節點的水位求出邊中點的流速。求網格中心點i的水位時,根據連續方程(1),仍然沿用原A2D的求法,在三角形控制體A(圖7a)中根據格林公式求得:式中:l為繞A一周的正向曲線,即l的方向為逆時針,A始終在其左邊。求網格節點m的水位時,也以格林公式(13)求解,不過控制體A變為包圍節點m的多邊形(圖7b)。控制體A有可能完全包圍節點m,也有可能如圖7b一般缺失部分角度,無論哪種情形,l均為繞A一周的正向曲線。在圖7b的情形下,僅需計算j4-i3-j3-i2-j2-i1-j1上的線積分,因為j1-m-j4上線積分等于0無通量進出。每一小段的流速為該小段所連接邊中點的流速。求邊中點j的流速時,對正壓梯度力項的算法作了修改。在大圓湖模型中,無斜壓梯度力項,平流項也非常小可忽略不計,故對精度影響最大的項正是正壓梯度力項gD∂ζ∂x和gD∂ζ∂y。TSNS配置算法中,動量方程僅正壓梯度力項較原A2D模式有改變,其它項均維持不變。TSNS配置算法在計算邊中點j的流速時,通過格林公式解出正壓梯度力項中的∂ζ∂x和∂ζ∂y,其中控制體A取為連接節點和中心點的菱形區域(見圖8),l為包圍A的正向曲線。由于節點和中心點的水位值模式均有計算,故沿i2-m2-i1-m1-i2的線積分易于求得。由于在以上計算方法中,同時計算了網格中心與網格節點兩套水位,在時間積分久后兩套的數值解可能存在分離的現象,所以將從算架構中的節點水位按照一定速度向主算架構中的中心點水位回歸。這里取每時間步0.05的回歸比例進行趨近。將上述TSNS算法在大圓湖模型中試驗,時間步長仍取5s,發現TSNS的精度與原A2D接近,穩定性更好(見圖9),第4.89d未出現不穩定抖動,水位場等值線也更平整。
4.4其它配置的計算結果除了TSNS配置外,本文還測試了一些其它配置,其中NSTS配置和ST配置這兩組試驗也得到了不錯的結果。在NSTS配置中,主算與從算架構與TSNS對換,NS主算、TS從算,回歸系數仍取0.05。其結果與TSNS略有不同(見圖10),精度接近TSNS,穩定性較好。而ST配置下也得到了穩定性尚可的結果(見圖11)。通過統計點A水位時間序列的平均誤差和均方根誤差來比較這3組試驗與原A2D模式的精度,可以發現TSNS配置的精度在所有4種配置中最好(見表1),同時誤差在模式4-5d時較之前時間段有所增大。由于模式最終應用時,以TS作為主算架構的TSNS配置更接近結構網格下的C網格,在守恒性、移動潮灘邊界處理等方面具有一定優勢和便利性,而NSTS配置和ST配置在各種邊界條件設計中存在一定的困難,故NSTS配置和ST配置不作重點介紹,具體算法細節在此省略,僅展示理想試驗下的結果,后續在實際海洋中的試算僅基于TSNS配置下進行。
5江浙沿海海域試算
5.1模式網格和基本設置基于改進了正壓梯度力算法的TSNS配置模式和原A2D模式,對東海區范圍內包含呂四測站、嵊山測站和定海測站的江浙沿海海域進行試算,以觀察模擬效果,并進行水位驗證,其中TSNS配置模式因配置變化對邊界條件進行了部分完善。模式的網格范圍見圖12,基于54坐標,包括了整個長江口、杭州灣和鄰近海區。東邊至124.5°E附近,北邊到34.3°N左右,南邊到28.4°N左右,長江上游邊界取在大通。長江口內、深水航道附近和島嶼附近的網格作了局部加密,最小網格分辨率可達100m,而口外網格則被放大,最大超過10000m。模式中深水航道的導堤和丁壩漲潮時淹沒、落潮時露出,是作為動邊界處理的。模式的時間步長統一取為1s。外海開邊界處的水位利用16個分潮(M2,S2,N2,K2,K1,O1,P1,Q1,MU2,NU2,T2,L2,2N2,J1,M1,OO1)的調和常數計算得出,而在長江口上游大通處則利用實測徑流量資料給出通量邊界條件。海表面的風場以每6h為1組、分辨率為0.5°×0.5°經緯度的氣象預報后處理結果給出底摩擦曼寧系數在模式中設為0.015。模式從2008年11月5日起計算,共計算30d。通過搜集的呂四站、嵊山站和定海站的實測水位資料,對原A2D模式和改進后的TSNS配置模式進行比對驗證。
5.2水位驗證結果模式共運行30d,輸出第10—30d的水位,與呂四站、嵊山站和定海站的實測資料進行比對(圖13—15,圖中上子圖的藍實線為原A2D計算水位,黑虛線為TSNS配置計算水位,紅點為實測水位;下子圖的黑實線為TSNS計算水位減去A2D計算水位的差,藍實線為0軸)。可以發現,TSNS配置模式的水位計算結果與原A2D模式(TS配置)的水位計算結果非常接近,兩者差異在3站均不到0.1m,同時模式計算結果與實測基本符合,基本反映了3個測站的水位變化規律,其中嵊山站和定海站存在較明顯的潮汐日不等現象,而定海站的振幅略偏大。由于此番改進主要目的是探討是模式本身的性能,故在江浙沿海海域進行試算并驗證時,底摩擦曼寧系數統一取了0.015,未針對不同海域進行分區設置,對驗證效果略有影響。而TSNS配置模式與原A2D模式雖然在算法上存在差異,但由于案例設定的物理環境相同,兩者受相同的邊界條件驅動,故兩者的水位計算結果與變化特征并沒有出現顯著差異。
6結論
本文通過對模式配置架構的改變,重新設計正壓梯度力項的算法,得到了3種比原二維A2D模式穩定性更高的算法,在大圓湖理想模型下模擬得到了較穩定的結果,同時精度與原A2D接近,均與解析解較為符合。其中TSNS配置算法中由于其主算架構TS配置更接近結構網格下的C網格,在守恒性、移動潮灘邊界處理等方面較好,故在完善了相應的邊界條件后,在江浙沿海海域進行了試算并與實測資料進行了對比驗證。水位驗證結果與實測基本符合,同時與原A2D模式計算水位之間無顯著 差異。原A2D的算法雖然與TSNS算法一樣也能在二維情形的實際海洋中穩定,且計算結果相近,但A2D的三維開發版本A3D在實際海洋計算中卻遭遇穩定性不佳的問題。大圓湖理想模型下更穩定的表現,有助于提升模式升級為三維后的穩定性,所以TSNS算法在穩定性方面的改進,為今后模式成功升級為三維打下基礎。經過改進后的TSNS配置模式可在風暴潮模擬等場合進行應用,同時由于其網格配置特性,可進一步優化邊界條件設定,根據需要靈活設置堤壩等特殊情形,具有一定的發展前景。
作者:陳昞睿 單位:國家海洋局東海預報中心