《計算流體力學(xué)》課程大作業(yè)

——基于渦量-流函數(shù)法的不可壓縮方腔驅(qū)動流問題數(shù)值模擬

1、引言和綜述

2、問題的提出，怎樣使用渦量-流函數(shù)方法建立差分格式

3、程序說明

4、計算結(jié)果和討論

5、結(jié)論

1引言

雖然不可壓縮流動的控制方程從形式上看更為簡單，但實際匕目前不可壓縮流動的數(shù)

值方法遠遠不如可壓縮流動的數(shù)值方法成熟。

考慮不可壓縮流動的N-S方程：

vu=o

?SU1(l.l)

—+V-(UU)=f——VP+vAU

ap

其中是運動粘性系數(shù)，認為是常數(shù)。將方程組寫成無量綱的形式：

v-u=o

<5111(1.2)

—+V-(UU)=f-VP+——AU

dtRe

其中Re是雷諾數(shù)。

從數(shù)學(xué)角度看，不可壓縮流動的控制方程中不含有密度對時間的偏導(dǎo)數(shù)項，方程表現(xiàn)出

橢圓?拋物組合型的特點；從物理意義上看，在不可壓縮流動中，壓力這?物理量的波動具有

無窮大的傳播速度，它瞬間傳遍全場，以使不可壓縮條件在任何時間、任何位置滿足，這就

是橢圓型方程的物理意義，這就造成不可壓縮的N-S方程不能使用比較成熟的發(fā)展型偏微

分方程的數(shù)值求解理論和方法。

如果將動量方程和連續(xù)性方程完全耦合求解，即使使用顯示的離散格式，也將會得到一

個剛性很強的、龐大的稀疏線性方程組，計算量巨大，更重要的問題是不易收斂。因此，實

際應(yīng)用中，通常都必須將連續(xù)方程和動量方程在一定程度上解耦。

目前，求解不可壓縮流動的方法主要有渦量

-流函數(shù)法,SIMPLE法及其衍生的改進方法,

有限元法，譜方法等，這些方法各有優(yōu)缺

點。其中渦量-流函數(shù)法是解決二維不可壓

縮流動的有效方法。作者本學(xué)期學(xué)習(xí)了研究

生計算流體課程，為了熟悉計算流體的基本

方法，選擇使用渦量-流函數(shù)法計算不可壓

縮方腔驅(qū)動流問題，并且對于不同雷諾數(shù)下

的解進行比較和分析，得出一些結(jié)論。

本文接下來的內(nèi)容安排為：第2節(jié)提出不可

壓縮方腔驅(qū)動流問題，并分析該問題怎樣使

用渦量一流函數(shù)方法建立差分格式.選擇邊

界條件。第3節(jié)介紹程序的結(jié)構(gòu)。第4節(jié)對

于不同雷諾數(shù)下的計算結(jié)果進行分析，并

與U.GHIA等人【1】的經(jīng)典結(jié)論進行對比,

評述本文所采用的計算方法。第五節(jié)給出結(jié)

論。

2問題的提出和分析

2.1經(jīng)典方腔驅(qū)動流問題

考慮如下圖所示的長度為I的正方形腔體，腔體上有一平板以速度U=1運動，其它三邊

為固壁條件。

Q75

>05

0?

025050.75

圖I.方腔驅(qū)動流示意圖

頂蓋方腔驅(qū)動流問題是個很經(jīng)典的問題，常常用于驗證不可壓縮流動數(shù)值方法的正確

性。U.GHIA等人于1982年發(fā)表的一篇文獻(見文獻【1】)計算了Re從100至旭的流動結(jié)

果，其結(jié)果得到廣泛的認同。

2.2渦量?流函數(shù)方法簡介

渦量一流函數(shù)法的基本思想是引入渦量和流函數(shù)：引入渦量，可以消去方程中的壓力

項，而引入流函數(shù)，可以使連續(xù)方程自然滿足。下面對該方法進行簡單推導(dǎo)：

考慮二維問題，將式(L2)寫成分量形式：

dudv?

—+—=0(1.3)

dtdt

dudududP1(d2ud2uy

(1.4)

dtdxdydxRc(dx~dy~>

dvcvdvdP1(d2vd2vy

—+14—+V—=

22(1.5)

dtdxdydyRe(dxdy>

式(L4)對求偏導(dǎo)數(shù)減去式(1.5)對求偏導(dǎo)數(shù)，考層到，推導(dǎo)出渦量滿足的方程為

(1.6)

然后引入流函數(shù)，定義為(1.7)

可見，連續(xù)性方程(1.3)自然成立。與的關(guān)系為

叫+嗎=8(18)

dx2dy2

式(1.6)~(1.8)構(gòu)成了一個封閉的方程組，由(1.6)計算出渦量，再由(1.8)式計算出流函數(shù),

利用(1.7)式計算出速度。這個方程組的特點是求解速度的時候完全不用考慮壓力項。若還需

要求解壓力場，則可以把式(1.4)對求偏導(dǎo)數(shù)，式式(1.5)對求偏導(dǎo)數(shù)，二者求和后整理得

到關(guān)于壓力的Poisson方程

Wdydx{dyj

以上推導(dǎo)出的渦量?流函數(shù)法在計算二維問題時很成功，但

是三維流動的流函數(shù)沒有直觀的物理意義，無法像二維流動

一樣直接定義，需要引入多個流函數(shù)，相應(yīng)解多個Poisson

方程，計算量很大，并不實用。對于本文的二維問題，該方法

就簡單易行。

2.3建立差分格式

2.3.1劃分網(wǎng)格

方腔驅(qū)動流的流動區(qū)域很簡單，均勻劃分為正方形的結(jié)構(gòu)網(wǎng)格即可，存儲網(wǎng)格時,X方

向使用標號i表示,y方向使用標號j表示,X和y方向的最大網(wǎng)格點標號分別為M和N。對

于Re小于等于1000的情況，使用100*10。網(wǎng)格,Re大于1000后的情況，使用256*256網(wǎng)格。

計算域如圖2所示：

圖2.100*100的均分網(wǎng)格

2.3.2建立差分方程

由于本題關(guān)注的是方控內(nèi)部的流動狀態(tài)，對于壓力分布不關(guān)心，因此不用建立壓力的差

分方程。

渦量的對流擴散方程U.6）使用FTCS格式離散得到：

需一%?J%廠%J1.

AriJ2\xiJ2Ay

（1.10）

=J_（①-+硯I」%+「2%+%_]

—Re'（Ax）2（Ay）2

該差分格式時間方向為I階精度，空間方向為2階精度。在（1.10）中，速度分量取的是n

時刻的值，已經(jīng)對方程進行了線性化處理。

流函數(shù)的Poisson方程中，二階導(dǎo)數(shù)都用中心差分離散：

這種中心差分可達到二階精度。

2.3.3設(shè)定邊界條件

（1）速度和流函數(shù)的邊界條件

由于沿著壁面是一條流線，所以流函數(shù)在邊界是常值，可以取為o：速度在邊界滿足無

滑移條件。

上邊界（）：；

下邊界（）：；

左邊界（）：；

右邊界（）：；

（2）渦量的邊界條件

根據(jù)渦晟的定義，在上下邊界，，所以；在左右邊界，，所以O(shè)；

左邊界（）：這里引入了虛擬網(wǎng)格點（-1J）,注意到，所以。

同理可得，右邊界（）：

下邊界（）：

上邊界（）：注意到，所以

下面考察構(gòu)造的這種渦量邊界條件的精度：

比如對于下邊界，流函數(shù)的Taylor展開為

di//Ay2d2i//

=匕。+△）'歹J竽乳

小、=科「by也+叱虹-空嗎,。'）

,?"。）2!（,0）3!a/|（,0）

則曾°+*+)一2/,°+%

0（與2）

dy(/.o)"

對于其他邊界的精度推導(dǎo)是類似的。所以構(gòu)

造的這種渦量邊界條件很有優(yōu)勢一一不僅

形式簡單，還能有2階精度。

3編程計算

3.1程序的結(jié)構(gòu)

程序流程圖如下圖所示:

圖3.流程圖

3.2程序說明

時間步長選擇o.(x)i,足以滿足穩(wěn)定性條件。

1、對于渦量場的計算：根據(jù)差分方程

端-%”"+1,廣“工"%+1-①。T1MilJ-2%+d]J峭+1-2%+%]

Z八，2Ax32AyRe(Ar)2(Ay)

該式中，速度分量取的是n時刻的值，已經(jīng)對方程進行了線性化處理。所以直接使用顯示法，

一步就可以將求出。

八，?+】八”,

%=電”+

Aj1J噥J-2dj+球T.j吠川-24j+%T(*J—以TJ“%—%T]

Rc(Ar)2(Ay)(“J2ArM2Ay

(1.12)

2.對于流函數(shù)場的計算:

流函數(shù)滿足泊松方程，差分形式為：

廠2匕；+仁：心-2K:+*工_

K?+亍=%

求解時要使用JACOBI迭代，由于，可以寫成

《產(chǎn)川=；(哨黃十屋：片+哨心+小，)-碗；

其中(k+l)表示第k+1次迭代，當時，認為精度達到，退出迭代。將此時的作為n+l時

刻的流函數(shù)場。

3.速度場的求解

由于，所以

2Ay'2Ar

4結(jié)果分析

圖4是Re=100,400,1000.3200,5000,7500,10000時的流線圖，每個雷諾數(shù)中，上圖是本文

計算的圖，下圖是文獻1中的結(jié)果。與文獻中圖形進行比較，流線圖很相似，計算結(jié)果是可

信的。

從圖4中看出,方腔驅(qū)動流中,Rc=100,400,10(H)時、只在左右兩個下角落出現(xiàn)二次渦,

Re=3200時，在靠近頂蓋(注意不是在頂蓋上)的左邊壁面又出現(xiàn)了一個二次渦，Re=75OO

時，右邊的F角落存在兩個二次渦，所有二次渦的強度都隨著Rc增大而增大。中央主渦的

中心當Re=100時偏向右上方，隨著Re增大逐漸移動到方腔的幾何中心上，當Re=3200也就

是當左邊壁面出現(xiàn)二次渦之后，主渦的中心幾乎保持不變?？梢灶A(yù)測，當Re接著增大，在現(xiàn)

有二次渦的地方會出現(xiàn)更多的二次渦。

(a)Rc=10()(b)Rc=4(X)(c)Re=l()()()

(d)Re=320()(e)Re=50(X)(DRe=7500

(d)Re=Ie4(右圖出自文獻I)

圖4.不同Rc時的流線圖

圖5顯示了不同Re數(shù)中，在通過方腔幾何中心的豎直線上速度u的分布情況，從圖5

中看出，隨著Rc增大，邊界層越來越薄,Re>5()00后，邊界層變薄的速率就很小了。高Rc數(shù)

時，方腔中部的速度分布幾乎是線性的。Re>3200時，在靠近y=l處,u的分布會山現(xiàn)一個小

凸起，這個現(xiàn)象也被之前的文獻提到過。

圖5.不同Re時過方腔幾何中心的豎直線上u的分布

圖6顯示了Re=100,400,1000,500040000時，文獻1中計算出的豎直中心線上u的分

布和本文計算的結(jié)果,其中，Re=100,400,1000時，本文使用100*100均分網(wǎng)格，文獻1使用

129*129優(yōu)化網(wǎng)格：Re=5000和10000時，本文使用257*257均分網(wǎng)格,文獻1使用257*257

優(yōu)化網(wǎng)格?？梢?，當Re比較低的時候，本文的計算是很準確的；當Re增加至5000,即使將

網(wǎng)格增加至257*257,仍然存在較大誤差，尤其是在靠近y=l頂蓋處，因為本文的網(wǎng)格是均

分的，即使加密，效果也遠遠不如文獻I中的優(yōu)化網(wǎng)格。

u

高Re時，140010003200500075010

計算時間較00000(1

長，如下表所

CPU時9166208810