偏微分方程组解法

1、偏微分方程组解法试求某厚度为10cm平壁原温度为20 C,现其两侧面分别维持在20 C和120 C,经过8秒后平壁温度分布，并分析温度分布随时间的变化直至温度分布稳定为止。t2t式中a为导温系数，cm2/s; a 2。解：模型转化为标准形式：2t2 t a初始条件为：t x, 020边界条件为：t 0,120, t 0.1,20函数：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

2、.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时刻温度分布图Word文档 plot(x,sol(2,:,1) 1s时刻，(因为本题sol第一行为0时亥U) hold onplot(x,sol(3,

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秒时的温度分布为:x/cm00.52631.05261.57892.10532.63163.1579t/ C120.0000112.5520105.165397.899490.810083.947777.3562x/cm3.68424.21054.73685.26325.78956.31586.8421t/ C71.071465.120259.520054.278449.393044.85184

4、0.6338x/cm7.36847.89478.42118.94749.473710.0000t/ C36.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.000092.227967.500747.576532.351120.00000 DWord文档不同时刻温度分布图将上图的视角转至 xt平面也得到本图，从本图可知当时间达到 15s时平壁的温度分布 已近稳定。Word文档某厚度为20cm钢板原温

5、度为20 C ,现将其置于1000 C的炉中加热，平壁导热系数为34.8W/m C ,对流传热系数 h 174W/m2 C ,导温系数为a 0.555 10 5 m2/s ;试分析温度分布随时间的变化及钢板表面温度达到I 500 C时所需的时间。2ta 2X解：模型转化为标准形式:2tX初始条件为:t x, 020边界条件为：t 0,0 (平壁中心坐标为 0,绝热)，h t 0.1,t 0.1,函数：pdefun1.m%偏微分方程(一维动态平壁两侧对流)function c,f,s=pdefun1(x,t,u,dudx)c=1/0.555e-5;f=dudx;s=0;icbun1.m%偏微分方

6、程初始条件(一维动态平壁两侧对流)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)Word文档%2160s的温度分布变

7、化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、2

8、160s时刻温度分布图plot(x,sol(19,:,1)hold onplot(x,sol(25,:,1)plot(x,sol(31,:,1)plot(x,sol(37,:,1)600s的温度分布变化Word文档不同时刻温度分布图Word文档0.2360s时为非由该图可知，当加热时间达到 2160s时，钢板的表面温度达到500 C,钢板中心温度为:uout,duoutdx=pdeval(0,x,sol(37,:,:),0*1e-2)371.2850 C有关传热量问题：(钢板的表面温度达到500 C时，所需的总传热量)Q QQo cV t0 t1 ， A exp0012Fo B对平板：(参见传

9、热学)sin 1sin 11 0.66818, A 21 1.07329, B 1 0.92723 HYPERLINK l bookmark67 o Current Document 1 sin 1 cos 112 _ 0 5352_ _一 Aexp ；Fo B 1.07329 e 0.92723 0.58270Word文档QocV t0 t1 0.5827 0.41726h 174 0.1a 0.555 10 5 360Bi 0.5 , Fo - z34.820.12根据以上两准数可知：该传热过程外对流及导热阻力相当；当加热时间小于正规阶段，加热时间大于 360s后进入正规阶段，从上图也可得

10、到该结论。1748Qo cV t0 t 1 0.120 1 0006.14486 1 08 J0.555 10 5Q 0.41726 6.14486 1082.564 108 JMatLab 解法：x=0:1:10*1e-2;t=0:60:2160;(表面温度达到 500 C时所需时间为 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结果：Q 2.47

11、97 108 JWord文档有一直径为40cm钢锭温度为20 C,将其置于900 C的炉中加热，平壁导热系数为34.8W/m C,对流传热系数h 174W/m2 C ,导温系数为a 0.695 10 5 m2/s；试分析温度分布随时间的变化及钢锭表面温度加热到750 C时所需的时间;长的圆柱。t 1 tax x x x初始条件为：t x, 020边界条件为：-t-02 0 （平壁中心坐标为 0,绝热），ht0.2, t x函数：钢锭可近似为无限t 0.2, xpdefun2.m icbun2.m bcfun2.m命令：x=0:1:20*1e-2;t=0:120:5520;(结合图5520s时，

12、钢锭表面温度达到750 C左右)sol=pdepe(1,pdefun2,icfun2,bcfun2,x,t); m=1(圆柱)Word文档mesh(x,t,sol)PDETOOL工具求解二维稳态与动态PDE :如图偏心环形空间表面温度为100 C ,外表面温度为20 C ;试给出其温度分布。如图导热物体，下表面温度为 热，试给出其温度分布。20 C,上部三角形截面处温度为100 C,其余各面绝Word文档有一砖砌的烟气通道，其截面形状如附图所示，边长为 1m,管道直径为0.5m。已知 外壁温分别为80 C、25 C,砖的导热系数为1.5 W/m C ,试确定该通道的温度分布、 距离任意相邻两直

13、角边各0.1m处的温度及每米长烟道上的散热量。Word文档距离相邻两直角边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);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.9063interp2(x,y,uxy, 0.4, 0.4)结果

14、：故其温度为29.9 Cn根据本题其四面对称，所以计算其中一面；传递的热量在表面处难以计算(圆形的表面)，但传递的热量必然通过外表面，其传热量只需将表面处各节点的一阶导数与节点间的面积、导热系数相乘即可得到传热量。命令：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(

15、dudx)*4结果：675.0683由以上命令得x方向一阶导数：026.172747.850174.575499.8289 120.8720 141.5859 160.9858 175.3053Word文档184.6430188.6017184.9701174.9577160.6397141.7009120.519597.641972.971950.269426.17510y方向一阶导数026.172750.264972.966397.6365 120.5146 141.6968 160.6367 174.9557184.9691188.6015184.6433175.3060160.9868

16、141.5868120.872299.827974.572047.848526.16960由以上边界处的一阶导数结合图可知在边界中心处的热通量最大；每米烟道的传热量为675W/m （根据传热学形状因子的计算方法得672.8W/m）同上题，烟气通道原温度为25 C ,零时刻开始外壁温分别维持在为80 C、25 C ,2t-2 ysoke菜单parameters选项的时间设置醇的导温系数为 0.00001 W/m C，试确定该通道温度分布随时间的变化及温度分布稳定 所需的时间。t2t- a -2X解：pdetool工具各菜单中的参数设置略，其中 为 0:100:8000。首先在稳态情况下（椭圆形）

17、求解（三角形网格经过两次细化），导出解u;而后解以上非稳态（抛物形）并导出解 u1 （三角形网格细化与稳态相同）和三角形网格参数p, e, t；命令:Word文档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.4511norm(u-u1(:,31)./u)0.8079norm(u-u1(:,51)./u)0.0770norm(u-u1(