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

Question

Ali Karim 2019-8-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/476530-how-to-define-stepchanges-variable-bc-for-different-timeintervalls-using-pdepe

回答： Ali Karim 2019-8-19

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;

5 个评论
显示 3更早的评论隐藏 3更早的评论

Ali Karim 2019-8-19

Thanks again. Could you give an example? I cannot hav eany if or for somewhere else in the code if that is what you mean?

Torsten 2019-8-19

Take a look at the last comment under

https://de.mathworks.com/matlabcentral/answers/451622-pdepe-boundary-condition-outputting-odd-results

请先登录，再进行评论。

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

Answer 1

Bill Greene 2019-8-19

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/476530-how-to-define-stepchanges-variable-bc-for-different-timeintervalls-using-pdepe#answer_388156

编辑：Bill Greene 2019-8-19

在 MATLAB Online 中打开

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