## Welcome to the World of Modelling and Simulation

### What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

### 2D Laplace Equation Solution by 5 Point Finite Difference Approximation

The temperature distribution of a rectangular plate is described by the following two dimensional (2D) Laplace equation:

Txx + Tyy = 0

The width (w), height (h), and thickness (t) of the plate are 10, 15, 1 cm, respectively. There is a heat source at the top edge, which is described as, T = 100 sin (πx / w) Celsius, and all other three edges are kept at 00 C. There is no heat flow along the plate thickness since the faces are insulated. Also, assume that there is no internal heat generation so that the problem becomes 2D temperature distribution.

Solve the Laplace equation with a 5 point finite difference approximation by MATLAB where Δx = Δy = 2.5 cm. The temperature distribution is defined by T(x,y)

The following MATLAB program develops a 5-point approximation stencil to solve the 2D Laplace equation.

MATLAB Program:

close all;
clc;
% Width, height & thickness of the thin plate (dimension, cm)
W=10.0;
h=15.0;
t=1.0;
% Gird size definition along x and y axes
xx=2.5;
yy=2.5;

imax=(W/xx)-1; % Array dimension along the width
jmax=(h/yy)-1; % Array dimension along the height
N_array=imax*jmax; % Total number of grid points

% Generating the Matrix, A
for j=1:N_array
for i=1:N_array
if i==j
a(i,j)=-4.0;
end
if j==(i+1)
a(i,j)=1.0;
end
if j==(i-1)
a(i,j)=1.0;
end
if j==(i+5)
a(i,j)=1.0;
end
if j==(i-5)
a(i,j)=1.0;
end
if i==(imax+2)
a(i,i+1)=0.0;
end
if i==(imax+3)
a(i,i-1)=0.0;
end
if i==(imax+2+imax+2)
a(i,i+1)=0.0;
end
if i==(imax+2+imax+2+1)
a(i,i-1)=0.0;
end
end
end

for i=1:imax
x(i)=(i)*xx;
Tb(i)=100*sin(pi*x(i)/W); % Top edge temperature (in Celcius)
end

% Generating the Matrix, b
k=imax+2;
for i=1:imax
b(k)=-Tb(i);
k=k+imax+2;
end

%Displaying A and b matrices
B=b';
T=GaussElimination(a,B); % Calling the function defined elsewhere
TT(:,1) = T(1:5);
TT(:,2) = T(6:10);
TT(:,3) = T(11:15);
TT = [TT;Tb];

imagesc([0:xx:xx*imax],[0:yy:yy*jmax],TT)
h=colorbar;
colormap jet;
ylabel(h,'Temperature, C')
set(gca,'YDir','normal')
xlabel('Width, cm')
ylabel('Height, cm')

The following function is created in MATLAB to solve the system of algebraic equations using Gauss elimination method:

function x = GaussElimination(A,b)
% Gauss Elimination without Pivoting
% x = GaussElimination(A,b)
% input:
% A = Coefficient Matrix
% b = Right Hand Side Vector
% output:
% x = Solution Vector

[m,n] = size(A);
if m~=n
error('Matrix A must be square');
end
nb = n+1;
ATM = [A b];

% Forward Elimination
for k = 1:n-1
for i = k+1:n
factor = ATM(i,k)/ATM(k,k);
ATM(i,k:nb) = ATM(i,k:nb)-factor*ATM(k,k:nb);
end
end

% Backward Substitution
x = zeros(n,1);
x(n) = ATM(n,nb)/ATM(n,n);
for i = n-1:-1:1
x(i) = (ATM(i,nb)-ATM(i,i+1:n)*x(i+1:n))/ATM(i,i);
end

disp('Solution Vector')
disp(x)
disp('Coefficient Matrix')
a1 = ATM(:,1:n);
disp(a1)
disp('The b Vector')
b1 = ATM(:,nb);
disp(b1)
disp('Verification of the Solution')
C = a1*x-b1;
disp(C)

Program Output:

Solution Vector
1.3046
3.3734
7.4183
15.8087
33.4596
1.8450
4.7707
10.4910
22.3568
47.3190
1.3046
3.3734
7.4183
15.8087
33.4596

Coefficient Matrix
Columns 1 through 11

-4.0000    1.0000         0         0         0    1.0000         0         0         0         0         0
0   -3.7500    1.0000         0         0    0.2500    1.0000         0         0         0         0
0         0   -3.7333    1.0000         0    0.0667    0.2667    1.0000         0         0         0
0         0         0   -3.7321    1.0000    0.0179    0.0714    0.2679    1.0000         0         0
0         0         0         0   -3.7321    0.0048    0.0191    0.0718    0.2679    1.0000         0
0         0         0         0         0   -3.7321    1.0718    0.0192    0.0051    0.0013    1.0000
0         0         0         0         0         0   -3.4050    1.0824    0.0220    0.0055    0.2872
0         0         0         0         0         0         0   -3.3673    1.0839    0.0210    0.0964
0         0         0         0         0         0         0         0   -3.3638    1.0786    0.0343
0         0         0         0         0         0         0         0         0   -3.3861    0.0124
0         0         0         0         0         0         0         0         0         0   -3.7047
0         0         0         0         0         0         0         0         0         0         0
0         0         0         0         0         0         0         0         0         0         0
0         0         0         0         0         0         0         0         0         0         0
0         0         0         0         0         0         0         0         0         0         0

Columns 12 through 15

0         0         0         0
0         0         0         0
0         0         0         0
0         0         0         0
0         0         0         0
0         0         0         0
1.0000         0         0         0
0.3179    1.0000         0         0
0.1088    0.3219    1.0000         0
0.0385    0.1094    0.3206    1.0000
1.0947    0.0323    0.0114    0.0037
-3.3489    1.1156    0.0393    0.0124
0   -3.2968    1.1193    0.0365
0         0   -3.2919    1.1072
0         0         0   -3.3318

The b Vector
0
0
0
0
-70.7107
-0.0907
-0.3887
-1.4838
-5.5569
-120.7386
-0.5983
-1.9828
-5.5408
-14.9919
-111.4802

Verification of the Solution
1.0e-13 *

-0.0022
0
0
-0.0711
0.1421
0.0053
-0.0078
-0.0644
-0.1421
0
-0.0011
-0.0044
0.0089
0.0355
0

The following plot shows the desired temperature distribution within the plate.