Can we solve 3 pde simultaneously using pdepe with respect to space and time??
2 views (last 30 days)
Show older comments
I have a system of 3 equations with 3 variables and corresponding boundary conditions. I am using pdepe to solve these. Is it possible? Can anyone help me correcting the code I have written? Sharing the equation in an image attached. 3rd equation is : dQp/dt=Keq*Cp*(lambda-(sigma+charge)*Qp)^charge-Qp*Csin^charge
%% Main function function sma_porediffusion %% globalising variables and constants global v x F Ac Dcoe ap Cin charge lambda epsilon k_eq sigma t m L dc N Run_time CV Csin; m=0;
%% specifying system parameters ap=.0023; %particle size epsilon_e= 0.7; %porosity epsilon_p=0.8; % L =5.2; %cm dc= 0.7; %cm F=0.5;%ml/min Ac=3.142*(dc/2)^2; %cm2 v=F/Ac;%cm/min Dcoe=4.43*10^-3;%cm^2/min charge=2.14; lambda= 1833; % capacity factor in mM Cin=0.0085; %mg/ml k_eq= 0.12; % equilibrium coefficient = kads/kdes sigma= 49; % steric hinderence factor N=50; % meshing intervals Run_time =1; % min Csin= 1; %mM mobile phase salt concentration % Csout= 400;%mM end salt concentration % dq/dt=keq*C*(lambda-(sigma+charge)*q)^charge-q*(Csalt^charge) % Area= pi*(dc^2)/4; % area of cross section cm^2
%% x=linspace(0,L,N); t=linspace(0,Run_time,N); sol=pdepe(m,@sma_pdecolumnpe,@sma_pdecolumnic,@sma_pdecolumnbc,x,t);
%% Extract the first solution component as conc. conc1 = sol(:,:,1); conc2=conc1./Cin; %surf(x,t,conc2)
%% A solution profile for conc. vs length. % figure, plot(x,conc1(end,:)) % title(strcat('Solution at t = ', num2str(t))) % xlabel('Distance x') % ylabel('Concentration (C)') % % %Plot conc. vs. time % figure, plot(t,conc1) % title('conc vs time: breakthrough conc1') % xlabel('Time (min)') % ylabel('Concentration (mg/ml)')
figure, plot(t,conc1(:,end)) title('conc vs time: breakthrough') xlabel('Time (min)') ylabel('Concentration (mg/ml)')
%% Inlet conditions function u0 = sma_pdecolumnic(x) if x==0 u0=[Cin;0;0]; else u0=[0;0;0]; end end %% describing equations function [c,f,s] = sma_pdecolumnpe(x,t,u,DuDx) %Defining variables %global v Cin charge lambda epsilon k_eq sigma k_a k_d Csin; c=[1;1;1]; f=[Dcoe;0;0].*DuDx; %dq/dt=adsorption coefficient*c(ionic capacity-(steric %factor+characteristic charge)*q)^ char charge - desorption coeff*q*Csaltin^ %char charge
sq=k_eq.*u(1).*(lambda-(sigma+charge).*u(2))^charge- u(2)*Csin^charge; %disp(sq); s=[-v*DuDx(1)-(1-epsilon_e)*k_eq*ap*(u(1)-u(2));(1-epsilon_p)*sq+k_eq*ap*(u(1)-u(2));sq]; end %% boundary conditions function [pl,ql,pr,qr] = sma_pdecolumnbc(xl,ul,xr,ur,t) %global Cin; pl = [ul(1)-Cin;0;0]; ql = [0;0;0]; pr = [0;0;0]; qr = [1;1;1]; end end
Thanks for your inputs
0 Comments
Answers (0)
See Also
Categories
Find more on Partial Differential Equation Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!