How to find mean temperature over area or volume from PDE solution?

I am using the PDE toolbox for a steady-state thermal analysis. The solver provides the temperature at each node and I can use evaluateHeatRate to get the heat flux through a face but I cannot find a way of getting the area- or volume-weighted mean temperature through the body. Is there a function that does this (rather than simply averaging the nodal values with no weighting function)?

(My body has 3 temperature boundaries. I want to solve using 2 independent temperature vectors, get the mean temperature for each and then define the body-mean temperature as Tbar = a*T1 + b*T2 + c*T3 so I can evaluate the pseudo-steady heat capacities dU/dT1, dU/dT2 etc).

Thanks

Roger

[darova]

Can you export data to MATLAB workspace (temperature, nodes and volumes)?

Answer 1

[Roger Moss]

Hi darova

Your comment made me think about the workspace variables - thanks for the nudge! Yes, of course one has nodal coordinates and temperatures from the PDE toolbox solution - the difficulty is the volumes. I tend to look at the "See also" lists when doing Help on functions which I am learning but these gave me no clues.

Then looking at the PDE toolbox functions I noticed area() and volume() and now I have some code that works, see attached code. For a 2D mesh, area() gives the area of each mesh element; one then needs to find the centroid coordinates for the element and interpolate the solution at these coordinates.

The example geometry I've used looks like this. There is convective heat transfer to the upper face only (heat input on left side, rejection on right).

My code works for 2D steady state meshes, either (x,y) or axisymmetric (r,z) (for which it converts the areas to volumes). I've only done a limited mount of testing but the answers seem sensible.

More generally, this might be a useful feature for inclusion in the toolbox, with options for whole mesh or per face output (and 3D equivalents).

Roger

[Roger Moss]

I've just realised one cannot easily view the attached code. Here is mean_temperature_example.m from above:

% How to find area or volume-weighted mean temperature for a 2D PDE model.
% RWM, April 2020

% uncomment/comment model type as required!
model=createpde('thermal','steadystate-axisymmetric');
%model=createpde('thermal','steadystate');

% Geometry is 2 rectangles: top-left, top-right
gd = [3 4, 0   0.3 0.3  0,  0     0    -0.5 -0.5
      3 4, 0.3 2   2  0.3,  0     0    -0.5 -0.5]';
    
[g, bt] = decsg(gd);
% delete one edge to leave 1 rectangles and 6 edges
[g,bt2] = csgdel(g,bt, 7);
geometryFromEdges(model, g);

% pdegplot(model,'FaceLabels','on') 
cond = 1; % thermal conductivity
thermalProperties(model,'Face', 1, 'ThermalConductivity', cond);

% pdegplot(model,'EdgeLabels','on','vertexlabels','on')
thermalBC(model,'Edge',3, 'ConvectionCoefficient',1, ...
    'AmbientTemperature',40); 
thermalBC(model,'Edge',4, 'ConvectionCoefficient', 2, ...
    'AmbientTemperature',0); 
generateMesh(model, 'GeometricOrder', 'quadratic'); 
% avoid 'linear', gives large heat imbalance

result = solve(model);

%-----------------------------------------------
% Calculating area or volume-weighted mean temperature
[A, AE] = area(model.Mesh); % AE is row vector

% use the nodal coordinates and the node indices for the corners of each
% element to assemble an array of element corner coordinates.
% Then find the centroid of each element
xyEM = zeros([2, size(result.Mesh.Elements)]); 
% 3D as [xORy, 3 or 6 nodes per element, element_number]
xyEM(:) = result.Mesh.Nodes(:, result.Mesh.Elements);

xyEc = squeeze(mean(xyEM, 2)); % 2*n, centroid (x,y) at element centres

% get temperatures at these centroid coordinates (returns a column vector):
T_nodes = interpolateTemperature(result, xyEc(1,:), xyEc(2,:)); 

% take weighted mean
if contains(model.AnalysisType, 'axisymmetric')
    % NB.  If using an axisymmetric model you cannot replace area() with 
    % volume() because volume() only works with 3D meshes.
    % area() returns the same element areas whether axi-symmetric or not so
    % you need to explicitly calculate element volumes:
    r = xyEc(1,:);  % axis of rotation is always x=0 line
    VE = AE .* (2*pi*r);  % volume of each element
    meanT = sum(T_nodes'.*VE)/sum(VE)
    labels = {'r','z'};
else
    meanT = sum(T_nodes'.*AE)/sum(AE)
    labels = {'x','y'};
end


%------------------------------------------------------
% plot temperature contours
T_levels = logspace(log10(min(result.Temperature)),...
                    log10(max(result.Temperature)), 30);
subplot(1,2,1)
pdeplot(model,'XYData',result.Temperature,'contour','on','levels', T_levels);
xlabel(labels{1})
ylabel(labels{2})
colormap(hsv(40))
title(model.AnalysisType)

% plot upper surface temperature distribution
subplot(1,2,2)
x = 0:0.01:2; 
Tintrp = interpolateTemperature(result,x, zeros(size(x))); 
plot(x, Tintrp)
xlabel(labels{1})
ylabel(['Temperature at ' labels{2} ' = 0'])
%-------------------------------------------------------
% Assess convergence:
Qn = [];
for j = 1:model.Geometry.NumEdges
    Qn(j) = evaluateHeatRate(result,'Edge',j);
    % unlike area(), this seems to know whether axisymmetric or not.
end
disp(['Fractional heat imbalance is ', num2str(2*sum(Qn)/sum(abs(Qn)))])

[darova]

Thank you for the answer

Answer 2

[Roger Moss]

I have now written a function to calculate mean temperature from the PDE solution. For good measure it also gives you volume, mass, centroid, heat capacity and internal energy, for the whole mesh or specified elements. It should work with steady or transient solutions, 2D, 2D-axisymmetric or 3D.

I have only tested it hastily so if you find bugs please let me know.

Enjoy!

Roger