Code covered by the BSD License

# Environmental Modeling

### Ekkehard Holzbecher (view profile)

20 Jul 2007 (Updated )

Introduction to Modelling and Simulation for Environmental Problems

slowsorp
```function slowsorp
% 1D transport - modelling with extensions for decay and slow sorption
%    using MATLAB pdepe
%
%   \$Ekkehard Holzbecher  \$Date: 2006/03/18 \$
%--------------------------------------------------------------------------
T = 16;                    % maximum time [s]
L = 8;                     % length [m]
D = 0.1;                   % diffusivity [m*m/s]
v = 0.5;                   % real fluid velocity [m/s]
lambdaf = 0;               % decay rate in fluid [1/s]
lambdas = 0;               % decay rate in solid [1/s]
theta = 0.2;               % porosity
rhob = 1200;               % porous medium bulk density [kg/m*m*m]
kappaf = 0.02;             % transition rate fluid to solid [1/s]
kappas =100000;            % transition rate solid to fluid [1/s]
c0f = 0.;                  % initial concentration in fluid [kg/m*m*m]
c0s = 0.;                  % initial concentration in solid [1]
cin = 1;                   % inflow concentration [kg/m*m*m]

M = 30;                    % number of timesteps  (>2)
N = 40;                    % number of nodes
%-------------------------- output parameters
gplot = 2;                 % =1: breakthrough curves; =2: profiles
gsurf = 0;                 % surface
gcont = 0;                 % =1: contours; =2: filled contours
ganim = 0;                 % animation of profiles; =1: single line; =2: all lines

t = linspace (T/M,T,M);    % time discretization
x = linspace (0,L,N);      % space discretization

%----------------------execution-------------------------------------------
options = odeset;
c = pdepe(0,@slowsorpde,@slowsorpic,@slowsorpbc,x,t,options,...
D,v,theta,rhob,kappaf,kappas,lambdaf,lambdas,[c0f;c0s],cin);

%---------------------- graphical output ----------------------------------
switch gplot
case 1
plot ([0 t],c(:,:,1))        % breakthrough curves
xlabel ('time'); ylabel ('concentration');
case 2
plot (x,c(:,:,1)','--')      % profiles
xlabel ('space'); ylabel ('concentration');
end
if gsurf                           % surface plot
figure; surf (x,[0 t],c(:,:,1));
xlabel ('space'); ylabel ('time'); zlabel('concentration');
end
if gcont figure; end
switch gcont
case 1
contour (x,[0 t],c(:,:,1))   % contours
grid on; xlabel ('space'); ylabel ('time');
case 2
contourf(x,[0 t],c(:,:,1))   % filled contours
colorbar; xlabel ('space'); ylabel ('time');
end
if (ganim)
[FileName,PathName] = uiputfile('*.mpg');
figure; if (ganim > 1) hold on; end
for j = 1:size(c,1)
axis manual;  plot (x,c(j,:,1),'r','LineWidth',2); YLim = [min(c0,cin) max(c0,cin)];
Anim(j) = getframe;
axis manual; plot (x,c(j,:,1),'b','LineWidth',2); YLim = [min(c0,cin) max(c0,cin)];
end
mpgwrite (Anim,colormap,[PathName '/' FileName]);     % mgwrite not standard MATLAB
movie (Anim,0);   % play animation
end

% --------------------------------------------------------------------------
function [c,f,s] = slowsorpde(x,t,u,DuDx,D,v,theta,rhob,kappaf,kappas,lambdaf,lambdas,c0,cin)
c = [1;1];
f = [D;0].*DuDx;
s = -[v;0].*DuDx - [lambdaf;lambdas].*u - ([kappaf -kappas]*u)*[1/theta;-1/rhob];

% --------------------------------------------------------------------------
function u0 = slowsorpic(x,D,v,theta,rhob,kappaf,kappas,lambdaf,lambdas,c0,cin)
u0 = c0;

% --------------------------------------------------------------------------
function [pl,ql,pr,qr] = slowsorpbc(xl,ul,xr,ur,t,D,v,theta,rhob,kappaf,kappas,lambdaf,lambdas,c0,cin)
pl = [ul(1)-cin;0];
ql = [0;1];
pr = [0;0];
qr = [1;1];```