可壓縮流動激波裝配在格心型有限體積法中的應(yīng)用

鄒東陽，劉君, *，鄒麗

1.大連理工大學(xué) 航空航天學(xué)院，大連 116024 2.大連理工大學(xué) 船舶工程學(xué)院，大連 116024

可壓縮流動激波裝配在格心型有限體積法中的應(yīng)用

鄒東陽1，劉君1, *，鄒麗2

1.大連理工大學(xué) 航空航天學(xué)院，大連 116024 2.大連理工大學(xué) 船舶工程學(xué)院，大連 116024

發(fā)展了一種基于格心型有限體積方法(FVM)的激波裝配算法。通過定義網(wǎng)格節(jié)點屬性可以靈活調(diào)用激波裝配和激波捕捉計算方法。在使用激波裝配方法時，激波節(jié)點運(yùn)動速度和下游運(yùn)動速度通過Rankine-Hugoniot(R-H)關(guān)系式獲得，同時采用非結(jié)構(gòu)動網(wǎng)格技術(shù)描述激波的運(yùn)動以及調(diào)整其他網(wǎng)格節(jié)點的位置。流過激波面元的通量為上游單元的基本通量，物理概念更加清晰，通量計算也更為準(zhǔn)確。在計算過程中，網(wǎng)格節(jié)點屬性可以發(fā)生變化，以此實現(xiàn)對帶有拓?fù)渥兓鲌龅拿枋?。?shù)值試驗表明：本文提出的計算方法不但具有較高的計算精度，同時能有效地避免由于捕捉激波而出現(xiàn)的數(shù)值問題。

激波裝配;非結(jié)構(gòu)動網(wǎng)格;有限體積法(FVM);計算流體力學(xué);可壓縮流動

當(dāng)流動速度大于聲速時，流動特性會發(fā)生很大的變化，非線性效應(yīng)也更加明顯。很多在低速流動研究中應(yīng)用良好的數(shù)值方法不再適用，需要發(fā)展出針對高超聲速流動更加有效的模擬手段[1]。關(guān)于高超聲速流動的數(shù)值研究在航空航天領(lǐng)域已歷經(jīng)了數(shù)十年，并且產(chǎn)生了大量的研究成果，極大地推進(jìn)了該領(lǐng)域的發(fā)展。但是對于激波，這一在高超聲速流動中出現(xiàn)的重要物理現(xiàn)象的模擬仍然是一個重要的挑戰(zhàn)。

計算激波一個最直接的辦法就是使用激波裝配法，即將激波面看做一個間斷，使其兩側(cè)的物理量滿足Rankine-Hugoniot(R-H)關(guān)系式。這種方式在計算流體力學(xué)初期作為一種重要的研究手段被廣泛應(yīng)用于高超聲速流動的數(shù)值計算中[2-3]。但是，激波裝配方法由于網(wǎng)格拓?fù)浣Y(jié)構(gòu)的限制很難應(yīng)用于具有內(nèi)嵌激波等復(fù)雜波系的實際問題中[4]。

從20世紀(jì)80年代起，基于捕捉算法的計算流體力學(xué)得到了迅速發(fā)展[5-12]，被廣泛用于航空航天、船舶、氣象、環(huán)境、生物醫(yī)學(xué)等領(lǐng)域。不可否認(rèn)捕捉算法的提出極大地促進(jìn)了計算流體力學(xué)的發(fā)展。但是，當(dāng)面對含有強(qiáng)激波的高超聲速流動時，這些在光滑區(qū)域應(yīng)用良好的計算方法，在激波附近最終還是降為一階，并且在強(qiáng)激波波后還會有非物理振蕩現(xiàn)象出現(xiàn)[13]。這些都是由于捕捉方法本身的設(shè)計缺陷所導(dǎo)致的。相對于激波陣面，波前流動為超聲速流動。按照超聲速流動特點，上游流動不應(yīng)受到下游流動的影響，但是對于“守恒格式”的激波捕捉方法來說卻不可能做到。處于激波附近的上游網(wǎng)格單元在計算時需要使用下游信息，否則具有守恒性的流動參數(shù)無法穿過激波。所以，從這點來說捕捉方法在構(gòu)造時就與激波的物理特點相矛盾[14]。

由于激波捕捉方法存在缺陷不得不重新審視傳統(tǒng)的激波裝配方法。文獻(xiàn)[15-17]將激波裝配方法同高精度計算方法相結(jié)合，利用激波裝配方法處理激波，在光滑流域使用高精度計算方法，得到全場一致的高精度計算結(jié)果。從某種意義上而言，這種做法為未來CFD發(fā)展指明了一個方向。但是在這些文獻(xiàn)中，激波裝配仍是在結(jié)構(gòu)網(wǎng)格上進(jìn)行的，所以也只是針對弓形激波進(jìn)行處理，并未涉及到較為復(fù)雜的內(nèi)嵌激波等。激波裝配的主要問題是因為拓?fù)潢P(guān)系造成的,傳統(tǒng)的激波裝配方法主要依賴的是結(jié)構(gòu)網(wǎng)格，這使得很難通過結(jié)構(gòu)的網(wǎng)格關(guān)系來描述具有非結(jié)構(gòu)特點的激波相交等復(fù)雜問題。自然而然，就會想到如果從底層改變使用的網(wǎng)格結(jié)構(gòu)可能會大大簡化原有裝配方法所遇到的困難。Paciorri等[18-22]將Moretti[23]的浮動激波裝配方法推廣到非結(jié)構(gòu)網(wǎng)格體系下，大大簡化了傳統(tǒng)激波裝配方法，并在很多問題上得以成功應(yīng)用。Paciorri等采用的是格點型有限體積，要求激波在背景網(wǎng)格上進(jìn)行滑動。激波運(yùn)動到新的位置時，通過挖洞處理將表征激波的曲線/面嵌入計算網(wǎng)格，開始進(jìn)行計算。激波運(yùn)動后挖洞區(qū)域由原始的計算網(wǎng)格填充，對應(yīng)節(jié)點參數(shù)通過插值獲得。與他們的處理方法不同，劉君等提出基于非結(jié)構(gòu)動網(wǎng)格技術(shù)的激波裝配方法，激波屬于網(wǎng)格的一部分，激波節(jié)點的移動帶動其他網(wǎng)格節(jié)點的運(yùn)動，計算過程中無需插值，該方法被成功用于求解許多定常/非定常流動問題[24-25]。

本文在非結(jié)構(gòu)激波裝配方法研究工作的基礎(chǔ)上，建立了一種激波裝配/捕捉混合算法，通過對網(wǎng)格節(jié)點定義，實現(xiàn)不同的求解算法之間的靈活調(diào)用。

1 計算方法

從1998年開始，本文作者所在的課題組一直在從事非結(jié)構(gòu)動網(wǎng)格技術(shù)相關(guān)的工作，并形成了一套具有特色的流體計算軟件[26]。本文在原來代碼的基礎(chǔ)上進(jìn)一步開展高超聲速流動計算研究，創(chuàng)新性地將激波裝配技術(shù)與原有CFD求解器有機(jī)結(jié)合起來。為了能夠更清晰地介紹本文的工作，新的計算軟件命名為MCFs。采用MCFs進(jìn)行流動模擬時，3種形式的數(shù)值解可能出現(xiàn)。① 捕捉解：在流場中沒有激波點被定義，整個流動完全采用捕捉方法進(jìn)行計算。② 裝配解：對流場中出現(xiàn)的激波都進(jìn)行了定義，激波相關(guān)參數(shù)完全由裝配方法確定。③ 混合解：顧名思義流場中出現(xiàn)的激波部分被定義了，部分激波參數(shù)采用裝配方法定義得到。

本節(jié)首先對原有的激波捕捉方法和非結(jié)構(gòu)動網(wǎng)格技術(shù)進(jìn)行簡單介紹，然后通過一個二維流動對激波裝配方法相關(guān)內(nèi)容進(jìn)行詳細(xì)說明。

1.1 激波捕捉方法

考慮無黏流動，控制方程為時間依賴的Euler方程。ALE(Arbitrary Lagrange-Euler)描述下積分形式的Euler方程可以寫為

(1)

式中：Ω為控制單元；?Ω為控制單元界面；xc為網(wǎng)格運(yùn)動速度；Q為守恒型變量向量；Fc為對流通量向量；n為界面外法線向量；V為控制體體積；S為控制體面元面積。并且有

(2)

采用格心型有限體積法(FVM)用于空間離散，時間推進(jìn)采用四步龍格-庫塔(R-K)方法。第i個單元的控制方程可以寫成如下半離散形式：

(3)

式中：Nf為第i個控制單元所包含的面元數(shù)量；Fk、nk和Sk分別為控制體第k個面元通量、外法線方向和面積。

1.2 非結(jié)構(gòu)動網(wǎng)格技術(shù)

在本文中，網(wǎng)格變形技術(shù)使用彈簧近似法[27]。當(dāng)邊界運(yùn)動后，移動內(nèi)部網(wǎng)格節(jié)點，以此來使整個系統(tǒng)的節(jié)點達(dá)到受力平衡，從而確定網(wǎng)格點的新位置。如圖1所示，考慮彈簧系統(tǒng)內(nèi)任意節(jié)點i,節(jié)點j為與之相連的節(jié)點。根據(jù)Hooke定律，節(jié)點i受力可以寫為

(4)

式中：Kij為節(jié)點i和j之間的彈簧倔強(qiáng)系數(shù)；Ni為節(jié)點i鄰居節(jié)點的個數(shù)。省略中間推導(dǎo)過程，網(wǎng)格點i新坐標(biāo)位置可以通過如下迭代方程進(jìn)行確定：

(5)

考慮折穿和扭轉(zhuǎn)效應(yīng)，彈簧倔強(qiáng)系數(shù)的計算公式可以寫為

(6)

式中：φ和ψ為彈簧剛度的修正因子,增加φ的值可以提高邊界處彈簧的剛度，這樣邊界節(jié)點的運(yùn)動引起的內(nèi)部網(wǎng)格變形會傳遞得更遠(yuǎn);β用來規(guī)定彈簧剛度，以防出現(xiàn)折穿現(xiàn)象。

當(dāng)邊界節(jié)點運(yùn)動位置過大時，單純地使用網(wǎng)格變形技術(shù)已經(jīng)無法保證計算網(wǎng)格具有較高的網(wǎng)格質(zhì)量，這時需要進(jìn)行網(wǎng)格重構(gòu)。重構(gòu)后的網(wǎng)格通過網(wǎng)格間信息傳遞技術(shù)對流場參數(shù)進(jìn)行插值?？梢宰C明，對于時空二階格式采用信息傳遞技術(shù)不增加插值誤差[28]。

圖1 彈簧模型示意圖Fig.1 Sketch of springs analogies

1.3 激波裝配方法

1.3.1 激波點標(biāo)識

為了直觀地對本文中的激波裝配方法進(jìn)行介紹，考慮如圖2所示的一個含有激波的二維流場。在使用MCFs計算程序進(jìn)行流動計算前，首先采用三角形單元對全場進(jìn)行離散，離散后的網(wǎng)格節(jié)點被標(biāo)記為普通節(jié)點(O)和激波節(jié)點(S)兩種。激波節(jié)點和普通節(jié)點相比，在激波節(jié)點上存在兩組或多組參數(shù)，而在普通節(jié)點上有且只有一組參數(shù)。

值得說明的是：對于定常問題，可以利用捕捉方法得到一個初場，然后通過一些相關(guān)辨識技術(shù)得到初始激波位置。將符合位置條件的網(wǎng)格節(jié)點標(biāo)記為激波點，經(jīng)過迭代收斂過程后，最終得到穩(wěn)定的流場。對于非定常問題，在計算過程中有時候存在激波的生成以及湮滅，所以在計算過程中需要使用更加準(zhǔn)確的激波辨識方法。

圖2 網(wǎng)格節(jié)點屬性定義Fig.2 Definition of properties of grid nodes

1.3.2 激波參數(shù)確定

由于本文采用的捕捉方法為格心型有限體積法，因此需要通過格心參數(shù)平均得到格點參數(shù)。對于普通節(jié)點來講，只需要將與該節(jié)點相鄰的所有單元進(jìn)行加權(quán)平均即可得到位于該節(jié)點的流動參數(shù)。而對于激波節(jié)點來說，由于其包含兩組流動參數(shù)，一組為上游流動參數(shù)，一組為下游流動參數(shù)。因此在通過單元格心參數(shù)來確定格點參數(shù)之前需要判斷相鄰單元位于激波上游還是下游。在激波節(jié)點的標(biāo)記完成后，將所有節(jié)點都屬于激波節(jié)點的面元標(biāo)記為激波面元。激波面元可以看做是激波陣面的離散，這些面元可以將流動區(qū)域分割成兩個部分，一部分為低壓區(qū)，位于激波上游；一部分為高壓區(qū)，位于激波下游。通過比較激波面元兩側(cè)單元熵sL和sR的大小來判斷標(biāo)識上下游單元。熵大的單元屬于下游單元，熵小的單元屬于上游單元。在上下游單元確定之后，激波節(jié)點上游參數(shù)Vu=[ρuuupu]T和下游參數(shù)Vd=[ρdudpd]T可以通過式(7)進(jìn)行確定：

(7)

χ={0Maτ k≤-1

1Maτ k>-1

(8)

圖3 激波節(jié)點參數(shù)確定Fig.3 Determination of shock nodes parameters

同時，也可以通過激波面元法向確定出激波點的法向n為

(9)

通過上述計算方法將單元格心參數(shù)平均到格點上，得到激波節(jié)點上的上下游參數(shù)Vu和Vd。對于上游流動區(qū)域來說，其位于間斷的上游，流動信息由上游單元通過間斷流向下游。在激波坐標(biāo)系下，上游流動為超聲速流動，因此處于下游的間斷不會對上游單元產(chǎn)生影響。也就是說，在激波上游使用捕捉方法獲得的流動信息是準(zhǔn)確的，即上游流動參數(shù)Vu不需要重新確定。由于下游流動區(qū)域位于間斷下游，熵、渦量以及向前運(yùn)動的聲波等信息都是從間斷傳播過來的，所以使用捕捉方法得到的下游區(qū)域的流動計算結(jié)果都不是正確的，即下游流動參數(shù)Vd需要重新確定。

(10)

式中：au為上游區(qū)域的聲速。根據(jù)激波上游參數(shù)，可以得到

(11)

(12)

將式(11)和式(12)合成一個關(guān)于Mau,rel的方程：

(13)

可以看出，式(13)左側(cè)隨Mau,rel線性變化。因此采用牛頓公式可以很容易地求解出來流馬赫數(shù)在激波法向的相對分量Mau,rel。Mau,rel一旦確定，其他4個未知參數(shù)都可以相繼得到。

1.3.3 激波通量計算

在格心型有限體積方法中，需要求解各個單元界面的流動通量。對于普通面元來說，通過各種通量計算方法在單元界面上進(jìn)行通量分解，進(jìn)而求出各個單元的通量值。而對于激波面元來說，流過激波面元的通量計算方法要更為簡單。由前文分析可知，對于激波上游而言，其流動相對于激波面為超聲速的。激波間斷位于單元下游，可以視為上游單元的超聲速出口，因此從上游流過激波面元的通量可以表示為

(14)

1.3.4 網(wǎng)格點運(yùn)動

1.3.2節(jié)在確定激波節(jié)點下游流動參數(shù)的同時也確定了激波點的運(yùn)動速度。本文采用非結(jié)構(gòu)動網(wǎng)格技術(shù)描述激波節(jié)點的運(yùn)動。根據(jù)計算的激波點速度w以及激波點法向n可以確定在此時間步內(nèi)節(jié)點運(yùn)動的位移為w·n·Δt。如圖4所示，在激波節(jié)點的運(yùn)動確定后，通過彈簧近似方法確定其他普通網(wǎng)格節(jié)點在新時刻的位置，從而得到新時刻的計算網(wǎng)格。

Solid line: grids at T=t; Dashed line: grids at T=t+Δt圖4 網(wǎng)格節(jié)點運(yùn)動示意圖Fig.4 Sketch of grid nodes motion

2 算例測試

本文通過將提出的計算方法應(yīng)用到3個二維算例中來驗證該方法的可行性。同時，通過比較使用激波裝配和不使用激波裝配的計算結(jié)果來對該方法在計算激波時的特性進(jìn)行說明。

2.1 超聲速圓柱繞流

考慮馬赫數(shù)Ma=20的高超聲速流動經(jīng)過一個半徑為L的無限長圓柱的問題。如圖5所示，計算域為包裹半圓柱的類扇形區(qū)域。作為一個含有激波的簡單問題，文獻(xiàn)[29]對該問題進(jìn)行過研究。因此有很多結(jié)果可以用來考核本文計算結(jié)果的好壞。

采用POINTWISE進(jìn)行網(wǎng)格劃分，網(wǎng)格類型選擇非結(jié)構(gòu)的三角形單元。離散時控制參數(shù)選擇Δ=0.1L,計算域離散后總共包含1 460個均勻分布的網(wǎng)格單元。

使用激波捕捉方法獲得激波裝配的初始流場，根據(jù)捕捉流場判斷出激波的大致位置，對相應(yīng)位置上的網(wǎng)格節(jié)點進(jìn)行標(biāo)記，形成計算所需的激波節(jié)點數(shù)據(jù)鏈表。圖6中Iter0給出了標(biāo)記出的初始激波。由于激波與計算網(wǎng)格的面元相匹配，所以給出的初始激波是一條不規(guī)則的曲線。在激波節(jié)點以及激波面元定義后，采用本文提出的裝配方法計算激波面元上的通量以及激波節(jié)點運(yùn)動速度，同時采用非結(jié)構(gòu)動網(wǎng)格技術(shù)調(diào)整其他普通網(wǎng)格節(jié)點的位置以避免計算網(wǎng)格由于激波節(jié)點的運(yùn)動而發(fā)生畸變。若干步迭代后，表征激波的曲線逐步光滑，波前后的流動逐步趨于穩(wěn)定。

圖5 計算區(qū)域Fig.5 Computational domain

圖6 激波收斂過程Fig.6 Process of shock wave convergence

圖7 計算的馬赫數(shù)云圖Fig.7 Computed Mach number iso-contours

圖8 激波壁面距離Fig.8 Shock wave-wall distance

圖9 歸一化的壁面壓力Fig.9 Normalized pressure at wall

圖7給出了收斂后的馬赫數(shù)云圖。圖8和圖9分別就本文計算得到的激波到壁面距離ds和壁面壓力p/p∞兩個參數(shù)隨角度Θ變化情況，與文獻(xiàn)[29]進(jìn)行了對比。從對比結(jié)果來看，本文中的激波裝配方法獲得的結(jié)果和文獻(xiàn)中的結(jié)果吻合得較好。

總的來說，采用裝配方法得到的流場結(jié)構(gòu)較為清晰。由于對強(qiáng)激波的處理，避免了激波附近數(shù)值振蕩引起整個流場的污染，所以采用很少的計算網(wǎng)格依然能夠獲得較為光滑合理的流場分布。

2.2 等截面通道內(nèi)激波的運(yùn)動

如圖10所示，一正激波流過一等截面的通道，激波運(yùn)動速度w=2.0。波前流動參數(shù)為(ρu,pu,vu)=(1.0,0.714 3,0),波后運(yùn)動參數(shù)為(ρd,pd,vd)=(2.667,3.214 3,1.25)。盡管此問題是一個簡單的一維問題，但是在計算過程中使用二維方式進(jìn)行模擬。如圖11所示，計算區(qū)域為[0,1.0]×[0,0.5]的矩形，初始激波位于x=0.25處，初始計算網(wǎng)格包括2 864個網(wǎng)格單元和1 535個網(wǎng)格節(jié)點。

圖10 流動參數(shù)設(shè)置Fig.10 Definition of flow parameters

圖11 初始計算網(wǎng)格Fig.11 Initial computational grids

根據(jù)激波初始位置將x=0.25處的所有網(wǎng)格節(jié)點(共26個)標(biāo)記為激波節(jié)點，由這些網(wǎng)格節(jié)點組成的所有面元(共25個)則被標(biāo)記為激波面元。對于這些被標(biāo)記為激波面元的通量按照前文所述的方法進(jìn)行求解，其他的面元通量采用van Leer分裂格式求解。在計算過程中由于標(biāo)記節(jié)點會進(jìn)行運(yùn)動，計算網(wǎng)格也隨之發(fā)生變形。當(dāng)網(wǎng)格質(zhì)量較差時，通過重構(gòu)方法來更新網(wǎng)格，所以計算網(wǎng)格數(shù)量并非完全不變。

圖12給出了y=0.2時沿x方向的密度分布，其中實線表示的是t=0時刻初始流場密度分布，點劃線表示的是t=0.125 5時刻的計算結(jié)果。從計算得到的密度分布可以看出：密度參數(shù)在激波處發(fā)生跳躍，激波厚度為零符合歐拉方程上描述的激波；在激波前后密度參數(shù)都為均值，沒有出現(xiàn)任何的非物理振蕩；激波從t=0時刻的初始位置x=0.25處運(yùn)動到t=0.125 5時刻的x=0.5處，激波速度為w=2.0,與理論速度保持一致，這說明該方法具有足夠高的計算精度。

圖12 沿x方向密度分布Fig.12 Density distribution along x direction

2.3 激波/渦的相互作用

圖13 計算區(qū)域和邊界條件Fig.13 Computational domain and boundary conditions

考慮一個激波和渦的相互作用問題[22, 30-31]。圖13給出了計算區(qū)域和相應(yīng)的邊界條件，其中計算區(qū)域為[0,2L]×[0,L]的矩形。在xs/L=0.7處有一個Mas=1.21的正激波。在t=0時刻，激波上游存在一個渦，渦核中心位于(xv/L,yv/L)=(0.5,0.5)處。在上游超聲速流動的作用下，渦勻速向位于下游的激波方向運(yùn)動，其運(yùn)動速度為|u∞|。在原點位于渦核中心的極坐標(biāo)下，渦只具有切向運(yùn)動速度。由渦引起的速度擾動場可以表示為

(15)

式中：τ=r/rc是一個無量綱的半徑參數(shù)，用來表示渦的影響范圍；可以通過ε、α、rc這些參數(shù)來控制渦的形狀以及影響強(qiáng)度。

在本算例中選擇ε=0.21,α=0.204,rc/L=0.05。這樣，渦的馬赫數(shù)為

(16)

式中：a∞為超聲速來流的聲速。

激波/渦的相交問題根據(jù)激波和渦的強(qiáng)度進(jìn)行分類。按照文獻(xiàn)[30]中定義，對于本文中的組合而言(Mas=1.21,Mav=0.3)，會出現(xiàn)帶有馬赫反射的強(qiáng)相交構(gòu)型出現(xiàn)。但是，在本文中由于缺乏相關(guān)的激波辨識模塊所以只考慮主激波。

初始計算網(wǎng)格共由40 000個三角形單元和20 301個網(wǎng)格節(jié)點組成。采用兩種方法對本算例進(jìn)行模擬：一種是利用傳統(tǒng)的激波捕捉方法進(jìn)行模擬，面元通量全部采用van Leer分裂方法進(jìn)行求解，記為S-C；另一種是對于標(biāo)記的激波面元采用本文的方法計算面元通量，其他普通面元同樣采用van Leer分裂格式進(jìn)行求解，記為S-F。在S-F方法中，初始激波位于x/L=0.7處，所以將x/L=0.7處的網(wǎng)格節(jié)點標(biāo)記為激波節(jié)點，對應(yīng)的面元標(biāo)記為激波面元。與S-C相比，由于節(jié)點的移動S-F的計算網(wǎng)格在計算過程中會隨計算時間T發(fā)生變形，如圖14所示，圖中實線表示激波所在位置。

在圖15中對S-C和S-F兩種方法的計算結(jié)果進(jìn)行了對比。從對比中可以看出，在捕捉法計算激波時需要若干網(wǎng)格才能描述激波，在這些網(wǎng)格內(nèi)熵值都比較高。過激波之后，捕捉得到的流場分布較不均勻。同激波裝配結(jié)果相比，沿著流動方向捕捉結(jié)果擾動影響區(qū)域更大。從圖中可以看出，在捕捉激波的下游有一條形狀較為規(guī)則的熵帶，而在裝配方法中不存在這樣的熵帶。通過分析， 認(rèn)為這條熵帶是由于捕捉激波所產(chǎn)生的數(shù)值誤差而造成的，在激波附近產(chǎn)生的誤差會沿著流動方向向下游傳播，引起下游流場均勻性發(fā)生變化，并且沒有出現(xiàn)衰減。由于本文中只對主激波進(jìn)行了裝配，對由于激波/渦相互作用產(chǎn)生的反射激波仍然采用捕捉方法計算，因此在激波裝配結(jié)果中馬赫桿下游也會出現(xiàn)一個較寬的激波區(qū)域，并且由于兩個反射激波距離較近，捕捉得到的激波無法將兩個激波分辨出來。

圖14 激波裝配法計算中隨時間推進(jìn)網(wǎng)格變形Fig.14 Grid deformations vs time evolution for shock-fitting methods

圖15 S-C和S-F方法得到的不同時刻熵的云圖Fig.15 Entropy iso-contours at different time instants obtained by S-C and S-F methods

3 結(jié) 論

發(fā)展了一種基于格心型有限體積法的激波裝配技術(shù)。通過定義網(wǎng)格節(jié)點屬性可以靈活調(diào)用激波裝配和激波捕捉計算方法。在使用激波裝配方法時，激波節(jié)點運(yùn)動速度通過R-H關(guān)系式獲得，同時采用非結(jié)構(gòu)動網(wǎng)格技術(shù)描述激波的運(yùn)動以及調(diào)整其他網(wǎng)格節(jié)點的位置。流過激波面元的通量為上游單元的基本通量，物理概念更加清晰，通量計算也更為準(zhǔn)確。數(shù)值試驗表明：本文提出的計算方法不但具有較高的計算精度，同時也能有效地避免由于捕捉激波而出現(xiàn)的數(shù)值問題。

[1] PIROZZOLI S. Numerical methods for high-speed flows[J]. Annual Review of Fluid Mechanics, 2011, 43(1):163-194.

[2] MORETTI G, SALAS M D. Numerical analysis of viscous one-dimensional flows[J]. Journal of Computational Physics, 1970, 5(3): 487-506.

[3] SALAS M D. Shock fitting method for complicated two-dimensional supersonic flows[J]. AIAA Journal, 1976, 14(5): 583-588.

[4] 賀國宏. 三階ENN格式及其在高超聲速黏性復(fù)雜流場求解中的應(yīng)用[D]. 綿陽：中國空氣動力研究與發(fā)展中心，1994： 6-11.

HE G H. Third-order ENN scheme and its application to the calculation of complex hypersonic viscous flows[D]. Mianyang: China Aerodynamics Research and Development Center, 1994: 6-11 (in Chinese).

[5] LAX P D. Hyperbolic system of conservation laws and the mathematical theory of shock waves[C]∥SIAM Regional series on Applied Mathematics, 1973.

[6] HARTEN A. High resolution schemes for hyperbolic conservative laws[J]. Journal of Computational Physics, 1983, 49(3): 357-393.

[7] OSHER S, CHAKRAVARTHY S. High resolution schemes and the entropy condition[J]. SIAM Journal on Numerical Analysis, 1984, 21(5): 955-984.

[8] CHAKRAVARTHY S, OSHER S. A new class of high accuracy TVD scheme for hyperbolic conservation laws: AIAA-1985-0363[R]. Reston, VA: AIAA, 1985.

[9] YEE H C. Upwind and symmetric shock-capturing schemes: NASA TM-89464[R]. Washington, D.C.: NASA, 1987.

[10] ROE P L. Generalized formulation of TVD Lax-Wendroff schemes: ICASE Report No. 84-53[R]. 1984.

[11] LIU X, DENG X G, MAO M L. High-order behaviors of weighted compact fifth-order nonlinear scheme[J]. AIAA Journal, 2007, 45(8): 2093-2097.

[12] QIU J X, SHU C W. On the construction, comparison and local characteristic decomposition for high order central WENO schemes[J]. Journal of Computational Physics, 2002, 183(1): 187-209.

[13] ARORA M, ROE P L. On post-shock oscillations due to shock capturing schemes in unsteady flows[J]. Journal of Computational Physics, 1997, 130 (1): 25-40.

[14] MORETTI G. Thirty-six years of shock-fitting[J]. Computers & Fluids, 2002, 31 (4-7): 719-723.

[15] RAWAT P S, ZHONG X. On high-order shock-fitting and front-tracking schemes for numerical simulation of shock-disturbance interactions[J]. Journal of Computational Physics, 2010, 229 (19): 6744-6780.

[16] PRAKASH A, PARSONS N, WANG X, et al. High-order shock-fitting methods for direct numerical simulation of hypersonic flow with chemical and thermal nonequilibrium[J]. Journal of Computational Physics, 2011, 230 (23): 8474-8507.

[17] KOPRIVA D A. Shock-fitted multidomain solution of supersonic flows[J]. Computer Methods in Applied Mechanics & Engineering, 1999, 175 (3-4): 383-394.

[18] PACIORRI R, BONFIGLIOLI A. A shock-fitting technique for 2D unstructured grids[J]. Computers & Fluids, 2009, 38 (3): 715-726.

[19] IVANOV M S, BONFIGLIOLI A, PACIORRI R, et al. Computation of weak steady shock reflections by means of an unstructured shock-fitting solver[J]. Shock Waves, 2010, 20(4): 271-284.

[20] PACIORRI R, BONFIGLIOLI A. Shock interaction computations on unstructured, two-dimensional grids using a shock-fitting technique[J]. Journal of Computational Physics, 2011, 230(8): 3155-3177.

[21] BONFIGLIOLI A, GROTTADAUREA M, PACIORRI R, et al. An unstructured, three-dimensional, shock-fitting solver for hypersonic flows[J]. Computers & Fluids, 2013, 73(6): 162-174.

[22] BONFIGLIOLI A, PACIORRI R, CAMPOLI L. Unsteady shock-fitting for unstructured grids[J]. International Journal for Numerical Methods in Fluids, 2016, 81(4): 245-261.

[23] MORETTI G, ABBETT M. A time-dependent computational method for blunt body flows[J]. AIAA Journal, 1966, 4(12): 2136-2141.

[24] 劉君， 鄒東陽， 徐春光. 基于非結(jié)構(gòu)動網(wǎng)格的非定常激波裝配方法[J]. 空氣動力學(xué)學(xué)報， 2015， 33(1)： 10-16.

LIU J, ZOU D Y, XU C G. An unsteady shock-fitting technique based on unstructured moving grids[J]. Acta Aerodynamica Sinica, 2015, 33(1): 10-16 (in Chinese).

[25] 劉君， 鄒東陽， 董海波. 動態(tài)間斷裝配法模擬斜激波壁面反射[J]. 航空學(xué)報， 2016, 37(3): 836-846.

LIU J, ZOU D Y, DONG H B. A moving discontinuity fitting technique to simulate shock waves impinged on a straight wall[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(3): 836-846 (in Chinese).

[26] 劉君， 徐春光， 白曉征. 有限體積方法和非結(jié)構(gòu)動網(wǎng)格[M]. 北京： 科學(xué)出版社， 2016.

LIU J, XU C G, BAI X Z. Finite volume methods and unstructured dynamic grids technique[M]. Beijing: Science Press, 2016 (in Chinese).

[27] BATINA J T. Unsteady Euler airfoil solutions using unstructured dynamic meshes[J]. AIAA Journal, 1990, 28(8): 1381-1388.

[28] LIU J, BAI X Z, GUO Z. A new method for transferring flow information among meshes[J]. Computational Fluid Dynamics Journal, 2007, 15(4): 509-514.

[29] LYUBIMOV A N, RUSANOV V V. Gas flow past blunt body, Part Ⅱ Table of the gasdynamic functions: NASA TT F-715[R]. Washington, D.C.: NASA, 1973.

[30] WEEKS T M, DOSANJH D S. Sound generation by shock-vortex interaction[J]. AIAA Journal, 1967, 5(4): 660-669.

[31] GRASSO F, PIROZZOLI S. Shock-wave-vortex interactions: Shock and vortex deformations, and sound production[J]. Theoretical and Computational Fluid Dynamics, 2000, 13(6): 421-456.

Applicationsofshock-fittingtechniqueforcompressibleflowincell-centeredfinitevolumemethods

ZOUDongyang1,LIUJun1, *,ZOULi2

1.SchoolofAeronauticsandAstronautics,DalianUniversityofTechnology,Dalian116024,China2.SchoolofNavalArchitectureandOceanEngineering,DalianUniversityofTechnology,Dalian116024,China

Ashock-fittingtechniqueforcell-centeredFiniteVolumeMethod(FVM)isdevelopedinthiswork.Itisflexibletoswitchamongshock-fittingandshock-capturingmethodsbychangingthenatureofgridnodes,whicharedefinedasshocknatureandcommonnature.Intheshock-fittingmethod,velocitiesofshocknodesanddownstreamstatesareobtainedbysolvingRankine-Hugoniot(R-H)relations.Theunstructureddynamicgridtechniqueisusedforshocktrackingandupdatingthepositionsofothercommonnodes.Thefluxacrossashockfaceequalsthebasicfluxofitsupstreamcell.Duringthecomputationalprocess,thenatureofthenodesisallowedtochange.Thus,itiseasiertoapplythismethodincomplexproblems,evenwithtopologicalchange.Thenumericalresultsshowtheproposedmethodisnotonlyofhighaccuracy,butalsoabletoavoidthetroublesinshock-capturing.

shock-fitting;unstructureddynamicgrid;FiniteVolumeMethod(FVM);computationalfluiddynamics;compressibleflow

2017-04-27；Revised2017-06-09；Accepted2017-06-19；Publishedonline2017-06-260907

URL：http://hkxb.buaa.edu.cn/CN/html/20171113.html

NationalNaturalScienceFoundationofChina(91541117)

.E-mailliujun65@dlut.edu.cn

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.121363

V211；O354.5

A

1000-6893(2017)11-121363-11

2017-04-27；退修日期2017-06-09；錄用日期2017-06-19；< class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2017-06-260907

http://hkxb.buaa.edu.cn/CN/html/20171113.html

國家自然科學(xué)基金(91541117)

.E-mailliujun65@dlut.edu.cn

鄒東陽，劉君，鄒麗．可壓縮流動激波裝配在格心型有限體積法中的應(yīng)用J． 航空學(xué)報,2017,38(11):121363.ZOUDY,LIUJ,ZOUL.Applicationsofshock-fittingtechniqueforcompressibleflowincell-centeredfinitevolumemethodsJ.ActaAeronauticaetAstronauticaSinica,2017,38(11):121363.

(責(zé)任編輯：李明敏)