首页/文章/ 详情

偏微分方程实例程序

过冷水

3年前浏览2300

过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，qq：927550334

过冷水在往期和大家分享了偏微分方程的数值解法，并且给出案例，本期过冷和大家分享更多实用案例，瞬态热传导、牛顿流体方程、固定床二维均相模型。

瞬态热传导

针对这样一个问题：

边界条件：

初始条件：

已知：

原方程可化为标准形式：

据此求解程序如下：

function pde
global h k rho Cp Tr
h=0.17;k=0.017;
T0=323;Tr=293;
rho=4.92*10^4;Cp=0.4605;
R=0.12;m=1;
rmesh=linspace(0,R,30);
tspan=0:0.001:0.1;
sol = pdepe(m,@pdefun,@icfun,@bcfun,rmesh,tspan);
surf( rmesh,tspan,sol)
xlabel('r'), ylabel('time'), zlabel('Temperature');
shading interp
colormap cool

function [c,f,s]= pdefun(x,t,u,du)
    c=rho*Cp;
    f=k*du;
    s=0;
end

function u=icfun(x)
    T0=323;
    u=T0
end
function [pl,ql,pr,qr]= bcfun(xl,ul,xr,ur,t)
    pl = ul;
    ql =0;
    pr =h*( ur-Tr);
    qr=k;
end
end

结果为：

牛顿流体运动方程

水平圆管中装有静止的牛顿流体﹐从某个时刻起外加压力△p，流体开始流动，以采用如下方程描述:

边界条件：

初始条件：

已知：

程序：

function pde
global rho miu deltP
rho=1000;
miu=1;
L=0.5;
dP=3*10^5;
deltP=dP/L;
R=0.002;
m=1;
rmesh=linspace(0,R,30);
tspan=0:0.001:0.1;
sol = pdepe(m,@pdefun,@icfun,@bcfun,rmesh,tspan);
surf( rmesh,tspan,sol)
shading interp
colormap cool
pp=pchip(sol(:,1),tspan);
SF=0.98*sol(end,1);
ST=ppval(pp,SF);
function [c,f,s]= pdefun(x,t,u,du)
    c=rho;
    f=miu*du;
    s=deltP;
end

function u=icfun(x)
    u=0
end
function [pl,ql,pr,qr]= bcfun(xl,ul,xr,ur,t)
    pl = ul;
    ql =0;
    pr =ur;
    qr=0;
end
end

结果为：

固定床均相模型

在固定床反应器中进行乙烯催化氧化制备环氧乙烷反应,如果忽略环氧乙烷的燃烧反应,则反应器内进行的反应如下。

采用二维拟均相反应器模型进行模拟,求产物环氧乙烷浓度和反应器温度沿反应器管长的变化。反应器模型可构建为：

初始条件:

边界条件:

其中D_er和λ_er分别为流体径向有效扩散和传热系数,数值分别为D_er=6×10^-⁴ m² /s，λ_er= 0.1 J/(m²· s.K); h_w为反应器器壁的导热系数,其值为250 J/(m²· s ·K);

源程序为：

dt=0.0508;
Tw=498;
L=12;
rmesh=linspace(0,dt,10);
tspan=linspace(0,L,40);
m=1;
sol = pdepe(m,@pdefun,@icfun,@bcfun,rmesh,tspan);
%求不同床层位置的平均温度和产物浓度
CEO=sol(:,:,3);
T=sol(:,:,4);
[row,col]=size(CEO);
for i=1:row
cpp=pchip(rmesh,CEO(i,:));
cav(i)=quadl(@ppval,rmesh(1),dt,[],[],cpp)./dt;
tpp=pchip(rmesh,T(i,:));
Tav(i)=quad(@ppval,rmesh(1),dt,[],[],tpp)/dt;
end
plot(tspan,cav,'Linewidth',2);
xlabel('Bed Length [m]','fontsize',16)
ylabel('EO Concentration [mol/m^3]','fontsize',16)
set(gca,'Fontsize',16)
figure
plot( tspan,Tav,'linewidth',2);
xlabel('Bed Length [m]','fontsize',16);
ylabel ( 'Bed Temperature [K]','fontsize',16)
set(gca,'Fontsize',16)

function C=icfun(r)
C=[224,14,0,498];
end

function [c,f,s]= pdefun(r,z,C,dCdr)
dt=0.0508;
Tw=498;
CEH=C(1);CO2=C(2);CEO=C(3);T=C(4);
us= 1.3;
rhouf = 6.06;
cp= 1160 ;
Uw=270;
Der= 6e-4;
lbr=0.1;
dH1=-210000;
dH2=-473000;
R=8.314 ;
k1= 35.2 * exp( -59860/R/T);
k2= 24700 * exp(-89791/R/T);
R1=810 *k1*CO2;
R2=2430* k2*CO2;
c=[us/Der;us/Der;us/Der;us*rhouf*cp/lbr];
f=[dCdr(1); dCdr(2); dCdr(3); dCdr(4)];
s=[-(2*R1 1/3*R2)/Der;
-(R1 R2)/Der;
2*R1/Der;
((-dH1*R1-dH2*R2)-4*Uw*(T-Tw)/dt)/lbr];
end

function [pl,ql,pr,qr]= bcfun(rl,ul,rr,ur,z)
dt=0.0508;
Tw=498;
hw =250;
lbr=0.6;
pl =[0;0;0;0];
ql=[1;1;1;1];
pr =[0;0:0;-hw*(ur-Tw)/lbr];
%pr =[0;0;0;0];
qr=[1;1;1;1];
end