One Dimensional Laplace Equation with Source term.

18 views (last 30 days)
Shahid Hasnain on 29 Aug 2021
Commented: Shahid Hasnain on 29 Aug 2021
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 on 29 Aug 2021
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.
.

Wan Ji on 29 Aug 2021
Edited: Wan Ji on 29 Aug 2021
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 on 29 Aug 2021
Thank you, friend. It is great now.