國立臺灣大學工學院土木工程學系 碩士論文
Department of Civil Engineering College of Engineering
National Taiwan University Master Thesis
應用形狀最佳化與拓樸最佳化於薄殼結構設計 Application of Shape Optimization and Topology
Optimization in Thin Shell Structure Design
黃懷寬
Huai-Kuan Huang
指導教授: 呂良正 教授 Advisor: Prof. Liang-Jenq Leu
中華民國 111 年 7 月
July, 2022
致謝
本論文能夠順利完成，首先要感謝的是恩師呂良正教授，最初在大四升研究 所要找指導教授時，看到老師展示碩士生研究內容時，便讓我萌生對最佳化的興 趣，進了研究所以後，碩一聽老師講解高等結構學時，也能夠感受到老師對教學 的用心，碩二則從每週 GM 中得到老師許多回饋，訓練我思考解決問題的方法，
非常感謝老師這兩年來照顧與教導，讓我從中學到許多待人處事的態度，學生在 此向老師致上最深的敬意。
論文口試期間，承蒙黃仲偉教授、宋裕祺教授以及郭世榮教授的詳加指教與 建議，讓本論文的內容可以更為完備，在此感謝各位教授的協助，也特別感謝黃 仲偉教授於百忙之中校閱本論文內容，指正出錯誤或是不足之處，讓我可以及時 修正。
感謝二博這兩年來的幫忙，無論是討論並提供我關於研究方向上的建議，還 是協助處理日常瑣碎的行政事務與經費報帳等，二博都給了我很大的幫助，謝謝 二博這些年來對於呂門的貢獻，祝二博在日本研究順利。感謝 R08 的學長姊鈺 翔、穎君以及啟聰，謝謝你們提供家產幫助我解答課業上的難題，也感謝鈺翔學 長在我研究初期協助我了解最佳化研究的程式，讓我的研究可以順利進行。感謝 與我同屆的俊廷以及牧軒，這兩年來在課業及研究上一同努力，彼此互相幫助也 互相鼓勵，給予我許多學習及研究的動力，感謝俊廷時常與我討論研究方向及內 容，幫助我做出許多研究上的決定；感謝牧軒經常在我寫程式遇到困難時提供我 解決辦法，幫助我克服許多研究上的難關。感謝學弟妹智傑、權恩、采霏、鈺庭 以及家成在平時以及口試期間的幫忙，謝謝智傑協助我論文校稿，幫我找出論文 中的錯字以及讓我的語句更為通順；謝謝權恩平時 GM 的提問，讓我從回答中反 思自己的研究內容；謝謝采霏在口試預演時幫我記錄老師的提問，並提供我簡報 上的建議；謝謝鈺庭於口試期間擔任紀錄，幫助我記錄各個口委的提問；謝謝家 成與我一同討論拓樸，給予我研究上不一樣的想法。這兩年與大家一起在研究室
共度的時光都會成為我寶貴的回憶，非常感謝大家。
黃懷寬謹致於 國立台灣大學土木研究大樓 511 室 2022 年 8 月
摘要
現今建築經常以美學為主要設計依據，其中便以薄殼結構結合流線型以及創 造開放空間感之設計較為常見，然而這些基於美觀之設計卻不一定是最有效率之 設計，因此本研究將結合建構自由曲面技術與最佳化理論，設計出兼具美觀與良 好力學行為之結構。
本研究採用 NURBS 做為自由曲面建構方法，透過控制點座標、節點向量以 及基底函數建立曲面形狀，再結合最佳化理論將控制點座標做為設計變數，以對 結構形狀進行調整，建構出在設定目標下具有最好表現之薄殼結構幾何形狀，並 進一步探討延伸應用於支承最佳化的方法。
除了形狀最佳化以外，本研究也會應用到結構最佳化中的拓樸最佳化，藉由 演算法迭代找出最有效率之材料分布位置，而本研究中所使用的拓樸最佳化演算 法為 BESO，將會以不同例題演示現有之 BESO 演算法的問題與不足之處，並針 對這些問題加以修正。
本研究會將經過形狀最佳化後之薄殼結構，應用拓樸最佳化方法將結構中無 效率處之材料移除，並將有效率處之材料保留，藉此將結構中之材料分配至最有 效率傳遞力量之位置，使整體結構在相同材料使用量的情況之下，勁度表現能夠 進一步提升。
本研究透過形狀最佳化與拓樸最佳化方法，應用於不同類型之薄殼結構，使 最終設計結果能夠具有良好之結構勁度表現，藉此作為未來工程師設計薄殼結構 時之參考依據。
關鍵字：自由曲面、結構最佳化、形狀最佳化、拓樸最佳化、薄殼結構、殼橋
Abstract
Today’s architecture is often designed based on aesthetics. Among them, it is more common to combine thin shell structures with streamlined designs and create a sense of open space. However, these designs based on aesthetics are not necessarily the most ef- ficient designs. Therefore, this study will combine free-form surface technique and opti- mization theory to design a structure considering both appearance and good mechanical behavior.
In this study, NURBS is used as a free-form surface construction method, and the surface shape is established through the coordinates of control points, node vectors and basis functions. Combined with optimization theory, the coordinates of control points are used as design variables to adjust the structure shape and construct the shape with the best objective function value. This study will also further explore the method of extending the application of shape optimization to support optimization.
In addition to shape optimization, this study will also apply topology optimization in structural optimization. The most efficient material distribution location is gradually found by algorithm. The topology algorithm used in this study is BESO. Different exam- ples will be used to demonstrate the problems of the existing BESO algorithm, and the solutions to these problems will be discussed.
In this study, the topology optimization will be applied to the shape optimization result of thin shell structure to remove the inefficient materials from the structure and keep the efficient parts. Therefore, the materials in the structure can be allocated to the most efficient position for force transmission, and the overall structure performance in terms of stiffness can be further improved with the same amount of materials used.
In this study, the shape optimization and topology optimization methods are applied to different kinds of thin-shell structures to achieve better stiffness performance. The final design will serve as a reference for future engineers when designing thin-shell structures.
Keywords: Free-form surface, Structural optimization, Shape optimization, Topology optimization, Thin shell structure, Shell bridge
目錄
Page
致謝 iii
摘要 v
Abstract vii
目錄 ix
圖目錄 xiii
表目錄 xv
第一章 緒論 1
1.1 研究動機 . . . 1
1.2 文獻回顧 . . . 2
1.3 研究內容 . . . 3
第二章 形狀最佳化與拓樸最佳化方法 5 2.1 前言 . . . 5
2.2 自由曲面建構方法 . . . 5
2.2.1 NURBS . . . 6
2.3 最佳化問題描述 . . . 8
2.4 結構最佳化 . . . 9
2.4.1 形狀最佳化 . . . 10
2.4.2 拓樸最佳化 . . . 11
2.5 結構最佳化分析方法介紹 . . . 11
2.5.1 直接法 (Direct method) . . . 12
2.7 雙向結構最佳化演進法 (Bi-directional Evolutionary Structure Opti-
mization，BESO) . . . 15
2.7.1 方法流程 . . . 15
2.7.2 篩選投影 . . . 16
2.7.3 穩定策略 . . . 17
2.7.4 移除準則 . . . 18
2.7.5 收斂條件 . . . 18
2.8 小結 . . . 19
第三章 形狀最佳化案例探討 21 3.1 前言 . . . 21
3.2 形狀最佳化之程式架構 . . . 21
3.3 平板屋頂案例 . . . 24
3.4 殼橋案例 . . . 26
3.4.1 支承最佳化 . . . 30
3.5 小結 . . . 32
第四章 拓樸最佳化案例探討 33 4.1 前言 . . . 33
4.2 懸臂梁案例 . . . 33
4.3 簡支梁案例 . . . 35
4.3.1 不同初始拓樸形狀 . . . 37
4.3.2 修改穩定策略 . . . 40
4.3.3 限制每步迭代元素改變數量 . . . 42
4.4 上承橋梁案例 . . . 45
4.4.1 收斂結果之判斷 . . . 47
4.4.2 材料集中化 . . . 49
4.5 小結 . . . 52
第五章 拓樸最佳化於薄殼結構之應用 53 5.1 前言 . . . 53
5.2 方法驗證 . . . 53
5.3 方形曲面屋頂案例 . . . 56
5.3.1 殼結構最佳材料配置 . . . 57
5.3.2 加勁層最佳材料配置 . . . 58
5.4 殼橋案例 . . . 63
5.5 小結 . . . 67
第六章 結論與未來展望 69 6.1 結論 . . . 69
6.2 未來展望 . . . 70
參考文獻 71
圖目錄
2.1 曲線與控制點示意圖 (連嘉玟 (2017)) . . . 5
2.2 NURBS 基底函數遞迴關係 (連嘉玟 (2017)) . . . 7
2.3 傳統設計方法 (a) 與結構最佳化設計方法 (b) 之比較 (Arora (2017)) . . 10
2.4 最佳化方法之分類 (Arora (2017)) . . . 11
2.5 棋盤化效應示意圖 . . . 17
3.1 結構最佳化主程式架構 . . . 21
3.2 Scipy 模組功能簡圖 (林享樑 (2018)) . . . 22
3.3 本研究曲面形狀最佳化程式架構 . . . 23
3.4 [例題 3-1] 結構示意圖 . . . 24
3.5 [例題 3-1] 文獻形狀最佳化結果 . . . 24
3.6 [例題 3-1] 控制點分布圖 (單位: 公尺) . . . 25
3.7 Musmeci Bridge 實際結構攝影圖 . . . 26
3.8 殼橋案例結構示意圖 . . . 27
3.9 [例題 3-2] 殼橋案例控制點分布圖 (單位: 公尺) . . . 28
3.10 [例題 3-2] 初始形狀與最佳化結果側視圖比對 . . . 29
3.11 文獻殼橋案例形狀 . . . 29
3.12 [例題 3-3] 殼橋案例支承最佳化控制點分布圖 (單位: 公尺) . . . 31
3.13 [例題 3-2] 與 [例題 3-3] 最佳化結果側視圖比對 . . . 32
4.1 [例題 4-1] 懸臂梁案例結構示意圖 . . . 33
4.2 [例題 4-1] 懸臂梁案例拓樸過程收斂圖 . . . 34
4.3 [例題 4-2] 簡支梁例題結構示意圖 . . . 35
4.4 [例題 4-2] 簡支梁案例拓樸過程收斂圖 . . . 36
4.5 [例題 4-3] 簡支梁案例由軟到硬拓樸過程收斂圖 . . . 38
4.6 限制的方法流程圖 . . . 43
4.7 [例題 4-5] 簡支梁案例限制的方法拓樸過程收斂圖 . . . 43
4.8 [例題 4-6] 上承橋梁案例結構示意圖 . . . 45
4.9 [例題 4-6] 上承橋梁案例不限制的方法拓樸過程收斂圖 . . . 46
4.10 [例題 4-6] 上承橋梁案例限制的方法拓樸過程收斂圖 . . . 47
4.11 材料集中化示意圖，黑色為實心元素，白色為空心元素 . . . 50
4.12 [例題 4-7] 採用不限制的方法拓樸過程收斂圖 . . . 51
5.1 [例題 5-1] 平板案例結構示意圖 . . . 54
5.2 [例題 5-1] 平板案例控制點分布圖 (單位: 公尺) . . . 54
5.3 [例題 5-1] 平板案例最佳化結果俯視圖 . . . 55
5.4 [例題 5-1] 文獻最佳化結果 . . . 56
5.5 [例題 5-2] 方形曲面屋頂案例殼結構最佳材料配置俯視圖 . . . 58
5.6 [例題 5-3] 固定 Vf 而改變 Vbottom下層殼拓樸最佳化結果俯視圖 . . . 61
5.7 [例題 5-3] 固定 Vbottom而改變 Vf 下層殼拓樸最佳化結果俯視圖 . . . 62
5.8 [例題 5-4] Vbottom為 0.3 而 Vf 為 0.5 下層殼拓樸最佳化結果 . . . 64
5.9 [例題 5-4] Vbottom為 0.5 而 Vf 為 0.5 下層殼拓樸最佳化結果 . . . 64
5.10 [例題 5-4] Vbottom為 0.7 而 Vf 為 0.5 下層殼拓樸最佳化結果 . . . 65
5.11 [例題 5-4] Vbottom為 0.5 而 Vf 為 0.3 下層殼拓樸最佳化結果 . . . 66
5.12 [例題 5-4] Vbottom為 0.5 而 Vf 為 0.7 下層殼拓樸最佳化結果 . . . 66
表目錄
3.1 [例題 3-1] 初始形狀與最佳化結果比較 . . . 25
3.2 [例題 3-2] 初始形狀與最佳化結果比較 . . . 28
3.3 [例題 3-2] 固接與鉸接比較 . . . 30
3.4 [例題 3-2] 與 [例題 3-3] 最佳化結果之比較 . . . 31
4.1 [例題 4-1] 懸臂梁案例拓樸最佳化結果 . . . 34
4.2 [例題 4-2] 簡支梁案例拓樸最佳化結果 . . . 36
4.3 [例題 4-3] 不同初始拓樸形狀最佳化結果比較 . . . 38
4.4 [例題 4-3] 簡支梁案例由軟到硬拓樸過程 . . . 39
4.5 [例題 4-4] 簡支梁案例修改穩定策略拓樸最佳化結果 . . . 40
4.6 [例題 4-4] 簡支梁案例修改穩定策略拓樸過程收斂圖 . . . 41
4.7 [例題 4-3] 簡支梁案例由軟到硬第 36 步至第 39 步拓樸形狀 . . . 42
4.8 [例題 4-5] 限制與不限制的方法最佳化結果比較 . . . 43
4.9 [例題 4-5] 簡支梁案例限制的方法拓樸過程 . . . 44
4.10 [例題 4-6] 上承橋梁案例不限制的方法最佳化結果 . . . 46
4.11 [例題 4-6] 上承橋梁案例限制的方法最佳化結果 . . . 46
4.12 [例題 4-6] 限制的方法拓樸過程 . . . 48
4.13 [例題 4-7] 採用不限制的方法最佳化結果 . . . 50
4.14 [例題 4-7] 採用限制的方法最佳化結果 . . . 51
5.1 [例題 5-1] 平板案例最佳化結果 . . . 55
5.2 [例題 5-2] 方形曲面屋頂案例殼結構最佳材料配置最佳化結果 . . . . 57
5.3 [例題 5-3] 固定 Vf 而改變 Vbottom下層殼拓樸最佳化結果 . . . 60
5.4 [例題 5-3] 固定 Vbottom並改變 Vf 下層殼拓樸最佳化結果 . . . 60
5.5 [例題 5-4] 固定 Vf 而改變 Vbottom拓樸最佳化結果 . . . 63
5.6 [例題 5-4] 固定 Vbottom而改變 Vf 拓樸最佳化結果 . . . 65
第一章 緒論
1.1 研究動機
現今建築設計常以美觀、流線型以及創造開放空間感為設計依據，其中薄殼 結構便有別於傳統方正的梁柱結構系統，以自由曲面設計出流線型的的薄殼結 構，然而，基於美學的設計卻未必是最有效率的設計，若要設計兼具美觀與優良 的力學行為之結構，則需要更多的工具提供輔助，因此本研究以最佳化理論結合 自由曲面，建構出在設定目標下具有最好表現之薄殼結構幾何形狀。
本研究將會延續本研究團隊近年對於薄殼結構形狀最佳化之方法，嘗試應用 於不同類型之薄殼結構，並探討支承位置最佳化之方法，而本研究在建構出薄殼 結構之最佳幾何形狀後，會再運用最佳化理論中的拓樸最佳化，進一步將薄殼結 構中材料無效率處移除，僅保留有效率處，藉此找出結構中最有效率之材料配 置，期望藉由應用形狀最佳化與拓樸最佳化方法提供不同的設計概念，作為未來 工程師的參考依據。
1.2 文獻回顧
結構最佳化主要可以分為形狀最佳化、拓樸最佳化以及尺寸最佳化三種類 型，而針對薄殼結構進行形狀最佳化及拓樸最佳化為本研究主要的研究內容，其 中形狀最佳化部分藉由參數化的自由曲面方法來控制結構形狀，常見的方法包括 貝茲曲線Bézier (1968)、均勻 B 樣條 (Uniform B-spline) 以及現今許多商業繪圖軟 體所使用的非均勻有理 B 樣條 (NURBS) Versprille (1975)，而本研究主要是使用 NURBS 來建構自由曲面，後續章節中將會對其進行詳細介紹並以不同例題進行演 示。
運用自由曲面技術結合最佳化方法對薄殼結構進行設計，其發展開始於近十 年，最早是由Kegl and Brank (2006) 提出完整的設計流程，利用貝茲曲線作為曲 面形狀的建構方法，並與最佳化理論結合進行設計，再將此方法應用於桁架結構 上。而 Espath et al. (2011) 則將 NURBS 與最佳化理論結合，並且分別進行線性與 非線性分析，針對不同初始形狀之平板結構進行最佳化設計，最終分別得到線性 與非線性兩種不同之最佳化設計成果。透過此方式將最佳化方法應用於薄殼結構 設計。
在拓樸最佳化方法中，常見的演算法有結構最佳化演進法 (ESO) Xie and Steven (1993)、雙向最佳化演進法 (BESO) Querin et al. (1998) 以及固體等向性懲罰 函數法 (SIMP) Rozvany et al. (1992)，其中 ESO 方法是基於結構承受外力時，並非 所有元素均勻分配應力，利用演算法找出處於無效率狀態之元素並加以剔除，雖 然 ESO 方法簡單有效，但其過程只能逐步移除而不能填回，導致容易求解失敗，
而本研究中所使用的 BESO 最初則延伸自 ESO 方法，對空洞處估算敏感度因子並 提供填回的動作，Huang and Xie (2007) 則對相關研究進行整理，並且提出以平滑 化與投影策略為主之 BESO 架構，成為現行之 BESO 方法。
透過上述文獻的概念，本研究將運用最佳化方法對薄殼結構進行設計，找到 在設定的目標之下，具有最佳表現之結構。而本研究所參考之薄殼結構除了有較 為常見之薄殼屋頂以外，同時也參考位於義大利一座特殊之殼橋 Musmeci Bridge，
利用殼結構對橋面板提供支承，而Briseghella et al. (2016) 則參考 Musmeci Bridge 並用最佳化方法設計一座人行殼橋。
在 最 佳 化 之 研 究 工 具 方 面， 本 研 究 團 隊 過 去 曾 分 別 使 用 Fortran (王 建 凱
(2005))、C++ (李宗豪 (2005)) 以及 MATLAB (施可葳 (2013)) 等多種不同程式語 言，結合商業有限元素分析軟體 ABAQUS 以開發結構最佳化之工具，然而以上 所述之研究工具皆需要讀取以及改寫 ABAQUS 的輸入及輸出檔，建模時的效率 便因此而較低，同時也造成分析所需的時間增加，因此連嘉玟 (2017) 參考Zuo and Xie (2015) 所提出的利用 Python 結合 ABAQUS 進行拓樸最佳化的程式概念，利用 ABAQUS 所提供的 Python 程式語言 API 介面 (Application Programming Interface)，
並且擴充最佳化演算法和自由曲面建模技術至程式中，而林享樑 (2018) 將其程式 架構延伸，將不同的建構曲面方法導入，並且結合 SIMP 演算法進行曲面拓樸最 佳化，並於實際薄殼結構應用 CFD 分析，簡孟笙 (2019) 則開發生成主應力線之演 算法，並擴充至結構最佳化設計中，林家萱 (2020) 探討結合形狀、支承位置以及 厚度的多層次最佳化方法，孫鈺翔 (2021) 則以理論解驗證主應力線加勁之版殼結 構分析結果之正確性。本研究將以本研究團隊近年之程式概念，沿用以 NURBS 曲面建構薄殼結構，搭配 Python 語言之其他相關擴充功能，控制 ABAQUS 進行 有限元素分析，以完成薄殼結構最佳化設計。
1.3 研究內容
第一章闡述本研究之動機及背景，並回顧過去相關的研究文獻。
第二章介紹本研究建構自由曲面的方法，並且介紹本研究於形狀最佳化與拓 樸最佳化所使用之演算法。
第三章先介紹本研究形狀最佳化之程式架構，再以不同薄殼結構案例演示形 狀最佳化之效果，並探討其相關之延伸應用。
第四章以不同平面案例演示拓樸最佳化演算法之特點以及不足之處，並探討 可能的改進方法。
第五章則結合第三章的形狀最佳化結果與第四章的拓樸最佳化方法，探討拓 樸最佳化於薄殼結構之應用。
第六章則為結論與未來展望，總結本研究之內容並提出未來可以進一部探討 之研究方向。
第二章 形狀最佳化與拓樸最佳化方法
2.1 前言
本章節首先將介紹自由曲面之建構方法，藉由參數化曲面來描述幾何形狀，
接著介紹對於最佳化問題之描述，以建立最佳化問題之數學模型 (Mathematical Model)，並且介紹不同結構最佳化之類型，最後則分別介紹本研究於形狀最佳化 以及拓樸最佳化中所使用的演算法。
2.2 自由曲面建構方法
隨著電腦科技的發展以及現代工程對美學要求的提升，許多設計採用平滑 而流線的形狀，因此而發展出許多不同的幾何建模方法，其中較為常見的方法 為貝茲曲線 (Bézier Curve)、B 樣條 (B-spline) 和非均勻有理 B 樣條 (NURBS)，在 這之中又以 NURBS 最為常見，許多不同的商業繪圖軟體所採用的建模方法便為 NURBS，使用者只需要調整少量的控制點座標便可以建構出整個曲線或曲面形 狀，方便使用者對於模型形狀進行調整及修改，曲線與控制點之示意圖可以參考 下圖 2.1 。
圖 2.1: 曲線與控制點示意圖 (連嘉玟 (2017))
由本研究團隊林享樑 (2018) 的研究結論，建議採用 NURBS 作為建構曲面的 方法，可以避免全域效應 (Global Effect) 於調整形狀過程中產生，並可以調整局部 曲面以進行更細微之設計，因此本研究亦選擇使用 NURBS 作為建構曲面的方法，
並將於以下小節針對 NURBS 進行詳細說明，介紹其調整模型的參數以及基底函 數。
2.2.1 NURBS
NURBS(Non-Uniform Rational B-Splines) 曲線之主要概念源自於 B-spline，其 不同之處在於 NURBS 引入對於個別控制點之權重，因此具有更多調整幾何形狀 的自由度。NURBS 曲線的數學定義如下:
C(u) =
∑n i=1
Ni,p(u)wiPi
∑n i=1
Ni,p(u)wi
, u∈ [a, b] (2.1)
其中，Pi 為控制點 (Control Point)，而 Ni,p 則是基底函數:
Ni,0(u) =
1, if ui ≤ u < ui+1
0, otherwise.
(2.2)
Ni,p(u) = u− ui
ui+p− ui
Ni,p−1(u) + ui+p+1− u ui+p+1− ui+1
Ni+1,p−1(u) (2.3) U ={a, . . . , a, up+1, . . . , um−p−1, b, . . . , b}1×(m+1) (2.4)
Ni,p 由式2.4之節點向量 U 決定，wi 為權重，p 為曲線的階數，n 為控制點數 量，m + 1 為節點數量，Ni,0(u)為起始方程，在半開區間 [ui, ui+1)的範圍內其值 為 1，其餘則為 0，Ni,p(u)則是透過兩個 p− 1 階之基底函數線性組合而成，基底 函數具有區間以及遞迴性如下圖 2.2 所示。
圖 2.2: NURBS 基底函數遞迴關係 (連嘉玟 (2017))
若是要由 NURBS 曲線擴展至 NURBS 曲面，則需要定義兩個方向的節點向 量，分別為 u 方向 p 階以及 v 方向 q 階之節點向量，其定義如下:
S(u, v) =
∑n i=1
∑m j=1
Ni,p(u)Nj,q(v)wi,jPi,j
∑n i=1
∑m j=1
Ni,p(u)Nj,q(v)wi,j
, 0≤ u, v ≤ 1 (2.5)
U ={0, . . . , 0, up+1, . . . , ur−p−1, 1, . . . , 1}1×(r+1) (2.6) V ={0, . . . , 0, vq+1, . . . , vs−q−1, 1, . . . , 1}1×(s+1) (2.7)
其中 Pi,j 為雙向的控制點網 (Control Point Net)，Ni,p 為由式 (2.6) 節點向量 U 定義的基底函數，Nj,q 為由式 (2.7) 節點向量 V 定義的基底函數，wi,j 為權重，p 與 q 分別為曲面在 u 向以及 v 向之階數，n 與 m 則分別為 u 向以及 v 向之控制點 數量，r = n + p + 1，r + 1 為 u 向節點數量，s = m + q + 1，s + 1 為 v 向節點數 量。
透過雙向控制點網便可以建立平滑的 NURBS 曲面形狀，而本研究便透過改 變控制點座標來調整曲面形狀，結合後續小節所介紹的最佳化方法，達成在特定 目標之下具有最佳表現之曲面形狀。
2.3 最佳化問題描述
為了將牽涉較多複雜數學條件的最佳化問題描述清楚，需要將最佳化之目標 與限制條件描述為標準形式的方程組，而Arora (2017) 將最佳化問題的標準數學描 述式定義如下:
尋找一組 n 維之設計變數 (Design Variables x = (x1, x2, . . . , xn)) ，在給定設計 變數之上限與下限的情形之下，並滿足等式限制式 (Equality Constraints) 以及不等 式限制式 (Inequality Constraints)，使目標函數 (Objective Function) 能夠具有最大值 (Maximize) 或是最小值 (Minimize)。上述最佳化問題之定義可以表達成數學描述 式如下:
Min. or Max. :
f (x) = f (x1, x2, . . . , xn) (2.8)
Subject to :
hi(x) = hi(x1, x2, . . . , xn) = 0, i = 1, 2, . . . , p (2.9) gi(x) = gi(x1, x2, . . . , xn)≤ 0, j = 1, 2, . . . , q (2.10) (xmin)k ≤ xk ≤ (xmax)k, k = 1, 2, . . . , n (2.11)
其中 f (x) 是最佳化之目標函數，h(x) 為等式限制式，g(x) 為不等式限制式，
p與 q 則分別為等式與不等式限制式之數量。在 n 個設計變數當中，(xmin)k 和 (xmax)k分別為第 k 個設計變數 xk之上限及下限。
一般可以依限制條件之有無對最佳化問題進行分類，有限制條件之最佳化問 題，稱為束制最佳化問題 (Constrained Optimization Problem)；而沒有限制條件之 最佳化問題，則稱為無束制最佳化問題 (Unconstrained Optimization Problem)。
除了依照限制條件之有無進行分類以外，最佳化問題也可以依照目標函數及 限制條件是否皆為設計變數 x 之線性方程式做分類，若皆為線性方程式者稱為線
性規劃問題 (Linear Programming Problem)；反之，若有任何一個方程式為設計變 數之非線性方程式，則歸類為非線性規劃問題 (Nonlinear Programming Problem)。
由於一般工程問題除了設計變數有範圍限制以外，還需要符合實際物理限 制條件，因此大多屬於束制最佳化問題。此外，大多數結構最佳化問題都歸類 為非線性規劃問題，因此各方程式的解析式 (Analytic Form) 難以直接求出，而 且目標函數是否為平滑連續的函數型式亦不容易確認，導致求解一般結構最佳 化問題時，並不適合使用過往的變分法 (Variational Method) 求極值的方法 (Kamat (1993))。目前求解非線性的最佳化問題時，會根據目標函數之特性選擇適合的數 值方法 (Numerical Method)，而大多數的演算法則需要經過多次迭代計算 (Iterative Calculation) 以漸漸求得最佳解 (Optimal Solution)。
2.4 結構最佳化
下圖 2.3 為傳統設計方法與結構最佳化設計方法之比較，兩者之間最主要的 不同之處在於更新設計變數的方法，傳統設計方法是設計者依照經驗來判斷並更 新設計變數，而結構最佳化方法則是以最佳化演算法之理論為基礎，對設計變數 進行迭代並更新直到滿足收斂條件為止。相較於傳統設計方法需要人為以經驗判 斷來更新，結構最佳化方法則是利用電腦進行迭代，因此可以將人為判斷之需求 省略，能夠以更有效率且穩定之方法進行結構設計。
一般求解最佳化設計問題之流程可以參考下圖 2.3(b) ，首先要先將欲求解之 問題建立為最佳化問題，並依據問題設定初始設計以及建立分析模型，再對最佳 化問題之設計變數、目標函數以及限制條件做出定義。接著進行系統分析，結構 最佳化之系統分析包括有限元素分析以及最佳化分析兩個部分。每一步迭代的最 佳化分析結果將提供一組新的設計變數，並將其更新至分析模型以重覆進行系統 分析。當最佳化問題滿足收斂條件時，則視為達到最佳化設計。在整個結構最佳 化的過程中，每一步迭代過程皆需要進行有限元素分析以及最佳化分析，因此是 結構最佳化設計中最為困難且重要的部分。
圖 2.3: 傳統設計方法 (a) 與結構最佳化設計方法 (b) 之比較 (Arora (2017))
依照處理類型之不同可以將結構最佳化設計分為以下三種：尺寸最佳化 (sizing optimization)、形狀最佳化 (shape optimization) 以及拓樸最佳化 (topology optimization)。以上三種結構最佳化類型之設計流程大略相似，都需要先定義初始 模型，並重複進行有限元素分析以及最佳化分析直到滿足收斂條件，主要的不同 之處在於設計變數以及使用之演算法。本研究中主要使用形狀最佳化以及拓樸最 佳化作為結構最佳化之設計方法，並將於以下小節中進一步做介紹。
2.4.1 形狀最佳化
形狀最佳化是藉由改變結構之幾何形狀以達到目標函數的最佳解，最小化結 構順從度通常為其目標函數。
形狀最佳化可以應用於不少結構系統中，常見的如桁架結構系統，透過改變 節點座標使桁架系統在相同的受力情況之下能夠具有較大之勁度。
而本研究團隊近幾年的研究內容，則是將曲面形狀最佳化應用於板殼結構，
運用第 2.2 節中所介紹的自由曲面建模方法，透過控制點及基底函數、節點向量 等控制參數構建出曲面形狀，藉由改變控制點座標來調整曲面的幾何形狀，以達 到單位體積勁度最大化，因此控制點座標即為最佳化之設計變數，同時也會搭配 位移或應力限制式進行分析。
2.4.2 拓樸最佳化
拓樸最佳化透過演算法得到最有效率之材料分布位置，使結構在相同的材料 用量之下具有最大的勁度，讓材料集中在最有效率之位置，也就是力量傳遞的重 要路徑上。
拓樸最佳化之目標函數通常為最小化結構順從度，設計變數則為元素密度，
因此在拓樸最佳化問題中，設計變數的數量往往數以千計，為了增加計算效率 以處理大量數據，拓樸最佳化演算法與一般最佳化演算法有所不同，需要利用 獨特的設計變數更新準則，常見的拓樸最佳化演算法有結構最佳化演進法 ESO (Evolutionary Structural Optimization)、雙向結構最佳化演進法 BESO (Bi-directional Evolutionary Structural Optimization)、元素交換法 EEM (Element Exchange Method) 以及固體等向性懲罰函數法 SIMP (Solid Isotropic Materal with Penalization)，其中，
本研究中所採用的拓樸最佳化方法為雙向結構最佳化演進法 BESO，並將於本章 後續小節中詳細介紹。
2.5 結構最佳化分析方法介紹
結構最佳化開始發展於 1980 年代，其分析方法種類眾多，本節參考Arora (2017) 將非線性規劃問題分類為以下兩大數值方法：直接法 (direct method) 以及間 接法 (indirect method)，其分類方式可以參考下圖 2.4 。
圖 2.4: 最佳化方法之分類 (Arora (2017))
其中，直接法的設計變數數量通常較少，更新設計變數的準則也較為簡單，
設計變數之間比較不會彼此互相影響；而間接法設計變數數量則通常較多，通常 各個方法更新設計變數的方式都會有所不同，因此需要藉由不同的更新設計變數 準則，將所需計算的參數數量減少，可以大幅節省計算時間，因此適用於大型的 變數矩陣運算。以下小節將針對本研究於形狀最佳化中所使用的直接法進行詳細 介紹。
2.5.1 直接法 (Direct method)
直接法首先要給予設計變數初始數值 (initial design variable)，並利用各種數學 規劃法以計算每一步迭代過程中對於設計變數之微小擾動量 ∆x(k)，其中，此微小 擾動量可以分為方向 (direction) 以及步長 (step size) 兩個部分，直接法便藉由此機 制來勘查在設計變數上下限範圍內目標函數之起伏變化。最佳化迭代過程便是透 過尋找出第 k 次迭代時之方向以及沿著該方向之步長，計算而得到下次迭代的新 設計變數 x(k+1)，如下式 (2.12) 所示:
x(k+1)= x(k)+ ∆x(k) (2.12)
設計變數組隨著每一步迭代而漸漸接近最佳解，當設計變數不再有新的改變 量時，便滿足收斂條件並且視為得到最佳解。
數學規劃法可以應用於許多不同最佳化問題，而對於一般的工程問題來說，
由於設計範圍較廣且限制條件較多，因此相較於直接搜索的演算方式，循序式的 求解方法會更加受到青睞，其中，循序線性規劃法 (Sequential Linear Programming, SLP) 以及循序二次規劃法 (Sequential Quadratic Programming, SQP) 便為兩種常見 循序式的數學規劃法，以上兩種方法都將較為複雜之非線性限制條件拆解成多組 計算相對單純的子問題，兩者不同之處在於前者是拆解為線性規劃子問題，後者 則是拆解為二次規劃子問題。
非束制最佳化問題與束制最佳化問題之迭代過程邏輯相同，兩者最大的不同 之處在於束制最佳化問題在迭代過程中，在決定搜索方向和步長的同時，也需要 將限制條件之約束納入考量，也就是可行之方向以及可行之步長。束制最佳化問
題會利用懲罰函數法來處理違背限制條件之情形，以最小化目標函數之問題為 例，倘若某一組設計變數違背限制條件時，將會藉由懲罰函數 (penalty function) 來增加目標函數的數值，以作為對該組設計變數之懲罰，一般常用動態的懲罰函 數 (dynamic penalty function) 對違背限制條件之目標函數值進行調整，也就是對於 超過限制條件越遠者，所給予的懲罰函數也會越大，藉以將此組設計變數之貢獻 排除；而無束制最佳化問題則不會受到限制條件的約束，因此可以直接對於方向 以及步長進行計算，讓目標函數最小化。
以梯度法作為基礎的數學規劃法，與勘查地形的方式相似，但是當最佳化問 題的設計變數數量增加時，直接法求解所需要的迭代次數往往會較多，並且當接 近最佳解時，收斂的速度經常會趨於緩慢 (Kirsch (2012))。
本研究於曲面形狀最佳化中所選用之演算法為數學規劃法中的序列最小二次 數學規劃法 SLSQP (Sequential Least Squares Programming optimization algorithm)，
也就是循序二次規劃法 (SQP) 的一種，以擬牛頓法 (Quasi-Newton Method) 為原 理，並將牛頓法於計算二次微分時效率較低的問題改良。
2.6 結 構 最 佳 化 演 進 法 (Evolutionary Structure Opti- mization，ESO)
本研究中所使用的拓樸最佳化演算法為雙向結構最佳化演進法，其概念則是 源自於結構最佳化演進法 (ESO)，因此將先於本節對 ESO 方法之原理以及流程進 行介紹。
ESO 方法之概念原理為當結構承受外力作用時，由於結構中每個元素所承受 的應力大小會有所不同，位於力量傳遞重要路徑之元素會承受較大之應力，而所 承受應力較小之元素即視為無效率之元素，因此利用演算法逐步將其移除，僅保 留有效率之元素，使結構可以漸漸演進至最佳拓樸形狀，ESO 之方法流程可以整 理如下:
Step 1 : 將初始設計領域切分成有限元素網格。
Step 2 : 進行有限元素分析。
Step 3 : 根據拓樸最佳化之目標函數，計算每個元素的敏感度因子並進行排序。
Step 4 : 依據排序結果，決定本次迭代中的無效率元素並加以移除。
Step 5 : 重複步驟 2 至步驟 4 直到滿足收斂條件，即可得到最佳拓樸形狀。
上述步驟中若目標函數為最小化結構順從度，則以元素之結構順從度 (也就是 兩倍應變能) 為其敏感度因子 αe，移除元素的方式則是以移除比率 RR (Removal Ratio) 作為移除門檻，並透過以下式 (2.13) 將敏感度因子小於移除門檻之無效率 元素移除，並保留敏感度因子大於移除門檻之元素。
αi
αoverall < RR (2.13)
式中 αi 為第 i 個元素之敏感度因子，αoverall 則為整體結構之敏感度因子，當 元素符合上式之條件時便將其移除。
移除無效率之元素後，為了進行下一步迭代，需要透過演進速率 ER (Evolu- tional Ratio) 對移除比率 RR 進行更新，再用新的移除比率 RR 來移除下一步無效 率之元素，直到達到預設的收斂條件為止，移除比率 RR 的更新方式如下:
RRnew = RRold+ ER (2.14)
ESO 常見的收斂條件則有以下三種: 最大演進次數 IT ERmax、最大移除比率 RRmax 以及最大移除率 RVmax ，其中，移除率 (Removed Volume Ratio) 為目前結 構中實心元素個數佔原有元素個數的百分比。當拓樸演進至滿足所選用的停止條 件時，即視為達到最佳拓樸形狀。
2.7 雙向結構最佳化演進法 (Bi-directional Evolutionary Structure Optimization，BESO)
第 2.6 節中所介紹的 ESO 方法雖然簡單有效，但由於其演進過程只能單向性 的移除元素，因此若是遇到某次演進步能量分布變化過大時，將無法透過填回元 素進行補償，容易造成拓樸結果失敗，因此雙向結構最佳化演進法 BESO 便被提 出以改進此問題，提供已移除之元素填回的動作。
BESO 最初利用位移外插 (Displacement Extrapolation) 方法估算相接於現有元 素邊界但不存在元素網格處之敏感度因子，作為衡量元素填回之標準，而 BESO 作者Huang and Xie (2007) 整理並提出以平滑化與投影策略為概念之方法成為 現行 BESO 之架構；Zuo and Xie (2015) 則將 BESO 方法以 Python 程式語言連接 ABAQUS 進行拓樸最佳化，並透過將元素之楊氏模數乘以 0.001 以模擬空心元素 之行為，以達到軟移除的效果，而本研究便以 Zuo and Xie 的方法作為主要程式架 構，將於本節詳細介紹。
2.7.1 方法流程
BESO 之方法流程大致與 ESO 相似，兩者主要不同之處在於對元素敏感度因 子之處理、收斂條件以及初始材料分布不一定填滿設計領域，BESO 之方法流程 可以整理如下:
Step 1 : 將初始設計領域離散為有限元素網格。
Step 2 : 設定初始材料分布，建議但不一定為全填滿設計領域。
Step 3 : 根據題目選取適當之篩選半徑 rmin、容許誤差 τ 以及考量收斂之步數 N，並 用篩選半徑計算各元素之權重。
Step 4 : 進行有限元素分析。
Step 5 : 根據拓樸最佳化之目標函數，計算每個元素的敏感度因子。
Step 6 : 計算出敏感度因子後，再配合步驟 3 算出之權重進行篩選投影，再以穩定策
Step 7 : 依據排序結果，以移除準則決定本次迭代中的無效率元素加以移除，並將有 效率之元素填回。
Step 8 : 重複步驟 4 至步驟 7 直到滿足收斂停止條件，即可得到最佳拓樸形狀。
上述步驟中所提到的篩選投影、穩定策略、移除準則以及收斂停止條件將會 於以下小節中進一步詳細介紹。
2.7.2 篩選投影
篩選投影為將篩選半徑範圍內的元素原始敏感度因子進行加權平均，可以對 空網格處元素敏感度因子進行估算，同時也具有平滑化元素敏感度因子之效果，
可以避免棋盤化效應之產生，如下圖 2.5 所示。針對元素 i 之篩選投影計算方式如 下式 (2.15) 以及式 (2.16):
αi =
∑M j=1
w(rij)αej
∑M j=1
w(rij)
(2.15)
w(rij) =
rmin− rij , rij < rmin
0 , rij ≥ rmin
(2.16)
式中 αi為元素 i 之敏感度因子，以元素 i 為中心，j 為篩選半徑 rmin 範圍內 之元素 (共有 M 個)，αje為元素 j 之原始敏感度因子，rij 為元素 j 中心至元素 i 中 心之距離，w(rij)則為線性權重函數。
(a) 棋盤化效應之結構
(b) 經過篩選投影平滑化後之結構 圖 2.5: 棋盤化效應示意圖
篩選半徑之大小會影響拓樸最佳化結果之桿件細緻程度，採用較小之篩選半 徑可以保留較為細小之桿件，然而拓樸結果可能會過於複雜；而採用較大之篩選 半徑則可以將較為細小之桿件消除，只保留較大之桿件，然而過大之篩選半徑則 可能會干擾拓樸形狀，得到較差之結果，篩選半徑之選擇主要取決於欲設計之桿 件尺寸。
2.7.3 穩定策略
由於移除及填回元素時，忽略了同時變更多個元素的交互影響，因此容易形 成數個重複拓樸形式不斷交替出現之情形，造成最後會不停迭代而無法收斂，因 此需要利用穩定策略來降低此情形發生之可能性，其方法為將此步計算出的敏感 度因子與前一步之敏感度因子進行平均，計算方式如下式 (2.17):
αi = αki + αki−1
(2.17)
式中 k 為目前演進之步數，αik與 αik−1 分別為元素 i 此步與前步之敏感度因 子，而平均後之 αi即為元素 i 於本次排序所使用之敏感度因子。
2.7.4 移除準則
由於 BESO 方法具有填回元素之自由，並且其初始設計也不一定為全設計領 域填滿實心元素，因此 BESO 之移除準則會與原始 ESO 略有不同。在移除及填回 元素之前，需要先計算每步演進之體積，而本研究採用固定體積比例的方法更新 每步演進之體積，此方法是由目前之體積改變一定比例作為下一步之體積，能夠 確保每次演進的體積改變量為相對較小之更動，以避免影響最終拓樸結果，其更 新體積方式可以表示如下:
Vi =
M ax(Vi−1(1− ER), Vtarget), when Vi−1 > Vtarget M in(Vi−1(1 + ER), Vtarget), when Vi−1 ≤ Vtarget
(2.18)
其中 Vi 代表第 i 次演進之體積，ER 為固定體積改變比例，一般大多設定為 0.02，Vtarget 則為目標體積。
算出此步演進之體積後，先將所有元素之敏感度因子進行排序，並取敏感度 因子較大者作為實心元素，取到實心元素之體積總和等於此演進步之體積為止，
其餘元素則為空心元素，透過此方法達到移除及填回元素之目的。
2.7.5 收斂條件
與 ESO 方法不同，由於 BESO 方法可以對元素進行移除及填回，所以達到目 標體積後可以繼續進行迭代以求得最佳拓樸形狀，因此其停止條件並非達到目標 體積便立即停止，而是改成以目標函數是否趨近於收斂作為判斷標準，其數學式 可以表示如下:
Error = |∑N
i=1(Ck−i+1− Ck−N−i+1)|
∑N
i=1Ck−i+1 ≤ τ (2.19) 其中 C 為系統之結構順從度，τ 為容許誤差，N 為考量收斂之步數，k 為目 前演進步數，當誤差小於容許值即視為達到收斂。
2.8 小結
本章節主要介紹本研究中所使用的方法以及理論，包括建構自由曲面所使用 的 NURBS，接著對最佳化問題進行描述，並將結構最佳化分為尺寸最佳化、形 狀最佳化與拓樸最佳化三個類型，最後再詳細介紹本研究於形狀最佳化所使用之 SLSQP 以及拓樸最佳化所使用之 BESO 和其前身 ESO 演算法，而本研究後續章 節將會以不同案例對本章所介紹之方法進行探討，以了解其相關延伸應用以及現 有方法不足之處，並進一步探討可能的改進方法，以提供未來設計薄殼結構時作 為參考。
第三章 形狀最佳化案例探討
3.1 前言
本章節主要將利用第二章所介紹的形狀最佳化方法進行案例探討，首先將介 紹本研究團隊形狀最佳化之程式架構，透過 Python 程式語言控制有限元素分析軟 體 ABAQUS 進行建模與分析，並將分析結果回傳給演算法以進行形狀最佳化；接 著再以不同案例演示薄殼結構形狀最佳化之效果，最後再進一步探討由形狀最佳 化延伸應用於支承最佳化之方法。
3.2 形狀最佳化之程式架構
ABAQUS 為一套有限元素分析軟體，常用於機械、土木等領域之結構分 析，並且具有 API (Application Programming Interface) 介面，提供使用者以 Python 對 ABAQUS 進 行 控 制， 而 本 研 究 團 隊 自 連 嘉 玟 (2017) 開 始 便 採 用 Python 與 ABAQUS 連結的方式進行最佳化，如下圖 3.1 所示，不同於過去使用其他程式語 言時，需要對輸入及輸出檔進行改寫及讀取，因此可以提高分析效率。
圖 3.1: 結構最佳化主程式架構
此外，Python 也提供使用者引入第三方模組庫 (third-party library) 以對程式進 行擴充，例如本研究中所使用的自由曲面工具 NURBS 以及最佳化演算法 SLSQP，
便是引入自 Python 第三方模組庫，分別為 nurbs 模組以及最佳化工具 Scipy 模組 庫中的 optimize 模組 (如下圖 3.2 所示)。
圖 3.2: Scipy 模組功能簡圖 (林享樑 (2018))
本研究形狀最佳化之程式架構可以參考下圖 3.3 ，其中包含模型建立函式、
模型性質與邊界條件設定函式、有限元素分析設定函式以及最佳化分析函式，四 個函式於最佳化迴圈中重複執行，直到達到收斂條件。
其中，模型建立函式是依照第 2.2 節所介紹的方法，透過設定控制點座標、
階數及節點向量等參數，並經由計算得到 NURBS 曲面，再將計算之曲面離散成 點座標，並將點座標一一輸入 ABAQUS 成為元素之節點座標，完成後即可得到 符合 NURBS 曲面形狀之模型。最佳化迭代過程中每當設計變數更新時，便透過 此函式改變模型形狀。
模型性質與邊界條件設定函式可以調整 ABAQUS 模型之參數，包含材料性 質、邊界條件、斷面設定、載重形式以及綁定約束等，可以利用此函式根據最佳 化問題之設定，對模型參數進行調整。
模型形狀建立完成並設定好參數後，再透過有限元素分析設定函式控制 ABAQUS 進行分析，並讀取分析結果以進行後續處理，如應變能、位移及應力 等，作為最佳化演算法之目標函數及限制式。
最佳化分析函式則將有限元素分析設定函式所回傳之目標函數及限制式，代 入由 optimize 模組所引入之最佳化演算法 SLSQP 以進行最佳化分析，最後再將演 算法所更新之設計變數代回模型建立函式以進行下一次迭代。
圖 3.3: 本研究曲面形狀最佳化程式架構
3.3 平板屋頂案例
[例題 3-1]
本節參考 Espath et al. (2011) 之例題，初始形狀為邊長 10 公尺，厚度 30 公 分之正方形平板屋頂，材料參數楊氏模數為 30 GPa，柏松比為 0.3，密度為 2500 kg/m3，邊界條件為四個角落的鉸支承，承受的載重為自重，本例題之結構示意 圖可以參考下圖 3.4 ，而文獻之形狀最佳化結果則可以參考下圖 3.5 。
本例題使用 64 個控制點描述曲面幾何形狀，由於整體結構為 1/8 對稱，並且 支承所在的角落紅色控制點於最佳化過程中維持固定不動，因此獨立的設計變數 為 9 個控制點的 z 座標，控制點之分布可以參考下圖 3.6 ，其中藍色控制點代表 z 座標為設計變數之控制點，其餘黑色控制點座標則可以透過對稱得到。
圖 3.4: [例題 3-1] 結構示意圖
圖 3.5: [例題 3-1] 文獻形狀最佳化結果
圖 3.6: [例題 3-1] 控制點分布圖 (單位: 公尺)
在厚度不變的情形下，形狀變化後表面積必然大於等於初始平板，因此體積 也會大於等於初始平板，然而體積增加會導致自重增加，造成總載重與初始平板 不同而難以直接比較最佳化形狀之效果，因此為了方便比較初始形狀與最佳化形 狀，最佳化過程中會藉由調整厚度使體積維持不變，最佳化結果與初始形狀之比 較如下表 3.1 :
表 3.1: [例題 3-1] 初始形狀與最佳化結果比較
Initial Optimization Result
Shape
C (N× m) 13686 214.28
V (m3) 30.0 30.0
CV (N× m4) 410588 6428.41
umax(m) 0.026 0.00035 uavg (m) 0.019 0.00030
由上表 3.1 可知，在相同材料使用量 (體積一樣) 的情況下，經過形狀最佳化 改變形狀後，整體結構的結構順從度約為初始平版的 0.016 倍，也就是結構勁度 提升了 64 倍，同時，結構的最大位移與平均位移也分別為初始平版的 0.013 倍與 0.016 倍，顯示形狀最佳化對於提升結構勁度具有良好的效果，比較本研究最佳化 形狀與參考文獻之形狀最佳化結果，也可以發現兩者形狀相似，顯示本研究形狀 最佳化方法具有可信度。
3.4 殼橋案例
本節殼橋案例所參考的是Briseghella et al. (2016) 的例題，而文獻之例題則是 參考自位於義大利的 Musmeci Bridge，實際結構照片可以參考下圖 3.7 ，Musmeci Bridge 為一座特殊的殼橋，其設計為將殼結構作為橋梁之下部結構以提供橋面板 支承，由於設計時需要兼顧美觀以及提供人行之目的，因此需要於橋面板與下部 殼結構之間保留足夠空間供行人行走。
圖 3.7: Musmeci Bridge 實際結構攝影圖
[例題 3-2]
本例題主要參考Briseghella et al. (2016) 的例題，設計一個由殼結構所支撐的 橋，支承位置則設在下方殼結構的四個角落與上方橋面板兩端，殼再於跨距 1/3 與 2/3 的位置與上方橋面板相接，以提供橋面板支承 (如下圖 3.8 所示)，所有支承 皆為鉸支承，板與殼相接處則為固接，橋梁跨距為 40 公尺，寬度為 6 公尺，高度 為 6 公尺，橋面板厚 1 公尺，殼厚 0.3 公尺，橋面板與殼結構採用相同的混凝土材 料，楊氏模數為 30 GPa，柏松比為 0.3，密度為 2500 kg/m3，載重則為自重以及 參考公路橋梁設計規範的人行道活載重並取保守值 3920 N/m2加載於橋面板上。
圖 3.8: 殼橋案例結構示意圖
形狀最佳化之設計領域為下方殼結構，上方橋面板於最佳化過程中固定不 變。為了確保建模時殼與橋面板有確實相接在一起，需要將相接處殼結構的元 素節點固定在與橋面板相同高度處，因此建模時需要將殼結構的形狀分成三個 NURBS 曲面來表達，並將殼結構與橋面板相接處設為 NURBS 曲面之角落，如此 便可以藉由固定角落控制點使殼結構與橋面板能確實相接。由於對稱的關係，第 一個與第三個 NURBS 曲面形狀相同，此外，橋梁沿軸向的兩側亦為對稱，因此 獨立的設計變數為 8 個控制點 z 座標，控制點分布與 NURBS 曲面劃分如下圖 3.9 所示。
圖 3.9: [例題 3-2] 殼橋案例控制點分布圖 (單位: 公尺)
本例題與平板屋頂案例相同，為了方便比較初始形狀與最佳化形狀，最佳化 過程中也會藉由調整殼結構厚度使總體體積維持不變，另外，由於本例題所參考 之殼橋有美觀上的需求，因此會限制控制點最大高度需要比橋面版低 3 公尺。本 例題之最佳化結果與初始形狀之比較如下表 3.2 ，初始形狀與最佳化形狀之側視 圖則可以參考下圖 3.10 ，其中，圖 3.10 與表 3.2 中之形狀皆只截取殼結構之形狀 而隱藏橋面板，以方便觀察最佳化結果:
表 3.2: [例題 3-2] 初始形狀與最佳化結果比較
Initial Optimization Result
Shape
C (N× m) 4.52× 105 1.99× 105
V (m3) 316.6 316.6
CV (N× m4) 1.43× 108 6.31× 107
umax(m) 0.092 0.038
uavg (m) 0.057 0.023
(a) 初始形狀側視圖
(b) 最佳化形狀側視圖
圖 3.10: [例題 3-2] 初始形狀與最佳化結果側視圖比對
圖 3.11: 文獻殼橋案例形狀
由上表 3.2 可知，在相同材料使用量 (體積一樣) 的情況下，經過形狀最佳化 後，整體結構的結構順從度約為初始形狀的 0.44 倍，也就是結構勁度提升為 2.27 倍，同時，結構的最大位移與平均位移也分別為初始形狀的 0.41 倍與 0.40 倍，顯 示形狀最佳化能提升整體結構之表現，比較最佳化形狀與參考文獻結果，也可以 發現兩者形狀相似。
表 3.3: [例題 3-2] 固接與鉸接比較
固接 鉸接
C (N× m) 99618.5 100476 殼角落 σmax(N/m2) 5.30×107 5.31×107 板與殼相接 σmax(N/m2) 1.32×107 1.16×107 殼角落 σmin(N/m2) -8.04×107 -8.05×107 殼角落支承反力 (N) 4.63×106 4.63×106
另外，觀察表3.3中的結果可以發現，若是將殼結構與橋面板相接處由固接改 為鉸接，雖然可以將殼結構相接處的最大張應力由 13.2 MPa 降低為 11.6 MPa，然 而降低後的應力值仍然超過混凝土抗拉強度，同時，殼結構的張應力最大值皆位 於角落接地處，並且會由 53.0 MPa 提高為 53.1 MPa，而整體結構的結構順從度 也會提高為 100476 N× m，以上結果顯示將殼結構與橋面板相接處由固接改為鉸 接，不僅無法解決相接位置應力集中問題，同時還會使結構其他區域應力提高，
整體結構的表現也較差，因此後續殼橋案例維持將殼結構與橋面板相接處設為固 接。
3.4.1 支承最佳化
[例題 3-3]
由例題 3-2 可以發現，NURBS 曲面角落控制點可以直接代表殼結構與橋面 板相接的所在位置，雖然其高度必須限制與橋面版相同而不能作為設計變數，然 而，若是將控制點 x 座標設為設計變數，則代表可以找出殼結構與橋面板沿著軸 向的最佳相接位置，達到支承最佳化的效果，因此可以同時進行形狀最佳化與支 承最佳化。
本例題與例題 3-2 相同，為一跨距 40 公尺之殼橋，不同之處在於設計變數新 增控制點 x 座標，如下圖 3.12 中所示，本例題位於殼結構角落的紅色控制點仍然 維持不變，而殼結構與橋面板相接處的紫色控制點之 x 座標則為新增的獨立設計 變數，其餘控制點之 x 座標則平均分布於紅色與紫色控制點之間，在減少設計變 數的同時，也可以避免控制點分布過於集中。
圖 3.12: [例題 3-3] 殼橋案例支承最佳化控制點分布圖 (單位: 公尺)
本例題之形狀與支承最佳化結果和例題 3-2 形狀最佳化結果之比較可以參考 下表 3.4 ，而本例題與例題 3-2 最佳化形狀之側視圖比對則可以參考下圖 3.13 ， 與例題 3-2 相同，圖 3.13 與表 3.4 中之最佳化形狀只截取殼結構之形狀而隱藏橋 面板，以方便觀察最佳化結果:
表 3.4: [例題 3-2] 與 [例題 3-3] 最佳化結果之比較
Only Shape Optimization Result Shape and Support Optimization Result
Shape
x (m) 6.67 9.11
C (N× m) 1.99× 105 1.72× 105
V (m3) 316.6 316.6
CV (N× m4) 6.31× 107 5.46× 107
umax(m) 0.038 0.034
uavg (m) 0.023 0.019
(a) [例題 3-2] 形狀最佳化結果側視圖
(b) [例題 3-3] 形狀與支承最佳化結果側視圖
圖 3.13: [例題 3-2] 與 [例題 3-3] 最佳化結果側視圖比對
由上表 3.4 可知，在相同材料使用量 (體積一樣) 的情況下，本例題經過形狀 與支承最佳化改變形狀與支承位置後，整體結構的結構順從度約只有形狀最佳化 的 0.86 倍，也就是結構勁度提升了 16%，同時，結構的最大位移與平均位移也分 別比只有形狀最佳化的降低 11% 與 17%，顯示形狀與支承最佳化能進一步提升整 體結構之表現。另外，比較本例題與例題 3-2 之側視圖，也可以發現本例題最佳 化形狀之殼結構與橋面版相接位置更靠近兩端，顯示原始設計支承位於跨距 1/3 與 2/3 的位置並非最佳支承位置。
3.5 小結
本章節先介紹了本研究形狀最佳化之程式架構，再分別以平板屋頂案例以及 殼橋案例演示，不僅可以成功找到具有較佳表現之結構形狀，最佳化結果也與參 考文獻相似，同時本章再進一步將形狀最佳化的方法延伸，探討應用於支承最佳 化的方法，例題也顯示同時進行形狀最佳化與支承最佳化之結果，會比只有進行 形狀最佳化之結果具有更好的表現。
第四章 拓樸最佳化案例探討
4.1 前言
本章將使用第二章介紹的拓樸最佳化方法進行案例探討，分別以簡支梁、懸 臂梁以及上承橋梁等平面案例演示，探討拓樸最佳化演算法之特點，並嘗試對現 有演算法之缺點進行改進，期望能更穩定且可靠地找出最佳拓樸形狀，使材料分 布於設計領域中最有效率處。
4.2 懸臂梁案例
[例題 4-1]
本例題之結構示意圖可以參考下圖 4.1 ，設計領域為長 16 公分、深 10 公分、
厚度 0.1 公分的懸臂梁，懸臂梁自由端中點承受向下集中載重 10 kN，設計領域切 分為邊長 0.2 cm 之元素，實心元素之楊氏模數為 200 GPa，柏松比為 0.3，並將楊 氏模數乘以 0.001 來代表空心元素，以達到對元素軟移除之目的。
圖 4.1: [例題 4-1] 懸臂梁案例結構示意圖
本例題 BESO 之參數為篩選半徑 rmin=1.0cm，目標體積則設為 40%，本例題 之初始拓樸形狀與拓樸最佳化結果之比較可以參考下表 4.1 ，表中之體積數值皆 會以相對於總體積的百分比顯示，而拓樸收斂過程則可以參考下圖 4.2 。
表 4.1: [例題 4-1] 懸臂梁案例拓樸最佳化結果
Initial Optimization Result
Shape
C (kN× mm) 12.04 27.82
V (%) 100.0 40.0
CV (kN× mm) 12.04 11.13
(a) 結構順從度收斂圖 (單位: kN× mm) (b) 體積收斂圖 (無單位)
圖 4.2: [例題 4-1] 懸臂梁案例拓樸過程收斂圖
由上表 4.1 可知，經過拓樸最佳化移除 60% 的體積後，結構順從度提升為初 始形狀的 2.31 倍，雖然這代表拓樸結果的整體結構勁度下降，然而結構順從度與 體積之乘積則下降了約 7.7%，代表單位體積之結構勁度上升，顯示拓樸最佳化可 以找出無效率元素的所在位置並加以移除，使材料分布於最有效率處。
本例題之拓樸最佳化結果顯示，現有之演算法對於簡單的懸臂梁例題，可以 順利找到最佳拓樸形狀，然而，若是進一步應用於更為複雜之問題則仍然有不足 之處，因此以下將以例題演示現有演算法之缺點，並且嘗試對其進行改進。
4.3 簡支梁案例
本節先以簡單的簡支梁例題為範例，將拓樸最佳化結果作為對照，並嘗試改 變初始拓樸形狀，以測試目前之演算法對於不同初始拓樸形狀，是否可以順利地 收斂至一致的解。
[例題 4-2]
本例題之設計領域為長 12 公尺、深 6 公尺、厚度 0.1 公尺的簡支梁，簡支梁 跨距中點下緣承受向下集中載重 30 kN，跨距四分之一與四分之三位置的下緣則 承受向下集中載重 15 kN，設計領域切分為邊長 0.1 公尺之方形元素，實心元素的 楊氏模數為 210 GPa，柏松比為 0.3，並將楊氏模數乘以 0.001 來代表空心元素。
本例題結構示意圖可以參考下圖 4.3 ，BESO 參數為篩選半徑 rmin=0.3 公尺，目標 體積則為 30%，本例題之初始拓樸形狀與拓樸最佳化結果之比較可參考下表 4.2 ， 拓樸收斂過程則可以參考下圖 4.4 。
圖 4.3: [例題 4-2] 簡支梁例題結構示意圖
表 4.2: [例題 4-2] 簡支梁案例拓樸最佳化結果
Initial Optimization Result
Shape
C (kN× mm) 1.725 3.072
V (%) 100.0 30.0
CV (kN× mm) 1.725 0.922
(a) 結構順從度收斂圖 (單位: kN× mm) (b) 體積收斂圖 (無單位)
圖 4.4: [例題 4-2] 簡支梁案例拓樸過程收斂圖
本例題之結果與懸臂梁例題相似，經過拓樸最佳化移除 70% 的體積後，雖 然結構順從度比初始形狀提升了 78%，然而結構順從度與體積之乘積則降低了約 47%，代表單位體積之結構勁度上升，顯示拓樸最佳化可以找出無效率元素的所 在位置並加以移除，使材料能夠分布於設計領域中最有效率處。
4.3.1 不同初始拓樸形狀
由於 ESO 與 BESO 方法皆是由前一步之分析結果來決定下一步之演進方向，
因此最佳化結果可能會受到初始拓樸形狀影響，而 BESO 又具有將移除之元素填 回的自由，因此初始拓樸形狀不一定要從全設計領域填滿硬材料開始，亦可以改 為從全設計領域填滿軟材料開始，以測試演算法是否能收斂至一至的解。然而初 始設計領域全填滿軟材料代表初始體積為 0，導致式 (2.18) 的更新體積方式將不 適用，因此需要對式 (2.18) 進行修改，使其亦可以適用於初始體積為 0 的情形，
修改方式如下:
Vi =
M ax(Vi−1(1− ER), Vtarget), when Vi−1 > Vtarget M in(1− (1 − Vi−1)(1− ER), Vtarget), when Vi−1 ≤ Vtarget
(4.1)
修改內容主要是針對前一次迭代體積小於目標體積的更新體積方式進行修 改，由原本將前一次迭代實心元素體積增加固定比例，作為下一次迭代實心元素 體積，改成將前一次迭代空心元素體積減少固定比例，作為下一次迭代空心元素 體積，再以總體積減去空心元素體積便可以得到實心元素體積，透過此修改方法 便可以適用於初始體積為 0 的情形。
另外，由於由軟到硬的拓樸過程中，空心元素體積會逐漸減少，因此以上所 採用之修改方法還有一個優點，也就是當體積逐漸接近目標體積而趨近於收斂 時，可以使每次迭代之體積改變量也會逐漸減少，避免趨近於收斂時出現過大之 改變，同時保持最初幾步迭代不至於演進過慢。以下將以例題演示不同初始拓樸 形狀對最佳化的過程與結果之影響。
[例題 4-3]
本例題之邊界條件、載重形式以及材料參數皆與例題 4-2 相同，並且將初始 設計領域改成全填滿軟材料，也就是將拓樸過程由硬到軟改成由軟到硬，以測試 因為不同初始拓樸形狀所產生的差異。
