王圣業(yè) 王光學(xué) 董義道 鄧小剛

1)(國防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院,長沙 410073)2)(中山大學(xué)物理學(xué)院,廣州 510275)

基于雷諾應(yīng)力模型的高精度分離渦模擬方法?

王圣業(yè)1)王光學(xué)2)董義道1)鄧小剛1)?

1)(國防科學(xué)技術(shù)大學(xué)航天科學(xué)與工程學(xué)院,長沙 410073)2)(中山大學(xué)物理學(xué)院,廣州 510275)

湍流流動(dòng),雷諾應(yīng)力模型,分離渦模擬,加權(quán)緊致非線性格式

1 引 言

準(zhǔn)確預(yù)測翼型和機(jī)翼在接近甚至超過失速迎角時(shí)的氣動(dòng)特性,已成為現(xiàn)代飛行器設(shè)計(jì)的重要方面.然而,對于高雷諾數(shù)大迎角條件下的分離流動(dòng),精確的計(jì)算流體力學(xué)(computational fluid dynamics,CFD)模擬仍十分困難[1].

任何流動(dòng)模擬的結(jié)果都依賴于所代表的物理模型的準(zhǔn)度以及求解相關(guān)方程的數(shù)值方法的精度[2].為了精確預(yù)測流動(dòng)物理,主要是湍流,最好是直接求解所有尺度的湍流脈動(dòng)甚至是最小尺度,即直接數(shù)值模擬;或求解大尺度湍流而模型化最小尺度,即大渦數(shù)值模擬(large eddy simulation,LES).然而,在真實(shí)飛行雷諾數(shù)下,即使附著邊界層內(nèi)的最大尺度脈動(dòng)也會(huì)變得很小,這將使計(jì)算花費(fèi)巨增而難以承擔(dān).因此,基于雷諾平均納維-斯托克斯方程(Reynolds average Navier-Stokes,RANS)的湍流模型方法仍然是工業(yè)應(yīng)用CFD的中堅(jiān)力量[3].

混合雷諾平均/大渦模擬(hybrid RANS/LES)方法結(jié)合了大渦模擬方法在分離流動(dòng)區(qū)域的高分辨率與雷諾平均方法在附著邊界層中的高效率,在大迎角分離模擬中受到了廣泛的關(guān)注[4?6].而在hybrid RANS/LES方法中,Spalart等[7]于1997年提出的分離渦模擬(detached eddy simulation,DES)方法由于其簡單、高效等特點(diǎn)得到了廣泛的發(fā)展和應(yīng)用[8].然而,最初版本的DES方法基于一方程Spalart-Allmaras(SA)模型,在大規(guī)模分離流動(dòng)支配的算例中表現(xiàn)良好,但對于由逆壓梯度造成的初始尾緣分離等問題卻很難準(zhǔn)確估計(jì).因此,為了更準(zhǔn)確地模擬翼型和機(jī)翼的失速分離問題,尤其在失速迎角附近,DES的RANS部分需要采用更先進(jìn)的湍流模型[9].2001年,Strelets[10]提出了基于二方程剪切壓力傳輸(shear stress transport,SST)模型的DES方法,對翼型最大升力迎角附近時(shí)的分離問題有一定的提高.然而,由于線性渦粘模型自身的限制,其效果仍不理想.2006年,Greschner等[11,12]發(fā)展了基于EASM LLk-ε模型(一類非線性渦粘模型)的DES方法,并在串聯(lián)圓柱-翼型算例中與傳統(tǒng)k-ε-DES方法進(jìn)行了對比.結(jié)果表明,EASM-DES方法能夠更精細(xì)地分辨脫落渦結(jié)構(gòu),并且預(yù)測的氣動(dòng)力更接近試驗(yàn)值.但由于EASM模型仍基于Boussinesq假設(shè)框架,其對雷諾應(yīng)力的求解仍然不直接,因此模擬效果對算例的敏感性較大.

雷諾應(yīng)力模型(Reynolds stress model,RSM)是更高階矩封閉的湍流模型,采用六個(gè)方程求解雷諾應(yīng)力項(xiàng)和一個(gè)方程求解尺度約束變量(例如湍動(dòng)能耗散率ε或比耗散率ω),能夠直接估計(jì)雷諾應(yīng)力各向異性、流線彎曲效應(yīng)等[13].2011年,Probst等[9]基于ε型RSM模型,發(fā)展了εh-DES類方法,并在HGR-01翼型上取得了很好的結(jié)果.但由于ε型RSM模型魯棒性較差,其適應(yīng)性將不如ω型RSM模型[14].

近年來,德國宇航局的Eisfeld等[13,15,16]發(fā)展了一種魯棒的RSM模型,Speziale-Sarkar-Gatski/Launder-Reece-Rodi(SSG/LRR)-ω模型,并在一系列復(fù)雜的航空分離流動(dòng)中證明了其優(yōu)勢.美國航空航天局的Rumsey等[2]在湍流模型資源(turbulence modeling resource,TMR)網(wǎng)站上將SSG/LRR-ω模型設(shè)為最推薦的RSM,并開展了廣泛的驗(yàn)證與確認(rèn)工作.國內(nèi),董義道等[17]基于TMR網(wǎng)站,成功開展了SSG/LRR-ω模型的初步應(yīng)用研究,但在應(yīng)用過程中發(fā)現(xiàn)基于RANS方法的SSG/LRR-ω模型對翼型大迎角狀態(tài)的模擬仍不理想.本文基于上述工作,在SSG/LRR-ω模型上發(fā)展了一類DES方法,并在典型航空算例(17?,45?及60?迎角的NACA0012翼型;12?迎角的NACA4412翼型;24.6?迎角的鈍前緣三角翼)中進(jìn)行了驗(yàn)證.為了進(jìn)行對比研究,本文還同時(shí)采用了傳統(tǒng)基于線性渦粘模型的分離渦模擬方法,SST-DES方法.

對于DES方法,在RANS部分采用更先進(jìn)的湍流模型,而在分離區(qū)LES部分則對數(shù)值精度的要求更高.相比傳統(tǒng)二階方法,高階方法在采用較少計(jì)算花費(fèi)獲得更低誤差方面具有更大的潛力[18?20].由于在計(jì)算效率和計(jì)算精度上的優(yōu)勢,高階方法受到了廣泛關(guān)注,并被推薦用于LES或hybrid RANS/LES方法中[21].本文中對RANS方程和模型方程的離散均采用了高階精度的加權(quán)緊致非線性格式(weighted compact nonlinear scheme,WCNS)[22,23],并且在計(jì)算網(wǎng)格導(dǎo)數(shù)時(shí)采用了滿足自由流保持的對稱網(wǎng)格導(dǎo)數(shù)算法(symmetrical conservative metric method,SCMM)[24].

2 湍流模型

2.1 SSG/LRR-ω模型

對Navier-Stokes方程進(jìn)行Favre平均后出現(xiàn)雷諾應(yīng)力項(xiàng).傳統(tǒng)線性渦粘模型,如SA模型和SST模型,對雷諾應(yīng)力項(xiàng)的求解均基于Boussinesq假設(shè),而RSM則直接對其求解.

其中t,xk(k=1,2,3)分別為時(shí)間和空間坐標(biāo)分量;,分別表示Favre平均的密度、速度分量和雷諾應(yīng)力張量,而全雷諾應(yīng)力不加特殊說明,下面出現(xiàn)的變量均為Favre平均后的量.(1)式右端第一項(xiàng)為生成項(xiàng),可準(zhǔn)確求解:

右端第三項(xiàng)為耗散項(xiàng),

其中ε表示各向同性耗散率,需要通過額外的湍流尺度方程求解得到,δij為Kronecker符號.右端第四項(xiàng)為輸運(yùn)項(xiàng),

其中各向異性張量

表1 SSG/LRR-ω模型封閉系數(shù)中SSG和LRR部分的貢獻(xiàn)Table 1.Values of closure coefficients for the SSG and the LRR contributions to the SSG/LRR-ω redistribution term.

其他各項(xiàng)定義如下:

為了封閉上述模型,需要額外引入尺度方程,SSG/LRR-ω模型借鑒了Menter的思路,對ε方程和ω方程進(jìn)行了混合[15]:

方程(9)的變量為比耗散率ω.各向同性耗散率通過下式計(jì)算:

混合函數(shù)定義如下:

其中,dw為距壁面法向距離,其他系數(shù)具體見表2.

表2 ω方程系數(shù)中ε和ω部分的貢獻(xiàn)Table 2.Values of coefficients of ω-equation corresponding to the ε and ω parts.

2.2 SSG/LRR-DES類方法

DES方法由Spalart等[7]提出,將LES方法和SA模型結(jié)合,并通過比較當(dāng)?shù)鼐W(wǎng)格尺度與壁面距離實(shí)現(xiàn)兩種方法的自動(dòng)切換.其后,Strelet[10]借鑒Spalart的思想,通過比較當(dāng)?shù)鼐W(wǎng)格尺度與湍流長度尺度,將DES方法引入SST模型.本文采用的SSG/LRR-ω模型與SST模型均為基于比耗散率ω的湍流模型,因此在構(gòu)造SSG/LRR-DES方法時(shí)將主要借鑒SST-DES的構(gòu)造思路.

首先定義網(wǎng)格長度尺度

和湍流長度尺度

然后得到SSG/LRR-DES的限制器

對雷諾應(yīng)力方程耗散項(xiàng)中的各向同性耗散率(10)式進(jìn)行修正:

這里CDES為DES常數(shù),通過Menter的混合函數(shù)(12)式得到

原始DES方法的限制器在復(fù)雜網(wǎng)格上處理RANS和LES的轉(zhuǎn)換過程中過于生硬,會(huì)造成?；瘧?yīng)力損耗(modeled stress depletion,MSD)等問題[8].2006年,Spalart等[26]借鑒Menter的SST模型構(gòu)造思想,采用“延遲LES函數(shù)”改善了原始版本的MSD問題,發(fā)展出了延遲分離渦模擬(delayed DES,DDES)方法.2008年,Shur等[27]將DDES方法和壁面模型大渦模擬(wall-modeled LES,WMLES)方法結(jié)合,發(fā)展了強(qiáng)化延遲分離渦模擬(improved DDES,IDDES)方法,并在非定常分離流動(dòng)中得到了廣泛的應(yīng)用.本文的計(jì)算中均采用IDDES形式,其構(gòu)造過程如下.

首先采用WMLES中網(wǎng)格尺度的加權(quán)思想,定義新的IDDES限制器:

其中,fd為原DDES中的轉(zhuǎn)換函數(shù),fd=1?tanh(8rd)3,fB為原WMLES中的轉(zhuǎn)換函數(shù),fB=min(2exp(?9α2),1),這里α=0.25?dw/hmax;fe為壁面模擬的控制函數(shù),其作用是保證在壁面附近網(wǎng)格分辨率滿足LES要求時(shí),忽略過渡區(qū)雷諾應(yīng)力的影響.

其中fe1定義為

fe2定義為

(22)式中rd,rdt和rdl均與文獻(xiàn)[27]中相同.

另外,在計(jì)算當(dāng)?shù)鼐W(wǎng)格尺度時(shí),Shur等[27]還引進(jìn)了壁面距離的影響,新的Lg為

其中hmax=max(?x,?y,?z),hwn為壁面垂直方向的網(wǎng)格步長.IDDES中的系數(shù)為Cw=0.15,ct=1.87和cl=5.0.

3 數(shù)值方法

3.1 高精度數(shù)值方法

WCNS格式由Deng等[22]在2000年提出.其后,不同學(xué)者[20,28,29]發(fā)展了多種形式的WCNS格式,并在廣泛的流動(dòng)問題中證明了該格式的優(yōu)勢.本文采用的是WCNS系列格式中一種典型的五階顯式離散格式WCNS-E-5[23].WCNS-E-5格式由于其低耗散、高魯棒和優(yōu)秀的自由流與渦保持特性,被廣泛應(yīng)用于各種實(shí)際流動(dòng)問題的高精度數(shù)值模擬中.其中岡敦殿等[30]將WCNS-E-5格式應(yīng)用于平板圓臺(tái)突起物繞流的LES中,并和采用平面激光散射技術(shù)的試驗(yàn)結(jié)果進(jìn)行了對比,證明了格式精細(xì)模擬湍流流動(dòng)的可行性.

另外本文在計(jì)算網(wǎng)格導(dǎo)數(shù)及雅克比時(shí),采用了滿足幾何守恒律的對稱守恒網(wǎng)格導(dǎo)數(shù)算法[24],有利于提高高精度有限差分方法的魯棒性并減小數(shù)值誤差.本文中時(shí)間推進(jìn)均采用子迭代步基于LU-SGS(lower-upper symmetric-Gauss-Seidal)方法的雙時(shí)間步法.

3.2 對雷諾應(yīng)力方程的高精度離散

較差的數(shù)值穩(wěn)定性是限制RSM使用的障礙,尤其是結(jié)合高精度數(shù)值方法時(shí).在迭代的過程中,雷諾應(yīng)力項(xiàng)很可能不滿足Schumann[31]提出的現(xiàn)實(shí)性限制.Chassaing等[32]提出一種魯棒的隱式格式,依賴于使用當(dāng)?shù)仉p時(shí)間步技術(shù)和顯式增強(qiáng)現(xiàn)實(shí)性限制.Yossef[33]發(fā)展了一種無條件正收斂隱式格式和一種現(xiàn)實(shí)性限制,保證了任意時(shí)間步內(nèi)雷諾正應(yīng)力項(xiàng)為正值.本文中,在雷諾應(yīng)力方程的迭代過程中添加如下限制條件:

除了時(shí)間迭代方面,也應(yīng)同時(shí)關(guān)注雷諾應(yīng)力方程的空間離散.對于加權(quán)型有限差分格式,非線性權(quán)是影響數(shù)值穩(wěn)定性的主要方面.需要指出的是,在對雷諾應(yīng)力方程采用WCNS類格式離散時(shí),加權(quán)公式中的小量εIS應(yīng)取為10?18,而非文獻(xiàn)[23]中對Navier-Stokes方程離散時(shí)的10?6.

4 計(jì)算結(jié)果

4.1 NACA0012翼型大迎角氣動(dòng)特性模擬

NACA0012翼型大迎角分離算例是推動(dòng)DES發(fā)展的經(jīng)典算例.2001年,Strelets等[7]在提出SST-DES類方法時(shí),就在該算例上與原始的SADES類方法進(jìn)行了對比,但發(fā)現(xiàn)并無本質(zhì)的差異.2016年,Yang和Zha[34]結(jié)合高精度WENO(weighted essentially nonoscillatory)格式,采用改進(jìn)的SA-IDDES方法對該算例進(jìn)行了模擬,但發(fā)現(xiàn)在中等迎角下SA-IDDES方法還沒有SA-URANS方法準(zhǔn)確.本文為了更好地對SSG/LRR-IDDES方法進(jìn)行對比研究,計(jì)算條件和網(wǎng)格均參考Yang和Zha的工作[34],并同時(shí)采用了SST-IDDES方法.

圖1 NACA0012翼型三維粗糙網(wǎng)格Fig.1.Three-dimensional coarse grid of the NACA0012 airfoil.

NACA0012翼型基于弦長c的雷諾數(shù)為1.3×106,基于自由流速度的馬赫數(shù)為0.5.在0?—90?迎角范圍內(nèi),選取17?(大范圍邊界層分離),45?和60?(大規(guī)模脫體分離)三個(gè)典型狀態(tài)進(jìn)行模擬.由于該算例在高雷諾數(shù)(Re>1.0×105)下受雷諾數(shù)影響很小,因此升力和阻力系數(shù)的參考值選取Re=2.0×106的試驗(yàn)值[34,35].本文計(jì)算采用粗和密兩套O型網(wǎng)格,網(wǎng)格單元數(shù)分別為192×102×30和288×102×30.第一層網(wǎng)格小于1;遠(yuǎn)場長度約為100c;展向長度為1.0c.圖1為NACA0012翼型三維網(wǎng)格示意圖.時(shí)間推進(jìn)采用雙時(shí)間步方法,周期為T=c/U∞,時(shí)間步長0.01T.時(shí)間平均統(tǒng)計(jì)從50T時(shí)開始,持續(xù)50T.

圖2展示了3個(gè)典型迎角下統(tǒng)計(jì)時(shí)間的平均升、阻力系數(shù).在17?迎角時(shí),翼型處于失速狀態(tài),所有方法得到的阻力系數(shù)相近且與試驗(yàn)吻合較好;但對于升力系數(shù),SSG/LRR-IDDES在粗網(wǎng)格上略低于試驗(yàn)值,而在密網(wǎng)格上與試驗(yàn)值吻合很好;而其他方法均明顯高于試驗(yàn)值.Yang和Zha[34]在17?迎角的計(jì)算中,得出了WENO+SA-IDDES結(jié)果明顯差于WENO+SA-URANS的現(xiàn)象,但并未給出解釋.而本文采用的線性渦粘模型、SST模型并未出現(xiàn)該現(xiàn)象.在45?和60?迎角時(shí),翼型處于過失速狀態(tài),所有IDDES方法的結(jié)果均明顯優(yōu)于URANS方法.同時(shí)對比兩類URANS方法,SSG/LRR-ω模型略優(yōu)于SST模型,而該現(xiàn)象與相關(guān)湍流模型在分離渦處對雷諾應(yīng)力預(yù)測的能力有關(guān).

圖3展示了17?迎角時(shí)升、阻力系數(shù)在100T內(nèi)的變化過程,其中AoA表示迎角.可以看到在25T后,SST-URANS和SSG/LRR-URANS得到的升、阻力系數(shù)出現(xiàn)類似簡諧振蕩的發(fā)展過程.而SSTIDDES和SSG/LRR-IDDES方法得到的升、阻力系數(shù)均為無周期振蕩,表明其均能夠描述分離湍流的隨機(jī)性.對比試驗(yàn)值,SST-IDDES得到的升力系數(shù)明顯偏高,并且隨網(wǎng)格加密無明顯改善,反映了基于線性渦粘模型的DES方法在翼型最大升力迎角附近模擬的局限.

經(jīng)典的渦拉伸原理[36]表明,三維性是湍流最本質(zhì)的特性之一.圖4給出的粗網(wǎng)格上100T時(shí)刻Q判據(jù)為0的三維等值面圖中,SST-URANS預(yù)測的渦脫落過程明顯為二維過程;SSG/LRRURANS能夠預(yù)測出上翼面附近較小尺度的渦,但整個(gè)渦脫落過程仍未表現(xiàn)出明顯的三維性;而SST-IDDES和SSG/LRR-IDDES均得到了高度混亂的三維渦結(jié)構(gòu),但SSG/LRR-IDDES對前緣渦拉伸-彎曲-破裂的過程則模擬的更為精細(xì).

圖5對比了不同方法在粗網(wǎng)格上100T時(shí)刻0.5展長處的展向速度,其中SST-URANS得到的展向流動(dòng)十分微弱,也印證了其為二維渦脫落過程;SSG/LRR-URANS得到的展向流動(dòng)稍強(qiáng)于SSTURANS,但也未表現(xiàn)出明顯的三維性;而SSTIDDES和SSG/LRR-IDDES均預(yù)測出了明顯的展向流動(dòng).

圖2 (網(wǎng)刊彩色)0?—90?范圍內(nèi)3個(gè)典型迎角下的升阻力結(jié)果對比Fig.2.(color online)Comparisons of lift and drag coefficients at 3 typical attack angles of among 0?–90?.

圖3 (網(wǎng)刊彩色)17?迎角時(shí)升、阻力系數(shù)變化過程Fig.3.(color online)Lift and drag coe ff cient history at attack angle of 17?.

圖4 (網(wǎng)刊彩色)17?迎角下100T時(shí)刻Q判據(jù)為0的等值面圖,顏色由馬赫數(shù)標(biāo)識(shí)Fig.4.(color online)Iso-surface of the Q-criterion=0 at 100T and attack angle of 17?,colored by Mach number.

圖5 (網(wǎng)刊彩色)17?迎角下100T時(shí)刻,0.5展長處展向速度云圖Fig.5.(color online)Spanwise velocity contours of 50%span at 100T and attack angle of 17?.

圖6展示了45?迎角時(shí)升、阻力系數(shù)在100T內(nèi)的變化過程.可以看到在25T后,SST-URANS得到的升、阻力系數(shù)出現(xiàn)類似簡諧振蕩的發(fā)展過程,而SSG/LRR-URANS并未得到有規(guī)律的簡諧振蕩過程,這與RSM能夠估計(jì)雷諾應(yīng)力的各向異性有關(guān).而SST-IDDES和SSG/LRR-IDDES方法得到的升、阻力系數(shù)均在試驗(yàn)值附近無規(guī)律振蕩.

圖6 (網(wǎng)刊彩色)45?迎角時(shí)升、阻力系數(shù)變化過程Fig.6.(color online)Lift and drag coe ff cient history at attack angle of 45?.

圖7 (網(wǎng)刊彩色)45?迎角下100T時(shí)刻Q判據(jù)為0的等值面圖,顏色由馬赫數(shù)標(biāo)識(shí)Fig.7.(color online)Iso-surface of the Q-criterion=0 at 100T and attack angle of 45?,colored by Mach number.

來流流過機(jī)翼后形成脫體的分離湍流,并將產(chǎn)生明顯的展向流動(dòng).圖7給出了粗網(wǎng)格上45?迎角,100T時(shí)刻,Q判據(jù)為0的三維等值面圖.SSTURANS僅預(yù)測出了大尺度的展向渦,其渦脫落過程仍為二維過程;SSG/LRR-URANS能夠預(yù)測出少量的流向和法向渦,但仍未表現(xiàn)出明顯的三維過程;而SST-IDDES和SSG/LRR-IDDES得到了高度混亂的大小渦結(jié)構(gòu),預(yù)測出了明顯的三維效應(yīng). 圖8對比了不同方法在粗網(wǎng)格上100T時(shí)刻0.5展長處的展向速度,更加清晰地表明各方法對于分離湍流的預(yù)測能力. 這也是圖3中SST-IDDES和SSG/LRR-IDDES的結(jié)果優(yōu)于SSG/LRR-URANS、更優(yōu)于SST-URANS的原因.

圖8 (網(wǎng)刊彩色)45?迎角下100T時(shí)刻,0.5展長處展向速度云圖Fig.8.(color online)Spanwise velocity contours of 50%span at 100T and attack angle of 45?.

圖9 (網(wǎng)刊彩色)60?迎角時(shí)升、阻力系數(shù)收斂過程Fig.9.(color online)Lift and drag coe ff cient history at attack angle of 60?.

圖10 (網(wǎng)刊彩色)60?迎角下100T時(shí)刻Q判據(jù)為0的等值面圖,顏色由馬赫數(shù)標(biāo)識(shí)Fig.10.(color online)Iso-surface of the Q-criterion=0 at 100T and attack angle of 60?,colored by Mach number.

圖11 (網(wǎng)刊彩色)60?迎角下100T時(shí)刻,0.5展長處展向速度云圖Fig.11.(color online)Spanwise velocity contours of 50%span at 100T and attack angle of 60?.

圖9展示了60?迎角時(shí)升、阻力系數(shù)的變化過程.與45?迎角時(shí)的結(jié)果類似,在25T后,SSTURANS得到的升、阻力系數(shù)出現(xiàn)明顯的周期性;而SSG/LRR-URANS，SST-IDDES和SSG/LRRIDDES呈現(xiàn)出無規(guī)律振蕩,但SSG/LRR-URANS得到的升、阻力系數(shù)平均值明顯高于試驗(yàn)值.圖10給出了60?迎角100T時(shí)刻下,Q判據(jù)為0的三維等值面圖.SST-URANS預(yù)測的渦脫落過程仍為二維過程;SSG/LRR-URANS得到的結(jié)果也仍未表現(xiàn)出明顯的三維過程;而SST-IDDES和SSG/LRRIDDES得到了高度混亂的大小渦結(jié)構(gòu),預(yù)測出了明顯的三維效應(yīng).圖11給出了60?迎角100T時(shí)刻0.5展長處的展向速度云圖,結(jié)果也與45?迎角時(shí)類似.

4.2 NACA4412翼型尾緣分離模擬

NACA4412翼型是經(jīng)典的低速湍流驗(yàn)證算例[37],有多種條件下的詳細(xì)試驗(yàn)數(shù)據(jù)進(jìn)行對比.本文選擇的是以Wadcock[38]實(shí)施的最大升力構(gòu)型(12?迎角)試驗(yàn)為參考的算例.試驗(yàn)條件為:翼型弦長0.9 m,雷諾數(shù)1.64×106,自由流速度29.1 m/s以及馬赫數(shù)0.085.該條件下翼型尾緣開始出現(xiàn)分離,對任何hybrid RANS/LES方法都是巨大的挑戰(zhàn)[39].

本文計(jì)算的條件設(shè)置與試驗(yàn)相同.采用粗和密兩套O型網(wǎng)格,展向拉伸0.1 m,網(wǎng)格單元數(shù)分別為169×123×30和249×123×30.第一層網(wǎng)格均小于1.圖12為NACA4412翼型三維網(wǎng)格示意圖.時(shí)間推進(jìn)采用雙時(shí)間步法,周期為T=c/U∞,時(shí)間步長0.01T.時(shí)間平均統(tǒng)計(jì)從40T時(shí)開始,持續(xù)40T.

表3列出了平均升力系數(shù)和分離位置的計(jì)算值與試驗(yàn)值.可以看到,SST-URANS、SSG/LRRURANS和SST-IDDES得到的結(jié)果均與試驗(yàn)值偏差較大,并且隨著網(wǎng)格加密未有改善. 而SSG/LRR-IDDES的結(jié)果與試驗(yàn)吻合較好.

圖12 NACA4412翼型三維粗網(wǎng)格Fig.12.Three-dimensional coarse grid of the NACA4412 airfoil.

表3 NACA4412翼型平均升力系數(shù)和分離位置的計(jì)算值與試驗(yàn)值對比Table 3.Comparation of computional results of lift coefficient and location of separation for NACA4412 airfoil case.

圖13 (網(wǎng)刊彩色)NACA4412翼型80T時(shí)刻,Q判據(jù)為0等值面圖,顏色由馬赫數(shù)標(biāo)識(shí)Fig.13.(color online)Iso-surface of the Q-criterion=0 at 80T for NACA4412 airfoil case,colored by Mach number.

圖13給出了粗網(wǎng)格上80T時(shí)刻Q判據(jù)為0的等值面圖.在翼型尾緣處,流動(dòng)開始分離并產(chǎn)生湍流尾跡渦.SSG/LRR-IDDES成功模擬了該過程,并且分離位置與試驗(yàn)吻合.圖14給出了該網(wǎng)格上得到的平均壓力系數(shù)分布和試驗(yàn)值的對比.SSG/LRR-IDDES在前緣得到的吸力值明顯低于另三種方法,并且與試驗(yàn)值更接近,這也是其得到更準(zhǔn)確升力系數(shù)的原因.在翼型尾緣處,由于SSG/LRR-IDDES成功模擬了尾跡渦結(jié)構(gòu),其得到的壓力系數(shù)也更接近試驗(yàn)值.圖15展示了翼型尾緣(x/c=0.952)和尾跡(x/c=1.282)兩個(gè)典型站位處的流向速度分布對比.在x/c=0.952處,SSG/LRR-IDDES成功捕捉到了逆向速度,計(jì)算結(jié)果與試驗(yàn)吻合.表明SSG/LRR-IDDES能在邊界層附近更好地感受逆壓梯度的影響,對于弱非定常流動(dòng)較傳統(tǒng)SST-IDDES方法有一定提高.而僅采用URANS方式會(huì)引入過多湍流黏性,使弱非定常運(yùn)動(dòng)的發(fā)展受到抑制.在x/c=1.282處,SSG/LRR-IDDES得到的流向速度略大于試驗(yàn)值,且與其他三種方法無明顯區(qū)別.

圖14 (網(wǎng)刊彩色)NACA4412翼型表面壓力系數(shù)分布Fig.14.(color online)Distribution of pressure coefficient on the surface of NACA4412 airfoil.

圖15 (網(wǎng)刊彩色)NACA4412翼型流向速度分布比較Fig.15.(color online)Comparison of streamwise velocity pro files for NACA4412 airfoil.

4.3 鈍前緣三角翼渦破裂模擬

現(xiàn)代戰(zhàn)斗機(jī)和導(dǎo)彈多采用三角翼布局,以獲得良好的飛行品質(zhì)和機(jī)動(dòng)性能.然而在大迎角下,三角翼會(huì)產(chǎn)生渦破裂現(xiàn)象對氣動(dòng)特性造成影響.本文采用的計(jì)算模型為65?后掠三角翼,是NASA Langley中心為研究該外形的雷諾數(shù)、馬赫數(shù)影響而完成的試驗(yàn)?zāi)Ｐ蚚40].該模型分為四部分:前緣、平板部分、后緣及整流罩.前緣部分提供4種可替換外形,本文選用中等半徑鈍前緣外形.來流條件為:Ma=0.85,Re=6×106.迎角選擇24.6?,為渦破裂現(xiàn)象發(fā)生的臨界迎角.

本文采用的計(jì)算網(wǎng)格自主生成,網(wǎng)格單元數(shù)約6百萬,網(wǎng)格結(jié)構(gòu)為多塊對接網(wǎng)格.圖16給出了三角翼的表面網(wǎng)格.網(wǎng)格拓?fù)洳捎肅-H型,以避免翼尖奇性軸的產(chǎn)生.計(jì)算區(qū)域的遠(yuǎn)場邊界取為50倍根弦長.壁面的第一排網(wǎng)格達(dá)到了10?6弦長,網(wǎng)格在背風(fēng)區(qū)和剪切層附近均進(jìn)行了適當(dāng)?shù)募用?以保證分離區(qū)、附面層內(nèi)和剪切層的數(shù)值模擬精度.時(shí)間推進(jìn)采用雙時(shí)間步法,周期為T=cref/U∞,其中cref為參考?xì)鈩?dòng)弦長,時(shí)間步長取0.01T.時(shí)間平均統(tǒng)計(jì)從10T時(shí)開始,持續(xù)20T.

圖16 鈍前緣三角翼計(jì)算網(wǎng)格Fig.16. Computational mesh for the blunt-edge deltawing.

圖17展示了鈍前緣三角翼不同站位處壓力分布與試驗(yàn)[41]的對比,其中η=2z/bl,bl為當(dāng)?shù)卣归L,cr為翼根弦長.回顧鈍前緣三角翼分離的基本特性.相比傳統(tǒng)尖前緣,鈍前緣三角翼在前緣處有轉(zhuǎn)捩過程[42],但本文計(jì)算的馬赫數(shù)較高,迎角也較大,因此轉(zhuǎn)捩過程很短,分離渦為全湍流狀態(tài).主渦旋轉(zhuǎn)產(chǎn)生展向流動(dòng),在旋渦下部、三角翼的上翼面加速,出現(xiàn)負(fù)壓力峰值,該點(diǎn)沿展向到翼邊緣是逆壓梯度,將誘導(dǎo)邊界層分離,產(chǎn)生二次分離和二次旋渦,在主渦負(fù)壓力峰值外側(cè)又出現(xiàn)二次負(fù)壓力峰值.觀察x/cr=0.4站位的壓力分布,SST-IDDES,SSG/LRR-URANS和SSG/LRR-IDDES均成功捕捉到了二次渦結(jié)構(gòu),而其中SSG/LRR-IDDES與試驗(yàn)值吻合最好.

當(dāng)迎角大到一定程度時(shí),主渦開始破裂.本文選擇的24.6?迎角,為渦破裂現(xiàn)象發(fā)生的臨界迎角,而x/cr=0.6是該迎角下渦破裂發(fā)生前的臨界站位.因此該站位的模擬對CFD方法是個(gè)挑戰(zhàn).觀察本文計(jì)算值與試驗(yàn)值的對比,四種方法均產(chǎn)生了不同程度的偏差.其中SSG/LRR-IDDES在吸力峰處與試驗(yàn)吻合很好,但在翼根處吻合較差.x/cr=0.6站位的翼根處臨近整流罩頭部(存在速度駐點(diǎn)),壓力系數(shù)由負(fù)值迅速變?yōu)榻咏?從而形成較大的壓力梯度,因此給數(shù)值模擬帶來很大困難.

圖17 (網(wǎng)刊彩色)鈍前緣三角翼不同站位處壓力分布對比Fig.17.(color online)Comparisons of surface pressure with experiments for blunt-edge deltawing at typical stations.

圖18 (網(wǎng)刊彩色)三角翼20T時(shí)刻Q判據(jù)等值面圖,顏色由壓力系數(shù)標(biāo)識(shí)Fig.18.(color online)Iso-surface of the Q-criterion at 20T for deltawing case,colored by pressure coefficient.

圖19 (網(wǎng)刊彩色)三角翼20T時(shí)刻流線圖,顏色由馬赫數(shù)標(biāo)識(shí)Fig.19.(color online)Streamlines at 20T for deltawing case,colored by Mach number.

在x/cr=0.8處,渦破裂已發(fā)生,上翼面吸力峰消失.SST-IDDES,SSG/LRR-URANS和SSG/LRR-IDDES均成功地預(yù)測了該現(xiàn)象,而SST-URANS仍得到了主渦及二次渦結(jié)構(gòu).結(jié)合其他站位上SST-URANS得到的壓力分布,可以看出在24.6?迎角下,其并未預(yù)測到渦破裂現(xiàn)象,即推遲了渦破裂的發(fā)生.而該現(xiàn)象的產(chǎn)生明顯由SSTURANS在分離區(qū)過大的湍流黏性導(dǎo)致.由于來流為跨聲速,前緣渦引起的上翼面流動(dòng)加速,形成了局部超音速區(qū)域并產(chǎn)生了激波.對于跨聲速流動(dòng),Rumsey[14]在ONERO M6機(jī)翼和NASA CRM翼身組合構(gòu)型中均表明,SSG/LRR-ω模型能更好地預(yù)測激波附近雷諾應(yīng)力的劇烈變化,對激波誘導(dǎo)分離的模擬能力明顯優(yōu)于SST等渦粘模型.本文發(fā)展的SSG/LRR-IDDES方法繼承了RSM的優(yōu)勢,在非定常區(qū)(x/cr=0.8和x/cr=0.95)得到了優(yōu)于SST-IDDES方法的結(jié)果.

圖18和圖19分別展示了20T時(shí)刻三角翼Q判據(jù)等值面圖和流線圖.SST-URANS得到的仍為定常狀態(tài)的前緣分離渦,在24.6?迎角時(shí)未發(fā)生渦破裂.SSG/LRR-URANS雖然預(yù)測到了渦破裂,但渦破裂后的流動(dòng)未表現(xiàn)出明顯的非定常特性.表明對渦破裂后形成的非定常運(yùn)動(dòng),需采用諸如DES,LES等的尺度模擬方法.SST-IDDES和SSG/LRR-IDDES均較好地模擬了渦破裂后的非定常流動(dòng)狀態(tài),但就精細(xì)度而言,SSG/LRRIDDES要更為優(yōu)秀.

5 結(jié) 論

本文借鑒SST-IDDES方法的構(gòu)造方式,在SSG/LRR-ωRSM上發(fā)展了SSG/LRR-IDDES分離渦模擬方法.通過結(jié)合高精度WCNS格式在NACA系列翼型及鈍前緣三角翼算例上進(jìn)行了驗(yàn)證,并和傳統(tǒng)SST-IDDES方法以及SSG/LRRURANS和SST-URANS方法進(jìn)行了對比.

為說明SSG/LRR-IDDES方法在翼型大迎角氣動(dòng)力模擬方面的提升,選取NACA0012翼型17?,45?和60?三個(gè)典型狀態(tài)進(jìn)行模擬.在17?迎角時(shí),翼型邊界層分離擴(kuò)大到前緣附近,升力系數(shù)接近最小.此時(shí),采用SST-URANS和SSG/LRR-URANS均會(huì)過高估計(jì)升力,而采用傳統(tǒng)基于線性渦粘模型的SST-IDDES方法也并未明顯改善.這說明對于傳統(tǒng)基于線性渦粘模型的IDDES方法,由于邊界層附近RANS部分本身對逆壓梯度引起的分離的模擬能力較差,使得整個(gè)方法較URANS并無明顯提高.而基于RSM,可以更準(zhǔn)確地估計(jì)逆壓梯度區(qū)的雷諾應(yīng)力變化,使得結(jié)合IDDES時(shí)得到了很好的結(jié)果.到45?和60?迎角時(shí),翼型產(chǎn)生大范圍脫體渦,升力系數(shù)回升.此狀態(tài)是傳統(tǒng)IDDES方法模擬的優(yōu)勢方面,而SSG/LRR-IDDES也繼承了該特性,并較SSG/LRR-URANS有較大提高.

為說明SSG/LRR-IDDES方法在翼型最大升力迎角附近的模擬能力,選取NACA4412翼型12?狀態(tài)進(jìn)行模擬.此時(shí),翼型尾跡由于逆壓梯度影響開始分離并產(chǎn)生湍流尾跡渦.采用SST-URANS和SSG/LRR-URANS均會(huì)過高估計(jì)升力,并延遲了尾跡分離;而采用SST-IDDES也未明顯改善.但是采用SSG/LRR-IDDES方法給出了與試驗(yàn)吻合很好的氣動(dòng)力特性和準(zhǔn)確的分離過程,表明該方法對模擬翼型失速迎角附近特性的能力提高.

為進(jìn)一步說明SSG/LRR-IDDES方法在三維機(jī)翼失速迎角附近的模擬能力,選取鈍前緣三角翼24.6?狀態(tài)進(jìn)行模擬.該迎角下,三角翼主渦由于激波誘導(dǎo)的影響,在x/cr=0.6站位后發(fā)生非定常渦破裂現(xiàn)象.采用SST-URANS方法無法成功預(yù)測渦破裂現(xiàn)象;而采用SSG/LRR-URANS方法雖然預(yù)測到了渦破裂,但渦破裂后的流動(dòng)未表現(xiàn)出明顯的非定常特性.采用SST-IDDES方法和SSG/LRR-IDDES方法均較好地模擬了渦破裂后的非定常流動(dòng)狀態(tài),但在表面壓力分布、流場精細(xì)度等方面,SSG/LRR-IDDES方法更加優(yōu)秀.

另外,本文采用的WCNS-E-5格式能夠在粗糙網(wǎng)格上獲得較高的保真度,體現(xiàn)了高精度數(shù)值方法在模擬分離流動(dòng)中的效率優(yōu)勢.但在加密網(wǎng)格計(jì)算時(shí),也出現(xiàn)了網(wǎng)格收斂性差的現(xiàn)象,如NACA0012翼型60?迎角時(shí),SSG/LRR-IDDES在粗網(wǎng)格下得到的平均氣動(dòng)力系數(shù)略優(yōu)于密網(wǎng)格.該現(xiàn)象的產(chǎn)生與IDDES方法基于網(wǎng)格長度尺度來調(diào)整湍流黏性有關(guān),也是DES類方法需要克服的問題之一[8].下一步,將繼續(xù)結(jié)合高精度格式對SSG/LRR-IDDES方法在更廣泛的流動(dòng)中進(jìn)行驗(yàn)證并提高網(wǎng)格收斂性.

[1]Slotnick J,Khodadoust A,Alonso J,Darmofal D,Gropp W,Lurie E,Mavriplis D 2014CFD Vision 2030 Study:A Path to Revolutionary Computational Aerosciences(Washington,DC:Langley Research Center,NASA)Tech.Rep.NASA/CR-2014-218178

[2]Eisfeld B,Rumsey C,Togiti V 2016AIAA J.54 1524

[3]Rumsey C 201452nd Aerospace Sciences MeetingNational Harbor,Maryland,January 13–17,2014 AIAA 2014-0201

[4]Tucker P 2006Int.J.Numer.Meth.Fluids51 261

[5]Richez F,Pape A,Costes M 2015AIAA J.53 3157

[6]Xu G,Jiang X,Liu G 2016Acta Mech.Sin.32 588

[7]Spalart P,Jou W H,Strelets M,Allmaras S 1997Comments on the Feasibility of LES for Wings,and on Hybrid RANS/LES Approach(Columbus:Greyden Press)

[8]Spalart P 2009Annu.Rev.Fluid Mech.41 181

[9]Probst A,Radespiel R,Knopp T 201120st AIAA Computational Fluid Dynamics ConferenceHonolulu,Hawaii,June 27–30,2011 AIAA 2011-3206

[10]Strelets M 200139th AIAA Aerospace Sciences Meeting and ExhibitReno,NV,8-11 January 2001,AIAA 2001-0879

[11]Greschner B,Thiele F,Gurr A,Casalino D,Jacob M200612th AIAA/CEAS Aeroacoustics ConferenceCambridge,Massachusetts,May 8–10,2006 AIAA 2006-2628

[12]Greschner B,Thiele F,Jacob M,Casalino D 2008Comput.Fluids37 402

[13]Cécora R D,Radespiel R,Eisfeld B,Probst A 2015AIAA J.53 739

[14]Rumsey C 2015 in Eisfeld B(ed.)Differential ReynoldsStress Modeling for Separating Flows in Industrial Aerodynamics(Springer Tracts Mechanical Engineering)p19

[15]Eisfeld B,Brodersen O 200523rd AIAA Applied Aerodynamics ConferenceToronto,Ontario Canada,June 6–9,2005 AIAA 2005-4727

[16]Togiti V,Eisfeld B,Brodersen O 2014J.Aircraft51 1331

[17]Dong Y D,Wang D F,Wang G X,Deng X G 2016J.National Univ.Defense Technol.38 46(in Chinese)[董義道,王東方,王光學(xué),鄧小剛2016國防科技大學(xué)學(xué)報(bào)38 46]

[18]Shu C 2003Int.J.Comput.Fluid D17 107

[19]Wang Z,Fidkowski K,Abgrall R,Bassi F,Caraeni D,Cary A,Deconinck H,Hartmann R,Hillewaert K,Huynh H,Kroll N,May G,Persson P O,van Leer B,Visbal M 2013Int.J.Numer.Meth.Fluids72 811

[20]Wang S,Deng X,Wang G,Xu D,Wang D 2016Int.J.Comput.Fluid D30 469

[21]Georgiadis N,Rizzetta D,Fureby C 2010AIAA J.48 1772

[22]Deng X,Zhang H 2000J.Comput.Phys.165 22

[23]Deng X,Liu X,Mao M,Zhang H 200517th AIAA Computational Fluid Dynamics ConferenceToronto,Ontario Canada,June 6–9,2005 AIAA 2005-5246

[24]Deng X,Mao M,Tu G,Liu H,Zhang H 2011J.Comput.Phys.230 1100

[25]Bellot G,Corrsin S 1971J.Fluid Mech.48 273

[26]Spalart P,Deck S,Shur M,Squires K,Strelets M,Travin A 2006Theor.Comput.Fluid Dyn.20 181

[27]Shur M,Spalart P,Strelets M,Travin A 2008Int.J.Heat Fluid Fl.29 1638

[28]Nonomura T,Fujii K 2009J.Comput.Phys.228 3533

[29]Liu H,Ma Y,Yan Z,Mao M,Deng X 20148th International Conference on Computational Fluid DynamicsChengdu,China,July 14–18,2014 ICCFD8-2014-0082

[30]Gang D D,Yi S H,Zhao Y F 2015Acta Phys.Sin.64 054705(in Chinese)[岡敦殿,易仕和,趙云飛2015物理學(xué)報(bào)64 054705]

[31]Schumann U 1977Phys.Fluids20 721

[32]Chassaing J,Gerolymos G,Vallet I 2003AIAA J.41 763

[33]Yossef Y 2014J.Comput.Phys.276 635

[34]Yang Y,Zha G 201646th AIAA Fluid Dynamics ConferenceWashington,D.C.,USA,June 13–17,2016 AIAA 2016-3185

[35]Shur M,Spalart P,Strelets M,Travin A 1999Proceedings of the 4th International Symposium on Engineering Turbulence Modelling and MeasurementsCorsica,France,May 24–26,1999 p669

[36]Chen M Z 2002Fundamentals of Viscous Fliud Dynamics(Beijing:Higher Education Press)p239(in Chinese)[陳懋章 2002 粘性流體動(dòng)力學(xué)基礎(chǔ) (北京:高等教育出版社)第239頁]

[37]Chen Y,Guo L D,Peng Q,Chen Z Q,Liu W H 2015Acta Phys.Sin.64 134701(in Chinese)[陳勇,郭隆德,彭強(qiáng),陳志強(qiáng),劉衛(wèi)紅2015物理學(xué)報(bào)64 134701]

[38]Wadcock A 1987Investigation of Low Speed Turbulent Scparatcd Flow Around Airfoils(Washington,DC:Ames Research Center,NASA)Tech.Rep.NASA-CR-177450

[39]Roy R,Stoellinger M 201553rd AIAA Aerospace Sciences MeetingKissimmee,Florida,January 5–9,2015 AIAA 2015-1982

[40]Luckring M,Hummel D 2013Aerosp.Sci.Technol.24 77

[41]Chu J,Luckring M 1996Experimental Surface Pressure Data Obtained on 65 deg Delta Wing Across Reynolds Number and Mach Number Ranges.Vol.3:Medium-Radius Leading Edge(Washington,DC:Ames Research Center,NASA)NASA-TM-4645-Vol-3

[42]Luckring M 2013Aerosp.Sci.Technol.24 10

High-order detached-eddy simulation method based on a Reynolds-stress background model?

Wang Sheng-Ye1)Wang Guang-Xue2)Dong Yi-Dao1)Deng Xiao-Gang1)?

1)(College of Aerospace Science and Engineering,National University of Defense Technology,Changsha 410073,China)
2)(School of Physics,Sun Yat-sen University,Guangzhou 510275,China)

20 March 2017;revised manuscript

14 May 2017)

Referring to the construction of shear stress transport-improved delayed detached-eddy simulation(SST-IDDES)method,a variant of IDDES method based on the Speziale-Sarkar-Gatski/Launder–Reece–Rodi(SSG/LRR)-ωReynoldsstress model(RSM)as Reynolds-averaged Navier-Stokes(RANS)background model,is proposed.Through combining high-order weighted compact nonlinear scheme(WCNS),the SSG/LRR-IDDES method is applied to three aeronautic cases and compared with traditional methods:SST-unsteady Reynolds-averaged Navier-Stokes(URANS),SSG/LRRURANS,and SST-IDDES.To verify the SSG/LRR-IDDES method in simulating airfoil stalled flow,NACA0012 airfoil is adopted separately at attack angles of 17?,45?and 60?.At the attack angle of 17?,SST-URANS,SSG/LRR-URANS,and SST-IDDES methods each predict a higher lift coefficient than the experimental data,while the SSG/LRR-IDDES method obtains a better lift coefficient result and a higher fidelity vortical flow structure.It indicates that the RSM can improve the prediction of RANS-mode for pressure-induced separations on airfoil surfaces in detached-eddy simulation.At the attack angles of 45?and 60?,the SSG/LRR-IDDES method captures the massively separated flow with threedimensional vortical structures and obtains a good result,which is the same as that from the traditional SST-IDDES method.To indicate the improvement of the SSG/LRR-IDDES method in simulating airfoil trailing edge separation,NACA4412 airfoil is adopted.At the attack angle of 12?(maximum lift),the trailing edge separation is mainly induced by pressure gradient.The SSG/LRR-IDDES method can predict the separation process reasonably and obtains a good lift coefficient and location of separation compared with experimental results.However,none of other methods can predict trailing edge separation.It confirms that when RSM is adopted as RANS background model in detached-eddy simulation,the ability to predict pressure-induced separation on airfoil surface is improved.For further verifying the SSG/LRR-IDDES method for simulating three-dimensional separated flow,blunt-edge deltawing at the attack angle of 24.6?is adopted.At this attack angle,the primary vortex will break,which is difficult to predict by using the SSTURANS method.For the SSG/LRR-URANS method,it predicts the vortex breakdown successfully,but the breakdown process does not show any significant unsteady characteristic.The SST-IDDES and the SSG/LRR-IDDES methods both predict a significant unsteady vortex breakdown.But in terms of the accuracy of surface pressure and the fidelity of unsteady flow,the result obtained by the SSG/LRR-IDDES method is better than by the SST-IDDES method.

turbulence flows,Reynolds stress model,detached eddy simulation,weighted compact nonlinear scheme

(2017年3月20日收到;2017年5月14日收到修改稿)

基于Speziale-Sarkar-Gatski/Launder-Reece-Rodi(SSG/LRR)-ω雷諾應(yīng)力模型發(fā)展了一類分離渦模擬方法,結(jié)合高精度加權(quán)緊致非線性格式在典型翼型及三角翼算例中進(jìn)行了驗(yàn)證,并和傳統(tǒng)基于線性渦粘模型的分離渦模擬方法進(jìn)行了對比.結(jié)果表明:基于SSG/LRR-ω模型的分離渦模擬方法,提高了原雷諾應(yīng)力模型對非定常分離湍流的模擬能力;同時(shí)相比于傳統(tǒng)基于線性渦粘模型的分離渦模擬方法,尤其是在翼型最大升力迎角和三角翼渦破裂迎角附近,該方法在平均氣動(dòng)力預(yù)測的準(zhǔn)確度、分離湍流模擬的精細(xì)度等方面更加優(yōu)秀.

10.7498/aps.66.184701

?國防科學(xué)技術(shù)大學(xué)科研計(jì)劃(批準(zhǔn)號:ZDYYJCYJ20140101)資助的課題.

?通信作者.E-mail:xgdeng2000@vip.sina.com

感謝中山大學(xué)國家超級計(jì)算廣州中心在計(jì)算資源方面提供的幫助.

PACS:47.27.–i,47.27.em,47.27.ep,47.11.BcDOI:10.7498/aps.66.184701

*Project supported by the Foundation of the National University of Defense Technology of China(Grant No.ZDYYJCYJ20140101).

?Corresponding author.E-mail:xgdeng2000@vip.sina.com