一垂直風機之三維模擬

國

國

國

國 立

立

立

立 交

交

交

交 通

通

通 大

通

大

大 學

大

學

學

學

機

機

機

機 械

械

械 工

械

工

工

工 程

程

程

程 學

學

學

學 系

系

系

系

碩

碩

碩

碩 士

士

士 論

士

論

論 文

論

文

文

文

一垂直風機之三維模擬

3D Simulation of a Vertical-Axis Wind

Turbine

研

研 究

研

研

究

究

究 生

生

生：

生

：陳

：

：

陳

陳 虹

陳

虹

虹

虹 汝

汝

汝

汝

指

指

指

指 導

導

導

導 教

教

教 授

教

授

授

授：

：崔

：

：

崔

崔

崔 燕

燕

燕 勇

燕

勇

勇

勇 教授

教授

教授

教授

中

中

中

中 華

華

華 民

華

民

民

民 國

國

國 101 年

國

年

年

年 6 月

月

月

月

一垂直風機之三維模擬

一垂直風機之三維模擬

一垂直風機之三維模擬

一垂直風機之三維模擬

3D Simulation of a Vertical-Axis Wind Turbine

研 究 生：陳虹汝 Student：Hung-Ju Chen 指導教授：崔燕勇 Advisor：Yeng-Yung Tsui 國 立 交 通 大 學 機 械 工 程 學 系 碩 士 論 文 A Thesis

Submitted to Department of Mechanical Engineering

College of Engineering National Chiao Tung University

in partial Fulfillment of the Requirements for the Degree of

Master in

Mechanical Engineering June 2012

Hsinchu, Taiwan, Republic of China

I

一垂直風機之三維模擬

研究生：陳虹汝 指導教授：崔燕勇 博士 國立交通大學機械工程學系﹙研究所﹚碩士班 摘 摘 摘 摘 要要要要 本研究係針對一結合升力型與阻力型風機特性的新型式垂直風機進行 氣動力、性能與流場結構的研究。此種風機之葉輪具有三個 S 型葉片與兩 個翼形葉片，並可將數個葉輪堆疊進而得到較佳的性能。為瞭解其氣動力 特性，利用計算流體力學進行數值模擬，數值方法採用有限體積法，並利 用移動網格處理風機葉片的旋轉，搭配 RNG k−ε紊流模式進行三維與二維 暫態模擬。除數值模擬外，另進行實驗予以驗證模擬之準確性。 三維模擬結果顯示此風機之扭矩係數與風機端速比呈線性遞減關係， 而功率係數則與端速比呈二次式關係。與單層葉輪的實驗結果比較可發現， 二維暫態的計算高估了扭矩的輸出，而三維暫態模擬則與實驗結果相當吻 合。從流場中可發現，自風機後緣流出的氣流，在二維計算中形成一週期 性變化的渦流，此渦流逐漸往下游逸散，三維計算中則呈現一穩定的流場， 這表示在三維模擬中，葉輪上下的流場有穩定風機尾流之效果，此三維效 應是二維模擬無法反應的。此外，模擬結果亦顯示，將葉輪堆疊而成的多 層風機較同高度之單層風機僅具稍高的輸出功率，但扭矩隨時間的變化則 隨風機層數的增加而減少，對於穩定風機旋轉有相當大的助益。

II

3D Simulation of a Vertical-Axis Wind Turbine

Student：Hung-Ju Chen Advisor：Dr. Yeng-Yung Tsui

Department of Mechanical Engineering National Chiao Tung University

ABSTRACT

A vertical-axis wind turbine which combines the merits of the lift type and drag type wind turbines is studied in this project. It consist of three S-type and two airfoil-type turbine blades. Several stages of these turbine blades can be stacked together to achieve better performance. The aerodynamic characteristics and performance are investigated by using computational fluid dynamic. The numerical scheme is based on the finite volume method. The moving mesh technique is adopted to tackle the rotation of the turbine wheel.With use of the RNG k−ε _{model to characterize turbulence, unsteady calculations are}

undertaken. Experiments are also conducted to assess the accuracy of the numerical simulation.

Two-dimensional and three-dimensional, unsteady simulations are conducted. It is shown that the resulting moment coefficient becomes a linearly decreasing function of the tip speed ratio and the power coefficient is a quadratic function of the tip speed ratio. Comparison with experimental measurements for a one-stage wind turbine indicates that good agreement is obtained by the 3-D unsteady simulation. However, the torques and powers are over-predicted by the 2-D unsteady calculations. It can be detected from the results that after the flow passing the wind turbine, a periodic shedding vortex is formed in 2-D calculations while a relativity steady flow is observed in 3-D calculations. This implies that the flows beyond the upper and lower sides of the wind turbine have an effect of stabilizing the flow in the near wake of the turbine. This 3-D effect cannot be accounted for in 2-D calculations.

Comparing with the single-stage wind turbine, only slightly higher power performance can be obtained by stacking several stages of turbines together. However, the variation of the torque with respect to time becomes smoother as the number of stage increase. This character helps stabilize operation of the wind turbine.

III

誌謝

誌謝

誌謝

誌謝

感謝崔燕勇教授在研究所期間的諄諄教誨，除了課業上的知識外，還有 許多做事的方法與態度。也感謝金大仁教授實驗室的蘇宏明、邱奕豪、黃 振瑋與其他同學，在實驗上幫助了我許多。 另外感謝林仕文學長在學習過程中的指導，感謝林子翔學長在研究上的 指導，使我的研究之路相對清楚，感謝陪伴我度過研究所歲月的所有學長、 同學與學弟，有你們一起打拼，使這段歲月更加充實。 也要感謝父母與家人對我的支持與鼓勵，使我能一路走到現在，最後更 要感謝每一個幫助過我的人，有你們的幫助才能完成。

IV

目錄

目錄

目錄

目錄

摘 摘 摘 摘 要要要要 ... I ABSTRACT ... II 誌謝 誌謝 誌謝 誌謝 ... III 目錄 目錄 目錄 目錄 ... IV 表目錄 表目錄 表目錄 表目錄 ... VIII 圖目錄 圖目錄 圖目錄 圖目錄 ... IX 符號表 符號表 符號表 符號表 ... XIV 第一章 第一章 第一章 第一章 緒論緒論緒論 ... 1 緒論 1.1. 前言前言前言前言 ... 1 1.2. 探討風機效能依據探討風機效能依據探討風機效能依據探討風機效能依據 ... 3 1.2.1. 功率係數Cp ... 3

1.2.2. 端速比(Tip Speed Ratio, TSR) ... 3

1.3. 文獻回顧文獻回顧文獻回顧文獻回顧 ... 4

1.3.1. 水平軸風力發電機的研究 ... 4

1.3.2. 垂直軸升力式風力發電機的研究 ... 6

1.3.3. 垂直軸阻力式風力發電機的研究 ... 8

V 1.4. 研究內容研究內容研究內容研究內容 ... 15 第二章 第二章 第二章 第二章 數學模式數學模式數學模式 ... 17 數學模式 2.1. 基本假設基本假設基本假設基本假設 ... 17 2.2. 統御方程式統御方程式統御方程式統御方程式 ... 17 2.2.1. 旋轉區 ... 17 2.2.2. 靜止區 ... 17 2.2.3. RNG k-
Model ... 18 2.2.3.1. 紊流黏滯係數t ... 19

2.2.3.2. Inverse Effective Prandtl Numbers, α ... 19

2.2.3.3. 項 ... 19 2.3. 邊界條件邊界條件邊界條件邊界條件 ... 20 2.3.1. 進口邊界條件 ... 20 2.3.2. 出口邊界條件(outflow) ... 21 2.3.3. 對稱邊界條件(Symmetry Condition) ... 21 2.3.4. 壁面邊界條件 ... 21 第三章 第三章 第三章 第三章 數值方法數值方法數值方法 ... 22 數值方法 3.1. 簡介簡介簡介簡介 ... 22 3.2. 傳輸方程式離散化傳輸方程式離散化傳輸方程式離散化傳輸方程式離散化 ... 22 3.2.1. 非穩態項(Unsteady Term) ... 23 3.2.2. 對流項(Convection Term) ... 23

VI 3.2.3. 擴散項(Diffusion Term) ... 24 3.2.4. 網格梯度計算 ... 24 3.2.5. 源項(Source Term) ... 25 3.2.6. 線性代數方程式 ... 25 3.3. 壓力與速度之耦合關係壓力與速度之耦合關係壓力與速度之耦合關係壓力與速度之耦合關係 ... 25 3.3.1. 面上的質量流率 ... 25 3.3.2. 修正壓力方程式 ... 26 3.4. 解題步驟解題步驟解題步驟解題步驟 ... 27 第四章 第四章 第四章 第四章 實驗設備與方法實驗設備與方法實驗設備與方法 ... 28 實驗設備與方法 4.1. 實驗設備實驗設備實驗設備實驗設備 ... 28 4.1.1. 工業用風扇 ... 28 4.1.1. 風管 ... 28 4.1.2. 風速計 ... 29 4.1.3. 單層垂直軸風機之葉輪 ... 29 4.1.4. 扭力計 ... 29 4.1.5. 訊號擷取器 ... 30 4.1.6. 轉速計 ... 30 4.2. 實驗方法與步驟實驗方法與步驟實驗方法與步驟實驗方法與步驟 ... 30 第五章 第五章 第五章 第五章 結果與討論結果與討論結果與討論 ... 32 結果與討論

VII 5.1. 單層風機單層風機單層風機單層風機 ... 32 5.1.1. 網格數的影響 ... 32 5.1.2. 流場結構 ... 33 5.1.3. 轉速對扭矩的影響 ... 34 5.1.4. 入口風速對扭矩的影響 ... 34 5.1.5. 葉輪直徑D對扭矩的影響 ... 34 5.1.6. 平均扭矩係數 與端速比TSR之關係 ... 35 5.1.7. 功率係數Cp 與端速比TSR之關係... 35 5.1.8. 3D與2D暫態模擬之比較 ... 36 5.2. 實際實驗實際實驗實際實驗實際實驗 ... 38 5.2.1. 風速量測 ... 38 5.2.2. 模擬與實驗結果比較 ... 39 5.3. 多層風機多層風機多層風機多層風機 ... 40 5.3.1. 個別階層與總和之扭矩 ... 40 5.3.2. 風機層數與性能之比較 ... 42 5.3.3. 流場結構 ... 43 第六章 第六章 第六章 第六章 結論結論結論 ... 45 結論 參考文獻 參考文獻 參考文獻 參考文獻 ... 47

VIII

表目錄

表目錄

表目錄

表目錄

表 5. 1 不同實驗條件之平均入口風速 ... 39 表 5. 2 單層風機各係數比較表 ... 40 表 5. 3 不同層數風機的最大 Cp 值與最佳 TSR 值 ... 42

IX

圖目錄

圖目錄

圖目錄

圖目錄

圖 1. 1 水平軸升力型風車【3】 ... 53 圖 1. 2 水平軸阻力型風車【4】 ... 53 圖 1. 3 垂直軸阻力型(Savonius)風機【5】 ... 54 圖 1. 4 垂直軸升力型(Darrieus)風機【6】 ... 54 圖 1. 5 垂直軸升力型(H 型)風機【6】 ... 55 圖 1. 6 Darrieus 和 Savonius 混合型風機【7】 ... 55 圖 1. 7 改良型垂直軸風機 ... 56 圖 1. 8 垂直軸多層風機 ... 56 圖 1. 9 各種風機之 Cp 與 TSR 關係圖【9】 ... 57 圖 1. 10 斜葉片風機【26】 ... 57 圖 1. 11 上游方向角示意圖【27】 ... 58 圖 1. 12 Savonius 風機葉片重疊比例示意圖【30】 ... 58

圖 1. 13 Zephyr Vertical axis Wind Turbine (ZVWT) 【32】 ... 59

圖 1. 14 設有導流管之 Savonius 風機示意圖【34】 ... 59

圖 1. 15 螺旋式 Savonius 風機【36】 ... 60

圖 1. 16 改良式垂直軸風機示意圖【37】 ... 60

X 圖 1. 18 Savonius 風機前方設有一組擋板設計圖【41】 ... 61 圖 1. 19 Savonius 風機前方放置一圓弧型擋板示意圖【43】 ... 62 圖 1. 20 Savonius 風機前方裝設兩塊板子之示意圖【44】 ... 62 圖 1. 21 Savonius 風機在 Darrieus 風機上方之複合風機【45】 ... 63 圖 1. 22 Savonius 風機在 Darrieus 風機內部之複合風機【46】 ... 63 圖 1. 23 複合風機之 Cm-TSR 曲線【46】 ... 63 圖 1. 24 2D 擬暫態、暫態與實驗之 Cm-TSR 關係圖【47】 ... 64 圖 2. 1 流場示意圖 ... 64 圖 3. 1 ∆ri定義 ... 65 圖 3. 2 位置向量示意圖 ... 65 圖 3. 3 解題流程圖 ... 66 圖 4. 1 工業用風扇 ... 67 圖 4. 2 風管 ... 67 圖 4. 3 風管尺寸 ... 68 圖 4. 4 Lutron AM-4200 風速計 ... 68 圖 4. 5 單層風機葉片 ... 69 圖 4. 6 HBM T22 扭力計(50 N-m) ... 69 圖 4. 7 兩種不同尺寸之 coupling ... 70

XI 圖 4. 8 Rinstrum R320 訊號擷取器 ... 70 圖 4. 9 OMEGA HHT13 數位式轉速計 ... 71 圖 4. 10 無風管實驗配置示意圖 ... 71 圖 4. 11 有風管實驗配置示意圖 ... 72 圖 5. 1 XY 平面計算區域示意圖 ... 72 圖 5. 2 XZ 平面計算區域示意圖 ... 73 圖 5. 3 U0 = 8 m/s， N  60rpm，扭矩與時間關係圖 ... 73 圖 5. 4 U0 = 8 m/s， N  60rpm，扭矩與網格數關係圖 ... 74 圖 5. 5 XY 平面之計算網格 ... 74 圖 5. 6 XZ 平面之計算網格 ... 75 圖 5. 7 U0 = 8 m/s， N  120rpm，Cm 與時間關係圖 ... 75 圖 5. 8 在 10~10.5 秒，Cm 為極值時 3D 模擬的流線圖 ... 76 圖 5. 9 入口風速 U0= 8, 16m⁄s 時，扭矩與風機轉速 N 之關係圖 ... 77 圖 5. 10 風機轉速 N= 60, 120rpm 時，扭矩與入口風速 U0之關係圖 ... 77 圖 5. 11 入口風速 U0 = 8 m/s，風機轉速 N = 90, 120 rpm 時，扭矩與葉輪直 徑 D 之關係圖 ... 78 圖 5. 12 Cm-TSR 關係圖 ... 78 圖 5. 13 Cp-TSR 關係圖... 79 圖 5. 14 在 23~23.5 秒，Cm 為極值時 2D 模擬的流線圖 ... 80

XII 圖 5. 15 3D 與 2D 模擬結果在入口風速 U0 = 8m/s，轉速 N = 120rpm 下單一 週期之 Cm 與時間關係圖 ... 81 圖 5. 16 3D 與 2D 模擬結果在 U0 = 8m/s 下，扭矩與轉速 N 之關係圖... 81 圖 5. 17 3D 與 2D 模擬結果的 Cm-TSR 關係圖 ... 82 圖 5. 18 3D 與 2D 模擬結果的 Cp-TSR 關係圖... 82 圖 5. 19 風速量測示意圖 ... 83 圖 5. 20 0.5HP 風扇距離風機 0.7 公尺情況下的等風速線圖 ... 83 圖 5. 21 0.5HP 風扇和風管在距離風機 0.7 公尺情況下的等風速線圖 ... 84 圖 5. 22 模擬與實驗結果的 Cm-TSR 關係圖 ... 84 圖 5. 23 模擬與實驗結果的 Cp-TSR 關係圖 ... 85 圖 5. 24 二層風機在轉速 N = 60 rpm，入口風速 U0 = 8 m/s 的扭矩與時間關 係圖 ... 85 圖 5. 25 三層風機在轉速 N = 60 rpm，入口風速 U0 = 8 m/s 的扭矩與時間關 係圖 ... 86 圖 5. 26 四層風機在轉速 N = 60 rpm，入口風速 U0 = 8 m/s 的扭矩與時間關 係圖 ... 86 圖 5. 27 各種風機在轉速 N = 60 rpm，入口風速 U0 = 8 m/s 的總和扭矩與時 間關係比較 ... 87 圖 5. 28 三層與四層風機轉速 N = 60 rpm，入口風速 U0 = 8 m/s 的總和扭矩

XIII 與時間關係比較 ... 87 圖 5. 29 單層與多層風機的 Cm-TSR 比較圖 ... 88 圖 5. 30 單層與多層風機的 Cp-TSR 比較圖 ... 88 圖 5. 31 在 10~10.5 秒，Cm 為極值時，二層風機第一階的流線圖... 89 圖 5. 32 在 10~10.5 秒，Cm 為極值時，二層風機第二階的流線圖... 90 圖 5. 33 在 10~10.5 秒，Cm 為極值時，三層風機第一階的流線圖... 91 圖 5. 34 在 10~10.5 秒，Cm 為極值時，三層風機第二階的流線圖... 92 圖 5. 35 在 10~10.5 秒，Cm 為極值時，三層風機第三階的流線圖... 93 圖 5. 36 在 10~10.5 秒，Cm 為極值時，四層風機第一階的流線圖... 95 圖 5. 37 在 10~10.5 秒，Cm 為極值時，四層風機第二階的流線圖... 96 圖 5. 38 在 10~10.5 秒，Cm 為極值時，四層風機第三階的流線圖... 98 圖 5. 39 在 10~10.5 秒，Cm 為極值時，四層風機第四階的流線圖... 99

XIV

符號

符號

符號

符號表

表

表

表

符號 符號 符號 符號 定義定義定義定義 ,  關係係數 A 風機受風面積

α 鬆弛因子、Inverse effective Prandtl number

Cm 扭矩係數  平均扭矩係數 Cp 功率係數 | 最大功率係數 D 葉輪直徑  紊流生成項 H 葉輪高度 I 紊流強度 J 係數矩陣 k 紊流動能  ! 質量流率 N 風機轉速 P 風機功率

XV p 壓力 "# 位置向量 S 源項 $%& 平均應變率張量 TSR 端速比 '$|()*%+,+ 最佳端速比 U0 入口風速 Uin 接觸風機葉片前風速 Uout 流出風機葉片後風速 Utip 葉片尖端速度 -.# 流速 -.#/ 網格移動速度 W 權重函數 Γ 擴散係數 ε 紊流消散率 μ 層流黏滯係數 3 紊流黏滯係數 4!! 有效黏滯係數 ρ 空氣密度

XVI τ 風機扭矩 ϕ 變數 ω ..# 風機角速度 上標 上標 上標 上標 n 此一 time step 的已知量 n+1 下一 time step 的未知量 ‘, * 修正量 下標 下標 下標 下標 i, j 表 x, y 方向 f 面上的值 c0, p 主格點 ci, nb 相鄰格點

1

第一章

第一章

第一章

第一章 緒論

緒論

緒論

緒論

1.1.

前言

前言

前言

前言

人類使用風能已有相當長一段時間，好比在使用石化燃料以前，人們使 用風帆推動船隻、利用風車研磨穀物或汲水，然而，隨著石化時代的來臨， 風能也漸漸淡出人類的生活。但在電力剛踏入人類生活的早期，風能也不是 完全沒有貢獻，丹麥科學家 Poul la Cour 據信是第一個將風車與發電機連接 藉以發電的人【1】，但隨著柴油引擎與蒸氣渦輪機的發明，風能便逐漸退出 發電的舞台。然而，隨著能源短缺與環保意識高漲，風能作為低污染的再生 能源再次獲得重視。 風能來自於太陽能，雖由太陽輻射能轉換成風能的比例極低，但估算全 球風能總量相當是全球總發電量的 8 倍，約為地表可開發之水能資源的 10 倍，而目前已開發的僅是其中一小部份，若 1%的風能可供利用，則可減少 全球 3%的能源消耗，產生全球 8~9%的總發電量。【2】 傳統風力發電機有依旋轉軸方向與空氣作用力兩種分類法。依旋轉軸方 向可分為水平式(horizontal axis wind turbine, HAWT)與垂直式(vertical axis

wind turbine, VAWT)兩種。而依照空氣作用力可分為阻力型與升力型兩類。

由以上兩種分類方式，風機大致可分為以下幾類：

2 常見的風力發電機(如圖 1.1)【3】即屬此型，葉片數通常為 2 至 3 片，利 用升力使葉片轉動對轉子造成扭矩，當靜止時，葉片受風吹而產生阻力即可 自發性的轉動，因此不需額外提供能量便可自行啟動，故一般商用型風力發 電機以此型為主。

2.

水平軸阻力型風機

形如圖 1.2【4】，典型例子為古荷蘭風車，以農業用途居多，葉片數較多， 主要靠流體的阻力轉動，效率較低，轉速也較慢，適用於高扭矩之工作。

3.

垂直軸阻力型風機

一般稱為 Savonius 風機(如圖 1.3)【5】，葉片設計為 S 形，當氣流衝擊葉 片時，因兩邊形狀差異產生壓力差，進而產生旋轉力矩。效率雖較差，但任 何風向皆可運轉。此外，由於風機表面積較大，亦有在表面裝設太陽能板藉 此增加發電量的設計。

4.

垂直軸升力型風機

圖 1.4【6】的 Darrieus 風機為此類代表，而圖 1.5【6】的 H 型風機亦屬 此類。葉片截面類似機翼之流線型，當氣體流經葉片表面時，因兩翼面間的 壓力差而產生升力，進而產生力矩，使用機轉動。在高轉速下，因葉片角度 固定，會產生動力失速(stall)的問題，此時運轉效率將降低，此外由於 Darrieus 風機完全靠升力轉動，故無法自行啟動，但此問題可藉由與 Savonius 風機組 合改善(如圖 1.6【7】)。

3 本文係探討依改良式垂直式風力發電機(圖 1.7)，結合阻力型與升力型風 機特點，且可將數個葉輪堆疊(如圖 1.8)來提昇啟動效能。風機內有中央 3 個 S 型葉片，作用類似阻力型風機，而外側有兩個似翼形結構，在不同風向下 有升力效應產生，類似升力型風機，藉此提高效率及自啟動性。

1.2.

探討風機效能依據

探討風機效能依據

探討風機效能依據

探討風機效能依據

1.2.1. 功率係數 Cp 1919 年由 Betz【3】提出，定義如下： Cp  ₁ ; 2 <=>?@ 其中： ; 風機功率，;  AB，(A為風機扭矩，B為角速度) < 空氣密度 => 風速 @ 風機受風面積 Cp 值越大表示風機擷取越多風能。根據貝茲極限(Betz limit)，當接觸風機葉 片前風速=_%C為流出風機葉片後風速=_(,*的 3 倍(DEF DGHI  3)時，Cp 有最大值 59%。

1.2.2. 端速比(Tip Speed Ratio, TSR)

假設接觸風機葉片前風速為=_>，葉片尖端速度為=_*%)，則端速比定義為 TSR  DIEN

4

功率係數 Cp 與端速比 TSR 為探討風機效率的重要依據。圖 1.9【9】 為各型風機之端速比 TSR 與功率係數 Cp 關係圖，由此圖可知，當一理想 風機 Cp 值高時，必具有高 TSR 值。升力型風機(High speed two-blade, Modern

three blade, Darrieus rotor 等)便具有高 TSR 值，高 Cp 值之特性，工作範圍

約在 TSR = 3 ~ 7 之間，最佳 Cp 值則在 35%到 45%不等，適用於一般商業 發電。阻力型風機(American multi-blade, Savonius rotor 等)則是轉速低，Cp 值亦小， TSR 值多在 2 以下，最佳 Cp 值也不超出 30%多少，由此可知阻 力型風機較適用於高扭矩之工作。

1.3.

文獻回顧

文獻回顧

文獻回顧

文獻回顧

1.3.1. 水平軸風力發電機的研究

水平軸風機性能預測最常見的簡化方法是 Blade Element Momentum

method (BEM) 【10】，假設 2D 流場，由已知空氣動力學參數，升力係數_P與 阻力係數_Q求得外力，藉此求解動量方程式，因此需要準確的空氣動力學參 數才能得到較準確的結果。 而為探討葉片末端所產生的尾流，進而發展出 actuator disc【11, 12】及 actuator line【13, 14】等方法，前者於葉片後方流場建立圓盤型計算區，流經 旋轉中的葉片因而受影響的流體將通過此計算區，藉此探討其性質，後者則 是使用螺線作為計算區域，以減少計算量。

5

但若要探討葉片幾何尺寸的影響，仍需要 3D 計算方法【15, 16】，才有 能力獲得較精確的結果。

2007 年，Troldborg et al.【17】利用 large-eddy simulations(LES)以及 actuator

line 法模擬入口為紊流之流場，並與入口為均勻流的結果在平均風速相同下

比較，發現流場中尾流的大幅度運動主要由入口流中較大的渦流控制。

2008 年，Kasmi 與 Masson【18】在標準 k-ε model 的方程式中增加一能

量傳遞項，藉此解決標準 k-ε model 紊流動能過高的問題，並利用三葉片之

Nibe-B 630-kW、 Danwin 180-kW 與二葉片之 NASA/DOE Mod-0A 100-kW

風機驗證，結果顯示作者提供的方法較標準 k-ε model 準確。

2011 年，Fernando et al.【19】針對大氣邊界層(atmospheric boundary layer)

及其對風機與流場的影響進行模擬，為了提高空間與時間的解析度，使用

large-eddy simulation(LES)作為紊流模型，由於 LES 在大氣邊界層的精度與次

網格(subgrid-scale, SGS)紊流通量及風機引起的力有關，因此使用 tuning-free

Lagrangian scale-dependent dynamic models 處理 SGS 通量，而風機引起的力，

例如升力與阻力，則是使用 actuator-disk models 與 actuator-line models 處理。 將模擬結果與在風洞中對位於邊界層流中的三葉片水平軸風機的縮小模型所 測得的結果比較發現兩者的特性頗為相似。

水平軸升力型風機常利用扭曲葉片提昇效能，但由於製作困難度同時提 高，使得扭曲葉片的製造成本較高，因此在 2008 年，Lanzafame 與 Messina

6 【20】提出以分為兩個不同 pitch angle 區塊的風機葉片，兩個區塊之間則以 翼尖小翼(winglet)連接，希望達成與扭曲葉片相似的性能提昇效果與較低的 製造成本。為了最佳化此種葉片的配置，作者提出一由 BEM 衍生的數值方 法，並將模擬結果與實驗比較，結果顯示，此種葉片雖因葉尖漩渦而有部份 能量損失，但其效能仍較一般無扭曲葉片的風機高。 2007 年，Wang et al.【21】將一小型風機設置在導管中，藉以提昇低風 速下的風機效率，使用計算流體力學進行評估並在風洞中實驗，藉以驗證數 值模擬的準確性，由模擬結果可知，在存在導管的條件下，風速較沒有導管 的情況增加了 1.5 倍，功率輸出亦增加了 2.2 倍，而比較功率輸出曲線可發現， 數值模擬與風洞實驗的結果相當接近。 1.3.2. 垂直軸升力式風力發電機的研究

2010 年，Mukinović et al.【22】針對雙直 NACA-4418 葉片之 H 型風機

進行 2D 與 3D 數值模擬，並利用該結果改善 Double Multiple Stream Tube

(DMST) model。由 ANSYS-CFX 使用 k-ω SST model 的計算結果與 DMST model 的比較可以發現在端速比小於 2 下，兩者的 Cp 相當接近，但當端速比

較大時便有較大的差異，此差異主要源於 DMST 沒有考慮動態失速效應

(dynamic stall effect)，比較兩者的正規化扭矩與旋轉角度之關係可發現，在接

近0°時，計算流體力學預測值有較大的負向扭矩，兩者差異源於 DMST 使用

7

model 中。此外，作者亦使用 DLR 研發之 TAU，使用 Spalart-Almaras turbulence

model 進行 2D 與 3D 模擬並比較，由兩者之扭矩與旋轉角度關係可發現，主

要差異為扭矩最大值之大小，2D 計算得到的扭矩最大值較 3D 為大，由 3D 壓力差的不連續性分佈可發現，造成此結果的主因為有限的葉片高度(span)。

2010 年，Gupta 與 Biswas【23】利用移動網格與標準 k-ε model 在 2D

穩態下模擬一扭曲的三葉片 H-Darrieus 風機，並進行實際實驗，藉比較升力 和阻力係數與攻角之關係式驗證此一計算流體力學模型，並利用此一有效之 計算流體力學模型求取一最佳扭曲角度並評估其自啟動能力。 2011 年，Yamada et al.【24】針對一 H 型風機進行實驗，研究葉片厚度 與曲度對風機之瞬時扭矩與功率的影響，結果顯示在端速比相對較小的情況 下，具較小曲度與較大厚度的葉片之平均扭矩與功率增加較多，此外厚度對 對稱葉片的影響較非對稱葉片小。 2011 年，Kjellin et al.【25】在兩組固定轉速下改變風速對一 H 型風機進 行實驗，實驗結果顯示 Cp 最大值 0.29 發生在端速比為 3.3 的情況下，此外， 由於 Cp 曲線在接近最大值的附近相當平緩，將使系統易於達到最大 Cp 值。 2011 年，Armstrong 與 Tullis【26】針對螺旋葉片之高成本缺點提出以斜 葉片取代的設計(如圖 1.10【26】)，斜葉片具有螺旋葉片之減少負載變動的 優點以及可與直葉片匹敵的成本之優勢。藉由風洞實驗得知，斜葉片與直葉

8

片相似，皆對角度敏感(pitch-sensitive)，並有最佳角度S1.5°與S5.5°。在葉片 速度比─也就是最大端速比─為 2.1 時，有最佳 Cp 值為 0.27。

2012 年，Greenblatt et al.【27】在 H 型風機葉片之前緣裝設電漿介質放

電制動器(dielectric barrier discharge plasma actuator)作為控制流體分離之用， 而此制動器僅作用於上游一半方向角(如圖 1.11 【27】)。經由低速風洞實驗 得知，扣除制動器所需能量後，Cp 可提昇至 38%，此外，根據既有資料推論 當風機放大為 5 倍時，電漿所需能量將佔輸出功率之 3.3%，放大為 10 倍時 則是 1.7%。 2009 年，Deglaire et al.【28】提出適用於任意形狀與運動的垂直軸升力 式風機單一葉片的 2D 暫態解析模型，可用於描述葉片附近的非黏性壓力分 佈及作用於葉片上的外力，此模型係由共形映射(conformal mapping

techniques)與洛朗級數分解(Laurent’s series decomposition)所衍生，分析結果

與實驗頗為吻合，在相同精準度的條件下，計算時間較早期的分析模型為短。

1.3.3. 垂直軸阻力式風力發電機的研究

2011 年，Dobrev 與 Massouh【29】利用 k-ω model 與 DES(detached eddy simulation)針對 Savonius 風機進行 2D 與 3D 模擬，並與 PIV 實驗所得真實流

場比較，結果顯示 3D DES/k-ω model 結果最接近實驗數據。

2012 年，Akwa et al.【30】研究 Savonius 風機葉片重疊比例(如圖 1.12

9

Eddy Viscosity Turbulence Model k-ω SST，進行 2D 模擬，定義重疊比例為重

疊長度與葉片直徑之比例。結果顯示最佳重疊比例為 0.15，此時 Cp 最大值 為 0.3161。

2011 年，Shigetomi et al.【31】利用 PIV 對同一平面上的兩個 Savonius

風機之不同配置進行研究，發現兩種增加功率的交互作用，一是彎曲主流方 向使下游風機增加轉動，另一是由渦流脫離與週期性的壓力變化所造成的局 部流的週期性耦合。

2011 年，Pope et al.【32】使用套裝軟體 FLUENT 6.3.26，在 2D 穩態與

標準 k-ε model 及 MRF 模擬一 Savonius 變化型風機(如圖 1.13【32】)之效能， 並使用此數據及無因次參數提供預測此風機在幾何尺寸與操作參數變化後對 應 Cp 的方法，結果顯示預測值與實際進行數值計算之結果相當接近。

2010 年，D’Alessandro et al.【33】針對一葉片部份重疊之 Savonius 風機

進行研究，利用流固耦合探討流場與葉片的交互作用，流場使用 Direct

Numerical Simulation (DNS)以及 RANS 求解，葉片則是視為剛體，使用 Second

Cardinal Equation of Dynamic 求解。結果顯示數值計算所得的性能與相同狀

況下之實驗結果有相當的一致性。

2007 年，Irabu 與 Roy【34】為降低逆向扭矩並提高功率，在 Savonius

風機外側附加一導流管(如圖 1.14【34】)。結果顯示，在有導流管情況下， 二葉片風機可提昇 1.23 倍輸出功率，三葉片風機則是 1.5 倍。

10 2008 年，Saha et al.【35】利用風洞實驗驗證風機層數(1~3 層)、葉片數 量(2 或 3 片)與葉片扭曲與否(扭曲葉片與半圓葉片)對 Savonius 風機的影響， 另針對二層三葉片的 Savonius 風機進行單向閥設置與研究。結果顯示，不論 風機層數多寡及葉片形狀，二葉片的最大 Cp 值均較三葉片為大，這主要是 因為當有三片葉片時，流體由前方葉片反射至後方葉片，使後方葉片產生逆 向扭矩，因而降低其輸出功率；葉片形狀方面則是扭曲葉片較半圓葉片具有 較高的最大 Cp 值，這樣的現象是由合力作用點的位置不同所致，對於半圓 葉片，合力與葉片中心垂直，而對於扭曲葉片，施力點較靠近葉片尖端，力 臂較大，因此扭矩與輸出功率也較大；二層風機較其餘兩種具有較高的最大 Cp 值；因為單向閥可以有效的降低返回葉片所受的逆向扭矩，因此有單向閥 的二層三葉片風機較沒有的組合有較佳的效能。

2009 年，Zhao et al.【36】探討一具有螺旋狀扭曲葉片之 Savonius 風機(如

圖 1.15【36】)，並對不同幾何尺寸進行試驗。運用線性 k-ε model 進行 3D 模擬。結果顯示，對此設計而言，二葉片之 Cp 值較三葉片為佳，且對於靜 止之風機在任意角度其靜扭矩(static torque)為正值。

2009 年，Kamoji et al.【37】利用風洞測試一改良式 Savonius 風機(如圖

1.16【37】)，探討有無中央軸、葉片重疊比例、葉片周角大小、展弦比與雷

諾數對風機效能的影響。結果顯示，在實驗組合中，無中心軸 Savonius 風機 具有最高的 Cp 值。

11

而由於傳統的 Savonius 風機的靜扭矩(static torque)過大，以致在特定角 度無法自行啟動，因此 Kamoji et al.【38】同年(2009)亦提出以扭曲葉片降低 Savonius 風機之靜扭矩，在風洞中針對旋轉 90 度的扭曲葉片之 Savonius 風 機進行靜扭矩係數、扭矩係數與功率係數的量測，此外，並對中央軸的有無 與葉片重疊比例(0.0、0.1 與 0.16)進行探討。結果顯示，具扭曲葉片的 Savonius 風機在各個角度之靜扭矩係數皆為正值；在無中央軸，葉片重疊比例為零的 扭曲葉片 Savonius 風機具有最大功率係數 0.174，並與傳統 Savonius 風機的 功率係數近乎相同；而當中央軸存在時，其功率係數較同條件下的其餘組合 為小；此外，扭曲葉片的 Savonius 風機對雷諾數敏感，當雷諾數增加時，最 大功率係數亦增加。

2010 年，Mohamed et al.【39】探討在二葉片與三葉片 Savonius 風機前

方放置一擋板降低逆向扭矩的效果(如圖 1.17【39】)。此研究使用套裝軟體

ANSYS-Fluent 進行模擬，紊流模型為 realizable k-ε model。結果顯示，二葉

片風機 Cp 值提昇至 27.3%，三葉片風機則是 27.5%，但考慮造成本與提昇效 果，二葉片風機較理想。

據此，Mohamed et al.【40】於 2011 年，利用相同的方法進行擋板之大 小和位置與葉片幾何形狀之最佳化，並對效能提昇與自啟動能力進行研究。 結果顯示，在端速比為 0.7 時較標準 Savonius 風機有 38.9%的 Cp 相對提昇，

12

在整個操作範圍則至少有 30%的增加，而在各角度之靜扭矩(static torque)皆 為正，且足以自啟動。

2008 年，Altan 與 Atilgan【41】討論 Savonius 風機前方放有兩片擋板的

情況(如圖 1.18【41】)，進行實際實驗並模擬在固定角度(45°、60°、90°)下 靜止風機之靜扭矩(static torque)。數值模擬採用 Fluent 6.0 作為模擬工具，使

用標準 k-ε model。結果顯示擋板有助於提昇扭矩，此外當擋板越長效果越好，

最大靜扭矩發生在風機角度為60°，擋板角度分別為α  15°、β  45°時，模

擬與實驗數據亦大致吻合。

為降低 Savonius 風機扭矩的震盪與提昇啟動效能，Hayashi et al.【42】 於 2005 年在風洞中進行各階層相位差 120 度之三層 Savonius 風機與導流板 的研究。結果顯示三層 Savonius 風機在旋轉一週中的靜扭矩(static torque)與 扭矩曲線皆較單層 Savonius 風機來的平順以及當有導流板存在時，啟動效能 與扭矩係數皆有提昇。

2009 年，Hu 和 Tong【43】在一四葉片 Savonius 風機前方設有一圓弧擋

板(如圖 1.19【43】)，藉以減少逆向扭矩。使用之紊流模型為 RNG k-ε model

並以 SIMPLE 法求解 2D 流場。結果顯示當擋板角度約為30°時將有最大扭矩。

2010 年，Mohamed 與 Thévenin【44】在 Savonius 風機前方裝設兩塊板

子(如圖 1.20【44】)，此兩塊板子分別為降低逆向扭矩的擋板與提昇扭矩的 導流板，藉以達到在都市等低風速區域所需的低啟動扭矩與高輸出功率，而

13

作者的主要針對此兩塊板子的形狀、位置和彎曲與否進行最佳化，以取得最 佳組合。利用數值模擬進行研究，使用 ANSYS Fluent 搭配 in-house

optimization libaray(OPAL)達到自動最佳化的效果，以 Cp 值作為評估最佳化 成果的唯一參數，使用評估方式則為 Evolutionary Algorithms，紊流模型則是 使用 realizable k S ε model。結果顯示，在有板子的情況下，其 Cp 值較沒有 的增加且可操作之 TSR 範圍亦增加，其中，最佳組合的 Savonius 風機的輸出 功率較無板子的增加了 48%。 1.3.4. 垂直軸複合式風力發電機的研究

2008 年，Gupta et al.【45】提出三葉片 Savonius 風機結合三葉片 Darrieus

風機的設計，其中，Savonius 風機在 Darrieus 風機上方，見圖 1.21【45】，藉 以得到較 Savonius 風機為高的效率與較 Darrieus 風機為佳的啟動效能，變數 為 Savonius 風機部份之葉片重疊比例。經由次音速風洞實驗可知，當葉片重 疊比例上升時，Cp 值將隨之下降，其中，最大 Cp 為 51%，此時葉片重疊比 例為 0。

2009 年，Elmabrok【46】針對結合 Darrieus 與 Savonius 風機的複合風機

進行理論研究，風機外型見圖 1.22【46】，此風機為 Savonius 風機在內部，

Darrieus 在外部，利用 multiple stream tube theory 預測 Darrieus 風機的扭矩係

數 Cm，又利用作者於 2003 年所提出的針對 Savonius 風機的新理論方法預測

14 此複合風機的扭矩係數 Cm。除理論研究外，作者又在次音速風洞中進行實 驗，藉以量測此風機的功率與扭矩。量測結果與理論解的比較見圖 1.23【46】， 實驗條件為葉片重疊比例為 0.25，Savonius 風機葉片為半圓形。由圖中可以 發現，在端速比較小時(TSR > 1.0)，主要由 Savonius 風機掌控，Cm(圖中表 示為 CQ)值隨 TSR 值增加而減少；當 TSR 值介在 1.0 到 2.5 時，則由 Darrieus 風機主導，Cm 值隨 TSR 值增加而增加；當 TSR > 2.5 時，兩者的影響力相 差不大，因此 Cm 值隨 TSR 增加而減少。此外，由圖中亦可發現理論解與實 驗結果有一定程度的相似性。

2011 年，林子翔【47】利用多重參考座標系統(Multiple Reference Frames,

MRF)處理葉片轉動，紊流模型則使用 RNG k-ε model，每隔 30°或 10°進行 一次計算，模擬在固定角度的一複合式風機之效能與 2D 流場，由於除葉片 角度為 0°時以外，扭矩係數 Cm 並不會達到一穩定值，故使用暫態的算法， 求取固定角度之 Cm 值隨時間變化關係，再取一段時間後的 Cm 加以平均， 視為該角度之扭矩係數，將此種計算方法稱之為擬暫態。 接著討論風機轉速、入口風速與風機直徑對輸出扭矩的影響，結果顯示 扭矩與風機轉速呈一線性遞減關係，與入口風速呈二次式關係，與風機直徑 呈三次式關係。再將以上有因次參數無因次化，得到無因次參數 Cm 與端速 比 TSR 呈線性遞減關係，再將 Cm 轉換為功率係數 Cp，得到 Cp 與端速比呈 二次式關係。

15 最後與使用移動網格計算之暫態 2D 以及實驗結果比較，三種方法得到 的參數皆有相同的關係，但值各有不同，圖 1.24【47】為擬暫態(fixed rotor)、 暫態(rotating)與實驗(exp)之 Cm-TSR 關係圖，由此圖可發現趨勢相同而值不 同的現象，在相同端速比的條件下，Cm 值由大至小分別為暫態、實驗、擬 暫態，造成這些差異的原因有很多，使實驗結果較小的主因為實驗設備過於 簡陋，尤其是實驗時，葉輪相當容易晃動，使輸出扭矩不準確，而擬暫態主 要是因忽略流體慣性，所以輸出扭矩較小。因此認為 MRF 僅適用於定性分 析，在定量分析之能力上稍嫌不足。

1.4.

研究內容

研究內容

研究內容

研究內容

本研究使用計算流體力學(Computational Fluid Dynamics, CFD)針對單層 和多層風機進行 3D 流場之模擬與性能預測，模擬工具則是使用套裝軟體 ANSYS Fluent 13.0，並利用移動網格處理風機葉片的旋轉，將風機葉片旋轉 範圍內設為旋轉區，外圍則定為靜止區。紊流模型則是使用 RNG k-ε Model 配合標準壁函數模擬風力發電系統之流場。而本文將探討單層風機在不同風 速、轉速、葉輪直徑下產生的扭矩與其關係，再將這些有因次參數轉換為扭 矩係數 Cm、功率係數 Cp 與端速比 TSR 等無因次參數，並比較各無因次參 數間的關係，此外，另針對風機層數進行探討，比較單層與多層風機差異及 最大功率係數與風機層數的關係。

16 除以上 3D 流場模擬外，亦使用與風機中央平面相同的計算區域進行 2D 暫態計算，藉以比較 3D 與 2D 模擬之差異。並與林子翔【47】所作 2D 擬暫 態模擬結果比較，藉以探討暫態與擬暫態模擬之差異性。 此外，為了驗證數值模擬的準確性，針對單層風機進行實際實驗，而 由於缺少可進行全尺寸實驗的大型風洞，故採用以工業風扇作為風源的簡 化實驗，除工業風扇直接吹向葉輪的組合外，又有在工業風扇與葉輪間設 置一裝有蜂巢結構之風管的組合，在有風管的情況下，除了平均入口風速 有所改變外，入口風速的均勻度也隨之改善，將實驗所測得的風機扭矩與 轉速的關係無因次化後可得扭矩係數 Cm 與端速比 TSR 之關係，加以計算 後又可得功率係數 Cp 與端速比 TSR 的關係，將由實驗所得的此兩個關係 與模擬結果比較，藉以驗證模擬的準確性。

17

第二章

第二章

第二章

第二章 數學模式

數學模式

數學模式

數學模式

2.1.

基本假設

基本假設

基本假設

基本假設

工作流體為空氣，假設流場為等溫狀態且忽略重力的影響。

2.2.

統御方程式

統御方程式

統御方程式

統御方程式

因分析垂直式風力發電機時需將計算區域分為兩個區塊(如圖 2.1)，葉輪 掃過之旋轉區域與外圍靜止區域。旋轉區域以給定的角速度B..#旋轉，靜止區 則不旋轉，也就是B..#  0，以下分別介紹兩區域之統御方程式 2.2.1. 旋轉區 連續方程式 (2.1) 動量方程式 (2.2) 其中，-.#_/為網格移動速度，由給定的角速度B..#計算而得，-.#_/  B..# Z "#。 2.2.2. 靜止區 在靜止區中，由角速度B..#  0，可知網格移動速度-.#_/  0，因此靜止區之 統御方程式如下： 連續方程式 (2.3) 動量方程式 (2.4)

(

_g

)

0

V V

d

dV

u

u

dA

dt

∫

ρ

+

∫

∂

ρ

−

=

i

(

g

)

V V V V

d

udV

u u

u

dA

u dA

SdV

dt

∫

ρ

+

∫

∂

ρ

−

= Γ∇

∫

+

∫

i

i

0

V V

d

dV

u dA

dt

∫

ρ

+

∫

∂

ρ

=

i

V V V V

d

udV

uu dA

u dA

SdV

dt

∫

ρ

+

∫

∂

ρ

= Γ∇

∫

+

∫

i

i

18

2.2.3. RNG k-[ Model

RNG k-ε model 是利用 Renormalization group (RNG)方法由瞬時 Navier-Stokes 方程式推導得到，與標準 k-ε model 有以下幾個不同點：

在ε方程式中增加一_\項，藉此增加應變流的準確性。

對旋轉流體的準確性更佳。

使用公式求解 turbulent Prandtl numbers，而標準 k-ε Model 使用固定常數。
標準 k-ε model 僅適用於高雷諾數流場，而 RNG k-ε model 因計算有效 黏滯係數的公式包含低雷諾數效應，若近牆區域處理得當便可適用於低 雷諾數流場。 k 方程式 (2.5) ε方程式 (2.6) 為紊流生成項：   Sρ- `-]^-_^<sub>`b</sub>a] &  2cdd$%& `-a] `b& (2.7)

$%&代表平均應變率張量(mean rate-of-strain tensor)：

$%&  1_{2 e}`-<sub>`b</sub>a] & f `-a_ `b%g (2.8) h\與i\係利用 RNG 方法所求得之常數，分別為h\  1.42與i\  1.68， 而$與$_\則為源項。 ( ) ( _i) ( _eff ) _{k k} i j j k k ku G S t ρ x ρ x αµ x ρε ∂ ₊ ∂ ₌ ∂ ∂ _{+ − +} ∂ ∂ ∂ ∂ 2 1 2 ( ) ( _i) ( _eff ) _k i j j u C G C R S t x x x ε k ε k ε ε ε ε ε ρε ρε αµ ∂ ₊ ∂ ₌ ∂ ∂ _{+ − − +} ∂ ∂ ∂ ∂

19 2.2.3.1. 紊流黏滯係數₃ RNG k-ε model 計算紊流黏滯係數使用的公式為： k l<im √o  1.72 qr sqr? _{S 1 f t}kqr (2.9) 其中： qr uvww u 、t x 100 (2.9)式適用於低雷諾數與靠近壁面之流場，在高雷諾數流場之情形下， 此式為： 3  <um i  (2.10) 其中，µeff = +µ µt 。此外，以 RNG 方法可求得u  0.0845，與在標準 k-ε model 中由經驗設定的 0.09 頗為接近。

2.2.3.2. Inverse Effective Prandtl Numbers, α

α為 k 與

ε

項的 inverse effective Prandtl numbers，而 RNG k-ε model 計算 α的公式如下： y_zz S 1.3929 > S 1.3929y >.{?ih y_zz f 2.3929 >f 2.3929y >.?{|} _ 4!! (2.11) 其中z_>  1.0。若在高雷諾數流場中( u uvww ≪ 1)，α x 1.393。 2.2.3.3. _\項 此項為 RNG k-ε Model 與標準 k-ε Model 最大的差異，表示為：

20 \  u< ?_{€1 S  ⁄ ‚>} 1 f ƒ?  i m (2.12) 其中η  $m ⁄ 、_>  4.38、β  0.012、S  s2$_%&$_%&。 將(2.12)代入(2.6)可得ε方程式最後形式為： ∂ ∂t €<‚ f`b`%€<-%‚  ` `b&lz4!! ` `b&o f h\  m S i\∗ < i m (2.13) 而_i\∗  _i\ f‡ˆ‰Š€h‹‰ ‰⁄ ‚O hŒ‰Š (2.14)

2.3.

邊界條件

邊界條件

邊界條件

邊界條件

討論的流場為一開放區域，包含進出口、對稱邊界與固體壁面，其中固 體壁面為轉動中的葉輪構成之移動壁面。 2.3.1. 進口邊界條件 定義進口為一水平方向固定流速之均勻流，而進口之 k 與ε值定義如下： k 3_{2 €=}>Ž‚i (2.15) ε  ρum i   * ‹h (2.16) 其中： U> 進口速度 I 紊流強度(turbulence intensity)，假設為 10% *

21 2.3.2. 出口邊界條件(outflow) 利用外推(extrapolates)計算流場內部出口的性質，使各性質滿足在垂直出 口方向的梯度為零的條件，此條件是基於完全發展流而訂定的，若是該梯度 不可忽略的位置，例如有回流發生的位置，就不適合假設為出口。此外，另 針對質量守恆進行訂正。 2.3.3. 對稱邊界條件(Symmetry Condition) 為了減少計算範圍，使用對稱邊界條件(Symmetry Condition)，此類邊界 條件假設對稱邊界面上垂直方向速度為零，其餘流場性質則假設梯度為零。 2.3.4. 壁面邊界條件 由於無滑移邊界條件，故壁面上流體速度等於固體壁面移動速度，而

22

第三章

第三章

第三章

第三章 數值方法

數值方法

數值方法

數值方法

3.1.

簡介

簡介

簡介

簡介

本研究使用 ANSYS Fluent 13.0 作為模擬工具，並搭配 Gambit 建立網格。 為了取得較準確的結果，將計算區域分為葉輪掃過之旋轉區與此區以外的靜 止區兩區域，紊流模型則使用 RNG k-ε model，求解 3D 暫態流場。 本章節主要介紹 ANSYS Fluent 13.0 所使用的數值方法，3.2 節處理動量 方程式離散化，3.3 節為壓力與速度耦合關係式，3.4 節則介紹 SIMPLE 演算 法解題步驟。

3.2.

傳輸方程式離散化

傳輸方程式離散化

傳輸方程式離散化

傳輸方程式離散化

在數學模式章節中所提到的連續方程式、動量方程式及 k 和ε方程式無論 在旋轉區或是靜止區均可用下列通式表示： ’ ∂<“_{∂t k”} • f – <“€-.# S -.#/‚ ∙ k@#  – ؙ“ ∙ k@# f ’ $š ˜k” (3.1) 對於靜止區-.#_/  0。

由於使用 pressure-based solver，ANSYS FLUENT 使用隱式(implicit)離散， 故(3.1)式可改為： ’ ∂<“_{∂t k”} • f – < CŒh_“CŒh_›-.#CŒh_{S -.#/}CŒh_{œ ∙ k@#  – Γ} ˜CŒh™“CŒh∙ k@# f ’ $˜CŒhk” š (3.2) 使用有限體積法對(3.1)式進行離散化，以下分別對非穩態項、對流項、 擴散項以及源項分別詳述：

23 3.2.1. 非穩態項(Unsteady Term) ’ ∂<“_{∂t k”} •  ∂<“ ∂t ∆V (3.3) 而ž˜ ž3可表示為： ∂“ ∂t  F€“‚ (3.4)

使用 Second Order Implicit 對(3.4)式進行離散化： 3“CŒh_{S 4“}C_{f “}C‹h 2∆t  F€“CŒh‚ (3.5) “CŒh _{ ›4“}C_{S “}C‹h_{f 2∆ ¡€“‚œ 3⁄ (3.6)} “CŒh_、“C_與“C‹h_{分別為t  t f ∆t, t, t S ∆t時的“值，除了“}CŒh_為未知數 外，“C與“C‹h均已知。 3.2.2. 對流項(Convection Term) 對於對流項利用中點定理(Midpoint rule)將面積分轉為差分形式： – <“€-.# S -.#/‚ ∙ k@# x ¢ <!“!£-.#!S -.#_/!¤ ∙ @#! !  ¢  !“! ! (3.7)

其中 _!為質量流率(mass flow rate)。

對於空間採二階上風差分法(Second Order Upwind Scheme)處理，藉以確 保較高的精度：

“!  “ f ™“ ∙ "# (3.8)

其中：

24 ™“ 上游網格中心之梯度 "# 由上游網格中心至面上中點的位置向量 3.2.3. 擴散項(Diffusion Term) 與擴散項類似，利用中點定理將面積分轉為差分形式： – ؙ“ ∙ k@# x ¢›Γ˜™“œ_d ∙ @#! ! (3.9) 3.2.4. 網格梯度計算

使用 Least Squares Cell-Based 求取計算對流項與擴散項所需的™“，定義 為沿向量"_%的 c0 與 ci 網格間的“值差，數學型式如下： €™“‚_¥¦∙ ∆"§  “¨% S “¨> (3.10) ∆"§定義為網格 c0 中心至網格 ci 中心之距離(如圖 3.1)。 將所有相鄰點代入(3.10)式，則可整理為： © ª«€™“‚¥¦  Δ“ (3.11) © ª«為只與幾何形狀有關之係數矩陣。 使用 Least Squares 求解網格梯度(™“_>  “_­®̂ f “_°±̂ f “_²m³)即為利用權重 【48】分解矩陣，解法為將權重與差異向量(difference vector, Δ“  “_¨%S “_¨>) 相乘，形式如下： €“_­‚_¨>  ¢ ´­ %>∙ €“¨% S “¨>‚ C %µh (3.12) ›“°œ_¨>  ¢ ´°%>∙ €“¨% S “¨>‚ C %µh (3.13)

25 €“_²‚_¨>  ¢ ´² %>∙ €“¨% S “¨>‚ C %µh (3.14)

Least Squares Cell-Based 的優勢在於使用非結構網格時，其精度可與

node-based 的結果匹敵，兩者結果皆較 cell-based 好許多，但 Least Squares

Cell-Based 之計算量較 node-based 少。 3.2.5. 源項(Source Term) ’ $˜k” š x $˜∆” (3.15) 3.2.6. 線性代數方程式 將非穩態項、對流項、擴散項及源項合併與整理後可得離散化之傳輸代 數方程式： )“  ¢ C¶“C¶ C¶ f  (3.16) 其中下標 p 者為被計算點，nb 為相鄰網格。 為使疊代過程較穩定，引入一鬆弛因子α (under-relaxation factor)： “  “(·¸ f z∆“ (3.17) 則(3.15)式可修改為： )“ z  ¢ _C¶ C¶“C¶ f  f1 S zz )“(·¸ (3.18)

3.3.

壓力與速度之耦合關係

壓力與速度之耦合關係

壓力與速度之耦合關係

壓力與速度之耦合關係

3.3.1. 面上的質量流率

26

為滿足連續方程式，必須計算面上速度，ANSYS Fluent 引用 Rhie【49】 所提出的方法求取面上質量流率：  !  €<q....#@‚C !  ª!@! (3.19) 其中ª_!為通過面 f 的質量通量(mass flux)，計算法如下： ª!  <!),¨>q_C,¨> f ),¨hqC,¨h ),¨> f ),¨h f k!¹›¨>f €™‚¨>∙ "0....#œ S ›¨h f €™‚¨h ∙ "1....#œº  ª»! f k!€¨>S ¨h‚ (3.20) 其中_¨為壓力，q_C為垂直方向之速度，下標 c0、c1 為與面相鄰兩網格，r0 與 r1 則為網格中心至面中心之位置向量(如圖 3.2 所示)，k_!則為 之函數。 ₎ 3.3.2. 修正壓力方程式

在模擬過程中的每一個 time step 裡使用 SIMPLE 演算法進行計算，此演 算法利用速度與壓力間的關係強制質量守恆成立並求解壓力場。若利用一猜 測的壓力∗解動量方程式，則由(3.20)式所得的面上通量ª_!∗形式如下： ª!∗  ª»!∗ f k!€¨> S ¨h‚ (3.21) 但上式並不滿足連續方程式，因此需加入修正項ª_!^，修正後之面上通量形 式如下： ª!  ª!∗ f ª!^ (3.22) 使其滿足連續方程式，而 SIMPLE 演算法假定ª_!^形如： ª_!^ _{ k} !€¨>^ S ¨h^ ‚ (3.23)

27 p’為網格之壓力修正項。 SIMPLE 演算法將通量修正項((3.22)式至(3.23)式)代入離散之連續方程式， 藉此得到網格壓力修正項的離散方程式： )′  ¢ C¶C¶^ C¶ f  (3.24) 其中源項 b 為網格淨流率：   ¢ ª!∗@! ½w¾¿vÀ ! (3.25)

ANSYS Fluent 使用代數多重網格(Algebraic Multigrid, AMG)法求解壓力

修正方程式(3.24)式後可得壓力修正項 p’，用於修正網格壓力與面上通量：   ∗ _{f z} )′ (3.26) ª!  ª!∗ f k!€¨>^ S ¨h^ ‚ (3.27) 其中z₎為鬆弛因子，且修正後滿足連續方程式。

3.4.

解題步驟

解題步驟

解題步驟

解題步驟

解題步驟之流程圖如圖 3.3。 在實際計算上，若單一 time step 中計算次數達到某一限值(目前設為 20)， 但仍未滿足收斂條件，將直接進入下一個 time step，而每一 time step 之初始 值為上一 time step 的最終值。

28

第四章

第四章

第四章

第四章 實驗設備與方法

實驗設備與方法

實驗設備與方法

實驗設備與方法

本研究除利用數值模擬進行探討外，亦針對單層風機進行簡易實驗，藉 以驗證數值模擬的準確性。實驗部份由交通大學複合材料結構實驗室協助完 成，除扭力計、訊號擷取卡、轉速計與風管外之器材皆由該實驗室提供。本 章主要針對實驗設備與實驗步驟進行介紹，實驗結果將在第五章進行討論。

4.1.

實驗設備

實驗設備

實驗設備

實驗設備

4.1.1. 工業用風扇 本實驗為一開放流場，使用工業用風扇(圖 4.1)作為風源，此風扇為單相 6 極四葉片，功率為 0.5HP 的風扇，中心高 0.65m，葉片直徑約 0.59m，風速 分佈並不均勻，其中，中心與邊緣區域之風速較弱。 4.1.1. 風管 為處理風扇之不均勻的風速分佈，在風扇與風機之間放置一裝設有鋁製 蜂巢板的風管(圖 4.2)，藉以使風速分佈較為平均。風管尺寸如圖 4.3 所示， 風管截面積約為0.8m Z 0.8m，深度為 0.9m，為與風扇高度配合，風管下緣 高度為 0.52m。

29 4.1.2. 風速計 為了量測風機前的風速，使用 Lutron 公司所製造的 AM-4200 風速計(圖 4.4)進行量測，此風速計的量測範圍為0.8~30.0  Â⁄ ，精準度為Ã3%。 4.1.3. 單層垂直軸風機之葉輪 圖 4.5 為本研究所探討的垂直軸風機之葉輪的單層實體【50】，此葉輪由 三個 S 型葉片、兩個似翼形葉片、上頂板與下底板所組成，葉輪直徑 1m， 高度 0.28m。使用樹脂與玻璃纖維布之複合材料製成，此種材料的優點為抗 腐蝕、耐磨耗且重量輕，密度約為鐵的二分之一。 4.1.4. 扭力計 本研究使用德國 HBM 公司製造的 50 N-m 型 T22 扭力計(圖 4.6)量測扭 力值，此型號扭力計可量測範圍為Ã50 N-m，精度為 0.25 N-m。扭力計之兩 端分別連接測試端與驅動端，在本研究中分別為垂直軸風機之葉輪與發電機， 由於扭力計、葉輪與發電機之軸徑各不相同，故使用 coupling 連接(圖 4.7)。 扭力計所測得的扭力值經轉換後輸出直流電壓訊號，其值介於Ã5V之間，當 直流電壓訊號為正時，表示輸入扭力計的扭力為順時針方向，負值則表示逆 時針方向，5V 表示扭力計可量測範圍的最大值，在此型扭力計為 50 N-m。

30 4.1.5. 訊號擷取器 本研究所使用的訊號擷取器為 Rinstrum 公司製造的 R320 擷取器(圖 4.8)， 扭力計所輸出的訊號經由訊號線傳至訊號擷取器，將直流電壓訊號轉換為扭 力值，最後將扭力值顯示於液晶螢幕上。 4.1.6. 轉速計 本研究使用 OMEGA 公司所生產的 HHT13 轉速計(圖 4.9)，此型號的轉 速計具有非接觸量測的功能，原理係利用雷射反射，將反光貼紙貼在軸心上， 將轉速計發出紅光雷射照射到反光貼紙上，經反射後由轉速計接收並計算轉 速，此種測量方法的量測範圍為 5~200,000 RPM。

4.2.

實驗方法與步驟

實驗方法與步驟

實驗方法與步驟

實驗方法與步驟

在本研究中，風機性能試驗分為兩種(圖 4.10、圖 4.11)，兩種實驗方式使 用的器材與連接方式及實驗步驟大致相同，兩者間的最大差異在於風管的有 無。法一，如圖 4.10 所示，工業用風扇所吹出的風直接與風機接觸，而法二， 如圖 4.11 所示，在工業用風扇與風機間另設置風管，使風機迎風面上的風速 分佈較為均勻，而與法 B 比較，由於法 A 的工業用風扇直接與風機接觸，因 此法 A 在風機迎風面上的平均風速較大。除此之外，並與不同的風源與風機 間距離(0.6, 0.7, 0.8m)進行搭配，藉以獲得不同的入口風速。

31 其餘的實驗設備則是以相同的方式連接，葉輪的旋轉軸與扭力計連接， 扭力計的另一端與發電機連接，扭力計的訊號則是利用訊號線傳輸至訊號擷 取器，再將扭力值以液晶螢幕顯示，風機轉速則是使用轉速計以紅光雷射進 行非接觸式的量測。而發電機另與可變電阻連接，此可變電阻即為發電機的 負載，藉由改變負載，也就是電阻值的大小使風機轉速產生變化，藉以得到 在不同輸出扭矩值條件下的風機轉速。

32

第五章

第五章

第五章

第五章 結果與討論

結果與討論

結果與討論

結果與討論

考慮一個 3D 流場，如圖 5.1(XY 平面)、圖 5.2(XZ 平面)所示，有多層與 單層之分，每以葉輪包含三個 S 型葉片與兩個似翼形葉片。葉輪直徑 D 為 1 m， 不論層數多寡，葉輪高度 H 皆為 1 m，R 為葉輪直徑的 25 倍(R = 25D)，L 為 葉輪直徑的 50 倍(L = 50D)，旋轉區域為各方向的葉輪尺寸向外延伸 0.1 倍 (1.1D、1.1H)。2D 的計算範圍與 3D 流場中與葉輪相交的 XY 平面相同，也 就是圖 5.1 所表示的計算區域。 左側假設為一流速為 U0的均勻流入口，右側則假設為出口，其餘邊界皆 設為對稱邊界。

5.1.

單層風機

單層風機

單層風機

單層風機

5.1.1. 網格數的影響 圖 5.3 為入口風速 U0 = 8 m/s，風機轉速 N = 60 rpm 下，不同網格所得 的扭矩與時間在 9 秒到 10 秒間的關係圖，由圖中可以發現，四種網格所得 的扭矩與時間的曲線的趨勢相當相似，網格數為 238,180 與 406,794 的兩組 網格其扭矩較其餘兩組網格(網格數為 563,560 與 879,713)所求得的略小， 此現象同樣可在平均扭矩(自第五秒起累計平均)與網格數的關係圖(圖 5.4) 中發現，由圖中可以發現，在網格數從 238,180 增加到 563,560 時，平均扭 矩有明顯的增加，但當網格數由 563,560 增加到 879,713 時，兩者的平均扭

33 矩相差約 1.8%，變化量減少，而為了減少計算時間，採用網格數為 563,560 的網格。 因此，以下主要針對網格數為 563,560 的網格進行討論，此網格的 XY 平面之計算網格約 8040 個，形式如圖 5.5 所示，Z 方向分割為 70 格，XZ 方向的網格如圖 5.6 所示。 5.1.2. 流場結構 圖 5.7 為在入口風速 U0 = 8 m/s，風機轉速 N = 120 rpm 的情況下的扭矩 係數(Moment coefficient, Cm)與時間之關係圖，Cm 定義如下： 2 0 1 4 Cm U AD τ ρ = _(5.1) τ表示風機產生之扭矩(n-m)，A 代表風機受風面積(m2_{)，A = DH，D 為風機} 直徑(m)。每週期平均為針對每一週期(0.5 秒)的 Cm 個別平均，自第五秒末 累計平均則是從第五秒末(第 10 個週期)起，對 Cm 值累計平均。 圖 5.8 則為該條件下 10 到 10.5 秒間當 Cm 為極值時葉輪中央平面的流 線圖，從圖中可發現葉輪中心部份的流體隨著葉輪以順時針方向旋轉，自 左上緣流入葉輪內的流體則在右側下方流出，並與流過葉輪下方的流體混 和，進而形成小渦流，在整個週期中，此流場結構並不隨著葉片位置改變 而有太大的變化。

34 5.1.3. 轉速對扭矩的影響 圖 5.9 為入口風速 U0 = 8, 16 m/s，風機轉速 N = 15, 30, 60, 90, 120, 150, 180, 210 rpm 的條件下所產生的扭矩與風機轉速關係圖，由此圖可知不論入 口風速為何，扭矩與轉速幾乎呈一線性遞減關係。此外，在入口風速 U0 = 8 m/s，風機轉速 N 約為 120 rpm 時，扭矩近乎為零，這代表在入口風速 U0 = 8 m/s 時，此風機的自由旋轉速度約為 120 rpm。 5.1.4. 入口風速對扭矩的影響 圖 5.10 為風機轉速 N = 60, 120 rpm，入口風速 U0 = 4, 8, 12, 16, 18, 20 m/s 時的扭矩與入口風速的關係圖，由圖中可以發現在入口風速較小時，以 風機轉速 N = 60 rpm 的情況為例，約是入口風速 U0 < 16 m/s，扭矩隨入口 風速而增加，但增量隨之減少，且當 U0 > 16 m/s，當入口風速增加時，扭 矩反隨之下降，很明顯的，扭矩與速度呈非線性關係，之後將證明兩者呈 二次關係。 5.1.5. 葉輪直徑 D 對扭矩的影響 圖 5.11 為入口風速 U0 = 8 m/s，風機轉速 N = 90, 120 rpm，葉輪直徑為 0. 5, 1.0, 1.5, 2.0 m 的條件下所產生的扭矩與葉輪直徑關係圖，觀察此圖可 發現，在葉輪直徑大於等於 1.0 m 的情況下，隨著葉輪直徑增加，扭矩值隨 之下降，但兩者間並不呈線性關係，之後將證明兩者呈三次式關係。

35

5.1.6. 平均扭矩係數ÅÆ 與端速比 TSR 之關係

端速比(Tip Speed Ratio, TSR)為討論風機性能所需的一個重要無因次 參數，定義為葉片尖端速度與入口風速的比值： 0 60 N D TSR U π = (5.2) 將扭矩與風機轉速、入口風速及葉輪直徑的關係無因次化可得無因次化 的平均扭矩係數Cm與端速比 TSR 的關係，也就是圖 5.12，由此圖可發現平 均扭矩係數Cm與端速比 TSR 呈一線性遞減關係： Cm= −aTSR b+ (5.3) 上式中的 a, b 皆為正值，3D 模擬的結果可得 a = 0.572，b = 0.513，若將Cm 與 TSR 的定義帶入，則可得有因次參數間的關係： 2 2 0 0 240 4 a b U AD N U AD π τ = − ρ + ρ (5.4) 風機受風面積 A = DH，H 為風機高度，因此(5.4)式又可改寫為： 3 2 2 0 0 240 4 a b U D HN U D H π τ = − ρ + ρ (5.5) 由上式可證明前述現象：扭矩與風機轉速呈線性遞減關係、入口風速呈二次 式關係、與葉輪直徑呈三次式關係。 5.1.7. 功率係數 Cp 與端速比 TSR 之關係 功率係數(Power coefficient, Cp)為評估風機效能所需的一個重要無因次 參數，定義如下：

36 3 0 1 2 P Cp U A ρ = _(5.6) (5.6)式中，P 代表風機之功率，P=τω，ω=2π 60N，也就是風機的角速度， 比較 Cm 與 Cp 的定義可得兩者的關係： Cp =Cm TSR⋅ (5.7) 將Cm與 TSR 的關係，(5.3)式帶入(5.7)式便可得 Cp 與 TSR 的關係式：

(

)

2 Cp= − ⋅a TSR + ⋅b TSR (5.8) 由上式可知 Cp 與 TSR 呈二次式關係(見圖 5.13)，利用(5.8)式可求取最佳 TSR 之值： 2 optimum b TSR a = (5.9) 當 TSR 為最佳 TSR 值時，此時的 Cp 最大值為: 2 4 MAX b Cp a = (5.10) 將由模擬結果所得之 a 與 b 值帶入(5.9)式與(5.10)式，便可得到最佳 TSR = 0.448，最大 Cp 值為 0.115。 5.1.8. 3D 與 2D 暫態模擬之比較 為評估模擬的維度對結果的影響，另進行 2D 計算，其網格(8040)與 3D 計算中單一平面所使用的相同。 圖 5.14 為入口風速 U0 = 8m/s，轉速 N = 120rpm，2D 暫態模擬在 23 到 23.5 秒間的流線圖，與前面 3D 暫態模擬(圖 5.8)的結果比較，可發現有很大

37 的差異，圖 5.8 中流體自葉輪流出後與流經葉輪下方外圍的氣流匯合，呈一 近似穩態的流場。而圖 5.14 中可看出氣流自葉輪後方流出後即刻形成渦流， 此渦流隨著氣流往下游逸散，整體而言，呈現一週期性的非穩態結構。此 二者間的差異性可能歸因於在實際三維流場中，風機上下外圍的氣流較為 穩定，因而壓抑近風機之不穩定渦流的形成，而在 2D 流場中無此機制，因 此呈現不穩的週期性型態。 圖 5.15 則為兩者在入口風速 U0 = 8m/s，轉速 N = 120rpm 條件下單一 週期之 Cm 與時間關係圖，由圖中可以發現兩者皆與 sin 函數有類似的趨勢， 但極值的大小與出現時間有些許的不同，而 2D 模擬所觀察到的較不穩定流 場可能有利於較高極值，此外，2D 模擬結果在 Cm 上升至快接近零時出現 一小段持平的現象，此現象在 3D 模擬的結果中並沒有出現，因此推論這些 不同點可能是 2D 模擬無法完全反應實際是三維的流場性質所致。 圖 5.16 為在固定入口風速 U0= 8m/s 下，2D 和 3D 模擬的平均扭矩與風 機轉速 N 的關係圖，而該圖可發現不論是 3D 或 2D 模擬皆有相同的關係， 也就是扭矩與轉速呈線性遞減關係，除此之外，不論轉速為何，2D 模擬的 平均扭矩皆較 3D 的高。 將圖 5.16 的結果無因次化後可得到扭矩係數 Cm 與端速比 TSR 之關係， 如圖 5.17 所示，將 Cm 轉換為 Cp 則可得圖 5.18。由圖 5.17 可發現，不論是 2D 或 3D 模擬，Cm 與 TSR 呈現一線性遞減關係，且不論 TSR 值為多少，

38 2D 計算所得的 Cm 值均較 3D 的高，此現象與圖 5.15 單一週期的結果相似， 可能同樣與 2D 模擬無法完全反應三維流場性質有關，因此 2D 模擬所計算出 的 Cm 值較相同 TSR 值的 3D 計算結果為高。

5.2.

實際實驗

實際實驗

實際實驗

實際實驗

5.2.1. 風速量測 使用風源為 0.5HP 的工業風扇，量測條件共有四組，兩組不含風管， 風扇與葉輪距離分別為 0.6 與 0.7 公尺(0.5HP 0.6m、0.5HP 0.7m)，另外兩組 則是包含風管，風管與葉輪間的距離為 0.7 與 0.8 公尺(0.5HP +風管 0.7m、 0.5HP +風管 0.8m)。 風速量測的示意圖如圖 5.19 所示，量測風速的平面為包含葉輪與風源 的最近點，與風向垂直的平面，量測範圍則是以葉輪為基準，由水平方向 距葉輪左側外緣 10 公分，垂直方向距葉輪下緣 7 公分處作為起始點，接著 在水平與垂直方向以 10 公分為單位取點，水平方向共取 13 個點，垂直方 向則取 5 個點，共取 65 個取樣點，對個別取樣點分別量取風速，便可得到 入口風速分佈，圖 5.20 即為 0.5HP 風扇在距離風機 0.7 公尺(0.5HP 0.7m)條 件下的風速分佈圖，而圖 5.21 則為 0.5HP 風扇和風管在距離風機 0.7 公尺 (0.5HP+風管 0.7m)情況下的等風速線圖。比較兩圖可發現，在有風管的條 件下，於 0.4 m < x < 0.8 m 的位置，風速約在 5 ~ 7 m/s，而在沒有風管的情