清華大學(xué)傳熱大作業(yè)-無限大平壁簡單數(shù)值計算報告(共6頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上傳熱大作業(yè)報告* 熱動* *一、 大作業(yè)題目一厚度為0.1m的無限大平壁，兩側(cè)均為對流換熱邊界條件，初始時 兩側(cè)流體溫度與壁內(nèi)溫度一致，tf1=tf2=t0=5 ；已知兩側(cè)對流換熱系數(shù)分別為h1=11 W/m2K、h2=23W/m2K, 平壁材料的導(dǎo)熱系數(shù)l=0.43W/mK，導(dǎo)溫系數(shù)a=0.3437×10-6 m2/s。如果一側(cè)的環(huán)境溫度tf1突然升高為50并維持不變，計算在其它參數(shù)不變的條件下，平壁內(nèi)溫度分布及兩側(cè)壁面熱流密度隨時間的變化規(guī)律（用圖形表示）。要求：將全部計算內(nèi)容（包括網(wǎng)格的劃分、節(jié)點方程組、計算框圖、程序及計算結(jié)果）用A4紙打印。二、 網(wǎng)

2、格劃分如圖，將無限大平板作為一維處理，本題為一維非穩(wěn)態(tài)導(dǎo)熱問題，對流換熱邊界條件。l 空間網(wǎng)格劃分：平板總厚度為delta=0.1m，定義空間步長為 dx=0.005m，則距離份數(shù)為N=delta/dx=20份。定義xn為以0為首項，以dx為公差的等差數(shù)列，尾項為delta=0.1m，共有N+1項，則xn中的每一項即表示一個沿平板厚度方向中的劃分點。l 時間網(wǎng)格劃分：設(shè)總時間長度為T= s，定義時間步長為dtao=20s，則時間份數(shù)為M=5000份。定義 taom是以0為首項，以dtao為公差的等差數(shù)列，尾項為T=s，共有M+1項,則taom中每一項即表示一個時刻。三、 計算框圖l 程序中的各

3、個變量的名稱及意義：1. 題設(shè)中各個常數(shù)lambda=0.43 導(dǎo)熱系數(shù)；a=0.3437e-6 熱擴散率；h1=11 邊界對流換熱系數(shù)；h2=23 邊界對流換熱系數(shù)2；t0=5 初始溫度；tf1=50 初始流體溫度；tf2=5 初始流體溫度2；delta=0.1 總距離長度（無限大平板厚度）；2. 網(wǎng)格劃分所設(shè)的變量T= 總時間長度（在T時間內(nèi)考慮本問題）；dtao=20 定義時間步長；dx=0.005 定義距離步長；M=floor(T/dtao) 時間份數(shù)=總時間/時間步長（向下取整）；N=floor(delta/dx 距離份數(shù)=總厚度/距離步長(向下取整）；tao=0:dtao:T 定義

4、時間劃分單元（以0為首項，以dtao為公差的等差數(shù)列，尾項為T），共有M+1項；x=0:dx:delta 定義距離劃分單元（以0為首項，以dx為公差的等差數(shù)列，尾項為delta），共有N+1項；3. 判定穩(wěn)定性的準(zhǔn)則數(shù)Bi1=h1*dx/lambda 邊界節(jié)點網(wǎng)格畢渥數(shù)；Bi2=h2*dx/lambda 邊界節(jié)點網(wǎng)格畢渥數(shù)2；Fo=a*dtao/dx2 傅里葉數(shù)；l 程序計算框圖開始輸入delta,T,dtao,dx,M,N,tao,xFo>1/(2*Bi1+2)&&Fo>1/(2*Bi2+2)NoYes建立t(M+1,N+1)溫度矩陣，令t(1,:)=t0,令m=

5、1m=m+1打印“不穩(wěn)定”t(m,1)=2*Fo*(t(m-1,2)+Bi1*tf1)+(1-2*Bi1*Fo-2*Fo)*t(m-1,1)t(m,N+1)=2*Fo*(t(m-1,N)+Bi2*tf2)+(1-2*Bi2*Fo-2*Fo)*t(m-1,N+1)t(m,n)=Fo*(t(m-1,n-1)+t(m-1,n+1)+(1-2*Fo)*t(m-1,n)m>M+1NoYes輸出溫度矩陣t(M+1,N+1)和相應(yīng)圖象停機四、 程序代碼本程序在MATLAB R2008a中運行通過，以下是源代碼（%后為注釋）：lambda=0.43;%導(dǎo)熱系數(shù)a=0.3437e-6;%熱擴散率h1=11

6、;%邊界對流換熱系數(shù)h2=23;%邊界對流換熱系數(shù)2t0=5;%初始溫度tf1=50;%初始流體溫度tf2=5;%初始流體溫度2delta=0.1;%總距離長度（無限大平板厚度）T=;%總時間長度（在T時間內(nèi)考慮本問題）dtao=20;%定義時間步長dx=0.005;%定義距離步長M=floor(T/dtao);%時間份數(shù)=總時間/時間步長（向下取整）N=floor(delta/dx);%距離份數(shù)=總厚度/距離步長(向下取整）tao=0:dtao:T;%定義時間劃分單元（以0為首項，以dtao為公差的等差數(shù)列，尾項為T），共有M+1項x=0:dx:delta;%定義距離劃分單元（以0為首項，以

7、dx為公差的等差數(shù)列，尾項為delta），共有N+1項Bi1=h1*dx/lambda;%邊界節(jié)點網(wǎng)格畢渥數(shù)Bi2=h2*dx/lambda;%邊界節(jié)點網(wǎng)格畢渥數(shù)2Fo=a*dtao/dx2;%傅里葉數(shù)if Fo>1/(2*Bi1+2)&&Fo>1/(2*Bi2+2)%判斷穩(wěn)定性，不穩(wěn)定則顯示畢渥數(shù)、傅里葉數(shù) disp('不穩(wěn)定'); disp(Bi1); disp(Bi2); disp(Fo); disp(1/(2*Bi1+2); disp(1/(2*Bi2+2);else%若穩(wěn)定，則進行迭代計算 t=zeros(M+1,N+1);%建立一個(M+

8、1)*(N+1)的溫度矩陣，M+1為時間節(jié)點個數(shù)，N+1為空間節(jié)點個數(shù)，以便進行迭代計算 q1=zeros(M+1,1);%根據(jù)題目要求算兩壁面處熱流密度 q2=zeros(M+1,1); t(1,:)=t0;%初始溫度均為t0=5 for m=2:M+1%m=1時是初值上一行已計算出，則從m=2一直計算到m=M+1，m對應(yīng)的時刻是tao=（m-1)dtao t(m,1)=2*Fo*(t(m-1,2)+Bi1*tf1)+(1-2*Bi1*Fo-2*Fo)*t(m-1,1);%首先計算一邊界這個時刻溫度 t(m,N+1)=2*Fo*(t(m-1,N)+Bi2*tf2)+(1-2*Bi2*Fo-2

9、*Fo)*t(m-1,N+1);%再計算另一邊界這個時刻的溫度 q1(m)=h1*(tf1-t(m,1); q2(m)=h2*(t(m,N+1)-tf2); for n=2:N%然后計算內(nèi)部，n=1和n=N+1時是邊界節(jié)點溫度，上面兩行已經(jīng)計算出，n對應(yīng)的坐標(biāo)是x=(n-1)*dx t(m,n)=Fo*(t(m-1,n-1)+t(m-1,n+1)+(1-2*Fo)*t(m-1,n); end end %以下是畫圖 figure plot(x,t(1,:),x,t(11,:),x,t(21,:),x,t(51,:),x,t(101,:),x,t(1001,:),x,t(5001,:); lege

10、nd('t=0s','t=200s','t=400s','t=1000s','t=2000s','t=20000s','t=',0); title('一定時間下溫度隨距離的分布','fontsize',12,'fontweight','bold','fontname','楷體'); axis(0,0.1,0,40); figure plot(tao,t(:,1),tao,t(:,6),ta

11、o,t(:,11),tao,t(:,16),tao,t(:,21); legend('x=0','x=0.025','x=0.05','x=0.075','0.1',0); title('一定位置處溫度隨時間的分布','fontsize',12,'fontweight','bold','fontname','楷體'); axis(0,0,40); figure mesh(x,tao,t); title('溫度隨時間和空間的分布','fontsize',12,'fontweight','bold','fontname','楷體'); figure plot(tao,q1,tao,q2); legend('q1','q2'); title('兩壁面熱流密度隨時間變化曲線','fontsize',12,'fontweight','bold','fon