王 亮，許厚謙，吳 偉，薛 銳

（南京理工大學(xué)能源與動力工程學(xué)院，江蘇南京 210094）

一種新的計算非穩(wěn)態(tài)流動的自適應(yīng)無網(wǎng)格方法

王 亮＊，許厚謙，吳 偉，薛 銳

（南京理工大學(xué)能源與動力工程學(xué)院，江蘇南京 210094）

基于帶權(quán)點填充布點方法，結(jié)合指定流場參數(shù)梯度的分布，提出了一種新的無網(wǎng)格自適應(yīng)方法計算非穩(wěn)態(tài)問題。該方法的自適應(yīng)探測器E由節(jié)點的權(quán)重和參數(shù)梯度組成。使用自適應(yīng)探測器識別出加密或者粗化的區(qū)域，形成加密空腔和粗化空腔，然后重新填充布點。新節(jié)點的權(quán)重由自適應(yīng)探測器E結(jié)合該處的梯度大小通過預(yù)測－校正迭代算法計算得到。另外，由于參數(shù)梯度在激波附近波動范圍比較大，所以新節(jié)點的權(quán)重存在最大值和最小值。首先對預(yù)先設(shè)置壓強梯度的流場區(qū)域進行自適應(yīng)布點，并與傳統(tǒng)自適應(yīng)布點結(jié)果進行對比，點云分布圖顯示所提自適應(yīng)方法重新生成的點云結(jié)構(gòu)疏密分布更加合理。進一步，將此自適應(yīng)方法應(yīng)用于Riemann問題和激波碰撞圓柱問題，計算結(jié)果表明該方法可以在節(jié)點數(shù)目較低的情況下顯著提高運動激波的分辨率。在激波碰撞圓柱問題中，比較了使用自適應(yīng)算法和非自適應(yīng)算法得到相當(dāng)?shù)慕Y(jié)果所使用時間，前者是后者的20.6%，因此該自適應(yīng)算法在計算效率方面也具有較大的優(yōu)勢。

無網(wǎng)格方法；自適應(yīng)方法；激波；點云

0 引 言

在計算流體力學(xué)中，無網(wǎng)格方法越來越受到關(guān)注［14］，這主要是由于無網(wǎng)格方法空間離散相對靈活，易于處理邊界較為復(fù)雜或存在復(fù)雜邊界變形的流動問題。無網(wǎng)格方法空間離散是將計算區(qū)域劃分為一系列的節(jié)點，節(jié)點之間無線段相連，避免了多邊形單元的生成，減少了存儲計算信息所占用的內(nèi)存空間。然而無網(wǎng)格方法也存在一定的缺陷：首先是其精度相對于網(wǎng)格方法較低，當(dāng)節(jié)點間距過疏時，計算結(jié)果誤差較大甚至無法分辨激波間斷；另外是計算效率較低，主要是由于形函數(shù)的計算涉及多次矩陣求逆、矩陣相乘等運算，而且每個離散單元與周圍單元間的通量計算個數(shù)（二維空間，一般大于等于6）也大于網(wǎng)格方法（二維空間，一般為3或者4）。所以，研究無網(wǎng)格快速算法是現(xiàn)在無網(wǎng)格方法研究重點之一。

自適應(yīng)技術(shù)是一種被廣泛使用于提高計算效率的方法，特別是在網(wǎng)格方法中［5－7］。在無網(wǎng)格方法中，大量學(xué)者［8－12］也對自適應(yīng)技術(shù)做作了較多的研究：馬志華［13］和蒲賽虎［14］根據(jù)流場中壓強梯度的大小，在原有節(jié)點的基礎(chǔ)上進行快速的局部加密或者粗化；馬志華［15］隨后又基于移動節(jié)點的思想提出了點云自適應(yīng)移動技術(shù)，該方法將流場中的離散點向梯度較大的區(qū)域聚集，主要適用于穩(wěn)態(tài)問題；倪國喜等［16］根據(jù)流場參數(shù)的分布，結(jié)合Voronoi圖方法對整個計算區(qū)域重新布點，清晰展現(xiàn)了點云基于參數(shù)大小的自適應(yīng)分布；Jun等［17］采用再生核粒子法（RKPM）自身具有多尺度分辨率的特性，將解分為高階模態(tài)和低階模態(tài)，在大梯度的區(qū)域利用高階模態(tài)過濾器進行精細處理；Liu等［18］結(jié)合事后誤差估計和Voronoi圖對流場進行加密。

本文基于梯度準則和帶權(quán)點填充布點技術(shù)［19］提出了一種新的無網(wǎng)格自適應(yīng)方法。該方法采用流動參數(shù)梯度大小和節(jié)點權(quán)重相結(jié)合的自適應(yīng)探測器對計算流場進行探測，形成加密空腔或粗化空腔。點云加密或粗化是在完全摒棄原有點云的基礎(chǔ)上重新生成的，充分地體現(xiàn)了參數(shù)的不均勻分布對點云結(jié)構(gòu)的影響，簡單靈活，穩(wěn)態(tài)和非穩(wěn)態(tài)問題均適用。為了驗證所提無網(wǎng)格自適應(yīng)方法在自適應(yīng)區(qū)域內(nèi)點云分布的合理性，本文首先對預(yù)先設(shè)置壓強梯度的流場進行自適應(yīng)布點，然后使用該方法分別計算非穩(wěn)態(tài)的黎曼問題和激波碰撞圓柱問題。

1 控制方程

本文采用無黏Euler控制方程，在直角坐標系中可表示為：

式中，U是守恒變量，F(xiàn)及G為無粘通量，分別為：

式中，ρ為密度，u、v為運動速度X和Y方向的分量，p為氣體壓強，E為單位質(zhì)量總能。狀態(tài)方程采用剛性氣體狀態(tài)方程。

式中，P∞和γ分別是壓力擴展系數(shù)和氣體比熱比，氣體為空氣時，P∞＝0，γ＝1.4。

2 數(shù)值方法

本文采用最小二乘方法逼近空間導(dǎo)數(shù)，基函數(shù)采用一階線性函數(shù)，其具體形式如下：

以流場中某一點i為例，其衛(wèi)星點個數(shù)為N，此點云中N＋1個點均滿足式（3）。所以可以通過一定的運算得到如下矛盾方程組：

則，控制方程（1）中的通量項可以表示為：

HLLC是一種近似求黎曼解的方法，由Toro等［20］對Harten及Lax等［21］提出的HLL方法進行改進所得。吳偉［19］將該方法引入到無網(wǎng)格方法中，計算結(jié)果中激波的分辨率比較高。所以本文采用文獻［19］中的HLLC格式計算數(shù)值通量。

3 自適應(yīng)技術(shù)

3.1 帶權(quán)點填充布點技術(shù)概述

帶權(quán)點填充布點技術(shù)是由吳偉［19］提出的一種新的無網(wǎng)格布點方法。該方法的主旨是賦予每一離散節(jié)點一個實數(shù)的權(quán)重，其物理意義可以理解為以該點為圓心，權(quán)重大小為半徑的圓球面。帶權(quán)點填充布點就是將具有相同或不同半徑的圓球面填充到計算區(qū)域中。在實際操作中，為了加強該方法的靈活可用性，兩個相鄰的虛擬圓球面之間并不要求嚴格意義上的相切，即允許部分的相交或相離。其具體的布點過程分為以下幾步：

1）根據(jù)實際情況，對邊界進行均勻或按照某一規(guī)則非均勻布點，每一個邊界點i的權(quán)重ri＝為點i與相鄰節(jié)點j的距離，N為相鄰點的個數(shù)；

2）將邊界作為初始陣面，向內(nèi)部填充帶權(quán)點。以陣面邊AB為例，其最佳推進帶權(quán)點C的權(quán)重取兩端點A和B權(quán)重的平均，通過節(jié)點C所在的球面與節(jié)點A和B所在的球面相切的原理，得到帶權(quán)點C的坐標。然后在C點周圍尋找現(xiàn)有的帶權(quán)點，形成候選點鏈表，經(jīng)過穿透、交叉等幾何判斷，刪除無效的帶權(quán)點。最后根據(jù)現(xiàn)有帶權(quán)點與理想推進點的重疊系數(shù)大小確定最佳推進點。具體的計算方法見文獻［19］。

3）填充布點結(jié)束后對內(nèi)部節(jié)點進行Laplace光順處理，使得流場節(jié)點分布更加合理。

3.2 自適應(yīng)方法

在以上帶權(quán)點填充布點技術(shù)的基礎(chǔ)上，本文提出了一種新的自適應(yīng)布點方法。該方法的自適應(yīng)探測器采用參數(shù)梯度大小與節(jié)點權(quán)重相結(jié)合的方法確定，具體的表達式為：

1）計算流場中每一節(jié)點的探測器E的取值，并對整個區(qū)域內(nèi)的探測器值取平均，得到E0，即該流場的自適應(yīng)基準值。

其中，M為計算域內(nèi)節(jié)點個數(shù)。

2）對整個區(qū)域內(nèi)的點進行加密或粗化判斷：當(dāng)ri＞1.1αrmax，Ei＞E0時，將點i標注為加密待刪除點；當(dāng)ri＜0.85βrmax，Ei＜0.1E0時，將點i標注為粗化待刪除點。然后分別形成加密空腔和粗化空腔。

3）對形成的空腔進行帶權(quán)點填充布點。新填充點C的權(quán)重采用以下公式計算。

4）填充布點結(jié)束后對新增加的內(nèi)部節(jié)點進行Laplace光順處理，使得流場節(jié)點分布更加合理。

4 數(shù)值算例

4.1 布點示例

本文首先對預(yù)先設(shè)置壓強梯度值的流場進行自適應(yīng)布點，以檢驗采用本文自適應(yīng)方法重構(gòu)點云分布的合理性。布點區(qū)域為一邊長為4的正方形，中心點位于原點（0，0）。其中，本文假設(shè)在以中心點為圓心，半徑為1的圓形區(qū)域內(nèi)壓強梯度值為外部區(qū)域的壓強梯度值均為零。布點的最大權(quán)重為0.1，最小權(quán)重為0.0125，E0＝5。圖1（a）給出了相應(yīng)的布點結(jié)果，圖1（b）給出了采用傳統(tǒng)自適應(yīng)方法［14］的布點結(jié)果。由于傳統(tǒng)方法是在原有節(jié)點的基礎(chǔ)上在兩節(jié)點中心處添加新節(jié)點或者刪除新增加的點，所以要進行多次自適應(yīng)布點才能使得大梯度區(qū)域點的權(quán)重達到0.0125。通過觀察可以發(fā)現(xiàn)傳統(tǒng)方法的節(jié)點分布呈現(xiàn)階梯型，節(jié)點的自適應(yīng)相對簡單，無法準確地與梯度的分布相吻合，冗余節(jié)點的產(chǎn)生是無法避免的，而采用本文的方法進行自適應(yīng)布點時，節(jié)點疏密過渡更加自然，與梯度的分布情況更加吻合，能夠最大程度的減少或避免冗余節(jié)點的產(chǎn)生。

圖1 自適應(yīng)布點結(jié)果Fig.1 Adaptive point distributions

4.2 Riemann問題

Riemann問題是經(jīng)典的非穩(wěn)態(tài)激波管問題。激波管內(nèi)放置一薄膜，薄膜兩側(cè)分別充滿不同狀態(tài)的空氣。本文采用二維軸對稱模型進行求解。空氣的初始狀態(tài)值為：

0.0≤x≤0.5：

ρ＝1.0，u＝0.0，p＝1.0

0.5＜x≤1.0：

ρ＝0.125，u＝0.0，p＝0.1

t＝0時刻，薄膜破裂，由于密度壓力的不同，空氣開始流動。

初始時刻，計算域內(nèi)均勻分布的無網(wǎng)格節(jié)點數(shù)為1234，如圖2（a）所示，采用非自適應(yīng)算法得到的計算結(jié)果如圖3所示，激波的捕捉精度較低，厚度較大，無法達到計算要求。所以本文采用自適應(yīng)無網(wǎng)格方法對計算域內(nèi)點云進行實時加密及粗化。對于該算例，自適應(yīng)探測器公式中的流場參數(shù)梯度取密度梯度。圖2（b）給出了t＝0.2時刻計算域內(nèi)自適應(yīng)點云分布圖，在密度變化較為劇烈的區(qū)域，點云分布的很密，并且左側(cè)的膨脹波區(qū)域點云分布由疏到密，右側(cè)強激波附近點云分布密集，對應(yīng)的激波管內(nèi)密度及壓力分布曲線如圖3所示。比較自適應(yīng)結(jié)果和非自適應(yīng)結(jié)果，可以發(fā)現(xiàn)激波的分辨率得到了很大的改善，而且自適應(yīng)結(jié)果與解析解吻合的較好。該算例說明本文的自適應(yīng)方法可以準確捕捉非定常流動中的簡單運動激波。

圖2 計算域內(nèi)t＝0.0和t＝0.2時刻節(jié)點分布Fig.2 Gridless points distributions att＝0.0andt＝0.2

圖3 在t＝0.2時刻的壓力與密度分布曲線Fig.3 Curves of pressure and density att＝0.2

4.3 瞬時激波碰撞圓柱

在以上計算的基礎(chǔ)上，本算例給出了自適應(yīng)無網(wǎng)格方法在二維非穩(wěn)態(tài)無黏圓柱繞流問題中的應(yīng)用。計算區(qū)域為長度L＝6.0，高度是H＝3.0的管道，圓柱位于管道中心處，半徑為1.0。管道內(nèi)初始激波的位置在圓柱左側(cè)0.45處，激波左側(cè)流場狀態(tài)為：ρ＝2.667，u＝1.479，v＝0.0，p＝4.5；激波右側(cè)流場狀態(tài)為：ρ＝1.0，u＝0.0，v＝0.0，p＝1.0。計算采用的初始無網(wǎng)格節(jié)點數(shù)為4703，均勻分布。

圖4給出了隨著時間的推進，計算域內(nèi)的壓力等值線圖以及對應(yīng)無網(wǎng)格點云分布的變化。觀察發(fā)現(xiàn)點云動態(tài)分布合理，在激波附近一直比較密，激波運動過后的區(qū)域點云被及時粗化，激波厚度始終較小，由此可以說明本文的自適應(yīng)算法可以實時捕捉較為復(fù)雜的運動激波。圖5給出了不采用自適應(yīng)方法，計算域內(nèi)較密集地均勻分布76389個無網(wǎng)格節(jié)點時，計算得到的t＝1.2時刻壓力等值線圖。比較發(fā)現(xiàn)兩種方式得到的激波位置及強度吻合的很好，激波分辨率很高。比較兩種算法所使用的時間，自適應(yīng)算法所使用的時間只占直接采用較密節(jié)點計算的20.6%，所以本文自適應(yīng)算法也可以有效提高計算效率。本文計算效率低于文獻［14］，這是由于本文自適應(yīng)點云重構(gòu)需要耗費較長的時間。然而，比較本文和文獻［14］中t＝1.2時刻自適應(yīng)點云分布圖，在激波附近本文點云分布更加緊密，其他區(qū)域也相對比較稀疏，說明本文方法可以更加準確地探測到計算域內(nèi)需要自適應(yīng)的區(qū)域，并對其進行重新布點，使得參數(shù)變化較為劇烈的區(qū)域內(nèi)的點云分布更加密集。

圖4 在t＝0.9、1.2、1.5時刻的自適應(yīng)結(jié)果（左側(cè)是壓力等值線圖，右側(cè)是自適應(yīng)無網(wǎng)格節(jié)點分布圖）Fig.4 Adaptive results att＝0.9，1.2，1.5（Left：pressure contours，Right：adaptive gridless points distributions）

圖5 節(jié)點總數(shù)為76389時，直接計算得到的t＝1.2時刻壓力等值線Fig.5 Pressure contours att＝1.2 when the number of gridless points is 76389

5 結(jié) 論

在帶權(quán)點填充布點的基礎(chǔ)上，本文結(jié)合流場中指定流場參數(shù)的梯度分布規(guī)律發(fā)展出了一種新的自適應(yīng)無網(wǎng)格方法。該方法在自適應(yīng)區(qū)域重新生成點云，簡單，靈活，實用，能夠準確求解非定常問題。為驗證自適應(yīng)后節(jié)點分布的魯棒性，首先針對預(yù)先設(shè)置壓強梯度分布的計算域進行了自適應(yīng)布點，得到的無網(wǎng)格自適應(yīng)點云分布圖顯示這種方法能夠根據(jù)梯度的大小對計算域合理地進行加密或粗化。隨后，非穩(wěn)態(tài)的Rimeann問題和瞬時激波碰撞圓柱問題的計算結(jié)果顯示本文方法能夠準確識別激波所在的區(qū)域，并能實時捕捉運動激波，而且點云的自適應(yīng)結(jié)果合理，進一步展示了所提自適應(yīng)方法在非穩(wěn)態(tài)問題計算中的有效性。在瞬時激波碰撞圓柱問題中，本文比較了自適應(yīng)計算和非自適應(yīng)計算得到相當(dāng)?shù)慕Y(jié)果所耗費的時間，結(jié)果顯示自適應(yīng)算法的計算時間比較小，所以本文的算法是可以較好地提高無網(wǎng)格計算效率的。本文后續(xù)的工作將考慮改進點云重構(gòu)方法以進一步提高自適應(yīng)計算效率并將該自適應(yīng)方法應(yīng)用于求解存在復(fù)雜邊界變形的多介質(zhì)流動問題。

［1］Sun Yingdan，Wang Gang，Ye Zhengyin.Extending the AUSM＋－up scheme into the gridless method［J］.Acta Aerodynamica Sinica，2005，4（4）：511－515.（in Chinese）孫迎丹，王剛，葉正寅，等.AUSM＋－up格式在無網(wǎng)格算法中的推廣［J］.空氣動力學(xué)學(xué)報，2005，4（4）：511－515.

［2］Hong Luo，Joseph D Baum，Rainald Lhner.A hybrid cartesian grid and gridless method for compressible flows［J］.Journal of Computational Physics，2006，214（2）：618－632.

［3］Ma Zhihua，Chen Hongquan，et al.A new type hybrid algorithm using local gridless techniques［J］.Acta Aerodynamica Sinica，2008，26（03）：344－348.（in Chinese）馬志華，陳紅全，等.基于局部無網(wǎng)格的混合算法研究［J］.空氣動力學(xué)學(xué)報，2008，26（03）：344－348.

［4］Wu Wei，Xu Houqian，Wang Liang，et al.Numerical simulation for the external combustion of base－bleed projectile using gridless method［J］.Journal of Harbin Institute of Technology，2014，21（3）：95－100.

［5］Kang Hongwen，Gu Xiangqian，Liu Chongjian.Study on weigh functions of adaptive grid technique［J］.Acta Mechanica Sinica，2002，34（5）：790－795.（in Chinese）康紅文，谷湘潛，柳崇健.自適應(yīng)網(wǎng)格技術(shù)中的權(quán)函數(shù)問題研究［J］.力學(xué)學(xué)報，2002，34（5）：790－795.

［6］Han Zhirong，Lu Zhiliang，et al.Grid adaption for separation flow［J］.Acta Aerodynamica Sinica，2012，1（1）：86－89.（in Chinese）韓志熔，陸志良，郭同慶，等.一種用于分離流動的網(wǎng)格自適應(yīng)算法［J］.空氣動力學(xué)學(xué)報，2012，1（1）：86－89.

［7］Xu Heyong，Ye Zhengyin，Zhang Weiwei.Numerical simulation of hypersonic flow based on unstructured adaptive grid method［J］.Engineering Mechanics，2012，29（3）：226－229.（in Chinese）許和勇，葉正寅，張偉偉.基于非結(jié)構(gòu)自適應(yīng)網(wǎng)格技術(shù)的高超聲速流動數(shù)值模擬［J］.工程力學(xué)，2012，29（3）：226－229.

［8］Lee，Sangho，Hyo Jinkim，Sukky Jun.Two scale meshfree method for the adaptivity of 3－D stress concentration problems［J］.Computational Mechanics，2000，26（4）：376－87.

［9］You Y，Chen J S，Lu H.Filters，reproducing kernel，and adaptive meshfree method［J］.Computational Mechanics，2003，31（3－4）：316－326.

［10］Lee C K，Zhou C E.On error estimation and adaptive refinement for element free Galerkin method：Part II：Adaptive refinement［J］.Computers and Structures，2004 82（4－5）：429－443.

［11］Rossi R，Alves M K.Recovery based error estimation and adaptivity applied to a modified element－free Galerkin method［J］.Computational Mechanics，2004，33（3）：194－205.

［12］Ma Z H，Hong Q C.Mesh free adaptive algorithm for solving Euler equations on structured grid points［J］.Transactions of Nanjing University of Aeronautics and Astronautics，2005，22（4）：271－275.

［13］Ma Z H，Chen H Q.Simulations of transonic inviscid flows over airfoilsusing meshfree adaptive algorithm［J］.Modern Physics Letters B，2005，19（28－29）：63－66.

［14］Pu S H，Chen H Q.Gridless adaptive method for simulating flows with moving shocks［J］.Aeronautical Computing Technique，2012，42（5）：4－8.（in Chinese）蒲賽虎，陳紅全.模擬存在運動激波流場的無網(wǎng)格自適應(yīng)算法［J］.航空計算技術(shù)，2012，42（5）：4－8.

［15］Ma Z，Chen H，Zhou C.A study of point moving adaptivity in gridless method［J］.Computer Methods in Applied Mechanics＆Engineering，2008，197（21）：1926－1937.

［16］Ni G X，Wang R L，Lin Z.Adaptive distribution of particles in a meshfree method［J］.Chinese Journal of Computational Physics，2006，23（4）：419－424.（in Chinese）倪國喜，王瑞利，林忠.無網(wǎng)格方法中粒子分布與自適應(yīng)研究［J］.計算物理，2006，23（4）：419－24.

［17］Jun S，Im S.Multiple－scale meshfree adaptivity for the simulation of adiabatic shear band formation［J］.Computational Mechanics，2000，25（2－3）：257－266.

［18］Liu G R，Tu Z H.An adaptive procedure based on background cells for meshless methods［J］.Computer Methods in Applied Mechanics ＆Engineering，2002，191（s17－18）：1923－1943.

［19］Wu W.Gridless method for numerical simulation of complex chemical reaction flow field and its application［D］.Nanjing：Nanjing University of Science ＆Technology，2015.（in Chinese）吳偉.復(fù)雜化學(xué)反應(yīng)流場數(shù)值模擬的無網(wǎng)格方法及應(yīng)用［D］.南京：南京理工大學(xué)，2015.

［20］Toro E F，Spruce M，Speares W.Restoration of the contact surface in the HLL－Riemann solver［J］.Shock Waves Journal，1994，4：25－34

［21］Harten A，Lax P D，Leer B V.On upstream differencing and Godunov－type schemes for hyperbolic conservation laws［J］.SIAM Revision，1983，25：35－6.

A novel gridless adaptive method for simulating unsteady flows

Wang Liang＊，Xu Houqian，Wu Wei，Xue Rui

（School of Energy and Power Engineering，Nanjing University of Science ＆Technology，Nanjing 210094，China）

Based on the weighted－point filling strategy，a novel gridless adaptive method is developed to capture shock waves with high resolutions in unsteady flows.The adaptive indicatorEis the function of the local specified parameter gradients and the weight of points.The refining and coarsening regions are identified by the indicator，and then form refining cavities and coarsening cavities which are filled with new points.Using the predictor－corrector iterative algorithm，the weight of the new point is obtained by the adaptive indicatorEand the local gradient.To prevent the violent change of the point weights with the parameter gradients，the maximum and the minimum are stipulated respectively in the filling process.The proposed adaptive method is firstly used to fill points into the flow field which is set with designated pressure gradient in advance.By comparison with the traditional adaptive result，the sparse and dense distributions of clouds of points are more satisfactory and more reasonable.Further studies are the calculations of Riemann flow and incident shock past a cylinder.The results show that the adaptive method can capture the moving shock waves successfully in the unsteady flow fields with scanty points.For achieving equivalent pressure contours of the flow of the incident shock past a cylinder，the time spent using the presented adaptive method is 20.6%of the conventional unadaptive method，which suggests that it has a definite advantage on the computational efficiency.

gridless method；adaptive strategy；shocks；clouds of points

V211.3

Adoi：10.7638／kqdlxxb－2015.0185

0258－1825（2016）04－0497－06

2015－10－10；

2015－11－26

王亮＊（1989－），男，安徽阜陽人，博士研究生，計算流體力學(xué).E－mail：310082216＠njust.edu.cn

王亮，許厚謙，吳偉，等.一種新的計算非穩(wěn)態(tài)流動的自適應(yīng)無網(wǎng)格方法［J］.空氣動力學(xué)學(xué)報，2016，34（4）：497－502.

10.7638／kqdlxxb－2015.0185 Wang L，Xu H Q，Wu W，et al.A novel gridless adaptive method for simulating unsteady flows［J］.Acta Aerodynamica Sinica，2016，34（4）：497－502.