一新型垂直式風機性能定性分析

國

立 交 通 大 學

機

械 工 程 學 系

碩 士 論 文

一新型垂直式風機性能定性分析

Qualitative Analysis of the Performance of a

Vertical-Axis Wind Turbine

研 究 生：林 子 翔

指

導 教 授：崔 燕 勇 教授

一新型垂直式風機性能定性分析

Qualitative Analysis of the Performance of a

Vertical-Axis Wind Turbine

研 究 生：林子翔 Student：Tzu-Hsiang Lin 指導教授：崔燕勇 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 of Science

in

Mechanical Engineering July 2011

Hsinchu, Taiwan, Republic of China

一新型垂直式風機性能定性分析

研究生：林子翔 指導教授：崔燕勇 博士 國立交通大學機械工程學系﹙研究所﹚碩士班 摘 要 本研究係發展一新型垂直式風機，此風機具有升力型及阻力型風機的 特性，並可將數個風機葉輪堆疊以得到較大的輸出功率。為瞭解其氣動特 性及風場結構，主要利用計算流力的方法做數值模擬，此對應進行並簡易 的實驗方法量測，以驗證數值分析的可靠度。數值方法為有限體積法，為 處理風場葉片的旋轉，採用多重參考座標系統，並假設在”擬暫態”下，將 風機固定在數個預先設定的旋轉角度位置，以非穩態的方式計算，如此可 簡化全暫態計算的複雜性。 結果顯示此垂直式風機的扭矩與風機轉速成線性遞減的關係，而此扭 矩與風速成二次式關係，並與風機直徑成三次關係式，將這些參數經過無 因次化後，可以得到扭矩係數與風機葉尖速度比間呈線性遞減的關係，進 而可推導出功率係數與葉尖速度比成二次式關係，並可找出最大功率係數 及最佳的葉尖速度比值。上述結果是在固定旋轉位置的擬暫態假設下得到， 在全暫態的計算(風機隨時間連續轉動)中，上述無因次參數間的關係仍然 適用，兩相比較，顯示此簡化模擬在定性上有相當的可信度，但在定量上 兩者有相當程度的差異，而全暫態計算所得的扭矩係數及功率係數高出許 多。另外實驗所得數據介於兩者之間，由此可知擬暫態分析因流體慣性被 忽略而低估扭力輸出，由於實驗設備過於簡陋，及過程有許多不嚴謹之處， 造成實驗值的不準確。

Qualitative Analysis of the Performance of a

Vertical-Axis Wind Turbine

Student：Tzu-Hsiang Lin Advisor：Dr. Yeng-Yung Tsui

Department of Mechanical Engineering National Chiao Tung University

ABSTRACT

A new type of vertical-axis wind turbine is developed in this project. This new design combines the merits of the lift type and drag type of wind turbines. It also possesses a feature by which several turbine wheels can be stacked together to get higher output power. The method of computational fluid dynamics is mainly used to investigate the flow structure and aerodynamic characteristics, also conducted is the experimental work to validate the numerical simulation. The numerical scheme is based on a finite volume method. The multiple reference frames (MRF) is adopted to tackle the rotation of the turbine wheel. To simplify the complicated problems, the quasi-unsteady state is assumed so that the wind turbine is fixed at a number of specified angular positions and unsteady computations are undertaken.

It is obtained from the simulation that the resulted torque is related to the angular speed of the wind turbine in a linear decreasing fashion, a quadratic function of the wind velocity and a cubic function of the diameter of the wind turbine. After non-dimensionalization, the resulting moment coefficient becomes a linearly decreasing function of the tip speed ratio. It can be derived from this relationship to show that the power coefficient is a quadratic function of the tip speed ratio and, furthermore, to find the maximum power coefficient and the optimum tip speed ratio. It is also shown that by using the fully unsteady model, in which the wind turbine continuously rotates in the simulation, the above correlations are still valid. However, comparing with the quasi-unsteady calculations, the resulting moment coefficient and power coefficient are much higher. The experimental results lies in between of the fully unsteady and quasi-unsteady predictions. For the quasi-unsteady on calculations, the effect of inertia is not property accounted for. As for experiments, the experimental rigs are too rough and there are a lot of uncertainties about the measurements, cause the experimental inaccuracy.

誌謝

感謝我的指導教授崔燕勇老師，您教導我的不只是課業上的知識，還有 許多做事方法與態度，並且增加了我的抗壓性，雖然常常無法達到您的要 求，但您還是很有耐心的指導我，謝謝！在此也要感謝金大仁老師實驗室 的蘇宏明、邱奕豪、黃振瑋與其他同學，在實驗上幫助了我許多。 另外感謝胡育昌與林仕文學長，總是細心地照顧實驗室的學弟妹們。感 謝陳信宏學長的照顧，雖然只有短短的一年。感謝有情有義的郭大慶，你 一定是捨不得我們才留下來的吧！感謝同學黃義政和賴胤男，碩一時很怕 meeting 一起被飆的畫面還記憶猶新！感謝學弟們黃裕堂、吳奉起和鄭凱鴻， 大家一起加油，要不然老師會在你後面非常火！最後感謝學妹陳虹汝幫助 我的論文與計畫，妳很優秀，加油！ 老人球友與重訓夥伴們，我要感謝你們幫助我走出戶外當個陽光的宅宅！ 感謝綠豆、蘇學長、邱哥、古必、Kevin 哥、型男阿康等，你們讓我的視野 更廣闊，對未來有不一樣的期許。感謝聖焱和熟哥，碩班生活因為你們而 更加充實，在實驗室壓力大的時候聽憤怒青年熟哥開砲就對了！ 媽，我畢業了！感謝我的爸媽在我的成長過程中一路支持我，相信我的 選擇與包容我的任性，讓我照著自己的步調完成我想要做的事情。最後感 謝每一個幫助過我的人，陪伴我度過這難忘的兩年，謝謝！

目錄 摘要...i ABSTRACT...ii 誌謝...iii 目錄...iv 圖目錄...vii 符號表...x 第一章 緒論 ... 1 1.1 前言 ... 1 1.2 文獻回顧 ... 4 1.3 研究內容 ... 9 第二章 數學模式 ... 10 2.1 基本假設 ... 10 2.2 多重參考座標系之統御方程式 ... 10 2.2.1 旋轉區域之統御方程式 ... 11 2.2.2 靜止區域之統御方程式 ... 11 2.2.2 均化方程式 ... 11 2.2.3 RNG k-ε Model ... 12 2.3 邊界條件 ... 16 2.3.1 進口邊界條件 ... 16

2.3.2 出口邊界條件 ... 16 2.3.3 對稱邊界條件 ... 17 2.3.4 壁面邊界條件 ... 17 第三章 數值方法 ... 18 3.1 簡介 ... 18 3.2 傳輸方程式離散化 ... 19 3.2.1 非穩態項(Unsteady Term) ... 19 3.2.2 對流項(Convection Term) ... 19 3.2.3 擴散項(Diffusion Term) ... 20 3.2.4 網格梯度的計算 ... 21 3.2.5 源項(Source Term) ... 21 3.2.6 線性代數方程式 ... 21 3.3 壓力與速度之偶合關係式 ... 22 3.3.1 計算面上的質量流率 ... 22 3.3.2 壓力修正方程式 ... 23 3.4 解題步驟 ... 24 第四章 實驗設備與方法 ... 25 4.1 實驗設備 ... 25 4.1.1 工業用風扇 ... 25

4.1.2 風速計 ... 25 4.1.3 垂直軸風力機葉片 ... 25 4.1.4 扭力計 ... 26 4.1.5 訊號擷取卡 ... 26 4.1.6 轉速計 ... 26 4.2 實驗方法與步驟 ... 27 第五章 結果與討論 ... 28 5.1 流場參數及計算網格 ... 28 5.2 流場結構 ... 28 5.3 網格數的影響 ... 29 5.4 葉片數目的探討 ... 30 5.5 轉速 N 對扭矩的影響 ... 30 5.6 入口風速 U0對扭矩的影響 ... 31 5.7 風輪直徑 D 對扭矩的影響 ... 31 5.8 平均扭矩Cm與葉尖速度比TSR 間的關係 ... 31 5.9 功率係數 Cp與葉片速度比TSR 的關係 ... 32 5.10 擬暫態與全暫態的比較 ... 33 5.11 模擬與實驗數據的比較 ... 35 第六章 結論 ... 38

圖目錄 圖1.1 水平式風力發電機【2】 ... 42 圖1.2 水平軸風車【3】 ... 42 圖1.3 垂直式阻力型(Savonius)風力機【4】 ... 43 圖1.4 垂直式 Darrieus 風力機【4】 ... 43 圖1.5 垂直式 H 型風力機【4】 ... 44 圖1.6 Darrieus 和 Savonius 混和型風機【5】 ... 44 圖1.7 改良式垂直風力發電機 ... 45 圖1.8 堆疊式垂直風力發電機 ... 45 圖1.9 各種風機之 Cp 與 TSR 關係圖【7】 ... 46 圖1.10 包覆 guide-box 之垂直軸風機示意圖【19】 ... 46 圖1.11 帶有擾流板之垂直型風機示意圖【20】 ... 47 圖1.12 螺旋式 Savonius 風車【22】 ... 47 圖1.13 垂直軸風力發電機示意圖【23】 ... 48 圖1.14 改良式垂直風力機示意圖【24】 ... 48 圖2.1 風機旋轉位置角度α ... 49 圖2.2 流場示意圖 ... 49 圖3.1 主網格與鄰近網格示意圖 ... 50 圖3.2 面上相鄰網格示意圖 ... 50

圖4.1 工業用風扇 ... 51 圖4.2 Lutron AM-4200 風速計 ... 51 圖4.3 單層風力機葉片 ... 52 圖4.4 HBM T22 扭力計(50 N-m) ... 52 圖4.5 兩種不同尺寸之 coupling ... 53 圖4.6 rinstrum R320 訊號擷取器 ... 53 圖4.7 OMEGA HHT13 數位式轉速計 ... 54 圖4.8 實驗設備示意圖 ... 55 圖5.1 流場計算示意圖 ... 56 圖5.2 固定位置角度α ... 56 圖5.3 計算網格 ... 57 圖5.4 Cm 在各角度位置 α 隨時間的變化 ... 58 圖5.5 α ₌90o_{時各時間點的流場 ... 59} 圖5.6 不同網格數Cm 值的分布 ... 60 圖5.7 不同網格數的Cm 值 ... 60 圖5.8 不同葉片數Cm 值的分布 ... 61 圖5.9 不同轉速時的扭矩 ... 61 圖5.10 不同入口風速時的扭矩 ... 62 圖5.11 不同風輪直徑時的扭矩 ... 62

圖5.12 Cm 值與 TSR 的關係 ... 63 圖5.13 Cp值與TSR 的關係 ... 63 圖5.14 全暫態計算初期 Cm 值隨時間的變化(陳虹汝提供) ... 64 圖5.15 全暫態計算後期 Cm 值隨時間的變化(陳虹汝提供) ... 64 圖5.16 不同分析方法所得 Cm 值隨風機旋轉角度的變化(N=60rpm) ... 65 圖5.17 不同分析方法所得 Cm 值隨風機旋轉角度的變化(N=90rpm) ... 65 圖5.18 不同分析方法所得 Cm 值隨風機旋轉角度的變化(N=120rpm) .... 66 圖5.19 比較全暫態及擬暫態計算模式於不同轉速時所得扭矩 ... 66 圖5.20 比較不同分析方法的 Cp值 ... 67 圖5.21 在不同位置量測風速之示意圖 ... 67 圖5.22 風扇徑向風速分布圖 ... 68 圖5.23 模擬與實驗數據的Cm 值 ... 68 圖5.24 模擬與實驗數據的 Cp值 ... 69

符號表

符號 定義 a 無因次參數間的關係係數 a 鬆弛因子 α 旋轉位置角度 A 風機受風面積 b 無因次參數間的關係係數 Cm 扭矩係數 Cm 平均扭矩係數 Cp 風車轉子效率係數 D 風機直徑

ε

紊流耗散率 C f F 對流通量 D f F 擴散通量 H 風機高度 I 紊流強度；單位張量 f J 面通量 k 紊流動能 f m 面上質量流率

N 風機轉速

ρ

密度 P 風機所產生的功率 P 壓力 r K 位置向量 R 風機半徑 S_φ 源項

τ

風機所產生的扭矩

τ

應力張量 TSR 端速比

μ

層流黏滯係數 t μ 紊流黏滯係數 0 U 入口風速 v 速度 r v 相對速度

ω

角速度 W 權重函數 上標 T 轉置矩陣

n 此一時間的已知數 n+1 下一時間的未知數 ′ 修正量 下標

φ

變數 f 面上的值 c0, p 主格點 c1, nb 相鄰格點 i 表x, y 方向

第一章 緒論

1.1 前言

隨著世界經濟發展，各國對能源的需求是逐年增加，而石化能源也逐漸 耗竭，伴隨著環保意識逐漸高漲，尋找再生性替代能源已被人類所重視。再 生性能源有許多形式，包含太陽能、風能、地熱、潮汐、生質能、水力等， 因其具有取之不盡、用之不竭的特性，且對環境的污染極低，而這些能源又 可稱為綠色能源。 風能的來源即為太陽能，太陽輻射能量約20 億分之一投射到地表上，而 其中轉換成風能的比例極低，據估算全球風能總量，相當於全世界總發電量 的 8 倍，約為地表可開發之水能資源的 10 倍，而目前已被開發的僅是微不 足道的一小部分，如果1%風能被利用，可減少全球 3%的能源消耗，可產生 世界 8~9%的總電量。【1】而風力發電擁有較高效率、低成本、無地形限制 等優點，故目前在綠色能源中為最可行的一項選擇。 傳統風力發電機有兩種分類方式。第一種依照其旋轉軸方向決定，可分 為水平式(horizontal axis wind turbine, HAWT)與垂直式(vertical axis wind turbine, VAWT)。第二種則是依照空氣作用力來決定，可分為阻力型和升力型 兩大類。由以上兩種分類方式，風力發電系統的發展大致如下：

1. 水平軸升力型風機：

其葉片數通常為 2 至 3 片，利用升力推動葉片對轉子造成扭矩，且其靜 止時葉片受到的阻力即可轉動，不需外力協助啟動，故一般商用型風力 發電機還是以此型為主。 2. 水平軸阻力型風機： 此類型如圖1.2【3】，為一典型的古荷蘭風車，通常是以農業用途居 多，因其葉片數較多，主要是靠流體的阻力來轉動，其效率較低，速度 也較慢，適用於高扭矩之工作。 3. 垂直軸阻力型風機： 一般稱為Savonius 風機(如圖 1.3)【4】，其葉片設計成 S 形，當氣流 沖擊葉片時，由於兩邊形狀的差異產生不同的壓力，而產生旋轉力矩。 雖其效率較差，但優點是不需要固定風向的限制，即任何風向皆可使其 運轉。因此風機表面積較大，也有應用在其表面上裝置太陽能電池增加 發電量。 4. 垂直軸升力型風機： 圖1.4【4】為此型風機的代表，即 Darrieus 風車，而圖 1.5【4】是 H 型風機。由於此類型風力機其完全靠升力轉動，故需要一外力協助啟動。 而此風機在高轉速下，因葉片角度固定，而會產生動力失速(stall)，故其 運轉效率也較低。若將Darrieus 風機和 Savonius 風機結合可解決其無法 自行啟動的問題(如圖 1.6)【5】。

本文所探討的為一組改良式的垂直式風力發電系統(如圖 1.7)，其結合了 阻力型與升力型風機之特點，此風機可用數個風輪堆疊分開(如圖 1.8)。由於 風機內有數個彎曲S 型葉片形成數個流道，則其作用類似於阻力型風機，而 其因流道的緣故，在不同的風吹角度下，會有升力效應的產生，如此可提高 功率輸出，改善其效能。 探討風機效能之依據： 1919 年，Betz【2】提出風車轉子效率係數 Cp，其定義如下： 3 0 1 2 P P C U A

ρ

= 其中： 風機所產生的功率 P=

τω

(τ 為產生之扭矩，

ω

代表角速度) 空氣密度

ρ

風速 U₀ 風機受風面積 A 當風機Cp值愈高，表示風機截取愈多的風能。根據貝茲極限(Betz limit)， Cp 最大值為 59%，假設流體還沒接觸到風力發電機葉片以前的風速為U0， 而流出風力發電機葉片後的速度為U₁，則在 1 0 1 3 U U = 時，會有最大 Cp值59%。

端速比(Tip Speed Ratio, TSR)λ：假設流體還沒接觸到風力發電機葉片以

前的風速為U₀，葉片尖端的速度為Utip，則端速比定義為

0

U

風車轉子效率係數Cp與端速比λ 為討論一風機效能最重要之依據。圖 1.9

【7】為各種型式風車之端速比 λ 與風車轉子效率係數 Cp關係圖。由此可知

針對一理想風車，若要有高Cp值，則必須要有高端速比λ。升力型風機(Modern

three blade, Darrieus rotor 等)則具有高端速比 λ、高 Cp值特性，其工作範圍λ

約落在 3 至 7 之間，最佳 Cp值在 35%至 45%不等，適合做為一般商業發電

用途。阻力型風機(American multi-blade, Savonius rotor 等)其轉速低、Cp值小，

一般端速比 λ 不超過 2，最佳 Cp值也在 30%以下，其中 Savonius 型風機 Cp 值則約20%左右，故阻力型風機較適合高扭矩之工作。

1.2 文獻回顧

早期因為電腦運算速度較慢及記憶單位容量較小，所以在模擬風力發電 機的研究裡，發展出不同的簡化方法，而隨著電腦運算速度快速成長，現今 對此類流場的模擬也愈來愈完整。

對水平式風力機的設計，最廣泛的方法是利用Blade Element Momentum

(BEM) method【8】，此方法假設二維流場，從已知的機翼升力係數(CL)及阻

力係數(CD)推得外力，再加入動量方程式中求解，但此方法必須要有準確的

機翼動力資訊，才會得到較高的準確性與性能預測。而為了探討葉片末端所

產生的尾流，發展出actuator disc【9,10】及 actuator line【11,12】等方法，

前者針對葉片後流場建立圓盤型計算區域，而經由轉動葉片影響後的流體通 過此區域，進而探討其流場物理性質；後者則是使用螺線作為計算區域，以

節省計算量。而要探討葉片幾何形狀的影響，則需利用與物體邊界吻合網格

的全三維計算方法【13,14】。而在垂直式風力機方面，也有一些簡化的方法，

不過要探討其流場物理性質和葉片形狀的影響，還是需要利用與物體邊界吻 合網格的全三維計算方法。

水平軸式風力發電機的研究：

2002 年，Sørensen 和 Shen【11】使用 actuator line 法簡化直徑 40 米的 Nordtank NKT500 型水平式風機，此風機有三個 LM19.1 葉片，正常功率輸

出為500kW。結果發現風速愈高(Wind Speed=12m/s)，模擬與實驗上的誤差

會來到5%，推測原因可能是機翼動力資訊不夠精準所致，而風速超過 12m/s

後，流場會產生動力失速現象，使輸出功率成長趨於平緩甚至下降。

2006 年，Mandas et al.【14】探討一水平式風機，比較傳統 BEM 方法與

全三維計算方法的結果。全三維計算方法採用 MRF 求解，並使用不同的紊

流模式(Spalart-Allmaras Model 與 SST Model)。而結果發現在低風速時，BEM 與全三維計算的動力輸出趨勢是一致的，但風速提高後卻有很大的差異，推 測為高轉速時機翼產生的動力失速影響所致，故在高轉速時適用全三維的計 算方法。

2006 年，Hu et al.【15】利用 BEM method 探討一水平式風力機葉片(NREL S809)的動力失速(stall-delay)現象，利用 RNG k-ε Model 做紊流分析，且使用

葉片來的好。而葉片在高轉速時，最大功輸出並不高，推測是動力失速現象 導致。

垂直軸式升力型風力發電機的研究：

2007年，Simão Ferreira et al.【16】針對二維的Darrieus風機做一完整的 探討。此風機所採用的葉片為NACA 0012型號，由於Darrieus風機屬於一升 力型風機，故在葉片上也會有動力失速現象產生。而此研究也做了四種不 同紊流模式的比較，分別為Spalart-Allmaras Model, Standard k-ε Model和 large eddy models (LES and DES Model)，結果發現DES Model的結果最接近

實驗數據。URANS models由於對大渦流的計算不夠準確導致整體準確性較

低。LES的結果略低於DES，是在靠近牆的區域不夠準確造成的。

2008年，Antheaume et al.【17】利用Vortex models方法簡化一垂直式 Darrieus風車，使用NACA0012型號之葉片，紊流模式為k-ε Model並和實驗 作對照。而此研究除了討論不同工作流體的效率，也把多個風力發電機串 接在一起討論整體平均效率。結果發現要使平均效率提升，除了縮短風力 發電機之間的距離，也可以增加多組風力發電機。另外工作流體為空氣的 效果遠小於水。而在堆疊多個風機時，二維的整體效率的增加效果比三維 的情況來的更好。 2009 年，Howell et al.【18】討論一垂直型升力風機，使用 CFD 模擬 2D 與3D 流場性質，並和實驗數據做比較。紊流模式採用 RNG k-ε Model，分

析不同葉片數與 2D、3D 流場的差異。模擬結果顯示雙葉片與三葉片動力

輸出趨勢相似，但三葉片僅需較低的TSR。而在 2D 的動力輸出也較 3D 與

實驗數據高，造成此現象可能是在葉片尖端的數個大型渦流影響其升力結 果導致。

垂直軸式阻力型風力發電機的研究：

2007 年，Irabu 和 Roy【19】運用 guide-box 包覆一 Savonius 風機(圖 1.10) 【19】，藉此降低逆向扭矩，以達到提升功率輸出目的。結果顯示在兩葉片

風機可提升1.23 倍的功率輸出，而對於三葉片風機則提升約 1.5 倍。此外，

兩葉片型風機抽電效率略高於三葉片型，在啟動時更為明顯。

2008 年，Altan 和 Atılgan【20】設計一組擾流板置於 Savonius 風車前緣

(如圖 1.11)【20】，測試在固定的靜止風機角度下(45o、60o、90o)擾流板對產

生扭矩的影響。此研究使用Fluent 6.0 做為模擬工具，紊流模式為標準 k-ε

Model，並與實驗比較。結論指出擾流板愈長效果愈好，最大扭矩發生在風

機角度 60o與擾流板角度 β=15o、α=45o時，擾流板對提升扭矩有明顯的效

果，且模擬與實驗數據大致吻合。

2008 年，Saha et al.【21】利用一風洞探討兩種不同葉片形狀之 Savonius 風車轉子，包括半圓形與扭曲葉片，其葉片上設計一單向進氣閥，藉此提 升扭矩。此實驗探討不同葉片形狀、葉片數量、堆疊層數與進氣閥的影響。

而葉片數增加，功率輸出隨之下降，對於Savonius 型風車兩葉片最為理想。

在葉片型狀方面，扭曲葉片Cp值也較一般半圓形葉片高。最後為進氣閥對

功率輸出的影響，對於三葉片 Savonius 型風車，裝置一進氣閥可以有效增

加其扭矩輸出。

2009 年，Zhao et al.【22】探討一螺旋式 Savonius 風車轉子(如圖 1.12)【22】，

且對其幾何外型做不同設計之測試。而其運用線性 k-ε Model 及全三維計算 方法得到完整的流場物理性質，結果發現對此螺旋式Savonius 風車轉子而言， 二葉片的風車轉子效率係數Cp較三葉片好，而在任何角度下，風對靜止轉子 的轉矩都是同方向的。 2009 年，Hu 和 Tong【23】探討一改良式的阻力型垂直風機，由於此風 機傳統型效率並不高，故在迎風面處加一圓弧形檔板(如圖 1.13)【23】，以 減少反方向之轉矩。紊流模式為RNG k-ε Model 並用 SIMPLE 演算法求解 二維流場。結果顯示檔板角度約30o時提升轉矩會有最好的效果。

2009 年，Kamoji et al.【24】使用風洞實驗測試一改良型 Savonius 風車(圖 1.14)【24】，探討有無中心軸對功率輸出的影響，另外也測試了葉片重疊比 例、展弦比、葉片圓弧角度與雷諾數的影響。最後結果顯示此阻力型垂直 風機無中心軸所產生的功率輸出較有中心軸高。

2010 年，Mohamed et al.【25】探討二葉片與三葉片 Savonius 風機，並

模式為 realizable k-ε Model，另外針對風機的旋轉使用滑移網格(Sliding Mesh)。結果顯示無論擋板角度為何，皆有助於扭矩的提升。對於二葉片風 機Cp值可提升27.3%，三葉片風機則是 27.5%。而無論在 Cp值、製造費用、 製造難度與複雜性皆為二葉片風機較理想。

1.3 研究內容

一般風機氣動特性的研究需使用大型風洞進行實驗【26,27】，除了經由 實驗量測與觀察，個人電腦的運算速度提高，可以更加的快速準確地進行模 擬。在此主要採用計算流體力學(Computational Fluid Dynamics, CFD)的方法 來探討，並以簡易的實驗做為輔助。

本文使用商用軟體ANSYS Fluent 做為模擬工具，對於風機葉片轉動的處

理，將採用多重參考座標系統(Multiple Reference Frames, MRF)，在風機葉片

旋轉範圍內為旋轉區，外圍則為靜止區。而紊流流場以RNG k-ε 模式配合壁 函數去模擬風力發電系統之流場。本文將探討風機在不同角度、風速、轉速、 直徑下時所產生的扭矩，進而探討其間的關係。 由於缺乏大型風洞實驗，且尚在初期開發階段，有諸多問題急需改善， 故本實驗僅測試單層風輪，實驗結果將用於對數值計算的結果，作一粗略的 評估。

第二章 數學模式

2.1 基本假設

所分析的垂直式風力發電機之工作流體為空氣，對此流場我們作以下幾 項假設： １. 流場擬似暫態 原先我們假設流場為一擬穩態(quasi-steady)系統，定義一風機旋轉 位置角度α，為 S 型葉片尖端至風機中心與水平軸的夾角，固定不同的 α (圖 2.1)，即利用 snap shot 的方式將流場視為穩態。但模擬後發現， 在流場變化較劇烈的某些旋轉位置角度α 下無法得到穩態解，故後續模 擬採用擬似暫態方式求解。 ２. 流場為二維不可壓縮流且忽略重力的影響。 由於本研究實驗部分測試單層風機效能(如圖 1.7)，故模擬將其簡化 為二維流場。 ３. 流場為等溫狀態。

2.2 多重參考座標系之統御方程式

分析垂直式風力發電機必須將計算區域分為兩個部分(如圖 2.2)，一為風 輪掃過的旋轉區域，外圍部分則為靜止區域。旋轉區域內給定一角速度ωG， 靜止區則ωG =0，以下就兩區域統御方程式分別介紹：

2.2.1 旋轉區域之統御方程式 連續方程式 ∇ ⋅

ρ

vr =0 G (2.1) 動量方程式 v

( ) ( )

v vr v p F t

ρ

ρ

ρ ω

τ

∂ + ∇ ⋅ + × = −∇ + ∇ ⋅ + ∂ G G G JG G JG (2.2) 其中： 相對速度 vr = − ×v ω r G G JG G (2.3) 在此 rG為位置向量 應力張量

(

)

2 3 T v v vI

τ μ

= ⎡_⎢ ∇ + ∇ − ∇ ⋅ ⎤_⎥ ⎣ ⎦ G G G (2.4) 在此I 為單位張量 在(2.2)式中

ρ ω

( )

JG G×v 項包含因旋轉座標產生的體力(body force)。 2.2.2 靜止區域之統御方程式 靜止區域中，因角速度ωG =0，故(2.3)式為vr =v G G ，動量方程式中體力項

( )

v

ρ ω

JG G× 則為零。因此靜止區域之統御方程式如下： 連續方程式 ∇ ⋅

ρ

vG=0 (2.5) 動量方程式 v

( )

vv p F t

ρ

ρ

τ

∂ + ∇ ⋅ = −∇ + ∇ ⋅ + ∂ G GG JG (2.6) 以下介紹之統御方程式皆以靜止區域統御方程式表示： 2.2.2 均化方程式 因流場為一非穩態紊流流場，故透過ensemble averaging 過程，形成一 均化方程式。其定義如下：

對於一瞬時速度u ，可分為平均項_i u 與變動項i ' i u 。 ' i i i u = +u u (2.7) 將(2.7)代入(2.5)(2.6)中並進行時均化過程後可得到以下形式： 連續方程式

( )

_i 0 i u x

ρ

∂ = ∂ (2.8) 動量方程式

( )

_i

(

_{i j}

)

i u u u t

ρ

x

ρ

∂ ∂ + = ∂ ∂ (2.9)

(

' '

)

j i i j i i j i j u p u u u x x

μ

x x x

ρ

⎡ ⎛ ∂ ⎞⎤ ∂ ∂ ∂ ∂ − + ⎢ ⎜_⎜ + ⎟_⎟⎥+ − ∂ ∂ _{⎢⎣ ⎝}∂ ∂ _⎠⎥⎦ ∂ 上式中

(

−

ρ

u u_i' '_j

)

為雷諾應力，此項需要利用紊流模式來求解。 2.2.3 RNG k-ε Model 根據Boussinesq(1877)的假設，雷諾應力與平均應變率呈線性關係： ' ' ₂ 2 3 i j t ij ij u u D k

ρ

μ

ρ δ

− = − (2.10) 其中： 平均應變率 1 2 i j ij j i u u D x x ⎛_{∂ ∂} ⎞ = ⎜_⎜ + ⎟_⎟ ∂ ∂ ⎝ ⎠ (2.11) 紊流黏滯係數 μ_t 紊流動能 k 紊流耗散率

ε

常數 C_μ =0.0845

(2.10)式中等號右邊加入 2 3

ρ δ

k ij − 項是為了當i=j 時，能夠使等式兩邊成立。 將(2.10)代入(2.9)中，動量方程式可重新整理為

( )

(

)

2 3 j i i i j eff ij i i i j i u p u u u u k t

ρ

x

ρ

x x

μ

x x

ρ δ

⎡ ⎛ ∂ ⎞ ⎤ ∂ ₊ ∂ _{= −}∂ ₊ ∂ ∂ _{+ −} ⎢ ⎜_⎜ ⎟_⎟ ⎥ ∂ ∂ ∂ ∂ _{⎢⎣ ⎝}∂ ∂ _{⎠ ⎥⎦} (2.12) 其中：

μ

_eff = +

μ μ

_t (2.13)

μ

為層流黏滯係數 t μ 為紊流黏滯係數

RNG k-ε Model 是利用 Re-Normalisation Group (RNG) 方法由瞬時 Navier-Stokes 方程式推導得到，與標準線性 k-ε Model 比較有以下特點：

• 在ε 方程式中增加了一 R_ε 項以增加應變流的準確性。

• 對於旋轉流體有更佳的準確性，較適合解有渦流產生的流場。

• RNG 理論使用公式求解紊流普朗特常數(turbulent Prandtl numbers)，

而標準k-ε Model 則是固定常數。

• 標準k-ε Model 僅適用於高雷諾數流場，而 RNG k-ε Model 可適用於

低雷諾數流場。

準確性與可信度。k 與 ε 的方程式如下： k 方程式

( )

(

_i

)

_{k eff k} i i j k k ku G t

ρ

x

ρ

x

α μ

x

ρε

⎛ ⎞ ∂ ₊ ∂ ₌ ∂ ∂ _{+ −} ⎜ ⎟ ⎜ ⎟ ∂ ∂ ∂ _⎝ ∂ _⎠ (2.14) ε 方程式

( )

(

_i

)

i u t

ρε

x

ρε

∂ ∂ + = ∂ ∂ (2.15) 2 1 2 eff k i j C G C R x ε x ε k ε k ε

ε

ε

ε

α μ

ρ

⎛ ⎞ ∂ ∂ _{+ − −} ⎜ ⎟ ⎜ ⎟ ∂ _⎝ ∂ _⎠ k G 為紊流生成項

( )

' ' i 2 i k i j t ij j j u u G u u S x x

ρ

∂

μ

∂ = − = ∂ ∂ (2.16) 1 2 j i ij j i u u S x x ⎡_∂ ∂ ⎤ = _⎢ + _⎥ ∂ ∂ ⎢ ⎥ ⎣ ⎦ (2.17)

此處S 代表平均應變率張量(mean rate-of-strain tensor) _ij

紊流黏滯係數μ_t RNG 理論計算紊流黏滯係數所使用公式為： 2 3 ˆ ˆ 1.72 ˆ 1 _v k v d dv v C

ρ

εμ

⎛ ⎞ = ⎜ ⎟ ⎜ ⎟ _{− +} ⎝ ⎠ (2.18) 其中： ˆ eff v

μ

μ

= 、C_v ≈100

(2.18)式可適用於低雷諾數與靠近壁面之流場，在高雷諾數流場此式則為 2 t k C_μ

μ

ρ

ε

= (2.19)

Inverse Effective Prandtl Numbers α_k、α_ε

在RNG 理論中計算α_k與α_ε公式如下： 0.6321 0.3679 0 0 1.3929 2.3929 1.3929 2.3929 mol eff

α

α

μ

α

α

μ

− + = − + (2.20) 其中α₀ =1.0 若在高雷諾數流場中( mol 1 eff

μ

μ

 )，則αk =αε ≈1.393。 而RNG k-ε Model 與 k-ε Model 最主要的差異為 ε 方程式中的 R_ε 項，其 表示為： 3 2 0 3 1 1 C R k μ ε

η

ρη

_{η ε}

βη

⎛ ₋ ⎞ ⎜ ⎟ ⎝ ⎠ = + (2.21) 其中

η

Sk

ε

≡ 、S ≡ 2S S_{ij ij} 、η0 =4.38、

β

=0.012 將(2.21)代回(2.15)整理後可得到 ε 方程式最後形式為：

( )

(

)

* 2 1 2 i eff k i i j u C G C t x x ε x ε k ε k

ε

ε

ε

ρε

ρε

⎛

α μ

⎞

ρ

∂ ₊ ∂ ₌ ∂ ∂ _{+ −} ⎜ ⎟ ⎜ ⎟ ∂ ∂ ∂ _⎝ ∂ _⎠ (2.22) 在此 3 * 0 2 2 3 1 1 C C C μ ε ε

η

η

_η

βη

⎛ ₋ ⎞ ⎜ ⎟ ⎝ ⎠ ≡ + + (2.23)

在應變率小的流場中(η η< )，₀ * 2 C _ε會大於C_2ε ，此時C₂*_ε ≈2.0，與標準 k-ε Model 中C_2ε =1.92相近。而在η η> 時₀ C₂*_ε也小於C_2ε ，這會使RNG k-ε Model 比標準 k-ε Model 具有更低的紊流黏滯係數。另外在 RNG 理論裡， (2.15)式中常數為C_1ε =1.42、C_2ε =1.68。

2.3 邊界條件

本文欲探討之流場為一開放區域(如圖 2.2)，包含進出口、對稱邊界與固 體壁面，其中固體壁面為一轉動中的風輪所構成的移動壁面。以下將分別 探討之。 2.3.1 進口邊界條件 定義進口的速度分佈為一水平方向固定流速之均勻流，對於進口的 k 與ε 值定義如下： 2 3 ( ) 2 in k = U I ； 1 2 t k C_μ

μ

ε ρ

μ μ

− ⎛ ⎞ = _{⎜ ⎟} ⎝ ⎠ 其中： in U 進口速度，即圖2.2 中 inlet 處。 I 紊流強度(turbulence intensity)，在此假設為 10%。 t

μ

μ

紊流黏滯係數比例(Turbulent Viscosity Ratio)，在此假設為 100。

2.3.2 出口邊界條件

其適用於出口壓力未知之內部流場與外部流場。

2.3.3 對稱邊界條件

ANSYS Fluent 內提供的對稱邊界條件(Symmetry Condition)假設以下兩 種情況：

• 對稱邊界面上垂直方向速度為零

• 垂直方向流場各性質梯度為零

2.3.4 壁面邊界條件

由於流場為非滑移流，故在壁面上流體的速度會等於固體壁面的移動速

度。且RNG k-ε Model 在壁面上將使用標準壁函數(Standard Wall Function )

第三章 數值方法

3.1 簡介

ANSYS Fluent 為目前市占率最大的計算流體力學分析軟體，其分析結果 快速準確，故本研究使用 ANSYS Fluent 做為模擬工具。此軟體提供了流場 分析、熱傳、質傳、化學反應及燃燒等模擬。目前版本(13.0.0)已將建立物理 模型(DesignModeler)、網格生成與切割(Meshing)與後處理整合到 Fluent 模組 裡。使用者僅需繪製物理模型與進行網格建立分割工作，再以有限體積法 (Finite Volume Method, FVM)求解，讓使用者在短時間內可以快速了解流場各 種變化，以縮短開發產品所需時間。

雖然ANSYS Fluent 13.0 版已將許多功能整合至模組裡，但本研究網格建

立與分割工作採用另一套軟體Gambit，此軟體可支援 SolidWorks 與 Pro/E 圖

檔，並可讓使用者自訂網格參數等，其自由度相對於 ANSYS Fluent 模組較 高，且可將建立好的網格匯至ANSYS Fluent 求解。 本研究計算範圍分別為風輪掃過的區域和此區外的二維流場。由於此流 場是非穩態系統，需要大量的電腦計算時間，為了節省計算量，我們將其 簡化假設為擬暫態系統。紊流模式則為RNG k-ε 模式。 本章目的為介紹ANSYS Fluent 所使用的數值方法。3.2 節處理動量方 程式離散化，3.3 節為壓力與速度之偶合關係式，在 3.4 節中介紹 SIMPLE 演算法解題步驟。

3.2 傳輸方程式離散化

根據第二章的連續方程式、動量方程式及k 和 ε 方程式，可用下列的通 式來表示：

( )

vr S t

ρφ

ρφ

φ

φ

φ ∂ + ∇ ⋅ = ∇ ⋅ Γ ∇ + ∂ G (3.1) 接下來以有限體積積分法對(3.1)式進行離散化，以下就非穩態項、對流項、 擴散項及源項分別詳述： 3.2.1 非穩態項(Unsteady Term) 對非穩態項直接對其取體積分：

( )

V _t dV _t V

ρφ

ρφ

∂ ∂ = ∂ ∂

∫

(3.2) 上式等號右邊 t

φ

∂ ∂ 可表示為：

( )

F t

φ

_φ

∂ = ∂ (3.3)

對(3.3)式離散化 ANSYS Fluent 預設採用 First Order Implicit：

( )

1 1 n n n F t

φ

+

φ

_φ

+ − ₌ Δ (3.4)

( )

1 1 n n n tF φ + ₌φ _{+ Δ} φ + (3.5) n

φ

為這一個時間的已知數、為

φ

n+1_{下一個時間的未知數。} 3.2.2 對流項(Convection Term)

ρφ

∇ ⋅ G

將體積分轉化為面積分，再以中點定理(Midpoint rule)將面積分轉為差分形 式： r r V∇⋅

ρφ

v dV =

ρφ

v d A⋅

∫

G

v

∫

G JK (3.6) C f r _{f f} rf _{f f f} f f f v d A v A m F

ρφ

⋅ ≈

∑

ρ φ

⋅ =

∑

φ

=

∑

∫

G JK G JK 

v

(3.7) 其中： f

m 質量流率(mass flow rate)

C f

F 對流通量(flux)

在 求 解 的 過 程 中 ， 對 於 對 流 項 的 處 理 ， 是 採 用 二 階 上 風 差 分 法 (Second-Order Upwind Scheme)作運算，可保持較高準確度：

(

)

C f f f f F =m

φ

=m

φ

+ ∇ ⋅

φ

rK (3.8) 其中：

φ

上游網格中心之值

φ

∇ 上游網格中心之梯度 r K 上游網格中心至面上中點之位置向量 3.2.3 擴散項(Diffusion Term) 同樣對擴散項∇ ⋅ Γ ∇ 進行有限體積積分，藉由高斯定理轉換體積分成_φ

φ

面積分，再由中點定理將面積分轉化為差分形式，由此可得： V∇ ⋅ Γ ∇φ

φ

dV = Γ ∇ ⋅φ

φ

d A

∫

v

∫

JG (3.9)

(

)

D f _f f f f d A S F φ

φ

φ

φ

Γ ∇ ⋅ ≈

∑

Γ ∇ ⋅ =

∑

∫

JG JK

v

(3.10)

3.2.4 網格梯度的計算

為了計算對流項與擴散項中的∇

φ

，ANSYS Fluent 13.0 使用 Least

Squares Cell-Based 方法，其定義如下：

( )

∇φ _c0⋅ Δ =ri

(

φci −φc0

)

(3.11) 在此Δ 為主格點至相鄰格點中心之距離(如圖 3.1)，可將上式整理為： r_i

[ ]

J

( )

∇φ _c0 = Δφ (3.12) 其中：

[ ]

J 代表由幾何形狀導出方程式之係數矩陣

(

ci c0

)

φ

φ

φ

Δ = − 0 xiˆ yˆj

φ φ

φ

∇ = + (3.13) (3.13)式中使用權重函數 W【28】可得各分量：

( )

₀ 0

(

0

)

1 n x x _c i ci c i W

φ

φ

φ

= =

∑

⋅ − (3.14)

( )

₀ 0

(

0

)

1 n y y _c i ci c i W

φ

φ

φ

= =

∑

⋅ − (3.15) 3.2.5 源項(Source Term) 對源項 S_φ可直接對其進行體積分： V S dVφ ≈ S Vφ

∫

(3.16) 3.2.6 線性代數方程式 將非穩態項、對流項、擴散項及源項合併後，可得離散化之傳輸代數方

程式： p nb nb nb a φ =

∑

a φ +b (3.17) 其中下標p 為主格點，下標 nb 為相鄰格點 。 為了使疊代過程較穩定，引入一鬆弛因子 a (under-relaxation factor)： old a φ φ= + Δφ (3.18) 其中下標 old 代表前一次疊代之值。其對動量代數方程式之修正如下式： 1 p nb nb p old nb a _a a b a a a

φ

φ

−

φ

=

∑

+ + (3.19)

3.3 壓力與速度之偶合關係式

3.3.1 計算面上的質量流率

為滿足連續方程式，須計算面上的速度。ANSYS Fluent 引用 Rhie【29】

所提出方法求得面上質量流率：

(

)

f n f f f m =

ρ

v AJJG =J A (3.20) 在此

( )

(

)

(

( )

)

(

)

, 0 , 0 , 1 , 1 0 ₀ 0 1 ₁ 1 , 0 , 1 p c n c p c n c f f f c c c c p c p c a v a v J d p p r p p r a a ρ + = + + ∇ ⋅ − + ∇ ⋅ + JG JG (3.21)

(

0 1

)

ˆ f f c c J d p p = + − 其中下標c0,c1 為面上相鄰兩格點，rJG₀與rJG₁為網格中心至面中心位置向量(見 圖3.2)，而d 則為_f a 的函數。 _p

3.3.2 壓力修正方程式 根據 Patankar 所提的 SIMPLE 法則，利用前次疊代的壓力 p∗，可求解 動量方程式，並可得到面上通量： J*_f = Jˆ*_f +d_f

(

p_c*₀ − p_c*₁

)

(3.22) 而上式並不滿足連續方程式，因此需再加入修正項J'_f： * ' f f f J =J +J (3.23) 其中J'_f =d_f

(

p_c'₀ − p_c'₁

)

， p 為壓力修正量。 ' 令修正後滿足連續方程式： 0 face N f f f J A =

∑

(3.24) 將(3.22)與(3.23)代入(3.24)整理後得到壓力修正方程式： _{p nb nb}' nb a p′ =

∑

a p +b (3.25) 其中源項b 為網格淨質量流率： * face N f f f b=

∑

J A (3.26)

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

力修正方程式(3.25)後可得壓力修正量 'p ，並修正面通量J ： _f

p= p*+a p_p ' (3.27)

J_f = J*_f +d_f

(

p_c'₀ − p_c'₁

)

(3.28)

3.4 解題步驟

在此介紹解題步驟為一個 t+ 下的流程： I. 給予初始速度v 、壓力0 p 、0 k 與0

ε

0分佈及紊流黏性

μ

_t0。 II. 求解動量方程式。 III. 計算質量流率m ，由壓力修正方程式解得 p′ 。 _f IV. 修正速度、壓力及質量流率。 V. 求解 k 與 ε 方程式，並計算新的紊流黏性μt ∗ 。 VI. 將求得的新值作為初始值，重覆步驟 II~V，直到獲得一個 t+ 下收斂解。

第四章 實驗設備與方法

本研究除了使用模擬探討外，也與實驗結果做比較。實驗部分由交通大 學複合材料結構實驗室協助完成，在此除了扭力計、訊號擷取卡、轉速計 外之儀器皆由該實驗室所提供。本章目的為介紹實驗設備與實驗步驟，實 驗結果將在第五章與模擬結果討論。

4.1 實驗設備

4.1.1 工業用風扇 本實驗為一開放性流場，採用工業用風扇做為風源，圖4.1 為該風扇之 影像圖，為一單相6 極之四葉片型風扇，其葉片直徑約 0.59m，中心位置與 邊緣風速較弱，靠近二分之一半徑處風速最強。 4.1.2 風速計 為了測量風扇所產生之風速，採用Lutron 公司生產之 AM-4200 風速計 (圖 4.2)量測。其量測範圍為 0.8 ~ 30.0 m s ，精確度 3%± 。 4.1.3 垂直軸風力機葉片 圖4.3 為本研究所探討的垂直式風機葉片實體，由 3 個 S 型葉片與 2 個 翼型葉片組成，包含上頂板與下底板。風機直徑為 1m，高度為 0.28m。其 製作採用樹脂與玻璃纖維布之複合材料，優點抗腐蝕耐磨耗，且重量輕， 約為鐵的二分之一。

4.1.4 扭力計 本研究量測扭力值使用德國HBM 公司製造之 T22 扭力計(圖 4.4)，此型 號扭力計可量測扭力值最小值為0.5%，在此選用 50 N-m 型 T22，即量測範 圍為0.25~50 N-m。本設備需接被測試端與驅動端，在本研究則為風力機葉 片與發電機。由於扭力計、風機葉片、發電機軸徑皆不同，故須使用coupling 連接(圖 4.5)。此外量得扭力值經轉換後輸出為直流電壓 5V± ，並由訊號擷 取設備接收顯示。 4.1.5 訊號擷取卡 本設備即用來接收扭力計測得之數據，在此採用 Rinstrum 公司製造之 擷取卡，型號為R320(圖 4.6)。由扭力計輸出之直流電壓 5V± 進入擷取卡後 顯示於液晶螢幕上，其對應扭力計量測範圍，即+5V 為50 N-m， 5V− 則為 -50 N-m。 4.1.6 轉速計 本實驗使用美國OMEGA 公司製造之數位式轉速計，型號為 HHT13(圖 4.7)。量測範圍為5 ~ 200,000 rpm ，誤差 0.01%± 。其作用原理係利用紅外 線反射，將反光貼紙貼在軸心上，由轉速計照射紅外線至反光貼紙上，反 射後便可接收計算轉速。本實驗將轉速計置於風機外 3 公尺處，以避免破 壞流場結構。

4.2 實驗方法與步驟

圖 4.8 為 4.1 節介紹之實驗設備連接後示意圖。測量方法首先將扭力計 校正歸零，再開啟風扇使風機葉片運轉，靜待數分鐘後等風機旋轉穩定後， 記錄扭力數值與轉速值。此外測量不同條件皆須等待風機旋轉穩定後再採 樣數據，以降低誤差值的產生。 本實驗測量工業用風扇距風機0.7 米時的輸出扭矩，首先量測無負載時 的自由轉速，即沒有連接發電機，僅靠風速推動之風機轉速。其次再串接 不同大小電阻(2.5, 5, 10, 20, 30, 40, 50 Ω)，記錄各種負載下的風機轉速與 輸出扭力值，並將在後面章節與模擬數據做一討論。

第五章 結果與討論

5.1 流場參數及計算網格

此流場目前考慮二維部分，計算範圍如圖 5.1 所示，工作流體為空氣， 入口速度為一均勻流U₀，其中風輪直徑 1m，L 則為 50 倍之風輪直徑，旋 轉區域為 1.1 倍風機直徑(虛線範圍內)，固定位置角度

α

(圖 5.2)，每 30o做 一次計算。

α

= 時，計算網格如圖 5.3 所示，其中包含了風輪近處之網格，0o 紅線區域內則為旋轉區。 本研究測試了不同流場的參數，先探討網格數量對計算結果之影響，接 著比較不同葉片數的扭矩輸出，其中入口速度(U₀ =4, 8, 12, 16 m s/ )，風 機轉速(N =30, 60, 90, 120 RPM )，風機直徑(D=0.5, 1, 1.5, 2 m)，並縮小

α

間隔至10o，最後再將擬暫態結果與全暫態計算結果做一比較。

5.2 流場結構

由於先前擬穩態的假設中將風輪固定在某些特定旋轉角度位置，然而除 了當風輪葉片約略平行於風向時方可得到穩態解，在其他角度流場呈非穩 態變化，故使用暫態求解。在入口風速U₀ =8 m s、轉速N =60 RPM 下， 圖5.4 顯示扭矩係數(Moment Coefficient) Cm 隨時間的變化圖，Cm 定義如 下： 2 0 1 2 Cm U AR

τ

ρ

= (5.1)

其中

τ

代表產生之扭矩(n-m)，A 代表風機受風面積(m2)，R(=D/2)代表風機 之半徑(m)。 從圖中可以看出只有在角度位置

α

= 時 Cm 到達穩定值，當0o

α

=₆₀o_、 90o_、120o_{時，Cm 成近乎週期性的變化，在}

α

=₁₅₀o_{時其改變則無固定振幅。} 在後面的結果，係將時間 10sec 以後到約 20sec 間的 Cm 值取其平均(Cm) 做為代表。 上述的週期及近乎週期的非穩態現象是由於當風機角度偏離 0o 時，其

葉片相對風向成鈍體形態(bluff body)，因而產生類似 Karman vortex street

的不穩定流場，此現象可由圖 5.5 得到證明，此圖顯示在

α

=90o_{(垂直於風} 向)時，在一週期變化中四個時間點於風機附近流場的變化情形，在頭一個 時間點，緊鄰風機後方的渦流在風輪上方形成，並以順時鐘方向旋轉，在 第二個時間點，此順時針渦流往下游傳遞，並逐漸脫離風輪，此時 Cm 值 增加至最高點，之後在第三個時間點，新的渦流於下方形成，其轉動方向 與上述渦流相反，呈逆時針方向旋轉，在週期末了，此渦流也逐漸脫離風 輪，此時Cm 值降到最低。

5.3 網格數的影響

圖 5.6 為入口風速U₀ =8 m s、轉速N =60 RPM 時，使用不同網格數 於各角度Cm 值的分布，由圖可知網格數在 8003 至 42180 時各角度 Cm 值 分布趨勢相同，但網格數在8003 至 23108 時平均扭矩係數(Cm)值有很大的

變化(見圖 5.7)，而在網格數 23108 與 42180 時的Cm值相差約1%，故之後 計算選用網格數23108，以減少實際計算所需時間。

5.4 葉片數目的探討

圖5.8 為單 S 型葉片(1s)、三 S 型葉片(3s)、單 S 型加兩翼型葉片(1s+2a) 與五葉片風機(3s+2a)各角度 Cm 值的分布(入口風速U₀ =8 m s、轉速 90 N = RPM )，在

α

=30o時三葉片、五葉片及單 S 加翼形葉片風機皆有較 大 Cm 值，而在

α

=90o_至120o_{產生較低 Cm 值。單 S 型葉片風機 Cm 值分} 布趨勢與前者不同，Cm 值變化也較小，其平均扭矩輸出五葉片風機最高 (Cm=0.0051)，其次為單 S 型加兩翼型葉片風機(Cm= −0.0076)與單 S 型葉 片風機(Cm= −0.0325)，而三 S 型葉片風機輸出為最低(Cm= −0.0515)，此 結果顯示本研究所採用的五葉片風機(3s+2a)有較佳的性能。

5.5 轉速 N 對扭矩的影響

圖5.9 是U₀ =8 m s及 D=1m 時，轉速N =30, 60, 90, 120 RPM 所產生 的扭矩，很明顯的，其間為一線性遞減的關係，值得注意的是在N =90 RPM 時扭矩值約為零，此表示在風速U₀ =8 /m s時，此風輪的自由旋轉速度約為 90RPM。

5.6 入口風速 U

0

對扭矩的影響

圖 5.10 顯示在N =90 RPM 及 D=1m 時，入口風速U₀ =4,8,12,16 /m s的 扭力值，可以看出扭力隨著風速遞增，其間呈現一非線性關係，在後面我 們將證明其間是一個二次的關係式。

5.7 風輪直徑 D 對扭矩的影響

圖5.11 顯示在U₀ =8 m s時，於兩種轉速下 N=45 及 90RPM，改變風輪 直徑 D 對扭力的影響，可以看出存在一最佳直徑值，此時扭力最大，在後 面我們將證明直徑與扭力間存在一個三次的關係式。

5.8 平均扭矩Cm與葉尖速度比 TSR 間的關係

葉尖速度比(TSR)為一重要無因次參數，其大小將決定風機的性能及效 率，其定義為葉片尖端速度與風速的比值： 0 R TSR U

ω

= (5.2) 在此角速度 2 60 N

π

ω

= 將上面圖 5.9 至圖 5.11 的結果無因次化後顯示於圖 5.12，很明顯的，無因 次的扭矩係數Cm及葉尖速度比 TSR 呈現一線性遞減關係： Cm= −aTSR b+ (5.3) 上式中的ａ,b 值為正值，在本計算中其值約為 a=0.4、b=0.24，將Cm及TSR 的定義代入上式，可得有因次參數間的關係：

2 2 0 0 240 4 a b A D NU ADU τ = − ρ π + ρ (5.4) 風機截面積A 可寫成A=DH，H 是風機高度，如此可得： 3 2 2 0 0 240 4 a b H D NU HD U τ = − ρ π + ρ (5.5) 上式可證明前面現象的結果：τ 與 N 間成一線性關係，τ 與 U0間成二次式 的關係，而τ 與 D 間成三次式關係。

5.9 功率係數 C

p

與葉片速度比

TSR 的關係

風機的效能可由功率係數Cp判定，其定義如下： 3 0 1 2 P P C U A

ρ

= (5.6) 上式 P 代表風機所產生的功率 P=

τω

，此處 2 60 N

π

ω

= 代表角速度，將(5.1) 式與(5.6)式比較可得到： P m C =C ⋅TSR (5.7) 將Cm與 TSR 的關係式(5.3)代入可得： 2 ( ) P C = −a TSR + ⋅b TSR (5.8) 如此得到Cp與TSR 間成二次關係式，將此式微分一次並令其為零可得 最佳TSR 值： 1 | 2 optimum b TSR a = (5.9) 此時Cp值為最大

2 max 1 | 4 P b C a = (5.10) 上面結果顯示於圖 5.13，將 a=0.4, b=0.24 代入可得最佳 TSR=0.3 及最大 Cp=0.036。

5.10 擬暫態與全暫態的比較

為了評估擬暫態模擬的準確性，另外也與全暫態(即風輪隨時間旋轉)的 計算做比較，全暫態模擬部分由本實驗室陳虹汝同學提供。圖5.14 與圖 5.15 為入口風速U₀ =8 m s、轉速N =90 RPM 下 Cm 值隨時間的分布，Cm 值 隨時間高低起伏代表風機實際轉動下不同角度產生的扭矩。可以看出在全 暫態計算初期(0 至 5 秒)Cm 值還未達到穩定，隨著時間增加流場漸趨於穩 定，後期(28 至 30 秒)便達到一穩定的週期性變化。 在擬暫態分析中，旋轉角度位置間隔α 為每 30o做一次模擬，由於間隔 角度不夠細小，無法準確預測峰值，而造成扭力過小，故縮小 α 間隔至每 10o做一次模擬，並將不同分析方法做一比較。圖 5.16 及圖 5.17 比較在半 個旋轉週期下，擬暫態(間隔 10o與 30o)及全暫態在不同旋轉角度位置的 Cm 值變化。圖5.16 的轉速 60rpm，可看出擬暫態結果最大 Cm 值落在

α

=30o_以 前，而其他位置角度下擬暫態兩者的線型相近，而在全暫態模擬下 Cm 值 變化較大，且最大Cm 值發生在約

α

=50o_{。圖5.17 為轉速 90rpm 的結果，} 同樣地擬暫態結果最大 Cm 值發生在

α

=30o_{以前，其於角度下線型也非常}

α

=

其餘的計算在此位置附近，Cm 都呈現一小山峰型態的變化，因此兩者有較 大的差異。最後圖5.18 則為轉速 120rpm 的結果，在此轉速下全暫態模擬之 最小Cm 值也接近

α

=120o_{，且Cm 值類似一 sin 函數分布。不過總體而言，} 擬暫態與全暫態模擬在定性上有類似的結果。 圖 5.19 顯示不同計算方法在固定風速U₀ =8 m s下、不同轉速 N 的平 均扭矩輸出，固定位置角度計算的時間步長為0.01 秒，全暫態計算則為 1e-4 秒，由圖可看出將間隔位置角度 α 縮小後扭力值略為提升，但全暫態計算 的扭力值遠超過擬暫態兩者的結果，這可能部分歸因於在擬暫態模擬中流 體的慣性被忽略。另外從

τ

= 的轉速值來看，全暫態模擬中的自由旋轉速0 度約為 240RPM，遠高於擬暫態中的 90RPM，不過從線型來看，除了在低 轉速下，全暫態所得的

τ

與 N 亦近似線性，因此擬暫態的模擬在定性上， 仍有相當程度的可靠度。 利用圖 5.19 可以得到擬暫態(10ｏ_{)與全暫態計算下Cm 與 TSR 間的關係} 係數分別為a=0.48、b= 0.32 及 a=0.326、b= 0.518。圖 5.20 將擬暫態與全暫 態 Cp值做一比較，很明顯的可以看出，全暫態分析的功率係數高出許多。 由前述結果與各種風機之 Cp 與 TSR 關係圖(圖 1.9)比較可知，全暫態分析 之特性與一般 Savonius 型風機類似，且最佳 TSR 及最大功率係數 Cp也與 之相近。

5.11 模擬與實驗數據的比較

在此為了求得流經風機的平均風速，在距風扇 0.7m 處量測不同位置的 風速(如圖 5.21)，其中包含了中心位置(r=0m)、約 r=0.15m 處、邊緣處(r=0.3m)、 外圍(r=0.4m)，每隔 45o量測，將測得數值平均後可得徑向風速分布圖(如圖 5.22)，在此以一二次曲線近似： 2 114 31 6 u= − r + r+ (5.11) 將上式積分後可得平均風速U 如下： ₀ 2 0 0 0 2 R u rdrd U R π

θ

π

=

∫ ∫

(5.12) 其積分上限R=0.4m，經多次實驗量測後，得到平均風速U 約為 5.46 m s 。 ₀ 實驗所得數據經由本節推導過程可得到Cm與 Cp值，圖 5.23 為模擬與 實驗結果之Cm值比較。實驗中風扇與風機之距離即為 0.7m，由此圖可得知， 實驗結果介於擬暫態模擬與全暫態模擬之間，且Cm值皆隨著 TSR 增加而 下降，推測擬暫態模擬低估了扭力輸出，但實驗結果的Cm值與全暫態模擬 有一段落差，這部分有可能為實驗上的功率損失造成的。 圖5.24 則為模擬與實驗結果之 Cp值比較，在此可發現實驗數據最佳Cp 值接近 0.11，整體而言 Cp值皆高於擬暫態分析，但與Cm值相同地遠低於 全暫態模擬。 由(5.9)式與(5.10)式可以整理出不同分析方法與實驗結果的 a、b 係數、 TSR 及最大功率係數 C 為：

固定位置角度α(每 10o) 全暫態分析 實驗結果 a 0.48 0.326 0.62 b 0.32 0.518 0.53 |_optimum TSR _{0.33 0.8 0.43} max | P C 0.053 0.215 0.113 從上表可得知擬暫態與實驗結果的最佳端速比TSR 與最大功率係數 Cp皆低 於全暫態分析，如前述擬暫態分析可能忽略流體慣性與 α 間隔不夠細小導 致低估扭力輸出。而本型式風力發電機葉片尚在初期開發階段中，實驗上 有許多部分急需改善，以下幾點為目前所面臨到的問題，可能為實驗數據 與全暫態模擬結果誤差之主因： • 風力機葉片外型較一般常見葉片複雜，故在製造上較不易，且風機葉 片尺寸大，在精度方面也不足。 • 單層風機葉片模型其中心軸過細，導致承受風吹時中心軸無法有效支 撐葉片重量，容易晃動造成能量損失。 • 承接上點，風力機塔架為單一尺寸，且塔架為針對四層風機設計，無 法與單層風機有效匹配。 • 在風機中心軸與塔架連接部分，無法使用滾珠軸承以降低磨耗，且塔 架多處腐蝕生鏽，導致在實際承受風速下嚴重搖晃。

• 由於風力機尺寸過大，欠缺大型風洞做測試，僅能使用大型工業用風 扇做為風源，而產生的入口風速為帶有旋轉的紊流，且缺乏整流罩將 其風源整流為一均勻流，故無法有效模擬實際在自然風速下的狀況。

第六章 結論

從第五章的結果可歸納出以下幾點結論： 1. 利用固定風機位置的擬暫態模擬計算，其結果顯示平均扭矩與風機轉速 間成線性遞減關係，與風速間成二次的遞增關係，而與風機直徑間成三 次式關係。 2. 將入口風速、風機轉速、風機直徑與輸出扭矩等無因次化，可得端速比 TSR 與平均扭矩係數Cm，上述有因次參數的三種關係式可進而合成單 一的Cm與 TSR 間的線性遞減關係。 3. 由Cm與 TSR 的關係可推出功率係數 Cp與TSR 間呈一二次式關係，進 而可得到相對應的最大Cp值與最佳TSR 值。 4. 擬暫態分析的功率輸出較全暫態分析低，主要可歸因於模擬中流體的慣 性被忽略，並與位置角度α 間隔不夠細小有關。 5. 擬暫態的模擬雖然在定量上與全暫態的模擬及實驗結果有相當的差異， 但上述的Cm及 Cp與TSR 間的線性及二次的關係式仍然適用。 6. 模擬結果與實驗數據間的差異，一重要原因可歸咎於實驗的不準度，由 於不具有大型風洞，只使用風扇作為風源，其風速分佈不均勻且不穩定， 另外實驗裝置精度不足，導致風機整體搖晃振動，喪失了旋轉能量，這 些都是造成實驗誤差的重要因素。

References

1. 馬振基。風能與風力發電技術，臺北市，五南出版公司，2009, 10。 2. Walker J. F., Wind Energy Technology, John Wiley & sons, 1997.

3. 張希良。風力發電技術，新文京開發出版公司，2007。

4. Hau E., Windkraftanlagen: Grundlagen, Technik, Einsatz, Wirtscha- ftlichkeit, Springer, 2003.

5. 美 國 密 蘇 里 州 立 大 學 網 頁 ：

http://extension.missouri.edu/publications/DisplayPub.aspx?P=G1981

6. Patel M. R., Wind and Solar Power Systems, CRC Press LLC, 1999.

7. 工業技術研究院/能源與資源研究所。風力機=Wind Turbine，臺北市， 經濟部能源委員會，1991。

8. Glauert H., “Airplane propellers”, Aerodynamic theory, by Durand W. F., Dover, New York, 1963.

9. Sørensen J. N., Shen W. Z., Munduate X., “Analysis of wake states by a full-field actuator disc model”, Wind Energy, Vol. 1, No. 2, pp. 73-88, December 1998.

10. Mikkelsen R., “Actuator disc methods applied to wind turbines”, Technical University of Denmark, MEK-FM-PHD, 2003.

11. Sørensen J. N., Shen W. Z., “Numerical modeling of wind turbine wakes”, Journal of Fluids Engineering, Vol. 124, No. 2, pp. 393-399, 2002.

12. Troldborg N., Sørensen J. N., Mikkelsen R., “Numerical simulations of wake characteristics of a wind turbine in uniform inflow”, Wind Energy, Vol. 13, No. 1, pp. 86-99, January 2010.

performance of a multi-megawatt wind turbine”, Journal of Physics: Conference Series, Vol 75: 012007, 2007.

14. Mandas N.,  Cambuli F., Carcangiu C. E., “Numerical prediction of

horizontal axis wind turbine flow”, EWEC-2006, Athens, 2006.

15. Hu D., Hua O., Du Z., “A study on stall-delay for horizontal axis wind

turbine”, Renewable Energy, Vol 31, No 6, pp. 821-836, May 2006.

16. Simão Ferreira C. J., Bijl H., van Bussel G., van Kuik G., “Simulating

Dynamic Stall in a 2D VAWT: Modeling strategy, verification and validation with Particle Image Velocimetry data”, Journal of Physics: Conference Series, Vol 75: 012023, 2007.

17. Antheaume S., Maître T.,  Achard J. L., “Hydraulic Darrieus turbines

efficiency for free fluid flow conditions versus power farms conditions”, Renewable Energy, Vol 33, No 10, pp. 2186-2198, October 2008.

18. Howell R. Qin N., Edwards J., Durrani N., “Wind tunnel and numerical

study of a small vertical axis wind turbine”, Renewable Energy, Vol 35, No 2, pp. 412-422, February 2010.

19. Irabu K., Roy J. Na., “Characteristics of wind power on Savonius rotor using a guide-box tunnel”, Experimental Thermal and Fluid Science, Vol 32, pp. 580-586, 2007.

20. Altan B. D., Atılgan M., "An experimental and numerical study on the improvement of the performance of Savonius wind rotor", Energy Conversion and Management, Vol 49, No 16, pp. 3425-3432, 2008.

21. Saha U.K., Thotla S., Maity D., "Optimum design configuration of Savonius rotor through wind tunnel experiments", Journal of Wind Engineering and Industrial Aerodynamics, Vol 96, pp. 1359-1375, 2008.

the performance of Savonius Rotor Based on Numerical Study”, International conference on sustainable power generation and supply, Supergen’09, 2009, Nanjing, 2009.

23. Hu Y. H., Tong Z. M., “The Influence of Windshield on Aero-dynamic Performance of VAWT”, 2009 International Conference on Energy and Environment Technology, iceet, Vol. 1, pp. 893-896, Guilin, China, October 2009.

24. Kamoji M.A., Kedare S.B., Prabhu S.V., “Experimental investigations on single stage modified Savonius rotor”, Applied Energy, Vol 86, pp.1064-1073, 2009.

25. Mohamed M.H., Janiga G., Pap E., Thévenin D., "Optimization of Savonius turbines using an obstacle shielding the returning blade", Renewable Energy, Vol 35, pp. 2618-2626, 2010.

26. Simms D.,  Schreck S.,  Hand M., Fingersh L. J.,  NREL unsteady

aerodynamics experiment in the NASA Ames wind tunnel: a comparison of predictions to measurements, NREL/TP-500-29494, NREL, Golden, Co, 2001.

27. Ebert P. R., Wood D. H., “The near wake of a model horizontal-axis wind turbine-I. experimental arrangements and initial results”, Renewable Energy, Vol. 12, No. 3, pp. 225-243, 1997.

28. Anderson W. K., Bonhus D. L., "An Implicit Upwind Algorithm for Computing Turbulent Flows on Unstructured Grids", Computers Fluids, Vol 23, No 1, pp. 1-21, 1994.

29. Rhie C. M., Chow W. L., "Numerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge Separation", AIAA Journal, Vol 21, No 11, pp.

圖1.1 水平式風力發電機【2】

圖1.3 垂直式阻力型(Savonius)風力機【4】

圖 1.5 垂直式 H 型風力機【4】

圖1.7 改良式垂直風力發電機

圖1.9 各種風機之 Cp 與 TSR 關係圖【7】

圖1.11 帶有擾流板之垂直型風機示意圖【20】

圖1.13 垂直軸風力發電機示意圖【23】

圖2.1 圖 風機旋轉 2.2 流場 轉位置角度 場示意圖 度α

圖3.1 主網格與鄰近網格示意圖 3.2 面上相鄰網格示意圖

c0

c1

r

1 r JG 0 rJG f A JJG

c0

c1

r

i

圖4.1 工業用風扇

圖 4.3 單層風力機葉片

圖4.5 兩種不同尺寸之 coupling

圖4.8 實驗設備示意圖 扭力計 發電機 風輪 可變電阻 顯示器 風扇

U

₀

U

圖5.1 5.2 1 流場計 2 固定位 計算示意圖 位置角度α 圖 α