PDEPE function - BCFUN error

2023 11 20
1 个回答
回答已采纳
更新时间：2023 11 21
3 次查看（30 天）

0 个投票

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

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

 采纳的回答

Torsten
Torsten 2023-11-20
编辑：Torsten 2023-11-20

0 个投票

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 个评论
显示 12更早的评论 隐藏 12更早的评论

Sanley Guerrier
Sanley Guerrier 2023-11-20
Yes, hc, env(:,1), Hs, and temp(:,6) are vectors of same length.
Sanley Guerrier
Sanley Guerrier 2023-11-20
I tried size(ptop), it didn't work.
Torsten
Torsten 2023-11-20
Yes, hc, env(:,1), Hs, and temp(:,6) are vectors of same length.
ptop = interp1(t_mesh,hc,t)*(Ttop-interp1(t_mesh,env(:,1),t))+ eps*sigma*Ttop^4-A*interp1(t_mesh,Hs,t);
qtop = k
if you want to set the condition
interp1(t_mesh,hc,t)*(Ttop-interp1(t_mesh,env(:,1),t))+ eps*sigma*Ttop^4-A*interp1(t_mesh,Hs,t) + k*alpha*dT/dz = 0
at the top.
Sanley Guerrier
Sanley Guerrier 2023-11-20
Thank you,
This command is working now but the output is only one row. Something is still not working.
ptop = interp1(t_mesh,hc,t)*(Ttop-interp1(t_mesh,env(:,1),t))+ eps*sigma*Ttop^4-A*interp1(t_mesh,Hs,t);
qtop = k
Torsten
Torsten 2023-11-20
编辑：Torsten 2023-11-20
This command is working now but the output is only one row.
What do you mean ? ptop must be one single value, namely you boundary condition evaluated at time t.
Sanley Guerrier
Sanley Guerrier 2023-11-20
I meant the solution T should be a matrix length (z_mesh X t_mesh) but I run the code it only returned a row vector T.
Sanley Guerrier
Sanley Guerrier 2023-11-21
移动：Torsten 2023-11-21
When I used this top boundary condition it works fine. But when I switched to other one it didn't work.
function [ptop, qtop, pbot, qbot] = bcfun(ztop, Ttop, zbot, Tbot, t)
ptop = Ttop - interp1(t_mesh, temp(:,1), t, 'linear');
qtop = 0;
pbot = Tbot - interp1(t_mesh, temp(:,6), t, 'linear');
qbot = 0;
end
Torsten
Torsten 2023-11-21
I meant the solution T should be a matrix length (z_mesh X t_mesh) but I run the code it only returned a row vector T.
Don't you get a message that the integration stopped because the stepsize became too small or something similar ? If integration didn't proceed till the second output time, the only thing that pdepe returns are the initial conditions for T.
Sanley Guerrier
Sanley Guerrier 2023-11-21
Warning: Failure at t=1.958356e+00. Unable to meet integration tolerances without reducing the step size below
the smallest value allowed (3.552714e-15) at time t.
> In ode15s (line 725)
In pdepe (line 291)
In run_model (line 53)
In driver (line 58)
Warning: Time integration has failed. Solution is available at requested time points up to t=1.958333e+00.
> In pdepe (line 305)
In run_model (line 53)
In driver (line 58)
Elapsed time is 18.348141 seconds.
Torsten
Torsten 2023-11-21
编辑：Torsten 2023-11-21
Yes, that's what I meant - my guess is that t_mesh(2) > 1.958356, isn't it ? So you have an explanation why you only get one row for T.
Concerning the error I cannot help. It seems the reason is the boundary condition at the top you supply. Check whether it's correct. Note that it reads ptop + qtop*f = 0 where f in your case is equal to alpha*dT/dz. Thus you have ptop + k*alpha*dT/dz = 0 as boundary condition.
Sanley Guerrier
Sanley Guerrier 2023-11-21
Thank you, I got it.
Sanley Guerrier
Sanley Guerrier 2023-11-21
I forgot to convert hc to m/hr because every was in m/hr.
Do you know why when I used qtop = -k works but qtop = k doesn't?
Torsten
Torsten 2023-11-21
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 2023-11-21
That makes sence.
Thanks a lot for helping out.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Structural Mechanics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by