采用有限體積法的自然對流換熱拓?fù)鋬?yōu)化數(shù)值方法

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

關(guān)鍵詞："自然對流；拓?fù)鋬?yōu)化；有限體積法；強(qiáng)化散熱；流體固體統(tǒng)一描述

中圖分類號："TH122"文獻(xiàn)標(biāo)志碼：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

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

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

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

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

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

在傳統(tǒng)的力學(xué)結(jié)構(gòu)拓?fù)鋬?yōu)化領(lǐng)域，有限元方法是被廣泛使用的數(shù)值求解方法。當(dāng)問題拓展至自然對流拓?fù)鋬?yōu)化時，這種慣性也得以保留。在上述關(guān)于自然對流散熱器的拓?fù)鋬?yōu)化設(shè)計中，幾乎都采用了有限單元法作為數(shù)值求解方法。另外一方面，在流體數(shù)值分析領(lǐng)域，特別是自然對流問題，作為一種熱流雙向強(qiáng)耦合的復(fù)雜問題，其求解難度較高。有限體積法因其自身優(yōu)勢而得到了廣泛應(yīng)用，Bari等"［17］根據(jù)大量調(diào)研指出，目前有限體積法（FVM）是自然對流問題分析中最常用的數(shù)值方法。然而，目前基于FVM進(jìn)行自然對流拓?fù)鋬?yōu)化的工作還較少。對于基于FVM的熱-流耦合拓?fù)鋬?yōu)化問題，仍存在諸多問題，特別是缺少基本的優(yōu)化框架。

針對上述問題，本文旨在建立基于FVM的自然對流換熱器的拓?fù)鋬?yōu)化設(shè)計方法框架，以此指導(dǎo)散熱器設(shè)計，本文工作也有助于開發(fā)基于有限體積法的拓?fù)鋬?yōu)化計算工具。因此，本文首先對換熱器的自然對流換熱機(jī)理進(jìn)行研究，建立其熱流耦合模型，并推導(dǎo)數(shù)值解法，以獲得對自然對流換熱問題的數(shù)值分析能力。在此基礎(chǔ)上，提出能強(qiáng)化自然對流的拓?fù)鋬?yōu)化方法，并構(gòu)建完整的拓?fù)鋬?yōu)化流程，從而獲得對自然對流的優(yōu)化能力，最后通過3組算例對所提出的方法進(jìn)行了驗證。

1"流體固體統(tǒng)一描述的自然對流模型及求解

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

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

1.1"流體固體統(tǒng)一描述的自然對流換熱模型

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

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

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

考慮浮力效應(yīng)的動量守恒方程被改寫為

能量守恒方程為

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

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

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

si=－μκu（4）

式中：μ為流體材料的動態(tài)黏度；κ為多孔介質(zhì)域內(nèi)的有效滲透率，表達(dá)式為

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

∞，x∈Ωs（5）

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

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

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

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

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

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

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

2"基于有限體積法的自然對流拓?fù)鋬?yōu)化

根據(jù)上節(jié)流體固體統(tǒng)一描述的自然對流控制方程，以及有限體積法對自然對流進(jìn)行的描述及求解，本節(jié)提出了自然對流拓?fù)鋬?yōu)化方法。熱量的傳遞與物質(zhì)的傳輸是自然對流的主要特征，為了強(qiáng)化散熱過程或強(qiáng)化流動過程，需要對散熱器進(jìn)行合理的材料分布與結(jié)構(gòu)設(shè)計。

本節(jié)在基于FVM 的控制方程離散方法中，引入了基于結(jié)構(gòu)的離散化描述方法，搭建了結(jié)構(gòu)的描述框架并構(gòu)建了流體與固體間的材料屬性插值模型；提出了拓?fù)鋬?yōu)化數(shù)學(xué)模型，明確了設(shè)計變量、目標(biāo)函數(shù)以及約束條件等，并提出了自然對流拓?fù)鋬?yōu)化方法。

2.1"自然對流拓?fù)鋬?yōu)化方法

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

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

由于設(shè)計變量的引入，尋找最優(yōu)材料分布的結(jié)構(gòu)設(shè)計問題，變成了尋找合理的設(shè)計變量場的數(shù)值優(yōu)化問題。約束函數(shù)使用體積約束，以此建立了拓?fù)鋬?yōu)化的數(shù)學(xué)模型

式中：f"obj是自然對流散熱器拓?fù)鋬?yōu)化的目標(biāo)函數(shù)；γ是設(shè)計變量場向量；物理場u、P、T是計算域中速度矢量、壓力及溫度場的數(shù)值分布；γn是節(jié)點n對應(yīng)的設(shè)計變量；Vn是節(jié)點n對應(yīng)的控制體積（二維模型為面積）；V"Ωd是設(shè)計區(qū)域Ωd的體積（二維模型為面積）；Nd是設(shè)計變量的總數(shù)；r是固體體積上限。

強(qiáng)化散熱通過對固體結(jié)構(gòu)進(jìn)行優(yōu)化來借助自然對流，盡可能地降低指定區(qū)域的溫度，以達(dá)到散熱效果。借鑒文獻(xiàn)［12-14］的思路，采用指定區(qū)域的溫度作為目標(biāo)函數(shù)，此時目標(biāo)函數(shù)f"obj可寫為

式中：nk為指定區(qū)域的離散單元數(shù)；Ti是其中第i個離散單元節(jié)點上的溫度；nh為計算域內(nèi)的總單元數(shù)。

強(qiáng)化物質(zhì)傳輸通過對固體結(jié)構(gòu)進(jìn)行優(yōu)化來借助自然對流，盡可能地增強(qiáng)指定區(qū)域的流體流動速率，以達(dá)到增強(qiáng)物質(zhì)傳輸效果，此時目標(biāo)函數(shù)f"obj寫為

式中：v"xi是其中第i個離散單元節(jié)點上指定方向（如x方向）的速度，該速度有正負(fù)之分以強(qiáng)調(diào)強(qiáng)化流動的方向。對于強(qiáng)化物質(zhì)傳輸，目標(biāo)函數(shù)也可以選擇能耗或者關(guān)心區(qū)域的壓降作為目標(biāo)函數(shù)。

自動微分方法"［21］已經(jīng)在很多領(lǐng)域如深度學(xué)習(xí)中得到應(yīng)用，該方法的基本思想在于對任意復(fù)雜函數(shù)的處理：將其視作一系列基本可為微分算子的組合，并借助鏈?zhǔn)椒▌t來計算函數(shù)導(dǎo)數(shù)。使用該方法求導(dǎo)時，相關(guān)基本算子僅被進(jìn)行符號轉(zhuǎn)換，因而可以保證求導(dǎo)的精度，同時也可以避免在解析法求導(dǎo)中的可能存在的表達(dá)式膨脹問題。在具體的拓?fù)鋬?yōu)化方法框架構(gòu)建中，本文借助自動微分工具CoDiPack"［21］來計算目標(biāo)函數(shù)、約束函數(shù)對設(shè)計變量的靈敏度。

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

2.2"自然對流拓?fù)鋬?yōu)化方法流程

結(jié)合第一節(jié)所提基于FVM的自然對流求解方法，以及本節(jié)所提的結(jié)構(gòu)描述模型及整套的拓?fù)鋬?yōu)化方法，針對一個具體的強(qiáng)化自然對流的結(jié)構(gòu)設(shè)計問題，建立了自然對流拓?fù)鋬?yōu)化流程，主要步驟如下。

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

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

（3）更新設(shè)計模型。更新設(shè)計變量場，并進(jìn)行濾波操作以避免棋盤格現(xiàn)象，同時更新材料屬性。

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

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

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

3"自然對流拓?fù)鋬?yōu)化求解方法驗證

首先通過二維浮力流問題，驗證第一節(jié)所提出的熱流耦合模型及數(shù)值分析方法的正確性，并選用商業(yè)軟件COMSOL進(jìn)行了對比分析。隨后，以強(qiáng)化物質(zhì)傳輸為優(yōu)化目標(biāo)，實施了自然對流拓?fù)鋬?yōu)化設(shè)計，驗證了提出的拓?fù)鋬?yōu)化方法的有效性。

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

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

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

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

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

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

3.2"強(qiáng)化散熱效果的自然對流拓?fù)鋬?yōu)化（散熱器）

為驗證本文方法在強(qiáng)化散熱拓?fù)鋬?yōu)化方面的有效性，以二維流固耦合自然對流的強(qiáng)化散熱拓?fù)鋬?yōu)化為算例進(jìn)行驗證。如圖5所示，在在封閉的矩形方腔內(nèi)，存在著流體區(qū)域和設(shè)計區(qū)域，設(shè)計區(qū)域內(nèi)的材料屬性未定。底部有一熱源進(jìn)行加熱，固體將熱量傳導(dǎo)至周圍的空氣，使之產(chǎn)生自然對流，設(shè)計目標(biāo)為最小化熱源處的平均溫度。

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

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

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

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

圖7（a）和（b）展示了拓?fù)鋬?yōu)化前后的速度場與溫度場結(jié)果對比。可以看出，優(yōu)化結(jié)果呈樹狀分布，有效地增大了固體與周圍流體的接觸面積，這使得熱量可以更有效地以熱傳導(dǎo)形式傳出，進(jìn)而對熱源處進(jìn)行散熱。流體流動速度有所降低，其原因是當(dāng)熱流密度載荷為1.0×105W/m2時，優(yōu)化結(jié)果的格拉曉夫數(shù)為5×103，此時結(jié)構(gòu)對流換熱強(qiáng)度較低，熱擴(kuò)散為主要的傳熱形式。

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

3.3"強(qiáng)化物質(zhì)傳輸?shù)淖匀粚α魍負(fù)鋬?yōu)化（微泵）

針對強(qiáng)化物質(zhì)傳輸?shù)耐負(fù)鋬?yōu)化問題，目標(biāo)函數(shù)為式（11），采用微泵算例進(jìn)行驗證。如圖8所示，有一矩形腔體，其由2個方塊組成。上半部分為一倒U型的流道，流道內(nèi)充斥著空氣并與固體區(qū)域換熱，下半部分灰色為設(shè)計域。

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

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

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

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

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

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

4"結(jié)"論

本文建立了基于有限體積法的自然對流換熱器拓?fù)鋬?yōu)化方法，并對方法進(jìn)行了驗證。

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

（2）建立了基于有限體積法的自然對流拓?fù)鋬?yōu)化方法。通過構(gòu)建流體與固體間的材料屬性插值模型，將結(jié)構(gòu)優(yōu)化設(shè)計問題轉(zhuǎn)換為設(shè)計變量場尋優(yōu)的數(shù)學(xué)問題；針對自然對流拓?fù)鋬?yōu)化中如何增強(qiáng)散熱能力與增強(qiáng)物質(zhì)傳輸能力的問題，提出了拓?fù)鋬?yōu)化的數(shù)學(xué)模型，明確了設(shè)計變量、目標(biāo)函數(shù)以及約束條件等，提出了高效的自然對流拓?fù)鋬?yōu)化方法框架。

（3）為驗證提出的自然對流求解方法，分別用提出的求解方法和商業(yè)軟件，對二維浮力流問題進(jìn)行分析求解，驗證了所提自然對流換熱模型及求解方法在二維問題分析上的合理性。以散熱器和微泵為例進(jìn)行自然對流拓?fù)鋬?yōu)化設(shè)計，驗證了所提出的拓?fù)鋬?yōu)化方法在優(yōu)化自然對流在強(qiáng)化散熱和強(qiáng)化物質(zhì)傳輸方面的可行性和有效性。

參考文獻(xiàn)：

［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］田沿杰， 林曉玲， 汪林， 等. 國外艦艇一體化反應(yīng)堆技術(shù)應(yīng)用 ［J］. 艦船科學(xué)技術(shù)， 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］陽祥， 陶文銓. 高瑞利數(shù)下封閉腔內(nèi)自然對流的數(shù)值模擬 ［J］. 西安交通大學(xué)學(xué)報， 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］丁昊， 陶文銓. 水平熱板上液態(tài)金屬自然對流的數(shù)值模擬 ［J］. 工程熱物理學(xué)報， 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］田東東， 王會， 刁永發(fā)， 等. 金屬泡沫銅/石蠟復(fù)合相變材料融化傳熱特性的實驗研究 ［J］. 西安交通大學(xué)學(xué)報， 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］解巖， 歐陽潔， 周文， 等. 自然對流換熱的高精度非結(jié)構(gòu)網(wǎng)格有限體積法 ［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］"李含靈， 藍(lán)代彥， 張顯明， 等. 自然對流散熱齒的拓?fù)鋬?yōu)化 ［J］. 工程熱物理學(xué)報， 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］"吳璇， 陳群. 基于不動點迭代的自然對流熱沉拓?fù)鋬?yōu)化 ［J］. 工程熱物理學(xué)報， 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］. 計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報， 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.

（編輯"陶晴）