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

國立交通大學

機械工程研究所

碩 士 論 文

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

Nonlinear Instability of Plane Poiseuille Flow with

Modulated Pressure Gradient

研究生：廖啟豪

指導教授：楊文美 博士

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

Nonlinear Instability of Plane Poiseuille Flow with Modulated

Pressure Gradient

研 究 生：廖 啟 豪 Student：Chi-Hao Liao 指導教授：楊 文 美 Adviser：Wen-Mei Yang 國 立 交 通 大 學 機 械 工 程 學 系 碩 士 論 文 A Thesis

Submitted to Institute of Mechanical Engineering College of Engineering

National Chiao Tung University in Partial Fulfillment of The Requirements

for the Degree of Master of Science

in

Mechanical Engineering July 2006

Hsinchu, Taiwan, Republic in China

摘 要

本 文 以 非 線 性 方 法 探 討 調 制 平 面 普 修 流 的 不 穩 定 現 象 ， 首 先 建 立 流 線 擾 動 方 程 式，然 後 使 用 頻 譜 法 (Spectral method) 以 雙 重 級 數 展 開 來 模 擬 流 線 函 數 ， 基 底 函 數 在 x 方 向 採 用 複 數 形 式 的 Fourier 級 數，在 y 方 向 選 用 具 有 快 速 收 歛 特 性 及 高 精 準 度 的 Chebyshev 多 項 式 。 在 數 值 方 法 部 分 ， 線 性 項 中 具 有 基 態 (basic state) 速 度 者 及 非 線 性 項 採 用 Adam-Bashforth Method 離 散 ， 其 餘 則 採 用 Crank-Nicolon Method 離 散 ， 然 後 將 雙 重 級 數 代 入 離 散 化 的 流 線 擾 動 方 程 式 後 ， 在 x方 向 以

Galerkin method 處 理 ， y方 向 以 定 置 法 (Collocation method)

處 理 ， 使 方 程 式 得 以 轉 化 成 代 數 方 程 式 並 可 逐 步 求 解 。 找 出 不 同 調 制 條 件 下 的 臨 界 雷 諾 數 ， 並 將 數 據 繪 製 成 偏 移 量 (threshold shift)圖，藉 此 探 討 調 制 作 用 對 流 場 臨 界 雷 諾 數 的 影 響。利 用 時 間 級 數 圖、FFT 頻 譜 圖、相 平 面 圖 及 Poincare 映 射 圖 探 討 超 越 臨 界 雷 諾 數 之 後 的 調 制 流 場 行 為 。 結 果 發 現 ， 在 低 頻 時 調 制 會 造 成 明 顯 的 負 偏 移 ， 而 在 中 至 高 頻 間 調 制 會 造 成 微 小 的 正 偏 移 ， 而 調 制 強 度 會 放 大 偏 移 的 效 果 。 對 超 越 臨 界 雷 諾 數 之 後 的 調 制 流 場 而 言 ， 隨 著 雷 諾 數 增 加 流 場 會 從 原 本 的 單 一 週 期 性 運 動 ， 逐 漸 演 變 至 最 終 的 混 沌 (Chaos) 狀 態 ， 而 調 制 作 用 所 造 成 的 影 響 為 ： 增 加 調 制 頻 率 將 造 成 流 場 變 動 加 快 提 前 進 入 混 沌 ； 增 加 調 制 振 幅 會 使 得 流 場 容 易 分 配 能 量 至 與 調 制 頻 率 無 關 的 頻 率 ， 導 致 流 場 提 前 進 入 混 沌 狀 態 。

誌謝

本論文可以順利完成，在此感謝指導教授楊文美博士的悉心指導 與督促，特此致上感激與謝意。 承蒙口試委員胡毓仁博士、盧定昶博士與游明輝博士於口試期 間，給予論文指正與寶貴意見，使得本論文更為嚴謹與完善，在此亦 表達由衷的感謝。 本研究得以順利完成，感謝同窗好友秋傑，以及學長豪傑、和學 弟文山、武伸於課業上及生活上的互相提攜，使得研究生涯更為充實。 最後，僅以本論文成果獻給我親愛的家人，尤其是我的母親，感 謝你們的關心與鼓勵。

目錄

摘要………..…………i 誌謝………..…………...ii 目錄………iii 表目錄……….... v 圖目錄………vi 符號說明………..viii 第一章 緒論……….1 1. 1 文獻回顧……….………..1 1. 2 研究目的……..………...4 第二章 數學模式……….5 2. 1 統御方程式………...5 2. 2 無因次化………...6 2. 3 基態解………...6 2. 4 擾動方程式之建立………...7 第三章 數值方法………...10 3. 1 頻譜法………...10 3.1.1 雙重級數展開………..………... 10 3.1.2 離散方程式………...………11 3.1.3 代數方程式………...………11 3. 2 數值解的分析法……...……….……….13 3.2.1 時間級數與FFT頻譜分析...……….…………. ……..…13 3.2.2 相平面圖………...14 3.2.3 Poincare映射圖………...15

第四章 結果與討論………...16 4. 1 非調制流場的穩定性……….…………16 4. 2 調制作用對臨界雷諾數的影響………..……..…….…....17 4. 3 臨界雷諾數後之調制流場行為………..…….…..17 4.3.1 調制頻率的影響………17 4.3.2 調制振幅的影響………20 結論………...22 參考文獻………...23

表 目 錄

表一 調制振幅(a)ε =1和(b)ε =0.5時，在不同調制頻率下所對應的

圖 目 錄

圖一 調制普修流物理模式示意圖…………...27 圖二 當M =3，α =1.020545時，針對非調制流場(ε =0)不同N所得 到的臨界雷諾數。………. 28 圖三 當N =43，α =1.020545時，針對非調制流場(ε =0)不同M 所 得到的臨界雷諾數。………..29 圖四 偏移量圖………..30 圖五 當ε =1，ω =0.1，Re=6000時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...31 圖六 當ε =1，ω =0.1，Re=6500時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。……….. 33 圖七 當ε =1，ω =0.1，Re=8500時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...35 圖八 當ε =1，ω =0.1，Re=10000時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...37 圖九 當ε =1，ω =1，Re=6000時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...39 圖十 當ε =1，ω =1，Re=6500時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...41 圖十一 當ε =1，ω =1，Re=8500時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...43 圖十二 當ε =1，ω =1，Re=10000時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...45 圖十三 當ε =1，ω =10，Re=6000時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...47 圖十四 當ε =1，ω =10，Re=6500時的(a)時間級數圖、(b)頻譜圖、

(c)相平面圖及(d)Poincare 映射圖。………...49 圖十五 當ε =1，ω =10，Re=8500時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...51 圖十六 當ε =1，ω =10，Re=10000時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...53 圖十七 當ε =0.5，ω =1，Re=6500時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...55 圖十八 當ε =0.5，ω =1，Re=8500時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...57 圖十九 當ε =0.5，ω =1，Re=10000時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...59 圖二十 當ε =0.5，ω =1，Re=11000時的(a)時間級數圖、(b)頻譜圖、 (c)相平面圖及(d)Poincare 映射圖。………...61

符 號 說 明

符 號

mn a 時間係數 H 平行板之間距 N M , 流線函數展開項之項數 P 壓力 Re 雷諾數 c Re 臨界雷諾數 0 Re 非調制臨界雷諾數 t 時間 n T n階Chebyshev多項式 v u, 速度分量 0 u 非調制時流道中央處最大速度 y x, 平面座標 α 波數 Δ 2 2 2 2 y x ∂ ∂ + ∂ ∂ 或偏移量

(

Rec−Re0

)

/Re0 2 Δ 4 4 2 2 4 4 4 2 y y x x ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ε 調制振幅 μ 黏滯係數 ρ 密度 φ y方向展開函數 ψ 流線函數 ω 調制頻率

上標說明

基態

‘ 擾動量

第 一 章

緒 論

1.1 文獻回顧

典型的普修流(Poiseuille flow)為藉由一固定壓力梯度驅使流體 流動之流場。由於普修流是流體力學中最基本的問題之ㄧ，在考慮穩 定層流流動的前提下，經由Navier-Stokes 方程式可得速度的單一函數 基態解，故當學者們一旦發展出新的實驗或數值方法，普修流便成為 驗證新方法正確性的最佳流場形式。 依照流體所流經空間的幾何特性，普修流大致可被區分為兩種類 型，分別是流經等截面圓管普修(Circular Poiseuille flow)流及流經兩 固定平行板間的平面普修流(Plane Poiseuille flow)。關於普修流穩定性 的研究及其沿革也可依此劃分為兩個研究方向介紹如下。 1.圓管普修流 可追朔至 19 世紀末，Osborne Reynolds 以圓導管做了普修流的實 驗，發現層流在雷諾數為13000 時發生不穩定。後人更進一步考慮實 驗物表面的差異來作研究，發現在表面非常粗糙時臨界(critical)雷諾 數約為2000，表面非常平滑時臨界雷諾數可達 40000。 Sexl(1927) 首先將圓導管普修流穩定性的問題理論化，推導出相 關的方程式。之後Sexl 和 Lessen,Fox, Bhat 以及 Liu(1964)更深入探討 複雜的非對稱管流的問題。在實驗以及理論對照下，他們發現了普修 流在受到極微小的擾動下呈穩定，但在一般的擾動下則呈現不穩定的 狀態。

2.平面普修流－Orr-Sommerfeld 方程式

W. Orr & A. Sommerfeld(1907)兩人針對平行通道(channel)內之普 修流提出流線函數的線性擾動方程式用以計算流場的穩定性，稱之為

Orr-Sommerfeld 方程式。時至 1950-1960 年代，隨著計算機科技初步 發展，學者試著以各種數值方法針對 Orr-Sommerfeld 方程式作數值 解，藉此探討線性普修流的問題求得發生不穩定的臨界雷諾數。 Shen(1955)以 Lin(1954)所發表的近似分析(asymsptotic analysis)法得 到 Re₀ =5360 ， α =1.05 。Thomas(1953) 率 先 以 有 限 元 素 法 得 到

5780

Re₀ = ，α =1.026；Betchov & Criminale(1967)隨後也以有限元素法 得到Re₀ =5360，α =1.05。 Orszag(1971)以 Chebyshev 多項式搭配 QR 矩陣特徵值演算法得 到臨界雷諾數發生於Re₀ =5772.22，α =1.02056。他以方程式中 Im(λ) 之正負值判斷解的穩定與否。當雷諾數不大於 5772.22 時，Im(λ)皆 為負值，一旦超過5772.22 則 Im(λ)當中會有正值出現。Orszag 的結 果不但證明Chebyshev 多項式具有高度精確及快速收斂的特性，同時 可節省大量CPU 時間，相當適合用來處理流體運動方面的問題。 Fortin 與 Jardak(1994)針對普修流的穩定性分析進行了完整的研 究。他們將問題分細分為一維及二維、線性及非線性，首先解出最基 本 的 線 性 Orr-Sommerfeld 方 程 式 ， 得 到 臨 界 雷 諾 數 發 生 於 020545 . 1 = α ，Re₀ =5772.22，與 Orszag 的結果吻合。接著將非線性加 入Orr-Sommerfeld 方程式中，採用有限元素法直接解出整個非線性的 擾動方程式，同時將完整的流場繪出。結果顯示在低雷諾數時流場呈 現穩定的流動，外加擾動在一定的時間內達到收斂，當雷諾數提高到 臨界值時，時間級數圖顯示擾動並不隨時間收斂，而是在一定的範圍 內震盪。進一步利用 FFT(Fast Fourier Transform)頻譜圖及相平面圖 (Phase plane)分析得知流動行為僅有一組基頻與其倍數頻率的組合， 為單一週期性的運動，稱為週期性流場。為了瞭解流場超過臨界雷諾 數時的行為，再繼續提高雷諾數，隨著雷諾數的增加，流場接著出現 具有第二組頻率的擬週期性(quasi-periodic)運動，直到雷諾數達到 12900 時，流場已產生混沌現象。

洪英棋(2005)採用頻譜法(Spectral method)將流線函數以雙重級 數展開，在水平方向以複數形式的 Fourier 級數展開，垂直方向以 Chebyshev 多項式展開，而展開係數則為時間的函數，藉此模擬流場 的變化，得到與Fortin 一致的結果。 近年來，調制形式的問題已開始受到重視並逐漸發展。藉著調制 作用導致流場提前發生不穩定，促使流場中的熱量或質量加速傳遞， 例如醫學上可藉由調制作用使得混合配方的時間顯著減少；相對的藉 著調制作用使得流場穩定化，機械潤滑更為穩定順利，或是使印表機 噴墨更為均勻。關於調制問題的研究，調制旋轉圓柱間之流場堪為代 表。

Donnelly et al. (1962)和 Donnelly(1964)以實驗的方式研究流場受 到調制作用的影響，探討內圓柱在不同的調制頻率及調制振幅轉速下 對整個流場穩定性的影響。為了解調制作用對流場臨界雷諾數的影 響，文中定義偏移量(Threshold shift)Δ=

(

Rec−Re0

)

/Reo，其中Rec為調

制作用下的臨界雷諾數，Re₀為定轉速下的臨界雷諾數。結果發現流 場在低頻時有明顯的負向偏移，中、高頻時有微小的正向偏移。

Kuhlman(1984)，Kuhlman，Roth 和 Lucke(1989)將 Navier-Stokes 方程式轉換為 Lorenz 方程組，並以有限差分法求解，藉以探討外圓 柱靜止，內圓柱調制轉速旋轉的問題。結果發現調制作用會使得基態 流場(basic flow)趨向不穩定，而當一個週期中，在內圓柱的調制轉速 出 現 超 臨 界 的 情 形 下 ， 在 較 大 的 調 制 振 幅 時 ， 則 會 有 次 調 諧 流 (subharmonic)的產生。 Ganske，Gebhartdt 和 Grossmann(1994)以數值方法分析內圓柱調 制轉動，外圓柱靜止條件下的穩定性。結果發現在低調制頻率時會產 生較不穩定的情形，而在高頻時臨界雷諾數則相等或略高於非調制 值。另外，在相同的調制頻率下，發現較大的調制振幅會有較低的臨 界雷諾數，亦即流場較不穩定。

李鎰清(2001)藉由一改良的 Lorenz 模型，輔以 Runge-Kutta 數值 方法對調制旋轉圓柱間流場求解，探討調制作用對流場的影響。發 現在外管固定，內管調制轉動時，低頻流場會產生負向偏移，中、 高頻流場出現正向偏移；在外管調制轉動，內管等速轉動時，低頻 流場產生正向偏移，而隨著調制頻率上升，最終偏移量會趨近於零。

1.2 研究目的

本文嘗試以頻譜法對調制普修流場的流線擾動方程式求其數值 解後，藉流線函數定義得到流場速度。而後尋找不同調制條件下的臨 界雷諾數，繪製成偏移量圖以探討調制作用對流場臨界雷諾數的影 響。並配合時間級數圖、FFT 頻譜圖、相平面圖及 Poincare 映射圖等 分析方式，探討超越臨界雷諾數後的調制流場行為。

第 二 章

數學模式

本文中研究的調制普修流之物理模式如圖一所示。兩平行板其間 充滿黏性流體，相距2H 長度。吾人定義座標在水平方向為x，並假 設x方向為無限長，垂直方向為y，水平方向施加一調制壓力梯度： ) cos 1 ( / x c t P ∂ = +ε ω ∂ − 其中，ε為調制強度，ω為調制頻率。

2.1 統御方程式

為了簡化問題，吾人作如下假設： 1.系統為二維。 2.考慮不可壓縮之牛頓流體。 3.等溫流場，μ =const.。 則統御方程式可表示如下： 連續方程式(continuity equation)： 0 = ∂ ∂ + ∂ ∂ y v x u (2.1) 動量方程式(momentum equations)： ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 y u x u x P y u v x u u t u _μ ρ (2.2a) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 y v x v y P y v v x v u t v _μ ρ (2.2b) 上式中μ為黏滯係數(viscosity)，ρ為流體密度(density)。

邊界條件(boundary conditions) ：

(

x,±H,t

) (

=v x,±H,t

)

=0 u (2.3)

2.2 無因次化

吾人在此定義無因次參數如下： H x x∗ = H y y∗ = 0 u u u∗ = 0 u v v∗ = 2 0 u P P ρ = ∗ H t u t∗ = 0 ν H u₀ Re= 0 u H ω ω∗ = 上式中u₀為ε =0時流道中央處最大速度。將各無因次參數代入統御方 程式並除去上標後，得到無因次化的統御方程式如下： 連續方程式(continuity equation)： 0 = ∂ ∂ + ∂ ∂ y v x u _(2.4) 動量方程式(momentun equations)： ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 Re 1 y u x u x P y u v x u u t u _(2.5a) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 Re 1 y v x v y P y v v x v u t v _(2.5b) 其中邊界條件(boundary conditions)為：

(

x,±1,t

) (

=v x,±1,t

)

=0 u (2.6)

2.3 基態解

考慮流體在基態(basic state)時為一維週期性運動，u =u( ty, )、 0 v= ，則(2.5a)可簡化為：

(

)

2₂ Re 1 cos 1 Re 2 y u t t u ∂ ∂ + + = ∂ ∂ _{ε ω (2.7)} 邊界條件轉變成：

(

± t1,

)

=0 u (2.8) 令基態速度u

( )

y,t =us

( )

y +up

( )

y,t ，其中us為與時間無關之穩態 解，而up為與時間有關之暫態解。將其代入(2.7)中並利用邊界條件可 得： 2 1 y us = − (2.9a) ( ) ( ) ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + + − ⎟ ⎠ ⎞ ⎜ ⎝ ⎛− × = ω β_β+₊β_{β −}−_ββ₊+_ββ ω ε i t y i_iy y_ii y p e e e e e i real u 1 Re 2 _(2.9b) 其中， 2 Reω β = 。

2.4 擾動方程式之建立

在基態上加入擾動量後的流場速度及壓力表示如下：

( )

y t u

(

x y t

)

u u = , + ′ , , (2.10a)

(

x y t

)

v v=0+ ′ , , (2.10b)

( )

x t P

(

x y t

)

P P= , + ′ , , (2.10c) 然後，將(2.10)代入無因次的統御方程式中即可到速度的擾動方程式 如下： 連續方程式： 0 = ∂ ′ ∂ + ∂ ′ ∂ y v x u _(2.11)

動量方程式： = ∂ ′ ∂ ′ + ∂ ∂ ′ + ∂ ′ ∂ ′ + ∂ ′ ∂ + ∂ ′ ∂ y u v y u v x u u x u u t u ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∂ ′ ∂ + ∂ ′ ∂ + ∂ ′ ∂ − 2₂ 2 ₂ Re 1 y u x u x P _(2.12a) = ∂ ′ ∂ ′ + ∂ ′ ∂ ′ + ∂ ′ ∂ + ∂ ′ ∂ y v v x v u x v u t v ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∂ ′ ∂ + ∂ ′ ∂ + ∂ ′ ∂ − 2_{2 2} Re 1 y v x v y P _(2.12b) 其邊界條件為：

(

,±1,

)

= ′

(

,±1,

)

=0 ′ x t v x t u (2.13) 利用流線函數和速度分量之關係，即 y u ∂ ∂ = ′ ψ 、 x v ∂ ∂ − = ′ ψ ，則(2.12) 可表示成： = ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ ∂ 2 2 2 2 2 y x y u x y x y y x u y t ψ ψ ψ ψ ψ ψ ψ ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∂ ∂ + ∂ ∂ ∂ + ∂ ′ ∂ − 3₂ 3₃ Re 1 y y x x P ψ ψ _(2.14a) = ∂ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ − ∂ ∂ − ∂ ∂ ∂ − y x x x y x u x t ψ ψ ψ ψ ψ ψ 2 2 2 2 2 2 ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ∂ ∂ ∂ + ∂ ∂ − ∂ ′ ∂ − 3₃ 3 ₂ Re 1 y x x y P ψ ψ (2.14b) 將(2.14a)對y偏微分減去(2.14b) 對x偏微分，消去壓力項，得流 線函數的非線性擾動方程式如下： ψ ψ ψ ψ ψ ψ ψ 2ψ 2 2 Re1 Δ = Δ ∂ ∂ ∂ ∂ − Δ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ − Δ ∂ ∂ + Δ ∂ ∂ y x x y x y u x u t (2.15) 其邊界條件為：

(

, 1,

)

(

, 1,

)

(

, 1,

)

=0 ∂ ± ∂ = ∂ ± ∂ = ± y t x x t x t x ψ ψ ψ (2.16)

(2.15)中 2₂ 2₂ y x ∂ ∂ + ∂ ∂ = Δ ， 2 4₄ 2₂ 4 ₂ 4₄ y y x x ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ = Δ 。 在(2.15)式中若考慮調制振幅ε =0(無調制作用)的情況，並忽略左 式的第四和第五兩非線性項，即為Orr-Sommerfeld 方程式。換言之， Orr-Sommerfeld 方程式可視為(2.15)式的一個特例。

第 三 章

數值方法

在上一章中推導出的流線函數擾動方程式，因其具有非線性項， 無法直接求得其解析解，需藉由數值方法的使用方能處理。本文中採 用頻譜法處理此問題。

3.1 頻譜法(Spectral method)

吾人首先以雙重級數展開模擬流線函數，同時配合邊界條件決定 各空間方向的展開函數。然後將擾動方程式離散化，再將雙重級數展 開式代入經重新整理劃分的離散方程式。最後則是使用 Galerkin method 處理x方向及Collocation method 處理y方向得到代數方程式。

3.1.1 雙重級數展開

吾人以雙重級數展開來模擬流線函數，因為流場在x方向呈現週 期性，故採用複數形式的 Fourier 級數展開，在y方向則採用符合邊 界條件的基底函數展開，此函數由具有快速收歛特性及高精準度的 Chebyshev 多項式所構成，至於展開係數則為時間的函數。可將流線 函數表示如下：

(

)

( ) ( )

i mx n m n mn t y e a t y x φ α ψ , , =

∑∑

(3.1) 其中基底函數表示如下：

( )

⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ − − + − − − ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − + − = 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 3 1 3 1 1 4 4 φ 而Tn

( )

y 為第一型n階 Chebyshev 多項式，其定義如下：

T_n

( )

y =cos

(

n⋅cos−1 y

)

, n=0,1,2... 其中y的定義域為

[ ]

−1,1 。

3.1.2 離散方程式

針對流線函數的擾動方程式(2.16)加以離散化(discretization)。在 其線性項部份，凡是本身具有基態速度u者採用Adam-Bashforth Method 離散，即 2 3 ₋ −1 = f k f k f ，其餘採用Crank-Nicolon Method 離散， 即 2 1 k k f f f = + + ；非線性項則一律採用Adam-Bashforth Method 進行離 散。經推導出離散化後的擾動方程式如下所示： Δ +1− Δ2 +1 = Re 2 k k dt ψ ψ 1 1 2 2 3 − _Δ − ∂ ∂ + Δ ∂ ∂ − Δ k k k k x u dt x u dt _{ψ ψ} ψ k k k k k dt x y u dt x y u dt ψ ψ 1 2_ψ 2 1 2 2 2 Re 2 2 2 3 _{+ Δ} ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ + − − 1 1 2 2 3 − _Δ − ∂ ∂ ∂ ∂ + Δ ∂ ∂ ∂ ∂ − k k k k x y dt x y dt ψ _ψ ψ _ψ 1 1 2 2 3 − _Δ − ∂ ∂ ∂ ∂ − Δ ∂ ∂ ∂ ∂ + k k k k y x dt y x dt ψ _ψ ψ _ψ (3.2)

3.1.3 代數方程式

將模擬流線函數的雙重級數(3.1)代入(3.2)當中，然後在x方向以

Galerkin method 處理，y方向以定置法(Collocation method)處理，而

定置點的取法為 s N N s y_s , 1,2, , 1 cos = _L + = π 。如此可將統御方程式轉化 為代數方程式，以矩陣形式表示如下：

k+1 = k,k−1 b AX (3.3) 其中，係數矩陣A為一與時間無關之方陣，大小為

(

)

2 N M × ，矩陣X的 元素為待求的時間係數群，其和矩陣b的大小皆為

(

M×N

)

。此矩陣等 式的組成如下： ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ MM A A A L M O M M L L 0 0 0 0 0 0 22 11 = ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ +1 2 1 k M a a a M 1 2 1 2 1 − ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ k M k M LII LII LII LI LI LI M M 1 2 1 2 1 − ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + k M k M H H H G G G M M 其中，係數矩陣A及矩陣LI、LII、G、H內的元素通式表示如下： n

( )

s n

( )

s n

( )

s n

( )

s mm sn m y dt y m dt y y m A =−α2 2φ +φ ″ − α4 4φ + α2 2φ ″ Re Re 2 n( )

( )

ys dt 4 Re 2 φ − (3.4a)

( )

( )

_n

( )

_s n s n n n s n mm s m y dt y y m LI α2 2φ φ α4 4φ Re 2

∑

∑

∑

+ ″ + − =

( )

( )

( )

(

) ( )

n s n k s n s n s n n y y t y u m i dt y dt y m dt_{α φ}

_∑

_φ

_∑

_{α φ}

∑

″ + + ∂ _∂ − 2 2 4 2 ₂, 2 3 Re 2 Re

(

) ( )

(

s k

) ( )

n s n s n k s n y t y u m i dt y t y u m i dt ₋ ″ +

∑

α φ

∑

α , φ 2 3 , 2 3 3 3 _(3.4b)

(

) ( )

(

s k

) ( )

n s n s n k s n mm s i mu y t y dt y t y u m i dt LII =−

∑

α3 3 ₋₁ φ +

∑

α , ₋₁ φ ″ 2 , 2

∑

(

) ( )

∂ ∂ − − n s n k s y y t y u m i dt _{α φ} 2 1 2 _, 2 (3.4c)

( ) ( )

( ) ( )

_⎥ ⎦ ⎤ ⎢ ⎣ ⎡ _′ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ − =

∑∑

∑∑

p q s q k pq m n s n k mn m s i ma t y p a t y dt G α3 φ 2 φ 2 3

( ) ( )

( )

( )

_{( )}

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ +

∑∑

∑∑

p q s q k pq m n s n k mn t y a t y ma i dt 3 2 3 _{α φ φ}

( ) ( )

( ) ( )

_⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ′ +

∑∑

∑∑

p q s q k pq m n s n k mn t y p a t y a i dt _α3 _φ 3 _φ 2 3

( ) ( )

( ) ( )

_⎥ ⎦ ⎤ ⎢ ⎣ ⎡ _″ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ′ −

∑∑

∑∑

p q s q k pq m n s n k mn t y pa t y a i dt _{α φ φ} 2 3 _(3.4d)

( ) ( )

( ) ( )

_⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ′ − =

∑∑

₋

∑∑

₋ p q s q k pq m n s n k mn m s i a t y p a t y dt H α3 ₁ φ 3 ₁ φ 2

( ) ( )

( ) ( )

_⎥ ⎦ ⎤ ⎢ ⎣ ⎡ _″ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ′ +

∑∑

₋

∑∑

₋ p q s q k pq m n s n k mn t y pa t y a i dt _{α φ φ} 1 1 2

( ) ( )

( ) ( )

_⎥ ⎦ ⎤ ⎢ ⎣ ⎡ _′ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ +

∑∑

₋

∑∑

₋ p q s q k pq m n s n k mn t y p a t y ma i dt _{α φ φ} 1 2 1 3 2

( ) ( )

( )

( )

( )

_⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ −

∑∑

₋

∑∑

₋ p q s q k pq m n s n k mn t y a t y ma i dt 3 1 1 2 α φ φ (3.4e) (3.4b)和(3.4c)為線性項，(3.4d)和(3.4e)為非線性項，其中m=1,2,...,M ， N n s, =1,2,..., 。而未知數最新一次的值，可由前兩次的解求得。 將求得之時間係數代入模擬流線函數的雙重級數中，配合流線函 數的定義則可計算出速度分量。

3.2 數值解的分析方法

3.2.1 時間級數與 FFT 頻譜分析

對於動力系統可選擇其內部與時間相關的函數 f(t)或是其動力 變數的時間級數(Time series)來表示其隨時間演變的過程，本文選取 擾動在y方向的速度分量v′作為研究依據。由於任何時間函數 f (t)通

常可以視為許多週期性分量之疊代所組成，為得知時間函數 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 plane)的觀念則是來 自於相空間。例如一質點運動於一維座標下運動，其位置(x)與速度 ) v ( 則是可表示其狀態之相空間，而此空間為一平面；而當質點運動 於三維座標下時，則其需要六個維度的相空間來描述其瞬時狀態，其 中包含三個位置分量及三個速度分量，而相空間的架構有許多種方 式，例如動量可以取代速度分量。對於複雜參數的動力系統，其相平 面通常極為紊亂，此時須藉由其他方法才較為容易判讀系統狀態。至

於本文則是採用水平向及鉛直向的擾動速度作為表示流場狀態之相 空間。

3.2.3 Poincare 映射圖

Poincare 映射圖是一種用來簡化相平面圖的方法，藉由固定週期 下，對相平面的軌跡取一投影來取樣觀察(strobe)，而在此吾人是以壓 力梯度的調制週期來做取樣，使得吾人可以清楚的觀察動力系統的演 變過程。 由 Poincare 映射圖吾人可以得知動力系統頻率(ω₀)對於取樣頻 率(ω_s)的比例關係。當動力系統的自然頻率為(p q)ωs時(p q為有理 數)，則於 Poincare 映射圖上會顯現q個映射點於一個封閉軌道之上， 而剩下[q−(p+1)]個點則會與先前的映射點位置重複，是故可以忽略 這些映射點。對一個消散(Dissipative)系統而言，其映射點會移動至適 當的吸子(Attractor)而形成一點，而當取樣頻率與自然頻率成有理數 關係時，則會形成散佈在封閉軌道上的點群，若完全不相關時，其映 射圖上的點甚至會填滿而形成一封閉曲線。

第 四 章

結果與討論

在本章中先探討在非調制(ε=0)的條件下流場的穩定性，其次討論 調制作用對流場臨界雷諾數的影響，最後則藉時間級數圖、頻譜圖、 相平面圖及Poincare 映射圖探討超越臨界雷諾數後之調制流場行為。

4.1 非調制(ε=0)流場的穩定性

Fortin 與 Jardak(1994)以有限元素法及洪英棋(2005)以頻譜法對於 非調制的普修流場在波數(wave number)α=1.020545 時求解，皆獲得臨 界雷諾數 Re0=5772.22 的精確結果。在前一章當中流線函數的 x 與 y 方向分別以複數形式Fouries series 及 Chebyshev polynomials 作雙重級 數展開，其展開項數(M×N)將影響數值計算結果的準確度，因此吾人 必須藉改變(M×N)的項數去尋找準確的臨界雷諾數。理論上，展開項 數取的愈多其結果之準確度也愈高，但項數太多將會耗費大量的運算 時間，因此取最少的展開項數而能達到相當的準確度才是最佳的選 擇。 先考慮展開項數(3×5)項，發現臨界雷諾數 Re0=573，然後逐漸提 高 N 項，所得到的臨界雷諾數隨之改變如圖二所示，一直到(3×43) 項，此時可發現 Re0=5664，再增加至(3×45)項仍是 Re0=5664，顯示 增加N 項已無法改善準確度，須朝增加M 項著手。 增 加 M 項 使 展 開 項 數 由 (3×43) 項 增 加 至 (7×43) 項 得 到 Re0=5772.22，再增加至(9×43)仍是 Re0=5772.22，如圖三所示，因此 (7×43)項為最佳的展開項數。

4.2 調制作用對臨界雷諾數的影響

為了解在不同的調制條件下，流場臨界雷諾數的變化情形，在此 定義偏移量(Threshold shift)Δ =(Rec−Re0)/Re0，其中Rec為調制作 用下的臨界雷諾數，正偏移量代表調制臨界雷諾數大於非調制時的臨 界雷諾數(即 5772.22)，負偏移量代表調制臨界雷諾數小於非調制時的 臨界雷諾數。 在不同調制振幅及頻率時的臨界雷諾數列於表 1，若以偏移量表 示如圖四。由圖四可知流場在頻率ω=0.01 至頻率 ω=1 之間具有負的 偏移量，且頻率越低偏移量越小，亦即臨界雷諾數越小，流場越提前 發生不穩定。而在頻率ω=1 至頻率 ω=10 之間則具有正的偏移量，亦 即此時流場相較非調制的結果有較高的臨界雷諾數，即流場延後發生 不穩定。至於調制振幅的作用則是放大提前或延後發生不穩定，亦即 原本不穩定的流場變得更不穩定，原本穩定的流場則更加穩定。

4.3 臨界雷諾數後之調制流場行為

4.3.1 調制頻率的影響

圖五為當調制振幅 ε=1，調制頻率 ω=0.1，雷諾數 Re=6000 時的 (a)時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。為更 容 易觀 察到 調制頻率 與頻 譜頻率 的 關 係 ， 在 此 以 角 頻 率(angular frequency)的形式來表示頻譜頻率，即角頻率形式等於頻率(frequency) 形式乘上2π。由圖五(a)發現流場呈週期性擾動，經頻譜分析由圖五 (b)得知其頻譜頻率為單一頻率 0.1(即 ω)，表示流場為單一週期運動。 觀察圖五(c)呈現一封閉曲線，顯示於圖五(d)則為一點。圖六為當 ε=1，ω=0.1，Re=6500 時的(a)時間級數圖、(b)頻譜圖、(c)相平面圖 及(d) Poincare 映射圖。由圖六(a)發現流場仍呈週期性擾動，但振幅

隨著時間略有變化，由圖六(b)得知除原來的頻率又另外產生兩個強度 微弱的頻率分別是0.075(即 3ω/4)和 0.125(即 5ω/4)。圖六(c)中相平面 為一由許多封閉曲線所疊合而成的圓環，顯示於圖六(d)則為一線段。 圖七為當 ε=1，ω=0.1，Re=8500 時的(a)時間級數圖、(b)頻譜圖、(c) 相平面圖及(d) Poincare 映射圖。由圖七(a)發現流場雖仍呈週期性擾 動，但波形與先前相較已有明顯的不同，經頻譜分析由圖七(b)發現又 增加一新的頻率0.15(即 3ω/2)。圖七(c)相平面呈現由許多封閉曲線所 疊合而成的兩個圓環，投影於圖七(d) 則為兩條線段。圖八為當 ε=1， ω=0.1，Re=10000 時的(a)時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。由圖八(a)發現流場呈現混亂狀態，再經由圖八(b) 觀察亦發現流場不具週期性，此時流場已進入混沌(chaos)狀態。圖八 (c)中相平面由無數雜亂的封閉曲線所組成，顯示於圖八(d)則為無數 的映射點。 圖九為當調制振幅ε=1，調制頻率 ω=1，雷諾數 Re=6000 時的(a) 時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。圖九(a) 顯示擾動呈單一週期變化，由圖九(b)觀察頻譜頻率為單一頻率 1(即 ω)，代表單一週期(periodic)運動。圖九(c)中的相平面為一封閉曲線， 降維後顯示於圖九(d)為一映射點。圖十為當調制振幅 ε=1，調制頻率 ω=1，雷諾數 Re=6500 時的(a)時間級數圖、(b)頻譜圖、(c)相平面圖 及(d) Poincare 映射圖。由圖十(a)發現流場雖不是單一週期運動，但 仍可觀察出具有週期性，經由頻譜分析由圖十(b)得知頻譜頻率有 1(即 ω)和 1.5(即 3ω/2)兩個，表示流場為擬週期(quasi-periodic)運動。圖十 (c)中相平面為兩條規則性的封閉曲線，顯示於圖十(d)為兩映射點。 圖十一為當調制振幅 ε=1，調制頻率 ω=1，雷諾數 Re=8500 時的(a) 時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。由圖十 一(a)發現其與圖十(a)相似，唯一差異是振幅較大，經頻譜分析由圖 十一(b)觀察同樣是得到 1 和 1.5 兩個頻率，所以此時流場仍處於擬週

期運動。圖十一(c)中的相平面同樣為兩條規則性的封閉曲線，顯示於 圖十一(d)同樣為兩映射點。圖十二為當調制振幅 ε=1，調制頻率 ω=1， 雷諾數 Re=10000 時的(a)時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。由圖十二(a)及圖十二(b)皆無法觀察任何週期性，代 表流場已進入混沌狀態。圖十二(c)中呈現由無數曲線填滿的複雜圖 形，顯示於圖十二(d)中為無數的映射點。 圖十三為當調制振幅 ε=1，調制頻率 ω=10，雷諾數 Re=6000 時 的(a)時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。由 圖十三(a)發現流場為單一週期運動，經頻譜分析由圖十三(b)得知此 單一頻率為6.8。圖十三(c)中相平面為一封閉曲線，顯示於圖十三(d) 中為一映射點。圖十四為當調制振幅 ε=1，調制頻率 ω=10，雷諾數 Re=6500 時的(a)時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。由圖十四(a)發現流場呈週期性擾動，經頻譜分析由圖十四(b) 得知除 6.8 之外又產生另一個主要頻率 10.7，此時流場為擬週期運 動。圖十四(c)中相平面為兩封閉曲線，顯示於圖十四(d)中為兩映射 點。圖十五為當調制振幅 ε=1，調制頻率 ω=10，雷諾數 Re=8500 時 的(a)時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。由 圖十五(a)發現流場已呈混亂的狀態，進一步作頻譜分析由圖十五(b) 觀察到在緊鄰兩個主要頻率的位置生成一些不相干的頻率，代表流場 此時產生輕微的混沌現象。圖十五(c)中相平面已呈現由無數曲線填滿 的圖形，顯示於圖十五(d)中為無數的映射點。圖十六為當調制振幅 ε=1，調制頻率 ω=10，雷諾數 Re=10000 時的(a)時間級數圖、(b)頻譜 圖、(c)相平面圖及(d) Poincare 映射圖。由圖十六(a)發現流場擾動變 得更加混亂，經頻譜分析由圖十六(b)得知流場已不具週期性，表示流 場已完全進入混沌狀態。圖十六(c)中相平面依然是由無數曲線填滿的 圖形，且比先前的圖形顯得更加濃密，顯示於圖十六(d)中也依然是無 數的映射點。

經由上述的討論可得到在超越臨界雷諾數後調制頻率對流場的 影響，可歸納如下: 在低頻(ω=0.1)及中頻(ω=1)時，流場在進入混沌 之前其週期性與調制頻率存在若干比例關係，於低頻時存在有 ω、 3ω/4、5ω/4 及 3ω/2 的關係；於中頻時存在有 ω 和 3ω/2 的關係。但在 高頻(ω=10)時，此關係則不復存在，此現象顯示在高調制頻率時，壓 力梯度變動得太快導致流場中流體來不及充分吸收和傳遞能量所 致。另外，隨著調制頻率的增加，相對造成流場變動得愈快，導致流 場愈容易進入混沌的狀態。選取中頻(圖九~十二)與高頻(圖十三~十六) 作比較則可觀察到此現象，當雷諾數為 8500 時，高頻流場已開始進 入混沌狀態，而中頻流場則仍處於擬週期運動。

4.3.2 調制振幅的影響

圖十七為當調制振幅 ε=0.5，調制頻率 ω=1，雷諾數 Re=6500 時 的(a)時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。由 圖十七(a)發現擾動呈單一週期變化，經頻譜分析由圖十七(b)得知產 生單一頻率 1(即 ω)，代表單一週期運動。圖十七(c)的相平面為一封 閉曲線，顯示於圖十七(d)為一映射點。圖十八為當調制振幅 ε=0.5， 調制頻率 ω=1，雷諾數 Re=8500 時的(a)時間級數圖、(b)頻譜圖、(c) 相平面圖及(d) Poincare 映射圖。圖十八(a)顯示流場呈週期性擾動， 經頻譜分析由圖十八(b)得知產生兩個頻率分別是 1(即 ω)和 1.5(即 3ω/2)，此時流場為擬週期運動。圖十八(c)中相平面為兩封閉曲線， 顯示於圖十八(d)為兩映射點。圖十九為當調制振幅 ε=0.5，調制頻率 ω=1，雷諾數 Re=10000 時的(a)時間級數圖、(b)頻譜圖、(c)相平面圖 及(d) Poincare 映射圖。觀察圖十九(a)發現流場同樣呈週期性擾動， 經頻譜分析由圖十九(b)得知同樣產生 1 和 1.5 兩個頻率，流場仍為擬 週期運動。圖十九(c)中相平面同樣為兩封閉曲線，顯示於圖十九(d) 同樣為兩映射點。圖二十為當調制振幅 ε=0.5，調制頻率 ω=1，雷諾

數 Re=11000 時的(a)時間級數圖、(b)頻譜圖、(c)相平面圖及(d) Poincare 映射圖。由圖二十(a)發現流場擾動呈現混亂，由圖二十(b) 得知頻譜分析的結果，流場進入混沌狀態。圖二十(c)中相平面為無數 的封閉曲線，投影於圖二十(d)為無數的映射點。 選取在中頻時較大調制振幅 ε=1 (圖九~十二)和較小調制振幅 ε=0.5 (圖十七~二十)作比較可觀察到在超越臨界雷諾數後調制振幅 對流場的影響。在雷諾數為 6500 時，較大調制振幅已進入擬週期運 動，較小調制振幅仍處於單一週期運動，當雷諾數增加至10000 時， 較大調制振幅已進入混沌狀態，較小調制振幅仍處於擬週期運動。由 此得知增加調制振幅會使得流場容易分配能量至各種與調制頻率無 關的頻率分量，導致流場容易進入混沌狀態。

第 五 章

結 論

本文中藉由頻譜法以 Chebyshev 多項式及複數形式的 Fourier 級 數對流線函數作雙重級數展開來模擬調制平面普修流並求其數值 解，研究調制作用對流場產生的影響。綜合前一章所得到的結果，吾 人得到以下的結論。 1. 在非調制(ε=0)條件下，當展開項數為(7×43)項，dt=0.02 時，得到 α=1.020545，Re0=5772.22 的精確值，與 Fortin 等人（1994）及洪 英棋(2005)的結果一致。 2. 就調制作用對流場臨界雷諾數的影響而言，流場在低頻時，有明 顯的負偏移量，而在中頻至高頻之間，則有些微的正偏移量。至 於調制振幅則具有放大偏移的效應。 3. 調制流場在臨界雷諾數之後的行為，低頻(ω=0.1)流場原先為單一 週期運動(即 ω)，隨著雷諾數的增加，流場先進入由 ω、3ω/4 及 5ω/4 等頻率疊合而成的週期運動，然後進入由 ω、3ω/4、5ω/4 及 3ω/2 疊合而成的週期運動，最終進入混沌狀態。至於中頻(ω=1)及高頻 (ω=10)流場原先同樣為單一週期運動，隨著雷諾數的增加，流場進 入擬週期運動，最終處於混沌的狀態。 4. 就調制作用對超越臨界雷諾數後之流場的影響而言，調制頻率的 增加，將會使得流場的週期減短，導致流場的變化過程加快;而增 加調制強度，不但會使流場的振盪變大，也容易使流場產生與調 制頻率無關的振盪，因此增加調制頻率與調制強度都將會使得流 場容易形成混沌。

參 考 文 獻

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, 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.

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

Donnelly, R. J. 1964 Experiments on the stability of viscous flow between rotating cylinders III. Enhancement of stability by modulation. Proc. Roy. Son. London A 281, 130-139.

Fortin A. and Jardak M., 1994 Old and new results on the two-Dimensional Poiseuille flow. J. Computational Physics. 115, 455. Ganske, A., Gebhardt, T., and Grossmann, S. 1994 Taylor-Couette flow

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.

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.

Kuhlmann, H., Roth, D. and Lucke, M. 1989 Taylor vortex flow under harmonic modulation of the driving force. Phys. Rev. A 39, 745-762. Lin, C. C., 1955 The Theory of Hydrodynamic Stability. Cambridge

Monographs On Mechanics And Applied Mathematics.

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

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

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. Journal of Applied Mechanics. 45, 13-18.

李鎰清 2001 圓柱之流場-一種改良 Lorenz 模型. 國立交通大學碩 士論文.

ω _{0.01 0.02 0.05 0.1 0.2 0.5 1 2 5 10 20 50 100} c Re _{5380 5382 5382 5386 5456 5620 5768 5864 5844 5776 5772 5772 5772} (a) ω _{0.01 0.02 0.05 0.1 0.2 0.5 1 2 5 10 20 50 100} c Re _{5642 5644 5644 5646 5674 5728 5762 5812 5806 5774 5772 5772 5772} (b) 表 1 調制振幅(a)ε=1 和(b)ε=0.5 時，在不同調制頻率下所對應的臨界雷諾數。

圖一 調制普修流物理模式示意圖 ) cos 1 ( t c x P _{ε ω} + = ∂ ∂ − _2H x y

圖二 當M =3，α =1.020545時，針對非調制流場(ε =0)不同N 所得到

圖三 當N =43，α =1.020545時，針對非調制流場(ε =0)不同M 所得到

(a)

(b)

圖五 當ε =1，ω =0.1，Re=6000時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖六 當ε =1，ω =0.1，Re=6500時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖七 當ε =1，ω =0.1，Re=8500時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖八 當ε =1，ω =0.1，Re=10000時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖九 當ε =1，ω =1，Re=6000時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十 當ε =1，ω =1，Re=6500時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十一 當ε =1，ω =1，Re=8500時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十二 當ε =1，ω =1，Re=10000時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十三 當ε =1，ω =10，Re=6000時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十四 當ε =1，ω =10，Re=6500時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十五 當ε =1，ω =10，Re=8500時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十六 當ε =1，ω =10，Re=10000時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十七 當ε =0.5，ω =1，Re=6500時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十八 當ε =0.5，ω =1，Re=8500時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖十九 當ε =0.5，ω =1，Re=10000時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續

(a)

(b)

圖二十 當ε =0.5，ω =1，Re=11000時的(a)時間級數圖、(b)頻譜圖、

(c)

(d) 續