# Crank Nicolson scheme using Finite Volume for 1D Transient Problem.

4 views (last 30 days)
Shahid Hasnain on 11 Sep 2021
Edited: darova on 13 Sep 2021
Dear Community,
I have attached discretized scheme (Problem_1_1 and Problem_1_2). According to complete MATLAB code which is working. I am getting suspicious about results when I compare them with books that are related to the Fully Implicit scheme. A few days earlier I posted my question. Any best idea which can do to satisfy?
Your help will be highly appreciated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
%% Sort out Inputs
Theta = 0.5; % Crank - Nicolson scheme
% Number of Control Volumes
N = 5;
% Domain length
L = 0.02; %[m]
% Grid Spacing (Wall thickness)
h = L/(N);
% Diffusivity (Thermal conductivity)
k = 10; %[W/m-K]
% Specific Heat Capacity
Row_c = 10e6; %[J/m^(3)-K]
% The center of the first cell is at dx/2 & the last cell is at L-dx/2.
x = h/2 : h : L-(h/2);
% Diffusivity (Heat)
alpha = k/Row_c;
% Simulation Time
t_Final = 40;
t_steps = 8;
% Discrete Time Steps
dt = t_Final/(t_steps);
% Left Surface Temperature
dTdt = 0; %[\circ c]
% Right Surface Temperature
T_b = 0; %[\circ c]
% Parameteric Setup
lambda = (alpha.*dt)/(h^2)
ap_0 = Row_c * h/(dt);
dt_Check = Row_c *(h^2)/k;
% Initializing Variable
% C*T^(n+1) = BB*T^(n)
%
% At Node # 1 dpdx=0 is adjusted during simplification of scheme
%
% (1+lamda/2)T_P = (1/2)*lamda*T_E +(T_P)_0
%
% At Node # 5
%
% (1+2*lamda)T_P = lamda*T_W + lamda*T_b+(T_P)_0
%
% At Node Interior Nodes
%
% -(lamda/2)*T_W + (1+lamda)T_P -(lamda/2)*T_E = (lamda/2)*T_W_0 + (1-lamda)T_P_0 +(lamda/2)*T_E_0
%
% Unknowns at time level n
T_Old = zeros(N,1);
% Unknowns at time level n+1
T_New = zeros(N,1);
% Initial Value (Initial Condition)
T_Old(:) = 200;
C = zeros(N,N);
D = zeros(N,1);
for t=1:dt:t_Final
for i = 1: N
if i > 1 && i < N
C(i,i) = (1 + lambda);
BB(i,i) = (1-lambda);
C(i,i+1) = -0.5*lambda;
BB(i,i+1) = 0.5*(lambda);
C(i,i-1) = -0.5*lambda;
BB(i,i-1) = 0.5*(lambda);
D(i)=T_Old(i);
elseif i == 1
% For 1st boundary node
C(i,i) = (1 + 0.5*lambda);
C(i,i+1) = -0.5*lambda;
D(i) = T_Old(i);
else
% For Last boundary node
C(i,i) = (1 + 2*lambda);
C(i,i-1) = -lambda;
BB(i,i) = 0;
D(i) = T_b*lambda + T_Old(i);
end
end
Temp = inv(C)*BB;
T_New = Temp*T_Old;
T_Old = T_New;
end
T_New

### Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!