Sto provando a risolvere un problema di carbocementazione usando ODE45, in particolr modo sto usando il metodo delle linee. Lo scopo dell'esercizio è quello di ricavare dopo un determinato istante di tempo T il grado di assorbimento del carbonio. Non riesco proprio a capire come risolvere perchè mi sembra che sia inserito tutto correttamente. qualcuno mi può aiutare?
clear
clc
%cabocemntazione
%defnizione parametri
L=1; Temp=1273; D0=20*10^-6; Q=142000; R=8.314; T=10;
nx=100; nt=100;
%controllo parametri
if any([L Temp T D0 nx-2 nt-2])<=0
error('controlla i parametri')
end
%inizializzazione
D=D0*exp(-Q/(R*Temp)); %diffusività
dx=L/nx; %ho discretizzato dx perchè mi sefve nella formula di p
p=D/dx^2;
x=linspace(0,L,nx+1);%non ho fatto il tempo non mi serve
A=[p*ones(nx,1) -2*p*ones(nx,1) p*ones(nx,1)];
AA=spdiags(A, -1:1, nx,nx);
A(end,end-1)=2*p;
a=sparse(nx,1);
%condizione iniziale
c=GetPhi(nx);
[t,y]=ode45(@system,[0 T], c(2:end),[],A,a,p);
c(2:end)=y(end,:);
c(1)=Getg1;
plot(x,c,'b*');
axis([0 L 0 Getg1]);
xlabel('x'); ylabel('%w di C');
legend('MOL');
grid
pause(0.1)
%---------subfunction----------
function phi=GetPhi(nx)
c0=0.2;
phi=c0*ones(nx+1,1);
end
function g1=Getg1
cs=1;
g1=cs;
end
function g2=Getg2
g2=0;
end
function cdot=system(t,c,A,a,p);
a(1)=p*Getg1;
a(end)=p*Getg2;
cdot=A.*c+a ;
end
clear
clc
%cabocemntazione
%defnizione parametri
L=1; Temp=1273; D0=20*10^-6; Q=142000; R=8.314; T=10;
nx=100; nt=100;
%controllo parametri
if any([L Temp T D0 nx-2 nt-2])<=0
error('controlla i parametri')
end
%inizializzazione
D=D0*exp(-Q/(R*Temp)); %diffusività
dx=L/nx; %ho discretizzato dx perchè mi sefve nella formula di p
p=D/dx^2;
x=linspace(0,L,nx+1);%non ho fatto il tempo non mi serve
A=[p*ones(nx,1) -2*p*ones(nx,1) p*ones(nx,1)];
AA=spdiags(A, -1:1, nx,nx);
A(end,end-1)=2*p;
a=sparse(nx,1);
%condizione iniziale
c=GetPhi(nx);
[t,y]=ode45(@system,[0 T], c(2:end),[],A,a,p);
c(2:end)=y(end,:);
c(1)=Getg1;
plot(x,c,'b*');
axis([0 L 0 Getg1]);
xlabel('x'); ylabel('%w di C');
legend('MOL');
grid
pause(0.1)
%---------subfunction----------
function phi=GetPhi(nx)
c0=0.2;
phi=c0*ones(nx+1,1);
end
function g1=Getg1
cs=1;
g1=cs;
end
function g2=Getg2
g2=0;
end
function cdot=system(t,c,A,a,p);
a(1)=p*Getg1;
a(end)=p*Getg2;
cdot=A.*c+[p*Getg1; a; p*Getg2] ;
end
