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.
11 ビュー (過去 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 で Ordinary Differential Equations についてさらに検索
製品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!