肋片散热分析—计算传热学课程设计

中国石油大学(华东)

储建学院热能与动力工程系

《计算传热学程序设计》

设计报告

学生姓名：龚波

学号：08123217

专业班级：热能与动力工程08－2班

指导教师：黄善波

2011年 7 月 5 日

1 设计题目

在工程实际中，往往需要增加（对流）传热量，应用比较广泛的较为有效的一种方法就是增加换热面积，即采用肋片—在材料消耗量增加较少的条件下能较多地增大换热面积。在一些换热设备中，肋片得到了广泛地应用，如制冷装置的冷凝器、散热器、空气加热器等等。

1.1 设计题目

某等截面圆柱形直肋，设肋端是绝热的。试分析在一定的金属消耗量下，为使肋片的散热量达到最大时所需要的肋片尺寸，并分析肋片的材料、表面传热系数对该尺寸的影响。

1.2 已知参数

为了求得数值结果和利用结果进行分析，现给定题目相关已知量，包括肋片材料导热系数λ=λοk(T)=400(1+0.0035T),肋基温度T w=95℃，肋表度黑度ε=0.80，周围空气温度T f=20℃，环境辐射温度T s=15℃，肋表面空气的表面换热系数h c=8W/(m2?℃)。

2 物理与数学模型

2.1 物理模型

发生在肋片的导热过程严格地说是多维的。如图1所示，暴露于恒温流体的圆柱肋片（肋高为L，直径为D）。由于圆柱直肋各处受热均匀，再加上肋片通常是由金属材料制成的，导热系数比较大，可以想象肋片内温度将仅沿肋高方向发生明显变化，再直径方向上变化相比很小。因此，假设该圆柱直肋在同一截面上温度相同，则该问题可转化为等截面直肋一维稳态导热问题。

t/d x=0

图1 圆柱肋片物理模型图

2.2 数学模型

以肋基为坐标原点，圆柱肋片厚度方向为坐标正方向，建立坐标系如图2所示。 基于上述物理模型，则该问题的数学模型可描述如下： ()()44

0c f b s d dT A

U h T T T T dx dx λεσ????--+-= ?????

（1-a ）

左右两侧相应的边界条件分别是第一类边界条件和第二类边界条件，分别描述如下：

左边界

w x T

T == （1-b ）

右边界

0x L

dT dx

== （1-c

）

图2 圆柱肋片数学模型图

3 数值处理与程序设计

3.1数学模型无量纲化

为了使数值计算结果具有更普遍的意义，将上述数学模型无量纲化。为此定义

x x L

=

，

f w f

T T T T θ-=- （2）

控制方程无量纲化后，方程整理为

()

4

4

2

2

2

3

2

10f s

c b w f w f w

f T T

d dk d U L h T T k d A T T T T

d x d x

θθθεσθθλ?

????????

?

????+-?+-+

-= ? ??? ? ? ?--???????

???

???

?

（3）

定义 1d dk k k d θ

=

，L x U L S A

=

，c c h L

Bi ο

λ=

，

()

3

b w f

R T T

L

N ο

εσλ-=

，

f f w f

T T T θ=-，

s s w f

T T T θ=- （4）

将上述定义带入式（3）中，整理得：

()

{

}

2

2

4

4

20Lx d c R f

s S d d k Bi N k d x d x θ

θθθθθ??

??+-++-= ?????

?? （5-a ） 左边界

1x θ

== （5-b)

右边界

1

0x d d x

θ==

（5-c ）

3.2 试射法的形式

令

1y θ=，

12dy d y d x

d x

θ== （6）

则有试射法形式模型

12dy y d x

= （7-a ）

()

{

}

4

2

4

2211Lx d c R f

s S dy k y Bi y N y k

d x

θθ??=-+

++-???

?

(7-b)

左边界

1112

1P y Q y W += （7-c ）

其中，P 1=1,Q 1=0,W 1=1

右边界

2122

2P y Q y W += （7-d ）

其中，P 2=0,Q 2=1,W 2=0

3.3 程序编写

圆柱直肋一维稳态导热数学模型是二阶常微分两点边值问题，可以采用试射法求解。其基本思想是将边值问题转换为初值问题求解。

3.3.1 设计特点

在主程序外设置全局变量，为使在调用各子程序时，不会因实参与形参的作用范围而无法编译、运行程序。

在主程序头部，对参数赋值，对体积和肋高赋值应注意范围和两者的关联性。此处赋值V=0.00002m3，L=0.5m，保证程序结果为最大传热量，而且保证了足够的计算空间又不至于过分浪费系统资源。

利用循环实现计算最大传热量的过程，首先调用肋高函数得到按线性规律递减的肋高，再调用shoot函数计算相应肋高时的肋基温度梯度，调用热量函数求解热量Q[g]，输出各个肋高下的肋基温度梯度和热量，为了便于了解热量随肋高的变化关系。比较各肋高下的热量值，将最大热量值对应下标保留。然后，输出最大热量Q[max]和相应的肋高LG[max]，再根据几何关系求解圆柱肋片的面积A，半径r和此时的最佳长径比CJB （肋高与半径的比值）。再次调用shoot函数，求解最大传热量时圆柱肋片的温度分布和温度梯度。

求出最大传热量使用后，对程序进行验证，用户只需根据实际情况对热量函数RL，用户子程序的相关参数进行设置，不需要对验证程序进行操作，即可对程序结果进行验证。本程序在无辐射和导热率为定值时，即λ=C（常数），N R=0时，验证程序自动执行。

本程序采用的试射法考虑了物性的变化，辐射的影响，且对模型进行了无量纲化，因此具有普遍的适用性。

3.3.2 程序流程

先给程序中相关参数赋值，给定材料体积，利用试射法计算各个肋高是的肋基温度和温度梯度，根据温度梯度求肋片相应肋高的传热量，比较各个传热量值确定最大传热量，最后输出最大传热量对应结果，如果初参数满足验证程序的条件，执行验证程序并输出验证程序的结果，程序结束。程序流程图如下。

图3 程序流程图

4 模型与程序的验证

4.1 模型验证

为了方便利用解析解验证程序，将本题简化为常物性、无辐射等截面直肋一维稳态导热模型。已知肋片材料导热系数λ=100 W/(m ?℃),肋基温度T w =95℃，周围空气温度T f =20℃，肋表面空气的表面换热系数h c =8W/(m 2?℃)。建立坐标系，列出其控制方程式及定解条件：

2

2

2d m dx

θθ= (8-a)

00,;θθ===-w f x t t ,

0θ==d x H dx

(8-b)

其中过余温度f t t θ=-，

m =

式（8-a ）是一个二阶线性微分方程，由两边界条件可求出精确解为

()()

ch m x h ch m h θθ-????= （9）

4.2 程序验证

将式（9）中参数换算成无量纲形式，然后编程，计算出每个节点温度的解析解（验证程序见附录）和数值解（验证源程序见附录），进行比较，如表格1。

表1 λ=100、无辐射圆柱直肋无量纲温度值数值解和分析解

x 数值解 理论解 百分误差/%

0 1 1 0 0.05

0.966899 0.966899 0

0.1 0.935906 0.935905 0.00010685

0.15 0.906951 0.906951 0 0.2 0.879973 0.879973 0 0.25 0.854913 0.854912 0.00011697

0.3 0.831715 0.831715 0 0.35 0.81033 0.810329 0.00012341 0.4 0.79071 0.79071 0 0.45 0.772814 0.772813 0.0001294 0.5 0.756602 0.756601 0.00013217 0.55 0.742038 0.742037 0.00013476 0.6 0.729091 0.72909 0.00013716 0.65 0.717733 0.717732 0.00013933 0.7 0.707939 0.707938 0.00014126 0.75 0.699687 0.699686 0.00014292 0.8 0.69296 0.692959 0.00014431 0.85 0.687744 0.687742 0.00029081 0.9 0.684026 0.684024 0.00029239 0.95 0.681798 0.681797 0.00014667 1

0.681056

0.681055

0.00014683

00.20.40.60.81

1.20

0.3

0.60.9

1.2

无量纲位置

无量纲温度

图4 圆柱直肋无量纲温度分布曲线

由上述图表可知圆柱肋片分析解和数值解相差不大，二者吻合较好，可以说明所编

制的数值解法的程序是正确的。

5 计算结果与分析

5.1 肋高与热量的关系

材料的导热率λ=400(1+0.0035T)，肋基温度T w=95℃，肋表度黑度ε=0.80，周围空气温度T f=20℃，环境辐射温度T s=15℃，肋表面空气的表面换热系数h c=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

从表中结果易看出在肋高为0.50m和0.49m的时候，肋基温度梯度为正值，传热量为负值，与实际情况不符。这是因为肋高增加，一定的耗材下，肋片直径变小，对流换热量处理成广义热源已不合适，即不能作为一维稳态导热模型看待。此时，增大肋片体积、增加导热率或减小对流换热系数，又能满足模型使用调节。忽略表中前两行的数据绘图见图5。

图5 不同长度下肋片传热量曲线

由图可知传热量随着肋高先增后减，传热量最大在肋高L=0.23m取得。因为Φ=-λAy0[1]，在肋基，温度始终为t w，即导热系数不变，肋基温度梯度y0[1]为负且和圆柱截面积A随肋高增加而变小，所以存在最佳肋高使传热量最大。

5.2 表面换热系数的影响

材料的导热率λ=400(1+0.0035T)，肋表度黑度ε=0.80，圆柱肋片在不同表面换热系数h时，为使传热量最大，相应最佳肋高LG[max]和最佳长径比CJB（肋高与半径的比值）见表3。

表3 不同表面传热量下的LG[max]和CJB

由表易知，随着表面传热系数的增加，最佳肋高是逐渐减小的。表面传热系数的增加，传热量增加，由Φ=-λAy0 [1]知，需要增大肋基导热面积，所以最佳肋高减小。5.3 材料导热率的影响

无辐射，肋表面空气的表面换热系数h c=8W/(m2?℃)时，不同导热系数λ时，为使传热量最大，相应最佳肋高LG[max]和最佳长径比CJB（肋高与半径的比值）见表4。

表4 不同传热系数下的LG[max]和CJB

由表易知，随着导热系数的增加，最佳肋高是逐渐增大的。因为导热系数变大，传热量增加，由Φ=h c A(t-t f )知，需增加圆柱侧面积以加强换热，所以最佳肋高增加。

6 结论

在肋基，温度始终维持不变，即导热系数不变，肋基温度梯度为负且和圆柱截面积A 随肋高增加而变小，由傅里叶公式可知存在最佳肋高使肋片传热量最大，在题目已知条件下，当肋高

L=0.23m 时取得最大散热量Φ=15.0432W ；表面传热系数的增加，传热量增加，由傅里叶公式知，需要增大肋基导热面积，所以最佳肋高减小；导热系数变大，传热量增加，由对流换热公式知，需增加圆柱侧面积以加强换热，所以最佳肋高增加。

参考文献

[1] 黄善波，刘中良.计算传热学基础.中国石油大学（华东）热能与动力工程系，2009 [2] 杨世铭，陶文铨.传热学（第四版）.高等教育出版社，2007

附录1 主要程序

表5 程序列表

已知参数赋值

//输入圆柱肋体积V 、任意给定肋高L 及肋高变化步长bc V=0.00002; L=0.5;

bc=0.01;

void fct(int N,double x,double y[],double f[]) //函数子程序，用户根据具体条件进行修改

{

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.868 667,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;

}

求最大传热量

//调用各个函数求最大传热量

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("%f ",LG[g]); //输出长度

fprintf(fp,"%f ",LG[g]);

printf("%f ",y0[1]); //输出肋基出温度梯度 fprintf(fp,"%f ",y0[1]);

printf("%6.4f \n",Q[g]); //输出对应传热量

fprintf(fp,"%6.4f \n",Q[g]);

//保存最大传热量

if(Q[k]

}

//输出最大传热量，并输出对应圆肋的尺寸

max=k;

LG[max]=LeiGao(L,bc,max);

LGZ=LG[max];

A=V/LG[max];

r=sqrt(A/3.14);

CJB=LG[max]/r;

printf("最佳肋高LG[max]=%f\n",LGZ);

fprintf(fp,"最佳肋高LG[max]=%f\n",LGZ);

printf("最佳面积A=%f\n",A);

fprintf(fp,"最佳面积A=%f\n",A);

printf("最佳半径r=%f\n",r);

fprintf(fp,"最佳半径r=%f\n",r);

printf("最佳长径比CJB=%f\n",CJB);

fprintf(fp,"最佳长径比CJB=%f\n",CJB);

x=xa;

//利用试射法确定m

shoot(N,NS,x,h,M1,M2,Eps,y0,y);

Q[max]=RL(V, y0);//求最大传热量

printf("最大热量Q[max]=%6.4f ",Q[max]);

fprintf(fp,"最大热量Q[max]=%6.4f ",Q[max]);

printf("\n\n\n");

fprintf(fp,"\n\n\n");

//输出最大传热量时的温度分分布

printf("输出最大传热量时的温度分布\n");//显示在屏幕上

fprintf(fp,"输出最大传热量时的温度分布\n");//保存到文件中

//输出表头

printf(" x y1 y2\n");//显示在屏幕上

fprintf(fp," x y1 y2 \n");//保存到文件中 //输出x=a时的结果

printf("%6.4f ",xa);

fprintf(fp,"%6.4f ",xa);

for(i=0;i

{

printf("%10.6f ",y0[i]);

fprintf(fp,"%10.6f ",y0[i]);

}

printf("\n");

fprintf(fp,"\n");

x=xa;

//调用R-K方法计算并输出后续各点的值

for(i=0;i

y[i]=y0[i];

for(j=0;j

{

rungek(N,&x,h,y); //根据求出的m解决问题

printf("%6.4f ",x);

fprintf(fp,"%6.4f ",x);

for(i=0;i

{

printf("%9.6f ",y[i]);

fprintf(fp,"%9.6f ",y[i]);

}

printf("\n");

fprintf(fp,"\n");

}

验证程序

//令lmd=常数，无辐射，验证程序的正确性

if(kd==0.0&&Nr==0.0)

{

printf("\n\n令lmd=100，无辐射，验证程序理论温度\n");

fprintf(fp,"\n\n令lmd=100，无辐射，验证程序理论温度\n");

MM=sqrt(0.16/r);

MMH=MM*LGZ;

printf(" x y1 y2\n");//显示在屏幕上 fprintf(fp," x y1 y2 \n");//保存到文件中 x=0.0000;

for(i=0;h*i<=1.0;i++) { x=h*i;

printf("%6.4f ",x); fprintf(fp,"%6.4f ",x);

LLWD=(exp(MMH*(x-1))+exp(MMH*(1-x)))/(exp(MMH)+exp(-1.0*MMH));//温度分布

printf("%9.6f ",LLWD); fprintf(fp,"%9.6f ",LLWD); printf("\n"); fprintf(fp,"\n"); }

}

附录2 数学模型的无量纲化过程推导

针对式（1）进行无量纲化处理，为此定义

x x L

=

，

f w f

T T T T θ-=- （10）

其中Tf 、Tw 均为常数，假定λ=λοk(T)，则无量纲化过程如下。

2

2d dT d T d dT A A A dx dx dx dx dx λλλ??=+ ??? 2

2

2

d T d dT A A dx dx dx λλ??

=+ ???

（11）

()()

w f

f

w f d T T

T

T T dT d dx

L

d x

d xL

θθ??-+-?

?

=

=

（12-a ）

2

2

22

2

w f T T d T d dT d dx

dx dx L d x

θ-??== ??? (12-b)

()()w f w f f d k d dk

dx

T T d d T T T οο

λλλθ

θ=

=

-??-+??

（12-c ）

()

2

2

2

2

2

2

w f

w f w f T T T T A d dk d A L

T T d L

d x d x

ο

λθθλ

θ

--??

+

?-??

()(){}

4

4

0c w f b w f f s U h T T T T T T θεσθ??--+-+-=??

(13)

整理

()

4

4

2

2

2

3

2

10f

s

c b w f w f w

f T T

d dk d U L h T T k d A T T T T

d x d x θθθεσθθλ?

????????

?

????+-?+-+

-= ? ??? ? ? ?--???????

???

???

?

(14) 定义 1d dk k k d θ

=

，L x U L S A

=

，c c h L

Bi ο

λ=

，

()

3

b w f

R T T

L

N ο

εσλ-=

，f f w f

T T T θ=

-，s s w f

T T T θ=

- (15)

将上述定义带入式（3）中，整理得：

()

{

}

2

2

4

4

20Lx d c R f

s S d d k Bi N k d x d x θ

θθθθθ??

??+-++-= ?????

?? (16-a) 左边界

1x θ

== (16-b)

右边界

1

0x d d x

θ== (16c)

热能与动力工程系

《计算传热学程序设计》成绩考核表

教师签字：——此页单独占一页！