1-D Transient Conduction - Implicit Method

51 views (last 30 days)
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')

Accepted Answer

Ravi
Ravi on 28 Dec 2023
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

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!