Why does evaluateHeatRate in PDE Solver give different answers depending on mesh size?

I have a simple half-space 2D thermal model with a top surface at 0, a specified heat flux at the bottom, insulated sides, and a hot rectangular body on the left size. I am interested in the total thermal energy lost from the body in a given time, and therefor sum the results of evaluateHeatRate for the top, right, and bottom edges of the hot rectangle. If I use the default mesh generation I get one (reasonable) number, but if I make the mesh significantly finer by setting Hmax, the calculated heat lost goes up by a factor of two. It keeps rising as I made the mesh finer and finer. Any suggestions as to why this is happening appreciated.

[Ravi Kumar]

Hi Allen,

Would you be able to share the code. Heat rate is computed based on the gradients of temperature as post processing quantity, you can check which mesh gives the accurate gradients.

Regards,

Ravi

[Asad (Mehrzad) Khoddam]

In FEM solutions (PDE Solver uses this method) for very fine meshes, the stiffness matrix will be ill-conditioned and becuase of the limiti on the the floating point digtis, we will have a numerical instability and the solution will not be accurate.

[Allen]

Here is a stripped-down script that illustrates the problem. Two mesh sizes are controlled by choosing either generateMesh statement (lines 45-46). The integral at the end (trapz) is the total heat flow out of the box at the last time step. It varies by almost a factor of two depending on which mesh is chosen, and would be even larger if a smaller mesh were chosen.

%% Thermal Model Illustrating Effect of Mesh Size on Calculated Heat Flow
% Define geometry and plot it to show edge designations.
 
width = 30000; % width of model in m
depth = 50000; % depth of model in m
% define edges of hot body
right = 10000; top = -5000; bottom = -30000; % m
%% 
% Define rectangle and plot it to show edge designations.
 
R1 = [3 4  0 width width 0 0  0 -depth -depth]';
R2 = [3 4  0 right right 0 top top bottom bottom]';
gd = [R1 R2];
sf = 'R1+R2';
ns = [82 82;49 50];
g = decsg(gd,sf,ns); % I have no idea what these last 2 lines do
 
Tmodel = createpde('thermal','transient');
geometryFromEdges(Tmodel,g);
pdegplot(Tmodel,'EdgeLabels','on',"FaceLabels","on"); axis equal;
%% 
% Set boundary conditions.
 
thermalBC(Tmodel,'Edge',1,'Temperature',0);
thermalBC(Tmodel,'Edge',2,'HeatFlux',0);
thermalBC(Tmodel,'Edge',3,'HeatFlux',1.3e6); % heat flux in J a-1 m-2
thermalBC(Tmodel,'Edge',7,'HeatFlux',0);
thermalBC(Tmodel,'Edge',8,'HeatFlux',0);
thermalBC(Tmodel,'Edge',9,'HeatFlux',0);
%% 
% Set initial conditions
 
T0 = @(location)(-0.02)*location.y; % this sets the initial T distribution
thermalIC(Tmodel,T0,'Face',1);
thermalIC(Tmodel,850,'Face',2); % set rectangular block to 850
%% 
% Set up properties of differential equation. I am using time units of years 
% (in thermal conductivity and time).
 
thermalProperties(Tmodel,'ThermalConductivity',8e7,'MassDensity',2500,'SpecificHeat',1000);
%% 
% Generate and plot mesh. Add ",'Hmax',x" to set maximum mesh size to x, 'Hmin' 
% to set minimum
 
generateMesh(Tmodel); % let Matlab set mesh
% generateMesh(Tmodel,'Hmax',500); % alternative with max mesh length 2000 m
 
pdeplot(Tmodel);
%% 
% Solve the equations and plot the results.
 
tlist = 1:1000:10000; % solution times in years
Tresults=solve(Tmodel,tlist); % solve PDE at listed times
T = Tresults.Temperature; % extract temperature
pdeplot(Tmodel,'XYData',T(:,end),'Contour','on','Levels',[0:100:1000],'ColorMap','parula');
caxis([0 1000]);
axis equal
%% 
% Evaluate heat flux out of box
 
Q4 = -evaluateHeatRate(Tresults,'Edge',4); % negative sign to make heat flowing out of box positive
Q5 = -evaluateHeatRate(Tresults,'Edge',5);
Q6 = -evaluateHeatRate(Tresults,'Edge',6);
Q = Q4+Q5+Q6; % J/a
 
disp('Values and sum of heat rate')
Q
trapz(tlist,Q)

[Ravi Kumar]

Thanks for the code. I think your refined mesh results are (relatively) more accurate then the default mesh one. You are also calculating heat rate at an interface that has significant jump in the temperature, for example, about 800 on the top edge of the hot body. Neither the default mesh nor the mesh with Hmax set to 500 is sufficient to resolve the gradients at these interfaces. Any small error in gradients would be amplified by the large thermal conductivity and length scale of the problem domain. Ideally you need very fine mesh only around these boundary and coarse mesh away from it to compute gradient accurately. It is not easy to generate such mesh in PDE Toolbox.

Unrelated to your question, I also notice that you convert thermal conductivity to account for your time-scale. Should you also be making such adjustments to specific heat, as Joule does have dimension of time embedded in it?

[Allen]

Thank you for the comments. I have found a more stable solution to the heat loss problem by using interpolateTemperature across the hot body to determine the temperature loss (and thus heat loss). This value, too, increases as the mesh size decreases, but at a much slower rate. evaluateHeatFlux numbers oscillate more and more as the mesh size goes down, in keeping with Asad's comment above about stiffness.

[Allen]

P.S. I have switched the calculation to seconds for consistency, but the answers are the same. The J terms cancel out as they are shared by k and c, but better to keep everything consistent.