中 華 大 學

中 華 大 學 碩 士 論 文

題目：脈衝震波對衝擊平板熱傳效應之數值研究

系 所 別：機械與航太工程研究所 學號姓名：M09508009 陳 永 軒 指導教授：鄭 藏 勝 博 士

中華民國 九十七 年 七 月

中文摘要

本文探討超音速脈衝噴流系統的物理現象，脈衝噴流與衝擊流 不同的地方在於脈衝噴流流出震波管時為震波波形，可應用在超音速 飛行器的引擎排氣管和火箭發射器等等。

當脈衝噴流從震波管流出衝擊加熱平板後，其震波與平板之交互 作用、反射震波對流場結構及平板熱傳效應的影響都是值得深入探討 的物理問題，因此，吾人對流場無平板、絕熱平板及加熱平板三種脈 衝噴流之流場結構進行數值模擬，對於脈衝噴流之流場結構，可觀察 到脈衝噴流發展大致可分為三個階段，第一階段為噴流流出震波管後 在出口端產生震波繞射，第二階段為形成一非穩定之馬赫圓盤，第三 階段為形成第一個震波囊結構。

為了瞭解此脈衝噴流衝擊加熱平板後對平板熱傳效應之影響，可 探討紐塞爾數的變化情形，而影響紐塞爾數的主要參數為溫度梯度，

當平板至震波管距離增加時溫度梯度也增加，紐塞爾數分佈現象會呈 現平板溫度越大其塞爾數越大的情形。在固定壁面溫度的情況下，平 均紐塞爾數皆隨著震波管至平板距離的增加而增加，當壁面溫度T_W = 2T₁時，其平均紐塞爾數增加的速度比其它壁面溫度要快。

致 謝

誠 摯 的 感 謝 我 的 指 導 教 授 鄭藏 勝 博 士 ， 令 我 在 數 值 模 擬 部 分 這 塊 領 域 上 瞭 解 到 其 要 點 ， 且 不 時 的 討 論 並 指 點 我 正 確 的 方 向 ， 使 我 在 兩 年 中 獲 益 匪 淺 ， 老 師 對 學 問 的 嚴 謹 更 是 我 學 習 的 指 標 。 本 論 文 的 完 成 另 外 亦 得 感 謝 實 驗 室 內 學 弟 柳 文 祥 的 協 助 ， 亦 得 感 謝 從 大 學 直 到 研 究 所 都 在 一 起 的 朋 友 們 ， 盧 韋 廷 、 高 迺 迪 、 曹 鈞 鴻 、 侯 忠 慶 、 葉 建 宏 、 蕭 世 昌 、 陳 炯 漢 ， 同 輩 分 的 高 祥 恩 、 還 有 研 究 所 甲 班 熱 流 組 的 同 學 及 學 弟 ， 給 我 學 術 上 指 導 的 學 長 ， 黃 博 丞 、 李 儀 威 、 陳 家 政 、 劉 光 倫 、 張 遼 、 周 中 祇 、 張 立 嚴 、 柯 志 偉 ， 還 有 在 我 疲 累 時 ， 陪 我 一 起 打 球 的 中 華 大 學 排 球 隊 ， 裡 面 有 很 多 我 所 珍 惜 的 回 憶 ， 兩 年 過 去 了 ， 實 驗 室 裡 共 同 的 生 活 點 滴 ， 佔 了 日 常 生 活 中 的 大 部 分 ， 感 謝 你 們 。 最後，將本文獻給我最敬愛的家人還有我的雙親，感謝您無怨無悔的 養育與無時無刻的關懷照顧，還有經濟上與精神上的支持，讓我能專 注於課業研究中，謝謝我所認識的所有人。

目 錄

中文摘要... i

致 謝 ... ii

目 錄 ... iii

圖目錄 ... v

符號說明... vii

希臘符號...xiii

第一章 緒論... 1

1-1 前言 ... 1

1-2 文獻回顧 ... 2

1-3 研究動機 ... 6

第二章 物理問題 ... 7

2-1 物理模式 ... 7

2-2 基本假設 ... 7

第三章 數值方法 ... 9

3-1 統御方程式 ... 9

3-2 時間積分 ... 12

3-3 空間差分 ... 13

3-4.1 起始條件 ... 18

3-4.2 邊界條件 ... 18

3-5 局部平均紐賽爾數及紐賽爾數 ... 20

第四章 結果與討論 ... 21

4-1 無平板之脈衝噴流結構 ... 21

4-2 格點測試 ... 23

4-3 絕熱平板之脈衝噴流結構 ... 24

4-4 平板熱傳效應之影響 ... 26

4-5 相同震波管至平板距離及不同平板溫度對紐塞爾數之影響 ... 27

第五章 結論... 29

參考文獻... 31

圖目錄

圖 2-1 二維軸對稱脈衝噴流衝擊平板模擬區域示意圖... 33 圖4-1 在無平板時之脈衝噴流流場示意圖[12] ... 34 圖4-2 在 t = 250μs 之密度與壓力等位線圖。(a)PE/P1 = 2.2；(b) PE/P1 = 3.61；(c) PE/P1 = 5。 ... 35 圖4-3 在不同震波管出口壓力下之非穩定馬赫圓盤的直徑( DM )和位

置( XM )圖。 ... 36 圖4-4 格點測試圖。Grid1: ΔX = Δr = 0.025D；Grid2: ΔX = Δr =

0.00125D；Grid3: ΔX = Δr = 0.001D；Grid4: ΔX = Δr = 0.0075D。 ... 37 圖4-5 在震波管至平板距離 x/D = 2.5、時間 t = 4.7 ~ 14 且噴嘴壓力

比PE/P1 = 5 時，震波衝擊絕熱平板後之流場結構數值紋影圖。

... 38 圖4-6 在相同噴嘴出口壓力比( PE/P1 = 5 )，不同平板至震波管距離

( x/D = 2.0 ~ 4.0 )，相同無因次時間( t = 10 )，相同平板溫度( TW

= 8T1 )時之溫度等位線圖。(a) x/D = 2.0；(b) x/D = 2.5；(c) x/D

= 3.0；(d) x/D = 3.5；(e) x/D = 4.0。 ... 39 圖4-7 在相同噴嘴出口壓力比( PE/P1 = 5 )，平板至震波管距離為 x/D = 3.4 ~ 3.5，相同無因次時間( t = 10 )，不同平板溫度( TW = 2T1 ~

2T1；(b) TW = 4T1；(c) TW = 6T1；(d) TW = 8T1；(e) TW = 10T1。 ... 40 圖4-8 在相同噴嘴壓力比( PE/P1 = 5 )下，相同平板至震波管距離，不

同平板溫度( TW = 2T1 ~ 10 T1 )下之局部紐塞爾數分佈圖。(a) x/D = 2.0；(b) x/D = 2.5。... 41 圖4-9 在相同噴嘴壓力比( PE/P1 = 5 )下，不同平板至震波管距離(x/D

= 2.0~4.0)，相同平板溫度( TW = 2T1 )之局部紐塞爾數分佈圖。

... 44 圖4-10 在相同噴嘴壓力比( PE/P1 = 5 )下，改變平板至震波管距離及

平板溫度( TW = 2T1 ~ 10T1 )之平均紐塞爾數分佈圖。 ... 45

符號說明

D ： 震波管直徑 K ： 熱傳係數 Ms : 震波馬赫數

Nu ： 紐塞爾數(Nusselt Number) P ： 氣體壓力

Pr ： 普蘭特數(Prandtl Number) R ： 震波管半徑

Ra ： 瑞理數(Rayleigh Number) S ： 黏滯係數之 Sutherland 常數 S1 ： 熱傳係數之 Sutherland 常數 T1 ： 參考溫度

T2 ： 震波進入之溫度 TW ： 平板溫度

u ： x 方向之速度分量 v ： r 方向之速度分量

希臘符號

γ ： 比熱值，γ =1.4

μ ： 黏滯係數

第一章 緒論

1-1 前言

衝擊流( Impinging Jet )系統具有極高的熱傳與質傳速率，利用高 速噴出的流體衝擊物件表面，可以達到快速加熱、冷卻或乾燥的目 的，因此，在工業上被廣泛的應用於熱金屬或電子元件的冷卻、紙張 或紡織品的乾燥、渦輪葉片或引擎燃燒室的散熱等等。此外，在火箭 和飛彈發射時，由噴嘴噴出的高溫氣體直接衝擊地表或防護材，也是 衝擊流的一種例子。而一般的衝擊流系統大多在探討次音速下的物理 現象，但在超音速下脈衝噴流( Supersonic Pulse Jet )系統的物理現象 則較少人研究。

本文主要是探討超音速脈衝噴流系統的物理現象，脈衝噴流與 衝擊流不同的地方在於脈衝噴流流出震波管時為震波波形，而此脈衝 震波具有高速高溫高壓的特性，可應用在超音速飛行器的引擎排氣管 和火箭發射器等等，因此，當脈衝噴流衝擊平板時，震波管至平板的 距離以及平板在不同溫度下，對震波結構與平板熱傳率之影響是值得 深入探討的問題。

1-2 文獻回顧

為了分析脈衝噴流流出震波管後非理想膨脹震波撞擊平板後反 射，震波與平板的交互作用、反射震波對震波流場的影響以及平板在 不同溫度下的熱傳效應，將參考文獻分為熱傳和震波兩部分探討：第 一部分為熱傳分析的研究，第二部分為震波模擬的研究。

熱傳分析部分：Heyerichs 及 Pollard [1]以紊流模式進行衝擊流 之熱傳研究，其所使用之 Wilcox 紊流模式包含 k-ε 和 k-ω 兩種，且 使用控制體積法和 QUICK differencing 法求出平均瑞理數和能量方 程，其紊流模式預測紐塞爾數跟實驗相比誤差不超過2%，且收斂時 間快 24%，證實 k-ω 模式可以預測出較複雜的熱傳現象。Rahimi 及 Owen [2]以軸對稱不足膨脹噴流進行熱傳分析的實驗，其噴嘴至衝擊 物件表面間隔 3、6 和 10 個單位噴嘴直徑下，以震波管壓力比 3.04 和 5.08 的紐塞爾數來做比較，發現當物件表面溫度較高時膨脹區域 內徑向的紐塞爾數會變高，表示其熱傳速率較好。Lee 等人[3]以漩渦 噴流衝擊平板表面進行熱傳的實驗分析，發現在噴流具有漩渦時之平 均紐塞爾數大於沒有漩渦時之平均紐塞爾數，其紐塞爾數最大值發生 在漩渦數為S=0.21 且噴嘴至平板距離為 L/d = 2 時，而紐塞爾數最小 值發生在漩渦數為S=0.77 且噴嘴至平板距離為 L/d = 2 時。Chung 及 Luo [4]以數值方法對非穩定之衝擊流進行熱傳分析，在空間上使用六

階精度的有限差分法，而時間上使用三階精度的 Runge-Kutta 積分 法，研究發現當平板溫度不變且雷諾數越高時，其紐塞爾數越高，而 在漩渦處的紐塞爾數會比其他位置高。Ramanujachari 等人 [5]以超音 速衝擊流衝擊垂直平板孔洞進行熱傳的實驗分析，為了探討在平板孔 洞的熱傳變化，改變空氣速度 60~210m/s、雷諾數 3×10⁴~1×10⁵和噴 嘴至平板距離87.5~375mm 做分析，觀察衝擊流在平板孔洞的紐塞爾 數的變化，發現愈靠近中心點的孔洞其熱傳速率愈好。Juhany 及 Hant[6]進行超音速薄膜冷卻的實驗研究，空氣和氦氣以馬赫數 Ms = 1.3 和 Ms = 2.2 被注入在一馬赫數 Ms = 2.4 的氣流裡，測量平板上垂 直點壓力跟溫度的分佈，發現空氣在馬赫數Ms = 1.3 時其冷卻效果最 好。

震波模擬部分：Jiang 等人[7]以數值模擬和實驗方法研究震波由 震 波 管 進 入 膨 脹 管 時 的 震 波 結 構 變 化 ， 數 值 方 法 是 採 用 Dispersion-Controlled 法，而流場結構為軸對稱型式，其主要在探討 不同馬赫數Ms = 1.3 與 Ms = 1.5 下，震波與剪切層的交互作用以及反 射震波與渦卷( Vortex )所形成的流場結構，發現當馬赫數越大時，即 膨脹係數越大，其震波流出震波管時剪切層也越大，同時反射震波和 馬赫反射震波會跟震波結合，使震波流場結構產生多處渦卷。Teng

型震波聚焦在圓柱管內的研究，而環型震波有兩個噴流出口分別在管 的上方跟下方，兩噴流同時流出且在管中間區域結合，發現環型震波 會有兩個聚焦點並且有聚焦、繞射和反射三種震波結構，且環形震波 在聚焦點可以在百萬分之一秒達到高溫高壓。Liang 等人[9]以 WENO 數值方法解尤拉方程式探討震波衝擊平板後，反射震波與渦卷的交互 作用，計算之馬赫數範圍介於1.2 ≦ Ms ≦ 2.0，且以 90 度及 180 度 繞射為例子，利用數值紋影效果的繪圖方法來觀察震波結構，研究結 果發現當震波馬赫數大於或等於1.4 時，震波管下游呈現一停滯反射 震波，而當震波馬赫數小於1.4 時，反射震波最終將反射回震波管內。

Sun 及 Takayama[10]利用數值方法進行震波 90 度繞射的模擬研究，

計算結果發現使用較密網格解尤拉方程式時，其剪切層會產生與實驗 結果不符合的小渦卷，若使用較疏網格雖不會在剪切層產生小渦卷，

但無法預測實驗所觀察到的第二個渦卷( Secondary Vortex )。而解層 流 Navier-Stokes 方程式雖可預測第二個渦卷卻無法抑制剪切層小渦 卷的產生，但若將方程式加上 k-ε 紊流模式則可得到與實驗相符的結 果 。Kim 及 Park [11]利用 Modified Low-Diffusion Flux-Splitting Scheme( MLDFSS )數值方法研究超音速噴流衝擊平板後之振盪行 為，計算變化的參數包括噴嘴至平板距離( H/D = 2~4 )及噴嘴壓力比 ( PE/P1 =2~5 )，研究結果發現無論變換噴嘴至平板的距離或變換噴嘴

壓 力 比 ， 平 板 表 面 壓 力 之 振 盪 頻 率 皆 具 有 層 級 行 為( Staging behavior )。Ishii 等人[12]以實驗和數值方法進行脈衝噴流的研究，其 主要在探討不同震波管壓力比下，脈衝噴流流出後的震波結構，以及 馬赫圓盤直徑與軸向位置隨時間之變化情形，研究結果發現在各種不 同震波管壓力比下，馬赫圓盤皆呈現不穩定現象，也就是馬赫圓盤的 直徑及軸相位置隨著時間而變化，同時數值模擬亦準確的計算此結 果。

由以上文獻得知衝擊流多數是在探討次音速下衝擊流對平板熱 傳效應的影響，甚少探討在超音速狀態且具有震波之噴流對熱傳效應 的影響，而具有震波之噴流則大多著重在噴流結構的研究，甚少探討 脈衝噴流衝擊平板後之震波與平板的交互作用，以及平板的熱傳效應 和反射震波對震波流場的影響。本文則是以二維軸對稱納維爾-史托 克方程式( Navier-Stokes equation )來模擬超音速脈衝噴流的震波流場 結構，改變震波管至平板的距離和改變不同的平板溫度，探討此震波 流場的結構並分析其平板的熱傳效應。

1-3 研究動機

由以上文獻得知大部分對震波的研究，多在探討震波跟渦旋的 交互作用和不同震波管出口壓力比對震波流場的影響，甚少討論當脈 衝噴流衝擊平板時，平板與震波管間之距離以及平板在不同溫度下，

對震波結構與平板熱傳速率之影響，如果能對衝擊平板的熱傳效應有 正確的分析，則對脈衝噴流系統物理現象的瞭解將會有極大幫助。

第二章 物理問題

2-1 物理模式

本 文 是 以 二 維 軸 對 稱 納 維 爾-史脫克方程式 ( Navier-Stokes equation )模擬脈衝噴流衝擊平板後之震波結構，模擬區域示意圖如圖 2-1 所示。其中為R震波管半徑，L為震波管出口到平板的距離，PE、 TE及 ρ^Ｅ分別為震波管出口之壓力、溫度及密度，Ｐ₁、Ｔ１及 ρ^１分別 為大氣之壓力、溫度及密度，本研究探討當震波管至平板距離Ｌ = 2D~4D 及平板溫度為絕熱與 TW = 2T1~10T1時，對震波與平板的交互 作用、反射震波與渦卷( Vortex )交互作用及平板熱傳率之影響，其中 T1 = 289K。

2-2 基本假設

統御方程式主要是以納維爾-史脫克方程式為依據，對流場的特 性，作了下列幾點的假設：

(1) 流體為理想氣體。

(2) 為一黏性流體。

(3) 為一可壓縮流場。

震波之關係式如下所示：

) 1 1(

1 2 ²

1

2 +

+ +

= M_s

P P

γ

γ ( 2-1 )

1 ) 1 2 (

1

1

2 − +

= +

P M_s P

γ

γ ( 2-2 )

⎟⎟

⎟⎟

⎠

⎞

⎜⎜

⎜⎜

⎝

⎛

− + +

− + +

=

1 2 1 2

1 2 1 2

1 1 1

1 1

P P P P P

P T T

γ γ γ γ

( 2-3 )

1 2 1 2

1 2

1 1

) 1( 1 1

P P P P

− + +

− + +

= γ γ

γ γ ρ

ρ ( 2-4 )

上式下標定義如下：

1：表示震波前方之氣體狀態 2：表示震波後方之氣體狀態

第三章 數值方法

3-1 統御方程式

本研究使用二維納維爾-史脫克方程式( Navier-Stokes equation ) 為统御方程式，且假設為理想氣體狀態，流體為可壓縮且具黏性之流 體，其守恆方程式如下所示：

v v

v H

r F x H E r F x E t

Q r r r

r r r r

∂ + +∂

∂

= ∂

∂ + +∂

∂ +∂

∂

∂ ( 3-1 )

其中，Q^r為守恆變數定義為：

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

= e

v Q u

ρ ρ ρ

r

其中 ρ 為密度，u、v 分別為 x 與 r 方向速度分量，而 e 為單位 體積之能量，可由理想氣體方程式得到：

(

^{2 2}

2 1

1 u v

e P + +

= − ρ

γ

)

( 3-2 ) 而Er Fr

、 為非黏性通量，Hr 為非黏性通量之源項定義如下：

( )

⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎢

⎣

⎡

+

= +

v P e

uv P u u

E ρ

ρ ρ r 2

( )

⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎢

⎣

⎡

+

= +

v P e

P v uv u

F ₂

ρ ρ ρ r

( )

⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎢

⎣

⎡

+

=

v P e v uv v

H r

ρ ρ ρ ρ

2

r 1

而Er_v Fr_v

、 為黏性通量，Hr_v則為黏性通量之源項定義如下：

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

=

x xr xx

Ev

β τ τ 0

r

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

=

r rr xr

Fv

β τ τ 0 r

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢

⎣

⎡

⎟⎠

⎜ ⎞

⎝

⎛

∂

− ∂

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂

− ∂

−

⎟⎠

⎜ ⎞

⎝

⎛

∂

− ∂

⎟⎠

⎜ ⎞

⎝

− ⎛

−

⎟⎠

⎜ ⎞

⎝

⎛

∂

− ∂

=

r uv x r r

v r r r

v

r v r r r

v r v x r

H r

r rr xr

v

μ μ μ

β

μ μ τ

τ

μ τ

θθ

Re 3 2 Re

3 2 Re

3 2

Re 3 2 Re

3 2 Re 3 2 0

1

2 2

r

上式中 τ 為剪應力及 β 為熱通量。

⎟⎠

⎜ ⎞

⎝

⎛

∂

−∂

∂

= ∂

r v x u

xx 2

Re 3 2 μ

τ

⎟⎠

⎜ ⎞

⎝

⎛

∂ + ∂

∂

−∂

= r

v x u

rr 2

Re 3 2 μ

τ

⎟⎠

⎜ ⎞

⎝

⎛ +

∂ +∂

∂

− ∂

= r

v r v x

u ) 2

Re ( 3 2 μ

τ_θθ

⎟⎠

⎜ ⎞

⎝

⎛

∂ +∂

∂

= ∂

= x

v r u

rx

xr Re

τ μ τ

x T v M

u

e xr

xx

x ∂

∂

− − +

= ₂

) 1 Pr(

Re γ τ μ

τ β

r T v M

u

e rr

xr

r ∂

∂

− − +

= ₂

) 1 Pr(

Re γ τ μ

τ β

在方程式( 3-1 )式中之無因次化參數分別如下:

D x x

ˆ

= ˆ ， D r r

ˆ

= ˆ ，

Ue

D t t

ˆ ˆ

= ˆ ，

ρe

ρ ρ ˆ

= ˆ

， Te

T T ˆ

= ˆ Ue

u u ˆ

= ˆ ， Ue

v v ˆ

= ˆ ，

μe

μ μ ˆ

= ˆ

ˆ2

ˆ ˆ

e eU P P

= ρ ，

k C_p

= μ

， ₂ Pr ˆ ˆ

ˆ

e eU e e

= ρ ，

e e

eu D

μ

= ρ Re

上式參數中有上標符號“＾”代表有因次之物理量，而下標 e 則表

示震波管出口處之物理狀態，其中 μ 為無因次黏滯係數，係經由

Sutherland law 計算如下：

) 4 . 110 (

) 4 . 110 1

2 (

3

e e

T T

T T

+

= +

μ ( 3-3 )

3-2 時間積分

本 研 究 中 對 統 御 方 程 式( 3-1 ) 式 之 時 間 積 分 ， 使 用 四 階 Runge-Kutta scheme：

^{v v} H_v r F x H E r F x E t

Q r r r

r r r r

∂ + +∂

∂

= ∂

∂ + +∂

∂ +∂

∂

∂

( )

ⁿ

n tLQ

Q

Qr r r

Δ +

= 2

) 1

1 (

( )

( )¹

) 2 (

2 1 tLQ Q

Qr rⁿ r

Δ +

=

( )

( )²

) 3 (

2 1 tLQ Q

Qr rⁿ r

Δ +

=

( ) ( ) ( )

(

^{1 2 3}

) ( )

^{( )3}

1

6 2 1

3

1 Q Q Q Q tLQ

Qrⁿ rⁿ r r r r

Δ + + +

+

−

+ =

( )

^Q

dt L Q dr r

=

H r H

F x F E

r F E

x E Q

L ^{v v} _v

j i j i j

i j i

r r r

r r r

r r

r + −

∂ +∂

∂ +∂

⎥⎦

⎢ ⎤

⎣

⎡ −

−Δ

⎥⎦

⎢ ⎤

⎣

⎡ −

−Δ

= + − + −

2 , 1 2 , 1 2,

, 1 2 1

1 ) 1

( ( 3-4 )

3-3 空間差分

對於納維爾-史脫克方程式中對流項( convective terms )之空間差 分的近似值計算，本文是採用WENO scheme ( Weighted Essentially Non-oscillatory scheme )，源自 1987 年由 Engquist 等人[13]提出 ENO scheme 在數值不連續處容易出現失真，所以 Liu 等人[14]與 Jiang 及 Shu[15]在 1994 及 1995 年分別對於 ENO 的差分方式進行平滑化 ( Smoothness )，其方法為在對流項Er(Q)

及Fr(Q)

進行空間差分時，對 於計算差分時所需的格點皆分配一個小小的重量( weight )ω_k於其 上，以進行平滑化，並藉以由補差的方式使得精密度能夠達到( 2r-1 ) 階，本文計算時選則r = 3，所以在空間上的精度可以到達 5 階。

為了求得方程式( 3-4 )中的

j i

E , 2 +1

r 的高階近似，因此Jiang 及 Shu

[15]定義：

∫

₋^+ΔΔ

( )

= Δ ²

2

) 1 (

x x

x xh d

Q x

E^r ξ ζ

( ) [ ( ) ( ) ]

x x x x h

x h Q

E Δ

−Δ Δ −

= + 2 2

r

接 著 對 全 流 場 進 行 通 量 分 離( flux splitting )， 本 文 係 採 用 Lax-Friedrichs ( LF ) scheme 其Er(Q)之通量分離式如下：

) ( ) ( )

(Q E Q E Q

Er = r⁺ + r⁻ 1 r r

其中α =maxEr'(Q)

也 就 是 必 須 再 差 分 以 得 到 左 邊 邊 界及 右 邊 邊 界 的 數 值 通 量 ( numerical flux )。

) ( )

( )

( ,

2 , 1

2 , 1

2

1 Q E Q E Q

E

j i j

i j

i

− + +

+

+ = r + r

r ( 3-5 )

在計算 ( )

2,

1 Q

E

j i

+ +

r 之正負數值通量時，採用下式以得到其近似值：

∑

⁻

= + + − + +

+

+ = ¹ ⋅ ⋅

0 , 1, ,

2,

1 ^r ( ,..., )

k

j k i s j

r k i s r k s j k

i

E l E

l q

E^r ω

∑

⁻

= − + −+ +

−

+ = ¹ ⋅ ⋅

0 , 1, ,

2,

1 ^r ( ,..., )

k

j k i s j

r k i s r k s j k

i

E l E

l q

E^r ω

) ,...,

( _{1, 1,}

,s k s i r j s i r j

k⁺ =ω l ⋅E_{− +} l ⋅E_{+ −}

ω

) ,...,

( _{2, ,}

,s k s i r j s i r j

k⁻ =ω l ⋅E_{− +} l ⋅E₊

ω

(

^{k=0,1,...,r−1}

)

± +2 i 1

Er 可寫為：

所以

∑

=

± +

±

+ = ^m ⋅

s

j s i j

i

E r

E

1 ,

2 , 1

2 1

r r

上述的 及 為特徵矩陣。

而 為一個對於計算差分時所需的格點分配之一個小重量，其

定義如下：

ls r_s

± s

ωk ,

±

±

±

± ±

+

= +

2 1 0

, α α α

ω_ks α^k ( k = 0, 1, 2 )

(

0

)

²

0

1

+ +

= + ε IS α

(

1

)

²

1

6

+ +

= + ε IS α

( )

²

2 + = 3 α

對於上式的ε其存在的意義在於為了避免分母為零導致程式發 散，所以必須假設一極小量值，Jing 及 Shu[15]建議ε =1×10⁻⁶。

( ) (

2, 1, ,

)

²

2 , , 1 ,

2

0 4 3

4 2 1

12

13 _{+ +}

− +

− +

+

− +

−

+ = E_{i j} − E_{i j} +E_ij + E_{i j} − E_{i j} + E_ij

IS r r r r r r

( ) (

1, 1,

)

²

2 , 1 ,

, 1

1 4

2 1 12

13 ₊

+ + + −

+ + +−

+ = E_{i j} − E_{i j} +E_{i j} + E_{i j} −E_{i j}

IS r r r r r

( ) (

, 1, 2,

)

²

2 , 2 ,

1 ,

2 3 4 3

4 2 1

12

13 +

+ +

+ +

+ + +

+ +

+ = E_{i j} − E_{i j} +E_{i j} + E_ij − E_{i j} + E_{i j}

IS r r r r r r

(

0

)

²

0

3

−

−

= + ε IS α

(

1

)

²

1

6

−

−

= + ε IS α

(

2

)

²

2

1

−

−

= + ε IS α

( ) (

1, , 1,

)

²

2 , 1 ,

, 1

0 4 3

4 2 1

12

13 ₋

− +

−−

−+

−

−−

− = E_{i j} − E_{i j} +E_{i j} + E_{i j} − E_ij + E_{i j}

IS r r r r r r

( ) (

, ,

)

²

2 , 2 ,

1 ,

1 4

2 1 12

13 _{− − −}

− +

− +

− = E_{i j} − E_{i j} +E_{i j} + E_ij −E_{i j}

IS r r r r r

( ) (

1, 2, 3,

)

²

2 , 3 ,

2 ,

1

2 3 4 3

4 2 1

12

13 ₋

− +

− +

− +

− +

− +

− = E_i+ _j − E_{i j} +E_{i j} + E_{i j} − E_{i j} + E_{i j}

IS r r r r r r

所以( 3-4 )式可改寫為：

( ) ( )

( )

6 5 2

6

2 5

6

11 7

2

, 2 ,

1 ,

2

, 1 ,

, 1 1

, ,

1 ,

2 0 2,

1

++ ++

+ +

+ + +

+

− + +

+

− +

− + +

+

− + +

+ +

+ − +

= −

j i j i j i

j i j i j i j

i j i j i j

i

E E

E Q

E E

E Q E

E E

E Q

r r

r

r r

r r

r r r

( 3-6 )

( ) ( )

( )

6 5 2

6

2 5

6

11 7

2

, 2 ,

1 ,

2

, 1 ,

, 1 1

, ,

1 ,

2 0 2,

1

++ ++

+ +

+ + +

+

− + +

+

− +

− + +

+

− + +

+ +

+ − +

= −

j i j i j i

j i j i j i j

i j i j i j

i

E E

E Q

E E

E Q E

E E

E Q

r r

r

r r

r r

r r r

( 3-7 )

r

而得到Ev_i,j的值。至於Fr(Q)

之通量分離步驟與方法和Er(Q)

類似，因此 不再重覆敘述。

對於納維爾-史脫克方程式中黏滯項( viscous terms )之空間差分 的近似值計算，本文是採四階中央差分( forth-order centered difference scheme )[16]，其方程式定義如下：

( ) ( )

[

2 , 8 , 8 ( , ) ( 2 , )

]

12

1 u x h r u x h r u x h r u x h r h

dx

du = − − − + + − +

( ) ( )

[

, 2 8 , 8 ( , ) ( , 2 )

]

12

1 u x r h u x r h u x r h u x r h h

dr

du = − − − + + − +

) , ( 30 ) , ( 16 ) , 2 ( 12 [

1

2 2

2

r x u r h x u r h x h u

dx u

d = − − + − −

)]

, 2 ( ) , (

16u x+h r −u x+ h r +

) , ( 30 ) , ( 16 ) 2 , ( 12 [

1

2 2

2

r x u h r x u h r x h u

dr u

d = − − + − −

)]

2 , ( ) , (

16u x r+h −u x r+ h +

) , ( 10 ) , ( 10 ) , ( 10 24 [

1

2 2

h r h x u h r h x u h r h x h u

dxdr u

d = − − − − + − + −

+10u(x+h,r+h)−u(x−2h,r−h)+u(x−2h,r+h)

+u(x+2h,r−h)−u(x+2h,r+h)−u(x−h,r−2h) )]

2 , ( ) 2 , ( ) 2 ,

(x h r h u x h r h u x h r h

u − + + + − + + +

+

)]

, 2 ( ) , ( 8 ) , ( 8 ) , 2 ( 12 [

1 x h r x h j x h j x h j

h dx

dμ = μ − − μ − + μ + −μ +

)]

, 2 ( ) , ( 8 ) , ( 8 ) , 2 ( 12 [

1 T x h r T x h j T x h j T x h j h

dx

dT = − − − + + − +

因此統御方程式中 及

dr F dr_v dx

E dr_v

可寫成下式：

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢

⎣

⎡

+ − + −

+ +

+

+ +

+

+

− +

+

−

=

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎣

⎡

=

2 2

2 2

2

2 2

2 2 2 2

2

) 1 Pr(

) 1 Pr(

1 )

( ) (

)]

3( 2 2

[ )]

3( 2 2

[ 0

Re 1 0

dx T d M dx

dT dx d M dx

vd dx

dv dx ud dx

du

dxdr u d dx

v d dr

du dx dv dx d

dxdr v d dx

u d dx

u d dr

dv dx du dx

du dx d

dx d

dx d

dx d

dx E d

e e

xr xr

xx xx

x xr xx

v

γ μ μ

γ τ τ

τ τ μ μ μ μ

β τ τ r

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎢

⎣

⎡

+ − + −

+ +

+

+

− +

+

−

+ +

+

=

⎥⎥

⎥⎥

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎢⎢

⎢⎢

⎣

⎡

=

2 2

2 2

2 2 2

2 2 2 2 2

) 1 Pr(

) 1 Pr(

1

)]

3( 2 2

[ )]

3( 2 2 [

) (

) (

0

Re 1 0

dr T d M dr

dT dr d M dr

vd dr

dv dr ud dr

du

dr v d dxdr

u d dr

v d dr

dv dx du dr

dv dr d

dr u d dxdr

v d dr

du dx dv dr d

dr d

dr d

dr d

dr F d

e e

rr rr

xr xr

r rr xr

v

γ μ μ

γ τ τ

τ τ μ μ μ μ

β τ τ r

3-4 起始條件與邊界條件

二維軸對稱脈衝噴流衝擊平板模擬區域示意圖如圖2-1 所示。

3-4.1 起始條件：

本研究之起始條件設定是參考文獻[12]之實驗條件，所有的氣體 性質皆以震波管出口之震波性質為無因次化依據。因此，震波管出口 之無因次震波壓力、溫度、密度及速度分別為0.6867、1、1 及 1，而 大氣之無因次壓力、溫度、密度及速度分別為0.1373、0.6314、0.3168 及0。

3-4.2 邊界條件：

在入口處之邊界：其值為起始條件所設定之狀態，並保持固定 值。

壁面邊界：採用無滑動邊界 ( non-slip boundry condition )，其壁 面上速度為零，平板壁面溫度除了考慮熱傳效應時變化其溫度外，其 餘皆假設為絕熱狀態。

自由邊界：在此是使用黎曼不變量邊界條件( Riemann invariant boundary condition )，此邊界法是假設一二維流場在邊界時可視為近 似一維自由流，因此根據特徵線理論，沿特徵線R⁺及R⁻各有一個黎 曼不變量，其可表示為

1 2 + −

=

+ c

V

R _n

1 2

− −

=

−

r V c

R _n

其中 C 為音速，V_n =V⋅nr，亦即表示垂直於流場邊界的速度分量，

以向外為正，R⁺及R⁻相對應的特徵速度為：

c V_n −

1 = λ

c V_n +

2 = λ

此時已有兩個方程式，但是為了求解四個變數值，所以還需要 另外兩個方程式，此時可以選擇熵( entropy )的關係式 _⎟⎟

⎠

⎜⎜ ⎞

⎝

= ⎛ ρ

S ln p 及切線

速度( V_t )為另外兩個關係式，來求解四個特徵值。

當超音速出流時，λ₁及λ₂皆大於 ，所以所有邊界值皆由流場內 部插分而得，當次音速出流時，

0

1<0

λ ，λ₂ >0，此時R⁻由自由流條件 而來，R⁺、 及熵( S )則由流場內部給定，如此將可求得各變數之邊 界值。

Vt

3-5 局部紐賽爾數及平均紐賽爾數

為了瞭解不同加熱溫度對壁面熱傳效應之影響，本文使用局部 紐賽爾數( Local Nusselt number )以及平均紐賽爾數( Average Nusselt number )來觀測震波對流場熱傳效應之影響，其定義如下：

L

w

x

T T T x D

Nu

∂

∂

−

= −

) ) (

(

1

( 3-8 )

∫

= ^L

avg

Nu x dx

Nu L

0 ( )

1 ( 3-9 )

第四章 結果與討論

本研究主要在探討脈衝噴流從震波管流出衝擊加熱平板之後，

其震波與平板之交互作用、反射震波對流場結構及平板熱傳效應的影 響。惟在未進行上述研究之前，吾人先進行流場無平板時，脈衝噴流 在不同的出口壓力下從震波管流出後之流場結構模擬，計算震波馬赫 圓盤直徑( DM )及位置( XM )隨時間之變化情形，並將計算結果與實驗 數據[12]和其他研究者之模擬結果[11]相比較，以驗證程式的正確性。

4-1 無平板之脈衝噴流結構

吾人模擬之脈衝噴流係以 Ishii 等人[12]之實驗條件為參考依 據，震波管直徑為D = 2cm，震波管壁厚為 1cm，震波管出口與大氣 之壓力比為PE/P1 = 2.2、3.61 及 5.0，震波馬赫數為 Ms = 1.0、1.0 及 1.02，大氣溫度為 T1 = 289K。由於 Ishii 等人[12]在其數值模擬研究中 提到，計算區域若採用 r = 4R 及 x = 15R 便足以模擬脈衝噴流流出震 波管後至t=1000μs 的噴流發展( Jet evolution )，同時若無因次化之網 格 大 小 採 用 Δx = Δy = 0.025 ， 則 可 準 確 模 擬 噴 流 邊 界 ( Jet Boundary )、膨脹波( Expansion Wave )、震波( Shock Wave )及渦卷結 構( Vortical Structure )等現象。因此，對於無平板之脈衝噴流結構模 擬，本研究採用之計算區域為 r = 4R 及 x = 15R，網格大小為 Δx = Δy

= 0.025。

Ishii 等人[12]在實驗過程中觀察到脈衝噴流發展大致可分為三 個階段，第一階段為噴流流出震波管後在出口端產生震波繞射( Shock Diffraction )，第二階段為形成一非穩定之馬赫圓盤，第三階段為形成 第一個震波囊結構。而在第二階段時之噴流結構示意圖如圖 4-1 所 示。圖中DM及XM表示馬赫圓盤的直徑及位置，而( XV, YV )代表第 一個渦卷(Vortex)的位置，同時圓桶狀震波( Barrel Shock )、噴流邊界 ( Jet Boundary )、第一個渦卷( 1^st Vortex )、平滑線( Slip Line )、反射 震波( Reflected Shock )、第二個渦卷( 2^nd Vortex )及渦卷誘發震波 ( Vortex-induced Shock )都在這個階段形成。圖 4-2 為 t = 250μs 時，三 種不同震波管壓力比為PE/P1 = 2.2、PE/P1 = 3.61 及 PE/P1 = 5 時的密度 與壓力等位線圖。從圖4-2(a)-(c)中可觀察到馬赫圓盤位置約在 x/R = 2.5、x/R = 3.8 及 x/R = 4.2 處，而渦卷中心約在 x/R = 3.8、x/R = 4.2 及x/R = 5.1 的位置，且有噴流邊界、膨脹波、震波及渦卷等等的流 場結構皆已成形。

為了進一步驗證計算程式的正確性，吾人將三種不同壓力比下所 計算之馬赫圓盤直徑及位置隨時間變化的情形分別為 Ishii 等人[12]

之實驗與計算結果，以及和Kim 及 Park[11]的計算結果相互比較，如 圖4-3 所示。由圖中可見，當出口壓力比為 P_E/P1 = 2.2 時，馬赫圓盤

直徑最大値( DM ≒ 1cm )大約發生在 t = 95μs 附近，而其軸向位置隨 時間而往下游移動，當t = 200μs 時馬赫圓盤軸向位置幾乎保持在 XM

= 2.4cm 的位置處。當出口壓力比升高至 3.61 及 5 時，馬赫圓盤直徑 隨著時間的增加而增大，其最大馬赫圓盤直徑發生在約t = 100μs 時，

然後隨著時間的增加而逐漸縮小，當t = 300μs 時，馬赫圓盤直徑幾 乎不再變化。而馬赫圓盤的位置則隨時間的增加而往下游移動，當t = 250μs 時馬赫圓盤軸向位置幾乎保持在固定位置，不再往下游移動。

由圖4-3 的計算結果可見，本研究所發展的程式不但可準確的預測馬 赫圓盤直徑與軸向位置的不穩定現象，其結果與實驗及其他數值模擬 結果相比較亦非常吻合，足見本程式在邊界條件的設定及程式的撰寫 正確性。

4-2 格點測試

格點數的選擇對流場計算的結果影響很大，格點數過少時，會 因為格點間距離過大，無法精確的分析流場間的物理量，造成計算上 的誤差，格點數過多時，雖然可以較精確描述流場內的物理量，但會 增加計算時間，基於精確度和時效性的考量，在進行數値模擬時須做 格點測試。

吾人為了瞭解計算脈衝噴流衝擊平板時所需的網格大小，選擇

震波管至平板距離x/D = 2.5，震波管出口壓力比 PE/P1 = 5 進行格點 測試，並觀測在r/D = 0 處平板表面之壓力振盪情形，如圖 4-4 所示。

所使用的四種網格大小分別為：Grid1 : Δx = Δr = 0.025D；Grid2 : Δx = Δr = 0.00125D；Grid3 : Δx = Δr = 0.001D；Grid4 : Δx = Δr = 0.0075D，

其中 D 表示震波管的直徑，Grid1、Grid2、Grid3 和 Grid4 四種格點 數的最大壓力，分別為 P = 4.439P1 、P = 6.170P1 、P = 7.826P1 及 P = 7.423P1，其中P1 = Pa，Grid1 和 Grid2 的壓力相差較大，而 Grid2、

Grid3 和 Grid4 壓力相差較小，Grid2 和 Grid3 的壓力誤差為 21%，Grid3 和Grid4 的壓力誤差為 5.2%，為了減少其計算時間採用 Grid3 為本研 究之計算格點，以x/D = 2.5 為例，當 Δx = Δr = 0.01D 時格點數為 300(x)×250®。

4-3 絕熱平板之脈衝噴流結構

吾 人 為 了 探 討 脈 衝 噴 流 衝 擊 平 板 後之 震 波 流 場 結 構 ， 以 圖 4-5(a)-(f)中的流場現象來做名稱上的定義，圖 4-5(a)-(f)為震波管至平 板距離 x/D = 2.5、無因次時間 t=4.7~14 且噴嘴壓力比 P_E/P₁ = 5 時，

震波衝擊絕熱平板後之流場結構數值紋影圖，其中震波現象包括有起 始震波( Primary Shock Wave )、第一個反射震波( Reflected Shock 1 )、

噴流邊界( Jet Boundary )、圓桶狀震波( Barrel Shock )、馬赫圓盤

( Mach Disk )、平滑線( Slip Line )、反射震波( Reflected Shock )、第 一個渦旋( Vortex )、第二個反射震波( Reflected Shock 2 )、渦卷誘發 震波( Vortex-induced Shock )、第二個渦旋( 2^nd Vortex )和三匯合點 ( T )，而三匯合點包含反射震波、馬赫圓盤和平滑線，上述反射震波 及反射震波與渦卷交互作用現象在 Ishii 等人[12]的實驗和數值模擬 中無法觀察到。

在圖 4-5(a)-(c)中，震波管出口端產生震波繞射現象，其起始震 波會在特定位置，即馬赫數Ms = 1 時產生馬赫圓盤，而後在渦卷處 形成一個震波囊結構，形成震波囊結構後，起始震波持續往前直至衝 擊平板後反射，其第一個反射震波( R1 )會因為反射後跟震波結構相互 影響而改變其波形，會往回行進直至撞擊馬赫圓盤後再度往平板反 射，此第一個反射震波會在平板和馬赫圓盤之間進行往復的反射運 動，直至能量消散為止，而衝擊平板後不會因為震波結構而改變其波 形的反射震波，稱為第二個反射震波( R2 )，在圖 4-5(d)-(f)中，第一 個反射震波會和震波囊結構結合，跟圖4-2 相比，震波不再是單方向 的行進，而是在一個區域內往復的傳遞能量，在傳遞的過程中，震波 結構會受到震波破裂、結合及分離等等現象的影響而改變其形狀大 小，在圖4-5(d)中，第一個反射震波和馬赫圓盤後的反射震波結合而

所形成的位置可用來觀察馬赫圓盤產生的位置及大小，在圖 4-5(f) 中，震波囊衝擊平板後反射，使原先的震波結構受到反射震波的交互 作用影響而改變形狀，包括噴流邊界和圓桶狀震波都受到反射震波的 撞擊而變形，而馬赫圓盤也會因為受到衝擊而縮小，並在此時產生馬 赫圓盤最小值。

4-4 平板熱傳效應之影響

圖 4-6(a)-(e)為在相同噴嘴出口壓力比( P_E/P₁ = 5 )下，不同震波 管至平板距離( x/D = 2.0 ~ 4.0 )，相同無因次時間( t = 10 )，相同平板 溫度( T_W = 8T₁ )時之密度等位線圖，由圖中可見，起始震波已衝擊加 熱平板後反射且在震波管至平板距離x/D = 2.0 及 2.5 的情況下，震波 囊結構已有破裂現象，而在震波管至平板距離為x/D = 3.0、3.5 及 4.0 的情況下，震波囊結構仍然完整。圖 4-7 為在相同噴嘴出口壓力比 ( P_E/P₁ = 5 )，震波管至平板距離為 x/D = 2.5，相同無因次時間( t = 10 )，不同平板溫度( T_W = 2T₁ ~ 10T₁ )時之溫度等位線圖，圖 4-7(a)-(e) 之平板溫度分別為T_W = 2T₁、4T₁、6T₁、8T₁及10T₁，將圖4-7 和圖 4-5(e)的震波結構相比較，會發現其平板溫度對震波結構的影響不 大，為了瞭解加熱平板上之熱傳效應，選取接近平板位置( x/D = 3.4 ~ 3.5 )的溫度等位線圖來觀察其近壁面之溫度分佈，如圖 4-8 所示，由

圖中可見，當平板溫度逐漸增高時其等溫線也隨著增密，較密的等溫 線其溫度梯度也較高，因此其熱傳效果也會較好。

4-5 相同震波管至平板 距 離及不同平板溫度對紐塞爾數之 影響

吾人從( 3-8 )式中瞭解紐塞爾數的定義，當平板壁面與大氣的溫 度差越大則紐塞爾數越小，反之亦然，而紐塞爾數增加時會使熱對流 效率提高，此亦表示其熱傳效率好，為了瞭解此脈衝噴流衝擊加熱平 板( 即不同平板溫度 )後對平板熱傳效應之影響，可探討平板上紐塞 爾數的變化情形，由於所有溫度參數皆已無因次化，因此影響紐塞爾 數大小的主要參數為溫度梯度。

圖 4-9(a)-(e)為在相同噴嘴壓力比( PE/P1 = 5 )下，相同震波管至 平板距離，不同平板溫度( TW = 2T1 ~ 10T1 )下之平板局部紐塞爾數分 佈圖，在圖 4-9(a)中紐塞爾數的分佈情形為當平板溫度增加時其紐塞 爾數越小，惟平板溫度 TW = 2T1時紐塞爾數分佈不同，其紐塞爾數約 從噴流中心位置r/D = 0 時逐漸增加，至 r/D = 0.2 後紐塞爾數才維持 一平緩數值，直至r/D = 0.8 時紐塞爾數才又快速增加，而在平板溫度 TW = 4T1 ~ 10T1時則紐塞爾數呈現平緩增加的趨勢。在圖4-9(b)中紐 塞 爾 數 的 分 佈 情 形 為 當 平 板 溫 度 增 加 時 其 紐 塞 爾 數 越 大 ， 在 圖

4-9(c)、(d)及(e)中，平板溫度 TW = 4T1 ~ 10T1的紐塞爾數已不再有大 幅度的變化，由於平板溫度TW = 2T1的紐塞爾數變化較平板溫度TW = 4T1 ~ 10T1不同，因此特別探討平板溫度 TW = 2T1時平板上之紐塞爾 數變化。圖4-10(a)-(e)為在相同噴嘴壓力比( PE/P1 = 5 )下，不同震波 管至平板距離( x/D = 2.0 ~ 4.0 )，相同平板溫度( TW = 2T1 )之平板局部 紐塞爾數分佈圖，其紐塞爾數都會先增加後減少而後再增加，在震波 管至平板距離x/D ＜ 3.0 時，其紐塞爾數約在 r/D = 0.5 後增加較劇 烈，而在震波管至平板距離x/D ≧ 3.0 時，其紐塞爾數約在 r/D = 0.5 後增加較平緩。

平板熱傳效應之影響也可由平均紐塞爾數得知，而圖 4-11 為在 相同噴嘴壓力比( PE/P1 = 5 )下，改變震波管至平板距離及平板溫度 ( TW = 2T1 ~ 10T1 )之平均紐塞爾數分佈圖，由圖中可見，在固定震波 管至平板的距離之下，平板壁面溫度 TW = 2T1時其平均紐塞爾數皆比 其它壁面溫度要低，當x/D ≧ 3 之後，壁面溫度為 TW = 4T1 ~ 10T1

之平均紐塞爾數值非常相近，此結果顯示，當震波管至平板距離x/D

≧ 3 之後，平板溫度若大於 4T1，則其壁面溫度再增高也對平均紐塞 爾數無影響。在固定壁面溫度的情況下，平均紐塞爾數皆隨著震波管 至平板距離的增加而增加，當壁面溫度TW = 2T1時，其平均紐塞爾數 增加的速度比其它壁面溫度要快。

第五章 結論

本文利用 WENO 數值方法成功模擬出在二維軸對稱流場內之震 波結構與物理現象，主要在探討脈衝噴流衝擊加熱平板之後，其震波 與平板之交互作用、反射震波對流場結構及平板熱傳效應的影響，為 了瞭解上述之震波流場結構，吾人對流場進行無平板、絕熱平板及加 熱平板三種脈衝噴流之數值模擬。

對於無平板之脈衝噴流流場結構之發展大致上可分為三個階 段，第一階段為噴流流出震波管後在出口端產生震波繞射，第二階段 為形成一非穩定之馬赫圓盤，第三階段為形成第一個震波囊結構，上 述震波現象在Ishii 等人[12]的實驗和數值模擬中亦可觀察到，對於絕 熱平板之脈衝噴流結構之發展，重點為震波跟平板之交互作用和平板 上熱傳效應之影響。震波由震波管流出後，其起始震波會在特定位 置，即馬赫數M_s =1時產生馬赫圓盤，起始震波會持續往前直至衝擊 平板後反射，第一個反射震波會在平板和馬赫圓盤之間進行往復的反 射運動，直至能量消散為止，第一個反射震波和馬赫圓盤後的反射震 波互相作用而有小渦卷產生，此小渦卷稱為第二個渦卷，對於加熱平 板之脈衝噴流流場結構之模擬，會發現其震波結構跟絕熱平板之震波 結構差異性不大，不同點在於當平板溫度逐漸增高時其等溫線也隨著

為了瞭解此脈衝噴流衝擊加熱平板( 即不同平板溫度 )後對平 板熱傳效應之影響，可探討平板上紐塞爾數的變化情形，由於所有溫 度參數皆已無因次化，因此影響紐塞爾數大小的主要參數為溫度梯 度。在圖 4-9(a)中紐塞爾數的分佈情形為當平板溫度增加時其紐塞爾 數越小，而在平板溫度TW = 4T1 ~ 10T1時則紐塞爾數呈現平緩增加的 趨勢，在圖 4-9(b)中紐塞爾數的分佈情形為當平板溫度增加時其紐塞 爾數越大，圖4-10(a)-(e)，其紐塞爾數都會先增加後減少而後再增加，

在震波管至平板距離x/D ＜ 3.0 時，其紐塞爾數約在 r/D = 0.5 後增 加較劇烈，而在震波管至平板距離x/D ≧ 3.0 時，其紐塞爾數約在 r/D

= 0.5 後增加較平緩。當震波管至平板距離 x/D ≧ 3 之後，平板溫度 若大於4T1，則其壁面溫度再增高也對平均紐塞爾數無影響。在固定 壁面溫度的情況下，平均紐塞爾數皆隨著震波管至平板距離的增加而 增加，當壁面溫度TW = 2T1時，其平均紐塞爾數增加的速度比其它壁 面溫度要快。