傳熱學(xué)編程(顯式和隱式格式)

PAGEPAGE9傳熱學(xué)編程熱動102班學(xué)號：2010031294姓名：唐旭第一部分隱式格式求解C語言程序求解#include<stdio.h>#defineN11#defineK121floatabsf(floatx){if(x<0)x=(-1)*x;returnx;}main(){ inti,j,l; intn;/*迭代次數(shù)*/ 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);/*當(dāng)計算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);/*邊界點溫度用熱平衡法推導(dǎo)出公式*/if(absf(t[i][j+1]-tmp)>max)max=absf(t[i][j+1]-tmp); } } if(max<=eeee)break;} /*輸出溫度分布，其中l(wèi)控制輸出值的排列;這個結(jié)果是橫軸為x，縱軸為τ的直角坐標(biāo)下從左上角開始依次的*/ printf("\n經(jīng)數(shù)值離散計算的溫度分布為(隱式差分格式)：\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; } } }輸出結(jié)果利用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）平板中溫度的瞬態(tài)分布圖線第二部分顯式格式求解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; /*循環(huán)開始，每次計算一個時刻*/ 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];/*當(dāng)計算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;/*邊界點溫度用熱平衡法推導(dǎo)出公式*/ } } } /*輸出溫度分布，其中l(wèi)控制輸出值的排列;這個結(jié)果是橫軸為x，縱軸為τ的直角坐標(biāo)下從左上角開始依次的*/ printf("\n經(jīng)數(shù)值離散計算的溫度分布為（顯式差分格式）：\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; } } }輸出結(jié)果利用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,