buondary condition derivative equal zero PDE
2 views (last 30 days)
Show older comments
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')
3 Comments
Accepted Answer
Torsten
on 22 Mar 2022
Edited: Torsten
on 23 Mar 2022
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 Comments
Torsten
on 24 Mar 2022
Edited: Torsten
on 24 Mar 2022
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,:))
More Answers (0)
See Also
Categories
Find more on General PDEs in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!