中国石油大学(华东)
储建学院热能与动力工程系
《计算传热学程序设计》
设计报告
学生姓名：龚波 学 号：08123217
专业班级：热能与动力工程08－2班 指导教师：黄善波
2011年 7 月 5 日
1 设计题目
在工程实际中，往往需要增加（对流）传热量，应用比较广泛的较为有效的一种方法就是增加换热面积，即采用肋片—在材料消耗量增加较少的条件下能较多地增大换热面积。在一些换热设备中，肋片得到了广泛地应用，如制冷装置的冷凝器、散热器、空气加热器等等。 1.1 设计题目
某等截面圆柱形直肋，设肋端是绝热的。试分析在一定的金属消耗量下，为使肋片的散热量达到最大时所需要的肋片尺寸，并分析肋片的材料、表面传热系数对该尺寸的影响。 1.2 已知参数
为了求得数值结果和利用结果进行分析，现给定题目相关已知量，包括肋片材料导热系数λ=λοk(T)=400(1+0.0035T),肋基温度Tw=95℃，肋表度黑度ε=0.80，周围空气温度Tf=20℃，环境辐射温度Ts=15℃，肋表面空气的表面换热系数hc=8W/(m2?℃)。
2 物理与数学模型
2.1 物理模型
发生在肋片的导热过程严格地说是多维的。如图1所示，暴露于恒温流体的圆柱肋片（肋高为L，直径为D）。由于圆柱直肋各处受热均匀，再加上肋片通常是由金属材料制成的，导热系数比较大，可以想象肋片内温度将仅沿肋高方向发生明显变化，再直径方向上变化相比很小。因此，假设该圆柱直肋在同一截面上温度相同，则该问题可转化为等截面直肋一维稳态导热问题。
tfεtsDtwdt/dx=0tfH
图1 圆柱肋片物理模型图
2.2 数学模型
1
以肋基为坐标原点，圆柱肋片厚度方向为坐标正方向，建立坐标系如图2所示。 基于上述物理模型，则该问题的数学模型可描述如下： Ad?dT?44hc?T?Tf????b?T?Ts???0 ????U???dx?dx? （1-a）
左右两侧相应的边界条件分别是第一类边界条件和第二类边界条件，分别描述如下：
左边界
T右边界
x?0?Tw （1-b）
dTdxx?L?0 （1-c）
图2 圆柱肋片数学模型图
3 数值处理与程序设计
3.1数学模型无量纲化
为了使数值计算结果具有更普遍的意义，将上述数学模型无量纲化。为此定义
x?xL??，
T?TfTw?Tf （2）
控制方程无量纲化后，方程整理为
?d?1dk?d??UL?????hc????b?Tw?Tf??2kd??dx?A??dx?222?3??Tf????Tw?Tf??????Ts?????T?Tf??w4????4??????0 （3） ????2
定义 kd?1dkkd?，SLx?ULA，Bic?hcL?? ，
NR???b?Tw?T???f3L?f?，
TfTw?Tf?s?，
TsTw?Tf （4）
将上述定义带入式（3）中，整理得：
SLx?d???k?Bic??NR?????fd??2??k?dx?dxd?22??44 ??s??0 （5-a）???左边界
?右边界
d?dxx?1x?0?1 （5-b)
?0
（5-c）
3.2 试射法的形式
令
??y1，
d?dx?dy1dx?y2 （6）
则有试射法形式模型
dy1dxdy2dx??kdy2?2?y2 （7-a）
SLxk?Bicy1?NR??y1??f???44??s? (7-b)
???左边界
P1y1?Q1y?2W 1 （7-c）
其中，P1=1,Q1=0,W1=1
右边界
P2y1?Q2y?2W 2 （7-d）
其中，P2=0,Q2=1,W2=0
3.3 程序编写
圆柱直肋一维稳态导热数学模型是二阶常微分两点边值问题，可以采用试射法求解。其基本思想是将边值问题转换为初值问题求解。
3
3.3.1 设计特点
在主程序外设置全局变量，为使在调用各子程序时，不会因实参与形参的作用范围而无法编译、运行程序。
在主程序头部，对参数赋值，对体积和肋高赋值应注意范围和两者的关联性。此处赋值V=0.00002m3，L=0.5m，保证程序结果为最大传热量，而且保证了足够的计算空间又不至于过分浪费系统资源。
利用循环实现计算最大传热量的过程，首先调用肋高函数得到按线性规律递减的肋高，再调用shoot函数计算相应肋高时的肋基温度梯度，调用热量函数求解热量Q[g]，输出各个肋高下的肋基温度梯度和热量，为了便于了解热量随肋高的变化关系。比较各肋高下的热量值，将最大热量值对应下标保留。然后，输出最大热量Q[max]和相应的肋高LG[max]，再根据几何关系求解圆柱肋片的面积A，半径r和此时的最佳长径比CJB（肋高与半径的比值）。再次调用shoot函数，求解最大传热量时圆柱肋片的温度分布和温度梯度。
求出最大传热量使用后，对程序进行验证，用户只需根据实际情况对热量函数RL，用户子程序的相关参数进行设置，不需要对验证程序进行操作，即可对程序结果进行验证。本程序在无辐射和导热率为定值时，即λ=C（常数），NR=0时，验证程序自动执行。
本程序采用的试射法考虑了物性的变化，辐射的影响，且对模型进行了无量纲化，因此具有普遍的适用性。 3.3.2 程序流程
先给程序中相关参数赋值，给定材料体积，利用试射法计算各个肋高是的肋基温度和温度梯度，根据温度梯度求肋片相应肋高的传热量，比较各个传热量值确定最大传热量，最后输出最大传热量对应结果，如果初参数满足验证程序的条件，执行验证程序并输出验证程序的结果，程序结束。程序流程图如下。
4
开始 输出最大热量下的结果 给程序中相关参数赋值 保存结果至文件 给定材料体积 判断初参数 计算各个肋高下肋基温度和温度梯度 Y 执行验证程序并输出结果 N Y 程序结束 N 计算热量、比较并求出最大热量
图3 程序流程图
4 模型与程序的验证
4.1 模型验证
为了方便利用解析解验证程序，将本题简化为常物性、无辐射等截面直肋一维稳态导热模型。已知肋片材料导热系数λ=100 W/(m?℃),肋基温度Tw=95℃，周围空气温度Tf=20℃，肋表面空气的表面换热系数hc=8W/(m2?℃)。建立坐标系，列出其控制方程式及定解条件：
d?dx22?m? (8-a)
d?dx2x?0,???0?tw?tf; x?H,?0
(8-b)
其中过余温度??t?tf，m?hp/??Ac?为一常量。
式（8-a）是一个二阶线性微分方程，由两边界条件可求出精确解为
???0ch??m?x?h???ch?mh? （9）
4.2 程序验证
5
将式（9）中参数换算成无量纲形式，然后编程，计算出每个节点温度的解析解（验证程序见附录）和数值解（验证源程序见附录），进行比较，如表格1。
表1 λ=100、无辐射圆柱直肋无量纲温度值数值解和分析解
x 0 0.05
0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1
数值解 1 0.966899 0.935906 0.906951 0.879973 0.854913 0.831715 0.81033 0.79071 0.772814 0.756602 0.742038 0.729091 0.717733 0.707939 0.699687 0.69296 0.687744 0.684026 0.681798 0.681056
理论解 1 0.966899 0.935905 0.906951 0.879973 0.854912 0.831715 0.810329 0.79071 0.772813 0.756601 0.742037 0.72909 0.717732 0.707938 0.699686 0.692959 0.687742 0.684024 0.681797 0.681055
百分误差/%
0 0
0.00010685
0 0 0.00011697
0 0.00012341 0 0.0001294 0.00013217 0.00013476 0.00013716 0.00013933 0.00014126 0.00014292 0.00014431 0.00029081 0.00029239 0.00014667 0.00014683
1.21无量纲温度0.80.60.40.2000.30.6无量纲位置0.91.2
图4 圆柱直肋无量纲温度分布曲线
数值解理论解由上述图表可知圆柱肋片分析解和数值解相差不大，二者吻合较好，可以说明所编
6
制的数值解法的程序是正确的。
5 计算结果与分析
5.1 肋高与热量的关系
材料的导热率λ=400(1+0.0035T)，肋基温度Tw=95℃，肋表度黑度?=0.80，周围空气温度Tf=20℃，环境辐射温度Ts=15℃，肋表面空气的表面换热系数hc=8W/(m2?℃)时，圆柱肋片肋基无量纲温度梯度和传热量见表2。
表2 不同长度下肋片的传热量 任意长度L/m 肋基温度梯度 热量?/W
0.5 10 -54.9246 0.49 10 -57.1893 0.48 -1.98619 11.8371 0.47 -1.92857 11.988 0.46 -1.87095 12.141 0.45 -1.81331 12.2957 0.4 -1.52509 13.0883 0.39 -1.46747 13.2479 0.36 -1.29491 13.7196 0.35 -1.23757 13.8721 0.34 -1.18038 14.0208 0.33 -1.12339 14.1647 0.32 -1.06663 14.3028 0.27 -0.78888 14.859 0.26 -0.73514 14.9324 0.25 -0.68225 14.9889 0.24 -0.63034 15.0266 0.23 -0.57955 15.0432 0.22 -0.53001 15.0366 0.21 -0.48189 15.0044 0.2 -0.43533 14.9441 0.19 -0.39049 14.8529 0.18 -0.34756 14.7295 0.17 -0.30666 14.5704 0.1 -0.08934 12.2674 0.07 -0.03744 10.4921 0.06 -0.0256 9.7631 0.05 -0.01629 8.9465 0.04 -0.00935 8.0207 0.03 -0.00456 6.9549 0.02 -0.00166 5.6963 0.01
-0.00029
3.9285
7
从表中结果易看出在肋高为0.50m和0.49m的时候，肋基温度梯度为正值，传热量为负值，与实际情况不符。这是因为肋高增加，一定的耗材下，肋片直径变小，对流换热量处理成广义热源已不合适，即不能作为一维稳态导热模型看待。此时，增大肋片体积、增加导热率或减小对流换热系数，又能满足模型使用调节。忽略表中前两行的数据绘图见图5。
16传热量/W1412108642000.20.40.6实际肋高/m
图5 不同长度下肋片传热量曲线
由图可知传热量随着肋高先增后减，传热量最大在肋高L=0.23m取得。因为
?=-λAy0 [1]，在肋基，温度始终为tw，即导热系数不变，肋基温度梯度y0[1]为负且和
圆柱截面积A随肋高增加而变小，所以存在最佳肋高使传热量最大。 5.2 表面换热系数的影响
材料的导热率λ=400(1+0.0035T)，肋表度黑度ε=0.80，圆柱肋片在不同表面换热系数h时，为使传热量最大，相应最佳肋高LG[max]和最佳长径比CJB（肋高与半径的比值）见表3。
表3 不同表面传热量下的LG[max]和CJB
h/(W/m2?℃) LG[max]/m CJB 4 0.15 23.019014 5 0.14 20.755915 7 0.13 18.572264 9 0.13 18.572264 10 0.12 16.471066 12 0.12 16.471066 由表易知，随着表面传热系数的增加，最佳肋高是逐渐减小的。表面传热系数的增加，传热量增加，由?=-λAy0 [1]知，需要增大肋基导热面积，所以最佳肋高减小。 5.3 材料导热率的影响
无辐射，肋表面空气的表面换热系数hc=8W/(m2?℃)时，不同导热系数λ时，为使传热量最大，相应最佳肋高LG[max]和最佳长径比CJB（肋高与半径的比值）见表4。
8
表4 不同传热系数下的LG[max]和CJB
λ/(W/mk) LG[max]/m CJB 100 0.18 30.25928 300 0.28 58.706592 400 0.31 68.389963 600 0.36 85.586167 700 0.39 96.504316 800 0.41 104.022099 由表易知，随着导热系数的增加，最佳肋高是逐渐增大的。因为导热系数变大，传热量增加，由?=hcA(t-tf)知，需增加圆柱侧面积以加强换热，所以最佳肋高增加。
6 结论
在肋基，温度始终维持不变，即导热系数不变，肋基温度梯度为负且和圆柱截面积A随肋高增加而变小，由傅里叶公式可知存在最佳肋高使肋片传热量最大，在题目已知条件下，当肋高L=0.23m时取得最大散热量?=15.0432W；表面传热系数的增加，传热量增加，由傅里叶公式知，需要增大肋基导热面积，所以最佳肋高减小；导热系数变大，传热量增加，由对流换热公式知，需增加圆柱侧面积以加强换热，所以最佳肋高增加。
参考文献
[1] 黄善波，刘中良.计算传热学基础.中国石油大学（华东）热能与动力工程系，2009 [2] 杨世铭，陶文铨.传热学（第四版）.高等教育出版社，2007
附录1 主要程序
表5 程序列表
序号 1 程序名称 全功能程序 程序功能 对应图表 按题目要求，输出各个肋高下的传热量，求解最大传图3、4、5，表热量以及此时肋片的尺寸，输出温度和温度梯度的分布，验证程序 1、2、3、4 已知参数赋值
//输入圆柱肋体积V、任意给定肋高L及肋高变化步长bc V=0.00002; L=0.5;
bc=0.01;
void fct(int N,double x,double y[],double f[]) //函数子程序，用户根据具体条件进行修改
9
{
double Kd,Slx,Bic;
Kd=0.2625/(2.026025+0.2625*y[0]); Slx=2*LGZ*sqrt(3.14*LGZ/V); Bic=0.08*LGZ/4.0; Nr=0.000191*LGZ/4.0; f[0]=y[1]; //f[0]=dy1/dx
f[1]=-Kd*y[1]*y[1]+Slx*(Bic*y[0]+Nr*(pow((y[0]+3.908667),4)-pow(3.868667,4))); return; }
void pqw1(double Y,double *P,double *Q,double *W) //左边界处的第三类边界条件(x=xa)
//P1*y1+Q1*y2=W1-用户应根据具体条件进行修改 {
*P=1.0; *Q=0.0; *W=1; }
void pqw2(double Y,double *P,double *Q,double *W) //右边界处的第三类边界条件(x=xb)
//P2*y1+Q2*y2=W2-用户应根据具体条件进行修改 { *P=0; *Q=1; *W=0; return; } 求最大传热量
10
//调用各个函数求最大传热量 for(g=0;L-g*bc>0;g++) {
//用肋高函数求肋高
LG[k]=LeiGao(L,bc,k);
LGZ=LG[k];
x=xa;
//利用试射法确定m
shoot(N,NS,x,h,M1,M2,Eps,y0,y);
//求热量
Q[k]=RL(V, y0); LG[g]=LeiGao(L,bc,g);
LGZ=LG[g];
x=xa;
//利用试射法确定m
shoot(N,NS,x,h,M1,M2,Eps,y0,y); Q[g]=RL(V, y0);
//输出任意长度及所对应的热量
printf(\输出长度 fprintf(fp,\
printf(\输出肋基出温度梯度 fprintf(fp,\
printf(\输出对应传热量 fprintf(fp,\ //保存最大传热量
if(Q[k] } //输出最大传热量，并输出对应圆肋的尺寸 max=k; LG[max]=LeiGao(L,bc,max); 11 LGZ=LG[max]; A=V/LG[max]; r=sqrt(A/3.14); CJB=LG[max]/r; printf(\最佳肋高LG[max]=%f\\n\ fprintf(fp,\最佳肋高LG[max]=%f\\n\ printf(\最佳面积A=%f\\n\ fprintf(fp,\最佳面积A=%f\\n\ printf(\最佳半径r=%f\\n\ fprintf(fp,\最佳半径r=%f\\n\ printf(\最佳长径比CJB=%f\\n\ fprintf(fp,\最佳长径比CJB=%f\\n\ x=xa; //利用试射法确定m shoot(N,NS,x,h,M1,M2,Eps,y0,y); Q[max]=RL(V, y0);//求最大传热量 printf(\最大热量Q[max]=%6.4f \ fprintf(fp,\最大热量Q[max]=%6.4f \ printf(\ fprintf(fp,\ //输出最大传热量时的温度分分布 printf(\输出最大传热量时的温度分布\\n\显示在屏幕上 fprintf(fp,\输出最大传热量时的温度分布\\n\保存到文件中 //输出表头 printf(\显示在屏幕上 fprintf(fp,\保存到文件中 //输出x=a时的结果 printf(\ fprintf(fp,\ for(i=0;i 12 { printf(\ fprintf(fp,\ } printf(\ fprintf(fp,\ x=xa; //调用R-K方法计算并输出后续各点的值 for(i=0;i rungek(N,&x,h,y); //根据求出的m解决问题 printf(\ fprintf(fp,\ for(i=0;i printf(\ fprintf(fp,\ } printf(\ fprintf(fp,\ } 验证程序 //令lmd=常数，无辐射，验证程序的正确性 if(kd==0.0&&Nr==0.0) { printf(\令lmd=100，无辐射，验证程序理论温度\\n\ fprintf(fp,\令lmd=100，无辐射，验证程序理论温度\\n\ MM=sqrt(0.16/r); 13 MMH=MM*LGZ; printf(\显示在屏幕上 fprintf(fp,\保存到文件中 x=0.0000; for(i=0;h*i<=1.0;i++) { x=h*i; printf(\ fprintf(fp,\ LLWD=(exp(MMH*(x-1))+exp(MMH*(1-x)))/(exp(MMH)+exp(-1.0*MMH));//温度分布 printf(\ fprintf(fp,\ printf(\ fprintf(fp,\ } } 附录2 数学模型的无量纲化过程推导 针对式（1）进行无量纲化处理，为此定义 x?xL??， T?TfTw?Tf （10） 其中Tf、Tw均为常数，假定λ=λοk(T)，则无量纲化过程如下。 2d?dT?dTd?dTdTd??dT??A??A?A ?A??A???? （11）22dx?dx?dxdxdxdxdx?dx?22 dTdx?d???Tw?T?dxLf??Tf?????Tw?Tfd?Ldx （12-a） dTdx222d?dT?Tw?Tfd?????22dx?dx?Ldx (12-b) 14 d?dx?d???k?d???Tw?Tf??Tf???2???dkTw?Tfd? （12-c） 2A?Tw?Tfd?L2dx2?A??dk?Tw?TfL2?2Tw?Tfd??d?????dx?4 4?Uhc??Tw?Tf????b???Tw?Tf??Tf??Ts?0 ???? (13) 整理 ?d?1dk?d??UL?????hc????b?Tw?Tf??2kd?A?dx??dx??222?3??Tf????Tw?Tf??????Ts?????T?Tf??w4????4??????0 (14) ????定义 kd?1dkkd?，SLx?ULA，Bic?hcL?? ， NR???b?Tw?T???f3L，?f?TfTw?Tf，?s?TsTw?Tf (15) 将上述定义带入式（3）中，整理得： SLx?d???k?Bic??NR?????fd??2??k?dx?dxd?22??44??s??0 (16-a) ???左边界 ?右边界 d?dxx?1x?0?1 (16-b) ?0 (16c) 15 热能与动力工程系 《计算传热学程序设计》成绩考核表 指标 1 2 3 4 5 6 考核内容 模型、方法和计算结果的可靠性（含程序考核） 讨论、分析的充分性、详实性 报告格式的规范性 报告内容的阐述、问答问题的情况 平时的表现 合计 分值 30 25 20 10 15 100 得分 教师签字： ——此页单独占一页！ 16