《计算传热学程序设计.doc》由会员分享，可在线阅读，更多相关《计算传热学程序设计.doc（29页珍藏版）》请在三一文库上搜索。
1、中国石油大学（华东）储运与建筑工程学院热能与动力工程系计算传热学程序设计设计报告学生姓名： 学 号: 专业班级： 指导教师2012年7月7日1、设计题目有一房屋的砖墙厚 8=0.3 m , 2=0.85 W/(m ), p=1.05 X 06 J/( m3 K)，室内温度 Tfi保持20r不变，表面传热系数hi=6W/(m2C)。开始时墙的温度处于稳定状态，内墙 表面温度Twi为15C寒潮入侵后，室外温度Tf2下降为-10C，外墙的表面传热系数为 35W/(m2 C)。试分析寒潮入侵后多少时间内墙壁面方可感受到外界气温的变化。1.1已知参数壁厚，墙壁导热系数，密度与比热容的乘积，室内和寒潮入侵
2、后室外空气温度，室 内空气和外墙的表面传热系数，开始时稳定状态下的内墙表面温度。1.2求解寒潮入侵多少时间后内墙壁面可感受到外界气温的变化？2物理与数学模型2.1物理模型该墙面为常物性，可以假设：(1)其为无限大平面，(2)只有在厚度方向传热, 没有纵向传热，则该问题转化为一维常物性无限大平面非稳态导热问题。2.2数学模型以墙外表面为坐标原点，沿厚度方向为坐标正方向，建立坐标系。基于上述模型， 取其在x方向上的微元作为研究对象，则该问题的数学模型可描述如下:T)(1a)初始条件:(1b)在两侧相应的边界条件是第三类边界条件，分别由傅立叶定律可描述如下:左边界:|x 0h2 (Tx 0 Tf2)
3、( 1c)X右边界：xhi(TxTfi)(1d)X3数值处理与程序设计3.1数值处理采用外点法用均匀网格对求解区域进行离散化，得到的网格系统如图2所示。一共使用了 0N-1共N个节点。dl图2墙壁内的网格划分节点间距s为:(2a)此例中墙壁导热系数为常值，无源项。则可采用有限体积法对控制方程离散化，得 到离散方程为：apTpaETEawTwb式中:0apaEawap(2b)0c xaE，aw， ap(2c)xx. %-0bapTp(2d)其中的上标“0表示此为上一时刻的值，分别为节点所在控制容积左右边界上的导热系数，由于墙壁导热系数不变，故都等于人At为时间步长。由元体能量平衡法可以得知左右边
4、界节点的离散方程分别为：左边界节点：一: ' " - - ' _ '二 一'（3）右边界节点：r 二1. -亍八（4）离散方程的详细推导过程见附录。3.2程序设计由物理模型可以知道本问题为一维导热问题，一维导热问题的离散方程在取遍所有节点后形成的是三对角的代数方程组，采用追赶法进行求解。程序构成和方法：程序由主程序和一个子程序构成。主程序进行变量定义和各已知参 数的输入，以及左右边界节点和内部节点控制方程的输入；子程序tdma实现追赶法用来计算每个节点新的温度。Thomas算法求解过程分为两步：消元和回代。消元是从系数 矩阵的第二行起，逐一将每一行的非
5、零元素消去一个，使原来的三元方程化为二元方程。 消元进行到最后一行时，二元方程就化为一元方程，直接得到最后一个未知数的值。然 后逐一往前回代，由各二元方程求出其它未知解。程序特点：该程序有很强的适应性，一维常物性非稳态平壁导热问题都可以使用此 程序，只要适当更改边值条件即可。还可以进行修改解决非常物性问题。程序中对输出节点，最大输出量都进行了控制，对计算结果的分析有很大帮助。而且 Thoms算法的优点需要内存小，工作量小，程序设计简单。程序流程图：首先对变量赋值，然后由初始条件建立初始温度场，接着从左边界， 内部节点，到右边界进行迭代，直到满足精度要求为止，最后输出结果，程序结束。程 序流程如
6、下图3。4、模型与程序验证4.1模型本题简化为厚度为2 =0.3m的一维非稳态模型如图4所示，初始温度为15C，在 其中间建立坐标系，左两边为对流换热，且换热系数相同都为 h=25 W/（m2 -C），且流体 温度Tf=-1OC对于x 0,列出其导热微分方程式及定解条件：Q(0x0)(5)图3程序流程图aT(x,O) To(0 x )cT(x,)x 0x引入过余温度:h T( , ) TT(x,)(6)(7)(8)T(x, ) T(9)图4 一维导热简化模型2 -C)直接根据公式得到解析解如下:(,)0.Cnexp(n 12nFO)COS( n )(10)式中，Fo ay ,，系数Cn应该使上
7、述无穷级数在0是满足初始条件，由傅里叶级数理论可得：(11)(12)2si n cos n sinn是超越方程的根，称为特征根。Bitan n ,n 1,2,n其中Bi4.2程序验证(1) 由模型可以得到相关信息然后进行编程，同等时间下计算出中心处温度的解 析解和数值解进行比较，数据记录在表 1o然后计算出相对误差，作图5,观察数值解 与分析解的比较曲线。由图表中可以发现，平壁中心不同时刻温度值的分析解和数值解相差不是很大，二 者吻合的比较好，可以说明所编制的数值解法的程序是正确的。 相对误差先增大后减小, 增大的原因是此时温度接近零度，相对误差的基数比较小，所以造成相对误差较大，但 是此时的
8、绝对误差并不大，在合理范围内，所以除去个别点外，都满足误差小于百分之1。可以验证所编数值解法的程序是正确的（2）空间步长对墙内壁的温度影响如图 6及表2。在程序编写过程中用网格节点 数对空间步长进行控制，为了观察空间步长对墙内壁温度的影响，表中选择了三个不同 的空间步长，分别为选取51,101,201个网格节点，则相应的空间步长为0.006,0.003, 0.0015。根据不同步长时温度的变化曲线可以看出，空间步长对内墙壁的影响不大，当空间步长控制在合理范围时可以忽略空间步长的影响。表1分析解与数值解比较时间（h）分析解（C）数值解（C）相对误差（%）0151500.55614.85414.8
9、11-0.289481.11113.46313.425-0.282261.66711.30611.3-0.053072.2229.0669.0810.1654532.7786.9746.9990.3584743.3335.0825.1120.5903193.8893.3933.4250.9431184.4441.891.9231.7460325. 000.5540.5886.1371845.556-0.631-0.598-5.229796.111-1.684-1.651-1.959626.667-2.618-2.586-1.222317.222-3.447-3.417-0.870327.778
10、-4.184-4.154-0.717028.333-4.837-4.809-0.578878.889-5.417-5.391-0.479979.444-5.932-5.907-0.42144图5平壁中心不同时刻的数值解和分析解表2空间步长对温度影响数据时间(h)N=51N=101N=2010.00015.00015.00015.0000.50015.00015.00015.0001.00014.99814.99814.9981.50014.98014.98114.9812.00014.92414.92414.9242.50014.82014.82014.8203.00014.67514.675
11、14.6763.50014.50114.50214.5024.00014.31014.31114.3114.50014.11114.11214.1125.00013.91113.91113.9125.50013.71513.71513.7156.00013.52513.52513.5256.50013.34313.34413.3447.00013.17113.17213.1727.50013.00913.01013.0108.00012.85812.85812.8588.50012.71612.71612.7169.00012.58312.58312.5839.50012.46012.4601
12、2.4609.91712.36312.36312.36414.Q13.512+512,01 JIIllJIJJI2 3456789IO11时间<h>图6空间步长对墙内壁温度影响(3) 时间步长对温度的影响如图7和表3,根据图中曲线可以看出时间步长选择 50s, 100s, 200s时基本重合，对墙内壁温度影响不大15.5-«-T=100s-*-T=200s24G时间(h)图7时间步长对温度的影响S10015015.0000150.694150.69415.0000.667151.38914.9881.38914.9871.33314.9872.08314.9122.083
13、14.910214.9182.77814.7462.77814.7442.66714.773.47214.5133.47214.5113.33314.5584.16714.2454.16714.244414.3084.86113.9664.86113.9664.66714.0445.55613.6925.55613.6935.33313.7816.2513.4316.2513.433613.5276.94413.1886.94413.1906.66713.2887.63912.9647.63912.9667.33313.0668.33312.768.33312.762812.8629.0281
14、2.5749.02812.5768.66712.675表3时间步长对温度的影响时间(h)T=50 (s)时间(h)T=100 (s)时间(h)T=200 (s)b图8墙内壁温度随时间的变化曲线5计算结果与分析5.1墙内壁温度分析根据题目中要求，计算寒潮入侵多长时间后内墙壁可以感受到外界气温的变化，通过建模，方程离散化，最终通过程序求解方程，得到图8和表4。由图可以看出，开始阶段，内墙壁温不变，随着时间的进一步深入，内壁温度开始降低，当很长时间后，温 度变化基本趋于平缓，直到再次平衡。根据图 8就可以得到墙内壁温度开始发生变化的 时间。表4墙内壁温度随时间变化数据表时间（h）温度C）时间（h）温
15、度（C）时间（h）温度（C）0.000156.66713.28513.33311.7630.556157.22213.09813.88911.6911.11114.9977.77812.92414.44411.6261.66714.9678.33312.76215.00011.5652.22214.8838.88912.61215.55611.512.77814.7449.44412.47316.11111.4583.33314.56210.00012.34516.66711.4113.88914.35410.55612.22717.22211.3684.44414.13311.11112.1
16、1817.77811.3295.00013.91111.66712.01818.33311.2925.55613.69312.22211.92618.88911.2596.11113.48412.77811.84119.44411.2285.2导热系数入对墙内壁温度的影响墙的导热系数对内表面的影响，在图 9和表5中发现，导热系数对内壁温度影响比 较大，2=1.2时，温度下降的趋势会更快，要比 2=0.85时下降快的多，下降速度更快， 更短时间内达到稳态，因此导热系数越大温度扩散越快，导热系数越小则温度变化越慢, 需要更长时间到达稳态，但是这时对于要求恒温的空间有好处，受波动影响更小。表5导热系
17、数对墙内壁温度的影响时刻（h）入=0.85入=1.0入=1.201515150.5561515151.11114.99714.9914.9711.66714.96714.92414.8332.22214.88314.7714.5642.77814.74414.54214.2073.33314.56214.26613.8083.88914.35413.96813.3984.44414.13313.66612.999513.91113.37112.625.55613.69313.08812.2676.11113.48412.82311.9426.66713.28512.57511.6447.222
18、13.09812.34611.3727.77812.92412.13611.1258.33312.76211.94210.9018.88912.61211.76410.6989.44412.47311.60210.5140246时刻h)8 103 21 111图9导热系数对墙内壁温度的影响5.3墙外换热系数h的影响墙外表面传热系数对温度分布的影响，如图10和表6影响不大。151412 101 1 1on4211816图10墙外换热系数对温度影响表6墙外换热系数对温度影响时刻h=20h=35h=50(h)W/(m2 C)W/(m2 C)W/(m2 C)01515151.11114.99814.9
19、9714.9962.22214.91414.88314.8653.33314.66214.56214.5074.44414.31114.13314.0435.55613.93713.69313.5736.66713.5813.28513.1447.77813.25712.92412.7688.88912.97212.61212.4461012.72412.34512.17311.11112.5112.11811.94212.22212.32511.92611.74813.33312.16711.76311.58514.44412.03111.62611.44815.55611.91511.51
20、11.33416.66711.81511.41111.23717.77811.7311.32911.15618.88911.65711.25911.0885.4墙内壁传热热系数的影响由图11可以看出，墙内表面传热系数对内表面温度影响较大，当传热系数比较小时受墙外流体温度影响明显，传热系数越大受墙外流体温度影响越小，当到了极限大的 情况下，内表面温度则等于墙内流体温度，不再受墙外流体温度影响彳iII1ii*II024 S 8101214161820时刻(h)0注窓图11墙内表面传热系数对温度影响表7墙内表面传热系数对温度影响时刻h=2h=6h=10(h)W/(m2 :C ) W/(m 2 :C
21、) W/(m 2 :C )0.00015.00015.00015.0001.11114.99414.99714.9992.22214.78914.88314.9553.33314.16914.56214.8374.44413.29214.13314.6885.55612.33213.69314.5406.66711.39013.28514.4087.77810.50712.92414.2958.8899.70012.61214.20110.0008.97212.34514.12311.1118.31812.11814.05912.2227.73511.92614.00713.3337.2141
22、1.76313.96414.4446.74911.62613.92915.5566.33611.51013.90016.6675.96811.41113.87617.7785.64011.32913.85718.8895.34811.25913.8415.5墙壁厚度S对内壁温度的影响改变墙壁的厚度3,内壁温度将发生变化。在表8和图12中可以发现，不同厚度的 墙壁对外界温度的感应快慢是不一样的，随着墙壁厚度的增加，感应到外界温度变化的 时间越来越长，且温度变化越来越慢，墙壁越厚要越长的时间才能达到新的平衡16024& K 101214 lb 1«20时何（町图12 墙壁厚度对内
23、壁温度的影响表8墙壁厚度对内壁温度的影响时间（h）0.15m0.3m0.45m01515151.11113.86414.997152.22211.35414.88314.9993.3339.56814.56214.9894.4448.43614.13314.9555.5567.72713.69314.8946.6677.28513.28514.817.7787.00912.92414.7138.8896.83612.61214.608106.72812.34514.50111.1116.66112.11814.39712.2226.61911.92614.29613.3336.59211.76
24、314.214.4446.57611.62614.1115.5566.56611.5114.02616.6676.55911.41113.94917.7786.55511.32913.87718.8896.55311.25913.8116结论(1) 整个墙壁经历了以下变化过程：首先外壁直接与室外冷空气接触，温度变化很快,随着时间的推移，墙内各点的温度也开始变化，并影响到右边内墙壁的温度，慢慢降低(2) 对于墙内壁的温度随时间的变化，变化趋势总是由快到慢，最后重新达到稳态。 当改变墙壁厚度的时候，随着墙壁厚度的增加，对于外界温度的感应是越来越慢。这对 于一些对温度有要求的地反很重要，根据温度的变
25、化作出适当的调整措施。参考文献1 黄善波 刘中良编著 计算传热学基础2 杨世铭陶文铨编著传热学(第四版)北京高等教育出版社3 王元明编数学物理方程与特殊函数(第三版)高等教育出版社附录1数学模型的离散过程推导2 P 彳PP 彳p11P1 P WPXXX_ 丄 _ 丄e2X P X X ex wX X按照Taylor级数展开法，温度对时间的偏导有向前差分格式，中心差分格式和向后差分格式，使用向后差分格式向后差分格式：pp Pp 1Pp Pp 1PaPp E11pp WxxexexwPxwxpp 1p E11pp W3)得XX由axP2Pxexecxxw0PxPxExWPaIpPaEEa3WWba
26、p aEawap3e0，aPx空 b aPTp05)(7)边界条件处理：右边界根据元体能量平衡法处理:如图1:其中(8)Bhr (TfrTn 1)， n 2(9)hr (TfrTN 1)CX N 1为常数a c ,xx2（11）将式（12-b）（12-c）入式（12-a）以得到：c xhrN 1N 2hTc xT0hr 1 frTN 1（2xx2同理可以得到左边界点离散方程:(hi:x)T。hTflT-c xT0（13）x2x2因此，离散方程为：apTpaETEawTwb（14）式中，apaEaW0aP（15）aE，aw，a0c xp（16）xxb0 0apTp（17）式中，上标“ 0表示上一
27、时刻值，为热扩散.系数，为时间步长。由元体能量平衡法可以分别求出左右边界节点的离散方程,左边界为：00(?h2 )T。一T1h2Tf2ap0_ T 0（18）2xx2右边界为：00ap(寸h1)T N 1-Tn 2 h1Tf1ap 丁01 N 1（19）2xx22程序清单序号程序名称程序功能计算结果（1）设计程序用数值方法求解 温度 表1、4,图5、8 改变空间步长：表2,图6 改变时间步长：表3,图7 改变导热系数：表5,图9 改变墙外换热系数：表6,图10 改变墙内换热系数：表7,图11 改变墙壁厚度：表8,图12（2）验证程序用分析解求解温 度表1,图52.1设计程序主要程序段（1）初始
28、时刻的输入dx=dlt/(N-1); 空间步长apO=dc*dx/dtao;系数 apOfor(i=0;i<N;i+)TOi=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw; 初始温度Twi=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw;Tw0i=Tr-h1*(Tf1-Tr)*(dlt-i*dx)/lw;/ 假设温度场边界条件及中间节点程序输入左边界节点（墙外）C0=lw/dx;B0=ap0/2.+h2+lw/dx;A0=0.0;f0=h2*Tf2+ap0*T00/2.0;/内部节点 for(i=1;i<N-1;i+)各系数Ai=lw/dx;Ci=le/dx;
29、Bi=Ai+Ci+ap0;fi=ap0*T0i;/右边界节点AN-1=lw/dx;BN-1=h1+lw/dx+ap0/2.;CN-1=0.0;fN-1=ap0*T0N-1/2.+h1*Tf1(3) 追赶法程序tdma(N,A,B,C,f,Tw);/计算迭代间的最大误差Er=0.0;for(i=0;i<N;i+)Err=fabs(Twi-Tw0i); 各节点的误差 if(Err>Er)Er=Err;/求各节点中的最大误差lter=lter+1; 迭代次数if(Er>Eps)/不满足精度，重新迭代计算for(i=0;i<N;i+)Tw0i=Twi;goto loop;2.2
30、模型验证程序：利用公式求分析解输出内部时刻各点温度for(j=1;Tao<Taomax*3600;j+)Tao=Tao+dtao;/计 算时刻Fo=lmda*Tao/(CP*dlt*dlt);for(i=0;i<N;i+)a=0.0;b=Pi/2.0;th=0.0;/初始值定义n=i*dx/dlt;doa0=a;bO=b;doc=(a0+b0)/2.;fa=aO*ta n(aO)-Bi;fc=c*ta n(c)-Bi;if(fa*fc<O.O)bO=c;elseaO=c;while(fabs(aO-bO)>1.e-5);u=(aO+bO)/2.; 特征根Cn=2*sin(u)/(u+cos(u)*sin(u);/ 系数Ex=exp(-u*u*Fo)*cos(u* n);th=th+C n*Ex;a=a+Pi;b=b+Pi;while (fabs(C n*Ex)>1.e-5);Ti=th*thO+Tf;热能与动力工程系计算传热学程序设计成绩考核表指标考核内容分值得分1模型、方法和计算结果的可靠性（含程序考核）402讨论、分析的充分性、详实性203报告格式的规范性204答辩情况105平时表现106合计100教师签字：