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

传热学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 #include #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

/*循环开始，每次计算一个时刻*/

for(j=0;j

for(i=0;i

/*下面对每一个时刻进行迭代求解对应的温度分布，公式按传热学课本P178页公式*/ cha=1;

while(cha>0.001) {

for(i=0;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

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

{

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文件结果算的后不会立即退出，方便观看*/

}