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

1、精選優(yōu)質(zhì)文檔-傾情為你奉上計算流體力學(xué)課程大作業(yè)基于渦量-流函數(shù)法的不可壓縮方腔驅(qū)動流問題數(shù)值模擬張伊哲 航博1011、 引言和綜述2、 問題的提出，怎樣使用渦量-流函數(shù)方法建立差分格式3、 程序說明4、 計算結(jié)果和討論5、 結(jié)論1引言雖然不可壓縮流動的控制方程從形式上看更為簡單，但實際上，目前不可壓縮流動的數(shù)值方法遠(yuǎn)遠(yuǎn)不如可壓縮流動的數(shù)值方法成熟?？紤]不可壓縮流動的N-S方程：(1.1)其中是運(yùn)動粘性系數(shù)，認(rèn)為是常數(shù)。將方程組寫成無量綱的形式： (1.2)其中Re是雷諾數(shù)。從數(shù)學(xué)角度看，不可壓縮流動的控制方程中不含有密度對時間的偏導(dǎo)數(shù)項，方程表現(xiàn)出橢圓-拋物組合型的特點(diǎn)；從物理意義上看，在不

2、可壓縮流動中，壓力這一物理量的波動具有無窮大的傳播速度，它瞬間傳遍全場，以使不可壓縮條件在任何時間、任何位置滿足，這就是橢圓型方程的物理意義。這就造成不可壓縮的N-S方程不能使用比較成熟的發(fā)展型偏微分方程的數(shù)值求解理論和方法。如果將動量方程和連續(xù)性方程完全耦合求解，即使使用顯示的離散格式，也將會得到一個剛性很強(qiáng)的、龐大的稀疏線性方程組，計算量巨大，更重要的問題是不易收斂。因此，實際應(yīng)用中，通常都必須將連續(xù)方程和動量方程在一定程度上解耦。目前，求解不可壓縮流動的方法主要有渦量-流函數(shù)法，SIMPLE法及其衍生的改進(jìn)方法，有限元法，譜方法等，這些方法各有優(yōu)缺點(diǎn)。其中渦量-流函數(shù)法是解決二維不可壓縮

3、流動的有效方法。作者本學(xué)期學(xué)習(xí)了研究生計算流體課程，為了熟悉計算流體的基本方法，選擇使用渦量-流函數(shù)法計算不可壓縮方腔驅(qū)動流問題，并且對于不同雷諾數(shù)下的解進(jìn)行比較和分析，得出一些結(jié)論。本文接下來的內(nèi)容安排為：第2節(jié)提出不可壓縮方腔驅(qū)動流問題，并分析該問題怎樣使用渦量-流函數(shù)方法建立差分格式、選擇邊界條件。第3節(jié)介紹程序的結(jié)構(gòu)。第4節(jié)對于不同雷諾數(shù)下的計算結(jié)果進(jìn)行分析，并且與U.GHIA等人【1】的經(jīng)典結(jié)論進(jìn)行對比，評述本文所采用的計算方法。第五節(jié)給出結(jié)論。2問題的提出和分析2.1經(jīng)典方腔驅(qū)動流問題考慮如下圖所示的長度為1的正方形腔體，腔體上有一平板以速度U=1運(yùn)動，其它三邊為固壁條件。圖1.方

4、腔驅(qū)動流示意圖 頂蓋方腔驅(qū)動流問題是個很經(jīng)典的問題，常常用于驗證不可壓縮流動數(shù)值方法的正確性。U.GHIA等人于1982年發(fā)表的一篇文獻(xiàn)（見文獻(xiàn)【1】）計算了Re從100到的流動結(jié)果，其結(jié)果得到廣泛的認(rèn)同。2.2渦量-流函數(shù)方法簡介渦量-流函數(shù)法的基本思想是引入渦量和流函數(shù)：引入渦量，可以消去方程中的壓力項，而引入流函數(shù)，可以使連續(xù)方程自然滿足。下面對該方法進(jìn)行簡單推導(dǎo)：考慮二維問題，將式(1.2)寫成分量形式：式(1.4)對求偏導(dǎo)數(shù)減去式(1.5)對求偏導(dǎo)數(shù)，考慮到，推導(dǎo)出渦量滿足的方程為 (1.6)然后引入流函數(shù)，定義為 (1.7)可見，連續(xù)性方程(1.3)自然成立。與的關(guān)系為 (1.8)

5、式(1.6)(1.8)構(gòu)成了一個封閉的方程組，由(1.6)計算出渦量，再由(1.8)式計算出流函數(shù)，利用(1.7)式計算出速度。這個方程組的特點(diǎn)是求解速度的時候完全不用考慮壓力項。若還需要求解壓力場，則可以把式(1.4)對求偏導(dǎo)數(shù)，式式(1.5)對求偏導(dǎo)數(shù)，二者求和后整理得到關(guān)于壓力的Poisson方程 (1.9)以上推導(dǎo)出的渦量-流函數(shù)法在計算二維問題時很成功，但是三維流動的流函數(shù)沒有直觀的物理意義，無法像二維流動一樣直接定義，需要引入多個流函數(shù)，相應(yīng)解多個Poisson方程，計算量很大，并不實用。對于本文的二維問題，該方法就簡單易行。2.3建立差分格式2.3.1劃分網(wǎng)格方腔驅(qū)動流的流動區(qū)域

6、很簡單，均勻劃分為正方形的結(jié)構(gòu)網(wǎng)格即可，存儲網(wǎng)格時，x方向使用標(biāo)號i表示，y方向使用標(biāo)號j表示，x和y方向的最大網(wǎng)格點(diǎn)標(biāo)號分別為M和N。對于Re小于等于1000的情況，使用100*100網(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)心，因此不用建立壓力的差分方程。渦量的對流擴(kuò)散方程(1.6)使用FTCS格式離散得到： (1.10) 該差分格式時間方向為1階精度，空間方向為2階精度。在(1.10)中，速度分量取的是n時刻的值，已經(jīng)對方程進(jìn)行了線性化處理。流函數(shù)的Po

7、isson方程中，二階導(dǎo)數(shù)都用中心差分離散： (1.11)這種中心差分可達(dá)到二階精度。2.3.3設(shè)定邊界條件（1）速度和流函數(shù)的邊界條件由于沿著壁面是一條流線，所以流函數(shù)在邊界是常值，可以取為0；速度在邊界滿足無滑移條件。上邊界（）： ；下邊界（）： ；左邊界（）：；右邊界（）：；（2）渦量的邊界條件根據(jù)渦量的定義，在上下邊界，所以；在左右邊界，所以。； 左邊界（）：這里引入了虛擬網(wǎng)格點(diǎn)（-1,j），注意到，所以。同理可得，右邊界（）：下邊界（）：上邊界（）：注意到，所以下面考察構(gòu)造的這種渦量邊界條件的精度：比如對于下邊界，流函數(shù)的Taylor展開為則對于其他邊界的精度推導(dǎo)是類似的。所以構(gòu)造的

8、這種渦量邊界條件很有優(yōu)勢不僅形式簡單，還能有2階精度。3編程計算3.1程序的結(jié)構(gòu)程序流程圖如下圖所示：設(shè)定速度和流函數(shù)的邊界條件設(shè)定渦量、流函數(shù)和速度的初始值主程序開始Cfd_driven.f90.定義變量調(diào)用compute_node劃分網(wǎng)格，計算節(jié)點(diǎn)坐標(biāo)根據(jù)流函數(shù)的值更新渦量邊界條件，調(diào)用函數(shù)compute_omiga計算渦量場根據(jù)計算出的渦量場，調(diào)用函數(shù)compute_fai計算流函數(shù)場根據(jù)計算出的流函數(shù)場，調(diào)用函數(shù)compute_speed計算速度場是否收斂？（根據(jù)相鄰兩次速度場之1-范數(shù)）輸出結(jié)果程序結(jié)束是否圖3.流程圖3.2程序說明時間步長選擇0.001，足以滿足穩(wěn)定性條件。1、對于

9、渦量場的計算：根據(jù)差分方程該式中，速度分量取的是n時刻的值，已經(jīng)對方程進(jìn)行了線性化處理。所以直接使用顯示法，一步就可以將求出。（1.12）2、對于流函數(shù)場的計算：流函數(shù)滿足泊松方程，差分形式為：求解時要使用JACOBI迭代，由于，可以寫成其中（k+1）表示第k+1次迭代，當(dāng)時，認(rèn)為精度達(dá)到，退出迭代。將此時的作為n+1時刻的流函數(shù)場。3、速度場的求解由于，所以4結(jié)果分析圖4是Re=100,400,1000,3200,5000,7500,10000時的流線圖,每個雷諾數(shù)中，上圖是本文計算的圖，下圖是文獻(xiàn)1中的結(jié)果。與文獻(xiàn)中圖形進(jìn)行比較，流線圖很相似，計算結(jié)果是可信的。從圖4中看出，方腔驅(qū)動流中，

10、Re=100，400, 1000時，只在左右兩個下角落出現(xiàn)二次渦，Re=3200時，在靠近頂蓋（注意不是在頂蓋上）的左邊壁面又出現(xiàn)了一個二次渦，Re=7500時，右邊的下角落存在兩個二次渦, 所有二次渦的強(qiáng)度都隨著Re增大而增大。中央主渦的中心當(dāng)Re=100時偏向右上方，隨著Re增大逐漸移動到方腔的幾何中心上，當(dāng)Re=3200也就是當(dāng)左邊壁面出現(xiàn)二次渦之后，主渦的中心幾乎保持不變?？梢灶A(yù)測，當(dāng)Re接著增大，在現(xiàn)有二次渦的地方會出現(xiàn)更多的二次渦。 (a) Re = 100 (b) Re=400 (c)Re=1000 (d) Re = 3200 (e) Re=5000 (f)Re=7500 (d)

11、 Re = 1e4（右圖出自文獻(xiàn)1）圖4.不同 Re時的流線圖圖5顯示了不同Re數(shù)中，在通過方腔幾何中心的豎直線上速度u的分布情況，從圖5中看出，隨著Re增大，邊界層越來越薄，Re5000后，邊界層變薄的速率就很小了。高Re數(shù)時，方腔中部的速度分布幾乎是線性的。Re3200時，在靠近y=1處，u的分布會出現(xiàn)一個小凸起，這個現(xiàn)象也被之前的文獻(xiàn)提到過。圖5.不同Re時過方腔幾何中心的豎直線上u的分布圖6顯示了Re=100，400，1000，5000,10000時，文獻(xiàn)1中計算出的豎直中心線上u的分布和本文計算的結(jié)果，其中，Re=100，400，1000時，本文使用100*100均分網(wǎng)格，文獻(xiàn)1使用

12、129*129優(yōu)化網(wǎng)格；Re=5000和10000時，本文使用257*257均分網(wǎng)格，文獻(xiàn)1使用257*257優(yōu)化網(wǎng)格?？梢?，當(dāng)Re比較低的時候，本文的計算是很準(zhǔn)確的；當(dāng)Re增加至5000，即使將網(wǎng)格增加至257*257，仍然存在較大誤差，尤其是在靠近y=1頂蓋處，因為本文的網(wǎng)格是均分的，即使加密，效果也遠(yuǎn)遠(yuǎn)不如文獻(xiàn)1中的優(yōu)化網(wǎng)格。圖6. 不同Re時豎直中心線上u的分布對比高Re時，計算時間較長，如下表所示。Re100400100032005000750010000CPU時間（s）90166208810940110515095結(jié)論本文通過渦量-流函數(shù)法對不同雷諾數(shù)下的方腔驅(qū)動流進(jìn)行了模擬，得到

13、結(jié)論如下：對于方腔驅(qū)動流：1、隨著Re增加，二次渦逐個出現(xiàn)。2、隨著Re增加，主渦的中心逐漸移動到方腔的幾何中心，Re3200后，主渦的位置幾乎不變。3、隨著Re增加，邊界層越來越薄。高Re數(shù)時，方腔中部的速度分布幾乎是線性的，在靠近y=1處，u的分布會出現(xiàn)一個小凸起。對于本文使用的算法1、Re3200時，計算結(jié)果精良，耗時短。2、由于網(wǎng)格是均分，高Re時計算結(jié)果在邊界處有較大誤差。3、可以改進(jìn)計算流函數(shù)泊松方程時的方法，比如使用共軛梯度法，使得收斂加快，節(jié)省計算時間。參考文獻(xiàn)：1GHIA U, GHIA K N, SHIN C T . High-Re solutions for incompressible flow using the Navier-Stokes Eq