Calculate total heat loss by conduction given temperature and depth profile vectors

Hi all,

I need to solve this problem in Matlab. I have a hot body placed over cold ground, this having a thermal conductivity k = 1.

I need to estimate the total heat loss by conduction into the substratum at each timestep.

For each of these timesteps I have a temperature profile and the associated vertical profile at 1cm intervals (see .mat files attached).

How can I estimate the total heat loss (w/m2) at each time interval (i.e., the heat loss for every profile)?

Any help woud be grately appreciated!

Answer 1

[Torsten]

Here is another way to determine the temporal change of heat content in the substratum.

rho*cp*dT/dt = d/dx (k*dT/dx)

Integrating with respect to x gives

d/dt (integral_{x=-20}^{x=0} (rho*cp*T) dx ) = k*dT/dx @h=0 - k*dT/dx @h=-20

Above, we tried to approximate the right-hand side. Here, we try to approximate the left-hand side.

close all
clc
load("time_in_days.mat")
load("Depth_Profiles.mat")
load("Temperature_Profiles.mat")
num_intervals = size(Temperature_Profiles, 2);
depth_index = find(Depth_Profiles <= -20, 1, 'last');
rhocp = 1.0;    % Specify rho*cp of the substratum [J/(m^3*K)]
for i = 1:num_intervals
    temperature_profile = Temperature_Profiles(:, i);
    energy_content(i) = trapz(Depth_Profiles(depth_index:end),temperature_profile(depth_index:end))*rhocp; %[J/m^2]
end
figure(1)
plot(time_in_days,energy_content) %[J/m^2]

change_of_energy_content = gradient(energy_content)./gradient(time_in_days*24*3600);%[W/m^2]
figure(2)
plot(time_in_days,change_of_energy_content*1e3) %[mW/m^2]

[Simone]

Hi @Torsten, thanks a lot for your time and support.

I have tried to go through your appraoch, but I am still a bit lost.

Just to clarify, if needed, what I am after is to find a single w/m2 value (M_cond) for each curve, that should tell me how much heat I lost (or I am loosing) from the hot body into the substratum at that specific time. For instance:

Now, reading a few books, I found:

As such, if ∆h must be calculated up to a depth where temperature difference applies, I came up with a code (see below). I have also tried to apply the code you kindly provided, as well as a modified version of it (to take into account the 'depth where temperature difference applies').

However, all these results are not consistent, and I don't really know which one is most likely the right output:

Heat Loss by Conduction (Mcond) [Book Approach] after 60 Seconds is                           5.689e+04 W/m²
Heat Loss by Conduction (Mcond) [Book Approach] after 1 Year is                               80.5295   W/m²
Heat Loss by Conduction (Mcond) [Torsten Approach] after 60 Seconds is                        3.373e+04 W/m²
Heat Loss by Conduction (Mcond) [Torsten Approach] after 1 Year is                            161.9705  W/m²
Heat Loss by Conduction (Mcond) [Torsten Approach with variable depth] after 60 Seconds is    1.781e+05 W/m²
Heat Loss by Conduction (Mcond) [Torsten Approach with variable depth] after 1 Year is        524.761   W/m²

Could you kindly provide a bit more guidance?

Thank you very much indeed.

clear 
close all
clc
load("time_in_days.mat")
load("Depth_Profiles.mat")
load("Temperature_Profiles.mat")

T_Ambient = 285;          % Kelvin
T_Substratum = T_Ambient; % Kelvin
k = 1.5;                  % Thermal conductivity in W/m/K
rho = 2700;               % Density in kg/m^3
Cp = 900;                 % Specific heat capacity in J/kg/K
rhocp = rho*Cp;           % Specify rho*cp of the substratum [J/(m^3*K)]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% Method 1 (extracted from book [see attached figure]) %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

First_Curve_Affected_Depth = find(Temperature_Profiles(:,1)>T_Substratum + 0.01); 
First_Curve_Affected_Depth = First_Curve_Affected_Depth(1);
Last_Curve_Affected_Depth = find(Temperature_Profiles(:,end)>T_Substratum + 0.01); 
Last_Curve_Affected_Depth = Last_Curve_Affected_Depth(1);
Qcond_First_Curve = k * (...
    (Temperature_Profiles(end,1) - Temperature_Profiles(First_Curve_Affected_Depth,1))...
    /...
    ((Depth_Profiles(end) - Depth_Profiles(First_Curve_Affected_Depth)))...
    );
Qcond_Last_Curve = k * (...
    (Temperature_Profiles(end,1) - Temperature_Profiles(Last_Curve_Affected_Depth,1))...
    /...
    ((Depth_Profiles(end) - Depth_Profiles(Last_Curve_Affected_Depth)))...
    );
disp(strcat(['Heat Loss by Conduction (Mcond) [Book Approach] after ', ...
    num2str(round(seconds(days(time_in_days(1))))), ' Seconds is                           ', ...
    sprintf('%.3e', Qcond_First_Curve), ' W/m²']))
disp(strcat(['Heat Loss by Conduction (Mcond) [Book Approach] after ', ...
    num2str(round(years(days(time_in_days(end))))), ' Year is                               ', ...
    num2str(Qcond_Last_Curve), '   W/m²']))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% Method 2 (Torsten) %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load("time_in_days.mat")
load("Depth_Profiles.mat")
load("Temperature_Profiles.mat")

num_intervals = size(Temperature_Profiles, 2);
depth_index = find(Depth_Profiles <= -20, 1, 'last');
for i = 1:num_intervals
    temperature_profile = Temperature_Profiles(:, i);
    energy_content(i) = trapz(Depth_Profiles(depth_index:end),temperature_profile(depth_index:end))*rhocp; %[J/m^2]
end
change_of_energy_content = gradient(energy_content)./gradient(time_in_days*24*3600);%[W/m^2]
disp(strcat(['Heat Loss by Conduction (Mcond) [Torsten Approach] after ', ...
    num2str(round(seconds(days(time_in_days(1))))), ' Seconds is                        ', ...
    sprintf('%.3e', change_of_energy_content(1)), ' W/m²']))
disp(strcat(['Heat Loss by Conduction (Mcond) [Torsten Approach] after ', ...
    num2str(round(years(days(time_in_days(end))))), ' Year is                            ', ...
    num2str(change_of_energy_content(end)), '  W/m²']))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% Method 3 (Torsten with variable depth) %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load("time_in_days.mat")
load("Depth_Profiles.mat")
load("Temperature_Profiles.mat")

rho = 2700;                % Density in kg/m^3
Cp = 900;                  % Specific heat capacity in J/kg/K
rhocp = rho*Cp;            % Specify rho*cp of the substratum [J/(m^3*K)]
num_intervals = size(Temperature_Profiles, 2);
for i = 1:num_intervals
    temperature_profile = Temperature_Profiles(:, i);
    Depth_Affected_Idx = find(abs(temperature_profile - (T_Ambient)) <= 1e-2); % Add tolerance to avoid numerical artefacts moving affected depth far behind i.e. 285.00000000001    
    Depth_Affected_Idx = Depth_Affected_Idx(end)+1;    
    %     depth_index = find(Depth_Profiles <= -20, 1, 'last');
    energy_content(i) = trapz(Depth_Profiles(Depth_Affected_Idx:end),temperature_profile(Depth_Affected_Idx:end))*rhocp; %[J/m^2]
end
change_of_energy_content = gradient(energy_content)./gradient(time_in_days*24*3600);%[W/m^2]
disp(strcat(['Heat Loss by Conduction (Mcond) [Torsten Approach with variable depth] after ', ...
    num2str(round(seconds(days(time_in_days(1))))), ' Seconds is    ', ...
    sprintf('%.3e', change_of_energy_content(1)), ' W/m²']))
disp(strcat(['Heat Loss by Conduction (Mcond) [Torsten Approach with variable depth] after ', ...
    num2str(round(years(days(time_in_days(end))))), ' Year is        ', ...
    num2str(change_of_energy_content(end)), '   W/m²']))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% P L O T %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0.01 0.05 0.99 0.95]);
hold on
plot(Temperature_Profiles(:,1),Depth_Profiles,'b-',LineWidth=2)
plot(Temperature_Profiles(First_Curve_Affected_Depth,1),...
    Depth_Profiles(First_Curve_Affected_Depth),'bo',MarkerSize=10,...
    LineWidth=2)
plot(Temperature_Profiles(:,end),Depth_Profiles,'r-',LineWidth=2)
plot(Temperature_Profiles(Last_Curve_Affected_Depth,end),...
    Depth_Profiles(Last_Curve_Affected_Depth),'ro',MarkerSize=10,...
    LineWidth=2)
set(gca,'YLim',[-25 0])
xlabel('Temperature (K)', 'FontSize', 20, Interpreter='latex');
ylabel('Depth (m b.g.l.)', 'FontSize', 20, Interpreter='latex');
legend(...
    strcat(['Vertical T Profile After: ',...
    num2str(round(seconds(days(time_in_days(1))))),' Seconds']),...
    strcat(['Equilibrium with T Substratum (+0.01) reached at: ',...
    num2str(Depth_Profiles(First_Curve_Affected_Depth)),' Meters']),...
    strcat(['Vertical T Profile After: ',...
    num2str(round(years(days(time_in_days(end))))),' Year']),...
    strcat(['Equilibrium with T Substratum (+0.01) reached at: ',...
    num2str(Depth_Profiles(Last_Curve_Affected_Depth)),' Meters']),...
    'location','SE','interpreter','latex','fontsize',12)
axis square
grid on
set(gca,'box','on','Color', 'white', 'XColor', 'black', 'YColor', 'black');

text(0.25, 0.55, ['$M_{\mathrm{cond}}$ after ', ...
    num2str(round(seconds(days(time_in_days(1))))), ' Seconds is ', ...
    sprintf('%.3e', Qcond_First_Curve), ' W/m$^2$'], 'Units', 'normalized',...
    'Interpreter', 'latex', 'FontSize', 16, 'Color', 'blue');


text(0.25, 0.45, ['$M_{\mathrm{cond}}$ after ', ...
    num2str(round(years(days(time_in_days(end))))), ' Year is ', ...
    num2str(Qcond_Last_Curve), ' W/m$^2$'], 'Units', 'normalized',...
    'Interpreter', 'latex', 'FontSize', 16, 'Color', 'red');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Torsten]

I gave you two ways to compute the heat input into the substratum. One way is to compute the heat flux over the boundaries of your domain, the other is to compute the temporal change of energy content within the domain. What specific question do you have ?

The book's approach is the same as mine - the temperature difference is

deltaT/deltah = (temperature_profile(end) - temperature_profile(end-1)) / 0.01;

[Simone]

Hi @Torsten, I wonder which is the "best" way of calculating it given my specific problem. I mean, all I need to know is the total and "istantaneous" heat (W/m2) lost into the substratum at a given specific time (i.e., t1). Altough you kindly provided 2 different approaches, and I "found" another one by myself, these outputs are not consistent between them. So the final question is: which approach yelds more realistic W/m2 values?

[Simone]

The book's approach is the same as mine - the temperature difference is

deltaT/deltah = (temperature_profile(end) - temperature_profile(end-1)) / 0.01;

If so, why do they yeld different results? Can you spot anything wrong with my code?

Thanks a lot for your patience

[Torsten]

These are my two suggestions:

load("time_in_days.mat")
load("Depth_Profiles.mat")
load("Temperature_Profiles.mat")
num_intervals = size(Temperature_Profiles, 2);
depth_index = 1;%find(Depth_Profiles <= -20, 1, 'last');
rhocp = 2700*900;    % Specify rho*cp of the substratum [J/(m^3*K)]
for i = 1:num_intervals
    temperature_profile = Temperature_Profiles(:, i);
    energy_content(i) = trapz(Depth_Profiles(depth_index:end),temperature_profile(depth_index:end))*rhocp; %[J/m^2]
end
%figure(1)
%plot(time_in_days,energy_content) %[J/m^2]
change_of_energy_content = gradient(energy_content)./gradient(time_in_days*24*3600);%[W/m^2]
figure(1)
plot(time_in_days(1:10),change_of_energy_content(1:10)) %[W/m^2]

change_of_energy_content(10)

ans = 1.5132e+03

k = 1; % thermal conductivity
num_intervals = size(Temperature_Profiles, 2);
total_heat_loss = zeros(1, num_intervals);
for i = 1:num_intervals
    temperature_profile = Temperature_Profiles(:, i);
    depth_index = 1;%find(Depth_Profiles <= -20, 1, 'last');
    dT_dh_surface = (temperature_profile(end) - temperature_profile(end-1)) / 0.01;
    dT_dh_Depth = (temperature_profile(1+depth_index) - temperature_profile(1+depth_index - 1)) / 0.01;
    total_heat_loss(i) = k * (dT_dh_surface - dT_dh_Depth);
end
figure(2)
plot(time_in_days(1:10),total_heat_loss(1:10))

total_heat_loss(10)

ans = 1.6875e+03

[Torsten]

If so, why do they yeld different results? Can you spot anything wrong with my code?

You have to fix a domain over the course of your computation and compute the gradients at its boundaries for all times. The difference will give you the inflow of energy into this part of your domain. I chose the complete depth, you may want to choose a smaller part. But it must remain fixed - otherwise you cannot compare the results.

[Simone]

Thanks a lot @Torsten for your step by step eplanation and support, it was superb.

[Torsten]

Here is a second-order approximation of the temperature gradients at the boundaries, but it doesn't seem to change that much:

dT_dh_surface = (1.5*temperature_profile(end) - 2*temperature_profile(end-1) + 0.5*temperature_profile(end-2)) / 0.01;
dT_dh_Depth = (-1.5*temperature_profile(2+depth_index) + 2*temperature_profile(1+depth_index) - 0.5*temperature_profile(1+depth_index - 1)) / 0.01;

For reference:

https://en.wikipedia.org/wiki/Finite_difference_coefficient

Answer 2

[Torsten]

How can I estimate the total heat loss (w/m2) at each time interval (i.e., the heat loss for every profile)?

Total heat flow into the domain is (k*dT/dh @h=0) - (k*dT/dh @h=-20). Maybe you only want to evaluate this at one end.

To estimate the gradient, use a finite difference approximation.

[Simone]

Hi @Torsten, thanks for getting back!

Maybe you only want to evaluate this at one end. Sorry, can you provide more details on this? What does it mean?

To estimate the gradient, use a finite difference approximation: I came up with this code, but the output does not make any sense... :

close all
clc
k = 1; % thermal conductivity
num_intervals = size(Temperature_Profiles, 2);
total_heat_loss = zeros(1, num_intervals);
for i = 1:num_intervals
    temperature_profile = Temperature_Profiles(:, i);
    depth_index = find(Depth_Profiles <= -20, 1, 'last');
    dT_dh_surface = (temperature_profile(2) - temperature_profile(1)) / 0.01;
    dT_dh_Depth = (temperature_profile(depth_index) - temperature_profile(depth_index - 1)) / 0.01;
    total_heat_loss(i) = k * (dT_dh_surface - dT_dh_Depth);
end
figure
plot(time_in_days,total_heat_loss)
set(gca, 'yScale', 'log');
set(gca, 'xScale', 'log');

% I have attached time_in_days vector

Could you kindly provide further assistance?

[Torsten]

Maybe you only want to evaluate this at one end. Sorry, can you provide more details on this? What does it mean?

It means that at h = 0, you have a heat flow into the domain and at h = -20, you have a heat flow out of the domain. Maybe you only want to compute the flow into the domain (k*dT/dh @h=0) or only the flow out of the domain (-k*dT/dh @h=-20) and not the net flow (k*dT/dh @h=0) - (k*dT/dh @h=-20).

To estimate the gradient, use a finite difference approximation: I came up with this code, but the output does not make any sense... :

You could use difference quotients of higher order that involve more h-points than only two to approximate the temperature gradient. But if it's correct that the distance between two adjacent h-points is 0.01 m, your code looks fine in principle. But you have to use

dT_dh_surface = (temperature_profile(end) - temperature_profile(end-1)) / 0.01;

instead of

dT_dh_surface = (temperature_profile(2) - temperature_profile(1)) / 0.01;