环形容器内双扩散对流的数值模拟_石科峰

第22卷第2期2005年3月中国科学院研究生院学报Journal of the Graduate School of the Chinese Academy of Sciences Vol .22March No .22005＊国家自然科学基金项目(50336040)资助文章编号:1002-1175(2005)02-0152-05环形容器内双扩散对流的数值模拟＊石科峰　卢文强(中国科学院研究生院物理科学学院,北京100049)(2004年6月9日收稿;2004年7月9日收修改稿)Shi KF ,Lu WQ .Double diffusive convection in an annular vertical cavity .Journal of the G raduate School of the Chinese Academy of Sciences ,2005,22(2):152～156摘　要　通过数值模拟,研究了温度和浓度梯度为水平方向、内壁和外壁上的温度和浓度保持定值的竖直环形容器内的双扩散对流结构,其中浮升力之比N =G r S G r T =-1,G r S 和G r T 分别为溶质格拉晓夫数和热格拉晓夫数.首先对环形容器与矩形容器内的双扩散对流状况进行了比较;然后着重分析了在双扩散对流情况下竖直环形容器内壁面处的热(溶质)边界层以及平均传热(质)率Nu (Sh )随Le 数的变化情况.结果表明,溶质边界层的厚度随着Le 数的增大逐渐减小,而热边界层随Le 数的变化不大;与Sh 数一直随Le 数的增大呈明显的上升趋势不同,当Le >10后Nu 数随Le 数的变化趋于平缓.关键词　双扩散对流,晶体生长,边界层,有限单元法中图分类号　O5511　引言双扩散对流是由于流体当中热量和溶质具有不同扩散率而产生的一种复杂的对流传热传质过程,这种物理现象存在于很多领域当中,如海洋学、地质学以及晶体生长[1].在多组元晶体生长过程当中,液相区内的热量和溶质的传输是通过扩散和对流来完成的.生长界面的推移所产生的温度梯度和浓度梯度决定着热量和溶质的扩散传输过程;流体中由于温度梯度和浓度梯度的存在而产生的浮升力则决定着热量和溶质的对流传输过程[2].这两种传输过程联合作用下流体的流动状况就是一种双扩散对流,它决定着相变界面处液相区的温度和浓度分布,是能否获得高品质的单晶的关键.因此研究双扩散对流,能够为以后指导制备高品质的单晶奠定基础.目前,一些国外学者已经开展了双扩散对流的研究,但大多数文献均是针对矩形容器内的双扩散对流,如文献[3～7]数值模拟了不同边界条件下矩形容器内的双扩散对流结构及其波动情况,文献[8～10]则通过实验研究分析了矩形容器内多组元混合物凝固过程中的双扩散对流现象.与此同时,针对环形和圆柱形容器内的双扩散对流的研究却比较少,M .Marcoux [11]仅仅对竖直环形多孔介质中的双扩散对流进行了探讨.本文通过有限单元法求解圆柱坐标系下的涡量、流函数、能量和组分控制方程,数值模拟了竖直环形容器内的双扩散对流结构,着重探讨和分析了竖直环形容器内壁面处热(溶质)边界层和平均传热(质)率随Le 数的变化情况.图1　环形容器示意图与坐标系2　控制方程及数值方法2.1　控制方程我们选取的物理模型和坐标系如图1所示.环形容器的内外壁面的半径分别为r 1和r 2,径向厚度L =r 2-r 1,高度为H ,宽高比A =L H .假定系统内部不存在内部热源、黏性耗散、化学反应以及热辐射现象,并且B oussinesq 假设成立,流体密度ρ(T ,C )=ρ0(1-βT (T -T 0)-βS (C -C 0)),其中ρ0=ρ(T 0,C 0)为温度T 0和浓度C 0时的密度,βT =(-1 ρ0)( ρ T )C 和βS =(1- ρ0)( ρ C )T 为热膨胀系数和溶质膨胀系数,其它热物性为常数.同时假定βT >0及βC <0(即密度值与温度值成反比,与浓度值成正比).在自然对流传热传质中,压力常常不是主要的求解变量,因此采用涡量流函数方程特别合适.流函数和涡量的定义如下:u r =1r ψ z ,u z =-1r ψ r ,ω= u r z - u z r.(1) 得到无量纲化的涡量方程、流函数方程、能量方程及组分方程如下: ω t +1r ψ z ω r -1r ψ r ω z -ωr 2 ψ z = 2ψ r 2+1r ω r + 2ω z 2-ωr2-G r T T r -Gr S C r , 2ψ r 2+ 2ψ z2-1r ψ r =r ω, T t +1r ψ z T r -1r ψ r T z =1Pr 2T r 2+1r T r + 2T z2, C t +1r ψ z C r -1r ψ r C z =1PrLe 2C r 2+1r C r + 2C z2.(2)其中,热格拉晓夫数G r T =g βT ΔTH 3υ2;溶质格拉晓夫数G r S =g βS ΔCH 3υ2;普朗特数Pr =υ α;刘易斯数Le =α D ;D 、α、υ、g 分别为溶质扩散系数、热扩散系数、黏性系数和重力加速度.在上述控制方程中,采用的无量纲化标准为:尺度用H ;速度用υ H ;时间用H 2 υ;温度用ΔT =T 2-T 1;浓度用ΔC =C 2-C 1;涡量用υ H 2;流函数用H υ.相应的边界条件如下:H =0∶ T z = C z =0,ψ=0,H =1∶ T z = C z=0,ψ=0,r =r 1=0.5∶T =C =1,ψ=0,r =r 2=1.5∶T =C =0,ψ=0.(3) 通过无量纲化的努谢尔特数(Nusselt number )和舍伍德数(Sher wood number )可以定义内外壁面上的平均传热和传质率:Nu =∫10 T r r =r 1,r 2d z ,Sh =∫10 C r r =r 1,r 2d z .(4)2.2　数值方法本文采用有限单元法来求解上述控制方程.计算单元为4节点的矩形单元,在单元内构造双线性的153　第2期石科峰,卢文强:环形容器内双扩散对流的数值模拟插值函数=b1+b2ξ+b3η+b4ξη.在求解控制方程的过程中,首先求解能量方程和组分方程,得到计算节点上的温度梯度和浓度梯度,然后将其带入到涡量方程当中求解涡量值,最后通过涡量求解流函数方程得到流函数,由流函数可以求得各个节点处的速度值,如此迭代.其中在求解涡量方程和流函数方程的过程中,需要引入一个确定壁面涡量的内迭代过程.迭代过程均采用超松弛(SOR)迭代,其中壁面涡量的确定采用Thom公式[12],圆柱坐标系下的Thom公式经推导为:ω1=2(ψ2-ψ1)rδ2.(5)2.3　计算程序的验证为了验证计算方法的正确性,我们采用扩大内径的方法来模拟矩形容器内的双扩散对流结构.因为当容器厚度与内径相比变得很小时,即r1L1,这时内外壁面的面积已经相当,容器的圆弧效应大大削弱,就可以用来模拟截面前后均趋于无穷的矩形容器内部的双扩散对流状况.如图2(a)所示.取r1= 100时计算所得的流函数、温度和浓度场(Pr=1、Le=21、G r T=6.5×104)与K.Ghorayeb[5]关于矩形容器内双扩散对流的结构为中心对称的描述一致;同时我们得到的容器竖直壁面上Nu和Sh的值分别为3.634和11.305,K.Ghora yeb的结果为3.664和10.77,二者符合得较好.图2　矩形容器(a)和环形容器(b)内的流函数线、温度场与浓度场模拟(Pr=1,Le=21,Gr T=6.5×104)表1　计算中心点处的流函数值、温度值与浓度值随网格划分的变化(Le=20,G rT=1×104)网格划分流函数ψ温度T浓度C21×21-1.36280.27780.266851×51-1.31160.28740.307871×71-1.30700.28730.30942.4　计算网格的确定为了确定计算网格,我们分别取了21×21、51×51、71×71的网格划分.表1显示了容器中心处的的流函数、温度和浓度计算值与网格节点数的关系.可以看出51×51网格划分下的解与71×71网格划分下的解相差不超过1%,同时考虑到求解所消耗的时间,我们选取了51×51作为基本网格划分.3　结果与分析在双扩散对流问题当中,容器内部流体的流动状况、传热量和传质量取决于以下5个参数,它们是宽高比A、G r T、P r、Le以及浮升力比N=G r S G r T.在这些参数当中,以Le以及N对双扩散对流的状况影响最大[4].试验和数值分析均表明,在Le数较大的情况下多组分流体在壁面处存在一个很薄的溶质边界层,但它并不存在于整个壁面边界处.在下面的研究分析中,我们给定竖直环形容器的内径r1=0.5, A=1.0,N=-1.0,首先模拟了环形容器内的双扩散对流结构,并对比分析了它与矩形容器内的双扩散对流结构的差别;随后我们分别对环形容器内壁面处的热边界层和溶质边界层随Le数的变化关系,以及内壁面处的传热(质)率Nu(Sh)随Le数的变化关系进行了分析.3.1　环形容器内的双扩散对流结构图2(b)给出了Pr=1.0,G r T=6.5×104,Le=21时,竖直环形容器内部相应的流函数场、温度场及浓度场.通过和图2(a)所示的矩形内相应的场相比较,可以看到,环形容器内双扩散对流的结构已经不再具有如矩形容器内的那样的中心对称结构,并且靠近内壁面处的温度及浓度的径向梯度比外壁面处的明显要大;同时,内壁面处的等流函数线所覆盖的高度,以及热边界层和溶质边界层所涉及的高度比外壁面处的要大.这主要是由于内壁面的面积小于外壁面造成的.154中国科学院研究生院学报第22卷3.2　内壁面处热边界层和溶质边界层随Le 数的变化情况图3　不同Le 数情况下的温度边界层(a )和溶质边界层(b )(Pr =7.0,G r T =1×104)图3表明了Pr =7.0,G r T =1×104条件下内壁面处的热边界层(a )和溶质边界层(b )随Le 数的变化关系(为了研究方便,图3中是对r =0.5～0.7、z=0.0～0.4区域的局部放大).从图中可以看出,热边界层随Le 数的变化并不是特别明显;而溶质边界层的厚度明显随着Le 数的增大而逐渐减小.在其他参数不变的情况下,Le 数的增大表明溶质扩散系数D 的减小,因此从壁面处开始的溶质扩散波及的区域变小,边界层变薄.在较大Le 数(Le >2)的情况下,热扩散系数比溶质扩散系数大,因而溶质边界层比热边界层薄,相对较薄的溶质边界层的变化并不会对热边界层产生明显的影响.3.3　内壁面处的平均传热率(Nu )和传质率(Sh )随Le 数的变化情况图4(a )显示内壁面处的平均传热率(Nu )在Le<10的情况下随Le 数的变化比较明显,但当Le >10以后的Nu 变化就趋于平缓;图4(b )则显示内壁面处平均传质率(Sh )随着Le 数的增大呈稳定的增长趋势.我们可以从边界层厚度变化的角度来解释产生图4(a )和图4(b )两种不同变化规律的原因.当Le <10时,溶质边界层的变化对靠近壁面处流体的流动状况影响相当明显,从图5可以看到,Le <10时内壁面处竖直方向上的流速u z 随Le 数变化相当显著,图4　环形容器内壁面处平均传热率(Nu )和传质率(Sh )随Le 数的变化关系(Pr =7.0,Gr T =1×104)图5　z =0.2处的流速u z 分布随Le 数的变化关系(Pr =7.0,G r T =1×104)对流作用得到加强,热边界层的变薄导致了Nu 数的增大;当Le >10以后,溶质边界层已经变得相当薄,在它的外边浓度梯度很小,由温度梯度产生的浮升力占着主导作用,这时候溶质边界层的变化已经对壁面处的流体流动状况影响很小,如图5所示.因此,在热边界层变化不大的情况下,Nu 数的变化也相当平缓;与此相反,当Le 数增大时,随着溶质边界层的厚度逐渐变薄,壁面处的浓度梯度逐渐变大,于是Sh 数则呈一直上升的趋势.4　结论本文利用有限单元法求解通过涡量-流函数法建立的控制方程,数值模拟了具有径向温度梯度和浓度梯度,且浮升力方向相反的竖直环形容器内的双扩散对流的结构.分析表明:(1)与矩形容器相比,环形容器内的双扩散对流不再具有中心对称的结构,内壁面处的边界层效应要比外壁面处的明显;(2)溶质边界层的厚度随Le 数的增大逐渐变薄,但热边界层的变化并不显著;(3)当Le <10时,随着Le 数的155　第2期石科峰,卢文强:环形容器内双扩散对流的数值模拟156中国科学院研究生院学报第22卷增大Nu数上升得比较显著,Le>10后则趋于平缓,但Sh数则随Le数的增大一直呈稳定上升的趋势.References[1]　Turner JS.Double-diffusive phenomena.Annu.Re v.Fluid Me ch.,1974,6:37～56[2]　Zhang KC,Zhang LH.Science and Technology of Crystal Growth.(Second Edition).Beijing:Science Press,1997.401～442(in Chinese)[3]　Murty VD.A finite element solution of double diffus ive convection.Int.C omm.He at M as s Tr ansfe r,1988,15(2):165～177[4]　Nishimura T,Waka matsu M,Morega A M.Oscillatory double-diffusive convection in a rectangular encl os ure with combined horiz ontal temperatureand c oncentration gradients.Inter nati onal Jour nal of Heat and Mass Tr ansfe r,1998,41(11):1601～1611[5]　Ghorayeb K,Khall ouf H,Mojtabi A.Ons et of oscill atory flo ws in double-diffusive convection.Inter nati onal Jour nal of Heat and Mas s T rans fer,1999,42(4):629～643[6]　Sundaravadivelu K,Kandaswa my P.D ouble di ffusive nonlinear convection in a square cavity.Fluid Dynamics Res earc h,2000,27(5):291～303[7]　Bennacer R,Mohamad AA,Akrour D.Transient natural convection in an enclosure with horiz ontal temperature and vertical sol utal gradients.Int.J.The rm.Sci.,2001,40(10):899～910[8]　Nishimura T,Imoto T,M i yashita H.Occurrence and development of double-diffus ive convection during s olidification of binary s yste m.Inte rnationalJour nal of Heat and Mass Tr ansfe r,1994,37(10):1455～1464[9]　Chen FL.Formation of double-diffusive la yer in the directional solidification of binary s olution.Journal of C rys tal Gr owth,1997,179(1-2):277～286[10]　Ghenai C,Mudunuri A,Lin CX,et al.D ouble-diffusive convection during solidificati on of a metal analog s ystem(N HCL-H2O)in a differentially4 heated cavity.Experi me ntal T her mal and Flui d Sc ienc e,2003,28(1):23～35[11]　Marcoux M,Charrier-Mojtabi MC,Azaiez M.Double-diffus ive convection in an annul ar vertical porous layer.Inter nati onal Journal of Heat andMas s Trans fer,1999,42(13):2313～2325[12]　Tao WQ.Numerical Heat Trans fer.Second Edition.Xi'an:X i'an Jiaotong Univers ity Press,2001.309～314(in Chines e)附中文参考文献[2]　张克从,张乐主编.晶体生长科学与技术.上册.第二版,北京:科学出版社,1997.401～442[12]　陶文铨编著.数值传热学.第二版.西安:西安交通大学出版社,2001.309～314Dou ble Diffusive Convection in an Annular Vertical CavitySHI Ke-Feng　LU W en-Qiang(C ollege of Physic s Science,G raduate School,C hine se Academy of Sci ence s,Beijing100049,China)A bstract　A numerical study is presented of the structure of double diffusive convection in an annual cavity with horizontal heat and mass gradients due to constant temperatures and concentrations along the vertical side-walls, respectively.The buoyancy ratio N=G r S G r T is equal to-1,where G r S and G r T are the solutal and thermal Grashof numbers,respectively.The str uctures of double diffusive c onvection in rectangular and annular cavity were compared.Numerical results are reported for the effect of the Le numbers on the thermal(solute)boundary layer and the avera ge rate of heat(mass)transfer Nu(Sh).The results show that the thickness of solutal boundar y layer decreases when the Lewis Number increases,but the thermal boundar y is not significantly influenced when the Le wis number varies;contrasting with the steady increase of the Sher wood number,the Nusselt number varies slo wly when the Lewis number is higher than10.Key words　double diffusive convection,crystal growth,boundary layer,finite ele ment method。