耦合BEM與CFD方法計算浮體式風機於規則波中運動之研究

(1)

國立臺灣大學工學院工程科學及海洋工程學系碩士論文

Department of Engineering Science and Ocean Engineering College of Engineering

National Taiwan University Master Thesis

耦合 BEM 與 CFD 方法計算浮體式風機於規則波中運動之研究

A Coupled Blade Element Momentum – Computational Fluid Dynamics Model for Evaluating Floating Wind Turbine

Motion under Linear Waves

李仲凱 Chung-Kai Lee

指導教授：郭真祥博士

Advisor: Jen-Shiang Kouh, Dr. Ing.

中華民國 103 年 7 月

July, 2014

(2)

誌謝

承蒙恩師郭真祥博士於碩士班期間的指導，本論文才得以順利完成。感謝老師這段期間的教導與照顧，使學生於碩士班期間瞭解到不論研究或處理事情都該按部就班，遇到問題才可追本溯源進而解決。在此也感謝口試委員蔡進發教授、

蔡國忠教授、趙修武教授對本論文提出的建議與指導，使本論文更加完整。

感謝實驗室的宗岳、淳宇兩位博班學長在學業、生活上提出各自的看法，討論時深刻的學習到處理事情的態度與解決問題的方式。感謝引棋、柏延、育瑋、

子豪、佾臻、怡凱、冠碩、育吟學長姐的照顧、陪伴、傾聽與鼓勵，讓我在實驗室的生活更加充實、多方面的交流帶來許多歡樂。還要感謝子靖、旻玫兩位同學，

雖然已先拋棄你們，但心會與你們同在的。再次謝謝 CAD 大家庭為我的碩士生涯增添許多美好的回憶，讓我體會到求學中不僅要努力學習也不忘適當的放鬆娛樂。

最後更要感謝父母、叔叔與其他家人們，你們的支持與無所不在的陪伴讓我更有動力、更堅信能夠完成碩士班學業，最後也畫下了完美的句點。進入職涯後不會辜負大家對我的栽培與期望。

(3)

I

摘要

浮體式風機所搭配之浮體平台對其受風與波浪耦合作用下之運動程度扮演重要的角色，在離岸風場位址選擇受陸上地形及風切影響較小之浮體平台有利於減少風機之發電量損失。本研究以 NREL 5MW 風機當作目標風機，並分別搭配 OC3-Hywind 之 Spar 型浮體平台與 OC4 DeepCwind 之半潛型浮體平台為計算模擬之對象。研究方法為雷諾平均納維爾史托克(RANS)求解浮體平台受水動力作用之流場並選用適當的紊流模型及額外搭配 BEM 程式估算風機之氣動力性能。以 BEM 進行目標風機之氣動力性能計算結果與 NREL 公佈之資料相互驗證，轉子功率計算結果與資料趨勢相符，額定風速(11.4m/s)下之差異約為 7.64%。BEM 程式於縱搖運動下修正之計算結果顯示考慮真實風機之控制情形其平均功率會下降，且風機推力與功率變動週期與縱搖角速度週期相同，與參考文獻以 RANS 計算之結論相符。耦合 BEM 與 CFD 軟體計算浮體式風機於波高 4m 週期 10s 之規則波中進行縱搖運動之數值模擬，並比較 Spar 型與半潛型浮體式風機之實際發電平均功率耗損率以半潛型浮體式風機之耗損率較小，此結果與文獻亦相同，顯示耦合 BEM 與 CFD 軟體計算浮體式風機受風、波耦合作用下之可行性。藉由不同的波浪參數探討以發電量為目標時 Spar 型浮體式風機與半潛型浮體式風機之適用範圍，由數值模擬結果顯示當波浪週期小於 8.4s 之線性波，將選則 Spar 型浮體式風機；當波浪週期大於 10s 與波浪近似於靜水狀態則選擇半潛型浮體式風機能夠減少功率損失。

關鍵字：葉片元素動量理論、計算流體力學、離岸浮體式風機、Spar 型浮體平台、

半潛型浮體平台

(4)

II

Abstract

The offshore wind turbine will be influenced by the couple of aerodynamics and hydrodynamics.The selection of floating platform will directly affect the generating power of the wind turbine. The floating system of the OC3-Hywind and the OC4 DeepCwind is applied in this paper. Through the Blade Element Momentum theory evaluate the power curve of NREL-5MW wind turbine and compare to the NREL data.

It shows two of them have similar power curve, but they have about 7.64% difference in the rating velocity condition. When we consider the real wind turbine control system in the pitch motion, the average power is no longer that high anymore. Within the motion of turbine, it causes the power of turbine have the same period as the motion of turbine.

In this paper, coupling the Blade Element Momentum (BEM) and computational Fluid Dynamics model (CFD) to evaluate two type of floating wind turbine, spar and semisubmersible wind turbine. Under the wave height 4m and period 10s, it shows that the semisubmersible type floating wind turbine has lower power loss than the spar type of wind turbine. This simulation case is verified with the reference.

According to this research present a suitable wave condition to reduce the power loss for both of the floating wind turbine. If the wave period is less than 8.4s, so the spar type floating wind turbine has been recommended. But if the wave period is larger than 10 s and the wave is just like calm sea state, the semisubmersible floating wind turbine will be better.

Keywords：Blade Element momentum theory, CFD, Offshore wind turbine, Spar platform, Semisubmersible platform

(5)

III

圖目錄

圖 1-1 全球再生能源統計[1] ... 1

圖 1-2 歐洲市場風能之新裝置容量預測圖 ... 2

圖 1-3 全球離岸式風機裝置容量 ... 2

圖 1-4 丹麥 Vindeby 離岸風場[4] ... 3

圖 1-5 未來 40 年風力發電裝置容量預測[5] ... 4

圖 1-6 淺水區域之基礎型式[6] ... 5

圖 1-7 過渡水域基礎應用型式(左為三腳架，右為套管式)[7] ... 5

圖 1-8 深水區域浮體平台型式[6] ... 6

圖 1-9 離岸風場之離岸距離與水深分佈[8] ... 7

圖 1-10 浮體式風機穩定型式示意圖[9] ... 7

圖 1-11 浮體平台穩定策略分析圖[9] ... 8

圖 2-1 紊流無因次化壁面距離與速度分佈圖[21] ... 16

圖 2-2 計算網格之 y^值 ... 18

圖 2-3 體積分率說明 ... 19

圖 2-4 軸向動量理論示意圖 ... 21

圖 2-5 旋轉跡流示意圖 ... 24

圖 2-6 轉子圓盤環狀面積力矩 ... 25

圖 2-7 葉片元素理論示意圖 ... 26

圖 2-8 二維翼型速力圖 ... 26

圖 2-9 浮體式風機六自由度運動示意圖 ... 33

圖 2-10 浮體式風機世界座標與局部座標系 ... 34

圖 2-11 懸鍊線示意圖[19] ... 36

圖 2-12 Spar 型浮體平台與錨鍊系統示意圖 ... 38

(9)

VII

圖 2-13 繫錨點和繫錨力說明[19] ... 39

圖 2-14 浮體式風機座標系統 ... 41

圖 2-15 BEM 對縱搖運動修正之參數示意圖 ... 42

圖 3-1 NREL 5MW 風機葉片幾何參數定義 ... 47

圖 3-2 NREL 各剖面翼型幾何 ... 47

圖 3-3 DU40 翼型之氣動力性能 ... 48

圖 3-8 NACA64 翼型之氣動力性能 ... 50

圖 3-9 BEM 程式計算單一葉片剖面之流程圖 ... 51

圖 3-10 NREL 5-MW 風機之功率曲線與推力曲線[29] ... 53

圖 3-11 BEM 與 NREL 之功率比較 ... 54

圖 3-12 BEM 與 NREL 之轉子推力比較 ... 54

圖 3-13 額定風速下翼剖面之攻角分佈圖 ... 55

圖 3-14 OC3-Hywind 浮體式風機示意圖[32] ... 57

圖 3-15 Spar 浮體平台幾何模型 ... 58

圖 3-16 DeepCwind Semisubmersible 浮體平台之幾何配置[33]... 60

圖 3-17 半潛型浮體平台模型示意圖 ... 60

圖 3-18 OC4 DeepCwind Semisubmersible 錨鍊系統之繫錨鍊配置圖 ... 61

圖 3-19 圓形數值模擬之流場範圍與邊界條件示意圖[34] ... 62

圖 3-20 圓形前方局部加密區域 ... 64

圖 3-21 圓形之阻力係數與雷諾數之關係[35] ... 65

圖 3-22 圓形邊界網格獨立性測試 ... 66

(10)

VIII

圖 3-23 滑移壁面(Slip Wall)邊界驗證 ... 67

圖 3-24 壓力出口之全壓驗證 ... 67

圖 3-25 一維波浪方程特徵線影響區域示意圖 ... 68

圖 3-26 自由液面區域網格佈置圖 ... 69

圖 3-27 純波浪驗證觀測點位置(0,150) ... 70

圖 3-28 Wave0 觀測點高度與時間之關係圖 ... 70

圖 4-1 BEM 對縱搖運動修正之參數示意圖 ... 73

圖 4-2 不同縱搖角速度振幅之風機性能計算 ... 74

圖 4-3 BEM 計算風機於 0.005sin(² ) p 10 ω  πt 之性能 ... 75

圖 4-4 OC3-Hywind 錨鍊系統之繫錨鍊配置圖 ... 76

圖 4-5 Spar 浮體平台之錨鍊示意圖 ... 77

圖 4-6 Spar 型浮體式風機錨鍊 1 於初始狀態之參數值 ... 78

圖 4-7 OC4 DeepCwind Semisubmersible 錨鍊系統之繫錨鍊配置圖 ... 80

圖 4-8 半潛型浮體式風機錨鍊 1 於初始狀態之參數值 ... 81

圖 4-9 BEM 與 RANS 耦合計算流程圖 ... 83

圖 4-10 浮體式風機流場示意圖 ... 83

圖 4-11 浮體式風機之網格佈置策略 ... 85

圖 4-12 半潛型風機計算邊界條件驗證圖 ... 86

圖 4-13 半潛型風機計算觀測點高度隨時間變化圖 ... 86

圖 4-14 半潛型風機計算壓力出口邊界條件驗證圖 ... 87

圖 4-15 Spar 型浮體風機於H4m，T 10 s 之縱搖角速度與時間之關係圖 ... 88

(11)

IX

圖 4-16 Spar 型浮體式風機於H5.2m，T10s 之之縱搖角速度與時間之關係圖 ... 89 圖 4-17 Spar 型浮體風機於 H 3.6 m，T8.4s 之之縱搖角速度與時間之關係圖 ... 90 圖 4-18 半潛型浮體式風機於H4m，T 10 s 縱搖角速度與時間之關係圖 ... 92 圖 4-19 半潛型浮體式風機於H5.2m，T 10 s 縱搖角速度與時間之關係圖 ... 93 圖 4-20 半潛型浮體式風機於H3.6m，T8.4s 縱搖角速度與時間之關係圖.... 94 圖 4-21 Spar 型浮體式風機於波高 3.6m，週期 8.4s 力矩分析 ... 99 圖 4-22 半潛型浮體式風機於波高 3.6m，週期 8.4s 力矩分析 ... 99 圖 4-23 Spar 型浮體風機於波高 5.2m 週期 10s 下縱搖角速度振幅與頻率關係圖 ... 101 圖 4- 24 Spar 型浮體風機於波高 3.6m 週期 8.4s 下縱搖角速度振幅與頻率關係圖 ... 101 圖 4-25 半潛型浮體風機於波高 5.2m 週期 10s 下縱搖角速度振幅與頻率關係圖 ... 102 圖 4-26 半潛型浮體風機於波高 3.6m 週期 8.4s 下縱搖角速度振幅與頻率關係圖 ... 102 圖 4-27 浮體式風機適用分佈圖 ... 104

(12)

X

表目錄

表 3-1 風機等級分類[31] ... 44

表 3-2 NREL-5MW 離岸風機基本參數[29] ... 45

表 3-3 NREL-5MW 葉片結構參數[29] ... 45

表 3-4 NREL-5MW 葉片幾何參數[29] ... 46

表 3-5 風速與葉片節距角之關係 ... 53

表 3-6 BEM 與 NREL 功率之比較 ... 55

表 3-7 葉片半徑 0.7R 之翼型雷諾數 ... 56

表 3-8 OC3-Hywind 浮體平台結構參數[32] ... 57

表 3-9 OC3-Hywind 錨鍊系統參數[32] ... 58

表 3-10 DeepCwind Semisubmersible 之浮體平台幾何與結構參數 ... 59

表 3-11 OC4 DeepCwind Semisubmersible 浮體式風機之錨鍊系統 ... 61

表 3-12 圓形邊界網格獨立性測試 ... 66

表 3-13 純波浪的時間步進間隔 ... 69

表 4-1 縱搖運動下之計算組數 ... 73

表 4-2 OC3-Hywind 錨鍊系統參數[32] ... 76

表 4-3 Spar 型浮體式風機於初始狀態之錨鍊系統狀態 ... 79

表 4-4 OC4 DeepCwind Semisubmersible 浮體式風機之錨鍊系統 ... 80

表 4-5 半潛型浮體式風機於初始狀態之錨鍊系統狀態 ... 81

表 4-6 浮體式風力機計算之波浪參數 ... 82

表 4-7 浮體式風機計算域空間 ... 84

表 4-8 邊界條件設定 ... 85

表 4-9 Spar 型浮體風機於H4m，T 10 s 之計算結果 ... 88

表 4-10 Spar 型浮體風機於H5.2m，T10s 之計算結果 ... 89

(13)

XI

表 4-11 Spar 型浮體風機於 H 3.6 m，T8.4s 之計算結果 ... 91

表 4-12 半潛型浮體式風機於H4m，T 10 s 之計算結果 ... 92

表 4-13 半潛型浮體式風機於H5.2m，T 10 s 之計算結果 ... 93

表 4-14 半潛型浮體式風機於H3.6m，T8.4s 之計算結果... 95

表 4-15 浮體式風機於縱搖運動下模擬結果比較 ... 96

表 4-16 各波浪參數之模擬結果比較 ... 98

(14)

XII

符號表

符號單位描述符號單位描述

P Pa 壓力 C _D - 阻力係數

 Pa-s 流體動力黏滯係數 c m 翼型弦長

ρ kg/m³ 密度 V m/s 剛體運動速度向量

T N 推力 ω rad/s 剛體旋轉速度向量

U m/s 流速 H m 波高

P atm Pa 大氣壓力 T s 週期

m kg/s 質量流率 λ m 波長

a - 軸向誘導速度因子 Φ degree 繫錨點連結角度

C T - 推力係數 T _H N 懸鍊線水平力

P w Watt/s 功率 T _V N 懸鍊線垂直力

C P - 功率係數 U _in m/s 相對入流速度向量 Q m²/s² 扭矩 ω_p rad/s 縱搖角速度

Ω rpm 轉子轉速 _θ degree 葉片方位角

I kg-m² 轉動慣量 ψ _wt degree 風機方位角 ω rpm 旋向誘導轉速 c _w m/s 波浪群速度

'

a - 旋向誘導速度因子 θp degree 縱搖角度

φ degree 入流角 θ_p degree 平均縱搖角度

α degree 攻角 δp - 功率增減率

r m 翼斷面之半徑位置 δT - 推力增減率

(15)

XIII

β degree 總節距角 δp - 平均功率增減率

θ T degree 扭轉角 δT - 平均推力增減率

,0

θp degree 節距角 δp_r - 以額定功率為上限

之功率增減率

C L - 升力係數 δp_r - 以額定功率為上限

之平均功率增減率

(16)

1

第1章緒論

1.1 研究背景

從工業革命開始，科技逐漸進步使得全球對於能源的需求日益提升，並且使用傳統石化能源與火力發電會排放大量的二氧化碳，進而導致溫室效應與全球暖化之現象，故為了減少二氧化碳之排放，聯合國於 1997 年邀約參加國簽訂《京都議定書》(Kyoto Protocol)，目的為將大氣中的溫室氣體含量穩定於適當之水平，防止氣候產生劇烈變化對人類造成威脅，所以為了因應協定各國逐漸投入再生能源與替代能源的研究，如太陽能、潮汐能、波浪發電、風力發電、核能…等。其中風力發電之建造成本逐漸降低，使風力發電為目前發展最佳之再生能源。由 21 世紀再生能源政策網(Renewable Energy Policy Network for the 21^st Century)可知全球風力發電機之裝機容量於 2006 年至 2011 年間每年之容量成長率均高於 20%，並且至 2012 年全球總容量已達 282GW[1]，且約佔了總再生能源發電量之 59%，如圖 1-1 所示。

圖 1-1 全球再生能源統計[1]

根據 DNV．GL 之市場估計，風能主要市場如英國、德國、丹麥等歐洲各國每年總裝置容量將逐年上升，顯示了風能之未來發展，如圖 1-2 所示。

(17)

2

圖 1-2 歐洲市場風能之新裝置容量預測圖

雖然風力發電之發展快速，但全球陸上風場有限，因此為了繼續增加風力發電之發電量，必須將風場擴展至海洋區域發展離岸式風機。離岸式風機依照其裝置之水深可分成固定式風機與浮體式風機兩種，其中固定式風機均設於近岸且淺水區，並以打樁方式將平台固定於海床上，適用水深約為 85m 以內；而浮體式風機則較不受水深之影響，可裝置於水深 100m 以上之深水海域。全球發展離岸風機之國家以丹麥、德國、英國與荷蘭為先驅，並已有許多離岸風場進行營運。至 2014 年底，全球之離岸風場裝機容量約為 4.7GW，如圖 1-3 所示，其中歐洲占 68%的總離岸風力裝置量。

圖 1-3 全球離岸式風機裝置容量

台灣西部沿岸與澎湖地區擁有豐富之風能，因此適合發展風力發電，但台灣西岸適合開發風場之區域已逐漸飽和，故台灣當局亦開始將風能推廣至離岸發電。

經濟部能源局於 101 年 7 月 3 日公告實施「風力發電離岸系統示範獎勵辦法」，

(18)

3

輔助二家民營業者與一家國營業者建立三個示範風場，並要求業者於民國 104 年完成示範機組之建置、測試與竣工，並取得電業執照開始商業運轉，以協助達成政府所提出的千架風力發電機的政策目標。

目前固定式風機之打樁技術已相當成熟，因此現階段離岸風電研究著重於深水區之浮體式風機，浮體式風機可安裝於水深 100m 以上之區域，因此風場可架設範圍更廣大，可提升風力發電的發電量。浮體式風機結構上為風機與浮體平台之結合，故有別於固定式風機，因此需更多的分析與測試。

1.2 離岸型風力發電介紹與浮體式風機介紹

1970 年代麻州大學的 Heronemus 教授提出浮體式風機的概念[2]，但當時再生能源尚未有發展之急迫性，因此風力發電技術尚未成熟，故此概念不受重視。直到 1990 年，爆發第三次石油危機，歐洲各國企業投入風機產業，使風力發電開始蓬勃發展，技術亦逐漸成熟，因此離岸風電之概念再度被提起。1990 年”World Wind”

公司建立了第一座離岸式風機，該風機位於瑞典北方離岸 250m、水深 7m 之海域；

隔年，位於丹麥的 Vindeby 離岸風場正式完工，共裝置 11 座 0.45MW 發電量之風機，總裝置發電量為 4.95MW[3]。

圖 1-4 丹麥 Vindeby 離岸風場[4]

(19)

4

隨後歐洲各國陸續訂定相關發展計畫，離岸風力發電開始快速發展，並且推廣至世界各地，根據國際能源署(International Energy Agency, IEA)評估，未來 40 年離岸風力發電將呈現大幅成長，如圖 1-5 所示。

圖 1-5 未來 40 年風力發電裝置容量預測[5]

由於離岸風電有別於以往陸上風機，並且擁有下面幾個特點：

(1) 陸上風能易受周遭建築物與地形影響，海上則無，因此風能更為穩定。

(2) 海域廣闊，能開發較大型的風場。

由於上述之優點，各國紛紛開始建置離岸型風場。而離岸型風機體需搭配平台，平台的基礎型式依水深可分成下面三類：

(1) 淺水區域(Shallow water)基礎。

(2) 過渡區域(Transition water)基礎。

(3) 深水區域(Deep water)基礎。

其中淺水區域適用之水深約為 0m~50m，基礎型式有兩種：重力式(Gravity Base) 與單樁式(Monopile)。單樁式基礎以打樁方式將風機固定於海床上，故建置時須與海事工程進行搭配，而重力式則利用自身之重量將之進行故定，淺水域型式如圖 1-6 所示。

(20)

5

單樁式基礎其缺點為基礎打入海床後，由於海流流經基樁時於基樁後方會產生渦流，因此易產生流砂進而掏空基座之打樁深度，故須在海床之外圍鋪砂石以減輕掏空之情況。而重力式基礎之缺點為易受海床成分的影響，並且對波浪方向性相當敏感。

圖 1-6 淺水區域之基礎型式[6]

過渡水域適用之水深約為 30m~100m，型式如圖 1-7 所示，左圖為三角架型式 (Tripod)，右圖為多管式(Jacket)。

圖 1-7 過渡水域基礎應用型式(左為三腳架，右為套管式)[7]

設計多管式基礎時，必須考慮其自然頻率問題，由於容易受到碎波的影響，

故須將自然頻率設計至較低頻之區域，但其基礎於海床上較不需考慮流砂的問題。

(21)

6

表 1-1 淺水區域與過渡區域基礎型式比較

Type Typical water depth Typical Size Typical Weight Gravity ~20m  30m 1000~5000t Monopile ~35m  4~5m 600~700t

Tripod ~50m ~35x60m 1000t

Jacket ~60m ~25x60m 700~900t

深水區域適用於超過 60m 之深水域，應用的基礎型式為浮體平台並使用錨鍊系統固定浮體式風機之運動範圍，如圖 1-8 所示。共分四類：Semisubmersible 型、

Barge 型、Spar 型與 Tension-leg platform(TLP)型。圖中兩種 Spar 型的差異在於錨鍊系統型式，一為張錨式(右一)，另一為懸垂式(左三)。

圖 1-8 深水區域浮體平台型式[6]

目前各國致力於離岸風場之開發，但由於海岸線有限，因此淺水與近岸區域已逐漸開發完成，故近年來離岸風場的發展趨勢為離岸漸遠、水域較深之區域擴展，圖 1-9 為歐洲風能協會(European Wind Energy Association, EWEA)統計至 2012 年離岸風場之離岸距離與水深分佈，藍色為已建置完成，紅色為建造中，綠色為審批中之項目。從審批中的項目分佈來看，可知離岸風場之開發區域之水深和離岸距離正逐漸增加。

(22)

7

圖 1-9 離岸風場之離岸距離與水深分佈[8]

浮體式風機是由兩個部分組成，第一部分為與傳統固定式風力機主體相同，

包括轉子葉片(Rotor)、輪轂(Hub)、機艙(Nacelle)以及塔柱(Tower)，第二部分是浮體平台，用於提供浮力以支撐風力機，另外搭配適當的錨鍊系統用於固定風力機，

避免風機受潮水作用漂走。美國國家能源實驗室 (National Renewable Energy Laboratory ,NREL)根據使浮體風機達成穩定之策略以及錨鍊系統分為三類[9]，如圖 1-10 所示。

圖 1-10 浮體式風機穩定型式示意圖[9]

(23)

8

圖 1-10 由左至右分別為壓載平衡式(Ballast Stabilized)、繫纜線平衡式(Mooring Line Stabilized)以及浮力平衡式(Buoyancy Stabilized)，下面將針對此三種平衡模式進行說明：

(1) 壓載平衡式(Ballast Stabilized)：此平衡方式為利用壓載使重心降低至浮心以下的位置，以產生較大之扶正力臂與慣性回復力，因其幾何較為細長因此可得較高之慣性矩以減少縱搖與橫搖運動的幅度，並利用懸垂式錨鍊(catenary)加以固定系統，Spar 型浮體式風機為此種型式。

(2) 繫纜線平衡式(Mooring Line Stabilized)：此種方式之錨鍊系統是垂直固定於海床上以平衡浮體平台，並利用較高張力的錨鍊加以固定。Tension Leg platform (TLP)為此種系統，目前大型海上鑽油平台多也使用此種型式。

(3) 浮力平衡式(Buoyancy Stabilized)：平衡方式為利用較大的水線面積與平均分佈之浮力以產生扶正力矩，再使用懸垂式錨鍊固定結構物。Barge 型浮體平台具代表性，Semisubmersible 型也類似此種型式。

一般而言，設計浮體平台時均考量到浮力、錨鍊、壓載三種模式，但各模式所佔比例不同，因此產生不同型式的浮體平台，策略分析如圖 1-11 所示。

圖 1-11 浮體平台穩定策略分析圖[9]

1.3 文獻回顧

繼 1970 年代麻州大學的 Heronemus 教授提出浮體式風機概念[2]及 1990 年”World Wind”公司建立了第一座離岸式風機後，初期大多關注於近岸式風力發電，

(24)

9

直至 2000 年左右浮體式風力機之相關研究才逐漸發表。浮體式風機為風力機主體搭配浮體平台，因此對於要選擇何種浮體平台成為重要的研究課題。起初，浮體平台的設計概念均來自海上鑽油平台，各國紛紛提出不同的浮體平台型式。1998 年英國的 Tong 提出搭配 Spar 型浮體平台的”FLOAT”，發電容量為 1.4MW；同年英國的 Henderson 和 Patel 提出 Semisubmersible 型式的浮體平台設計概念[10]；2002 年 Lagerwey the Windmaster 公司提出 Dutch Tri-floater 的浮體式風機設計[11]，亦為 Semisubmersible 型的浮體平台；2004 年 NREL 的 Musial 等人提出 TLP 型浮體平台之設計，上方承載 5MW 風機，並比較 TLP 型與 Dutch Tri-floater 型平台[12]；

同年，日本的 T. Kogaki 提出 Spar 型浮體平台的設計概念[13]；2005 年 NREL 根據浮體式風機的固定與平衡方式將浮體平台分為三類，並且分析個別的優缺點[9]。

在 2000 年以後，浮體式風機的分析工具更加完善，2007 年 Jonkman 等人提出耦合時域分析方法，並對 5MW 發電量之 Barge 型浮體式風機進行分析，此耦合時域之分析方法考量了空氣動力、水動力、結構、錨鍊以及控制部分[14]。

除了數值分析工具外，許多模型試驗也被進行，如 2006 年 Nielson 等人進行 Hywind 浮體式風力機的模型試驗[13]；Tomoaki 於日本五島列島進行 1:10 之 Spar 型浮體式風機模型試驗，他們的模型試驗結果與模擬結果進行比較後驗證了這些分析的可靠性[15]。

1.4 研究目的與方法

浮體式風機由於遭受風、波浪與海流與繫纜等作用，因此須考慮氣動力、水動力與繫纜力三部分，本研究基於 IEC61400-3 設計負荷條件(Design Loading Case ,DLC)2.4 進行計算，其中僅考慮空氣、波與繫纜三種受力形式，且將風與波浪設為同向。

風機之風力計算使用葉片元素動量理論(Blade Element Momentum, BEM)進行求解轉子之推力與發電量。BEM 方法為常用的簡易於評估風機性能之方法，如 2007 年 NREL 開發的 FAST 模擬程式之氣動力的基礎亦是使用 BEM 方法，並加入了

(25)

10

Prandtl 的葉尖損失(Tip loss)[16]與輪轂損失(Hub loss)，以及 Beddoes-Leishman 的動態失速(dynamic stall)模型[14]；2012 年 R. Malki 等人使用 BEM 方法與數值模擬軟體進行潮流發電機之性能評估[17]；同年，Céline Faudot 等人使用 BEM 耦合波浪力方法計算潮流發電機葉片之受力[18]。顯示 BEM 方法於計算風或流動力效能，

被廣為應用。並且 BEM 方法相較於 CFD 軟體求解風機氣動力性能極具時間之經濟性。

水動力過去則利用勢流理論計算，忽略流體的黏滯性，並且假設流場為無旋性與不可壓縮之理想流體，將統御方程式簡化為拉普拉斯(Laplace)方程，再藉由邊界元素法(Boundary Element Method)或稱小板法(Panel Method)進行離散，並於小板上分佈渦流(Vortex)、源流(Source)、偶流(dipole)或匯流(Sink)進行流場計算，此種方法在過去常應用於船舶與海洋結構物之計算。但以勢流理論求解波浪力負荷時，

為簡化問題的求解多必須假設物體為靜止或物體運動相對於波浪極小，然而由研究可知當浮體式風機受外在負荷時可能產生較大之旋轉角度，故在某些特定情況下並不適用。水動力問題亦可使用 Morison 方程求解，但限用於波浪之波長大於圓柱直徑 5 倍的情況，且無法直接估算阻力。因此本研究使用計算流體力學求解浮體平台受波浪作用下之黏性流問題，其控制方程式為雷諾平均化那維爾-史托克 (Reynolds-Averaged Navier-Stokes, RANS)方程式，求解時先對空間域進行離散，離散方法為有限體積法(Finite Volume Method)，並搭配適當之紊流模型求解黏性流流場，此方法不僅可解析波浪力，在阻力部分亦可直接求解，並選用適當之網格佈置可增加其準確性。

本研究耦合 BEM 方法與計算流體力學(CFD)之套裝軟體 Star CCM+以模擬容量 5MW 之 Spar 型浮體式風機與 Semisubmersible 型浮體式風機，並搭配懸垂式錨鍊系統作為研究對象，比較兩種型式之運動型態差異與其適用海域在波浪中之運動與風機性能。

(26)

11

1.5 本文架構

本論文共分五個章節，各章節內容如下述。

第 1 章緒論：介紹研究背景、離岸風機與浮體式風機、文獻回顧以及研究目的與方法。

第 2 章理論基礎：針對計算流體力學使用之理論與統御方程式作概略性之介紹，

並簡介風機效能計算方法、對風機於縱搖運動下進行修正之推導。

第 3 章浮體式風機幾何與數值驗證：進行目標風機效能計算之驗證以及浮體平台之幾何、圓形阻力模擬與驗證、模擬流場範圍與邊界條件、純波浪傳遞驗證、網格佈置策略。

第 4 章浮體式風機數值計算：內容包括計算風機於週期性縱搖角度變化下之性能與耦合 BEM 與 CFD 方法進行浮體式風機之數值模擬，並與文獻驗證方法可行性，

接著探討浮體平台運動之差異。

第 5 章結論：針對第 4 章模擬結果進行浮體式風機於不同波浪參數下之運動特性與功率耗損交叉分析。

(27)

12

第2章理論基礎

2.1 計算流體力學

本節介紹流體力學之理論基礎-統御方程式，包括連續方程式、動量守恆方程式，

以及計算流體力學之紊流模型、壁面函數、數值方法以及自由液面。

2.1.1 統御方程式

使用數值方法求解流場問題時，主要為求解流體力學的控制方程式。流體流動時必須遵守物理定律，此定律即為建立控制方程的依據。最基本的兩個統御方程式為連續方程式以及動量守恆方程式。

連續方程式的基本物理意義為在單位時間內質量流入控制體積內與流出的差值即為該控制體積內部增加的質量，方程式如下：

(2-1) 其中為密度，v = ev_i _i i1, 2,3，e₁,e₂,e₃為卡氏座標 x,y,z 軸的單位向量。

動量守恆方程式的物理意義為控制體積內動量隨時間的變化量等於作用在上之外力總和。

動量守恆方程式如下：

2

( 2^j)

i i

j i i

j i j

v v P v

ρ v g F

t x x x

        

      i j, 1, 2,3 (2-2)

其中 P 為壓力，為流體動力黏滯係數，F 為作用在控制體積上的外力。本研究_i

僅考慮重力加速度場並定義垂直方向向上為正，故g₁g₂ 0，g₃ 9.81 /m s²。本研究之數值方法是藉由有限體積法將空間離散，再求解三維且不可壓縮之雷諾平均化納維爾-史托克方程式(Reynolds Averaged Navier-Stokes 簡稱 RANS)。

(28)

13

其統御方程式中之連續方程式與動量守恆方程式以雷諾平均法處理過程如下：

將速度及壓力分成平均量及變動量：

v v v' (2-3)

  '

P P P (2-4)

接著將其速度及壓力採取時間平均量的觀念，可知變動量對時間平均後其值為零，

並將(2-3)與(2-4)代入連續方程式與動量守恆方程式，並假設為不可壓縮流，因此連續方程式與動量守恆方程式以雷諾平均法處理後表示如下：

連續方程式：

, 0

vi i  (2-5)

動量方程式：

2

( ' ' )

( ^j) ^j ⁱ

i i

j i i

j i j j

v v v

v v P

ρ v g F

t x x x x

 

        

       i, j 1, 2,3 (2-6)

式(2-5)、式(2-6)即為 RANS 方程式，其中式(2-6)等號右邊第三項稱為雷諾應力項 (Reynolds stresses)。

2.1.2 紊流模型

使用雷諾平均法處理統御方程式時會產生待解未知數 v' ，即為雷諾應力項。

為了求解此非線性項，會利用紊流模型進行求解。紊流模型(Turbulence Model)是以半經驗的方式建立紊流流場中流體分子運動的模式，在許多工程問題上，都已被證實具有一定的正確性與可信度，故已被廣泛使用。

本研究以數值方法處理風機受水動力作用之性能，採用之方法沿用文獻[19]

使用之紊流模型，該文獻使用 Star CCM+提供的 SST k紊流模型計算海上浮體風力發電機之浮體運動以及風力發電機性能，在風力發電機性能計算以及浮體運

(29)

14

動及二維圓柱阻力均顯示不錯的準確性。故本研究選用 SST k紊流模型作為計算之基礎。

有關 k紊流模型之細節可參考文獻[20]，相關之數學關係式如下：

紊流動能方程式：

( ) ( _j ) _ij ⁱ ^* ( _k _t)

j j j j

v k

ρk ρv k τ β ρωk + μ+σ μ

t x x x x

 



        

       (2-7)

紊流動能消散率方程式：



1



2

( ) ( _j ) _ij ⁱ ² ( _{ω t}) 2 1 _ω 1

j t j j j j j

v ω k ω

ρω ρv ω τ βρω + μ+σ μ F ρσ

t x x x x x x

 



            

          

(2-8)

其中μ_t為渦旋黏滯係數

1

1 2

max( ; )

t

μ a k

a F

   (2-9)

y

  



v ，常數a₁0.31

混合函數 F₂為壁面距離 y 之函數

2 2 tanh(arg )2

F  (2-10)

2 2

arg max(2 ;500 ) 0.09

k

ωy y ω

  (2-11)

式(2-8)中的 ₍ _{ω t}₎

j j

μ+σ μ ω

x x

 

 

 

    為正交擴散項，此為將kε模型轉換為 k ω 模型 時的產生項，而函數 F1考量邊界層區域的 k ω 方程之係數以及自由剪切層和自由 流區域的 k ε 方程之係數，函數 F₁在無滑移壁面上時為 1，在邊界層內則趨近於 1，邊界層邊緣時則為 0，混合函數 F1公式如下：

4 1 tanh(arg )1

F  (2-12)

(30)

15

2

1 * 2 2

4

arg min max ;500 ; ^ω

kω

ρσ k k

β ωy y ω CD y

   

     (2-13)

10 2

max(2 1 ;10 )

k ω

j j

k ω

CD ρσ

x x

 

  

   (2-14)

式(2-7)中 SST k紊流模型常數

* 0.09

β  0.41 k

另外 SST k紊流模型中的係數σ σ β γ_k, _ω, , 可以表示，此紊流模型為 Standard k與 Standard k ε 兩者特性之結合，故可將分解為_及_，其中_代表 Standard k模型的係數，_為 Standard k ε 模型轉換後的係數，如下所示：

1 1

{σ σ β γ_ω, _k, , } F _ (1 F) _

   

    (2-15)

第一組係數_：

1 1 1

2 1

1 1

0.5, 0.85, 0.075 0.553

ω k

ω

* *

σ σ β

β k

γ = σ

β β

  

  (2-16)

第二組係數_：

2 2

1.0, 0.856, 0.0828 0.440

ω2 k

2 * ω2 *

σ σ β

β k

γ = σ

β β

  

  (2-17)

2.1.3 壁面函數

當流體流經無滑移邊界(No-slip condition)時，因受黏性效應，所以流體在靠近固體邊界時流速會急速衰減至與邊界相等的速度，造成壁面區域的流速具有相當

(31)

16

大的梯度。為了解析流體的流速在此區域的快速變化，因此在固體邊界附近的網格必須佈置得較為密集。

壁面函數法(Wall Function)是由於解析壁面剪應力 u τ μ

y

 

 之速度梯度時須使用極細密網格才會達到準確性。若使用半經驗速度分佈函數(Velocity profile)即能降低總網格數。

圖 2-1 紊流無因次化壁面距離與速度分佈圖[21]

壁面函數相關公式如下：

*

u u u

  (2-18)

yρu*

y μ

   (2-19)

其中 u⁺為無因次化的速度參數，u 為流場平均流速，u^*為摩擦速度，y⁺為無因次後的壁面距離，物理意義為局部網格的雷諾數，y為壁面與壁面法線方向第一個網格大小，μ為流體黏滯係數，ρ為流體密度。由實驗結果得知，當0y^ 5時u^與

y^的區域稱為黏性次層流，關係如下：

(32)

17

*

u yu

u y

u ν

    (2-20)

當30 y^ 800，u^與y^為對數關係，公式如下：

1 ln(y ) 5

u K

    (2-21)

其中 K 為 Von Karman 常數，K0.41，當5 y^ 30 視為過渡區(Transition region)，

此範圍內沒有公式可以表示。

求解不同的局部網格雷諾數 y^時有其對應之璧面函數方法，可分為下列三種：

(1) High y^ Wall treatment：壁面之第一層層狀網格大小需佈置於圖 2-1 之對數區域，即y^ 30。

(2) Low y^ Wall treatment：此區域之第一層層狀網格大小需極為精細， y^需小於 1 才可適用。

(3) All y^ Wall treatment：此方法處理壁面函數時，會依照每個網格的 y^選擇合適之壁面函數進行求解，但此種方法相較於(1)、(2)其精準度較差。

本研究選用 All y^ Wall treatment 作為初步設定，二維圓型計算結果發現前 方之停滯點與分離後之區域易求解出較低之 y^，而自由液面區域之流體速度較快

以相同之層狀網格設定易出現較高之 y^，不適用(1)、(2)之設定，故選用 All y^ Wall treatment 進行計算。圖 2-2 為本研究經由網格獨立性測試後所採用的網格設 定對應之 Spar 型浮體平台 y^值分佈，浮體平台於自由液面區域之 y^範圍為 200 y^  500，而自由夜面以下之 y^較低，範圍為 0 y^ 200，因此為了求解 自由液面之流場特性需選用 All y^ Wall treatment 壁面函數方法求解黏性效應。

(33)

18

圖 2-2 計算網格之 y^值

2.1.4 數值方法

本研究使用有限體積法將計算空間進行離散，原理為將空間域劃分為許多為微小的體積網格，網格的幾何中心稱為單元中心，計算時將流場的變數置於每個控制體積的單元中心上。數值求解為在空間離散建立後對統御方程式進行空間離散，進而建立求解的代數聯立方程式，並且搭配邊界條件求解流場中的變數。求解壓力與速度耦合時採用 SIMPLE(Semi-Implicit Method for Pressure-Linked Equations)演算法[22]，為壓力耦合方程半隱式解法，其原理：先假設壓力場代入離散後之動量方程式，求得速度向量，再代入連續方程式中求解壓力修正方程式進而求解修正過後的壓力值，最後再將壓力場代入動量方程式得到新的速度場，

如此進行迭代直到壓力與速度達到收斂。

(34)

19

2.1.5 雙相流數值方法

本研究採用流體體積法[23] (Volume of Fluid Method, 簡稱 VOF)進行空氣與水兩者的流場計算，此法為利用一種等效流體來取代兩種流體，等效流體的物理 性質(如密度與黏性)，為兩種流體的體積分數(Volume Fraction,c)內插求得，

本研究的自由液面位置定義在c0.5處，如圖 2-3 所示。VOF 法提供經由固定網格追蹤明顯流體介面方法，其基礎在於流體體積分數，它包含於每一個控制體積中，控制體積可以為空氣(c0)、水(c1)以及空氣與水混合體(0 c 1)，以圖 2-3 說明，上為空氣，下則為水。

Vw

c V (2-22)

(1 )

w a

ρ cρ  c ρ (2-23)

(1 )

w a

μ cμ  c μ (2-24) 守恆方程式：

Dc 0

Dt  (2-25)

圖 2-3 體積分率說明

(35)

20

2.2 風力發電機理論

本節介紹風力發電機的理論基礎，包括動量理論、葉片元素理論以及葉片元素動量理論，將在以下分別介紹。

2.2.1 動量理論

動量理論為能量守恆觀點，上游與下游的動能差即為整個控制體積內受到的作用力，下面將介紹軸向動量理論與角動量理論。

(1) 軸向動量理論

軸向動量理論為假設有一個空氣的流管，風機對流場的影響集中於此一流管中，流管邊界為一流線，此流線為在風機作用下所形成之自由流線 (free-streamline)。

當流體流經轉子平面(Rotor Plane / Actuator disk)時，會對轉子平面作功，因此動能會降低，基於質量守恆，流進系統之質量會與流出系統之質量相同，進而產生漸擴之流管，如圖 2-4 所示。

軸向動量理論基於下列假設進行推導：

a. 流體為無黏性(μ0)

b. 流體不可壓縮( ρ = constant ) c. 流場為非旋流( U 0)

d. 無限葉片數，將葉片的掃風面積視為圓盤(Actuator disk) e. 圓盤內均勻受力

f. 圓盤前後速度連續U₂ U₃

(36)

21

圖 2-4 軸向動量理論示意圖

軸向動量理論藉由推導出流入流管的平均流速U₁與流出流管的平均流速U₄ 差異，估算風機轉子平面所受到的推力。

圓盤所受之推力T為前後兩端的壓力差P₂與P₃乘以圓盤面積A，因此可得下列式子：



2 3

 

2 3



²

T  P P A P P 4D (2-26) T 為轉子平面所受推力，P₂與P₃分別為轉子平面前、後方壓力，A為轉子平面面積，D為轉子直徑。

由於軸向動量理論的八個假設，因此可將伯努力方程(Bernoulli equation)應用在位置(1)至(2)與位置(3)至(4)的流線上，公式如下：

2 2

1 2 2

1 1

2 ¹ 2

P ρU P  ρU (2-27)

2 2

3 3 4 4

1 1

2 2

P  ρU P  ρU (2-28) 其中ρ為流體密度，並假設控制體積上下游邊界與風機位置距離無窮遠，代表上下游邊界不再受到風機影響，則壓力會等於一大氣壓，公式如下：

1 4 atm

P  P P (2-29)

(37)

22

將式(2-27)與式(2-28)相減後代入圓盤前後速度連續並相等(U₂U₃)之關係可得：



² ²



2 3 4

1

2 ¹

P P  ρ U U (2-30) 從能量守恆定律，控制體積流出與流入的動量差等於控制體積內所受作用力 (推力)可得下列公式：

1 1 4 4

mU m U T (2-31)

因為質量守恆定律mm₁m₂m₄，並且將式(2-26)代入式(2-31)可將式(2-31) 改寫如下：

1 4 ( 2 3)

mU mU  P P A (2-32) 其中，A 代表圓盤截面積， m 為流經圓盤的質量流率，公式如下：

mρAU2 (2-33) 將式(2-30)與式(2-32)聯立可得下列公式：

 

2 1 4

1

U  2 U U (2-34) 此時提出軸向速度誘導因子(Axial Induction Factor)a，定義為流體入流速度與流體在圓盤的軸向速度之降低比率，公式如下：

1 2

1

U U

a U

  (2-35)

將式(2-35)代入式(2-34)可得下游流速U₄與軸向速度誘導因子a的關係：

4 1(1 2 )

U U  a (2-36)

利用上游流速與下游流速之關係式(2-36)，代入式(2-31)可得圓盤所受之總推力，推力公式如下：

2 1

1 4 (1 )

T 2 ρAU a a (2-37)

(38)

23

推力係數C_T為：

2 1

1 2

T

C T

ρAU

 (2-38)

將推力公式(2-37)代入推力係數公式，可得推力係數與軸向速度誘導因子的關係，關係如下：

4 (1 )

CT  a a (2-39) 同理可得圓盤所受的功率，所受功率為單位時間內流入流管與流出流管的動能差：

2 2 3 2

w 1 4 1

1 1

P ( ) 4 (1 a)

2m U U 2ρAU a

    (2-40)

功率係數C_P為：

w 3

P 1 2

P

1

C

ρAU

 (2-41)

功率 P 使用式(2-40)代入，可得：

4 (1 a)2

CP  a  (2-42)

當 ¹,1

a3 時，C_P有極值，分別為最大C_P與最小C_P，最大C_P即為 Betz Limit，

為貝茲(Betz)於 1926 年所提出，其值如下：

,max

16 0.5926

P 27

C   (2-43)

此極限代表風機能夠從空氣動能中擷取的最大比例為 59.26%。

(2) 角動量理論

前述的軸向動量理論中，並不考慮轉子旋轉時所造成的旋轉跡流(Rotational wake)，角動量理論是以軸向動量理論為基礎，進一步加入轉子旋轉時對流體產生

(39)

24

的角速度，並且求出氣流對圓盤造成的旋轉力矩。

當轉子以角速度Ω做旋轉運動時，在轉子平面前方並無角速度，因此只有軸向速度U₂U₁(1a)，U₁為流體流入流管的速度，但是當流體流過轉子平面後其角速度會增加至Ω ω ，其中Ω為轉子角速度，ω為轉子平面前後的變化量。

圖 2-5 旋轉跡流示意圖

根據角動量隨時間的變化為物體所受之力矩Q之關係可得下列公式：

dMa

Q dt (2-44)

其中M 為角動量_a M_a Iω，而 I 為轉動慣量I mr²，因此可將力矩整理成：

2

( ) 2

dIω d mr ω

Q mr ω

dt dt

   (2-45)

其中 r 為半徑。如圖 2-5，葉片單位半徑產生環狀面積的力矩dQ為：

dQdmr ω2 (2-46)

因質量流率dmρU dA₂ ，其中 dA 為環狀面積 2πrdr ，dr 為單位半徑，因此可 得：

22

dm ρU πrdr (2-47) 將式(2-47)代入式(2-46)可整理得單位環狀面積下所受之力矩dQ為：

2

2 2

dQ ρU r πrdr (2-48)

(40)

25

圖 2-6 轉子圓盤環狀面積力矩

此時定義旋向速度誘導因子(angular velocity induction factor)：

' 2Ω

a  ω (2-49)

利用式(2-49)旋向誘導速度與旋向誘導因子之關係、式(2-35)轉子平面軸向速度與軸向誘導因子可將每一環狀面積所得之力矩dQ整理得：

2 1

4 '(1 )1 Ω 2

dQ a a 2ρU r πrdr (2-50) 進一步可求出環狀面積上的功率dP_w為：

2 3

w 1

P Ω 4 (1 ) 'Ω

d dQ  πρU a a r dr (2-51)

2.2.2 葉片元素理論(Blade Element Theory)

葉片元素理論是將風機葉片幾何沿著翼展(Span)方向分成 n 個截面，每個截面的受力是由此截面翼型的升力係數與阻力係數決定，而且假設每個截面不會互相影響，在這些前提下可求出每個截面平行於轉子平面的切向力以及垂直於轉子平面的正向力，並且可藉由切向力與截面幾何所在半徑位置的乘積求得力矩進而得到功率。示意圖如圖 2-7，風機葉片半徑為R，截面所在半徑位置為 r ，截面弦長 為c，截面長為 dr 。

(41)

26

R r

c

dr

圖 2-7 葉片元素理論示意圖

每個葉片截面都有相對風速U_rel，此相對風速是風速、轉子轉速之切線速度、

軸向誘導速度、切線誘導速度之向量所組成，其中風速為U₁，軸向誘導速度為U a₁ ，轉子轉速之切線速度與流體之切線誘導速度兩者的向量合為：

Ω ( ) Ω Ω ' Ω (1 ')

2

r ω r r ra  r a (2-52)

所以轉子的總切線速度為Ω (1r a')，總軸向速度為U₁(1a)。因此葉片截面之速度、

力分量如圖 2-8，其中β為截面的總節距角(section pitch angle)，為整支葉片的節距角θ 與截面的扭轉角_p_,0 θ_T之總和。φ為截面之相對風速與轉子平面之夾角，α為攻角。

圖 2-8 二維翼型速力圖

(42)

27

由圖 2-8 可知相對入流角φ為：

1(1 ) 1

tan Ω(1 ) 1 ' ^r

U a a

φ λ

r a' a

 

 

  (2-53)

其中λ_r的定義為：

1

Ω

r

λ U

 r (2-54)

由於可得相對入流角φ，因此相對速度U_rel以及斷面對應之攻角α為：

1(1 ) sin( )

rel

U a

U φ

  (2-55)

α φ β  (2-56)

已知攻角α後，可得此斷面之升力係數C_L與阻力係數C_D，再由升阻力係數之 定義可算得此斷面於截面長 dr 下之升力 dL 、阻力 dD 為：

1 2

2 ^rel ^L

dL ρU cC dr (2-57)

1 2

2 ^rel ^D

dD ρU cC dr (2-58)

其中c為弦長。利用截面升阻力與相對入流角φ之關係可得截面於垂直轉子平面之正向力dF_N為：

cos( ) sin( )

dFN dL φ dD φ (2-59) 也可得平行轉子平面之切線方向力dF_T 為：

sin( ) cos( )

dFT dL φ dD φ (2-60) 由於需考慮風機葉片總數為B，並將式(2-57)與式(2-58)代入式(2-59)可整理得此葉片截面半徑上的總推力：

 

2

D

1 cos( ) C sin( )

N 2 rel L

dF Bc ρU C φ  φ dr (2-61) 同理，葉片半徑上總切線方向力為：

(43)

28

 

1 2

sin( ) C cos( )

T 2 rel L D

dF Bc ρU C φ  φ dr (2-62) 葉片上之力矩dQ為切線方向力dF_T 與截面半徑位置 r 之乘積：

dQdF rT (2-63) 將式(2-62)代入式(2-63)可將斷面所受力矩整理成：

 

1 2

sin( ) C cos( )

2 ^rel ^L ^D

dQBc ρU C φ  φ rdr (2-64)

2.2.3 葉片元素動量理論(Blade Element Momentum Theory)

在 2.2.1 與 2.2.2 項分別介紹了動量理論與葉片元素理論，將這兩種理論結合即為葉片元素動量理論(BEM)。BEM 之目的為求出斷面幾何對流體產生的軸向速度誘導因子(Axial induction factor)以及旋向速度誘導因子(Angular induction factor) 也就是a與 'a 的關係式。另外，從式(2-53)可知a、 'a 與φ三者是耦合的，所以在 BEM 中必須利用迭代法找出a、 'a 與φ。

在此將前兩小節的公式整理如下：

整理軸向動量理論推導出推力的公式(2-37)：

2

4 (1 ) 1

dT  a a ρU πrdr (2-65) 整理角動量理論推導出力矩的公式(2-50)：

2

4 '(1 ) 1Ω

dQ a a ρU r πrdr (2-66) 整理葉片元素理論剖面正向力公式(2-62)以及力矩公式(2-64)：

 

2 2

1 2

(1 )

' cos( ) C sin( )

N sin L D

U a

dF σ πρ C φ φ dr

φ

   (2-67)

 

2 2

1 2

(1 )

' sin( ) C cos( )

sin ^L ^D

U a

dQ σ πρ C φ φ rdr

φ

   (2-68)

其中 'σ 為局部弦周比(Local solidity)

(44)

29 ' 2 σ Bc

 π (2-69)

為了求出a與 'a 的關係式，令軸向動量理論之推力與葉片元素理論之正向力相 等，即式(2-65)與式(2-67)相等得到：

 

2

cos( ) sin( ) 1 ' 4sin

L D

C φ C φ

a σ

a φ

 

 (2-70)

接下來令角動量理論之扭力與葉片元素理論之力矩相等，即式(2-66)與式(2-68) 相等可得：

 

^D

2

sin( ) C cos( )

' '

1 4 sin

L r

C φ φ

a σ

a λ φ

 

 (2-71)

整理式(2-70)可得軸向速度誘導因子a成：

 

2

1 1 4sin

' _Lcos( ) _Dsin( )

a φ

σ C φ C φ



 

(2-72)

將軸向速度誘導因子a之公式(2-72)代入式(2-71)即可整理旋向速度誘導因子 '

a 為：

   

 

' 1

2sin cos

' _Lcos( ) _Dsin( ) 1

a φ φ

σ C φ C φ



 

(2-73)

在葉片元素動量理論中為了找出φ角，因此必須給a與 'a 的起始值，再利用迭 代方法找出最適當的a、 'a 與φ。當迭代已達到收斂時，可經由式(2-67)、式(2-68) 得到斷面的推力與力矩，每個葉片斷面均迭代完成後可經由積分方式計算整體風機所承受之總推力T與功率P_w。

h

R

T 



r dT ^(2-74) Pw Ω

h

R

r dQ





^(2-75)

其中，r_h為輪轂半徑，積分完成後即可評估葉片效能。

耦合BEM與CFD方法計算浮體式風機於規則波中運動之研究

國立臺灣大學工學院工程科學及海洋工程學系 碩士論文

National Taiwan University Master Thesis

耦合 BEM 與 CFD 方法計算浮體式風機於規則波中運動之研究

A Coupled Blade Element Momentum – Computational Fluid Dynamics Model for Evaluating Floating Wind Turbine

Motion under Linear Waves

李仲凱 Chung-Kai Lee

指導教授：郭真祥 博士

Advisor: Jen-Shiang Kouh, Dr. Ing.

中華民國 103 年 7 月

July, 2014

誌謝

摘要

Abstract

目錄

圖目錄

表目錄

符號表

第1章 緒論

1.1 研究背景

1.2 離岸型風力發電介紹與浮體式風機介紹

1.3 文獻回顧

1.4 研究目的與方法

1.5 本文架構

第2章 理論基礎

2.1 計算流體力學





2.2 風力發電機理論



 







 

 

 

 

 

 

 

 

 

   

 





國立臺灣大學工學院工程科學及海洋工程學系碩士論文

指導教授：郭真祥博士

第1章緒論

第2章理論基礎