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

Yash Toshniwal

### Yash Toshniwal (view profile)

さんによって質問されました 2019 年 6 月 26 日

### Torsten (view profile)

さんによって コメントされました 2019 年 6 月 26 日
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

#### 0 件のコメント

サインイン to comment.

## 1 件の回答

### Torsten (view profile)

2019 年 6 月 26 日

Use MATLAB's "pdepe".
As boundary condition for q, set
pl = 0, ql = 1, pr = 0, qr = 1.

Yash Toshniwal

### Yash Toshniwal (view profile)

2019 年 6 月 26 日
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

### Yash Toshniwal (view profile)

2019 年 6 月 26 日
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

### Torsten (view profile)

2019 年 6 月 26 日
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.