低亞聲速和跨聲速矩形柱繞流的大渦模擬研究

陳建平，黃偉希，許春曉

（清華大學(xué) 航天航空學(xué)院，北京 100084）

低亞聲速和跨聲速矩形柱繞流的大渦模擬研究

陳建平，黃偉希，許春曉

（清華大學(xué) 航天航空學(xué)院，北京 100084）

采用基于非結(jié)構(gòu)網(wǎng)格的有限體／有限元混合格式和大渦模擬的方法求解可壓縮Navier-Stokes方程，研究了不同長(zhǎng)寬比矩形柱低亞聲速和跨聲速繞流的流動(dòng)特性。在雷諾數(shù)為22000時(shí)，對(duì)來(lái)流馬赫數(shù)等于0．1和0．75，截面長(zhǎng)寬比分別為1∶1、2∶1、3∶1和4∶1的矩形柱繞流進(jìn)行了大渦模擬，以研究長(zhǎng)寬比和壓縮性對(duì)矩形柱繞流流場(chǎng)的影響。馬赫數(shù)為0．1時(shí)，Strouhal數(shù)隨著長(zhǎng)寬比的增大先降低再增大然后再降低；長(zhǎng)寬比為3∶1和4∶1時(shí)會(huì)有流動(dòng)的再附產(chǎn)生；柱體上表面的三維特性在長(zhǎng)寬比大時(shí)更明顯。馬赫數(shù)為0．75時(shí)，Strouhal數(shù)隨著長(zhǎng)寬比的增大逐漸減?。煌牧髅}動(dòng)和渦脫落受到抑制；方柱的近尾跡區(qū)域，有兩種形成機(jī)制不同的局部超聲速區(qū)。

湍流；低亞聲速；跨聲速；大渦模擬；矩形柱

0 引 言

鈍體繞流包含豐富的物理現(xiàn)象，例如在低馬赫數(shù)范圍內(nèi)存在流動(dòng)的分離和再附、剪切層的失穩(wěn)和轉(zhuǎn)捩、周期性的渦脫落以及尾跡的發(fā)展；跨聲速范圍內(nèi)會(huì)出現(xiàn)激波和渦、剪切層以及尾跡等的相互作用，尾跡中還會(huì)形成局部超聲速區(qū)。對(duì)這些典型流動(dòng)現(xiàn)象的研究不僅具有明確的科學(xué)意義，也有很高的工程價(jià)值。

作為兩種基本的鈍體幾何形狀，圓柱和矩形柱的繞流問(wèn)題得到了大量的實(shí)驗(yàn)和數(shù)值研究。對(duì)于圓柱情形，不可壓縮和可壓縮繞流兩方面的研究都比較充分，但針對(duì)矩形柱繞流的研究目前還僅限于不可壓縮或者低馬赫數(shù)的情形。Lyn等［1］利用水洞實(shí)驗(yàn)對(duì)雷諾數(shù)為21400的方柱繞流的流場(chǎng)做了全面的測(cè)量，為數(shù)值模擬方法的驗(yàn)證提供了參考。Grigoriadis等［2］采用浸沒(méi)邊界方法對(duì)方柱繞流進(jìn)行了數(shù)值模擬，得到了比較準(zhǔn)確的Strouhal數(shù)以及不同位置處的速度剖面。謝志剛等［3］采用有限體／有限元混合方法對(duì)雷諾數(shù)為22000、馬赫數(shù)為0．1的方柱繞流進(jìn)行了大渦模擬，亞格子模式采用渦粘和渦擴(kuò)散模式，并且將低馬赫數(shù)預(yù)處理方法應(yīng)用到了非定常流動(dòng)的數(shù)值模擬中，數(shù)值模擬得到的回流區(qū)長(zhǎng)度、恢復(fù)速度以及脈動(dòng)量較以往的數(shù)值模擬結(jié)果都有較大改進(jìn)。鄧小兵等［4］利用虛擬壓縮方法對(duì)三維方柱不可壓縮繞流進(jìn)行了大渦模擬，將雷諾數(shù)為22000時(shí)的主渦脫落頻率、分離區(qū)長(zhǎng)度、平均升力和平均阻力與實(shí)驗(yàn)數(shù)據(jù)做了對(duì)比，得到了較為滿意的結(jié)果。Nakaguchi等［5］實(shí)驗(yàn)研究了不同雷諾數(shù)和不同長(zhǎng)寬比的不可壓矩形柱繞流的流動(dòng)特性，結(jié)果顯示在長(zhǎng)寬比為0．6時(shí)，阻力系數(shù)達(dá)到最大值，長(zhǎng)寬比為2．8時(shí)渦脫落頻率不連續(xù)，該實(shí)驗(yàn)結(jié)果被后續(xù)的一系列相關(guān)研究［6-10］所證實(shí)。Bruno等［11］采用大渦模擬方法研究了長(zhǎng)寬比為5∶1的矩形柱繞流特性，觀察到了柱體上表面發(fā)生流動(dòng)分離和再附以及遠(yuǎn)離前緣點(diǎn)的上表面渦結(jié)構(gòu)的三維特性。

下面介紹一下目前關(guān)于圓柱跨聲速繞流的實(shí)驗(yàn)和數(shù)值研究。Murthy＆Rose［12］，Macha［13］和Rodriguez［14］實(shí)驗(yàn)研究了雷諾數(shù)為105量級(jí)的圓柱跨聲速繞流，結(jié)果表明流場(chǎng)中存在一些激波結(jié)構(gòu)，跨聲速范圍內(nèi)隨著馬赫數(shù)的增大阻力上升，馬赫數(shù)高于0．9時(shí)周期性渦脫落消失。Miserda＆Leal［15］采用有限體積法計(jì)算了二維可壓縮的Navier-Stokes方程，對(duì)來(lái)流馬赫數(shù)為0．8的圓柱繞流進(jìn)行了數(shù)值模擬，主要研究了圓柱的受力和圓柱周?chē)牧鲌?chǎng)結(jié)構(gòu)。由于實(shí)驗(yàn)時(shí)非定常的數(shù)據(jù)測(cè)量存在一定的困難，而數(shù)值研究大都求解歐拉方程，沒(méi)有流場(chǎng)結(jié)構(gòu)的細(xì)致分析，因此還有很多流動(dòng)物理機(jī)制有待探索。Xu等16采用大渦模擬方法求解了三維的可壓縮N-S方程，研究了不同馬赫數(shù)下圓柱跨聲速繞流的流動(dòng)特性，其研究表明，在跨聲速范圍存在一個(gè)臨界馬赫數(shù)（該馬赫數(shù)約為0．9），當(dāng)來(lái)流馬赫數(shù)小于臨界馬赫數(shù)時(shí)，流動(dòng)為非定常狀態(tài)；當(dāng)來(lái)流馬赫數(shù)大于臨界馬赫數(shù)時(shí)，流動(dòng)為準(zhǔn)定常狀態(tài)。通過(guò)對(duì)這兩種不同流動(dòng)狀態(tài)的流動(dòng)特性和相關(guān)物理機(jī)理的分析，研究了局部超聲速區(qū)的形成機(jī)制、湍流剪切層的演化特性和不穩(wěn)定性以及流場(chǎng)的湍流特性和渦結(jié)構(gòu)等。矩形柱與圓柱繞流的不同點(diǎn)在于矩形柱是固定點(diǎn)分離，針對(duì)矩形柱跨聲速繞流的流動(dòng)特性還有待研究。

本文中，我們采用基于非結(jié)構(gòu)網(wǎng)格的有限體／有限元混合格式和大渦模擬方法，對(duì)低亞聲速和跨聲速條件下的可壓縮矩形柱繞流開(kāi)展研究。我們主要關(guān)注矩形柱繞流的統(tǒng)計(jì)量、流動(dòng)的分離和再附以及三維渦結(jié)構(gòu)等，分別對(duì)低亞聲速和跨聲速情況下的流動(dòng)特性作了分析，在跨聲速情況下我們還考察局部超聲速的形成及其物理機(jī)理。

1 控制方程和數(shù)值方法

1.1 控制方程

我們考慮可壓縮流動(dòng)的連續(xù)方程、動(dòng)量方程，即Navier-Stokes（N-S）方程和能量方程，形式如下：

式中ui表示速度分量，p＝ρRT表示壓強(qiáng)，其中ρ為密度，R為氣體常數(shù)，T為溫度，σij＝2μ（Sij－Skkδij／3）為粘性應(yīng)力張量，其中μ為動(dòng)力粘性系數(shù)，Sij＝（?ui／?xj＋?uj／?xi）／2為變形率張量，E＝CvT＋uiui／2表示總能量，其中Cv為定容比熱容，qi＝－κ?T／?xi為熱傳導(dǎo)項(xiàng)，其中κ為熱擴(kuò)散系數(shù)。

采用密度加權(quán)的Favre濾波方法，由式（1）～式（3）可得：

式中CI為模式常數(shù)。亞格子熱傳導(dǎo)Qi采用渦擴(kuò)散模式：

式中kt為亞格子熱傳導(dǎo)系數(shù)，采用下式確定：

1.2 數(shù)值方法

本文采用基于非結(jié)構(gòu)網(wǎng)格的有限體／有限元混合方法，即采用有限體積法計(jì)算對(duì)流通量、采用有限單元法計(jì)算粘性通量，在時(shí)間方向上利用四階龍格-庫(kù)塔方法進(jìn)行推進(jìn)。為此我們需要構(gòu)造兩套重合的網(wǎng)格，即有限體網(wǎng)格和有限元網(wǎng)格，如圖1所示，有限元單元為三角形（二維）或四面體（三維），在此基礎(chǔ)上以每個(gè)節(jié)點(diǎn)為中心構(gòu)造有限體單元，兩套網(wǎng)格均覆蓋整個(gè)計(jì)算域。

我們利用有限單元的形函數(shù)，可由單元節(jié)點(diǎn)上物理量的值得到其在單元內(nèi)的分布及空間導(dǎo)數(shù)，進(jìn)而求得有限體單元邊界上的粘性通量。另一方面，對(duì)流通量的計(jì)算采用了Roe格式，其優(yōu)點(diǎn)是能夠捕捉激波并且不需要迭代。由于流動(dòng)參數(shù)值存儲(chǔ)在節(jié)點(diǎn)上（即有限體單元的中心），計(jì)算通量需要邊界上的值，因此需要在有限體單元內(nèi)進(jìn)行數(shù)據(jù)重構(gòu)。這里我們采用了MUSCL插值，使得Roe格式具有空間二階精度。為解決低亞聲速情況下Roe格式引起的剛性問(wèn)題，需要對(duì)數(shù)值方法進(jìn)行預(yù)處理。上述方法在文獻(xiàn)［3］中有詳細(xì)的解釋。此外，對(duì)于可壓縮流動(dòng)，Roe格式可能導(dǎo)致違反熵條件的非物理解，在激波附近引起數(shù)值不穩(wěn)定。因此，在馬赫數(shù)較高時(shí)，我們做了兩方面的改進(jìn)來(lái)抑制非物理解的產(chǎn)生。一是進(jìn)行熵修正，即對(duì)對(duì)流項(xiàng)系數(shù)矩陣的特征值絕對(duì)值過(guò)小的情況進(jìn)行限制；二是加通量限制器，本文中我們采用了基于非結(jié)構(gòu)網(wǎng)格的WBAP限制器［18］。

圖1 網(wǎng)格劃分示意圖Fig.1 Schematic of unstructured mesh

2 結(jié)果與分析

2.1 不同長(zhǎng)寬比矩形柱繞流的統(tǒng)計(jì)特性研究

本文的算例中，我們采用來(lái)流馬赫數(shù)Ma＝0．1和Ma＝0．75，基于矩形柱高度和來(lái)流速度的雷諾數(shù)Re＝22000，矩形柱長(zhǎng)寬比分別取1∶1，2∶1，3∶1和4∶1。當(dāng)Ma＝0．1時(shí)，計(jì)算域入口采用Steger-Warming通量條件，即利用矢通量分裂法計(jì)算從域外流入域內(nèi)的對(duì)流通量，上下邊界采用對(duì)稱邊界條件，出口采用基于特征線的無(wú)反射邊界條件，展向采用周期邊界條件，矩形柱表面采用絕熱不可滑移條件。當(dāng)Ma＝0．75時(shí)，計(jì)算域入口、出口、上下邊界都采用遠(yuǎn)場(chǎng)邊界條件，展向采用周期邊界條件，矩形柱表面采用絕熱不可滑移邊界條件。圖2為計(jì)算域、坐標(biāo)系和網(wǎng)格分布示意圖，矩形柱截面的寬度和來(lái)流速度分別用D和U∞表示。Ma＝0．1時(shí)計(jì)算域流向（x）長(zhǎng)度取為24D，矩形柱前端面距入口6．5D，計(jì)算域橫向（y）寬度H＝14D，展向（z）寬度為4D，計(jì)算時(shí)間步長(zhǎng)取0．003D／U∞。Ma＝0．75時(shí)計(jì)算域流向和橫向長(zhǎng)度均取36D，計(jì)算時(shí)間步長(zhǎng)取為0．006D／U∞，其余不變。整個(gè)計(jì)算域采用非結(jié)構(gòu)網(wǎng)格離散，柱體周?chē)木W(wǎng)格最密，尺寸約為0．01D，尾跡中次之，遠(yuǎn)離柱體的地方網(wǎng)格最疏，如圖2所示。Ma＝0．1時(shí)計(jì)算域總節(jié)點(diǎn)數(shù)為117158，總單元數(shù)為655662；Ma＝0．75時(shí)總節(jié)點(diǎn)數(shù)為274424，總單元數(shù)為1571242。隨著矩形柱長(zhǎng)寬比的變化，總節(jié)點(diǎn)數(shù)和總單元數(shù)略有不同。

圖2 計(jì)算域、坐標(biāo)系和網(wǎng)格示意圖Fig.2 Computational domain，coordinates and mesh

圖3和圖4分別顯示了矩形柱的阻力系數(shù)和升力系數(shù)脈動(dòng)均方根隨長(zhǎng)寬比的變化，并和現(xiàn)有的低亞聲速結(jié)果做了比較［10，19］。從圖中可以看出，低亞聲速和跨聲速時(shí)阻力系數(shù)都隨著長(zhǎng)寬比的變大而逐漸降低，跨聲速時(shí)升力系數(shù)脈動(dòng)均方根比低亞聲速時(shí)小很多，說(shuō)明湍流脈動(dòng)受到抑制，可以從圖11顯示的瞬時(shí)渦結(jié)構(gòu)圖像看出。

圖3 矩形柱的阻力系數(shù)隨長(zhǎng)寬比的變化Fig.3 Drag coefficient around rectangular cylinders with diferent aspect ratios

圖4 矩形柱的升力系數(shù)脈動(dòng)均方根隨長(zhǎng)寬比的變化Fig.4 Lift coefficient RMS of rectangular cylinders with different aspect ratios

流動(dòng)的特征頻率通常用Strouhal數(shù)反映，其定義為：

式中U∞無(wú)窮遠(yuǎn)處的來(lái)流速度，f為流動(dòng)的頻率，可以利用升力系數(shù)或流場(chǎng)中某點(diǎn)的壓強(qiáng)信號(hào)的頻譜得到。圖5顯示了Ma＝0．1時(shí)矩形柱繞流的Strouhal數(shù)隨長(zhǎng)寬比的變化，及其與文獻(xiàn)結(jié)果的比較［9-10，19-20］，這里流動(dòng)頻率采用升力系數(shù)的頻譜得到?？梢钥闯?，長(zhǎng)寬比從1∶1增大到2∶1時(shí)，Strouhal數(shù)降低，而長(zhǎng)寬比從2∶1變?yōu)?∶1時(shí)，Strouhal數(shù)有一個(gè)跳躍，出現(xiàn)兩個(gè)主要的流動(dòng)頻率，之后長(zhǎng)寬比再增大時(shí)，兩個(gè)Strouhal數(shù)都逐漸減小。這種變化與矩形柱表面流體的分離與再附相關(guān)。當(dāng)長(zhǎng)寬比為1∶1和2∶1時(shí)，流體從前緣分離之后在柱體的上下表面沒(méi)有產(chǎn)生再附，因而長(zhǎng)寬比越大，渦脫落頻率越??；在長(zhǎng)寬比介于2∶1和3∶1之間時(shí)，流體分離之后產(chǎn)生再附，然后再分離，此時(shí)流動(dòng)的兩個(gè)頻率分別對(duì)應(yīng)從前緣分離和從再附點(diǎn)分離的流體產(chǎn)生渦脫落的頻率，后者的值較大；之后隨著長(zhǎng)寬比的增大，再附點(diǎn)與矩形柱后緣角點(diǎn)的距離變大，渦脫落的頻率也隨之變小。

圖6表示Ma＝0．75時(shí)矩形柱的Strouhal數(shù)隨長(zhǎng)寬比的變化，這里流動(dòng)頻率分別采用升力系數(shù)和特定點(diǎn)的壓強(qiáng)信號(hào)的頻譜得到，以顯示流場(chǎng)的結(jié)構(gòu)特征。由圖可見(jiàn)，Strouhal數(shù)逐漸減小，其中采用升力系數(shù)和剪切層中的壓強(qiáng)信號(hào)得到的Strouhal數(shù)相同，而矩形柱后中心線（y＝0）上的壓強(qiáng)信號(hào)則有兩個(gè)不同的頻率，分別為升力系數(shù)對(duì)應(yīng)的Strouhal數(shù)的一半和兩倍。由于柱體上下方的渦交替脫落，使得柱體后方中心線上流體的波動(dòng)頻率是渦脫落頻率的兩倍，而另一個(gè)小的頻率表明存在長(zhǎng)周期的運(yùn)動(dòng)，需要進(jìn)一步研究。另外，當(dāng)長(zhǎng)寬比為4∶1時(shí)，只有一個(gè)大的主導(dǎo)頻率，小頻率所對(duì)應(yīng)的長(zhǎng)周期流動(dòng)結(jié)構(gòu)很弱，因此在圖中沒(méi)有顯示。

圖5 Ma＝0.1時(shí)矩形柱繞流的Strouhal數(shù)隨長(zhǎng)寬比的變化Fig.5 Strouhal around rectangular cylinders with different aspect ratios at Ma＝0.1

圖6 Ma＝0.75時(shí)矩形柱繞流的Strouhal數(shù)隨長(zhǎng)寬比的變化Fig.6 Strouhal around rectangular cylinders with diferent aspect ratios at Ma＝0.75

圖7顯示了不同長(zhǎng)寬比的矩形柱后中心線上流向平均速度沿x方向的分布，原點(diǎn)對(duì)應(yīng)矩形柱后駐點(diǎn)。圖中負(fù)的極值點(diǎn)表示回流區(qū)的強(qiáng)度，隨x增大的漸進(jìn)值反映了恢復(fù)速度的大小。可以看出，Ma＝0．1時(shí)，長(zhǎng)寬比為2∶1的回流最強(qiáng)，與前角點(diǎn)的分離流在柱體表面發(fā)生再附有關(guān)。圖中也顯示了已有的方柱繞流實(shí)驗(yàn)結(jié)果［1］作比較，可以看出回流區(qū)符合較好，而恢復(fù)速度有一定程度的高估，這也是方柱繞流問(wèn)題數(shù)值模擬中普遍存在的問(wèn)題，具體原因還有待進(jìn)一步分析。Ma＝0．75時(shí)隨著長(zhǎng)寬比的增大，矩形柱后回流區(qū)速度峰值單調(diào)變小，而出口恢復(fù)速度逐漸變大。

壓力系數(shù)Cp的定義為：

圖7 不同長(zhǎng)寬比矩形柱后中心線上流向平均速度的分布Fig.7 Mean streamwise velocity U along the center line behind the rectangular cylinder with different aspect ratios at Ma

圖8 平均壓力系數(shù)沿著不同長(zhǎng)寬比的矩形柱表面的分布Fig.8 Mean pressure distribution on the surface of the rectangular cylinder with different aspect ratios at Ma

圖9 矩形柱后中心線上脈動(dòng)速度均方根u′和v′的分布Fig.9 RMS of velocity fluctuations u′and v′on the center line behind the cylinder

圖9顯示了不同長(zhǎng)寬比矩形柱后中心線上流向速度脈動(dòng)均方根和橫向速度脈動(dòng)均方根沿x方向的分布?？梢钥闯鲭S著長(zhǎng)寬比的增加，除了低亞聲速情況下靠近柱體區(qū)域變化比較復(fù)雜，總體上流向和橫向的速度脈動(dòng)均方根都逐漸減小。需要指出的是，這里速度脈動(dòng)為瞬時(shí)值與當(dāng)?shù)貢r(shí)均值之差，沒(méi)有考慮周期性渦脫落的相位平均，所以圖9顯示的脈動(dòng)速度均方根可以反映渦的強(qiáng)度，但不能反映湍流脈動(dòng)強(qiáng)度。因此，盡管在跨聲速流動(dòng)中湍流脈動(dòng)受到抑制（如圖11所示），從圖9（a，b）和圖9（c，d）的對(duì)比中可以看到跨聲速和低亞聲速流動(dòng)的脈動(dòng)速度均方根為同一量級(jí)，結(jié)果與圖3顯示的升力系數(shù)均方根有較大差異。另外，由于柱體后方的展向渦結(jié)構(gòu)會(huì)引起較大的橫向速度，圖9（b，d）顯示橫向速度脈動(dòng)均方根v′的極大值位于2倍柱體寬度處，其值對(duì)應(yīng)于柱體后方渦脫落時(shí)離開(kāi)柱體的距離，也大致相當(dāng)于回流區(qū)的長(zhǎng)度，從后面的圖10和11可以看出。

2.2 不同長(zhǎng)寬比矩形柱繞流的流場(chǎng)特性研究

圖10 不同長(zhǎng)寬比矩形柱平均流動(dòng)的流線以及壓強(qiáng)云圖Fig.10 Streamline and pressure of mean flow around the rectangular cylinder with different aspect ratios

為了清楚顯示不同長(zhǎng)寬比矩形柱繞流是否再附，我們給出了平均流場(chǎng)的流線和壓強(qiáng)分布，如圖10所示?？梢钥闯?，Ma＝0.1時(shí)長(zhǎng)寬比為3∶1和4∶1的矩形柱有明顯的再附點(diǎn)，流體從前緣點(diǎn)分離之后在柱體的上下表面發(fā)生再附，然后再與柱體分離。長(zhǎng)寬比為1∶1和2∶1的矩形柱則沒(méi)有流動(dòng)的再附現(xiàn)象發(fā)生。Ma＝0.75時(shí)，流體從前緣分離之后并不會(huì)在柱體上下表面再附，并且隨著長(zhǎng)寬比的增大，尾跡里的低壓區(qū)壓強(qiáng)逐漸增大，從而使得阻力系數(shù)減?。ㄒ?jiàn)圖3）。

采用Q準(zhǔn)則可以對(duì)流場(chǎng)的渦結(jié)構(gòu)進(jìn)行識(shí)別和顯示，其定義如下：

圖11顯示了不同長(zhǎng)寬比矩形柱繞流的瞬時(shí)流場(chǎng)的Q等值面，Ma＝0．1時(shí)Q＝1．0，Ma＝0．75時(shí)Q＝0．1?？梢钥闯?，Ma＝0．1時(shí)矩形柱尾跡里存在復(fù)雜的三維渦結(jié)構(gòu)，而隨著長(zhǎng)寬比的增大尾跡里的渦逐漸變?nèi)?。Ma＝0．75時(shí)，渦結(jié)構(gòu)的三維特性明顯減弱，尾跡中存在強(qiáng)的展向渦形成的渦街結(jié)構(gòu)，近尾跡區(qū)連接兩展向渦的流向渦結(jié)構(gòu)。隨著長(zhǎng)寬比的增大，渦街中兩展向渦間的流向距離逐漸變大，橫向距離逐漸減小，當(dāng)長(zhǎng)寬比為4∶1時(shí)，渦基本在排列一條直線上。

圖11 不同長(zhǎng)寬比矩形柱繞流的Q等值面Fig.11 Q contours around the rectangular cylinder with different aspect ratios

盡管跨聲速時(shí)鈍體繞流的渦結(jié)構(gòu)較為簡(jiǎn)單，但近尾跡區(qū)中也存在復(fù)雜的流動(dòng)現(xiàn)象。Xu［16］指出圓柱跨聲速繞流問(wèn)題中，圓柱近尾跡區(qū)中存在兩種形成機(jī)制不同的局部超聲速區(qū)，分別為渦渦相互作用和渦與激波相互作用所形成。本文中我們考察了方柱跨聲速繞流的近尾跡區(qū)的流動(dòng)特征，發(fā)現(xiàn)也存在兩種形成機(jī)制不同的局部超聲速區(qū)。圖12給出了瞬時(shí)Ma＝1的等值面，圖中有兩塊單獨(dú)的區(qū)域，分別記為A和B，其中A處是回流區(qū)中局部超聲速區(qū)，B處為尾跡中形成的局部超聲速區(qū)。

圖12 瞬時(shí)Ma＝1的等值面Fig.12 Instantaneous Ma＝1 contours

圖13 Ma＝0.75方柱升力系數(shù)曲線表示的時(shí)間點(diǎn)位置Fig.13 Time points on the lift coefficient curve at Ma＝0.75

為了觀察局部超聲速的形成過(guò)程，我們?cè)贛a＝0.75時(shí)方柱升力系數(shù)的時(shí)間演化曲線上選取了同一個(gè)渦脫落周期內(nèi)的6個(gè)不同時(shí)間點(diǎn)a～f（見(jiàn)圖13），各時(shí)刻對(duì)應(yīng)的馬赫數(shù)等值線如圖14所示。可以看到，a時(shí)刻回流區(qū)內(nèi)馬赫數(shù)都小于1，超聲速區(qū)并沒(méi)有形成；b時(shí)刻在回流區(qū)中出現(xiàn)了M＞1的一塊小的區(qū)域，即局部超聲速區(qū)；c時(shí)刻局部超聲速區(qū)變大并且向下游移動(dòng)，需要指出的是流動(dòng)結(jié)構(gòu)的遷移速度和當(dāng)?shù)亓鲃?dòng)速度不同，盡管局部超聲速區(qū)位于回流區(qū)，仍然具有正的遷移速度；d時(shí)刻局部超聲速區(qū)繼續(xù)變大，并且要與上方剪切層中的超聲速流體匯合，同時(shí)可以看到下方超聲速剪切層從柱體脫落；e時(shí)刻前面所形成的局部超聲速區(qū)已經(jīng)完全與上方剪切層中的超聲速流體匯合，并且在方柱后不再有M＞1的區(qū)域；f時(shí)刻出現(xiàn)新的M＞1的區(qū)域，在隨后的演化中將與下方剪切層中的超聲速流體匯合，同時(shí)伴隨著上方剪切層的脫落。因此，回流區(qū)中的局部超聲速區(qū)形成后逐漸向下游遷移，然后與剪切層超聲速流體匯合在一起，并逐漸與方柱分離。這里我們可以推斷尾跡中形成的局部超聲速區(qū)是由于超聲速剪切層脫落而形成的。

圖14 一個(gè)周期內(nèi)不同時(shí)刻的馬赫數(shù)等值線分布（實(shí)線：M＞1，虛線：M＜1，時(shí)間間隔為0.3D／U∞）Fig.14 Instantaneous iso-lines of Mach number with solid lines Ma＞1 anddashed lines Ma＜1，the time interval is 0.3D／U∞

為進(jìn)一步分析回流區(qū)中的局部超聲速區(qū)的成因，我們考察了b，e兩個(gè)時(shí)刻（見(jiàn)圖13）的展向渦量和Ma＝1的等值面，如圖15所示?？梢钥闯觯瑑蓚€(gè)反向旋轉(zhuǎn)的渦的相互作用，使它們中間的流體加速，進(jìn)而形成局部超聲速區(qū)。由于對(duì)稱性，回流區(qū)中的局部超聲速區(qū)在每個(gè)渦脫落周期內(nèi)形成兩次，分別在方柱下方的渦向上運(yùn)動(dòng)和上方的渦向下運(yùn)動(dòng)時(shí)形成。

圖15 展向渦量云圖和瞬時(shí)Ma＝1的等值面Fig.15 Instantaneous spanwise vorticity and Ma＝1 contours at the instants b and e

3 結(jié) 論

本文采用有限體／有限元混合格式和大渦模擬方法求解可壓縮N-S方程，研究了不同長(zhǎng)寬比矩形柱低亞聲速和跨聲速繞流的流動(dòng)特性。低亞聲速范圍內(nèi)采用低馬赫數(shù)預(yù)處理與MUSCL插值重構(gòu)，跨聲速范圍內(nèi)采用熵修正和WBAP通量限制器，時(shí)間推進(jìn)格式采用四階龍格-庫(kù)塔方法。雷諾數(shù)為22000，來(lái)流馬赫數(shù)分別取0．1和0．75，截面長(zhǎng)寬比為1∶1、2∶1、3∶1和4∶1。研究表明：

馬赫數(shù)為0．1時(shí)，Strouhal數(shù)隨著長(zhǎng)寬比的增大先降低再增大然后再降低；長(zhǎng)寬比為3∶1和4∶1時(shí)柱體表面會(huì)有流動(dòng)的再附產(chǎn)生；渦結(jié)構(gòu)具有明顯的三維特性。

馬赫數(shù)為0．75時(shí)，Strouhal數(shù)隨著長(zhǎng)寬比的增大逐漸減??；湍流脈動(dòng)以及渦脫落受到抑制；方柱的近尾跡區(qū)域中有兩種形成機(jī)制不同的局部超聲速區(qū)。

［1］LYN D A，EIVAN S，RODI W，et al．A laser-doppler velocimetry study of the en semble-averaged characteristics of the turbulent wake of a square cylinder［J］．Journal of Eluid Mechanics，1995，304：285-319．DOI：10．1017／S0022112095004435

［2］GRIGORIADIS D G E，BARTZIS J G，GOULAS A．LES of the flow past a rectangular cylinder using the immersed boundary concept［J］．International Journal for Numerical Methods in Eluids，2003，41：615-632．DOI：10．1002／fld．458

［3］謝志剛，許春曉，崔桂香，等．方柱繞流大渦模擬［J］．計(jì)算物理，2007，24（2）：171-180．

［4］鄧小兵，張涵信，李沁．三維方柱不可壓縮繞流的大渦模擬計(jì)算［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2008，26（2）：167-172．

［5］NAKAGUCHI H，HASHIMOTO K，MUTO S．An experimental study on aerodynamicdrag of rectangular cylinders［J］．Journal of Japan Society of Aeronautic Space Science，1968，16：1-5．

［6］LAROSEG L，AUTEUIL A D．Experiments on 2D rectangular prisms at high Reynolds numbers in a pressurized wind tunnel［J］．Journal of Wind Engineering and Industrial Aerodynamics，2008，96：923-933．DOI：10．1016／j．jweia．2007．06．018

［7］SH ADARAM A，EARD M and ROSTAMY N．Experimental study of near wake flow behind a rectangular cylinder［J］．A-merican Journal of Applied Sciences，2008，5（8）：917-926．

［8］TAYLOR I，VEZZA M．Prediction of unsteady flow around square and rectangular section cylinders using adiscrete vortex method［J］．Journal of Wind Engineering and Industrial Aerod ynamics，1999，82：247-269．DOI：10．1016／S0167-6105（99）00038-0

［9］SHIMADA K and ISHIH ARA T．Application of a modifiedk-ε model to the prediction of aerodynamic characteristics of rectangular cross-section cylinders［J］．Journal of Eluid and Structures，2002，16：465-485．DOI：10．1006／jfls．2001．0433

［10］YU D，KAREEM A．Parametric study of flow around rectangular prisms using LES［J］．Journal of Wind Engineering and Industrial Aerodynamics，1998，77-78：653-662．DOI：10．1016／S0167-6105（98）00180-9

［11］BRUNO L，ERANSOS D，COSTE N，et al．3D flow around a rectangular cylinder：a computational study［J］．Journal of Wind Engineering and Industrial Aerodynamics，2010，98：263-276．DOI：10．1016／j．jweia．2009．10．005

［12］MURTHY V S，ROSE W C．Detailed measurements on a circular cylinder in cross flow［J］．AIAA Journal，1978，16（6）：549-550．DOI：10．2514／3．60930

［13］MACHA J M．Drag of circular cylinders at transonic Mach numbers［J］．Journal of Aircraft，1977，14（6）：605-607．DOI：10．2514／3．58828

［14］RODRIGUEZ O．The circular cylinder in subsonic and transonic flow［J］．AIAA Journal，1984，22（12）：1713-1718．DOI：10．2514／3．8842

［15］MISERDA R E B，LEAL R G．Numerical simulation on the unsteady aerodynamic forces over a circular cylinder in transonic flow［R］．AIAA 2006-1408．

［16］XU C Y，CHEN L W，LU X Y．Effect of Mach number on transonic flow past a circular cylinder［J］．Chinese Science Bulletin，2009，54（11）：1886-1893．DOI：10．1007／s11434-009-0325-x

［17］YOSHIZAWA A，Statistical theory for compressible turbulent shear flows，with the application to subgrid modeling［J］．Phys．Eluids，1986，29：2152-2164．DOI：10．1063／1．865552

［18］LIW A，REN Y X，LEI G D，et al．The multi-dimensional limiters for solving hyperbolic conservation laws on unstructured grids［J］．Journal of Computational Physics，2011，230：7775-7795．DOI：10．1016／j．jcp．2011．06．018

［19］TAMURA T，ITO Y．Aerodynamic characteristics and flow structures around a rectangular cylinder with a section of variousdepth／breath ratios［J］．Journal of Structural and Construction Engineering（Transactions of Architectural Institute of Japan），1996，486：153-162．

［20］NORBERG C．Elow around rectangular cylinders：pressure forces and wake frequencies［J］．Journal of wind Engineering and Industrial Aerodynamics，1993，49：187-196．DOI：10．1016／0167-6105（93）90014-E

Large eddy simulation of flow around a rectangular cylinder at low subsonic and transonic speeds

CHEN Jianping，H UANG Weixi，XU Chunxiao
（School of Aerospace，Tsinghua University，Beijing 100084，China）

Elow characteristics around a rectangular cylinder with different aspect ratios at low subsonic and transonic speeds were studied by large eddy simulation．AtRe＝22000，large eddy simulations were performed for flows around a rectangular cylinder with aspect ratios of 1∶1，2∶1，3∶1 and 4∶1 atMa＝0．1 and 0．75 respectively，to study the effects of the aspect ratio and the compressibility on the flow characteristics．AtMa＝0．1，the Strouhal numberdecreases first and then increases and finallydeceases with the increase of the aspect ratio．Elow reattachment is observed for the cases with aspect ratios of 3∶1 and 4∶1．Complex three-dimensional flow structures were observed in the wake．AtMa＝0．75，the Strouhal number and the turbulence intensitiesdecrease with the increase of the aspect ratio．It is also observed that there are two kinds of mechanisms about the formation of the local supersonic zones in the near wake region．

turbulent flow；low subsonic flow；transonic flow；LES；rectangular cylinder

V211.3

Adoi：10.7638／kqdlxxb-2013.0029

0258-1825（2014）06-0791-09

2013-03-11；

2013-12-24

國(guó)家自然科學(xué)基金（11472154，11322221，11132005）；空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放課題（SKLA20120106）

陳建平（1988-），男，碩士研究生，研究方向：湍流大渦模擬．

黃偉希，男，副教授，博士．E-mail:hwx＠tsinghua.edu.cn

陳建平，黃偉希，許春曉．低亞聲速和跨聲速矩形柱繞流的大渦模擬研究［J］．空氣動(dòng)力學(xué)學(xué)報(bào)，2014，32（6）：791-799．

10．7638／kqdlxxb-2013．0029 CHEN J P，HUANG W X，XU C X．Large eddy simulation of flow around a rectangular cylinder at low subsonic and transonic speeds［J］．ACTA Aerodynamica Sinica，2014，32（6）：791-799．