传热学编程(显式和隐式格式)

PAGEPAGE9传热学编程热动102班学号：2010031294姓名：唐旭第一部分隐式格式求解C语言程序求解#include<stdio.h>#defineN11#defineK121floatabsf(floatx){if(x<0)x=(-1)*x;returnx;}main(){ inti,j,l; intn;/*迭代次数*/ floata,x,y,Fo,Bi,tmp,max,eeee=0.0000001; floatt[N][K]; y=15;/*y代表Δτ取时间步长为15秒*/ x=0.1/(N-1); a=1.39/100000; Fo=(a*y)/(x*x); Bi=1163*x/50; /*赋予初场温度为80摄氏度*/ for(i=0;i<N;i++) for(j=0;j<K;j++) t[i][j]=80;/*迭代计算*/for(n=0;;n++){max=0; for(j=0;j<K-1;j++) { for(i=0;i<N-1;i++) {tmp=t[i][j+1]; if(i==0) t[i][j+1]=(Fo*(t[i+1][j+1]+t[i+1][j+1])+t[i][j])/(1+2*Fo);/*当计算t[0]时，要用到t[-1],其中t[-1]=t[1]的(对称分布)*/ else t[i][j+1]=(Fo*(t[i-1][j+1]+t[i+1][j+1])+t[i][j])/(1+2*Fo); t[N-1][j+1]=(t[N-1][j]+2*Fo*(t[N-2][j+1]+Bi*300))/(1+2*Fo+2*Fo*Bi);/*边界点温度用热平衡法推导出公式*/if(absf(t[i][j+1]-tmp)>max)max=absf(t[i][j+1]-tmp); } } if(max<=eeee)break;} /*输出温度分布，其中l控制输出值的排列;这个结果是横轴为x，纵轴为τ的直角坐标下从左上角开始依次的*/ printf("\n经数值离散计算的温度分布为(隐式差分格式)：\n"); l=0; for(j=K-1;j>=0;j-=20) for(i=N-1;i>=0;i--) { printf("%6.2f",t[i][j]); l=l+1; if(l==N) { printf("\n"); l=0; } } }输出结果利用MATLAB软件绘图（1）MATLAB绘图语言x=0.0:0.01:0.1;y30=[294.91293.76292.69291.71290.83290.07289.43288.93288.57288.35

288.27];y25=[291.46289.53287.73286.08284.61283.33282.27281.42280.81280.44

280.32];y20=[285.67282.43279.41276.64274.18272.03270.24268.82267.80267.18

266.97]y15=[275.95270.51265.44260.80256.66253.06250.05247.68245.95244.91

244.56]y10=[259.64250.51242.00234.22227.26221.23216.18212.19209.30207.55

206.97]y05=[231.95216.57202.28189.27177.70167.69159.36152.79148.06145.20

144.24]y=[y30;y25;y20;y15;y10;y05];plot(x,y)>>grid>>text(0.1,288.27,'30min')>>text(0.1,280.32,'25min')>>text(0.1,266.97,'20min')>>text(0.1,244.56,'15min')>>text(0.1,206.97,'10min')>>text(0.1,144.24,'5min')（2）平板中温度的瞬态分布图线第二部分显式格式求解C语言程序求解#include<stdio.h>#defineN11#defineK6001main(){ inti,j,l; floata,x,y,Fo,Bi; floatt[N][K]; y=0.3;/*y代表Δτ*/ x=0.1/(N-1); a=1.39/100000; Fo=(a*y)/(x*x); Bi=1163*x/50; printf("\n显式格式条件:"); printf("\n1、Fo=%3.3f<0.5\t",Fo); printf("\t2、1-2Fo*Bi-2Fo=%4.6f>0\n\n",1-2*Fo*Bi-2*Fo); /*时刻为零时，赋予初场温度*/ for(i=0;i<N;i++) t[i][0]=80; /*循环开始，每次计算一个时刻*/ for(j=0;j<K-1;j++) { for(i=0;i<N;i++) { 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[1]的(对称分布)*/ 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-1][j]*(1-2*Fo*Bi-2*Fo)+2*Fo*t[N-2][j]+2*Fo*Bi*300;/*边界点温度用热平衡法推导出公式*/ } } } /*输出温度分布，其中l控制输出值的排列;这个结果是横轴为x，纵轴为τ的直角坐标下从左上角开始依次的*/ printf("\n经数值离散计算的温度分布为（显式差分格式）：\n"); l=0; for(j=K-1;j>=0;j-=1000) for(i=N-1;i>=0;i--) { printf("%6.2f",t[i][j]); l=l+1; if(l==N) { printf("\n"); l=0; } } }输出结果利用MATLAB软件绘图（1）MATLAB绘图语言>>x=0.0:0.01:0.1;y30=[295.12294.01292.98292.04291.20290.47289.86289.38289.03288.82288.75]y25=[291.75289.88288.15286.56285.13283.90282.87282.05281.46281.10280.98]y20=[286.06282.91279.97277.28274.88272.79271.05269.67268.67268.07267.87]y15=[276.45271.12266.15261.61257.55254.02251.08248.75247.06246.04245.70]y10=[260.20251.19242.79235.12228.26222.31217.33213.39210.54208.82208.24]y05=[232.58217.33203.14190.20178.65168.65160.31153.72148.96146.09145.12]y=[y30;y25;y20;y15;y10;y05];>>grid>>text(0.1,288.75,'30min')>>text(0.1,280.98,'25min')>>text(0.1,276.87,'20min')>>plot(x,