賁安慶, 竇華書(shū)

(浙江理工大學(xué)機(jī)械與自動(dòng)控制學(xué)院, 310018 杭州)

可壓縮機(jī)翼繞流的數(shù)值模擬及其穩(wěn)定性分析

賁安慶, 竇華書(shū)

(浙江理工大學(xué)機(jī)械與自動(dòng)控制學(xué)院, 310018 杭州)

運(yùn)用數(shù)值模擬的方法分別模擬了馬赫數(shù)為0.5、攻角為3和8的可壓縮的機(jī)翼繞流流動(dòng),同時(shí)研究了馬赫數(shù)為0.75、攻角為1.5,5和8的具有激波的可壓縮流動(dòng)機(jī)翼繞流,模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)符合良好。采用能量梯度方法分析了流體流動(dòng)的穩(wěn)定性,研究發(fā)現(xiàn):Spalart-Allmaras湍流模型能夠準(zhǔn)確反映出可壓縮機(jī)翼繞流流動(dòng)的流場(chǎng)特性;對(duì)于未產(chǎn)生激波的可壓縮機(jī)翼繞流,背風(fēng)的一側(cè)首先發(fā)生失穩(wěn),且在機(jī)翼前端的上緣首先發(fā)生失穩(wěn);對(duì)于具有激波的機(jī)翼繞流,激波處的能量梯度最大,首先發(fā)生失穩(wěn)。

可壓縮; 數(shù)值模擬; 機(jī)翼繞流; 能量梯度法; 穩(wěn)定性

0 引 言

繞流流動(dòng)現(xiàn)象是流體力學(xué)的經(jīng)典研究課題之一[1],也是眾多理論分析、數(shù)值模擬以及實(shí)驗(yàn)研究的對(duì)象[2-3]。它廣泛存在于自然界中,如河水繞過(guò)橋墩、風(fēng)吹過(guò)建筑物和空氣繞過(guò)飛機(jī)等;還大量出現(xiàn)在實(shí)際問(wèn)題中,如水流對(duì)渡槽槽墩、橋梁、海洋鉆井平臺(tái)支柱、海底輸運(yùn)管線、樁基碼頭等[1,4-5]。盡管對(duì)這一現(xiàn)象的研究已經(jīng)有一個(gè)多世紀(jì),但是直到現(xiàn)在仍是流體力學(xué)中的一項(xiàng)艱巨挑戰(zhàn)[6]。

眾所周知,層流繞流時(shí),流體產(chǎn)生的摩擦阻力相對(duì)于湍流繞流要小的多。在機(jī)翼繞流中,層流繞流極大地減少了能源消耗,因而機(jī)翼繞流現(xiàn)象一直受到國(guó)內(nèi)外研究學(xué)者的廣泛關(guān)注。由此產(chǎn)生的層流控制技術(shù)(LFC),其主要目的是通過(guò)各種手段調(diào)控機(jī)翼周圍的流體使其處于或者保持為層流狀態(tài)[7]。而實(shí)現(xiàn)流體控制技術(shù)的前提,是首先判斷出流體所處的狀態(tài),并找出流場(chǎng)中最容易首先發(fā)生失穩(wěn)的位置或者流體失穩(wěn)后最不穩(wěn)定的位置,從而加以調(diào)控。目前已有研究機(jī)翼穩(wěn)定性的方法主要有線性穩(wěn)定性分析方法、en等方法。然而這些方法都主要被用來(lái)預(yù)測(cè)湍流轉(zhuǎn)捩,且預(yù)測(cè)結(jié)果與實(shí)驗(yàn)結(jié)果存在較大誤差,而且當(dāng)流體處于完全湍流狀態(tài)時(shí),目前還沒(méi)有有效的數(shù)值方法判斷并找出流場(chǎng)中的最不穩(wěn)定位置。

本文運(yùn)用數(shù)值模擬的方法,分別模擬了Ma=0.5,0.75和不同攻角工況下的可壓縮氣體的繞流現(xiàn)象,并將計(jì)算結(jié)果與相關(guān)文獻(xiàn)進(jìn)行了比較,然后用能量梯度法判斷了流場(chǎng)的穩(wěn)定性。

1 能量梯度理論

近年來(lái),竇華書(shū)等[8-10]提出了一種基于牛頓力學(xué)、與N-S方程相兼容的用于分析流動(dòng)失穩(wěn)和轉(zhuǎn)捩問(wèn)題的新理論——能量梯度理論。該理論指出,粘性流體流動(dòng)的不穩(wěn)定性取決于展向獲得的能量(ΔE)與沿流線方向由于摩擦而引起的能量損失(ΔH)的相對(duì)大小(粒子在其平衡位置處的振動(dòng)示意圖如圖1所示,符號(hào)含義及具體內(nèi)容見(jiàn)文獻(xiàn)[8-10]);展向獲得的比較大的能量會(huì)放大擾動(dòng),沿流線方向損失的能量則會(huì)吸收擾動(dòng),使流動(dòng)趨于穩(wěn)定;對(duì)于任一給定的擾動(dòng),湍流轉(zhuǎn)捩與否取決于這兩者的相對(duì)大小;當(dāng)展向獲得的能量達(dá)到臨界值時(shí),層流無(wú)法平復(fù)此擾動(dòng),此時(shí)流動(dòng)有可能失穩(wěn);最終,當(dāng)展向的能量梯度與沿流線方向損失的能量梯度之比足夠大的時(shí)候,有可能引發(fā)湍流。

由上可知,流動(dòng)失穩(wěn)的準(zhǔn)則可以表示為:

(1)

(2)

2 幾何模型及數(shù)值方法

2.1 幾何模型

本文選取的計(jì)算域如圖2所示,機(jī)翼的弦長(zhǎng)為L(zhǎng)=1 m,機(jī)翼上游距機(jī)翼前端距離為15L,下游距機(jī)翼尾部為15L。

2.2 網(wǎng)格生成

本文運(yùn)用C型網(wǎng)格生成法劃分劃分機(jī)翼周圍網(wǎng)格,計(jì)算域網(wǎng)格和機(jī)翼周圍網(wǎng)格如圖3所示。其中在機(jī)翼周圍由于流體的流動(dòng)變化較為劇烈,對(duì)其周圍網(wǎng)格進(jìn)行加密處理,具體如圖3和圖4所示。最終,網(wǎng)格的最小間距為2.69×10-2,網(wǎng)格總數(shù)為11.5萬(wàn)。

2.3 控制方程

由于本文計(jì)算的工況來(lái)流馬赫數(shù)均大于0.3,要考慮氣體的可壓縮性,故其控制方程為:

(3)

(4)

(5)

由于此種工況下,流動(dòng)已為完全湍流,本文選用的湍流模型為Spalart-Allmaras湍流模型,其輸運(yùn)方程為:

(6)

2.4 邊界條件

邊界條件設(shè)置如圖2中所示,機(jī)翼采用壁面無(wú)滑移條件,機(jī)翼周圍的計(jì)算區(qū)域設(shè)為壓力遠(yuǎn)場(chǎng)邊界條件。

3 結(jié)果及討論

3.1 計(jì)算結(jié)果與文獻(xiàn)數(shù)值結(jié)果對(duì)照

分別將Ma=0.5、α=3和Ma=0.75、α=1.5兩種工況下計(jì)算得到的壓力系數(shù)與文獻(xiàn)中Jameson、Mavriplis[11]和Jameson[12]的數(shù)值結(jié)果作對(duì)比,結(jié)果如圖5和圖6所示。從圖中可以看出,本文的計(jì)算結(jié)果與文獻(xiàn)中的計(jì)算結(jié)果吻合良好。而后在此基礎(chǔ)上又分別計(jì)算了Ma=0.5、α=8和Ma=0.75、α=5或8三種工況,其壓力系數(shù)分布如圖6所示。

3.2 結(jié)果分析及討論

圖7和圖8分別給出了Ma=0.5 不同攻角工況下計(jì)算得到的流線圖、速度云圖、壓力云圖和K值云圖。在圖7中,由流線圖可知,機(jī)翼周圍流動(dòng)并沒(méi)有發(fā)生流動(dòng)分離,然而由于此時(shí)攻角α=3并不為零,流體流動(dòng)整體呈不對(duì)稱狀態(tài)。從速度云圖可知,由于攻角的原因,機(jī)翼的上表面的速度較大,而此處對(duì)應(yīng)的壓力較小。對(duì)應(yīng)圖7(d)中K值云圖可知,K值較大的地方集中分布在機(jī)翼的上表面,根據(jù)能量梯度理論可知,機(jī)翼的上表面更容易首先發(fā)生失穩(wěn)。且通過(guò)研究K值的大小分布得到,在X=0.1 m附近K值最大,根據(jù)能量梯度理論——K值最大的地方最先發(fā)生失穩(wěn)可知,在X=0.1 m附近流動(dòng)最先發(fā)生失穩(wěn),這與前人的研究成果相一致[7,13]。由圖8(a)流線圖和圖8(b)速度云圖可知,此時(shí)在機(jī)翼尾部區(qū)域初步發(fā)生了流動(dòng)的分離。對(duì)應(yīng)圖8(d)K值云圖可知,隨著攻角的增大,機(jī)翼表面的不穩(wěn)定區(qū)域向左移動(dòng),且在流動(dòng)分離的邊緣處K值較大。且研究K值的大小分布發(fā)現(xiàn),K值最大值發(fā)生在X=0.05 m處,根據(jù)能量梯度理論可判斷,此處將最先發(fā)生流動(dòng)的失穩(wěn),這與前人的研究結(jié)果[7]基本一致。

圖9—圖11給出了Ma=0.75對(duì)應(yīng)不同攻角工況下計(jì)算得到的流線圖、速度云圖、壓力云圖和K值云圖。隨著對(duì)應(yīng)的機(jī)翼的攻角不斷增大,揭示了具有激波的機(jī)翼繞流流動(dòng)的分離過(guò)程。根據(jù)上圖中的(b)速度云圖和(c)壓力云圖可知,此時(shí)在機(jī)翼的上表面已經(jīng)形成了激波。在攻角由1.5°增大到8°的過(guò)程中,流動(dòng)發(fā)生了分離,且在攻角α=8時(shí),由于流動(dòng)分離,在機(jī)翼的后方形成了分離渦。由圖9—圖11的(b)圖可以發(fā)現(xiàn),隨著攻角的增大,激波形成的區(qū)域不斷的向機(jī)翼的前端移動(dòng)。分別對(duì)比上圖中各個(gè)工況下的(b)速度云圖和(d)K值云圖可以發(fā)現(xiàn),較大的出現(xiàn)在激波位置處,這是由于激波處速度過(guò)渡較大更容易引起流動(dòng)的不穩(wěn)定性而導(dǎo)致的。而此處的K值最大,說(shuō)明此處的流動(dòng)將最先發(fā)生失穩(wěn),這與前人研究成果一致[14-15]。

4 結(jié) 論

本文運(yùn)用計(jì)算流體力學(xué)的方法成功的模擬了馬赫數(shù)Ma=0.5和Ma=0.75不同攻角工況下來(lái)流為可壓縮氣體的機(jī)翼繞流流動(dòng)。其中選用Couple算法計(jì)算了速度和壓力耦合,用有限體積法離散動(dòng)量方程,計(jì)算結(jié)果與文獻(xiàn)的數(shù)值結(jié)果吻合的很好;用能量梯度方法研究了流場(chǎng)的穩(wěn)定性,最終得出如下結(jié)論:

a) Spalart-Allmaras湍流模型適合于計(jì)算可壓縮機(jī)翼流動(dòng)。

b) 對(duì)于未產(chǎn)生激波的可壓縮機(jī)翼繞流,機(jī)翼的背風(fēng)面?zhèn)?即機(jī)翼前端的上緣,更容易首先發(fā)生失穩(wěn)。

c) 對(duì)于具有激波的機(jī)翼繞流,激波處的K值最大,最容易首先發(fā)生失穩(wěn)。

d) 能量梯度方法非常適于研究流體流動(dòng)的穩(wěn)定性問(wèn)題,可以充分的揭示其失穩(wěn)機(jī)理。

[1] 李文春. 三維氣固兩相圓柱繞流的直接數(shù)值模擬[D]. 杭州: 浙江大學(xué), 2005.

[2] Mair W A. Bluff bodies and vortex shedding[J]. Journal of Fluid Mechanics, 1971, 45: 209-224.

[3] Mujumdar A S, Douglas W J M. Vortex shedding from slender cylinders of various cross sections[J]. Journal of Fluids Engineer, 1973, 95( 3): 474-476.

[4] 楊 爍, 吳寶山. 二維圓柱繞流數(shù)值模擬[J]. 中國(guó)造船, 2008, 48(B11): 533-540.

[5] Rahman M M, Karim M M, Alim M A. Numerical investigation of unsteady flow past a circular cylinder using 2-D finite volume method[J]. Journal of Naval Architecture and Marine Engineering, 2007, 4(1): 27-42.

[6] Rajani B N, Kandasamy A, Majumdar S. Numerical simulation of laminar flow past a circular cylinder[J]. Applied Mathematical Modeling, 2009, 33: 1228-1247.

[7] Dagenhart J R, William S S. Crossflow stability and transition experiments in swept-wing flow[R/OL]. USA, Hampton: NASA Langley Research Center[1999-07-01]. http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19990056592.pdf.

[8] Dou H S. Mechanism of flow instability and transition to turbulence[J]. International Journal of Non-Linear Mechanics, 2006, 5(41): 512-517.

[9] Dou H S. Physics of flow instability and turbulent transition in shear flows[J]. International Journal of Physical Science, 2011, 6(6): 1411-1425.

[10] Dou H S, Khoo B C, Yeo K S. Instability of Taylor-Couette flow between concentric rotating cylinders[J]. International Journal of Thermal Science, 2008, 47: 1422-1435.

[11] Jameson A, Mavriplis D. Finite volume solution of the two-dimensional Euler equation on a regular triangular mesh[J]. AIAA, 1985, 21: 578-592.

[12] Jameson A. Numerical computation of transonic flows with shock waves[J]. Symposium Transsonicum, Gottingen, 1976, 56: 384-414.

[13] Johansen J. Prediction of laminar/turbulent transition in airfoil flows[J]. Journal of Aircraft, 1999, 36(4): 731-734.

[14] Wurz W, Wagner S. Experimental investigations of transition development in attached boundary layers and laminar separation bubbles[M]// K?rner H, Hilbig R. New Results in Numerical and experimental Fluid Mechanics. Vieweg+ Teubner Verlag, 1997: 413-420.

[15] Krimmelbein N, Radespiel R, Nebel C. Numerical aspects of transition prediction for three-dimensional configurations[C]//35th AIAA Fluid Dynamics Conference and Exhibit, Toronto, AIAA-2005-4764. 2005: 106.

(責(zé)任編輯： 康 鋒)

Numerical Simulation of Compressible Flow Around the Airfoiland Its Stability Analysis

BEN An-qing, DOU Hua-shu

(School of Mechanical Engineering and Automation, Zhejiang Sci-Tech University,Hangzhou 310018, China)

The compressible flows around an airfoil with the mach number of 0.5 and angles of attack of 3 and 8 were simulated with numerical simulation. Meanwhile, the compressible flows around an airfoil with shock wave, the mach number of 0.75 and angles of attack of 1.5, 5 and 8 were studied. The simulation results well accorded with experimental data. Then the stability of fluid flow was analyzed with the energy gradient theory. It is found that Spalart-Allmaras turbulence model can accurately reflect flow characteristics of fluid flow. For compressible flow around an airfoil without shock wave, instability first appears at the leeward side and the upper edge of the front end of airfoil; for flow around an airfoil with shock wave, energy gradient at the shock wave reaches the maximum and instability happens first here.

compressible; numerical simulation; flow around an airfoil; energy gradient theory; stability

1673- 3851 (2015) 05- 0675- 07

2014-12-10

國(guó)家自然科學(xué)基金項(xiàng)目(51276142,51106141);浙江理工大學(xué)研究生創(chuàng)新項(xiàng)目(11130032661215)

賁安慶(1989-),男,山東臨沂人,碩士研究生,主要從事湍流的數(shù)值模擬及穩(wěn)定性方面的研究。

竇華書(shū),E-mail:huashudou@yahoo.com

O355

A