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

國立交通大學

機械工程學系

碩 士 論 文

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

Linear Instability of Plane Poiseuille Flow with Modulated

Pressure Gradient

研 究 生：邱文山

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

Linear Instability of Plane Poiseuille Flow with Modulated Pressure Gradient 研 究 生：邱 文 山 Student：Wen-Sung Chiou 指導教授：楊 文 美 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 2007

Hsinchu, Taiwan, Republic in China

摘 要

為了解脈動對流動行為的影響，本文以線性方法探討調制平面普 修流的不穩定現象，首先建立流線擾動方程式，依據 Normal Mode Analysis 及 Floquet 定理（Coddington and Levinsin[1955]）可將擾動

量以 Fourier - Chebyshev 雙重級數展開來模擬流線函數，並利用 Galerkin method 與定置法（Collocation method） 將擾動方程式轉化

成代數形式之特徵值問題求數值解，探討不同調制振幅與頻率下，特 徵值對流場穩定性的影響，並尋找調制振幅、頻率與臨界雷諾數之間 的關係。 在不同調制參數求解過程中，依特徵值虛部為零特性可找出相對 應的臨界曲線與臨界雷諾數，最後將結果整理為臨界雷諾數與頻率關 係圖，圖中顯示當調制振幅愈大、調制頻率愈低時會導致流場愈不穩 定，在高頻率時則會趨於無調制，綜合而言調制作用對平面普修流具 有使流場提前發生不穩定的功用。

誌 謝

首先感謝指導教授楊文美老師兩年來的悉心指導與諄諄教誨，使 本論文能順利的完成。恩師的嚴謹教學態度，讓我們培養出獨立思考 的能力與認真的處事態度，使學生受益良多，浩瀚師恩，非筆墨難以 形容，永銘於心，謹誌於此，以表由衷感謝與敬意。 口試期間承蒙傅武雄老師、游明輝教授及胡毓仁博士等不吝賜教 ，並提供諸多寶貴意見，使本文更加完備，深表謝意。 感謝實驗室學長豪傑，同學武伸以及學弟梓源、榮豐及國元在生 活上的熱情幫忙與大力支持，使得兩年的研究生活更為充實平順，並 且感謝在交通大學求學期間所有教導過我的師長，使我在學識的領域 上有所成長與突破，在此謹致上我最誠摯的謝意。 最後謹以本論文獻給我的父母及親友，感謝你們對我的生活無止 境的關懷與照顧，並在精神上全力的支持與鼓勵，此恩此情銘記於 心。

目 錄

摘要………....….….i 誌謝………...…..ii 目錄………..…..iii 表目錄………....….v 圖目錄………...……….…vi 符號說明………..viii 第一章 緒論………...1 1-1 前言………....1 1-2 文獻回顧……….………...2 1-3 研究目的……..………...5 第二章 數學模式………...….….6 2-1 統御方程式………...…...6 2-2 無因次化………...7 2-3 基態解………...8 2-4 擾動方程式之建立………...9 第三章 數值方法………...11 3-1 Fourier-Chebyshev 雙重級數展開……….11 3-2 建立代數特徵方程式………..12 3-3 求解特徵值………..15 第四章 結果與討論………...16

4-1 調制作用與基態(Basic state)關係………..16

4-2 展開項數的影響………..17

4-2-1 無調制(

ε

=0)狀態下展開項數 ... 17

4-2-2 調制狀態下展開項數 ... 17

4-2-3 調制狀態下極值區域展開項數結果 ... 18

4-3 臨界曲線(Marginal stability curve)……….18

4-3-1 無調制狀態下之臨界曲線 ... 18 4-3-2 調制狀態下之臨界曲線 ... 19 4-3-3 調制參數對臨界曲線之影響 ... 19 4-4 流線函數擾動量隨時間變化關係………..20 4-5 調制頻率、振幅範圍與實際心跳參數………..20 4-6 中頻區域………..21 第五章 結論………...23 參考文獻………...25

表 目 錄

表一 無調制結果與 Orszag(1971) 項數結果比較 ... 28 表二 當 ε=1 時調制頻率對 M 項數關係 ... 29 表三 不同 ω、ε 所對應之臨界雷諾數與項數(M)關係 ... 30 表四 不同 ω、ε 所對應之臨界波數(α)值 ... 31 表五 各種動物實際心跳頻率整理(每分鐘計) ... 32

圖 目 錄

圖一 調制平面普修流物理模型示意圖 ... 33 圖二 相同調制頻率下(ω =0.0001)不同調制振幅對基態之影響 (a) ε=1 (b) ε=0.5 (c) ε=0.25 (d) ε=0 ... 34 圖三 相同調制振幅下(ε=1)不同調制頻率對基態之影響 (a) 0.0001 = ω (b)ω=0.001(c)ω=0.01(d)ω =0.1 ... 36 圖四 無調制作用時 M=1，不同 N 項數與臨界雷諾數關係測試 圖 ... 38 圖五 當ω=0.00001、ε=1 時，N=22 不同取項 M 與誤差 ... 39 圖六 無調制平面普修流(ε=0)α與 Re 關係圖 ... 40 圖七 當ε=1 時不同調制頻率下α與 Re 關係圖 (a)ω =0.00001(b) 0.0001 = ω (c)ω=0.001(d)ω =0.0036(e)ω=0.01(f)ω =0.1 ... 41 圖八 當 ε=1 時不同ω範圍對臨界曲線影響圖 (a)0.0003≤ω ≤ 001 . 0 (b) 0.002 ≤ω ≤0.005 _(c) 0.01≤ω ≤0.03 _{... 44} 圖九 相同調制頻率下(ω=0.01)不同調制振幅對臨界曲線之影 響 ... 46 圖十 臨界雷諾數與調制頻率圖 ... 47 圖十一 α與ω關係圖 ... 48 圖十二 當ω =0.0001、ε=1、Re_c = 4005，討論流場中(x,y)=(0.5,0.5) 處流線函數擾動量隨時間變化關係圖(a) Re =3980 (b)

圖十三 當ω =0.00316、ε=1、Re_c =4003.7 ，討論流場中(x,y) = (0.5,0.5) 處 流 線 函 數 擾 動 量 隨 時 間 變 化 關 係 圖 (a) 0 5 39 Re = (b)Re = 4003.7(c)Re = 4050 _{... 50} 圖十四 當ω=0.01、ε=1、Re_c =3055.9，討論流場中(x,y)=(0.5,0.5) 處 流 線 函 數 擾 動 量 隨 時 間 變 化 關 係 圖(a) Re =2900 (b) 9 . 3055 Re = (c)Re =3150 _{... 51} 圖十五 本文結果與 Straatman(2002)結果對照 ... 52

符 號 說 明

符 號

A、B、X 矩陣 c 特徵值 C 常數 f 頻率 H 平行板之間距 N , M 流線函數展開項之項數 M1 展開範圍 P 壓力 Re 雷諾數 c Re 臨界雷諾數 0 Re 非調制臨界雷諾數 t 時間 n T n階Chebyshev多項式 v u, 速度分量 0 u 非調制時流道中央處最大速度 x,y,z 平面座標 α 波數(wave number) ∆ 2 2 2 2 y x ∂ ∂ + ∂ ∂ 4 4 4 2∂ ₊ ∂ + ∂

ε 調制振幅 µ 黏滯係數 ρ 密度 φ y方向展開函數 ψ 流線函數擾動量 ω 調制頻率

上標說明

基態 ´ 擾動量

第 一 章

緒 論

1-1 前言

自然界中存在許多調制液動系統，例如血液循環系統中的流動、 大氣層在白天的加熱與夜晚的冷卻、以及一些生物化學系統 . . . 等。 近年來，調制形式的流場研究逐漸受到重視，其原因為調制作用具有 使流場提前或延後發生不穩定的特性，所以當我們了解調制作用對流 場的影響後，則可有效將其特性加以運用，例如:醫學上的配方混合， 可藉由調制作用使配方混合之流場提前發生不穩定，運用此特性讓流 場中的熱量與質量快速傳遞，可大幅減少混合時間，在工業上也可應 用在印表機的噴霧流，利用調制特性產生更穩定的流場使印表機噴出 的噴墨更為均勻。但調制性流場問題較一般非調制性流場複雜許多， 所以過去學者們會對物理或數學上提出適當假設藉以簡化問題，有許 多基本或代表性的問題被相繼提出，例如：普修流（Poiseille Flow）、 庫耶流（Couette Flow）. . .等。本文中將以平面普修流為模型之基本 流場，利用線性理論探討調制作用對平面普修流流場之影響，在下一 節將介紹前人所作之相關研究可幫助吾人進一步了解平面普修流研 究的演變。

1-2 文獻回顧

普修流是以固定壓力梯度驅使流體流動的流場，在層流狀態下可

經由 Navier-Stokes 方程式得到速度的單一函數基態解。普修流可依

照流體所流經空間的幾何特性大致劃分為兩種類型，一、流經等截面 的圓管普修流(Circular Poiseuille flow)。二、流經兩固定平行板間的 平面普修流(Plane Poiseuille flow)。對平面普修流而言過去前人將其 分為調制與非調制型式的問題加以討論，可整理如下：

1.非調制平面普修流－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(λ)皆

果不但證明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 一致的結果。 2.調制平面普修流

在無調制方面其結果與 Tomas(1953)、Lin(1954)一致，調制方面則提 出以下現象：a.在很小調制振幅下流體會趨於穩定，大調制振幅則會 使流體不穩定；b.不同調制振幅、頻率會導致波數與雷諾數改變，使 臨 界 曲 線(marginal stability curve) 產 生 移 動 ； c. 相 較 於 無 調 制 (Re_c ≈5772.22) Rec愈大則會導致臨界曲線內部涵蓋區域越狹窄。 Herbert(1972)利用很小調制振幅描述平面普修流模型在邊界附 近能量轉換過程，經能量法分析後，發現擾動剪應力波除以調制剪應 力波大於1/2 時，調制會抑制能量轉換為擾動形式，使擾動量呈週期 性遞減，流場呈現穩定狀態。 Hall(1975)以數學形式討論高頻調制下之平面普修流模型，探討 波 數 、 臨 界 雷 諾 數 關 係 並 推 導 出 以 下 數 學 型 式 ：α=1.02h−1、 )] ( 10 709 . 2 1 [ 4 . 5774 Re 2 7 3 5 2 − ₊ − ∈ − = σ O σ c 、 0.2639[1 4.8 10 ( )] 3 6 2 5 2 − ₊ − ∈ − = σ O σ c 其中σ ω 2/ν h = : 頻率參數、∈ : 振幅參數(∈= ，無調制)、 h : 1/20 模型寬度、ν:動黏滯係數(kinematic viscosity)，當σ →∞，α =1.02、 4 . 5774 Re₀ = 、c₀ =0.2639結果趨於無調制狀態。 Davis(1976)針對各種調制形式問題做整理，其中在平面普修流方 面探討Grosch & Salwen、Herbert . . .等人所做之研究，提出 a.當極低 頻時，ReM ~ Re₀/(1+∆)_{；b.當極高頻時， Re (Re Re )/Re 0} 0 0 ≈ − = M M δ (i.e. ReM =Re0)。

Straatman(2002)沿用 Grosch & Salwen 所提出之調制平面普修

流模型，將調制頻率無因次化定義引用Wo=ω / Re 數(Womersley

number)作修改，使用大振幅求解臨界雷諾數、波數、與調制參數間 之關係，結果顯示較高Wo No. 與 Davis(1976)、Hall(1975)一致趨於

無調制，但在較低Wo No. 處作者提出可能具有奇異點存在，結果則 與 Davis(1967)的低頻現象有所不同。 廖啟豪（2006）採用頻譜法對調制普修流場的流線擾動方程式求 其數值解後，藉流線函數定義得到流場速度。而後尋找不同調制條件 下的臨界雷諾數，繪製成偏移量圖以探討調制作用對流場臨界雷諾數 的影響。並配合時間級數圖、FFT 頻譜圖、相平面圖及 Poincare 映射 圖等分析方式，探討超越臨界雷諾數後的調制流場行為。 3.心臟實驗 Nerem(1972)等人使用探針量測狗的胸腔大動脈，並記錄瞬時與 時間-平均速度分佈，實驗結果顯示擾動發生與雷諾數峰值及頻率參 數有關。

1-3 研究目的

吾人欲了解脈動現象對流場之影響，本文嚐試使用平面普修流模 型，透過壓力調制作用模擬脈動效應，藉以探討脈動行為對流場穩定 性之影響。在以下章節中將推導出流線函數線性擾動方程式，並將其 轉化為代數形式之特徵值問題，最後利用數值方式求解，探討調制參 數與臨界雷諾數之關係。

第 二 章

數 學 模 式

考慮上下相距2H之固定平行板間其內部充滿黏性流體如圖一， 在y 軸方向，原點座標固定於 1/2 流場寬度處；在 x 軸方向假設無限 長並施加調制壓力梯度： ) C 1 ( / x C os t P ∂ = +ε ω ∂ − 其中，ε為調制振幅、ω為調制頻率、 C 為常數。

2-1 統御方程式

經由下列假設條件簡化問題： 1.系統為二維、無 z 方向分量（w=0） 2.不可壓縮流體（ ρ =const.） 3.假設流場為等溫，不受溫度影響（T=const.） 則統御方程式可表示如下： 連續方程式(continuity equation)： 0 = ∂ ∂ + ∂ ∂ y v x u (2.1-1a) 動量方程式(momentum equations)：

⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 y u x u x P y u v x u u t u _µ ρ (2.1-1b) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 y v x v y P y v v x v u t v _µ ρ (2.1-1c) 上式中 µ 為黏滯係數(viscosity)， ρ 為流體密度(density)。 邊界條件(boundary conditions) ：

(

x,±H,t

) (

=v x,±H,t

)

=0 u (2.1-2)

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.2-1a) 動量方程式(momentun equations)：

⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 Re 1 y u x u x P y u v x u u t u (2.2-1b) ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ + ∂ ∂ − = ∂ ∂ + ∂ ∂ + ∂ ∂ 2 2 2 2 Re 1 y v x v y P y v v x v u t v (2.2-1c) 其中邊界條件(boundary conditions)為：

(

x,±1,t

) (

=v x,±1,t

)

=0 u (2.2-2)

2-3 基態解

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

(

)

2 2 Re 1 C 1 Re 2 y u t os t u ∂ ∂ + + = ∂ ∂ _{ε ω} (2.3-1) 邊界條件為：

(

± t1,

)

=0 u (2.3-2) 令基態速度u

( )

y,t =us

( )

y +up

( )

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

2-4 擾動方程式之建立

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

( )

y t u

(

x y t

)

u u = , + ′ , , (2.4-1a)

(

x y t

)

v v=0+ ′ , , (2.4-1b)

( )

x t P

(

x y t

)

P P = , + ′ , , (2.4-1c) 將(2.4-1)代入無因次的統御方程式（2.2-1）中，可到速度的擾動方程 式如下： 連續方程式： 0 = ∂ ′ ∂ + ∂ ′ ∂ y v x u (2.4-2) 動量方程式： = ∂ ′ ∂ ′ + ∂ ∂ ′ + ∂ ′ ∂ ′ + ∂ ′ ∂ + ∂ ′ ∂ y u v y u v x u u x u u t u ) ( Re 1 2 2 2 2 y u x u x P ∂ ′ ∂ + ∂ ′ ∂ + ∂ ′ ∂ − (2.4-3a) = ∂ ′ ∂ ′ + ∂ ′ ∂ ′ + ∂ ′ ∂ + ∂ ′ ∂ y v v x v u x v u t v _{( )} Re 1 2 2 2 y v x v y P ∂ ′ ∂ + ∂ ′ ∂ + ∂ ′ ∂ − (2.4-3b) 其邊界條件為：

(

,±1,

)

= ′

(

,±1,

)

=0 ′ x t v x t u (2.4-4) 消除非線性項並引用流線函數和速度分量之關係，即 y u ∂ ∂ = ′ ψ 、

= ∂ ∂ ∂ ∂ − ∂ ∂ ∂ + ∂ ∂ ∂ y u x y x u y t ψ ψ ψ 2 2 ) ( Re 1 3 3 2 3 y y x x P ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ − ψ ψ _(2.4-5a) = ∂ ∂ − ∂ ∂ ∂ − 2 2₂ x u x t ψ ψ _{( )} Re 1 2 3 3 3 y x x y P ∂ ∂ ∂ + ∂ ∂ − ∂ ∂ − ψ ψ _(2.4-5b) 為消去壓力擾動P′將(2.4-5a)對y變數偏微分減去(2.4-5b) 對 x 變數偏微分得流線函數的線性擾動方程式如下： ) ( Re 1 ) ( ) ( 2 2 2 ψ ψ ψ ψ = ∆ ∂ ∂ ∂ ∂ − ∆ ∂ ∂ + ∆ ∂ ∂ x y u x u t (2.4-6) 其中∆及_∆2_{表空間運算子，} 2 2 2 2 y x ∂ ∂ + ∂ ∂ = ∆ 、 ₄ 4 2 2 4 4 4 2 2 y y x x ∂ ∂ + ∂ ∂ ∂ + ∂ ∂ = ∆ 其邊界條件為：

(

, 1,

)

(

, 1,

)

(

, 1,

)

=0 ∂ ± ∂ − = ∂ ± ∂ = ± x t x y t x t x ψ ψ ψ (2.4-7) 吾人使用線性理論求得式(2.4-6) 流線函數的線性擾動方程式， 在下一章節將介紹數值求解方法並計算不同調制参數下對臨界雷諾 數之影響。

第 三 章

數 值 方 法

3-1 Fourier-Chebyshev 雙重級數展開

由於擾動方程式為線性偏微分方程式，且方程式中有時間週期函

數的係數，依據Normal Modes Analysis，可將擾動量表示為：

) ( ) , ( ) , , ( i x ct e t y g t y x = α −

ψ

(3.1) 其中α 為 x 方向之波數、c=cr +ici為波速、g( ty, )為時間週期函數， 依據Floquet 定理(Coddington and Levinsin [1955])未知函數g( ty, )之 週期與基態流相同，因此吾人可以使用時間及空間雙重級數加以展開 ，其中時間部份以複數富氏級數（Complex Fourier Series）展開；空 間部分則以第一類n 階 Chebyshev polynomials（Fox and Parker [1968]） 展開，其定義為： ) ( ) ( 1 y Cos n Cos y T_n = ⋅ − (3.1-1) 經由以上過程則可將g( ty, )表示成以下形式：

∑ ∑

− = + = − = M1 1 M 3 N 4 ) ( ) , ( m n t im n mn y e a t y g φ ω (3.1-2) 其中φ_n( y)為基底函數型式如下：

⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ − − + − − − − + − = n odd T n T n T n even T n T n y n n n 1 2 2 3 2 2 0 2 2 2 1 3 3 1 3 1 ) 1 4 ( 4 T ) ( φ _{、n=4,5,6... (3.1-3)} 上式已滿足式（2.4-7）所有邊界條件。在式(3.1-2)中M、 N 為展開項 數，M1=(M-1)/2、m=-M1~M1、n=4~ N+3，展開項數規則為 1,3,5... M= 、N=1,2,3...例如：M = 、3 N= 則3 m=−1,0,1、n=4,5,6。 將式（3.1）帶入式（2.4-6）可整理成下式 ] ) )( Re[( ) ( ₂ 2 2 2 2 2 2 g y u g D c u i g D ∂ ∂ − − − = −α α α D g t( ) Re 2 _−α2 ∂ ∂ + (3.1-4) 其中 y D ∂ ∂ = ， 若在無調制情況下(ε=0），方程式中週期函數的係數、 速度分佈均與時間無關（g =g(y)、u =us），則式（3.1-4）變成下式： ] ) )( Re[( ) ( ₂ 2 2 2 2 2 2 g y u g D c u i g D ∂ ∂ − − − = −α α α (3.1-5) 上式即為Orr-Sommerfeld equation

3-2 建立代數特徵方程式

首先利用 Galerkin Method 消去時間項，再利用 Collocation

g y u D u i D t D ) Re ( ) Re[ ( ) ]} {( ₂ 2 2 2 2 2 2 2 2 ∂ ∂ − − − − ∂ ∂ − −α α α α =c iα Re(α2 −D2)g (3.2-1) 將式式（3.1-2）帶入上式可得 mn t im n m n a e y y u D u i D t D α α α α φ − ω ∂ ∂ − − − − ∂ ∂ − −

∑∑

{( 2 2)2 Re ( 2 2) Re[ ( 2 2) 2₂]} ( ) mn t im n m n a e y D i c

∑∑

α α − φ − ω = Re( 2 2) ( ) _(3.2-2) 首先利用 Galerkin Method，將（3.2-2）乘上 ip t e ω ，p=-M1~M1並對 每個週期進行積分，在積分前可將基態流與時間有關部份轉變為以下 型式： t i p t i p p u e u e u = ₁ ω + ₂ −ω 其中 } / )] ( Im[ )] ( [ { 2 1 1 realU y U y i u_p = − 、 { [ ( )] Im[ ( )]/ } 2 1 2 realU y U y i u_p = + ] 1 )[ Re 2 ( ) ( _{( )} ) ( β β β β β β β β ω ε i i y i y y i y e e e e i y U _{+ − +} + − + + + − − = 其次利用Collocation Method ，將空間項分別以 ) 1 N 2 I ( + =Cos π y_I ，I=1、 2、3 …N 帶入；經過上述運算後可得代數特徵方程式，表示為 BX c AX = 其中 c 為特徵值、A、B、X 與其矩陣內部單元表示如下：

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = − − M , M , 1 M , M , M , 1 M , 2 , 3 , 3 , 2 , 2 , 2 , 1 , 2 , 2 , 1 , 1 , 1 , 0 . 0 . . . 0 . . 0 . 0 . 0 J I J I J I J I J I J I J I J I J I A A A A A A A A A A ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = M , M , 2 , 2 , 1 , 1 , 0 . . 0 0 . . . . . . 0 0 . . 0 J I J I J I B B B B _、 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = J J J a a a X , M , 2 , 1 . . 其中： ) ( ] ) ( ) )( ( Re[ 2 2 1 2 2 2 1 1 -, , J I I p I p k k J I y y y u D y u i A α α φ ω π ∂ ∂ + − = ) Re( 2 ) ( 2 { 2 2 2 2 2 , , _ω α π α π _{− + −} = D ip D A_Ik_Jk 2 Re[ ( )( ) (₂ )]} ( ) 2 2 2 I J I s I s y y y u D y u iα α φ ω π ∂ ∂ + − + ) ( ] ) ( ) )( ( Re[ 2 2 2 2 2 2 2 1 , , J I I p I p k k J I y y y u D y u i A α α φ ω π ∂ ∂ + − = + ) ( ) Re( 2 2 2 , , J I k k J I i D y B α α φ ω π ₋ = ) 1 2 / ) 1 M (( − + − = k p 、k =1 , 2 ... M、I,J =1 , 2... N 以上所述的矩陣單元，使用Galerkin Method 處理時可將時間係數 分成三種情形運算 i m t e−( −1)ω、e−imωt及e−i(m+1)ωt，對每個週期進行積分後， 則可得到三種不同形式：

∫

⋅ − − ⋅ =_⎩⎨⎧ = _≠= ₊+ 1 , 0 1 , ) 1 ( p m p m A A dt e e A i m ωt ipωt

∫

⋅ − ⋅ =_⎩⎨⎧ = _≠= p m p m A A dt e e A im t ip t , 0 , ω ω

∫

⋅ − + ⋅ =_⎩⎨⎧ = _≠= ₋− 1 , 0 1 , ) 1 ( p m p m A A dt e e A i m ωt ipωt

3-3 求解特徵值

上一節中，吾人得到代數特徵方程式AX =c BX型式，其中A、 B均為複數矩陣型式，吾人在此可利用 Fortran IMSL 程式庫之中的

EIGZC 副程式或 MATLAB 程式語言中的 eig 指令求解廣義特徵值 c， 因為c=cr +ici在本文中吾人利用ci特性作為判斷基態流場穩定與否 之依據，其特性為當ci <0時擾動隨時間衰減，流場會趨於穩定，ci >0

時擾動隨時間成長，流場趨於不穩定，ci趨近於零時則介於兩者之間，

第 四 章

結 果 與 討 論

本章將從最基本的基態解及展開項數開始探討，接著討論無調制 (ε =0)與調制作用下臨界曲線的變化，最後將結果綜合歸納並與前人 所提及之現象加以比較。

4-1 調制作用與基態(Basic state)關係

在本文第二章數學模式中經由統御方程式吾人可求解出流場基 態解，在調制作用下基態解相較於非調制作用下增加了調制振盪部分

( )

y t up , 的影響，當調制振幅與頻率改變時流場行為受壓力調制影響而 產生不同變化，其速度(u)與位置(y)變化可由圖二與圖三描述，圖中 不同線型分別代表 T/6、2T/6、3T/6、4T/6、5T/6、T 其中 T 為調制 基態解之週期(T =2π / ω)，由圖可知壓力調制導致流速在單一週期 內產生前後振盪行為。圖二(a) ~ (d)為相同調制頻率(ω =0.0001)下調 制振幅(ε =1~0)對流場的影響，由圖可知當調制振幅愈大流速振盪範 圍也愈大如圖二(a)，當調制振幅逐漸減小時振盪範圍也逐漸縮小如圖 二(c)，圖二(d)則為無調制作用此時速度分佈與時間無關，最大速度 為1。圖三(a) ~ (d)為相同調制振幅(ε =1)下調制頻率(ω=0.0001~0.1) 對流場的影響，由圖可知調制頻率愈小流場振盪範圍愈大如圖三(a)，

當調制頻率逐漸變大時，振盪範圍逐漸縮小，最後流場趨於無調制狀 態如圖三(d)，此結果與 Straatman(2002)之基態解行為一致。

4-2 展開項數的影響

4-2-1 無調制(

ε

=

0

)狀態下展開項數

由Orszag(1971)、Fortin(1994) . . . 等人的結果可知在無調制平面 普修流中臨界雷諾數(Re₀ =5772.22)、臨界波數(α =1.020545)的精確 結果。在前一章中吾人將流線函數以雙重級數形式展開，其展開項數 (M,N)將影響數值計算結果的準確度，因此吾人將藉由改變展開項數 來尋找準確的臨界雷諾數。 在無調制狀態下(ε =0)取 M = 1、N=12 , 13 , 14 . . .項做項數展開 測試其結果如圖四所示，由圖可知當N = 22 項時結果趨於Re₀。另一 方面，測試結果顯示不論M 取 1 , 3 , 5 . . .均不影響測試結果，所以在 無調制狀況下只受展開項數 N 的多寡影響，展開項數越多則越趨近 於 Re₀，由以上結果本文將固定取 N=22 項做為往後計算之基準。表 一為Orszag(1971)與本文在臨界值條件下項數(N=20 , 26. . . )所對應 之波速(c=cr +ici)比較，結果相當一致。

4-2-2 調制狀態下展開項數

經過上一節討論N 取 22 項，用來測試不同調制參數下應取之 M

至3 項，相同情況下ω減小時則須使用較多項數才具有足夠的計算精 度，例如ω=0.00001為保持誤差小於 5 % 則 M 必須取 13 項以上， 表中err(%)= (R_M -R_M-2)/R_M ×100%為誤差定義其中R_M為取 M 項展 開求得之雷諾數，R_M-2為取M-2 項求得之雷諾數。

4-2-3 調制狀態下極值區域展開項數結果

由Davis (1976)的研究可知，當調制流場在極高或極低頻率時臨 界雷諾數會趨於一固定值，在極高頻率時會趨於無調制狀態Rec ≈Re₀， 極低頻率時會有Rec ≈Re₀/(1+

ε

)關係，以ε =1為例極高與極低頻時會 分別趨於 5772.22 與 2886.11，在此嘗試逼近極低頻時理論值所需 M 項數，其結果如圖五所示，橫軸為項數(M)縱軸為誤差(error % )其中 誤差定義error(%)= (Rec_M -2886)/2886 ×100%其中Rec_M 為取 M 項之

臨界雷諾數，由圖可知 M = 17 以後趨於平緩誤差小於 2%，結果與

Davis 所述相同。

4-3 臨界曲線(Marginal stability curve)

4-3-1 無調制狀態下之臨界曲線

在平面普修流穩定性分析中，無調制情況下（ε =0）所求解之方 程式為著名之 Orr-Sommerfeld 方程式，其中判定流場穩定或不穩定 的重要數值結果可整理成圖六，在圖六中縱軸為雷諾數(Re)橫軸為波 數(α)，圖中符號為計算所得之值，實線為符號數值經由曲線擬合後 求得，此實線即稱為臨界曲線，由圖中可知臨界曲線之最低點為臨界

雷諾數(Re_c ≈5772.22)與臨界波數(

α

_c ≈1.02)，此結果與前人一致，在 本文中此為調制平面普修流之特例。

4-3-2 調制狀態下之臨界曲線

在調制情況下臨界曲線會因調制參數不同而產生偏移。圖七(a) ~ (f)為固定調制振幅(ε =1)下不同ω所對應之臨界曲線，由圖可知調制 作用下之臨界曲線形狀與無調制作用時相似，臨界波數範圍分佈於 0.9 ~ 1.2 間，在低頻與高頻時會趨於 1.02。臨界雷諾數則在 2886.11 ~ 5772.22 間，低頻與高頻分別趨近 2886.11 與 5772.22。

4-3-3 調制參數對臨界曲線之影響

圖八為固定調制振幅(ε =1)下對不同調制頻率計算臨界曲線的變 化，並將其結果依ω範圍加以疊加成圖八(a)、(b)、(c)，圖中橫軸為 波數(α )縱軸為雷諾數( eR )，虛線部分為無調制狀態，由圖八(a)可知 在ω =3×10-4 ~10-3區間，臨界曲線位置隨調制頻率增加而上升；由圖 八(b)可知在_{ω =2×10}-3 _~5×10-3_{區間，臨界曲線位置隨調制頻率增加} 而下降；由圖八(c)可知在ω=0.01~0.03區間，臨界曲線位置隨調制頻 率增加而上升；綜合圖八結果可知臨界雷諾數愈小臨界曲線包含區域 則愈大，各臨界曲線均在無調制臨界曲線下。圖九為ω =0.01時不同ε 對臨界曲線之影響，由圖中得知當ε上升會使臨界雷諾數下降，又因 0.01 = ω 落於非極高或極低頻區域，臨界波數非固定為 1.02，由圖可 ε

其中橫軸為調制頻率(ω)縱軸為臨界雷諾數(Rec)，圖中標明不同調制 振幅(ε =1、0.5、0.3、0.1)與相對應理論值(2886、3848、4440、5247)， 圖中可知當調制振幅減小時曲線向上提升，振幅越小則趨於無調制； 在調制頻率極大與極小時臨界雷諾數分別趨於理論值；圖十所對應之 相關數據整理於表三。圖十一為臨界波數與頻率關係圖，相關數據整 理於表四，圖中可知極高與低頻區域臨界波數趨於 1.02，中頻區域則 依不同頻率而變化，此區域可看出臨界曲線偏移現象較明顯。

4-4 流線函數擾動量隨時間變化關係

依據第三章數值方法擾動量的定義，吾人在此將計算後之特徵值 代入，選取流場中一點對流線函數擾動量隨時間變化做分析，結果如 圖十二 ~ 十四，圖中橫軸為時間，縱軸為擾動量(ψ )，由表三可知 在ω=0.001、ε =1所對應之臨界雷諾數為4005，圖十二(a)、(b)、(c) 分別代表小於臨界雷諾數、等於臨界雷諾數、大於臨界雷諾數擾動量 隨時間的變化，由圖十二(a)可知擾動量隨時間呈週期性衰減，圖十二 (b)擾動量隨時間呈週期性振盪，圖十二(c)擾動量隨時間呈週期性成 長。在圖十三、十四則代表不同調制參數下之結果，趨勢均與圖十二 相同，經由以上分析得知圖面結果與數學模式推導一致。

4-5 調制頻率、振幅範圍與實際心跳參數

動物的血管大致可區分為動脈、靜脈與微血管三種，其直徑大小

因部位及功能而不同。由表五可概估動物實際心跳頻率約1 ~ 3 HZ， 一般而言血管直徑大小約

γ

×10-3m，γ 為尺寸因子常數依部位而不同 (人主動脈直徑約 52. ~ 53. cm，γ = 25 ~ 35)，流速由 Nerem(1972) 心臟實驗得知大約0 ~ 200cm/s (人動脈約 1 m / s；微血管 > 0.3 mm / s)，故由 2-2 節無因次定義_ω* _{=ω H/u0可概估大主動脈所對應之無因} 次頻率約落於圖十0.025 ~ 0.105 之間，同理代入不同

γ

值可估出其他 對應之頻率範圍。其次由 Nerem 實驗速度與時間關係圖可概估心跳 振幅約 0.6，經過以上分析可知，在此調制範圍內不穩定皆會較無調 制作用時提前發生。

4-6 中頻區域

Davis(1976) 將各種不同調制型式問題結果做整理，大致上可分 為平行剪應力流(Parallel shear flow)、溫度不穩定(Thermal instability) 與離心不穩定(Centrifugal instability)三類，平面普修流屬於平行剪應 力流的一種，但相較於其他調制形式問題(如溫度不穩定)對於非極大 或極小調制頻率部分目前仍無較明確的了解與認識，本文圖十所示結 果在極高或極低頻部分均趨於理論值，但在ω=10-4~10-2區間所產生 之跳起現象卻與其他已知調制形式的趨勢不同，在 Straatman(2002) 結果中特定頻率(Wo=1.107,_ω= Wo2/Re_{)附近也出現未知現象無法求} 得 數 值 解 ， 作 者 猜 測 可 能 有 奇 異 點 存 在 於 區 間 內 ， 圖 十 五 為 Straatman(2002)與本文結果之比較，由圖可知高頻區域與本文相同且

現象或許與本文圖十之跳起現象有關，雖然吾人欲找出原因但從各相 關參數中探索與檢驗下仍無法找出明確關係以及對應之物理現象，或 許應該利用其他方法或數值模型才能更進一步了解未知的調制現 象。

第 五 章

結 論

本文透過數值解探討調制參數對平面普修流之穩定性影響，綜合 上一章結果吾人可得以下結論： 1. 在無調制（

ε

=0 ）作用時展開項數受項數N 影響，當展開項數為 (1×22)時得到Re₀ =5772.22、

α

c =1.02056經比較特徵值後其結果 與Orszag(1971)、Fortin(1994) . . . 等人一致。 2. 數值計算結果包含 Nerem(1972)心臟實驗概估之實際調制振幅與 頻率範圍，且調制作用會導致流場更加不穩定。 3. 調制作用下展開項數 N=22 項，ω越小 M 項數則須越大，且臨界 雷諾數均低於無調制(Re₀ =5772.22)，與 Straatman(2002)一致。 4. 調制作用下流場在極高或極低頻率時臨界雷諾數會趨於一固定值 ， 極 高 頻 時 會 趨 於 無 調 制 狀 態 Rec ≈Re₀ ， 極 低 頻 時 會 有 ) 1 /( Re Re_c ≈ ₀ +

ε

關係，結果與 Davis (1976)一致。 5. 調制作用會導致Rec、αc改變，使臨界曲線產生偏移，其次若Rec 愈大則所涵蓋不穩定區域則愈小這些現象與行為均與 Grosch & Salwen(1968) 一致。 6. 由Rec與ω關係圖及擾動量(

ψ

)與時間關係圖可知在固定ω、ε、α

當Re≈Rec時擾動量隨時間不發散亦不衰減，流場呈週期性振盪；

參 考 文 獻

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

Chen Wen-Yau and Chen, C. F. 1999 Effect of gravity modulation on the stability of convection in a vertical slot. J. Fluid Mech. 395, 327.

Criminale, W. O., Jackson, T. L. & Joslin, R. D. 2003 Theory and Computation in Hydrodynamic Stability.

Davis, Stephen H. 1976 The stability of time-periodic flow. Annu. Rev.

Fluid Mech. 7, 57.

Drazin, P. G. & Reid, W. H. 1981 Hydrodynamic Stability. Cambridge University Press.

Fortin A. and Jardak M., 1994 Old and new results on the two - Dimensional Poiseuille flow. J. Computational Physics. 115,455. Fedele Francesco, Hitt Darren and Prabhu Rachakonda, 2003

Revisiting the stability of pulsatele pipe flow. European Journal

of Mechanics B/Fluids. 24, 237.

Grosch Chester E. and Salwen Harold, 1968 The syability of steady and time-dependent plane Poiseuille flow. J. Fluid Mech. 34, 177.

Herbert, D. M. 1972 The energy balance in modulated plane Poiseuille flow. J. Fluid Mech. 56, 73.

Lin, C. C., 1955 The Theory of Hydrodynamic Stability. Cambridge Monographs on Mechanics and Applied Mathematics.

Mcintire, L. V. & Lin, C. H. 1972 Finite amplitude instability of second-order fluids in plane Poiseuille flow. J. Fluid Mech. 52, 273.

Nishioka, M., and Ichikawa, Y. 1975 An experimental investigation of the stability of plane Poiseuille flow. J. Fluid Mech. 72, 731.

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

Pekeris, C. L. & Shkoller , B. 1967 Stability of plane Poiseuille flow to periodic disturbances of finite amplitude in the vicinity of the neutral curve. J. Fluid Mech. 29, 31.

Nerem, R. M., Seed, W. A. and Wood, N. B. 1972 An experimental study of the velocity distribution and transition to turbulence in the aorta, J. Fluid Mech. 52, 137.

Sarpkaya, T. 1966 Experimental determination of the critical Reynolds number for pulsating Poiseuille flow, J. Basic Eng. 88, 589.

Straatman, A. G., Khayat, R. E., and Haj-Qasem, E. 2002 On the hydrodynamic stability of pulsatile flow in a plane channel.

Physics of Fluids. vol. 14, No. 6.

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

91, 780.

Wood, N. B. 1999 Aspects of fluid dynamics applied to the large arteries, J. Theor. Biol. 199, 137.

楊文美，「黏彈性流體在調節重力場下的穩定」，行政院國家科學委

楊慶泰，「圓柱間調制庫耶流之穩定性」，國立交通大學機械工程研 究所碩士論文，民國84 年。 李鎰清，「圓柱之流場-一種改良 Lorenz 模型」，國立交通大學機械 工程研究所碩士論文，民國94 年。 洪英棋，「平面普修流之非線性不穩定」，國立交通大學機械工程研 究所碩士論文，民國94 年。 廖啟豪，「調制平面普修流之非線性不穩定」，國立交通大學機械工 程研究所碩士論文，民國95 年。

表一 無調制結果與 Orszag(1971) 項數結果比較 N cr + i ci 20 0.26400173436 + 0.000000050160i 0.26400165486 + 0.000000066646i 26 0.26400173957 - 0.000000003021i 0.26400166008 + 0.000000013465i 38 0.26400173958 - 0.00000000302i 0.26400166008 + 0.00000001347i 44 0.26400173958 - 0.00000000301i 0.26400166012 + 0.00000001352i 50 0.26400173944 - 0.00000000306i 0.26400165986 + 0.00000001344i Orszag(1971) 本文計算結果

表二 當 ε=1 時調制頻率與 M 項數關係 M\ω 0.00001 0.0001 0.001 0.01 0.1 1 3 3381.4 3394.1 4387.1 3159.6 5761.3 5772.3 5 3093.5 3128.76 4158.29 3055.9 5761.27 5772.3 7 3000.4 3011.97 4074.52 3055.6 --- --- 9 2958.6 2969.95 4035.22 --- --- --- 11 2936.3 2947.45 4015.13 --- --- --- 13 2922.9 2934 4005 --- --- --- 15 2914.2 --- --- --- --- --- 17 2908.3 --- --- --- --- --- 19 2904.1 --- --- --- --- --- 21 2901.0 --- --- --- --- --- M 項數選擇 (條件 : err < 5 %) 13 13 13 5 3 3

表三 不同 ω、ε 所對應之臨界雷諾數與項數(M)關係 N=22\M 13 13 13 13 13 13 13 11 5 3 3 3 3 ω 0.00001 0.0000316 0.0001 0.000316 0.001 0.002 0.00316 0.00562 0.01 0.0316 0.1 0.316 1 ε=1.0 2923 2924 2934 3036 4005 4898 4003.7 3180.2 3055.9 5441 5761 5772.3 5772.3 ε=0.5 3880.8 3882.3 3894 4047.3 4937.6 5503.3 5008.3 4671.8 5103.6 5693.8 5769.5 5772.3 5772.3 ε=0.3 4466.2 4467.9 4484.4 4632.1 5299.1 5648.2 5381.8 5330.9 5547.4 5744.1 5771.3 5772.3 5772.3 ε=0.1 5259.6 5260.8 5271.8 5360 5634.6 5747.3 5705.3 5722.7 5747.9 5769.1 5772.1 5772.3 5772.3 ε=0.0 5772.3 5772.3 5772.3 5772.3 5772.3 5772.3 5772.3 5772.3 5772.3 5772.3 5772.3 5772.3 5772.3 30

表四 不同 ω、ε 所對應之臨界波數(α)值 ω 0.00001 0.0000316 0.0001 0.000316 0.001 0.002 0.00316 0.00562 0.01 0.0316 0.1 0.316 1 ε=1.0 1.02 1.02 1.02 1.02 1.03 1.02 1.04 1.07 1.1 1.03 1.02 1.02 1.02 ε=0.5 1.02 1.02 1.02 1.02 1.02 1.02 1.02 1.04 1.04 1.02 1.02 1.02 1.02 ε=0.3 1.02 1.02 1.02 1.02 1.02 1.02 1.02 1.03 1.03 1.02 1.02 1.02 1.02 ε=0.1 1.02 1.02 1.02 1.02 1.02 1.02 1.02 1.02 1.02 1.02 1.02 1.02 1.02 31

表五 各種動物實際心跳頻率整理(每分鐘計) 動物 心跳頻率 小白鼠 175 兔 150 ~ 210 貓 102 ~ 140 狗 90 ~ 100 綿羊 70 ~ 110 人 72 牛、馬 40 象 25

 

2H −∂P/∂x =C (1+

ε

Cos

ω

t) y

x

  (a) (b) 圖二 相同調制頻率下(ω=0.0001)不同調制振幅對基態之影響 (a) ε=1 (b) ε=0.5 (c) ε=0.25 (d) ε=0 0 0.5 1 1.5 2 2.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u y T/6 2T/6 3T/6 4T/6 5T/6 T 0 0.5 1 1.5 2 2.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u y T/6 2T/6 3T/6 4T/6 5T/6 T

  (c) (d) 續 0 0.5 1 1.5 2 2.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u y T/6 2T/6 3T/6 4T/6 5T/6 T 0 0.5 1 1.5 2 2.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u y ε=0

  (a) (b) 圖三 相同調制振幅下(ε=1)不同調制頻率對基態之影響 (a) ω=0.0001 (b) ω=0.001 (c) ω=0.01 (d) ω=0.1 0 0.5 1 1.5 2 2.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u y T/6 2T/6 3T/6 4T/6 5T/6 T 0 0.5 1 1.5 2 2.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u y T/6 2T/6 3T/6 4T/6 5T/6 T

  (c) (d) 續 0 0.5 1 1.5 2 2.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u y T/6 2T/6 3T/6 4T/6 5T/6 T 0 0.5 1 1.5 2 2.5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u y T/6 2T/6 3T/6 4T/6 5T/6 T

  圖四 無調制作用時 M=1，不同 N 項數與臨界雷諾數關係測試圖 10 12 14 16 18 20 22 24 26 28 30 5740 5750 5760 5770 5780 5790 5800 5810 5820 5830

N

Re

c

  圖五 當ω=0.00001、ε=1 時，N=22 不同取項 M 與誤差關係圖 3 5 7 9 11 13 15 17 19 21 -5 0 5 10 15 20

M

error(%)

  圖六 無調制平面普修流(ε=0)α與 Re 關係圖 0.9 0.95 1 1.05 1.1 5600 5800 6000 6200 6400 6600 6800 7000 7200 7400

Re

α

  (a) (b) 圖七 當ε=1 時不同調制頻率下α與 Re 關係圖(a) ω=0.00001 (b) 0.9 0.95 1 1.05 1.1 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 Re α 0.9 0.95 1 1.05 1.1 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 Re α

  (c) (d) 續 0.9 0.95 1 1.05 1.1 3800 4000 4200 4400 4600 4800 5000 5200 5400 5600 Re α 0.9 0.95 1 1.05 1.1 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 Re α

  (e) (f) 續 0.9 0.95 1 1.05 1.1 1.15 1.2 3000 3500 4000 4500 5000 5500 Re α 0.9 0.95 1 1.05 1.1 5600 5800 6000 6200 6400 6600 6800 7000 7200 Re α

  (a) (b) 圖八 當ε=1 時不同ω範圍對臨界曲線影響圖(a) 0.0003≤ω≤0.001  (b) 0.002 ≤ω≤0.005 (c) 0.01≤ω≤0.03 0.9 0.95 1 1.05 1.1 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 Re α Steady ω=0.0003 ω=0.0005 ω=0.0007 ω=0.001 0.9 0.95 1 1.05 1.1 1.15 3000 3500 4000 4500 5000 5500 6000 6500 7000 Re α Steady ω=0.002 ω=0.003 ω=0.004 ω=0.005

  (c) 續 0.9 0.95 1 1.05 1.1 1.15 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 Re α Steady ω=0.01 ω=0.015 ω=0.017 ω=0.02 ω=0.03

  圖九 相同調制頻率下(ω=0.01)不同調制振幅對臨界曲線之影響 0.9 0.95 1 1.05 1.1 1.15 1.2 3000 3500 4000 4500 5000 5500 6000 6500 7000 7500 Re α Steay ε = 1 ε = 0.7 ε = 0.5 ε = 0.3

 

 

  (a) (b) (c) 圖十二 當ω=0.001、ε=1、Rec =4005，討論流場中(x,y)=(0.5,0.5)處流 線函數擾動量隨時間變化關係圖(a)Re=3980(b)Re=4005(c) 025 4 Re= ψ ψ ψ

  (a) (b) (c) 圖十三 當ω=0.00316、ε=1、Rec = 4003.7，討論流場中(x,y)=(0.5,0.5) 處 流 線 函 數 擾 動 量 隨 時 間 變 化 關 係 圖 (a) Re=3950 (b) 7 . 4003 Re= (c)Re=4050 ψ ψ_ψ′ ψ

  (a) (b) (c) 圖十四 當ω=0.01、ε=1、Rec =3055.9，討論流場中(x,y)=(0.5,0.5)處流 線函數擾動量隨時間變化關係圖(a) = = ψ ψ ψ