Pdepe matlab convection diffusion

10 views (last 30 days)
Tasneem Na
Tasneem Na on 27 Jul 2015
Commented: Torsten on 28 Jul 2015
I want to solve this pde with initial and boundary conditions given. Tried Matlab's pdepe, but does not work satisfactorily. Maybe the boundary conditions is creating problem for me. I also used this isotherm equation for equilibrium: q = K*Cp^(1/n). This is convection-diffusion equation but i could not find any write ups that addresses solving this type of equation properly.
sol = pdepe(m,@ParticleDiffusionpde,@ParticleDiffusionic,@ParticleDiffusionbc,x,t);
% Extract the first solution component as u.
u = sol(:,:,:);
function [c,f,s] = ParticleDiffusionpde(x,t,u,DuDx)
global Ds
c = 1/Ds;
f = DuDx;
s = 0;
function u0 = ParticleDiffusionic(x)
global qo
u0 = qo;
function [pl,ql,pr,qr] = ParticleDiffusionbc(xl,ul,xr,ur,t,x)
global Ds K n
global Amo Gc kf rhop
global uavg
global dr R nr
sum = 0;
for i = 1:1:nr-1
r1 = (i-1)*dr; % radius at i
r2 = i * dr; % radius at i+1
r1 = double(r1); % convert to double precision
r2 = double(r2);
sum = sum + (dr / 2 * (r1*ul+ r2*ur));
end;
uavg = 3/R^3 * sum;
ql = 1;
pl = 0;
qr = 1;
pr = -((kf/(Ds.*rhop)).*(Amo - Gc.*uavg - ((double(ur/K)).^2).^(n/2) ));`
dq(r,t)/dt = Ds( d2q(r,t)/dr2 + (2/r)*dq(r,t)/dr )
q(r, t=0) = 0
dq(r=0, t)/dr = 0
dq(r=dp/2, t)/dr = (kf/Ds*rhop) [C(t) - Cp(at r = dp/2)]
q = solid phase concentration of trace compound in a particle with radius dp/2
C = bulk liquid concentration of trace compound
Cp = trace compound concentration at particle surface
  2 Comments
Walter Roberson
Walter Roberson on 27 Jul 2015
I moved the text above the code to make the question easier to read.
Torsten
Torsten on 28 Jul 2015
It's not obvious for me how
pr = -((kf/(Ds.*rhop)).*(Amo - Gc.*uavg - ((double(ur/K)).^2).^(n/2) ));
corresponds with
dq(r=dp/2, t)/dr = (kf/Ds*rhop) [C(t) - Cp(at r = dp/2)]
Further, I don't understand what you try to average in the for-loop. You refer to radii r1 and r2, but you use the values of u at the boundary. Moreover, averaging would be choosing r^2*u instead of r*u in spherical coordinates.
Best wishes
Torsten.

Sign in to comment.

Answers (0)

Products

Community Treasure Hunt

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

Start Hunting!