基于DISTMESH的三角形網(wǎng)格自適應(yīng)生成方法與應(yīng)用

張靖婷王憲業(yè)李 鋮葛建忠*

(1.華東師范大學(xué) 河口海岸國家重點(diǎn)實(shí)驗(yàn)室,上海201100；2.上海市海洋環(huán)境預(yù)報中心,上海200062)

海洋數(shù)值計(jì)算的基礎(chǔ)與核心是在時間和空間上對流體動力學(xué)方程進(jìn)行離散。在計(jì)算區(qū)域中,水平方向空間離散而成的平面網(wǎng)格,主要分為結(jié)構(gòu)、無結(jié)構(gòu)和混合網(wǎng)格[1]。當(dāng)前主流的海洋數(shù)值模式中,DELFT3D[2],ROMS[3],POM[4]等模式采用了結(jié)構(gòu)網(wǎng)格,常用正交或非正交曲線網(wǎng)格,TELEMAC[5],FVCOM[6]等模式采用了無結(jié)構(gòu)網(wǎng)格,網(wǎng)格單元為三角形,DELFT-FM[7]模式采用了混合網(wǎng)格。結(jié)構(gòu)網(wǎng)格計(jì)算高效,在簡單區(qū)域能充分發(fā)揮優(yōu)勢,但其對復(fù)雜區(qū)域的幾何擬合度較差。無結(jié)構(gòu)網(wǎng)格能較好地擬合不規(guī)則邊界,方便局部加密,且能迅速改變分辨率,適用范圍更廣,尤其針對復(fù)雜區(qū)域,無結(jié)構(gòu)網(wǎng)格更具優(yōu)勢。

海洋數(shù)值計(jì)算時,經(jīng)常需要離散邊界不規(guī)則、地形多變的復(fù)雜區(qū)域。這些區(qū)域的網(wǎng)格生成,不僅應(yīng)保證各個單元的質(zhì)量,還應(yīng)合理控制分辨率,在岸界、地形顯著變化處生成小單元,保證離散的幾何擬合度,在其余部分生成大單元,避免網(wǎng)格規(guī)模的不必要擴(kuò)大,從而提高計(jì)算效率；同時,網(wǎng)格應(yīng)使局部曲線光滑和整體疏密過渡自然。在保證上述質(zhì)量要求的前提下,應(yīng)盡可能縮短網(wǎng)格生成所用的時間,從而提高整體效率。

生成三角形網(wǎng)格的方法較多,采用商業(yè)軟件(如SMS[8])、開源軟件(如TRIANGLE[9]和EASYMESH[10])等均可快速生成三角形網(wǎng)格,TRIANGLE,GMESH[11],DISTMESH[12]等三角形網(wǎng)格生成方法的一項(xiàng)對比測試[13]表明,針對區(qū)域的均勻剖分,隨著區(qū)域復(fù)雜性的提高和網(wǎng)格規(guī)模的擴(kuò)大,TRIANGLE生成網(wǎng)格速度最快,DISTMESH質(zhì)量表現(xiàn)最優(yōu)。但隨著區(qū)域復(fù)雜性的提高和網(wǎng)格規(guī)模的擴(kuò)大,各種已有軟件及開源庫所生成的網(wǎng)格質(zhì)量問題凸顯,如出現(xiàn)鈍角三角形、相鄰單元面積比過大等。這些問題的解決只能依賴人工手動優(yōu)化,并且,其局部加密也需要較多的人工控制。利用人工網(wǎng)格調(diào)整手段提高網(wǎng)格質(zhì)量相當(dāng)耗時,需進(jìn)行大量的調(diào)試和優(yōu)化,難以提高網(wǎng)格生成效率,開展網(wǎng)格自適應(yīng)生成研究可有效緩解網(wǎng)格質(zhì)量和生成效率的矛盾。

針對海洋數(shù)值計(jì)算中常見的邊界不規(guī)則、地形多變的復(fù)雜區(qū)域,目前的網(wǎng)格自適應(yīng)生成方法缺乏實(shí)用性,難以在保證質(zhì)量和效率的同時,滿足實(shí)際問題的復(fù)雜分辨率要求?；贒ISTMESH方法,本研究提出一種普適性的、能精確擬合復(fù)雜岸線和地形的三角形網(wǎng)格自適應(yīng)生成方法,以進(jìn)行海洋動力學(xué)數(shù)值模擬空間離散,并將該方法應(yīng)用于理想潮汐汊道和崇明東灘復(fù)雜潮溝系統(tǒng)的網(wǎng)格生成,其中崇明東灘的網(wǎng)格已被用于FVCOM模型流場計(jì)算。

1 DISTMESH方法

DISTMESH方法適用于單連通和復(fù)連通的離散區(qū)域,前者僅包含外邊界(指岸線與開邊界),后者包含外邊界和內(nèi)邊界(如島嶼、水工建筑物邊界)。DISTMESH方法通過計(jì)算機(jī)迭代,一定時間內(nèi)使網(wǎng)格滿足質(zhì)量和分辨率要求。其網(wǎng)格自適應(yīng)生成過程如圖1所示。每一次迭代過程,包括3個步驟：1)對給定點(diǎn)集,利用Delaunay[14]三角剖分算法進(jìn)行網(wǎng)格剖分,確定拓?fù)浣Y(jié)構(gòu)；2)利用距離函數(shù)確定區(qū)域內(nèi)的點(diǎn),即保留上一次迭代后留在區(qū)域內(nèi)的點(diǎn)集；3)確定保留點(diǎn)集的位置及其拓?fù)浣Y(jié)構(gòu)后,由分辨率函數(shù)確定“力-位移”關(guān)系,從而將點(diǎn)移動至新位置,也即給定下一次迭代過程的初始點(diǎn)集。其中,距離函數(shù)帶正、負(fù)號,分別與區(qū)域外、內(nèi)相對應(yīng),通過正負(fù)判斷確定保留點(diǎn)集,距離函數(shù)的絕對值為平面上各點(diǎn)到邊界的最近距離；Delaunay三角剖分算法具有空外接圓和最大化最小角特性,因此,對給定點(diǎn)集,能給出質(zhì)量最優(yōu)的拓?fù)浣Y(jié)構(gòu)；分辨率函數(shù)給出模型計(jì)算所需網(wǎng)格中各點(diǎn)的分辨率,提出了對目標(biāo)網(wǎng)格分辨率即單元的大小和分布的要求；“力-位移”關(guān)系由分辨率函數(shù)確定,由網(wǎng)格現(xiàn)有分辨率與目標(biāo)分辨率之間存在的差異構(gòu)造出對三角形邊的拉“力”,“力”使邊長“變形”,即點(diǎn)發(fā)生“位移”。在不斷迭代中,“力-位移”關(guān)系的建立與求解最終使網(wǎng)格分辨率符合要求。

圖1 DISTMESH方法網(wǎng)格自適應(yīng)生成過程Fig.1 The adaptive mesh generating process of DISTMESH

1.1 “力-位移”關(guān)系的建立與求解

給定網(wǎng)格點(diǎn)的當(dāng)前位置和拓?fù)浣Y(jié)構(gòu)后,需求解網(wǎng)格點(diǎn)移動的新位置。將三角形網(wǎng)格類比為結(jié)構(gòu)力學(xué)中的二維桁架結(jié)構(gòu),三角形邊比作桁架桿件,三角形頂點(diǎn)比作桁架節(jié)點(diǎn),力使桿件變形即節(jié)點(diǎn)移動。類比僅為使過程便于理解,以下討論和方程不涉及實(shí)際的物理意義和量綱。由節(jié)點(diǎn)的當(dāng)前位置pn和力f(pn),解桁架結(jié)構(gòu)靜力平衡方程,得到節(jié)點(diǎn)的移動位置：

定義桁架上的力：

式中,l為點(diǎn)處于當(dāng)前位置pn時的已知現(xiàn)有邊長,l0為點(diǎn)處于新位置pn+1時的期望邊長；為使點(diǎn)擴(kuò)散,定義f恒為正；當(dāng)l接近l0時,f接近0。期望邊長計(jì)算式：

式中,rb為由分辨率函數(shù)給出的現(xiàn)有三角形邊中點(diǎn)的目標(biāo)分辨率；由于f恒為正,因此設(shè)置系數(shù)z=1.2[12]；代表現(xiàn)有邊長和目標(biāo)邊長的比例關(guān)系。

1.2 分辨率函數(shù)及評價指標(biāo)

實(shí)際問題的復(fù)雜各異性產(chǎn)生了對目標(biāo)網(wǎng)格的特定分辨率要求。針對海洋數(shù)值計(jì)算網(wǎng)格的生成,構(gòu)建全局分辨率函數(shù)和局部加密函數(shù),分別對整個區(qū)域和局部區(qū)域進(jìn)行分辨率控制。全局分辨率函數(shù)確定網(wǎng)格單元分布的宏觀要求,即近岸生成小單元,隨離岸距離變大,生成單元擴(kuò)大；局部加密函數(shù)在地形突變處(如存在陸坡、深槽等),修正全局分辨率函數(shù)的結(jié)果(如分辨率過低),以更準(zhǔn)確地描述其復(fù)雜的流體動力變化過程,減小模型計(jì)算誤差、避免模型計(jì)算溢出中斷。為了確保分辨率函數(shù)控制的網(wǎng)格迭代能持續(xù)進(jìn)行至生成網(wǎng)格質(zhì)量良好,適用于模型計(jì)算,提出分辨率函數(shù)合理性的評價指標(biāo)。網(wǎng)格迭代至質(zhì)量滿足要求且繼續(xù)迭代的質(zhì)量提升效果不大時,為了合理中止其迭代,提出網(wǎng)格質(zhì)量評價指標(biāo),對網(wǎng)格質(zhì)量進(jìn)行綜合評估,作為網(wǎng)格迭代的中止判據(jù)。

1.2.1 全局分辨率函數(shù)

全局分辨率函數(shù)中,自變量為網(wǎng)格點(diǎn)與其距離最近的岸線點(diǎn)之間的距離：

式中,(x,y)、(xc,yc)分別為網(wǎng)格點(diǎn)、最近岸線點(diǎn)的位置坐標(biāo),dm為距離上限。因變量為網(wǎng)格點(diǎn)分辨率：

1.2.2 局部加密函數(shù)

局部加密函數(shù)中,自變量為加密區(qū)域網(wǎng)格點(diǎn)與其距離最近的加密點(diǎn)之間的距離：

式中,(x,y)、(xr,yr)分別為網(wǎng)格點(diǎn)、最近加密點(diǎn)的位置坐標(biāo),d4為加密區(qū)域距離上限。因變量為加密區(qū)域網(wǎng)格點(diǎn)分辨率：

1.2.3 分辨率函數(shù)合理性評價指標(biāo)

曲線越光滑、疏密過渡越均勻的網(wǎng)格質(zhì)量越好,但通常其網(wǎng)格規(guī)模較大且相鄰單元大小更接近,網(wǎng)格迭代和模型計(jì)算效率較低；反之,網(wǎng)格規(guī)模的減小和相鄰單元大小變化較大時,迭代中網(wǎng)格質(zhì)量將難以提升,模型計(jì)算誤差由此產(chǎn)生。應(yīng)綜合考慮網(wǎng)格質(zhì)量和整體效率構(gòu)建合理的分辨率函數(shù),引入評價指標(biāo)：

式中,r為分辨率函數(shù)設(shè)計(jì)的目標(biāo)網(wǎng)格中的分辨率集合,eac指目標(biāo)網(wǎng)格中相鄰三角形面積改變的理論最大值。eac越接近0,目標(biāo)網(wǎng)格中疏密過渡越均勻。認(rèn)為eac∈[0,0.5]時,網(wǎng)格迭代持續(xù)進(jìn)行至生成的網(wǎng)格質(zhì)量良好。但eac越小,整體效率越低,因此,實(shí)際應(yīng)用中,eac值不應(yīng)過小。

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

網(wǎng)格質(zhì)量,包括網(wǎng)格均勻性和單元質(zhì)量。整個網(wǎng)格的現(xiàn)有分辨率越接近目標(biāo)分辨率,網(wǎng)格均勻性越好。對均勻性的質(zhì)量評價,采用三角形外接圓半徑與圓心處目標(biāo)分辨率比值的標(biāo)準(zhǔn)偏差,其值越小,質(zhì)量越好,認(rèn)為其值小于0.05時,均勻性質(zhì)量滿足要求。單個三角形的質(zhì)量：

式中,rin為三角形內(nèi)切圓半徑,rout為三角形外接圓半徑,a,b,c為三角形邊長；Q值越大,質(zhì)量越好；Q=1,即三角形為等邊三角形時,質(zhì)量最好,Q＞0.5時,三角形質(zhì)量滿足要求。

采用以下質(zhì)量評價指標(biāo)對網(wǎng)格質(zhì)量進(jìn)行綜合評估。網(wǎng)格均勻性質(zhì)量Quni描述整個網(wǎng)格的現(xiàn)有分辨率接近目標(biāo)分辨率的程度,其值越接近0,質(zhì)量越好,表達(dá)式：

單元平均質(zhì)量Qmean指網(wǎng)格中所有三角形質(zhì)量的平均值,其值越接近1,質(zhì)量越好,表達(dá)式：

單元最小質(zhì)量Qmin指網(wǎng)格中質(zhì)量最差的一個三角形的質(zhì)量,其值越接近1,質(zhì)量越好,表達(dá)式：

人工調(diào)整比例Q05指網(wǎng)格中質(zhì)量小于0.5即需人工調(diào)整的三角形個數(shù)占所有三角形總個數(shù)的比例,其值越接近0,質(zhì)量越好,表達(dá)式：

式(10)～式(13)中,(rout)i為第i個三角形的外接圓半徑,(r)i為第i個三角形的外接圓圓心處的目標(biāo)分辨率,Qi為第i個三角形的質(zhì)量,N為網(wǎng)格中所有三角形的總個數(shù),n為網(wǎng)格中需人工調(diào)整的三角形個數(shù)。

2 應(yīng)用實(shí)例

應(yīng)用本文方法時,首先參考eac構(gòu)建出合理的分辨率函數(shù),之后全局分辨率函數(shù)對整個區(qū)域的分辨率進(jìn)行宏觀控制,需要局部加密時,再結(jié)合局部加密函數(shù)對局部區(qū)域分辨率進(jìn)行控制,并保證加密區(qū)與附近區(qū)域分辨率的自然過渡。合理的分辨率函數(shù)保證網(wǎng)格迭代的進(jìn)行,但在網(wǎng)格質(zhì)量達(dá)到良好且提升極其緩慢的情況下,迭代可能仍在進(jìn)行,因此將網(wǎng)格質(zhì)量評價指標(biāo)作為網(wǎng)格中止迭代的判據(jù),選取可用網(wǎng)格和保證整體效率。將本文方法應(yīng)用于理想潮汐汊道和崇明東灘復(fù)雜潮溝系統(tǒng),其網(wǎng)格生成分別利用了全局分辨率函數(shù)和全局分辨率函數(shù)、局部加密函數(shù)相結(jié)合的手段。

2.1 理想潮汐汊道

理想潮汐汊道是近岸地貌動力學(xué)數(shù)值模擬中的典型理想模型[16],其為長15 km(X方向),寬14 km(Y方向)的矩形區(qū)域,北部為開放海域,南部為潟湖,由寬2 km的潮汐通道連通,外邊界由北側(cè)開邊界和其余岸線邊界組成。為探討本文方法在不同的全局分辨率函數(shù)下生成同一區(qū)域的網(wǎng)格的情況,在全局分辨率函數(shù)合理的前提下,對理想潮汐汊道構(gòu)建3個參數(shù)設(shè)置不同的全局分辨率函數(shù)type 1,type 2和type 3(表1),其中,預(yù)設(shè)岸線點(diǎn)分辨率rc均為200 m。

由eac(表1)判斷,在3個分辨率函數(shù)下,網(wǎng)格迭代均能持續(xù)進(jìn)行至生成的網(wǎng)格質(zhì)量良好,且網(wǎng)格中疏密過渡均勻程度：type 3＞type 2＞type 1。網(wǎng)格點(diǎn)與其最近岸線點(diǎn)的實(shí)際最大離岸距離約為7 km(圖3a),在此范圍內(nèi),網(wǎng)格點(diǎn)分辨率隨離岸距離的3種函數(shù)變化關(guān)系(圖2)：隨離岸距離的增大,type 1的分辨率迅速增大,之后保持較為平緩的增長；type 2的分辨率始終較為勻速的增長；type 3的分辨率經(jīng)歷了由緩慢增長至較快增長的變化。這代表了空間分布特征各異的3種網(wǎng)格分辨率要求,其單元的大小和分布上的差異表現(xiàn)為：type 1的網(wǎng)格高分辨率區(qū)的范圍較小,迅速過渡到低分辨率區(qū),過渡區(qū)的相鄰單元面積變化較快(圖3b)；type 2的網(wǎng)格相鄰單元的面積變化速率較平穩(wěn)(圖3c)；type 3的網(wǎng)格高分辨率區(qū)的范圍較大,相鄰單元面積變化緩慢,過渡到低分辨率區(qū)后,相鄰單元面積變化較快(圖3d)。

表1 理想潮汐汊道全局分辨率函數(shù)參數(shù)和評價Table 1 The parameters and evaluation of global resolution function in ideal tidal inlet

在各函數(shù)下迭代生成網(wǎng)格(圖4),對同一分辨率函數(shù),迭代10,50,150次的網(wǎng)格分辨率均能體現(xiàn)相似的空間分布特征。type 1函數(shù)下,向岸小范圍區(qū)域內(nèi)分布著小單元,離岸大范圍區(qū)域內(nèi)分布著大單元,兩者具有較為明顯的分界,分界處迅速由小單元過渡到大單元(圖4a,4b,4c)；type 2函數(shù)下,整個區(qū)域中相鄰單元面積變化平穩(wěn),疏密過渡自然(圖4d,4e,4f)；type 3函數(shù)下,向岸大范圍區(qū)域內(nèi)分布著小單元,離岸小范圍區(qū)域內(nèi)分布著大單元,疏密過渡較均勻(圖4g,4h,4i)。

各函數(shù)下的網(wǎng)格迭代隨時間的變化(表2、圖5)均表現(xiàn)出較為相似的規(guī)律。隨迭代次數(shù)的增大,Qmean不斷增大,前期增速迅速,迭代150次時,各函數(shù)下的迭代總時間均小于0.25 min,Qmean均大于0.960(表2),之后增速變緩,迭代約400次時均至0.975左右(圖5)；Q05保持減小趨勢,迭代150次時,Q05均小于0.08%(表2),即存在約5個需人工質(zhì)量調(diào)整的三角形,迭代約500次之后均為0(圖5)；由于迭代隨機(jī)性的影響,前期Qmin和Quni存在較大波動變化(表2、圖5),但總體上,Qmin呈增大趨勢,迭代約400次時,其均值0.7左右,Quni呈減小趨勢,迭代約400次時,其均值0.02左右(圖5)；整體網(wǎng)格質(zhì)量呈升高趨勢。

圖2 理想潮汐汊道各函數(shù)下分辨率隨離岸距離變化關(guān)系Fig.2 Relationship between resolution and offshore distance under each function in ideal tidal inlet

圖3 理想潮汐汊道的距離和3個分辨率函數(shù)值Fig.3 Distance and three resolution function values in ideal tidal inlet

圖4 理想潮汐汊道各函數(shù)下迭代10,50和150次的網(wǎng)格Fig.4 The meshes after 10,50,and 150 iterations under each function in ideal tidal inlet

在全局分辨率函數(shù)合理的前提下,本文方法在不同的全局分辨率函數(shù)下生成理想潮汐汊道的網(wǎng)格,其結(jié)果顯示,eac能作為全局分辨率函數(shù)合理性的有效判據(jù)；本文方法對不同的分辨率要求均表現(xiàn)出良好的適應(yīng)性,在相對較短的時間內(nèi),獲得了質(zhì)量較好的網(wǎng)格。網(wǎng)格迭代次數(shù)相同時,迭代總時間、節(jié)點(diǎn)數(shù)、單元數(shù)和Qmean基本遵循：type 3＞type 2＞type 1(表2、圖5),表明函數(shù)相應(yīng)的網(wǎng)格節(jié)點(diǎn)、單元數(shù)越多,迭代總時間越長,其網(wǎng)格質(zhì)量越高；但各函數(shù)下Qmean均達(dá)到較高水平,其間的差距非常微小。因此使用本文方法時,可根據(jù)實(shí)際問題,在合理的前提下自由地構(gòu)建合適的分辨率函數(shù)。

表2 理想潮汐汊道各函數(shù)下迭代10,50,150次時的情況Table 2 The evaluation after 10,50,and 150 iterations under each function in ideal tidal inlet

圖5 理想潮汐汊道各函數(shù)下網(wǎng)格質(zhì)量隨迭代次數(shù)變化的過程Fig.5 The mesh quality evolution under each function in ideal tidal inlet

2.2 崇明東灘潮溝系統(tǒng)

潮溝作為潮間帶和相鄰水域進(jìn)行物質(zhì)和能量交換的重要通道,對光灘和鹽沼的沖淤過程和地貌形態(tài)起著非常重要的控制作用[17],同時對生物遷徙、魚類洄游等的生態(tài)環(huán)境產(chǎn)生重要影響。崇明島位于長江口入海口,崇明東灘位于崇明島的東端,灘面遼闊,這里發(fā)育了數(shù)量眾多、地形復(fù)雜的網(wǎng)絡(luò)狀潮溝(圖6c)。崇明東灘數(shù)值模擬研究需要精確地擬合其復(fù)雜的潮溝地形,大量的網(wǎng)格加密工作由此產(chǎn)生。為探討本文方法擬合復(fù)雜地形生成高分辨率網(wǎng)格的情況,在分辨率函數(shù)合理的前提下,對崇明東灘構(gòu)建全局分辨率函數(shù)和局部加密函數(shù),全局分辨率函數(shù)參數(shù)為dm=30 000,ym=1 000,fr=10,fd=30,s1=3,s2=6,此時eac=0.091。將岸線點(diǎn)分辨率預(yù)設(shè)在40～500 m范圍內(nèi),其中,近潮溝岸線的分辨率在40～50 m范圍內(nèi),潮溝加密點(diǎn)的最高分辨率為6 m,由此取h0=7.5 m。取s=1.14,并取h=0.53h0=3.975 m。在d3=3h=11.925m時,eac=0.488,分辨率變化速率接近限值,保持d3處的分辨率變化速率至d4=170h=675.750 m,此處分辨率接近附近非加密區(qū)域分辨率。

由eac判斷,網(wǎng)格迭代能有效進(jìn)行并使生成的網(wǎng)格質(zhì)量良好,高eac值的局部加密函數(shù)要求潮溝與其附近的分辨率差異明顯(圖7d),對潮溝精確模擬的同時避免了網(wǎng)格規(guī)模劇增。網(wǎng)格點(diǎn)與其最近岸線點(diǎn)的實(shí)際最大離岸距離約為30 km(圖7a)。分辨率函數(shù)對單元的大小和分布的控制要求為：整個區(qū)域中相鄰單元面積變化較平緩,近距高分辨率區(qū)的范圍非常小,過渡到遠(yuǎn)距低分辨率區(qū)后,相鄰單元面積變化加快(圖7b)；加密區(qū)域中,加速增長區(qū)范圍非常小,迅速過渡到勻速增長區(qū),其過渡處的相鄰單元面積變化非常快,并在勻速增長區(qū),保持該面積變化速度(圖7c)。

圖6 長江、崇明東灘和潮溝區(qū)域水深Fig.6 Water depth of the Changjiang River,Chongming Dongtan and tidal creek area

圖7 崇明東灘的距離和分辨率函數(shù)值Fig.7 Distance and resolution function values in Chongming Dongtan

迭代網(wǎng)格并對網(wǎng)格質(zhì)量(圖8)進(jìn)行評估。隨迭代次數(shù)增大,Qmean增大,Q05減小,受迭代隨機(jī)性的影響,Qmin和Quni有波動變化,但分別具有增大和減小趨勢,整體上,網(wǎng)格質(zhì)量呈增高趨勢；迭代105次(迭代總時間127 min)時,Qmean達(dá)到了0.960,Qmin接近0；迭代120次(迭代總時間128 min)時,Quni為0.020,此時Qmin仍較小,但Q05接近0；迭代330次(迭代總時間141 min)時,Qmin達(dá)到了0.655 6,即不存在需人工調(diào)整的三角形。生成的網(wǎng)格(迭代330次)在非加密區(qū)域中的相鄰單元面積變化平緩,疏密過渡自然(圖9a)；潮溝能明顯的從附近區(qū)域中區(qū)分出來,潮溝與其附近的疏密過渡自然(圖9b,9c)。

圖8 崇明東灘網(wǎng)格質(zhì)量隨迭代次數(shù)變化的過程Fig.8 The mesh quality evolution under each function in Chongming Dongtan

圖9 崇明東灘生成網(wǎng)格圖Fig.9 The mesh of Chongming Dongtan

崇明東灘數(shù)值模擬研究需在長江口大環(huán)境中考慮,而潮溝處的高分辨率網(wǎng)格限制了模型時間步長,考慮到模型計(jì)算效率,本文采用模型嵌套網(wǎng)格技術(shù),將崇明東灘小區(qū)(圖6b)模型嵌套于丁平興等已經(jīng)發(fā)展非常成熟的長江口大區(qū)(圖6a)模型[18]中。崇明東灘網(wǎng)格迭代至330次(迭代總時間141 min)時不存在需人工調(diào)整的三角形,選取該網(wǎng)格進(jìn)行處理后用于模型計(jì)算。計(jì)算網(wǎng)格有101 247個節(jié)點(diǎn)和201 711個三角形,最高分辨率為7 m。小區(qū)模型垂向上采用σ坐標(biāo),均勻分11層,采用0.1 s的時間步長,計(jì)算22 d。在大區(qū)模型中,長江口上界取至潮流界江陰以上200 km左右的大通,錢塘江上界取至老鹽倉附近。外海開邊界東至124°30′E,北至34°30′N,南至28°30′N。整個區(qū)域覆蓋了長江各個入海口。網(wǎng)格分辨率在外海陸架海域?yàn)?.5～5.0 km,在長江口則精細(xì)到200 m左右。更詳細(xì)的介紹可參考[18]。大區(qū)模型垂向上采用σ坐標(biāo),均勻分11層,采用1 s的時間步長,計(jì)算26 d。流場在外海通過水位作為模型的驅(qū)動,采用8個主要分潮M2,S2,K1,O1,N2,K2,P1和Q1。由于水位、流速的響應(yīng)較快,模型采用冷啟動方式,即設(shè)置初始場為0。

模型計(jì)算流場的最大淹沒時刻(圖10a)淹沒區(qū)域可到達(dá)離岸150 m處,東灘外的開闊海域呈現(xiàn)顯著的落潮流特征。在崇明東灘的潮溝系統(tǒng)(圖10b)中,潮溝與灘面這兩種地形結(jié)構(gòu)的流場特征存在顯著差異,由于水流漫灘距離短,灘面流主要發(fā)生在離岸較遠(yuǎn)的灘面下緣,其流速較小,最大流速約為0.2 m/s；由于灘面水深變化緩慢,流速的空間差異也較小。潮溝處(圖10c)的流則主要受到潮溝地形的約束,基本沿著潮溝方向流動,其流速較大,最大流速可達(dá)0.5 m/s；潮灘中流速的空間差異較大,主要與水深有關(guān),一般水深越深,流速越大。

圖10 崇明東灘嵌套模型流場Fig.10 The hydrodynamic field in Chongming Dongtan by the nesting model

在分辨率函數(shù)合理的前提下,采用本文方法通過全局分辨率函數(shù)和局部加密函數(shù)相結(jié)合的手段,生成崇明東灘復(fù)雜潮溝系統(tǒng)的高分辨率網(wǎng)格并進(jìn)行模型計(jì)算,其結(jié)果顯示,本文方法能在相對較短的時間內(nèi)完成大量的網(wǎng)格加密工作擬合出復(fù)雜的地形,生成網(wǎng)格的分辨率較好地滿足了分辨率函數(shù)的要求,其網(wǎng)格質(zhì)量較好。利用質(zhì)量評價指標(biāo)能為中止網(wǎng)格迭代選取出計(jì)算網(wǎng)格提供可靠參考,保證了迭代效率。由于網(wǎng)格能在流速較大且空間變化顯著的潮溝系統(tǒng)中進(jìn)行加密,而在流速較小且空間變化較小的灘面上采用較少的網(wǎng)格,在保證精確模擬流場的同時,減小了模型網(wǎng)格規(guī)模,提高了計(jì)算效率。對比其他方法,如在灘面上采用和潮溝相同的空間分辨率,則網(wǎng)格節(jié)點(diǎn)數(shù)約為220 000,增加2.2倍,這將顯著影響模型的計(jì)算效率,而將人工網(wǎng)格調(diào)整手段運(yùn)用在極其復(fù)雜的潮溝系統(tǒng)中,調(diào)整難度較大,需要大量的生成調(diào)試過程。在實(shí)際復(fù)雜問題應(yīng)用中,本文方法實(shí)現(xiàn)了對網(wǎng)格單元的大小和分布的全面精確的控制,有效解決了對網(wǎng)格分辨率特定的控制需求,高效地獲得了高精度、高質(zhì)量的網(wǎng)格。

3 結(jié) 語

針對海洋數(shù)值模型的空間離散,基于DISTMESH,本研究提出了一種能精確擬合復(fù)雜岸線、地形的三角形網(wǎng)格自適應(yīng)生成方法。該方法保留了DISTMESH方法生成網(wǎng)格質(zhì)量高、速度快的優(yōu)點(diǎn),并將DISTMESH方法中分辨率函數(shù)構(gòu)建為適應(yīng)于海洋數(shù)值模型的全局分辨率函數(shù)和局部加密函數(shù),可以對網(wǎng)格生成的空間布局進(jìn)行多種類型的預(yù)先設(shè)置,獲得了對分辨率全面的控制；eac評價指標(biāo)確定分辨率函數(shù)的合理性,保證迭代的有效開始；多個質(zhì)量評價指標(biāo)對網(wǎng)格質(zhì)量進(jìn)行評估,合理中止迭代,以提高網(wǎng)格生成效率。將該方法應(yīng)用于理想和實(shí)際案例的網(wǎng)格生成,案例表明,本文方法能滿足理想和實(shí)際條件下海洋數(shù)值模型的空間離散要求,高效地生成高精度、高質(zhì)量的三角形網(wǎng)格。在實(shí)際復(fù)雜案例的崇明東灘潮溝系統(tǒng)的應(yīng)用中,能很好地刻畫高復(fù)雜度的潮溝網(wǎng)格,并保證了較高的網(wǎng)格質(zhì)量,從而保證了數(shù)值模型的在此復(fù)雜邊界條件下的可靠計(jì)算。