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.
7 次查看(过去 30 天)
显示 更早的评论
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 %%%%%%
1 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!