偏微分方程组解法

偏微分方程组解法

某厚度为10cm平壁原温度为20C，现其两侧面分别维持在20C和120C，试求经过8秒后平壁温度分布，并分析温度分布随时间的变化直至温度分布稳定为止。

t2ta2 x式中a为导温系数，cm2/s；a2。 解：

模型转化为标准形式：

1t2t ax2初始条件为：

tx,020

边界条件为：

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.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时刻） hold on

plot(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秒时温度分布 x

sol(9,:,1)

经过8秒时的温度分布为： x/cm 0 0.5263 1.0526 1.5789 2.1053 t/C 120.0000 112.5520 105.1653 97.8994 90.8100 x/cm 3.6842 4.2105 4.7368 5.2632 5.7895 t/C 71.0714 65.1202 59.5200 54.2784 49.3930 x/cm 7.3684 7.8947 8.4211 8.9474 9.4737 t/C 36.7095 33.0419 29.5877 26.2982 23.1207 或者求第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

Word文档

2.6316 83.9477 6.3158 44.8518 10.0000 20.0000 3.1579 77.3562 6.8421 40.6338

`

不同时刻温度分布图

将上图的视角转至xt平面也得到本图，从本图可知当时间达到15s时平壁的温度分布已近稳定。

Word文档

`

某厚度为20cm钢板原温度为20C，现将其置于1000C的炉中加热，平壁导热系数为34.8W/mC，对流传热系数h174W/m2C，导温系数为a0.555105m2/s；试分析温度分布随时间的变化及钢板表面温度达到500C时所需的时间。

t2ta2 x 解：

模型转化为标准形式：

1t2t2 ax初始条件为：

tx,020

边界条件为：

t0,t0.1,，ht0.1,t 0（平壁中心坐标为0，绝热）xx

函数： 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)

Word文档

`

%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 on

plot(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 on

plot(x,sol(25,:,1)) plot(x,sol(31,:,1)) plot(x,sol(37,:,1))

600s的温度分布变化

Word文档

`

2160s的温度分布变化

不同时刻温度分布图

Word文档

`

a0.5551053601740.10.2 Bi0.5，Fo2234.80.1h根据以上两准数可知：该传热过程外对流及导热阻力相当；当加热时间小于360s时为非

正规阶段，加热时间大于360s后进入正规阶段，从上图也可得到该结论。

由该图可知，当加热时间达到2160s时，钢板的表面温度达到500C，钢板中心温度为： [uout,duoutdx]=pdeval(0,x,sol(37,:,:),[0]*1e-2) 371.2850C

有关传热量问题：（钢板的表面温度达到500C时，所需的总传热量）

QQ1Aexp12FoB ，Q0cVt0t00对平板：（参见传热学）

10.66818，A2sin1sin11.07329，B0.92723

1sin1cos11Aexp12FoB1.07329e0.53520.927230.5827 0QQ110.58270.41726 Q0cVt0t0Word文档

`

Q0cVt0t174810.12010006.1448610J 50.55510Q0.417266.144861082.564108J

MatLab解法： 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:36

q=q+174*((sol(i,11,1)+sol(i+1,11,1))/2-1000)*60;（所需的热量均是从表面传递进入钢板的） end q

结果：Q2.4797108J

Word文档

`

有一直径为40cm钢锭温度为20C，将其置于900C的炉中加热，平壁导热系数为

34.8W/mC，对流传热系数h174W/m2C，导温系数为a0.695105m2/s；试分

析温度分布随时间的变化及钢锭表面温度加热到750C时所需的时间；钢锭可近似为无限长的圆柱。

t1tax xxx初始条件为：

tx,020

边界条件为：

t0,t0.2,，ht0.2,t 0（平壁中心坐标为0，绝热）

xx函数：

pdefun2.m icbun2.m bcfun2.m

命令：

x=[0:1:20]*1e-2;

t=[0:120:5520]; (结合图5520s时，钢锭表面温度达到750C左右) sol=pdepe(1,pdefun2,icfun2,bcfun2,x,t); m=1(圆柱) mesh(x,t,sol)

Word文档

`

PDETOOL工具求解二维稳态与动态PDE：

如图偏心环形空间表面温度为100C，外表面温度为20C；试给出其温度分布。

如图导热物体，下表面温度为20C，上部三角形截面处温度为100C，其余各面绝热，试给出其温度分布。

Word文档

`

有一砖砌的烟气通道，其截面形状如附图所示，边长为1m，管道直径为0.5m。已知外壁温分别为80C、25C，砖的导热系数为1.5W/mC，试确定该通道的温度分布、距离任意相邻两直角边各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);

interp2(x,y,uxy,-0.4,-0.4) 结果：29.9010 interp2(x,y,uxy,-0.4, 0.4) 结果：29.9004 interp2(x,y,uxy, 0.4,-0.4) 结果：29.9015 interp2(x,y,uxy, 0.4, 0.4) 结果：29.9063 故其温度为29.9C

t n根据本题其四面对称，所以计算其中一面；传递的热量在表面处难以计算（圆形的表面），但传递的热量必然通过外表面，其传热量只需将表面处各节点的一阶导数与节点间的面积、导热系数相乘即可得到传热量。 命令：

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.0801

dudx=( 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

QAWord文档

`

184.6430 188.6017 184.9701 174.9577 160.6397 141.7009 120.5195 97.6419 72.9719 50.2694 26.1751 0 y方向一阶导数

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C，零时刻开始外壁温分别维持在为80C、25C，砖的导温系数为0.00001W/mC，试确定该通道温度分布随时间的变化及温度分布稳定所需的时间。

2t2ttax2y2  解：pdetool工具各菜单中的参数设置略，其中solve菜单parameters选项的时间设置为0:100:8000。

首先在稳态情况下（椭圆形）求解（三角形网格经过两次细化），导出解u；而后解以上非稳态（抛物形）并导出解u1（三角形网格细化与稳态相同）和三角形网格参数p，e，t；

命令：

Word文档

`

pdemesh(p,e,t,u1(:,1)) 0s pdemesh(p,e,t,u1(:,2)) 100s pdemesh(p,e,t,u1(:,4)) 300s pdemesh(p,e,t,u1(:,7)) 600s pdemesh(p,e,t,u1(:,11)) 1000s pdemesh(p,e,t,u1(:,21)) 2000s pdemesh(p,e,t,u1(:,31)) 3000s pdemesh(p,e,t,u1(:,51)) 5000s pdemesh(p,e,t,u1(:,81)) 8000s

稳定时间（以某时刻的温度分布与稳态温度分布的相对误差来判断） norm((u-u1(:,21))./u) 结果：2.4511 norm((u-u1(:,31))./u) 0.8079 norm((u-u1(:,51))./u) 0.0770 norm((u-u1(:,81))./u) 0.0047

可见在8000s以后的温度分布已基本没有变化。

Word文档

`

u1

u0

h=1 r=0 c=1 a=0 f=1 0.25

0.2

0.15 0.1

0.05

0 10.5

0

-0.5-0.5 -1-1

误差plot Height3D中选择user entry : u-(1-x.^2-y.^2)/4

-4 x 10 1.5Word文档

10.50-0.510.50 `

u1-x-y

u0

h=1 r=0 c=1 a=0 f=1-x-y

需采用solve parameters中use nonlinear solver方可收敛

0.08

0.06

0.04

0.02

0

-0.02 -0.04 -0.06 -0.08-1 -0.50-0.50.5 1-1

Word文档

100.5 `

Word文档