非定常流場的分區(qū)計算方法

傅翔， 李高華， 王福新

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

0 引言

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

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

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

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

1 數(shù)值計算方法

1.1 基本控制方程

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

(1)

(2)

(3)

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

式(1)和式(2)是描述渦量和速度關系的渦量運動學方程。式(3)為渦量傳輸方程，它描述了流場中渦量的演化方式，即流場中的渦量變化率是由渦的傳輸、拉伸旋轉以及耗散決定的。聯(lián)立渦量傳輸方程和式(1)、式(2)即可求得流場中每個時刻的渦量分布。

1.2 非定常氣動力計算方法

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

?RSvdR

(4)

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

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

?RSvdR

(5)

(6)

1.3 計算步驟

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

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

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

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

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

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

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

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

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

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

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

9)回到第一步，進入下一個時刻。

3 計算結果分析

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

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

圖1 整個網(wǎng)格

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

(7)

其中αm是翼型的最大攻角，k是縮減頻率，c是翼型弦長，U∞是無窮處來流速度。

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

(8)

為了得到強烈的非定常流場，將翼型最大攻角設為αm=20°，將縮減頻率設為k=2.9。

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

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

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

表1 一個震蕩周期T所需要的計算時間

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

分區(qū)計算方法計算出的一個震蕩周期內(nèi)翼型附近的瞬時渦量云圖，如圖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 一個震蕩周期內(nèi)翼型附近的瞬時渦量云圖

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

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

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

4 總結

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

圖4 一個震蕩周期內(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.