1D HEAT CONDUCTION

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;

[Gwendolyn]

@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;

Answer 1

[Torsten]

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).