印刷板式微通道換熱器流動與傳熱特性的理論模型

余智強，吳建澤，任亞濤，何明鍵，于喜奎，齊宏

（1 哈爾濱工業大學能源科學與工程學院，黑龍江 哈爾濱 150001； 2 空天熱物理工業和信息化部重點實驗室，黑龍江 哈爾濱 150001； 3 中國航空工業集團公司沈陽飛機設計研究所，遼寧 沈陽 110035）

引 言

隨著世界能源資源減少、能源成本增加和環境問題日益突出，節約能源成為可持續發展目標的關鍵之一。根據美國能源信息管理局（EIA）的數據，到2040 年時世界能源消耗將增加56%。無論是化石能源還是可再生能源，都需要更新、開發更加高效的能源技術與設備[1]。在創新應用中，高效、緊湊型印刷板式換熱器（PCHE）在熱管理系統整體效率中發揮著至關重要的作用[2]。

1981 年，Tuckerman 等[3]首先提出了微通道散熱器的概念，并預測其能以1000 W/m2的速率散去熱量。微通道是通過微變形、微銑削、激光加工、化學刻蝕等技術，在單層金屬薄板上加工出的微小的流道[4]。由于流道加工過程與印刷電路板類似，此類換熱器被稱為“ 印刷板式微通道換熱器”（PCHE）[4-5]。相比于傳統換熱器，PCHE 具有傳熱系數高、體積小、緊湊度高、質量輕等特點，在生物工程、微加工系統、微泵、微熱管等領域有廣泛應用[6]。在“微尺度效應”[7]的作用下，傳統換熱器特性的計算方法變得不再適用，很多研究者針對這一問題展開研究，建立相關的準則關聯式來預測PCHE 的流動與傳熱特性。

在對微通道內流動特性的研究中，Muzychka等[8]對非圓管內流體力學特性進行了詳細的分析，提出以橫截面積的平方根作為特征尺寸比水力直徑更合適，同時構建新的模型并簡化了非圓管f-Re的預測，基于實驗數據獲得了很好的相關性。Hun[9]采用Heatric UK 制造的PCHE 開展實驗研究，提出Fanning 摩擦系數和Nusselt 數與平均Reynolds 數-Re呈線性關系，指出Re的函數依賴于通道中尖角彎頭處的形狀損失；通過數值模擬和實驗測試數據對Fanning 摩擦系數關聯式進行修正，給出在不同角度、節距和水利直徑下關聯式的修正因子。Chen[10]通過實驗研究了換熱器兩側壓損，提出壓降主要由兩部分產生：換熱器核芯的壓降，分配裝置（如入口/出口和流量分配集管）的壓降，在Hun[9]研究的基礎上根據換熱器內部動量變化、摩擦作用、突擴、突縮等因素的影響提出了PCHE 壓力損失的計算模型，對比發現實驗測量結果與模型之間差異是由于通道變化造成的，并且在層流區域壓降強烈依賴于流道的幾何形狀和尺寸。Chen 等[11]采用以Sharifi 等[12]和Chen[10]為基準的動態模型，對縮小尺度的直流道PCHE 的穩態和瞬態行為進行了數值預測，通過實驗驗證了該動力學模型在He/He條件下預測直通道PCHE性能的適用性。Bahrami等[13]基于橢圓通道提出了一種新的近似解，用于測定任意截面單連通微通道中充分發展的層流單相流的壓降，該模型通過矩形、梯形和三角形微通道的實驗數據進行了驗證，還與不同橫截面形狀（超橢圓、梯形、正弦、菱形、圓形扇形等）通道的數值結果進行了比較，預測壓降誤差在8%以內。

在對微通道內換熱特性的研究中，Kim 等[14]采用光刻技術加工出兩個不同結構的PCHE，以水作為冷、熱介質研究不同Reynolds數下的換熱特性，提出用修正的威爾遜圖法計算熱側和冷側傳熱系數的關聯式，與實驗結果相比呈現出了很好的精度，誤差在7%以內。Park 等[15]采用離散化方法得到了沿PCHE 徑向的局部Nusselt 數，并采用不同的數據簡化方法與PCHE 實驗結果進行了比較，發現在設計PCHE 時，基于入口和出口點平均值的典型數據簡化方法不適用于性能發生顯著變化的傳熱過程，離散化方法是反映PCHE 特性變化的有效數據縮減方法。Yoon 等[16]針對高溫氣冷堆中PCHE 展開研究，發現根據實驗數據得出換熱器的效率是被高估的，在評估高溫換熱器性能時引入熱損失部分，提高計算精度。Mortean 等[17]針對方形通道的交錯流換熱器展開研究，基于前人研究提出了一個數學模型來預測換熱器的傳熱特性，對比分析實驗數據和數學模型理論計算結果，發現Stephan 等[18]的關聯式在模型中與實驗數據比較時表現出最好的結果，特別是在低Reynolds 數時。Mortean 等[19]對不銹鋼PCHE 的傳熱和壓降進行了完整的分析，提出了一個預測換熱器壓降的模型，以水為兩側工作流體進行了72 次測試，發現Muzychka 等[8]提出的關聯式與實驗結果吻合較好；提出了一個壓降模型，計算了設備芯、集管等部件的貢獻，得出進口和出口集管對總壓降的貢獻最大，占50%以上。Sarmiento等[2,20]在已有研究的基礎上，以文獻中1482個數據點作為數據庫，使用Churchill 等[21]提出的方法建立了一個適用于大多數常見非圓形管道的模型，模型計算結果與實驗數據的最大誤差小于30%；針對半圓形微通道提出了一種新的方法來考慮通道間翅片幾何形狀的變化，同時引入Yoon 等[16]提出的熱損失概念，建立了PCHE 換熱特性的預測模型，并將模型結果、實驗結果、數值結果進行比較，結果表明三者吻合較好，綜合誤差小于8%。

由上述研究可以發現在流動與傳熱特性準則關聯式的研究中，前者的研究主要針對摩擦特性與流動狀態之間的關系，目的是獲得求解Fanning 摩擦系數f的關聯式，并通過摩擦系數間接得到壓損；而后者的研究主要針對傳熱特性與通道結構、流動狀態之間的關系，通常是獲得求解Nusselt 數的計算關聯式，然后間接計算對流傳熱系數、換熱量等特性參數。然而，大量的研究都是針對特定工質、結構、工況等條件展開，得到的計算關聯式大多只適用于特定情況，更換其中某一條件可能導致其不再適用。

因此，本文基于目前已有的研究，篩選出適用于層流換熱的流動與傳熱特性關聯式，并建立相關計算模型，通過文獻中已公布的不同結構、工質、工況等條件下的實驗測試結果與計算模型中不同關聯式所得經驗解進行對比，獲得其中適用性最高的計算關聯式組成計算模型。最后，提出一種PCHE換熱核芯結構，對其進行兩種工質（航空燃油RP-3與航空潤滑油4050）對流換熱的數值仿真與經驗計算，對比數值解與經驗解之間的精度，進一步驗證計算模型的適用性。

1 理論模型

1.1 幾何結構

對于PCHE換熱核芯（圖1），定義其整體結構尺寸為長L、寬W、高H；由Nc層冷板和Nh層熱板層疊而成，每層板含有通道數M；冷、熱通道的高為b、寬為w、柵格厚e、基底厚a，冷熱流道結構相同，使用下角標c、h區分冷側、熱側。對于本文PCHE幾何模型所采用的半圓形流道，使用Sarmiento等[20]和Liu等[22]提出的理論將其等效為矩形結構，具體方法將在后文中詳細敘述。

圖1 PCHE結構示意圖[2]Fig.1 Structure of PCHE[2]

基于PCHE 的幾何參數，根據文獻[23]，緊湊度α是傳熱面積與體積的比值：

式中，Atotal為總傳熱面積，m2；Vtotal為換熱器核心體積，m3；P為濕周，m。

孔隙度σ為自由流面積Afree與前緣面積的比值，表達式為[23]：

由式（3）可知，換熱器的緊湊度α與孔隙度σ成正比，與當量直徑d成反比?？梢?，為了提高緊湊度α（即更小的換熱體積內具有更大的傳熱面積），就必須減少當量直徑或增加孔隙率。

1.2 傳熱特性理論計算

對換熱器內的換熱過程進行簡化處理來確定傳熱路徑，以便進一步分析和討論，有如下假設：微通道尺寸較小，忽略通道內的輻射和自然對流；固體具有恒定特性，流體不可壓縮、性質無突變，且能以平均溫度為基礎表示；層流流動下恒定傳熱速率、穩態傳熱；流體均勻分布于各通道之間，具有均勻的流速剖面。

基于上述假設，熱交換器的總傳熱系數UA（或總熱阻1/UA）可以表示為[23]：

式中，Rc為對流熱阻；Rd為污垢熱阻；Rwall為基底導熱熱阻；R″d為污垢系數；ηo表示表面總效率；Awall表示平均壁面面積；h表示材料的傳熱系數；ks表示材料熱導率。換熱器中傳熱特性計算簡化為翅片結構，如圖2所示。

圖2 傳熱計算簡化示意圖Fig.2 Simplified schematic diagram of heat transfer calculation

PCHE 的應用中，由于其通道孔徑很小，需要保持工質潔凈，不能有摻混雜質、局部結渣等問題造成微通道堵塞。在計算過程中，可以忽略污垢熱阻的影響，式(4)可以簡化為：

式中，χ表示效率參數。

式中，φ表示通道的寬高比，定義為通道最大邊與最小邊的比值；Red表示以當量直徑為特征尺寸的Reynolds數，表達式為：

式中，ρ為流體密度；V為流速；μ為動力黏度；m?為流體質量流量。

上述表達式中，壁面熱流和通道徑向溫度被認為是常數，是目前緊湊型熱交換器的較好近似。經過上述過程可推導得到關于總傳熱系數UA（或總熱阻1/UA）的求解式，其中重點需要求解的冷、熱流體的對流傳熱系數（或Nusselt數），該值需要通過大量實驗數據推導計算獲得，許多研究者根據不同流動狀態（層流區、過渡區、湍流區）提出了不同的Nusselt 數關聯式。針對微通道內層流具有代表性的關聯式如表1所示。

表 1 微通道內層流換熱Nusselt數關聯式Table 1 Nusselt number correlation of laminar heat transfer in microchannels

表1中8個Nu關聯式均以當量直徑作為特征尺寸，除此之外部分學者研究發現以截面面積的平方根A作為特征尺寸擬合得到的Nusselt數關聯式模型同樣可以取得極佳的計算精度，部分學者研究結果如下。

Muzychka等[8,30]評估了幾種特征尺寸，對橢圓和矩形通道進行了比較，發現圓形通道和矩形通道的解基本上已經坍縮為單一值，并且在小長徑比的極限下，橢圓通道和矩形通道的計算結果也趨于一致，可以清楚地看出，對于層流數據的無量綱化，橫截面積的平方根比水力直徑更合適。Muzychka 等基于截面面積的平方根A作為特征尺寸，提出Nusselt數關聯式模型：

層流狀態下的Fanning摩擦系數表達式為：

此外，常數C1和C3是熱邊界系數，而C2和C4取決于Nusselt數是局部值和平均值；f(Pr)是一個依賴Prandtl數的因子；γ是形狀參數（表2）。

表2 關聯式系數［20，30］Table 2 Correlation coefficient［20，30］

Sarmiento 等[20]使用Churchill等[21,31]提出的方法，建立了一個適用于大多數常見的非圓形管道的模型，該模型利用層流和紊流的換熱表達式生成了一個過渡方程，表達式為：

求得冷、熱側流體Nu后，間接推導對流傳熱系數h，然后再計算得到總傳熱系數UA，最后計算總換熱量等性能參數區分與轉換?；诒疚牡难芯績热?，采用LMTD 法、ε-NTU 法構建計算模型，介紹如下。

換熱量是通過對數平均溫差ΔTLM來確定，適用于熱物性不發生較大變化且工況遠離臨界條件的情況[15,38-39]，當熱物性會發生劇烈變化時，需要采用離散化方法[15,40-41]。本文屬于前者，總熱導率可表示為：

對于校正因子F的計算，不同流動形式的換熱器取值有所不同，完全逆流的換熱器的F值取1。對于混合交錯流換熱器，校正因子F的定義式為：

式中，C表示換熱工質的熱容量；Cmin表示冷、熱工質熱容量的最小值。

恒壓下，工質比定壓熱容cp變化不劇烈時可視為常數，在忽略環境熱量散失時總傳熱量Q可表示為：

式中，C*表示冷、熱工質的熱容量的最大值與最小值之比。不同流道布置形式效能ε表達式有所不同，本文主要涉及準逆流、非混合交錯流。

逆流換熱器效能ε的表達式[42]為：

根據Kuppan[37]、Shah 等[43]提出的解析表達式，非混合交錯流換熱器的效能ε為：

在換熱器運行時，由于環境的影響其熱損失可能很大。Yoon 等[16]考慮周圍環境熱損失提出了模型：

式中，λ為熱損失比；T∞為環境溫度；Taverage,c和Taverage,h分別為冷流體和熱流體的平均溫度。

1.3 流動特性理論計算

對于微通道換熱器的流動特性的計算，本文基于現有文獻中提出的相關實驗關聯式估算的換熱器整體的Fanning 摩擦系數，然后通過摩擦系數間接求解得到換熱器各個關鍵位置產生的壓損。在進行壓降計算前，進行如下假設：

（1）入口流量穩定；

（2）流體在每個通道中均勻分布；

（3）所有通道的幾何形狀相同；

（4）重力效應被忽略，換熱器是水平定向的；

（5）流道通道的橫截面為半圓形。

圖3 為PCHE 在實驗測量過程中可能會存在的壓損分布示意圖。由圖可得，當前PCHE 中的壓降大致可以分為三部分：①換熱器核芯中的流動壓降Δpcore；②通道入口突縮和出口突擴相關的壓降Δpcontraction、Δpexpansion；③進口和出口流量分配集管中的壓降Δpheader。

圖3 壓損分布示意圖Fig.3 Schematic diagram of pressure loss distribution

其中，對于封頭內部的壓損暫沒有相關實驗數據，因此暫時忽略其影響，只針對換熱核芯內的流動特性展開分析與討論。所以PCHE 熱側和冷側總壓降可計算為：

相關文獻中PCHE 各部分的壓力損失理論計算方法大致如下。

（1）換熱器核芯中的流動壓降

換熱器核心中的壓力損失由兩部分組成：摩擦產生的壓力損失Δpfriction、動量變化產生的壓力損失Δpacceleration和通道折角產生的壓力損失Δpangle。其中，摩擦損失同時考慮了表面摩擦和形狀阻力效應，由流動面積變化引起的內部收縮和膨脹產生的任何額外損失（如有）也被納入換熱核芯摩擦損失項，表達式為：

式中，V表示具體體積；Vin和Vout分別在入口和出口溫度和壓力下進行評估。下角標in、out分別表示入口、出口；m代表平均。

Fanning 摩擦系數f基于有效壁剪切應力計算，由于通道內部的剪切應力無法通過實驗等方法直接測量，使得對于微通道內Fanning 摩擦系數f的計算關系式沒有統一的、普適性的實驗關聯式。許多學者基于現有的理論進行研究，致力于通過實驗測量來擬合實驗關聯式來預測通道內部的摩擦損失，本文基于現有的研究內容，整理出一些應用相對較廣的關于層流流動的實驗關聯式，如表3所示。

表3 微通道Fanning摩擦系數f的關系式Table 3 The Fanning friction coefficient f of the microchannel

由于換熱核芯內的動量變化產生的壓縮歸因于流體被加熱或冷卻，導致換熱核芯內的動量率變化或流動加速或減速效應。其中，正值表示流量加速時的壓力降低，負值表示流量減速時的壓力升高。由此產生的壓降由式（35）給出：

流體經過一個彎曲通道時，必然會在流體上作用一個向內的徑向力，外壁面壓力增大，在內壁面上有壓損。在流動慣性的作用下，使得彎管和折管內側往往會產生旋渦，造成局部壓力損失，如圖4所示。

圖4 通道內部折角Fig.4 Internal angle of channel

根據Darcy-Weisbach 公式，折管內局部阻力系數ζangle表達式為[46]：

（2）入口突縮和出口突擴相關的壓降Δpcontraction、Δpexpansion

Kays等[47]分析評估了換熱器核心區域的壓降收縮和膨脹產生的影響，人為地將換熱器的入口段的壓降分成兩個部分：①流動截面積發生變化導致壓力損失，此過程不考慮摩擦作用；②突縮段的不可逆自由膨脹引起的壓損（可以解釋為：流體流過突然收縮斷面，產生流動邊界層分離，導致隨著斷面下游的速度分布發生變化，然后動量速率也發生變化，進而引起相應應力的變化，從而產生壓力損失）。入口段壓損僅因面積變化引起的壓降由式（39）給出：

式中，G是換熱器微通道（或管束）內質量速度，即核芯內質量流動速度；σ是核芯自由流動面積和橫截面積之比，即孔隙率，也稱收縮比。

對換熱核芯入口壓降的第二個貢獻是與突然收縮和動量率變化后的自由膨脹相關的損失，收縮損失系數Kc乘以換熱核芯進口處的動態速度頭，可算出不可逆壓降分量：

式中，Kc為收縮損失系數，是收縮比σ、Reynolds數Re和流動截面幾何形狀的函數。因此，換熱核芯入口壓降為：

對換熱核芯入口壓降的第二個貢獻是與突然收縮和動量率變化后的自由膨脹相關的損失，收縮損失系數Kc乘以換熱核芯進口處的動態速度頭，可算出不可逆壓降分量。

類似地，出口段的壓損也可以分成兩個部分：①流動截面積發生變化而引起的壓損，此過程不考慮摩擦作用，其表達形式與入口側壓損相同；②突擴段的不可逆自由膨脹、動量變化等現象引起的壓損[47]。出口段由于截面面積增加導致減速引起的壓損升高為：

不可逆自由膨脹和突然膨脹后動量率變化相關的壓力損失，由出口損失系數Ke乘以換熱核芯出口處的動態速度頭得出：

式中，Ke是突擴（或出口）壓力損失系數，為膨脹比1/σ、Reynolds 數Re和流動截面幾何形狀的函數。因此，換熱核芯出口處的壓損為：

Kc和Ke是收縮段和擴展段幾何結構的函數，在某些情況下還與管內Reynolds 數有關，原因在于管內速度的變化會引起動量變化，同時入口和出口的動量也會發生變化。在實際應用時，Kc和Ke中已經包含了速度分布變化與壓力變化的影響，Kays 等[47]通過分析大量的實驗結果證實了經驗解的可靠性。Lee[48]提出了計算該系數的模型，表達式為：

在換熱器的流動與傳熱特性計算過程中，通過上述實驗關聯式計算平均Nusselt 數與Fanning 摩擦系數，從而間接確定冷、熱側的對流傳熱系數h，然后通過進一步確定換熱器整體的總傳熱系數UA（或總熱阻1/UA），最后計算得到PCHE 的總換熱量、效能、壓損等特性參數。

2 精度分析

2.1 計算模型的構建

基于上述關于PCHE 中傳熱特性的理論計算，本文通過整理分析設計了一個用于預測換熱器傳熱特性的計算模型，并且基于Matlab 語言編寫了一個計算換熱器傳熱特性的計算程序，框圖如圖5所示。

圖5 程序框圖Fig.5 Flow chart

由圖5 可以看出傳熱特性計算過程，本文建立計算模型的基本流程為：

（1）輸入數據。工質入口溫度Tin、壓力p、流量qm以及換熱器整體的幾何尺寸（H，w，L，a，e，等）；

（2）提出一個效能εas的假設值作為初始值，估算出口溫度Tout；

（3）基于入口、出口溫度確定特征溫度，進一步計算流動工質的物性參數；

（5）基于平均Nusselt 數的計算關聯式估算冷、熱側換熱工質的Nu；

（6）通過Nusselt 數確定冷側、熱側對流傳熱系數h；

（7）計算換熱器的總傳熱系數UA（或總熱阻1/UA）、傳熱單元NTU和總換熱量Q；

（8）確定換熱器的效能ε；

（9）判斷冷側入口溫度和出口溫度是否符合要求，即Tc,out-Tc,in＞ 0，符合則進行下一步計算，否則計算停止并報錯；

（10）判斷換熱器效能ε與效能假設值εas之差是否小于10-6，若滿足，進入下一步；若不滿足誤差要求，則令εas=ε，然后基于效能假設值更新出口溫度，重復步驟(3) ～ (10)直到滿足誤差要求；

（11）輸出計算結果，完成計算。

基于上述計算模型，可以通過理論計算的方式獲得PCHE 的傳熱特性，但現有的關聯式種類較多，無法直接判斷哪一種關聯式可以很好地適配本文所研究的問題，需要通過實驗數據來進行進一步的篩選，本文通過已有文獻中的實驗數據，來確定精確度、適用性相對較高的傳熱計算關聯式。

2.2 實驗數據庫構建

為了驗證關于PCHE 的10 種傳熱經驗關聯式以及6 種流動經驗關聯式的精確度與適用性，本文基于文獻[2,17,19-20,45,49-50]中與緊湊型換熱器相關的實驗數據構建實驗數據庫，一共選用了426 個數據點組成實驗數據庫（表4），其中包含逆流換熱、交錯流換熱兩種換熱方式，換熱工質有氣/氣換熱、氣/液換熱和液/液換熱三種組合。

表4 數據庫相關組成及研究內容Table 4 Composition of database and research content

基于文獻中所給出微通道的幾何條件與邊界條件，采用本文所建立的計算模型計算其相關工況下的傳熱特性獲得經驗解，然后將經驗解與實驗數據庫中的實驗解進行對比分析。在計算過程中，考慮到換熱過程中由于溫度和壓力對于換熱工質熱物性的影響，采用開源熱物性數據庫CoolProp 進行相關工質物性的獲取[51]來表征換熱工質在相應工況下的熱物性變化。

由于本文所采用相關數據的工況遠離臨界條件，預計工質不會出現熱物性的劇烈變化。此外，為評估計算模型計算結果的精度，使用均方根誤差（RMSE）來評估二者數據間的差異。RMSE 的計算公式為[2]：

式中，Dk和Yk分別表示理論值和實驗值；Ni表示可供比較的數據總數。

2.3 關聯式精度分析與篩選

基于本文構建的計算模型，可以獲得根據數據庫文獻所給出的緊湊型換熱器幾何特征與邊界條件計算得到換熱器的傳熱特性的經驗解，其中包括出口溫度Tout、效能ε、總傳熱系數UA（或總熱阻1/UA）以及總換熱量Q。下面進行計算模型所得經驗解與數據庫中的實驗解之間精度的對比分析。

圖6 顯示了換熱器出口溫度Tout的經驗解和數據庫中實驗解的對比。從圖6(a)、(b)中可以看出，所有的經驗關聯式對于出口溫度Tout的計算中絕大多數數據點的誤差維持在10%以內，不同的關聯式都表現出較好的計算精度。從整體來看，經驗解與實驗解的吻合度較好，尤其在高溫段時精度明顯較好，進一步對比低溫區段[圖6(c)、(d)]，可以看到盡管不同關聯式計算所得經驗解與實驗解之間仍然存在不同的差異，但整體預測精度依然較好。相 對 來 說，Stephan 等[18]、Muzychka 等[8]和Sarmiento等[20]所提出關聯式計算所得經驗解更接近文獻中的實驗值。

圖6 出口溫度Tout對比分析Fig.6 Comparative analysis of outlet temperature Tout

圖7 給出了計算模型對換熱器總傳熱系數UA（或總熱阻1/UA）與總效能ε的計算結果和實驗測試結果的對比。從圖7(a)中對于總傳熱系數UA（或總熱阻1/UA）的計算結果可以看到，Sarmiento 等[20]的關聯式在低Re時的計算結果與實驗結果吻合較好，計算結果的誤差基本維持在10%以內；但隨Re增大，即Re＞2300 后，計算結果明顯低于實驗測量值，誤差增大到15%甚至更大。這種現象表明，該關聯式在一定程度上低估了在高Re下對流換熱的傳熱系數。類似地，其他關聯式也表現出類似的情況。其中，Lee等[26]提出關聯式的計算結果相比實驗解明顯偏小，低估了對流傳熱系數；與之相反，Peng等[28]提出的關聯式的計算結果遠遠高估了對流傳熱系數UA的值，最大誤差超過了30%；Muzychka 等[8]關聯式在高Re下的計算結果的精度則更差一些。類似地，對比計算模型對于換熱器效能ε的計算結果的誤差，盡管Stephan 等[18]、Sarmiento 等[20]的關聯式都存在一定的計算偏差，不過僅有少數的數據點計算誤差達到30%，排除這些誤差較大的點之后，二者的計算精度相比其余關聯式的表現來看依然較好。其余幾種關聯式的表現大同小異，呈現出類似的誤差分布。進一步對比上述對出口溫度Tout經驗解與實驗解，從10種計算關聯式中篩選出精度相對 最 高 的3 種 關 聯 式，即Muzychka 等[8]、Stephan等[18]、Sarmiento等[20]的關聯式。

圖7 總傳熱系數UA（或總熱阻1/UA）與效能ε 對比分析Fig.7 Comparative analysis of total heat transfer coefficient UA (or entire thermal resistance 1/UA) and efficiency ε

圖8 給出了本計算模型對換熱器換熱量Q的計算結果與實驗測量值的對比。從圖8(a)中可以看出，10 種關聯式整體的均表現出相對較好的計算精度，雖然其中有個別數據點的誤差較大，不過絕大部分數據點的計算誤差保持在20%以內。選取數據庫中局部的一些數據點進行對比分析[圖8(b)～(d)]，可以看出一些關聯式在換熱量和較低時精度明顯變差，比如Lee 等[26]、Grigull 等[27]提出的關聯式。與之相反，Stephan等[18]、Sarmiento等[20]提出的關聯式盡管在一部分數據點上計算存在較大的誤差，但計算精度仍在可接受范圍，無論是整體數據點還是局部數據點都表現出較好的計算精度。

圖8 總換熱量Q對比分析Fig.8 Comparative analysis of the total heat exchange amount Q

經過上述對比分析，可以得出本文建立的計算模型可以較好地預測換熱器的傳熱特性，其中，Stephan 等[18]、Sarmiento 等[20]提出的關聯式很好地預測了不同工況下換熱器傳熱特性。前者雖然是基于圓形通道的實驗擬合而得到經驗公式，但對于非圓形通道的實驗預測依然保持很好的預測精度；后者盡管在較低Re時會過分估計對流換熱量，但隨著Re的增加，誤差逐漸減小，計算結果的精度也同步提高。因此，通過對比實驗測量值與計算模型的理論值，可以初步確定計算模型采用Stephan 等[18]或者Sarmiento 等[20]提出的計算關聯式來預測PCHE 的傳熱特性。

圖9 所示為模型計算結果與數值仿真結果的對比分析，可以看到除了Fried 等[44]提出的關聯式誤差較大，Shah 等[24]的計算誤差稍大，其余關聯式均和數值模擬結果吻合很好，誤差基本保持在10%以內。從圖中可以看到，預測模型與數值仿真結果不僅變化趨勢保持一致，數值也基本吻合，由此可以進一步說明前面針對實驗測量值偏大的分析是正確的。

圖9 換熱器Fanning摩擦系數f與數值仿真值對比分析Fig.9 Comparison and analysis of Fanning friction coefficient f of heat exchanger with numerical simulation

圖10 顯示了本文計算模型對換熱器Fanning 摩擦系數f的計算結果與實驗測量結果的對比，可以看到除了Fried 等[44]提出的關聯式預測結果很差，Red＜ 2000 時的計算結果遠遠偏離了其余三種關聯式的計算結果，而2000 ＜Red＜4000 時與其余關聯式又有所重合，但下降速度則較低。其余關聯式的預測曲線相比實驗測試結果的變化趨勢基本一致，但數值之間仍有一定差距，其中原因可能在于通道內部轉折彎角處出現二次流，增大了壓力損失，同時換熱器的加工過程中出現了部分通道的畸變導致換熱器內部出現局部的收縮和膨脹導致壓力損失，從而使得實驗測量值大于模型計算值。相比之下，Kim 等[14]、Muzychka 等[8]提出的關聯式所得到的經驗解與數值解更為接近，不僅誤差小而且趨勢吻合較好，二者對于實驗解的預測盡管有一定的誤差，但是整體的變化趨勢十分接近。因此，通過對比實驗測量值、數值仿真值與計算模型的理論值，可以初步確定計算模型采用Kim 等[14]或Muzychka 等[8]提出的計算關聯式來預測PCHE 的流動特性。

圖10 換熱器Fanning摩擦系數f與實驗測量值對比分析Fig.10 Comparison and analysis of Fanning friction coefficient f of heat exchanger with experimental measurement

綜上所述，通過對比實驗測量值與計算模型的理論值，可以初步確定本計算模型采用Stephan等[18]、Sarmiento 等[20]提出的關聯式來計算PCHE 的傳熱特性，采用Kim 等[14]、Muzychka 等[8]提出的關聯式計算流動特性。

3 仿真分析對比

3.1 邊界條件與仿真模型

前文依據本文建立的計算模型，分析了不同Nusselt 數的經驗關聯式對于微通道換熱器傳熱特性的預測情況，對比了經驗解與實驗解之間的精度分布，篩選出2 種計算關聯式對于層流區域內對流換熱計算精度和適用性表現較好。接下來，本文針對一種印刷板式微通道換熱器（PCHE）的換熱核芯進行仿真分析，該結構由冷側、熱側換熱板結構組成，如圖11所示?；谏鲜鲇嬎隳Ｐ鸵约昂Y選出的經驗關聯式對該換熱核芯進行傳熱特性的求解，對比分析經驗解與數值解的誤差。

圖11 換熱核芯結構示意圖Fig.11 Structure of heat transfer core

從圖11中可以看出，該換熱核芯為兩種換熱工質的對流換熱，在換熱板的中間部分為主要的換熱形式，即準逆流；在換熱板的兩端（即入口、出口）為交錯流。其中，兩種換熱工質分別為航空燃油RP-3以及航空潤滑油4050，在數值仿真與模型計算時，前者的熱物性采用Deng 等[52-54]和Cheng 等[55]的實驗測試結果來表征，后者的物性依據楊九高[56]給出的物性關系式來表征，兩種換熱工質物性均隨溫度變化，采用多項式擬合，擬合多項式如表5所示。

表5 航空燃油RP-3擬合多項式［52-55］Table 5 Fitting polynomial of aviation fuel RP-3［52-55］

航空潤滑油4050的物性關系式為：密度ρ

式中，ρ表示密度，kg/m3；cp表示定壓比熱容，J/(kg·K)；λ表示熱導率，W/(m·K)；υ表示運動黏度，mm2/s；T表示溫度，K。

換熱板片通常采用化學刻蝕的方法進行加工，加工得到的通道橫截面為接近半圓結構的橢圓形，在理論研究時會對該結構進行等效處理，本文采用Sarmiento 等[20]和Liu 等[22]提出的理論為基礎，將半圓形截面微通道等效為一種矩形結構，如圖12所示。

等效后的矩形翅片的平均厚度x′、y′可以根據截面積來進行估計，二者的表達式為：

式中，R表示近似圓形通道的半徑，m；AT1、AT2分別表示半圓形翅片的冷側截面積、熱側截面積，m2；其余參數見圖12，詳細推導過程見文獻[2]。

圖12 半圓形微通道截面示意圖[2]Fig.12 Semicircular microchannel cross section[2]

由于該PCHE 換熱核芯結構中換熱板片數量較多，整體數值模擬顯然是十分困難的，將會耗費大量的計算資源和時間?？梢园l現換熱核芯結構呈現出周期性分布，即核芯部分以一層冷側換熱板與一層熱側換熱板組成一個周期單元周期性分布。因此，抽取其中一個換熱單元作為仿真對象，PCHE換熱核芯的周期換熱單元以及換熱單元邊界類型如圖13所示。

圖13 換熱單元及邊界類型Fig.13 Heat transfer unit and boundary type

其中，換熱板片的結構參數見表6，換熱單元的邊界類型有：

表6 換熱板片幾何尺寸Table 6 Heat exchange plate geometry

（1）入口邊界 邊界類型為速度入口（表7），并且背壓為3 MPa；

（2）出口邊界 出口邊界類型為壓力出口，背壓為3 MPa；

（3）周期邊界 換熱單元上下兩個面；

（4）絕熱邊界 通道出口處固體材料的截面；

（5）恒溫邊界 通道入口處固體材料的截面，溫度等于來流溫度；

（6）對流邊界 外側，對流傳熱系數為8 W/(m2·K)，環境溫度為298.15 K。

對換熱單元模型進行網格無關性驗證分析，驗證結果如圖14所示。隨著模型網格數量的增加，換熱單元出口速度不斷下降逐漸趨向于一定值（即0.22 m/s），但網格數量超過2200 萬個后，模擬結果基本保持恒定，以此確定本模型數值模型時的網格數量為≥2200萬個。

圖14 網格無關性驗證Fig.14 Grid independence verification

確定換熱單元以及其邊界條件后，基于入口Reynolds 數均低于2300，求解模型選擇層流模型（Laminar）。為保證求解精度，采用雙精度求解計算，求解算法采用Coupled 算法的二階迎風格式進行求解，殘差設置為1×10-10，迭代步數設置在2000步左右。如果模擬時出現回流現象，則先使用Simple 算法的一階迎風格式進行求解，在計算基本收斂且回流現象消失，即為高階求解建立初場，然后再使用Coupled 算法的高階形式進行求解來保證求解精度。

3.2 數值解與經驗解精度分析

基于建立的PCHE 周期換熱單元，針對相應的工況開展數值仿真和解析計算的結果對比分析研究。工況安排見表7，其中熱側換熱工質航空潤滑油4050的入口溫度和速度是一個恒定值，相對發生變化的只有冷側換熱工質航空燃油RP-3，其中入口流量的范圍為1.0～5.0 m3/h，入口溫度的范圍40～80℃），合計45 組種工況開展相應的數值仿真和經驗公式計算。

表7 工況分布Table 7 Working condition

在進行解析計算時則依賴于本文所建立的計算模型來預測換熱單元的傳熱特性，本文建立的計算模型最終確定采用Stephan 等[18]、Sarmiento 等[20]提出的關聯式來計算PCHE 的傳熱特性，在進行解析計算時將同時采用上述候選關聯式進行計算，并對比不同關聯式的計算模型的經驗解與數值仿真結果的數值解之間的誤差分布。

圖15給出了換熱單元冷、熱側換熱工質出口溫度的數值解與兩種關聯式的經驗解之間的對比結果，從圖中可以看出，兩種關聯式對換熱工質出口溫度的預測均存在一定的偏差，不過偏差基本都保持在10%以內，進一步證明本文研究結論的正確性，側面說明兩種關聯式對于換熱器的傳熱特性預測具有一定的指導意義。對比而言，雖然Stephan等[18]提出的關聯式僅適用于層流，但其無論是對于冷側工質，還是熱側工質的求解結果的連續性相對更好；相反，Sarmiento等[20]提出的關聯式同時覆蓋層流、過渡區以及紊流，雖然其適用范圍更廣，但在本文所研究內容中其適用性則相對較差，可以看出其求解結果不穩定，對于熱側工質出口溫度的預測結果尤為明顯。因此，對于本文所研究的相關內容，Stephan 等[18]提出的關聯式顯然更具實用性。進一步分析，由Stephan 等[18]提出的關聯式對于冷、熱側換熱工質的求解結果與數值仿真結果的對比可以看到經驗解相較于數值解明顯偏低[圖15(a)]，經驗解相較于數值解明顯偏高[圖15(b)]。初步分析，其原因在于數值仿真計算時，默認求解過程是理想的，對于很多不可控因素的影響進行了忽略，比如環境溫度、加工精度等，然而Stephan 等[18]提出的關聯式是基于實驗測試結果進行擬合而獲得，所以使得計算關聯式的求解結果相較于數值仿真的理想結果偏低。

圖15 出口溫度的數值解與經驗解對比Fig.15 Numerical solution of outlet temperature compared with empirical solution

在確定傳熱特性的計算關聯式之后，在此基礎上對比分析采用Kim 等[14]、Muzychka等[8]提出的關聯式計算流動特性的經驗解與數值仿真結果之間的偏差，圖16 給出了換熱器壓損與Fanning 摩擦系數的經驗解與數值解的對比結果，由于本節中僅有冷側換熱工質的溫度和流速發生變化，而熱側換熱工質保持恒定，因此只對冷側換熱工質的流動特性進行了對比分析。從圖16(a)中可以看到，兩種關聯式對于換熱器冷側壓損pcost的計算結果與數值仿真結果吻合度表現很好，二者誤差均在10%以內，只有在整體壓損較?。ɑ蛉肟诹魉佥^?。r存在相對較大的偏差，不過整體而言二者的吻合度、穩定性的表現均很好。進一步比較可以發現，Kim 等[14]提出的關聯式與仿真結果的吻合度相較于Muzychka 等[8]的模型更高一點，不過二者的誤差量級基本一致。圖16(b)中給出了冷側Fanning 摩擦系數的經驗解與數值解之間的對比分析，其結果表現與壓損pcost的表現基本一致，整體誤差均維持在10%以內，并且隨著流速的增加兩種關聯式的預測結果與數值結果越接近，相對而言Kim 等[14]提出關聯式的表現略勝一籌。

圖16 pcost 和 f數值解與經驗解對比Fig.16 Comparison of numerical and empirical solutions of pcost and f

綜上所述，通過對比采用Stephan 等[18]、Sarmiento等[20]提出的關聯式對于傳熱特性的經驗解和采用Kim 等[14]、Muzychka等[8]提出的關聯式對于流動特性的經驗解與數值仿真結果所得的數值解之間的精度，可以得出，Stephan 等[18]與Kim 等[14]提出的關聯式對于本文換熱單元的傳熱與流動特性的預測精度相對最高。

4 結 論

本文從理論計算的角度出發，基于現有的研究成果建立了PCHE 的經驗解的計算模型，對比分析現有研究中提出的傳熱與流動特性的計算關聯式，篩選出精度相對最高的關聯式，并針對一種PCHE結構流動與傳熱特性進行了數值仿真與計算模型求解，對比分析了數值解與經驗解，得出以下結論。

（1）基于現有研究成果，以PCHE 整體流動與換熱特性的理論計算出發選取目前研究文獻中所給出的關于層流Nusselt 數以及Fanning 摩擦系數的計算關聯式，建立了較為全面的用于預測PCHE 流動與換熱特性的計算模型。

（2）基于現有文獻中關于層流的實驗研究構建了實驗數據庫，對比分析具有代表性的關聯式的經驗解與實驗解的精度，并確定采用Stephan 等[18]、Sarmiento 等[20]提出的關聯式來計算PCHE 的傳熱特性，采用Kim 等[14]、Muzychka 等[8]提出的關聯式計算流動特性。

（3）基于一種PCHE 核芯結構，選取換熱單元作為研究對象進行數值仿真和計算模型計算，然后對比分析了流動與傳熱特性的數值解與經驗解之間的誤差分布，絕大多數點的誤差在10%以內，驗證了本文提出計算模型在理論求解中的準確性與可信度。此外，Stephan 等[18]與Kim 等[14]提出的關聯式對于本文換熱單元的傳熱與流動特性的預測精度相對最高。