倾斜多孔介质方腔内纳米流体自然对流的格子Boltzmann 方法模拟*

张贝豪 郑林

(南京理工大学能源与动力工程学院, 南京 210094)

1 引 言

近几十年来人们在电力电子设备、太阳能集热器、核反应堆等领域的强化传热方面开展了广泛的研究, 而自然对流是该研究领域所涉及到的经典流动与传热问题之一. 在这些设备当中, 大多采用水、冷冻液和乙二醇等传统换热工质来进行强化换热, 与之相比, 在传统的纯液体中添加纳米颗粒,所得到的纳米流体具有较高的导热系数, 并且能满足一些特殊条件下的传热与冷却要求, 因而其强化传热特性引起了研究者的广泛关注[1].

国内外许多学者已经对方腔内纳米流体自然对流现象进行了研究. 例如: 李新芳和朱冬生[2]采用CFD 方法对影响纳米流体强化传热的因素进行了分析, 结果表明, 随着纳米颗粒体积分数的增大,流体的能量传输得到强化. 张晶等[3]采用4 种不同的粘度模型计算纳米流体粘度, 并对二维U 型腔内Al2O3-H2O 纳米流体的自然对流换热进行了数值模拟. Khanafer 等[4]数值模拟研究了二维方形腔内的纳米流体自然对流换热, 结果表明, 随着纳米颗粒体积分数的增加, 热效率增强. Hatami[5]采用有限元方法研究了不同纳米流体(TiO2和Al2O3)在加热肋片的矩形腔内的自然对流问题,结果表明, 随着体积分数的增大, 虽然Al2O3-H2O 纳米流体在热肋片处Nuave数缓慢增大, 但是TiO2-H2O 纳米流体在热肋片处Nuave数先增大后减小. Chen 和Du[6], Selimefendigil 和Oztop[7],Hatami[8]研究了纳米流体在方腔内的自然对流换热, 结果表明, 随着Ra 数的增加, 传热速率得到提高. Jahanshahi 等[9]采用瞬态热线法和哈密顿模型得到了不同体积分数下的流体导热系数, 分别研究了SiO2-H2O 纳米流体在热壁面处不同Ra 数下Nuave数随体积分数增加的变化趋势, 结果表明采用哈密顿模型计算得到的Nuave数呈现轻微减小的趋势.

多孔介质中的自然对流与换热问题普遍存在于自然界和工业应用中, 得到了学者们的广泛研究. 目前, 对于封闭多孔介质腔体内自然对流的研究大致可以分为两大类: 一类为通过水平多孔介质层的传热问题, 特别是水平多孔层存在高低温壁面, 即经典的Benard 流动问题[10], 另一类是具有两个不同温度的垂直壁面以及两个绝热的水平壁面的多孔介质传热问题[11]. Lan 和Prakash 等[12]采用Darcy-Brinkman-Forchheimer (DBF) 模型探究了不同达西数、Ra数和毕渥数下达西和非达西流体的流动机制. Kaluri 和Basak[13]研究了二维多孔介质方腔内纳米流体的自然对流换热问题.Alsabery 等[14]研究了倾斜方腔多孔层内纳米流体的自然对流问题, 结果表明施加较小的倾斜角可以显著提高壁面的Nuave数. Toosi 和Siavashi[15]对二维多孔介质腔内纳米流体自然对流现象进行了数值模拟, 结果表明同时采用多孔介质和纳米流体可以获得最佳的传热性能. 相比于均匀受热的边界条件, 近年来, 学者们对于壁面处于非均匀温度加热的研究不断增多. Wu 等[16]采用非达西模型的多孔介质控制方程, 分析了左右壁面均为正弦温度分布的自然对流传热特性. Basak 等[17]对比了多孔介质方腔内均匀边界条件和非均匀边界条件影响下自然对流强弱的变化情况, 研究指出, 底部中心加热壁面处为非均匀温度分布时, 换热效率更高. 前面的研究主要分析纯流体的自然对流问题,由于纳米流体具有强化传热特征, 考虑纳米颗粒影响的多孔介质方腔内自然对流的研究同样受到一些学者的关注. 比如: Ghasemi 和Siavashi[18]采用LBM 方法研究了不同线性温度分布的边界条件对多孔介质内纳米流体换热效率的影响. 此外,Sivasankaran[19]利用交替方向隐式的有限差分方法分析了下壁面为非均匀温度分布的倾斜多孔介质方腔内Cu-H2O 纳米流体的自然对流问题, 发现在Ra= 103时Nuave数随纳米颗粒体积分数的增加而减少.

综上所述, 以往研究表明纳米流体在多孔介质方腔内具有强化作用, 但是对于不同边界条件和计算模型, 其流动和换热过程表现出不同结论, 有的结果表现为强化换热效果[20], 有的则表现为削弱换热作用[21]. 值得注意的是, 非均匀温度边界条件对倾斜多孔介质方腔内的纳米流体流动与传热过程产生了重要影响. 与此同时, 由于Al2O3纳米颗粒在分散系中化学性质较为稳定且价格低廉, 故被广泛应用于工业强化传热领域中, 因此, 本文对非均匀温度边界条件下倾斜多孔介质方腔内Al2O3-H2O 纳米流体自然对流进行数值模拟, 分析了ϵ,Ra数,ϕ和g等参数对纳米流体流动传热特性的影响.

2 问题描述

2.1 物理模型

本文主要研究多孔介质方腔内纳米流体的自然对流问题. 如图1 所示, 方腔内纳米流体受到重力场的作用, 其下壁面与水平方向的夹角为g, 左壁面温度分布为T=sin(πy/H) , 而右壁面则保持低温Tc, 上下壁面绝热且不可渗透, 流体在固体壁面无滑移. 由于多孔介质内复杂的物理结构, 本文将通过以下4 种假设来简化物理问题: 1) 假设多孔介质为均质且各向同性; 2) 粘性耗散在能量方程中可以忽略不计; 3) 该流体被认为是不可压缩的牛顿流体; 4) 可以采用标准的 Boussinesq 假设.

图1 物理模型示意图Fig. 1. Schematic diagram of the physical model.

2.2 无量纲控制方程

基于上述基本假设, 本文采用DBF 模型, 对多孔介质方腔内纳米流体的自然对流问题模型的无因次形式的连续性方程、动量方程和能量方程进行描述[22]:

其中,X和Y是无因次横向纵向坐标,U和V是沿X和Y方向的无因次速度,Ra是瑞利数,Pr是普朗特数,Da是达西数,T是无因次流体温度,P是无因次压力,r是流体密度,µ是流体粘度,k是流体导热系数,f是纳米颗粒体积分数,g是方腔旋转角度,ϵ是多孔介质孔隙率. 下标f代表基液, 下标nf代表纳米流体, 下标m代表考虑多孔介质的有效物理量.

相应地, 无因次边界条件由以下表达式来确定:

上壁面:Y= 1

下壁面:Y= 0

左壁面:X= 0

右壁面:X= 1

2.3 纳米流体的物性参数

本文所研究的纳米流体是在基液水中加入一定体积分数的Al2O3纳米颗粒而形成的悬浮液. 假设纳米流体中的水和Al2O3处于热平衡状态, 并且纳米颗粒为球形结构且均匀分布在基液中, 不存在纳米颗粒的团聚和沉淀现象. 表1 所示为基液、Al2O3纳米颗粒和多孔介质的物性参数.

利用这些已知的热物性参数, 可以得到不同体积分数以及不同孔隙率下系统内的热物性参数[25,26],具体计算公式如表2 所示:

为描述热壁面处自然对流强度, 本研究组分析了壁面处局部Nu数和Nuave数, 其计算公式分别为[26]

表1 H2O, Al2O3 和玻璃纤维的热物理性质Table 1. Thermophysical properties of water,Al2O3 and glass fibers.

表2 纳米流体的热物性参数计算公式Table 2. Calculation formula for thermodynamic properties of nanofluids.

3 网格选择和程序验证

与传统的数值方法相比, 格子玻尔兹曼方法(lattice Boltzmann method, LBM)的介观性质和独特的计算特点使其具有演化过程清晰、计算效率高、易实现等优点[27], 已经被成功地应用到传热传质、湍流、多相流、多孔介质等复杂流体系统上[28−32].因此, 本文采用LBM 开展倾斜多孔介质方腔内纳米流体自然对流的数值模拟研究. 为了验证程序的可靠性, 首先对经典的多孔介质方腔自然对流问题进行网格无关性验证, 其中多孔介质方腔的上下壁面为绝热壁面, 左壁面为高温壁面, 右壁面为低温壁面, 不考虑倾斜角变化, 流体在固体壁面处无滑移. 在该模拟研究中, 保持Da= 107,ϵ= 0.9999,Ra= 106,Pr= 6.2 不变, 并采用四种均匀网格(80 × 80, 100 × 100, 120 × 120, 140 × 140)进行网格无关性验证, 得到了不同网格数下热壁面处的Nuave数. 与文献[33]对比的结果如表3 所示,随着网格数的增大,Nuave数的相对误差逐渐减小,当网格数为120 × 120 时,Nuave数的相对误差为0.83%, 综合考虑结果准确性以及计算效率, 本文所有算例将选用120 × 120 作为计算网格数.

为了进一步验证所采用LBM 的可行性, 对二维方腔内的自然对流问题进行数值模拟, 其中, 方腔上下壁面为绝热壁面, 左壁面为高温壁面, 右壁面为低温壁面, 流体在固体壁面处无滑移. 首先验证忽略多孔介质的影响, 即ϵ= 0.9999, 并在模拟时取Pr= 0.72 和Da= 107, 计算热壁面处的Nuave数, 并与已有文献[33]的结果进行对比. 从表4 可知, 本文的计算结果与文献[33]的计算结果吻合得非常好.

表3 不同网格数与文献[33]的Nuave 数比较Table 3. Comparison of Nuave number with literature[33] in different grids number.

表4 本文与文献[33]的Nuave 数值结果的比较Table 4. Comparison of Nuave number with previous literature[33].

其次, 考虑多孔介质的影响, 采用DBF 模型进行数值模拟计算. 在模拟过程中保持ϵ= 0.6,Pr= 1.0 和Da= 10–2不变, 同样, 在表5 中对比了不同Ra数情况下热壁面处的Nuave数. 结果表明, 本文的结果与文献[28]的吻合得比较好, 误差大约在3%以内.

表5 本文与文献[34]的Nuave 数值结果的对比Table 5. Comparison of Nuave number with previous literature[34].

4 结果与讨论

4.1 孔隙率 ϵ 对方腔内流动换热的影响

本文数值模拟了多孔介质方腔内Al2O3-H2O 纳米流体的自然对流问题, 分析了不同物理参数对纳米流体换热效率的影响. 在模拟过程中,无量纲参数的选取为Da= 0.01,f= 0.01,Ra=105,Pr= 6.2,ϵ= 0.9 和g= 0°.

由于孔隙率与多孔介质渗透率有关, 它的大小会影响多孔介质方腔内纳米流体的流动状态和传热效率. 如图2 所示, 当ϵ= 0.3 时, 流场中心处形成了一个不规则的椭圆形流胞, 方腔中心区域内等温线与水平壁面存在一定倾斜, 并且左右壁面温度边界层较厚, 因此其换热方式以导热为主; 随着孔隙率的增大, 当ϵ= 0.5 时, 椭圆形主流胞变得更加扁平, 温度边界层逐渐被挤压到左右壁面, 其对流换热强度增大; 当ϵ= 0.7 时, 流体受到的约束阻力进一步变小, 使得其对流换热区域扩大; 当ϵ= 0.9 时, 即在高孔隙率条件下, 流胞在方腔中心处转变为扁平的椭圆形流胞, 并且主流胞中心存在分裂成两个小流胞的趋势, 方腔中心区域存在大量与水平壁面平行的等温线, 因此方腔内对流换热占据主导作用. 值得注意的是, 由于左侧垂直壁面为正弦温度分布, 在方腔左上角位置存在壁面吸收热的现象, 导致在左上角处形成了一个相反于主对流胞的逆时针小漩涡.

为了分析孔隙率对流体流速的影响, 本研究组还比较了不同ϵ下热壁面处局部V和上壁面处局部U的分布图象. 如图3 所示, 随着孔隙率的增大, 热壁面附近的竖直流速明显增大; 同样地, 上壁面附近的水平速度也存在相同的增长趋势. 由于为正弦温度分布, 热浮升力在热壁面中点位置最强, 因此竖直速度的峰值在Y = 0.5 附近; 与此同时, 由于左上角小流胞的阻碍作用, 上壁面附近水平速度的峰值出现在X = 0.8 附近. 值得注意的是, 由于方腔的左上角存在一个逆时针旋转的小流胞, 因此, 在其附近的竖直和水平速度均为负值.综上所述, 随着孔隙率的增大, 流体受到多孔介质的阻力减弱导致流速增强, 提升了热壁面处的换热效率.

图2 不同 ϵ 下温度场和流场的分布图像 (a) ϵ = 0.3; (b) ϵ = 0.5; (c) ϵ = 0.7; (d) ϵ = 0.9Fig. 2. Streamlines, isotherms contours for different ϵ : (a) ϵ= 0.3; (b) ϵ= 0.5; (c) ϵ= 0.7; (d) ϵ = 0.9.

图3 (a)不同 ϵ 下X = 0 处的竖直速度分布; (b) Y = 1 处的水平速度分布图Fig. 3. (a) Vertical velocity distribution at X = 0; (b) horizontal velocity distribution at Y = 1 for different ϵ .

图4 (a)不同 ϵ 下热壁面处Nuave 数分布曲线; (b)热壁面处局部Nu 数分布曲线Fig. 4. (a) At the heated wall Nuave number; (b) local Nu number for different ϵ .

如图4(a)所示, 随着孔隙率的增大, Nuave数呈现上升趋势, 然而这种增长幅度逐渐减弱, 这一变化规律与以往均匀温度分布所得的结论相符[35].这是由于在较高孔隙率的情况下, 方腔内部多孔介质的阻碍作用较小, 孔隙率对方腔内传热效率的提升不明显. 为了进一步分析孔隙率对热壁面附近局部换热特性的影响, 本研究组还分析了热壁面处的局部Nu 数随孔隙率变化的情况, 如图4(b)所示,热壁面0.2 < Y < 0.6 为主换热区域, 不同孔隙率之间的变化差异并不明显. 与此同时, 在热壁面0.6 < Y < 0.8, 随着 ϵ 的增大, 局部Nu 数逐渐减小, 这是由于位于方腔中心的主流胞对流逐渐增强, 流胞形态逐步转变为偏平椭圆形, 使得温度边界层沿着主流胞流动方向扩散, 因此, 在该区域内温度边界层逐渐变厚, 局部换热效率下降. 但是在0 < Y < 0.2 和0.9 < Y < 1 靠近上下壁面区域附近, 其变化趋势恰好相反, 这是由于最大温度在热壁面中点位置, 向上或向下壁面延伸时热壁面温度逐渐减小, 在上下壁面局部形成了左侧低温右侧高温的现象, 因此在方腔左上角和左下角出现了热壁面向方腔内部吸热的现象, 并且形成了逆时针方向的小漩涡. 随着孔隙率的变大, 主流胞逐渐扩展至几乎整个方腔区域, 上下两个小流胞逐渐被挤压变小, 向左换热的区域逐渐被压缩并且换热效率变弱, 局部Nu 数的绝对值逐渐变小主导了Nuave数的发展趋势, 因此, 随着孔隙率的增大, 流胞形态发生转变, 减少了热壁面从方腔内部的吸热量, 进一步促进了方腔内的对流换热效率, 使得热壁面处的Nuave数不断增大.

4.2 Ra 数对方腔内自然对流的影响

在自然对流问题中, Ra 数是表征自然对流综合强弱的重要无量纲准则数, 为了研究Ra 数对方腔内自然对流的影响, 保持 ϵ =0.9, Da = 0.01, f =0.01, Pr = 6.2 不变, Ra 数在103—106变化. 如图5所示, 当Ra = 103时, 对流换热现象几乎可以忽略不计, 圆滑的主流胞集中在靠近热壁面区域内,大量等温线呈现上下对称并垂直于水平方向, 温度梯度大部分集中在方腔左半侧区域内, 因此该状态下流场内主要以导热为主; 当Ra = 104时, 方腔中心处形成了一个椭圆形状的对流胞, 由于浮升力的增强, 流体在进入热壁面下侧时促进了对流换热, 而流体在进入热壁面上侧时等温线将沿着流动方向扩散, 导致局部换热能力减弱, 此时方腔内换热方式仍然以导热为主; 当Ra = 105时, 方腔中心处存在大量水平方向的等温线, 方腔中心的流胞开始出现向多个流胞分裂的趋势, 因此对流区域内主要以对流换热为主; 当Ra = 106时, 此时方腔中心处的椭圆形流胞分裂成2 个不规则小漩涡, 主流胞占据了几乎整个多孔介质方腔, 边界层被挤压变薄, 并且边界层附近存在很多沿流体流动方向的等温线. 值得注意的是, 随着Ra 数的不断增大, 方腔左上方的对流胞在逐渐增大, 其中心位置不断向右侧移动, 挤占主流胞的空间.

本研究组还比较了不同Ra 数下热壁面处竖直速度和上壁面处水平速度分布图像. 如图6(a)和图6(b)所示, 在低Ra 数时热壁面处的流速十分微弱, 在高Ra 数时热壁面处流体流速得到了显著提升, 这意味着Ra 数的增加将直接增强方腔内自然对流的强度, 腔体内流胞环流变得越来越剧烈,流动速率加快, 温度边界层变薄, 强化了壁面处的传热能力, 使得方腔内部温度变得越来越均匀. 随着Ra 数的增大, 在主换热区域流速增大的同时,方腔左上角小漩涡的流动速率也得到增大, 并且速度为负的区域明显扩大.

如图7(a)所示, 均匀温度分布模型表明热壁面的Nuave数随Ra 数增大而增大[20], 但是随着Ra 数的增大, 当Ra = 104时, 热壁面Nuave数出现了微弱减小现象. 为了进一步解释这一现象, 本研究组分析了局部Nu 数的分布曲线. 如图7(b)所示, 随着瑞利数的增大, 局部Nu 数的峰值位置集中在0.2 < Y < 0.5 附近, 这是由于热浮升力的增强使得流速加快, 下侧热边界层变薄, 并且在高瑞利数时该换热区域占据主导地位, 因此, 在Ra =105和106时, 局部Nu 数产生了显著的增长趋势.

图5 不同Ra 下温度场和流场的分布图像 (a) Ra = 103; (b) Ra = 104; (c) Ra = 105; (d) Ra = 106Fig. 5. Streamlines, isotherms contours for different Ra number: (a) Ra = 103; (b) Ra = 104; (c) Ra = 105; (d) Ra = 106.

图6 (a) 不同Ra 下X = 0 处的竖直速度分布; (b) Y = 1 处的水平速度分布图Fig. 6. (a) Vertical velocity distribution at X = 0; (b) horizontal velocity distribution at Y = 1 for different ϵ .

但是, 在0.6 < Y < 0.8, 当Ra = 103和104时局部会出现Nu 数被削弱的现象, 究其原因, 一方面,是由于流胞逐渐转变为扁平的椭圆形, 对左上侧的热边界层挤压作用减弱; 另一方面, 随着瑞利数的增大, 中心流胞流速增大并转变为椭圆形状, 由于流体流速的增强, 该区间内等温线沿流胞主流区域流动方向扩散, 使得该局部温度边界层变厚. 此外,由于热壁面为非均匀边界条件, 在上下壁面附近处出现了局部Nu 数小于0 的情况, 随着瑞利数的增加, 中心区域的换热效率增强, 同时上壁面附近向左的热通量也得到加强, 而下壁面受到主流胞挤压作用向左侧的换热效率微弱下降.

通过上述分析发现, 当Ra = 104时, 热壁面上侧局部边界层沿流体流动方向扩散, 局部换热效率下降显著, 并且低温区域吸热现象随Ra 数增大而增大, 导致左壁面的Nuave数出现了减小现象.与此同时, 在高Ra 数情况下, 由于方腔内主换热区域换热效率增长更加明显, 其壁面Nuave数增长显著, 受非均匀温度边界层影响较小.

图7 (a)不同下Ra 热壁面处Nuave 数; (b)热壁面处局部Nu 数分布曲线Fig. 7. (a) At the heated wall Nuave number; (b) local Nu number for different Ra.

4.3 倾斜角g 对流体流动换热的影响

当多孔介质方腔逆时针旋转一定角度时, 流体的流动形态和换热强度将会发生明显变化. 为了分析 倾 斜 角g 的 影 响,保 持 ϵ = 0.9, Da = 0.01,ϕ=0.01 , Ra = 105, Pr = 6.2 不变, 倾斜角g 在0° ≤ g ≤ 120°变化. 从图8 可知, 当g ≠ 0°时, 方腔中心流胞的流线更加圆滑, 主流胞也更加向中心集中, 并且方腔中心为单一流胞, 未出现分裂成多个流胞的趋势. 从温度场可知, 当方腔施加倾角时,在方腔中心形成了大量沿流胞主流速方向的弯曲等温线, 但是其温度边界层的变化在宏观云图中变化不明显.

因此, 本研究组计算了水平中心线的局部温度分布以及热壁面处Vave/Nuave的对比图, 其中,Vave为沿Y 方向的平均速度. 如图9(a)所示, 与g = 0°相比, 热壁面附近温度下降更为迅速, 这意味着在有倾斜角的情况下热壁面中点附近的温度边界层变薄, 换热效率得到一定加强. 深入分析热壁面整体的换热效率. 如图9(b)所示, 当g =40°时, 发现壁面处的换热效率得到了较为显著的提升; 当g = 80°时, Nuave数存在削弱换热现象;当g = 120°时, 其强化传热作用较为微弱, 这一变化规律与以往文献[36]所得结论一致. 由于倾斜角的变化使得热浮升力的方向发生改变, 流体在X 方向上存在浮升力的驱动, 影响了热壁面处Vave的大小, 并且壁面处局部V 会直接影响温度梯度以及壁面处传热强度, 因此Vave/Nuave的变化趋势基本相同, 并呈现上下波动的趋势.

为了进一步探究g 对流体传热的影响, 本研究组分析了热壁面处局部V 的分布曲线, 由于g =120°时旋涡流动方向为逆时针, 本文选取它的相反数作为参考. 从图10(a)可以看出, 当g = 40°时,其局部V 的峰值存在最大值; 当g = 80°时, 其局部V 的峰值最小, 并且绝大多数区域小于g = 0°时的数值; 当g = 120°时, 速度峰值数值与g = 0°时几乎接近, 因此g = 40°时强化效果最为明显,而g = 80°时则产生削弱现象. 最后, 为了进一步细致分析g = 120°时局部Nu 数的分布规律, 本研究组给出了不同倾角下局部Nu 数的分布曲线. 如图10(b)所示, 当g = 0°, 40°和80°时, 非均匀温度边界条件在中点以上位置出现了与主换热区域相反趋势的分布情况, 而最大局部Nu 数存在于热壁面中点以下区域. 这是由于当流体流动方向为顺时针时, 中点以下的流体受到浮升力和主流包挤压作用的影响, 对流换热作用最显著. 由图10(b)可知,主换热区域内的换热强度由强到弱依次为40°,0°和80°. 当g = 120°时, 流胞旋转方向的转变使得热壁面最大局部Nu 数出现在壁面中点以上的位置, 而局部Nu 数的峰值介于0°和40°之间.

对方腔施加倾角后, 热浮升力的有效作用发生改变, 使得壁面的竖直速度发生改变, 而壁面Vave与Nuave数关联性增强. 当g = 40°时, 通过提升主换热区域的竖直速度可以提高该区域内局部Nu 数, 使得整体壁面换热效率增强. 类似地, 当g = 80°时, 壁面竖直速度的削弱对流体传热效率的影响较为明显, 导致壁面换热效率下降显著. 这一规律与文献[36]所得结论相符, 表明主换热区域对壁面Nuave数影响较为明显, 壁面低温区域从方腔内吸热的强度受倾斜角的影响较小, 这是由于主换热区域和主流区域近似重叠, 因此, 随倾斜角的变化, 均匀或非均匀边界条件下主换热区域内Nuave数强化或削弱作用都将起到主导作用.

图8 不同g 下温度场和流场的分布图像 (a) g = 0°; (b) g = 40°; (c) g = 80°; (d) g = 120°Fig. 8. Streamlines, isotherms contours for different g number: (a) g = 0°; (b) g = 40°; (c) g = 80°; (d) g = 120°.

图9 (a)不同g 时Y = 0.5 处局部温度分布曲线; (b)热壁面处Vave/Nuave 的分布曲线Fig. 9. (a) Local temperature distribution along the Y = 0.5; (b) average velocity in the y direction & Nuave number at the heated wall in different g.

图10 (a)不同g 下热壁面处局部竖直速度V; (b)热壁面处局部Nu 数的分布曲线Fig. 10. (a) Local velocity in the y direction; (b) local Nuave number at the heated wall in different g.

本研究组还研究了不同倾斜角下孔隙率的增加对方腔内强化传热的影响. 如图11(a)所示,在施加倾斜角的情况下, 随着孔隙率的增大, 当ϵ = 0.9 时, Nuave数出现轻微减小; 而没有施加倾角时, Nuave数随 ϵ 的增大而增大. 当g = 40°时,Nuave数的增长幅度最大, 当 ϵ = 0.7 时, Nuave数达到 峰 值,与g = 0°的 情 况 相 比, Nuave数 在 ϵ =0.7 时提升了17.72%.

当方腔不施加倾角时, 壁面处的Nuave数随孔隙率增加而增加, 但是在施加倾斜角条件时, 方腔热壁面的Nuave数出现先增加后减小的趋势, 为了解释这一现象, 以g = 40°和0°时孔隙率为0.7 和0.9 的情况为例进行分析. 如图11(b)所示, 当g =0°时, 在主换热区域内(0.2 < Y < 0.6), 孔隙率增大, 局部Nu 数减小的幅度较小, 而位于上壁面附近区域内, 局部Nu 数的增长成为主导. 局部Nu 数是由导热系数比值km/kf和热壁面温度梯度所组成, 当g = 0°时, 左壁面温度梯度的增加幅度相对于km/kf减小幅度更加显著, 因此, 随着孔隙率的增加, Nuave数呈现出单调递增的趋势. 当g =40°时, 方腔内对换热能力起主导作用的依然位于主换热区域(0.2 < Y < 0.6), 在此区域内, 随着孔隙率的增大, 局部Nu 数减小的幅度更为明显, 从而导致了孔隙率为0.9 时Nuave数的轻微减小. 这是由于随着孔隙率的进一步增大, 多孔介质导热对方腔传热效率的贡献逐渐变小, 即km随孔隙率增大而减小, 而施加不同倾斜角后, 左壁面处的温度梯度在0.7 < ϵ < 0.9 增长的幅度减缓, 倾斜角的变化使得km/kf的减小幅度相比于左壁面温度梯度的增加幅度更加显著, 因此加热壁面的局部Nu 数出现了减小的情况. 值得注意的是, 当g = 80°与g = 0°进行对比时可知, 热壁面处的Nuave数在ϵ =0.9 时产生了削弱的作用. 这可能是由于在低孔隙率的情况下, 沿X 方向的作用力可以更有效地促进流场的对流换热, 而随着孔隙率的增加, 方腔内流速加快, 在热壁面沿Y 方向上热浮升力不再起到主导作用, 从而使得热壁面的换热效率低于g = 0°时的换热效率. 通过分析发现, 在高孔隙率时对方腔施加不同的倾角, 随孔隙率增大, 壁面温度梯度增长幅度以及多孔介质对流体流速的影响减弱, 使得热壁面处的传热能力得到抑制, 并且当ϵ =0.7 时左壁面处的Nuave数达到峰值, 这表明在高孔隙率时, 壁面处的Nuave数受到km/kf减小的影响将更加明显.

图11 (a)随着 ϵ 的增加不同g 时热壁面Nuave 数分布曲线; (b)当g = 0°, 40°时, 不同 ϵ 下局部Nu 数的分布曲线Fig. 11. (a) Variation of Nuave number as a function of ϵ in different g at the heated wall; (b) when g = 0°, 40°, variation of local Nu number at the heated wall in different ϵ .

图12 (a)随着f 的增加不同g 下热壁面Nuave 数分布曲线; (b)当g = 0°, 40°时, 不同f 下局部Nu 数的分布曲线Fig. 12. (a) Variation of Nuave number as a function of f in different g at the heated wall; (b) when g = 0°, 40°, variation of local Nu number at the heated wall in different f.

最后, 本研究组分析了不同g 下随着纳米颗粒体积分数的增加, 方腔内部的换热机理. 如图12(a)所示, 随着纳米颗粒体积分数的增加, 当g = 40°,80°和120°时, 热壁面处的Nuave数存在微弱增大的趋势, 且当g = 40°时, Nuave数最大, 强化换热现象最明显, 而当g = 0°时, Nuave数反而出现了微弱削弱的趋势, 且 ϕ=0.04 较 ϕ =0 时的Nuave数微弱减小了0.32%. Ho 等[20]采用同样的导热系数和粘度公式进行方腔内数值模拟, 发现随着Al2O3-water 纳米流体体积分数的增加, 热壁面的Nuave数缓慢增加. 导致本文的结果与文献[26]的差异主要由以下两个方面原因: 一方面, 由于本文采用非均匀温度边界, 在增强高温壁面换热效率的同时, 由于靠近上下两个壁面处存在从方腔内部吸热的现象, 这部分区域也随之被加强并在g = 0°时占据主导趋势; 另一方面, 在g = 0°时, 随体积分数的增加流体粘度变大, 热壁面的平均温度梯度缓慢减小, 这一减小幅度相比km/kf的增长幅度更为显著, 导致Nuave数出现缓慢下降的现象.

但是, 当对方腔施加任意倾角时, 热壁面Nuave数缓慢增长, 为了进一步分析这种不同变化趋势的原因, 在图12(b)中, 以g = 0°和40°, 纳米颗粒体积分数为0%和4%的情况为例进行分析,一方面, 当g = 40°时, 在主换热区(0.2 < Y < 0.6)附近, 随着孔隙率的增大, 局部Nu 数随之增长;另一方面, 虽然随体积分数的增大流体, 有效粘度逐渐变大导致流体流速降低, 对流强度缓慢下降,但是施加不同的倾斜角后, 壁面温度梯度呈现微弱增加趋势, 而km/kf增长幅度固定不变, 使得在g = 40°时Nuave数单调递增趋势更加明显. 这说明通过对方腔施加一定倾斜角, 改变了壁面温度梯度变化趋势, 进而改变了其Nuave数. 因此, 热壁面的换热效率取决于g 和 ϕ 对壁面温度梯度的影响.

5 结 论

纳米流体作为良好的换热媒介被广泛地应用到强化传热领域, 很多因素会影响其强化传热特性, 其中包括多孔介质孔隙率、纳米流体体积分数、方腔倾斜角度以及边界条件等因素. 本文以非均匀温度边界条件为例, 采用LBM 系统研究方腔内Al2O3-H2O 纳米流体的对流传热特性, 重点讨论了不同g 对腔内换热效率的作用, 以及 ϵ和 ϕ 对不同倾斜角下多孔介质方腔内流动换热的影响. 通过对非均温度边界条件的情况研究表明, Ra 数、孔隙率和纳米颗粒体积分数并不像均匀温度边界条件那样始终起到强化传热的效果.

主要结论如下:

1)在低Ra 数时, 流体受非均匀温度边界层低温区域影响出现吸热现象, 热壁面处的Nuave数出现轻微减小的趋势, 这与以往采用均匀温度分布的文献[20]所得结论存在差异; 在高Ra 数时, 由于瑞利数综合表现为热浮升力的强度, 方腔内主换热区域换热效率增长更加明显, 热壁面处Nuave数增大最为显著, 受到非均匀温度边界层影响较小.

2)随着g 的增加, 方腔内热浮升力随之发生改变, 通过对壁面流速的改变, 对其壁面强化换热效果产生不同影响, 当g = 40°时, 左壁面的Nuave数增长幅度较为明显; 而当g = 80°时, 左壁面的Nuave数反而产生削弱作用. 这一规律与文献[36]所得结论相符, 这是由于主换热区域和主流区域近似重叠, 而壁面Vave与 Nuave数关联性较强, 因此壁面Nuave数受主换热区域影响较为明显, 而壁面吸热区域受到倾斜角的影响较小.

3)均匀温度边界条件下的研究表明, 壁面处的Nuave数随孔隙率增加而增加[12]; 而非均匀温度边界条件下的结果表明, 在高孔隙率对方腔施加不同倾角时, 增大孔隙率对壁面温度梯度的促进作用减弱, 因此, 孔隙率的增加对热壁面处的传热能力存在抑制作用.

4)对于均匀温度分布边界条件下的研究表明,热壁面的Nuave数随Al2O3-H2O 纳米流体体积分数增加而缓慢增加[20], 而本文非均匀温度边界条件下的研究表明, 当g= 0°时, 随ϕ的增大流体受粘性增长的影响较为明显, 热壁面温度梯度缓慢下降,Nuave数存在微弱减小的趋势. 而施加倾斜角后, 热壁面温度梯度受流体粘性增长的影响减弱,导致随ϕ的增大热壁面温度梯度微弱增大, 使Nuave数呈现微小增长趋势.

通过上述总结归纳发现, 在本文所研究的非均匀温度边界物理模型下, 纳米颗粒体积分数所产生的强化换热作用不再明显. 给定Ra数时, 要有效改善纳米流体流动换热效率, 就需要利用多孔介质对有效导热系数的提升, 以及倾斜角对系统内的扰动作用, 采用合适的孔隙率和方腔倾斜角度对方腔进行干预.