摘要：針對(duì)傳統(tǒng)方塊型金屬泡沫散熱器散熱效率不高，難以滿足高功率電子器件散熱需求的問題，發(fā)展了基于格子Boltzmann和水平集方法的三維拓?fù)鋬?yōu)化方法。首先，采用格子Boltzmann方法求解多孔介質(zhì)內(nèi)的流動(dòng)和傳熱過程；然后，采用伴隨格子Boltzmann方法計(jì)算梯度；最后，根據(jù)反應(yīng)-擴(kuò)散方程更新水平集函數(shù)。通過對(duì)底部放置發(fā)熱元件方腔的自然對(duì)流過程開展拓?fù)鋬?yōu)化，得到特征格拉曉夫數(shù)為2.4×102～1.2×105的4種優(yōu)化散熱器結(jié)構(gòu)，定量評(píng)估了散熱器的強(qiáng)化換熱性能并分析相關(guān)機(jī)理。研究結(jié)果表明：隨著格拉曉夫數(shù)的增大，散熱主導(dǎo)機(jī)制由熱傳導(dǎo)轉(zhuǎn)變?yōu)闊釋?duì)流，優(yōu)化結(jié)構(gòu)由伸展的樹枝狀向收縮的花朵狀結(jié)構(gòu)轉(zhuǎn)變，以留出更多空間使得流動(dòng)渦旋充分發(fā)展；與傳統(tǒng)的金屬泡沫方塊結(jié)構(gòu)和開槽的金屬泡沫塊結(jié)構(gòu)相比，所提優(yōu)化結(jié)構(gòu)表現(xiàn)出優(yōu)異的散熱性能，其分支結(jié)構(gòu)和中空結(jié)構(gòu)能夠優(yōu)化方腔內(nèi)的流動(dòng)，從而將散熱效率提升了29.7%以上，表明了所提拓?fù)鋬?yōu)化方法的有效性。研究工作可為新型金屬泡沫散熱器的設(shè)計(jì)提供理論指導(dǎo)。

關(guān)鍵詞：散熱器設(shè)計(jì)；強(qiáng)化換熱；伴隨格子Boltzmann方法；金屬泡沫；拓?fù)鋬?yōu)化

中圖分類號(hào)：TK414.2"文獻(xiàn)標(biāo)志碼：A

DOI：10.7652/xjtuxb202410007"文章編號(hào)：0253-987X（2024）10-0072-11

Design of Metal Foam Heat Sinks Guided by the Topology Optimization of Natural Convection

LUO Jiwang， CHEN Li， ZHENG Xinjian， YANG Qirui， TAO Wenquan

（Key Laboratory of Thermo-Fluid Science and Engineering of MOE， Xi’an Jiaotong University， Xi’an 710049， China）

Abstract：To overcome the challenges of low heat dissipation efficiency in traditional block-type metal foam heat sinks and their inability to meet the cooling requirements of high-power electronic devices， a three-dimensional topology optimization method based on the lattice Boltzmann and level set methods is developed. Firstly， the lattice Boltzmann method is employed to simulate the flow and heat transfer processes within the porous medium. Then， the adjoint lattice Boltzmann method is used to calculate the gradients. Finally， the level set function is updated based on the reaction-diffusion equation. By conducting topology optimization of the natural convection process in a square cavity with heat-generating elements at the bottom， four optimized heat sink structures with characteristic Grashof numbers ranging from 2.4×102 to 1.2×105 are obtained. The enhanced heat transfer performance of the heat sinks is quantitatively evaluated， and the underlying mechanisms are analyzed. The research findings indicate that as the Grashof number increases， the dominant heat transfer mechanism shifts from conduction to convection. The optimized structures transform from extended branching patterns to compact flower-like structures， creating additional space for the development of more flow vortices. Compared with traditional block-shaped metal foam structures and slotted metal foam block structures， the proposed optimized structures demonstrate superior heat dissipation performance. The branch structures and hollow structures optimize the flow within the cavity， resulting in an improvement in heat dissipation efficiency of over 29.7%. These results highlight the effectiveness of the proposed topology optimization method and provide theoretical guidance for the design of novel metal foam heat sinks.

Keywords：heat sink design; enhanced heat transfer; adjoint lattice Boltzmann method; metal foam; topology optimization

日益增長的芯片功率密度對(duì)當(dāng)前熱管理技術(shù)提出了極大挑戰(zhàn)，使用多孔金屬泡沫強(qiáng)化換熱正在受到越來越多研究人員的關(guān)注。金屬泡沫是一類包含高導(dǎo)熱系數(shù)金屬骨架和大量連通孔隙的多孔介質(zhì)，具有低密度、高比表面積及高導(dǎo)熱系數(shù)等優(yōu)點(diǎn)，被廣泛應(yīng)用于強(qiáng)化單相強(qiáng)制對(duì)流換熱、單相自然對(duì)流換熱、池沸騰及流動(dòng)沸騰換熱等換熱過程［1］。早在21世紀(jì)初，金屬泡沫就被應(yīng)用于強(qiáng)化電子器件自然對(duì)流冷卻相關(guān)方向。Mahdi等［2］發(fā)現(xiàn)將鋁金屬泡沫作為CPU散熱器，可在傳統(tǒng)翅片散熱器基礎(chǔ)上進(jìn)一步將總熱阻降低70%。此后，研究人員深入研究了影響金屬泡沫散熱性能的主要因素，如泡沫孔密度、樣品高度及長-寬-高比、擺放傾角、開槽寬度等，并優(yōu)化了金屬泡沫結(jié)構(gòu)強(qiáng)化散熱效果［1］。然而，在上述研究中，參數(shù)化研究所獲得的優(yōu)化參數(shù)或結(jié)構(gòu)并不能使金屬泡沫的拓?fù)浣Y(jié)構(gòu)發(fā)生根本改變，因而難以獲得能使散熱性能大幅提升的散熱器結(jié)構(gòu)。

拓?fù)鋬?yōu)化是獲得高性能優(yōu)化結(jié)構(gòu)設(shè)計(jì)的一種可靠方法，其基本思想是求解給定設(shè)計(jì)域內(nèi)的材料最優(yōu)分布問題。拓?fù)鋬?yōu)化最早起源于結(jié)構(gòu)力學(xué)，經(jīng)過30余年的快速發(fā)展，已被應(yīng)用于流體力學(xué)及傳熱學(xué)［3-5］。一般而言，拓?fù)鋬?yōu)化方法根據(jù)其結(jié)構(gòu)表征可以分為密度法和水平集法兩類，前者使用與多孔介質(zhì)類似的具有中間物性值的過渡材料代替不連續(xù)界面［6］，后者則引入高維水平集函數(shù)（LSF）的0等值線來表征材料的分布［7］。2003年，Borrvall和Petersson［8］首次將拓?fù)鋬?yōu)化應(yīng)用于求解以降低耗散功率為目標(biāo)的Stokes流動(dòng)優(yōu)化問題，開創(chuàng)了流動(dòng)問題拓?fù)鋬?yōu)化研究的先河。隨后，Yoon［9］和Dede［10］幾乎同時(shí)將基于密度法的拓?fù)鋬?yōu)化拓展到流動(dòng)傳熱領(lǐng)域，用以獲得高效強(qiáng)制對(duì)流換熱器結(jié)構(gòu)。由于自然對(duì)流過程的強(qiáng)耦合和強(qiáng)非線性特點(diǎn)，關(guān)于自然對(duì)流的拓?fù)鋬?yōu)化研究起步較晚。2014年，Alexandersen等［11］首次開展了二維自然對(duì)流過程的拓?fù)鋬?yōu)化研究，獲得了不同格拉曉夫數(shù)下的散熱器結(jié)構(gòu)。隨后，他們進(jìn)一步開展了相關(guān)的三維拓?fù)鋬?yōu)化研究，獲得了可以兼顧提升散熱性能和降低材料消耗的散熱器結(jié)構(gòu)［12］。Coffin等［13］和Li等［14］采用基于水平集方法的拓?fù)鋬?yōu)化方法，獲得了二維及三維的優(yōu)化散熱器結(jié)構(gòu)。需要指出的是，上述研究獲得的優(yōu)化散熱器結(jié)構(gòu)均由固體金屬而非金屬泡沫制成。考慮到金屬泡沫在增加比表面積、減少散熱器質(zhì)量、降低金屬材料消耗成本等方面的巨大優(yōu)勢，基于金屬泡沫強(qiáng)化自然對(duì)流的拓?fù)鋬?yōu)化研究亟待開展。

在以往的拓?fù)鋬?yōu)化研究中，一般使用直接離散求解Navier-Stokes方程的宏觀方法，如有限元法（FEM）及有限體積法（FVM），進(jìn)行數(shù)值模擬及相關(guān)優(yōu)化研究，僅有少部分研究選擇基于氣體動(dòng)理論的格子Boltzmann方法（LBM）。以往的研究表明，LBM法具有計(jì)算效率高、并行性好、處理復(fù)雜界面能力強(qiáng)等優(yōu)點(diǎn)，因此被廣泛應(yīng)用于多孔介質(zhì)流及多相流等問題［15］。在拓?fù)鋬?yōu)化中，除了直接對(duì)所關(guān)心的物理問題進(jìn)行仿真外，還需求解目標(biāo)函數(shù)關(guān)于設(shè)計(jì)變量的梯度以驅(qū)動(dòng)優(yōu)化算法尋得優(yōu)化解，因此通常還需求解額外的伴隨方程以簡化梯度計(jì)算。在以LBM為基礎(chǔ)的拓?fù)鋬?yōu)化方法中，這種求解伴隨方程策略也被稱為伴隨LBM（ALBM）。2007年，Pingen等［16］率先將LBM和ALBM法應(yīng)用于流體流動(dòng)拓?fù)鋬?yōu)化，成功復(fù)現(xiàn)了文獻(xiàn)［8］的結(jié)果。接著，他們將提出的方法應(yīng)用于非牛頓流體流動(dòng)以及流動(dòng)換熱過程的拓?fù)鋬?yōu)化。然而，在他們提出的方法中，必須要組裝整個(gè)系統(tǒng)的Jacobian矩陣并求其逆矩陣，實(shí)際上與LBM法計(jì)算簡潔的特性背道而馳，由于優(yōu)化算法的求解效率大大降低，所以難以被應(yīng)用于大尺度拓?fù)鋬?yōu)化求解。2013年，Krause等［17］提出了基于連續(xù)Boltzmann方程的ALBM法，推導(dǎo)得到的伴隨方程與原始Boltzmann方程具有相似的形式，可以進(jìn)一步采用類似LBM法的離散格式求解，從而很好地保持了LBM法的諸多優(yōu)點(diǎn)。隨后，Yaji等［18］基于連續(xù)Boltzmann方程導(dǎo)出了ALBM方程并建立了相應(yīng)的伴隨邊界格式。后來，Yaji等［19］又提出了另一種基于速度離散Boltzmann方程的ALBM法，可以引入常用的邊界格式（如Zou-He格式）推導(dǎo)伴隨方程。在Yaji等工作的基礎(chǔ)上，Dugast等［20］提出了一種新的ALBM推導(dǎo)策略，并引入多松弛時(shí)間（MRT）LBM法來提高伴隨計(jì)算的穩(wěn)定性。近期，Luo等［21-22］和Tanabe等［23］分別通過發(fā)展相應(yīng)的ALBM模型，獲得了自然對(duì)流過程的二維及三維拓?fù)鋬?yōu)化結(jié)構(gòu)，再次證明了LBM法對(duì)于復(fù)雜多物理場拓?fù)鋬?yōu)化的適用性和獨(dú)特優(yōu)勢。

在Luo等［24］近期的研究中，首次提出了自然對(duì)流系統(tǒng)內(nèi)基于LBM的金屬泡沫散熱器拓?fù)鋬?yōu)化方法，并獲得了多個(gè)參數(shù)組合下的二維拓?fù)鋬?yōu)化結(jié)構(gòu)?；诖?，本文針對(duì)金屬泡沫散熱器的三維拓?fù)鋬?yōu)化，采用基于ALBM及水平集法的拓?fù)鋬?yōu)化方法，獲得了不同格拉曉夫數(shù)下金屬泡沫散熱器的最優(yōu)結(jié)構(gòu)，分析了強(qiáng)化傳熱百分比及相關(guān)強(qiáng)化傳熱機(jī)理，并與傳統(tǒng)金屬及開槽金屬泡沫方塊進(jìn)行對(duì)比。該研究工作為金屬泡沫散熱器結(jié)構(gòu)優(yōu)化設(shè)計(jì)提供了新的技術(shù)途徑，也為新一代輕質(zhì)高效多孔介質(zhì)結(jié)構(gòu)優(yōu)化設(shè)計(jì)提供了理論指導(dǎo)。

1"格子Boltzmann模型

1.1"流動(dòng)LBM模型

本文采用Guo等［25］發(fā)展的多孔介質(zhì)流動(dòng)格子Boltzmann模型，模擬孔隙及多孔介質(zhì)內(nèi)的流動(dòng)。在LBM模擬中，對(duì)碰撞項(xiàng)采用BGK碰撞模型近似后，分布函數(shù)的演化表達(dá)式可寫為

fit+ci·Δfi=－1τf（fi－feqi）+Fi（1）

式中：ci為表征離散空間方向的向量，也被稱為離散速度；fi為空間x點(diǎn)處、t時(shí)刻、ci方向上的速度分布函數(shù)，對(duì)于三維流動(dòng)，通常采用D3Q19模型所規(guī)定的19個(gè)離散速度表示

c0c1c2c3c4c5c6c7c8c9c10c11c12c13c14c15c16c17c18=

01－100001－11－11－11－100000001－10011－1－100001－11－1000001－1000011－1－111－1－1（2）

τf為與速度分布函數(shù)f及流體有效黏度νe相關(guān)的無量綱松弛時(shí)間

νe=13（τf－0.5）（3）

feqi為平衡態(tài)分布函數(shù)，對(duì)于多孔介質(zhì)流動(dòng)可表示為

feqi=ωiρ1+3（ci·u）+92ε（ci·u）2－32εu2（4）

其中，ε為多孔介質(zhì)的當(dāng)?shù)乜紫堵?；ρ、u分別為流體密度和速度矢量，可由分布函數(shù)fi的零階矩和一階矩表示為

ρ=∑8i=0fi；"ρu=∑8i=0cifi（5）

ωi為權(quán)重系數(shù)，D3Q19模型的ωi可寫為

ωi=1/3，i=0

1/18，i=1～6

1/36，i=7～18（6）

Fi為離散外力項(xiàng)，可表示如下

Fi=ωiρ

3ci·F+1－12τf·9ε（ci·u）（ci·F）－3ε（u·F）（7）

其中，F(xiàn)為外力項(xiàng)向量，對(duì)于本文中研究的多孔介質(zhì)內(nèi)的自然對(duì)流過程，可表示為Darcy阻力項(xiàng)、Forchheimer阻力項(xiàng)以及Boussinesq假設(shè)確定的浮升力項(xiàng)三者之和

F=－ενfku－εFεkuu－εβg（T－T0）（8）

其中，νf為流體黏度；k為滲透率；Fε為形狀因子；β為體積膨脹系數(shù)；g為重力加速度向量；T為溫度；T0為參考溫度。

對(duì)于本文求解的封閉腔內(nèi)自然對(duì)流過程，四周及上、下壁均為無滑移壁面，因此采用標(biāo)準(zhǔn)反彈格式

fi，（ci·nlt;0）=fi，（ci·ngt;0）（9）

1.2"傳熱LBM模型

通過求解速度分布函數(shù)f獲得密度和速度后，引入溫度分布函數(shù)g求解溫度場。溫度分布函數(shù)的演化方程［26］與式（1）相類似

git+ci·Δgi=－1τg（gi－geqi）+ω′iS（10）

式中：ci仍為離散速度，對(duì)于求解能量守恒方程這類被動(dòng)標(biāo)量方程，可采用更少的離散速度，這里采用D3Q7模型，即僅采用式（2）的前7個(gè)離散速度；S為熱源項(xiàng)；τg為與溫度分布函數(shù)g相關(guān)的無量綱松弛時(shí)間，可由有效熱擴(kuò)散系數(shù)ae=σ0（τg－0.5）/3確定，其中σ0為有效比熱與流體比熱之比的參考值；geqi為平衡態(tài)分布函數(shù)，可寫為

geqi=T（σ－σ0）+ω′iT（σ0+4ci·u），i=0

ω′iT（σ0+4ci·u），i=1～6（11）

其中，σ為有效比熱與流體比熱之比；ω′i為D3Q7模型的權(quán)重系數(shù)，取值如下

ω′i=1/4，i=0

1/8，i=1～6（12）

溫度T可通過對(duì)g求和獲得，表示為

T=∑6i=0gi（13）

至于熱邊界條件，除底面為絕熱壁面外，四周及頂壁均為溫度為0的恒溫壁面。絕熱壁面處采用與式（9）類似的反彈邊界，恒溫壁面處采用如下邊界格式（以頂壁為例）

g6=T0－g0－g1－g2－g3－g4－g5（14）

求解1.2小節(jié)中的流動(dòng)模型及1.3小節(jié)中的傳熱模型，實(shí)際上相當(dāng)于求解如下質(zhì)量守恒方程、多孔介質(zhì)流動(dòng)量守恒方程及多孔介質(zhì)流能量守恒方程［24］

Δ·u=0

ut+（u·Δ）uε=－1ρΔ（εp）+Δ·（νeΔu）+F

σTt+Δ·（Tu）=Δ·（aeΔT）+S （15）

式中：p為壓力。

將式（15）所述方程無量綱化，可得到控制多孔介質(zhì)自然對(duì)流的無量綱數(shù)，分別為達(dá)西數(shù)Da、普朗特?cái)?shù)Pr、格拉曉夫數(shù)Gr以及熱擴(kuò)散系數(shù)比γ，定義為

Da=kh2； Pr=νfaf

Gr=βg（Ts－T0）h3ν2f； γ=aeaf（16）

式中：h為特征高度；Ts為發(fā)熱元件的平均溫度；af為流體熱擴(kuò)散系數(shù)。

2"拓?fù)鋬?yōu)化模型

2.1"優(yōu)化問題定義

一般地，拓?fù)鋬?yōu)化問題可視為給定設(shè)計(jì)域內(nèi)的材料分布優(yōu)化問題。對(duì)于本文中待解的金屬泡沫散熱器拓?fù)鋬?yōu)化問題，需要優(yōu)化金屬泡沫的結(jié)構(gòu)設(shè)計(jì)，其中金屬泡沫作為單一多孔材料應(yīng)具有恒定的孔隙率和恒定的有效輸運(yùn)系數(shù)。換言之，我們需要確定設(shè)計(jì)域內(nèi)金屬泡沫和純流體兩種材料在空間中的最優(yōu)分布，因此可采用如下特征方程表述

χ（x）=0，x為金屬泡沫點(diǎn)

1，x為純流體點(diǎn)（17）

由于上述特征函數(shù)高度不連續(xù)，導(dǎo)致所定義的優(yōu)化問題是病態(tài)的，所以進(jìn)一步采用連續(xù)的水平集函數(shù)φ（x）的0等值面表征金屬泡沫與流體的界面，寫為

0lt;φ（x）≤1，x為純流體點(diǎn)

φ（x）=0，x在界面處

－1≤φ（x）lt;1，x為金屬泡沫點(diǎn)（18）

由此，特征方程與水平集函數(shù)間的映射關(guān)系可表示如下

χφ（x）=

0，－1≤φ（x）lt;0

1，0≤φ（x）≤1（19）

式中：χφ為設(shè)計(jì)變量。

根據(jù)上述材料分布的定義，整個(gè)拓?fù)鋬?yōu)化問題［24］可表述為

min χφ（J+R）=12V20∫X0T2dX+12κ∫X1Δφ2dX

s.t."G=∫X1（1－χφ）dX－Vmax∫X1dX≤0

fit+ci·Δfi=－1τf（fi－feqi）+Fi

git+ci·Δgi=－1τg（gi－geqi）+ω′iS（20）

式中：J為目標(biāo)函數(shù)，定義為發(fā)熱區(qū)域X0平均溫度的平方；V0為X0的體積；X1為設(shè)計(jì)域；R為Tikhonov正則項(xiàng)；κ為正則化系數(shù)；G為體積約束；Vmax為給定金屬泡沫體積分?jǐn)?shù)的上限。

需要指出的是，式（1）和式（10）作為優(yōu)化問題的等式約束，隱式地規(guī)定了設(shè)計(jì)變量χφ與狀態(tài)變量f、g之間的函數(shù)關(guān)系，該函數(shù)關(guān)系可通過拓展表示如下

fit+ci·Δfi=－1τffi－feqi［ρ（f），u（f），ε（χφ）］+

Fi［ρ（f），u（f），T（g），ε（χφ），k（χφ）］（21）

git+ci·Δgi=－1τg（χφ）gi－geqi［u（f），T（g）］+ω′iS（22）

其中，孔隙率ε、滲透率k和松弛時(shí)間τg與設(shè)計(jì)變量χφ間的關(guān)系可采用材料物性理性近似格式表示

1ε（χφ）=1εf+1εp－1εfq1（1－χφ）q1+χφ

1k（χφ）=1kf+1kp－1kfq1（1－χφ）q1+χφ

1τg（χφ）=1τg，f+1τg，p－1τg，fq2（1－χφ）q2+χφ（23）

其中，下標(biāo)f、p分別代表純流體點(diǎn)和多孔介質(zhì)點(diǎn)；q1和q2為不同插值項(xiàng)的懲罰系數(shù)。

2.2"反應(yīng)-擴(kuò)散方程

在優(yōu)化迭代過程中，水平集函數(shù)及其對(duì)應(yīng)的特征函數(shù)不斷更新，從初始結(jié)構(gòu)逐步演化為優(yōu)化結(jié)構(gòu)。這一演化是由最速下降法驅(qū)動(dòng)的，對(duì)式（20）中定義的目標(biāo)函數(shù)求導(dǎo)，可得到如下反應(yīng)-擴(kuò)散方程形式的結(jié)構(gòu)演化方程

φt=－K（ΔJ－κΔ2φ）（24）

式中：K為步長；ΔJ為目標(biāo)函數(shù)J關(guān)于水平集函數(shù)φ（x）的梯度。

本研究中，采用伴隨方法求解得到梯度，即通過求解額外的輔助變量（伴隨變量）簡化梯度的計(jì)算。具體來說，通過拉格朗日乘子法構(gòu)造如下拉格朗日函數(shù)

L（f*i，fi，g*i，gi，χφ，λ）=J（gi）+

R1（f*i，fi，gi，χφ）+R2（g*i，gi，fi，χφ）+λG（χφ）（25）

式中：λ為體積約束的拉格朗日乘子；R1和R2均為約束項(xiàng)，定義為伴隨變量f*i及g*i與相應(yīng)等式約束的內(nèi)積，表達(dá)式如下

R1=∫tf0∫X∑18i=0f*ifit+ci·Δfi+1τf（fi－feqi）－FidXdt（26）

R2=∫tf0∫X∑6i=0g*igit+ci·Δgi+1τg（gi－geqi）－ω′iSdXdt（27）

式中：tf為達(dá)到穩(wěn)態(tài)所需的時(shí)間；X為計(jì)算域。

根據(jù)之前的研究［18］，若伴隨變量滿足如下形式的伴隨方程

Lfδf=R1fδf+R2fδf=0（28）

Lgδg=Jgδg+R1gδg+R2gδg=0（29）

那么梯度ΔJ就可簡化為

ΔJ=dJdφ=Lχφ=R1χφ+R2χφ+λGχφ（30）

下面將對(duì)伴隨方程和梯度表達(dá)式的具體形式做詳細(xì)介紹。

2.3"伴隨模型和梯度計(jì)算

由于本文中使用的LBM模型與我們之前工作中的LBM模型是一致的，僅在離散速度的具體取值上有差別，所以相應(yīng)的伴隨方程的推導(dǎo)過程也完全一致。本文將不再詳細(xì)展開相應(yīng)的推導(dǎo)過程，僅對(duì)伴隨方程本身及其求解過程作詳細(xì)介紹，對(duì)推導(dǎo)過程感興趣的讀者可以參考文獻(xiàn)［24］。將式（28）和（29）展開，逐項(xiàng)推導(dǎo)可得到對(duì)應(yīng)的伴隨方程，分別表示如下

－f*it－ci·Δf*i+1τff*i－∑18j=0f*jfjeqfi－

∑18j=0f*jFjfi－1τg∑6k=0g*kgeqkfi=0（31）

－g*it－ci·Δg*i+1τgg*i－∑6k=0g*kgeqkgi－

∑18j=0f*jFjgi+TV20t=tf=0（32）

上述伴隨方程的邊界條件可由式（28）和（29）展開得到。例如，對(duì)于無滑移壁面的反彈格式，相應(yīng)的伴隨邊界條件為

f*i，（ci·ngt;0）=f*i，（ci·nlt;0）（33）

對(duì)于式（14）所示的恒溫壁面邊界格式，伴隨邊界格式如下

g*i，（ci·ngt;0）=－g*i，（ci·nlt;0）（34）

此外，伴隨方程的初始條件也可一并得到

f*it=tf=0 ； g*it=tf=0（35）

上述伴隨方程、邊界條件以及初始條件，一并構(gòu)成了完整的伴隨模型。在時(shí)間和空間上采用與LBM法類似的一階顯式格式離散伴隨方程式（31）及（32），可得到對(duì)應(yīng)的伴隨格子Boltzmann方程。進(jìn)一步通過依次實(shí)施伴隨碰撞、伴隨遷移和伴隨邊界條件，可求解得到伴隨分布函數(shù)f*和g*，這種方法即被稱為ALBM法。

進(jìn)一步展開式（30），得到梯度的表達(dá)式

ΔJ=－λ+∑6k=0g*k（gk－geqk）（1/τg）χφ

+∑18j=0f*j－1τffeqjεεχφ－Fjεεχφ+Fjkkχφ（36）

上式中的各個(gè)偏導(dǎo)數(shù)項(xiàng)可通過式（21）～（22）的函數(shù)關(guān)系，以及第1節(jié)中的相應(yīng)表達(dá)式計(jì)算得到其顯式表達(dá)式，無須采用迭代算法求解。由此發(fā)現(xiàn)，在使用伴隨方法后，梯度ΔJ實(shí)際上是相應(yīng)的分布函數(shù)f、g，以及伴隨分布函數(shù)f*、g*的顯式表達(dá)式。一旦通過LBM法求解得到了f、g，再通過ALBM法求解得到f*和g*，就可以快捷地計(jì)算得到梯度ΔJ。式（36）中的拉格朗日乘子λ用于準(zhǔn)確實(shí)施體積約束［27］，計(jì)算公式可寫為

λ=∫X1ΔJ－κΔ2φdX∫X1dXeG（37）

在獲得梯度后，即可根據(jù)“反應(yīng)-擴(kuò)散”方程更新設(shè)計(jì)域內(nèi)的結(jié)構(gòu)設(shè)計(jì)，而后進(jìn)入下一步優(yōu)化迭代?？傮w的優(yōu)化算法流程如圖1所示，并可概述如下：①首先給定設(shè)計(jì)域內(nèi)的初始結(jié)構(gòu)；②通過求解式（1）和（10）求解自然對(duì)流傳熱過程，得到f和g的解；③通過求解式（31）和（32）求解伴隨方程，得到f*和g*的解；④通過式（36）基于當(dāng)前f、g以及f*、g*計(jì)算梯度；⑤根據(jù)式（24）的反應(yīng)-擴(kuò)散方程更新水平集函數(shù)，基于更新后的水平集函數(shù)使用式（19）更新結(jié)構(gòu)；⑥重復(fù)步驟②～⑤直至優(yōu)化結(jié)構(gòu)及目標(biāo)函數(shù)收斂，得到最終的優(yōu)化結(jié)構(gòu)。

由于本文是文獻(xiàn)［24］二維金屬泡沫拓?fù)鋬?yōu)化工作的直接拓展，相關(guān)數(shù)值模型和優(yōu)化算法的驗(yàn)證已在其中展示，有興趣的讀者可以參考文獻(xiàn)［24］。

3"結(jié)果與討論

本文中，自然對(duì)流拓?fù)鋬?yōu)化問題的設(shè)置如圖2所示，放置尺寸為16×16×8（格子單位）的發(fā)熱元件（固體加熱域）于底部絕熱壁面的正中央，令其具有S=0.1（格子單位）的恒定均勻體熱源。加熱域正上方為112×112×80的設(shè)計(jì)域，內(nèi)可放置由金屬泡沫材料制成的散熱器，用以冷卻發(fā)熱元件。計(jì)算域整體尺寸為160×160×160，除底面外，四周及頂部均為T0=0℃的冷壁面。計(jì)算域內(nèi)填充的空氣作為換熱介質(zhì)，在加熱域與冷壁面間溫差的作用下，受熱體積膨脹，密度減小而產(chǎn)生浮力，將在計(jì)算域內(nèi)形成自然對(duì)流，加速發(fā)熱元件與冷壁面間的熱交換。上述幾何、網(wǎng)格及工況設(shè)置均參考文獻(xiàn)［22］，相關(guān)網(wǎng)格獨(dú)立性的討論也已在其中說明。需要指出的是，由于開放邊界下要求設(shè)計(jì)域外有更大的外部空間［28］，本文采用與設(shè)計(jì)域存在適當(dāng)距離的低溫壁面代替外部大空間以減少計(jì)算量。

上述計(jì)算域內(nèi)自然對(duì)流過程的控制方程由式（15）給出。由式（15）可知，控制上述自然對(duì)流過程的無量綱數(shù)分別為達(dá)西數(shù)Da、普朗特?cái)?shù)Pr、格拉曉夫數(shù)Gr和熱擴(kuò)散系數(shù)比γ。本文中，由于介質(zhì)為空氣，普朗特?cái)?shù)Pr=0.692；采用文獻(xiàn)［29］中的金屬泡沫樣品，達(dá)西數(shù)Da=1.38×10－5；熱擴(kuò)散系數(shù)比γ=114.1。格拉曉夫數(shù)Gr定義為浮升力與黏性力之比，可用于確定自然對(duì)流過程中傳熱的主導(dǎo)機(jī)制。當(dāng)Gr較小時(shí)，浮升力也較小，計(jì)算域內(nèi)流體流速較低，對(duì)流傳熱較弱，發(fā)熱元件與冷壁面間主要依靠熱傳導(dǎo)進(jìn)行換熱；當(dāng)Gr較大時(shí)，較大的浮升力使得流體快速在發(fā)熱元件與冷壁面間循環(huán)流動(dòng)，形成多個(gè)大的渦旋，強(qiáng)化了對(duì)流傳熱，熱對(duì)流取代熱傳導(dǎo)成為傳熱主導(dǎo)機(jī)制。通過改變體積膨脹系數(shù)β，對(duì)多個(gè)Gr下的自然對(duì)流過程進(jìn)行拓?fù)鋬?yōu)化，得到不同換熱主導(dǎo)機(jī)制下的優(yōu)化散熱器結(jié)構(gòu)。

值得指出的是，由于本文考慮的自然對(duì)流過程相對(duì)應(yīng)的Gr在102～105之間，因而可獲得相應(yīng)的穩(wěn)態(tài)層流解，加之問題的對(duì)稱性，實(shí)際采用的計(jì)算區(qū)域僅為圖2的1/4，且在計(jì)算域內(nèi)部的切割面上使用對(duì)稱邊界條件。為了避免設(shè)計(jì)域內(nèi)全部填充金屬泡沫的平凡解，設(shè)置金屬泡沫填充的體積分?jǐn)?shù)上限Vmax為0.1。

通過調(diào)節(jié)體積膨脹系數(shù)β，可以改變對(duì)流換熱強(qiáng)度，進(jìn)而可以在優(yōu)化完成后通過式（16）后驗(yàn)計(jì)算相應(yīng)的Gr。為呈現(xiàn)不同傳熱主導(dǎo)機(jī)制下的優(yōu)化結(jié)構(gòu)，β分別取1×10－9、1×10－8、1×10－7和1×10－6，對(duì)應(yīng)的Gr分別為2.4×102、2.3×103、1.7×104和1.2×105。為了評(píng)估本文所提優(yōu)化散熱器結(jié)構(gòu)對(duì)強(qiáng)化自然對(duì)流傳熱的效果，以2種傳統(tǒng)結(jié)構(gòu)作為參考結(jié)構(gòu)進(jìn)行對(duì)比分析。參考結(jié)構(gòu)1為常用的金屬泡沫方塊結(jié)構(gòu)，參考結(jié)構(gòu)2為開槽的金屬泡沫塊結(jié)構(gòu)。為了保證比較的公平性，方塊結(jié)構(gòu)的尺寸為56×56×32，開槽金屬泡沫塊的翅片尺寸為6×66×33，這使得2種參考結(jié)構(gòu)的金屬泡沫填充體積與優(yōu)化結(jié)構(gòu)相等。在后續(xù)進(jìn)行的參考結(jié)構(gòu)自然對(duì)流仿真中，由于特征溫差變化，導(dǎo)致實(shí)際的Gr并不完全與上述4個(gè)特征Gr相等，但這種變化不影響定性分析，因此為了分析方便，本文始終選取上述4個(gè)特征Gr進(jìn)行后續(xù)討論。

采用第1、2節(jié)中發(fā)展的拓?fù)鋬?yōu)化方法對(duì)設(shè)計(jì)域內(nèi)金屬泡沫材料的分布進(jìn)行拓?fù)鋬?yōu)化，得到不同Gr下優(yōu)化后的散熱器結(jié)構(gòu)。圖3、圖4分別給出了優(yōu)化結(jié)構(gòu)及2種參考結(jié)構(gòu)的速度和溫度分布對(duì)比。優(yōu)化結(jié)構(gòu)及參考結(jié)構(gòu)2在計(jì)算域內(nèi)沿著z方向中心線上的速度和溫度分布由圖5給出。對(duì)應(yīng)的Ts列于表1。

由圖3～圖5及表1可知，隨著Gr增加，計(jì)算域內(nèi)的流體流速顯著增加，發(fā)熱元件及散熱器本身的溫度均有明顯下降。觀察各個(gè)Gr下的優(yōu)化散熱器結(jié)構(gòu)可以發(fā)現(xiàn)，當(dāng)Gr低至2.4×102時(shí)，優(yōu)化結(jié)構(gòu)呈現(xiàn)復(fù)雜的樹枝狀，由指向計(jì)算域?qū)蔷€方向的主分支和一系列次級(jí)分支構(gòu)成，可將發(fā)熱元件產(chǎn)生的熱量傳導(dǎo)至周圍冷壁面。當(dāng)Gr從2.4×102增加至2.3×103時(shí)，優(yōu)化結(jié)構(gòu)變化不明顯，僅一部分細(xì)小分支結(jié)構(gòu)呈現(xiàn)萎縮趨勢，此時(shí)雖然流速有所增加，但仍主要靠導(dǎo)熱傳熱，上述2種情況對(duì)應(yīng)的目標(biāo)函數(shù)值相差不大。然而，當(dāng)Gr從2.3×103增加至1.7×104時(shí)，優(yōu)化結(jié)構(gòu)明顯收縮，形成彎曲而非筆直的分支，以留出周圍空間使得計(jì)算域內(nèi)各個(gè)大渦旋得以發(fā)展，此時(shí)計(jì)算域內(nèi)速度明顯上升，目標(biāo)函數(shù)值明顯下降，對(duì)流傳熱在總傳熱中占據(jù)明顯比重。當(dāng)Gr進(jìn)一步從1.7×104增加至1.2×105時(shí)，優(yōu)化結(jié)構(gòu)也進(jìn)一步呈現(xiàn)收縮趨勢以加速流體流動(dòng)，計(jì)算域內(nèi)流體速度上升，目標(biāo)值下降，此時(shí)傳熱完全轉(zhuǎn)變?yōu)閷?duì)流主導(dǎo)。綜上，隨著Gr增加，傳熱隨著空氣對(duì)流的加強(qiáng)得以不斷強(qiáng)化，對(duì)于優(yōu)化結(jié)構(gòu)而言，需要不斷收縮以適應(yīng)計(jì)算域內(nèi)流體渦旋發(fā)展。由此可知，所得優(yōu)化結(jié)構(gòu)在不同的傳熱主導(dǎo)機(jī)制下呈現(xiàn)出完全不同的形態(tài)特征，且具有鮮明的物理意義。值得指出的是，雖然本文中散熱器是由可滲透的金屬泡沫制成，但從圖5（a）中z=10～20時(shí)的低流速可知，由于本文中金屬泡沫滲透率較低，實(shí)際上流體難以穿透大片區(qū)域的金屬泡沫。盡管如此，金屬泡沫高導(dǎo)熱系數(shù)帶來的強(qiáng)化傳熱效果依然很好。由圖3～圖4還可知，對(duì)于給定的參考結(jié)構(gòu)，計(jì)算域內(nèi)的流體流速同樣隨著Gr增加而不斷增加，傳熱主導(dǎo)機(jī)制也逐漸由導(dǎo)熱傳熱向?qū)α鱾鳠徂D(zhuǎn)換，因而散熱效果明顯增強(qiáng)。進(jìn)一步由表1及圖5（b）可知，參考結(jié)構(gòu)2的散熱性能始終比參考結(jié)構(gòu)1要好，但始終比優(yōu)化結(jié)構(gòu)差。

為了定量評(píng)估優(yōu)化散熱器結(jié)構(gòu)的強(qiáng)化傳熱效果，引入冷壁面的平均努塞爾數(shù)Nuave［22］，定義如下

Nuave=h∫X0SdXkf（Ts－T0）A（38）

式中：kf為流體導(dǎo)熱系數(shù)；A為冷壁面的總面積。

由上式可知，在優(yōu)化問題中設(shè)置最小化平均溫度Ts，本質(zhì)上等價(jià)于最大程度強(qiáng)化換熱（最大化Nuave）。將式（38）統(tǒng)計(jì)得到的Nuave列于表2，其中E表示優(yōu)化結(jié)構(gòu)相對(duì)于參考結(jié)構(gòu)2的Nuave的提升百分比，用于定量評(píng)估優(yōu)化結(jié)構(gòu)的強(qiáng)化傳熱效果。由表2可見，相較于參考結(jié)構(gòu)2，優(yōu)化結(jié)構(gòu)可將傳熱效率至少提高29.7%，相較于參考結(jié)構(gòu)1的提升幅度則更大。這種強(qiáng)化傳熱效果可以解釋如下，當(dāng)Gr較低時(shí)，優(yōu)化結(jié)構(gòu)主要通過形成復(fù)雜的樹枝狀結(jié)構(gòu)，盡可能地將高導(dǎo)熱系數(shù)的金屬泡沫分支伸向冷壁面，以降低發(fā)熱元件與冷壁面間的導(dǎo)熱熱阻。此時(shí)，由于簡單的金屬泡沫方塊和開槽金屬泡沫塊在空間中的延伸寬度和高度均不如優(yōu)化結(jié)構(gòu)，因此散熱性能不如優(yōu)化結(jié)構(gòu)。但隨著Gr增加，優(yōu)化結(jié)構(gòu)不斷收縮。實(shí)際上，當(dāng)Gr增加至1.2×105時(shí)，優(yōu)化結(jié)構(gòu)在空間中的延伸寬度已經(jīng)和開槽金屬泡沫基本相當(dāng)。此時(shí)，優(yōu)化結(jié)構(gòu)實(shí)現(xiàn)強(qiáng)化傳熱的機(jī)理不再是簡單地依靠拓寬散熱器空間延伸寬度和高度來強(qiáng)化導(dǎo)熱，而是依靠分支結(jié)構(gòu)和中心空心結(jié)構(gòu)的配合來優(yōu)化流動(dòng)進(jìn)而強(qiáng)化對(duì)流傳熱。實(shí)際上，相較于整塊的金屬泡沫方塊，開槽金屬泡沫塊具有更好散熱效果的原因就在于開槽金屬泡沫的泡沫塊間存在可供流體受熱后加速的純流體區(qū)，這使得單位時(shí)間內(nèi)有更多流體能更充分地與金屬泡沫換熱，獲得更強(qiáng)的換熱效果。而優(yōu)化散熱器結(jié)構(gòu)則在開槽金屬泡沫的基礎(chǔ)上進(jìn)一步優(yōu)化了金屬泡沫和純流體區(qū)的分布，從而帶來更好的流體加速效果，也就獲得了更強(qiáng)的散熱能力。這種優(yōu)化的主要特征為：在樹枝狀或者花朵狀優(yōu)化結(jié)構(gòu)中，由金屬泡沫樹枝或花瓣包圍而形成了熱源正上方的空心區(qū)域，使得與樹枝或花瓣接觸后受熱加速的流體能在其中暢通無阻地流動(dòng)，從而加快了流體在發(fā)熱元件與冷壁面間循環(huán)的速度，提升了換熱效率。這種周圍密、中心疏的翅片分布已被證實(shí)能有效強(qiáng)化自然對(duì)流傳熱［14］，而在優(yōu)化結(jié)構(gòu)中也存在此種類似結(jié)構(gòu)，進(jìn)一步驗(yàn)證了優(yōu)化所得結(jié)構(gòu)的合理性。此外，由表2可見，強(qiáng)化傳熱提升百分比并不隨Gr的增加而單調(diào)增加，而是在高Gr時(shí)有所回落。這是因?yàn)殡S著Gr增加，參考結(jié)構(gòu)和優(yōu)化結(jié)構(gòu)的換熱均會(huì)增強(qiáng)，而強(qiáng)化傳熱提升百分比能否增加取決于二者增強(qiáng)的相對(duì)比例，因而其數(shù)值未必隨Gr單調(diào)增加。

值得指出的是，本文聚焦單一金屬泡沫材料散熱器的結(jié)構(gòu)設(shè)計(jì)，但在自然對(duì)流系統(tǒng)的多孔介質(zhì)散熱器拓?fù)鋬?yōu)化中，多孔介質(zhì)材料的底層微結(jié)構(gòu)（單胞結(jié)構(gòu)）可選用不同形式，相關(guān)的有效輸運(yùn)系數(shù)等可以通過孔尺度模擬獲得。關(guān)于不同單胞結(jié)構(gòu)、不同底層孔隙率及不同多孔材料體積填充率等對(duì)優(yōu)化結(jié)構(gòu)的影響，感興趣的讀者可以參考課題組近期的多尺度拓?fù)鋬?yōu)化工作［30］。

4"結(jié)論與展望

本文發(fā)展了基于格子Boltzmann方法和水平集方法的三維拓?fù)鋬?yōu)化方法，用于優(yōu)化由金屬泡沫制成的散熱器結(jié)構(gòu)。首先，介紹了用于求解多孔介質(zhì)內(nèi)流動(dòng)和傳熱的格子Boltzmann模型及相關(guān)邊界格式；然后，介紹了基于反應(yīng)-擴(kuò)散方程的水平集方法及相關(guān)的材料物性插值格式；最后，給出了相應(yīng)的伴隨格子Boltzmann模型及梯度計(jì)算式。采用上述發(fā)展的拓?fù)鋬?yōu)化方法，對(duì)底部放置發(fā)熱元件的方腔內(nèi)自然對(duì)流過程進(jìn)行拓?fù)鋬?yōu)化，獲得了4個(gè)特征Gr下的優(yōu)化散熱器結(jié)構(gòu)，結(jié)合2類參考結(jié)構(gòu)定量評(píng)估了優(yōu)化結(jié)構(gòu)的強(qiáng)化傳熱性能，分析了相關(guān)的強(qiáng)化換熱機(jī)理，得到以下主要結(jié)論。

（1）采用的格子Boltzmann方法具有計(jì)算效率高、并行性好、處理復(fù)雜界面能力強(qiáng)等優(yōu)點(diǎn)，基于此發(fā)展的伴隨格子Boltzmann模型同樣具備上述優(yōu)點(diǎn)，這使得針對(duì)復(fù)雜多物理場的大規(guī)模三維拓?fù)鋬?yōu)化結(jié)構(gòu)設(shè)計(jì)成為可能。基于此，本文實(shí)現(xiàn)了金屬泡沫散熱器的三維拓?fù)鋬?yōu)化結(jié)構(gòu)設(shè)計(jì)。

（2）觀察不同Gr下的優(yōu)化結(jié)構(gòu)可以發(fā)現(xiàn)，隨著Gr增加，計(jì)算域內(nèi)的流體流速顯著增加，發(fā)熱元件及散熱器本身的溫度均有明顯下降，散熱主導(dǎo)機(jī)制由熱傳導(dǎo)轉(zhuǎn)變?yōu)闊釋?duì)流，優(yōu)化結(jié)構(gòu)隨之由伸展的樹枝狀向收縮的花朵狀結(jié)構(gòu)變化，以留出更多空間使得流動(dòng)渦旋的發(fā)展更加充分。

（3）選用金屬泡沫方塊結(jié)構(gòu)和帶有開槽的金屬泡沫塊結(jié)構(gòu)作為參考結(jié)構(gòu)，與優(yōu)化設(shè)計(jì)結(jié)構(gòu)進(jìn)行對(duì)比研究后發(fā)現(xiàn)，拓?fù)鋬?yōu)化通過形成分支結(jié)構(gòu)和中心空心結(jié)構(gòu)的配合以優(yōu)化流動(dòng)進(jìn)而強(qiáng)化對(duì)流傳熱。相較于2種參考結(jié)構(gòu)，拓?fù)鋬?yōu)化結(jié)構(gòu)可使散熱效率至少提升29.7%，表現(xiàn)出了優(yōu)異的散熱性能。

本文為金屬泡沫散熱器的優(yōu)化設(shè)計(jì)提供了新方法和新見解，也為新一代輕質(zhì)高效多孔介質(zhì)結(jié)構(gòu)優(yōu)化設(shè)計(jì)提供了理論指導(dǎo)。后續(xù)研究將集中在更加復(fù)雜的多尺度多物理場耦合輸運(yùn)過程的多孔介質(zhì)結(jié)構(gòu)優(yōu)化理論方法和應(yīng)用實(shí)踐等方面，如考慮電池內(nèi)流動(dòng)-擴(kuò)散-反應(yīng)的多孔電極結(jié)構(gòu)拓?fù)鋬?yōu)化。

參考文獻(xiàn)：

［1］ZHAO Changying. Review on thermal transport in high porosity cellular metal foams with open cells ［J］. International Journal of Heat and Mass Transfer， 2012， 55（13/14）： 3618-3632.

［2］MAHDI H， LOPEZ P J R， FUENTES A， et al. Thermal performance of aluminium-foam CPU heat exchangers ［J］. International Journal of Energy Research， 2006， 30（11）： 851-860.

［3］朱繼宏， 周涵， 王創(chuàng)， 等. 面向增材制造的拓?fù)鋬?yōu)化技術(shù)發(fā)展現(xiàn)狀與未來 ［J］. 航空制造技術(shù)， 2020， 63（10）： 24-38.

ZHU Jihong， ZHOU Han， WANG Chuang， et al. Status and future of topology optimization for additive manufacturing ［J］. Aeronautical Manufacturing Technology， 2020， 63（10）： 24-38.

［4］吳璇， 陳群. 基于不動(dòng)點(diǎn)迭代的自然對(duì)流熱沉拓?fù)鋬?yōu)化 ［J］. 工程熱物理學(xué)報(bào)， 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.

［5］吳璇， 陳群. 基于偽密度的對(duì)流換熱熱沉結(jié)構(gòu)拓?fù)鋬?yōu)化 ［J］. 工程熱物理學(xué)報(bào)， 2022， 43（3）： 763-769.

WU Xuan， CHEN Qun. Topology optimization of convection heat sink based on pseudo-density method ［J］. Journal of Engineering Thermophysics， 2022， 43（3）： 763-769.

［6］董馨， 劉小民. 微流控芯片通道結(jié)構(gòu)的拓?fù)鋬?yōu)化研究 ［J］. 西安交通大學(xué)學(xué)報(bào)， 2018， 52（6）： 143-149.

DONG Xin， LIU Xiaomin. Topology optimization of the channel structure in microfluidic chips ［J］. Journal of Xi’an Jiaotong University， 2018， 52（6）： 143-149.

［7］董馨， 劉小民. 流體拓?fù)鋬?yōu)化的參數(shù)選擇與分析 ［J］. 西安交通大學(xué)學(xué)報(bào)， 2018， 52（8）： 132-138.

DONG Xin， LIU Xiaomin. Selection and analysis of parameters for fluid topology optimization ［J］. Journal of Xi’an Jiaotong University， 2018， 52（8）： 132-138.

［8］BORRVALL T， PETERSSON J. Topology optimization of fluids in stokes flow ［J］. International Journal for Numerical Methods in Fluids， 2003， 41（1）： 77-107.

［9］YOON G H. Topological design of heat dissipating structure with forced convective heat transfer ［J］. Journal of Mechanical Science and Technology， 2010， 24（6）： 1225-1233.

［10］DEDE E M. Optimization and design of a multipass branching microchannel heat sink for electronics cooling ［J］. Journal of Electronic Packaging， 2012， 134（4）： 041001.

［11］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.

［12］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.

［13］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.

［14］LI Hao， KONDOH T， JOLIVET P， et al. Optimum design and thermal modeling for 2D and 3D natural convection problems incorporating level set-based topology optimization with body-fitted mesh ［J］. International Journal for Numerical Methods in Engineering， 2022， 123（9）： 1954-1990.

［15］CHEN Li， HE An， ZHAO Jianlin， et al. Pore-scale modeling of complex transport phenomena in porous media ［J］. Progress in Energy and Combustion Science， 2022， 88： 100968.

［16］PINGEN G， EVGRAFOV A， MAUTE K. Topology optimization of flow domains using the lattice Boltzmann method ［J］. Structural and Multidisciplinary Optimization， 2007， 34（6）： 507-524.

［17］KRAUSE M J， THTER G， HEUVELINE V. Adjoint-based fluid flow control and optimisation with lattice Boltzmann methods ［J］. Computers amp; Mathematics with Applications， 2013， 65（6）： 945-960.

［18］YAJI K， YAMADA T， YOSHINO M， et al. Topology optimization using the lattice Boltzmann method incorporating level set boundary expressions ［J］. Journal of Computational Physics， 2014， 274： 158-181.

［19］YAJI K， YAMADA T， YOSHINO M， et al. Topology optimization in thermal-fluid flow using the lattice Boltzmann method ［J］. Journal of Computational Physics， 2016， 307： 355-377.

［20］DUGAST F， FAVENNEC Y， JOSSET C. Reactive fluid flow topology optimization with the multi-relaxation time lattice Boltzmann method and a level-set function ［J］. Journal of Computational Physics， 2020， 409： 109252.

［21］LUO Jiwang， CHEN Li， HE An， et al. Topology optimization of convective heat transfer by the lattice Boltzmann method ［J］. International Journal for Numerical Methods in Fluids， 2023， 95（3）： 421-452.

［22］LUO Jiwang， CHEN Li， KE Hanbing， et al. Three-dimensional topology optimization of natural convection using double multiple-relaxation-time lattice Boltzmann method ［J］. Applied Thermal Engineering， 2024， 236： 121732.

［23］TANABE Y， YAJI K， USHIJIMA K. Topology optimization using the lattice Boltzmann method for unsteady natural convection problems ［J］. Structural and Multidisciplinary Optimization， 2023， 66（5）： 103.

［24］LUO Jiwang， CHEN Li， XIA Yang， et al. Topology optimization of natural convection using porous metal foam based on the adjoint lattice Boltzmann method and level set method ［J］. Computers amp; Fluids， 2023， 265： 106007.

［25］GUO Zhaoli， ZHAO Tianshou. Lattice Boltzmann model for incompressible flows through porous media ［J］. Physical Review E， 2002， 66（3）： 036304.

［26］CHEN Sheng， YANG Bo， ZHENG Chuguang. Simulation of double diffusive convection in fluid-saturated porous media by lattice Boltzmann method ［J］. International Journal of Heat and Mass Transfer， 2017， 108： 1501-1510.

［27］曲東越， 張海兵， 徐建安， 等. 基于改進(jìn)水平集結(jié)構(gòu)拓?fù)鋬?yōu)化方法的理論研究 ［J］. 應(yīng)用力學(xué)學(xué)報(bào)， 2019， 36（4）： 895-900.

QU Dongyue， ZHANG Haibing， XU Jianan， et al. Theoretical research of structure topology optimization based on the level set method ［J］. Chinese Journal of Applied Mechanics， 2019， 36（4）： 895-900.

［28］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.

［29］BHATTACHARYA A， MAHAJAN R L. Metal foam and finned metal foam heat sinks for electronics cooling in buoyancy-induced convection ［J］. Journal of Electronic Packaging， 2006， 128（3）： 259-266.

［30］LUO Jiwang， CHEN Li， XIA Yang， et al. Three-dimensional multi-scale topology optimization of porous heat sink with predetermined unit cells for natural convection heat transfer ［J］. International Journal of Heat and Mass Transfer， 2024， 225： 125398.

（編輯"李慧敏）