偏微分方程组解法

《偏微分方程组解法》由会员分享，可在线阅读，更多相关《偏微分方程组解法（6页珍藏版）》请在装配图网上搜索。

1、-偏微分方程组解法 *厚度为10cm平壁原温度为20，现其两侧面分别维持在20和120，试求经过8秒后平壁温度分布，并分析温度分布随时间的变化直至温度分布稳定为止。式中为导温系数，；。 解： 模型转化为标准形式：初始条件为：边界条件为：，函数：pdefun.m%偏微分方程(一维动态传热)function c,f,s=pdefun(*,t,u,dud*)c=1/2e-4;f=dud*;s=0;icbun.m%偏微分方程初始条件(一维动态传热)function u0=icbun(*)u0=20;bcfun.m%偏微分方程边界条件(一维动态传热)function pl,ql,pr,qr=bcfun(

2、*l,ul,*r,ur,t)pl=ul-120;ql=0;pr=ur-20;qr=0;命令：*=linspace(0,10,20)*1e-2;t=linspace(0,15,16);sol=pdepe(0,pdefun,icfun,bcfun,*,t);mesh(*,t,sol(:,:,1) %温度与时间和空间位置的关系图%画1、2、4、6、8、15s时刻温度分布图plot(*,sol(2,:,1) 1s时刻，因为此题sol第一行为0时刻hold onplot(*,sol(3,:,1)plot(*,sol(5,:,1)plot(*,sol(7,:,1)plot(*,sol(9,:,1)plot

3、(*,sol(16,:,1)计算结果：%第8秒时温度分布*sol(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.587

4、726.298223.120720.0000或者求第8秒时，*=0,2,4,6,8,10cm处的温度uout,duoutd*=pdeval(0,*,sol(9,:,:),0,2,4,6,8,10*1e-2)120.0000 92.2279 67.5007 47.5765 32.3511 20.0000不同时刻温度分布图 将上图的视角转至*t平面也得到本图，从本图可知当时间到达15s时平壁的温度分布已近稳定。 *厚度为20cm钢板原温度为20，现将其置于1000的炉中加热，平壁导热系数为，对流传热系数，导温系数为；试分析温度分布随时间的变化及钢板外表温度到达500时所需的时间。 解： 模型转化为

5、标准形式：初始条件为：边界条件为：平壁中心坐标为0，绝热，函数：pdefun1.m%偏微分方程(一维动态平壁两侧对流)function c,f,s=pdefun1(*,t,u,dud*)c=1/0.555e-5;f=dud*;s=0;icbun1.m%偏微分方程初始条件(一维动态平壁两侧对流)function u0=icbun1(*)u0=20;bcfun1.m%偏微分方程边界条件(一维动态平壁两侧对流)function pl,ql,pr,qr=bcfun1(*l,ul,*r,ur,t)pl=0;ql=1;pr=174*(ur-1000);qr=34.8; %平壁两侧置于同一流体中具有对流传热

6、,平壁中心为绝热命令：%600s的温度分布变化*=linspace(0,10,20)*1e-2;t=0:60:600;sol=pdepe(0,pdefun1,icfun1,bcfun1,*,t);mesh(*,t,sol)%2160s的温度分布变化t=0:60:2160;sol=pdepe(0,pdefun1,icfun1,bcfun1,*,t);mesh(*,t,sol)%60、120、180、240、300、360、420s时刻温度分布图plot(*,sol(2,:,1) 60s(t网格为0:60:2160,其时间间隔为0,60,120,180,2160,第二点为60s)hold onpl

7、ot(*,sol(3,:,1)plot(*,sol(4,:,1)plot(*,sol(5,:,1)plot(*,sol(6,:,1)plot(*,sol(7,:,1)plot(*,sol(8,:,1)%1080、1440、1800、2160s时刻温度分布图plot(*,sol(19,:,1)hold onplot(*,sol(25,:,1)plot(*,sol(31,:,1)plot(*,sol(37,:,1)600s的温度分布变化2160s的温度分布变化不同时刻温度分布图，根据以上两准数可知：该传热过程外对流及导热阻力相当；当加热时间小于360s时为非正规阶段，加热时间大于360s后进入正规

8、阶段，从上图也可得到该结论。由该图可知，当加热时间到达2160s时，钢板的外表温度到达500，钢板中心温度为：uout,duoutd*=pdeval(0,*,sol(37,:,:),0*1e-2)371.2850有关传热量问题：钢板的外表温度到达500时，所需的总传热量，对平板：参见传热学，MatLab解法：*=0:1:10*1e-2;t=0:60:2160; 外表温度到达500时所需时间为2160ssol=pdepe(0,pdefun1,icfun1,bcfun1,*,t);size(sol)length(t)q=0;for i=1:36q=q+174*(sol(i,11,1)+sol(i+

9、1,11,1)/2-1000)*60;所需的热量均是从外表传递进入钢板的endq结果： 有一直径为40cm钢锭温度为20，将其置于900的炉中加热，平壁导热系数为，对流传热系数，导温系数为；试分析温度分布随时间的变化及钢锭外表温度加热到750时所需的时间；钢锭可近似为无限长的圆柱。初始条件为：边界条件为：平壁中心坐标为0，绝热，函数：pdefun2.micbun2.mbcfun2.m命令：*=0:1:20*1e-2;t=0:120:5520; (结合图5520s时，钢锭外表温度到达750左右)sol=pdepe(1,pdefun2,icfun2,bcfun2,*,t); m=1(圆柱)mesh

10、(*,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处；命令：*=-0.5:0.05:0.

11、5;y=-0.5:0.05:0.5;u*y=tri2grid(p,t,u,*,y);interp2(*,y,u*y,-0.4,-0.4) 结果：29.9010interp2(*,y,u*y,-0.4,0.4) 结果：29.9004interp2(*,y,u*y, 0.4,-0.4) 结果：29.9015interp2(*,y,u*y, 0.4,0.4) 结果：29.9063故其温度为29.9根据此题其四面对称，所以计算其中一面；传递的热量在外表处难以计算圆形的外表，但传递的热量必然通过外外表，其传热量只需将外表处各节点的一阶导数与节点间的面积、导热系数相乘即可得到传热量。命令：*=-0.5:0

12、.05:0.5;y=-0.5:0.05:0.5;u*y=tri2grid(p,t,u,*,y);dud*=( u*y(:,2)-u*y(:,1)/(*(2)-*(1) *方向一阶导数q=1.5*0.05*sum(dud*)*4 结果：675.0801dud*=( u*y(2,:)-u*y(1,:)/(y(2)-y(1) y方向一阶导数q=1.5*0.05*sum(dud*)*4 结果：675.0683由以上命令得*方向一阶导数：0 26.1727 47.8501 74.5754 99.8289 120.8720 141.5859160.9858 175.3053 184.6430 188.60

13、17 184.9701 174.9577 160.6397141.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根据传

14、热学形状因子的计算方法得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,

15、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(:,81)./u) 0.0047可见在8000s以后的温度分布已根本没有变化。h=1 r=0 c=1 a=0 f=1误差plot Height3D中选择user entry : u-(1-*.2-y.2)/4-*-yh=1 r=0 c=1 a=0 f=1-*-y需采用solve parameters中use nonlinear solver方可收敛. z.