平面普修流之非線性不穩定

(1)

國立交通大學

機械工程學系

碩士論文

平面普修流之非線性不穩定

Nonlinear Instability of Plane Poiseuille Flow

研究生：洪英棋

指導教授：楊文美博士

(2)

平面普修流之非線性不穩定

Nonlinear Instability of Plane Poiseuille Flow

研究生：洪英棋 Student：Ying-Chi Hung 指導教授：楊文美 Advisor：Wen-Mei Yang 國立交通大學機械工程學系碩士論文 A Thesis

Submitted to Instuitute of Mechanical Engineering College of Engineering

National Chiao Tung University in Partial Fulfillment of The Requirments

for the Degree of Master of Science

in

Mechanical Engineering July 2005

Hsinchu, Taiwan, Republic in China

(3)

摘要

普修流是流體力學中最基本也是最重要的問題，當學者們發展出新的實驗方法或是數值方法，普修流為驗證新方法正確性的依據。本文以數值方法探討普修流在非線性擾動作用下的行為，在數值方法上，吾人先以複數形式的 Fourier 級 數對 x 方向展開，在 y 方向則選用具有高精確度及收斂快速 特性的契比希夫多項式 (Chebyshev polynomials) 來對空間座標作雙重級數展開，並以 Collocation Method 與 Galerkin Approximation 將方程式轉換為代數方程式，並將代數方程式中時間非線性項利用 Adam- Bashforth Method 作離散化，而線性項則採用 Crank-Nicolson Method 處理。本文中以改變雷諾數的大小來觀察流場隨時間的變化，結果發現，流場在低雷諾數時擾動會隨著時間呈現快速收斂，提高雷諾數，流場擾動行為變為劇烈，收斂至穩定值的時間也隨之增加。而在臨界雷諾數 Re =5772.22時，擾動不再c 隨著時間收斂至定值，呈現單一週期性的狀態。當超過臨界雷諾數在 Re=6000~10000範圍內，流場會隨著雷諾數的增加使擾動的行為變為劇烈，此時會產生第二個不相干的頻率，呈現擬週期性運動。隨著雷諾數提升至 11000時，流場擾動的更加劇烈，且呈現出微弱的混沌現象，當雷諾數為 12900時，流場擾動混亂無法判讀，此時流場混沌現象更加劇烈。

(4)

誌謝

本論文可以順利完成，在此感謝指導教授楊文美博士的悉心指導與督促，特此致上感激與謝意。承蒙口試委員胡毓仁博士、崔燕勇博士與傅武雄博士於口試期間，給予論文指正與寶貴意見，使得本論文更為嚴謹與完善，在此亦表達由衷的感謝。本研究得以順利完成，感謝同窗好友友約與弘仁，以及實驗室學長勝文、豪傑和學弟啟豪、秋傑於課業上及生活上的互相提攜，使得研究生涯更為充實。最後，僅以本論文成果獻給我敬愛的父母、親愛的家人以及默默支持我的女友靜宜，感謝你們的關心與鼓勵。

(5)

iii

目錄

摘要………..…………i 誌謝……….ii 目錄………iii 圖目錄……….v 符號說明………...vii 第一章緒論……….1 1. 1 文獻回顧……….………..1 1. 2 研究目的……..………...3 第二章數學模式……….4 2. 1 統御方程式………...4 2. 2 無因次化………...5 2. 3 基態解………...5 2. 4 擾動方程式之建立………...6 第三章數值方法……….9 3. 1 頻譜法………...9 3.1.1 雙重級數展開………..………...9 3.1.2 離散方程式………...………10 3.1.3 代數方程式………...………10 3. 2 數值解的分析法……...……….……….13 3.2.1 時間級數與FFT頻譜分析……….…………. ………13 3.2.2 相平面圖………...13 3.2.3 Poincare映射圖……….14 第四章結果與討論………...15

(6)

4. 1 線性擾動方程式及準確度分析……….15 4. 2 非線性擾動方程式及高於臨界雷諾數之流場行為………….16 4. 3 由時間級數與 FFT 頻譜圖探討流場行為……….16 4. 4 由相平面圖與 Poincare 映射圖探討流場的行為…...17 第五章結論………...20 參考文獻………...21

(7)

圖目錄

圖一普修流模型………24 圖二當 M=3，dt=0.5，α＝1.020545 時不同 N 所得到的臨 界雷諾數………25 圖三 Rec對 M&N，dt=0.5，α＝1.020545 時所得到的臨界 雷諾數………26 圖四當 dt=0.5，α＝1.020545 時，速度隨時間變化圖 Re= (a)1000 (b) 4000(c) 5772.22 (d) 6000………27 圖五當 N=31，dt=0.5，α＝1.020545 時不同 M 所得到的 臨界雷諾數………...29 圖六加入非線性項後 dt=0.5，α＝1.020545 時，速度隨時 間變化圖，Re=(a)1000 (b) 4000...30 圖七加入非線性項後,雷諾數 Re=5772.22，時間 dt=0.2，波數α ＝1.020545 的(a)時間級數圖(b)頻譜圖...31 圖八加入非線性項後,雷諾數 Re=6000，時間 dt=0.2，波數 α＝1.020545 的(a)時間級數圖(b)頻譜圖…...32 圖九加入非線性項後,雷諾數 Re=7000，時間 dt=0.2，波數 α＝1.020545 的(a)時間級數圖(b)頻譜圖…...33 圖十加入非線性項後,雷諾數 Re=8000，時間 dt=0.2，波數 α＝1.020545 的(a)時間級數圖(b)頻譜圖…...34 圖十一加入非線性項後,雷諾數 Re=9000，時間 dt=0.1，波數 α＝1.020545 的(a)時間級數圖(b)頻譜圖…...35 圖十二加入非線性項後,雷諾數 Re=10000，時間 dt=0.1，波數 α＝1.020545 的(a)時間級數圖(b)頻譜圖…...36 圖十三加入非線性項後,雷諾數 Re=11000，時間 dt=0.1，波數 α＝1.020545 的(a)時間級數圖(b)頻譜圖…...37

(8)

圖十四加入非線性項後,雷諾數 Re=12900，時間 dt=0.1，波數 α＝1.020545 的(a)時間級數圖(b)頻譜圖…...38 圖十五雷諾數 Re=5772.22，時間 dt=0.2，波數α ＝1.020545 的(a)相平面圖(b)Poincare 映射圖...39 圖十六雷諾數 Re=6000，時間 dt=0.2，波數α ＝1.020545 的 (a)相平面圖(b)Poincare 映射圖...40 圖十七雷諾數 Re=7000，時間 dt=0.2，波數α ＝1.020545 的 (a)相平面圖(b)Poincare 映射圖...41 圖十八雷諾數 Re=8000，時間 dt=0.2，波數α ＝1.020545 的 (a)相平面圖(b)Poincare 映射圖...42 圖十九雷諾數 Re=9000，時間 dt=0.1，波數α ＝1.020545 的 (a)相平面圖(b)Poincare 映射圖...43 圖二十雷諾數 Re=10000，時間 dt=0.1，波數α ＝1.020545 的(a)相平面圖(b)Poincare 映射圖...44 圖二十一雷諾數 Re=11000，時間 dt=0.1，波數α＝1.020545 的(a)相平面圖(b)Poincare 映射圖...45 圖二十二雷諾數 Re=12900，時間 dt=0.1，波數α ＝1.020545 的(a)相平面圖(b)Poincare 映射圖...46

(9)

vii

符號說明

符號

mn a 時間函數 c 常數 f 頻率 H 平行板之間距 N M , 流線函數展開項之項數 p 壓力 Re 雷諾數 t 時間 n T n階chebyshev多項式 v u, 速度分量 y x, 平面座標 φ y方向展開函數 ψ 流線函數(stream function) µ 黏滯係數 ρ 密度 α 波數(wave number) ∆ 2 2 2 2 x y ∂ ₊ ∂ ∂ ∂ 2 ∆ 4 4 4 4 2 2 4 2 x x y y ∂ ₊ ∂ ₊ ∂ ∂ ∂ ∂ ∂

上標說明

- 基態解 ´ 擾動量

(10)

k 第k個時間間隔

下標說明

(11)

第一章

緒論

1.1 文獻回顧

普修流(Poiseuille flow)是一固定壓力梯度使得流體流動之場。普修流是流體力學中最基本也是最重要的問題之一，當學者們發展出新的實驗方法或是數值方法，普修流這樣基本的問題成為了驗證新方法正確性的依據。普修流流場中，高低不同的雷諾數產生不同的物理現象，當雷諾數超過臨界雷諾數時就發生紊流(turbulence)。若新的實驗或數值方法的結果能夠與物理現象吻合，表示方法為可行，學者們再以同樣方法為基礎將研究擴展到更複雜的問題上，由此可見其重要性。早在 19 世紀末 Osborne Reynolds 以圓管做了普修流的實驗，他發現層流在雷諾數為 13000 時發生不穩定，而後人發現在實驗物表面非常粗糙時臨界(critical)雷諾數約為 2000，而若表面為平滑則臨界雷諾數可達 40000。隨後 W. Orr & A. Sommerfeld(1907)針對平行通道 (channel)內之普修流提出一計算流體穩定性之 Orr-Sommerfeld 方程式。

Sexl(1927)將在圓管普修流穩定性的問題理論化，他針對軸對稱之圓管流推導出相關之方程式。之後 Sexl 和 Lessen,Fox, Bhat 以及 Liu(1964)探討了較為複雜的非對稱管流的問題。在實驗以及理論對照下，他們發現了普修流在受到極微小的擾動下呈穩定，但在一般的擾動下則呈現不穩定的狀態。時至 1950-1960 年代，計算機科技初步發展，學者在試著以各種數值方法探討線性普修流的問題，針對 Orr-Sommerfeld 方程式找出其發生不穩定的臨界雷諾數。Shen(1955)使用由 Lin(1954)所發表的近

(12)

似分析 (asymsptotic analysis) 法得到 Re_c =5360 ， α ＝ 1.05 ， Lin並且建立了中性曲線(neutral curve)；Thomas(1953)率先以有限

差分法得到Re_c =5780，α ＝1.026；Betchov & Criminale(1967)也以

有限元素法得到Rec=5767，α ＝1.02。Grosch & Salwen(1968)利用

Chandrasekhar -Reid(1961)正交函數得到Rec=5750，α ＝1.025。 Orszag(1971)以多項的 Chebyshev 展開式搭配 QR 矩陣特徵值演算法得到該式臨界雷諾數發生於Re_c=5772.22，α ＝1.02056。他以 Orr-Sommerfeld 方程式中 Im(λ)之正負值判斷解的穩定與否。當 Re =5772.22 時，Im(λ)全為負值，雷諾數超過此值 Im(λ)會有正值出現。Orszag 的結果具重大的參考價值，證明 Chebyshev 展開式的精確性及快速收斂的特性，並可節省大量 CPU time，適合用來求解流體運動的問題。

對於非線性 Navier Stokes 方程式的探討最先為 Meksyn & Stuart(1951)，而 Dowell(1969)接續了此研究並更系統化地探討了非線性項加入對方程式的影響。Dowell 對流場的水平方向採 Fourier 級數展開，而垂直方向採用符合邊界條件的兩組餘弦函數(cosine)展開，並使用了多項展開模式(spectral method)試著精確地描述流場現象。在不含非線性項的情形，Re=80,000，α ＝2 的條件下，Dowell 使用了四十項 y 方向奇函數展開得到與 Thomas 使用 100 個有限差分點(finite -difference points)相同的特徵值值，說明 spectral method 的可行性。在包含非線性項的同樣條件，y 方向取十六項展開的情形，發現第一階展開項係數對時間解呈 limit cycle，而第二階係數則收斂至一固定值。Dowell 另外發現垂直方向若有奇函數的展開能引發非線性的擾動現象，若僅使用偶函數則無法引發此現象。Dowell 以當時計算機速度無法得到非常理想的結果，並建議在先由線性方程式計算出合適的 eigenfunction 再代入完整的方程式中，能夠減少使用的展開項數。

(13)

Fortin 與 Jardak(1994)對於普修流的穩定性分析做了一個完整的統合。他們將問題分為一維及二維、線性及非線性，首先以數值方法解出線性的 Orr-Sommerfeld 方程式，得到臨界雷諾數發生於α ＝ 1.020545，Rec=5772.22，與 Orszag 的結果吻合，而後採用疏密不同的網格，利用有限元素法直接解出非線性的擾動方程式，並且將完整的流場繪出。結果顯示在低雷諾數(Re=525)時流場呈現穩定的流動，擾動在一定的時間內達到收斂：在雷諾數提高到臨界值時，時間級數圖顯示擾動並不隨時間收斂，利用 FFT(Fast Fourier Transform)及相平面圖(Phase plane)得知流動行為僅有一組基頻與其倍數頻率的組合，為單一週期性的運動，稱為週期性流場。當雷諾數提高至 6000 時，出現第二組頻率，流場轉變為擬週期性(quasi-periodic)的運動。當雷諾數為 8000 至 10000 時流場依然為擬週期性的運動，不同之處僅在於流動行為更加劇烈。雷諾數提高至 11000 時，流場呈現微弱混沌現象(weakly chaos)，即流場已產生兩組以上之頻率。當雷諾數達到 12900 時，流場已產生混沌現象，無法由 FFT 及相平面圖解讀其現象。Fortin 的結果與前人吻合，並在超過臨界雷諾數後有完整的分析，為近年研究普修流非線性現象的重要指標。

1.2 研究目的

普修流的穩定性問題早在 1970 年代左右已深入地被探討，學者們從實驗、理論、數值方法切入此問題在不同雷諾數下，是否符合實際物理現象。而吾人在本文欲嘗試先以 Chebyshev polynomials 及複數形式的 Fourier 級數對空間座標作雙重級數展開來模擬平面流場在低於以及超過臨界雷諾數的物理現象。

(14)

第二章

數學模式

本文中研究的普修流物理模式如圖一所示。兩平行板距為2H ， 其間充滿黏性流體，定義座標在水平方向為 x，並假設 x 方向為無限 長，垂直方向為 y，水平方向存在一固定壓力梯度(−∂P ∂x=c)。吾人以此為基礎，建立本章之數學模式。

2.1 統御方程式

考慮不可壓縮之牛頓流體，其直角座標(x,y)的統御方程式如下： 連續方程式(continuity equation)： 0 u v x y ∂ ₊ ∂ ₌ ∂ ∂ ( 2 . 1 ) 動量方程式(momentum equations)： 2 2 2 2 ( u u u v u) P ( u u) t x y x x y

ρ

∂ + ∂ + ∂ = −∂ +

µ

∂ +∂ ∂ ∂ ∂ ∂ ∂ ∂ (2.2a) 2 2 2 2 ( v u v v v) P ( v v) t x y y x y

ρ

∂ + ∂ + ∂ = −∂ +

µ

∂ +∂ ∂ ∂ ∂ ∂ ∂ ∂ (2.2b) 上式中P為壓力項， µ為黏滯係數(viscosity)，ρ為流體密度(density)。邊界條件(boundary conditions) ： 0 ) , , ( ) , , (x ±H t =v x ±H t = u (2.3)

(15)

2.2 無因次化

定義無因次參數如下： x* x H = * y y H = * u u u_∞ = * v v u_∞ = * 2 P P u

ρ

∞ = * u t t H ∞ = Re u H v ∞ = 上次中u∞為流道中央處最大速度。將統御方程式無因次化，並除去上標後，得以下之形式：連續方程式(continuity equation)： 0 u v x y ∂ ₊∂ ₌ ∂ ∂ (2.4) 動量方程式(momentun equations)： 2 2 2 2 1 ( ) Re u u u P u u u v t x y x x y ∂ ₊ ∂ ₊ ∂ _{= −}∂ ₊ ∂ ₊∂ ∂ ∂ ∂ ∂ ∂ ∂ (2.5a) 2 2 2 2 1 ( ) Re v v v P v v u v t x y y x y ∂ ∂ ∂ ∂ ∂ ∂ + + = − + + ∂ ∂ ∂ ∂ ∂ ∂ (2.5b) 邊界條件(boundary conditions) ： 0 ) , 1 , ( ) , 1 , (x ± t =v x ± t = u (2.6)

2.3 基態解

考慮流體在基態(basic state)為一維穩態運動，u =u y( )、v=0，基態速度及壓力梯度可以下式表示： 2 ( ) 1 u y = − y (2.7)

(16)

2 ( ) Re P x = − x_(2.8)

2.4 擾動方程式之建立

本文欲探討雷諾數對於一維基態造成不穩定的問題，在基態加上擾動量如下所示： ( ) ( , , ) u=u y +u x y t′ (2.9a) 0 ( , , ) v= +v x y t′ (2.9b) ( ) ( , , ) P = P x +P x y t′ (2.9c) 並代入統御方程式中，可得到擾動方程式。連續方程式(continuity equation)： u v 0 x y ′ ′ ∂ ₊∂ ₌ ∂ ∂ (2.10) 動量方程式(momentun equations)： y u v dy y u d v x u u x u u t u ∂ ′ ∂ ′ + ′ + ∂ ′ ∂ ′ + ∂ ′ ∂ + ∂ ′ ∂ ( ) ) ( Re 1 2 2 2 2 y u x u x P ∂ ′ ∂ + ∂ ′ ∂ + ∂ ′ ∂ − = _(2.11a) y v v x v u x v u t v ∂ ′ ∂ ′ + ∂ ′ ∂ ′ + ∂ ′ ∂ + ∂ ′ ∂ ) ( Re 1 2 2 2 2 y v x v y P ∂ ′ ∂ + ∂ ′ ∂ + ∂ ′ ∂ − = _(2.11b) 方程式(2.11)為描述擾動量之動量方程式，其邊界條件如下：

(17)

u′(x,±1,t)=v′(x,±1,t)=0 (2.12) 引用流線函數(stream function)

ψ

( , , )x y t 和速度分量之關係，即 y u ∂ ∂ = ′

ψ

_、 x v ∂ ∂ − = ′

ψ

_{，則(2.11)可表示為：} 2 2 2 2 2 u u t y x y y x y x y x y

ψ

∂ ₊ ∂ ₊ ∂ ∂ ₋ ∂ ∂ ₋ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ 3 3 2 3

1 Re

P

x

x y

y

ψ

′

⎛

⎞

∂

= −

+

_⎜

+

_⎟

∂

_⎝

∂ ∂

∂

_⎠

(2.13a) 2 2 2 2 2 2 u t x x y x x y x

ψ

∂ ∂ ∂ ∂ ∂ ∂ − − − + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ 3 3 3 2

1 Re

P

y

x

y x

ψ

′

⎛

⎞

∂

= −

−

_⎜

+

_⎟

∂

_⎝

∂

∂ ∂

_⎠

(2.13b) 將(2.13a)對y偏微分減去(2.13b)對x偏微分以消去 P′ ，可得以下的方程式：

x

u

x

u

t

∂

″

−

∆

∂

+

∆

∂

_ψ

ψ

2

Re

1 ∆

=

∆

∂

−

∆

∂

+

y

x

y

(2.14) 上式中 2 2

dy

u

d

u

″

=

_， 2 2 2 2 x y ∂ ∂ ∆ = + ∂ ∂ ， 4 4 4 2 4 2 2 4 2 x x y y ∂ ∂ ∂ ∆ = + + ∂ ∂ ∂ ∂ 其邊界條件為:

0 )

,

1 ,

(

)

,

1 ,

(

)

,

1 ,

(

=

∂

±

∂

=

∂

±

∂

=

±

y

t

x

t

x

t

x

ψ

_(2.15) (2.14)式中若消去非線性項即為 Orr-Sommerfeld 方程式:

(18)

ψ

2

Re

1 ∆

=

∂

″

−

∆

∂

+

∆

∂

x

u

x

u

t

(2.16)

對於 Orr-Sommerfeld 方程式，Orszag(1971)以多項的 Chebyshev 展開式搭配 QR 矩陣特徵值演算法得到該式臨界雷諾數發生於 c Re =5772.22，α＝1.02056，而 Fortin 與 Jardak(1994)對於普修流的穩定性分析做了一個完整的統合，以數值方法解出線性的 Orr-Sommerfeld 方程式，得到臨界雷諾數發生於α ＝1.020545，Rec =5772.22，與 Orszag 的結果吻合。

(19)

第三章

數值方法

在前一章方程式 (2.14) 中除了前三項和等號右邊項為線性，其餘二項皆為非線性，無法直接求得其解析解，所以使用數值方法來求解此類問題。

3.1 頻譜法 (Spectral method)

3.1.1 雙重級數展開

在數值求解中吾人將流線函數以雙重級數展開，來模擬 流場的變化，由於流場在 x方向以週期性出現，因此以複數形 式的 Fourier級數展開，在 y方向則選用具有高精確度及收斂快 速特性的契比希夫多項式 (Chebyshev polynomials)，而展開之係數為時間的函數，則可將流線函數表示為： mx i m n n mn t y e a t y x

φ

α

ψ

( , , ) =

∑∑

( ) ( ) _{( 3 . 1 )} 上式中

φ

n是滿足速度邊界條件的基底函數，表示如下：

⎪

⎩

⎪⎪

⎨

⎧

−

+

−

+

−

=

n

odd

T

n

T

n

T

n

even

T

n

T

n

T

y

n n n 1 2 2 2 3 2 2 0 2 2 2

1

3

1

3

1 )

1

4 (

4 )

(

φ

) ( y T_n 為第一型 n階契比希夫多項式 (Chebyshev polynomials) ，其定義如下：

(20)

,...

2 ,

1 ,

0 ,

)

cos

cos(

)

(

y

=

n

⋅

−1

y

n

=

T

_n 其中 y之定義域為 [-1,1]。

3.1.2 離散方程式

在將方程式 (2.14)加以離散化 (discretization)過程中，線性項部份，除時間項採用 I mp l i ci t M et h o d外，其它均使用 Crank-Nicolson Method處理，即 2 1 k k f f f = + + ;非線性項則採用 Adam-Bashforth Method處理，即 2 3 − −1 = f k f k f 。則未知數之最新一次的值，可由前兩次的結果求得。經推導出離散化的流線方程式，如下所示： 1 2 1 1 1

Re

2

+ + + +

₋

∆

_∆

∂

″

∆

−

∆

∂

∆

+

∆

k k k k

t

x

u

t

x

u

t

_ψ

ψ

_ψ

ψ

k k k k

t

x

u

t

x

u

t

_ψ

ψ

_ψ

ψ

2

Re

2

2 ∆

∆

+

∂

″

∆

+

∆

∂

∆

−

∆

=

1 1

2

3

−

_∆

−

∂

∆

+

∆

∂

∆

−

k k k

x

y

t

x

y

t

ψ

_ψ

ψ

_ψ

1 1

2

3

−

_∆

−

∂

∆

−

∆

∂

∆

+

k k k k

y

x

t

y

x

t

ψ

_ψ

ψ

_ψ

(3.2)

3.1.3 代數方程式

將 (3.1)式模擬流線函數的多項式代入 (3.2)當中，在 x方向 以 Galerkin method處理 ,y方向以定置法 (Collocation methed)

處理而定置點取為 i N N i y_i , 1, 2, 3,... 1 cos = + =

π

_{，如此可將}

(21)

統御方程式轉化為代數方程式，以矩陣表示如下： 1 , 1 − +

₌

kk k

F

AX

( 3 . 3 ) 其中，矩陣 A的大小為 2 ) (M ×N ,矩陣 X及 F的大小為(M×N)，其組成如下：

⎥

⎦

⎤

⎢

⎣

⎡

mm

A

L

M

O

M

L

0

22 11 1 2 1 +

⎥

⎦

⎤

⎢

⎣

⎡

k m

X

M

1 2 1 2 1 2 1 −

⎥

⎦

⎤

⎢

⎣

⎡

+

⎥

⎦

⎤

⎢

⎣

⎡

+

⎥

⎦

⎤

⎢

⎣

⎡

=

k m k m k m

H

G

F

M

其中，矩陣 A內各矩陣及矩陣 F、 G、 H內各矩陣表示如下， ) ( ) ( 2 2 ) ( 2 ) ( 2 2 2 3 3 i j i i j i j mm ij u y i m y t y y m A α φ α π φ α π φ α α π ₊ ″ ₋ ∆ − = ) ( ) ( 2 2 ) ( ) ( 2 2 i j i i j i u y i m y t y m i y u t _α _φ α π φ α α π ∆ ″ ₋ ∆ ″ + ) ( Re 2 ) ( Re 2 2 4 4 2 2 i j i j m y t y m t ₊ ∆ ″ ∆ − α φ α π φ α α π ) ( Re 2 2 "" i j y t _φ α π ∆ − (3.4)

∑

+ ″ − = n i n k mn n i n k mn m i m a t y a t y F 2 2 2 ( ) ( ) 2 ( )φ ( ) α π φ α α π

∑

∆ − ∆ ″ + n i n k mn i n i n k mn i u y i ma t y t y t a m i y u t ) ( ) ( ) ( 2 2 ) ( ) ( ) ( 2 2 3 3 α φ α π φ α α π

∑

∆ ″ + ∆ + n i n k mn n i n k mn i m a t y t y t ma i y u t ) ( ) ( Re 2 2 ) ( ) ( ) ( 2 2 α4 4 φ α π φ α α π ) ( ) ( Re 2 2 ) ( ) ( Re 2 2 2 "" i n k mn n n i n k mn a t y t y t a m t _φ α π φ α α π ∆ + ″ ∆ −

∑

(3.5)

(22)

⎢ ⎣ ⎡ ∆ = ) (

∑∑

( ) ( ) ) 2 3 ( m n mx i i n k mn m i i ma t y e t G α φ α

∑∑

− ′ + ⋅ m n mx i i n k mn p q px i i q k pq t y e i ma t y e a p ( ) ( ) ) ( ( ) ( ) ) ( α2 2 φ α α φ α ) ) ( ) ( ( ) ) ( ' " ) ( (

∑∑

+

∑∑

′ ⋅ m n mx i i n k mn p q px i i q k pq t y e a t y e a φ α φ α

∑∑

−

∑∑

′ ⋅ p q m n mx i i n k mn px i i q k pq t y e a t y e a p iα3 3 ( )φ ( ) α ( ( )φ ( ) α ) ⎥ ⎦ ⎤ ″ ⋅(

∑∑

( ) ( ) ) p q px i i q k pq t y e pa iα φ α (3.6) ⎢ ⎣ ⎡ ′ ∆ = ) (

∑∑

( ₋ ) ( ) ) 2 ( ₁ m n mx i i n k mn m i a t y e t H φ α

∑∑

− +

∑∑

′ ⋅ ₋ ₋ p q m n mx i i n k mn px i i q k pq t y e a t y e a p i ( ) ( ) ) ( ( ) ( ) ) ( α3 3 ₁ φ α ₁ φ α ) ) ( ) ( ) ) ( ) ( (

∑∑

₋₁ ″ +

∑∑

₋₁ ⋅ m n mx i i n k mn p q px i i q k pq t y e i ma t y e pa iα φ α α φ α

∑∑

− − ′ − − ⋅ m n mx i i n k mn p q px i i q k pq t y e i ma t y e a p ( ) ( ) ) ( ( ) ( ) ) ( α2 2 ₁ φ α α ₁ φ α ⎥ ⎦ ⎤ ⋅(

∑∑

( ₋1) "'( ) ) p q px i i q k pq t y e a φ α (3.7) (3.4)、 (3.5)是表示為線性項部份， (3.6)、 (3.7)則是非線性項部份，其中 m=1,…,M， i,j=1,…,N，如此我們便可以先將線性項部份之解求出後再將非線性項加入，而未知數之最新一次的值，可由前兩次的結果求得。將求得之對時間系數代入流線函數 (3.1)中，根據流線函數的定義可計算出速度分量

x

v

y

u

∂

−

=

∂

=

ψ

,

ψ

(23)

3.2 數值解的分析方法

3.2.1 時間級數與 FFT頻譜分析

動力系統可藉由其與時間有關的函數 f(t)或是其動力變數的時間級數 (Time series)來表示其隨時間演變的過程。而任何時間函數 f (t)通常可以表示為許多週期性分量的疊代組成，而求取其中各週期性分量的相對強度則稱為頻譜分析 (Spectral analysis)，而頻譜分析為研究混沌的一個重要方法。在此藉由傅立葉轉換 (Fourier transform)來求取各分量頻譜強度，將欲求的時間函數轉換為一頻率函數，並分離其分量強度以描繪我們所需之頻譜強度圖。當一個時間函數 f(t)為一週期性函數時，該頻譜圖將會出現一個頻率，而當時間級數為一非週期性函數時，該頻譜圖則會以振盪的連續頻率，而在圖上為連續曲線。

吾人在此採用快速傅立葉轉換 (Fast Fourier transform)(簡稱 FFT)，藉以求取我們所需之頻譜強度圖，該演算法是由 70 年代的學者所提出，其為離散傅立業轉換 (Discrete Fourier transform)所演化而來，但其運算速率更為有效率，可較快獲取我們所要的轉換函數。我們可藉由轉換後之頻譜圖來判讀動力系統的運動形態，若為一個同頻 (Synchronous)的週期性振盪，則產生一個頻率；而當次調諧 (Subharmonic)的週期性運動時，除了原來的頻率外將產生第二個頻率，而對於進入混沌狀態的系統，其頻譜圖則是無法判讀的連續不規則圖形。

3.2.2 相平面圖

Baker 和 Gollub (1990)研究中提出，動力系統中的相空間 (Phase space)是一種數學上正交座標系統下的空間，用以確立動力系統中相同瞬時狀態下之各時間變數，而相平面 (Phase

(24)

plane)的觀念則是來自於相空間。例如一質點運動於一維座標下運動，其位置(x)與速度 (v)則是可表示其狀態之相空間，而此空間為一平面；而當質點運動於三維座標下時，則其需要六個維度的相空間來描述其瞬時狀態，其中包含三個位置分量及三個速度分量，而相空間的架構有許多種方式，例如動量可以取代速度分量。對於複雜參數的動力系統，其相平面通常極為紊亂，此時須藉由其他方法才較為容易判讀系統狀態。

3.2.3 Poincare 映射圖

Poincare 映射圖是一種用來簡化相平面圖的方法，藉由固定週期下，對相平面的軌跡取一投影來取樣觀察 (strobe)，而在此我們是以內普修流之非線性週期來做取樣，使得我們可以清楚的觀察動力系統的演變過程。由 Poincare 映射圖我們可以得知動力系統自然頻率(ω₀) 對於取樣頻率 (ωs)的比例關係。當動力系統的自然頻率為 s ) q p ( ω 時 (p q為有理數 )，則於 Poincare 映射圖上會顯現q個映射點於一個封閉軌道之上，而剩下[q−(p+1)]個點則會與先前的映射點位置重複，是故可以忽略這些映射點。對一個消散 (Dissipative) 系統而言，其映射點會移動至適當的吸子 (Attractor)而形成一點，而當取樣頻率與自然頻率成有理數關係時，則會形成散佈在封閉軌道上的點群，若完全不相關時，其映射圖上的點甚至會填滿而形成一封閉曲線。

(25)

第四章

結果與討論

本章中吾人將探討平面 Poiseuille 流在低於以及超過臨界雷諾數的物理現象，在線性擾動方程式(2.16)和非線性擾動方程式(2.14)部份我們改變雷諾數去分析對流場所產生的影響。

4.1 線性擾動方程式及準確度分析

首先吾人考慮平面普修流的線性擾動方程式，在前一章數值方法 中 x 方向與 y 方向以複數形式 Fouries series 及 Chebyshev polynomials

雙重級數展開的項數(M×N)將影響數值計算結果的準確度，因此我們必需改變(M×N)的項數去尋找準確的臨界雷諾數。理論上，展開項數取的愈多其結果之準確度愈高，但項數取太多將耗費大量的時間在運算上，因此取最少的展開項數而能達到相當準確度才是最佳的選擇。考慮展開項數為(3× 項，5) α ＝1.020545 時發現臨界雷諾數發生在 667，然後提高 N 項，一直到(3×31)項，這時可發現臨界雷諾數為 5772.22，且提高 N 項時臨界雷諾數會上下振盪如圖二的結果所示，當取到(3×31)或(3×33)項時臨界雷諾數都為 5772.22。吾人另外再取展開項數為(5×5) 項，α ＝ 1.02054 這個結果顯示當展開項數為 ) 31 31 ( × 、(33×33)項的臨界雷諾數也發生在 5772.22 如圖三，展開項數(3×31)及(31×31)項所運算出來的結果都和 Orszag(1971)吻合，為了運算時間考量以下取至(3×31)項即可。流場在低雷諾數時，速度隨時間會呈現快速收斂如圖四(a)，與 Fortin(1994)的結果一致。提高雷諾數，流場振盪變得更劇烈，收斂至穩定值的時間也隨之增加如圖四(b)，在物理意義上提高雷諾數代表擾動使流場的流動情形變得更劇烈，流場需要更長的時間達到穩定

(26)

狀態，因此時間函數收斂的時間增加與實際的物理現象符合。當雷諾數達到 5772.22 時流場流動的行為並不會隨著時間收斂或發散而以週期性運動發生如圖四(c)，當雷諾數高於臨界雷諾數時則時間對於速度函數會產生發散行為如圖四(d)所示。

4.2 非線性擾動方程式及高於臨界雷諾數之流場行為

為了解超過臨界雷諾數後流場的行為，吾人將非線性項加入，而展開項數取(3×31)項，所得到的臨界雷諾數為 5835。再提高 M 項試著尋找臨界雷諾數，在(7×31)項時得到臨界雷諾數為 5772.22，如圖五。而流場在低雷諾數時，速度隨時間會呈現快速收斂如圖六(a)。當提高雷諾數，流場振盪變得更劇烈，收斂至穩定值的時間也隨之增加如圖六(b)。以下結果為用在(7×31)項探討超過臨界雷諾數流場的行為，包含單一週期性(periodic)、擬週期性(quasi-periodic)、混沌 (chaos) 現象，吾人將結果以時間級數圖配合 FFT 來獲取頻譜圖，而後將時間級數解投影作成相平面圖與 Poincare 映射圖，藉以瞭解流場之行為。

4.3 由時間級數與 FFT 頻譜圖探討流場行為

圖七為加入非線性項後，雷諾數 Re=5772.22，時間 dt=0.2，波數 α ＝1.020545 的時間級數圖與頻譜圖，時間級數圖顯示擾動並不會隨著時間收斂，呈現單一週期性的波動變化，吾人從 FFT 頻譜圖上看到頻率為單一個頻率 f1=64.1，確為單一週期性運動型態。圖八為雷諾數 Re=6000，時間 dt=0.2，波數α ＝1.020545 的時間級數圖與頻譜圖，在時間級數圖上的波動無法判斷出是什麼形態的運動，然而從 FFT 頻譜圖可看出有二個不相干頻率出現 f₁=64.1，

(27)

2 f =96.4，因此我們可判斷此時流場為擬週期性運動(quasi-periodic)。圖九、圖十分別為雷諾數 Re=7000、Re=8000，時間 dt=0.2，波數α ＝1.020545 的時間級數圖與頻譜圖，吾人在時間級數圖上看到波動振盪大小增加，藉由 FFT 頻譜圖看到有二個不相干頻率出現 1 f =64.1，f₂=96.4，所以我們可判定為擬週期性運動。圖十一、圖十二分別為雷諾數 Re=9000、Re=10000，時間 dt=0.1，波數α ＝1.020545 的時間級數圖與頻譜圖，由於在時間級數圖上看到波動更加劇烈，藉由 FFT 頻譜圖上看到二個不相干頻率出現 1 f =64.1，f₂=96.4，只是第二個頻率的強度增加，因此我們可判定為擬週期性運動。圖十三為雷諾數 Re=11000，時間 dt=0.1，波數α ＝1.020545 的時間級數圖與頻譜圖，從在時間級數圖上看到的波動開始呈現微弱的混亂，吾人無法判讀是為擬週期性運動或者是混沌現象，藉由 FFT 頻譜圖亦無法找到其規則性，這時流場開始產生微弱的混沌現象。圖十四為雷諾數 Re=12900，時間 dt=0.1，波數α ＝1.020545 的時間級數圖與頻譜圖，從在時間級數圖上看到的波動混亂無法判讀，藉由 FFT 頻譜圖亦無法找到其規則性，這時流場的混沌現象逐漸增強。

4.4 由相平面圖與 Poincare 映射圖探討流場的行為

圖十五為雷諾數 Re=5772.22，時間 dt=0.2，波數α ＝1.020545 的相平面圖與 Poincare 映射圖，我們將圖七與圖十五互相比較，可發現在圖七(b)中頻譜圖顯示流場當時是處於單一週期性運動狀態，在圖十五(b)中顯示於 Poincare 映射圖上只有一點，而相平面圖也呈現一封閉曲線如圖十五(a)，所以我們可判斷為單一週期性運動狀態。圖十六為雷諾數 Re=6000，時間 dt=0.2，波數α ＝1.020545 的相

(28)

平面圖與 Poincare 映射圖，我們將圖八與圖十六互相比較，可發現在圖八(b)中頻譜圖顯示流場是處於擬週期性運動狀態，而在圖十六(b) 中顯示於 Poincare 映射圖上有二點，而相平面圖呈現規則性二條的封閉曲線如圖十六(a)。因此我們可判斷其為擬週期性運動狀態。圖十七、圖十八為雷諾數 Re=7000、Re=8000，時間 dt=0.2，波數α ＝1.020545 的相平面圖與 Poincare 映射圖，我們將圖九與圖十七、圖十與圖十八互相比較，可發現在圖九(b)、圖十(b)中頻譜圖顯示流場是處於擬週期性運動狀態，而在圖十七(b)、圖十八(b)中顯示於 Poincare 映射圖上有二點，而相平面圖呈現規則性二條的封閉曲線如圖十七(a)、圖十八(a)。因此我們還是可判斷其為擬週期性運動狀態。圖十九、圖二十為雷諾數 Re=9000、Re=10000，時間 dt=0.1，波數α ＝1.020545 的相平面圖與 Poincare 映射圖，我們將圖十一與圖十九、圖十二與圖二十互相比較，可發現在圖十一、圖十二中，時間級數圖看到的波動更加劇烈，頻譜圖顯示流場是處於擬週期性運動狀態，而在圖十九(b)、圖二十(b)中顯示於 Poincare 映射圖上有二點，而相平面圖呈現規則性二條的封閉曲線如圖十九(a)、圖二十(a)。因此我們還是可判斷其為擬週期性運動狀態。圖二十一為雷諾數 Re=11000，時間 dt=0.1，波數α ＝1.020545 的相平面圖與 Poincare 映射圖，我們將圖十三與圖二十一互相比較，可發現圖十三中時間級數圖上看到的波動開始呈現微弱的混亂，頻譜圖上的頻率出現些為的雜亂，而圖二十一(b)中顯示 Poincare 映射圖則形成無數點的狀態，而相平面圖則由雜亂的封閉曲線所組成如圖二十一(a)。此時流場已經進入微弱的混沌狀態。圖二十二為雷諾數 Re=12900，時間 dt=0.1，波數α ＝1.020545 的相平面圖與 Poincare 映射圖，我們將圖十四與圖二十二互相比較，可發現圖十四(b)中頻譜圖上的頻率雜亂無法判斷，而圖二十二(b)

(29)

中顯示 Poincare 映射圖則形成無數點的狀態，而相平面圖則由更雜亂的封閉曲線所組成如圖二十二(a)，此時已無週期性可言流場已經進入混沌狀態。

(30)

第五章

結論

本文以 Chebyshev polynomials 及複數形式的 Fourier 級數對空間座標作雙重級數展開來模擬平面普修流在低於以及超過臨界雷諾數的物理現象，藉由改變雷諾數的大小來觀察流場形態，根據第四章的結果可歸納以下幾點結論： 1．流場在低雷諾數時，速度隨時間會呈現快速收斂，當提高雷諾數時，時間函數收斂至固定值的時間增加，且時間函數振盪的幅度較低雷諾數時大，與 Fortin(1994)的結果一致。 2．在臨界雷諾數(Re =5772.22)時，所觀察到的流場行為顯示，時c 間函數並不會隨著時間而收斂，呈現單一週期性的狀態，與 Orszag 及 Fortin 的結果一致。當雷諾數至 6000~10000 範圍內，流場流動行為愈來愈劇烈，並開始產生第二個不相干的頻率呈現擬週期性運動。 3．當雷諾數為 11000 時，流場擾動的情形更加劇烈，且呈現出微弱的混沌現象，當臨界雷諾數為 12900 時，我們看到流場擾動混亂而無法判讀，這時流場的混沌現象更加劇烈。

(31)

參考文獻

Baker, G. L. and Gollub, J. P. 1900 Chaotic dynamics an introduction Cambridge university press.

Barenghi, C. F. and Jones, C. A. 1989 Modulated Taylor-Couette flow. J. Fluid Mech. 208, 127-160.

Betchov, R. and Criminale, W. O. Jr., 1967 Stability of parallel flows.Academic Press.

Chandrasekhar, S. 1961 Hydrodynamic and Hydromagnetic Stability. 637-642.

Chock, D. P. & Schechter, R. S. 1972 Critical Reynolds number of the Orr-Sommerfeld equation. The Physics of Fluids. 16-2, 329.

Davey, A. and Drazin, P. G. 1969 The stability of Poiseuille flow in a pipe. J. Fluid Mech. 36, 209-218.

Davis, S. H. 1976 The stability of time-periodic flows. Ann. Rev. Fluid Mech. 8, 57-74.

Fortin, A. and Jardak, M. 1994 Old and new results on the two-Dimensional Poiseuille flow.J.Computational Physics. 115, 455. Grosch, C. E. and Salwen, H. 1968 The stability of steady and

time-dependent plane Poiseuille flow. J. Fluid Mech. 34, 177.

George, W. D. and Hellums, J. D. 1972 Hydrodynamic stability In Plane Poiseuille flow with finite amplitude disturbances. J. Fluid Mech. 51-4, 687.

Georgescu, A. 1985 Hydrodynamic stability theory. Martinus Nijhoff Publishers.

(32)

Jimenez, J. 1990 Transtion to turbulence in two-dimensional Poiseuille flow. J. Fluid Mech. 218, 265.

Kuhlmann, H. 1984 Model for Taylor-Couette flow. Phy. Rev. A 32, 1703-1707.

Lessen, M. Fox, J. A. Bhat, W. V. & Liu, T. Y. 1964 Stability of Hagen-Poiseuille Flow. Phys Fluids. 7, 1384.

Lin, C. C. 1954 The Theory of Hydrodynamic Stability. Cambridge Monographs On Mechanics And Applied Mathematics.

Nobutake Itoh, 1977 Nonlinear Stability of parallel flows with subcritical Reynolds numbers. Part 1. An asymptotic theory valid for small amplitude disturbance. J. Fluid Mech. 82-3, 455.

Orszag, S. A. 1971 Accurate solution of the Orr-Sommerfeld stability equation. J. Fluid Mech. 50, 689.

Reid, W. H. and Bart, S. 2002 On the spectral problem for Poiseuille flow in a circular pipe. Fluid Dynamics Research 33, Issue:1-2.

Shen, S. F. 1955 Calaulated amplified oscillations in plane Poiseuille and Blasius flow. J. Aero. Sci. 21, 62.

Soibelman, I. and Merion, D. 1991 Finite-amplitude bifurcations in plane Poiseuille flow: two-dimensional Hopf bifurcation. J. Fluid Mech. 229, 389.

Stewartson, K. 1971 A non-Linear instability theory for a wave system in plane Poiseuille flow,J. Fluid Mech. 48, 529.

Thomas, L. H. 1953 The stability of plane Poiseuille flow. Phy. Rev. 91, 780.

Wolf, L. Jr. and Lavan, Z. and Nielsen, H. J. 1978 Numerical Computation of the Stability of Plane Poiseuille Flow. J..Applied

(33)

李鎰清 2001 圓柱之流場-一種改良 Lorenz 模型. 國立交通大學碩士論文.

林紀元 2004 水平同心圓柱間自然對流之分歧現象.國立交通大學碩士論文.

(34)

圖一普修流模型

dP

c

dX

−

=

X

Y

2H

(35)

N 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 Rec 0 1000 2000 3000 4000 5000 6000 7000 8000 圖二 當 M=3，dt=0.5，α＝1.020545 時不同 N 所得到的臨界雷諾數

(36)

M&N 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 Rec 0 1000 2000 3000 4000 5000 6000 7000 圖三 Rec對M &N ，dt=0.5，α＝1.020545 時所得到的臨界雷諾數

(37)

(a)

(b)

圖四當 dt=0.5，α ＝1.020545 時，速度隨時間變化圖，Re=

(38)

(c)

(39)

M 1 3 5 7 9 11 Rec 5720 5740 5760 5780 5800 5820 5840 圖五當 N=31，dt=0.5，α＝1.020545 時不同 M 所得到的臨界雷諾 數

(40)

(a)

(b)

圖六加入非線性項後 dt=0.5，α ＝1.020545 時，速度隨時間變

(41)

(a)

(b)

圖七加入非線性項後,雷諾數 Re=5772.22，時間 dt=0.2，波數α ＝

(42)

(a)

(b)

圖八加入非線性項後,雷諾數 Re=6000，時間 dt=0.2，波數α ＝

(43)

(a)

(b)

圖九加入非線性項後,雷諾數 Re=7000，時間 dt=0.2，波數α ＝

(44)

(a)

(b)

圖十加入非線性項後,雷諾數 Re=8000，時間 dt=0.2，波數α ＝

(45)

(a)

(b)

圖十一加入非線性項後,雷諾數 Re=9000，時間 dt=0.1，波數α ＝

(46)

(a)

(b)

圖十二加入非線性項後,雷諾數 Re=10000，時間 dt=0.1，波數α ＝

(47)

(a)

(b)

圖十三加入非線性項後,雷諾數 Re=11000，時間 dt=0.1，波數α ＝

(48)

(a)

(b)

圖十四加入非線性項後,雷諾數 Re=12900，時間 dt=0.1，波數α ＝

(49)

(a)

(b)

圖十五雷諾數 Re=5772.22，時間 dt=0.2，波數α＝1.020545

(50)

(a)

(b)

圖十六雷諾數 Re=6000，時間 dt=0.2，波數α＝1.020545

(51)

(a)

(b)

圖十七雷諾數 Re=7000，時間 dt=0.2，波數α＝1.020545

(52)

(a)

(b)

圖十八雷諾數 Re=8000，時間 dt=0.2，波數α＝1.020545

(53)

(a)

(b)

圖十九雷諾數 Re=9000，時間 dt=0.1，波數α＝1.020545

(54)

(a)

(b)

圖二十雷諾數 Re=10000，時間 dt=0.1，波數α ＝1.020545

(55)

(a)

(b)

圖二十一雷諾數 Re=11000，時間 dt=0.1，波數α ＝1.020545

(56)

(a)

(b)

圖二十二雷諾數 Re=12900，時間 dt=0.1，波數α ＝1.020545

平面普修流之非線性不穩定

國立交通大學

機械工程學系

碩 士 論 文

平面普修流之非線性不穩定

Nonlinear Instability of Plane Poiseuille Flow

研 究 生：洪英棋

指導教授：楊文美 博士

平面普修流之非線性不穩定

Nonlinear Instability of Plane Poiseuille Flow

摘 要

誌謝

目錄

圖 目 錄

符 號 說 明

符 號

上標說明

下標說明

第 一 章

緒 論

1.1 文獻回顧

1.2 研究目的

第 二 章

數學模式

2.1 統御方程式

ρ

µ

ρ

µ

2.2 無因次化

ρ

2.3 基態解

2.4 擾動方程式之建立

ψ

ψ

ψ

ψ

ψ

ψ

ψ

ψ

ψ

ψ

1

Re

P

x

x y

y

ψ

ψ

′

⎛

⎞

∂

∂

∂

= −

+

⎜

+

⎟

∂

⎝

∂ ∂

∂

⎠

ψ

ψ

ψ

ψ

ψ

ψ

1

Re

P

y

x

y x

ψ

碩士論文

研究生：洪英棋

指導教授：楊文美博士

摘要

圖目錄

符號說明

符號

第一章

緒論

第二章

_⎜

_⎟

_⎝

_⎠

_⎜

_⎟

_⎝

_⎠

_ψ

_ψ