流熱力耦合的高性能結構拓撲優化設計方法

摘要：拓撲優化和增材制造技術的快速發展為高性能復雜裝備提供了高效的產品設計和制造方法。目前高性能結構拓撲優化只考慮熱力耦合或者流熱耦合的設計，且大多局限于簡單結構，未考慮流熱力三場共同作用下的設計，限制了結構性能的提升。針對流熱力多物理場工況下的高性能復雜結構設計這一挑戰，提出了一種流熱力耦合拓撲優化方法，以提高復雜結構的承溫能力。首先引入流場、溫度場和結構位移場的控制方程，對計算域的流固材料進行統一表征；然后以最小化平均溫度為目標，以流動能量耗散和結構柔度為約束，建立了流熱力耦合的拓撲優化模型，并結合變分法和拉格朗日函數開展了設計變量的靈敏度分析；最后將所建立的拓撲優化模型應用于渦輪的結構設計，得到了散熱性能良好、流道分布合理的可增材制造結構。

關鍵詞：拓撲優化；變密度法；多物理場；高性能結構

中圖分類號：TH122；TP391.7

DOI：10.3969/j.issn.1004132X.2024.03.011

0引言

高性能結構設計一直以來都是非常具有挑戰性的工作。作為一種先進的結構設計方法，拓撲優化（topologyoptimization，TO）因其無限的設計自由度和正向設計的理念而成為高性能結構的理想設計工具［1-2］。近年來，增材制造（additivemanufacturing，AM）因其可以實現傳統工藝難以生產的復雜幾何結構而引起了廣泛關注［3-4］。相比于車削、鍛造等傳統減材加工工藝，增材制造通過逐層打印的方式完成由計算機輔助設計（computeraideddesign，CAD）模型到三維結構的制造，極大地縮短了產品開發周期，減少了材料成本。拓撲優化與增材制造的結合為高性能復雜裝備提供了更高效的產品設計和制造方法，具有顯著的研究意義和社會價值［5］。

拓撲優化最初是針對結構力學問題而開發的，現已廣泛應用于由單一物理場控制的一系列工程問題，例如最小化結構柔度設計［6］、最佳熱傳導結構的設計［7-8］、流體流動設計［9］、波傳播、聲學設計［10-11］和生物醫學設計［12-13］等。在結構柔度設計方面，FERRARI等［14］給出了針對結構柔度拓撲優化的緊湊型MATLAB代碼，討論了加速策略的使用，大大減少了總體計算時間。WU等［15］提出一種針對shell-infill結構的柔度拓撲優化方法，該方法能夠生成同時具有外部殼體和內部填充材料的多孔結構。WANG等［16］結合柔度比例濾波器提出了一種新的密度變量更新方案，避免了獲取靈敏度信息時涉及的數值推導和計算問題，提高了收斂速度。KRIEGESMANN等［17］考慮載荷和空間變化的彈性模量的不確定性，提出了一種穩健的拓撲優化方法，該方法利用一階二階矩（firstordersecondmoment，FOSM）方法來估計柔度的平均值和方差，降低了計算成本。在熱力耦合結構設計方面，RODRIGUES等［18］首先使用均勻化方法建立了最小化結構柔度的熱彈性拓撲優化模型。YANG等［19］研究了機械力和熱負載共同作用下的動態柔度拓撲優化。高效散熱裝置［20］、特定溫度條件下的最大剛度結構［21］以及受溫度約束的結構［22］等設計先后被提出。TAKEZAWA等［23］設計了強度和熱傳導約束下的熱彈性結構。DEATON等［24］提出了一種基于應力的p范數聚合函數的熱結構設計。ZHENG等［25］考慮到材料屬性、載荷等相關的混合不確定性，實現了穩健性熱力耦合拓撲優化設計。

受限于流固耦合邊界的識別與表征以及流動問題的大規模計算量，涉及耦合三維流場的數值計算非常具有挑戰性。在流固耦合拓撲優化設計方面，YOON［26］研究了考慮流固邊界變化的拓撲優化問題，將結構平衡方程和流體N-S方程同時擴展到了流固域，避免了邊界處數據的傳輸及迭代求解。PICELLI等［27］提出了一種針對流固耦合系統頻率響應的拓撲優化方法，該方法可應用于汽車和飛機乘客艙的降噪以及水下結構的振動控制等。LI等［28］使用反應擴散方程來更新水平集函數，并利用貼體自適應網格進行流固邊界捕捉，研究了三維流固耦合問題。在流熱耦合設計方面，YOON等［29］在N-S方程中引入了摩擦力項，并忽略了流體內部的浮力和粘黏耗散，以最小化熱順應性為目標設計了利用強制對流傳熱的二維散熱結構。YU等［30］和YAJI等［31］分別使用移動組件法（movingmorphablecomponent，MMC）和水平集法研究了流動傳熱拓撲優化問題。ZHAO等［32］研究了考慮傳熱和結構承載的冷卻直管二維截面設計問題，并提出了一種由固體、液體和孔洞組成的多相材料拓撲優化方法。

上述高性能結構拓撲優化工作只考慮熱力耦合、流固耦合及流熱耦合的設計，且大多研究都局限于簡單幾何模型下的小規模問題，未考慮流熱力三場共同作用下的設計。然而，工程上以渦輪機械為代表的很多復雜裝備處于流動傳熱結構應力共同作用下的工況，因此，為實現此類復雜裝備的高性能結構設計，研究流熱力耦合的拓撲優化更具有實際意義。本文基于變密度法，提出一種高性能結構流熱力耦合拓撲優化方法，以提升復雜結構的性能。首先采用Brinkman懲罰方式對流體域和固體域進行建模，引入流場、溫度場和結構位移場的控制方程。其次以平均溫度為目標，流動能量耗散和結構柔度為約束，構建了流熱力耦合拓撲優化模型。為實現單元設計變量的靈敏度計算，通過引入弱形式的拉格朗日函數和變分法推導了連續形式的伴隨方程。使用基于有限體積法離散的OpenFOAM軟件對控制方程和伴隨方程進行求解，拓撲優化問題的求解策略為移動漸近線法（methodofmovingasymptotes，MMA），并利用兩種數值過濾算法穩定優化歷程。最后將建立的拓撲優化模型應用于渦輪的結構設計，得到了散熱性能良好、流道分布合理的可增材制造結構。需要強調的是，本文的研究著重于復雜結構流熱力耦合拓撲優化算法的理論及其實現，具體的增材制造過程及可制造性等方面的內容不在本文的研究范疇。

1流熱力多物理場控制方程

考慮圖1所示的多物理場問題模型，整體計算域Ω由固體域Ωs和流體域Ωf組成，計算域邊界包含流體入口及出口邊界Γin、Γout，固定溫度邊界ΓT0，熱通量邊界Γq，固定位移邊界ΓD0，應力邊界Γf，其余邊界為應力自由的絕熱邊界。該問題中的各物理場分別由N-S方程、能量方程和線彈性平衡方程進行描述，各物理場之間的耦合情況如下：

（1）流熱耦合。主要考慮流場對溫度場的影響，即流體的速度變化通過對流和擴散過程影響溫度場的分布，忽略溫度場對流速場的影響。

（2）流力（流固）耦合。主要考慮設計域中固體材料部分對流速場的阻礙作用，忽略流速場引起的結構變形。

（3）流熱力耦合。主要考慮彈性變形引起的結構應力，忽略溫度變化產生的熱應力。流固材料屬性均為各向同性，且不隨溫度發生變化。

1.1流動控制方程

本文將流體視為穩態不可壓縮充分發展的層流狀態，流動控制方程由N-S方程、連續性方程及相應邊界條件表示：

1.2能量方程

計算域的溫度場由能量方程進行控制，本文采用穩態熱傳導方程作為能量控制方程，能量方程通過對流項與N-S方程進行弱耦合。能量方程及其邊界條件定義如下：

1.3線彈性平衡方程

計算域固體材料的位移場由線彈性平衡方程進行控制，張量形式的線彈性平衡可表示為

其中，D代表結構位移場；fb為固體體積力項；I為單位矩陣；在固定位移邊界上指定位移D0；在應力邊界上指定面力fN；其余邊界為自由邊界。μ和λ為材料的拉梅常數，它們與彈性模量E和泊松比ν的關系分別為

2流熱力拓撲優化模型

2.1優化模型

本文使用基于密度法的框架進行結構的流熱力耦合拓撲優化設計，其基本思想是在設計域Ω內用設計變量場γ來描述域內流固兩相材料的分布：γ=1代表流體材料，γ=0代表固體材料。在優化過程中，兩相材料的物理屬性由標量場γ插值描述，從而將設計變量引入多物理場控制方程。

散熱性能和降溫效果通常是影響受熱結構綜合性能的重要因素，是本文優化的主要目標。本文以結構的平均溫度作為目標函數來衡量降溫散熱效果，平均溫度目標函數Ψ的定義如下：

Ψ=∫ΩTdΩ/∫ΩdΩ（6）

一般情況下，作為冷卻介質的流體通過冷卻通道的壓力損失越大表明結構的散熱性能越好，但實際上流體的流動損失受到流動效率要求的限制。此外，在涉及流場的優化設計時若不進行流場能量耗散的約束往往會產生大量懸浮固體，不具備制造條件，因此，本文對流場的能量耗散進行約束，流動能量耗散約束函數定義如下：

設計域內流體材料的分布在提高結構散熱性能的同時必定會影響結構的剛度，結構剛度不足會導致結構承載性能的降低，因此，在優化時引入結構柔度約束函數以保證結構剛度。柔度約束函數定義如下：

式中，為流體材料體積分數上限；為流動能量耗散約束上限；為結構柔度約束上限。

2.2材料插值方案

在基于達西滲流理論的流體拓撲優化框架下，BORRVALL等［9］首先采取基于有理數近似的插值策略（rationalapproximationofmaterialproperties，RAMP），在流體拓撲優化問題中得到了很好的應用。因此，本文采取RAMP策略對各物理場屬性進行插值描述。

流體和固體材料的流動屬性通過對式（2）中的材料流阻系數插值進行表征，插值格式為

其中，αmax為材料最大流阻系數，即固體材料的流阻系數；qα為控制懲罰強度的常數，本文取0.01；當設計變量γ=1時，α（γ）=0，流動控制方程退化為N-S方程；當γ=0時，α（γ）=αmax，此時流動阻力達到最大值，流速趨近于0，表示此處為固體材料。

流體和固體材料的傳熱屬性通過對材料熱導率插值進行表征，插值格式為

式中，kf、ks分別為流體區域和固體區域的熱導率；qk的作用與qα一致，本文取qk=0.01。

流體和固體材料的力學屬性通過對材料的彈性模量插值進行表征，插值格式為

其中，Emax表示固體材料的彈性模量，Emin表示流體材料即孔洞的彈性模量，通常取Emin=0.001來避免數值求解過程中出現的矩陣奇異問題。qE的作用與qα和qk一致，本文取qE=0.01。

3靈敏度計算方法

通過計算目標函數對設計變量的靈敏度從而更新設計變量的分布是求解拓撲優化問題的主要求解策略。由于本文設計變量γ的變化會引起狀態場量（速度、溫度、位移）的變化，而狀態場量對設計變量的導數隱含在控制方程中無法直接獲得，這給靈敏度分析帶來了困難。針對這一難點，本文采用基于拉格朗日乘子法和連續伴隨法的思想來求設計變量的靈敏度。與常見的離散伴隨法敏度分析相比，連續伴隨法需要手動推導連續形式的伴隨方程，進而求解伴隨方程獲得伴隨變量，最終靈敏度可由伴隨變量經過簡單的計算獲得。

首先引入弱形式下的拉格朗日方程，即在整個求解域上將狀態方程與拉格朗日乘子（伴隨變量）相乘后增廣至目標函數（或約束函數）上。拉格朗日函數L定義如下：

其中，Φ為待求靈敏度的目標或約束函數；θ為拉格朗日乘子，即伴隨變量；u^、p^、T^、D^分別為伴隨速度、伴隨壓力、伴隨溫度和伴隨位移；R代表狀態控制方程，即式（1）～式（4）。狀態控制方程始終滿足R≡0，因此對于某個目標或者約束函數來說始終有L≡Φ，求解目標或者約束函數的靈敏度即可轉換為求解式（13）的靈敏度。令式（13）對狀態變量的變分為0可求得一組伴隨變量。通過引入伴隨變量，靈敏度分析時只需對式（13）中與設計變量相關的項直接求導，而不必考慮其余狀態量的中間效應。

3.1能量耗散靈敏度

求能量耗散約束函數J的靈敏度時，式（13）中Φ=J。由于能量方程與流動方程是解耦求解的，且能量耗散約束與變量無關，故能量耗散約束的伴隨變量θ=（u^，p^，0，0）。為獲得能量耗散約束對設計變量的靈敏度，需引入拉格朗日函數的總變分：

為使任意擾動量δp和δu下，上式均成立，必須滿足其中的積分項為0。令上式中的體積分項為0，得到伴隨方程；令上式中的邊界積分項為0，得到伴隨邊界條件。最終得到的伴隨N-S方程及邊界條件：

需要指出的是，在求解伴隨方程的過程中認為原始方程的量是已知的，即在求解伴隨變量（u^，p^）時，原始場量（u，p）是已知量。至此，將式（21）的解u^i代入式（16）即可得到能量耗散約束的靈敏度。

3.2平均溫度靈敏度

目標函數平均溫度Ψ的靈敏度推導過程與3.1節一致，需要注意的是該函數同時與u、p、T相關，因此需要求解一組伴隨N-S方程和伴隨能量方程以獲得伴隨變量（u^，p^，T^）。平均溫度的靈敏度表達式為

3.3柔度靈敏度

柔度約束函數C的靈敏度表達式為

4優化求解策略

4.1優化算法流程

本文使用基于有限體積法離散的OpenFOAM軟件作為多物理場問題和拓撲優化問題的求解器，拓撲優化求解算法采用移動漸近線法。優化算法的流程圖見圖2，具體執行步驟如下：

（1）對設計變量進行初始化，根據式（10）～式（12）更新材料的流阻系數、熱導率和彈性模量。

（2）分別求解式（1）、式（3）、式（4）得到當前迭代步下的流速場u、溫度場T和位移場D，進而計算式（6）～式（8）得到當前迭代步下的目標函數和各約束函數的值，在求解溫度場T時使用了流速場u的數據。

（3）求解式（21）得到能量耗散約束函數的伴隨變量，求解式（23）、式（24）得到平均溫度目標函數的伴隨變量。根據式（16）、式（22）和式（25）分別計算能量耗散約束、平均溫度目標和結構柔度約束的靈敏度。

（4）根據得到的靈敏度信息調用MMA算法更新設計變量。

（5）判斷當前迭代步下設計變量是否完成收斂或迭代次數是否達到上限，若滿足條件則拓撲優化流程結束，否則使用當前迭代步更新的設計變量重新更新材料屬性，繼續迭代求解直至滿足程序結束條件。

4.2數值過濾方法

基于密度法的拓撲優化流程一般存在棋盤格效應、中間密度單元過多等數值不穩定現象。為抑制上述數值不穩定現象，本文使用基于Helmholtz型偏微分方程的PDE過濾［35］和基于雙曲正切函數的投影［36］對密度場進行過濾，經兩層過濾后可得到邊界清晰的拓撲結果。

首先對設計變量場γ使用PDE過濾以消除棋盤格效應，PDE過濾函數如下：

5數值算例

以燃氣輪機為代表的渦輪機械是典型的多物理場工況作用下的復雜結構。本文以渦輪盤體和渦輪葉片作為研究對象，采用提出的流熱力拓撲優化模型對其進行結構設計。冷卻氣流由渦輪軸心位置處引入渦輪內部，通過內部冷卻通道截面的對流換熱帶走溫度，從而提高渦輪的承溫能力。渦輪材料及冷卻空氣的相關屬性見表1。

5.1渦輪盤體設計

給定渦輪盤體幾何模型如圖5所示，尺寸單

位為mm。為減少計算量，取渦輪盤體的一個旋轉周期作為計算域，渦輪盤體計算域及邊界情況如圖6a所示。為保證優化設計時空氣流動能夠充分發展，入口和出口處都有一段延長導管。中間區域為拓撲優化設計域，入口和出口處為純流動狀態的流體材料。計算域離散為6×105個六面體網格，計算時采用24核并行計算。

計算模型邊界條件設置如下：指定冷卻通道入口處的法向速度為0.05m/s，溫度為20℃，冷卻通道出口處為零壓力絕熱邊界，其余壁面為絕熱無滑移邊界。熱載荷由體熱源引入，給定熱源密度為10kW/m3；力載荷由葉片表面均布壓力引入，給定壓力大小為0.1MPa，同時固定轉軸處位移為0。

為了得到流熱力耦合拓撲優化各項約束的參考值，首先對圖6b所示參考設計進行模擬。該結構的流體材料體積分數為1.5%，經仿真分析，圖6b所示結構的能量耗散J0=1.0×10-5kg·m2/s3，結構柔度C0=6×10-10J，平均溫度為130.7℃。為提高結構的承溫能力，在參考設計的基礎上，設置許用流體材料分數為10%，能量耗散約束為=1.5J0，結構柔度約束為=1.1C0，設計變量的初始值為0.5。拓撲優化收斂準則設置為最大迭代次數400或設計變量的增量‖γk+1－γk‖lt;10－3。

由于設計流場的多物理場問題存在極強的非線性，為保證數值計算穩定，將初始最大流阻系數αmax設置為108，之后隨迭代次數的增加而逐漸增大至1011。同時，在優化前期對能量耗散約束和結構柔度約束也進行放松：在初始迭代階段，設置能量耗散約束上限為10J0，結構柔度約束上限為2C0。隨著優化過程的進行，控制這兩個約束上限分別線性減小至和。渦輪盤體的拓撲優化迭代過程如圖7所示，目標及約束迭代曲線如圖8所示。

在優化前期，流阻系數較小，能量耗散值也較小；同時由于流體材料的大面積存在，結構柔度較高；由于此時對固體材料的流速懲罰不足，固體材料內部仍存在不合實際的對流換熱，故結構的平均溫度較低。隨著優化過程的進行，流阻系數逐漸增大，能量耗散約束上限和結構柔度約束上限逐漸減小，流體材料由分散狀態逐漸收斂為完整的管道；同時固體材料內部的流速逐漸懲罰為零，對流換熱只能發生在流體管道的表面，而且固體材料的質量熱容遠小于流體材料的質量熱容，因此結構的平均溫度有所升高。優化迭代至300步左右時，設計域域內流體形成了貫通入口和出口的管道，此時管道內流體形成連續的流動，流場能量耗散瞬時減小。最終結構拓撲變化逐漸穩定，平均溫度為120℃。圖9為盤體優化前后的溫度場云圖，可以看出優化后的溫度場整體上低于參考設計值，渦輪盤體最終設計結果如圖10所示。表2給出了渦輪盤體優化前后的性能對比，與參考設計相比，拓撲優化后結構的平均溫度降低了8.2%。

5.2渦輪葉片設計

渦輪葉片是將燃氣熱能轉化為渦輪旋轉動能的核心部件，其工作溫度已接近甚至超過了材料的熔點，高效的渦輪葉片冷卻通道是保證渦輪穩定工作的關鍵。本文對圖11a所示渦輪葉片模型進行設計，該模型含有一個入口和一個出口，計算域離散為3×105個六面體網格，計算時采用24核并行計算。

渦輪葉片計算模型邊界條件設置如下：指定冷卻通道入口處的法向速度為0.05m/s，溫度為20℃，冷卻通道出口處為零壓力絕熱邊界，其余壁面為絕熱無滑移邊界。熱載荷由體熱源引入，給定熱源密度為10kW/m3；力載荷由葉片表面均布壓力引入，給定壓力大小為0.1MPa，同時固定葉根處位移為0。

首先對圖11b所示參考設計進行模擬，以獲得流熱力耦合拓撲優化各項約束的參考值。參考設計的平均溫度為69.7℃，流體材料體積分數為1%，流動能量耗散J0=5×10-7kg·m2/s3，結構柔度C0=2×10-10J。由于渦輪葉片處受熱較多且承受的力學載荷相對較少，故可以適當增加流體材料占比。在參考設計的基礎上設置許用流體材料分數為40%，能量耗散約束為=12J0，結構柔度約束為=1.7C0，設計變量的初始值為0.5。拓撲優化收斂準則設置為最大迭代次數400或設計變量的增量‖γk+1－γk‖lt;10－3。本算例采取與5.1節一致的放松策略，將初始最大流阻系數αmax設置為108，之后隨迭代次數的增加而逐漸增大至1011。同時，在初始迭代階段，

設置能量耗散約束上限為25J0，結構柔度約束上限為3C0。隨著優化過程的進行，控制這兩個約束上限分別線性減小至和。

圖12和圖13給出了渦輪葉片優化設計的迭代過程和目標及約束函數的迭代曲線。從優化結果可以看出，當流體材料許用體積和流動能量耗散許用值增加后，生成的最優構型更加復雜。由于大量增加了流體材料用量，即使在優化初期對能量耗散和結構柔度約束進行了大量放松，仍難以找到穩定的優化方向，優化前期出現了少許振蕩。隨著流阻系數的增加，能量耗散逐步增大；能量耗散約束上限及結構柔度約束上限隨迭代步數的增加逐漸縮緊到許用值，最終在許用能量耗散和結構柔度約束下，平均溫度目標收斂至54.5℃。圖14給出了葉片優化前后的溫度場分布情況，可以看出優化后葉片的溫度明顯低于參考設計，葉片流道截面形狀如圖15所示。表3給出了葉片優化前后的性能對比，與參考設計相比，拓撲優化后葉片結構的平均溫度降低了21.5%。

6結論

（1）本文針對流熱力多物理場工況下的高性能復雜結構設計這一挑戰，提出了一種復雜結構流熱力耦合拓撲優化方法，為高性能結構的設計和制造提供了新的方案。

（2）基于變分法和拉格朗日乘子法建立了多物理場拓撲優化問題的連續伴隨敏度分析框架，推導了伴隨狀態方程及其邊界條件。使用移動漸近線法對拓撲優化問題進行求解，結合Helmholtz型偏微分方程過濾和雙曲正切方程投影，改善了拓撲優化過程的不穩定現象，得到邊界清晰的可增材制造拓撲結構。

（3）將所提的流熱力拓撲優化模型應用于實際工程裝備的設計，分別對渦輪盤體和渦輪葉片進行結構設計，得到了散熱性能良好、流道分布合理的高性能結構。通過比較發現，本文方法可以獲得比參考設計承溫能力更好的結構。因此，考慮流熱力耦合的拓撲優化設計有其實際必要性。

參考文獻：

［1］LIUJikai，GAYNORAT，CHENShikui，etal.CurrentandFutureTrendsinTopologyOptimizationforAdditiveManufacturing［J］.StructuralandMultidisciplinaryOptimization，2018，57（6）：2457-2483.

［2］劉成穎，譚鋒，王立平，等.面向機床整機動態性能的立柱結構優化設計研究［J］.機械工程學報，2016，52（3）：161-168.

LIUChengying，TANFeng，WANGLiping，etal.ResearchonOptimizationofColumnStructureDesignforDynamicPerformanceofMachineTool［J］.JournalofMechanicalEngineering，2016，52（3）：161-168.

［3］DUPLESSISA，BROECKHOVENC，YADROITSAVAI，etal.BeautifulandFunctional：aReviewofBiomimeticDesigninAdditiveManufacturing［J］.AdditiveManufacturing，2019，27：408-427.

［4］劉偉，李能，周標，等.復雜結構與高性能材料增材制造技術進展［J］.機械工程學報，2019，55（20）：128-151.

LIUWei，LINeng，ZHOUBiao，etal.ProgressinAdditiveManufacturingonComplexStructuresandHigh-performanceMaterials［J］.JournalofMechanicalEngineering，2019，55（20）：128-151.

［5］朱繼宏，周涵，王創，等.面向增材制造的拓撲優化技術發展現狀與未來［J］.航空制造技術，2020，63（10）：24-38.

ZHUJihong，ZHOUHan，WANGChuang，etal.StatusandFutureofTopologyOptimizationforAdditiveManufacturing［J］.AeronauticalManufacturingTechnology，2020，63（10）：24-38.

［6］SIGMUNDO，MAUTEK.TopologyOptimizationApproaches［J］.StructuralandMultidisciplinaryOptimization，2013，48（6）：1031-1055.

［7］IGAA，NISHIWAKIS，IZUIK，etal.TopologyOptimizationforThermalConductorsConsideringDesign-dependentEffects，IncludingHeatConductionandConvection［J］.InternationalJournalofHeatandMassTransfer，2009，52（11/12）：2721-2732.

［8］SUBRAMANIAMV，DBOUKT，HARIONJL.TopologyOptimizationofConductiveHeatTransferDevices：anExperimentalInvestigation［J］.AppliedThermalEngineering，2018，131：390-411.

［9］BORRVALLT，PETERSSONJ.TopologyOptimizationofFluidsinStokesFlow［J］.InternationalJournalforNumericalMethodsinFluids，2003，41（1）：77-107.

［10］DHRINGMB，JENSENJS，SIGMUNDO.AcousticDesignbyTopologyOptimization［J］.JournalofSoundandVibration，2008，317（3/4/5）：557-575.

［11］RUPPCJ，EVGRAFOVA，MAUTEK，etal.DesignofPhononicMaterials/StructuresforSurfaceWaveDevicesUsingTopologyOptimization［J］.StructuralandMultidisciplinaryOptimization，2007，34（2）：111-121.

［12］SUTRADHARA，PAULINOGH，MILLERMJ，etal.TopologicalOptimizationforDesigningPatient-specificLargeCraniofacialSegmentalBoneReplacements［J］.ProceedingsoftheNationalAcademyofSciencesoftheUnitedStatesofAmerica，2010，107（30）：13222-13227.

［13］PARKJ，LEEDamian，SUTRADHARA.TopologyOptimizationofFixedCompleteDentureFramework［J］.InternationalJournalforNumericalMethodsinBiomedicalEngineering，2019，35（6）：e3193.

［14］FERRARIF，SIGMUNDO.ANewGeneration99LineMatlabCodeforComplianceTopologyOptimizationandItsExtensionto3D［J］.StructuralandMultidisciplinaryOptimization，2020，62（4）：2211-2228.

［15］WUJun，CLAUSENA，SIGMUNDO.MinimumComplianceTopologyOptimizationofShell–infillCompositesforAdditiveManufacturing［J］.ComputerMethodsinAppliedMechanicsandEngineering，2017，326：358-375.

［16］WANGHui，CHENGWenming，DURun，etal.ImprovedProportionalTopologyOptimizationAlgorithmforSolvingMinimumComplianceProblem［J］.StructuralandMultidisciplinaryOptimization，2020，62（2）：475-493.

［17］KRIEGESMANNB，LDEKERJK.RobustComplianceTopologyOptimizationUsingtheFirst-orderSecond-momentMethod［J］.StructuralandMultidisciplinaryOptimization，2019，60（1）：269-286.

［18］RODRIGUESH，FERNANDESP.AMaterialBasedModelforTopologyOptimizationofThermoelasticStructures［J］.InternationalJournalforNumericalMethodsinEngineering，1995，38（12）：1951-1965.

［19］YANGXiongwei，LIYueming.TopologyOptimizationtoMinimizetheDynamicComplianceofaBi-materialPlateinaThermalEnvironment［J］.StructuralandMultidisciplinaryOptimization，2013，47（3）：399-408.

［20］PARKJ，NGUYENTH，SHAHJJ，etal.ConceptualDesignofEfficientHeatConductorsUsingMulti-materialTopologyOptimization［J］.EngineeringOptimization，2019，51（5）：796-814.

［21］DEATONJD，GRANDHIRV.StiffeningofRestrainedThermalStructuresviaTopologyOptimization［J］.StructuralandMultidisciplinaryOptimization，2013，48（4）：731-745.

［22］ZHUXuefeng，ZHAOChao，WANGXuan，etal.Temperature-constrainedTopologyOptimizationofThermo-mechanicalCoupledProblems［J］.EngineeringOptimization，2019，51（10）：1687-1709.

［23］TAKEZAWAA，YOONGH，JEONGSH，etal.StructuralTopologyOptimizationwithStrengthandHeatConductionConstraints［J］.ComputerMethodsinAppliedMechanicsandEngineering，2014，276：341-361.

［24］DEATONJD，GRANDHIRV.Stress-basedDesignofThermalStructuresviaTopologyOptimization［J］.StructuralandMultidisciplinaryOptimization，2016，53（2）：253-270.

［25］ZHENGJing，CHENHong，JIANGChao.RobustTopologyOptimizationforStructuresunderThermo-mechanicalLoadingsConsideringHybridUncertainties［J］.StructuralandMultidisciplinaryOptimization，2022，65（1）：29.

［26］YOONGH.TopologyOptimizationforStationaryFluid-structureInteractionProblemsUsingaNewMonolithicFormulation［J］.InternationalJournalforNumericalMethodsinEngineering，2010，82（5）：591-616.

［27］PICELLIR，VICENTEWM，PAVANELLOR，etal.EvolutionaryTopologyOptimizationforNaturalFrequencyMaximizationProblemsConsideringAcoustic-structureInteraction［J］.FiniteElementsinAnalysisandDesign，2015，106（C）：56-64.

［28］LIHao，KONDOHT，JOLIVETP，etal.Three-dimensionalTopologyOptimizationofaFluid–StructureSystemUsingBody-fittedMeshAdaptionBasedontheLevel-setMethod［J］.AppliedMathematicalModelling，2022，101：276-308.

［29］YOONGH.TopologicalDesignofHeatDissipatingStructurewithForcedConvectiveHeatTransfer［J］.JournalofMechanicalScienceandTechnology，2010，24（6）：1225-1233.

［30］YUMinghao，RUANShilun，WANGXinyu，etal.TopologyOptimizationofThermal–fluidProblemUsingtheMMC-basedApproach［J］.StructuralandMultidisciplinaryOptimization，2019，60（1）：151-165.

［31］YAJIK，YAMADAT，KUBOS，etal.ATopologyOptimizationMethodforaCoupledThermal–fluidProblemUsingLevelSetBoundaryExpressions［J］.InternationalJournalofHeatandMassTransfer，2015，81：878-888.

［32］ZHAOXi，ZHOUMingdong，LIUYichang，etal.TopologyOptimizationofChannelCoolingStructuresConsideringThermomechanicalBehavior［J］.StructuralandMultidisciplinaryOptimization，2019，59（2）：613-632.

［33］PANTONRL.IncompressibleFlow［M］.Hoboken：Wiley，2013.

［34］于明豪.液冷散熱器流道與熱源分布拓撲優化研究［D］.大連：大連理工大學，2021.

YUMinghao.ResearchonTopologyOptimizationofCoolingChannelsandHeatSourcesDistributionforLiquid-cooledHeatSink［D］.Dalian：DalianUniversityofTechnology，2021.

［35］LAZAROVBS，SIGMUNDO.FiltersinTopologyOptimizationBasedonHelmholtz-typeDifferentialEquations［J］.InternationalJournalforNumericalMethodsinEngineering，2011，86（6）：765-781.

［36］XUShengli，CAIYuanwu，CHENGGengdong.VolumePreservingNonlinearDensityFilterBasedonHeavisideFunctions［J］.StructuralandMultidisciplinaryOptimization，2010，41（4）：495-505.