How to define stepchanges, variable BC for different timeintervalls using pdepe?
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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 Commenti
Torsten
il 19 Ago 2019
Take a look at the last comment under
https://de.mathworks.com/matlabcentral/answers/451622-pdepe-boundary-condition-outputting-odd-results
Risposta accettata
Bill Greene
il 19 Ago 2019
Modificato: Bill Greene
il 19 Ago 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
0 Commenti
Più risposte (1)
Vedere anche
Categorie
Scopri di più su Stability Analysis in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!