張慶河，王容基，李文俊，冉國全

(1. 天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072；2. 中交第三航務(wù)工程勘察設(shè)計院有限公司，上海 200032)

隨著計算機(jī)技術(shù)的發(fā)展，數(shù)值模擬方法在海洋動力環(huán)境的研究中得到了越來越廣泛的應(yīng)用．海洋動力如潮流等的模擬，一般通過有限差分法、有限元法和有限體積法等傳統(tǒng)空間離散方法，求解淺水方程而實現(xiàn)．近年來，間斷有限元(discontinuous Galerkin，DG)方法逐漸受到人們的關(guān)注．DG 的基本思想是在多項式空間內(nèi)尋找近似解，不強(qiáng)制物理量在相鄰單元的重合節(jié)點上相同，同時根據(jù)守恒性采用合適的數(shù)值通量實現(xiàn)相鄰單元的耦合[1]．DG 方法兼具有限體積法和有限元法的優(yōu)點，具有高階精度和局部守恒性，適用多種自適應(yīng)網(wǎng)格加密算法，同時還具有并行計算效率高的優(yōu)勢[2]．因此，利用間斷有限元方法開發(fā)新型海洋動力數(shù)值模型，成為近年來海洋與海岸動力數(shù)值模擬的一個重要方向．目前DG 海洋動力數(shù)值模擬方面比較典型的工作有：Dawson 等[3]將DG 方法應(yīng)用于大范圍潮流和風(fēng)暴潮模擬，源項和對流項采用標(biāo)準(zhǔn)的DG 格式，擴(kuò)散項則采用降階處理的局部DG 方法，通過求解球坐標(biāo)系下的二維淺水方程，模擬了強(qiáng)颶風(fēng)Ike 作用下的沿海潮波運動，計算結(jié)果表明該模型具有較高精度；Kubatko 等[4]建立了具有h-p 自適應(yīng)的DG 淺水模型，求解沿海復(fù)雜水體運動，模擬結(jié)果和ADCIRC 有限元模型進(jìn)行了對比，結(jié)果表明DG模擬結(jié)果更精確，并且在大規(guī)模計算時有更高的并行計算效率；Karna 等[5]在SLIM 間斷有限元水動力模型基礎(chǔ)上，提出一種全隱式干濕方法，保證了水深非負(fù)、質(zhì)量守恒和數(shù)值穩(wěn)定的性質(zhì)，模型在數(shù)值實驗和Scheldt 河口的實際應(yīng)用中均得到驗證．

除了上面提到的優(yōu)點外，間斷有限元方法也存在一定的缺點，如求解變量較多、計算量偏大等．為了減小計算量，李龍翔等[6]提出利用無積分節(jié)點間斷有限元求解對流擴(kuò)散方程，李文俊等[7]則進(jìn)一步在直角坐標(biāo)系下建立了無積分節(jié)點間斷有限元二維淺水方程模型．

在研究大范圍海域水流運動時，由于需考慮地球曲率和科氏加速度隨緯度的變化，直角坐標(biāo)下的二維淺水方程不再適用，通常需要求解球坐標(biāo)下的二維淺水方程．普通經(jīng)緯度球坐標(biāo)系下的數(shù)值模型往往存在兩方面的問題[8]：一是坐標(biāo)奇異，模型中某些項趨于無窮大；二是經(jīng)緯網(wǎng)格在極點附近網(wǎng)格輻合，通常有3 種解決方法，即改變格點劃分并轉(zhuǎn)換坐標(biāo)系統(tǒng)[9-10]、引入濾波器[11]和采用隱式時間積分方法[12]，第一種方法是比較可行的，從網(wǎng)格劃分上解決了網(wǎng)格輻合問題，通過轉(zhuǎn)化到平面坐標(biāo)避免了坐標(biāo)奇異，間斷有限元數(shù)值格式也容易實現(xiàn)．為此本文將采用極射赤面投影將球面坐標(biāo)轉(zhuǎn)換為平面坐標(biāo)的方法求解球坐標(biāo)下的控制方程，以建立可求解大范圍海域水流運動的無積分間斷有限元數(shù)值模型．

1 極射赤面投影坐標(biāo)系二維間斷有限元水動力模型的建立

1.1 極射赤面投影坐標(biāo)系二維淺水方程

原始球坐標(biāo)下的二維淺水方程為

式中：h 為總水深；t 為時間；η為相對靜水面的波動水深；R 為地球平均半徑；λ為經(jīng)度，λ∈[0,2π)；φ為緯度；u、v 分別為λ、φ方向上水深平均速度分量；f 為科氏力分量，f=2ωsinφ，ω為地球自轉(zhuǎn)角速度；K 為拖曳力系數(shù)，Cd為底摩阻系數(shù)；g 為重力加速度．動量方程(2)中右端第1 項為底坡源項，第2 項為科氏力項，第3 項為由徑向和緯度方向相對運動產(chǎn)生的附加項，第4 項為底摩阻項．

本文采用的極射赤面投影法[13]如圖1 所示，轉(zhuǎn)換后的正交坐標(biāo)x、y 和原經(jīng)緯坐標(biāo)系下的正交坐標(biāo)λ、φ的映射關(guān)系為

式中m 為映射參數(shù)，通過幾何相似得到，即圖1 中NA′ 與OB 的比值，即

其中α為常量，對北半球投影變換時α=1 ，對南半球投影時α=?1．

圖1 極射赤面投影法示意Fig.1 Sketch of stereographic projection method

為保證數(shù)值模型的和諧性質(zhì)，本模型采用Bollermann 等[14]的方法，對變換后的底坡源項進(jìn)行改造，使底坡源項與壓力項在靜水條件下平衡，即采用式(5)對其進(jìn)行拆分：

極射赤面投影坐標(biāo)系下的速度場分量U 、V 與原經(jīng)緯坐標(biāo)系下的速度分量u、v 的關(guān)系為

最后，將式(6)代入經(jīng)緯坐標(biāo)下的底摩阻項，可得到極射赤面投影坐標(biāo)下淺水方程的控制方程

其中，Cf為赤平投影坐標(biāo)下的科氏力系數(shù)，

1.2 無積分節(jié)點間斷有限元數(shù)值離散

采用節(jié)點間斷有限元方法進(jìn)行空間離散，令Ωk為求解域非重疊劃分后的第k 個單元， Ω?k為其邊界，n 為邊界處的單位外法向量．在定義計算域Ωk上最高不超過p 階的多項式空間Vp(Ωk)上，選擇基函數(shù) φi( x )組成的線性組合作為近似解逼近U，其中為基函數(shù)系數(shù)，使得殘差方程(8)達(dá)到最小．

在節(jié)點間斷有限元方法中，插值基函數(shù)和試驗函數(shù)均采用Lagrange 插值基函數(shù)定義，將近似解代入殘差方程連續(xù)積分計算后，得到方程強(qiáng)形式，為

通量項的體積分和面積分計算均采用無積分方法[6]，利用節(jié)點基函數(shù)插值性質(zhì)，將單元內(nèi)通量項用基函數(shù)近似為

由于物理通量F 在邊界是不連續(xù)的，有限體積法中成熟的數(shù)值通量方法可以容易地應(yīng)用于間斷有限元方法中，故引入法向通量數(shù)值通量F*[15]．進(jìn)而可得到未知系數(shù)Uk關(guān)于時間的常微分方程

隨后采用顯式四階五步龍格庫塔方法時間遞進(jìn)求解，最終得到不同時刻具有高階時間和空間精度的數(shù)值解．

1.3 邊界條件處理

采用虛擬邊界方法施加邊界條件，在可滑移固壁邊界中，水質(zhì)點沿邊界切向速度不變，法向流速為零，水位為零梯度，令分別代表虛擬邊界節(jié)點和邊界內(nèi)部節(jié)點的物理量，虛擬邊界物理量為

2 模型驗證

為了驗證本文建立的極射赤面投影坐標(biāo)系下二維間斷有限元數(shù)值模型的合理性，下面采用渤海海域兩個不同尺度的算例對潮波和潮流運動進(jìn)行模擬，并與實測數(shù)據(jù)進(jìn)行比較．

2.1 渤海海域潮波模擬

2.1.1 模擬區(qū)域及計算配置

渤海潮波運動數(shù)值模擬區(qū)域包括整個渤海海域，其 模 擬 范 圍 為：117.58° ～121.50°E ，37.14° ～40.91°N．圖2 顯示了渤海海域及模擬和實測潮汐分量進(jìn)行比較的渤海驗潮站分布情況．模型采用比三角形網(wǎng)格計算效率更高的非結(jié)構(gòu)化任意四邊形網(wǎng)格，利用有限元網(wǎng)格生成器Gmesh 生成，如圖3 所示，整個計算區(qū)域共有5 974 個四邊形單元和5 776 個網(wǎng)格節(jié)點，最大網(wǎng)格和最小網(wǎng)格尺度分別約為0.11°和0.02°，渤海水深數(shù)據(jù)源于 2005 年中國近海 1∶250 000 海圖．模擬時東邊界設(shè)為開邊界，其他為陸地邊界，開邊界潮位從Chinatide[16]潮汐系統(tǒng)提取．模擬區(qū)域底摩阻系數(shù)統(tǒng)一取0.002，暫不考慮河流及風(fēng)場強(qiáng)迫的影響．模型計算時間步長為10 s，模擬時間段為2006 年6 月29 日至2006 年8 月2 日，共35 d.

圖2 渤海驗潮站分布Fig.2 Tidal stations in the Bohai Sea

圖3 渤海計算域水平非結(jié)構(gòu)四邊形網(wǎng)格投影圖Fig.3 Distribution of the arbitrary quadrilateral meshes after being projected in the Bohai Sea

圖4 M2、K1 分潮模擬結(jié)果檢驗Fig.4 Model validation of M2 and K1 constituents

表1 M2 和K1 分潮實測模擬調(diào)和常數(shù)比較Tab.1 Comparison between simulated and measured results of the tidal harmonic constants of M2 and K1

2.1.2 調(diào)和常數(shù)驗證

根據(jù)準(zhǔn)調(diào)和分析方法，采用T-Tide[17]對30 d 的潮位計算結(jié)果進(jìn)行分析，得到計算區(qū)域8 個潮位站主要分潮M2和K1的調(diào)和常數(shù)，如表1 所示．圖4 為M2和K1分潮8 個驗證點模擬值與實測值的比較．模擬結(jié)果與渤海沿岸8 個驗潮站M2分潮和K1分潮觀測資料進(jìn)行了對比，得到M2分潮振幅絕對平均偏差為8.60 cm，相關(guān)系數(shù)為0.938，遲角絕對平均偏差為10.75°，相關(guān)系數(shù)為0.990；K1分潮振幅絕對平均偏差為1.78 cm，相關(guān)系數(shù)為0.975，遲角絕對平均偏差為10.50°，相關(guān)系數(shù)為0.997．總的來看，不同測站M2和K1兩個分潮的振幅和遲角與實測結(jié)果接近，模型對渤海潮波系統(tǒng)模擬結(jié)果較為準(zhǔn)確．個別站點，如塘沽的振幅和遲角偏差較大，可能與近岸水深地形刻畫不夠細(xì)致有關(guān)．

圖5 為M2和K1分潮的同潮圖，模擬結(jié)果表明M2分潮在渤海海域分布2 個無潮點，分別位于黃河口附近(119°03′E，38°11′N)和秦皇島外海(120°06′E，39°48′N)，K1分潮無潮點位于渤海海峽(120°42′E，38°16′N)，與沈育疆[18]、Fang[19]和朱學(xué)明等[20]研究中的無潮點位置接近．從等振幅線和等遲角線分布來看，本文模擬結(jié)果與已有研究成果都基本一致，能較好展現(xiàn)渤海潮波傳播與分布特征．

圖5 M2、K1 分潮同潮圖Fig.5 Cotidal chart of M2 and K1 constituents in the Bohai Sea

2.2 渤海灣局部區(qū)域潮流過程驗證

2.2.1 模擬區(qū)域及計算配置

為模擬渤海灣東北部局部區(qū)域的潮流運動，對計算區(qū)域加密，整個計算區(qū)域共有19 436 個四邊形單元和18 809 個網(wǎng)格節(jié)點，網(wǎng)格最大和最小分辨率分別約為0.05°和0.000 6°，相當(dāng)于5 000 m 和60 m，如圖6 所示．圖7 顯示擬進(jìn)行模擬和實測水位、流速和流向的渤海驗潮站分布情況，其中T1 為潮位測站，V1～V6 為潮流測站．區(qū)域地形數(shù)據(jù)取自1∶50 000海圖(2011 年)，測站附近選用實測地形資料，模擬區(qū)域底摩阻系數(shù)統(tǒng)一取0.002，模型計算時間步長為0.8 s，模擬時間段為2012 年10 月8 日至2012 年10月20 日．

圖6 渤海灣加密計算域水平網(wǎng)格投影分布Fig.6 Distribution of refined meshes after being projected in the Bohai Bay

圖7 潮流驗證點空間分布圖Fig.7 Observation stations for model verification in the Bohai Bay

2.2.2 潮位和流速、流向驗證

圖8 顯示了2012 年10 月10 日至2012 年10 月20 日實測潮位過程和模擬潮位過程的比較情況．可以看出，T1 潮位過程驗證點的模擬結(jié)果與實測水位過程比較吻合，相位基本一致，高潮和低潮位模擬都較好，模擬潮位最大誤差為0.17 m．

圖9 顯示了2012 年10 月18 日至19 日大潮期間實測和模擬的流速、流向的比較情況．總體來看，各測點模擬流速、流向與實測結(jié)果吻合較好，部分時刻兩者的差距可能和地形刻畫不夠完全準(zhǔn)確以及底摩阻系數(shù)簡單取為常數(shù)有關(guān)．

圖8 T1驗潮站潮位驗證Fig.8 Tide level verification in station T1

圖9 V1～V6測站潮流驗證Fig.9 Tidal current verification in stations V1—V6

總結(jié)上面兩種不同區(qū)域尺度的模擬與實測結(jié)果的比較可以得知，采用極射赤面投影將球面坐標(biāo)轉(zhuǎn)換為平面坐標(biāo)的方法建立的二維無積分間斷有限元數(shù)值模型，可以較好地模擬大范圍海域和局部海域的潮波傳播和潮流運動．

3 結(jié)論與展望

為了利用間斷有限元方法實現(xiàn)大范圍海域的二維水動力模擬，采用極射赤面投影將球面坐標(biāo)轉(zhuǎn)換為平面坐標(biāo)的方法求解球坐標(biāo)下的淺水流動控制方程，建立了基于無積分節(jié)點間斷有限元方法的二維水動力數(shù)值模型，并將建立的模型運用到整個渤海海域潮波和渤海灣局部潮流模擬．模型模擬結(jié)果與實測結(jié)果的比較表明，所建立的模型不僅可以較好地模擬大范圍海域的潮波傳播，而且通過非結(jié)構(gòu)化四邊形網(wǎng)格的局部加密，也可以較好地模擬復(fù)雜地形附近水位、流速和流向變化規(guī)律．今后我們將進(jìn)一步加入波流耦合等功能，將該模型進(jìn)一步推廣應(yīng)用于大范圍海域風(fēng)暴潮等問題的模擬．