二维椭圆型流动传热通用程序.docx

1、二维椭圆型流动传热通用程序二维椭圆型流动传热通用程序变量表及算例说明（本材料仅供教学参考）西安交通大学CFD&NHT/EHT 研究中心陶文铨教授2002/10/15一、 FORTRAN 变量表 3二、 关于程序的主要说明 6(1)二维椭圆型流动和传热问题通用计算机算法方面的特点 6(2)各程序的主要功能 7(3)三种坐标系统 8(4)网格系统与节点命名方法 9十一个例题的已知条件与求解内容 12例题1直角坐标中二维稳态无内热源的导热 12例题2空心圆柱内的稳态热传导 12例题3正方形管道内的充分发展对流换热 13例题4内壁上有直肋的环形通道内充分发展对流换热 13例题5给定流场条件下温度场的计

2、算 14例题6二维突扩通道中的流动与换热 14例题7方形通道内的复杂充分流动 15例题8旋转圆盘上的冲击流动 15例题9轴对称燃烧内的瞬间燃烧过程 16例题10带中心射流的通道内的紊流换热 16例题11有质量源的流动问题 17FORTRAN 变量表List of FORTRAN VariablesACOFQuantity calculated by subroutine DIFLOW to give the combined convection and diffusion effect.AIM (I, J)The coefficient On .AIP (I, J)The coefficie

3、nt aE .AJM (I, J)The coefficient Os .AJP (I, J)The coefficient On .AP (I, J) tThe coefficient ap ; also Sp in GAMSOR.APTThe unsteady term / 匚 t .AREALocal variable, usually the area of a C.V . face.ARHOaLocal variable, (area)( ).ARX (J)The area of the main C.V . face normal to the x direction.ARXJ (

4、J)The part of ARX (J) that overlaps on the C.V . for V (I, J).ARXJP (J)The part of ARX (J) that overlaps on the C.V . for V (I, J+1).BLBLCCoefficients used in the block correction.BLMBLPCON (I, J)The constant term b in the discrimination equation; also stands for SC in GAMSOR.DENOMTemporary storage.

5、DIFFDiffusion conductance D.DTThe time step t .DU (I, J)de Influencing U (I, J).DV (I, J)d n Influencing V (I, J).F (I, J, NF)Various .FLTemporary storage leading to FLOW.FLMTemporary storage leading to FLOW.FLOWMass flow rate through a C.V . face.FLPTemporary storage leading to FLOW.FV (J)Interpola

6、tion factors which give the mass flowFVP (J)Vr At a main grid point, I, J as FV (J) V(I, J)+FVP (I, J)：? vr (i,j+1)FX (I)Interpolation factors, which give the interface.FXM (I)Density RHOM (at the location of U (I, J) as FX (I) - RHO (I, J)+FXM (I)RHO (I-1, J).FY (J)Interpolation factors, which give

7、 the interface.FYM (J)Density RHOM (at the location of V (I, J) as FY (J) RHO (I, J)+FYM (J)RHO (I, J-1).GAM (I, J)The diffusion coefficient -.IIndex denoting the position in x.IBEGTemporary values used in PRINT.IENDIFSTThe first value of I for which the print-out is arranged; used in PRINT.IITempor

8、ary index.IPREFThe value of I for the grid point, which is used as a reference for pressure.1STThe first internal-point value of I.ISTFIST-1; used in SOLVE.ITERA counter for iterations.IT1Temporary values used in SOLVE.IT2JIndex denoting the position in y.JFLTemporary index used in PRINT.JFSTSimilar

9、 to ILST.JJTemporary indexJLSTSimilar to ILSTJPREFSimilar to IPREF.JSTThe first internal-point value of J.JSTFJST-1; used in SOLVE.JT1Temporary values used in SOLVE.JT2LASTThe maximum number of iterations allowed by the user.LBLK (NF)When. TRUE. The block correction for F (I, J, NF) is used.LPRINT (

10、NF)When. TRUE. , F (I, J, NF) is printed.LSOLVE (NF)When. TRUE. , We solve for F (I, J, NF).LSTOPWhen. TRUE. Computation stops.L1The value of I for the last grid location in the x direction.L2(L1-1).L3(L1-2).MODEIndex for the coordinate system; =1 for Xy, =2 forX , =3 for二.M1The value of I for the l

11、ast grid location in the y direction.M2(M1-1).M3(M1-2).NTemporary storage for NF.NF6Index denoting a particular .NFMAXThe largest value of NF for which storage is assigned.NGAMNFMAX+3; GAM (I, J) can be considered as F (I, J, NGAM).NPNFMAX+1; P (I, J) can be considered as F (I, J, NP).NRHONFMAX+2; R

12、HO (I, J) can be considered as F (I, J, NRHO).NTIMES (NF)The number of repetitions of the sweeps in SOLVE for the variable F (I, J, NF).P (I, J)The pressure p.PC (I, J)1The pressure correction p .PREFThe pressure at the reference point.PT (I) or PT (J)QT (I) or QT (J)Transformed coefficients in the

13、TDMA.R (J)The radius r for a main grid point I, J.REL1.0-RELAX (NF).RELAX (NF)Relaxation factor for F (I, J, NF).RHO (I, J)pThe density .RHOCONThe value of for a constant-density problem.RMN (J)The value of radius r for the location to which V (I, J) refers.SMAXIThe largest absolute value of the “ m

14、ass source ” PusequaAiohne1SSUM1The algebraic sum of all the “ mass sources p iietfisition.SX (J)Scale factor for the x direction at the main grid locations Y (J).SXMN (J)Scale factor for the x direction at interface locations YV(J).TEMPTemporary storage.TIMETime t for unsteady problems.TITLE (NF)Al

15、phabetic title for F (I, J, NF).U (I, J)The x-direction velocity u.V (I, J)The y-direction velocity v.VOLVolume of the C.V.X (I)The values of x at grid points.XCV (I)The x-direction widths of main C.V .XCVI (I)The part of XCV (I) that overlaps on the C.V . for U (I, J).XCVIP (I)The part of XCV (I) t

16、hat overlaps on the C.V . for U (I+1,J)XCVS (I)The x-direction width of the staggered C.V . for U (I, J).XDIF (I)The difference X (I)-X (I-1).XLThe x-direction length of the calculation domain.XU (I)The locations of the C.V. faces; i.e. the location of U (I, J).丫 (J)The values of y at grid points.YC

17、V (J)The y-direction widths of main C.V .YCVR (J)The area r - y for a main C.V .YCVRS (J)The area r : y for the C.V. for v (I, J).YCVS (J)The y-direction width of the staggered C.V . for V (I, J).YDIF (J)The difference Y (J)-Y (J-1).YLThe y-direction length of the calculation domain.YV (J)The locati

18、ons of the C.V. faces; i.e. the location of V (I, J).二关于程序的主要说明(1) 二维椭圆型流动和传热问题通用计算机程序算法方面的特点1、 采用原始变量法，即以速度 U、V及压力P作为直接求解的变量2、 守恒型的差分格式，离散方程系对守恒型的控制方程通过对控制容积作积分而得出的，无论网格疏密程度如何，均满足在计算区域内守恒的条件；3、 采用区域离散化方法 B，即先定控制体界面、再定节点位置4、 采用交叉网格，速度 U、V与其他变量分别存储于三套网格系统中；5、 不同的项在空间离散化过程中去不同的型线假设， 源项采用局部线性化方法；扩散一对流项

19、采用乘方格式(但很容易转化为中心差分、迎风差分或混合格式) ；街面上的扩散系数采用调和平均法，而密度与流速则用线性插值；6、 不稳态问题采用全隐格式，以保证在任何时间步长下均可获得具有物理意义的解；7、 边界条件采用附加源项法处理；8、 耦合的流速与压力采用 SIMPLE算法来求解；9、 迭代式的求解方法，对非线性问题，整个求解过程具有迭代性质；对于代数方程也采用迭代法求解；10、 采用交替方向先迭代法求解代数方程并补以块修正技术以促进收敛。（ 2） 各方程块的主要功能1、 主程序 MAIN : 规定整个计算过程的流程，决定是否停止计算。2、 子程序 SETUP1: 设置与网格系统有关的，在计

20、算过程中保持不变的几何参数及设置U、V、P P、p、CON（Sc）、AP（Sp）的初值。3、 子程序SETUP2 : 1.建立每一轮迭代中变量 U、V及P和其他变量离散方程的系数 as、 aw、a”、as、ap及b项。所有变量公用一套数组逐个求解，逐个确定各自目标的系数，以节省内存。在建立系数过程中调用 DIFLOW子程序；2调用SOLVE子程序求解代数方程； 3 迭代次数指标及时间步长增值；4、 子程序 DIFLOW ：计算乘方格式中的系数 A（|p|）;5、 子程序 SOLVE: 采用交替方向线迭代法并辅以块修正技术求解代数方程；6、 子程序 GRID ：设置为建立网格系统所必须的量，包括

21、（1） 设置 x,y 方向求解区域的宽度 XL 及 YL（2） 设置 x,y 方向上各自的节点数 L1 及 M1（3） 设置 x,y 方向上控制体积的界面位置， XU（I）,I=2,L1;YV（J）,J=2,M1 。如为均分网格可确定子程序 UGRID, 如为 非均分网格需由用户在 GRID 中把 XU（I）,YV（J） 一一设置好。（4） 规定坐标系， MODE=1 ， 2， 3 分别相应于直角、圆柱轴对称 及极坐标（5） 当 MODE 不等于 1 时，规定径向起始点半径 R（I）7、 子程序 START: 设置初值，包括（1 ）、对不稳态问题给出初始条件（2） 、对稳态问题给出迭代求解的假

22、定值（3） 、已知的边界值也可在此块中，在设置初值或假定值时一并送入，但对边界条件 随时间而异的不稳态不能在此块中赋值，因在整个计算过程中此会仅执行一次。8、 子程序 DENSE:规定密度场，对于密度为常数的问题，此块中可不设任何语句，但必 须保留 ENTRY 及 RETUIRN 语句9、 子程序 BOUND: 规定边界条件，包括（1） 设置各变量的边界条件（2） 对各特殊问题所需的量，如 Nu, fRe 等等可在此块中计算，但也可在 OUTPUT 中设置10、 子程序 OUTPUT ：打印输出，包括（1） 每做一轮迭代，输出一行信息，以观察收敛情形；（2） 调用 PRINT 子程序，实现二维

23、物理量场的输出；（3） 用户所需输出的其他特殊变量11 、 子程序 GAMSOR ：设置扩散系数及源项，包括（1） 设置控制方程中规定的源项 Sc,Sp（2） 对第二，三类边界条件设置与边界相邻的控制体中的附加源项， Sc,ad及Sp,ad（3） 内节点及边界节点规定扩散系数，扩散采用附加源项法时令边界扩散系数为零。12、 子程序 UGRID: 为均分网格设置界面位置13、 子程序 PRINT: 输出二维物理场(3)二种坐标系统 直角坐标系，Z方向厚度为1;M0DE=1。x,l轴对称圆柱坐标，计算对0 =1弧度 的区域进行；M0DE=2 ; Y(J)可以 从任何起点算起，但 R(J)须从对称

24、轴开始；R(1)维计算区的最小半径。两个区域XL,YL相等，但 R不 同极坐标系，Z方向厚度为1; M0DE=3;计算区域的0角须小于2n ; Y(J)可从任何点起算，但 R(J)比从中心点开始；R(1)为计算区域边界的最小半径。MODEX(I)Y(J)R(J)SX(J)1xy1.01.02xyr1.030yrr网格系统与节点命名方法4V(2,M1)U(L1,M2)：U(2,M2)Fi * * i T 1 I 丨 p(L1i,M1) * “F*jL ! i I L !XCV(I)(c)x方向的主要几何参数ARXJP(J)M1 y(M1) yv(M1)XCVS(3)M2U(2,J)XDIF(3)

25、(d)速度 U(3,J)的控制容ARX(J)ARXJ(J)YDIF(J) ycvs(j)ycv(J-1)yv(3)2 y(2)1 y(i) yv(2)V(I,J) | I1 1 I , J 3 y(3)(e)y方向的主要几何参边界压力P(1,2)、P(2,1)由内点值境外推而得；参考压力点p(1,1)按下式计算：P(1,1)=P(1,2)+P(2,1)-P ( 2,2)变量内点下标起始值变量1STJST,P,P22U32V23三、一个例题的已知条件与求解内容例题1直角坐标中二维稳态无内热源的导热已知：求解区域如图 1所示，四个边界上的温度由下式决定； T=x+y+xy求：该区域中的温度分布。图

26、1例题2空心圆柱内的稳态热传导已知：求解区域如图 2所示。左边界为给定温度，右边界为对流，上边界为绝热，下边界为 给定热流。整个计算区域内都有源项。求：该区域中的温度分布。计算方法：用附加源项法处理第 2, 3类边界条件例题3正方形管道内的充分发展对流换热已知：常物性流体在均匀壁温的正方形截面管道内作充分发展的对流换热（如图 3）。求：截面上的速度分布、温度分布、 fRe及Nu。计算方法：充分发展对流换热问题的处理特点例题4内壁上有直肋的环形通道内充分发展对流换热已知：常物性流体在带直肋的环形通道内作充分发展层流对流换热，内表面为均匀壁温（周 向）,轴向则呈线性变化，外表面绝热。 Ri=1,R

27、2=2, a =15o。求：截面上速度分布、温度分布及 fRe、Nu。计算方法：耦合问题的一种处理方法。例题5给定流场条件下温度场的计算已知：流体流过一个直角的两表面，速度场为 u=Ax,V=-AY,a=10,T in=500,Tw=100求计算区域中温度的分布。TinT out绝热例题6二维突扩通道中的流动与换热已知：由两平行平板组成的突扩通道，尺寸如图所示。流动为层流。入口流速均匀， Vin=100,入口温度均匀，Tin=100,通道壁温均匀，Tw=300。流体Pr=0.7,分子粘性 尸1，密度按下式变化，密度分布及压力场。计算方法：出口边界条件的处理图6例题7方形通道内的复杂充分流动已知

28、：常物性流体在如下图示方形管道内作充分发展的流动与换热 （重力项中的密度采用Boussines 假设）.p3=10：T仁0,T2=1,PR=0.7,尸1.0 dp/dz=const 计算中取-3000,进入充分发展时 右边进入的热量等于左边导出的热量 求截面上速度 U、V的分布，Z方向的分量 W的分布，温度分布及压力场。计算方法：如何用一个二维程序来计算三个速度分量。例题9轴对称燃烧内的瞬间燃烧过程已知：燃料与空气如图 9所示进入轴对称燃烧室，设：1、 仅有三种组分：燃料，氧气及燃烧产物2、 氧气与燃料的扩散系数相等3、 反应在瞬间内完成4、 对各个组分Cp均相同U1=2O,U2=5O,尸 1,Pr=0.7 燃料焓 Hfu=3X 1O【Vo=O.2U2 求：燃料室中的速度场、温度场与浓度场。例题10带中心射流的通道内的紊流换热已知：如图10，带一般中心射流的流体进入一平行板通道。流动为紊流， 卩=106T=100V=10例题11已知：一深为0.2米的湖泊尺寸如图11所示。湖泊中有一源与汇。湖面上与空气发生换热。Tair=100, h=2。Uin=100,Tin=500,湖的四个边均视为绝热。源处的入口强度为 4。湖深方向当作均匀处理（即按 X Y坐标中的二维问题处理）。求：湖泊中的速度、温度分布计算方法：质量源项方法的应用。图11