二維導(dǎo)熱物體溫度場的數(shù)值模擬

1、二維導(dǎo)熱物體溫度場的數(shù)值模擬作者： 學(xué) 號： 學(xué)院(系)：能源與動力工程學(xué)院專業(yè)：能源動力系統(tǒng)及自動化班級： 二維導(dǎo)熱物體溫度場的數(shù)值模擬 一：物理問題 有一個用磚砌成的長方形截面的冷空氣通道，其截面尺寸和示意圖如圖1-1所示，假設(shè)在垂直紙面方向上冷空氣及磚墻的溫度變化很小，可以近似地予以忽略。在下列兩種情況下試計算：（1）磚墻橫截面上的溫度分布；（2）垂直于紙面方向的每米長度上通過磚墻的導(dǎo)熱量。 第一種情況：內(nèi)外壁分布均勻地維持在0及30； 第二種情況：內(nèi)外表面均為第三類邊界條件，且已知： 磚墻的導(dǎo)熱系數(shù) 二：數(shù)學(xué)描述 該結(jié)構(gòu)的導(dǎo)熱問題可以作為二維問題處理，并且其截面如圖1-1所示，由于對稱

2、性，僅研究其1/4部分即可。 其網(wǎng)絡(luò)節(jié)點(diǎn)劃分如圖 f a c （m，n） b = n e m d 上述問題為二維矩形域內(nèi)的穩(wěn)態(tài)、無內(nèi)熱源、常物性的導(dǎo)熱問題，對于這樣的物理問題，我們知道，描寫其的微分方程即控制方程，就是導(dǎo)熱微分方程： 第一類邊界條件：內(nèi)外壁分布均勻地維持在0及30；=30=0 第三類邊界條件：內(nèi)外表面均為第三類邊界條件，且已知： 磚墻的導(dǎo)熱系數(shù)三：方程的離散如上圖所示，用一系列與坐標(biāo)軸平行的網(wǎng)絡(luò)線把求解區(qū)域劃分成許多子區(qū)域，以網(wǎng)格線的交點(diǎn)作為需要確定溫度值的空間位置，即節(jié)點(diǎn)，節(jié)點(diǎn)的位置已該點(diǎn)在兩個方向上的標(biāo)號m、n來表示。每一個節(jié)點(diǎn)都可以看成是以它為中心的小區(qū)域的代表，如上（m

3、，n）：對于（m，n）為內(nèi)節(jié)點(diǎn)時：由熱平衡法可以得到，當(dāng)=時： 1 對于（m，n）為邊界節(jié)點(diǎn)時：l 恒溫邊界只需特殊考慮位于絕熱平直邊界上的節(jié)點(diǎn)：l 對流邊界分為角點(diǎn)、絕熱邊界點(diǎn)和對流邊界點(diǎn)。1.絕熱邊界點(diǎn)：2.對流邊界點(diǎn)：3.外角點(diǎn)：4.內(nèi)角點(diǎn)：四：編程思路及流程圖 開始 輸入已知參數(shù)說明邊界條件取定初始試探值ta（i, j）=0計算新的內(nèi)節(jié)點(diǎn)和邊界點(diǎn)溫度t(i, j)t=ta比較所有節(jié)點(diǎn)|ta(i, j)- t(i, j)|<epsilo ？否是計算內(nèi)外邊界散熱量及熱平衡偏差輸出內(nèi)外邊界散熱量及熱平衡偏差結(jié)束畫出溫度場模擬圖 matlab程序如下function=wendu()t=z

4、eros(12,16);tf=zeros(12,16);q1x=0;q1y=0;q1=0;q2x=0;q2y=0;q2=0;n=0;for i=1:12 t(i,1)=30;endfor j=1:16 t(1,j)=30;endt0=t;t=diedai1(t);for i=2:5 for j=2:16 while t(i,j)-t0(i,j)>=0.00001; t0=t; t=diedai1(t); end endendfor i=6:12 for j=2:5 while t(i,j)-t0(i,j)>=0.00001; t0=t; t=diedai1(t); end ende

5、ndt0=t for i=2:11 q1x=q1x+0.53*(t(i,1)-t(i,2);endq1x=q1x+0.53*(t(12,1)-t(12,2)/2;for j=2:15 q1y=q1y+0.53*(t(1,j)-t(2,j);endq1y=q1y+0.53*(t(1,16)-t(2,16)/2;q1=(q1x+q1y)*4for i=6:11 q2x=q2x+0.53*(t(i,5)-t(i,6);endq2x=q2x+0.53*(t(12,5)-t(12,6)/2;for j=6:15 q2y=q2y+0.53*(t(5,j)-t(6,j);endq2y=q2y+0.53*(t

6、(5,16)-t(6,16)/2;q2=(q2x+q2y)*4n=2*abs(q1-q2)/(q1+q2) t0=tf;tf=diedai2(tf);for i=2:5 for j=2:16 while tf(i,j)-t0(i,j)>=0.00001; t0=tf; tf=diedai2(tf); end endendfor i=6:12 for j=2:5 while tf(i,j)-t0(i,j)>=0.00001; t0=tf; tf=diedai2(tf); end endendt0=tf q1x=0;q1y=0;q1=0;q2x=0;q2y=0;q2=0;n=0;for

7、 i=1:11 q1x=q1x+10*0.1*(30-tf(i,1);endq1x=q1x+10*0.05*(30-tf(12,1);for j=2:15 q1y=q1y+10*0.1*(30-tf(1,j);endq1y=q1y+10*0.05*(30-tf(1,16);for i=6:11 q2x=q2x+4*0.1*(tf(i,6)-10);endq2x=q2x+4*0.05*(tf(12,6)-10);for j=7:15 q2y=q2y+4*0.1*(tf(6,j)-10);endq2y=q2y+4*0.05*(tf(6,16)-10);q1=(q1x+q1y)*4;q2=(q2x+

8、q2y)*4;n=2*abs(q1-q2)/(q1+q2) for i=7:12 for j=7:16 tf(i,j)=10; endendsubplot(211);pcolor(t)shading interp;colormap(hot)hold oncontour(t,3,'k')hold offcolorbark=caxis;subplot(212);pcolor(tf)shading interp;colormap(hot)hold oncontour(tf,3,'k')hold offcaxis(k)colorbar function t1=dieda

9、i1(t) for i=2:5 for j=2:15 t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4; end t(i,16)=(2*t(i,15)+t(i-1,16)+t(i+1,16)/4;endfor i=6:11 for j=2:5 t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4; endendfor j=2:5 t(12,j)=(2*t(11,j)+t(12,j-1)+t(12,j+1)/4;endt1=t; function t1=diedai2(t)t(1,1)=(t(1,2)+t(2,1)+

10、2*10*0.1*30/0.53)/(2*(10*0.1/0.53+1);for j=2:15 t(1,j)=(2*t(2,j)+t(1,j-1)+t(1,j+1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);endt(1,16)=(2*t(2,16)+2*t(1,15)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);for i=2:11 t(i,1)=(2*t(i,2)+t(i-1,1)+t(i+1,1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);endt(12,1)=(2*t(12,2)+2*t(11

11、,1)+2*10*0.1*30/0.53)/(2*(10*0.1/0.53+2);for i=2:5 for j=2:15 t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4; end t(i,16)=(2*t(i,15)+t(i-1,16)+t(i+1,16)/4;endt(6,6)=(2*(t(6,5)+t(5,6)+t(7,6)+t(6,7)+2*4*0.1*10/0.53)/(2*(4*0.1/0.53+3);for j=7:15 t(6,j)=(2*t(5,j)+t(6,j-1)+t(6,j+1)+2*4*0.1*10/0.53)/(2*(4*

12、0.1/0.53+2);endt(6,16)=(2*t(5,16)+2*t(6,15)+2*4*0.1*10/0.53)/(2*(4*0.1/0.53+2);for i=7:11 t(i,6)=(2*t(i,5)+t(i-1,6)+t(i+1,6)+2*10*0.1*4/0.53)/(2*(4*0.1/0.53+2);endt(12,6)=(2*t(12,5)+2*t(11,6)+2*4*0.1*10/0.53)/(2*(4*0.1/0.53+2);for i=6:11 for j=2:5 t(i,j)=(t(i,j-1)+t(i,j+1)+t(i-1,j)+t(i+1,j)/4; endendfor j=2:5 t(12,j)=(2*t(11