1-D Transient Heat Conduction With No Heat Generation [FDM] [Crank-Nicholson]

3 views (last 30 days)
I am currently coding a 1 D Transient Heat Conduction using Crank Nicholson method and I would like an expert opinion as to the accuracy of the result
%1-D Transient Heat Conduction With No Heat Generation [FDM][CN]
clear all
clc
clf
%% Variable Declaration
n = 21; %number of nodes
L = 1; %length of domain
A = sparse(n,n); %initilizing Space
B = sparse(n,n);
Bx = zeros(n,1);
Ta = zeros(n,1);
Tb = zeros(n,1);
dx = L/(n-1); %domain element
dt = 20; %time step
tmax = 4000; %total Time steps (s)
t = 0:dt:tmax;
Tl = 100; %temperature at left face °C
Tr = 100; %temperature at right face °C
x = linspace(0,L,n); %linearly spaced vectors x direction
alpha = 1e-4; %thermal diffusivity (m^2/s)
r = alpha * dt /(2 *dx^2); %for stability, must be 0.5 or less
%% Set Up Matrix
for i = 2 : n-1
A(i,i-1) = -r;
A(i,i ) = (1+2*r); %Implicit Matrix
A(i,i+1) = -r;
end
A(1 ,1 ) = 1;
A(n ,n ) = 1;
for j = 2 : n-1
B(j,j-1) = r;
B(j,j ) = (1-2*r); %Explicit Matrix
B(j,j+1) = r;
end
B(1 ,1 ) = 1;
B(n ,n ) = 1;
%% Boundry Condition
Bx(1,1) = Tl; %Left Wall (Dirichlet conditions)
Bx(n,1) = Tr; %Right Wall(Dirichlet conditions)
%% Solution
for k = 2 : length(t) %time steps
Bx(2 : n-1) = Tb(2 : n-1);
Tb = B * Bx;
fprintf('Time t=%d\n',k-1);
end
Ta = A\Tb; %Solve CN Matrix
%% Plot
pos1 = [0.1 0.17 0.4 0.7];
subplot('Position',pos1)
plot(x,Ta);
title('Nodes Vs Temperature');
xlabel('Nodes (i)');
ylabel('Temperature °C');
grid on
grid minor
a=0:L;
d=0:0;
pos2 = [0.51 0.45 0.45 0.1];
subplot('Position',pos2)
ax = gca;
imagesc(a,d,Ta', 'Parent', ax);
ax.YAxis.Visible = 'off';
title('Temperature Gradient');
xlabel('Nodes (i)');
h = colorbar;
ylabel(h, 'Temperature °C')
colormap jet
  4 Comments
Leroy Coelho
Leroy Coelho on 3 Jul 2020
-r\cdot{T}_{i-1,j+1} + (1+2\cdot{r})\cdot{T}_{i,j+1} -r\cdot{T}_{i+1,j+1}=r\cdot{T}_{i-1,j} + (1-2\cdot{r})\cdot{T}_{i,j} +r\cdot{T}_{i+1,j}
r= \frac {alpha\cdot{dt}}{2\cdot{dx^2}}

Sign in to comment.

Accepted Answer

Leroy Coelho
Leroy Coelho on 7 Jul 2020
Code Update and cross checked with other sources for accuracy
%% Solution
for k = 2 : length(t) %time steps
Tb = B * Bx;
Ta = A\Tb; %Solve CN Matrix
Bx(2 : n-1) = Ta(2 : n-1);
fprintf('Time t=%d\n',k-1);
end

More Answers (1)

Bjorn Gustavsson
Bjorn Gustavsson on 3 Jul 2020
My suggestion is that you add a line:
Tall(:,k) = Tb;
at the end of your "time-step"-loop. Then after you're done with it all you can simply calculate the time-derivative and second-order spatial derivative of your solution. Then you can check how well those satisfy the heat-equation by simply calculating the required gradients:
subplot(2,2,1)
imagesc(t,x,Tall)
colorbar
subplot(2,2,2)
imagesc(t,x,dT1)
caxis([-1 1]/100)
caxis([0 .2])
colorbar
[d2T21,d2T22] = gradient(dT2,t,x);
subplot(2,2,3)
imagesc(t,x,d2T22)
colorbar
subplot(2,2,4)
imagesc(t,x,(dT1-alpha*d2T22)./Tall)
caxis([-1 1]/100)
colorbar
Here I've not bothered with calculating the gradients properly centred for your comparison - but that surely is your job...
HTH

Tags

Community Treasure Hunt

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

Start Hunting!