PDEPE function - BCFUN error
2 views (last 30 days)
Show older comments
Sanley Guerrier
on 20 Nov 2023
Commented: Sanley Guerrier
on 21 Nov 2023
Hello all,
Can someone help me with the BCFUN error?
Thank you!
Error using pdepe
Unexpected output of BCFUN. For this problem BCFUN must return four column vectors of length 1.
Error in run_model (line 53)
sol = pdepe(0, @pdefun, @icfun, @bcfun, z_mesh, t_mesh);
Error in driver (line 58)
[z_mesh, t_mesh, T] = run_model(alpha, z_sensor, time, temp, env, hc, k, eps, sigma, A, Hs);
function [z_mesh, t_mesh, T] = run_model(alpha, z_sensor, time, t_mesh, temp, env, hc, k, eps, sigma, A,Hs)
% Create the time mesh using the measurement times [days].
% t_mesh is shifted to start at 0.
t_mesh = convertTo(time, 'datenum') - convertTo(time(1), 'datenum');
% Create the space mesh using even spacing between the top and
% bottom sensors, plus the locations of the intermediate sensors.
z_mesh = unique([linspace(z_sensor(1), z_sensor(end), 100), z_sensor]);
% Solve the heat equation.
sol = pdepe(0, @pdefun, @icfun, @bcfun, z_mesh, t_mesh);
% Extract the first solution component as T.
T = sol(:,:,1);
%------------------------------------------------
% pdefun
%------------------------------------------------
function [c, f, s] = pdefun(z, t, T, dTdz)
c = 1;
f = alpha * dTdz;
s = 0;
end
%------------------------------------------------
% icfun
%
% Notes:
% o The inital temperature profile is given by
% interpolating the initial field data.
%------------------------------------------------
function T0 = icfun(z)
T0 = interp1(z_sensor, temp(1, :), z, 'linear');
end
%------------------------------------------------
% bcfun
%
% Notes:
% o The boundary conditions through time are
% given by the flux and temperature at the bottom sensors.
%------------------------------------------------
function [ptop, qtop, pbot, qbot] = bcfun(ztop, Ttop, zbot, Tbot,t)
% hc and Hs are vectors
ptop = hc.*(Ttop-env(:,1))+ eps*sigma*Ttop^4-A*Hs;
qtop = k;
pbot = Tbot - interp1(t_mesh, temp(:,6), t, 'linear');
qbot = 0;
end
end
0 Comments
Accepted Answer
Torsten
on 20 Nov 2023
Edited: Torsten
on 20 Nov 2023
Use
function [ptop, qtop, pbot, qbot] = bcfun(ztop, Ttop, zbot, Tbot,t)
% hc and Hs are vectors
ptop = hc.*(Ttop-env(:,1))+ eps*sigma*Ttop^4-A*Hs;
qtop = k;
pbot = Tbot - interp1(t_mesh, temp(:,6), t, 'linear');
qbot = 0;
size(ptop)
end
and see what the size of ptop is. It should be 1x1, but I doubt it because I guess that env(:,1) is a vector.
14 Comments
Torsten
on 21 Nov 2023
You must decide why the boundary condition is
ptop - k*alpha*dT/dz = 0
and not
ptop + k*alpha*dT/dz = 0
It's equivalent to an inversion of the direction of the heat flux at the top.
More Answers (0)
See Also
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!