二维局部受热腔体内电热对流问题模拟和分析1)

刘镇涛 肖 莉 和 琨 汪 垒 ,2)

* (中国地质大学(武汉)数学与物理学院,武汉 430074)

† (中国地质大学(武汉)数学科学中心,武汉 430074)

引言

电流体动力学(electro-hydrodynamics,EHD)属于电磁学和流体动力学等多学科交叉领域,它描述了带电粒子在库仑力驱动下引起的电对流[1].作为EHD 的基本课题之一,电热对流(electro-thermo convection,ETC)考虑电场、流场和温度场等多个物理场耦合及带电粒子与流体介质之间的相互作用,也是EHD 的重要研究课题之一[2-3].近年来,随着能源问题的日益突出,节能环保等相关问题越来越受到广泛的关注,科研人员对于电热对流强化传热应用逐渐重视[4-5].ETC 与许多其他传统强化传热方法相比,具有节约资源、提高传热效率、结构简单和安全等优点,在强化传热领域受到许多研究人员的关注[6].

近年来众多学者对ETC 问题开展了相关研究,早期对ETC 的研究主要是对简化的理论模型和实验数据进行分析.Atten[7]的研究结果表明电热对流可以显著提高传热效率,传热效率提升了15 倍.Mccluskey 和Atten[8]的研究结果显示,电热对流对传热的影响比纯自然对流要高一个数量级.Wu 等[9]研究了电热对流中亚临界或超临界的分叉类型和有限幅度稳定性标准,并发现了超临界点和亚临界点的取值很大程度上取决于普朗特数Pr和无量纲离子迁移率M.随后,Traore 等[10]提出一种求解ETC 控制方程的数值模拟方案,该方案能够很好解决ETC 问题的多物理场耦合问题,基于该模型,作者还对二维腔体中的ETC 进行了数值研究.Dantchi 等[11]采用有限体积法求解ETC 方程,模拟中考虑了二维腔体内从下壁面电极和从侧壁面电极注入电荷两种情况,结合相关无量纲参数进行分析,研究结果表明电场增强了传热,对于足够大的T值,努塞特数与瑞利数无关.Luo 等[12]设计了格子Boltzmann 模型来求解ETC 方程,使用4 个不同的分布函数分别计算流场、电势、电荷密度分布和温度场,还进行了多尺度分析来将离散的格子Boltzmann 方程恢复到宏观控制方程,并通过对几种问题进行验证,最后得到的结果与解析解等其他数值结果一致.除此之外,Lu 等[13]采用了高分辨率逆风方案(high resolution upwind scheme)探究了方腔内介电液体中电热对流的流动结构和传热现象,通过研究相关无量纲参数如热瑞利数、无量纲离子迁移率以及普朗特数,最后发现了多种流动结构以及传热规律.Hu 等[14]采用了格子Boltzamnn 方法研究了不同几何结构下的二维电热对流问题,探讨了瑞利数、电瑞利数、长宽比和普朗特数等无量纲参数对流场流线、等温线、电荷密度分布和界面的平均努塞特数的影响,结果表明增加电瑞利数可以提升传热效率并改变了流场的结构.Wang 等[15]提出了一种谱元法(spectral element method)的新型数值方法,在验证该数值方法能有效的解决电热对流问题后,又模拟了正弦温度边界方腔内的二维电热对流问题,研究结果表明,与均匀温度边界相比,正弦温度边界能提高传热效率.还有一些关于电热对流的拓展研究.Chand 等[16]研究了两平行电极板之间的纳米流体在电场影响下的热不稳定性.Srivastava[17]对多孔介质中的二元流体电热对流进行了稳定性分析.

上述研究在很大程度上加深了人们对电热对流在传热过程中的内在机理.但是应当指出的是,上述研究,模型大多数是水平放置的两平行电极板[7,9,12,18]或者是局部单极注入电荷[19-20],这些研究中较少考虑温度分布不均匀的情况,实际上在许多研究领域如对自然对流的研究,非均匀温度边界问题一直是热门的研究课题.Poulikakos[21]研究了二维方腔左侧壁面上半部分加热,下半部分冷却,其他壁面隔热的情况,结果表明,当加热区域位于冷却区域下方时,整体的传热增强.Sathiyamoorthy等[22]研究了二维方腔两侧壁面的温度是具有线性变化的情况,结果表明,当壁面两侧的温度呈线性变化时,方腔底部的努塞特数会变得十分不稳定.此外,Bilgen 和Yedder[23]研究了二维矩形腔体垂直壁面的温度呈正弦分布而其他壁面绝热的情况.另外,考虑加热器长度和位置的情况,Oztop 和Abu-Nada[24]研究了部分加热的二维腔体中纳米流体的自然对流,发现了加热器的位置会影响流线结构和温度分布.最近,Wang 等[25]研究了局部加热的二维方腔中纳米流体的自然对流,结果表明传热效率随着加热器的长度增加而降低.从这些文献中不难看出,不论是在电热对流中考虑局部单极注入电荷还是在自然对流中考虑非均匀温度边界情况下,都会对传热效率产生影响,并产生复杂的流动结构.而电热对流本质是电对流和自然对流耦合的问题,因此对于局部单极注入电荷和局部加热相结合的机制需要进行更深入的研究,以便于理解非均匀温度边界的电热对流的流体流动和传热特性.

本文采用了格子Boltzmann 方法(LBM)求解控制方程组,其中包括Navier−Stokes 方程,电荷密度守恒方程,电势的泊松方程以及温度方程,另外,为了减少计算时间以及提升计算效率,本文使用NVIDIA图形处理平台,编写计算统一设备架构程序.研究了在局部受热条件下,方腔内流体流动和系统传热效率的变化情况,并为研究其他非均匀温度边界的电热对流问题提供参考数据.

1 问题描述

1.1 物理模型和控制方程

如图1 所示,介电液体被限制在长度为H的二维方腔中.左侧壁面是加热的电极板,温度值固定为θ0,电势保持为 φ0,长度为h,电极板的中心点到下壁面的距离用 δ 来表示.右壁面温度固定为 θ1,电势保持为 φ1,且 θ0>θ1,φ0>φ1,上下壁面是绝热且绝缘的.为方便模拟,假设方腔中的流体是不可压缩的牛顿流体,流体的密度变化是满足标准的Boussinesq假设的,自由电荷的产生只发生在介电液体与电极界面之间的电化学反应,并且注入的电荷密度取恒定值为q0.另外,忽略了黏性热耗散以及焦耳热的影响[26].

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

基于上述假设,二维局部受热腔体内电热对流问题可以通过连续方程、动量方程、电势的泊松方程、电场的定义方程、电荷密度守恒方程和温度方程共同描述[27]

在上述方程组中,u,E,g,q,θ,θref,φ 分别为流体速度、电场强度、重力加速度、电荷密度、温度、参考温度、电势.p,υ,β,ε ,K,D,χ 分别代表压力、运动学黏性、热膨胀系数、介电常数、离子迁移率、电荷扩散系数以及热扩散系数.为了更方便的研究电热对流问题,引入以下无量纲参数[28]

其中,Ra为热瑞利数,表示浮升力与黏性力的比值.T为电瑞利数,表示库仑力与黏性力的比值.对于热瑞利数Ra以及电瑞利数T而言,在实际实验中可依据温度差和电压差灵活设置[29],为了系统探究其影响,本文将之分别约束在 0 ≤T≤1200,1.0×104≤Ra≤1.0×105区间.C为电荷注入强度,罗康[30]的理论研究和数值模拟研究中,一般取C=10.0作为其强注入的代表性值,为此本文中所涉及的电荷注入强度C在未加强调的情况下,其值均为10.0.无量纲离子迁移率M是所谓的水动力学迁移率与真实迁移率的比值.对于带有注入粒子的介电液体而言,实验显示其值一般大于等于3[3],例如甲醇(H+,M=4.1)、氯苯(C l−,M=4.9)、乙醇(H+,M=8.8;C l−,M=26.5)、硝基苯(C l−,M=22)、碳酸丙烯(C l−,M=51),鉴于此,为了更加清晰地认识无量纲离子迁移率M对于电热对流的影响,本文所考察的无量纲离子迁移率M介于5.0 和50.0 之间.α 为无量纲电荷扩散系数,其值一般介于 1 0−3和1 0−4之间[28],在本文中参考了文献[31],将其值固定为 1 .0×10−4.Pr为动量边界层与热边界层之比,反映了流体物理性质对热传递的影响,本文重点关注一类液体形态的介电介质[3],因此将其值设置在5.0 和50.0 之间.另外,还考虑了电极板的长度h和电极板的中心点到下壁面的距离δ,分别取值为 0 .5H≤h≤H,0 .25H≤δ ≤0.75H.对于温度边界、电势边界、电荷密度边界以及速度边界,它们的数学描述如下,左侧壁面电极板边界条件如下

而对于速度边界而言,方腔的4 个壁面均为无滑移边界

2 数值方法

在数值方法选择上,本文采用LBM 进行模拟,与有限体积法和有限元方法等传统方法不同,LBM 在处理复杂边界问题时具有高效和简单的优点[32],又因为LBM 具有天然并行性[33],所以使得并行计算成为可能[34-35].经过数十年的发展,LBM 的应用领域从单相问题延伸到多相流[36]、多孔介质[37-38]、在电流体动力学[9,12,30]等诸多领域,并受到了许多研究人员的广泛关注.在本文中,拟并采用4 个演化方程分别求解Navier−Stokes 方程,电势的泊松方程,电荷密度守恒方程和温度方程.接下来将依次介绍.

2.1 流场格子Boltzmann 模型

在假设流体为不可压牛顿流体的前提下,采用的是Guo 等[39]提出的格子Boltzmann 模型进行求解,该模型的演化方程可以表示如下

其中外力F=qE+gβ(θ −θref).τf为流场的无量纲松弛时间,它与 υ 的关系如下

最后计算宏观速度u以及压力项p,由以下公式计算得到

2.2 电势场格子Boltzmann 模型

本质上而言,电势控制方程(3)为一类典型的泊松方程.在这里采用下述演化方程[41]

其中,ζ 为扩散系数,取 ζ 等于0.5 以保持数值的稳定性.此外,源项R=q/ε.平衡态分布函数定义如下

因为电势的平衡态分布函数为线性的平衡态形式,为加快运算效率,采用D2Q5 的格子模型[41],该模型的权系数取值如下

而电场E则可以通过式(4)和电势的一阶分布函数计算得到[12],具体计算表达式如下

2.3 电荷密度的格子Boltzmann 模型

从文献[30]中可知电荷的传输机制有三种,第一种是在电场作用下,电荷会发生漂移;第二种是电荷在流体速度的影响下运动;第三种则是电荷的扩散.而电荷密度守恒方程属于强对流方程,对于该方程的求解在算法中起着重要作用,因此为了能很好的求解该方程,使用罗康[30]提出的格子Boltzmann模型来求解,其演化方程表示如下

上 述等式中的权重 ωi,离散速度ci都与流场相同.

2.4 温度场格子Boltzmann 模型

本文中的温度方程是忽略了黏性热耗散和焦耳热后的简化方程,从实现的角度和精度的满足两个方面,简化的温度方程能很好地应用于各种传热问题中[42].因此采用文献[12]中的方法,温度的演化方程如下

3 边界处理和程序流程

在格子Boltzmann 方法中,边界条件的处理一般采用Guo 等[44]提出的非平衡态外推格式,所得的分布函数精度为二阶,具有很好的数值稳定性.另外,程序的具体步骤为:(1)初始化为系统平衡态,此时电荷运输和温度传递处于平衡状态,流场处于静止状态;(2)进行碰撞流动,并通过式(19)计算速度u;(3)计算电势的演化方程(21),结合式(25)和式(26)得到电势 φ 和电场强度E;(4)计算电荷密度分布的演化方程(27),通过式(30)得到电荷密度q;(5)计算温度的演化方程(31),并计算式(34) 得到温度 θ ;(6) 计算式(17) 得到结果并代入流场的演化方程(13)中,依次重复步骤(2)～ 步骤(6)直到系统稳定并收敛为止.最后,定义二维问题中的局部努塞特数Nuy和平均努塞特数Nuav[34]以及最大的流体速度Vmax分别用来表示传热效率和流动强度,通过如下公式计算得到

4 数值验证

如图2 所示,对系统处于平衡态时的电荷密度q和x方向的电场强度Ex的数值解与解析解进行了比较.其中,平衡态电荷密度q和电场强度Ex的解析解表示如下

图2 电荷密度 q,水平方向电场强度 E x 的数值解与解析解对比Fig.2 Comparisons of the charge density q and horizontal electric field strength E x between the analytical solutions and numerical solutions

上式中的参数a和b的取值与电荷注入强度C有关[12],在验证中取C=10.

结果表明,通过当前LBM 代码计算得到的数值解 与解析解吻合较好.

5 数值模拟结果与讨论

5.1 库仑力和浮升力作用下的电热对流问题探究

本节中,对二维局部受热腔体结构中的电热对流问题进行了模拟,并初步对该结构下系统的传热效率、流体流动以及电热对流的分岔结构进行了分析.实验中将采用 3 00×300的网格分辨率,其可满足网格无关性验证.如图3(a)所示,首先观察了在热对流很弱时(热瑞利数Ra=400),电对流主导流体流动的情况,其他参数分别为C=10,M=10,Pr=10,h=0.5H,δ=0.5H.图3(a)中描述了随着电瑞利数T的增加和降低,相应传热效率Nuav的变化情况,图3(a) 中还给出了相关分岔点处的电荷密度分布图.对于电瑞利数T的值逐渐增加的情况,当电瑞利数T低于临界值[45](Tc≈174,Pr=10,M=10) 时,Nuav的变化不是很大,在前半段几乎保持水平状态,并且电荷密度分布比较均匀,这都表明流体没有形成强烈的对流,库仑力无法有效的激发流体流动,因此传热效率变化速度缓慢.然而当电瑞利数T超过临界值Tc以后,Nuav的值会突然向上跳跃式增加,并且电荷分布图中会出现空电荷区,空电荷区的形成与电荷运输过程中的漂移和对流机制之间的竞争有关[7,9,12,46].从电荷密度守恒方程(5)中的对流速度(KE+u)可以知道,一旦该区域的最大速度超过了电荷的漂移速度KE,强对流结构就会阻碍扩散到这些区域的电荷,结果就出现了空电荷区,并且会发生传热增强的现象.对于逐渐减少电瑞利数T的情况,当从电瑞利数T=300 开始逐渐减少时,Nuav并不能沿着增长时的曲线完全返回,这意味着强大的对流可以维持很大范围内的电瑞利数T值.但是当电瑞利数T接近148 附近时,Nuav会突然下降,可以观察到电瑞利数T=148 和T=174之间会形成所谓的回滞区,并且如预期中那样,空电荷区消失,Nuav会沿着曲线原路返回.接着观察在电对流很弱时(电瑞利数T=80),热对流主导流体流动的情况,图3(b)中随着Ra的增加和降低,往返的Nuav共用相同的演化曲线,表现出超临界特征(Pr=10,M=10).另外,Nuav变化过程中并没有出现回滞区,这说明流体流动表现出规律且单一的机制.在图3(b)中,还插入了热瑞利数Ra=800时的电荷密度分布和温度分布图,通过观察发现它们并没有形成空电荷区,电荷密度分布与温度分布相似,说明了流体运动十分弱,速度无法超过电荷的漂移速度.这也说明当电瑞利数T的取值很低时,只凭借浮升力的作用无法产生电热对流中分岔结构.通过对两种不同情况下的电热对流的稳定性分析,能够对电热对流中许多现象的产生有一定的初步认识,接下来会探讨在库仑力和浮升力共同作用下的电热对流问题.

图3 (a) 热瑞利数 Ra=400,变化电瑞利数 T 的分岔结构和(b) 电瑞利数T=80,变化热瑞利数Ra 的分岔结构 (C=10,M=10,Pr=10,h=0.5H,δ=0.5H)Fig.3 (a) The distributions of the maximum speed Vmax vs.electric Rayleigh number T .(b) The distributions of the maximum speed Vmax vs.Rayleigh number R a (C=10,M=10,P r=10,h=0.5H,δ=0.5H)

固定参数Ra=1×104,C=10,M=10,Pr=10,h=0.5H,δ=0.5H,图4 描绘了电瑞利数T在100～1200 范围内,不同强度的库仑力作用下的系统传热效率和流体流动的变化情况.可以看出,由于库仑力和浮升力之间的相互竞争,每条曲线都对应着不同的流动状态.随着电瑞利数T的增加,库仑力变大,Nuav和Vmax显著增加,系统的传热效率和流动强度都会变大.值得注意的是,当电瑞利数T较小(例如电瑞利数T=100)时,系统始终处于稳定的对流状态,此时的传热效率和流动强度的增强是有限的.接着当T增加时,则会在演化的最后出现一种不稳定的对流状态(T=200),并且Nuav会围绕某一值波动.但是当电瑞利数T的值增大时(5 00 ≤T≤800),系统再次出现稳定的对流状态(电瑞利数T=500)以及周期性流动状态(电瑞利数T=600,T=800时),主要原因是库仑力逐渐增大后,会在对流中占据主导地位.随着电瑞利数T的增加(当电瑞利数T=1200时),流动变得十分不稳定,此时的库仑力急剧增加并由稳定或周期性的对流状态变成混沌状.

图4 在不同电瑞利数T 下,(a) 平均努塞特数 N uav 和(b) 最大速度 Vmax 随时间变化的曲线 (R a=1×104,C=10,M=10,P r=10,h=0.5H,δ=0.5H)Fig.4 The distributions of (a) the average Nusselt number N uav vs.time and (b) the maximum speed Vmax vs.t ime under different electric Rayleigh number T (R a=1×104,C=10 ,M=10,P r=10,h=0.5H,δ=0.5H)

为了直观地理解电瑞利数T对电热对流的影响,如图5 所示,给出了稳定和非稳定情况下的电荷密度分布、温度分布和流线分布图.图5(a)～ 图5(e)分别为电瑞利数T=0,2 00,5 00,8 00,1 200时,其余无量纲参数为Ra=1×104,C=10,M=1 0,Pr=10,h=0.5H,δ=0.5H.可以看到,当电瑞利数T=0时,表示为纯自然对流的情况,方腔内只有一个沿顺时针方向旋转的漩涡,且上下对称.当电瑞利数T=200 时,流线分布图与电瑞利数T=0时相似,都出现了一个漩涡,虽然在方腔顶部出现了两个次级涡,但是库仑力还比较弱,对流动和传热产生的贡献很小,电荷只能在浮升力的影响下朝方腔顶部移动,下方的电荷密度则会变小.当电瑞利数T=500时,流线变得十分复杂,不仅由一个漩涡分裂成两个漩涡,而且方腔的4 个角落均出现了4 个小涡,此时的库仑力急剧增加,对流动和传热效率的影响大于浮升力,因此引导电荷逐渐摆脱浮升力的束缚,进而形成了上下对称的空电荷区,此时温度的分布受到了库仑力的影响,温度的热羽流与空电荷区形状相似.当电瑞利数T继续增加到800 和1200 时,流线中的上下两个漩涡都分别超过了分界线,流体流动的平衡再次被打破,由于强大的库伦力导致了对流变得十分不稳定,空电荷区变大,传热效率显而易见的增加.

图5 电瑞利数(a) T=0,(b) T=200,(c) T=500,(d) T=800,(e) T=1200 时的电荷密度分布(上)、温度分布(中)和流线分布(下)图(R a=1×104,C=10 ,M=10,P r=10,h=0.5H,δ=0.5H)Fig.5 The distributions of the charge density (top),temperature (middle) and streamlines (bottom) under different electrical Rayleigh number T.(a) T=0,(b) T=200,(c) T=500,(d) T=800,(e) T=1200 (R a=1×104,C=10,M=10,P r=10,h=0.5H,δ=0.5H)

最后,在图6(a) 中,给出了在热瑞利数Ra=1×104,5 ×104,1 ×105情况下的平均努塞特数Nuav随电瑞利数T变化的曲线,并固定其他无量纲参数C=10,M=10,Pr=10,h=0.5H,δ=0.5H.通过上述分析,知道电瑞利数T是存在一个临界阈值Tc,当电瑞利数T超过该值时,Nuav会急剧增加.在图6(a)中,可以发现类似的现象,热瑞利数Ra=1×105时,此时电瑞利数T的临界阈值在400～ 500 之间,电瑞利数T的取值超过500 以后,Nuav会发生明显的变化.

图6 (a) 探究不同电瑞利数T 下的平均努塞特数 N uav 的变化过程和(b) 探究不同热瑞利数 R a 下的平均努塞特数 N uav 的变化过程(C=10,M=10,P r=10,h=0.5H,δ=0.5H)Fig.6 (a) The distributions of the average Nusselt number N uav vs.electrical Rayleigh number T under different Rayleigh number R a .(b) The distributions of the average Nusselt number N uav vs.Rayleigh number R a under different electrical Rayleigh number T (C=10,M=10,P r=10,h=0.5H,δ=0.5H)

需要指出的是,热瑞利数Ra=5×104时,电瑞利数T在 3 00～400之间出现的凹点是因为此区间内的库仑力和浮升力对腔体下半区域的对流影响相当,加之二者所致的涡流方向在此区域又相反[27],如此二者相互竞争尤为激烈,腔体内的总体能量减少(主要体现在动能上),进而导致传热效率相对电瑞利数T≤300时略有降低,从而出现凹点.而当热瑞利数Ra=104,电瑞利数T=300～400时,因此时浮升力相对库仑力的影响较小,系统内库仑力影响绝对占优,因此方腔内的传热效率随着电瑞利数T的增加而增加,从而也不会出现上述凹点.当电瑞利数T小于临界阈值Tc时,浮升力的作用会增强系统的传热效率,但是当电瑞利数T大于临界阈值Tc以后,Nuav的上升是由于库仑力增强导致的.当电瑞利数T足够大时,因为库仑力在流动中发挥主导作用,反观浮升力的影响可以忽略,所以传热增强不再取决于热瑞利数Ra,可以看到3 条曲线最后收敛为一个点.在图6(b)中,探究了不同电瑞利数T值情况下,Nuav随着热瑞利数Ra的演化曲线,当电瑞利数T的取值大于500 以后,曲线的斜率是逐渐下降的,这是库仑力在主导热传递引起的,浮升力的影响在逐渐减弱,特别是当电瑞利数T=800时,曲线几乎是水平的,Nuav的变化与热瑞利数Ra的取值无关.可以注意到,与纯自然对流问题(电瑞利数T=0)情况相比,Nuav显著增加,当热瑞利数Ra=1×104时,在电瑞利数T的取值为500 的情况下,较于纯自然对流传热大约增强了 8 4.31%,说明加入电场后会强化传热.在本节中,探究了电热对流中的基本特征,并加深了对电热对流中库仑力和浮升力共同作用下流体流动结构、温度分布以及电荷的运输机制的理解.接下来,会对其 他无量纲参数进行研究.

4. B短语辨析。A.为……担心；B.为……感到自豪；C.害怕；D.为……感到遗憾。联系前一句描述，可知此处指的是我为你感到自豪。故选B。

5.2 电极板位置影响探究

本节中,研究了电极板位置的影响.图7 描绘了电极板位置对在不同电瑞利数T下的传热效率的影响,其余各无量纲参数为Ra=1×104,C=10,Pr=10,M=10,h=0.5H.其中,当电瑞利数T=800,热瑞利数Ra=1×104时,相比于浮升力,库仑力影响绝对占优;将加热壁面从下半区域 (δ=0.25H)更改为上半区域 (δ=0.75H)后,浮升力变化带来的影响基本可以忽略,此时系统传热效率主要由库仑力决定,又因两种情况下流动结构以及库仑力分布等基本镜面对称,从而T=800时,上下两种位置对应的传热效率基本重合.另外,当加热壁面处于中间位置(δ=0.5H) 时,上下两部分的速度分布相对更加均匀,流体流动速度比另外两个位置下的要大,对流强度相对更强,温度边界层也相对更薄,因而此位置下的传热效率一般优于另外两种情况.另外,在 0 ≤T≤225 区间,电极板在上半区域时(δ=0.75H)的Nuav小于在下半区域时(δ=0.25H)的Nuav.在图7中画出了电瑞利数T=200 时,A和B两点处的平均努塞特数变化曲线图,可以看出B点处的对流十分不稳定,这也说明库仑力和浮升力竞争相较于A点处更加激烈,并在极大程度上影响了系统的传热效率,因此 δ=0.75H时的传热效率低于 δ=0.25H时的传热效率,但是随着电瑞利数T的取值继续增加,库仑力在对流中的影响增强,当电瑞利数T=800时,库仑力占据绝对优势,此时系统的传热效率与热瑞利数Ra无关,所以两个位置重合.最后,在电瑞利数T=200 到T=300之间,3 条曲线分别出现了一个临界阈值,这与图6(a)中热瑞利数Ra=1×104时出现的临界阈值范围相同.

图7 电极板在不同位置下,随着电瑞利数 T 增加,平均努塞特数Nuav 的变化曲线 (R a=1×104,C=10,P r=10,M=10,h=0.5H)Fig.7 The distributions of the average Nusselt number N uav vs.electrical Rayleigh number T under different heating locations δ (R a=1×104,C=10,P r=10,M=10,h=0.5H)

如图8 所示,画出了当电瑞利数T=250时,δ=0.25H和 δ=0.75H时的电荷密度分布、温度分布和流线分布图.在这两种情况下,电荷密度分布图中都出现了空电荷区,根据在5.1 节中的分析,可以知道,此时的库仑力增大并与浮升力相互竞争,通过观察和分析两种情况下的流线分布(图8(c1)和图8(c2)),当电极板在下半区域时,由于库仑力增强,导致了下部分对流减弱,并且图8(c1)下部分的漩涡比图8(c2)的漩涡小,因此系统的传热效率也降低了.所以从图7 中可以发现电极板在上半区域的Nuav高于在下半区域时的Nuav,并且在之后的发展中一直保持(即使差异会越来越小).

图8 电极板在上半区域和下半区域时的电荷密度分布图(上)、温度分布图(中)和流线分布图(下) (T=250,R a=1×104,C=10,Pr=10,M=10,h=0.5H)Fig.8 The distributions of the charge density (top),temperature(middle) and streamlines (bottom) as the electrode plate is at the top and bottom positions (T=250,R a=1×104,C=10 ,P r=10,M=10,h=0.5H)

图9 给出了电瑞利数T=800时,电极板在下半区域图9(a),中间图9(b),上半区域图9(c)位置时的电荷密度分布、温度分布和流线分布图.可以看到在电荷分布图中都形成了一个很明显的空电荷区,此时的库仑力在流动过程中占据主导地位.当电极板在上半区域和下半区域时,库仑力分布不均匀,从而导致速度大小分布不均匀.电极板在下半区域时,下半部分的速度大小会比上半部分的小,因此方腔下半部分的热量转移相比上半部分的热量转移要少,从温度分布图中的热羽的形状也可以发现靠近壁面的温度边界层变厚.无论是电极板在上半区域还是在下半区域,与电极板在中间时相比,都会使得靠近壁面一侧的温度梯度降低,从而传热效率变低.通过以上分析说明了当电极板放在中间位置时,传热效率最佳.另外,电极板在上半区域和下半区域时的漩涡结构类似,这也说明了为什么传热效率会随着电瑞利数T的取值增加而最终接近.通过对比图8 和图9 中的电荷密度分布,可以发现电荷分布会受到库仑力和浮升力的影响,当电瑞利数T的取值很小时,库仑力较弱,左壁面的电荷会由于浮升力作用向上壁面运动,而在右壁面一定范围内的电荷向下运动后,接着会沿着下壁面向左壁面流动,从而形成一个顺时针的对流.虽然库仑力始终向右并会对下部分一定区域的对流造成阻碍,但是电瑞利数T的取值很小,浮升力大于库仑力,因此流动的方向仍然和在浮升力影响下的方向一致.但是随着电瑞利数T的取值增加,库仑力逐渐增加,流动的方向逐渐受到库仑力的影响,最终的对流方向会被库仑力主导,这也可以从图9 中的温度分布图中的热羽形状和电荷分布形状相似看出.

图9 电极板在(a)下半区域,(b)中间,(c)上半区域,电荷密度分布(左)、温度分布(中)和流线分布(右)图 (T=800,R a=1×104,C=10,Pr=10,M=10,h=0.5H)Fig.9 The distributions of the charge density (left),temperature (middle) and streamlines (right) as the electrode plate is at the top,middle and bottom positions (T=800,R a=1×104,C=10,P r=10,M=10,h=0.5H)

5.3 电极板长度影响探究

本节中,考虑了电极板长度h对传热的影响,相应的数据结果如图10 所示,其余各无量纲参数为Ra=1×104,C=10,Pr=10,M=10,δ=0.5H.从图10中可以看出,对于给定的电瑞利数T,传热效率会随着电极板的长度增加而降低,并且3 种情况下电瑞利数T的临界阈值都在 2 00～300之间,在大于临界阈值以后,它们的变化趋势十分相似,再结合上一节中的临界阈值,发现电极板长度h和位置 δ 对临界阈值的范围影响不大.

图10 电极板长度不同时,随着电瑞利数 T 增加,平均努塞特数Nuav 的变化曲线 (R a=1×104,C=10,P r=10,M=10,δ=0.5H)Fig.10 The distributions of the average Nusselt number N uav vs.electrical Rayleigh number T under different length of electrode plate h(R a=1×104,C=10 ,P r=10,M=10,δ=0.5H)

分别画出了电瑞利数T=250(图11)和电瑞利数T=300(图12)时的电荷密度分布、温度分布和流线分布图.在图11 中,随着电极板的长度增加,温度分布图中的热羽流变小至一半并类似于纯自然对流(电瑞利数T=0),空电荷区向下移动并逐渐远离电极板,从流线分布图中可以看出下方的一级漩涡最终退化为次级漩涡,下方的对流由有序逐渐变为无序,从而整体的传热效率逐渐下降.又如图12 所示,当电瑞利数T=300时,此时的库仑力在与浮升力的竞争中占据了优势,并且在方腔的中心区域形成了空电荷区.但是随着电极板两端的增加,电场力减弱,热壁面的速度下降,无法将更多的热量转移到介电液体中,同时热壁面的温度边界层变厚,温度梯度降低,因此传热效率会随着电极板的长度增加而降低.

图11 不同长度下,(a) h=0.5H,(b) h=0.75H,(c) h=H,电荷密度分布图(左)、温度分布图(中)和流线分布图(右) (T=250,R a=1×104,C=10 ,P r=10,M=10,δ=0.5H)Fig.11 The distributions of the charge density (left),temperature (middle) and streamlines (right) under different length of electrode plate.(a) h=0.5H,(b) h=0.75H,(c) h=H (T=250,R a=1×104,C=10,P r=10,M=10,δ=0.5H)

图12 不同长度下,(a) h=0.5H,(b) h=0.75H,(c) h=H,电荷密度分布图(左)、温度分布图(中)和流线分布图(右) (T=300,R a=1×104,C=10 ,P r=10,M=10,δ=0.5H)Fig.12 The distributions of the charge density (left),temperature (middle) and streamlines (right) under different length of electrode plate.(a) h=0.5H,(b) h=0.75H,(c) h=H (T=300,R a =1×104,C=10,P r=10,M=10,δ=0.5H)

5.4 无量纲离子迁移率和普朗特数的联合影响

如图13 所示,还研究了无量纲离子迁移率M和普朗特数Pr的联合影响,其余无量纲参数为Ra=1×104,T=300,C=10,δ=0.5H,h=0.5H.从图13 可以看出在给定Pr的情况下,随着M的增加,平均努塞特数Nuav会从一开始的急剧下降转变为趋势平缓的情况.M越小,Nuav会随着Pr增加而增加,当M很大时(M=50),Nuav几乎保持不变.事实上,如文献[9]指出,M的大小决定了静电方程组与Navier-Stokes 方程组的耦合紧密程度,M越大,两个方程组的耦合程度越弱,进而库仑力对流场的作

用减弱.另外,不同的M和Pr数的流体流动状态也存在一定差异,主要分为稳态和非稳态,具体可见图13中的标记.

图13 固定普朗特数Pr,在不同无量纲离子迁移率M 下,平均努塞特数 N uav 的变化曲线 (R a=1×104,T=300,C=10,δ=0.5H,h=0.5H)Fig.13 The distributions of the average Nusselt number N uav vs.dimensionless ion mobility M under different Prandtl number P r (R a=1×104,T=300,C=10,δ=0.5H,h=0.5H)

最后,取图13 中Pr=50的曲线肘部处和端点处的值即M=25,M=30和M=50,画出这3 点处的电荷密度分布、温度分布和流线分布图,如图14所示,当M=25和M=30时,空电荷区会变得不明显,还会偏移中心向下且远离电极板,这都说明库仑力的影响在减弱,而且流线分布图中的漩涡结构的差异更加明显,由两个反向旋转的漩涡变成了一个漩涡结构,这不仅减弱了对流而且破坏了有序的对流结构,因此降低了系统的传热效率.当M=50时,差异会变得十分明显,虽然离子迁移率增加导致大量电荷注入到方腔中,但是两方程组的耦合程度减弱,库仑力减弱导致系统的对流减弱.所以当M越大时,系统的传热效率反而会减小.

图14 不同M 值下的电荷密度分布(左)、温度分布(中)和流线分布(右)图.(a) M=25,(b) M=30,(c) M=50 (P r=50,R a=1×104,T=300,C=10,δ=0.5H,h=0.5H)Fig.14 The distributions of the charge density (left),temperature (middle) and streamlines (right) under different dimensionless ion mobility M .(a) M=25,(b) M=30,(c) M=50 (P r=50,R a=1×104,T=300,C=10,δ=0.5H,h=0.5H)

6 结论

本文采用了格子Boltzmann 方法对二维局部受热腔体结构中电热对流问题进行模拟,首先对系统的传热效率、流体流动以及电热对流的分岔结构进行了分析,对在库仑力和浮升力共同作用下流体的流动、温度的分布以及电荷的运输机制有了初步理解,然后探究了电极板的位置和长度的影响,最后探讨了无量纲离子迁移率M和普朗特数Pr的联合影响.根据以上研究中呈现的数值结果,主要结论如下:

(1) 固定热瑞利数Ra,平均努塞特数Nuav会随着电瑞利数T的增加而增加,这主要是由于库仑力增强导致的.当电瑞利数T超过临界阈值(Tc,Pr=10,M=10)后,Nuav会急剧增加,并且热瑞利数Ra越大,临界阈值也会增加.通过稳定性分析可以发现,电瑞利数T会存在一个亚临界点Tf,然而热瑞利数Ra不会出现分岔结构,会表现为超临界特征.另外,当电瑞利数T足够大时,库仑力占据绝对优势,浮升力对传热的影响减弱,Nuav的值不再依赖于热瑞利数Ra.

(2) 传热效率很大程度上受到电极板位置的影响,并且在中间位置时的传热效率最佳,这主要是由于电极板在中间位置时的对流最强导致的.

(3) 传热效率会随着电极板的长度增加而降低,这主要是由于电极板长度增加以后,方腔底部的对流减弱,从而导致传热效率降低.

(4) 当M越大时,系统的传热效率越低,这是因为随着M的增加,库仑力会对流场的影响减弱,从而导致对流减弱和传热效率降低.