1D HEAT CONDUCTION

12 views (last 30 days)
Giovanni Karahoxha
Giovanni Karahoxha on 5 Dec 2023
Edited: Gwendolyn on 8 Dec 2023
I succesfully created the following script for 1D conduction using the explicit discretization thanks to the natural way to write in on paper using matrix operations. I'm trying to write down (using pen and paper!) the matrix operations for the implicit discretization (the problem is the same for Crank Nicolson)... Obviously a change on the time step doesn't work... e.g. fourier * (T(i-1, k+1) -...)
%% Finite differences - Explicit discretization
close all; clear; clc;
%% Problem's input
alpha = input('Thermal diffusivity [mc/s] = ');
L = input('Lenght of the bar [m] = ');
Tleft = input('Imposed temperature on the left [K] = ');
Tright = input('Imposed temperature on the right [K] = ');
T0 = input('Temperature before the perturbation [K] = ');
m = input('How many nodes for along x? ');
n = input('How many time steps? ');
tmax = input('How long do you want to see? ');
%deltat = tmax / n;
deltax = L / m;
%fourier = alpha*deltat/deltax^2 % Fron page 222 Mills: Fourier has to be minor than 0.5 to avoid divergent oscillation
% That means 0.5*deltax^2/alfa = deltat to avoid divergent oscillation
% Fix fourier 0.4
fourier = 0.4;
deltat = ( fourier * deltax^2 ) / alpha;
% At the beginning, we have a matrix of the temperature with the quite
% temperature everywhere along x!
T = ones(m, n+1)*T0;
% Live plot
figure;
h = mesh(linspace(0, L, m), linspace(0, tmax, n+1), T');
xlabel('Position [m]');
ylabel('Time [s]');
zlabel('Temperature [K]');
title('Temperature in the bar',LineWidth=19);
%% Iterative calculations
for k = 1:n
% For step k, it is calculated the temperature at node i for time step
% k+1
for i = 2:m-1
T(i, k+1) = T(i, k) + fourier * (T(i-1, k) - 2*T(i, k) + T(i+1, k));
end
% Border conditions
T(1, k+1) = Tleft;
T(end, k+1) = Tright;
% Updated mesh plot
set(h, 'ZData', T');
title(['Temperature distribution at time t = ', num2str(k * deltat)]);
drawnow;
end
%% Final plot
figure(1)
plot(linspace(0, L, m), T(:, end));
xlabel('Position [m]');
ylabel('Temperature [K]');
title(['Temperature at time t = ', num2str(tmax,2)]);
pause;
  1 Comment
Gwendolyn
Gwendolyn on 7 Dec 2023
Edited: Gwendolyn on 8 Dec 2023
@basket randomBoundary conditions need to be incorporated into the system of equations before solving.
# Modify A matrix for boundary conditions
A[1, 1] = 1;
A[m, m] = 1;
# Modify the right-hand side vector for boundary conditions
b = T_prev + fourier * np.concatenate((Tleft, np.zeros(m-2), Tright))
b[1] += fourier * Tleft;
b[m] += fourier * Tright;

Sign in to comment.

Answers (1)

Torsten
Torsten on 5 Dec 2023
Moved: Torsten on 5 Dec 2023
In each time step, you have to solve the linear system of equations
T(1,k+1) = Tleft
T(i, k+1) - deltat/deltax^2*alpha* (T(i-1, k+1) - 2*T(i, k+1) + T(i+1, k+1)) = T(i, k)
T(end,k+1) = Tright
Put this in matrix form
A * T^(k+1) = b(T^k)
and solve for T^(k+1) given T^(k).

Categories

Find more on Biotech and Pharmaceutical in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!