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

一維非穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維非穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解一維非穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解V:1.0精細整理，僅供參考一維非穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解日期：20xx年X月計算傳熱學(xué)程序報告題目：一維非穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解姓名：學(xué)號：學(xué)院：能源與動力工程學(xué)院專業(yè)：工程熱物理日期：2014年5月25日一維非穩(wěn)態(tài)導(dǎo)熱問題數(shù)值解求解下列熱傳導(dǎo)問題：1.方程離散化對方程進行控制體積分得到：非穩(wěn)態(tài)項：選取T隨x階梯式變化，有擴散項：選取一階導(dǎo)數(shù)隨時間做顯示變化，有進一步取T隨x呈分段線性變化，有，整理可以得到總的離散方程為：2.計算空間和時間步長取空間步長為：h=L/N網(wǎng)格Fourier數(shù)為：（小于時穩(wěn)定）時間步長為：建立溫度矩陣與邊界條件T=ones(N+1,M+1)T(:,1)=Ti(初始條件溫度都為0）T(1,:)=To(邊界條件x=0處溫度為1）T(N+1,:)=Te(邊界條件x=L處溫度為0）差分法求解溫度由離散方程可得到：轉(zhuǎn)化為相應(yīng)的溫度矩陣形式：輸入界面考慮到方程的變量，采用inputdlg函數(shù)設(shè)置5個輸入變量，對這5個變量設(shè)置了默認值，如圖1所示。在計算中可以改變不同的數(shù)值，得到不同的結(jié)果，特別注意穩(wěn)定條件的臨界值是。根據(jù)設(shè)置的默認值，得到的計算結(jié)果如圖2所示。圖1matlab變量輸入界面圖2默認值的計算結(jié)果結(jié)果分析根據(jù)上面的分析，給出了程序的輸入界面，以及默認值狀態(tài)下的數(shù)值解。可以通過改變不同的輸入值，得到需要的分析結(jié)果，總結(jié)出了下面4點結(jié)論：取F0=，得到一維非穩(wěn)態(tài)導(dǎo)熱結(jié)果如下圖所示圖2F0=時一維非穩(wěn)態(tài)導(dǎo)熱從圖中可以看出，對于長度L=1的細桿，初始時刻t=0時溫度為0，邊界條件x=0時，T=1,邊界條件x=1時，T=0。隨著時間的增加，溫度從x=0通過導(dǎo)熱的形式傳遞到x=1，不同時刻不同位置桿的溫度都不同，并且隨著時間的增加，桿的溫度也逐漸增加。取F0=，可以得到不同位置的溫度響應(yīng)曲線，如下圖所示圖3F0=時不同x位置處的溫度響應(yīng)圖中紅色曲線代表x=位置的溫度瞬態(tài)響應(yīng)，黑色曲線代表x=位置的溫度瞬態(tài)響應(yīng)，藍色曲線代表x=位置的溫度瞬態(tài)響應(yīng)。從圖中可以看出，隨著x的增加，曲線與x軸的交點值越大，溫度開始傳遞到該位置的所需的時間越長。隨著x的增加，溫度響應(yīng)曲線的變化速率越慢，最終的達到的溫度也越低。取F0=，得到不同位置的溫度響應(yīng)曲線如下圖所示圖4F0=時不同x位置處的溫度響應(yīng)圖中三條曲線分別是x=，x=，x=位置的溫度瞬態(tài)響應(yīng)。與圖3的F0=進行對比，兩種情況下的F0值不同，F(xiàn)0值越大表明熱擴散系數(shù)的值越大。從圖中可以看出熱擴散系數(shù)對于導(dǎo)熱的影響，F(xiàn)0=時，與F0=相比較，各位置開始響應(yīng)時所需的時間較長，而且各位置響應(yīng)曲線的變化速率較小，最終的達到的溫度也較低，說明了熱擴散系數(shù)越小，熱傳導(dǎo)越慢，傳遞效率越低。取F0=，得到非穩(wěn)定的數(shù)值解如圖所示圖5F0=時一維非穩(wěn)態(tài)導(dǎo)熱圖6F0=時不同x位置處的溫度響應(yīng)從圖中可以看出，對于顯示格式的離散方程，并不是所有的F0值都能得到有意義的解，必須要求F0<時才能得到穩(wěn)定的數(shù)值解，當(dāng)F0>時，會出現(xiàn)物理上不真實的解。附件：（matlab程序）functionheat_conduction()%一維齊次熱傳導(dǎo)方程%設(shè)置輸入界面options={'空間桿長L','空間點數(shù)N','時間點數(shù)M','擴散系數(shù)a','穩(wěn)定條件的值Fo(臨界值',};topic='一維非穩(wěn)態(tài)導(dǎo)熱';%標(biāo)題欄顯示lines=1;%輸入行為1行def={'1','100','1000','1',''};%默認值輸入f=inputdlg(options,topic,lines,def);%輸入框設(shè)置L=eval(f{1});%設(shè)置輸入值N=eval(f{2});M=eval(f{3});a=eval(f{4});Fo=eval(f{5});%Fo的值必須小于，小于波動%計算空間步長與時間步長h=L/N;%空間步長x1=0:h:L;x=x1';n=Fo*h^2/a;%時間步長tm=n*M;%傳導(dǎo)總時間t1=0:n:tm;t=t1';%計算初始條件與邊界條件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邊界條件%利用差分法求解溫度矩陣Tfork=1:Mm=2;whilem<=N;T(m,k+1)=Fo*(T(m+1,k)+T(m-1,k)-2*T(m,k))+T(m,k);m=m+1;endend%將時間空間的一維坐標(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)熱');%圖像名稱xlabel('長度x');%x軸名稱ylabel('時間t');%y軸名稱zlabel('溫度T');%z軸名稱subplot(2,2,2);A=T(11,:);%取矩陣第11列的值plot(A,'r');%二維曲線繪制legend('A=');%顯示函數(shù)名稱title('x=瞬態(tài)響應(yīng)');xlabel('時間t');ylabel('溫度T');axis([0100001]);%坐標(biāo)軸數(shù)值范圍subplot(2,2,3);B=T(21,:);%取矩陣第21列plot(B,'k');legend('B=');title('x=瞬態(tài)響應(yīng)');xlabel('時間t');ylabel('溫度T');axis([0100001]);subplot