=10
。
四:
编程思绪及步骤图
结束
否
图3-1
五、程序及运行结果
第三类边界条件:
1、试验程序(C语言):
//1.cpp:
定义控制台应用程序入口点。
//
#include"stdafx.h"
#include
#include
int_tmain(intargc,_TCHAR*argv[])
{
inti,j,l;
floatdt=1.0,dx=0.1,dy=0.1;
floatt[13][17],a[13][17];
floatq1=0,q2=0,q=0,e;
floatlmd=0.53,h1=10.33,h2=3.93,t1=30,t2=10,ep=1.0e-7;
/*打印出题目*/
printf("\t\t\t二维稳态导热问题\t\t");
printf("\n\t\t\t\t\t\t----陈振兴装备\n");
printf("\n题目:
二维导热物体温度场电模拟试验\n");
printf("\n矩形区域,l1=2.2;l2=3;l3=2;l4=1.2,假设区域内无内热源,导热系数为常熟,内外表面均为第三类边界条件且已知t1=30;t2=10;h1=10.33;h2=3.93;LMD=0.53;求该矩形区域内温度分布及垂直于纸面方向单位长度上经过墙体导热量。
\n");
/*各节点上温度值*/
{
for(j=0;j<17;j++)
{
t[0][j]=30.0;
}
for(i=1;i<13;i++)
{
t[i][0]=30.0;
}
for(i=7;i<13;i++)
{
t[i][7]=10.0;
}
for(j=8;j<17;j++)
{
t[7][j]=10.0;
}
for(i=1;i<7;i++)
for(j=1;j<17;j++)
t[i][j]=20,a[i][j]=0;
}
{
for(i=7;i<13;i++)
for(j=1;j<7;j++)
t[i][j]=20,a[i][j]=0;
}
while(dt>=ep)
{
{
for(i=1;i<7;i++)
for(j=1;j<17;j++)
a[i][j]=t[i][j];
}
{
for(i=7;i<13;i++)
for(j=1;j<7;j++)
a[i][j]=t[i][j];
}
{
for(i=6;i<12;i++)
for(j=2;j<6;j++)
t[i][j]=(t[i-1][j]+t[i+1][j]+t[i][j-1]+t[i][j+1])/4;
}
{
for(i=2;i<6;i++)
for(j=2;j<16;j++)
t[i][j]=(t[i-1][j]+t[i+1][j]+t[i][j-1]+t[i][j+1])/4;
}
{
for(j=2;j<6;j++)
t[12][j]=(t[12][j-1]+t[12][j+1]+2*t[11][j])/4;
}
{
for(i=2;i<6;i++)
t[i][16]=(t[i-1][16]+t[i+1][16]+2*t[i][15])/4;
}
{
for(i=2;i<12;i++)
t[i][1]=(dx*h1*t1/lmd+(t[i+1][1]+t[i-1][1])/2+t[i][2])/(2+dx*h1/lmd);
}
{
for(j=2;j<16;j++)
t[1][j]=(dy*h1*t1/lmd+(t[1][j+1]+t[1][j-1])/2+t[2][j])/(2+dy*h1/lmd);
}
{
for(i=7;i<12;i++)
t[i][6]=(dx*h2*t2/lmd+(t[i+1][6]+t[i-1][6])/2+t[i][5])/(2+dx*h2/lmd);
}
{
for(j=7;j<16;j++)
t[6][j]=(dy*h2*t2/lmd+(t[6][j+1]+t[6][j-1])/2+t[5][j])/(2+dy*h2/lmd);
}
t[1][1]=(h1*dx*t1+lmd*(t[2][1]+t[1][2])/2)/(lmd+h1*dy);
t[1][16]=(h1*dx*t1+lmd*(t[1][15]+t[2][16])/2)/(lmd+h1*dy);
t[6][16]=(h2*dx*t2+lmd*(t[6][15]+t[5][16])/2)/(lmd+h2*dy);
t[12][1]=(h1*dx*t1+lmd*(t[12][2]+t[11][1])/2)/(lmd+h1*dy);
t[12][6]=(h2*dx*t2+lmd*(t[12][5]+t[11][6])/2)/(lmd+h2*dy);
t[6][6]=(h2*dy*t2+lmd*(t[5][6]+t[6][5]+t[7][6]/2+t[6][7]/2))/(3*lmd+h2*dx);
{
for(i=1;i<7;i++)
{
for(j=1;j<17;j++)
dt=dt+abs(t[i][j]-a[i][j]);
dt=dt/(6*16);
}
for(i=7;i<13;i++)
{
for(j=1;j<7;j++)
dt=dt+abs(t[i][j]-a[i][j]);
dt=dt/(6*6);
}
}
}
printf("温度分布为:
\t\t\t\t\t\t\t\t\t");
l=0;
{
for(i=1;i<7;i++)
{
for(j=1;j<17;j++)
printf("%3.1f",t[i][j]);
l=l+1;
if(l==16)
{
printf("\n");
l=0;
}
}
l=0;
for(i=7;i<13;i++)
{
for(j=1;j<7;j++)
{printf("%3.1f",t[i][j]);
l=l+1;
if(l==6)
{
printf("\n");
l=0;
}
}
}
}
{
for(j=2;j<16;j++)
q1=q1+(30-t[1][j])*h1*dx;
for(i=2;i<12;i++)
q1=q1+(30-t[i][1])*h1*dy;
for(j=7;j<17;j++)
q2=q2+(t[6][j]-10)*h2*dx;
for(i=7;i<12;i++)
q2=q2+(t[i][6]-10)*h2*dy;
q1=q1+h1*(dx/2*(30-t[1][16])+dy/2*(30-t[12][1])+dx*(30-t[1][1]));
q2=q2+h2*(dx/2*(t[6][16]-10)+dy/2*(t[12][6]-10)+dx*(t[7][7]-10));
}
q=(q1+q2)/2;
e=abs((q2-q1)/q);
printf("单位长度上1/4墙体导热量为:
%4.2fW,偏差为:
%3.2f",q,e);
getchar();getchar();
return0;
}
运行结果图:
图3-2
试验算得导热量为97.62W,与数值模拟偏差为(26.73*4W-97.62W)/(26.73*4)W*100%=8.7%
2、数值模拟程序(matlab):
z=[29.929.729.529.329.129.028.828.728.728.628.628.528.528.528.528.2
29.729.128.528.027.426.926.526.226.025.925.825.725.625.625.625.5
29.528.527.626.625.724.824.123.623.323.122.922.822.822.722.722.7
29.328.026.625.223.822.521.620.920.520.220.019.919.919.819.819.8
29.127.425.723.821.919.918.717.917.517.217.117.016.916.916.916.9
29.026.924.822.519.916.715.214.614.314.214.114.014.014.014.013.9
28.826.524.121.518.615.20000000000;
28.726.223.520.817.814.60000000000;
28.625.923.220.317.314.20000000000;
28.525.722.820.017.014.00000000000;
28.325.422.619.716.713.60000000000;
27.225.122.519.516.312.50000000000];
v=[182226];
[xx,yy]=meshgrid(y,x);
surf(xx,yy,z);colorbar;xlabel('x');ylabel('y');zlabel('z');az=0;el=-90;view(az,el);
shadinginterp;
axistight;
figure,contour(xx,yy,z,v);
gridon
数值模拟图:
图3-3
图3-4
第一类边界条件
1、试验程序(C语言):
#include"stdafx.h"
#include
#include
int_tmain(intargc,_TCHAR*argv[])
{
inti,j,l;
floatdt=1.0,dx=0.1,dy=0.1;
floatt[12][16],a[12][16];
floatq1=0,q2=0,q=0,e;
floatlmd=0.53,t1=30,t2=0,ep=1.0e-7;
/*打印出题目*/
printf("\t\t\t二维稳态导热问题\t\t");
printf("\n\t\t\t\t\t\t----陈振兴装备\n");
printf("\n题目:
二维导热物体温度场电模拟试验\n");
printf("\n矩形区域,l1=2.2;l2=3;l3=2;l4=1.2,假设区域内无内热源,导热系数为常熟,内外表面均为第一类边界条件且已知t1=30;t2=0;LMD=0.53;求该矩形区域内温度分布及垂直于纸面方向单位长度上经过墙体导热量。
\n");
/*各节点上温度值*/
{
for(j=0;j<16;j++)
{
t[0][j]=30.0;
}
for(i=1;i<12;i++)
{
t[i][0]=30.0;
}
for(i=5;i<12;i++)
{
t[i][5]=0.0;
}
for(j=6;j<16;j++)
{
t[5][j]=0.0;
}
for(i=1;i<5;i++)
for(j=1;j<15;j++)
t[i][j]=20,a[i][j]=0;
}
{
for(i=5;i<11;i++)
for(j=1;j<5;j++)
t[i][j]=20,a[i][j]=0;
}
while(dt>=ep)
{
{
for(i=1;i<5;i++)
for(j=1;j<15;j++)
a[i][j]=t[i][j];
}
{
for(i=5;i<11;i++)
for(j=1;j<5;j++)
a[i][j]=t[i][j];
}
{
for(i=5;i<11;i++)
for(j=1;j<5;j++)
t[i][j]=(t[i-1][j]+t[i+1][j]+t[i][j-1]+t[i][j+1])/4;
}
{
for(i=1;i<5;i++)
for(j=1;j<15;j++)
t[i][j]=(t[i-1][j]+t[i+1][j]+t[i][j-1]+t[i][j+1])/4;
}
{
for(j=1;j<5;j++)
t[11][j]=(t[11][j-1]+t[11][j+1]+2*t[10][j])/4;
}
{
for(i=1;i<5;i++)
t[i][15]=(t[i-1][15]+t[i+1][15]+2*t[i][14])/4;
}
{
for(i=1;i<5;i++)
{
for(j=1;j<15;j++)
dt=dt+abs(t[i][j]-a[i][j]);
dt=dt/(4*14);
}
for(i=5;i<11;i++)
{
for(j=1;j<5;j++)
dt=dt+abs(t[i][j]-a[i][j]);
dt=dt/(4*6);
}
}
}
printf("温度分布为:
\t\t\t\t\t\t\t\t\t");
l=0;
{
for(i=0;i<6;i++)
{
for(j=0;j<16;j++)
printf("%4.1f",t[i][j]);
l=l+1;
if(l==16)
{
printf("\n");
l=0;
}
}
l=0;
for(i=6;i<12;i++)
{
for(j=0;j<6;j++)
{printf("%4.1f",t[i][j]);
l=l+1;
if(l==6)
{
printf("\n");
l=0;
}
}
}
}
{
for(j=1;j<15;j++)
q1=q1+(30-t[1][j])*lmd;
for(i=1;i<11;i++)
q1=q1+(30-t[i][1])*lmd;
for(j=6;j<15;j++)
q2=q2+t[4][j]*lmd;
for(i=6;i<11;i++)
q2=q2+t[i][4]*lmd;
q1=q1+lmd*((30-t[1][15])/2+(30-t[11][1])/2);
q2=q2+lmd*(t[4][15]/2+t[11][4]/2+t[4][4]);
}
q=(q1+q2)/2;
e=abs((q2-q1)/q);
printf("单位长度上墙体导热量为:
%4.2fW,偏差为:
%3.2f",4*q,e);
getchar();getchar();
return0;
}
运行结果图:
2、数值模拟程序(matlab):
z=[30.030.030.030.030.030.030.030.030.030.030.030.030.030.030.030.0;
30.029.028.127.126.225.524.924.524.324.224.124.124.024.024.024.0;
30.028.126.124.222.320.719.618.918.518.318.218.118.018.018.018.0;
30.027.124.221.218.115.513.913.012.512.312.212.112.012.012.012.0;
30.026.222.318.113.69.17.46.76.46.26.16.16.06.06.06.0;
30.025.520.715.59.10.00.00.00.00.00.00.00.00.00.00.0;
30.024.919.613.97.40.00.00.00.00.00.00.00.00.00.00.0;
30.024.518.913.06.70.00.00.00.00.00.00.00.00.00.00.0;
30.024.318.512.66.40.00.00.00.00.00.00.00.00.00.00.0;
30.024.218.312.36.20.00.00.00.00.00.00.00.00.00.00.0;
30.024.118.212.26.10.00.00.00.00.00.00.00.00.00.00.0;
30.024.118.212.26.10.00.00.00.00.00.00.00.00.00.00.0];
v=[121824];
x=1:
1:
12;
y=1:
1:
16;
[xx,yy]=meshgrid(y,x);
surf(xx,yy,z);colorbar;xlabel('x');ylabel('y');zlabel('z');az=0;el=-90;view(az,el);
shadinginterp;
axistight;
figure,contour(xx,yy,z,v);
gridon
六、试验感想
这次花了很长时间来完成传热学数值模拟试验,总来说还是收获很大。
首先从试验本身来看,结合书本上相关数值求解讲解,以及在网上搜索部分例题,我对数值计算中起决定性作用划分网格、求各节点对应温度值计算方法有了比较清楚认识。
数值求解在处理实际问题中确实能够发挥很大作用,能最靠近真实情况将实际问题解用数学模型表示出来,从而加深对这些问题了解,同时也能加强对二维稳态导热问题等相关知识掌握。
其次从具体实施角度来看,我在编改程序过程中也再次具体认识和掌握了很多C语言以及matlab相关知识,对迭代求解思绪有了比较清楚地认识。
而且也认识到编写程序所需要严谨与认真是又快又好处理问题关键。
不然浪费了时间精力还做了无用功,得不偿失。
七、参考文件
[1]、西安交通大学等编,传热学试验指导书,热与流体试验中心
[2]、周振红等主编,Fortran90/95高级程序设计,黄河水利出版社,
[3]、杨世铭,陶文铨编著,传热学,高等教育出版社,