buondary condition derivative equal zero PDE
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Edoardo Bertolotti
le 21 Mar 2022
Modifié(e) : Torsten
le 24 Mar 2022
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 commentaires
Réponse acceptée
Torsten
le 22 Mar 2022
Modifié(e) : Torsten
le 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 commentaires
Torsten
le 24 Mar 2022
Modifié(e) : Torsten
le 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,:))
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Geometry and Mesh dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!