finite element method – Boundary Value of Function DIsagrees with Neumann Boundary Heat Equation

I’m trying to solve a time-dependent 2D heat equation with a source, initial specified heat distribution and heat flux loss at the boundary which goes with the temperature difference. Thus far, I’ve managed to get the prior two parts working but the source term seems to make the derivative at the boundary disagree with the value specified by the Neumann conditions with no error message. My code is shown below.

The first part is setting up the mesh for the FEM:

dist = 0.0025; (*0.0021*)
RTSpacing = 0.001;  (*0.002*)
DResistor = <|
  "Lx" -> 1.55*10^-3`, "Ly" -> 0.2*10^-3`, "Lz" -> 2.2*10^-3|>;
DThermistor = <|"Lx" -> 3.2*10^-3`, "Ly" -> 0.45*10^-3`, 
   "Lz" -> 4.5*10^-3|>;
(CapitalOmega)Resistor = 
 CreateRect({0, RTSpacing*0.5` + getLy(DResistor)}, 
  DResistor); (*Can also do rounding radius*)


(CapitalOmega)Thermistor = 
 CreateRect({0, -RTSpacing*0.5` - getLy(DThermistor)}, DThermistor);
(CapitalOmega)Air = (CapitalOmega)Disk(
  1); (*Original is 0.2*10^-2*)
(CapitalOmega)Plot = 
(CapitalOmega)Square(dist);
TotalMesh = 
  CombineDomaintoMesh((CapitalOmega)Air, {(CapitalOmega)Resistor, 
(CapitalOmega)Thermistor});
PlotMesh = 
  CombineDomaintoMesh((CapitalOmega)Plot, {(CapitalOmega)Resistor, 
(CapitalOmega)Thermistor});
PlotMesh("Wireframe")

Which produces a mesh that looks like this:
Wireframe plot of the mesh used to generate the solution and plot

CombineDomainintoMesh is defined as shown below. It takes a list of domains for the innershapes and has worked fine in numerous examples thus far:

CombineDomaintoMesh(outershape_, innershapes_) := 
 Module({innercoords, coords, outermesh, BMesh},
  outermesh = ToBoundaryMesh(outershape);
  innercoords = Flatten(ToBoundaryMesh(#)((1)) & /@ innershapes, 1);
  coords = Join(outermesh((1)), innercoords);
  BMesh = 
   ToBoundaryMesh("Coordinates" -> coords, 
    "BoundaryElements" -> outermesh((3)), 
    "PointElements" -> outermesh((4)));
  ToElementMesh(BMesh)
  )

Then some functions are defined to handle the piecewise values of the properties like the source and the heat capacities etc.

PropertyFunction(x_, y_, 
   properties_) := (properties((1)) - properties((3)))*
    RectPiecewise(x, 
     y, (CapitalOmega)Resistor) + (properties((2)) - 
      properties((3)))*
    RectPiecewise(x, y, (CapitalOmega)Thermistor) + 
   properties((3))*DiskPiecewise(x, y, (CapitalOmega)Air);
HeatSource(x_, y_, S_) := 
  S*RectPiecewise(x, y, (CapitalOmega)Resistor) ;
(*Properties are in order of resistor, thermistor, paste*)

Then we have the initial conditions and material properties being substituted.

T0 = 294;
ThermConduct = {30, 30, 0.0259};
(*0.0259*)
VolHeatCapacity = {2.43*10^6, 2.43*10^6, 1210};
(*1210*)
VolumeResistor = getVolume(DResistor);
Voltage = 13;
Resistance = 82;
SourceStrength = Voltage^2/(Resistance*VolumeResistor);
(*Variables and parameters for heat equation*)

vars = {T(t, x, y), t, {x, y}};
pars = <|"ThermalConductivity" -> 
    PropertyFunction(x, y, ThermConduct)*IdentityMatrix(2),
   "SpecificHeatCapacity" -> PropertyFunction(x, y, VolHeatCapacity),
   "HeatSource" -> HeatSource(x, y, SourceStrength)|>;

Then we define the Neumann Boundary Conditions below and plug this into the NDSolveValue to get the solution.

(*Defining Neumann Boundary Conditions*)
HeatTransferCoeff = 1000;
(CapitalGamma) = 
 HeatTransferValue({x == dist, y == dist, y == -dist, x == -dist}, 
  vars, pars, <|"HeatTransferCoefficient" -> HeatTransferCoeff, 
   "AmbientTemperature" -> T0|>)

pde = {{HeatTransferPDEComponent(vars, pars) == (CapitalGamma)}, 
   T(0, x, y) == T0};
thermsol3 = 
 Refine(NDSolveValue(pde, 
   T, {x, y} (Element) PlotMesh, {t, 0.0, 1000}), {Element(T, 
    PositiveReals)})

This produces a reasonable looking solution but taking the difference between the Neumann condition and the actual derivative at the boundary reveals that the difference is significant.

Plot(Head(D(thermsol3(t, x, y), y))(k, 0, -dist) - 
  HeatTransferCoeff*(294 - thermsol3(k, 0, -dist)), {k, 0, 
  1000})

Plot of the difference in set value vs. actual derivative value over time

The actual slope seems to just be vastly larger than the value that should be set at the boundary. As I mentioned earlier, this is not an issue without the source term and for a given domain, one can predict where the offset will be shifted to. I’ve read the documentation for the NeumannBoundaryValue (https://reference.wolfram.com/language/ref/NeumannValue.html) and it seems like the source term should be treated as their “f” parameter and effectively ignored as part of the boundary conditions where it is zero anyway. Is this something to do with the mesh not being fine enough around the boundary? It also seems as though the boundary is specifying the heat flux not the temperaure gradient. However, introducing the thermal conductivity to convert this to a temperature gradient doesn’t rectify it either.

Any help or insight on specifying the NeumannValues would be greatly appreciated, thanks!