傅翔， 李高華， 王福新

(上海交通大學(xué) 航空航天學(xué)院，上海 200240)

0 引言

計(jì)算流體力學(xué)方法即指用數(shù)值手段來(lái)求解描述流場(chǎng)運(yùn)動(dòng)的N-S方程，它是獲取流場(chǎng)詳細(xì)數(shù)據(jù)的最簡(jiǎn)單直接的方法，被廣泛用于航空航天設(shè)計(jì)、汽車外形設(shè)計(jì)、化學(xué)工業(yè)處理、生物醫(yī)學(xué)工程以及仿生學(xué)研究等領(lǐng)域[1]。

計(jì)算流體力學(xué)方法起源于20世初60年代，Lax等人提出了N-S方程差分逼近的穩(wěn)定性理論[2]，為計(jì)算流體力學(xué)發(fā)展打下了良好的基礎(chǔ)。1971年，美國(guó)的弗魯姆和哈洛利用計(jì)算機(jī)第一次成功的計(jì)算出了二維長(zhǎng)方柱體的繞流，并對(duì)卡門渦街的形成做了細(xì)致的分析[3]。此后，隨著一大批先進(jìn)數(shù)值計(jì)算理論的提出[4,5]以及電子計(jì)算機(jī)的快速發(fā)展，計(jì)算流體力學(xué)方法進(jìn)入了飛速發(fā)展的時(shí)代。到20世紀(jì)80年代之后，計(jì)算流體力學(xué)方法在網(wǎng)格生成、方程離散、湍流理論以及迭代方法等方面已經(jīng)基本完善，在科學(xué)研究以及工程應(yīng)用方面的應(yīng)用也越來(lái)越多[6]。

即便計(jì)算流體力學(xué)方法在最近幾十年取得了巨大的進(jìn)步，但依然還存在著一個(gè)很大的缺點(diǎn)，就是需要耗費(fèi)非常多的計(jì)算時(shí)間。這一缺點(diǎn)在面對(duì)復(fù)雜的非定常流場(chǎng)計(jì)算時(shí)得到了放大。以計(jì)算旋翼非定常流場(chǎng)為例，通常計(jì)算一個(gè)完整的旋轉(zhuǎn)周期就需要幾個(gè)月的時(shí)間[7]，這還是在采用大規(guī)模并行處理的前提下。

為了減少計(jì)算時(shí)間，提高研究效率，本文基于速度-渦形式的N-S方程，提出了新的分區(qū)計(jì)算方法。該方法將整個(gè)流場(chǎng)計(jì)算域分為勢(shì)流域、邊界層域和N-S方程域3個(gè)區(qū)域，由于勢(shì)流域占據(jù)了絕大部分流場(chǎng)且不需要求解渦量輸運(yùn)方程，所以能夠在保證精度的情況下極大的減少計(jì)算時(shí)間。進(jìn)一步的，用該方法計(jì)算了震蕩翼型所產(chǎn)生的非定常流場(chǎng)，計(jì)算結(jié)果清晰的顯示了前緣渦和尾緣渦的脫落過(guò)程。

1 數(shù)值計(jì)算方法

1.1 基本控制方程

流體的運(yùn)動(dòng)是由N-S方程決定的，N-S方程可以通過(guò)變換寫成很多種形式，比如速度-壓力形式、速度-壓力形式以及渦量-流函數(shù)形式等。其中速度-渦形式的N-S方程可表示為[8]式(1)—式(3)。

(1)

(2)

(3)

其中v和ω分別是速度和渦量，ν是流體的動(dòng)粘性系數(shù)。

式(1)和式(2)是描述渦量和速度關(guān)系的渦量運(yùn)動(dòng)學(xué)方程。式(3)為渦量傳輸方程，它描述了流場(chǎng)中渦量的演化方式，即流場(chǎng)中的渦量變化率是由渦的傳輸、拉伸旋轉(zhuǎn)以及耗散決定的。聯(lián)立渦量傳輸方程和式(1)、式(2)即可求得流場(chǎng)中每個(gè)時(shí)刻的渦量分布。

1.2 非定常氣動(dòng)力計(jì)算方法

由于速度-渦形式的N-S方程沒(méi)有出現(xiàn)壓力項(xiàng)，所以直接采用渦量矩定理[9]進(jìn)行面積分來(lái)計(jì)算物體所受到的非定常氣動(dòng)力為式(4)。

?RSvdR

(4)

其中ρ是流體密度，r是矢徑。式(4)表明物體所受氣動(dòng)力與流場(chǎng)內(nèi)所有渦量的渦量矩的時(shí)間變化率相關(guān)，即只要知道流場(chǎng)內(nèi)的渦量分布，就可以得到相應(yīng)的氣動(dòng)力。

將式(4)分別向x和y方向投影即可得物體的阻力和升力為式(5)、(6)。

?RSvdR

(5)

(6)

1.3 計(jì)算步驟

分區(qū)計(jì)算方法的核心在于將整個(gè)流場(chǎng)分為勢(shì)流域、邊界層域和N-S方程域。在每個(gè)時(shí)刻，該方法將渦量為0的流場(chǎng)區(qū)域定義為勢(shì)流域，將物體表面的區(qū)域定義為邊界層域，將脫落渦區(qū)域設(shè)為N-S方程域。

由于勢(shì)流域內(nèi)不需要求解渦量輸運(yùn)方程而且勢(shì)流域占據(jù)了絕大部分流場(chǎng)，所以這樣的設(shè)置能夠極大的減少計(jì)算量。

分區(qū)計(jì)算方法的具體求解步驟如下：

1)在n+1時(shí)刻，更新物體邊界上的速度。

2)設(shè)置迭代步數(shù)k=1。

3)在N-S方程域內(nèi)求解渦量輸運(yùn)方程，在邊界層域內(nèi)求解邊界層方程，得到n+1時(shí)刻流體內(nèi)部的渦量預(yù)估值。

4)利用第三步求出的n+1時(shí)刻流體內(nèi)部的渦量預(yù)估值和物體邊界上的速度求出n+1時(shí)刻物體邊界上的渦量預(yù)估值。

5)更新迭代步數(shù)，k=k+1。

6)重復(fù)第三步到第五步，直到算出的n+1時(shí)刻物體邊界上的渦量預(yù)估值滿足收斂準(zhǔn)則。

7)用式(2)求出全域的速度場(chǎng)。

8)利用式(5)和式(6)計(jì)算阻力和升力。

9)回到第一步，進(jìn)入下一個(gè)時(shí)刻。

3 計(jì)算結(jié)果分析

用該方法計(jì)算了俯仰震蕩翼型繞流流場(chǎng)。翼型為NACA0012對(duì)稱翼型，雷諾數(shù)設(shè)為1×106。用商業(yè)軟件生成繞翼型的O型網(wǎng)格，網(wǎng)格外邊界的半徑為翼型長(zhǎng)度的25倍，整個(gè)網(wǎng)格及翼型附近的網(wǎng)格結(jié)構(gòu)，如圖1所示。

(a)及翼型附近網(wǎng)格(b)示意圖

圖1 整個(gè)網(wǎng)格

為了避免每個(gè)時(shí)間步都要生成網(wǎng)格，故在翼型運(yùn)動(dòng)的過(guò)程中，網(wǎng)格跟隨翼型一起運(yùn)動(dòng)。將轉(zhuǎn)動(dòng)軸取在四分之一弦長(zhǎng)處，翼型攻角的變化規(guī)律為式(7)。

(7)

其中αm是翼型的最大攻角，k是縮減頻率，c是翼型弦長(zhǎng)，U∞是無(wú)窮處來(lái)流速度。

設(shè)f為震蕩頻率，則k和f存在以下關(guān)系式為式(8)。

(8)

為了得到強(qiáng)烈的非定常流場(chǎng)，將翼型最大攻角設(shè)為αm=20°，將縮減頻率設(shè)為k=2.9。

由于在翼型震蕩的過(guò)程中，渦量主要集中在翼型背風(fēng)方向，所以將翼型背風(fēng)附近區(qū)域設(shè)置為N-S方程域，將翼型迎風(fēng)附近區(qū)域設(shè)置為邊界層域，將其它部分設(shè)為勢(shì)流域，如圖2所示。

圖2 勢(shì)流域、邊界層域及N-S方程域示意圖

分別用本文提出的分區(qū)計(jì)算方法和計(jì)算流體力學(xué)商業(yè)軟件Fluent計(jì)算一個(gè)周期T的震蕩翼型非定常流場(chǎng)，將計(jì)算所耗用時(shí)間,如表1所示。

表1 一個(gè)震蕩周期T所需要的計(jì)算時(shí)間

由表1可知，在相同的計(jì)算條件下，分區(qū)計(jì)算方法比Fluent能夠減少三分之一的計(jì)算時(shí)間，大大提高了計(jì)算效率。

分區(qū)計(jì)算方法計(jì)算出的一個(gè)震蕩周期內(nèi)翼型附近的瞬時(shí)渦量云圖，如圖3所示。

(a)t=0T(b)t=0.143T(c)t=0.286T(d)t=0.429T(e)t=0.571T(f)t=0.714T(g)t=0.857T(h)t=T

圖3 一個(gè)震蕩周期內(nèi)翼型附近的瞬時(shí)渦量云圖

從圖3中可以看出分區(qū)計(jì)算方法準(zhǔn)確的計(jì)算出了前緣渦和尾緣渦的脫落過(guò)程。在俯仰運(yùn)動(dòng)的初始時(shí)刻(t=0T)，翼型正在由0攻角位置向上抬起，一個(gè)逆時(shí)針?lè)较虻奈簿墱u逐漸在翼型尾緣的右后方形成，同時(shí)，一個(gè)順時(shí)針?lè)较虻那熬墱u也在翼型上表面形成。當(dāng)t=0.143T時(shí)，尾緣渦吸收了從尾緣脫落的渦量并逐漸變大，而前緣渦則保持了原來(lái)的大小并沿著翼型表面向下游運(yùn)動(dòng)。當(dāng)t=0.286T時(shí)，翼型的攻角幾乎已經(jīng)達(dá)到最大，尾緣渦開(kāi)始從尾緣脫落，前緣渦則運(yùn)動(dòng)到翼型尾緣位置。當(dāng)t=0.429T時(shí)，翼型開(kāi)始向下低頭，已脫落的尾緣渦順著來(lái)流向下游運(yùn)動(dòng)，而前緣渦則與新生成的順時(shí)針?lè)较虻奈簿墱u融合在一起。當(dāng)t=0.571T時(shí)，新生成的尾緣渦逐漸變大，同時(shí)一個(gè)新的逆時(shí)針?lè)较虻那熬墱u又在翼型下表面形成。當(dāng)t=0.714T時(shí)，翼型幾乎已達(dá)到負(fù)的最大攻角，新生成的尾緣渦繼續(xù)從尾緣吸收渦量并逐漸變大，而新生成的前緣渦則沿著翼型下表面向下游運(yùn)動(dòng)。當(dāng)t=0.857T時(shí)，翼型又開(kāi)始做抬頭運(yùn)動(dòng)，新生成的尾緣渦從尾緣處脫落并順著來(lái)流向下游運(yùn)動(dòng)。如此反復(fù)，就形成了以對(duì)渦形式存在的反卡門渦街。由于設(shè)置的縮減頻率比較高，流動(dòng)的不穩(wěn)定性也比較大，所以脫落的渦街并不是水平的，而是向上方偏離的。

翼型在一個(gè)震蕩周期內(nèi)所受到的升力和阻力，如圖4所示。

從圖4中可以看出，翼型所受的升力在翼型抬頭的過(guò)程中緩慢上升，至尾緣渦從尾緣脫落的時(shí)刻達(dá)到最大。在此之后，升力開(kāi)始逐漸下降，并在新生尾緣渦從尾緣脫落的時(shí)刻達(dá)到負(fù)的最大值。最后，由于翼型重新開(kāi)始抬頭運(yùn)動(dòng)，升力又逐漸上升。而翼型所受到的非定常阻力均為負(fù)值且呈現(xiàn)出正弦函數(shù)分布，即說(shuō)明了在高頻率震蕩情況下，翼型受到的是推力。這與魚類通過(guò)快速擺動(dòng)尾巴向前游動(dòng)現(xiàn)象吻合。

4 總結(jié)

本文基于渦-速度形式的N-S方程，提出一種分區(qū)計(jì)算非定常流場(chǎng)的方法，將整個(gè)計(jì)算流場(chǎng)劃分為勢(shì)流域、邊界層域以及N-S方程域三個(gè)部分。利用該方法計(jì)算了高頻俯仰震蕩翼型的非定常流場(chǎng)，發(fā)現(xiàn)該方法比商業(yè)計(jì)算軟件Fluent節(jié)省近三分之一的計(jì)算時(shí)間。此外，利用該方法計(jì)算出的非定常流場(chǎng)，對(duì)高頻震蕩翼型反卡門渦街的形成以及非定常力的演化過(guò)程進(jìn)行了詳細(xì)的研究。

圖4 一個(gè)震蕩周期內(nèi)翼型的升力和阻力(N)

[1] Versteeg H K, Malalasekera W. An introduction to computational fluid dynamics-The finite volume method[M]. Singapore: World Scientific, 1995.

[2] Reddy S C, Trefethen L N. Lax-stability of fully discrete spectral methods via stability regions and pseudo-eigenvalues [J]. Computer Methods in Applied Mechanics & Engineering, 1990, 80(1-3):147-164.

[3] Anderson, John David. Computational fluid dynamics: the basics with applications[M]. Beijing: Tsinghua Press, 2002.

[4] Blazek J, Blazek J. Computational Fluid Dynamics: Principles and Applications, Second Edition[J]. Computational Fluid Dynamics Principles & Applications, 2001, 55(2):1-4.

[5] Oberkampf W L, Trucano T G. Verification and validation in computational fluid dynamics[J]. Progress in Aerospace Sciences, 2002, 38(3):209-272.

[6] Al-Baali A G, Farid M M. Fundamentals of Computational Fluid Dynamics[J]. 2002, 55(4):33-44.

[7] Gupta R, Biswas A. Computational fluid dynamics analysis of a twisted three-bladed H-Darrieus rotor[J]. Journal of Renewable & Sustainable Energy, 2010, 2(4):691-694.

[8] Davies C, Carpenter P W. A Novel Velocity-Vorticity Formulation of the Navier-Stokes Equations with Applications to Boundary Layer Disturbance Evolution[J]. Journal of Computational Physics, 2001, 172(1):119-165.

[9] Wu J C. Theory for aerodynamic force and moment in viscous flows[J]. AIAA Journal, 1981, 19(4): 432-441.