# How to define stepchanges, variable BC for different timeintervalls using pdepe?

9 views (last 30 days)
Ali Karim on 17 Aug 2019
Answered: Ali Karim on 19 Aug 2019
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;
Torsten on 19 Aug 2019
Take a look at the last comment under

Bill Greene on 19 Aug 2019
Edited: Bill Greene on 19 Aug 2019
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

Ali Karim on 19 Aug 2019
Thank you very much. I see the problem with my code now. Thank you Torsten and Bill Greene.