【正文】 。在正規(guī)階段可以用無窮級數(shù)的第一項作為近似解，但其計算涉及到求解隱式的超越方程，過程十分繁瑣；故工程上利用哈斯勒圖進(jìn)行計算，但結(jié)果誤差很大；而對非正規(guī)階段，則必須依次計算無窮級數(shù)中的各項直至達(dá)到精度為止，計算量大。若采用MATLAB求解該問題則十分快捷方便。根據(jù)拋物型微分方程的標(biāo)準(zhǔn)式： （9）式中u=u(t,x)與t=t(τ,x)對應(yīng)，m為形狀因子，m=0,1,2分別為對應(yīng)無限大平板、柱形和球形；function[c,f,s]=pdefun(x,t,u,dudx)c=1/。f=dudx。s=0。function u0=icbun(x)u0=20。function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)pl=0,ql=1。pr=174*(ur1000)。qr=。而后對空間和時間劃分網(wǎng)格，應(yīng)用pdepe命令求解該偏微分方程，其解結(jié)果存放在sol中，命令如下：x=linspace(0,20)。t=[0:60:2160]。sol=pdepe(0,@pdefun,@icfun,@bcfun,x,t)。mesh(x,t,sol)對一個偏微分方程sol為mnl矩陣，其中第一維與時間網(wǎng)格對應(yīng)，第二維與空間網(wǎng)格對應(yīng)；用mesh(x,t,sol)命令可以繪出空間、時間與偏微分方程解（對本問題來講就是不同時刻不同位置的溫度）的關(guān)系如圖4所示。而用plot(x,sol(m,:,l))可以畫出在60m60(60為時間網(wǎng)格步長)時刻導(dǎo)熱物體的溫度分布，如圖5所示。 以上時間網(wǎng)格設(shè)置為t=[0:60:2160],所以在圖4中x方向的各條線與不同的時刻0、60s、120s、……2100s、2160s共37個時刻對應(yīng)；從圖中可知當(dāng)傳熱時間達(dá)到500℃；而60s60時刻x0處的溫度可用[uout,duoutdx]=pdeval(0,x,sol(m+l,:,:),x0)命令求解，其中uout為溫度，duoutdx為對應(yīng)的溫度梯度；如加熱時間達(dá)到2160s時，執(zhí)行命令[uout,duoutdx]=pdeval(0,x,sol(37,:,:),0)℃。而從圖5可知鋼板加熱時間小于240s時，其溫度分布曲線隨時間而變化，故該階段傳熱的非正規(guī)階段；當(dāng)加熱時間大于240s后，溫度分布形狀不再變化，進(jìn)入了正規(guī)階段。將鋼板的表面溫度加熱到500℃整個傳熱過程為動態(tài)的，但是鋼板得到的熱量均是通過表面?zhèn)鬟f到內(nèi)部，所以只要將每一時間網(wǎng)格內(nèi)的表面?zhèn)鳠崃肯嗉蛹纯傻玫娇倐鳠崃?，而每一時間網(wǎng)格內(nèi)的傳熱量只要時間間隔足夠短就可用下式計算： （10）其中，m為時間網(wǎng)格數(shù)，n應(yīng)為鋼板表面處的網(wǎng)格節(jié)點，由于空間網(wǎng)格劃分為x=linspace(0,20)共有21個節(jié)點，節(jié)點n=21即代表鋼板表面；由以下程序即可求得總傳熱量Q=108J。q=0。for m=1:36 %36為時間網(wǎng)格數(shù)（不含零時刻）q=q+174*((sol(m,21,1)+sol(m+1,21,1))/21000)*60。end對于二維動態(tài)導(dǎo)熱，同樣可用PEDTOOL工具箱求解。如有一磚砌變長為1m煙氣通道。已知煙氣通道原溫度為25℃，零時刻開始內(nèi)外壁溫分別維持在80℃、25℃，℃，℃,試確定:(1) 該通道的溫度分布隨時間的變化。(2) ； （3）傳熱穩(wěn)定后每米長煙道上的散熱量。本問題為二維動態(tài)導(dǎo)熱過程，其微分方程為： （11）通過solve菜單parameters選項可對時間網(wǎng)格及初始條件進(jìn)行設(shè)置，時間網(wǎng)格設(shè)置為0：100：8000即從零時刻開始計算到8000s，時間間隔為100s；初始溫度設(shè)置為25℃。將PDETOOL計算的網(wǎng)格數(shù)據(jù)p,e,t和因變量即溫度結(jié)果（保存于u矩陣中）導(dǎo)出，其中u為nm矩陣,其中第一維與網(wǎng)格節(jié)點對應(yīng),第二維與時間網(wǎng)格對應(yīng)。由命令pdemesh(p,e,t,u(:m))即可繪制出100m100s時刻的溫度分布圖，如圖6所示。其中x、y為煙道界面的坐標(biāo)，t為溫度；圖6中分別是t=0s、100s、300s、600s、1000s、2000s、5000s、8000s時刻的溫度分布圖，從圖中可看出溫度分布的變化情況，可見在5000s之后溫度分布已不再隨時間變化。若將8000s的溫度分布視為穩(wěn)定的溫度分布，%，而5000s后的溫度分布完全可以認(rèn)為進(jìn)入穩(wěn)定傳熱階段。那么要求其中任一點的溫度，則必須將三角形網(wǎng)格轉(zhuǎn)化為直角坐標(biāo)網(wǎng)格，由于煙道為邊長1m的方形管道，在利用PDETOOL求解時是以坐標(biāo)零點為中心，所以其求解域為x,y∈[,]。因此利用以下程序及網(wǎng)格轉(zhuǎn)化命令tri2grid將三角形網(wǎng)格p,t上解u轉(zhuǎn)化為直角網(wǎng)格x,y上的解uxy。x=[::]。y=[::]。uxy=tri2grid(p,t,u(:,81),x,y)。而后利用二維插值函數(shù)interp2求出對應(yīng)位置的溫度，（，），所以插值命令為interp2(x,y,uxy,),℃。傳熱穩(wěn)定后的散熱量若根據(jù)內(nèi)外兩面的溫度差來計算，由于其傳熱面和熱傳遞路徑復(fù)雜無法直接計算，傳熱學(xué)中其計算采用形狀因子法，針對不同形狀的導(dǎo)熱物體其形狀因子也不相同，在應(yīng)用中具有一定的局限性；而利用MATLAB數(shù)值模擬則克根據(jù)溫度分布直接求得傳熱量，因為傳遞的熱量必然通過外表面，其傳熱量只需將表面處各節(jié)點的一階倒數(shù)與節(jié)點間的傳熱面積、導(dǎo)熱系數(shù)相乘即可得到傳熱量。 （12）根據(jù)本題由于煙道外表面的四面對稱，計算其中一面再乘以4倍即為總傳熱量；命令如下：dudx=(uxy(2,:)uxy(1,:))/(x(2)x(1))q=**sum(dudx)*4根據(jù)表面處溫度分布求得的差商近似代替表面處的每一網(wǎng)格內(nèi)的溫度梯度dudx，㎡，由此計算出傳熱量為675W/m。4．小結(jié)傳熱學(xué)是一門基于理論計算和實驗的傳統(tǒng)學(xué)科，以往對傳熱學(xué)的研究中，實驗是一種最常用也最重要的手段。但是，由于傳熱學(xué)學(xué)科的自身參數(shù)多、傳熱情況復(fù)雜的特點，實驗的設(shè)計和進(jìn)行存在許多困難。相似原理的提出，雖對實驗研究的簡化起到了決定性的作用，但計算公式和修正系數(shù)的確定，任需要進(jìn)行大量實驗及煩瑣的計算。MATLAB等計算軟件的廣泛應(yīng)用為計算傳熱學(xué)的發(fā)展和傳熱學(xué)的研究提供了十分便利的條件。與以往的研究模式相比，實踐證明，使用MATLAB能方便高效地對已建立的模型及復(fù)雜的方程組進(jìn)行數(shù)值計算和模擬仿真，同時其強大的圖形句柄功能為使結(jié)果可視化、直觀化和形象化提供了方便。也可以對模擬結(jié)果用于實驗進(jìn)行驗證，這就大大縮短了研究時間，減少了實驗次數(shù)，這也為提高研究的經(jīng)濟性提供了有效的途徑。另一方面，MATLAB軟件在傳熱學(xué)研究中已經(jīng)得到了廣泛的應(yīng)用，同時促進(jìn)了本學(xué)科的發(fā)展。隨著MATLAB自身的不斷發(fā)展和完善下，筆者認(rèn)為，其在傳熱學(xué)的理論研究和工程應(yīng)用方面都將發(fā)揮更大的作用。參考文獻(xiàn)[1]許國良，王曉墨，鄔田華，陳維漢. 工程傳熱學(xué). 北京：中國電力出版社，[2]（德）夸特羅尼，（德）色拉瑞著；. 北京：清華大學(xué)出版社，：Scientific Computing with MATLAB(國外計算機科學(xué)經(jīng)典教材)[3]魏巍. MATLAB應(yīng)用數(shù)學(xué)工具箱技術(shù)手冊. 北京：國防工業(yè)出版社，.[4]祝用華，趙慧，浦健，張奕. ，：63～67.[5]，：986～990.[6]許波，. 北京：清華大學(xué)出版社，2000.[7]李燦，高彥棟，（自然科學(xué)版），.[8]姚東，王愛民，馮峰，：.[9]（第2版）.西安：西安交通大學(xué)出版社，2