平面普修流之動態模型

國立交通大學

機械工程研究所

碩 士 論 文

平面普修流之動態模型

Dynamic Models of Plane Poiseuille Flow

研 究 生：古 勤 暐

指導教授：楊 文 美 博士

平面普修流之動態模型

Dynamic Models of Plane Poiseuille Flow

研 究 生：古 勤 暐 Student：Ching-Wei Ku 指導教授：楊 文 美 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 2004

Hsinchu, Taiwan, Republic in China

摘 要

普修流是流體力學中最基本也是最重要的問題，當學者們發展出 新的實驗方法或是數值方法，普修流這樣基本的問題成為了驗證新方 法正確性的依據。本文以數值方法探討平面普修流流場在不同雷諾數 的行為，文中以三組 Lorenz 模型簡化統御方程式為聯立非線性常微 分方程組，採用 Runge-Kutta 數值方法求解，再利用時間級數圖觀察 流場的行為，並藉由求取線性化方程組的特徵值，及以牛頓法求取穩 態解驗證吾人所建立的模型。 結果發現，三組模型在小於臨界雷諾數時符合物理現象，在超越 臨界雷諾數後沒有非零解產生。吾人求取線性化方程組的特徵值，亦 發現超越臨界雷諾數沒有不穩定的特徵值出現。利用牛頓法求取模型 的穩態解，得到皆為零解，因此驗證此模型未能夠得到臨界雷諾數之 後的合理現象。推測原因，可能是因為展開項數不足，或使用之展開 函數不適當所造成。

誌謝

這篇論文的誕生，首先要感謝我的指導教授－楊文美博士在這一 段期間給予勤暐的指導。小至論文中的排版、符號安排，大如論文中 的程式，老師都能夠找到盲點並提供指正，使得吾人獲易良多，並且 在習得流體力學深入的學問之外，更了解了做學問的方法。 承蒙口試委員盧定昶博士以及游明輝博士於口試期間給予本論 文的寶貴意見及指教，使得本文內容得以更加完善及正確，在此由衷 的表示感謝。 以及感謝鎰清學長，即使是在大夜班後的早晨仍然親自來到實驗 室和我們討論程式，你的鼓舞與指導一直是我們完成論文的一大動 力。完成論文的現在，雖然已經來不及和你分享我的喜悅和感動，勤 暐還是會一直記得和感謝你。 感謝盈嘉在撰寫論文時給我的支持，在煩悶時妳總能傾聽我的心 情，讓我能夠在遇到挫折時再站起來繼續努力。感謝余伯、志凱、紀 元、朝任學長，以及同窗盈立、勝文、宏文、柏霖、永賢、昱志、家 宏、家棟、碩一同學宏仁、英棋、友約、自皓，大學同學朝雯、秋華、 國昌，在功課以及生活上陪我渡過這段難忘的日子，你們為我的生命 添上更多的色彩。 最後僅將本論文獻給我親愛的父母及家人，有你們的照顧與關心 我才能夠順利完成本文，感謝你們。

目 錄

摘要……….i 誌謝………ii 目錄………...………iii 表目錄………...……….……v 圖目錄………..……….vi 符號說明………...…..………….…….vii 第一章 緒論………..………1 1.1 文獻回顧……….……….1 1.2 研究目的……….……….5 第二章 數學模式………..6 2.1 統御方程式………..……6 2.2 無因次化………..………7 2 . 3 基 態 解 … … … . . 7 2.4 擾動方程式之建立……….8 2.5 Lorenz 模型的建立………10 第三章 數值方法……….………...16 3.1 Runge-Kutta Method………...………..………16 3.2 Newton’s Nethod………...……….18 3.3 時間級數圖………...……….19 3.4 線性化與求取特徵值………....19

第四章 結果與討論………20 4.1 以時間級數圖觀察模型行為……….………20 4.2 以線性化模型判斷解的正確性……….………21 4.3 牛頓法的驗證……….………22 4.4 超過臨界值無法獲得非零解的可能原因………23 第五章 結論………25 參考文獻………26 附錄………39

表 目 錄

表一 線性化之 Model(1)，各時間函數在不同雷諾數之特徵值……29

表二 線性化之 Model(2)，各時間函數在不同雷諾數之特徵值……30

圖 目 錄

圖一 基本物理模式示意圖………...……….32

圖二 Model(1)，函數 E 在不同雷諾數之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

………...……….………...……….33

圖三 Model(1)，函數 F 在不同雷諾數之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

………...……….………...……….34

圖四 Model(2)，函數 F 在不同雷諾數之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

………...……….………...……….35

圖五 Model(2)，函數 G 在不同雷諾數之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

………...……….………...……….36

圖六 Model(3)，函數 E 在不同雷諾數之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

………...……….………...……….37

圖七 Model(3)，函數 F 在不同雷諾數之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

符號說明

符 號 , ... A B J _時間函數 c 常數 i f Runge-Kutta 法中 Lorenz 變數之函數 g 重力(y 方向) H 平行板之間距 ij k 運算變數 m,n 流線函數展開項之項數 P 壓力 Re 雷諾數 t 時間 , u v _速度分量 x,y 平面座標 φ _{y 方向展開函數} ψ _{流線函數(Stream Function)} µ _黏滯係數 ρ 密度 α 波常數(wave number) ∆ 2 2 2 2 x y ∂ ₊ ∂ ∂ ∂ 2 ∆ 4 4 4 2 2 2 4 4 x x y y ∂ ₊ ∂ ₊ ∂ ∂ ∂ ∂ ∂ 上 標 基態解 ' 擾動量

第 一 章

緒 論

1. 1 文獻回顧

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

Sexl(1927a)將在圓管普修流穩定性的問題理論化，他針對軸對稱 之圓管流推導出相關之方程式，而 Corcos 和 Sellars(1959)以及 Gill (1965)的論文深入地探討了 Sexl 的問題。之後 Sexl(1927b)和 Lessen, Fox, Bhat 以及 Liu(1964)探討了較為複雜的非對稱管流的問題。在實 驗以及理論對照下，他們發現了普修流在受到極微小的擾動下呈穩 定，但在一般的擾動下則呈現不穩定的狀態。

時至 1950-1960 年代，計算機科技初步發展，學者在試著以各種 數值方法探討線性普修流的問題，針對 Orr-Sommerfeld Equation 找出 其發生不穩定的臨界雷諾數。首先 Shen(1954)使用由 Lin(1955)所發 表 的 近 似 分 析 (asymsptotic analysis)法 得 到Rc=5360，αc＝1.05，

Lin並且建立了中性曲線(neutral curve)；Thomas 率先以有限元素法 得到Rc =5780，αc＝1.026；Betchov & Criminale(1967)也以有限元素

法得到R_c=5767，αc＝1.02。Grosch & Salwen(1968)利用 Chandrasekhar

-Reid(1961)正交函數得到Rc=5750，αc＝1.025。他們並探討了在流場 加入水平方向調制壓力梯度對臨界雷諾數的影響。結果發現在調制振 幅大過於 0.105 時( 為該文中定義之最大速度)，會造成劇烈的減 穩定(destabilization)效果，並推斷為外加振盪作用與流體共振的結果。 Kerczek(1981)利用 Chebyshev 展開式探討同樣的問題，發現在臨界點 附近，在調制頻率為中高頻時造成流場穩定化(stabilization)的作用， 在較低的調制頻率範圍則造成輕微的簡穩定，與 Grosch & Salwen 相 異，作者推測是由於計算方式的不同所造成。

0

U U₀

Orszag(1971)以多項的 Chebyshev 展開式搭配 QR 矩陣特徵值演 算法得到該式臨界雷諾數發生於Rc=5772.22，αc＝1.02056。他以 Orr-

Sommerfeld Equation 式中 Im(λ)之正負值判斷解的穩定與否。當 Re

=5772.22 時，Im(λ)為負值，雷諾數超過此值所得到 Im(λ)全為正。

Orszag 的結果具重大的參考價值，證明 Chebyshev 展開式的精確性、 及快速收斂的特性，並可節省大量 CPU time，適合用來求解流體運 動的問題。

隔年 Chock & Schechter(1972)假設 Orr-Sommerfeld Equation 式中 (neutral stability)的情形下，將該式置換為四條線性聯立方程式， 並使之符合邊界條件，利用四階的 Runge-Kutta 法求得 0 i c = c R=5772.225±

0.001，αc＝1.0205±0.001，與 Orszag 的結果相近。Chock & Schechter

的計算方式雖然較為精簡，但文中提及該方式的缺陷在雷諾數到達 15000 以上時，其解已無參考價值，作者建議採用 Chebyshev 展開式 或有限元素法或求解較為恰當。

Davis(1976)廣泛地探討周期性流場(time - periodic flow)的穩定 包含了史托克層(Stokes layer)、帶有調制壓力梯度的普修流、以及振 盪邊界層等此類形式之問題。Davis 發現改變液動系統的振福和頻率 造成增穩定或減穩定的現象，能應用在工程上，減穩定能夠增加系統 的熱、質、能量的傳遞；增穩定能夠提高系統的效率。

對於非線性 Navier Stokes 方程式的探討(2.16 式)最先為 Meksyn & Stuart(1951)，而 Dowell(1969)接續了此研究並與更系統化地探討了 非線性項的加入對方程式的影響。Dowell 對流場的水平方向採 Fourier 展開，而垂直方向採用符合邊界條件的兩組餘弦函數(cosine)展開， 並使用了多項展開模式(spectral method)試著精確地描述流場現象。在

不含非線性項的情形，Re=80,000，α_c＝2 的條件下，Dowell 使用了

四十項 y 方向奇函數展開得到與 Thomas 使用 100 個有限差分點(finite -difference points)相同的特徵值值，說明 spectral method 的可行性。 在包含非線性項的同樣條件，y 方向取十六項展開的情形，發現第一 階展開項係數對時間解呈 limit cycle，而第二階係數則收斂至一固定 值。Dowell 另外發現垂直方向若有奇(odd)函數的展開能引發非線性 的擾動現象，若僅使用偶函數則無法引發此現象。Dowell 以當時計 算機速度無法得到非常理想的結果，並建議在先由線性方程式計算出 合適的 eigenfunction 再代入完整的方程式中，能夠減少使用的展開項 數。 Fortin 與 Jardak(1994)對於普修流的穩定性分析做了一個完整的

統合。他們將問題分為一維及二維、線性及非線性，首先以數值方法 解 出 線 性 的 Orr-Sommerfeld equation， 得 到 臨 界 雷 諾 數 發 生 於 1.020545 α= ， ，與 Orszag 的結果吻合，而後採用疏密不同 的網格，利用有限元素法直接解出非線性的擾動方程式，並且將完整 的流場繪出。結果顯示在低雷諾數(Re=525)時流場呈現穩定的流動， 擾動在一定的時間內達到收斂：在雷諾數提高到臨界值時，時間級數 圖顯示擾動並不隨時間收斂，利用 FFT(Fast Fourier Transform)及相平 面圖(Phase plane)得知流動行為僅有一組基頻與其倍數頻率的組合， 為單一週期性的運動，稱為週期性(periodic)流場。當雷諾數提高至 6000 時，出現第二組頻率，流場轉變為擬週期性(quasi-periodic)的運 動。當雷諾數為 8000 至 10000 時流場依然為擬週期性的運動，不同 之處僅在於流動行為更加劇烈。雷諾數提高至 11000 時，流場呈現微 弱混沌現象(weakly chaos)，即流場已產生兩組以上之頻率。當雷諾數 達到 12900 時，流場已產生混沌現象，無法由 FFT 及相平面圖解讀 其現象。Fortin 的結果與前人吻合，並在超過臨界雷諾數後有完整的 分析，為近年研究普修流非線性現象的重要指標。 R_c =5772.22 有關動力系統的研究，氣象學家 Lorenz(1963)首先以簡單的三角 函數對動量以及能量方程式做傅立葉展開，得到一組三式的非線性聯 立方程式，模擬了流體在自由邊界(free-free boundary)的運動情形。 方程式中 Rayleigh Number 超過臨界值時會得到難以解讀的混沌現 象，引起人們對於簡單行為背後的複雜現象的興趣，開啟後世研究動 力系統及混沌學的先河。 周煒超(2001)與李鎰清(2001)分別以四項以及九項展開式模擬兩 同心圓柱間流場在調制作用下的動力模型。他們發現在外管固定，內 管做調制轉動情形下，當運動系統的雷諾數大於非調制臨界雷諾數

時，流場會先做擬週期性振盪運動，隨著雷諾數的增大，運動系統會 由擬週期性振盪轉變成與調制無關的週期性振盪，最後進入混沌狀 態。而流場中調制頻率的增加，使得流場運動形態轉變較快，或是調 制振幅較大時，流場內吸子的擴張及收縮效應較大，皆易使得流場進 入混沌的狀態之中。而周的結果與 Kuhlmann et al.及線性理論有些許 誤差，李的結果則更貼近 Kuhlmann。由 Branghi 和 Jones(1989)文獻 中發現這是由於在動力模型簡化運動方程過程時，展開項數的多寡會 影響結果的準確性。

1.2 研究目的

本文之主要研究目的為試圖利用 Lorenz 方程組探討二維平板普 修流的穩定性。普修流的穩定性問題早在 1970 年代左右已深入地被 探討，學者們從實驗、理論、數值方法切入此問題在不同雷諾數下， 是否符合實際物理現象。在參閱相關文獻後發現數學及物理學家大多 直接以較為複雜的數值方法解出擾動方程式並加以探討，而吾人在本 文欲嘗試以較為精簡之 Lorenz 方程組模擬出平面流場在低於以及超 過臨界雷諾數的物理現象，期望能夠得到與 Fortin 一文中相近的結 果。 第二章首先建立運動方程組，並將之無因次化，再加入擾動項以 建立擾動方程組。吾人嘗試在平板之 y 方向以正交及非正交兩種不同 函數，x 方向以傅立葉作展開推導成六項及七項 Lorenz 方程組，利用 第三章所介紹的 Runge-Kutta method 求解時間函數，在第四章討論所 得到的結果，並增加展開項數至十項，探討展開項數多寡是否造成不 同的解。

第 二 章

數 學 模 式

本文中研究的普修流物理模式如圖一所示。兩平行板距為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) 上式中 為壓力項，g 為重力加速度(gravitational acceleration)，P µ 為黏滯係數(viscosity)，ρ為流體密度(density)。 邊界條件(boundary conditions) ： ( , , ) ( , , ) 0 u x −H t =u x H t = (2.3)

2.2 無因次化

吾人在此定義無因次參數如下： x* x H = * y y H = * u u u∞ = * v v v_∞ = * 2 P P u

ρ

∞ = * u t t H ∞ = Re u H v ∞ = 將統御方程式無因次化，並除去上標後，為以下之形式： 連續方程式(continuity equation)： u v 0 x y ∂ ₊∂ ₌ ∂ ∂ (2.4) 動量方程式(momentun equations)： 2 2 2 1 ( Re u u u P u u v t x y x x y ∂ ₊ ∂ ₊ ∂ _{= −}∂ ₊ ∂ ₊ ∂ ∂ ∂ ∂ ∂ ∂ ∂ 2) u (2.5a) 2 2 2 1 ( Re v v v P v u v t x y y x y ∂ ₊ ∂ ₊ ∂ _{= −}∂ ₊ ∂ ₊∂ ∂ ∂ ∂ ∂ ∂ ∂ 2) v (2.5b) 邊界條件(boundary conditions) ： ( , 1, ) ( ,1, ) 0 u x − t =u x t = (2.6)

2.3 基態解

考慮流體在基態(basic state)為一維穩態運動，u=u y( )、 ，流 道中央處速度 0 v= u_∞為最大速度。基態速度及壓力梯度可以下式表示： 2 ( ) 1 u y = −y (2.7)

2 ( ) Re P x = − x _(2.8) 如此即解得u x y t( , , )之基態解。

2.4 擾動方程式之建立

吾人在本文中欲探討雷諾數大小對於系統造成不穩定性的問 題，分別於 、 、 加上擾動量，並代入統御方程式中，以得到本 節中之擾動方程式，並以此擾動方程式作為後續之統御方程式。 u v P 分別在基態上加上擾動量u x y t′( , , )、v x y t′( , , )、P x y t′( , , )： ( ) ( , , ) 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) 將三個速度分量分別代入(2.4)、(2.5)，並將此式減去(2.4)、(2.5)可得： 連續方程式(continuity equation)： u v 0 x y ′ ′ ∂ ₊∂ ₌ ∂ ∂ (2.10) 動量方程式(momentun equations)： 2 2 2 2 1 ( ) Re u u u u u u v v t x x y P u u x x y ′ ′ ′ ∂ ∂ _′∂ _′∂ _′∂ + + + + ∂ ∂ ∂ ∂ ∂ ′ ′ ′ ∂ ∂ ∂ = − + + ∂ ∂ ∂ u y ′ (2.11a) 2 2 2 2 1 ( ) Re v v v u u v t x x v y P v v y x ′ ′ ′ ∂ ∂ _′∂ _′∂ + + + ∂ ∂ ∂ ∂ ′ ′ ∂ ∂ ∂ = − + + ∂ ∂ ∂y ′ ′ (2.11b)

方程式(2.11)為描述擾動量之動量方程式，其邊界條件如下： 1, 1, 0 y= − y = u′= =v′ P′= _(2.12) 引用流線函數(stream function)

ψ

( , , )x y t ，和速度分量之關係 u y

ψ

∂ = ∂ 、v x

ψ

∂ = − ∂ ，將流線函數代入(2.11)中，可得： 2 2 2 2 3 3 2 3 1 Re u u t y x y y x y x y x y P x x y y 2

ψ

ψ

ψ

ψ

ψ

ψ ψ

ψ

ψ

∂ ₊ ∂ ₊∂ ∂ ₋∂ ∂ ₋∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ′ ⎛ ⎞ ∂ ∂ ∂ = − + _⎜ + _⎟ ∂ _⎝∂ ∂ ∂ _⎠ (2.13a) 2 2 2 2 2 3 3 3 2 1 Re u t x x y x x y x P y x y x 2

ψ

ψ

ψ ψ

ψ

ψ

ψ

ψ

∂ ∂ ∂ ∂ ∂ ∂ − − − + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ′ ⎛ ⎞ ∂ ∂ ∂ = − − _⎜ + _⎟ ∂ _⎝ ∂ ∂ ∂ _⎠ (2.13b) 將(2.13a)對 偏微分減去(2.13b)對 偏微分以消去y x P′，可得以下的方 程式： 2 1 Re u u t x x y x x y

ψ

ψ

ψ

ψ

ψ

ψ

ψ

ψ

″ ∂ _{∆ +} ∂ _{∆ −} ∂ ₊ ∂ ∂ _{∆ −}∂ ∂ _∆ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = ∆ (2.14) 其中 2 2 2 2 x y ∂ ∂ ∆ = + ∂ ∂ ， 4 4 2 4 2 2 2 4 4 x x y y ∂ ∂ ∂ ∆ = + + ∂ ∂ ∂ ∂ 其邊界條件為: 0

ψ

= at y = −1 1and y = (2.15) (2.14)式消去非線性項可得到 Orr-Sommerfeld equation: 2 2 2 1 Re u u t x y x

ψ

ψ

ψ

∂ ∂ ∂ ∂ ∆ + ∆ − = ∆ ∂ ∂ ∂ ∂

ψ

(2.16)

2.5 Lorenz 模型之建立

2.5.1 Legendre 正交函數展開 為符合邊界條件以及物理現象，將流線函數之 y 方向選擇以 Legendre 正交函數展開: 2 1 ( ) [ ( 1) ] 2 ! l l l l l d P y y l dy = − _(2.17) 其中l=0,1, 2...，代表 Legendre 函數之階數。 令一數列w y( )為 0 0 1 1 2 2 3 3 0 ( ) _{l l} ... l w y B P B P B P B P B P ∞ = =

∑

= + + + _(2.18) 代入邊界條件 ( 1) (1) ( 1) (1) w − =w = Dw − = Dw (2.19) 可得以下式： 0 1 2 3 4 ( 1) ... 0 w − =B − +B B −B +B = (2.20a) 0 1 2 3 4 (1) ... 0 w = B +B +B +B +B = (2.20b) 0 1 2 3 4 ( 1) 3 6 10 15 ... 0 Dw − = B − B + B − B + B = (2.20c) 0 1 2 3 4 (1) 3 6 10 15 ... 0 Dw = B + B + B + B + B = (2.20d) 化簡B B B B B0, 1, 2, 3, 4可得一新函數φl 4 ( ) _{l l} l w y B

φ

∞ = =

∑

_(2.21) 其中l = 4, 5, 6, 7...∞，在此為簡化方程式僅取φ_l前兩項 2 4 35 ( 1 8 y

φ

= − 2 ) _(2.22a)

2 5 63 ( 1 8 y y

φ

= − 2 ) _(2.22b) 而 x 方向則選擇以 in x e α 展開，其中l 、n為展開項之項數： 1 ( , , ) ( ) ( ) n in x ln l n l x y t A t y e α

ψ

=∞ ∞

φ

=−∞ = =

∑ ∑

_(2.23) 若以六項簡化，則 ( , , ) 4( ) 5( ) 1 4 2 5 i x i x x y t A y B y A eα A e α

ψ

_≈

φ

₊

φ

₊

φ

₊

φ

− 4 5 4 4 5 5 ( ) ( ) ( ) ( ) ( ) ( ) sin

( ) ( ) cos ( ) ( ) sin ( ) ( ) cos

A t y B t y C t y x D t y x E t y x F t y x

φ

φ

φ

α

φ

α

φ

α

φ

= + + + + +

α

(2.24) 其中

cos x

α

及sin x

α

代表水平方向 x 之週期性，其波長為2 /

π α

，

α

為 波常數(wave number)。 在此使用 Galerkin 法將方程式中之空間函數消除，將（2.24）式 代入（2.14）式，再乘上

φ

4、

φ

5、

φ

4sin x

α

、

φ

4cos x

α

、

φ

5sin x

α

、

φ

5cos x

α

而後對六式對

x

及y 積分， 方向積分範圍為x −

π α

/ 到

π α

/ ， y 方向 積分範圍為 0 到 1，可得到以下六式(Model 1)： 21 11A A = (4ED 4FC + ) 22

α

α

Re − −  _(2.25a) 45B B = 2Re −  _(2.25b) 3 3 2 3 4 2 1 C = (260D +286D 630AF 286( +3) C 2457AF +1260BD 9009BD 286 Re C C 1716 9009 ) Re Re

α

α

α

α

α

α

α

α

− − − − − − 

α

(2.25c)

3 3 2 3 4 2 1 D = ( 260C 286C 1260BC 286( +3) D +9009BC +2457AE +630AE 286 Re D D 1716 9009 ) Re Re

α

α

α

α

α

α

α

α

− − − − − − 

α

(2.25d) 3 2 4 2 3 1 E = ( 180F 1638F 5103BF 234(11+ ) E E 11375AD +1750AD +234 +5148 Re Re E +57915 ) Re

α

α

α

α

α

α

α

α

− − − − −  (2.25e) 3 2 3 4 1

F = ( 180E 1638E 5103BE 234(11+ ) F F 11375AC +1750AC 234 5148 Re Re F 57915 ) Re

α

α

α

α

2

α

α

α

− − − − − − − 

α

(2.25f) 其中A B C D E F, , , , ,     為對時間微分之 Lorenz 變數 (2.25)式既為吾人所建立之 Lorenz 模型，其中包含六個非線性常 微分方程式。此組聯立方程式包含非線性項，無法直接解出，必須以 數值方法求解。為探討較多展開項對結果的影響，在附錄以相同展開 方式加上四項，以十項展開，並以相同數值方法求解，結果在第四章 探討。展開項之型式如下： 4 5 1 4 2 2 2 5 3 4 4 5 ( , , ) ( ) ( ) i x i x i x i x x y t A y B y A e A e A e A e α α α α

ψ

φ

φ

φ

φ

−

φ

φ

− ≈ + + + + + 4 5 4 4 5 5 4 4 5 5 ( ) ( ) ( ) ( ) ( ) ( )sin

( ) ( )cos ( ) ( )sin ( ) ( )cos ( ) ( )sin 2 ( ) ( )cos 2 ( ) ( )sin 2 ( ) ( )2cos A t y B t y C t y x D t y x E t y x F t y x G t y x H t y x I t y J t y x x

φ

φ

φ

α

φ

α

φ

α

φ

α

φ

α

φ

α

φ

φ

α

= + + + + + + + + +

α

(2.26)

2.5.2 Power Series 展開法 本文嘗試以不同展開方式探討普修流之現象。在管流之 方向以 第二種方式做展開。 y Power Series 之型式如下: 2 3 0 1 2 3 0 ... p p p a x a a x a x a x ∞ = = + + + +

∑

(2.27) 為一非正交函數。使用如前節之方法使之符合邊界條件，可得通式如 下: 2 1 ( ) 2 2 p p p p y y

φ

= ⎡ − + − ⎤

⎢⎣ ⎥⎦ for even l (2.28a)

3 1 2 ( ) 2 2 p p p p y y

φ

=⎡ − − + − ⎤ ⎢⎣ y⎥⎦ 2 ) 2 ) ) for odd l (2.28b) 取首三項，可得到 2 4 (y 1

φ

= − (2.29a) 2 5 y y( 1

φ

= − (2.29b) 6 2 6 (y 3y 2

φ

= − + (2.29c) 以相同於前節展開方式做展開，並以七項簡化如下： 4 5 6 1 4 2 5 ( , , )x y t A ( )y B ( )y C ( )y A ei xα A e i xα

ψ

_≈

φ

₊

φ

₊

φ

₊

φ

₊

φ

− 4 5 6 4 4 5 5 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) sin ( ) ( ) cos ( ) ( ) sin ( ) ( ) cos A t y B t y C t y D t y x E t y x F t y x G t y x

φ

φ

φ

φ

α

φ

α

φ

α

φ

= + + + + + +

α

(2.30) 將上式代入（2.14）式並使用 Galerkin 法將方程式中空間函數消 除。由於φ4、φ5、φ6並非完全之正交函數，使用 Galerkin 法後無法得 到如同(2.25)式一般在各方程式中左手邊僅有單一時間相關項，而有 一項以上之時間相關項，為使方程式能夠使用 Runge-Kutta 數值方法

運算，作如下之運算: 將原方程式化為矩陣型式如下：

[ ]

1 2 7 ( , ..., ) ( , ..., ) ( , ..., ) A R A B G R A B G B S R A B G G ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢₌ ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦   # #  ⎥ (2.31) 其中 S 矩陣為時間相關項之係數，R₁至R₇為等式右手邊相關項。將 S 矩陣作逆運算，並乘上等式右手邊之矩陣如下：

[ ]

1 1 2 7 ( , ..., ) ( , ..., ) ( , ..., ) A R A B G R A B G B S R A B G G − ⎡ ⎤ _{⎡ ⎤} ⎢ ⎥ _{⎢ ⎥} ⎢ ⎥ _{= ⎢ ⎥} ⎢ ⎥ _⎢ ⎢ ⎥ _{⎢ ⎥} ⎢ ⎥ ⎢ ⎥ _{⎣ ⎦} ⎣ ⎦   # #  ⎥ (2.32) 即可得到如同(2.25)式之型式。 運算後可得以下七式(Model 2): 1 C A = ( 608FE +608GD +40755 +6435 ) 286 −

α

α

Re Re  A _(2.33a) 45B B = 2Re −  _(2.33b) 1 C C = ( 64FE +64GD +6435 +1287 ) 78

α

α

Re Re − ₋  A _(2.33c) 3 2 3 3 3 2 4 1

D = ( 260E 286E +312AG 286(3+ ) +144CG +80AG +936CG 160BE D D D +1144BE 9009 +1716 +286 ) Re Re Re

α

α

α

α

α

α

α

α

α

α

α

− _{− −} − +  (2.33d)

3 3 2 3 3 4 2 1 E = (260D +286D 936FC +160BD 286(3+ ) 1144BD 144FC 80AF 312AF E E E 286 +1716 +9009 ) Re Re Re

α

α

α

α

α

α

α

α

α

α

− ₋ − − − − + 

_α

(2.33e) 3 3 2 3 4 2 1 F = ( 20G 182G +144EC 26(11+ )

+80AE 520AE 600EC 72BG

F F F +26 +572 +6435 ) Re Re Re

α

α

α

α

α

α

α

α

α

− _{− −} − − − 

α

(2.33f) 3 2 3 3 4 2 1 G = ( 20F 182F 600CD 26(11+ ) +80AD 520AD 72BF +144CD G G G 26 572 6435 ) Re Re Re

α

α

α

α

α

α

α

α

α

− − − − − − − − 

α

(2.33g) 在方程式左手邊各包含單一時間相關係數，此型式即可進行相關數值 運算。

第 三 章

數 值 方 法

吾人於第二章所建立的非線性方程組(2.25)式為六個一階非線性 常微分方程式所構成，無法直接求得其解，需使用數值方法以求得其 解。本章將介紹吾人所用以求解方程組之數值方法，和驗證數值解所 使用的 Newton’s method，以及觀察時間函數的方法。

3.1 Runge-Kutta Method

Runge-Kutta method 為人們用以求解微分方程組的解法之一，是 一種有效率且誤差小的高準確度數值法，為德國科學家 Runge 以及 Kutta 所創。吾人利用此法以求解在此所建立之非線性方程組。此法 之優點在於，相較於同樣用以求解微分方程組之 Euler 法，其運算速 度較快且可減少 round-off 誤差。此法為自發性之方法，在程式設計 上較為簡易，而隨著的階數不同，所推導的形式也會隨著不同的推導 過 程 有 所 不 同 。 在 此 以 六 項 [Model(1)] 展 開 為 例 ， 使 用 四 階 之 Runge-Kutta method 求取數值解： A、B、C、D、E、F 之函數形式如下： A = f A B C D E F t1( , , , , , , ) B = f2( , , ,A B C D E F t, , , ) C = f A B C D E F t3( , , , , , , ) D = f4( , , ,A B C D E F t, , , ) E = f A B C D E F t5( , , , , , , ) F = f6( , , ,A B C D E F t, , , )

運算方式如下: 11 1( , , , , , , ) k =hf A B C D E F t k21 =hf A B C D E F t2( , , , , , , ) k₃₁=hf A B C D E F t₃( , , , , , , ) k41 =hf A B C D E F t4( , , , , , , ) k₅₁=hf A B C D E F t₅( , , , , , , ) k61 =hf A B C D E F t6( , , , , , , ) k₁₂ =hf A k₁( + ₁₁/ 2,B+k₂₁/ 2,C+k₃₁/ 2,D+k₄₁/ 2,E+k₅₁/ 2,F+k₆₁/ 2,t+h/ 2) 22 2( 11/ 2, 21/ 2, 31/ 2, 41/ 2, 51/ 2, 61/ 2, / 2) k =hf A k+ B+k C+k D+k E+k F+k t+h 32 3( 11/ 2, 21/ 2, 31/ 2, 41/ 2, 51/ 2, 61/ 2, / 2) k =hf A k+ B+k C+k D+k E+k F+k t+h 42 4( 11/ 2, 21/ 2, 31/ 2, 41/ 2, 51/ 2, 61/ 2, / 2) k =hf A k+ B+k C+k D+k E+k F+k t+h 52 5( 11/ 2, 21/ 2, 31/ 2, 41/ 2, 51/ 2, 61/ 2, / 2) k =hf A k+ B+k C+k D+k E+k F+k t+h k₆₂ =hf A k₆( + ₁₁/ 2,B+k₂₁/ 2,C+k₃₁/ 2,D+k₄₁/ 2,E+k₅₁/ 2,F+k₆₁/ 2,t+h/2) 2) 2) 2) 2) ) ) ) ) ) ) 13 1( 12/ 2, 22/ 2, 32/ 2, 42/ 2, 52/ 2, 62/ 2, / 2) k =hf A k+ B+k C+k D+k E+k F+k t+h 23 2( 12/ 2, 22/ 2, 32/ 2, 42/ 2, 52/ 2, 62/ 2, / k =hf A k+ B+k C+k D+k E+k F+k t+h 33 3( 12/ 2, 22/ 2, 32/ 2, 42/ 2, 52/ 2, 62/ 2, / k =hf A k+ B+k C+k D+k E+k F+k t+h 43 4( 12/ 2, 22/ 2, 32/ 2, 42/ 2, 52/ 2, 62/ 2, / k =hf A k+ B+k C+k D+k E+k F+k t+h 53 5( 12/ 2, 22/ 2, 32/ 2, 42/ 2, 52/ 2, 62/ 2, / k =hf A k+ B+k C+k D+k E+k F+k t+h 63 6( 12/ 2, 22/ 2, 32/ 2, 42/ 2, 52/ 2, 62/ 2, / 2) k =hf A k+ B+k C+k D+k E+k F+k t+h k₁₄ =hf A k₁( + ₁₃,B+k₂₃,C+k₃₃,D+k₄₃,E+k₅₃,F+k₆₃,t+h k₂₄ =hf A k₂( + ₁₃,B+k₂₃,C+k₃₃,D+k₄₃,E+k₅₃,F+k₆₃,t+h k₃₄ =hf A k₃( + ₁₃,B+k₂₃,C+k₃₃,D+k₄₃,E+k₅₃,F+k₆₃,t+h k₄₄ =hf A k₄( + ₁₃,B+k₂₃,C+k₃₃,D+k₄₃,E+k₅₃,F+k₆₃,t+h k₅₄ =hf A k₅( + ₁₃,B+k₂₃,C+k₃₃,D+k₄₃,E+k₅₃,F+k₆₃,t+h k₆₄ =hf A k₆( + ₁₃,B+k₂₃,C+k₃₃,D+k₄₃,E+k₅₃,F+k₆₃,t+h 其中kij為運算變數， 為時間間隔，可藉由 之改變來控制精確度。 h h 解可獲得如下：

A1,n+1 = A1,n +(k11+2k12+2k13+k14) /6 6 6 6 6 6 B1,n+1 =B1,n +(k21+2k22+2k23+k24) / C1,n+1 =C1,n +(k31+2k32+2k33+k34) / D_1,n+1 = D_1,n +(k₃₁+2k₃₂+2k₃₃+k₃₄) / E_1,n+1 = E_1,n +(k₃₁+2k₃₂+2k₃₃+k₃₄) / F_1,n+1 = F_1,n +(k₃₁+2k₃₂+2k₃₃+k₃₄) /

此即為 Lorenz 方程式之解。對於十項展開[Model(3)]以及 Power Series 所得到之方程式(2.33) [Model(2)]同樣以此法求解，不同之處僅為方程 式多寡。

3.2 Newton’s Method

吾人採用 Newton’s Method 以驗證由 Runge-Kutta 求得之 steady state 解。以六項展開式為例， f A B C D E F t_n( , , , , , , )=0為六個非線性方程 式所組成的系統，n=1~6，Newton’s Method 之疊代方法如下： 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [ , , , , , [ , , , , , ] k k k k k k k k k k k k k k k k k k k k k k k k A A B B C C J A B C D E F D D E E F F f A B C D E F − − − − − − − − − − − − − − − − − − − ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢₌ ⎥₋ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ i ] 其中k=1,2,3...，為疊代次數，J為 Jacobian 矩陣，其定義如下：

1 1 1 2 2 2 6 6 6 ... ... ... f f f A B F f f f J A B F f f f A B F ∂ ∂ ∂ ⎡ ⎤ ⎢_{∂ ∂ ∂} ⎥ ⎢ ⎥ ∂ ∂ ∂ ⎢ ⎥ ⎢ ⎥ = ∂ ∂ ∂ ⎢ ⎥ ⎢ ⎥ ⎢_{∂ ∂ ∂} ⎥ ⎢_{∂ ∂ ∂} ⎥ ⎣ ⎦ 給定適當的初始值開始疊代，Newton’s Method 即能夠在較少的時間 格點求得 steady state 之解。

3.3 時間級數圖

動力系統可藉由其時間函數 f(t)，或是其動力變數之時間級數 (Time series)來表示此系統隨時間演變的過程。流場在低臨界雷諾數 時，其時間對於時間函數的作圖會呈現快速收斂情形；在雷諾數提 高，流場流動行為變為劇烈，收斂至穩定值的時間也隨之增加。在流 場雷諾數高於臨界雷諾數時，時間級數圖一般而言則不會發生收斂的 情形，取而代之為週期性或擬週期性，甚至非週期性的流場。

3.4 線性化與求取特徵值

將 Model(1)、(2)、(3)中的非線性項消去，並將等式右手邊的時 間函數化為矩陣型式。吾人採用 Matlab 中所內建的”eig”指令求取線 性化方程組矩陣之特徵值，以判斷臨界雷諾數值。

第 四 章

結 果 與 討 論

4.1 以時間級數圖觀察模型行為

吾人將三組模型以 Runge-Kutta 方式求解，將每組模型取出兩個 時間函數對時間作圖，得到圖形如圖(二)至圖(七)所示。各圖中分別 取五組雷諾數，並將波常數固定在α=1，比較模型在不同雷諾數所表 現的行為。 Model(1)中取 E、F 兩時間函數作圖，如圖(二)、(三)所示，在 雷諾數為 500 時，時間函數在短時間即達到一固定收斂值，與Fortin( 1994)的結果符合。當雷諾數提高至 5000 時，時間函數收斂至固定值 的時間增加，且時間函數振盪的幅度較低雷諾數時大。在物理意義上 提高雷諾數代表流場的流動情形變得更劇烈，流場需要更長時間達到 穩定，因此時間函數收斂的時間增加與實際物理意義符合。 吾人將雷諾數提高至 Orszag 由線性理論所得到的臨界雷諾數值 ( =5772.22)之上，即 Re=6000，試圖觀察流場是否出現與 Fortin 文 中符合的擬週期性的流動情形。而所得到的結果顯示時間函數僅較雷 諾數為 5000 時需要更長的時間達到收斂值，所得到時間函數的解與 Fortin 所解出之 limit cycle 解不符合，意即流場在經過足夠長的時間 依舊會達到穩定，並未產生不穩定。此與實際的物理現象不符合，吾 人再將雷諾數提高至 10000，觀察是否為模型的缺陷導致無法得到合 理的物理現象。結果顯示雷諾數為 10000 時流場依然為穩定，僅發現 達到穩定的時間向後延長。為了數學上的探討，再將雷諾數提高至 50000，結果同樣發現時間函數在足夠的時間內趨向定值，而收斂時 c R

間再向後延長。 在 Model(1)無法得到預期的結果後，吾人嘗試以同樣展開方式並 增加展開項數至十項，得到 Model(3)，如附錄所示，探究是否因為展 開項數數量不足，而無法得到合理的物理現象。在展開項中取出 E、 F 兩項時間級數作圖，如圖(六)、(七)所示。所得到圖形與圖(二)、(三) 雷同，低雷諾數的條件下時間函數在較短的時間達到收斂，而隨著雷 諾數增加，時間函數收斂的時間亦加長，並未產生不穩定流場，顯示 以 Legendre 正交函數展開無法得到所預期能與物理現象符合的結果。

吾人嘗試不同的 y 方向展開式 Power Series 做展開[Model(2)]， 所得到方程式為(2.33)式。取出時間函數 F、G 作圖，如圖(四)、(五) 所示。結果與 Model(1)、(3)相似，改變 y 方向展開式尚無法得到與 物理現象符合的結果。

4.2 以線性化模型判斷解的正確性

吾人以 Runge-Kutta 法無法得到預期的結果，為驗證所建立的模 型在數學上的正確性，吾人將三組模型線性化，意即將非線性項去 除，並將方程式各時間函數化為矩陣型式，利用 Matlab 內建的 ”eig”指令求取矩陣的特徵值值，並改變雷諾數以觀察線性模型發生不 穩定解的臨界雷諾數所在。在第一章已介紹 Orszag(1971)利用此法得 到線性方程組最精確的解。吾人利用該方法判斷結果是否與 4.1 節符 合。

表(一)、(二)、(三)分別代表線性化之 Model(1)、(2)、(3)在不同 雷諾數下所求得之特徵值值。在表(一)中，雷諾數為 10 時可發現 求 得 特 徵 值 之 實 部 值 皆 為 負 值 ， 代 表 各 時 間 函 數 隨 時 間 衰 減。提高雷諾數至 100、1000、5000 等，此時發現特徵值之實部值仍 為負值，但其值逐漸增大接近 0，代表時間函數隨時間衰減情形隨著 雷諾數增大而減緩，意即流場的黏滯力影響逐漸小於慣性力，因而需 要更長的時間達到穩定，與物理現象符合。 將雷諾數提高至 10000，此值已遠高於 Orszag 所解得之臨界雷諾 數，發現特徵值實部仍為負值，並未轉變成表示流場已進入不穩定解 的正值。單純以數學上的觀點，將雷諾數提高至 50000、100000 甚至 更高，觀察特徵值實部值是否變為正。結果即使在雷諾數 500000 時， 此值僅趨近於 0，代表在有限的時間內流場仍會趨向穩定。 觀察表(二)、(三)解之型態，與表(一)的趨勢相似，並未發現特徵 值實部為正值的情形。此結果雖然與物理現象不符合，在模型中無法 找到臨界雷諾數，發現流場轉變為不穩定之處。但三組結果與吾人在 4.1 所得到的結果符合，驗證了這三組模型無法在超過臨界雷諾數之 後得到非零解。

4.3 牛頓法的驗證

吾人將方程組中時間函數去除，以牛頓法求解模型的 steady-state 解。可發現 Model(1)、(2)、(3)在不同雷諾數的情形下，皆得到 0 解， 並無其他組非零解產生(如 Lorenz Equation 在提高瑞里數產生不同解 的現象)。此結果與 4.1 所得到的結果一致，可驗證三組模型皆僅有一 組解產生。

4.4 超過臨界值無法獲得非零解的可能原因

吾人所建立的三組模型所得到結果皆顯示，在臨界雷諾數以下流 場能在相對較短時間內達到穩定，提高流場雷諾數則需要相對較長時 間才能達到穩定，與實際物理現象以及前人文獻結果符合。但在超過 臨界雷諾數的情形下，則無法以模型得到符合物理現象的解，亦與前 人文獻結果不符。吾人利用 4.2 以及 4.3 兩種不同驗證方法證明 4.1 節所得到結果的正確性，推斷無法以模型得到合理現象的解釋如下： 1.由使用 Legendre 函數做 y 方向展開所建立的 Model(1)以及Model(3) 所得到相近的結果顯示，在少量(6、10 項)展開項數的情形下皆無法 得到合理的結果。由於再增加ｘ方向展開項運算更加複雜，本文中 尚未完成繼續增加展開項數的研究。而依 Dowell 僅使用一組正、餘 弦函數的 x 方向展開模式加上大量不同的 y 方向展開能夠得到非線 性現象，x 方向展開項數應為足夠。吾人亦曾嘗試僅使用正弦或是 餘弦單一函數對 x 方向做傅立葉展開，然而發現會在正交積分的過 程會導致基態流場消去，因此一組模式必須同時使用正、餘弦函數 做 x 方向展開。

2.使用 Power Series 做 y 方向展開所建立的 Model(2)亦無法描述正確 的物理現象。吾人在建立三組模型時，發現 y 方向的展開模式若選 擇為偶函數，則在正交積分的過程會導致非線性項相互消去，形成 線性的方成組，因此欲探討非線性的問題，必須搭配奇函數做展開， 與 Dowell 的研究吻合。而吾人使用 Power Series 前兩項與 Legendre 函數相似，僅係數不同，所加上的第三項 y 方向函數非奇非偶，由 結果顯示並無特別作用。而 Model(1)~(3)雖然包含非線性項，卻並 未發生非線性的作用。單就數學上與其他非線性聯立方程式比較，

在某些條件下非線性作用未發生是合理的，原因可能來自於模型建 立的方式，包含 y 方向展開項數的不足以及展開函數的選擇。此方 向的研究，可以嘗試在在 y 方向使用不同的函數，或是 y 方向使用 更多展開項數，而在 x 方向同樣使用傅立葉展開，或許能在超越臨 界雷諾數的條件下，引發擾動方程式中非線性項的作用。 3.吾人嘗試利用與 Lorenz 相同的方式建立動力模型，但無法與 Lo -renz 一般，僅利用極少展開項即描繪出合理的物理現象。由前人文 獻可知，使用 Chebyshev 以及傅立葉函數對擾動方程式做雙重展 開，以及有限元素法仍為解普修流問題的較佳方式。利用動力模型 的方法雖可簡化計算過程，但準確性及是否能正確描繪物理現象在 本文沒有正面的答案。因此最後建議，使用如第二點的方式，期盼 能夠使此模型符合物理現象。

第 五 章

結 論

本文藉由 Lorenz 模型，輔以 Runge-Kutta 數值方法求解，再利用 時間級數圖觀察平面普修流流場的行為，綜合前一章的結果以及驗證 方式，吾人得到以下的結論： 1. Model(1)、(2)、(3)在低雷諾數以及臨界雷諾數( =5772.22)以下的 流場皆能與 Fortin(1994)的結果吻合，在較短的時間皆能收斂至固 定解，達到穩定的流場。 c R 2. Model(1)、(2)、(3)在超過臨界雷諾數皆未能得到符合物理現象的 非零解，在求得線性化模型的特徵值皆為負值，以及牛頓法解得 穩態解皆為零解的驗證下，證明模型的數學上的正確性。單就數 學上而言，某些情形下非線性作用未發生是合理的，原因可能來 自模型建立的方式。 3. 超過臨界值無法得到合理現象的原因可能為展開項數不足以及展 開函數的選擇不適合，在 x 方向展開方式及項數與 Dowell 相似甚 至更多的合理條件下，本文的改進可嘗試以不同的 y 函數做展開， 以及增加 y 方向展開項數，或許能在臨界雷諾數以上引發擾動方 程式中非線性項的作用。

參 考 文 獻

Baker, G. L. and Gollub, J. P. 1990 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, O. C., 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.

Donnelly, R. J., Reif, F. and Suhl, H. 1962 Enhancement of hydrodynamic stability by modulation. Phys. Rev. Lett. 9, 363-365.

Fortin A. and Jardak M., 1994 Old and new results on the two-dimensional Poiseuille flow. J. Computational Physics. 115, 455

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

Grosch, C. E. and Salwen, H., 1968 The stability of steady and time-dependent plane Poiseuille flow. J. Fluid Mech, 34, 177

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

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

Lessen, M., Fox, J. A., Bhat W. V. & Liu, T. Y. 1964 Phys. Fluids, 7, 1384.

Lin, C. C., 1955 The Theory of Hydrodynamic Stability. Cambridge

Monographs On Mechanics And Applied Mathematics

Lorenz, E. N. 1963 Deterministic nonperiodic flow. J. Atmos. Sci. 20, 130-141.

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 Ng,Bart S. 2002 On the spectral problem for Poiseuille flow in a circular pipe. Fluid Dynamics Research 33, Issue: 1-2

Shen, S. F. 1954 Calculated amplified oscillations in plane Poiseuille and Blasius flows. 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. Phys. Rev. 91, 780

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

Applied Mechanics. 45, 13-18

李鎰清 2001 圓柱之流場-一種改良 Lorenz 模型.

國立交通大學

碩士論文

.

周煒超 2001 旋轉圓柱間流場之動力模型.

國立交通大學碩士論

Re=10 Re=100 Re=1000 Re=5000

A –2.2542 + 0.6474i –0.2254 + 0.6474i –0.0225 + 0.6474i –0.0045 + 0.6474i B –2.2542 – 0.6474i –0.2254 – 0.6474i –0.0225 – 0.6474i –0.0045 – 0.6474i C –0.9625 + 0.4773i –0.0963 + 0.4773i –0.0096 + 0.4773i –0.0019 + 0.4773i D –0.9625 – 0.4773i –0.0963 – 0.4773i –0.0096 – 0.4773i –0.0019 – 0.4773i E –1.05 –0.105 –0.0105 –0.0021

F –2.25 –0.225 –0.0225 –0.0045

Re=10000 Re=50000 Re=100000 Re=500000 A –0.0023 + 0.6474i –0.0005 + 0.6474i –0.0002 + 0.6474i –0.0000 + 0.6474i B –0.0023 – 0.6474i –0.0005 – 0.6474i –0.0002 – 0.6474i –0.0000 – 0.6474i C –0.0010 + 0.4773i –0.0002 + 0.4773i –0.0001 + 0.4773i –0.0000 + 0.4773i D –0.0010 – 0.4773i –0.0002 – 0.4773i –0.0001 – 0.4773i –0.0000 – 0.4773i E –0.001 –0.0002 –0.0001 0 F –0.0022 –0.0004 –0.0002 0

Re=10 Re=100 Re=1000 Re=5000

A –0.9875 –0.0988 –0.0099 –0.002 B –5.0125 –0.5012 –0.0501 –0.01 C –2.2542 + 0.6474i –0.2254 + 0.6474i –0.0225 + 0.6474i –0.0045 + 0.6474i D –2.2542 – 0.6474i –0.2254 – 0.6474i –0.0225 – 0.6474i –0.0045 – 0.6474i E –0.9625 + 0.4773i –0.0962 + 0.4773i –0.0096 + 0.4773i –0.0019 + 0.4773i F –0.9625 – 0.4773i –0.0962 – 0.4773i –0.0096 – 0.4773i –0.0019 – 0.4773i G –2.25 –0.225 –0.0225 –0.0045

Re=10000 Re=50000 Re=100000 Re=500000 A –0.001 –0.0002 –0.0001 0

B –0.005 –0.001 –0.0005 –0.0001 C –0.0023 + 0.6474i –0.0005 + 0.6474i –0.0002 + 0.6474i –0.0000 + 0.6474i D –0.0023 – 0.6474i –0.0005 – 0.6474i –0.0002 – 0.6474i –0.0000 – 0.6474i E –0.0010 + 0.4773i –0.0002 + 0.4773i –0.0001 + 0.4773i –0.0000 + 0.4773i F –0.0010 – 0.4773i –0.0002 – 0.4773i –0.0001 – 0.4773i –0.0000 – 0.4773i G –0.0022 –0.0004 –0.0002 0

Re=10 Re=100 Re=1000 Re=5000

A –2.3433 + 1.3436i –0.2343 + 1.3436i –0.0234 + 1.3436i –0.0047 + 1.3436i B –2.3433 – 1.3436i –0.2343 – 1.3436i –0.0234 – 1.3436i –0.0047 – 1.3436i C –0.9625 + 0.4773i –0.0963 + 0.4773i –0.0096 + 0.4773i –0.0019 + 0.4773i D –0.9625 – 0.4773i –0.0963 – 0.4773i –0.0096 – 0.4773i –0.0019 – 0.4773i E –2.2542 + 0.6474i –0.2254 + 0.6474i –0.0225 + 0.6474i –0.0045 + 0.6474i F –2.2542 – 0.6474i –0.2254 – 0.6474i –0.0225 – 0.6474i –0.0045 – 0.6474i G –1.0214 + 1.3247i –0.1021 + 1.3247i –0.0102 + 1.3247i –0.0020 + 1.3247i H –1.0214 – 1.3247i –0.1021 – 1.3247i –0.0102 – 1.3247i –0.0020 – 1.3247i I –1.05 –0.105 –0.0105 –0.0021

J –2.25 –0.225 –0.0225 –0.0045

Re=10000 Re=50000 Re=100000 Re=500000 A –0.0023 + 1.3436i –0.0005 + 1.3436i –0.0002 + 1.3436i –0.0000 + 1.3436i B –0.0023 – 1.3436i –0.0005 – 1.3436i –0.0002 – 1.3436i –0.0000 – 1.3436i C –0.0010 + 0.4773i –0.0002 + 0.4773i –0.0001 + 0.4773i –0.0000 + 0.4773i D –0.0010 – 0.4773i –0.0002 – 0.4773i –0.0001 – 0.4773i –0.0000 – 0.4773i E –0.0023 + 0.6474i –0.0005 + 0.6474i –0.0002 + 0.6474i –0.0000 + 0.6474i F –0.0023 – 0.6474i –0.0005 – 0.6474i –0.0002 – 0.6474i –0.0000 – 0.6474i G –0.0010 + 1.3247i –0.0002 + 1.3247i –0.0001 + 1.3247i –0.0000 + 1.3247i H –0.0010 – 1.3247i –0.0002 – 1.3247i –0.0001 – 1.3247i –0.0000 – 1.3247i I –0.001 –0.0002 –0.0001 0

J –0.0022 –0.0004 –0.0002 0

(A) (B)

(C) (D)

(E)

圖二 Model(1)，函數 E 之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

(A) (B)

(C) (D)

(E)

圖三 Model(1)，函數 F 之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

(A) (B)

(C) (D)

(E)

圖四 Model(2)，函數 F 之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

(A) (B)

(C) (D)

(E)

圖五 Model(2)，函數 G 之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

(A) (B)

(C) (D)

(E)

圖六 Model(3)，函數 E 之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

(A) (B)

(C) (D)

(E)

圖七 Model(3)，函數 F 之時間級數圖 (A)Re=500 (B)Re=5000 (C)Re=6000 (D)Re=10000 (E)Re=50000

附 錄

吾人利用 Legendre 函數建立之改良 Lorenz 方程組為十個一階非 線性常微方程式(Model 3)如下： 21 11A A = ( 4FC +8IH +4ED 8JG + ) 22 α α α α Re − − − 45B B = 2Re − 3 3 3 2 3 3 2 4 1 C = (572D +520D +2520 BD 18018BD +5670 FH 572(3+ ) 15561FH +4095JD +5670 EG+4095IC 15561EG 1260 AF C C C 4914AF 18018 3432 572 ) Re Re Re α α α α α α α α α α α α α α α − − − − − − − − 3 3 2 3 3 2 4 1 D = (520C +572C 4914AE 1260AE +2520BC 572(3+ ) 18018BC 15561FG +5670FG 4095JC 5670EH +15561EH D D D +4095ID 18018 +3432 +572 ) Re Re Re α α α α α α 3 α α α α α α α α − − − − − − − + α 3 3 2 3 3 4 2 1 E = (3276F +360 F 3500AD +22750AD +10206BF 468( +11) 15750HD +11375HD 15750CG +5103JF +11375GC +5103IE E E E 468 10296 115830 ) Re Re Re α α α α α α α α α α α α α − − − − − − α 3 3 2 3 3 4 2 1

F = (3276E +360 E+22750AC +10206BE 3500AC 468( +11) +11375GD 11375HC 5103EJ 15750GD +15750HC +5103IF F F F +468 +10296 +115830 ) Re Re Re α α α α α α α α α α α α α − ₋ − − − α

3 2 3 3 4 2 1 G = ( 5733EC +5733FD 2080H +18018BH 286(3+4 ) 572H 10080BH +5040AJ +4914AJ G G G +4576 +6864 9009 ) Re Re Re α α α α α α α α α α − − − − − + α 3 2 3 3 4 2 1

H = (4914AI +5040AI +5733FC +5733ED +18018BG 572G 286(3+4 ) H H H 2080G 10080BG 4576 6864 9009 ) Re Re Re α α α α α α α α α α − − − − − − α 2 2 2 2 3 2 3 4 2 1 I = (11375D 5103E +5103F 11375C +2880J +6552J 468(4 +11) I I +20412BJ +45500AH 28000AH 7488 41184 115830 ) Re Re Re α α α α α I α α α α α α α − − − − − − 3 2 3 4 2 1

J = (1440I +3276I +11375CD +5103FE +10206BI +22750AG 234(4 +11) J J J 14000AG +3744 +20592 +57915 ) Re Re Re α α α α α α α α α − − α