buondary condition derivative equal zero PDE
2 次查看(过去 30 天)
显示 更早的评论
Hello, I am trying to solve a P.D.E. problem with explicit and implicit method (finite difference methods)
So I am building a grid for the Temperature profile through space and time.
How can I set that the derivative is zero at radius = 0?
*there was an error in the previous code, coordinates were wrong, like it was pointed out correctly in the comments
clear all
clc
rho=8900; %[kg/m^3] density
c=15; %[J/m*s*C] conducibility
Cp=600; %[m] specific heat
diffus= c/(rho*Cp); %[m^2/s] diffusivity
R=0.1; %[m] %radius
t_start = 0.0;
t_final = 20;
time_steps = 1000;
space_steps = 30;
r = (linspace(0.0,R,space_steps)).';
dr = r(2)-r(1); % vs dr = R/space_steps; %[m]
dt= 0.5*dr^2/diffus; %vs dt = t/time_steps; %[sec]
A=diffus*dt/dr^2;
x=linspace(0,R,space_steps); %discretization
LL=length(x);
time=linspace(t_start,t_final,time_steps);%discretization
TT=length(time);
% T start (T=1000C) u=temperature
for i = 1:LL+1
x(i) =(i-1)*dx;
u(i,1) = 1000 + 273.15; %Kelvin
end
% T boundary (r=0 T=dT/dr=0 - r=R 25C)
for k=1:TT+1
u(1,k) = ????
u(LL+1,k) = 25+ 273.15; %Kelvin
time(k) = (k-1)*dt;
end
% Explicit method
for k=1:TT % Time
for i=2:LL % Space
u(i,k+1) =u(i,k) + 0.5*A*(u(i-1,k)+u(i+1,k)-2.*u(i,k));
end
end
mesh(x,time,u')
title('Temperatures: explicit method','interpreter','latex')
xlabel('r [m]')
ylabel('time [sec]','interpreter','latex')
zlabel('Temperature','interpreter','latex')
采纳的回答
Torsten
2022-3-22
编辑:Torsten
2022-3-23
rho = 8900;
cp = 600;
D = 15;
a = D/(rho*cp);
R = 0.05;
uR = 25 + 273.15;
u0 = 1000 + 273.15;
rstart = 0.0;
rend = R;
nr = 30;
r = (linspace(rstart,rend,nr)).';
dr = r(2)-r(1);
tstart = 0.0;
tend = 1000.0;
dt = 0.5*dr^2/a;
t = tstart:dt:tend;
nt = numel(t);
A = a * dt/dr^2;
u = zeros(nr,nt);
u(:,1) = u0;
u(nr,:) = uR;
for it = 1:nt-1
u(2:nr-1,it+1) = u(2:nr-1,it) + A*( ...
(1 + dr./(2*r(2:nr-1))).*u(3:nr,it) - ...
2*u(2:nr-1,it) + ...
(1 - dr./(2*r(2:nr-1))).*u(1:nr-2,it));
u(1,it+1) = u(2,it+1);
end
plot(r,[u(:,1),u(:,round(numel(t)/20)),u(:,round(numel(t)/10)) ,u(:,round(numel(t)/5)),u(:,numel(t))]);
2 个评论
Torsten
2022-3-24
编辑:Torsten
2022-3-24
rho = 8900;
cp = 600;
D = 15;
a = D/(rho*cp);
R = 0.1;
uR = 25 + 273.15;
u0 = 1000 + 273.15;
rstart = 0.0;
rend = R;
nr = 241;
r = (linspace(rstart,rend,nr)).';
dr = r(2)-r(1);
tstart = 0.0;
tend = 1000.0;
dt = 0.5*dr^2/a;
t = tstart:dt:tend;
nt = numel(t);
A = a * dt/dr^2;
u = zeros(nr,nt);
u(:,1) = u0;
u(nr,:) = uR;
for it = 1:nt-1
u(2:nr-1,it+1) = u(2:nr-1,it) + A*( ...
(1 + dr./(2*r(2:nr-1))).*u(3:nr,it) - ...
2*u(2:nr-1,it) + ...
(1 - dr./(2*r(2:nr-1))).*u(1:nr-2,it));
u(1,it+1) = u(2,it+1);
end
figure(1)
plot(r,[u(:,1),u(:,round(numel(t)/20)),u(:,round(numel(t)/10)) ,u(:,round(numel(t)/5)),u(:,numel(t))]);
% Post processing
r_query = 0.005;
ir_query = r_query/dr + 1; % nr above was adjusted so that r_query is a grid point (i.e. r_query/dr is integer)
% Of course, 2d interpolation with interp2 is also an
% option
T_query = 873.15;
t_query = interp1(u(ir_query,:),t,T_query);
t_query
figure(2)
plot(t,u(ir_query,:))
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 General PDEs 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!