1 view (last 30 days)

Show older comments

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!