PDE Toolbox Tolerance Error - Non-Constant Non-Linear Neumann Condition

Hello, I have been trying the to model a 1-d diffusion problem, with the following governing equation:

The problem has a non-linear neumman boundary condition where heat flux is caused by convection and radiation.

This boundary condition is futher complicated by the fact that the temperature of the surrounding fluid changes over time:

I have tried numerous methods but to no avail - however, following tutorials, I feel I have gotten very close using the PDE toolbox.

%% PROPERTIES AND VARIABLES

% MATERIAL PROPERTIES

k = 1.3; % thermal conductivity, W/(m-K)

rho = 2705; % density, kg/m^3

C_p = 0.91*1000; % specific heat, J/(kg-K)

alpha = k/(rho*C_p); % thermal diffusivty

T_ambient = 20;

%% PDE MODEL SETUP

% Create a generic PDE analysis model of system size 1.

numberOfPDE = 1;

model = createpde(numberOfPDE);

% DEFINE GEOMETERY

radius = 0.18; % meter

height = radius; % meter

surface_area = 2*pi*radius;

g = decsg([3 4 0 0 radius radius -height/2 height/2 height/2 -height/2]');

geometryFromEdges(model,g); % Convert into PDE geometry.

% Plot of geometry with boundary ID

figure;

pdegplot(model,'EdgeLabels','on');

title 'Geometry With Edge Labels Displayed';

%% BOUNDARY CONDITIONS

applyBoundaryCondition(model, 'neumann', 'edge', 3, 'g', @mygfun)

%% MODEL COOEFICIENTS AND INITIAL CONDITIONS

% EQUATION OF FORM:

% m*d^2u/dt^2 + d*du/dt - div(c*grad(u)) + a*u = f

c = 1; % 2nd Order Space Derrivative Coefficient

d = -1/alpha; % 1st Order Time Derrivative Coefficient

m = 0; % 2nd Order Time Derrivative Coefficient

a = 0; % Solution Multiplier Coefficient

f = 0; % Source Coefficient

specifyCoefficients(model,'m',m,'d',d,'c',c,'a',a,'f',f);

% Initial Conditions.

setInitialConditions(model,T_ambient);

%% MESH

hmax = .01; % element size

msh = generateMesh(model,'Hmax',hmax);

figure;

pdeplot(model);

axis equal

title 'Plate With Triangular Element Mesh'

xlabel 'X-coordinate, meters'

ylabel 'Y-coordinate, meters'

%% SOLVE

p = msh.Nodes;

endTime = 5000;

tlist = 0:50:endTime;

numNodes = size(p,2);

u0(1:numNodes) = T_ambient; % Define initial condition

% Tolerances

model.SolverOptions.RelativeTolerance = 1.0e-3;

model.SolverOptions.AbsoluteTolerance = 1.0e-4;

R = solvepde(model,tlist);

u = R.NodalSolution;

The boundary condition uses a function (mygfunc) to calculate the heat flux:

function funOut = mygfun (location,state)

sigma = 5.670373e-8; % Stefan-Boltzmann constant, W/(m^2-K^4)

phi = 1; % Configuration factor.

e_m = 0.9; % Surface emissivity.

e_f = 1; % Fire emissivity.

alpha_convection = 25;

radius = 0.18; % meter

height = radius; % meter

surface_area = 2*pi*radius;

k = 1.3; % thermal conductivity, W/(m-K)

T_gas = 20 + 345*log10(8*(state.time/60) + 1);

% Room Temperature is defined by the standard fire curve, C.

% It is a function of time.

% Divided by 60 to convert from seconds to minutes.

% HEAT FLUX

h_c = alpha_convection*(T_gas - state.u);

h_r = phi*e_m*e_f*sigma((T_gas + 273)^4 - (state.u + 273)^4);

h_total = h_c + h_r;

heatFlux = h_total*surface_area;

funOut = heatFlux/(-k);

end

However when I run the solver I recieve the following tollerance error:

I am really unfamiliar with the PDE toolbox and no expert at matlab - and ive been stuck on this for quite a while so I would really really really appreciate any help I can get.

Best,

Victor

Answer 1

[victor noorani]

As it takes 14 seconds to diverge I ran it for the first 10 and looked at the results.

This doesn't make sense as heat flux should be consistent along that entire edge. But maybe it indicates something is wrong with my BC.