一維非穩(wěn)態(tài)導(dǎo)熱問(wèn)題的數(shù)值解

1、 計(jì)算傳熱學(xué)程序報(bào)告 題目：一維非穩(wěn)態(tài)導(dǎo)熱問(wèn)題的數(shù)值解 姓名： 學(xué)號(hào)： 學(xué)院：能源與動(dòng)力工程學(xué)院 專(zhuān)業(yè)：工程熱物理 日期：2014年5月25日一維非穩(wěn)態(tài)導(dǎo)熱問(wèn)題數(shù)值解求解下列熱傳導(dǎo)問(wèn)題：1.方程離散化對(duì)方程進(jìn)行控制體積分得到： 非穩(wěn)態(tài)項(xiàng)：選取T隨x階梯式變化，有 擴(kuò)散項(xiàng)：選取一階導(dǎo)數(shù)隨時(shí)間做顯示變化，有 進(jìn)一步取T隨x呈分段線性變化，有 ， 整理可以得到總的離散方程為： 2.計(jì)算空間和時(shí)間步長(zhǎng)取空間步長(zhǎng)為： h=L/N網(wǎng)格Fourier數(shù)為： （小于0.5時(shí)穩(wěn)定）時(shí)間步長(zhǎng)為： 3. 建立溫度矩陣與邊界條件 T=ones(N+1,M+1) T(:,1)=Ti (初始條件溫度都為0） T(1,:)

2、=To (邊界條件x=0處溫度為1） T(N+1,:)=Te (邊界條件x=L處溫度為0）4. 差分法求解溫度由離散方程可得到： 轉(zhuǎn)化為相應(yīng)的溫度矩陣形式： 5. 輸入界面考慮到方程的變量，采用inputdlg函數(shù)設(shè)置5個(gè)輸入變量，對(duì)這5個(gè)變量設(shè)置了默認(rèn)值，如圖1所示。在計(jì)算中可以改變不同的數(shù)值，得到不同的結(jié)果，特別注意穩(wěn)定條件的臨界值是0.5。根據(jù)設(shè)置的默認(rèn)值，得到的計(jì)算結(jié)果如圖2所示。 圖1 matlab變量輸入界面 圖2 默認(rèn)值的計(jì)算結(jié)果6. 結(jié)果分析 根據(jù)上面的分析，給出了程序的輸入界面，以及默認(rèn)值狀態(tài)下的數(shù)值解?？梢酝ㄟ^(guò)改變不同的輸入值，得到需要的分析結(jié)果，總結(jié)出了下面4點(diǎn)結(jié)論： （

3、1） 取F0=0.48，得到一維非穩(wěn)態(tài)導(dǎo)熱結(jié)果如下圖所示 圖2 F0=0.48時(shí)一維非穩(wěn)態(tài)導(dǎo)熱從圖中可以看出，對(duì)于長(zhǎng)度L=1的細(xì)桿，初始時(shí)刻t=0時(shí)溫度為0，邊界條件x=0時(shí)，T=1,邊界條件x=1時(shí)，T=0。隨著時(shí)間的增加，溫度從x=0通過(guò)導(dǎo)熱的形式傳遞到x=1，不同時(shí)刻不同位置桿的溫度都不同，并且隨著時(shí)間的增加，桿的溫度也逐漸增加。（2） 取F0=0.48，可以得到不同位置的溫度響應(yīng)曲線，如下圖所示 圖3 F0=0.48時(shí)不同x位置處的溫度響應(yīng)圖中紅色曲線代表x=0.1位置的溫度瞬態(tài)響應(yīng)，黑色曲線代表x=0.2位置的溫度瞬態(tài)響應(yīng)，藍(lán)色曲線代表x=0.4位置的溫度瞬態(tài)響應(yīng)。從圖中可以看出，隨

4、著x的增加，曲線與x軸的交點(diǎn)值越大，溫度開(kāi)始傳遞到該位置的所需的時(shí)間越長(zhǎng)。隨著x的增加，溫度響應(yīng)曲線的變化速率越慢，最終的達(dá)到的溫度也越低。（3） 取F0=0.25，得到不同位置的溫度響應(yīng)曲線如下圖所示 圖4 F0=0.25時(shí)不同x位置處的溫度響應(yīng)圖中三條曲線分別是x=0.1，x=0.2，x=0.4位置的溫度瞬態(tài)響應(yīng)。與圖3的F0=0.48進(jìn)行對(duì)比，兩種情況下的F0值不同，F(xiàn)0值越大表明熱擴(kuò)散系數(shù)的值越大。從圖中可以看出熱擴(kuò)散系數(shù)對(duì)于導(dǎo)熱的影響，F(xiàn)0=0.25時(shí)，與F0=0.48相比較，各位置開(kāi)始響應(yīng)時(shí)所需的時(shí)間較長(zhǎng)，而且各位置響應(yīng)曲線的變化速率較小，最終的達(dá)到的溫度也較低，說(shuō)明了熱擴(kuò)散系數(shù)越

5、小，熱傳導(dǎo)越慢，傳遞效率越低。（4） 取F0=0.51，得到非穩(wěn)定的數(shù)值解如圖所示 圖5 F0=0.51時(shí)一維非穩(wěn)態(tài)導(dǎo)熱 圖6 F0=0.51時(shí)不同x位置處的溫度響應(yīng)從圖中可以看出，對(duì)于顯示格式的離散方程，并不是所有的F0值都能得到有意義的解，必須要求F0<0.5時(shí)才能得到穩(wěn)定的數(shù)值解，當(dāng)F0>0.5時(shí)，會(huì)出現(xiàn)物理上不真實(shí)的解。 附件：（matlab程序）function heat_conduction() %一維齊次熱傳導(dǎo)方程%設(shè)置輸入界面options='空間桿長(zhǎng)L','空間點(diǎn)數(shù)N' ,'時(shí)間點(diǎn)數(shù)M','擴(kuò)散系數(shù)a'

6、;,'穩(wěn)定條件的值Fo(臨界值0.5)',;topic='一維非穩(wěn)態(tài)導(dǎo)熱'%標(biāo)題欄顯示lines=1;%輸入行為1行def='1','100','1000','1','0.48'%默認(rèn)值輸入f=inputdlg(options,topic,lines,def);%輸入框設(shè)置L=eval(f1);%設(shè)置輸入值N=eval(f2);M=eval(f3);a=eval(f4);Fo=eval(f5);%Fo的值必須小于0.5，小于0.5波動(dòng)%計(jì)算空間步長(zhǎng)與時(shí)間步長(zhǎng)h=L/N;%空間步長(zhǎng)x1

7、=0:h:L;x=x1'n=Fo*h2/a;%時(shí)間步長(zhǎng)tm=n*M;%傳導(dǎo)總時(shí)間t1=0:n:tm;t=t1'%計(jì)算初始條件與邊界條件Ti=x.*0;%初始條件To=1+t.*0;%x=0的邊界條件Te=t.*0;%x=L的邊界條件%建立溫度矩陣TT=ones(N+1,M+1);T(:,1)=Ti;%第一列為初始條件T(1,:)=To;%第一行為x=0邊界條件T(N+1,:)=Te;%最后一行為x=L邊界條件%利用差分法求解溫度矩陣Tfor k=1:M m=2; while m<=N; T(m,k+1)=Fo*(T(m+1,k)+T(m-1,k)-2*T(m,k)+T(m

8、,k); m=m+1; endend%將時(shí)間空間的一維坐標(biāo)轉(zhuǎn)化為二維坐標(biāo)Y,X=meshgrid(t1,x);%根據(jù)溫度矩陣T繪圖subplot(2,2,1);mesh(X,Y,T);%三維圖繪制view(1,-1,1);%調(diào)整視圖角度title('非穩(wěn)態(tài)導(dǎo)熱');%圖像名稱(chēng)xlabel('長(zhǎng)度x');%x軸名稱(chēng)ylabel('時(shí)間t');%y軸名稱(chēng)zlabel('溫度T');%z軸名稱(chēng)subplot(2,2,2);A=T(11,:);%取矩陣第11列的值plot(A,'r');%二維曲線繪制legend('

9、;A=0.1');%顯示函數(shù)名稱(chēng)title('x=0.1瞬態(tài)響應(yīng)');xlabel('時(shí)間t');ylabel('溫度T');axis(0 1000 0 1);%坐標(biāo)軸數(shù)值范圍subplot(2,2,3);B=T(21,:);%取矩陣第21列plot(B,'k');legend('B=0.2');title('x=0.2瞬態(tài)響應(yīng)');xlabel('時(shí)間t');ylabel('溫度T');axis(0 1000 0 1);subplot(2,2,4);C=T(41,:);%取矩陣第41列plot(A,'r');hold on;%多條曲線繪制plot(B,'k');plot(C);ho