when i use this code without reaction terms ( * react(j,i)*)the curve and result is reasonable(fig 1) but with applying reaction terms(fig 2) it goes wrong. i can not find where is my mistake. i would be grateful if anybody could help me.
problem at solving coupled PDE-ODE with reaction between 4 species
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
I have solved a time dependent coupled PDE-ODE system for 4 compound where there are some rection between them. my probem is that when i add reaction terms in my model, i got wrong result.the first coumpound should decrease and three others should increase but first one remains fixed and others decrease. that is very strange for me. i am sure i wrote correct expression for reactions but the result is reverse. this is my code:
clc
clear
format long
tf=60*1*60;
nt=100;
dt=tf/nt;
tspan = 0:dt:tf;
nx = 50;
%%%%%%%parameters%%%%
Dab = 11.3*(1e-6);%
ux=0.5;
L=2e-3*(5/95);
Cin=[ 286.4879091 3.508114893 3.769381628 9.247436195]*1e-3*2.95;
yout=[233.9457273 5.289110609 26.08221981 31.552135]*1e-3*2.95;
K=[11.55507758 1.841304056 561.4994587 60.69688482 1.959738497];
Ncmp=4;
for i=1:8*nx
y0(1,i)=0;
end
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'InitialStep', 0.0001);
[t,y] = ode15s(@fun,tspan,y0,options,nx,Ncmp,L,Dab,Cin,ux,K);
[s1,s2] = size(t);
plot(t/60,y(1:s1,2*nx-1))
hold on
plot(t/60,y(1:s1,4*nx-1))
hold on
plot(t/60,y(1:s1,6*nx-1))
hold on
plot(t/60,y(1:s1,8*nx-1))
hold on
plot(25,yout,'r*')
function dydt=fun(~,y,nx,Ncmp,L,Dab,Cin,ux,K)
kg=0.66;
area=4;
eps=0.65;
rho=4.23*1e6;
B=rho;
E=1; %<---E can be changed
dx=L/(nx+1);
dydt=zeros(8*nx,1);
k1=K(1);
k1a=K(2);
k2=K(3);
k3=K(4);
k4=K(5);
for j=1:Ncmp
for i=1:2*nx
C(j,i)=y(i+2*nx*(j-1)); % concentration of j in node i,
end
end
%%%%%%%%%%%Reactions%%%%%%%%
for i=2:2:2*nx
rd1(i)=k1*k1a*C(1,i)/(1+k1a*C(1,i));
rd2(i)=k2*C(2,i);
rd3(i)=k3*C(3,i);
rd4(i)=k4*C(4,i);
react(1,i)=-rd1(i);
react(2,i)=rd1(i)-rd2(i);
react(3,i)=rd2(i)-rd3(i);
react(4,i)=rd3(i)-rd4(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
for j=1:Ncmp
i=1;
C0(j)=Cin(j);
dcdx(j,i)=(C(j,i+2)-C0(j))/2/dx;
d2cdx2(j,i)=(C(j,i+2)-2*C(j,i)+C0(j))/dx/dx;
dcdT(j,i)= Dab*d2cdx2(j,i)-ux*dcdx(j,i)-kg*area*((1-eps)/eps)*(B*C(j,i)-C(j,i+1));
dcdT(j,i+1)=kg*area*(B*C(j,i)-C(j,i+1))+E*react(j,i+1); %<---E can be changed
for i = 3:2:2*nx-3
dcdx(j,i)=(C(j,i+2)-C(j,i-2))/2/dx;
d2cdx2(j,i)=(C(j,i+2)-2*C(j,i)+C(j,i-2))/dx/dx;
dcdT(j,i)= Dab*d2cdx2(j,i)-ux*dcdx(j,i)-kg*area*((1-eps)/eps)*(B*C(j,i)-C(j,i+1));
dcdT(j,i+1)=kg*area*(B*C(j,i)-C(j,i+1))+E*react(j,i+1); %<---E can be changed
end
i = 2*nx-1;
C(j,i+2)=C(j,i)-(C(j,i-2)-C(j,i))/3;
dcdx(j,i)=(C(j,i+2)-C(j,i-2))/2/dx;
d2cdx2(j,i)=(C(j,i+2)-2*C(j,i)+C(j,i-2))/dx/dx;
dcdT(j,i)= Dab*d2cdx2(j,i)-ux*dcdx(j,i)-kg*area*((1-eps)/eps)*(B*C(j,i)-C(j,i+1));
dcdT(j,i+1)=kg*area*(B*C(j,i)-C(j,i+1))+E*react(j,i+1); %<---E can be changed
end
for j=1:Ncmp
for i = 1:2*nx
dydt(i+2*nx*(j-1)) = dcdT(j,i);
end
end
end
10 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!