One Dimensional Laplace Equation with Source term.

Hi, Community,

Need some help to solve 1 Dimensional Laplace equation (Figure attached) by using the Finite-Difference Iterative scheme (I used Jacobi Iterative Method). MATLAB Code is working. When I compare it with analytical, Difference (Error) is significant which is not acceptable. If it is possible, a better idea to improve results, if it is possible? Your idea will be highly appreciated.

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

%

% This code solves the steady 1D diffuion equation using FDM for the

% unknown Temperature. Analytical solution comparison is done with

% Numerical. Example # 2 on Versteeg 2E book Pages 135-139

%

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

clear all

close all

clc

%% Sort out Inputs

tic

% Number of Control Volumes

N = 5;

% Domain length

L = 0.02; %[m]

% Grid Spacing

h = L/(N);

% Diffusivity (Thermal conductivity)

k = 0.5; %[W/m-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]

% Left Boundary Condition

T_a = 100; %[\circ c]

% Right Boundary Condition

T_b = 200; %[\circ c]

% Initialising Temperature

T=zeros(N,1);

% Left Boundary Condition

T(1)= 100; %[\circ c]

% Right Boundary Condition

T(N) = 200; %[\circ c]

% Interior of Domain

index_x = 2:N-1;

%% Formation of Coefficient Matrices, Numerical Solution

% Tolerence

eps = 1.e-6;

% Iteration

It = 0;

T_old = T;

Error = 2*eps;

% Calculation

while (Error > eps)

It = It +1;

T(index_x)= 0.5*(T(index_x-1) + T(index_x+1)+ (q*h*h)/(k));

Error = max(abs(T(:)-T_old(:)));

if any(isnan(T(:))) || any(isinf(T(:)))

fprintf(' Iteration Diverge \n');

return;

end

T_old = T;

end

T_Numeric = T;

%% Analytical Results

X_G = x;

T_Analytical=((T_b-T_a)/L + (q/(2*k))*( L - X_G )).*(X_G) + T_a;

%% Error Analysis

M_Error(N)=0;

for i=1:N

% This x is only use to compare results on Versteeg book

k=find(abs(X_G - x(i))<=1e-3);

% No need to find k (You can but no Need)

% Absolute Error

M_Error(i) = abs(T_Numeric(i)-T_Analytical(i));

% Relative Error

Relative_Error(i) = abs(T_Numeric(i)-T_Analytical(i))/abs(T_Numeric(i));

% Percentage Error

Percentage_Error(i) = Relative_Error(i).*100;

end

%% Post Processing

% Data Table Comparison of Different Results

Nodes = 1 : N;

fprintf('Node Grid_{Loc.} T_{Num.} T_{An.} Error\n\n');

fprintf('====================================================\n');

for i = 1: N;

fprintf('%3.0f %3.5f %3.5f %3.5f %3.5f\n',[Nodes(i); X_G(i); T_Numeric(i); T_Analytical(i); Percentage_Error(i)]);

end

% Graphical Representation

plot(100*x,T_Numeric,'--','MarkerSize',5,'MarkerFaceColor','b');

% hold on

% plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','b');

hold on

plot(100*X_G,T_Analytical,'r','LineWidth',2);

% hold on

% plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','r');

% xlim([0 2]);

% ylim([50 300]);

%Boundary Points Inclusion

xlabel('Distance (m)');

ylabel('Temperature [\circ C]');

legend('T_{Numeric}','T_{Analytical}');

title ('Numerical & Analytical Results Comparison');

grid on;

[Star Strider]

I ran the code just to see what it does and what the problem could be.

Perhaps if you could post ‘Example # 2 on Versteeg 2E book Pages 135-139’ (without violating any copyright) it would be easier to see what the problem is. (I am not aware of that book. Since I do not have it in my library, I cannot refer to it.)

%% Sort out Inputs
tic
% Number of Control Volumes
N = 5;
% Domain length
L = 0.02; %[m]
% Grid Spacing
h = L/(N);
% Diffusivity (Thermal conductivity)
k = 0.5;  %[W/m-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]       
% Left Boundary Condition
T_a = 100;    %[\circ c] 
% Right Boundary Condition
T_b = 200;    %[\circ c] 
% Initialising Temperature
T=zeros(N,1);
% Left Boundary Condition
T(1)= 100;    %[\circ c] 
% Right Boundary Condition
T(N) = 200;    %[\circ c] 
% Interior of Domain
index_x = 2:N-1;
%% Formation of Coefficient Matrices, Numerical Solution
% Tolerence
eps = 1.e-6;
% Iteration 
It = 0;
T_old = T;
Error = 2*eps;
% Calculation
while   (Error > eps)
        It = It +1; 
        T(index_x)= 0.5*(T(index_x-1) + T(index_x+1)+ (q*h*h)/(k));
        Error = max(abs(T(:)-T_old(:)));
     if any(isnan(T(:))) || any(isinf(T(:)))
        fprintf(' Iteration Diverge \n');
        return;
    end
    T_old = T;
end
T_Numeric = T;
%% Analytical Results
X_G = x;
T_Analytical=((T_b-T_a)/L + (q/(2*k))*( L - X_G )).*(X_G) + T_a;
%% Error Analysis
M_Error(N)=0;
 for i=1:N
  % This x is only use to compare results on Versteeg book
  k=find(abs(X_G - x(i))<=1e-3);
  % No need to find k (You can but no Need)
  % Absolute Error
  M_Error(i) = abs(T_Numeric(i)-T_Analytical(i));
  % Relative Error
  Relative_Error(i) = abs(T_Numeric(i)-T_Analytical(i))/abs(T_Numeric(i));
  % Percentage Error
  Percentage_Error(i) = Relative_Error(i).*100;
 end
%% Post Processing
% Data Table Comparison of Different Results
Nodes = 1 : N; 
fprintf('Node  Grid_{Loc.}  T_{Num.}     T_{An.}       Error\n\n');

Node  Grid_{Loc.}  T_{Num.}     T_{An.}       Error

fprintf('====================================================\n');

====================================================

for i = 1: N;
fprintf('%3.0f   %3.5f     %3.5f     %3.5f    %3.5f\n',[Nodes(i); X_G(i); T_Numeric(i); T_Analytical(i); Percentage_Error(i)]);
end

  1   0.00200     100.00000     146.00000    46.00000
  2   0.00600     173.00000     214.00000    23.69942
  3   0.01000     214.00000     250.00000    16.82243
  4   0.01400     223.00000     254.00000    13.90135
  5   0.01800     200.00000     226.00000    13.00000

% Graphical Representation
plot(100*x,T_Numeric,'--','MarkerSize',5,'MarkerFaceColor','b');
% hold on
% plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','b');
hold on
plot(100*X_G,T_Analytical,'r','LineWidth',2);
% hold on 
% plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','r');
% xlim([0  2]);
% ylim([50 300]);
%Boundary Points Inclusion 
xlabel('Distance (m)');
ylabel('Temperature [\circ C]');
legend('T_{Numeric}','T_{Analytical}', 'Location','best');
title ('Numerical & Analytical Results Comparison'); 
grid on;

figure
plot(T_Analytical, T_Numeric, '+')
hold on
p = polyfit(T_Analytical, T_Numeric, 1)

p = 1×2

1.1319 -64.7447

pv = polyval(p, T_Analytical);
plot(T_Analytical, pv, '--k')
hold off
grid

I added the last plot simply to see if there is any consistent relationship.

.

Answer 1

[Wan Ji]

Hi, friend

I have found your boudary conditions are applied at the centre of control volumes of both ends rather than on the boundary of the control volumes viz (x=0, and x=L). You need firstly tell the differences between the control volumes and the edges (here indicate the nodes) that we are interested in.

So you shoud give an array of temperature with size 2*N+1. And the segment size should be h/2. Then it works.

%
% This code solves the steady 1D diffuion equation using FDM for the
% unknown Temperature. Analytical solution comparison is done with
% Numerical. Example # 2 on Versteeg 2E book Pages 135-139
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc
%% Sort out Inputs
tic
% Number of Control Volumes
N = 5;
% Domain length
L = 0.02; %[m]
% Grid Spacing
h = L/(N);
% Diffusivity (Thermal conductivity)
k = 0.5;  %[W/m-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]       
% Left Boundary Condition
T_a = 100;    %[\circ c] 
% Right Boundary Condition
T_b = 200;    %[\circ c] 
% Initialising Temperature
Nh = 2*N + 1; %######## the  number of node with temperature 
T=zeros(Nh,1); %########
% Left Boundary Condition
T(1)= T_a;    %[\circ c] 
% Right Boundary Condition
T(Nh) = T_b;    %[\circ c] 
% Interior of Domain
index_x = 2:Nh-1;
%% Formation of Coefficient Matrices, Numerical Solution
% Tolerence
eps = 1.e-6;
% Iteration 
It = 0;
T_old = T;
Error = 2*eps;
% Calculation
while   (Error > eps)
        It = It +1; 
        T(index_x)= 0.5*(T(index_x-1) + T(index_x+1)+ (q*h*h/4)/(k)); %#### here divided by 4, cause h/2*h/2
        Error = max(abs(T(:)-T_old(:)));
     if any(isnan(T(:))) || any(isinf(T(:)))
        fprintf(' Iteration Diverge \n');
        return;
    end
    T_old = T;
end
T_Numeric = T(2:2:end); %%##### choose the temperature result at the control volume centre
%% Analytical Results
X_G = x;
T_Analytical=((T_b-T_a)/L + (q/(2*k))*( L - X_G )).*(X_G) + T_a;
%% Error Analysis
M_Error(N)=0;
 for i=1:N
  % This x is only use to compare results on Versteeg book
  k=find(abs(X_G - x(i))<=1e-3);
  % No need to find k (You can but no Need)
  % Absolute Error
  M_Error(i) = abs(T_Numeric(i)-T_Analytical(i));
  % Relative Error
  Relative_Error(i) = abs(T_Numeric(i)-T_Analytical(i))/abs(T_Numeric(i));
  % Percentage Error
  Percentage_Error(i) = Relative_Error(i).*100;
 end
%% Post Processing
% Data Table Comparison of Different Results
Nodes = 1 : N; 
fprintf('Node  Grid_{Loc.}  T_{Num.}     T_{An.}       Error\n\n');
fprintf('====================================================\n');
for i = 1: N;
fprintf('%3.0f   %3.5f     %3.5f     %3.5f    %3.5f\n',[Nodes(i); X_G(i); T_Numeric(i); T_Analytical(i); Percentage_Error(i)]);
end
% Graphical Representation
plot(100*x,T_Numeric,'--','MarkerSize',5,'MarkerFaceColor','b');
% hold on
% plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','b');
hold on
plot(100*X_G,T_Analytical,'r','LineWidth',2);
% hold on 
% plot([0, 100.*L], [T_a, T_b], 'd', 'LineWidth',3,'MarkerFaceColor','r');
% xlim([0  2]);
% ylim([50 300]);
%Boundary Points Inclusion 
xlabel('Distance (m)');
ylabel('Temperature [\circ C]');
legend('T_{Numeric}','T_{Analytical}');
title ('Numerical & Analytical Results Comparison'); 
grid on;

The result printed here:

Node  Grid_{Loc.}  T_{Num.}     T_{An.}       Error

====================================================
  1   0.00200     146.00000     146.00000    0.00000
  2   0.00600     213.99999     214.00000    0.00000
  3   0.01000     249.99999     250.00000    0.00000
  4   0.01400     253.99999     254.00000    0.00000
  5   0.01800     226.00000     226.00000    0.00000
  

And you see the plot curves

The two curves are completely coincident.

[Shahid Hasnain]

Thank you, friend. It is great now.