clc
clear all
dx=0.01; %x方向的步长
dt=1; %t方向的步长
r1=dt/(dx^2);r2=dt/dx;%计算r的值
x=0:dx:1.5; %得到x的序列
t=0:dt:60; %得到t的序列
M=length(x)-1;
N=length(t)-1;
%相关数据
cp=4178;p=1000;WH=0.5;m=0.05;uin=45;%a=0.015;v=1;
u0=20;u1=37;
u=ones(N+1,M+1);
u=u*u1;
% q=(u-u0)*(4.692+1.934+3.804)/WH;
% Phi(1,:)=100; %设置初值条件：Phi(x,0)=100
% Phi(2:N+1,1)=0; %设置边界条件：Phi(0,t)=0
% Phi(2:N+1,M+1)=0; %设置边界条件：Phi(1,t)=0
%j=1;i=1;
%根据差分方程，计算Phi的数值解
a=3e-7;v=5e-5;
u(:,1)=uin;u11=u;u22=u;
for j=1:N
for i=1:M
q=(u11-u0)*(4.692+1.934+3.804)/WH;
u11(j+1,i+1)=(1-2*r1*a)*u11(j,i+1)+(a*r1+v*r2/2)*u11(j,i)+(a*r1-v*r2/2)+q(j,i+1)*dt/(cp*p);
% Phi(j+1,i)=(1-2*r)*Phi(j,i)+r*(Phi(j,i-1)+Phi(j,i+1));
end
end
for j=1:N
for i=1:M
u22(j+1,i+1)=u22(j,i)-m*dt*(uin-u22(j,i))/(p*WH*dx);
end
end
u=u11-u+u22;
for j=1:N
u(j,M+1)=2*u(j,M)-u(j,M-1);
end
% q(j,i)=(u(j,i)-u0)*(4.692+1.934+3.804)/WH;
%flag=(min(min(u)')>=37&&max(max(u)')<=60&&max(max(u)')>=44);
u=u(2:N+1,2:M+1);
%image(u)
pcolor(u)
colorbar;
% [x,t]=meshgrid(x,t);
% mesh(u) %绘制(x,t,Phi)的三维图
xlabel('x')
ylabel('t')
%zlabel('\u(x,t)')
title('扩散方程的数值模拟')
评论0