碟形天線之流固耦合分析與最佳化

全文

(1)國立交通大學機械工程學系碩士論文. 碟形天線之流固耦合分析與最佳化. Fluid Structure Interaction Analysis and Optimization on Dish Antenna. 研究生：洪智偉指導教授：洪景華. 教授. 中華民國九十四年六月.

(2) 碟形天線之流固耦合分析與最佳化研究生：洪智偉. 指導教授：洪景華教授. 國立交通大學機械工程學系摘要. 設計一個新的碟形天線，在量產前為達到規格的要求，必需要先經過風洞測試，以確保其在自然的風吹雨打的環境中，不會產生永久的變形而影響其接收訊號的品質。傳統上，碟形天線的結構強度設計，譬如碟盤的厚度，必需使用試誤法來反覆的在風洞中測試，一旦結構強度不夠而測試失敗則必需重來，但若為了強化結構又使用過多的材料導致生產成本無法降低，經過多次的測試不但耗費金錢也耗費時間。本研究嘗試以流固耦合的方式，將計算流體力學分析軟體與有限元素應力分析軟體整合，對碟形天線中，碟盤以及支撐桿這兩部份做強度分析，並將此分析與風洞實驗進行驗証，最後再與最佳化軟體連結，對碟盤以及連結器做節省使用材料的最佳化設計。. I.

(3) Fluid Structure Interaction Analysis and Optimization On Dish antenna Student：ChihWei Hung. Adviser：Dr. Chinghua Hung. Department of Mechanical Engineering National Chiao Tung University. ABSTRACT. Before mass-producing a new type of dish antenna, this dish antenna must go through the survival test of wind tunnel to make sure it can work well under the natural wind condition. Traditionally, the design of the strength of the dish antenna uses the method of trial and error according to the wind tunnel tests. This procedure costs a lot of time and money. The main objective of this research is to utilize the numerical scheme of fluid-structure interaction, which integrates CFD and FEM, to analyze the strength of the dish and the connecting pipe. After verified by the results from wind tunnel test, this approach was further incorporated with optimization software to reduce the cost of the antenna dish and the connector.. II.

(4) 誌. 謝. 首先要感謝的是指導教授洪景華博士，其亦師亦友的身份，除了在學業上對學生的提攜，並且不吝嗇的分享生活上的知識與經驗，讓我在二年的研究生涯中，除了對機械這塊領域有更深一層的認識之外，也對自己未來的方向有進一步的思考。在二年研究所的路程中，感謝洪榮崇學長對我的照顧，還要感謝的有宇中、中興與政成學長，共同渡過兩年甘苦時光的琇晶與維德，已經畢業的建溢、雅雯和奇忠學長姊，以及嘉偉、銘傑與麒禎學弟，很高興在研究生涯能夠認識你們，也感謝清大動機的朋友們，在新竹陪伴我渡過了六年。同時也感謝幫助我論文中實驗部份以及模擬困難的朋友們，沒有你們，我不可能可以如期完成這份研究。最後我要感謝我的父母，在我的求學過程一直大力的支持，你們的關心，點滴在心頭，謝謝你們。. 謝謝一路陪伴我走來的所有人。. III.

(5) 目錄中文摘要………………………………………………………………..Ⅰ 英文摘要…………..………………………………………………….. Ⅱ 誌謝……………………………………………………………………..Ⅲ 目錄………………………………………………………………....… Ⅵ 表目錄……………..………………………………………………….. Ⅶ 圖目錄…………………………………………………………………..Ⅷ 第一章緒論...............................................................................................1 1.1 前言 ..................................................................................................1 1.2 研究動機與目標..............................................................................2 1.3 碟形天線的規格..............................................................................2 1.4 研究方法 ..........................................................................................3 1.5 文獻回顧 ..........................................................................................4 第二章計算流體力學分析與有限元素應力分析整合 ..........................9 2.1 前言 ..................................................................................................9 2.2 軟體簡介 ..........................................................................................9 2.2.1 前處理器I-DEAS簡介 .............................................................9 2.2.2 計算流體力學分析軟體STAR-CD簡介 ...............................11 2.2.3 有限元素分析軟體ABAQUS簡介........................................11 2.3 材料性質 ........................................................................................12 2.4 整合流程說明 ................................................................................12 2.4.1 整合流程 ................................................................................12 2.4.2 I-DEAS建立模型 ....................................................................13 2.4.3 STAR-CD之壓力計算.............................................................13 IV.

(6) 2.4.4 C++語言程式進行壓力轉檔 ..................................................14 2.4.5 ABAQUS之負荷輸入 .............................................................15 2.5 分析結果 ........................................................................................15 2.5.1 流固耦合分析結果 ................................................................15 2.5.2 其它討論－驗証假設條件與簡化模型 ................................16 第三章風洞實驗.....................................................................................18 3.1 實驗設備與操作環境....................................................................32 3.2 實驗規劃 ........................................................................................32 3.3 實驗結果與討論............................................................................33 第四章最佳化系統規劃 ........................................................................36 4.1 前言 ................................................................................................39 4.2 軟體簡介 ........................................................................................39 4.3 有限元素分析與最佳化整合流程................................................40 4.4 最佳化系統規劃............................................................................41 4.5 分析結果 ........................................................................................43 第五章結論.............................................................................................49 5.1 結論 ................................................................................................49 5.2 未來展望 ........................................................................................50 附錄A. ALE法進行流固耦合分析.......................................................54. A.1. 前言...........................................................................................54. A.2. ALE法之網格描述與座標系統...............................................54. A.3. ALE法分析情況簡介...............................................................55. A.4. 軟體簡介...................................................................................56 V.

(7) A.4.1. LS-DYNA.........................................................................56. A.4.2. 前處理器FEMB................................................................56. A.4.3. 後處理器LS-Pre/Post .......................................................57. A.5. 模型建立與邊界條件設定 ......................................................57. A.6. 分析結果與討論 ......................................................................58. 附錄 B 壓力轉換C++程式語言原始碼…………………………………..65 附錄 C MOST INPUT C 語言程式原始碼………..…...………………68 附錄 D MOST 資料輸入檔內容………..……………...………………73. VI.

(8) 表目錄表 2. 1. 碟盤的材料性質 .......................................................................18. 表 2. 2. 空氣的材料性質 .......................................................................18. 表 2. 3. 碟盤編號 1030 的節點各方向應變值………………………..19. 表 2. 4. 支撐桿編號 1680 的節點各方向應變值 .................................19. 表 2. 5. 比較有無接收器碟盤節點 1030 各方向應變值 .....................19. 表 3. 1. 碟盤上的應變規在風速 40m/s下所量測到的應變值 ............35. 表 3. 2. 支撐桿上的應變規在風速 40m/s下所量測到的應變值 ........35. 表 3. 3. 碟盤與支撐桿模擬分析結果與實驗結果比較 .......................35. 表 4. 1. 原始尺寸與最佳化結果比較表 ...............................................44. 表 4. 2. 目標函數－原始尺寸與最佳化結果比較表 ...........................44. 表A. 1. LS-DYNA碟盤及支撐桿上量測點各方向的應變值..............60. 表A. 2. LS-DYNA碟盤及支撐桿模擬與實驗的比較表......................60. VII.

(9) 圖目錄圖 1. 1. 歐洲地區電視分佈一覽 .........................................................6. 圖 1. 2. W公司－碟形天線AU2..........................................................7. 圖 1. 3. 碟形天線AU2 的背視圖.........................................................8. 圖 1. 4. 碟形天線AU2 的側視圖.........................................................8. 圖 2. 1. 拉伸試驗試片尺寸圖 ...........................................................20. 圖 2. 2. MTS萬能拉伸試驗機 ...........................................................20. 圖 2. 3. 應力應變曲線圖 ...................................................................21. 圖 2. 4. 應力應變曲線圖－線性範圍 ...............................................21. 圖 2. 5. 整合流程圖 ...........................................................................22. 圖 2. 6. I-DEAS輸出之碟盤模型 ......................................................23. 圖 2. 7. I-DEAS輸出之支撐桿模型 ..................................................24. 圖 2. 8. STAR-CD建立之模型...........................................................25. 圖 2. 9. C++言語程式搜尋壓力示意圖 ............................................26. 圖 2. 10. STAR-CD之壓力分佈結果－正面視角 ............................27. 圖 2. 11. STAR-CD之壓力分佈結果－斜視角.................................27. 圖 2. 12. ABAQUS碟盤von Mises等效應力分佈等高線圖 ............28. 圖 2. 13. ABAQUS支撐桿von Mises等效應力分佈等高線圖 ........29. 圖 2. 14. STAR-CD速度分佈圖＿傾斜角＝0∘（只有碟盤） ......30. 圖 2. 15. STAR-CD速度分佈圖＿傾斜角＝2.5∘（只有碟盤） ...30. 圖 2. 16. 傾斜角＝0∘的碟盤von Mises等效應力分佈等高線圖 ..31. 圖 2. 17 傾斜角＝2.5∘的碟盤von Mises等效應力分佈等高線圖 .31 圖 3. 1. 碟形天線在風洞中架設的背視圖 .......................................36. 圖 3. 2. 碟形天線在風洞中架設的側視圖 .......................................36. VIII.

(10) 圖 3. 3. 應變規在碟盤上的位置圖 ...................................................37. 圖 3. 4. 應變規在支撐桿上的位置圖 ...............................................38. 圖 4. 1. 連結器－實體圖與CAD圖...................................................45. 圖 4. 2. MOST流程和使用者介面 ....................................................46. 圖 4. 3. 有限元素分析與最佳化系統之整合流程作業 ...................47. 圖 4. 4. 最佳化－目標函數與疊代次數關係圖 ...............................48. 圖A. 1. LS-DYNA建立之模型與邊界條件設定..............................61. 圖A. 2. LS-DYNA碟盤von Mises等效應力分佈等高線圖 .............62. 圖A. 3. LS-DYNA碟盤x-y方向的應變分佈等高線圖 ....................62. 圖A. 4. LS-DYNA支撐桿的應變分佈等高線圖..............................63. 圖A. 5 STAR-CD速度流場分佈圖...................................................64 圖A. 6. LS-DYNA分析結果速度分佈圖..........................................64. IX.

(11) 第一章. 緒論. 1.1 前言科技日新月異，人們接受資訊的管道也愈來愈多樣化，從平面的報章雜誌演進到聲光效果十足的電視節目，單純的只有地方上的電視台節目，已經滿足不了追求娛樂、知識的現代人。直播衛星提供了更多節目選擇。根據市場統計，台灣的境外與境內衛星電視頻道供應商，在一百一十家以上，然而有線電視系統可以容納的頻道數，約七十個左右，直播衛星所提供的視訊服務，在兩百個頻道以上。[1] 直播衛星利用距離地球表面 3 萬 6 千公里的同步衛星，使用 Ku 頻或 C 頻的頻率，進行訊號的傳送，由於傳送的是數位訊號，因此直播衛星擁有絕對清晰的數位畫面，與幾近完美的立體音效，直播衛星的收視戶利用碟形天線接收訊號，可以收看到比有線電視系統所提供更好的影音品質。另外有線電視必須透過線纜，才能進行訊號傳輸，因此偏僻的地區居民，常因為有線電視線纜沒有在當地舖設，而無法訂閱有線電視。直播衛星以無線方式播出訊號，訂戶只需購買 18 吋大小的碟形天線，就可接收衛星電視訊號，這對於偏遠地區的居民來說，確實是較有線電視更好的選擇。目前日本、美國、加拿大和澳洲為主要推廣的數位衛星電視國家。在日本藉由衛星傳遞數位訊號，以克服日本山脈綿延的地形；在美國、加拿大和澳洲，由於地廣人稀的先天條件，使的有線電視的發展在環境上就受到了限制，但對於直播衛星來說這是一個很好的發展優勢。 1.

(12) 因為直播衛星有著節目品質豐富、傳遞訊號不受地域限制等優點，以歐洲來說，依據 SES ASTRA Satellite Monitor 的調查，歐洲的電視用戶目前僅有 21.4%是使用衛星電視（如圖 1.1），所以對直播衛星市場抱持著樂觀的人士認為，未來直播衛星的佔有的比率會逐漸升高。因此，如何針對收看直播衛星節目的碟形天線，有效的降低其生產成本，是一個非常值得研究的課題。. 1.2 研究動機與目標設計一個新的碟形天線，在量產前為了達到規範的要求，必需要先經過風洞測試，以確保其在自然的風吹的環境中，不會產生永久的變形而影響其接收訊號的品質。傳統上，碟形天線的結構強度設計，譬如碟盤的厚度，必需使用試誤法來反覆的在風洞中測試，一旦結構強度不夠而測試失敗則必需重來，但若為了強化結構又使用過多的材料導致生產成本無法降低，經過多次的測試不但耗費金錢也耗費時間。本研究的目標，在於利用現有的分析模擬軟體，建立出碟形天線受風壓產生影響的模擬分析流程，將得到的模擬分析結果提供設計者作為參考，以期望減少碟形天線在量產前時間與金錢的花費，並降低生產的成本同時滿足強度剛性的要求。. 1.3 碟形天線的規格一般碟形天線在販售時需註明其在多少風速的操作下，可正常接收訊號、不產生永久變形、不會損壞分解。通常來說此規格的制定是依據下列幾個因素：. 2.

(13) (1). 碟盤的大小。. (2). 接收的頻率（Ku 頻或 C 頻）。. (3). 產品的壽命與售價（公司的政策）。. (4). 用於接收或是發射信號。如果用於發射時，設計上所能承受之風速將越高，各個風速之定義內容也將不同。. (5). 產品視各個地區的不同需求去做改變。. 本研究以園區某 W 公司銷往美國的碟形天線 AU2（如圖 1.2、圖 1.3、圖 1.4，各部位的名稱如圖上所標示）作為模擬分析的目標，該公司根據美國地區的需求訂出 AU2 的規格為：在風速 20m/s 的環境下可正常接收訊號，在風速 40m/s 下天線不會產生永久的變形，在風速 60m/s 不會損壞分解。. 1.4 研究方法為了分析碟形天線是否產生了永久變形，以確定天線在經過風吹之後還是能夠正常接收訊號，因此本研究的模擬目標先設定為，在風速 40m/s 環境下的碟形天線，碟盤與支撐桿的應力應變分析。本研究嘗試以流固耦合（Fluid Structure Interaction；FSI）的方式，將計算流體力學（Computational Fluid Dynamics；CFD）分析軟體與有限元素（Finite Element Method；FEM）應力分析軟體整合。此方法的先前假設條件為不考慮碟形天線的盤面在受到風壓之後的變形，是否會影響流場的改變；因此將對碟形天線的分析的過程分成二個部分：第一個部份是利用計算流體力學分析軟體，分析出碟形天線的固定盤面在風洞試驗中碟盤. 3.

(14) 前後的壓力差；第二部份將流場分析結果，轉至相對應的有限元素分析中的元素上來做衛星碟盤結構的應力應變分析，接著計算出碟盤與前方接收器受到風壓所產生的受力，並將此受力施加給下方支撐桿的部份做應力應變分析。經由風洞實驗，量得碟盤與支撐桿上的應變值，與模擬分析的結果做校正後，再進行最佳化的流程；最佳化所針對的目標函數為碟形天線的材料成本，而成本的計算則包含了碟盤的體積、連接器的大小，而設計變數為碟盤的厚度以及碟盤上固定點的位置。另外在附錄一裡，本研究也曾嘗試過另一種流固耦合的方法，利用 Arbitrary Lagrangian-Eulerian（ALE）方法。此方法最大好處在於只需要建立單一模型、使用單一求解核心，就可以求得碟盤上的應力應變分析結果，並不需要進行兩種軟體的整合。但由於模擬結果與風洞實驗結果比較誤差過大，故僅放在附錄 A 中以供參考。. 1.5 文獻回顧關於整合計算流體力學分析軟體與有限元素應力分析軟體的方法，過去已有許多人使用來應用在不同的問題上。如愛爾蘭的 Penrose 等[2]，整合了商業軟體 CFX-5 和 ANSYS 來分析動脈血管，在實驗的驗証下証明此方法有良好的分析結果；而瑞士的 Filsinger 等[3]，整合了 ABB Turbo Systems, Ltd. 所開發的計算流體力學軟體與 ABAQUS，針對渦輪增壓器的振動行為做分析，可得到渦輪葉片上的負載分佈結果。. 4.

(15) 關於碟形天線的研究雖多，但絕大部份都是著重於通信、接收訊號方面的問題，至於對其結構受風壓分析的研究，目前還少有人提出相關討論。因此本研究希望能建立起一套碟形天線在受到風壓而產生影響的分析流程，提供給碟形天線設計者作參考，以期望減少時間與金錢的花費，有效的降低生產的成本。. 5.

(16) 圖 1. 1. 歐洲地區電視分佈一覽. （資料來源：SES ASTRA Satellite Monitor，2004/6）. 6.

(17) 碟盤. 碟盤上的四個固定點. 接收器. 支撐桿. 圖 1. 2. W 公司－碟形天線 AU2. 7.

(18) ︵連結器︶連結碟盤上四個固定點與支撐桿. 圖 1. 3. 碟形天線 AU2 的背視圖. 圖 1.4. 碟形天線 AU2 的側視圖. 8.

(19) 第二章. 計算流體力學分析與有限元素應力分析整合. 2.1 前言此方法的假設條件為，並不考慮碟形天線在受到風壓之後所產生的變形，是否會影響流場的改變，因此可以將此分析拆成兩個部份：計算流體力學分析與有限元素應力分析。目前市面上在處理計算流體力學分析方面比較有名的軟體有 FLOW-3D、STAR-CD、FLUENT 等，而有限元素分析軟體有 ABAQUS、ANSYS、LS-DYNA 等。在本研究中以 I-DEAS 作為前處理器，STAR-CD 作為計算流體力學分析的軟體， ABAQUS 作有限元素分析軟體，並以自行撰寫的 C++語言程式作為中間轉換壓力分佈的介面。. 2.2 軟體簡介 2.2.1. 前處理器 I-DEAS 簡介. 電腦數值分析模擬的第一步驟便是建立正確的幾何模型。在絕大多數的工程問題上，幾何形狀複雜且不規則，所以需要一套強大的前處理器（Preprocessor）輔助建立模型的外觀及分割網格。本研究採用專業前處理器軟體 I-DEAS（Integrated Design Engineering Analysis Software）來建立幾何模型，因為 I-DEAS 與其他前處理器比較，具以下優點： (1) 幾何外型建立：參數化的繪圖方式，便於建立外型及修改，適合用來建立有限元素分析模型。 (2) 親切的使用者介面：具圖形界面的常用繪圖指令，讓 9.

(20) 使用者更容易上手，並有選單式功能表調整內部參數。 (3) 優秀的網格建立功能：不但具有自動建立網格（Free mesh）的功能，也能在簡單的模型上建立規則的網格（Mapped mesh），此外，在建立自動網格時可作局部加密，建立規則網格時也可以權重的方式調整密度。 (4) 定義邊界方便：I-DEAS 除了可以在元素或節點上給定限制或負載外，亦可將之施加於幾何外型上（Geometry based），例如對一平面施加壓力或限制條件。利用這種方式亦可以將我們需要平面上的元素篩選出。 (5) 與其他有限元素分析軟體介面的整合性高：可將建立完成的模型輸出成其他有限元素分析軟體的輸入檔，並可載入其他有限元素軟體產生的檔案。. I-DEAS 軟體是由多個應用程式模組（Application）所組成，這些應用程式的模組包括了設計 (Design) ､繪圖（ Drafting ）､模擬（ Simulation ）､測試（ Test ）､製造（Manufacturing）､管理（Management）和幾何模型轉換（Geometry translators）等，是一套完整且功能大的電腦輔助設計､分析與製造的軟體。本研究主要將使用設計與模擬兩個模組，建立好模型之後，分別輸出 IGES 檔，由 STAR-CD 進行流場分析，以及輸. 10.

(21) 出 ABAQUS 輸入檔，加上 STAR-CD 分析得到的壓力分佈，最後利用 ABAQUS/STANDARD 進行有限元素分析。. 2.2.2. 計算流體力學分析軟體 STAR-CD 簡介. STAR-CD 為一泛用型計算流體力學軟體，針對流場結構提供了精確的分析並作定性及定量上的探討，並針對特殊問題提供各式模型方便使用者設定，如旋轉機械，燃燒行為，兩相傳輸等。並可進行暫態，靜態的分析以及移動網格等設定。 STAR-CD 亦為使用者提供了一標準使用流程 STARGIDE，從導入 CAD-產生網格-設定流體性質-求解方式 -後處理結果輸出，使用者只要依循此設定流程即可輕易完成分析，縮短了摸索軟體的時間。. 2.2.3. 有限元素分析軟體 ABAQUS 簡介. ABAQUS 是由 ABAQUS 公司發展，是一套功能強大的工程模擬分析軟體，無論是在處理線性（Linear）或非線性（Nonlinear）問題，都有優秀的表現。使用者只需提供工程資料，如結構幾何、材料行為、邊界條件及施予的負載， ABAQUS 會自動連續的調整負載增量（Increment）與收斂公差，使模擬結果更為精確，更有效率。針對本研究所使用的模組如下： (1). ABAQUS/STANDARD：標準有限元素求解主程式，處理一般線性或非線性問題。. 11.

(22) (2). ABAQUS/CAE：ABAQUS 的視窗前後處理軟體。由於前處理器建立幾何模型功能不理想，本研究只採用其後處理功能，觀看應力、應變等數值分析結果。. 2.3 材料性質在建立模型前必需先確定材料係數。材料由 W 公司所提供，為該公司所生產碟形天線 AU2 的碟盤以及其下方支撐桿的材料。材料的各項係數則利用拉伸試驗取得。拉伸試驗由交大機械複合材料結構實驗室協助完成。試片的規範參考美國 ASTM 標準金屬拉伸試驗 E8M（試片尺寸如圖 2.1）；拉伸試驗在 MTS 萬能拉伸試驗機上（如圖 2.2）進行。實驗的結果整理成應力應變曲線圖（如圖 2.3～圖 2.4），可得到楊氏係數（Young’s modulus）、抗拉強度（Tensile strength）和降伏強度（Yield stress），加上一般材料手冊上所查到鋼材的浦松比（Poisson’s ratio），所得到的材料性質列於表 2.1。. 2.4 整合流程說明 2.4.1. 整合流程此方法的假設條件為，並不考慮碟形天線在受到風壓之. 後的變形，是否會影響流場的改變。因此將對碟形天線的分析的過程分成計算流體力學分析與有限元素分析兩部份。模型的建立是以 W 公司所生產的碟形天線 AU2 為參考樣本。整合的流程如圖 2.5 所示，而流程詳細的文字說明將依序在接下來的章節作介紹。. 12.

(23) 2.4.2. I-DEAS 建立模型使用 I-DEAS Design Application 建立出碟形天線的碟. (1). 盤、前方接收器以及支撐桿部份的近似模型，輸出 IGES 檔，供 STAR-CD 滙入。使用 I-DEAS Simulation Application 功能，對碟盤建立. (2). 網格、節點、元素和設定邊界條件（如圖 2.6），輸出 ABAQUS 輸入（Input）檔。由於盤面的尺寸與厚度之比大於 20，這種情況下適合選擇使用殼元素（Shell element）來建立網格。之後可將 STAR-CD 所分析的盤碟上壓力值輸入後，再進行碟盤部份的有限元素應力分析。滙入由 W 公司所提供支撐桿的 IGES 檔，使用 I-DEAS. (3). Simulation Application 功能建立網格、節點、實體元素和設定邊界條件（如圖 2.7），其中支撐桿上端部份的元素設為殼元素並設為剛體（Rigid body），為連結器（連接碟盤中四個固定點以及支撐桿的部份）簡化之模型，以作為之後承接碟盤受力用途。最後輸出支撐桿 ABAQUS 輸入檔。分析完碟盤的有限元素應力分析後，可得到碟盤所受到的和力，將此受力與前方接收器的受力，一併施加至支撐桿的剛體元素上，最後進行有限元素應力分析。. 2.4.3. STAR-CD 之壓力計算. 此部份的模擬由交大機械熱流應用實驗室完成。模擬的. 13.

(24) 環境為天線在風洞中，以正面承受 40m/s 的風速，而分析完成後所需的資料為風在碟盤前後產生的壓力差。 (1). 將入風口的速度設為 40m/s。. (2). 模擬的分析過程假設為均勻穩態狀況；紊流模式訂為 k-ε模式。. (3). 將 IGES 檔匯入的模型設為固體，並建立流體網格。. (4). 為了保持風速到盤面之前還能維持 40m/s，除了入風口、出風口及地面，其餘壁面皆設為無摩擦的性質。. (5). 工作流體（空氣，Air）的性質，如表 2.2 所示。. (6). 設定完成的模型的側視圖如圖 2.8 所示。. (7). 分析完成後，將盤面前、後的流體網格壓力值輸出，得到的資料格式為（x，y，z，P）；其中 x、y、 z 為該網格的中心座標值，P 則為該網格的壓力值。. 2.4.4. C++語言程式進行壓力轉檔以 C++語言程式進行以下各步驟的工作，原始程式碼如附錄二。 (1). 讀取 STAR-CD 輸出結果，包含盤面前、後網格的中心座標值，以及該網格的壓力值。. (2). 讀取 ABAQUS 輸入檔裡所有殼元素的節點座標值；並計算出所有殼元素的中心座標值。. (3). 設定一適當的搜尋半徑，以該元素的中心為圓心，半徑內所有網格上壓力的平均值給定為該元. 14.

(25) 素的壓力值。圖 2.9 為示意圖。 (4). 2.4.5. 輸出所有有限元素上的壓力值。. ABAQUS 之負荷輸入 (1). 將 C++ 語言程式輸出的壓力值輸入至碟盤的 ABAQUS 輸入檔。. (2). 給定材料係數後，啟動 ABAQUS/STANDARD 分析求解，以 ABAQUS/CAE 觀看碟盤部份分析完成後的結果。. (3). 計算出碟盤受到風壓產生的受力，將此受力與前方接收器之受力一併施加到支撐桿的剛體元素上，以 ABAQUS/STANDARD 模組分析求解，再以 ABAQUS/CAE 觀看支撐桿部份分析完成後的結果。. 2.5 分析結果 2.5.1. 流固耦合分析結果. 圖 2.10、圖 2.11 為 STAR-CD 分析結果的壓力分佈圖，可以看出風速 40m/s 所產生的壓力主要是在碟盤面上以及前方接收器的正面上。從碟盤的 von Mises 等效應力分佈等高線圖（如圖 2.12）來看，其中紅色圈起來的為邊界條件設定為固定的元素，可觀察出應力集中的部份為中間固定點的外側兩邊。圖 2.12 中箭頭所指編號 1030 的節點，為盤面上應力集中. 15.

(26) 的節點，其各方向的應變值如表 2.3 所示。碟盤的受力和為 357.9 N，而由圖 2.10 可估算出前方接收器在 z 方向上的受力大約為 26 N，將此二受力施加到支撐桿上。分析完成後，從支撐桿的 von Mises 應力分佈等高線圖（如圖 2.13），可以得到在箭頭所指編號 1680 的節點，是支撐桿上除邊界限制外應力最大值的節點，其各個方向的應變值如表 2.4 所示。在下一章中，將以本次模擬分析的結果與風洞實驗的結果來進行比較，以驗証並校正模擬模式的準確性。. 2.5.2. 其它討論－驗証假設條件與簡化模型. 在前言中有提到，此方法的假設條件為：不考慮碟形天線在受到風壓之後的變形，是否會影響流場的改變。驗証的方法為：假設天線受到風壓而產生變形，使得原本面向正前方的碟盤造成了一個 2.5∘的傾斜角，此傾斜角是假設支撐桿上固定螺絲的滑移與支撐桿本身的彈性變形所造成的傾角。從 STAR-CD 的速度分析圖（如圖 2.14、2.15）來看，傾斜角對流場有造成些微的影響，但差異不大。將兩者不同傾斜角所產生的壓力差資料輸入到 ABAQUS 做應力分析後，從傾斜角為 0∘和 2.5∘的碟盤碟盤 von Mises 等效應力分佈等高線圖（如圖 2.16 和圖 2.17），可以看出兩者的應力分佈趨勢與數值大小皆非常的接近，因此“不考慮碟形天線在受. 16.

(27) 到風壓之後的變形，是否會影響流場的改變＂的這個假設，是合理的。合理的減化模型可以達到節省模擬分析時間的目的，因此嘗試將前方接受器與下方支撐桿忽略，模型只留下碟盤部份，觀察這樣的流場對碟盤應力應變分析的影響。表 2.5 為比較流場分析中，是否有考慮前方接收器的影響，碟盤上編號 1030 的節點各方向的應變值。由表可以看出兩者分析的結果相差不多，因此可以考慮在後續的碟形天線研究中，有關流場分析的部份，為了達到節省分析時間的目的，可忽略接收器的影響，而僅分析碟盤的部份。. 17.

(28) 表 2. 1. 碟盤的材料性質. 係數. 對應值. 單位. 楊氏係數（Young's modulus）. 202000. MPa. 浦松比（Poisson's ratio）. 0.29. 降伏強度（Yield stress）. 160. MPa. 抗拉強度（Tensile strength）. 290. MPa. 表 2. 2. 空氣的材料性質. 係數. 對應值. 單位. 密度（Density）. 1.205. Kg/m3. 分子黏滯係數（Molecular Viscosity） 1.81*10-5 N-s/m2. 18.

(29) 表 2. 3. 碟盤編號 1030 的節點各方向應變值. Node 1030. E11 -86. E22 565. E12 339 （單位：micro-strain）. 表 2. 4. Node 1680. 支撐桿編號 1680 的節點各方向應變值. E11 -113. E22 62. E33 454. E12 9. E13 -41. E23 634. （單位：micro-strain）. 表 2. 5. 比較有無接收器碟盤節點 1030 各方向應變值 Node 1030 有接收器無接收器相差. E11. E22. E12. -86. 565. 339. -99. 602. 362. 15%. 6.50%. 6.80%. （單位：micro-strain）. 19.

(30) 單位：mm. 圖 2. 1. 拉伸試驗試片尺寸圖. 圖 2. 2. MTS 萬能拉伸試驗機 20.

(31) 應力應變曲線圖應力（MPa） 350 300 250 200 150 100 50 0 0.0. 0.1. 0.2. 0.3. 0.4. 0.5. 0.0002. 0.00025. 應變. 圖 2. 3. 應力應變曲線圖. 應力應變曲線圖＿線性範圍應力（GPa） 0.045 y = 202.52x + 0.0011. 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0 0. 0.00005. 0.0001. 0.00015 應變. 圖 2. 4. 應力應變曲線圖－線性範圍. 21.

(32) Original CAD File Import STAR-CD. I-DEAS to create model For Dish. Data of Pressure C++ Transfer. Dish’s ABAQUS Input file For Pipe ABAQUS/Standard Total Force Of Dish ABAQUS Output. Pipe’s ABAQUS Input file. ABAQUS/Standard. ABAQUS Output. 圖 2. 5. 整合流程圖. 22.

(33) 圖 2. 6. I-DEAS 輸出之碟盤模型. 23.

(34) 圖 2. 7. I-DEAS 輸出之支撐桿模型註：深色網格為剛體元素. 24.

(35) 圖 2. 8. STAR-CD 建立之模型. 25.

(36) P為流體網格壓力. 圖 2. 9. C++言語程式搜尋壓力示意圖. 26.

(37) 單位：Pa. 圖 2. 10. 圖 2. 11. STAR-CD 之壓力分佈結果－正面視角. STAR-CD 之壓力分佈結果－斜視角 27.

(38) 單位：MPa. Node. 1030. 圖 2. 12. ABAQUS 碟盤 von Mises 等效應力分佈等高線圖. 28.

(39) Node. 1680. 圖 2. 13. ABAQUS 支撐桿 von Mises 等效應力分佈等高線圖. 29.

(40) 圖 2. 14. 圖 2. 15. STAR-CD 速度分佈圖＿傾斜角＝0∘（只有碟盤）. STAR-CD 速度分佈圖＿傾斜角＝2.5∘（只有碟盤）. 30.

(41) 圖 2. 16. 傾斜角＝0∘的碟盤 von Mises 等效應力分佈等高線圖. 圖 2. 17 傾斜角＝2.5∘的碟盤 von Mises 等效應力分佈等高線圖 31.

(42) 第三章. 風洞實驗. 3.1 實驗設備與操作環境風洞實驗係由委託中山科學研究院－航空研究所的低速風洞組來進行實驗。目的在與數值分析結果進行比較，以驗証並校正模擬模式的準確性。碟形天線架設在風洞中的情況如圖 3.1、圖 3.2 所示；其中圖 3.2 左下方的顯示數字為當時風洞中的風速。. 3.2 實驗規劃為了驗証第二章中所模擬出來的結果，選擇量測碟形天線受風壓產生的應變來做比較。參考先前模擬分析的結果，由圖 2.12 碟盤的 von Mises 等效應力分佈等高線圖，可以看出在中間固定點的外側 45°的地方有最大的應力值；由圖 2.13 支撐桿的 von Mises 等效應力分佈等高線圖來看，雖然最大的應力值產生在固定的邊界元素上，但該處並不容易貼上應變規，故選擇貼在箭頭處會較適當。因此實驗時，應變規所貼的位置，在碟盤的部份，應變規貼在盤面中間固定點的斜外側，且其量測應變值的方向為 45°（如圖 3.3）；在支撐桿的部份，貼在靠近固定座的桿子中間部位（如圖 3.4）。另外應變規配置的方式採用全橋式，這種配置可提供最高的敏感度及最少的誤差，因為全橋式配置提供最大的輸出，雜訊的影響相對也較小。實驗開始的十分鐘，風洞的風速由 0m/s 慢慢提升到 40m/s， 32.

(43) 等風速到達 40m/s 後五分鐘，才開始讀取應變值。. 3.3 實驗結果與討論盤面上左右各貼一個應變規，支撐桿上貼一個，實驗共進行兩次，因此可得到四個盤面和二個支撐桿的數據，如表 3.1、表 3.2 所示。由碟盤應變規量測的結果，表 3.1（單位為 micro-strain）來看，數據 2 為 395 與數據 4 為 660 與其它二個數據 580、435 相比，有明顯的差距；支撐桿上應變規量測二次的結果為 1156、 1045，而這二次實驗的結果差異性不大。為了將模擬分析的結果與實驗結果做比較，必需先算出編號 1030 節點的 45°方向應變值。將表 2.3 的結果代入公式： ε x ' = ε x l 2 + ε y m 2 + ε z n 2 + γ xy lm + γ yz mn + γ zx nl. （2.1）. 其中ε、γ為各座標方向的應變值； l、m、n 為新座標與原座標的夾角 cosine 值，以碟盤為例： l = m =1. 2 ，n=0；. 表 3.3 為模擬分析結果與實驗結果的比較，在碟盤部份，模擬結果與實驗結果最高值與最低值的誤差在-38%～3.5%，可以看出實驗時量測碟盤上應變值時的差異性很大，但模擬分析的結果為 409micro-strain，仍落於實驗結果 395～660micro-strain 的區間；支撐桿的部份，量測二次的結果非常接近，模擬結果與實驗結果的最高值與最低值的誤差在-50.3%～-45%。而模擬產生誤差的原因可能有下面幾個： (1). 建立的碟盤模型為近似的幾何形狀，因此碟盤上的應. 33.

(44) 力應變分析和真實上會有所出入，連帶影響計算施加到支撐桿上的受力。 (2). 簡化模擬的結果。在流場分析方面，只考慮碟盤受風壓的影響，而忽略了其它部份的影響；應力分析時是將碟盤固定在空間中，完全沒有位移，但實際上碟盤應該是與支撐桿連結，支撐桿的變形會造成碟盤的位移。. (3). 在風洞啟動後，支撐桿上的固定螺絲因為風吹而造成滑移的現像，使得碟盤的重心改變，也影響了在支撐桿上應變規所量取到的數值。因為在停風後有發現螺絲滑移的現像，使得應變規讀取到的數值沒有歸零，表示在實驗時所讀取到的數值已有了誤差。. (4). 實驗時風洞中的風速，是否有提供穩定的流速，也是造成模擬與實驗值有差異的可能原因。. (5). 由碟盤上應變規所量取的四個數值，其變動幅度很大，推測原因是碟盤在受測時，盤面的角度沒有完全校正到正面風向的位置，因此造成盤面上兩個應變規量測值的差異性。. 34.

(45) 表 3. 1. 碟盤上的應變規在風速 40m/s 下所量測到的應變值 1. 580. 2. 395. 3. 435. 4. 660 單位：micro-strain. 表 3. 2. 支撐桿上的應變規在風速 40m/s 下所量測到的應變值 1. 1156. 2. 1045 單位：micro-strain. 表 3. 3. 碟盤與支撐桿模擬分析結果與實驗結果比較碟盤. 支撐桿. 409. 575. 517.5. 1100.5. -20.9%. -47.8%. 660. 1156. -38%. -50.3%. 395. 1045. 3.5%. -45%. 模擬實驗平均值與平均值的誤差實驗最高值與最高值的誤差實驗最低值與最低值的誤差. 35.

(46) 圖 3. 1. 碟形天線在風洞中架設的背視圖. 圖 3. 2. 碟形天線在風洞中架設的側視圖. 36.

(47) 圖 3. 3. 應變規在碟盤上的位置圖 37.

(48) 圖 3. 4. 應變規在支撐桿上的位置圖. 38.

(49) 第四章. 最佳化系統規劃. 4.1 前言在確認了碟形天線流固耦合分析的準確性之後，接下來要說明的碟形天線的最佳化系統規劃；因此在本章中，要利用“減少材料的成本＂來達到降低生產成本這個最佳化的目的。由於支撐桿的分析結果與實驗結果誤差過大，因此在本章中，最佳化的分析僅考慮碟盤的部份。最佳化的目標函數中，其中一項是利用減少碟盤厚度來減少其本身的體積，另一項則藉由調整固定點的位置來縮小連結器（如圖 4.1）的體積。最佳化的方法則是利用連接最佳化軟體 MOST 與有限元素分析軟體 ABAQUS 來完成最佳化的目標。. 4.2 軟體簡介 MOST（Multifunctional Optimization System Tool ）為本校機械系曾錦煥教授以 C 語言所發展出來的一套最佳化軟體，被用來解連續和不連續變數的多目標最佳化問題[4]。因為本研究同時使用到連續以及不連續的設計變數，所以採用了 MOST。 (1) 連續變數：這模組使用 Sequential Quadratic Programming （SQP）method 作為單目標最佳化處理器（Single Objective Optimizer ），處理連續設計變數（ Continuous Design Variables）問題。這模組已經由許多最佳化問題驗証，MOST 可處理大量的最佳化問題，卻仍然具有優越的收斂性。[5] (2) 不連續變數：Branch and Bound Method（BBM）[6]是廣為所知並使用的不連續變數最佳化問題解決法。這套演算法將 39.

(50) 轉換不連續設計變數空間（Discontinuous Design Space）為連續，而在連續設計變數空間（Continuous Design Space）中再以上述的 SQP 模組求得最佳設計值。 (3) 多目標函數（Multi-objective）：在工程應用上常存在多目標最佳化問題，而多目標函數正好可解決這類問題。MOST 提供 Compromise Programming 、 Goal Programming 和 Surrogate Worth Trade-off method 等提供決策者找到最折衷答案。下面五個副程式（Subroutine）是 MOST 主要項目[4]： (1) cusermf / USERMF：目標函數計算副程式。 (2) cusercf / USERCF：限制方程式（Constraint Functions）計算副程式。 (3) cusermg / USERMG：目標函數微分方程計算副程式。 (4) cusercg / USERCG：限制方程式微分方程計算副程式。 (5) cuserou / USEROU：提供使用者依使用者需求格式輸出資料。圖 4.2 為 MOST 流程和使用者介面說明。. 4.3 有限元素分析與最佳化整合流程整合的流程如圖 4.2 所示，主要的想法為利用反覆修改有限元素分析輸入檔的內容，來進行最佳化的動作。詳細的說明如下： (1) 準備好先前使用的碟盤 ABAQUS 輸入檔。 (2) 進行有限元素分析，並產生分析結果輸出檔。. 40.

(51) (3) 讀取分析結果輸出檔中最佳化所需的資訊（如應力值等），進行最佳化的步驟。 (4) 由 MOST 核對目標函數，若收斂則為最佳化數值，並停止程式；若不收斂則執行下一步驟。 (5) MOST 計算獲得新的設計變數，並更新原先的 ABAQUS 輸入檔（如厚度、節點位置等），重複步驟(2)。. 4.4 最佳化系統規劃最佳化的目標函數為降低材料的成本，其中考慮到兩個項目，其中一項為減少碟盤的材料，而另一項為減少連結器的材料。由於碟盤的曲面是經由設計的特殊曲面，因此能夠針對硬盤部份，能夠改變的設計變數就只有碟盤的厚度tdish，因此碟盤的體積為碟盤面積Adish × tdish；連結器的部份如圖 4.1 所示，假設調整固定點後所會變動的只有虛線所圍起來的部份，其面積為A，並且與R（連結器中心點到固定點的長度）所圍成的圓形面積AR為正比關係。因此若調整固定點位置後新的R”所圍成的面積為AR”，則連結器調整後的面積為A × ( AR” / AR )，因此連結器虛線部份的體積等於A × ( AR” / AR ) × tc，tc為連結器的厚度。另外由 W 公司所提供的資料，碟盤部份所使用的板材的厚度有幾個規格可供選擇，分別為 0.4mm、0.5mm，以上厚度的價錢為每公斤 32 元，以及 0.6mm、0.65mm、0.7mm、0.8mm，這幾個厚度的價錢為每公斤 30 元；而連結器的所使用的板材厚. 41.

(52) 度固定為 1.5mm，其價錢為每公斤 28 元。此外在嘗試減少碟盤厚度，調整固定點位置的同時，以限制條件要求碟盤在承受 40m/s 風速的情況時，不致於發生碟盤超過降伏應力的情形而導致了塑性變形。由於碟盤是成對稱形狀，因此在控制固定點位置時，只需考慮一個固定點的位置。在此目的下，定義下列條件：目標函數：f = Min { Adish × tdish × D × P1 + A × ( AR” / AR ) × tc × D × P2}；限制條件： σ ele ≤ σ Y / 1.5 ；設計變數：非連續變數：碟盤厚度tdish（單位：mm）；連續變數：固定點座標 (x , y ) （單位：mm）。其中： Adish = 碟盤的面積（單位：mm2）； tdish = 碟盤的厚度（單位：mm）； D = 板材的材料密度（單位：kg/mm3）； P1 = 碟盤板材的價格（單位：元）； A = 圖 4.1 虛線所指的部份的面積（單位：mm2）； AR = 連結器中心點到固定點的長度所圍成的圓形面積（單位：mm2）； AR” = 調整後的長度所圍成的圓形面積（單位：mm2）；. tc = 連結器的厚度（單位：mm）； P2 = 連結器板材的價格（單位：元）；. 42.

(53) σ ele = 碟盤上所有元素的有效應力值（單位：Mpa）； σ Y = 碟盤的降伏應力值（單位：Mpa）；. 1.5 = 安全係數. 4.5 分析結果最佳化的結果如圖 4.4 所示，原始尺寸的目標函數為 39.29 元，而經過最佳化流程之後得到的最佳化數值為 32.37 元。整理最佳化結果於表 4.1、表 4.2。對碟盤厚度來說，由原來的 0.7mm 變成為 0.6mm，碟盤部份的成本減少了約 14.28 %；而連結器部份，調整固定點之後，原來的面積等比例縮小了 21.01%，因此成本同樣的也降低了 21.01%。調整固定點後，且厚度為 0.6mm 的碟盤進行有限元素應力分析後，由輸出結果檔可以得知，碟盤上的最大應力值為 105.8Mpa，乘上安全係數 1.5 後為 158.7Mpa，並沒有超過降伏應力 160Mpa，因此此最佳化結果是可以作為參考的。 W 公司所生產的碟形天線 AU2 型，在去年銷往美國的數量達百萬台以上，因此在量產 AU2 之前，若能配合本研究中所提出的最佳化流程，預估此一型的碟型天線可節省的材料成本約新台幣七百萬元。. 43.

(54) 表 4. 1. 原始尺寸與最佳化結果比較表. 原始尺寸最佳化結果. 厚度. x. y. 46. -30. 0.7. 34.18. -40.29. 0.6 單位：mm. 表 4. 2. 目標函數－原始尺寸與最佳化結果比較表目標函數碟盤部份連結器部份原始尺寸. 35.34 元. 5.6 元. 最佳化結果. 30.29 元. 4.42 元. 14%. 21%. 減少比例. 44.

(55) R. 圖 4. 1. 連結器－實體圖與 CAD 圖 45.

(56) INPUT DATA. MOST INP.C Read input data Multiobjective MOOD.C Multi-objective Optimizer. User Interaction. BRAND.C Noncontinuous Optimizer. FIXMF FIXMG FIXCF FIXCG. SQP.C Single Objective Continuous. xxxxMF xxxxMG xxxxCF xxxxCG. User Supplied Subroutines Describing a Specific Optimization Problem cusermf cusercf. Multi-objective. Have fix-DV. cusermg. Design Variable Optimizer. cusercg No fix-DV. Single-objective. cuserou. OUTPUT. xxxx : GOAL, COMP, SWT, OBJ. 圖 4. 2. MOST 流程和使用者介面. 46.

(57) 圖 4. 3. 有限元素分析與最佳化系統之整合流程作業. 47.

(58) 單位：新台幣. 40.00 39.29. 39.00 38.00 37.00 36.00 35.00 34.00. 32.37. 33.00 32.00 31.00 30.00 1. 2. 圖 4. 4. 3. 4. 5. 6. 7. 8. 9. 10. 11. 最佳化－目標函數與疊代次數關係圖. 48. 12. 13.

(59) 第五章. 結論. 5.1 結論在台灣，碟形天線收看衛星電視節目的市場一直沒辦法拓展開，但這並不代表碟形天線並沒有國際的市場。以 W 公司所銷往美國的碟形天線 AU2 為例，在 2004 年的出貨量達百萬台以上，因此如何有效的降低生產成本就顯得更為重要。本研究提出以流固耦合的方法，針對碟形天線在受到一定風速下的情況做應力應變分析，以自行撰寫的程式將計算流體力學分析軟體所分析出的壓力結果，轉換至有限元素分析軟體。而分析結果在與風洞實驗結果比較後，碟盤的應變值為 409micro-strain，與實驗結果的誤差在-38%～3.5%；而支撐桿的應變值為 575micro-strain，與實驗結果的誤差在-50.3%～-45%。但由於支撐桿部份誤差過大，後續最佳化分析只考慮與碟盤相關的部份。在確認此流固耦合分析的準確性之後，對碟盤以及連結器這兩個部份，利用減少材料的成本來降低生產成本，而結果顯示，碟形天線 AU2 在經由最佳化流程後，目標函數由 39.29 元降低至 32.37 元，節省了約 7 元。以 AU2 去年的出貨量達一百萬台以上來計算，估計若當初在天線量產前能配合此研究所提出的最佳化分析流程，可節省約新台幣 700 萬元。在本研究中雖然是以 W 公司的碟形天線 AU2 為分析目標，建立起流固耦合分析的流程，但此流程同時也適用於一般外型的碟形天線。未來在設計新型的碟形天線時，可以利用本研究中所提出的分析流程，來減少設計時間，在滿足天線強度要求 49.

(60) 的同時，並降低生產的成本。. 5.2 未來展望由於支撐桿的模擬分析與實驗結果誤差過大，因此在本論文中並無針對支撐桿的部份做最佳化的分析；若能夠修正產生誤差的因素，可以針對支撐桿的壁厚、管半徑等來做最佳化的分析。關於風洞實驗部份，由於初次實驗時只量測單一方向的應變值，有欠周全，因此建議下次實驗時可以利用菊花型應變計（strain rosette），來量測三個方向的應變值，並與模擬分析的結果做比較。另外對於碟盤盤面的方向，必需將其精度的要求要提高，以減少實驗時兩邊對稱位置上的應變規，所讀取數值的差異。在本研究中主要探討的部份皆為碟形天線強度分析的部份，僅針對天線不要產生永久變形的情況，不過剛性也是需要考慮的一部份。碟盤的變形和支撐桿的偏擺都會影響接收訊號的能力，若 W 公司能進一步提供碟盤變形量以及支撐桿偏擺對接收訊號影響的規範，在針對碟形天線分析時也必需將上述兩種情況加以考慮。本模擬分析時所討論的分析情況，是假設碟形天線受到正面方向均勻強度的風吹，但在真實使用的情況，自然界的風是有可能從各個方向吹來的，雖說直觀上，天線在承受正面風向會產生最大的應力，但側向風所造成碟盤以及接收器的位移的分析，需加入到剛性分析下。因此對碟形天線在不同方向風吹. 50.

(61) 下的分析，也是必需加以考慮的一個分析情況。另外還有模態分析的部份，對碟形天線來說，其中一項規範是天線的第一共振頻率必需在規範值以上，不過對此分析，本論文中所使用簡化之模型便不適用，因此必需建立較完整的模型，較複雜的邊界條件設定是其困難之處。本研究在附錄 A 中提出另一種流固耦合的方法－ALE 法，其優點為只需建立單一模型，使用單一求解核心，就能求得碟形天線在受風吹時的應力應變分析。但誤差過大且與原廠代理商溝通後，還是未能找出誤差的主因。因此未來若能找出分析產生誤差的原因，建立正確的邊界條件以及各項係數，利用此法就可以免除整合計算流體力學分析軟體與有限元素分析軟體的麻煩，同時可以節省租用兩套軟體的費用。. 51.

(62) 參考文獻 [1]. 彭玉賢，「直播衛星在台灣及亞太地區的市場潛力」，公共電視台研究發展部，http://www.pts.org.tw/~rnd/p4/991025.htm，民國 88 年。. [2]. J. M. T. Penrose, D. R. Hose, C. J. Staples, I. S. Hamill, I. P. Jones and D. Sweeney, “Fluid Structure Interactions: Coupling of CFD and FE,” CAD-FEM User’s Meeting International FEM Technologietage, 2000.. [3]. D. Filsinger, J. Szwedowicz and O. Schafer, “Approach to Unidirectional. Coupled. CFD–FEM. Analysis. of. Axial. Turbocharger Turbine Blades,” Journal of Turbomachinary, Vol. 124, pp. 125-131, January 2002. [4]. C. H. Tseng, “MOST 1.1 Manual,” Technical Report No. AODL-96-01, National Chiao Tung University, January 1996.. [5]. C. H. Tseng and J. S. Arora, “On Implementation of Compational Algorithms of Optimal Design 1: Preliminary Investigation,” International Journal for Numerical Methods in Engineering, Vol. 26, pp. 1365-1384, 1988.. [6]. J. S. Arora, Introduction to Optimum Design, Chapter 15 “Discrete Variable Optimum Design Concepts and Methods,” ELSEVIER, 2004.. [7]. 童寶鴻，「振動圓柱對渠道內熱流場影響之研究－第二章ALE 法之運動學理論與數學模式」，國立交通大學機械工程學系博士. [8]. 論文，民國92年六月。 J. Donea, S. Giuliani and J.P. Halleux, “An Arbitrary 52.

(63) Lagrangian-Eulerian Finite Element Method for Transient Dynamics Fluid-Structure Interactions,” Computer Methods in Applied Mechanics and Engineering, Vol. 33, pp. 689-723, 1982.. 53.

(64) 附錄A. ALE 法進行流固耦合分析. A.1 前言 Arbitrary Lagrangian-Eulerian （ ALE ）是結合拉格朗日（Lagrangian）和尤拉（Eulerian）座標系統所衍生出的座標系統，能有效處理移動邊界的問題，並且在計算的過程中控制網格的變形。此方法最大好處在於只需要建立單一模型、使用單一求解核心，就可以求得碟盤上的應力應變分析結果，並不需要進行兩種軟體的整合。. A.2 ALE 法之網格描述與座標系統[7] 現今描述流體運動的座標系統可分為兩種，一為拉格朗日座標系統，另一為尤拉座標系統。利用拉格朗日座標系統描述流場時，觀察流場的焦點集中在各個流體質點上；在數值計算的過程中，計算網格將隨流體質點的運動而移動。反之，尤拉座標系統觀測流場的焦點固定在空間中，計算網格也固定在計算區域中，不隨流體質點運動。 ALE 的概念是綜合拉格朗日與尤拉兩種座標系統，將觀測流場的焦點定在新定義的參考定義域上，此參考定義域以任意速度 w 移動。根據 Donea 等[8]的整理，ALE 座標描述系統與拉格朗日、尤拉座標系統間的關係，由 w 決定： (1). 當 w=0 時，參考定義域與計算網格皆靜止不動，此時 ALE 描述流場的觀點與尤拉座標系統相同。. (2). 當 w=v 時，v 為流體質點的運動速度。此時考定義域的移動速度與流體質點速度相同，ALE 描述流場的觀 54.

(65) 點則與拉格朗日座標系統相同。 (3). 當 w≠v≠0 時，此時參考定義域描述流場的觀點即為 ALE 法。. A.3 ALE 法分析情況簡介對於使用 ALE 法來分析流固耦合問題，可以將問題分成三種情況分析： (1). 用基於拉格朗日格式的有限元素方法模擬結構的變形和應力；. (2). 用基於純尤拉格式的有限體積方法描述材料（包括氣體和液體）流動；. (3). 對通過流體與固體介面傳遞相互作用的流體－結構耦合分析，採用基於混合的拉格朗日格式和純尤拉格式的有限元素與有限體積技術，完成全耦合的流體結構相互作用計算。. 在第二章中，處理流體-結構相互作用問題時，是將流體產生的力作為「預先確定」的載荷作用到結構上進行分析，而 ALE 法則不然，它是通過直接耦合結構網格（拉格朗日網格）和流體材料網格（尤拉網格）間的回應自動地、精確地算出每一時間步流－固介面處的物理性質。在這個過程中，一方面，尤拉材料流動引起的壓力載荷通過耦合演算法自動作用到結構的有限元素網格上，在這種壓力作用下，結構的有限元素網格將發生變形，結構的變形也反過來影響尤拉材料的流動和壓力值，這種結構變形和流體載荷間的. 55.

(66) 相互影響得到完全耦合的流體－結構回應。由於原先使用的 ABAQUS/STANDARD 並無提供解決流固耦合的功能，因此本研究選擇使用 LS-DYNA 有限元素軟體，其功能在下節中簡介。. A.4 軟體簡介 A.4.1. LS-DYNA. LS-DYNA 為泛用型的有限元素軟體，由美國 Livermore Software Technology Corp.（簡稱 LSTC）所開發。與其它軟體最大不同之處，強調「 One Model ， One Code ， Multi-Result」，只需建立一次有限元素模型，利用 LS-DYNA 核心求解程式，即可求解各式不同的物理現象及多階段分析，最重要的是其中包含了以 ALE 法處理流固耦合的問題；且內建專用前處理器 FEMB 與後處理器 LS-Pre/Post，相容性極高，使用方便。. A.4.2. 前處理器 FEMB. FEMB（Finite Element Model Builder）為美國知名軟體公司 ETA 所開發，專為 LS-DYNA 量身訂作的前處理程式，有以下特點：人性化的視窗介面、學習時間短。直覺性的操作模式。 CAD 功能。自動網格功能、網格修正調整。. 56.

(67) 豐富的 LS-DYNA Contact 以及材料設定、支援最新版本 LS-DYNA 核心。與 CAD 軟體整合度高。不需要額外指令輸入。直接輸出 LS-DYNA 專有 Input deck 檔案。由於 FEMB 為 LS-DYNA 內建專用的前處理器，親切的使用介面，可完全處理建構模型、網格切割、邊界條件設定等前置作業，不需再另外手動修改 input 檔等等優點，因此選擇其作為前處理器。. A.4.3. 後處理器 LS-Pre/Post. LS-Pre/Post 為 LSTC 原廠所設計的後處理程式。 LS-Pre/Post 具備絕佳的數值處理能力可直接讀取 LS-DYNA 的計算結果，進行數值資料的彙整及二次運算。可顯示板殼厚度、輸出各項力學資料：應力、應變、塑性應變、溫度、位移、速度、加速度、成型極限圖、能量分佈及動態視覺呈現。與另一個由 ETA 所開發的後處理器 PostGL 相比， LS-Pre/Post 在觀看動態結果的速度更快，並且提供比 PostGL 更多樣性的功能，故選擇其作為觀看 LS-DYNA 分析結果的後處理器。. A.5 模型建立與邊界條件設定在這裡所考慮的模型只有碟盤與支撐桿兩個，模擬出來. 57.

(68) 的結果直接與風洞實驗做比較。模擬所設定的條件，皆與第二章中相同。以下以條列式說明： (1). 以 FEMB 將匯入的碟盤與支撐桿的 IGES 檔建立網格和設定邊界條件，並將兩個物件之間以剛體元素連接起來，剛體元素視為連結器的簡化模型；由於模型為對稱，且網格節點數目龐大，因此將模型簡化為原來的 1/2。材料係數如表 2.1。. (2). 建立流體材料 Euler 網格，其材料性質為空氣，如表 2.2。. (3). 將流體的 Equation of State（EOS）給定為 Linear Polynomial Model 的形式，並將其參數設定為理想氣體；EOS 為決定流體施加到結構上壓力的方程式。. (4). 完成的模型如圖 A.1，方形網格即為流體材料 Euler 網格，其中圖 A.1 圈起來的節點為流體流入的邊界，視為風洞中風的流速，因此節點的速度設為 40m/s。. (5). 啟動 LS-DYNA 分析，分析結果以 LS-Pre/Post 觀看。. A.6 分析結果與討論由於 LS-DYNA 分析是以暫態方式進行分析，並不像 STAR-CD 是以穩態方式進行，取模擬結果中，出現數值比較穩定的一段時間的階段做為分析的結果。圖 A.2 為 LS-DYNA 碟盤 von Mises 等效應力分佈等高線圖，與圖 2.12 ABAQUS 的結果比較，應力集中的部份同樣都是位於四個固定點上，最大的應力產生在中間兩個固定點的外側。. 58.

(69) 但不同的是在圖 2.12 中，最大應力的節點，其值約為 103MPa，而圖 A.2 中最大的應力值約為 50MPa。圖 A.3 為 LS-DYNA 碟盤 x-y 方向的應變分佈等高線圖（座標方向如圖上所標示），圖 A.4 為 LS-DYNA 支撐桿的應變分佈等高線圖，其讀取數值的節點位置如圖 A.3、A.4 所示，數值結果如表 A.1 所列。經由式 2.1 轉換後的結果與風洞實驗結果比較，以及誤差比例列於表 A.2，在碟盤上模擬分析與實驗結果的誤差為 77%，在支撐桿上的誤差為 60.7%。由於模型只考慮了碟盤以及支撐桿這兩個部份，忽略了前方接收器，簡化的模型有可能是造成分析誤差的原因。另外再比較 STAR-CD 與 LS-DYNA 的流場分析結果，將兩個分析流程的模型皆簡化為只剩下碟盤以方便比較；圖 A.5 為 STAR-CD 的速度流場分佈圖，而圖 A.6 為 LS-DYNA 分析結果的速度分佈圖，觀察兩套軟體的流場速度分析結果並不盡相同，尤其在 STAR-CD 中，碟盤前方的有一個明顯的流體減速區，但在 LS-DYNA 中並不明顯，因此推測 LS-DYNA 分析誤差的可能原因之一，是產生在流場分析的部份。. 59.

(70) 表 A. 1. 500255 5001552. LS-DYNA 碟盤及支撐桿上量測點各方向的應變值 E11 -165 -287. E22 84 113. E33 31 380. E12 -157 10. E13 6 370. E23 14 26. 單位：micro-strain. 表 A. 2. LS-DYNA 碟盤及支撐桿模擬與實驗的比較表. 模擬實驗模擬與實驗的誤差(%). 碟盤. 支撐桿. 119. 432. 517.5. 1100.5. -77.0. -60.7 單位：micro-strain. 60.

(71) 40m/s. Ground 圖 A. 1. LS-DYNA 建立之模型與邊界條件設定. 61.

(72) 單位：GPa. 圖 A. 2. LS-DYNA 碟盤 von Mises 等效應力分佈等高線圖. Node 500255. 圖 A. 3. LS-DYNA 碟盤 x-y 方向的應變分佈等高線圖 62.

(73) Node 5001552. 圖 A. 4. LS-DYNA 支撐桿的應變分佈等高線圖. 63.

(74) 單位：m/s. 圖 A. 5. 圖 A. 6. STAR-CD 速度流場分佈圖. LS-DYNA 分析結果速度分佈圖. 64.

(75) 附錄B. 壓力轉換 C++程式語言原始碼. //2004.10.07 #include <iostream.h> #include <iomanip.h> #include <fstream.h> #include <math.h> double node[1297][4];//節點資料 int i,j,p,q; int element[1264][5];//該元素由哪幾個 node 組成 double center[1264][5];//該元素的中心點位置,及前後壓力 double front[1170][4];//前壓力 double back[1166][4];//後壓力 double x,y,z; int count=0; double R=15.; void main() { ifstream input("node.txt"); while(!input.eof()) { for(i=0;i<1297;i++) for(j=0; j<4; j++) { input>>node[i][j]; } } input.close();//讀取 node 資料 ifstream input1("element.txt"); while(!input1.eof()) { for(i=0;i<1264;i++) for(j=0; j<5; j++) { input1>>element[i][j]; } } input1.close();//讀取 element 資料 ifstream input2("front.txt"); while(!input2.eof()) { for(i=0;i<1170;i++) { for(j=0; j<4; j++) { input2>>front[i][j]; cout<<front[i][j]; } 65.

(76) cout<<endl; } } input2.close();//讀取前方壓力值和網格中心點位置 ifstream input3("back.txt"); while(!input3.eof()) { for(i=0;i<1166;i++) for(j=0; j<4; j++) { input3>>back[i][j]; } } input3.close();//讀取後方壓力值和網格中心點位置 //calculate center 計算出網格中心點位置 for(i=0; i<1264; i++) { for(p=0; p<1297; p++) { for(q=1; q<5; q++) { if(node[p][0]==element[i][q]) { center[i][0]+=node[p][1]; center[i][1]+=node[p][2]; center[i][2]+=node[p][3]; count++; } } } center[i][0]=center[i][0]/4; center[i][1]=center[i][1]/4; center[i][2]=center[i][2]/4; cout<<count; count=0; } //print the center[][] information ofstream output3; output3.open("center.txt"); for ( i=0; i<1264; i++) { output3<<center[i][0]<<setw(15)<<center[i][1]<<setw(15)<<center[i][2]<<endl; } output3.close(); //進行搜尋並運算壓力值 count=0; for(i=0; i<1264; i++) { for(j=0; j<1170; j++) { 66.

(77) x=center[i][0]-front[j][0]; y=center[i][1]-front[j][1]; z=center[i][2]-front[j][2]; if(x*x+y*y+z*z<R*R) { center[i][3]+=front[j][3]; count++; } } center[i][3]=center[i][3]/count; cout<<endl<<count; count=0; } count=0; for(i=0; i<1264; i++) { for(j=0; j<1166; j++) { x=center[i][0]-back[j][0]; y=center[i][1]-back[j][1]; z=center[i][2]-back[j][2]; if(x*x+y*y+z*z<R*R) { center[i][4]+=back[j][3]; count++; } } center[i][4]=center[i][4]/count; count=0; } //輸出壓力值 ofstream output; output.open("Pfront.txt"); for ( i=0; i<1264; i++) { output<<center[i][3]<<endl; } output.close(); ofstream output1; output1.open("Pback.txt"); for ( i=0; i<1264; i++) { output1<<center[i][4]<<endl; } output1.close(); }//void main(). 67.

(78) 附錄C. MOST Input C 語言程式原始碼. #include <math.h> #include <stdio.h> #include <stdlib.h> #include <iostream.h> cusermf (x,f,nv,nobj,kk) int nv,nobj,kk; double x[],*f; { int ierr =0; double w1=28;/*the price of steel plate*/ double w2=30;/*the price of stell pipe*/ double pi=3.1415926; double Aconnect=12000;/*connector's area*/ double Aradius=pi*(46*46+82.5*82.5);/*原來半徑所圍成的面積*/ double Anew=pi*(pow(x[0],2)+pow(fabs(x[1]-x[3])/2,2));/*新的半徑所圍成的面積*/ double Adish=214366;/*Area of dish*/ double d=7850.0/1000000000.0;/* kg/m^3 */ if(x[4]<0.6) w2=32; *f = Aconnect*(Anew/Aradius)*d*w1+Adish*x[4]*d*w2; nobj = nobj; /* never used */ nv = nv; /* never used */ kk = kk ; /* never used */ return ierr; } cusercf (x,g,nv,neql,nc,kk) int nv,neql,nc,kk; double x[],g[]; { int ierr=0; double mises=0;/*frequency*/ FILE *openfile; int i,j,k,m; double center[1264][3]; int element[1264][5]; double node[1297][4]; char s[99]; double check;/*use for check*/ double xx,yy,y;/*distance between two nodes*/ int number;/*save node's number*/ int fix[4][4]; x[2]=x[0];/*let the x[2]=x[0]*/ /*renew the dish_for_abaqus.inp*/ system("del dish_for_abaqus.inp"); system("copy backup.inp dish_for_abaqus.inp"); /*read the center element data*/ 68.

(79) openfile=fopen("center.txt","r"); for(i=0; i<1264; i++) { for(j=0; j<3; j++) { fscanf(openfile,"%s",s); center[i][j]=atof(s); } } fclose(openfile); /*read the element data*/ openfile=fopen("element.txt","r"); for(i=0; i<1264; i++) { for(j=0; j<5; j++) { fscanf(openfile,"%s",s); element[i][j]=atof(s); } } fclose(openfile); /*read the node data*/ openfile=fopen("node.txt","r"); for(i=0; i<1297; i++) { for(j=0; j<4; j++) { fscanf(openfile,"%s",s); node[i][j]=atof(s); } } fclose(openfile); /*search for the nearest element with {x[0] and x[1]} and {-1*x[0] and x[1]}*/ /*with {x[2] and x[3]} and {-1*x[2] and x[3]}*/ for(k=0; k<2; k++) { for(j=0; j<2; j++ ) { check=1000000;/*use for check*/ for(i=0; i<1264; i++) { if(j==0) { xx=center[i][0]-x[0+2*k]; y=center[i][1]-x[1+2*k]; } if(j==1) { xx=center[i][0]-(-1)*x[0+2*k]; y=center[i][1]-x[1+2*k]; 69.

(80) } if(xx*xx+y*y<check) { number=i; check=xx*xx+y*y; } } /*re write the new position*/ for(i=0; i<4; i++) fix[j+2*k][i]=element[number][i+1]; /*find the node in this element near the (1000,1000) or (-1000,-1000) */ if(j==0) xx=1000; if(j==1) xx=-1000; yy=1000; for(i=0; i<1297; i++) if(node[i][0]==fix[j+2*k][0]) check=(xx-node[i][1])*(xx-node[i][1])+(yy-node[i][2])*(yy-node[i][2]); number=0; for(i=1; i<4; i++) { for(m=0; m<1297; m++) if(node[m][0]==fix[j+2*k][i]) y=(xx-node[m][1])*(xx-node[m][1])+(yy-node[m][2])*(yy-node[m][2]); if(y<check) { check=y; number=i; } } /*nodes are arrayed by clockwise*/ for(i=0; i<4; i++) { openfile=fopen("dish_for_abaqus.inp","r+"); while(check!=fix[j+2*k][number]) { fscanf(openfile,"%s",s); check=atof(s); } fseek(openfile,0,1); for(m=0; m<1297; m++) if(node[m][0]==fix[j+2*k][number]) y=node[m][3]; if(i==0) { if(j==0) check=1; 70.

(81) if(j==1) check=-1; fprintf(openfile," fprintf(openfile," fprintf(openfile," } if(i==1) { if(j==0) check=1; if(j==1) check=-1; fprintf(openfile," fprintf(openfile," fprintf(openfile," } if(i==2) { if(j==0) check=1; if(j==1) check=-1; fprintf(openfile," fprintf(openfile," fprintf(openfile," } if(i==3) { if(j==0) check=1; if(j==1) check=-1; fprintf(openfile," fprintf(openfile," fprintf(openfile," } fclose(openfile); number++; if(number==4) number=number-4;. ,%f",check*x[0+2*k]); ,%f",x[1+2*k]); ,%f",y);. ,%f",check*x[0+2*k]+7*j); ,%f",x[1+2*k]+7*(j-1) ); ,%f",y);. ,%f",check*x[0+2*k]-7*check); ,%f",x[1+2*k]-7); ,%f",y);. ,%f",check*(x[0+2*k]+7*(j-1))); ,%f",x[1+2*k]-7*j); ,%f",y);. } } } /*renew the boundary condition */ openfile=fopen("dish_for_abaqus.inp","r+"); while(s[0]!='*'||s[1]!='N'||s[2]!='S'||s[13]!='0') fscanf(openfile,"%s",s); fseek(openfile,0,1); fprintf(openfile,"\n"); for(i=0; i<4; i++) 71.

(82) for(j=0; j<4; j++) fprintf(openfile,"%d,",fix[i][j]); fclose(openfile); /*key in the thickness */ openfile=fopen("dish_for_abaqus.inp","r+"); while(s[0]!='M'||s[1]!='A'||s[2]!='T'||s[13]!='R') fscanf(openfile,"%s",s); fseek(openfile,0,1); fprintf(openfile,"\n"); fprintf(openfile,"%f,",x[4]); fclose(openfile); /*run ABAQUS*/ system("abaqus job=dish_for_abaqus int"); system("del dish_backup.dat"); system("rename dish_for_abaqus.dat dish_backup.dat"); /*del the file we never used*/ system("del dish_for_abaqus.com"); system("del dish_for_abaqus.fil"); system("del dish_for_abaqus.msg"); system("del dish_for_abaqus.mdl"); system("del dish_for_abaqus.odb"); system("del dish_for_abaqus.par"); system("del dish_for_abaqus.pes"); system("del dish_for_abaqus.pmg"); system("del dish_for_abaqus.prt"); system("del dish_for_abaqus.sta"); system("del dish_for_abaqus.stt"); /*read the mises*/ openfile=fopen("dish_backup.dat","r"); for(i=0; i<4; i++) { while(s[0]!='M'||s[1]!='A'||s[2]!='X'||s[6]!='M') fscanf(openfile,"%s",s); fscanf(openfile,"%s",s); } /*get the Mises stress from dish_backup.dat*/ mises=atof(s); fclose(openfile); g[0] =mises-160/1.5; g[1] =625-((x[0]-x[2])*(x[0]-x[2])+(x[1]-x[3])*(x[1]-x[3])); nc = nc; /* never used */ neql = neql; /* never used */ kk = kk; /* never used */ nv = nv; /* never used */ return ierr; } 以下略. 72.

(83) 附錄D MOST 輸入資料檔內容 tit = Optimization of the Fix position and thickness of Dish nv =5 nobj = 1 neql = 0 niql = 2 ntrs = 1000 ipr =1 iact = 5 igrad = 1 del = 1.0000e-02 acs = 1.0000e-02 acv = 1.0000e-01 act = 1.0000e-10 x[1] = 46 25 100 x[2] = -30 -185 100 x[3] = 1 1 1 x[4] = -185 -185 -185 x[5] = 0.7 0.4 0.8 tdv[5] = discrete 6 set[5] = 0.4 0.5 0.6 0.65 0.7 0.8. 73.

(84)