I'm try to define the stepchange of bc in pdepe which pl varies with time. The following is my code.
Can anyone help? plzzzz…
T = 16; % maximum time [s]
L = 8; % length [m]
D = 0.1; % diffusivity [m*m/s]
v = 0.5; % real fluid velocity [m/s]
theta = 0.2; % porosity
rhob = 1200; % porous medium bulk density [kg/m*m*m]
kappaf = 0.00; % transition rate fluid to solid [1/s]
kappas =0.1; % transition rate solid to fluid [1/s]
lambdaf = 0; % decay rate in fluid [1/s]
lambdas = 0; % decay rate in solid [1/s]
c0f = 1; % initial concentration in fluid [kg/m*m*m]
c0s = 0.1; % initial concentration in solid [1]
M = 100; % number of timesteps (>2)
N = 40; % number of nodes
%-------------------------- output parameters------------------------------
gplot = 1; % =1: breakthrough curves; =2: profiles 1
t = linspace (T/M,T,M); % time discretization
x = linspace (0,L,N); % space discretization
cin = 10; % inflow concentration [kg/m*m*m]
bctimes=[0,5,5.1,T];bcVals=[cin,cin,0,0];%----------------------execution-------------------------------------------
bcFunc=@(xl,ul,xr,ur,t) slowsorpbc(xl,ul,xr,ur,t,bctimes,bcVals);options = odeset;c = pdepe(0,@slowsorpde,@slowsorpic,bcFunc,x,t,options,D,v,theta,rhob,kappaf,kappas,lambdaf,lambdas,[c0f;c0s]);%---------------------- graphical output ----------------------------------
Y=c(:,N,1);switch gplot case 1 plot (t,c(:,N,1)) % breakthrough curves
xlabel ('time'); ylabel ('concentration'); case 2 plot (x,c(:,:,2)','--') % profiles
xlabel ('space'); ylabel ('concentration');end% --------------------------------------------------------------------------
function [c,f,s] = slowsorpde(x,t,u,DuDx,D,v,theta,rhob,kappaf,kappas,lambdaf,lambdas,c0)c = [1;1];f = [D;0].*DuDx;s = -[v;0].*DuDx - [lambdaf;lambdas].*u - ([kappaf,-kappas]*u)*[1/theta;-1/rhob];end% --------------------------------------------------------------------------function u0 = slowsorpic(x,D,v,theta,rhob,kappaf,kappas,lambdaf,lambdas,c0)u0 = c0;end % --------------------------------------------------------------------------function [pl,ql,pr,qr] = slowsorpbc(xl,ul,xr,ur,t,bctimes,bcVals)pl = [ul(1)-interp1(bctimes,bcVals,t);0];ql = [0;1];pr = [0;0];qr = [1;1];end
Best Answer