華中科技大學 碩士學位論文 散熱結構拓撲優化 姓名：李朕 申請學位級別：碩士 專業：機械電子工程 指導教師：熊蔡華 20080501 華 中 科 技 大 學 碩 士 學 位 論 文 I 摘 要 摘 要 拓撲優化方法在連續體結構、MEMS 機構和散熱結構等結構優化設計方面得到 了廣泛應用。拓撲優化至今仍是結構優化中最具挑戰性的研究領域。本文針對散熱 結構拓撲優化問題，提出了改進的密度懲罰法，并進行了理論和數值實驗研究。 文中首先介紹了原始密度懲罰法基本理論，研究了密度懲罰法中的材料模型和 有限元熱傳導模型。接著討論了拓撲優化中溫度場的有限元求解格式及二維和三維 單元導熱矩陣的構造方法，對如何選擇插值函數、確定形函數待定系數以及構建物 理矩陣進行了深入的研究。運用密度懲罰法和熱傳導分析模型研究了散熱結構拓撲 優化模型。隨后針對散熱結構拓撲優化中普遍存在的數值不穩定性問題進行了討論， 并研究了消除數值不穩定性的方法。 接下來提出了改進的密度懲罰法模型，研究了拓撲優化中各向異性導熱材料拓 撲優化與方向角優化問題，提出了基于溫度梯度的方向角優化準則算法。研究了含 各向異性導熱材料區域溫度場有限元求解模型和數值算法。通過引入旋轉矩陣研究 了材料角度對拓撲優化目標函數的影響。 最后研究了散熱結構拓撲優化設計的數值算法和處理技巧，并開發了三維散熱 結構拓撲優化和各向異性導熱材料拓撲-角度混合優化的計算程序，針對一系列典型 的散熱結構拓撲優化算例進行了數值實驗。 關鍵詞：關鍵詞： 拓撲優化 散熱結構 SIMP 各向異性材料 華 中 科 技 大 學 碩 士 學 位 論 文 II Abstract Topology optimization has been widely used in continuum structures, MEMS and heat radiation structures design. Topology optimization is one of the most challenging tasks in structural design. In this thesis, a modified density penalization method has been proposed. Numerical study has been carried out to investigate the efficiency of this method. First, a fundamental theory of SIMP method is introduced in this thesis. Moreover, the material model and finite element method (FEM) in density penalization method are discussed. Key steps of 2-dimenation and 3-dimenation element matrix formulation, including interpolating function selection and determination of the coefficients of shape function are implemented. The SIMP method and the heat conduction analyze model are integrated to formulate topology optimization model of heat radiation structures. Then numerical instability in topology optimization methods such as checkerboards, mesh-dependence and local minima are discussed and corresponding methods to avoid those problems are listed in the study. A modified density penalization model is also proposed and the issue of anisotropic material in heat conduction structure topology optimization is studied. Besides, a criteria method based on temperature gradient for the optimization of material orientation is proposed. Finally, a computer program for topological optimization of heat conduction structures is developed using SIMP. Furthermore, numerical experiments for a series of representative examples are implemented, which shows that this method is feasible and effective. Keywords: Topological Optimization, Heat Radiation Structures, SIMP, anisotropic material 獨創性聲明 本人聲明所呈交的學位論文是我個人在導師指導下進行的研究工作及取得的 研究成果。盡我所知，除了文中已經標明引用的內容外，本論文不包含任何其他人 或集體已經發表或撰寫過的研究成果。對本文的研究做出貢獻的個人和集體，均已 在文中以明確方式標明。本人完全意識到本聲明的法律結果由本人承擔。 學位論文作者簽名： 李朕 日期： 2008 年 5 月 24 日 學位論文版權使用授權書 本學位論文作者完全了解學校有關保留、使用學位論文的規定，即：學校有權 保留并向國家有關部門或機構送交論文的復印件和電子版，允許論文被查閱和借 閱。本人授權華中科技大學可以將本學位論文的全部或部分內容編入有關數據庫進 行檢索，可以采用影印、縮印或掃描等復制手段保存和匯編本學位論文。 保密，在_年解密后適用本授權書。 本論文屬于 不保密?。 （請在以上方框內打“” ） 學位論文作者簽名：李朕 指導教師簽名：熊蔡華 日期： 2008 年 5 月 24 日 日期：2008 年 5 月 24 日 華 中 科 技 大 學 碩 士 學 位 論 文 1 1 緒論緒論 1.1 課題的來源課題的來源 本學位論文得到以下項目的聯合資助： 國家自然科學基金重點項目“機電表面功能結構及相關熱物理問題的基礎研究” （批準號： 50436010） 國家重點基礎研究發展計劃（973）項目“數字化制造基礎研究” 。 （批準號： 2005CB724100） 1.2 課題的目的和意義課題的目的和意義 隨著微電子領域對產品性能的無限追求，芯片集成度不斷提高，帶來致命的高 熱流密度問題，電子器件的冷卻問題越來越突出。英特爾公司負責芯片設計的首席 執行官帕特-蓋爾欣格曾指出， 如果芯片耗能和散熱的問題得不到解決， 到 2005 年芯 片上集成了 2 億個晶體管時，就會熱得像“核反應堆”，2010 年時會達到火箭發射時 高溫氣體噴射的水平，而到 2015 年就會與太陽的表面一樣熱。目前面積 2 2cm芯片 上的功耗可達125W。為了解決散熱問題，芯片制造商在發布新一代產品的時候不約 而同的采用巨型的風冷散熱器，甚至用水冷系統進行散熱1。但是在量產型號中，由 于風冷系統可靠性高，成本低廉等優點，仍然是處理器散熱的首選。除了微電子領 域，在化工等領域同樣存在高熱流密度問題。隨著現代制造技術的不斷進步，有效 散熱空間日趨減小，許多場合散熱空間是封閉或半封閉的，要求散熱系統能在短時 間內通過極高的熱流密度。因此，具有高熱流密度機電器件功能結構表面散熱問題 已成為工程熱物理、機械等學科急待解決的關鍵科學和技術問題。 本論文旨在根據復雜表面結構的熱功能需求，優化熱功能結構的拓撲和形狀， 為微電子及化工等領域強化傳熱結構的優化設計提供理論依據。 本論文提出了改進的密度懲罰法。提出了基于溫度梯度的材料方向角優化準則 算法，解決了各向異性導熱材料拓撲與角度混合優化問題。運用該方法開展了大量 的數值實驗，取得了較好的效果，對散熱結構的優化設計具有較好的指導意義。 華 中 科 技 大 學 碩 士 學 位 論 文 2 1.3 拓撲優化基本描述拓撲優化基本描述 由于人類掌握的資源有限，最大化地利用現有原材料一直是人們追求的目標， 也是許多研究領域不變的話題。早期的產品設計并沒有涉及優化的概念，而是依賴 設計人員的經驗進行設計。在設計過程中，通過簡單地試算和調整有時可以得到局 部最優的結果。但是設計余量往往留得過大，造成了浪費。在六十年代末隨著計算 機技術的發展， 尺度優化已經應用于工程結構的優化設計中， 該方法得到了A. Miele、 E.J.Haug和J.Taylor等人的發展2。隨著有限元技術逐漸成熟，七十年代形狀優化蓬 勃發展。進入八十年代，伴隨著小型化高性能計算機的出現和數值計算方法走向成 熟，拓撲優化進入了活躍發展時期。 從工程設計的角度看，結構優化設計大致有三個層次：拓撲優化、形狀優化、 尺寸優化，其中拓撲優化是結構設計的核心。 結構拓撲優化的研究始于桁架結構優化。Maxwell早在1854年就進行了應力約 束下最小桁架的基本拓撲分析3。 在隨后的一百年間受理論和工具的限制，拓撲優化 發展緩慢。直到20世紀60年代隨著計算技術的發展，Dom等人把數值方法引入到 該領域，拓撲優化的研究才重新活躍起來4。 由于受到理論限制，因此早期的結構優化研究主要是針對設計域形狀和拓撲固 定的尺寸優化問題展開的。后來隨著結構邊界優化問題的提出，形狀優化方法應運 而生3。 當今在航空和汽車制造業，尺寸和形狀優化技術已經廣泛用于結構和零部件 設計。形狀優化方法也常用于電磁、電化學和聲學零件的設計。目前已有許多成功 的算法可以處理形狀優化問題5。 然而， 采用邊界優化等形狀優化方法進行產品設計， 設計質量依賴于人為給定的初始設計。如果初始拓撲和邊界條件定義不當，得到的 設計結果往往不能令人滿意。采用拓撲優化可以避免這個問題。Bendse3 對拓撲優 化理論進行了開創性研究，實現了拓撲優化與形狀優化的結合，為高效高質量的快 速原型設計提供了一個非常有效的理論工具。 結構拓撲優化包括離散結構拓撲優化和連續體結構拓撲優化兩大類。目前國內 外研究主要針對連續體結構。因為離散結構優化問題的目標函數和約束條件不連續、 不可微，可行域退化為不連通的集合。普通的求解方法難以處理這類優化問題。經 過近30年的發展，連續體結構優化基本理論已經比較成熟，目前正向工程熱物理、 光學和MEMS等多學科交叉領域發展。 華 中 科 技 大 學 碩 士 學 位 論 文 3 拓撲優化按模型不同可以分為兩大類：微結構法和宏觀法。 微結構法又稱為材料法，把材料看作有限個微小單元的集合，以不同的材料分 布代表拓撲，主要代表是均勻化方法和密度懲罰法模型。 宏觀法又叫幾何法，主要包括進化法、泡泡法和近年來發展迅速的水平集拓撲 優化方法。幾何法用設計域中的圓和孔等基本幾何元素代表拓撲，比如進化法中采 用移動的邊界來表示結構的變化，水平集方法采用移動的零水平集界面來描述拓撲 結構。 1.4 國內外發展現狀國內外發展現狀 盡管結構拓撲優化設計難度大，但它所帶來的潛在經濟效益卻是尺寸優化和形 狀優化難以達到的。因此，許多學者致力于拓撲優化的研究。目前研究較多的結構 拓撲優化方法主要有均勻化方法、變密度法、進化法和水平集方法。 均勻化方法均勻化方法 復合材料的大量應用降低了產品制造成本，但是也增加了產品設計的復雜性。 均勻化方法（homogenization）通過構造等價的微結構材料模型來描述復合材料，以 微結構單胞的幾何尺寸及方位角為設計變量，把拓撲優化問題轉化為數學規劃問題， 再通過某種優化算法得到使結構的某種性能指標達到最優時設計區域各處的材料分 布，從而得到最優的拓撲結構。該方法為復合材料產品設計提供了一個很好的工具。 均勻化方法的理論基礎產生于上世紀70年代末。 利用Bensoussan等人發展起來的 基于攝動理論的周期性結構分析方法，Bendse3建立了材料微結構尺寸與材料宏觀 彈性特性之間的關系，提出了均勻化拓撲優化方法。由于均勻化方法具有較嚴謹的 數學基礎，因而已成為拓撲優化研究中的重要方法。國內大連理工大學劉書田和程 耿東等 6在復合材料彈性剛度張量預測、熱膨脹系數張量預測、導熱系數張量預測、 局部應力和位移分布分析等方面開展了均勻化方法的應用研究。 目前，均勻化方法已能成功處理多工況的二維、三維連續體結構拓撲優化問題。 變密度法變密度法 另一種基于微結構法的拓撲優化方法是以固體各向同性材料懲罰模型(SIMP: Solid Isotropic Material with Penalization)為代表的變密度法。該模型假設材料為各向 華 中 科 技 大 學 碩 士 學 位 論 文 4 同性材料，僅以密度為設計變量，不需要引入細微結構和附加的均勻化過程。程序 實現簡單，計算效率高7。Sigmund7發表了一個簡練的拓撲優化程序，引起了拓撲 優化研究者的廣泛興趣。從中派生出了許多優化程序用于解決不同領域中的拓撲優 化問題，例如柔性體結構優化、光學通路優化等。 進化法進化法 Xie和Steven8, 9提出了漸進結構優化方法（ESO：Evolutionary Structural Optimizatio） 。ESO的基本思想很簡單，即通過逐步去除設計域中無效或低效材料獲 得最優結構的一種結構優化方法。該方法在最小柔順性設計、特定自振頻率等拓撲 優化設計中獲得了極大成功。該算法通用性好，不僅能解決尺寸優化，還可同時實 現形狀與拓撲優化。ESO法的局限在于刪除的材料不能被恢復，而雙向漸進結構法 (BESO：BiDirectional Evolutionary Structural Optimization)彌補了ESO法的不足。 BESO法在優化過程中既可以刪除材料單元，又可以添加材料單元 10。Young11成 功地把雙向漸進結構法應用到三維拓撲結構優化中。 水平集方法水平集方法 Osher 和Sethian提出了水平集方法(Level Set Method)12，該方法可以有效描述 界面運動問題。Wiegmann13首先把水平集方法應用到彈性體結構的邊界設計中，采 用移動的邊界來表示設計空間中材料的增加或刪除。隨后Allair等14把經典的形狀 靈敏度算法結合到水平集方法中，實現了水平集方法在彈性體結構拓撲與形狀構優 化的整合。Wang 等15提出了結構的向量水平集表示方法。Seung-Hyun Ha16, 17等將 水平集方法應用于熱傳導問題的形狀和拓撲優化，取得了較好效果。國內也在水平 集拓撲優化領域展開了廣泛的研究。大連理工大學梅玉林18提出了改進的水平集拓 撲優化算法，該算法適用于一般目標函數、多材料、多約束和多載荷工況的結構拓 撲優。賈海朋19基于進化算法的思想提出了ESO 插孔的水平集拓撲優化算法，解決 了水平集方法對任意結構、任意初始邊界的拓撲優化問題。莊春剛20研究了基于水 平集方法的多載荷散熱結構拓撲優化問題。華中科技大學蔣良杰21在SIMP法和水 平集方法散熱結構拓撲優化設計方面展開了研究，開發了MATLAB和FEMLAB下 的拓撲優化計算程序。 結構優化中的優化算法結構優化中的優化算法 華 中 科 技 大 學 碩 士 學 位 論 文 5 結構優化問題的解決最終歸結為求解數學優化問題。為了提高優化的經濟性， 需要選擇時間復雜度與空間復雜度都相對較低的優化算法。對于所有優化問題都可 歸結為在給定條件(例如約束條件)下求目標函數的最優值問題。實際工程優化問題 中，優化問題的約束條件和目標函數不僅是非線性的，而且通常是隱式函數，直接 求解非常困難，所以優化算法的選用至關重要。對于不同層次的優化問題需要選用 不同的優化算法。目前在拓撲優化中運用較多的優化算法有OC法(Optimality Criteria)、SLP法(Sequential Linear Programming)和MMA法(Method of Moving Asymptotes)等。 OC法是基于直覺的準則法，把數學中最優解應滿足的Kuhn-Tucker條件作為最 優結構應滿足的準則。用優化準則來更新設計變量和拉格朗日乘子3。采用準則法不 需要計算目標函數的梯度等額外信息，所以計算量不大。在一般情況下，所需迭代 次數少，目標函數收斂較快。但其缺點也明顯，首先它只能保證特定情況和參數下 解的收斂性；其次由于采用準則法獲得的結果僅僅是人們所“期望”看到的結果， 所以不能保證獲得最優解；最后對不同類型的約束、變量、目標函數需導出不同的 優化準則，通用性差。 MMA法即移動漸進線法，最早由Svanberg22 提出。MMA法用一顯式的線性 凸函數來近似代替隱式的目標和約束函數。由事先確定的左、右漸進點和原函數在 各點的導數符號來確定每一步迭代的近似函數。當左、右漸進點分別趨近負無窮大 和正無窮大時，MMA法就等同于用SLP近似。該方法的優點在于它是全局收斂的， 所以MMA法能從理論上保證解的存在性，另外MMA法對初值不敏感，迭代過程 穩定。缺點是計算效率較低。 SLP法即序列線性規劃法。該法通用性好，但收斂速度慢，對初值敏感。數學 規劃法是以規劃論為基礎，理論嚴謹，適用面廣，且收斂性有保證。但由于計算量 大，對拓撲優化這一類含多個變量的優化問題不適合。20世紀70年代中期以后，結 構優化設計中的規劃法吸收了準則法的優點，根據力學特性進行了某些改進，如顯 式逼近、變量連接、選擇有效約束、引入倒數變量、采用對偶求解技術等，使計算 效率得到了顯著提高。 其它一些新興的優化算法也逐漸引起了拓撲優化研究者的興趣。比如遺傳算法 （GA：Genetic Algorithms） 、粒子群算法（PSO：Particle Swarm Optimization）和模 擬退火算法等。遺傳算法屬于半隨機算法，它具有自適應性，不需要描述問題的全 華 中 科 技 大 學 碩 士 學 位 論 文 6 部特點，搜索過程不受優化函數的連續性約束，也不需要導數等額外信息。但這種 方法存在結構分析重復次數過多、收斂速度慢等缺點，不適合求解大規模的工程優 化問題23。粒子群算法在圖像識別中體現出極高的效率，但在拓撲優化中的應用還 有待研究。 1.5 結構優化進展實例結構優化進展實例 下面簡要介紹國內外結構優化設計方法應用的一些實例，他們所開展的研究工 作對本論文有很好的借鑒意義： 由丹麥技術大學（Technical University of Denmark）機械工程學院 Ole Sigmund，Martin Bendse 等人領導的研究隊伍很早就在拓撲優化方面展開了研究。 早在上世紀 90 年代初期，該學院就在與丹麥政府合作的一個小型衛星支架優化項目 中實現了多約束下的三維結構拓撲優化。 圖 1.1 小型衛星支撐結構優化 圖1.1 a）為支架的設計要求，需要在60x60x80cm的空間中安置4臺攝像機、 望遠鏡頭、推進器、電池和電子通信設備等，設備總重量小于80kg，支撐結構的總 重量限制為12kg并需在衛星發射時能承受數G的加速度。整個設計域為288000個 立方體有限單元。最終優化結果如圖1.1 b）所示。 另一個典型的優化算例是三維橋梁的優化算例，如圖 1.2。 華 中 科 技 大 學 碩 士 學 位 論 文 7 圖 1.2 三維橋梁 從文獻和例子來看，目前結構拓撲優化設計已經取得了廣泛的理論研究成果， 比如上述衛星支架優化項目和三維橋梁優化項目。在應用方面也取得了一定成果， 目前已經出現了一些成熟的拓撲優化商業軟件。但在散熱結構優化方面的研究還比 較滯后，遠遠不能滿足復雜的工程需求。大部分研究還停留在數值實驗階段。需要 跳出這一局限，為散熱結構的工程設計提供高效的設計方法。本論文的散熱結構拓 撲優化正是在這種背景下提出的。主要解決三維散熱結構拓撲優化和面向各向異性 材料的拓撲-方向角混合優化問題。在本文第二章和第三章給出了詳細的介紹。 1.6 本文研究內容與章節安排本文研究內容與章節安排 本文針對散熱結構的拓撲優化問題，利用密度懲罰法進行了理論和數值實驗研 究。提出了改進的密度懲罰法，適用于各向異性材料拓撲-角度混合優化。利用 MATLAB開發了二維和三維散熱結構拓撲優化的數值算法和計算程序。全文內容具 體安排如下： 第一章為緒論，闡述了課題來源和項目研究的背景和意義，并詳細介紹了拓撲 優化問題的描述和各種拓撲優化方法的研究現狀。 第二章首先介紹了拓撲優化的數學模型和基于SIMP模型的密度懲罰法理論。 在 此基礎上分析了密度懲罰法散熱結構拓撲優化模型，為后續的數值實驗建立了理論 基礎。 第三章研究了拓撲優化中各向異性導熱材料方向角優化問題，提出了改進的密 華 中 科 技 大 學 碩 士 學 位 論 文 8 度懲罰法模型，建立了拓撲和方向角混合優化模型。提出了基于溫度梯度的方向角 優化準則算法。 第四章研究了散熱結構拓撲優化設計的數值算法和處理技巧， 并在MATLAB中 進行了三維散熱結構拓撲優化和各向異性導熱材料拓撲優化的幾種典型算例的數值 實驗。 第五章總結全文，并對下一步研究工作重點進行了分析和展望。 華 中 科 技 大 學 碩 士 學 位 論 文 9 2 密度懲罰法拓撲優化理論密度懲罰法拓撲優化理論 拓撲一詞最早源于希臘詞“topos” ，意思是地點、空間或區域。拓撲學研究區域 中物體的變形，比如四邊形、八面體在空間中的形變。拓撲區域的定義為： 3 所有 子集 （包括線和點的集合） 。 拓撲優化研究的問題是尋找 3 中最優拓撲使目標函數 （剛 度、自振頻率等）最小。從數學的觀點來講，所有的拓撲變形都可以視為變換或是 映射。規定拓撲的變換和逆變換都是具有唯一性的連續變換。這就為拓撲優化解的 連續性提供了保證。 研究拓撲優化問題離不開形狀優化的討論。在拓撲學中，區域的變形和大小并 不影響拓撲，而區域的增減對拓撲產生影響。聯系到工程設計中，設計的總體布置 往往比其設計尺寸更能決定設計的優劣。所以，單純的形狀優化很難使設計達到最 優。在產品的早期設計階段引入拓撲優化是提高設計質量非常有效的手段。從另一 個方面來講，單純的拓撲優化僅僅在設計域上生成一些孔洞，這些含有孔洞的拓撲 對設計者而言不具有實用意義，因此拓撲優化需要結合形狀優化生成有意義的產品 形狀。Qing24 結合形狀優化和拓撲優化，研究了熱傳導結構優化問題。Zhou25使用 商業拓撲優化代碼研究了形狀、拓撲和尺寸三者優化的整合問題。Zhou提出的方法 在優化問題中可以組合不同設計變量，為設計者提供了極大的自由度。Rafael26等學 者開發了ANSYS下的拓撲優化算法，并利用該代碼成功進行了馬達轉子的優化。近 二十年來，隨著計算機運算速度和存儲能力的提高以及有限元等分析工具的成熟， 采用計算機求解拓撲優化問題逐漸成為一個熱門領域。 2.1 優化模型優化模型 由于拓撲優化的復雜性，只有少量的問題可以用解析的方法人工進行求解。為 了滿足工程中的復雜要求，需要借助計算機用數值方法求解拓撲優化問題。因此解 決拓撲優化的挑戰首先在于如何在計算機中用離散模型描述拓撲；其次是如何建立 一個可供計算機求解的優化模型。正如緒論中所述，材料法假設結構的宏觀特性由 微尺度的單胞結構（類似自然界中組成材料的分子和原子）特性來決定。這種方法 雖然可以提供自然界中所有材料的精確數學描述，然而它表示一個普通零件所需要 華 中 科 技 大 學 碩 士 學 位 論 文 10 的存儲空間和計算量已經大大超出了現有計算機的能力，因而必須用簡化的材料數 學模型及有限個變量來描述拓撲結構。微結構法的思想就是把拓撲優化轉化為有限 個參數的尺度優化（sizing problem）問題。 2.1.1 材料模型材料模型 基于材料微結構的拓撲優化起源于均勻化方法。均勻化方法中材料模型分為兩 大類：一種是多層材料模型，另一種是含有孔洞的微結構。前者的均勻化方程可以 通過解析的方法求解，而后者的單胞模型只能用數值方法求解。為了方便計算機求 解拓撲優化問題，一般選擇第二種模型。 先考慮二維均勻化模型。設計域屬于歐拉空間（ 2 ） ，其中的材料單元微結 構特性可以由對稱矩陣( , , ) ij Ea b來描述。 圖 2.1 均勻化方法材料模型 a和b分別表示第e個單元材料在x和y方向的尺度。以上是最早的微結構材料 描述方法。 在 SIMP 法中簡化了材料模型， 僅用密度x一個變量來表示單元材料性質， 把優化問題歸結為固定域上的 0-1 整數規劃問題。 對單元e， 密度 e x 取 0 代表沒有材 料，取 1 代表有材料。有材料的區域對結構的剛度有貢獻，而無材料的區域為空白 區域（即不存在結構，不傳遞力） 。采用這樣的處理可以在把原問題轉化為一個固定 設計域上的優化問題。從而避免了由于結構變化導致的網格重新劃分和設計變量改 變的繁瑣過程。SIMP 法中的材料模型可以用密度函數來表示 3： 0 ( )( ) p E xxE= 1p (2.1) e a b 1 1 x y o 華 中 科 技 大 學 碩 士 學 位 論 文 11 其中( ),0( )1x xx是密度函數， 0 E是給定的固體材料彈性模量。密度 函數上的指數p是懲罰因子。因為拓撲優化關心的是材料的有無，所以密度函數選 擇離散函數(2.2)式： ( ) 1 /0 s s xif x xif = (2.2) Bendse在研究中發現，如果直接采用(2.2)式的離散密度函數，不引入含中間密 度的復合材料，無法保證拓撲優化中解的存在性。為了得到可行解，需要對原問題 進行一定的調整。通過引入密度處于0和1之間的復合材料把密度函數( )x松弛化 為一個連續函數，再通過加上懲罰因子使占有相同容量的中間密度材料比原有材料 的剛度低。在優化過程中( )x向0或1趨近減少密度中間值。 圖2.2表示不同懲罰因子下的材料密度函數曲線。 從上到下依次是1p =到5的密 度變化曲線。從圖上可以看出，懲罰因子越大，中間密度越趨向于0=（空洞） 。 通過這樣的處理使連續變量優化結果趨近于離散變量優化結果。但是過大的懲罰因 子會使拓撲優化求解過程出現震蕩，所以p值需要小心選擇。 圖 2.2 不同懲罰因子 p 下的密度函數曲線 SIMP法的最大問題是，設計結果不僅僅顯式的依賴于懲罰因子p，而且還受到 有限元網格的影響。只有通過增加約束等方法來改良函數性態才能獲得最優解。 Sigmund和Petersson27研究了存在于拓撲優化中的數值不穩定性問題，并介紹了多 種解決方法。對拓撲優化中數值不穩定性的研究放在第本章第二節。 華 中 科 技 大 學 碩 士 學 位 論 文 12 由于SIMP法具有思想簡單， 易于實現等優點， 它的應用范圍從單純的彈性體結 構拓撲優化拓展到了其它領域，例如MEMS、導電、傳熱和特征頻率等領域。本文 主要研究基于熱傳導的拓撲優化問題。 在傳熱學中， 材料的性質用導熱率來描述。 將(2.1)式中的彈性模量替換成熱傳 導系數就建立了材料導熱率與密度間的函數關系： 0 ( )( ) p xx= (2.3) 式中 0 是給定的導熱率。(2.3)式表示了SIMP法導熱材料數學模型。該模型簡 單，因此局限也大，只能描述各向同性材料。在下一章中，結合均勻化方法中的部 分理論，提出了改進的密度懲罰法模型，使該模型適用于各向異性材料的拓撲優化。 2.1.2 熱傳導模型熱傳導模型 許多與時間無關的物理系統，即穩態系統中的物理問題都可以用橢圓偏微分方 程來描述。比如彈性體結構拓撲優化求解的物理問題是彈性力學問題，本質上是求 解一個拉普拉斯方程。正如第一章所述，通過改變邊界條件和載荷，彈性力學拓撲 優化模型可以擴展到其它研究領域中，如導電、熱應力、熱傳導等問題的拓撲優化 中。因為它們的本質都是求解類似的偏微分方程。在研究熱傳導結構的拓撲優化問 題之前，首先研究熱傳導模型。 熱傳導過程滿足傅立葉定律： gradTq = - (2.4) q為熱流密度矢量，gradT是空間某點的溫度梯度向量，是熱傳導系數。對于 一維導熱問題，直接對(2.4)式積分即可求出熱量Q。對于多維熱傳導問題，其熱傳 導控制方程為： ()()() xyzvp TTTT qc xxyyyzt += (2.5) p c是導熱材料的比熱容)/(KkgJ，代表材料密度 3 /mkg， xyz ，分別是 材料沿xyz， ，方向的導熱系數()/Wm Ki， v q為物體的內熱源強度/W kg 熱傳導問題的邊界條件有三類： 華 中 科 技 大 學 碩 士 學 位 論 文 13 第一類Dirichlet邊界條件，給定邊界 1 上的溫度值： ( )()T xyztT t=， ， ， (2.6) 第二類Neumann邊界條件，給定邊界 2 上的熱流密度： ( ) xxyyyy TTT q t xyy += nnn (2.7) 第三類邊界條件為 3 上的混合邊界條件： () xxyyyy TTT h TT xyy += nnn (2.8) 圖 2.3表示 1 和 2 作為邊界條件的熱傳導區域： 1 2 圖 2.3 熱傳導區域示意圖 對于穩態熱傳導，溫度對時間的變化率為0： 0 T t = (2.9) 由于各向同性導熱材料各個方向導熱率相同 xyz =， 故公式(2.5)可以簡 化為： 222 222 ()0 v TTT q xyz += (2.10) 如果給定問題的初始條件： 華 中 科 技 大 學 碩 士 學 位 論 文 14 ()0()T xyztTxyz=， ， ， ， (2.11) 對于二維穩態熱傳導， 在滿足邊界條件(2.6)(2.8)以及初始條件(2.11)式的溫度場 使以下泛函I取極小值，即： () () 123 0 , 1 2d 2 min xyv BC T IC T TT Iq T xy =+ 2 2 (2.12) 通過變分求解就可以求出設計域中的溫度場。離散化(2.12)式，設計域即被劃 分為e個單元。為了方便網格劃分選擇矩形網格，每個單元被劃分為由四個節點組成 的正方形。單元可以表示為圖2.4的形式： x y o 圖 2.4 單元形式以及節點編號 要形成單元剛度矩陣，首先要選擇單元插值函數。4節點單元可以選擇雙線性插 值函數。通過插值函數，在節點溫度已知的情況下，單元內任意點的溫度可以表示 成節點溫度的插值形式： ()() 0123 , e Tx yaa xa ya xyf x y=+= (2.13) 在節點處有： (), iii T x yT= (1,2,.,4)i = (2.14) 把 1 T到 4 T代入(2.13)式中得到： 011111 122222 233333 344444 1 1 A 1 1 aTxyx y aTxyx y aTxyx y aTxyx y = e t (2.15) 通過(2.15)式可以求解出(2.13)中的待定系數 ： 1 A = e t (2.16) 華 中 科 技 大 學 碩 士 學 位 論 文 15 把溫度插值函數由(2.13)的坐標形式重寫為節點溫度的形式得到： ()() 1 ,( , ) eee Tx yx yx y = AttfN (2.17) ( , )N x y稱為單元形函數。 溫度沿x和y方向的導數可以表示為： () () 0 0123 1 2 0123 3 010 001 a T aa xa ya xy ay xx C Ta x aa xa ya xy yy a + = + (2.18) 將(2.16)式代入(2.18)式，得到： 1 ee T x T y = = C AttB (2.19) B為單元的幾何矩陣。 將(2.19)式代入(2.12)式，并求變分極值0 e I = t ，得到單元溫度的求解公式： eee tiK= Q (2.20) 其中 e K為單元導熱矩陣： d T T exy xxyy =+ NNNN K (2.21) 將 e K寫成矩陣的形式： d T e = BK (2.22) 為單元熱導率矩陣： 0 0 x y = (2.23) 把(2.22)式的面積分轉化為二次積分得到單元導熱矩陣的完整形式： 24 12 d d xy T e xy x y= BK (2.24) 華 中 科 技 大 學 碩 士 學 位 論 文 16 把單元導熱矩陣和右端載荷向量疊加以后得到求解溫度場的有限元格式： K=TQ (2.25) 式中 1 N e e KK = =是組裝后的總體傳熱矩陣， 1 T N e e= =t為疊加的總體節點溫 度向量， 1 N e e QQ = =為疊加后的總體右端向量。 (2.25)式是一個代數方程，給定初始以及邊界條件后，利用計算機可以方便地求 解出設計域的溫度場。 三維散熱結構拓撲優化中的單元傳熱矩陣構建與二維的類似，采用正6面體單 元，每個單元具有8節點，節點排列順序見圖2.5。 T2 T1 T4 T5 T6 T8 T7 T3 圖 2.5 正 8 面體單元節點排列 對三維正6面體8節點單元，選擇(2.26)式的插值函數： ()() 01234567 , , , e Tx y zaa xa ya za xya xza yza xyzf x y z=+= (2.26) 將8個節點坐標代入(2.26)式得到： 1111111 11 111 1 2222222222222 3333333333333 4444444444444 5555555555555 6666666666666 77 8 1 1 1 1 1 1 1 Txyzx yx zy zx y z Txyzx yx zy zx y z Txyzx yx zy zx y z Txyzx yx zy zx y z Txyzx yx zy zx y z Txyzx yx zy zx y z xyT T = e t 0 1 2 3 4 5 777777777776 888888888888 7 A 1 a a a a a a zx yx zy zx y za xyzx yx zy zx y za = (2.27) 華 中 科 技 大 學 碩 士 學 位 論 文 17 通過(2.15)式可以求解出(2.27)中的待定系數 。 溫度沿x、y和z方向導數可以表示為： () () () 01234567 01234567 01234567 01000 00100 00010 T aa xa ya za xya xza yza xyz xx T aa xa ya za xya xza yza xyz yy T aa xa ya za xya xza yza xyz zz yzyz xzxz xyxy + += + 0 1 2 3 4 5 6 7 a a a a C a a a a = (2.28) (2.24)式中的單元導熱矩陣變為三階矩陣： 00 00 00 x y z = (2.29) 將(2.22)的體積分化為三次積分： 468 111 d d d zxy T e zxy x y z= BK (2.30) 通過以上推導，得到了拓撲優化中三維溫度場的求解算法。在第3章中還將介 紹含各向異性導熱材料區域中的溫度場有限元求解算法。 2.1.3密度法優化模型密度法優化模型 最早引入拓撲優化研究的是彈性體結構最小化柔度問題。經過近三十年的發展， 此類優化問題的研究已經非常成熟。一些公司把此類拓撲優化模塊整合到了商業軟 件中。如Altair公司在Hyperworks中提供了Optistruct拓撲優化軟件包。 彈性體結構拓撲優化的數學模型可以表示為(2.31)式的形式： 華 中 科 技 大 學 碩 士 學 位 論 文 18 :( ) . .:( ) Minimize C xf udx stx dxV = = i (2.31) 其中f為力載荷，u為位移。目標函數( )C x為結構的平均柔度。約束( )x dxV = 代 表體積限制。優化的目標是在設計域中給定材料體積約束下，尋找最合適的材料分 布使得目標函數值最小。 為了方便用數值方法求解優化問題，需要把設計域離散成含n個單元的區域： 1,., 1 :() . .:,0,1 n n n e e Minimize C xx stxV x = = (2.32) 若nV,.,1，對于集合 n x1 , 0，滿足( )x dxV = 非空，則(2.31)式有解28， 但不能保證解的唯一性。式(2.32)是密度法拓撲優化的一般格式，引入帶懲罰因子的 密度函數(2.1)后得到SIMP模型如 (2.33)式： 1 ,., 1 :() . .:,0,1 pp n n n e e Minimize C xx stxV x = = (2.33) 拓撲優化中選擇合適的目標函數是問題的第一步，也是最為關鍵的一步。彈性 體結構優化中通常選擇結構的柔度作為目標函數。散熱結構研究的物理問題與彈性 體結構存在一定的相似性，目標函數可以通過類似的推導得到。Dems29, 30研究了邊 界問題中瞬態溫度場的靈敏度分析問題，為散熱結構拓撲優化提供了理論依據。其 思想可以推廣到熱傳導結構的拓撲優化中。一般散熱器的設計目標是要在給定的散 熱空間中保證工作器件的正常運行，即最高溫度最小。但是以最高溫度點作為設計 變量的函數，具有不連續性（最高溫度點所在的點的空間位置會隨著設計變量的變 化而改變） 。這會對求解造成困難19。可根據程新廣31, 32提出的最小熱傳遞勢容耗散 原理定義一個性態良好的目標函數。對于穩態熱傳導問題，傳入區域的熱量為： += s nds qqdQ (2.34) 從(2.34)式可知，熱量由兩部分組成。qd 代表區域內的熱源產生的熱量； n s q ds 項代表從邊界S出入區域的熱量。在非穩態對流情況下，將各點溫度與熱流密 華 中 科 技 大 學 碩 士 學 位 論 文 19 度的乘積積分稱為散熱弱度： n s CqTdq Tds =+ (2.35) 在實體內部和邊界上傳入的熱量一定的條件下，溫度越低，則散熱弱度越小， 整個設計域上的平均溫度越低。根據傅立葉定律把q的表達式代入(2.35)式，得到： () ijijij s iij TTT CKnhT TdsKd xxx =+ (2.36) 在對流換熱邊界上： )( fi i ij TTn x T K= (2.37) 把(2.37)式代入(2.36)式得到： 2 ij s ij TT ChT dsKd xx =+ (2.38) 對于穩態情況，不考慮對流邊界，將(2.38)式的前一項去掉得到： ij ij TT CKd xx = (2.39) (2.39)式也叫做熱量傳遞勢容，它的量綱是焦耳開爾文（ J K） ，它代表了物體 向周圍介質導熱的能力。隨著熱量從高溫傳向低溫，熱量傳遞中勢容存在耗散，當 其熱量傳遞勢容耗散為最小時，熱量傳遞效率最高。散熱結構拓撲優化的目標就是 使勢容耗散最小。經過離散化處理后，(2.39)式可以化成： 1 ( ) n T eee e