偏微分方程组解法.doc

偏微分方程组解法 某厚度为10cm平壁原温度为20，现其两侧面分别维持在20和120，试求经过8秒后平壁内温度分布，并分析温度分布随时间的变化直至温度分布稳定为止。式中为导温系数，；。 解： 模型转化为标准形式：初始条件为：边界条件为：，函数：pdefun.m%偏微分方程(一维动态传热)function c,f,s=pdefun(x,t,u,dudx)c=1/2e-4;f=dudx;s=0;icbun.m%偏微分方程初始条件(一维动态传热)function u0=icbun(x)u0=20;bcfun.m%偏微分方程边界条件(一维动态传热)function pl,ql,pr,qr=bcfun(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) %温度与时间和空间位置的关系图%画1、2、4、6、8、15s时刻温度分布图plot(x,sol(2,:,1) 1s时刻，（因为本题sol第一行为0时刻）hold onplot(x,sol(3,:,1)plot(x,sol(5,:,1)plot(x,sol(7,:,1)plot(x,sol(9,:,1)plot(x,sol(16,:,1)计算结果：%第8秒时温度分布xsol(9,:,1)经过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.587726.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不同时刻温度分布图 将上图的视角转至xt平面也得到本图，从本图可知当时间达到15s时平壁内的温度分布已近稳定。 某厚度为20cm钢板原温度为20，现将其置于1000的炉中加热，平壁导热系数为，对流传热系数，导温系数为；试分析温度分布随时间的变化及钢板表面温度达到500时所需的时间。 解： 模型转化为标准形式：初始条件为：边界条件为：（平壁中心坐标为0，绝热），函数：pdefun1.m%偏微分方程(一维动态平壁两侧对流)function c,f,s=pdefun1(x,t,u,dudx)c=1/0.555e-5;f=dudx;s=0;icbun1.m%偏微分方程初始条件(一维动态平壁两侧对流)function u0=icbun1(x)u0=20;bcfun1.m%偏微分方程边界条件(一维动态平壁两侧对流)function pl,ql,pr,qr=bcfun1(xl,ul,xr,ur,t)pl=0;ql=1;pr=174*(ur-1000);qr=34.8; %平壁两侧置于同一流体中具有对流传热,平壁中心为绝热命令：%600s内的温度分布变化x=linspace(0,10,20)*1e-2;t=0:60:600;sol=pdepe(0,pdefun1,icfun1,bcfun1,x,t);mesh(x,t,sol)%2160s内的温度分布变化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网格为0:60:2160,其时间间隔为0,60,120,180,2160,第二点为60s)hold 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内的温度分布变化2160s内的温度分布变化不同时刻温度分布图，根据以上两准数可知：该传热过程内外对流及导热阻力相当；当加热时间小于360s时为非正规阶段，加热时间大于360s后进入正规阶段，从上图也可得到该结论。由该图可知，当加热时间达到2160s时，钢板的表面温度达到500，钢板中心温度为：uout,duoutdx=pdeval(0,x,sol(37,:,:),0*1e-2)371.2850有关传热量问题：（钢板的表面温度达到500时，所需的总传热量），对平板：（参见传热学），MatLab解法：x=0:1:10*1e-2;t=0:60:2160; （表面温度达到500时所需时间为2160s）sol=pdepe(0,pdefun1,icfun1,bcfun1,x,t);size(sol)length(t)q=0;for i=1:36q=q+174*(sol(i,11,1)+sol(i+1,11,1)/2-1000)*60;（所需的热量均是从表面传递进入钢板的）endq结果： 有一直径为40cm钢锭温度为20，将其置于900的炉中加热，平壁导热系数为，对流传热系数，导温系数为；试分析温度分布随时间的变化及钢锭表面温度加热到750时所需的时间；钢锭可近似为无限长的圆柱。初始条件为：边界条件为：（平壁中心坐标为0，绝热），函数：pdefun2.micbun2.mbcfun2.m命令：x=0:1:20*1e-2;t=0:120:5520; (结合图5520s时，钢锭表面温度达到750左右)sol=pdepe(1,pdefun2,icfun2,bcfun2,x,t); m=1(圆柱)mesh(x,t,sol)PDETOOL工具求解二维稳态与动态PDE： 如图偏心环形空间内表面温度为100，外表面温度为20；试给出其温度分布。 如图导热物体，下表面温度为20，上部三角形截面处温度为100，其余各面绝热，试给出其温度分布。 有一砖砌的烟气通道，其截面形状如附图所示，边长为1m，内管道直径为0.5m。已知内外壁温分别为80、25，砖的导热系数为1.5，试确定该通道的温度分布、距离任意相邻两直角边各0.1m处的温度及每米长烟道上的散热量。 解：距离相邻两直角边0.1m处的温度，即以上图中坐标（0.4,0.4）、（-0.4,0.4）、（0.4,-0.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) 结果：29.9010interp2(x,y,uxy,-0.4, 0.4) 结果：29.9004interp2(x,y,uxy, 0.4,-0.4) 结果：29.9015interp2(x,y,uxy, 0.4, 0.4) 结果：29.9063故其温度为29.9根据本题其四面对称，所以计算其中一面；传递的热量在内表面处难以计算（圆形的表面），但传递的热量必然通过外表面，其传热量只需将表面处各节点的一阶导数与节点间的面积、导热系数相乘即可得到传热量。命令：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方向一阶导数q=1.5*0.05*sum(dudx)*4 结果：675.0801dudx=( uxy(2,:)-uxy(1,:)/(y(2)-y(1) y方向一阶导数q=1.5*0.05*sum(dudx)*4 结果：675.0683由以上命令得x方向一阶导数：0 26.1727 47.8501 74.5754 99.8289 120.8720 141.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方向一阶导数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由以上边界处的一阶导数结合图可知在边界中心处的热通量最大；每米烟道的传热量为675W/m（根据传热学形状因子的计算方法得672.8W/m） 同上题，烟气通道原温度为25，零时刻开始内外壁温分别维持在为80、25，砖的导温系数为0.00001，试确定该通道温度分布随时间的变化及温度分布稳定所需的时间。 解：pdetool工具各菜单中的参数设置略，其中solve菜单parameters选项的时间设置为0:100:8000。 首先在稳态情况下（椭圆形）求解（三角形网格经过两次细化），导出解u；而后解以上非稳态（抛物形）并导出解u1（三角形网格细化与稳态相同）和三角形网格参数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稳定时间（以某时刻的温度分布与稳态温度分布的相对误差来判断）norm(u-u1(:,21)./u) 结果：2.4511n