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).
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];
TT =
padarray(TT,[1 0],'pre');
TT =
padarray(TT, [0 1]);
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
No comments:
Post a Comment