基因演算法於二維空蝕水翼逆向設計及最佳設計之應用(3/3)

行政院國家科學委員會專題研究計畫 成果報告

基因演算法於二維空蝕水翼逆向設計及最佳設計之應用

(3/3)

計畫類別： 個別型計畫 計畫編號： NSC91-2611-E-002-004- 執行期間： 91 年 08 月 01 日至 92 年 07 月 31 日 執行單位： 國立臺灣大學工程科學及海洋工程學系暨研究所 計畫主持人： 郭真祥 計畫參與人員： 劉佳龍、簡鴻斌 報告類型： 完整報告 處理方式： 本計畫可公開查詢

中 華 民 國 92 年 11 月 6 日

行政院國家科學委員會補助專題研究計畫成果報告

基因演算法於二維空蝕水翼逆向設計及最佳設計之應用

Application of Genetic Algorithms to Inverse Design and

Optimal Design of Two-Dimensional Cavitating Hydrofoils

計畫類別：■個別型計畫 □整合型計畫

計畫編號：NSC 91-2611-E-002-004

執行期間：89 年 08 月 01 日至 92 年 07 月 31 日

計畫主持人：

郭真祥

計劃參與人員： 劉佳龍

簡鴻斌

執行單位：國立台灣大學工程科學及海洋工程學系

中華民國

92年07月31日

摘 要

二維翼形廣泛地被應用在航空器、船舶或者泵等流體機械由來已久，相 關的研究或設計已累積了大量而可靠的系列翼形資料，例如 NACA 翼形系 列、Eppler 翼形系列等等。在船舶方面，船形設計高速化，螺槳設計上必需 更能減少振動、增加效率、靜音、抗空化等。基於上述原因，本研究主要目 的在於建立一套有效的設計工具以期能夠改變現有的翼形，使之達到性能上 的要求。 在研究方式上，分成兩部分進行，第一部分將二維翼形設計問題分成壓 力分佈最佳化及逆向設計二步驟在勢流場中進行，首先求取最佳的壓力分 佈，將設計目標透過限制條件型態呈現，經由基因演算法搜索符合設計條件 的最佳壓力分佈。壓力分佈的表示方式採用 B 木條曲線(B-spline)來表示。藉 著控制點的移動，使曲線能夠產生局部變形，同時導入基因演算法操作控制 點，讓壓力分佈達到最佳化。其次，將前一步驟所得最佳壓力分佈設為目標 壓力分佈，由該壓力分佈來反算對應的水翼形狀。而第二部分則將勢流場中 的計算延伸到黏性流場，其中以自行建立完成之二維黏性流計算方法出發， 導入以流體體積法為觀念之二相流的數值模型，利用等效流體模型，考慮流 場中每一個控制體積內的兩種流體的分佈，使計算方法得以考慮流場中兩種 不同的流體，再利用壓力、速度和密度之間的偶合關係使流場之計算可趨於 穩定，其中經由疊代方式進行計算，直到片狀空化流場的形狀穩定為止。整 個期末報告以各個部分之進行方法與獲得結果分開撰寫，方便陳述及閱讀。 關鍵詞：二維空蝕水翼、幾何最佳化設計、基因演算法、黏性流計算、片狀 空化流場、二相流法。

Abstract

Marine applications involving attached bubbles on propellers and underwater vehicles are obvious examples of cavitation phenomena. For naval architecture issues, cavitation is always a very important topic relating to high-speed ships. However, the high-speed propulsors including water jet system, super-cavitating propeller and surface-piercing propeller are facing critical challenge from cavitation problems. In this research work we divide the method into two parts. In the first part, we try to improve the pressure distribution of known hydrofoil sections by GAs. And then we link the potential flow calculation method for two-dimensional hydrofoils and GAs to solve the inverse problems with a specified pressure distribution. In the second part, for simulating sheet cavitation flows, the adopted method is starting from a two-dimensional flow code capable of computing turbulent hydrofoil flows. In order to take cavitation region into account, a two-phase flow scheme is further implemented in this code. The proposed approach treats gas and liquid phases as an effective fluid with density varied in space. The flow field is computed in both phases with vapor pressure recovered inside the cavity via a pressure-velocity-density coupling scheme. The computation is iteratively conducted, until the shape of the sheet cavitations become stable. This research is to simulate steady sheet cavitation flows of a two-dimensional hydrofoil using a two-phase approach. The size and shape of the sheet cavitation, and the lift and drag coefficient of a cavitating hydrofoil are then predicted.

Keywords: GAs、Potential flow、Cavitation flow、Two phase flow scheme、 Pressure-velocity-density scheme。

目 錄

摘要

……….….

Abstract

………

目錄

………...………... I

表格目錄

………. IV

圖例目錄

………... V

第一部分

第一章 文獻回顧

1. 二維翼型設計歷史………. 1

第二章 基因演算法與依據可行性行為基礎的逞罰方式

2. 前言……….………. 1

2.1 簡介……….…..……. 1

2.2 由演化論看基因演算法………. 2

2.3 基因演算法流程與運算子………. 2

2.4 編碼與解碼………. 3

2.5 基因運算(Genetic Operator)………. 4

2.6 選擇運算(Selection Operator)…..………. 7

2.7 基因參數………. 8

2.8 結論………. 9

第三章 小板法

3. 前言………. 12

3.1 統御方程式………..………. 12

3.2 邊界條件………. 13

3.3 相關係數………. 15

3.4 數值方法與離散化………. 16

第四章 二階段設計法

4. 前言……….………. 19

4.1 壓力分佈設計………. 19

4.1.1 壓力分佈曲線的建立………. 19

4.1.2 以基因演算法進行壓力分佈設計………. 26

4.2 二維翼形的逆向設計………. 27

第五章 結果與討論

5. 前言………. 28

5.1 二維翼形逆向設計之驗證………..………. 28

5.2 二階段設計之應用例…..………. 31

第六章 結論

第二部分

第一章 研究目的及文獻回顧

第二章 理論基礎

2.1 統御方程式………..…..……. 39

2.1.1 Reynolds –Averaged Navier-Stokes 方程式…..……. 39

2.1.2 連續方程式及體積比率方程式………. 39

2.2 空化模型………. 40

2.3 紊流模型………. 40

2.4 邊界條件………. 41

2.4.1 入口條件………. 42

2.4.2 出口條件………. 42

2.4.3 壁面函數………. 42

第三章 數值方法

3.1 有限體積法………. 43

3.1.1 交錯放置………. 43

3.1.2 非交錯放置………. 44

3.1.3 有限體積法………. 44

3.2 壓力修正方程式………. 47

3.3 數值網格產生及流場邊界條件………. 48

3.4 流場求解流程及參數設定………. 51

3.4.1 流場求解流場………. 51

3.4.2 計算流程………. 51

3.4.3 參數設定………. 52

第四章 計算結果與討論

4.1 空化長度之定義………. 53

4.2 實驗結果比較……...………. 55

4.3 理論結果比較………. 56

4.4 密度變化………. 57

4.4.1 攻角改變………. 57

4.4.2 空化係數改變………. 57

第五章 結論與未來展望

5.1 結論………. 59

5.2 未來展望……...……….. 59

參考文獻

62

表格目錄

第一部分

表

4.1 壓力分佈設計之基因參數………..……….. 27

表

4.2 翼形逆向設計之基因參數……… 27

表

2.1

k-ε

模式的經驗常數

………....……… 41

表

_{3.1 積分形式的數值近似……… 46}

表

_{3.2 方程式離散後各項之結果……… 46}

表

3.3 二維翼形流場邊界條件……… 50

第二部分

表

4.1 攻角為 2.3 空化係數為 0.734 之比較表格……….. 55

表

4.2 攻角為 4.0 空化係數為 1.0 之比較表格……… 55

圖例目錄

第一部分

圖

_{2.1 基因演算法流程圖……… 3}

圖

2.2 解碼與編碼關係圖……… 4

圖

2.3 單點交換……… 5

圖

2.4 兩點交換……….……... 6

圖

2.5 均勻交換………..…. 6

圖

2.6 突變………..……… 7

圖

2.7 輪盤選擇示意圖……….. 8

圖

3.1 二維翼形勢流場示意圖……… 13

圖

3.2 速度勢與邊界 S 的關係圖………..……… 14

圖

4.1 壓力分佈建構圖……….……..……… 21

圖

4.2 壓力分佈最佳化流程圖…….………..……… 24

圖

4.3 二維翼形建構圖………..………. 26

圖

_{4.4 二維翼形反向設計流程圖……… 26}

圖

5.1 NACA

4412 與初始翼形之比較……….…..……..……… 28

圖

5.2 NACA

4412 與初始翼形壓力分佈之比較………. 29

圖

5.3 NACA

4412 與逆向設計結果之比較…….………..………..… 29

圖

5.4 NACA

4412 與逆向設計結果壓力分佈之比較……….… 30

圖

5.5 NACA

2410 與逆向設計結果之比較……….…..…….. 30

圖

5.6 NACA

2410 與逆向設計結果壓力分佈之比較……….... 31

圖

5.7 NACA

4412

壓力分佈與滿足設計條件之壓力分佈之比較…. 32

圖

5.8 NACA

4412 與限制條件下設計之翼形之比較………. 32

圖

5.9 逆向求取 NACA 4412 翼形之適應值對世代數圖………. 33

圖

5.10 逆向求取 NACA 2410 翼形之適應值對世代數圖………. 33

圖

5.11 壓力分佈設計之適應值對世代數圖……… 34

圖

5.12 逆向求取滿足設計條件翼形之適應值對世代數圖……… 34

圖

6.1 二維翼形最佳化流程圖……… 36

第二部分

圖

3.1 質量守恆時的控制體積（左）

，u 方向動量守恆時的控制體積

（中）

，v 方向動量守恆時的控制體積（右）………

43

圖

3.2 速度與壓力皆放置於同一網格內…………..………44

圖

3.3 網格點與單元中心之關係………45

圖

3.4 單元節點 P 與四周單元節點(W，E，N，S)關係………….……...47

圖

3.5 網格分佈形式………..…. 48

圖

3.6 靠近壁面的放大圖……… 49

圖

3.7 二維流場計算的範圍.……….. 50

圖

3.8 流場求解流程……… 51

圖

4.1 NACA

66(mod)攻角 4 度的 Cp 圖……….. 53

圖

4.2 放大 A 點的 Volume Fraction 圖………. 54

圖

4.3 Volume Fraction 等值線圖……… 54

圖

4.4 攻角為 2.3 空化係數為 0.734 之比較圖……….. 55

圖

4.5 攻角為 4.0 空化係數為 1.0 之比較圖……… 55

圖

4.6 攻角為 4.0 空化係數為 0.84 之比較圖……… 56

圖

4.7 攻角為 4.0 空化係數為 0.91 之比較圖……… 56

圖

4.8 NACA 0012 改變

α

的

ρ

m

圖………..……..……… 57

圖

4.9 NACA 0012 改變

σ

的

ρ

m

圖………. 58

圖

5.1 （上）本文之計算網格，（下）加密後之計算網格………… 60

第一部分

第一章 文獻回顧

1.

二維翼形設計歷史

關於二維翼形設計相關研究，在國內外均有論文發表。在逆向設計方面 1945 年 Lighthill[1]以保角轉換方法（conformal mapping）設計滿足壓力分 佈之二維翼形。1987 年 Giles、Drela[2]以有限體積法(finite volume method) 求解尤拉方程式（Euler equation），再以牛頓法決定滿足壓力分佈之二維翼 形。1996 年 De Flaco、Del Balio、Della Chioppa 及 Tarantino[3]以發展策 略（evolution strategies, ESs）與基因演算法（genetic algorithms, GAs） 結合而成的培育者基因演算法（breeder genetic algorithms, BGAs）對已知 的壓力分佈進行反向設計。1996 年國內辛敬業教授[4]採用積分翼表面對應 之速度場，從而得到速度勢。再以邊界元素法求解特殊邊界條件之格林函數 （green function）求得設計之二維翼面。

在翼形最佳化方面，1990 年 van den Dam、van Eanmond、及 Slooff[5][6] 以勢流理論描述壓力分佈，再以梯度方法求解最佳化壓力分佈，從而得到最 佳化翼形。1995 年 De Flaco[7]等人以平行化基因演算法，對升力係數、阻 力 係 數 等 結 合 成 的 適 應 函 數 ， 求 得 最 佳 化 的 翼 形 。 同 年 Obayashi 及 Takanashi[8]基因演算法結合逆向設計求得最佳化翼形。在國內則有潘大知 教授及邱永裕碩士[9]以平行化基因演算法以直接設計法對二維翼形進行最 佳化設計。

第二章 基因演算法與依據可行性為基礎的懲罰方式

2.

前言

本章將分成兩部分討論。第一部份介紹基因演算法，第二部分則介紹「依 據可行性為基礎的懲罰方式」。

第一部分 基因演算法

2.1.

簡介

基因演算法是 1975 年 John Holland 與其同事及學生於密西根大學發展 出的一種搜尋方式。基因演算法是建構於遺傳變異(genetics)與自然界天擇 (natural selection)機制的演化型搜尋方式。在計算過程中，問題的解藉由 0 與 1 編成的二進位碼 (binary code) 模擬基因組成染色體(chromosome)，並 且每一組二進位碼可視為單一的人工生物個體。數個個體則組成生物族群 (population)。生物族群在選擇(selection)與基因操作(genetic operation)如 互換(crossover)，突變(mutation)下產生下一個世代(generation)，代代逐 漸演化，而得到最佳解。

2.2.

由演化論看基因演算法

1859 年達爾文(Charles Darwin) 出版《物種起源》的初版，正式向公 眾推介其理念。其中最重要的概念為「物競天擇，適者生存」（此概念又名 為天擇說）。其中敘述在資源有限的環境下，生態族體中的每一個生物個體， 為了生存，而必須相互競爭。在競爭的過程中，適應程度較差的個體所爭取 到的資源，相對少於適應程度較佳的個體。適應程度不佳者也較容易被所處 環境中的選擇機制淘汰，成為競爭的失敗者。適應程度佳的個體因為能獲得 較多的資源，其繁殖子代的機會較多，有利於其本身留下自己的基因。 1950 年代經由遺傳學的進步，世人得以瞭解遺傳機制在演化中扮演的角 色：在繁衍的過程中染色體上攜帶不同的基因組合。透過有性生殖的方式， 染色體隨機交換，讓個體間隨機交換基因組合，產生新的子代。另外，基因 突變的機制，可以使基因組合多樣化，增加適應環境的機會。 因此，由以上敘述可得知兩點： (一) 、在選擇運算下，基因演算法的計算結果總是朝向更適合目標的方 向演化。 (二) 、在遺傳運算下，遺傳訊息會隨機交換，有利於提供基因演算法多 元化的搜尋方向。

2.3.

基因演算法流程與運算子

基因演算法包含以下步驟：首先，隨機產生初始族群。之後，族群中的

個體由適應函數評估其適應能力，適應能力較佳者配與較高的機率，使之能 在選擇機制下生存，也可以獲得較多的互換(crossover)機會。之後，由選擇 (selection)機制進行淘汰，在此一階段適應能力不佳的個體會被剔除。接著 由基因運算重組基因使組合多樣化，產生新的子代。新的世代產生後，再次 評估(Evaluate)、選擇、基因運算、產生新世代，重複循環直至最佳解達到 穩定的狀態。流程如圖 2.1 所示。 隨機產生一組解的群體 計算每個解的適應度 選擇適應度較高的解 進行交換運算 進行突變運算 程式結束 程式開始 NO NO 滿足解的個數設定數 ? YES 滿足演化世代的數目 ? YES 圖 2.1 基因演算法流程圖

2.4.

編碼與解碼

在基因演算法中，使用二進位碼(如 101001010110…)模擬自然界中生 物的基因排列。 如圖 2.2，基因演算法行演算的過程中，交換、突變等基因運算是在二 進位編碼空間（Binary Space）進行運算。評估與選擇等選擇運算則是在數 字解空間（Solution Space）中進行運算，這兩個空間中的連結的維繫，就

要藉由編碼(encoding)與解碼(decoding)來達成。 圖 2.2 解碼與編碼關係圖 針對編碼的工作，需要先決定以二進位表示一設計參數值所需之位元數 目，假設個體上第 i 個變數 xi位於 [a_i , b_i] 中，且，計算精確度須達到小數 後 n 位，使用 mi 個位元。mi可表示為 : 1 2 10 ) ( 2 i−1 _{< − ×} n _< mi ₋ i i m a b

(2.1) 如此可在 [ai , bi] 上畫出

2

mi 個等分點，而第 c 個等分點所代表的數 值為 : (2.2)

2.5.

基因運算(Genetic Operator)

基因運算是模仿生物在繁衍過程時，進行基因組合重組的過程。基因運 算包含交換與突變。 交換(crossover) 交換過程是由兩個親代進行基因交換產生兩個子代。首先在配對池中選 擇兩個適應程度較高的個體作為親代。再以亂數決定一到數個交換點，將二

decoding

coding space

Genetic

Operations

solution space

Evalution and

selection

encoding

1 2 − − × + = i m i i i i a b p a x

進位碼切割成數段。然後將切割出的片段進行交換，便可以產生新子代。依 照交換點的個數，交換的進行方式可分為下列三種：

a. 單點交換(One Point Crossover)

此一交換方式隨機且均勻地選擇一個小於二進位碼字串長度的正整數作 為交換點。對交換點的右側染色體片段進行交換。

b. 兩點交換(Two Point Crossover)

類似單點交換，隨機且均勻地選擇兩個交換點，並互換交換點中間的片 段。 1 2 3 4 5 6 1 0 0 1 0 1 親代1 親代 2 1 2 3 4 5 6 子代1 子代 2 單點交換 0 1 0 0 1 1 1 0 0 0 1 1 0 1 0 1 0 1 位置 圖 2.3 單點交換(圖上虛線部份為交換點)

1 2 3 4 5 6 1 0 0 1 0 1 親代 1 親代 2 1 2 3 4 5 6 子代 1 子代 2 單點交換 0 1 0 0 1 1 1 0 0 0 1 1 0 1 0 1 0 1 位置 圖 1.3 單點交換(圖上虛線部份為交換點) 兩 兩 4 2 c. 均勻交換(Uniform Crossover) 隨機且均勻地選擇多個交換點，並交換第偶數個或第奇數個片段。圖 2.5 為在偶數片段處交換的結果。 1 2 3 4 5 6 1 0 0 1 0 1 親代 1 親代 2 1 2 3 4 5 6 子代 1 子代 2 單點交換 0 1 0 0 1 1 1 0 0 0 1 1 0 1 0 1 0 1 位置 圖 1.3 單點交換(圖上虛線部份為交換點) 均勻 5 均勻 2

突變(Mutation)

在交換作用的過程中，子代的搜索方向會逐步地靠近，因而有可能掉入 局部最佳解之情況。突變的作用在於彌補交換有可能產生的缺失。突變的作 用機制由隨機決定開啟。在這一階段，基因演算法會在 0 到 1 之間產生一個 亂數，當亂數小於突變率則進行突變。突變發生時會隨機選擇某一個體的二 進位編碼上的某個位元，將值由 1 轉變成 0 或由 0 變成 1。 0 0 0 0 0 0 突變 0 1 0 0 0 0 親代 1 0 0 0 0 0 0 0 1 0 0 0 0 圖 2.6 突變

2.6.

選擇運算(Selection Operator)

選擇運算中，根據「物競天擇、適者生存」。針對個體的適應程度，給 予不同的機率。適應程度越高的個體可分配到越大的機率生存。適應程度低 的個體則越容易被淘汰。在基因演算法中，適應程度越高的個體也越容易被 選入交換池中成為親代，與其他個體產生下一代。 常用的選擇方式有輪盤選擇、競賽選擇等。

a. 輪盤選擇（Roulette Wheel Selection）

轉盤中每一槽的大小，是相對應于每一個個體的適應值大小。在一個世 代中，每一個個體的適應值總和即為輪盤每一槽的面積總和。進行選擇時就 如同丟飛鏢到輪盤上，面積較大的槽就較被容易被射中。若 xκ_{表示第 k 個個} 體，則輪盤選擇的進行方式如下： 1.計算第 k 個個體之適應函數值

₍ k₎ x eval 2.計算總適應值 3.計算個體對應的機率 4.計算個體對應的累積機率 ) ( _ 1

∑

− = = pop size k k x eval sum F sum F x eval p k k _ ) ( =

∑

= k j k p q

5.隨機在[0,1]間產生亂數 r 6.若 qk-1 < r < qk ，選擇第 k 個個體 圖 2.7 輪盤選擇示意圖 b. 競賽選擇（Tournament Selection） 每次在所有個體中任意選出兩個個體，而將適應值較高的個體留下，再 次重複上述程序便可選出另一個適應值較高之個體，藉由兩個被選出之個體 進行交換。 重置(Replacement)與菁英(Elitist)策略 重置是將親代中適應度較高的個體，經由選擇操作直接放入子代中，參 與下一世代的演化。重置的作用在於確保較優良的個體能夠保留下來，有助 於基因演算法在計算過程較於穩定。而菁英策略是直接將親代最佳的個體直 接放入子代，可確定最佳解不會在演算過程遺失。

2.7.

基因參數

基因參數是用來描述基因演算法的系統參數。這些參數的設定會影響到 基因演算法在計算上的表現。藉由對基因參數的瞭解我們將更能夠清楚的瞭 解基因演算法運作的情況。 a. 族群大小(Population Size , N) 族群大小會影響基因演算法的表現與效率。族群小演化壓力較大，演化 的速度較快。但族群小無法提供足夠關於搜尋的空間的資訊有可會使基因演 算法收斂在局部最佳解。大的族群有助於搜尋整個空間，提供較多的資訊。 第k個體所佔的機率

但族群龐大會使收斂速度降低，甚至有可能會面臨無法接受的速度。 b. 交換率(Crossover Rate , Cr) 交換率是交換進行的頻率，在每一代中有 N×Cr 個個體會參與交換。若交 換率越高，新的基因組合會越快出現在族群中。但交換率太高反而會使優良 的個體還沒發揮作用就被取代。交換率過低則會使基因演算法搜尋腳步停滯 不前。 c. 突變率(Mutation Rate,Mu) 突變率是交換突變進行的頻率。假設每個個體的二進位碼字串長度為 L，每一代會發生約 N×L×Mu 次突變。突變有益於增加族群的變動性。然而太 高的突變率會使搜尋變得完全隨機化。一般突變率設為 0.01 至 0.03 之間。

2.8.

結論

基因演算法的特點可歸類為以下幾點： 1. 基因演算法僅對解集合進行運算而不是針對問題的函數運 算，因此不需要提供有關梯度的資訊即可求取最佳解。 2. 基因演算法屬於多路徑隨機的搜尋方式而非決定性的搜尋 方式。

第二部分 依據可行性為基礎的懲罰方式

以基因演算法解有限制條件的最佳化問題時，最常用來處理限制條件的 方法為將限制條件轉換成懲罰函數(penalty function)的方式。假設某個限制 條件的最佳化問題為： minimize

_f(x)

(2.3) subject to n i x x x K k x h J j x g u i i l i k j ,..., 2 , 1 , ,... 2 , 1 , 0 ) ( ,..., 2 , 1 , 0 ) ( = ≤ ≤ = = = ≥

其中 _(x) f 代表欲求最小值的目標函數，g_j(x)代表第 j 個不等式限制條 件， _(x) h_k 是第 k 個等式限制條件，xi則為第 i 個變數，變數的範圍限制在區 間[ l i x ，x_iu]內。透過懲罰函數則可以將目標函數改寫成：

∑

+ = + = J K j j j g x R x f x F 1 2 ) ( ) ( ) (

(2.4) 其中 _(x) F 代表適應函數，其定義為目標函數與限制條件衍生的懲罰項的 總和。在這裡 j R 代表懲罰參數，而 _{g (x)} j 代表的是懲罰函數值，當 ) (x g_j 大 於 0 時傳回 0，當 ) (x gj 小於 0 時傳回函數值的絕對值。另外 k 個等式形的限 制條件可以藉由下式轉換成不等式。 K k x h x gJ+k( )≡δ − k( ) ≥0 =1,2,...,

(2.5) δ 代表一個很小的正整數，例如設為 0.000001。如此不等式條件就會增 加至 J+K 個，如此即可將等式限制條件併入適應函數中。 然而要使用懲罰函數求解帶有限制條件的最佳化問題，會遭遇到的困難 在於如何決定懲罰參數的大小。懲罰參數太小會導致懲罰函數沒有作用，懲 罰參數太大卻也可能使某些有潛力的不適用解在還沒發揮作用時就被淘汰。 要決定懲罰參數的大小，常常需要進行許多次的試誤。在累積大量的計算經 驗後，方能決定較為合適的懲罰參數。在這過程中往往耗費調相當的時間成 本。 Deb 於文獻[10]中提出「依據可行性為基礎的懲罰方式」。這個方法主 要由三個規則組成： (1.) 任何適用(feasible)的解都比所有不適用(infeasilbe)的解好。 (2.) 在兩個適用解中，其目標函數值較好的為佳。 (3.) 在兩個不適用解中，受到懲罰較少的解比較好。 在這個方法下，將不再需要懲罰函參數。第一條規則僅判斷解是否有受 到懲罰。第二條規則僅比較目標函數。而第三條規則僅比較接受懲罰的程度。

如此便不必同時比較目標函數值與懲罰函數值而可以將適用解與不適用解清 楚的區分出來。以可行性為基礎的懲罰方式之適應函數可定義如下： (2.6) 上式中 max f 為適用解中最差的目標函數值，因此所有不適用解的適應函 數值皆較適用解的適應函數值低。      + = ≥ =

∑

= otherwise x g f m j for x g if x f x F m j j j 1 max ( ) ,..., 2 , 1 0 ) ( ) ( ) (

第三章 小板法

3.

前言

吾人採用小板法(panel method)與基因演算法結合，發展二維翼形逆向 設計的工具。在本章中，吾人將解釋小板法基礎理論與數值方法之應用。 小板法乃是以數值方法求解物體表面上的源流或偶流之分佈強度，使擾 動速度勢能滿足邊界條件的方法。小板法的基本觀念為利用有限個線段近似 物體。每一線段稱為一小板（panel），然後在每一小板上分佈未知強度之源 流及偶流，並要求每一小板中心滿足其邊界條件。 在求解源流及偶流強度時，可將積微分式之邊界條方程式離散化，而得 到一組與源流或偶流分佈強度有關之線性聯立方程式。再利用數值方法求解 此聯立方程式組，則可求得源流或偶流分佈強度。物體表面上之速度與壓力 分佈即可依序計算。

3.1.

統御方程式

假設一穩定(steady)均勻流Ur_∞流向一二維翼形，假設流體為非黏性

(inviscid) 、 非 旋 性 (irrotational) 、 密 度 均 勻 (uniform) 且 為 不 可 壓 縮 (incompressible)流體，如圖 3.1 所示，在二維翼形周圍的流場中必存在一 總 速 度 勢 _Φ (velocity potential) 函 數 滿 足 拉 普 拉 斯 方 程 式 (Laplace’s equation)： 0 2_Φ= ∇

(3.1) 其中總速度勢_Φ可表示成入流速度勢(inflow potential)_φ∞與誘導速度勢 (perturbation potential)_φ之和： φ φ φ + = ⋅ + = Φ _∞ Ur_∞ pr

(3.2) 其中_Ur_∞為入流速度向量，_{U∞ =[U∞x U∞y]=[U∞cosα U∞sinα]}，_pr_{( yx, )}為 速度場中任一點之位置向量。又總速度_Vr與總速度勢_Φ的關係式如下： Φ ∇ = Vr

(3.3)

y x

α

∞

U

圖 3.1 二維翼形勢流場示意圖

3.2.

邊界條件

除了上述的關係式之外，尚需滿足下列邊界條件： 1. 無窮遠處條件：無窮遠處的流體不受干擾，流體速度須 趨近於入流速度。 ∞ ∞ → Φ ∇ Ur

(3.4) 2. 物體表面邊界條件（Neumann）：即流體不穿透物體表面。 0 = ⋅ = Φ = ⋅ Φ ∇ V n n nr r r ∂ ∂

_(3.5) nr為物體邊界上向外之單位法向量 3. 庫塔條件(kutta condition)：在水翼的尾緣(trailing edge)處，速度為有限值，此條件之目的在限制本方法之 解為唯一解。 ∞ < e t Vr_.

(3.6) 流場中一點 _pr_{( yx, )}之誘導速度勢_φ，如圖3.2 所示，其邊界值問題可由下 列積分方程式求解： ds q p G n q n q ds n q p G q q p S q q q S ] ( , ) ) ( ) ( [ ) , ( )] ( ) ( [ ) ( r r r r r r r r r

_∫

_∫

∂ ∂ − ∂ ∂ − ∂ ∂ − = − − φ φ φ φ φ

(3.7)

其中 pr

：

為外流場之任一場點 qr：為物體表面之任一點 ) (qr φ ：為外流場之誘導速度勢 ) (qr − φ ：為內流場之誘導速度勢 ) , ( qp G r r ：為二維時之格林函數(Green function) ) , ( ln 2 1 ) , (p q r p q G r r r r π = q p q p r(r,r)= r−r ，表示qr點到pr點之距離 若令µ(qr)=φ(qr)−φ−(qr) q q n q n q q ∂ ∂ − ∂ ∂ = ( ) −( ) ) (r φ r φ r σ 則_µ(qr)可視為_qr點之偶流強度，而_σ(qr)可視為_qr點之源流強度 為滿足物體表面邊界條件，令物體表面法線向量速度為零： 0 = ⋅ = ∂ ∂ + ∂ ∂ = ∂ Φ ∂ _∞ n V n n n r r φ φ

_(3.8) 由於在一物體表面法向分佈偶流時，其內外表面之法向速度為連續，故： n n ∂ ∂ = ∂ ∂φ φ−

_(3.9) 圖 3.2 速度勢與邊界 S 的關係圖 ∞ U − φ S ∞ S nr nr p r qr φ ⋅ ⋅

事實上內流場針對求解之問題並無任何意義，因此可以假設內流場之誘 導速度勢為： ∞ − _=−φ φ

(3.10) 則可將式(3.9)改寫成： n n n ∂ ∂ − = ∂ ∂ = ∂ ∂φ φ− φ∞

_(3.11) 由式(3.11)可得 0 ) ( ₌ ∂ + ∂ − _∞ n φ φ

_(3.12) 將式(3.11)代入式(3.7)，且因滿足式(3.8)之邊界方程式，則可將式(3.7) 改寫成： ds n G q p S q

∫

+ _∂∂ = ( ( ) _∞) ) ( φ φ φ r r

(3.13) 由式(3.13)可知，此時在物體表面只須分佈偶流即可，而其強度為： ) ( ) ( ) ( ) ( ) ( q q q q q r r r r r Φ = + = − =_{φ φ}− _{φ φ∞} µ

_(3.14) 此式即表示物體表面之偶流強度等於為其外表面之總速度勢。最後由式 (3.10)、式(3.13)及式(3.14)可得物體內表面之總速度勢為零之方程式如下： 0 ) ( = ∂ ∂ + = Φ− _∞

∫

ds n G qr µ φ

(3.15) 式(3.15)即為 Dirichlet 邊界條件。因此由上述之推導，若將內流場之誘 導速度勢設為_−φ∞，則可將物體表面邊界條件轉換為 Dirichlet 內邊界條件。 藉由式(3.15)加上其他必要之邊界條件，可求得表面分佈之偶流強度(即物體 外表面之總速度勢)，而由此外表面總速度勢之梯度，即可求出物體表面速度 及壓力等物理量。

3.3.

相關係數

利用上述之邊界條件求得翼形上之偶流強度(即速度勢)，繼而由速度勢 沿翼形上切線方向之梯度求得切向速度 t Vr ，代入伯努力(Bernoulli)方程式以 計算壓力 P，而各點壓力係數 CP為: 2 2 1 ( ) 2 / 1 _{∞ ∞} ∞ _{= −} − ≡ U V U P P C_P t r ρ

(3.16) 又相關之無因次性能係數，升力係數 CL及阻力係數 CD可依下列關係式 計算： s d U U t C C U L C S P L

∫

∞ ∞ ∞ ⋅ − = ≡ ( ) 2 / 1 2 r r ρ (3.17)

∫

∞ ∞ ∞ ⋅ − = ≡ S P D ds U U n C C U D C ( ) 2 / 1 2 r r ρ (3.18)

3.4.

數值方法與離散化

考慮無限水深之情況，一均勻流_Ur_∞流經一可產生升力之二維翼，採用奇 異點分佈(以速度勢為基礎的定值小板法)方法，來求解式(3.2)式之總速度 勢，在翼形表面及跡流上任一點_qr分佈偶流，則流場內一點 _pr之總速度勢可 由下式表示： w q S w B q S _n r ds _n r ds q p U p w w B ) ln ( 2 ) ln ( 2 ) ( ) ( ∂ ∂ + ∂ ∂ + ⋅ = Φ _∞

∫

∫

π µ π µ r r r r _(3.19) 其中 ) (qr µ ：為翼形表面分佈的偶流強度 w µ ：為跡流分佈的偶流強度 ) , ( yx pr ：擬計算誘導速度勢之場點 ) , ( yx qr ：奇異點所在之點 q p r= r−r ：為點qr至點pr之距離 q n ∂ ∂ _：表示 qr點處之翼形表面法向導數

B S 、Sw為翼形表面及尾部跡流 又根據 Dirichlet 邊界條件，物體內部的總速度勢為 0。 0 ) ( = Φ− _{p (3.20)} 所以翼形內部之速度勢可寫成： w q S w B q S _n r ds _n r ds q p U p w w B ) ln ( 2 ) ln ( 2 ) ( 0 ) ( ∂ ∂ + ∂ ∂ + ⋅ = = Φ− _∞

∫

∫

π µ π µ r r r r _(3.21) 以 數 值 方 法 求 解 式 (3.12) ， 將 翼 形 表 面 離 散 為 D N 個 小 板 j B S ， D N j=1~ ，整個跡流視為一個小板S_w，由翼形尾緣延伸至無窮遠處，因此共 有 D N +1 個小板分佈偶流，每個小板上偶流均為定值，因此式(3.21)可離散 化如下式所示： w i w N j j i j i i D p U p _, 12 , 2 0 ) ( α π µ α π µ ₊ + ⋅ = = Φ

∑

= ∞ − r r r _(3.22) 其中

∫

_∂∂ = j B j j S B q j i r ds n ln ) ( , α w S q w i r ds n w w

∫

_∂∂ = ( ln ) , α 另外，庫塔條件式(3.6)式可表示為： 1 1 µ µ µ_w =Φ_ND −Φ = _ND − (3.23) 求解 ND +1 個未知數，需 ND +1 個方程式。可由式（3.22）及（3.23）解之。 解得翼形表面分佈的偶流強度 _D N µ µ µ₁, ₂,... 及 w µ 之後，再進一步求解翼 形表面的切線速度。 假設速度勢沿著翼形表面的小板呈二次分佈，則可表示成:

(3.24) (3.25) 由上式第 i 個小板的梯度為 又第 i-1 個及第 i+1 個小板速度勢可表為: (3.26) (3.27) 其中li為第i個小板的長度。 令 (3.28) (3.29) 由式(3.26) 、式(3.27) 、式(3.28)及式(3.29)可求得: (3.30) 在第 1 個小板及最後一個小板_Vt則分別可簡化成 DB 及 DA。如此切線速 度即可求得。再由壓力係數定義即可求得相關係數。 t V t cs b ds d = ⋅ Φ = + = Φ ₂ b ds d s i = Φ =0 2 ) (s =Φ_i+bs+cs Φ 2 / ) ( 1 2 1 i a a a i i i = Φ +bl + cl l = − l +l Φ _{− −} 2 / ) ( ₁ 2 1 i b b b i i i =Φ +bl +cl l = l +l Φ+ + a i i l DA=(Φ₊₁−Φ )/ b i i l DB=(Φ₋₁−Φ )/ t a b b a V l l l DB l DA b = − ⋅ − ⋅ =

第四章 二階段設計法

4.

前言

在二維翼形設計問題中，我們經常遭遇到兩類設計情況，分別為直接設 計問題(Direct problem) 和逆向設計問題(Inverse problem)。依據文獻[1]， 直接設計問題指的是由給定的二維翼形幾何形狀計算其流場的特徵如流速及 壓力等。直接設計法通常採用最佳化演算法與流體動力解析程式結合。演算 過程中疊代二維翼形的幾何來對給定的目標函數進行最佳化（例如求最小的 阻力係數）。在逆向設計問題中採用設計變數當作描述壓力分佈的參數，求 得滿足某些要求條件下(例如最小負壓力係數必須大於水之蒸氣壓)的壓力分 佈，而非直接描述二維翼形的幾何形狀。本研究採取二階段逆向設計的架構， 對已知的二維翼形進行局部的修改。第一階段中，將對於設計變數如升力係 數、空化係數等之要求，利用限制條件來表示，再透過最佳化的方法，求得 滿足限制條件的壓力分佈。第二階段的工作則負責將第一階段得到的壓力分 佈，作為目標壓力分佈，以流體動力解析程式與基因演算法結合，求得對應 壓力分佈的二維翼形幾何形狀。

4.1.

壓力分佈設計

流體機械如螺槳、飛機翼等在流體中工作時，均利用其二維翼形之上下 表面壓力分佈的不同所產生的壓力差，來提供此流體機械所需的推力、吸力 或升力。壓力分佈曲線的形狀、包含的面積、壓力係數的最小值都會造成流 體機械性能上的差異。正因為有此一特性，吾人可藉由設計壓力分佈，使對 應的二維翼形更加符合性能上的要求。

4.1.1.

壓力分佈曲線的建立

在以基因演算法進行壓力分佈設計前，首先需建立壓力分部曲線的表示 方式。從設計的角度上來看，吾人希望曲線的表示方法能使用較少控制點而 準確地表示。同時也希望在使用上有較高的自由度，可以方便地對壓力分佈 曲線形狀進行局部的調整。基於上述考量，在本研究中壓力分佈採用 B 木條 曲線表示。 B 木條曲線的定義與特性

一條由 n+1 個控制點組成的 k 階 B 木條曲線(B-spline curve)，其定義如 下: 1 2 , ) ( ) ( 1 _{min max} 1 , + ≤ ≤ ≤ ≤ =

∑

+ = n k t t t t N D t P n i k i i

(4.1) 上式中_P(t)為曲線上任一點的位置，是參數_t的函數。Di為曲線第 i 個 控制點的位置。 _{( )} , t

N_ik 是階數(order)k 的 B 木條曲線基底函數（B-spline Basis

function）。定義為：

(4.2a)

(4.2b)

i

v 是 knot vector 的元素。在本文中使用 Open Uniform knot vector，其

計算規則如下： (4.3) 式(4.1)中 B 木條曲線只有在一部分範圍的參數值下才不是零。在範圍外 的 B 木條曲線基底函數皆為零。B 木條曲線基底函數的這一項特性使得 B 木 條曲線的控制點僅能影響一定範圍內的曲線。所以吾人可以藉由調整部分控 制點，改變某一部分曲線而不影響其他部分的曲線段。B 木條曲線還可以經 由改變階數(order)調整曲線的平順程度。當階數越高，曲線會越平滑。反之， 階數數越低，曲線則越接近控制點的連線。 除了以上兩點，B 木條曲線的起點與終點分別是第一點控制點及最後一 點控制點。而端點的斜率也分別由控制點的前兩點 1 B 、B₂與後兩點B_n、B_n+1 決定。 1 1 , 1 1 1 , , ) ( ) ( ) ( ) ( ) ( + + − + + − + − − − + − − = i k i k i k i i k i k i i k i v v t N t v v v t N v t t N    ≤ < = + otherwise v t v if t Ni i i 0 1 ) ( 1 1 , 1 2 2 1 1 1 0 + + ≤ ≤ + + − = + ≤ ≤ + − = ≤ ≤ = k n i n k n v n i k k i v k i v i i i

以 B 木條曲線描述壓力分佈曲線 利用 B 木條曲線，表示壓力分佈，其中將翼形上下表面之壓力分佈分開 (以導緣及尾緣為分界點)，並各以一條 B 木條曲線表示，如圖 4.1 所示，上 方的曲線對應翼形上表面的壓力分佈，下方的曲線則對應下表面的壓力分 佈。兩條曲線由左到右分別由 11 個控制點組成，在假設無攻角狀況下，第一 點的控制點座標固定在(0,1)，以及最後一點控制點的 x 座標固定為 1。同時 上下兩條曲線最後一點設為同一點，使之滿足水流流離尾緣時上下翼面的速 度相同(亦即壓力相同)。改變其餘控制點的 xy 座標值，來對整個壓力分佈曲 線之形狀進行調整

。

0 0.2 0.4 0.6 0.8 1 x/C 0.8 0.4 0 -0.4 -0.8 Cp

Pressure distribution curve Control points 圖 4.1 壓力分佈建構圖

4.1.2.

以基因演算法進行壓力分佈設計

本階段的主要工作，在於修改已知二維翼形的壓力分佈使翼形能夠更符 合工作環境的限制與需求。因此吾人參考文獻[8]，設定壓力分佈設計問題如 下 minimize : f(Α)=CL−CL−specified (4.4)

subject to: Cp.u及 Cp.l分別表示翼形上表面及下表面的壓力係數，而 Cp*表示負的空 化係數。S 表示翼形上表面尾緣壓力分佈曲線斜率最大值。Ninflection-u 及 Ninflection-l分別代表上下表面壓力分佈曲線上反曲點的個數。另外 A 為壓力分 佈曲線控制點，下標 1 至 11 為翼型下表面壓力分佈曲線控制點，下標 12 至 22 為上表面壓力分佈曲線控制點。 ] [Α1 Α2 Α3 Α22 = ... Α 其中向量Αi包含個控制點座標值，即       = Α i i i y x _。 此最佳化問題定義為要求升力係數達到給定值，並且符合限制條件。升 力係數可由下式求得： (4.5) 其中限制條件 1 g 為確保上下兩條壓力係數曲線包圍的面積為正，也就是 升力方向向上。限制條件 2 g 為限制上方的壓力曲線最小值不會小於負的空化 係數，此一限制乃考慮翼形上表面發生空化現象的避免。限制條件 3 g 為參考 文獻[8]所得，限制尾緣壓力梯度的變化，當尾緣壓力梯度變化較小時，可以 減少流離的發生。然而文獻[8]為設計飛機翼斷面，可能與水翼上有所不同， 所以吾人將尾緣壓力分佈曲線斜率最大值以 S 表示。另外加上限制條件 4 g 、 5 g 控制壓力分佈曲線的平滑度。限制條件g₄限制二維翼形上表面的壓力分佈 9 . 0 1 . 0 0 2 ) ( . 5 9 . 0 1 . 0 0 0 ) ( . 4 0 ) ( . 3 0 ) ( . 2 1 0 0 ) ( . 1 inf 5 inf 4 1 . 3 * 2 . . 1 < < ≥ − = < < ≥ − = ≥ − = ≥ − = < < > − = − − → C x in N g C x in N g dx dC S g C C g c x in C C g l lection u lection c x u p p peak suction p u p l p Α Α Α Α Α

∫

− = = ∞ 1 0 . . 2 ) ( ) ( 2 1 C x d C C C U L C_{L pl pu} ρ

曲線中段不可以有反曲點發生。限制條件_g5限制二維翼形下表面的壓力分佈 曲線反曲點不可多於兩個以上。限制條件 4 g 、g₅為吾人觀察 NACA 系列等翼 形所得之結論。 在最佳化過程，給定初始壓力曲線幾何控制點 A0 。壓力曲線幾何控制點 A 的搜尋範圍如下式所示： ∆Χ Α Α= 0 +

_(4.6) 其中∆Χ代表搜尋範圍。 上述的壓力分佈表示方式與基因演算法結合，每一個壓力分佈扣除固定 的控制點座標則對應到 38 個變數。我們將變數變化的範圍設定為-0.05 xi0 <

△xi < 0.05xi0及 -0.1(yi0-1) <△yi < 0.1(yi0-1)。

其中 y 方向座標值 yi0-1，表示 Cp值以 1 為基準值，當沒有速度擾動時， 壓力係數應為 1。 以基因演算法進行壓力分佈最佳化步驟如下: 步驟一、產生初始族群。 步驟二、計算族群中個體的的適應值。首先將個體由二進位碼解碼還原成 B 木條曲線控制點組。接著呼叫 B 木條曲線副程式將控制點繪成壓力 分佈曲線。由式（4.5）計算壓力分佈曲線所包圍的面積可得到升力 係數。然後檢查每一個體是否違反限制條件。若無違反限制條件，則 升力係數數值即為適應含數值。否則以「依據可行性為基礎的懲罰方 式」對違反限制條件的個體計算其適應含數值。       ∆ ∆ ∆ ∆ ∆ ∆ = 22 22 2 1 2 1 ... ... y x y y x x ∆X

步驟三、進行菁英策略。將具有最高適應含數值的個體無條件地複製到新 一代的族群中。 步驟四、應用基因操作產生下一代的個體。依據每一個體的適應含數值， 以競賽選擇方式選擇一對親代個體，經由交換和突變產生一對子代個 體。 步驟五、若所需的下一代個體數未達族群大小，回到步驟四。否則將此新 的一代當做下次演化的世代。 步驟六、若未達最大世代數，回到步驟四。 步驟七、將最後一代中，具由最高適應函數值的個體當作最佳解輸出。 上述全部步驟之流程圖如圖 4.2 所示。參數設定如表 4.1 所示。 圖 4.2 壓力分佈最佳化流程圖 隨機產生一組解的群體 計算每個解的適應度 選擇適應度較高的解 進行交換運算 進行突變運算 程式結束 程式開始 NO NO 滿足解的個數設定數 ? YES 滿足演化世代的數目 ? YES

4.2.

二維翼形的逆向設計

當壓力分佈設計階段完成後，接下來的工作就是由已知的壓力分佈求得 對應的二維翼形。在這一階段，本研究採用勢流計算程式計算程式結合基因 演算法，進行二維翼形的逆向設計的工作。 在本研究中，二維翼形幾何的表示亦採用上下表面兩條 B 木條曲線（型 態為 open uniform B 木條曲線）來表示翼形，如圖 4.3。在每一條曲線上分 佈 11 個控制點。控制點的第一點及最後一點分別固定在座標（0,0）及(1,0) 上以確定翼形是封閉的形狀。並且將控制點第二點 x 座標固定為 0，確保在 導緣在兩條曲線連接處斜率連續。變化控制點第二點的 y 軸座標，及其餘控 制點的 xy 座標值，來對整個翼形進行形狀變化。 運用基因演算法在給定一個壓力分佈之下計算二維翼形的幾何形狀。最 佳化問題之目標函數如下式所示： minimize

(4.7) 其中， _{( ,( ) )} i P C x C Β 為由 B 木條曲線所表示之翼形所算出的壓力分佈， ) ) (( _i Pt C x C 為目標翼形的壓力分佈。而矩陣 B 包含所有二維翼形的幾何控制 點，其中向量Βi則表示個控制點之倒置向量，下標 1 至 11 為翼型下表面控 制點，下標 12 至 22 為上表面控制點。 在最佳化過程，首先給定一目標壓力分佈 Cpt及初始水翼幾何控制點 B0。 二維翼形翼幾何控制點 B 的搜尋範圍如下式所示： ∆Χ Β Β= 0+

_(4.8) 其中_∆Χ代表搜尋範圍，其數值介於 [-0.04C, 0.04C] 之間，而 C 代表 翼展弦長(Chord length)。搜尋範圍為文獻[3]所建議的大小，並且搜尋範圍 大小亦超過翼形最大厚度的一半，因此吾人亦採用相同的搜尋範圍。

∑

_ − _ = 2 ) ) (( ) ) ( , ( i Pt i P C x C C x C f Β ] ... [Β1 Β2 Β3 Β22 = Β       ∆ ∆ ∆ ∆ ∆ ∆ = 11 11 2 1 2 1 ... ... y x y y x x ∆X

0 0.2 0.4 0.6 0.8 1 x/C -0.04 0 0.04 0.08 0.12 y/C Foil section Control points 圖 4.3 二維翼形建構圖 將所有控制點的位置向量 B 代入 B 木條曲線程式中產生二維翼形。再將 產生的二維翼形連結流體動力分析程式，產生每組水翼幾何對應的壓力分 佈，經由基因演算法將翼形幾何的控制點逐步修正使二維翼壓力分佈能與目 標壓力分佈一致。處理過程之流程如圖 4.4 所示。參數設定如表 4.2 所示。 圖 4.4 二維翼形反向設計流程圖 隨機產生一組解的群體 計算每個解的適應度 選擇適應度較高的解 進行交換運算 進行突變運算 程式結束 程式開始 NO NO 滿足解的個數設定數 ? YES 滿足演化世代的數目 ? YES

參數說明 參數值 編碼位元 10 族群大小 100 最大世代數 1000 交換率 0.8 突變率 0.02 選擇方式 競爭選擇 交換方式 均勻交換 表 4.1 壓力分佈設計之基因參數 參數說明 參數值 編碼位元 14 族群大小 100 最大世代數 4000 交換率 0.8 突變率 0.02 選擇方式 輪盤選擇 交換方式 均勻交換 表 4.2 翼形逆向設計之基因參數

第五章 結果與討論

5.

前言

本章將就第四章中所提的壓力分佈設計與二維翼形逆向設計之結果作討 論，並再最後結合兩者進行一個設計例子。

5.1.

二維翼形逆向設計之驗證

在本節中，將以 NACA4412 及 NACA2410 作為方法驗證之對象，並討 論以基因演算法進行逆向設計之限制。 首先，圖(5.1)、(5.2)、(5.3)、(5.4)所示為 NACA4412 為逆向設計目 標的驗證。其中(5.1)圖所示為 NACA4412 原有之幾何形狀(以虛線表示)， 以及用於逆向設計之初始翼形(以實線表示)，後者是以 NACA4412 原有形狀 為樣本，以手動方式調整而決定 B 木條曲線控制點。 0 0.2 0.4 0.6 0.8 1 x/C -0.04 0 0.04 0.08 0.12 y/C

Target foil section Initial foil section

圖 5.1 NACA4412 與初始翼形之比較

圖(5.2)為 NACA4412 之壓力分佈與手動取得之翼形的壓力分佈的比 較。可以看出兩者間有相當的差距。圖(5.3)為藉由基因演算法調整初始手動 取得控制點後之翼形與 NACA4412 形狀之比較。

0 0.2 0.4 0.6 0.8 1 x/C 0.8 0.4 0 -0.4 -0.8 -1.2 Cp

Target pressure distribution Initial pressure distribution

圖 5.2 NACA4412 與初始翼形壓力分佈之比較 0 0.2 0.4 0.6 0.8 1 x/C -0.04 0 0.04 0.08 0.12 y/C

Target foil section Cal. foil section

圖 5.3 NACA4412 與逆向設計結果之比較

圖(5.4)為基因演算法調整後翼形之壓力分佈與 NACA4412 壓力分佈之 比較，可以看出在導緣之部分已有明顯之改善，但在尾緣仍有差異存在。演 化過程如圖(5.9)所示。

0 0.2 0.4 0.6 0.8 1 x/C 0.8 0.4 0 -0.4 -0.8 Cp

Target pressure distribution Cal. pressure distribution

圖 5.4 NACA4412 與逆向設計結果壓力分佈之比較 接著，圖(5.5)及圖(5.6)為以上述基因演算法方式取得之 NACA4412 翼 形 為 初 始 翼 形 逆 向 設 計 取 得 NACA2410 之 計 算 結 果 。 圖 (5.5) 為 NACA2410 與初始翼形(及 NACA4412)及逆向設計結果之比較。圖(5.6)則 為上述三者之壓力分佈比較。可以明顯看出控制點調整的結果。演化過程如 圖(5.10)所示。 0 0.2 0.4 0.6 0.8 1 x/C -0.04 0 0.04 0.08 0.12 y/C

Target foil section Initial foil section Cal. foil section

0 0.2 0.4 0.6 0.8 1 x/C 0.8 0.4 0 -0.4 -0.8 Cp

Target pressure distribution Initial pressure distribution Cal. pressure distribution

圖 5.6 NACA2410 與逆向設計結果壓力分佈之比較

5.2.

二階段設計之應用例

在本節將展示以二階段設計法，嘗試針對 NACA4412 壓力分佈，針對設 計條件進行搜尋之實例。 本例中，假定之升力係數與 NACA4412 相同為 0.48，而空化係數為 0.7。 尾緣壓力分部曲線斜率(式(4.4)中 g3式子裡之 S)限制與文獻[8]相同設為 2.3。初始之壓力分佈為採用手動取得 NACA4412 壓力分佈曲線。 在圖(5.7)可以大致判斷滿足設計條件之壓力分佈面積與 NACA4412 相 同，而翼形之表面壓力係數最小值則滿足空化係數大於-0.7 的要求。演化過 程如圖(5.11)所示。 圖(5.8)展示以 NACA4412 為初始翼形逆向求得滿足設計條件之翼形。 演化過程如圖(5.12)所示

0 0.2 0.4 0.6 0.8 1 x/C 0.8 0.4 0 -0.4 -0.8 Cp

Cal. pressure distribution Target pressure distribution Initial pressure distribution

圖 5.7 NACA4412 壓力分佈與滿足設計條件之壓力分佈之比較 0 0.2 0.4 0.6 0.8 1 x/C -0.04 0 0.04 0.08 0.12

y/C _{Cal. foil section} Intial foil section

0 1000 2000 3000 4000 Generation 0 10 20 30 40 50 F itnes s Best Average 圖 5.9 逆向求取 NACA4412 翼形之適應值對世代數圖 0 1000 2000 3000 4000 Generation 0 10 20 30 40 50 F itnes s Best Average 圖 5.10 逆向求取 NACA2410 翼形之適應值對世代數圖

0 200 400 600 800 1000 Generation 10 20 30 40 50 F itnes s Best Average 圖 5.11 壓力分佈設計之適應值對世代數圖 0 1000 2000 3000 4000 Generation 0 10 20 30 40 50 F itnes s Best Average 圖 5.12 逆向求取滿足設計條件翼形之適應值對世代數圖

第六章 結論

本研究採用二階段設計法，第一階段採用基因演算法針對現有的二維翼 形進行壓力分佈的改變，使壓力分佈能滿足性能上的要求。第二階段藉由基 因演算法結合小板法將上一階段得到的壓力分佈反向求取對應之翼形。 逆向設計方法是一種強而有力的設計工具。用這些工具可以求出給定壓 力分佈對應的二維翼形，但是在使用這一些設計方法時，使用者必須先將他 的設計目標適當地轉換成壓力分佈，並且同時滿足流體動力性能上的要求。 這對於二維翼形的設計者造成某種程度的不方便。在本研究中展示以最佳化 的方式，採用基因演算法有效地將設計目標轉換成滿足所設定限制條件之壓 力分佈，如此可提供設計者某一程度的便利性。 然而在研究進行中累積的經驗告訴吾人，二維翼形逆向設計及最佳設計 分成二階段進行仍有一些的問題待釐清。其中問題在於壓力分佈雖然可以提 供吾人關於翼形升力的特性，但是卻無法提供阻力及相關的資訊。這個問題 在翼形進行最佳化時，形成評估的盲點。因為能夠滿足如式(4.4)及其限制條 件所表示之最佳化問題的壓力分佈曲線並非唯一，然而對於這些滿足求解條 件的壓力曲線並沒有可以作為評估的條件。為了彌補這個缺點，吾人擬參考 文獻[8]中採用一體求解的方式，即在整個最佳化問題中，加入了阻力最小的 觀點，作為最佳化問題的依據，如此，修改後的最佳化問題定義可寫成: minimize : f(Α)=Cd (6.1) subject to: 其中前述反向設計之程式，亦須整合於求解過程中，以便由壓力分佈求 9 . 0 1 . 0 0 2 ) ( . 6 9 . 0 1 . 0 0 0 ) ( . 5 0 ) ( . 4 0 ) ( . 3 1 0 0 ) ( . 2 ) ( . 1 inf 6 inf 5 1 / . 4 * 3 . . 2 1 < < ≥ − = < < ≥ − = ≥ − = ≥ − = < < > − = − = − − → − c x in N g c x in N g dx dC S g C C g x in C C g C C g l lection u lection c x u p p peak suction p u p l p specified L L Α Α Α Α Α Α

取對應之翼形及所屬之阻力。

阻 力 的 預 測 可 採 用 Squire-Young relation[11] 計 算 。 Squire-Young relation 是由邊界層理論與勢流速度推導而得的經驗公式。由已知的壓力分 佈，吾人可以使用 Squire-Young relation 計算阻力係數，進而以阻力係數作 為評估的標準。 吾人重新構思的最佳化過程如圖(6.1)，如此，可使整個二維翼形最佳化 過程更加完整。 圖 6.1 二維翼形最佳化流程圖 隨機產生一組解的群體 計算每個解的適應度 選擇適應度較高的解 進行交換運算 進行突變運算 程式結束 程式開始 NO NO 滿足解的個數設定數 ? YES 滿足演化世代的數目 ? YES 翼形逆向設計 是否滿足給定升力係數 ? YES Cd=Cd max NO

第二部分

第一章 研究目的及文獻回顧

因為速度或攻角的改變，導致水翼或螺槳某些部分的壓力低於水的飽和 蒸汽壓時，就會發生空化現象。一般來說，空化問題是在高速作動下之流體 機械最希望避免的不利現象，它不僅造成機械效率衰減，同時其所伴隨之振 動、噪音與結構侵蝕問題，亦對於流體機械的可靠度產生負面的影響，因此 預測空化發生的條件，便成為水翼或螺槳葉片設計必須考量的因素。本研究 的目的在於，雖然可以實驗方式分析空化現象，但基於經濟與時間的考量， 若能使用計算方式進行研究，在工程上精度要求的條件下，不僅可以大幅縮 短產品開發時程，同時可以進行產品性能之最佳化，提高產品的市場價值； 另一方面，藉著有系統的研究來與世界之發展接軌，以提昇國內在這個領域 的成果與知識。 空化現象在造船工程上的重要性是眾所皆知，特別是在提高船速所使用 的高速推進系統上如：噴水推進器、超空化螺槳、穿水式螺槳等皆面臨了空 化問題的挑戰。同時該現象的複雜度及困難度高，加上空化流場屬於紊流流 場，更加阻礙吾人對於其物理現象的了解。目前為止，大部分的研究仍然以 實驗方式進行，但近年來，由於計算流體力學的進步與發展，已有學者使用 計算方式研究此類問題。早期的研究是將重點放在勢流理論的計算，如： Kinnas and Fine [12]， Kinnas[13]， Kinnas [14]， Fine and Kinnas[15]， Kinnas and Mazel [16]，至於國內方面，則有陳建宏教授[17]，他分別處理二維超空 蝕、三維局部空蝕水翼數值計算以及二維超空蝕水翼之空蝕係數與翼形幾何 參數關係之回歸分析，台大郭真祥教授發展完成處理二維水翼數值計算之高 階小板法[18]，亦完成了二維局部空蝕及超空蝕水翼的計算[19]。但是隨著數 值方法以及模型的進步，現在的研究已經轉移到黏性流場的模擬計算，以觀 察空化區現象為出發點的計算方法大致可以分成兩類，其一是界面追蹤法 （Interface tracking scheme），此法是將空化區當成一固體邊界，整個計算空 間只考慮在液體的情況下，如：Chen and Heister [20]、Deshpande et al. [21]； 其二是二相流法（Two phase flow scheme），此法是不預設空化區的大小， 整個計算網格考慮液體與氣體都存在的情況，如：Delannoy and Kuney [22]、

Chen and Heister [23]，到目前為止，假使能夠使用一較好的空化模型，較能 模擬片狀空化流場。

本文是以Reynolds-Averaged Navier-Stokes equation (RANS) 為統御方程

式（governing equation），利用數值方法加以離散，同時將計算空間離散成 網格點，再配合 κ - ε 紊流模型來計算雷諾應力，以空化模型（cavitation model）加上速度、密度與壓力間的偶合關係來計算二維翼形之空化流場，並 進一步解析相關的流體動力特性及預測空化初始是否發生，以驗證數值計算 結 果 之 適 用 性 ； 另 外 針 對 二 維 翼 形 之 四 個 不 同 的 幾 何 類 型 ， 分 別 為 NACA0012、NACA4412、NACA66（mod）、Eppler，以及數值參數的變化， 探討空化係數的影響和無空化情況下的動力特性。

第二章 理論基礎

2.1. 統御方程式

2.1.1. Reynolds-Averaged Navier-Stokes 方程式

本 文 是 以 計 算 空 化 流 場 為 出 發 點 ， 在 此 情 況 下 ， 假 設 一 非 穩 定 （unsteady）、黏性（viscid）、可壓縮（compressible）流體，以攻角α 流向 一二維翼。雖然Navier-Stokes 方程式可以計算層流以及紊流的流場，然而卻 不適合直接去計算紊流流場，起因於無法計算在高雷諾數下，速度及壓力隨 著時間改變的微小擾動，因此為了計算黏性空化流場，便將Navier-Stokes 方 程式修正為Reynolds-Averaged Navier-Stokes（RANS）方程式如下：

( ) ( )

u uu P

[ ] [

u u

]

t m +∇⋅ m =−∇ +∇⋅ m∇ − ∇ m∇⋅ ∂ ∂ _{ρ ρ µ µ} 3 2 _（2.1） 其中等號左邊第一項為隨時間改變的動量變化，第二項為非線性且對流 場影響甚大的對流傳導項（convection term），等號右邊第一項為壓力梯度變 化項，第二項為流體黏性效應所造成的擴散項（diffusion term），最後一項 則是為了計算空化效應所加入的壓縮項（compressibility effect）。

2.1.2. 連續方程式及體積比率方程式

對於空化流場可壓縮流體而言，由於密度會隨著不同時間而改變，所以 連續方程式（continuity equation）必須修正為

( )

=0 ⋅ ∇ + ∂ ∂ u t m m ρ ρ _（2.2） 式（2.2）中的混合密度ρm定義如下： ) 1 ( l v l l m ρ α ρ α ρ = + − （2.3） 而在發生空化時，每一個計算網格流體所佔的體積就會因此改變，體積比率 方程式就是計算在同一個網格內，流體和氣體所佔的比率：

) ( −₊ + =       ⋅ ∇ + ∂ ∂ m m u t l l & & α α _（2.4） 其中 − m& 是汽化率， m&+ 表示凝結率，α l 即是液體所佔的體積比率。

2.2. 空化模型

空化區內密度會隨著時間而改變，因此空化模型便需要引用混合密度的 觀念，而最常用的方法是以傳輸方程模型（transport equation-based model） （TEM）為主，如：Singhal et al.[24]，Merkle et al.[25]，Kunz et al.[26]，Senocak and Shyy[27]，上述這些不同的空化模型，主要的區別是在於源流項（source term），亦即 −

m& 和 m&+ 的不同，然而在質量守恆的條件下，不同的模型都

有不同的經驗參數去調整，本文是使用Senocak and Shyy[27]來考慮空化區的 現象，如下：

[

]

∞ ∞ − ₌ − t U p p MIN C m l l v l v dest ) 2 1 ( , 0 2 ρ ρ α ρ & （2.5） ∞ + ₌ − t C m l l l v prod ρ α α ρ 2(1 ) & （2.6） 除此之外，空化模型忽略了表面張力和浮力效應所造成的影響，而在式 （2.5）、（2.6）中的 Cdest和 Cprod分別是經驗常數，控制著空泡的崩潰和生 成。

2.3. 紊流模型

紊流模型最主要的目的就是求解推導RANS 過程中所增加的雷諾應力，

一般是使用Reynolds stress model（RSM），然而 RSM 在計算上增加了六項 雷諾應力，需要用六個方程式去解，無形中增加了許多的計算時間，同時所 要界定的係數就相對增加，但是此模式在某些較為複雜流場處預測紊流的效 果不錯，因此仍有其應用價值。另一方面，便有了代數方程求解紊流參數的 方法出現，其中最有名的是包氏渦流黏性模式（Boussinesq Eddy Viscosity Model）[28]，本文使用此模式來計算雷諾應力，假設關係如下：

k u _ij t ij µ δ τ 3 2 + ∇ − = r （2.7） ε ρ µ mCµk2 t= （2.8） 其中µt為紊流黏度（Turbulence viscosity），C_µ為一經驗常數，約為0.09，δ_ij 在 i = j 為 1 及 i≠j 為 0，k為紊流動能（Turbulence Kinetic Energy），ε為紊 流消散率（Turbulence Dissipation），本文利用雙方程式模式中之 k-ε 模式 （k-ε model）其中依照 Launder and Spalding[29]之建議來計算 k 和ε，如下式：

(

)

ε σ µ µ ρ ρ + −         ∇       + ⋅ ∇ =       ⋅ ∇ + ∂ ∂ k k t m mk k u k P t （2.9）

(

)

k c k P c u t k t m m 2 2 1 ε ε ε σ µ µ ε ρ ε ρ ε ε ε − +         ∇       + ⋅ ∇ =       ⋅ ∇ + ∂ ∂ （2.10） 其中σ_k、σ_ε、c_ε1、c_ε2經由實驗和計算得來，如下表： MODEL Cµ Cε1 Cε2 σk σε k-ε 0.09 1.44 1.92 1.0 1.3 表 2-1 k-ε 模式的經驗常數

P_k為紊流動能製造項（Production of Turbulent Kinetic Energy），定義為：

j i ij k x U P ∂ ∂ − = τ （2.11） 由於流體在壁面為不滑移之狀態（no-slip condition），所以靠近壁面流場的 速度梯度變化很大，此處局部雷諾數很小，而紊流模式是以高雷諾數的流體 為主，因此對於靠近壁面的流場需要做額外的處理，在下一節會另有說明。

2.4. 邊界條件