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

传热学C程序源二维稳态导热的数值计算2.1物理问题一矩形区域，其边长L=W=1，假设区域内无内热源，导热系数为常数，三个边温度为T1=0，一个边温度为T2=1，求该矩形区域内的温度分布。2.2数学描述对上述问题的微分方程及其边界条件为：2222TT0xyx=0，T=T1=0x=1，T=T1=0y=0，T=T1=0y=1，T=T2=1该问题的解析解：112121(1)sinnnnshyTTnLxnTTnLshWL2.3数值离散2.3.1区域离散区域离散x方向总节点数为N，y方向总节点数为M，区域内任一节点用I,j表示。2.3.2方程的离散对于图中所有的内部节点方程可写为：2222,,0ijijttxy用I,j节点的二阶中心差分代替上式中的二阶导数，得：+1,,-1,,+1,,-1222+2+0ijijijijijijTTTTTTxy上式整理成迭代形式：22,1,-1,,1,-12222+2()2()ijijijijijyxTTTTTxyxy(i=2,3……,N-1)，(j=2,3……,M-1)补充四个边界上的第一类边界条件得：1,1jTT(j=1,2,3……,M),1NjTT(j=1,2,3……,M),1ijTT(i=1,2,3……,N),2iMTT(i=1,2,3……,N)传热学C程序源之一维非稳态导热的数值计算#includestdio.h#includemath.h#defineN10#defineK11main(){inti,j,l;floatcha;floata,x,y,Fo,Bi;floatt[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.1f0.5\t,Fo);printf(\t2、1-2Fo*Bi-2Fo=%4.2f0\n\n,1-2*Fo*Bi-2*Fo);/*时刻为零时，赋予初场温度*/for(i=0;iN;i++)t[i][0]=1000;/*循环开始，每次计算一个时刻*/for(j=0;jK-1;j++){for(i=0;iN;i++)b[i][j]=t[i][j];/*下面对每一个时刻进行迭代求解对应的温度分布，公式按传热学课本P178页公式*/cha=1;while(cha0.001){for(i=0;iN-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]的(对称分布)*/elset[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;iN;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;iN;i++){if(t[i][j]999.99)printf(%6.1f,t[i][j]);elseprintf(%6.2f,t[i][j]);l=l+1;if(l==N){printf(\n);l=0;}}getchar();/*为了是生成的exe文件结果算的后不会立即退出，方便观看*/}