芯片級節流制冷器換熱器結構多目標優化研究

童 欣，邱 杰，李家鵬，槐 陽，陳俊元，謝坤圓

（昆明物理研究所，昆明 650223）

0 引言

近年來，芯片行業發展迅速，小型化和集成化成為主要發展趨勢，廣泛應用于生物、光學和IT等領域[1]。隨著芯片體積的縮小，適配的制冷系統也需相應縮小。芯片級節流制冷器的體積比常規節流制冷器小一個數量級左右，最早由美國斯坦福大學Little等[2-4]于20世紀70年代發明。首次采用了微加工技術，在材料表面加工出微米級的細微槽道，構成了高、低壓換熱器以及節流區域，并采用鍵合工藝進行組裝密封。

芯片級節流制冷器在國外已經較為成熟，美國MMR公司推出了一系列制冷量為25～250 mW的產品，主要應用于半導體霍爾效應和塞貝克效應等檢測儀器、CCD陣列芯片以及光學器件中[5]。荷蘭TWENTE 大學 Lerou等[6-9]、Derking等[10-11]和 Haishan等[12-15]也對芯片級節流制冷器進行了廣泛深入的研究，提出了由點陣組成的高、低壓側換熱器結構，如圖1[15]所示。此種結構的MMR可根據不同工況需求調整尺寸和制冷量，以雙級制冷方式達到28 K制冷溫度。MMR能夠與微型熱電制冷器進行耦合，在實現相同制冷量的同時有效降低所需的氣源壓力。

圖1 TWENTE大學芯片級節流制冷器Fig.1 MMR made by TWENTE university

圖1中的逆流換熱器的流動和換熱性能對制冷器整體性能起到至關重要的作用。文獻[8]對圖2所示的不同形狀的換熱器進行了實驗研究，發現菱形和正弦型換熱器擁有最佳的流動性能，然而綜合考慮制造加工難度以及換熱性能，最終采用了圓形點陣的方案，圓形點陣換熱器的進一步優化未見報道。由于MMR結構微小，制造難度較大，數值模擬是優化其性能的有效方法。

圖2 MMR不同形狀的換熱器結構Fig.2 Different shape of heat exchangers of MMR

由于MMR最大尺寸為60 mm，最小尺寸僅為1 μm，差距巨大，無法對整體MMR進行建模。本文對高、低壓側換熱器的若干柱體組成的換熱單元進行參數化建模，使模型的所有幾何參數均為自變量結構參數的函數，生成結構計算實驗點。然后使用CFD-Fluent進行數值計算，得到該結構的流動和換熱性能數據。最后采用CCD優化算法、Kriging響應面以及多目標遺傳算法優化方法對MMR換熱器的努塞爾數Nu以及摩擦因子f進行優化，分別得到高壓側和低壓側換熱器的理論最佳結構，并對最佳結構進行了數值計算驗證。

1 MMR換熱器的計算模型及數值模擬方法

1.1 計算模型的建立

為了對計算結果進行有效驗證，高、低壓換熱器計算模型的建立和幾何規律采用文獻[15]中的結構，如圖3（a）所示。計算模型采用參數化建模，模型自變量結構參數為柱體底部圓直徑Dpillar以及沿流動方向相鄰柱體間距Sspace，當自變量結構參數發生改變時模型也會隨之發生改變，如圖3（b）（c）所示。

圖3 MMR換熱器的計算模型Fig.3 Calculation model of MMR heat exchanger

1.2 基本方程和數值方法

由于MMR換熱器尺寸微小，流體在換熱器內無法充分發展為湍流，因此采用層流模型進行計算，基本方程包括質量、動量和能量方程，表達式如下。

連續性方程：

式中：v為速度矢量，m·s-1；T為靜態溫度，K；ρ為流體密度，kg·m-3；p為靜壓，Pa；μ為流體分子黏度，Pa·s；cp為流體比熱，J·kg-1·K-1；λ為流體導熱系數，W·m-1·K-1；?為湍動能耗散項。在層流模型中?可忽略不計。

計算模型物性參數和邊界條件設置為：高、低壓側換熱器流體均采用氮氣作為工質，物性在Refprop中查取，換熱器中流動計算過程按照工質物性恒定進行；高壓側換熱器進口氮氣設置為300 K以及8 MPa，低壓側換熱器進口氮氣按中波紅外探測器的典型工作溫度97 K以及壓力0.6 MPa進行設置；因模型大小會隨自變量結構參數改變，計算模型進口處采用速度進口邊界條件以保證在模型變換時進口雷諾數Re變化不大，高壓側換熱器進口速度為0.57 m/s，低壓側換熱器進口速度為0.9 m/s；模型出口采用背壓為0的壓力出口邊界條件；高壓側換熱器采用上接觸面和柱體表面定壁溫，低壓側換熱器采用下接觸面和柱體表面定壁溫邊界條件，以模擬實際MMR工作時高、低壓側換熱器通過中間制冷器片進行導熱換熱的狀態；壁溫與流體溫差設置為5 K，即高壓側換熱壁面為295 K，低壓側換熱壁面為102 K，其余壁面均按絕熱面處理；所有壁面均設置為無滑移標準壁面函數，壓力-速度耦合采用SIMPLE算法，動量和能量方程均為二階迎風格式，各項殘差均設置為1×10-6。采用以上邊界條件計算得到的高壓側換熱器中部平面溫度分布云圖如圖4所示（模型結構參數見3.1節中的初始結構），模型最右邊溫度出現分布不均區域，原因為數值計算時出口處壓力計算值與背壓存在少量不匹配引發的回流現象，該區域面積較小對換熱器整體性能評價影響不大。

圖4 高壓側換熱器中部溫度分布云圖Fig.4 Temperature contour of mid plane in high-pressure heat exchanger

1.3 計算模型的驗證

計算模型在進行優化算法前需進行網格無關性驗證。高、低壓側換熱器計算模型網格均為四面體網格，采用MESH生成，高壓側換熱器生成網格數量為60萬、100萬、200萬和360萬，低壓側換熱器生成網格數量為110萬、160萬、210萬和440萬。不同網格模型的Nu和f采用與文獻[14]中相同的計算方法，并與在相同Re下該文獻提出的關聯式計算值Nucor、fcor進行對比，如圖5所示。

圖5 高、低壓側換熱器計算模型網格無關性驗證Fig.5 Grid independence verification of calculation model of high and low pressure heat exchanger

由計算結果可見，隨著網格數量增加，高壓以及低壓側換熱器模型Nu和f的變化均趨于平緩，表明當網格數量足夠密集時，網格變化對計算結果影響顯著減小，而過于密集的網格會導致計算速率大幅下降。高壓側換熱器模型計算結果隨網格數量變化更顯著，應選擇相對較多的網格，而低壓側換熱器模型計算結果隨網格數量變化不大，可選擇相對較少的網格。綜合考慮計算準確性以及計算速率，高壓側選擇200萬網格模型，低壓側選擇160萬網格模型，兩模型的Nu和f與文獻[14]中關聯式的計算偏差分別為：高壓側4.0%、3.8%；低壓側5.1%、1.5%，能夠滿足計算精度要求。

網格無關性驗證后還需對模型在各個工況下的計算精度進行驗證，文獻[14]高低壓側換熱器性能關聯式適用范圍為Re≤300，因此分別對高低壓側換熱器模型工況在Re=100～300進行計算驗證，結果如圖6所示。計算結果與關聯式符合良好，高壓側換熱器Nu的最大誤差5.2%、平均誤差2.3%，f最大誤差3.5%、平均誤差2.3%；低壓側換熱器Nu最大誤差6.1%、平均誤差5.3%，f最大誤差1.8%、平均誤差1.5%。造成以上計算誤差的原因可能是網格劃分以及工質物性設置的不同。

圖6 高、低壓換熱器計算模型驗證Fig.6 Calculation model verification of high and low pressure heat exchanger

2 數值計算

2.1 計算模型幾何參數

高、低壓側換熱器幾何參數及進口條件設置如表1所列，其中柱體底部直徑Dpillar以及沿流動方向相鄰柱體間距Sspace為模型自變量。自變量變化范圍需在一定合理范圍內，否則會導致在高壓工況下制冷器片鍵合失效，造成結構損壞。高壓側換熱器柱體為交叉排布以提高換熱能力，低壓側換熱器柱體為垂直排布以減小流動阻力，降低蒸發腔背壓從而降低制冷溫度，兩模型換熱器柱體排列形式如圖3（b）（c）。

表1 高、低壓側換熱器計算模型幾何參數及進口條件Tab.1 Geometric parameters and inlet conditions of calculation model of high and low pressure heat exchanger

2.2 數據處理

本文采用努塞爾數Nu和摩擦因子f來表征高低壓側換熱器的換熱和流動性能，Nu越高換熱性能越好，f越低換熱器阻力越小，流動性能也越好。模型的換熱性能參數Nu計算公式如下：

式中：Dh為水力直徑，m；Vflow為流通體積，m3；h為換熱系數，W·m-2·K-1；Q為總換熱量（為計算模型中所有換熱壁面的換熱量總和），W；Atrans為換熱面積（柱體表面積+上或下換熱面），m2；Δtm為對數平均溫差，K；Tin、Tout、Tw分別為進口、出口以及壁面平均溫度，K。

式（8）～（10）主要計算f以及Re。

式中：Δp為模型壓降，Pa；pin、pout為進口和出口平均壓力，Pa；v為進口流速，m·s-1；L為換熱段長度（8、16列柱體排列總長度），m；μ為流體分子黏度，Pa·s。

3 計算結果分析與討論

3.1 響應面方法和多目標參數優化理論

在得到模型初始結構的換熱和壓力性能（高壓側換熱器Dpillar=0.28 mm、Sspace=0.03 mm；低壓側換熱器Dpillar=0.13 mm、Sspace=0.025 mm）后對高低壓模型的Dpillar和Sspace進行改變，變化范圍如表1所列。通過計算點設計試驗計算量較小、精度高和預測性好的中心組合設計CCD（Central Composite Designs）方法[16]，從而得到計算點用于建立響應面。該方法更加強調變量邊界處的取值計算，相比于傳統的正交試驗法所需的計算點更少、計算速率更快，更加有利于體現影響因素之間的交互作用。

響應面采用Kriging插值方法[17]，是一種通過已知點來預測未知觀察點的插值方法，利用方差的變化來表達空間的變化，可以保證由空間分布得到的預測值的誤差最小[18]。

利用Kriging響應面中的數據進行換熱器結構多目標優化，在一定的設計參數范圍內尋找同時滿足Nu最大化和f最小化的換熱器結構[19-21]。

3.2 響應面計算結果及分析

CCD優化算法及Kriging響應面生成的高、低壓側換熱器的Nu、f關于結構變量Dpillar、Sspace的響應面如圖7和圖9所示，從響應面中能夠大致分析出Nu、f關于結構參數的變化趨勢。

圖7 高壓側換熱器性能-結構響應面Fig.7 High-pressure heat exchanger performance-structure response surfaces

由圖7（a）可見，高壓側換熱器中Dpillar較小區域的Nu較高，在此區域較大的Sspace也能在一定范圍內提高Nu。這是因為MMR換熱器流道十分微小，工質流體呈層流流動，較小的Dpillar能夠使得流體更易繞柱體環流流動，較大的Sspace使得層流能夠更充分地發展，從而提升換熱性能。值得注意的是在Sspace=0.03 mm附近Nu變化規律產生改變，此時Nu受Dpillar變化的影響更大并且在Dpillar=0.28 mm附近Nu值快速減小，圖7（a）出現凹陷區域。這可能是因為Sspace在0.03 mm附近時層流流型發生較顯著改變，此時Dpillar改變對換熱器換熱性能影響更大。在Sspace=0.03 mm時，Dpillar=0.24～0.28 mm區間Nu隨Dpillar的增加減小較快，而在Dpillar=0.28～0.32 mm區間Nu隨Dpillar的增加而減小的斜率較低。為了更好地對此現象進行分析，分別對Sspace=0.03 mm，Dpillar=0.24、0.26、0.28、0.30、0.32 mm五個結構進行單獨研究，提取每個模型流道中部截面進口區域流體速度分布云圖，如圖8所示。可見隨著Dpillar增加流體流速也隨之增加，并且在柱體后方逐漸出現渦旋（圖中方框所示），渦旋處流體滯留，在增加流體阻力的同時也對換熱產生不利影響，在Dpillar=0.24～0.28 mm階段渦旋逐漸出現并充分發展，Dpillar＞0.28 mm后渦旋無顯著變化。因此Sspace=0.03 mm時Nu隨Dpillar增大而變化的速率先快后慢。

圖8 高壓側換熱器不同Dpillar速度分布云圖Fig.8 Velocity contours of high-pressure heat exchanger with different Dpillar

由圖 7（b）可見，較小的 Dpillar和 Sspace使f顯著降低，當柱體更小、柱體間距更大時換熱通道對流體的阻力也將減小。與之類似，在低壓側換熱器中，如圖9（a）（b）所示，在低壓側換熱器中隨著 Dpillar減小Nu顯著增加，并且受Sspace影響不大。這是因為低壓換熱器的柱體排列行間距較大，為2倍Dpillar+Sspace，因此Sspace的改變對層流流型無顯著影響。而較小的Dpillar和較大的Sspace也有利于降低換熱器的摩擦因子f，并且相較于高壓側換熱器，低壓側換熱器中Sspace對f的影響更為顯著，原因是由于低壓側換熱器柱體陣列呈垂直排布，且柱體行間距大，換熱器的流動阻力主要來自平行于流動方向的柱體間隙即Sspace。

圖9 低壓側換熱器性能-結構響應面Fig.9 Low-pressure heat exchanger performance-structure response surfaces

通過響應面計算能夠得出MMR換熱器換熱和流動性能關于幾何結構的變化規律，對換熱器優化起到一定指導作用。并且根據響應面數據可同時對兩個目標函數進行多目標優化。

3.3 多目標優化

換熱器優化目標函數應為：努塞爾數Nu盡可能大、摩擦因子f盡可能小，這兩個目標函數在某些時候存在沖突，即換熱較佳的結構通常會伴有較大流動阻力。因此優化過程需首先利用響應面中數據獲得Nu和f的所有非支配解的解集，即帕累托前沿（Pareto front），響應面中不存在Nu和f性能同時優于帕累托前沿中解的結構點。

自變量結構參數變化如表1所列，其中高壓側換熱器Dpillar=0.28 mm、Sspace=0.03 mm以及低壓側換熱器Dpillar=0.13 mm、Sspace=0.025 mm為初始結構。利用響應面得到換熱器的帕累托前沿如圖10所示。

由圖10可見，在所有優化結構中，更高Nu的換熱器結構也伴隨著更高的f，為在帕累托前沿中選取更加合理的優化結構點，采用權重系數k來評價Nu和f在優化過程中的相對重要性[22]。其計算如下：

圖10 計算模型優化算法帕累托前沿Fig.10 Pareto optimal fronts for calculation models

函數F（x）反映了在Nu和f的帕累托前沿變化內各備選優化點性能參數與前沿中最優性能參數Numax和fmin的偏差大小，F（x）值最小時即為Nu和f在當前權重系數下同時接近最佳性能的結構參數點。

換熱器分別選取權重系數k=0.1、0.2～0.9范圍內得到的優化結構參數和性能以及與原始結構對比如表2所列。

由表2可見，取不同權重系數k時，高、低壓側換熱器性能較原始結構均得到了優化。根據式（11）定義，k值越高優化越傾向于強化Nu性能，反之則傾向于強化f性能。對于MMR制冷器而言，較高換熱性能有利于提升降溫速率，但MMR一般為非金屬材料，無法通過不斷提高供氣壓力來提升制冷性能，因此為保證制冷器有足夠的節流制冷效應，流動阻力也不宜過大。根據MMR工作特性，對換熱器換熱和流動性能進行均衡考慮，因此選擇k=0.5時的結構參數作為優化結構，并加以進一步驗證和分析。

表2 換熱器不同權重系數優化結果Tab.2 Optimization results of heat exchanger with different weighting factor

對k=0.5時的優化結構進行CFD重新計算，誤差在1%以內，證明了多目標參數優化的準確性。本文選取高壓側換熱器的最佳結構點為Dpillar=0.24 mm、Sspace=0.033 mm，較原始結構Nu提升了5.8%，f下降了25.7%。低壓側換熱器最佳結構點為Dpillar=0.10 mm、Sspace=0.032 mm，較原始結構Nu提升了14.5%，f下降了2.3%。為進一步分析最佳結構與原始結構流動和換熱的區別，分別選取兩結構中部進口部分平面的速度和溫度分布如圖11、12所示。

由圖11（a）可見，高壓側換熱器優化結構的流體流速相對較慢，并且柱體后方的渦旋更少（圖中方框所示），有利于減小流動阻力和強化換熱。圖11（b）進一步證明該觀點，在方框所示區域優化結構的流體溫度相對更低，說明具有更高換熱效率。

圖11 高壓側換熱器速度和溫度分布云圖Fig.11 Velocity and temperature contours of high-pressure heat exchanger

由圖12可見，與高壓側換熱器類似，低壓側換熱器優化結構中流體流速相對較低，并且流體溫度升高相同值所需的流動距離變短（圖（b）中虛線所示）。

圖12 低壓側換熱器速度和溫度分布云圖Fig.12 Velocity and temperature contours of low-pressure heat exchanger

綜上所述，高、低壓側換熱器的優化結構較原始結構在換熱和流動性能上得到提升。進一步對比優化結構與原始結構參數可得出以下結論：

對于高壓側換熱器而言，較小的Dpillar能夠提高換熱性能，降低流動阻力；Sspace值過大時會使得換熱器柱體陣列過于稀疏，從而影響換熱性能。

對于低壓側換熱器而言，較小的Dpillar能提升換熱性能，對流動性能優化效果不明顯。這是由于低壓換熱器的柱體陣列垂直排布以及較大行間距已經在很大程度上降低了流動阻力，進一步優化空間有限。該優化結構分別位于模型變量中Dpillar的下限以及Sspace的上限，本文未拓寬自變量取值范圍做進一步計算，這是因為在實際MMR換熱器中低壓換熱器的Dpillar不能過小，Sspace也不能過大，否則會降低換熱通道表面鍵合強度以及使制冷器承受更大應力，在高工作壓力下更易損壞。因此計算所得優化結構應該在MMR換熱器結構參數合理變化范圍內。

通過多目標參數優化分別得到了高、低壓側換熱器在層流工況下的最佳結構點，并且優化后換熱器的換熱和壓力性能均得到了提升。

4 結論

本文采用多目標參數優化數值模擬的方法，對結構參數影響MMR換熱器的傳熱和流動性能進行了研究，并對換熱器結構進行了優化，得到結論：

（1）參數化建模以及多目標參數優化可利用較少的計算資源實現MMR換熱器的多參數變量計算以及多目標函數優化，為MMR的設計提供參考。

（2）對于高、低壓側換熱器而言，較小的柱體底部直徑Dpillar有利于提升換熱性能及降低流動阻力。因此設計中在保證結構強度的前提下適當降低Dpillar以提升換熱性能。

（3）對于高壓側換熱器，適當增加柱體間距Sspace使層流更充分地發展以提高換熱和壓力性能，Sspace存在最佳值。對于低壓側換熱器，較大的Sspace僅對換熱性能存在有利影響，對壓力性能影響不大。

（4）利用多目標參數優化方法得到了高、低壓側換熱器在Re=150、210附近工況以及結構參數在表1所列變化范圍內的流動、換熱性能最佳結構，分別為：Dpillar=0.24 mm、Sspace=0.033 mm以及Dpillar=0.10 mm、Sspace=0.032 mm；Nu分別提升了 5.8%、14.5%，f分別降低了25.7%、2.3%。