偏微分方程組解法

1、偏微分方程組解法 某厚度為10cm平壁原溫度為20，現(xiàn)其兩側(cè)面分別維持在20和120，試求經(jīng)過8秒后平壁內(nèi)溫度分布，并分析溫度分布隨時間的變化直至溫度分布穩(wěn)定為止。式中為導(dǎo)溫系數(shù)，；。 解： 模型轉(zhuǎn)化為標(biāo)準(zhǔn)形式：初始條件為：邊界條件為：，函數(shù)：pdefun.m%偏微分方程(一維動態(tài)傳熱)function c,f,s=pdefun(x,t,u,dudx)c=1/2e-4;f=dudx;s=0;icbun.m%偏微分方程初始條件(一維動態(tài)傳熱)function u0=icbun(x)u0=20;bcfun.m%偏微分方程邊界條件(一維動態(tài)傳熱)function pl,ql,pr,qr=bcfun(

2、xl,ul,xr,ur,t)pl=ul-120;ql=0;pr=ur-20;qr=0;命令：x=linspace(0,10,20)*1e-2;t=linspace(0,15,16);sol=pdepe(0,pdefun,icfun,bcfun,x,t);mesh(x,t,sol(:,:,1) %溫度與時間和空間位置的關(guān)系圖%畫1、2、4、6、8、15s時刻溫度分布圖plot(x,sol(2,:,1) 1s時刻，（因?yàn)楸绢}sol第一行為0時刻）hold onplot(x,sol(3,:,1)plot(x,sol(5,:,1)plot(x,sol(7,:,1)plot(x,sol(9,:,1)pl

3、ot(x,sol(16,:,1)計算結(jié)果：%第8秒時溫度分布xsol(9,:,1)經(jīng)過8秒時的溫度分布為：/cm00.52631.05261.57892.10532.63163.1579/120.0000112.5520105.165397.899490.810083.947777.3562/cm3.68424.21054.73685.26325.78956.31586.8421/71.071465.120259.520054.278449.393044.851840.6338/cm7.36847.89478.42118.94749.473710.0000/36.709533.041929.5

4、87726.298223.120720.0000或者求第8秒時，x=0,2,4,6,8,10cm處的溫度uout,duoutdx=pdeval(0,x,sol(9,:,:),0,2,4,6,8,10*1e-2)120.0000 92.2279 67.5007 47.5765 32.3511 20.0000不同時刻溫度分布圖 將上圖的視角轉(zhuǎn)至xt平面也得到本圖，從本圖可知當(dāng)時間達(dá)到15s時平壁內(nèi)的溫度分布已近穩(wěn)定。 某厚度為20cm鋼板原溫度為20，現(xiàn)將其置于1000的爐中加熱，平壁導(dǎo)熱系數(shù)為，對流傳熱系數(shù)，導(dǎo)溫系數(shù)為；試分析溫度分布隨時間的變化及鋼板表面溫度達(dá)到500時所需的時間。 解： 模型

5、轉(zhuǎn)化為標(biāo)準(zhǔn)形式：初始條件為：邊界條件為：（平壁中心坐標(biāo)為0，絕熱），函數(shù)：pdefun1.m%偏微分方程(一維動態(tài)平壁兩側(cè)對流)function c,f,s=pdefun1(x,t,u,dudx)c=1/0.555e-5;f=dudx;s=0;icbun1.m%偏微分方程初始條件(一維動態(tài)平壁兩側(cè)對流)function u0=icbun1(x)u0=20;bcfun1.m%偏微分方程邊界條件(一維動態(tài)平壁兩側(cè)對流)function pl,ql,pr,qr=bcfun1(xl,ul,xr,ur,t)pl=0;ql=1;pr=174*(ur-1000);qr=34.8; %平壁兩側(cè)置于同一流體中具

6、有對流傳熱,平壁中心為絕熱命令：%600s內(nèi)的溫度分布變化x=linspace(0,10,20)*1e-2;t=0:60:600;sol=pdepe(0,pdefun1,icfun1,bcfun1,x,t);mesh(x,t,sol)%2160s內(nèi)的溫度分布變化t=0:60:2160;sol=pdepe(0,pdefun1,icfun1,bcfun1,x,t);mesh(x,t,sol)%60、120、180、240、300、360、420s時刻溫度分布圖plot(x,sol(2,:,1) 60s(t網(wǎng)格為0:60:2160,其時間間隔為0,60,120,180,2160,第二點(diǎn)為60s)ho

7、ld onplot(x,sol(3,:,1)plot(x,sol(4,:,1)plot(x,sol(5,:,1)plot(x,sol(6,:,1)plot(x,sol(7,:,1)plot(x,sol(8,:,1)%1080、1440、1800、2160s時刻溫度分布圖plot(x,sol(19,:,1)hold onplot(x,sol(25,:,1)plot(x,sol(31,:,1)plot(x,sol(37,:,1)600s內(nèi)的溫度分布變化2160s內(nèi)的溫度分布變化不同時刻溫度分布圖，根據(jù)以上兩準(zhǔn)數(shù)可知：該傳熱過程內(nèi)外對流及導(dǎo)熱阻力相當(dāng)；當(dāng)加熱時間小于360s時為非正規(guī)階段，加熱時間大

8、于360s后進(jìn)入正規(guī)階段，從上圖也可得到該結(jié)論。由該圖可知，當(dāng)加熱時間達(dá)到2160s時，鋼板的表面溫度達(dá)到500，鋼板中心溫度為：uout,duoutdx=pdeval(0,x,sol(37,:,:),0*1e-2)371.2850有關(guān)傳熱量問題：（鋼板的表面溫度達(dá)到500時，所需的總傳熱量），對平板：（參見傳熱學(xué)），MatLab解法：x=0:1:10*1e-2;t=0:60:2160; （表面溫度達(dá)到500時所需時間為2160s）sol=pdepe(0,pdefun1,icfun1,bcfun1,x,t);size(sol)length(t)q=0;for i=1:36q=q+174*(so

9、l(i,11,1)+sol(i+1,11,1)/2-1000)*60;（所需的熱量均是從表面?zhèn)鬟f進(jìn)入鋼板的）endq結(jié)果： 有一直徑為40cm鋼錠溫度為20，將其置于900的爐中加熱，平壁導(dǎo)熱系數(shù)為，對流傳熱系數(shù)，導(dǎo)溫系數(shù)為；試分析溫度分布隨時間的變化及鋼錠表面溫度加熱到750時所需的時間；鋼錠可近似為無限長的圓柱。初始條件為：邊界條件為：（平壁中心坐標(biāo)為0，絕熱），函數(shù)：pdefun2.micbun2.mbcfun2.m命令：x=0:1:20*1e-2;t=0:120:5520; (結(jié)合圖5520s時，鋼錠表面溫度達(dá)到750左右)sol=pdepe(1,pdefun2,icfun2,bcfu

10、n2,x,t); m=1(圓柱)mesh(x,t,sol)PDETOOL工具求解二維穩(wěn)態(tài)與動態(tài)PDE： 如圖偏心環(huán)形空間內(nèi)表面溫度為100，外表面溫度為20；試給出其溫度分布。 如圖導(dǎo)熱物體，下表面溫度為20，上部三角形截面處溫度為100，其余各面絕熱，試給出其溫度分布。 有一磚砌的煙氣通道，其截面形狀如附圖所示，邊長為1m，內(nèi)管道直徑為0.5m。已知內(nèi)外壁溫分別為80、25，磚的導(dǎo)熱系數(shù)為1.5，試確定該通道的溫度分布、距離任意相鄰兩直角邊各0.1m處的溫度及每米長煙道上的散熱量。 解：距離相鄰兩直角邊0.1m處的溫度，即以上圖中坐標(biāo)（0.4,0.4）、（-0.4,0.4）、（0.4,-0.

11、4）或（-0.4,-0.4）處；命令：x=-0.5:0.05:0.5;y=-0.5:0.05:0.5;uxy=tri2grid(p,t,u,x,y);interp2(x,y,uxy,-0.4,-0.4) 結(jié)果：29.9010interp2(x,y,uxy,-0.4, 0.4) 結(jié)果：29.9004interp2(x,y,uxy, 0.4,-0.4) 結(jié)果：29.9015interp2(x,y,uxy, 0.4, 0.4) 結(jié)果：29.9063故其溫度為29.9根據(jù)本題其四面對稱，所以計算其中一面；傳遞的熱量在內(nèi)表面處難以計算（圓形的表面），但傳遞的熱量必然通過外表面，其傳熱量只需將表面處各節(jié)點(diǎn)

12、的一階導(dǎo)數(shù)與節(jié)點(diǎn)間的面積、導(dǎo)熱系數(shù)相乘即可得到傳熱量。命令：x=-0.5:0.05:0.5;y=-0.5:0.05:0.5;uxy=tri2grid(p,t,u,x,y);dudx=( uxy(:,2)-uxy(:,1)/(x(2)-x(1) x方向一階導(dǎo)數(shù)q=1.5*0.05*sum(dudx)*4 結(jié)果：675.0801dudx=( uxy(2,:)-uxy(1,:)/(y(2)-y(1) y方向一階導(dǎo)數(shù)q=1.5*0.05*sum(dudx)*4 結(jié)果：675.0683由以上命令得x方向一階導(dǎo)數(shù)：0 26.1727 47.8501 74.5754 99.8289 120.8720 141

13、.5859 160.9858 175.3053 184.6430 188.6017 184.9701 174.9577 160.6397 141.7009 120.5195 97.6419 72.9719 50.2694 26.1751 0y方向一階導(dǎo)數(shù)0 26.1727 50.2649 72.9663 97.6365 120.5146 141.6968 160.6367 174.9557 184.9691 188.6015 184.6433 175.3060 160.9868 141.5868 120.8722 99.8279 74.5720 47.8485 26.1696 0由以上邊界處的

14、一階導(dǎo)數(shù)結(jié)合圖可知在邊界中心處的熱通量最大；每米煙道的傳熱量為675W/m（根據(jù)傳熱學(xué)形狀因子的計算方法得672.8W/m） 同上題，煙氣通道原溫度為25，零時刻開始內(nèi)外壁溫分別維持在為80、25，磚的導(dǎo)溫系數(shù)為0.00001，試確定該通道溫度分布隨時間的變化及溫度分布穩(wěn)定所需的時間。 解：pdetool工具各菜單中的參數(shù)設(shè)置略，其中solve菜單parameters選項(xiàng)的時間設(shè)置為0:100:8000。 首先在穩(wěn)態(tài)情況下（橢圓形）求解（三角形網(wǎng)格經(jīng)過兩次細(xì)化），導(dǎo)出解u；而后解以上非穩(wěn)態(tài)（拋物形）并導(dǎo)出解u1（三角形網(wǎng)格細(xì)化與穩(wěn)態(tài)相同）和三角形網(wǎng)格參數(shù)p，e，t；命令：pdemesh(p,e,t,u1(:,1) 0spdemesh(p,e,t,u1(:,2) 100spdemesh(p,e,t,u1(:,4) 300spdemesh(p,e,t,u1(:,7) 600spdemesh(p,e,t,u1(:,11) 1000spdemesh(p,e,t,u1(:,21) 2000spdemesh(p,e,t,u1(:,31) 3000spdemesh(p,e,t,u1(:,51) 5000spdemesh(p,e,t,u1(:,81) 8000s穩(wěn)定時間（以某時刻的溫度分布與穩(wěn)態(tài)溫度分布的相對誤差來判斷）norm(u-u1(:,21)./u) 結(jié)果：2.4511n