Let us assume that there is a temperature distribution across a rectangular plate which is defined by the following unsteady one-dimensional convection-diffusion equation.
Tt + uTx = αTxx
T(x, 0) = 100 x / L; 0 ≤ x ≤ L
And, the boundary conditions are,
T(0, t) = 0; T(L, t) = 100;
Assume, u = 0.1 cm/s. Develop a MATLAB program that solves this problem with the Backward Time Centered Space (BTCS) Approach. Consider, Δx = 0.1 and ΔT = 0.5. Plot the results for various times steps, such as, 0.5, 2.5, 5, 10, and 50 seconds. The Forward Time Centered Space (FTCS) approach has been shown here.
Backward Time Centered Space (BTCS) Approach:
The following MATLAB program implements the BTCS approach to solve the unsteady one-dimensional convection-diffusion equation using the input parameters and boundary conditions.
%
Backward Time Centered Space (BTCS) Approach
close all;
clc;
clc;
%
Input Properties from Problem
L=1.0; %
(dimension in cm)
alpha=0.01;
u=0.10;
Tt=50.0; % Maximum simulation time
%
Mesh/Grid size calculation
%
Gird size definition along x and y axes
dx=0.1;
dt=0.5;
c=u*dt/dx;
d=alpha*(dt/dx^2);
%
No of grid points
imax=(L/dx)+1;
% Array dimension along the width
tmax=Tt/dt;
% Array dimension along the height
ndim=imax-2;
%
Coefficient Matrix Generation for Gauss-Seidel Method
for
j=1:imax-2
for
i=1:imax-2
if i==j
a(i,j)=1.0+2.0*d;end
if j==i+1
a(i,j)=c/2.0-d;end
if j==i-1
a(i,j)=-(c/2.0+d);end
end
end
for i=1:imax
x(i)=(i-1)*dx;
T(i,1)=100*x(i)/L;
end
for t=1:tmax
time(t+1)=t*dt;
for
i=1:imax-2
b(i)=T(i+1,t);
end
b(1)=b(1)+(c/2.0+d)*T(1,t);
b(imax-2)=b(imax-2)-(c/2.0-d)*T(imax-2,t);
% Calling Gauss-Seidel
function defined another 'm-file'
T(2:ndim+1,t+1)=GaussSeidelMethod(a,b);
T(1,t+1)=T(1,1);
T(imax,t+1)=T(imax,1);
end
hold on
plot(x,T(:,1),'LineWidth',2)
plot(x,T(:,2),'-p','LineWidth',2)
plot(x,T(:,6),'-*','LineWidth',2)
plot(x,T(:,11),'-s','LineWidth',2)
plot(x,T(:,21),'--','LineWidth',2)
plot(x,T(:,100),'^-','LineWidth',2)
xlabel('x, cm')
ylabel('Temperature,
\circ C')
legend('t=0s','t=0.5s','t=2.5s','t=5s', 't=10s', 't=50s')
To solve the algebraic equations in BTCS
method, the following MATLAB function develops Gauss-Seidel algorithm to solve
the equations.
function x =
GaussSeidelMethod(A,b,ET,maxiteration)
% x
= GaussSeidelMethod(A,b): Gauss Seidel without relaxation
%
Input:
% A
= Coefficient Matrix
% b
= Right Side Vector
%
ET = Error Tolerance as Stop Criterion (default in program = 0.00001%)
%
maxiteration = Maximum Iterations (default in program = 100)
%
Output:
% x
= Solution Vector
%
The following logic determines if the given number of inputs are less
%
than required
if nargin
< 2
error('Minimum 2 Input Arguments Required'),
end
%
The following condition checks the maximum number of iterations;
%
If not given, this program considers maxiteration = 100
if nargin
< 4 || isempty(maxiteration)
maxiteration = 100;
end
%
The following loop is for error tolerance for solution;
%
If none mentioned, this program considers it as 0.00001
if
nargin<3||isempty(ET)
ET=0.00001;
end
%
Gauss-Seidel Iteration Approach
[m,n] =
size(A);
if m~=n
error('A Must be a Square
Matrix');
end
C = A;
format short
for i = 1:n
C(i,i) =
0;
x(i) = 0;
end
x = x';
for i = 1:n
C(i,1:n)
= C(i,1:n)/A(i,i);
end
for i = 1:n
d(i) =
b(i)/A(i,i);
end
iteration
= 0;
while (1)
xold =
x; %
Value of x is being updated simultaneously
for i = 1:n
x(i) =
d(i)-C(i,:)*x;
if x(i) ~=
0
r =
abs(b-A*xold); %
Equation for Residual
end
end
iteration
= iteration+1; %
Number of iterations updated
%End
Loop Condition
if
max(abs(b-A*xold)) <= ET | iteration >= maxiteration
break,
end
end
Program Output:
No comments:
Post a Comment