1 D Transient by FVM Explicit (Working Code) but Implicit Issue with Time.

Hi, Community,

Need some help to solve 1 D Unsteady Diffusion Equation by Finite Volume (Fully Implicit) Scheme . MATLAB Code is working for

a- Explicit FVM (I compared results with attached book)

b- Finite Volume (Fully Implicit) Scheme is also working but results are different when I compared with book results (Pages 266 - 269). Also I put t_Final=2 to compare my coefficient matrix with book which prefectly match. So when I change t_Final=40. Reults are different from book.

Am I missing something? My question is to improve my results in case of Finite Volume (Fully Implicit) Scheme. Thank you in advance. Your idea will be highly appreciated.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% This code solves the unsteady 1D diffuion equation using FVM for the

% unknown Temperature.

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

close all

clc

%% Sort out Inputs

% 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]

% Uniform Heat Generation, A Source Term

q = 1000e3; % [W/m^(3)]

% 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);

% Cross sectional area of the 1D domain

A = 1; %[m^2]

% Diffusivity (Heat)

alpha = k/Row_c;

% Simulation Time

t_Final = 40;

% Discrete Time Steps

dt = 2;

% Left Surface Temperature

T_a = 0; %[\circ c]

% Right Surface Temperature

T_b = 0; %[\circ c]

% Parameteric Setup

lambda = (alpha.*dt)/(h^2);

%% Initializing Variable

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

%% Formation of Coefficient Matrices

% Time loop

for n = 1:t_Final/dt;

% First Node

T_New(1) = lambda*T_a + (1-lambda)*T_Old(1) + lambda*T_Old(2);

% Interior Nodes

for i = 2: N-1

T_New(i) = lambda*T_Old(i-1) + (1-2*lambda)*T_Old(i) + lambda*T_Old(i+1);

end

% Last Node

T_New(N) = lambda*T_Old(N-1) + (1-3*lambda)*T_Old(N) + 2*lambda*T_b;

% Update calculations

T_Old = T_New;

end

T_New

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Fully Implicit Working Code

%%%%%%%%%%%%%%%%%%%%%%%%%

%% Sort out Inputs

% 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]

% Uniform Heat Generation, A Source Term

q = 1000e3; % [W/m^(3)]

% 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);

% Cross sectional area of the 1D domain

A = 1; %[m^2]

% Diffusivity (Heat)

alpha = k/Row_c;

% Simulation Time

t_Final = 2;

% Discrete Time Steps

dt = 2;

% Left Surface Temperature

T_a = 0; %[\circ c]

% Right Surface Temperature

T_b = 0; %[\circ c]

% Parameteric Setup

lambda = (alpha.*dt)/(h^2);

%% Initializing Variable

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

% % Time loop

for n = 1:t_Final/dt

C = zeros(N,N);

D = zeros(N,1);

for i = 1: N

if i > 1 && i < N

C(i,i) = (1 + 2*lambda).*T_Old(i);

C(i,i+1) = -lambda.*T_Old(i);

C(i,i-1) = -lambda.*T_Old(i);

D(i)=T_Old(i);

elseif i == 1

% For 1st boundary node

C(i,i) = (1 + lambda).*T_Old(i);

C(i,i+1) = -lambda.*T_Old(i);

D(i) = T_Old(i);

else

% For Last boundary node

C(i,i) = (1 + 3*lambda).*T_Old(i);

C(i,i-1) = -lambda.*T_Old(i);

D(i) = 2*T_b*lambda + T_Old(i);

end

end

T_New = (inv(C)*D);

% Update calculations

T_Old = T_New

end

C