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

鄒東陽，劉君, *，鄒麗

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

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

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

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

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

激波裝配;非結構動網格;有限體積法(FVM);計算流體力學;可壓縮流動

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

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

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

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

本文在非結構激波裝配方法研究工作的基礎上，建立了一種激波裝配/捕捉混合算法，通過對網格節點定義，實現不同的求解算法之間的靈活調用。

1 計算方法

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

本節首先對原有的激波捕捉方法和非結構動網格技術進行簡單介紹，然后通過一個二維流動對激波裝配方法相關內容進行詳細說明。

1.1 激波捕捉方法

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

(1)

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

(2)

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

(3)

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

1.2 非結構動網格技術

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

(4)

式中：Kij為節點i和j之間的彈簧倔強系數；Ni為節點i鄰居節點的個數。省略中間推導過程，網格點i新坐標位置可以通過如下迭代方程進行確定：

(5)

考慮折穿和扭轉效應，彈簧倔強系數的計算公式可以寫為

(6)

式中：φ和ψ為彈簧剛度的修正因子,增加φ的值可以提高邊界處彈簧的剛度，這樣邊界節點的運動引起的內部網格變形會傳遞得更遠;β用來規定彈簧剛度，以防出現折穿現象。

當邊界節點運動位置過大時，單純地使用網格變形技術已經無法保證計算網格具有較高的網格質量，這時需要進行網格重構。重構后的網格通過網格間信息傳遞技術對流場參數進行插值。可以證明，對于時空二階格式采用信息傳遞技術不增加插值誤差[28]。

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

1.3 激波裝配方法

1.3.1 激波點標識

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

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

圖2 網格節點屬性定義Fig.2 Definition of properties of grid nodes

1.3.2 激波參數確定

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

(7)

χ={0Maτ k≤-1

1Maτ k>-1

(8)

圖3 激波節點參數確定Fig.3 Determination of shock nodes parameters

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

(9)

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

(10)

式中：au為上游區域的聲速。根據激波上游參數，可以得到

(11)

(12)

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

(13)

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

1.3.3 激波通量計算

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

(14)

1.3.4 網格點運動

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

Solid line: grids at T=t; Dashed line: grids at T=t+Δt圖4 網格節點運動示意圖Fig.4 Sketch of grid nodes motion

2 算例測試

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

2.1 超聲速圓柱繞流

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

采用POINTWISE進行網格劃分，網格類型選擇非結構的三角形單元。離散時控制參數選擇Δ=0.1L,計算域離散后總共包含1 460個均勻分布的網格單元。

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

圖5 計算區域Fig.5 Computational domain

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

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

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

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

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

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

2.2 等截面通道內激波的運動

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

圖10 流動參數設置Fig.10 Definition of flow parameters

圖11 初始計算網格Fig.11 Initial computational grids

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

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

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

2.3 激波/渦的相互作用

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

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

(15)

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

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

(16)

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

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

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

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

圖14 激波裝配法計算中隨時間推進網格變形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 結 論

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

[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格式及其在高超聲速黏性復雜流場求解中的應用[D]. 綿陽：中國空氣動力研究與發展中心，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] 劉君， 鄒東陽， 徐春光. 基于非結構動網格的非定常激波裝配方法[J]. 空氣動力學學報， 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] 劉君， 鄒東陽， 董海波. 動態間斷裝配法模擬斜激波壁面反射[J]. 航空學報， 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] 劉君， 徐春光， 白曉征. 有限體積方法和非結構動網格[M]. 北京： 科學出版社， 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">網絡出版時間

時間：2017-06-260907

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

國家自然科學基金(91541117)

.E-mailliujun65@dlut.edu.cn

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

(責任編輯：李明敏)