index calculations

Hi, I need to do index calculations. I've a vertical column divided into thin layers, for which properties are considered constant. The temperature T(i + 1) and the heat flow Q(i + 1), at the bottom of each layer, are determined from the temperature T(i) and the heat flow Q(i) at the top of each layer by

T(i+1) = T(i) + (Q(i).*Dz(i))./k(i) - (A(i).*Dz(i).^2)./2.*k(i)

and

Q(i+1) = Q(i) - A(i).*Dz(i)

where A(i) and k(i) are known for each layer and Dz(i) is the layer thickness and set to 1. Somebody can suggests me how implement in matlab this operation? thanks Gianluca

[Jan]

This looks like valid Matlab already. What exactly is the problem?

[gianluca]

Hi Jan,

I would implement a .m file to compute one-dimensional steady-state conductive geotherm using a boot-strapping method (Chapman, 2011), which requires thermal conductivity (k), internal heat source (A), surface temperature (T at i=1) and surface heat flow (Q at i = 1) as input. I think of having to implement a for command but I do not know how link the expression of T and Q and compute the temperatures at each interval.

Answer 1

[Jan]

%define depth range and intervals [m]
z = (0:100:10000)';   % No addition [ and ] around vectors

%define variables
% k = zeros(size(z));   % No pre-allocation required
% A = zeros(size(z));   % No pre-allocation required
Q = zeros(size(z));
T = zeros(size(z));

To = 18;        %surface temperature [°C]
Qo = 0.07;      %surface heat flow [W/m2]

k = 2 + 0.00005 .* z;             % Spaces around operators increase readability
A = 2.5e-6 .* exp(-0.0001 .* z);  % 2.5e-6 is faster than 2.5*10^-6

diffz = diff(z);
Q(1) = Qo;                        % Once outside the loop  
for i = 1:(length(z)-1)
  Q(i+1) = Q(i) - A(i) .* diffz(i);
end
% Faster, but needs more memory:
% Q = cumsum([Qo; -A(1:end-1) .* diffz])  

T(1) = To;                         % Once outside the loop
for i = 1:(length(z)-1)
  T(i+1) = T(i) + (Q(i) .* diffz(i)) ./ ...
           k(i) - (A(i) .* diffz(i) .^ 2) ./ 2 .* k(i);
end
plot(z,T)

The spaces around the operators are not required, but this can help to avoid ambiguities, compare:

2.*3
2. *3
2 .* 3

[gianluca]

thanks Jan, your advices will be useful also in the future.

Answer 2

[gianluca]

This script, to compute one-dimensional steady-state conductive geotherm with internal heat source, works. Any advice to improve it is welcome. Thanks in advance

%define depth range and intervals [m]
 z = [0:100:10000]';

%define variables
k = zeros(size(z));
A = zeros(size(z));
Q = zeros(size(z));
T = zeros(size(z));

%define boundary conditions
To = 18;        %surface temperature [°C]
Qo = 0.07;      %surface heat flow [W/m2]

%compute thermal conductivity [W/(m K)] and heat production [W/m3] vertical profile
k = 2 + 0.00005.*z;       %e.g. thermal conductivity increases linearly with depth
A = 2.5.*10.^-6.*exp(-0.0001.*z);       %e.g. heat production decreases exponentially with depth

%compute the heat flow Q at the bottom of each interval
for i = 1:(length(z)-1)
Q(1) = Qo;
Q(i+1) = Q(i) - A(i).*(z(i+1) - z(i));
end

%compute geotherm
for i = 1:(length(z)-1)
T(1) = To;
T(i+1) = T(i) + (Q(i).*(z(i+1) - z(i)))./k(i) - (A(i).*(z(i+1) - z(i)).^2)./2.*k(i);
end
plot(z,T)