传热学上机C程序源答案之一维非稳态导热的数值计算

二维稳态导热的数值计算
2.1物理问题
一矩形区域，其边长L=W=1，假设区域内无内热源，导热系数为常数，三个边温度为T1=0，一个边温度为T2=1，求该矩形区域内的温度分布。

2.2 数学描述 对上述问题的微分方程及其边界条件为：2222T T 0x y
∂∂+=∂∂ x=0，T=T 1=0
x=1，T=T 1=0
y=0，T=T 1=0
y=1，T=T 2=1 该问题的解析解：112121(1)sin n n n sh y T T n L x n T T n L sh W L ππππ∞=⎛⎫⋅ ⎪---⎛⎫⎝⎭=⋅ ⎪-⎛⎫⎝⎭⋅ ⎪⎝⎭
∑ 2.3数值离散
2.3.1区域离散
区域离散x 方向总节点数为N ，y 方向总节点数为M ，区域内任一节点用I,j 表示。

2.3.2方程的离散 对于图中所有的内部节点方程可写为：2222,,0i j i j
t t x y ⎛⎫⎛⎫∂∂+= ⎪ ⎪∂∂⎝⎭⎝⎭ 用I,j 节点的二阶中心差分代替上式中的二阶导数，得：
+1,,-1,,+1,,-1222+2+0i j i j i j
i j i j i j T T T T T T x y --+=
上式整理成迭代形式：()()22
,1,-1,,1,-12222+2()
2()i j i j i j i j i j y x T T T T T x y x y ++=++++ (i=2,3……,N -1)，(j=2,3……,M -1)
补充四个边界上的第一类边界条件得：1,1j T T = (j=1,2,3……,M)
,1N j T T = (j=1,2,3……,M)
,1i j T T = (i=1,2,3……,N)
,2i M T T = (i=1,2,3……,N)
#include<stdio.h>
#include<math.h>
#define N 10
#define K 11
main()
{
int i,j,l;
float cha;
float a,x,y,Fo,Bi;
float t[N][K],b[N][K];
/*打印出题目*/
printf("\t\t\t一维非稳态导热问题\t\t");
printf("\n\t\t\t\t\t\t----何鹏举\n");
printf("\n题目：补充材料练习题三\n");
y=1;/*y代表Δτ*/
x=0.05/(N-1);
a=34.89/(7800*712);
Fo=(a*y)/(x*x);
Bi=233*x/34.89;
printf("\n显示格式条件:");
printf("\n1、Fo=%3.1f<0.5\t",Fo);
printf("\t2、1-2Fo*Bi-2Fo=%4.2f>0\n\n",1-2*Fo*Bi-2*Fo);
/*时刻为零时，赋予初场温度*/
for(i=0;i<N;i++)
t[i][0]=1000;
/*循环开始，每次计算一个时刻*/
for(j=0;j<K-1;j++)
{
for(i=0;i<N;i++)
b[i][j]=t[i][j];
/*下面对每一个时刻进行迭代求解对应的温度分布，公式按传热学课本P178页公式*/ cha=1;
while(cha>0.001)
{
for(i=0;i<N-1;i++)
{
if(i==0)
t[i][j+1]=Fo*(t[i+1][j]+t[i+1][j])+(1-2*Fo)*t[i][j];/*当计算t[0]时，要用到t[-1],其中t[-1]=t[2]的(对称分布)*/
else
t[i][j+1]=Fo*(t[i+1][j]+t[i-1][j])+(1-2*Fo)*t[i][j];
t[N-1][j+1]=t[N-2][j]*(1-2*Fo*Bi-2*Fo)+2*Fo*t[N-1][j]+2*Fo*Bi*20;/*边界点温度用热平衡法推导出公式*/
}
cha=0;
for(i=0;i<N;i++)
cha=cha+abs(t[i][j]-b[i][j]);
cha=cha/N;
}
}
/*输出温度分布，其中l控制输出值的排列;这个结果是横轴为x，纵轴为τ的直角坐标下从左上角开始依次的*/
printf("\n经数值离散计算的温度分布为：\n");
l=0;
for(j=K-1;j>=0;j--)
for(i=0;i<N;i++)
{
if(t[i][j]>999.99)
printf("%6.1f ",t[i][j]);
else
printf("%6.2f ",t[i][j]);
l=l+1;
if(l==N)
{
printf("\n");
l=0;
}
}
getchar();/*为了是生成的exe文件结果算的后不会立即退出，方便观看*/
}。