buondary condition derivative equal zero PDE

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]

In which coordinate system did you set up your PDE ?

If it's in cylindrical or spherical coordinates (r direction), your substitution of the conduction term is wrong.

[Edoardo Bertolotti]

yes, cylindrical

[Edoardo Bertolotti]

rho=8900; %[kg/m^3] %density

r=0.05; %[m] %radius

c=15; %[J/m*s*C] %conducibility

Cp=600; %[m] %specific heat

diffus= c/(rho*Cp); %[cm^2/s] %diffusivity

a=(1/0.001)*(R/diffus)+(1/h^2)

Answer 1

[Torsten]

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))]); 

[Edoardo Bertolotti]

Hello,

honestly I cannot understand the loop, but graph-wise, seems to make sense

mesh(r,t,u')

title('Temperatures: explicit method','interpreter','latex')

xlabel('r [m]')

ylabel('time [sec]','interpreter','latex')

zlabel('Temperature','interpreter','latex')

I did a mistake in the data: the radius is 0.1 meters. And the problem ask when I reach 873.15K at 0.005 meters.

So by setting this new radius, I get

T = 873.23K at r=0.0034 at t=605.3270 (coordinates 2,287 of matrix "u")

T = 872.4383K at r=0.0069 at t=603.2105 (coordinates 3,286 of matrix "u")

Does it make sense?

[Torsten]

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,:))