傳熱學(xué)大作業(yè)

傳熱學(xué)第四章大作業(yè)二維穩(wěn)態(tài)導(dǎo)熱問題的數(shù)值解法姓名：班級：學(xué)號：第一題：如圖所示，一個無限長矩形柱體，其橫截面的邊長分別為和，常物性。該問題可視為二維穩(wěn)態(tài)導(dǎo)熱問題，邊界條件如圖中所示，其中=0.6m，=0.4m，=60℃，=20℃，。編寫程序求解二維導(dǎo)熱方程。繪制x=/2和y=/2處的溫度場，并與解析解進(jìn)行比較。已知矩形內(nèi)的溫度場的解析解為：解：(1)建立控制方程及定解條件控制方程：定解條件：（2）區(qū)域離散化（確立節(jié)點）將矩形區(qū)域分為M*N個網(wǎng)格，其中x方向上的步長;y方向的步長。設(shè)節(jié)點為（m,n）。（3）建立節(jié)點離散方程對節(jié)點（m,n）有：內(nèi)節(jié)點:化簡得邊界節(jié)點：（4）編程求解，程序見附錄取M=N=50得到矩形區(qū)域各節(jié)點溫度（見附件一），為方便在這里僅給出M=N=10時溫度分布數(shù)據(jù)，如下表：606060606060606060606060.4486560.8287961.0919261.213461.1883861.0274160.7531460.39773606060.9277461.7095862.2468162.491962.4371862.1051761.5423660.81427606061.4716862.6990863.5315463.9036663.8092363.285562.4050261.26912606062.1241863.8641265.0208665.5227665.3720864.6243963.381261.78318606062.9474565.2883666.8012267.4319167.2002766.1835164.5150962.37953606064.040367.0806368.9736469.7243169.377468.0322465.856663.08431606065.5775969.3828171.6537872.5045271.998270.2500867.4632863.92776606067.9025872.3668974.9651375.8878875.1709172.9295569.4027464.945536066.1803471.7557176.1803479.021138079.0211376.1803471.7557166.1803460M=N=50時矩形區(qū)域節(jié)點溫度分布圖如圖：M=N=50時x=/2處數(shù)值解與解析解溫度場分布如下：M=N=10時將數(shù)值解與理論解進(jìn)行比較，如下表所示：y00.4/90.4/9*20.4/9*30.4/9*40.4/9*50.4/9*60.4/9*70.4/9*80.4數(shù)值解6061.213462.491963.9036665.5227667.4319169.7243172.5045275.8878880理論解6061.1744762.4128263.7824465.3578167.2246469.4844972.2602975.7030380相對誤差00.064%0.127%0.190%0.252%0.308%0.345%0.338%0.244%0M=N=50時y=/2處數(shù)值解與解析解溫度場分布如下圖：M=N=10時將數(shù)值解與理論解進(jìn)行比較，如下表所示：x00.6/90.6/9*20.6/9*30.6/9*40.6/9*50.6/9*60.6/9*70.6/9*80.6數(shù)值解6062.1241863.8641265.0208665.5227665.3720864.6243963.381261.7831860理論解6062.1372464.016765.4116966.1539566.1539565.4116964.016762.1372460相對誤差00.021%0.238%0.598%0.954%1.182%1.204%0.993%0.569%0第二題：將第一題中處的邊界條件變?yōu)?，其他條件不變。編寫程序求解二維導(dǎo)熱方程并計算從y=0處導(dǎo)入的熱量。當(dāng)時，該二維導(dǎo)熱問題可簡化為一維導(dǎo)熱問題。在一維的近似下，試計算從y=0處導(dǎo)入的熱量，并比較不同L2/L1下的比值。由該問題的解析解可知：L2/L10.0070.010.050.080.10.99870.99120.9560.930.912解：編程得到M=N=50溫度分布數(shù)據(jù)見附件二，得到溫度分布圖如下：這里也僅給出M=N=10時溫度分布數(shù)據(jù)，見下表：606060606060606060606058.9248458.0307557.4169257.1105157.1105157.4169258.0307558.92484606057.7692255.9369354.697254.0848354.0848354.697255.9369357.76922606056.4364853.5797551.6986550.78750.78751.6986553.5797556.43648606054.7896250.7889648.2692447.0839947.0839948.2692450.7889654.78962606052.605147.3399744.2467242.854242.854244.2467247.3399752.6051606049.4740242.925739.4683438.0055138.0055139.4683442.925749.47402606044.5750837.1376833.8034832.5066732.5066733.8034837.1376844.57508606036.1261529.5260127.2331326.4314826.4314827.2331329.5260136.126156020202020202020202020從y=0處導(dǎo)入的熱量可以近似看作從y=0向y=0.4/9處傳遞的熱量，y=0.4/9處的溫度分別為：x00.6/90.6/9*20.6/9*30.6/9*40.6/9*50.6/9*60.6/9*70.6/9*80.6T6058.9248458.0307557.4169257.1105157.1105157.4169258.0307558.9248460則整理得當(dāng)時，該二維導(dǎo)熱問題可簡化為一維導(dǎo)熱問題。當(dāng)時，令M=N=100，得=1129571.429看作一維導(dǎo)熱時不同下比值如下表：L2/L10.0070.010.050.080.11129571.42978954015281292890729081142857.143800000160000100000800000.98840.98690.95510.92890.91140.99870.99120.9560.930.912附錄Matlab程序：第一題：functionchuanredazuoyeclearclcL1=0.6;%矩形長度L2=0.4;%矩形寬度Tw1=60;Tw2=20;globalMN%設(shè)置網(wǎng)格數(shù)M=input('請輸入將區(qū)間[0,L1]等分的個數(shù)M:');N=input('請輸入將區(qū)間[0,L2]等分的個數(shù)N:');XDIF=L1/M;%x方向上的步長YDIF=L2/N;%y方向上的步長axis([0,M,0,N]);grid%設(shè)置網(wǎng)格U=initial(M,N,Tw1,Tw2,L1,XDIF);AIM=YDIF/XDIF;AIP=YDIF/XDIF;AJM=XDIF/YDIF;AJP=XDIF/YDIF;AP=2*YDIF/XDIF+2*XDIF/YDIF;CON=0;%離散方程系數(shù)T=rechuandao(U,M,N,AIM,AIP,AJM,AJP,AP);A=flipud(T);mesh(A);title('溫度分布')xlabel('x')ylabel('y')zlabel('T')functionU=rechuandao(U,M,N,AIM,AIP,AJM,AJP,AP)while1temp=U;fori=2:M-1forj=2:N-1U(i,j)=(AIM*U(i,j-1)+AIP*U(i,j+1)+AJM*U(i-1,j)+AJP*U(i+1,j))/AP;%設(shè)置迭代條件endendeps=abs(U-temp);ifmax(max(eps))<1e-8break;%限制迭代次數(shù)endendfunctionU=initial(M,N,Tw1,Tw2,L1,XDIF)U=zeros(M,N);%賦溫度矩陣初值U(:,1)=Tw1;U(:,M)=Tw1;U(N,:)=Tw1;fori=1:MU(1,i)=Tw1+Tw2*sin(pi*i*XDIF/L1);end%初始和邊界條件的設(shè)定第二題：functiongaibianjietiaojianclearclcL1=0.6;%矩形長度L2=0.4;%矩形寬度Tw1=60;Tw2=20;globalMN%設(shè)置網(wǎng)格數(shù)M=input('請輸入將區(qū)間[0,L1]等分的個數(shù)M:');N=input('請輸入將區(qū)間[0,L2]等分的個數(shù)N:');XDIF=L1/M;%x方向上的步長YDIF=L2/N;%y方向上的步長axis([0,M,0,N]);grid%設(shè)置網(wǎng)格U=initial(M,N,Tw1,Tw2);AIM=YDIF/XDIF;AIP=YDIF/XDIF;AJM=XDIF/YDIF;AJP=XDIF/YDIF;AP=2*YDIF/XDIF+2*XDIF/YDIF;CON=0;%離散方程系數(shù)T=rechuandao(U,M,N,AIM,AIP,AJM,AJP,AP);A=flipud(T);mesh(A);title('改變邊界條件后溫度分布')xlabel('x')ylabel('y')zlabel('T')functionU=rechuandao(U,M,N,AIM,AIP,AJM,AJP,AP)while1temp=U;fori=2:M-1forj=2:N-1U(i,j)=(AIM*U(i,j-1)+AIP*U(i,j+1)+AJM*U(i-1,j)+AJP*U(i+1,j))/AP;%設(shè)置迭代條件endendeps=abs(U-t