采用有限體積法的自然對流換熱拓撲優化數值方法

摘要："針對基于有限體積法的熱-流耦合數值分析缺乏拓撲優化框架問題，建立了采用有限體積法的自然對流換熱拓撲優化方法。建立了流體固體統一描述的自然對流換熱數學模型，提出了基于有限體積法的自然對流求解方法；通過引入基于結構的離散化描述方法，搭建了結構的描述框架并構建了流體與固體間的材料屬性插值模型；針對自然對流中增強散熱能力與增強物質傳輸能力問題，提出了拓撲優化的數學模型，明確了設計變量、目標函數以及約束條件，建立了自然對流拓撲優化方法框架；通過3個算例，驗證了提出的數學模型和求解方法，并與商業軟件COMSOL進行了對比。結果表明：提出的自然對流求解方法與COMSOL計算結果相比，壓力、速度、溫度的數值誤差均小于3%；對于增強散熱能力問題，采用提出的拓撲優化方法獲得的最高溫度優于COMSOL結果。所提基于有限體積法的自然對流拓撲優化方法的可行性和有效性得到了驗證。

關鍵詞："自然對流；拓撲優化；有限體積法；強化散熱；流體固體統一描述

中圖分類號："TH122"文獻標志碼：A

DOI："10·7652/xjtuxb202408011"文章編號：0253-987X（2024）08-0103-11

A Numerical Method for Topology Optimization of Natural Convection

Heat Transfer Based on Finite Volume Method

DU Fei1， TIAN Zhenxiong1， LIU Honglei"2，3， GUO Shuzhe3， GUO Junkang3， LI Baotong3

（1. School of Astronautics， Northwestern Polytechnical University， Xi’an 710072， China;

2. State Key Laboratory of Fluid Power and Mechatronic Systems， Zhejiang University， Hangzhou 310027， China;

3. School of Mechanical Engineering， Xi’an Jiaotong University， Xi’an 710049， China）

Abstract："A topology optimization method for natural convection heat transfer based on the finite volume method （FVM） is established to address the lack of a topology optimization framework in numerical analysis of thermal-fluid coupling based on the FVM. A fluid-solid unified mathematical model of natural convection heat transfer is first established， and a solution algorithm based on the FVM is proposed. By introducing a structure-based discrete description method， a structure description framework and a material attribute interpolation model between the fluid and the solid are constructed. In addition， mathematical models of topology optimization for heat dissipation enhancement and material transport enhancement are proposed， and the design variables， objective functions and constraints are clarified. A framework of natural convection topology optimization algorithm is established. Finally， three examples are given to verify the proposed mathematical models and algorithm， and a comparison is made with the commercial software COMSOL. The results show that the numerical differences of pressure， speed and temperature between the natural convection algorithm proposed in this paper and COMSOL are less than 3%. For enhancing heat dissipation capacity， the maximum temperature obtained by the proposed topology optimization method is superior to the COMSOL optimization result. The results verify the feasibility and effectiveness of the proposed natural convection topology optimization method based on the FVM.

Keywords："natural convection; topology optimization; finite volume method; strengthen heat dissipation; fluid solid unified description

自然對流換熱以流體中的溫度梯度或密度梯度作為驅動力驅使流體對流換熱，散熱過程無需外圍設備，故此類散熱器結構簡單、成本低、可靠性高、壽命長且無噪聲。自然對流廣泛應用在核反應堆非動能循環、電子設備散熱等方面"［1］。以核反應堆為例，自然對流是一種可靠的非動能循環方法，在核反應堆主泵停止工作時，僅依靠冷熱段密度差驅動冷卻劑循環，導出堆芯余熱，從而保證反應堆的安全。對于核潛艇反應堆而言，核反應堆依靠自然對流實現冷卻水循環能力是衡量核潛艇動力性能的重要指標，是實現固有安全性的關鍵"［2］。

隨著核反應堆、電子設備等對于性能要求的不斷提高，對于自然對流換熱能力的要求也不斷提高。自然對流換熱能力已成為設備性能提升的瓶頸。提高自然對流換熱能力，除了研制更高導熱系數的散熱材料外，還在于設計更加高效的散熱器結構。在傳統自然對流散熱結構"［3］的設計過程中，設計人員的工程經驗往往起到了主導作用，這種設計缺乏系統的熱力學分析和優化理論的指導。在當前機電設備尺寸小型化以及功率大型化的發展趨勢下，這種依靠經驗的設計模式，顯得越來越力不從心。拓撲優化是在結構中材料布局形式未知的情況下，依據對物理場的評估與分析，尋求結構中材料的最優分布形式。自然對流換熱器的設計問題，本質上是最佳材料分布問題，故拓撲優化是實現換熱器設計的有效途徑。

自然對流換熱的數值模擬，是對換熱器進行優化的基礎。陽祥和陶文銓"［4］采用有限差分法對高瑞利數下二維封閉方腔內自然對流進行了數值模擬，揭示了不同物質發生的流動與換熱特性。隨后，丁昊和陶文銓"［5］采用類似的離散方法研究了不同瑞利數下液態鉛鉍合金的自然對流換熱特性。田東東等"［6］采用實驗的方法研究了不同厚度金屬泡沫銅對換熱和對流強度的影響。有限體積法因其存儲空間小、穩定性強、便于推廣到三維以及能適用于復雜的求解區域，被廣泛應用于流體流動、傳熱的研究。針對自然對流換熱問題，解巖等"［7］提出了一種對流項離散格式修正方法，用于非結構網絡劃分。

為借助拓撲優化方法進行自然對流換熱器的設計，首先需建立換熱器對流換熱過程的物理模型。由于對流換熱問題的復雜性，在工程中常使用牛頓冷卻定律（Newton’s law of cooling，NLC）將自然對流問題簡化為導熱問題，并賦予對流邊界條件，通過選擇合適的對流系數來近似的模擬自然對流換熱過程"［8］。Burns等"［9］使用具有恒定對流系數的模型，對散熱器進行拓撲優化設計，并得到了有著樹狀分叉的散熱器結構。然而，基于NLC 模型的自然對流拓撲優化無法將流體的流動納入分析中，而只能用對流系數來衡量其對換熱器的影響，這樣的設計結果并不能完全的反應流體的流動與換熱特性。

針對上述問題，Alexandersen等"［10］于2014 年首次提出了使用全階自然對流模型進行的散熱器的拓撲優化設計。該模型基于不可壓縮的穩態層流假設，通過耦合 Navier-Stokes 方程（動量守恒）與對流-擴散方程（能量守恒）來完成對換熱及對流過程的描述。其中，耦合的方程使用有限元方法（finite element method，FEM）離散并使用變密度法（SIMP）來進行優化。Alexdandersen等"［11］進一步提出基于變密度法的散熱器設計方法，在該方法基礎上，設計了用于大功率LED 燈的新型散熱器"［12］。類似的，Li等"［13］針對自然對流散熱器設計，采用拓撲優化方法研究了加熱功率、熱源尺寸、材料導熱系數等對于散熱器設計的影響規律。李含靈等"［14］采用基于變密度的自然對流拓撲優化方法，用于三維散熱齒的優化設計。Coffin 和Maute"［15］基于水平集方法進行二維和三維自然對流散熱器的拓撲優化設計，使用全階的Navier-Stokes 方程，其物理模型使用拓展有限元完成求解工作。吳璇和陳群"［16］引入了Brinkman摩擦阻力項，利用變分法推導出自然對流過程的歐拉方程，得到了不動點迭代方法。

在傳統的力學結構拓撲優化領域，有限元方法是被廣泛使用的數值求解方法。當問題拓展至自然對流拓撲優化時，這種慣性也得以保留。在上述關于自然對流散熱器的拓撲優化設計中，幾乎都采用了有限單元法作為數值求解方法。另外一方面，在流體數值分析領域，特別是自然對流問題，作為一種熱流雙向強耦合的復雜問題，其求解難度較高。有限體積法因其自身優勢而得到了廣泛應用，Bari等"［17］根據大量調研指出，目前有限體積法（FVM）是自然對流問題分析中最常用的數值方法。然而，目前基于FVM進行自然對流拓撲優化的工作還較少。對于基于FVM的熱-流耦合拓撲優化問題，仍存在諸多問題，特別是缺少基本的優化框架。

針對上述問題，本文旨在建立基于FVM的自然對流換熱器的拓撲優化設計方法框架，以此指導散熱器設計，本文工作也有助于開發基于有限體積法的拓撲優化計算工具。因此，本文首先對換熱器的自然對流換熱機理進行研究，建立其熱流耦合模型，并推導數值解法，以獲得對自然對流換熱問題的數值分析能力。在此基礎上，提出能強化自然對流的拓撲優化方法，并構建完整的拓撲優化流程，從而獲得對自然對流的優化能力，最后通過3組算例對所提出的方法進行了驗證。

1"流體固體統一描述的自然對流模型及求解

盡管自然對流現象廣泛存在，但其作用機理卻較為復雜：熱源的存在會影響整個空間的溫度分布情況，并造成空間中流體的密度差異，在重力場的作用下浮力產生并驅使流體流動，而流體流動又會影響溫度的分布。這是一種傳熱與流動的雙向強耦合問題。此外，當形狀待定的固體散熱器被引入時，整個對流換熱過程將更加復雜。

假設流體為不可壓縮層流，流體所受浮力采用Boussinesq假設，建立流體的熱流耦合模型。為了將固體納入熱流耦合分析過程，基于多孔介質假設對流體模型改造，建立了固體流體統一描述的熱流耦合模型。在數值分析方法上，選用有限體積法在二維情況下離散模型，推導自然對流物理場的迭代求解格式。

1.1"流體固體統一描述的自然對流換熱模型

對于流體的自然對流換熱問題，其控制方程包含質量守恒方程、動量守恒方程和能量守恒方程。假設冷卻介質流動狀態為層流，質量守恒方程和能量守恒方程中的流體特性不隨溫度和壓力變化，動量守恒方程中除體積力外其余各項中的流體特性不隨溫度和壓力變化。進一步假設流體不可壓縮，且換熱過程十分充分，已經達到穩態。在穩態狀態下，流體的特性及問題物理場（速度場、溫度場、壓力場等）分布均不隨時間變化。

綜上可知，在流體性質不變、不可壓縮、等溫和穩定流動以及忽略黏性耗散的假設下，同時為了在流體不可壓縮假設下考慮浮力的影響，引入Boussinesq假設。此時，質量方程為

[XZ（180#]Δ[XZ）]·u=0（1）

考慮浮力效應的動量守恒方程被改寫為

能量守恒方程為

式中：u、T分別為流體速度矢量、溫度；P為包含重力壓頭的流體壓力場；ρ為流體密度；cp為流體比定壓熱容；k是流體導熱系數；ρ0為參考流體密度；T0為參考流體溫度；β為熱體積膨脹系數；si是動量源項；sT是體積熱源項。

在上述控制方程中，僅考慮了流體的流動換熱模型，而未將流動與換熱特性不同的“固體”納入考慮。為能將任意形狀的固體納入熱流耦合分析過程，引入多孔介質的概念，整個問題領域被視為一個多孔介質區域，流體在其中流動時不考慮固液邊界效應，流速分布均勻且受一定的摩擦阻力。摩擦阻力會阻礙流體流動并吸收部分動量，這種阻礙用區域的滲透率來衡量。如果多孔介質的滲透率很低，那么這里的流體由于吸收了大部分的動量而幾乎是靜止的，從流動特性的角度來看，流體可以看作是固相，反之視為液相。根據多孔介質中固相與液相比例的變化來表示此處流固分布的變化。

為此，在動量守恒方程（式（2））中，引入動量吸收項中的si，基于Brinkman 摩擦項來量化多孔介質對流體流動阻礙，即

si=－μκu（4）

式中：μ為流體材料的動態黏度；κ為多孔介質域內的有效滲透率，表達式為

κ（x）=[HL（2：1，Z;2，Z]0，x∈Ωf

∞，x∈Ωs（5）

式中：Ωf為流體區域；Ωs為固體區域。同時，在能量守恒方程中，根據流體與固體材料屬性的不同重新定義熱導率k。

1.2"基于有限體積法的求解方法

考慮到自然對流具有雙向耦合性且計算量大，為提高自然對流熱沉拓撲優化效率，本文選擇計算復雜度較低且在對流換熱求解方面較為成熟的有限體積法作為數值方法。為了提高自然對流問題的數值求解效率，將自然對流控制方程重構為偽可壓縮流體格式"［18］，并引入參數向量和矩陣

式中：U是狀態變量的向量；Fc是對流通量的矩陣；Fv是黏性通量的矩陣；Q是源項；δ是偽壓縮系數"［19］；u*=［u1， u2， u3］T是速度矢量；ei （i=1，2，3）是對應方向的單位方向向量。

采用外節點有限體積單元對自然對流區域進行離散，將連接節點所屬各子區域的中心和邊界中點作界面線構建對偶剖分。二維情況下，其離散方式分別如圖1所示。圖中，Ωi是節點i的控制體積；節點j是節點i的任意相鄰節點；n"ij是單元外部法向量，位于節點i和節點j的控制體積的共同邊界S"ij處。

應用上述參數向量和矩陣，可將自然對流控制方程寫為統一形式。將自然對流控制方程對控制體積Ωi作積分，并使用高斯散度定理可得

式中：N（i）為節點i的相鄰節點集；[AKF～]c"ij和[AKF～]v"ij分別為對流通量Fc"［20］和黏性通量Fv在邊界S"ij的中點處沿其外法線方向n"ij的投影分量；ΔS"ij為控制體積Ωi和Ωj共同界面S"ij的面積（二維模型為長度）；[JB（|]Ωi[JB）|]為控制體積（二維模型為面積）；Ri（U）為式中穩態項的殘差。采用隱式歐拉方法對控制方程進行時間離散，對于穩態項在第（n+1）時間層上的殘差，對時間tn進行一階泰勒展開，并考慮守恒向量的時間相關性，可得守恒向量的迭代求解格式

2"基于有限體積法的自然對流拓撲優化

根據上節流體固體統一描述的自然對流控制方程，以及有限體積法對自然對流進行的描述及求解，本節提出了自然對流拓撲優化方法。熱量的傳遞與物質的傳輸是自然對流的主要特征，為了強化散熱過程或強化流動過程，需要對散熱器進行合理的材料分布與結構設計。

本節在基于FVM 的控制方程離散方法中，引入了基于結構的離散化描述方法，搭建了結構的描述框架并構建了流體與固體間的材料屬性插值模型；提出了拓撲優化數學模型，明確了設計變量、目標函數以及約束條件等，并提出了自然對流拓撲優化方法。

2.1"自然對流拓撲優化方法

在上節所提多孔介質域的概念下，引入了導熱率和滲透率兩種材料屬性，以此來區分固體與液體。這在物理模型上很好地解決了流體、固體的差異，但在實際求解中，計算的最小單元是以有限體網格做區分的。此時，材料屬性是基于網格分配的。

為此，對區域x，引入設計變量場γ（x）來描述整個區域內的二元（固體、液體）材料分布情況。對區域內的任意一點xi，用γ（xi）=0表示固體材料，用γ（xi）=1表示流體材料。對于0lt;γ（xi）lt;1之間非固體也非流體的部分，參考前面多孔介質的概念進一步的定義。在這一部分中，xi點處的材料被認為是多孔材料，其性質介于固體和流體之間，γ（xi）代表流體體積占比。這樣，材料的導熱系數、滲透率均可利用設計變量進行插值。

由于設計變量的引入，尋找最優材料分布的結構設計問題，變成了尋找合理的設計變量場的數值優化問題。約束函數使用體積約束，以此建立了拓撲優化的數學模型

式中：f"obj是自然對流散熱器拓撲優化的目標函數；γ是設計變量場向量；物理場u、P、T是計算域中速度矢量、壓力及溫度場的數值分布；γn是節點n對應的設計變量；Vn是節點n對應的控制體積（二維模型為面積）；V"Ωd是設計區域Ωd的體積（二維模型為面積）；Nd是設計變量的總數；r是固體體積上限。

強化散熱通過對固體結構進行優化來借助自然對流，盡可能地降低指定區域的溫度，以達到散熱效果。借鑒文獻［12-14］的思路，采用指定區域的溫度作為目標函數，此時目標函數f"obj可寫為

式中：nk為指定區域的離散單元數；Ti是其中第i個離散單元節點上的溫度；nh為計算域內的總單元數。

強化物質傳輸通過對固體結構進行優化來借助自然對流，盡可能地增強指定區域的流體流動速率，以達到增強物質傳輸效果，此時目標函數f"obj寫為

式中：v"xi是其中第i個離散單元節點上指定方向（如x方向）的速度，該速度有正負之分以強調強化流動的方向。對于強化物質傳輸，目標函數也可以選擇能耗或者關心區域的壓降作為目標函數。

自動微分方法"［21］已經在很多領域如深度學習中得到應用，該方法的基本思想在于對任意復雜函數的處理：將其視作一系列基本可為微分算子的組合，并借助鏈式法則來計算函數導數。使用該方法求導時，相關基本算子僅被進行符號轉換，因而可以保證求導的精度，同時也可以避免在解析法求導中的可能存在的表達式膨脹問題。在具體的拓撲優化方法框架構建中，本文借助自動微分工具CoDiPack"［21］來計算目標函數、約束函數對設計變量的靈敏度。

為了確保拓撲優化問題解的存在并避免形成棋盤圖案，必須對設計施加一些限制。常見的方法是對靈敏度和設計變量應用濾波器"［22］。本文同時采用了敏度過濾"［23］和密度過濾方法。

2.2"自然對流拓撲優化方法流程

結合第一節所提基于FVM的自然對流求解方法，以及本節所提的結構描述模型及整套的拓撲優化方法，針對一個具體的強化自然對流的結構設計問題，建立了自然對流拓撲優化流程，主要步驟如下。

（1）初始化。確定待優化的模型，針對當前設計模型，確定優化過程中用于初始化和判定的相關參數：設計變量場γ=［γ1，γ2，…，γn］的初始值，允許使用的最大固體體積占比r，優化循環迭代次數，目標函數收斂條件等。

（2）設立自然對流求解模型。針對上述設計模型，確定問題的尺寸、材料屬性、邊界條件等，并將計算域進行有限體網格劃分。

（3）更新設計模型。更新設計變量場，并進行濾波操作以避免棋盤格現象，同時更新材料屬性。

（4）數值分析與評估。將更新后的材料屬性嵌入自然對流控制方程中，進行迭代求解。根據設定的迭代求解次數及收斂條件，得到當前設計變量場下的物理場分布情況。根據設計需求，基于式（10）或（11）計算當前目標函數f"obj以及約束函數。基于自動微分法，分析目標函數、約束函數對設計變量的靈敏度。

（5）優化更新。使用基于梯度"［24］的移動漸近線方法（method of moving asymptotes，MMA）"［25］來對設計變量進行優化。其中，輸入為設計變量場、目標函數值與約束函數值、相關的靈敏度向量等。輸出為新的設計變量場，同時對設計變量場再次密度濾波。更新后，對結果進行判斷，包括是否達到迭代次數上限以及目標函數變換量已經小于特定值。若不滿足條件，則回到步驟（3）中，重新更新設計模型，并重復上述過程；若滿足條件，則停止迭代，并輸出當前設計變量場結果作為最終設計結果。

在優化迭代過程當中，多次涉及到對相關數據（如設計變量場、目標函數靈敏度、約束函數靈敏度）的濾波操作，且整個過程數據流動較為復雜。圖2展示了在具體實現中，相關數據的處理過程。設計變量場數據讀入自然對流求解器后，進行分析，分別得到目標函數、約束函數及其靈敏度，在此處對靈敏度進行一次濾波。將濾波后的靈敏度、目標函數、約束函數、上一代優化后的設計變量場，一起輸入到MMA優化器，得到更新后的設計變量場，并存儲。然后，此處對設計變量場進行一次濾波。根據收斂進行判定，重復進行上述過程，或者輸出結果。

3"自然對流拓撲優化求解方法驗證

首先通過二維浮力流問題，驗證第一節所提出的熱流耦合模型及數值分析方法的正確性，并選用商業軟件COMSOL進行了對比分析。隨后，以強化物質傳輸為優化目標，實施了自然對流拓撲優化設計，驗證了提出的拓撲優化方法的有效性。

3.1"自然對流模型及求解方法的驗證

二維浮力流問題為典型的自然對流問題。如圖3所示算例，封閉方形區域內充斥著流體，其上下壁為絕熱，左壁為高溫壁，壁溫Th設置為20℃，右壁為低溫壁，壁溫Tc設置為10℃，矩形腔的所有邊界均滿足無滑移邊界條件。考慮重力影響，重力加速度g方向向下，數值設置為9.8m/s2。方腔長度尺寸D設置為0.1m。流體在左側壁受熱后，膨脹并浮升，并在右側放熱后下沉，繼而形成自然對流循環。此處，流體選用空氣，材料參數如下：流體密度為"1.2886kg/m3，流體動力黏度為1.716×105Pa·s，流體比熱容為"1004.7J/（kg·K），流體體積熱膨脹系數為0.00347K"－1，流體導熱系數為0.0257W/（m·K）。

COMSOL軟件具備出色的多物理場耦合分析能力，適合用于進行自然對流問題的仿真分析。通過與商用COMSOL軟件的二維求解結果對比，可以驗證本文所提方法在二維自然對流問題分析上的合理性及求解精度。在COMSOL分析模型中，選擇二維問題，設置物理場為層流，問題狀態設置為穩態，層流節點為不可壓縮流動，并選擇包含重力選項。在多物理場節點中，選擇布辛涅斯克近似。

圖4（a）為采用COMSOL高精度網格求解結果。在網格劃分中，將整個區域使用普通的三角形網格劃分，并對邊界層加密，設置網格密度為超細化，包含16 992個域單元和600個邊界元。圖4（b）為采用COMSOL普通網格求解的結果，其與圖4（a）主要區別在于采用正四邊形網格進行劃分，并刪去對邊界層的處理。圖4（c）為采用本文所提FVM方法進行求解的結果，使用正四邊形有限體網格對問題域進行劃分，高度和寬度方向的單元數為150。

從圖4（a）的結果來看：壓力出現明顯分層現象，x方向最大速度出現在左上角與右下角，其值為0.027m·s"－1，y方向最大速度出現在兩側壁板中間，其值為0.048m·s"－1；溫度分布也出現明顯的分層現象，上層為高溫空氣，下層為低溫空氣，中層溫度位于二者之間。這與實際情況相符。從圖4（b）求解結果來看：壓力分層與圖4（a）類似，x方向最大速度出現在左上角與右下角，其值為0.028m·s"－1，y方向最大速度出現在兩側壁板中間，其值為0.048m·s"－1；溫度分層與前者類似。從圖4（c）求解結果來看，x方向最大速度出現在左上角與右下角，其值為"0.025m·s"－1，y方向最大速度出現在兩側壁板中間，其值為0.048m·s"－1；溫度分層與前者類似，壓力場的分層相同。

對比圖4中3種求解方法可知，（a）中基于邊界加密化的三角形網格分析結果更接近于流體自然對流換熱的真實值，而（b）、（c）中使用了更為簡單的網格。比較圖4中3種方法得到的結果可知，壓力場、速度場與溫度場的分布趨勢十分接近，各類誤差均小于3%，由此證明本文方法在求解流體自然對流問題時能提供較好的計算精度。

3.2"強化散熱效果的自然對流拓撲優化（散熱器）

為驗證本文方法在強化散熱拓撲優化方面的有效性，以二維流固耦合自然對流的強化散熱拓撲優化為算例進行驗證。如圖5所示，在在封閉的矩形方腔內，存在著流體區域和設計區域，設計區域內的材料屬性未定。底部有一熱源進行加熱，固體將熱量傳導至周圍的空氣，使之產生自然對流，設計目標為最小化熱源處的平均溫度。

矩形方腔尺寸[HJ2.07mm]為H×W=8cm×14cm，設計區域尺寸為h×w=5cm×8cm。矩形腔上邊界及左右邊界溫度恒定為Tw=15℃，除熱源以外的其余下邊界絕熱，而且矩形腔的所有邊界均滿足無滑移邊界條件。熱源設置為線熱源，其尺寸B=0.4cm，熱流密度q0=1.0×105W/m2。材料參數同3.1小節。此外，固體導熱系數為2.57W/（ m·K），流體反滲透率為0，固體反滲透率為1.0×105Pa·s/m2。該算例的優化函數、設計變量可描述為

在計算中，固體導熱系數為2.57W/（m·K），流體反滲透率為0，固體反滲透率為1.0×105Pa·s/m2。整個計算區域被離散為140×160個有限體胞體，其尺寸為H/160。相應地，總自由度數為90804，設計域由79×99個有限體單元組成，且優化設計變量數為8000。濾波半徑被設為有限體胞體尺寸的2.5倍，即R"min=H/64。自然對流熱傳遞的拓撲優化數學模型如式（9）所示，目標函數如式（10）所示。

熱沉固體材料體積占比小于50%，初始的材料分布采用均勻分布，即設計變量γn的初始值為0.5。設計過程包括優化的大循環和每一步優化中的仿真循環。在仿真循環中，當物理場的殘差小于10"－12量級或者相對于初始值下降了８個數量級時，認為仿真達到收斂條件。仿真循環的最大迭代次數設定為3500。對于優化循環，收斂條件設定為2次毗連迭代中，設計變量的最大變化值小于0.01。

設計迭代過程如圖6所示，平均每次迭代用時為725.5s。可以看出，最終設計結果類似于樹狀，均勻向四周展開分叉。初始布局下目標函數為116.85℃，優化后目標函數值為87.85℃，溫度降低了約30℃。經過98次迭代，最終收斂。

圖7（a）和（b）展示了拓撲優化前后的速度場與溫度場結果對比。可以看出，優化結果呈樹狀分布，有效地增大了固體與周圍流體的接觸面積，這使得熱量可以更有效地以熱傳導形式傳出，進而對熱源處進行散熱。流體流動速度有所降低，其原因是當熱流密度載荷為1.0×105W/m2時，優化結果的格拉曉夫數為5×103，此時結構對流換熱強度較低，熱擴散為主要的傳熱形式。

為對上述拓撲優化結果進行進一步驗證，采用COMSOL軟件對該算例進行了計算，采用的參數一致，拓撲優化計算結果如圖7（c）所示。可以看出，采用COMSOL的計算結果與本文方法的結果較為接近，但是略有不同，特別是其優化后的最大溫度為90.53℃，略高于本文方法得到的目標函數值87.85℃。上述對比證明了本文方法的有效性。

3.3"強化物質傳輸的自然對流拓撲優化（微泵）

針對強化物質傳輸的拓撲優化問題，目標函數為式（11），采用微泵算例進行驗證。如圖8所示，有一矩形腔體，其由2個方塊組成。上半部分為一倒U型的流道，流道內充斥著空氣并與固體區域換熱，下半部分灰色為設計域。

設計域左側為高溫壁，溫度為25℃，右側為低溫壁，溫度為15℃，其余壁面為絕熱壁，所有壁面均滿足無滑移條件。優化目標為增強黃線處流體向右的流速。腔體整體尺寸為H×L=8mm×4mm；設計域尺寸為L×L=4mm×4mm；U型流道寬度"W=0.8mm。材料參數同3.1小節。此外，固體導熱系數為2.57W/（ m·K），流體反滲透率為0，固體反滲透率為1.0×105Pa·s/m2。該算例的優化函數、設計變量可描述為

在進行實際的設計前，首先基于FVM對問題進行離散。整個計算域適用50×100個有限體單元離散，其自由度為20000；設計域由50×50個有限體單元組成，設計變量數為2 500。濾波半徑被設為有限體單元尺寸的2.5倍，即H/40。在該問題中，自然對流熱傳遞的拓撲優化數學模型如式（9）所示，設計目標為最大化目標函數計算面上的速度，且指定速度方向向右為正，目標函數采用式（11）。

設計域中固體材料體積占比上限為50%，初始的材料分布采用均勻分布，即設計變量γn的初始值為0.5。在仿真循環中，當物理場的殘差小于10"－12 數量級或者相對于初始值下降了8個數量級時，認為仿真達到收斂條件。仿真循環的最大迭代次數設定為3500。對于優化循環，收斂條件被設定為兩次毗連迭代中，設計變量的最大變化值小于0.01。為便于后續計算，設計結果被映射為一黑白的二元結果。

基于本文方法對該問題進行尋優，優化迭代過程如圖9所示。在40次迭代后，問題收斂并結束優化。問題初始目標函數值為0.012m/s，最終目標函數值為0.08m/s，且整個區域變得聯通。可以看出，最終結構為一正U型流道，與上部相通，形成了一個完整的對流通道。不足的是，最終設計結果中仍存在較多的中間密度單元，因此對最終優化結構進行后處理。首先，為滿足材料用比的約束前提，將迭代結果的密度矩陣進行過濾，得到滿足50%固體材料體積占比條件的結構。其次，以減少流道突變、增強其連通性為設計需求，對流道結構進行光滑化處理，增強其結構可描述性，并且這樣更符合實際加工技術要求性。

初始設計布局、最終設計結果，及其各自物理場的分布如圖10所示。可以看出：在最終的優化結果中，設計區域出現了近似U型的流道；在左側高溫壁的驅動下，流體膨脹，并向上浮動；在右側低溫壁的作用下，

流體收縮而下沉；U型流道的存在使得整個區域聯通，并在整個區域形成促使流體循環流動，流道內最大流體速度由0.022m/s提升至"0.049m/s。該系統無需外界動力部件，就能自發完成物質輸送，故稱之為自然對流微泵。該算例的實施，證明了本文方法對于提升強化自然對流物質輸送能力的效果。

4"結"論

本文建立了基于有限體積法的自然對流換熱器拓撲優化方法，并對方法進行了驗證。

（1）建立了自然對流換熱的物理和數學模型，提出了基于有限體積法的高效自然對流求解方法。通過引入常物性假設、Boussinesq假設等，建立了層流自然對流換熱模型。為了同時考慮流體與固體介質，基于多孔介質概念，對流體與固體引入了滲透率和導熱率等參數，據此修正了控制方程，得到了能統一描述流體與固體特性的控制方程。基于有限體積法對該方程進行求解，并推導了相應的迭代求解公式。

（2）建立了基于有限體積法的自然對流拓撲優化方法。通過構建流體與固體間的材料屬性插值模型，將結構優化設計問題轉換為設計變量場尋優的數學問題；針對自然對流拓撲優化中如何增強散熱能力與增強物質傳輸能力的問題，提出了拓撲優化的數學模型，明確了設計變量、目標函數以及約束條件等，提出了高效的自然對流拓撲優化方法框架。

（3）為驗證提出的自然對流求解方法，分別用提出的求解方法和商業軟件，對二維浮力流問題進行分析求解，驗證了所提自然對流換熱模型及求解方法在二維問題分析上的合理性。以散熱器和微泵為例進行自然對流拓撲優化設計，驗證了所提出的拓撲優化方法在優化自然對流在強化散熱和強化物質傳輸方面的可行性和有效性。

參考文獻：

［1］LEE Y J， KIM S J. Thermal optimization of the pin-fin heat sink with variable fin density cooled by natural convection ［J］. Applied Thermal Engineering， 2021， 190： 116692.

［2］田沿杰， 林曉玲， 汪林， 等. 國外艦艇一體化反應堆技術應用 ［J］. 艦船科學技術， 2019， 41（2）： 154-157.

TIAN Yanjie， LIN Xiaoling， WANG Lin， et al. Research on the application of foreign warships’ integrated reactor technology ［J］. Ship Science and Technology， 2019， 41（2）： 154-157.

［3］JOO Y， KIM S J. Comparison of thermal performance between plate-fin and pin-fin heat sinks in natural convection ［J］. International Journal of Heat and Mass Transfer， 2015， 83： 345-356.

［4］陽祥， 陶文銓. 高瑞利數下封閉腔內自然對流的數值模擬 ［J］. 西安交通大學學報， 2014， 48（5）： 27-31.

YANG Xiang， TAO Wenquan. Numerical simulations for natural convection with high Rayleigh number in a tall rectangular cavity ［J］. Journal of Xi’an Jiaotong University， 2014， 48（5）： 27-31.

［5］丁昊， 陶文銓. 水平熱板上液態金屬自然對流的數值模擬 ［J］. 工程熱物理學報， 2021， 42（8）： 2035-2039.

DING Hao， TAO Wenquan. Numerical simulation of natural convection of liquid metal above a heated horizontal plate ［J］. Journal of Engineering Thermophysics， 2021， 42（8）： 2035-2039.

［6］田東東， 王會， 刁永發， 等. 金屬泡沫銅/石蠟復合相變材料融化傳熱特性的實驗研究 ［J］. 西安交通大學學報， 2020， 54（3）： 188-196.

TIAN Dongdong， WANG Hui， DIAO Yongfa， et al. Experimental study on the melting heat transfer characteristics of paraffin/copper foam composite ［J］. Journal of Xi’an Jiaotong University， 2020， 54（3）： 188-196.

［7］解巖， 歐陽潔， 周文， 等. 自然對流換熱的高精度非結構網格有限體積法 ［J］. 計算物理， 2013， 30（3）： 337-345.

XIE Yan， OUYANG Jie， ZHOU Wen， et al. A high accuracy unstructured grid finite volume method for natural convection heat transfer ［J］. Chinese Journal of Computational Physics， 2013， 30（3）： 337-345.

［8］ALEXANDERSEN J. Topology optimisation for coupled convection problems ［D］. Copenhagen： Technical University of Denmark， 2013.

［9］BRUNS T E. Topology optimization of convection-dominated， steady-state heat transfer problems ［J］. International Journal of Heat and Mass Transfer， 2007， 50（15/16）： 2859-2873.

［10］"ALEXANDERSEN J， AAGE N， ANDREASEN C S， et al. Topology optimisation for natural convection problems ［J］. International Journal for Numerical Methods in Fluids， 2014， 76（10）： 699-721.

［11］"ALEXANDERSEN J， SIGMUND O， AAGE N. Large scale three-dimensional topology optimisation of heat sinks cooled by natural convection ［J］. International Journal of Heat and Mass Transfer， 2016， 100： 876-891.

［12］"ALEXANDERSEN J， SIGMUND O， MEYER K E， et al. Design of passive coolers for light-emitting diode lamps using topology optimisation ［J］. International Journal of Heat and Mass Transfer， 2018， 122： 138-149.

［13］"LI Hanling， LAN Daiyan， ZHANG Xianming， et al. Investigation of the parameter-dependence of topology-optimized heat sinks in natural convection ［J］. Heat Transfer Engineering， 2022， 43（11）： 922-936.

［14］"李含靈， 藍代彥， 張顯明， 等. 自然對流散熱齒的拓撲優化 ［J］. 工程熱物理學報， 2022， 43（5）： 1357-1361.

LI Hanling， LAN Daiyan， ZHANG Xianming， et al. Topology optimization of the natural convection fin heat sink ［J］. Journal of Engineering Thermophysics， 2022， 43（5）： 1357-1361.

［15］"COFFIN P， MAUTE K. A level-set method for steady-state and transient natural convection problems ［J］. Structural and Multidisciplinary Optimization， 2016， 53（5）： 1047-1067.

［16］"吳璇， 陳群. 基于不動點迭代的自然對流熱沉拓撲優化 ［J］. 工程熱物理學報， 2020， 41（9）： 2241-2246.

WU Xuan， CHEN Qun. Topology optimization of natural convection heat sink based on fixed point method ［J］. Journal of Engineering Thermophysics， 2020， 41（9）： 2241-2246.

［17］"BARI A， ZARCO-PERNIA E， GARCA DE MARA J M. A review on natural convection in enclosures for engineering applications： the particular case of the parallelogrammic diode cavity ［J］. Applied Thermal Engineering， 2014， 63（1）： 304-322.

［18］"ROE P L. Approximate Riemann solvers， parameter vectors， and difference schemes ［J］. Journal of Computational Physics， 1981， 43（2）： 357-372.

［19］"CHORIN A J. A numerical method for solving incompressible viscous flow problems ［J］. Journal of Computational Physics， 1967， 2（1）： 12-26.

［20］"WEISS J M， MARUSZEWSKI J P， SMITH W A. Implicit solution of the Navier-Stokes equations on unstructured meshes ［C］//13th Computational Fluid "Dynamics Conference. Reston， VA， USA： AIAA， 1997： AIAA 1997-2103.

［21］"SAGEBAUM M， ALBRING T， GAUGER N R. High-performance derivative computations using CoDiPack ［J］. ACM Transactions on Mathematical Software， 2019， 45（4）： 38.

［22］"SIGMUND O. Morphology-based black and white filters for topology optimization ［J］. Structural and Multidisciplinary Optimization， 2007， 33（4）： 401-424.

［23］"龍凱， 傅曉錦. 考慮密度梯度的敏度過濾方法 ［J］. 計算機輔助設計與圖形學學報， 2014， 26（4）： 669-674.

LONG Kai， FU Xiaojin. Sensitivity filtering method considering density gradient ［J］. Journal of Computer-Aided Design amp; Computer Graphics， 2014， 26（4）： 669-674.

［24］"GRIEWANK A， WALTHER A. Evaluating derivatives： principles and techniques of algorithmic differentiation ［M］. 2nd ed. Philadelphia， PA， United States： Society for Industrial and Applied Mathematics， 2008.

［25］"SVANBERG K. The method of moving asymptotes： a new method for structural optimization ［J］. International Journal for Numerical Methods in Engineering， 1987， 24（2）： 359-373.

（編輯"陶晴）