HI, I am trying to solve two Simulta PDEs for diffusion through packed bed using pdepe. I have understood the syntax but the code isn't getting solved. Can anyone help me correct the code? The problem might be with identifying the values of c,f,s.

1 次查看（过去 30 天）

simopera6 2015-9-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/245974-hi-i-am-trying-to-solve-two-simulta-pdes-for-diffusion-through-packed-bed-using-pdepe-i-have-unde

编辑： Torsten 2015-10-8

Equation 1: Dx*D2C/Dx2-v*DC/DX-Dq/Dt=Dc/Dt Equation 2: Dq/Dt= K1*[C*(qm-q)-Kd*q] Dx= diffusion coefficient v= velocity q= concentration in solid phase (function of c) qm= maxium value of colid phase concentration K1= adsorption constant kd= adsorption coeff. I have intial and boundary conditions which I can write.

Please see the following code.

function pdecolumn %defining variables from Biochemical engineering 49 221-228 global v x Dcoe Cin qm k1 kd t m L dc; m=0; v=1;%cm/min Dcoe= 1.6;% L =2.50; %cm dc=1.6; %cm Cin=3; %mg/ml qm=113.6; %mg/ml k1=0.035;%ml/mg.min kd=0.008;%mg/ml

x=linspace(0,L,101); t=linspace(0,12000,101);

sol=pdepe(m,@pdecolumnpe,@pdecolumnic,@pdecolumnbc,x,t); tend= 12000; % Extract the first solution component as conc. conc1 = sol(:,:,1); conc2=sol(:,:,1)./Cin; surf(x,t,conc2) % A solution profile for conc. vs length. figure, plot(x,conc1(end,:)) title(strcat('Solution at t = ', num2str(tend))) xlabel('Distance x') ylabel('Concentration (C)')

%Plot conc. vs. time figure, plot(t,conc1) title('conc vs time: breakthrough') xlabel('Time (s)') ylabel('Concentration (C)')

figure, plot(t,conc2(end,:)) title('conc vs time: breakthrough') xlabel('Time (s)') ylabel('Concentration (C)') end

%% % main file function [C,f,s] = pdecolumnpe(x,t,u,DuDx) %Defining variables global v Dcoe qm k1 kd;

C=[1;1]; f=[Dcoe.*DuDx(1)-(v.*u(1));0];

s1= -DuDt(2); s2= k1.*(u(1).*(qm-u(2)))-kd.*u(2); s=[s1;s2]; end

%% % Initial conditions function u0 = pdecolumnic(x) global Cin; u0=[Cin;0]; end

%% %Boundary conditions function [pl,ql,pr,qr] = pdecolumnbc(xl,ul,xr,ur,t) global Cin; pl = [ul(1)-Cin; ul(2)]; ql = [0;0]; pr = [0;0]; qr = [1;1] ; end %%%%%%