結論

12 雷諾平均化納維爾-史托克方程式(Reynolds Averaged Navier-Stokes 簡稱 RANS)。

13

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

2.1.2 紊流模型

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

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

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

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

14

15

0.5, 0.85, 0.075 0.553

1.0, 0.856, 0.0828 0.440

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^的區域稱為黏性次層流，關係如下：

17

18

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

2.1.4 數值方法

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

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

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 體積分率說明

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₃

21

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

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

23

24

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 。

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 二維翼型速力圖

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

28

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 與φ三者是耦合的，所以在

29

30 限葉片的功率會高出許多，因此使用以升力面方法(Lifting Surface Method)推導出 葉尖損失的數學模型[16]，此模型產生之修正項F介於 0 到 1 之間。

31

(2) Glauert Correction

Glauert Correction 是一個經驗公式用來針對不合理的軸向速度誘導因子a做 修正，因為在軸向動量理論中推導出式(2-36)，當a0.5 時，U₄變成負值代表產生 回流，但實際上並不會發生回流現象，因此當a0.5時，氣流流過轉子後，實際 流場會相較於軸向動量理論中預測的流場還要複雜，所以上一小節求得的a就要 進行修正。

Glauert Correction [24]經驗公式為：

 

32 一個六自由度的剛體運動(6-DOF Rigid Body)，此六自由度為三軸的平移運動和三 軸旋轉運動，座標系統將中心點設於浮體式風機之中心，依照右手定則定義入流 風方向為 X’軸正向，Z’軸方向與重力方向相反，其中沿 X’-Y’-Z’三軸方向之平移 運動依序為：縱移(Surge)、橫移(Sway)及起伏(Heave)，另繞 X’-Y’-Z’三軸之旋轉 運動分別為：橫搖(Roll)、縱搖(Pitch)及平擺(Yaw)。如圖 2-9 所示。

33

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

本研究模擬浮體式風機考慮了縱搖，故套用使用之軟體內部提供的六自由度 模型求解器(Six Degree of Freedom solver, 6-DOF solver)。此求解器可模擬剛體運動，

並藉由計算物體在流場中受到水動力、重力或使用者設定之外力(如 BEM 法計算 之氣動力及繫纜力)等，所構成之合力與合力矩，再進一步求解剛體運動方程式，

以計算物體下個時間所在的位置。剛體運動相關方程式如下所示[25]：

剛體平移方程式：

md dtV 

f (2-90)

剛體轉動方程式：

M d M

dtω  

ω ω n (2-91)

其中，m為剛體質量，f 為作用在剛體上的合力向量， V 為剛體質心之平移 速度向量，M 為慣性矩，ω為剛體質心之轉動速度向量， n 為剛體所受之合力矩 向量。

34

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

圖 2-10 為浮體式風機所用之參考座標系，其中使用二個卡式座標系統，說明 如下：

(1) XYZ 座標(紫色)：為固定之世界座標系統，原點座標為(0, 0, 0)，依右手定則 定義 X 軸向與入流方向相反，Z 軸向與重力方向相反。此座標系用來描述流 體之流動速度向量。

(2) X’Y’Z’座標(綠色)：為浮體式風機之局部座標系統，座標原點設置於浮體式風 機整體之重心上，用來描述浮體式風機之六自由度運動，此座標系相對於世 界座標系會隨著風機不同的姿態而改變，但對於計算域是固定的。

選用軟體內之六自由度模型求解器的方式為旋轉或平移整個計算域來模擬剛 體運動，因此在內部不需要重新生成網格，所以可省略重新長網格之計算時間，

流場計算域的運動是參考剛體的局部座標系統，而流場中的流體則持續以世界座 標系為參考座標保持流動，因此當計算域在運動時空間中同一位置的流體之物理 性質(如：速度、壓力)是保持不變的。

35

30，並據平面前進波理論在無限水深(Deep water)之條件下可推導 出四個重要關係式，公式如下[26]：

本研究以美國國家再生能源實驗室(National Renewable Energy Laboratory)單 機容量為 5MW 搭載 Spar 型浮體平台與半潛型(Semisubmersible)浮體平台為研究對 象，此兩種浮體式平台均搭配懸垂式錨鍊。錨鍊系統以位移-力的方式計算，即為 利用錨鍊的繫錨點(錨鍊與浮體式風機的接點)位移來計算繫錨力之變化量以及繫

36

錨力大小，而浮體式風機運動後受錨鍊之作用力是以 JAVA 副程式進行計算。計 算繫錨力之步驟如下：

(1) 計算浮體式風機初始狀態時懸鍊線於浮體式風機繫錨點上之作用力(簡稱初始 繫錨力)，包含繫錨點的水平力T_H⁰與垂向力T_V⁰。

(2) 當浮體式風機運動時，利用繫錨點的位移計算繫錨力的變化量dT 與_Hⁱ dT 。 _Vⁱ (3) 由時間步 i 時的繫錨力T 、_Hⁱ T 加上繫錨點位移過後求出之_Vⁱ dT 、_Hⁱ dT ，即可_Vⁱ

求得T_H^i1、T_V^i1。

(4) 將求得之繫錨力代入浮體式風機之剛體運動方程式，即可計算下一時刻之運 動狀態。

(5) 重複進行(2)~(4)之過程，直到運動達成收斂。

繫錨點與繫錨力相關資訊如下圖所示：

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

其中 B 為懸鍊線固定於海床處，P 為繫錨點，Lc 為懸鍊線長度(m)，H為懸鍊

其中 B 為懸鍊線固定於海床處，P 為繫錨點，Lc 為懸鍊線長度(m)，H為懸鍊

在文檔中 耦合BEM與CFD方法計算浮體式風機於規則波中運動之研究 (頁 26-123)