Axial Dispersion Model, How to solve the following set of equations

here dq/dt = k*(q*-q)

q* = HC

H= henry constant

Initial conditions

at t=0

z=0 to l

q=0

c=0

boundary condition

at t>0

z=0

c=2

How can i solve the following equation, i am stuck with the above problem

Answer 1

[Torsten]

Use MATLAB's "pdepe".

As boundary condition for q, set

pl = 0, ql = 1, pr = 0, qr = 1.

[Yash Toshniwal]

L=10; %Width of plate
dx=0.01; % Grid Space
x=0:dx:10; %Space vector
dt=0.003; %Time step
t_final=30;
t=0:dt:t_final; % Time vector
n=floor(L/dx)+1; %Number of spacial nodes
m=floor(t_final/dt)+1; %Number of time nodes
k=0.1;
H=2;
ph=2.33;
u=1.061; %Convection velocity
c=(u*dt)/(dx); %Convection Number
alpha=0.006; %Thermal diffusivity
d=(alpha*dt)/(dx^2); %Diffusion Coeff
Pe=(u*L/alpha); %Peclet Number
T=zeros(m,n);
nod=numel(x);
T(1,:)=0; %At t=0
qstar=H.*T;
q=qstar.*T;
T(:,1)=2;%1st surface
qstar=H.*T;
dq=k*(qstar-q);
for j=1:m-1
   for  i=2:n-1 
       T(j+1,i)=T(j,i)-(c/2)*(T(j,i+1)-T(j,i-1))+d*(T(j,i+1)-2*T(j,i)+T(j,i-1))-(dq(j,i).*ph);
        qstar(j+1,i)=H*T(j,i);
        q(j+1,i)=qstar(j,i)*T(j,i);
        dq(j+1,i)=k*(qstar(j,i)-q(j,i));
   end
end
C1=T(500,:); %C at 2 min
C2=T(1000,:); %C at 6 min
C3=T(2000,:); %C at 10 min
C4=T(2500,:); %C at 30 min
plot(x,C1,'r');
hold on
plot(x,C2,'b');
hold on
plot(x,C3,'g');
hold on
plot(x,C4,'m');
hold off
xlabel('x(cm)');
ylabel('Concentration(C)');
title('Concentration profile ');
legend('c=2min','c=6min','c=10min');
xlim([ 0 10]);
ylim([0 2]);

The following is my matlab code, can you please suggest me where i am going wrong?

[Yash Toshniwal]

function pdex4

m = 2;

x = [0:0.01:10];

t = [0:0.003:30];

solution = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);

c1=solution;

figure(x,t,c1);

end

function [u,f,s,p] = pdex4pde(x,t,c,DcDx)

H=2;

k=0.1;

u = 1;

C=2;

f = 0.006.* DcDx;

s = -1.061;

q=0;

p= k*((H*C)-q);

end

function [c0,q] = pdex4ic(x)

c0=0;

H=2;

c=2;

q=H*c*c0;

end

function [pl,ql,pr,qr] = pdex4bc(xl,cl,xr,cr,t)

pl = cl;

ql = 2;

pr = 0;

qr = 0;

end

%------------------------------------

This is how i implemented the pdepe function, as i have a minimal idea about how to proceed using the solver

[Torsten]

m = 2;

Why ? Do you work in a sphere ?

function [u,f,s,p] = pdex4pde(x,t,c,DcDx)

Why 4 return parameters ? pdepe needs 3.

u = 1;

f = 0.006.* DcDx;

s = -1.061;

u, f and s must be vectors of length 2 since you have equations for C and q.

function [c0,q] = pdex4ic(x)

Why 2 return parameters ? pdepe needs one array of length 2.

pl = cl;

ql = 2;

pr = 0;

qr = 0;

These settings don't make sense at all.

My advice:

Read the documentation of "pdepe".

[Hasti]

Hi Torsten,

In another question here

https://se.mathworks.com/matlabcentral/answers/371313-error-in-solving-system-of-two-reaction-diffusion-equations

I see that you commented ""pdepe" is not designed to solve mixtures of partial differential equations and ordinary differential equations (the equation for u is a partial differential equation, the equation for v is an ordinary differential equation)."

In my idea these two cases are similar, so I am wondering if such a system can be solved with pdepe?

Appreciate your reply.