池沸腾加热器换热优化设计——基于Boltzmann方法研究并选取加热器最优参数

Optimized Heat Exchange Design for Pool Boiling Heaters—Studying and Selecting Optimal Parameters for Heaters Based on the Boltzmann Method

摘要: 本文采用了相变格子Boltzmann模型研究加热器高度、长度、表面润湿性和表面结构对池沸腾换热性能的影响。我们发现加热器高度较低时汽液相变速度更大，因此加热器高度越低热流密度越大。另一方面，加热器长度越短，换热性能越好；其次，在低过热度下全疏水加热器换热性能最好，而在高过热度下全亲水换热器性能则更好；最后，我们还研究了加热器表面微结构的影响，发现锥形表面加热器的换热性能要优于光滑表面加热器。

Abstract: This study uses the phase-change lattice Boltzmann model to investigate the influence of heater height, length, surface wettability and surface structure on the pool boiling heat transfer performance. We found that when the heater height is lower, the vapor-liquid phase change rate is higher, leading to a greaster heat flux density. Additionally, shorter heater lengths exhibit better heat transfer performance. Moreover, at low degrees of superheat, fully hydrophobic heaters demonstrate the best heat transfer performance, whereas at high degrees of superheat, fully hydrophilic heaters perform better. Lastly, we examined the impact of heater surface microstructures and found that the conical surface heater exhibits better heat transfer performance compared to the smooth surface heater.

1. 引言

随着科学技术水平的不断发展，高温高热流设备在热能动力、机械制造、石油化工、航天航空等众多领域的应用也愈加广泛。而在工业中最常使用的沸腾方式之一便是大容器饱和沸腾。在大容器饱和沸腾中，随着大量气泡的产生、增长和脱离会带走大量气化潜热并扰动池中的流体从而强化换热。因此大容器饱和沸腾有着换热能力强，换热均匀等优点。然而，沸腾换热性能会受到临界热流密度的限制，一旦热流密度超过临界热流密度，将会导致设备烧毁或者换热性能剧烈下降。因此，众多学者致力于研究影响临界热流密度的因素。而换热器作为直接影响换热过程的设备，其特性对沸腾换热性能有着直接的影响[1]-[3]。

一些学者发现换热器尺寸是影响沸腾换热器性能的重要因素。Mody [4]等实验研究了两种矩形加热器尺寸(10 × 10 mm和34.5 × 32 mm)下三种电介质液体(PP1、PP1C、FC-87)的沸腾换热过程，发现加热器尺寸的增加使得PP1、PP1C、FC-87沸腾的临界热流密度(critical heat flux, CHF)分别提高了53%、92%、103%。为了进一步揭示CHF与加热器尺寸之间的关系，Lee [5]等人实验研究了长度变化范围为3~16 mm，宽度变化范围为15~40 mm的13种不同矩形加热器情况下水和电子氟化液的沸腾过程，他们发现了CHF值随着加热器宽度增加而降低，并发现了影响CHF大小的关键因素为沸腾过程不同尺寸加热器导致的干斑大小和流体再润湿速度不同。Lee [6]等人使用宽度为10、30和50 mm、长度为10、30和100 mm的不锈钢加热器评估加热器尺寸对CHF的影响，发现与加热器宽度相比，加热器长度对CHF的影响更显著，且随着加热器长度的增加，CHF趋于减少。当加热器宽度超过30 mm后，宽度的增加对于CHF的影响不显著。当加热器长度超过100 mm后，长度的增加对于CHF的影响也不显著。Kam [7]等研究了长方体不锈钢平板加热器的尺寸和平板与水平表面的夹角对CHF的影响，研究得出对于一定角度范围，在100至200 mm范围内长度对CHF值无显著影响。而在40至50 mm范围内随着宽度增加，CHF值的总体趋势趋于下降，并随着角度的增加趋于收敛。

除了矩形加热器，Pattanayak [8]等研究了圆柱形加热器管径和管长对沸腾的影响，发现管径较小时形成汽穴的尺寸也较小，使得CHF幅值增加。另一方面，他们还发现随着加热管的长度增加，成核气泡的尺寸也会增加，导致在加热器表面频繁发生气泡合并，从而使得较长的加热器更易烧毁。Sang [9]等实验对比了涂有Al₂O₃纳米颗粒的加热器表面和光滑方形加热器表面的沸腾过程，发现无论是否有涂层，随着加热器尺寸的增加，都会导致临界热流密度下降。他们把临界热流密度减小的原因归结为加热器尺寸增加引起供给冷却液体的电阻路径变长。

上述研究均是在单个加热器基础上进行的，与单加热器情况相比，采取多加热器的方式可以控制气泡成核位点、气泡之间的相互作用以及加热的均匀性，因此不同数量的加热器对沸腾过程带来的影响也不尽相同[10]。Baltis等[11]研究了垂直排列的多个加热器的成核点位在流动沸腾中受到的影响，发现上游加热器的额外对流会使下游加热器的气泡成核频率增加，气泡分离直径增大，但该影响会随着加热器间距和下游加热器初始气泡成核频率的增加而减小。此外，上游加热器的成核气泡对下游加热器处的气泡成核具有抑制作用，这种抑制作用会随着下游气泡成核频率的增加而减弱。Balch等[12]指出采用多加热器可以提高成核点数量，并研究了五个加热器的情况下FC-72的沸腾换热过程，发现了与单个加热器的情况相比，采用多加热器的方式可以延迟膜态沸腾的开始。Sun [13]等人研究了垂直管两侧各存在多个加热器时的沸腾过程，发现管道下方产生的气泡在脱离加热器后会消失在周围流体中，管道中部产生的气泡脱离加热器后倾向于与周围气泡合并，而管道上方产生的气泡在脱离加热器后会以孤立气泡的状态继续运动。Yadav等[14]研究了流动沸腾过程中单个加热器、水平方向上均分分布的两个加热器以及三个加热器情况下空泡份额的分布，观察到低入口过冷度条件下的局部空隙率高于高入口过冷度条件下的局部空泡份额，并且随着加热器数量的增加，最大空泡份额也相应增加。

为了进一步提高沸腾换热性能，一些学者建议采用粗糙表面，如在加热器表面增加粗糙度或者翅片，并探索了表面粗糙结构对沸腾换热性能的影响。Eldesouki等人[15]实验研究不同圆形水平粗糙钢加热面上的沸腾过程和纳米颗粒粗糙表面的沸腾过程，重点分析了三种表面粗糙度(0.14微米，1.35微米，2.87微米)对沸腾换热性能的影响，发现粗糙表面可以提高124%换热性能。Mario等人[16]研究了粗糙铜管表面的沸腾换热过程，发现与光滑表面相比，粗糙表面的沸腾换热系数可以提高1.5倍。Cong-Yao等人[17]研究了六种粗糙表面上R134a的沸腾过程，发现当表面粗糙度从0.133增加到5.742时，其热通量可以从40 kW/m²增加到280 kW/m²。在基础加热表面增加翅片同样可以改变表面微结构并增加换热面积，进而强化换热。Jaswal等人[18]提出在加热器表面增加矩形和梯形翅片提高沸腾换热性能，并发现了矩形和梯形翅片表面的热通量值分别提高了52.3%和101.5%。Ma等人[19]以FC-72为工质，实验研究了均质和分形阵列微针翅片表面的池沸腾传热特性，得出与均匀的微针翅片表面相比，传热面积比较小的阵列微针翅片表面可以达到相似甚至更高的CHF。Mahyar Ghazvini等人[20]则在池沸腾设施中，对规则翅片和改良空心翅片两种翅片结构的传热性能进行了实验研究，得出改进翅片的传热系数比普通翅片高约20%。

前人的研究工作表明加热器尺寸、数量以及加热器表面形态对沸腾换热性能有非常大的影响。然而，由于沸腾过程的复杂性，目前对于该问题的研究还不充分。如加热器的尺寸与传热性能之间有什么规律；加热器的个数越多传热性能是否越好；润湿性和加热面的改变对传热性能具体有着怎样的影响等。基于此，本文将通过Gong和Cheng的改进格子Boltzmann模型[21]在保持边界条件不变的情况下，依次改变加热器的高度，长度，个数，润湿性和加热器形状进一步研究加热器形态对沸腾换热性能的影响，以便于人们更清晰的认识沸腾换热机理，并进一步开发传热性能高的加热器形态。

2. 池沸腾传热模型设计

2.1. 物理问题和边界条件

本课题研究大容器饱和沸腾过程，定义特征长度l₀、特征速度u₀和特征时间t₀如下：

$l_0=\sqrt{\frac{\sigma }{g\left(p_l-p_v\right)}}{u}_0=\sqrt{g{l}_0};t_0=\frac{l_0}{t_0}$

其中l₀为特征长度，σ为表面张力，ρ_l为液相密度，ρ_υ为气相密度。

如图1所示，物理计算区域是一个二维矩形区域，宽度为L_x，长度为L_y。底部红色区域为加热器，长为L_H，高为H。在数值模拟中，通过改变加热器的高度H、长度L_H和数量来进行优化设计。

Figure 1. Schematic diagram of the compute area

图1. 计算区域示意图

加热面底部采用定温加热，底部温度T_W视为常数。加热器两端面为绝热边界。计算区域上边界为对流边界，底部为绝热边界，计算域两侧边界为周期边界。计算域内部填充饱和水。

2.2. 数学模型

本文采用Gong和Cheng提出的格子Boltzmann两相流相变传热模型[22]来进行模拟实验。在玻尔兹曼模型中，密度分布函数和温度分布函数用来模拟相变过程。其中密度分布的演化方程由下式推出：

$f_i\left(x+c_i{\delta }_t,t+{\delta }_t\right)-f_i\left(x,t\right)=-\frac{1}{\tau }\left[f_i\left(x,t\right)-f_i{}^{\left(eq\right)}\left(x,t\right)\right]+\Delta f_i\left(x,t\right)$

其中$i=0,1,2,\cdots$ ；$f_i\left(x,t\right)$ 为位置x和时间t处速度为c_i的密度分布函数。τ为弛豫时间。$f_i{}^{eq}\left(x,t\right)$ 是平衡态分布函数，具体形式为：

$f_i{}^{eq}={\omega }_i\rho \left[1+\frac{c_i+u}{c_s^2}+\frac{{\left(c_i\cdot u\right)}^2}{2c_s^4}-\frac{u^2}{2c_s^2}\right]$ ,

${\omega }_i$ 为加权系数，$c_s=c/\sqrt{3}$ 为模型常数，$c={\delta }_x/{\delta }_t$ 。${\delta }_x$ 是网格步长，${\delta }_t$ 是时间步长。

式(1)中的$\Delta f_i\left(x,t\right)$ 为力项，其具体关系式可表达为：

$\Delta f_i\left(x,t\right)=f_i{}^{\left(eq\right)}\left[\rho \left(x,t\right),u+\Delta u\right]-f_i{}^{\left(eq\right)}\left[\rho \left(x,t\right),u\right]$ ,

其中$\Delta u=F{\delta }_t/\rho$ ，F可以由下式表达：

$F=F_{int}\left(x\right)+F_s\left(x\right)+F_g\left(x\right)$ ,

其中$F_{int}\left(x\right)$ 表示对应x处位置粒子间的相互作用力，$F_s\left(x\right)$ 为对应x处位置固体表面与流体之间的相互作用力，$F_g\left(x\right)$ 为对应x处位置的所受引力。而$F_{int}\left(x\right)$ 表达式为：

$F_{int}\left(x\right)=-\beta c_0\psi \left(x\right)G\nabla \psi \left(x\right)-\left(1-\beta \right)c_{0}G\nabla {\psi }^2\left(x\right)/2$

其中$\beta$ 是加权因子，c₀是由晶格结构决定的常数，G为表示粒子间相互作用力强度常数，$\psi \left(x\right)$ 为有效密度。

$F_{int}\left(x\right)=-\beta \psi \left(x\right){\displaystyle {\sum }_{x^{\prime }}G\left(x,x^{\prime }\right)\psi \left(x^{\prime }\right)\left(x^{\prime }-x\right)}-\frac{1-\beta }{2}{\displaystyle {\sum }_{x^{\prime }}G\left(x,x^{\prime }\right){\psi }^2\left(x^{\prime }\right)\left(x^{\prime }-x\right)}$ ,

其中$G\left(x,x^{\prime }\right)$ 反映了相邻流体粒子间的相互作用强度，可由下式得出：

$G\left(x,x^{\prime }\right)=\{\begin{array}{l}{g}_1,\left|x-x^{\prime }\right|=1\\ g_2,\left|x-x^{\prime }\right|=\sqrt{2}\\ 0,其他\end{array}$ ,

在D2Q9模型中g₁ = 2g，g₂ = g/2。$\psi \left(x\right)$ 即有效密度可表示为：

$\psi \left(x\right)=\sqrt{\frac{2\left(p-\rho c_s^2\right)}{c_{0}g}}$ ,

在D2Q9模型中，取c₀为6.0。P可由状态方程得到，本文使用P-R状态方程：

$p=\frac{\rho RT}{1-b\rho }-\frac{\alpha {\rho }^2\epsilon \left(T\right)}{1+2b\rho -b^2{\rho }^2}$ ,

$\epsilon \left(T\right)={\left[1+\left(0.37464+1.54226\omega -0.26992{\omega }^2\right)\left(1-\sqrt{\frac{T}{T_{cr}}}\right)\right]}^2$ ,

$\omega$ 为无因次因子，$a=0.45724R^2{T}_{cr}^2/p_{cr}$ ，$b=0.0778R{T}_{cr}/p_{cr}$ ，其中T_cr表示临界温度，p_cr表示临界压力，固体与流体之间的相互作用力F_s为：

$F_s\left(x\right)=-\phi \left(x\right){\displaystyle {\sum }_i{G}_s}{\omega }_{i}s\left(x+e_i{\delta }_t\right)\cdot e_i{\delta }_t$

其中G_s为调节接触角的流固相互作用强度。S(x)为指标函数，当x为固体时，S(x) = 1，当x为流体状态时，S(x) = 0。

重力F_g可由重力加速度g，流体平均密度${\rho }_{ave}$ 表示：

$F_g\left(x\right)=g\cdot \left[\rho \left(x\right)\cdot {\rho }_{ave}\right]$

本文所采用的是D2Q9模型，粒子在平面内的运动方向一共可以描述为九种不同的方向，e₁~e₈如图2所示往不同方向运动，而e₀则是在原点保持不动。具体方向如下

$e_i\{\begin{array}{l}\left(0,0\right),i=0\\ \left(\pm 1,0\right)c,\left(0,\pm 1\right)c,i=1~4\\ \left(\pm 1,\pm 1\right)c,i=5~8\end{array}$

Figure 2. Spatial distribution map of D2Q9

图2. D2Q9空间分布图

3. 数值结果及分析

3.1. 加热器的高度

本小节研究加热器高度对池沸腾换热性能的影响。这里保证加热器长度不变，即加热器长度为200，加热器高度分别为30、50、70、90。

图3给出了不同加热器高度下得到的饱和沸腾曲线。从图中我们可以发现当壁面多热度较低(0.94 < T_w/T_c 0.95)时，各加热器换热性能近似。这是因为各加热器均处于自然对流阶段，系统吸热量不高，因此各加热器之间换热性能差异并不明显。而对于中壁面过热度(0.95 < T_w/T_c 1.00)高度为30时加热器换热性能最好。以上现象可以由图4所示的各高度下加热器在T_w/T_c = 0.97、0.99、1.01时的气泡图解释。如图4(a)所示，当加热器为30时，成核点位多，气泡脱离迅速，从而可以更好地换热。在高过热度(1.00 < T_w/T_c 1.02)时，高度为30的加热器换热性能急剧下降，而高度为90格子的加热器换热性能相对较好。这是因为如图4(a)所示高度为30的加热器在该过热度下表面形成气膜，气泡脱离困难，而高度为90的加热器仍处于核态沸腾(如图4(d)所示)，换热量持续增加。同时根据图中的结果，我们还能发现T_w/T_c = 0.99时，不同高度加热器率达到热流密度最大值时对应的壁面过热度差距很小，如高度为30，50，70和90的加热器分别在T_w/T_c = 0.99，1.00，1.00，1.00时达到热流密度最大点，即临界热流密度。

同时根据图4(b)中高度为50格子的加热器表面气泡运动状态可知，与加热器高度30时的情况一致，当T_w/T_c =1.01时加热器表面形成稳定的气膜。而根据图4(c)和图4(d)可知，随着加热器高度的增加，气泡脱离速度变慢，表面为不同温度下加热器表面上气泡的面积更大，这也解释了如图3中所示的在高过热度时热流密度随着高度增加而减小的现象。

Figure 3. Heat absorption curves for various heater heights at different degrees of overheating

图3. 各加热器高度在不同过热度下吸热量曲

(a)

(b)

(c)

(d)

Figure 4. The bubble states of heaters at different heights under different wall superheat conditions, from left to right: T_w/T_c = 0.97, 0.99 and 1.01. (a) At a height of 30 grids; (b) At a height of 50 grids; (c) At a height of 70 grids; (d) At a height of 90 grids

图4. 不同高度加热器在不同壁面过热度时的气泡状态，从左到右依次为T_w/T_c = 0.97、0.99和1.01。(a) 高度为30格子；(b) 高度为50格子；(c) 高度为70格子；(d) 高度为90格子

3.2. 加热器的表面润湿性

本小节研究加热器表面润湿性对池沸腾换热性能的影响。为了使研究结果更具有一般性，这里同时考虑了亲水和疏水情况，并选择两个加热器进行研究，其示意图如图5所示，其中左测加热器表面的接触角为θ₁，右测加热器表面接触角为θ₂。根据θ₁和θ₂的不同组合，我们考虑了三种不同情况，分别标识为Case 1，Case 2和Case 3。在Case1中，θ₁ = θ₂ = 110˚，Case2中，θ₁ = θ₂ = 70˚，而Case 3中，θ₁ = 70˚，θ₂ = 110˚，关于不同情况的详细参与见表1。在数值模拟中，加热器高度为20，长度为100，两加热器间距为20。

Figure 5. Schematic diagram illustrating the contact angle settings of two heaters (Red represents the heaters, while blue represents the liquid droplets)

图5. 两个加热器的接触角设置示意图(红色为加热器，蓝色为液滴)

Table 1. The contact angle settings of the two heaters

表1. 两个加热器的接触角设置

样本序号	接触角设置	润湿性
Case1	θ1 = 110˚，θ2 = 110˚	全疏水
Case2	θ1 = 70˚，θ2 = 70˚	全亲水
Case3	θ1 = 70˚，θ2 = 110˚	亲疏混合

图6给出了Case 1、Case 2和Case 3得到的饱和沸腾曲线。从图中结果可以发现当过热度较低时，如T_w/T_c 0.95时，Case 1即全疏水加热器情况得到的热流密度最大，Case 2即全亲水加热器得到的热流密度热度最小，而Case 3即一个加热器表面疏水而另一个加热器表面亲水对应的热流密度在疏水情况和亲水情况之间。通过观察图7所示的在沸腾过程中气泡成核和生长的演化过程可以发现，在低过热度时沸腾初期Case 1情况两个疏水表面都发生了成核现象，Case 2情况两个亲水表面都没有成核现象发生，而Case 3中成核只发生在疏水样本表面。这意味着疏水表面有利于气泡成核。再进一步观察气泡生长状态图可以发现表面为疏水性的加热器产生的气泡体积要大于亲水性加热器产生的气泡，这意味着疏水表面相变速率更快。以上结果解释了图6中所示的低多热度时疏水性表面热流密度的原因。

Figure 6. The heat absorption curves for different wetting properties between the heaters (Case 1: fully hydrophobic, Case 2: fully hydrophilic, Case 3: mixed wettability)

图6. 加热器之间不同湿润性的吸热量曲线(Case 1为全疏水，Case 2为全亲水，Case 3为亲疏混合。)

(Case1) (Case2) (Case3)

(a)

(Case1) (Case2) (Case3)

(b)

Figure 7. Samples with three different wettabilities: (a) nucleation state; (b) growth state (T_w/T_c = 0.95)

图7. 三种润湿性样本。(a) 成核状态；(b) 生长状态(T_w/T_c = 0.95)

随着壁面过热度的增加，热流密度变化趋势发生了改变，如在0.95 < T_w/T_c 1.00时，Case 2即全亲水性加热器所达到的热通量最大而Case 1即全疏水表面影响的热通量最小，这与低过热度时得到的结论相反。为了揭示引起这一现象的机理，图8展示了T_w/T_c = 0.99时三种润湿性样本的气泡图，我们可以发现：在该过热度下，Case 1表面产生很多气泡并连成一片，覆盖了部分加热面；而Case 2和Case 3加热器表面的气泡独立生长。这与图6中热流密度的变化趋势一致。除此之外根据图8在t^* = 89.33时不通过热度(T_w/T_c = 0.99)下的气泡状态图，我们还发现表面为亲疏混合性样本的气泡发生了偏离。这是因为两个加热器亲疏水性的不同而相互影响，从而发生气泡上升状态的不稳定。

(Case1) (Case2) (Case3)

Figure 8. The bubble state diagrams of the three wettabilities samples when T_w/T_c = 0.99, t^* = 89.33

图8. T_w/T_c = 0.99, t^* = 89.33时三种润湿性样本的气泡状态图

3.3. 加热器的长度与表面形状

本小节研究加热器长度与表面粗糙结构对池沸腾换热性能的影响，这里粗糙结构选择为锥形结构。在数值模拟中，加热器的高度为30，长度分别取100、150、200、250、300格子。

图9给出了光滑表面加热器和锥形表面加热器得到的饱和沸腾曲线。从图中结果可以发现，长度为100的加热器热通量都远高于其他长度的临界热流密度。这是因为随着长度的增加，加热面积也随之增加，从而使更多的汽化核心被激活，使气泡发生合并，直径增大，频率下降，进而影响换热性能。其次，我们从图中还可以发现不论加热器表面形状如何，在低过热度情况下加热器的换热性能相近，而过热度较高时长度较短的加热器换热性能优于其他加热器。

(a) (b)

Figure 9. The saturation boiling curves obtained for smooth surface heaters (a) and conical surfaces (b) at different lengths

图9. 光滑表面加热器(a)和锥形表面(b)在不同长度下得到的饱和沸腾曲线

(a) (b)

Figure 10. Bubble state diagrams for two sets of heaters with a length of 200. (a) Conical surface: T_w/T_c =0.99; (b) Smooth surface: T_w/T_c = 0.99

图10. 两组长度为200的加热器气泡状态图。(a) 锥形：T_w/T_c = 0.99；(b) 光滑：T_w/T_c = 0.99

再者，相同长度下锥形表面加热器的吸热量要明显高于光滑表面加热器。为了解释这一现象，图10展现了两组加热器在长度为100时的气泡状态图。如图所示，三角形加热器表面生成的气泡受到的扰动更大，且生成气泡的体积更小，数量更多，脱离速度更快。最后，图9中的数据还表明锥形表面加热器的临界热流密度与莱登弗罗斯点相比光滑表面加热器都要在较小的过热度达到，即整体沸腾曲线发生了右移。由以上结果可知，从加热器的长度来说，不论表面形状如何变化，过热度较低时换热性能相近，过热度较高时长度越短换热性能越好。从加热器的形状来说，表面形状不规则的加热器在相同的长度下换热性能更好。

4. 结论

本文基于气液相变格子Boltzmann模型研究了加热器的高度、长度、个数、润湿性和表面结构对池沸腾过程中气泡运动过程和大饱和沸腾曲线的影响，主要结论如下：

(1) 低过热度下不同加热器高度得到的热量密度差距较小，而在高过热度时，加热器高度较小时对应的池沸腾换热性能最好。

(2) 低过热度下全疏水性表面加热器性能较好，而高过热度下全亲水性表面加热器性能较好。

(3) 不论加热器表面结构如何，低过热度下各加热器换热性能近似，当过热度较高时加热器长度越短，换热性能越好。

(4) 加热器表面有微结构存在时气泡在成核和上升过程中受到的扰动更大，且与光滑表面相比，锥形加热器得到的大容器饱和沸腾曲线向右偏移。

基金项目

本研究由国家自然科学基金(51976128, 52376068)及上海市浦江人才(22PJD047)资助。

NOTES

^*通讯作者。