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

3 次查看（过去 30 天）
Yi Peng Biao
Yi Peng Biao 2020-8-17
回答： Bill Greene 2020-9-3
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

请先登录，再回答此问题。

回答（2 个）

Rohit Pappu
Rohit Pappu 2020-9-2
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 2020-9-3
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.

类别

在 Help Center 和 File Exchange 中查找有关 Just for fun 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by