肋片散热数值计算2016年12月目录一、题目------------------------------------------3二、数值计算--------------------------------------4 （1）网格划分-----------------------------------4 （2）节点方程-----------------------------------5 （3）计算方式-----------------------------------6 （4）计算结果-----------------------------------6 （5）温度分布云图-------------------------------7 （6）误差分析-----------------------------------10三、结论------------------------------------------10四、程序------------------------------------------11五、参考文献--------------------------------------15一、题目肋片优化问题考虑三种不同形状的肋片，如图所示。
材料均为硬铝，热导率为，肋根半厚度为4mm，肋高为25mm。
对于梯形肋和圆弧边肋，最右端的平面部分半厚度为1mm，且圆弧在最右端的切线为水平线。
肋根温度即227℃，肋外流体温度即27℃，表面对流换热系数为。
试编程求解每种肋片的温度分布及散热量，并讨论肋片形状对散热量、材料需求量的优劣。
散热量17115 W 15605 W 14726 W截面面积200 mm2 125 mm2 101 mm2q/A 85.6 124.8 145.8肋顶温度367 K（94℃）355 K（82℃）342K（69℃）二、数值计算鉴于肋片对称，因此只研究上半部分即可，肋半厚度处可以按绝热处理。
（1）网格划分如图所示，用方格逼近曲面边界（三种形状的肋片都这样处理），即用图中蓝色网格边界来替代红色实际边界。
网格在边界处的取法用以下规则：对于第i列节点，如果其中的第j个节点在红线上方或者恰好在红线上，并且第j-1个节点在红线下方或者正好在红线上，那么就取第j个节点作为第i列上的边界节点，最后把所有边界节点连成锯齿状得到用方格逼近的边界（蓝色边界）。
网格间距delta在运行程序时输入。
（建议值0.1mm，如果输入的值过大，那么计算结果不精确；输入值过小则程序需要很长时间才能算出结果）（2）、节点方程内部节点肋根换热面传热量Q 按照肋根部的导热量计算，忽略根部y方向上的温度梯度（根据程序的运行结果来看，这种近似是合理的，对结果的影响很小），由求得的温度数据，在根部用(t(j,1)-t(j,2))/delta作为温度梯度，（由于程序的原因，这里t(j,1)表示第1列第j个温度值）再用温度梯度乘上微元面积delta*1（为计算方便，肋宽取为1m），然后把所有的值求和，就得到从肋根部导出的热量。
因为只研究上半部分，求出的散热量只有一半，再乘2即可。
具体到程序中所使用的节点方程，一共有一下11类（3）、计算方式简单迭代法，矩形肋片允许误差取为10^（-6），梯形肋片和圆弧形肋片的允许误差取为10^（-4）。
（4）、计算结果肋片形式允许误差e网格间距散热量Q 肋端温度截面面积AQ/A矩形肋片10^（-6）0.1mm 17394 96.2 200 86.97梯形肋片10^（-4）0.1mm 15107 75.6 125 199.93圆弧形肋10^（-4）0.1mm 14758 65.1 101 146.12 片（5）、温度分布云图矩形肋片温度分布梯形肋片温度分布圆弧形肋片温度分布（6）、误差分析由于用矩形边界代替曲面边界，导致截面周长增加，所以在程序中引入修正项L/L0，其中L为实际对流换热边界长度，L0为网格对流换热边界长度，将求得的热量乘上该修正项，能够减小误差。
修正步骤已经写入程序。
三、结论经过反复测试，发现对于矩形肋片，只要设置的允许误差e足够小，迭代求得的肋端温度值基本不受网格间距大小的影响，但是求得的散热量对网格间距比较敏感，当网格间距设为0.1mm，允许误差e=0.000001时，结果已经基本准确。
对于变截面肋片，即梯形和圆弧形肋片，测试发现求得的散热量Q和肋端温度对e和间距delta都很敏感；发现设置同样的e值，即e=0.0001，delta=0.1mm时，求得肋端温度75.7℃；delta=0.05mm 时，求得的肋端温度反而变小了，68℃。
但是可以预见，当e和delta取的都足够小时，求得的散热量和肋端温度将会越来越接近，但是鉴于将e和delta同时取很小，程序将运行很长时间，所以不再尝试。
对比计算结果可见，其实三种截面的肋片散热量差别不是很大，但是其Q/A值却差很大，也就是说矩形截面肋片用料比较多，而圆弧形、梯形肋片比较省材。
但是也看到，虽然圆弧形肋片省材，但其肋端温度较低，也就是说肋端的散热温差较小，这其实也是一种材料利用不充分的体现。
本题目属于给定换热系数h和导热系数λ的情况，而我们知道，肋片效率ηf=th(ml)*m,m=,因此肋片的效率取决于肋片的周长与面积比。
四、程序矩形肋片程序juxingleipian.mfunction juxingleipianh=2800;%对流换热系数lamda=187;%导热系数sprintf('输入网格间距delta，单位为mm,建议值0.1或0.2(其中0.2运行时间较短)') delta=input('delta=');%输入网格间距，单位mmdelta=0.001*delta;xnum=0.025/delta;%x轴划分数ynum=0.008/delta;%y轴划分数tf=300-273;%流体温度t0=500-273;%肋根温度x=zeros(ynum+1,xnum+1);t=zeros(ynum+1,xnum+1);w=zeros(ynum+1,xnum+1);c=0; sprintf('输入迭代允许误差e，(若delta输入0.2，则e的建议值为10^-6)当两次迭代之间的误差小于该数时，停止迭代')e=input('e=');%输入允许误差for i=1:ynum+1;j=1:xnum+1;t(i,j)=100;end;%任意假定一组初始温度值y=1;while y==1forj=2:xnum;x(1,j)=(2*t(2,j)+t(1,j+1)+t(1,j-1)+2*h*delta*tf/lamda)/(4+2*h*delta/lamda);end; %节点方程forj=2:xnum;x(ynum+1,j)=(2*t(ynum,j)+t(ynum+1,j+1)+t(ynum+1,j-1)+2*h*delta*tf/lamda) /(4+2*h*delta/lamda);end;%节点方程for i=1:ynum+1;x(i,1)=t0;end;%节点方程fori=2:ynum;x(i,xnum+1)=(2*t(i,xnum)+t(i+1,xnum+1)+t(i-1,xnum+1)+2*h*delta*tf/lamda)/ (4+2*h*delta/lamda);end%节点方程for i=2:ynum;j=2:xnum;x(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j))/4;end;%节点方程x(1,xnum+1)=(2*h*delta*tf/lamda+t(1,xnum)+t(2,xnum+1))/(2+2*h*delta/lamda);%节点方程x(ynum+1,xnum+1)=(2*h*delta*tf/lamda+t(ynum,xnum+1)+t(ynum+1,xnum))/(2+2*h*d elta/lamda);%节点方程for i=1:ynum+1;j=1:xnum+1;w(i,j)=abs(x(i,j)-t(i,j));endif (max(max(w))<=e)%判断两次迭代的误差是否小于允许值y=0;endt=x;c=c+1;endfor i=1:ynum+1;tidu(i)=(t(i,1)-t(i,2))/delta;end;%求肋根部温度梯度Q=lamda*sum(tidu(1,:))*delta%温度梯度与微元面积乘积求和，得到散热量temp=t(ynum/2+1,xnum+1)a=linspace(0,25,xnum+1);b=linspace(0,8,ynum+1);[aa,bb]=meshgrid(a,b);figuremesh(aa,bb,t);%温度分布图figurecontourf(a,b,t,50);shading flatend梯形肋片&&圆弧形肋片程序laddershaped.m注：本程序改变红色字的公式即可分别计算梯形肋片和圆弧形肋片function laddershapedh=2800;%对流换热系数r=634/6;%半径rc=0;lamda=187;%导热系数tf=300-273;%流体温度t0=500-273;%肋根温度sprintf('输入网格间距delta，单位为mm')delta=input('delta=');%输入网格间距，单位mmsprintf('输入设定的误差值，当两次迭代的误差小于该值时停止迭代');wucha=input('wucha=');%输入设定的误差值delta=0.001*delta;xnum=0.025/delta;%x轴划分数ynum=0.008/delta;%y轴划分数t=zeros(ynum/2+1,xnum+1);%设初始温度场为100℃T=zeros(ynum/2+1,xnum+1);%下面要确定第i列对应的行数i=1;while (i<=xnum+1)%求出每一列的数据个数y(i)=(-3*(i-1)*delta*1000+100)/25;%梯形肋片；肋片形状可以任意更改，只需写出不同截面的方程即可%y(i)=((- 3*(i-1)*delta*1000 - 242)^(1/2)*(3*(i-1)*delta*1000 - 392)^(1/2))/3 + 320/3;%圆弧形肋片for j=1:ynum/2+1;y1(j)=(j-1)*delta*1000-y(i);y2(j)=y1(j)-delta*1000;if(y1(j)>=0 && y2(j)<=0)z(i)=j;endendi=i+1;%求下一列数据的个数end%下面是节点方程p=1;while(p==1)%从这里开始循环for j=1:ynum/2+1;T(j,1)=t0;end;%根部边界条件，第1类for i=2:xnum;T(1,i)=(t(1,i-1)+t(1,i+1)+2*t(2,i))/4;end;%第2类T(1,xnum+1)=(t(1,xnum)+t(2,xnum+1)+h*delta*0.5*tf/lamda)/(h*delta*0.5/lamda+2);%第3类forj=2:0.001/delta;i=xnum+1;T(j,i)=(2*t(j,i-1)+t(j+1,i)+t(j-1,i)+2*h*delta*tf/lamda)/(h*delta* 2/lamda+4);end;%第12类i=2;while (i>=2 && i<=xnum)if(z(i)==z(i-1) && z(i)==z(i+1))T(z(i),i)=(2*t(z(i)-1,i)+t(z(i),i+1)+t(z(i),i-1)+2*h*delta*tf/lamda)/(h*delta*2/lamda+4);%第7类else if(z(i)==z(i-1) && z(i)>z(i+1))T(z(i),i)=(t(z(i)-1,i)+t(z(i),i-1)+2*h*delta*tf/lamda)/(2+2*h*delta/lamda);%第6类T(z(i)-1,i)=(t(z(i),i)+t(z(i)-1,i+1)+2*(t(z(i)-2,i)+t(z(i)-1,i-1))+2*h*delta*tf/lamda)/(h*delta* 2/lamda+6);%第5类else if (z(i)==z(i+1) && z(i-1)>z(i))T(z(i),i)=(2*t(z(i)-1,i)+t(z(i),i+1)+t(z(i),i-1)+2*h*delta*tf/lamda)/(h*delta*2/lamda+4);%第8类endendendi=i+1;endi=xnum+1;T(z(i),i)=(t(z(i)-1,i)+t(z(i),i-1)+2*h*delta*tf/lamda)/(2+2*h*delta/lamda);%第4类i=2;j=2;k=2;l=2;while(i<=xnum)%纯导热部分的节点方程if(z(i)==z(i-1) && z(i)==z(i+1))while(j<z(i))T(j,i)=0.25*(t(j+1,i)+t(j-1,i)+t(j,i+1)+t(j,i-1));%第9类j=j+1;end;else if(z(i)==z(i-1) && z(i)>z(i+1))while(k<z(i)-1)T(k,i)=0.25*(t(k+1,i)+t(k-1,i)+t(k,i+1)+t(k,i-1));%第10类k=k+1;end;else if(z(i)==z(i+1) && z(i-1)>z(i))while(l<z(i))T(l,i)=0.25*(t(l+1,i)+t(l-1,i)+t(l,i+1)+t(l,i-1));%第11类l=l+1;end;endendendi=i+1;j=2;k=2;l=2;endi=1;j=1;while(i<=xnum+1)while(j<=z(i))e(j,i)=abs(T(j,i)-t(j,i));j=j+1;endi=i+1;j=1;endt=T;%把新的一组温度值赋值给tif max(max(e))<wucha;p=0;end;c=c+1;endfor i=1:ynum/2+1;tidu(i)=(t(i,1)-t(i,2))/delta;end;%求肋根部温度梯度Q=lamda*sum(tidu(1,:))*delta%温度梯度与微元面积乘积求和，得到散热量%Qreal=2*Q/(25+3+1)*(25^2+3^2)^0.5;%梯形肋片的散热量Qreal=2*Q*25.54859743/(25+3+1);%圆弧形肋片的散热量temp=t(1,xnum+1); %求肋端部温度tempQreala=linspace(0,25,xnum+1);b=linspace(0,4,ynum/2+1);[aa,bb]=meshgrid(a,b);figuremesh(aa,bb,t)%温度分布图figurecontourf(a,b,t,50)shading flatsprintf('参考值--矩形肋片：17115W，94；梯形肋片：15605W，82度；圆弧形肋片：14726W，69度')end五、参考文献《传热学（第五版）》.章熙民，任泽霈，梅飞鸣. 中国建筑工业出版社.。
