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

3 次查看（过去 30 天）
victor noorani
victor noorani 2022-7-12
回答： victor noorani 2022-7-12
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

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

回答（1 个）

victor noorani
victor noorani 2022-7-12
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.

类别

在 Help Center 和 File Exchange 中查找有关 Heat Transfer 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by