Clear Filters
Clear Filters

PDEPE function - BCFUN error

3 views (last 30 days)
Sanley Guerrier
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

Accepted Answer

Torsten
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
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.
Sanley Guerrier
Sanley Guerrier on 21 Nov 2023
That makes sence.
Thanks a lot for helping out.

Sign in to comment.

More Answers (0)

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!