MATLAB Answers

Mesh problem with PDEPE

6 views (last 30 days)
Gianpio Caringella
Gianpio Caringella on 8 May 2020
Hi guys,
I'm using pdepe in order to solve the system of PDEs associated to a chromatographic colum considering the Lump Kinetic Model for the mass transfer between the two phases along the columwn.
The code I wrote works only when the grid point are a multiple of ten, like 10,100,1000 while in the other cases I have the error : "Error using pdepe (line 293)Spatialdiscretization has failed.Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative."
I report here the equations I used
And the equation for q* is :
Here is the code,I think that maybe there are some issues with the boundary condtions related to the behaviour of q ( in the code u(2)) :
clear all
clc
%Coefficients
D=1e-5;
v=0.1; %That's u in the equations above
k=100;
eps=0.4;
H=0.85;
L=1;
t_in=1;
q_sat=1;
t_start=0;
t_end=10;
ic=0;
c_in=0.5;
N_x=50; % space points
N_t=100; % time points
t_span = linspace(t_start,t_end,N_t);
x_mesh = linspace(0,L,N_x);
% Solver
pdefun = @(x,t,u,dudx) (L_isotherm(x,t,u,dudx,D,v,q_sat,k,H,eps));
pdebc = @(xl,ul,xr,ur,t) (boundary(xl,ul,xr,ur,t,v,c_in,t_in));
pdeinc=@(x) initial(x,ic);
m=0; % Problemy simmetry / 0 per una "slab"
sol = pdepe(m,pdefun,pdeinc,pdebc,x_mesh,t_span); %
%Printing the results
figure()
for i=1:length(t_span)
plot(x_mesh,sol(i,:,1));
pause(1)
end
%Langmuir isotherm u(1)=c(x,t) u(2)= q(x,t) c = concentration liquid
%phase, q concentration in solid phase
function [c,f,s] = L_isotherm(x,t,u,dudx,D,v,q_sat,k,H,eps)
%
c = [1 ; 1];
f = [D*dudx(1) - v * u(1) ; 0 ]; %flux term
s = [-k / eps * ( q_sat * H * u(1) / ( 1 + H * u(1) ) - u(2) ) ; k /(1- eps )* ( q_sat * H * u(1) / ( 1 + H * u(1) ) - u(2))];
end
%Initial conditions
function u0 = initial(x,ic)
u0=[ic;ic];
end
%Boundary conditions for a rectangular pulse
function [pl,ql,pr,qr] = boundary(xl,ul,xr,ur,t,v,c_in,t_in)
pl=[ - (c_in) * ( t>0 & t<t_in); ur(2) ];
ql=[-1/v ; 0];
qr=[1/v ; 0];
pr=[ ur(1); ur(2)];
end

  2 Comments

Bill Greene
Bill Greene on 9 May 2020
For your left BC you have:
pl=[ - (c_in) * ( t>0 & t<t_in); ur(2) ];
Did you intend?
pl=[ - (c_in) * ( t>0 & t<t_in); ul(2) ];
Gianpio Caringella
Gianpio Caringella on 9 May 2020
Hi,
yeah I've noticed that error this morning and know the things work better.
I was wondering if, instead of putting ul(2) = ur(2)=0 at boundaries , the analytical solution of the equation dq/dt=k/(1-eps)*(q*-q) can be used since in the points it becomes an ODE.

Sign in to comment.

Answers (0)

Tags

Products


Release

R2019a