2. 墨尔本皇家理工大学 航空与机械工程学院, 维多利亚 3083, 澳大利亚
2. School of Aerospace, Mechanical and Manufacturing Engineering, RMIT University, Victoria 3083, Australia
核态沸腾凭借其极高的传热效率在工业领域有着非常广泛的应用，尤其在微电子领域的芯片冷却及核反应堆的熔融物容器冷却系统等方面有着广阔的应用前景。但是，临界热通量(critical heat flux,CHF)的存在却大大限制了核态沸腾的传热强度。多年来，如何提高CHF一直是传热学界的研究热点之一。强化沸腾传热的方法有很多，其中向基液中加入纳米颗粒制成纳米流体以其高效、简便等特点成为最有前景的强化沸腾换热的方法之一。大量实验表明，纳米流体的CHF相对于其基液的提高幅度可高达200%，这为核态沸腾强化换热提供了十分广阔的前景。但是，作为一种全新的换热介质，纳米流体的核态沸腾具有很多不同于纯工质的特性，其中最为显著的就是纳米颗粒沉积及其引起一系列沸腾特征的变化[1]。目前，人们对其内在机理还缺乏足够的认识，相应的理论预测模型尚无建立，这也大大阻碍了纳米流体作为换热介质在工业系统中的应用。
事实上，纳米流体核态沸腾也具有纯工质核态沸腾的一般规律，已在沸腾模拟中得到广泛应用的壁面热通量拆分(heat flux partitioning，HFP)模型[2]理论上同样适用于纳米流体，只是由于其特有的纳米颗粒沉积现象[1]对汽泡的成核、生长及脱离特征造成了影响，因此需要对HFP模型中的有关汽泡动力学参数进行重新定义。
因此，本文基于现有实验数据对纳米流体核态沸腾的相关物理参数进行了重新定义，并以之修正HFP模型。数值计算结果与实验数据对比表明，修正后的HFP模型对于纳米流体具有较好的预测精度。
1 壁面热通量拆分模型 1.1 模型描述根据HFP模型的基本思路(见文[2])，核态沸腾壁面上的热量通过对流、蒸发和激冷3种方式传入流体，即总热量可按如下方式拆分：
| $\begin{array}{l} q = \frac{{\rm{\pi }}}{6}d_{\rm{b}}^3w{\rho _{\rm{g}}}f{N_a}{h_{{\rm{fg}}}} + \\ \frac{2}{{\sqrt {\rm{\pi }} }}f{A_{\rm{q}}}\sqrt {{t_w}{\lambda _{\rm{l}}}{\rho _{\rm{l}}}{c_{{\rm{p,l}}}}} \left( {{T_{\rm{W}}} - {T_{\rm{l}}}} \right) + \\ {A_{\rm{c}}}St{\rho _{\rm{l}}}{c_{{\rm{p,l}}}}{u_{\rm{l}}}\left( {{T_{\rm{W}}} - {T_{\rm{l}}}} \right). \end{array}$ | (1) |
由于纳米包覆层的存在，纳米流体中加热壁面上的沸腾汽泡表现出不同于纯工质汽泡的脱离特征。实验观测发现，纳米流体核态沸腾的汽泡脱离直径不仅比纯工质大，而且还会随着壁面过热度的增大而变化[3]。但由于问题的复杂性，即使对于纯工质的核态沸腾，目前也还没有一个公式或理论模型能够将这些因素系统地加以考虑，汽泡脱离直径的计算还主要依赖经验方法。对于纳米流体这一新型的换热介质，目前尚未有可靠的公式来计算其核态沸腾汽泡脱离直径。
在机理理解不足的情况下，为了增加数学模型结果的准确性并同时反映纳米流体核态沸腾的现象特征，本文对Gerardi等[3]的实验数据进行了拟合，将汽泡脱离直径dbW拟合为壁面过热度ΔTsup的函数，
| $\begin{array}{l} {d_{{\rm{bW}}}} = - 1.91 \times {10^{ - 3}} + 4.21125 \times {10^{ - 4}}\Delta {T_{\sup }} - \\ 1.70945 \times {10^{ - 5}}\Delta T_{\sup }^2 + 2.03938 \times {10^{ - 7}}\Delta T_{\sup }^3. \end{array}$ | (2) |
出于分析和对比的目的，本文还分别采用了Tolubinsky-Kostanchuk(见文[2])以及Phan等[4]提出的汽泡脱离直径公式(式(3)及(4))进行了计算。其中：Tolubinsky和Kostanchuk将汽泡脱离直径整理为液体过冷度ΔTsub的函数； Phan等的公式是在经典的Fritz公式的基础上发展而来，加入了液体接触角θ的影响。
| ${d_{{\rm{bW}}}} = \min \left[{0.0006\exp \left( {\frac{{\Delta {T_{\sup }}}}{{45}}} \right),0.00014} \right].$ | (3) |
| ${d_{{\rm{bW}}}} = {C_{{\rm{bW}}}}\frac{{2 + 3\cos \theta - {{\cos }^3}\theta }}{4}\sqrt {\frac{\sigma }{{g\left( {{\rho _{\rm{l}}} - {\rho _{\rm{g}}}} \right)}}} .$ | (4) |
研究发现，纳米包覆层是促使活化核心密度发生变化的主要因素[5]。多孔型纳米包覆层不仅会显著地改变壁面的表面性质，还会改变加热壁面的表面形态。这是活化核心密度变化的两个主要诱因。
为了突出考虑以上两个因素对于纳米流体核态沸腾活化核心密度的影响，Ganapathy等[6]在 Benjamin-Balakrishnan公式[7]的基础上纳入液体接触角θ及壁面粗糙度变化对活化核心密度的影响,
| $\begin{array}{l} {N_{{\rm{a,nf}}}} = 218.8\frac{1}{\gamma }P{r^3}\Delta T_{\sup }^3\left[{\left( {14.5 - 4.5\left( {\frac{{{R_{\rm{a}}}p}}{\sigma }} \right) + } \right.} \right.\\ 0.4\left. {{{\left( {\frac{{{R_{\rm{a}}}p}}{\sigma }} \right)}^2}} \right){\left( {\frac{{1 - \cos \theta }}{{1 - \cos {\theta ^ * }}}} \right)^{ - 3}}{\left. {{{\left( {\frac{{{R_{\rm{a}}}}}{{{D_{\rm{p}}}}}} \right)}^{ - 0.5}}} \right]^{ - 0.4}}. \end{array}$ | (5) |
同样，出于进一步比较和分析的目的，本文还分别将Benjamin-Balakrishnan[7]活化核心密度公式(式(6))和Wang-Dhir[8]公式(式(7))纳入模型进行了对比计算:
| $\begin{array}{l} {N_a} = 218.8\frac{1}{\gamma }P{r^3}\Delta T_{\sup }^3\\ \left[{\left( {14.5 - 4.5\left( {\frac{{{R_{\rm{a}}}p}}{\sigma }} \right) + } \right.} \right.0.4{\left. {{{\left( {\frac{{{R_{\rm{a}}}p}}{\sigma }} \right)}^2}} \right]^{ - 0.4}}. \end{array}$ | (6) |
| ${N_a} = 7.81 \times {10^{ - 29}}{C_n}\left( {1 - \cos \theta } \right){\left( {\frac{{2\sigma {T_{sat}}}}{{{\rho _v}{h_{fg}}\Delta {T_{\sup }}}}} \right)^{ - 6}}.$ | (7) |
本文以Gerardi等[3]的水基氧化硅纳米流体池沸腾实验为基础进行方程的求解和分析。数值计算的原型为一直径200 mm、高200 mm的圆柱状容器。直径为50 mm的加热表面位于容器底部中心处。由于两相流场具有轴对称性，本文将其简化为100 mm×200 mm的2维计算域。计算域采用矩形结构化网格划分，网格测试表明当网格密度到达 75×50时，计算结果具有网格无关性。加热表面上的热流密度为50~400 kW/m2。容器的其余表面根据实验条件假定为绝热表面。通过将容器顶部定义为气相的质量汇项可近似模拟液体表面的气体排放过程。
容器内流体的流动过程采用双流体模型(见文[2])进行模拟。本文建立起来的HFP模型采用迭代二分法求解后，作为边界条件纳入到双流体模型中。最后，通过商用计算流体力学软件CFX-4进行求解。当液相的残差小于1×10－4时，判定计算收敛。
2.2 模拟结果及讨论 2.2.1 汽泡脱离直径本文作者所在课题组前期[9]分别对式(3)及(4)进行了分析。根据Kim等[5]的实验数据，取纯水静态接触角为79°，纳米流体静态接触角为22°。将式(3)及式(4)纳入模型后其计算结果与Gerardi等[3]的实验数据对比如图1所示。
图1表明，虽然式(3)的计算结果在较低壁面过热度时与实验数据吻合较好，但这一公式并没有反映出纯水或纳米流体沸腾汽泡脱离直径随着热流密度增大的变化规律。式(4)反映出汽泡脱离直径随液体接触角减小而增大，这与大多数实验观察在规律上一致，但却没能作到定量预测。显然，式(3)和式(4)的计算结果均与实验数据均存在较大偏差，在纳米流体核态沸腾系统中，都不能作为计算汽泡脱离直径的封闭方程，因此本文根据Gerardi等的实验数据将汽泡脱离直径拟合为壁面过热度的函数。对于纳米流体核态沸腾系统，本文拟合的汽泡脱离直径计算式(2)与Gerardi的实验数据对比如图2所示。
|
| 图 2 汽泡直径公式(式(2))计算结果与实验数据(文[2])对比 |
图2表明，出当壁面过热度在10~35 K时，式(2)与实验数据有较好的吻合度。但是，同样也要指出，虽然式(2)与Gerardi等[3]的实验数据吻合较好，但该公式只考虑了壁面过热度对汽泡脱离直径的影响，并没有考虑对汽泡脱离直径事实上存在影响的诸多其他因素，包括液体静态接触角、流动速度及浮力等。因此，式(2)也有其一定的局限性。如何从机理角度计算汽泡脱离直径仍然是未来的研究重点之一。
2.2.2 活化核心密度影响活化核心密度的因素很多，Wang-Dhir公式(式(7))和Benjamin-Balakrishnan公式(式(6))分别考虑了液体接触角和壁面材料及粗糙度的影响，分别将它们整合入HFP模型后，其计算结果与实验数据的对比如图3所示。
|
| 图 3 各活化核心密度计算式与实验数据的对比 |
尽管式(6)和式(7)的精度均不能令人满意，但相比之下，式(6)比式(7)与实验数据吻合更好。正因为式(6)考虑了更多的影响因素，因此被认为更有可能从机理上描述汽泡成核。式(5)正是在式(6)基础上加入更多影响因素后提出的。式(5)不仅通过润湿性能参数(1-cosθ)/(1-cosθ*)来反映液体接触角对活化核心密度的影响，还通过Ra/Dp考虑了表面粗糙度及纳米颗粒尺寸的影响。将式(5)作为封闭方程纳入HFP模型后其计算结果与文[6]的对比如图4和5所示。
|
| 图 4 静态接触角对活化核心密度预测的影响 |
图4进一步探究了壁面静态接触角的变化对纳米流体活化核心密度的影响。计算条件为加热壁面平均粗糙度200 nm，纳米颗粒直径为50 nm，液体接触角θ依次为22°、50°、60°、70°。图4表明当静态接触角为22°时，模型计算结果比Ganapathy等的预测结果偏小； 当静态接触角在50°~70°范围内时，计算结果与Ganapathy等的预测结果较为符合。图4还表明，对于给定的静态接触角，活化核心密度随着壁面过热度的增大而增加； 对于给定的壁面过热度，活化核心密度随着接触角的降低而减小。这与Kim等[5]及Geradi等[3]的实验观测一致。
图5展示了表面粗糙度及纳米颗粒尺寸相对大小对活化核心密度的影响。计算条件为：壁面粗糙度200 nm，静态接触角60°，纳米颗粒的直径依次为20、50、100、200、300、400、500 nm。由图5可以看出，在给定的壁面过热度及粗糙度条件下，活化核心密度随着纳米颗粒直径的增加而减小。但由于实验数据的缺乏，活化核心密度与表面粗糙度及纳米颗粒尺寸的关系未来还需要更多的实验来验证。
|
| 图 5 纳米颗粒直径对活化核心密度预测的影响 |
可见，Ganapathy等的公式(式(5))虽然能从趋势上反映出液体静态接触角以及纳米颗粒尺寸对活化核心密度的影响，却没能够作到定量预测。 未来还需要对Ganapathy公式进一步改进和完善，同时对活化核心密度与纳米颗粒尺寸的关系进行进一步验证。
2.2.3 传热系数由于纳米流体核态沸腾机理研究的缺乏，当前还不能准确地判断纳米流体传热系数(heat transfer coefficient，HTC)的变化趋势。Das等发现纳米流体的沸腾传热系数和其基液相比会降低，然而这一结论却与Wen和Ding等的实验结果恰好相反。Kole和Dey的实验数据甚至表明纳米流体的传热系数会随纳米颗粒含量的增加先升高再降低(见文[10])。尽管实验观测到的纳米流体核态沸腾传热系数变化差别很大，但目前普遍认为这些变化仍然是由于纳米包覆层导致的，因此仍然可以从壁面形态及润湿性的角度着手分析。为了研究液体接触角以及表面粗糙度和纳米颗粒尺寸对纳米流体沸腾传热系数的影响，本文对前述HFP模型计算结果进行了分析，并与Ganapathy等的实验数据[6]进行了对比，如图6和7所示。
|
| 图 6 静态接触角对传热系数的影响 |
|
| 图 7 纳米颗粒直径对传热系数的影响 |
图6表明，在相同壁面过热度下，预测的沸腾传热系数随着液体接触角的降低而降低。这与大多数实验测试结果一致，因为沸腾传热系数的大小主要受活化核心密度的影响，当减小的接触角导致活化核心密度减小时，纳米流体的传热性能就会随之降低。
图7表明，在相同壁面过热度下，沸腾传热系数随着纳米颗粒直径的增加而减小。事实上，纳米颗粒尺寸对纳米流体沸腾传热性能的影响规律十分复杂。本文对Das等[11]的研究成果与Ganapathy等的实验数据[6]进行了汇总分析，如图8所示。图8中参数ψ的定义为平均表面粗糙度与纳米颗粒尺寸的相对大小，ψ=Ra/Dp。
从图8可以看出，当ψ位于1.0左右时，纳米流体核态沸腾传热弱化，反之当ψ远离1.0时，纳米流体核态沸腾传热均得到强化。Narayan等[12]认为，当纳米颗粒尺寸与壁面微结构尺寸大小相当时，纳米颗粒能够沉积到这些微结构中，导致原本存在的活化核心消失，导致沸腾传热性能下降； 反之，当两者尺寸相差较大时，沉积到壁面上的纳米颗粒能够创造出更多的活化核心，从而强化沸腾传热。但是，式(5)中的(Ra/Dp)-0.5项并没有反映出这一规律。这也是今后研究所需要解决的问题之一。
|
| 图 8 纳米流体传热系数的提升与表面-颗粒相互影响参数的实验数据汇总图 |
纳米颗粒在加热壁面上沉积形成纳米包覆层进而导致壁面的表面性质及形态发生改变是纳米流体核态沸腾不同于纯工质沸腾的特有现象。多孔型的纳米包覆层不仅会提高壁面的润湿性能，还会显著地改变加热壁面的表面形态。目前普遍认为，这正是纳米流体不同于纯工质沸腾传热特征的根本原因。当基于经典HFP模型建立纳米流体核态沸腾传热预测模型时，以上现象也是需要重点考虑的因素。本文基于现有实验数据建立了纳米流体核态沸腾的汽泡动力学参数计算公式，以之修正已有的HFP模型并根据文献实验条件进行了数值计算。本文研究表明:
1) 对加热壁面润湿性能的提高以及壁面粗糙度变化的准确描述是建立纳米流体适用的HFP模型的关键所在。以往研究者提出的计算公式虽然能够在趋势上反映纳米流体核态沸腾汽泡动力学参数的变化特征，但却没有能够实现定量预测。
2) 鉴于纳米包覆层对于汽泡成核、生长及脱离规律的影响，在未来研究中，应当将液体接触角及加热壁面粗糙度的表征作为建立纳米流体核态沸腾传热理论模型的重点。
| [1] | 李祥东, 屠基元. 纳米流体核态沸腾的机理探讨及数值模拟 [J]. 工程热物理学报, 2013, 34(6): 1096-1100.LI Xiangdong, TU Jiyuan. Mechanistic modelling of nucleate boiling of nanofluids [J]. Journal of Engineering Thermophysics, 2013, 34(6): 1096-1100. (in Chinese) |
| [2] | Yeoh G, Tu J. Modelling Subcooled Boiling Flows [M]. New York: Nova Science Publishers, 2009: 5-36. |
| [3] | Gerardi C, Buongiorno J, Hu L, et al. Infrared thermometry study of nanofluid pool boiling phenomena [J]. Nanoscale Research Letters, 2011, 6(1): 1-17. |
| [4] | Phan H, Caney N, Marty P, et al. Surface wettability control by nanocoating: The effects on pool boiling heat transfer and nucleation mechanism [J]. International Journal of Heat and Mass Transfer, 2009, 52(23): 5459-5471. |
| [5] | Kim S, Bang I, Buongiorno J, et al. Surface wettability change during pool boiling of nanofluids and its effect on critical heat flux [J]. International Journal of Heat and Mass Transfer, 2007, 50(19/20): 4105-4116. |
| [6] | Ganapathy H, Sajith V. Semi-analytical model for pool boiling of nanofluids [J]. International Journal of Heat and Mass Transfer, 2013, 57(1): 32-47. |
| [7] | Benjamin R, Balakrishnan A. Nucleation site density in pool boiling of saturated pure liquids: Effect of surface microroughness and surface and liquid physical properties [J]. Experimental Thermal and Fluid Science, 1997, 15(1): 32-42. |
| [8] | Wang C, Dhir V. Effect of surface wettability on active nucleation site density during pool boiling of water on a vertical surface [J]. Journal of Heat Transfer, 1993, 115(3): 659-669. |
| [9] | Li X, Li K, Tu J, et al. On two-fluid modeling of nucleate boiling of dilute nanofluids [J]. International Journal of Heat and Mass Transfer, 2014, 69: 443-450. |
| [10] | Wu J, Zhao J. A review of nanofluid heat transfer and critical heat flux enhancement: Research gap to engineering application [J]. Progress in Nuclear Energy, 2013, 66: 13-24. |
| [11] | Das S, Narayan G, Baby A. Survey on nucleate pool boiling of nanofluids: The effect of particle size relative to roughness [J]. Journal of Nanoparticle Research, 2008, 10(7): 1099-1108. |
| [12] | Narayan G, Baby A, Das S. Mechanism of enhancement/deterioration of boiling heat transfer using stable nanoparticle suspensions over vertical tubes [J]. Journal of Applied Physics, 2007, 102(7): 074317-1-074317-7. |
| [13] | Wen D, Ding Y. Experimental investigation into the pool boiling heat transfer of aqueous based alumina nanofluids [J]. Journal of Nanoparticle Research, 2005, 7(2/3): 265-274. |
| [14] | Das S, Putra N, Roetzel W. Pool boiling nano-fluids on horizontal narrow tubes [J]. International Journal of Multiphase Flow, 2003, 29: 1237-1247. |
| [15] | Bang I, Chang S. Boiling heat transfer performance and phenomena of Al2O3-water nano-fluids from a plain surface in a pool [J]. International Journal of Heat and Mass Transfer, 2005, 48: 2407-2419. |