I have a steady state heat transfer problem trying to model a thin film heater wedged between two significantly thicker slabs. The slab thicknesses may be as much as three orders of magnitude larger than the heater. I want temperature resolution within the slabs and a temperature of the heater but I don’t need temperature resolution within the film per se. I tried treating the entire thin film as dissipating energy using an internal heat source command of this form [internalHeatSource(thermalmodel,6000,"Face",2)]. My checks on model validity included looking for higher temperature within the heater compared to the two bounding slabs using a command of this form [pdeplot(thermalmodel,"XYData",thermalresults.Temperature(:,end), ... "Contour","on",... "ColorMap","jet")]. I also checked energy conservation in the steady state by comparing the known amount of power dissipated within the heater to the energy rate leaving the boundaries of the system using commands of the form [Qn1= evaluateHeatRate(thermalresults, "Edge",1)] for each edge and summing up these energy rates..
This approach was not successful. I did not see significant temperature increases that I would expect for heater power levels that I used. Conservation of energy was not observed. I plotted the mesh results and in some or all cases would get strange looking results. In some or all cases Matlab returned the following feedback: “Warning: Found elements with poor shape quality.” I tried manipulating the mesh size to get finer resolution using commands of the form [mesh= generateMesh(thermalmodel, 'Hmax', 1e-6, 'Hmin', 1e-7)]. I was not successful.
I feel like the problem is related to the mismatch in mesh resolution that may be required in the heater compared to the larger slabs, but I don’t know if that’s true and I don’t know how to handle it if that is true.
The code is included below.
Any help would be greatly appreciated.
john
%% Create a thermal model for transient analysis
thermalmodel = createpde("thermal","steadystate");
% Define the geometry components
Rfilmbottom = [3,4,[-5.0, 5.0, 5.0, -5.0, 0.0, 0.0, 0.124375,0.124375]/1000]';
Rfilmhtr = [3,4,[-5.0, 5.0, 5.0, -5.0, 0.124375,0.124375,0.125625,0.125625]/1000]';
Rfilmtop= [3,4,[-5.0, 5.0, 5.0, -5.0, 0.125625,0.125625, 0.25, 0.25]/1000]';
% Combine the geometry data
gdm = [ Rfilmbottom Rfilmhtr Rfilmtop ];
ns = char('Rfilmbottom', 'Rfilmhtr', 'Rfilmtop');
g = decsg(gdm, 'Rfilmbottom + Rfilmhtr + Rfilmtop', ns');
% Create geometry from edges
geometryFromEdges(thermalmodel, g);
% Create a larger figure window
figure('Position', [100, 100, 1200, 800]);
% Figure 1
% Plot with edge labels
pdegplot(thermalmodel, "EdgeLabels", "on", "FaceLabels", "on");
%Figure 2
figure
% Adjust limits and aspect ratio for better visibility
% Overall view
pdegplot(thermalmodel, "EdgeLabels", "on", "FaceLabels", "on");
xlim([-5.0e-3 5.0e-3]);
ylim([0.0e-3 0.5e-3])
daspect([50 1 1]);
% Around heater
%pdegplot(thermalmodel, "EdgeLabels", "on", "FaceLabels", "on");
%xlim([-5.0e-3 5.0e-3]);
%ylim([1.1e-4 1.4e-4])
%daspect([500 1 1]);
% Properties assigned assuming corresponding order as per above
% Substrate
%thermalProperties(thermalmodel,"ThermalConductivity",0.050);
%Film
thermalProperties(thermalmodel,"ThermalConductivity",0.025);
% Boundary conditions
% Edge 1 is bottom edge of film
thermalBC(thermalmodel,"Edge",1,"Temperature",30);
% Edge 2 is top edge of film
thermalBC(thermalmodel,"Edge",2,"Temperature",30);
% Edge 3 is bottom edge of heater film- no BC for continuity
%thermalBC(thermalmodel,"Edge",3,"HeatFlux",0);
% Edge 4 is top edge of heater film- no BC for continuity
%thermalBC(thermalmodel,"Edge",4,"HeatFlux",0);
% Edge 5 is right edge of bottom film- adiabatic
thermalBC(thermalmodel,"Edge",5,"HeatFlux",0);
% Edge 6 is right edge of heater film- adiabatic
thermalBC(thermalmodel,"Edge",6,"HeatFlux",0);
% Edge 7 is right edge of top film- adiabatic
thermalBC(thermalmodel,"Edge",7,"HeatFlux",0);
% Edge 8 is left edge of bottom film- adiabatic
thermalBC(thermalmodel,"Edge",8,"HeatFlux",0);
% Edge 9 is left edge of heater film- adiabatic
thermalBC(thermalmodel,"Edge",9,"HeatFlux",0);
% Edge 10 is left edge of top film- adiabatic
thermalBC(thermalmodel,"Edge",10,"HeatFlux",0);
internalHeatSource(thermalmodel,6000,"Face",2);
mesh= generateMesh(thermalmodel);
%mesh= generateMesh(thermalmodel, 'Hmax', 5e-8, 'Hmin', 5e-9);
%mesh= generateMesh(thermalmodel, 'Hmax', 1e-6, 'Hmin', 1e-7);
% Figure 5
figure
pdemesh(thermalmodel)
pdeplot(mesh)
%xlim([-1.0e-3 1.0e-3]);
%ylim([-2.0e-3 1.0e-3])
%daspect([1 1 1]);
xlim([-5.0e-3 5.0e-3]);
ylim([0.0e-3 0.5e-3])
daspect([50 1 1]);
title("Mesh with Refined Elements");
%Warning: Found elements with poor shape quality.
thermalresults = solve(thermalmodel)
T= thermalresults.Temperature;
msh= thermalresults.Mesh;
[qx,qy] = evaluateHeatFlux(thermalresults);
%figure
%pdeplot(thermalmodel,"XYData",thermalresults.Temperature(:,end), ...
% "Contour","on",...
% "FlowData",[qx(:,end),qy(:,end)], ...
%"ColorMap","jet")
%title("Temperature (C) & Heat flux fields (W/m2)")
%xlabel("Width (m)")
%ylabel("Height (m)")
% Figure 6
figure
pdeplot(thermalmodel,"XYData",thermalresults.Temperature(:,end), ...
"Contour","on",...
"ColorMap","jet")
title("Temperature field (C)")
xlabel("Width (m)")
ylabel("Height (m)")
xlim([-1.0e-3 1.0e-3]);
ylim([-2.0e-3 1.0e-3])
daspect([1 1 1]);
% Figure 7
figure
pdeplot(thermalmodel,"FlowData",[qx(:,end), qy(:,end)])
title("Heat flux field (W/m2)")
xlabel("Width (m)")
ylabel("Height (m)")
xlim([-1.0e-3 1.0e-3]);
ylim([-2.0e-3 1.0e-3])
daspect([1 1 1]);
%isp(qx(:,end))
%disp(qy(:,end))
Qn1= evaluateHeatRate(thermalresults, "Edge",1)
Qn2= evaluateHeatRate(thermalresults, "Edge",2)
Qn3= evaluateHeatRate(thermalresults, "Edge",3)
Qn4= evaluateHeatRate(thermalresults, "Edge",4)
Qn5= evaluateHeatRate(thermalresults, "Edge",5)
Qn6= evaluateHeatRate(thermalresults, "Edge",6)
Qn7= evaluateHeatRate(thermalresults, "Edge",7)
Qn8= evaluateHeatRate(thermalresults, "Edge",8)
Qn9= evaluateHeatRate(thermalresults, "Edge",9)
Qn10= evaluateHeatRate(thermalresults, "Edge",10)
Edgetotalwatts= Qn1+Qn2+Qn4+Qn5+Qn6+Qn7+Qn8+Qn9+Qn10
