馬赫數(shù)變化對(duì)平板上方柱繞流特性的影響

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

（南京理工大學(xué)瞬態(tài)物理重點(diǎn)實(shí)驗(yàn)室，南京 210094）

0 引言

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

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

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

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

1 計(jì)算方法與模型

1.1 計(jì)算方法

三維可壓N－S方程經(jīng)過(guò)Faver濾波，略去非線性項(xiàng)后，得到以下形式：連續(xù)性方程：

動(dòng)量方程：

能量方程：

上述數(shù)學(xué)模型在無(wú)量綱化后，采用有限體積法進(jìn)行離散。時(shí)間推進(jìn)采用三階精度的Runge－Kutta法，對(duì)流項(xiàng)采用三階WENO 格式，而粘性項(xiàng)則用中心差分格式。計(jì)算過(guò)程中還對(duì)網(wǎng)格進(jìn)行自適應(yīng)加密來(lái)提高時(shí)間與空間的求解效率。計(jì)算過(guò)程中還對(duì)高壓力梯度區(qū)域采取自適應(yīng)加密［13］來(lái)提高對(duì)激波的分辨率。

1.2 三維計(jì)算模型

選用平板上方柱繞流的計(jì)算模型為文獻(xiàn)［6］中的FD－03風(fēng)洞實(shí)驗(yàn)?zāi)Ｐ?，具體結(jié)構(gòu)與計(jì)算域如圖1 所示。計(jì)算域的長(zhǎng)×寬×高大小為33cm×24cm×10cm，方柱根部與來(lái)流入口距離為16.8cm，方柱截面大小為2.5cm×2.5cm，高度為3cm。分別取來(lái)流馬赫數(shù)Ma＝2.9與5.0。

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

2 計(jì)算結(jié)果與討論

圖2為Ma＝5時(shí)，方柱的實(shí)驗(yàn)紋影（圖2（a））與本文的計(jì)算結(jié)果（圖2（b））比較?？芍瑑烧叻浅Ｎ呛?。對(duì)于高超聲速方柱繞流（Ma＝5.0），方柱前端的二維弓形激波貼近柱體迎風(fēng)面，形成激波層。同時(shí)，計(jì)算能夠清晰揭示流場(chǎng)中分離激波和三維弓形激波碰撞后所形成的類(lèi)“λ”波系結(jié)構(gòu)，以及柱體繞流過(guò)程中的復(fù)雜波系結(jié)構(gòu)。

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

圖2 Ma＝5.0時(shí)，方柱繞流的實(shí)驗(yàn)紋影（a）與本文的計(jì)算結(jié)果（b）比較Fig.2 Comparison of experimental［6］and numerical schlieren at Ma＝5.0

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

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

圖4 方柱上游平板對(duì)稱(chēng)線壓力分布Fig.4 The pressure distribution along the symmetry axis of the flat

圖5 方柱前緣表面對(duì)稱(chēng)線壓力分布Fig.5 The pressure distribution on the leading edge of cylinder

圖6和7分別為Ma＝5.0與2.9時(shí)，平板表面（xoy面）及柱體對(duì)稱(chēng)面（xoz面）上的流線分布。可知，流體經(jīng)過(guò)激波或遇到方柱障礙后出現(xiàn)了逆壓梯度，并導(dǎo)致分離，在柱體上游的平板表面形成了典型的分離與再附的拓?fù)浣Y(jié)構(gòu)。圖中S 表示分離點(diǎn)，A表示再附點(diǎn)，S2、A2、S3、S4分別表示二次、三次與四次分離和再附點(diǎn)。以上各圖每個(gè)分離點(diǎn)對(duì)應(yīng)一條分離線，而再附點(diǎn)和再附線相對(duì)應(yīng)。表面摩擦線收攏于分離線，并從再附線發(fā)散，均符合鞍點(diǎn)和結(jié)點(diǎn)交替出現(xiàn)的拓?fù)湟?guī)律［14］。另外，圖6與圖7 的xoz面還顯示了不同時(shí)刻方柱前端分離區(qū)流場(chǎng)的形成與發(fā)展。

圖6 Ma＝5.0時(shí)，平板表面及對(duì)稱(chēng)面流線圖Fig.6 Streamlines on both the surface of the flat and the symmetry plane at Ma＝5.0

圖7 Ma＝2.9時(shí)，平板表面及對(duì)稱(chēng)面流線圖Fig.7 Streamlines on both the surface of the flat and the symmetry plane at Ma＝2.9

障礙物繞流分離區(qū)中的旋渦數(shù)量和旋渦結(jié)構(gòu)一直是被關(guān)注且尚未完全解決的問(wèn)題。Sedney［15］認(rèn)為Reynolds數(shù)的不同會(huì)導(dǎo)致分離區(qū)內(nèi)旋渦個(gè)數(shù)的變化，當(dāng)二次分離現(xiàn)象被發(fā)現(xiàn)以后，Sedney又認(rèn)為主渦將分叉并以偶數(shù)方式增加。后來(lái)Norman提出了噴流迷宮模型，認(rèn)為主渦可以分叉為兩個(gè)同向的旋渦，但不會(huì)產(chǎn)生二次分離，由此產(chǎn)生了奇數(shù)個(gè)旋渦。鄧小剛等［11］模擬方舵激波誘導(dǎo)湍流邊界層分離流場(chǎng)結(jié)構(gòu)時(shí)，得到流動(dòng)發(fā)展經(jīng)過(guò)2 渦－4渦－6渦－4渦的發(fā)展過(guò)程并將繼續(xù)變化。

圖6與圖7中所描述的方柱分離區(qū)的流場(chǎng)結(jié)構(gòu)表明了分離區(qū)流動(dòng)的高瞬態(tài)特性。t＝0.00025s時(shí)，馬蹄渦結(jié)構(gòu)開(kāi)始形成，隨著分離區(qū)的往上游發(fā)展，主渦與邊界層作用，渦管變小，形成新的逆壓梯度，并造成二次分離，產(chǎn)生兩個(gè)順時(shí)針渦及其中間的逆時(shí)針?lè)聪驕u（見(jiàn)圖6（b）），與此同時(shí)在方柱前緣的角落渦以及逆時(shí)針渦的影響下，與其相鄰的兩個(gè)順時(shí)針渦被迫抬升，迫使壁面流動(dòng)再次分離，最終形成圖6（c）、圖6（d）中所顯示多渦結(jié)構(gòu)。結(jié)合x(chóng)oy和xoz平面可以看到每一個(gè)馬蹄渦渦管都對(duì)應(yīng)于一條分離線和一條再附線，隨著時(shí)間的發(fā)展，分離區(qū)中的旋渦數(shù)目由2個(gè)到4個(gè)到6個(gè)再到8個(gè)，結(jié)構(gòu)也越來(lái)越復(fù)雜，進(jìn)而表現(xiàn)為表面流線在保持拓?fù)浣Y(jié)構(gòu)的同時(shí)越來(lái)越扭曲。

對(duì)比圖7與圖6，兩圖中各渦的形成相似，但由于來(lái)流馬赫數(shù)變小，激波邊界層干擾明顯減弱，表現(xiàn)為分離區(qū)長(zhǎng)度變小，分離區(qū)中渦系結(jié)構(gòu)變得相對(duì)簡(jiǎn)單。從圖7（d）中看到，來(lái)流馬赫數(shù)較小時(shí)，分離區(qū)的渦被抬升，因而僅在平板表面留下一條分離線和一條再附線。

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

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

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

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

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

3 結(jié)論

本文通過(guò)大渦模擬方法，應(yīng)用WENO 格式以及自適應(yīng)網(wǎng)格加密技術(shù)，對(duì)來(lái)流馬赫數(shù)，Ma＝2.9，5.0條件下的三維方柱超聲速繞流進(jìn)行了數(shù)值研究，與相關(guān)實(shí)驗(yàn)結(jié)果相符。數(shù)值結(jié)果表明，兩種來(lái)流條件下，柱體前端平板上的壓力分布具有相同趨勢(shì)，一般隨與柱體距離變小而增大，而最高峰值點(diǎn)靠近柱體表面，但Ma＝2.9的壓力最高峰值下降較多。另外，對(duì)于Ma＝5.0來(lái)流，激波邊界層干擾明顯，表現(xiàn)為分離區(qū)長(zhǎng)度較大，分離區(qū)中渦系結(jié)構(gòu)復(fù)雜，分離區(qū)中的旋渦數(shù)目由2個(gè)到4個(gè)到6個(gè)再到8個(gè)，拓?fù)浣Y(jié)構(gòu)扭曲，其繞流尾渦形態(tài)同樣較為復(fù)雜，柱體尾渦大幅擺動(dòng)，呈現(xiàn)外擴(kuò)形態(tài)，而Ma＝2.9的尾渦線則顯得狹長(zhǎng)。

［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］.航空學(xué)報(bào)，1996，17（S1）：2－7.

［6］李素循.激波與邊界層主導(dǎo)的復(fù)雜流動(dòng)［M］.北京：科學(xué)出版社，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］鄧小剛，張涵信.數(shù)值研究平板方舵激波－湍流邊界層干擾［J］.力學(xué)學(xué)報(bào)，1993，25（6）：651－657.

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

［11］李艷麗，李素循.高超聲速繞鈍舵層流干擾流場(chǎng)特性研究［J］.宇航學(xué)報(bào)，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］張涵信.分離流與旋渦運(yùn)動(dòng)的結(jié)構(gòu)分析［M］.北京：國(guó)防工業(yè)出版社，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.