# problem at solving coupled PDE-ODE with reaction between 4 species

1 view (last 30 days)
Moji on 9 Nov 2018
Edited: Moji on 9 Nov 2018
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
Moji on 9 Nov 2018
in this line :
dcdT(j,i+1)=kg*area*(B*C(j,i)-C(j,i+1))+E*react(j,i+1);
i changed react(j,i+1) to react(j,i). and also i changed the for loop in below lines. i used for i=1:2:2*nx instead of for i=2:2:2*nx:
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);