盧 付 強(qiáng),宋 志 堯,雷 夢 玲,丁 凱 孟

(1.南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023;2.南京師范大學(xué)大規(guī)模復(fù)雜系統(tǒng)數(shù)值模擬江蘇省重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023;3.江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心，江蘇 南京 210023)

基于Laplace方程的球面正交曲線格網(wǎng)生成方法

盧 付 強(qiáng)1,2,3,宋 志 堯1,2,3,雷 夢 玲1,3,丁 凱 孟1,3

(1.南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023;2.南京師范大學(xué)大規(guī)模復(fù)雜系統(tǒng)數(shù)值模擬江蘇省重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023;3.江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心，江蘇 南京 210023)

為克服基于傳統(tǒng)經(jīng)緯度格網(wǎng)的全球海洋及大氣數(shù)值模式產(chǎn)生的極點(diǎn)問題，同時顧及全球格網(wǎng)編碼便捷及信息查詢,該文將立方體球心投影的格網(wǎng)經(jīng)球極投影轉(zhuǎn)換到二維平面區(qū)域，再應(yīng)用Laplace方程數(shù)值求解生成正交曲線格網(wǎng)；然后采用球極逆投影,得到球面上一種新的正交曲線格網(wǎng)。通過與等角球心投影、保角變換兩種球面格網(wǎng)生成方法的結(jié)果相比較,從格網(wǎng)兩邊的長度、單元的面積、正交性等指標(biāo)進(jìn)行統(tǒng)計(jì)分析，對所生成格網(wǎng)的質(zhì)量進(jìn)行評價。結(jié)果表明,等角球心投影法生成的格網(wǎng)單元大小較均勻，但正交性不理想；Laplace方程法生成的正交曲線格網(wǎng)稍優(yōu)于保角變換法，可作為應(yīng)用有限差分方法和有限體積方法的全球海洋及大氣大尺度數(shù)值模擬的基礎(chǔ)格網(wǎng)。

Laplace方程；球極投影；正交曲線格網(wǎng)；Neumann正交化

0 引言

隨著全球氣候逐漸變暖，極端天氣現(xiàn)象發(fā)生的頻率呈升高趨勢,需對全球氣候變化進(jìn)行及時預(yù)測，以便采取可能的措施減小惡劣天氣產(chǎn)生的危害。而這些預(yù)測在很大程度上依賴于全球大氣模式、海洋模式的發(fā)展，因此需要構(gòu)建一個理想的全球大氣、海洋數(shù)值模型都適用的網(wǎng)格，以便進(jìn)一步實(shí)現(xiàn)相互耦合的數(shù)值模擬。建立全球大氣、海洋數(shù)值模型的第一步就需要選擇一個好的全球離散網(wǎng)格。目前在大氣、海洋數(shù)值模式中所采用的網(wǎng)格分為兩大類[1]：一類是有結(jié)構(gòu)網(wǎng)格，網(wǎng)格的排列結(jié)構(gòu)縱橫有序，相鄰節(jié)點(diǎn)或網(wǎng)格單元的數(shù)目是固定的，易于結(jié)點(diǎn)(單元)的編號和布置變量，便于方程的離散和計(jì)算編程，如傳統(tǒng)或退化的經(jīng)緯度格網(wǎng)[2,3]、立方體球心投影網(wǎng)格[4,5]、陰陽交疊網(wǎng)格[6,7]等。另一類是無結(jié)構(gòu)網(wǎng)格，即由任意三角形或六邊形組成的網(wǎng)格，如基于正二十面體產(chǎn)生的球面三角網(wǎng)格[8,9]、球面正六邊形網(wǎng)格[10,11]。此類網(wǎng)格在全球大氣或海洋運(yùn)動數(shù)值模型中得到了廣泛應(yīng)用，其主要優(yōu)點(diǎn)是：整體上網(wǎng)格大小較均勻。但在這類網(wǎng)格上對數(shù)據(jù)進(jìn)行操作，需要預(yù)先對每個網(wǎng)格建立索引、拓?fù)潢P(guān)系表，并且在這類網(wǎng)格上建立數(shù)值模型的求解方法通常較復(fù)雜。

正交曲線網(wǎng)格屬于有結(jié)構(gòu)網(wǎng)格，具有很多優(yōu)點(diǎn)：首先，網(wǎng)格的編碼與領(lǐng)域查詢比較方便,模型方程組在正交曲線坐標(biāo)系中的表達(dá)形式比較簡潔；其次，便于數(shù)值有限差分方法或者有限體積方法的設(shè)計(jì)和實(shí)現(xiàn)等。許多大氣數(shù)值模式建立在經(jīng)緯度網(wǎng)格上，但是經(jīng)緯度網(wǎng)格在趨于兩個極點(diǎn)時經(jīng)線收斂，除了譜變換方法，在采用基于格點(diǎn)的方法進(jìn)行數(shù)值計(jì)算時會引發(fā)極點(diǎn)問題[12]，同時也產(chǎn)生了解決這些問題的方法：高緯度濾波、半隱式半拉格朗日方法等。但是這些方法在全球海洋模式中卻難以應(yīng)用[13]，所以必須考慮采用一種新的正交網(wǎng)格，采用數(shù)值求解Laplace方程組，結(jié)合球極投影方法，可生成球面正交曲線網(wǎng)格，既利于全球地理信息的表達(dá)，也適用于地球流體大尺度運(yùn)動的數(shù)值模擬，同時也是生成三維球體準(zhǔn)正交的曲面六面體網(wǎng)格的基礎(chǔ)。

1 球面四邊形網(wǎng)格生成方法

1.1 立方體球心投影方法

考慮球的內(nèi)接立方體，分別在立方體的每個表面上建立局部二維笛卡爾坐標(biāo)，然后做一個伸縮變換，將表面上每個點(diǎn)變換到外接球面上。假設(shè)球半徑為R，可知立方體的邊長2a與半徑有關(guān)系式：

(1)

選取立方體中心為三維笛卡爾坐標(biāo)系的原點(diǎn)，坐標(biāo)軸X、Y、Z分別與立方體的前面、右面、頂面垂直(圖1)，其中左子圖為立方體各面的展開圖，右子圖為立方體與外接球面在三維笛卡爾坐標(biāo)系中的位置關(guān)系。記(x,y)為立方體每個表面上的局部二維笛卡爾坐標(biāo)，有-a≤x,y≤a。

圖1 立方體6個表面的展開圖及其外接球面在三維笛卡爾坐標(biāo)系中的位置

Fig.1 Layout of an open cube and the cube and its circum-sphere under 3D absolute Cartesian coordinates

從立方體表面變換到球面上有兩種方法：

(1)立方體球心等距投影網(wǎng)格。直接對局部坐標(biāo)進(jìn)行等距網(wǎng)格剖分(圖2)，其中顯示了在立方體表面X=a上的局部二維平面笛卡爾坐標(biāo)(x,y)與半徑為R的球面上點(diǎn)(X,Y,Z)之間的關(guān)系：

(2)

(2)立方體球心等角投影網(wǎng)格。首先對局部笛卡爾坐標(biāo)做一個變換：

(3)

圖2 立方體表面X=a上的點(diǎn)經(jīng)過球心投影到球面上

Fig.2 Gnomonic projection of point on the cube surfaceX=ato the sphere surface

1.2 保角變換方法

應(yīng)用保角變換將網(wǎng)格點(diǎn)在各復(fù)平面與球面之間變換時，網(wǎng)格的角度保持不變，并且除去8個角點(diǎn)，映射在每個網(wǎng)格點(diǎn)上是光滑的，保角變換方法的主要思想是先采用廣義球極投影將球面上的點(diǎn)映射到復(fù)平面上，然后在復(fù)平面之間采用Taylor 級數(shù)展開式實(shí)現(xiàn)一系列復(fù)變換及逆變換[14]。保角變換生成的網(wǎng)格總體上正交性很好并且網(wǎng)格的長寬比接近1，即形狀與正方形相近。但在8個奇點(diǎn)(角點(diǎn))處，網(wǎng)格偏離正交性固定在30°，在角點(diǎn)附近的網(wǎng)格正交性也受到影響，正交性有所降低。盡管不如球面經(jīng)緯網(wǎng)格在極點(diǎn)處收斂那么嚴(yán)重,嚴(yán)格的正交性要求將導(dǎo)致網(wǎng)格在8個奇點(diǎn)處收斂,可以對網(wǎng)格在局部進(jìn)行拉伸以改善這種狀況。

2 球面正交曲線網(wǎng)格生成

采用某種投影方法將球面上的網(wǎng)格變換到二維平面區(qū)域，然后借鑒平面上數(shù)值網(wǎng)格生成方法生成正交曲線網(wǎng)格，最后采用逆變換，將生成的網(wǎng)格變換到球面上。這個過程的主要思路可由圖3描述。在一般情形下，采用投影變換后，網(wǎng)格的幾何形狀、面積、角度等屬性都會發(fā)生變化。在此選取常用的3種地圖投影變換方法：圓錐、圓柱及赤平投影[18,19]，這些變換是保持正交性質(zhì)的，即網(wǎng)格在球面與平面之間變換時，格網(wǎng)的角度不發(fā)生變化。

圖3 算法的總體路線

Fig.3 The process of the presented algorithm

2.1 地圖投影變換

由于立方體6個表面上的網(wǎng)格是對稱的，所以只需在一個面上生成網(wǎng)格，然后采取合適的旋轉(zhuǎn)變換就可得到覆蓋整個球面的網(wǎng)格。在此選取赤道區(qū)域所在的一個球面區(qū)域作為基準(zhǔn)面，圖4顯示了立方體球面前表面上初始網(wǎng)格分別經(jīng)過Lambert正形圓錐投影、Mercator投影、赤平投影變換后的網(wǎng)格。立方體球面前表面上的格網(wǎng)關(guān)于格網(wǎng)中心具有四重對稱性，從圖4可看出：經(jīng)過圓錐或圓柱投影，格網(wǎng)的四重對稱性有損失；而經(jīng)過赤平投影變換，格網(wǎng)的對稱性得到保持，所以采用赤平投影或球極投影。

圖4 經(jīng)過圓柱、圓錐以及赤平投影變換之后,立方體球面前表面上格網(wǎng)點(diǎn)的分布

Fig.4 Distribution of gird points on front face of cubed sphere under Mercator projection,Lambert conformal conic projection and stereographic projection

2.2 Laplace微分方程方法

Laplace微分方程是一種橢圓形偏微分方程，橢圓形網(wǎng)格生成是通過求解擬線性橢圓形偏微分方程組的一種網(wǎng)格生成方法，可以對網(wǎng)格的光滑性及正交性進(jìn)行控制。在生成球面網(wǎng)格時，初始網(wǎng)格取立方體等距球心投影網(wǎng)格，將在球面上物理網(wǎng)格的每個表面通過球極投影變換到二維的平面空間，然后對網(wǎng)格施加所需的限制條件，在投影平面上求解橢圓型偏微分方程組。在大部分基于有限差分方法的地球流體動力學(xué)數(shù)值模擬中，網(wǎng)格的正交性及網(wǎng)格尺寸的均勻性是最重要的，對網(wǎng)格施加正交性有Neumann條件限制和Dirichlet條件限制[20]兩種方法。Neumann條件限制：網(wǎng)格點(diǎn)可以沿著每個面的邊界滑動，而Dirichlet條件限制則是固定在邊界上的點(diǎn)，通過調(diào)整網(wǎng)格內(nèi)點(diǎn)并且保持網(wǎng)格尺寸的準(zhǔn)一致性。Neumann條件限制通過要求未知量對兩個參數(shù)的偏導(dǎo)數(shù)的乘積等于零來計(jì)算每個網(wǎng)格點(diǎn)沿網(wǎng)格邊界的偏移量，但這種做法將導(dǎo)致網(wǎng)格單元在角點(diǎn)鄰近區(qū)域趨于收斂，網(wǎng)格單元的面積偏小，可在局部對網(wǎng)格進(jìn)行拉伸來改善這種狀況。

假定一個二維初始變換W:(ξ,η)→(x,y),將計(jì)算空間[1,M]×[1,N]變換到物理空間Ω?R2,其中正整數(shù)M、N表示沿兩個方向的網(wǎng)格點(diǎn)數(shù)目，網(wǎng)格線索引為ξ=1,2,…,M;η=1,2,…,N。初始變換常由一般的方法如代數(shù)方法中的超限插值技術(shù)得到。對于給定的初始變換，二維平面正交曲線網(wǎng)格的生成方法多采用數(shù)值求解下列擬線性偏微分方程組[21]:

g22(xξξ+Pxξ)-2g12xξη+g11(xηη+Qxη)=0
g22(yξξ+Pyξ)-2g12yξη+g11(yηη+Qyη)=0

(4)

其中：

(5)

(6)

(7)

其中:(x,y)為原二維物理區(qū)域中的網(wǎng)格點(diǎn)坐標(biāo)，(ξ,η)為變換區(qū)域中網(wǎng)格點(diǎn)坐標(biāo)。網(wǎng)格的度量項(xiàng)g11、g22分別表示網(wǎng)格沿兩個方向的長度，g12可用于表示網(wǎng)格單元的角度大小。P、Q表示網(wǎng)格點(diǎn)疏密的控制函數(shù)，可通過改變P、Q對網(wǎng)格點(diǎn)的分布進(jìn)行調(diào)整，若P=Q=0，則產(chǎn)生步長均勻的網(wǎng)格。

為了求解擬線性橢圓形方程組(4)，需給出必要的邊界條件，此處設(shè)置控制函數(shù)P=Q=0，取Neumann正交邊界條件，允許邊界點(diǎn)通過自由滑動達(dá)到邊界上的正交。采用有限差分方法求解方程組(4)，其中網(wǎng)格的度量項(xiàng)g11、g12、g22可以根據(jù)網(wǎng)格點(diǎn)坐標(biāo)表達(dá)式預(yù)先計(jì)算出來，或者采用差分近似方法，系數(shù)中含有待求解函數(shù)的非線性項(xiàng)，采用前一次的迭代值計(jì)算這些系數(shù)，只對二階偏導(dǎo)數(shù)以及與控制函數(shù)有關(guān)的導(dǎo)數(shù)項(xiàng)進(jìn)行離散，采用中心差分近似，選取映射到計(jì)算空間中的離散步長Δξ=Δη=1，則:

(8)

(9)

采用超松弛迭代方法(SOR)求解上述方程組：

(10)

其中wi,j為松弛因子，應(yīng)滿足條件：0

2.3 Neumann正交邊界條件

在求解區(qū)域邊界上應(yīng)該滿足正交性條件：

XξXη=0,在ξ=1,M

(11)

XξXη=0在η=1,N

(12)

其中X表示網(wǎng)格點(diǎn)的二維坐標(biāo)，采用有限差分方法對式(11)、式(12)進(jìn)行近似求解，在迭代求解方程組時，每次循環(huán)中邊界條件式都要進(jìn)行更新。為了防止邊界點(diǎn)滑動距離過大，超過其鄰近點(diǎn)，需要設(shè)置臨界值，特別是在鈍角區(qū)域附近容易出現(xiàn)這種情形。合適的臨界值對邊界每個點(diǎn)的滑動距離起限制作用，一般要求在每次迭代過程中距離邊界上最近兩個鄰近點(diǎn)的距離至多減少一半。

3 結(jié)果分析與比較

3.1 網(wǎng)格質(zhì)量評價指標(biāo)

網(wǎng)格的4個指標(biāo)有：網(wǎng)格分別沿參數(shù)ξ方向與η方向的長度hξ,hη網(wǎng)格單元的面積大小area以及角度的正交性偏差DO：

(13)

(14)

3.2 計(jì)算結(jié)果分析

選取立方體的頂面Z=a作為基礎(chǔ)，在其上建立局部二維笛卡爾坐標(biāo)系(x,y)，然后變換到半徑為R的球面上。取a=1，采用北極投影，將立方體球心等距投影網(wǎng)格映射到與北極點(diǎn)相切的二維平面上，數(shù)值求解橢圓形微分方程組(4)，采用Neumann正交邊界條件(式(11)、式(12))，控制函數(shù)取值P=Q=0；然后再采用北極逆投影得到球面上的正交網(wǎng)格。圖5給出了在球極投影下，各網(wǎng)格點(diǎn)處地圖放大系數(shù)的等值線圖，可以看出投影系數(shù)接近1 ，即網(wǎng)格經(jīng)過球極投影后，幾何度量變化不是很大。

目前，也可采用等角或等距球心投影方法以及保角映射方法生成立方體球面網(wǎng)格，但等角或等距球心投影產(chǎn)生的網(wǎng)格不具有正交性，而保角變換方法或求解橢圓形方程產(chǎn)生的網(wǎng)格正交性較好。下面給出上述方法生成網(wǎng)格的度量性質(zhì)的比較結(jié)果。

選取網(wǎng)格單元數(shù)分別為45×45、90×90、180×180的情形時進(jìn)行計(jì)算，相對應(yīng)的空間分辨率大約為2°、1°、30°。表1給出了不同網(wǎng)格數(shù)時，分別采用等角投影方法、保角變換、Laplace方程方法生成網(wǎng)格的統(tǒng)計(jì)結(jié)果，包括網(wǎng)格單元最大邊與最短邊的長度、單元角度的正交性最大偏差以及平均偏差、網(wǎng)格單元面積的最大值和最小值。

圖5 各網(wǎng)格點(diǎn)處的地圖放大因子等值線

Fig.5 Map factor of each grid point after elliptic-smoothing under north polar stereographic projection

如果頂面的網(wǎng)格已經(jīng)生成，則采用合適的旋轉(zhuǎn)變換生成其他5個表面上的正交網(wǎng)格。圖6顯示了網(wǎng)格數(shù)在45×45×6時，采用等角投影方法、保角變換、 Laplace方程方法生成球面網(wǎng)格的正交性偏差及網(wǎng)格單元面積的變化?？煽闯霰＝亲儞Q、Laplace方程方法生成的球面網(wǎng)格非常相似，絕大多數(shù)網(wǎng)格的正交性偏差不超過10°，只是接近8個角點(diǎn)的網(wǎng)格單元的正交性偏差增大了一些，最大偏差發(fā)生在8個角點(diǎn)處。另外，從表1也可看出，保角變換與Laplace方程方法生成網(wǎng)格的統(tǒng)計(jì)量差別相當(dāng)小，基本在同一個數(shù)量級。綜合起來，在相同情形下，等角投影方法生成球面網(wǎng)格單元面積的大小較均勻，但正交性偏差較大；而保角變換、Laplace方程方法生成球面網(wǎng)格單元的面積總體上不如等角投影方法的均勻，尤其是在8個角點(diǎn)附近區(qū)域網(wǎng)格單元的面積偏小。

圖 7 給出了網(wǎng)格逐次分半加密時，分別采用三種方法，F(xiàn)ortran編程的CPU消耗時間隨網(wǎng)格數(shù)目的變化，可以看出在同等格網(wǎng)數(shù)目的條件下，求解Laplace方程組耗時最多，等角球心投影方法耗時最少，保角變換方法耗時稍多于等角球心投影方法。即便如此，生成空間分辨率為幾分的格網(wǎng)最多時間耗費(fèi)不過10 s，并且在發(fā)展全球大氣或海洋數(shù)值模式過程中，網(wǎng)格生成是一次性的，這個時間是可接受的。

4 結(jié)語

基于立方體球心投影網(wǎng)格，結(jié)合球極投影技術(shù)，采用Laplace方程組在球面上生成了一種廣義正交曲線網(wǎng)格，并且對網(wǎng)格的質(zhì)量進(jìn)行分析評價。通過與等角球心投影、保角變換方法的結(jié)果進(jìn)行比較，表明等角球心投影生成的網(wǎng)格單元大小較均勻，但網(wǎng)格不具有正交性；而保角變換方法、數(shù)值求解微分方程組生成的網(wǎng)格的正交性質(zhì)較好，但在8個角點(diǎn)鄰域的網(wǎng)格單元面積偏小。在應(yīng)用于實(shí)際大氣或海洋大尺度運(yùn)動的具體數(shù)值模型中，這個缺陷可以通過對網(wǎng)格進(jìn)行拉伸，或采用隱式時間積分方法來彌補(bǔ)。因此,這種廣義正交曲線網(wǎng)格可用于地理空間數(shù)據(jù)的表達(dá)以及全球大氣或海洋運(yùn)動的數(shù)值模擬。

表1 在不同網(wǎng)格數(shù)情形下三種方法生成網(wǎng)格的統(tǒng)計(jì)結(jié)果

Table 1 Quality statistics of grids generated by three different methods with different cell-numbers

45×4590×90180×180等角投影方法保角變換方法Laplace方程方法網(wǎng)格邊長單元面積正交性偏差(°)網(wǎng)格邊長單元面積正交性偏差(°)網(wǎng)格邊長單元面積正交性偏差(°)最大值3.4906E-021.7453E-028.7266E-03最小值2.4684E-021.2342E-026.1707E-03均值3.2367E-021.6212E-028.1132E-03方差6.9822E-061.6732E-064.0929E-07最大值1.2000E-033.0648E-047.7482E-05最小值8.7658E-042.1682E-045.3487E-05均值1.0343E-032.5857E-046.4642E-05方差8.5667E-095.364E-0103.3550E-011最大值28.814229.421529.7109均值8.06318.06568.0661最大值5.1088E-022.5973E-021.3097E-02最小值1.7561E-027.1262E-032.8604E-03均值3.2365E-021.6188E-028.1008E-03方差3.2343E-058.1043E-062.0272E-06最大值1.3519E-033.4438E-048.7732E-05最小值4.2396E-046.9325E-051.1316E-05均值1.0343E-032.5856E-046.4643E-05方差2.2336E-081.3998E-098.759E-011最大值12.922512.930312.9310均值0.45220.22610.1131最大值3.6401E-021.8622E-021.0004E-02最小值6.5321E-032.5311E-038.9817E-04均值3.1804E-021.5916E-027.9523E-03方差1.7828E-054.7045E-061.5216E-06最大值1.3000E-033.4978E-041.0225E-04最小值7.5305E-051.1096E-051.8763E-06均值1.0343E-032.5857E-046.4643E-05方差5.2956E-083.7047E-093.137E-010最大值12.688912.564712.5055均值0.73910.35000.2063

圖6 不同投影生成網(wǎng)格的角度正交性偏差與網(wǎng)格單元面積大小 圖7 三種方法的計(jì)算效率隨網(wǎng)格數(shù)目的變化

Fig.6 Angle deviation orthogonality and cell area on equi-angular sphere grids based on different method Fig.7 The efficiency of three methods under different grid-numbers

[1] STANIFORTH A,THUBURN J.Horizontal grids for global weather and climate prediction models:A review[J].Quarterly Journal of the Royal Meteorological Society,2012,138:1-26.

[2] HOLLOWAY L,SPELMAN M J,MANABE J.Latitude-longitude grid suitable for numerical time integration of a global atmospheric model[J].Monthly Weather Review,1973,101(1):69-78.

[3] FADEEV R Y.Algorithm for reduced grid generation on a sphere for a global finite difference atmospheric model[J].Computational Mathematics and Mathematical Physics,2013,53(2):237-252.

[4] RONCHI C,IACONO R,PAOLUCCI P S.The "Cubed Sphere":A new method for the solution of partial differential equations in spherical geometry[J].Journal of Computational Physics,1996,124(1):93-114.

[5] MCGREGOR J L.Semi-Lagrangian advection on a cubic gnomonic projection of the sphere[J].Atmosphere-Ocean,1997,35(1):153-169.

[6] KAGEYAMA A,SATO T."Yin-Yang grid":An overset grid in spherical geometry[J].Geochem.Geophys.Geosyst.,2004,5,Q09005,doi:10.1029/2004GC000734.

[7] 艾細(xì)根,劉宇迪,崔新東,等.基于球面陰陽網(wǎng)格數(shù)值計(jì)算研究進(jìn)展[J].氣象科技,2014,42(1):13-22.

[8] TOMITA H,SATOH M,GOTO K.An optimization of the Icosahedral Grid modified by Spring Dynamics[J].Journal of Computational Physics,2002,183(1):307-331.

[9] IGA S,TOMITA H.Improved smoothness and homogeneity of icosahedral grids using the spring dynamics method[J].Journal of Computational Physics,2014,258(1):208-226.

[10] MIURA H,KIMOTO M.A comparison of grid quality of optimized spherical hexagonal-pentagonal geodesic grids[J].Monthly Weather Review,2005,133(10):2817-2833.

[11] 童曉沖,奔進(jìn),汪瀅.利用數(shù)值投影變換構(gòu)建全球六邊形離散格網(wǎng)[J].測繪學(xué)報,2013,42(2):268-276.

[12] NAUGHTON M,COURTIER P. A pole problem in the reduced Gaussian grid[J].Quarterly Journal of the Royal Meteorological Society,1994,120:1389-1407.

[13] ADCROFT A,CAMPIN J M,HILL C,et al.Implementation of an atmosphere-ocean general circulation model on the expanded spherical cube[J].Monthly Weather Review,2004,132:2845-2863.

[15] NAIR R D,THOMAS S J,LOFT R D.A discontinuous Galerkin transport scheme on the cubed sphere[J].Monthly Weather Review,2004,133(4):814-828.

[16] LU F,SONG Z,FANG X.Generation of Orthogonal Curvilinear Grids on the Sphere Surface Based on Laplace-Beltrami Equations[C].IOP Conference Series:Earth and Environmental Science,2014,(19) 012012.

[17] PUTMAN W M,LIN S J.Finite-volume transport on various cubed-sphere grids[J].Journal of Computational Physics,2007,227(1):55-78.

[18] 沈桐立,田永祥,葛孝貞.數(shù)值天氣預(yù)報[M].北京:氣象出版社,2003.60-71.

[19] 王美玲,付夢印.地圖投影與坐標(biāo)變換[M].北京:電子工業(yè)出版社,2014.51-125.

[20] THOMPSON J F,SONI B K,WEATHERILL N P.Handbook of Grid Generation[M].CRC Press,1999,ISBN 0-8493-2687-7.

[21] SPEKREIJSE S P. Elliptic grid generation based on Laplace equations and algebraic transformations[J].Journal of Computational Physics,1995,118(1):38-61.

Generation Method of Orthogonal Curvilinear Grids on Sphere Surface Based on Laplace Equations

LU Fu-qiang1,2,3,SONG Zhi-yao1,2,3,LEI Meng-ling1,3,DING Kai-meng1,3

(1.Key Laboratory of Virtual Geographic Environment(Ministry of Education),Nanjing Normal University，Nanjing 210023;2.Jiangsu Key Laboratory for Numerical Simulation of Large Scale Complex Systems,Nanjing Normal University,Nanjing 210023;3.Jiangsu Center for Collaborative Innovation in Geographical Information Resource Development and Application,Nanjing 210023,China)

In order to overcome the pole-problem induced by grid-based methods on traditional longitude-latitude grid for global atmosphere and ocean circulation models,a new kind of orthogonal curvilinear grids based on the gnomonic cube grids on the sphere was presented.Firstly,the gnomonic cube grids was on a two-dimensional plane by stereographic projection,then Laplace equations are solved to generate quasi-orthogonal curvilinear grids on the projection plane.Finally,the grids are transformed back to the sphere surface by the inverse of stereographic projection.The produced grid configuration has quasi-uniform cell-size on the whole sphere surface except the neighbourhood of eight major corner vertices.Moreover,another two methods which include equi-angular gnomonic projection and conformal mapping were also introduced to generate spherical grids.The resulting comparisons between three methods would be presented at last.Some quantities such as the grid length along two directions,the angle deviation from orthogonality and the area of the cell were used to evaluate the quality of the grids generated by three methods.The results show that the equi-angular gnomonic projection produces grids of the most uniform size,but are non-orthogonal,while elliptic-smoothing combined with stereographic projection method produces quasi-orthogonal curvilinear grids that are much like that of conformal mapping method.The small grid cells in the neighbourhood of eight corner vertices could be improved by stretching method,or this small defect could be circumvented by implicit time integration method,which demonstrate the kind of grid newly produced is fit to be a model grid on which the finite difference method or finite volume method can be implemented for numerical simulation of global atmosphere or ocean dynamics on large scale.

Laplace equations;stereographic projection;orthogonal curvilinear grids;Neumann orthogonal boundary condition

2015-06-18；

2015-10-13

國家自然科學(xué)基金資助項(xiàng)目(41306078)

盧付強(qiáng)(1984-),男,博士,研究方向?yàn)榈貓D學(xué)與地理信息系統(tǒng)。E-mail:lufuqiang2000@163.com

10.3969/j.issn.1672-0504.2016.01.004

P208

A

1672-0504(2016)01-0018-06