馬赫數變化對平板上方柱繞流特性的影響

薛大文，陳志華，孫曉暉，韓珺禮

（南京理工大學瞬態物理重點實驗室，南京 210094）

0 引言

方柱作為飛行器表面凸起的典型外形，其繞流特性與對飛行器的影響一直受到研究與關注，特別是對超聲速與高超聲速方柱繞流及其所包含的激波與邊界層相互作用等的研究，不僅具有理論價值，而且還有實際的工程應用背景。

自20世紀50年代開始，人們進行了大量的風洞試驗和飛行實驗以研究超聲速流動問題。因受實驗與計算條件的限制，其流場細節與分離機理等仍沒完全揭示清楚。早在1968年，Young等［1］實驗研究了超聲速鈍舵的激波邊界層干擾，并測得平板表面壓力分布。Dolling與Bogdonoff［2］同樣對前緣為半圓柱的鈍舵超聲速繞流進行了實驗測量，并給出了鈍舵前緣的壓力雙峰型分布以及平面表面的壓力。隨后，Shigeru等［3］研究了超聲速條件下尖舵和鈍舵繞流的激波－邊界層干擾，實驗所得的壓力分布以及油流圖譜顯示了干擾區的不穩定性，同時實驗利用彩色油流技術證實了二次分離的存在。Barberis等［4］研究比較了鈍舵后掠角對干擾流場影響，實驗所得的一系列紋影圖和油流圖譜顯示了后掠角的增大可以有效減小干擾強度和干擾范圍。王世芬等［5］對超聲速圓柱繞流進行了系列實驗研究，得到了物面平均壓力和熱流分布。而李素循［6］對高超聲速與超聲速下，有高度變化的圓柱與方柱進行了風洞實驗研究，得到了比較詳細的平板表面的壓力分布、物面油流圖譜和紋影照片。

隨著計算機技術的發展，數值模擬在此領域的應用越來越廣泛。Hung［7］用MacCormack格式先后模擬了平板上圓形鈍舵和方舵的超聲速流場，計算結果與實驗結果比較吻合，但是由于MacCormack 格式的固有缺點，結果對激波分辨率不高。Chaussee等［8］利用雷諾平均方法先后模擬了圓錐、再入式飛行器、航天飛機等繞流流場，計算了升阻力與壓力分布。同樣鄧小剛［9］基于三維雷諾平均方法，利用NND 格式對方舵超聲速繞流激波邊界層干擾進行了模擬，結果描述了弓形激波和分離激波碰撞后所形成的“λ”波系結構并分析了流場的弱非定常性。馬漢東等［10］同樣基于雷諾平均法與修正的ROE 方法，對高超聲速且壁面有傳熱條件下的高超聲速圓柱繞流進行了數值模擬，得到了上游分離區的馬蹄形五渦結構及二次分離激波。李艷麗［11］也用雷諾平均方程模擬了高超聲速鈍舵層流干擾，并與風洞試驗數據做了對比，研究了干擾流場的流動特性以及前緣后掠角對流場的影響，隨著后掠角的減小，干擾流場區域擴大，壓力載荷也逐漸增加。Shen等［12］同樣基于雷諾平均方程，利用WENO 格式對有壓縮拐角的柱／錐體高超聲速流動的激波邊界層干擾進行了研究，對比討論了WENO 格式階數對模擬精度的影響，與實驗結果的對比顯示了5階WENO 格式的精度明顯高于3階格式。

雖然前人對超聲速繞流問題做了不少數值研究，但大多基于雷諾平均方程，只能得到各流場參數的統計平均量，丟失了包含在脈動運動中的大量有重要意義的信息。大渦模擬方法則通過體積平均，對大尺度運動直接計算，而對于比網格尺度小的小尺度運動對大尺度運動的影響則通過建立模型來模擬。因而既具有較直接模擬少得多的計算量，同時還可保證模擬流場的精確性。因而本文就采用近年來發展較快的大渦模擬方法，同時利用WENO 格式來模擬前人研究較少的方柱繞流問題，結合采用自適應網格加密技術（AMR，Adaptive Mesh Refinement）對方柱超聲速繞流及其激波－邊界層干擾流場進行模擬計算，進一步研究和討論來流馬赫數對其流場結構的影響。

1 計算方法與模型

1.1 計算方法

三維可壓N－S方程經過Faver濾波，略去非線性項后，得到以下形式：連續性方程：

動量方程：

能量方程：

上述數學模型在無量綱化后，采用有限體積法進行離散。時間推進采用三階精度的Runge－Kutta法，對流項采用三階WENO 格式，而粘性項則用中心差分格式。計算過程中還對網格進行自適應加密來提高時間與空間的求解效率。計算過程中還對高壓力梯度區域采取自適應加密［13］來提高對激波的分辨率。

1.2 三維計算模型

選用平板上方柱繞流的計算模型為文獻［6］中的FD－03風洞實驗模型，具體結構與計算域如圖1 所示。計算域的長×寬×高大小為33cm×24cm×10cm，方柱根部與來流入口距離為16.8cm，方柱截面大小為2.5cm×2.5cm，高度為3cm。分別取來流馬赫數Ma＝2.9與5.0。

圖1 計算模型圖（單位：cm）Fig.1 Computational model（unit：cm）

2 計算結果與討論

圖2為Ma＝5時，方柱的實驗紋影（圖2（a））與本文的計算結果（圖2（b））比較。可知，兩者非常吻合。對于高超聲速方柱繞流（Ma＝5.0），方柱前端的二維弓形激波貼近柱體迎風面，形成激波層。同時，計算能夠清晰揭示流場中分離激波和三維弓形激波碰撞后所形成的類“λ”波系結構，以及柱體繞流過程中的復雜波系結構。

圖3為不同馬赫數（Ma＝2.9，5.0）條件下，三維方柱繞流的壓力等勢分布的切片與平板壓力分布。兩種情況下的柱體前端壓力分布總體趨勢相同，即方柱頂端前半部分與弓形激波之間以及方柱前端底部存在高壓區。同時，分離激波與方柱頂端弓形激波的相交位置基本相同。但受來流影響，高超聲速來流的弓形激波彎曲程度較低馬赫數大。而平板上壓力分布具有相同趨勢。相關討論參見本節后面。

圖2 Ma＝5.0時，方柱繞流的實驗紋影（a）與本文的計算結果（b）比較Fig.2 Comparison of experimental［6］and numerical schlieren at Ma＝5.0

圖3 馬赫數Ma＝2.9，5.0時，方柱三維繞流的壓力切片與平板壓分布圖Fig.3 Three－dimensional pressure distribution of the cylinder at Ma＝2.9and 5.0

圖4為兩種馬赫數條件下的方柱上游平板對稱線上的壓力計算值與Ma＝5.0時的壓力實驗值的比較。圖中以方柱迎風面處作為橫坐標原點，同時取x／d（其中d為方柱寬度）作為橫坐標，以p／p∞為縱坐標。可知，對于Ma＝5.0，數值模擬結果和實驗基本相符，但在x＝－1.0～－2.75區間內相差相對較大，這主要是由于此區域的旋渦的非定常發展與變化造成的，因此，此區域的壓力一直處于脈動狀態（參見圖6與圖7）。另外，柱體上游平板對稱線上的壓力同樣與測試相符，并呈雙峰值分布，最高峰值點靠近柱體表面。而在距柱體較遠的上游處，x≤－2.75，計算值與實驗值完全相符。一般地，壓力隨著與柱體的距離變小而增大，約在x／d＝－1.8處達到第一峰值。此位置點正好產生逆壓梯度而引起流動的二次分離。且在第一峰值和最高峰值之間又有極小值出現，此處應是渦核下部高速低壓區。隨后由于臨近柱面，壓力值迅速上升，由于柱體前緣根部角落渦的存在，在靠近柱面處壓力又由極大值又開始下降。來流馬赫數Ma＝2.9時，壓力雙峰的基本特征和Ma＝5.0時相似，且第一峰值相差不大，壓力抬升點也大致相同，但是第二峰值之間的差異較大，隨著馬赫數變小，數值下降較多。圖5 為方柱前緣對稱線上的壓力分布。由于分離激波位置剛好略高于方柱頂部，因此在z／d約為1.4處形成壓力極大值。同樣，計算值和實驗數據吻合良好。這在Ma＝5.0和Ma＝2.9兩種不同馬赫數下，柱體前端壓力分布趨勢基本一致。但隨著馬赫數減小，壓力峰值有所下降。

圖4 方柱上游平板對稱線壓力分布Fig.4 The pressure distribution along the symmetry axis of the flat

圖5 方柱前緣表面對稱線壓力分布Fig.5 The pressure distribution on the leading edge of cylinder

圖6和7分別為Ma＝5.0與2.9時，平板表面（xoy面）及柱體對稱面（xoz面）上的流線分布。可知，流體經過激波或遇到方柱障礙后出現了逆壓梯度，并導致分離，在柱體上游的平板表面形成了典型的分離與再附的拓撲結構。圖中S 表示分離點，A表示再附點，S2、A2、S3、S4分別表示二次、三次與四次分離和再附點。以上各圖每個分離點對應一條分離線，而再附點和再附線相對應。表面摩擦線收攏于分離線，并從再附線發散，均符合鞍點和結點交替出現的拓撲規律［14］。另外，圖6與圖7 的xoz面還顯示了不同時刻方柱前端分離區流場的形成與發展。

圖6 Ma＝5.0時，平板表面及對稱面流線圖Fig.6 Streamlines on both the surface of the flat and the symmetry plane at Ma＝5.0

圖7 Ma＝2.9時，平板表面及對稱面流線圖Fig.7 Streamlines on both the surface of the flat and the symmetry plane at Ma＝2.9

障礙物繞流分離區中的旋渦數量和旋渦結構一直是被關注且尚未完全解決的問題。Sedney［15］認為Reynolds數的不同會導致分離區內旋渦個數的變化，當二次分離現象被發現以后，Sedney又認為主渦將分叉并以偶數方式增加。后來Norman提出了噴流迷宮模型，認為主渦可以分叉為兩個同向的旋渦，但不會產生二次分離，由此產生了奇數個旋渦。鄧小剛等［11］模擬方舵激波誘導湍流邊界層分離流場結構時，得到流動發展經過2 渦－4渦－6渦－4渦的發展過程并將繼續變化。

圖6與圖7中所描述的方柱分離區的流場結構表明了分離區流動的高瞬態特性。t＝0.00025s時，馬蹄渦結構開始形成，隨著分離區的往上游發展，主渦與邊界層作用，渦管變小，形成新的逆壓梯度，并造成二次分離，產生兩個順時針渦及其中間的逆時針反向渦（見圖6（b）），與此同時在方柱前緣的角落渦以及逆時針渦的影響下，與其相鄰的兩個順時針渦被迫抬升，迫使壁面流動再次分離，最終形成圖6（c）、圖6（d）中所顯示多渦結構。結合xoy和xoz平面可以看到每一個馬蹄渦渦管都對應于一條分離線和一條再附線，隨著時間的發展，分離區中的旋渦數目由2個到4個到6個再到8個，結構也越來越復雜，進而表現為表面流線在保持拓撲結構的同時越來越扭曲。

對比圖7與圖6，兩圖中各渦的形成相似，但由于來流馬赫數變小，激波邊界層干擾明顯減弱，表現為分離區長度變小，分離區中渦系結構變得相對簡單。從圖7（d）中看到，來流馬赫數較小時，分離區的渦被抬升，因而僅在平板表面留下一條分離線和一條再附線。

為了能形象描繪方柱上游分離區流場結構的三維特性，我們結合流線與壓力等勢分布，對此區域內的主要特征進行了顯示。圖8為t＝0.001s時，不同馬赫數（Ma＝2.9，5.0）的方柱繞流分離區及附近的空間流線與壓力等勢分布。圖8（a）顯示了分離區中四個順時針方向渦管的空間流線形態，而圖8（b）則刻畫了兩個順時針方向渦形態，揭示了分離流場中各點纏繞著向下游蜿蜒，形成橫向的翼展狀結構。各流線分布趨勢大致相同，在方柱阻礙和逆時針渦作用下，從上游往下游，流線位置開始抬升。圖中三個平面上的彩色代表壓力等勢分布，其中紅色代表壓力值最高區，藍色則為最低區。結合圖6和圖7，可知每條分離線位置對應一個壓力峰值，弓形激波與柱面的迎風面及頂端前半部分則對應流場中的壓力極大值，而平板表面拓撲結構中的各條流線正是空間結構的各渦所留下的蜿蜒扭曲的“足跡”。

方柱繞流的尾渦形態同樣被關注。圖9 為t＝0.001s時，不同馬赫數（Ma＝2.9，5.0）的柱體下游尾渦空間流線圖。流線以壓力值著色。對于高超聲速方柱繞流（Ma＝5.0，圖9.（a）），貼近柱體背風面的渦對由其下部流線以螺旋形態往上抬升同時往下游發展，形成的近壁“羊角狀”尾渦對，另一對則從柱體兩側下部繞至柱體背風面后從下部往上抬升以螺旋形態沿著渦軸線往外散開并一直持續下去，從而形成遠壁“喇叭狀”尾渦對。而對于超聲速方柱繞流（Ma＝2.9，圖9（b）），柱體尾部可以清晰地看到由下往上螺旋發展的“羊角狀”渦對，而往下游的“喇叭狀”渦對則不明顯。相比可知，馬赫數較大時，柱體尾渦大幅擺動，呈現外擴形態，而小馬赫數的尾渦線則顯得狹長。圖10為對應時刻兩種馬赫數條件下，方柱尾流在平板表面上的流線分布。同樣清晰顯示了兩種情況下尾流的異同。在靠近柱體壁面流線形態大致相同。但在遠離方柱的下游，對于Ma＝5.0時，由于其螺旋形態“喇叭狀”渦對的影響，在平板表面留下了曲折的再附線，而Ma＝2.9時的再附線卻顯得直接明了和有序。

圖8 t＝0.001s時刻柱體上游的空間流線形態Fig.8 Three－dimensional streamlines in front of the cylinder at time t＝0.001s

圖9 t＝0.001s時刻，方柱尾流的空間形態Fig.9 Wake behind the cylinder at time t＝0.001s

圖10 t＝0.001s時刻，方柱尾流的平板表面流線分布Fig.10 Streamlines of the wake on the flat at time t＝0.001s

3 結論

本文通過大渦模擬方法，應用WENO 格式以及自適應網格加密技術，對來流馬赫數，Ma＝2.9，5.0條件下的三維方柱超聲速繞流進行了數值研究，與相關實驗結果相符。數值結果表明，兩種來流條件下，柱體前端平板上的壓力分布具有相同趨勢，一般隨與柱體距離變小而增大，而最高峰值點靠近柱體表面，但Ma＝2.9的壓力最高峰值下降較多。另外，對于Ma＝5.0來流，激波邊界層干擾明顯，表現為分離區長度較大，分離區中渦系結構復雜，分離區中的旋渦數目由2個到4個到6個再到8個，拓撲結構扭曲，其繞流尾渦形態同樣較為復雜，柱體尾渦大幅擺動，呈現外擴形態，而Ma＝2.9的尾渦線則顯得狹長。

［1］YOUNG F L，KAUFMAN L G，ROBERT H K.Experimental investigation of interactions between blunt fin shock waves and adjacent boundary layers at Mach numbers 3and 5［R］.OHIO：Aerospace Research Labs Wright－Patterson，1968.

［2］DOLLING D S，BOGDONOFF S M.Blunt fin－induced shock wave／turbulent boundary－layer interaction［J］.AIAAJournal，1982，20（12）：1674－1680.

［3］SHIGERU A，SYOZO M，SATOSHI O.Intensive studies on unsteady secondary separations and oscillating shock waves in three－dimensional shock waves／turbulent boundary layer interaction region induced by sharp／blunt fins［R］.AIAA－93－2939.

［4］BARBERIS D，MOLTON P.Shock wave－turbulent boundary layer interaction in a three dimensional flow［R］.AIAA－95－0227.

［5］王世芬，王宇.鈍緣舵高超音速湍流分離特性［J］.航空學報，1996，17（S1）：2－7.

［6］李素循.激波與邊界層主導的復雜流動［M］.北京：科學出版社，2007：68－80.

［7］HUNG C M，BUNING D G.Simulation of blunt－fininduced shock wave and turbulent boundary－layer interaction［J］.JournalofFluidMechanics，1985，154：163－185.

［8］DENNY S，CHAUSSE E.High－speed flow calculations past 3－D configurations based on the Reynolds averaged Navier－Stokes equations［R］.NASA－100082，1988.

［9］鄧小剛，張涵信.數值研究平板方舵激波－湍流邊界層干擾［J］.力學學報，1993，25（6）：651－657.

［10］馬漢東，李素循，吳禮義.高超聲速繞平板上直立圓柱流動特性研究［J］.宇航學報，2000，21（1）：1－5.

［11］李艷麗，李素循.高超聲速繞鈍舵層流干擾流場特性研究［J］.宇航學報，2007，28（6）：1472－1478.

［12］SHEN Y Q，ZHA G C，MANUEL A H.Simulation of hypersonic shock wave／boundary layer interaction using high order WENO scheme［R］.AIAA 2010－1047.

［13］BERGER M，COLELLA P.Local adaptive mesh refinement for shock hydrodynamics［J］.JournalofComputer＆Physics，1988，82：64－84.

［14］張涵信.分離流與旋渦運動的結構分析［M］.北京：國防工業出版社，2005：18－22.

［15］SEDNEY R.A survey of the effects of small protuberance in boundary layer flows［J］.AIAAJournal，1973，2（6）：782－792.