1-D Transient Conduction - Implicit Method

Hi,

I am writing this code which is intended to find the minimum required depth of water pipes required for them to not reach 0 degrees after 60 days, using the finite-difference method, with a soil temperature starting at 20 degrees and air temperature of -15 degrees.

While the code is not resulting in errors, the Temperature results received are not what is expected as they are expected to increase gradually to 20 degrees, as depth increases and decrease to -15 degrees as time increases.

Hoping someone with more experience can help me out!

clear all
p=2050;     % density(kg/m^3)
k=0.52;     % conductivity(W/m-K)
c=1840;     % specific heat capacity(J/kg-K)
T_ini=20;    % initial temperature(K)
T_inf=-15;    % external temperature(K)
alpha=k/(p*c);      % Let thermal diffusivity be alpha
M=30;      % number of time steps
N=30;
t=60*60*60*24; % time (s)
DELTA_t=t/M;  % time step duration(s)
time=1:DELTA_t:t;

L=6;
DELTA_x=L/N;

Ti=zeros(N,M);

Ti(1,1:M)=T_ini;
Ti(2:N-1,1)=T_inf;


a = zeros(N,1);
b = zeros(N,1);
c = zeros(N,1);
d = zeros(N,1);

a(1)= 1;
b(1)= 0;
c(1)= 0;
d(1)= T_ini;
a(N)= 1;
b(N)= 0;
c(N)= 0;
d(N)= T_inf;

Ti(1,1:M)=T_inf;
Ti(2:N,1:M)=T_ini;

for i=2:N-1
    for j=1:M-1
        a(i)=2*((alpha*DELTA_t)/DELTA_x^2)+1;
        b(i)=(alpha*DELTA_t)/DELTA_x^2;
        c(i)=b(i);
        d(i)=-c(i)*Ti(i-1,j+1)+a(i)*Ti(i,j+1)-b(i)*Ti(i+1,j+1); 
    end
end
% Creating the tri-diagonal matrix
Ai=zeros(N);
for i=1:N
    for j=1:N
        if i==j+1
            Ai(i,j)=-b(i);
        elseif i==j
            Ai(i,j)=a(i);
        elseif i==j-1
            Ai(i,j)=-c(i);
        end
    end
end

Ai_inv=inv(Ai);
for i=2:N-1
    for j=1:M-1
        Ti(i,j+1)=(-c(i)*Ti(i-1,j+1)-b(i)*Ti(i+1,j+1))/a(i);
        Ti(i:N,1)=d(i);
        Ti(i:N,j+1)=Ti(i:N,j).*(Ai_inv(i:N,j));
    end
end

Z=N*(0.1/L);

plot(Ti(Z,:),time, 'b',Ti(N,:),time, 'r')

xlabel('Temperature')
ylabel('Time')
grid on
set(gca,'fontsize',16)
legend('0.1','xm')

Answer 1

[Ravi]

Hi Gabriel Quattromani,

I understand that you are facing some issues in implementing the code for transient conduction.

Please find the below modifications.

• Initial conditions:

The initial temperature profile should be “T_ini” throughout the soil depth, except at the surface which is exposed to the air temperature “T_inf”. The original code incorrectly set the initial temperatures for the interior nodes to “T_inf”.

• Boundary conditions:

The original code incorrectly set the boundary condition at the bottom of the domain to “T_inf”. The modification ensures that the bottom boundary condition is insulated or set to the initial soil temperature “T_ini”.

• Matrix assembly and time-stepping:

The tridiagonal matrix should be assembled once outside the time loop, and at each time step, the system “Ai * Ti(:,j+1) = d” should be solved to update the temperature profile. The modified code correctly assembles the tridiagonal matrix and performs the time-stepping loop to solve for the temperature at each time step.

Here is the complete code.

clear all 
p=2050;     % density(kg/m^3) 
k=0.52;     % conductivity(W/m-K) 
c=1840;     % specific heat capacity(J/kg-K) 
T_ini=20;    % initial temperature(C) 
T_inf=-15;    % external temperature(C) 
alpha=k/(p*c);      % thermal diffusivity 
M=30;      % number of time steps 
N=30;      % number of spatial steps 
t=60*60*60*24; % time (s) 
DELTA_t=t/M;  % time step duration(s) 
L=6;         % depth (m) 
DELTA_x=L/N; % spatial step size(m) 

% Initialize temperature matrix 
Ti=zeros(N+1,M+1); 
Ti(:,1)=T_ini;  % Initial temperature profile 
Ti(1,:)=T_inf;  % Surface temperature 

% Coefficient matrix 
a = (1 + 2*alpha*DELTA_t/DELTA_x^2) * ones(N+1, 1); 
b = (-alpha*DELTA_t/DELTA_x^2) * ones(N+1, 1); 
c = b; 

% Modify coefficients for boundary conditions 
a(1) = 1; 
b(1) = 0; 
c(1) = 0; 
a(N+1) = 1; 
b(N+1) = 0; 
c(N+1) = 0; 

% Assemble the tridiagonal matrix 
Ai = diag(a) + diag(b(2:end), -1) + diag(c(1:end-1), 1); 

% Time-stepping 
for j=1:M 
    d = Ti(:,j);  % Right-hand side 
    d(1) = T_inf;  % Boundary condition at the surface 
    d(N+1) = T_ini;  % Boundary condition at the bottom 
    Ti(:,j+1) = Ai\d;  % Solve the system for the next time step 
end 

% Plotting 
depths = (0:N)*DELTA_x; 
times = (0:M)*DELTA_t; 
% Plot temperature versus depth for the last time step 
figure; 
plot(depths, Ti(:,end), 'b'); 
xlabel('Depth (m)'); 
ylabel('Temperature (C)'); 
grid on; 
set(gca,'fontsize',16); 
title('Temperature profile after 60 days'); 

Hope this solution helps.

Thanks,

Ravi Chandra