I am trying to implement a radiation boundary condition using the new solver for the heat transfer applicaition for the PDE Toolbox. I reduced the problem to just a square domain with fixed temperture on the bottom and radiation from the remaining edges. The temperture only increases slightly in the wrong direction. Please let me know what I am missing! Here is the code:
% Given
Lx = 0.1; % m, maximum length in x-direction
Ly = 0.1; % m, maximum length in y-direction
Hmax = Lx/10; % m, maximum mesh spacing
kt = 16; % W/m-K, thermal conductivity
epsilon = 0.9; % emissivity
T_inf = 10; % K, temperature of surroundings
T_b = 100; % K, temperture of base
% Create PDE thermal model container
thermalmodel = createpde('thermal','steadystate');
% Create 2D geometry and append geometry to model
R1 = [3;4;0;Lx;Lx;0;0;0;Ly;Ly];
gd = R1; % geometry description matrix
sf = 'R1'; % set formula
ns = char('R1')'; % name space matrix
dl = decsg(gd,sf,ns); % decomposed geometry matrix
pg = geometryFromEdges(thermalmodel,dl);
% Set boundary conditions and thermal properties
BC1 = thermalBC(thermalmodel,'Edge',1,'Temperature',T_b);
BC2 = thermalBC(thermalmodel,'Edge',2:4,'Emissivity',epsilon,...
'AmbientTemperature',T_inf);
thermalmodel.StefanBoltzmannConstant = 5.670373E-8; % W/m^2-K^4
mtl = thermalProperties(thermalmodel,'ThermalConductivity', kt);
% Generate mesh, solve problem, and plot results
msh = generateMesh(thermalmodel,'Hmax',Hmax);
thermalmodel.SolverOptions.ResidualTolerance = 1.0e-05;
results = solve(thermalmodel);
pdeplot(thermalmodel,'XYData',results.Temperature,'ColorMap','jet');
colorbar
