Step Function as Input Boundary Condition to MATLAB PDEPE solver

2 次查看（过去 30 天）
Siamack
Siamack 2018-3-19
评论： Bill Greene 2018-3-22
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;

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

采纳的回答

Bill Greene
Bill Greene 2018-3-22
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 个评论
Siamack
Siamack 2018-3-22
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 2018-3-22
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by