模拟辐射传热的离散坐标法的改进

资源类型:pdf 资源大小:845.00KB 文档分类:工业技术 上传者:周韦臣

相关文档

批量下载下列文档

文档信息

【作者】 杨占春  武文斐  李义科 

【关键词】辐射传热 离散坐标法 假散射 分段积分 

【出版日期】2005-05-25

【摘要】 针对离散坐标法模拟求解炉内辐射换热的问题,在对辐射传递方程采用控制容积法离散时,按辐射方向上各界面之间的投影关系划分区域,对辐射传递方程在各子控制体上沿辐射传输方向进行积分,求解各子控制体的辐射强度,再取容积平均值作为最终节点的辐射强度的新离散坐标方法,定义为分段积分离散坐标法,并对三维矩型燃气炉内辐射换热进行了数值模拟,与传统离散坐标方法及区域法进行了比较,比较结果表明新离散坐标法较好的解决了假散射问题,提高了离散坐标法的精度.

【刊名】计算物理

全文阅读

0 引言随着计算机技术和计算数学的迅速发展,辐射换热数值计算作为传热学的重要组成部分,近20多年得到了空前发展.离散坐标法的主要优点是将辐射传递方程转化成微分表达式,便于同一般输运方程耦合求解,计算工作量小,而且近几年出现了几种新的离散坐标取法,如SRAPN[1]法和T[2]N法,使得离散坐标法的精度大幅提高,可以和区域法媲美,因此离散坐标法将是今后一段时期内发展和应用较多的一种辐射换热模型.然而,其本身仍存在难以克服的缺点,如射线影响、假散射、对差分格式敏感等,笔者认为假散射的根源在于对离散方程进行积分时,是按坐标方向进行的,没有按传播方向进行,从表1中可以看出,在很多情况下,W面根本影响不到E面,而是对T面和N面产生影响,传统做法忽略了这一点,W面根本和E面没有关系,却用W面来推算E面,相反,S面和B面必然对E面产生影响,却没有考虑进来,为此,本文提出一种按各面之间在传播方向上的投影关系,将微元控制体划分区域,分段对辐射传递方程进行积分的方法,定义为分段积分离散坐标法.1 辐射传递方程的坐标离散离散坐标法是基于对辐射强度的方向变化进行离散,通过求解覆盖整个4π立体角空间的一套离散方向上的辐射传递方程而得到问题的解,在三维直角坐标系中,对于吸收、辐射、各向同性散射介质的辐射传递方程可写为μIx+ξIy+ηIz=-(kα+ks)I+kαIb+ks4π∫4πIΦ(Ω,Ω′)dΩ′,(1)漫灰表面的边界条件为I=εIb,w+ρwπ∫n·s′≤0I|n·s|dΩ′.(2)在离散坐标法中,积分项用数值积分代替,写为μmImx+ξmImy+ηmImz=-(kα+ks)Im+kαIb+ks4π∑m′ωm′m,m′Im′,(3)式中的方向余弦μm,ξm,ηm和权值ωm根据离散坐标取法确定,这里不做阐述,参见相关文献.对于漫-灰体边界壁面,相应的边界条件为May,2005 Im=εwIb,w+ρwπ∑m′,n·s≤0ωm′|n·s|dΩ′.(4)  先按传播方向上各界面之间的投影关系划分区域,以方向余弦μm,ξm,ηm>0为例,参见图1(a),(b),(c).其中,ABCD面称为控制体的B界面、EFGH为T界面、ACEG为W界面、BDFH为E界面、ABEF为S界面、CDGH为N界面.分3种情况:(a)从A点引射线交于T面P1点;(b)从A点引射线交于N面P2点;(c)从A点引射线交于E面P3点.然后在各子控制体上分别对(3)式积分,如第1种情况下的VWT,m子控制体W T,μm(-AWT,mIWT,m)+ηm(ATW,mITW,m)=-βm,pVWTIWTP,m+Sm,pVWT,(5)其中 βm,p=kα+ks-ks4πωmm,m, Sm,p=kαIb,p+ks4π∑m′≠mωm′m,m′, AW=ΔxΔy, AWT=AW-12ξmηmΔz2,  ATW=μmηmAW-12ξmηmΔz2,  VWT=12μmηmΔyΔz2-13μmξmη2mΔz3.令IWTP,m=fITW,m+(1-f)IW,m,即ITW,m=[IWTP,m-(1-f)IW,m] f,(6)将(6)式代入(5)式有IWTP,m=(μmAWT,mfIW,m+ηm(1-f)ATW,mIW,m+fVWT,mSm,p) (ηmATW,m+fβm,pVWT,m).(7)将(7)式代入(6)式可推得ITW,m.对其他各子控制体以同样方式进行积分,下面仅给出相应子节点的计算式,和要用到的相应子控制体的界面面积和体积,过程从略.  W-N:IWNP,m=(μmAWN,mfIW,m+ξm(1-f)ANW,mIW,m+fVWN,mSm,p) (ξmANW,m+fβm,pVWN,m),(8)AWN=12ξmηmΔz2,  ANW=12μmηmΔz2,  VWN=16μmξmη2mΔz3;  W-E: 无;  S-T:ISTP,m=(ξmAST,mfIS,m+ηm(1-f)ATS,mIS,m+fVST,mSm,p) (ηmATS,m+fβm,pVST,m),(9)AS=ΔxΔz, AST=AS-12μmηmΔz2, ATS=ξmηmAS-12μmηmΔz2, VST=12μmηmΔxΔz2-13μmξmη2mΔz3;  S-N: 无;  S-E:ISEP,m=(ξmASE,mfIS,m+μm(1-f)AES,mIS,m+fVSE,mSm,p) (μmAES,m+fβm,pVSE,m),(10)ASE=12μmηmΔz2, AES=12ξmηmΔz2, VSE=16μmξmη2mΔz3;  B-T:IBTP,m=(ηmABT,mfIB,m+ηm(1-f)ATB,mIB,m+fVBT,mSm,p) (ηmATB,m+fβm,pVBT,m),(11)AB=ΔxΔy, ABT=ATB=AB-ξmηmAS-μmηmAW+μmξmηmΔz2, VBT=ABTΔz;  B-N:IBNP,m=(ηmABN,mfIB,m+ξm(1-f)ANB,mIB,m+fVBN,mSm,p) (ξmANB,m+fβm,pVBN,m),(12)ABN=ATS, ANB=AS-ANW, VBN=AWNΔx-2VWN;  B-E:IBEP,m=(ηmABE,mfIB,m+μm(1-f)AEB,mIB,m+fVBE,mSm,p) (μmAEB,m+fβm,pVBE,m),(13)ABE=AB-ABN-ABT, AEB=AW-AES, VBE=ASEΔy-2VSE.完成各子控制体的求解后,我们对整个控制体ABCDEFGH取容积平均值得IP,m=(IWTP,mVWT,m+IWNP,mVWN,m+ISTP,mVST,m+ISEP,mVSE,m+IBTP,mVBT,m+IBNP,mVBN,m+IBEP,mVBE,m) V总.(14)对各下游界面取面积平均值有(a)从A点引射线交于T面P1点(a)TheradialfromAtoP1onsurfaceT(b)从A点引射线交于N面P2点(b)TheradialfromAtoP2onsurfaceN(c)从A点引射线交于E面P3点(c)TheradialfromAtoP3onsurfaceE图1 按传播方向上各界面之间的投影关系划分区域(μm,ξm,ηm>0)Fig.1 Thetinydominatedbodyisplottedintodistrictsaccordingtotheprojectionrelationbetweentheinterfacesintheradiativeheat transferingdirection(μm,ξm,ηm>0)IT,m=(ITW,mATW+ITS,mATS+ITB,mATB) AT,(15)IN,m=(INW,mANW+INB,mANB) AN,(16)IE,m=(IES,mAES+IEB,mAEB) AE.(17)表1中另外两种情况的推导与此类似,当μm,ξm,ηm取其他值时也完全可以仿照上面的过程推导,这里不再赘述.上面的工作看起来很烦琐,但一旦完成,辐射传递方程的求解就转化为对式(6)~(17)的求解,可以借助计算机程序来实现了.首先假设温度场,对离散方程(6)~(17)采用逐点推进的方法求解,对每一个离散方向,都从其上游边界开始计算,利用式(7)~(17)计算节点处的辐射强度,再由一系列与式(6)形式一样的式子计算下游界面的辐射强度,并作为下一节点的上游界面辐射强度,按此从一个边界逐点推进到另一个边界.每一个离散方向都求解完后,将辐射强度值代入能量平衡方程,可求出温度分布,代替上一次假设的温度场,再求解辐射强度,如此反复求解,直到满足迭代精度,求解完毕.2 三维矩形炉膛内辐射传递过程的数值求解为了验证此方法正确性,本文采用诸多文献都使用的算例:求解一三维矩形模型炉膛.计算区域划分成7×7×13的网格,矩形炉膛的有关参数如下:尺寸及物性参数:Lx=Ly=2m, Lz=4m, kα=0 5m-1, ks=0, q=5 0kW·m-3;边界条件:z=0:T=1200K,ε=0 85;  z=Z:T=400K,ε=0 7;  其他边界:T=900K,ε=0 7.新方法采用阶梯格式(即f=1 0)计算,传统方法分别采用阶梯格式和指数格式进行计算,坐标的离散均采用TN(N=6)法,并与区域法的计算结果进行了比较.图2(a)给出了在y=1 0,z等于0 4,2 0,3 6处沿x方向介质温度的分布规律,图2(b)给出了在y=1 0,z等于0 0,4 0处沿x方向壁面热流密度的分布规律,从图2中可以看出,采用同样的阶梯差分格式,新方法比传统方法的结果更接近区域法结果,但比传统方法采用指数格式计算的结果稍差,除了上部区域误差较大,其他部位与区域法结果吻合的较好.表1给出了区域法计算的各点温度值及在阶梯差分格式下新旧方法的结果,并计算了误差,在z=0 4处,平均误差减少了0 19%;在z=2 0处,传统方法计算结果偏高.平均误差为0 16%,新方法计算结果偏低,平均误差-0 25%,平均误差增加了0 09%;在z=3 6处,平均误差减少了1 05%.从总体来看,在相同差分格式下,计算精度有所提高.图2 新离散坐标法与传统离散坐标法及区域法结果比较Fig 2 AcomparisonofthenewDOM,thetraditionalDOM,andthezonemethod —区域法[3~5],—新离散坐标法,传统离散坐标法指数格式,·传统离散坐标法阶梯格式( —zonemethod[3~5],—newDOMinladderformat,traditionalDOMinexponentformat,·traditionalDOMinladderformat)表1 各点温度及误差表Table1 Temperaturesanderrors坐标值区域法结果 K传统离散坐标法结果 K新方法结果 K传统方法误差 %新方法误差 %(0 2,1 0,0 4)1039.91028.81029.4-1.07-1.00(0 6,1 0,0 4)1056.41041.81044.2-1.38-1.15(1 0,1 0,0 4)1061.91045.01048.2-1.59-1.29(1 4,1 0,0 4)1056.41041.81045.1-1.38-1.04(1 8,1 0,0 4)1039.91028.81028.8-1.07-1.07(0 2,1 0,2 0)946.5948.2945.30.18-0.13(0 6,1 0,2 0)950.8952.6947.90.19-0.31(1 0,1 0,2 0)953.1953.9949.00.08-0.43(1 4,1 0,2 0)950.8952.6948.30.19-0.26(1 8,1 0,2 0)946.5948.2945.40.18-0.12(0 2,1 0,3 6)880.6888.2882.10.860.17(0 6,1 0,3 6)872.4884.5873.71.390.15(1 0,1 0,3 6)870.8883.6872.11.470.15(1 4,1 0,3 6)872.4884.5873.71.390.15(1 8,1 0,3 6)880.6888.2881.10.860.06传统方法平均误差为-1 30%,新方法平均误差为-1 11%,误差减少值为0 19%传统方法平均误差为0 16%,新方法平均误差为-0 25%,误差减少值为0 09%传统方法平均误差为1 19%,新方法平均误差为0 14%,误差减少值为1 05%3 结论与分析①该方法的应用,消除了假散射问题,提高了计算精度.②不用细化网格,节约机时,运行效率提高.③离散坐标法常用的3种差分格式中(阶梯格式、菱形格式、指数格式),阶梯格式精度最差,指数格式精度最高,而本文采用分段积分法计算后,阶梯格式计算误差有所降低,在高温区,与区域法吻合很好,基本上达到了指数格式计算的精度,说明新方法在消除假散射问题上的有效性

1 2

问答

我要提问