複合材料葉片振動行為之研究

國 立 交 通 大 學

工學院專班精密與自動化工程學程

碩 士 論 文

複合材料葉片振動行為之研究

A study on vibration behavior of composite sandwich wind blade

研 究 生 : 唐 榕 崧

指導教授 : 金 大 仁 教授

複合材料葉片振動行為之研究

A study on vibration behavior of composite sandwich wind blade

研 究 生: 唐 榕 崧

Student：Jung-Sung Tang

指導教授: 金 大 仁 博士 Advisor：Tai-Yan Kam

國 立 交 通 大 學

工學院專班精密與自動化工程學程

碩 士 論 文

A Thesis

Submitted to Degree Program of Automation and Precision Engineering

College of Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of

Master of Science

in

Automation and Precision Engineering

February 2009

Hsinchu, Taiwan, Republic of China

複合材料葉片振動行為之研究

研究生：唐榕崧 指導教授：金大仁 博士 國立交通大學工學院專班精密與自動化工程學程 摘要 本文重點分兩部份，第一部份介紹 Wilson 所推導之葉片元素動量理 論，引用其理論說明葉片幾合外型計算方式及演算流程，並計算風力作用 於葉片上的氣動載荷，將瞬時氣動載荷施於葉片上分析葉片振動行為。 第二部份介紹複合材料三明治結構風力葉片自然頻率的實驗量測方 法、葉片受瞬時激振力下的振動試驗方法，以及提供簡單的阻尼量測實驗 及估算方式，並利用有限元素分析軟體 ANSYS 進行分析；本文先以鋁板比 較振動實驗數據與有限元素分析結果，驗證有限元素分析可行性，同法分 析葉片振動行為，並以實驗驗證之，最後分析葉片受瞬時陣風作用下的應 力分佈情形，本研究結果將有助於葉片破壞、失效分析、及可靠度評估。

A study on vibration behavior of composite sandwich wind

blade

Student : Jung-Sung Tang Advisor : Dr. Tai-Yan Kam

Degree Program of Automation and Precision Engineering National Chiao Tung University

ABSTRACT

This work consists of two parts. The first part reviews the methods for determining the wind loads on wind blades, in particular, Wilson's blade element momentum theory. Based on Wilson's method, an algorithm is presented for determining a wind blade's geometry and aerodynamic load. The distribution of instantaneously applied aerodynamic load on the wind blade is analyzed and used for the wind blade's vibration analysis.

Part two deals with the experimental and theoretical studies of a composite sandwich type wind blade. The wind blade's natural frequencies and its response under impulsive loads are measured. Damping of the wind blade is

approximated using the measured response. A finite element model of the wind blade is constructed using the FEM package ANSYS. Validity of the FEM code for vibration analysis is first demonstrated by comparing the theoretical results with the experimental data for an aluminum plate. The finite element analysis of the composite wind blade is then analyzed. The experimental results have also validated the accuracy of the finite element model. The composite wind blade subjected to instantaneously applied wind load is analyzed to study the stress distribution in the blade. The results obtained in the study will be useful for the failure analysis of the blade.

誌謝 研究所求學過程中，首先要感謝指導教授 金大仁博士對我在研究方法 上的指導與輔正，感謝中山科學研究院各級長官，在我進修這一段時間諸 多關懷與問候，特別是組裡的同仁協助學生解決理論基礎上的各方問題， 使得對於不熟悉領域的觀念有深入的認識與了解，從而解決研究上的許多 問題。 求學過程中感謝研究室學長在實驗上協助，完成了架設振動量測實驗 平台，並建立了實驗數據分析軟體架構，感謝專班的同學在求學過程中一 起成長，互相學習，並要感謝交大游泳隊巫教練以及全體泳隊同學，在求 學過程中增加了歡樂氣份，增添了豐富的色彩，感謝體育室的王教練以及 健身房的學長和有氧課的同學，強健了我的體魄，並維持健康的身體。 榕崧 2009.02 于桃園

目錄 中文提要 ……… i 英文提要 ……… ii 誌謝 ……… iii 目錄 ……… iv 表目錄 ……… vi 圖目錄 ……… vii 一、 緒論……… 01 1.1 前言……… 01 1.2 文獻回顧……… 02 1.3 研究方法……… 04 二、 基本概念及原理介紹……… 06 2.1 風力機種類……… 06 2.2 葉片空氣動力學理論……… 07 2.2.1 渦流效應……… 08 2.2.2 葉尖損失效應……… 09 2.2.3 二維翼型理論……… 09 2.2.4 葉片元素動量理論……… 12 2.2.5 貝茲（Batz）理論……… 13 2.3 葉片幾何外型介紹……… 14 2.3.1 翼剖面/翼型（airfoil）外型的選擇……… 14 2.3.2 旋轉最大外徑……… 17 2.3.3 選擇翼尖流速比率……… 17 2.3.4 葉片數量……… 18 2.3.5 葉片寬度……… 19 2.3.6 葉片扭轉角……… 19 2.3.7 徑向位置 r 處效率因子最佳化……… 21 三、 葉片幾何模型計算及建立……… 22 3.1 葉片幾何模型計算流程………. 22 3.2 葉片氣動性能計算流程……… 24 3.3 實例計算……… 26 3.4 三維立體模組建立……… 30 四、 有限元素分析……… 32 4.1 自然頻率、自然模態分析……… 32 4.1.1 ANSYS 模型建立步驟……… 32 4.2 暫態振動分析……… 34

4.2.1 ANSYS 模型建立步驟……… 34 4.2.2 ANSYS 模擬分析中各參數的取得……… 36 4.3 有限元素分析驗證……… 37 五、 葉片振動試驗……… 39 5.1 自然頻率量測實驗……… 39 5.2 阻尼量測實驗……… 40 5.3 應變量測實驗……… 41 5.4 有限元素分析與實驗結果討論……… 42 六、 結論與未來研究方向……… 47 6.1 結論……… 47 6.2 未來研究方向……… 48 參考文獻 ……… 50 附錄一 效率因子最佳化程式碼……… 54 附錄二 葉片幾何參數限性化程式碼……… 57 附錄三 葉片三維立體模組自動化繪製程式碼……… 58

表 目 錄 表 1 各種能源的發電成本……… 60 表 2 旋轉最大外徑尺寸對應輸出功率……… 60 表 3 各站位參數計算結果……… 61 表 4 各站位參數線性化計算結果……… 62 表 5 材料常數設定……… 62 表 6 鋁板在無拘束條件下的自然頻率實驗與分析結果比較…… 63 表 7 鋁板在單邊拘束條件下的自然頻率實驗與分析結果比較… 63 表 8 加速度與應變之實驗與分析結果比較……… 64 表 9 自然頻率實驗與分析結果比較……… 64 表 10 葉片單邊夾持下加速度之實驗與分析結果比較……… 65 表 11 P01 與 P02 應變結果之實驗與分析比較……… 65

圖 目 錄 圖 1 荷蘭式風車………. 66 圖 2 多葉片式風車………. 66 圖 3 徑流吸力式………. 67 圖 4 螺旋槳式風車……… 67 圖 5 達利亞斯式（Darrieus）風車……….. 68 圖 6 S 型式（Savonius）風車………. 68 圖 7 槳式（Paddle）風車………. ……… 69 圖 8 不同風力機型式 Cp對 λ 曲線……… 69 圖 9 環境風對風力機發生的相互關係……… 70 圖 10 葉片幾何……… 70 圖 11 氣流速度三角圖……… 71 圖 12 葉片受風作用力分解圖……… 71 圖 13 葉片受風作用力簡圖……… 72 圖 14 相對風造成葉片壓力分佈示意圖……… 72 圖 15 葉片承受作用力分佈示意圖……… 73 圖 16 Batz 理論模型……… 73 圖 17 翼剖面的幾何定義……… 74 圖 18 CP與λ的關係……… 74 圖 19 λ與CP、B的關係……… 75 圖 20 α -CL、α -CD的曲線……… 75 圖 21 葉片模型計算流程圖……… 76 圖 22 翼剖面使用站位說明……… 77 圖 23 各型號之翼型曲線……… 77 圖 24 分段的情況說明……… 77 圖 25 尋找最佳攻角示意圖……… 78 圖 26 C_r對站位及β_r對站位曲線……… 78 圖 27 葉片俯視圖……… 79 圖 28 各元素的坐標系統示意圖……… 79 圖 29 力對時間圖的資料點……… 79 圖 30 阻尼量測試驗之加速度及位移對時間曲線……… 80 圖 31 鋁板無拘束條件下自然頻率量測試驗……… 80 圖 32 鋁板在無拘束條件下的自然頻率量測結果……… 81 圖 33 鋁板單邊拘束條件下自然頻率量測試驗……… 81 圖 34 鋁板在單邊拘束條件下的自然頻率量測結果……… 82 圖 35 加速度及位移對時間變化曲線……… 82 圖 36 實驗與分析之加速度對時間變化曲線……… 83

圖 37 實驗與分析之應變對時間變化曲線……… 83 圖 38 實驗與分析之加速度趨勢線比較……… 84 圖 39 實驗與分析之應變趨勢線比較……… 84 圖 40 複合材料葉片示意圖……… 85 圖 41 自然頻率量測試驗台架……… 85 圖 42 自然頻率量測結果……… 86 圖 43 阻尼量測實驗示意圖……… 86 圖 44 應變規黏貼示意圖……… 87 圖 45 有限元素模型……… 87 圖 46 各階自然模態……… 88 圖 47 自然頻率實驗與分析結果比較曲線……… 89 圖 48 邊界挾持下葉片的自然頻率……… 89 圖 49 實驗與 ANSYS 分析結果比較曲線……… 90 圖 50 應變實驗量測結果……… 90 圖 51 網格細緻處理前、後比對圖……… 91 圖 52 應變量測點 P01 實驗與分析比較結果……… 91 圖 53 應變量測點 P02 實驗與分析比較結果……… 92 圖 54 應變量測點 P01 實驗與分析之應變趨勢線比較……… 92 圖 55 應變量測點 P02 實驗與分析之應變趨勢線比較……… 93 圖 56 風速及葉片受力與時間關係圖……… 93 圖 57 不同時間下葉片受力對徑向位置之變化曲線……… 94 圖 58 瞬時陣風作用下葉片徑向方向應變結果……… 94 圖 59 瞬時陣風作用下葉片 Von Mises 應變結果……… 95 圖 60 瞬時陣風作用下葉片位移及加速度對時間變化曲線……… 95 圖 61 應力分析網格細致處理示意圖……… 96 圖 62 12ms 時的葉片迎風面應力分佈……… 96 圖 63 12ms 時的葉片背風面應力分佈……… 97 圖 64 12ms 時整體葉片變形情形……… 97

一、 緒論 1.1 前言 根據電力公司資料顯示，台灣每年電力需求逐漸上昇，科技的進步帶 來能量需求的大幅提升，如何提升足夠的能量供應，而又必須符合環保的 需要（二氧化碳與溫室效應問題），兼顧經濟成長、環境保護的永續發展， 這已不僅僅只是台灣單一地區的問題，我國目前能源供應有97％來自石化 燃料；全球暖化速度加快、環保意識逐漸受到重視，國際間已在1997 年簽 訂京都協議書，提出對於各國降低二氧化碳的排放量要求，並在2005 年二 月十六日正式生效，不僅僅於此，國內的能源使用一直以來面臨著許多挑 戰，擴大核能的使用決策，將面對核廢料處理的問題，核電廠除役後，如 何尋求替代能源更是一大難題，對於集中式電力供應系統，在容易發生地 震、颱風等自然災害的國內而言，大大降低了區域電力供應的自主性與彈 性；石化燃料不僅是環保問題，國際能源的價格波動，以及能源的來源是 否確保，皆深深影響供電的穩定性，全球經過了1974 年及 1980 年兩次石 油危機，已帶來深深的影響及警惕，於2004 年底國際油價已衝到每桶原油 50 美元，2006 年更是飆到每桶原油 78.65 美元，2007 年底更是逼進每桶原 油100 美元，其他石化燃料也是同步上漲，有限的資源逐漸減少價格日漸 高漲是必然的趨勢，可見在生能源的重要性、急迫性[1]~[10]。 在生能源的利用已在國際間受到相當大重視，風力發電亦屬於再生能

源種類之ㄧ，根據資料顯示，風力發電為全球增長最快的能源，每年增長 超過30％。全球風力發電在 2003 年約達到 3200 萬千瓦，大約相等於 32 座 標準的核電廠，近幾年來全球風電累計增長率，一直保持在30％以上，顯 示其在市場上的成長及為快速，技術的成熟已經促使許多國家積極開發益 於環境的風能，可見風力發電是決不可忽視的。 1.2 文獻回顧 關於風力機葉片氣動理論方面，Rankine 在 1865 年[11]用線性動量理論 形成了螺旋槳流場的簡單模型。德國空氣動力科學家Betz[12]於 1919 年推 導出風力機風能利用系數所能達到的最大值，然而實際狀況下無法達到該 理想情形，Glauert[13]考慮了 Schmitz 提出的渦流效應影響，把動量理論結 合葉片元素理論來分析螺旋槳、風車周圍的流動並將Rankine 一維流發展成 為有旋轉作用的二維流。1974 年時 Wilson 和 Lissaman[14]將葉片元素動量 理論應用於風力機，考慮葉尖損失效應並引用Prandtl 推導的葉尖損失因 子，形成了經典的葉片元素動量理論理論（Blade Element Momentum， BEM），此理論簡單且易於使用。從 Wilson 和 Lissaman 之后，葉片元素動 量理論理論被廣泛用于風力機的設計[15]和性能計算[16]，並且用來確定風 力機的動態載荷[17]，同時葉片元素動量理論還不斷地被進一步改進和完善 [18]。

風力機葉片之自然頻率、模態及振動行為在數值計算方法研究的演化 上，傳統採用旋翼機葉片經典理論，自1958 年 Houbolt 和 Brooks 推導出了 擺振-扭轉偶合的運動微分方程[19]，1978 年 Wendell 考慮了氣動載荷提出 風力機葉片氣動彈性穩定性問題[20]，同時期 Kottapalli 等人用非偶合的旋 轉效應模態對該問題進行研究，得出了葉片動態響應[21]，Chopre[22][23]、 Miller 和 Dugundji[24]等開始研究氣動彈性響應問題，近年來 Chaviaropoulos 以簡單的彈簧支承的剛性葉片模型研究葉片的擺振響應[25]，Thomsen 和 Petersen 用葉片的中心鉸接模型研究葉片的擺振穩定性[26]，Hansen 研究單 獨葉片的擺振、扭轉的耦合運動[27]，由於風力機葉片的幾何構型相當複 雜，近年來葉片所使用的材料也由金屬轉變為複合材料，又普遍以三明治 結構為主，1987 年時 David 用有限元素方法分析葉片動態問題[28]，2000 年時許多中國科學家[29][30]開始投入有限元素分析應用於風力機葉片上。 在研究三明治板的文獻方面，Mau [31]提出多層一階平板理論及其他學 者陸續提出的相似理論[32~35] ，Kant[36、37]基於高階位移場模式發展出 的等參單元有限元素法來做三明治結構的彎曲分析；對於其動態特性也就 是自然頻率(Natural Frequencies)及其相應的模態(Mode Shape)的分析方 式，在許多有限元素法或結構動態分析的文獻[38-47] 裡諸多介紹，1971 年 Bathe K.J.提出的子空間迭代法(Subspace Iteration Method)[48]。

1.3 研究方法 根據工研院2006 年的資料，各種能源的發電成本如表一[6]，傳統能源 會因為油價的高漲反應在能源成本上，在生能源的必要性由此可見，台灣 目前風力發電機的零組件全部都來自於國外進口，其中葉片屬於消耗性零 件，長期不斷的運作下，可能遭受昆蟲撞擊、電擊、疲勞、火災等因素而 損壞，如國人有能力自行自製、生產，則可確保風力機持續運轉，減少進 口等待的時間，開發自製及維修能力，除了能大幅降低成本，更能為國內 帶來大量工作機會，目前世界各大風力機廠家對於風力機零組件中，自製 的項目裡一定有控制系統以及葉片模組，可見葉片對於風力發電產業的重 要性，國內可以先從中小型葉片開始發展，培養能夠自製及維修的能力， 並通過國際認證後，可外銷其他發展中國家，未來可持續研發大型葉片， 以提供各種不同的需求。隨著風力發電市場的快速發展，葉片的需求急速 增長，除了新增設的風力機外，汰舊換新以及損壞維修等，都是良好的機 會，可見葉片的商機以及發展潛力相當龐大。 第四章將介紹有限元素分析軟體ANSYS 的使用及相關設定，並使用鋁 板進行動態試驗與ANSYS 分析結果比較，驗證 ANSYS 分析的可行性，第 五章分別介紹葉片於無拘束情形下自然振動情形，提供自然頻率量測的實 驗方法，比較分析與實驗結果以確定模組的正確性，從結果中了解前五階

的自然頻率範圍以及振動模態行為，提供設計時避免與風力機葉片轉速達 到共振；接著介紹葉片單邊受拘束條件下，受一瞬時激振力的暫態分析， 了解葉片受瞬時激振力下所產生的變形行為、及應力、應變分佈情形，並 說明葉片暫態試驗方法，將實際試驗量測結果與分析數據做比較，找出葉 片的阻尼比，了解其振動行為。最後引用葉片元素動量理論計算氣動載荷， 將氣動載荷對時間的變化施加於葉片進行分析，了解葉片受瞬時陣風作用 下的振動行為。 本文介紹之風力機葉片採用複合材料三明治結構，此材料是由高強度 非等向性材料的面層(face sheet)與等向性低密度材料的夾心層(core shee)所 構成其材料性質在層板與層板間是不連續的，因而使得層板在受力狀況之 下，應力的分佈和變形會有很大的差異，而每一層間的情形也相當不同。 幾何外型以機械繪圖軟體SolidWorks 繪製，中繼檔案交換格式利用 Parasolid 作為與ANSYS 軟體交換介面，有限元素模型的網格在芯材採用高階體元素 (solid 186)，面層採用高階殼元素(shell 99)，分析自然模態（Modal）時，計 算特徵植的演算法用子空間迭代法（Subspace），分析暫態（transient）的 演算法使用全矩陣法（Full）。

二、 基本概念及原理介紹 2.1 風力機種類 風力機的分類有許多種，通常可以用容量大小、轉子配置方式、力學 型態、旋轉速度等等做為分類的方式，本文採用轉子配置方式型態的不同 做一個簡單的介紹，轉子不同的配置方式造成氣流的前進方向有所不同， 可分成(1)水平軸式、以及(2)垂直軸式[49]~[59]。 (1)水平軸式風力機：其主要的特徵就是葉片轉子圍繞著一個水平軸旋轉， 且氣流延著水平軸流動，典型的代表有荷蘭式風車（圖1）、多葉片式風 車（圖2）、徑流吸力式（圖 3）、螺旋槳式等等；目前較有純熟技術的 則為螺旋槳式，螺旋槳式風車（圖4）也已應用於大型風力發電機上，目 前有良好的成效，而在螺旋槳式中又可分做上風型和下風型，上風型葉片 位於塔的前方，必須裝置方向控制功能，通常採用設計垂直安定面控制方 向，下風型葉片位於塔的後方，可免裝方向控制功能，但容易造成低頻噪 音，對人體產生有感壓力，因而不常見使用。 (2)垂直軸式風力機：主要特徵為風機轉軸與風向垂直的風力機稱之，其主 要優點為可以接受任何方向來的風，無需控制方向的功能機構，且整體結 構簡單，維修容易，因為此型式之風力機的發電機、齒輪箱可裝設於地面， 易於維護，常見的形式有達利亞斯（Darrieus 圖 5）、S 型（Savonius 圖 6）、

以及槳式（Paddle 圖 7）等三大類，也有將其中二類混用的混合型；此種 型態最常使用的為達利亞斯與S 型混合型態，目前也是各國努力研發的主 力。 以上分類除了S 型在氣動力學屬性上為阻力型（抗力型），其餘都屬於升力 型（楊力型），之間的優缺點除了以上探討的機構複雜度、附屬機電設備設 計難易性、噪音問題、維修問題等等，最重要的取決項目就是做功的效率， 從圖8 可以得知，上風型的三葉片設計具有較高的效率。 2.2 葉片空氣動力學理論 上風型螺旋槳式風力機葉片設計通常採用Glauert法和Wilson法，Wilson 法主要考慮了動量理論、渦流效應、及葉尖損失效應，在以上的理論基礎 上並假定葉片各個徑向斷面之間互相獨立、忽略有限葉片數對氣流的週期 性影響、忽略翼型阻力[11~18]；本文將分別以實際狀態及理想狀態討論葉 片設計，實際狀態所指的就是基於Wilson理論所推導出的，而理想狀態下指 的就是基於Betz理論下並忽略渦流效應、及葉尖損失效應等；關於動量理 論、渦流效應、葉尖損失效應、及Betz理論如后說明[60~62]。 首先介紹其相關參數定義及原理，簡單的說就是葉片受風的作用力而 產生動作旋轉，因旋轉而產生軸功，再將軸功利用電磁產生電能，環境產 生的風對風力機產生的相互關係，以圖9的簡圖來說明。葉片的幾何定義如

圖10所示，D為風力機旋轉最大外徑（單位m）、R為葉片的最大半徑長度

（單位m）、r為葉片某個徑向長度函數（單位m）、dr為葉片某段徑向半徑 （單位m）、C為翼弦長（Chord、單位m）、B為葉片數量、φ為傾斜角（單 位°）、α 為葉形攻角（angle of attack、單位°）、β為葉片扭轉角（blade twist angle、單位°）。 2.2.1 渦流效應 渦流效應主要是假設葉片為一相同半徑之圓盤，不考慮黏性效應，當 均勻環境風速V₁（單位m/s）通過一掃掠面積時，其速度值會減小，在圓盤 上（葉片處）的軸向風速V（單位 m/s）必定小於V₁，通過圓盤後的環境風 速V₂（單位m/s）亦會降低，引入一個軸向干擾因子a值，以計算式表示如 下： 1 a)V -(1 V= (2-1) 1 2 (1-2a)V V = (2-2) Schmitz 提出當考慮有限長的葉片在旋轉時，其翼尖部分會產生翼尖渦 流（wine tip vortex），在與轂（hub）搭接處也會產生渦流區域，因為渦流 的存在造成葉片的能量損耗，因此實際上通過葉尖的流線為一個螺旋線， 葉片的上游及下游所受到的氣流速度三角圖可以圖11 表示之。

rev/s 或 rpm），Schmitz 理論引入一個徑向的誘導因子b值，因此可以得知在 葉片處的徑向風速及氣流對葉片的旋轉角速度可以下式表示之： r b) (1 r 2 U_r ⎟ = + Ω ⎠ ⎞ ⎜ ⎝ ⎛ +Ω = ω (2-3) Ω + = + Ω (1 b) 2 ω _(2-4) 2.2.2 葉尖損失效應 在葉片元素動量理論中假設了葉片為無窮多片，然而這與風力機的葉 片設計有很大的差異，葉片產生的翼尖渦流大大的影響整體葉片效率， Prandtl 針對此一情形提出修正方式，推導出葉尖損失因子（符號F）： ) arcos(e 2 F -f π = (2-5) φ Rsin r -R 2 B f = (2-6) 2.2.3 二維翼型理論 根據二維翼型理論整個葉片受風作用力的情形如圖12 說明相互關係， 以力圖方式簡化表示如圖13，其中U為徑向風速（單位m/s）、W為相對風 速（單位m/s），之間的關係式為： ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = Ω + = R r V ) b (1 r ) b (1 U_{r r r} λ ₁ (2-7) φ sin V W= (2-8) 下標r表示在半徑r 處的特徵值，b為渦流效應產生的徑向的誘導 因子，理想狀態下可假設無此影響，而λ則為翼尖流速比率（Tip speed

ratio），物理意義為最大半徑 R 處的徑向速度與環境風速 V1的比值。 環境的風通過葉片時對於葉片有軸向風及徑向風兩種方向，將兩個方 向的風合成相對風探討葉片受力情形，相對風造成葉片不同作標位置上的 壓力分佈，因壓力差在葉片上氣動中心上產生升力以及阻力，在討論葉片 結構時，利用升力及阻力的合向量淨力，經由作標轉換分解成結構承受的 作用力（Thrust、符號FT、單位N）、以及產生旋轉的作用力（DrivingForce、 符號FD、單位N），實際上這兩種力事分佈在葉片表面上的，如示意圖 15， 為了便於計算，假設二者的力作用於氣動中心並以下列公式求得： r r r r

(r)(Thrust) L cosφ D sinφ

FT = + (2-9)

r r r r

(r)(DrivingForce) L sinφ D cosφ

FD = − (2-10) 由數學式可以看出，影響最主要的參數為升力（Lift、符號L、單位 N）、 阻力（Drag、符號D、單位N），相關定義及計算式如下： dr C 2 1 r 2 r r C W L = _L ρ (2-11) dr C 2 1 r 2 r r C W D = _D ρ (2-12) 其中CL為升力係數（Lift coefficient）、CD為阻力係數（Drag coefficient），

L C 及C_D可經由翼剖面（Airfoil）外型的參考資料中查圖或查表找到。 有了各向分力，可大致估算出風力機的性能軸向推力（符號T、單位 N）、扭矩（Shaft torque、符號M、單位N-m）：

∫

=BRF_T(r) T (2-13)

∫

= R 0 D(r) rF B M (2-14) 為了方便探討其性能，定義無因次的系數C_T推力係數、C_M扭矩係數： 2 2 1 r T(r) _{0.5 dr} T π ρV C = (2-15) 3 3 1 r M(r) _{0.5 dr} M π ρV C = (2-16) 風機產生的功（Mechanical power、符號P、單位W）分別以實際狀態 下與理想狀態下討論：

(

1-a

)

FVdr b r 4 3 2 _{r r r 1} r = ρπ Ω P 實際狀態下 (2-17-1) dr C ) 1 ( 1 2 1 r r ) ( ) ( 2 r r 3 1 ) ( r ρ λ λ λ r D r L r L C C V C P = + − 理想狀態下 (2-17-2) 得到功之後可求出一個重要的參數指標，功率係數（Power coefficient、 符號C_p）： w p _P P C = (2-18) r 3 r r r r 2 p(r) b (1-a )F d 8 C λ λ λ = 實際狀態下 (2-19-1) 2 3 1 r p(r) _{0.5 V r} P C π ρ = 理想狀態下 (2-19-2) ) arcos(e 2 F -fr r =_π (2-20) φ Rsin r -R 2 B fr = (2-21) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = Ω = R r V r 1 r λ λ (2-22) 其中λ_r則為在r 處的流速比率（speed ratio）， F為Prandtl 修正動量理 論後推導出的葉尖損失因子（Tip-lose factor）。

2.2.4 葉片元素動量理論 1865 年 Rankine 用線性動量理論形成了螺旋槳的簡單模型，由於理論 過於簡化與實際風力機葉片狀況不相符合，為了更能有效的討論葉片上的 性能參數，產生了葉片元素動量理論，主要基於每一個單位長度葉片元素 討論，且葉片元素之間不互相影響，並作用於葉片元素不同位置之力都相 同，亦即不考慮三維效應及假設葉片為無窮多片。 根據葉片元素動量理論，並導入質量守恆定律、能量守恆定律、及白 努力定律，葉片上的軸向推力、扭矩可以下式表達之：

(

V -V

)

4 V a (1-a )Frdr m T 2 _{r r r} 2 1 r = & = πρ 1 (2-23) dr r )F a -(1 b V 4 r m M 3 r r r 1 2 r = &Ω = πρ Ω (2-24) 式中的m_& 為空氣單位時間的質流率，將上列二式與葉片受力計算的軸 向推力、扭矩整理後可以分別得到a、b的線性關係式：

(

)

(

)

(

L() r D(r) r

)

r 2 r 2 r r r r r r 2 T(r) r sin C cos C rsin 8 BC a -1 F a F a -1 rdr )F a -(1 a V 4 BF T 1 φ φ φ π πρ + = ⇒ = = (2-25)

(

L(r) r D(r) r

)

r r r r 3 r r r 1 D(r) r cos C -sin C rsin2 4 BC b 1 F b dr r )F a -(1 b V 4 r BF M φ φ φ π πρ = + ⇒ Ω = = (2-26) 整理以上兩式之後可以得到a與b的關係式又稱能量方程式：

) a F -(1 a 1) (b b 2 _{r r r} r r r + λ = (2-27) 翼弦長度、功及功率係數的計算式則為：

(

)

(

)

r L(r) r 2 2 r r r r r r BC r cos sin a -1 F a -1 F a 8 C φ φ π = (2-28)

(

1-a

)

FVdr b r 4 3 2 _{r r r 1} r = ρπ Ω P (2-29) r 3 r r r r 2 p(r) b (1-a )F d 8 C λ λ λ = (2-30) 2.2.5 貝茲（Betz）理論 Betz 理論首先做了以下的假設：(1)氣流通過葉片的流動假設為一個單 元流管。(2)葉片無錐角、傾角、偏角，看似一個平面轉盤。(3)葉片旋轉時 無摩擦阻力，且無渦流效應影響，氣流通過前跟通過後的靜壓相等 （P₁ =P₂）。(4)推力均勻的分佈在葉片上。(5)通過的氣流為均勻、穩定、且 不可壓縮。 相關參數如圖16 說明，在 Batz 理論基礎下引用白努力定律（Bemoulli） 得出一個關係式如下： a 2 1 2 1 V P 2 1 P V 2 1_{ρ + = ρ + (2-31)} b 2 2 2 2 V P 2 1 P V 2 1_{ρ + = ρ + (2-32)} 因此求解出三者速度V₁、V₂、與V的關係式如下： 2 V V V₌ 1+ 2 _(2-33) 根據動量理論找出V₁、V₂之間的關係：

a) -(1 V V= ₁ (2-34) 2a) -(1 V V₂ = ₁ (2-35) 由上式中可以判斷，當 2 1 a= 時，V₂ =0，此與事實不相符合；當 2 1 a> 時， 0 V₂ < ，此也與實際情形不相符合；因此，在正常狀態下軸向干擾因子 2 1 a< 。 將風機功率計算是以環境風速V₁改寫成為下式： a)A -a(1 V 2 2 V -2 V m P 3 1 2 2 2 1 ρ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = & (2-36) 最大風機功率P發生在 0 da dP = 時，

(

)

0 3a 4a -1 AV 2 da dP 3 2 1 = + = ρ _(2-37) 風機功率P的極值發生在 3 1 a= 或a=1時，而前述已導出在正常狀態下軸 向干擾因子 2 1 a< ，因此最大風機功率發生在 3 1 a= ，帶入功率係數可以求出： 0.593 27 16 AV 0.5 P C ₂ 1 p ≅ = = ρ _(2-38) 功率係數C_p =0.593也就是著名的貝茲極限（Betz limit）。 2.3 葉片幾何外型介紹 2.3.1 翼剖面/翼型（airfoil）外型的選擇 翼剖面的幾何外型，決定了葉片所能產生的升立及阻力的大小，傳統 的翼型為美國太空總署所命名的NACA 系列，以 4、5、或 6 位數字表達翼

剖面的幾何外型，例如NACA 4 4 1 5（圖 17）[63]： 1.最大弧高為 4%弦長 2.最大弧高發生位置距前緣 0.4 弦長 3.厚度為 15%的弦長 由於傳統的翼型已漸漸無法滿足使用的需求，因而各國陸續從80 年 代，開發出各種新翼型[64]~[68]：

1.SERI 翼型 from National Renewable Energy Laboratory ,美國

此類翼型主要是針對各種不同直徑的風力機所做的翼型設計，主要 是由 Tangler 和 Somers 所設計的，應用於直徑 10 至 30 公尺的風力機上 葉根（ 0.4

Rr = ）部份採用SERI S807、中間部份採用 SERI S805A、葉尖

（ 0.95

Rr = ）部分採用SERI S806A，應用於 21 至 35 公尺的風力機上葉

根部份採用 SERI S814、中間部份採用 SERI S813、葉尖部分採用 SERI S812，在大尺寸（直徑 36 公尺以上）風力機方面設計出在葉根部份的 SERI S818、中間部份的 SERI S817、葉尖部分的 SERI S816。

2.NREL 翼型 from National Renewable Energy Laboratory ,美國

此類翼型同SERI 一樣是由美國可再生能源實驗室所研發的翼型， 主要分為薄型和厚型，分別應用在中型葉片的設計和大型葉片設計上， 此類翼型能有效降低由灰塵、昆蟲殘骸等對葉面造成的表面粗度增加而 造成的性能減小。

3.RISΦ-A 和 RISΦ-B 翼型 from RISΦ National Laboratory ,丹麥

此兩種是由丹麥RISΦ 國家實驗室所研究設計的，RISΦ-A 主要的設 計有七種，最佳攻角在 10°，其幾何特徵為前緣呈現尖銳，因此具有對 前緣粗糙度的不敏感性，兩種翼型表示方法通常採用 RISΦ-X-xx，最後 兩位xx 表示最大厚度與弦長之比值。

4.FFA-W 翼型 from Flygtekniska Forsoksanstalten Aeronautical Research Institute of Sweden ,瑞典 此型是由瑞典航空研究所研發設計的，主要有三種系列，編號以 FFA-W1-xxx、FFA-W2-xxx、FFA-W3-xxx 表示，FFA-W1 的薄翼型對於 灰塵、昆蟲殘骸造成的前緣粗糙具有良好的性能，即使在粗糙的前緣下 仍具有較高的升阻比，而FFA-W2 系列的翼型對於平滑和粗糙的表面都 具有良好的性能。

5.DU 翼型 from Delft University ,荷蘭

此類翼型的編號方式為DU xx-W-xxx，前面數字 xx 代表發展的年 份，後面數字xxx 代表最大厚度與弦長之比值。

6.AH 翼型 from Institute for Areodynamics and Gasdynamics of the Univer- sity of Stuttgart ,德國

此類翼型的編號方式為AH xx-W-xxx，前面數字 xx 代表發展的年 份，後面數字xxx 代表最大厚度與弦長之比值。

傳統的翼型設計從葉片根部到葉片尖部皆採用相同的翼剖面，雖然可 以減少製作上的困難以及設計計算上的複雜，但是性能上卻未必能達到理 想狀態，在近期研發的翼型中可以見到翼剖面在徑向方向是有所變化的， 不僅風力機專用的翼型大大提升了性能，在徑向方向做翼型變化更是解決 了空氣動力上失速的問題以及改善了雷諾數對葉片的影響。 2.3.2 旋轉最大外徑（符號D、單位m） 選定尺寸是相當重要的，不僅決定了結構設計、輸出功率、耗資成本 外，甚至大大影響可行性的分析，依據各國的參考資料大概可以統計出尺 寸對應輸出功率的範圍為何，茲將表列說明如表2[69] 通常也可以使用計算來選定旋轉最大外徑尺寸，依據所需要輸出的功 率（符號P_e、單位W）來估算： m g p 2 3 1 e C 4 D V 2 1 P ρ π _⎟⎟ η η ⎠ ⎞ ⎜⎜ ⎝ ⎛ = (2-39) 其中η_g為發動機產生的損失效率，η_m為傳動時損失等等，C_p在此處以 0.4 估算。 2.3.3 選擇翼尖流速比率λ 本節在說明一次λ的定義，以數學式來解釋可以說[葉片頂端的速度] / [環境風速]：

1 R V U = λ (2-40) λ在葉片設計中為最重要的參數之ㄧ，其所影響的參數眾多，茲將介紹 幾個比較重要的影響。 1.λ影響了葉片的轉速（符號N、單位 rpm），通常以下面的數學式估算 之： π π λ Ω = =60 30 N 1 D V _(2-41) 2.λ設計的好壞影響了功率係數（CP），太小的λ無法產生足夠的功率， 過大的λ亦可能造成翼尖失速而產生效率下降，CP與λ的關係可參考圖 18。 3.λ越高則渦流效應造成的誘導速度影響就會減少，徑向誘導因子_b的影 響也就降低了。 2.3.4 葉片數量（符號B） 決定葉片的數量基本上與λ值有直接的關係，通常使用下式計算： 2 80 B λ = (2-42) 除了以公式計算之外，尚有一個方法可以參考，上一小節提到λ與C_P之 間的關係，其實在仔細討論其中之後，可以加入B一起探討，三者之間的關 係，可用圖19[70]來表示，圖中可以看出，不同的葉片數量的最高C_P所對 應的λ值均不相同。

2.3.5 葉片寬度（翼弦長度）（符號C、單位m） 當選擇了適當的翼剖面後，就可依照所需要的翼弦長度而做等比率縮 放，在此提供一經驗公式來幫助尺寸設計，分別考慮實際狀態以及理想狀 態下，可以下列公式計算：

(

)

(

)

r L(r) r 2 2 r r r r r r BC r cos sin a -1 F a -1 F a 8 C φ φ π = 實際狀態 (2-43-1) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = r R B R 9 16 C_{r 2} λ π 理想狀態 (2-43-2) 式中所有參數值都已在上述步驟中求得，代入後葉片寬度則可計算出 來，設計者不一定要完全一樣的採用計算出來的值，實際設計時仍要考慮， 實際受力情形以及製造加工可行性，必要時可將各徑向分段計算出的值， 做一個線性化的處理。 2.3.6 葉片扭轉角（符號β、單位度°） 在前節的理論分析中已經說明了α、β、及φ三者的關係，以數學式表 示則為： α φ β = − (2-44) 由此可知此三者是不可分開來討論的，攻角絕對性的影響C_L、C_D，而β 設計緊繫著α，當φ值計算出來後，只要設計出良好的α值，即可計算出β所 需的值。 分別以實際狀態下與理想狀態下討論φ值的數學計算式：

r r r r _{1 b} a -1 1 tan + = λ φ 實際狀態 (2-45-1) r R 3 2 tan _r λ φ = 理想狀態 (2-45-1) 從翼剖面的設計中可以得到α、CL、CD三者之間的關係，圖20 為典型 的曲線。 由定義中我們知道CL帶給風力機正面的效應，CD帶給風力機負面的效 應，因而必須在CL與CD之間找到ㄧ個α的平衡值，當CL與CD相差的倍數最 高時，所對應的攻角我們稱之為最佳攻角αopt。 從翼型的設計資料中可以看出升力係數CL與阻力係數CD不僅僅與攻角 α 相關，由於流體是粘性的，因此必須將雷諾數（Reynolds number）考慮 進去，葉片上的雷諾數在每個徑向位置有所不同，且每個雷諾數所對應的 最佳攻角值也不一樣，雷諾數可用下列公式估算之： νr r r C W Re = (2-46) ν 即為氣體的動黏滯係數（Kinematic viscosity、單位 m2/s），在常溫常 壓（1atm、20℃）下為 1.51x10-5 m2/s，從數學式可以明顯看出雷諾數與徑 向位置關係，因此在選取最佳攻角時必須先將各徑向位置的雷諾數求出， 在根據翼型參考數據，找出各個雷諾數所對應的升力、阻力係數與攻角的 曲線，在適用的雷諾數區線上找出最大升阻北所對應的攻角值，則稱該值 為各徑向位置的最佳攻角。 當最佳攻角值設計後，β值在徑向位置因α 及φ值的影響而產生非線性

的變化，經過計算後可能造成加工上的困難、製作上的不便，因而與翼展 長同樣的可將計算的結果線性化後，在進行優化設計，不過值得注意的一 點是，改變計算出的β值勢必影響了原先α 值的設計，必須反算α 值確定線 性化後CL與CD的比值是否仍在設計規格範圍內。 2.3.7 徑向位置r處效率因子最佳化 由於風通過葉片時一定會有損失，因而在軸向方向以一個a值作為修 正，葉片在轉動時會產生渦流現象，因而在徑向方向會引發一個誘導速度， 以b值作為修正，以理想狀態下可用Batz 理論推倒出 3 1 a= ，並且沒有b值的 影響，然而a、b確實存在且影響甚大，並且隨著徑向方向而變化，整理前 幾節的數學式，將約束非線性最佳化規劃如下： ［設計變數（design variable）］：

[

r r

]

i a b x = (2-47) ［目標函數（objective function）］： Max.

( )

3 r r r r 2 r r p(r) 8 _{b (1-a )F} d C x G λ λ λ = ：　 (2-48) ［限制條件（constrain）］：經動量理論可以推導出a、b的關係式 ) a F -(1 a 1) (b b 2 _{r r r} r r r + λ = (2-49) ［設計變數上、下限（up、low boundary）］： b 0 0.5 a 0 x x x U i i L i ≤ ≤ ≤ ≤ ≤ (2-50)

三、葉片幾何模型計算及建立 3.1 葉片幾何模型計算流程. 整體計算流程可以規劃出以下步驟（流程圖如圖21）： 1. 根據所需的輸出功率P_e，選定適合的額定風速V₁。 2. 選定適合的翼型。 3. 估算旋轉所需最大外徑D，利用計算式： m g p 2 3 1 e C 4 D V 2 1 P ρ π _⎟⎟ η η ⎠ ⎞ ⎜⎜ ⎝ ⎛ = (3-1) 其中P_e、V₁由(1)中得知，空氣密度為已知ρ=1.225㎏/m3，損失效率 g η 、η_m分別假設為η_g =0.9、η_m =0.85，功率係數依經驗設定C_p =0.4。 4. 依據選用的葉片數量B，設計適用的翼尖流速比率λ。 5. 將葉片在徑向位置等比率分段，大約每 50 ㎜一段。 6. 各徑向位置之效率因子最佳化求解a、b值， ［設計變數］：

[

r r

]

i a b x = (3-2) ［目標函數］： min.

( )

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ₃ r r r r 2 r r p(r) _- 8 _{b (1-a )F} d C : x G λ λ λ (3-3) ［限制條件］：

) a F -(1 a 1) (b b 2 _{r r r} r r r + λ = (3-4) ［設計變數上、下限］： b 0 0.5 a 0 x x x U i i L i ≤ ≤ ≤ ≤ ≤ (3-5) 式中的λ於(4)中確認，各逕向的λ_r由公式： ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = Ω = R r V r 1 r λ λ (3-6) 計算之，葉尖損失因子F及參數f經由下式計算得知： ) arcos(e 2 F -fr r =_π (3-7) r r _Rsin r -R 2 B f φ = (3-8) 各徑向位置的傾斜角φ_r，由下式計算： r r r r _{1 b} a -1 1 tan + = λ φ (3-9) 7. 計算各徑向位置之葉片寬度C_r。

(

)

(

)

r L(r) r 2 2 r r r r r r BC r cos sin a -1 F a -1 F a 8 C φ φ π = (3-10) 其中C_L(r)在(2)中確定，B於(4)中選定，a_r、F_r於(6)中計算得知。 8. 計算各徑向位置最佳攻角α_r。 9. 計算各徑向位置之葉片扭轉角β_r。 r r r φ α β = − (3-11) 10. 線性化徑向位置之葉片寬度C_r、及葉片扭轉角β_r。 11. 完成葉片模型計算。

3.2 葉片氣動性能計算流程 經過葉片寬度C_r、扭轉角β_r線性化計算後，以及額定風速V₁不一定一 直維持在設計狀態下，為了求出在各種額定風速，以及修正後的C_r、β_r下 的各徑向位置的氣動性能，必須先將在偏離設計情形下的干涉因子a、b 值，以迭代方式求出，求解過程順序如下： 1. 首先給定初始a、b值，利用下列式子求出傾斜角φ_r： r r r r b 1 a -1 1 tan + = λ φ (3-12) 2. 將求得的傾斜角φr，以及各徑向位置線性化設計的葉片寬度Cr，代入 下列二式，分別求出a、b值：

(

L(r) r D(r) r

)

r 2 r r r _{C cos C sin} rsin 8 BC a -1 a _{φ φ} φ π + = (3-13)

(

L(r) r D(r) r

)

r r r r _{C sin -C cos} rsin2 4 BC b -1 b _{φ φ} φ π = (3-14) L C 、C_D可於翼型資料中根據設計攻角找出相對的值，攻角α 利用線性化 計算後的葉片扭轉角β_r以下式求出： r r r φ α β = − (3-15) 3. 將第(2)步驟中的a₍₂₎、b₍₂₎代入第(1)步驟，依序迭代計算n次 4. 檢驗是否收歛： a (n) (n) 1) (n a a -a ε = + _(3-16)

b (n) (n) 1) (n b b -b ε = + _(3-17) 當 -6 a =10 ε 、 -6 b =10 ε 判定為收斂。 根據上述方法求得的各徑向位置干涉因子後，可求出相對徑向位置 的軸向推力、扭矩、功等，利用下面三式進行計算：

(

V -V

)

4 V a (1-a )Frdr m T 2 _{r r r} 2 1 r = & = πρ 1 (3-18) dr r )F a -(1 b V 4 r m M 3 r r r 1 2 r = &Ω = πρ Ω (3-19)

(

1-a

)

FVdr b r 4 M 3 2 _{r r r 1} r r =Ω = ρπ Ω P (3-20) 軸向推力係數、扭矩係數、功能係數： 2 2 1 r T(r) _{0.5 dr} T π ρV C = (3-21) 3 3 1 r M(r) _{0.5 dr} M π ρV C = (3-22) r 3 r r r r 2 r p(r) b (1-a )F d 8 C λ λ λ = (3-23) 葉片的總軸向推力、總扭矩、總功可經由積分求得：

∫

=R 0 r dT T (3-24)

∫

=R 0 r dM M (3-25)

∫

=R 0 r dP P (3-26) 依據不同的翼尖流速比率λ，可求出λ對T、M、P的關係；也可以計 算在某個額定轉速N下，不同風速V₁所對應的T、M、P關係；或是相同風

速V₁下，不同額定轉速N所對應的T、M、P關係；或者依所需要觀察的相 關性能曲線，做相互處理比對，以繪出各種性能曲線。 3.3 實例計算 1. 設計所需的輸出功率P_e為1 KW，依據台灣的風能情況，設定額定風 速V₁為10 m/s。 2. 翼型選擇方面參考市售機型＂Whisper H40＂所使用的 NREL S822、 S823，以上兩型翼型為較早期開發出來的，適合應用在 3 至 10 公尺 的風力機；本計劃預計將葉片長度設計在3 公尺以內，依據美國可再 生能源實驗室所研發的翼型NREL 的研究報告說明，在 1 至 3 公尺以 內的小型風力機葉片適合使用的型號為S833、S834、S835，此三型 分別使用的區域如圖22 說明。各型號之翼型曲線如圖 23 所示。 3. 計算旋轉所需最大外徑D： m g p 2 3 1 e ₄ C D V 2 1 P ρ π _⎟⎟ η η ⎠ ⎞ ⎜⎜ ⎝ ⎛ = (3-27) 其中P_e =1000 W、V₁ =10 m/s、ρ=1.225 ㎏/m3、η_g =0.9、η_m =0.85 0.4 C_p = 。 85 . 0 9 . 0 4 . 0 4 D 10 .225 1 2 1 000 1 2 3 _{× ×} ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ × × = π (3-28) 6064 . 2 85 . 0 9 . 0 4 . 0 10 225 . 1 4 2 1000 2 1 3 ⎟_⎠ = ⎞ ⎜ ⎝ ⎛ × × × × × × × = ⇒ π D (3-29) 2.6 公尺。

4. 風力機的型式採用一般大型風力發動機所設計的水平軸式上風型發 力機，偶數葉片易造成旋轉機構產生共振，葉片數量過多也將造成成 本提高、製作困難、品質控制不易，一般來說葉片數量3 能產生足夠 的功能係數CP，且不會造成成本的過度提升，最有幫助的是能從三葉 片的設計過程中吸取經驗，以提昇未來對於大型風力機設計能力；翼 尖流速比率的選擇在3 葉片、額定風速V₁ =10 m/s 的情況下宜選擇 6 = λ 。 5. 葉片最大外徑為 2.6 公尺，考慮葉片中心軸榖（Hub）以及組合介面 需求，因此，單一葉片長度為1 公尺，吾人在徑向方向等分為 20 份， 大約每50 ㎜一段，分段的情況說明如圖 24 所示。 6. 本文最佳化程式採用市售套裝軟體 Matlab 內建的最佳化程式庫中的 fmincon 指令（程式摘要如附錄 main.m

、

myfun.m

、

mycon.m），在

Matlab by fmincon 使用的方法是是序列二次規劃法（SQP）[71] [72]， 其基本理論就是將約束化問題（Constrain Question）轉化為求解一系 列的二次規劃問題，也就是每次迭代步驟中解決一個二次規劃子問 題，迭代時，收斂經由使用近似牛頓法（Quasi-Newton）得到的 Langrange 函數構成的赫世矩陣（Hessian Matrix）來保證，從而轉化 二次規劃子問題，其解產生搜索過程的搜索方向，以數學式表示，當 給定

(x) g f(x) ) L(x, m _i 1 i i

∑

= + = λ λ (3-30) 通過線性化非線性約束條件得到二次規劃問題如下： d ) f(x d H d 2 1 R d min _T k k T n +∇ ∈ (3-31) 0 (x) g d (x) g T _i i + = ∇ i=1,...,m_e (3-32) 0 (x) g d (x) g T _i i + ≤ ∇ i=m_e+1,...,m (3-33) （其中H矩陣的更新，於fmincon中採用BFGS方法） 此問題再通過二次規劃解，形成新的迭代方程 k k k 1 k x d x ₊ = +α (3-34) 此法僅在可行區域（Feasible Domain）求解，因此，求解一個非 線性約束化問題比求解無約束化問題，所迭代的次數更少。對於各徑 向位置之效率因子a、b值最佳化求解過程規劃如下： ［設計變數］：

[

r r

]

i a b x = (3-35) ［目標函數］： min.

( )

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ₃ r r r r 2 r r p(r) _- 8 _{b (1-a )F} d C : x G λ λ λ (3-36) ［限制條件］： ) a F -(1 a 1) (b b 2 _{r r r} r r r + λ = (3-37) ［設計變數上、下限］：

b 0 0.5 a 0 x x x U i i L i ≤ ≤ ≤ ≤ ≤ (3-38) 7. 經過最佳化求解後可以得出各徑向位置的a_r、F_r等值，葉片寬度C_r利 用下式計算：

(

)

(

)

r L(r) r 2 2 r r r r r r BC r cos sin a -1 F a -1 F a 8 C φ φ π = (3-39) 其中C_L在NREL 公司翼型的報告中可以查得。 8. 計算各徑向位置最佳攻角α_r。 最佳攻角值αopt產生在當升阻比最大時，尋找的方法為，在座標原點 做”升力對阻力曲線”的切線，切點處所對應的升阻比為最大值，由該 切點處尋找”升力對攻角曲線”所對應的攻角值，即為最佳攻角，如圖 25 說明。 各徑向位置的雷諾數不相同，利用公式計算之： νr r r C W Re = (3-40) 將計算得到的雷諾數以內插或外插，即可求出各徑向位置的最佳攻 角。 9. 利用最佳化後得出的a_r、b_r值計算出各徑向位置的傾斜角φ_r，以及上 一步驟獲得的各徑向位置的最佳攻角αr，計算各徑向位置之葉片扭轉 角β_r，如下式： r r r φ α β = − (3-41)

經程式演算後，結果如表3，步驟 6 至步驟 9 計算過程的 Matlab 程 式碼如附錄一。 10. 線性化徑向位置之葉片寬度C_r、及葉片扭轉角β_r。 以徑向位置r分別當作葉片寬度C_r、及葉片扭轉角β_r的函數，假設 為四次多項試，如下式所列：

( )

c c c c c r r a x b x c x d x e C = 4 + 3+ 2 + 1 + (3-42)

( )

β β β β β βr r =a x +b x +c x +d x +e 1 2 3 4 _(3-43) 本文利用Matlab內建指令polyfit（程式摘要如附錄二）分別求出兩 式的a、b、c、d、e係數，考慮葉片形狀分佈、組裝界面處等，忽 略葉柄部分（ R R r 2 . 0 ~ 0 = ）及葉尖部分（ R R r 1 ~ 9 . 0 = ）的線性化處理， 因此僅針對 R R r 9 . 0 ~ 25 . 0 = 做分析計算，結果如下所列：

( )

₌₋0.56 4 ₊1.684 3 ₋1.403 2 ₋0.168 1₊0.642 x x x x r C_r (3-44)

( )

₌24.625 4 ₋104.793 3 ₊169.159 2 ₊₋131.973 1 ₊48.626 x x x x r r β (3-45) 依據上式結果將r帶入，重新求出線性化後的C_r、及β_r的值，茲將 結果與比較詳列於表4，C_r對站位及β_r對站位曲線如圖26。 11. 完成葉片設計。 3.4 三維立體模組建立 本文利用三維機械繪圖軟體SolidWorks 繪製三維葉片模組，翼型曲線

經由NREL 報告中獲得 c x 、 c y 資訊，將前節求出的C_r帶入，獲得各徑向位 置的xr、yr座標，整理繪圖時所需的相關數值資料xr、yr、βr於同一檔案

airfoil.txt，利用 SolidWorks 內建的巨集功能，以及附加程式 API 語言的 「Open For Input As 打開文件、讀入指令」、「boolstatus =

Part.Extension.SelectByID2 在平面繪製草圖指令」、「Part.SketchSpline 繪 製不規則曲線指令」、「Part.Extension.RotateOrCopy 旋轉不規則曲線指 令」、「boolstatus = Part.Extension.SelectByID2 疊層拉伸指令」等等，規劃 一自動畫繪製葉片程式檔案（程式摘要如附錄三），完成葉片自動化繪製 後，再依據結合介面的需要設計葉柄處，即完成葉片設計，葉片完成示意 圖如圖27。

四、有限元素分析 4.1 自然頻率、自然模態分析 隨著運動速度的提高，對機器或結構動態特性的了解要求日益迫切， 所謂的動態特性指的就是自然頻率(Natural Frequencies)及其相應的振型 (Mode Shape)，本質上固有頻率適度顯示結構體的剛性行為；其與複合材料 板受力產生的位移與應變，成相對應的比例關係。在許多有限元素法或結 構動態分析的文獻裡諸多介紹，結構可能因共振而產生過大之變位或應 力，超過材料之容許極限而產生破壞；或是振動中所產生之反復應力，雖 然其應力值並未超過材料之容許極限，但卻可能因為反復次數過多而造成 材料疲勞破壞；另外振動亦會伴隨產生噪音、結構的顫動，使得使用者產 生不舒服的感覺。對於結構量測低階自然頻率遠比量測其位移與應變容 易，通常低階自然頻率在工程設計較受重視，因低階自然頻率容易受外界 激發而引起結構共振。 4.1.1 ANSYS模型建立步驟 在元素使用方面芯材採用編號solid186，面層的纖維/樹酯採用編號 shell99，分析自然模態採用無拘束下邊界條件，計算特徵植所採用的演算法 選用子空間迭代法（Subspace），實際操作步驟如下段說明。

1. Preprocessor → Element type：選擇心材 solid186，面層 shell99。 2. Preprocessor → Real constant：設定面層實體常數，包含排列角度、厚

度等等。

3. Preprocessor → Material Props → Material Models：設定各材料性質。 4. Preprocessor → Modeling：先使用三維繪圖軟體建立模組，儲存成

ANSYS 可以讀取的格式（如 Parasolid、IGS、Sat 等），以 ANSYS 匯 入葉片三維立體模組。

5. Utility Menu → WorkPlane → Local Coordinate Systems → Create Local CS：建立元素座標，由於纖維材料具有方向性，需先定義座標系統。 6. Preprocessor → MeshTool：選擇元素參數、材料性質、各元素之尺寸

大小，並分割元素。 準備進行模態分析：

7. Solution → Analysis Type → New Analysis：選擇分析型態，自然頻率 模態分析點選“Modal”。

8. Solution → Analysis Type →Analysis Options：設定使用的演算法及要 分析的模態個數，選用Subspace 進行模態分析。

9. Solution → Solve → Current Ls：求解。

10.General Postproc →Results Summary：列出所有的自然頻率。

自然頻率。

12. General Postproc →Plot Results →Contour Plot →Nodal Solu：再選 Nodal Solution →Displacement vector sum，圖示出各節點的位移。 以上步驟可以得到葉片在無拘束條件下的自然頻率及模態。 4.2 暫態振動分析 風力葉片在受風的作用下會產生一個正向的結構承受力以及周向的旋 轉力，葉片必須在風的作用下有足夠的強度，風力與時間的關係是隨機變 化的，有時風力在短時間之內迅速增加二到三倍，甚至數倍，因而暫態分 析將有助於了解葉片在瞬時陣風作用下的位移改變及應力、應變分佈情 形，本節將討論葉片在受瞬時激振力情形下的振動行為分析方式。 4.2.1ANSYS模型建立步驟 暫態振動分析的前處理步驟與自然頻率分析一樣，可參閱4.1.1節的1-6 步驟，由於暫態分析相當耗費計算機系統資源，通常建立的網格比較粗， 用以節省計算時間及系統資源，為了能有效節省並善用資源，且能準確分 析量測點的結果，本文在實驗量測位置的節點上進行網格細分，並提供一 些節省資源的方法使用[73]。 準備進行暫態振動分析：

細分網格，將原先較為粗的網格在需要分析的資料點處進行細分， 可得到較為精確的解。

8.Preprocessor→NumberingCtrls→Compress Numbers→Nodes、Elements： 壓縮節點、元素編號，使用此指令可以避面編號的跳號，此功能可 將已設定的所有相關資料轉換。

9. Preprocessor→Numbering Ctrls→Element Reorder→Reorder by XYZ：將 元素內部編號重新排列，可以使得資料結構完整，計算時更有效率。 10. Solution→Analysis Type→New Analysis：選擇暫態分析點選“Tran-

sient”。

11. Solution → Analysis Type →Sol’n Controls：設定暫態的步進時間、結 尾時間、阻尼值等等。

12. Solution → Define Loads →Apply→Structurl：設定挾持邊界、受風下 葉片所承受的力、以及葉片轉動的速度。

13. Preprocessor→Archive Model→Write：將有限元素分析資料寫到文字 檔案中，本指令可以寫出乾淨的分析模型。

14. Utility Menu>File>Clear & Start New：清除所有資料，並重新開啟一 個乾淨的新檔案。

15. Preprocessor→Archive Model→Read：將之前寫出的乾淨模型讀入。 16. Solution → Solve → Current Ls：求解。

4.2.2ANSYS模擬分析中各參數的取得 1.材料常數的給定 葉片製作所使用的材料包含面層的純樹酯、玻璃纖維/樹酯，芯材的 PS（聚苯乙烯）發泡，材料常數經實驗求得表列於表 5，並依序於 ANSYS 中設定編號為MAT=1、2、3。 2. 元素座標的設定 由於纖維材料具有方向性，必須在各面層的網格中定義坐標系統， 使得在實體常數中設定的角度有所意義，依據不同位置的面層定意坐標 系統，葉片的部分於徑向方向（翼展方向）設定為X，葉片旋轉方向（翼 弦方向）為Y，葉片表面的法向量方向為 Z，而葉柄的部分為圓柱幾何， 軸向設定為 X，弧面切線方向設定為 Y，表面的法向量方向為 Z，設定 完成後在產生元素時即可有規律且整齊的元素，如圖28 所示，各元素的 坐標系統都整齊的排列。 3.力的給定 Kistler 敲擊槌透過訊號分析儀將訊號存錄於個人電腦，利用個人電 腦中的軟體將存錄的力對時間圖取數個資料點紀錄下來，如圖29 所示， 在ANSYS 軟體中使用表單（table）的功能 parameters → Array Parameters →Define/Edit，將紀錄的力對時間資料點對應輸入新增的表單裡，儲存 表單內容，在進行力的設定時 Solution → Define Loads →Apply→ Structurl→Force/Moment 其中的 Apply as 選項裡選擇 Existing table，點選 前一步驟所儲存的表單，完成暫態力的給定。

圖30 為阻尼量測試驗之加速度及位移對時間曲線，加速度為量測葉 片接近尖端50mm 處位置的量測數據，位移根據加速度量測數據對時間 積分後求得；由於 α-damping 對低頻影響較大對高頻沒什麼影響而 β-damping 對高頻影響較大對低頻沒什麼影響，因此，低頻下可以忽略 β-damping，α-damping 則利用第一階頻率計算；圖中得知在 27.4ms 至 274.7ms 內共振動六次，利用下式計算之， t m f Δ = (4-1) 得知第一階頻率f=6/(274.7-24.4)=24.26Hz，本文使用單自由度公式來計 算α-damping，

(

)

m V V_{n n m} π ξ 2 ln + = (4-2) 其中

ξ

為此共振頻率之阻尼比，Vn為第 n 次的位移，Vn+m為第 n+m 次的 位移，計算後得到

ξ

=0.0189，利用下式並忽略低頻時的 β-damping，可 以求得系統α-damping， Ω ≅ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ Ω + Ω = 2 2 2 α β α ξ (4-3) 其中α 為與質量矩陣有關的阻尼比，Ω為模態之角自然頻率，計算後求 得 α= 2x(24.26x2π)x 0.0189=5.76。

進入求解器中進行設定Solution → Analysis Type →Sol’n Controls， 選擇後會出現設定的視窗，點選Transient 設定，在 Damping Coefficients 裡將計算出的α-damping 輸入，Full Transient Option 選擇 Stepped loading。

4.3 有限元素分析驗證

由於複合材料的行為較為複雜，本文先以均勻等向性之材料鋁合金， 驗證分析之可行性，鋁合金採用Al 6061，楊式系數為 6.92e10 Pa，柏松比 為0.33，密度的實測值為 2669.3kg/m3，幾何形狀為一[300 X 29.74 X 4.05]

mm3的矩形平板，首先以第 4.1.1 節介紹的方法進行無拘束條件下的自然頻 率分析，網格使用solid 186 六面體元素，實驗方面如圖 31 所示，鋁板以橡 皮筋與試驗台架接合，加速規放置於鋁板中心處，圖32 量測結果曲線，第 一階到第五階自然頻率依序為235.5、648.5、1266.3、1355.1、1662.3，ANSYS 分析結果依序為235.5、649.4、1273.4、1346.0、1671.8，相關結果如表 6， 其間之最大誤差量不到1%，可以確定模態分析是可行的。 其次驗證在單邊夾持下的自然頻率，實驗方面如圖33 所示，夾持鋁板 一端32.3mm 處，加速規置於中間量測數據，應變規於中間另一面量測應 變，比對實驗與分析在單邊拘束條件下自然頻率的誤差量，實驗量測結果 曲線如圖34，與分析結果比較列於表 7，由表中可看出最大的誤差量為 1.08%，因此可判斷 ANSYS 分析拘束條件下模態是準確的。 驗證了模態分析準確性後，開始進行暫態振動分析，利用4.2.2 節中的 第4 項說明的方法計算鋁板的阻尼，圖 35 為鋁板所量測的加速度對時間曲 線，以及利用加速度對時間積分後計算的位移對時間曲線圖，第一階頻率 f=46.6Hz，模態之角自然頻率 Ω=292.8，位移經 9 次振動後從 V9=9.450e-5 衰減至V18=7.396e-5，計算後得知

ξ

=0.00433，因而得到 α= 2.54。 使用第4.2.1 節所述之方法計算鋁板的加速度與應變對時間的變化，並 與實驗結果比較，圖36 為加速度對時間的比較結果，圖 37 為應變對時間 的比較結果，為了判斷方便，茲將每個波峯及波谷以曲線連接起來，並將 其線性化相互比對，加速度趨勢線如圖38，應變趨勢線如圖 39，由圖中可 以判斷出兩者的變化趨勢相似，波型也近似， 分別在不同時間點下取加速 度及應變的結果值互相進行比較，比較結果詳列於表8，加速度最大誤差在 5.03%，應變最大誤差在 5.99%，可見暫態振動分析之方法是可行且準確的。

五、葉片振動試驗 本文試驗所使用之葉片為一實際量測全長為 998 mm 的複合材料葉 片，單隻葉片實測重量為750g，採用樹酯轉注模塑成型法（Resin Transfer Molding，RTM）製做，RTM 製作方式為目前風力葉片最常使用的方法之 ㄧ，具有成型效率高、不易產生揮發物質、產品表面較光滑、無污染環境 等優點，主要步驟是先將芯材包覆纖維後放入模具中，模具封閉後先行抽 真空，約達0.01tour 時，將樹酯注射進入模具，纖維完全凈潤樹酯後並確定 模具內空氣皆以排出，將模具完全封閉置於室溫熟化。葉片結構部分為三 明治構型，芯材使用PS 發泡，發泡密度達 13.6 kg/m3，面層使用玻璃纖維/ 樹酯雙軸針織纖維包覆，採用的雙軸針織纖維為兩個方向為正交且互相垂 直針織布，每層為0.2mm 共鋪疊五層，葉柄部分使用多層玻璃纖維/樹酯纏 繞，複合材料葉片示意如圖40。 5.1 自然頻率量測實驗 本實驗主要量測葉片在無拘束條件下的自然頻率，將葉片以橡皮筋連 接於試驗台架上，加速規分別黏貼於三處不同位置接收訊號，以敲擊槌激 振，將量測資料儲存後，記錄下其自然頻率。 1. 基本設備有: (1) 訊號分析儀 (2) Kistler 敲擊槌 (3) Isotron Acceleronmeter 加速規 (4) 個人電腦 (5) 葉片 (6) 試驗台架 2.實驗程序

(1) 將敲擊槌及加速規的訊號端組裝於訊號分析儀上。 (2) 將訊號分析儀與個人電腦連結後，打開訊號分析儀之電源，並開啟 電腦應用程式。 (3) 把葉片吊掛於試驗台架上，並將加速規分別於三處不同位置端黏於 葉片上，如圖41 所示。 (4) 於應用程式中開啟訊號接收鈕，以敲擊槌敲及葉片不同位置約三至 五次。 (5) 於應用程式中按下停止紀錄鈕，並呼叫出加速度對頻率曲線，將各 皆頻率記錄下來，量測結果如圖42。 5.2 阻尼量測實驗 本實驗將葉片底部端挾持，葉尖端放置加速規量測，敲擊後第一階振 動行為乃上下擺振，根據量測出的波形對時間圖可以求出第一階頻率，依 據單自由度公式求其阻尼值。 1. 基本設備有: (1) 訊號分析儀 (2) Kistler 敲擊槌 (3) Isotron Acceleronmeter 加速規 (4) 個人電腦 (5) 葉片 (6) 試驗台架 2.實驗程序 (1) 將敲擊槌及加速規的訊號端組裝於訊號分析儀上。 (2) 將訊號分析儀與個人電腦連結後，打開訊號分析儀之電源，並開 啟電腦應用程式。

尖 50mm 位置上。 (4) 於應用程式中開啟訊號接收鈕，以敲擊槌敲擊距離葉柄 545mm 處 位置，如圖 43 所示。 (5) 於應用程式中按下停止紀錄鈕，並呼叫出加速度對頻率曲線，將 各第ㄧ階波型的波峰加速度值及時間記錄下來。 5.3 應變量測實驗 本實驗將葉片底部端挾持，距離葉柄端135mm（編號 P01）和 225mm （編號 P02）處黏貼應變規，如圖 44 所示，應變規坐標沿著葉柄向葉尖之 方向，敲擊葉片後量測其應變對時間之變化，並與ANSYS 分析結果比較。 1. 基本設備有: (1) Ni cDAQ-9172 資料擷取系統 (2) Ni 9237 模組（用於接收應變規數據） (3) Ni 9234 模組（用於接收敲擊槌數據） (4) KYOWA KFG-2-120-C1 應變規 (5) Kistler 敲擊槌 (6) 葉片 (7) 試驗台架 (8) 個人電腦 2.實驗程序 (1) 將 National Instruments 公司的資料擷取系統模組與 Ni 9237 模組 及 Ni 9234 模組裝。 (2) 將應變規接於 Ni 9237 模組，Kistler 敲擊槌接於 Ni 9234 模組，並 將 Ni cDAQ-9172 資料擷取系統與個人電腦連結（模組組裝示意圖 如圖 47）後，打開所有儀器電源，並開啟電腦應用程式，程式軟 體使用 NI 公司所提供之 Ni DAQmax。

(3) 試驗台架將葉片底端挾持固定住，並量測葉片是否確實水平且穩 固。 (4) 於應用程式中開啟訊號接收鈕，以敲擊槌敲擊距離葉柄 545mm 處 位置。 (5) 於應用程式中按下停止紀錄鈕，並呼叫出各點應變對時間曲線， 以及敲擊槌施的力對時間曲線，將各應變、力、及時間資料記錄下 來。 5.4有限元素分析與實驗結果討論 葉片實際量測重量為750g，本文建立的有限元素模型（圖45）重量亦 為750g，暫態分析需要龐大的計算資源，為了能有效節省資源，減少計算 時間，建立的格點沒有很細緻，網格收斂約在千分之ㄧ，共有3245個元素、 4033節點。無拘束條件下分析自然頻率所得的第一階至第五階自然頻率分 別為95.3 Hz、214.6 Hz、247.8 Hz、348.9 Hz、387.5Hz，各階自然模態綜整 如圖46，由模態圖可看出自然振動行為有彎曲振動以及扭轉振動，其中彎 曲振動又包含了翼展方向（flap- wise）和翼弦方向（edge-wise）。 無拘束條件下自然頻率分析結果與實驗結果相比較第一階到第五階的 誤差量約在4.2%以內，詳細結果資料如表9，互相比較實驗結果與分析值的 曲線（如圖47），可以看出有相同的趨勢；由於葉片採用複合材料並以手工 製作，實際葉片與分析模組會有所差異性，如樹酯含量的分佈情形、纖維

分佈密度、纖維排列方向、幾何形狀製作時產生的變異等等，造成葉片分 析與實驗的誤差量較鋁板的為大。 進行暫態分析時先驗證單邊挾持下的自然頻率，分析結果顯示第一階 至第三階分別為27 Hz、50 Hz、110Hz，根據實驗結果顯示（實驗數據如圖 48 所示），測量出的第一階頻率為 25Hz、第三階頻率為 109Hz，比較二者 結果相當近似，得到初步可靠性的驗證。其中第一階的彎曲振動自然頻率 與發電機葉片旋轉頻率息息相關，ㄧ般工程設計上第一階自然頻率不能與 發電機葉片旋轉頻率整數倍重合，必免產生共振，對於三葉片的發電機設 計要求上，第一階自然頻率必須大於發電機葉片旋轉頻率的320%，以本文 葉片範例來看，葉片可應用在旋轉頻率為7.813Hz（亦等於 468.75rpm）以 下的發電機。續而進行暫態分析，選擇距離葉尖50mm 處做為一參考點， 比較分析結果與實驗數據在參考點位置上的加速度對時間的變化，將兩條 結果曲線相互比較於圖49，圖中觀察出兩者曲線相當近似，個別比較幾個 時間特徵點（如圖中所標示之位置），時間資料點之差異值也都在6.6%以 內，阻尼影響振動的收斂行為也近似，因此，可以確定暫態分析模組的準 確性以及阻尼值的適用性。 使用暫態分析比較葉片受瞬時激振力時，應變對時間的變化趨勢，圖 50為應變實驗量測結果，應變量測的兩點（P01、P02）位置與施力為同一 平面，點P01距離夾持端葉柄的位置較近，從實驗結果判斷，兩點應變經敲 擊後立刻同步產生正應變，點P01初始最大應變為1.263e-4大於點P02初始最 大應變為3.379e-5，初步判斷實驗結果符合理論物理現象；進行有線元素分 析時，必須先確認元素的方向性是否與應變規黏貼的方向一致，於總域座 標中葉柄向葉尖之徑向方向（翼展方向）定義為X軸，複合材料的各個元素