数值传热大作业

1放置竖直孤立平板的二维围场内的空气流动与换热的数值分析（西安交通大学能源与动力工程学院，710049，西安）摘要：针对内部放置孤立平板的二维围场内的空气流动与换热问题，在稳态、常物性和壁面温度以及孤立平板温度恒定的条件下，采用SIMPLER算法，对围场内部的空气进行了流动与换热的数值模拟计算。在瑞利数Ra=10000时，计算得到了二维围场内的流线、等温线以及热线。关键词：SIMPLER算法、孤立平板、流线、等温线、热线Abstract:Inordertoinvestigatethefluxionandheattransferofairina2Dsquareenclosurewithanisolatedplate.SIMPLERalgorithmwasadoptedbasedontheReylondsconservationequationsofthesteady-stateconstantpropertylaminarflowandaconstanttemperatureoftheisolatedplateandtheinnerwallsoftheenclosurecondition.SlovefluidvelocityandtemperaturefieldsintheenclosureforRa=10000,anddrawthediagramsofstreamlines,isothermsandheatlines.Keywords:SIMPLERalgorithm;isolatedplate;streamlines;isotherms;heatlines2主要符号表aR瑞利数rP普朗特数空气运动粘度m2/sg重力加速度kg.m/s2k空气导热系数W/(m℃)空气体膨胀系数1/℃pc空气比热容J/kg.℃空气密度kg/(m3s)hT金属板温度℃cT围场壁面温度℃T温差℃3一、引言封闭空腔内孤立物体自然对流换热是一个重要的研究课题，从某种角度讲，大空间自然对流是封闭腔内孤立物体自然对流的一个特例。文献[1,2]应用原始变量法求解了圆形空腔内存在同心的孤立圆弧状或八边形物体时的自然对流换热问题，文献[3]研究了封闭腔内竖直孤立物体的自然对流。作为这一类问题的继续，本文着重研究了二维围场内竖直放置的孤立平板的自然对流，旨在计算分析得出在瑞利数等于10000时，二维围场内的流线、等温线以及热线。二、物理问题与数值描写本文研究的二维围场与孤立平板的布置如图1所示，正方形围场的壁面恒定温度为280cT℃，边长17Hmm，内部竖直放置的孤立平板温度恒定为180hT℃。围场内其余地方填充着空气。图1、二维围场平面示意图数值计算是在如下基本假设的前提下进行的：1、在整个计算过程中假设空气的物性是恒定不变的；2、空气的粘性耗散忽略不计；3、围场外表面假设是绝热的，与外界没有热交换；4、空气流动以及换热处于稳态、层流状态下。所研究的问题的控制方程如下：4连续性方程：0uvxy(1.1)动量方程：X方向：2222()uvpuuuvxyxxy(1.2)Y方向：2222()vvpvvuvxyyxy(1.3)能量方程：2222()TTTTuvkxyxy(1.4)其中边界条件：围场壁面：(0,)(,0)(,)(,)cTyTxTHyTxHT(1.5)(0,)(,0)(,)(,)0uyuxuHyuxH(1.6)(0,)(,0)(,)(,)0vyvxvHyvxH(1.7)孤立平板：温度恒定hTT四周边界速度0uv(1.8)三、数值求解方法此次数值模拟计算采用SIMPLER算法，并用区域离散方法B以及交错网格。为保证迭代计算过程中孤立平板处的速度恒为零（或与主流区域相比要小若干数量级的小数），采用了以下方法[4]：1、在每一层次的迭代计算前令孤立平板的速度为零，以保证孤立平板中的节点对流体区中的节点速度起滞止的影响；2、在求解速度的代数方程前令孤立平板各速度离散方程主对角元的系数为一很大值（此次计算过程中，该值取为1030），以保证预估值**0uv；3、计算压力修正值时，使孤立平板各速度修正值计算公式的系数（即,endd等）取一个近似于零的值，如10-38，以使孤立平板中各速度修正值也为零。四、计算结果与讨论在瑞利数10000aR条件下，计算得出速度场和温度场，并利用tecplot对计算数据进行处理得出流线、等温线、热线图如下：5图2、流线从流线图我们可以看出，在孤立平板的左右两侧的流场中出现了许多封闭的流线，说明空气在上述区域形成了自然循环流动。图3、等温线6对温度场中的等温线进行观察，我们可以很容易的得出在孤立平板的四个角处的等温线比较密集，说明四个顶角处的温度梯度比较大，而且围场壁面处的等温线也比较密集，温度梯度也比较大。图4、热线对热线进行分析，我们可以粗略的看出热流密度的方向是从二维围场墙壁指向孤立平板，这和围场墙壁温度（280℃）大于孤立平板的温度（180℃）的情况是相符合的，也验证了热量是从高温物体向低温物体单向传递的事实。在孤立平板的右侧，由于空间足够大，热线形成了一个自我闭合的同心曲线簇。五、结论1、靠近围板壁面和孤立平板的区域，由于温差比较大，对流效果比较明显，导致速度梯度比较大，流线密集。2、围板壁面和孤立平板四个顶角等处的温度梯度比较大，等温线比较密集；3、热流密度的方向是由高温物体指向低温物体，即热量传递方向是由高温物体向低温物体单向传递的。。参考文献[1]杨沫、陶文铨、陈钟欣，“缝隙度对封闭圆内开缝圆环自然对流换热的影响，”工程物理学报，Vol.12,No.2,1991,pp.193-196.[2]杨沫、陶文铨、陈钟欣，“开缝八边形离相封闭母线自然对流换热的数值计算，”中国电机工程学报，Vol.11,No.4,1991,pp.58-64.[3]王秋旺、王育清、陶文铨、杨沫，“几何位置对封闭方腔内水平孤立平板自然对流换热的影响，”工程物理学报，Vol.15,No.2,1994,pp.195-199.[4]陶文铨编著，数值传热学（第二版），2001，西安交通大学出版社，244-244.1附录源程序CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCMODULEUSER_LC******************************************************************INTEGER*4I,J,M00,M0,N0,N00REAL*8TIN,TW,VIN,VOUT,AMU,COND,BETA,TREF,FLOWIN,FL,GRA,TMC******************************************************************ENDMODULECCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCSUBROUTINEUSERC*******************************************************************USESTART_LUSEUSER_LIMPLICITNONEC******************************************************************C-----------------------PROBLEMSIX-----------------------CLaminarfluidflowandheattransferinatwo-dimensionalC---------------------suddenexpansion--------------------C******************************************************************CENTRYGRIDTITLE(1)='VEL_U'TITLE(2)='VEL_V'TITLE(4)='TEMP.'TITLE(11)='PRESSURE'TITLE(8)='EX'TITLE(9)='EY'RELAX(1)=0.8RELAX(2)=0.8LSOLVE(1)=.TRUE.LSOLVE(4)=.TRUE.LPRINT(1)=.TRUE.LPRINT(2)=.TRUE.LPRINT(4)=.TRUE.LPRINT(11)=.TRUE.2LPRINT(8)=.TRUE.LPRINT(9)=.TRUE.LAST=200XL=0.017YL=0.017L1=63M1=63N0=20N00=28M0=18M00=44CALLUGRIDRETURNCENTRYSTARTTREF=20DO100J=1,M1DO101I=1,L1U(I,J)=0.0V(I,J)=0.0T(I,J)=0.0T(1,J)=280.0T(L1,J)=280.0T(I,1)=280.0T(I,M1)=280.0U(1,J)=0.0U(L1,J)=0.0V(I,1)=0.0V(I,M1)=0.0101ENDDO100ENDDODO120J=M0,M00DO121I=N0,N00T(I,J)=180.0U(I,J)=0.0V(I,J)=0.0121ENDDO120ENDDOAMU=3.472E-5CPCON=1034.0COND=3.931E-2GRA=9.8BETA=3.676E-33RHO=0.746RETURNCENTRYDENSERETURNCENTRYBOUNDDO34I=1,L1DO35J=1,M1IF(I==1)THENDTX(I,J)=(T(I+1,J)-T(I,J))*L1/XLELSEIF(I==L1)THENDTX(I,J)=(T(I,J)-T(I-1,J))*L1/XLELSEIF(N0.LE.I.AND.I.LE.N00)THENDTX(I,J)=0.0ELSEDTX(I,J)=(T(I+1,j)-T(I-1,J))*L1/XL/2.0ENDIF35ENDDO34ENDDODO44J=1,M1DO45I=1,L1IF(J==1)THENDTY(I,J)=(T(I,J+1)-T(I,J))*M1/YLELSEIF(J==M1)THENDTY(I,J)=(T(I,J)-T(I,J-1))*M1/YLELSEIF(M0.LE.J.AND.J.LE.M00)THENDTY(I,J)=0.0ELSEDTY(I,J)=(T(I,j+1)-T(I,J-1))*M1/YL/2.0ENDIF45ENDDO44ENDDODO24I=1,L1DO25J=1,M1EX(I,J)=0.746*U(I,J)*CPCON*T(I,J)-COND*DTX(I,J)25ENDDO24ENDDODO26J=1,M1DO27I=1,L1EY(I,J)=0.746*V(I,J)*CPCON*T(I,J)-COND*DTY(I,J)27ENDDO26ENDDO4RETURNENTRYOUTPUTIF(ITER==0)THENPRINT401WRITE(8,401)401FORMAT('ITER',7X,'V(4,7)',19X,'T(4,7)')ELSEPRINT403,ITER,V(4,7),T(4,7)WRITE(8,403)ITER,V(4,7),T(4,7)403FORMAT(I6,1P4E15.3)ENDIFIF(ITER==LAST)THENCALLPRINTENDIFRETURNENTRYGAMSORDO500J=1,M1DO501I=1,L1GAM(I,J)=AMUIF(NF==4)THENGAM(I,J)=COND/CPCONIF(M0.LE.J.AND.J.LE