一维非稳态导热的数值计算

传热学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,,-1
2
2
2+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文件结果算的后不会立即退出，方便观看*/
}。