MATLAB Answers

Step Function as Input Boundary Condition to MATLAB PDEPE solver

14 views (last 30 days)
Siamack
Siamack on 19 Mar 2018
Commented: Bill Greene on 22 Mar 2018
Dear Matlab Users I have the following code to find the transient solution to the 1D heat conduction equation. The code below intends to use a rectangular heat pulse as the initial temperature but i receive an error ( Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t). Can anyone tell me what is issue with the boundary condition as I set it and if there is a better way to deal with the matter.
function out = parabolic_silicon(~)
global rho cp k T_i T_o
L = 10; % [cm]
k = 1.3; % [W/cmK]
rho = 2.33; % [g/cm^3]
cp = 0.7; % [J/gK]
T_i = 295; % [K]
T_o = T_i + 5.2e-3; % [K]
t_end= 1; % [s]
m = 0;
x = linspace(0,L,500);
t = linspace(0,t_end,500);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
Temperature = sol(:,:,1);
L = strcat('t = ',num2str(t(:)),' Seconds');
plot(x,sol);
% legend(L,'location','northeast');
% % --------------------------------------------------------------
% OUTPUT
out = {x,t,sol};
% % --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global rho cp k
c = rho*cp;
f = k*DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 295;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global T_i T_o
pl = ul - (T_i + (T_o*(heaviside(t) - heaviside(t - 0.5))));
ql = 0;
pr = ur - T_i;
qr = 0;

  0 Comments

Sign in to comment.

Accepted Answer

Bill Greene
Bill Greene on 22 Mar 2018
I suggest you smooth the sharp corners in your boundary condition profile using a Heaviside function similar to the code below. Note, you can change the sharpness of the corners by varying the r variable.
function k=heavisideSmooth(t)
k1 = 0; k2 = 1; c = 0;
r = 1e3;
ecr = exp(c*r);
erx = exp(r*t);
k = (k1*ecr + k2*erx)./(ecr+erx);
end

  2 Comments

Siamack
Siamack on 22 Mar 2018
Thanks alot Bill, I used the smooth function and it worked. Just one more question and that of: your suffestion was based on the error message regarding the integration tolerance?
Regards Siamack
Bill Greene
Bill Greene on 22 Mar 2018
The error message indicates the ODE solver wasn't able to take even a tiny, first step. So I assumed that was due to the sharp discontinuity in the BC at t=0.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!