李付鵬，汪繼文，歐 莽

（1.安徽大學(xué)計(jì)算智能與信號(hào)處理教育部重點(diǎn)實(shí)驗(yàn)室，安徽合肥 230039；2.安徽大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，安徽合肥 230039）

文章編號(hào)：1001?246X（2015）01?0033?07

一種光滑粒子流體動(dòng)力學(xué)數(shù)據(jù)的后處理方法

李付鵬1，2，汪繼文1，2，歐 莽2

（1.安徽大學(xué)計(jì)算智能與信號(hào)處理教育部重點(diǎn)實(shí)驗(yàn)室，安徽合肥 230039；
2.安徽大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，安徽合肥 230039）

提出一種光滑流體動(dòng)力學(xué)數(shù)據(jù)后處理的三角網(wǎng)格化方法.首先確定流場中每個(gè)粒子作用域內(nèi)的配對(duì)粒子集；再以當(dāng)前粒子的配對(duì)粒子集為限定條件進(jìn)行Delaunay三角剖分，可得到由粒子作為節(jié)點(diǎn)的三角形單元集；最后對(duì)“虛假”單元做過濾處理.算例表明復(fù)雜邊界條件和嚴(yán)重變形的流體區(qū)域上的SPH計(jì)算結(jié)果都可用該方法處理.

光滑粒子方法；非凸區(qū)域；Delaunay三角化；后處理

0 引言

近年來，光滑粒子流體動(dòng)力學(xué)（smoothed particle hydrodynamics，簡稱為SPH）方法［1］逐漸成為解決物理問題的一個(gè)有效工具，與傳統(tǒng)的基于網(wǎng)格的數(shù)值方法，如有限差分、有限元、有限體積方法不同，這是一種典型的無網(wǎng)格拉格朗日數(shù)值方法，它利用核函數(shù)對(duì)物理問題進(jìn)行近似處理，用離散的粒子來描述宏觀連續(xù)分布微觀仍為粒子的流體，而每個(gè)粒子則攜帶了其所在位置的流體的各種性質(zhì)，如質(zhì)量、密度、速度、能量等，算法節(jié)點(diǎn)可以任意分布，并非有某種網(wǎng)格聯(lián)系在一起，這種自適應(yīng)的特點(diǎn)使得該方法在模擬流體劇烈變形或間斷問題時(shí)具有顯著優(yōu)勢，目前方法的應(yīng)用已擴(kuò)展到氣體動(dòng)力學(xué)、不可壓縮、爆炸、固體力學(xué)和彈性體等多個(gè)領(lǐng)域.然而，與該算法應(yīng)用快速發(fā)展所不相適應(yīng)的是，方法的數(shù)據(jù)后處理技術(shù)一直沒有太大的發(fā)展，甚至可以說是一大難題［2］.如果將計(jì)算后的結(jié)果直接圖形化，得到的是大量無組織的離散粒子，無法表現(xiàn)出流體的性態(tài)，更無法對(duì)流體的等值線（面）、流線圖等分析，現(xiàn)有的比較成熟的基于網(wǎng)格的可視化技術(shù)也無法直接利用.目前，就我們查閱到的文獻(xiàn)而言，有基于Splash、Paraview等軟件開發(fā)的SPH后處理程序［3］，但這類軟件較為專業(yè)，通用性不好，同時(shí)在流體力學(xué)方面的處理能力有限.

文獻(xiàn)［4－5］提出了一種處理SPH數(shù)據(jù)后處理的方法，該方法先利用Delaunay三角剖分技術(shù)對(duì)離散數(shù)據(jù)三角化，然后根據(jù)每個(gè)單元中三個(gè)節(jié)點(diǎn)是否彼此為粒子作用對(duì)，決定是否將該單元從單元集中刪除，從算例來看，這種方法是可行的.本文根據(jù)SPH算法的特點(diǎn)，對(duì)上述方法進(jìn)行了改進(jìn)，提出了一種更為簡便的SPH數(shù)據(jù)后處理方法，以粒子支持域（影響域）為限定條件進(jìn)行限定三角剖分（constrained Delaunay triangulation，簡稱為CDT或conforming Delaunay triangulation，簡稱為RDT），與文獻(xiàn)［4－5］相比，方法減少了所生成的多余的三角形網(wǎng)格的數(shù)量，刪除的空白的三角形網(wǎng)格也隨之減少，提高了算法運(yùn)算的效率.這樣做另一個(gè)優(yōu)點(diǎn)在于三角網(wǎng)格化輸出的圖形能夠清晰地體現(xiàn)各粒子之間的影響關(guān)系，凡與某個(gè)節(jié)點(diǎn)（粒子）相關(guān)聯(lián)的其它節(jié)點(diǎn)（粒子）均為數(shù)據(jù)輸出時(shí)刻與該粒子有相互影響的粒子，而且方法更加簡便易行，一般情況下不需要過濾“虛假”單元，簡化了算法流程，增強(qiáng)了算法的適應(yīng)性，驗(yàn)證結(jié)果表明方法穩(wěn)定性好，凸區(qū)域和嚴(yán)重變形的非凸區(qū)域上的SPH計(jì)算結(jié)果都可用該方法處理.

首先介紹SPH方法中與數(shù)據(jù)三角化相關(guān)的幾個(gè)概念以及數(shù)據(jù)后處理的特點(diǎn)；接下來闡述SPH數(shù)據(jù)三角化的算法思想和步驟；最后通過算例驗(yàn)證，得出結(jié)論.

1 SPH粒子支持域和粒子配對(duì)

SPH方法近似偏微分方程的過程包含兩個(gè)步驟：核近似和粒子近似，核近似是通過對(duì)場函數(shù)及權(quán)函數(shù)積分實(shí)現(xiàn)，粒子近似是在一個(gè)有限區(qū)域內(nèi)對(duì)所有粒子進(jìn)行加權(quán)求和實(shí)現(xiàn).對(duì)于場函數(shù)f（x），核近似和粒子近似可以分別表示為

其中變量x是位置矢量，變量h被稱為光滑長度，W（x-x′，h）即為光滑核函數(shù)，變量h用來確定核函數(shù)的支持域，是核寬度的一種度量，核近似算子習(xí)慣用角括?。迹緲?biāo)記.場函數(shù)的導(dǎo)數(shù)可以采取與上述類似的方法求出，將場函數(shù)與其導(dǎo)數(shù)的粒子近似公式與實(shí)際應(yīng)用的方程相結(jié)合，便可以求解特定的物理問題.

從以上公式可以看出，光滑核函數(shù)在SPH近似方法中起著重要作用，它決定了函數(shù)表達(dá)式的精度和計(jì)算效率.在核函數(shù)的表達(dá)式中有兩個(gè)參數(shù)，分別為粒子的位置矢量和光滑長度，對(duì)于某個(gè)特定的粒子i而言，其周圍粒子與該粒子的距離｜xi－x′｜＞κh時(shí)，W（xi－x′，h）＝0，式中κ是與點(diǎn)x處光滑函數(shù)相關(guān)的常數(shù)，此時(shí)場函數(shù)f（x）＝0，就是說，粒子i周圍的粒子對(duì)粒子i沒有影響；僅當(dāng)" ｜x-x′｜≤κh時(shí)，W（xi－x′，h）≠0，粒子i周圍的粒子對(duì)粒子i產(chǎn)生影響．SPH定義｜ x-x′｜≤κh所確定的以粒子i位置為中心的以κh為半徑的區(qū)域?yàn)榱Ｗ觟支持域．粒子i支持域內(nèi)的所有點(diǎn)的信息都被用來決定粒子i的信息，支持域外的信息則對(duì)該點(diǎn)沒有影響．支持域與粒子的光滑長度有關(guān)，光滑長度h與因子κ的乘積決定了該粒子的支持域．對(duì)于兩個(gè)光滑長度相同的配對(duì)粒子i和j來說，粒子i和j彼此均在對(duì)方的支持域內(nèi).

給出了支持域的概念，可以很容易地確定粒子的配對(duì)粒子，對(duì)于給定的粒子，其支持域內(nèi)的所有粒子均為該粒子的配對(duì)粒子.如果整個(gè)算法過程中光滑長度為常數(shù)，配對(duì)粒子互相影響，彼此互為配對(duì)粒子，而對(duì)于可變光滑長度的方法，配對(duì)粒子可能是單向的，即一個(gè)粒子是另一個(gè)粒子的配對(duì)粒子；反之未必成立，也就是說可能出現(xiàn)其中一個(gè)粒子對(duì)另一個(gè)粒子發(fā)生作用，但不是相互作用的，這違背牛頓第三運(yùn)動(dòng)定律，因此必須設(shè)法保證粒子間相互作用的對(duì)稱性，否則會(huì)導(dǎo)致流場不守恒.

SPH支持域內(nèi)粒子搜索算法常用的有全配對(duì)搜索法、鏈表搜索法、樹形搜索法.

2 SPH數(shù)據(jù)后處理

對(duì)于基于網(wǎng)格的數(shù)值方法，生成網(wǎng)格的目的在于為數(shù)值計(jì)算提供一個(gè)高質(zhì)量的網(wǎng)格，以便提高數(shù)值計(jì)算的精度，網(wǎng)格主要的目的是用于數(shù)值計(jì)算，盡管也包括數(shù)據(jù)結(jié)果的顯示和數(shù)值分析；而無網(wǎng)格算法生成網(wǎng)格的目的僅僅在于數(shù)據(jù)結(jié)果的顯示和數(shù)值分析.

兩類數(shù)值方法的數(shù)據(jù)處理的目的不同，生成網(wǎng)格質(zhì)量的要求也不相同.基于網(wǎng)格的數(shù)值方法要求三角剖分必須滿足Delaunay優(yōu)化準(zhǔn)則，對(duì)于流場變量梯度變化較為劇烈的局部區(qū)域還需要進(jìn)行加密處理，這樣才能得到較好質(zhì)量的計(jì)算網(wǎng)格；而SPH算法雖然同樣要求三角剖分滿足Delaunay優(yōu)化準(zhǔn)則，但還有一個(gè)更為重要的限定條件，就是三角化后的三角形單元的三條邊應(yīng)當(dāng)是配對(duì)邊，即構(gòu)成該三角形的三個(gè)頂點(diǎn)（粒子）在互相的作用域之內(nèi)，只有彼此相互作用或影響的粒子對(duì)相連（即構(gòu)成三角形的一條邊）才有意義，否則會(huì)導(dǎo)致流場的“非物質(zhì)”區(qū)域被剖分的現(xiàn)象，兩個(gè)粒子無互相作用，就表明它們之間沒有物質(zhì)接觸，即沒有質(zhì)量.如果一個(gè)三角形單元中有一對(duì)節(jié)點(diǎn)（由粒子構(gòu)成的三角形頂點(diǎn)）之間沒有物質(zhì)接觸，就表明這個(gè)單元內(nèi)部有間斷，因此從連續(xù)介質(zhì)力學(xué)的觀點(diǎn)出發(fā)，該單元就是沒有真實(shí)的流場，因此不是配對(duì)的粒子就沒有必要生成三角形單元；因此SPH數(shù)據(jù)的三角剖分屬于限定三角剖分.

常見的限定三角剖分有CDT和RDT兩種類型.CDT是指剖分域的邊界或內(nèi)部限定邊在三角剖分結(jié)果中出現(xiàn)并且不允許被細(xì)分，既不允許在剖分域中加入額外的點(diǎn)；RDT允許在域中加入點(diǎn)從而使得域的所有邊界存在于剖分結(jié)果中.對(duì)SPH粒子點(diǎn)數(shù)據(jù)集可以用CDT或RDT方法三角化，但意義不同，SPH粒子點(diǎn)具有描述所求解的問題應(yīng)具有的各物理量（質(zhì)量、密度、壓力等），因此如果數(shù)據(jù)集用于數(shù)值計(jì)算，就必須按照CDT的方法剖分，圖形化的數(shù)據(jù)保持原有場域的各物理量不變；如果數(shù)據(jù)集用于視覺效果（比如游戲等），就可以按照RDT方法三角化，加入些許虛擬點(diǎn)，增強(qiáng)視覺效果.兩個(gè)不同的方法必然產(chǎn)生不同質(zhì)量的三角化的網(wǎng)格，CDT方法由于不允許加點(diǎn)，可能造成SPH數(shù)據(jù)集邊界附近的三角形單元不滿足Delaunay優(yōu)化準(zhǔn)則.RDT允許加入額外的點(diǎn)，SPH數(shù)據(jù)集域內(nèi)的所有單元均符合Delaunay優(yōu)化準(zhǔn)則.

3 SPH粒子三角化算法

算法以SPH粒子支持域內(nèi)的配對(duì)粒子為約束條件，以Delaunay三角剖分算法為基礎(chǔ)生成二維三角化網(wǎng)格，剖分時(shí)Delaunay優(yōu)化準(zhǔn)則（最小角最大化或空外接圓特性）檢查及其換邊操作僅限定在配對(duì)粒子集合中.

設(shè)SPH數(shù)據(jù)輸出時(shí)刻，SPH粒子的集合為P，三角化后的粒子集合為Q，算法具體描述如下：

1）首先根據(jù)粒子的支持域概念調(diào)用搜索算法檢索粒子集合P，確定每個(gè)粒子的配對(duì)集合；

2）對(duì)配對(duì)后的粒子集合按照粒子的坐標(biāo)排序，生成有序粒子點(diǎn)集；

3）從有序粒子點(diǎn)集中依序確定當(dāng)前粒子，調(diào)用Delaunay算法對(duì)當(dāng)前粒子及其配對(duì)集合中的粒子三角剖分，將參與三角剖分后的粒子依序加入到Q；

4）若當(dāng)前粒子的配對(duì)集合中的元素（粒子）都已三角化，則從已經(jīng)三角化的集合Q中依序確定新的當(dāng)前粒子；

5）依次檢測并刪除當(dāng)前粒子的配對(duì)粒子集合中與集合Q有交集的元素（粒子），對(duì)當(dāng)前粒子配對(duì)集合中存留的元素（粒子）三角化處理，并將三角化后的元素（粒子）添加到Q集合中；

6）若當(dāng)前粒子配對(duì)集合中的元素處理完畢，從Q集合中依序選取新的當(dāng)前粒子，重復(fù)步驟5，直至遍歷Q集合中的所有元素，算法結(jié)束.

需要說明的是：結(jié)束時(shí)，可能依然存在P中的部分元素（粒子）沒有被三角化的現(xiàn)象，原因主要有以下兩個(gè)方面：這部分粒子沒有與其他的粒子相互作用，多表現(xiàn)為由于流體碰撞飛濺出去的粒子；另一種可能是這部分粒子與其它粒子沒有物質(zhì)接觸，這部分粒子構(gòu)成的流體與其它流體之間存在著間斷或不連通的現(xiàn)象.此外，上述算法中最初參與三角剖分的粒子的選取很重要，如果最初參與剖分的是孤立粒子，由于孤立粒子周圍沒有配對(duì)粒子，將導(dǎo)致算法無法進(jìn)行下去.這時(shí)要將該粒子直接刪除，重新依序選擇新的參與剖分的粒子.

另一點(diǎn)需要說明的是：后處理的精細(xì)化問題.由于對(duì)流體局部細(xì)節(jié)表現(xiàn)的需要或者視覺效果的需要，需要將粒子占據(jù)的區(qū)域離散的更細(xì)，從而使得物理量場的變化能被更精細(xì)地反映出來.表現(xiàn)細(xì)節(jié)可以采取粒子分裂技術(shù)［6］，基本的原則是在對(duì)粒子及新增粒子的物理量重新分配時(shí)，要遵循質(zhì)量守恒和動(dòng)量守恒；如果是為了視覺的需要（比如游戲等），可考慮加入些許虛擬點(diǎn)，增強(qiáng)視覺效果，上述兩類處理都屬于RDT類型的三角剖分，由于允許加入額外的點(diǎn)，SPH數(shù)據(jù)集域內(nèi)的所有單元均符合Delaunay優(yōu)化準(zhǔn)則.后處理的精細(xì)化的方法［7］可以很容易地加入上述算法之中.

4 網(wǎng)格的質(zhì)量

在基于網(wǎng)格技術(shù)的數(shù)值方法中，都非常重視剖分網(wǎng)格的質(zhì)量，因?yàn)閿?shù)值方法的計(jì)算精度都會(huì)不同程度的受到網(wǎng)格質(zhì)量的影響，一般要求計(jì)算網(wǎng)格不能存在很小或很大的夾角.在二維網(wǎng)格中，就三角形剖分而言，一般用三角單位的外接圓半徑與內(nèi)切圓半徑的比值或者用三角形單位的外接圓半徑與最短邊長的比值來度量網(wǎng)格單元的質(zhì)量［8］，質(zhì)量越好的網(wǎng)格剖分，三角形單元越接近于正三角形.

而在無網(wǎng)格方法中，網(wǎng)格化的目的不是用于數(shù)值計(jì)算，而是用于數(shù)據(jù)后處理，所以網(wǎng)格化的圖形只要準(zhǔn)確反應(yīng)數(shù)值計(jì)算的結(jié)果就行了，對(duì)剖分的質(zhì)量要求也與基于網(wǎng)格化的數(shù)值方法不一樣，網(wǎng)格化后的圖形一方面要能夠表現(xiàn)出流場的運(yùn)動(dòng)的特征，另一方面又要能夠體現(xiàn)出粒子之間的相互關(guān)系，這是最重要的兩個(gè)目的.

5 算例

5.1 算例1

算法對(duì)文獻(xiàn)［4］描述的算例進(jìn)行了模擬.算例的計(jì)算域和初始粒子分布如圖1所示，初始時(shí)刻所有粒子靜止，沒有外力作用，圓形液滴以某一初始速度勻速向下運(yùn)動(dòng)，圖1－4顯示了流場中的粒子在初始時(shí)刻、碰撞、融合、飛濺等時(shí)刻的分布和三角化后對(duì)應(yīng)的各時(shí)刻的網(wǎng)格曲面圖.沒有對(duì)三角化后的網(wǎng)格作特殊處理，飛濺的粒子能夠被算法有效的“過濾”，沒有多余的單元，與文獻(xiàn)［4］相比，不需要再對(duì)三角化后的網(wǎng)格再次“過濾”刪除多余的單元.

圖1 初始時(shí)刻粒子分布和網(wǎng)格曲面Fig.1 At t＝0.0s particle distributions and element sets derived by Delaunay triangulation

圖2 2.5×10－3s時(shí)刻粒子分布和網(wǎng)格曲面Fig.2 At t＝2.5×10－3s particle distributions and element sets derived by Delaunay triangulation

圖3 5.0×10－3s時(shí)刻粒子分布和網(wǎng)格曲面Fig.3 At t＝5.0×10－3s particle distributions and element sets derived by Delaunay triangulation

圖4 12.5×10－3s時(shí)刻粒子分布和網(wǎng)格曲面Fig.4 At t＝12.5×10－3s particle distributions and element sets derived by Delaunay triangulation

5.2 算例2

二維圓潰壩問題.問題的求解區(qū)域是一個(gè)（0，50 m）×（0，50 m）的正方形，在正方形的中心是一個(gè)半徑為11 m的圓柱面形壩體.初始時(shí)，壩內(nèi)水深10 m，壩外水深1 m，且水是靜止的.研究在壩體突然潰決（完全消失）后水的徑向運(yùn)動(dòng)的規(guī)律.這一問題可以用來檢驗(yàn)方法保持對(duì)稱性的性能.SPH方法輸出的數(shù)據(jù)進(jìn)行網(wǎng)格化處理后再導(dǎo)入TECPLOT軟件，得到時(shí)間為0.69 s時(shí)的水流自由表面圖和水深等高線圖，分別如圖5和圖6所示.從圖5可知，本文所使用的SPH方法具有較好的對(duì)稱性，也較好地顯示了壩決后0.69 s時(shí)的水流形態(tài)，但對(duì)圖6做進(jìn)一步分析可知，在流場的局部區(qū)域，水深等高線與基于網(wǎng)格的方法［9－10］相比略顯雜亂，這說明本文所使用的SPH方法還需要進(jìn)一步改進(jìn).

圖5 水流自由表面Fig.5 Water surface profile after breaking of a circular dam

圖6 水深等高線Fig.6 Contours ofwater height after breaking of a circular dam

5.3 算例3

二維部分潰壩問題.問題的求解區(qū)域是一個(gè)（0，200m）×（0，200m）的正方形，正中有一條大壩將水庫和下游一分為二，上游水深10m，下游水深5m.大壩在距離水庫一邊95m處，突然有一段75m長的壩體潰決，要求模擬這部分壩體倒塌后水流的情況，該問題可以用來檢驗(yàn)方法對(duì)于復(fù)雜邊界的適應(yīng)性.問題的平面圖見圖7.

SPH計(jì)算的粒子分布見圖8.采用虛粒子來處理固定邊界條件.把網(wǎng)格化后的數(shù)據(jù)導(dǎo)入TECPLOT軟件，可得到壩體潰決后7.2s時(shí)的水流自由表面圖和水深等高線圖，如圖9、10所示.從圖可知，水流的整體形態(tài)和等高線的分布都與基于網(wǎng)格的數(shù)值方法［9，11］得出的結(jié)果基本保持一致，但在流體局部細(xì)節(jié)的表現(xiàn)上，本文所使用的SPH方法與傳統(tǒng)的數(shù)值方法計(jì)算的結(jié)果還有一定的差異，如下游壩體潰決的兩側(cè)所形成的漩渦沒有傳統(tǒng)數(shù)值方法明顯，邊界附近的計(jì)算結(jié)果也需要進(jìn)一步改善.這說明經(jīng)過網(wǎng)格化方法處理后的數(shù)據(jù)，能夠較為準(zhǔn)確地表現(xiàn)出所使用的數(shù)值方法的性能，同時(shí)也能發(fā)現(xiàn)數(shù)值方法所存在的問題.

圖7 二維部分潰壩問題Fig.7 Problem domain for partial dam break test

圖8 二維部分潰壩問題的粒子分布Fig.8 Particle distributions in partial dam break test

圖9 水流自由表面Fig.9 Water surface profile after breaking of dam

圖10 水深等高線Fig.10 Contours ofwater elevation in partial dam break test

上述給出的算例表明，本文提出的SPH數(shù)據(jù)后處理的方法，在處理二維流體力學(xué)數(shù)值問題時(shí)具有一定的可行性，對(duì)于較為復(fù)雜的邊界問題也能較好地處理.SPH數(shù)據(jù)通過網(wǎng)格化方法處理后，可以像基于網(wǎng)格的方法一樣通過TECPLOT流體力學(xué)軟件進(jìn)行數(shù)值分析，數(shù)值分析的結(jié)果準(zhǔn)確地反映了所使用的SPH算法的性能和存在的問題.

6 結(jié)論

提出一種光滑粒子流體動(dòng)力學(xué)的二維數(shù)據(jù)后處理的方法，以粒子支持域?yàn)橄薅l件對(duì)離散粒子三角剖分，對(duì)于物質(zhì)連通的流場區(qū)域，能夠自動(dòng)過濾飛濺的孤粒子，網(wǎng)格化后的數(shù)據(jù)一般不需要作特殊的處理，可以直接導(dǎo)入TECPLOT軟件進(jìn)行數(shù)值分析，算法簡單易于實(shí)現(xiàn).需要指出的是，本文不適合直接處理粒子飛濺較為嚴(yán)重的數(shù)值問題，這些飛濺的粒子可能會(huì)對(duì)流場的特性有較大影響，這類問題可采取粒子分裂方法，將粒子占據(jù)的區(qū)域離散成更多的粒子，盡可能地達(dá)到保留流場特性的目的.還需要指出的是，本文網(wǎng)格化后的數(shù)據(jù)表現(xiàn)的僅僅是流場的整體特征，無法表現(xiàn)SPH方法所特有的粒子飛濺、卷曲等細(xì)節(jié)特征.如何將網(wǎng)格化的方法與離散粒子相結(jié)合，圖形化后的數(shù)據(jù)既能體現(xiàn)流場的整體特征，又能體現(xiàn)流場飛濺、卷曲等細(xì)節(jié)特征.這些是我們下一步的工作.

［1］ Randles PW，Libersky L D.Smoothed particle hydrodynamics：Some recent improvements and applications［J］.Computer Methods in Applied Mechanics and Engineering，1996，139（1）：375－408.

［2］ Massidda，Luca.ARMANDO，a SPH code for CERN?Some theory，a short tutorial，the code description and some examples ［R］.European Organization for Nuclear Research，CERN?AB?Note?2008?039?ATB，2008：1－36.

［3］ Biddiscombe J，Graham D，Maruzewski P，et al.Visualization and analysis of SPH data［J］.ERCOFTAC Bulletin，SPH Special Edition，2008，76：9－12.

［4］ Zheng Jun，Yu Kaiping，Wei Yingjie，Zhang Jiazhong.A general study on post?processing of smoothed particle hydrodynamics ［J］.Chinese Journal of Computational Physics，2011，28（2）：213－218.

［5］ 鄭俊，于開平，張嘉鐘，魏英杰.非凸區(qū)域上SPH計(jì)算結(jié)果后處理方法研究［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2011，43（3）：12 －18.

［6］ Vacondio R，Rogers B D，Stansby P K.Accurate particle splitting for smoothed particle hydrodynamics in shallow water with shock capturing［J］.International Journal for Numerical Methods in Fluids，2012，69（8）：1377－1410.

［7］ Feldman J，Bonet J.Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems［J］. International Journal for Numerical Methods in Engineering，2007，72（3）：295－324.

［8］ 李海生.Delaunay三角剖分理論及可視化應(yīng)用研究［M］.哈爾濱：哈爾濱工業(yè)大學(xué)出版社，2010：36－38.

［9］ Anastasiou K，Chan C T.Solution of the 2D shallow water equations using the finite volumemethod on unstructured triangularmeshes［J］.International Journal for Numerical Methods in Fluids，1997，24（11）：1225－1245.

［10］ Alcrudo F，Garcia?Navarro P.A high?resolution Godunov?type scheme in finite volumes for the 2D shallow?water equations［J］. International Journal for Numerical Methods in Fluids，1993，16（6）：489－505.

［11］ Wang Jiwen，Liu Ruxun.Finite volume methods for solving the problem of discontinuous solution［J］.Chinese Journal of Computational Physics，2001，18（2）：97－105.

Article ID：1001?246X（2015）01?0040?11

Abstract： A numerical model，which solves equations with spectral element method and resolves atmospheric near space，is developed.Model validations are performed.The spectral element model with atmospheric near space resolved（SEMANS）is integrated for 10 years under configuration of81 local elements in each projected face of cubed spherewith spectral approximation of8?th degree Gauss?Lobatto?Legendre interpolation polynomial and 66 vertical layers with top layer pressure 4.5×10－6hPa.Through verification against ERA?Interim reanalysis dataset from ECMWF（European Center for Medium?Range Weather Forecasts）and COSPAR（Committee on Space Research）international reference atmosphere 1986，it is found that SEMANS reproduces features of 2 wave?lengths pattern along zonal circle in the northern hemisphere and 1 wave?length pattern along zonal circle in the southern hemisphere at30 hPa；SEMANS reproduces zones of low temperature at about 100 hPa，0.001 hPa and high temperature at about 1 hPa，above0.000 1 hPa；SEMANSalso reproducesmain features of zonalmean zonalwind in January and July below 0.001 hPa level. Key words： spectral element；near space；atmospheric model；numerical simulation

0 Introduction

The atmosphere is divided into a number of layers associated with different features in temperature structure.The structure is determined primarily by the absorption of solar radiation［1］. Actions from military and commerce have extended greatly their scope in atmosphere in the vertical. Tomeet the demand for exploration and application，scientists proposed the concept of near space. Near space is the region of Earth′s atmosphere that lies between 20 and 100 km above sea level，encompassing the stratosphere，mesosphere，and the lower thermosphere.It is above where airliners fly in the sky but below orbiting satellites in the space.Satellites in the space cost high with long deployment cycles and are difficult to replenish after damage.Airliners in the sky can be attacked easily，have a low ability for survival and are also difficult to replenish after damage.Flight vehicles in the near space have high efficiency?cost ratio，goodmobility，little difficulty in payload technique and are easy for replenishment and maintenance.Since the distance from flight vehicles in the near space to earth surface is only about 1/20 to 1/10 of that from low orbit satellites to earth surface，sensors on those near space flight vehicles can detect low power transmissions that cannot be caught by satellites and it is easier to implementhigh resolution earth observation.Near space craft provides a rare air carrier for military applications such as airborne early warning，reconnaissance and monitoring.Itwill be an inevitable trend for early warning and tracking of cruisemissile to use near space craft.In addition，fighter air craft can fly up to 20 km；strategic ballistic missile flight canreach height of dozens of km；space airplanes can fly with a speed of12～25 Mach in an altitude of 30～100 km.These actions in near space demand meteorological supportwhich cannot be supplied by traditional numericalmodels and it becomes urgent and necessary to develop atmosphericmodels with the near space resolved.Under the context of global environment change，the atmospheric near space is subject to long?term changes due to bothman?made and solar variability effects.Man?made changes in the ozone layer and in the concentrations of other trace gases are expected to causemajor changes in the temperature structure within these regions，which will be most noticeable at high altitudes.These changes influence the atmospheric circulation，including possible effects on tropospheric climate，and may influence space operations and radio?wave communications［1］.Since much remains to be learned about this region，models that include many of the important effects related to dynamics，chemistry，and energetics，are important tools for carrying out related researches.

In community of geophysical computation，three numerical methods，finite difference method （FDM），spectral method（SM），and finite element method（FEM）are most frequently used. FDM，which is implemented through approximating differential quotientwith difference quotient，is efficient regard to programming and computational cost.The weakness of FDM is its low speed of convergence.SM，which ismore efficient than FDM［2］，is a high order，p?type weighted residual method.Since the basis functions of SM are global，it is usually applied to regular geometrical problem with smooth solutions.FEM，which is flexible for geometry and converges slowly as the number of the elements increased，is a low order，h?type weighted residual method.Spectral elementmethod（SEM），proposed by Patera，combine advantages of Galerkin spectralmethodswith those of finite elementmethods［3］.In SEM physical domain is broken up into several elements，which is usually triangular or quadrilateral［4－5］.Within each elementa spectral representation based on high order interpolation polynomials is used.Convergence is either obtained by increasing the degree of polynomials or by increasing the number of elements.A shallow water equation case showed that SEM provides not only good spational resolution and geometrical flexibity，but also affordability for long?range simulations［6］.SEM had been used for atmospheric dynamic core to perform three?dimensional climate modeling for low and middle atmosphere［7］.There have been many works on stratosphere simulation with numericalmodels discretized by FDM and SM.And it has been found that results from low?top models（model top below the stratopause）have very weak stratospheric variability on daily and interannual time scales，which leads to strongly underestimated frequency ofmajor sudden stratospheric warming events［8－9］.There are few atmosphericmodelswith model top above 100 km.Of which，the whole atmosphere community atmosphere model （WACCM）［10］and the whole atmospheremodel（WAM）［11］are themost popular ones.WACCM is an extension，upward through the thermosphere，of the national Center for Atmospheric Research （NCAR）community atmosphere model（CAM），which is an atmospheric component of the community earth system model（CESM）.WAM is also an extension，upward through the thermosphere，of the National Oceanic and Atmospheric Administration（NOAA）global forecast system（GFS）model as partof the integrated dynamics in the Earth's atmosphere（IDEA）project. TheWACCM and WAM are discretized with finite volumemethod（FVM）and SM，respectively.

In this work，we develop a spectral element model with atmospheric near space resolved （SEMANS）.A brief description of themodel is given in Section 1.Design of numerical experiment and datasets used for model validation are given in Section 2.Analysis on simulation results and model validation are shown in Section 3.In the last section，conclusions and discussions are given.

1 M odel discretization and physical processes

1.1 M odel equations

With static approximation，the control equations of SEAMNS are composed of a set of equations，which can be derived from their common forms appeared in text books of atmospheric dynamic.We have equation ofmomentum

equation of continuity

equation of temperature

and equation of constituent related towater（water vapor，cloud water，rain water，icewater and so on）

With appropriate choices for A and B，s is identical toσat the lowestmodel layer and p at high model layers.s is a hybrid ofσand p formiddlemodel layers.

The diagnostic equation for the geopotential is

whereφsurfis surface geopotential.

Tomake the equations complete，two diagnostic equations forandω

are added.Eqs.（1）－（8）make up themodel equations.

1.2 Spatial discretization

Themodel equations are discretized with SEM in the horizontal.To solve the polar problem inspherical coordinate，the earth sphere surface is projected onto its inner contact cubic and numerical computations are implemented in the projected local coordinates on elements of six faces.A 2?step projection（denoted by P）is needed to transform from spherical coordinate to projected local coordinates on elements（Step 1：From spherical coordinate to projected coordinates on six faces；Step 2：From projected coordinates on six faces to local coordinates on elements of six faces）.

Fig.1 A cubed sphere（There are 81 elements on each projected face.）

The projected faces are divided into quadrilateral elements （Fig.1）.The quadrilateral element has been broadly used in such works of Taylor et al［12］，Giraldo et al［13］and Dennis et al［14］. Similar to that in Refs.［12－14］， Gauss?Lobatto?Legendre（GLL）gridding scheme is applied with SEMANS to facilitate the integration satisfying the Gaussian quadrature rule. Differences between these works lie in that，the configuration of the order of polynomial and the number of elements is different；themodel equations are different；the physics in thesemodels are different.Works of Refs.［12－14］deal with problems of single layer or multi?layers mainly within troposphere with simple physical process parameterization schemes adopted，whereas SEMANS can simulate the atmosphere up to the low thermosphere with complicated physical process parameterization schemes used.

The SEM is a kind of weighted residualmethod.Themodel equations are written in Galerkin form and then transformed to formula expressed in local coordinates on elements whose domain of definition is［－1：1；－1：1］.

The variables are interpolated as

Here，αaremodel variable such as u，v，T and qi，αijare expansion coefficients ofαand they vary with time，x and y are local coordinates on elements.The interpolation functions hiare Legendre cardinal functions.

In SEMANS，the spatial integrations are approximated by Gauss?Lobatto integration.Choosing Legendre cardinal functions as weight functions，the integration equations are discretized to equations of expansion coefficients of?α/?t．The solutions can be expressed with projection operator P as

where RHS means other items on the right hand side of equation.If we exchange the sequence of projection and temporal integration，the computation can be divided into 2 steps：

The first step is implemented within local elements and no data exchange between elements is needed.The computation in the second step needs data exchange between elements.

Themodel equations are discretized with FDM in the vertical.stop＝s1/2and sK＋1/2＝1．Thediscretization operator of?/?s isδs（X）k＝（Xk＋1/2－Xk－1/2）/（sk＋1/2－sk－1/2）and?/?s isδs（X）k＝（2πkΔsk）－1［（π）k＋1/2（Xk＋1－Xk）＋（π）k－1/2（Xk－Xk－1）］whereΔsk＝sk＋1/2－sk－1/2．

1.3 Physical processes

Zhang and McFarlane scheme［15］is used for deep convection.Park and Bretherton scheme［16］is used for shallow convection.A nonlocal diffusion scheme［17］is used for boundary layer parameterization.The near surface constant flux over land is calculated with Monin?Obukhov similarity theory and those over ocean water and sea ice are calculated with bulk formula.Below 60 km in altitude，solar radiation is computed withδ?Eddington approximation［18］and thermal radiation is computed with Ramanathan et al scheme［19］.Over 100 km，radiation heating rates from TIME?GCM（thermosphere ionospheremesosphere electrodynamics general circulationmodel）［20］are used. For height between 60 km and 100 km，weighted average of TIME?GCM results and radiation transfer output is used and weight coefficients are dependent onmodel layer.No chemical processes are considered at present.

1.4 Temporal integration scheme and boundary conditions

Leap?frog scheme is used for temperal integration.At boundaries of top and bottom，no flux conditions in which（π）1/2＝（π）K＋1/2＝0，are adopted.

2 Numerical experiments and datasets

Themodel atmosphere is divided into 66 layers with a pressure at the top 4.5×10－6hPa （Fig.2）.Each projected face is divided into 81 local elements and GLL interpolation polynomial of degree 8 is used to discretize the equation in each local element（Fig.3）.

Fig.2 Vertical coordinate and model layers

Fig.3 Gauss?Legendre?Lobatto interpolation polynomialswith degree 8

The data used for low boundary conditions andmodel initial values are all from the ERA?Interim dataset［21］.Dataset used for model validation are ERA?Interim dataset（from 1985 to 1994）and COSPAR（Committee on Space Research）International Reference Atmosphere1986（CIRA?86）［22］.

Initial values are from the atmospheric state at0 hour，0 minute，0 second of GMT（Greenwich Mean Time）on January 1st1989.Formodel layers above 1 hPa，values of variables such as zonal wind，meridional wind and temperature are identical to those of the highest model layer below 1hPa.Digital filter initialization［23］is used to reduce the initial imbalances and spurious gravity wave amplitude.The sea surface temperature，sea ice concentration，soil temperature and soilmoisture are all 12 monthly means of 1989，which are interpolated to model time needed within a year and used as bottom conditions.Land cover fraction and model terrain are interpolated from U.S. Geological Survey（USGS）dataset.

SEMANS has been integrated for 10 years and monthly mean model results are saved for analysis.During the integration，globalmean total energy（sum of sensible heat，potential energy，latent heat for phase change between cloud water and cloud ice，latent heat for phase change between vapor and cloud ice）and airmass（represented by surface air pressure）aremonitored every step.There is no obvious strange value found in theses variables.The globalmeanmonthly averaged surface air pressure is shown in Fig.4.It can be seen that the simulated seasonal variation of global mean airmass is reasonable and there is no obvious trend longer than a year.

Fig.4 Globalmean monthly averaged surface air pressure（Unit of vertical axis is Pa and unit of horizontal axis is year.）

3 Simulation results

Since atmospheric data of observation is scarce above altitude of 10 hPa level，simulation results are verified for stratosphere below altitude of 10 hPa level only.SEMANS can reproduce the features of 2 wave?lengths pattern along zonal circle in the northern hemisphere and 1 wave?length pattern along zonal circle in the southern hemisphere（Fig.5）.The simulated strength of disturbance in the northern hemisphere is weaker than that of ERA?Interim and the simulated strength of disturbence in the southern hemisphere is similar to that of ERA?Interim.Compared to result from ERA?Interim，the simulated extent of lower value area near the equator enclosed by 23 800 gpm is smaller.

SEMANS reproduced main features of zonal mean zonal wind in January.The positions of westerly belt，easterly belt and large wind speed centers in January simulated by themodel coincide with those of ERA?Interim well.The simulated intensity ofwesterly jet in northern hemisphere below 100 hPa level and above 30 hPa level isweaker than those of ERA?Interim.The simulated intensity of easterly near the equator is also weaker than those of ERA?Interim（Figs.6（a），6（b））.

SEMANS also reproduced main features of zonalmean zonal wind in July.The positions of westerly belt，easterly beltand largewind speed centers in July simulated by themodel coincidewiththose of ERA?Interim well.The simulated intensity ofwesterly jet in southern hemisphere above 30 hPa level is stronger than that of ERA?Interim.The span of simulated easterly belt below 100 hPa near the equator is smaller than that of ERA?Interim（Figs.6（c），6（d））.

Fig.5 Distribution of geopotential height at30 hPa from（a）SEMANSand（b）ERA?Interim for 10?year average （Contour interval is 200 gpm.）

Fig.6 Distribution of zonalmean zonal wind from SEMANS（left）and ERA?Interim（right）in January（upper）and July （lower）for 10?year average（Contour intervals are 5 m·s－1for（a），（b）and 10 m·s－1for（c），（d）.）

The lower part（0～120 km）of CIRA?86，consisting of themonthlymean values of temperature and zonalwind with almost global coverage（80°N～80°S），is used for evaluation of SEMANS.As illustrated in Fig.7，SEMANS reproduces the zones of low temperature atabout100 hPa，0.001 hPa and high temperature atabout1 hPa；the characteristic ofhigh temperature above0.0001 hPa is alsoreproduced.But the simulated temperature at about 0.001 hPa is higher than that of CIRA?86 in January（Figs.7（b），7（a））.The high temperature center of 1 hPa level at North Pole is not reproduced by SEMANS in July（Figs.7（d），7（c））.

Fig.7 Distribution of zonalmean temperature from CIRA?86（left）and SEMANS（right）in January（upper）and July（lower）（Unit of vertical axis is hPa and contour interval is 20 K.）

Compared to CIRA?86，SEMANS reproduces the main features of distribution of zonal mean zonal wind below 0.001 hPa level，whereas there are much more discrepancies above that level （Figs.8（b），8（a）and 8（d），8（c））.These discrepancies are at least partly due to effect of exclusion of gravity wave in SEMANS.

4 Conclusion and discussion

A spectral elementmodel with atmospheric near space resolved was developed and preliminary work onmodel validation was done.Due to scarcity of observation data over10 hPa level，simulation results are verified against ERA?Interim reanalysis dataset for atmospheric stratosphere below 10 hPa level and reference atmosphere CIRA?86 for high levels.SEMANS reproduced the features of 2 wave?lengths pattern along zonal circle in the northern hemisphere and 1 wave?length pattern along zonal circle in the southern hemisphere at 30 hPa.SEMANS also reproduced themain features of zonalmean zonalwind in January and July.There areminute discrepancies between model results and ERA?Interim reanalysis dataset.Since there are differences even between reanalysis datasets from different sources，the discrepancies found here is thought to be acceptable.SEMANS reproduces the zones of low temperature at about100 hPa，0.001 hPa and high temperature at about 1 hPa and above 0.0001 hPa.But there is a spurious high temperature band around 0.2 hPa level.SEMANS reproduces the main features of distribution of zonalmean zonal wind below 0.001 hPa level，whereas there aremuch more discrepancies above that level.

Fig.8 Distribution of zonalmean zonalwind from CIRA?86（left）and SEMANS（right）in January（upper）and July（lower）（Unit of vertical axis is hPa and contour interval is 20 m·s－1.）

The chemical processes and gravity wave drag effects are not considered in SEMANSatpresent. To improve the performance of SEMANS，morework such as consideration of chemical processes and gravity wave drag effects and fine tuning ofmodel parameters needs to be done.

The advantage of SEM over FEM and FDM is in the exponential convergence.In practice，to improvemodel resolution，one can either keep the order of polynomial（denoted by N）fixed and increase the number of elements（denoted by K），or keep K fixed and increase N.In addition，one can reposition elements，keeping both N and K fixed，or change N in different elements.In simulation with spectralmodel，such cases as negative terrain height and negative specific humidity may occur due to discontinuity.This is ameliorated but can notbe avoided in SEM.In recent years，numerical model with FDM attracts many researchers'attention again due to advantages in implementation of parallelization and computation property（The amount of computation doesn't increase dramatically with enhancement of resolution）.

Acknowledgments：The work is sponsored by the Natural Science Foundation of China（Grant No.41276190）.The ERA?Interim reanalysis datasets are downloaded from the ECMWF data server website.We thank the two anonymous reviews for their suggestions.

［1］ Geller M A，Brasseur G P，Evans JV，etal.Upper?atmosphere and near?earth space research entering the twenty?first century［M］∥Board on Atmospheric Sciences and Climate Commission on Geosciences，Environment，and Resources，National Research Council.The atmospheric sciences entering the twenty?first century.Washington DC，USA：National Academy Press，1998：199－271.

［2］ Haidvogel D B，Robinson A R，Schulman E E.The accuracy，efficiency，and stability of three numerical models with application to open ocean problems［J］.JComp Phys，1980，34：1－53.

［3］ Patera A T.A spectral elementmethod for fluid dynamics：Laminar flow in a channel expansion［J］.J Comput Phys，1984，54：468－488.

［4］ Liu Xiying.Introduction of physical processes into a dynamic framework of spectral element and numerical experiment of medium range weather prediction［C］∥Zhou Liandi.Proceedings of the 20th National Seminar of Hydrodynamics，Beijing，China：Ocean Press，2007：224－230.

［5］ Liu Xiying.Numerical simulation of typhoon movement with a regional spectral element barotropic atmospheric model［J］.Chinese Journal of Computational Physics，2011，28：35－40.

［6］ Ma Hong.A spectral element basin model for the shallow water equations［J］.JComp Phys，1993，109：133－149.

［7］ Thomas S，Loft R.The NCAR spectral element climate dynamical core：Semi?implicit Eulerian formulation ［J］.Journal of Scientific Computing，2005，25：307－322.

［8］ Charlton?Perez A J，et al.On the lack of stratospheric dynamical variability in low?top versions of the CMIP5 models［J］.JGeophys Res Atmos，2013，118：2494－2505.

［9］ Austin J，Horowitz LW，Daniel Schwarzkopf M，etal.Stratospheric ozone and temperature simulated from the preindustrial era to the present day［J］.JClimate，2013，26：3528－3543.

［10］ Brakebusch M，Randall C E，Kinnison D E，et al.Evaluation of whole atmosphere community climate model simulations of ozone during Arctic winter 2004－2005［J］.JGeophys Res Atmos，2013，118：2673 －2688.

［11］ Wang H，F(xiàn)uller?Rowell T J，Akmaev R A，et al.First simulations with a whole atmosphere data assimilation and forecast system：The January 2009 major sudden stratospheric warming［J］.JGeophys Res，2011，116：A12321.

［12］ Taylor M，Tribbia J，Iskandarani M.The spectral elementmethod for the shallow water equations on the sphere［J］.JComput Phys，1997，130：92－108.

［13］ Giraldo F X，Rosmond T E.A scalable spectral element Eulerian atmosphericmodel（SEE?AM）for NWP：Dynamical core tests［J］.Mon Wea Rev，2004，132：133－153.

［14］ Dennis J，et al.CAM?SE：A scalable spectral element dynamical core for the Community Atmosphere Model［J］.Int JHigh Perform Comput Appl，2012，26：74－89.

［15］ Zhang G J，McFarlane N A.Sensitivity of climate simulations to the parameterization of cumulus convection in the Canadian Climate Centre general circulation model［J］.Atmosphere?Ocean，1995，33：407－446. ［16］ Park S，Bretherton C S.The university ofwashington shallow convection and moist turbulence schemes and their impact on climate simulationswith the community atmospheremodel［J］.JClimate，2009，22：3449 －3469.

［17］ Holtslag A A M，Boville B A.Local versus nonlocal boundary?layer diffusion in a global climate model ［J］.JClimate，1993，6：1825－1842.

［18］ Briegleb B P.Delta?Eddington approximation for solar radiation in the NCAR Community Climate Model ［J］.JGeophys Res，1992，97：7603－7612.

［19］ Ramanathan V，Downey P.A nonisothermal emissivity and absorptivity formulation for water vapor［J］.J Geophys Res，1986，91：8649－8666.

［20］ Marsh D，Roble R.TIME?GCM simulations of lower?thermospheric nitric oxide seen by the halogen occultation experiment［J］.JAtmos Solar Terr Phys，2002，64：889－895.

［21］ Dee D P，Uppala SM，Simmons A J，etal.The ERA?Interim reanalysis：Configuration and performance of the data assimilation system［J］.Quart JR Meteorol Soc，2011，137：553－59.

［22］ Fleming E L，Chandra S，Shoeberl M R，et al.Monthly mean global climatology of temperature，wind，geopotential height and pressure for 0～120 km［R］.Technical Memorandum 100697，National Aeronautics and Space Administration，Washington D.C.，1988：1－56.

［23］ Chen M，Huang X Y.Digital filter initialization for MM5［J］.Mon Wea Rev，2006，134：1222－1236.

A Post?processing M ethod for Smoothed Particle Hydrodynam ics Data

LIFupeng1，2，WANG Jiwen1，2，OU Mang2
（1.Key Laboratory of Intelligent Computing＆Signal Processing，Anhui University，Hefei 230039，China；
2.Department ofComputer Science and Engineering，Anhui University，Hefei 230039，China）

A post processing method for smoothed particle hydrodynamics data is presented on basis of Delaunay triangulation. Pairing particle set in scope of each particle in flow field is confirmed，and then discretized by Delaunay triangulation on condition of present pairing particle set，which leads to a set of triangle elements with particles as nodes.At last non material units are filtered. Examples indicate that SPH resultwith complex boundary conditions and severe deformation of fluid region can be treated.

smoothed particlemethod；non convex region；Delaunay triangulation；post processing

A Spectral Element M odel w ith Atmospheric Near Space Resolved（SEMANS）

LIU Xiying，LIU Chunyan，YAO Shanshan

（College ofMeteorology and Oceanography，PLAUniversity ofScience and Technology，Nanjing 211101，China）

P435 Document code：A

O302

A

2014－01－24；

2014－06－18

安徽高校省級(jí)自然科學(xué)研究（KJ2013A009）資助項(xiàng)目

李付鵬（1974－），男，安徽阜陽，博士生，從事計(jì)算流體力學(xué)、計(jì)算機(jī)流體動(dòng)畫模擬研究，E?mail：lifupenganda＠163.com

Received date： 2014－01－24；Revised date： 2014－06－18

Received date：2013－10－21；Revised date：2014－08－29

Foundation item s：Supported by National Natural Science Foundation of China（41276190）

Biography：Liu Xiying（1970－），male，PhD，associate professor，major in atmospheric dynamics and numerical simulation，E?mail：liuxynj＠gmail.com