• 沒有找到結果。

中 華 大 學

N/A
N/A
Protected

Academic year: 2022

Share "中 華 大 學"

Copied!
74
0
0

加載中.... (立即查看全文)

全文

(1)

中 華 大 學 碩 士 論 文

採用大尺度渦流模擬法於正方形孔穴內混合 對流之暫態與紊流場分析

Large eddy simulation of mixed convection in transient laminar and turbulent flows in a

square cavity

系 所 別:機械工程碩士班

學號姓名: M09908013 沈富第 指導教授: 楊 一 龍 博士

中華民國 101 年 九 月

(2)

中文摘要

本研究採用大尺度渦流於正方形孔穴內混合對流之暫態與紊流場分析,次格點模式為 Smagorinsky model,統御方程式採用有限體積法求解,計算之結果其雷諾數分別介 於 400 至 10,000 以及理查森數介於 0.1 至 0.44 之間,層流暫態流場之發展過程與文 獻中實驗結果相當吻合,至於紊流熱傳之發展與最終之變動量也與前人結果很接近,

至於非穩定熱分層相較穩定熱分層會有 21%以上之熱傳增益,最後在頻率響應分析上 可發現邊界層內流速與壓力之反應一致皆有基頻之倍頻現像,但在溫度傳遞上則偏向 基頻以內之低頻為主,因此溫度與速度之跳動呈現不一致之行為。

關鍵詞:大尺度渦流、混合對流、過渡流場、驅動腔

(3)

英文摘要

This work presents large eddy simulation of mixed convection in transient,

two-dimensional laminar and turbulent flows in a square cavity.The Smagorinsky model is employed for the sub-grid treatment. The simulations are based on the finite volume solution of the conservation equations, and are presented for Reynolds and Richardson numbers ranging from 400 to 10,000 and from 0.1 to 0.44, respectively. The agreement in solutions for the transient laminar flow is very well with those in the literature. Results for the development of turbulent flow and final turbulent fluctuation are compared well the previous work. In general, the results show a 21% increase of heat transfer from unstable stratification to stable one. The power spectral density of velocity and pressure moves to high frequency regime with significant resonant by the base frequency. The power spectral density of temperature dominates within the base frequency. The results indicate that the temperature and the velocity components fluctuate separately.

Keywords:LES,、Mixed convection, Transient flow, Driven cavity

(4)

誌謝

本研究得以順利完成,首先誠摯的感謝指導教授 楊一龍

本論文承蒙

博士,由於老師的細心 的提攜與指導,讓我在碩士期間了解的人生理念與實務領域的深奧,不時的討論並指 點我正確的方向,使我在這些期間獲益匪淺。

傅武雄博士、鄭藏勝

在我求學的日子裡,感謝碩士班的經勝、冠儒等學長,謝謝你們過去所做之研究 成果,提供了我學習基礎。也感謝同學

博士,對本文提供了許多寶貴的意見,因為有各 位委員的建議與指正,使得本論文更加完善。

子儀、子慶與碩翰以及學弟冠霖與佑達

最後,我要感謝我的雙親以及家人,感謝你們在我求學的過程中,默默的支持與 鼓勵,才能讓我無後顧之憂的完成學業,謹以本文獻給你們。

在研究 過程的協助與勉勵,並且陪我走過在實驗室裡共同的生活點滴。

(5)

目錄

中文摘要 ... i

Abstract ... ii

致謝 ... iii

目錄 ... iv

圖目錄 ... vi

符號定義 ... vii

第一章 緒論 ... 1

1.1 前言 ... 1

1.2 文獻回顧 ... 1

1.3 採用方法 ... 5

1.4 文章安排 ... 5

第二章 物理問題描述 ... 7

2.1 混合紊流問題與邊界條件 ... 7

2.2 大尺度窩流數值模擬 ... 7

2.3 基本假設 ... 8

2.4 無因次參數 ... 9

第三章 數值方法 ... 11

3.1 Navier-Stokes 方程式 ... 11

3.2 時間積分 ... 12

3.3 流通量計算 ... 12

3.4 穩定性分析 ... 13

3.5 人工黏滯 ... 17

3.6 固體壁面 ... 18

(6)

4.1 不同理查森數之影響 ... 19

4.2 暫態流場下直接數值分析與大尺度窩流之比較 ... 21

4.3 穩定熱分層與非穩定熱分層之影響 ... 22

第五章 結論 ... 26

5.1 結論 ... 26

5.2 未來展望 ... 27

參考文獻 ... 28

(7)

圖目錄

圖 2-1 邊界條件 ... 31

圖 4-1、不同位置下瞬時速度與溫度變化過程(Re=400、Ri=0.1) ... 32

圖 4-2、不同位置下瞬時速度與溫度變化過程(Re=400、Ri=1) ... 33

圖 4-3、不同位置下瞬時速度與溫度變化過程(Re=400、Ri=10) ... 34

圖 4-4 頂面與底面平均紐賽數之變化過程 ... 35

圖 4-5 不同時間下之流線分布圖(Re=400、Ri=0.1) ... 36

圖 4-6 不同時間下之流線分布圖(Re=400、Ri=1) ... 37

圖 4-7 不同時間下之流線分布圖(Re=400、Ri=10) ... 38

圖 4-8 不同時間下之溫度分布圖(Re=400、Ri=0.1) ... 39

圖 4-9 不同時間下之溫度分布圖(Re=400、Ri=1) ... 40

圖 4-10 不同時間下之溫度分布圖(Re=400、Ri=10) ... 41

圖 4-11 不同時間下之瞬時水平速度變化過程(Re=4000、Ri=0.44) ... 42

圖 4-12 不同時間下之瞬時溫度變化過程(Re=4000、Ri=0.44) ... 43

圖 4-13 不同時間下之流線分布圖(Re=4000、Ri=0.44) ... 44

圖 4-14 不同時間下之溫度分布圖(Re=4000、Ri=0.44) ... 45

圖 4-15 不同時間下之溫度分布圖,穩定熱分層(Re=10000、Ri=0.1) ... 47

圖 4-16 不同時間下之溫度分布圖,非穩定熱分層(Re=10000、Ri=0.1) ... 49

圖 4-17 時間下之溫度分布圖 ... 50

圖 4-18 中央與半腰處時間平均之溫度分布 ... 50

圖 4-19 穩定熱分層之底部時間平均紐賽爾數分布 ... 51

圖 4-20 非穩定熱分層之底部時間平均紐賽爾數分布 ... 51

圖 4-21 時間平均之流線分布圖 ... 52

圖 4-22 中央與半腰處時間平均之速度分布圖 ... 52

(8)

圖 4-24 四個邊界層 5mm 處之瞬時水平速度變化過程 ... 54

圖 4-25 四個邊界層 5mm 處之瞬時垂直速度變化過程 ... 55

圖 4-26 四個邊界層 5mm 處之瞬時溫度變化過程 ... 56

圖 4-27 四個邊界層 5mm 處之瞬時壓力變化過程 ... 57

圖 4-28 四個邊界層 5mm 處瞬時水平速度之頻率響應圖 ... 58

圖 4-29 四個邊界層 5mm 處瞬時垂直速度之頻率響應圖 ... 59

圖 4-30 四個邊界層 5mm 處瞬時溫度之頻率響應圖 ... 60

圖 4-31 四個邊界層 5mm 處瞬時壓力.之頻率響應圖 ... 61

圖 4-32 中央與半腰處 Utop V ' 10 , Utop U ' 10 分布圖 ... 62

圖 4-33中央與半腰處500 2' '

U

top

V U

分布圖 ... 62

圖 4-34 中央與半腰處

T T

∆ ' 分布圖 ... 63

圖 4-35 中央與半腰處

T U T U

top∆ ' ' 分布圖 ... 63

圖 4-36 中央與半腰處

T

U

V

T

U

top2∆ ' ' ' 分布圖 ... 64

(9)

符號定義

c : 局部音速

c

c

F

E ,

: 流速流通量

a

a

F

E ,

: 壓力流通量

v

v

F

E ,

: 擴散流通量

H

i : 內插函數 H : 總焓

k : 局部熱傳導係數

K0 : 入口熱傳導係數 LH : 加熱長度

Cp : 壓力係數, 2 0

0

2

1 U

P C

p

P

ρ

= −

Cf : 摩擦係數, 2

2

0

1 U

C

f w

ρ

= τ

M0 : 入口馬赫數,

RT

0

M U

= γ

Nu : 平均紐塞數,

( ) Nu ds

x x Nu

1

x 0

=

Nu : 局部紐塞數,

x T Tk Nu hk

= ∆

0

N : 耦合之點數 P : 壓力

Pr : 普朗特數, vα

ρ

: 密度

Re : 雷諾數,

H v

U

(10)

T : 溫度

T : 頂部與底部壁面之溫差

u

: 為

x

方向之速度

v

: 為

y

方向之速度

µ

: 黏滯係數

n : 壁面之法線方向

τ

ij : 流體剪應力

Γ

: 預置矩陣 ∆ : 時間步階

t

(11)

第一章 緒論

1.1 前言

傳統正方形孔穴內混合對流穩態流動之分析已經發展相當完備,本研究則希望記 錄驅動腔內從靜止速度開始上方驅動至流場接近穩定狀態下之暫態過程,以了解隨時 間過程中流場發展與熱傳現象。透過不同理查森數之設定,來探討驅動腔內浮力效應 對流場中溫度、速度與熱傳之影響。並針對紊流場之發展,使用大尺度渦流分析穩定 熱分層與非穩定熱分層之暫態過程,觀察體內渦流發展的複雜現象,以協助紊流模型 之建立。

1.2 文獻回顧

Prasad 與 Koseff[1]於 1998 年進行了一系列的實驗在頂部方型驅動腔內,腔體 深度與寬度為 150mm,雷諾數為 Re=3200 到 Re=10000 之範圍,翼展長寬比為 0.25:1 至 1:1 之流速量測。流場可視化使用聚苯乙烯珠與二維多普勒測速儀以了解腔體內動 量傳輸之過程。文獻中,觀察了在不同雷諾數與翼展長寬比下,沿著水平與垂直中心 線上的平均剖面速度與均方根剖面速度、剖面速度跳動。

Desai 與 Vafai [2]於 1995 年利用一圓筒帶動數形成對外環凹槽產生一個開放式 環型腔內流場,並控制外環每一凹槽的熱通量,可針對浮力效應與剪切流場進行熱傳 實驗與數值模擬之研究。由於凹槽很淺,外環半徑較大,因此可用於近似一系列水平 孔穴之混合對流實驗,由於每個孔穴之局部表面溫度皆有測量,可以用於對流之計 算,實驗中歸納出紐塞數與瑞利數之關聯。最後也利用數值模擬,以有限元素法進行 三維時間平均(雷諾)穩態流體運動和傳熱方程來解釋浮力引起的流動的基本物理特 性,至於平均紐塞數之比較上,實驗與數值模擬兩者結果差異不大。

(12)

浮力效應下紊流熱通量之模擬沿用流場之代數紊流模型,研究中採用一個混合時間尺 度,可針對流速項與紊流能量項做各別之時間尺度設定,以免除固定紊流普朗特常數 之設定。因此浮力效應,平均速度、紊流動能以及消散率和溫度變化以及其消散率皆 能以顯式方程式表示。沿續先前之近壁面的能量方程式處理方法,能透過近壁面的溫 度變化與消散率以及速度分布即可獲得溫度分布。擴充此一模型至浮力效應下之紊流 流場分析。模擬中進行底部加熱橫向完全發展流通道之分析以及垂直加熱管之紊流流 動,模擬結果與實驗數據結果一致。由於這些模擬結果,能夠了解在不可壓縮紊流場 下具浮力與無浮力效應下的基本物理現象。

Versteegh. 與 Nieuwstadt[4]於 1999 年以直接數值模擬兩個不同垂直加熱牆其 瑞利數介於 5.4x105~5.0×106之自然對流流場,其模擬結果與先前 Dafa'Alla 與 Betts 之實驗結果相比較,有很好的一致性。模擬進行了紊流的平均溫度與平均剖面 速度的統計,並與 George 與 Capp 所提出之自然對流垂直加熱表面紊流邊界層所歸 納之速度與溫度邊界層分布進行比對,其中在近牆的內層和外層的中心區域,直接數 值模擬結果與George 與 Capp 在平均剖面溫度分布一致,但在內層的速度上,平均 剖面速度的統計與 George 與 Capp 結果有些差異。至於雷諾應力、溫度變異與溫度 熱通量之比較上,也顯示速度變異量之尺度差異較大。

Peng 與 Davidson[5]於 2000 年在不同兩個加熱牆且有浮力效應流場的孔穴內 進行大渦流模擬數值研究。流體特性為穩定的熱分層與低紊流強度。流場模擬的時間 平均結果與實驗相當吻合,但在紊流統計上則模擬結果與實驗有些許的誤差。並在左 右熱傳面附近,所產生之高速渦捲,流經頂面與底面時並無法與外部結構融合,因此 雖然在看似均溫之區域其垂直方向之網格需增加以處理此一特殊渦流結構。

Ampofo, 與 Karayiannis[6] 於 2000 年對低紊流強度下自然對流空氣孔穴內進 行實驗研究。此孔穴高為 0.75M,寬為 0.75M 和深為 1.5M 的二維流場。孔穴內左右加 熱面溫度分別為 50 與 10 度,其瑞利數為 1.58x109,在實驗中於孔穴內不同位置同時 記錄了局部的溫度與速度分布,以及平均值和跳動量。也提供了局部的紐塞爾數,壁

(13)

面上的剪切應力以及紊流動能,擴散率和溫度變化,由於實驗數值有很高的精度,可 提供計算流體力學進行直接數值模擬之比對。

Tian 與 Karayiannis[7]於 2000 年對沿續先前低紊流強度下自然對流空氣孔穴 內之實驗結果。瑞利數為 1.58×105紀錄量包括:溫度的跳動量與速度 U、V 的跳動量 以及雷諾應力,流體流動受壁面剪切流屬於非等向性之擾動,在低紊流區的基頻介於 0.1-0.2 赫茲。其中流速之功率頻譜會往高頻方向延伸而溫度則屬低頻現象,且邊界 層內溫度與速度波動分布並非高斯分布,分析結果顯示,溫度和速度波動是不一致的。

Aounallah 等人[8]於 2007 年採用不同紊流模式針對下左右不同溫差之封閉孔穴 內之自然對流現象進行數值模擬。並針對不同的傾斜角度與波形振幅,利用各種紊流 模式進行熱傳比較。其中低雷諾數紊流模型(K-ω、K-ε、K-ω-SST)與較疏網格 DNS 模 擬結果與 Ampofo 與 Karayiannis 的實驗數據做比較,其結果顯示 K- -SST 模型與 實驗數據最為吻合。雖然 DNS 流場與溫度場的模擬與實驗平均結果也相當吻合,但紊 流的統計上差異較大,實驗整體熱通量比模擬的熱通量要高。此外,透過不同振幅的 波動壁,會改變局部的熱傳趨勢,在高瑞利數下,採用波形壁之熱傳量比正方形孔穴 熱傳量要高,文獻中並針對瑞利數在 109~1012範圍下對紐塞數之做一關聯式。

Ji 等人[9]2007 研究了暫態混合對流在一個驅動腔體。此腔體頂部溫度高於底部 的穩定分層的結構,左右壁面都為絕熱,在上壁面啟動後提供一由左向右的水平速 度。實驗進行不同的格拉斯霍夫數(Gr)與雷諾數(Re)對穩態層流流場之影響。在穩定 狀態下,流場會形成一個大循環結構在腔體中央,其溫度分布均勻。當 Ri 增加時,

在穩定狀態下,原本的大迴流結構會被擠壓到頂部,形成一狹小迴流區域,造成溫度 分布在腔體的下方,形成穩定的熱分層。實驗也針對非穩態紊流場下,記錄溫度在初 期的發展上,於中央頂面與底面附近微弱的波動現象。

Angeli 等人[10]於 2007 年在方形孔穴內加入一個熱水平圓柱進行了數值模擬。

採用二維假設,觀察數值結果呈現出穩定對稱、穩定非對稱、非穩態週期與非穩態擾 動等四種流動現象,研究中利用瑞利數與直徑與高度比進行分類,其熱傳分析結果也

(14)

與其他的文獻作一比較,呈現相當一致之結果。

Yilmaz 與 Fraser[11]於 2007 年在垂直通道中,由一個等溫度加熱牆與一垂直面 玻璃牆組成一二維流道進行多普勒測速儀(LDA)記錄,並在通道上方出口測量速度 與溫度分布,並以商業軟體採用在低雷諾數 K-ω紊流模型進行數值與實驗比較,數值 方法,能相當準確的預測熱傳與出口流量,也能捕捉剖面速度、溫度與紊流動能分佈,

並將熱傳量與流量結果歸納出對應之方程式,並透過圖形與其他文獻做比較。

Wu 等人[12]於 2008 年針對方形空氣孔穴內左熱右冷之自然對流場中將頂面加熱 改變其流場結構,實驗採用格拉夫數(Grashof number)1.3×108,且上壁無因次溫度 由 1.4 至 2.3。量測結果顯示,流線受上方穩定溫度分層之結構影響,在到達上壁面 角落前提早分離。因此孔穴內形成三個區域:分層的核心區域、浮力效應下流線上升 之區域、與在分離上方穩定溫度分層區域。當頂面越熱其高度分層的流線結構區域越 大,使得在上壁溫度更穩定的增加。這個現象相似 Kulkarni 等人[16]在自然對流熱 分層中恆溫垂直壁面與 Chen 與 Eichorn[17]在垂直表面自然對流熱分層流體結構的 熱傳分析。

Sleiti [13]於 2008 年以數值研究冷熱空氣受浮力於一交換開口之暫態流動現 象,其上部腔體為冷空氣,下部腔體為熱空氣,因冷熱在不穩定流場模式下,此一水 平交換開口之長寬比對流體交換具有相當大之影響。結果討問包含 L/D=1、0.5、2.0。

L 為開口長度、D 為狹縫寬度。結果顯示,在 L/D=0.5 時流量交換會比在 L/D=2.0 時 增加許多。針對不穩定模式下的流動頻率與瞬態的衝擊現象在上游和相對應的下游腔 體做了記錄並與文獻實驗數據做一比較。

Trias 等人[14]於 2010 年以直接數值模擬在寬高比為 4 之左右壁加熱腔,進行 瑞利數 6.4x108, 2x109, 1x1010, 3x1010, 與 1x1011等五個條件進行時間平均流場分析 與比較,其瑞利數涵蓋弱紊流至完全發展之紊流。其流場明顯觀察在高瑞利數時邊界 層的過渡區會往上游移動,並提早產生分離,因此會造成頂部和底部之迴流區增加,

萎縮中央溫度分層區域。因此熱分層值明顯大於 1。

(15)

Santos 等人[15]2010 年研究大渦流模擬二維方形孔穴內紊流與層流之混合對流 暫態分析。採用 Smagorinsky 模型方式處理。使用有限元素法解 Navier-Stokes 方程 式,雷諾數為 400 到 10000,理查森數為 0.1 到 0.44 之範圍。紊流與層流模擬結果與 文獻相比非常接近,最大誤差量僅 5%。結果顯示,因浮力造成之分層結構對流場與 熱傳有很大的影響。

1.3 採用方法

本研究延續先前有限體積法之數值方法,將二階之速度與溫度微分更改為四階計 算,由先前學長[18]之研究可知擴散項若採高階能免除人工黏滯之需求,能提供較精 確之數值解答,因此本論文將二階之有限體積流場解答,擴充至四階等向性元素之微 分法於擴散項之計算。另外,本研究之結果需與文獻不可壓縮流數據比對,因此採用 低馬赫數設定之混合對流流場,因此需運用預置矩陣之方式,以縮小可壓流場之特徵 值差異過大以避免收斂困難之發生,另外也加入時間項,並與預置矩陣整合以符合暫 態流場時間精度之需求。研究旨在開發一暫態分析程式,能應用於如同先前文獻所探 討之現象,以掌握流場與熱傳之瞬時資訊,以協助未來在其它工程上之應用,另外也 透過直接數值模擬之跳動統計,了解速度與溫度在邊界層內傳遞之行為。

1.4 文章安排

本研究論文分為五個章節,架構簡介如下:

第一章:為緒論部份,說明本論文先前有關文獻與採用方法。

第二章:介紹本研究中所探討之物理現象與無單位係數定義。

第三章:本章介紹統御方程式與數值方法並對時間積分、有限體積法、人工黏滯、邊 界條件的設定並加以說明。

第四章:為模擬結果驗證,將數值模擬答案與參考文獻之結果進行比對,並探討不同

(16)

第五章:為結論及未來展望,歸納研究中所得之成果並說明未來發展之方向。

(17)

第二章 物理問題描述 2.1 混合對流問題與邊界條件

本論文探討的問題是上下壁具有溫度差,而左右垂直壁面為絕熱,且上壁面向右 等速移動之正方形腔體,由於浮力與剪力之作用,對腔內空氣產生之二維混合對流現 象。當底面熱頂面冷為非穩定熱分層,相反則為穩定熱分層,如圖 2 所示。其溫度與 速度之邊界條件處理如同圖 2 所示,但對可壓縮流程式需加入第四個條件以符合方程 式數目,在此程式以動量方程式外差其壁面壓力,並於每次疊代後平移每一點之壓力 值以符合孔穴內中央處壓力為一定值。

2.2 大尺度渦流數值模擬

採用大尺度渦流是基於 Boussinesq 之渦流黏滯假設,其紊流應力在局部下,生成 與消散維持平衡,當不可壓縮或近似不可壓縮流下其紊流應力可寫成:

 

 

∂ + ∂

= ∂

i j j

i sgs

ij

x

V x

ν V

τ

(2.1)

其中

ν

sgs為動態渦流黏滯係數,當使用 Smagorinsky 之次格點模式為,其動態渦流黏 滯係數定義如下:

ν

sgs

= C

s2

2

S

(2.2)

其中

C

s為 Smagorinsky 常數;

為次格點之特徵長度,

S

為流場過濾後之應變率,

目前計算所採用之特徵長度以下式處理

2 2

=1

=

i

x

i (2.3) 而應變率之計算以 box 過濾進行

ij ij

S S

S = 2

(2.4)

(18)

其中

S

ij

= 2 1 ( ∂ V

i

x

j

+ ∂ V

j

x

i

)

,至於 Smagorinsky 常數一般介於 01.1~0.2,利 用三維等向性紊流之半理論方式推導。至於紊流細下之傳熱量亦可表示為:

i sgs

i

x

q T

= α ∂

(2.5)

其中

α

sgs為渦流熱擴散係數,可利用擾動普朗特係數結合動態渦流黏滯係數計算如 下式所示

C S

sgs s sgs

2 2

Pr ∆

α =

(2.6)

其中∆為網格之特徵長度, S 為流場過濾後之應變率,

C 為紊流常數,為紊流之普

S 朗特數(Prandtl number)。在本研究中,=0.9 、格拉斯霍夫數為 7.04×106、雷諾數 為 4000、普朗克數為 6、理查森數為 0.44 為固定值,分析在紊流常數為 0 與 0.1 時,

其流場暫態下流場之速度、溫度變化,以及流場結構與溫度分佈。

2.3 基本假設

在本論文中,我們採用二維 Navier-Stokes 方程式描述孔穴內空氣之流動與熱傳,

並做了以下之假設 1.流場為層流。

2.空氣為理想氣體

( P

=

ρ RT )

3.黏滯係數

( ) µ

、隨溫度變化之關係可透過 Sutherland’s Law 近似。

4.熱傳導係數

( ) k 由普朗特數計算 k

=

µ C

p Pr,其中Pr=6

5.大尺度渦流計算時,Smagorinsky 常數

C

s=0.1,紊流普朗特係數

Pr

sgs=0.9。

(19)

2.4 無因次參數

由於傳統不可壓縮流之流體膨脹係數為一定值,但在可壓縮流下,氣體於孔穴中

每一點溫度皆不一致,因此採用冷熱壁面之平均溫度

( )

2

0

bot

top

T

T

=

T

+ 定義流場之

平均膨脹係數,同理流體黏滯係數與熱傳導係數於葛拉斯赫夫數計算時採用平均溫度 (

T

0)處理,其葛拉斯赫夫數(Grashof number)定義如下:

2 0 0

3

ν T

TH Gr g

=

(2.7)

其中,

g

為重力加速度,

T = T

bot

T

top在底面與頂面之溫差,

H

正方形孔穴之深

度,

ν

0為空氣在平均溫度下之動黏滯係數。而上蓋驅動速度之大小以無因次之馬赫 數表示如下:

top top

RT M U

= γ

(2.8) 而孔穴內流動之雷諾數( Reynolds number ) 亦採用平均溫度( T0 )處理,其定義如下:

Re ν

0

H U

top

=

(2.9)

此外也可透過理查森數(Richardson)結合雷諾數與葛拉斯赫夫數如下:

Re

2

Ri = Gr

(2.10)

至於分析壁面熱傳導量之多寡上,皆採用無因次之紐塞數( Nusselt number)表示,其定 義如下:

n T TK Nu HK

= ∆

0

(2.11) 其中

n

為壁面之法線方向。至於平均紐塞數可由壁面積分計算之

(20)

ds L Nu

Nu 1

L

0

=

(2.12)

最後本研究中在記錄流場發展過程不同的時間下之流速與溫度分佈。皆採用無因次之 時間尺度,其定義如下:

H

t

*

= tU

top (2.13)

(21)

第三章 數值方法 3-1 Navier-Stokes 方程式

本研究之數值方法求解二維非穩態 Navier-Stokes 方程式,其守恆方程式可表示 如下:

=0

− +

+

S

y F x E y F x E t

U

c c v v

(3-1)

其中

U

E

c

F

c

E

v

F

v

S

如附錄 A 所示。透過虛擬時間項以及預置矩陣 (Preconditioned)之導入可以維持每一物理時間下之時間精度,並在每一格以其局部之 時間步階值(τ )進行疊代以加速收斂,方程式(3-1)可改寫成

=0

− + +

∂ +

Γ∂

S

y F x E y F x E t

U

p

U

c c v v

τ (3-2)

其中

U

p =

[ p

,

u

,

v

,

T ]

T,守恆變數 U 與預置變數 Up間之間可以透過變數矩陣轉換 U=MUp。而Γ為預置矩陣,最早由 Choi 與 Merkle [19]所提出以改善對流項之特徵 值以避免低馬赫數之收斂困難狀況,而後續 Weiss 與 Smith [20]之研究。建議

[ ]

[ ]

[ ] ( )

1 0 0

1 0

1 0

1 p

RT T

u RT u T

v RT v T

H RT u v C H T

ρ

ρ ρ

ρ ρ

ρ ρ ρ

Θ + −

 

 Θ + − 

 

Γ =  Θ + − 

 

 

Θ + −

  

 

(3-3)

其中Θ =1Vref2 −1c2C為聲速,本文參考 Edwards and Liou [21]所提出的參考速度

Vref ,其定義如下

V

ref2 =min

c

2, max

( V

2,

KV

max2

)

(3-4) 其中

V

max為正方形孔穴內之最大流速而

K

為一常數。

(22)

3-2 時間積分

在本研究中在物理時間項採用二階 backward 差分,並配合預置變數之積分,方 程式(3-2)可以展開為下列式

2 0 4 ) (

3 1

=

− +

∆ +

+

∆ + +

∂ Γ∂

y S F x E y F x E t

U U U

M U

U

c c v v

n n p

k p

τ

(3-5)

其中

t

U M

p

∆ 2

3 項可與預置矩陣項合併成為新的預置矩陣 M

t + ∆ Γ

=

Γ 2

' 3 τ

,最後虛擬時 間項則以三階 Runge-Kutta 方式進行整體時間積分以獲得非穩態流場解答

) ( ) 3 ( 2 3 2 3 1

) ( ) 4 ( 1 4 1 4 3

) ( ) (

2 1 ' 2

1

1 1 ' 1

2

1 ' 1

U R U

U U

U R U

U U

U R U

U

n n

n

n n

+

Γ

∆ + +

=

Γ

∆ + +

=

Γ

∆ +

=

τ τ τ

(3-6)

其中∆ 為虛擬時間,而

τ

S

y F x

E y

F x

E t

U U

R U

c c v v

n n

k

  +

 

∂ + ∂

∂ + ∂

 

 

∂ + ∂

− ∂

∆ +

− −

=

2 4

3

1

(3-7) 至於局部所需之疊代數目,亦即收斂與否乃透過(3-7)式之殘值量來做判斷,以確保積 分後之 Uk變數能符合所設定之暫態精度(10-7)。

3-3 流通量計算

二維流場於右側截面之流通量有

j , 2 i+1

E

j 2, i+1

F ,將採用平均值計算,則如下:

( )

(

i,j i 1,j

)

j 2, i 1

j , 1 i j , j i

2, i 1

2 1 2 1

+ + + +

+

=

+

=

F F F

E E E

(3-8)

因此,該截面積上之總流通量為

j y

x

F

i

An

An E

Flux

2, j 1

2,

i+1

+

+

=

(3-9)

(23)

其中,

A

為該流通量之投影面積。

若以四邊形之網格而言,該計算體積之總流通量為

+ +

=

+ +

x y

side

An F

An E

Flux

j 2, i 1 j

2, i 1 4

+ +

=

x y

side

An F

An E

Flux

j 2, i 1 j

2, i 1 4

+ +

=

+ +

x y

side

An F

An E

Flux

2 j 1 , 2 i

j 1 , 4 i

y x

side

An F

An E

Flux

2 j 1 , 2 i

j 1 ,

4

=

i

+

(3-10)

3-4 穩定性分析

在使用顯式法時,未避免發散採用 von Neumann 法來作分析。考慮統御方程式如 下:

c c v v

E F E F

U

t x y x y

∂ ∂ ∂ ∂

∂ + + = +

∂ ∂ ∂ ∂ ∂

在 x-y 座標下,所採用的計算格點為非均勻網格,為便於分析則將物理 x-y 座標系統 轉換為計算均勻網格點

ξη

座標系統,其結果如下:

c c v v

E F E F

J U

t ξ η ξ η

′ ′ ′ ′

∂ ∂ ∂ ∂

∂ + + = +

∂ ∂ ∂ ∂ ∂

(3-11)

其中

U,F ,G ,F ,G , J

c

′ ′ ′ ′

c v v 如附錄 B-1。(3-7)式為非線性方程式,將其推導為線性方程 式:



 

∂ + ∂

∂ + ∂

∂ + ∂



 

∂ + ∂

− ∂

∂ =

η η ξ

η ξ ξ

T Q S Q

S Q A Q

A Q t

Q 2

2 1 2 2 2

2 1 2

1 (3-12)

(24)

而 Q A F

Q ,

A1 Ec 2 c

∂ ′

∂ =

∂ ′

= 其中

Q, A

1

, A

2

, S

1

, S

2

, T

1如附錄 B-2。將中央差分法應用於

上式可得:

( ) ( )

 

 

− −

− + ∆

∆ + + −

∆ + + −

 

 

∆ + −

− −

∆ =

+

+ +

+

+

+

+

+ +

η η

ξ

η ξ

η ξ

2 Q Q

2 Q Q

2 T

Q 2Q S Q

Q 2Q S Q

2 Q A Q

2 Q A Q

t Q Q

n 1 j 1, i n

1 j 1, i n

1 j 1, i n

1 j 1, 1 i

2 n

1 j i, n

j i, n

1 j i, 2

2 n

j 1, i n

j i, n

j 1, i 1

n 1 j i, n

1 j i, 2 n

j 1, - i n

j 1, i 1 n

j i, 1 n

j

J

i,

(3-13)

使用 von Neumann 法:令

Q

ni,j

= e

at

e

ikxξ

e

ikyη代入(3-9)式且等號兩邊各除以Q ,則可ni,j

( ) ( )

 

 

 − − −

+

+ + −

+ + −

 

 

 − + −

− =

η η

ξ

η ξ

2 Δ 2 Δ

2 Δ T

Δη S 2

Δξ S 2

2 Δ 2 Δ A

Δt A

Δ Δ Δ

Δ Δ Δ

1

2 Δ Δ

2 2 Δ Δ

1

Δ Δ

2 Δ Δ

1 t

η η ik

ik ξ ik η η ik

ik ξ ik

η η ik

ξ ik ξ ik

ik

η η ik

ξ ik ξ ik

a ik

y y

x y

y x

y x y

x

y x y

x

e e e

e e e

e e

e e

e e

e e

I J e

其中 I 為單位矩陣,

k

x

, k

y為誤差相角。

再應用三角函數關係

2 e sin e

2 , e cos e

i i

i

iθ θ θ θ

θ

θ

= + = 可得

(25)

( ) [ ] ( ) [ ]

) sin

) Δ sin

T

1 ) 2S cos

1 ) 2S cos

) A sin

) A sin

Δt

1

2 y 2 2

1

2 1

t

η η ξ

ξ

η η ξ ξ

η η ξ ξ

∆ ∆

∆ ∆ +

∆ ∆ +



 

 ∆

+ ∆

∆ ∆

− =

y x

x

y x

a

(k (k

(k (k

(k i (k

I i J e

化簡上式可得:

[ ]

( ) ( )

[ ]

) sin ) sin T

1 cos

2S 1

cos t 2S

) sin A ) sin t A

I

y 1

y 2

1

y 2

1 t

(k (k

k J k

i (k J (k

e

x x

x a

− +

∆ − +

∆ +

=

(3-14)

上式

e

at為擴大矩陣

t Z) P(

G

G, J

= ∆

將(3-10)式對流項與擴散項兩部份探討,以求

得最小時間步階,設定Z=ZI +ZR則:

[ ]

( ) ( )

[ 2S cos 1 2S cos 1 T sin ) sin ) ]

Z

) sin A ) sin A Z

y 1

y 2

1 R

y 2

1 I

k ( k

( k

k

k ( k

(

x x

x

− +

=

+

=

為符合穩定性條件

λ

G ≤1則:

1. 考慮

t Z ) 1 (

I

 ≤

 

 ∆

J

,亦即 t

( )

Z ) 1

( I max  ≤

 

 ∆

λ

J

Z

I如附錄 B-2 所示,則

Z

I的四個 特徵值為:

2

2

b

a c U , U ,

U ± +

k

x

, k

y

π 2

時,可得最大之特徵值為:

( ) Z

I max

= MAX [ U ± c a

2

+ b

2

]

λ

因此對流項時間步階

∆ t

c為:

( )

2

t = 2 ×

J

(26)

其中2 2為四階 Runge-Kutta 之容許值。

2.

t Z ) 1 (

R

 ≤

 

 ∆

J

,則

( ) ( )

[ ]

( ) ( )

[ 2S cos 1 2S cos 1 T sin ) sin ) ]

t

) sin ) sin T 1 cos 2S 1 cos t 2S

t Z

y 1

y 2

1

y 1

y 2

1 R

k ( k ( k

J k

k ( k ( k

J k J

x x

x x

− +

∆ −

− +

∆ −

∆ =

當kx,ky

π

時,cos

k

x −1,cos

k

y −1有最大值 2,而當

k

x

, k

y

π

2時,

) sin ),

sin ( k

x

( k

y 有最大值 1,因此

[

1 2 1

]

R t 4S 4S T

t Z ≤ ∆ + +

J J

由於T 的特徵值非常複雜,因此假設格點是近乎直交的,則1 T 的四個特徵值如下: 1

( ) ( ) ( )

( ) ( )

ξ

ρ η

µ λ µ

λ λ

λ

x 3 y

T

0 T T

T

t l 1 4

1 3 1 2 1 1

J

= +

=

=

=

因此

( )

ξ

ρ

η

µ

µ

y x

T1 3l t

J

= +

S

1的四個特徵值為

( ) ( ) ( ) ( )

1

( ) (

2 2

)

4

1 3 1 2 1 1

x 1 y

S

0 S S

S

η

ρ

η

λ

λ λ

λ

− +

=

=

=

=

RJ γ K

因此

( ) [ (

2 2

) (

2 2

) ]

2

1

4 1 y x y x

S 4 S

4

η η ξ ξ

ρ

+ + +

= −

+ RJ

γ K

所以

( ) ( ) [ ( ) ( ) ]

 

 

 + + − + + +

≤ ∆

l t 2 2 2 2

R

4 1 y x y x

x 3 y

Z t t

ξ ξ η η ξ

η ρ

ρ µ µ

RJ γ K J

J J

則擴散項的時間步階為:

(27)

R

v

Z

t = J

(3-16)

因此時間步階取兩者中最小之值。

t = MIN ( ∆ t

c

, ∆ t

v

)

(3-17)

3-5 人工黏滯

在數值計算時,若加上人工黏滯可避免數值震盪。本論文使用 Jameson 的人工 黏滯[21]。考慮一維模擬之方程式

= 0 +

x

t

cu

u

(3-18)

透過人工黏滯之導入

xxxx xx

x

t

cu v xu v x u

u + =

2

∆ −

4

3 (3-19)

而(3-14)式之解為

u = e

ikxiωt,將其帶入(3-15)式,則可得

( )

(

v2 xk2 v4 x3k4

)

kct x

i

e e

e

ω

=

+ (3-20)

此時,當

v

2

v

4之值皆大於零時,則可使數值震盪現象呈現指數衰退。如果 在方程式中加入

v

3

x

2

u

xxx,則將產生

e

ix2k3v3t的項,此項只會修改解答相位的改 變,對於避免數值的震盪並無實質上幫助。由於人工黏滯

u

xx 項的加入,產生一階的 誤差,而

u

xxxx項則造成三階的誤差,因此

v

2不可取的過大,以免影響數值結果的 精確度,本研究只使用四階微分以避免數值震盪的現象,其值皆固定在 0.001。最後,

加上人工黏滯的方程式為

( )

x,y

y , x y ,

x

Flux AV

dt

V dU + ∑ =

(3-21)

( )

(28)

3-6 固體壁面

考慮黏滯性流時,由於在固體壁面上之無滑動條件,所以流通量之計算可表示成:

( ) ( )

( )

0

0

wall xx xy

xy wall yy

P dy

E F dS

P dx

τ τ

τ τ

 

  − −  

   

+ =  

 − + − 

   

 

 

∫ ∫

(3-22)

其中

P

wall為壁面上之靜壓力,

τ

ij 為壁面上之剪應力。

(29)

第四章 結果與討論 4.1 不同理查森數之影響

本節所探討之問題,共有三組參數進行數值模擬,透過 Ji 等人[9]之實驗數據作 層流暫態流場之驗證。其中雷諾數為 400、普朗特數為 6,理查森數分別為 Ri=0.1、

Ri=1 與 Ri=10 時之混合對流現象,並記錄下腔體中之溫度與速度之發展過程以及二 維瞬時溫度與速度分布,作為程式開發之依據。

腔體內之暫態溫度與速度發展過程,取在腔體中央斷面(X=0.5)於垂直高度位在 Y=0.27、Y=0.48 與 Y=0.93 等三處進行記錄,並與文獻實驗結果進行比較。圖 4-1 為 理查森數為 0.1 時之暫態溫度與速度變化過程。由圖 4-1 中六張圖形可觀察出目前數 值結果與實驗數據非常吻合,但在 0.48 處之水平速度於 t*>10 之後目前結果量測之 差異相當明顯,由流速與溫度初期之發展主要受到頂面剪切流帶動之影響,浮力效應 不明顯,且水平速度在 t*>50 之後幾乎已沒有太大變化,顯示流場結構已趨於穩定,

但反觀上中下三處之溫度變化,即使 t*>100 後仍有明顯之改變,顯示熱傳與浮力效 應之傳遞較為遲緩。

圖 4-2 為理查森數為 1 時之暫態溫度與速度變化過程。由這些圖可觀察出目前數 值結果與實驗數據仍相當雷同,但在 0.48 處之溫度發展於 t*>50 之後與量測結果差 異較大,另外在 0.93 處初期速度之發展於 t*<25 期間有相當落差,由流速與溫度初 期之發展雖然受到頂面剪切流帶動之影響,但因浮力效應增加會抑制剪切迴流之發 展,雖流速與先前 Ri=0.1 弱,但其發展延遲許多,因此水平速度在 t*>200 之後仍持 續在改變,但反觀上中下三處之溫度變化,t*>150 後變化已趨緩,顯示浮力效應已 發展完成但對流項仍在持續作用。

圖 4-3 為理查森數為 10 時之暫態溫度與速度變化過程。由這些圖可觀察出目前 數值結果與實驗數據仍接近,但在 0.48 處之溫度發展於 t*>50 之後與量測結果差異 較大,另外在 0.93 處速度之發展皆低於實驗值,由於浮力效應更強,因此剪切迴流

(30)

t*>30 後已達穩定之低流動訴度甚至呈現靜止之狀態,由於對流項很弱,因此溫度發 展初期看不出迴流帶動之現象,以熱擴散為主,其發展也變得很緩慢,t*>200 溫度 仍在持續改變,顯示熱平衡需非常漫長。

圖 4-4 為頂面與底面之平均紐塞數之變化過程,由此圖之比較可明顯觀察出,

Ri=0.1 之熱傳量遠高於 Ri=1 與 Ri=10 之流場,頂面與底面之平均紐塞數已趨於一 致,顯示熱平衡已發展接近完成,對 Ri=0.1 之頂面平均熱傳量初期雖有剪切迴流之 攪拌,但迴流並沒有貫穿至底面因此效果不足,因此仍有賴熱擴散持續進行熱平衡,

至於 Ri=10 之頂面平均熱傳量初期發展顯示更多之迴流攪拌,但時間更短效果更弱,

顯示迴流幾乎沒有貫穿孔穴,無法將冷熱空氣做一快速交換,因此在頂面與底面之平 均紐塞數之發展上呈兩條平行線之模式在發展,距離熱平衡將要耗費相當多之計算時 間。

圖 4-5 至 4-7 為在不同的理查森數下,腔體暫態之流線分布,與 Ji 文獻之數值 結果對照相當一致。其中 Ri=0.1 時其流場結構受到由右向左的速度帶動下,腔體之 右上方形成順向渦流結構,隨時間發展下,渦流結構逐漸發展至腔體中央並分布在腔 體的大部分區域,其中 t*=5 之結果與 t*=200 之流動已非常接近,顯示對流效應發展 迅速。但在圖 4-6,Ri=1 時驅動之渦流初期被限制在右上方處,且下方有另一反轉渦 流阻礙頂面與底面之傳熱,且隨時間發展,浮力效應較大之情況下,t*=5 時近底面 處又形成第三個渦流之結構,在 t*=200 時主迴流有加大,略為擠壓另外兩個渦流形 成略為偏置之造型,至於圖 4-7,Ri=10 下之流線分布顯示主渦流之貫穿能力很弱,於 t*=1.4 之後皆一直處在上方四分之三處,因穩態分層之效應孔穴內不斷產生逆向與 順向渦流,在 t*=200 時已發展出五層之迴流區塊,因此對流效應被減弱,且流速也 很低,因此個別渦流對局部均溫性也無法產生效果,所以要達到熱平衡需要長時間的 擴散。

圖 4-8 至 4-10 為在不同的理查森數下,腔體暫態之溫度分布,與 Ji 文獻之數值 結果對照也是相當一致。其中圖 4-8Ri=0.1 時其溫度分布受到順向渦流之帶動,在

(31)

t*=5 時已有一半的冷熱空氣進行攪扮,在 t*=200 時其溫度場已非常接近穩態分布。

圖 4-9,Ri=1 之溫度初期分布結構與先前之流線分布一致,顯示在 t*<5 期間主要受渦 流之影響,但在 t*=200 時溫度場與三個渦流結構不一致,顯示底面溫度之發展主要 受擴散作用而非對流效應。圖 4-10,Ri=10 之溫度分布結構與先前之流線分布差異頗 大,僅在右上角有一渦流結構外,其餘區域皆維持相當水平之熱分層,當 t*=200 時 頂部溫度已趨於均溫,但中央至底面仍維持一穩態熱分層之結構。

4.2 暫態流場下直接數值分析與大尺度渦流之比較

使用雷諾數為 4000、格拉斯霍夫數為 7.04×106、普朗特數為 6 時,觀察腔體內 之暫態溫度與速度發展過程,取在腔體中央斷面(X=0.5)於垂直高度位在 Y=0.28、

Y=0.5 與 Y=0.94 等三處進行記錄,並與文獻實驗結果進行比較。分別以直接數值分 析與採用大尺度渦流模擬法之紊流模型,其 Smagorinsky 常數為 0.1、紊流之普朗特 數為 0.9,兩組結果進行暫態流場分別計算。

圖 4-11 為不同位置下之瞬時水平速度變化過程,瞬時水平速度在發展初期僅在 t*=20 與 t*=30 才受到順向渦流之影響。比較直接數值分析與採用 Y=0.94 處有大振 幅跳動現象,顯示出主迴流之影響,但在 Y=0.5 與 Y=0.28 處分別在大尺度渦流模擬 在 t*<50 之前期差異很小,但伴隨渦流之發展與流場之攪動,大尺度渦流模擬之流動 其波動量有減小之趨勢,且流速較快趨於時間平均值。

圖 4-12 為不同位置下之瞬時溫度變化過程。比較 Ji 等人之實驗數據,目前模 擬結果前段有些許落差,不過 觀察實驗數據可發現其起始溫度並非模擬所設定之線性 分布,尤其是 Y=0.94 之啟始值並非 0.94 而是落在 0.78 左右,同樣 Y=0.28 之啟始值 並非 0.28 而是落在 0.22 左右,但溫度在後續之發展上,數值與 與實驗結果逐漸吻合。

比較大尺度渦流與直接數值分析之結果,發現採用大尺度渦流分析瞬時溫度之振幅較 小,在瞬時溫度分布上與實驗結果也較為接近,反觀直接數值分析在 Y=0.94 與 Y=0.28

(32)

分析,有略好之準確性。

圖 4-13 不同時間下之瞬時流線分布圖,比較大尺度渦流之結果相較直接數值分 析,在初期 t*=20 差異很小,但在 t*>40 之後,其主迴流擴展較快,且反轉渦流較小 且被擠壓至更靠近壁面,在 t*=200 時其流場較趨於時間平均之流動,反觀直接數值 分析其主渦流距底面仍有相當距離,且右下角與左上角之反轉渦流仍未形成。

圖 4-14 為不同時間下之溫度分布圖,大尺度渦流之結果相較直接數值分析,在 初期 t*=20 差異很小,但在 t*>40 之後,其主迴流擴展較快,造成溫度混合更快速,

因此溫度梯度較高且複雜,但到達最後 t*=200 時,大尺度渦流之溫度分布則較平順 且熱邊界層更靠近壁面,顯示出比直接數值分析較快接近時間平均之溫度分布。

4.3 穩定熱分層與非穩定熱分層之影響

由於大尺度渦流在前面紊流場分析上較接近實驗結果。因此在雷諾數為 10000、

格拉斯霍夫數為 10×107、普朗特數為 6、理查森數為 0.1 下,將以大尺度渦流法處理。

在此紊流場下針對頂面熱與冷(stable)以及頂面冷與底面熱(unstable)兩種條件進 行暫態分析。

圖 4-15 為穩定熱分層下於 t*=5.6、t*=9.4、t*=13.1、t*=16.9、t*=30.0、

t*=40.9、t*=60.8、t*=72.2 等八個時段之瞬時溫度分布。其中圖左側為目前結果,

右側為 Santos 與 Piccoli 文獻之結果,比較此一系列圖形,目前分析之結果發展比 較迅速,且當後段時不均勻的冷熱氣泡於孔穴內分布較多且明顯。在 t*=5.6 時,受 到頂面剪切流帶動,在右上角形成一主要渦流;t*=9.4 渦流逐漸往下移動並開始出 現不規則之捲繞現象;t*=13.1 主渦流分裂成兩個捲曲結構且向上漂移而非貫穿下方 冷流體;t*=16.9 至 60.8 時渦流開始侵蝕底部冷流體分層,且渦流不斷捲繞與分裂 形成許多更細更小之渦流,有助於熱流體與冷流體混合效果。當到達 t*=72.2 時,混 合過程接近完成,但尚未達到完全混合。

圖 4-16 為非穩定熱分層下於 t*=5.6、t*=9.4、t*=13.1、t*=16.9、t*=30.0、

參考文獻

相關文件

If that circle is formed into a square so that the circumference of the original circle and the perimeter of the square are exactly the same, the sides of a pyramid constructed on

Which of the following is used to report the crime of damaging the Great Wall according to the passage.

The case where all the ρ s are equal to identity shows that this is not true in general (in this case the irreducible representations are lines, and we have an infinity of ways

Miroslav Fiedler, Praha, Algebraic connectivity of graphs, Czechoslovak Mathematical Journal 23 (98) 1973,

The min-max and the max-min k-split problem are defined similarly except that the objectives are to minimize the maximum subgraph, and to maximize the minimum subgraph respectively..

Experiment a little with the Hello program. It will say that it has no clue what you mean by ouch. The exact wording of the error message is dependent on the compiler, but it might

and Dagtekin, I., “Mixed convection in two-sided lid-driven differentially heated square cavity,” International Journal of Heat and Mass Transfer, Vol.47, 2004, pp. M.,

In flow field boarding up for the smooth wall surface, gets down the board is the Sinusoidal wavy wall, but under board wavy wall oscillation amplitude A=2.54(mm), the wave length λ