How to define a step function as input BC in pdepe?

1 view (last 30 days)
I'm trying to define a stepchange of inflow concentration(cin) in my system. The following is my code, in line 27 cin should be from 10 to 0 at t=50, so I defined cin = 10.*heaviside(50-t), but there is always an error :
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
Error in ADSmainsolver>slowsorpbc (line 53)
pl = [ul(1)-cin;0];
Error in pdepe (line 250)
[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});
Error in ADSmainsolver (line 30)
c = pdepe(0,@transpdes,@slowsorpic,@slowsorpbc,x,t,options,D,v,theta,rhob,F,kd,k,c0,cin);
Here's my code, could anyone help me,plzzzzzz :
%ADSmainsolver
clear all
clc
close all
% numerical PDE solution for transport (a+d) and nonideal sorption
%--------------------------------------------------------------------------
T = 156.5; % maximum time [min]
T1 = 1.5; % minimum time [min]
L = 15; % length [m]
D = 0.4; % diffusivity [cm^2/min]
v = 1; % real fluid velocity [cm/min]
theta = 0.4; % porosity
rhob = 1.6; % porous medium bulk density [g/cm^3]
c0f = 0; % initial concentration in fluid [ug/L]
c0s = 0; % initial concentration in solid [ug/kg]
c0=[c0f;c0s];
F = 1; %fraction of ideal partition
kd = 1; %partition coefficient [kg/L]
k = 0; % sorption rate limit const. [1/min]
M = 32; % number of timesteps (>2)
N = 100; % number of nodes
%-------------------------- output parameters------------------------------
gplot = 1; % =1: breakthrough curves; =2: profiles 1
t = linspace (T1,T,M); % time discretization
x = linspace (0,L,N); % space discretization
cin = 10.*heaviside(50-t);% inflow concentration [ug/L]
%----------------------solver-------------------------------------------
options = odeset;
c = pdepe(0,@transpdes,@slowsorpic,@slowsorpbc,x,t,options,D,v,theta,rhob,F,kd,k,c0,cin);
Y=c(:,N,1);
%---------------------- graphical output ----------------------------------
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] = transpdes(x,t,u,DuDx,D,v,theta,rhob,F,kd,k,c0,cin)
c = [1+rhob*F*kd/theta;1];
f = [D;0].*DuDx;
s = -[v;0].*DuDx - ([k*(1-F)*kd,-k]*u)*[rhob/theta;-1];
end
function u0 = slowsorpic(x,D,v,theta,rhob,F,kd,k,c0,cin)
u0 = c0;
end
% --------------------------------------------------------------------------
function [pl,ql,pr,qr] = slowsorpbc(xl,ul,xr,ur,t,D,v,theta,rhob,F,kd,k,c0,cin)
pl = [ul(1)-cin;0];
ql = [0;1];
pr = [0;0];
qr = [1;1];
end

Answers (2)

Rohit Pappu
Rohit Pappu on 2 Sep 2020
There’s an error in Line 27 because, cin is of size 1x32 and MATLAB cannot vertically concatanate it with a scalar.
To define a unit step function, without using Heaviside function :
unitstep = 10*ones(size(t));
unitstep(t>=50) = 0;

Bill Greene
Bill Greene on 3 Sep 2020
Here is the way you want to define such a BC:
if(t>=50)
pl = [ul(1);0];
else
pl = [ul(1)-10;0];
end
One thing you have to remember is that boundary condition function may be called for any value of time between the starting and ending times-- not just the times where you have requested output.

Categories

Find more on Fractals in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!