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
