MATLAB: How to define stepchanges, variable BC for different timeintervalls using pdepe

different bc for different time inervallspdepestepchange at bc

I am trying to define stepchanges of relative humidity at the boundary of a slab varying from 0.9 to 0.3 every 12 hours without any progress. The code is working without any issues for only one single BC but not otherwise. Any idea how I can get it? Thanks!
function parabolicmoisture
global  tp a_v tend v_s delta exi  
L=0.15;     
k=0.028;     
rho=156;
cp=990;
mu=4.25;
D=25*10^-6;
delta=D/mu;
v_s=17.28e-3;
exi=5/0.65;
a_v=(delta*v_s)/exi;
m = 0;
tp=24*3600;
tend=2*tp; 
x = linspace(0,L,1000);
t = linspace(0,tend,2000);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
RH = sol(:,:,1);
% A solution profile can also be illuminating.
plot(x,RH(end,:),'r--','LineWidth',2)
hold on
plot(x,RH(3*end/4,:),'b.-','LineWidth',1)
hold on
plot(x,RH(2*end/4,:),'k','LineWidth',1)
hold on
plot(x,RH(1*end/4,:),'*')
title(strcat('Solution at t = ', num2str(tend)))
legend('48h','36h','24h','12h')
xlabel('Distance x')
ylabel('RH (%)')
 
% %Plot surface temperature vs. time
% figure, plot(t,RH(:,1))
% title('Surface RH')
% xlabel('Time (s)')
% ylabel('RH (%)')
% --------------------------------------------------------------


 
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global a_v v_s delta exi
c=1/a_v;
f = 1*DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 0;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global  tend 
% 
   if  t<(tend/4)
         pl = ul-0.9;
   elseif  t>=(tend/4) || t<tend/2
          pl = ul-0.3;
   elseif  t>=tend/2  || t<3*tend/4
         pl = ul-0.9;
   elseif t>=3*tend/4
       pl = ul-0.3;
   end
 
ql = 0;
pr = ur;
qr = 0;

Best Answer

The reason that your code didn't show a change in the BC with time is due to the programming error that Torsten pointed out. If you had simply made the corrections he suggested (as I did) pdepe would have failed at the time of the first BC change. The problem is the sharp step change in BC value at the transition times. If you simply allow the BC to change over a small, but non-zero, time interval, pdepe will be able to obtain a solution. In the code I show below, I implemented this strategy and arbitrarily picked 10 seconds as the BC transition time. I also used the matlab function interp1 to substantially simplify the coding.
function parabolicmoisture
global  tp a_v tend v_s delta exi  
L=0.15;     
k=0.028;     
rho=156;
cp=990;
mu=4.25;
D=25*10^-6;
delta=D/mu;
v_s=17.28e-3;
exi=5/0.65;
a_v=(delta*v_s)/exi;
m = 0;
tp=24*3600;
tend=2*tp;
d=10; % small time over which bc changes occur
bcTimes=[0 tend/4 tend/4+d tend/2 tend/2+d 3*tend/4 3*tend/4+d tend];
bcVals=[.9 .9 .3 .3 .9 .9 .3 .3];
x = linspace(0,L,100);
t = linspace(0,tend,200);
% test the bc table
%figure; plot(t, interp1(bcTimes, bcVals, t)); return;
bcFunc=@(xl,ul,xr,ur,t) pdex1bc(xl,ul,xr,ur,t,bcTimes,bcVals);
sol = pdepe(m,@pdex1pde,@pdex1ic,bcFunc,x,t);
RH = sol(:,:,1);
figure; plot(t/tend, RH(:,1)); grid;
title('Solution at left end as a function of normalized time.');
% A solution profile can also be illuminating.
plot(x,RH(end,:),'r--','LineWidth',2)
hold on
plot(x,RH(3*end/4,:),'b.-','LineWidth',1)
hold on
plot(x,RH(2*end/4,:),'k','LineWidth',1)
hold on
plot(x,RH(1*end/4,:),'*')
title(strcat('Solution at t = ', num2str(tend)))
legend('48h','36h','24h','12h')
xlabel('Distance x')
ylabel('RH (%)')
 
% %Plot surface temperature vs. time
% figure, plot(t,RH(:,1))
% title('Surface RH')
% xlabel('Time (s)')
% ylabel('RH (%)')
% --------------------------------------------------------------


end
 
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global a_v v_s delta exi
c=1/a_v;
f = 1*DuDx;
s = 0;
end
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 0;
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t,bcTimes,bcVals)
pl = ul-interp1(bcTimes, bcVals, t);
ql = 0;
pr = ur;
qr = 0;
end
Related Question